做综合能源系统运行优化这项工作，我接触最多的一条技术线就是冷热电联供（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(i+1), 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场景里是很顺手的一套工具，但真正让结果有价值的地方，永远是模型是否贴合物理规律、结果是否让人看得懂用得上。这个道理放在任何优化问题上都成立。
