基于POA的冷热电联供系统日前优化调度与MATLAB实现

冷热电联供（CCHP）系统的日前调度，本质上是一个多变量、多约束、强耦合的非线性优化问题。既要盯紧电、热、冷三条能量流的实时平衡，又要跟随分时电价灵活调整购电策略，再加上蓄能装置的状态约束，问题维度随便就能破百。这两年我一直在用元启发式算法折腾这类问题，实测了几种主流算法之后，POA（鹈鹕优化算法）是其中收敛表现最稳的一个。这篇文章就围绕一套我实际在运行的MATLAB代码展开，从算法原理、系统建模、代码架构到调参排坑，把整个流程掰开揉碎讲清楚，希望对正在做综合能源优化或者准备入坑智能算法调度的朋友有点参考价值。

1. 先把POA算法讲透：鹈鹕捕鱼和寻优有什么关系

1.1 从生物学行为到数学公式

POA的全称是Pelican Optimization Algorithm，鹈鹕优化算法，2022年提出。核心思想是模拟鹈鹕在觅食过程中的两块招牌动作：高空俯冲逼近猎物，以及贴着水面滑翔并用翅膀拍水驱赶鱼群。这两套动作在算法里被拆成两个阶段，分别对应全局搜索和局部开发。

第一阶段叫探索阶段，位置更新公式长这样：

X_new = X + rand × (P - I × X)

其中P是当前迭代中随机生成的猎物位置，I是一个随机整数，取1或2。可以这样理解：鹈鹕朝着自己锁定的猎物方向飞行，如果I=1，更新步长适中，偏向局部修正；如果I=2，步长直接翻倍，个体更容易跳到搜索空间较远的区域，保证种群不会过早聚成一团。这个阶段全局搜索能力强，尤其适合冷热电联供这种可行域分散的问题。

第二阶段叫开发阶段，公式是：

X_new = X + 0.2 × (1 - t / T_max) × (2 × rand - 1) × X

这里t是当前迭代次数，T_max是最大迭代次数。系数0.2×(1-t/T_max)会随迭代进行从0.2附近逐渐衰减到0。这意味着迭代前期扰动幅度大，算法还在大范围搜索的尾巴上；迭代后期扰动幅度越来越小，相当于鹈鹕在水面小范围精确修正位置，搜索精度逐步提升。这个自适应衰减的思路和模拟退火的降温过程异曲同工，但实现起来要简单得多。

还有一点容易被忽略：POA每次位置更新后不是无条件采纳的。更新完要重新计算适应度，如果新位置的适应度更差，就留在原地。这种“贪心接受”策略保证了种群整体质量不会退化，也是POA在许多标准测试函数上表现稳定的重要原因。

1.2 跟粒子群、遗传算法放一起比，优势在哪

我拿同一个冷热电联供的24小时调度算例做过横向对比，统一使用相同的种群规模和迭代次数，最终结果如下（相对值以粒子群结果为基准）：

这只是一次测试的相对结果，不代表绝对优劣，但趋势已经很明显：POA在冷热电联供这类“维度不算高、约束一堆”的问题上，收敛速度确实比粒子群快一大截，最终解的质量也更好。它真正的优势可以总结成三条。

第一，参数少到几乎没有。粒子群要调惯性权重、两个学习因子，遗传算法要调交叉概率、变异概率、选择策略，POA只需要设置种群规模和迭代次数，调参成本极低。对工程项目来说，少一个参数就少一个坑。

第二，贪心更新策略天然适合带约束的问题。冷热电联供的可行域往往是破碎的，很多元启发式算法更新完很容易跑飞，POA因为每次更新都做“更好才替换”，配合罚函数能比较稳地往可行域方向收敛。

第三，代码实现非常简洁。整个算法核心不到60行MATLAB，不依赖任何额外工具箱，这对实际落地来说太重要了。后文我会把代码完整贴出来。

2. 冷热电联供系统建模：这步做不好，后面全是白费

2.1 系统结构和能量流

冷热电联供系统的典型结构是这样的：天然气进入燃气轮机（或者内燃机）发电，发电后的高温烟气进入余热回收装置，产生蒸汽或热水，一部分直接对外供热，一部分驱动吸收式制冷机制冷。余热不够的时候，燃气锅炉补燃；冷量不够的时候，电制冷机补充。整套系统和电网并网运行，电不够就买电，电多了可以选择卖给电网（取决于当地政策）。

从能量流角度看，系统里有电、热、冷三种能量形式，设备之间存在强耦合。这就决定了它的优化模型不会像单设备调度那样简单——某个时段燃气轮机多发电，烟气余热也会同步增多，制冷机的输入热量和供热侧的可用热量都会跟着变化，决策变量牵一发而动全身。

我常用的设备模型包括以下几种：

燃气轮机的发电功率与燃料消耗关系，用近似线性模型表示：

F_MT(t) = a × P_MT(t) + b

余热回收量跟随发电功率变化：

Q_HR(t) = η_recovery × Q_exhaust(t)

吸收式制冷机的制冷量和输入热量挂钩：

Q_AC(t) = COP_AC × Q_hot_in(t)

电制冷机就是用电换冷：

Q_EC(t) = COP_EC × P_EC(t)

燃气锅炉相对简单：

Q_GB(t) = η_GB × F_GB(t)

储能设备（蓄热罐、蓄电池）可以统一用一个SOC递推模型：

SOC(t+1) = SOC(t) + η_ch × P_ch(t) - P_dis(t) / η_dis

当然，具体项目里这些模型要根据实际设备参数标定。工程上常见做法是查设备铭牌数据，拿额定工况点做拟合；如果厂家给了效率曲线，就直接按曲线取值。千万不要拍脑袋填系数，一套错误参数会让优化结果完全失去物理意义。

2.2 能量平衡约束和运行约束

冷热电联供的约束条件分三类，缺一不可。

第一类是能量平衡约束，这是硬约束中的硬约束：

电平衡方程：

P_MT(t) + P_grid_buy(t) = P_load(t) + P_EC(t) + P_ch(t) - P_dis(t)

热平衡方程：

Q_HR(t) + Q_GB(t) = Q_heat_demand(t) + Q_AC_in(t)

冷平衡方程：

Q_AC(t) + Q_EC(t) = Q_cooling_demand(t)

需要注意，热平衡里不只包含热负荷，还要算上吸收式制冷机消耗的热量。很多新手第一次建模，只写了热负荷，漏了制冷机热输入，结果能量怎么都对不上，排查半天才发现是这里的问题。

第二类是设备出力上下限约束。每台设备都有最小出力和最大出力限制，燃气轮机还可能有爬坡率约束。储能设备要限制SOC范围，防止过充过放，一般把SOC限制在0.1到0.9之间。

第三类是运行状态约束。如果设备带启停状态变量（0-1变量），问题就变成了混合整数优化，纯连续优化算法处理起来要想办法。我的做法是把启停变量单独拿出来，用二进制编码去处理，或者更省事的做法是：在目标函数里加入启停成本惩罚项，让算法自己权衡“多运行一小时”和“多一次启停”之间的经济性，从而隐式决定设备的运行状态。

2.3 目标函数怎么定

冷热电联供的调度目标，工程上最常用的是运行成本最小化。运行成本包括三块：

C_total = C_fuel + C_grid + C_OM

燃料成本是天然气消耗量乘气价，购电成本是各时段购电量乘分时电价，运行维护成本通常简化成各设备出力乘以维护系数。

如果项目有碳排放考核，还要加上碳税或者碳排放成本。我做过一个多目标版本，把运行成本和碳排放量同时放进目标函数，用线性加权法求解。这里有个关键点要提一下：权重系数不能想当然地随便定。建议先分别跑一次单目标优化，摸清两个目标函数的大致量级范围，再根据量级去定权重，否则很可能出现某个目标在整个优化过程中几乎不起作用的情况。

还有一套思路是加碳排放约束，而不是把碳排放放进目标函数。比如规定全天碳排放总量不能超过某个限额，这样优化目标仍是经济性最优，碳排放作为约束条件存在。这种方式更适合那些已经有明确碳排放指标的项目。

3. MATLAB代码架构：从数据准备到算法主循环

3.1 文件结构和模块划分

整套代码我习惯分成五个部分：主程序、数据初始化、目标函数、约束检查、POA算法核心。文件过多不好维护，五个部分刚好够用。

主程序负责加载负荷数据、设定算法参数（种群规模、迭代次数）、初始化种群、调用POA主循环、输出结果。数据初始化模块把24小时的负荷曲线、分时电价、设备参数整理成结构体，方便后续所有函数调用。

我一般用一个结构体把数据包起来，比如param.P_MT_max、param.eta_recovery这样，可读性好，也不容易记混。所有参数放在代码开头，集中修改，不用跑到函数内部去翻。

具体来说，数据初始化里要准备的基础数据包括：

• 24小时电负荷曲线、热负荷曲线、冷负荷曲线（单位：kW）
• 分时购电价（单位：元/kWh）
• 天然气价格（单位：元/m³）
• 燃气轮机额定功率、发电效率、余热回收效率
• 吸收式制冷机COP、电制冷机COP
• 燃气锅炉效率
• 储能容量、充放效率、SOC上下限

这些数据直接决定了优化结果是否贴近实际。我曾经遇到过负荷曲线单位不一致的问题——设计院给的数据是MW，我代码里默认是kW，算出来的方案偏差了三个数量级，排查了很久才找到原因。建议在数据初始化模块里做好单位统一，并且加一行注释标明所有数据的单位。

3.2 决策变量编码：维度怎么设计

冷热电联供的决策变量，以24小时调度为例，每个时段需要决定的通常是：燃气轮机出力、燃气锅炉出力、吸收式制冷机热输入量、电制冷机出力、蓄能装置充放功率，可能还有电网购电功率。

把所有设备24个时段的出力全部串成一个行向量，维度就是：

dim = 设备类型数 × 24

我举个例子。假设系统有5类可控设备，则决策变量维度就是5×24=120。每个变量归一化到[0,1]区间，在目标函数里再映射回实际物理范围。这么做的好处是POA的初始化和位置更新都不用处理复杂的边界问题，只要在映射时处理好约束就行。

x_actual = lb + x_norm × (ub - lb)

在实际处理中，我还会把每类设备的变量排成连续的一段，而不是交叉排列。比如第1到24维是燃气轮机出力，第25到48维是燃气锅炉出力，这样在解析变量时逻辑清晰，不容易出错。代码里我会在目标函数开头写好注释块，标明变量段落的划分方式。

3.3 约束处理：罚函数还是别的办法

冷热电联供的约束处理，我实测过三种方案，各有适用场景。

第一种是罚函数法，也是我用得最多的。把违反约束的程度作为惩罚项加到适应度函数里。优点是实现简单，代码改动最小；缺点是罚系数难调。罚系数太小，约束形同虚设，算法会大量产出一看就不合理的方案；罚系数太大，可行解和不可行解的适应度差异被拉得过大，算法会绕开那些“差一点就可行”的优质区域，错过潜在的最优解。

第二种是直接修缮法。位置更新完之后，逐个检查约束，如果能量不平衡，按优先级调整设备出力。这个方法在单设备问题上很好用，但在冷热电联供这种多能量强耦合的系统里，修缮逻辑写起来非常复杂，容易引入新的错误。我只在早期做简单案例时用过，后面基本放弃了。

第三种是可行性法则，又被称为Deb法则。在比较两个解时遵循三条规则：可行解永远优于不可行解；两个不可行解中约束违反总量小的更优；两个可行解中目标函数值小的更优。这个方案不需要罚系数，理论上更干净，但需要修改POA算法内部的比较逻辑，实现成本比罚函数法高一些。

工程上我推荐先用罚函数法快速验证算法逻辑，等模型稳定了再切换到可行性法则做精细调优。大部分项目跑下来，罚函数调好了精度也完全够用。

4. 核心代码逐段拆解：目标函数、约束处理和POA实现

4.1 POA主循环代码

上面提到的POA主循环，在MATLAB里的实现大概是这样的：

function [best_pos, best_score, history] = POA_main(obj_func, dim, lb, ub, N, MaxIter) % 初始化种群 X = rand(N, dim) .* (ub - lb) + lb; fitness = zeros(N, 1); for i = 1:N fitness(i) = obj_func(X(i, :)); end [best_score, idx] = min(fitness); best_pos = X(idx, :); history = zeros(MaxIter, 1); for t = 1:MaxIter for i = 1:N % 第一阶段：逼近猎物（探索） I = randi([1, 2]); P = rand(1, dim) .* (ub - lb) + lb; X_new = X(i, :) + rand(1, dim) .* (P - I * X(i, :)); X_new = boundary_check(X_new, lb, ub); if obj_func(X_new) < fitness(i) X(i, :) = X_new; fitness(i) = obj_func(X_new); end % 第二阶段：水面滑翔（开发） r = 0.2 * (1 - t / MaxIter); X_new = X(i, :) + r * (2 * rand(1, dim) - 1) .* X(i, :); X_new = boundary_check(X_new, lb, ub); if obj_func(X_new) < fitness(i) X(i, :) = X_new; fitness(i) = obj_func(X_new); end end % 更新全局最优 for i = 1:N if fitness(i) < best_score best_score = fitness(i); best_pos = X(i, :); end end history(t) = best_score; end end

写代码时有几个细节值得注意。边界检查函数我建议单独写一个，用min(max())把越界变量拉回边界，或者用随机重置到边界内。如果直接用固定边界拉回，当边界比较窄时，大量个体会堆在边界上，种群多样性会快速下降。

第一阶段里的P是随机生成的猎物位置，每只鹈鹕“看到”的猎物都不一样，因此种群探索方向足够分散。这一点跟粒子群所有粒子都朝全局最优飞有很大区别，也是POA不容易早熟的关键原因。

4.2 目标函数和约束检查代码

目标函数的框架大概是这样的：

function cost = obj_func(x) global param load_data price_data % 解析决策变量：1-24 燃气轮机出力，25-48 燃气锅炉出力， % 49-72 吸收式制冷机热输入，73-96 电制冷机出力 P_MT = x(1:24); Q_GB = x(25:48); Q_AC_in = x(49:72); P_EC = x(73:96); % 电平衡：计算电网购电 P_load_total = load_data.elec + P_EC; P_grid = P_load_total - P_MT; % 电网购电不能为负，负的表示售电，按售电竞价结算 P_buy = max(P_grid, 0); P_sell = max(-P_grid, 0); % 热平衡：余热 + 锅炉补燃 >= 热负荷 + 制冷机热输入 Q_HR = param.eta_recovery * (param.a * P_MT + param.b); heat_excess = Q_HR + Q_GB - load_data.heat - Q_AC_in; % 允许弃热，但不能缺热 heat_vio = sum(max(0, -heat_excess)); % 冷平衡 Q_cool = param.COP_AC * Q_AC_in + param.COP_EC * P_EC; cool_vio = sum(max(0, load_data.cool - Q_cool)); % 储能SOC约束 SOC_vio = check_soc(x, param); % 单独函数处理 % 运行成本 cost_fuel = param.gas_price * ... (sum(param.a * P_MT + param.b) + sum(Q_GB / param.eta_GB)); cost_grid = sum(price_data.elec_buy .* P_buy) - ... sum(price_data.elec_sell .* P_sell); cost_om = sum(param.k_MT .* P_MT + param.k_EC .* P_EC + ... param.k_GB .* Q_GB); cost = cost_fuel + cost_grid + cost_om; % 罚函数 penalty = param.penalty_coef * (heat_vio + cool_vio + SOC_vio); cost = cost + penalty; end

写目标函数的时候，能量平衡约束怎么判断要特别想清楚。上面这个例子中热平衡我用了max(0, -heat_excess)来计算违反量，因为工程上确实存在“弃热”操作，也就是说产热量大于需求量是可以接受的，但是需求量大于产热量是绝对不容许的。如果用abs来算违反量，反而会把合理的弃热方案当成不可行解，影响优化方向。

冷平衡那边逻辑反过来：冷量不够就是违反约束，冷量多了算设备选型过大或者策略不够优。不同能量流的约束逻辑不一样，不能一概用abs。

罚系数的量级要和目标函数匹配合。假设正常运行成本是几万块钱的量级，罚系数至少要设到几百上千，才能把违反约束的方案压下去。如果罚系数设得太大，算法会躲开那些“只有轻微违反但潜力巨大”的区域。我的习惯是先跑一次不带约束的优化，看看目标函数的大致范围，再按这个范围定罚系数。

4.3 运行结果分析与可视化

跑完程序之后，我一般会输出三样东西：收敛曲线、各设备出力时序图、能量平衡校验表。

收敛曲线用来判断算法有没有收敛。正常情况下应该是一条前期快速下降、后期趋于平稳的曲线。如果曲线到后期还在大幅波动，说明参数没调好或者罚函数有问题；如果前期下降非常缓慢，大概率是初始化范围不对，或者边界映射写错了。

设备出力时序图是把24小时的最优出力画出来，直观反映调度策略是否合理。例如在分时电价高的时段，燃气轮机应该多发；电价低的时段，可以适当多买电、少自发电。储能设备应该是低谷充电、高峰放电。如果画出来的曲线不符合这个逻辑，即使目标函数值很好看，也要怀疑代码逻辑出了问题。

figure; subplot(2,2,1); plot(1:24, P_MT_opt, '-o', 'LineWidth', 1.5); xlabel('时段/h'); ylabel('功率/kW'); title('燃气轮机出力'); grid on; subplot(2,2,2); plot(1:24, P_buy_opt, '-s', 'LineWidth', 1.5); xlabel('时段/h'); ylabel('功率/kW'); title('电网购电'); grid on; subplot(2,2,3); plot(1:24, Q_AC_opt, '-^', 'LineWidth', 1.5); xlabel('时段/h'); ylabel('功率/kW'); title('吸收式制冷输入'); grid on; subplot(2,2,4); plot(1:24, P_EC_opt, '-d', 'LineWidth', 1.5); xlabel('时段/h'); ylabel('功率/kW'); title('电制冷出力'); grid on;

这一步看着简单，但真的能发现问题。有一次我画完图才发现，算法给出的购电曲线在凌晨电价低谷期是零，燃气轮机却满负荷运行。单纯看目标函数值是完全正常的，但画图一看就知道不对——凌晨电价那么低，应该买电而不是让设备满发。后来排查发现是维护成本系数设置不合理，导致算法宁愿多用设备也不买便宜的电。如果不画图，这个问题很难被察觉。

能量平衡校验表则是在优化结束后，把每个时段的电力、热力、冷力平衡情况列出来，人工核对一遍。我会专门写一个check_balance函数，输出每个时段的供给量和需求量差值，超过容忍阈值就告警。这套校验逻辑在代码交付时非常加分，能帮甲方快速建立信任。

5. 运行调参与踩坑记录：从新手到落地的实战经验

5.1 收敛慢、结果差，先排查这几个典型问题

我在多个项目里踩过的坑不少，挑几个典型的列出来，遇到了可以照着排查。

第一个坑是决策变量顺序不一致。目标函数里变量解析顺序和初始化时的定义一旦不一致，程序不会报错，但结果会非常离谱。比如决策变量是5类设备、每类24个变量，我习惯在代码开头用大块注释把变量排列顺序写清楚，每次改动都对照着查一遍。

第二个坑是储能SOC超出安全范围。SOC是一个递推量，如果充放电效率的符号写错，SOC会越走越偏，甚至跑到负值。我在递推后加了一个clamp函数强制把SOC限制在[0,1]内，同时在完成一次完整的24小时递推后，检查最终SOC是否回到初始值附近。如果偏离太远，说明储能策略存在问题，可能是充放电功率跨度设置太大。

第三个坑是罚函数的系数大小问题。我做过一个对比实验，目标函数值量级在10^4到10^5之间，罚系数取10的时候约束违反在结果里经常出现；取10^6之后，约束倒是满足了，但算法卡在可行域边缘找不到更优解。最后我的做法变成了两步走：先用较小的罚系数快速迭代200代，让个体进入可行域的大致区域，再改用较大的罚系数跑200代，把解拉到严格可行的范围内。这个方法比固定一个罚系数从头跑到尾收敛效果更好。

第四个坑是边界设置的问题。有些设备在特定工况下上限和下限非常接近，比如燃气轮机的技术最小出力是30%额定功率，如果优化程序里上限设成了额定功率、下限设成了0，算法就会在0到30%之间产生大量物理上不可行的解。建议所有设备的出力和约束都按真实物理下限设置。

5.2 算法参数怎么调最省心

POA需要设置的参数真的很少，但种群规模和迭代次数还是值得说两句。

种群规模我建议30到50之间。太小了，探索能力不足，容易陷入局部最优；太大了，每次迭代计算目标函数的次数太多。冷热电联供的目标函数本身不算复杂，但如果有储能递推、多设备耦合，几百个个体的计算量累积起来还是很可观的。我一般先取40，观察收敛曲线情况再增减。

迭代次数500到1000之间足够。POA收敛速度快，通常200代左右就能看到明显的平台期，再多迭代边际收益很低。画一下收敛曲线，看到曲线变平了就可以停了，没必要为了保险硬跑1000代。

还有一个容易被忽略的参数是边界映射方式。我在使用归一化决策变量时，映射回实际范围用的是：

x_actual = lb + x_norm × (ub - lb)

这里如果lb和ub设置得不好，比如某个时段燃气轮机的上限比下限还小，程序会直接崩溃或者产生一堆无效解。建议在数据初始化模块里对所有设备的上下限做一次校验，确保ub大于lb，并且把启动过程中的出力限制也算进去。

5.3 工程落地还想多说几句

如果这套代码只是做研究验证，前面提到的基本够了。但要真正落地到实际项目，还有两件事必须补上。

第一件事是模型校准。前面用到的设备模型都是简化版本，真实设备的热电比、部分负载效率、启停损耗都是非线性的。工程上建议拿一台真实机组的历史运行数据，对设备模型参数做最小二乘拟合，然后再交给算法使用。我接手过某个项目，厂家给的额定效率是38%，实际运行下来只有35%左右，如果用额定参数做优化，调度结果和实际运行偏差会非常大。

第二件事是鲁棒性验证。调度结果对负荷预测误差相当敏感。分时电价预测错了，或者第二天气温突变导致冷负荷暴涨，优化出来的方案就可能完全不可行。我一般会在求解完成后做一次敏感性分析，看关键设备出力在负荷波动正负10%范围内的变化情况，确认方案不会因为一点偏差就崩溃。

另外必须提醒一句：POA是单次求解算法，每次跑的结果会略有差异。我在正式出结果前，一般会把同一个算例连续跑5次，取目标函数值最小的那次，同时对比多次结果的差异程度。如果差异太大，说明算法没有充分收敛，需要适当增加种群规模或者迭代次数。

最后再分享一个小技巧：调试的时候别直接跑24小时整段优化，先跑3小时或者6小时的子问题，把每个变量的变化打印出来，配合手算验证约束是否满足。等子问题全部通过，再扩展到全时段。这个方法帮我省下的排查时间至少都按天算了。

从算法原理到系统建模，从代码架构到实际调参，这一套流程跑下来，我认为POA和冷热电联供调度的搭配是值得尝试的。相比传统粒子群，它的收敛速度和解的质量都有明显优势；相比精确求解方法，它又不受问题规模和非线性程度的束缚。如果你手里正好有一个CCHP优化项目，不妨用我这里提供的框架先跑通一版，再根据实际数据进行迭代调整。能源调度的魅力就在于，模型永远可以更精细，算法永远有改进空间，而收益隐藏在每个细节的打磨之中。