Python学习9

大家好，这里是七七，本次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)
)

1. 导入相关的库：使用import语句导入matplotlib.pyplot库用于绘图，导入numpy库用于进行数值计算，导入pyswarms库用于粒子群优化算法的求解。

2. 设置中文显示：通过设置plt.rcParams的参数，将字体设置为中文，解决中文显示乱码的问题。

3. 定义列表list_1：该列表包含了一个分类的名称，如花叶类、花菜类等。

4. 定义预测进价和预测销量：分别定义了predit_buy和predit_sale两个二维数组，表示了不同分类的商品的预测进价和预测销量。

5. 定义折扣和耗损率：分别定义了predit_omega和predit_gama两个一维数组，表示了不同分类的商品的折扣和耗损率。

6. 定义天数：使用day变量定义了天数，表示了需要进行优化的天数。

7. 初始化位置和粒子数：使用ini_pos和n_particles两个变量分别定义了初始位置和粒子数。

8. 定义维度和边界：使用n_dimensions、lower_bound和upper_bound定义了粒子的维度和边界。

9. 定义粒子范围：使用bounds变量定义了粒子的范围，其中bounds是一个元组，包含了lower_bound和upper_bound。

10. 初始化给定的初始值粒子：使用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函数，计算修正后的销量。