2025年美赛B题加急版论文首发+B题补充数据集

import numpy as np
import matplotlib.pyplot as plt

from pymoo.core.problem import ElementwiseProblem
from pymoo.algorithms.moo.nsga2 import NSGA2
from pymoo.termination import get_termination
from pymoo.optimize import minimize

def multiobj_sustainable_tourism():
    """
    多目标优化示例：可持续旅游 (Python + pymoo)
    对标 MATLAB gamultiobj 示例。
    """

    # 1. 定义基准值 (以 2021 年为例)
    RP_baseline = 31973.0       # 居住人口
    MI_baseline = 45148.0       # 个人收入中位数
    TR_baseline = 1271412500.0  # 收入总值

    lb = np.array([0.9*RP_baseline, 0.9*MI_baseline, 0.9*TR_baseline])
    ub = np.array([1.1*RP_baseline, 1.1*MI_baseline, 1.1*TR_baseline])

    # 2. 定义多目标问题 (继承自 ElementwiseProblem)
    class TourismProblem(ElementwiseProblem):
        def __init__(self):
            super().__init__(
                n_var=3,           # 决策变量维度
                n_obj=2,           # 目标函数个数
                n_ieq_constr=0,    # 不等式约束个数
                n_eq_constr=1,     # 等式约束个数
                xl=lb,             # 决策变量下界
                xu=ub,             # 决策变量上界
                type_var=np.float64
            )
            # 一些基准值放这里，便于在 _evaluate 中调用
            self.TR_baseline = TR_baseline
            self.baselineEnv = 100.0
            self.alpha = 1e-9

        def _evaluate(self, x, out, *args, **kwargs):
            """
            x: [RP, MI, TR]
            目标1: f1 = 0.006816 * TR - 9.3361e+06 (相当于 -TCT, 需要最小化)
            目标2: f2 = 100 + 1e-9*(TR - TR_baseline) (环境影响)
            约束:  TR - [5.8475e+08 - 6131.7*RP + 19581*MI] = 0
            """
            RP, MI, TR = x[0], x[1], x[2]

            # ---- 目标函数 ----
            f1 = 0.006816 * TR - 9.3361e+06
            f2 = self.baselineEnv + self.alpha*(TR - self.TR_baseline)

            # ---- 等式约束 (h = 0) ----
            #   TR - (5.8475e+08 - 6131.7*RP + 19581*MI) = 0
            model2_val = 5.8475e+08 - 6131.7*RP + 19581.0*MI
            h = TR - model2_val

            # 给 pymoo 输出
            out["F"] = [f1, f2]    # 两个目标
            out["H"] = [h]         # 1 个等式约束

    # 初始化问题
    problem = TourismProblem()

    # 3. 设置算法参数，相当于 gamultiobj 的选项
    algorithm = NSGA2(
        pop_size=100,          # 种群规模
        # 可根据需要增加其他参数设置
    )

    # 4. 设置终止准则 (最大迭代代数)
    termination = get_termination("n_gen", 100)

    # 5. 开始优化
    res = minimize(
        problem,
        algorithm,
        termination,
        seed=1,
        save_history=True,
        verbose=True     # 在控制台输出迭代信息
    )

    # 6. 结果分析
    X = res.X   # 帕累托前沿上的决策变量集合 (可能有多行)
    F = res.F   # 对应的目标函数值 [f1, f2]

    print("=== 帕累托解 (局部查看) ===")
    for i in range(len(X)):
        RP_sol, MI_sol, TR_sol = X[i]
        f1_val, f2_val = F[i]
        print(f"Solution {i+1}:  RP={RP_sol:.2f}, MI={MI_sol:.2f}, TR={TR_sol:.2f}, "
              f"Obj1={f1_val:.2f}, Obj2={f2_val:.6f}")

    print(f"\nExit message: {res.algorithm.verbose}")

    # 7. 绘制帕累托前沿
    # f1 在横轴，f2 在纵轴
    plt.figure(figsize=(6,5))
    plt.scatter(F[:, 0], F[:, 1], c='red', marker='o')
    plt.xlabel('f1 = 0.006816*TR - 9.3361e+06  (minimize)')
    plt.ylabel('f2 = EnvImpact             (minimize)')
    plt.title('Pareto Front (NSGA2 Result)')
    plt.grid(True)
    plt.show()


# 入口
if __name__ == "__main__":
    multiobj_sustainable_tourism()