梯级水光互补调度模型复现：基于随机优化的可消纳电量期望最大化

研究能源调度、做水库优化、搞新能源并网的人，看到“EI复现”加上“梯级水光互补”“最大化可消纳电量期望”这一串关键词，基本就知道说的是哪类问题了：一条流域上的若干个梯级水电站，搭配光伏电站做日前短期优化调度，目标是让整个系统在光伏出力不确定的情况下，上网电量的期望值尽可能大。我最近花了接近两周时间，把这个模型用Python完整复现了一遍，从目标函数、约束条件到场景生成、求解验证全部跑通，过程中踩了不少坑，也把论文里略过的一些关键细节补全了。这篇博文就把我的完整思路和代码实现过程写下来，给正在复现水光互补调度、梯级水电短期优化或者随机优化模型的朋友做个参考。

1. 梯级水光互补到底补的是什么：先搞懂工程问题再谈建模

1.1 水电站的调节能力决定互补上限

光伏出力的特点，做新能源的人都很熟悉：白天中午大发，傍晚归零，夜间完全没有出力，而且受云层影响波动剧烈。没有水库配合的情况下，电网为了接纳光伏，要么在午间被迫弃光，要么让火电机组快速压出力来腾空间。水电不一样，水库本身就是天然的储能设施，可以把白天的水先存住，等到晚上或者负荷高峰再放出来发电。

但单座水电站和梯级水电站的互补能力完全不是一个量级。单站调节，只能决定“我自己这个水库什么时候放水”；梯级电站则涉及更复杂的耦合：上游电站放出来的水，经过一段河道的流动时间，会变成下游电站的入库流量，下游电站可以再发一次电。也就是说，上游发过电的水，到下游还能产生额外的电量收益。这是梯级比单站多出来的“隐藏资源”，也是建模时必须考虑上下游水力联系的直接原因。

另外还有一个经常被忽略的地方：水流从上游到下游存在时间延迟。如果上游水库在t时段放水，这批水可能要到t+τ时段才进入下游水库。τ的大小取决于河道长度和流速，可能是0.5小时，也可能是几个小时。复现论文模型时如果忽略这个时滞，水量平衡就错了，下游电站的出力过程也会整体偏移，最后算出来的消纳电量期望值看着很高，实际上根本不可能实现。

1.2 “最大化可消纳电量期望”到底在优化什么

“可消纳电量”不是“发电量”，这是很多初学者第一个理解偏差。发电量只关心机组发出来多少电，可消纳电量关心的是电网实际接纳了多少电。一个电站发出来的电如果没有负荷、没有外送通道消纳，那这部分电量就是弃电，不但不能算作收益，反而可能带来调峰压力。

所以模型里必然存在一个电网消纳能力约束：每个时段整个系统的上网功率不能超过电网允许的上限。这个上限可能来自负荷曲线，也可能来自联络线传输容量。在这个约束之下，目标函数要最大化的是水电和光伏加起来真正上网的电量期望值。

再来理解“期望”。光伏明天的出力是多少，今天不可能完全确定，只能基于预测给出一个概率分布。传统的确定性调度只拿一条预测曲线来优化，实际出力一旦偏离预测，很可能出现两种情况：一种是实际出力小于预测，系统浪费了备用容量；另一种是实际出力大于预测，超出消纳边界被迫紧急弃光。随机期望模型用多个可能场景来描述未来的不确定性分布，每个场景是一个可能的光伏出力序列，带一个发生概率，目标函数变成所有场景下上网电量的概率加权和。这样做出来的调度计划，不再只是对一条曲线最优，而是对一簇可能情况整体最优。

1.3 短期调度管的是什么时间尺度

题目标注了“短期优化调度”，这个“短期”一般指日前调度，常见的是把一天分成96个时段（每15分钟一个点）或者24个时段（每小时一个点）。为什么要做这么细的时间粒度？因为水光互补的主要矛盾就发生在一天之内：光伏正午大发、早晚为零，负荷早晚有高峰，水电需要在一天内多次调整出力来配合。如果时间步长太粗，比如按天或者按周调度，光伏出力的日变化特征就完全被磨平了，互补效果根本体现不出来。

模型的时间尺度还决定了要使用什么数据。径流预测、光伏出力预测、负荷预测都是24小时级别的；水库的起调水位和末水位则是从更长期的调度方案中传递下来的边界条件。短期调度模型的输出，是各电站各时段的出力计划、发电流量计划和弃电计划，这些计划最终要下发给现场执行。

2. 目标函数与约束的拆解：可消纳电量期望怎么算出来

2.1 目标函数的数学表达

用最简洁的方式描述这个问题，就是把所有场景、所有时段的上网电量加起来，按场景概率加权：

max Σ_{s∈Ω} π_s · Σ_t Δt · ( Σ_i P^H_{i,t,s} + Σ_j P^PV_{j,t,s} )

变量含义如下：

• Ω是场景集合，s是某个光伏出力场景；
• π_s是场景s的概率，所有概率加起来等于1；
• Δt是时段长度，功率乘时间就是电量；
• P^H_{i,t,s}是第i个水电站t时段s场景的上网出力；
• P^PV_{j,t,s}是第j个光伏电站在t时段s场景的上网出力。

注意“上网出力”和“实际消纳”在这里是同一个东西，因为凡是不能上网的部分都已经作为弃电被排除了。模型只优化能被电网消纳的电量，目标值本身就是“可消纳电量期望”。

2.2 约束条件：水量平衡、出力限制和消纳边界

第一组核心约束是水量平衡，这是梯级水电站模型的灵魂：

V_{i,t+1,s} = V_{i,t,s} + ( I_{i,t,s} + Q^up_{i-1,t-τ,s} - Q^tur_{i,t,s} - Q^spill_{i,t,s} ) · Δt

每一项的含义：

• V_{i,t,s}是水库i在t时段s场景的库容；
• I_{i,t,s}是水库i本地的天然来水（区间径流）；
• Q^up_{i-1,t-τ,s}是上游水库经过τ时段水流延迟后到达本水库的出库流量；
• Q^tur_{i,t,s}是本水库的发电流量；
• Q^spill_{i,t,s}是本水库的弃水流量。

这个约束的本质是水量守恒：库容变化等于入库减出库。如果上游出库流量没有经过时滞就直接计入下游入库，梯级电站之间就没有真正“梯级”起来，模型退化成几个独立水库的并联调度。

第二组约束是水电出力与发电水量的关系。物理上，水电出力等于水的重力势能转化：P = η · ρ · g · Q · H，其中H是发电水头。水头会随库容变化，所以这个关系是非线性的。论文里常见的处理有三种：一是用平均水头或者额定水头近似；二是把出力曲线按流量分段线性化；三是在模型里引入水头动态方程做迭代逼近。复现时你需要先看清楚论文用哪一种，如果论文没有明确说明，工程上最常用的是分段线性化，Gurobi的addGenConstrPWL可以直接完成。

第三组约束是各类运行边界：

• 库容上下限：V_min ≤ V_{i,t,s} ≤ V_max，下限通常由死水位决定，上限由正常蓄水位决定；
• 发电流量上限：0 ≤ Q^tur_{i,t,s} ≤ Q_max，受机组过流能力限制；
• 生态流量要求：下泄流量不能小于某个最小值，这个不能省，否则模型会为了发电把河道放干；
• 光伏消纳上界：0 ≤ P^PV_{j,t,s} ≤ P^PV_max_{j,t,s}，P^PV_max是场景s下光伏预测出力，模型可以选择弃掉一部分光伏；
• 电网消纳能力约束：Σ_i P^H + Σ_j P^PV ≤ C_t，C_t是时段t的系统消纳上限。

2.3 为什么期望目标不能简化成平均场景优化

我刚开始复现时想过一个偷懒的办法：把各场景的光伏出力按概率加权平均成一条曲线，然后跑确定性模型。这样确实可以省掉场景维度，求解速度飞快。但结果和完整随机模型对不上，原因主要有两个。

第一，库容约束是对每个场景分别成立的，也就是说“无论明天实际来多少光伏，我制定的一组放水计划都必须让水库不越限”。如果用平均场景优化，单独看光伏高的场景，水库可能被迫大量弃水或者突破出力上限；单独看光伏低的场景，水电出力可能不够，消纳电量偏低。期望模型在优化时自动权衡了所有场景的计划可行性，平均场景模型做不到这一点。

第二，弃电决策是分场景的。光伏大发场景下可能弃光，光伏小发场景下可能不弃，这种与场景相关的决策直接进入目标函数；如果先平均再优化，相当于把高发场景的弃光压力和平发场景的出力空间强行揉在一起，最终调度计划在现场根本执行不了。

3. Python代码落地的关键环节：变量定义、约束实现与求解框架

3.1 为什么选Gurobi而不是自己写优化算法

复现这类论文模型，第一选择就是用Gurobi这样的商用求解器，配合Python接口。原因很现实：模型变量数量动辄几万个、约束几万个，自己写梯度投影、拉格朗日松弛或者启发式算法，光调参就能耗掉两周，而且无法保证全局最优。Gurobi学术版可以免费申请license，工业界也能接受这个工具链。

有的读者可能会问：用Pyomo框架不行吗？Pyomo本身是建模语言，底层还是调Gurobi或者开源求解器。对小模型没关系，但Pyomo的建模开销在超大规模场景下会拖慢求解速度；直接用gurobipy手写约束，反而更直观可控。复现论文时，我建议直接gurobipy，思路最顺。

环境方面没有特殊要求，Python 3.8以上版本，安装numpy、pandas、gurobipy、matplotlib就够了。如果还没装好Python环境和Gurobi，先把基础环境搞定再回来看后面的代码。

3.2 整体代码结构怎么设计

我的目录结构分五个模块，职责非常清晰：

• load_data.py：读取电站参数、径流序列、光伏场景、电网边界等输入数据；
• generate_scenarios.py：生成光伏出力场景并计算概率，必要时做场景削减；
• build_model.py：核心建模，定义变量、目标函数和全部约束；
• run_case.py：主程序，调用上面三个模块求解并保存结果；
• analyze_results.py：解析输出，绘制出力曲线、库容过程、弃电统计表。

这种结构的好处是：换数据不用改模型，换模型不用动数据读取；做算例对比的时候，只要改run_case.py里的参数组合就行。复现论文时你会反复调试，模块分离绝对是加分项。

3.3 核心建模代码示例

下面给出模型构建的核心代码骨架，以一个小算例为例：2级梯级水电站，1座光伏电站，24个时段，20个场景。变量规模大概在几千个，Gurobi几秒钟就能求解。

import numpy as np from gurobipy import Model, GRB, quicksum def build_model(data, scenarios): m = Model("cascade_hydro_pv") I = range(data["n_hydro"]) # 水电站索引 J = range(data["n_pv"]) # 光伏场索引 T = range(data["T"]) # 时段索引 S = range(len(scenarios["prob"])) # 场景索引 dt = data["dt"] # 时段长度（小时） # 决策变量 V = m.addVars(I, T, S, lb=data["v_min"], ub=data["v_max"], name="Volume") Q_tur = m.addVars(I, T, S, lb=0, name="TurbineFlow") Q_spill = m.addVars(I, T, S, lb=0, name="SpillFlow") P_h = m.addVars(I, T, S, lb=0, name="HydroPower") P_pv = m.addVars(J, T, S, lb=0, name="PVUsed") # 目标函数：最大化可消纳电量期望 objective = quicksum( scenarios["prob"][s] * dt * ( quicksum(P_h[i, t, s] for i in I) + quicksum(P_pv[j, t, s] for j in J) ) for s in S for t in T ) m.setObjective(objective, GRB.MAXIMIZE) # 光伏出力上界：每个场景的预测出力 for s in S: for j in J: for t in T: m.addConstr(P_pv[j, t, s] <= scenarios["pv_max"][s][j][t]) # 电网消纳能力约束 for s in S: for t in T: m.addConstr( quicksum(P_h[i, t, s] for i in I) + quicksum(P_pv[j, t, s] for j in J) <= data["grid_limit"][t] ) # 水量平衡约束（含梯级时滞） tau = data["flow_delay"] # 上游到下游的时滞时段数 for i in I: for s in S: for t in list(T)[1:]: inflow = data["inflow"][i][t][s] if i > 0: # +上游出库（发电流量+弃水）经过时滞到达本库 t_up = t - tau[i] if t_up >= 0: inflow += Q_tur[i - 1, t_up, s] + Q_spill[i - 1, t_up, s] m.addConstr( V[i, t, s] == V[i, t - 1, s] + dt * ( inflow - Q_tur[i, t, s] - Q_spill[i, t, s] ) ) # 水电出力上限：简化线性出力系数K for i in I: for s in S: for t in T: m.addConstr(P_h[i, t, s] <= data["p_coef"][i] * Q_tur[i, t, s]) m.addConstr(P_h[i, t, s] <= data["p_cap"][i]) return m, V, Q_tur, Q_spill, P_h, P_pv

这段代码已经把目标函数、水量平衡、光伏上界、消纳约束和出力上限都搭起来了。实际复现时，如果论文模型里水电出力用的是水头动态方程，就把p_coef * Q_tur这一段替换成分段线性化约束。

Gurobi的分段线性化可以这样写：

# flow_points和power_points是流量-出力分段点坐标 for i in I: for s in S: for t in T: m.addGenConstrPWL( Q_tur[i, t, s], P_h[i, t, s], flow_points[i], power_points[i], name="pwl_hydro")

addGenConstrPWL会自动把非线性函数线性化，求解器会引入辅助变量，不影响整体模型结构。

4. 场景生成、削减与数据口径：模型能不能跑出结果看这一步

4.1 输入数据清单

复现这个模型，最花时间的不是代码本身，而是数据准备。需要的输入数据可以分为四类：

如果是从网上找的开源数据，经常需要花很多时间做单位换算和口径统一。比如径流单位可能是m³/s，也可能直接给m³；库容单位可能是万m³也可能是亿m³；时间粒度可能是15分钟也可能是1小时，这些不统一，模型一跑就是数值灾难。

4.2 光伏场景生成：从历史曲线到概率场景

光伏场景生成最常用的方法是基于历史出力数据聚类。把过去一段时间内每天96点的出力曲线作为样本，用K-means聚成K类，每个类的中心曲线作为一个典型场景，类内样本占比作为该场景概率。代码很简单：

from sklearn.cluster import KMeans import numpy as np def generate_pv_scenarios(history_curves, n_clusters=20, seed=42): # history_curves: shape (n_samples, T)，每行是一条日出力曲线 kmeans = KMeans(n_clusters=n_clusters, random_state=seed, n_init=10) labels = kmeans.fit_predict(history_curves) prob = np.bincount(labels) / len(labels) centers = kmeans.cluster_centers_ return prob, centers

聚类方法的缺点是场景曲线是类中心，不是真实观测曲线，极端天气情况会被磨平。更贴近论文复现需求的是随机规划里经典的“后向场景削减法”：先采样大量场景，然后反复合并距离最近的两个场景，把小概率场景的概率叠加到大概率场景上，直到数量降到目标值。这样保留的是真实场景，概率分布也更合理。

不管用哪种方法，场景概率加起来必须等于1。这个我在调试时犯过低级错误，归一化漏了，目标值直接虚高了好几万MWh，排查了半天才发现。

4.3 最容易翻车的四个数据口径问题

第一是功率和电量的换算。模型目标函数里Δt的单位必须是小时。如果调度时段是15分钟，Δt就是0.25，不能拿0去乘；否则功率MW和电量MWh直接混用，目标值错得离谱。

第二是流量和库容的单位一致性。入库流量单位如果是m³/s，乘以Δt秒数就是m³。但库容上下限如果是万m³，中间就差了10^4倍。复现时我遇到过库容平衡约束一直报不可行，最后发现是一个数值干掉了三个数量级。

第三是梯级时滞的处理层级。时滞参数τ对应的是时段数，如果调度时段是15分钟而水流延迟是40分钟，τ就是2或者3，必须取整并单独处理。有的论文直接把时滞忽略掉，复现时你需要注意，如果忽略时滞导致的结果异常，往往不是求解器的问题，而是模型结构比论文简化太多。

第四是光伏场景和径流场景的关系。径流也有随机性，但短期调度的径流误差通常比光伏小很多；多数论文把径流作为确定性输入，只对光伏做随机场景。复现时要确认清楚你的模型里哪些是随机变量、哪些是确定性参数，别把两套随机过程混在一起。

5. 求解与调试中的实测经验：收敛性、数值尺度与常见报错

5.1 场景数到底取多少：计算代价和精度之间的平衡

场景数越多，不确定性刻画越精细，求解规模线性增长，求解时间却往往超线性增长。我自己跑小算例的经验是：2级电站、1个光伏场、24时段、20个场景的MILP模型，Gurobi十几秒能求解；场景加到50个，需要三到四分钟；如果再加机组启停的二进制变量，可能直接飙升到半小时以上。

所以在做算例对比的时候，不建议一上来就跑50个场景。先从5个场景把模型调通，确认约束没问题，再逐步加场景数，同时观察目标值变化。如果场景数从20增加到40，目标值变化不到2%，说明20个场景已经足够；如果目标值还在显著变化，再提高场景数。这个技巧能帮你把调试时间压缩一大半。

5.2 数值病态：量纲不一致是最大隐性杀手

复现时最让我头疼的不是模型不可行，而是求解器给出的结果表面正常、实际物理经济性不合逻辑。后来才发现是数值尺度问题。目标函数里，库容可能是几亿m³，出力几千MW，概率零点几，最终目标值是几万MWh，各数量级差出10^8，Gurobi算起来非常吃力，甚至会警告Numerical trouble。

解决方案很直接：把所有变量统一到同一数量级。库容单位改成万m³，流量单位用万m³/时段或者m³/s，出力用MW，目标系数自然就落在合理区间。如果求解器还是数值困难，就把目标函数整体除以10000，再不行就把Gurobi的NumericFocus参数设为3，让求解器花更多精力做数值净化。

5.3 不可行问题的定位思路

模型跑出来显示infeasible或者unbounded，第一反应不要改参数瞎试，按照下面这个思路排查：

1. 先单独固定光伏出力场景为0，去掉消纳约束，只保留水量平衡和库容约束，看模型是否可行。如果这都不可行，问题出在水库本身的初始库容、入库流量和库容上下限之间的冲突。
2. 加入光伏场景后如果不可行，检查是不是光伏预测出力在某些时段高于电网消纳上限，而模型又加了一条“必须用完光伏”的约束。正常情况下模型应该是可以弃光的，如果论文模型不允许弃光，就会在午后时段无解。
3. 加入梯级时滞后不可行，检查上游出库流量经过时滞的那一项是否超出了下游库容空余。汛期尤其是这样，上游集中放水，下游水库装不下，必须允许下游弃水才能保住可行解。
4. 还定位不了，直接用Gurobi的computeIIS()提取不可行约束集合，一条一条看是哪些约束打架，这是最权威的定位手段。

5.4 MIP求解的参数经验

这类模型如果只有连续变量，其实是线性规划，Gurobi求解很快。但一旦涉及机组启停、分段线性化或者整数决策，就变成MILP，需要关注MIP Gap。我一般把MIPGap设为0.1%，再设一个TimeLimit兜底，保证不管能不能证明最优，都能在可接受时间内拿到一个工程可用解。这两个参数非常关键：

m.Params.MIPGap = 0.001 m.Params.TimeLimit = 600 m.Params.Threads = 4

Threads设成4或者8可以让Gurobi并行求解，多核机器上提速明显。如果模型本身就很大，MemoryLimit也要关注，防止内存溢出。

6. 结果怎么看、对比怎么做：复现完只是起点

6.1 四个必须检查的结果现象

拿到优化结果后，不要急着画图贴进报告，先做四步物理合理性检查。

第一，看弃光时段。弃光应该集中在光伏大发的中午时段；如果出现清晨或傍晚弃光，模型肯定有问题，大概率是水电出力下限或者电网消纳约束写错了。

第二，看弃水时段。弃水应该出现在入库流量大、或者水库接近满库需要被迫加大泄流的时段。如果丰水期弃水量反而低于枯水期，那水量平衡可能没考虑径流随机性。

第三，看库容过程线。正常结果里，库容通常不会长期贴着上限或者下限运行。如果一张图里库容线全程压在上限，说明约束太紧，模型无非是在把水“放在水库里保底线”；如果全程压在下限，说明水电根本没发挥调节作用，光伏直接把水电的发电空间挤没了。

第四，看各场景的目标值离散程度。20个场景下，各场景的消纳电量应该有一个合理的方差：光伏大发场景消纳电量偏高，光伏小发场景偏低。如果所有场景的目标值几乎一致，说明场景生成出了问题，概率全部集中到了平均曲线上。

6.2 设置几个关键对比实验

复现论文时，一定要设计对照组，否则你根本不知道模型的价值在哪。我建议至少做三组对比。

第一组：确定性模型对照随机模型。把光伏场景概率集中到一条预测曲线上，跑确定性模型，对比随机模型的期望消纳电量。这一步能直接验证随机优化是否真的提升了期望值，也是论文里最常展示的图表。

第二组：纯水电模型对照水光互补模型。关掉光伏电站，只做梯级水电优化，再打开光伏做水光联合优化。两组结果的差值就是光伏接入带来的增量消纳空间，也是“互补价值”的量化体现。

第三组：不同消纳上限的敏感性分析。把电网消纳上限从某个基准值逐渐放松，观察可消纳电量期望和弃光率的变化曲线。这样的结果放到论文讨论部分很有说服力，也能帮你理解模型在哪个边界条件下开始“瓶颈凸显”。

整理成一个简单的对比表格式：

实际数值取决于你的数据和参数，但趋势上随机模型的结果应该不低于确定性模型，水光互补模型的总消纳电量会明显高于单独水电模型。

6.3 复现之后还能往哪些方向扩展

模型能跑通之后，如果你要做自己的研究，有几条比较自然的扩展路线。

一是加入风电和多能互补，把光伏场景生成扩展成风光联合场景，需要考虑风、光出力之间的相关性，场景生成难度会明显增加。二是把目标从单目标扩展到多目标，比如同时考虑消纳电量期望最大化、水库生态流量保障和发电经济效益，用约束法或者加权法处理。三是把确定性场景换成分布鲁棒，不精确知道场景概率，只给一个概率分布的模糊集合，优化最坏情况下的期望值，这对消纳边界附近的风险控制特别有意义。四是加入机组组合层，在出力计划之上再考虑水电机组的启停状态和最小运行时间，模型复杂度上一个台阶，但更贴近真实电厂运行约束。

我个人建议，扩展方向不要贪多，先把你复现的这个基础模型跑扎实，再挑一个方向做深做透。

最后分享一个我在整个复现过程中体会最深的点：复现EI论文的最大价值，不在于你拿到了一套能跑的代码，而在于把论文里没写出来的那些隐藏假设和数据口径，一个个用自己的代码验证了一遍。我可能花了三天时间在场景概率归一化这种看似无关紧要的细节上排错，但也正因为经历过这些，才真正理解了期望值模型和确定性模型之间的本质差异。如果你也在复现这个题目，不用急着追求一步到位，先把小算例跑通，再逐步加场景、加约束，最后你会发现，模型的规律慢慢都会变得符合物理直觉，那才是真正“复现”成功的时刻。