企业官网建设流程全解析

1. 项目背景：为什么梯级水光互补调度值得复现

最近在复现一个EI期刊里的短期优化调度模型：梯级水光互补系统最大化可消纳电量期望，用Python代码实现。这个题目听起来长，但拆开看就是几个关键词的组合：梯级水电站群、光伏电站、短期调度（比如日前24小时）、目标函数是“期望可消纳电量最大化”。做电力系统优化的人应该不陌生，这是新能源消纳问题的一个典型场景。

先说清楚这项目能解决什么问题。传统调度中，水电和光伏各自单独调度，容易出现光伏大发时段水电站没有主动压出力，或者光伏出力波动导致弃光；反过来，如果水电站死板地按计划出力，又可能占用了光伏的消纳空间。梯级水光互补调度的核心思想，是让梯级水电站通过调节自身出库流量，去“补偿”光伏出力的随机波动，使得整个系统在满足负荷/外送通道限制的前提下，尽量多地吸收水电和光伏电量。这里的“最大化可消纳电量期望”，就是指在考虑光伏出力随机性的情况下，让系统能够消耗掉的（而不是弃掉/限掉的）电量，在期望意义下最大。

这个项目适合谁看？如果你在做水电优化调度、新能源消纳、或者刚接触Python建模优化，都挺合适。我会从模型原理讲起，给出数据构造、优化代码框架、画图结果分析，最后把我在复现中踩过的坑列出来。整个过程不依赖商业求解器，用Scipy和Numpy就能跑通，方便你在一台普通笔记本上完成复现。

2. 模型核心拆解：目标函数、随机性与梯级约束的数学表达

2.1 目标函数：从“期望”到可计算的场景均值

“最大化可消纳电量期望”这句描述里，“期望”二字意味着光伏出力不是确定的，而是随机的。在短期调度中，我们通常在日前做决策，但光伏出力在未来某一小时的实际值是不确定的。期望最大化，指的是在所有可能的光伏场景下，目标函数的平均值最大。

数学上可以写成：

[ \max \ E\left[\sum_{t=1}^{T} A_t(\omega)\right] ]

其中 (A_t(\omega)) 表示在随机场景 (\omega) 下，时段 (t) 的可消纳电量。可消纳电量直观理解就是系统能够接纳的电力总量，如果总有功（水电+光伏）超过消纳上限（比如负荷需求或外送通道容量），超出部分就会被削减，所以可消纳电量是：

[ A_t = \min(D_t,\ P_h(t) + P_pv(t,\omega)) ]

其中 (D_t) 是时段 (t) 的消纳能力（如负荷）， (P_h) 是梯级水电总出力， (P_pv) 是光伏出力。

期望怎么算？实践中最常用的是蒙特卡洛场景法：生成 (S) 个光伏出力场景，把目标函数变成所有场景的平均值：

[ \max \frac{1}{S}\sum_{s=1}^{S}\sum_{t=1}^{T} \min(D_t,\ P_h(t)+P_pv(t,s)) ]

注意一个关键点：优化变量是水电站的出库流量（或者出力）序列，这个序列是日前确定的，不能随着每个光伏场景改变。所以目标函数中，(P_h(t)) 在所有场景下是同一套决策变量，而 (P_pv(t,s)) 随场景变化。这样才体现“鲁棒/期望”的意义——我们不知道明天光伏到底是多少，但根据下游负荷需求和光伏概率分布，先定好水电的调度计划，让平均消纳电量最大。

2.2 梯级水电站的物理约束长什么样

梯级水电站和单站最大的不同，在于上下游之间有水量联系。上游电站的出库流量（包括发电流量和弃水流量），经过一定的滞时（比如1~2小时），会成为下游电站的入库流量的一部分。这种水力耦合关系如果不建模，调度计划根本没有物理可实现性。

拿两级梯级电站距离，t时段的水量平衡方程：

[ V_{i,t+1} = V_{i,t} + (Q_{in,i,t} - Q_{out,i,t} - S_{i,t}) \cdot \Delta t ]

• (V_{i,t})：第i级水库在第t时段末的库容；
• (Q_{in,i,t})：入库流量，对上游水库来说就是天然径流，对下游水库来说天然径流加上上游水库出库流量滞后；
• (Q_{out,i,t})：发电流量（用于水轮机组发电的流量）；
• (S_{i,t})：弃水流量（不发电但需要下泄的流量）。

对于上游水库 (i=1)，(Q_{in,1,t} = I_{1,t})（天然来水）；对于下游水库 (i=2)，(Q_{in,2,t} = I_{2,t} + Q_{out,1,t-\tau} + S_{1,t-\tau})，这里 (\tau) 是水流滞时。如果时间尺度是小时，滞时通常取0~2小时。

水电站出力方程是另一个关键。常用的简化表达：

[ P_{h,i,t} = 9.81 \cdot \eta_i \cdot H_i(t) \cdot Q_{out,i,t} / 1000 ]

单位：(H) 是水头（米），(Q) 是发电流量（立方米/秒），(\eta) 是综合效率。如果水头变化不大，可以用平均水头 (H_{avg}) 替代，把9.81etaH/1000打包成常数 (k_i)，于是出力线性化为：

[ P_{h,i,t} = k_i \cdot Q_{out,i,t} ]

这样一来，优化模型就能变成一个带线性/非线性约束的规划问题，后面用Scipy求解方便很多。

还有一组基本约束：

• 库容上下限：(V_{i,\min} \le V_{i,t} \le V_{i,\max})
• 发电流量上下限：(Q_{i,\min} \le Q_{out,i,t} \le Q_{i,\max})
• 弃水流量非负：(S_{i,t} \ge 0)
• 出力上限：(P_{h,i,t} \le P_{i,\max})
• 末端库容要求（调度周期末水位不低于某个值）：(V_{i,T+1} \ge V_{i,end})

这些约束尽量都要放进复现代码里，否则跑出来的结果只能算玩具。但如果某些非线性约束导致难收敛，可以做线性化处理。后面我会给一个简化但逻辑完整的版本。

2.3 光伏随机性的处理：场景生成与等价转化

光伏出力的随机性来源于云层遮挡、温度变化、辐照度波动等。在日前调度中，我们通常有光伏出力的预测曲线（比如基于数值天气预报得到），但真实出力会在预测值附近波动。常见的建模方法是假设光伏出力服从某种概率分布，然后通过随机抽样生成场景。

最简单的场景生成方式是对预测值叠加正态分布随机扰动：

[ P_pv(t,s) = \max\left(0,\ \bar{P}{pv}(t) + \epsilon{t,s}\right) ]

其中 (\epsilon_{t,s} \sim N(0, \sigma_t^2))，(\sigma_t) 反映预测误差。也可以使用拉丁超立方抽样，保持场景多样性。注意 (\epsilon) 不能太大，否则会出现负的光伏出力，需要截断。

生成后的 (S) 个场景 (P_pv(t,s)) 放进目标函数平均，就把随机期望问题转化为一个确定性的场景平均问题。理论上场景数越多，期望估计越准，但计算量线性上升。一般可以先取50个场景，后续可以下降。

这里还有一个技巧：如果用Scipy的optimize，目标函数里要求导或差分会很慢。建议把所有场景的循环写成向量化运算，用Numpy的二维数组一次算完，速度快很多。

2.4 简化假设与合理性

在复现论文时，完全还原论文所有细节往往不现实。论文里可能还考虑了水头动态变化、机组组合、最小运行时间、爬坡速率等。复现第一步应该搭一个“骨架”模型，验证核心目标函数和梯级水量平衡逻辑是否跑得通，然后逐步添加复杂度。

我在本项目中做了一些简化：

• 用平均水头系数代替水头-库容非线性关系；
• 不考虑机组启停，把每个水电站看作一台等效机组；
• 不设置爬坡约束，后续可以扩展；
• 弃水流量作为决策变量，但只在需要维持库容安全时才会出现；
• 假设负荷曲线已知，光伏场景通过随机扰动生成。

这些简化不影响对核心机制的验证。当你拿到了基准结果，再一步步把论文里的复杂约束加上去，就是完整复现的路径。

3. Python实现方案与代码架构

3.1 环境准备与工具选型

Python生态里做优化调度，工具选择因人而异。这里我说一下我的选择：

• Numpy：所有向量和矩阵运算；
• Scipy.optimize.minimize：SLSQP算法，支持等式和不等式约束；
• Pandas：处理时间序列数据；
• Matplotlib：画调度结果曲线。

为什么不直接上Gurobi或Cplex？因为很多读者没有商业求解器授权，而且线性化后的模型用SLSQP也能解。但如果你之后要做大规模或精确求解，建议把模型转成线性规划，用Gurobi的Python接口来解。

安装环境也很简单，Anaconda自带Numpy、Pandas、Matplotlib，只用额外确认Scipy版本。我的测试环境是Python 3.10，Scipy 1.9.1，没有任何兼容问题。

3.2 数据结构设计

建模前先把所有参数和时间序列组织好。推荐用Python字典来存储水电站参数，比如：

plants = { 'up': { 'k': 0.070, # 出力系数 MW/(m3/s) 'Vmax': 800, # 库容上限 10^4 m3 'Vmin': 150, 'V0': 500, 'V_end': 400, # 调度周期末约束 'Qmax': 120, # 最大发电流量 m3/s 'Qmin': 10, 'Smax': 50, # 最大弃水流量 'inflow': inflow_up, # T维数组，天然来水 }, 'down': { 'k': 0.045, 'Vmax': 1200, 'Vmin': 300, 'V0': 800, 'V_end': 600, 'Qmax': 200, 'Qmin': 20, 'Smax': 80, 'inflow': inflow_down, 'lag': 1, # 上游到下游的滞时 } }

时间尺度取T=24小时。库容单位统一为“万立方米”，流量单位“立方米/秒”，二者相乘再乘时间步长（秒）得到库容变化，要仔细换算。一小时是3600秒，所以：

[ \Delta V = (Q_{in} - Q_{out} - S) \times 3600 / 10^4 \quad \text{(单位：万m3)} ]

这个单位换算很容易出错，后面我会专门讲。

负荷曲线和光伏预测曲线也需要准备好：

load = np.array([...]) # 24个时段，单位MW pv_forecast = np.array([...]) # 24个时段，单位MW

然后生成场景矩阵，形状是(S, T)：

S = 20 noise_sigma = 0.5 * pv_forecast # 简单取预测值的50%作为标准差 scenarios = np.zeros((S, T)) for s in range(S): scenarios[s] = np.maximum(0, pv_forecast + np.random.normal(0, noise_sigma))

实际还可以用拉丁超立方抽样，但在这个复现里够用了。

3.3 核心优化模型建模（scipy minimize + SLSQP）

决策变量需要仔细定义。假设有两个水电站，每个时段有发电流量 (Q_{out}) 和弃水流量 (S)。那么总变量数 = 24 * 2 * 2 = 96。变量排列顺序可以是：

• x[0:24]：上游发电流量
• x[24:48]：上游弃水流量
• x[48:72]：下游发电流量
• x[72:96]：下游弃水流量

目标函数要接收决策变量，并返回负期望可消纳电量（因为minimize是求最小）。

def objective(x): q_up = x[0:24] s_up = x[24:48] q_down = x[48:72] s_down = x[72:96] p_up = k_up * q_up p_down = k_down * q_down # 这里用向量化计算所有场景的消纳电量 # p_total 形状 (S, T)，需要广播 p_total = p_up + p_down + scenarios # scenarios (S,T) absorb = np.minimum(load, p_total) # 场景下每个时段可消纳电量 return -np.mean(absorb.sum(axis=1))

约束条件分两类。

第一类是等式约束，即水量平衡。水量平衡是时序递推的，所以得手写每个时段的递推公式，不能直接用向量化：

def water_balance_eq(x): q_up = x[0:24]; s_up = x[24:48] q_down = x[48:72]; s_down = x[72:96] qout_up = q_up + s_up qout_down = q_down + s_down V_up = np.zeros(25) V_down = np.zeros(25) V_up[0] = plants['up']['V0'] V_down[0] = plants['down']['V0'] for t in range(24): # 上游水量平衡 V_up[t+1] = V_up[t] + (plants['up']['inflow'][t] - qout_up[t]) * 3600 / 1e4 # 下游水量平衡：上游出库考虑滞时 inflow_down_t = plants['down']['inflow'][t] if t - plants['down']['lag'] >= 0: inflow_down_t += qout_up[t - plants['down']['lag']] V_down[t+1] = V_down[t] + (inflow_down_t - qout_down[t]) * 3600 / 1e4 # 把所有时段的库容约束偏差都返回到一个数组? SLSQP需要多个约束 # 这里仅返回最终末库容约束和所有库容越限量之和 ...

但是SLSQP要求等式约束是形如cons_eq['fun'] = 0。水量平衡本身是递推的，我们更应该把它表达成一组等式约束，每个时段一个等式。做法是让约束函数返回一个数组，长度是2*24：

def water_balance(x): q_up = x[0:24]; s_up = x[24:48] q_down = x[48:72]; s_down = x[72:96] V_up = np.zeros(25); V_down = np.zeros(25) V_up[0] = V0_up; V_down[0] = V0_down for t in range(24): V_up[t+1] = V_up[t] + (inflow_up[t] - q_up[t] - s_up[t]) * 3600 / 1e4 inflow_d = inflow_down[t] if t-lag >=0: inflow_d += q_up[t-lag] + s_up[t-lag] V_down[t+1] = V_down[t] + (inflow_d - q_down[t] - s_down[t]) * 3600 / 1e4 # 返回每个水库每个时段的等式约束值：V[t+1] 应等于递推值，所以这里的偏差为零？但这样检查不出来 # 其实我们是在递推，V是该时段的库容，它由上一个时段决定，这里等式约束已经自动满足。 # 真正的约束是库容上下限不等式，可以在边界约束中体现。 # 末库容约束可以单独写等式约束。 cons = [] # 末库容 cons.append(V_up[-1] - V_end_up) cons.append(V_down[-1] - V_end_down) return cons

这里要说明一下：水量平衡其实是通过递推公式“硬编码”在代码中，而不是作为寻优变量来改变，因此不需要写成等式约束。变量只有流量，库容由流量序列通过递推唯一确定。这样处理的好处是大大减少变量维度，坏处是目标函数里必须按顺序递推。SLSQP支持这种隐式约束。

库容上下限可以作为不等式约束，在每个时段检查库容是否越限：

def volume_bounds(x): # 递推库容，返回V_up[1:] - Vmax, Vmin - V_up[1:] 等。 ... v_up = ... (长度为24) v_down = ... c1 = v_up - Vmax_up c2 = Vmin_up - v_up c3 = v_down - Vmax_down c4 = Vmin_down - v_down return np.concatenate([c1, c2, c3, c4])

SLSQP的不等式约束要求fun(x) >= 0。所以我们让返回值为Vmax - V和V - Vmin都大于零。

最终约束定义：

cons = [ {'type': 'ineq', 'fun': volume_bounds}, {'type': 'eq', 'fun': end_volume_constraint}, {'type': 'ineq', 'fun': lambda x: plants['up']['Qmax'] - x[0:24]}, ... ]

边界约束直接放在bounds里：

bounds = [] for _ in range(24): bounds.append((Qmin_up, Qmax_up)) bounds.append((0, Smax_up)) for _ in range(24): bounds.append((Qmin_down, Qmax_down)) bounds.append((0, Smax_down))

这样整个优化问题就构建好了。求解一行代码：

res = minimize(objective, x0, method='SLSQP', bounds=bounds, constraints=cons, options={'maxiter': 500})

x0可以用所有发电流量取中间值，弃水为零。

3.4 关键代码片段与参数计算

这里给出一段可运行的核心循环框架，你可以直接复制再改参数：

import numpy as np from scipy.optimize import minimize import matplotlib.pyplot as plt T = 24 S = 20 np.random.seed(42) # 假设负荷曲线 load = np.array([50,48,46,45,46,50,55,60,65,68,70,72,71,70,68,66,65,64,63,60,58,55,53,50], dtype=float) # 光伏预测（白天多，晚上零） pv_forecast = np.zeros(T) pv_forecast[8:18] = [20,35,50,65,70,72,68,55,35,15] # 生成场景 sigma = 0.2 * pv_forecast + 2 scenarios = np.zeros((S, T)) for s in range(S): scenarios[s] = np.maximum(0, pv_forecast + np.random.normal(0, sigma)) # 水电站参数 k_up=0.070; k_down=0.045 V0_up=500; Vmax_up=800; Vmin_up=150; V_end_up=400 V0_down=800; Vmax_down=1200; Vmin_down=300; V_end_down=600 Qmin_up=10; Qmax_up=120; Qmin_down=20; Qmax_down=200 Smax_up=50; Smax_down=80 lag=1 inflow_up=np.array([80]*T) inflow_down=np.array([20]*T) def update_volume(x): q_up = x[0:24]; s_up = x[24:48] q_down = x[48:72]; s_down = x[72:96] V_up=[V0_up]; V_down=[V0_down] for t in range(T): v_up = V_up[-1] + (inflow_up[t] - q_up[t] - s_up[t]) * 3600 / 1e4 V_up.append(v_up) inflow = inflow_down[t] if t-lag >= 0: inflow += q_up[t-lag] + s_up[t-lag] v_down = V_down[-1] + (inflow - q_down[t] - s_down[t]) * 3600 / 1e4 V_down.append(v_down) return np.array(V_up[1:]), np.array(V_down[1:]) def objective(x): q_up = x[0:24]; q_down = x[48:72] p_up = k_up * q_up p_down = k_down * q_down p_total = p_up + p_down + scenarios absorb = np.minimum(load, p_total) return -np.mean(absorb.sum(axis=1)) def volume_cons(x): V_up, V_down = update_volume(x) return np.concatenate([Vmax_up - V_up, V_up - Vmin_up, Vmax_down - V_down, V_down - Vmin_down]) def end_cons(x): V_up, V_down = update_volume(x) return np.array([V_up[-1] - V_end_up, V_down[-1] - V_end_down]) cons = [ {'type':'ineq','fun': volume_cons}, {'type':'eq','fun': end_cons}, {'type':'ineq','fun': lambda x: Qmax_up - x[0:24]}, {'type':'ineq','fun': lambda x: x[0:24] - Qmin_up}, {'type':'ineq','fun': lambda x: Qmax_down - x[48:72]}, {'type':'ineq','fun': lambda x: x[48:72] - Qmin_down}, {'type':'ineq','fun': lambda x: Smax_up - x[24:48]}, {'type':'ineq','fun': lambda x: x[24:48]}, {'type':'ineq','fun': lambda x: Smax_down - x[72:96]}, {'type':'ineq','fun': lambda x: x[72:96]}, ] bounds = [] for _ in range(T): bounds.append((Qmin_up, Qmax_up)) bounds.append((0, Smax_up)) for _ in range(T): bounds.append((Qmin_down, Qmax_down)) bounds.append((0, Smax_down)) x0 = np.zeros(4*T) x0[0:24] = 60 x0[48:72] = 80 res = minimize(objective, x0, method='SLSQP', bounds=bounds, constraints=cons, options={'maxiter': 800, 'ftol': 1e-8}) print('成功：', res.success, res.message) print('最佳期望可消纳电量：', -res.fun)

这段代码我在本地跑过，可以收敛。但参数不同可能会失败，后面会讲调试方法。

4. 实操过程：从数据到结果的完整复现步骤

4.1 构造仿真场景：水电站系数与光伏场景

第一步是先确认单位统一。很多复现跑不出结果，都是因为流量单位用了立方米每秒，而水库库容单位用了亿立方米，导致水量平衡数值差好几个数量级。上面的代码里我把库容单位定为万立方米，流量乘以3600秒再除以1e4，换算成万立方米，这样V的数量级在数百左右，和目标值匹配。

水电站出力系数k的计算：假设水头为50米，效率0.85，则 (k = 9.81 \times 0.85 \times 50 / 1000 = 0.4169) MW/(m3/s)。在实际测试中，为了让出力曲线接近负荷，我调低到0.070和0.045，因为上游水头和流量偏小。你可以根据实际情况调整。

光伏场景生成时，噪声的标准差如果取预测值的20%，白天波动可能超过10MW，晚上预测值为0，再加常数2MW，体现晚上可能有的出力误差。场景数量先取20，运行时间比较友好。

4.2 优化求解与结果提取

SLSQP对初始值很敏感，如果初值离可行域太远，可能无法收敛。我的经验是：发电流量初值取最大和最小流量之间的中间值，弃水取0；不用随机初值。如果还是无法收敛，可以先关闭末端库容等式约束，得到结果后再逐步加入。

求解结束后，检查res.success。如果success为False，不要急着改参数，先打印res.message，一般会提示“Singular matrix C”或者“Inequality constraints incompatible”。前者通常是目标函数或约束在初始点不可微，后者是约束之间矛盾。下面我会单独讲常见问题。

4.3 结果可视化与调度曲线解读

结果提取的关键是把决策变量转换成方便看的pandas表格：

q_up = res.x[0:24]; s_up = res.x[24:48] q_down = res.x[48:72]; s_down = res.x[72:96] V_up, V_down = update_volume(res.x) p_up = k_up * q_up p_down = k_down * q_down

然后画4个子图：

1. 水电出力与负荷曲线；
2. 库容变化曲线；
3. 发电流量与弃水流量；
4. 某个有代表性场景下的实际总出力（水电+光伏场景均值）与可消纳电量对比。

画图用Matplotlib两行代码即可：

plt.figure(figsize=(12,8)) plt.subplot(2,2,1) plt.plot(range(24), p_up, label='up hydro') plt.plot(range(24), p_down, label='down hydro') plt.plot(range(24), load, 'k--', label='load') plt.legend() plt.subplot(2,2,2) plt.plot(range(24), V_up, label='up reservoir') plt.plot(range(24), V_down, label='down reservoir') plt.legend() # ... plt.tight_layout() plt.savefig('schedule.png', dpi=150)

观察结果时重点看两点：总出力曲线是否尽量贴住负荷曲线（说明消纳能力强）；库容曲线是否滑落过猛导致末端库容恰好压在下限上。如果末端库容落在下限上，说明可用库容已经用到极限，可能最优解受末水位约束限制。

5. 常见问题与避坑指南

5.1 求解器不收敛怎么办

这是复现里遇到最多的坑。我踩过几次后，总结出下面这张排查表：

另外SLSQP对目标函数光滑性有要求，如果目标函数里用了np.minimum(load, p_total)，在负荷交叉点有不可导，但问题不大。如果严格有问题，可以用平滑化函数近似。

5.2 单位与量纲错误

单位错误会让你得到荒谬的库容值。我在水量平衡里把3600/1e4写错过多次。建议每一步都print库容序列，检查数值变化幅度是否合理。比如流量落差为50m3/s持续1小时，库容变化是50*3600/10000=18万立方米，对于最小库容150万立方米来说，变化比例约12%，属于正常范围。如果库容变化是几千，那一定是除以1e6而不是1e4。

5.3 约束冲突与不可行解

当库容上限收紧时，可能没有任何解决方案能让所有时段的流量都在允许范围内，同时满足末端库容。此时不要硬调求解器，而是放宽一个条件。比如把Vmin下调、把Qmax调大、或者把V_end调低。做模型实验时，我们往往需要先跑一个宽松基准，再逐步收紧。

5.4 随机场景数量怎么选

场景数太少，期望估计方差大；场景数太多，目标函数计算会拖慢优化。我用20个场景在24时段、96个变量的问题上，单次目标函数计算约在毫秒级，SLSQP迭代数百次也能接受。如果模型加了机组组合和爬坡约束，场景数要降到5个，或者做场景削减（如k-means聚类）。这个项目是复现论文，通常论文会写清楚场景削减方法，复现时先不做也影响不大。

5.5 性能与编码技巧

目标函数里的np.minimum(load, p_total)会自动在(S,T)矩阵上广播，比for循环快很多。水量平衡递推则无法避免循环，所以24个时段的循环还好。如果T扩展到96个时段，建议用Numba或者Cython加速递推部分，否则SLSQP每次迭代都调用递推函数，整体时间会明显拉长。

还有一个小技巧：先固定一个场景（S=1）测试模型能否求解、结果是否有物理意义，再放大场景数。这能迅速把建模错误从随机性影响中分离出来。

6. 复现后的扩展空间与个人心得

把基础模型跑通之后，这个项目的扩展空间其实很大。我自己复现时，第一版也是上面这种平均水头+场景均值模型，跑通后分别加了水头动态、机组组合、爬坡约束，发现目标函数和调度曲线都会有一定变化，但核心逻辑不变。如果追求EI论文复现精度，一定要把论文里的公式和参数仔细对照，尤其是库容-水位曲线、水头-出力关系，这些数据通常来自真实电站，编码时不能照抄示意数据。

另一个值得尝试的扩展是把目标函数从期望最大化改成鲁棒最大化，即用最坏场景下的消纳电量来衡量壁垒，这样调度策略会更保守，更贴近实际运行人员对风险的偏好。代码改动不大，只需要把目标函数里的np.mean改成np.min。

最后说一点个人体会。这种调度优化模型虽然是学术复现，但想真正用于生产，还需要解决中长期来水预报、实时反馈控制、不确定集选取等问题。但作为理解“梯级水光互补”问题的入门代码，它足够让你看清决策变量、随机场景、水力耦合三者之间的相互作用。我建议你先把这个骨架代码跑通，然后尝试修改负荷曲线或者光伏场景分布，观察调度曲线如何自适应变化。只有亲手改过参数，你才算真正“吃透”这个模型。