碳交易机制下考虑需求响应的综合能源系统优化运行

碳交易机制下考虑需求响应的综合能源系统优化运行
参考文献：《碳交易机制下考虑需求响应的综合能源系统优化运行》
摘 要:综合能源系统是实现“双碳”目标的有效途径,为进一步挖掘其需求侧可调节潜力对碳减排的作用,提出了一种碳交易机制下考虑需求响应的综合能源系统优化运行模型。首先,根据负荷响应特性将需求响应分为价格型和替代型 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  24
1800  1760  1750  1750  1725  1800  1820  2300  2400  2700  2500  2050  2100  2080  2100  2120  2140  2200  2300  2400  2300  1900  2000  2020
1520  1500  1600  1650  1700  1900  1850  2100  2400  2350  2300  2400  2410  2430  2440  2450  2455  2460  2100  1950  2000  1890  1820  1750
0  0  0  0  0  204  429  599  770  1248  1200  1450  1350  1555  1476  1064  1071  795  325  0  0  0  0  0
1600  1800  1650  1640  1630  1630  1800  1750  1400  1300  1250  1500  1450  1600  1550  1520  1500  1550  1300  1350  1380  1400  1450  1500
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.78  0.78
0.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.68
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.37  0.37
0.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.36
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  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,:); %向电网售电价
%% 供应测定义机组变量
%CHP
P_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);
end
for 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; %储电容量约束
end
ES_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);
end
for 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; %储热容量约束
end
HS_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(i)<=P_GT(i)*h_GT*r_ORC*(1-a)];
    
    C = [C, 0<= B_grid(i)<= B_grid_sign*1500];
   C = [C, 0<= S_grid(i)<=(1-B_grid_sign)*1500]; %外部电网联络线约束
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:24
C_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:24
C_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:24
Q_carbon=Q_carbon+(PP*e_GT*P_GT(i)+h_GT*P_GT(i)+P_GB(i));                              
end
value(Q_carbon)
%% huatu
x=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 on
stairs(x,pe_a,'--b')
hold on
stairs(x,ph_b,'-r')
hold on
stairs(x,ph_a,'--r')
title('价格曲线');
legend('初始电价','分时电价','初始热价','分时热价');

figure(4)
plot(x,P_PV,'-m')
hold on
plot(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 on
bar(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 on
bar(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)<P(j)<max(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)<P(j)< max(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)<P(j)< max(P)%平
                Z(i,j)=-0.1;
           else%平 
                Z(i,j)=-0.1;
           end
         end  
    end
end  

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;
    end
end

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)<P2H*ph_a(i)  %转换后价格成本较低则替换
      Prl_e(i)=0;
    else
      Prl_e(i)=e_W4*load_e(i);
    end
end
for i=1:24
    if pe_a(i)/P2H<ph_a(i) %转换后价格成本较低则替换
      Prl_h(i)=0;
    else
      Prl_h(i)=h_W4*load_h(i);
    end
end