做水电调度的朋友应该都有体会,光伏大规模并网之后,原来那套"按经验排计划"的方式越来越不好使了。光伏出力跟着天气走,早上一波爬坡、傍晚一波陡降,中间还可能被云层切出好几个"凹坑"。水电虽然响应快,但也不是想怎么调就怎么调——库容有限、来水不确定、下游还有生态流量要求。在这种两头都受约束的情况下,怎么让水电机组和光伏电站协同出力,既保证发电收益、又让并网功率尽量平滑、还能提高新能源消纳率,就成了一个典型的多目标优化问题。
我最早接触这个课题是在一个小型梯级水电站群,装机不大但光伏配了不少。最开始我们试过人工经验调度,也试过线性加权把多个目标揉成一个单目标去优化,效果都不理想。后来换了基于非支配排序遗传算法(NSGA-II)的多目标优化调度方案,一次性生成一组Pareto最优调度方案,让运行人员根据当天的实际情况挑方案,这才算把问题解决得比较舒服。这篇文章就把我完整跑通的思路和Python代码实现分享出来,内容包括问题建模、NSGA-II算法原理、手写算法核心代码、结果解读和工程调参经验,适合正在做新能源调度、电力系统优化或者刚接触多目标优化的朋友参考。
1. 水光互补为什么值得做:先看清调度问题的底层矛盾
1.1 光伏出力天生"不听话"
光伏出力的本质是跟光照强度走的,而光照强度受昼夜、云层、季节、纬度多方面影响,表现在出力曲线上就是:夜间出力为零、白天出现一个峰值、阴雨天剧烈波动。下面这组典型日光伏出力数据(装机30MW,归一化后)很能说明问题:
| 时段 | 出力/MW | 时段 | 出力/MW |
|---|---|---|---|
| 00:00-04:00 | 0 | 12:00-13:00 | 28.5 |
| 05:00-06:00 | 0-3 | 14:00-15:00 | 22.5 |
| 07:00-08:00 | 9-18 | 16:00-17:00 | 15 |
| 09:00-10:00 | 22.5-27 | 18:00-19:00 | 6-3 |
| 11:00-12:00 | 28.5-30 | 20:00-24:00 | 0 |
注意17:00到18:00这个时段,出力会从十几兆瓦快速掉到接近零,如果水电不能及时顶上,整个并网功率就会出现一个巨大的下坡,对电网频率支撑非常不友好。更让人头疼的是,这个曲线并不是固定不变的,一片云飘过来,10分钟内出力可能掉一半。所以说光伏出力天生"不听话",这是水光互补调度要解决的第一层矛盾。
1.2 水电是理想的补偿电源,但也不是无限可调
水电之所以被当作光伏的天然搭档,核心原因有两个:一是响应速度快,水电机组从接到指令到完成出力调整通常只需要1-5分钟,比火电快一个量级;二是调节范围大,在满足最小技术出力的前提下,水电可以大幅调整负荷。但"理想"这个词得打上引号,水电的调节能力是受到硬约束的:
- 库容约束:水库能存的水有限,多发意味着后半夜或者明天可能没水可发。
- 来水约束:天然来水不确定,枯水期和丰水期差别非常大。
- 生态流量约束:下游不能断流,最小下泄流量是死约束。
- 震动区约束:水轮机组在部分负荷区运行会有振动,实际可调区间不是连续的。
这就意味着,水电不能简单地"光伏多我就少发,光伏少我就多发",必须统筹考虑一整天的水量平衡。你看,第二层矛盾也出来了:水电想补偿光伏,但自己有"体力上限"。
1.3 互补调度到底在"补"什么
从运行角度看,水光互补调度的核心诉求有三个:时间互补、功率互补、电量互补。
时间互补好理解,白天光伏出力大,水电可以压低甚至停机蓄水;晚上光伏归零,水电再顶上来保证供电;凌晨和傍晚这种光伏快速变化的时段,水电起爬坡补偿作用。功率互补是指联合出力的总功率曲线不能大起大落,爬坡率要控制在合理范围之内。电量互补则是从更长时间尺度(周、月、年)上考虑,不能让水电因为补偿光伏而过度放水,导致后续时段无水可用。
说白了,水光互补调度就是要在"水电怎么用"上做文章:每个时段水电该发多少?是在白天多存水还是多放水?这些都是决策变量。而判断一个方案好不好,不能只看收益,还要看电网稳定性和新能源利用率,这就自然引出了多目标优化的需求。
2. 多目标优化与NSGA-II:为什么不能只追求发电量最大
2.1 多目标之间的"按起葫芦浮起瓢"
初学者最容易犯的错误是把多目标调度的所有目标用权重揉成一个数,然后用单一目标优化去解。比如"总收益最大"加上"波动最小"就变成 min(收益 + 100×波动)。这个做法不是不行,但有一个致命问题:你凭什么定这个100?不同目标的量纲不同,权重本质上反映了决策者的偏好,但实际运行中这个偏好是时变的——丰水期更在乎消纳,枯水期更在乎收益,极端天气更在乎平稳。靠人工拍脑袋定权重,很难保证每次结果都合理。
更现实的问题是,多个目标之间往往是冲突的。拿水光互补调度来说:
- 想发电收益最大,最好让水电在电价高的晚高峰和早高峰开足马力;
- 想并网出力最平稳,最好让水电出力曲线跟光伏出力曲线形成镜像补偿,全天保持一条平线;
- 想让光伏消纳率最高,水电就得随时给光伏"让路",光伏一出力水电立刻压低。
这三个愿望不可能同时完全满足。比如电价高峰恰好是光伏出力不足的傍晚,这时水电开足马力去追收益,联合出力就会从低谷猛冲上去,波动自然变大。你选择了收益最大化,就必然牺牲一部分平稳性;你选择了绝对平稳,就必然损失一部分高峰期的发电收入。
这就是多目标优化的核心特征:不存在一个解能让所有目标同时达到最优,只存在一组互不支配的折中解,需要调度员根据实际场景从中挑选。
2.2 Pareto最优解的直觉理解
讲清楚"互不支配",一个二维例子就够了。假设现在只有两个目标——收益最大化(目标1)和并网波动最小化(目标2),有两个候选调度方案A和B:
- A方案:收益100万元,波动指标50
- B方案:收益95万元,波动指标40
那么A和B谁更好?说不清。A收益高但波动大,B波动小但收益低,两个方案在不同的维度上各自领先,这时候我们就说A和B互不支配,它们都属于Pareto最优解。但如果还有一个方案C:收益90万元,波动指标60,那C就被A支配了——A在收益和波动两个指标上都比C好。
把所有互不支配的解放在一起,就形成了Pareto前沿。NSGA-II这类多目标进化算法做的事情,就是通过不断进化,让解集尽可能逼近真实的Pareto前沿,并且在整条前沿上分布得均匀分散——既要有高收益低平滑度的方案,也要有低收益高平滑度的方案,还要有中间状态。
2.3 NSGA-II的三个核心机制
NSGA-II能在众多多目标优化算法里成为最常用的基线方法,靠的是三个精心设计的机制:
快速非支配排序。把所有个体按照支配关系分成一层一层的"前沿等级":第一层是当前种群中所有不被任何其他个体支配的个体;去掉第一层之后,剩下的个体里再找出互不支配的,作为第二层;依此类推。这个排序结果决定了每个个体的"优劣等级",等级越低(越靠前)说明这个解越优秀。
拥挤度距离。同一层内,个体之间的优劣怎么比较?答案是比较拥挤度——个体在目标空间里周围其他个体密集的程度。拥挤度大的个体说明它周围比较空旷,保留它有利于维持解的多样性;拥挤度小的个体周围挤满了同类,淘汰它也不会让前沿丢失太多信息。这样算法在保留前沿多样性的同时,也能收敛到比较完整的Pareto形状。
精英保留策略。每次进化时,把父代和子代合并成一个大的种群,统一做非支配排序和拥挤度排序,然后从好的开始依次挑选出下一代(种群大小不变)。父代中的优秀解不会因为一次交叉变异就丢失,保证了算法的收敛性不会退化。
这三个机制环环相扣,保证了NSGA-II既能收敛到前沿,又能在前沿上分布得均匀。
2.4 为什么选NSGA-II而不是别的算法
有人会问,现在MOEA/D、NSGA-III、SPEA2这些算法也很多,为什么我推荐NSGA-II?我的理由是:NSGA-II实现简单、参数少、稳定性好,是理解多目标进化算法最好的切入点。MOEA/D需要你预先设计权重向量分布,NSGA-III在高维问题上优势明显但代码复杂度高,SPEA2需要维护额外的档案集合。对于水光互补调度这种目标数量在2-4个、实时性要求又不像日内滚动那么极端的问题,NSGA-II的性价比是最高的。而且它的核心算子(锦标赛选择、模拟二进制交叉、多项式变异)理解了之后,往其他算法迁移非常容易。
3. 调度问题的数学模型:目标函数、约束与编码
3.1 目标函数怎么定
建立模型时我采用一个典型的日调度场景:调度周期为24小时,时间间隔取1小时,水电站装机50MW,光伏电站装机30MW,需要为水电机组确定全天24个时段的出力计划。我选了三个目标,分别对应经济效益、电网安全和新能源消纳:
目标1:发电收益最大化
$$\max f_1 = \sum_{t=1}^{24} \left( P^{H}_t + P^{PV}_t \right) \times c_t$$
其中 $P^{H}_t$ 是水电在t时段出力,$P^{PV}_t$ 是光伏出力(由预测曲线给定),$c_t$ 是分时电价。考虑午间光伏大发时段电价相对较低,而晚高峰电价较高,这个目标会驱动水电向高峰时段集中。
目标2:并网出力波动最小化
$$\min f_2 = \sum_{t=1}^{23} \left| \left( P^{H}{t+1} + P^{PV}{t+1} \right) - \left( P^{H}_t + P^{PV}_t \right) \right|$$
这个目标直接衡量联合出力曲线的总爬坡量,数值越小,说明并网功率越平稳,对电网频率和电压的冲击就越小。因为光伏曲线是固定的,水电曲线就成了决定波动大小的关键。
目标3:光伏消纳率最大化(弃光量最小化)
$$\max f_3 = \sum_{t=1}^{24} \left( P^{PV}{t} - P^{curtail}{t} \right) / \sum_{t=1}^{24} P^{PV}_{t}$$
严格来说这个目标需要引入弃光变量并联合约束求解。在简化版模型里,我假设水光联合出力在满足并网上下限的前提下以全额消纳为优先,弃光只在联合出力超过并网上限时发生。这样模型既能反映实际物理过程,又不至于太复杂。
3.2 约束条件有哪些
约束是调度模型的骨架,没有约束的优化结果在现实中根本无法执行。我至少考虑了以下几类:
- 水电出力上下限:$P^{H}{min} \le P^{H}t \le P^{H}{max}$,本文取 $P^{H}{min}=5$MW,$P^{H}_{max}=50$MW。这里要注意,水电机组不能在低于最小技术出力的状态下稳定运行。
- 水电爬坡约束:$|P^{H}_{t+1} - P^{H}_t| \le R_H$,取 $R_H = 15$MW/h。水电站虽然有快速调节能力,但上下游水位变动、引水系统水击等问题要求出力变化速率必须受限。
- 日发电量/水量平衡约束:$\sum_{t=1}^{24} P^{H}t \cdot \Delta t = E^{H}{target}$,$\Delta t=1$h。这相当于说一天之内水库可用的水量对应一个固定的发电量。如果水库在白天过度放水,后半夜就没有足够的水量维持晚高峰出力。
- 联合并网功率上限:$P^{H}t + P^{PV}t \le P^{grid}{max}$,取 $P^{grid}{max} = 75$MW。当光伏大发而水电又不能瞬间压低时,就有弃光风险。
在NSGA-II的框架里,这些约束我统一用罚函数处理:对越界的个体,在目标函数里加上惩罚项,让它在非支配排序中处于劣势。罚函数系数是一个需要谨慎调节的变量——太小约束失效,太大又压缩了搜索空间,后面第6章我会专门讲。
3.3 决策变量与编码
决策变量就是水电机组24个时段的出力向量 $\mathbf{x} = [P^{H}_1, P^{H}2, \ldots, P^{H}{24}]$,实数编码,维度24。为什么不用二进制编码?因为调度变量是连续量,二进制编码在解码、交叉、变异时都要做额外转换,而且精度受编码长度限制。实数编码配合模拟二进制交叉(SBX)和多项式变异,可以保证变量在连续空间里平滑演化,Gen操作后得到的新个体天然在取值范围内,只需要再做一次越界钳位就行。
每个调度方案在算法里就是一条"染色体"(长度为24的浮点数组),种群则由100条染色体构成。目标函数接收到一条染色体后,结合光伏预测曲线、电价曲线和负荷信息,计算出对应的三个目标值。这个过程完全解耦,目标函数内部不需要知道任何遗传操作细节,方便以后替换更精细的水电模型。
4. Python完整实现:从数据准备到Pareto前沿输出
4.1 环境准备与基础数据构造
实现只用numpy和matplotlib两个库,不需要额外装DEAP或pymoo,这样便于读者理解算法内部机制。安装环境时建议用Python 3.8以上版本,pip安装numpy和matplotlib即可:
pip install numpy matplotlib基础数据包括光伏预测出力、分时电价、水电参数这三部分。为了让结果可复现,我给随机种子设了一个固定的值。下面是数据构造代码:
import numpy as np import matplotlib.pyplot as plt # 固定随机种子,保证结果可复现 np.random.seed(42) T = 24 # 调度时段数,1h一个时段 # 光伏预测出力曲线(30MW装机,归一化后乘装机容量) pv_norm = np.array([ 0.00, 0.00, 0.00, 0.00, 0.03, 0.15, 0.35, 0.60, 0.85, 0.95, 1.00, 0.92, 0.75, 0.58, 0.40, 0.22, 0.10, 0.02, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00 ]) PV_CAP = 30.0 # MW pv_power = pv_norm * PV_CAP # 分时电价(元/kWh),平段、峰段、谷段 price = np.array([ 0.38, 0.35, 0.32, 0.32, 0.35, 0.50, 0.75, 0.75, 0.70, 0.55, 0.42, 0.42, 0.42, 0.42, 0.55, 0.70, 0.75, 0.72, 0.60, 0.50, 0.42, 0.38, 0.35, 0.35 ]) # 水电参数 H_MIN = 5.0 # MW 最小技术出力 H_MAX = 50.0 # MW 最大出力 RAMP = 15.0 # MW/h 爬坡限制 E_TARGET = 600.0 # MWh 日发电量目标 GRID_MAX = 75.0 # MW 并网功率上限 # 水电日发电量对应水量平衡:所有时段出力之和应接近E_TARGET这里面E_TARGET=600MWh是我人为给定的日发电量目标。实际工程中这个值来自水库调度图或中长期计划的分解结果,代表当天"能用的水"折算出的电量总额。
4.2 目标函数与约束判定代码
目标函数接收一条24维的调度曲线,返回三个目标值。注意:约定所有目标统一为最小化,所以收益和消纳率都要加负号。
def evaluate(x, pv_power=pv_power, price=price): """ 计算一个调度方案的目标函数值 x: 水电24时段出力向量 [P_H_1, ..., P_H_24] (MW) 返回: [亏损(负收益), 总爬坡量, 弃光惩罚] """ total = x + pv_power # 联合出力 # 目标1:发电收益最大化 -> 最小化负收益 revenue = np.sum(total * price) # 单位:MW * 元/kWh = 千元 f1 = -revenue # 目标2:并网出力波动最小化 f2 = np.sum(np.abs(np.diff(total))) # 目标3:弃光量最小化 # 联合出力超过并网上限的部分视为弃光(简化处理) curtail = np.maximum(total - GRID_MAX, 0.0) f3 = np.sum(curtail) return np.array([f1, f2, f3]) def constraint_violation(x): """ 返回约束越界的总惩罚量,用于罚函数 包括:水电上下限、爬坡约束、日发电量偏差 """ viol = 0.0 # 上下限越界 viol += np.sum(np.maximum(x - H_MAX, 0.0)) + np.sum(np.maximum(H_MIN - x, 0.0)) # 爬坡越界 ramp_diff = np.abs(np.diff(x)) - RAMP viol += np.sum(np.maximum(ramp_diff, 0.0)) # 日发电量偏差 viol += 0.1 * abs(np.sum(x) - E_TARGET) # 系数0.1缩放 return viol注意constraint_violation返回的是标量,后续会以罚函数形式叠加到目标值上。把"日发电量偏差"单独算而不是直接约束每个时段,是为了保留水电调度的灵活性——只要一天总水量不超,具体各时段怎么分配由算法自己去权衡。
4.3 NSGA-II核心算子的Python实现
首先是快速非支配排序。为了效率,我按种群中每个个体与其他个体的支配关系来构建分层,时间复杂度是O(MN²),M是目标数,N是种群大小。对本文这种100×3的小规模问题完全够用:
def dominates(a, b): """判断a是否支配b:a在所有目标上不差于b,且至少一个目标严格更优""" return np.all(a <= b) and np.any(a < b) def fast_non_dominated_sort(values): """ 快速非支配排序,返回前沿分层列表 values: (pop_size, n_obj) 目标值矩阵,所有目标越小越好 """ n = values.shape[0] dominated_count = np.zeros(n) dominate_list = [[] for _ in range(n)] fronts = [[]] for i in range(n): for j in range(n): if i == j: continue if dominates(values[i], values[j]): dominate_list[i].append(j) elif dominates(values[j], values[i]): dominated_count[i] += 1 if dominated_count[i] == 0: fronts[0].append(i) k = 0 while fronts[k]: next_front = [] for p in fronts[k]: for q in dominate_list[p]: dominated_count[q] -= 1 if dominated_count[q] == 0: next_front.append(q) k += 1 fronts.append(next_front) return fronts[:-1]然后是拥挤度距离计算。核心思想是在每个目标维度上,把一个前沿内的个体按目标值排序,然后计算每个个体前后相邻个体之间的归一化距离之和,边界个体直接赋无穷大:
def crowding_distance(values, front): """ 计算某个前沿内所有个体的拥挤度距离 values: (pop_size, n_obj) 目标值矩阵 front: 当前前沿的个体索引列表 """ if len(front) <= 2: return {p: float('inf') for p in front} dist = {p: 0.0 for p in front} n_obj = values.shape[1] for obj in range(n_obj): front_sorted = sorted(front, key=lambda p: values[p, obj]) obj_min = values[front_sorted[0], obj] obj_max = values[front_sorted[-1], obj] norm = obj_max - obj_min if norm < 1e-12: continue dist[front_sorted[0]] = float('inf') dist[front_sorted[-1]] = float('inf') for k in range(1, len(front_sorted) - 1): dist[front_sorted[k]] += (values[front_sorted[k+1], obj] - values[front_sorted[k-1], obj]) / norm return dist锦标赛选择、模拟二进制交叉(SBX)和多项式变异是遗传操作的三板斧。锦标赛选择的思路是随机抽两个个体,优先选非支配等级低的;等级相同选拥挤度大的。SBX交叉模拟二进制编码交叉的分布特性,让后代有较大概率落在两个父代之间附近,也有一定概率跳得更远:
def tournament_selection(pop, ranks, dists, k=2): """锦标赛选择:随机取k个个体,返回最优的一个""" idx = np.random.choice(len(pop), k, replace=False) best = idx[0] for i in idx[1:]: if ranks[i] < ranks[best]: best = i elif ranks[i] == ranks[best] and dists[i] > dists[best]: best = i return pop[best].copy() def sbx_crossover(p1, p2, eta_c=20): """模拟二进制交叉,输入两个父代向量,返回两个子代向量""" u = np.random.rand(len(p1)) beta = np.empty_like(u) idx = u <= 0.5 beta[idx] = (2 * u[idx]) ** (1 / (eta_c + 1)) beta[~idx] = (1 / (2 * (1 - u[~idx]))) ** (1 / (eta_c + 1)) c1 = 0.5 * ((1 + beta) * p1 + (1 - beta) * p2) c2 = 0.5 * ((1 - beta) * p1 + (1 + beta) * p2) return c1, c2 def polynomial_mutation(child, eta_m=20): """多项式变异,变异后做边界钳位""" u = np.random.rand(len(child)) delta = np.empty_like(u) idx = u < 0.5 delta[idx] = (2 * u[idx]) ** (1 / (eta_m + 1)) - 1 delta[~idx] = 1 - (2 * (1 - u[~idx])) ** (1 / (eta_m + 1)) child = child + delta * (H_MAX - H_MIN) return np.clip(child, H_MIN, H_MAX)4.4 主循环与可视化输出
主循环是整个算法的心脏:生成父代 -> 产生子代 -> 合并排序 -> 挑出下一代。每一代的流程完全一致,循环MAX_GEN次后返回最终种群:
POP_SIZE = 100 MAX_GEN = 200 CX_PROB = 0.9 MUT_PROB = 1.0 / 24 PENALTY_COEF = 10.0 # 约束惩罚系数 # 初始化种群:均匀随机生成水电出力 pop = np.random.uniform(H_MIN, H_MAX, (POP_SIZE, T)) for gen in range(MAX_GEN): # 1. 评价所有个体目标值(含罚函数) vals = np.zeros((len(pop), 3)) for i, ind in enumerate(pop): f = evaluate(ind) viol = constraint_violation(ind) vals[i] = f + PENALTY_COEF * viol # 2. 非支配排序 + 拥挤度距离 fronts = fast_non_dominated_sort(vals) ranks = np.zeros(len(pop), dtype=int) dists = np.zeros(len(pop)) for k, front in enumerate(fronts): for p in front: ranks[p] = k d = crowding_distance(vals, front) for p in front: dists[p] = d[p] # 3. 锦标赛选择父代 parents = [] for _ in range(POP_SIZE): parents.append(tournament_selection(pop, ranks, dists)) # 4. 交叉 + 变异生成子代 offspring = [] for i in range(0, POP_SIZE, 2): p1, p2 = parents[i], parents[i+1] if np.random.rand() < CX_PROB: c1, c2 = sbx_crossover(p1, p2) else: c1, c2 = p1.copy(), p2.copy() if np.random.rand() < MUT_PROB: c1 = polynomial_mutation(c1) if np.random.rand() < MUT_PROB: c2 = polynomial_mutation(c2) offspring.append(c1) offspring.append(c2) # 5. 父代+子代合并,精英选择 combined = np.vstack([pop, np.array(offspring[:POP_SIZE])]) combined_vals = np.zeros((len(combined), 3)) for i, ind in enumerate(combined): f = evaluate(ind) viol = constraint_violation(ind) combined_vals[i] = f + PENALTY_COEF * viol fronts = fast_non_dominated_sort(combined_vals) combined_ranks = np.zeros(len(combined), dtype=int) combined_dists = np.zeros(len(combined)) for k, front in enumerate(fronts): for p in front: combined_ranks[p] = k d = crowding_distance(combined_vals, front) for p in front: combined_dists[p] = d[p] # 按 (rank, -distance) 排序后取前POP_SIZE个 sort_key = [(combined_ranks[i], -combined_dists[i]) for i in range(len(combined))] order = np.argsort(np.array(sort_key, dtype=[('r', int), ('d', float)])) pop = combined[order[:POP_SIZE]] vals = combined_vals[order[:POP_SIZE]] # 最终结果 final_vals = vals final_front = [i for i, r in enumerate(combined_ranks[order[:POP_SIZE]]) if r == 0]可视化部分比较简单,三维目标空间直接画散点图,也可以两两配对画二维投影。我一般同时画两幅图:一幅是种群整体的分布,一幅是最终Pareto前沿的曲线:
# 绘制Pareto前沿:三维目标空间 fig = plt.figure(figsize=(10, 7)) ax = fig.add_subplot(111, projection='3d') ax.scatter(final_vals[final_front, 0], final_vals[final_front, 1], final_vals[final_front, 2], c='crimson', s=30, alpha=0.8) ax.set_xlabel('负收益 (越小越好)') ax.set_ylabel('总爬坡量 (越小越好)') ax.set_zlabel('弃光量 (越小越好)') plt.tight_layout() plt.savefig('pareto_front_3d.png', dpi=150) plt.show()实际上三维图对浏览者不太友好,我更常用的是三张二维投影子图,分别在收益-波动、收益-弃光、波动-弃光三个平面上看Pareto分布。
5. 结果解读与算法参数调优
5.1 Pareto前沿怎么看
跑完200代,种群会稳定在一个三维的Pareto前沿面上。我第一次跑出来的时候其实有点懵,三维散点图虽然好看,但脱离调度场景很难直接指导运行。后来我的习惯是:先在二维投影里找规律,再回到调度曲线里验证物理含义。
举个例子,在"负收益-总爬坡量"二维投影里,前沿通常是一条单调递减的曲线——横轴负收益越小(收益越高),纵轴爬坡量就越大。这说明"多赚钱"和"平稳出力"确实存在直接冲突。在实际项目里,运行人员会先画一条水平线表示"今天允许的最大爬坡量",然后去看这条线跟Pareto前沿的交点对应的收益值,再从交点附近挑一个调度方案。如果当天电网运行比较稳定、对爬坡不太敏感,就可以往收益大的方向挑;如果遇到特殊运行方式,就把平稳性放在更优先的位置。
另外,从Pareto前沿上还可以反推出水电的调度行为规律。我观察过一个有趣的现象:前沿左端(重收益)的方案,水电出力曲线几乎全部集中在18:00-21:00晚高峰,白天的光伏时段压到最小技术出力甚至接近停机;前沿右端(重平稳)的方案,水电曲线跟光伏曲线形成近乎镜像的反向关系,光伏升水电降、光伏降水电阻。这两种极端方案对应的物理含义都说得通,验证了模型没有跑偏。
5.2 关键参数如何设置
NSGA-II参数不多,但每一个都直接影响结果质量。我把自己实测下来比较可靠的参数范围整理成了表格:
| 参数 | 推荐范围 | 说明 |
|---|---|---|
| 种群大小 | 50-200 | 本文取100。太小前沿不完整,太大单代计算慢 |
| 最大迭代次数 | 100-500 | 本文取200,已经能收敛,500代以后基本无变化 |
| 交叉概率 | 0.8-0.95 | 本文取0.9,交叉是产生新解的主要手段 |
| 变异概率 | 1/维度-1/5 | 本文取1/24,相当于平均每条染色体有一个基因变异 |
| SBX分布指数ηc | 15-30 | 越大后代越接近父代,本文取20 |
| 多项式变异ηm | 15-30 | 越大变异幅度越小,本文取20 |
| 罚函数系数 | 1-100 | 需要调试,太小会越界,太大会丢失搜索方向 |
一个很重要的经验:种群大小比迭代次数更关键。我开始用50个个体跑300代,Pareto前沿总是有缺口;后来改成100个个体跑200代,前沿一下子就完整了。原因是种群多样性不足时,即使给再多代,种群也很难跳出局部区域,交叉算子没有足够的基因素材可用。
5.3 与单目标优化的对比
为了验证多目标优化的价值,我把同一个调度问题用线性加权法转成单目标也跑了一遍。加权系数取"收益权重0.5、爬坡权重0.3、弃光权重0.2"(预先归一化),然后用粒子群算法搜索最优解。结果很有意思:加权法找到的那个"最优解",放到NSGA-II的Pareto前沿上一查,发现它被一条前沿上的解支配——收益比它高、爬坡比它小、弃光也不比它多。
这个现象说明了多目标进化算法的一个核心优势:它不需要事先知道权重,就能通过种群进化的方式找到一条覆盖整个偏好空间的解集。你事后从Pareto前沿上挑出来的解,往往比你事前拍脑袋定权重得到的最优解更接近真实最优。当然,如果决策者目标非常明确、偏好长期固定不变,加权法也有它的用武之地,毕竟计算量小、实现简单。但对调度决策支持这类场景,Pareto解集的"一次计算、反复挑选"特性,明显更适合实际运行。
6. 工程实现中的经验教训:我踩过的坑
6.1 约束处理不能只靠罚函数
上一章的代码里我用的是一个固定惩罚系数,实际用的时候踩了坑。罚函数系数太小,种群里的个体会大量越界,非支配排序里"越界的烂解"会占据前沿位置;罚函数系数太大,又等于给优化目标加了一个大的常偏置,个体会过度"躲避"约束边界,导致最优解永远落在可行域的内部而非边界上。而真正的Pareto最优往往就骑在几个等式约束边界上(比如日发电量刚好等于E_TARGET)。
后来我的解决方法是动态罚函数:开始时惩罚系数给小一点(比如1.0),让算法能在可行域和不可行域交界处自由探索;随着代数增加线性增大到10-20,逐步把种群压进可行域。这个思路有点像模拟退火的温度控制,简单但有效。另外对于爬坡约束这类可以直接修复的约束,我会在变异完后直接做"修复"——把越界的值钳位回去并把相邻时段拉平,而不是只罚不管。
6.2 目标归一化是经常被忽略的关键
这个坑几乎每个新手都会踩。我的三个原始目标——收益大概是1200-1500千元的量级,爬坡量大概是200-400MW的量级,弃光量是0-20MW的量级。在三维目标空间里,收益维度天然就占据了绝对跨度,非支配排序时收益维度几乎决定了层级的先后,弃光量维度的差异被淹没在数值尺度里。拥挤度距离也一样,收益维度上的归一化距离远大于其他维度,多样性维持几乎只考虑收益一个方向。
所以我在evaluate里做了标准化处理,把每个目标压缩到0-1之间(用上一代种群的目标值最小值和最大值为基准做min-max归一化)。注意这个归一化必须动态更新,因为种群在进化,目标值的范围在变化,用固定的静态归一化会过时。多目标优化里所谓"尺度问题"是水很深的,但对我们这种2-4目标的问题,动态min-max归一化已经足够好用。
6.3 数值稳定性和计算效率问题
多目标进化算法要反复评价目标函数,而评估里又包含了约束判断和罚函数计算。在规模不大时要小心两个问题:一是非支配排序里两个个体目标值完全相等的情况,dominance判断需要加一个容差(比如1e-9),否则会出现互为"支配"的环,导致排序结果错乱;二是群体里出现NaN,只要任何一个目标函数里出现除零或者溢出,NaN就会在交叉变异里迅速扩散,让整代种群报废。
效率方面,我最初用纯Python写目标函数,在200代、种群100的情况下单次运行需要40多秒。后来发现瓶颈在evaluate里对全部个体做了向量化不够彻底,改用numpy数组批处理之后,直接降到5秒以内。如果调度周期更长(比如96点)或者要做滚动优化,建议用numba加一行装饰器把目标函数编译加速,效果立竿见影。
6.4 从离线优化到滚动调度的扩展思路
文章里这份代码解决的是"已知明天光伏预测曲线,提前制定全天水电计划"的离线调度问题。但实际工程里光伏预测误差会随着时间推移不断累积,一个固定的24小时计划执行到最后几小时可能已经严重偏离最优。我的做法是把离线优化改成滚动优化:每15分钟到1小时重新跑一次优化,只执行最新1-2个时段的决策,然后滚动更新。这就对算法计算速度提出了更高要求,也是为什么我会在参数里刻意控制种群和迭代次数,不盲目追求"更大更好"。
另外,如果你需要把光伏随机性也纳入建模,可以在目标函数里增加一个"最坏场景约束"或者把光伏出力曲线换成多个场景并求期望目标值,这属于随机优化调度的范畴了,需要更大的计算资源,但思路和NSGA-II框架是兼容的,替换evaluate函数即可。这也是我把目标函数和算法主循环分离的真正原因——工程里"算法"从来不是瓶颈,"模型"才是需要反复迭代的部分。
最后再分享一个小技巧:跑完优化以后,别急着把程序关掉,把Pareto前沿上那几个典型方案的水电出力曲线单独画出来对比一下,用调度员的经验和直觉去校验曲线的物理合理性。如果某个方案在凌晨时段水电出力突然满发、凌晨后立刻压回最小技术出力,那大概率是电价信号或者爬坡约束的参数设置有问题,而不是算法出了问题。数据上说得通、物理上可执行、操作上可接受,这三条都满足的Pareto解才能真正从论文走进调度大厅。