请在场景1代码的基础上修改场景2,优化变量全部改成场景1的,请给出修改后的场景2代码,即场景1补充CCS后的代码:%% 场景1:考虑碳交易+P2G
clc
clear
%% 加载电负荷、热负荷
load date.mat;
P_load_e = date(1, :); % 总电负荷
H_h = date(2, :); % 高品位热负荷
H_l = date(3, :); % 低品位热负荷
H_load_h = H_h + H_l; % 总热负荷
% 成本系数和能源价格
c_coal = [0.0004818, 2.695, 52, 0.0004751, 0.3536, 0.001004] * 0.1471;
p_coal = 800; % 标煤价格 (元/吨)
p_gas = 3.2; % 天然气价格 (元/m³)
p_buy_e = 1000 * [0.17, 0.17, 0.17, 0.17, 0.17, 0.17, 0.17, 0.49, 0.83, 0.83, 0.83, 0.83, 0.49, 0.49, 0.49, 0.49, 0.49, 0.83, 0.83, 0.83, 0.83, 0.49, 0.49, 0.49]; % 电价系数
V_buy_h = 5 * [10, 10.4, 11.3, 11.7, 12.4, 13.2, 14.1, 14.6, 15.1, 15.3, 14.7, 14.4, 14, 13.8, 13.3, 12.7, 12.2, 12.3, 11.9, 11.8, 11.4, 11.1, 10.4, 10.2]; % 热价基数
C_sell_h = 180; % 售热价格 (元/MWh)
% 气象数据(辐照度 kW/m²,风速 m/s)
P_Solar = [0,0,0,0,0,0.09,1.15,2.82,4.54,6.20,7.44,8.01,7.93,7.09,5.93,4.57,2.97,1.27,0.15,0,0,0,0,0];
P_Wind = [7.14,7.20,7.22,6.89,7.18,7.07,6.43,5.97,6.33,6.56,6.63,6.70,6.65,6.92,6.86,6.85,6.89,7.01,7.02,6.97,7.05,7.06,7.03,7.10];
%% 新能源模型参数
% 光伏模型 (单位: MW)
S_PV = 40000 * 3.75; % 光伏面积 (m²)
k_eta_PV = 0.157; % 光电转换效率
P_PV_max = zeros(1, 24);
for i = 1:24
P_PV_max(i) = (S_PV * P_Solar(i) * k_eta_PV) / 1000;
end
% 风电模型 (单位: MW)
Q_S_WT = 80000 * 4.8; % 单机额定容量 (kW)
v_cut_in = 3; v_rated = 12; v_cut_out = 20; % 风速参数
P_WT_max = zeros(1, 24);
for i = 1:24
v = P_Wind(i);
if v >= v_cut_in && v <= v_rated
P_WT_max(i) = Q_S_WT * (v - v_cut_in) / (v_rated - v_cut_in) / 1000;
elseif v > v_rated && v <= v_cut_out
P_WT_max(i) = Q_S_WT / 1000;
else
P_WT_max(i) = 0;
end
end
%% 定义决策变量
P_th_h = sdpvar(1, 24, 'full'); % 火电机组热出力
P_th_e = sdpvar(1, 24, 'full'); % 火电机组电出力
z_th_1 = binvar(1, 24); % 机组1状态
z_th_2 = binvar(1, 24); % 机组2状态
P_e_eb = sdpvar(1, 24, 'full'); % 电锅炉电功率
P_g_gb = sdpvar(1, 24, 'full'); % 燃气锅炉燃气功率
P_g_chp = sdpvar(1, 24, 'full'); % CHP机组燃气功率
B_cr_h = binvar(1, 24); % 储热状态 (0:储热, 1:放热)
P_loss_e = sdpvar(1, 24, 'full'); % 失负荷功率
m_ccs_st = sdpvar(1, 24, 'full'); % CCS解吸量
m_ccs_ab = sdpvar(1, 24, 'full'); % CCS吸收量
P_e_p2g = sdpvar(1, 24, 'full'); % P2G用电功率
P_h2_chp = sdpvar(1, 24, 'full'); % CHP掺氢功率
level_cool = binvar(3, 24); % COAP冷却级别
Z_coap_mt = sdpvar(1, 24, 'full'); % 活性炭质量
P_e_coap = sdpvar(1, 24, 'full'); % COAP能耗
B_coap_x1 = sdpvar(2, 24, 'full'); % COAP辅助变量
P_pv_e = sdpvar(1, 24, 'full'); % 光伏上网功率
P_wt_e = sdpvar(1, 24, 'full'); % 风电上网功率
P_ess_e = sdpvar(1, 24, 'full'); % 电化学储能净功率
P_ess_cha = sdpvar(1, 24, 'full'); % 储能充电功率
P_ess_dis = sdpvar(1, 24, 'full'); % 储能放电功率
H_tes_dis = sdpvar(1, 24, 'full'); % 储热罐释热功率
H_tes_cha = sdpvar(1, 24, 'full'); % 储热罐储热功率
SOC_HSS = sdpvar(1,24);
%% 系统参数
SOC_HSS_max = 80; % 储热罐最大容量 (MWh)
SOC_ESS_max = 40; % 储电能最大容量 (MWh)
P_eb_max = 10; % 电锅炉最大功率 (MW)
P_gb_max = 40; % 燃气锅炉最大功率 (MW)
k_eb = 0.99; % 电锅炉效率
k_gb = 0.81; % 燃气锅炉效率
K_chp_h = 0.47; % CHP机组气-热比
K_chp_e = 0.35; % CHP机组气-电比
%% 约束条件
n = 24;
Cons = [];
for i = 1:n
% 火电热功率约束
Cons = [Cons, (0 <= P_th_h(1,i)) & (P_th_h(1,i) <= 2*98.18)];
Cons = [Cons, implies(P_th_h(1,i) >= 2*57.27, z_th_1(1,i) == 1)];
Cons = [Cons, z_th_1(1,i) + z_th_2(1,i) == 1];
% 火电电功率约束
Cons = [Cons, (2*43.37 <= P_th_e(1,i)) & (P_th_e(1,i) <= 2*135)];
% 火电热-电耦合约束
Cons = [Cons, implies(z_th_2(1,i), (P_th_e(1,i) <= 2*135 - 0.25*P_th_h(1,i)) & (P_th_e(1,i) >= 2*54 - 0.1857*P_th_h(1,i)))];
Cons = [Cons, implies(z_th_1(1,i), (P_th_e(1,i) <= 2*135 - 0.25*P_th_h(1,i)) & (P_th_e(1,i) >= 1.64*P_th_h(1,i) - 2*50.56))];
% CHP机组约束
Cons = [Cons, (0 <= P_g_chp(1,i)) & (P_g_chp(1,i) <= 2*50)];
% 电锅炉约束
Cons = [Cons, (0 <= P_e_eb(1,i)) & (P_e_eb(1,i) <= P_eb_max)];
% 失负荷约束
Cons = [Cons, P_loss_e(1,i) >= 0];
% 储热系统功率约束
Cons = [Cons, (0 <= H_tes_cha(1,i)) & (H_tes_cha(1,i) <= 20)];
Cons = [Cons, (0 <= H_tes_dis(1,i)) & (H_tes_dis(1,i) <= 20)];
% 储热容量约束 (动态)
if i == 1
initial_H = 0.3 * SOC_HSS_max;
else
initial_H = 0;
end
Cons = [Cons, (0 <= initial_H + sum(0.98*H_tes_cha(1,1:i) + k_eb*P_e_eb(1,1:i) - H_tes_dis(1,1:i)) <= SOC_HSS_max)];
% 储热状态互斥
Cons = [Cons, implies(B_cr_h(i) == 0, H_tes_dis(1,i) == 0)];
Cons = [Cons, implies(B_cr_h(i) == 1, H_tes_cha(1,i) == 0)];
% 电功率平衡
Cons = [Cons, (P_th_e(1,i) + K_chp_e*P_g_chp(1,i) + P_wt_e(1,i) + P_pv_e(1,i) + P_loss_e(1,i) - P_e_eb(1,i) - P_e_p2g(1,i) - P_e_coap(1,i) + P_ess_e(1,i) - (5 + 0.269*m_ccs_st(1,i))) == P_load_e(i)];
% 热功率平衡
Cons = [Cons, P_th_h(1,i) + 0.98*H_tes_dis(1,i) - H_tes_cha(1,i) + K_chp_h*P_g_chp(1,i) + k_gb*P_g_gb(1,i) == H_load_h(i)];
% 燃气锅炉约束
Cons = [Cons, (0 <= P_g_gb(1,i)) & (P_g_gb(1,i) <= P_gb_max)];
% CCS约束 (当前禁用)
Cons = [Cons, m_ccs_st(1,i) == 0];
Cons = [Cons, m_ccs_ab(1,i) == 0];
% P2G约束
Cons = [Cons, (0 <= P_e_p2g(1,i)) & (P_e_p2g(1,i) <= 30)];
Cons = [Cons, (0 <= P_h2_chp(1,i)) & (P_h2_chp(1,i) <= P_e_p2g(1,i))];
Cons = [Cons, (P_g_chp(1,i) + P_g_gb(1,i))*0.003539 >= 5*P_h2_chp(1,i)*0.01083]; % 掺氢比约束
% COAP冷却级别 (固定使用第三级)
Cons = [Cons, level_cool(1,i) + level_cool(2,i) + level_cool(3,i) == 1];
Cons = [Cons, level_cool(3,i) == 1];
Cons = [Cons, P_e_coap(1,i) == 0]; % 所有级别能耗均为0
% 活性炭约束
Cons = [Cons, (0 <= Z_coap_mt(1,i)) & (Z_coap_mt(1,i) <= 20)];
% 电化学储能约束
Cons = [Cons, (-20 <= P_ess_e(1,i)) & (P_ess_e(1,i) <= 20)];
Cons = [Cons, (-20 <= P_ess_cha(1,i)) & (P_ess_cha(1,i) <= 0)];
Cons = [Cons, (0 <= P_ess_dis(1,i)) & (P_ess_dis(1,i) <= 20)];
Cons = [Cons, implies(P_ess_e(1,i) >= 0, (P_ess_dis(1,i) == P_ess_e(1,i)) & (P_ess_cha(1,i) == 0))];
Cons = [Cons, implies(P_ess_e(1,i) <= 0, (P_ess_dis(1,i) == 0) & (P_ess_cha(1,i) == P_ess_e(1,i)))];
% 储能容量约束 (动态)
if i == 1
initial_ESS_SOC = 0.3 * SOC_ESS_max;
else
initial_ESS_SOC = 0;
end
Cons = [Cons, (0 <= initial_ESS_SOC + sum(P_ess_dis(1,1:i))/0.95 + 0.95*sum(P_ess_cha(1,1:i)) <= SOC_ESS_max)];
% COAP辅助变量约束
Cons = [Cons, B_coap_x1(1,i) == (P_th_e(1,i) + P_th_h(1,i)) * 6.48];
Cons = [Cons, B_coap_x1(2,i) == (P_th_e(1,i) + P_th_h(1,i)) * 2.88];
% 新能源出力约束
Cons = [Cons, (0 <= P_pv_e(1,i)) & (P_pv_e(1,i) <= P_PV_max(i))];
Cons = [Cons, (0 <= P_wt_e(1,i)) & (P_wt_e(1,i) <= P_WT_max(i))];
end
% 燃气轮机爬坡约束
for i = 2:n
Cons = [Cons, (2*-28.6 <= P_g_chp(1,i) - P_g_chp(1,i-1)) & (P_g_chp(1,i) - P_g_chp(1,i-1) <= 2*28.6)];
end
% 储热和储能周期平衡约束
Cons = [Cons, (-1 <= sum(0.98*H_tes_cha(1,1:24) + k_eb*P_e_eb(1,1:24) - H_tes_dis(1,1:24)) <= 1)];
Cons = [Cons, (-1 <= sum(P_ess_dis(1,1:24))/0.95 + 0.95*sum(P_ess_cha(1,1:24)) <= 1)];
%% 燃料购能成本
Fun = 0;
for i = 1:n
% 火电机组燃料成本
C_th_c = p_coal * (c_coal(3) + c_coal(2)*P_th_e(i) + c_coal(5)*P_th_h(i) + c_coal(1)*P_th_e(i)^2 + c_coal(6)*P_th_e(i)*P_th_h(i) + c_coal(4)*P_th_h(i)^2);
% CHP机组和燃气锅炉燃气成本
C_chp_c = p_gas * (P_g_chp(i) / 0.01083);
C_gb_c = p_gas * (P_g_gb(i) / 0.01083);
% 环保代价成本
C_cb_c = 6.4 * B_coap_x1(1,i) + 8 * B_coap_x1(2,i) + 3.16;
% 碳排放计算
C_mc = 2.57 * (C_th_c/p_coal) + 0.054*0.98*3.67*(P_g_chp(i) + P_g_gb(i) - P_h2_chp(i));
% 碳配额计算
C_Ec = 0.69135*(P_th_e(i) + K_chp_e*P_g_chp(i)) + 0.3*(P_th_h(i) + K_chp_h*P_g_chp(i) + k_gb*P_g_gb(i));
% 碳交易成本
C_ct_c = 205 * (C_mc - C_Ec);
% 失负荷惩罚成本
C_buy_c = 5 * p_buy_e(i) * P_loss_e(i);
% 售热售电收益
I_sale = p_buy_e(i) * (P_th_e(i) + K_chp_e*P_g_chp(i) + P_wt_e(i) + P_pv_e(i) - P_e_eb(i) - P_e_p2g(i) - P_e_coap(i) + P_ess_e(i)) + C_sell_h * (P_th_h(i) + 0.98*H_tes_dis(i) - H_tes_cha(i) + K_chp_h*P_g_chp(i) + k_gb*P_g_gb(i));
% P2G产气收益
I_ch4_c = p_gas * (((P_e_p2g(i)*0.85 - P_h2_chp(i)) * 0.75) / 0.01083);
% 总成本
C0_all = C_th_c + C_gb_c + C_chp_c - p_gas*P_h2_chp(i)/(0.01083) + C_cb_c + C_buy_c + C_ct_c - I_ch4_c;
% 目标函数累加
Fun = Fun + (I_sale - C0_all);
end
%% 设备投资和运维成本
% 火电建设成本
C_Ld = (4589000*((135*2)*(1+0.07)) + 100) / (365*30);
C_cs = 0;
for i=1:n
% 累加燃料相关成本
fuel_cost_i = C_th_c + C_gb_c + C_chp_c - p_gas*P_h2_chp(i)/0.01083 - I_ch4_c;
C_cs = C_cs + fuel_cost_i;
end
% 火电运维成本
C_com = (0.1536/(1-0.1536)) * (C_Ld + C_cs);
% 光伏成本(建设+运维)
C_pv = (6000*1000*200)/(365*25) + 12*sum(P_PV_max(1:24));
% 风电成本(建设+运维)
C_wt = (8000*1000*440)/(365*20) + 20*sum(P_WT_max(1:24));
% 电锅炉成本
C_eb = (3200*1000*P_eb_max)/(365*20) + 2*sum(P_e_eb(1,1:24));
% 储热罐成本
C_yan = (300*1000*SOC_HSS_max)/(365*5) + 1.5*sum(H_tes_cha(1,1:24));
% 电化学储能成本
C_ess = (2000*1000*SOC_ESS_max)/(365*5);
% 煤电容量电价补偿
I_c_e = (330*1000*(370))/365;
% 总设备成本
C_q_om = C_Ld + C_com + C_pv + C_wt + C_eb + C_yan + C_ess - I_c_e;
%% 弃风弃光惩罚
P_v_a = P_PV_max - P_pv_e; % 光伏弃光量
P_w_a = P_WT_max - P_wt_e; % 风电弃风量
F_a_em = 257.396 * sum(P_w_a + P_v_a); % 惩罚成本
%% 目标函数
Fun = Fun - C_q_om - F_a_em; % 总收益减去设备成本和弃能惩罚
Fun = -Fun; % 转换为最小化问题
%% CPLEX优化求解
ops = sdpsettings('solver', 'cplex', 'verbose', 2);
ops.cplex.timelimit = 120; % 设置求解时间限制为120秒
optimize(Cons, Fun, ops);
%% 结果可视化
figure(1)
x = 1:24;
hold on;
P_pv_e_val = value(P_pv_e);
P_wt_e_val = value(P_wt_e);
% 初始化弃能区域标记
has_pv_curtail = false;
has_wt_curtail = false;
% 绘制弃能区域
for i = 1:24
% 光伏弃光区域
if P_PV_max(i) > P_pv_e_val(i)
patch([x(i)-0.4, x(i)+0.4, x(i)+0.4, x(i)-0.4], ...
[P_pv_e_val(i), P_pv_e_val(i), P_PV_max(i), P_PV_max(i)], ...
[0.97, 0.75, 0.40], 'EdgeColor', 'none', 'FaceAlpha', 0.5);
has_pv_curtail = true;
end
% 风电弃风区域
if P_WT_max(i) > P_wt_e_val(i)
patch([x(i)-0.4, x(i)+0.4, x(i)+0.4, x(i)-0.4], ...
[P_wt_e_val(i), P_wt_e_val(i), P_WT_max(i), P_WT_max(i)], ...
[0.40, 0.70, 0.95], 'EdgeColor', 'none', 'FaceAlpha', 0.5);
has_wt_curtail = true;
end
end
% 绘制实际出力曲线
p_pv_actual = plot(x, P_pv_e_val, '-o', 'Color', [0.90, 0.55, 0.20], 'LineWidth', 1, 'MarkerSize', 5);
p_wt_actual = plot(x, P_wt_e_val, '-o', 'Color', [0.20, 0.50, 0.80], 'LineWidth', 1, 'MarkerSize', 5);
% 绘制最大出力曲线
p_pv_max = plot(x, P_PV_max, '-s', 'Color', [0.80, 0.35, 0.10], 'LineWidth', 1);
p_wt_max = plot(x, P_WT_max, '-s', 'Color', [0.10, 0.30, 0.60], 'LineWidth', 1);
% 创建图例项
legend_items = [p_pv_actual, p_pv_max];
legend_labels = {'Actual PV Output', 'Maximum PV Output'};
if has_pv_curtail
legend_items(end+1) = patch(NaN, NaN, [0.97, 0.75, 0.40], 'FaceAlpha', 0.5);
legend_labels{end+1} = 'Curtailed PV Power';
end
legend_items(end+1:end+2) = [p_wt_actual, p_wt_max];
legend_labels(end+1:end+2) = {'Actual WT Output', 'Maximum WT Outpu'};
if has_wt_curtail
legend_items(end+1) = patch(NaN, NaN, [0.40, 0.70, 0.95], 'FaceAlpha', 0.5);
legend_labels{end+1} = 'Curtailed WT Power';
end
legend(legend_items, legend_labels);
grid on;
xlabel('time (h)');
ylabel('Power (MW)');
title('Output Curve of New Energy Units');
xlim([0, 25]);
ylim([-1 260]);
set(gca, 'XTick', 0:5:25);
box on;
hold off;
figure(2)
% Electric Power Balance
subplot(1,3,1);
hold on;
P = [P_g_chp*K_chp_e; P_th_e; P_pv_e; P_wt_e; -P_e_eb; -P_e_p2g; P_ess_cha; P_ess_dis]';
e = bar(P, 'stacked');
hold on;
plot(1:24, P_load_e, '-g*');
grid on;
colors_e = [
0.7 0.2 0.7 % CHP (purple)
0.9 0.4 0.2 % Thermal (orange)
1 0.8 0.2 % PV (yellow)
0.5 0.8 0.2 % WT (green)
0.6 0.6 0.6 % EB (gray)
0.8 0.3 0.3 % P2G (red)
0.2 0.8 0.8 % ESS Charge (cyan)
0.3 0.7 0.9 % ESS Discharge (light blue)
];
hold on;
for i = 1:size(colors_e,1)
e(i).FaceColor = colors_e(i,:);
e(i).EdgeColor = 'none';
end
hold on;
legend_items = {'CHP','CTP','PV','WT','EB','P2G','ESS(cha)','ESS(dis)','Electric Load'};
legend(legend_items);
xlabel('Time (h)');
ylabel('Power (MW)');
axis([0,25,-200,700]);
xticks(0:2:24);
title('Electric Power Balance');
% Thermal Power Balance
subplot(1,3,2);
H1 = [K_chp_h*P_g_chp; P_th_h; k_gb*P_g_gb; 0.98*H_tes_dis; -(H_tes_cha)]';
h = bar(H1, 'stacked');
hold on;
plot(1:24, H_load_h, '-r*');
grid on;
colors_h = [
0.7 0.2 0.7 % CHP heat
0.9 0.4 0.2 % Thermal heat
0.3 0.6 1 % GB heat
0.5 0.8 0.2 % TES discharge
0.8 0.3 0.3 % TES charge
];
for i = 1:size(colors_h,1)
h(i).FaceColor = colors_h(i,:);
h(i).EdgeColor = 'none';
end
hold on;
legend('CHP','CTP','GB','TES(dis)','TES(cha)','Thermal Load');
xlabel('Time (h)');
ylabel('Power (MW)');
axis([0,25,-50,290]);
xticks(0:2:24);
title('Thermal Power Balance');
% Nature Gas Balance
subplot(1,3,3);
V_chp1 = (P_g_chp(1,:));
V_gb = (P_g_gb);
V_ch4 = (((P_e_p2g*0.85-P_h2_chp)*0.75));
V_buy_g = V_chp1 + V_gb - V_ch4;
V = [V_chp1; V_gb; -V_ch4]';
g = bar(V, 'stacked');
hold on;
plot(1:24, value(V_buy_g), '-b*');
grid on;
colors_g = [
0.7 0.2 0.7 % CHP
0.3 0.6 1 % GB
0.8 0.3 0.3 % P2G
];
hold on;
for i = 1:size(colors_g,1)
g(i).FaceColor = colors_g(i,:);
g(i).EdgeColor = 'none';
end
legend('CHP','GB','P2G','Gas Purchase');
xlabel('Time (h)');
ylabel('Gas Flow (MW)');
axis([0,25,-50,200]);
xticks(0:2:24);
title('Nature Gas Balance');
figure(3);
E_SOC_elec(1)=SOC_ESS_max*0.3;
H_SOC_heat(1)=SOC_HSS_max*0.3;
for i=2:25
E_SOC_elec(i)=0.3*SOC_ESS_max+sum(P_ess_dis(1,1:i-1))/0.95+0.95*sum(P_ess_cha(1,1:i-1));
H_SOC_heat(i)=0.3*SOC_HSS_max+sum(0.98*H_tes_cha(1,1:i-1)+k_eb*P_e_eb(1,1:i-1)-H_tes_dis(1,1:i-1));
end
% 子图1:电储能状态变化
subplot(2,1,1);
P1 = [P_ess_cha; P_ess_dis]';
bar(P1, 'stacked');
hold on
plot(0:24, E_SOC_elec, 'g*-');
xlabel('Time (h)');
ylabel('ESS (MWh)');
title('Electric Energy Storage State Variation');
grid on;
ylim([-30, 60]);
xticks(0:2:24);
% 子图2:热储能状态变化
subplot(2,1,2);
H1 = [-H_tes_cha; H_tes_dis]';
bar(H1, 'stacked');
hold on
plot(0:24, H_SOC_heat, 'r*-');
xlabel('Time (h)');
ylabel('TES (MWh)');
title('Thermal Energy Storage State Variation');
grid on;
ylim([-30, SOC_HSS_max]);
xticks(0:2:24);
% 火电机组运行轨迹
figure(4);
% 定义四边形可行域顶点坐标
x_vertex = [0, 114.54, 196.36, 0];
y_vertex = [108, 86.73, 220.91, 270];
% 填充可行域(SCI素雅配色)
fill(x_vertex, y_vertex,[0.8 0.9 0.8], 'EdgeColor', 'none');
hold on;
% 绘制运行点轨迹(SCI专业配色,优化标记尺寸)
plot(value(P_th_h), value(P_th_e), 'o-', 'Color', [0.80, 0.35, 0.10], 'MarkerSize', 6);
xlabel('Heat Output (MW)');
ylabel('Power Output (MW)');
title('Thermal Power Unit Operation Trajectory');
legend('Feasible Region', 'Operation Points');
grid on;
axis([-10, 210, 70, 280]);
figure(5)
for i=2:25
m_th_so2(i-1)=(P_th_e(1,i-1)+P_th_h(1,i-1))*6.48;
m_th_nox(i-1)=(P_th_e(1,i-1)+P_th_h(1,i-1))*2.88;
m_coap_so2(i-1)=-B_coap_x1(1,i-1)+m_th_so2(i-1);
m_coap_nox(i-1)=-B_coap_x1(2,i-1)+m_th_nox(i-1);
end
level_cool=zeros(24,1);
color_so2_gen = [0.80, 0.35, 0.10];
color_so2_ads = [0.90, 0.55, 0.20];
color_nox_gen = [0.10, 0.30, 0.60];
color_nox_ads = [0.20, 0.50, 0.80];
color_level = [0.8 0.9 0.8];
color_power = [0.10, 0.50, 0.20];
subplot(2,1,1);
yyaxis left;
plot(1:24, value(m_th_so2), '-', 'Color', color_so2_gen, 'DisplayName', 'SO₂ Generation Amount');
hold on;
plot(1:24, value(m_coap_so2), '--', 'Color', color_so2_ads, 'DisplayName', 'SO₂ Adsorption Capacity');
ylabel('SO₂ mass (t)');
yyaxis right;
plot(1:24, value(m_th_nox), '-', 'Color', color_nox_gen, 'DisplayName', 'NOₓ Generation Amount');
plot(1:24, value(m_coap_nox), '--', 'Color', color_nox_ads, 'DisplayName', 'NOₓ Adsorption Capacity');
ylabel('NOₓ mass (t)');
title('Pollutant Gas Treatment');
xlabel('Time (h)');
legend('Location', 'best');
grid on;
xticks(0:2:24);
subplot(2,1,2);
yyaxis left;
bar(1:24, level_cool, 0.4, 'FaceColor', color_level, 'EdgeColor', 'none', 'DisplayName', 'Operation Level (Not activated)');
ylabel('Cooling Stage');
ylim([0, 5]);
yyaxis right;
plot(1:24, value(P_e_coap), '-*', 'Color', color_power, 'MarkerSize', 6, 'DisplayName', 'COAP Power Consumption');
ylabel('Power (MW)');
title('COAP Operation Status');
xlabel('Time (h)');
legend('Location', 'best');
grid on;
xticks(0:2:24);
%% 场景2:考虑碳交易+P2G+CCS
%% 初始化
close all
clear
clc
warning off
load date.mat
%加载电负荷、热负荷
c=[0.0004818,2.695,52,0.0004751,0.3536,0.001004]*0.1471;
p_coal=800;
p_gas=3.2;
%煤耗特性热电联产/火电(注:除1000,将千瓦时转化为兆瓦时)/一吨标煤价格
a=1000*[0.17,0.17,0.17,0.17,0.17,0.17,0.17,0.49,0.83,0.83,0.83,0.83,0.49,0.49,0.49,0.49,0.49,0.83,0.83,0.83,0.83,0.49,0.49,0.49];
V1=5*[10,10.4,11.3,11.7,12.4,13.2,14.1,14.6,15.1,15.3,14.7,14.4,14,13.8,13.3,12.7,12.2,12.3,11.9,11.8,11.4,11.1,10.4,10.2];
b=180;
%售电售热价格
pi(1)=1;
% 气象数据(辐照度 kW/m²,风速 m/s)
P_Solar = [0,0,0,0,0,0.09,1.15,2.82,4.54,6.20,7.44,8.01,7.93,7.09,5.93,4.57,2.97,1.27,0.15,0,0,0,0,0];
P_Wind = [7.14,7.20,7.22,6.89,7.18,7.07,6.43,5.97,6.33,6.56,6.63,6.70,6.65,6.92,6.86,6.85,6.89,7.01,7.02,6.97,7.05,7.06,7.03,7.10];
% 光伏模型 (单位: MW)
S_PV = 40000*3.75; % 光伏面积 (m²)
k_eta_PV = 0.157; % 光电转换效率
P_PV_max = zeros(1,24);
for i = 1:24
P_PV_max(i) = (S_PV * P_Solar(i) * k_eta_PV) / 1000;
end
% 风电模型 (单位: MW)
Q_S_WT = 80000*4.8; % 单机额定容量 (kW)
v_cut_in = 3; v_rated = 12; v_cut_out = 20; % 风速参数
P_WT_max = zeros(1,24);
for i = 1:24
v = P_Wind(i);
if v >= v_cut_in && v <= v_rated
P_WT_max(i) = Q_S_WT * (v - v_cut_in) / (v_rated - v_cut_in) / 1000;
elseif v > v_rated && v <= v_cut_out
P_WT_max(i) = Q_S_WT / 1000;
else
P_WT_max(i) = 0;
end
end
nt=1;
%场景生成、削减,得到风光出力场景/生成nt个场景
%% 定义变量
H_th_h=sdpvar(1,24,'full');
%CHP机组热出力
P_th_e=sdpvar(1,24,'full');
%火电机组电出力
H_tes_dis=sdpvar(1,24,'full');
%储热罐储能热出力
H_tes_cha=sdpvar(1,24,'full');
%储热罐储热量
P_eb_e=sdpvar(1,24,'full');
%电锅炉出力
P_gb_g=sdpvar(1,24,'full');
%电锅炉出力
P_chp=sdpvar(1,24,'full');
%CHP机组出力
z_chp=binvar(1,24);
z_chp2=binvar(1,24);
%CHP机组热电状态
cr=binvar(1,24);
%储热充放能状态,0充,1放;
P_loss_e=sdpvar(1,24,'full');
%失负荷,5倍价格购电
m_ccs_st=sdpvar(1,24,'full');
%ccs解吸塔处理量
m_ccs_ab=sdpvar(1,24,'full');
%ccs吸收处理量
P_p2g_e=sdpvar(1,24,'full');
%P_g2p功率
P_chp_H2=sdpvar(1,24,'full');
% Y_h2=sdpvar(1,24,'full'); B=Y_h2/(1-Y_h2)
% 燃气掺氢比
level_cool=binvar(3,24);
% COAP冷却塔温度,进行cool级冷却(一级20、二级0、三级-20)
coap_mt=sdpvar(1,24,'full');
% 活性炭质量
P_coap_e=sdpvar(1,24,'full');
% coap能耗
P_ess_e=sdpvar(1,24,'full');
% ES功率
P_pv_e=sdpvar(1,24,'full');
%光伏上网功率
P_wt_e=sdpvar(1,24,'full');
%风电上网功率
P_ess_cha=sdpvar(1,24,'full');
%储能充电功率
P_ess_dis=sdpvar(1,24,'full');
%储能放电功率
coap_x1=sdpvar(2,24,'full');
%coap辅助变量
P_load_e=date(1,:);
H_h=date(2,:);
H_l=date(3,:);
H_load_h=H_h+H_l;
SOC_HSS_max=80;
SOC_ESS_max=40;
E_SOC_elec(1)=SOC_ESS_max*0.3;
H_SOC_heat(1)=SOC_HSS_max*0.3;
Pebmax=10;
Pgbmax=40;
k_eb=0.99;
k_gb=0.8;
K_chp_H=0.4;
K_chp_P=0.35;
k_ccs_co2=1/(25*1.977*0.001);
v_ccs0=6000;
%% 约束条件
n=24;
C=[];
for i=1:n
C=[C,(0<=H_th_h(1,i))&(H_th_h(1,i)<=2*98.18)];
% 火电热功率约束
C=[C,implies(H_th_h(1,i)>=2*57.27,z_chp(1,i)==1)];
C=[C,z_chp(1,i)+z_chp2(1,i)==1];
C=[C,implies(z_chp2(1,i),(P_th_e(1,i)<=2*135-0.25*H_th_h(1,i))&(P_th_e(1,i)>=2*54-0.1857*H_th_h(1,i)))];
C=[C,implies(z_chp(1,i),(P_th_e(1,i)<=2*135-0.25*H_th_h(1,i))&(P_th_e(1,i)>=1.64*H_th_h(1,i)-2*50.56))];
% 火电热电耦合特性约束
C=[C,(2*43.37<=P_th_e(1,i))&(P_th_e(1,i)<=2*135)];
% 火电电功率约束
C=[C,(0<=P_chp(1,i))&(P_chp(1,i)<=2*50)];
% CHP机组总功率约束
C=[C,(0<=P_eb_e(1,i))&(P_eb_e(1,i)<=Pebmax)];
% 电锅炉出力约束
C=[C,P_loss_e(1,i)>=0];
% 失负荷惩罚为正,发电量不够5倍惩罚
C=[C,(0<=H_tes_cha(1,i))&(H_tes_cha(1,i)<=20)];
% 储热储能功率约束
C=[C,(0<=H_tes_dis(1,i))&(H_tes_dis(1,i))<=20];
% 储热释能功率约束
C=[C,(0<=0.3*SOC_HSS_max+sum(0.98*H_tes_cha(1,1:i)+k_eb*P_eb_e(1,1:i)-H_tes_dis(1,1:i)))&(0.3*SOC_HSS_max+sum(0.98*H_tes_cha(1,1:i)+k_eb*P_eb_e(1,1:i)-H_tes_dis(1,1:i))<=SOC_HSS_max)];
% 储热容量约束
C=[C,(P_th_e(1,i)+K_chp_P*P_chp(1,i)+P_wt_e(1,i)+P_pv_e(1,i)+P_loss_e(1,i)-P_eb_e(1,i)-P_p2g_e(1,i)-P_coap_e(1,i)+P_ess_e(1,i)-(5+0.269*m_ccs_st(1,i)))==P_load_e(i)];
% 电功率平衡约束
C=[C,implies(cr(i)==0,H_tes_dis(1,i)==0)];
C=[C,implies(cr(i)==1,(H_tes_cha(1,i)==0))];
% 储热不能同时充放热
C=[C,H_th_h(1,i)+0.98*H_tes_dis(1,i)-H_tes_cha(1,i)+K_chp_H*P_chp(1,i)+k_gb*P_gb_g(1,i)==H_load_h(i)];
% 热功率平衡约束
C=[C,(0<=P_gb_g(1,i))&(P_gb_g(1,i)<=Pgbmax)];
% 燃气锅炉功率约束
C=[C,(0<=m_ccs_st(1,i))&(m_ccs_st(1,i)<=120)];
% CCS再生塔最大碳捕集量限制
C=[C,(0.1*0.9*2.57*((0.4104*P_th_e(1,i)+0.06618*H_th_h(1,i)+7.092))<=m_ccs_ab(1,i))&(m_ccs_ab(1,i)<=0.9*2.57*((0.4104*P_th_e(1,i)+0.06618*H_th_h(1,i)+7.092)))];
% CCS吸收塔吸收量限制
C=[C,(0<=(v_ccs0+sum(k_ccs_co2*m_ccs_ab(1,1:i)-k_ccs_co2*m_ccs_st(1,1:i))))&(((v_ccs0+sum(k_ccs_co2*m_ccs_ab(1,1:i)-k_ccs_co2*m_ccs_st(1,1:i))))<=12000)];
% ccs富液罐体积约束
C=[C,(0<=(v_ccs0+sum(-k_ccs_co2*m_ccs_ab(1,1:i)+k_ccs_co2*m_ccs_st(1,1:i))))&(((v_ccs0+sum(-k_ccs_co2*m_ccs_ab(1,1:i)+k_ccs_co2*m_ccs_st(1,1:i))))<=12000)];
% ccs贫液罐体积约束
C=[C,(0<=P_p2g_e(1,i))&(P_p2g_e(1,i)<=30)];
% P2G功率约束
C=[C,(0<=P_chp_H2(1,i))&(P_chp_H2(1,i)<=P_p2g_e(1,i))];
% 掺氢能量约束
C=[C,(P_chp(1,i)+P_gb_g(1,i))*0.003539>=5*P_chp_H2(1,i)*0.01083];
% CHP机组掺氢比约束
C=[C,level_cool(1,i)+level_cool(2,i)+level_cool(3,i)==1];
C=[C,level_cool(3,i)==1];
% 冷却级别选择
C=[C,implies(level_cool(1,i)==1,P_coap_e(1,i)==0)];
C=[C,implies(level_cool(2,i)==1,P_coap_e(1,i)==0)];
C=[C,implies(level_cool(3,i)==1,P_coap_e(1,i)==0)];
% COAP能耗
C=[C,(0<=coap_mt(1,i))&(coap_mt(1,i)<=20)];
% 活性炭处理能力约束
C=[C,(-20<=P_ess_e(1,i))&(P_ess_e(1,i)<=20)];
% 电化学储能功率约束
C=[C,(-20<=P_ess_cha(1,i))&(P_ess_cha(1,i)<=0)];
C=[C,(0<=P_ess_dis(1,i))&(P_ess_dis(1,i)<=20)];
C=[C,implies(P_ess_e(1,i)>=0,(P_ess_dis(1,i)==P_ess_e(1,i))&(P_ess_cha(1,i)==0))];
C=[C,implies(P_ess_e(1,i)<=0,(P_ess_dis(1,i)==0)&(P_ess_cha(1,i)==P_ess_e(1,i)))];
C=[C,(0<=0.3*SOC_ESS_max+sum(P_ess_dis(1,1:i))/0.95+0.95*sum(P_ess_cha(1,1:i)))&(0.3*SOC_ESS_max+sum(P_ess_dis(1,1:i))/0.95+0.95*sum(P_ess_cha(1,1:i))<=SOC_ESS_max)];
% 电化学储能容量约束
C=[C,implies(level_cool(1,i)==1,coap_x1(1,i)==(P_th_e(1,i)+H_th_h(1,i))*6.48-0)];
C=[C,implies(level_cool(2,i)==1,coap_x1(1,i)==(P_th_e(1,i)+H_th_h(1,i))*6.48-0)];
C=[C,implies(level_cool(3,i)==1,coap_x1(1,i)==(P_th_e(1,i)+H_th_h(1,i))*6.48-0)];
C=[C,implies(level_cool(1,i)==1,coap_x1(2,i)==(P_th_e(1,i)+H_th_h(1,i))*2.88-0)];
C=[C,implies(level_cool(2,i)==1,coap_x1(2,i)==(P_th_e(1,i)+H_th_h(1,i))*2.88-0)];
C=[C,implies(level_cool(3,i)==1,coap_x1(2,i)==(P_th_e(1,i)+H_th_h(1,i))*2.88-0)];
%新能源出力约束
C=[C,(0<=P_pv_e(1,i))&(P_pv_e(1,i)<=P_PV_max(i))];
C=[C,(0<=P_wt_e(1,i))&(P_wt_e(1,i)<=P_WT_max(i))];
end
for i=2:n
C=[C,(2*-28.6<=P_chp(1,i)-P_chp(1,i-1))&(P_chp(1,i)-P_chp(1,i-1)<=2*28.6)];
% 燃气轮机爬坡约束
end
C=[C,(-1<=sum(0.98*H_tes_cha(1,1:24)+k_eb*P_eb_e(1,1:24)-H_tes_dis(1,1:24))<=1)];
C=[C,(-1<=sum(P_ess_dis(1,1:24))/0.95+0.95*sum(P_ess_cha(1,1:24))<=1)];
% C=[C,(sum(P_es(1,1:24))==0)];
%% 目标函数
F=0;
I_sale0=0;
C00=0;
punish0=0;
C_cs=0;
C_Ct0=0;
mc0=0;
Ec0=0;
CB0=0;
c_ccs0=0;
for k=1:nt
for i=1:n
C_th1=p_coal*(c(3)+c(2)*P_th_e(1,i)+c(5)*H_th_h(1,i)+c(1)*P_th_e(1,i)^2+c(6)*P_th_e(1,i)*H_th_h(1,i)+c(4)*H_th_h(1,i)^2);
% 火电机组燃料成本
I_ch4=p_gas*(((P_p2g_e(1,i)*0.85-P_chp_H2(1,i))*0.75)/0.01083);
% P2G产气收益
C_chp1=p_gas*(P_chp(1,i)/0.01083);
C_gb=p_gas*(P_gb_g(1,i)/0.01083);
% CHP机组燃料成本
m_so2=level_cool(1,i)*coap_mt(1,i)*22.8761+level_cool(2,i)*coap_mt(1,i)*55.7041+level_cool(3,i)*coap_mt(1,i)*100.008;
m_nox=level_cool(1,i)*coap_mt(1,i)*32.333+level_cool(2,i)*coap_mt(1,i)*116.349+level_cool(3,i)*coap_mt(1,i)*169.142;
% 污染气体吸收量
QQ=P_th_e(1,i)+H_th_h(1,i);
C_cb=6.4*(coap_x1(1,i))+8*(coap_x1(2,i))+1+2+0.16;
CB(i)=C_cb;
CB0=CB0+pi(k)*C_cb;
% 环保代价
mc=2.57*((C_th1)/p_coal)+0.054*0.98*3.67*(P_chp(1,i)+P_gb_g(1,i)-P_chp_H2(1,i));
% 碳排放
Ec=0.69135*(P_th_e(1,i)+K_chp_P*P_chp(1,i))+0.3*(H_th_h(1,i)+K_chp_H*P_chp(1,i)+k_gb*P_gb_g(1,i));
% 碳配额
C_ct=205*(mc-Ec-m_ccs_ab(1,i));
% 碳交易
C_buy=5*a(i)*P_loss_e(1,i);
H_eb=k_eb*P_eb_e(1,i);
H_gb=k_gb*P_gb_g(1,i);
C_cs=C_cs+pi(k)*(C_th1+C_gb+C_chp1-p_gas*P_chp_H2(1,i)/(0.01083)-I_ch4);
%燃料成本
Pccs=5+0.269*m_ccs_st(1,i);
% 碳捕集能耗
c_ccs=50*(m_ccs_st(1,i)-(0.054*0.98*3.67*((P_p2g_e(1,i)*0.85-P_chp_H2(1,i))*0.75)));
c_ccs0=c_ccs0+pi(k)*c_ccs;
% 碳封存成本
I_sale=a(i)*(P_th_e(1,i)+K_chp_P*P_chp(1,i)+P_wt_e(1,i)-P_p2g_e(1,i)+P_pv_e(1,i)-P_eb_e(1,i)-P_coap_e(1,i)-Pccs+P_ess_e(1,i))+b*(H_th_h(1,i)+0.95*H_tes_dis(1,i)-H_tes_cha(1,i)+K_chp_H*P_chp(1,i)+H_gb);
%售热售电收益
C0=C_th1+C_gb+C_chp1-p_gas*P_chp_H2(1,i)/(0.01083)+C_cb+C_buy+C_ct+c_ccs-I_ch4;
% 成本
I_sale0=I_sale0+pi(k)*I_sale;
C00=C00+pi(k)*C0;
C_Ct0=C_Ct0+pi(k)*C_ct;
mc0=mc0+pi(k)*mc;
Ec0=Ec0+pi(k)*Ec;
F=F+pi(k)*(I_sale-C0);
end
end
%%其他成本/固定收入
C_L=(4589000*((135*2)*(1+0.07+0.12))+100)/(365*60);
%火力发电建设成本
C_com=(0.1536/(1-0.1536))*(C_L+C_cs);
%火力发电运维成本
C_pv=(6000*1000*200)/(365*25)+12*sum(P_PV_max(1:24));
%光伏成本(建设加运维)
C_w=(8000*1000*440)/(365*20)+20*sum(P_WT_max(1:24));
%风电成本(建设加运维)
C_eb=(3200*1000*Pebmax)/(365*20)+2*sum(P_eb_e(1:24));
%电加热成本
C_yan=(300*1000*SOC_HSS_max)/(365*5)+1.5*sum(sum(H_tes_cha));
%储热罐成本
C_es=(2000*1000*SOC_ESS_max)/(365*5);
%电化学储能成本
I_c=(330*1000*(370))/365;
%煤电容量电价补偿
C_q=C_L+C_com+C_pv+C_w+C_eb+C_yan+C_es-I_c;
C_zhejiu=C_L+(6000*1000*200)/(365*25)+(8000*1000*440)/(365*20)+(3200*1000*Pebmax)/(365*20)+(300*1000*SOC_HSS_max)/(365*5)+(2000*1000*SOC_ESS_max)/(365*5);
C_yunwei=C_com+12*sum(P_PV_max(1:24))+20*sum(P_WT_max(1:24))+2*sum(P_eb_e(1:24))+1.5*sum(sum(H_tes_cha));
P_v_a=P_PV_max-P_pv_e;
P_w_a=P_WT_max-P_wt_e;
F_a=257.396*sum(P_w_a+P_v_a);
%% 计算
F=F-C_q-F_a;
F=-F;
ops = sdpsettings('solver','cplex', 'verbose', 2);
options.cplex.MaxTime =120;
ops.cplex.timelimit = 30;
aa=optimize(C,F,ops);
%% 结果出具、绘图
value(-F);
P_th_e=double(P_th_e);
H_th_h=double(H_th_h);
H_tes_dis=double(H_tes_dis);
H_tes_cha=double(H_tes_cha);
P_eb_e=double(P_eb_e);
P_chp=double(P_chp);
P_loss_e=double(P_loss_e);
P_gb_g=double(P_gb_g);
H_gb=k_gb*P_gb_g;
C_Ct0=double(C_Ct0);
mc0=double(mc0);
Ec0=double(Ec0);
m_ccs_st=double(m_ccs_st);
m_ccs_ab=double(m_ccs_ab);
Pccs=5+0.269*m_ccs_st;
P_p2g_e=double(P_p2g_e);
P_chp_H2=double(P_chp_H2);
P_coap_e=double(P_coap_e);
level_cool=double(level_cool);
coap_mt=double(coap_mt);
P_ess_e=double(P_ess_e);
CB=double(CB);
coap_x1=double(coap_x1);
P_pv_e=double(P_pv_e);
P_wt_e=double(P_wt_e);
c_ccs0=double(c_ccs0);
P_ess_cha=double(P_ess_cha);
P_ess_dis=double(P_ess_dis);
M_CO2_P2G=(0.54*0.98*3.67*((P_p2g_e*0.85-P_chp_H2)*0.75));
M_CCS_CO2=m_ccs_st;
for i=2:25
E_SOC_elec(i)=0.3*SOC_ESS_max+sum(P_ess_dis(1,1:i-1))/0.95+0.95*sum(P_ess_cha(1,1:i-1));
H_SOC_heat(i)=0.3*SOC_HSS_max+sum(0.98*H_tes_cha(1,1:i-1)+k_eb*P_eb_e(1,1:i-1)-H_tes_dis(1,1:i-1));
% m_coap_so2(i-1)=cool(1,i-1)*mt(1,i-1)*22.8761+cool(2,i-1)*mt(1,i-1)*55.7041+cool(3,i-1)*mt(1,i-1)*100.008;
% m_coap_nox(i-1)=cool(1,i-1)*mt(1,i-1)*32.333+cool(2,i-1)*mt(1,i-1)*116.349+cool(3,i-1)*mt(1,i-1)*169.142;
m_th_so2(i-1)=(P_th_e(1,i-1)+H_th_h(1,i-1))*6.48;
m_th_nox(i-1)=(P_th_e(1,i-1)+H_th_h(1,i-1))*2.88;
m_coap_so2(i-1)=-coap_x1(1,i-1)+m_th_so2(i-1);
m_coap_nox(i-1)=-coap_x1(2,i-1)+m_th_nox(i-1);
end
%% 结果可视化
% 新能源出力
figure(1);
set(groot, 'DefaultAxesFontName', 'Times New Roman');
set(groot, 'DefaultTextFontName', 'Times New Roman');
set(groot, 'DefaultLineLineWidth', 1.5);
x = 1:24;
for i = 1:24
if P_PV_max(i)-P_pv_e(i)>0
PC1(i)=patch([x(i)-0.4, x(i)+0.4, x(i)+0.4, x(i)-0.4], [P_pv_e(i), P_pv_e(i), P_PV_max(i), P_PV_max(i)], 'b');
set(PC1(i), 'FaceColor', [0.97, 0.75, 0.40] , 'EdgeColor', 'none');
bj=i;
end
if P_WT_max(i)-P_wt_e(i)>0
PC4(i)=patch([x(i)-0.4, x(i)+0.4, x(i)+0.4, x(i)-0.4], [P_wt_e(i), P_wt_e(i), P_WT_max(i), P_WT_max(i)], 'b');
set(PC4(i), 'FaceColor', [0.40, 0.70, 0.95] , 'EdgeColor', 'none');
bj1=i;
end
end
hold on
PC2=plot(x,P_pv_e,'-rv');
set(PC2, 'Color', [0.90, 0.55, 0.20] ,'LineWidth',1);
hold on
PC3=plot(x,P_PV_max,'-ro');
set(PC3, 'Color', [0.80, 0.35, 0.10] ,'LineWidth',1);
hold on
PC5=plot(x,P_wt_e,'-rv');
set(PC5, 'Color', [0.20, 0.50, 0.80],'LineWidth',1);
hold on
PC6=plot(x,P_WT_max,'-ro');
set(PC6, 'Color', [0.10, 0.30, 0.60] ,'LineWidth',1);
grid on;
legend([PC1(bj),PC2,PC3,PC4(bj1),PC5,PC6], 'Curtailed PV Power','Actual PV Output','Maximum PV Output','Curtailed WT Power','Actual WT Output','Maximum WT Output');
xlabel('time (h)');
ylabel('Power (MW)');
title('Output Curve of New Energy Units');
ylim([-1 240]);
% 功率平衡图绘制
figure(2);
% 全局样式设置(SCI期刊标配)
set(groot, 'DefaultAxesFontName', 'Times New Roman');
set(groot, 'DefaultTextFontName', 'Times New Roman');
set(groot, 'DefaultLineLineWidth', 1.5);
% set(gcf, 'Color', 'white'); % 背景设为白色
set(gcf, 'Position', [10 10 1500 350]); % 调整图窗尺寸适配排版
% 子图1:电负荷平衡图
subplot(1,3,1);
P_v_a = P_PV_max - P_pv_e;
P_w_a = P_WT_max - P_wt_e;
P = [P_chp*K_chp_P; P_th_e; P_pv_e; P_wt_e; -P_eb_e; -Pccs; -P_p2g_e; -P_coap_e; P_ess_cha; P_ess_dis]';
e = bar(P, 'stacked');
hold on;
plot(1:24, P_load_e, '-b*', 'LineWidth', 1, 'MarkerSize', 6); % 优化标记尺寸
grid on;
% 电负荷平衡图配色(SCI专业色系,保留原色系逻辑并优化对比度)
colors_e = [
0.7 0.2 0.7 % CHP机组 - 紫色
0.9 0.4 0.2 % 火电机组 - 橙色
1 0.8 0.2 % 光伏机组 - 黄色
0.5 0.8 0.2 % 风电机组 - 绿色
0.6 0.6 0.6 % 电锅炉 - 灰色
0.2 0.6 0.8 % CCS能耗 - 蓝色
0.8 0.3 0.3 % P2G能耗 - 红色
0.1 0.1 0.6 % COAP能耗 - 深蓝
0.2 0.8 0.8 % 储能充电 - 青色
0.3 0.7 0.9 % 储能放电 - 蓝绿色
];
for i = 1:size(colors_e,1)
e(i).FaceColor = colors_e(i,:);
e(i).EdgeColor = 'none'; % 去除边框更素雅
end
% 图例与标签(英文显示,优化图例位置)
legend('CHP','CTP','PV','WT','EB','CCS','P2G','COAP','ESS Charging','ESS Discharging','Electric Load', 'Location', 'best', 'FontSize', 9);
xlabel('Time Interval');
ylabel('Power (MW)');
axis([0,25,-200,650]);
set(gca, 'FontSize', 10);
set(gca, 'GridLineStyle', ':', 'GridColor', [0.7,0.7,0.7]); % 优化网格样式
xticks(0:2:24); % X轴刻度优化,避免拥挤
title('Electric Power Balance Diagram');
% 子图2:热负荷平衡图
subplot(1,3,2);
H1 = [K_chp_H*P_chp; H_th_h; H_gb; 0.98*H_tes_dis; -(H_tes_cha)]';
h = bar(H1, 'stacked');
hold on;
plot(1:24, H_load_h, '-r*', 'LineWidth', 1, 'MarkerSize', 6);
grid on;
% 热负荷平衡图配色
colors_h = [
0.7 0.2 0.7 % CHP供热 - 紫色
0.9 0.4 0.2 % 火电供热 - 橙色
0.3 0.6 1 % 燃气锅炉 - 蓝色
0.5 0.8 0.2 % 储热罐放热 - 绿色
0.8 0.3 0.3 % 储热罐充热 - 红色
];
for i = 1:size(colors_h,1)
h(i).FaceColor = colors_h(i,:);
h(i).EdgeColor = 'none';
end
legend('CHP','CTP','GB','TES Discharge','TES Charging','Thermal Load', 'Location', 'best');
xlabel('Time Interval');
ylabel('Power (MW)');
axis([0,25,-50,290]);
set(gca, 'FontSize', 10);
set(gca, 'GridLineStyle', ':', 'GridColor', [0.7,0.7,0.7]);
xticks(0:2:24);
title('Thermal Power Balance Diagram');
% 子图3:气平衡图
subplot(1,3,3);
V_chp1 = (P_chp(1,:));
V_gb = (P_gb_g);
V_ch4 = (((P_p2g_e*0.85-P_chp_H2)*0.75));
V_buy = V_chp1 + V1 + V_gb - V_ch4;
V = [-V_chp1; -V_gb; V_ch4; V_buy]';
g = bar(V, 'stacked');
hold on;
plot(1:24, V1, '-m*', 'MarkerSize', 6);
grid on;
% 气平衡图配色
colors_g = [
0.7 0.2 0.7 % CHP - 紫色
0.3 0.6 1 % 燃气锅炉 - 蓝色
0.8 0.3 0.3 % P2G - 红色
0.5 0.8 0.2 % 气负荷 - 绿色
];
for i = 1:size(colors_g,1)
g(i).FaceColor = colors_g(i,:);
g(i).EdgeColor = 'none';
end
legend('CHP','GB','P2G','Gas Purchase','Gas Load', 'Location', 'best', 'FontSize', 9);
xlabel('Time Interval');
ylabel('Gas Volume (MW equivalent)'); % 补充气维度标注,更规范
axis([0,25,-150,270]);
set(gca);
set(gca, 'GridLineStyle', ':', 'GridColor', [0.7,0.7,0.7]);
xticks(0:2:24);
title('Gas Power Balance Diagram');
% 储能状态变化
figure(3);
% 全局样式设置(SCI期刊标配)
set(groot, 'DefaultAxesFontName', 'Times New Roman');
set(groot, 'DefaultTextFontName', 'Times New Roman');
set(groot, 'DefaultLineLineWidth', 1.5);
% set(gcf, 'Color', 'white'); % 背景设为白色
% 子图1:电储能状态变化
subplot(2,1,1);
% 绘制电储能曲线(SCI专业配色,优化标记尺寸)
plot(0:24, E_SOC_elec, '*-', 'Color', [0.10, 0.30, 0.60], 'MarkerSize', 6);
xlabel('Time (h)');
ylabel('Stored Electric Energy (MWh)');
title('Electric Energy Storage State Variation');
grid on;
% 优化网格样式(浅灰色虚线,不喧宾夺主)
set(gca, 'GridLineStyle', ':', 'GridColor', [0.7, 0.7, 0.7]);
set(gca, 'LineWidth', 1); % 坐标轴线条宽度
ylim([0, SOC_ESS_max]);
set(gca, 'FontSize', 10);
xticks(0:2:24); % X轴刻度优化,避免拥挤
% 子图2:热储能状态变化
subplot(2,1,2);
% 绘制热储能曲线(SCI专业配色,优化标记尺寸)
plot(0:24, H_SOC_heat, '*-', 'Color', [0.80, 0.35, 0.10], 'MarkerSize', 6);
xlabel('Time (h)');
ylabel('Stored Thermal Energy (MWh)');
title('Thermal Energy Storage State Variation');
grid on;
set(gca, 'GridLineStyle', ':', 'GridColor', [0.7, 0.7, 0.7]);
set(gca, 'LineWidth', 1);
ylim([0, SOC_HSS_max]);
set(gca, 'FontSize', 10);
xticks(0:2:24);
% 火电机组运行轨迹
figure(4);
set(groot, 'DefaultAxesFontName', 'Times New Roman');
set(groot, 'DefaultTextFontName', 'Times New Roman');
set(groot, 'DefaultLineLineWidth', 1.5);
% set(gcf, 'Color', 'white'); % 背景设为白色
% 定义四边形可行域顶点坐标
x_vertex = [0, 114.54, 196.36, 0];
y_vertex = [108, 86.73, 220.91, 270];
% 填充可行域(SCI素雅配色)
fill(x_vertex, y_vertex,[0.8 0.9 0.8], 'EdgeColor', 'none');
hold on;
% 绘制运行点轨迹(SCI专业配色,优化标记尺寸)
plot(H_th_h, P_th_e, '*-', 'Color', [0.80, 0.35, 0.10], 'MarkerSize', 6);
xlabel('Heat Output (MW)');
ylabel('Power Output (MW)');
title('Thermal Power Unit Operation Trajectory');
legend('Feasible Region', 'Operation Points', 'Location', 'best');
grid on;
set(gca, 'GridLineStyle', ':', 'GridColor', [0.7, 0.7, 0.7]);
set(gca, 'LineWidth', 1);
axis([-10, 210, 70, 280]);
% P2G-CCS系统分析
figure(5);
set(groot, 'DefaultAxesFontName', 'Times New Roman'); % 统一字体为SCI标配
set(groot, 'DefaultTextFontName', 'Times New Roman');
set(groot, 'DefaultLineLineWidth', 1.5); % 线条宽度适配SCI印刷
% set(gcf, 'Color', 'white'); % 背景设为白色
subplot(2,1,1);
yyaxis left;
plot(1:24, P_p2g_e, 'b-*', 'MarkerSize', 6); % 绘制电功率输入曲线
ylabel('Electric Power Input (MW)');
ylim([-1 50]);
yyaxis right;
plot(1:24, V_ch4, 'r-o', 'MarkerSize', 6); % 绘制天然气产量曲线
ylabel('Natural Gas Production (10^4 m^3)');
ylim([-1 50]);
title('P2G System Performance');
xlabel('Time (h)');
legend('Power Consumption', 'Gas Production', 'Location', 'best');
grid on;
set(gca, 'GridLineStyle', ':', 'GridColor', [0.7,0.7,0.7]); % 网格样式优化
set(gca, 'LineWidth', 1);
xticks(0:2:24); % X轴刻度优化
ax1 = gca;
% 左侧Y轴
ax1.YAxis(1).Color = 'k';
ax1.YAxis(1).TickLabelColor = 'k';
ax1.YAxis(1).Label.Color = 'k';
% 右侧Y轴
ax1.YAxis(2).Color = 'k';
ax1.YAxis(2).TickLabelColor = 'k';
ax1.YAxis(2).Label.Color = 'k';
% 子图2:碳源利用分析
subplot(2,1,2);
carbon_flow = [M_CCS_CO2;M_CO2_P2G];
bar(carbon_flow', 1.2, 'FaceColor', 'flat'); % 绘制碳流量柱状图
colormap([0.10,0.30,0.60; 0.80,0.35,0.10]);
title('Carbon Source Utilization Analysis');
xlabel('Time (h)');
ylabel('CO_2 (ton)');
legend('CCS Capture Amount', 'P2G Consumption Amount', 'Location', 'best');
ylim([-1 170]);
grid on;
set(gca, 'GridLineStyle', ':', 'GridColor', [0.7,0.7,0.7]);
set(gca, 'LineWidth', 1);
xticks(0:2:24);
% COAP处理污染物系统分析
figure(6)
set(groot, 'DefaultAxesFontName', 'Times New Roman'); % 统一字体为SCI标配
set(groot, 'DefaultTextFontName', 'Times New Roman');
set(groot, 'DefaultLineLineWidth', 1.5); % 线条宽度适配SCI印刷
% set(gcf, 'Color', 'white'); % 背景设为白色
level_cool=zeros(24,1);
color_so2_gen = [0.80, 0.35, 0.10];
color_so2_ads = [0.90, 0.55, 0.20];
color_nox_gen = [0.10, 0.30, 0.60];
color_nox_ads = [0.20, 0.50, 0.80];
color_level = [0.85, 0.85, 0.85];
color_power = [0.10, 0.50, 0.20];
subplot(2,1,1);
yyaxis left;
plot(1:24, m_th_so2, '-', 'Color', color_so2_gen, 'DisplayName', 'SO₂ Generation Amount');
hold on;
plot(1:24, m_coap_so2, '--', 'Color', color_so2_ads, 'DisplayName', 'SO₂ Adsorption Capacity');
ylabel('SO₂ mass (t)');
ylim([0, max([m_th_so2, m_coap_so2,m_th_nox, m_coap_nox])*1.2]);
yyaxis right;
plot(1:24, m_th_nox, '-', 'Color', color_nox_gen, 'DisplayName', 'NOₓ Generation Amount');
plot(1:24, m_coap_nox, '--', 'Color', color_nox_ads, 'DisplayName', 'NOₓ Adsorption Capacity');
ylabel('NOₓ mass (t)');
ylim([0, max([m_th_so2, m_coap_so2,m_th_nox, m_coap_nox])*1.2]);
title('Pollutant Gas Treatment');
xlabel('Time (h)');
legend('Location', 'best');
grid on;
set(gca, 'GridLineStyle', ':', 'GridColor', [0.7,0.7,0.7]);
set(gca, 'LineWidth', 1);
xticks(0:2:24);
% 设置子图1纵坐标刻度和字体为黑色
ax1 = gca;
ax1.YAxis(1).Color = 'k';
ax1.YAxis(1).TickLabelColor = 'k';
ax1.YAxis(1).Label.Color = 'k';
ax1.YAxis(2).Color = 'k';
ax1.YAxis(2).TickLabelColor = 'k';
ax1.YAxis(2).Label.Color = 'k';
subplot(2,1,2);
yyaxis left;
bar(1:24, level_cool, 0.4, 'FaceColor', color_level, 'EdgeColor', 'none', 'DisplayName', 'Operation Level');
ylabel('Cooling Stage');
ylim([0, 5]);
yyaxis right;
plot(1:24, P_coap_e, '-*', 'Color', color_power, 'MarkerSize', 6, 'DisplayName', 'COAP Power Consumption');
ylabel('Power (MW)');
ylim([0, max(P_coap_e+100)*1.3]);
title('COAP Operation Status');
xlabel('Time (h)');
legend('Location', 'best');
grid on;
set(gca, 'GridLineStyle', ':', 'GridColor', [0.7,0.7,0.7]);
set(gca, 'LineWidth', 1);
xticks(0:2:24);
% 设置子图2纵坐标刻度和字体为黑色
ax2 = gca;
ax2.YAxis(1).Color = 'k';
ax2.YAxis(1).TickLabelColor = 'k';
ax2.YAxis(1).Label.Color = 'k';
ax2.YAxis(2).Color = 'k';
ax2.YAxis(2).TickLabelColor = 'k';
ax2.YAxis(2).Label.Color = 'k';
C00=double(C00);
F1=[double(I_sale0)-double(punish0)+sum(V1*p_gas/0.01083),I_c,-double(C_cs)-sum(V1*p_gas/0.01083),-double(C_zhejiu),-double(C_yunwei)-c_ccs0,-C_Ct0,-sum(double(CB0)),-double(F_a),-5*sum(P_loss_e.*a)];
sum(F1)
disp([' 总收益: ', num2str(sum(F1)/(1000*7.06), '%.2f'), ' 千美元']);
disp('========== 运行收益分解 ==========');
disp([' 售能收益: ', num2str(F1(1)/(1000*7.06), '%.2f'), ' 千美元']);
disp([' 容量电价补偿收益: ', num2str(F1(2)/(1000*7.06), '%.2f'), ' 千美元']);
disp([' 燃料成本: ', num2str(F1(3)/(1000*7.06), '%.2f'), ' 千美元']);
disp([' 设备折旧成本: ', num2str(F1(4)/(1000*7.06), '%.2f'), ' 千美元']);
disp([' 运维成本: ', num2str(F1(5)/(1000*7.06), '%.2f'), ' 千美元']);
disp([' 碳交易收益/成本: ', num2str(F1(6)/(1000*7.06), '%.2f'), ' 千美元']);
disp([' 环保代价处理成本: ', num2str(F1(7)/(1000*7.06), '%.2f'), ' 千美元']);
disp([' 弃电惩罚成本: ', num2str(F1(8)/(1000*7.06), '%.2f'), ' 千美元']);
disp([' 失负荷惩罚成本: ', num2str(F1(9)/(1000*7.06), '%.2f'), ' 千美元']);
disp([' 弃风率: ',num2str(100*(1-sum(P_wt_e)/sum(P_WT_max)), '%.2f'),'%']);
disp([' 弃光率: ',num2str(100*(1-sum(P_pv_e)/sum(P_PV_max)), '%.2f'),'%']);
disp([' 失负荷电量: ',num2str(sum(P_loss_e), '%.2f'),'MW']);
% 成本与收益可视化
figure(7);
conv = 1/(1000*7.06); % 转换系数:千美元
threshold = 1e-3; % 近0阈值,绝对值小于该值的项目过滤
item_names = {
'Energy Sales Revenue'; % 1:售能收益
'Capacity Tariff Revenue'; % 2:容量电价补偿收益
'Fuel Cost'; % 3:燃料成本
'Equipment Depreciation Cost'; % 4:设备折旧成本
'O&M Cost'; % 5:运维成本
'Carbon Trading Revenue/Cost'; % 6:碳交易收益/成本
'Environmental Treatment Cost'; % 7:环保代价处理成本
'Curtailed Power Penalty Cost'; % 8:弃电惩罚成本
'Load Loss Penalty Cost' % 9:失负荷惩罚成本
};
item_values = F1 * conv; % 转换为千美元
% 收益项目(正数值,过滤近0)
revenue_mask = (item_values > 0) & (abs(item_values) > threshold);
revenue_names = item_names(revenue_mask);
revenue_values = item_values(revenue_mask);
% 成本项目(负数值,过滤近0,取绝对值用于饼图)
cost_mask = (item_values < 0) & (abs(item_values) > threshold);
cost_names = item_names(cost_mask);
cost_values = abs(item_values(cost_mask));
set(groot, 'DefaultAxesFontName', 'Times New Roman');
set(groot, 'DefaultTextFontName', 'Times New Roman');
set(groot, 'DefaultLineLineWidth', 1.5);
pie_colors = [0.80,0.35,0.10; 0.10,0.30,0.60; 0.10,0.50,0.20; 0.70,0.20,0.70; 0.60,0.60,0.60; 0.20,0.60,0.80; 0.30,0.70,0.90];
bar_colors = [0.10,0.30,0.60; 0.80,0.35,0.10; 0.10,0.50,0.20];
subplot(1,2,1);
explode = zeros(1, length(cost_values));
min_sector_idx = find(cost_values / sum(cost_values) < 0.02);
if ~isempty(min_sector_idx)
explode(min_sector_idx) = 0.5; % 极小扇区爆炸
else
explode(end) = 1;
end
h_pie = pie(cost_values, explode);
pie_patches = findobj(h_pie, 'Type', 'patch'); % 获取饼图扇区对象
for i = 1:length(pie_patches)
pie_patches(i).FaceColor = pie_colors(i,:);
end
legend(cost_names, 'Location', 'best', 'FontSize', 8);
title('Cost Breakdown (Thousand USD)', 'FontSize', 12);
set(gca, 'FontSize', 10);
t1 = title('Cost Breakdown (Thousand USD)', 'FontSize', 12);
pos = t1.Position;
t1.Position = [pos(1), pos(2) + 0.15, pos(3)]; % 增大Y坐标实现标题上移
subplot(1,2,2);
h = barh(revenue_values);
set(gca, 'YTick', 1:length(revenue_names));
set(gca, 'YTickLabel', revenue_names);
title('Revenue Breakdown (Thousand USD)', 'FontSize', 12);
xlabel('Revenue (Thousand USD)', 'FontSize', 10);
grid on;
set(gca, 'FontSize', 10);
for i = 1:length(revenue_values)
text(revenue_values(i), i, sprintf('%.2f', revenue_values(i)), ...
'HorizontalAlignment', 'left', 'VerticalAlignment', 'middle', ...
'FontSize', 9, 'FontName', 'Times New Roman');
end
set(gcf, 'Position', [100 100 1200 500]); % 调整图窗大小
set(gca, 'FontSize', 10);
最新发布