1主要内容
该程序参考文献《梯级水光互补系统最大化可消纳电量期望短期优化调度模型》，构建了以最大化整体可消纳电量期望为目标的梯级水光互补系统短期优化调度模型。此模型以机组作为最小调度单元，对电站约束、机组约束以及电网约束展开精细化建模。通过在电站和时段之间合理调配梯级负荷，深度挖掘梯级水电在电网供电支撑与光伏互补协调方面的双重作用，进而有效提升互补系统的整体消纳水平。
在模型求解环节，运用分段线性逼近、引入0-1整数变量、发电水头离散等一系列线性化方法与建模技巧，对模型中的非线性约束进行处理，将原始模型成功转换为混合整数线性规划问题。随后，在Matlab环境中借助CPLEX工具实现求解，最终实现了对包含4个水电站、15 台机组以及2个光伏群所构建的互补系统的优化调度。（源代码在机组对应关系等方面存在问题，本次进行更正后发布）
目标函数：
以梯级水光互补系统的可消纳电量期望最大为目标，程序未考虑每种场景下的概率，因此未考虑光伏的不确定性，这部分功能不难实现，期待下次将这部分进行完善。
约束条件：
模型实现的约束包括：
电站约束 机组约束 电网约束
上述约束为约束大类，具体约束内容可参考原文献。
线性化处理：
流程示意：
步骤 1：读取基础数据并设置计算条件。包括区间流量、梯级发电计划、光伏预测出力、光伏历史预测与实际出力、分区断面约束、爬坡能力等。
步骤 2：模型转换处理。采用 2.1 节所述模型转换方法，对非线性约束进行线性化处理。
步骤 3：光伏出力场景构建。根据计划日光伏预测出力以及 2.2 节所述方法构建光伏出力场景。
步骤 4：模型求解。将目标函数与转化后的约束结合构成的 MILP 模型，在Matlab环境中，调用 CPLEX 求解类，实现模型求解。
步骤 5：结果输出。输出互补系统整体可消纳电量期望值，不同组合场景下的电站出力、机组出力、机组开停机、出库流量、水库水位等结果信息。
2部分代码
%% 清除内存空间
clc
clear
close all
warning off
yalmip('clear')
%% 参数设置
S = 1; % 光伏出力的场景数
I = 4; % 水电站数目
J = 2; % 光伏电厂数目
G = 4; % 约束断面数
T = 96; % 调度期总时段数
dt = 0.15; % 调度时间间隔
tau = 1; % 电站与其上游电站的水流滞时
Z_up_max = [1140 970 837 760]; % 坝前水位上限
Z_up_min = [1071 936 814 709]; % 坝前水位下限
Z_up_begin = [1076 950 822 720]; % 调度期初始水位
Z_up_end = [1076 950 822 720]; % 调度期末控制水位
dZ = [0.5 0.1 0.3 0.2]; % 允许的调度期末水位偏差
Q_max = [3866 11142 15956 18360]; % 出库流量上限
Q_min = [866 1142 5956 8360]; % 出库流量下限
P_hydro_max = [600 695 600 1250]; % 水电站的出力上限
P_hydro_min = [100 125 100 150]; % 水电站的出力下限
Ni = [3 4 3 5]; % 电站i所包含的机组总数
Y_i = ... % 各水电站包含的机组集合
{[1 2 3],[4 5 6 7],[8 9 10],[11 12 13 14 15]};
N = sum(Ni); % 该区域机组总数
P_in_max = ... % 第n台机组的出力上限
[200 200 200 190 190 190 125 200 200 200 250 250 250 250 250];
P_in_min = ... % 第n台机组的出力下限
[45 45 45 40 40 40 30 0 0 0 50 50 50 70 70];
Q_p_int_max = ... % 电站i的发电流量上限
[490.5 490.5 490.5 632.2 632.2 632.2 632.2 994.5 994.5 994.5 1250 1250 1250 1250 1250];
Q_p_int_min = 0; % 电站i的发电流量下限
P_1_3nk_max = ... % 电站1-3的机组第k个振动区的出力上限
[90 90 90 110 110 110 80 130 130 130];
P_1_3nk_min = ... % 电站1-3的机组第k个振动区的出力下限
[45 45 45 40 40 40 30 0 0 0];
P_4nk_max = ... % 电站4的机组第k个振动区的出力上限
[70 70 70 80 80;90 90 90 170 170];
P_4nk_min = ... % 电站4的机组第k个振动区的出力下限
[50 50 50 70 70;80 80 80 140 140];
T_on_in = 8; % 电站 i 的第 n 台机组的最小开机持续时段数
T_off_in = 8; % 电站 i 的第 n 台机组的最小停机持续时段数
M_on_in = 8; % 电站 i 的第 n 台机组的最大开机次数
dP_in = 60; % 电站 i的第 n 台机组的爬坡能力为60MW/15 min
te = 4; % 机组在一轮出力升降过程中需持续的最少时段数
% ai = ... % 电站 i 的水头损失系数
% [9 9 9 7 7 7 7 8 8 8 4 4 4 4 4]*1e-8;
% bi = ... % 电站 i 的水头损失常数
% [12 12 12 6 6 6 6 15 15 15 9 9 9 9 9]*1e-5;
ai = ... % 电站 i 的水头损失系数
[9 8 4]*1e-8;
bi = ... % 电站 i 的水头损失常数
[12 15 9]*1e-5;
P_plan = ... % 梯级水电发电计划
[2200 2200 2200 2200 2100 2100 2100 2100 2050 2050 2050 2050 2000 2000 1995 1995 ...
1995 1900 1900 1900 1995 1995 2100 2100 2205 2205 2310 2310 2415 2520 2625 2650 ...
2700 2700 2700 2750 2750 2850 2850 2850 2750 2750 2750 2750 2660 2660 2660 2570 ...
2570 2570 2565 2565 2565 2565 2660 2660 2660 2660 2755 2755 2755 2755 2780 2780 ...
2780 2780 2800 2800 2755 2755 2755 2755 2755 2800 2800 2800 2800 2850 2850 2850 ...
2850 2755 2755 2700 2700 2700 2700 2600 2600 2600 2520 2520 2520 2415 2415 2415];
e = 0.02; % 允许偏差
L_gt = [900 2000 800 600]; % 第g个约束断面的负荷容量上限
C_gt = [900 1500 750 500]; % 第g个约束断面的输电容量上限
M_ab = 12; % 功率调整（向上和向下）时段数上限
N_pv = [550 1200]; % 光伏场 j 的装机容量
P_pv0 = ... % 光伏场j的预测出力
[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 14 31 68 98 139 175 126 209 229 247 260 274 288 304 309 309 293 288 329 332 325 337 337 299 277 302 206 189 237 164 135 87 67 96 89 44 56 88 40 62 39 9 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 25 30 73 84 99 355 446 209 524 313 358 382 521 436 775 744 825 658 706 803 797 703 492 458 526 329 439 408 498 258 337 139 177 164 190 113 73 62 80 41 171 69 23 5 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0];
P_pv = ... % 光伏场j的实际出力
[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 15 32 74 94 133 159 119 203 207 268 251 260 304 330 304 315 267 265 321 346 304 331 317 305 299 275 227 193 229 151 130 88 61 96 84 47 61 83 40 68 39 9 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 26 32 75 83 102 364 490 213 547 286 328 384 525 438 720 711 749 635 659 729 820 641 488 492 476 326 404 438 451 256 342 138 162 157 192 122 72 63 85 39 159 76 24 5 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ];
dP_pv = abs(P_pv - P_pv0); % 光伏预测偏差
tau_int = 0.005; % 出力曲线的系数
%% 决策变量
P_hydro = sdpvar(I,T); % 水电站i在时段t的出力
P_w = sdpvar(G,T); % 第g个约束断面在时段 t 的弃电量
V_it = sdpvar(I,T); % 水电站 i 在时段 t 末的库容
I_it = sdpvar(I,T); % 水电站i在时段t的入库流量
Q_it = sdpvar(I,T); % 电站 i 在时段 t的出库流量
R_it = sdpvar(I,T); % 电站i - 1和电站i之间在时段t的区间流量
Q_p_it = sdpvar(I,T); % 水电站i在时段t的发电流量
Q_d_it = sdpvar(I,T); % 电站 i 在时段 t 的弃水流量
Z_up_it = sdpvar(I,T); % 电站i所在水库在时段t的坝前水位
Z_down_it = sdpvar(I,T); % 水电站i在时段t的入库流量
P_int = sdpvar(N,T); % 第n台机组在时段t的出力
u_int = binvar(N,T); % 第n台机组在时段t的开停机状态
Q_p_int = sdpvar(N,T); % 第n台机组时段 t 的发电流量
y_on_int = binvar(N,T); % 第n台机组在时段t的启动操作变量
y_off_int = binvar(N,T); % 第n台机组在时段t的停机操作变量
H_int = sdpvar(N,T); % 第 n 台机组在时段 t的发电水头
H_loss_int = sdpvar(N,T); % 第 n 台机组在时段 t的水头损失
P_sgt = sdpvar(G,T); % 约束断面g在t时段所有机组总出力
theta_1_3ntk = binvar(10,T,2); % 水电站1-3的机组振动区指示变量
theta_4ntk = binvar(5,T,3); % 水电站4的机组振动区指示变量
alpha_int = binvar(N,T); % 1表示时段 t + 1 功率向下调节
beta_int = binvar(N,T); % 1表示时段 t + 1 功率向上调节
%% 目标函数
F = sum(sum(P_hydro) + sum(P_pv) - sum(P_w));
%% 约束条件
cons = [];
for t = 2:T
% 式2
cons = [cons , V_it(:,t) == V_it(:,t-1) + 3600*(I_it(:,t) - Q_it(:,t))*dt];
cons = [cons , abs(Z_up_it(:,t) - Z_up_it(:,t-1)) <= [0.005;0.02;0.05;0.01]];
cons = [cons , abs(V_it(:,t) - V_it(:,t-1)) <= 0.5];
% 式3
cons = [cons , I_it(1 , t) == R_it(1 , t)];
for i = 2:I
cons = [cons , I_it(i , t) == Q_it(i - 1 , t - 1) + R_it(i , t)];
end
end
cons = [cons , V_it(:,1) == [35;7.2;22;18]];
% 式4
cons = [cons , Q_it == Q_p_it + Q_d_it];
% 式5
cons = [cons , Z_up_max'*ones(1,T) >= Z_up_it >= Z_up_min'*ones(1,T)];
% 式6
cons = [cons , Z_up_it(:,1) == Z_up_begin'];
3程序结果
下载链接见发布者备注信息！
