大家好,这里是七七,本次的Python学习的例子是NSGAII优化算法的实现代码
目录
import numpy as np
import pandas as pd
from tool import *
import deap as ea
address=r'data/q3.csv'
data=pd.read_csv(address)
data_40_rows=data.head(40)
cost_pre=data_40_rows['批发价格(元/千克)'].tolist() #预测成本
sale_pre=data_40_rows['平均销量'].tolist() #预测销量
for i in range(40):
if sale_pre[i]<=2.5:
sale_pre[i]=2.5
mean_sale=data_40_rows['平均单价']
alph=[0.1]*40 #损耗率
discount=[0.8]*40 #打折
ini_data=np.array([i]*30+[0]*10+[0.7]*40+sale_pre)
need_list=data_40_rows['Q'].tolist()
weight=[0.5]*40
def Modify(sale_price):
sale=[]
for i in range(40):
ss=(1-(sale_price[i]-mean_sale[i])/(mean_sale[i])*weight[i])+sale_pre[i]
sale.append(ss)
return sale
def Obj_function(X):
Objv=[]
Cv=[]
for i in range (len(X)):
x=X[i]
#目标一:最大化利润
decide=x[:40]
profit=x[40:80]
nums=x[80:120]
#计算成本
cost=[]
for i in range(len(cost_pre)):
cost.append(decide[i]*cost_pre[i])
#计算售价
sale_price=[]
for i in range(len(cost_pre)):
s=cost[i]*(1+profit[i])
sale_price.append(s)
#计算销量
sale_modify=Modify(sale_price)
#计算总销售额
w1=[]
for i in range(len(sale_modify)):
bad=(sale_modify[i]*alph[i])*(discount[i]*sale_price[i])
good=(sale_modify[i]*(1-alph[i]))*(sale_price[i])
w1.append(bad+good)
sum_sale=sum(w1)
#计算总成本
w2=[]
for i in range(len(cost)):
w2.append(nums[i]*cost[i])
sum_cost=sum(w2)
penalty=0
for i in range(40):
if nums[i]<2.5:
penalty-=200
f1=sum_sale-sum_cost+penalty
#目标二:最大化需求
satif_need=[]
for i in range(len(decide)):
satif_need.append(decide[i]*need_list[i])
f2=sum(satif_need)/sum(need_list)
#确定评价函数值
x_Objv=np.hstack([f1,f2])
#定义约束
list_cv=[sum(decide)-33,27-sum(decide)]
x_Cv=np.array(list_cv)
Objv.append(x_Objv)
Cv.append(x_Cv)
Objv=np.array(Objv)
Cv=np.array(Cv)
return Objv,Cv
problem=ea.Problem(
name='NSGAII_for_q3',
M=2,
maxormins=[-1,-1],
Dim=120,
varTypes=[1]*40+[0]*80,
lb=[0]*40+[0]*40+ data_40_rows['最小销量'].tolist(),
ub=[1]*40+[1.2]*40+ data_40_rows['最大销量'].tolist(),
evalVars=Obj_function
)
#混合编码
Encodings=['BG','RI']
Field1=ea.crtfld(Encodings[0],problem.varTypes[:40],ranges=np.array([problem.lb[:40],problem.ub[:40]]))
Field2=ea.crtfld(Encodings[1],problem.varTypes[:40],ranges=np.array([problem.lb[:40],problem.ub[:40]]))
Fields=[Field1,Field2]
population=ea.PsyPopulation(Encodings=Encodings,Fields=Fields,NIND=100)
algorithm=ea.moea_psy_NSGA2_archive_templet(
problem,
population,
MAXGEN=300,
logTras=20,
prophetPop=ini_data,
maxTrappedCount=20,
)
res=ea.optimize(algorithm,seed=1,verbose=True,drawing=1,
outputMsg=True,drawLog=True)
print(f"最优解是:{res['Vars'][0]}")
print(f"最优解值是:{res['Vars'][0]}")
print(f"hv结果是:{res['Vars'][0]}")
ini_data=np.array([1]*30+[0]*10+[0.7]*40+sale_pre)
ini_data
?是一个变量,它是一个由 Numpy 库创建的一维数组。这个数组包含了四个部分的数据,分别是值为 1 的 30 个元素、值为 0 的 10 个元素、值为 0.7 的 40 个元素和?sale_pre
?中存储的元素。
x_Objv=np.hstack([f1,f2])
x_Objv
?是一个 numpy 数组,它通过?np.hstack()
?函数将两个数组?f1
?和?f2
?水平堆叠在一起。
这意味着?f1
?和?f2
?必须具有相同的维度,并且它们将在水平方向上进行堆叠,形成一个新的?x_Objv
?数组。换句话说,x_Objv
?的每一列将从?f1
?和?f2
?中对应的列中获取元素。
problem=ea.Problem(
name='NSGAII_for_q3',
M=2,
maxormins=[-1,-1],
Dim=120,
varTypes=[1]*40+[0]*80,
lb=[0]*40+[0]*40+ data_40_rows['最小销量'].tolist(),
ub=[1]*40+[1.2]*40+ data_40_rows['最大销量'].tolist(),
evalVars=Obj_function
)
Problem
?是在?ea
?模块中创建一个优化问题的类的对象。通过使用该对象,你可以定义和处理多目标优化问题。
在?Problem
?函数调用中,传递了一些参数来设置问题的属性和限制条件:
name='NSGAII_for_q3'
:设置问题的名称为?NSGAII_for_q3
。M=2
:表示问题有两个目标函数。maxormins=[-1,-1]
:表示两个目标函数都是最小化目标。Dim=120
:设置问题的维度为 120。varTypes=[1]*40+[0]*80
:指定了每个变量的类型,其中前 40 个变量是连续变量(对应?1
),后 80 个变量是离散变量(对应?0
)。lb=[0]*40+[0]*40+ data_40_rows['最小销量'].tolist()
:设置了问题的变量下界(lower bound)。前 40 个变量的下界都为 0,接下来的 40 个变量仍为 0,最后 40 个变量是从?data_40_rows
?中取出的一列名为?最小销量
?的数据转换为列表。ub=[1]*40+[1.2]*40+ data_40_rows['最大销量'].tolist()
:设置了问题的变量上界(upper bound)。前 40 个变量的上界都为 1,接下来的 40 个变量都为 1.2,最后 40 个变量是从?data_40_rows
?中取出的一列名为?最大销量
?的数据转换为列表。evalVars=Obj_function
:指定了用于评估变量的目标函数。
通过调用?Problem
?函数,你创建了一个名为?NSGAII_for_q3
?的优化问题对象,并设置了相关的属性和限制条件。
Encodings=['BG','RI']
Field1=ea.crtfld(Encodings[0],problem.varTypes[:40],ranges=np.array([problem.lb[:40],problem.ub[:40]]))
Field2=ea.crtfld(Encodings[1],problem.varTypes[:40],ranges=np.array([problem.lb[:40],problem.ub[:40]]))
Fields=[Field1,Field2]
population=ea.PsyPopulation(Encodings=Encodings,Fields=Fields,NIND=100)
首先定义了一个名为?Encodings
?的列表,其中包含了两个编码类型:“BG” 和 “RI”。
接下来,使用?ea.crtfld()
?函数创建了两个字段对象?Field1
?和?Field2
。这两个字段对象分别使用了?Encodings
?列表中的第一个和第二个编码类型。
ea.crtfld()
?函数的参数如下:
- 第一个参数是编码类型,可以是 “BG”(二元编码)或 “RI”(实数编码)。
- 第二个参数是一个列表,用于指定每个变量的类型。在这里,使用?
problem.varTypes[:40]
?取出了问题对象中前 40 个变量的类型。- 第三个参数是一个范围 (
ranges
) 数组,用于指定每个变量的取值范围。在这里,使用了?np.array([problem.lb[:40], problem.ub[:40]])
,取出了问题对象中前 40 个变量的下界(problem.lb[:40]
)和上界(problem.ub[:40]
)。
然后,将?Field1
?和?Field2
?添加到一个列表?Fields
?中。
最后,通过调用?ea.PsyPopulation()
?函数创建了一个名为?population
?的粒子种群对象。ea.PsyPopulation()
?函数的参数如下:
Encodings
:表示个体的编码类型列表。Fields
:表示个体所使用的字段列表。NIND
:表示种群中的个体数。
在这里,将之前定义的?Encodings
?列表和?Fields
?列表作为参数传递给?ea.PsyPopulation()
?函数,并设置了种群中的个体数为 100。
通过以上代码,你创建了一个带有两种编码的粒子种群,每个编码类型对应一个字段对象。?
algorithm=ea.moea_psy_NSGA2_archive_templet(
problem,
population,
MAXGEN=300,
logTras=20,
prophetPop=ini_data,
maxTrappedCount=20,
)
res=ea.optimize(algorithm,seed=1,verbose=True,drawing=1,
outputMsg=True,drawLog=True)
首先使用?ea.moea_psy_NSGA2_archive_templet()
?函数创建了一个名为?algorithm
?的多目标优化算法对象。这个算法是基于 NSGA-II 算法的带存档的粒子群优化算法模板。
ea.moea_psy_NSGA2_archive_templet()
?函数的参数如下:
- 第一个参数?
problem
:表示要解决的优化问题对象。- 第二个参数?
population
:表示种群对象。MAXGEN
:表示算法的最大进化代数。logTras
:表示在每个进化代数之后,将群体全部划分成多少个轨迹段。prophetPop
:预知种群,表示具有先验知识或已知最优解的种群,用于指导群体的进化。maxTrappedCount
:表示算法最大停滞计数值,用于提前终止算法。
然后,将?algorithm
?对象传递给?ea.optimize()
?函数,以进行优化操作。ea.optimize()
?函数的参数如下:
- 第一个参数?
algorithm
:表示要使用的优化算法对象。seed
:表示随机数种子,用于控制算法的随机性。verbose
:表示是否显示详细信息。drawing
:表示是否绘制图形。outputMsg
:表示是否输出算法的运行信息。drawLog
:表示是否绘制优化过程的日志。
通过以上代码,你使用了 NSGA-II 算法的带存档的粒子群优化算法模板来进行多目标优化。