[matlab] 3.粒子群优化算法

%% I. 清空环境
clc
clear all

%% II. 绘制目标函数曲线图
x = -5:0.01:5;
y = 2.1*(1-x+5*x.^3).*exp(-x.^2/2)-x.^2;
figure
plot(x, y)
hold on

%% III. 参数初始化 这里默认惯性因子为1
c1 = 1.49445;   %加速常数 取Eberhart的参考值
c2 = 1.49445;   %加速常数 取Eberhart的参考值

maxgen = 30;   % 进化次数  
sizepop = 100;   %种群规模

Vmax = 0.5;  %速度步长上限
Vmin = -0.5; %速度步长下限
popmax = 5; % x的上限
popmin = -5; % x的下限

%% IV. 产生初始粒子和速度
for i = 1:sizepop
    % 随机产生一个种群
    pop(i,:) = -10*rands(1)+5;    %产生[-5,5]的初始种群
    V(i,:) = 0.5 * rands(1);  %初始化[-0.5,0.5]速度
    % 利用fun.m子函数，计算适应度并存储到fitness(i)中
    fitness(i) = fun(pop(i,:));   
end

%% V. 个体极值和群体极值
[bestfitness, bestindex] = max(fitness);
zbest = pop(bestindex,:);   %全局最佳
gbest = pop;    %个体最佳
fitnessgbest = fitness;   %个体最佳适应度值
fitnesszbest = bestfitness;   %全局最佳适应度值

%% VI. 迭代寻优
for i = 1:maxgen
    
    for j = 1:sizepop
        % 速度更新
        V(j,:) = V(j,:) + c1*rand*(gbest(j,:) - pop(j,:)) + c2*rand*(zbest - pop(j,:));
        V(j,V(j,:)>Vmax) = Vmax;
        V(j,V(j,:) Vmin;
        
        % 种群更新
        pop(j,:) = pop(j,:) + V(j,:);
        pop(j,pop(j,:)>popmax) = popmax;
        pop(j,pop(j,:) popmin;
        
        % 适应度值更新
        fitness(j) = fun(pop(j,:)); 
    end
    
    for j = 1:sizepop    
        % 个体最优更新
        if fitness(j) > fitnessgbest(j)
            gbest(j,:) = pop(j,:);
            fitnessgbest(j) = fitness(j);
        end
        
        % 群体最优更新
        if fitness(j) > fitnesszbest
            zbest = pop(j,:);
            fitnesszbest = fitness(j);
        end
    end 
    yy(i) = fitnesszbest;          
end

%% VII. 输出结果并绘图
[fitnesszbest zbest];
plot(zbest, fitnesszbest,‘r.‘,‘MarkerSize‘,10)
plot(zbest, fitnesszbest,‘ro‘,‘MarkerSize‘,16)
x_text=[‘x=‘,num2str(zbest)]; % x横坐标转换为字符串
y_text=[‘y=‘,num2str(fitnesszbest)];  % y横坐标转换为字符串
max_text=char(‘全局最优‘,x_text,y_text);  % 生成标志最大值点的字符串
text(zbest+0.3, fitnesszbest-1.4,max_text) %图上绘出全局最优的数据值
figure
plot(yy)
title(‘最优个体适应度‘,‘fontsize‘,12);
xlabel(‘进化代数‘,‘fontsize‘,12);ylabel(‘适应度‘,‘fontsize‘,12);

function y = fun(x)
% 函数用于计算粒子适应度值
%x           input           输入粒子 
%y           output          粒子适应度值 
y = 2.1*(1-x+5*x.^3).*exp(-x.^2/2)-x.^2;