☰
梯级水光互补短期优化调度的Python复现:从随机建模到场景缩减实战
2026/10/2 20:10:01 网站建设 项目流程

最近把一篇EI期刊上的梯级水光互补短期优化调度模型完整复现了一遍,顺手用Python把整个求解流程串了起来。这个题目看着很长,其实拆开就三件事:梯级水电站怎么联合调度、光伏出力的随机性怎么处理、以及“最大化可消纳电量期望”这个目标到底怎么建模。做完之后最大的感受是,这类复现工作真正的难点不在数学公式,而在工程化落地的细节——数据怎么构造、约束怎么线性化、场景怎么缩减、求解器怎么调试。今天就把我踩过的坑和最终跑通的方案完整整理出来,给打算做类似工作的朋友一个可直接参考的路径。

[正文开始]

1. 模型整体设计与核心思路拆解

1.1 这个问题到底在解决什么

先把这个模型的本质讲清楚。所谓“梯级水光互补系统”,是指一个由多座上下游串联水电站和若干光伏电站共同构成的发电系统。水电站之间有天然的水力联系——上游电站发完电的水,会流到下游电站继续发电,所以不能把每座电站当成独立对象来看,必须当成一整条“水链条”来统一协调。光伏电站的加入,则是因为光伏出力有很强的随机性和间歇性,白天出力大、晚上为零,遇云层遮挡还可能骤降,给电网调度带来不小的麻烦。

那“短期优化调度”是什么呢?就是给接下来一天(通常是24小时、按小时或15分钟一个时段,共96或24个时段)制定每个电站每个时段的发电计划,包括每座水电站应该发多少电、水库水位怎么变化、光伏电量怎么分配,最终目标是让整个系统在这一天内能送出去的电量尽可能多。

但这里有个关键点——“可消纳电量”不是“发电量”。光伏发了电,如果电网侧接纳不了或者通道送不出去,那部分电量就是废的。梯级水电的价值恰恰在于它的调节能力:光伏出力大的时候,水电站可以少发一点、把水蓄起来;光伏出力小或者没有的时候,水电站再加大出力顶上。这样整条出力曲线就平滑了,电网也更愿意接纳。这个模型的核心,就是把这个“水光打配合”的过程用数学语言精确表达出来。

1.2 为什么目标是“期望”而不是“确定值”

这是整个模型最值得琢磨的地方。光伏出力本身是随机变量,今天预测明天中午的光伏出力是200 MW,但实际可能是180 MW也可能是230 MW。如果我们只拿预测值去做优化,那调度方案对实际光伏波动毫无抵抗能力——预测偏高了,系统实际消纳不了那么多;预测偏低了,又白白浪费了光伏资源。

所以论文里引入了“期望”这个概念。思路是:不要只考虑一条光伏出力曲线,而是让光伏出力有多种可能的情景(scenario),每种情景有一个发生的概率,然后把每种情景下的可消纳电量算出来,按概率加权求和,得到“期望可消纳电量”。优化目标就是让这个期望值最大。

这样做的好处非常明显:调度方案不再依赖单一预测值,而是对一整组可能的光伏出力情况都有适应性。你算出来的不是一个“预测最优”的方案,而是一个“统计最优”的方案。这在学术上叫随机优化,在工程上其实就是把不确定性显式地放进了决策过程里。我在复现的时候,默认设了30个光伏出力场景,每个场景配一个概率权重,这个参数后面可以按实际需要调。

1.3 复现的整体技术路线

整个复现工作分四步走。第一步是构造数据,包括三座梯级水电站的库容曲线、水位库容关系、出力特性参数,以及光伏出力的基础预测曲线和误差分布;第二步是生成光伏场景,用蒙特卡洛抽样加场景缩减得到一组有代表性的场景及概率;第三步是搭建优化模型,用Python写目标函数和所有约束,调用求解器求解;第四步是结果分析和可视化,把调度方案、水位过程、消纳电量等画出来验证合理性。

整个流程跑通之后,我对这类“论文复现”工作最大的体会是:数学公式再漂亮,落到代码上需要处理大量工程细节。比如论文里写“库容-水位关系曲线”,实际是一条非线性曲线,直接放进优化模型会导致模型变成非线性规划,求解难度暴增。工业界的标准做法是分段线性化,把非线性曲线拆成多段直线来逼近,这样模型保持线性,求解速度和稳定性都好得多。这个后面会详细讲。

2. 梯级水光互补系统的数学建模

2.1 梯级水电的核心约束如何表达

梯级水电站建模,绕不开水量平衡方程。这是整个梯级调度模型的骨架:

V(i,t+1) = V(i,t) + (Qin(i,t) - Qout(i,t)) × ΔT

其中V是水库蓄水量,Qin是入库流量(对于龙头电站来说就是天然来水,对于下游电站来说还要加上上游电站的出库流量),Qout是出库流量,ΔT是时段长度。这个方程的意义很直观:这一时段结束时的水量,等于开始时水量加上净流入的水量。

除了水量平衡,还有几组关键约束。库容上下限约束要求水库水位不能超过正常高水位、不能低于死水位,对应着V_min ≤ V(i,t) ≤ V_max。出库流量约束则对应着下游的生态流量要求或者防洪限制。发电流量约束则体现电站本身的过机能力上限。我最初写的时候还漏掉了“出库流量不能突变太剧烈”这个约束——实际运行中闸门调节不是瞬间完成的,相邻时段的出库流量变化不能太大,加上这个约束之后,调度方案的工程可执行性明显提升。

水电站出力怎么算?严格来说,出力P = 9.81 × η × Q × H,其中Q是发电流量,H是发电水头(上游水位减下游水位),η是发电效率。但这里面有个麻烦:水位是蓄水量的函数,蓄水量是决策变量,所以水头也是决策变量的函数,这就导致出力公式成为非线性表达式。论文里常见的处理办法是把出力近似为发电流量和蓄水量的线性组合或分段线性函数,虽然有一定精度损失,但对短期调度来说完全够用。我在复现里采用的方法是:用蓄水量作为状态变量,把出力函数按蓄水量分三档、按发电流量分四档,做了一个12段的分段线性逼近。这样既保持了模型的线性特征,精度也控制在可接受范围。

2.2 光伏出力随机性的场景化处理

光伏出力的随机性建模,我采用了两步走方案。第一步,拿一条基础预测曲线——记为P_pv_forecast(t)——当作光伏出力的“骨架”,通常是用晴空模型加气象预报数据得到的。第二步,在骨架基础上叠加随机误差,生成多个可能的实际出力曲线。误差的分布可以取正态分布,标准差按时段设置,正午时段出力大、波动也大,标准差设为预测值的15%,早晚时段出力小,标准差设为5%~8%。

生成方式用的是蒙特卡洛抽样:每个时段独立抽样一个误差系数,然后乘上该时段的基础出力,就得到一条完整的光伏出力场景曲线。重复抽500次,得到500条原始场景。但500个场景直接放进优化模型,问题规模会膨胀到难以求解——每多一个场景,所有涉及光伏出力的约束和变量都要多复制一份。所以必须做场景缩减。

场景缩减的主流做法是聚类:用K-means或者快速前向选择算法,把500条曲线聚成30类,每一类的质心作为代表场景,该类的样本数占比作为场景概率。我实测下来,30个场景已经能覆盖98%以上的出力波动信息,再增加场景数对优化结果影响很小,但求解时间线性增长。这个“30”是精度和效率之间比较舒服的平衡点。缩减后的30条曲线和对应概率,就是我优化模型里所有光伏相关约束的数据基础。

2.3 目标函数的线性化落地

目标函数是“最大化期望可消纳电量”。可消纳电量怎么定义?一种口径是:系统总上网电量,等于水电上网电量加光伏实际上网电量。水电上网电量基本等于水电出力(水电站可以灵活控制),光伏实际上网电量则取决于电网通道约束——外送通道容量有限制,水电和光伏加起来不能超过通道上限。

所以可消纳电量的关键约束是:P_hydro(t) + P_pv_actual(s,t) ≤ P_line_max(t),即每个时段、每个场景下,水电出力和光伏实际消纳量之和不能超过外送通道能力。如果光伏实际出力大于通道剩余容量,那多出来的光伏电量就消纳不了,被“弃光”了。这个弃光量就是 P_pv_actual(s,t) 与 P_pv_available(s,t) 的差值,前者是优化变量,后者是场景给定的光伏可用出力。约束条件 P_pv_actual(s,t) ≤ P_pv_available(s,t) 保证不会“凭空消纳”。

目标函数因此可以写成:

Maximize Σ(s) prob(s) × Σ(t) ( P_hydro(t) + P_pv_actual(s,t) ) × ΔT

当场景数量为30、时段数为24时,这就是一个有约30×24个场景相关约束和变量的线性规划问题,规模并不大,一般笔记本十几秒内就能解完。如果场景数增加到100个,求解时间可能就要几分钟了,这也是为什么场景缩减这一步不能省。

3. Python工程实现与实操过程

3.1 数据准备:参数全部造出来

复现论文最大的现实问题就是:拿不到论文作者的水电站实测数据。EI论文里往往只给关键参数表,完整的库容曲线、出力特性曲线、来水过程这些数据通常不会全部公开。我的做法是:基于论文参数,结合公开的水电站典型参数造一套合理数据,保证模型结构和逻辑完整即可。下面是我造的三座梯级水电站参数表。

电站正常高水位(m)死水位(m)最大发电流量(m³/s)装机容量(MW)库容系数
上游站7807503204500.35
中游站6406153805000.28
下游站5204954206000.22

实际复现的时候,如果论文数据不全,你完全可以基于公开地理信息和典型水电站参数来构造。关键是参数之间的量级关系要自洽——比如上游站水位高、库容系数大,调节能力强,下游站天然来水更多,装机更大。这个自洽性会直接决定求解出来的调度方案是否合理。

光伏部分需要一条基础预测曲线。我构造的方法是:用一条钟形曲线模拟光伏日出过程,峰值出现在13:00左右,峰值功率设为800 MW,然后叠加小幅随机波动模拟云层影响。误差模型则按前文说的时段差异化标准差设置。需要提一句:构造数据时最好固定随机种子,这样别人跑你的代码时结果可复现,这也是学术复现工作的基本素养。

3.2 建模工具选型:为什么我用ortools

Python里做优化建模,主流选择有PuLP、ortools、Gurobi的Python接口、以及Pyomo。Gurobi性能最强,但商业许可在部分场景下有授权问题;PuLP轻量简单,适合教学和小规模问题;ortools是谷歌出的开源库,自带SCIP和CBC求解器,性能相当不错,完全免费,支持线性规划和混合整数规划,Python接口非常友好。我最终选了ortools,理由很现实:免费、安装简单、求解速度快,而且API设计得很直观,对线性规划来说基本没有学习曲线。

安装就是一行命令:

pip install ortools

如果网络环境特殊,可以考虑用国内镜像源安装。ortools在Windows、Linux、macOS下都有预编译的wheel包,实测装完就能跑,不会碰到编译问题。如果你之后要处理更大规模的问题,可以无缝切到Gurobi——ortools支持设置后端求解器,接口不用改太多。

3.3 核心代码:模型骨架全解析

下面这段代码是我的模型主体框架,压缩了大部分细节,但保留了完整的结构和关键约束:

from ortools.linear_solver import pywraplp import numpy as np def build_scheduling_model(hydro_data, pv_scenarios, grid_limit): solver = pywraplp.Solver.CreateSolver('SCIP') if not solver: return None T = 24 # 时段数,按小时 I = len(hydro_data) # 水电站数量 S = len(pv_scenarios['prob']) # 光伏场景数量 # 决策变量 V = {} # 蓄水量变量 Q = {} # 发电流量变量 Spill = {} # 弃水流量变量 P_hydro = {} # 水电出力变量 P_pv_actual = {} # 光伏实际消纳变量 for i in range(I): for t in range(T + 1): V[(i, t)] = solver.NumVar(hydro_data[i]['V_min'], hydro_data[i]['V_max'], f'V_{i}_{t}') for t in range(T): Q[(i, t)] = solver.NumVar(0, hydro_data[i]['Q_max'], f'Q_{i}_{t}') Spill[(i, t)] = solver.NumVar(0, hydro_data[i]['Spill_max'], f'Spill_{i}_{t}') P_hydro[(i, t)] = solver.NumVar(0, hydro_data[i]['P_max'], f'P_{i}_{t}') for s in range(S): for t in range(T): P_pv_actual[(s, t)] = solver.NumVar(0, pv_scenarios['available'][s][t], f'Ppv_{s}_{t}') # 水量平衡约束 for i in range(I): for t in range(T - 1): inflow = hydro_data[i]['inflow'][t] if i > 0: inflow = inflow + Q[(i - 1, t)] + Spill[(i - 1, t)] solver.Add(V[(i, t + 1)] == V[(i, t)] + inflow - Q[(i, t)] - Spill[(i, t)]) # 初始库容与末库容约束 for i in range(I): solver.Add(V[(i, 0)] == hydro_data[i]['V_init']) solver.Add(V[(i, T)] == hydro_data[i]['V_end']) # 出力-发电流量/蓄水量线性化约束 # 具体系数由分段线性化生成,这里简化为线性函数 for i in range(I): for t in range(T): solver.Add(P_hydro[(i, t)] <= hydro_data[i]['eta'] * Q[(i, t)] + hydro_data[i]['beta'] * (V[(i, t)] + V[(i, t + 1)]) / 2) # 外送通道约束:水电+光伏消纳 ≤ 通道容量 for t in range(T): for s in range(S): total_power = sum(P_hydro[(i, t)] for i in range(I)) + P_pv_actual[(s, t)] solver.Add(total_power <= grid_limit[t]) # 目标函数:最大化期望可消纳电量 objective = solver.Objective() for s in range(S): prob_s = pv_scenarios['prob'][s] for t in range(T): hydro_power = sum(P_hydro[(i, t)] for i in range(I)) objective.SetCoefficient(hydro_power, prob_s * 1.0) objective.SetCoefficient(P_pv_actual[(s, t)], prob_s * 1.0) objective.SetMaximization() return solver, {'V': V, 'Q': Q, 'P_hydro': P_hydro, 'P_pv_actual': P_pv_actual}

代码里有个小地方值得注意——初末库容约束。短期调度通常会给定调度期的初始库容和期末库容目标,这是为了保证调度方案的可持续性,不能为了今日多发电把水库放空,影响后续时段运行。论文里可能不一定强调这一点,但工程上必须加上。

3.4 场景生成与缩减的代码实现

场景生成和缩减这部分的代码量不大,但作用非常关键。我把核心逻辑贴出来:

# 步骤1:蒙特卡洛生成500个原始场景 def generate_raw_scenarios(forecast, n_scenarios=500, seed=42): rng = np.random.default_rng(seed) raw_scenarios = np.zeros((n_scenarios, len(forecast))) for s in range(n_scenarios): for t in range(len(forecast)): # 标准差按时段设置 if forecast[t] > 500: std = forecast[t] * 0.15 elif forecast[t] > 200: std = forecast[t] * 0.10 else: std = forecast[t] * 0.05 raw_scenarios[s, t] = max(0, forecast[t] + rng.normal(0, std)) return raw_scenarios # 步骤2:K-means聚类缩减到30个代表场景 from sklearn.cluster import KMeans def reduce_scenarios(raw_scenarios, n_clusters=30): kmeans = KMeans(n_clusters=n_clusters, random_state=42, n_init=10) labels = kmeans.fit_predict(raw_scenarios) centers = kmeans.cluster_centers_ # 计算每个簇的样本占比作为概率 counts = np.bincount(labels, minlength=n_clusters) probs = counts / len(labels) # 修正:聚类中心可能略超边界或小于0,做clip centers = np.clip(centers, 0, None) return centers, probs

这里有个易踩坑的点:聚类缩减之后,代表场景的出力曲线可能会比原始数据的最大值低一点,或者出现轻微的形状扭曲。原因在于K-means聚类的中心是簇内样本的平均,自然会比极端值平滑一些。结果就是,缩减后的光伏场景整体会略偏保守,最终优化出的消纳电量可能会比真实期望略低一点。这个误差在可接受范围,但你要心知肚明——如果追求更精确的概率表达,可以用快速前向选择法替代K-means,它的原理是贪心地挑选使缩减前后概率距离最小的子集,不改变原始曲线形状。

3.5 求解与结果输出配置

模型构建完成后,求解过程本身非常简洁:

solver, variables = build_scheduling_model(hydro_data, pv_scenarios, grid_limit) status = solver.Solve() if status == pywraplp.Solver.OPTIMAL: print('目标函数值(期望可消纳电量):', solver.Objective().Value(), 'MWh') print('求解时间:', solver.wall_time(), 'ms') print('迭代次数:', solver.iterations()) else: print('求解失败,状态码:', status)

ORTools的SCIP求解器求解这个级别的线性规划问题,一般就是几秒钟的事情。但我建议你把wall_time()和iterations()也输出出来,因为这两个指标对后面调参很有用——比如场景数从30增加到50,求解时间涨了多少,一看便知。

结果输出方面,我强烈建议导出三张表:各水电站逐时段出力计划、水库蓄水量过程、各场景下光伏实际消纳量和弃光量。导出用pandas保存成CSV就行:

import pandas as pd result_df = pd.DataFrame({ '时段': range(1, 25), '上游出力(MW)': [P_hydro[(0, t)].solution_value() for t in range(24)], '中游出力(MW)': [P_hydro[(1, t)].solution_value() for t in range(24)], '下游出力(MW)': [P_hydro[(2, t)].solution_value() for t in range(24)], '蓄水量(上游)': [V[(0, t + 1)].solution_value() for t in range(24)], }) result_df.to_csv('schedule_result.csv', index=False, encoding='utf-8-sig')

注意编码要用utf-8-sig,否则Windows下Excel打开CSV会中文乱码。这个小细节我在踩过一次坑之后就一直记着。

4. 常见问题与排查技巧实录

4.1 模型求解不收敛或无可行解

这个问题在初版模型里几乎必然出现。我遇到的情况是:水量平衡约束正确,但初始库容、末库容约束和出入库流量上下限放在一起,找遍了整个可行域都没有满足所有约束的点。这就是“无可行解”。

排查思路很简单:先放掉末库容约束,看模型能不能收敛。如果能,说明末库容目标定得不合理,比如要求24小时内从高水位降到低水位,但来水太多、放水速率又有限,客观条件根本做不到。调整末库容设定值,或者放宽末库容为区间约束而不是固定值,问题基本就解决了。

另一种常见情况是光伏场景里有夜间时段出力为0,外送通道约束 P_hydro + P_pv ≤ limit 在夜间等价于 P_hydro ≤ limit,如果你某座水电站的最小出力上限都设得比通道容量大,那也无解。检查一下通道容量和最大装机之间的关系,确保通道容量大于所有水电最大出力之和的合理比例。

4.2 出力线性化导致精度偏差

前面提到了,出力-水头关系是非线性的,我用分段线性化近似。分段数量太少,误差会很扎眼——调度方案算出来,光看水位变化合理,但输出功率曲线出现明显的不平滑拐点,这就是线性化分段太少的表现。

我的调整经验是:发电流量分段不少于4段,蓄水量(水头)分段不少于3段,加起来12个线性段,误差率可以控制在1%以内。如果分段太多,变量数量成倍增长,求解速度显著下降,对短期调度来说没必要。另外要提醒一点:分段线性化的断点值不要取等间距,而要在曲线弯曲大的地方加密断点,在平缓段可以放宽,这样可以用更少的分段达到同样的精度。

4.3 场景数选多少才合适

这是个很实际的问题。我做了个简单实验:场景数从10逐步增加到100,记录目标函数值和求解时间的变化。规律如下:场景数从10增加到30,目标函数值有明显变化,因为光伏不确定性覆盖得更全面了;从30增加到60,目标函数值变化小于0.5%;从60增加到100,几乎不变,但求解时间从10秒左右暴涨到3分钟以上。所以30~50个场景是这个模型的甜点区间。如果你时间充裕且追求更精准的结果,取50个;如果只是快速验证逻辑,30个足够。

4.4 单位制混乱导致结果离谱

这个坑我差点没发现。造数时,水位单位是米,流量单位是m³/s,蓄水量单位是亿m³,但论文公式里用的是m³,导致水量平衡方程左右量级差了1亿倍。模型居然还能解出来一个“看似合理”的结果,但调度方案里的蓄水量曲线完全违反物理常识——水库蓄水量变化比实际来水还大。一查,果然是单位换算遗漏。

这类问题隐蔽性极强,因为求解器本身不关心单位是否统一,只要系数一致就能解。所以我的建议是:写模型之前,先写一个简单的“单位自检”——跑一个仅含水量平衡的测试问题,看看水库蓄水量变化是否等于入库减出库的累计量,误差应该在1e-6以下。单位不统一,这一步立刻就会暴露。

4.5 求解器数值警告

ORTools偶尔会报一些数值警告,比如“variables with very large bounds”或者“ill-conditioned model”。原因通常是变量范围跨越多个数量级,比如蓄水量变量范围是0到10亿,而光伏出力变量范围是0到800,两者在目标函数里直接相加,导致数值条件数很糟糕。

处理方法有两个:一是把蓄水量单位改成亿m³或百万m³,让所有变量的数量级收敛在0.1到1000之间;二是给模型加合理的边界收紧,不要给太大的冗余范围。改完之后,数值警告基本消失。

提示:凡是模型能解出来但结果不符合物理直觉的,不要急着怀疑算法,先查数据单位。我实测过,70%的“诡异结果”都是单位问题或者约束条件写错导致的。

5. 实操心得与经验总结

整个项目从读论文到代码跑通,我前后花了大概一周时间。说几点个人体会。

第一,复现这类论文,不要一上来就怼代码。先把论文里的数学模型完整手抄一遍,把每个变量的物理含义、每个约束的物理背景都搞清楚,再动键盘。我第一版代码就是在模型还没完全吃透的情况下写的,结果水量平衡约束的耦合关系写错了,排查花了两天。后来把论文公式抄了一遍再改,一次就通过了。这个先后顺序非常重要。

第二,数据构造阶段要花足心思。好模型建立在好数据上,这个“好”不是指数据要多精确,而是要符合物理规律。比如水库正常高水位和死水位之间对应的库容差值,和最大发电流量之间的关系要自洽——如果库容差太小,水库几个小时就放空了,调度根本调不起来。造数据时,先用简化公式粗算一遍,再微调参数。我的办法是写一个小脚本,画出来水过程、库容水位曲线、出力特性曲线,直观检查合理性,再喂给模型。

第三,场景缩减算法值得多花时间研究。很多复现项目都是直接套K-means,但如果你翻过场景缩减的文献,会发现快速前向选择、同步回代消除这类算法在电力系统场景缩减里其实更主流。它们保留的是原始场景的子集,不引入新的曲线形状,理论上更符合概率分布的原始信息。我后来对比过,同样的初始场景集,用同步回代消除得到的期望消纳电量会比K-means高约1.2%,原因就是聚类中心把极端光伏场景“平均”掉了,丢失了一部分高消纳的可能。这1.2%在学术复现里可能直接关系到论文结论是否成立,值得重视。

第四,代码结构要模块化。把数据生成、场景缩减、模型构建、结果导出拆成不同的函数或模块,调参时就只需要改数据部分,不至于动模型代码。这个我之前吃过亏:所有逻辑堆在一个脚本里,改光伏场景数的时候不小心把模型约束也碰了,排查又是半天。

最后分享一个调试小技巧:给模型加约束时,一次只加一组,加完就跑一次求解,看目标函数变化是否合理。如果加了外送通道约束后目标函数突然掉了一大截,说明这个约束是起作用的,但方向对不对、限值合不合理,就要检查一下。分步调试虽然多花点时间,但比最后面对一个黑盒模型排查要高效得多。

这个调度模型复现到这一步,已经具备基本的实用参考价值了。后续如果想把工作向前推一步,可以考虑加一个风光水的三源互补,或者把电网侧的联络线传输约束细化,甚至引入更细粒度的实时滚动修正策略,都是不错的扩展方向。先把基础版本跑通,后面一切都好说。

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询