clc
clear
%% 常量导入
% 电负荷
P_e_load= [1125.91468984867,1116.10893523230,1144.36212103951,1119.86779733002,1257.01203834709,1303.31364247203,1305.98753683134,1409.63778403649,1401.08325010268,1479.34743427615,1440.57455508449,1443.70276087454,1400.25137665787,1344.76042648101,1435.83768792681,1327.45950721548,1340.29160985833,1419.81291587739,1345.16739885315,1409.97900818528,1334.62530572066,1316.02226270621,1206.31538900082,1189.90502757710];
% 热负荷
P_h_load= [1251.03276366450,1272.94389557152,1307.94077358185,1305.48028523476,1360.36566863060,1259.66357074242,1243.00413398017,1336.52178221359,1120.05044137236,963.363623213736,956.352712044075,904.489802927447,890.210713727788,907.806731250927,826.780220676256,881.655248987816,1051.81754899827,1110.15277312622,1118.92034209325,1106.97368862136,1102.65695967542,1089.49837088430,1121.21697240827,1131.51942253980];
% 冷负荷
P_c_load= [560.667569473881,599.305129759678,601.402928695395,604.406267303822,639.083522597915,528.074031699563,614.040991490826,574.101828717455,655.498992194637,652.634828573160,770.093215102705,800.976839125784,820.413910851808,775.979118604715,741.892953141515,699.421969248419,737.469821850432,753.007801543912,711.675009348734,746.703455359319,688.281300673128,582.033015706014,531.873237458605,561.594694688039];
% 风、光预测出力
P_pv= [33.1720717442352,16.4592840977671,41.2310237491652,12.1089407662635,45.3684971776781,38.2739684412357,164.959887152324,368.854225425538,878.946323040853,1173.71785764822,1352.80584601190,1622.93596648321,1727.61945443386,1658.98413309822,1053.06269421854,553.860563968220,72.9396250912741,23.6706709972714,3.57075529813238,45.8151259079952,25.5899639890867,12.7284191464265,20.2283584573044,47.9582166141778];
P_wt= [1172.54500720798,1207.66430177826,1237.89928202681,1246.97942661048,1316.82079831985,1146.36054616807,1345.90082750747,1168.87444329790,925.205681677089,794.944485705988,772.871084354346,555.046179511916,431.010306522448,450.751740216153,475.425146714131,530.218067769728,614.392063479094,634.138967325128,762.646522250968,940.956234257272,1149.83755781112,1154.05536988983,1405.76536888628,1450.15849296105];
% 能源市场电力价格:
c_e_buy= [0.38*ones(1,7),0.68*ones(1,4),1.2*ones(1,3),0.68*ones(1,4),1.2*ones(1,4),0.38*ones(1,2)];
% 能源市场天然气价格:
c_g_buy= 3.7*ones(1,24);
%% 绿证-碳交易联合机制参数
lambda_G = 0.5; % 绿证单价 (元/本)
rho = 0.15; % 绿证配额系数
k_renew = 0.01; % 可再生能源-绿证换算系数 (本/kWh)
lambda_C = 0.268; % 碳交易单价 (元/kg)
beta = 0.5; % 碳排放配额转换系数 (kg/kW)
gamma = 0.2; % 燃气机组-碳排放换算系数 (kg/kWh)
k_offset = 1.0; % 绿证碳减排扣减系数
EF_grid = 0.8; % 基准排放因子 (kg/kWh)
EF_renew = 0.0; % 可再生能源排放因子 (kg/kWh)
%% 各变量及常量定义
P_g_CHP= sdpvar(1,24); % 输入CHP设备消耗天然气功率
P_CHP_e= sdpvar(1,24); % CHP设备的输出电功率
P_CHP_h= sdpvar(1,24); % CHP设备的输出热功率
P_g_GB= sdpvar(1,24); % 输入GB设备的天然气功率
P_GB_h= sdpvar(1,24); % GB设备输出的热功率
P_e_EB= sdpvar(1,24); % 输入EB设备的电功率
P_EB_h= sdpvar(1,24); % EB设备输出的热功率
P_e_EC= sdpvar(1,24); % 电制冷机EC的耗电功率
P_EC_c= sdpvar(1,24); % 电制冷机EC的制冷功率
P_h_AC=sdpvar(1,24); % 吸收式制冷机AC的耗热功率
P_AC_c=sdpvar(1,24); % 吸收式制冷机AC的制冷功率
% 电ES1、热ES2、冷ES3
P_ES1_cha=sdpvar(1,24);
P_ES2_cha=sdpvar(1,24);
P_ES3_cha=sdpvar(1,24);
% 充放功率
P_ES1_dis=sdpvar(1,24);
P_ES2_dis=sdpvar(1,24);
P_ES3_dis=sdpvar(1,24);
% 各储能的实时容量状态
S_1 = sdpvar(1,25);
S_2 = sdpvar(1,25);
S_3 = sdpvar(1,25);
% 充放标志二进制变量
B_ES1_cha=binvar(1,24);
B_ES2_cha=binvar(1,24);
B_ES3_cha=binvar(1,24);
% 充标志
B_ES1_dis=binvar(1,24);
B_ES2_dis=binvar(1,24);
B_ES3_dis=binvar(1,24);
% 能源市场购能约束
P_e_buy=sdpvar(1,24); % 购电功率
P_g_buy=sdpvar(1,24); % 购气功率
%% 约束条件
Constraints=[];
Constraints=[Constraints,
P_CHP_e==0.93*P_g_CHP, % CHP的气-电能量转换约束
P_CHP_h==0.73*P_g_CHP, % CHP的气-热能量转换约束
0<=P_g_CHP<=950, % CHP消耗的气功率上下限约束
-0.2*950<=P_g_CHP(2:24)-P_g_CHP(1:23)<=0.2*950, % CHP的爬坡约束(1-24时段)
];
Constraints=[Constraints,
P_GB_h==0.92*P_g_GB, % GB的气-热能量转换约束
0<=P_g_GB<=790, % GB的出力上下限约束
-0.2*800<=P_g_GB(2:24)-P_g_GB(1:23)<=0.2*800, % GB的爬坡约束(1-24时段)
];
Constraints=[Constraints,
P_EB_h==0.79*P_e_EB, % EB的电-热能量转换约束
0<=P_e_EB<=860, % EB的出力上下限约束
-0.2*860<=P_e_EB(2:24)-P_e_EB(1:23)<=0.2*860, % EB的爬坡约束(1-24时段)
];
Constraints=[Constraints,
P_EC_c==0.82*P_e_EC, % 电制冷机的电-冷的约束关系
0<=P_e_EC<=630, % 电制冷机上下限约束
-0.2*630<=P_e_EC(2:24)-P_e_EC(1:23)<=0.2*630, % EC的爬坡约束(1-24时段)
];
Constraints=[Constraints,
P_AC_c==0.63*P_h_AC, % 吸收式制冷机的热-冷的约束关系
0<=P_h_AC<=570, % 吸收式制冷机上下限约束
-0.2*570<=P_h_AC(2:24)-P_h_AC(1:23)<=0.2*570, % AC的爬坡约束(1-24时段)
];
Constraints=[Constraints,
0<=P_ES1_cha<=B_ES1_cha*450*2, % 储电设备的最大充电功率约束
0<=P_ES2_cha<=B_ES2_cha*500*2, % 储热设备的最大充热功率约束
0<=P_ES3_cha<=B_ES3_cha*350*2, % 储冷设备的最大充气功率约束
0<=P_ES1_dis<=B_ES1_dis*450*2, % 储电设备的最大放电功率约束
0<=P_ES2_dis<=B_ES2_dis*500*2, % 储热设备的最大放热功率约束
0<=P_ES3_dis<=B_ES3_dis*350*2, % 储冷设备的最大放气功率约束
S_1(1) == 0.3*450*2, % 初始状态
S_2(1) == 0.3*500*2,
S_3(1) == 0.3*350*2,
% 充放状态唯一
B_ES1_cha+B_ES1_dis<=1,
B_ES2_cha+B_ES2_dis<=1,
B_ES3_cha+B_ES3_dis<=1,
% 储能容量上下限约束
0.2*450*2 <= S_1(2:25) <= 0.8*450*2,
0.2*500*2 <= S_2(2:25) <= 0.8*500*2,
0.2*350*2 <= S_3(2:25) <= 0.8*350*2,
% 储能容量变化约束
S_1(2:25) == S_1(1:24) + 0.95*P_ES1_cha - P_ES1_dis/0.95,
S_2(2:25) == S_2(1:24) + 0.95*P_ES2_cha - P_ES2_dis/0.95,
S_3(2:25) == S_3(1:24) + 0.95*P_ES3_cha - P_ES3_dis/0.95
];
Constraints=[Constraints,
P_e_buy+P_wt+P_pv+P_CHP_e-P_e_EB-P_e_EC-P_ES1_cha+P_ES1_dis==P_e_load, % 电功率平衡约束
P_CHP_h+P_GB_h+P_EB_h-P_h_AC-P_ES2_cha+P_ES2_dis==P_h_load, % 热功率平衡约束
P_EC_c+P_AC_c-P_ES3_cha+P_ES3_dis==P_c_load, % 冷功率平衡约束
P_g_buy-P_g_GB-P_g_CHP==0, % 天然气功率平衡约束
0<=P_e_buy<=2000, % 购电功率约束
0<=P_g_buy<=2000, % 购气功率约束
];
%% 目标函数
% 绿证交易成本
N_GS = rho * sum(P_e_load);
N_GP = k_renew * sum(P_wt + P_pv);
C_GCT = lambda_G * (N_GS - N_GP);
% 碳交易成本
E_CS = beta * (1600); % CHP+GB装机容量
E_CP = gamma * (sum(P_g_CHP) + sum(P_g_GB));
E1 = 0.175 *[282.839969258212,1732.88303939672,1136.86541399768,1143.49431585863,1201.12679562320,1294.34994259177,366.005494600093,873.805241602914,693.509744882107,605.368456793456,365.639876711029,375.091247566023,238.464953397028,1031.72231573773,997.704198854051,1452.20246701193,2084.23506301809,1978.66973668847,1879.35968001180,1043.65763169695,1367.84703767313,1201.11363215881,936.277288668493,788.330252657417];
E2 = [36.6416896601208,156.279608587218,124.375098157906,67.7682313766277,178.894097016925,48.9147825719602,50.2727164336131,32.0706019551670,51.1095317009999,75.7459485300551,57.0564526390624,6.08591317072405,44.3251964509820,160.064016224196,97.1348111231326,126.590983521110,186.240471229658,121.232885825931,191.134132772061,106.614128200402,223.677927749615,190.936570017364,120.652086395541,59.6938460582076];
E3 = [35.2891913717701,85.4545836977783,64.8163201098496,15.6947419839914,87.4609204322101,30.5257635248543,34.1419884659322,12.6844110489195,18.7795150368931,74.8356339601018,2.15324583424892,5.38705563942302,40.4816180512034,127.440388458987,9.58839854274324,33.1505307220314,62.4570159047602,82.4053811220228,26.0999654857227,76.8930409954489,23.8802719669989,124.826185778421,59.6231164978258,46.5045936379491];
E4 = [21.6169753695504,84.5957827103781,34.2022802898281,7.52599580616450,70.0865991650976,6.95507957617122,17.0059295449738,11.4267832589820,10.7918590057747,63.2494452877671,1.59047413174750,3.15674476597780,9.98821284100896,84.9283416105617,0.800466488291493,20.7508990890867,41.2806247678819,60.1354797619188,23.2485994924254,75.5322819871507,18.3646237083536,72.5797473157448,55.3489178608123,26.9768667328856];
% 绿证-碳交易联合成本
E_offset = k_offset * N_GP * (EF_grid - EF_renew);
C_CET_all = lambda_C * (E_CS - E_CP - E_offset);
% 运维成本
C_om=0;
for t=1:24
C_om=C_om+P_g_CHP(t)*0.025+P_g_GB(t)*0.025+P_e_EB(t)*0.025+P_e_EC(t)*0.025+P_h_AC(t)*0.025+P_wt(t)*0.018+P_pv(t)*0.018;
end
% 购能成本
Cost_e=0;
Cost_g=0;
for t=1:24
Cost_e=Cost_e+c_e_buy(t)*P_e_buy(t); % 购电成本
end
for t=1:24
Cost_g=Cost_g+c_g_buy(t)*P_g_buy(t); % 购气成本
end
Cost_buy= Cost_e+Cost_g;
% 总目标函数
obj = Cost_buy + C_om + C_CET_all; % 购能成本+ 运维成本+ 联合交易成本
%% 模型求解
ops=sdpsettings('solver','cplex','verbose',0,'usex0',0);
ops.cplex.mip.tolerances.mipgap=1e-6;
result=optimize(Constraints,obj,ops);
%% 打印输出结果
disp([' 上网购能成本: ', num2str(value(Cost_buy)), ' 元']);
disp([' 机组运维成本: ', num2str(value(C_om)), ' 元']);
disp([' 联碳交易成本: ', num2str(value(C_CET_all)), ' 元']);
disp([' 总规划运行成本: ', num2str(value(obj)), ' 元']);
%% 绘图
figure(1)
x=1:24;
plot(x,P_e_load,'-go');
hold on
plot(x,P_h_load,'-r*');
hold on
plot(x,P_c_load,'-b+');
xlabel('Time/h');
ylabel('Power/kW');
grid on
title('冷、热、电负荷需求曲线');
legend('电负荷', '热负荷', '冷负荷');
grid on;
figure(2)
plot(P_pv,'-b*');
hold on
plot(P_wt,'-ro');
hold off
grid on
xlabel('时间/h');
ylabel('功率/kW');
title('风、光机组输出功率预测曲线');
legend('PV机组','WT机组');
grid on;
figure(3); % 电功率
Figure_E_in = [P_e_buy', P_CHP_e', P_wt', P_pv', P_ES1_dis']; % 产电量
bar(Figure_E_in, 'stacked');
hold on
Figure_E_out = [-P_ES1_cha', -P_e_EB', -P_e_EC']; % 耗电量
bar(Figure_E_out, 'stacked');
hold on
plot(P_e_load,'g-*','LineWidth',1);
applyhatch(gcf, '/\x.');
xlabel('时间/h');
ylabel('功率/kW');
title('电功率平衡结果');
legend('上网购电', 'CHP机组(电)', '风电机组', '光伏机组', '储能放电', '储能充电', '电锅炉', '电制冷机', '电负荷');
grid on;
axis([0 25,-2000,4000])
figure(4) % 热功率
Figure_H_in = [P_CHP_h',P_GB_h', P_EB_h', P_ES2_dis'];
bar(Figure_H_in, 'stacked');
hold on
Figure_H_out = [-P_ES2_cha', -P_h_AC'];
bar(Figure_H_out, 'stacked');
hold on
plot(P_h_load,'r-*','LineWidth',1);
xlabel('时间/h');
ylabel('功率/kW');
title('热功率平衡结果');
legend('CHP机组(热)', '燃气锅炉', '电锅炉', '储能放热', '储能充热', '吸收制冷机', '热负荷');
grid on;
axis([0 25,-500,2000])
figure(5) % 冷功率
Figure_C_in = [P_AC_c',P_EC_c', P_ES3_dis'];
bar(Figure_C_in, 'stacked');
hold on
Figure_C_out = [-P_ES3_cha'];
bar(Figure_C_out, 'stacked');
hold on
plot(P_c_load,'b-*','LineWidth',1);
xlabel('时间/h');
ylabel('功率/kW');
title('冷功率平衡结果');
legend('吸收制冷机','电制冷机', '储能放冷', '储能充冷', '冷负荷');
grid on;
axis([0 25,-200,1200])