粒子群算法 模拟退火算法

粒子群算法

(求解二元函数极值)

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 = 'max' #选择最大值还是最小值
x_limit = [-4,4]
y_limit = [-4,4]
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()

 

模拟退火算法

(解二元函数极值)

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

fun = '5 * sin(x * y) + x ** 2 + y ** 2'   #目标函数
choice = 'max'
x_limit = [-4,4]  #x取值范围
y_limit = [-4,4]  #y取值范围

T = 100  #初始温度为100
alpha = 0.95 #温度衰减系数
iter_time = 200 #最大迭代次数
max_repeat_time = 30 #最大重复次数
repeat_time = 0 #重复次数


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

def init_z(z = None,first = False): #随机产生新解,如果为初始值,则不需要在z周围生成
    if first:
        z = np.random.rand(2,)
        z[0] = z[0] * (x_limit[1] - x_limit[0]) + x_limit[0]
        z[1] = z[1] * (y_limit[1] - y_limit[0]) + y_limit[0]
        return z

    else:
        new_z = np.zeros((2,))
        new_z[0] = z[0] + (np.random.rand() - 0.5) * (x_limit[1] - x_limit[0]) * 0.1  #在z周围随机生成新解
        new_z[1] = z[1] + (np.random.rand() - 0.5) * (y_limit[1] - y_limit[0]) * 0.1  #在z周围随机生成新解

        while ((new_z[0] < x_limit[0]) or (new_z[0] > x_limit[1])):
            if new_z[0] < x_limit[0]:  #边缘检测
                r = np.random.rand()
                new_z[0] = r * x_limit[0] + (1 - r) * z[0]

            if new_z[0] > x_limit[0]:
                r = np.random.rand()
                new_z[0] = r * x_limit[1] + (1 - r) * z[1]

        while ((new_z[1] < y_limit[0]) or (new_z[1] > y_limit[1])):
            if new_z[1] < y_limit[0]:  #边缘检测
                r = np.random.rand()
                new_z[1] = r * y_limit[1] + (1 - r) * new_z[1]

            if new_z[1] > y_limit[0]:
                r = np.random.rand()
                new_z[1] = r * y_limit[0] + (1 - r) * new_z[1]

        return new_z


def update(z,T,count = 0):  #更新z

    if count >= iter_time:  #迭代次数达到最大,返回当前值
        return z

    new_z = init_z(z)
    f_new_z = f(fun,z[0],z[1])
    f_z = f(fun,z[0],z[1])

    if choice == 'max':
        if f_new_z < f_z:
            p = np.exp(-np.abs(f_new_z - f_z) / T)  #接受概率
            r = np.random.rand()
            if r < p:
                return new_z
            else:
                update(z,T,count = count + 1)  #不接受则重复操作
        else:
            return new_z
    if choice == 'min':
        if f_new_z > f_z:
            p = np.exp(-np.abs(f_new_z - f_z) / T)  #接受概率
            r = np.random.rand()
            if r < p:
                return new_z
            else:
                update(z,T,count = count + 1)

        else:
            return new_z


z = init_z(first = True)
global_best = z  #全局最优解

while T >= 0.000001:
    for i in range(100):  #为了保证搜索彻底,每个温度下都重复搜索多次
        z = update(z,T)  #新解产生
        f_z = f(fun,z[0],z[1])
        if (f_z - 1e-5) <= f(fun,global_best[0],global_best[1]) <= (f_z + 1e-5):  #浮点数精度问题
            repeat_time += 1
        else:
            repeat_time = 0
        if repeat_time >= max_repeat_time:  #连续多次找到同一个最优解
            break

        if choice == 'max':
            if f_z > f(fun,global_best[0],global_best[1]):
                global_best = z
        if choice == 'min':
            if f_z < f(fun,global_best[0],global_best[1]):
                global_best = z

    T = T * 0.95  #温度降低

print("最优解为: x = ",global_best[0],"y = ",global_best[1],"z = ",f(fun,global_best[0],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)
f_z = f(fun,x,y)
ax.plot_surface(x, y, f_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='blue')

plt.show()