做电力调度优化的同行多少都遇过这种场景系统里既有火电也有梯级水电火电煤耗曲线是非线性的水电又受来水、库容和上下游关系的多重限制两边纠缠在一起靠人工或者单目标优化很难给出让人满意的方案。我最近在MATLAB里用NSGA-III算法把梯级水电和火电机组的联合调度重新做了一遍目标就两个火电总燃料成本最小、梯级总耗水量最小。这篇文章把从问题建模到MATLAB代码落地的全过程完整拆开讲包括约束怎么处理、参考点怎么生成、Pareto前沿怎么选解以及我实际调试中踩过的不少坑希望能帮到正在做这类课题的朋友。1. 问题建模把目标、约束、变量先写清楚1.1 为什么水火要放在同一个框架里优化单纯做火电调度问题相对好解负荷给定之后按经济调度往下分配各台机组出力就行约束也主要是机组上下限和爬坡。但一旦把梯级水电加进来事情就变复杂了。水电出力的“燃料”是水而水的状态是动态的——上游放水之后要经过一段时间才变成下游水库的入库库容多了就可能弃水少了又满足不了下一个时段的负荷。这个问题的本质是在满足系统负荷的前提下火电和水电各自的出力应该如何分配。用多了火电成本高但水留下来了用多了水电省了燃料但水耗上升且梯级之间还会相互影响。两个目标天然冲突没法用一个简单的权重压成一个标量来优化——因为权重怎么定本身就是个难题。所以我选择了多目标优化的路子在同一个框架里同时优化成本和水耗最后得到一组互为折中的解让决策者根据实际来水情况、燃料价格去挑。1.2 两个目标函数的数学表达第一个目标火电总燃料成本最小化。火电机组的耗量特性通常用二次函数拟合也就是$$F_1 \sum_{i1}^{N_T} \sum_{t1}^{T} \left( a_i P_{i,t}^2 b_i P_{i,t} c_i \right)$$其中$P_{i,t}$ 是第 $i$ 台火电机组在第 $t$ 个时段的出力$a_i$、$b_i$、$c_i$ 是煤耗曲线系数。这三个系数不是拍脑袋定的它们来自机组的实际热力试验或历史运行数据拟合单位通常对应为元/MW²·h、元/MW·h 和元/h。真实项目里如果还要计算启停成本可以在目标函数后面加一个分段惩罚项我这次为了聚焦优化问题本身没有把机组启停作为离散变量考虑进去。第二个目标梯级水电总耗水量最小。水电站的出力和流量之间的关系可以用经典公式表示$$P 9.81 \cdot \eta \cdot Q \cdot H$$$Q$ 是发电流量$H$ 是净水头$\eta$ 是综合效率系数。在简化模型里如果假设净水头近似恒定那水电出力基本和发电流量成正比。这时第二目标可以写成$$F_2 \sum_{j1}^{N_H} \sum_{t1}^{T} Q_{j,t} \cdot \Delta t$$也就是整个调度周期内所有梯级水库的总放水量。两个目标放在一起看就有意思了如果我多放水发电火电就可以少出力$F_1$ 下降但 $F_2$ 上升反过来如果勒紧水龙头火电就得顶上$F_1$ 上升而 $F_2$ 下降。这就是典型的矛盾关系Pareto 前沿会是一条向左下凸的曲线或者说一条从左上到右下倾斜的曲线。看到这种结果曲线基本可以判断模型建对了要是两条目标同时变大变小那大概率是约束写错了。1.3 约束条件哪些坑不能踩多目标进化算法里约束不会被“硬性”保证它只是决定解是否可行因此在建模阶段就要把约束列全。我用到的约束大概分四类。功率平衡约束这是硬约束中的硬约束每个时段系统总出力必须等于负荷。$$\sum_{i1}^{N_T} P_{i,t} \sum_{j1}^{N_H} P_{H,j,t} P_{load,t}$$火电机组约束包括出力上下限$$P_{i,\min} \le P_{i,t} \le P_{i,\max}$$还有爬坡约束$$\left| P_{i,t} - P_{i,t-1} \right| \le R_i$$我曾经在初期版本里忽略了爬坡约束结果算法给出的解非常“漂亮”成本极低但相邻时段出力跳变几十上百兆瓦实际生产根本执行不了。这个约束必须加。梯级水电约束是另一类包括库容上下限、发电流量上下限、以及梯级水量平衡方程$$V_{j,t1} V_{j,t} \left( I_{j,t} Q_{up,j,t} - Q_{j,t} \right) \cdot \Delta t$$这里 $I_{j,t}$ 是当地自然来水$Q_{up,j,t}$ 是上游水库的出库流量。对于梯级第二个水库来说$Q_{up}$ 就是第一个水库的出力泄流。还有一类容易忽略的是水电出力和水头的耦合。如果模型里把净水头当成恒定值那出水流量和出力的关系就简单了但如果要考虑水位变化对水头的影响就得把库容映射成水位再用水位差算水头计算量会上一个台阶。我这次采用的是恒水头简化版本后续再扩展。1.4 决策变量怎么组织决策变量不是火电出力加水电出力那么简单因为水电出力是间接量真正的操作量是各水库的发电流量。我采用的编码方式是$$X \left[ P_{1,1}, \cdots, P_{N_T,T}, \ Q_{1,1}, \cdots, Q_{N_H,T} \right]$$也就是一个维度为 $N_T \times T N_H \times T$ 的实数向量。比如两台火电机组、两个梯级水库、24个时段维度就是 $2 \times 24 2 \times 24 96$。这在高维优化问题里算中等规模NSGA-III 处理起来没有压力。如果直接把水电出力当作决策变量那水量平衡就无法闭合成环梯级上下游关系也体现不出来所以我坚持用流量作为水侧的核心变量。2. NSGA-III算法为什么选它怎么理解2.1 NSGA-II的困境在哪里NSGA-II 是很多人的入门算法它的核心是两板斧非支配排序保证收敛性拥挤距离维持多样性。处理两个目标时拥挤距离确实够用解的分布也还行。但问题是一旦目标数到了三个以上拥挤距离就明显力不从心。原因在于“拥挤距离”只反映某个解周围有多少邻居当目标太多时空间体积指数膨胀距离计算变得不敏感解会很快堆到高维空间的一小片区域整个前沿面展开不了。NSGA-III 正是为了修复这个问题而诞生的。2.2 NSGA-III的核心流程NSGA-III 和 NSGA-II 的骨架几乎一样都是先做非支配排序再逐层选择。不同的是最后一步——从临界层选个体进入下一代时NSGA-III 不用拥挤距离而是用预先设定的参考点来引导选择。流程可以拆成五步对当前种群做非支配排序得到多个非支配层。从第一层开始依次选入下一代直到某一层放不下为止。对最后一个临界层里的所有个体做目标值归一化。把每个个体关联到距离最近的参考点。统计每个参考点被选中的个体数量优先保留那些“小生境”里个体少的参考点附近的解。第五步是整个算法的灵魂。它相当于告诉算法每个参考点代表一块区域区域里已经有太多解的就别加了优先往空点、少点的方向补人。这样解集就不会挤作一团。2.3 参考点是如何生成的参考点的生成方式是系统性的。假设目标数为 $M$每个目标方向等分为 $p$ 份那么生成的参考点数量为$$H C_{Mp-1}^{M-1}$$举个例子两个目标、每维等分 99 份参考点数是 100三个目标、每维等分 12 份参考点数就是 $C_{14}^{2} 91$。每一个参考点都在目标空间的单位单纯形上也就是说这些点的坐标之和总是 1。二维时参考点就是一根线段上的均匀刻度很好想象。三维时是一个正三角形上的网格点。更高维时这些点分布在单纯形面上彼此间距尽量均匀这就为解集多样性提供了“锚点”。2.4 和主流多目标算法的快速对比我干脆列个表格把几个主流算法的特点放一起方便你选型时直接对照。算法多样性机制适用目标数主要特点NSGA-II拥挤距离2个勉强3个经典实现简单3目标以上效果变差NSGA-III参考点小生境3个以上也适合2个解集分布均匀参数多一个参考点划分MOEA/D权重向量分解2-3个速度快但权重选取对结果影响较大SPEA2环境密度估计2个历史表现不错但计算复杂度偏高对于水火联合调度这种两到三个目标的问题NSGA-III 表现出了比 NSGA-II 更好的均匀性这是我在多次对比实验后的真实感受。如果你已经写了 NSGA-II 代码改造成 NSGA-III 的成本其实不大主要就是环境选择部分换成参考点逻辑。3. MATLAB代码实现一步一步拆开讲3.1 代码框架总览一套完整的代码我按功能拆成了以下几个文件main_nsga3主程序负责参数设置、循环和结果绘制initialize_population种群初始化evaluate_objectives目标函数评估non_dominated_sort非支配排序generate_reference_points生成参考点environmental_selection环境选择含归一化和关联sbx_crossover和polynomial_mutation交叉与变异算子主循环的逻辑是父代种群产生子代然后合并父代和子代做非支配排序和环境选择优选出下一代。这个“精英策略”保证了不错的最优解不会被丢掉。我用伪代码描述一下核心循环结构for gen 1:maxGen offspring []; for i 1:N/2 [p1, p2] tournament_selection(pop); [c1, c2] sbx_crossover(p1.x, p2.x, eta_c, lb, ub); c1 polynomial_mutation(c1, eta_m, lb, ub); c2 polynomial_mutation(c2, eta_m, lb, ub); offspring [offspring; gene(c1); gene(c2)]; end pop_all [pop; offspring]; fronts non_dominated_sort(pop_all); pop environmental_selection(pop_all, fronts, ref_points); end看到这个结构你就明白了NSGA-III 的框架不复杂难点都在那个environmental_selection里。3.2 种群初始化给算法一个好的起点我用一个结构体数组来存个体pop(i).x zeros(1, nvars); pop(i).f zeros(1, M); pop(i).cv 0; % 约束违反量 pop(i).rank 0;初始化时火电出力在各机组上下限内随机生成发电流量在水库允许范围内随机生成。但有一个环节很关键初始种群如果不做处理很多个体无法满足功率平衡后续罚函数计算会让整个搜索变得异常缓慢。我习惯在初始化时做一次“补差修正”——随机选一台火电机组把各时段的功率差额补上去for t 1:T diff load(t) - sum(P_H(:, t)) - sum(P_T(:, t)); % 选一台有余量的火电按差额修正 i randi(NT); P_T(i, t) P_T(i, t) diff; P_T(i, t) max(Pmin(i), min(Pmax(i), P_T(i, t))); end这样生成的初始解虽然不保证完全满足所有约束但至少功率平衡的方向不会歪得离谱收敛速度会快很多。3.3 目标函数和约束评估的细节评估函数返回两个值目标向量和约束违反量。核心代码如下function [f, cv] evaluate_objectives(x, data) % 解码决策变量 P_T reshape(x(1:data.NT*data.T), data.NT, data.T); Q reshape(x(data.NT*data.T1:end), data.NH, data.T); % 计算水电出力简化恒水头模型 P_H zeros(data.NH, data.T); for j 1:data.NH P_H(j, :) 9.81 * data.eta(j) * Q(j, :) * data.H(j) / 1000; end % 目标1火电燃料成本 f1 sum(data.a .* P_T.^2 data.b .* P_T data.c, all); % 目标2梯级总耗水量 f2 sum(Q, all) * data.dt; f [f1, f2]; % 约束违反量功率平衡 cv sum(abs(sum(P_T, 1) sum(P_H, 1) - data.load));注意我在约束违反量里暂时只写了功率平衡这只是最小版本。完整项目里还要把火电爬坡约束、梯级水量平衡、库容上下限全部叠进去。约束违反量算出来后在环境选择时直接作为一个排序维度违反量为零的解优先生。3.4 非支配排序和参考点生成的代码逻辑非支配排序在很多工具箱里都能找到现成实现但自己做一遍更能理解逻辑。基本思路是两层循环对种群中任意两个解比较目标向量看谁支配谁统计每个解被支配的次数和被它支配的解集合。function [front] non_dominated_sort(pop) N length(pop); n_dom zeros(N, 1); dominated_set cell(N, 1); front cell(1); for i 1:N for j i1:N if dominates(pop(i), pop(j)) dominated_set{i} [dominated_set{i}; j]; n_dom(j) n_dom(j) 1; elseif dominates(pop(j), pop(i)) dominated_set{j} [dominated_set{j}; i]; n_dom(i) n_dom(i) 1; end end end % n_dom 0 的个体属于第一前沿依次更新其余层的被支配计数 end参考点的生成相对独立两目标的时候最简单function ref generate_reference_points(M, p) if M 2 H p 1; ref zeros(H, M); for i 0:p ref(i1, :) [i/p, 1 - i/p]; end end % 三目标及以上需要用组合递推生成单纯形上的均匀点 end3.5 环境选择整个算法的核心所在环境选择是把 NSGA-II 改成 NSGA-III 的关键差异点。我把它拆成三步讲。第一步归一化。先找到当前种群每个目标的最小值构成理想点然后把所有目标值减去理想点。接着要算出每个目标对应的“截距”用来把不同目标拉到一个量纲上。这一步的目的是让火电成本和水耗这两个单位完全不同的目标在关联参考点时能够公平计算。第二步关联。把每个个体和所有参考点做内积计算找到离它最近的参考点。本质上就是计算个体到每条参考线之间的垂直距离距离最小的就是它的归属。第三步小生境保持。统计每个参考点被分配到的个体数量找到目前个体数最少的参考点然后从这个参考点关联的候选解中挑一个进入下一代。这样做会让算法持续往“有空间”的方向填充而不是把选择压力集中在少数区域。核心代码逻辑大致是function pop_new environmental_selection(pop_all, fronts, ref) pop_new []; layer 1; while length(pop_new) length(fronts{layer}) N pop_new [pop_new; fronts{layer}]; layer layer 1; end % 处理最后一个入不满的临界层 last_layer fronts{layer}; % 对当前已选解构造小生境再关联临界层解 niche_count prepare_niche(pop_new, ref); [chosen, niche_count] select_from_last_layer(last_layer, ref, niche_count, N - length(pop_new)); pop_new [pop_new; chosen]; end这段逻辑写对之后NSGA-III 的骨架就通了。3.6 主程序的输出与结果整理迭代结束后我把最终种群里非支配解集输出同时保存每一代的Pareto前沿快照方便看收敛过程。画图时直接用 scatter 画两个目标的散点图横轴是火电成本纵轴是梯级耗水量。如果算法收敛正常你会看到一条从左下向右上延伸的曲线——左下是省成本但费水的解右上是费成本但省水的解。除了画前沿我还习惯把每个个体对应的调度方案单独存下来这样后面选中的某个折中解可以直接取出逐时段的火电出力、水库放水流量和库容变化曲线用于决策支撑。4. 案例测试24时段联合调度4.1 测试系统参数我用一个简化系统来验证算法两台火电机组、两个串联水库、24个时段时段间隔1小时。系统参数如下参数机组1机组2a系数元/MW²·h0.0020.004b系数元/MW·h2015c系数元/h300200出力下限MW10080出力上限MW400350爬坡上限MW/h5060水库部分简化为两个库容区间都在 50~300 万立方米之间自然来水按典型的枯/丰中间值设置让这个问题既不过分简单也不至于太苛刻。日负荷曲线我用的是一条典型的工作日曲线凌晨负荷低白天出现两个高峰。负荷峰值约900MW低谷约500MW。4.2 算法参数设置NSGA-III 的参数经验值如下种群规模200最大迭代次数300参考点每维划分份数99对应100个参考点交叉概率0.9变异概率1/nvars这里约1%SBX分布指数20多项式变异分布指数20这些参数不需要精确到小数点先把大概范围定下来然后看收敛曲线是否平稳。如果200代以后Pareto前沿几乎不再变化说明迭代代数够了。如果还在明显移动就增加到400代。4.3 Pareto前沿解读与方案选择运行结束后典型的结果是一条平滑的凸曲线。曲线上的每一个点都代表一种调度策略。曲线左下角意味着火电出力被压得很低水电承担了大部分负荷水耗自然高曲线右上角则相反火电满发或接近满发水库很“省”但燃料成本很高。那么问题来了曲线上这么多点该选哪个作为最终调度方案实际项目中我不会只看一个点而是用两种方式选。一种是根据当天的来水预报和燃料价格换算一个“权衡比例”在线性化处理下选离决策线最近的解。另一种更稳妥选出几个典型方案比如最低成本方案、最低水耗方案和用模糊TOPSIS综合评估出的折中方案再让调度人员在这几个候选里做最终判断。4.4 和NSGA-II的对比实验我跑了同样的系统仅把环境选择从拥挤距离替换成参考点其余参数保持不变对比两个指标解集在目标空间中的分布均匀程度和最终超体积值。实际结果和预期一致NSGA-III 的Pareto前沿分布更均匀没有出现某些区段解特别密集、另一些区段却稀疏的情况。超体积值两块算法差异不算大但均匀性这项上 NSGA-III 优势明显尤其在考虑更多目标时会更突出。5. 常见问题与调试经验5.1 功率平衡约束老被破坏这是最容易出现的问题。原因有两个一是初始种群质量太差在随机初始化下很难产生满足功率平衡的个体二是罚函数权重设置不合理导致进化过程不把功率平衡当回事。我的处理方法是“修复惩罚”双管齐下。在评估任何个体之前先做一次补差修正把火电出力拉回负荷平衡的状态修正之后如果仍然越限再计算惩罚量。这个方法比纯罚函数稳健得多。5.2 梯级水量平衡的时序耦合梯级水库之间的耦合是最容易在代码里写错的点。第二个水库当天的入库有一部分来自第一个水库当天的放水。这个“当天到当天”的假设已经做了简化实际河流要考虑流达时间。代码实现时要注意更新的顺序先更新最上游水库再按顺序往下游推。千万不能并行更新否则下游用到的上游出库流量还是旧值算出来的库容变化就会混乱。我第一次写的时候就在这里栽了跟头整条水量曲线歪得离谱还以为算法出了问题排查很久才发现是更新顺序错了。5.3 收敛太慢或者解集太散解集太散通常是参考点划分粒度太大。两目标问题的话把 p 从 20 提到 99参考点数量会从 21 涨到 100多样性提升非常明显。收敛太慢则大概率是种群规模太小或者变异概率过低。我常用的调试手段是先输出每一代Pareto前沿的HV值并画收敛曲线如果曲线到后期还是大斜率下降说明还没收敛加大迭代代数如果前期就平坦、但解集分布很差那多半是参考点设置或选择压力出了问题。5.4 MATLAB运行太慢怎么提速NSGA-III 的瓶颈主要在环境选择里的归一化和关联计算。这部分如果写成双层循环种群200、100个参考点时还好一旦规模上去就明显卡顿。提速技巧有三个一是把目标向量计算全部向量化尽量避免逐时段循环二是参考点在每次迭代前一次性生成好不要重复计算三是关联时用矩阵运算替代逐个体内积。我在重构代码时做了这几处优化整体速度提升了差不多三倍。如果还要更快可以直接把评估函数丢进并行循环里跑反正评估是纯函数天然支持并行。6. 算法扩展的几种可行方向这套框架跑通之后扩展起来思路非常清晰。第一个方向是加第三个目标。比如把火电产生的污染物排放量作为独立目标加进去变成成本、水耗、排放三目标优化。这正是 NSGA-III 最有优势的场景因为三目标下它的参考点机制比 NSGA-II 的拥挤距离强太多。第二个方向是考虑不确定性。来水预报和负荷预测都不可能完全准确把不确定性的惩罚项或者场景集放进去问题就会变成一个随机优化或者鲁棒优化问题。NSGA-III 的框架不需要大改主要改评估函数和约束处理逻辑。第三个方向是加入水库弃水惩罚。当库容接近上限时即使有足够的发电能力也不得不多放水以避免风险这时把弃水流量放进目标里会得到一个更贴近实际运行的调度方案。以上这些扩展每一条都能单独写成一篇完整的文章。这里先把最核心的 NSGA-III 实现路径讲透后面再逐步展开。写这套代码让我最深的一点体会是算法框架从来不是难点难点在模型和约束的准确性。代码里的每个约束都对应一条物理规律漏掉一个算法给出的漂亮数字就可能是一堆不能执行的废解。所以不管你是做研究还是做工程第一步都应该是把目标、变量、约束的边界彻底想清楚再谈算法的事。后面如果你也打算跑一套自己的水火联合调度建议从最简单的两目标两水库版本起步跑通了再逐步加复杂度这条路最稳。