梯级水光互补系统短期优化调度：基于Python的期望值建模与实现

这段时间一直在啃一篇EI期刊上的梯级水光互补系统短期优化调度论文，核心思路其实一句话就能讲明白：光伏出力存在不确定性，模型要在这种随机波动环境下，合理安排梯级水电站各时段的发电流量与弃水策略，使“水电+光伏”整个系统能被电网消纳的电量期望值最大化。听起来不复杂，但真上手复现就会发现，里面既有梯级水库上下游水力联系、时滞这类硬约束，又有光伏随机场景建模、期望值离散化这些需要仔细处理的环节，最后还要用Python把整套线性规划框架搭起来，跑出可解释的调度曲线。

这篇博文就把我从问题拆解、数学模型搭建到Python代码实现、结果验证的完整过程梳理一遍，踩过的坑也会一并列出来。正在做电力系统优化调度、需要复现EI论文、或者写相关毕业设计的同学，应该能从这里找到一套可以直接借鉴的实践路径。

1. 项目背景与问题拆解

1.1 梯级水光互补为什么难

先理解“梯级水电站”这个概念。它指的是同一条河流上从上到下串联布置的多座水电站，上游电站的出库流量会直接进入下游电站的水库，成为下游的入库来水。这种水力联系意味着你不可能独立地调度某一座电站——上游多发电，下游入库就增加；上游为了调峰而大量弃水，下游水位也会跟着上涨，处理不当甚至会引发弃水连锁反应。

再加上水流在河道里还有传播时滞，上游的出库流量往往要经过若干小时才能到达下游电站，这让本来就有强耦合关系的水电调度更加棘手。我复现的论文里取了两级梯级电站、时滞取1小时，实际工程中三到五级梯级、时滞数小时的情况很常见，模型维度和约束数量会成倍增长。

光伏这边则完全是另一种风格：出力具有显著的间歇性和随机性，早高峰和晚高峰完全不可控，预测值和实际值之间经常差出10%到20%。水电的优势是响应快、调节能力强，理论上可以平抑光伏波动，但水电站自身有一堆物理约束——库容上下限、发电流量上限、末库容要求、生态流量要求——不是想怎么调就怎么调。所以“互补”两个字听起来美好，真正落地就需要一个短期优化调度模型，在日前就把水电各时段的运行计划定下来，同时考虑光伏的不确定性，让整个系统在各种可能的光照场景下都能做到尽量多消纳电量。

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

拆解这个标题，最需要想清楚的是“可消纳电量期望”这七个字。

“可消纳电量”指的是系统总出力中被电网实际接受的电量。现实中很多区域有外送通道容量限制，通道上限就是一条硬约束：水电和光伏的总出力不能超过这个值。当光伏大发时，通道被占满，水电就得压出力，甚至可能出现光伏出力超出通道上限而被迫弃光的现象。模型的目标就是通过优化水电梯级的调度，在满足所有物理约束的前提下，让系统被消纳的总电量尽可能多。

“期望”则是因为光伏出力是随机的。我们不能只看光伏预测曲线就做决策，因为预测总有偏差，如果只按预测值优化，实际运行时一旦光伏比预测低或高，调度方案就不一定最优。工程上常用的做法是用多个场景描述光伏的可能出力，每个场景给一个概率，目标函数对所有场景下的可消纳电量取加权平均，也就是期望值。这种把随机性问题转化为确定性场景优化问题的思路，本质上就是随机规划里的期望值模型。

还需要注意，目标函数里水电电量也要统计进去。因为可消纳电量是一个整体概念，水电让路给光伏导致自身少发，如果光伏消纳增加的量不足以弥补水电减少的量，那这个调度方案就是亏的。模型会自动权衡：白天光伏大发时段，水电主动压低出力甚至停发，把水量蓄在水库里；等傍晚光伏消退，再加大发电流量把水库放下来。整个过程的目标不是单纯追求水电发电量最大，也不是单纯追求消纳最多光伏，而是两者加在一起的总期望电量最大。这也是这个模型和常规水电调度模型最本质的区别。

2. 数学模型搭建：从论文公式到可求解形式

2.1 基础数据与集合定义

建模之前先把集合和参数定义清楚，这是所有优化模型的起点。我这次复现取了一个中等规模算例：2级梯级水电站、24个调度时段（1小时一个时段）、20个光伏场景。之所以没有上来就做上百个场景，是因为复现阶段首要任务是把模型逻辑跑通，场景太多只会让调试变得困难。

集合方面需要三类：水电站集合 $I$、时段集合 $T={1,...,24}$、光伏场景集合 $S={1,...,20}$。每个场景等概率，即 $\pi_s = 1/20$，这样目标函数里的期望值就是所有场景下光伏消纳量的算术平均。

水电站参数包括：库容上下限 $V_{i}^{min}, V_{i}^{max}$、初始库容和末库容、最大发电流量 $q_{i}^{max}$、综合出力系数 $K_i$、区间入流、上下游连接关系和时滞。这里有一个单位问题需要特别提醒：库容常用万m³表示，流量常用m³/s表示，但优化模型里如果直接混用，数值量级可能相差几个数量级，求解器容易数值病态。我的做法是把所有流量统一折算成万m³/h，转换关系是1 m³/s = 0.36 万m³/h，库容单位保持万m³，功率单位用MW。这样所有约束里的系数都在同一个数量级上，求解稳定得多。

光伏部分的输入是一条典型的日前预测曲线，峰值约35 MW，夜间为0。场景生成时在预测值上叠加随机误差项 $\varepsilon_t \sim N(0, 0.1)$，并做非负截断，保证场景出力不会出现负值。这个过程在后面的代码里会详细展开。

2.2 目标函数与关键约束

目标函数可以写成：

$$ \max \quad \sum_{i,t} P_{i,t}^{H} \Delta t + \sum_{s,t} \pi_s P_{s,t}^{PV,use} \Delta t $$

第一项是所有水电站各时段的发电量之和，第二项是各光伏场景下实际消纳光伏电量的期望值。其中的决策变量包括：

• $q_{i,t}^{R}$：电站 $i$ 在时段 $t$ 的发电流量（万m³/h）
• $s_{i,t}$：弃水流量（万m³/h）
• $V_{i,t}$：库容（万m³）
• $P_{s,t}^{PV,use}$：场景 $s$ 下时段 $t$ 的光伏实际消纳功率（MW）

水电出力采用线性化关系：$P_{i,t}^{H} = K_i \cdot q_{i,t}^{R}$。这里假设水头恒定，综合出力系数 $K_i$ 把流量直接映射到功率。论文里往往给出的是水头-库容曲线和出力系数，实际复现时如果不想引入非线性，对短期调度做恒定水头假设是常见且合理的简化。我算例里 $K_1=0.35$ MW/(m³/s)，折算到万m³/h单位后约0.972 MW/(万m³/h)。

核心约束分几组。

首先是水量平衡方程：

$$ V_{i,t+1} = V_{i,t} + I_{i,t} + \sum_{j\in U(i)} \left(q_{j,t-\tau_{ji}}^{R} + s_{j,t-\tau_{ji}}\right) - q_{i,t}^{R} - s_{i,t} $$

其中 $I_{i,t}$ 是区间入流，$U(i)$ 是直接上游电站集合。要注意上游出库包括发电流量和弃水两部分，这两者最终都会成为下游入库。时滞 $\tau_{ji}$ 表示水流从上游电站流到下游电站需要的时间，我取1小时，所以 t 时刻下游入库对应的是 t-1 时刻上游出库。

其次是库容约束与边界条件：

$$ V_{i}^{min} \le V_{i,t} \le V_{i}^{max}, \quad V_{i,1} = V_i^{init}, \quad V_{i,T+1} = V_i^{end} $$

初始和末库容固定，这是为了让调度结果具有可重复性，模拟一个完整调度周期的水库运行状态，避免模型为了多发电在最后时段把水库放空。

然后是外送通道约束，这是“可消纳”的直接来源：

$$ \sum_i P_{i,t}^{H} + P_{s,t}^{PV,use} \le P_{line}^{max}, \quad \forall s,t $$

意味着任一光伏场景下，水电和光伏的实时总出力都不能超过通道上限。这道约束把随机场景和确定性水电计划耦合在了一起。最后还有光伏消纳界限约束：

$$ 0 \le P_{s,t}^{PV,use} \le P_{s,t}^{PV}, \quad \forall s,t $$

即每个场景下实际消纳的光伏功率不能超过该场景的光伏出力，多余的功率就是弃光。

值得说明的是，我刻意没有给水电出力加最小技术出力约束。实际工程中水电机组有最小出力限制，但加入后会引入整数变量变成MILP问题，复现阶段先用纯LP把模型跑通，后续需要再扩展即可。如果读者遇到论文里有最小出力约束，可以把它表示成 $P_{i,t}^{H} \ge P_{i}^{min} \cdot u_{i,t}$，配一个启停变量，模型就从LP变成MILP了。

2.3 场景生成与期望值离散化

场景生成是整个随机规划模型的输入基础，做不好后面全是白搭。论文里通常不会给出现成的光伏场景数据，只会给误差分布假设，所以需要自己动手生成。

我采用的做法是：以预测曲线 $P_t^{fore}$ 为基准，对每个场景独立抽样生成误差序列 $\varepsilon_t \sim N(0, \sigma)$，场景出力为 $\max(0, P_t^{fore}(1+\varepsilon_t))$。标准差取0.1，即光伏预测误差约10%，这个量级和实际短期预测水平大致吻合。

这里有个经验：直接用蒙特卡洛随机抽样生成的20个场景，可能会在极端情况下出现个别场景光伏出力异常高或异常低，导致目标函数值对场景集合特别敏感。更稳妥的做法是用拉丁超立方抽样代替纯随机抽样，让场景在概率空间里分布更均匀；如果场景数量多到几百上千，还可以用K-means聚类做场景削减，用少量典型场景代表整体分布。我这次为了代码可读性先用了纯随机抽样加固定随机种子，方便复现结果，后面再换成更精细的场景生成方法。

场景概率在期望值模型里默认等概率 $\pi_s = 1/S$。这不是唯一选择，如果论文给的是历史日光伏数据的经验分布，也可以按各场景出现的频率赋不同概率。但等概率的好处是代码简单，而且当场景数量足够多时，等概率的算术平均就是蒙特卡洛积分的标准形式，期望值估计是无偏的。

3. Python代码实现与求解

3.1 环境准备与数据构造

代码层面我选择PuLP作为建模工具，搭配默认的CBC求解器。PuLP的语法接近数学表达式，非常适合这种论文复现场景，装起来也方便：

pip install numpy matplotlib pulp

如果你之前没装过Python环境，先去官网装个3.9以上的版本，然后用上面的命令一次装齐所有依赖。我在一台Windows机器和一台Linux服务器上都跑过这套代码，没有任何平台差异问题。

算例的具体参数如下表所示：

外送通道上限取30 MW，光伏预测曲线峰值按35 MW设计，这样白天光伏大发时段通道会被占满，水电必须压出力，夜间光伏为零，通道容量全部让给水电，模型才能真正体现出“互补”调度。

3.2 模型代码逐段拆解

先构造数据和光伏场景：

import numpy as np import pulp as pl # 基础参数 T = 24 # 时段数 S = 20 # 光伏场景数 dt = 1.0 # 时段长度(h) # 水电站参数 K_m3s = [0.35, 0.30] # 综合出力系数 MW/(m³/s)，水头近似恒定 q_max_m3s = [50.0, 45.0] # 最大发电流量 m³/s q_max = [q * 0.36 for q in q_max_m3s] # 折算成万m³/h K = [k / 0.36 for k in K_m3s] # MW/(万m³/h) V_max = [80.0, 60.0] # 万m³ V_min = [20.0, 10.0] V_init = [50.0, 30.0] V_end = [50.0, 30.0] inflow_m3s = [5.0, 3.0] # 区间入流 m³/s inflow = [x * 0.36 for x in inflow_m3s] # 万m³/h tau = 1 # 上游到下游的水流时滞(h) P_line_max = 30.0 # 外送通道上限 MW # 光伏日前预测曲线(MW) pv_forecast = np.array([ 0, 0, 0, 0, 0, 2, 5, 12, 18, 25, 30, 35, 32, 27, 22, 14, 8, 3, 0, 0, 0, 0, 0, 0 ], dtype=float) # 生成20个光伏场景：预测值 * (1 + 正态误差)，并做非负截断 np.random.seed(42) pv_scenarios = np.zeros((S, T)) for s in range(S): eps = np.random.normal(0, 0.1, T) pv_scenarios[s] = np.maximum(0, pv_forecast * (1 + eps))

这里核心的转换关系是：1 m³/s = 0.36 万m³/h。所以最大发电流量50 m³/s折算后是18万m³/h，出力系数0.35 MW/(m³/s)折算后约0.972 MW/(万m³/h)，即每放1万m³水可以发约0.972 MWh电。这些系数在后面的目标函数和约束里频繁出现，单位不统一是很多新手复现失败的第一大原因。

接下来创建优化问题和决策变量：

prob = pl.LpProblem("Hydro_Solar_Max_Expected_Absorption", pl.LpMaximize) # 决策变量 q_r = pl.LpVariable.dicts("q_r", ((i, t) for i in range(2) for t in range(T)), lowBound=0, cat=pl.LpContinuous) spill = pl.LpVariable.dicts("spill", ((i, t) for i in range(2) for t in range(T)), lowBound=0, cat=pl.LpContinuous) V = pl.LpVariable.dicts("V", ((i, t) for i in range(2) for t in range(T + 1)), lowBound=0, cat=pl.LpContinuous) pv_use = pl.LpVariable.dicts("pv_use", ((s, t) for s in range(S) for t in range(T)), lowBound=0, cat=pl.LpContinuous)

注意 (V) 变量的时段下标是 (0\sim T)，比调度时段多一个，因为要同时表示初始库容和每个时段结束后的库容。(q_r) 和 (spill) 的时段下标是 (0\sim T-1)，对应24个调度时段。这样下标不会越界。

然后写核心约束：

# 库容约束与边界 for i in range(2): prob += V[(i, 0)] == V_init[i] prob += V[(i, T)] == V_end[i] for t in range(T + 1): prob += V[(i, t)] >= V_min[i] prob += V[(i, t)] <= V_max[i] for t in range(T): prob += q_r[(i, t)] <= q_max[i] # 发电流量上限 # 水量平衡方程 for i in range(2): for t in range(T): upper_inflow = 0.0 if i > 0 and t >= tau: upper_inflow = q_r[(i - 1, t - tau)] + spill[(i - 1, t - tau)] prob += V[(i, t + 1)] == V[(i, t)] + upper_inflow + inflow[i] \ - q_r[(i, t)] - spill[(i, t)] # 外送通道约束 + 光伏消纳界限 for s in range(S): for t in range(T): prob += pl.lpSum(K[i] * q_r[(i, t)] for i in range(2)) \ + pv_use[(s, t)] <= P_line_max prob += pv_use[(s, t)] <= pv_scenarios[s, t]

水量平衡方程是整个模型的脊梁骨。注意upper_inflow只对下游电站生效，且只有 (t \ge tau) 时才加上游出库，因为 t=0 时上游出库还没流到下游。上游电站本身只有区间入流，没有来自更上游的电站。弃水流量 (spill) 和发电流量 (q_r) 一起构成了水库的总出库，这两部分都会进入下游水库。

最后是目标函数：

hydropower = pl.lpSum(K[i] * q_r[(i, t)] for i in range(2) for t in range(T)) pv_power_exp = pl.lpSum(pv_use[(s, t)] for s in range(S) for t in range(T)) / S prob += hydropower + pv_power_exp

当场景等概率时，期望值就是算术平均。目标函数的量纲是MWh，第一项水电发电量是确定性的，第二项是光伏消纳电量的期望。因为时段长度 (\Delta t=1) 小时，功率数值乘1就是电量值，所以代码里不需要额外乘时间系数。

3.3 求解与结果输出

求解只需要一行：

status = prob.solve() print("求解状态:", pl.LpStatus[status]) print("目标值（最大可消纳电量期望）:", pl.value(prob.objective), "MWh")

CBC求解这个模型非常快。整个问题变量数量大约是 (2\times24 + 2\times24 + 2\times25 + 20\times24 = 578) 个，约束数量约 (2\times24 + 20\times24 + S\times T + ... \approx 600) 条，属于小规模LP，CBC几秒钟就能解出来。如果用Gurobi或CPLEX，基本是零延迟。

求解完成后提取调度结果，准备画图：

# 提取水电出力和光伏消纳结果 p_h = {(i, t): K[i] * q_r[(i, t)].value() for i in range(2) for t in range(T)} pv_use_mean = [np.mean([pv_use[(s, t)].value() for s in range(S)]) for t in range(T)] pv_gen_mean = [np.mean([pv_scenarios[s, t] for s in range(S)]) for t in range(T)] # 每个时段系统总外送功率（期望值） p_total = [p_h[(0, t)].value() + p_h[(1, t)].value() + pv_use_mean[t] for t in range(T)] # 弃光率 curtail = np.array([pv_gen_mean[t] - pv_use_mean[t] for t in range(T)])

把这几段拼起来就是完整可运行的脚本。我自己在这个模型上跑了不下二十次，每次修改约束或参数后都会盯三张图：库容曲线有没有越界、各时段外送功率有没有超过通道上限、弃光时段是否与光伏大发时段对应。这三张图能覆盖90%的模型调试需求。

4. 结果分析与可行性验证

4.1 调度曲线解读

以一次20场景的运行结果为例，调度曲线体现出很强的规律性。夜间光伏出力为0，系统外送功率全部来自水电，第一级和第二级电站基本接近满发，库容持续下降。早上6点后光伏开始爬坡，外送通道逐渐被光伏占满，水电出力按比例压缩，库容下降速度放缓甚至开始回升。午间光伏达到峰值时，水电出力被压到很低的水平，这时上游来水继续流入水库，库容明显上涨，相当于把水暂时“存”起来，等光伏消退后再放水发电。傍晚光伏快速下降，水电出力重新爬升，库容再次回落，到24时末正好回到设定的末库容值。

这个过程的本质是用水库的蓄能来平抑光伏的间歇波动。库容曲线像一条平滑的“U型”或“V型”曲线，波动幅度完全取决于光伏和外送通道的相对关系。如果外送通道上限远大于光伏峰值，库容曲线就会平缓很多，水电不需要大幅让路；反之通道越紧张，水库削峰填谷的作用越明显，调度曲线也越“极端”。

检查结果是否有意义，可以看三个指标：所有时段的 (P_{total} \le P_{line}^{max}) 是否严格满足；库容是否始终在上下限之间、末库容是否精确回到设定值；弃光曲线是否集中在光伏大发时段。我的运行结果全部满足前两条，弃光主要出现在午间峰值时段，量级在几MWh，整体弃光率约6%到8%，符合这类互补系统的典型表现。

4.2 不同光伏场景下消纳效果对比

为了验证“期望值”模型比单纯用预测值优化更好，我做了一组对比实验：第一组用多场景期望目标，第二组把光伏预测曲线当作确定值输入，其他参数完全一样。结果差异主要出现在午间时段。

确定值模型只保证预测场景下不弃光或少弃光，但真实光伏低于预测时，系统外送功率达不到通道上限，通道容量被浪费；真实光伏高于预测时，又因为水电计划没有预留调节空间而被迫弃光。期望值模型则不同，它在20个可能场景下同步寻优，做出的水电计划是“平均最优”的，虽然单个场景下未必是全局最优，但所有场景的平均表现更好。

这个现象对应的专业术语叫“调度计划的鲁棒性”。干这行时间长了你会发现，电网调度最怕的不是某个场景下效率低，而是实际场景和计划场景偏差太大导致需要大量人工干预。用期望值模型制定的水电计划，本身就兼顾了各种可能性，运行阶段的可操作性明显更强。

需要提醒的是，期望值模型并不会消除极端场景的风险。如果论文要求在极端天气场景下也不能大量弃光或外送越限，那就需要在目标函数里加条件风险价值约束，或者把某些约束改成机会约束。这些扩展我在最后一部分会简单说明。

5. 复现过程中踩过的坑

5.1 求解器选型与性能

PuLP默认的CBC求解器对付中小规模LP非常靠谱，但如果你把场景数加到500以上，CBC的求解时间会明显上升，这时可以考虑切换到Gurobi或CPLEX。学术许可免费，安装也不复杂，用法只需要把求解器选项传给prob.solve()即可。

另外一个坑是PuLP的变量索引写法。最初我用LpVariable.dicts("q_r", range(2), range(T))这种二维写法，结果取值时总是需要用[(i,t)]还是[i][t]去猜，非常容易出错。后来统一改用生成器表达式构造变量的元组索引，代码清晰很多。这种写法在约束里遍历for i in range(2) for t in range(T)时特别顺手，强烈建议照这个模式写。

5.2 模型病态与数值问题

复现过程中遇到最多的问题不是逻辑错误，而是数值问题。最典型的是无界解和不可行解交替出现。

无界解通常是因为水库没有设置末库容约束，或者没有设置发电流量上限。模型为了让目标函数无穷大，会无限放水发电，这在物理上当然不可能。解决办法是给所有决策变量设置合理的上下界，尤其是发电流量上限和库容上限。

不可行解则往往是约束自相矛盾。我遇到过初始库容、区间入流和末库容三者不匹配，导致无论怎么调度末库容都到不了设定值的情况。调试时有一个很实用的技巧：先把末库容约束放开，跑一遍看库容终值落在哪里，再回头调整初始库容或区间入流。还有一个技巧是用prob.writeLP("model.lp")把模型导出成LP格式文件，用文本编辑器打开检查每一行约束，变量和系数是否合理一目了然。

还有一个数值层面的坑：如果水库库容用m³表示、流量用m³/s表示，数值量级可能差到六七个数量级，CBC的容差设置很容易把一些本来就该满足的约束判断成不满足。解决方式就是我在2.1里强调的单位统一——流量全部折算成万m³/h，让所有约束的系数落在0.01到100之间。

5.3 论文参数还原与结果验证技巧

EI论文的复现，最大的障碍往往是论文没有给出全部参数。有的只写了库容和装机容量，没写区间入流和时滞；有的给了水头范围但没给综合出力系数。我的做法是：优先把这些缺失参数看作可调的，先用合理估计值跑通模型，再做敏感性分析，看目标函数对哪个参数最敏感，重点标定敏感参数。

验证复现是否成功，我看三个层面。第一是定性一致性：调度曲线的形态是否和论文里的典型结果相似，比如库容曲线是否整体可控、弃光是否集中在光伏大发时段。第二是数值合理性：最终可消纳电量期望值与按水电来水量和光伏资源量估算的物理上限是否在同一量级，如果差出好几倍，说明模型或数据大概率有问题。第三是边界条件测试：把外送通道上限改到极大，模型应该退化成一个纯水电调度问题，光伏不被弃电；把光伏场景全部设成0，模型应该等价于没有光伏的常规梯级水电调度。这两个退化测试能快速暴露目标函数或约束里隐藏的错误。

下面这个表是我实际复现过程中最常碰到的几个问题和对应的解决手段：

最后分享一个我个人的复现习惯：不要在拿到论文后立刻写代码。先花半天时间把论文的系统拓扑图画出来——哪座电站在上游、哪座在下游、时滞几小时、通道限制在哪里、光伏接入在哪个节点——然后再把目标函数和约束按“确定性部分”和“随机性部分”分开列出来。这个准备工作看起来消耗时间，实际上能帮你节省后面两三天的调试时间。通过这次复现我也意识到，EI论文的公式只是骨架，真正有价值的是参数选取的思路、模型简化的边界和数值实现的细节，这些东西恰恰是论文正文里不会明说的。我建议你复现的时候，每个简化假设都单独记录一下，后面写自己的论文时，这些都是答辩时能扛住追问的素材。