再次优化以下代码,重构后输出:%% 日前调度优化模型
clc; clear;
%% 参数初始化
T = 24; % 调度周期
% 预测数据
P_prew = 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];
EleLoad = 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];
HeatLoad = 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; % 火电机组分段数
% 初始化分段边界
[h_min, h_max] = deal([200, 120, 100], [401, 456, 201]);
[k_min, k_max] = deal([-400, -455, -200], [-199, -119, -99]);
% 计算分段边界
hl1 = (h_max - h_min) / h_gn;
kl1 = (k_max - k_min) / k_gn;
gl1 = (P_G_max - P_G_min) / gn;
hl2 = arrayfun(@(min, max, step) min:step:max, h_min, h_max, hl1, 'UniformOutput', false);
kl2 = arrayfun(@(min, max, step) min:step:max, k_min, k_max, kl1, 'UniformOutput', false);
gl2 = arrayfun(@(min, max, step) min:step:max, P_G_min, P_G_max, gl1, 'UniformOutput', false);
%% 决策变量定义
% 电力相关
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'); % 响应后热负荷
% 分段线性化变量
% 火电机组分段
w = sdpvar(gn+1, T, 3); % 三维数组: [分段点, 时间, 机组]
z = binvar(gn, T, 3); % 三维二进制变量
P_G_line = sdpvar(3, T, 'full');
% 上/下分段变量
h_w = sdpvar(h_gn+1, T, 3);
h_z = binvar(h_gn, T, 3);
k_w = sdpvar(k_gn+1, T, 3);
k_z = binvar(k_gn, T, 3);
y1 = sdpvar(3, T, 'full');
y2 = sdpvar(3, T, 'full');
%% 约束构建
C = [];
% ===== 基础约束 =====
for t = 1:T
% 风电出力约束
C = [C, 0 <= P_w(t) <= P_prew(t)];
% 电锅炉约束
C = [C, sum(EB(:, t)) + P_w(t) <= P_prew(t)];
for i = 1:3
% 火电机组出力约束
C = [C, P_G_min(i) <= P_G(i, t) <= P_G_max(i)];
% 碳捕集系统约束
C = [C, E_G(i, t) == eg(i) * P_G(i, t)]; % 总碳排放
C = [C, E_total_co2(i, t) == E_CG(i, t) + 0.25 * E_bata * eg(i) * (y1(i, t) - y2(i, t))]; % 捕获总量
C = [C, 0 <= E_total_co2(i, t) <= P_yita(i) * E_bata * eg(i) * P_G_max(i)];
C = [C, P_B(i, t) == P_lamda(i) * E_total_co2(i, t)]; % 运行能耗
C = [C, P_G(i, t) == P_J(i, t) + P_D(i, t) + P_B(i, t)]; % 总功率平衡
C = [C, P_G_min(i) - P_lamda(i) * P_yita(i) * E_bata * eg(i) * P_G_max(i) - P_D(i, t) <= P_J(i, t) <= ...
P_G_max(i) - P_D(i, t)]; % 净出力范围
end
end
% 旋转备用约束
C = [C, min(sum(R_u), sum(P_G_max) - sum(P_G)) >= 0.01 * max(P_DE)];
% ===== 分段线性化约束 =====
% 通用分段约束函数
function C = add_segment_constraints(C, w, z, values, var, T, n_segments)
n_units = size(var, 1);
for i = 1:n_units
C = [C, var(i, :) == values{i} * squeeze(w(:, :, i))];
for t = 1:T
% 凸组合约束
C = [C, sum(w(:, t, i)) == 1, w(:, t, i) >= 0];
% 相邻约束
C = [C, w(1, t, i) <= z(1, t, i)];
for seg = 2:n_segments
C = [C, w(seg, t, i) <= z(seg-1, t, i) + z(seg, t, i)];
end
C = [C, w(n_segments+1, t, i) <= z(n_segments, t, i)];
% 激活一个分段
C = [C, sum(z(:, t, i)) == 1];
end
end
end
% 应用分段约束
C = add_segment_constraints(C, w, z, gl2, P_G_line, T, gn);
C = add_segment_constraints(C, h_w, h_z, hl2, y1, T, h_gn);
C = add_segment_constraints(C, k_w, k_z, kl2, y2, T, k_gn);
% 连接变量
for i = 1:3
C = [C, P_G(i, :) == gl2{i} * squeeze(w(:, :, i))];
C = [C, P_G_line(i, :) == (gl2{i}.^2) * squeeze(w(:, :, i))];
end
C = [C, y1 >= y2]; % 上下界关系约束
% ===== 机组运行约束 =====
% 爬坡约束
for t = 1:T-1
for i = 1:3
C = [C, -R_u(i) <= P_G(i, t+1) - P_G(i, t) <= R_u(i)];
end
for j = 1:2
C = [C, -50 <= GT_P(j, t+1) - GT_P(j, t) <= 50];
C = [C, -30 <= EB(j, t+1) - EB(j, t) <= 30];
end
end
% 燃气轮机约束
for t = 1:T
for j = 1:2
C = [C, 0 <= P_gas(j, t) <= 130];
C = [C, 0 <= GT_P(j, t) <= GT_P_max(j)];
C = [C, GT_P(j, t) == 0.98 * P_gas(j, t)];
C = [C, GT_H(j, t) == 0.6 * (GT_P(j, t) / 0.4)];
C = [C, 0 <= EB(j, t) <= P_prew(t) - P_w(t)];
C = [C, EB_H(j, t) == 0.98 * EB(j, t)];
end
end
% ===== 碳捕集系统动态约束 =====
V_PY0 = [1000, 1000, 1000]; % 贫液初始值
V_FY0 = [1000, 1000, 1000]; % 富液初始值
for i = 1:3
C = [C, V_PY0(i) == V_PY(i, T)]; % 终值平衡
C = [C, V_FY0(i) == V_FY(i, T)]; % 终值平衡
for t = 1:T
C = [C, V_CA(i, t) == (M_MEA * E_CG(i, t)) / (M_co2 * sita * CR * rou_R)]; % 溶液体积
C = [C, 0 <= V_FY(i, t) <= V_CR(i)];
C = [C, 0 <= V_PY(i, t) <= V_CR(i)];
C = [C, 0 <= (0.25 * eg(i) * (y1(i, t) - y2(i, t))) <= P_G(i, t)];
end
for t = 1:T-1
C = [C, V_FY(i, t+1) == V_FY(i, t) - V_CA(i, t+1)]; % 富液动态
C = [C, V_PY(i, t+1) == V_PY(i, t) + V_CA(i, t+1)]; % 贫液动态
end
end
% ===== 需求响应约束 =====
% 电负荷响应
for t = 1:T
if (t >= 11 && t <= 12) || (t >= 18 && t <= 20)
C = [C, -0.02 * EleLoad(t) <= P_cut(t) <= 0];
else
C = [C, P_cut(t) == 0];
end
end
% 热负荷响应
for t = 1:T
if (t >= 0 && t <= 3) || (t >= 11 && t <= 12)
C = [C, -0.02 * HeatLoad(t) <= H_cut(t) <= 0];
else
C = [C, H_cut(t) == 0];
end
end
% 响应平衡约束
C = [C, -0.08 * EleLoad <= P_tran <= 0.08 * EleLoad];
C = [C, sum(P_tran) == 0];
C = [C, P_DE == EleLoad + P_tran + P_cut];
C = [C, -0.08 * HeatLoad <= H_tran <= 0.08 * HeatLoad];
C = [C, sum(H_tran) == 0];
C = [C, H_DE == HeatLoad + H_tran + H_cut];
% ===== 系统平衡约束 =====
C = [C, P_w + sum(GT_P, 1) + sum(P_J, 1) == P_DE]; % 电力平衡
C = [C, sum(EB_H, 1) + sum(GT_H, 1) == H_DE]; % 热力平衡
%% 目标函数构建
% 组件成本计算
cost_DR_elec = 0.25 * gama_elec .* abs(P_tran) + 0.65 * gama_elec .* abs(P_cut);
cost_DR_heat = 0.25 * abs(H_tran) + 0.65 * abs(H_cut);
cost_gen = zeros(1, T);
cost_carbon = zeros(1, T);
cost_emission = zeros(1, T);
for i = 1:3
cost_gen = cost_gen + a(i) * P_G_line(i, :) + b(i) * P_G(i, :) + c(i);
cost_carbon = cost_carbon + 1.5 * K_R * (0.25 * E_bata * eg(i) * (y1(i, :) - y2(i, :)) + E_CG(i, :));
cost_emission = cost_emission + 12 * (E_G(i, :) - 0.25 * E_bata * eg(i) * (y1(i, :) - y2(i, :)) - E_CG(i, :) - 0.7 * P_G(i, :));
end
cost_gas = 50 * sum(P_gas, 1);
cost_EB = gama_elec .* sum(EB, 1);
cost_wind_curtail = 0.3 * gama_elec .* (P_prew - P_w - sum(EB, 1));
% 总目标函数
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(C, F, ops);
if result.problem == 0
disp(['最优值: ', num2str(value(F))]);
% 结果提取
[GT_P_val, P_w_val, P_G_val, EB_val] = deal(value(GT_P), value(P_w), value(P_G), value(EB));
[GT_H_val, EB_H_val, P_gas_val] = deal(value(GT_H), value(EB_H), value(P_gas));
[E_G_val, E_total_co2_val, E_CG_val] = deal(value(E_G), value(E_total_co2), value(E_CG));
[P_B_val, P_J_val, V_CA_val] = deal(value(P_B), value(P_J), value(V_CA));
[y1_val, y2_val] = deal(value(y1), value(y2));
[P_tran_val, P_cut_val, P_DE_val] = deal(value(P_tran), value(P_cut), value(P_DE));
[H_tran_val, H_cut_val, H_DE_val] = deal(value(H_tran), value(H_cut), value(H_DE));
%% 结果可视化
% 原始数据图
figure('Name', '原始数据');
plot(1:T, P_prew, 'r-o', 'LineWidth', 1.5, 'DisplayName', '风电出力');
hold on;
plot(1:T, EleLoad, 'g-*', 'LineWidth', 1.5, 'DisplayName', '电负荷');
plot(1:T, HeatLoad, 'b-^', 'LineWidth', 1.5, 'DisplayName', '热负荷');
legend show; xlabel('时间/h'); ylabel('功率/kW'); title('原始数据'); grid on;
% 电力平衡图
Plot_EleNet = [P_J_val(1, :); P_J_val(2, :); P_J_val(3, :);
P_w_val + sum(EB_val, 1);
GT_P_val(1, :); GT_P_val(2, :); sum(EB_val, 1)]';
figure('Name', '电功率平衡');
bar(1:T, Plot_EleNet, 'stacked');
hold on;
plot(1:T, EleLoad, 'k-', 'LineWidth', 2, 'DisplayName', '原始电负荷');
plot(1:T, P_DE_val, 'r-o', 'LineWidth', 1.5, 'DisplayName', '响应后电负荷');
legend([{'火电1','火电2','火电3','风电','燃气轮机1','燃气轮机2','电锅炉'}, {'原始电负荷','响应后电负荷'}]);
xlabel('时间/h'); ylabel('功率/kW'); title('电功率平衡'); grid on;
% 热力平衡图
Plot_HeatNet = [sum(GT_H_val, 1); sum(EB_H_val, 1)]';
figure('Name', '热功率平衡');
bar(1:T, Plot_HeatNet, 'stacked');
hold on;
plot(1:T, HeatLoad, 'k-', 'LineWidth', 2, 'DisplayName', '原始热负荷');
plot(1:T, H_DE_val, 'r-o', 'LineWidth', 1.5, 'DisplayName', '响应后热负荷');
legend('燃气轮机', '电锅炉', '原始热负荷', '响应后热负荷');
xlabel('时间/h'); ylabel('功率/kW'); title('热功率平衡'); grid on;
% 风电出力图
figure('Name', '风电功率');
plot(1:T, P_w_val + sum(EB_val, 1), 'g-o', 'LineWidth', 1.5, 'DisplayName', '实际利用');
hold on;
plot(1:T, P_prew, 'k-', 'LineWidth', 1.5, 'DisplayName', '预测出力');
legend show; xlabel('时间/h'); ylabel('功率/kW'); title('风电功率'); grid on;
% 碳捕集能耗图
figure('Name', '碳捕集能耗');
plot(1:T, sum(P_B_val, 1), 'g-^', 'LineWidth', 1.5);
xlabel('时间/h'); ylabel('功率/kW'); title('碳捕集能耗'); grid on;
% 火电机组净出力
figure('Name', '火电机组净出力');
plot(1:T, sum(P_J_val, 1), 'r-o', 'LineWidth', 1.5);
xlabel('时间/h'); ylabel('功率/kW'); title('火电机组净出力'); grid on;
% 二氧化碳捕集量
figure('Name', '碳捕集量');
plot(1:T, E_total_co2_val(1, :), 'r-', 'LineWidth', 1.5, 'DisplayName', '机组1');
hold on;
plot(1:T, E_total_co2_val(2, :), 'g-^', 'LineWidth', 1.5, 'DisplayName', '机组2');
plot(1:T, E_total_co2_val(3, :), 'b-*', 'LineWidth', 1.5, 'DisplayName', '机组3');
legend show; xlabel('时间/h'); ylabel('碳捕集量'); title('机组二氧化碳捕集量'); grid on;
else
error('优化失败: %s', result.info);
end