粒子群算法

import numpy as np
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from numpy import *

w = 0.5  #惯性因子[0,1]
c1 = 2  #学习因子
c2 = 2  #学习因子
partical_num = 100  #粒子数量
fun = '5 * sin(x * y) + x ** 2 + y ** 2'   #目标函数
choice = 'min' #选择最大值还是最小值
x_limit = [-5,5]
y_limit = [-5,5]
vx = [-5,5]  #x方向速度范围
vy = [-5,5]  #y方向速度范围
iter_time = 50 #迭代次数

global_best = np.zeros((2,))  #全局最优
local_best = np.zeros((2,partical_num))  #每个粒子的历史最优

def f(fun,x,y):  #目标函数
    return eval(fun)

def init_pos_and_velocity(partical_num):  #传入粒子总数,随机初始化位置和速度
    pos = np.random.rand(2,partical_num)
    pos[0,:] = pos[0,:] * (x_limit[1] - x_limit[0]) + x_limit[0]
    pos[1,:] = pos[1,:] * (y_limit[1] - y_limit[0]) + y_limit[0]

    velo = np.random.rand(2, partical_num)
    velo[0, :] = velo[0, :] * (vx[1] - vx[0]) + vx[0]
    velo[1, :] = velo[1, :] * (vy[1] - vy[0]) + vy[0]

    return pos,velo

def update_pos_and_velo(pos,velo):  #更新位置和速度
    #更新速度
    for i in range(partical_num):
        delta_vx = c1 * np.random.rand() * (local_best[0,i] - pos[0,i]) + c2 * np.random.rand() * (global_best[0] - pos[0,i])
        delta_vy = c1 * np.random.rand() * (local_best[1,i] - pos[1,i]) + c2 * np.random.rand() * (global_best[1] - pos[1,i])
        velo[0,i] = (w * velo[0,i] + delta_vx) if (vx[0] < (w * velo[0,i] + delta_vx) < vx[1]) else np.random.rand() * (vx[1] - vx[0]) + vx[0] #越界检测
        velo[1,i] = (w * velo[1,i] + delta_vy) if (vy[0] < (w * velo[1,i] + delta_vy) < vy[1]) else np.random.rand() * (vy[1] - vy[0]) + vy[0]

    for i in range(partical_num):
        pos[0,i] = (pos[0,i] + velo[0,i]) if x_limit[0] < (pos[0,i] + velo[0,i]) < x_limit[1] else (x_limit[1] - x_limit[0]) + x_limit[0] #越界检测
        pos[1,i] = (pos[1,i] + velo[1,i]) if y_limit[0] < (pos[1,i] + velo[1,i]) < y_limit[1] else (y_limit[1] - y_limit[0]) + y_limit[0]


    return pos,velo

def update_global_best(pos): #更新全局最优
    global global_best
    z = f(fun,pos[0,:],pos[1,:])
    if choice == 'max':
        if np.max(z) > f(fun,global_best[0],global_best[1]):
            global_best = pos[:,np.argmax(z)]
    else:
        if np.min(z) < f(fun, global_best[0], global_best[1]):
            global_best = pos[:, np.argmin(z)]

def update_local_best(pos,local_best):  #更新历史最优
    if choice == 'max':
        index = f(fun,pos[0,:],pos[1,:]) > f(fun,local_best[0,:],local_best[1,:])
        local_best[:,index] = pos[:,index]
    else:
        index = f(fun, pos[0, :], pos[1,:]) < f(fun, local_best[0,:], local_best[1,:])
        local_best[:, index] = pos[:, index]


pos,velo = init_pos_and_velocity(partical_num)  #初始化粒子群

update_local_best(pos,local_best)
update_global_best(pos)

for i in range(iter_time):
    pos,velo = update_pos_and_velo(pos,velo)
    update_local_best(pos,local_best)
    update_global_best(pos)

print('极值为: x = ',global_best[0],'y = ',global_best[1])

figure = plt.figure(figsize = (10,8),dpi = 80)
# 新建一个3d绘图对象
ax = Axes3D(figure)
# 定义x,y 轴名称
plt.xlabel("x")
plt.ylabel("y")
x = np.arange(x_limit[0],x_limit[1],(x_limit[1] - x_limit[0]) / 50)
y = np.arange(y_limit[0],y_limit[1],(y_limit[1] - y_limit[0]) / 50)
x,y = np.meshgrid(x,y)
z = f(fun,x,y)
ax.plot_surface(x,y,z, rstride=1, cstride=1, color = 'green',alpha = 0.5)
ax.scatter(global_best[0],global_best[1],f(fun,global_best[0],global_best[1]),color = 'red')
plt.show()