function [P_gg,fV,fPg,fQg,fPloss,fSwarm]=IAPSO()
%-----------基于免疫粒子群算法的无功优化程序------------------------
%-----------包含免疫记忆和调节,以及疫苗接种技术------------
%------------改进:分为多个子程序--------------------------------
t1 = clock;
%------给定算法参数----------------------------------------------
N=30;
M=10;
R=20;
w_min=0.4;
w_max=0.9;
iter_max=100;
iterdec=(w_max-w_min)./iter_max;
iter = 0;
%--------粒子群维数---------------
[gen_U, T, C_Q] = control;
Ng=size(gen_U,1);
Nt=size(T,1);
Nc=size(C_Q,1);
D=Ng+Nt+Nc;
%---------初始化粒子形成初始群体P0----------------
[x,Swarm,v,v_max]=inition(N,D);
%-------------计算当前粒子群体P0中粒子的适应度值--------
[fV,fPg,fQg,fPloss,fSwarm]=fitness(N,D,Swarm);
%------------生成免疫记忆粒子(抗体)------------
%-------------得到个体最优Pi----------------
PBest = x;
outPBest = Swarm;
fPBest = fSwarm;
%-------------得到全局最优Pg----------------
[fGBest, g] = min(fPBest);
GBest = PBest(g,:);
outGBest = outPBest(g,:);
%-------------------------开始迭代-----------------------------
while iter <iter_max
iter = iter+1;
time(iter)=iter;
%------------将全局最优位置作为疫苗V------------
bacterinx=GBest;
bacterinSwarm=outGBest;
bacterinv=v(g,:);
%-------------------修改惯性系数--------------------------------------
w_now = w_max - iterdec*iter;
%-------------------新粒子的生成----------------
%--------第一种:由PSO的更新公式产生N个新粒子----------------
[x1,v1,Swarm1]=renew(w_now,N,D,x,v,v_max,PBest,GBest);
%---------计算更新后粒子的适应值,并更新个体最优-------------
[fV,fPg,fQg,fPloss,fSwarm1]=fitness(N,D,Swarm1);
changeRows = fSwarm1 < fPBest;
fPBest(find(changeRows)) = fSwarm1(find(changeRows));
%------个体最优位置PBest的更新----------------------
PBest(find(changeRows), :) = x1(find(changeRows), :);
outPBest(find(changeRows), :) = Swarm1(find(changeRows), :);
%-------------得到个体最优Pi----------------
PBest1 = PBest;
outPBest1 = outPBest;
fPBest1 = fPBest;
%--------第二种:随机产生M个新粒子----------------
[x2,Swarm2,v2,v_max]=inition(M,D);
%---------计算随机产生的M个粒子的适应值,并更新个体最优-------------
[fV,fPg,fQg,fPloss,fSwarm2]=fitness(M,D,Swarm2);
%-------------得到个体最优Pi----------------
PBest2 = x2;
outPBest2 = Swarm2;
fPBest2 = fSwarm2;
%--------生成N+M个粒子的群体----------------
x=[x1;x2];
Swarm=[Swarm1;Swarm2];
v=[v1;v2];
PBest=[PBest1;PBest2];
outPBest=[outPBest1;outPBest2];
fPBest=[fPBest1,fPBest2];
fSwarm=[fSwarm1,fSwarm2];
%--------计算生成的N+M个粒子的选择概率--------
%--------基于粒子浓度的选择概率-----------
P=probabilitychroma(fSwarm,N,M);
%--------依概率大小选择N个粒子形成新的群体Ak---------
[x,v,Swarm,fSwarm,fPBest,PBest,outPBest]=choose(N,M,P,x,v,Swarm,fSwarm,fPBest,PBest,outPBest);
%----------保留群体Ak------------
Qx=x;
QSwarm=Swarm;
Qv=v;
QfPBest=fPBest;
QPBest=PBest;
QoutPBest=outPBest;
%-----------------接种疫苗,形成新的群体Rk-----------------------
k=1;
while k<=R
randSwarm=round(unifrnd(1,N));
randplace1=round(unifrnd(1,D));
randplace2=round(unifrnd(1,D));
randplace3=round(unifrnd(1,D));
x(randSwarm,randplace1)=bacterinx(randplace1);
x(randSwarm,randplace2)=bacterinx(randplace2);
x(randSwarm,randplace3)=bacterinx(randplace3);
Swarm(randSwarm,randplace1)=bacterinSwarm(randplace1);
Swarm(randSwarm,randplace2)=bacterinSwarm(randplace2);
Swarm(randSwarm,randplace3)=bacterinSwarm(randplace3);
v(randSwarm,randplace1)=bacterinv(randplace1);
v(randSwarm,randplace2)=bacterinv(randplace2);
v(randSwarm,randplace3)=bacterinv(randplace3);
k=k+1;
end
%---------计算接种后形成的群体Rk的适应值-------------
[fV,fPg,fQg,fPloss,fSwarm]=fitness(N,D,Swarm);
%------------------免疫选择---------------------
%---------适应值比接种前的差则放弃,否则保留-----------
changeRows = fSwarm < QfPBest;
Qx(find(changeRows), :) = x(find(changeRows), :);
Qv(find(changeRows), :) = v(find(changeRows), :);
QSwarm(find(changeRows), :) = Swarm(find(changeRows), :);
QfPBest(find(changeRows)) = fSwarm(find(changeRows));
QPBest(find(changeRows), :) = x(find(changeRows), :);
QoutPBest(find(changeRows), :) = Swarm(find(changeRows), :);
%-------------形成新的群体Pk----------------
x=Qx;
Swarm=QSwarm;
v=Qv;
fPBest=QfPBest;
PBest=QPBest;
outPBest= QoutPBest;
%----------得到全局最优值--------------
[fGBest, g] = min(fPBest);
GBest = PBest(g,:);
outGBest = outPBest(g,:);
%----------取得最优时的网损---------
[fV,fPg,fQg,fPloss,fSwarm]=fitness(1,D,outGBest);
P_gg(iter)=fPloss;
%------输出结果----------------------
disp('控制变量取值Best=:');
disp(outGBest);
disp('目标函数最优值fGBest=:');
disp(fGBest);
end
%-----------------------结束循环-----------------------------
et = etime(clock, t1);
%----------------最后取得的全局最优解,对应的控制变量、网损、发电机无功出力、节点电压的取值--------------
[fV,fPg,fQg,fPloss,fSwarm]=fitness(1,D,outGBest);
disp('控制变量取值Best=:');
disp(outGBest);
disp('目标函数最优值fGBest=:');
disp(fGBest);
disp('网损Ploss=:');
disp(fPloss);
%---------------画‘评价值’变化曲线--------------------------