“import numpy as np
import pickle
import time
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from pymoo.core.problem import ElementwiseProblem
from pymoo.algorithms.moo.nsga2 import NSGA2
from pymoo.core.callback import Callback
from pymoo.optimize import minimize
from pymoo.indicators.hv import Hypervolume
from pymoo.core.population import Population
from pymoo.core.variable import Real, Integer
from pymoo.core.mixed import MixedVariableMating, MixedVariableSampling, MixedVariableDuplicateElimination
from pymoo.operators.mutation.pm import PolynomialMutation
from pymoo.operators.crossover.sbx import SBX
class OptimProblem(ElementwiseProblem):
def __init__(self):
# 定义变量边界（48个变量，移除了x44/x47/x49）
self.vars = {
"x01": Real(bounds=(1538, 1870)), # S101-M
"x02": Real(bounds=(20148, 24544)), # S104-M
"x03": Real(bounds=(1827.5, 2179.5)), # S108-M
"x04": Real(bounds=(7863.51, 9538.54)), # S112-M
"x05": Real(bounds=(299.604, 364.467)), # S201-M
"x06": Real(bounds=(84837.5, 103471)), # S313-M
"x07": Real(bounds=(13.0821, 15.9591)), # S313-T
"x08": Real(bounds=(426.6, 521.4)), # S313-P
"x09": Real(bounds=(18857.3, 22984.1)), # H101-M
"x10": Real(bounds=(288.71, 350.558)), # H101-COT
"x11": Real(bounds=(18040.6, 21986.5)), # H102-M
"x12": Real(bounds=(540.196, 659.749)), # H102-COT
"x13": Real(bounds=(10632, 12973)), # H201-M
"x14": Real(bounds=(288.167, 351.655)), # H201-HOT
"x15": Real(bounds=(14517.8, 17590.1)), # H301-M
"x16": Real(bounds=(58.271, 70.9577)), # H301-HOT
"x17": Real(bounds=(72029.7, 87790.8)), # H302-M
"x18": Real(bounds=(57.225, 69.1627)), # H302-COT
"x19": Real(bounds=(6685.61, 8148.55)), # H303-M
"x20": Real(bounds=(-9.29859, -7.62888)), # H303-HOT
"x21": Real(bounds=(900.657, 1098.62)), # RE101-C
"x22": Real(bounds=(9.02238, 10.8986)), # RE101-H
"x23": Real(bounds=(1.35287, 1.64148)), # RE101-D
"x24": Real(bounds=(910.653, 1101.31)), # RE102-C
"x25": Real(bounds=(7.2783, 8.7985)), # RE102-H
"x26": Real(bounds=(0.722521, 0.876798)), # RE102-D
"x27": Real(bounds=(1351.32, 1647.75)), # RE103-C
"x28": Real(bounds=(7.22507, 8.77916)), # RE103-H
"x29": Real(bounds=(1.38167, 1.61998)), # RE103-D
"x30": Real(bounds=(1.80656, 2.19881)), # RE201-P
"x31": Real(bounds=(6.31061, 7.6632)), # RE201-L
"x32": Real(bounds=(0.453803, 0.549121)), # RE201-D
"x33": Real(bounds=(1.80579, 2.19799)), # F201-P
"x34": Real(bounds=(135.366, 164.98)), # F201-T
"x35": Real(bounds=(1.75542, 2.14253)), # AB301-P1
"x36": Integer(bounds=(11, 14)), # AB301-N
"x37": Real(bounds=(0.0930186, 0.113473)), # R301-P1
"x38": Real(bounds=(1.81906, 2.21977)), # R301-R
"x39": Real(bounds=(5889250, 6683890)), # R301-QN
"x40": Real(bounds=(0.757167, 0.923944)), # R302-P1
"x41": Real(bounds=(-71041.4, -58147.9)), # R302-Q1
"x42": Integer(bounds=(7, 8)), # R301-N
"x43": Real(bounds=(0.108235, 0.130973)), # V301-P
# "x44": Integer(bounds=(8, 9)), # R301-FEEDS
"x45": Integer(bounds=(7, 8)), # R302-N
"x46": Real(bounds=(0.772167, 0.941046)), # V302-P
# "x47": Integer(bounds=(8, 9)), # R302-FEEDS
"x48": Real(bounds=(1.78042, 2.17364)), # V303-P
# "x49": Integer(bounds=(12, 15)), # AB301-FEEDS
"x50": Real(bounds=(58, 71)), # F301-T
"x51": Real(bounds=(1.98292, 2.37614)), # F301-P
}
# 预加载模型（避免重复磁盘I/O）
self.models = {
'TAC': pickle.load(open("TAC.dat", "rb")),
'CO2R': pickle.load(open("CO2R.dat", "rb")),
'HD': pickle.load(open("HD.dat", "rb")),
'C1C2R': pickle.load(open("C1C2R.dat", "rb")),
'RE101Y': pickle.load(open("RE101Y.dat", "rb")),
'RE102Y': pickle.load(open("RE102Y.dat", "rb")),
'RE103Y': pickle.load(open("RE103Y.dat", "rb")),
'RE201Y': pickle.load(open("RE201Y.dat", "rb"))
}
super().__init__(vars=self.vars, n_obj=3, n_constr=5) # 修正约束数量
def _evaluate(self, x, out, *args, **kwargs):
# 1. 构建完整51维输入向量
input_vector = self._construct_input_vector(x)
# 2. 预测目标值
TAC, HD, C1C2R = self._predict_objectives(input_vector)
# 3. 预测约束值
CO2R, RE101Y, RE102Y, RE103Y, RE201Y = self._predict_constraints(input_vector)
# 4. 设置目标函数
out["F"] = np.array([
TAC, # 最小化TAC（原值）
# -CO2R, # 最小化负CO2纯度（相当于最大化CO2R）
HD, # 最小化HD（原值）
-C1C2R # 最小化负C1-C2回收（相当于最大化C1C2R）
])
# 5. 设置约束条件（g <= 0）
tolerance = 1e-6 # 约束容差
out["G"] = np.array([
0.90 - CO2R + tolerance, # 允许轻微违反约束
0.90 - RE101Y + tolerance,
0.90 - RE102Y + tolerance,
0.90 - RE103Y + tolerance,
0.90 - RE201Y + tolerance
])
def _construct_input_vector(self, x):
"""构建51维输入向量，动态计算x44/x47/x49"""
# 基础向量（按顺序提取48个变量）
base = np.array([x[f"x{k:02}"] for k in range(1, 52) if f"x{k:02}" in self.vars])
# 插入计算得到的变量（基于优化规则）
x44_val = x["x42"] + 1 # x44 = x42 + 1 R301-FEEDS
x47_val = x["x45"] + 1 # x47 = x45 + 1 R302-FEEDS
x49_val = x["x36"] + 1 # x49 = x36 + 1 AB301-FEEDS
# 在指定位置插入计算值
return np.insert(base, [43, 46, 48], [x44_val, x47_val, x49_val])
def _predict_objectives(self, X):
"""预测三个目标值"""
X = X.reshape(1, -1)
return (
self.models['TAC'].predict(X)[0],
# self.models['CO2R'].predict(X)[0],
self.models['HD'].predict(X)[0],
self.models['C1C2R'].predict(X)[0]
)
def _predict_constraints(self, X):
"""预测五个约束值"""
X = X.reshape(1, -1)
return (
self.models['CO2R'].predict(X)[0],
self.models['RE101Y'].predict(X)[0],
self.models['RE102Y'].predict(X)[0],
self.models['RE103Y'].predict(X)[0],
self.models['RE201Y'].predict(X)[0]
)
# 优化配置
def setup_algorithm(pop_size=200):
return NSGA2(
pop_size=pop_size,
sampling=MixedVariableSampling(),
crossover=SBX(prob=0.85, eta=20),# 降低交叉概率，增加分布指数
mutation=PolynomialMutation(prob=0.7, eta=30),# 调整变异概率
mating=MixedVariableMating(eliminate_duplicates=MixedVariableDuplicateElimination()),
eliminate_duplicates=MixedVariableDuplicateElimination()
)
# 运行优化
if __name__ == "__main__":
start_time = time.time()
problem = OptimProblem()
algorithm = setup_algorithm()
res = minimize(problem, algorithm, ('n_gen', 100), seed=42, save_history=True, verbose=True)
# 提取最优结果
opt_X = np.array([[sol.X[k] for k in problem.vars] for sol in res.opt])
opt_F = res.F
# 保存结果
np.savetxt("optimal_variables.csv", opt_X, delimiter=",", header=",".join(problem.vars.keys()), comments='')
np.savetxt("optimal_objectives.csv", opt_F, delimiter=",", header="-TAC,-HD,C1C2R", comments='')
# 可视化三个目标对
fig, ax = plt.subplots(1, 1, figsize=(8, 6))
scatter = ax.scatter(opt_F[:, 0], opt_F[:, 2], c=-opt_F[:, 1], cmap="viridis")
ax.set_title("TAC vs HD (color=C1C2R)")
ax.set_xlabel("TAC")
ax.set_ylabel("HD")
plt.colorbar(scatter, ax=ax)
plt.tight_layout()
plt.savefig("pareto_front_3metrics.png", dpi=900)
plt.show()
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111, projection='3d')
ax.scatter(opt_F[:, 0], -opt_F[:, 1], opt_F[:, 2], c='teal', s=40)
ax.set_xlabel('TAC')
ax.set_ylabel('C1C2R')
ax.set_zlabel('HD')
ax.set_title('3D Pareto Front: TAC vs CO2R vs HD')
plt.tight_layout()
plt.savefig("pareto_3d_TAC_C1C2R_HD.png", dpi=900)
plt.show()
# 生成所有二维目标对图像
plt.figure(figsize=(6,5))
plt.scatter(opt_F[:, 0], -opt_F[:, 1], c='darkblue')
plt.xlabel("TAC")
plt.ylabel("C1C2R")
plt.title("TAC vs C1C2R")
plt.tight_layout()
plt.savefig("pareto_TAC_vs_C1C2R.png", dpi=900)
plt.show()
plt.figure(figsize=(6,5))
plt.scatter(opt_F[:, 0], opt_F[:, 2], c='darkred')
plt.xlabel("TAC")
plt.ylabel("HD")
plt.title("TAC vs HD")
plt.tight_layout()
plt.savefig("pareto_TAC_vs_HD.png", dpi=900)
plt.show()
plt.figure(figsize=(6,5))
plt.scatter(-opt_F[:, 1], opt_F[:, 2], c='darkgreen')
plt.xlabel("C1C2R")
plt.ylabel("HD")
plt.title("C1C2R vs HD")
plt.tight_layout()
plt.savefig("pareto_C1C2R_vs_HD.png", dpi=900)
plt.show()
# 自动选最优解：距离理想点最近
# 仅从 C1C2R ≥ 0.80 的解中推荐最优解
valid_idx = np.where(-opt_F[:, 1] >= 0.80)[0]
if len(valid_idx) > 0:
subset = opt_F[valid_idx]
ideal_point = np.array([min(subset[:, 0]), max(subset[:, 1]), min(subset[:, 2])])
normalized = np.column_stack([subset[:, 0], -subset[:, 1], subset[:, 2]])
distances = np.linalg.norm(normalized - ideal_point, axis=1)
best_idx_in_subset = valid_idx[np.argmin(distances)]
print("\n==== 推荐解====")
print("索引:", best_idx_in_subset)
print("变量配置:")
for name, val in zip(problem.vars.keys(), opt_X[best_idx_in_subset]):
print(f" {name}: {val:.4f}")
print("对应目标函数值:")
print(f" TAC: {opt_F[best_idx_in_subset, 0]:.4f}")
print(f" C1C2R: {-opt_F[best_idx_in_subset, 1]:.4f}")
print(f" HD: {opt_F[best_idx_in_subset, 2]:.4f}")”这是我的代码，但是我不想规定C1C2R一定要大于多少了，请为我修改嗲吗
扫一扫
1466
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
