粒子群算法的简单逻辑框架

大致的按照此计算步骤实现：

就对应上面的标准的PSO流程图和伪代码进行编程：

例一：求：f=x^2+y^2;   -10<x,y<10   求取min
例二：求：f = ((sin(x^2+y^2))^2-cos(x^2+y^2)+1)/((1+0.1*(x^2+y^2))^2)-0.7;   最小值

clc;                       %清屏
clear all;                 %清除所有变量
close all;                 %清图

N = 1000;                   %粒子群个数
D = 2;                     %粒子维数
T = 200;                   %最大迭代次数
c1 = 1.5;                  %学习因子 1.5
c2 = 1.5;                  %学习因子 1.5
w = 0.8;                   %惯性权重
Xmax = 10;                 %位置最大值
Xmin = -10;                %位置最小值
Vmax = 10;                 %速度最大值
Vmin = -10;                %速度最小值
%%%%%%%%%%%初始化种群个体(限定位置和速度)%%%%%%%%%%%%%
% x =rand(N,D) * (Xmax-Xmin)+Xmin;        %初始化位置
% v =rand(N,D) * (Vmax-Vmin)+Vmin;        %初始化速度
x =rand(N,D) * (Xmax-Xmin)+Xmin;        %初始化位置
v =rand(N,D) * (Vmax-Vmin)+Vmin;        %初始化速度
%%%%%%%%%%%随机初始化每个粒子以及评估每个粒子得到全局最优%%%%%%%%%%%%%%%
pbest_value = ones(N,1);
pbest=ones(N,2);
for i = 1:N
    pbest(i,:)=x(i,:);
    pbest_value(i) = fitness(x(i,1),x(i,2));
end
[gbest_value,gbest_index]=min(pbest_value);
gbest=x(gbest_index,:);
%%%%%%%%%%%%%%%%%%开始迭代%%%%%%%%%%%%%%%%%
for i = 1:T
 for j = 1:N
 %%%%%%%%%%%%更新位置和速度值%%%%%%%%%%%%%%%
        v(j,:) = w*v(j,:)+c1*rand*(pbest(j,:)-x(j,:))+c2*rand*(gbest-x(j,:));
        x(j,:) = x(j,:)+v(j,:);
 %%%%%%%%%%%%%%边界条件处理%%%%%%%%%%%%%%%%%
   for ii=1:1:D
         if (v(j,ii) > Vmax)
                v(j,ii) = Vmax;
         end
          if (v(j,ii) < Vmin)
                v(j,ii) = Vmin;
         end
         if (x(j,ii) > Xmax)
                x(j,ii) = Xmax ;
         end
         if (x(j,ii) < Xmin)
                x(j,ii) = Xmin;
         end
   end
 %%%%%%%%%更新个体最优位置和最优值%%%%%%%%%%%%%
        if (fitness(x(j,1),x(j,2)) < pbest_value(j))
            pbest(j,:) = x(j,:);
            pbest_value(j) = fitness(x(j,1),x(j,2));
        end
 %%%%%%%%%%更新全局最优位置和最优值%%%%%%%%%%%%
        if (pbest_value(j) < gbest_value)
            gbest_value = pbest_value(j);
            gbest = pbest(j,:);
        end
 end  
 gb(i)=gbest_value;
end
 plot(gb)

% function fitvalue=fitness(x,y)
% result=x^2+y^2;
% fitvalue=result;
% end
function y = fitness(x,y)
    y = ((sin(x^2+y^2))^2-cos(x^2+y^2)+1)/((1+0.1*(x^2+y^2))^2)-0.7;
end