检查和整理以下代码:%% 场景3:综合能源系统优化调度(含碳交易+CCS+P2G)
clc; clear;
%% ================== 基础数据定义 ==================
% 电负荷 (MW)
P_Load_e = [217.10, 247.90, 267.08, 203.38, 193.86, 227.46, 274.50, 310.90, 330.36, 345.48, 400.78, 346.60, 319.30, 305.58, 359.76, 366.34, 346.04, 367.04, 417.72, 433.40, 407.36, 284.30, 320.28, 333.86];
% 热负荷 (MW)
H_Load_h = [183.41, 206.28, 209.52, 219.91, 226.59, 213.21, 226.84, 203.28, 169.32, 152.46, 146.45, 172.79, 183.65, 199.35, 174.41, 160.08, 152.46, 147.38, 129.13, 142.30, 134.67, 155.00, 142.30, 127.05];
% 气象数据
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,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,7.10];
% 能源价格
C_p_coal = 780; % 标煤价格 (元/吨)
C_p_gas = 3.2; % 燃气价格 (元/m³)
% 分时电价 (元/kWh -> 元/MWh)
C_buy_e = [0.17*ones(1,8), 0.49*ones(1,10), 0.83*ones(1,6)] * 1000;
% 碳交易参数
C_p_carbon = 230; % 碳交易价格 (元/吨CO2)
C_ccs_fixed = 50; % CCS固定运行成本 (元/吨CO2)
% P2G系统参数
P_P2G_max = 80; % P2G最大电功率输入 (MW)
k_eta_P2G = 0.65; % P2G电转气效率
k_co2_P2G = 1.96; % P2G单位产气量CO2消耗 (kg/m³)
C_fixed_p2g = 1.2e6; % P2G年固定成本 (元/年)
%% ================== 新能源机组模型 ==================
% 光伏模型
S_PV = 40000; % 光伏面积 (m²)
k_eta_PV = 0.157; % 光电转换效率
P_PV_max = (S_PV * P_Solar * k_eta_PV) / 1000;
% 风电模型
Q_S_WT = 80000; % 额定容量 (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;
end
end
%% ================== 设备参数设置 ==================
% 储能系统参数
SOC_ESS_max = 75; % 电储能容量 (MWh)
SOC_HSS_max = 50; % 储热罐容量 (MWh)
P_ess_max = 30; % 蓄电池最大充放电功率 (MW)
H_hss_max = 20; % 储热罐最大充放热功率 (MW)
% 机组运行参数
P_CHP_max = 130; % CHP机组最大功率
P_GB_max = 100; % 燃气锅炉最大功率
P_EB_max = 120; % 电锅炉最大功率
G_heat_v = 0.01083; % 燃气热值 (MWh/m³)
% 能量效率参数
k_eta_CHP_e = 0.42; k_eta_CHP_h = 0.29; k_eta_gb = 0.82; k_eta_eb = 0.97;
k_eta_ess_cha = 0.95; k_eta_ess_dis = 0.93; k_eta_hss_cha = 0.95; k_eta_hss_dis = 0.93;
D_self_elec = 0.0001; D_self_heat = 0.0002;
% 系统运行参数
N_Hours = 24; % 时段数
Ramp_elec_th = 20; % 火电电出力爬坡率 (MW/h)
Ramp_heat_th = 15; % 火电热出力爬坡率 (MW/h)
% 火电机组煤耗系数
k_c_xs = [0.0004818, 2.695, 52, 0.0004751, 0.3536, 0.001004];
% CCS系统参数
k_ccs_co2 = 20.23; v_ccs_int = 6000; CCS_ab_max = 120;
P_CCS_base = 2.51; P_CCS_ton = 0.269;
% 初始状态参数
P_TH_int = 50.73; H_TH_int = 30.11;
%% ================== 决策变量定义 ==================
% 基础设备变量
P_TH_e = sdpvar(1, N_Hours); H_TH_h = sdpvar(1, N_Hours);
P_CHP_e = sdpvar(1, N_Hours); H_CHP_h = sdpvar(1, N_Hours);
H_GB_h = sdpvar(1, N_Hours); H_EB_h = sdpvar(1, N_Hours);
P_g_CHP = sdpvar(1, N_Hours); P_g_GB = sdpvar(1, N_Hours);
P_e_EB = sdpvar(1, N_Hours); P_e_CCS = sdpvar(1, N_Hours);
P_buy_e = sdpvar(1, N_Hours); P_pv_e = sdpvar(1, N_Hours);
P_wt_e = sdpvar(1, N_Hours);
% 储能变量
P_ESS_cha = sdpvar(1, N_Hours); P_ESS_dis = sdpvar(1, N_Hours);
H_HSS_cha = sdpvar(1, N_Hours); H_HSS_dis = sdpvar(1, N_Hours);
E_SOC_elec = sdpvar(1, N_Hours+1); H_SOC_heat = sdpvar(1, N_Hours+1);
% CCS系统变量
M_CCS_ab = sdpvar(1, N_Hours); M_CCS_st = sdpvar(1, N_Hours);
T_v_rich = sdpvar(1, N_Hours+1); T_v_lean = sdpvar(1, N_Hours+1);
% P2G系统变量
P_e_P2G = sdpvar(1, N_Hours); V_g_P2G = sdpvar(1, N_Hours);
M_co2_P2G = sdpvar(1, N_Hours); Gas_buy = sdpvar(1, N_Hours);
% 状态变量
B_ess_s = binvar(1, N_Hours); B_hss_s = binvar(1, N_Hours);
%% ================== 约束条件 ==================
Cons = [];
% 初始状态
E_SOC_int = 0.53 * SOC_ESS_max; H_SOC_int = 0.47 * SOC_HSS_max;
T_v_rich(1) = v_ccs_int; T_v_lean(1) = v_ccs_int;
Cons = [Cons, E_SOC_elec(1) == E_SOC_int, H_SOC_heat(1) == H_SOC_int];
for i = 1:N_Hours
% CHP机组约束
Cons = [Cons, 0 <= P_g_CHP(i) <= P_CHP_max, P_CHP_e(i) == k_eta_CHP_e * P_g_CHP(i), H_CHP_h(i) == k_eta_CHP_h * P_g_CHP(i)];
% 燃气锅炉约束
Cons = [Cons, 0 <= P_g_GB(i) <= P_GB_max, H_GB_h(i) == k_eta_gb * P_g_GB(i)];
Cons = [Cons, 0 <= P_e_EB(i) <= P_EB_max, H_EB_h(i) == k_eta_eb * P_e_EB(i)];
% P2G约束
Cons = [Cons, 0 <= P_e_P2G(i) <= P_P2G_max];
Cons = [Cons, V_g_P2G(i) == k_eta_P2G * P_e_P2G(i) / G_heat_v];
Cons = [Cons, M_co2_P2G(i) == (k_co2_P2G * V_g_P2G(i)) / 1000];
Cons = [Cons, Gas_buy(i) + V_g_P2G(i) == (P_g_CHP(i) + P_g_GB(i)) / G_heat_v];
% 火电机组约束
Cons = [Cons, 50.73 <= P_TH_e(i) <= 135.54, 30.11 <= H_TH_h(i) <= 98.18];
Cons = [Cons, P_TH_e(i) >= 1.64*H_TH_h(i) - 50.56, P_TH_e(i) <= 135.54 - 0.25*H_TH_h(i)];
% 储能约束
Cons = [Cons, 0 <= P_ESS_cha(i) <= P_ess_max * (1 - B_ess_s(i)), 0 <= P_ESS_dis(i) <= P_ess_max * B_ess_s(i)];
Cons = [Cons, 0 <= H_HSS_cha(i) <= H_hss_max * (1 - B_hss_s(i)), 0 <= H_HSS_dis(i) <= H_hss_max * B_hss_s(i)];
% CCS约束
Cons = [Cons, 0 <= M_CCS_ab(i) <= CCS_ab_max, 0 <= M_CCS_st(i) <= CCS_ab_max];
Cons = [Cons, T_v_rich(i+1) == T_v_rich(i) + k_ccs_co2*M_CCS_ab(i) - k_ccs_co2*M_CCS_st(i), 0 <= T_v_rich(i+1) <= 12000];
Cons = [Cons, T_v_lean(i+1) == T_v_lean(i) - k_ccs_co2*M_CCS_ab(i) + k_ccs_co2*M_CCS_st(i), 0 <= T_v_lean(i+1) <= 12000];
Cons = [Cons, P_e_CCS(i) == P_CCS_base + P_CCS_ton * M_CCS_st(i)];
% 能量平衡
Cons = [Cons, P_TH_e(i) + P_CHP_e(i) + P_wt_e(i) + P_pv_e(i) + P_buy_e(i) + P_ESS_dis(i) == P_Load_e(i) + P_ESS_cha(i) + P_e_EB(i) + P_e_CCS(i) + P_e_P2G(i)];
Cons = [Cons, H_TH_h(i) + H_CHP_h(i) + H_GB_h(i) + H_EB_h(i) + H_HSS_dis(i) == H_Load_h(i) + H_HSS_cha(i)];
% 新能源出力限制
Cons = [Cons, 0 <= P_pv_e(i) <= P_PV_max(i), 0 <= P_wt_e(i) <= P_WT_max(i), 0 <= P_buy_e(i) <= 200];
% 储能状态更新
Cons = [Cons, E_SOC_elec(i+1) == (E_SOC_elec(i) - D_self_elec * SOC_ESS_max) + k_eta_ess_cha*P_ESS_cha(i) - P_ESS_dis(i)/k_eta_ess_dis];
Cons = [Cons, 0.17*SOC_ESS_max <= E_SOC_elec(i+1) <= 0.93*SOC_ESS_max];
Cons = [Cons, H_SOC_heat(i+1) == (H_SOC_heat(i) - D_self_heat * SOC_HSS_max) + k_eta_hss_cha*H_HSS_cha(i) - H_HSS_dis(i)/k_eta_hss_dis];
Cons = [Cons, 0.23*SOC_HSS_max <= H_SOC_heat(i+1) <= 0.95*SOC_HSS_max];
end
% 火电机组爬坡约束
Cons = [Cons, P_TH_e(1) >= P_TH_int - Ramp_elec_th, P_TH_e(1) <= P_TH_int + Ramp_elec_th];
Cons = [Cons, H_TH_h(1) >= H_TH_int - Ramp_heat_th, H_TH_h(1) <= H_TH_int + Ramp_heat_th];
for i = 2:N_Hours
Cons = [Cons, -Ramp_elec_th <= P_TH_e(i)-P_TH_e(i-1) <= Ramp_elec_th];
Cons = [Cons, -Ramp_heat_th <= H_TH_h(i)-H_TH_h(i-1) <= Ramp_heat_th];
end
% 储能结束状态约束
Cons = [Cons, 0.45*SOC_HSS_max <= H_SOC_heat(end) <= 0.55*SOC_HSS_max];
Cons = [Cons, 0.45*SOC_ESS_max <= E_SOC_elec(end) <= 0.55*SOC_ESS_max];
% CCS结束状态约束
Cons = [Cons, v_ccs_int - 1000 <= T_v_rich(end) <= v_ccs_int + 1000];
Cons = [Cons, v_ccs_int - 1000 <= T_v_lean(end) <= v_ccs_int + 1000];
%% ================== 目标函数 ==================
F1_fuel_cost = 0; F2_carbon_cost = 0; F3_penalty_cost = 0;
for i = 1:N_Hours
% 火电煤耗计算
Coal_th_kg = k_c_xs(1)*P_TH_e(i)^2 + k_c_xs(2)*P_TH_e(i) + k_c_xs(3)*H_TH_h(i)^2 + k_c_xs(4)*H_TH_h(i) + k_c_xs(5)*P_TH_e(i)*H_TH_h(i) + k_c_xs(6);
Coal_th_ton = Coal_th_kg / 1000;
% 总燃料成本计算
F1_fuel_cost = F1_fuel_cost + (C_p_coal * Coal_th_ton + C_p_gas * Gas_buy(i));
% 碳排放计算
Coal_emission = 2.67 * Coal_th_ton;
Gas_emission = (P_g_CHP(i) + P_g_GB(i)) / G_heat_v * 1.96 / 1000;
Total_emission = Coal_emission + Gas_emission;
% 碳配额计算
Carbon_quota_per = 0.57;
Carbon_quota = Carbon_quota_per * P_Load_e(i);
% 碳交易成本
F2_carbon_cost = F2_carbon_cost + C_p_carbon * (Total_emission - Carbon_quota - M_CCS_ab(i) - M_co2_P2G(i));
% 购电成本计算
F3_penalty_cost = F3_penalty_cost + (C_buy_e(i) * P_buy_e(i));
end
% 设备固定成本 (按年分摊)
days_per_year = 365;
C_fixed_th = (4.589e6 * 135.54) / days_per_year;
C_fixed_chp = (3.2e6 * 130) / days_per_year;
C_fixed_gb = (700e3 * P_GB_max) / days_per_year;
C_fixed_eb = (800e3 * P_EB_max) / days_per_year;
C_fixed_ess = (2e6 * SOC_ESS_max) / days_per_year;
%% 场景3:综合能源系统优化调度(含碳交易+CCS+P2G)
clc
clear
%% ================== 基础数据定义 ==================
% 电负荷 (MW)
P_Load_e = [217.10, 247.90, 267.08, 203.38, 193.86, 227.46, 274.50, 310.90, 330.36, 345.48, 400.78, 346.60, 319.30, 305.58, 359.76, 366.34, 346.04, 367.04, 417.72, 433.40, 407.36, 284.30, 320.28, 333.86];
% 热负荷 (MW)
H_Load_h = [183.41, 206.28, 209.52, 219.91, 226.59, 213.21, 226.84, 203.28, 169.32, 152.46, 146.45, 172.79, 183.65, 199.35, 174.41, 160.08, 152.46, 147.38, 129.13, 142.30, 134.67, 155.00, 142.30, 127.05];
% 气象数据(辐照度 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,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,7.10];
% 能源价格参数
C_p_coal = 780; % 标煤价格 (元/吨)
C_p_gas = 3.2; % 燃气价格 (元/m³)
% 分时电价 (元/kWh -> 元/MWh)
C_buy_e = [0.17*ones(1,8), 0.49*ones(1,10), 0.83*ones(1,6)] * 1000;
% 碳交易参数
C_p_carbon = 230; % 碳交易价格 (元/吨CO2)
C_ccs_fixed = 50; % CCS固定运行成本 (元/吨CO2)
%% ================== 新能源机组模型 ==================
% 光伏模型 (单位: MW)
S_PV = 40000; % 光伏面积 (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; % 额定容量 (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
%% ================== 设备参数设置 ==================
% 储能系统参数
SOC_ESS_max = 75; % 电储能容量 (MWh)
SOC_HSS_max = 50; % 储热罐容量 (MWh)
P_ess_max = 30; % 蓄电池最大充放电功率 (MW)
H_hss_max = 20; % 储热罐最大充放热功率 (MW)
% 机组运行参数
P_CHP_max = 130; % CHP机组最大功率
P_GB_max = 100; % 燃气锅炉最大功率
P_EB_max = 120; % 电锅炉最大功率
G_heat_v = 0.01083; % 燃气热值 (MWh/m³)
% 能量效率参数
k_eta_CHP_e = 0.42; % CHP气电转换效率
k_eta_CHP_h = 0.29; % CHP气热转换效率
k_eta_gb = 0.82; % 燃气锅炉效率
k_eta_eb = 0.97; % 电锅炉效率
k_eta_ess_cha = 0.95; % 蓄电池充电效率
k_eta_ess_dis = 0.93; % 蓄电池放电效率
k_eta_hss_cha = 0.95; % 储热罐充电效率
k_eta_hss_dis = 0.93; % 储热罐放电效率
D_self_elec = 0.0001; % 电储能自放电率 (占容量的比例)
D_self_heat = 0.0002; % 热储能自放电率 (占容量的比例)
% 系统运行参数
N_Hours = 24; % 时段数
Ramp_elec_th = 20; % 火电电出力爬坡率 (MW/h)
Ramp_heat_th = 15; % 火电热出力爬坡率 (MW/h)
% 火电机组煤耗系数
k_c_xs = [0.0004818, 2.695, 52, 0.0004751, 0.3536, 0.001004];
% CCS系统参数
k_ccs_co2 = 20.23; % CO2处理系数
v_ccs_int = 6000; % CCS初始体积 (m³)
CCS_ab_max = 120; % CCS最大吸收量 (吨/小时)
P_CCS_base = 2.51; % CCS基础能耗 (MW)
P_CCS_ton = 0.269; % CCS单位处理量能耗 (MW/吨)
% 初始状态参数
P_TH_int = 50.73; % 火电机组初始电出力
H_TH_int = 30.11; % 火电机组初始热出力
%% ================== 决策变量定义 ==================
% 定义决策变量
P_TH_e = sdpvar(1, N_Hours); % 火电机组电出力
H_TH_h = sdpvar(1, N_Hours); % 火电机组热出力
P_CHP_e = sdpvar(1, N_Hours); % CHP机组电出力
H_CHP_h = sdpvar(1, N_Hours); % CHP机组热出力
H_GB_h = sdpvar(1, N_Hours); % 燃气锅炉热出力
H_EB_h = sdpvar(1, N_Hours); % 电锅炉热出力
P_g_CHP = sdpvar(1, N_Hours); % CHP机组消耗燃气
P_g_GB = sdpvar(1, N_Hours); % 燃气锅炉耗气
P_e_EB = sdpvar(1, N_Hours); % 电锅炉耗电
P_e_CCS = sdpvar(1, N_Hours); % CCS系统耗电
P_buy_e = sdpvar(1, N_Hours); % 上网购电
P_pv_e = sdpvar(1, N_Hours); % 光伏实际出力
P_wt_e = sdpvar(1, N_Hours); % 风电实际出力
% 储能变量
P_ESS_cha = sdpvar(1, N_Hours); % 电储能充电
P_ESS_dis = sdpvar(1, N_Hours); % 电储能放电
H_HSS_cha = sdpvar(1, N_Hours); % 储热罐充电功率
H_HSS_dis = sdpvar(1, N_Hours); % 储热罐放电功率
E_SOC_elec = sdpvar(1, N_Hours+1); % 储电状态
H_SOC_heat = sdpvar(1, N_Hours+1); % 储热状态
% CCS系统变量
M_CCS_ab = sdpvar(1, N_Hours); % CCS吸收处理量 (吨)
M_CCS_st = sdpvar(1, N_Hours); % CCS解吸处理量 (吨)
T_v_rich = sdpvar(1, N_Hours+1); % 富液罐体积 (m³)
T_v_lean = sdpvar(1, N_Hours+1); % 贫液罐体积 (m³)
% 状态变量
B_ess_s = binvar(1, N_Hours); % 储电状态 (0=充电,1=放电)
B_hss_s = binvar(1, N_Hours); % 储热状态 (0=充电,1=放电)
%% ================== 约束条件 ==================
Cons = [];
% 设置合理初始状态
E_SOC_int = 0.53 * SOC_ESS_max; % 初始电储能状态
H_SOC_int = 0.47 * SOC_HSS_max; % 初始热储能状态
T_v_rich(1) = v_ccs_int; % CCS初始富液罐体积
T_v_lean(1) = v_ccs_int; % CCS初始贫液罐体积
Cons = [Cons, E_SOC_elec(1) == E_SOC_int];
Cons = [Cons, H_SOC_heat(1) == H_SOC_int];
for i = 1:N_Hours
% ==== 设备运行约束 ====
% CHP机组约束
Cons = [Cons, 0 <= P_g_CHP(i) <= P_CHP_max];
Cons = [Cons, P_CHP_e(i) == k_eta_CHP_e * P_g_CHP(i)];
Cons = [Cons, H_CHP_h(i) == k_eta_CHP_h * P_g_CHP(i)];
% 燃气锅炉约束
Cons = [Cons, 0 <= P_g_GB(i) <= P_GB_max];
Cons = [Cons, H_GB_h(i) == k_eta_gb * P_g_GB(i)];
Cons = [Cons, 0 <= P_e_EB(i) <= P_EB_max];
Cons = [Cons, H_EB_h(i) == k_eta_eb * P_e_EB(i)];
% ==== 火电机组约束 (热-电联产模式) ====
% 电出力范围
Cons = [Cons, 50.73 <= P_TH_e(i) <= 135.54];
% 热出力范围
Cons = [Cons, 30.11 <= H_TH_h(i) <= 98.18];
% 热-电耦合约束
Cons = [Cons, P_TH_e(i) >= 1.64*H_TH_h(i) - 50.56];
Cons = [Cons, P_TH_e(i) <= 135.54 - 0.25*H_TH_h(i)];
% ==== 储能约束 ====
Cons = [Cons, 0 <= P_ESS_cha(i) <= P_ess_max * (1 - B_ess_s(i))];
Cons = [Cons, 0 <= P_ESS_dis(i) <= P_ess_max * B_ess_s(i)];
Cons = [Cons, 0 <= H_HSS_cha(i) <= H_hss_max * (1 - B_hss_s(i))];
Cons = [Cons, 0 <= H_HSS_dis(i) <= H_hss_max * B_hss_s(i)];
% ==== CCS约束 ====
% 吸收处理量约束
Cons = [Cons, 0 <= M_CCS_ab(i) <= CCS_ab_max];
% 解吸处理量约束
Cons = [Cons, 0 <= M_CCS_st(i) <= CCS_ab_max];
% 富液罐体积约束
Cons = [Cons, T_v_rich(i+1) == T_v_rich(i) + k_ccs_co2*M_CCS_ab(i) - k_ccs_co2*M_CCS_st(i)];
Cons = [Cons, 0 <= T_v_rich(i+1) <= 12000];
% 贫液罐体积约束
Cons = [Cons, T_v_lean(i+1) == T_v_lean(i) - k_ccs_co2*M_CCS_ab(i) + k_ccs_co2*M_CCS_st(i)];
Cons = [Cons, 0 <= T_v_lean(i+1) <= 12000];
% CCS能耗约束
Cons = [Cons, P_e_CCS(i) == P_CCS_base + P_CCS_ton * M_CCS_st(i)];
% ==== 各类能量平衡 ====
% 电功率平衡 (增加CCS能耗)
Cons = [Cons,
P_TH_e(i) + P_CHP_e(i) + P_wt_e(i) + P_pv_e(i) + P_buy_e(i) + P_ESS_dis(i) == P_Load_e(i) + P_ESS_cha(i) + P_e_EB(i) + P_e_CCS(i)
];
% 热功率平衡
Cons = [Cons,
H_TH_h(i) + H_CHP_h(i) + H_GB_h(i) + H_EB_h(i) + H_HSS_dis(i) == H_Load_h(i) + H_HSS_cha(i)
];
% 新能源出力限制
Cons = [Cons, 0 <= P_pv_e(i) <= P_PV_max(i)];
Cons = [Cons, 0 <= P_wt_e(i) <= P_WT_max(i)];
Cons = [Cons, 0 <= P_buy_e(i) <= 200];
% ==== 储能状态更新 ====
Cons = [Cons, E_SOC_elec(i+1) == (E_SOC_elec(i) - D_self_elec * SOC_ESS_max) + k_eta_ess_cha*P_ESS_cha(i) - P_ESS_dis(i)/k_eta_ess_dis];
Cons = [Cons, 0.17*SOC_ESS_max <= E_SOC_elec(i+1) <= 0.93*SOC_ESS_max];
Cons = [Cons, H_SOC_heat(i+1) == (H_SOC_heat(i) - D_self_heat * SOC_HSS_max) + k_eta_hss_cha*H_HSS_cha(i) - H_HSS_dis(i)/k_eta_hss_dis];
Cons = [Cons, 0.23*SOC_HSS_max <= H_SOC_heat(i+1) <= 0.95*SOC_HSS_max];
end
% ===== 火电机组爬坡约束 =====
% 初始时段基于实际初始值
Cons = [Cons, P_TH_e(1) >= P_TH_int - Ramp_elec_th];
Cons = [Cons, P_TH_e(1) <= P_TH_int + Ramp_elec_th];
Cons = [Cons, H_TH_h(1) >= H_TH_int - Ramp_heat_th];
Cons = [Cons, H_TH_h(1) <= H_TH_int + Ramp_heat_th];
% 后续时段爬坡率
for i = 2:N_Hours
Cons = [Cons, -Ramp_elec_th <= P_TH_e(i)-P_TH_e(i-1) <= Ramp_elec_th];
Cons = [Cons, -Ramp_heat_th <= H_TH_h(i)-H_TH_h(i-1) <= Ramp_heat_th];
end
% ===== 储能结束状态约束 =====
Cons = [Cons, 0.45*SOC_HSS_max <= H_SOC_heat(end) <= 0.55*SOC_HSS_max];
Cons = [Cons, 0.45*SOC_ESS_max <= E_SOC_elec(end) <= 0.55*SOC_ESS_max];
% ===== CCS结束状态约束 =====
Cons = [Cons, v_ccs_int - 1000 <= T_v_rich(end) <= v_ccs_int + 1000];
Cons = [Cons, v_ccs_int - 1000 <= T_v_lean(end) <= v_ccs_int + 1000];
%% ================== 目标函数 ==================
F1_fuel_cost = 0;
F2_carbon_cost = 0;
F3_penalty_cost = 0;
for i = 1:N_Hours
% 火电煤耗计算
Coal_th_kg(i) = k_c_xs(1)*P_TH_e(i)^2 + k_c_xs(2)*P_TH_e(i) + k_c_xs(3)*H_TH_h(i)^2 + k_c_xs(4)*H_TH_h(i) + k_c_xs(5)*P_TH_e(i)*H_TH_h(i) + k_c_xs(6);
Coal_th_ton(i) = Coal_th_kg(i) / 1000;
% 燃气消耗量 (m³)
Gas_volume(i) = (P_g_CHP(i) + P_g_GB(i)) / G_heat_v;
% 总燃料成本计算
F1_fuel_cost = F1_fuel_cost + (C_p_coal * Coal_th_ton(i) + C_p_gas * Gas_volume(i));
end
for i = 1:N_Hours
% 碳排放计算(煤炭碳排、天然气碳排)
Coal_emission(i) = 2.67 * Coal_th_ton(i); % 煤炭排放 (吨CO2)
Gas_emission(i) = Gas_volume(i) * 1.96 / 1000; % 燃气排放 (吨CO2)
Total_emission(i) = Coal_emission(i) + Gas_emission(i);
% 碳配额计算 (基于电负荷计算)
Carbon_quota_per = 0.57; % 每MWh电量分配0.57吨CO2配额
Carbon_quota(i) = Carbon_quota_per * P_Load_e(i);
% 碳交易成本(考虑CCS捕集量)
F2_carbon_cost = F2_carbon_cost + C_p_carbon * (Total_emission(i) - Carbon_quota(i) - M_CCS_ab(i));
end
for i = 1:N_Hours
% 购电成本计算
F3_penalty_cost = F3_penalty_cost + (C_buy_e(i) * P_buy_e(i));
end
% 设备固定成本 (按年分摊)
days_per_year = 365;
C_fixed_th = (4.589e6 * 135.54) / days_per_year; % 火电机组固定成本
C_fixed_chp = (3.2e6 * 130) / days_per_year; % CHP机组固定成本
C_fixed_gb = (700e3 * P_GB_max) / days_per_year; % 燃气锅炉固定成本
C_fixed_eb = (800e3 * P_EB_max) / days_per_year; % 电锅炉固定成本
C_fixed_ess = (2e6 * SOC_ESS_max) / (days_per_year); % 电储能固定成本
C_fixed_hss = (1.5e6 * SOC_HSS_max) / (days_per_year); % 热储能固定成本
C_fixed_ccs = (3.5e6 * C_ccs_fixed) / (days_per_year); % CCS系统固定成本
F4_Com_fixed = C_fixed_th + C_fixed_chp + C_fixed_gb + C_fixed_eb + C_fixed_ess + C_fixed_hss + C_fixed_ccs;
% 总目标函数
Objective = F1_fuel_cost + F2_carbon_cost + F3_penalty_cost + F4_Com_fixed;
%% ================== CPLEX求解模型 ==================
ops = sdpsettings('solver','cplex', 'verbose', 1);
sol = optimize(Cons, Objective, ops);
if sol.problem == 0
disp('==================== 优化成功 ====================');
disp(['>>>总成本: ', num2str(value(Objective)/10000, '%.2f'), ' 万元']);
% 成本分解
cost_fuel = value(F1_fuel_cost);
cost_carbon = value(F2_carbon_cost);
cost_penalty = value(F3_penalty_cost);
cost_fixed = value(F4_Com_fixed);
cost_total = value(Objective);
disp('========== 运行成本分解 ==========');
disp([' 燃料成本: ', num2str(cost_fuel/10000, '%.2f'), ' 万元']);
disp([' 购电成本: ', num2str(cost_penalty/10000, '%.2f'), ' 万元']);
disp([' 固定成本: ', num2str(cost_fixed/10000, '%.2f'), ' 万元']);
disp([' 碳交易成本: ', num2str(cost_carbon/10000, '%.2f'), ' 万元']);
% CCS捕集量统计
Total_emission = value(Total_emission);
M_CCS_ab = value(M_CCS_ab);
disp(['>>>总碳排放量: ', num2str(sum(Total_emission)), ' 吨']);
disp(['>>>总碳捕集量: ', num2str(sum(M_CCS_ab)), ' 吨']);
else
disp('========== 优化失败 ==========');
disp(['错误代码: ', num2str(sol.problem)]);
disp(['错误信息: ', sol.info]);
yalmiperror(sol.problem);
end
%% ================== 结果可视化 ==================
% 提取变量值
P_buy_e = value(P_buy_e);
P_ESS_cha = value(P_ESS_cha);
P_ESS_dis = value(P_ESS_dis);
P_e_EB = value(P_e_EB);
P_pv_e = value(P_pv_e);
P_wt_e = value(P_wt_e);
P_g_CHP = value(P_g_CHP);
P_g_GB = value(P_g_GB);
M_CCS_ab = value(M_CCS_ab);
P_e_CCS = value(P_e_CCS);
H_HSS_cha = value(H_HSS_cha);
H_HSS_dis = value(H_HSS_dis);
E_SOC_elec = value(E_SOC_elec);
H_SOC_heat = value(H_SOC_heat);
P_TH_e = value(P_TH_e);
H_TH_h = value(H_TH_h);
P_CHP_e = k_eta_CHP_e * P_g_CHP;
H_CHP_h = k_eta_CHP_h * P_g_CHP;
H_GB_h = k_eta_gb * P_g_GB;
H_EB_h = k_eta_eb * P_e_EB;
% 新能源出力曲线
figure(1);
x = 1:24;
plot(x, P_pv_e, '-r*');
hold on;
plot(x, P_wt_e, '-b*');
grid on;
legend('光伏出力', '风电出力');
xlabel('时间 (小时)');
ylabel('出力 (MW)');
title('新能源机组出力曲线');
% 电功率平衡图
figure(2);
Figure_E_in = [P_buy_e', P_pv_e', P_wt_e', P_TH_e', P_CHP_e', P_ESS_dis'];
Figure_E_out = [-P_ESS_cha', -P_e_EB', -P_e_CCS'];
e_in = bar(1:24, Figure_E_in, 'stacked');
hold on;
e_out = bar(1:24, Figure_E_out, 'stacked');
hold on;
plot(x, P_Load_e, 'g-*');
grid on;
colors_in = {
[0.2 0.6 0.8] % 上网购电 - 蓝色
[1 0.8 0.2] % 光伏机组 - 黄色
[0.5 0.8 0.2] % 风电机组 - 绿色
[0.9 0.4 0.2] % 火电机组 - 橙色
[0.7 0.2 0.7] % CHP机组 - 紫色
[0.1 0.1 0.6] % 储能放电 - 深蓝
};
for i = 1:length(e_in)
set(e_in(i), 'FaceColor', colors_in{i}, 'EdgeColor', 'none');
end
colors_out = {
[0.6 0.6 0.6] % 储能充电 - 灰色
[0.3 0.3 0.3] % CCS能耗 - 深灰
[0.8 0.3 0.3] % 电锅炉 - 红色
};
for i = 1:length(e_out)
set(e_out(i), 'FaceColor', colors_out{i}, 'EdgeColor', 'none');
end
legend('上网购电', '光伏机组', '风电机组', '火电机组', 'CHP机组', '储能放电', '储能充电', 'CCS能耗', '电锅炉', '电负荷');
xlabel('时间 (小时)');
ylabel('功率 (MW)');
title('电功率平衡图');
% 热功率平衡图
figure(3);
Figure_H_in = [H_TH_h', H_CHP_h', H_GB_h', H_EB_h', H_HSS_dis'];
Figure_H_out = [-H_HSS_cha'];
h_in = bar(1:24, Figure_H_in, 'stacked');
hold on;
h_out = bar(1:24, Figure_H_out, 'stacked');
hold on;
plot(x, H_Load_h, 'r-*');
grid on;
colors_heat = {
[0.9 0.4 0.2] % 火电供热 - 橙色
[0.7 0.2 0.7] % CHP供热 - 紫色
[0.3 0.6 1] % 燃气锅炉 - 蓝色
[0.8 0.3 0.3] % 电锅炉 - 红色
[0.1 0.1 0.6] % 储热罐放热 - 深蓝
[0.5 0.8 0.2] % 储热罐充电 - 绿色
};
for i = 1:length(h_in)
h_in(i).FaceColor = colors_heat{i};
h_in(i).EdgeColor = 'none';
end
for i = 1:length(h_out)
h_out(i).FaceColor = colors_heat{length(h_in)+i};
h_out(i).EdgeColor = 'none';
end
legend('火电机组', 'CHP机组', '燃气锅炉', '电锅炉', '储罐放热', '储罐充热', '热负荷');
xlabel('时间 (小时)');
ylabel('功率 (MW)');
title('热功率平衡图');
% 储能状态变化
figure(4);
subplot(2,1,1);
plot(0:24, E_SOC_elec, 'b-*');
title('电储能状态变化');
xlabel('时间 (小时)');
ylabel('储电量 (MWh)');
grid on;
ylim([0 SOC_ESS_max]);
subplot(2,1,2);
plot(0:24, H_SOC_heat, 'r-*');
title('热储能状态变化');
xlabel('时间 (小时)');
ylabel('储热量 (MWh)');
grid on;
ylim([0 SOC_HSS_max]);
% 火电机组运行点
figure(5);
H_range = linspace(30.11, 98.18, 100);
P_min_bound = 1.64*H_range - 50.56;
P_max_bound = 135.54 - 0.25*H_range;
fill([H_range, fliplr(H_range)], [P_min_bound, fliplr(P_max_bound)], [0.8 0.9 0.8]);
hold on;
plot(H_TH_h, P_TH_e, 'r*-');
xlabel('热出力 (MW)');
ylabel('电出力 (MW)');
title('火电机组运行轨迹');
legend('可行域', '运行点');
grid on;
axis([20,100,-20,160]);
% CCS运行状态
figure(6);
subplot(2,1,1);
bar_data = [M_CCS_ab', M_CCS_st'];
h = bar(1:24, bar_data, 'stacked');
set(h(1), 'FaceColor', [0.8 0.2 0.2], 'EdgeColor', 'none');
set(h(2), 'FaceColor', [0.9 0.6 0.1], 'EdgeColor', 'none');
title('CCS碳捕集与再生量');
xlabel('时间 (小时)');
ylabel('CO2量 (吨)');
legend('吸收量', '解吸量');
grid on;
subplot(2,1,2);
plot(1:24, value(T_v_rich(1:24)), 'b-*');
hold on;
plot(1:24, value(T_v_lean(1:24)), 'r-*');
ylim([-10 13000]);
title('CCS富液罐和贫液罐体积');
xlabel('时间 (小时)');
ylabel('体积 (m³)');
legend('富液罐', '贫液罐');
grid on;