粒子群算法PSO

clc;
clear;
close all;
 
f= @(x) - (x - 10) .^ 2 + x .* sin(x) .* cos(2 * x) - 5 * x .* sin(3 * x) ; % 适应度函数表达式(求这个函数的最大值)  
figure(1);
fplot(f, [0 20], 'b-');                 % 画出初始图像 
 
d = 1;                           % 空间维数(该例子是1维)  
N = 15;                          % 初始种群个数         
x_limit = [0, 20];               % 设置位置限制  
v_limit = [-1, 1];               % 设置速度限制    
 
x = x_limit(1) + (x_limit(2) - x_limit(1)) * rand(N, d); %初始每个粒子的位置  
v = rand(N, d);                  % 初始每个粒子的速度  
 
pbest = x;                       % 初始化每个个体的历史最佳位置 
gbest = zeros(1, d);             % 初始化种群的历史最佳位置 ,创建全0数组 
 
fp_best = zeros(N, 1);           % 初始化每个个体的历史最佳适应度 为 0  
fg_best = -inf;                  % 初始化种群历史最佳适应度 为 负无穷  
 
iter = 50;                       % 最大迭代次数 
 
c1 = 1.6;                        % 自我学习因子  
c2 = 1.8;                        % 群体学习因子 
w=0.4;                           % 惯性学习因子
hold on;
plot(x, f(x), 'ro');
title('初始状态图');  
 
record = zeros(iter, 1);        % 记录器(用于记录 fg_best 的变化过程)
figure(2);
i = 1; 
while i <= iter  
    fx = f(x) ;                 % 计算个体当前适应度 
    % 个体寻优
    for j = 1:N        
        if fp_best(j) < fx(j)   % 第j个个体的历史（最优）适应度小于当前的适应度则更新
            fp_best(j) = fx(j); % 更新个体历史最佳适应度
            pbest(j) = x(j);    % 更新个体历史最佳位置 
        end   
    end
    % 群体寻优
    if fg_best < max(fp_best)   % 第i次迭代，种群历史最优适应值，小于当前迭代次数的群体最优适应值则更新
        [fg_best,ind_max] = max(fp_best);     % 更新群体历史最佳适应度
        gbest = pbest(ind_max);               % 更新群体历史最佳位置，其中 ind_max 是最大值所在的下标
        % 注：将上面的两个式子换成下面两个是不可以的
        % [gbest, ind_max] = max(pbest);      % 更新群体历史最佳位置，其中 ind_max 是最大值所在的下标
        % fg_best = fp_best(ind_max);         % 更新群体历史最佳适应度   
    end  
    
    v = v * w + c1 * rand() * (pbest - x) + c2 * rand() * (repmat(gbest, N, 1) - x); % 速度更新
    % 注： repmat(A,r1,r2):可将矩阵 扩充 为每个单位为A的r1*r2的矩阵
    
    % 边界速度处理  
    v(v > v_limit(2)) = v_limit(2);
    v(v < v_limit(1)) = v_limit(1);  
    
    x = x + v;% 位置更新  
    
    % 边界位置处理  
    x(x > x_limit(2)) = x_limit(2);  
    x(x < x_limit(1)) = x_limit(1);  
    
    record(i) = fg_best;%最大值记录  
    
    % 画动态展示图
    zuo_x = 0 : 0.01 : 20;  
    plot(zuo_x, f(zuo_x), 'b-', x, f(x), 'ro');
    title('状态位置变化')  
    pause(0.1)  
    i = i + 1;
%     if mod(i,10) == 0   % 显示进度
%         i
%     end
end  
figure(3);
plot(record);
title('收敛过程'); 
zuo_x = 0 : 0.01 : 20;  
 
figure(4);
plot(zuo_x, f(zuo_x), 'b-', x, f(x), 'ro');
title('最终状态图')
 
disp(['最佳适应度：',num2str(fg_best)]);  
disp(['最佳粒子的位置x：',num2str(gbest)]);

% 主调用函数，以下新建一个文件，命名为f_xy
% function z = f_xy(x, y)
% % 适应度函数
% z = - x.^2 - y.^2 + 10*cos(2*pi*x) + 10*cos(2*pi*y) + 100;
% end
 
clc;
clear;
close all;
 
[ZuoBiao_x, ZuoBiao_y] = meshgrid(-10:0.1:10,-5:0.1:5);   % 产生二维坐标
ZuoBiao_z = f_xy(ZuoBiao_x, ZuoBiao_y);
figure(1);
s = mesh(ZuoBiao_x, ZuoBiao_y, ZuoBiao_z);      % 画网格曲面图
s.FaceColor = 'flat';    % 修改网格图的属性
 
% 初始化种群 
N = 100;                        % 初始种群个数  
D = 2;                          % 空间维数  
iter = 50;                      % 迭代次数       
x_limit = [-10, 10; -5, 5];     % 位置限制  
v_limit = [ -2,  2; -1, 1];     % 速度限制  
 
x = zeros(N, D);
for i = 1:D 
    x(:,i) = x_limit(i, 1) + (x_limit(i, 2) - x_limit(i, 1)) * rand(N, 1);%初始种群的位置  
end  
v(:,1) = rands(N, 1) * 1;       % 初始种群x方向的速度 
v(:,2) = rands(N, 1) * 2;       % 初始种群y方向的速度 
 
p_best = x;                     % (初始化)每个个体的历史最佳位置  
f_best = zeros(1, D);           % (初始化)种群的历史最佳位置  
 
fp_best = zeros(N, 1) - inf;    % (初始化)每个个体的历史最佳适应度为负无穷  
fg_best = -inf;                 % (初始化)种群历史最佳适应度为负无穷
 
w = 0.8;                        % 惯性权重
c1 = 0.5;                       % 自我学习因子  
c2 = 0.5;                       % 群体学习因子 
 
figure(2);
s = mesh(ZuoBiao_x, ZuoBiao_y, ZuoBiao_z);    % 画网格曲面图
s.FaceColor = 'flat';                         % 修改网格图的属性
hold on;
plot3(x(:,1), x(:,2),f_xy(x(:,1), x(:,2)), 'ro');
title('初始状态图');  
 
figure(3); 
i = 1;  
record = zeros(iter, 1);                 % 记录器  
while i <= iter  
    fx = f_xy(x(:,1), x(:,2));       % 个体当前适应度     
    for j = 1:N        
        if fp_best(j) < fx(j)        % 记录最大值
            fp_best(j) = fx(j);      % 更新个体历史最佳适应度  
            p_best(j,:) = x(j,:);    % 更新个体历史最佳位置  
        end   
    end  
    if fg_best < max(fp_best)  
        [fg_best, ind_max] = max(fp_best);    % 更新群体历史最佳适应度  
        f_best = p_best(ind_max, :);          % 更新群体历史最佳位置  
    end  
    
    v = v * w + c1 * rand() * (p_best - x) + c2 * rand() * (repmat(f_best, N, 1) - x); % 速度更新
    
    % 速度处理
    for t=1:N
        for k=1:D
            if v(t,k) > v_limit(k,2)      % 超速处理
                v(t,k) = v_limit(k,2);
            elseif v(t,k) < v_limit(k,1)  % 慢速处理
                v(t,k) = v_limit(k,1);
            end   
        end
    end
 
    x = x + v;            % 位置更新
 
    % 边界处理
    for t=1:N
        for k=1:D
            if x(t,k) > x_limit(k,2)      % 超过边界上限
                x(t,k) = x_limit(k,2);
            elseif x(t,k) < x_limit(k,1)  % 超过边界下限
                x(t,k) = x_limit(k,1);
            end
        end
    end
    
    record(i) = fg_best;   % 最大值记录  
    i = i + 1;
    if mod(i,10)  == 0
        i                  % 收敛进度输出
    end
    pause(0.1) 
    
end
 
figure(4)
plot(record);
title('收敛过程');
 
% 画最终状态图
figure(5);
[ZuoBiao_x, ZuoBiao_y] = meshgrid(-10:0.1:10,-5:0.1:5);   % 产生二维坐标 
s = mesh(ZuoBiao_x, ZuoBiao_y, ZuoBiao_z);      % 画网格曲面图
s.FaceColor = 'flat'; % 修改网格图的属性
hold on;
scatter3(x(:,1), x(:,2), f_xy(x(:,1), x(:,2)), 'ro');
title('最终状态图')
 
disp(['最大值：',num2str(fg_best)]);
disp(['变量取值：',num2str(f_best)]);
 
% 下面一段用来输出 粒子群最佳适应度的排行榜
[socre,ind] = sort(fp_best, 'descend')  % 降序排序
DaAn = zeros(3,2);
for i=1:N
    row = ind(i);
    DaAn(i,:) = p_best(row,:);
end