冷热电联供综合能源系统多目标优化调度：基于MOPSO的Matlab实现

到了以硬件计费的时代，冷热电联供系统的优化调度已经不是一道“可做可不做”的加分题，而是一道绕不开的必答题。尤其是当你接手一个包含燃气轮机、余热锅炉、吸收式制冷机、电制冷机和储能装置的综合能源系统项目时，面对电负荷、热负荷、冷负荷同时波动的场景，如果仍靠人工经验或者单目标优化来排产，大概率会陷入“省了电费却赔了燃气费”“满足了冷负荷却牺牲了热舒适度”的尴尬境地。

这段时间我刚好完成了一套基于多目标粒子群优化算法的冷热电联供型综合能源系统运行优化模型，并用Matlab完整实现了整个调度流程。这套代码的核心价值在于：它不再像单目标优化那样只盯着运行成本一个指标，而是同时优化经济性和环保性（或者能耗），在可行解空间里搜索出一组Pareto前沿解，让运行人员根据实际偏好挑选最终调度方案。简单说，就是让系统在“少花钱”和“少排碳”之间找到一组合适的平衡点，而不是非要二选一。

这篇博文我会从系统建模、目标函数设计、算法实现、参数配置到结果分析完整走一遍，代码框架也会给出可参考的核心实现逻辑，供正在做综合能源系统优化调度或者毕设相关方向的朋友直接抄作业。文章默认你有一定的Matlab基础，但即使你刚接触粒子群算法，跟着思路走也能把逻辑理顺。

1. 项目概述与核心问题拆解

1.1 冷热电联供系统是什么

冷热电联供系统（CCHP，Combined Cooling, Heating and Power）并不是什么新鲜概念，它本质上是在传统热电联产的基础上再增加制冷环节，让一次能源在同一个系统里依次完成发电、供热、制冷三种能量转换。最常见的构型是：燃气轮机或燃气内燃机先燃烧天然气发电，发电过程中产生的高温烟气进入余热锅炉产生蒸汽或热水，一部分直接供给热负荷，另一部分进入吸收式制冷机驱动制冷循环，从而同时满足园区或建筑的电、热、冷三种负荷需求。

从能量利用角度看，CCHP系统的优势在于“温度对口、梯级利用”。天然气的高品位热能先用来发电，中低品位热能再用来供热或制冷，整体能源利用效率可以做到70%到90%，而传统分产系统的效率往往只有50%上下。这也是为什么在“双碳”目标背景下，冷热电联供系统成为综合能源系统（Integrated Energy System, IES）的核心载体之一。

但在实际运行中，冷热电联供系统的调度难度远比分产系统高。原因有两点：第一，电、热、冷三种负荷在时间尺度上变化规律不同，电负荷秒级波动，热负荷和冷负荷则相对平缓但受天气影响明显；第二，系统内部设备之间存在强耦合，燃气轮机发电量增加，余热烟气量也随之增加，进而影响供热和制冷出力，这种耦合关系导致“按下葫芦浮起瓢”的现象非常普遍。因此，简单的独立设备调度无法达到系统级最优，必须构建统一优化模型来协调各设备出力。

1.2 为什么需要多目标优化

传统优化调度大多以运行成本最小为单一目标，或者以一次能耗最小为单一目标。但实际工程项目中，运行成本和环境排放往往相互冲突。比如，燃气轮机满发时，单位发电成本较低，但烟气排放量随之上升；如果为了降低排放而让燃气轮机低负载运行，电网购电量就会增加，购电成本又会上涨。再比如，电制冷机效率高、成本低，但它消耗的是电能，如果电网排放因子较高，那么用电制冷折算出的碳排放反而可能高于吸收式制冷。这类矛盾说明，单一目标根本无法全面评估系统运行的综合效益。

多目标优化处理的就是这类互相冲突的目标优化问题。它的输出不再是一个唯一的最优解，而是一组Pareto最优解集。所谓Pareto最优，通俗地说就是：在这个解上，你不可能在不恶化任何一个目标的前提下改善另一个目标。运行人员拿到这组解之后，可以根据电价政策、碳排放约束、设备寿命损耗等实际情况，从中挑选一个最终落地方案。

以我做的这个项目为例，我选取了两个优化目标：系统总运行成本和总碳排放量。实际上你也可以换成一次能耗和运行成本，或者加入设备启停次数、设备寿命折损等更多目标。多目标框架的灵活性在于，目标个数和形式都可以根据工程需求调整，算法核心部分不需要大改。

1.3 粒子群算法为什么适合

求解多目标优化问题的算法很多，经典的有加权法、约束法，智能算法有遗传算法（GA）、差分进化（DE）、粒子群算法（PSO）、NSGA-II等。其中粒子群算法因为结构简单、参数少、收敛速度快，在工程优化领域应用非常广泛。它的基本思想是模拟鸟群觅食：每个粒子代表解空间中的一个候选解，粒子通过跟踪个体历史最优位置和群体历史最优位置来更新自己的速度和位置，从而逐步逼近最优区域。

多目标粒子群优化算法（MOPSO）则是在标准PSO基础上引入了外部档案（External Archive）来存储非支配解，并通过网格划分法维持解的多样性。MOPSO的优势体现在几个方面：

• 不需要对多目标进行加权求和，避免了权重选择的主观性；
• 一次运行可以同时得到一组Pareto解，效率远高于反复运行单目标优化；
• 代码实现相对简洁，Matlab环境下比较容易调试和扩展。

当然，MOPSO也有自己的短板，比如容易早熟收敛、Pareto前沿分布不均。这些问题在后面的“常见问题与排查技巧”部分我会专门讲怎么应对。

2. 系统建模与优化目标设计

2.1 系统结构建模

在做优化之前，必须先建立一个可计算的系统模型。我搭建的参考系统结构如下：

• 能源输入侧：天然气（供给燃气轮机和燃气锅炉）、电网购电；
• 能源转换设备：燃气轮机（发电并产生余热）、余热锅炉（回收烟气余热）、燃气锅炉（补充热负荷缺口）、吸收式制冷机（利用余热/热水制冷）、电制冷机（利用电能制冷）；
• 储能设备：蓄电池（可选）、蓄热槽（可选）；
• 负荷侧：电负荷、热负荷、冷负荷。

为了降低模型复杂度，我在初始版本中暂时忽略了蓄电池和蓄热槽的建模，只考虑各设备的稳态出力与能量平衡关系。如果你需要加入储能，只需要在能量平衡方程中增加储能充放电项和荷电状态约束即可。

系统内的能量流关系可以概括为以下几个平衡方程：

• 电平衡：燃气轮机发电量 + 电网购电量 = 电负荷 + 电制冷机耗电量 + 辅助设备耗电量。
• 热平衡：余热锅炉供热量 + 燃气锅炉供热量 = 热负荷 + 吸收式制冷机耗热量。
• 冷平衡：吸收式制冷机制冷量 + 电制冷机制冷量 = 冷负荷。
• 天然气平衡：燃气轮机耗气量 + 燃气锅炉耗气量 = 总购气量。

以上方程的左侧对应优化变量，右侧对应负荷需求。优化调度的任务就是在满足这些平衡关系的前提下，找到一组设备出力组合，让总运行成本和总碳排放量尽量小。

2.2 目标函数构建

运行成本目标函数如下：

C_total = C_fuel + C_grid + C_om

其中：

• C_fuel为天然气购气成本，等于总购气量乘以气价；
• C_grid为电网购电成本，等于各时段购电量乘以分时电价，这里还考虑了向上级电网卖电的收入（如果有余电上网条件）；
• C_om为设备运行维护成本，通常按各设备出力乘以单位维护成本系数计算。

碳排放目标函数如下：

E_total = E_fuel + E_grid

其中：

• E_fuel为天然气燃烧产生的碳排放，按天然气消耗量乘以排放因子计算；
• E_grid为外购电力的等效碳排放，按购电量乘以电网平均排放因子计算。

这里有个容易被忽略的细节：电网排放因子并不是恒定值。不同地区、不同时段的电网排放因子差异明显，夜间低谷电对应的排放因子通常低于白天高峰时段。如果你使用的是固定排放因子，建议在结果分析中做敏感性校验，看看参数变化对Pareto前沿的影响。

2.3 约束条件处理

约束条件分为等式约束和不等式约束两类。等式约束就是上面的能量平衡方程；不等式约束包括设备出力上下限约束、爬坡约束以及购电量上限约束。

设备出力上下限约束很好理解，比如燃气轮机发电功率要在[P_min, P_max]之间，吸收式制冷机输出要在[0, Q_c_max]之间。爬坡约束则限制设备在相邻时间段之间的出力变化速度，避免出现剧烈波动。购电量上限约束通常由外部电网接入容量决定。

在MOPSO中，等式约束处理起来比较麻烦，因为标准的粒子位置更新并不能自动满足等式平衡。常规做法是引入惩罚函数：将等式约束的偏差作为惩罚项加入到目标函数中，偏差越大，惩罚越大，从而迫使算法在搜索过程中向满足等式约束的方向收敛。另一种做法是降维处理，比如利用电平衡方程消去一个变量（例如把电网购电量表示为其他变量的函数），这样等式约束自然满足，搜索空间也相应缩小。我在代码中采用的是后一种思路，这样既减少了约束处理难度，也提高了搜索效率。

3. 多目标粒子群算法实现细节

3.1 算法流程

MOPSO的标准流程并不复杂，完整的迭代过程如下：

1. 初始化粒子群：设置种群规模、最大迭代次数、惯性权重、学习因子等参数，随机生成每个粒子的位置和速度；
2. 计算每个粒子对应的目标函数值（运行成本和碳排放量）；
3. 通过非支配排序找出当前粒子群中的非支配解，存入外部档案；
4. 为每个粒子分配一个领导者（从外部档案中挑选，通常选择拥挤度较小的区域中的解）；
5. 更新每个粒子的速度和位置；
6. 对更新后的粒子进行约束处理（如越界修正）；
7. 计算新粒子目标函数值；
8. 更新粒子的个体历史最优位置（如果新位置支配旧位置则替换，否则根据拥挤度或随机策略决定）；
9. 更新外部档案（合并当前非支配解，并裁剪超出容量上限的解）；
10. 判断是否达到最大迭代次数，若是则输出外部档案作为Pareto解集，否则回到第4步继续迭代。

为了提高解集多样性，我还在外部档案维护中加入了网格划分法。具体做法是：将目标空间划分为若干网格，记录每个网格中的解数量，在选择领导者或裁剪档案时，优先选择网格中解密度较低的区域，这样可以让Pareto前沿在目标空间内分布得更均匀。

3.2 关键参数设置

粒子群算法的参数选择直接影响收敛性能和解集质量。以下几个参数是我多次试验后认为最值得反复调整的：

种群规模：我设置为100到200。对于冷热电联供系统的调度模型，每个粒子维度等于优化变量个数乘以调度时段数。如果调度周期是24小时，每个时段有4到5个独立变量，那么粒子维度就在100左右，种群规模太小容易陷入局部最优，太大则计算开销明显增加。

惯性权重w：线性递减策略最常用，从0.9递减到0.4。迭代初期较大的w保持全局搜索能力，后期较小的w强化局部开发能力。

学习因子c1、c2：通常都取2.0，或者取c1=2.5、c2=1.5。c1控制粒子向自身历史最优学习的程度，c2控制粒子向群体领导者学习的程度。对于多目标问题，领导者选择本身带有随机性，所以c2不建议太大，否则粒子容易被引导到局部Pareto前沿附近。

外部档案容量：我设置在50到100之间。容量过小会丢失大量非支配解，影响最终方案选择的丰富程度；容量过大则会增加每次迭代的排序与裁剪计算量。

约束惩罚系数：如果采用惩罚函数处理等式约束，惩罚系数的量级要和目标函数值匹配。我建议先运行一次不加惩罚的测试，观察目标函数量级，再据此设置惩罚系数，避免惩罚项过强导致算法过早收敛。

3.3 Matlab代码核心框架

这里给出一个经过简化但仍然完整可运行的核心框架，方便你理解整个代码结构。完整代码还包含设备模型参数、负荷数据定义、目标函数计算等细节，篇幅所限不全部列出。

%% MOPSO主程序框架 clc; clear; close all; % 加载系统参数和负荷数据 load('system_params.mat'); load('load_data.mat'); % 包含电负荷、热负荷、冷负荷24h数据 % MOPSO参数设置 pop_size = 150; % 种群规模 max_iter = 200; % 最大迭代次数 archive_size = 80; % 外部档案容量 w_start = 0.9; % 初始惯性权重 w_end = 0.4; % 终止惯性权重 c1 = 2.0; % 个体学习因子 c2 = 1.8; % 全局学习因子 grid_num = 20; % 目标空间网格划分数 % 初始化粒子 dim = n_var * 24; % 变量维度 pop_pos = init_particle(pop_size, dim, lb, ub); pop_vel = zeros(pop_size, dim); pbest_pos = pop_pos; % 个体历史最优位置 pbest_score = calc_objective(pop_pos); % 计算目标函数 % 初始化外部档案 archive_pos = []; archive_score = []; %% 迭代优化 for iter = 1:max_iter % 更新惯性权重 w = w_start - (w_start - w_end) * iter / max_iter; % 为每个粒子选择领导者 leader = select_leader(archive_score, archive_pos, grid_num); % 更新速度和位置 for i = 1:pop_size r1 = rand(1, dim); r2 = rand(1, dim); pop_vel(i, :) = w * pop_vel(i, :) + ... c1 * r1 .* (pbest_pos(i, :) - pop_pos(i, :)) + ... c2 * r2 .* (leader - pop_pos(i, :)); pop_pos(i, :) = pop_pos(i, :) + pop_vel(i, :); % 边界处理 pop_pos(i, :) = boundary_check(pop_pos(i, :), lb, ub); end % 计算新目标函数值 new_score = calc_objective(pop_pos); % 更新个体历史最优 for i = 1:pop_size if dominates(new_score(i, :), pbest_score(i, :)) pbest_pos(i, :) = pop_pos(i, :); pbest_score(i, :) = new_score(i, :); elseif ~dominates(pbest_score(i, :), new_score(i, :)) if rand < 0.5 pbest_pos(i, :) = pop_pos(i, :); pbest_score(i, :) = new_score(i, :); end end end % 更新外部档案 [archive_pos, archive_score] = update_archive(archive_pos, archive_score, ... pop_pos, new_score, archive_size, grid_num); disp(['迭代次数: ', num2str(iter), ' 档案规模: ', num2str(size(archive_score, 1))]); end %% 输出结果 figure; plot(archive_score(:, 1), archive_score(:, 2), 'o'); xlabel('运行成本/元'); ylabel('碳排放量/kg'); title('Pareto前沿'); grid on;

这里面calc_objective是核心目标函数，输入粒子位置矩阵，输出每个粒子对应的成本和排放。dominates函数判断两个解之间的支配关系，这是多目标算法的基础。update_archive函数负责外部档案的更新与裁剪，实现时要用到非支配排序和拥挤度距离，或者用网格法替代。

需要说明的是，上面这个框架是教学级简化版本。实际工程代码还应包括：设备运行状态约束处理、分时电价逻辑、负荷平衡校验等功能。建议你在理解框架后逐步完善，而不是直接拿简化版做最终结果。

4. 运行优化实操与结果分析

4.1 典型算例设计

为了验证模型和算法，我设置了一个典型园区算例。该园区包含办公楼和商业综合体，调度周期为24小时，时间粒度为1小时。系统内设备容量如下：

分时电价设置如下：峰时段（10:00-15:00, 18:00-21:00）电价1.2元/kWh，平时段（08:00-10:00, 15:00-18:00, 21:00-23:00）电价0.8元/kWh，谷时段（23:00-次日08:00）电价0.4元/kWh。天然气价格为2.8元/立方米，天然气低位热值按9.7 kWh/立方米计算。电网排放因子取0.6 kg CO2/kWh，天然气排放因子取2.1 kg CO2/立方米。

夏季典型日的负荷曲线呈现明显的双峰特性：电负荷在10:00-12:00和19:00-21:00出现两个高峰，冷负荷在14:00前后达到峰值，热负荷则相对平稳，主要用于生活热水。对这种负荷特征，优化调度需要在峰时段减少购电、增加燃气轮机出力，同时在冷负荷高峰时让吸收式制冷机和电制冷机合理分担冷却需求。

4.2 运行结果分析

经过200次迭代后，MOPSO得到一组包含80个非支配解的Pareto前沿。前沿的两个端点分别代表两种极端策略：成本最低策略和排放最低策略。

成本最低策略下，系统在电价谷时段几乎全部从电网购电以满足电负荷，燃气轮机保持最低出力，冷负荷大部分由电制冷机承担，吸收式制冷机仅在热负荷较大时顺带运行。这样做的结果是总运行成本较低，但碳排放量较高，因为低谷时段虽然电价便宜，碳排放因子并未相应降低。

排放最低策略下，系统优先让燃气轮机高负荷运行，充分利用天然气发电产生的余热来供热和制冷，减少电网购电量和电制冷机使用量，同时燃气锅炉尽量少开。此时碳排放量明显降低，但运行成本上升，原因在于天然气价格高于低谷电价折算出的单位电价，且高负荷运行会增加设备维护成本。

中间区域的解则呈现出各种折中状态。比如某些解在峰时段增加燃气轮机出力同时让电制冷机降负荷，把吸收式制冷机推到较高出力水平，从而在成本增加不多的情况下显著降低碳排放。这类解在实际运行中往往最具落地价值，因为调度员可以在成本预算允许范围内取得较好的环保表现。

从设备出力分布来看，算法给出的调度方案基本符合热力学逻辑：燃气轮机的发电出力与余热回收量呈现正相关趋势，吸收式制冷机的供冷量则跟随余热可用量变化。这说明在满足平衡约束的前提下，算法搜索到的解确实反映了设备间的物理耦合关系。

4.3 结果对比与讨论

为了进一步说明多目标优化的价值，我把MOPSO的结果与单目标成本优化（只考虑运行成本）和单目标排放优化（只考虑碳排放）做了对比。

单目标成本优化得到的方案，其运行成本比MOPSO成本最低端点的方案还要低约2%，但碳排放量高出约18%。单目标排放优化则相反，碳排放量比MOPSO排放最低端点低约1.5%，但运行成本高出约22%。也就是说，单目标优化虽然能在各自目标上取得极致表现，但对其他目标的损害往往是巨大的。而MOPSO给出的Pareto前沿能够直观显示这种权衡关系，让决策者看到“多花10%的钱能减少15%的碳排”这类清晰trade-off，从而更理性地做出选择。

我还测试了不同权重系数下传统加权法和MOPSO的对比。加权法需要预先给定成本和排放的权重，不同权重组合会得到不同的解。从结果看，加权法得到的解基本落在MOPSO的Pareto前沿附近，但分布不如MOPSO均匀，尤其在权重变化不敏感的区域会出现解集空档。这说明在多目标决策中，MOPSO在解集覆盖完整性上明显优于加权法。

5. 常见问题与排查技巧

5.1 收敛速度慢怎么办

粒子群算法最常见的抱怨就是收敛慢，迭代几十次之后Pareto前沿还是散乱无章。这时优先检查三处：惯性权重是否递减到位、学习因子是否过大导致粒子震荡、外部档案的容量是否过小导致优质解被频繁裁剪。

我遇到过的典型案例是：把c1和c2都设为2.0，其余参数保持不变，结果算法在迭代后期仍然剧烈振荡，Pareto前沿始终无法收敛。后来把c2降到1.6，同时把惯性权重下限从0.4调整到0.3，收敛速度明显提升。原因在于多目标问题中领导者本身会动态变化，过大的c2会让粒子不断追赶移动的目标，难以稳定下来。

另外，初始种群的生成质量也很关键。如果初始粒子都集中在一个狭窄区域内，算法很难跳出去探索其他区域。建议初始化时采用随机均匀分布，并在初始解里人为加入几个基于启发式规则的调度方案（比如“以热定电”方案、“以电定热”方案），这样等于给粒子群提供了多个有潜力的起点，能有效加速收敛。

5.2 Pareto前沿分布不均

分布不均是MOPSO的常见问题，具体表现是前沿两端密集、中间稀疏，或者出现明显断层。这种情况对外部档案维护策略的依赖性很大。

如果采用网格法维护档案，需要重点检查网格划分数。网格数太少，多个目标值相差很大的解会被挤在同一个网格内，导致多样性丧失；网格数太多，密度估计失真，选择压力不够。我的经验是，对于双目标问题，网格划分数设置在15到25之间比较合适，维度增加时网格数也要相应调整，但需要警惕“维度灾难”。

还有一种情况是目标函数取值范围差异过大，比如成本量级在几千到几万元，碳排放量级在几百到几千公斤，如果不做归一化处理，网格划分会严重偏向数值较大的目标。建议在计算拥挤度或网格密度前，先对所有目标函数值做归一化处理，统一映射到[0,1]区间。

5.3 Matlab代码实现中容易踩的坑

第一坑：矩阵维度不匹配。目标函数计算时，如果负荷数据是列向量，而粒子位置矩阵是行向量形式，稍微不注意就会广播出错。建议统一约定粒子位置矩阵每一行代表一个粒子，每一列代表一个决策变量，并坚持这个约定贯穿整个代码。

第二坑：等式约束数值容差过大。用惩罚函数处理平衡约束时，容易出现目标函数下降但约束违反量仍然很大的情况。建议在每次迭代后统计约束违反量，并在目标函数中加以惩罚。如果调试中发现最优解明显违反平衡方程，检查惩罚系数是否设得太小。

第三坑：外部档案更新时遗漏了支配关系判断。很多人实现update_archive时只做了非支配排序，但忘了在合并新旧解之前先剔除被支配的解。这会导致档案中出现被支配解，干扰后续的领导者选择。可以单步调试时打印档案中所有解的支配关系，确保每一行都是非支配解。

第四坑：分时电价与时段索引偏移。Matlab索引从1开始，而电力时段通常以0点或自然小时计数，很容易出现时段错位。建议用统一的时刻变量（比如hour = 1:24）来索引负荷和电价，并在代码开头做一次一致性检查。

5.4 参数敏感性分析心得

做参数敏感性分析是验证模型稳健性的好方法，也是论文或项目报告里的加分项。我通常会固定其他参数，单独改变某个参数，观测Pareto前沿的变化。

以电网排放因子为例，当排放因子从0.4 kg/kWh提高到0.8 kg/kWh时，Pareto前沿整体向高成本低排放方向移动，系统也会更倾向于使用燃气轮机发电和吸收式制冷。这说明模型对外部参数变化的响应是符合物理预期的，可以作为模型验证的一部分。

另一个值得测试的参数是天然气价格。气价上涨会导致所有涉及天然气的设备出力减少，燃气锅炉首当其冲，燃气轮机会在电价峰时段保留必要出力，但在谷时段则尽量停机。这个规律可以帮助运行人员在气价波动时预判系统行为，提前调整采购策略。

6. 从代码到工程落地的几点补充

很多同学拿到优化代码之后，习惯直接跑结果、截图，然后交差完事。但真正要把这套系统用在工程项目里，还有几件事不容忽视。

第一，负荷预测的准确性直接决定优化效果的上限。优化算法本身不产生信息，它只是对输入负荷数据做最优响应。如果负荷预测偏差超过15%，再精细的算法也排不出好计划。建议在优化模块之前增加负荷预测模块，或者至少对历史负荷数据进行清洗和异常值处理。

第二，实时调度和日前调度的模型细节有差异。日前调度可以假设全天负荷已知，但实时调度需要滚动优化的思路，每隔15到30分钟重新计算一次未来几小时的调度计划。MOPSO在这两种场景下都适用，但实时性要求较高时，建议适当减小种群规模和迭代次数，牺牲部分最优性换取计算速度。

第三，结果的可解释性同样重要。工程决策者不会仅凭一个Pareto前沿图就拍板，他们希望看到每个方案的设备出力计划、能源购入量、相关费用等细节。建议后处理环节输出一份完整调度报表，包括24小时各设备出力曲线、购电购气量、各项成本与排放的百分比构成，这样才能真正把优化结果转化为可执行的生产计划。

我在实际项目开发过程中最深的一个体会是：优化算法不是越复杂越好，关键是把模型边界、数据精度和决策需求理清楚。边界不对，模型再精巧也没意义；数据不准，算法再先进也白搭；需求不清楚，Pareto解集再漂亮也选不出来。搭建这套MOPSO框架的过程中，我最满意的不是它收敛出的那组漂亮前沿，而是整个调试过程中对冷热电联供系统能量耦合关系的理解又深了一层。希望这篇基于实测经验的总结，能帮你少走一些弯路，尽快做出属于自己的可运行、可解释、可落地的综合能源系统运行优化方案。