修改和优化以下代码：clc
clear
%% ===================== 数据导入 =====================
% 风电预测功率
P_wt_pre = [419.04,374.32,392.12,398.68,460.74,405.60,235.60,222.69,260.45,277.39,183.57,120.09,81.92,127.76,163.73,195.31,209.13,267.56,374.10,407.58,421.39,468.05,494.70,479.29];
% 电负荷数据
P_Load_e = [483.652,556.596,567.088,572.472,559.540,528.948,564.380,597.140,621.736,556.064,708.320,714.772,604.900,522.044,484.160,491.404,576.340,708.172,787.652,674.048,595.500,578.148,514.456,550.844];
% 热负荷数据
P_Load_h = [358.680,358.420,370.790,369.580,344.640,321.890,352.310,326.540,315.740,265.470,232.010,227.020,296.580,292.170,267.080,273.990,288.170,270.460,227.610,345.550,322.250,342.400,347.970,370.340];
%% ===================== 参数设置 =====================
% ===== 火电机组参数 =====
P_MT_max = [400,455,200]; % 机组最大出力 (MW)
P_MT_min = [200,120,100]; % 机组最小出力 (MW)
Rate_u = [200; 200; 100]; % 机组爬坡率 (MW/h)
a = [0.00048,0.00031,0.0002]; % 燃料成本二次项系数 (元/MW²)
b = [16.2,17.3,16.6]; % 燃料成本一次项系数 (元/MW)
c = [1000,970,700]; % 燃料成本常数项 (元)
C_gas_price = 50; % 天然气价格 (元/单位)
% ===== 碳捕集参数 =====
C_k_eg = [0.91,0.95,0.98]; % 机组碳排放系数 (t_CO2/MWh)
P_lamda = [0.18,0.18,0.18]; % 单位捕碳量能耗 (MWh/t_CO2)
P_yita = [1.05,1.05,1.05]; % 再生塔效率系数 (关键参数)
P_GD_eg = [15,15,10]; % 碳捕集固定能耗 (MW)
C_cb_price = 120; % 碳交易价格 (元/t_CO2)
f_quota_rate = 0.65; % 免费配额比例
E_bata_c = 0.81; % 碳捕集的效率
% ===== 溶剂系统参数 =====
M_MEA = 61.08; % MEA摩尔质量 (g/mol)
M_co2 = 44; % CO2摩尔质量 (g/mol)
C_s_CR = 30; % 醇胺溶液浓度 (%)
C_rou_gL = 1010; % 溶液密度 (g/L)
B_V_CR = [2500,2500,2500]; % 溶液储罐容量 (L)
K_k_R = 0.95; % 溶剂成本系数 (元/L)
C_fai = 1.21; % 溶剂运行损耗系数
% ===== 需求响应参数 =====
DR_e_tra_limit = 0.12; % 电负荷转移比例上限
DR_h_tra_limit = 0.12; % 热负荷转移比例上限
DR_e_cut_limit = 0.03; % 电负荷削减比例上限
DR_h_cut_limit = 0.03; % 热负荷削减比例上限
% ===== 高峰时段定义 =====
e_peak_hours = [10,11,12,18,19,20]; % 电负荷高峰时段
h_peak_hours = [1,2,3,21,22,23]; % 热负荷高峰时段
% ===== 新增参数 =====
C_wt_om = 15; % 风电单位运维成本 [元/MWh]
%% ===================== 决策变量定义 =====================
% --- 能源生产变量 ---
P_CHP_e = sdpvar(1,24,'full'); % CHP机组电出力 (MW)
P_MT_e = sdpvar(3,24,'full'); % 火电机组出力 (MW)
P_wt_e = sdpvar(1,24,'full'); % 风电实际出力 (MW)
P_e_EB = sdpvar(1,24,'full'); % 电锅炉耗电量 (MW)
P_e_CCS = sdpvar(3,24,'full'); % 碳捕集耗电量 (MW)
P_CHP_h = sdpvar(1,24,'full'); % CHP机组热出力 (MW)
P_EB_h = sdpvar(1,24,'full'); % 电锅炉热出力 (MW)
% --- 天然气购买 ---
P_g_CHP = sdpvar(1,24,'full'); % 天然气消耗量 (单位)
% --- 碳捕集系统 ---
E_Amt_co2 = sdpvar(3,24,'full'); % 机组碳排放量 (t_CO2)
E_co2_CCS = sdpvar(3,24,'full'); % 捕获的CO2量 (t_CO2)
F_CO2_aux = sdpvar(3,24,'full'); % 超额碳排放量 (t_CO2)
% --- 溶剂系统 ---
B_V_CA = sdpvar(3,24,'full'); % 溶剂循环量 (L)
B_V_FY = sdpvar(3,24,'full'); % 富液体积 (L)
B_V_PY = sdpvar(3,24,'full'); % 贫液体积 (L)
% --- 需求响应 ---
P_tra_e_pos = sdpvar(1,24,'full'); % 电负荷正向转移量
P_tra_e_neg = sdpvar(1,24,'full'); % 电负荷负向转移量
P_cut_e = sdpvar(1,24,'full'); % 可削减电负荷
P_DR_Load_e = sdpvar(1,24,'full'); % 需求响应后电负荷
H_tra_h_pos = sdpvar(1,24,'full'); % 热负荷正向转移量
H_tra_h_neg = sdpvar(1,24,'full'); % 热负荷负向转移量
H_cut_h = sdpvar(1,24,'full'); % 可削减热负荷
P_DR_Load_h = sdpvar(1,24,'full'); % 需求响应后热负荷
%% ===================== 约束条件 =====================
Cons = [];
% ===== 能源平衡约束 =====
for t = 1:24
% 电功率平衡
Cons = [Cons,
P_CHP_e(t) + P_wt_e(t) + sum(P_MT_e(:, t)) == P_DR_Load_e(t) + P_e_EB(t) + sum(P_e_CCS(:, t))
];
% 热功率平衡
Cons = [Cons,
P_CHP_h(t) + P_EB_h(t) == P_DR_Load_h(t) % 产热=用热
];
end
% ===== 机组运行约束 =====
for t = 1:24
% 风电约束：实际出力不超过预测值
Cons = [Cons,
0 <= P_wt_e(t) <= P_wt_pre(t)
];
% CHP机组约束：效率模型
Cons = [Cons,
0 <= P_g_CHP(t) <= 850, % 天然气消耗范围
P_CHP_e(t) == 0.48 * P_g_CHP(t), % 电效率48%
P_CHP_h(t) == 0.32 * P_g_CHP(t) % 热效率32%
];
% 电锅炉约束：电热转换效率
Cons = [Cons,
P_EB_h(t) == 0.97 * P_e_EB(t) % 热效率97%
];
% 火电机组约束 (含碳捕集系统)
for i = 1:3
% 总碳排放计算：发电量×排放系数
Cons = [Cons,
E_Amt_co2(i, t) == P_MT_e(i, t) * C_k_eg(i)
];
% 碳捕集过程 (物理模型)
Cons = [Cons,
% 碳捕集量 = 捕集效率 × 待捕集量
E_co2_CCS(i, t) == E_bata_c * E_Amt_co2(i, t),
% 捕集量上下限约束
0 <= E_co2_CCS(i, t) <= P_yita(i) * E_bata_c * C_k_eg(i) * P_MT_max(i),
E_co2_CCS(i, t) <= E_Amt_co2(i, t), % 捕集量不超过总排放
% 碳捕集能耗 = 单位捕集能耗 × 碳捕集量 + 固定能耗
P_e_CCS(i, t) == P_lamda(i) * E_co2_CCS(i, t) + P_GD_eg(i),
% 火电机组出力范围
P_MT_min(i) <= P_MT_e(i, t) <= P_MT_max(i),
];
end
end
% ===== 机组爬坡约束 =====
for t = 1:23
% CHP机组爬坡：
Cons = [Cons,
-150<= P_CHP_e(t+1) - P_CHP_e(t) <= 150
];
% 电锅炉爬坡：
Cons = [Cons,
-130 <= P_e_EB(t+1) - P_e_EB(t) <= 130
];
% 火电机组爬坡：
for i = 1:3
Cons = [Cons,
-Rate_u(i) <= P_MT_e(i, t+1) - P_MT_e(i, t) <= Rate_u(i)
];
end
end
% ===== 溶剂系统动态约束 =====
B_V_PY0 = [1000; 1000; 1000]; % 贫液初始量 (L)
B_V_FY0 = [1000; 1000; 1000]; % 富液初始量 (L)
for i = 1:3
% t=1时段特殊处理：初始状态
Cons = [Cons,
B_V_FY(i, 1) == B_V_FY0(i) + B_V_CA(i, 1), % 富液=初始+新增
B_V_PY(i, 1) == B_V_PY0(i) - B_V_CA(i, 1) % 贫液=初始-消耗
];
% t=2-24时段：溶剂流动平衡
for t = 2:24
Cons = [Cons,
B_V_FY(i, t) == B_V_FY(i, t-1) + B_V_CA(i, t), % 富液积累
B_V_PY(i, t) == B_V_PY(i, t-1) - B_V_CA(i, t) % 贫液消耗
];
end
% 溶液容量约束：储罐限制
Cons = [Cons,
0 <= B_V_FY(i, :) <= B_V_CR(i), % 富液容量
0 <= B_V_PY(i, :) <= B_V_CR(i) % 贫液容量
];
% 周期平衡约束 (允许±50L偏差)
Cons = [Cons,
abs(B_V_PY(i, 24) - B_V_PY0(i)) <= 50, % 贫液周期平衡
abs(B_V_FY(i, 24) - B_V_FY0(i)) <= 50 % 富液周期平衡
];
% ===== 新增溶剂循环量约束 =====
% 溶剂循环量与CO2捕集量的化学平衡
for t = 1:24
% 计算单位CO2捕集所需溶剂体积 [L/t_CO2]
solvent_per_co2 = (M_MEA / M_co2) / (C_s_CR/100 * C_rou_gL) * 1e6;
Cons = [Cons,
B_V_CA(i, t) == E_co2_CCS(i, t) * solvent_per_co2
];
end
end
% ===== 碳交易约束 =====
for i = 1:3
for t = 1:24
% 实际排放量 = 总碳排放 - 碳捕集量
Actual_emission = E_Amt_co2(i,t) - E_co2_CCS(i,t);
% 免费碳配额 = 配额率 × 排放系数 × 发电量
Free_quota = f_quota_rate * C_k_eg(i) * P_MT_e(i,t);
% 超额排放量计算
Cons = [Cons,
F_CO2_aux(i,t) >= Actual_emission - Free_quota, % 超额排放部分
F_CO2_aux(i,t) >= 0 % 非负约束
];
end
end
% ===== 需求响应约束 =====
% 转移负荷限制：每小时转移量上限
for t = 1:24
Cons = [Cons,
% 电负荷转移量不超过比例限制
0 <= P_tra_e_pos(t) <= DR_e_tra_limit * P_Load_e(t),
0 <= P_tra_e_neg(t) <= DR_e_tra_limit * P_Load_e(t),
% 热负荷转移量不超过比例限制
0 <= H_tra_h_pos(t) <= DR_h_tra_limit * P_Load_h(t),
0 <= H_tra_h_neg(t) <= DR_h_tra_limit * P_Load_h(t)
];
end
% 转移总量平衡 (允许1%偏差)
Cons = [Cons,
abs(sum(P_tra_e_pos) - sum(P_tra_e_neg)) <= 0.01*sum(P_Load_e),
abs(sum(H_tra_h_pos) - sum(H_tra_h_neg)) <= 0.01*sum(P_Load_h)
];
% 可削减负荷约束：仅高峰时段允许
for t = 1:24
% 电负荷削减约束
if ismember(t, e_peak_hours)
Cons = [Cons,
-DR_e_cut_limit*P_Load_e(t) <= P_cut_e(t) <= 0
];
else
Cons = [Cons,
P_cut_e(t) == 0 % 非高峰时段不能削减
];
end
% 热负荷削减约束
if ismember(t, h_peak_hours)
Cons = [Cons,
-DR_h_cut_limit*P_Load_h(t) <= H_cut_h(t) <= 0
];
else
Cons = [Cons,
H_cut_h(t) == 0 % 非高峰时段不能削减
];
end
end
% 响应后负荷计算
Cons = [Cons,
% 电负荷 = 原始负荷 + 净转移量 + 削减量
P_DR_Load_e == P_Load_e + (P_tra_e_pos - P_tra_e_neg) + P_cut_e,
% 热负荷 = 原始负荷 + 净转移量 + 削减量
P_DR_Load_h == P_Load_h + (H_tra_h_pos - H_tra_h_neg) + H_cut_h
];
%% ===================== 目标函数 =====================
% 总目标：最小化运行成本 = 燃料成本 + 燃气成本 + 环境成本
% 1. 燃料成本 (火电机组二次函数)
F_1 = 0;
for i = 1:3
for t = 1:24
F_1 = F_1 + a(i)*P_MT_e(i,t)^2 + b(i)*P_MT_e(i,t) + c(i);
end
end
% 2. 燃气成本 (CHP机组)
F_2 = C_gas_price * sum(P_g_CHP);
% 3. 环境成本 (溶剂成本 + 碳交易成本)
F_3 = 0;
for i = 1:3
for t = 1:24
% 溶剂成本 = 成本系数 × 损耗系数 × 溶液量
F_3 = F_3 + K_k_R * C_fai * B_V_CA(i,t);
% 碳交易成本 = 碳价 × 超额排放量
F_3 = F_3 + C_cb_price * F_CO2_aux(i,t);
end
end
% 总成本
F_all = F_1 + F_2 + F_3 ;
最新发布
博客记录了日常琐事，提到课程无趣，当天搞定上网、听歌和 flash，晚上和爱人逛超市很累，扬华网挂了打算明天再看，还计划列个国庆计划。
扫一扫
764
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
