搞电力调度和运筹优化的朋友对“梯级水光互补系统”这几个字应该不陌生。梯级水电站一条链串下来加上一片光伏共同外送电力这是典型的多能互补场景。最近我在复现一篇关于该系统的EI论文核心目标很直接在短期日内尺度上通过优化调度让电网尽可能多地消纳水光发电量而且要把光伏出力不确定性带来的风险装进模型。整个复现过程用的是Python跑通后发现这类模型最大的难点其实不在数学表达而在把物理约束、随机场景和求解器三者揉在一起。这篇博文就把我的复现思路、模型搭建、Python代码骨架以及踩过的坑完整记录下来适合正在做电力系统优化调度、论文复现或者课程设计的同学参考。1. 先把问题搞清楚梯级水光互补短期调度到底在做什么1.1 梯级水光互补系统的物理结构理解模型之前先得看懂物理对象。所谓梯级水电站就是一条河流上从上到下建了多座水电站上一级的出库流量会进入下一级的入库。比如上游A站发完电的水经过一段时间滞后变成下游B站的来水。这种水力联系让调度问题天然带上了“时序耦合”的特征不能像独立电站那样各管各的。在这个基础上再叠加光伏电站就出现了一个有意思的互补特性光伏出力集中在白天和水电站的调节能力形成搭配。水电站可以快速调节出力用来填补光伏波动带来的缺口光伏则在白天多发让水电站可以少放水、多存水留到晚上或阴天再发力。整体外送时还要考虑联络线传输功率上限、弃电惩罚、水库水位约束等一系列现实边界。如果只是做“多发电”问题就退化成了简单的水库调度。而论文里讲的是“最大化可消纳电量期望”就意味着电网侧可能吃不下那么多电。光伏大发的中午如果外送通道容量不够就会有弃光水库水位太高又遇上来水偏大可能被迫弃水。调度模型的任务就是在各种约束之间找平衡点让最终被电网真正消纳的电量期望值最大。1.2 “最大化可消纳电量期望”这句话怎么读这句话的关键词有两个一个是“可消纳电量”一个是“期望”。可消纳电量不是发电量发出来但被弃掉的部分不算有效贡献。目标函数里真正要最大化的是“上网电量”或者“负荷侧被满足的电量”。在送出通道受限的场景里这个区别极其关键它直接影响优化模型的结构。再说“期望”。光伏出力有随机性预报不可能完全准。今天中午的辐照度到底是1000瓦每平方米还是600模型不能拍脑袋定死。处理不确定性的标准做法是场景法把历史数据或预测误差分布抽样成若干个光伏出力场景每个场景算一次电量最后按概率加权平均得到的就是“期望可消纳电量”。这样调度方案不是只对某一个确定性预测最优而是对所有可能出现的情况整体上最优。所以这个模型的本质是一个含随机变量的短期调度优化问题。目标函数带期望算子决策变量是水电站出库流量、出力计划和光伏消纳计划约束覆盖水力、电力、通道多个维度。1.3 短期优化调度模型的边界与目标短期调度一般指未来24小时到72小时的调度计划时间颗粒度通常取1小时或15分钟。相比中长期水库调度短期模型更重视出力过程、通道约束和机组启停等运行细节相比实时调度它又多了一层“提前安排”的属性。模型求解结果可以作为实际运行的参考计划也可以嵌入滚动优化框架里每过一段时间用最新预测重新计算。复现这类EI论文最重要的不是把公式抄一遍而是搞清楚作者的决策变量集合、目标函数构成、约束条件覆盖范围以及不确定性怎么建模。这四个要素一旦确定模型框架就清晰了。Python代码只是把框架翻译成机器能算的形式。我复现时用的方法是先把论文里的符号表全部列出来再一项项对应到数组和矩阵里最后用线性规划或者混合整数线性规划求解。2. 数学模型目标函数与约束体系怎么写清楚2.1 目标函数期望可消纳电量最大化先给出一个简化的数学表达便于后面对照代码。设调度周期T为24个时段水电站编号为ii1,2,...,I光伏电场编号为j。决策变量包括水电站i在t时段的出力P_h(i,t)、光伏电场j在t时段的出力P_pv(j,t)最大值受场景限制、弃电量P_curt(t)、弃水量S_spill(i,t)等。目标函数可以写成最大化 Z Σ_t [ Σ_i P_h(i,t) Σ_j P_pv(j,t) ] × Δt − λ × Σ_t [ P_curt(t) Σ_i S_spill(i,t) ] × Δt这里的λ是惩罚系数用来让优化器尽量不要弃水弃光。严格地说“可消纳电量最大化”应该直接写成最大化Σ(P_h P_pv_used)而P_pv_used不能超过光伏理论出力。用惩罚项的好处是让模型在“少发”和“多发但被弃”之间自动做取舍比硬约束更灵活。如果采用场景法目标函数就要写成加权和形式最大化 Σ_s p_s × Σ_t [ Σ_i P_h(i,t,s) Σ_j P_pv_used(j,t,s) ] × Δt其中p_s是场景s的概率所有场景概率加起来等于1。注意场景引入之后水电站出力这类“可调度”变量通常仍然要求不随场景变化即所谓“非预期性”约束否则模型会拿未来信息作弊。这一点很多新手容易忽略。2.2 约束条件水量、库容、出力与通道约束条件是小水电联合优化的重头戏。我按类别列一个清单后面代码也是按这个顺序写约束的。水量平衡约束是梯级耦合的关键。水库i在t1时段的库容等于t时段库容加上入库流量减去出库流量。入库流量包含上游水库出库流量的滞后项和区间入流。公式写出来是V(i,t1) V(i,t) [ R_in(i,t) Σ_{u∈上游} Q_out(u,t−τ) − Q_out(i,t) ] × Δt这里τ是水流滞后时间区间入流R_in是外生给定值。这条约束把梯级电站串成了一条链写代码时要用矩阵把上下游关系表达出来不然很容易张冠李戴。库容约束很简单每个水库的水位或库容要保持在最小和最大之间。如果论文里用的是水位还要通过库容曲线转换成库容。复现时建议统一用库容立方米而非水位因为库容变化才是线性的。出库流量约束同样有上界和下界。下界可能是生态流量要求上界是机组最大过流能力。另外水电站出力不是独立的它和发电流量、水头有关。精确表达式是非线性的P_h η × ρ × g × Q × H。线性化处理通常有两种一是假定水头变化不大直接用“出力发电系数×发电流量”的线性关系二是分成多个水头区间分别建立线性约束。简化复现时用前者就够但要在论文里注明这是近似。联络线通道约束是造成弃电的根源之一。所有电站的出力之和减去本地负荷再扣除弃电不能超过联络线传输能力。这种约束让模型具备了“外送瓶颈”的物理含义。最后还有最小技术出力约束火电场景里常见水光系统里对应的是水电机组振动区、光伏逆变器最低功率等限制。如果论文没写可以不加但这会让结果偏乐观。2.3 随机性建模场景法场景法在代码里体现为对光伏出力数组的批量扩展。假设我们预测24小时光伏出力曲线是μ_pv(t)预测误差服从一定分布那么可以生成N个场景P_pv(j,t,s) max(0, μ_pv(t) ε(t,s))其中ε(t,s)采样自正态分布或历史误差分布。为了保证期望计算的合理性场景数不能太少我一般生成200到500个场景。如果模型规模大还可以用场景缩减技术比如K-means聚类挑选代表性场景。但复现论文时先别做花活直接用Monte Carlo抽样加固定随机种子就能对齐大多数论文结果。还有一个细节光伏场景生成时要设置上限装机容量超出部分截断掉否则可能出现等于1.2倍装机的离谱场景最终优化结果也会被带偏。我在复现时吃过这个亏后面调试章节会详细说。3. Python代码实现核心模块与求解思路3.1 数据准备和场景批量生成我把整个项目拆成四个模块数据模块、模型构建模块、求解模块、结果分析模块。数据模块负责读取水库参数、光伏预测曲线、负荷预测曲线并生成随机场景。核心数据结构用字典或dataclass保存都可以。我的习惯是这样组织# water_reservoir.py dataclass class Reservoir: name: str v_min: float # 最小库容 m^3 v_max: float # 最大库容 m^3 v0: float # 初始库容 q_min: float # 最小出库 m^3/s q_max: float # 最大出库 m^3/s k: float # 发电系数线性化 upstream: list # 上游水库编号 delay: int # 水流滞后时间 h光伏场景生成部分我用的方法比较朴素但很稳读取预测曲线设定预测误差标准差比如预测值的20%用numpy随机抽样生成N个场景。注意把随机种子固定住方便别人复现你的结果。import numpy as np def generate_pv_scenarios(prediction, capacity, n_scenes, seed42): rng np.random.default_rng(seed) n_hours len(prediction) scenarios np.zeros((n_scenes, n_hours)) for s in range(n_scenes): # 简单乘性误差模型 scale rng.normal(1.0, 0.15, n_hours) power prediction * scale power np.clip(power, 0, capacity) scenarios[s, :] power return scenarios这里用了乘性误差比加性误差更贴合光照强度的物理特征——晴天和云遮天的波动幅度是成比例的。每个场景概率相等即1/N。如果论文里给的是离散概率分布就把p_s数组作为参数传进来。3.2 用scipy.optimize.milp搭建线性规划骨架模型规模不大时用scipy.optimize.linprog就够了。但很多EI论文里的模型带二进制变量比如机组启停、弃电状态切换这时候就必须用scipy.optimize.milp。milp函数在SciPy 1.9.0之后才正式可用所以先确认环境版本。我在代码里优先用milp是因为它支持整数变量而且接口比linprog更贴近“决策变量数组约束矩阵”的标准写法。调用方式是这样的from scipy.optimize import milp, LinearConstraint, Bounds # 定义约束矩阵 A右侧 ub constraints LinearConstraint(A, lb, ub) bounds Bounds(lb_decision, ub_decision) result milp(c, constraintsconstraints, integralityintegrality, boundsbounds)这里最花时间的是组装约束矩阵A。模型决策变量一多手动拼矩阵往往要到几百列非常容易错位。我的经验是提前用一个字典记录每个变量名对应的列索引再按约束类型循环填矩阵。这样即使后期加了变量索引也不会乱。3.3 核心代码解读决策变量、目标、约束下面给一个简化但完整的模型骨架。假设只有1个水库和1个光伏电站24小时周期N个场景。先定义决策变量的排列方式# 决策变量顺序 # 时段 t: 0..23 # 场景 s: 0..N-1 # 每个 (s,t) 有 # 水电站出力 P_h 索引 base 0 # 光伏消纳 P_pv 索引 base 1 # 弃电量 P_curt 索引 base 2 # 水库出库 Q_out 索引 base 3 # 库容 V 索引 base 4列数等于场景数 × 时段数 × 5。如果场景数200时段数24就有24000列矩阵规模还算可控。约束按行追加水量平衡约束、库容上下限、出库上下限、光伏消纳上限、通道约束。目标函数系数c要写成最小化形式所以把“负的可消纳电量”作为目标向量。弃电和弃水用很小的惩罚系数确保只在必须时才发生。写起来是这样c np.zeros(n_vars) # 期望可消纳电量最大化 最小化 -E[P_h P_pv_used] # 因为每个场景概率相等直接对场景求和再除以 N for s in range(N): for t in range(T): idx_h var_index(P_h, s, t) idx_pv var_index(P_pv, s, t) idx_curt var_index(P_curt, s, t) c[idx_h] -1.0 / N c[idx_pv] -1.0 / N c[idx_curt] 0.01 / N # 弃电惩罚水量平衡约束是最容易写错的地方。对场景s、时段t、水库i而言库容关系是V[i, t1] V[i, t] inflow[i, t] Q_out_upstream[i, t] - Q_out[i, t]注意把单位统一。库容用立方米流量用立方米每秒时段长度用秒。我一般先把水位库容曲线拟合成线性函数然后把所有约束写到同一个矩阵里。约束矩阵的行数很多但列数固定构建过程可以用列表逐行append最后再转成稀疏矩阵。对大规模模型推荐用scipy.sparse保存否则内存会爆。4. 实操过程从零跑通一个简化算例4.1 算例说明与参数设置我复现时用的是一个简化算例一座上游水库加一座下游水库光伏电站额定容量300兆瓦联络线外送上限400兆瓦本地负荷曲线取典型夏季日负荷。光伏预测曲线用一条峰值在中午12点的钟形曲线预测误差标准差取预测值的15%。场景数取100随机种子固定为42。水库参数设置如下参数上游水库下游水库最小库容百万m³8050最大库容百万m³500300初始库容百万m³300150最小出库m³/s00最大出库m³/s600800发电系数MW·s/m³0.0050.004滞后时间h10之所以把滞后时间设为1小时是为了体现梯级耦合的效果上游放水后1小时才进入下游。这样调度计划不会出现“下游先看到未来水再决定现在出力”的穿越现象。4.2 结果输出与合理性检查跑完优化后我习惯先画三张图第一张是水电站出力与光伏消纳的时序曲线第二张是水库库容变化曲线第三张是弃电量柱状图。如果结果合理光伏大发的中午时段水电站出力会明显压低给光伏让路傍晚光伏出力掉下去水电站再把库容里的水放出来发电填晚高峰。还有两个关键数值要核验。第一弃电量必须是“被迫”产生的也就是通道容量或者库容约束被卡住时才会出现不会平白无故弃电。第二水库末库容要落在允许范围内最好还和初始库容差别不大这样模型才有滚动调度的实用价值。如果发现弃电量全是0先别高兴。那可能是惩罚系数设太小模型宁愿完全不弃也不去多发电反过来如果惩罚系数设太大弃电量又会异常高。我用的惩罚系数不超过目标函数中电量系数的2%实测下来比较稳。4.3 求解器选型与性能调优用scipy自带的HIGHS求解器处理30000列以内的问题没有问题但规模一大瓶颈就来了。我在200个场景、24个时段、5个电站的模型里变量数已经超过15万约束数超过8万这时内存占用明显上升求解时间也从几秒跳到了几分钟。这时候有几个优化思路特别值得尝试。第一是用场景缩减技术把200个场景聚类成20个代表场景变量数直接砍掉90%解的质量损失在1%以内。第二是约束矩阵用稀疏格式存储scipy.sparse.lil_matrix逐行填充再转成csr_matrix传给milp。第三是如果论文允许把部分整数变量松弛成连续变量MILP就变成了LP求解速度快几个量级。我这里要特别强调一个经验不要依赖默认参数。milp函数里有一个options参数可以设置求解时间上限和gap容忍度。复现时间有限时我会设置options { time_limit: 300, mip_rel_gap: 0.01, }这样即使问题规模很大也能在5分钟内拿到一个gap不超过1%的可行解。EI论文复现的核心是验证模型逻辑不是追求毫秒级收敛设个合理阈值就行。5. 复现路上最值得记录的踩坑清单5.1 线性化陷阱水头-出力关系怎么处理很多论文会用分段线性函数近似水电站出力特性出力发电系数×发电流量×净水头。如果把水头当成常数模型是线性的求解快但精度差。如果必须考虑水头变化又不想引入非线性我有一个折中方案把库容范围分成K个区间每个区间单独用线性关系然后引入二进制变量选择所在区间。这个做法的缺点是需要额外增加约束代码量比较可观。复现时如果论文没重点讨论水头变化建议直接用常数水头。我在第一个版本里为了“忠实复现”强行加了水头分段结果求解时间从几十秒涨到半小时还经常不收敛。后来对比发现对常规日调度来说水头变化对可消纳电量期望的影响通常小于3%不值得为此付出这么大的计算代价。5.2 别把梯级上下游关系写错这个坑说出来有点丢人但确实发生了。我在组约束矩阵时上游出库对下游入库的贡献项方向写反了导致下游水库莫名其妙多了几倍水。排查了三天最后是画了一张“水量从哪里来到哪里去”的流图文字版才发现的。给大家一个自查技巧随便挑一个时段手动计算某水库的库容变化和代码输出的结果对比一下。具体做法是取一个简单场景比如光伏出力全为0水电按最大出库放水然后看下游水库库容是否按预期增涨。如果一小时内库容增量和上游放水量×3600秒相等说明时序耦合关系写对了一半。另一半是滞后项的索引偏移也要单独核对。5.3 代码工程配置的几点小建议复现这类项目Python环境配置的坑比算法本身还多。我的建议是先用一个干净的虚拟环境不要一股脑装一堆库。核心依赖只有numpy、scipy、pandas和matplotlib。如果要用场景缩减加scikit-learn。有个问题值得提醒scipy.optimize.milp是较新版本才有的函数如果你用的是老环境import会直接报错。我第一次跑就遇到这个情况解决方案是升级scipypip install --upgrade scipy如果网络条件限制太多可以先把whl包下载好再离线安装。和pulp、cvxpy这些库相比scipy的接口比较底层但足够应付标准的线性约束问题。cvxpy语法更优雅不过当变量规模特别大时我用下来的感觉是scipy更直接一些多花一点组装矩阵的时间换来更少的性能损耗。6. 复现EI论文的通用方法论6.1 一篇论文拿到手先拆这几个部分复现EI论文的第一步不是写代码而是拆解论文结构。我拿到文章后会分成三块读一是问题描述和符号定义二是数学建模公式三是算例设置与结果图。符号定义这一块尤其关键。EI论文的符号往往不统一有的用P表示出力有的用N表示机组数还有的用π表示场景概率。我的习惯是花两小时把所有符号整理成对照表再用统一的命名映射到代码变量。这一步很枯燥但省下的调试时间远超花掉的时间。算例设置也必须看仔细。很多复现结果对不上不是因为模型写错而是因为初始水位、机组极限、通道容量这些边界条件和论文不一样。论文正文没写的参数我会先按工程常识补上然后在记录里注明“假设值”避免误导后人。6.2 从论文公式到Python代码的翻译技巧公式翻译成代码核心是矩阵化思考。看到求和符号先想清楚“这题是求矩阵里哪一轴的和”看到下标t和s就想到二维数组的轴顺序。水电站梯级调度里最常用的是三维数组维度分别是场景、时段、电站编号写变量索引时要统一顺序。我的翻译三板斧是第一把公式里的所有变量名列出来去重后分配索引第二把所有等式约束用“左边减右边等于0”的格式整理第三把不等式约束统一写成“左边减去右边≤0”。这样做完之后组装矩阵就是一个机械化过程。还有一点特别实用先写一个小规模确定性模型比如把场景数降到1验证模型和论文结果是不是对得上。验证通过后再把这个模型封装成函数传入场景数组。这样能把“随机性实现有问题”和“模型本身有问题”两类错误分开。6.3 我习惯的复现顺序和验证节奏复现顺序上我推荐六个步骤参数准备、确定性模型搭建、单场景求解、多场景扩展、结果可视化、敏感性分析。这六步每步都有独立的验证点不要跳步。单场景求解这步尤为划算。只有1个场景时期望就等于那个场景的值优化结果应该非常接近人工经验判断。比如光伏中午大发、负荷不高的时候水电站出库应该减少晚高峰负荷上来水电站出库应该增加。如果单场景结果都说不通就别急着往多场景推。多场景扩展后重点看两点一是各个场景的弃电位置是否合理二是调度方案的鲁棒性是否体现。情景里阴天光伏弱的日子水电机组应该早点发力晚上不至于把库容放空。如果所有场景下的出力计划都一模一样说明不确定性没有真正进入模型大概率是某处下标写错了。敏感性分析是最后一步也是论文复现最有价值的一环。我会把联络线容量从300逐步调到500兆瓦观察可消纳电量期望的变化曲线。通道容量低时曲线饱和得很厉害说明瓶颈在外送通道容量提升到一定程度后曲线增长放缓说明瓶颈转移到了发电能力。这条曲线画出来模型的整个物理逻辑就通了。最后分享一个我自己的小习惯每复现一篇论文我都会把代码里的关键参数抽到一个配置文件里用dataclass或者json统一管理。这样面对不同算例时改参数十分方便也方便日后继续在工程里使用。这个模型的后续扩展方向也很多比如加入储能、考虑机组组合、把确定性优化转成鲁棒优化只要基础框架搭好往上加东西会非常顺手。