% 改進的快速粒子群優化算法 (APSO):
function apso
Lb=[0.1 0.1 0.1 0.1]; %下邊界
Ub=[2.0 10.0 10.0 2.0]; %上邊界
% 預設參數
para=[25 150 0.95]; %[粒子數，疊代次數，gama參數]
% APSO 優化求解函數
[gbest,fmin]=pso_mincon(@cost,@constraint,Lb,Ub,para);
% 輸出結果
Bestsolution=gbest % 全局最優個體
fmin
%% 目标函數
function f=cost(x)
f=1.10471*x(1)^2*x(2)+0.04811*x(3)*x(4)*(14.0+x(2));
% 非線性限制
function [g,geq]=constraint(x)
% 不等式限制條件
Q=6000*(14+x(2)/2);
D=sqrt(x(2)^2/4+(x(1)+x(3))^2/4);
J=2*(x(1)*x(2)*sqrt(2)*(x(2)^2/12+(x(1)+x(3))^2/4));
alpha=6000/(sqrt(2)*x(1)*x(2));
beta=Q*D/J;
tau=sqrt(alpha^2+2*alpha*beta*x(2)/(2*D)+beta^2);
sigma=504000/(x(4)*x(3)^2);
delta=65856000/(30*10^6*x(4)*x(3)^3);
F=4.013*(30*10^6)/196*sqrt(x(3)^2*x(4)^6/36)*(1-x(3)*sqrt(30/48)/28);
g(1)=tau-13600;
g(2)=sigma-30000;
g(3)=x(1)-x(4);
g(4)=0.10471*x(1)^2+0.04811*x(3)*x(4)*(14+x(2))-5.0;
g(5)=0.125-x(1);
g(6)=delta-0.25;
g(7)=6000-F;
% 如果沒有等式限制，則置geq=[]；
geq=[];
%% APSO Solver
function [gbest,fbest]=pso_mincon(fhandle,fnonlin,Lb,Ub,para)
if nargin<=4,
para=[20 150 0.95];
end
n=para(1);% 粒子種群大小
time=para(2); %時間步長，疊代次數
gamma=para(3); %gama參數
scale=abs(Ub-Lb); %取值區間
% 驗證限制條件是否合乎條件
if abs(length(Lb)-length(Ub))>0,
disp('Constraints must have equal size');
return
end
alpha=0.2; % alpha=[0,1]粒子随機衰減因子
beta=0.5; % 收斂速度(0->1)=(slow->fast);
% 初始化粒子群
best=init_pso(n,Lb,Ub);
fbest=1.0e+100;
% 疊代開始
for t=1:time,
%尋找全局最優個體
for i=1:n,
fval=Fun(fhandle,fnonlin,best(i,:));
% 更新最有個體
if fval<=fbest,
gbest=best(i,:);
fbest=fval;
end
end
% 随機性衰減因子
alpha=newPara(alpha,gamma);
% 更新粒子位置
best=pso_move(best,gbest,alpha,beta,Lb,Ub);
% 結果顯示
str=strcat('Best estimates: gbest=',num2str(gbest));
str=strcat(str,' iteration='); str=strcat(str,num2str(t));
disp(str);
fitness1(t)=fbest;
plot(fitness1,'r','Linewidth',2)
grid on
hold on
title('适應度')
end
% 初始化粒子函數
function [guess]=init_pso(n,Lb,Ub)
ndim=length(Lb);
for i=1:n,
guess(i,1:ndim)=Lb+rand(1,ndim).*(Ub-Lb);
end
%更新所有的粒子 toward (xo,yo)
function ns=pso_move(best,gbest,alpha,beta,Lb,Ub)
% 增加粒子在上下邊界區間内的随機性
n=size(best,1); ndim=size(best,2);
scale=(Ub-Lb);
for i=1:n,
ns(i,:)=best(i,:)+beta*(gbest-best(i,:))+alpha.*randn(1,ndim).*scale;
end
ns=findrange(ns,Lb,Ub);
% 邊界函數
function ns=findrange(ns,Lb,Ub)
n=length(ns);
for i=1:n,
% 下邊界限制
ns_tmp=ns(i,:);
I=ns_tmp
ns_tmp(I)=Lb(I);
% 上邊界限制
J=ns_tmp>Ub;
ns_tmp(J)=Ub(J);
%更新粒子
ns(i,:)=ns_tmp;
end
% 随機性衰減因子
function alpha=newPara(alpha,gamma);
alpha=alpha*gamma;
% 帶限制的d維目标函數的求解
function z=Fun(fhandle,fnonlin,u)
% 目标
z=fhandle(u);
z=z+getconstraints(fnonlin,u); % 非線性限制
function Z=getconstraints(fnonlin,u)
% 罰常數 >> 1
PEN=10^15;
lam=PEN; lameq=PEN;
Z=0;
% 非線性限制
[g,geq]=fnonlin(u);
%通過不等式限制建立罰函數
for k=1:length(g),
Z=Z+ lam*g(k)^2*getH(g(k));
end
% 等式條件限制
for k=1:length(geq),
Z=Z+lameq*geq(k)^2*geteqH(geq(k));
end
% Test if inequalities
function H=getH(g)
if g<=0,
H=0;
else
H=1;
end
% Test if equalities hold
function H=geteqH(g)
if g==0,
H=0;
else
H=1;
end