大家好,这里是七七,本次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]))
代码1
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之间。
总的来说,这段代码是为粒子群优化算法做准备的,包括定义了粒子的维度、边界、位置范围和初始位置等。
代码2
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,并返回相应值;否则,直接返回计算出的售价。
代码3
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,即每个粒子的适应度值列表。需要注意的是,这里返回的是适应度值的负值,因为粒子群算法通常是追求最小化问题,而此处是最大化问题。
代码4
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
库绘制收敛曲线,用于判断算法是否已经找到最优解。
代码5
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函数,计算修正后的销量。