搞过综合能源优化的人都知道冷热电联供系统的运行优化看着简单真正把模型建起来、把算法调通细节多得能写一本书。最近我把这套基于多目标粒子群算法的冷热电联供综合能源系统运行优化代码重新整理了一遍从目标函数设计到Pareto前沿提取踩了不少坑也攒了不少经验。这篇文章就把完整思路、建模方法、算法原理和代码实现一次性讲透项目里需要做综合能源系统调度、园区冷热电联供方案对比或者准备把多目标粒子群算法落到实际工程场景的同学可以直接照着复现。1. 问题拆解与建模思路1.1 冷热电联供系统的物理架构先把系统架构说清楚。典型的冷热电联供综合能源系统核心设备包括燃气轮机或者燃气内燃机、余热锅炉、吸收式制冷机、电制冷机再加上蓄电池和蓄热/蓄冷储罐。能源输入端是天然气和电网购电输出端是电负荷、热负荷和冷负荷。燃气轮机烧天然气发电发电后的高温烟气进余热锅炉回收热能一部分直接供热另一部分驱动吸收式制冷机供冷。电负荷不够时从电网买电冷负荷高峰期电制冷机补充制冷能力储能环节用来平抑负荷波动。这套系统的耦合关系是优化建模的重点。电、热、冷三条能量流不是独立的燃气轮机发电越多余热越多可供的热和冷也越多但燃气消耗成本和碳排放同时上升。如果单纯追求经济性可能倾向于多用燃气轮机多发电阻止购电如果追求环保性又要尽量减少天然气消耗。两条目标往两个方向拉扯单目标优化根本压不住这种矛盾必须上多目标。我自己做建模时习惯用能量枢纽的方式来表达这套系统。把所有设备看成能量转换节点输入侧是电网和天然气输出侧是电、热、冷三类负荷中间每一台设备都是一个转换环节。这样写成的模型不依赖具体工程场景换一个设备配置只需要改参数表不用重写约束方程后面扩展到含氢储能、含电动汽车充放电的复杂系统也方便。1.2 设备模型与能量耦合关系设备模型是优化的物理基础每个设备我都是用“输入能量乘以转换效率等于输出能量”这个基本思路建模的。燃气轮机的模型稍微特殊一点除了发电效率还得关注热电比——也就是发电的同时产生多少可利用余热。我用的是定热电比模型实际工程中燃气轮机的热电比会随负荷率变化定热电比模型在部分负荷工况下有一定误差但作为运行优化层面的研究够用了。想做得更精细可以用变热电比曲线不过那样非线性约束会增加不少粒子群求解难度也会上升。余热锅炉的作用是把燃气轮机的排烟余热转化成供热蒸汽或热水模型相对简单就是一个输入余热量乘以回收效率。吸收式制冷机用余热驱动制冷量和输入热量的比值就是制冷系数。电制冷机就是一个耗电设备制冷量除以耗电量得到能效比。储能设备我统一用荷电状态来建模蓄电池用SOC表示蓄热和蓄冷罐用储罐容量比例表示充放能时有上下限约束还要注意储能不能同时充和放这个约束不少初学者会漏掉。整个优化过程是典型的多时间尺度调度问题我采用的是日前调度模式把一天24小时作为优化周期每个小时一个决策点一共24个时段。每个时段要决策的是燃气轮机出力、电网购电量、吸收式制冷机供冷量、电制冷机供冷量、储能设备充放功率。决策变量个数大概是5乘以24等于120个多目标粒子群算法跑这种规模的优化问题完全没有压力。1.3 约束条件的工程细节约束条件是模型里最容易出问题的地方。我梳理了一下这个系统需要满足的约束大致有四类功率平衡约束、设备出力上下限约束、爬坡约束和储能约束。功率平衡约束是整个模型的主心骨。电功率平衡要求每一时刻的燃气轮机发电、电网购电、蓄电池放电的总和等于电负荷、电制冷机耗电、蓄电池充电之和。热功率平衡是余热供热量加蓄热放热等于热负荷加蓄热充热。冷功率平衡是吸收式制冷加电制冷加蓄冷放冷等于冷负荷加蓄冷充冷。这三条约束每个时段都要严格满足如果计算结果中有某一时刻功率失衡算法迭代再多轮也没意义。设备出力上下限约束比较直接每台设备都有最大最小出力边界。爬坡约束是限制设备出力在相邻时段的变化速率主要是燃气轮机的爬坡速率限制电网购电的爬坡我一般也加上避免优化结果出现相邻时段购电量剧烈波动的不现实方案。储能约束包括充放能功率限制和容量边界蓄电池、蓄热罐、蓄冷罐的SOC都要维持在一个合理区间内。这里有个经常被忽略的问题储能设备的充放能会有损耗我一般用充放效率分别建模充电效率和放电效率取0.9左右这样更贴近实际。2. 目标函数设计与多目标博弈2.1 经济性目标运行成本怎么算才合理经济性目标我定义为系统日运行总成本包含三个部分购电成本、购气成本和设备运维成本。购电成本用分时电价乘以每小时购电量来算分时电价要分峰、平、谷三个时段设定峰时段价格大约是谷时段的3倍左右这个结构会直接影响优化结果中电网购电量的时间分布。购气成本是天然气单价乘以每小时燃气消耗量燃气消耗量和燃气轮机发电量的关系由发电效率决定公式是燃气消耗量等于发电量除以发电效率再除以天然气低位热值再乘一个单位换算系数。设备运维成本按照各设备出力乘以对应单位运维系数来算。运维系数不要拍脑袋定参考实际工程数据一般是燃气轮机每千瓦时几分钱储能设备因为充放电循环寿命损耗会高一些。把三个部分加起来得到总运行成本作为第一个目标函数。这里要注意量纲和数量级的匹配购电和购气成本是万元级别运维成本可能是几千元级别如果目标函数之间数量级差太大后面多目标优化的Pareto前沿会畸形最好做归一化或者在目标函数里乘以合适的系数把数量级拉齐。2.2 环保性目标碳排放折算是关键环保目标我选择了系统运行过程中的二氧化碳排放量。碳排放来源主要是两个电网购电对应的间接排放和天然气燃烧的直接排放。电网购电的碳排放用区域电网平均碳排放因子来折算天然气燃烧的排放量用单位热值含碳量乘以燃烧效率来算。这两种排放源的重要性在不同场景下权重差异很大如果电网排放因子高优化结果会倾向于多用天然气多发电阻止高碳电力的购入如果天然气价格高而排放因子低结果又不一样。这种目标间的博弈正是多目标优化的价值所在。除了碳排放指标有的项目会在第三个目标里加入一次能源利用率或者火用效率我这次没把它设成独立目标而是放在约束条件里要求系统一次能源利用率不低于某个阈值。这样做的考虑是目标函数数量增加会显著加大Pareto前沿求解难度两个目标时前沿是一条清晰曲线三个目标时前沿变成曲面粒子群的收敛性和解的分布性都不太容易控制。初学者我建议先从双目标做起两个目标的方案做透了再加第三个目标。2.3 为什么非要多目标单目标加权不香吗这是我最常被问到的问题。既然有经济和环保两个目标把两个目标加权求和转成单目标再做粒子群优化不就行了技术上确实可行但有个重要缺陷权重系数很难定得合理。经济目标和环保目标量纲不同、数量级不同权重系数的物理意义模糊实际操作中需要反复试凑而且一组权重只能得到一个解。真实工程中决策者希望看到的是经济性和环保性不同权衡下的多个候选方案然后在方案之间做决策多目标算法一次运行可以输出一整个Pareto前沿这种“方案套餐”的价值是单目标加权给不了的。从算法角度说多目标粒子群算法和单目标粒子群算法的核心结构也不同多目标版要有非支配排序、外部档案维护、leader选取策略这些额外机制。这些机制正是让Pareto前沿保持完整性和均匀性的关键后面代码实现部分会详细说。3. 多目标粒子群算法核心原理与参数选择3.1 从标准PSO到多目标多了哪些东西标准粒子群算法的逻辑很简单一群粒子在解空间里飞行每个粒子记住自己的历史最优位置同时跟随群体最优位置速度更新公式是经典的速度等于惯性项加个体认知项加社会认知项。这个机制在单目标问题里很直接因为“最优”就是那个函数值最小的点。到了多目标里“最优”不再是单一概念两个解A和BA的经济性好但环保性差B的环保性好但经济性差到底谁更优这时就要引入Pareto支配的概念。我用自己的理解来解释Pareto支配解A支配解B意味着A在所有目标上都优于或者等于B并且至少有一个目标严格优于B。如果A和B互不支配它们就属于同一个非支配层级。所有不被其他任何解支配的解组成第一层非支配解集也就是Pareto前沿。多目标粒子群算法的目标就是把粒子群推向这个前沿同时保证前沿上的解分布尽量均匀。所以多目标粒子群算法相比标准粒子群增加了三个关键模块非支配排序用来给所有解分层外部档案用来存储当前找到的Pareto前沿解leader选取策略决定粒子跟随哪个Pareto解飞行。这三个模块的优劣直接决定算法效果我调代码的经验是外部档案的维护和leader选取策略对最终前沿质量的影响比惯性权重和学习因子参数的调整影响更大。3.2 外部档案与网格密度法外部档案是算法的“成果留存区”粒子每次迭代找到的非支配解都会尝试存入外部档案如果新解支配档案里的某个旧解旧解要被删除。为了防止档案无限膨胀一般要给档案设置容量上限比如50到100个解。档案满了以后新解怎么处理就需要一个密度估计机制——在解的密集区域删除一些离得近的解保留稀疏区域的解让前沿分布更均匀。密度估计我推荐用自适应网格法比拥挤距离法实现简单效果好。思路是把目标空间划分成若干网格统计每个网格里解的个数解所在的网格里解越多说明这个区域越密集。新增解时如果它落在稀疏网格优先加入删除时从最密集网格里删。这种方法计算量小而且网格的边界会根据当前解集的范围动态调整不会因为解的分布偏移导致网格失效。这个网格密度机制直接用于leader选取也有不错的效果粒子从网格中选择一个含解较少的网格再从这个网格里随机选一个解作为自己的领导粒子。这样引导粒子朝前沿的稀疏区域飞避免了所有粒子挤在一起导致前沿分布不均的问题。3.3 参数选择的实践经验多目标粒子群算法的参数比标准PSO更多一些我把自己调好后一直沿用的参数写在这里供参考参数取值说明种群规模80太大计算慢太小前沿不完整最大迭代次数150复杂场景可以加到200惯性权重w0.9线性降到0.4前期探索、后期收敛个体学习因子c11.5保持个体历史经验社会学习因子c21.5跟随leader的力度外部档案容量60与种群规模匹配网格划分数10×10目标数为2时的常用配置惯性权重线性递减是我强烈推荐的做法。前期w大粒子速度快搜索范围广有利于探索整个解空间后期w小粒子速度减慢步长变小有利于在Pareto前沿附近精细搜索。我一开始用的是固定w0.8结果收敛速度慢差不多的迭代次数下前沿质量明显不如线性递减的方案。c1和c2我习惯取相等的值让粒子的“个人经验”和“群体经验”权重保持平衡需要特殊说明的是这里的群体经验在多目标里指的不是整个种群的最优而是每个粒子各自选定的leader。有一个调参细节值得说输出Pareto前沿的时候从外部档案里取解不直接取种群的解。因为外部档案经过密度维护解的分布性和支配性都比种群里的解更干净。种群每迭代一次外部档案都同步更新到迭代结束时候选集的质量已经足够好了。4. 代码实现与关键步骤解析4.1 主程序架构与初始化我用的MATLAB版本是通用的R2020a以上的功能整个优化程序分成四个模块主程序文件、目标函数文件、约束处理文件和MOPSO算法核心文件。主程序负责加载数据、调用算法、输出结果。目标函数文件接收一个决策变量向量返回经济成本和碳排放两个目标值。约束处理文件负责判断一个解是否满足所有约束。MOPSO核心文件就是算法本体包含粒子更新、非支配排序、外部档案维护三块逻辑。初始化部分要做的第一件事是配置设备参数和负荷数据。设备参数包括燃气轮机容量和效率、余热锅炉效率、吸收式制冷机COP、电制冷机EER、储能容量和充放功率上限等。负荷数据是24小时的电、热、冷负荷曲线我一般从典型日的实测数据或者模拟数据获取。然后是电价、气价、碳排放因子的参数设定这些数据对优化结果的影响非常直接。粒子的初始化有个容易被忽略的细节初代粒子的位置必须满足约束否则后续更新会一直在不可行域里打转。我的办法是采用“约束修正法”来处理随机生成初始解后检查约束如果不满足就对决策变量在边界范围内重新采样直到满足约束为止。初始化之后做非支配排序把初始Pareto前沿存入外部档案然后开始迭代。4.2 粒子更新与Pareto前沿提取标准粒子群速度更新公式我就不重复推导了直接给多目标版的核心逻辑。每个粒子迭代时按照更新公式计算新速度和位置但leader不是全局最优而是根据网格密度从外部档案中挑选出来的稀疏区域解。为了让种群在前沿分布均匀性上更好我还加了一个小概率变异操作以0.1的概率对粒子的某个决策变量做随机扰动这个操作在典型测试函数上能有效防止种群早熟收敛到前沿的局部区域。来看主循环代码框架for iter 1:maxIter w wMax - (wMax - wMin) * iter / maxIter; for i 1:popSize % 从外部档案中选择稀疏区域的解作为leader leader selectLeader(archive, gridInfo); % 速度更新 v(i,:) w * v(i,:) ... c1 * rand * (pbest(i,:) - x(i,:)) ... c2 * rand * (leader - x(i,:)); % 边界处理 v(i,:) max(min(v(i,:), vMax), vMin); x(i,:) x(i,:) v(i,:); % 边界约束修正 x(i,:) boundFix(x(i,:)); % 处理等式约束用分摊法或惩罚法使功率平衡成立 x(i,:) balanceFix(x(i,:), loadData); % 计算新解的两个目标值 [cost, emission] objFunction(x(i,:), sysParam); newObj [cost, emission]; % 更新个体历史最优 if dominates(newObj, pbestObj(i,:)) ... || (isequal(newObj, pbestObj(i,:)) rand 0.5) pbest(i,:) x(i,:); pbestObj(i,:) newObj; end end % 合并种群与外部档案做非支配排序更新外部档案 [archive, gridInfo] updateArchive(archive, x, objValues, gridNum); end这段框架基本就是MOPSO算法的核心骨架具体工程代码里还要加状态向量记录、每代Pareto前沿的保存逻辑、绘图函数的调用等内容。实际跑的时候建议在外部档案更新函数里加一个计数器记录每代被成功加入到档案中的新解数量。如果连续几十代这个数量都是0说明算法已经收敛了可以提前终止迭代省时间。4.3 折中解选取模糊隶属度法算法跑完得到的是一个Pareto前沿上面有几十个非支配解但工程上最终只需要一个运行方案。怎么从几十个解里挑一个最合适的我用的方法是模糊隶属度法。对Pareto前沿上的每个解先分别计算它在每个目标上的隶属度隶属度函数对成本类目标取越小越好某个解的某个目标隶属度越接近1说明这个目标越优。然后把每个解的所有目标隶属度求和作为这个解的满意度指标取满意度最大的解作为折中最优解。实现这个逻辑的代码很简短% objSet是Pareto前沿上所有解的目标矩阵每行一个解每列一个目标 objMin min(objSet, [], 1); objMax max(objSet, [], 1); numObj size(objSet, 2); numSol size(objSet, 1); satisfaction zeros(numSol, 1); for i 1:numSol muSum 0; for j 1:numObj if objMax(j) objMin(j) mu 1; else mu (objMax(j) - objSet(i,j)) / (objMax(j) - objMin(j)); end muSum muSum mu; end satisfaction(i) muSum; end [~, bestIdx] max(satisfaction); bestSolution repmatSol(bestIdx, :);模糊隶属度法相当于一种无偏好的决策方式它不偏向经济也不偏向环保而是选一个综合满意度最高的解。如果决策者明确倾向于经济性可以给经济目标的隶属度乘以权重因子再求和这样折中解会朝经济性更好的方向偏移。5. 仿真结果与典型场景分析5.1 算例配置与Pareto前沿形态我用某园区的典型日负荷数据做了仿真电负荷峰值约4000kW热负荷峰值约2500kW冷负荷峰值约1500kW负荷曲线有明显峰谷差。设备参数上燃气轮机额定发电功率设为2000kW发电效率0.35余热回收效率0.5吸收式制冷机COP为1.2电制冷机EER为3.5蓄电池容量2000kWh蓄热罐容量500kWh分时电价峰时段1.2元/kWh、平时段0.7元/kWh、谷时段0.35元/kWh。跑完多目标粒子群算法后Pareto前沿呈现出一条平滑下降的形状横轴是运行成本纵轴是碳排放量。前沿最左端解运行成本最低对应的排放量最高最右端解排放量最低但成本明显上升。整条前沿展示了经济和环保两个目标之间的完整权衡关系这是单目标优化完全无法输出的信息。比较有意思的是折中解对应的调度方案。折中解下燃气轮机在午间电价较高时段保持较高出力谷时段降低出力、优先用电网购电蓄电和蓄冷主要从谷时段工作。储能设备的作用很明显——在谷电价时把电能储存起来峰时段释放既压低成本又减少峰时段从电网购电的间接排放。热负荷主要由余热回收提供冷负荷优先由吸收式制冷机承担电制冷机只在冷负荷尖峰时补充出力。5.2 不同场景下的调度策略对比我把天然气价格从2.2元/立方米调到3.5元/立方米再跑一遍Pareto前沿明显右移所有解的成本都上升了同时前沿上解的分布也发生变化高气价下优化方案更倾向于减少燃气轮机出力、增加购电。如果此时谷电价时段充分蓄电、峰时段用储能顶峰整体系统运行成本仍然可控说明储能设备对能源价格的波动有缓冲作用。电价结构的影响同样显著。当峰谷电价差拉大时储能的经济价值更加凸显折中解中蓄电池在谷时段充电量明显增加。如果峰谷价差缩小到一定程度储能充放电循环的损耗成本可能高于电价套利收益优化结果中储能几乎不动作。所以在做方案报告时分时电价结构的设定对结论影响很大至少要跑两三组不同的价差条件给决策者参考。5.3 结果可视化与工程落地建议Pareto前沿的可视化用散点图即可将I代到最终代的非支配解集叠加在同一张图上可以清晰看到解的分布逐渐从分散走向前沿面收敛过程一目了然。折中解的功率平衡情况用堆叠面积图展示可以直观看到每个小时的电力来源构成和负荷去向构成。热、冷系统的能量流动建议用桑基图表现能讲清楚每一股能量的来源与去向。做工程报告时我还会补充一个“方案对比表”把经济最优解、折中最优解、环保最优解对应的运行成本、碳排放量和主要设备的出力结果放在一起横向对比。这张表能帮助非技术背景的决策者快速理解不同方案的价值差异比放一堆公式直观得多。曾经有个项目汇报就靠这张表促成了方案的最终确定实际效果比放十张曲线图都好。6. 常见问题与排查技巧实录6.1 粒子越界与迭代发散迭代过程里最常遇到的问题是粒子位置越界。速度更新公式中有随机数参与粒子一步可能飞出决策空间的边界如果不处理后面计算目标函数时会得到异常结果甚至报错。我处理边界的策略是速度边界和位置边界双重限制速度超过上限就截断到上限位置超出上下界就拉回到边界上。需要注意的是拉回边界之后速度也要相应调整可以令越界处的速度反方向这样粒子会被弹回可行域内部不会贴死在边界上无法动弹。还有一类问题是目标函数值出现NaN导致整个算法崩溃根源一般是约束处理中有除以接近于零的中间量。比如某个时段热负荷刚好很小吸收式制冷机功率下限约束和平衡约束会产生冲突算出来的中间流量接近零除零错误就来了。解决办法是在所有除法运算前加一个极小值保护abs(x) 1e-6时直接令分母等于1e-6。这个小改动解决了我至少十次仿真中的崩溃问题。6.2 前沿分布不均与收敛过快Pareto前沿分布不均是最常见的优化质量问题表现为前沿两端稀疏、中间密集或者干脆解全都堆在一个角落。我遇到的情况多数是leader选取策略没效果粒子都跟着同一个leader飞。检查外部档案的网格参数如果网格划分数太少密度估计分辨率不够稀疏区域和密集区域区分不出来leader选取就退化成随机选择。我把网格划分数从5×5调成10×10后前沿分布明显改善。收敛过快的问题也很好识别迭代到二三十代时外部档案几乎不再更新前沿的形状也固定了。这种情况通常是惯性权重下降太快或者变异概率太小。把线性递减的终值从0.3提高到0.45同时把变异概率从0.05提高到0.15粒子在后期的探索能力就恢复了。我在多组算例上验证过调整这两个参数后前沿质量显著提升。6.3 结果不合理时的排查路径仿真结果偶尔会出现一些“反常识”的调度方案比如明明电价很高优化结果却在大幅购电。遇到这种情况不要急着改算法先排查模型逻辑。我会按这个顺序排查先看功率平衡约束的代码确认每个时段的电平衡确实满足再看储能设备的充放能方向约束看蓄电池是否被允许同时充放电——如果有这种物理上不可能的状态优化结果就会变得夸张最后看目标函数代码确认成本计算方向没有写反比如有的版本会把碳排放量写成负数结果环保目标变成了排放越多越好。另外提醒一个量纲坑有时候仿真结果看起来“巧合”地漂亮其实是目标函数里漏了单位换算导致碳排放量和成本在数值上出现了不真实的同向变化。我在目标函数文件开头用注释把每个变量的单位标注清楚迭代过程中定期把某个已知方案的中间量打印出来用手工核算验证。这一步可以省去很多无效的调试时间。6.4 运行太慢的优化措施仿真程序跑一轮要十几分钟如果要做参数敏感性分析可能要跑几十上百轮速度太慢就会影响效率。我采用两种手段来提速。第一种是向量化优化把粒子群循环内部能矩阵化的部分尽量矩阵化尤其是目标函数计算和边界处理用向量化写法可以把单轮迭代时间压缩到原来的三分之一。第二种是提前终止机制设置外部档案连续多代不更新时停止迭代我设的是连续15代不更新就终止。实际测试中大部分算例在80到100代时就已经收敛强制跑满150代浪费了不少时间提前终止机制把整体耗时大约压缩了40%。7. 一些实操后的体会整套程序从模型搭建到算法调通我前前后后改了很多版本。印象最深的一版修改是把单目标的加权粒子群全部推翻重写成多目标结构那次重写之后很多以前讲不通的问题突然就清楚了。加权单目标得到的只是一个方案工程师做完仿真往往说不清为什么选择这个解——决策者问一句“环保再苛刻一点会多花多少钱”单目标解法完全答不上来多目标前沿可以很自然地回答这个问题。最后分享一个我一直在坚持的小习惯每次跑完仿真把目标函数和约束相关的参数集中放在一个文件里并写清楚版本号和数据来源。项目做久了参数会经常变动没有版本管理的话同一个模型跑出来的结果可能完全不同而你已经说不清是哪版参数了。这个习惯帮我避免了好几次返工。这套冷热电联供运行优化的思路不仅适用于以燃气轮机为核心的系统。微电网、园区综合能源、含电动汽车充放电的楼宇能源系统本质上都是多能源耦合下的多目标权衡问题把设备模型换一换、目标函数改一改同样的算法框架完全可以平移过去。后续我打算在这个基础上加入不确定场景的随机优化让模型对负荷预测偏差更有鲁棒性到时候再和大家分享。