% 改進的快速粒子群優化算法 (APSO):

function apso

Lb=[6 50??0.032??0.040 0.007 0.08 300 ]; %下邊界

Ub=[8 64??0.038??0.050 0.010 0.13 350]; %上邊界

% 默認參數

para=[40 100 0.95]; %[粒子數，迭代次數，gama參數]

% APSO 優化求解函數

[gbest,fmin]=pso_mincon(@result,@constraint,Lb,Ub,para);

% 輸出結果

Bestsolution=gbest % 全局最優個體

fmin

%% 目標函數

function f=result(x)

f=-((1.5*10^(-3)*sin(x(1)*pi()/180))/(((x(7)/((7.5398*10^8)*x(3)*x(5)*sin(pi()*x(1)/180)))*(99.591*x(3))*(1-exp(-3.9*10^(-3)/x(3)))*(sqrt(1+9*(x(6))^2))*(3*10^(-7)+1.2582*10^(-6)*(sin(pi()*x(1)/180)/x(7))^0.3 ))*(54.3914/(25.2340*x(7)*x(6)*x(3)/(sin(x(1)*pi()/180))+32.4657))))

% 非線性約束

function [g,geq]=constraint(x)

% 不等式限制條件

A=tan(x(2)*pi()/180);

B=x(4)*sin(x(1)*pi()/180)-0.1*x(6)*x(3);

C=0.1*x(4)*sin(x(1)*pi()/180)+x(6)*x(3);

D=sin(x(1)*pi()/180)*2*pi()*x(5)*x(3);

R=x(7)/D;

Q=tan(x(1)*pi()/180);

F=B/C

S=x(3)/x(5)

t=54.407/((x(7)*x(6)*x(3)*25.234/sin(x(1)*pi()/180))+32.4657)

g(1)=F-A;

g(2)=R-2.5;

g(3)=x(6)-Q;

g(4)=2.5-S;

g(5)=S-4;

g(6)=0.15-t;

g(7)=t-0.3;

g(8)=50-x(7);

g(9)=x(7)-350;

% 如果沒有等式約束，則置geq=[]；

geq=[];

%%??APSO Solver

function [gbest,fbest]=pso_mincon(fhandle,fnonlin,Lb,Ub,para)

if nargin<=7,? ?? ?%用來判定輸入變量的個數

para=[40 200 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