使用基础粒子群(PSO)算法求解一元及二元方程的Python代码

最近在看清风的数学建模视频，其中有两道题：

求一元函数的最值问题

1. 题目
   求函数 y = 11 s i n ( x ) + 7 c o s ( 5 x ) y =11sin(x)+7cos(5x) y=11sin(x)+7cos(5x)在 x ∈ [ − 3 , 3 ] x∈[-3,3] x∈[−3,3]内的最大值。
2. 流程图
   
3. 代码实现

# 第一步，绘制函数图像
import numpy as np
import matplotlib.pyplot as plt

def func(x):
    return 11*np.sin(x)+7*np.cos(5*x)

x0 = np.linspace(-3,3,1000)
y0 = func(x0)
fig,ax = plt.subplots()
ax.plot(x0,y0)
ax.set_title('y = 11sin(x)+7cos(5x)')

# 第二步，设置粒子群算法的参数
n = 10 # 粒子数量
narvs = 1 # 变量个数
c1 = 2 # 个体学习因子
c2 = 2 # 社会学习因子
w = 0.9 # 惯性权重
K = 50 # 迭代次数
vmax = (3-(-3))*0.2 # 粒子的最大速度
x_lb = -3 # x的下界
x_ub = 3 # x的上界

# 第三步，初始化粒子
x = x_lb + (x_ub-x_lb)*np.random.rand(n)
v = -vmax + 2*vmax*np.random.rand(n)

# 第四步，计算适应度
fit = func(x) # 计算每个粒子的适应度
pbest = x # 初始化这n个例子迄今为止找到的最佳位置
ind = np.argmax(fit) # 找到适应度最大的那个粒子的下标
gbest = x[ind]
gbest_total = np.zeros(K)
    
# 第五步，更新粒子速度和位置
for j in range(K): # 外层循环，共K次
    v = w*v + c1*np.random.rand(n)*(pbest-x) + c2*np.random.rand(n)*(gbest-x)
    v[np.where(v<-vmax)] = -vmax # 速度小于-vmax的元素赋值为-vmax
    v[np.where(v>vmax)] = vmax # 速度大于vmax的元素赋值为vmax
    x = x + v # 更新第i个粒子的位置
    x[np.where(x<x_lb)] = x_lb
    x[np.where(x>x_ub)] = x_ub
    
    # 第六步，重新计算适应度并找到最优粒子
    fit = func(x) # 重新计算n个粒子的适应度
    for k in range(n): # 更新第k个粒子迄今为止找到的最佳位置
        if fit[k]>func(pbest[k]):
            pbest[k] = x[k]
    if np.max(fit)>func(gbest): # 更新所有粒子迄今找到的最佳位置
        gbest = x[np.argmax(fit)]
    gbest_total[j] = func(gbest)

ax.scatter(x,fit,c='r',marker='x')
ax.scatter(gbest,func(gbest),c='r')
ax.plot(gbest*np.ones(10),np.linspace(-20,func(gbest),10),'--')
ax.plot(np.linspace(-3,gbest,10),func(gbest)*np.ones(10),'--')
ax.set_ylim(-20,20)
ax.autoscale(axis='x',tight=True)
fig,ax = plt.subplots()
ax.plot(np.arange(K),gbest_total)
plt.rcParams['font.sans-serif']=['SimHei'] #用来正常显示中文标签
ax.set_title('算法找到的最大值随优化次数的变化曲线')
ax.set_xlabel('循环次数')
ax.autoscale(axis='x',tight=True)
print('找到的最优解为：',func(gbest))    

4. 运算结果

• 打印输出最优解：找到的最优解为： 17.492618154736594
• 显示最终粒子所在的位置：
  
  其中红色叉号表示最后一次循环结束后粒子所在的位置；红色圆圈(用橙色和绿色的线标识横纵坐标值)表示算法寻到的最优解。
• 循环过程中最优解的变化曲线：
  

求二元函数的最值问题

1. 题目
   求函数 y = x 1 2 + x 2 2 − x 1 x 2 − 10 x 1 − 4 x 2 + 60 y =x_1^2+x_2^2-x_1x_2-10x_1-4x_2+60 y=x12​+x22​−x1​x2​−10x1​−4x2​+60在 x 1 , x 2 ∈ [ − 15 , 15 ] x_1,x_2∈[-15,15] x1​,x2​∈[−15,15]内的最小值。
2. 流程图
   
3. 代码实现

# 第一步，绘制函数图像
import numpy as np
import matplotlib.pyplot as plt
import mpl_toolkits.mplot3d

def func(x,y):
    return x**2+y**2-x*y-10*x-4*y+60

x0 = np.linspace(-15,15,100)
y0 = np.linspace(-15,15,100)
x0,y0 = np.meshgrid(x0,y0)
z0 = func(x0,y0)
fig = plt.figure(constrained_layout=True)
ax = fig.add_subplot(projection='3d')
ax.plot_surface(x0,y0,z0,cmap=plt.cm.viridis,alpha=0.7)
#ax.plot_wireframe(x0,y0,z0) # 另一种绘图方式
ax.set_title('$y = x_1^2+x_2^2-x_1x_2-10x_1-4x_2+60$')

# 第二步，设置粒子群算法的参数
n = 30 # 粒子数量
narvs = 2 # 变量个数
c1 = 2 # 个体学习因子
c2 = 2 # 社会学习因子
w = 0.9 # 惯性权重
K = 40 # 迭代次数
vxmax = np.array([(15-(-15))*0.2,(15-(-15))*0.2]) # 粒子在x方向的最大速度
x_lb = np.array([-15,-15]) # x和y的下界
x_ub = np.array([15,15]) # x和y的上界

# 第三步，初始化粒子   
x = x_lb + (x_ub-x_lb)*np.random.rand(n,narvs)
v = -vxmax + 2*vxmax*np.random.rand(n,narvs)

# 第四步，计算适应度
fit = func(x[:,0],x[:,1]) # 计算每个粒子的适应度
pbest = x # 初始化这n个例子迄今为止找到的最佳位置
ind = np.argmax(fit) # 找到适应度最大的那个粒子的下标
gbest = x[ind,:]
gbest_total = np.zeros(K)

# 第五步，更新粒子速度和位置
for j in range(K): # 外层循环，共K次
    for p in range(n):
        v[p,:] = w*v[p,:] + c1*np.random.rand(narvs)*(pbest[p,:]-x[p,:]) + c2*np.random.rand(narvs)*(gbest-x[p,:])
    loc_v = np.where(v<-vxmax)
    v[loc_v] = -vxmax[loc_v[1]] # 速度小于-vmax的元素赋值为-vmax
    loc_v = np.where(v>vxmax)
    v[loc_v] = vxmax[loc_v[1]] # 速度大于vmax的元素赋值为vmax
    x = x + v # 更新第i个粒子的位置
    loc_x = np.where(x<x_lb)
    x[loc_x] = x_lb[loc_x[1]]
    loc_x = np.where(x>x_ub)
    x[loc_x] = x_ub[loc_x[1]]
    
    # 第六步，重新计算适应度并找到最优粒子
    fit = func(x[:,0],x[:,1]) # 重新计算n个粒子的适应度
    for k in range(n): # 更新第k个粒子迄今为止找到的最佳位置
        if fit[k]<func(pbest[k,0],pbest[k,1]):
            pbest[k,:] = x[k,:]
    if np.min(fit)<func(gbest[0],gbest[1]): # 更新所有粒子迄今找到的最佳位置
        gbest = x[np.argmin(fit),:]
    gbest_total[j] = func(gbest[0],gbest[1])

ax.scatter(x[:,0],x[:,1],fit,c='r',marker='x')
fig,ax = plt.subplots()
ax.plot(np.arange(K),gbest_total)
plt.rcParams['font.sans-serif']=['SimHei'] #用来正常显示中文标签
plt.rcParams['axes.unicode_minus'] = False #用来显示负号
ax.set_title('算法找到的最小值随优化次数的变化曲线')
ax.set_xlabel('循环次数')
ax.autoscale(axis='x',tight=True)
print('找到的最优解为：',func(gbest[0],gbest[1]))

4. 运算结果

• 打印输出最优解：找到的最优解为： 8.031080661489156
• 显示最终粒子所在的位置：
  

图中红色叉号表示最后一次循环结束后粒子所在的位置。

• 循环过程中最优解的变化曲线：
  
  可以看出，当循环次数高于50之后（我最开始设置的循环次数为100），最优解也不会改变，因此后续是无效循环，可以据此适当调整循环次数。