水光互补多目标优化调度：基于NSGA-II的Python实现

前阵子在做一个水电与新能源联合调度的仿真项目，最让我头疼的一个问题恰好就是水光互补优化调度本身：水电和光伏明明在时间特征上非常互补，但把它们放进同一个模型里以后，目标却互相打架。单纯追求总发电量最大，调度出来的出力曲线波动剧烈，电网侧根本不敢接；只追求出力平稳，又白白丢掉大量光伏电量。后来我才意识到，问题的根源在于"目标个数"——这类调度本质上是一个典型的多目标优化问题，必须同时处理发电经济性、出力平稳性、弃光率等多个互相冲突的指标。于是我把目光转向非支配排序遗传算法，也就是常说的 NSGA-II，用 Python 从零实现了一套完整的水光互补多目标调度代码。这篇文章就把这整套东西摊开来讲：数学模型怎么建、NSGA-II 的核心机制是什么、Python 代码怎么组织、最终怎么从 Pareto 前沿里挑出真正能落地的调度方案。

1. 为什么水光互补调度必须拆成多目标来做

1.1 水电与光伏的出力特性恰好错位

很多人在接触水光互补时，第一反应是"这俩加起来不就行了吗"。实际上没那么简单，你得先理解两条出力的天然时间曲线。

光伏出力呈现出极强的昼夜周期性和天气波动性：白天十点到下午两三点是发电高峰，傍晚开始快速下跌，晚上出力直接归零；遇到云层遮挡，出力还能在半小时内上下波动 30% 以上。水电则恰恰相反，它可以通过水库调节把能量"攒起来"，在任意时刻灵活释放，尤其是具有日调节或季调节能力的水电站，响应速度可以做到分钟级。

这就是互补的基础：光伏出力高涨时，水电可以降低出力、蓄水；光伏出力骤降时，水电快速顶上。理想情况下，两者的合成出力会变成一条相对平稳、总量又可观的曲线。但注意，这只是理想情况——现实中，水电站有自己的水量平衡约束，水库不能无限蓄水，机组出力也有上下限，光伏预测还存在误差。所有这些约束加起来，使得"如何安排水电每一时段的出力"变成一个有边界条件的优化问题，而不是简单的加减法。

1.2 单目标模型顾此失彼的实测表现

我在项目里先后跑了三组试验：第一组只优化总发电量，第二组只优化出力平稳性，第三组用 NSGA-II 做多目标优化。下表是三组试验在同一个典型日下的模拟结果，数据是我用本地光伏典型出力曲线和固定来水序列跑出来的，主要看趋势，不追求精确：

单目标模型最大的问题是它把决策者逼到了墙角：当你把发电量作为唯一目标，算法一定会让水电机组拼命顶满，库容紧张的时候就出现弃水，光伏午间出力高峰时又可能出现外送通道拥堵，逼着你弃光；当你把平稳性作为唯一目标，算法又会倾向于让水电始终保持在一个低而稳定的出力，光伏高峰时期不再压水，这时弃光率飚升，总发电量也跟着下降。这个两难不是参数调不好，而是模型本身只给了一个极端方向。

1.3 "互补"的本质是把时间错位变成调度红利

所以，"互补"这个概念真正有价值的地方在于：它给了调度一个额外自由度。水电出力的每个时点安排，就相当于在时间轴上重新分配能量。光伏出力高的时候少发水电、蓄水，相当于把水能"存"到更晚的时段；光伏出力低的时候多放水发电，弥补出力低谷。这样既保证了总量，又削平了波动。

但问题也随之而来：水库的调节能力不是无限的。你不可能为了削峰填谷，让水电在不该出力的时候一点不发——那样可能造成水库水位顶到上限导致被迫弃水；也不可能让水电在低谷时段疯狂出力——那样很快会把库容放空，后面无米下锅。于是，调度的本质变成：在水量平衡、库容限制、机组爬坡限制三条硬约束之下，去手机版找一条权衡"发电量和波动性"的出力序列。这从数学上讲就是一个带约束的多目标优化问题，适合用基于 Pareto 支配关系的进化算法来求解。

2. 优化模型的建立：变量、目标与约束条件

2.1 调度时段与决策变量设计

建立模型第一步是确定时间尺度和决策变量。我在代码里默认以一天 24 个小时为调度周期，每小时一个时段，也就是 T=24；如果你手头的数据是 15 分钟一个点，把 T 改成 96 就行，代码骨架不需要大改。

决策变量选什么？在水光互补系统中，光伏出力是"看天吃饭"的外部输入，不能作为决策变量；真正能主动调节的只有水电机组的出力序列。所以我规定：

• 决策变量：水电机组每个时段的有功出力 $P_h(t)$，$t=1,2,\dots,T$
• 输入数据：光伏预测出力 $P_{pv}(t)$、水库来水序列 $I(t)$、初始库容 $V_0$
• 实值编码的染色体长度为 T，即一个个体就是一条完整的"日内水电出力曲线"

这里有个很容易被新手忽略的点：如果系统中包含多座梯级水电站，决策变量就变成"电站数量 × T"，每座电站都要有自己的出力序列。为了讲清楚核心思路，文章里我按单座水电站 + 光伏电站来叙述，多电站只需在此基础上把变量维度撑开。

2.2 两个核心目标函数

我把核心目标压缩成两个：系统总发电量最大，以及出力过程最平稳。这两个目标在物理意义上直接反映了"经济性"和"电网友好性"。

目标一，发电量最大：

$$F_1 = \sum_{t=1}^{T} [P_h(t) + P_{pv}(t)] \cdot \Delta t$$

目标二，出力波动最小。这里有几种取法，最常用的是相邻时段出力差的平方和，或者负方向取方差：

$$F_2 = \sum_{t=1}^{T-1} [P_{total}(t+1) - P_{total}(t)]^2$$

其中 $P_{total}(t) = P_h(t) + P_{pv}(t) - P_{curtail}(t)$，$P_{curtail}(t)$ 是该时段弃光功率。由于弃光只可能发生在光伏出力过大时，模型中我会用约束惩罚来限制它：如果水电出力下限导致总出力和光伏出力之和超过外送通道上限，就记一笔弃光惩罚。

实际项目中还可以加入第三个目标，比如"弃光量最小"或者"水库水位变化最小"，目标数量增加到三个以后 Pareto 前沿就从二维曲线变成三维曲面，可视化和选方案都会更复杂。所以我的建议是：先用两个核心目标把系统跑通，后续再按业务需求加维度。

2.3 必须满足的物理约束清单

多目标进化算法处理约束的方式不像线性规划那样干净，通常会通过罚函数或者可行性修复来完成。但不管用哪种方式，约束本身必须列清楚：

其中水量平衡是最容易出错的地方。很多初学者把发电流量直接当成出力除以一个固定转换系数 $\eta$，但忽略了一个事实：如果水库已经达到库容上限，而来水还在持续流入，那么多出来的水量必须以弃水 $S(t)$ 的形式排出，这部分能量没有发电。于是胡克定律一样，反推出来的库容变化和实际过程经常对不上。

2.4 为什么这个模型天然适合 NSGA-II

如果你尝试过用传统加权法解这个模型，你会发现它有两个致命问题。第一，目标函数中的波动项包含相邻时段差，整体问题是高度非线性的；第二，可行域被多条不等式约束切割得支离破碎，不是一个简单的凸集。这种情况下，加权法每次只能从一个权重组合中解得一个点，你想得到整条 Pareto 前沿就得反复试权重，而且一旦前沿非凸，加权法在凹段上的解永远找不全。

非支配排序遗传算法却完全不同。它采用种群搜索，每一代同时维护一组候选解，通过非支配排序把这些解划分成若干层次，推动种群整体向 Pareto 前沿逼近。它不要求模型可微、不要求凸性，只需要你能写清楚每个个体对应的目标函数值，本质上把优化问题变成"评估 + 进化"两个模块。这种特性让它在工程配置复杂的调度问题里非常受欢迎。这也是我选择它而不是粒子群或传统线性规划的原因：模型改一个约束，代码改动量很小，算法本身不需要跟着重构。

3. NSGA-II 的核心机制：非支配排序、拥挤度与精英保留

3.1 "支配"是什么，非支配前端怎么排

要理解 NSGA-II，先得理解多目标优化里"谁更优"的判断标准。两个目标的情况很好类比：甲方案发电量高而波动也大，乙方案发电量低但波动也小，这两个方案谁也不敢说完全压制对方。只有当甲方案的每个目标都不差于乙方案，并且至少有一个目标严格优于乙方案时，我们才说甲支配乙。

非支配排序做的事情，就是在一个种群内部把所有个体按"支配关系"分层：

• 第 1 层：种群中不被任何其他个体支配的个体，它们组成了当前最好的前端。
• 第 2 层：去掉第 1 层后，剩余个体中不被剩余个体支配的个体。
• 依此类推，得到一组从前到后的前端序列 $F_1, F_2, \dots$。

这个排序保证了算法优先保留"综合能力强"的个体。代码实现上，最直接的是两两比较的快速非支配排序，时间复杂度 $O(MN^2)$，其中 $M$ 是目标数量，$N$ 是种群规模。对于调度问题来说，目标数通常不超过 3，种群规模几百，这个复杂度完全可以接受。

3.2 拥挤度距离：替算法守住解的多样性

只做非支配排序还不够——如果算法只盯着前端最靠前的少量个体，种群很快会朝着 Pareto 前沿的某一小段集中，这就是所谓的早熟收敛。为了保持群体多样性，NSGA-II 引入了拥挤度距离这个概念。

拥挤度的大致计算思路是：对某一层前端中的每个个体，分别按每一个目标值排序，然后计算该个体与相邻两个个体在目标空间上的距离之和。距离越大，说明这个解在种群中的分布越"稀疏"，越值得保留。

你可以这么理解：同样是第 1 层的优秀候选，一个方案在 Pareto 前沿的左边，一个在中间，一个在右边，三个方案各自代表了不同的权衡取向。如果只保留其中一个，调度者就失去了选择余地。拥挤度排序的目的就是让前沿上的解尽可能铺开，保证最终交到使用者手里的不是"单一答案"，而是一整条可供挑选的曲线。

3.3 从父代到子代：锦标赛选择、SBX 交叉与多项式变异

NSGA-II 的进化流程和普通遗传算法没有本质区别，但它的选择算子是基于"非支配序号 + 拥挤度距离"进行的。常见做法是二元锦标赛选择：每次随机抽两个个体，优先比较非支配层级，层级更低的胜出；如果层级相同，则比较拥挤度距离，距离更大的胜出。这样做的目的是让胜出的个体既优秀又"身边空旷"。

交叉算子方面，实数编码的调度问题最适合用模拟二进制交叉（SBX）。SBX 的核心思想是让子代在父代附近产生，且两个子代关于父代中心对称分布；它的分布靠参数 $\eta_c$ 控制，$\eta_c$ 越大，子代越接近父代。变异算子则用多项式变异（Polynomial Mutation），它让个体以一定概率随机扰动某个基因位，扰动幅度由参数 $\eta_m$ 控制。这两个算子加在一起，既保证了算法在局部精细搜索，又保留了跳出局部最优的随机性。

3.4 精英保留策略为什么能防止退化

早期多目标遗传算法有个通病：每一代进化过程中，最好的解可能会因为交叉变异被破坏，造成"优秀基因丢失"，算法收敛过程来回震荡。NSGA-II 的解法是把父代种群和子代种群合并成一个两倍大小的候选池，然后在合并后的池子里重新做非支配排序和拥挤度排序，最后选出前 $N$ 个个体作为下一代种群。

这个精英保留机制带来一个非常实际的好处：无论子代多糟糕，父代中已有的最优非支配解一定有机会被保留到下一代。收到的效果就是算法每代的表现不会再倒退，收敛曲线十分稳定，你对运行结果也有了更大的确定性。这在跑调度工程时尤其重要——你不会希望夜间跑一个通宵优化，第二天早上看结果发现第 50 代比第 200 代还好，那说明代码写崩了。

4. Python 代码实现：从数据准备到算法主循环

4.1 模块划分与依赖安装

写算法前先把环境准备好。这篇文章的代码只依赖三个基础库：

pip install numpy pandas matplotlib

如果你的机器上还没装 Python，或者对 numpy 的安装方式不太熟，直接按常规做法：到官网装 Python 3.8 以上版本，然后在命令行执行上面的 pip 命令即可。pandas 主要用于读取光伏出力和来水的数据文件，matplotlib 用于最后把 Pareto 前沿画出来，算法核心只用到 numpy。

代码我按五个模块组织，结构如下：

• data_prep.py：生成或读取光伏出力序列、来水序列，打包成输入字典。
• problem.py：定义目标评估函数evaluate_individual和约束惩罚逻辑。
• nsga2.py：实现非支配排序、拥挤度计算、锦标赛选择、SBX 交叉、多项式变异。
• main.py：组装初始化、进化循环、精英保留。
• plot_results.py：可视化 Pareto 前沿和典型出力曲线。

这样做的好处是，后续如果要把单水电站改成梯级电站群，只需要改problem.py里的决策变量维度和水量平衡逻辑，nsga2.py一套进化机制完全不需要动。

4.2 光伏出力与来水数据准备

在我们这个例子里，光伏出力用一条典型日曲线按如下方式生成，方便复现；如果你有真实运行数据或预测系统输出，直接用pandas.read_csv读进来替换即可。

import numpy as np import pandas as pd np.random.seed(42) T = 24 # 调度时段数，每小时一个点 # 典型晴天的光伏出力曲线，6点到18点有出力 pv_power = np.zeros(T) for t in range(T): if 6 <= t <= 18: pv_power[t] = np.sin(np.pi * (t - 6) / 12) * 40 # 单位 MW，峰值 40 # 来水序列：可以是一条平稳径流 inflow = np.full(T, 12.0) # 单位 m3/s # 水电站与水库参数 V0 = 180.0 # 初始库容，万 m3 Vmin = 100.0 Vmax = 260.0 Ph_min = 5.0 # 水电最小出力 MW Ph_max = 50.0 # 水电最大出力 MW ramp_max = 12.0 # 最大爬坡速率 MW/h eta = 8.0 # 出力-发电流量转换系数

注意，eta的单位和数值取决于实际水电站的水头效率特性。我这里用一个简化常数表示"每立方米每秒发电流量能转换成多少兆瓦出力"，真正做工程时应该用库水位和发电流量的二维特性曲线查表。

4.3 目标评估函数

目标评估是整个算法的核心，也是所有约束发挥作用的地方。我在 Python 中定义了两个目标，并对约束违反量施加二次罚函数。这样保持了进化的连续性，梯度不光滑也无关紧要，因为 NSGA-II 不依赖梯度。

def evaluate_individual(P_h, pv_power, inflow, params): T = len(pv_power) V0 = params["V0"] Vmin = params["Vmin"] Vmax = params["Vmax"] eta = params["eta"] Ph_min = params["Ph_min"] Ph_max = params["Ph_max"] ramp_max = params["ramp_max"] # 通过水量平衡推算库容过程 V = np.zeros(T) spill = np.zeros(T) V[0] = V0 + inflow[0] - P_h[0] / eta if V[0] > Vmax: spill[0] = V[0] - Vmax V[0] = Vmax elif V[0] < Vmin: V[0] = Vmin # 简化处理：低于下限按修复到下限 for t in range(1, T): V[t] = V[t-1] + inflow[t] - P_h[t] / eta if V[t] > Vmax: spill[t] = V[t] - Vmax V[t] = Vmax elif V[t] < Vmin: V[t] = Vmin # 目标1：总发电量，取反方向 total_power = P_h + pv_power - spill * 0 # 弃水不发电 f1 = -np.sum(total_power) # 目标2：相邻时段出力波动 diff = np.diff(total_power) f2 = np.sum(diff ** 2) # 约束惩罚项 penalty = 0.0 # 出力上下限惩罚 penalty += np.sum((np.maximum(P_h - Ph_max, 0)) ** 2) penalty += np.sum((np.maximum(Ph_min - P_h, 0)) ** 2) # 爬坡约束惩罚 if T > 1: ramp_violation = np.maximum(np.abs(np.diff(P_h)) - ramp_max, 0) penalty += np.sum(ramp_violation ** 2) # 库容惩罚 penalty += np.sum((np.maximum(V - Vmax, 0)) ** 2) penalty += np.sum((np.maximum(Vmin - V, 0)) ** 2) penalty_weight = 1000.0 return f1 + penalty_weight * penalty, f2 + penalty_weight * penalty

代码里有几个处理细节值得解释。

第一，目标1取负号是因为我把 NSGA-II 统一写成小化问题；如果不取负，代码层面就需要单独区分最小化和最大化，容易出错。第二，水量平衡在每一时刻强行把库容钳制在上下限之间，超限部分记为弃水。这里我用的是一种"修复式"处理——库容超限后直接把水位拉回边界。实践表明它比单纯罚函数更容易让种群收敛到可行区域。第三，罚函数系数我取了 1000，一个"足够大到能淘汰明显违法个体"的经验值。具体数值可以根据目标量纲调整，但建议把它设为目标量级最大项的 10 到 100 倍。

4.4 非支配排序与拥挤度计算

nsga2.py里的非支配排序我采用最直观的两两比较实现。虽然理论上还可以用更高级的排序算法，但调度问题的种群规模在 200~300 左右，瓶颈不在排序，没必要过度优化。

def dominates(a, b): # 最小化问题：a 支配 b 当且仅当 a 全不劣且至少一个严格优于 b return all(x <= y for x, y in zip(a, b)) and any(x < y for x, y in zip(a, b)) def non_dominated_sort(fitness_values): N = len(fitness_values) S = [[] for _ in range(N)] n = [0] * N fronts = [[]] for i in range(N): for j in range(N): if i == j: continue if dominates(fitness_values[i], fitness_values[j]): S[i].append(j) elif dominates(fitness_values[j], fitness_values[i]): n[i] += 1 if n[i] == 0: fronts[0].append(i) k = 0 while fronts[k]: next_front = [] for i in fronts[k]: for j in S[i]: n[j] -= 1 if n[j] == 0: next_front.append(j) k += 1 if next_front: fronts.append(next_front) else: break return fronts

拥挤度距离的计算我按标准流程实现：每个目标方向单独处理，边界个体距离设为无穷大，中间个体累加归一化后的相邻距离。

def crowding_distance(front, fitness_values): distance = np.zeros(len(front)) num_obj = len(fitness_values[0]) for m in range(num_obj): front_sorted = sorted(front, key=lambda i: fitness_values[i][m]) distance[front_sorted[0]] = np.inf distance[front_sorted[-1]] = np.inf fmin = fitness_values[front_sorted[0]][m] fmax = fitness_values[front_sorted[-1]][m] if fmax == fmin: continue for k in range(1, len(front) - 1): distance[front_sorted[k]] += ( fitness_values[front_sorted[k+1]][m] - fitness_values[front_sorted[k-1]][m] ) / (fmax - fmin) return distance

实现中一个容易踩的坑是if fmax == fmin的分支。很多初版代码没处理同一个目标方向上所有个体值完全一样的情况，导致分母为 0，程序直接崩掉。我在写第一个版本时就在这个位置吃过亏，后来养成了习惯：凡是归一化计算，先检查分母。

4.5 遗传算子与精英选择

交叉和变异的实现直接沿用经典 SBX 与多项式变异。多项式变异的扰动公式里有一个随机指数项，实现时我习惯把均匀分布变量 u 压到极小值以上，避免对数爆炸：

def sbx_crossover(p1, p2, eta_c=15.0): child1 = np.empty_like(p1) child2 = np.empty_like(p2) rand = np.random.random(len(p1)) beta = np.where(rand <= 0.5, np.power(2.0 * rand, 1.0 / (eta_c + 1.0)), np.power(2.0 * (1.0 - rand), -1.0 / (eta_c + 1.0))) child1 = 0.5 * ((1 + beta) * p1 + (1 - beta) * p2) child2 = 0.5 * ((1 - beta) * p1 + (1 + beta) * p2) return np.clip(child1, 0, 1), np.clip(child2, 0, 1) def polynomial_mutation(child, eta_m=20.0, p_m=0.1): for i in range(len(child)): if np.random.random() < p_m: u = np.random.random() if u < 1e-6: u = 1e-6 if u <= 0.5: delta = np.power(2 * u, 1.0 / (eta_m + 1.0)) - 1.0 else: delta = 1.0 - np.power(2 * (1 - u), 1.0 / (eta_m + 1.0)) child[i] = np.clip(child[i] + delta, 0, 1) return child

这里有个工程细节：我先把个体统一编码到 [0,1] 区间，评估时才映射到真实的出力上下限。这么做的好处是交叉和变异算子不用关心每个基因位的物理量纲，也不用处理水电出力上下限不一样的问题；等到调用evaluate_individual时再做一个P_h = lb + child * (ub - lb)的映射就行。

精英选择逻辑如下：

def environmental_selection(parents, offspring, fitness_parents, fitness_offspring): combined = parents + offspring combined_fitness = fitness_parents + fitness_offspring pop_size = len(parents) fronts = non_dominated_sort(combined_fitness) new_pop = [] new_fitness = [] for front in fronts: if len(new_pop) + len(front) <= pop_size: for idx in front: new_pop.append(combined[idx]) new_fitness.append(combined_fitness[idx]) else: dist = crowding_distance(front, combined_fitness) order = sorted(front, key=lambda idx: -dist[idx]) remain = pop_size - len(new_pop) for idx in order[:remain]: new_pop.append(combined[idx]) new_fitness.append(combined_fitness[idx]) break return new_pop, new_fitness

4.6 主循环与 Pareto 前沿输出

主循环负责把上述函数串起来：初始化种群 → 评估 → 进化迭代 → 每代输出前端的解。我贴一个能直接跑通的最小版本：

def run_nsga2(params, pv_power, inflow, pop_size=100, max_gen=200, pc=0.9, pm=0.1): T = len(pv_power) # 初始种群 pop = [np.random.rand(T) for _ in range(pop_size)] for gen in range(max_gen): # 映射到物理变量并评估 fitness = [] for ind in pop: P_h = params["Ph_min"] + ind * (params["Ph_max"] - params["Ph_min"]) fitness.append(evaluate_individual(P_h, pv_power, inflow, params)) # 生成子代：锦标赛 + 交叉 + 变异 offspring = [] while len(offspring) < pop_size: # 锦标赛选择两个父代 idx1 = np.random.randint(pop_size) idx2 = np.random.randint(pop_size) parent1 = pop[idx1] if fitness[idx1][0] + fitness[idx1][1] <= fitness[idx2][0] + fitness[idx2][1] else pop[idx2] idx1 = np.random.randint(pop_size) idx2 = np.random.randint(pop_size) parent2 = pop[idx1] if fitness[idx1][0] + fitness[idx1][1] <= fitness[idx2][0] + fitness[idx2][1] else pop[idx2] if np.random.random() < pc: c1, c2 = sbx_crossover(parent1, parent2) else: c1, c2 = parent1.copy(), parent2.copy() c1 = polynomial_mutation(c1, p_m=pm) c2 = polynomial_mutation(c2, p_m=pm) offspring.extend([c1, c2]) offspring = offspring[:pop_size] # 精英保留 pop, fitness = environmental_selection(pop, offspring, fitness, [evaluate_individual( params["Ph_min"] + ind * (params["Ph_max"] - params["Ph_min"]), pv_power, inflow, params) for ind in offspring]) # 返回最终第1层 Pareto 前沿 fronts = non_dominated_sort(fitness) pareto_front_idx = fronts[0] pareto_pop = [pop[i] for i in pareto_front_idx] pareto_fitness = [fitness[i] for i in pareto_front_idx] return pareto_pop, pareto_fitness

这里需要说明，主循环里的锦标赛选择我做了简化处理，直接按两个目标的和来比较。更严格的做法应该是：先比较非支配层级，再比较拥挤度；但在这个小规模示例里，简化的目的是让代码更短、更容易理解。完整工程实现时请用non_dominated_sort+crowding_distance组装锦标赛选择算子，效果会更好。

5. 结果解读：如何从 Pareto 前沿挑选落地调度方案

5.1 典型 Pareto 前沿长什么样

跑完 200 代以后，代码会返回一组非支配解。把它们的目标函数值画在二维平面上，就会出现一条从左上到右下的弯曲线。横轴我习惯取"发电量取反后的绝对值"（即总发电量），纵轴取波动指标。左下角代表高发电量低波动，右上角代表低发电量高波动，两个目标互为牺牲。

我在项目里跑出来的前沿大致呈 L 形下降。在发电量偏小的区域，曲线比较陡峭：稍微增加一点出力波动，就能换来明显的发电量提升。在发电量偏大的区域，曲线开始平缓：这时候再想增加发电量，需要付出的波动代价非常大，性价比很低。

这个形状本身就是调度信息的体现。它告诉你：你的水光互补系统最"划算"的运行区间在哪些地方，以及系统的瓶颈在哪里。如果前沿在某一端突然截断，说明你触碰到了约束边界，比如水库库容不够了，或者光伏装机容量限制了总出力。

5.2 肘点选择和折中策略

面对一整条 Pareto 前沿，工程上最常用的选点方法是"肘点法"——优先选择曲率最大、两个目标权衡收益最均衡的那个点。计算方式是把前沿两端连接成一条直线，然后找到距离这条直线最远的点。

def select_elbow_point(pareto_fitness): f1 = np.array([fit[0] for fit in pareto_fitness]) f2 = np.array([fit[1] for fit in pareto_fitness]) # 归一化 f1_norm = (f1 - f1.min()) / (f1.max() - f1.min()) f2_norm = (f2 - f2.min()) / (f2.max() - f2.min()) idx_start = 0 idx_end = len(f1_norm) - 1 A = f1_norm[idx_end] - f1_norm[idx_start] B = f2_norm[idx_end] - f2_norm[idx_start] distances = np.abs(A * f1_norm - B * f2_norm + f2_norm[idx_start] * f1_norm[idx_end] - f2_norm[idx_end] * f1_norm[idx_start]) / np.sqrt(A**2 + B**2) elbow_idx = np.argmax(distances) return elbow_idx

肘点选择的哲学是"不做极端偏好判断"：它默认决策者没有非常明确的倾向，只是想找一个各方面都能接受的平衡点。实际工程中，调度员往往还有自己的偏好——比如午间负荷高峰期间希望出力更高，比如外部输电通道对波动率有硬性指标。这时就应该把这些偏好转成约束或者权重，去 Pareto 前沿上重新筛选。比如电网要求波动指标低于某阈值，你只需要把前沿上所有满足条件的非支配解挑出来，再在其中选发电量最高的。

5.3 用超体积和分布间距评价算法质量

选完方案，还有一个问题经常被忽略：怎么判断这次算法运行的质量好不好？我在调参时主要看两个指标。

第一个是超体积指标（Hypervolume，HV），它度量 Pareto 前沿与一个参考点围成的目标空间体积。HV 越大，说明这一组解在整体上距离理想界越近、覆盖面越广。第二个是分布间距指标（Spacing），它统计相邻非支配解之间的欧氏距离波动情况。Spacing 越小，说明前沿上的点分布越均匀。

这些指标不需要每次跑都算，但在你调整种群大小、迭代代数、交叉变异概率时非常有用。我通常的做法是：固定随机种子，分别用不同参数组合跑 10 次，比较 HV 均值和方差，选择均值高且方差小的参数。你可以把 hv 的计算简化理解为在目标空间里填充小立方体，不必引入复杂的外部库；如果追求严谨，可以装pymoo这类现成库来算，不过自己写一遍几十行代码也够用。

6. 调参经验与常见的坑

6.1 约束惩罚容易踩的两个雷

第一个雷是罚函数权重过小。如果权重设置得太温和，算法会认为违反库容约束也不算什么大不了的事，最终的"最优解"可能根本不在可行域里。我建议先用一个很大的权重把种群压入可行域，再考虑目标函数的精细区分。判断方法很简单：画一张图，把每一代种群中违反约束的个体数量可视化，如果到 100 代以后还有大量个体整天踩线，说明罚函数力度不够。

第二个雷是库容修复逻辑写错顺序。在我列出的evaluate_individual里，如果库容先超上限就强制截断，但这时弃水没有正确地从水量平衡中体现，后续时段的库容过程就会失真。正确顺序必须是：先计算天然来水和发电用水后的库容，再判断是否超限，超限部分记为弃水，最后让库容回到界内。顺序反了，跑出来的调度方案可能让你白白扔掉大量可发的水量。

6.2 种群大小与迭代代数怎么定

这个问题没有标准答案，但我可以给一个经过试验的经验区间。对于决策变量维度为 24、目标数为 2 的水光互补问题，种群规模取 100 到 150，迭代代数取 200 代左右，通常能得到质量不错的 Pareto 前沿。如果调度周期改成 96 时段，决策变量维度变成 96，搜索空间指数级增大，种群规模建议提升到 200 以上，迭代代数增加到 300 到 500。

判断是否收敛有一个简单方法：每隔 20 代记录一次当前第 1 层前端的 HV 值，画一条 HV 随代数变化的曲线。当曲线进入平台期、后续涨幅小于 1% 时，说明算法已经稳定。我实测下来，24 时段问题通常在 100 代以后进入平台期，但为了保险，一般多跑几十代。

6.3 交叉变异参数与随机种子

SBX 的 $\eta_c$ 我习惯取 15，多项式变异的 $\eta_m$ 取 20，变异概率 $p_m$ 取 0.1。这个组合在多个调度实例里表现稳定，但不是唯一正确解。需要提醒的是，变异概率不要设太高——调度问题中决策变量之间存在强的时间连续性，某个时段出力的突变往往造成水量平衡连锁反应，过高变异率会让种群在后期仍然剧烈震荡，Pareto 前沿难以收敛。想追求更稳的收敛，可以尝试让变异概率在迭代后期动态衰减，从 0.15 线性降到 0.03。

随机种子也是工程中必须注意的。我跑算法时总是固定一个种子数组做多组对比，否则你根本无法判断某一组参数效果好是参数本身的功劳还是随机运气。推荐的做法是：同一组参数下跑 5 个不同种子，取 HV 均值和标准差来评价。

6.4 光照数据颗粒度与调度场景选择

最后想多说一句关于光伏数据的事。很多初版实现喜欢直接用"最大出力"这样的静态数据，但实际调度中光伏预测出力是随时段变化的，而且不同天气类型差异极大。典型晴天、多云天、阴雨天生成三条场景曲线，分别跑一遍 NSGA-II，你会发现 Pareto 前沿的形状完全不同：晴天时水电几乎不需要承担午间出力，前沿整体偏向左下角；阴雨天时水电几乎全程兜底，两条目标曲线都变差。

这就是为什么在实际工程里，多目标优化的价值不只是"算出一条曲线"，而是"面对不同光伏场景都能提前给出可行的调度预案"。调试阶段我用典型的晴天数据，上线前一定要替换成多场景数据做鲁棒性验证，否则调度方案经不起天气变化的考验。

整条代码跑下来以后，我自己最大的体会是：NSGA-II 并不神秘，它只是把"怎么权衡"这个问题从中选取一个答案，变成了生成一大批候选答案；真正考验功力的地方，反而在模型的目标函数设计、约束处理和结果解读。把这三件事想明白，代码层面的实现就是水到渠成的事情。希望这套方法能帮你少走一些弯路，尤其是水量平衡和罚函数那部分，值得多花时间验证。