基于多能互补的热电联供型微网优化运行（matlab代码）

%%  赋值
MT=intvar(1,24,'full');
Pgrid=intvar(1,24,'full');  %  购 、 售
Pnas=intvar(2,24,'full'); % 充、 放
H=intvar(1,24,'full');%锅炉
Q=intvar(1,24,'full');%可中断电负荷
Hti=intvar(1,24,'full');%充热
Hto=intvar(1,24,'full');%放热
xx=intvar(1,24,'full');%%%%%时平负荷量
yy=intvar(1,24,'full');
%% 0-1赋值
I_Hti=binvar(1,24,'full');%充热
I_Hto=binvar(1,24,'full');%放热
I_MT=binvar(1,24,'full');
I_Pnas=binvar(2,24,'full');% 1运行 0停止
I_Q=binvar(1,24,'full');
%%   目标函数
for i=1:24%发电成本
    Cf(1,i)=Cgas(i)*(aF(1)*MT(1,i)+bF(1)*I_MT(1,i));
end
for i=1:24%%%%余热回收
    H_cycle(1,i)=aH(1)*MT(1,i)+bH(1)*I_MT(1,i);
end
for k=1:24 %PCC交互成本 % 1-5，23-24 谷 % 6-12，19-22 峰 % 13-18 平 
    if k>=1&&k<7
        Cgrid(1,k)=Pgrid(1,k).*buy(3);
    elseif k>=7&&k<13
        Cgrid(1,k)=Pgrid(1,k).*buy(1);
    elseif k>=13&&k<19
        Cgrid(1,k)=Pgrid(1,k).*buy(2);
    elseif k>=19&&k<23
        Cgrid(1,k)=Pgrid(1,k).*buy(1);
    else
        Cgrid(1,k)=Pgrid(1,k).*buy(3);
    end
end
for k=1:24 %需求响应单位成本
   if k>=1&&k<7
        bu_q(1,k)=0.9*buy(1);
        bu_x(1,k)=0.5*buy(1);
        bu_p(1,k)=0.6*Cgas(1,k);
    elseif k>=7&&k<13
        bu_q(1,k)=0.9*buy(1);
        bu_x(1,k)=0.5*buy(1);
        bu_p(1,k)=0.6*Cgas(1,k);
    elseif k>=13&&k<19
        bu_q(1,k)=0.9*buy(1);
        bu_x(1,k)=0.5*buy(1);
        bu_p(1,k)=0.6*Cgas(1,k);
    elseif k>=19&&k<23
        bu_q(1,k)=0.9*buy(1);
        bu_x(1,k)=0.5*buy(1);
        bu_p(1,k)=0.6*Cgas(1,k);
    else
        bu_q(1,k)=0.8*buy(1);
        bu_x(1,k)=buy(1)/2;
        bu_p(1,k)=0.6*Cgas(k);
    end
end
for k=1:24  %% 切除负荷成本
        if k>=7&&k<=12
          Ck1(1,k)=(Q(1,k).*bu_q(1,k))+xx(1,k)*bu_x(1,k);
    elseif  k>=19&&k<=20
          Ck1(1,k)=(Q(1,k).*bu_q(1,k))+xx(1,k)*bu_x(1,k);
       else
          Ck1(1,k)=(Q(1,k).*bu_q(1,k))+xx(1,k)*bu_x(1,k);
        end
end
for i=1:24%锅炉成本
    Ch(1,i)=Cgas(i)*(H(1,i))/LHV;
end
F=0;%目标函数
mm=3.1;
for k=1:24 %1.8
     F=F+Cf(1,k)+Cgrid(1,k)+Ch(1,k)+(Pnas(1,k)+Pnas(2,k))*0.024;
end
for k=1:24 %SOC值
    SOC(k)=(500+sum(Pnas(1,1:k).*I_Pnas(1,1:k)-(Pnas(2,1:k)).*I_Pnas(2,1:k)))/1000;   
end
begin=500;
% for i=1:24%%热储能
%     L(1,i)=begin*h_n+h_charge*Hti(1,i)-Hto(1,i);%%%热储能容量
%     begin=L(1,i);
% end
%%  约束条件
constraints=[];
%% 状态约束
for k=1:24  %Pgrid状态 %Pnas状态
    constraints=[constraints,I_Pnas(1,k)+I_Pnas(2,k)<=1];
%     constraints=[constraints,I_Hti(1,k)+I_Hto(1,k)<=1];
end
constraints=[constraints,sum(I_Pnas(1,1:24)+I_Pnas(2,1:24))<=14];
%% 上下限约束
for k=1:24     
    constraints=[constraints,25.*I_MT(1,k)<=MT(1,k)<=145.*I_MT(1,k)]; 
    
    constraints=[constraints,Pgrid_min<=Pgrid(1,k)<=Pgrid_max]; 
    constraints=[constraints,Pnas_min.*I_Pnas(1,k)<=Pnas(1,k)<=Pnas_max.*I_Pnas(1,k)]; 
    constraints=[constraints,Pnas_min.*I_Pnas(2,k)<=Pnas(2,k)<=Pnas_max.*I_Pnas(2,k)];    
end
%%% MT爬坡率
for i=1:23
    constraints=[constraints,-55<=MT(1,i+1)-MT(1,i)<=65];
end
%   PCC深度限制
 for k=1:23
     constraints=[constraints,-90<=Pgrid(1,k+1)-Pgrid(1,k)<=90];
 end
 %%荷电状态
 for k=1:24
      constraints=[constraints,SOC_min<=300+sum(Pnas(1,1:k)-Pnas(2,1:k))<=SOC_max];      
 end
 %%储能充放电深度限制
 for k=1:23
     constraints=[constraints,-50<=Pnas(1,k+1)-Pnas(2,k+1)-Pnas(1,k)+Pnas(2,k)<=50];
 end
     constraints=[constraints,sum(Pnas(1,1:24))==sum(Pnas(2,1:24))];
%%%锅炉上下限、爬坡率
for i=1:24
    constraints=[constraints,30<=H(1,i)<=H_max];
end
for i=1:23   
   constraints=[constraints,-90<=H(1,i+1)-H(1,i)<=90];
end
%%功率平衡
for k=1:24
constraints=[constraints,MT(1,k)-Pnas(1,k)+Pnas(2,k)+Pgrid(1,k)+Pwind(1,k)==load(k)+L5(1,k)+L6(1,k)+L7(1,k)];
end