干这一行最绕不开的就是综合能源系统的调度问题电、热、气多种能源耦合在一起目标还不止一个——你要省钱还要减排两个目标往往还打架。传统做法是把两个目标加权成一个单目标去优化但权重怎么定不同场景下权重该不该变这些问题在实际工程里很难回答。所以我直接把目光投向多目标进化算法最后选定NSGA-II非支配排序遗传算法在Matlab里把整个综合能源优化调度模型跑通了。这篇文章就把整条链路拆开讲为什么综合能源调度适合用NSGA-II目标函数和约束条件怎么建模Matlab代码怎么组织才不乱以及我在调参和踩坑过程中积累的一些经验。内容偏工程向有仿真需求的电力专业研究生、做园区能源调度的工程师或者刚接触多目标优化的新手都能在这篇文章里找到可以“抄作业”的部分。1. 综合能源优化调度的问题拆解1.1 从“电热气”耦合说起先别急着写代码得先把问题本身想清楚。综合能源系统IES一般会包含微型燃气轮机MT、燃气锅炉GB、电锅炉EB、储能电池ESS、光伏PV、风电WT以及上级电网的购电通道。电、热、气三种能量流在系统内部互相转换、互相制约。举个例子微型燃气轮机烧天然气发电发完电之后的高温烟气还能通过余热回收装置制热——这就是“热电联产”的基本逻辑。但这个逻辑本身就自带冲突你若想让燃气轮机多发电以降低购电成本它的余热也会同步增加如果此时热负荷并不高多余的热量就只能废弃造成能源浪费。反过来如果你只盯着“满足热负荷”来调度又可能让燃气轮机在电价低谷期也满负荷运行经济性很差。电锅炉更是直接把两个系统彻底绑死功率给大了热满足了电费账单也跟着飞涨。储能装置则像一个时间维度上的调节器在电价低谷充电、高峰放电但也加剧了调度模型的状态耦合。所以综合能源调度本质上是一个“既要、又要、还要”的问题——要经济、要低碳、还要满足随时波动的电热负荷。单目标优化在这里说不通因为你根本找不到一个放之四海而皆准的权重比例。1.2 目标函数设计成本和碳排放我把目标函数拆成两个运行总成本和碳排放量。这里为简化模型把弃风弃光惩罚也折算进成本里碳排单独作为一个目标函数。第一个目标运行成本F1包括购电费用从上级电网买电的费用按分时电价计算注意峰谷平各时段电价差异很大燃料费用燃气轮机和燃气锅炉消耗天然气的费用燃气轮机燃料成本是其输出电功率的二次函数更贴近实际特性设备运维费用各设备单位出力对应的维护成本储能和燃气轮机这两块的运维单价通常偏高弃风弃光惩罚为保证结果不过度牺牲可再生能源消纳加入一个惩罚项数值上等于弃风弃光量乘以一个不小的惩罚系数。第二个目标碳排放量F2包括购电等效碳排放电网输入的每度电对应的平均排放因子这个因子在不同区域电网差异很大算的时候要留意用当地数据而非网上的通用值天然气燃烧排放燃气轮机和燃气锅炉烧天然气直接产生的二氧化碳。两个目标的数学形式展开后都是一个时段的量乘上调度周期步长取1小时再累加到总调度周期24小时上。写成代码就是一个函数输入是决策变量构成的向量输出是[F1, F2]这两个数。1.3 约束条件电、热、储能一个都不能少约束条件决定了模型的物理可实现性。我在模型里考虑了四类约束电功率平衡约束各个电源出力加储能充放电功率、再减去电锅炉消耗的功率必须等于电负荷热功率平衡约束燃气轮机余热回收、燃气锅炉产热、电锅炉产热三者的总和必须等于热负荷设备出力上下限约束每台设备都有出力范围不能越界储能SOC状态约束电池的荷电状态按时间递推受充放电功率上限和容量约束限制。对应到NSGA-II里面这些问题都归约到一个点上决策变量的编码和约束处理方式。我的做法是决策变量取每个时段各设备的出力值约束用惩罚函数法整合进目标函数里。具体怎么设计下一部分展开讲。2. NSGA-II核心机制与Matlab实现思路2.1 非支配排序和拥挤度算法为什么“能打”NSGA-II的核心机制有三板斧快速非支配排序、拥挤度距离计算、精英保留策略。快速非支配排序解决的是“谁比谁好”的问题。两个目标的情况下一个解如果在成本和碳排放上都优于另一个解那它就能支配后者。所有不被任何解支配的解排在第一前沿去掉它们之后剩下的解再重新判断支配关系得到第二前沿、第三前沿……这样一层层剥开。调度问题里的Pareto前沿就是第一前沿那些解构成的曲线。拥挤度距离解决的是“前沿上解太挤”的问题。光有排序还不够如果某种区域堆了一大堆放不下而另一边只有零星几个解那Pareto前沿就畸形了。拥挤度距离把某个解沿着每个目标方向上前后的两个解的欧氏距离算一遍再归一化求和距离越大说明周围越“空旷”在淘汰的时候优先保留。精英保留策略则是把父代和子代合并成一个2N大小的种群先按前沿层级排排完再按拥挤度排取前N个作为下一代。父代里表现好的个体不会被后代替换冲掉保证算法不会“退化”。这三个机制加在一起让NSGA-II在综合能源这种典型的多目标、非线性、带约束问题上表现很稳不需要额外调节参数太多就能得到一个分布性不错的解集。2.2 决策变量编码与种群数据结构设计在Matlab里实现NSGA-II最核心的一个设计决定是种群的数据结构。我建议直接用结构体数组struct array而不是用元胞数组或者classdef类。原因有三个结构体数组在修改单个字段时不会影响其他字段语义清晰排序、选择、淘汰的时候用一个字段如rank来辅助筛选非常方便访问pop(i).dec和pop(i).obj在代码里可读性极高后续给别人看代码也更容易理解。决策变量编码是这样的假设调度周期是24小时每个时段有6个可调度设备/通道——燃气轮机输出电功率、燃气锅炉产热量、电锅炉耗电量、储能充电功率、储能放电功率、电网购电量。那么一个个体即一个调度方案的决策向量长度就是6*24144。种群初始化的时候在这个144维空间里按照上下界随机撒点就是初始解集。这里有个细节因为储能充电和放电在同一个时段不应该同时为正我通过设置决策变量的逻辑关系来规避而不是额外加约束。具体做法是把“净充电功率”当作一个决策变量正值表示充电负值表示放电在计算目标函数时再把它拆成充电功率和放电功率并各自带上对应的效率系数。这样既保证了物理合理性又少了一个变量的维度。2.3 为什么不用权重法而要用多目标进化算法这其实是大部分初学者最容易问的问题“我把碳排目标乘以一个系数跟成本加在一起用粒子群或者单目标遗传算法直接优化不就行了”这样做在单次求解上没问题但实际工程里你会遇到两个麻烦。一是权重系数没有先验依据每个系统的最佳权重都不一样必须反复试。二是单次运行只能给你一个解你无法给决策者提供一个“成本最低但碳排略高”和“碳排最低但成本略高”的备选方案集。而用NSGA-II跑一轮直接送上来一整个Pareto前沿让决策者在上面挑。对于综合能源调度这种涉及经济部门和环保部门利益平衡的场景多解输出比单解输出实用得多。当然多目标进化算法的代价是计算量大。种群规模200迭代300轮目标函数每评估一次就要跑一遍24小时的潮流平衡计算整体耗时在普通PC上大概几十秒到一两分钟。这个量级对于日前调度这种“提前一天算一次、不追求实时”的场景完全可以接受。3. Matlab代码实现与关键参数配置3.1 环境与代码结构我用的是Matlab R2022b实测R2019b之后的版本跑这份代码没有任何问题。整套代码不需要额外装任何工具箱NSGA-II是自己纯手写的不依赖全局优化工具箱也不用YALMIP或CPLEX。代码按功能拆成以下几个文件结构上比较清晰main_NSGAII_IES.m主程序负责读取参数、初始化种群、调用进化循环、绘图initPop.m种群初始化解析上下界约束生成第一代解集nonDominatedSort.m快速非支配排序返回每个个体的rank和所在前沿编号crowdingDistance.m计算每个个体的拥挤度距离selection.m二元锦标赛选择算子sbxCrossover.m模拟二进制交叉polynomialMutation.m多项式变异environmentSelect.m精英保留环境选择合并父子种群再裁剪回NevaluateObj.m目标函数评估包含设备模型、约束判断和惩罚函数。主循环大约长这样for gen 1:maxGen % 生成子代 offspringPop []; while length(offspringPop) popSize p1 selection(pop, rank, crowdDist); p2 selection(pop, rank, crowdDist); [c1, c2] sbxCrossover(p1.dec, p2.dec, etaC, lb, ub); c1 polynomialMutation(c1, etaM, mutateProb, lb, ub); c2 polynomialMutation(c2, etaM, mutateProb, lb, ub); c1.obj evaluateObj(c1.dec, sysParams); c2.obj evaluateObj(c2.dec, sysParams); offspringPop [offspringPop, c1, c2]; end % 合并父代和子代 combinedPop [pop, offspringPop]; [rank, crowdDist] nonDominatedSortAndDistance(combinedPop); % 精英保留选择前popSize个 pop environmentSelect(combinedPop, rank, crowdDist, popSize); end3.2 几段核心函数的实现细节快速非支配排序是NSGA-II里最绕的一段逻辑但Matlab实现起来并不复杂。关键是两重循环里要同时统计“支配别人的解集合”和“被支配计数”function [rankSet, frontCell] nonDominatedSort(popObj) N size(popObj, 1); dominatedCount zeros(N, 1); dominateSet cell(N, 1); rank zeros(N, 1); frontCell {}; for i 1:N for j 1:N if i j, continue; end if all(popObj(i,:) popObj(j,:)) any(popObj(i,:) popObj(j,:)) dominateSet{i}(end1) j; elseif all(popObj(j,:) popObj(i,:)) any(popObj(j,:) popObj(i,:)) dominatedCount(i) dominatedCount(i) 1; end end end front find(dominatedCount 0); frontIdx 1; while ~isempty(front) frontCell{frontIdx} front; nextFront []; for k 1:length(front) i front(k); rank(i) frontIdx; for m 1:length(dominateSet{i}) j dominateSet{i}(m); dominatedCount(j) dominatedCount(j) - 1; if dominatedCount(j) 0 nextFront(end1) j; end end end front nextFront; frontIdx frontIdx 1; end end拥挤度距离计算比较直接核心就是按每个目标方向排序后用相邻个体的目标差除以当前目标的范围累积到原个体上。首尾两个点的拥挤度直接设成无穷大保证边缘个体被保留。function dist crowdingDistance(popObj, frontIdx) m size(popObj, 2); n length(frontIdx); dist zeros(1, n); for k 1:m [~, order] sort(popObj(frontIdx, k)); dist(order(1)) inf; dist(order(end)) inf; fmin popObj(frontIdx(order(1)), k); fmax popObj(frontIdx(order(end)), k); if fmax - fmin 1e-9 continue; end for j 2:n-1 dist(order(j)) dist(order(j)) ... (popObj(frontIdx(order(j1)), k) - popObj(frontIdx(order(j-1)), k)) / (fmax - fmin); end end end在这个函数里fmax-fmin的归一化必不可少。成本目标动辄几万元到几十万元碳排放目标可能只有几吨到几十吨量级差了好几个数量级。如果不归一化量纲大的目标会直接淹没另一个目标的信息。3.3 目标函数里的设备模型与约束处理目标函数是整个调度的灵魂也是最容易出错的地方。我给出一个简化的设备模型版本方便你按自己的系统替换燃气轮机发电效率取30%热电比取1.3意味着发1单位电的同时产生1.3单位的热量热量经过余热回收装置后有个0.9的回收效率燃气锅炉热效率取85%天然气低热值取9.7 kWh/m³气价取2.5元/m³电锅炉的电热转换效率取95%储能电池充放电效率分别取95%和92%容量500 kWhSOC限值 [0.1, 0.9]最大充放电功率100 kW电网购电电价按峰谷平三档峰段1.2元/kWh、平段0.75元/kWh、谷段0.4元/kWh电网排放因子取0.581 tCO2/MWh天然气排放因子取0.2 tCO2/MWh。每评估一个个体先根据决策变量算出每个时段的购电量、燃料消耗量然后算电功率平衡偏差和热功率平衡偏差再加上SOC越限的惩罚。惩罚项直接加到两个目标函数里用足够大的惩罚系数来保证进化过程会优先淘汰那些违反物理约束的个体。目标函数的骨架代码如下function [obj, balancePenalty] evaluateObj(dec, sysParams) Hs 24; P_MT dec(1:Hs); H_GB dec(Hs1:2*Hs); P_EB dec(2*Hs1:3*Hs); P_batt dec(3*Hs1:4*Hs); % 正值充电负值放电 P_grid dec(4*Hs1:5*Hs); % 电动率平衡 elecBalance P_MT P_grid sysParams.P_pv sysParams.P_wt ... - P_EB - P_batt - sysParams.P_load; % 热功率平衡 heatBalance sysParams.heatRatio * P_MT * sysParams.HR_eff ... H_GB sysParams.EB_eff * P_EB - sysParams.H_load; balancePenalty sum(max(abs(elecBalance)-1e-3, 0).^2 ... max(abs(heatBalance)-1e-3, 0).^2); % 成本计算和碳排放计算... obj [cost sysParams.penaltyCoef * balancePenalty, ... emission sysParams.penaltyCoef * balancePenalty]; end这里有个我在实际调试中踩过的坑惩罚系数太小进化算法发现违反约束还能混进下一代最终解在功率平衡上有明显偏差惩罚系数太大数值上又会把成本、碳排的真实差异“淹没”目标函数全被惩罚项主导Pareto前沿退化成一个点。我建议惩罚系数设置成设备最大成本量级的100倍以上并且把偏差平方而不是偏差本身作为惩罚项这样既能区分小偏差和大偏差又不会让罚函数过于“陡峭”。3.4 参数整定与运行配置NSGA-II本身参数不多但每个参数对结果影响不小。我给出实测比较稳妥的一组组合popSize 200种群太小容易早熟太大计算时间成倍增加maxGen 300多数算例在150代左右Pareto前沿基本收敛300代够安全etaC 20SBX交叉的分布指数越大生成的子代越接近父代etaM 20多项式变异的分布指数同上mutateProb 1 / nVar变异概率设置为决策变量维数的倒数实测这个经验公式比固定0.1效果好得多锦标赛选择的参赛个体数取2。跑一次24小时调度的耗时大概在40~90秒之间取决于具体的目标函数复杂度和电脑性能。建议初始试验阶段先跑maxGen 100验证代码没bug再上完整300代。4. 常见问题与排查技巧实录4.1 初始种群全是不可行解很多人第一次跑通代码后一看初始种群发现绝大多数个体在电功率平衡或热功率平衡上偏差很大Pareto前沿压根不存在。这种情况几乎都是决策变量初始化的取值范围设置问题。设备上下限设得太宽随机生成的解在平衡约束附近完全碰不上“可行域”。解决办法有两个一是上下限先设窄比如燃气轮机上下限就按 [0.3*rated, rated] 设给随机初始化一点引导二是在种群初始化里面直接加入一部分“启发式个体”——手动构造几个调度方案比如燃气轮机全天满发的方案、储能恒功率充放的方案、全部购电的方案把这些已知可行的方案混入初始种群。实测混入启发式个体后收敛速度会明显加快。4.2 Pareto前沿分布不均或者聚成一团一种是前沿整体堆在一小块区域。常见原因是拥挤度距离的归一化方向出了问题两个目标的量纲差距太大导致拥挤度计算形同虚设。检查crowdingDistance里是否对每个目标分别做了归一化没有的话补上。另一种情况是前沿出现“断带”——某一段区间几乎没有解。这通常是变异算子跨不过局部不可行区造成的可以适当增大etaM或者提高变异概率。还有一个思路把调度问题的约束做一些等价变换比如把电功率平衡里的松弛变量购电量显式表达成其他变量的函数从而从模型层面消掉一部分平衡等式。4.3 结果每次跑都不一样这是种群算法的天性不是bug。但如果你要在论文或工程报告里汇报结果必须保证可复现。方法很简单在主程序最前面固定随机种子rng(42);固定种子后在同一台机器上同一份代码跑10次结果完全一致。另外我建议项目里多试几个种子比如42、7、2024取多次运行的综合结果来画Pareto前沿而不是只依赖单次运行。原因是一个种子只能代表随机搜索的一次路径结果可能带有偶然性。4.4 从Pareto前沿里怎么选最终调度方案多目标优化的经典难题就是“解太多选哪个”。我处理这个问题用的是模糊隶属度方法。先把前沿上每个解的每个目标做归一化然后对每个目标定义一个满意度函数——成本越低满意度越高碳排越低满意度越高。每个解的总体满意度取两个目标满意度的最小值。最后选总体满意度最大的那个解作为折中方案。这个方法的好处是无偏、可解释性强。你也可以用TOPSIS法原理上类似。实际用下来模糊隶属度选出的方案在成本和碳排之间往往能取得一个肉眼可见的平衡单独画出来看也是一个符合工程直觉的调度结果。5. 实验结果分析与调度方案的工程启示5.1 典型仿真结果怎么看我用自己设的一组建模参数跑完300代之后把200个末端个体画成散点图。横轴是运行成本元/天纵轴是碳排放量吨/天。得到的Pareto前沿是一条自左上到右下倾斜的凸曲线。左下端代表“最环保但不便宜”的方案燃气轮机尽量少开全靠电网购电和燃气锅炉供热碳排低但购电成本高昂。右上端代表“最省钱但不绿色”的方案燃气轮机几乎满发利用气价相对电价便宜的优势压低成本但天然气燃烧带来的碳排放居高不下。折中方案落在曲线中段燃气轮机的出力维持在60%~80%额定功率区间储能以“谷充峰放”的方式运行。这个结果本身就说明综合能源调度的“绿色”和“经济”之间存在一个清晰的物理约束边界NSGA-II的作用就是把这整条边界找出来。5.2 储能对Pareto前沿的影响我做过一组对比实验其他参数不变把储能容量从0改成500 kWh再看Pareto前沿位置的变化。结果显示装储能之后整条前沿明显向右下移动——在同样的碳排放水平下成本可以降低约6%~10%在同样的成本预算下碳排放也能降低不少。这背后的逻辑很直观储能增加了系统的调度自由度让调节能力更强相当于把“多花钱减碳”和“少花钱多排碳”这两个极端之间的可调节范围扩大了。在工程上这意味着多目标优化不止是算法问题也是设备配置方案评估的一个有力工具。你完全可以把这个仿真代码改造成“容量配置优化”的雏形——把储能容量、燃气轮机台数也变成决策变量跑一次就能得到设备和运行策略的联合Pareto前沿。5.3 从仿真到工程落地的一些思考算出来的Pareto前沿漂亮不代表能直接搬到实际系统里。这里有几个工程问题仿真中常常被忽略一是设备模型参数的不确定性燃气轮机效率随环境温度变化、储能电池效率随SOC变化仿真里用的定值在工程里会产生偏差二是突发工况约束比如某条线路检修导致通道容量下降、极端天气下光伏出力骤降都意味着你只能在日前调度方案的基础上做滚动修正。有人问NSGA-II能不能做实时调度理论上可以但由于每轮进化都涉及种群评估计算开销较大现阶段更适合做日前调度、设备容量规划这类“离线决策”场景。实时调度一般用基于NSGA-II离线生成规则库再换用MPC等在线算法执行。但无论如何做多目标调度项目的最大价值不在于“算出最优解”而在于把复杂问题显式化——你把成本、碳排、设备限制、能量平衡全部翻译成了可计算的目标和约束然后机器帮你把整条权衡曲线都算出来。决策者要做的是在这条曲线上根据政策、市场和风险偏好选一个点。这个过程比任何“拍脑袋定权重”都科学得多。最后说个细节方面的经验仿真写代码时建议把evaluateObj里的设备参数集中放在一个结构体sysParams里不要散落在各个函数的输入参数里。我一开始图省事直接把数值硬编码在目标函数里结果换一组电价数据时改到怀疑人生。把这些参数提炼成一个结构体后换算例、做敏感性分析都只用改一处。看似只是代码整洁问题实际能省下大把调试时间。
