粒子群算法

%%%粒子群伪代码
begin
    Initialize;                        % 包括初始化粒子群数、粒子初始速度和位置
    [x, xd] = judge(x, pop_size);      % 调用 judge 函数，初始化第一次值
    
    for num = 2: 最大迭代次数
        wk = wmax - num * (wmax - wmin) / max_gen;   % 计算惯性权重
        r1 = rand();                   % 随机产生加速权重
        r2 = rand();                   % 随机产生加速权重
        
        PSO 算法
        迭代求 vk, xk;
        
        While 判断 vk 是否满足条件
            再次重新生成加速权重系数 r1, r2;
            
            PSO 算法
            再次迭代求 vk, xk 数值
        end
        
        调用 [x, xd] = judge(x, pop_size); 重新计算目标函数数值
        
        判断并重新生成 pj 数值;
        判断并重新生成 pjd 数值;
        如果迭代前数值 > 迭代后的数值
            累加迭代次数
        end
    end
    
    输出随机数种子、进度、最优迭代次数，每个函数的数值和目标函数的数值
    用 ASCII 保存粒子位移的数值
    用 ASCII 保存粒子速度的数值
end

function [xm, fv] = PSO(fitness, N, c1, c2, w, M, D)
% 定义初始化条件
% c1 学习因子1
% c2 学习因子2
% w 惯性权重
% M 最大迭代次数
% D 搜索空间维数，自变量个数
%xm目标函数取最小值时的自变量
%fv目标函数的最小值
% N 初始种群个体数量
% 初始化粒子群的位置和速度（可以在这里限定位置和速度的范围）

format long;
for i = 1:N
    for j = 1:D
        x(i,j) = randn;  % 随机初始化位置
        v(i,j) = randn;  % 随机初始化速度
    end
end

% 计算各个粒子的适应度，并初始化 Pi 和 Pg
for i = 1:N
    p(i) = fitness(x(i,:));
    y(i,:) = x(i,:);
end

pg = x(N,:);  % pg 为全局最优
for i = 1:(N-1)
    if fitness(x(i,:)) < fitness(pg)
        pg = x(i,:);
    end
end

% 进入主要循环，根据公式依次迭代，直到满足精度要求
for t = 1:M
    for i = 1:N
        % 更新速度、位移
        v(i,:) = w*v(i,:) + c1*rand*(y(i,:)-x(i,:)) + c2*rand*(pg-x(i,:));
        x(i,:) = x(i,:) + v(i,:);
        if fitness(x(i,:)) < p(i)
            p(i) = fitness(x(i,:));
            y(i,:) = x(i,:);
        end
        if p(i) < fitness(pg)
            pg = y(i,:);
        end
    end
    Pbest(t) = fitness(pg);
end

% 最后给出计算结果
disp('*******************************************');
disp('目标函数取最小值时的自变量：');
xm = pg;
disp(xm);
disp('目标函数的最小值为：');
fv = fitness(pg);
disp(fv);
disp('*******************************************');
end

两个适应度函数

% 定义GW函数
function y = GW(x)
    [row, col] = size(x); % 获取 x 的大小
    if row > 1
        error('输入参数错误'); % 如果输入的 x 不是行向量，抛出错误
    end
    y1 = 1/4000 * sum(x.^2); % 计算 y1 = (1/4000) * x 各元素平方和
    y2 = 1; % 初始化 y2 为 1
    for h = 1:col
        y2 = y2 * cos(x(h)/sqrt(h)); % 逐元素更新 y2，乘以 cos(x(h)/sqrt(h))
    end
    y = y1 - y2 + 1; % 计算最终的 y 值
    y = - y; % 取 y 的相反数，返回
end

% 主脚本
clear all
clc

x = [-7:0.05:7];
y = x;
[X,Y] = meshgrid(x,y); %二维网格绘图
[row,col] = size(X);

for l = 1:col
    for h = 1:row
        z(h,l) = GW([X(h,l),Y(h,l)]);
    end
end

surf(X,Y,z); %三维曲面图
shading interp  %使曲面图的着色平滑过渡。

% 定义RA函数
function y = RA(x)
    [row, col] = size(x);
    if row > 1
        error('输入的参数错误');
    end
    y = sum(x.^2 - 10 * cos(2 * pi * x) + 10);
    y = - y;
end

% 主脚本
clear all
clc

x = [-5:0.1:5];
y = x;
[X,Y] = meshgrid(x,y);
[row,col] = size(X);

for l = 1:col
    for h = 1:row
        z(h,l) = RA([X(h,l),Y(h,l)]);
    end
end

surf(X,Y,z);
shading interp