from pymoo.core.problem import ElementwiseProblem
from pymoo.algorithms.moo.nsga2 import NSGA2
from pymoo.optimize import minimize
import numpy as np
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
import time
import pickle
start = time.time()
class MyProblem(ElementwiseProblem):
def __init__(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
}
super().__init__(vars=vars, n_obj=4, n_constr=7) # 目标数和约束数
def _evaluate(self, x, out, *args, **kwargs):
# 提取 51 个特征
x = np.array([x[f"x{k:02}"] for k in range(1, 52)])
# 获取预测输出
output = self.XGB_model(x)
f1 = -output[0] # 最小化TAC
f2 = output[1] # 最大化塔顶CO2纯度
f3 = -output[2] # 最小化能耗
f4 = output[3] # 最大化C1-C2回收
# 获取约束条件
g_values = self.Constraint_model(x)
# 提取整数类型变量
x36 = int(x[35]) # 由于索引从 0 开始,x36 在第36个位置是 x[35]
x49 = int(x[48]) # x49 对应位置为 x[48]
x42 = int(x[41]) # x42 对应位置为 x[41]
x44 = int(x[43]) # x44 对应位置为 x[43]
x45 = int(x[44]) # x45 对应位置为 x[44]
x47 = int(x[46]) # x47 对应位置为 x[46]
# 约束计算
g1 = 0.90 - g_values[0] # 当预测值<0.9时g1>0(违反约束)
g2 = 0.90 - g_values[1]
g3 = 0.90 - g_values[2]
g4 = 0.90 - g_values[3]
# 根据规则修改 x49, x44 和 x47
x49 = x36 + 1 # x49 = x36 + 1
x44 = x42 + 1 # x44 = x42 + 1
x47 = x45 + 1 # x47 = x45 + 1
# 目标和约束输出
out["F"] = np.column_stack([f1, f2, f3, f4]) # 目标
out["G"] = np.column_stack([g1, g2, g3, g4, x49, x44, x47]) # 约束
def XGB_model(self, Xtest):
Xtest = Xtest.reshape(1, -1)
# 加载XGBoost模型
etr1 = pickle.load(open("TAC.dat", "rb"))
etr2 = pickle.load(open("CO2R.dat", "rb"))
etr3 = pickle.load(open("HD.dat", "rb"))
etr4 = pickle.load(open("C1C2R.dat", "rb"))
# 获取预测值
TAC_test_predict = etr1.predict(Xtest)
CO2R_test_predict = etr2.predict(Xtest)
HD_test_predict = etr3.predict(Xtest)
C1C2R_test_predict = etr4.predict(Xtest)
return TAC_test_predict, CO2R_test_predict, HD_test_predict, C1C2R_test_predict
def Constraint_model(self, Xtest):
Xtest = Xtest.reshape(1, -1)
# 加载约束模型
etr5 = pickle.load(open("RE101Y.dat", "rb"))
etr6 = pickle.load(open("RE102Y.dat", "rb"))
etr7 = pickle.load(open("RE103Y.dat", "rb"))
etr8 = pickle.load(open("RE201Y.dat", "rb"))
# 获取预测值
RE101Y_test_predict = etr5.predict(Xtest)
RE102Y_test_predict = etr6.predict(Xtest)
RE103Y_test_predict = etr7.predict(Xtest)
RE201Y_test_predict = etr8.predict(Xtest)
return RE101Y_test_predict, RE102Y_test_predict, RE103Y_test_predict, RE201Y_test_predict
# 优化算法设置
algorithm = NSGA2(
pop_size=100,
n_offsprings=10,
sampling=MixedVariableSampling(),
crossover=SBX(prob=0.9, eta=15),
mutation=PolynomialMutation(prob=1.0, eta=20),
mating=MixedVariableMating(eliminate_duplicates=MixedVariableDuplicateElimination()),
eliminate_duplicates=MixedVariableDuplicateElimination(),
)
# 最优化问题求解
problem = MyProblem()
res = minimize(problem,
algorithm,
('n_gen', 100),
seed=10,
save_history=True,
verbose=True)
# 获取结果
labels = list(res.opt[0].X.keys())
bounds = np.array([problem.vars[name].bounds for name in labels]).T
X = np.array([[sol.X[name] for name in labels] for sol in res.opt])
F = res.F # 目标
# 合并所有可行解
all_pop = Population()
for algorithm in res.history:
all_pop = Population.merge(all_pop, algorithm.off)
B = all_pop.get("F") # 可行解
# 输出运行时间
end = time.time()
print(end - start)
最新发布