%標(biāo)準(zhǔn)粒群優(yōu)化算法程序
% 2007.1.9 By jxy
%測試函數(shù)：f(x,y)=100(x^2-y)^2+(1-x)^2, -2.048<x,y<2.048
%求解函數(shù)最小值

global popsize; %種群規(guī)模
%global popnum; %種群數(shù)量
global pop; %種群
%global c0; %速度慣性系數(shù),為0—1的隨機(jī)數(shù)
global c1; %個體最優(yōu)導(dǎo)向系數(shù)
global c2; %全局最優(yōu)導(dǎo)向系數(shù)
global gbest_x; %全局最優(yōu)解x軸坐標(biāo)
global gbest_y; %全局最優(yōu)解y軸坐標(biāo)
global best_fitness; %最優(yōu)解
global best_in_history; %最優(yōu)解變化軌跡
global x_min; %x的下限
global x_max; %x的上限
global y_min; %y的下限
global y_max; %y的上限
global gen; %迭代次數(shù)
global exetime; %當(dāng)前迭代次數(shù)
global max_velocity; %最大速度

initial; %初始化

for exetime=1:gen
outputdata; %實(shí)時輸出結(jié)果
adapting; %計(jì)算適應(yīng)值
errorcompute(); %計(jì)算當(dāng)前種群適值標(biāo)準(zhǔn)差
updatepop; %更新粒子位置
pause(0.01);
end

clear i;
clear exetime;
clear x_max;
clear x_min;
clear y_min;
clear y_max;

%程序初始化

gen=100; %設(shè)置進(jìn)化代數(shù)
popsize=30; %設(shè)置種群規(guī)模大小
best_in_history(gen)=inf; %初始化全局歷史最優(yōu)解
best_in_history(=inf; %初始化全局歷史最優(yōu)解
max_velocity=0.3; %最大速度限制
best_fitness=inf;
%popnum=1; %設(shè)置種群數(shù)量

pop(popsize,8)=0; %初始化種群,創(chuàng)建popsize行5列的0矩陣
%種群數(shù)組第1列為x軸坐標(biāo)，第2列為y軸坐標(biāo)，第3列為x軸速度分量，第4列為y軸速度分量
%第5列為個體最優(yōu)位置的x軸坐標(biāo)，第6列為個體最優(yōu)位置的y軸坐標(biāo)
%第7列為個體最優(yōu)適值，第8列為當(dāng)前個體適應(yīng)值

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); %初始狀態(tài)下個體最優(yōu)值等于初始位置
pop(i,6)=pop(i,2); %初始狀態(tài)下個體最優(yōu)值等于初始位置
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); %全局最優(yōu)初始值為種群第一個粒子的位置
gbest_y=pop(1,2);

%適值計(jì)算
% 測試函數(shù)為f(x,y)=100(x^2-y)^2+(1-x)^2, -2.048<x,y<2.048

%計(jì)算適應(yīng)值并賦值
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) %若當(dāng)前適應(yīng)值優(yōu)于個體最優(yōu)值，則進(jìn)行個體最優(yōu)信息的更新
pop(i,7)=pop(i,8); %適值更新
pop(i,5:6)=pop(i,1:2); %位置坐標(biāo)更新
end
end

%計(jì)算完適應(yīng)值后尋找當(dāng)前全局最優(yōu)位置并記錄其坐標(biāo)
if best_fitness>min(pop(:,7))
best_fitness=min(pop(:,7)); %全局最優(yōu)值
gbest_x=pop(find(pop(:,7)==min(pop(:,7))),1); %全局最優(yōu)粒子的位置 
gbest_y=pop(find(pop(:,7)==min(pop(:,7))),2);
end

best_in_history(exetime)=best_fitness; %記錄當(dāng)前全局最優(yōu)

%實(shí)時輸出結(jié)果

%輸出當(dāng)前種群中粒子位置
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%標(biāo)準(zhǔn)粒群優(yōu)化算法程序
% 2007.1.9 By jxy
%測試函數(shù)：f(x,y)=100(x^2-y)^2+(1-x)^2, -2.048<x,y<2.048
%求解函數(shù)最小值

global popsize; %種群規(guī)模
%global popnum; %種群數(shù)量
global pop; %種群
%global c0; %速度慣性系數(shù),為0—1的隨機(jī)數(shù)
global c1; %個體最優(yōu)導(dǎo)向系數(shù)
global c2; %全局最優(yōu)導(dǎo)向系數(shù)
global gbest_x; %全局最優(yōu)解x軸坐標(biāo)
global gbest_y; %全局最優(yōu)解y軸坐標(biāo)
global best_fitness; %最優(yōu)解
global best_in_history; %最優(yōu)解變化軌跡
global x_min; %x的下限
global x_max; %x的上限
global y_min; %y的下限
global y_max; %y的上限
global gen; %迭代次數(shù)
global exetime; %當(dāng)前迭代次數(shù)
global max_velocity; %最大速度

initial; %初始化

for exetime=1:gen
outputdata; %實(shí)時輸出結(jié)果
adapting; %計(jì)算適應(yīng)值
errorcompute(); %計(jì)算當(dāng)前種群適值標(biāo)準(zhǔn)差
updatepop; %更新粒子位置
pause(0.01);
end

clear i;
clear exetime;
clear x_max;
clear x_min;
clear y_min;
clear y_max;

%程序初始化

gen=100; %設(shè)置進(jìn)化代數(shù)
popsize=30; %設(shè)置種群規(guī)模大小
best_in_history(gen)=inf; %初始化全局歷史最優(yōu)解
best_in_history(=inf; %初始化全局歷史最優(yōu)解
max_velocity=0.3; %最大速度限制
best_fitness=inf;
%popnum=1; %設(shè)置種群數(shù)量

pop(popsize,8)=0; %初始化種群,創(chuàng)建popsize行5列的0矩陣
%種群數(shù)組第1列為x軸坐標(biāo)，第2列為y軸坐標(biāo)，第3列為x軸速度分量，第4列為y軸速度分量
%第5列為個體最優(yōu)位置的x軸坐標(biāo)，第6列為個體最優(yōu)位置的y軸坐標(biāo)
%第7列為個體最優(yōu)適值，第8列為當(dāng)前個體適應(yīng)值

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); %初始狀態(tài)下個體最優(yōu)值等于初始位置
pop(i,6)=pop(i,2); %初始狀態(tài)下個體最優(yōu)值等于初始位置
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); %全局最優(yōu)初始值為種群第一個粒子的位置
gbest_y=pop(1,2);

%適值計(jì)算
% 測試函數(shù)為f(x,y)=100(x^2-y)^2+(1-x)^2, -2.048<x,y<2.048

%計(jì)算適應(yīng)值并賦值
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) %若當(dāng)前適應(yīng)值優(yōu)于個體最優(yōu)值，則進(jìn)行個體最優(yōu)信息的更新
pop(i,7)=pop(i,8); %適值更新
pop(i,5:6)=pop(i,1:2); %位置坐標(biāo)更新
end
end

%計(jì)算完適應(yīng)值后尋找當(dāng)前全局最優(yōu)位置并記錄其坐標(biāo)
if best_fitness>min(pop(:,7))
best_fitness=min(pop(:,7)); %全局最優(yōu)值
gbest_x=pop(find(pop(:,7)==min(pop(:,7))),1); %全局最優(yōu)粒子的位置 
gbest_y=pop(find(pop(:,7)==min(pop(:,7))),2);
end

best_in_history(exetime)=best_fitness; %記錄當(dāng)前全局最優(yōu)

%實(shí)時輸出結(jié)果

%輸出當(dāng)前種群中粒子位置
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