%参数记录
clear all;
clc;
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)*D+6)<0.5*(max(-SBatnom*RBatdis,
