在综合能源系统优化调度这个方向上NSGA-II几乎是绕不开的名字。只要涉及多目标——比如既要省钱又要减排、既要响应快又要设备寿命长——单目标优化那套加权求和的路子就很难兼顾而NSGA-II这种基于非支配排序的多目标进化算法可以直接在目标空间里找出一组互不支配的Pareto解让决策者根据当天的电价、天气和负荷情况去挑最合适的调度方案。这篇文章我把我自己用Matlab实现NSGA-II求解综合能源优化调度问题的完整过程拆开讲从算法原理、系统建模到代码实现和调参避坑一次讲透。先说清楚这个东西解决什么问题。传统能源系统里电、热、气、冷是分开规划、分开调度的各管各的。综合能源系统的核心思路是打破这些壁垒把光伏、风电、微型燃气轮机、储能、电锅炉、吸收式制冷等设备耦合在一起通过协调各设备的出力让整个系统在满足电热冷负荷的前提下运行成本最低、碳排放最少、可再生能源消纳率最高。但这个协调过程非常复杂因为三个目标之间往往是冲突的——多烧气确实能补上出力缺口但成本和碳排放同时就上去了储能多放电能省钱但频繁充放会缩短电池寿命。这种多个矛盾目标同时优化的问题就是典型的多目标优化问题也是NSGA-II的主场。1. 为什么综合能源调度非要用NSGA-II1.1 多目标问题用单目标算法是强行凑合我刚接触综合能源优化调度那会儿最先想到的方案其实特别朴素把成本、碳排放、可再生能源利用率这三个目标加权重求和变成一个目标然后扔给遗传算法或者粒子群去跑。这么做确实简单但用起来会发现几个特别难受的问题。第一个问题是权重极其难定。成本权重调大一点算法就会牺牲减排目标碳排放量嗖嗖往上涨把减排权重调上去运行成本又压不下来。每次工况一变比如电价峰谷时段变化了、光照条件变了原来的权重组合就失效了得重新试。第二个问题是加权法本质上只是在凸的Pareto前沿上能找到较好解而综合能源系统这个优化问题目标函数和约束条件都是非线性的可行域也不一定是凸的加权法找出来的解很可能漏掉大量真正有意义的折中方案。NSGA-II这类多目标进化算法就不存在这个问题。它本质上是在同时优化多个目标函数而不是把多个目标揉成一个。算法最终输出的不是单个最优解而是一整组Pareto最优解集。比如某一天预测明天上午光伏出力很强、电价处于平价时段那Pareto前沿上碳排放最低的那个解会倾向于让光伏满发、储能充电燃气轮机停机而运行成本最低的那个解可能会在下午电价尖峰时段让储能放电、燃气轮机顶上去。决策者晚上看到这一组解能根据明天的实际生产计划、上级下达的减排指标或者售电价格预期挑一个最合适的调度策略执行整个决策过程灵活得多。1.2 NSGA-II的三大核心机制决定了它就是干这个的说实话多目标进化算法不止NSGA-II一种MOEA/D、SPEA2、PESA-II这些也各有特点。但NSGA-II在综合能源调度场景里能够成为事实标准靠的是它那三个设计得非常精巧的机制。第一个是非支配排序。种群里的每个个体在调度问题里就是一个完整的设备出力方案都会被拿出来两两比较。如果方案A在所有目标上都优于或等于方案B而且至少在一个目标上严格优于B那A就支配B。把所有不被任何其他方案支配的个体挑出来就是第一层Pareto前沿去掉它们之后剩下的个体里再挑一轮就是第二层以此类推。这样做的好处是不需要任何权重信息纯靠支配关系就能给群体分层每一层内部的方案在目标空间里都是谁也压不过谁的平级方案。第二个是拥挤度距离。分层解决的是解的优先级问题但光分层不够——第一层可能有几百个方案怎么从中选出下一代保留的个体如果随机保留Pareto前沿上的解会分布得非常不均匀某些区域挤成一团另一些区域大片空白。NSGA-II的做法是计算每个解在目标空间里与相邻解的欧氏距离距离越大说明这个解周围越空旷越应该被保留下来这样最终得到的Pareto前沿分布就非常均匀决策者能看到各个折中程度下的完整方案图谱。第三个是精英保留策略。在生成下一代种群时父代和子代合并成一个大种群先做非支配排序然后按照层级从低到高依次放入下一代直到某一层放不下为止这一层再按拥挤度从大到小挑选。这个设计的妙处在于优秀的父代个体不会被进化过程冲掉算法不会出现调着调着最优解丢了的尴尬情况。对综合能源调度这种计算耗时较长的场景每一步迭代的效率都要用在刀刃上精英保留策略保证了收敛性是有保障的。2. 先建模再写码综合能源系统怎么用数学描述2.1 典型系统架构与设备模型在动笔写NSGA-II代码之前最关键的一步其实是把综合能源系统的物理模型建好。这个模型的质量直接决定了优化结果的可靠程度。我参照一个比较典型的园区级综合能源系统来做介绍光伏和风电作为可再生能源电源微型燃气轮机作为可控电源电储能负责削峰填谷电锅炉和吸收式制冷机负责热负荷和冷负荷的供给同时系统可以与上级电网购售电。设备模型是优化的基础每个设备都需要用数学公式描述其输入输出关系。光伏出力的简化模型是P_pv(t) η_pv * S_pv * I(t)其中η_pv是光伏板光电转换效率S_pv是光伏板总面积I(t)是t时刻的太阳辐照强度。这是理想模型实际做的时候还要考虑温度修正系数、逆变器效率、光伏板老化衰减等。风电出力模型类似用的是风速的三次方关系P_wt(t) 0.5 * ρ * A * Cp * v(t)^3ρ是空气密度A是风轮扫掠面积Cp是风能利用系数v(t)是风速。这里的Cp在实际中不可能超过贝兹极限约0.593工程上取0.4左右比较合理。微型燃气轮机的模型要复杂一些。它的燃料成本与输出电功率呈非线性关系通常用二次函数拟合C_mt(t) a * P_mt(t)^2 b * P_mt(t) c同时它的热电联产特性决定了发电的同时会产出一定热量H_mt(t) P_mt(t) * η_hr / η_eη_hr是热回收效率η_e是发电效率。这个热电耦合关系是综合能源系统的精髓——电和热一起算才能体现联产的优势。电储能系统用SOC荷电状态描述SOC(t1) SOC(t) (P_ch(t) * η_ch - P_dis(t) / η_dis) * Δt / E_cap其中η_ch和η_dis分别是充放电效率E_cap是电池额定容量Δt是调度时段间隔。P_ch和P_dis不能同时大于零这组互斥约束在编码时需要特别处理。2.2 三个优化目标和它们的矛盾关系我选的三个优化目标很典型综合运行成本最低、碳排放量最小、可再生能源弃用率最低。这三个目标在数学上都写成全天24小时的累计值。综合运行成本C_total包括四部分向电网购电费用峰谷分时电价、燃气轮机的燃料成本、设备运行维护成本减去向电网售电的收入如果有余电。写成数学形式C_total Σ[ P_buy(t)*price_buy(t) - P_sell(t)*price_sell(t) C_fuel(t) C_om(t) ]注意这里的price_buy(t)是分时电价峰时段和平时段价格可以差出三四倍这对优化结果影响非常大。我实际测试过如果忽略分时电价储能系统基本不会在低谷充电、高峰放电整个调度方案的经济性就会差很多。碳排放量E_co2主要来自两方面外购电力的间接排放和燃气轮机的直接排放。E_co2 Σ[ P_buy(t)*EF_grid P_mt(t)*EF_gas ]EF_grid是电网排放因子不同地区的电网清洁程度不同取值差别很大北方火电占比高的地区这个值能到0.8以上南方水电风电占比高的地区可能只有0.4左右。在做区域碳排放核算时这个参数需要查当地电网的最新数据不能凭感觉填。可再生能源弃用率R_curtail定义为弃风弃光电量占可再生能源理论可发电量的比例目标是最小化这个比例也就是尽量把光伏风电全消纳掉R_curtail 1 - Σ(P_pv_used(t) P_wt_used(t)) / Σ(P_pv_avail(t) P_wt_avail(t))这三个目标的经济含义和物理含义完全不同互相之间还有内在冲突。最典型的就是降成本和减碳排之间的矛盾在夜间低谷电价时段燃气轮机发电的燃料成本可能比从电网买电还便宜但排碳一定比电网买电高假设电网里有一部分风光水核的清洁电反过来在光伏大发的时段为了消纳光伏、降低弃用率储能在低价时充电再在高价时放电虽然能赚峰谷价差但电池充放电的损耗和折旧会增加运维成本。这种矛盾关系用单目标加权法很容易被权重抹平而NSGA-II可以把这些矛盾完整地呈现在Pareto前沿上。2.3 约束条件的处理方式约束条件决定了优化解是否物理可行它的重要性甚至超过目标函数。我整理了一份综合能源调度中必须满足的约束清单约束类型数学表达物理含义处理策略电功率平衡P_grid P_pv P_wt P_mt P_dis P_load P_ch P_eb任意时刻发电与用电必须平衡等式约束用罚函数处理热功率平衡H_mt H_eb H_load H_ac热源输出必须满足热负荷等式约束用罚函数处理冷功率平衡Q_ac Q_load制冷量匹配冷负荷等式约束机组出力上下限P_mt_min ≤ P_mt ≤ P_mt_max设备物理出力范围编码时强制限幅爬坡约束-R_down ≤ P_mt(t) - P_mt(t-1) ≤ R_up燃气轮机出力不能突变可行性修复储能SOC范围SOC_min ≤ SOC(t) ≤ SOC_max电池不能过充过放编码时强制限幅储能充放电约束0 ≤ P_ch ≤ P_ch_max0 ≤ P_dis ≤ P_dis_max充放电功率限制编码时处理互斥关系购售电约束0 ≤ P_buy ≤ P_buy_max0 ≤ P_sell ≤ P_sell_max联络线功率限制编码时强制限幅弃风弃光约束0 ≤ P_pv_used ≤ P_pv_avail实际消纳不能超过可用值编码时强制限幅电力平衡这个等式约束是最核心的它在任何时刻都不能被突破。我在代码里对等式约束采用动态罚函数策略把不平衡量乘以一个自适应系数加到目标函数里。这个系数不能一成不变迭代初期种群多样性好、约束违反程度大罚函数系数可以小一些让个体有更多探索空间迭代后期种群逐渐收敛罚函数系数要逐步增大促使个体往可行域内收缩。爬坡约束的处理要更小心因为它天然是跨时段耦合的。如果t时刻的出力定下来了那t1时刻的可行范围就被限死了这在遗传算法的变异操作里很容易被忽略。我的做法是在生成初始种群时就按照爬坡约束生成时序曲线后续的变异操作只对局部时段进行小幅扰动这样能保证大部分个体天然满足爬坡约束。3. NSGA-II算法原理拆解三段核心机制讲透3.1 快速非支配排序非支配排序的思路其实特别直观就像在班里评三好学生先找那些所有科目都不比任何人差、而且至少有一科比所有人都强的同学他们就是完全碾压的第一梯队剩下的同学里再来一轮找出第二梯队以此类推。在代码里每个个体i需要维护两个关键信息S_i表示被个体i支配的所有个体的集合n_i表示支配个体i的个体数量。第一层Pareto前沿就是所有n_i0的个体。找到它们之后遍历第一层每个个体的S_i集合把集合里每个个体的n_i减1减完发现某些个体的n_i变成0了它们就是第二层。这个过程循环下去直到所有个体都被分层为止。我在Matlab里实现这一段的时候最需要注意的是两层循环的复杂度问题。种群规模50个个体、3个目标函数两两比较一次就是2500次支配判断每次判断要做3次目标值比较排序一轮下来的计算量在3000代迭代里是可观的。实测下来用纯for循环实现非支配排序50个种群跑3000代这部分大约占整个算法耗时的15%到20%。如果想提速可以按照目标函数值对个体先排序再做比较能减少不少无效比较次数不过代码复杂度会上去建议先把功能的正确性做出来再优化速度。3.2 拥挤度距离计算第一层Pareto前沿上的解虽然都是精英但精英内部也有区别。如果10个精英个体挤在目标空间的一个角落里剩下的大片区域没有解那这10个精英的代表性就差很多决策者看不到完整的折中图谱。拥挤度距离就是用来量化这种拥挤程度的指标。计算方法是先对同一层的所有个体按照某个目标函数值排序把两端无穷大或者一个很大的数赋给边界解保证它们一定会被保留因为Pareto前沿的端点也是很有价值的参考方案。中间每个解的拥挤度距离就是这个解在前后两个相邻解在各个目标方向上的距离之和。目标函数值归一化之后这个距离是有量纲的归一化后比较才有意义。我在刚开始实现拥挤度时踩过一个坑忘记对每个目标函数做归一化就直接算欧氏距离。结果发现成本动辄几万块、碳排放几千公斤、弃电率才百分之几三个数量级完全不对等成本一个变量就把其他两个变量完全压下去了拥挤度距离几乎完全由成本决定Pareto前沿分布极其糟糕。后来把所有目标值都归一化到[0,1]区间再计算拥挤度分布立刻均匀多了。3.3 精英保留与锦标赛选择精英保留策略是整个NSGA-II收敛性的保证。每一代的进化流程是这样的父代种群P规模为N通过锦标赛选择、交叉和变异生成子代种群Q规模也是N然后把P和Q合并成规模2N的大种群R。对R做非支配排序从第一层开始把每一层整层放进下一代种群S直到某一次放完之后S的规模即将超过N。对这个临界的最后一层按照拥挤度距离从大到小排序只取前面几个个体填满S让S的规模恰好等于N。这样得到的新种群S既保留了父代中最优秀的个体也容纳了子代新探索出来的方案不会出现最优解在进化中丢失的情况。锦标赛选择则是在选择父母时从种群中随机抽两个个体优先挑Pareto层级更低的层级越小越好第一层是最优的如果层级相同就挑拥挤度距离更大的。这种先比层级再比分布的两级比较策略就是NSGA-II在收敛性和多样性之间做平衡的核心机制。4. Matlab代码实现完整流程与关键代码段4.1 整体代码框架与参数配置我用Matlab实现这套算法时整体的模块划分是这样的主脚本负责初始化参数、调用进化循环、输出结果种群初始化模块负责生成满足约束的初始调度方案目标函数计算模块负责把调度方案换算成成本、碳排放和弃电率NSGA-II核心模块负责非支配排序、拥挤度计算、选择交叉变异结果可视化模块负责画出Pareto前沿和各设备出力曲线。算法参数方面我常用的配置是种群规模N50最大迭代次数maxgen3000交叉概率Pc0.9变异概率Pm0.1交叉分布指数eta_c20变异分布指数eta_m20。这里特别要提交叉分布指数和变异分布指数它们是模拟二进制交叉SBX和多项式变异中的关键参数控制着子代与父代的相似程度。eta值越大子代越接近父代局部搜索能力强eta值越小子代偏离父代越远全局探索能力更强。对综合能源调度这种变量多、约束复杂的问题我建议eta_c取20左右、eta_m取20到30兼顾局部精细搜索和全局探索。决策变量的编码方式需要说一下。一天的调度周期按1小时划分共24个时段每个决策变量都是长度为24的向量。我选择编码的决策变量包括燃气轮机每个时刻的出力P_mt(1:24)、储能每个时刻的放电功率P_dis(1:24)、储能每个时刻的充电功率P_ch(1:24)、电锅炉每个时刻的耗电功率P_eb(1:24)、向电网购电功率P_buy(1:24)。这样共5个24维向量每个个体是一个5×24的矩阵。光伏和风电的出力按预测值直接代入不作为决策变量因为它们是不可控的。%% 参数初始化 pop_size 50; % 种群规模 maxgen 3000; % 最大迭代次数 n_var 5 * 24; % 决策变量维度5个可控设备24时段 lb zeros(n_var, 1); % 变量下限 ub ones(n_var, 1); % 变量上限归一化实际映射时还原 pc 0.9; % 交叉概率 pm 0.1; % 变异概率 eta_c 20; % 交叉分布指数 eta_m 20; % 变异分布指数4.2 目标函数与约束罚函数实现目标函数是所有优化算法的灵魂这一步写不准确后面所有结果都是空中楼阁。我在计算目标函数之前先把归一化的个体解码成实际的功率值然后逐个时段计算成本和碳排放最后汇总成三个目标值。function [f, constraint_violation] evaluate_objective(x) % 解码x是归一化的决策变量矩阵大小为5×24(实际为扁平向量) P_mt decode_Pmt(x(1:24)); % 燃气轮机出力 P_dis decode_Pdis(x(25:48)); % 储能放电 P_ch decode_Pch(x(49:72)); % 储能充电 P_eb decode_Peb(x(73:96)); % 电锅炉消耗 P_buy decode_Pbuy(x(97:120)); % 购电功率 % 外购功率、光伏、风电、负荷数据由外部数据文件读取 global data P_pv data.P_pv; % 光伏预测出力 P_wt data.P_wt; % 风电预测出力 P_load data.P_load; % 电负荷 H_load data.H_load; % 热负荷 % --- 目标1综合运行成本 --- price_buy data.price_buy; % 分时购电价 C_grid sum(P_buy .* price_buy) * data.dt; C_fuel sum(data.a * P_mt.^2 data.b * P_mt data.c) * data.dt; C_om sum(data.k_mt * P_mt data.k_es * (P_dis P_ch) data.k_eb * P_eb) * data.dt; f_cost C_grid C_fuel C_om; % --- 目标2碳排放量 --- E_co2 sum(P_buy) * data.EF_grid * data.dt sum(P_mt) * data.EF_gas * data.dt; % --- 目标3可再生能源弃用率 --- % 电平衡约束发电购电放电 负荷充电电锅炉 P_balance P_buy P_pv P_wt P_mt P_dis - P_load - P_ch - P_eb; % 逻辑优先消纳可再生弃用功率等于因系统约束无法消纳的部分 P_curtail max(0, P_pv P_wt - (P_load P_ch P_eb - P_mt - P_dis - P_buy)); R_curtail sum(P_curtail) / sum(P_pv P_wt 1e-6); % 加小量防止除零 f [f_cost, E_co2, R_curtail]; % --- 约束违反量 --- g []; % 等式约束电功率平衡忽略小误差 g [g, abs(P_balance)]; % 12×1 % 热功率平衡约束 H_eb data.eta_eb * P_eb; % 电锅炉产热 H_balance data.H_mt_coef * P_mt H_eb - H_load; g [g, abs(H_balance)]; % 不等式约束转化为g0表示违反 SOC compute_SOC(P_ch, P_dis); g_soc max(0, SOC - data.SOC_max) max(0, data.SOC_min - SOC); g [g, g_soc]; % 爬坡约束 ramp_violation max(0, abs(diff(P_mt)) - data.ramp_max); g [g, ramp_violation]; constraint_violation sum(max(0, g)); % 总违反量 end这里有个细节要特别留意代码里电功率平衡P_balance是24维向量在罚函数设计时我是用这个向量的绝对值之和作为约束违反量的的一部分。但单纯加罚函数有个问题罚函数系数太小约束就松弛了系数太大了Pareto前沿会被压得变形。我后来采用的方法是先把罚函数项放在次级目标里也就是用约束支配的思路——如果一个解是可行解它无论如何都优于所有不可行解如果两个解都不可行比谁违反约束更少。这个方法比单纯加罚函数靠谱得多但实现复杂度也高一些。我在快速原型阶段用的是动态罚函数效果已经够用了。4.3 非支配排序与拥挤度计算核心代码非支配排序的代码是NSGA-II最关键的部分。下面这段我精简化了数据结构的实现重点展示逻辑流程function [rank, front] non_dominated_sort(fitness) % fitness: N×M矩阵N个体数量M目标数量 N size(fitness, 1); dominates false(N, N); % 两两比较建立支配关系矩阵 for i 1:N for j 1:N if i ~ j % 检查i是否支配j if all(fitness(i,:) fitness(j,:)) any(fitness(i,:) fitness(j,:)) dominates(i,j) true; end end end end % 计算每个个体被支配的数量和它支配的个体列表 n_dominated sum(dominates, 1); % 支配i的个体数 S cell(N, 1); % S{i}保存被i支配的个体 for i 1:N S{i} find(dominates(i,:)); end % 第一层没有被任何个体支配的个体 front {}; front{1} find(n_dominated 0); rank zeros(N, 1); rank(front{1}) 1; % 逐层剥离 k 1; while ~isempty(front{k}) Q []; for i front{k} for j S{i} n_dominated(j) n_dominated(j) - 1; if n_dominated(j) 0 Q(end1) j; rank(j) k 1; end end end k k 1; front{k} Q; end front(end) []; % 删掉最后一个空的层级 end拥挤度距离的计算逻辑function distance crowding_distance(fitness, front_indices) % front_indices是同一层内个体的索引列表 n length(front_indices); m size(fitness, 2); distance zeros(n, 1); for i 1:m % 按第i个目标函数排序 [~, idx] sort(fitness(front_indices, i)); sorted_idx front_indices(idx); % 边界个体距离设为无穷 distance(idx(1)) inf; distance(idx(end)) inf; f_min fitness(sorted_idx(1), i); f_max fitness(sorted_idx(end), i); if f_max - f_min 1e-10 continue; % 目标函数值完全相同跳过 end % 中间个体的拥挤度 for j 2:n-1 distance(idx(j)) distance(idx(j)) ... (fitness(sorted_idx(j1), i) - fitness(sorted_idx(j-1), i)) / (f_max - f_min); end end end4.4 模拟二进制交叉和多项式变异NSGA-II采用实数编码交叉算子用模拟二进制交叉SBX变异算子用多项式变异。这两种算子最大的特点是子代和父代的相似度可以通过分布指数来控制这对多目标优化里平衡收敛性和多样性非常关键。SBX的核心思想是让子代在父代附近生成而且生成的分布概率与分布指数eta_c相关。具体实现时先生成一个均匀分布的随机数u然后按照概率分布函数反解出展开因子beta当u0.5时beta (2u)^(1/(eta_c1))否则beta (1/(2(1-u)))^(1/(eta_c1))。子代的生成公式是child1 0.5 * ((1beta)*parent1 (1-beta)*parent2) child2 0.5 * ((1-beta)*parent1 (1beta)*parent2)多项式变异则是以一定概率对个体的某个决策变量进行扰动扰动幅度的控制由eta_m决定。eta_m越大扰动幅度越小局部搜索越精细eta_m越小扰动越剧烈越容易跳出局部最优。function child polynomial_mutation(child, pm, eta_m, lb, ub) % 多项式变异 for i 1:length(child) if rand pm u rand; if u 0.5 delta (2*u)^(1/(eta_m1)) - 1; else delta 1 - (2*(1-u))^(1/(eta_m1)); end % 限制变异后的值在边界内 if child(i) delta*(ub(i)-lb(i)) ub(i) child(i) ub(i); elseif child(i) delta*(ub(i)-lb(i)) lb(i) child(i) lb(i); else child(i) child(i) delta*(ub(i)-lb(i)); end end end end4.5 主进化循环主循环是整个算法的主干它把这些子模块串成一条流水线。%% 主程序NSGA-II进化循环 % 初始化种群 population init_population(pop_size, n_var, lb, ub); % 计算初始种群的目标函数值 fitness zeros(pop_size, 3); cv zeros(pop_size, 1); for i 1:pop_size [fitness(i,:), cv(i)] evaluate_objective(population(i,:)); end % 进化主循环 for gen 1:maxgen % 1. 锦标赛选择产生父代 parent_idx tournament_selection(fitness, cv, pop_size); parent_pop population(parent_idx, :); % 2. 模拟二进制交叉 offspring zeros(pop_size, n_var); for i 1:2:pop_size if rand pc i1 pop_size [offspring(i,:), offspring(i1,:)] sbx_crossover(... parent_pop(i,:), parent_pop(i1,:), eta_c, lb, ub); else offspring(i,:) parent_pop(i,:); if i1 pop_size offspring(i1,:) parent_pop(i1,:); end end end % 3. 多项式变异 for i 1:pop_size offspring(i,:) polynomial_mutation(offspring(i,:), pm, eta_m, lb, ub); end % 4. 评估子代个体 offspring_fitness zeros(pop_size, 3); offspring_cv zeros(pop_size, 1); for i 1:pop_size [offspring_fitness(i,:), offspring_cv(i)] evaluate_objective(offspring(i,:)); end % 5. 合并父代和子代进行选择和精英保留 combine_pop [population; offspring]; combine_fitness [fitness; offspring_fitness]; combine_cv [cv; offspring_cv]; new_population []; new_fitness []; new_cv []; % 约束处理可行解优先 feasible_mask combine_cv 1e-6; if sum(feasible_mask) pop_size % 可行解多于种群规模只用可行解做选择 feasible_pop combine_pop(feasible_mask, :); feasible_fit combine_fitness(feasible_mask, :); feasible_cv combine_cv(feasible_mask, :); [rank, fronts] non_dominated_sort(feasible_fit); % 逐层填充新种群 remaining pop_size; for k 1:length(fronts) if remaining 0 break; end front_size length(fronts{k}); if front_size remaining % 整层都保留 new_population [new_population; feasible_pop(fronts{k}, :)]; new_fitness [new_fitness; feasible_fit(fronts{k}, :)]; new_cv [new_cv; feasible_cv(fronts{k}, :)]; remaining remaining - front_size; else % 最后一层按拥挤度挑选 cd crowding_distance(feasible_fit, fronts{k}); [~, cd_idx] sort(cd, descend); pick fronts{k}(cd_idx(1:remaining)); new_population [new_population; feasible_pop(pick, :)]; new_fitness [new_fitness; feasible_fit(pick, :)]; new_cv [new_cv; feasible_cv(pick, :)]; remaining 0; end end else % 可行解不足需要保留部分不可行解辅助搜索 % 详细实现略思路是对不可行解单独按约束违反量排序 end population new_population; fitness new_fitness; cv new_cv; % 记录当前代的最优Pareto前沿信息 if mod(gen, 100) 0 fprintf(Generation %d: Pareto front size %d\n, gen, size(fitness, 1)); end end % 输出最终Pareto最优解集 final_front_idx fitness(~any(isnan(fitness),2), :); save(pareto_results.mat, population, fitness, cv);这里有个很关键的细节想强调进化过程中某一代可行解数量可能少于种群规模甚至为零。这时候如果强行只保留可行解种群会严重萎缩进化没法继续进行。我在处理这种情况时会对不可行解按约束违反量进行单独排序把违反量小的一部分个体保留下来。它们虽然不可行但携带了可能有用的遗传信息保留下来说不定哪一代经过交叉重组后就变成可行解了。5. 结果可视化Pareto前沿和设备出力曲线怎么看5.1 Pareto前沿图的解读方法跑完3000代之后算法会输出一组Pareto最优解。我习惯先把三维的Pareto前沿画出来X轴是运行成本Y轴是碳排放量Z轴是弃电率每个散点代表一个调度方案。这个图最直观的价值在于展示三个目标之间的权衡关系。我拿一组实际跑出来的数据举例Pareto前沿上有47个解成本和碳排放之间表现出明显的正相关趋势但弃电率这个目标和两者的关系更复杂。当系统允许一定程度的弃电时储能的调度可以更从容成本和碳排放都能降下来一些但如果追求极低的弃电率储能和燃气轮机就要在光伏大发时段频繁调整出力成本和碳排放都会上升。这就是典型的三目标Pareto前沿的形态。在实际项目中我会把Pareto前沿图发给用户让用户根据自己的偏好挑选方案。比如有些用户对碳排放有硬性指标要求比如全年碳排放不得超过某个上限那就在Pareto前沿上筛选出碳排放低于该阈值的解集再从中挑成本最低的那一个。这种先展示全貌、再根据偏好决策的工作方式比直接给一个最优解要实用得多。5.2 典型调度方案的设备出力曲线分析挑选一个中间偏经济性的Pareto解画出它的24小时设备出力曲线能看到非常明显的规律性。在凌晨0点到6点的低谷电价时段储能系统以最大功率充电把便宜的电存起来燃气轮机保持最低出力运行或者直接停机冷热负荷一部分由电锅炉的蓄热特性来满足。从早上8点开始进入平段电价光伏出力逐渐爬升燃气轮机开始提高出力补充电负荷缺口。到了晚上18点到22点的峰段电价储能系统开始全力放电尽量少从电网买高价电燃气轮机也进入高出力状态。这套规律其实对应的就是低谷充电、高峰放电、光伏优先、气电调峰的经典运行策略算法自己摸索出来的结果跟工程经验高度吻合这本身就是对模型的验证。弃风弃光比较严重的场景值得单独看一下。某天预测有台风过境风电出力在凌晨达到满发而夜间负荷低谷期用电量低储能满了电锅炉也开到上限了风电还是有富余。这时Pareto前沿上有一些解选择降负荷运行燃气轮机宁可让燃气轮机停机也要尽量消纳风电。另一些解则选择弃掉一部分风电但换来了更低的成本和更少的碳排放。这个时候决策者就要权衡——今天的风电是零碳的多消纳风电虽然有弃电率低的面子价值但如果系统因为过度消纳风电而让储能频繁深度充放、加速电池老化综合来看是不是划算这种精细权衡的依赖正是NSGA-II输出的Pareto前沿的价值所在。6. 常见问题与Matlab实操避坑指南6.1 Matlab版本与工具箱兼容性我在多个Matlab版本上跑过这套代码R2018b之前和之后的版本在处理结构体数组和数组索引时的行为有一些细微差别R2021a及之后的版本在JIT编译加速上对循环的优化更好。如果你的代码在我上面给出的非支配排序双层循环上运行太慢可以尝试把循环改成向量化运算。但有一说一Matlab在纯for循环上的执行效率确实不如Python如果你要做超大种群或者超长迭代的实验可以考虑把这部分核心逻辑改写成mex文件或者直接用Python的numba加速性能能提升一个数量级。另外建议全程使用Matlab的profiler工具分析代码性能。我第一次跑这套代码时发现evaluate_objective函数占了总运行时间的60%以上因为里面反复调用了compute_SOC函数而SOC的计算是通过for循环逐时段递推的。后来我把这个函数向量化运行时间从原来的3.5秒一次种群评估降到了1.2秒3000代的总耗时从接近4个小时缩短到了一个半小时左右。6.2 收敛性差和Pareto前沿分布不均匀如果你的Pareto前沿跑了很多代还是挤在一个小区域或者分布很稀疏优先检查这几个方面先看种群规模。50个个体在三维目标空间里其实很稀疏如果Pareto前沿形状复杂个体数建议加到80到100。但种群规模上去了每代的计算量也上去了这是一个需要平衡的点。再看交叉变异参数交叉概率太小会导致种群多样性不足变异概率太大会让算法退化成随机搜索建议从Pc0.9、Pm0.1的经典配置开始调试。最后看目标函数归一化。如果三个目标函数的数量级差异过大拥挤度距离会被量级大的目标主导Pareto前沿的均匀性会受影响。我之前遇到过成本和碳排放相差三个数量级的情况归一化之后效果立竿见影。6.3 约束处理不当导致的不可行解堆积这是我在综合能源调度中遇到的最常见问题。很多初学者把等式约束直接写成硬约束比如强行要求电功率平衡方程严格等于零然后发现可行解的占比极低算法根本收敛不了。因为等式约束加上爬坡约束、SOC递推约束之后可行域被压缩得很厉害随机产生的个体几乎不可能正好落在可行域内。我的建议是在初期采用宽松的罚函数策略让不可行解以一定的比例参与进化给算法一个渐进逼近可行域的过程。到了算法后期罚函数系数逐渐增大最终收敛到可行解。另外一种更工程化的做法是直接构造可行解先按爬坡约束随机生成燃气轮机的出力曲线再用功率平衡方程反推储能的充放电功率这样产生的个体天然满足主要等式约束大大缩小了搜索范围。代价是决策变量的自由度和搜索空间变小了可能会漏掉某些极端但最优的解。两种方案各有利弊我通常在快速原型阶段用第一种在最终精调阶段用第二种。6.4 数据输入与时间尺度匹配这个坑特别隐蔽但一旦踩中结果就完全不对。我们做的是24小时调度时间粒度通常是1小时但要注意所有数据的时间尺度必须统一。我遇到过有人把光伏出力数据的时间粒度设成15分钟负荷数据是1小时储能SOC递推公式里Δt混用结果优化结果完全失真。光伏、风电、负荷、电价、设备参数所有输入数据必须严格对齐到同一时间粒度上并且DateTime序列要一致这是一个在写代码之前就要确认好的事情。7. 一点实战经验分享做了几个综合能源优化调度的实际项目之后我最大的体会是算法本身只是工具真正的瓶颈在于系统建模的准确性和边界条件的合理性。NSGA-II这种算法的高明之处在于它不需要人为设定权重能自己在一轮轮进化中探索出目标之间的权衡关系但这也意味着它非常依赖输入数据的质量。设备效率参数、分时电价、光伏风电预测曲线任何一个数据失真Pareto前沿都会偏离现实。最后分享一个调参的小技巧先用小种群规模比如20个个体和少代500代快速跑几轮确认代码没有bug、约束条件合理、Pareto前沿形状大致合理再逐步加大种群规模和迭代次数。我在快速验证阶段一般会在Matlab里加一个实时曲线绘制窗口每50代刷新一次Pareto前沿图。如果看到前沿面在迭代中持续向坐标原点方向移动说明算法在收敛如果前沿面的形状一直剧烈变化就要回头检查是不是目标函数写错了或者约束条件前后矛盾。这套调试思路比闷头跑大实验再回头看结果高效得多。
