梯级水光互补系统最大化可消纳电量期望短期优化调度模型这串词我在工位前盯了整整两天才敢说自己读懂了。论文要解决的问题其实一句话就能讲清楚在光伏出力预测不准、电网消纳能力有限的前提下怎么让上下游一串水电站配合光伏发电使未来24小时内能被电网实际吸收的电量达到最多。真正动手用Python复现这个EI模型时我发现论文里每个词都是坑。梯级意味着上下游电站之间水流有段时间滞后可消纳电量期望要求把光伏不确定性做成几十个场景再求平均短期优化调度则对求解效率提出了硬要求。这篇博文会把我从模型推导、代码实现到踩坑的经历完整过一遍。如果你是正在复现新能源调度方向论文的研究生或者工作中需要做水光互补优化的工程师照着这篇的思路能节省大量试错时间。1. 为什么复现这个模型梯级水光互补的核心矛盾1.1 可消纳电量期望比多发电更现实的目标很多刚接触这个方向的人会把目标理解成让系统发最多电其实这是个常见误区。电网要的不是机组能发多少而是能接收多少。光伏午间大发时如果负荷和输电断面吃不下多余的电只能弃掉水电也有类似情况上游一直放水但下游通道受限就会出现弃水。所以论文里强调的可消纳电量本质上衡量的是发电量与电网接纳能力的匹配程度。那为什么还要加期望两个字因为光伏出力不是一个确定值。气象预报有误差光伏实际出力可能比预测高也可能比预测低。如果把预测值当作真实值来做调度一旦实际偏差大调度方案就不够用了预测偏高会过度放弃光伏水电却占了通道预测偏低又会留下太多水电调节空间白白牺牲发电量。期望目标的做法是把各种可能的光伏出力都生成对应场景给每个场景一个概率最后把各场景的可消纳电量按概率加权作为方案的评价标准。这个目标能保证调度方案在各种可能的世界里平均表现最优而不是只在某个预测场景下最优。我打过一个比方像逛超市时货架库存有波动你的预算电网容量是固定的你要决定购物车各电源出力怎么装。如果照着预计库存拿货实际库存和预计不一致时你买到的需要的商品就少了。期望模型相当于同时考虑了库存偏高、库存偏低、库存正常三种情况选的购物车方案是三种情况平均收益最好的那一个。1.2 梯级水电的蓄与放互补的关键在于水量时序耦合水光互补听起来简单光伏大发电时水电少发、光伏没电时水电多发但一旦加上梯级两个字事情就变得复杂了。梯级电站是沿着同一条河流从上到下排列的多个水库上游电站放出去的水经过一段时间流到下游电站成为下游电站的来水。这意味着你的调度决策不能只看单个电站上游的一举一动都会影响下游未来几个小时的可用水量这个影响还有时间延迟。水库本身就是天然的储能装置。有调节能力的水库在光伏出力高的时段可以把水存住、压低水电出力给光伏让出消纳空间到了晚上光伏退出再加大放水发电。这个蓄放过程放在梯级系统里还需要考虑每个水库的库容上下的约束以及出库流量对下游入库流量的传递。上游如果为了短时补光伏缺口猛放水下游可能在几个小时后遭遇来水集中到顶被迫弃水等于把上游省下的电量又丢掉了。所以梯级水光互补的优化不是一个简单的功率分配问题而是一个带时序耦合的水量调度问题。理解了这个背景再看论文标题里的短期优化调度就知道模型要做的核心事情是什么在给定未来24小时的径流预报和光伏预测场景下求出每条河流每个时刻的蓄水量、发电流量、光伏消纳量使得期望可消纳电量最大。下面就把这个模型的数学骨架一层层拆开。2. 模型框架拆解目标函数、约束条件和不确定性建模2.1 目标函数把期望落到实处论文的目标函数可以写成这样max Σ_{s∈S} π_s · Σ_{t1..T} [ Σ_{i∈R} P_hydro(i,s,t) P_pv_used(s,t) ] · Δt其中s是光伏场景编号π_s是该场景对应的概率T是调度总时段通常24小时R是梯级水库集合。P_hydro是水电出力决策变量P_pv_used是光伏实际消纳量决策变量。这里有一个非常关键的理解点P_pv_used不是光伏预测出力而是决策变量。它的上界是场景s下时段t的光伏可用出力目标函数鼓励它尽可能大但约束会限制它不能超过电网消纳能力。换句话说模型允许弃光但弃多少、在哪里弃要由优化器来选择。很多论文复现者第一步就栽在这里——把P_pv_used写成等于预测值相当于把模型从优化光伏消纳变成了固定光伏消纳水电的调节作用直接被废掉一大半。后面我会单独用一节讲这个坑。目标函数里的Δt需要统一量纲功率用MW、时段按小时算乘起来就是MWh这样期望可消纳电量就是MWh级别的一个数值方便和论文结果对比。别小看这个单位统一好多复现代码表面上看能跑算出来的电量结果却大得离谱根源就是Δt混用了秒和小时。2.2 约束条件水量平衡是最容易错的一道关模型的基础约束是水量平衡方程每个水库每个场景每个时段都必须满足V(i,s,t1) V(i,s,t) Δt_water · ( I(i,s,t) Q_up(i,s,t-τ_i) - Q(i,s,t) - Q_spill(i,s,t) )每个量的含义我都标清楚V是水库蓄水量I是区间入流比如支流汇入和降雨Q_up是上游电站的出库流量经过时滞τ_i后到达本站Q是发电流量Q_spill是弃水流量。Δt_water如果以秒为单位流量用m³/s乘起来得到的是m³这样库容和流量的单位自洽。时滞τ_i是梯级模型的灵魂所在。上游出库到下游入库不是瞬时的水流在河道里流动需要时间不同河段可能是1小时到几小时不等。如果忽略这个延迟你会看到下游水位约束被毫无道理地打破因为上游今天放的水不可能立刻到下游你却当成已经到了。除水量平衡外还有几组常见约束库容上下限约束Vmin(i) ≤ V(i,s,t) ≤ Vmax(i)以及调度周期末库容约束短期调度通常固定或约束末水位否则模型会把水全部放光只追求电量不顾后续。发电流量约束Qmin(i) ≤ Q(i,s,t) ≤ Qmax(i)弃水流量非负。水电出力约束P_hydro(i,s,t) 与发电流量、水头相关。严格的水电出力是9.81·η·H·QH为发电水头Q为流量两者乘积形成非线性项。后面我会讲如何做线性化处理。光伏消纳约束0 ≤ P_pv_used(s,t) ≤ P_pv_avail(s,t)即消纳量不会超过该场景下能达到的光伏出力。电网消纳上限约束Σ_i P_hydro(i,s,t) P_pv_used(s,t) ≤ P_grid_limit(t)这条约束和期望目标共同定义了可消纳电量的物理边界。这里每条约束都要按场景展开因为每个场景下光伏不同相应的最优水量分配也会不同。场景之间唯一需要真正耦合的是水库的初末库容关系——起始库容是所有场景共用的末期库容也要所有场景共同满足这样调度计划才具备可执行性。2.3 不确定性场景从预测误差分布到场景集合有了目标函数和约束还缺场景从哪来。复现论文时通常分三步处理不确定性。第一步统计光伏预测误差。拿历史预测曲线和实际出力曲线做差得到误差序列然后拟合分布。不建议直接用标准正态分布因为正态分布两端无限延伸抽样会出现负的光伏出力明显不合理。更实用的做法是截断正态或者用beta分布拟合并限制在[0, P_installed]区间内。第二步生成带时序相关的抽样。光伏预测误差在相邻时段是有相关性的——如果10点实际比预测高11点大概率也偏高。独立逐时段抽样会得到锯齿状高频抖动的场景和真实偏差规律不符。最简单的处理是用历史误差序列求协方差矩阵做多元正态抽样生成完整的24小时误差曲线再加到预测曲线上得到每个场景的光伏出力曲线。第三步场景缩减。原始抽样可能抽500个甚至1000个场景但把这些所有场景全部放进模型变量和约束规模会大到求解器难以接受。通常用聚类方法把相似场景合并挑出50到100个代表性场景再按簇内样本数占比确定每个代表场景的概率。这一步后面代码部分会给出具体实现。3. Python代码实现求解器选型与建模避坑3.1 为什么我选Gurobi而不是纯开源方案先说结论复现这类随机优化调度模型我最推荐Gurobi。不是因为它最贵而是因为学术用户免费。用学校邮箱在官网申请学术license可以全功能使用没有变量数限制这对搞论文复现的人来说是零成本高收益。对比情况我列在下面求解器授权模式适用场景我的评价Gurobi闭源学术免费LP/MILP/QP大规模问题首选性能稳、API顺手CPLEX闭源学术免费与Gurobi同级也可以但建模层不太常用SCIP开源MILP小规模够用场景放大后偏慢HiGHS开源LP/QP、部分MILP纯LP不错MILP大规模吃力如果单位不允许安装商业软件可以考虑先用HiGHS跑小规模模型验证逻辑提交论文结果时再换Gurobi。但我要提醒一句场景数超过80、水库数超过3时HiGHS的求解时间会明显拉长而Gurobi通常还在秒级到分钟级。论文复现里你不是在写工业代码而是在和论文结果对照追求的是求解时间可接受、解的质量可复现Gurobi是最稳的选择。3.2 核心建模代码从变量、约束到目标函数下面是我整理出的代码骨架去掉了部分非必要细节重点保留和论文公式对应的结构。整体流程是定义参数 → 建变量 → 加约束 → 设目标 → 求解。import gurobipy as gp from gurobipy import GRB import numpy as np # ---------- 基础参数 ---------- N_RES 3 # 梯级水库个数 N_SCEN 50 # 光伏场景个数 T 24 # 调度时段小时 WATER_DT_SEC 3600.0 # 水量平衡用的时间单位秒 HOUR_DT 1.0 # 电量目标用的时间单位小时 tau [0, 1, 2] # 上游到本站的水流时滞小时 pi_s np.ones(N_SCEN) / N_SCEN # 场景概率均匀假设 # 水库参数示意值 V0 [120.0, 85.0, 60.0] # 初始库容百万m3 Vmin [60.0, 40.0, 30.0] Vmax [180.0, 110.0, 80.0] Qmax [250.0, 220.0, 200.0] # 最大发电流量 m3/s coef [0.65, 0.60, 0.55] # 简化出力系数 MW/(m3/s)隐含水头近似 # 时段数据示意实际从Excel/CSV读取 inflow np.random.rand(N_RES, N_SCEN, T) * 50 20 # 区间入流 m3/s pv_avail np.random.rand(N_SCEN, T) * 150 # 光伏可用出力 MW grid_limit np.full(T, 180.0) # 消纳上限 MW # ---------- 建立模型 ---------- model gp.Model(cascade_hydro_pv) # 决策变量 V model.addVars(N_RES, N_SCEN, T 1, lb0, nameV) # 库容 Q model.addVars(N_RES, N_SCEN, T, lb0, nameQ) # 发电流量 Qspill model.addVars(N_RES, N_SCEN, T, lb0, nameQspill) # 弃水流量 Ph model.addVars(N_RES, N_SCEN, T, lb0, namePh) # 水电出力 Ppv_use model.addVars(N_SCEN, T, lb0, namePpv_use) # 光伏消纳量 # ---------- 约束1水量平衡含时滞 ---------- for r in range(N_RES): for s in range(N_SCEN): for t in range(T): up 0.0 if r 0: idx t - tau[r] idx max(idx, 0) # t小于时滞时简化取0严格可用历史出库 up Q[r - 1, s, idx] * WATER_DT_SEC # 立方米 model.addConstr( V[r, s, t 1] - V[r, s, t] - ( inflow[r, s, t] * WATER_DT_SEC up - Q[r, s, t] * WATER_DT_SEC - Qspill[r, s, t] * WATER_DT_SEC ) 0, namefbalance_{r}_{s}_{t} )水量平衡这段是全文最需要细看的地方。我最初实现时把时滞索引写成了idx t - tau[r]没有加max(idx, 0)模型直接报索引越界。加了下限处理后又发现另一个问题t小于时滞时上游来水应该用调度开始前的历史出库值而不是0复现阶段我为了简化取0结果前几个时段的水位和论文对不上。后来改成读取历史实际出库记录作为边界才算完全对齐。继续加约束# ---------- 约束2库容上下限与末库容 ---------- for r in range(N_RES): for s in range(N_SCEN): for t in range(T 1): model.addConstr(V[r, s, t] Vmin[r], namefVmin_{r}_{s}_{t}) model.addConstr(V[r, s, t] Vmax[r], namefVmax_{r}_{s}_{t}) # 周期末库容回到设定值原文若用范围约束则改成不等式 model.addConstr(V[r, s, T] V0[r], namefVend_{r}_{s}) # ---------- 约束3发电流量上限 ---------- for r in range(N_RES): for s in range(N_SCEN): for t in range(T): model.addConstr(Q[r, s, t] Qmax[r], namefQmax_{r}_{s}_{t}) # ---------- 约束4水电出力简化线性化 ---------- for r in range(N_RES): for s in range(N_SCEN): for t in range(T): model.addConstr(Ph[r, s, t] coef[r] * Q[r, s, t], namefPh_{r}_{s}_{t}) # ---------- 约束5光伏消纳上界 ---------- for s in range(N_SCEN): for t in range(T): model.addConstr(Ppv_use[s, t] pv_avail[s, t], namefpv_{s}_{t}) # ---------- 约束6电网消纳上限 ---------- for s in range(N_SCEN): for t in range(T): model.addConstr( gp.quicksum(Ph[r, s, t] for r in range(N_RES)) Ppv_use[s, t] grid_limit[t], namefgrid_{s}_{t} ) # ---------- 目标函数期望可消纳电量最大化 ---------- obj gp.quicksum( pi_s[s] * HOUR_DT * ( gp.quicksum(Ph[r, s, t] for r in range(N_RES)) Ppv_use[s, t] ) for s in range(N_SCEN) for t in range(T) ) model.setObjective(obj, GRB.MAXIMIZE) model.optimize()跑完之后V.x、Q.x、Ppv_use.x就是各场景下的最优调度值。注意这里水量平衡里的库容单位是百万m³还是m³我代码里为了避免量纲混乱直接以m³为单位写实际工程数据如果是万m³需要先换算。代码写完后要逐行检查量纲我在调试时因为这个坑把库容约束放大了10000倍Gurobi报infeasible排查了半个小时才发现是V的单位和流量换算不一致。3.3 场景缩减用KMeans把500个场景压到50个场景生成完直接全部丢进模型大概率会卡。以500个场景、3个水库、24时段为例仅库容变量就有500×3×2537500个加上流量、出力、光伏变量总变量数轻松超过10万约束也是同等量级。虽然Gurobi硬吃也不是完全不行但求解时间会从几十秒变成几十分钟完全没必要。场景缩减是一个标准操作。from sklearn.cluster import KMeans # pv_scenes: shape (n_original_scenes, T)原始光伏场景曲线 kmeans KMeans(n_clusters50, random_state42).fit(pv_scenes) rep_scenes kmeans.cluster_centers_ # 50个代表场景 cluster_counts np.bincount(kmeans.labels_) weights cluster_counts / cluster_counts.sum() # 场景概率 # 用 rep_scenes 和 weights 替换原始场景其他约束、目标结构不变缩减后精度损失取决于原始场景的聚类程度。我的经验是预测误差分布平稳时500场景缩到50个期望目标值的变化通常在1%以内而求解时间可以从分钟级降到20秒级。如果论文里明确写了场景数那就按论文参数来但代码调试阶段建议先用50个场景跑通出结果后再放大验证。4. 完整的实现链路从数据准备到结果输出4.1 输入数据清单论文没给数据怎么破复现EI论文最难受的不是模型公式而是作者经常只给几张结果图具体流域参数、径流数据、预测误差参数一个都不给。我的应对方案是先把论文可用的参数抄下来比如水库个数、装机容量、预测误差标准差其余用公开典型数据补全。整理过一张数据清单照着准备不会漏数据项说明论文没给时的替代方案水库参数库容上下限、初始库容、最大发电流量、出力系数用该流域所在省份典型水电站公开参数近似径流与区间入流各水库各时段来水用历史水文站日平均流量按比例分摊到小时光伏装机与预测曲线24小时光伏可用出力用当地历史辐照度乘装机容量再乘转换效率预测误差分布均值、方差或协方差矩阵用装机容量的5%-15%作为标准差再调参对比电网消纳上限各时段允许接入的上限功率按负荷曲线或输电断面限额设定数据准备阶段有个原则先用看起来合理的数据把代码跑通再回头对照论文结果曲线微调参数。不要一上来就纠结某条径流曲线是否精确模型逻辑通了参数替换只是重跑一遍的事。4.2 迭代式开发先10个场景跑通再放大我强烈建议分三步走不要一步到位。第一步跑一个确定性模型即场景数设为1、光伏取预测值。这一步的目的是验证水量平衡、库容约束和出力约束没有写错。第二步跑10个场景验证随机逻辑——目标函数里场景概率是否生效Ppv_use是否合理。第三步再跑到50到100个场景出正式结果。每步都要加检查代码。我最常用的是把每个水库每个时段的库容打印出来肉眼扫一遍看是否在上下限之间再看看目标函数值和Gurobi求解日志里的约束数是否符合预期。如果模型报Infeasible优先检查三处末库容约束是不是与水量平衡矛盾、时滞索引是否越界、电网消纳上限是否小到所有电源同时开机都满足不了。4.3 结果曲线调度图和弃电率怎么画求解完成后先把结果写成CSV存档再画两张图。一张是各时段光伏消纳量、水电总出力、总出力与消纳上限的对比曲线另一张是各水库水位变化曲线。画图用matplotlib就够了import matplotlib.pyplot as plt pv_used_sol np.array([Ppv_use[s, t].x for s in range(N_SCEN) for t in range(T)]).reshape(N_SCEN, T) ph_sol np.array([Ph[r, s, t].x for r in range(N_RES) for s in range(N_SCEN) for t in range(T)]) ph_sol ph_sol.reshape(N_RES, N_SCEN, T) pv_avg pv_used_sol.mean(axis0) # 各时段光伏消纳期望 hydro_avg ph_sol.mean(axis1) # 各时段水电总出力期望 total_avg pv_avg hydro_avg plt.figure(figsize(10, 4)) plt.plot(range(T), pv_avg, labelPV used (expected)) plt.plot(range(T), hydro_avg, labelHydro (expected)) plt.plot(range(T), total_avg, labelTotal output (expected)) plt.plot(range(T), grid_limit, --, labelGrid limit) plt.xlabel(Hour) plt.ylabel(MW) plt.legend() plt.tight_layout() plt.show()指标计算就不能只看图了要按定义算期望可消纳电量是sum(pi[s] * sum(total[s,t] for t)) * HOUR_DT弃光率是1 - sum(Ppv_use) / sum(pv_avail)水能利用率可以定义为实际发电用水量与可发电用水量之比。算完这些指标再对照论文里给出的数量级比如弃光率是百分之几而不是百分之几十总消纳电量在几百MWh到几千MWh量级基本就能确认复现是否到位。5. 复现中踩过最深的三个坑5.1 水头-出力关系非线性模型直接做成了非凸论文里的水电出力公式常写成 P 9.81·η·H·QMWH是水头mQ是发电流量m³/s。问题出在H本身是库容的函数而库容又是变量于是出现变量×变量的乘积项整个约束变成非线性的。直接丢给Gurobi要么报unsupported要么把它当MIQP处理求解速度慢到怀疑人生。我在复现时采用了两段处理。第一阶段先用固定平均水头近似即把H当成常数P_hydro就变成了Q的线性函数整模型变成纯LP几十秒出解。这个近似在短期调度里是合理的因为一天内水库水位变化幅度通常不大水头变化对出力的影响属于二阶误差。第二阶段如果论文要求水头精确建模再对H-Q出力曲线做分段线性逼近把每个线性段用SOS2约束实现模型变成MILP速度会慢一些但仍是可解的。5.2 可消纳是决策变量不是预测值直接拿来用这个坑我在前面提过但值得再展开一次。最初实现时我把Ppv_use写死为预测值理由是光伏发多少我就消纳多少除非超过上限。跑出来的结果很离谱午间光伏预测80MW、电网上限100MW、水电可以发30MW的时候总出力变成110MW模型被迫弃掉10MW但水也没有存下来。为什么因为光伏出力被写死后水电只能被动跟着让路根本没有商量余地。改成决策变量后模型自己算出来最优解是水电22MW、光伏78MW总分刚好100MW电网容量被充分利用水库还多存了8MW水量的水晚上光伏归零时再放出来。这一改期望可消纳电量直接提升了几个百分点。判断你的模型写对了没有可以做一个简单测试把光伏可用出力设成零模型应该自动提高水电出力来补缺口再把电网上限设成很小模型应该自动降低光伏消纳而不是硬顶上限。如果这两个行为都正常说明决策变量的设计是合理的。5.3 论文公式和代码之间的落差用小案例反推验证EI论文的符号经常不统一同一篇里Q可能一会儿表示流量、一会儿表示电量时滞的定义也不写清楚。遇到这种情况我不建议直接对着大系统调试而是构造一个极简案例1个水库、1个时段、无光伏。这个案例你可以手算出最优解比如初始库容100、最大发电流量限制对应出力50MW那最优解就是发50MW末库容等于初始库容减去发电流量对应水量。把代码参数改成这个极简案例跑出来如果和手算一致说明水量平衡和目标是对的再逐步加入光伏、第二级水库、时滞每加一块都用行为是否符合物理直觉来校验。这样一层层叠上去基本能做到模型的每部分都有确凿依据。如果只是把大模型跑完看总结果一旦数字对不上你完全不知道是哪个约束错了排查成本反而更高。6. 验证与延伸随机模型到底带来多少提升6.1 确定性模型和随机模型的结果对比我复现时用同一套流域参数跑了两个模型一个只考虑光伏预测均值确定性一个按50个场景做随机优化下面是示意结果不是原论文数字只是展示对比方向指标确定性模型随机优化模型变化期望可消纳电量MWh428043150.8%弃光率8.7%6.1%-2.6个百分点水能利用率91.2%94.5%3.3个百分点提升看起来不大但关键不是数值而是来源。确定性模型在预测偏高时会浪费水电调节空间在预测偏低时又会占用光伏消纳通道随机模型看到全部场景后会在水库蓄放上做一个更稳妥的取舍——到底留多少库容应对光伏晚高峰到底放弃多少光伏给水电腾空间这些决策在全场景期望意义下是最优的。也正是这种牺牲个别场景、提升平均表现的特征让期望模型在实际运行中比单点预测更可靠。6.2 这个模型还能往哪里扩展代码结构搭好以后扩展方向其实很清晰。加储能只需要新增充放电变量和SOC荷电状态转移约束目标函数里加储能放电出力做现货市场优化把电量目标换成收益目标加入价格场景即可考虑日前-日内两阶段决策则要把模型改成两阶段随机规划第一阶段变量是水库蓄放计划第二阶段变量是光伏消纳和水电修正出力约束层级要重新组织。还有一个方向是分布鲁棒优化不假设光伏误差服从某个固定分布而是构造一个包含所有可能分布的模糊集目标改成最坏分布下的期望收益这在学术上更前沿代码也是在现有模型基础上加对偶约束。这次复现给我最大的体会是EI论文的模型骨架其实不难难的是把工程细节填对。目标函数就一行约束加起来也不超过十组但单位、时滞、决策变量自由度这些细节任何一个没想透求解器都会用一条infeasible或一个离谱数字来提醒你。最后分享一个实用技巧写这类随机优化模型先用10个场景跑通再逐步加到100个。如果10个场景都跑不动基本可以断定是约束规模写错了而不是求解器的问题用这个习惯调试能省下大把和求解日志死磕的时间。 
长河君 · 2026/10/6 4:52:27 · 阅读 1.2k 梯级水光互补可消纳电量期望最大化短期调度模型与Python实现
梯级水光互补系统最大化可消纳电量期望短期优化调度模型这串词我在工位前盯了整整两天才敢说自己读懂了。论文要解决的问题其实一句话就能讲清楚在光伏出力预测不准、电网消纳能力有限的前提下怎么让上下游一串水电站配合光伏发电使未来24小时内能被电网实际吸收的电量达到最多。真正动手用Python复现这个EI模型时我发现论文里每个词都是坑。梯级意味着上下游电站之间水流有段时间滞后可消纳电量期望要求把光伏不确定性做成几十个场景再求平均短期优化调度则对求解效率提出了硬要求。这篇博文会把我从模型推导、代码实现到踩坑的经历完整过一遍。如果你是正在复现新能源调度方向论文的研究生或者工作中需要做水光互补优化的工程师照着这篇的思路能节省大量试错时间。1. 为什么复现这个模型梯级水光互补的核心矛盾1.1 可消纳电量期望比多发电更现实的目标很多刚接触这个方向的人会把目标理解成让系统发最多电其实这是个常见误区。电网要的不是机组能发多少而是能接收多少。光伏午间大发时如果负荷和输电断面吃不下多余的电只能弃掉水电也有类似情况上游一直放水但下游通道受限就会出现弃水。所以论文里强调的可消纳电量本质上衡量的是发电量与电网接纳能力的匹配程度。那为什么还要加期望两个字因为光伏出力不是一个确定值。气象预报有误差光伏实际出力可能比预测高也可能比预测低。如果把预测值当作真实值来做调度一旦实际偏差大调度方案就不够用了预测偏高会过度放弃光伏水电却占了通道预测偏低又会留下太多水电调节空间白白牺牲发电量。期望目标的做法是把各种可能的光伏出力都生成对应场景给每个场景一个概率最后把各场景的可消纳电量按概率加权作为方案的评价标准。这个目标能保证调度方案在各种可能的世界里平均表现最优而不是只在某个预测场景下最优。我打过一个比方像逛超市时货架库存有波动你的预算电网容量是固定的你要决定购物车各电源出力怎么装。如果照着预计库存拿货实际库存和预计不一致时你买到的需要的商品就少了。期望模型相当于同时考虑了库存偏高、库存偏低、库存正常三种情况选的购物车方案是三种情况平均收益最好的那一个。1.2 梯级水电的蓄与放互补的关键在于水量时序耦合水光互补听起来简单光伏大发电时水电少发、光伏没电时水电多发但一旦加上梯级两个字事情就变得复杂了。梯级电站是沿着同一条河流从上到下排列的多个水库上游电站放出去的水经过一段时间流到下游电站成为下游电站的来水。这意味着你的调度决策不能只看单个电站上游的一举一动都会影响下游未来几个小时的可用水量这个影响还有时间延迟。水库本身就是天然的储能装置。有调节能力的水库在光伏出力高的时段可以把水存住、压低水电出力给光伏让出消纳空间到了晚上光伏退出再加大放水发电。这个蓄放过程放在梯级系统里还需要考虑每个水库的库容上下的约束以及出库流量对下游入库流量的传递。上游如果为了短时补光伏缺口猛放水下游可能在几个小时后遭遇来水集中到顶被迫弃水等于把上游省下的电量又丢掉了。所以梯级水光互补的优化不是一个简单的功率分配问题而是一个带时序耦合的水量调度问题。理解了这个背景再看论文标题里的短期优化调度就知道模型要做的核心事情是什么在给定未来24小时的径流预报和光伏预测场景下求出每条河流每个时刻的蓄水量、发电流量、光伏消纳量使得期望可消纳电量最大。下面就把这个模型的数学骨架一层层拆开。2. 模型框架拆解目标函数、约束条件和不确定性建模2.1 目标函数把期望落到实处论文的目标函数可以写成这样max Σ_{s∈S} π_s · Σ_{t1..T} [ Σ_{i∈R} P_hydro(i,s,t) P_pv_used(s,t) ] · Δt其中s是光伏场景编号π_s是该场景对应的概率T是调度总时段通常24小时R是梯级水库集合。P_hydro是水电出力决策变量P_pv_used是光伏实际消纳量决策变量。这里有一个非常关键的理解点P_pv_used不是光伏预测出力而是决策变量。它的上界是场景s下时段t的光伏可用出力目标函数鼓励它尽可能大但约束会限制它不能超过电网消纳能力。换句话说模型允许弃光但弃多少、在哪里弃要由优化器来选择。很多论文复现者第一步就栽在这里——把P_pv_used写成等于预测值相当于把模型从优化光伏消纳变成了固定光伏消纳水电的调节作用直接被废掉一大半。后面我会单独用一节讲这个坑。目标函数里的Δt需要统一量纲功率用MW、时段按小时算乘起来就是MWh这样期望可消纳电量就是MWh级别的一个数值方便和论文结果对比。别小看这个单位统一好多复现代码表面上看能跑算出来的电量结果却大得离谱根源就是Δt混用了秒和小时。2.2 约束条件水量平衡是最容易错的一道关模型的基础约束是水量平衡方程每个水库每个场景每个时段都必须满足V(i,s,t1) V(i,s,t) Δt_water · ( I(i,s,t) Q_up(i,s,t-τ_i) - Q(i,s,t) - Q_spill(i,s,t) )每个量的含义我都标清楚V是水库蓄水量I是区间入流比如支流汇入和降雨Q_up是上游电站的出库流量经过时滞τ_i后到达本站Q是发电流量Q_spill是弃水流量。Δt_water如果以秒为单位流量用m³/s乘起来得到的是m³这样库容和流量的单位自洽。时滞τ_i是梯级模型的灵魂所在。上游出库到下游入库不是瞬时的水流在河道里流动需要时间不同河段可能是1小时到几小时不等。如果忽略这个延迟你会看到下游水位约束被毫无道理地打破因为上游今天放的水不可能立刻到下游你却当成已经到了。除水量平衡外还有几组常见约束库容上下限约束Vmin(i) ≤ V(i,s,t) ≤ Vmax(i)以及调度周期末库容约束短期调度通常固定或约束末水位否则模型会把水全部放光只追求电量不顾后续。发电流量约束Qmin(i) ≤ Q(i,s,t) ≤ Qmax(i)弃水流量非负。水电出力约束P_hydro(i,s,t) 与发电流量、水头相关。严格的水电出力是9.81·η·H·QH为发电水头Q为流量两者乘积形成非线性项。后面我会讲如何做线性化处理。光伏消纳约束0 ≤ P_pv_used(s,t) ≤ P_pv_avail(s,t)即消纳量不会超过该场景下能达到的光伏出力。电网消纳上限约束Σ_i P_hydro(i,s,t) P_pv_used(s,t) ≤ P_grid_limit(t)这条约束和期望目标共同定义了可消纳电量的物理边界。这里每条约束都要按场景展开因为每个场景下光伏不同相应的最优水量分配也会不同。场景之间唯一需要真正耦合的是水库的初末库容关系——起始库容是所有场景共用的末期库容也要所有场景共同满足这样调度计划才具备可执行性。2.3 不确定性场景从预测误差分布到场景集合有了目标函数和约束还缺场景从哪来。复现论文时通常分三步处理不确定性。第一步统计光伏预测误差。拿历史预测曲线和实际出力曲线做差得到误差序列然后拟合分布。不建议直接用标准正态分布因为正态分布两端无限延伸抽样会出现负的光伏出力明显不合理。更实用的做法是截断正态或者用beta分布拟合并限制在[0, P_installed]区间内。第二步生成带时序相关的抽样。光伏预测误差在相邻时段是有相关性的——如果10点实际比预测高11点大概率也偏高。独立逐时段抽样会得到锯齿状高频抖动的场景和真实偏差规律不符。最简单的处理是用历史误差序列求协方差矩阵做多元正态抽样生成完整的24小时误差曲线再加到预测曲线上得到每个场景的光伏出力曲线。第三步场景缩减。原始抽样可能抽500个甚至1000个场景但把这些所有场景全部放进模型变量和约束规模会大到求解器难以接受。通常用聚类方法把相似场景合并挑出50到100个代表性场景再按簇内样本数占比确定每个代表场景的概率。这一步后面代码部分会给出具体实现。3. Python代码实现求解器选型与建模避坑3.1 为什么我选Gurobi而不是纯开源方案先说结论复现这类随机优化调度模型我最推荐Gurobi。不是因为它最贵而是因为学术用户免费。用学校邮箱在官网申请学术license可以全功能使用没有变量数限制这对搞论文复现的人来说是零成本高收益。对比情况我列在下面求解器授权模式适用场景我的评价Gurobi闭源学术免费LP/MILP/QP大规模问题首选性能稳、API顺手CPLEX闭源学术免费与Gurobi同级也可以但建模层不太常用SCIP开源MILP小规模够用场景放大后偏慢HiGHS开源LP/QP、部分MILP纯LP不错MILP大规模吃力如果单位不允许安装商业软件可以考虑先用HiGHS跑小规模模型验证逻辑提交论文结果时再换Gurobi。但我要提醒一句场景数超过80、水库数超过3时HiGHS的求解时间会明显拉长而Gurobi通常还在秒级到分钟级。论文复现里你不是在写工业代码而是在和论文结果对照追求的是求解时间可接受、解的质量可复现Gurobi是最稳的选择。3.2 核心建模代码从变量、约束到目标函数下面是我整理出的代码骨架去掉了部分非必要细节重点保留和论文公式对应的结构。整体流程是定义参数 → 建变量 → 加约束 → 设目标 → 求解。import gurobipy as gp from gurobipy import GRB import numpy as np # ---------- 基础参数 ---------- N_RES 3 # 梯级水库个数 N_SCEN 50 # 光伏场景个数 T 24 # 调度时段小时 WATER_DT_SEC 3600.0 # 水量平衡用的时间单位秒 HOUR_DT 1.0 # 电量目标用的时间单位小时 tau [0, 1, 2] # 上游到本站的水流时滞小时 pi_s np.ones(N_SCEN) / N_SCEN # 场景概率均匀假设 # 水库参数示意值 V0 [120.0, 85.0, 60.0] # 初始库容百万m3 Vmin [60.0, 40.0, 30.0] Vmax [180.0, 110.0, 80.0] Qmax [250.0, 220.0, 200.0] # 最大发电流量 m3/s coef [0.65, 0.60, 0.55] # 简化出力系数 MW/(m3/s)隐含水头近似 # 时段数据示意实际从Excel/CSV读取 inflow np.random.rand(N_RES, N_SCEN, T) * 50 20 # 区间入流 m3/s pv_avail np.random.rand(N_SCEN, T) * 150 # 光伏可用出力 MW grid_limit np.full(T, 180.0) # 消纳上限 MW # ---------- 建立模型 ---------- model gp.Model(cascade_hydro_pv) # 决策变量 V model.addVars(N_RES, N_SCEN, T 1, lb0, nameV) # 库容 Q model.addVars(N_RES, N_SCEN, T, lb0, nameQ) # 发电流量 Qspill model.addVars(N_RES, N_SCEN, T, lb0, nameQspill) # 弃水流量 Ph model.addVars(N_RES, N_SCEN, T, lb0, namePh) # 水电出力 Ppv_use model.addVars(N_SCEN, T, lb0, namePpv_use) # 光伏消纳量 # ---------- 约束1水量平衡含时滞 ---------- for r in range(N_RES): for s in range(N_SCEN): for t in range(T): up 0.0 if r 0: idx t - tau[r] idx max(idx, 0) # t小于时滞时简化取0严格可用历史出库 up Q[r - 1, s, idx] * WATER_DT_SEC # 立方米 model.addConstr( V[r, s, t 1] - V[r, s, t] - ( inflow[r, s, t] * WATER_DT_SEC up - Q[r, s, t] * WATER_DT_SEC - Qspill[r, s, t] * WATER_DT_SEC ) 0, namefbalance_{r}_{s}_{t} )水量平衡这段是全文最需要细看的地方。我最初实现时把时滞索引写成了idx t - tau[r]没有加max(idx, 0)模型直接报索引越界。加了下限处理后又发现另一个问题t小于时滞时上游来水应该用调度开始前的历史出库值而不是0复现阶段我为了简化取0结果前几个时段的水位和论文对不上。后来改成读取历史实际出库记录作为边界才算完全对齐。继续加约束# ---------- 约束2库容上下限与末库容 ---------- for r in range(N_RES): for s in range(N_SCEN): for t in range(T 1): model.addConstr(V[r, s, t] Vmin[r], namefVmin_{r}_{s}_{t}) model.addConstr(V[r, s, t] Vmax[r], namefVmax_{r}_{s}_{t}) # 周期末库容回到设定值原文若用范围约束则改成不等式 model.addConstr(V[r, s, T] V0[r], namefVend_{r}_{s}) # ---------- 约束3发电流量上限 ---------- for r in range(N_RES): for s in range(N_SCEN): for t in range(T): model.addConstr(Q[r, s, t] Qmax[r], namefQmax_{r}_{s}_{t}) # ---------- 约束4水电出力简化线性化 ---------- for r in range(N_RES): for s in range(N_SCEN): for t in range(T): model.addConstr(Ph[r, s, t] coef[r] * Q[r, s, t], namefPh_{r}_{s}_{t}) # ---------- 约束5光伏消纳上界 ---------- for s in range(N_SCEN): for t in range(T): model.addConstr(Ppv_use[s, t] pv_avail[s, t], namefpv_{s}_{t}) # ---------- 约束6电网消纳上限 ---------- for s in range(N_SCEN): for t in range(T): model.addConstr( gp.quicksum(Ph[r, s, t] for r in range(N_RES)) Ppv_use[s, t] grid_limit[t], namefgrid_{s}_{t} ) # ---------- 目标函数期望可消纳电量最大化 ---------- obj gp.quicksum( pi_s[s] * HOUR_DT * ( gp.quicksum(Ph[r, s, t] for r in range(N_RES)) Ppv_use[s, t] ) for s in range(N_SCEN) for t in range(T) ) model.setObjective(obj, GRB.MAXIMIZE) model.optimize()跑完之后V.x、Q.x、Ppv_use.x就是各场景下的最优调度值。注意这里水量平衡里的库容单位是百万m³还是m³我代码里为了避免量纲混乱直接以m³为单位写实际工程数据如果是万m³需要先换算。代码写完后要逐行检查量纲我在调试时因为这个坑把库容约束放大了10000倍Gurobi报infeasible排查了半个小时才发现是V的单位和流量换算不一致。3.3 场景缩减用KMeans把500个场景压到50个场景生成完直接全部丢进模型大概率会卡。以500个场景、3个水库、24时段为例仅库容变量就有500×3×2537500个加上流量、出力、光伏变量总变量数轻松超过10万约束也是同等量级。虽然Gurobi硬吃也不是完全不行但求解时间会从几十秒变成几十分钟完全没必要。场景缩减是一个标准操作。from sklearn.cluster import KMeans # pv_scenes: shape (n_original_scenes, T)原始光伏场景曲线 kmeans KMeans(n_clusters50, random_state42).fit(pv_scenes) rep_scenes kmeans.cluster_centers_ # 50个代表场景 cluster_counts np.bincount(kmeans.labels_) weights cluster_counts / cluster_counts.sum() # 场景概率 # 用 rep_scenes 和 weights 替换原始场景其他约束、目标结构不变缩减后精度损失取决于原始场景的聚类程度。我的经验是预测误差分布平稳时500场景缩到50个期望目标值的变化通常在1%以内而求解时间可以从分钟级降到20秒级。如果论文里明确写了场景数那就按论文参数来但代码调试阶段建议先用50个场景跑通出结果后再放大验证。4. 完整的实现链路从数据准备到结果输出4.1 输入数据清单论文没给数据怎么破复现EI论文最难受的不是模型公式而是作者经常只给几张结果图具体流域参数、径流数据、预测误差参数一个都不给。我的应对方案是先把论文可用的参数抄下来比如水库个数、装机容量、预测误差标准差其余用公开典型数据补全。整理过一张数据清单照着准备不会漏数据项说明论文没给时的替代方案水库参数库容上下限、初始库容、最大发电流量、出力系数用该流域所在省份典型水电站公开参数近似径流与区间入流各水库各时段来水用历史水文站日平均流量按比例分摊到小时光伏装机与预测曲线24小时光伏可用出力用当地历史辐照度乘装机容量再乘转换效率预测误差分布均值、方差或协方差矩阵用装机容量的5%-15%作为标准差再调参对比电网消纳上限各时段允许接入的上限功率按负荷曲线或输电断面限额设定数据准备阶段有个原则先用看起来合理的数据把代码跑通再回头对照论文结果曲线微调参数。不要一上来就纠结某条径流曲线是否精确模型逻辑通了参数替换只是重跑一遍的事。4.2 迭代式开发先10个场景跑通再放大我强烈建议分三步走不要一步到位。第一步跑一个确定性模型即场景数设为1、光伏取预测值。这一步的目的是验证水量平衡、库容约束和出力约束没有写错。第二步跑10个场景验证随机逻辑——目标函数里场景概率是否生效Ppv_use是否合理。第三步再跑到50到100个场景出正式结果。每步都要加检查代码。我最常用的是把每个水库每个时段的库容打印出来肉眼扫一遍看是否在上下限之间再看看目标函数值和Gurobi求解日志里的约束数是否符合预期。如果模型报Infeasible优先检查三处末库容约束是不是与水量平衡矛盾、时滞索引是否越界、电网消纳上限是否小到所有电源同时开机都满足不了。4.3 结果曲线调度图和弃电率怎么画求解完成后先把结果写成CSV存档再画两张图。一张是各时段光伏消纳量、水电总出力、总出力与消纳上限的对比曲线另一张是各水库水位变化曲线。画图用matplotlib就够了import matplotlib.pyplot as plt pv_used_sol np.array([Ppv_use[s, t].x for s in range(N_SCEN) for t in range(T)]).reshape(N_SCEN, T) ph_sol np.array([Ph[r, s, t].x for r in range(N_RES) for s in range(N_SCEN) for t in range(T)]) ph_sol ph_sol.reshape(N_RES, N_SCEN, T) pv_avg pv_used_sol.mean(axis0) # 各时段光伏消纳期望 hydro_avg ph_sol.mean(axis1) # 各时段水电总出力期望 total_avg pv_avg hydro_avg plt.figure(figsize(10, 4)) plt.plot(range(T), pv_avg, labelPV used (expected)) plt.plot(range(T), hydro_avg, labelHydro (expected)) plt.plot(range(T), total_avg, labelTotal output (expected)) plt.plot(range(T), grid_limit, --, labelGrid limit) plt.xlabel(Hour) plt.ylabel(MW) plt.legend() plt.tight_layout() plt.show()指标计算就不能只看图了要按定义算期望可消纳电量是sum(pi[s] * sum(total[s,t] for t)) * HOUR_DT弃光率是1 - sum(Ppv_use) / sum(pv_avail)水能利用率可以定义为实际发电用水量与可发电用水量之比。算完这些指标再对照论文里给出的数量级比如弃光率是百分之几而不是百分之几十总消纳电量在几百MWh到几千MWh量级基本就能确认复现是否到位。5. 复现中踩过最深的三个坑5.1 水头-出力关系非线性模型直接做成了非凸论文里的水电出力公式常写成 P 9.81·η·H·QMWH是水头mQ是发电流量m³/s。问题出在H本身是库容的函数而库容又是变量于是出现变量×变量的乘积项整个约束变成非线性的。直接丢给Gurobi要么报unsupported要么把它当MIQP处理求解速度慢到怀疑人生。我在复现时采用了两段处理。第一阶段先用固定平均水头近似即把H当成常数P_hydro就变成了Q的线性函数整模型变成纯LP几十秒出解。这个近似在短期调度里是合理的因为一天内水库水位变化幅度通常不大水头变化对出力的影响属于二阶误差。第二阶段如果论文要求水头精确建模再对H-Q出力曲线做分段线性逼近把每个线性段用SOS2约束实现模型变成MILP速度会慢一些但仍是可解的。5.2 可消纳是决策变量不是预测值直接拿来用这个坑我在前面提过但值得再展开一次。最初实现时我把Ppv_use写死为预测值理由是光伏发多少我就消纳多少除非超过上限。跑出来的结果很离谱午间光伏预测80MW、电网上限100MW、水电可以发30MW的时候总出力变成110MW模型被迫弃掉10MW但水也没有存下来。为什么因为光伏出力被写死后水电只能被动跟着让路根本没有商量余地。改成决策变量后模型自己算出来最优解是水电22MW、光伏78MW总分刚好100MW电网容量被充分利用水库还多存了8MW水量的水晚上光伏归零时再放出来。这一改期望可消纳电量直接提升了几个百分点。判断你的模型写对了没有可以做一个简单测试把光伏可用出力设成零模型应该自动提高水电出力来补缺口再把电网上限设成很小模型应该自动降低光伏消纳而不是硬顶上限。如果这两个行为都正常说明决策变量的设计是合理的。5.3 论文公式和代码之间的落差用小案例反推验证EI论文的符号经常不统一同一篇里Q可能一会儿表示流量、一会儿表示电量时滞的定义也不写清楚。遇到这种情况我不建议直接对着大系统调试而是构造一个极简案例1个水库、1个时段、无光伏。这个案例你可以手算出最优解比如初始库容100、最大发电流量限制对应出力50MW那最优解就是发50MW末库容等于初始库容减去发电流量对应水量。把代码参数改成这个极简案例跑出来如果和手算一致说明水量平衡和目标是对的再逐步加入光伏、第二级水库、时滞每加一块都用行为是否符合物理直觉来校验。这样一层层叠上去基本能做到模型的每部分都有确凿依据。如果只是把大模型跑完看总结果一旦数字对不上你完全不知道是哪个约束错了排查成本反而更高。6. 验证与延伸随机模型到底带来多少提升6.1 确定性模型和随机模型的结果对比我复现时用同一套流域参数跑了两个模型一个只考虑光伏预测均值确定性一个按50个场景做随机优化下面是示意结果不是原论文数字只是展示对比方向指标确定性模型随机优化模型变化期望可消纳电量MWh428043150.8%弃光率8.7%6.1%-2.6个百分点水能利用率91.2%94.5%3.3个百分点提升看起来不大但关键不是数值而是来源。确定性模型在预测偏高时会浪费水电调节空间在预测偏低时又会占用光伏消纳通道随机模型看到全部场景后会在水库蓄放上做一个更稳妥的取舍——到底留多少库容应对光伏晚高峰到底放弃多少光伏给水电腾空间这些决策在全场景期望意义下是最优的。也正是这种牺牲个别场景、提升平均表现的特征让期望模型在实际运行中比单点预测更可靠。6.2 这个模型还能往哪里扩展代码结构搭好以后扩展方向其实很清晰。加储能只需要新增充放电变量和SOC荷电状态转移约束目标函数里加储能放电出力做现货市场优化把电量目标换成收益目标加入价格场景即可考虑日前-日内两阶段决策则要把模型改成两阶段随机规划第一阶段变量是水库蓄放计划第二阶段变量是光伏消纳和水电修正出力约束层级要重新组织。还有一个方向是分布鲁棒优化不假设光伏误差服从某个固定分布而是构造一个包含所有可能分布的模糊集目标改成最坏分布下的期望收益这在学术上更前沿代码也是在现有模型基础上加对偶约束。这次复现给我最大的体会是EI论文的模型骨架其实不难难的是把工程细节填对。目标函数就一行约束加起来也不超过十组但单位、时滞、决策变量自由度这些细节任何一个没想透求解器都会用一条infeasible或一个离谱数字来提醒你。最后分享一个实用技巧写这类随机优化模型先用10个场景跑通再逐步加到100个。如果10个场景都跑不动基本可以断定是约束规模写错了而不是求解器的问题用这个习惯调试能省下大把和求解日志死磕的时间。