☰
基于NSGA-II的水光互补多目标优化调度Python实现
2026/10/13 11:32:32 网站建设 项目流程

先说明一下:这个项目核心就是把水光互补的优化调度问题,用非支配排序遗传算法(NSGA-II)跑多目标求解,最后用Python代码实现整套流程。说实话,这类题目在水电、新能源并网方向非常常见,但网上能直接跑通的完整代码没那么多,很多版本要么约束处理得太粗糙,要么帕累托前沿画出来丑得不忍直视。我自己做完一轮之后,把整套思路、代码结构、调参坑位都整理出来,希望能帮你少走点弯路。

为了让不知道背景的朋友也能跟上,我先把这个项目到底在做什么、为什么用NSGA-II、Python代码怎么落地讲清楚。后面章节会拆解每个设计细节,包含目标函数定义、约束处理、算法算子实现、结果怎么分析,以及我实测踩过的几个典型问题。如果你是做电力系统优化调度方向的研究生,或者刚接触多目标优化编程,这份内容可以直接按步骤参考。

1. 项目概述:水光互补调度到底在优化什么

1.1 这个模型解决什么问题

水光互补,字面意思就是水电站和光伏电站联合运行。光伏出力的特点是白天有、晚上没有,阴天一阵一阵地波动,完全看老天爷脸色。水电的优势是响应快、可调节能力强,尤其是有调节能力的水库电站,能根据负荷和光伏出力情况调整发电量。

但问题也出在这里:水电要配合光伏,就得频繁调整出力,这会带来机组运行、生态流量、发电效益之间的冲突。还有一块很现实的问题——白天光伏大发的时候,水电如果还满发,多余的电要么送不出去,要么造成弃光;晚上光伏出力归零,水电又得顶上去满足负荷。这一来一回,调度就不是简单“白天多学、晚上少学”的经验操作了,而是要考虑多个效益目标的协调问题。

这个项目就是把水光互补调度建模成一个多目标优化问题:决策变量是水电的逐时出力计划,目标之一是尽量增加光伏消纳(减少弃光,也可以理解成提高综合收益),目标之二是让水电出力过程尽量平稳、波动小。然后我用NSGA-II找出这两个目标之下的帕累托前沿,最终给出一批可供调度人员选择的方案,而不是只给一个“最优解”。

1.2 为什么必须用多目标

传统调度往往把多个目标加权成一个综合目标,比如“总运行成本最小化”。这个做法在工程上简单,但有一个致命缺陷:权重系数非常难定,而且调一次权重,解就变一次。多目标之间往往存在矛盾——你要是拼命消纳光伏,水电出力就会跟着光伏曲线大幅波动;你要是要求水电出力平稳,那某些时段就可能加大弃光。两个目标不可能同时达到最优,只能找一组“不坏”的折中方案。

这就像你想买一辆车,既要油耗低,又要动力猛,还要价格便宜。这三个指标天然冲突,加权打分最后买到的车大概率哪方面都平庸。你需要的是一批候选车型,每款在三个指标上各有侧重,然后根据实际偏好去选择。NSGA-II干的正是这件事:一次性找出多条“不互相支配”的方案,让调度人员根据日实际运行需要去挑。

我把单目标和多目标处理的区别整理成下面的表,方便新手理解:

处理方式核心思路典型问题适用情况
加权单目标法给每个目标乘以权重再求和权重选取主观,不同权重跑N次才能看趋势目标关系简单、运行规则固定的场景
约束转化法优先优化一个目标,另一个转成约束转化边界不好设,边界太严会无解目标优先关系明确的场景
多目标进化算法(NSGA-II)一次性搜索帕累托前沿计算量大一些,结果是一组解,决策成本转移给调度目标之间强冲突、需要灵活决策的场景

水光互补调度恰恰属于目标强冲突的场景,所以我毫不犹豫选了NSGA-II。

1.3 项目整体设计思路

整个项目分四层:

  • 数据层。准备负荷预测曲线、光伏预测出力曲线、水电参数(出力上下限、日发电量约束等)。我用的算例是24小时调度周期,时间分辨率为1小时,也就是一个个体长度为24。
  • 模型层。定义两个目标函数和一组约束条件。这里要特别小心,约束处理方式直接决定算法能不能收敛到可行域。
  • 算法层。用Python手写NSGA-II核心流程,包括种群初始化、快速非支配排序、拥挤度计算、锦标赛选择、模拟二进制交叉(SBX)、多项式变异、精英保留。
  • 展示层。把帕累托前沿画出来,再展示几个典型方案的出力曲线,方便从“一堆解”落到“实际调度怎么用”。

后面我会把每层的细节全部展开。

2. NSGA-II核心原理:从支配关系到帕累托前沿

2.1 理解“支配”比记算法名重要

很多人第一次接触NSGA-II,一头扎进代码里,结果没搞懂支配(dominance)是什么意思,代码看得云里雾里。这里我打个比方:

假设有三个方案A、B、C,每个方案算出来两个目标值:弃光电量和出力波动。A的弃光是2万千瓦时,波动是5;B的弃光是3万千瓦时,波动是4。你会发现A弃光更少但波动更大,B波动更小但弃光更多,这没法说A一定比B好,这就是非支配关系。但如果有个方案C,弃光1万千瓦时,波动3,比A、B都要好,那C就支配了A和B,直接淘汰掉A和B。

NSGA-II的第一步,就是把整个种群按照这种支配关系分层:最外层是没人能支配的个体,属于第一非支配层;第二层是除了第一层之外没人能支配的个体;以此类推。分层结束后,再配合拥挤度来保证解的多样性。这套机制保证了算法最后得到的不是一个“孤独的最优点”,而是一条分布均匀的帕累托前沿。

2.2 非支配排序和拥挤度具体做了什么

快速非支配排序的代码其实不复杂。每个个体要记录两件事:它支配了哪些个体,以及它被多少个体支配。先找出所有被支配数为0的个体放进第一层;然后把这一层支配的那些个体的“被支配数”减1,减到0的进入第二层,反复执行,直到所有个体分层完毕。复杂度是O(MN^2),M是目标个数,N是种群规模。目标数只有2到3个时,跑起来非常快。

拥挤度是为了防止算法把所有解都堆在帕累托前沿的某一个小区域。做法是对每一层里的个体,按某个目标函数值排序,最两端的个体拥挤度设为无穷大,中间个体的拥挤度等于相邻两个个体在这一目标上的差除以该层该目标的范围。目标函数个数多就多算几轮,最后把每个目标上的拥挤度加总。通俗理解就是:这个解周围同伴多不多?同伴越多,说明这块区域越拥挤,越应该被淘汰,因为留下来不利于保持解的分散性。

2.3 精英保留策略为什么能保住好解

NSGA-II的流程是父代种群做选择、交叉、变异,生成子代种群,然后父代和子代合并,总规模变成2N,再对这个大种群做非支配排序和拥挤度排序,取前N个作为下一代。这样做的好处是:就算父代里出现了某方面极好的解,交叉变异弄出的后代把父代全弄没了,好基因也不会直接丢掉,因为父代个体还保留在合并种群中参与竞争。

我在实际跑项目时感受很明显,如果不做精英保留,只拿子代当下一代,目标函数很容易越跑越差。原因很简单,交叉和变异有一定的破坏性,好解可能生出来的孩子反而不行。精英保留相当于给每一代的好方案买了一份“保险”。

3. Python代码实现:环境、建模、算子

3.1 用到的库与运行环境

我用的环境是Python 3.10,核心库只用到numpy和matplotlib。优化调度本身不太吃机器配置,24维变量、100个种群、跑200代,在普通笔记本上几十秒就能跑完。如果你还没装numpy和matplotlib,直接在终端执行:

pip install numpy matplotlib

这里有一个新手容易踩的坑:如果机器上有多个Python环境,装包的时候要确认pip属于当前正在使用的环境。实在不放心就用python -m pip install numpy matplotlib。我见过太多朋友在Jupyter里能import,但命令行脚本一跑就ModuleNotFoundError,十有八九是环境错位。这一类问题我放到第五部分再说,这里先继续讲核心内容。

3.2 问题建模:决策变量、目标函数、约束条件

这个项目里我把决策变量定义为水电站24小时的出力序列,用x表示,x[t]代表第t小时的水电出力,单位是MW。这个选择很直觉,因为水电是我们唯一能主动调度的对象。

目标函数我定义成两个,全部做最小化处理:

  • F1:弃光电量尽量小。公式可以写为sum(max(load[t] - x[t], 0))和sum(max(pv[t] - (load[t] - x[t]), 0))的组合,也就是缺电量和弃光量折算成一个综合“运行代价”。
  • F2:水电出力波动尽量小。我直接累加相邻时段水火电出力的差平方,也就是sum((x[t] - x[t-1])^2)。差平方比差绝对值更好用,因为平方能更明显地区分波动大的方案,让算法更积极地寻找平稳曲线。

约束条件有三个:

  • 水电出力上下限:x_low <= x[t] <= x_high,运行安全是底线。
  • 日发电量约束:sum(x[t]) * delta_t 不得低于某个最小值。这是考虑到水库有来水,日发电量不能太少,否则水库水位会异常偏高,甚至要溢流。
  • 时段供需平衡约束:x[t] + pv[t] 必须尽量接近负荷 load[t]。如果小于负荷就缺电,如果大于负荷就弃电,这两者都会进入F1的代价里。

至于光伏,我简化处理为预测值已知,并且允许弃光。这样不需要额外把光伏接入方式、逆变器控制之类的内容塞进模型,聚焦在算法本身。做研究或者写课程作业的话,这个简化完全够用。

这里补充一下我在建模时吃过亏的地方:一开始我把F1写成单纯的“最大化发电收益”,F2写成“波动最小”,结果发现两个目标的量纲差太远,帕累托前沿拐得非常厉害,算法退化到只优化F1的程度。后来改成把F1设计为“弃光+缺电代价最小化”,让两个目标都统一成越小越好,前沿形状就正常了许多。

3.3 核心算法代码与实现要点

我直接放核心代码。为了不过度堆砌,我把NSGA-II封装成了一个类,重点方法都保留。类结构如下,你可以按需补全读取数据的部分。

import numpy as np import matplotlib.pyplot as plt class NSGA2: def __init__(self, n_var=24, pop_size=100, generations=200, lb=None, ub=None, pv=None, load=None): self.n_var = n_var self.pop_size = pop_size self.generations = generations self.lb = lb if lb is not None else np.full(n_var, 50.0) # 水电出力下限 self.ub = ub if ub is not None else np.full(n_var, 300.0) # 水电出力上限 self.pv = pv # 光伏预测曲线,长度n_var self.load = load # 负荷曲线,长度n_var def init_population(self): pop = np.random.uniform(self.lb, self.ub, (self.pop_size, self.n_var)) return pop def evaluate(self, x): # 目标1:弃光缺电综合代价 diff = self.load - x - self.pv shortage = np.maximum(diff, 0).sum() # 缺电量 surplus = np.maximum(-diff, 0).sum() # 弃电量 f1 = shortage + 1.5 * surplus # 弃光惩罚高一点 # 目标2:水电出力波动平方和 f2 = np.sum(np.diff(x)**2) return [f1, f2] def dominates(self, a, b): # a支配b的条件:a所有目标不比b差,且至少一个严格更好 a_f1, a_f2 = self.evaluate(a) b_f1, b_f2 = self.evaluate(b) return (a_f1 <= b_f1 and a_f2 <= b_f2 and (a_f1 < b_f1 or a_f2 < b_f2)) def fast_non_dominated_sort(self, pop): # 返回每一层的个体索引列表 n = len(pop) domination_count = np.zeros(n, dtype=int) dominated_solutions = [[] for _ in range(n)] fronts = [[]] for i in range(n): for j in range(i+1, n): if self.dominates(pop[i], pop[j]): domination_count[j] += 1 dominated_solutions[i].append(j) elif self.dominates(pop[j], pop[i]): domination_count[i] += 1 dominated_solutions[j].append(i) for i in range(n): if domination_count[i] == 0: fronts[0].append(i) k = 0 while fronts[k]: next_front = [] for i in fronts[k]: for j in dominated_solutions[i]: domination_count[j] -= 1 if domination_count[j] == 0: next_front.append(j) k += 1 fronts.append(next_front) return [f for f in fronts if f] def crowding_distance(self, front_idx, pop): distance = dict.fromkeys(front_idx, 0.0) # 目标数量这里按2处理,如果有更多目标就循环处理多个目标 for m in range(2): # 按第m个目标排序 order = sorted(front_idx, key=lambda i: self.evaluate(pop[i])[m]) distance[order[0]] = float('inf') distance[order[-1]] = float('inf') fmin = self.evaluate(pop[order[0]])[m] fmax = self.evaluate(pop[order[-1]])[m] if fmax == fmin: continue for k in range(1, len(order)-1): distance[order[k]] += (self.evaluate(pop[order[k+1]])[m] - self.evaluate(pop[order[k-1]])[m]) / (fmax - fmin) return distance def selection(self, pop): idx1 = np.random.randint(0, len(pop)) idx2 = np.random.randint(0, len(pop)) # 通过排序后的锦标赛选择简化实现 # 这里只做简单二元锦标赛,完整实现需要携带支配关系信息 f1_1, f2_1 = self.evaluate(pop[idx1]) f1_2, f2_2 = self.evaluate(pop[idx2]) if self.dominates(pop[idx1], pop[idx2]): return pop[idx1].copy() elif self.dominates(pop[idx2], pop[idx1]): return pop[idx2].copy() else: return pop[idx1].copy()

这段代码里dominates方法和evaluate方法耦合在了一起,为了展示把支配比较简化为直接评估新旧个体。在实际的完整实现中,更高效的做法是提前把整个种群的目标函数值向量化计算,一次评估完所有个体的F1、F2,再做排序和拥挤度计算。我建议你用向量化方式实现,否则每跑一代要评估上万次函数,虽然本例依然能跑,但架构不干净,后续扩展到更复杂约束会很痛苦。

接下来是SBX交叉和多项式变异的实现。这两个算子是NSGA-II区别于普通遗传算法的关键。

def sbx_crossover(self, parent1, parent2, eta=20): child1 = np.zeros_like(parent1) child2 = np.zeros_like(parent2) for i in range(self.n_var): if np.random.rand() < 0.9: u = np.random.rand() if u <= 0.5: beta = (2 * u) ** (1 / (eta + 1)) else: beta = (1 / (2 * (1 - u))) ** (1 / (eta + 1)) child1[i] = 0.5 * ((1 + beta) * parent1[i] + (1 - beta) * parent2[i]) child2[i] = 0.5 * ((1 - beta) * parent1[i] + (1 + beta) * parent2[i]) else: child1[i] = parent1[i] child2[i] = parent2[i] return np.clip(child1, self.lb, self.ub), np.clip(child2, self.lb, self.ub) def polynomial_mutation(self, individual, pm=0.1, eta=20): mutated = individual.copy() for i in range(self.n_var): if np.random.rand() < pm: u = np.random.rand() if u < 0.5: delta = (2 * u) ** (1 / (eta + 1)) - 1 else: delta = 1 - (2 * (1 - u)) ** (1 / (eta + 1)) mutated[i] += delta * (self.ub[i] - self.lb[i]) return np.clip(mutated, self.lb, self.ub)

SBX交叉让子代能够根据父代差异进行有方向性的搜索,eta控制后代贴近父代的程度。eta越大,子代越集中在父代附近,收敛更精细但探索性弱;eta越小,子代探索范围越广但容易震荡。我建议先把eta设成20,这是NSGA-II论文里的经典取值,之后再看前沿分布调整到10到30之间。

多项式变异的作用是防止种群陷入局部最优。注意pm不要设太高,太高会让种群变成随机游走。变量维度是24时,pm取0.1左右已经足够,每个个体平均也就变异2到3个维度。

3.4 主循环:整合精英保留

def run(self, verbose=True): pop = self.init_population() for gen in range(self.generations): offspring = [] while len(offspring) < self.pop_size: p1 = self.selection(pop) p2 = self.selection(pop) c1, c2 = self.sbx_crossover(p1, p2) c1 = self.polynomial_mutation(c1) c2 = self.polynomial_mutation(c2) offspring.append(c1) if len(offspring) < self.pop_size: offspring.append(c2) combined = np.vstack([pop, np.array(offspring)]) fronts = self.fast_non_dominated_sort(combined) new_pop = [] for front in fronts: if len(new_pop) + len(front) <= self.pop_size: new_pop.extend(combined[front]) else: dist = self.crowding_distance(front, combined) sorted_front = sorted(front, key=lambda i: dist[i], reverse=True) need = self.pop_size - len(new_pop) new_pop.extend(combined[sorted_front[:need]]) break pop = np.array(new_pop) if verbose and (gen % 20 == 0 or gen == self.generations - 1): f1s = [self.evaluate(ind)[0] for ind in pop] f2s = [self.evaluate(ind)[1] for ind in pop] print(f"gen {gen}: F1 range [{min(f1s):.4f}, {max(f1s):.4f}], F2 range [{min(f2s):.4f}, {max(f2s):.4f}]") return pop

主循环的核心就是“合并-排序-优选”。每次生成子代后,新旧种群合并成2N规模,然后按非支配层一层一层往下一代里填,快填满时用拥挤度排序选拥挤度大的个体补进去。这一步就是所谓的精英保留,也是NSGA-II的名字里“非支配排序”的来源。

4. 案例实测与参数调整经验谈

4.1 一个能复现的24小时算例

我构造的算例是某个含水电和光伏的区域电网简化版。负荷曲线在白天有两个高峰,光伏曲线中午高、早晚低。水电出力下限设为80MW,上限设为280MW,日发电量约束设为平均出力不低于160MW。光伏峰值180MW,整体装机规模不大,主要用来模拟“白天光伏挤占水电空间,晚上水电重新顶上”的典型场景。

跑200代、种群规模100,得到帕累托前沿后我发现几个规律。当F2比较小,也就是水电出力平稳的时候,弃光和缺电很明显地集中在中午时段。因为中午光伏大发,负荷并没有同步大涨,水电又出于平稳考虑不肯快速压下来,多余电力只能弃掉。当F1比较小,也就是系统经济性优先的时候,水电出力曲线会死死贴着“负荷减光伏”的净负荷曲线走,波动值变大,但几乎不弃光不缺电。这两种方案都不算错,区别在调度偏好。

4.2 帕累托前沿图怎么看、怎么用

跑完之后把每个个体的两个目标值打点,x轴是F1,y轴是F2,左下角的点最理想。但实际前沿是一条向右下方倾斜的曲线。前沿左端代表强经济性方案,几乎不弃光;前沿右端代表强稳定性方案,水电几乎平稳,靠弃光换稳定。

实际工程中怎么选?如果当天系统备用容量充足、上网通道宽裕,可以选F1小一点的方案,多消纳光伏;如果当天水电有生态流量要求,或者机组不适合频繁调节,就选F2小一点的方案。这恰恰是多目标优化相比单目标最大的价值——不是替你做一个决定,而是给你看完整的抉择空间。

4.3 参数调优的几条经验

种群规模和迭代次数是我最先调整的两个参数。刚开始我用50个个体跑100代,前沿明显稀疏,中间段还有缺口。增加到150个个体、300代之后,前沿平滑了很多。继续加大到400代,边际改善就不明显了。对你来说,如果只做教学演示,100个个体、150代已经够用;如果是想发表在论文里,建议跑300代以上,并且至少重复5次独立实验取代表性前沿。

交叉概率的调整我比较谨慎。SBX交叉概率固定在0.9,不太会出问题。变异概率pm从0.05调整到0.2试过,pm=0.1时前沿分布最均匀。pm太小时,种群容易陷在某个前沿段里出不去;pm太大时,收敛变慢,前沿拖得很长但点很散。记住一个原则:探索和开发要平衡,探索由变异负责,开发由交叉和精英保留负责。

还有一类问题来自约束处理。由于日发电量约束的存在,很多随机生成个体完全不满足要求。我的做法是对可行不可行个体都得参与非支配排序,但当个体违反约束时,把目标函数值加上一个很大的罚项,这样既保证算法继续搜索,又不会让不可行解占据前沿位置。罚项系数也要调,太大会导致所有不可行解之间无法区分优劣,搜索效率下降;太小又会让破方案混进前沿。我调下来感觉罚项取最大目标值的1到2倍比较顺手。

5. 常见问题与排查技巧实录

5.1 等式约束难守,怎么办

如果你的模型里有严格的功率平衡等式约束,也就是水电加光伏必须等于负荷,直接塞进NSGA-II会让几乎所有的交叉变异后代都变成不可行解。这个问题我试过两种解法。

一种是把等式约束做松弛,允许实际出力和负荷之间有小偏差,偏差进目标函数里惩罚。这样可行域从一条线变成一个窄带,进化算法才能在里面操作。

另一种方法是变量重定义,把其中一个时段的水电出力改成其他时段算完之后的“平衡项”,也就是最后时刻出力由公式反算,再检查上下限。这个思路更巧妙,但会把问题从24维变量变成23维变量,而且反算出来的值可能远超约束范围,需要额外修复。

我的经验是:如果是首次跑通流程,用第一种松弛法;如果追求精度和工程可用性,用变量重定义加修复法。我在教学版本里选择了松弛法,因为代码简单,逻辑也透明。

5.2 种群过早收敛,前沿分布不均匀

如果你跑第一代的时候就发现所有点都挤在一个狭长区域,后面几代也基本不变,大概率是初始化出了问题。随机初始化如果只在一个很窄的范围里生成个体,种群多样性从一开始就没了。我用uniform(lb, ub)初始化,理论上范围够宽,但如果你把lb和ub设成同一个值,或者lb太接近ub,就会踩坑。

还有一种情况是选择压力太大。锦标赛选择如果不小心实现了“每次从全种群选最优”,很快整个种群就同质化了。标准做法是二元锦标赛,每次只随机抽两个比,赢的那个参与下一代,这能保留足够的多样性。

5.3 代码层面的“跑不动”问题

新手最常遇到的是:明明代码照着写,却拖得非常慢。主要因为非支配排序里频繁重复调用evaluate方法,而evaluate里又反复做数组差、条件筛选。最快的优化办法是把每个个体的目标函数提前算好缓存到列表里,排序和拥挤度都直接查缓存,只有当种群换代后才更新缓存。这个小优化能把运行时间压缩到原来的十分之一以下。

还有numpy版本太旧导致的运行问题,建议先把numpy升级到最新版。matplotlib画图时如果中文标签变成方框,确认一下系统有没有中文字体,或者在代码开头设置plt.rcParams['font.sans-serif'] = ['SimHei']。这些都是很小的点,但经常卡住人。

最后再分享一个我惯用的可视化技巧:在跑完帕累托前沿之后,别只画散点图,把前沿左端、中端、右端各挑一个方案,分别画出它们的水电出力曲线和光伏出力曲线叠在一起。这样一眼就能看出“经济性优先”和“稳定性优先”在实际调度曲线上长什么样,写报告或答辩的时候比单纯丢一张前沿图清晰得多。

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询