%标准粒群优化算法程序
%测试函数:f(x,y)=100(x^2-y)^2+(1-x)^2, -2.048<x,y<2.048
%求解函数最小值
global popsize; %种群规模
%global popnum; %种群数量
global pop; %种群
%global c0; %速度惯性系数,为0—1的随机数
global c1; %个体最优导向系数
global c2; %全局最优导向系数
global gbest_x; %全局最优解x轴坐标
global gbest_y; %全局最优解y轴坐标
global best_fitness; %最优解
global best_in_history; %最优解变化轨迹
global x_min; %x的下限
global x_max; %x的上限
global y_min; %y的下限
global y_max; %y的上限
global gen; %迭代次数
global exetime; %当前迭代次数
global max_velocity; %最大速度
initial; %初始化
for exetime=1:gen
outputdata; %实时输出结果
adapting; %计算适应值
errorcompute(); %计算当前种群适值标准差
updatepop; %更新粒子位置
pause(0.01);
end
clear i;
clear exetime;
clear x_max;
clear x_min;
clear y_min;
clear y_max;
%程序初始化
gen=100; %设置进化代数
popsize=30; %设置种群规模大小
best_in_history(gen)=inf; %初始化全局历史最优解
best_in_history(=inf; %初始化全局历史最优解
max_velocity=0.3; %最大速度限制
best_fitness=inf;
%popnum=1; %设置种群数量
pop(popsize,8)=0; %初始化种群,创建popsize行5列的0矩阵
%种群数组第1列为x轴坐标,第2列为y轴坐标,第3列为x轴速度分量,第4列为y轴速度分量
%第5列为个体最优位置的x轴坐标,第6列为个体最优位置的y轴坐标
%第7列为个体最优适值,第8列为当前个体适应值
for i=1:popsize
pop(i,1)=4*rand()-2; %初始化种群中的粒子位置,值为-2—2,步长为其速度
pop(i,2)=4*rand()-2; %初始化种群中的粒子位置,值为-2—2,步长为其速度
pop(i,5)=pop(i,1); %初始状态下个体最优值等于初始位置
pop(i,6)=pop(i,2); %初始状态下个体最优值等于初始位置
pop(i,3)=rand()*0.02-0.01; %初始化种群微粒速度,值为-0.01—0.01,间隔为0.0001
pop(i,4)=rand()*0.02-0.01; %初始化种群微粒速度,值为-0.01—0.01,间隔为0.0001
pop(i,7)=inf;
pop(i,8)=inf;
end
c1=2;
c2=2;
x_min=-2;
y_min=-2;
x_max=2;
y_max=2;
gbest_x=pop(1,1); %全局最优初始值为种群第一个粒子的位置
gbest_y=pop(1,2);
%适值计算
% 测试函数为f(x,y)=100(x^2-y)^2+(1-x)^2, -2.048<x,y<2.048
%计算适应值并赋值
for i=1:popsize
pop(i,8)=100*(pop(i,1)^2-pop(i,2))^2+(1-pop(i,1))^2;
if pop(i,7)>pop(i,8) %若当前适应值优于个体最优值,则进行个体最优信息的更新
pop(i,7)=pop(i,8); %适值更新
pop(i,5:6)=pop(i,1:2); %位置坐标更新
end
end
%计算完适应值后寻找当前全局最优位置并记录其坐标
if best_fitness>min(pop(:,7))
best_fitness=min(pop(:,7)); %全局最优值
gbest_x=pop(find(pop(:,7)==min(pop(:,7))),1); %全局最优粒子的位置
gbest_y=pop(find(pop(:,7)==min(pop(:,7))),2);
end
best_in_history(exetime)=best_fitness; %记录当前全局最优
%实时输出结果
%输出当前种群中粒子位置
subplot(1,2,1);
for i=1:popsize
plot(pop(i,1),pop(i,2),'b*');
hold on;
end
plot(gbest_x,gbest_y,'r.','markersize',20);axis([-2,2,-2,2]);
hold off;
subplot(1,2,2);
axis([0,gen,-0.00005,0.00005]);
if exetime-1>0
line([exetime-1,exetime],[best_in_history(exetime-1),best_fitness]);hold on;
end
%粒子群速度与位置更新
%更新粒子速度
for i=1:popsize
pop(i,3)=rand()*pop(i,3)+c1*rand()*(pop(i,5)-pop(i,1))+c2*rand()*(gbest_x-pop(i,1)); %更新速度
pop(i,4)=rand()*pop(i,4)+c1*rand()*(pop(i,6)-pop(i,2))+c2*rand()*(gbest_x-pop(i,2));
if abs(pop(i,3))>max_velocity
if pop(i,3)>0
pop(i,3)=max_velocity;
else
pop(i,3)=-max_velocity;
end
end
if abs(pop(i,4))>max_velocity
if pop(i,4)>0
pop(i,4)=max_velocity;
else
pop(i,4)=-max_velocity;
end
end
end
%更新粒子位置
for i=1:popsize
pop(i,1)=pop(i,1)+pop(i,3);
pop(i,2)=pop(i,2)+pop(i,4);
end
% A SIMPLE IMPLEMENTATION OF THE
% Particle Swarm Optimization IN MATLAB
function [xmin, fxmin, iter] = PSO()
% Initializing variables
success = 0; % Success flag
PopSize = 30; % Size of the swarm
MaxIt = 100; % Maximum number of iterations
iter = 0; % Iterations’counter
fevals = 0; % Function evaluations’ counter
c1 = 2; % PSO parameter C1
c2 = 2; % PSO parameter C2
w = 0.6; % inertia weight
% Objective Function
f = ^DeJong^;
dim = 10; % Dimension of the problem
upbnd = 10; % Upper bound for init. of the swarm
lwbnd = -5; % Lower bound for init. of the swarm
GM = 0; % Global minimum (used in the stopping criterion)
ErrGoal = 0.0001; % Desired accuracy
% Initializing swarm and velocities
popul = rand(dim, PopSize)*(upbnd-lwbnd) + lwbnd;
vel = rand(dim, PopSize);
for i = 1:PopSize,
fpopul(i) = feval(f, popul(:,i));
fevals = fevals + 1;
end
bestpos = popul;
fbestpos = fpopul;
% Finding best particle in initial population
[fbestpart,g] = min(fpopul);
lastbpf = fbestpart;
while (success == 0) & (iter < MaxIt),
iter = iter + 1;
% VELOCITY UPDATE
for i=1:PopSize,
A(:,i) = bestpos(:,g);
end
R1 = rand(dim, PopSize);
R2 = rand(dim, PopSize);
vel = w*vel + c1*R1.*(bestpos-popul) + c2*R2.*(A-popul);
% SWARMUPDATE
popul = popul + vel;
% Evaluate the new swarm
for i = 1:PopSize,
fpopul(i) = feval(f,popul(:, i));
fevals = fevals + 1;
end
% Updating the best position for each particle
changeColumns = fpopul < fbestpos;
fbestpos = fbestpos.*( 1-changeColumns) + fpopul.*changeColumns;
bestpos(:, find(changeColumns)) = popul(:, find(changeColumns));
% Updating index g
[fbestpart, g] = min(fbestpos);
currentTime = etime(clock,startTime);
% Checking stopping criterion
if abs(fbestpart-GM) <= ErrGoal
success = 1;
else
lastbpf = fbestpart;
end
end
% Output arguments
xmin = popul(:,g);
fxmin = fbestpos(g);
fprintf(^ The best vector is : \n^);
fprintf(^--- %g ^,xmin);
fprintf(^\n^);
%==========================================
function DeJong=DeJong(x)
DeJong = sum(x.^2);