（粒子群算法）的matlab及其python实现

粒子群算法的实现

1. 基本概念：

鸟群中有个体和群体，个体和群体的信息是可以互通的。个体在随机搜寻食物的过程中，只要跟踪离食物最近的群体，就能最有效地找到食物。

（1）粒子：优化问题的候选解，指鸟群中的独立个体；

（2）位置：候选解所在的位置，即鸟群个体的位置；

（3）速度：粒子的移动速度；

（4）适应度：评价粒子优劣的值，一般为优化目标函数的数值；

（5）个体极值：单个粒子迄今为止找到的最佳位置；

（6）群体极值：所有粒子迄今为止找到的最佳位置。

2.基本流程

（1）初始粒子；

（2）计算适应度值；

（在实际问题的解决中，构建目标函数是最重要的，也是最难的。）

（3）定义初始个体极值与群体极值；

（4）更新粒子位置与速度；（最经典的部分）

（5）更新个体极值和群体极值。

3、函数代码实现

（1）主函数：

（POSmain.m）

clear
clc

%绘制原图        图1目标函数的三维网格图
x1=-15:1:15;
x2=-15:1:15;
[x1,x2]=meshgrid(x1,x2);
y=x1.^2+x2.^2-x1.*x2-10.*x1-4.*x2+60;
mesh(x1,x2,y);    
hold on;

%%预设参数
n=100; %粒子群的规模
d=2; %变量个数
c1=2;
c2=2;
w=0.9;%权重一般设为0.9
K=50; %迭代次数

%%分布粒子的取值范围及速度的取值范围
x=-10+20*rand(n,d);  %产生一个n*d的随机矩阵,x中元素取值在[-10,10]中
v=-5+10*rand(n,d);   %产生一个n*d的随机矩阵,v中元素取值在[-5,5]中

%%计算适应度
fit=zeros(n,1);%创建一个n*1全为零的矩阵，实现每个个体适应度的初始化
for j=1:n
    fit(j)=fun(x(j,:));%x(j,:)代表了一个个体的位置。
end

pbest=x;%初始化每个个体的最好位置。
ind=find(min(fit)==fit);%find找到非零数的下标，此为找到最小适应度对应的下标。
gbest=x(ind,:);%初始化全体个体中的最好位置
%h=scatter3(x(:,1),x(:,2),fit,'o');  %图2 粒子的初始分布图，MATLAB三维散点图的绘制。

%%更新速度与位置
for i=1:K
    for j=1:n %每一个个体都得更新。
       v(j,:)=w*v(j,:) + c1*rand*(pbest(j,:)-x(j,:)) + c2*rand*(gbest-x(j,:));%rand是[0,1]随机数,速度更新公式
       
       v(j,find(v(j,:)<-5))=-5;%这里发现速度小于-5时取-5
       v(j,find(v(j,:)>5))=5;%这里发现速度大于5时取5
       
       x(j,:)=x(j,:)+0.5*v(j,:);%位置更新公式
       x(j,find(x(j,:)<-10))=-10;%这里发现位置小于-10时取-10
       x(j,find(x(j,:)>10))=10;%这里发现位置大于10时取10
       
       %重新计算适应度
       fit(j)=fun(x(j,:));
       if x(j,:)<fun(pbest(j,:))
           pbest(j,:)=x(j,:);
       end
       if fun(pbest(j,:))<fun(gbest)
           gbest=pbest(j,:);
       end
    end
    fitnessbest(i)=fun(gbest);
    pause(0.01);    %为了直观，每更新一次，暂停0.01秒
    h.XData=x(:,1);
    h.YData=x(:,2);
    h.ZData=fit;
end
 h=scatter3(x(:,1),x(:,2),fit,'o');  %图2 粒子的初始分布图，MATLAB三维散点图的绘制。
hold on;
figure(2)
plot(fitnessbest);
xlabel('迭代次数');

2、 目标函数：

 （fun.m）

function y=fun(x)

y=2*x(1)^2+x(2)^2-x(1)*x(2)-10*x(1)-4*x(2)+60;

（1）最后一次迭代的种群的图像

 

（2）迭代与适应度值的图像

4.python实现，适应度函数为f ( x , y ) = x^2 + y^2 + x ,求解f(x,y)的最小值点.

代码：

# 速度
# Vi+1 = w*Vi + c1 * r1 * (pbest_i - Xi) + c2 * r2 * (gbest_i - Xi)
# 位置
# Xi+1 = Xi + Vi+1
# vi, xi 分别表示粒子第i维的速度和位置
# pbest_i, gbest_i 分别表示某个粒子最好位置第i维的值、整个种群最好位置第i维的值

import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl

mpl.rcParams['font.sans-serif'] = ['SimHei']  # 指定默认字体
mpl.rcParams['axes.unicode_minus'] = False  # 解决保存图像是负号'-'显示为方块的问题

def fitness_func(X):
    """计算粒子的的适应度值，也就是目标函数值，X 的维度是 size * 2 """
    A = 10
    pi = np.pi
    x = X[:, 0]
    y = X[:, 1]
    return x**2+y**2+x

#更新速度

def velocity_update(V, X, pbest, gbest, c1, c2, w, max_val):
    """
    根据速度更新公式更新每个粒子的速度，有20个体，维度为2。
    :param V: 粒子当前的速度矩阵，20*2 的矩阵
    :param X: 粒子当前的位置矩阵，20*2 的矩阵
    :param pbest: 每个粒子历史最优位置，20*2 的矩阵
    :param gbest: 种群历史最优位置，1*2 的矩阵
    """
    size = X.shape[0]
    r1 = np.random.random((size, 1))
    r2 = np.random.random((size, 1))
    V = w*V+c1*r1*(pbest-X)+c2*r2*(gbest-X)
    # 防止越界处理
    V[V < -max_val] = -max_val
    V[V > max_val] = max_val
    return V

#更新粒子位置

def position_update(X, V):
    """
    根据公式更新粒子的位置
    :param X: 粒子当前的位置矩阵，维度是 20*2
    :param V: 粒子当前的速度举着，维度是 20*2
    """
    return X+V

#粒子群算法的主要流程

def pos():
    w = 1
    c1 = 2
    c2 = 2
    r1 = None
    r2 = None
    dim = 2
    size = 20
    iter_num = 1000 #迭代参数
    max_val = 0.5
    best_fitness = float(9e10)#初始化最优适应度值
    fitness_val_list = []#初始化适应度值
    # 初始化种群各个粒子的位置
    X = np.random.uniform(-5, 5, size=(size, dim))
    # 初始化各个粒子的速度
    V = np.random.uniform(-0.5, 0.5, size=(size, dim))
    # print(X)
    p_fitness = fitness_func(X)#初始化的个体的最优适应值
    g_fitness = p_fitness.min()#初始化全体个体的最优适应值
    fitness_val_list.append(g_fitness)

    # 初始化的个体最优位置和种群最优位置
    pbest = X
    gbest = X[p_fitness.argmin()]
    # 迭代计算
    for i in range(1, iter_num):
        V = velocity_update(V, X, pbest, gbest, c1, c2, w, max_val)
        X = position_update(X, V)
        p_fitness2 = fitness_func(X)
        g_fitness2 = p_fitness2.min()

        # 更新每个粒子的历史最优位置
        for j in range(size):
            if p_fitness[j] > p_fitness2[j]:
                pbest[j] = X[j]
                p_fitness[j] = p_fitness2[j]
            # 更新群体的最优位置
            if g_fitness > g_fitness2:
                gbest = X[p_fitness2.argmin()]
                g_fitness = g_fitness2
            # 记录最优迭代记录
            fitness_val_list.append(g_fitness)
            i += 1

    # 输出迭代结果
    print("最优值是：%.5f" % fitness_val_list[-1])
    print("最优解是：x=%.5f,y=%.5f" % (gbest[0], gbest[1]))

    # 绘图
    plt.plot(fitness_val_list, color='r')
    plt.title('迭代过程')
    plt.show()

pos()

运行结果：

5.总结
作为自学使用，参考两位大佬，万分感谢

Matlab粒子群算法（PSO）优化程序——经典实例_箫韵的博客-CSDN博客_粒子群优化算法matlab程序

粒子群优化算法（PSO）简介及MATLAB实现_Yancy的博客-CSDN博客_粒子群优化算法matlab程序

粒子群算法Python代码实现_大灰狼学编程的博客-CSDN博客_粒子群算法python代码