考虑柔性负荷的综合能源系统低碳经济优化调度

目录

1 主要内容

 2 代码

3 程序结果

4 程序链接 

---

1 主要内容

程序主要实现的是考虑柔性负荷的综合能源系统低碳经济优化调度，模型参考《考虑柔性负荷的综合能源系统低碳经济优化调度》，求解方法采用的是混合整数规划算法，通过matlab+cplex进行求解。

结合需求侧柔性负荷的可平移、可转移、可削减特性，构建了含风光储、燃气轮机、柔性负荷等在内的 IES 模型。综合考虑了系统运行成本和碳交易成本，建立了以总成本最低为优化目标的 IES 低碳经济调度模型，并考虑3种场景：

场景 1：考虑包括可平移、可转移、可削减的柔性电负荷和包括可平移、可削减的柔性热负荷参与系统优化调度的情况；

场景 2：只考虑包括可平移、可转移、可削减柔性电负荷，不考虑柔性热负荷参与系统优化调度的情况；

场景 3：不考虑柔性负荷参与系统优化调度的情况。

参考文献见下述两图。

 2 代码

%采用CPIEX求解某微网的运行优化情况，下层优化得出的微网向配电网购电或售电功率，以及各机组的出力%基于能源集线器概念，结合需求侧柔性负荷的可平移、可转移、可削减特性，构建了含风光储、燃气轮机、柔性负荷等%在内的 IES 模型。综合考虑了系统运行成本和碳交易成本，建立了以总成本最低为优化目标的 IES 低碳经济%调度模型，采用cplex求解器对算例进行求解。%场景1 考虑包括可平移、可转移、可削减的柔性电负荷和包括可平移、可削减的柔性热负荷参与系统优化调度的情况；clc;clear;close all;%读取数据 %电负荷、热负荷、光伏、风机、购电价、售电价e_load=[160  150  140  140  130  135  150  180  215  250  275  320  335  290  260  275  270  280  320  360  345  310  220  160];%电负荷h_load=[135  140 150 135 140 120 115 100 115 115 160 180 190 170 140 130 145 200 220 230 160 150 140 130];%热负荷ppv=[0 0 0  0 0  10 15 25 45 75 90 100 80 100 50  40 30 15 10 0 0 0 0 0  ];%光伏预测数据pwt=[60 65  70 75 80 85 90 100 125 150 130 110 100 120 125 130 140 160 180 200 175 160 155 150];%风机预测数据buy_price=[0.25  0.25 0.25 0.25 0.25 0.25 0.25 0.53 0.53 0.53 0.82 0.82 0.82 0.82 0.8 0.53 0.53 0.53 0.82 0.82 0.82 0.53 0.53 0.53];%购电价sell_price=[0.22 0.22 0.22 0.22 0.22 0.22 0.22 0.42 0.42 0.42 0.65 0.65 0.65 0.65 0.65 0.42 0.42 0.42 0.65 0.65 0.65 0.42 0.42 0.42];%售电价%需求响应数据Pcut=[10 10 10 10 10 10 15 15 25 50 50 50 50 50 50 50 50 50 50 50 40 40 15 10];%可削减电负荷Temp_Pcut=binvar(1,24,'full'); % 电负荷削减标志PPcut=sdpvar(1,24,'full');%电负荷消减量n1=zeros(1,1);%消减连续Hcut=[25 25 25 25 25 25 25 25 30 40 40 40 40 40 40 40 40 40 50 50 30 30 20 15];%可削减热负荷Temp_Hcut=binvar(1,24,'full'); % 热负荷削减标志HHcut=sdpvar(1,24,'full');%热负荷消减量n2=zeros(1,1);%消减连续​Ptran=[0 0 0 0 0 0 0 0 0 0 0 0 25 25 25 25 0 0 0 0 0 0 0 0 ];%可转移电负荷Temp_Ptran=binvar(1,24,'full'); % 可转移电负荷 转移标志PPtran=sdpvar(1,24,'full');%电负荷转移量​Pshift1=[0 0 0 0 0 0 0 0 0 0 0 25 25 0 0 0 0 0 0 0 0 0 0 0 ];%可平移电负荷1Temp_Pshift1=binvar(1,24,'full'); % 可平移电负荷1 平移标志PPshift1=sdpvar(1,24,'full');%可平移电负荷1量Pshift2=[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0  25 25 25 0 0 ];%可平移电负荷2Temp_Pshift2=binvar(1,24,'full'); % 可平移电负荷2 平移标志PPshift2=sdpvar(1,24,'full');%可平移电负荷2量Hshift=[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 45 45 45 0 0 0 0 ];%可平移热负荷Temp_Hshift=binvar(1,24,'full'); % 可平移热负荷 平移标志HHshift=sdpvar(1,24,'full');%可平移热负荷量​for i=1:24    Pfix(i)=e_load(i)-Pshift1(i)-Pshift2(i)-Ptran(i)-Pcut(i);%基础电负荷endfor i=1:24    Hfix(i)=h_load(i)-Hshift(i)-Hcut(i);%基础热负荷end​​%定义机组变量P_pv=sdpvar(1,24,'full');%光伏电输出功率P_wt=sdpvar(1,24,'full');%风机电输出功率P_mt=sdpvar(1,24,'full');%燃气轮机电输出功率P_GB=sdpvar(1,24,'full');%燃气锅炉输出热功率​Pbuy=sdpvar(1,24,'full');%从电网购电电量Psell=sdpvar(1,24,'full');%向电网售电电量Pnet=sdpvar(1,24,'full');%与电网交换功率Temp_net=binvar(1,24,'full'); % 购|售电标志​Pcharge=sdpvar(1,24,'full');%充电功率UPcharge=binvar(1,24,'full');%充电标志  Pdischarge=sdpvar(1,24,'full');%放电功率UPdischarge=binvar(1,24,'full');%放电标志  B=sdpvar(1,24,'full');%电储能余量​Hcharge=sdpvar(1,24,'full');%储热系统充热Hdischarge=sdpvar(1,24,'full');%储热系统放热UHcharge=binvar(1,24,'full'); %储热系统充热标志UHdischarge=binvar(1,24,'full'); %储热系统放热标志H=sdpvar(1,24,'full'); %热储能余量​​%储能参数%电储能参数E_storage_max=0.95*100;E_storage_min=0.4*100;e_loss=0.001;e_charge=0.9;e_discharge=0.9;%电储能容量/自损/充电/放电%热储能参数H_storage_max=0.95*100;H_storage_min=0.4*100;h_loss=0.001;h_charge=0.9;h_discharge=0.9;%热储能容量//自损/充热/放热%约束条件Constraints =[]; %% 电储能容量约束、SOC约束、充电约束、放电约束、充放电状态约束、爬坡约束B(1,1)=E_storage_min;%电储能初始 for t=2:25  %在一个周期内的充放电功率    Constraints=[Constraints,(B(mod(t-1,24)+1)==(B(mod(t-2,24)+1)*(1-e_loss)+(e_charge*Pcharge(mod(t-2,24)+1)-(1/e_discharge)*Pdischarge(mod(t-2,24)+1))))]; end% % %   %全周期净交换功率为零%     Constraints=[Constraints,B(1,24)==E_storage_min];%初始功率相等即可for i=1:24Constraints=[Constraints,E_storage_min<=B(1,i)<=E_storage_max];%容量约束限制end for i=1:24     Constraints=[Constraints,30*UPcharge(1,i)<=Pcharge(1,i)<=40*UPcharge(1,i)];%电储能充电约束     Constraints=[Constraints,30*UPdischarge(1,i)<=Pdischarge(1,i)<=40*UPdischarge(1,i)];%电储能放电约束 end %蓄电池充放电约束 for i=1:24     Constraints=[Constraints,UPcharge(1,i)+UPdischarge(1,i)<=1];   %不同时充放电  end   Constraints=[Constraints,sum(UPcharge(1,1:24))+sum(UPdischarge(1,1:24))==16];%使用寿命小于24​ %% 热储能容量约束、SOC约束、充热约束、放热约束、充放热状态约束H(1,1)=H_storage_min;%热储能初始 for t=2:25  %在一个周期内的充放热功率    Constraints=[Constraints,(H(mod(t-1,24)+1)==(H(mod(t-2,24)+1)*(1-h_loss)+(h_charge*Hcharge(mod(t-2,24)+1)-(1/h_discharge)*Hdischarge(mod(t-2,24)+1))))]; end% %  %全周期净交换功率为零%    Constraints=[Constraints,H(1,24)==H_storage_min];%初始功率相等即可for i=1:24Constraints=[Constraints,H_storage_min<=H(1,i)<=H_storage_max];%容量约束限制end for i=1:24     Constraints=[Constraints,5*UHcharge(1,i)<=Hcharge(1,i)<=30*UHcharge(1,i)];%热储能充电约束     Constraints=[Constraints,5*UHdischarge(1,i)<=Hdischarge(1,i)<=30*UHdischarge(1,i)];%热储能放电约束 end %蓄热池充放电约束 for i=1:24     Constraints=[Constraints,UHcharge(1,i)+UHdischarge(1,i)<=1];   %不同时充放热  end   Constraints=[Constraints,sum(UHcharge(1,1:24))+sum(UHdischarge(1,1:24))==16];%使用寿命小于24​%% 机组约束for i=1:24   Constraints = [Constraints,0<=P_pv(i)<=ppv(i)];%光伏上下限约束    Constraints = [Constraints,0<=P_wt(i)<=pwt(i)];%风机上下限约束   Constraints = [Constraints,0<=P_mt(i)<=65];%燃气轮机上下限约束   Constraints = [Constraints,0<=P_GB(i)<=160];%燃气锅炉上下限约束   Constraints = [Constraints, -160<=Pnet(i)<=160,0<=Pbuy(i)<=160, -160<=Psell(i)<=0]; %主网功率交换约束   Constraints = [Constraints, implies(Temp_net(i),[Pnet(i)>=0,Pbuy(i)==Pnet(i),Psell(i)==0])]; %购电情况约束   Constraints = [Constraints, implies(1-Temp_net(i),[Pnet(i)<=0,Psell(i)==Pnet(i),Pbuy(i)==0])]; %售电情况约束 end  %% 需求响应约束%% 可平移电负荷1量    Constraints= [Constraints,sum(Temp_Pshift1(1,1:24)) == 2,sum(Temp_Pshift1(1,5:21)) == 2];%可平移电负荷1 平移标志    for i=5:20 %时段区间为5~21-2+1   Constraints = [Constraints,sum(Temp_Pshift1(1,i:i+1)) >= 2*(Temp_Pshift1(1,i)-Temp_Pshift1(1,i-1))];%连续2个时段    end    for i=1:24       Constraints = [Constraints,PPshift1(1,i)== 25*Temp_Pshift1(1,i)];%可平移电负荷1量    end%% 可平移电负荷2量        Constraints = [Constraints,sum(Temp_Pshift2(1,1:24)) == 3,sum(Temp_Pshift2(1,7:23)) == 3];%可平移电负荷2 平移标志    for i=7:21 %时段区间为7~23-3+1    Constraints = [Constraints,sum(Temp_Pshift2(1,i:i+2)) >= 3*(Temp_Pshift2(1,i)-Temp_Pshift2(1,i-1)-Temp_Pshift2(1,i-2))];%连续3个时段    end       for i=1:24       Constraints = [Constraints,PPshift2(1,i)== 25*Temp_Pshift2(1,i)];%可平移电负荷2量       end%% 可平移热负荷量        Constraints = [Constraints,sum(Temp_Hshift(1,1:24)) == 3,sum(Temp_Hshift(1,5:21)) == 3];%可平移热负荷 平移标志​    for i=5:19%时段区间为5~21-3+1    Constraints = [Constraints,sum(Temp_Hshift(1,i:i+2)) >= 3*(Temp_Hshift(1,i)-Temp_Hshift(1,i-1))];%连续3个时段    end    for i=1:24       Constraints = [Constraints,HHshift(1,i)== 45*Temp_Hshift(1,i)];%可平移电负荷2量    end       %% 可转移电负荷(大于5自然会大于2)  for i=1:24      Constraints = [Constraints,Temp_Ptran(i)*8<=PPtran(i)<=Temp_Ptran(i)*26.7 ];%可转移电负荷  end      Constraints = [Constraints,sum(Temp_Ptran(1,1:24)) == 5,sum(Temp_Ptran(1,4:22)) ==5];%可转移电负荷            Constraints = [Constraints,sum(Temp_Ptran(1,1:24)) ==5];%可转移电负荷    for i=4:18 %时段区间为4~22-5+1    Constraints = [Constraints,sum(Temp_Ptran(1,i:i+4)) >= 5*(Temp_Ptran(1,i)-Temp_Ptran(1,i-1))];    end​%% 可削减电负荷​Constraints=[Constraints,sum(Temp_Pcut)==8,sum(Temp_Pcut(1,5:22))==8];Constraints=[Constraints,2<=n1<=5];    for i=5:22-n1+1 %时段区间为5~22-n1+1    Constraints = [Constraints,sum(Temp_Pcut(1,i:i+n1-1)) >= n1*(Temp_Pcut(1,i)-Temp_Pcut(1,i-1))];    endfor i=1:24       Constraints = [Constraints,0<=PPcut(1,i)<=Temp_Pcut(1,i)*0.9*Pcut(i)];%可消减电负荷end%% 可削减热负荷Constraints=[Constraints,sum(Temp_Hcut(1,1:24))==8,sum(Temp_Hcut(1,11:19))==8];Constraints=[Constraints,2<=n2<=5];    for i=11:19-n2+1 %时段区间为11~19-n2+1    Constraints = [Constraints,sum(Temp_Hcut(1,i:i+n1-1)) >= n1*(Temp_Hcut(1,i)-Temp_Hcut(1,i-1))];    endfor i=1:24       Constraints = [Constraints,Temp_Hcut(1,i)*0<=HHcut(1,i)<=Temp_Hcut(1,i)*0.9*Hcut(i)];%可消减热负荷end%% 电平衡   for i=1:24          Constraints = [Constraints,P_mt(i)+P_pv(i)+P_wt(i)+Pnet(i)-Pcharge(1,i)+Pdischarge(1,i)==Pfix(i)+Pcut(i)+PPshift1(i)+PPshift2(i)+PPtran(i)-PPcut(i)]; %电平衡约束   Constraints = [Constraints,P_GB(i)+0.83*P_mt(i)/0.45-Hcharge(1,i)+Hdischarge(1,i)==Hfix(i)+Hcut(i)+HHshift(i)-HHcut(i)]; %热平衡约束   end      %% 目标函数%% 从大电网的购电成本C_gridbuy=0;for i=1:24    C_gridbuy=C_gridbuy+Pbuy(i)*buy_price(i);end%% 向大电网的售电成本C_gridsell=0;for i=1:24    C_gridsell=C_gridsell+Psell(i)*sell_price(i);end%运行成本C_OM=0;for i=1:24 C_OM=C_OM+0.72*P_pv(i)+0.52*P_wt(i);%风机光伏运维成本end​%% 燃料成本C_fuel=0;for i=1:24 C_fuel=C_fuel+2.5*P_GB(i)/9.7+2.5*P_mt(i)/0.45/9.7;%耗气成本end%% 储能运行成本C_storge=0;for i=1:24 C_storge=C_storge+0.5*(Pcharge(i)+Pdischarge(i)+Hcharge(i)+Hdischarge(i));%储能运行成本end​%% 补偿成本C_L=0;for i=1:24    C_L=C_L+0.2*(PPshift1(i)+PPshift2(i))+0.1*HHshift(i)+0.3*PPtran(i)+0.4*PPcut(i)+0.2*HHcut(i);end%% 碳交易成本​Q_carbon=0;%碳排放量-碳配额量(克)for i=1:24    Q_carbon=Q_carbon+(((1303-798)*(Pbuy(i)+abs(Psell(i)))+(564.7-424)*(P_GB(i)/9.7+P_mt(i)/0.45/9.7)+...        (43-78)*P_wt(i)+(154.5-78)*P_pv(i)+91.3*(Pcharge(i)+Pdischarge(i))));end​E_v=sdpvar(1,5);%每段区间内的长度,分为5段,每段长度是2000lamda=0.15*10^(-3);%碳交易基价Constraints=[Constraints,   Q_carbon==sum(E_v),%总长度等于Q_carbon   0<=E_v(1:4)<=120000,%除了最后一段，每段区间长度小于等于120000g   0<=E_v(5),  ];%碳交易成本C_CO2=0;for v=1:5    C_CO2=C_CO2+(lamda+(v-1)*0.25*lamda)*E_v(v);end​​F= C_OM+C_fuel+C_gridbuy+C_gridsell+C_storge+C_L+C_CO2;ops = sdpsettings('solver','cplex', 'verbose', 2);%参数指定程序用cplex求解器optimize(Constraints,F,ops)

3 程序结果

 

 

 

  

上面是场景1出图12张，场景2和3也有全套的出图，学习必备佳品！