计及碳捕集电厂低碳特性的含风电电力系统源-荷多时间尺度调度方法

参考文献：《计及碳捕集电厂低碳特性的含风电电力系统源-荷多时间尺度调度方法》
摘要: 逐步提升风电等可再生能源发电占比，并对火电机组进行低碳化改造，辅之以多类需求侧资源，是实现能源电力碳达峰、碳中和目标的重要手段。首先，挖掘源荷两侧低碳资源并分析其低碳特性,源侧在碳捕集电厂中装设烟气旁路系统与溶液存储器,形成碳捕集电广综合灵活运行方式进而与风电协调配合:荷侧调用不同响应速度的价格型、激励型需求响应资源克服多时间尺度下碳捕集电厂综合灵活运行方式的局限，通过源荷资源协调优化，从而提高系统的低碳性能。其次，构建源荷协调的日前-日内-实时三阶段低碳经济调度模型,优化系统的负荷及旋转备用分配计划，并改善失负荷与弃风问题。最后，在改进的IEEE-39 节点系统中进行算例分析，结果表明本文调度方法能够利用源荷可调节资源的调度优势,实现电力系统低碳经济调度的目标。
关键词:碳捕集电厂;低碳特性:风电:需求侧资源:多时间尺度调度
1 综合灵活运行方式低碳特性分析
碳捕集电厂的灵活运行方式大致分为 3 类: 分流式运行方式、储液式运行方式与综合灵活运行方式.在碳捕集电厂同时引入烟气旁路系统与溶液存储器，即形成综合灵活运行方式。相较于储液式运行方式，碳捕集电广综合灵活运行方式可在某一时段根据系统需求主动排放 CO,到大气中，提高调度与决策的灵活性: 相较于分流式运行方式，碳捕集电厂综合灵活运行方式既可在负荷高峰将碳捕集能耗转移至负荷低谷，缓解负荷需求与捕碳需求的矛盾，又可在负荷低谷提高碳捕集能耗、降低碳捕集电厂净出力以消纳部分弃风。

上述碳捕集能耗的转移是通过调整溶液存储器的储液量实现的: 当某时段富液存储器储液量增多，贫液存储器储液量减少时 (即该时段本该捕集的 CO,只进行了吸收过程而并未实现解析过程)，碳捕集能耗降低，净出力提升，反之亦然。换言之，溶液存储器可通过灵活调整储液量实现碳捕集能耗的能量时移,从而挖掘系统低碳潜力。

碳捕集电广综合灵活运行方式内部能量流动情况如图 1 所示。

由上述分析可知，综合灵活运行方式碳捕集电厂相较于分流式碳捕集电厂的净出力范围更大，如图2 所示。

通过源荷多时间尺度低碳经济调度方法，可利用风电及负荷预测精度随时间尺度缩短而提高的特性制定更加精确的调度计划，有效应对实时阶段失负荷、弃风问题的同时最大限度地挖掘系统的低碳性能，降低系统成本。

2 算例概述
本文基于改进的 IEEE-39 节点系统进行算例仿真，将火电厂 G1、G2 改造为碳捕集电厂，火电参数详见附录 B表 B1: 分别在节点 9、19、22引入300MW 风电场，系统结构如图 5 所示。附录 A 图 A1 为风电及负荷日前、日内、实时预测曲线:设定 A类、B 类、C 类IDR 的调用容量分别不超过总负荷的 5%、5%、3%; 价格需求弹性矩阵 E与与峰谷平三时段的划分详见附录 B表 B2 与表 B3;其余参数详见附录 B 表 B4。本文算例采用 CPLEX 进行优化求解。

3 程序运行结果
1）电平衡

2）热平衡

4 程序

%% 实时调度                   
clc
clear
close all
warning off % 关闭版本不兼容提示

%% 数据导入 
% load A_IDR
A_IDR=[1.74589264293358,1.74589264293358,1.74589264293358,1.74589264293358,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-77.0159429369334,-77.0159429369334,-77.0159429369334,-77.0159429369334,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,73.6278967890000,73.6278967890000,73.6278967890000,73.6278967890000,1.64215350499977,1.64215350499977,1.64215350499977,1.64215350499977;
    -17.8950766400000,-17.8950766400000,-17.8950766400000,-17.8950766400000,-20.5940638400000,-20.5940638400000,-20.5940638400000,-20.5940638400000,-20.9823152000000,-20.9823152000000,-20.9823152000000,-20.9823152000000,-21.1814758400000,-21.1814758400000,-21.1814758400000,-21.1814758400000,-20.7030096000000,-20.7030096000000,-20.7030096000000,-20.7030096000000,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-19.0348424000000,-19.0348424000000,-19.0348424000000,-19.0348424000000,-20.3812280000000,-20.3812280000000,-20.3812280000000,-20.3812280000000;
    -20.0858924392000,-20.0858924392000,-20.0858924392000,-20.0858924392000,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,12.9925845112000,12.9925845112000,12.9925845112000,12.9925845112000,12.7131992288000,12.7131992288000,12.7131992288000,12.7131992288000,16.6082548632000,16.6082548632000,16.6082548632000,16.6082548632000,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,15.1456744432000,15.1456744432000,15.1456744432000,15.1456744432000,12.7460334904000,12.7460334904000,12.7460334904000,12.7460334904000,0,0,0,0,0,0,0,0,-9.89475745440003,-9.89475745440003,-9.89475745440003,-9.89475745440003,-19.4862098872000,-19.4862098872000,-19.4862098872000,-19.4862098872000,-20.7388867560000,-20.7388867560000,-20.7388867560000,-20.7388867560000;
    -5.02147310980000,-5.02147310980000,-5.02147310980000,-5.02147310980000,-5.01781581000000,-5.01781581000000,-5.01781581000000,-5.01781581000000,-5.19103884460000,-5.19103884460000,-5.19103884460000,-5.19103884460000,-5.17408227000000,-5.17408227000000,-5.17408227000000,-5.17408227000000,-4.82497634640000,-4.82497634640000,-4.82497634640000,-4.82497634640000,-4.50645875120000,-4.50645875120000,-4.50645875120000,-4.50645875120000,-4.93236797900000,-4.93236797900000,-4.93236797900000,-4.93236797900000,-4.57162519100000,-4.57162519100000,-4.57162519100000,-4.57162519100000,-4.42042063120000,-4.42042063120000,-4.42042063120000,-4.42042063120000,-3.71664747160000,-3.71664747160000,-3.71664747160000,-3.71664747160000,-3.24814612780000,-3.24814612780000,-3.24814612780000,-3.24814612780000,-3.17829980720000,-3.17829980720000,-3.17829980720000,-3.17829980720000,-4.15206371580000,-4.15206371580000,-4.15206371580000,-4.15206371580000,-4.09038655480000,-4.09038655480000,-4.09038655480000,-4.09038655480000,-3.73911266120000,-3.73911266120000,-3.73911266120000,-3.73911266120000,-3.83591721100000,-3.83591721100000,-3.83591721100000,-3.83591721100000,-4.03442780720000,-4.03442780720000,-4.03442780720000,-4.03442780720000,-3.78641861080000,-3.78641861080000,-3.78641861080000,-3.78641861080000,-3.18650837260000,-3.18650837260000,-3.18650837260000,-3.18650837260000,-4.83770668920000,-4.83770668920000,-4.83770668920000,-4.83770668920000,-4.51144597960000,-4.51144597960000,-4.51144597960000,-4.51144597960000,-4.79358702340000,-4.79358702340000,-4.79358702340000,-4.79358702340000,-4.87155247180000,-4.87155247180000,-4.87155247180000,-4.87155247180000,-5.18472168900000,-5.18472168900000,-5.18472168900000,-5.18472168900000];

A_IDR_P_tran=A_IDR(1,:); % 日前调度需求响应--电力转移
A_IDR_P_cut=A_IDR(2,:); % 日前调度需求响应--电力削减
A_IDR_H_tran=A_IDR(3,:); % 日前调度需求响应--热力转移
A_IDR_H_cut=A_IDR(4,:); % 日前调度需求响应--热力削减
% 调度周期
T=96;
% 风电预测功率
P_prew=(1/1.79)*1.4*[552.300186	589.654186	562.712186	516.618186	526.004186	526.726186	457.528186	493.628186	595.848186	551.768186	498.036186	531.324186	589.160186	561.002186	547.816186	583.650186	567.386186	648.326186	604.170186	528.398186	583.460186	570.730186	495.642186	500.848186	384.492186	331.520186	311.038186	316.586186	291.886186	313.356186	338.9000675	312.3131263	306.7788154	366.4881905	387.9301691	324.7722787	304.6930198	390.3218814	399.1100335	380.6716005	187.7772243	258.3049259	283.1119881	263.3386459	153.4867448	168.9772534	124.3690384	82.29158868	135.0483118	115.2749696	118.9181592	133.3518647	108.7491384	179.7774309	192.0697196	250.2773219	292.7997413	230.3927373	183.6431054	205.3353796	163.3413617	274.8340886	294.523999	246.0501095	185.8957647	294.2737035	297.221628	213.2892134	353.2321922	376.4886168	364.4886168	333.8486168	436.6086168	526.4086168	480.1949288	426.2349288	453.44	573.52	582.72	605.08	543.48	592.96	637.277248	627.442048	477.072768	658.608704	638.763456	603.706432	598.460992	696.1136	621.14752	675.78752	693.141184	674.432448	689.862784	679.153344];
% 电负荷
EleLoad=0.74*[1256.6848	1209.1268	1294.3448	1338.722	1306.5164	1391.4908	1376.0684	1303.0388	1394.3636	1417.724	1406.3084	1408.8032	1433.6	1431.1808	1220.94	1278.27	1308.636	1398.852	1497.762	1431.99	1268.568	1322.37	1312.164	1335.726	1392.93	1410.948	1481.256	1548.162	1610.658	1492.848	1558.872	1627.794	1563.786	1554.336	1287.09	1209.6	1393.812	1390.158	1543.374	1588.608	1685.124	1770.804	1848.294	1896.174	1881.81	1786.932	1622.754	1495.494	1464.624	1512.252	1547.154	1439.802	1332.198	1305.108	1316.196	1311.786	1339.044	1210.398	1190.868	1231.944	1282.806	1228.514	1148.252	1173.956	1254.792	1440.852	1476.72	1536.318	1730.106	1770.426	1857.114	1923.012	1949.85	1969.128	1964.214	1950.606	1808.1	1685.124	1498.14	1461.222	1344.098	1488.746	1422.064	1439.074	1426.978	1445.374	1425.466	1384.138	1354.052	1286.138	1488.746	1422.456	1414.336	1377.11	1349.768	1415.918];
% 热负荷
HeatLoad=(1/3.93)*0.7*[1478.665213	1409.599237	1375.159582	1361.999713	1388.879444	1408.572581	1418.559148	1418.185818	1452.532141	1457.198761	1468.585314	1462.238711	1475.771909	1452.438809	1425.185748	1401.665983	1399.799335	1354.439789	1292.187078	1231.521018	1178.041553	1265.02735	1288.080452	1285.747142	1416.785832	1384.586154	1378.239551	1417.532491	1357.98642	1283.3205	987.5912801	1056.50194	1190.31016	1240.87522	1158.89332	1184.80648	1110.16282	1043.31604	1011.66988	941.1539801	927.2801201	911.8010201	898.8444401	865.4783801	896.8952201	892.1941601	902.9722001	939.0901001	1120.59688	1165.5436	1119.1063	1113.02932	1156.14148	1148.22994	1112.68534	1068.88522	1082.64442	1049.62234	1059.94174	1024.5118	1080.23656	1076.79676	1141.57966	1200.62956	1191.45676	1132.52152	1004.4463	982.3169201	991.3750601	1062.901796	1148.020945	1145.50097	1092.414742	894.4984216	921.125442	1002.218144	1352.56823	1358.013378	1345.298549	1357.089735	1383.180466	1266.427336	1364.482582	1355.097642	1362.167174	1345.628357	1352.681646	1356.569756	1347.005506	1367.514373	1370.908669	1366.838104	1482.118512	1455.425446	1475.958574	1486.038473];
% 分时电价
gama_elec=100*[0.29	0.29	0.29	0.29	0.29	0.29	0.29	0.29	0.29	0.29	0.29	0.29	0.29	0.29	0.29	0.29	0.29	0.29	0.29	0.29	0.69	0.69	0.69	0.69	0.69	0.69	0.69	0.69	0.69	0.69	0.69	0.69	0.84	0.84	0.84	0.84	0.84	0.84	0.84	0.84	0.84	0.84	0.84	0.84	0.69	0.69	0.69	0.69	0.69	0.69	0.69	0.69	0.69	0.69	0.69	0.69	0.69	0.69	0.69	0.69	0.69	0.69	0.69	0.69	0.69	0.69	0.69	0.69	0.69	0.69	0.69	0.69	0.64	0.64	0.64	0.64	0.64	0.64	0.64	0.64	0.55	0.55	0.55	0.55	0.55	0.55	0.55	0.55	0.29	0.29	0.29	0.29	0.29	0.29	0.29	0.29];

%% 相关参数
a=[0.00048,0.00031,0.0002]; % 成本系数a
b=[16.2,17.3,16.6]; % 成本系数b
c=[1000,970,700]; % 成本系数c
eg=[0.91,0.95,0.98]; % 碳排放系数
GT_P_max=[150,80];
E_bata=0.81; % 碳捕集效率
P_yita=[1.05,1.05,1.05]; % 再生塔和压缩机最大工作状态系数
P_lamda=[0.269,0.269,0.269]; % 单位捕碳量能耗
P_G_max=[400,455,200]; % 机组最大出力
P_G_min=[200,120,100]; % 机组最小出力
P_D1=0.25*0.5*[20,20,15]; % 碳捕集固定能耗
P_D=ones(3,96);
for i=1:3
    P_D(i,:)=P_D1(i).*ones(1,96);  
end
R_u_intra=[50,50,25]; % 机组上爬坡
M_MEA=61.08; % M_MEA的摩尔质量
M_co2=44; % 二氧化碳的摩尔质量
sita=0.4; % 再生塔解析量
CR=30; % 醇胺溶液浓度(%)
rou_R=1.01; % 醇胺溶液密
V_CR=[2000,2000,2000]; % 溶液储液装置容量
K_R=1.17; % 乙醇胺溶剂成本系数
fai=1.5; % 溶剂运行损耗系数   
h_min=[200,120,100];
h_max=[401,456,201];
h_gn=5;  
hl1=(h_max-h_min)./h_gn;
hl2=zeros(3,h_gn+1);
k_min=[-400,-455,-200];
k_max=[-199,-119,-99];
k_gn=5; 
kl1=(k_max-k_min)./k_gn;
kl2=zeros(3,k_gn+1);
for i=1:3
hl2(i,:)=h_min(i):hl1(i):h_max(i);
end
for i=1:3
kl2(i,:)=k_min(i):kl1(i):k_max(i);
end

%% 决策变量
% 电力源出力
GT_P = sdpvar(2,96,'full'); % 燃气轮机电出力
P_w = sdpvar(1,96,'full'); % 风电机组出力
P_G = sdpvar(3,96,'full'); % 火电机组出力
EB=sdpvar(2,96,'full'); % 电锅炉出力
P_s_intra=sdpvar(1,96,'full'); % 失电负荷功率
% 热力源出力
GT_H = sdpvar(2,96,'full'); % 燃气轮机热出力
EB_H=sdpvar(2,96,'full'); % 电锅炉热出力
H_s_intra=sdpvar(1,96,'full'); % 失热负荷功率
% 天然气
P_gas=sdpvar(2,96,'full'); % 天然气需求
% 碳捕集相关
E_G=sdpvar(3,96,'full'); % 碳捕集机组产生的总碳排放
E_total_co2=sdpvar(3,96,'full'); % 机组捕获的总碳排放
E_CG=sdpvar(3,96,'full'); % 储液装置提供的待捕集二氧化碳量
P_B=sdpvar(3,96,'full'); % 机组运行能耗
P_J=sdpvar(3,96,'full'); % 机组净出力
V_CA=sdpvar(3,96,'full'); % 机组净出力
V_FY=sdpvar(3,96,'full'); % 富液体积
V_PY=sdpvar(3,96,'full'); % 贫液体积
% 电力需求响应
P_tran_intra=sdpvar(1,96,'full'); % 系统可转移电负荷
P_cut_intra=sdpvar(1,96,'full'); % 系统可削减电负荷
P_DE_intra=sdpvar(1,96,'full'); % 系统经过过需求响应后的电负荷
H_tran_intra=sdpvar(1,96,'full'); % 系统可转移热负荷
H_cut_intra=sdpvar(1,96,'full'); % 系统可削减热负荷
H_DE_intra=sdpvar(1,96,'full'); % 系统经过过需求响应后的热负荷
gn=5;  
P_G_line= sdpvar(3,96,'full'); % 火电机组出力
w1=sdpvar(gn+1,T,'full');
w2=sdpvar(gn+1,T,'full');
w3=sdpvar(gn+1,T,'full');
w4=sdpvar(gn+1,T,'full');
w5=sdpvar(gn+1,T,'full');
z1=binvar(gn, T,'full');
z2=binvar(gn, T,'full');
z3=binvar(gn, T,'full');
z4=binvar(gn, T,'full');
z5=binvar(gn, T,'full');
y1= sdpvar(3,96,'full'); 
h_w1=sdpvar(h_gn+1,T,'full');
h_w2=sdpvar(h_gn+1,T,'full');
h_w3=sdpvar(h_gn+1,T,'full');
h_w4=sdpvar(h_gn+1,T,'full');
h_w5=sdpvar(h_gn+1,T,'full');
h_z1=binvar(h_gn,T,'full');
h_z2=binvar(h_gn,T,'full');
h_z3=binvar(h_gn,T,'full');
h_z4=binvar(h_gn,T,'full');
h_z5=binvar(h_gn,T,'full');
y2= sdpvar(3,96,'full'); 
k_w1=sdpvar(k_gn+1,T,'full');
k_w2=sdpvar(k_gn+1,T,'full');
k_w3=sdpvar(k_gn+1,T,'full');
k_w4=sdpvar(k_gn+1,T,'full');
k_w5=sdpvar(k_gn+1,T,'full');
k_z1=binvar(k_gn, T,'full');
k_z2=binvar(k_gn, T,'full');
k_z3=binvar(k_gn, T,'full');
k_z4=binvar(k_gn, T,'full');
k_z5=binvar(k_gn, T,'full');

%% 约束条件
C = [];  %约束条件初始
for t=1:96
    for i=1:3
    C = [C,         
        0<=E_CG(i,t),
        0<=P_B(i,t),
        E_G(i,t)==eg(i)*P_G(i,t), % 碳捕集机组产生的总碳排放   
        E_total_co2(i,t)==E_CG(i,t)+0.25*E_bata*eg(i)*(y1(i,t)-y2(i,t)), % 机组捕获的二氧化碳总量
        0<=E_total_co2(i,t)<=P_yita*E_bata*eg(i)*P_G_max(i),   
        P_B(i,t)==P_lamda(i)*E_total_co2(i,t), % 机组运行能耗
        P_G(i,t)==P_J(i,t)+P_D(i,t)+P_B(i,t), % 机组输出总功率
        P_G_min(i)-P_lamda(i)*P_yita*E_bata*eg(i)*P_G_max(i)-P_D(i)<=P_J(i,t)<=P_G_max(i)-P_D(i), % 碳捕集电厂净出力范围
        0<= P_w(t)<= P_prew(t), % 风电出力区间约束 
        sum(EB(:,t))+P_w(t)<=P_prew(t);
        P_G_min(i)<= P_G(i,t)<=P_G_max(i), % 火电机组出力约束
        0<=H_s_intra(t)<=20;
        0<=P_s_intra(t)<=20;
        ];
    end
end
C=[C,min(sum(R_u_intra),sum(P_G_max)-sum(P_G))>=0.01*max(P_DE_intra),]; % 旋转备用约束
C = [C, y1(1,:)==hl2(1,:).^2*h_w1];
C = [C, y1(2,:)==hl2(2,:).^2*h_w2];
C = [C, y1(3,:)==hl2(3,:).^2*h_w3];
C = [C, h_w1(1,:)<=h_z1(1,:)];
for i=2:h_gn
    C = [C, h_w1(i,:)<=h_z1(i-1,:)+h_z1(i,:)];
end
C = [C, h_w1(h_gn+1,:)<=h_z1(h_gn,:)];
C = [C, sum(h_z1)==ones(1,T)];
C = [C, h_w2(1,:)<=h_z2(1,:)];
for i=2:h_gn
    C = [C, h_w2(i,:)<=h_z2(i-1,:)+h_z2(i,:)];
end
C = [C, h_w2(h_gn+1,:)<=h_z2(h_gn,:)];
C = [C, sum(h_z2)==ones(1,T)];
C = [C, h_w3(1,:)<=h_z3(1,:)];
for i=2:h_gn
    C = [C, h_w3(i,:)<=h_z3(i-1,:)+h_z3(i,:)];
end
C = [C, h_w3(h_gn+1,:)<=h_z3(h_gn,:)];
C = [C, sum(h_z3)==ones(1,T)];
C = [C, h_w4(1,:)<=h_z4(1,:)];
for i=2:h_gn
    C = [C, h_w4(i,:)<=h_z4(i-1,:)+h_z4(i,:)];
end
C = [C, h_w4(h_gn+1,:)<=h_z4(h_gn,:)];
C = [C, sum(h_z4)==ones(1,T)];
C = [C, h_w5(1,:)<=h_z5(1,:)];
for i=2:h_gn
    C = [C, h_w5(i,:)<=h_z5(i-1,:)+h_z5(i,:)];
end
C = [C, h_w5(h_gn+1,:)<=h_z5(h_gn,:)];

C = [C, sum(h_z3)==ones(1,T)];
C = [C, y2(1,:)==kl2(1,:).^2*k_w1];
C = [C, y2(2,:)==kl2(2,:).^2*k_w2];
C = [C, y2(3,:)==kl2(3,:).^2*k_w3];
C = [C, k_w1(1,:)<=k_z1(1,:)];
for i=2:k_gn
    C = [C, k_w1(i,:)<=k_z1(i-1,:)+k_z1(i,:)];
end
C = [C, k_w1(k_gn+1,:)<=k_z1(k_gn,:)];
C = [C, sum(k_z1)==ones(1,T)];
C = [C, k_w2(1,:)<=k_z2(1,:)];
for i=2:k_gn
    C = [C, k_w2(i,:)<=k_z2(i-1,:)+k_z2(i,:)];
end
C = [C, k_w2(k_gn+1,:)<=k_z2(k_gn,:)];
C = [C, sum(k_z2)==ones(1,T)];
C = [C, k_w3(1,:)<=k_z3(1,:)];
for i=2:k_gn
    C = [C, k_w3(i,:)<=k_z3(i-1,:)+k_z3(i,:)];
end
C = [C, k_w3(k_gn+1,:)<=k_z3(k_gn,:)];
C = [C, sum(k_z3)==ones(1,T)];
C = [C, k_w4(1,:)<=k_z4(1,:)];
for i=2:k_gn
    C = [C, k_w4(i,:)<=k_z4(i-1,:)+k_z4(i,:)];
end
C = [C, k_w4(k_gn+1,:)<=k_z4(k_gn,:)];
C = [C, sum(k_z4)==ones(1,T)];
C = [C, k_w5(1,:)<=k_z5(1,:)];
for i=2:k_gn
    C = [C, k_w5(i,:)<=k_z5(i-1,:)+k_z5(i,:)];
end
C = [C, k_w5(k_gn+1,:)<=k_z5(k_gn,:)];
C = [C, sum(k_z3)==ones(1,T),y1>=y2,];

%% 其他约束
% 爬坡约束
for  t = 1: 95
    for i=1:3
         C = [C,   
         P_G(i,t+1)-P_G(i,t)<=R_u_intra(i), % 火电机组上爬坡出力约束
        ];
    end
    for j=1:2
        C = [C,  
         GT_P(j,t+1)-GT_P(j,t)<=0.25*50, % 燃气轮机上爬坡出力约束
         EB(j,t+1)-EB(j,t)<=0.25*30, % 电锅炉上爬坡出力约束
         ];
    end
end
for t=1:96
    for i=1:2
        C=[C,
            0<=P_gas(i,t)<=130, % 燃气供给功率
            0<=GT_P(i,t)<=GT_P_max(i), % 燃气供给功率
            GT_P(i,t)==0.98*P_gas(i,t), % 燃气轮机电功率    
            GT_H(i,t)==0.6*(GT_P(i,t)/0.4), % 燃气轮机供热功率
            0<=EB(i,t)<=P_prew(t)-P_w(t), % 电锅炉出力
            EB_H(i,t)==0.98*EB(i,t), % 电锅炉热出力
          ];
    end
end
% 综合灵活运行方式碳捕集电厂约束
for t=1:96
    for i=1:3
        C = [C,
        V_CA(i,t)==0.25*((M_MEA*E_CG(i,t))/(M_co2*sita*CR*rou_R)), % 释放二氧化碳所需的溶液体积
        0<=V_FY(i,t)<=V_CR(i), % 富液体积
        0<=V_PY(i,t)<=V_CR(i), % 贫液体积
        0<=(0.25*eg(i)*(y1(i,t)-y2(i,t)))<=P_G(i,t),
        ];
    end
end 
V_PY0=[1000,1000,1000];
V_FY0=[1000,1000,1000];
for i=1:3
    C=[C,
    V_PY0(i)==V_PY(i,96),% 贫液初始平衡
    V_FY0(i)==V_FY(i,96),% 富液初始平衡
    ];
end
for t=1:95
    for i=1:3
        C = [C,
        V_FY(i,t+1)==V_FY(i,t)-V_CA(i,t+1), % 富液变化
        V_PY(i,t+1)==V_PY(i,t)+V_CA(i,t+1), % 贫液变化
        ];
    end
end
gl1=(P_G_max-P_G_min)./gn;
gl2=zeros(3,gn+1);
for i=1:3
gl2(i,:)=P_G_min(i):gl1(i):P_G_max(i);
end
C = [C, P_G_line(1,:)==gl2(1,:).^2*w1];
C = [C, P_G_line(2,:)==gl2(2,:).^2*w2];
C = [C, P_G_line(3,:)==gl2(3,:).^2*w3];
C = [C, w1(1,:)<=z1(1,:)];
for i=2:gn
    C = [C, w1(i,:)<=z1(i-1,:)+z1(i,:)];
end
C = [C, w1(gn+1,:)<=z1(gn,:)];
C = [C, sum(z1)==ones(1,T)];
C = [C, w2(1,:)<=z2(1,:)];
for i=2:gn
    C = [C, w2(i,:)<=z2(i-1,:)+z2(i,:)];
end
C = [C, w2(gn+1,:)<=z2(gn,:)];
C = [C, sum(z2)==ones(1,T)];
C = [C, w3(1,:)<=z3(1,:)];
for i=2:gn
    C = [C, w3(i,:)<=z3(i-1,:)+z3(i,:)];
end
C = [C, w3(gn+1,:)<=z3(gn,:)];
C = [C, sum(z3)==ones(1,T)];
C = [C, w4(1,:)<=z4(1,:)];
for i=2:gn
    C = [C, w4(i,:)<=z4(i-1,:)+z4(i,:)];
end
C = [C, w4(gn+1,:)<=z4(gn,:)];
C = [C, sum(z4)==ones(1,T)];
C = [C, w5(1,:)<=z5(1,:)];
for i=2:gn
    C = [C, w5(i,:)<=z5(i-1,:)+z5(i,:)];
end
C = [C, w5(gn+1,:)<=z5(gn,:)];
C = [C, sum(z3)==ones(1,T)];
C = [C, P_G(1,:)==gl2(1,:)*w1];
C = [C, P_G(2,:)==gl2(2,:)*w2];
C = [C, P_G(3,:)==gl2(3,:)*w3];
% 电力需求响应
for t=1:96
    if 44<=t<=48 
        C = [C,
            -0.02*(EleLoad(t)+A_IDR_P_tran(t)+A_IDR_P_cut(t))<=P_cut_intra(t)<=0  % 电负荷削减
            ];
    elseif 72<=t<=80
         C = [C,
            -0.02*(EleLoad(t)+A_IDR_P_tran(t)+A_IDR_P_cut(t))<=P_cut_intra(t)<=0  % 电负荷削减
            ];
    else
        C = [C,
            P_cut_intra(t)==0  % 电负荷不削减
            ];
    end
end
% 热力需求响应
for t=1:96
    if 0<=t<=12 
        C = [C,
            -0.02*(EleLoad(t)+A_IDR_H_tran(t)+A_IDR_H_cut(t))<=H_cut_intra(t)<=0  % 热负荷削减
            ];
    elseif 88<=t<=96
         C = [C,
            -0.02*(HeatLoad(t)+A_IDR_H_tran(t)+A_IDR_H_cut(t))<=H_cut_intra(t)<=0  % 热负荷削减
            ];
    else
        C = [C,
            H_cut_intra(t)==0  % 热负荷不削减
            ];
    end
end
C = [C,
     -0.08*(EleLoad+A_IDR_P_tran+A_IDR_P_cut)<=P_tran_intra<=0.08*(EleLoad+A_IDR_P_tran+A_IDR_P_cut), % 电负荷转移
     sum(P_tran_intra)==0, % 电负荷转移守恒
     P_DE_intra==(EleLoad+A_IDR_P_tran+A_IDR_P_cut)+P_tran_intra+P_cut_intra, % 需求响应后的电负荷
     -0.08*(HeatLoad+A_IDR_H_tran+A_IDR_H_cut)<=H_tran_intra<=0.08*(HeatLoad+A_IDR_H_tran+A_IDR_H_cut), % 热负荷转移
     sum(H_tran_intra)==0, % 热负荷转移守恒
     H_DE_intra==(HeatLoad+A_IDR_H_tran+A_IDR_H_cut)+H_tran_intra+H_cut_intra, % 需求响应后的热负荷
     ];
     C = [C,
         P_w+sum(GT_P)+sum(P_J)+P_s_intra-P_DE_intra==0, % 电力平衡
         sum(EB_H)+sum(GT_H)+H_s_intra-H_DE_intra==0, % 热力平衡
     ];
 
 %% 目标函数
 F=0;
 for t=1:96
     F1(t)=0.35*gama_elec(t)*abs(P_tran_intra(t))+0.7*gama_elec(t)*abs(P_cut_intra(t)); % 电能需求响应成本
     F2(t)=0.35*abs(H_tran_intra(t))+0.7*abs(H_cut_intra(t)); % 热能需求响应成本
     for i=1:3
         F3(i,t)=a(i)*P_G_line(i,t)+b(i)*P_G(i,t)+c(i); % 火电机组燃料费用
         F4(i,t)=0.25*E_bata*eg(i)*(y1(i,t)-y2(i,t))+E_CG(i,t); % 碳捕集量
         F5(i,t)=E_G(i,t)-0.25*E_bata*eg(i)*(y1(i,t)-y2(i,t))-E_CG(i,t)-0.7*P_G(i,t); % 碳排放
     end
     F6(t)=50*sum(P_gas(:,t)); % 燃气费用
     F7(t)=gama_elec(t)*sum(EB(:,t)); % 电锅炉运行成本
     F8(t)=0.3*gama_elec(t)*(P_prew(t)-P_w(t)-sum(EB(:,t))); % 弃风惩罚
     F9(t)=(H_s_intra(t)+P_s_intra(t)); % 失负荷惩罚
 end
 F=F+sum(F1)+sum(F2)+sum(sum(F3))+1.17*1.5*sum(sum(F4))+12*sum(sum(F5))+sum(F6)+sum(F7)+sum(F8)+142.857*sum(F9);
 F=0.25*F;
 
 %% 求解器配置
ops=sdpsettings('solver','cplex','verbose',2,'usex0',0);
ops.cplex.mip.tolerances.mipgap=1e-6;
result=optimize(C,F,ops);
F=value(F);
display(['最优规划结果 : ', num2str(F)]);

%% 数据获取
GT_P=value(GT_P); % 燃气轮机电出力
P_w=value(P_w); % 风电机组出力
P_G=value(P_G); % 火电机组出力
EB=value(EB); % 电锅炉出力
H_s_intra=value(H_s_intra);
P_s_intra=value(P_s_intra);
% 热力源出力
GT_H=value(GT_H) ; % 燃气轮机热出力
EB_H=value(EB_H); % 电锅炉热出力
% 天然气
P_gas=value(P_gas); % 天然气需求
% 碳捕集相关
E_G=value(E_G); % 碳捕集机组产生的总碳排放
E_total_co2=value(E_total_co2); % 机组捕获的总碳排放
E_CG=value(E_CG); % 储液装置提供的待捕集二氧化碳量
P_B=value(P_B); % 机组运行能耗
P_J=value(P_J); % 机组净出力
V_CA=value(V_CA); % 机组净出力
y1=value(y1);
y2=value(y2);
% 电力需求响应
P_tran_intra=value(P_tran_intra); % 系统可转移电负荷
P_cut_intra=value(P_cut_intra); % 系统可削减电负荷
P_DE_intra=value(P_DE_intra); % 系统经过过需求响应后的电负荷
% 热力需求响应
H_tran_intra=value(H_tran_intra); % 系统可转移热负荷
H_cut_intra=value(H_cut_intra); % 系统可削减热负荷
H_DE_intra=value(H_DE_intra); % 系统经过过需求响应后的热负荷

%% 数据分析与画图
% 电力系统
figure
for t=1:96
    Plot_EleNet(1,t)=(P_J(1,t));
    Plot_EleNet(2,t)=(P_J(2,t));
    Plot_EleNet(3,t)=(P_J(3,t));
    Plot_EleNet(4,t)=P_w(t)+sum(EB(:,t));
    Plot_EleNet(5,t)=GT_P(1,t);
    Plot_EleNet(6,t)=GT_P(2,t);
    Plot_EleNet(7,t)= P_s_intra(t);
    Plot_EleNet(8,t)=sum(EB(:,t));
end
Plot_EleNet=Plot_EleNet';
bar(Plot_EleNet,'stacked');
hold on
plot(EleLoad+A_IDR_P_tran+A_IDR_P_cut,'k-','LineWidth',1.5);
plot(P_DE_intra,'r-o','LineWidth',1.5);
xlabel('时间/h');
ylabel('功率/kW');
title('电功率平衡');
hold on
legend('火电机组1','火电机组2','火电机组3','风电','燃气轮机1','燃气轮机2','失负荷','电锅炉','电负荷','P-DE');


figure
% 热力系统
for t=1:96
    Plot_HeatNet(1,t)=sum(GT_H(:,t));
    Plot_HeatNet(2,t)=sum(EB_H(:,t));
    Plot_HeatNet(3,t)=H_s_intra(t);
end
Plot_HeatNet=Plot_HeatNet';
bar(Plot_HeatNet,'stacked');
hold on
plot(HeatLoad+A_IDR_H_tran+A_IDR_H_cut,'k-','LineWidth',1.5);
plot(H_DE_intra,'r-o','LineWidth',1.5);
xlabel('时间/h');
ylabel('功率/kW');
title('热功率平衡');
hold on
legend('燃气轮机','电锅炉','失负荷','热负荷','H-DE');