检查和优化重构以下代码:%% 日前调度优化模型
clc; clear;
%% 参数初始化
T = 24; % 调度周期
% 预测数据
P_wt_pre = 1.4 * [299.3168, 267.3737, 280.0854, 284.7727, 329.0996, 289.7108, 168.2844, 159.0641, ...
186.0346, 198.1329, 131.1193, 85.7753, 58.5152, 91.2576, 116.9506, 139.5097, ...
149.3775, 191.1110, 267.2125, 291.1269, 300.9949, 334.3191, 353.3572, 342.3515];
P_Load_e = 0.74 * [1209.1268, 1391.4908, 1417.7240, 1431.1808, 1398.8520, 1322.3700, 1410.9480, 1492.8480, ...
1554.3360, 1390.1580, 1770.8040, 1786.9320, 1512.2520, 1305.1080, 1210.3980, 1228.5140, ...
1440.8520, 1770.4260, 1969.1280, 1685.1240, 1488.7460, 1445.3740, 1286.1380, 1377.1100];
P_Load_h = 0.7 * [358.6767, 358.4154, 370.7885, 369.5773, 344.6412, 321.8899, 352.3120, 326.5447, ...
315.7443, 265.4748, 232.0104, 227.0214, 296.5760, 292.1705, 267.0795, 273.9941, ...
288.1734, 270.4585, 227.6077, 345.5505, 322.2461, 342.3991, 347.9680, 370.3373];
gama_elec = 100 * [0.29, 0.29, 0.29, 0.29, 0.29, 0.69, 0.69, 0.69, 0.84, 0.84, 0.84, 0.69, ...
0.69, 0.69, 0.69, 0.69, 0.69, 0.69, 0.64, 0.64, 0.55, 0.55, 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]; % 火电机组最小出力
R_u = 4 * [50, 50, 25]; % 机组爬坡率
% 碳捕集系统参数
P_D1 = 0.5 * [20, 20, 15]; % 碳捕集固定能耗
P_D = repmat(P_D1', 1, T); % 扩展至24小时
M_MEA = 61.08; % M_MEA摩尔质量
M_co2 = 44; % CO2摩尔质量
sita = 0.4; % 再生塔解析量
CR = 30; % 溶液浓度(%)
rou_R = 1.01; % 溶液密度
V_CR = [2000, 2000, 2000]; % 溶液储液装置容量
K_R = 1.17; % 溶剂成本系数
fai = 1.5; % 溶剂损耗系数
% 分段线性化参数
h_gn = 5; % 上分段数
k_gn = 5; % 下分段数
gn = 5; % 火电机组分段数
% 计算分段边界
hl1 = (P_G_max - P_G_min) / h_gn; % 上分段步长
kl1 = (P_G_min - P_G_max) / k_gn; % 下分段步长 (注意负号)
gl1 = (P_G_max - P_G_min) / gn; % 机组分段步长
%% 决策变量定义
% 电力相关
GT_P = sdpvar(2, T, 'full'); % 燃气轮机电出力
P_w = sdpvar(1, T, 'full'); % 风电出力
P_G = sdpvar(3, T, 'full'); % 火电机组总出力
P_J = sdpvar(3, T, 'full'); % 火电机组净出力
EB = sdpvar(2, T, 'full'); % 电锅炉出力
% 热力相关
GT_H = sdpvar(2, T, 'full'); % 燃气轮机热出力
EB_H = sdpvar(2, T, 'full'); % 电锅炉热出力
% 碳捕集系统
P_gas = sdpvar(2, T, 'full'); % 天然气需求
E_G = sdpvar(3, T, 'full'); % 总碳排放
E_total_co2 = sdpvar(3, T, 'full'); % 捕获的总碳排放
E_CG = sdpvar(3, T, 'full'); % 待捕集CO2量
P_B = sdpvar(3, T, 'full'); % 机组运行能耗
V_CA = sdpvar(3, T, 'full'); % 溶液体积
V_FY = sdpvar(3, T, 'full'); % 富液体积
V_PY = sdpvar(3, T, 'full'); % 贫液体积
% 需求响应
P_tran = sdpvar(1, T, 'full'); % 可转移电负荷
P_cut = sdpvar(1, T, 'full'); % 可削减电负荷
P_DE = sdpvar(1, T, 'full'); % 响应后电负荷
H_tran = sdpvar(1, T, 'full'); % 可转移热负荷
H_cut = sdpvar(1, T, 'full'); % 可削减热负荷
H_DE = sdpvar(1, T, 'full'); % 响应后热负荷
% 绝对值辅助变量 (用于线性化目标函数中的绝对值)
P_tran_pos = sdpvar(1, T, 'full'); % P_tran正部
P_tran_neg = sdpvar(1, T, 'full'); % P_tran负部
H_tran_pos = sdpvar(1, T, 'full'); % H_tran正部
H_tran_neg = sdpvar(1, T, 'full'); % H_tran负部
% 分段线性化变量
w = binvar(gn, T, 3); % 分段激活标志 (机组×时间×分段)
lambda = sdpvar(gn+1, T, 3);% 分段权重 (分段点×时间×机组)
%% 约束构建
Cons = [];
% ===== 基础约束 =====
% 风电出力约束
Cons = [Cons, 0 <= P_w <= P_wt_pre];
% 电锅炉约束
Cons = [Cons, sum(EB,1) + P_w <= P_wt_pre];
% 火电机组出力约束
for i = 1:3
Cons = [Cons, P_G_min(i) <= P_G(i,:) <= P_G_max(i)];
end
% 绝对值分解约束
Cons = [Cons, P_tran == P_tran_pos - P_tran_neg];
Cons = [Cons, P_tran_pos >= 0, P_tran_neg >= 0];
Cons = [Cons, H_tran == H_tran_pos - H_tran_neg];
Cons = [Cons, H_tran_pos >= 0, H_tran_neg >= 0];
% 碳捕集系统静态约束
for i = 1:3
for t = 1:T
% 总碳排放
Cons = [Cons, E_G(i,t) == eg(i) * P_G(i,t)];
% 捕集量约束
Cons = [Cons, 0 <= E_total_co2(i,t) <= P_yita(i) * E_bata * eg(i) * P_G_max(i)];
Cons = [Cons, E_total_co2(i,t) == E_CG(i,t) + 0.25 * E_bata * eg(i) * (P_G(i,t) - P_J(i,t))];
% 能耗平衡
Cons = [Cons, P_B(i,t) == P_lamda(i) * E_total_co2(i,t)];
Cons = [Cons, P_G(i,t) == P_J(i,t) + P_D(i,t) + P_B(i,t)];
% 净出力范围
P_J_min = P_G_min(i) - P_lamda(i) * P_yita(i) * E_bata * eg(i) * P_G_max(i) - P_D(i,t);
P_J_max = P_G_max(i) - P_D(i,t);
Cons = [Cons, P_J_min <= P_J(i,t) <= P_J_max];
end
end
% 旋转备用约束 (每小时)
for t = 1:T
Cons = [Cons, sum(R_u) >= 0.01 * P_DE(t)];
end
% ===== 分段线性化约束 =====
for i = 1:3 % 对每台机组
seg_points = linspace(P_G_min(i), P_G_max(i), gn+1); % 分段点
for t = 1:T % 每个时间段
% 凸组合约束
Cons = [Cons, sum(lambda(:,t,i)) == 1, lambda(:,t,i) >= 0];
% 功率表达式
Cons = [Cons, P_G(i,t) == seg_points * lambda(:,t,i)];
% 分段激活约束
for seg = 1:gn
Cons = [Cons, lambda(seg,t,i) <= w(seg,t,i)];
if seg < gn
Cons = [Cons, lambda(seg+1,t,i) <= w(seg,t,i) + w(seg+1,t,i)];
end
end
Cons = [Cons, sum(w(:,t,i)) == 1]; % 仅激活一个分段
end
end
% ===== 机组运行约束 =====
% 爬坡约束
for i = 1:3
for t = 1:T-1
ramp_down = -R_u(i);
ramp_up = R_u(i);
Cons = [Cons, ramp_down <= P_G(i,t+1) - P_G(i,t) <= ramp_up];
end
end
% 燃气轮机约束
for j = 1:2
% 出力范围
Cons = [Cons, 0 <= GT_P(j,:) <= GT_P_max(j)];
% 天然气转换
Cons = [Cons, P_gas(j,:) == GT_P(j,:) / 0.98];
Cons = [Cons, 0 <= P_gas(j,:) <= 130];
% 热电关系
Cons = [Cons, GT_H(j,:) == 0.6 * (GT_P(j,:) / 0.4)];
Cons = [Cons, EB_H(j,:) == 0.98 * EB(j,:)];
end
% ===== 碳捕集系统动态约束 =====
V_PY0 = [1000, 1000, 1000]; % 贫液初始值
V_FY0 = [1000, 1000, 1000]; % 富液初始值
for i = 1:3
% 初始条件
Cons = [Cons, V_FY(i,1) == V_FY0(i) - V_CA(i,1)];
Cons = [Cons, V_PY(i,1) == V_PY0(i) + V_CA(i,1)];
% 动态方程
for t = 2:T
Cons = [Cons, V_FY(i,t) == V_FY(i,t-1) - V_CA(i,t)];
Cons = [Cons, V_PY(i,t) == V_PY(i,t-1) + V_CA(i,t)];
end
% 溶液体积计算
for t = 1:T
Cons = [Cons, V_CA(i,t) == (M_MEA * E_CG(i,t)) / (M_co2 * sita * CR * rou_R)];
Cons = [Cons, 0 <= V_FY(i,t) <= V_CR(i)];
Cons = [Cons, 0 <= V_PY(i,t) <= V_CR(i)];
end
end
% ===== 需求响应约束 =====
% 电负荷响应
for t = 1:T
if ismember(t, [11,12,18,19,20]) % 高峰时段
Cons = [Cons, -0.02 * P_Load_e(t) <= P_cut(t) <= 0];
else
Cons = [Cons, P_cut(t) == 0];
end
end
% 热负荷响应
for t = 1:T
if ismember(t, [1,2,3,11,12]) % 高峰时段
Cons = [Cons, -0.02 * P_Load_h(t) <= H_cut(t) <= 0];
else
Cons = [Cons, H_cut(t) == 0];
end
end
% 响应平衡约束
Cons = [Cons, -0.08 * P_Load_e <= P_tran <= 0.08 * P_Load_e];
Cons = [Cons, sum(P_tran) == 0];
Cons = [Cons, P_DE == P_Load_e + P_tran + P_cut];
Cons = [Cons, -0.08 * P_Load_h <= H_tran <= 0.08 * P_Load_h];
Cons = [Cons, sum(H_tran) == 0];
Cons = [Cons, H_DE == P_Load_h + H_tran + H_cut];
% ===== 系统平衡约束 =====
Cons = [Cons, P_w + sum(GT_P,1) + sum(P_J,1) == P_DE]; % 电力平衡
Cons = [Cons, sum(EB_H,1) + sum(GT_H,1) == H_DE]; % 热力平衡
%% 目标函数构建
% 需求响应成本
cost_DR_elec = 0.25 * gama_elec .* (P_tran_pos + P_tran_neg) + ...
0.65 * gama_elec .* (-P_cut); % P_cut<=0
cost_DR_heat = 0.25 * (H_tran_pos + H_tran_neg) + ...
0.65 * (-H_cut); % H_cut<=0
% 发电成本
cost_gen = zeros(1,T);
for i = 1:3
cost_gen = cost_gen + a(i)*P_G(i,:).^2 + b(i)*P_G(i,:) + c(i);
end
% 碳捕集成本
cost_carbon = zeros(1,T);
for i = 1:3
cost_carbon = cost_carbon + 1.5 * K_R * E_total_co2(i,:);
end
% 排放成本
cost_emission = zeros(1,T);
for i = 1:3
net_emission = E_G(i,:) - E_total_co2(i,:);
cost_emission = cost_emission + 12 * net_emission;
end
% 其他成本
cost_gas = 50 * sum(P_gas,1);
cost_EB = gama_elec .* sum(EB,1);
wind_curtail = max(0, P_wt_pre - P_w - sum(EB,1)); % 弃风量
cost_wind_curtail = 0.3 * gama_elec .* wind_curtail;
% 总目标函数
F = sum(cost_DR_elec) + sum(cost_DR_heat) + sum(cost_gen) + ...
sum(cost_carbon) + sum(cost_emission) + sum(cost_gas) + ...
sum(cost_EB) + sum(cost_wind_curtail);
%% 模型求解
ops = sdpsettings('solver', 'cplex', 'verbose', 2, 'usex0', 0);
ops.cplex.mip.tolerances.mipgap = 1e-6;
result = optimize(Cons, F, ops);
%% 结果分析与可视化
if result.problem == 0
disp(['最优值: ', num2str(value(F))]);
% 结果提取
val = @(x) value(x); % 简写取值函数
[GT_P_val, P_w_val, P_G_val] = deal(val(GT_P), val(P_w), val(P_G));
[EB_val, GT_H_val, EB_H_val] = deal(val(EB), val(GT_H), val(EB_H));
[P_J_val, E_G_val, E_total_co2_val] = deal(val(P_J), val(E_G), val(E_total_co2));
[P_DE_val, H_DE_val] = deal(val(P_DE), val(H_DE));
% 可视化设置
t = 1:T;
colors = lines(7); % 获取7种区分度高的颜色
% 1. 系统功率平衡
figure('Name', '系统功率平衡', 'Position', [100, 100, 1200, 800]);
% 电功率平衡
subplot(2,1,1);
hold on;
area(t, [val(P_J(1,:)); val(P_J(2,:)); val(P_J(3,:));
val(P_w); val(GT_P(1,:)); val(GT_P(2,:))]', ...
'LineWidth', 1.5);
plot(t, P_Load_e, 'k--o', 'LineWidth', 2, 'MarkerSize', 8);
plot(t, P_DE_val, 'm-*', 'LineWidth', 2, 'MarkerSize', 8);
legend('火电1', '火电2', '火电3', '风电', '燃机1', '燃机2', ...
'原始电负荷', '响应后电负荷', 'Location', 'best');
title('电功率平衡');
xlabel('时间 (h)'); ylabel('功率 (kW)');
grid on; box on;
% 热功率平衡
subplot(2,1,2);
hold on;
area(t, [val(GT_H(1,:)); val(GT_H(2,:)); val(EB_H(1,:)); val(EB_H(2,:))]', ...
'LineWidth', 1.5);
plot(t, P_Load_h, 'k--o', 'LineWidth', 2, 'MarkerSize', 8);
plot(t, H_DE_val, 'm-*', 'LineWidth', 2, 'MarkerSize', 8);
legend('燃机1热', '燃机2热', '电锅炉1', '电锅炉2', ...
'原始热负荷', '响应后热负荷', 'Location', 'best');
title('热功率平衡');
xlabel('时间 (h)'); ylabel('功率 (kW)');
grid on; box on;
% 2. 成本分解
figure('Name', '成本分解', 'Position', [100, 100, 1000, 600]);
costs = [sum(val(cost_DR_elec)), sum(val(cost_DR_heat)), sum(val(cost_gen)), ...
sum(val(cost_carbon)), sum(val(cost_emission)), ...
sum(val(cost_gas)), sum(val(cost_EB)), sum(val(cost_wind_curtail))];
labels = {'需求响应(电)', '需求响应(热)', '发电成本', '碳捕集', ...
'碳排放', '天然气', '电锅炉', '弃风惩罚'};
pie(costs, labels);
title('总成本构成');
% 3. 碳捕集系统分析
figure('Name', '碳捕集分析', 'Position', [100, 100, 1200, 800]);
% CO2捕集量
subplot(2,1,1);
plot(t, val(E_total_co2(1,:)), 'r-o', 'LineWidth', 1.5); hold on;
plot(t, val(E_total_co2(2,:)), 'b-s', 'LineWidth', 1.5);
plot(t, val(E_total_co2(3,:)), 'g-^', 'LineWidth', 1.5);
legend('机组1', '机组2', '机组3', 'Location', 'best');
title('CO2捕集量');
xlabel('时间 (h)'); ylabel('捕集量 (kg)');
grid on; box on;
% 溶液体积变化
subplot(2,1,2);
plot(t, val(V_FY(1,:)), 'r--o', 'LineWidth', 1.5); hold on;
plot(t, val(V_PY(1,:)), 'r-s', 'LineWidth', 1.5);
plot(t, val(V_FY(2,:)), 'b--^', 'LineWidth', 1.5);
plot(t, val(V_PY(2,:)), 'b-d', 'LineWidth', 1.5);
legend('机组1富液', '机组1贫液', '机组2富液', '机组2贫液', 'Location', 'best');
title('溶液体积动态变化');
xlabel('时间 (h)'); ylabel('体积 (m³)');
grid on; box on;
else
error('优化失败: %s', result.info);
end
最新发布