看到这个标题，我第一反应是：“又是水电、又是光伏、又是期望、又是EI复现”，这套组合拳劝退了不少人。但把这个题拆开之后你会发现，它本质上是把一个很现实的工程问题——梯级水光互补系统短期优化调度——用“最大化可消纳电量期望”这个目标串起来，再用Python的混合整数线性规划（MILP）框架落地。我这次把自己复现同类EI论文的过程完整梳理了一遍，从数学模型、代码骨架到算例解读和踩坑记录都会讲到，希望能帮你少走一点弯路。适合正在做“双碳”方向研究的学生，也适合刚接触电力系统优化调度、想用Python把论文公式跑成可执行程序的工程师。

1. 项目背景与需求拆解

1.1 为什么梯级水电和光伏总被放在一起调度

先说“梯级水光互补”到底是个什么物理场景。它不像单一水电站那样简单：一条河流上往往有好几座水库，上游电站发完电，水会流到下游水库，下游电站再发一次。所以梯级电站之间在水量和时间上都存在强耦合关系。在此基础上再加入光伏电站，就形成了“水光互补”的格局。

互补的核心逻辑很好理解：光伏出力白天波动大，尤其是阴天、云层过境时可能短时间内从80%出力掉到20%；水电则相反，只要水库里有水、机组能调节，就可以快速增减出力。光伏骤降时，水电顶上；光伏大发时，水电压低出力甚至停机，把水蓄起来留到晚上再放。听起来很理想，但是电站运行约束不会让你这么轻松，水库不能随便蓄到漫坝，也不能把水位放到死水位以下，更不能让下游断流。上下游电站之间还存在流达时间，上游的出库流量要过几个小时甚至更久才能到下游水库。

所以这个问题的复杂度主要来自三层：时间上的水库蓄放联动、空间上的上下游梯级约束、以及光伏预测误差带来的不确定性。标题里提到的“短期优化调度”，一般就是提前24到96个时段做计划，时间颗粒度常用1小时或15分钟。

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

很多人一开始会误解，以为“可消纳电量”就是“电站总共发了多少电”。实际上不是。光伏和水电发出来，不代表能被电网接受或者被本地负荷用完。外送通道有容量上限，水电和光伏同时大发时，就只能弃掉一部分。弃掉的光伏叫“弃光”，从水库溢流出去不发电的水叫“弃水”。可消纳电量，就是最终真正被电网接纳、被负荷使用的那部分电量。

“期望”这个词更关键。光伏出力在调度计划制定时还不知道准确值，只能靠预测，预测又有误差。如果只给一个预测值，然后当成确定值去优化，那叫确定性调度。而“期望”意味着要让目标在多种可能的光伏场景下都表现得好，相当于求所有场景下可消纳电量的加权平均。这是随机规划里非常典型的处理方式，也是这类EI论文比传统确定性优化更有价值的地方。

在实际复现时，我会把目标函数写成“期望可消纳电量最大化”，具体包含水电出力、被消纳的光伏出力、以及弃光和弃水惩罚项。惩罚项非常重要，如果没有惩罚，模型可能为了追求电量数字好看，让弃光弃水肆无忌惮，结果算出来方案完全不可运行。

1.3 这个复现项目适合谁

我给这个项目定过位：它不是让你从零做一个水电站监控系统，而是把论文里的数学模型“翻译”成Python代码。适合三类人：

• 电力系统专业研究生，需要复现EI/SCI论文做对比实验或者毕设基础；
• 新能源并网、水电调度方向的算法工程师，想快速验证一种调度策略；
• 想系统学习Pyomo或者MILP建模的开发者，这个例子比课本上的运输问题更有工程代入感。

如果你之前只写过简单的线性规划，这个项目正好可以用来进阶：场景生成、带概率权重的目标函数、时间耦合约束、求解器调优，这些在一般教程里很少被放在同一个例子里讲清楚。

2. 数学模型拆解

2.1 目标函数的数学形式

目标函数我习惯先写清楚再动手写代码。假设有S个光伏出力场景，每个场景的概率为 (\pi_s)，收益率期望模型可以写成：

[ \max \sum_{s \in S} \pi_s \sum_{t \in T} \left( P_{h}(s,t) + P_{pv_use}(s,t) - \lambda_c C_{pv}(s,t) - \lambda_w C_{water}(s,t) \right) \Delta t ]

其中：

• (P_h(s,t))：场景s下t时段全部梯级水电站总出力；
• (P_{pv_use}(s,t))：场景s下t时段被真正消纳的光伏出力；
• (C_{pv}(s,t))：场景s下t时段的弃光量，等于光伏最大可发功率减去实际消纳功率；
• (C_{water}(s,t))：场景s下t时段的弃水量折算成对应可发电功率；
• (\lambda_c, \lambda_w)：弃光、弃水的惩罚系数，取值要远大于正常发电收益，这样模型才会主动避免弃电；
• (\Delta t)：时段时长，小时级模型就是1。

很多论文里还会把外送通道约束松弛成惩罚，或加上阻塞电价，但初版可以先把外送通道作为硬约束，后续再根据复现需要调整。

2.2 五类核心约束条件

我复现时遇到过无数次因为少了一条约束导致结果离谱的情况。短期优化调度里，下面五类约束基本是必备的：

第一，水量平衡约束。上游水库的出库流量（发电流量加弃水流量）加上天然入库，决定了水库下一时刻的蓄水量。梯级之间还需要考虑流达延迟：

[ V_{i,t+1} = V_{i,t} + Inflow_{i,t} + Q_{in,i,t-\tau} - Q_{h,i,t} - C_{water,i,t} ]

这里的 (Q_{in,i,t-\tau}) 就是上游电站的出库流量经过 (\tau) 个时段后到达下游水库的部分。 (\tau = 0) 可以简化成不考虑流达时间。

第二，库容边界约束。每座水库的蓄水量不能超过最大库容，也不能低于死库容对应的最低水位。短期调度还必须加“日末水位恢复”约束，比如 (V_{i,T} = V_{i,0})，不然模型会把水全部放空去发电，第二天没法正常运行。

第三，水电出力约束。水电站出力与发电流量和水头有关，严格说是非线性关系。复现时最常见的是“恒定水头线性化”，即把水头视为常数，用 (P_{h,i,t} = K_i Q_{h,i,t}) 近似，其中 (K_i) 是电站的综合出力系数。如果想更精确，可以做分段线性化，把水头-库容-出力曲面切成若干线性段，配合二进制变量构建MILP。

第四，光伏消纳约束。(0 \le P_{pv_use}(s,t) \le P_{pv_forecast}(s,t))，同时弃光量等于最大可发与实际消纳之差。这样模型不会出现“消纳了不存在的电”这种蠢事。

第五，系统通道约束。同一并网点下，所有电源总出力不能超过外送通道上限：

[ \sum_i P_{h,i,t} + P_{pv_use}(s,t) \le L_t ]

有些调度模型还会加水电爬坡约束或最小出力约束，初版可以先不加，但后期对比论文时必须补上。

2.3 光伏场景与期望值的工程简化

“期望”在代码里最简单粗暴的实现方式是样本平均逼近（Sample Average Approximation，SAA）。假设光伏预测曲线是 (P_{pv}^{forecast}(t))，我用历史预测误差分布去构造场景：

[ P_{pv}(s,t) = \max\left(0, P_{pv}^{forecast}(t) \cdot (1 + \epsilon_{s,t})\right) ]

其中 (\epsilon_{s,t}) 可以服从正态分布、Beta分布，或者直接用历史实测误差样本。场景数少的时候，直接用随机数据生成几十个场景就可以；如果场景数多到求解器受不了，再做场景削减，比如同步回代削减法或者k-means聚类。复现阶段，我建议先跑10到30个场景，验证模型能求解，再逐步加到100个场景看目标值变化趋势。

3. Python代码实现与关键环节

3.1 求解器选型：为什么不用遗传算法

做这类调度优化，最常见的错误是一上来就套遗传算法或粒子群。不是说启发式算法不行，而是可消纳电量期望最大化本质上是线性目标加线性约束，完全可以用MILP精确求解，结果有全局最优性保证。EI论文里大量用的也是商业求解器加线性化模型。

我这次用的是Pyomo建模，配合CBC开源求解器先跑通流程，再用Gurobi做正式实验。Pyomo的好处是模型定义与求解器解耦，你写一套代码，换求解器只需要改一行。环境准备大致如下：

pip install pyomo numpy pandas matplotlib

CBC可以用conda安装：

conda install -c conda-forge coincbc

如果你有Gurobi学术授权，SolverFactory('gurobi')直接可用，求解速度比CBC快一个量级，尤其场景数超过50时差距极其明显。

3.2 数据结构与参数化约定

写代码之前，我强烈建议把参数和数据分开。我用四个CSV文件组织输入数据：

• reservoir.csv：每座水库的最小库容、最大库容、初始库容、目标末库容、综合出力系数；
• pv_forecast.csv：每个时段的光伏预测最大出力；
• inflow.csv：每座水库各时段的天然入库流量；
• system.csv：各时段外送通道上限、系统负荷需求。

读入之后统一用字典存起来，方便后面在Pyomo里做索引。比如：

reservoir = pd.read_csv('reservoir.csv', index_col=0) inflow = pd.read_csv('inflow.csv', index_col=0) pv_base = pd.read_csv('pv_forecast.csv')['power'].values

参数化的好处是：换一个典型日、换一组电站参数，只需要改CSV，不需要动模型代码。这对后期做多组算例对比非常省心。

3.3 核心模型构建代码框架

下面给一个简化但可直接跑通思路的Pyomo代码骨架。假设只有24个时段、两座梯级水电站、S个光伏场景。

import numpy as np from pyomo.environ import * N_H = 2 # 水电站数量 T = 24 # 调度时段数 N_S = 20 # 光伏场景数 dt = 1 # 时段时长，单位小时 # 场景概率 pi_s = np.ones(N_S) / N_S # 生成光伏场景：预测值叠加随机误差 np.random.seed(42) pv_max = pv_base * (1 + np.random.normal(0, 0.12, size=(N_S, T))) pv_max = np.clip(pv_max, 0, None) # 水库参数（示意：两座电站） Vmin = {1: 100.0, 2: 80.0} Vmax = {1: 500.0, 2: 400.0} V0 = {1: 300.0, 2: 250.0} Vtarget = {1: 300.0, 2: 250.0} K_hydro = {1: 0.8, 2: 0.9} # 综合出力系数，简化处理 model = ConcreteModel() model.H = Set(initialize=[1, 2]) model.T = RangeSet(1, T) model.S = RangeSet(1, N_S) # 变量 model.V = Var(model.H, model.T, model.S, within=NonNegativeReals) # 库容 model.Qh = Var(model.H, model.T, model.S, within=NonNegativeReals) # 发电流量 model.Sp = Var(model.H, model.T, model.S, within=NonNegativeReals) # 弃水流量 model.Ph = Var(model.H, model.T, model.S, within=NonNegativeReals) # 水电出力 model.Ppv_use = Var(model.T, model.S, within=NonNegativeReals) # 光伏实际消纳 model.Cpv = Var(model.T, model.S, within=NonNegativeReals) # 弃光 # 目标函数：期望可消纳电量最大化 def obj_rule(m): expr = 0 for s in m.S: inner = sum(m.Ph[h, t, s] for h in m.H for t in m.T) * dt inner += sum(m.Ppv_use[t, s] for t in m.T) * dt inner -= 10.0 * sum(m.Cpv[t, s] for t in m.T) * dt inner -= 20.0 * sum(m.Sp[h, t, s] for h in m.H for t in m.T) * dt expr += pi_s[s-1] * inner return expr model.obj = Objective(rule=obj_rule, sense=maximize) # 水量平衡约束 def water_balance_rule(m, h, t, s): if t == 24: return Constraint.Skip inflow = 15.0 # 实际应从数据文件读取 # 上游来水简化处理：上游电站出库流量当天就到下游 upstream_in = 0 if h == 2: upstream_in = m.Qh[1, t, s] + m.Sp[1, t, s] return m.V[h, t+1, s] == m.V[h, t, s] + inflow + upstream_in \ - m.Qh[h, t, s] - m.Sp[h, t, s] model.water_balance = Constraint(model.H, model.T, model.S, rule=water_balance_rule) # 库容上下限 def vlimit_low_rule(m, h, t, s): return m.V[h, t, s] >= Vmin[h] def vlimit_up_rule(m, h, t, s): return m.V[h, t, s] <= Vmax[h] model.v_low = Constraint(model.H, model.T, model.S, rule=vlimit_low_rule) model.v_up = Constraint(model.H, model.T, model.S, rule=vlimit_up_rule) # 初始库容与末库容 def vinit_rule(m, h, s): return m.V[h, 1, s] == V0[h] def vend_rule(m, h, s): return m.V[h, 24, s] == Vtarget[h] model.v_init = Constraint(model.H, model.S, rule=vinit_rule) model.v_end = Constraint(model.H, model.S, rule=v_end_rule) # 水电出力与发电流量的线性关系 def power_hydro_rule(m, h, t, s): return m.Ph[h, t, s] == K_hydro[h] * m.Qh[h, t, s] model.power_hydro = Constraint(model.H, model.T, model.S, rule=power_hydro_rule) # 弃光定义：光伏可发 - 光伏实消 = 弃光 def curtail_rule(m, t, s): return m.Cpv[t, s] == pv_max[s-1, t-1] - m.Ppv_use[t, s] model.curtail_def = Constraint(model.T, model.S, rule=curtail_rule) # 外送通道约束 def channel_rule(m, t, s): return sum(m.Ph[h, t, s] for h in m.H) + m.Ppv_use[t, s] <= 1000.0 model.channel = Constraint(model.T, model.S, rule=channel_rule)

需要特别提醒的是：上面这种把所有决策变量都按场景展开的写法，严格来说是“场景独立决策”，每个方案只在对应场景下最优，没有考虑决策的不可预期性。复现初版可以这么跑，但写论文对比时，通常还需要加“非预期约束”，比如某些时段的库容决策不随场景变化，或者采用两阶段随机规划的写法。我建议先按全场景展开跑通，再加非预期约束，这样更容易定位问题。

3.4 求解与结果导出

模型构建完成后，求解代码非常简单：

solver = SolverFactory('cbc') results = solver.solve(model, tee=False) if results.solver.termination_condition == TerminationCondition.optimal: print('优化成功') else: print('求解异常:', results.solver.termination_condition)

提取结果时，我最常用的是value(model.V[h, t, s])，但要小心Pyomo变量对象在循环里的使用方式，建议用字典先缓存：

V_res = {(h, t, s): value(model.V[h, t, s]) for h in model.H for t in model.T for s in model.S} Ph_res = {(h, t, s): value(model.Ph[h, t, s]) for h in model.H for t in model.T for s in model.S} Ppv_res = {(t, s): value(model.Ppv_use[t, s]) for t in model.T for s in model.S}

得到结果后，用pandas转成DataFrame导出CSV，再交给matplotlib绘制功率曲线和库容曲线。这里有一个我经常踩的坑：value()拿到的是浮点数，但有些求解器返回的不是严格数值类型，直接参与pandas计算可能出问题，最好先用float()包一层。

4. 算例实验与结果分析

4.1 典型日数据与参数设置

为了让结果可对照，我设置了一个简化的典型日算例：两座梯级水电站，光伏电站容量800MW，外送通道上限1000MW，调度周期24小时。预测光伏曲线呈典型“单峰”形态，中午时段最大出力约720MW。

水库参数参考下表：

天然入库流量取常数，电站1每小时15万m³，电站2每小时10万m³。惩罚系数设为弃光10元/MWh、弃水20元/MWh，这个惩罚要高于正常情况下单位发电收益，但又不能高到让模型为了“惩罚”而故意压低正常出力。实际论文里会用上网电价作为基准，这里先随意给一个示意值。

4.2 功率曲线与水位过程怎么读

算例跑出来后，最直观的是看各时段总出力曲线。我通常会画一张包含三条曲线的图：光伏可发功率、光伏实际消纳功率、水电总出力。

光伏大发的中午时段，外送通道上限1000MW会被逼近，水电会自动压低出力，甚至可能出现停机蓄水。到了傍晚光伏快速下降的时段，水电出力会明显抬升，补足电力缺口。这就是“互补”在曲线上的直接体现。

库容曲线则需要分开看两座电站。上游电站因为来水条件更好，通常承担主要的调蓄功能，库容变化幅度会更大；下游电站除了自己发电，还需要消化上游的出库流量，库容变化往往受上游约束影响。如果发现某座电站的水位曲线一直贴着上限或下限运行，就要注意是不是约束设置太紧，或者惩罚系数不合理。

4.3 期望消纳电量与弃光率对比

场景数对目标值的影响非常值得观察。我做了一组对比实验，同样是24时段两库系统，CBC求解器，结果大致如下（示意值）：

可以看出，场景数从1增加到50时，期望可消纳电量明显上升、弃光率明显下降。这是因为确定性模型只针对一条预测曲线优化，遇到极端偏差场景时毫无准备；而随机优化方案虽然在某些“平均场景”下不是最优，但对各种可能出现的光伏场景都预留了水电调节空间，整体期望值自然更高。

但场景数继续从50增加到100时，增益明显变缓。这是随机规划的典型现象：样本均值逼近的误差随着样本数增加而收敛，但增长速度递减。实际工程中不一定非要跑到200个场景，要根据求解时间和精度需求权衡。

5. 复现路上的坑与排查经验

5.1 模型无可行解：按什么顺序排查

这个模型最常见的失败方式是求解器直接报“infeasible”，第一次遇到时很容易让人懵。我的排查顺序固定三步：

第一步，检查始末库容和天然来水是否自洽。如果初始库容等于目标末库容，而全天然来水又不够发电流量的最低需求，模型必崩。最简单的验证方法：在目标函数里暂时去掉惩罚项，把弃水变量和弃光变量固定为0，只让水量平衡约束参与，看是否存在解。

第二步，检查水量平衡约束中时间下标是否写错。比如V[h, t+1]在t=24时越界，或者上游来水错用了t+1时段的值。这类错误在Pyomo里有时候不会直接报错，而是产生一个变量索引错误。

第三步，检查外送通道约束是否过紧。通道上限设置过低时，水电被迫压出力，但水量平衡又要求发电流量至少满足某一下限，两个约束互相打架，就完全不可行。我一般先把通道上限放松50%，确认模型能解，再逐步收紧看边界在哪里。

5.2 求解时间爆炸：场景削减和线性化并施

随机场景数增加时，MILP规模增长很快。100个场景、24时段、2个电站就已经是上万变量，CBC可能跑几分钟甚至更久。这时候我常用的处理手段有三个：

一是场景削减，不要直接拿200个原始场景糊上去。先用同步回代法把相似场景合并，比如Gurobi或者scikit-learn的KMeans都能做，把200个削到20个，目标值损失往往不超过2%。二是检查有没有非必要的整数变量。如果水头用恒定线性近似，整个模型完全没有整数变量，只是一个大规模的线性规划，求解速度会快非常多；一旦分出控制水头的二进制变量，模型立刻变成MILP，难度上一个台阶。三是给求解器设定合理的MIP Gap。工程场景没必要非要0.00%最优，gap=0.5%时结果已经足够用于对比分析。

Gurobi里设置：

solver.options['MIPGap'] = 0.005

CBC里对应的是：

solver.options['ratio'] = 0.005

5.3 数值与论文对不齐：复现到底复现什么

这是EI复现最折磨人的环节。好不容易把代码跑通，结果和论文表里的数值对不上，差个百分之五以内还好说，差得多了就怀疑人生。

我的建议是：不要只盯随机种子。就算你把随机种子完全复现成论文一样的，对比时也可能因为以下因素对不上：

• 论文的预测误差场景分布参数没有公开，你用正态分布，他可能用的Beta分布；
• 弃光弃水惩罚系数不同，哪怕目标函数公式一样，数值结果也会差很多；
• 水头处理方式不同，恒定水头和分段线性水头算出来的发电流量边界完全两码事；
• 外送通道约束是硬约束还是罚函数，对弃光率影响非常大。

所以在复现阶段，更实际的目标是复现“方法结构和规律”，比如随机优化比确定性调度能提高多少消纳期望、场景数增加后目标值如何收敛、水库水位变化趋势是否符合物理直觉。这些规律和论文一致，就说明你的模型逻辑是对的。

5.4 惩罚系数怎么设，我的一点经验

惩罚系数是这类模型里最“玄学”的地方。设太小，模型不在乎弃光，外送通道一堵就让光伏弃着；设太大，模型会为了省惩罚而让水电机组满发，反过来造成大量弃水，同样不物理。

我自己常用的基准是上网电价。假设上网电价为0.4元/kWh，那么弃光惩罚可以取1.2到2倍上网电价，弃水惩罚取2倍左右。这样模型的逻辑是：如果光伏弃掉不如水电压低出力划算，水电机组自然会压低；如果水电为了避免弃水而多发电导致更低效，惩罚会及时拉住它。调参技巧是固定其他条件，只变一个惩罚系数，看弃光率和弃水率的变化曲线，直到两者处于一个比较均衡的水平。不要指望一次调好，这个环节通常要跑很多轮。

6. 复现完之后的几点个人体会

做完整套复现，我最深的感触是：随机优化最难的其实不是求解，而是把“期望”从一句黑话变成可计算的表达。当你在代码里第一次看到sum(pi_s * 目标量)跑完并输出一个收敛的期望值时，数学模型和工程实践之间的距离一下子就缩短了。

如果后续想继续扩展，我建议可以在现有模型上做三个方向的小改：一是把恒定水头换成分段线性水头，让水电出力更贴近实际；二是加入滚动时域调度，每4小时用最新预测重新求解一次；三是把目标函数从期望最大化改成“期望+条件风险价值（CVaR）”的加权形式，用来控制极端场景下的风险。每个方向的改动都不算大，但做完之后你对调度模型的理解会上一个台阶。这个模型虽然只是简化复现，但它已经足够让你掌握“梯级水电+随机新能源”这类问题的完整建模方法论。