如何在Pymoo中构建优先X0-X3的初始种群及代码修改

我的代码用于实现利润最大化和总面积最小化的多目标优化，其中X0-X3为免费使用且生产率更高的种植区域，X4-X7为有使用成本且生产率较低的区域。需要修改Pymoo代码，让初始种群优先使用X0-X3，仅在必要时使用X4-X7；同时需要了解Pymoo中种群个体的选择机制，以及如何编写代码实现初始种群的优先级设置。

NSGA-II的种群选择核心流程：

1. 初始种群生成：通过指定的sampling策略在变量上下限范围内生成初始个体，默认是均匀随机采样。
2. 锦标赛选择：从种群中随机挑选k个个体，通过支配关系和拥挤度比较选出最优个体作为父代，用于交叉变异。
3. 子代生成：父代通过交叉、变异操作产生新个体。
4. 合并与筛选：父代与子代合并后，通过非支配排序和拥挤度计算筛选出下一代种群，保留pop_size个最优个体。

初始种群的采样是优化起点，要实现优先级设置，核心是自定义采样策略，让初始个体倾向于优先填充X0-X3区域。

通过自定义Sampling类实现初始种群对X0-X3的优先使用：

• 大部分初始个体中，X4-X7初始值设为0，仅在X0-X3的最大产能无法满足需求约束时，少量分配X4-X7的面积。
• 保留小比例随机个体，保证种群多样性，避免陷入局部最优。

pip install -U pymoo==0.5.0
# 导入Pymoo相关库
import numpy as np
from pymoo.core.problem import ElementwiseProblem
from pymoo.algorithms.moo.nsga2 import NSGA2
from pymoo.core.sampling import Sampling
from pymoo.factory import get_crossover, get_mutation, get_termination
from pymoo.optimize import minimize
from pymoo.visualization.scatter import Scatter

# 初始参数（单位：R$/t）
S_preco = 2300

# 区域最大面积（单位：ha）
# MA - 1, TO - 2, PI - 3, BA - 4
# 新区域MA -5, TO -6, PI -7, BA -8
AreaMAX = [600000, 650000, 500000, 1500000, 309853, 930190, 14197, 55525]

# 生产率（单位：t/ha）
Produtividade = [3.206, 3.322, 3.377, 3.779, 3, 3.12, 3, 3.3]

# 成本参数
# 大豆种植成本（R$/t）
S_custo = [983, 966.5, 1051.5, 1040.167, 983, 966.5, 1051.5, 1040.167]
# 土地使用成本（R$/t）
Custo_terra = [0, 0, 0, 0, 7500, 9000, 10000, 20500]

# 大豆需求（单位：t）
S_demanda = 5056921 + 2332769 - (5056921/117887685 + 2332769/117887685)*18344777

# 其他产品价格与成本（固定值）
S_preco = 2500
F_preco = 2500
O_preco = 3500
B_preco = 3500

F_custo = 60
O_custo = 100
B_custo = 500

# 副产品需求（单位：t）
F_demanda = 1510028
O_demanda = 378913
B_demanda = 244031.02

# 自定义问题类
class MyProblem(ElementwiseProblem):
    def __init__(self):
        super().__init__(n_var=10,
                         n_obj=2,
                         n_constr=12,
                         xl=np.array([0,0,0,0,0,0,0,0,0,0]),
                         xu=np.array([AreaMAX[0],AreaMAX[1],AreaMAX[2],AreaMAX[3],AreaMAX[4],AreaMAX[5],AreaMAX[6],AreaMAX[7],1,1]))

    def _evaluate(self, X, out, *args, **kwargs):
        # 计算总产量
        total_prod = np.sum(X[:8] * Produtividade)
        # 利润最大化（取负转为最小化问题）
        f1 = -1*((total_prod * X[8] * S_preco) + 
                 (total_prod * (1 - X[8]) * (F_preco*0.4 + 0.2*(X[9]*O_preco + (1 - X[9])*B_preco))) - 
                 np.sum(X[:8] * (np.array(Produtividade)*np.array(S_custo) + np.array(Custo_terra))) - 
                 (total_prod * (1 - X[8]) * (F_custo + O_custo + 0.2*(1 - X[9])*B_custo)))/1e10
        # 总面积最小化
        f2 = np.sum(X[:8])/1e7

        # 约束条件
        g1 = X[0] - AreaMAX[0]
        g2 = X[1] - AreaMAX[1]
        g3 = X[2] - AreaMAX[2]
        g4 = X[3] - AreaMAX[3]
        g5 = X[4] - AreaMAX[4]
        g6 = X[5] - AreaMAX[5]
        g7 = X[6] - AreaMAX[6]
        g8 = X[7] - AreaMAX[7]
        g9 = -1*(total_prod * X[8] - S_demanda)
        g10 = -1*(total_prod * (1 - X[8])*0.4 - F_demanda)
        g11 = -1*(total_prod * (1 - X[8])*0.2*X[9] - O_demanda)
        g12 = -1*(total_prod * (1 - X[8])*0.2*(1 - X[9]) - B_demanda)

        out["F"] = [f1, f2]
        out["G"] = [g1, g2, g3, g4, g5 ,g6, g7, g8, g9, g10, g11, g12]

# 自定义采样器：优先使用X0-X3区域
class PrioritySampling(Sampling):
    def _do(self, problem, n_samples, **kwargs):
        # 初始化种群矩阵
        X = np.random.rand(n_samples, problem.n_var)
        # 变量上下限
        xl, xu = problem.xl, problem.xu

        # 80%的个体优先使用X0-X3，X4-X7初始为0
        priority_ratio = 0.8
        n_priority = int(n_samples * priority_ratio)
        
        # 处理优先种群
        for i in range(n_priority):
            # X0-X3随机分配（不超过最大值）
            X[i, 0] = np.random.uniform(xl[0], xu[0])
            X[i, 1] = np.random.uniform(xl[1], xu[1])
            X[i, 2] = np.random.uniform(xl[2], xu[2])
            X[i, 3] = np.random.uniform(xl[3], xu[3])
            # X4-X7初始为0
            X[i, 4:8] = 0
            # X8和X9随机分配
            X[i, 8] = np.random.uniform(xl[8], xu[8])
            X[i, 9] = np.random.uniform(xl[9], xu[9])

            # 检查总产量是否满足需求约束，若不满足则少量分配X4-X7
            total_prod = np.sum(X[i,:8] * Produtividade)
            # 计算所需最小产量
            if X[i,8] > 0:
                required_prod = max(S_demanda / X[i,8], 
                                   (F_demanda/0.4 + O_demanda/(0.2*X[i,9]) + B_demanda/(0.2*(1-X[i,9])))/3)
            else:
                required_prod = (F_demanda/0.4 + O_demanda/(0.2*X[i,9]) + B_demanda/(0.2*(1-X[i,9])))/3
            
            if total_prod < required_prod:
                deficit = required_prod - total_prod
                # 按X4-X7的最大产能比例分配补充面积
                max_prod_new = np.array(Produtividade[4:8]) * np.array(AreaMAX[4:8])
                total_max_new = np.sum(max_prod_new)
                if total_max_new > 0:
                    ratio = deficit / total_max_new
                    X[i,4:8] = np.array(AreaMAX[4:8]) * ratio
                    # 确保不超过区域最大面积
                    X[i,4:8] = np.minimum(X[i,4:8], np.array(AreaMAX[4:8]))

        # 剩余20%个体使用随机采样，保证多样性
        for i in range(n_priority, n_samples):
            X[i] = np.random.uniform(xl, xu)

        return X

# 初始化问题
problem = MyProblem()

# 配置NSGA-II算法，使用自定义采样器
algorithm = NSGA2(
    pop_size=1100,
    n_offsprings=10,
    sampling=PrioritySampling(),
    crossover=get_crossover("real_sbx", prob=0.9, eta=15),
    mutation=get_mutation("real_pm", eta=20),
    eliminate_duplicates=True
)

# 设置终止条件：1000代
termination = get_termination("n_gen", 1000)

# 运行优化
res = minimize(problem,
               algorithm,
               termination,
               seed=1,
               save_history=True,
               verbose=True)

# 处理结果：将f1转回利润最大化的正值
res.F[:,0] *= -1

# 绘制Pareto前沿
plot = Scatter(title = "目标空间")
plot.add(res.F, color="red")
plot.show()

# 保存种群变量到CSV
X = res.pop.get("X")
np.savetxt("NSGApop.csv", X, delimiter=",")

# 保存目标函数结果到CSV
F = res.pop.get("F")
np.savetxt("NSGAResult.csv", F, delimiter=",")

1. 自定义采样器PrioritySampling：
   • 80%的初始个体优先填充X0-X3区域，X4-X7初始为0；若X0-X3产能无法满足需求，则按X4-X7的最大产能比例分配补充面积。
   • 剩余20%个体使用随机采样，避免种群多样性不足。
2. 算法配置：将原随机采样替换为自定义的PrioritySampling()，确保初始种群符合优先级要求。
3. 目标函数简化：对总产量计算进行了整合，提升代码可读性。

内容的提问来源于stack exchange，提问作者Rafael Henrique