matlab-基于多目标算法的冷热电联供型综合能源系统运行优化-源码

clc; clear; close all; warning off; addpath(genpath(pwd)); tic %tic用来保存当前时间，而后使用toc来记录程序完成时间 %---燃气发电机参数---% PGmax=[1200,2800]; %燃气发电机最大电功率 PGmin=[300,300]; %燃气发电机最小电功率 PGnom=[1200,2800]; %燃气发电机额定电功率 Aa=8.935;Bb=33.157;Cc=-27.081;Dd=17.989; %这里的Aa,Bb,Cc,Dd分别对应四台不同的燃气发电机的发电效率 %-----------% Pgridmax=600; Pgridmin=400; % Pgridmin=0; % ----------------蓄电池组参数-------------------% SBatmax=180; %蓄电池保持稳定运行时的最大容量 SBatmin=40; %蓄电池保持稳定运行时的最小容量 SBatnom=200; %蓄电池在协调周期内的额定容量 RBatcha=0.20;%蓄电池组l在t时段的最大充电率 RBatdis=0.40;%蓄电池组l在t时段的最大放电率 dBat=0.04; %蓄电池组的自放电率 yBatcha=0.95;yBatdis=0.95; %充放电效率 %------------------------------------------------% % ---燃气锅炉参数---% Qbolnom=1000; %燃气锅炉额定功率 QHrsnom=1000; %热回收系统的额定功率 QAcnom=400; %吸收冷机的额定功率 QEcnom=300; %电制冷机的额定功率 yHrs=0.60; %热回收效率 Copac=0.80;%吸收制冷能效比 Copec=3; %电制冷能效比 % % QDtherm=[1300,1400,1360,1700,1600,2000,3100,3450,3900,4400,4300,4500,5400,5000,5100,4950,4600,4710,4400,4000,3650,3570,2900,2420]; %系统在一天中需要的热负荷 QDcool=[2350,2800,2850,3400,3800,4000,4700,4350,5100,7500,7600,8450,8700,8050,7700,7450,7200,6650,7300,5700,5100,4450,3600,3080]; %系统在一天中需要的冷负荷 w=0.7; c1=2; %加速因子 c2=2; N=100; %粒子个数 SBat=zeros(N,25);QAc=zeros(N,24); QEc=zeros(N,24);QHrs=zeros(N,24);yg=zeros(N,48); k=1;D=8; Tmax=1000; %最大迭代次数 %-----------% for i=1:N SBat(i,19)=63.5209380080000;%初始储能量，单位为kW.h, end for i=1:N %分时段初始化满足约束的初始值 %----------% for t=1:24 %-----------% %燃气发电机出力初始化 for j=1:2 X(i,(t-1)*D+j)=unifrnd(PGmin(j),PGmax(j)); V(i,(t-1)*D+j)=unifrnd(PGmin(j)-PGmax(j),PGmax(j)-PGmin(j)); end %锅炉提供的热功率用于热负荷部分量的初始化 for j=1:2 yg(i,(t-1)*2+j)=(Aa+Bb*(X(i,(t-1)*D+j)/PGnom(j))+Cc*(X(i,(t-1)*D+j)/PGnom(j))^2+Dd*(X(i,(t-1)*D+j)/PGnom(j))^3)/100; end QHrs(i,t)=(X(i,(t-1)*D+1)*(1-yg(i,(t-1)*2+1))+X(i,(t-1)*D+2)*(1-yg(i,(t-1)*2+2)))*yHrs; %热回收系统在单时段提供的总热功率 if QHrs(i,t)>1000 QHrs(i,t)=1000; end if QHrs(i,t)-QDtherm(t)>=0 X(i,(t-1)*D+3)=0; V(i,(t-1)*D+3)=0; else X(i,(t-1)*D+3)=unifrnd(abs(min(0,QHrs(i,t)-QDtherm(t))),QDtherm(t)); V(i,(t-1)*D+3)=unifrnd(abs(min(0,QHrs(i,t)-QDtherm(t)))-QDtherm(t),QDtherm(t)-abs(min(0,QHrs(i,t)-QDtherm(t)))); end %锅炉提供的热功率用于冷负荷部分量的初始化 X(i,(t-1)*D+4)=unifrnd(0,min((Qbolnom-X(i,(t-1)*D+3)),QDcool(t)/Copac)); V(i,(t-1)*D+4)=unifrnd(-min((Qbolnom-X(i,(t-1)*D+3)),QDcool(t)/Copac),min((Qbolnom-X(i,(t-1)*D+3)),QDcool(t)/Copac)); %微网在单个时段内与电网交互的电量初始化 X(i,(t-1)*D+5)=unifrnd(-Pgridmin,Pgridmax); V(i,(t-1)*D+5)=unifrnd(-Pgridmin-Pgridmax,Pgridmax+Pgridmin); %蓄电池组的充放电功率初始化 X(i,(t-1)*D+6)=unifrnd(max(-SBatnom*RBatdis,SBatmin-SBat(i,t)),min(SBatnom*RBatcha,SBatmax-SBat(i,t))); V(i,(t-1)*D+6)=unifrnd(max(-SBatnom*RBatdis,SBatmin-SBat(i,t))-min(SBatnom*RBatcha,SBatmax-SBat(i,t)),min(SBatnom*RBatcha,SBatmax-SBat(i,t))-max(-SBatnom*RBatdis,SBatmin-SBat(i,t))); %热回收系统用于提供冷功率的热功率初始化 for j=1:2 yg(i,(t-1)*2+j)=(Aa+Bb*(X(i,(t-1)*D+j)/PGnom(j))+Cc*(X(i,(t-1)*D+j)/PGnom(j))^2+Dd*(X(i,(t-1)*D+j)/PGnom(j))^3)/100; end QHrs(i,t)=(X(i,(t-1)*D+1)*(1-yg(i,(t-1)*2+1))+X(i,(t-1)*D+2)*(1-yg(i,(t-1)*2+2)))*yHrs; %热回收系统在单时段提供的总热功率 if QHrs(i,t)>QHrsnom QHrs(i,t)=QHrsnom; end X(i,(t-1)*D+7)=min(max(0,QHrs(i,t)-QDtherm(t)),QDcool(t)/Copac); V(i,(t-1)*D+7)=0; %电制冷机制冷所消耗的电功率初始化 QAc(i,t)=X(i,(t-1)*D+7)*Copac; %吸收制冷机在单时段提供的冷功率 X(i,(t-1)*D+8)=max(0,min((QDcool(t)-QAc(i,t))/Copec,2*QEcnom/Copec)); V(i,(t-1)*D+8)=0; X(i,(t-1)*D+4)=abs(min(0,(QAc(i,t)+ X(i,(t-1)*D+8)*Copec-QDcool(t)))); V(i,(t-1)*D+4)=0; %--------------------% if X(i,(t-1)*D+6)>=0 SBat(i,t+1)= SBat(i,t)*(1-dBat)+X(i,(t-1)*D+6)*yBatcha; %蓄电池组在单时段时的剩余容量 else SBat(i,t+1)= SBat(i,t)*(1-dBat)+X(i,(t-1)*D+6)/yBatdis; %蓄电池组在单时段时的剩余容量 end if SBat(i,t+1)>SBatmax SBat(i,t+1)=SBatmax; end if SBat(i,t+1)<SBatmin SBat(i,t+1)=SBatmin; end end end WEI=24*8; for i=1:N p(i)=fitness(X(i,:),WEI); y(i,:)=X(i,:); end pg = X(N,:); %Pg为全局最优 for i=1:(N-1) if fitness(X(i,:),WEI)<fitness(pg,WEI) pg=X(i,:); end end for K=1:Tmax if ~mod(K,10) %K与10作除法后的余数 fprintf('current iter:%d\n',K) end for i=1:N V(i,1:WEI)=w*V(i,1:WEI)+c1*rand*(y(i,1:WEI)-X(i,1:WEI))+c2*rand*(pg-X(i,1:WEI)); %粒子速度更新 %燃气发电机出力速度越限处理 %-----------% for t=1:24 %---------------% for j=1:2 if V(i,(t-1)*D+j)>0.5*(PGmax(j)-PGmin(j)) V(i,(t-1)*D+j)=0.5*(PGmax(j)-PGmin(j)); end if V(i,(t-1)*D+j)<0.5*(PGmin(j)-PGmax(j)) V(i,(t-1)*D+j)=0.5*(PGmin(j)-PGmax(j)); end end %锅炉提供的热功率用于热负荷部分量的速度越限处理 for j=1:2 yg(i,(t-1)*2+j)=(Aa+Bb*(X(i,(t-1)*D+j)/PGnom(j))+Cc*(X(i,(t-1)*D+j)/PGnom(j))^2+Dd*(X(i,(t-1)*D+j)/PGnom(j))^3)/100; end QHrs(i,t)=(X(i,(t-1)*D+1)*(1-yg(i,(t-1)*2+1))+X(i,(t-1)*D+2)*(1-yg(i,(t-1)*2+2)))*yHrs; %热回收系统在单时段提供的总热功率 if QHrs(i,t)>QHrsnom QHrs(i,t)=QHrsnom; end if QHrs(i,t)-QDtherm(t)>=0 if V(i,(t-1)*D+3)>0 V(i,(t-1)*D+3)=0; end if V(i,(t-1)*D+3)<0 V(i,(t-1)*D+3)=0; end else if V(i,(t-1)*D+3)>0.5*(QDtherm(t)-abs(min(0,QHrs(i,t)-QDtherm(t)))) V(i,(t-1)*D+3)=0.5*(QDtherm(t)-abs(min(0,QHrs(i,t)-QDtherm(t)))); end if V(i,(t-1)*D+3)<0.5*(abs(min(0,QHrs(i,t)-QDtherm(t)))-QDtherm(t)) V(i,(t-1)*D+3)=0.5*(abs(min(0,QHrs(i,t)-QDtherm(t)))-QDtherm(t)); end end %锅炉提供的热功率用于冷负荷部分量的速度越限处理 if V(i,(t-1)*D+4)>0 V(i,(t-1)*D+4)=0; end if V(i,(t-1)*D+4)<0 V(i,(t-1)*D+4)=0; end %微网在单个时段与电网交互的电量的速度越限处理 if V(i,(t-1)*D+5)>0.5*(Pgridmax+Pgridmin) V(i,(t-1)*D+5)=0.5*(Pgridmax+Pgridmin); end if V(i,(t-1)*D+5)<0.5*(-Pgridmin-Pgridmax) V(i,(t-1)*D+5)=0.5*(-Pgridmin-Pgridmax); end %蓄电池组的充放电功率的速度越限处理 if V(i,(t-1)*D+6)>0.5*(min(SBatnom*RBatcha,SBatmax-SBat(i,t))-max(-SBatnom*RBatdis,SBatmin-SBat(i,t))) V(i,(t-1)*D+6)=0.5*(min(SBatnom*RBatcha,SBatmax-SBat(i,t))-max(-SBatnom*RBatdis,SBatmin-SBat(i,t))); end if V(i,(t-1