基于粒子群算法的机组组合计算

% 问题模型 % （1）数学模型： % 以IEEE 3机6节点为工况模型 ，Load=850MW % （2）目标函数： % minF=∑(i=1，Ng)Fi(PGi)，Fi(PGi) =aiPi2+biPi+ci+Ei，考虑阀点(Valve-Point)效应 % （3）约束条件： % a.发电机组输出功率上下限约束，即PGi min<=PGi<= PGimax； % b.电力负荷平衡约束； % c.忽略网损，即∑(i=1，Ng)（PGi）= PGD； function pso %主程序 clear; clc; format long; xGbest=0; yPbest=[]; bound=[100 600;100 400;50 200;50 200]; %上下限 a=[0.00156;0.00194;0.00482;0.00482] ; b=[7.92;7.85;7.97;7.97]; c=[561;310;78;78]; g=[300;200;150;150]; h=[0.0315;0.042;0.063;0.063]; load=500; %负荷 popsize=20 ; %粒子数 MaxDD=100; %迭代次数 c1=1.4962; %学习因子 c2=1.4962; %学习因子 w=0.729; %惯性权重 dimsize=length(bound); vmax=vmaxinitial(dimsize,bound); %vmax产生过程 pop=initial(dimsize,popsize,bound,load); %初始种群 cost=obj(pop,dimsize,popsize,a,b,c,g,h,bound) ; %初始种群目标函数值 Pbest=pop(:,1:dimsize); %初始Pbest 和 Gbest [ggBest,xindex]=min(cost); xtemp=pop(xindex,1:dimsize); Gbest=xtemp; v=vinitial(popsize,dimsize,vmax); %初始速度种群 for m=1:MaxDD for t=1:popsize while(1) vv=vrenew(vmax,dimsize,Pbest,pop,Gbest,t,w,v,c1,c2); %进化速度 ppop=poprenew(vv,dimsize,t,pop,bound,load) ; %进化位置 if(ppop(t,1)<bound(1,2)&ppop(t,1)>bound(1,1)) break; end end ccost=hobjrenew(t,dimsize,a,ppop,b,c,g,h,bound) ; %进化目标函数值 pcost=hPestobj(t,dimsize,a,Pbest,b,c,g,h,bound) ; %Pbest的目标函数值 Gcost=hGbestobj(dimsize,a,Gbest,b,c,g,h,bound) ; %Gbest的目标函数值 if ccost<pcost %Pbest的更新 Pbest(t,1:dimsize)=ppop(t,1:dimsize); end if ccost<Gcost %Gbest的更新 Gbest=ppop(t,1:dimsize); GGcost=ccost; end pop(t,:)=ppop(t,:); v(t,dimsize:2)=vv(t,dimsize:2); end ggBest(m,:)=Gcost; xtemp(m,:)=Gbest; end yPbest=xtemp(MaxDD,:); xGbest=ggBest(MaxDD); yPbest,xGbest,pop plot(1:MaxDD,ggBest) function vmax=vmaxinitial(dimsize,bound) %子程序 for i=1:dimsize %vmax产生过程 vmax(i)=(bound(i,2)-bound(i,1))/20; end end function pop=initial(dimsize,popsize,bound,load) while (1) %初始种群 for i=dimsize:-1:2 pop(:,i)=rand(popsize,1)*(bound(i,2)-bound(i,1))+bound(i,1); end i=dimsize;s=zeros(popsize,1); while i>=2 s=s+pop(:,i); i=i-1; end pop(:,1)=load*ones(popsize,1)-s; if(pop(:,1)<ones(popsize,1)*bound(1,2)&pop(:,1)>ones(popsize,1)*bound(1,1)) break; end end end end function cost=obj(pop,dimsize,popsize,a,b,c,g,h,bound) t=1;cost=zeros(popsize,1); %初始种群目标函数值 while t<=dimsize cost=cost+a(t)*pop(:,t).^2+b(t)*pop(:,t)+c(t)+abs(g(t)*sin(h(t)*(pop(:,t)-bound(t,1)*ones(popsize,1)))); t=t+1; end end function v=vinitial(popsize,dimsize,vmax) %初始速度种群 vv=zeros(popsize,1);i=dimsize; while(1) while i>=2 v(:,i)=rand(popsize,1)*vmax(i); vv=vv+v(:,i); i=i-1; end v(:,1)=-vv; if(abs(v(:,1))<=ones(popsize,1)*vmax(1)) break; end end end function vv=vrenew(vmax,dimsize,Pbest,pop,Gbest,t,w,v,c1,c2) %进化速度 for dimindex=dimsize:-1:2 sub1=Pbest(t,dimindex)-pop(t,dimindex); sub2=Gbest(1,dimindex)-pop(t,dimindex); tempv(t,dimindex)=w*v(t,dimindex)+c1*rand*sub1+c2*rand*sub2; if tempv(t,dimindex)>vmax(dimindex) vv(t,dimindex)=vmax(dimindex); elseif tempv(t,dimindex)<-vmax(dimindex) vv(t,dimindex)=-vmax(dimindex); else vv(t,dimindex)=tempv(t,dimindex); end end end function ppop=poprenew(vv,dimsize,t,pop,bound,load) %进化位置 ppp=zeros; for dimindex=dimsize:-1:2 pp(t,dimindex)=pop(t,dimindex)+vv(t,dimindex); if pp(t,dimindex)>bound(dimindex,2) ppop(t,dimindex)=bound(dimindex,2); elseif pp(t,dimindex)<bound(dimindex,1) ppop(t,dimindex)=bound(dimindex,1); else ppop(t,dimindex)=pp(t,dimindex); end ppp=ppp+ppop(t,dimindex); end ppop(t,1)=load-ppp; end function ccost=hobjrenew(t,dimsize,a,ppop,b,c,g,h,bound) %进化目标函数值 ccost=zeros; for i=1:dimsize ccost=ccost+a(i)*ppop(t,i).^2+b(i)*ppop(t,i)+c(i)+abs(g(i)*sin(h(i)*(ppop(t,i)-bound(i,1)))); end end function pcost=hPestobj(t,dimsize,a,Pbest,b,c,g,h,bound) %Pbest的目标函数值 pcost=zeros; for i=1:dimsize pcost=pcost+a(i)*Pbest(t,i).^2+b(i)*Pbest(t,i)+c(i)+abs(g(i)*sin(h(i)*(Pbest(t,i)-bound(i,1)))); end end function Gcost=hGbestobj(dimsize,a,Gbest,b,c,g,h,bound) %Gbest的目标函数值 Gcost=zeros; for i=1:dimsize Gcost=Gcost+a(i)*Gbest(i).^2+b(i)*Gbest(i)+c(i)+abs(g(i)*sin(h(i)*(Gbest(i)-bound(i,1)))); end end