% 改进的快速粒子群优化算法 (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