大家好,这里是七七,本次Python学习专题的例子剖析已经接近尾声,再更新几期本专题就结束了,对于展示的例子,可以介绍的细节部分越来越少了。
今天来给大家介绍的是粒子群优化算法的实现代码。
import matplotlib.pyplot as plt
import numpy as np
import pyswarms as ps
plt.rcParams['font.sans-serif'] = [u'simHei']
plt.rcParams['axes.unicode_minus'] = False
list_1=['花叶类','花菜类','水生根茎类','茄类','辣椒类','食用菌']
predit_buy=[
[3.285864,3.2921748,3.2889733,3.285188,3.2851105,3.2964268,3.2876368] ,
[7.7414317,7.763459,7.814592,7.794937,7.747068,7.810813, 7.7633805] ,
[12.018661,11.912668,12.027704,11.941088,11.92054,12.118359,11.972251] ,
[4.5562034,4.601929,4.5483465,4.549116,4.532483,4.539543,4.601603],
[3.7067149,3.65774,3.6644902,3.6755412,3.6658049,3.6834998,3.6471841] ,
[4.015016,4.075036,4.025253,4.0783653,4.0397897,4.0211616,4.0569587] ]#预测进价
predit_sale=[
[190.75572,189.94437,189.2342,190.34938,189.08669,188.56415,190.10588],
[28.618061,28.872581,28.873682,28.74203,28.776909,28.86964,28.661997],
[26.962397,27.805391,27.65219,27.210875,27.88252,27.764929,27.24346],
[32.29913,31.795496,32.425,31.649815,31.683603,31.381622,31.967655],
[102.14867,101.67641,102.01936,102.55891,102.07538,102.73362,102.2766],
[56.307552,57.39569,58.130955,56.80816,56.4629,57.186737,57.953686]]#预测销量
predit_omega=[0.7]*6 #折扣
predit_gama=[0.1283,0.1551,0.1365,0.0668,0.0924,0.0945] #耗损率
day=6
ini_pos=[0.6,0.4,0.3,0.6,0.9,0.6,150,30,20,20,100,60]
ini_pos=np.array(ini_pos)
n_particles=1000
n_dimensions=12
lower_bound=np.array([0.3,0.2,0.2,0.3,0.3,0.3,10,0,0,10,10,20])
upper_bound=np.array([1.2,0.98,0.8,0.9,1.2,1.0,450,90,75,60,300,250])
bounds=(lower_bound,upper_bound)
weight_ini=0.3 #给定的初始值占总粒子的比例
pos_given=np.random.uniform(
low=lower_bound,high=upper_bound,size=(int(n_particles*weight_ini),n_dimensions)
)
pos_given=0.8*pos_given+0.2*ini_pos
pos_given=np.clip(pos_given,lower_bound,upper_bound)
pos_random=np.random.uniform(
low=lower_bound, high=upper_bound, size=(int(n_particles * (1-weight_ini)), n_dimensions)
)
Initial_pos=np.vstack((pos_given,pos_random))
def modify(y,x,idx):
if idx==0:
sale=-1.91*x+1.03*y+3.82
if sale<0:
return 0
else:
return sale
if idx==1:
sale=-0.25*x+0.94*y+1.07
if sale < 0:
return 0
else:
return sale
if idx == 2:
sale = -0.03 * x + 0.97 * y -0.98
if sale < 0:
return 0
else:
return sale
if idx == 3:
sale = -0.07 * x + 0.99 * y +1.70
if sale < 0:
return 0
else:
return sale
if idx == 4:
sale = -0.029 * x + 0.99 * y +2.54
if sale < 0:
return 0
else:
return sale
if idx == 5:
sale = -0.51 * x + 1.01 * y + 2.32
if sale < 0:
return 0
else:
return sale
#目标函数
def Objective_function(x):
profit_list=[]
for i in range(n_particles):
x_new=x[i]
profit=0
for idx in range(6):
alpha=x_new[idx] #利润率
beta=x_new[idx+6] #进货量
buy=predit_buy[idx][day] #预测进价
sale_lstm=predit_sale[idx][day] #预测销量
omega=predit_omega[idx] #折扣
gama=predit_gama[idx] #耗损率
sale_price_normal=buy*(1+alpha) #好货的售价
sale_price_discount= buy * (1 + alpha) *omega # 差货的售价
good=beta*(1-gama) #好的进货量
bad=beta*gama #差的进货量
sale_modify=modify(sale_lstm,sale_price_normal,idx) #预测的销量
w1=sale_modify*(1-gama)*sale_price_normal+sale_modify*gama*sale_price_discount
w2=beta*buy
if beta<=sale_modify:
profit+=(w1-w2)-(sale_modify-beta)*20
else:
profit+=(w1-w2)
profit_list.append(-profit)
return profit_list
options={'c1':0.5,'c2':0.5,'w':0.6} #个人 社会 继承
optimizer=ps.single.GlobalBestPSO(n_particles=n_particles,dimensions=n_dimensions,options=options,bounds=bounds,init_pos=Initial_pos)
best_position,best_cost=optimizer.optimize(Objective_function,iters=300,verbose=True)
fig,ax=plt.subplots(figsize=(8,6))
fig, ax = plt.subplots(figsize=(8, 6))
ax.plot(optimizer.cost_history)
ax.set_title("目标函数(-利润)变化曲线")
ax.set_xlabel("迭代次数")
ax.set_ylabel("目标函数(-利润)")
plt.show()
plt.show()
def pre_sale_mount(x1,x2):
list_sale=[]
for idx in range(6):
x=x1[idx]
y=x2[idx]
sale=modify(y,x,idx)
list_sale.append(sale)
return list_sale
def count_how_much(x,y):
list_how_much=[]
for i in range(6):
list_how_much.append((x[i]+1)*y[i])
return list_how_much
print("最大利润:",-(best_position))
print("进货量(变量):",best_cost[-6:].tolist())
sale_price=count_how_much( best_cost[-6:].tolist(),[row[day] for row in predit_buy])
print("销售价格:",sale_price)
print("利润率(变量):",best_cost[:6].tolist())
print("修正后的销量:",pre_sale_mount(sale_price,[row[day] for row in predit_buy]))
import matplotlib.pyplot as plt
import numpy as np
import pyswarms as ps
plt.rcParams['font.sans-serif'] = [u'simHei']
plt.rcParams['axes.unicode_minus'] = False
list_1=['花叶类','花菜类','水生根茎类','茄类','辣椒类','食用菌']
predit_buy=[
[3.285864,3.2921748,3.2889733,3.285188,3.2851105,3.2964268,3.2876368] ,
[7.7414317,7.763459,7.814592,7.794937,7.747068,7.810813, 7.7633805] ,
[12.018661,11.912668,12.027704,11.941088,11.92054,12.118359,11.972251] ,
[4.5562034,4.601929,4.5483465,4.549116,4.532483,4.539543,4.601603],
[3.7067149,3.65774,3.6644902,3.6755412,3.6658049,3.6834998,3.6471841] ,
[4.015016,4.075036,4.025253,4.0783653,4.0397897,4.0211616,4.0569587] ]#预测进价
predit_sale=[
[190.75572,189.94437,189.2342,190.34938,189.08669,188.56415,190.10588],
[28.618061,28.872581,28.873682,28.74203,28.776909,28.86964,28.661997],
[26.962397,27.805391,27.65219,27.210875,27.88252,27.764929,27.24346],
[32.29913,31.795496,32.425,31.649815,31.683603,31.381622,31.967655],
[102.14867,101.67641,102.01936,102.55891,102.07538,102.73362,102.2766],
[56.307552,57.39569,58.130955,56.80816,56.4629,57.186737,57.953686]]#预测销量
predit_omega=[0.7]*6 #折扣
predit_gama=[0.1283,0.1551,0.1365,0.0668,0.0924,0.0945] #耗损率
day=6
ini_pos=[0.6,0.4,0.3,0.6,0.9,0.6,150,30,20,20,100,60]
ini_pos=np.array(ini_pos)
n_particles=1000
n_dimensions=12
lower_bound=np.array([0.3,0.2,0.2,0.3,0.3,0.3,10,0,0,10,10,20])
upper_bound=np.array([1.2,0.98,0.8,0.9,1.2,1.0,450,90,75,60,300,250])
bounds=(lower_bound,upper_bound)
weight_ini=0.3 #给定的初始值占总粒子的比例
pos_given=np.random.uniform(
low=lower_bound,high=upper_bound,size=(int(n_particles*weight_ini),n_dimensions)
)
导入相关的库:使用import
语句导入matplotlib.pyplot
库用于绘图,导入numpy
库用于进行数值计算,导入pyswarms
库用于粒子群优化算法的求解。
设置中文显示:通过设置plt.rcParams
的参数,将字体设置为中文,解决中文显示乱码的问题。
定义列表list_1
:该列表包含了一个分类的名称,如花叶类、花菜类等。
定义预测进价和预测销量:分别定义了predit_buy
和predit_sale
两个二维数组,表示了不同分类的商品的预测进价和预测销量。
定义折扣和耗损率:分别定义了predit_omega
和predit_gama
两个一维数组,表示了不同分类的商品的折扣和耗损率。
定义天数:使用day
变量定义了天数,表示了需要进行优化的天数。
初始化位置和粒子数:使用ini_pos
和n_particles
两个变量分别定义了初始位置和粒子数。
定义维度和边界:使用n_dimensions
、lower_bound
和upper_bound
定义了粒子的维度和边界。
定义粒子范围:使用bounds
变量定义了粒子的范围,其中bounds
是一个元组,包含了lower_bound
和upper_bound
。
初始化给定的初始值粒子:使用pos_given
变量随机生成了一些粒子,数量为n_particles*weight_ini
,这些粒子的位置在lower_bound
和upper_bound
之间。
总的来说,这段代码是为粒子群优化算法做准备的,包括定义了粒子的维度、边界、位置范围和初始位置等。?
def modify(y,x,idx):
if idx==0:
sale=-1.91*x+1.03*y+3.82
if sale<0:
return 0
else:
return sale
if idx==1:
sale=-0.25*x+0.94*y+1.07
if sale < 0:
return 0
else:
return sale
if idx == 2:
sale = -0.03 * x + 0.97 * y -0.98
if sale < 0:
return 0
else:
return sale
if idx == 3:
sale = -0.07 * x + 0.99 * y +1.70
if sale < 0:
return 0
else:
return sale
if idx == 4:
sale = -0.029 * x + 0.99 * y +2.54
if sale < 0:
return 0
else:
return sale
if idx == 5:
sale = -0.51 * x + 1.01 * y + 2.32
if sale < 0:
return 0
else:
return sale
这段代码定义了一个函数modify
,用于更新每个商品的售价。
该函数接收三个参数,分别是商品的当前进价x
、商品的当前销量y
和商品的索引idx
。
对于不同的商品索引,函数分别计算出相应的售价,并返回该售价。
具体地,该函数使用如下代码计算销量和进价对售价的影响关系:
- 对于索引为0的商品,使用
sale = -1.91*x + 1.03*y + 3.82
计算售价;- 对于索引为1的商品,使用
sale = -0.25*x + 0.94*y + 1.07
计算售价;- 对于索引为2的商品,使用
sale = -0.03*x + 0.97*y - 0.98
计算售价;- 对于索引为3的商品,使用
sale = -0.07*x + 0.99*y + 1.70
计算售价;- 对于索引为4的商品,使用
sale = -0.029*x + 0.99*y + 2.54
计算售价;- 对于索引为5的商品,使用
sale = -0.51*x + 1.01*y + 2.32
计算售价。
最后,对于计算出的售价,如果小于0,则将其设为0,并返回相应值;否则,直接返回计算出的售价。
def Objective_function(x):
profit_list=[]
for i in range(n_particles):
x_new=x[i]
profit=0
for idx in range(6):
alpha=x_new[idx] #利润率
beta=x_new[idx+6] #进货量
buy=predit_buy[idx][day] #预测进价
sale_lstm=predit_sale[idx][day] #预测销量
omega=predit_omega[idx] #折扣
gama=predit_gama[idx] #耗损率
sale_price_normal=buy*(1+alpha) #好货的售价
sale_price_discount= buy * (1 + alpha) *omega # 差货的售价
good=beta*(1-gama) #好的进货量
bad=beta*gama #差的进货量
sale_modify=modify(sale_lstm,sale_price_normal,idx) #预测的销量
w1=sale_modify*(1-gama)*sale_price_normal+sale_modify*gama*sale_price_discount
w2=beta*buy
if beta<=sale_modify:
profit+=(w1-w2)-(sale_modify-beta)*20
else:
profit+=(w1-w2)
profit_list.append(-profit)
return profit_list
这段代码定义了一个目标函数Objective_function
,用于计算粒子群优化算法中粒子的适应度值。
该函数接收一个参数x
,表示一组粒子的位置。
在函数中,首先创建了一个空列表profit_list
,用于存储每个粒子的适应度值。
然后,通过一个循环遍历粒子群中的每个粒子。对于每个粒子,按照指定的规则计算适应度值。
具体地,通过另一个循环遍历6个商品的索引。对于每个商品,获取粒子位置中的两个参数alpha
和beta
,分别表示利润率和进货量。
然后,根据预测的进价、预测的销量、折扣率和耗损率,计算出好货的售价、差货的售价、好的进货量和差的进货量。
接着,根据给定的修改函数modify
,计算出修改后的销量,并结合其他参数计算出权重值w1
和w2
。
最后,根据进货量和销量的大小关系,计算出相应的适应度值,并将其加入到profit_list
中。
最终,函数返回profit_list
,即每个粒子的适应度值列表。需要注意的是,这里返回的是适应度值的负值,因为粒子群算法通常是追求最小化问题,而此处是最大化问题。
options={'c1':0.5,'c2':0.5,'w':0.6} #个人 社会 继承
optimizer=ps.single.GlobalBestPSO(n_particles=n_particles,dimensions=n_dimensions,options=options,bounds=bounds,init_pos=Initial_pos)
best_position,best_cost=optimizer.optimize(Objective_function,iters=300,verbose=True)
fig,ax=plt.subplots(figsize=(8,6))
fig, ax = plt.subplots(figsize=(8, 6))
ax.plot(optimizer.cost_history)
ax.set_title("目标函数(-利润)变化曲线")
ax.set_xlabel("迭代次数")
ax.set_ylabel("目标函数(-利润)")
plt.show()
plt.show()
这段代码是使用Python的pyswarms库进行粒子群优化算法的实现。
其中,首先定义了一个字典options
,包含了三个参数,分别是个人认知系数c1
、社会认知系数c2
和惯性权重w
。这些参数用于调整粒子群算法的搜索策略。
然后,定义了一个GlobalBestPSO
对象optimizer
,用于实现全局最优PSO算法。该对象接收了以下参数:
n_particles
:粒子群的个数;dimensions
:搜索空间的维度;options
:一个字典,包含了调整算法参数的信息;bounds
:一个元组,表示搜索空间每个维度的上下限;init_pos
:初始位置。
接着,调用optimizer.optimize
函数,开始进行优化计算。该函数用于最大化目标函数,期望找到最优的粒子位置。该函数需要接收两个参数:
Objective_function
:目标函数,用于评估粒子位置的适应度;iters
:迭代次数。
verbose=True
表示在迭代过程中显示详细信息。
最后,通过matplotlib
库绘制收敛曲线,用于判断算法是否已经找到最优解。
ef pre_sale_mount(x1,x2):
list_sale=[]
for idx in range(6):
x=x1[idx]
y=x2[idx]
sale=modify(y,x,idx)
list_sale.append(sale)
return list_sale
def count_how_much(x,y):
list_how_much=[]
for i in range(6):
list_how_much.append((x[i]+1)*y[i])
return list_how_much
print("最大利润:",-(best_position))
print("进货量(变量):",best_cost[-6:].tolist())
sale_price=count_how_much( best_cost[-6:].tolist(),[row[day] for row in predit_buy])
print("销售价格:",sale_price)
print("利润率(变量):",best_cost[:6].tolist())
print("修正后的销量:",pre_sale_mount(sale_price,[row[day] for row in predit_buy]))
这段代码首先定义了两个函数pre_sale_mount
和count_how_much
,然后输出了一些结果。
函数pre_sale_mount
接受两个参数x1
和x2
,分别表示两个列表。通过一个循环遍历6个商品的索引,取出对应位置的元素x
和y
,并调用modify
函数进行处理得到销量,并将结果添加到list_sale
列表中。最终返回list_sale
列表作为函数的输出。
函数count_how_much
接受两个参数x
和y
,同样是两个列表。通过一个循环遍历6个商品的索引,将每个商品的(x[i]+1)*y[i]
的结果添加到list_how_much
列表中。最终返回list_how_much
列表作为函数的输出。
接下来,代码打印了一些结果:
"最大利润:"后面紧跟着
-(best_position)
,即全局最优粒子位置的最大化目标函数的负值。这里使用负值表示了最大化利润。"进货量(变量):"后面紧跟着
best_cost[-6:].tolist()
,即全局最优粒子位置中最后6个维度的取值,表示优化后的进货量。"销售价格:"后面紧跟着
sale_price
,即通过调用count_how_much
函数计算出销售价格。"利润率(变量):"后面紧跟着
best_cost[:6].tolist()
,即全局最优粒子位置中前6个维度的取值,表示优化后的利润率。"修正后的销量:"后面紧跟着
pre_sale_mount(sale_price,[row[day] for row in predit_buy])
,即将销售价格和预测的购买量传入pre_sale_mount
函数,计算修正后的销量。