已按上述提出的修改意见进行修改，请再次检查修改后的代码：%% 日前调度模型
clc
clear
%% 数据导入
% 调度周期
T = 24;
% 风电预测功率
P_wt_pre = [419.043583980000,374.323177900000,392.119523040000,398.681756480000,460.739421520000,405.595056020000,235.598101760000,222.689675320000,260.448460160000,277.386108580000,183.566952380000,120.085357742000,81.9212982000000,127.760610824000,163.730879200000,195.313565580000,209.128520300000,267.555362200000,374.097494120000,407.577665040000,421.392893460000,468.046794740000,494.700020360000,479.292094960000];
% 电负荷
P_Load_e = 0.74*[1209.12680000000,1391.49080000000,1417.72400000000,1431.18080000000,1398.85200000000,1322.37000000000,1410.94800000000,1492.84800000000,1554.33600000000,1390.15800000000,1770.80400000000,1786.93200000000,1512.25200000000,1305.10800000000,1210.39800000000,1228.51400000000,1440.85200000000,1770.42600000000,1969.12800000000,1685.12400000000,1488.74600000000,1445.37400000000,1286.13800000000,1377.11000000000];
% 热负荷
P_Load_h = 0.72*[358.676650700000,358.415415000000,370.788488900000,369.577305000000,344.641167600000,321.889910800000,352.311998500000,326.544656500000,315.744330800000,265.474819400000,232.010437700000,227.021414800000,296.575979700000,292.170468200000,267.079475800000,273.994086500000,288.173414800000,270.458472200000,227.607740900000,345.550477800000,322.246141400000,342.399073100000,347.968033700000,370.337263500000];
% 分时电价
C_buy_e = 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];
%% 设置参数
P_CHP_max = [150,80]; % CHP机组处理（2台）
P_MT_max = [400,455,200]; % 火电机组最大出力
P_MT_min = [200,120,100]; % 火电机组最小出力
R_u = 4*[50,50,25]; % 火电机组上爬坡
a = [0.00048,0.00031,0.0002]; % 成本系数a
b = [16.2,17.3,16.6]; % 成本系数b
c = [1000,970,700]; % 成本系数c
C_eg = [0.91,0.95,0.98]; % 碳排放系数
E_bata = 0.81; % 碳捕集效率
P_yita = [1.05,1.05,1.05]; % 再生塔和压缩机最大工作状态系数
P_lamda = [0.269,0.269,0.269]; % 单位捕碳量能耗
P_D1 = 0.5*[20,20,15]; % 碳捕集固定能耗
P_D_gd = ones(3,24);
for i=1:3
P_D_gd(i,:)=P_D1(i).*ones(1,24);
end
M_MEA = 61.08; % M_MEA的摩尔质量
M_co2 = 44; % 二氧化碳的摩尔质量
C_sita = 0.4; % 再生塔解析量
C_CR = 30; % 醇胺溶液浓度(%)
C_rou_R = 1.01; % 醇胺溶液密度
V_CR = [2000,2000,2000]; % 溶液储液装置容量
K_R = 1.17; % 乙醇胺溶剂成本系数
C_fai = 1.5; % 溶剂运行损耗系数
% h区间分段
h_min = [200,120,100];
h_max = [401,456,201];
h_gn = 5;
hl_1 = (h_max-h_min)./h_gn;
hl_2 = zeros(3,h_gn+1);
% k区间分段
k_min = [-400,-455,-200];
k_max = [-199,-119,-99];
k_gn = 5;
kl_1 = (k_max-k_min)./k_gn;
kl_2 = zeros(3,k_gn+1);
for i=1:3
hl_2(i,:)=h_min(i):hl_1(i):h_max(i);
end
for i=1:3
kl_2(i,:)=k_min(i):kl_1(i):k_max(i);
end
% 热电比系数
heat_ratio = 0.6 / 0.4; % 燃气轮机热电比
%% 决策变量
% 电力源出力
P_CHP_e = sdpvar(2,24,'full'); % 燃气轮机电出力
P_wt_e = sdpvar(1,24,'full'); % 风电机组出力
P_MT_e = sdpvar(3,24,'full'); % 火电机组出力
P_e_EB = sdpvar(2,24,'full'); % 电锅炉耗电
% 热力源出力
P_CHP_h = sdpvar(2,24,'full'); % 燃气轮机热出力
P_EB_h = sdpvar(2,24,'full'); % 电锅炉热出力
% 天然气输入
P_g_CHP = sdpvar(2,24,'full'); % 天然气需求
% 碳捕集相关
E_G_all = sdpvar(3,24,'full'); % 碳捕集机组产生的总碳排放
E_total_co2 = sdpvar(3,24,'full'); % 机组捕获的总碳排放
E_CG_co2 = sdpvar(3,24,'full'); % 储液装置提供的待捕集二氧化碳量
P_B = sdpvar(3,24,'full'); % 机组运行能耗
P_J = sdpvar(3,24,'full'); % 机组的净出力
V_CA = sdpvar(3,24,'full');
V_FY = sdpvar(3,24,'full'); % 富液体积
V_PY = sdpvar(3,24,'full'); % 贫液体积
P_tra_e = sdpvar(1,24,'full'); % 可转移电负荷
P_cut_e = sdpvar(1,24,'full'); % 可削减电负荷
P_DE_Load_e = sdpvar(1,24,'full'); % 经过过需求响应后的电负荷
H_tra_h = sdpvar(1,24,'full'); % 可转移热负荷
H_cut_h = sdpvar(1,24,'full'); % 可削减热负荷
P_DH_Load_h = sdpvar(1,24,'full'); % 经过过需求响应后的热负荷
% 旋转备用变量
spi_reserve = sdpvar(1,24,'full');
% 正出力分段专用 (三维数组存储)
y1 = sdpvar(3, 24, 'full');
h_w = sdpvar(h_gn+1, 24, 3, 'full'); % 权重变量 [分段点×时间×机组]
h_z = binvar(h_gn, 24, 3, 'full'); % SOS2二进制变量 [分段×时间×机组]
% 负出力分段专用 (三维数组存储)
y2 = sdpvar(3, 24, 'full');
k_w = sdpvar(k_gn+1, 24, 3, 'full'); % 权重变量 [分段点×时间×机组]
k_z = binvar(k_gn, 24, 3, 'full'); % SOS2二进制变量 [分段×时间×机组]
%% 约束条件
Cons = []; % 初始化约束集合
% 1. 机组出力约束
for t = 1:T
% 电力平衡
Cons = [Cons,
sum(P_CHP_e(:, t)) + sum(P_MT_e(:, t)) + P_wt_e(t) == P_DE_Load_e(t) + sum(P_e_EB(:, t))];
% 热力平衡
Cons = [Cons,
sum(P_CHP_h(:, t)) + sum(P_EB_h(:, t)) == P_DH_Load_h(t)];
for i = 1:3 % 三个火电机组
% 碳捕集约束
Cons = [Cons,
E_G_all(i, t) == C_eg(i) * P_MT_e(i, t), % 总碳排放
0 <= E_CG_co2(i, t), % 待捕集CO2非负
0 <= P_B(i, t), % 碳捕集能耗非负
E_total_co2(i, t) == E_CG_co2(i, t) + 0.25 * E_bata * C_eg(i) * (y1(i, t) - y2(i, t)), % 捕获总碳量
0 <= E_total_co2(i, t) <= P_yita(i) * E_bata * C_eg(i) * P_MT_max(i), % 捕获量范围
P_B(i, t) == P_lamda(i) * E_total_co2(i, t), % 捕集能耗关系
P_MT_e(i, t) == P_D_gd(i, t) + P_J(i, t) + P_B(i, t), % 碳捕集电厂功率平衡
P_MT_min(i) <= P_MT_e(i, t) <= P_MT_max(i), % 火电机组出力范围
P_MT_min(i) - P_lamda(i)*P_yita(i)*E_bata*C_eg(i)*P_MT_max(i) - P_D_gd(i,t) <= P_J(i, t) <= P_MT_max(i) - P_D_gd(i, t) % 净出力范围
];
end
% 风电约束
Cons = [Cons,
0 <= P_wt_e(t) <= P_wt_pre(t), % 风电出力范围
sum(P_e_EB(:, t)) + P_wt_e(t) <= P_wt_pre(t) % 弃风约束
];
end
% 2. 旋转备用约束
for t = 1:T
Cons = [Cons,
spi_reserve(t) <= sum(R_u), % 爬坡能力限制
spi_reserve(t) <= sum(P_MT_max) - sum(P_MT_e(:, t)), % 可用容量限制
spi_reserve(t) >= 0.01 * P_DE_Load_e(t), % 最小备用要求
spi_reserve(t) >= 0 % 物理意义要求
];
end
% 3. 分段线性化约束
% 正出力分段 (y1)
for i = 1:3
for t = 1:T
% 二次函数近似
Cons = [Cons, y1(i, t) == (hl_2(i, :).^2) * h_w(:, t, i)];
% SOS2约束
Cons = [Cons,
sum(h_w(:, t, i)) == 1, % 权重和为1
h_w(:, t, i) >= 0, % 权重非负
h_w(1, t, i) <= h_z(1, t, i), % 首端点约束
h_w(end, t, i) <= h_z(end, t, i) % 末端点约束
];
% 中间点约束
for s = 2:h_gn
Cons = [Cons,
h_w(s, t, i) <= h_z(s-1, t, i) + h_z(s, t, i)
];
end
% 二进制变量约束
Cons = [Cons,
sum(h_z(:, t, i)) == 1, % 仅一个区间激活
h_z(:, t, i) >= 0 % 二进制变量非负
];
end
end
% 负出力分段 (y2)
for i = 1:3
for t = 1:T
% 二次函数近似
Cons = [Cons, y2(i, t) == (kl_2(i, :).^2) * k_w(:, t, i)];
% SOS2约束
Cons = [Cons,
sum(k_w(:, t, i)) == 1, % 权重和为1
k_w(:, t, i) >= 0, % 权重非负
k_w(1, t, i) <= k_z(1, t, i), % 首端点约束
k_w(end, t, i) <= k_z(end, t, i) % 末端点约束
];
% 中间点约束
for s = 2:k_gn
Cons = [Cons,
k_w(s, t, i) <= k_z(s-1, t, i) + k_z(s, t, i)
];
end
% 二进制变量约束
Cons = [Cons,
sum(k_z(:, t, i)) == 1, % 仅一个区间激活
k_z(:, t, i) >= 0 % 二进制变量非负
];
end
end
% 4. 爬坡约束
for t = 1:T-1
for i = 1:3 % 火电机组
Cons = [Cons,
-R_u(i) <= P_MT_e(i, t+1) - P_MT_e(i, t) <= R_u(i)];
end
for j = 1:2 % 燃气轮机和电锅炉
Cons = [Cons,
P_CHP_e(j, t+1) - P_CHP_e(j, t) <= 50,
P_e_EB(j, t+1) - P_e_EB(j, t) <= 30
];
end
end
% 5. 燃气系统约束
for t = 1:T
for j = 1:2
Cons = [Cons,
0 <= P_g_CHP(j, t) <= 130, % 供气限制
0 <= P_CHP_e(j, t) <= P_CHP_max(j), % 发电限制
P_CHP_e(j, t) == 0.98 * P_g_CHP(j, t), % 气-电转换
P_CHP_h(j, t) == heat_ratio * P_CHP_e(j, t), % 热电关系
0 <= P_e_EB(j, t) <= P_wt_pre(t) - P_wt_e(t), % 电锅炉限制
P_EB_h(j, t) == 0.98 * P_e_EB(j, t) % 电热转换
];
end
end
% 6. 碳捕集动态约束
V_PY0 = [1000, 1000, 1000]; % 贫液初始值
V_FY0 = [1000, 1000, 1000]; % 富液初始值
for i = 1:3
% 初始条件 (t=1)
Cons = [Cons,
V_CA(i, 1) == (M_MEA * E_CG_co2(i, 1)) / (M_co2 * C_sita * C_CR * C_rou_R),
V_PY(i, 1) == V_PY0(i) + V_CA(i, 1),
V_FY(i, 1) == V_FY0(i) - V_CA(i, 1)
];
% 后续时间点 (t=2 to T)
for t = 2:T
Cons = [Cons,
V_CA(i, t) == (M_MEA * E_CG_co2(i, t)) / (M_co2 * C_sita * C_CR * C_rou_R),
V_FY(i, t) == V_FY(i, t-1) - V_CA(i, t),
V_PY(i, t) == V_PY(i, t-1) + V_CA(i, t)
];
end
% 容量约束 (t=1 to T)
for t = 1:T
Cons = [Cons,
0 <= V_FY(i, t) <= V_CR(i), % 富液容量
0 <= V_PY(i, t) <= V_CR(i) % 贫液容量
];
end
% 周期平衡
Cons = [Cons,
V_PY0(i) == V_PY(i, T),
V_FY0(i) == V_FY(i, T)
];
end
% 7. 需求响应约束
% 电负荷响应
for t = 1:T
if (t >= 11 && t <= 12) || (t >= 18 && t <= 20)
Cons = [Cons, -0.02*P_Load_e(t) <= P_cut_e(t) <= 0];
else
Cons = [Cons, P_cut_e(t) == 0];
end
end
% 热负荷响应
for t = 1:T
if (t >= 1 && t <= 3) || (t >= 11 && t <= 12)
Cons = [Cons, -0.02*P_Load_h(t) <= H_cut_h(t) <= 0];
else
Cons = [Cons, H_cut_h(t) == 0];
end
end
% 响应平衡约束
Cons = [Cons,
-0.08*P_Load_e <= P_tra_e <= 0.08*P_Load_e, % 电转移范围
sum(P_tra_e) == 0, % 电转移守恒
P_DE_Load_e == P_Load_e + P_tra_e + P_cut_e, % 响应后电负荷
-0.08*P_Load_h <= H_tra_h <= 0.08*P_Load_h, % 热转移范围
sum(H_tra_h) == 0, % 热转移守恒
P_DH_Load_h == P_Load_h + H_tra_h + H_cut_h % 响应后热负荷
];
%% 目标函数
F=0;
for t=1:24
F1(t)=0.25*C_buy_e(t)*abs(P_tra_e(t))+0.65*C_buy_e(t)*abs(P_cut_e(t)); % 电能需求响应成本
F2(t)=0.25*abs(H_tra_h(t))+0.65*abs(H_cut_h(t)); % 热能需求响应成本
for i=1:3
F3(i,t)=a(i)*P_MT_line(i,t)+b(i)*P_MT_e(i,t)+c(i); % 火电机组燃料费用
F4(i,t)=0.25*E_bata*C_eg(i)*(y1(i,t)-y2(i,t))+E_CG_co2(i,t);
F5(i,t)=E_G_all(i,t)-0.25*E_bata*C_eg(i)*(y1(i,t)-y2(i,t))-E_CG_co2(i,t)-0.7*P_MT_e(i,t);
end
F6(t)=50*sum(P_g_CHP(:,t)); % 燃气费用
F7(t)=C_buy_e(t)*sum(P_e_EB(:,t)); % 电锅炉运行成本
F8(t)=0.3*C_buy_e(t)*(P_wt_pre(t)-P_wt_e(t)-sum(P_e_EB(:,t))); % 弃风惩罚
end
F=F+sum(F1)+sum(F2)+sum(sum(F3))+1.5*1.17*sum(sum(F4))+12*sum(sum(F5))+sum(F6)+sum(F7)+sum(F8);
%% 求解器配置
ops=sdpsettings('solver','cplex','verbose',2,'usex0',0);
ops.cplex.mip.tolerances.mipgap=1e-6;
result=optimize(Cons,F,ops);
F=value(F);
display(['最优规划结果 : ', num2str(F)]);
最新发布
博客介绍了单位K、M、G、T、P、E、Z，指出它们分别对应10的3次、6次、9次等幂次，属于信息技术领域中单位量级的知识内容。
扫一扫
2万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
