前阵子复现了一篇关于梯级水光互补系统短期优化调度的EI论文前后折腾了一周多才把代码完整跑通。说实话模型本身并不复杂真正费时间的反而是那些藏在公式背后的细节目标函数里的期望算子到底对谁取期望、水量平衡方程的时滞下标注没标注清楚、水电站出力曲线怎么线性化这些地方只要有一处理解偏了模型不是不可行就是结果违背常识。这篇文章是复现笔记也是踩坑记录给打算在新能源电力系统方向做EI论文复现、或者想动手实现水光互补调度模型的朋友。我会从标题拆解开始把数学模型、Python代码结构、算例结果分析和调试经验全部过一遍。先说明一点EI论文里很多参数和算例数据并不会完整公开我在复现时采用了自构造的测试算例参数量级参考公开流域数据和通用水电参数重点还原的是建模方法和求解逻辑。你手头如果拿到具体论文和原始数据把这套代码对应位置的参数替换掉就能复用。整体框架比较通用不挑特定的水库几何或光伏场站配置。1. 先把这个标题掰开揉碎梯级、水光互补、可消纳期望分别意味着什么很多同学拿到标题第一反应是搜代码这是最大的误区。标题里每个词都是信息量理解到位了后面建模就是顺水推舟的事。这一章我把“梯级水光互补系统最大化可消纳电量期望短期优化调度模型”拆成四个关键短语逐个讲清楚。1.1 梯级水电站为什么是“梯级”而不是单库梯级这个词描述的是多个水电站沿同一条河流自上而下串联布置的物理格局。上游水库放水发电之后水流不会消失而是经过一段时间流到下游水库成为下游水库的入库流量下游水库可以再次利用这部分水发电。这就是梯级电站“一水多用”的天然优势。梯级联合调度的核心耦合关系有两个一个是水量耦合上游的出库流量包含发电流量和弃水流量经延时后进入下游水库的水量平衡方程另一个是电量的时空耦合上游某个时段多放水发电可能造成下游未来几个时段入库流量增加进而影响下游的发电能力。所以梯级调度不能拆成单个水库独立优化否则下游水库的入库过程会失真优化结果在物理上不成立。短期优化调度里梯级的作用被进一步放大。光伏出力在一天内具有很强的波动性中午出力大、早晚出力小而水电可以通过调节水库放水来快速响应这种波动。上游水库可以在光伏大发时段适当蓄水、减少发电流量把水存到光伏出力下降的晚间再集中放水发电。这个“时间平移”能力就是梯级水光互补系统最核心的调节价值。1.2 “最大化可消纳电量期望”不是“最大发电量”这里最容易理解偏。题目里说的是“可消纳电量”不是“发电量”。可消纳意味着系统存在一个消纳上限这个上限可能来自外送输电通道的容量限制、电网负荷水平、电力系统安全稳定约束等。光伏发电若超出消纳空间超出部分只能丢弃也就是弃光。为什么要用“期望值”而不是直接用预测值因为光伏出力是随机的受云层、温度、日照强度影响预测值只是一个可能出现的确定值完全按预测值调度实际运行时大概率出现偏差。论文里普遍的做法是构造多个光伏出力场景每个场景对应一种可能出现的光照情况然后对每个场景下的可消纳电量取概率加权平均这个加权平均就是期望值。这里有一个建模层面的重要选择水电调度决策发生在光伏实际出力已知之前属于“现在决定”的变量而光伏消纳量是场景依赖的属于“看到场景后再调整”的变量。模型因此呈现两阶段随机优化的结构理解这一点目标函数和约束的代码写法就清晰了。1.3 短期调度的时间尺度与决策内容短期调度一般指未来24小时到72小时内的运行计划时段长度常用1小时。我复现时采用24个时段、每时段1小时的设定这是比较经典的做法。如果涉及日内滚动也可以压到15分钟一个时段但模型规模会成倍增加求解时间明显变长。调度决策的内容包括三块第一每个水库每个时段的发电流量和弃水流量这决定了水电出力过程第二每个时段光伏实际消纳功率和弃光功率第三由水量平衡推导出的水库库容变化过程。三块决策通过系统功率平衡约束耦合在一起形成完整的优化问题。一句话总结这个模型要干什么在给定光伏出力随机场景、来水预报和系统消纳上限的情况下通过优化梯级水库的放水计划和光伏消纳方案让整个系统在一天内期望可消纳电量达到最大同时保证水库运行约束安全。2. 数学模型拆解目标函数、约束条件与线性化处理这一章是复现的核心。论文里的模型通常写得比较紧凑但展开后其实就是一个目标函数加五大类约束。我会把每一部分的物理含义讲清楚并且给出可编程的数学形式。2.1 目标函数怎么表示“期望值最大化”目标函数的一般形式是最大化所有场景下系统可消纳电量的概率加权平均值论文里还会加上弃光惩罚项。可以用下面的公式表达[ \max \quad \sum_{s1}^{S} \pi_s \sum_{t1}^{T} \left[ \sum_{i1}^{N} P_{i,t}^{H} \Delta t P_{t,s}^{use} \Delta t - \lambda P_{t,s}^{shed} \Delta t \right] ]各符号的含义( S ) 是光伏出力场景数( \pi_s ) 是第 ( s ) 个场景的概率所有场景概率之和为1( P_{i,t}^{H} ) 是第 ( i ) 座水电站在时段 ( t ) 的出力( P_{t,s}^{use} ) 是场景 ( s ) 下时段 ( t ) 实际消纳的光伏功率( P_{t,s}^{shed} ) 是场景 ( s ) 下时段 ( t ) 的弃光功率( \lambda ) 是弃光惩罚系数取值通常远大于电价目的是让模型优先避免弃光。水电出力 ( P_{i,t}^{H} ) 不携带场景下标 ( s )因为水电调度计划需要在光伏随机性实现之前就确定下来这对应两阶段随机优化里的第一阶段决策。光伏消纳和弃光变量携带场景下标属于第二阶段决策。这个下标细节对编程影响很大写变量时一定要带对。弃光惩罚的本质是一个松弛机制。如果完全没有惩罚项模型在极端场景下可能通过大量弃光来换取水电发电量的增加这未必符合实际运行偏好。加上惩罚项后模型会自动在“多消纳光伏”和“多让水电发电”之间做权衡实际输出结果更贴近调度员的决策习惯。2.2 五大类约束的工程化表述我把复现中用到的约束整理成表格每一类都对应明确的物理意义。表格里的形式是编程时采用的最终形式比论文里更贴近代码实现。约束类型数学表达物理含义水量平衡( V_{i,t1} V_{i,t} (I_{i,t} Q_{i-1,t-\tau_i} S_{i-1,t-\tau_i} - Q_{i,t} - S_{i,t})\Delta t )上游来水、上游放水延时入库、本库出流共同决定库容变化库容边界( V_{i}^{min} \le V_{i,t} \le V_{i}^{max} )对应死水位和正常高水位短期调度还要加末库容约束出流边界( 0 \le Q_{i,t} \le Q_{i}^{max} ) ( 0 \le S_{i,t} \le S_{i}^{max} )机组过流能力和溢洪道泄流能力水电出力( P_{i,t}^{H} 9.81 \eta_i Q_{i,t}^{gen} H_{i,t} )水头、流量和机组效率决定发电功率光伏消纳( P_{t,s}^{use} P_{t,s}^{shed} P_{t,s}^{pv} )预测光伏分解为消纳和弃光两部分系统功率平衡( \sum_{i} P_{i,t}^{H} P_{t,s}^{use} \le P_{t}^{grid} )水电与光伏实际消纳之和不超过外送或负荷上限梯级水量平衡中那个 ( \tau_i ) 需要单独说明它是第 ( i-1 ) 座水库到第 ( i ) 座水库的水流传播时间。短期调度中如果流域范围不大( \tau_i ) 通常取0到2小时。我当时设置上游到下游的传播时间为1小时也就是上游 ( t ) 时段的出库最早在 ( t1 ) 时段才能进入下游水库的平衡方程。末库容约束是短期调度里特别容易被忽略的一类。如果只约束库容上下限模型大概率会把水库放到最高水位或者排到最低水位因为从单日收益角度看这往往是最优的但现实中这会导致第二天无法正常运行。常规做法是加 ( V_{i,T}^{min, end} \le V_{i,T} \le V_{i,T}^{max, end} )让调度周期结束时的库容回到一个合理区间。2.3 从论文公式到可求解模型线性化三板斧水电出力特性是模型里唯一的非线性来源。( P 9.81 \eta Q H ) 中发电流量 ( Q ) 是决策变量发电水头 ( H ) 又随库容和尾水位变化两个变量相乘就形成了双线性项。直接交给Gurobi这类求解器会报非凸非线性错误必须先线性化。我在复现中试过三种做法分别适用于不同精度要求第一种是固定水头法。对于调节能力不强、水位波动小的径流式电站可以把水头近似看作常数出力变成流量的线性函数。这是最简单的方式求解速度快但碰上日调节能力强的水库误差偏大。第二种是分段线性化。把发电流量分成若干段每段对应一个线性出力区间用SOS2约束保证模型在相邻分段之间插值。分段越多精度越高但变量数量随之增加。实际操作中我建议水电出力曲线分5到8段即可再多对结果影响很小求解时间却明显上升。第三种是McCormick包络法。对 ( z Q H ) 这个双线性项用四个线性不等式约束把它包围在凸包内。对于短期调度这种单峰曲线问题McCormick松弛的质量通常还不错。写出通用代码的示例是这样的# z Q * H 的 McCormick 线性化 # Q_lb, Q_ub 为流量边界H_lb, H_ub 为水头边界 m.addConstr(z Q_lb * H Q * H_lb - Q_lb * H_lb) m.addConstr(z Q_ub * H Q * H_ub - Q_ub * H_ub) m.addConstr(z Q_lb * H Q * H_ub - Q_lb * H_ub) m.addConstr(z Q_ub * H Q * H_lb - Q_ub * H_lb)这套线性化处理之后模型就变成混合整数线性规划或者线性规划Gurobi可以直接求解。如果模型里没有0-1整数变量线性规划求解速度非常快这也是短期调度能跑出高时效性的基础。3. Python代码架构从数据准备到Gurobi求解模型理解到位后代码实现就是水到渠成的事。我的代码结构分为四层数据定义层、场景生成层、模型构建层、结果输出层。这样分的好处是以后换算例只要改参数文件模型构建部分完全不用动。3.1 数据准备与单位统一这是复现中最容易被坑死的一步。论文里的库容单位经常写“亿立方米”来水单位却是“立方米每秒”两个单位直接相乘水量平衡方程会差出好几个数量级模型立刻不可行。我统一把所有流量换算成立方米每小时或者干脆在方程中乘以时间步长 ( \Delta t ) 来处理单位。具体做法是流量 ( m^3/s ) 乘以3600变成 ( m^3/h )再乘以时段长度1小时就得到该时段的体积。如果用半天或一小时之外的时段步长需要相应调整3600这个系数。我测试算例采用的参数如下供参考参数数值梯级水库数2水电站在线台数3光伏场站装机600 MW系统外送通道上限900 MW调度周期24 h水流传播时间1 h光伏场景数5削减后弃光惩罚系数200 元/MWh数据在代码里用字典和numpy数组组织。光伏预测出力是一个 ( 24 \times S ) 的二维数组每一列对应一条场景曲线来水预报是一维长度24的数组水库的参数库容上下限、最大出流、初始库容、末库容单独存成字典。这样主程序读起来非常清晰。3.2 光伏随机场景生成与削减光伏预测曲线一般只有一条但模型需要的是多个场景。我采用的做法是蒙特卡洛加场景削减。先生成200条场景在预测值基础上叠加均值为零、方差随时间变化的正态扰动中午光照强时方差大早晚方差小。然后在场景集上做削减。场景削减最实用的方案有两种同步回代消除法和KMeans聚类法。同步回代法更严谨它通过迭代删除与其他场景距离最近的场景并把概率转移给邻近场景最终保留的场景集在概率分布意义上最优。KMeans实现简单速度快在场景数不是特别大的情况下效果足够。我先后对比了两者最终用的是KMeans加概率归一的组合。from sklearn.cluster import KMeans import numpy as np # pv_samples: shape (N_scen, T)每行是一条光伏场景 k 5 kmeans KMeans(n_clustersk, random_state42).fit(pv_samples) centers kmeans.cluster_centers_ # 削减后的典型场景 labels kmeans.labels_ probs np.bincount(labels, minlengthk) / len(labels) # 各场景概率削减后得到5条典型场景曲线和对应的概率。这5条曲线进入目标函数时每个场景的权重就是 ( \pi_s )。需要特别提醒的是削减后的场景数量不要太多因为每个场景都会复制一份光伏消纳变量和弃光变量场景数从5增加到20模型规模大好几倍求解时间成倍增长但目标值改善往往很有限。3.3 模型构建核心代码模型构建我用的是Gurobi的Python接口。变量、约束、目标函数三部分分开写方便调试。核心代码如下注释里说明了每个变量和约束的作用import gurobipy as gp from gurobipy import GRB dt 3600 # 时段长度秒 T 24 S len(probs) N 2 m gp.Model(hydro_pv_short_term) # 阶段一变量水库运行计划不依赖场景 V m.addVars(N, T 1, lbV_min, ubV_max, nameV) # 库容 Q m.addVars(N, T, lb0, ubQ_max, nameQ) # 发电流量 SP m.addVars(N, T, lb0, ubSP_max, nameSP) # 弃水流量 Ph m.addVars(N, T, lb0, ubPh_max, namePh) # 水电出力 # 阶段二变量光伏消纳与弃光依赖场景 P_use m.addVars(T, S, lb0, nameP_use) # 实际消纳光伏 P_shed m.addVars(T, S, lb0, nameP_shed) # 弃光功率 # 水量平衡约束 for i in range(N): for t in range(T): inflow I[i, t] if i 0: tau travel_time[i] if t - tau 0: inflow Q[i - 1, t - tau] SP[i - 1, t - tau] m.addConstr(V[i, t 1] V[i, t] (inflow - Q[i, t] - SP[i, t]) * dt) # 初始库容与末库容 for i in range(N): m.addConstr(V[i, 0] V0[i]) m.addConstr(V[i, T] Vend_min[i]) m.addConstr(V[i, T] Vend_max[i]) # 光伏消纳与弃光 for t in range(T): for s in range(S): m.addConstr(P_use[t, s] P_shed[t, s] P_pv[t, s]) m.addConstr(Ph[0, t] Ph[1, t] P_use[t, s] P_grid[t]) # 目标函数最大化期望可消纳电量并加入弃光惩罚 obj gp.quicksum( probs[s] * (Ph[0, t] Ph[1, t] P_use[t, s] - penalty * P_shed[t, s]) for s in range(S) for t in range(T) ) m.setObjective(obj, GRB.MAXIMIZE) m.optimize()这个结构里有一个细节值得注意水量平衡中的时滞处理用了条件判断if t - tau 0。上游在调度期一开始放的水要经过 ( tau ) 小时才能到下游调度期前几小时下游应该收不到上游放水。如果忽略这个时间差把上游所有时段的出库都直接加入下游当天水量平衡会导致下游水量被重复计算。3.4 求解参数设置与结果导出Gurobi求解器默认参数在中小规模问题上表现不错但短期调度有时需要限制求解时间特别是在滚动优化场景下。我一般会加这样几个参数m.Params.TimeLimit 120 # 最长求解120秒 m.Params.MIPGap 0.02 # 允许2%的gap m.Params.Threads 4对于纯线性规划模型这组参数基本不会触发TimeLimit几秒内就能收敛如果模型因为加入整数变量变成了MIP这两个参数就是限制求解时长的保险丝。结果导出我习惯直接落盘成CSV包含每个时段的水库库容、发电流量、弃水流量、水电出力以及每个场景下的光伏消纳、弃光功率。落盘之后用pandas做汇总统计比直接从Gurobi变量对象里取值要方便得多。4. 算例结果分析看数字怎么验证模型的物理逻辑模型跑通只是第一步结果合理性验证才是复现真正收获的地方。我把算例结果整理成表格逐项分析了模型输出的物理逻辑这一步能帮你发现隐藏的建模错误。4.1 调度结果总览与时段级拆解下面是我在测试算例里得到的部分时段结果可以清晰地看到梯级水光互补的运行规律时段光伏预测 (MW)水电总出力 (MW)系统外送上限 (MW)弃光 (MW)上游库容变化08:002603806400缓降11:0045032070070缓升14:0051028071080缓升18:001805207000快速下降21:00605806400下降正午光伏出力达到峰值时水电出力明显压减上游水库库容上升相当于把水蓄起来傍晚光伏快速退坡时水电出力大幅增加水库放水发电补上缺口。这个“光伏高时蓄水、光伏低时放水”的模式正是梯级水光互补调度最理想的运行状态。如果算出来的库容过程是白天猛放水、晚上蓄水那大概率是目标函数或者约束写反了需要回头检查系统功率平衡方向。4.2 梯级水库调节过程怎么看判断梯级联动是否正确的关键指标是上下游水库的库容变化时序。模拟结果中下游水库在上午时段收到上游水库前一天傍晚放水的补充库容呈现缓慢上升趋势下午上游开始蓄水、放水量减少下游水库入库减少库容转为下降。这个过程与水流传播时间 ( \tau ) 的设置完全对应。还有一个验证点是弃水行为。梯级调度中如果上游来水较大而上游库容不足或下游无法承接会出现弃水。但模型应该优先利用水头高、耗水率低的水库发电而不是让水白白流走。测试算例里弃水主要集中在深夜来水高峰时段且下游库容接近上限时这正是物理上合理的弃水场景。4.3 场景削减对结果的影响我专门对比了三种方案只使用点预测、使用200条原始场景、使用削减后的5条场景。结果能清楚显示场景数量与优化质量的关系方案场景数目标值 (万kWh)求解时间 (s)弃光率点预测115200.812.4%原始场景200149638.510.1%削减场景515031.210.4%点预测方案目标值最高是因为模型把所有光伏预测都当作确定值敢于在系统功率平衡约束的边缘试探但一旦实际光伏低于预测就会大量弃水或弃光。200条原始场景的目标值最低因为极端场景拉低了期望值但也因此更保守、弃光率最低。削减后的5条场景在目标值和求解时间之间取得了很好的平衡。这个对比能回答复现时最常见的问题“场景削减有必要吗”答案是有必要但场景数不是越多越好。对短期调度来说5到10条典型场景就能覆盖主要的不确定性特征再多只增加计算负担。5. 复现路上的坑与调试经验这一章是全文最值钱的部分。模型和代码最终能跑通很大程度上是因为我把每个坑都详细记录并逐一解决了。每个坑背后都对应一条建模或编程上的通用经验换到别的调度模型同样适用。5.1 不可行模型先查单位再查约束冲突第一次求解直接报Model is infeasible我当时第一反应是约束写错了把水量平衡、库容上下限来回检查了好几遍全都没有问题。最后发现问题出在单位混用上游来水用的单位是 ( m^3/s ) 乘以 ( dt ) 得到了立方米而库容上下限用的却是论文里的亿立方米两个量级直接相差 ( 10^8 )模型当然无解。这个经历想说明的是排错要有顺序。遇到不可行第一步查所有参数的物理单位是否一致把所有流量统一成 ( m^3/s )把所有水量统一成 ( m^3 ) 或 ( 10^4 m^3 )。第二步查边界约束是否有交叉冲突比如库容下限过高加上初始库容过低模型连可行解都不存在。第三步才去查约束内部的等式和不等式方向。Gurobi提供一个很有用的工具模型不可行后可以调用m.computeIIS()找出不可行约束的最小集合。复现调试阶段用好这个功能能节省大量时间。5.2 场景削减要保留物理极值KMeans做场景削减有一个潜在风险聚类中心会被“均值化”极端场景的特征在削减后可能被抹平。比如光伏出力特别低、可能引起系统功率缺口的场景如果被聚类中心平均掉了模型得到的调度计划就缺乏对极端情况的预判能力。我的处理办法是在削减前先保留几条极端场景固定参与优化不允许聚类把它们的特征平滑掉。具体做法是把极端场景单独挑出来赋予一个固定概率然后再对剩余场景做KMeans。这样模型既保留了不确定性概率分布的典型代表又不丢失最危险的情况。5.3 分段线性化的分段点密度不是越密越好分段线性化处理水电出力特性时我一开始把发电流量分成了20段想知道精度能提升多少。结果模型直接变成了混合整数规划变量数暴增求解时间从几秒变成了十几分钟而目标值只变化了不到0.1%。后来我把分段数调到5到8段求解时间回到1秒以内精度损失完全可以接受。原因是短期调度中水电出力在大多数时段都在低流量段运行高流量段的分段点几乎用不到堆太多分段只会增加无意义的计算量。分段点的分布也值得注意流量较小的区间分段要密一些大流量区间可以稀疏一些这样拟合出来的出力曲线更贴近真实特性。5.4 求解精度与速度的权衡Gap和TimeLimit短期优化调度模型如果被投入到日内滚动运行中每次求解的时间窗口可能只有几分钟不可能无限等待最优解。这时候Gurobi的MIPGap参数就很有用。我把gap设置为2%求解器一旦找到足够接近最优的解就会提前停止大大缩短了求解时间。对于纯线性规划模型这个设置不产生什么影响因为线性规划求解本身很快。如果模型引入了0-1变量比如机组开停机状态MIPGap就非常关键了。实际运行时可以把gap从5%开始尝试观察目标值变化逐步收紧到1%或2%找到一个速度和精度都能接受的平衡点。5.5 随机种子必须固定这是我最后补的一个细节但重要性不亚于前面任何一个坑。场景生成和KMeans削减都有随机性如果每次运行结果都不一样后续分析根本没法做。我在所有使用随机数的地方都固定了种子包括numpy的np.random.seed(42)和KMeans的random_state42。这样同一套数据每次跑出来的结果完全一致也方便在校验模型时对比不同参数下的变化。另外一个相关的问题是论文复现时如果对方没有给随机种子你换一台机器或者换一个Python版本生成的场景都可能不同。所以我自己复现时一般不依赖远程随机数而是把生成好的场景集直接保存成CSV作为固定的输入文件。这样场景数据稳定后续所有调试都基于同一套数据展开避免因为场景变化来回折腾。最后再分享一个小经验。整个模型复现过程中最耗时间的不是最终写代码那一步而是想清楚“每个变量到底属于第一阶段还是第二阶段”。这个判断直接影响变量要不要带场景下标一旦搞反目标函数表面上还能求出值但结果会完全违背物理规律。我建议调试时先跑一个单场景版本的模型确认水电发电和光伏消纳的时序曲线符合直觉再加场景削减功能。单场景跑通多场景只是重复乘上概率权重而已出问题的概率会小很多。