做综合能源系统运行优化这项工作我接触最多的一条技术线就是冷热电联供CCHP。它把天然气发电后的余热回收起来再通过吸收式制冷和热交换把冷、热、电一块供出去能量梯级利用效率比单纯发电高不少。但真正落地的时候“怎么在一个24小时调度周期里分配各设备出力”这件事特别难既要省钱又要省一次能源还想压低碳排放这三个目标往往互相顶牛。多目标粒子群优化算法MOPSO是我在这个场景里用得最多、也最顺手的一套算法工具。它基于粒子群搜索框架叠加上Pareto支配关系和外部档案机制可以在一次运行里求出整条Pareto前沿给决策者留出充分的选择空间。这篇文章我把从系统建模、目标函数搭建、约束处理到Matlab里的算法主体实现再到算例结果取舍的完整流程全部拆开讲。适合刚接触CCHP优化、想真正把MOPSO用在Matlab仿真里的朋友参考。1. 冷热电联供系统到底长什么样1.1 设备构成与能量流典型CCHP系统包含燃气轮机、余热锅炉、吸收式制冷机、电制冷机、燃气锅炉、蓄电装置、蓄热装置这几个主要环节。天然气进入燃气轮机一部分能量变成电另一部分变成高温烟气和余热余热通过余热锅炉转化为蒸汽或热水供给吸收式制冷机制冷或者直接满足热负荷。当燃气轮机的余热不够用的时候燃气锅炉补上热量缺口。电力侧由燃气轮机发电和电网购电共同供给多余的电可以存入蓄电装置热侧如果产热大于负荷可以把热量存进蓄热罐。冷负荷由吸收式制冷机和电制冷机共同承担前者用的是热驱动后者用的是电驱动。这个结构决定了运行优化的核心逻辑设备之间能量耦合关系比较强某一台设备出力变化会同时牵扯到电平衡、热平衡、冷平衡所以不能把三套系统分开单独算必须整体建模。这也是CCHP运行优化比单纯微电网更麻烦的一个原因。1.2 运行优化是解什么数学问题以某个典型日为调度周期将一天分成若干时段惯例上按小时取24个时段。给定各时段的电负荷、热负荷、冷负荷以及天然气价格、分时电价、设备额定容量求解一组调度方案使多个优化目标达到平衡。决策变量一般是各时段每台可调设备的出力。我常用的变量清单包括燃气轮机发电功率电制冷机输入功率燃气锅炉供热量蓄电充电和放电功率蓄热充热和放热功率如果考虑设备启停状态还要加入0/1变量问题就从连续优化升级成了混合整数优化复杂度上了一个台阶。约束条件要覆盖的东西也不少。首先是三个能量平衡约束电平衡要求燃气轮机发电加电网购电加蓄电放电等于电负荷加电制冷耗电加蓄电充电热平衡要求余热回收加燃气锅炉供热加蓄热放热等于热负荷加吸收式制冷耗热加蓄热充热冷平衡要求吸收式制冷加电制冷等于冷负荷。其次是设备出力上下限、爬坡速率限制、蓄电蓄热的容量限制以及充放功率限制。2. 为什么是MOPSO而不是别的算法2.1 多目标问题的本质经济成本最小、一次能源消耗最小、碳排放最小这三个目标不是简单加权就能统一的。省钱的方案往往用更多低谷电但发电侧碳排放不低低碳方案倾向于让燃气轮机多发、用吸收式制冷结果购电少了、气耗高了成本又上去了。传统的加权法需要事先给定权重系数把多目标打成单目标。问题是权重怎么定本身就是主观的而且一次只能求一个解。想看清楚“多花多少钱能换多少减排量”就得反复改权重重新跑很多遍效率很低Pareto前沿非凸的时候还可能把一部分有效解漏掉。多目标优化换了一个思路不求单一最优解而是求一组Pareto最优解。一个方案在改进一步改善某个目标的同时不得不至少牺牲另一个目标时这个方案就是Pareto解。比如一个方案运行成本更低但碳排放更高另一个方案碳排放更低但成本更高两者都可能是Pareto解。最终怎么挑那是决策者的事算法只负责把候选方案集完整呈现出来。2.2 MOPSO的核心机制MOPSO的基本骨架是标准粒子群优化。每个粒子代表一组调度方案粒子有位置和速度在搜索空间里不断飞行靠个体历史最优pbest和全局引导者gbest调整飞行方向。跟标准PSO相比MOPSO有几处关键改动维护一个外部档案Archive保存当前已经找到的非劣解。粒子不求单一全局最优而是从档案里选择一个全局引导者。用Pareto支配关系来判断解的好坏档案里的解两两之间不互相支配。档案容量满了以后用拥挤度距离或者网格法裁剪留下分布更均匀的解避免所有粒子都堆在前沿的某个区域。这套机制保证了搜索过程始终沿着非劣方向推进同时保持解的分布性最终能得到一条覆盖较宽的Pareto前沿。2.3 为什么在CCHP场景里适用CCHP运行优化的数学模型一般是非线性、非凸的存在机组启停这类离散变量还有大量不等式约束。传统数学规划方法处理起来要么做线性化近似要么依赖专业求解器建摸门槛不低。粒子群算法不依赖梯度信息目标函数能算出来就行逻辑直白Matlab里把目标函数实现、约束处理一下算法主体几百行代码就能搭起来。这也是很多研究者和工程师把它当作快速验证工具的原因。它虽然不保证全局最优但配合合适的变异策略和档案更新机制在常规算例里往往能给出工程上足够好的调度方案。3. Matlab建模决策变量、目标函数和约束怎么布局做优化代码最忌讳上来就写算法。建模永远是第一位模型写不清楚算法再先进也算不出可信结果。3.1 决策变量的摆放方式我把24个时段的所有连续决策变量按顺序拼成一个长向量这样粒子群算法处理起来最方便。向量结构大概是这样的% 位置向量的切片规划 % 下标1:24 燃气轮机发电功率单位kW % 下标25:48 电制冷机输入功率单位kW % 下标49:72 燃气锅炉供热量单位kW % 下标73:96 蓄电放电功率单位kW % 下标97:120 蓄电充电功率单位kW % 下标121:144 蓄热放热功率单位kW % 下标145:168 蓄热充热功率单位kW为什么要把变量都放在一个向量里因为粒子群算法本身不关心变量结构粒子位置就是一个实时更新的向量。把变量统一拼成长向量以后适应度计算函数的唯一输入就是这个向量代码接口干净调试方便。后面如果想加入启停变量直接扩展向量长度就行。3.2 目标函数怎么写成代码我在项目里采用三个目标运行成本、一次能源消耗量、碳排放量。运行成本包括天然气费、购电费、设备运行维护费单位通常取万元一次能源消耗量用GJ或者吨标准煤表示碳排放量用吨表示。核心计算代码的骨架如下function obj evaluateObjective(x) % x: 决策变量向量 % obj: [cost, primaryEnergy, CO2] % 从x里切片 P_gt x(1:24); P_ec x(25:48); Q_gb x(49:72); P_dis x(73:96); P_ch x(97:120); H_dis x(121:144); H_ch x(145:168); % 燃气轮机发电消耗的天然气量单位m3/h fuel_gt P_gt ./ eta_gt * 3600 / LHV; % 购电功率由电平衡反推 P_buy E_load P_ec P_ch - P_gt - P_dis; % 运行成本 cost_fuel_gas fuel_gt .* price_gas Q_gb ./ eta_gb * 3600 / LHV .* price_gas; cost_buy max(0, P_buy) .* price_e; cost_om P_gt .* c_om_gt P_ec .* c_om_ec Q_gb .* c_om_gb; cost_total sum(cost_fuel_gas cost_buy cost_om) / 10000; % 换成万元 % 一次能源消耗量 primary_energy sum(fuel_gt * 3600 / LHV Q_gb / eta_gb * 3600 / LHV max(0, P_buy) * 0.3); % 碳排放量 carbon_total sum(fuel_gt * coef_carbon_gas max(0, P_buy) * coef_carbon_grid) / 1000; obj [cost_total, primary_energy, carbon_total]; end这段代码是概念示意实际跑的时候最需要注意的是量纲统一。电价单位是元/kWh购电功率单位是kW乘积再乘以时段长度才是费用。一次能源里电网电力的折算系数要看项目所在区域电网的供电煤耗各地差异很大。3.3 约束条件怎么落地我把约束处理分成两类。第一类是等式约束也就是电、热、冷三个能量平衡。这类约束不是靠罚函数硬扛而是利用边界变量反推。比如购电功率不是独立决策变量而是由电平衡公式算出来的P_buy E_load P_ec P_ch - P_gt - P_dis;这样只要把燃气轮机出力、电制冷功率、蓄电充放都作为粒子变量电平衡约束就自然被消解了。但这里有个问题反推出的P_buy可能为负值也就是系统从电网反向卖电。这种情况下要么允许卖电增加一个卖电收益项要么把它截断为零并计入惩罚。第二类是不等式约束主要是设备上下限、爬坡约束、蓄能容量。上下限直接在粒子位置更新后做越界截断蓄能容量约束我通常会通过时序累计来检查逐时刻看蓄电量和蓄热量是否落在可用范围内超出量计入罚函数。罚函数系数怎么选太小约束形同虚设太大目标函数被罚项淹没Pareto前沿会失真。我的经验是先跑一次不带约束的优化看目标函数量级然后让惩罚项量级比目标函数高2到3个数量级。比如成本在万元量级罚函数系数取10^2到10^3比较合适。注意要用相对量别把不同量纲的目标函数直接比较。4. MOPSO算法在Matlab里的完整实现建模完成后进入算法主体。这一部分我拆成四块讲初始化、档案更新、粒子更新、结果输出。4.1 初始化与Pareto排序粒子初始化遵循“随机但尽量可行”的原则。每个变量在上下限之间均匀随机取值。为了加快收敛、减少前期大量越界的情况我会把一部分粒子初始化为参考运行方案比如常见的“以电定热”策略即让燃气轮机发电功率直接跟踪电负荷其余部分完全随机。两者混合起来初期种群质量会好很多。初始化完成后计算全部粒子的三个目标值然后做第一次Pareto排序。Pareto排序我没有依赖外部工具箱手写一个支配函数就能用function flag dominates(obj1, obj2) % obj1支配obj2返回true否则false % 假设所有目标都是越小越好 flag all(obj1 obj2) any(obj1 obj2); end对所有粒子两两判断把不被任何其他粒子支配的个体放入外部档案。这里要注意初期档案往往塞满了粒子而目标值差异巨大所以第一次裁剪会很剧烈不用担心这是正常现象。4.2 外部档案维护与拥挤度距离档案容量一般设为50到200。容量太小容易丢失搜索多样性太大又会拖慢每步挑选领导者的速度。我最常用的值是100。档案满了以后用拥挤度距离排序裁剪。拥挤度距离的含义是某个解在目标空间的前后邻居围成的“矩形”边长总和越大的解说明它周围越空旷应该优先保留这样才能让Pareto前沿铺得更均匀。计算方式如下function cd crowdingDistance(objectives) % objectives: n行m列n是解数量m是目标数量 [n, m] size(objectives); cd zeros(n, 1); for j 1:m [~, idx] sort(objectives(:, j)); fmin objectives(idx(1), j); fmax objectives(idx(end), j); cd(idx(1)) inf; cd(idx(end)) inf; for i 2:n-1 if fmax ~ fmin cd(idx(i)) cd(idx(i)) ... (objectives(idx(i1), j) - objectives(idx(i-1), j)) / (fmax - fmin); end end end end裁剪时把拥挤度最小的粒子逐个淘汰直到档案长度满足上限。边界解拥挤度设为无穷大这个细节非常重要它保证了前沿的两个端点永远不会被剪掉。4.3 领导者选择与速度位置更新标准PSO的速度位置更新公式是v_new w * v_old c1 * r1 * (pbest - x) c2 * r2 * (gbest - x) x_new x_old v_new在MOPSO里pbest的更新不再按单目标比较而是用Pareto支配关系判断如果新粒子支配旧pbest则更新pbest。如果旧pbest支配新粒子则保留旧值。如果两者互不支配以一定概率我取0.5更新pbest给搜索增加多样性。gbest从档案里选但直接随机选一个会导致粒子在早期聚集到档案的某个局部区域。我用两阶段选择第一步从档案里随机抽出10个候选解第二步在这10个候选里挑拥挤度距离最大的那个。这样既保证引导解来自Pareto前沿又避免了总选同一个端点解。4.4 完整主循环框架与参数经验主循环框架如下for iter 1:maxIter w w_max - (w_max - w_min) * iter / maxIter; % 惯性线性递减 for i 1:popSize gbest leaderSelection(archive); v(i,:) w*v(i,:) c1*rand*(pbest(i,:) - x(i,:)) c2*rand*(gbest - x(i,:)); x(i,:) x(i,:) v(i,:); x(i,:) min(max(x(i,:), lb), ub); % 越界截断 obj evaluateObjective(x(i,:)); % 和pbest做Pareto支配比较 if dominates(obj, pbestObj(i,:)) pbest(i,:) x(i,:); pbestObj(i,:) obj; elseif rand 0.5 ~dominates(pbestObj(i,:), obj) pbest(i,:) x(i,:); pbestObj(i,:) obj; end end % 将本轮非劣解加入临时档案 tempArchive [archive; x, objEvaluated]; % Pareto排序剔除被支配解 archive paretoSort(tempArchive); % 超出容量则裁剪 archive trimArchive(archive, archiveSize); end惯性权重w从0.9线性降到0.4学习因子c1、c2我取1.5左右。种群规模常取100到200迭代次数300到500。常规24时段算例Matlab单次运行大概一分钟到几分钟能收敛。这个量级的参数基本够用。注意粒子速度上限对结果影响很大。变量值域在0到200范围内的话速度上限我习惯取变量范围的10%左右比如20。速度上限过大粒子在边界来回振荡过小只能做局部搜索很难跳出初值邻域。5. 算例结果与分析算法跑通之后最关键的步骤是读结果。很多人跑到一张图就收手其实怎么读图、怎么选解才是整条技术链里最有价值的部分。5.1 典型算例怎么设置我给出一个可复现的参考设置。24时段某个典型夏季日电负荷集中在上午和傍晚冷负荷在午后达到峰值热负荷很小主要是生活热水。能源价格使用峰谷电价结构高峰时段电价约1.1元/kWh平段0.7元/kWh低谷0.4元/kWh天然气价格取2.5元/m3。燃气轮机发电效率取32%余热回收效率取50%吸收式制冷COP约1.2电制冷COP约3.5燃气锅炉效率取90%。系统容量按实际项目场景设置就不写死了比如燃气轮机额定发电功率2000kW电制冷额定功率800kW蓄电容量2000kWh蓄热容量5000kWh。这些参数直接决定优化出来的Pareto前沿形状。前端设备参数不同前沿的曲率和拐点位置都会变。5.2 Pareto前沿怎么看运行结束后把档案里的解在目标空间画出来。最直观的是画成本-碳排放二维图第三个目标一次能源消耗可以用颜色映射或者气泡大小表示。典型结果如下Pareto前沿在左下角朝右上角延伸。成本低的方案倾向于在低谷时段多购电把燃气轮机压到低出力这样运行成本低但电网电对应的碳排放系数往往较高所以碳排放偏高碳排放低的方案倾向于让燃气轮机多发、用吸收式制冷承担冷负荷购电少了、天然气耗量上去了成本自然升高。前沿的弯曲程度体现了目标之间的冲突强度。如果前沿接近一条直线说明两个目标近似线性折中如果出现明显的拐点或凹陷段说明设备组合有不同运行模式在切换。曲线凹下去那一段里往往藏着工程上更有性价比的方案值得重点关注。5.3 怎么从一堆解里选运行方案Pareto前沿上的每个解都对应一种24小时设备出力计划但实际运行只能执行一个。我在项目里常用两种选解方式方式一是看曲线的工程可执行性。有些解虽然目标值很优但设备出力曲线剧烈波动蓄电充放频繁这种方案在真实系统里很难落地。优先选出力平滑、运行模式切换次数少的解。方式二是用理想点距离法。先算出三个目标分别在档案里的最小值组成一个“理想点”然后计算每个档案解到理想点的归一化欧氏距离选距离最小的作为折中最优解。再人工检查这个解的出力曲线是否平滑如果平滑就采用如果不平滑就在附近再挑一个。两种方式结合下来结果比较稳定而且给项目汇报时理由也站得住脚。6. 常见问题与排查方法这部分内容来自我实际跑Matlab代码时遇到的问题按出现频率从高到低排。6.1 结果一直不收敛Pareto前沿散乱原因大多是种群规模太小或者档案容量太小导致多样性丢失。把种群从50提到150迭代从200提到400通常就有明显改善。另一个隐蔽因素是领导者选择如果每次都选档案里拥挤度最小的解粒子会快速聚到局部区域前沿反而拉不长。我采用“随机初选拥挤度复选”之后前沿分布明显改善。6.2 约束被违反但解显示为可行这种情况常见于罚函数系数不够。检查方法很笨但很有效跑完以后单独把每个等式约束的最大偏差打印出来。如果电平衡最大偏差超过几十kW说明罚函数强度不足把系数提高一个数量级再看。同时记得提高迭代次数否则前期大量粒子越界搜索效率很低罚函数费用会压制所有粒子的正常飞行。6.3 优化耗时太长三个目标再加上24时段一天下来几万次适应度评估。如果用纯循环嵌套Matlab跑起来确实慢。我的改进方法有三个目标函数里尽量做向量化计算24时段的运算写成数组整体运算不要一个时段一个时段循环。Pareto排序本身是O(n^2)量级空间不大但可以减少每轮临时档案的规模。先对当前粒子全体做快速筛选只把本轮的非劣个体加入临时档案能省不少排序时间。如果模型已经验证正确迭代次数可以适当降到200档案容量保持100对常规结果的影响不大。6.4 优化结果是负值或者出力剧烈振荡主要原因有两个。一是边界处理用截断时速度方向和位置更新不匹配粒子反复撞边界造成相邻时段出力突变。解决方法是在位置更新后增加一个最小变化量检查如果某变量与其上一时刻的值差小于设定值就不做更新。二是变量切片顺序与适应度函数里的切片顺序不一致。这种低级错误最让人抓狂。我的习惯是写一个调试模式随机给一个合理初值逐步打印每个设备出力检查能量平衡是否大体成立。6.5 随机种子导致结果差异大MOPSO是随机算法每次跑出来的Pareto前沿会有细节差异这很正常。为了避免报告里的结果无法复现我在主程序开头固定随机数种子rng(42)相同参数下所有结果可以完全复现。如果要做参数敏感性分析就换不同种子跑若干遍再看统计特征。不要用单次结果下结论至少跑三到五次取分布范围。7. 我踩过的坑和后续可以扩展的方向模型和代码都稳定之后这套东西的扩展空间其实很宽。我把自己试过的和准备尝试的方向列一下。7.1 设备启停变量的混合整数处理连续变量版本有个不足它允许设备在很低的负载率运行工程上可能不合理。比如燃气轮机最小技术出力以下的区间往往被忽略掉或者设备频繁启停。要更贴近实际就得加入启停变量。MOPSO直接处理0/1变量不自然我建议用双层优化内层用粒子群算连续出力外层用另一轮粒子群或者枚举方式算启停方案两层循环用同一个目标函数。代价是计算量成倍增加但得到的运行曲线更符合工程常识。7.2 考虑风光出力和负荷的不确定性实际系统里风光出力、负荷预测、电价都带有不确定性。可以把预测误差建模成多个场景把目标函数改写成期望值或者条件风险价值。MOPSO里的外部档案和拥挤度机制仍然可以用只是评价函数变成对多个场景的平均评估运行时间会进一步上升。这时候给算法加并行计算就很有必要了Matlab的parfor可以直接用在粒子评估环节。7.3 结合机器学习做快速调度最近一个明显的趋势是用强化学习或图神经网络表示CCHP的能量流拓扑训练好的模型能快速输出接近最优的调度策略。进化算法在可解释性和硬约束处理上仍然有优势。两边不是替代关系更像互补。我也在尝试用MOPSO生成的高质量Pareto解作为监督样本训练一个代理模型把非线性的设备效率曲线给拟合出来这样既保留了Pareto语义又能让实时调度的响应速度大幅提升。7.4 工程落地时的二次校核最后提醒一句仿真结果再漂亮放到真实系统里也要经过调度执行层的二次校核。优化算出来的24小时计划本质上还是建议计划必须考虑设备响应时间、最小持续运行时间、检修停机窗口这些细约束。建模仿真阶段如果能把这些因素提前考虑进去后期调试和落地会少走很多弯路。我自己带项目时通常第一步让模型跑通并输出可解释的Pareto前沿图第二步做灵敏度分析第三步才进入细节约束的补充。顺序不能反。先把大局看清再抠细节否则很容易陷入参数调优的泥潭。多目标粒子群在CCHP场景里是很顺手的一套工具但真正让结果有价值的地方永远是模型是否贴合物理规律、结果是否让人看得懂用得上。这个道理放在任何优化问题上都成立。