微电网多目标调度，运用改进的多目标粒子群算法进行计算，里面包含经济等三个目标函数设置

function microgird_dispatch_2() %第一行 DG，第二行MT，第三行FC tic; w_start=0.9; w_end=0.4; c1=2; c2=2; Vmax=200; PopSize=20; MaxIter=300; iter=0; T=24; N=4; %可再生能源的种类数 % 输入原始数据,包括各时段负荷大小，发电机有功输出上下限，发电机耗量成本系数，各时段风电场预测的平均输出功率 pmax(1)=7;pmin(1)=0; pmax(2)=4;pmin(2)=0; pmax(3)=4;pmin(3)=0; % a=[1000,970,700,680,450,370,480,660,665,670]; % b=[16.19,17.26,16.60,16.50,19.70,22.26,27.74,25.92,27.27,27.79]; % c=[0.00048,0.00031,0.002,0.00211,0.00398,0.00712,0.00079,0.00413,0.00222,0.00173]; % 柴油发电机的耗量参数 a=0.0071;b=0.2333;c=0.4333; % 微型燃气轮机,P_mt表示燃气轮机发出的功率，Xl_mt表示燃气轮机的效率 % Xl_mt=0.0753*(P_mt/65)^3-0.3095*(P_mt/65)^2+0.1068; % 燃料电池 Xl_fc=0.4; % Cost_fc=Price_fc*P_fc/Xl_fc; xx=1:T; %负荷数据 % PL=[11,4,4,4,5,6,7.5,8,7.5,6,5,6,6.5,7,8,5,8.5,11,13,13,10,10,6.5,5]; PL=[3,4,4,4,5,6,6.5,7,7.5,8.5,9,10,10.5,10,9,8.5,9,10,11,11.5,10,9,5.5,5]; %PL=PL*10; plot(xx,PL,'-*r'); xlabel('t/h');ylabel('P/kW'); axis([1,24,0,15]); % 光伏和风力发电24小时的预测数据 PV=[0,0,0,0,0.0067,0.0133,0.01661,0.0298,0.043,0.0446,0.0529,0.05605,0.0543,0.0558,0.046,0.0394,0.0164,0.0131,0.0099,0.0049,0.0033,0,0,0]; WT=[0.2762,0.2762,0.3077,0.2824,0.3209,0.3275,0.2701,0.2824,0.2348,0.2521,0.2701,0.3275,0.3275,0.3479,0.3209,0.3077,0.3209,0.3342,0.3618,0.334,0.321,0.307,0.3,0.308]; PV_WT=PV(1:T)+WT(1:T); % figure(1); % plot(xx,PV,'-*r',xx,WT,'-b^'); % xlabel('t/h');ylabel('P/kW'); PGTTT=zeros(N,T); PGTTT(N,:)=PV_WT; TTT=T; %_________________________________________________________________________ %_________________________________________________________________________ for iii=1:TTT T=1; B=ones(N,T,PopSize); % B=zeros(N,T,PopSize); % for i=1:PopSize % B(:,:,i)=Y; %生成一个（N,T,PopSize）维的矩阵,生成了PopSize个11*24的矩阵，构成了A % end %随机初始化 Z=zeros(N,T,PopSize); Z(N,:,:)=PV_WT(iii); for jj=1:PopSize; X=zeros(N,T); for i=2:3 for j=1:T X(i,j)=rand(1)*(pmax(i)-pmin(i))+pmin(i); % X(i,j)=Y(i,j)*(rand(1)*(pmax(i)-pmin(i))+pmin(i)); end end X(N,:)=PV_WT(iii); % X(1,:)=PL(1:T)-sum(X(2:4,:)); X(1,:)=PL(iii)-sum(X(2:4,:)); k= X(1,:); % 对第一行平衡机组越限的调整，取其限值 change=k>pmax(1); k(find(change))=pmax(1); change=k<pmin(1); k(find(change))=pmin(1); X(1,:)=k; Z(:,:,jj)=X; end Z; M=Z; % 用在之后第一次迭代更新的时候，如果某个粒子不满足要求，则用初始化时的粒子替代 %适应值函数 Cost=zeros(PopSize,1); for jj=1:PopSize Cost(jj)=shiyingzhi_fadianchengben_2(Z(:,:,jj),iii); % Cost(jj)=shiyingzhi_pollution_2(Z(:,:,jj),iii); % Cost(jj)=shiyingzhi_zonghexiaoyi_2(Z(:,:,jj),iii); end V=rand(4,T,PopSize);%微粒速度随机初始化 %设定当前位置为粒子的最好位置，并记录其最好值。 PBest=Z; FPBest=Cost; %找到初始微粒群体的最好微粒 [Fgbest,r]=min(Cost); CF=Fgbest;%记录当前全局最优值 Best=Z(:,:,r);%用于保存最优粒子的位置 FBest=Fgbest; HengZB=1:MaxIter+1; ZongZB=zeros(1,MaxIter); %------------------------------------------------------------------------------------------------------- iter=0; %循环 while(iter<=MaxIter) iter=iter+1; % 更新惯性权重的值 w_now=((w_start-w_end)*(MaxIter-iter)/MaxIter)+w_end; A=Z(:,:,r); for i=1:PopSize A(:,:,i)=Z(:,:,r);%生成一个（11,24,PopSize）维的矩阵,生成了PopSize个11*24的矩阵，构成了A end %生成随机数 R1=rand(N,T,PopSize); R2=rand(N,T,PopSize); %速度更新 V=w_now*V+c1*R1.*(PBest-Z)+c2*R2.*(A-Z); %对进化后速度小于最小速度的微粒进行处理 changeRows=V>Vmax; V(find(changeRows))=Vmax; changeRows=V<-Vmax; V(find(changeRows))=-Vmax; %微粒位置进行更新 Z=Z+1.0*V.*B; %判断更新后的粒子是否满足约束条件，对不满足约束的粒子进行处理，使其满足约束，如果不满足则重新生成粒子的速度继而更新其位置，直到满足约束条件，如果重复 %次数超过规定的次数，则用原来可行的该粒子代替。 for i=1:PopSize for j=2:3 k=Z(j,:,i); change=k>pmax(j); k(find(change))=pmax(j); change=k<pmin(j); k(find(change))=pmin(j); Z(j,:,i)=k; end Z(N,:,i)=PV_WT(iii); % Z(1,:,i)=PL(1:T)-sum(Z(2:4,:,i)); Z(1,:,i)=PL(iii)-sum(Z(2:4,:,i)); k= Z(1,:,i); % 对第一行平衡机组越限的调整，取其限值 change=k>pmax(1); k(find(change))=pmax(1); change=k<pmin(1); k(find(change))=pmin(1); Z(1,:,i)=k; end %重新计算新位置的适应度值 Cost=zeros(PopSize,1); for jj=1:PopSize Cost(jj)=shiyingzhi_fadianchengben_2(Z(:,:,jj),iii); % Cost(jj)=shiyingzhi_pollution_2(Z(:,:,jj),iii); % Cost(jj)=shiyingzhi_zonghexiaoyi_2(Z(:,:,jj),iii); end %更新每个微粒的最好位置 d=Cost; P=d<FPBest; FPBest(find(P))=d(find(P));%适应值更换 PBest(:,:,find(P))=Z(:,:,find(P));%粒子位置更换 [Fgbest,g]=min(FPBest);%保存最好适应值 Best=PBest(:,:,g);%自己加的一句，因为如果下边的if语句如果不成立的话，Best就不更新了，不合逻辑 if Fgbest<CF%如果本次适应值好于上次则保存 [FBbest,h]=min(FPBest);%最好适应值为FBest Best=PBest(:,:,h); end CF=Fgbest; ZongZB(iter)=CF; end CF; Fgbest; Best; PGTTT(:,iii)=Best; cccc=sum(Best); figure(5); plot(HengZB,ZongZB); xlabel('迭代次数');ylabel('全局最优值/$'); end aaaaaa=sum(PGTTT); Cost=shiyingzhi_fadianchengben_all_2(PGTTT) % Cost=shiyingzhi_pollution_all_2(PGTTT) % Cost=shiyingzhi_zonghexiaoyi_all_2(PGTTT) %__________________________________________________________________________ %__________________________________________________________________________ % figure(2); % plot(HengZB,ZongZB); PGTTT % 输出最优的各类型微电源结果 figure(3); plot(xx,PGTTT(1,:),'-ksquare',xx,PGTTT(2,:),'-r*',xx,PGTTT(3,:),'-b^'); % 绘制24小时的各微电源的功率曲线图 xlabel('t/h');ylabel('P/kW'); axis([1,24,-1,8]); figure(4); plot(xx,sum(Best),':',xx,PL(1:T),'*',xx,sum(Best)-PL(1:T),'x'); toc;