这个项目是我从一篇EI期刊论文的复现工作里整理出来的。题目叫“梯级水光互补系统最大化可消纳电量期望短期优化调度模型”，听起来很长，但拆开看其实很清晰：对象是梯级水电站加上光伏电站，手段是短期优化调度，目标不是简单追求发电量最大，而是让系统在电网消纳能力有限的前提下，期望可消纳电量达到最大。复现过程中我踩了不少坑，包括光伏场景怎么生成、水库水量平衡约束怎么线性化、用Python求解线性规划时矩阵怎么拼不会崩，今天把这套模型的建模思路和Python实现过程完整写出来，给正在做水电、新能源调度、或者论文复现的同学一个可以直接落地的参考。
一句话总结：这篇博文会讲清楚“最大化可纳电量期望”为什么不能简单理解成多发电，梯级水库和光伏怎么配合，以及用Python怎么把含有随机场景的短期调度模型写成可求解的线性规划。
1. 项目要解决的问题与模型构建思路
1.1 为什么盯上“可消纳电量期望”这个指标
光伏出力天然带波动，而且是跟着天气走的。一个光伏电站一天的出力曲线，可能上午爬坡、中午尖峰、下午剧降，遇到云层遮挡还会在短时间内大幅摆动。如果系统里只有光伏，电网能消纳多少基本取决于当时的负荷和通道能力，光伏发多了只能弃光。梯级水电站加入之后，水库可以像一个巨大的“充电宝”一样调节出力，光伏出力大的时候少放水、光伏出力小的时候多放水，让整个系统的外送电量尽量平稳地贴着电网消纳上限走。
但问题在于：光伏出力是随机变量，不是一组确定值。同一时刻可能有多种出力场景，不同场景下最优调度策略并不一样。如果只用一组预测值做确定性优化，结果往往偏乐观，实际运行中遇到预测偏差时水库调整不过来。于是模型目标就变成了“最大化可消纳电量期望”——也就是在所有可能的光伏场景下，系统能够被电网真正消纳的电量按概率加权求和，在这个期望值最大化意义下安排水库各时段的发电流量与库容。
这里关键要理解“可消纳”不等于“总发电量”。总发电量再大，如果超过外送通道上限或者电网需求上限，多余部分只能丢弃。所以在模型里通常会引入一个外送能力上限C_t，系统出力超过C_t的部分不算数。目标函数因此带了一个“截断”效果，数学上可以通过线性化处理，后面会详细讲。
1.2 梯级水光互补相比独立光伏、水库调度的优势
单独一个水库配合光伏，调节能力受限于单库库容和入流条件。梯级电站不一样，上游水库放水后，水流经过一段时间到达下游水库，下游水库可以在时间上错峰利用这部分水资源。这种上下游之间的水量衔接，给调度模型增加了一个时间维度上的灵活性。
举一个简单的例子：如果上游水库在光伏低谷时段加大发电流量，这部分水不会立刻到达下游，而是在几个小时后进入下游水库。下游水库可以选择在这段时间存水，等下一次光伏低谷再放水发电。于是系统中相当于有多级“储能缓冲”，光伏波动可以被更好地平抑。
但梯级也带来了约束上的麻烦。上游的出库流量直接影响下游的入库流量，所以水量平衡不能按单库孤立建模。同时每个电站都有库容上下限、发电流量上下限、出力上下限，末水位通常也有要求，所有这些约束耦合在一起，调度问题的变量数和约束数会成倍增长。这也是为什么需要把问题写成线性规划形式，而不是靠人工试算去运行。Python生态里有现成的线性规划求解接口，变量上千规模时也能在秒级到分钟级内解决，这对短期调度来说完全够用。
1.3 短期调度的时间尺度与数据需求
短期优化调度一般指日前到日内时间尺度，常见时段长度为1小时，调度周期为24小时，也有的论文用15分钟或30分钟时段做日内滚动。这个项目以1小时时段、24小时周期为例，光伏场景用预测误差抽样生成，水库径流采用给定的入库流量预报序列。
数据准备是复现这类项目最费时间的一步，需要五类基础数据：
- 梯级水库参数：各水库正常蓄水位、死水位、对应库容、装机容量、出力系数、发电流量上下限、初始库容和末库容。
- 入库径流数据：各水库逐时段的天然入库流量预报值，单位常用立方米每秒。
- 光伏出力预测数据：未来24小时每个时段的光伏出力预测值，以及预测误差分布参数。
- 外送通道限制：每个时段系统能够向电网输送的最大功率。
- 场景数据：基于光伏预测误差采样生成的多个光伏出力场景，以及每个场景的概率权重。
有了这些数据，才能把目标函数和约束写成可计算的代数表达式。
2. 数学模型设计与约束条件梳理
2.1 目标函数：最大化可消纳电量期望的表达式
先给出符号定义。设调度周期时段数为T，梯级电站数量为I，光伏场景数为S。i表示电站序号，t表示时段，s表示场景。决策变量包括：各电站各时段发电流量q_{i,t}、库容V_{i,t}、弃水流量spill_{i,t}，以及各场景各时段系统的实际可消纳功率R_{s,t}。光伏场景出力记为PV_{s,t}，场景概率为p_s。
目标函数为：
[ \max \sum_{s=1}^{S} p_s \sum_{t=1}^{T} R_{s,t} \cdot \Delta t ]
其中Δt为时段长度，R_{s,t}满足：
[ R_{s,t} \le C_t ]
[ R_{s,t} \le PV_{s,t} + \sum_{i=1}^{I} P_{i,t} ]
两个约束一起限定了R_{s,t}必须同时小于外送通道上限C_t和系统总出力。因为目标是对R求最大化，所以优化过程会自动让R取两个上限中的较小值，也就是实际可消纳电量。这里没有引入0-1变量，逻辑上成立，因为R是连续变量。
水电站出力P_{i,t}用简化线性公式计算：
[ P_{i,t} = K_i \cdot H_i \cdot q_{i,t} ]
K_i为电站出力系数，H_i为平均发电水头，这里做了恒定水头近似。更精确的做法是让水头随库容变化，但那样会引入非线性，需要通过分段线性化处理。实际复现中，如果文献没有特别强调水头变化，先用平均水头是完全可接受的，后文会提怎么修正。
2.2 关键约束：水量平衡、水位库容、出力限制
水量平衡约束是梯级调度模型的核心，它把上下游电站串成了一个整体。对于每个电站i、每个时段t，库容变化等于入库流量减去出库流量：
[ V_{i,t+1} = V_{i,t} + (Qin_{i,t} + Qup_{i,t} - q_{i,t} - spill_{i,t}) \Delta t ]
Qin_{i,t}为电站i的区间天然入库流量（立方米每秒），Qup_{i,t}为上游电站的出库流量。需要注意单位换算：如果库容单位是立方米，流量单位是立方米每秒，则流量的水量要用流量乘以时段秒数，一小时就是3600秒，Δt这一项不能漏。
出库流量由发电流量和弃水流量两部分组成。实际工程中弃水不是常态，但建模阶段必须给出弃水变量的允许范围，否则水库库容达到上限而入库流量又很大时，模型会出现不可行解。对于上游电站，Qup_{i,t}等于上一个电站的q_{i-1,t}+spill_{i-1,t}，这就是梯级水量耦合。
库容上下限约束：
[ V_{i}^{min} \le V_{i,t} \le V_{i}^{max} ]
发电流量上下限：
[ q_{i}^{min} \le q_{i,t} \le q_{i}^{max} ]
弃水流量非负：
[ spill_{i,t} \ge 0 ]
此外，首末库容要固定。初库容V_{i,1}等于调度开始时的实际库容，末库容V_{i,T+1}要保持在设定值或允许的小范围内。末水位约束非常关键，如果不加，优化结果会把水库放空来追求发电量，这在实际调度中不可接受。处理末水位约束时建议用等式，但为了数值稳定性，可以直接让V_{i,T+1}等于一个预设值。
还有一类约束是出力爬坡限制，论文里有时会写，有时不写。梯级水库调节能力较强，短期调度中主要限制是机组出力变化速率。如果模型中没有爬坡限制，优化结果可能出现相邻时段发电流量剧烈跳变的情况，实际机组跟不上。我的经验是，即使原文献没给爬坡参数，也要在复现时加一个宽松的爬坡限幅，避免结果曲线出现锯齿状。
2.3 光伏随机性与场景削减的工程处理
光伏出力随机性通过场景描述。最简单的方法是把预测值当作期望，加上服从正态分布或Beta分布的预测误差，生成大量原始采样场景。比如生成500个场景，每个场景包含24个时段的光伏出力值。场景数量越多，随机性描述越精细，但优化问题规模线性增长，求解时间变长。
场景削减是平衡精度和求解速度的常用手段。常用方法有同步回代削减（SBR）和K-means聚类。复现时我用的是K-means，因为它实现简单、行为稳定。把500个24维场景向量聚类成20~30类，每类均值作为代表场景，该类场景数量占总样本数量的比例作为概率。削减后的场景要尽量保留原始样本的统计特征，我实测下来，光伏出力期望和标准差在削减前后偏差能控制在2%以内，完全够用。
场景生成的关键参数是预测误差标准差。通常用光伏预测误差的分段标准差，或者简化成全天统一标准差。复现项目直接用全天统一标准差，取值范围在预测出力的10%~15%之间。误差分布如果选正态分布，会生成负的光伏出力值，需要截断到0；如果选Beta分布，则天然落在[0,1]区间，更符合光伏出力占比特性，但参数估计复杂一点。
3. Python代码实现与关键环节拆解
3.1 整体框架与依赖库选型
我用的环境是Python 3.9，核心依赖库如下：
- numpy：所有数组运算和矩阵拼接。
- pandas：清洗和管理水电、光伏数据。
- scipy：线性规划求解器scipy.optimize.linprog。
- matplotlib：结果曲线绘制。
- sklearn：K-means场景削减（sklearn.cluster.KMeans）。
为什么用linprog而不是gurobi或cplex？因为复现场景不需要商业求解器，scipy自带的HiGHS求解器对中规模线性规划效率很高。这个模型的变量数量在1500左右，约束数量在3000左右，HiGHS通常几秒内就能求解完成，完全满足短期调度需求。如果你想扩展成混合整数二次规划，再考虑换成puLP+CBC或mip库。
整个实现流程分为四步：数据初始化、场景生成、优化建模、结果输出。下面按流程展示关键代码。
3.2 数据初始化与场景生成代码
先定义水电站数据类。为了简洁，我用字典存参数，也可以用pandas DataFrame。
import numpy as np # 三个梯级电站的简化参数 stations = [ { "name": "上库", "V0": 120.0, "Vmin": 80.0, "Vmax": 150.0, "Vend": 120.0, "qmin": 0.0, "qmax": 120.0, "K": 9.8, "H": 30.0, # K为出力系数，H为平均水头 "Qin": [80, 85, 90, 95, 100, 105, 110, 115, 120, 125, 130, 135, 140, 135, 130, 125, 120, 115, 110, 105, 100, 95, 90, 85] }, { "name": "中库", "V0": 200.0, "Vmin": 150.0, "Vmax": 250.0, "Vend": 200.0, "qmin": 0.0, "qmax": 150.0, "K": 9.8, "H": 25.0, "Qin": [50, 55, 60, 65, 70, 75, 80, 85, 90, 95, 100, 105, 110, 105, 100, 95, 90, 85, 80, 75, 70, 65, 60, 55] }, { "name": "下库", "V0": 180.0, "Vmin": 130.0, "Vmax": 220.0, "Vend": 180.0, "qmin": 0.0, "qmax": 180.0, "K": 9.8, "H": 20.0, "Qin": [40, 42, 45, 48, 50, 52, 55, 58, 60, 62, 65, 68, 70, 68, 65, 62, 60, 58, 55, 52, 50, 48, 45, 42] } ] T = 24 I = len(stations) # 外送通道上限：每个时段系统最大可消纳功率，单位MW C_limit = 800.0 * np.ones(T) # 光伏预测出力曲线（单位MW） pv_forecast = 100 + 80 * np.sin(np.linspace(0, np.pi, T)) + 20 * np.random.randn(T) pv_forecast = np.maximum(pv_forecast, 0)场景生成函数，这里用K-means做削减。原始采样场景数量设500个，削减后保留20个。
from sklearn.cluster import KMeans def generate_pv_scenarios(n_raw=500, n_reduced=20, seed=42): rng = np.random.default_rng(seed) raw_scenarios = [] for _ in range(n_raw): error = rng.normal(0, 0.12 * pv_forecast + 5) # 标准差取预测值的12%加5MW scenario = np.maximum(pv_forecast + error, 0) raw_scenarios.append(scenario) raw_scenarios = np.array(raw_scenarios) kmeans = KMeans(n_clusters=n_reduced, random_state=seed, n_init=10) labels = kmeans.fit_predict(raw_scenarios) reduced_scenarios = [] probabilities = [] for k in range(n_reduced): cluster_data = raw_scenarios[labels == k] reduced_scenarios.append(cluster_data.mean(axis=0)) probabilities.append(len(cluster_data) / n_raw) return np.array(reduced_scenarios), np.array(probabilities)这段代码的要点有两个。第一，误差标准差不能用单一常数，因为光伏出力为零的夜里，误差也接近零；我用了预测值的比例项加一个基础项，更接近实际。第二，K-means的n_init至少设为10，避免聚类结果陷入局部最优。场景削减完之后，reduced_scenarios的每一行就是一个代表场景，概率向量代表该场景发生的可能性。
3.3 优化问题构建与求解核心代码
决策变量顺序我用的是：
- q_flat：I×T个发电流量变量
- v_flat：I×T个库容变量（对应每个时段初库容，再加最后一个末库容会有I×T个？为简化，变量只存V_{i,1}到V_{i,T}，用W公式处理末库容方程）
- spill_flat：I×T个弃水变量
- r_flat：S×T个可消纳功率变量
这里要特别注意变量边界：q有上下限，v有上下限，spill下限0，r没有显式上下限（靠约束限制）。如果给r也设置上限C_t，可能导致模型可解但目标被限制，但约束本身已经包含R≤C_t，所以不需要再设边界。给r设一个较大的上界也可以，但建议只依赖约束，方便检查问题。
代码如下：
from scipy.optimize import linprog pv_scenarios, probs = generate_pv_scenarios() S = len(pv_scenarios) n_q = I * T n_v = I * T n_spill = I * T n_r = S * T # 方便索引的函数 def qi(i, t): return i * T + t # 第i个电站第t时段的发电流量索引 def vi(i, t): return n_q + i * T + t # 第i个电站第t时段的库容变量索引 def spi(i, t): return n_q + n_v + i * T + t def ri(s, t): return n_q + n_v + n_spill + s * T + t n_vars = n_q + n_v + n_spill + n_r # 目标系数：最大化期望可消纳电量，所以linprog取负 c = np.zeros(n_vars) for s in range(S): for t in range(T): c[ri(s, t)] = -probs[s] * 1.0 # 时段长度取1小时，所以Δt=1下面拼等式约束。第一组是水量平衡，第二组是首末库容方程。
A_eq = [] b_eq = [] # 1. 水量平衡：V(i,t+1) - V(i,t) = Qin(i,t) + Qup(i,t) - q(i,t) - spill(i,t) # 这里忽略水流时滞，上游出库直接进入下游 for i in range(I): q_up = 0 if i == 0 else f"q_{i-1}_t + spill_{i-1}_t" for t in range(T - 1): row = np.zeros(n_vars) row[vi(i, t + 1)] = 1.0 row[vi(i, t)] = -1.0 row[qi(i, t)] = 1.0 row[spi(i, t)] = 1.0 # 上游出库对下游是入流，在上游自己的水量平衡里已经扣掉， # 所以在下游模型中作为正项加入 if i > 0: row[qi(i - 1, t)] -= 1.0 row[spi(i - 1, t)] -= 1.0 A_eq.append(row) b_eq.append(stations[i]["Qin"][t]) # 2. 首时段库容固定：V(i,0) = V0 # 末时段库容固定：V(i,T) = Vend for i in range(I): row_init = np.zeros(n_vars) row_init[vi(i, 0)] = 1.0 A_eq.append(row_init) b_eq.append(stations[i]["V0"]) row_end = np.zeros(n_vars) row_end[vi(i, T - 1)] = 1.0 A_eq.append(row_end) b_eq.append(stations[i]["Vend"])注意这里我简化了水量平衡的写法，没有把矢量化的流量单位换算写进去。实际使用时，库容和流量单位必须统一。如果库容单位是万立方米，流量单位是立方米每秒，那每小时流量体积是3.6千立方米，即0.36万立方米。这个换算经常被忽略，一旦单位搞错，整个结果都不合理。
不等式约束分四组：库容上下限、发电流量上下限、弃水非负、可消纳功率与系统出力约束。
A_ub = [] b_ub = [] # 3. 发电流量上下限 for i in range(I): for t in range(T): row = np.zeros(n_vars) row[qi(i, t)] = 1.0 A_ub.append(row) b_ub.append(stations[i]["qmax"]) row = np.zeros(n_vars) row[qi(i, t)] = -1.0 A_ub.append(row) b_ub.append(-stations[i]["qmin"]) # 4. 库容上下限 for i in range(I): for t in range(T): row = np.zeros(n_vars) row[vi(i, t)] = 1.0 A_ub.append(row) b_ub.append(stations[i]["Vmax"]) row = np.zeros(n_vars) row[vi(i, t)] = -1.0 A_ub.append(row) b_ub.append(-stations[i]["Vmin"]) # 5. 弃水非负 for i in range(I): for t in range(T): row = np.zeros(n_vars) row[spi(i, t)] = -1.0 A_ub.append(row) b_ub.append(0.0) # 6. 可消纳功率约束：R(s,t) - sum_i P_i(t) <= PV(s,t) # 以及 R(s,t) <= C(t) for s in range(S): for t in range(T): row_r = np.zeros(n_vars) row_r[ri(s, t)] = 1.0 for i in range(I): # P_i(t) = K_i * H_i * q_i(t) row_r[qi(i, t)] -= stations[i]["K"] * stations[i]["H"] A_ub.append(row_r) b_ub.append(pv_scenarios[s, t]) row_c = np.zeros(n_vars) row_c[ri(s, t)] = 1.0 A_ub.append(row_c) b_ub.append(C_limit[t])然后调用linprog求解：
res = linprog( c, A_ub=np.array(A_ub), b_ub=np.array(b_ub), A_eq=np.array(A_eq), b_eq=np.array(b_eq), method="highs" ) if res.success: x = res.x q_opt = x[:n_q].reshape(I, T) v_opt = x[n_q:n_q + n_v].reshape(I, T) spill_opt = x[n_q + n_v:n_q + n_v + n_spill].reshape(I, T) r_opt = x[n_q + n_v + n_spill:].reshape(S, T) print("求解成功，目标函数值（期望可消纳电量）:", -res.fun) else: print("求解失败:", res.message)3.4 结果输出与可视化要点
计算完成后，至少要画三张图：各电站发电流量过程、库容过程、系统外送功率与光伏出力的对比。下面是简单的可视化代码。
import matplotlib.pyplot as plt time = np.arange(1, T + 1) plt.figure(figsize=(12, 4)) for i in range(I): plt.plot(time, q_opt[i], label=stations[i]["name"] + "发电流量") plt.xlabel("时段/h") plt.ylabel("流量/(m³/s)") plt.legend() plt.grid(True) plt.title("梯级电站发电流量优化结果") plt.show() # 期望可消纳功率 r_expect = np.zeros(T) for s in range(S): r_expect += probs[s] * r_opt[s] plt.figure(figsize=(12, 4)) plt.plot(time, r_expect, label="期望可消纳功率") plt.plot(time, pv_forecast, label="光伏预测出力") plt.plot(time, C_limit, label="外送通道上限", linestyle="--") plt.xlabel("时段/h") plt.ylabel("功率/MW") plt.legend() plt.grid(True) plt.title("系统可消纳功率与光伏出力对比") plt.show()可视化这一步容易被人忽略，但调试模型时，图比数值更能暴露问题。比如发电流量曲线如果出现相邻时段跳变异常，多半是爬坡约束缺失；库容曲线如果长期贴着上限，可能是弃水变量权重太低，模型选择让水库蓄满然后弃水，而不是充分利用水量。
4. 实操过程与参数调优实录
4.1 从文献参数到可运行数据的转换
文献里的水库参数往往是一张大表，列着正常蓄水位、死水位、总库容、调节库容、装机容量等。复现时不能照搬，需要整理成模型需要的形式。
最容易被卡住的是库容单位。很多文献的库容单位是亿立方米，而发电流量单位是立方米每秒，水量平衡方程要统一成同一个单位。我建议全部先换算成“立方米每秒·小时”这种水量单位，即把库容除以3600秒，得到等效的“立方米每秒·小时”数，然后和逐时段的流量数值直接相加减。这样就不容易出现单位错误。
还有一个常用技巧：如果文献只给了正常蓄水位和死水位对应库容，没有给库容曲线，先用两点线性插值估算每个库容对应的水位，不追求精确。短期调度中水库不一定会跨越很大库容范围，线性近似影响不大。
出力系数K_i需要根据装机容量反推。装机容量N_i、最大发电流量q_i^max、平均水头H_i的关系是N_i ≈ K_i * H_i * q_i^max，取K_i=9.8，可以解出合适的H_i，或者反过来给定水头，用实际额定数据校验。我在复现中发现，文献给的装机容量和水头往往对不上，最终微调K_i或H_i让最大出力等于装机容量，结果曲线更合理。
4.2 求解器选择与运行时长实测
初始版本我用的是method="simplex"，在小规模场景下能跑，但场景数增加到30以后求解时间会明显上升，还容易报数值警告。后来改成method="highs"，求解速度几乎提升了一个量级。
实测数据供参考：3个梯级电站、24个时段、20个光伏场景，变量数约1400，约束数约3200，HiGHS求解时间在0.15秒左右。把场景提升到50个，变量数约2000，求解时间也就0.3秒。这个速度对短期调度模型来说非常理想，完全可以把场景数加大到50甚至100，换取更稳定的期望值估计。
如果发现求解时间超过10秒，先检查你的矩阵构造方式是不是用了Python循环逐行append。变量数超过5000时，这种逐行构造方式会消耗大量时间在list和numpy转换上。我的建议是先用列表收集A_eq和b_eq，最后一次性np.array，不要反复在循环里调用np.vstack，那个内存开销非常可怕。
4.3 灵敏度分析与参数调优建议
复现论文时，除了跑出一个结果，还要做几条灵敏度曲线增强说服力。
一是外送通道上限C_t的影响。固定其他参数，改变C_t从600MW到1000MW，观察期望可消纳电量的变化。通道越宽，可消纳电量越大，但增长会有饱和趋势，因为水库发电流量有上限。这个曲线能直观体现“消纳瓶颈”在哪里。
二是末水位约束的影响。把末库容从预设值上下调10%，观察发电量变化。末库容越松，模型越能在高电价时段（光伏低、负荷高）多放水，但会影响下一个调度周期的运行。论文里常把这个写成“末水位对调度结果的影响分析”。
三是光伏预测误差标准差的影响。标准差从5%增加到20%，期望可消纳电量通常下降，因为系统需要预留更多调节能力应对不确定性。这一条曲线很有价值，可以用来说明随机优化的必要性。
实操建议：每调一个参数，固定随机种子seed，保证不同参数之间的对比不受场景随机性的干扰。场景生成函数里的seed=42一定要保留，否则每次运行场景都不一样，灵敏度曲线会变得毛糙。
5. 常见问题与排查技巧实录
5.1 不可行解与收敛困难的排查
这是我第一次复现时卡得最久的问题。linprog返回“The problem is infeasible”，但硬是看不出哪里矛盾。后来一步步排查，发现是末库容和初始库容设成同一个值，而中间入库流量总量不够，导致无论如何都回不到初始库容。解决办法是先不固定末库容，只加一个末库容≥某值的约束，跑通后再收紧成等式。
还有一种很隐蔽的矛盾：上游电站qmax过大，而下游库容小，上游突然放水导致下游库容越限。这时要在水量平衡约束里考虑水流时滞和传播时间，或者给下游库容加一个更大的Vmax。如果文献没有水流时滞，就默认上游出库在同一个时段到下游，这样耦合最紧，模型偏保守但也安全。
还有一个常用技巧：把目标函数系数全部乘以一个较大的缩放因子，比如100，让数值范围更接近1~100。linprog对目标系数尺度不太敏感，但对约束矩阵的尺度非常敏感。如果b_eq里有100、也有0.01，数值尺度差距过大，求解器容易出现数值病态。我的经验是所有约束的右侧尽量统一在一个量级，必要时把库容单位换算后就不会有太大差异。
5.2 运行时间过长怎么办
如果场景数量已经很大，比如200个，求解时间可能接近1秒到几秒，虽然还能接受，但如果你要做100组灵敏度分析，累积时间也不小。优化手段有三个：
- 降低场景数：用K-means削减到10~15个，期望值依然相对稳定。
- 压缩问题尺寸：把每个场景的可消纳功率R_{s,t}用解析方法代入。实际上因为R只出现在目标函数和两个上限约束里，可以消去R变量，让每个场景的约束只作用于水电站发电流量上，但代价是目标函数变成非线性凸函数，需要做分段线性化。复现时没必要，直接用R变量最简单。
- 使用稀疏矩阵：A_eq和A_ub转换成scipy.sparse.csr_matrix，可以显著降低内存占用。变量上5000之后建议这么做。
我亲测过，用HiGHS求解时，稀疏矩阵比稠密矩阵在2000变量规模下没有太大速度差异，所以代码里保持简单也是可以的。
5.3 结果水位越界、出力波动的处理
优化完成后检查库容曲线，偶尔会发现某个库容值正好卡在Vmax上，相邻时段出现尖角。这是因为线性水量平衡没有限制库容变化速率，模型可以让库容快速上升。工程上水库蓄泄水速率受闸门限制，短期模型通常忽略，但论文讨论里要提一句。
出力波动则更多来自光伏场景末端的极端值。当某个场景的光伏出力突然从低谷跳到高峰，系统可消纳功率R上限被PV抬高，但水电站出力调节不过来，整体外送曲线会跟着跳。处理办法是在目标函数里加一个小的二次惩罚？但线性规划不支持。更实用的方法是先对光伏场景做平滑滤波，或者对系统出力加爬坡约束。梯级水电短期调度一般会对“系统总出力变化速率”加以限制：
[
- \Delta P^{max} \le (PV_{s,t}+\sum_i P_{i,t}) - (PV_{s,t-1}+\sum_i P_{i,t-1}) \le \Delta P^{max} ]
但这样每个场景都会增加约束，场景多时约束规模跟着翻倍。实地场景数20左右时没问题。
最后说一个我在复现过程中特别有体会的细节：目标函数里的“期望”并不神秘，它本质上就是场景概率加权平均。很多人一开始纠结要不要用随机规划、鲁棒优化、机会约束等高级方法，其实EI论文里很多模型就是把光伏场景和线性规划结合，用期望值做目标，工程意义已经足够。把基础版本跑通，再往里面加风险偏好、备用容量、日内滚动，都是后话。
如果你也是第一次复现这类梯级水光互补调度模型，建议按这个顺序走：先把确定性模型（只有一个光伏预测场景）跑通，确认水量平衡和末库容逻辑正确；再加场景生成和期望目标；最后加复杂约束。不要一上来就把50个场景、爬坡约束、生态流量全塞进去，那样只会让排查问题变成灾难。模型本身不复杂，复杂的是数据和调试流程。把基础版本做稳，后面扩展完全水到渠成。
