扫码关注
参考文献：《碳交易机制下考虑需求响应的综合能源系统优化运行》摘 要:综合能源系统是实现“双碳”目标的有效途径,为进一步挖掘其需求侧可调节潜力对碳减排的作用,提出了一种碳交易机制下考虑需求响应的综合能源系统优化运行模型。首先,根据负荷响应特性将需求响应分为价格型和替代型 2 类,分别建立了基于价格弹性矩阵的价格型需求响应模型,及考虑用能侧电能和热能相互转换的替代型需求响应模型:其次,采用基准线法为系统无偿分配碳排放配额,并考虑燃气轮机和燃气锅炉的实际碳排放量,构建一种面向综合能源系统的碳交易机制;最后,以购能成本、碳交易成本及运维成本之和最小为目标函数,建立综合能源系统低碳优化运行模型,并通过4 类典型场景对所提模型的有效性进行了验证。通过对需求响应灵敏度、气轮机热分配比例和不同碳交易价格下系统的运行状态分析发现，合理分配价格型和替代型需求响应及燃气轮机产热比例有利于提高系统运行经济性,制定合理的碳交易价格可以实现系统经济性和低碳性协同。关键词:碳交易机制;需求响应;综合能源系统;优化运行1 背景
2019 年，电力行业二氧化碳排放占全国碳排放总量超过 40%。2020 年9 月,我国提出争取 2060 年前实现碳中和的目标,发展低碳电力迫在眉睫。目前,我国正在逐步推行碳交易市场，努力通过市场手段实现二氧化碳“零排放”。IEHS 能够实现电能、热能的互补协同,提高能源利用效率,满足用户多种能源梯级利用的同时保障持续可靠供能。本文构建了含需求响应的 IEHS 架构,如图1 所示。在该系统中,电能和气能分别由上级电网、气网供应，从上级气网购气用来供给热电联产装置(combined heat and power，CHP)和燃气锅炉(gasboiler，GB),剩余电能可出售给上级电网;能量耦合设备有 CHP、热泵(heat pump，HP) 和 GB,能实现电热能量双向流动;CHP 由燃气轮机(gas turbine,GT)、余热锅炉(waste heat boiler，WHB)和基于有机朗肯循环(organic Rankine cycle，ORC)的低温余热发电装置组成,运行方式为热电解耦,该运行方式能适应系统不同运行工况;HP 和 GB 消纳风电并承担部分热负荷。引入 DR 可以平抑负荷曲线波动,实现电热的交互耦合、削峰填谷并降低运行成本。
2 IEHS 优化运行模型
碳交易机制下考虑 DR 的 IEHS 优化运行模型旨在满足系统运行约束的前提下,实现整个网络经济性最佳。以购能成本、碳交易成本及运维成本之和最小为目标函数:
碳交易机制下考虑 DR 的 IEHS 优化运行约束有:风电出力约束、能量平衡约束、设备能量转换约束、储能设备约束和用户用电方式满意度约束。
本文所求问题为混合整数线性规划问题,首先分析价格型需求响应和替代型需求响应,得到需求响应后的负荷曲线;然后,引入碳交易机制，并将碳交易机制下的碳交易成本作为目标函数的组成部分;最后在满足风电出力约束、能量平衡约束、设备能量转换约束、储能设备约束和用户用电方式满意度约束的条件下,基于 MATLAB 平台调用 CPLEX 求解器求解求解流程如图 2 所示。
3 算例仿真
以北方某工业园区为研究对象，以 24 h 为一个运行周期,单位运行时间为 1 h。系统中已安装设备有由 GT、WHB 和基于 ORC 的低温余热发电装置组成的 CHP、HP、GB,参数见附表 A1;天然气价格为2.55 元/m’,分时电价见附表 A2;系统初始电、热负荷及风电预测出力见附图 A1；场景 3: 仅考虑需求响应。
4 程序运行结果
1）电负荷优化结果
2）热负荷优化
3）价格优化
4）电平衡
5）热平衡
5 程序
1）主程序
%% 碳交易机制下考虑需求响应的综合能源系统优化运行——魏震波%场景 3: 仅考虑需求响应;clc;clear;close all;% 程序初始化%% 读取数据shuju=[1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 241800 1760 1750 1750 1725 1800 1820 2300 2400 2700 2500 2050 2100 2080 2100 2120 2140 2200 2300 2400 2300 1900 2000 20201520 1500 1600 1650 1700 1900 1850 2100 2400 2350 2300 2400 2410 2430 2440 2450 2455 2460 2100 1950 2000 1890 1820 17500 0 0 0 0 204 429 599 770 1248 1200 1450 1350 1555 1476 1064 1071 795 325 0 0 0 0 01600 1800 1650 1640 1630 1630 1800 1750 1400 1300 1250 1500 1450 1600 1550 1520 1500 1550 1300 1350 1380 1400 1450 15000.78 0.78 0.78 0.78 0.78 0.78 0.78 0.78 0.78 0.78 0.78 0.78 0.78 0.78 0.78 0.78 0.78 0.78 0.78 0.78 0.78 0.78 0.78 0.780.35 0.35 0.35 0.35 0.35 0.35 0.35 0.68 0.68 1.09 1.09 1.09 0.68 0.68 0.68 0.68 0.68 0.68 0.68 1.09 1.09 1.09 0.68 0.680.37 0.37 0.37 0.37 0.37 0.37 0.37 0.37 0.37 0.37 0.37 0.37 0.37 0.37 0.37 0.37 0.37 0.37 0.37 0.37 0.37 0.37 0.37 0.370.25 0.25 0.25 0.25 0.25 0.25 0.25 0.36 0.36 0.36 0.58 0.58 0.58 0.58 0.36 0.36 0.36 0.25 0.25 0.58 0.58 0.58 0.36 0.360.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5]; %把一天划分为24小时load_e=shuju(2,:); %初始电负荷load_h=shuju(3,:); %初始热负荷P_PV=shuju(4,:); %光电预测P_WT=shuju(5,:); %风电预测pe_b=shuju(6,:); %需求响应前电价pe_a=shuju(7,:); %需求响应电价ph_b=shuju(8,:); %需求响应前热价ph_a=shuju(9,:); %需求响应热价
%% 需求侧定义变量Z=zeros(24,24); %需求弹性矩阵e_W1=0.5;e_W2=0.3;e_W3=0.15;e_W4=0.05;%约束：固定、可转移、可消减、可替代负荷占比50%，30%，15%，5% %这里进行4. 2. 2 需求响应灵敏度分析h_W1=0.5;h_W2=0.2;h_W3=0.2;h_W4=0.1;%约束：固定、可转移、可消减、可替代负荷占比50%，30%，15%，5% %这里进行4. 2. 2 需求响应灵敏度分析Psl_e=zeros(1,24);%转移电负荷量Pcl_e=zeros(1,24);%消减电负荷量Prl_e=zeros(1,24);%电负荷被替代量Psl_h=zeros(1,24);%转移热负荷量Pcl_h=zeros(1,24);%消减热负荷量Prl_h=zeros(1,24);%热负荷被替代量P2H=1.83; %电转热系数OP_load_e=zeros(1,24);%优化后的电负荷OP_load_h=zeros(1,24);%优化后的热负荷%% IES电网交互电价price_buy_grid=shuju(7,:); %向电网购电价price_sell_grid=shuju(10,:); %向电网售电价%% 供应测定义机组变量%CHPP_GT=sdpvar(1,24,'full');%燃气轮机输出功率e_GT=0.3;%燃气轮机供电效率h_GT=0.4;%燃气轮机供热效率P_WHB=sdpvar(1,24,'full');%余热锅炉输出功率r_WHB=0.80;%热回收效率P_ORC=sdpvar(1,24,'full');%ORC输出功率r_ORC=0.80;%ORC效率
P_GB=sdpvar(1,24,'full');%燃气锅炉输出功率h_GB=0.9;%燃气锅炉供热效率
P_HP=sdpvar(1,24,'full');%热泵输入功率COP_HP=4.4;%电制冷机冷系数
B_grid=sdpvar(1,24,'full');%购电电量 S_grid=sdpvar(1,24,'full');%售电电量 B_grid_sign=binvar(1,24,'full'); %购电标志
ES_char=sdpvar(1,24,'full');%储电系统充电ES_dischar=sdpvar(1,24,'full');%储电系统放电ES_char_sign=binvar(1,24,'full');%储电系统充电标志ES_max=400; ES_loss=0.01;ES_c_char=0.95;ES_c_discharge=0.9;%电储能最大容量；自损系数；充、放能效率
HS_char=sdpvar(1,24,'full');%储热系统充热HS_dischar=sdpvar(1,24,'full');%储热系统放热HS_char_sign=binvar(1,24,'full'); %储热系统充热标志HS_max=400; HS_loss=0.01;HS_c_char=0.95;HS_c_discharge=0.9;%热储能最大容量；自损系数；充、放能效率;原文0.8%% DR-需求侧响应优化Z_e=ElasticityMatrix(pe_a); %电价需求弹性矩阵Z_e_CL=diag(diag(Z_e)); %消减电负荷弹性矩阵,对角阵Z_e_SL=Z_e-Z_e_CL; %转移电负荷弹性矩阵
Z_h=ElasticityMatrix(ph_a); %热价需求弹性矩阵Z_h_CL=diag(diag(Z_h)); %消减热负荷弹性矩阵,对角阵Z_h_SL=Z_h-Z_h_CL; %转移热负荷弹性矩阵
%价格型需求响应[Psl_e,Pcl_e]=IBDR(Z_e_SL,Z_e_CL,load_e,pe_a,pe_b,e_W2,e_W3);%（转移电负荷弹性矩阵，削减电负荷弹性矩阵，初始电负荷，需求响应电价，需求响应前电价，可转移电负荷比例，可削减电负荷比例）[Psl_h,Pcl_h]=IBDR(Z_h_SL,Z_h_CL,load_h,ph_a,ph_b,h_W2,h_W3);%（转移热负荷弹性矩阵，削减热负荷弹性矩阵，初始热负荷，需求响应热价，需求响应前热价，可转移热负荷比例，可削减热负荷比例）%替代型需求响应[Prl_e,Prl_h]=RBDR(pe_a,ph_a,e_W4,h_W4,shuju);%（需求响应电价，需求响应热价，可替代电负荷占比，可替代热负荷占比）
OP_load_e=load_e Psl_e Pcl_e-Prl_e Prl_h/P2H;%优化后的电负荷（初始电负荷，转移电负荷，削减电负荷，电负荷被替代量，热负荷被替代量）OP_load_h=load_h Psl_h Pcl_h-Prl_h Prl_e*P2H;%优化后的热负荷（初始热负荷，转移热负荷，削减热负荷，热负荷被替代量，电负荷被替代量）%% IES供应侧储能约束 ES_start=80;HS_start=50; %电储能和热储能的初始能量for i=1:24 ES(1,i)=ES_start ES_char(1,i)*ES_c_char-ES_dischar(1,i)/ES_c_discharge; %储电初始容量约束 ES_start=ES(1,i);endfor i=1:23 ES(1,i 1)= ES(1,i)*(1-ES_loss) ES_char(1,i)*(1-ES_c_char)-ES_dischar(1,i)/ES_c_discharge; %储电容量约束endES_start=ES(1,24);
for i=1:24 EH(1,i)=HS_start HS_char(1,i)*(1-HS_c_char)-HS_dischar(1,i)/HS_c_discharge; %储热初始容量约束 HS_start=EH(1,i);endfor i=1:23 EH(1,i 1)= EH(1,i)*(1-HS_loss) HS_char(1,i)*HS_c_char-HS_dischar(1,i)/HS_c_discharge; %储热容量约束endHS_start=EH(1,24);
%% IES供应侧优化% 约束条件C=[];%%电储能设备运行约束 for i=1:24 %运行约束 C=[C,0<=es_char(1,i)<=250*es_char_sign(1,i)];<> C=[C,0<=es_dischar(1,i)<=250*(1-es_char_sign(1,i))];<> end for i=1:24 %余量约束 C=[C,0<=es(1,i)<=400];<> end %热储能设备运行约束 for i=1:24 %运行约束 C=[C,0<=hs_char(1,i)<=250*hs_char_sign(1,i)];<> C=[C,0<=hs_dischar(1,i)<=250*(1-hs_char_sign(1,i))];<> end for i=1:24 %余量约束 C=[C,0<=eh(1,i)<=400];<> end a=0.5; %这里进行4. 2. 3 GT 产热分配比例的影响%各个机组约束for i=1:24 C = [C,0<=p_gt(i)<=4000];%燃气轮机上下限约束<>% C = [C,0<=p_whb(i)<=1000];%余热锅炉上下限约束<> C = [C,0<=p_gb(i)<=1000];%燃气锅炉上下限约束> C = [C,0<=p_hp(i)<400];%热泵上下限约束<> C = [C,0<=p_orc(i)<=400];%orc上下限约束<> C = [C,P_GT(i)*h_GT*r_WHB*a<=p_whb(i)<=p_gt(i)*h_gt*r_whb*a];%余热回收分配公式，a为分配系数<> C = [C,P_GT(i)*h_GT*r_ORC*(1-a)<= p_orc=""> C = [C, 0<= b_grid=""> C = [C, 0<= s_grid="">end
%功率平衡约束for i=1:24 C = [C,B_grid(i)-S_grid(i) P_WT(i) P_PV(i) e_GT*P_GT(i) P_ORC(i)-P_HP(i)-ES_char(1,i) ES_dischar(1,i)==OP_load_e(i)]; %电平衡C = [C,P_WHB(i) P_GB(i) COP_HP*P_HP(i)-HS_char(1,i) HS_dischar(1,i)==OP_load_h(i)];%热平衡约束end
%% 目标函数%碳交易机制下考虑需求响应的综合能源系统以系统总收益最大为目标函数。（与原文不同）%收入，供应侧考虑用户侧需求响应时会给与经济补贴以鼓励社会责任Income=pe_a*OP_load_e' ph_a*OP_load_h' 6000;%包含系统运维成本、购售成本、碳交易成本，三部分构成成本% RIES运维成本GT=0.04;%燃气轮机单位运维成本WHB=0.025;%余热锅炉单位运维成本HP=0.025;%热泵单位运维成本PV=0.016;%光伏单位运维成本WT=0.018;%风机单位运维成本ES=0.018;%电储能单位运维成本HS=0.016;%热储能单位运维成本C_om=0;%运维成本for i=1:24C_om=C_om P_GT(i)*GT P_WHB(i)*WHB P_HP(i)*HP P_WT(i)*WT P_PV(i)*PV ES*(ES_char(1,i) ES_dischar(1,i)) HS*(HS_char(1,i) HS_dischar(1,i));end
H_gas=9.88;%天然气热值C_buy=0;%购能成本for i=1:24C_buy=C_buy B_grid(i)*price_buy_grid(i)-S_grid(i)*price_sell_grid(i) 2.55*(P_GT(i)/e_GT/H_gas P_GB(i)/h_GB/H_gas); end
% C_carbon_trade=0;%碳交易成本 PP=2.53;%电量的折算系数% for i=1:24% C_carbon_trade=C_carbon_trade 0.5*(0.6101-0.57)*(PP*e_GT*P_GT(i) h_GT*P_GT(i) P_GB(i)); %0.50yuan/t % end
%目标函数
f=C_om C_buy;op = sdpsettings('solver','cplex', 'verbose', 0);
optimize(C,f,op)CC=value(f) %总成本F=Income-CC%利润om=value(C_om);grid=value(C_buy);% car=value(C_carbon_trade)Q_carbon=0;%碳排放量 for i=1:24Q_carbon=Q_carbon (PP*e_GT*P_GT(i) h_GT*P_GT(i) P_GB(i)); endvalue(Q_carbon)%% huatux=1:24;figure(1)plot(x,OP_load_e,'-rs',x,load_e,'--bo');xlabel('时间/h');ylabel('电负荷/kW');title('需求响应前后电负荷曲线');legend('优化后电负荷','优化前电负荷');x=1:24;
figure(2)plot(x,OP_load_h,'-rs',x,load_h,'--bo');xlabel('时间/h');ylabel('热负荷/kW');title('需求响应前后热负荷曲线');legend('优化后热负荷','优化前热负荷');
figure(3)stairs(x,pe_b,'-b')hold onstairs(x,pe_a,'--b')hold onstairs(x,ph_b,'-r')hold onstairs(x,ph_a,'--r')title('价格曲线');legend('初始电价','分时电价','初始热价','分时热价');
figure(4)plot(x,P_PV,'-m')hold onplot(x,P_WT,'-c')title('风光预测');legend('光伏预测曲线','风机预测曲线');ee=value([B_grid;ES_dischar;(e_GT*P_GT P_ORC);P_WT;P_PV]);ee1=value([-ES_char;-P_HP;-S_grid]);
figure(5)bar(ee','stack');hold onbar(ee1','stack');plot(x,OP_load_e,'-gs');title('电负荷平衡');xlabel('时段');ylabel('功率/kW');legend('购电量','ES放电','CHP产电','风电出力','光伏出力','ES充电','HP耗电','卖电量','电负荷');hh=value([P_WHB;P_GB;HS_dischar;COP_HP*P_HP]);hh1=value([-HS_char]);
figure(6)bar(hh','stack');hold onbar(hh1','stack');plot(x,OP_load_h,'-rs');title('热负荷平衡');xlabel('时段');ylabel('功率/kW');legend('CHP产热','GB产热','HS放热','HP产热','HS充热','热负荷');
2）子程序
function [Z]=ElasticityMatrix(P)%%%总：用户需求弹性矩阵（分：根据文章分成SL型需求弹性矩阵，CL型需求弹性矩阵（对角阵））%时段 %谷 %平 %峰%低 -0.1 0.01 0.012%平 0.01 -0.1 0.016%峰 0.012 0.016 -0.1 %%数据来源《需求侧响应理论模型与应用研究_曾鸣》%构造需求弹性矩阵for i=1:24 if P(i)==min(P)%谷 for j=1:24 if P(j)==min(P)%低 Z(i,j)=-0.1; elseif P(j)==max(P)%峰 Z(i,j)=0.012; elseif min(P)
Z(i,j)=0.01; else%平 Z(i,j)=0.01; end end elseif P(i)==max(P)%峰 for j=1:24 if P(j)==min(P)%低 Z(i,j)=0.012; elseif P(j)==max(P)%峰 Z(i,j)=-0.1; elseif min(P)
Z(i,j)=0.016; else%平 Z(i,j)=0.016; end end else%平 for j=1:24 if P(j)==min(P)%低 Z(i,j)=0.01; elseif P(j)==max(P)%峰 Z(i,j)=0.016; elseif min(P)
Z(i,j)=-0.1; else%平 Z(i,j)=-0.1; end end endend
function [Psl,Pcl]=IBDR(Z_SL,Z_CL,load,p_a,p_b,W2,W3)Psl=zeros(1,24);%消减负荷量Pcl=zeros(1,24); %转移负荷量%价格型需求响应for i=1:24 sum1=0; for j=1:24 sum1=sum1 (Z_SL(i,j)*((p_a(j)-p_b(j))/p_b(j))); %可转移系数 end Psl(i)=W2*load(i)*sum1;end for i=1:24 sum2=0; for j=1:24 sum2=sum2 (Z_CL(i,j)*((p_a(j)-p_b(j))/p_b(j))); %可消减系数 end if sum2>0 sum2=0; %消减发生在价格上涨的时候 else Pcl(i)=W3*load(i)*sum2; endend
function [Prl_e,Prl_h]=RBDR(pe_a,ph_a,e_W4,h_W4,shuju)
%替代型需求响应%与文章不同，转换后成本较低方可发生替代P2H=1.83; %电转热系数%读取数据load_e=shuju(2,:); %初始电负荷load_h=shuju(3,:); %初始热负荷Prl_e=zeros(1,24);%电负荷被替代量Prl_h=zeros(1,24);%热负荷被替代量for i=1:24 if pe_a(i) Prl_e(i)=0; else Prl_e(i)=e_W4*load_e(i); endendfor i=1:24 if pe_a(i)/P2H Prl_h(i)=0; else Prl_h(i)=h_W4*load_h(i); endend
