【更正版】梯级水光互补系统最大化可消纳电量期望短期优化调度模型

目录

1 主要内容

目标函数：

约束条件：

线性化处理：

流程示意：

2 部分代码

3 程序结果

4 下载链接

---

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   00   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% 式2cons = [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];% 式3cons = [cons , I_it(1 , t) == R_it(1 , t)];for i = 2:Icons = [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 程序结果

• 4 下载链接