目录
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 ];%可平移电负荷1 Temp_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 ];%可平移电负荷2 Temp_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:24Pfix(i)=e_load(i)-Pshift1(i)-Pshift2(i)-Ptran(i)-Pcut(i);%基础电负荷 end for i=1:24Hfix(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:24 Constraints=[Constraints,E_storage_min<=B(1,i)<=E_storage_max];%容量约束限制 endfor i=1:24Constraints=[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:24Constraints=[Constraints,UPcharge(1,i)+UPdischarge(1,i)<=1]; %不同时充放电 endConstraints=[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:24 Constraints=[Constraints,H_storage_min<=H(1,i)<=H_storage_max];%容量约束限制 endfor i=1:24Constraints=[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:24Constraints=[Constraints,UHcharge(1,i)+UHdischarge(1,i)<=1]; %不同时充放热 endConstraints=[Constraints,sum(UHcharge(1,1:24))+sum(UHdischarge(1,1:24))==16];%使用寿命小于24 %% 机组约束 for i=1:24Constraints = [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+1Constraints = [Constraints,sum(Temp_Pshift1(1,i:i+1)) >= 2*(Temp_Pshift1(1,i)-Temp_Pshift1(1,i-1))];%连续2个时段endfor i=1:24Constraints = [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+1Constraints = [Constraints,sum(Temp_Pshift2(1,i:i+2)) >= 3*(Temp_Pshift2(1,i)-Temp_Pshift2(1,i-1)-Temp_Pshift2(1,i-2))];%连续3个时段endfor i=1:24Constraints = [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+1Constraints = [Constraints,sum(Temp_Hshift(1,i:i+2)) >= 3*(Temp_Hshift(1,i)-Temp_Hshift(1,i-1))];%连续3个时段endfor i=1:24Constraints = [Constraints,HHshift(1,i)== 45*Temp_Hshift(1,i)];%可平移电负荷2量end %% 可转移电负荷(大于5自然会大于2)for i=1:24Constraints = [Constraints,Temp_Ptran(i)*8<=PPtran(i)<=Temp_Ptran(i)*26.7 ];%可转移电负荷endConstraints = [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+1Constraints = [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+1Constraints = [Constraints,sum(Temp_Pcut(1,i:i+n1-1)) >= n1*(Temp_Pcut(1,i)-Temp_Pcut(1,i-1))];end for i=1:24Constraints = [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+1Constraints = [Constraints,sum(Temp_Hcut(1,i:i+n1-1)) >= n1*(Temp_Hcut(1,i)-Temp_Hcut(1,i-1))];end for i=1:24Constraints = [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:24C_gridbuy=C_gridbuy+Pbuy(i)*buy_price(i); end %% 向大电网的售电成本 C_gridsell=0; for i=1:24C_gridsell=C_gridsell+Psell(i)*sell_price(i); end %运行成本 C_OM=0; for i=1:24C_OM=C_OM+0.72*P_pv(i)+0.52*P_wt(i);%风机光伏运维成本 end %% 燃料成本 C_fuel=0; for i=1:24C_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:24C_storge=C_storge+0.5*(Pcharge(i)+Pdischarge(i)+Hcharge(i)+Hdischarge(i));%储能运行成本 end %% 补偿成本 C_L=0; for i=1:24C_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:24Q_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段,每段长度是2000 lamda=0.15*10^(-3);%碳交易基价 Constraints=[Constraints,Q_carbon==sum(E_v),%总长度等于Q_carbon0<=E_v(1:4)<=120000,%除了最后一段,每段区间长度小于等于120000g0<=E_v(5),]; %碳交易成本 C_CO2=0; for v=1:5C_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也有全套的出图,学习必备佳品!