粒子群 PSO matlab

% 参数初始化  赋予随机性能够寻解和防止局部最优解   让随机数据收敛能够找到最优解  
clear
clc
close all
%函数
figure(1);
lbx=-5;ubx=5; %函数自变量x范围【-5,5】
lby=-5;uby=5; %函数自变量y范围【-5,5】
ezmesh('x^2-2*sin(2*pi*x)+y^2-2*cos(2*pi*y)-15',[lbx,ubx,lby,uby],100);

% 参数初始化
tic %运行时间开始
p = 10;%解的个数(粒子个数)
n =2;%解的维度
E = 1;%期望解值
x = rand(p,n)*10-5;%初始解集 p*n
xmax = 5;%解值的上限
xmin = -5;%解值的下限
v = rand(p,n);%初始速度，速度理解为每次解变化范围 （粒子可变位置范围）
vmax = 2;%速度上限，解值变化上限
vmin = -2;%速度下限，解值变化下限
w = 0.6;%速度的惯性因子
c1 =2;%个体历史最优权重
c2 =2;%整体历史最优权重
maxnum = 100;%迭代次数
k=1;%记录迭代次数
fitness = inline('x(1)^2-2*sin(2*pi*x(1))+x(2)^2-2*cos(2*pi*x(2))-15','x');%定义适应度函数（目标函数，求最小值）
% 初始化
% 获取个体数据，及留存个体历史数据
for i=1:p
    for j =1:n %判断是否存在超出解值范围的解
        if x(i,j)>xmax
            x(i,j) =xmax;
        elseif x(i,j)<xmin
            x(i,:) =xmin;
        end
    end
    f(i) =fitness(x(i,:));%获取每个解在目标函数中对应的结果（适应度值）
end
hold on
plot3([x(:,1)],[x(:,2)],[f],'r*');%画出当前解的分布
pf = f; %保留每一个解对应的求解结果（适应度）
px = x; %保留每一个解
ff(1) = min(pf)%记录初始解中的最小值
% 迭代 
while (k<maxnum) % 判断是否进行迭代
    for i=1:p 
        for j =1:n %判断是否存在超出解值范围的解
            if x(i,j)>xmax
                x(i,j) =xmax;
            elseif x(i,j)<xmin
                x(i,:) =xmin;
            end
        end
        f(i) =fitness(x(i,:));%获取每个解在目标函数中对应的结果（适应度值）
        if f(i) < pf(i) %判断解更新后，是否出现了求解结果（适应度）最小值
             pf(i) = f(i);%记录该最小值对应解的求解结果（适应度，新的最小值）
             px(i,:) = x(i,:);%记录该最小值对应解
        end
    end

    [gf,t] = min(pf);%获取所有个体中历史求解结果最小值及对应解
    gx = px(t,:);%及对应解
    
    %进行随机 和 保持随机下能够收敛  累加v 保证随机性 选适应度最小值保证收敛 
    for i=1:p
        %根据该解上一次最小值对应的解，及历史所有解中具有最小值的解
        %和解的可变化范围，得到该解的变化量
        v(i,:) = w*v(i,:)+c1*rand*(px(i,:)-x(i,:))+c2*rand*(gx-x(i,:));
        for j=1:n%判断该解的变化量是否超出范围
            if v(i,j) > vmax
                v(i,j) = vmax;
            elseif  v(i,j) <vmin
                v(i,j) = vmin;
            end
        end
        x(i,:) = x(i,:)+v(i,:);%该解累加对应的变化量，得到新解 
    end
    
    h = get(gca, 'children'); % 获取坐标轴的children属性
    delete(h(1)); % 删除上次绘制的子图形
    hold on
    plot3([x(:,1)],[x(:,2)],[pf],'r*');%绘制新的解集
    pause(0.3);
    %     for i=1:p
    %         plot3([x(p-i+1,1)],[x(p-i+1,2)],[pf(p-i+1)],'r*');
    %     end

    k=k+1;% 迭代次数累加
    best_x = gx;% 记录目前最小结果对应的解值
    best_f = gf;% 记录目前的最小结果（最小适应度值）
    ff(k) = gf;%对每次最小结果（最小适应度值）数据进行保留
end

best_f %最终的最小结果（最小适应度值）
figure;
plot(ff)% 每次迭代中最小结果（最小适应度值）数据变化曲线
toc; %记录运行时间