看到这个标题我第一反应是“又是水电、又是光伏、又是期望、又是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,t1} 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 ApproximationSAA。假设光伏预测曲线是 (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 matplotlibCBC可以用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_col0) inflow pd.read_csv(inflow.csv, index_col0) 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, withinNonNegativeReals) # 库容 model.Qh Var(model.H, model.T, model.S, withinNonNegativeReals) # 发电流量 model.Sp Var(model.H, model.T, model.S, withinNonNegativeReals) # 弃水流量 model.Ph Var(model.H, model.T, model.S, withinNonNegativeReals) # 水电出力 model.Ppv_use Var(model.T, model.S, withinNonNegativeReals) # 光伏实际消纳 model.Cpv Var(model.T, model.S, withinNonNegativeReals) # 弃光 # 目标函数期望可消纳电量最大化 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(ruleobj_rule, sensemaximize) # 水量平衡约束 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, t1, 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, rulewater_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, rulevlimit_low_rule) model.v_up Constraint(model.H, model.T, model.S, rulevlimit_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, rulevinit_rule) model.v_end Constraint(model.H, model.S, rulev_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, rulepower_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, rulecurtail_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, rulechannel_rule)需要特别提醒的是上面这种把所有决策变量都按场景展开的写法严格来说是“场景独立决策”每个方案只在对应场景下最优没有考虑决策的不可预期性。复现初版可以这么跑但写论文对比时通常还需要加“非预期约束”比如某些时段的库容决策不随场景变化或者采用两阶段随机规划的写法。我建议先按全场景展开跑通再加非预期约束这样更容易定位问题。3.4 求解与结果导出模型构建完成后求解代码非常简单solver SolverFactory(cbc) results solver.solve(model, teeFalse) 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电站2最小库容 / 万m³10080最大库容 / 万m³500400初始库容 / 万m³300250目标末库容 / 万m³300250综合出力系数0.80.9天然入库流量取常数电站1每小时15万m³电站2每小时10万m³。惩罚系数设为弃光10元/MWh、弃水20元/MWh这个惩罚要高于正常情况下单位发电收益但又不能高到让模型为了“惩罚”而故意压低正常出力。实际论文里会用上网电价作为基准这里先随意给一个示意值。4.2 功率曲线与水位过程怎么读算例跑出来后最直观的是看各时段总出力曲线。我通常会画一张包含三条曲线的图光伏可发功率、光伏实际消纳功率、水电总出力。光伏大发的中午时段外送通道上限1000MW会被逼近水电会自动压低出力甚至可能出现停机蓄水。到了傍晚光伏快速下降的时段水电出力会明显抬升补足电力缺口。这就是“互补”在曲线上的直接体现。库容曲线则需要分开看两座电站。上游电站因为来水条件更好通常承担主要的调蓄功能库容变化幅度会更大下游电站除了自己发电还需要消化上游的出库流量库容变化往往受上游约束影响。如果发现某座电站的水位曲线一直贴着上限或下限运行就要注意是不是约束设置太紧或者惩罚系数不合理。4.3 期望消纳电量与弃光率对比场景数对目标值的影响非常值得观察。我做了一组对比实验同样是24时段两库系统CBC求解器结果大致如下示意值场景数期望可消纳电量 / MWh弃光率 / %求解时间 / s1确定性1685012.60.410173208.42.850177805.218.6100178904.887.3可以看出场景数从1增加到50时期望可消纳电量明显上升、弃光率明显下降。这是因为确定性模型只针对一条预测曲线优化遇到极端偏差场景时毫无准备而随机优化方案虽然在某些“平均场景”下不是最优但对各种可能出现的光伏场景都预留了水电调节空间整体期望值自然更高。但场景数继续从50增加到100时增益明显变缓。这是随机规划的典型现象样本均值逼近的误差随着样本数增加而收敛但增长速度递减。实际工程中不一定非要跑到200个场景要根据求解时间和精度需求权衡。5. 复现路上的坑与排查经验5.1 模型无可行解按什么顺序排查这个模型最常见的失败方式是求解器直接报“infeasible”第一次遇到时很容易让人懵。我的排查顺序固定三步第一步检查始末库容和天然来水是否自洽。如果初始库容等于目标末库容而全天然来水又不够发电流量的最低需求模型必崩。最简单的验证方法在目标函数里暂时去掉惩罚项把弃水变量和弃光变量固定为0只让水量平衡约束参与看是否存在解。第二步检查水量平衡约束中时间下标是否写错。比如V[h, t1]在t24时越界或者上游来水错用了t1时段的值。这类错误在Pyomo里有时候不会直接报错而是产生一个变量索引错误。第三步检查外送通道约束是否过紧。通道上限设置过低时水电被迫压出力但水量平衡又要求发电流量至少满足某一下限两个约束互相打架就完全不可行。我一般先把通道上限放松50%确认模型能解再逐步收紧看边界在哪里。5.2 求解时间爆炸场景削减和线性化并施随机场景数增加时MILP规模增长很快。100个场景、24时段、2个电站就已经是上万变量CBC可能跑几分钟甚至更久。这时候我常用的处理手段有三个一是场景削减不要直接拿200个原始场景糊上去。先用同步回代法把相似场景合并比如Gurobi或者scikit-learn的KMeans都能做把200个削到20个目标值损失往往不超过2%。二是检查有没有非必要的整数变量。如果水头用恒定线性近似整个模型完全没有整数变量只是一个大规模的线性规划求解速度会快非常多一旦分出控制水头的二进制变量模型立刻变成MILP难度上一个台阶。三是给求解器设定合理的MIP Gap。工程场景没必要非要0.00%最优gap0.5%时结果已经足够用于对比分析。Gurobi里设置solver.options[MIPGap] 0.005CBC里对应的是solver.options[ratio] 0.0055.3 数值与论文对不齐复现到底复现什么这是EI复现最折磨人的环节。好不容易把代码跑通结果和论文表里的数值对不上差个百分之五以内还好说差得多了就怀疑人生。我的建议是不要只盯随机种子。就算你把随机种子完全复现成论文一样的对比时也可能因为以下因素对不上论文的预测误差场景分布参数没有公开你用正态分布他可能用的Beta分布弃光弃水惩罚系数不同哪怕目标函数公式一样数值结果也会差很多水头处理方式不同恒定水头和分段线性水头算出来的发电流量边界完全两码事外送通道约束是硬约束还是罚函数对弃光率影响非常大。所以在复现阶段更实际的目标是复现“方法结构和规律”比如随机优化比确定性调度能提高多少消纳期望、场景数增加后目标值如何收敛、水库水位变化趋势是否符合物理直觉。这些规律和论文一致就说明你的模型逻辑是对的。5.4 惩罚系数怎么设我的一点经验惩罚系数是这类模型里最“玄学”的地方。设太小模型不在乎弃光外送通道一堵就让光伏弃着设太大模型会为了省惩罚而让水电机组满发反过来造成大量弃水同样不物理。我自己常用的基准是上网电价。假设上网电价为0.4元/kWh那么弃光惩罚可以取1.2到2倍上网电价弃水惩罚取2倍左右。这样模型的逻辑是如果光伏弃掉不如水电压低出力划算水电机组自然会压低如果水电为了避免弃水而多发电导致更低效惩罚会及时拉住它。调参技巧是固定其他条件只变一个惩罚系数看弃光率和弃水率的变化曲线直到两者处于一个比较均衡的水平。不要指望一次调好这个环节通常要跑很多轮。6. 复现完之后的几点个人体会做完整套复现我最深的感触是随机优化最难的其实不是求解而是把“期望”从一句黑话变成可计算的表达。当你在代码里第一次看到sum(pi_s * 目标量)跑完并输出一个收敛的期望值时数学模型和工程实践之间的距离一下子就缩短了。如果后续想继续扩展我建议可以在现有模型上做三个方向的小改一是把恒定水头换成分段线性水头让水电出力更贴近实际二是加入滚动时域调度每4小时用最新预测重新求解一次三是把目标函数从期望最大化改成“期望条件风险价值CVaR”的加权形式用来控制极端场景下的风险。每个方向的改动都不算大但做完之后你对调度模型的理解会上一个台阶。这个模型虽然只是简化复现但它已经足够让你掌握“梯级水电随机新能源”这类问题的完整建模方法论。