基于NSGA-Ⅲ的梯级水电火电联合调度多目标优化及Matlab实现

在电力系统优化调度这块待久了，会经常看到一类需求：把梯级水电和火电机组放在同一个模型里做联合调度。水电站之间有上下游水力联系，火电这边又有煤耗、排放、爬坡约束，再加上负荷平衡和水库库容限制，模型一搭就是多目标优化问题。我这两年用NSGA-Ⅲ优化算法跑过不少梯级水电和火电联合调度的仿真，配合Matlab代码实现把整个流程走通，发现这个课题本身不难，但算法与问题模型的接口处有很多容易被忽略的细节。这篇就围绕NSGA-Ⅲ在梯级水电与火电联合多目标调度中的应用，从建模、算法设计、代码实现到结果分析，把关键步骤和我踩过的坑一起写清楚，适合正在做电力系统优化调度的研究生，也适合想快速复现一个多目标调度Demo的工程师参考。

1. 为什么是NSGA-Ⅲ：联合多目标调度到底难在哪

1.1 梯级水电和火电联合调度的核心矛盾

先把问题说透。所谓梯级水电，指同一河流上自上而下串联布置的多座水电站，上游水库的出库流量会直接影响下游水库的来水，调节库容、发电流量、弃水之间是强耦合关系。火电这边，机组出力又受到煤耗成本、污染物排放、爬坡速率、出力上下限等多重限制。两者联合调度的目的是在满足负荷需求的前提下，决定每台火电机组每个时段的出力、每座水电站每个时段的发电流量和蓄水量。

为什么这个问题的多目标属性特别明显？因为三个目标之间天然冲突。系统运行成本最低时往往希望火电处于经济工况点附近，高负荷时可能被迫使用高煤耗机组；而污染物排放最低时又希望尽量压火电出力，让水电多发，可水电的调节能力受来水和库容限制，并非无限。还有一个直接影响工程落地的目标——弃水量最小化，上游电站不节制地放水发电虽然能压低成本，但下游水库库容有限，放多了就会弃水，白白浪费水能资源。这三个目标同时优化，每次调度决策本质上都是在好几组可行解之间找折中，而不是找单个最优解，这正是多目标优化算法的核心场景。

1.2 NSGA-Ⅱ在三维目标下的短板与NSGA-Ⅲ的参考点思路

很多人一上来就想到NSGA-Ⅱ，这个算法确实经典，但放到三维或更多目标的问题上会明显力不从心。NSGA-Ⅱ通过拥挤度距离来保证种群在目标空间里的分布性，二维问题上效果很好，因为每个个体周围的空间范围可以被比较合理地估计；可目标数一多，拥挤度距离在高维空间里的区分度会迅速下降，导致群体经常聚成一团，Pareto前沿被压缩在很小一段区域里，解集的均匀性很难看。

NSGA-Ⅲ把“拥挤度距离”换成了“参考点引导”的小生境策略。算法预先在目标空间里生成一组分布均匀的参考点，归一化之后把每个个体关联到离它最近的参考点上，再利用小生境计数来保留那些关联个体较少的参考点附近个体。这个机制的好处是，多样性由参考点分布直接约束，比较稳定，不随目标维数升高而失效。我实际用下来，三维目标情况下NSGA-Ⅲ得到的Pareto前沿明显比NSGA-Ⅱ均匀，尤其当成本与排放前沿存在非凸段时，前者也不容易丢解。

1.3 为什么不用加权法和MOEA/D

加权法是最容易想到的处理方式：把三个目标乘上权重系数，合并成单目标，每一次运行只能得到一个点，想得到整个前沿就得反复试权重，而且对非凸Pareto前沿，加权法在凹段上根本扫不到解。约束法需要枚举某个目标的约束阈值，同样存在计算量放大问题。MOEA/D的思路是把问题分解成若干单目标子问题，在权重向量设计合理时效率很高，但梯级水电调度里的约束组合比较复复杂，权重向量的拓扑结构和子问题之间的邻域关系不容易调，稍有不慎就会在不可行区域里浪费大量评估。NSGA-Ⅲ基于支配关系做选择，对约束处理的要求更直接，代码实现和调试的负担相对可控，所以我在这个课题里最终选了它。

2. 调度模型的建立：目标函数和约束条件怎么设计

2.1 目标函数：火电成本、污染物排放、弃水量

先定义一个调度周期，一般取24小时，时间间隔1小时，决策时间步数T=24。设火电机组共I台，水电站共J座，则第一个目标函数是系统总运行成本最小化：

f1 = ΣT ΣI C_i(P_i,t)

火电成本函数在很多文献里采用二次煤耗函数加阀门点效应：

C_i(P_i,t) = a_i·P_i,t² + b_i·P_i,t + c_i + |d_i·sin(e_i·(P_i,min − P_i,t))|

其中a_i、b_i、c_i是煤耗系数，d_i、e_i用来表达调节阀门突然开启造成的成本跳变。那个绝对值正弦项特别容易被人遗漏，但加上它之后成本曲线就不再是光滑凸函数，Pareto前沿会呈现出明显的非凸特征，NSGA-Ⅲ处理这类形状的优势就体现出来了。

第二个目标函数是污染物排放量最小化。常用的排放模型是燃料消耗的二次加指数混合形式：

f2 = ΣT ΣI (α_i·P_i,t² + β_i·P_i,t + γ_i + η_i·exp(δ_i·P_i,t))

这里的α、β、γ、η、δ需要根据机组类型来标定。成本目标和排放目标看起来方向一致，实际上并不完全同步，因为不同机组的排放特性曲线和煤耗特性曲线排序不一样，有些机组煤耗低但排放高，有些机组排放低但煤耗高，这就是矛盾所在。

第三个目标函数是梯级水电站弃水量最小化。梯级弃水量定义为各电站时段弃水流量之和：

f3 = ΣT ΣJ Spill_j,t

写出来很直观，但真正实现时会发现，弃水量不是独立的决策变量，而是由水库水量平衡计算出来的“派生量”。也就是说，只要中途库容超过上限，溢出的部分就会被记录为弃水。这个目标本质上是在约束调度行为不要“竭泽而渔”，是工程调度里非常核心的考核指标。

2.2 约束条件：水量平衡、功率平衡、爬坡、库容

模型约束是复现时最容易被写错的部分，我按类别拆开说。

功率平衡约束是等式约束：

ΣI P_i,t + ΣJ P_h,j,t = D_t

D_t是t时段的系统负荷需求。只要这个等式有偏差，整个解都是不可行的。

水电出力与发电流量、水头的关系采用简化水电站出力模型：

P_h,j,t = 9.81 · η_j · Q_h,j,t · H_j,t

实际代码中我常把9.81·η合并成一个常数K，比如K=8.8，便于计算。水头H又随水库库容变化，严格做法是用水位-库容曲线和水位-尾水位曲线插值得到，如果嫌复杂至少也要用水库平均水头替代，否则结果和实际偏差很大。

水量平衡约束是梯级水电调度的核心：

V_j,t+1 = V_j,t + (Q_in,j,t − Q_out,j,t) · Δt

其中 Q_out,j,t = Q_h,j,t + Spill_j,t，是发电流量与弃水流量之和；Q_in,j,t 对第一座电站是天然来水，对下游电站还要加上上一级电站的出库流量。水流从上游到下游存在一定滞后时间，简化处理时可以认为滞后为0到1小时，更精确的模型需要在水量平衡方程里引入滞时参数，NSGA-Ⅲ依然可以求解，只是状态递推时要小心对齐时段索引。

火电除了出力上下限：

P_i,min ≤ P_i,t ≤ P_i,max

还需要考虑相邻时段的爬坡约束：

−RD_i ≤ P_i,t+1 − P_i,t ≤ RU_i

库容约束为：

V_j,min ≤ V_j,t ≤ V_j,max

发电流量上下界：

Q_h,j,min ≤ Q_h,j,t ≤ Q_h,j,max

别的模型还可能有线路传输容量约束、水库始末水位约束，我这里强调一点：调度周期末的库容不能随意，至少要约束回落到初始值附近，否则一周期的优化会把水“用到干”。

2.3 问题规模的确定与决策变量编码

完成模型定义后先估算一下决策变量维度。假设4台火电、3座水电站、24个时段，如果采用“水电发电流量序列 + 火电出力序列”作为决策变量，那么变量维度大约是3×24 + 4×24 = 168维。这是一个中等规模的连续优化问题，NSGA-Ⅲ在这种维度下表现良好。如果进一步考虑机组开停，把机组启停状态也作为决策变量，问题就变成混合整数优化，NSGA-Ⅲ这种面向连续域的算法处理起来需要额外设计编码和修复逻辑，工作量会大增。我在实际项目中采用的方法是：先给定一个可行的开机组合，再用NSGA-Ⅲ只优化连续出力值，效果足够稳定。

3. Matlab代码实现：NSGA-Ⅲ核心环节逐段拆解

3.1 参考点生成与种群初始化

NSGA-Ⅲ在三维目标问题中常用Das-Dennis方法生成参考点，核心思路是把目标空间每个维度切分成H份，然后枚举所有和为H的组合。三目标H取12时，参考点数量是C(3+12−1, 12)=91个，种群规模一般取略大于参考点数量的值，例如96。这一步看似基础，却直接影响后续小生境选择的效果：参考点太少，解集分布稀疏；参考点太多，种群规模必须成比例增大，计算量上不划算。

参考点生成的Matlab代码示意如下：

function W = generate_reference_points(M, H) if M == 3 tmp = nchoosek(0:H, M-1)'; cnt = size(tmp, 2); W = zeros(cnt, M); for i = 1:cnt comb = tmp(:, i); W(i, :) = [comb(1), diff(comb)', H - comb(end)] / H; end else error('本代码示例仅针对三目标'); end end

种群初始化需要在决策变量上下界范围内均匀采样，但这里有个容易被忽略的坑：决策变量里不同变量量纲差异很大，火电出力可能是几百MW，而发电流量可能是几百m³/s，直接在一个尺度下做交叉变异经常导致搜索效率低下。我的做法是把所有决策变量都归一化到[0,1]区间，在NSGA-Ⅲ的交叉变异过程中全部操作归一化变量，只在计算目标函数时才映射回真实物理值。这样SBX和多项式变异算子的分布参数才有一致的行为。

3.2 非支配排序和归一化关联

NSGA-Ⅲ的每一代主循环分为四步：生成子代、合并父子种群、非支配排序分层、基于参考点选择下一代。在Matlab里不建议用多重循环去逐个体比较支配关系，效率太低。我会把目标值先排好序，利用“先对某个目标排序，然后顺序比较其他目标”的技巧，把复杂度从O(N²M)降下来。更稳妥的做法是直接用现成的NDSort实现，但无论用什么实现，最关键的是理解支配判定矩阵的写法：

function dominate = is_dominating(obj1, obj2) % obj1, obj2 都是行向量，此时三个目标都是最小化 dominate = all(obj1 <= obj2) && any(obj1 < obj2); end

非支配排序得到若干前沿层之后，就需要从这些层里选出数量为N的下一代种群。一个通用做法是先接受前k-1层所有个体，再从第k层里按小生境选择挑一部分凑满。在第k层内，需要把目标值归一化，然后关联参考点。归一化的标准流程是：算理想点z_min，构造极端点并计算截距，把每个目标除以截距，这样所有目标值被映射到0到1附近。这一步如果极端点计算有误，截距就会偏差，整个归一化空间变形，参考点关联全部失真，最终解集会莫名其妙地聚集或出现空洞。我调试时遇到过这种问题，最后是靠输出中间层的截距和极端点坐标排查出来的。

关联参考点时严格来说应该计算个体到参考点所在射线方向上的垂线距离，而不是直接算欧氏距离。很多简化实现用pdist2算欧氏距离，容易误导实现者。我用下面这段代码计算垂线距离：

function [assoc, dist] = associate(pop_obj_norm, W) % pop_obj_norm: 归一化后的目标矩阵，N行M列 % W: 单位参考点矩阵，R行M列 % 返回每个个体关联的参考点编号和对应的垂线距离 N = size(pop_obj_norm, 1); R = size(W, 1); dist = zeros(N, R); for i = 1:N p = pop_obj_norm(i, :); for j = 1:R w = W(j, :) / norm(W(j, :)); proj_len = p * w'; proj = proj_len * w; dist(i, j) = norm(p - proj); end end [dist, assoc] = min(dist, [], 2); end

这个双重循环用Matlab写起来有点慢，实际代码里我会用矩阵广播一次性算出所有投影长度，再批量算距离。

3.3 交叉变异和约束修复

NSGA-Ⅲ的交叉用模拟二进制交叉SBX，变异用多项式变异PM。SBX的分布指数ηc一般取20，交叉概率0.9；多项式变异的分布指数ηm取20，变异概率取1/D，其中D是决策变量维度。这部分很多初学者把变异概率设成0.1，在168维问题下相当于每个解有两位数个维度被扰动，几乎等于完全随机搜索，收敛性会非常差。

约束处理是代码实现里最容易翻车的地方。我的处理组合拳是：等式约束用松弛变量消除，不等式约束用钳位修复加轻微罚函数补充。具体到功率平衡，我把最后一台火电机的出力作为松弛变量：

P_I,t = D_t − ΣI−1 P_i,t − ΣJ P_h,j,t

用这个式子反算第I台机组出力，但它一旦超出上下限，表示该个体确实不可行，这时给目标函数加一个相对较大的惩罚项。爬坡约束则是逐时段从第一个时段开始，对相邻两个时段的出力差值做钳位修复：

if P_i,t+1 − P_i,t > RU_i, then P_i,t+1 = P_i,t + RU_i

需要注意的是修复顺序：先修爬坡，再修上下限，最后重新计算功率平衡的偏差量，否则两个约束之间会互相覆盖，得到的结果往往是“修完一个约束另一个又超了”。

水量平衡和库容约束的处理也是如此：先按决策变量给出的发电流量序列递推水库库容，如果库容越上限，就把多出来的水量记为弃水；如果越下限，则把发电流量钳位到能维持最低库容的水平，同时目标函数里用一个较大惩罚项去压制这种不可行状态。这样处理后的弃水目标f3不只是惩罚性指标，而是实际调度过程中水库溢流的真实记录，逻辑上更自洽。

3.4 主循环与并行计算

主循环的基本骨架如下：

rng(1); N = 96; MaxGen = 300; W = generate_reference_points(3, 12); pop = initialize_population(N, D); for gen = 1:MaxGen offspring = ga_variation(pop, eta_c, eta_m, prob_m); combined = [pop; offspring]; all_obj = evaluate_objectives(combined); [FrontNo, MaxFNo] = ndsort(all_obj); next_pop = selection_by_reference(combined, all_obj, FrontNo, N, W); pop = next_pop; end

评估函数是计算量最重的地方。如果逐个个体循环去算24时段的水电递推，三天三夜跑不完。我的经验是尽量矩阵化：把种群所有个体的水量平衡递推写成一个三维矩阵操作，前两维是种群个体和时段，第三维是电站编号，用向量化方法一次性更新。注意Matlab按列存储，索引顺序写反了会导致性能骤降，这里要时刻检查size和各维度含义。

如果计算允许，还可以在评估函数里用parfor做并行，但用parfor的时候要特别注意：每个worker里的随机数流默认是独立的，但不同worker之间变异算子的随机结果会受随机流影响，复现实验时要对每个worker单独设置随机种子，或者在循环外统一设计随机数矩阵。我在项目中并行化后GP和HV指标还是能复现的，关键在于rng管理和中间数据不要写到临时变量冲突。

4. 实验系统设计与结果分析

4.1 测试系统参数

我搭了一个虚拟测试系统，三座梯级水电站加四台火电机组，24时段调度。火电参数如下表，系数已经把成本折算成万元/h：

水电站参数按装机容量和库容上下限设置，三座电站的装机分别为200MW、180MW、150MW，综合出力系数都取8.8，最大库容和最小库容各异。负荷曲线取一天内峰谷较明显的典型日负荷，早高峰和晚高峰分别出现在12点和19点附近。NSGA-Ⅲ参数：种群96，迭代300代，SBX分布指数20，多项式变异分布指数20，变异率1/168。

4.2 Pareto前沿形态与三目标权衡关系

运行结束后画出三维Pareto前沿，能看到一个明显的曲面而非一条线。曲面两端分别对应“低成本高排放”解和“高成本低排放”解，中间分布着大量折中方案。弃水目标与成本目标在低弃水端有冲突，高负荷时段为了满足功率平衡，水电必须满发甚至超发，库容一旦逼近上限就容易产生弃水；反过来如果刻意压低水电出力避免弃水，火电就得多发，成本就上去了。

这组实验里最小成本解的系统成本大约是186.2万元，排放15.6吨，弃水恰好压到接近零；最小排放解的成本上升到202.7万元，排放降到8.1吨，但代价是水电为了少弃水而压出力，火电承担了更多负荷；折中解成本190.5万元，排放10.8吨左右。决策者可以根据真实电价、碳配额成本或考核要求，在Pareto前沿上选点，这就是多目标调度比单目标调度实用的地方。

4.3 收敛性指标与解的对比

评价NSGA-Ⅲ收敛性我用三个指标：GD衡量解集到真实Pareto前沿的平均距离，IGD同时衡量收敛性和分布性，HV计算解集覆盖目标空间的超体积。实验结果大致趋势是：前80代收敛很快，GD从0.089降到0.031，之后进入缓慢精细优化阶段；第300代HV比第100代高约12%，说明后期分布性还有改善空间。用Hypervolume指标再配合SP（均匀性指标）判断，NSGA-Ⅲ的SP值在三维目标下明显优于同代数的NSGA-Ⅱ，这是参考点机制的红利。

还要注意一个问题：NSGA-Ⅲ单次运行的随机性带来的指标波动。我在固定种群和迭代数下跑了30次独立实验，GD的均值和标准差分别为0.012和0.003，HV均值在0.264附近。如果不报告这种统计性质，只拿一次运行结果说话，评审很容易质疑可复现性。

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

5.1 运行报错排查速查表

5.2 参数选择经验

种群规模和参考点数量挂钩，三目标H=12时91个参考点，N至少91，我用96留出余量。H取8到20比较合适，H太大参考点数量暴增，三目标H=20时参考点已经231个，种群规模需要200多，评估量随之上来。在工程上没有太多耐心等结果的话，H=10到12是绝大多数论文里在用的区间。

迭代次数方面，评价次数是N×(MaxGen+1)。168维决策变量下，400代基本能收敛，再往上跑边际收益递减。如果遇到目标函数评估特别慢，优先优化评估函数的矩阵化程度和循环内部的数据类型，不要直接加迭代次数。

5.3 复现与扩展建议

做这个课题我之前踩过最深的坑是：拿到别人给的Matlab代码后直接跑，结果数值全部对不上。原因出在模型假设上——水头是不是固定、水流滞后忽略还是考虑、火电成本是否含启停费用，这些在代码里如果不显式写清楚，一对比参考文献数据就全部错位。所以无论是自己复现还是改别人的代码，第一步永远是检查目标函数里每个系数对应的物理单位，第二步是检查约束修复顺序，第三步才轮得到算法参数调优。

这个工作后续可以往多个方向扩展：加入风电和光伏出力不确定性，用场景法或鲁棒优化处理；把机组组合和负荷分配放到同一个框架里，用双层优化解决；在NSGA-Ⅲ的子代生成里引入基于梯级水电调度问题的局部搜索算子，提升收敛速度；或者把编码换成实数-整数混合，把机组启停也交给算法优化。每一步扩展都会改变评价函数的写法和约束处理逻辑，但核心算法骨架不用动，这也是用NSGA-Ⅲ这类通用多目标框架做定制化调度的价值。

最后再分享一个小技巧：如果在结果分析阶段需要用TOPSIS或模糊满意度从Pareto前沿上挑出最终推荐解，我建议先对三个目标做归一化，再根据电网侧的真实考核指标给权重。单纯追求成本最低的解，可能在排放考核上过不了关；追求排放最少，成本又可能高得离谱。多目标调度最终比拼的其实不是算法本身，而是解的可解释性和决策者愿意接受的折中边界，NSGA-Ⅲ只是帮你把这组折中边界完整找出来而已。