粒子学习算法一维实现

%基本粒子群优化算法（Particle Swarm Optimization）-----------
%名称：基本粒子群优化算法（PSO）
%作用：求解优化问题
%作者：周楠
%时间：2009-12-17-3:00
clear all;
clc;
%横坐标长度百分比
x_rate=0.2;
%纵坐标长度百分比
y_rate=1;
%字体显示位置横坐标长度百分比
x1_show=0.5;
x2_show=0.5;
%字体显示位置纵坐标长度百分比
y1_show=0.2;
y2_show=0.9;
%x1的范围[a1 b1];
a1=-1;
b1=2;
r1=rand;
r2=rand;
c1=2; %学习因子1
c2=2; %学习因子2
%惯性权重
%w=0.05;
w_0=0.9;
w_end=0.4;
MaxDT=200; %最大迭代次数
N=50; %初始化群体个体数目
eps=10^(-16); %设置精度(在已知最小值时候用)，即运算精度
rate=0.2; %如果步长超出边界，则缩短为原来的rate
%初始化种群的个体(位置和速度)
%rand是0-1的均匀分布
%初始位置变量生成
x1=((b1-a1)*rand(N))+a1;
%速度
v=zeros(N);
%某次迭代的最优解初始化
pbest=zeros(MaxDT);
pbest(1)=fx(x1(1));
%某次迭代的最优解的位置初始化
x_pbest=zeros(MaxDT);
x_pbest(1)=x1(1);
%gbest为全局最优解的初始化
gbest=0;
%draw_gbest为画出全局最优解的初始化
draw_gbest=0;
%gbest为全局最优解的位置初始化
x_gbest=x1(1);
%进入主要循环，按照公式依次迭代，直到满足精度要求
for t=1:MaxDT
for i=1:N
w(i)=((w_0-w_end)*(MaxDT-i))/(MaxDT)+w_end;
v(i)=w(i)*v(i)+c1*r1*(x_pbest(t)-x1(i))+c2*r2*(x_gbest-x1(i));
%当设定固定w=0.05;
%v(i)=w*v(i)+c1*r1*(x_pbest(t)-x1(i))+c2*r2*(x_gbest-x1(i));
x1(i)=x1(i)+v(i);
%对于步长太长，导致变量越界后的处理
while (((x1(i))>b1)||((x1(i))<a1))
% 方法一：如果步长超出边界，则缩短为原来的rate
% x1(i)=x1(i)-v(i);
% v(i)=v(i)*rate;
% x1(i)=x1(i)+v(i)
% 方法二：把如果步长超出边界，则取边界附近某 个值
if ((x1(i))>b1)
x1(i)=x1(i)-rate;
end
if ((x1(i))<a1)
x1(i)=x1(i)+rate;
end
end
%保存第t次迭代的局部最优解
if (fx(x1(i))>pbest(t))
pbest(t)=fx(x1(i));
x_pbest(t)=x1(i);
end
end
%保存全局最优解
if (pbest(t)>gbest)
gbest=pbest(t);
x_gbest=x_pbest(t);
end
end
%画出最优值
draw_gbest=zeros(MaxDT,1);
draw_gbest(1)=pbest(1);
for i=2:MaxDT
if (draw_gbest<pbest(i))
draw_gbest(i)=pbest(i);
else
draw_gbest(i)=draw_gbest(i-1);
end
end
%最后给出计算结果
disp('函数的全局最优位置为：')
Solution=x_gbest
disp('最后得到的全局优化极值为：')
Result=gbest
figure
x=(1:MaxDT)';
plot(x,pbest,'b:');
title('粒子群优化算法')
xlabel('迭代次数')
ylabel('优化极值')
aa=min(min(pbest));
bb=max(gbest);
string_cell{1}='/fontsize{10}/fontname{arial}红实线：到目前为止遇到的最好的解';
text((MaxDT+MaxDT*x_rate)*x1_show,((aa+aa*y_rate)+((bb+bb*y_rate)-(aa+aa*y_rate))*y1_show),string_cell,'color','m')
string_cell{1}='/fontsize{10}/fontname{arial}绿虚线：当前代中最好的个体';
text((MaxDT+MaxDT*x_rate)*x1_show,((aa+aa*y_rate)+((bb+bb*y_rate)-(aa+aa*y_rate))*y1_show*y1_show),string_cell,'color','m')
cell_string{4}=['/fontsize{15}/rm全局最优位置为= ' num2str(x_gbest,6)];
cell_string{5}='/fontsize{15}';
cell_string{6}=['/fontsize{15}/rm全局优化极值= ' num2str(gbest,6)];
cell_string{7}='/fontsize{15}';
cell_string{8}=['/fontsize{15}/rm迭代次数= ' num2str(MaxDT,6)];
text((MaxDT+MaxDT*x_rate)*x2_show,((aa+aa*y_rate)+((bb+bb*y_rate)-(aa+aa*y_rate))*y2_show),cell_string,'color','b','HorizontalAlignment','center')
hold on
plot(x,draw_gbest,'r-')
axis([0,fix(MaxDT+MaxDT*x_rate),fix(aa+aa*y_rate),fix(bb+bb*y_rate)])
hold off
grid on

function [ output_args ] =fx1Tox2(x)
%UNTITLED2 Summary of this function goes here
% Detailed explanation goes here
%for i=1:D
output_args=x(i)*sin(10*pi*x(i))+1.0;
%end
end