基于粒子群算法的冷热电联供系统优化调度Python实现

简介：面向能源动力与电气工程相关专业学生及研究人员，围绕微型燃气轮机冷热电联供系统的日前优化调度问题展开。程序选取夏季典型日，仅考虑电负荷与冷负荷两类需求，由燃气轮机出力、吸收式制冷与电制冷协同供能，并可向电网购电，以全天运行费用最小为目标函数，采用粒子群算法（PSO）求解最优调度方案。资源包共8个文件，含4个m脚本、2个jpg结果图及2个asv自动保存文件，所有代码逐行注释、结构简练，便于复现与二次开发，压缩包仅61KB，轻量易用。目前已有253人学习浏览，特别适合希望掌握PSO算法在综合能源系统优化中应用方法的入门及进阶学习者。通过运行源码与查看功率曲线图，可直观理解冷热电联供系统的建模过程、负荷平衡约束处理及粒子群寻优细节，是学习智慧能源优化调度的实用参考资料。

1. 微型燃气轮机冷热电联供系统优化的关键矛盾

微型燃气轮机冷热电联供系统的运行优化，难在电、热、冷三股能量互相耦合，调任何一个出力都会牵动另外两条平衡。很多团队初期按“以热定电”或“以电定热”的规则排程，一旦叠加上峰谷电价、蓄电池和蓄热罐，人工排出来的方案距离经济最优往往差出一大截。PSO算法把一天的设备出力编码成粒子位置，靠粒子之间的信息共享迭代搜索，不要求目标函数可导、不要求模型凸性，几百行代码就能形成一套可用的求解框架。这篇按建模、编码、实现、调参、验证的顺序写，适合调度工程师和刚转入综合能源方向的算法工程师直接上手。

2. CCHP系统建模与PSO算法的对接方式

2.1 能量流与设备组成

微型燃气轮机联供系统的骨架是一条能量梯级利用链：天然气在微燃机里燃烧发电，排出的高温烟气和缸套水通过余热回收产生热水或蒸汽，一部分直接供暖，一部分驱动溴化锂吸收式机组制冷。在这个基础上，系统还会配电制冷机、燃气锅炉、蓄电池和蓄热罐，用来补冷差、补热差、做跨时段平移。整个模型里，微燃机每发 1 kWh 电，同时产生约 1.6 kWh 的燃料余热（按 η_e=0.30、η_r=0.48 折算），这个“发电和产热没法拆开”的特性，是所有优化矛盾的根源。

提示：余热回收效率取的是相对燃料热值的系数，别和微燃机排烟温度简单等同。不同工况下效率有偏差，第一版模型用常数没问题，后续精修就按制造商部分负荷曲线插值。

下面是一组常见的小型系统参数，后面的示例代码直接按这套数值写：

这些参数是工程里比较典型的取值，实际项目要以设备铭牌和热力学核算为准。尤其微燃机的部分负荷效率，在 20% 负荷附近通常会比额定效率掉 3~5 个百分点，第一次做优化可以用常数，后续要精细再改成负荷的分段线性函数。

2.2 目标函数：运行成本怎么写

24 小时调度优化的目标，一般只写当天的可变成本，固定折旧和容量费不影响逐时排程，可以不进模型。成本项包括：微燃机燃料费和燃气锅炉补热燃料费，按天然气体积乘气价；从电网买电的费用，按分时电价结算，模型里不允许反送电时 P_grid 取 max(P_grid,0) 参与成本；微燃机、吸收式机组、电制冷机的运维费，按出力或制冷量乘一个小系数。

写成 numpy 风格的计算式很直观：

fuel_cost = (P_mgt / eta_e / LHV).sum() * gas_price boiler_cost = np.maximum(Q_boiler, 0).sum() / eff_boiler / LHV * gas_price grid_cost = np.maximum(P_grid, 0).dot(price) om_cost = 0.015 * P_mgt.sum() + 0.01 * Q_abs.sum() + 0.01 * Q_ec.sum() total_cost = fuel_cost + boiler_cost + grid_cost + om_cost

P_mgt 是每小时微燃机发电量，Q_boiler 是热平衡算出来的锅炉补热量，Q_abs/Q_ec 是吸收式和电制冷机的制冷量。LHV 按 9.7 kWh/Nm³ 折算，各气田的气价和热值不同，换项目时这两项必须校核。

2.3 三条能量平衡与一组不等式约束

模型的核心是三条小时级平衡：

• 电平衡：P_mgt + P_grid = E_load + Q_ec / COP_ec + P_bt
• 制冷平衡：Q_abs + Q_ec = C_load
• 热平衡：Q_rec + Q_boiler + P_ht = H_load + Q_abs / COP_abs

式中 P_bt 为蓄电池出力，正值放电、负值充电；P_ht 为蓄热罐出力，正值放热、负值蓄热；Q_rec 是余热回收热量，等于微燃机燃料输入乘 η_r。

不等式约束可以分成设备上下限、爬坡和储能三类：

PSO 对等式约束的处理，常见做法不是把残差写进罚函数硬压，而是在编码时做变量消元：制冷平衡直接解出 Q_ec = C_load - Q_abs；电平衡解出 P_grid；热平衡解出 Q_boiler。消完之后，真正需要惩罚的只剩“P_grid 不能小于 0、Q_boiler 不能小于 0、制冷量不能越限”这些不等式。

2.4 为什么选PSO而不是数学规划

CCHP 调度在数学上是混合整数非线性规划，效率曲线、储能累计和购售电切换都会带来非凸项。用商业求解器也能做，但工程团队没有许可证、换设备参数就要重新剪枝、再叠加重试和报表逻辑，维护成本不低。PSO 的优势在于它不关心目标函数是否光滑、是否凸，只要能把“一个解”编码成向量、能算出适应度，就可以开始迭代。实现上它对盒约束尤其友好：设备上下限直接变成粒子位置的上下界，其余约束进罚函数。

代价是结果有随机性，且不保证全局最优。因此行业内的常规操作是跑 10~20 次独立随机种子，取最优解，同时观察多次结果的标准差。如果多条收敛曲线峰值接近，说明解可信；如果每轮差异大，优先检查种群规模和迭代次数，而不是急着换算法。

3. 用Python实现PSO求解CCHP典型日优化

3.1 准备典型日负荷与分时电价

先把典型日电、热、冷负荷和分时电价写成数组。这里的负荷数据是示例值，工程中要替换为 SCADA 实测或日前预测序列。

import numpy as np T = 24 dt = 1.0 # 调度步长，单位小时 hours = np.arange(T) # 示例负荷曲线：工程中替换为实测或预测序列 E_load = 60 + 30 * np.sin(2 * np.pi * (hours - 7) / 24) # 电负荷 kW H_load = 50 * np.exp(-((hours - 9) ** 2) / 18) + 12 # 供暖热负荷 kW C_load = 45 * np.exp(-((hours - 14) ** 2) / 28) + 6 # 空调冷负荷 kW # 峰平谷三段电价，8-11 峰，18-21 峰，22-6 谷，其余平段 price = np.where((hours >= 8) & (hours <= 11), 0.98, np.where((hours >= 18) & (hours <= 21), 1.12, np.where((hours >= 22) | (hours <= 6), 0.32, 0.62)))

这里用电负荷按午前低、日落后小高峰的正弦形态近似，冷负荷峰值落在 14 时附近，热负荷峰值在 9 时前后。实际项目里这三个数组来自负荷预测接口，目标函数和 PSO 结构不用改。电价数组用了嵌套 np.where，多个条件之间必须加括号，否则 & 和 | 的优先级会出错。

3.2 决策变量编码与适应度函数

接着定义设备参数和适应度函数。决策变量只留了四组：P_mgt（微燃机发电）、Q_abs（吸收式制冷量）、P_bt（蓄电池出力）、P_ht（蓄热罐出力），其余变量通过平衡方程解出。

p_min_mgt, p_max_mgt = 20.0, 100.0 # kW q_max_abs = 80.0 p_max_ec = 80.0 p_bt_max = 40.0 p_ht_max = 30.0 e_bat_max = 100.0 h_tank_max = 60.0 ramp_max = 30.0 # kW/h 爬坡速率 eta_e, eta_r = 0.30, 0.48 cop_abs, cop_ec = 1.3, 3.0 eff_boiler = 0.90 eta_ch, eta_dch = 0.92, 0.95 # 电池充/放电效率 eta_ht_ch, eta_ht_dch = 0.95, 0.95 # 储热罐充/放热效率 LHV = 9.7 gas_price = 3.2 # 元/Nm3 soc0 = 0.5 def evaluate(x): P_mgt = np.clip(x[0:T], p_min_mgt, p_max_mgt) Q_abs = np.clip(x[T:2*T], 0.0, q_max_abs) P_bt = x[2*T:3*T] P_ht = x[3*T:4*T] P_fuel = P_mgt / eta_e Q_rec = P_fuel * eta_r Q_ec = C_load - Q_abs P_grid = E_load + Q_ec / cop_ec + P_bt - P_mgt Q_boiler = H_load + Q_abs / cop_abs - Q_rec - P_ht # 运行成本：微燃机燃气 + 锅炉补热 + 购电 + 运维 fuel_cost = (P_fuel / LHV).sum() * gas_price boiler_cost = np.maximum(Q_boiler, 0).sum() / eff_boiler / LHV * gas_price grid_cost = np.maximum(P_grid, 0).dot(price) om_cost = 0.015 * P_mgt.sum() + 0.01 * Q_abs.sum() + 0.01 * Q_ec.sum() cost = fuel_cost + boiler_cost + grid_cost + om_cost # 软约束惩罚 pen = 0.0 pen += 1e3 * np.sum(np.maximum(-P_grid, 0.0) ** 2) # 不允许反送电 pen += 1e3 * np.sum(np.maximum(-Q_boiler, 0.0) ** 2) # 热平衡供小于求 pen += 1e3 * np.sum(np.maximum(Q_ec - p_max_ec, 0) ** 2) pen += 1e3 * np.sum(np.maximum(-Q_ec, 0) ** 2) # 爬坡约束 dP = np.abs(np.diff(np.r_[P_mgt[0], P_mgt])) pen += 1e2 * np.sum(np.maximum(dP - ramp_max, 0) ** 2) # 蓄电池 SOC 滚动 soc = soc0 * e_bat_max for t in range(T): if P_bt[t] >= 0: soc -= P_bt[t] * dt / eta_dch # 放电 else: soc += (-P_bt[t]) * dt * eta_ch # 充电 if soc < 0.1 * e_bat_max or soc > 0.9 * e_bat_max: gap = min(soc - 0.1 * e_bat_max, 0.9 * e_bat_max - soc) pen += 1e4 * gap ** 2 pen += 1e2 * (soc - soc0 * e_bat_max) ** 2 # 终值回到初始附近 # 储热罐 SOC ss = 0.3 * h_tank_max for t in range(T): if P_ht[t] >= 0: ss -= P_ht[t] * dt / eta_ht_dch # 放热 else: ss += (-P_ht[t]) * dt * eta_ht_ch # 蓄热 if ss < 0.0 or ss > h_tank_max: gap = min(ss, h_tank_max - ss) pen += 1e4 * gap ** 2 pen += 1e2 * (ss - 0.3 * h_tank_max) ** 2 return cost + pen

这段代码的处理逻辑需要说明几点。第一，制冷平衡被直接消元了，Q_ec 由冷负荷减 Q_abs 得出，所以不需要额外惩罚 Q_abs + Q_ec = C_load 这个等式。第二，电平衡和热平衡分别解出 P_grid 和 Q_boiler，它们本身不参与罚函数，罚的是“解出来为负”也就是物理不可行的方向。第三，SOC 是逐时段累计的，不能只约束上下界，还要把终值拉回初始值附近，否则 PSO 会在最后一小时把电池放空来降低成本。

3.3 PSO迭代主循环

PSO 主循环按标准的速度-位置更新写，加上速度钳制和边界反射：

def pso(n_dim, bounds, obj_func, pop_size=40, max_iter=200, w_max=0.9, w_min=0.4, c1=2.0, c2=2.0, seed=0): rng = np.random.default_rng(seed) lb, ub = bounds[:, 0], bounds[:, 1] span = ub - lb pos = rng.uniform(lb, ub, size=(pop_size, n_dim)) vel = rng.uniform(-0.1 * span, 0.1 * span, size=(pop_size, n_dim)) pbest = pos.copy() pbest_val = np.array([obj_func(p) for p in pos]) g = np.argmin(pbest_val) gbest_val = pbest_val[g] gbest_pos = pbest[g].copy() best_hist = [gbest_val] for it in range(max_iter): w = w_max - (w_max - w_min) * it / max_iter r1 = rng.random((pop_size, n_dim)) r2 = rng.random((pop_size, n_dim)) vel = w * vel + c1 * r1 * (pbest - pos) + c2 * r2 * (gbest_pos - pos) v_max = 0.2 * span vel = np.clip(vel, -v_max, v_max) pos = pos + vel # 边界反射：不要只裁剪，反射保留更多探索性 for j in range(n_dim): low, up = lb[j], ub[j] pos[:, j] = np.where(pos[:, j] < low, 2 * low - pos[:, j], pos[:, j]) pos[:, j] = np.where(pos[:, j] > up, 2 * up - pos[:, j], pos[:, j]) pos = np.clip(pos, lb, ub) val = np.array([obj_func(p) for p in pos]) improved = val < pbest_val pbest[improved] = pos[improved] pbest_val[improved] = val[improved] g = np.argmin(pbest_val) if pbest_val[g] < gbest_val: gbest_val = pbest_val[g] gbest_pos = pbest[g].copy() best_hist.append(gbest_val) return gbest_pos, gbest_val, best_hist

速度更新用了三个可调量：惯性权重 w、个体学习因子 c1、群体学习因子 c2。w 从 0.9 线性降到 0.4，前期粒子飞得远、负责全局搜索，后期收缩做局部精修。c1 和 c2 都取 2.0 是老经验值，实际工程里对 CCHP 这种上百维的问题，c1 在 1.5~2.0、c2 在 1.8~2.5 的表现差别不大，后面第 4 章展开。注意我给 pso 函数加了 seed 参数，所有重复实验都固定随机种子，对比参数时才有可比性。

3.4 运行结果与调度策略解读

调用代码：

bounds = np.vstack([ np.tile([p_min_mgt, p_max_mgt], (T, 1)), np.tile([0.0, q_max_abs], (T, 1)), np.tile([-p_bt_max, p_bt_max], (T, 1)), np.tile([-p_ht_max, p_ht_max], (T, 1)), ]) best_x, best_cost, hist = pso(4 * T, bounds, evaluate, pop_size=40, max_iter=300)

跑完后 best_cost 会比规则调度低，收敛曲线的形态通常是前 30 代快速下降、100 代以后趋于平缓。把 best_x 按四个 T 段拆分，可以看到符合经济直觉的排程：

深夜电价低，电池选择充电而不是让微燃机多发电；午后冷负荷上来，吸收式制冷承担主要冷量，电制冷只补差额；晚高峰微燃机接近满发，电池和储热同时放电/放热来压低购电价。这些行为不是人工写进去的，而是目标函数和价格信号共同逼出来的。

4. PSO参数整定与收敛诊断的实战细节

4.1 惯性权重w的线性衰减

w 控制上一代速度有多少保留到下一代。w 接近 1 时粒子主要沿原方向飞，探索范围大；w 接近 0 时粒子迅速向 pbest/gbest 靠拢，收敛快但容易早熟。工程中最常见的是线性衰减：

w = w_max - (w_max - w_min) * it / max_iter

对应代码里 w_max=0.9、w_min=0.4。对 CCHP 这种 96 维变量的问题，衰减速度可以再慢一点，比如把 0.4 改成 0.5，代价是迭代次数要增加 50~100 代。如果发现收敛曲线后期还在大幅跳动，优先把 w_min 调高而不是调低，跳动能帮助逃离局部最优。

4.2 学习因子、种群规模与迭代次数的搭配

c1 让粒子倾向自己的历史最优，c2 让粒子倾向群体最优。c1 过大时粒子各飞各的，收敛慢；c2 过大时群体最优像磁铁，容易把粒子吸进局部极值。CCHP 调度目标里惩罚项占了一部分权重，适应度面比较崎岖，c1 和 c2 的差距不宜拉太大。给一张工程上常用的取值范围表：

种群规模对单次 eval 耗时的影响是线性的。evaluate 里做了 SOC 的逐时段循环，96 维变量乘 60 个粒子乘 300 代，纯 Python 大概要几十秒到一两分钟，可以接受。如果模型再加入机组启停等整数变量，单次 eval 变贵，建议先用 30 个粒子的试探性跑法调惩罚系数，最后再用 80 粒子出正式结果。

4.3 速度钳制、边界处理与早停判断

速度钳制限制每次位置跳动的幅度。v_max 取 0.2 倍变量跨度在多数问题里够用；当不同设备出力范围差异大时，要对每一维单独算 span，不能统一用同一个绝对值。

边界处理上，裁剪虽然简单，但会把一批粒子拍在边界上，种群多样性下降。反射处理把越过边界的粒子弹回可行域内，粒子保留向外的冲量，探索性更好。代码里用 np.where 做两次反射，最后再 clip，防止反射绕到另一侧。

早停判断要基于一段窗口而不是单次变化，否则粒子抖动会误触发。连续 10 代最优值相对变化小于 1e-4 时停，常见写法是：

if len(best_hist) > 60: delta = abs(best_hist[-1] - best_hist[-11]) if delta < 1e-4 * (abs(best_hist[-1]) + 1e-9): break

另外，正式结果不要只跑一次。用 5 个不同 seed 各跑一轮，观察最优成本的相对极差；极差超过 1% 就加迭代数或种群，这是判断 PSO 是否收敛到可信解最快的方法。

4.4 约束惩罚系数的自适应调节

罚函数方法最坑的是惩罚系数定多少。系数太小，粒子发现违反约束的成本更低，解会停在不可行域边缘；系数太大，罚项完全压过真实目标，粒子只关心“怎么不违规”，目标函数梯度信号被淹没。常见做法是先取 1e3 跑一轮，统计最优解里最大违规量，比如 SOC 越界 2 kWh、P_grid 出现 -5 kW，再把对应罚系数翻倍重跑。

如果懒得手动试，可以在 PSO 迭代过程中做简单自适应。每次迭代记录总体违规量 violation，用它调整罚权：

penalty = 1e3 while epoch % 50 == 0: if violation > 1e-3: penalty *= 1.5 else: penalty *= 0.95

这个方案的关键是把违规量和成本量纲对齐后再乘系数。SOC 的偏差单位是 kWh，能量不平的偏差单位也是 kWh，乘 1e3 相当于把每 kWh 违规定价为 1000 元，远高于正常电价和燃料成本，粒子就会优先把违规消除。注意等式约束的罚系数要比不等式高一个量级，因为等式残差不会像不等式那样自动归零。

5. 从单日优化到滚动优化与灵敏度分析

5.1 只执行第一时段的滚动优化

单日优化的解在预测准确时才可靠。到了实际运行，负荷预测和气温都在变，行业标准做法是滚动优化：保持 24 小时窗口不变，每一小时重新跑一次 PSO，只采用第一个时段的指令，下一时刻用更新的预测再优化。骨架代码如下：

horizon = 24 for now in range(0, 24 * 7, 1): load_forecast = get_forecast(now, horizon) x_opt, _, _ = pso(4 * horizon, bounds, evaluate, pop_size=40, max_iter=150) dispatch[now] = { 'P_mgt': x_opt[0], 'Q_abs': x_opt[T], 'P_bt': x_opt[2 * T], 'P_ht': x_opt[3 * T], }

滚动优化的好处是蓄电池 SOC 和历史误差会被下一轮优化自动修正，PSO 不用为预测偏差单独做反馈控制。代价是最优性差了一点，但换来了对扰动的鲁棒性。

5.2 电价和气价的灵敏度分析

用固定 seed 把天然气价格从 3.2 按 ±10% 步长扫一遍，观察总成本和微燃机日出力曲线。气价上涨时，PSO 会把部分发电量让给电网，同时把吸收式制冷向电制冷转移，这是模型对气价信号的最直接响应；电价低谷段充电、高峰段放电的动作则几乎不受气价影响。灵敏度分析的输出是一张“输入变化 % → 总成本变化 %”的表，比单点最优解更有决策价值。

5.3 三分钟验收清单

换任何一组新负荷数据，都按这套项目检查：

• 最优解处电、冷、热三条平衡残差小于 1e-3；
• P_grid 序列无负值，Q_boiler 无负值；
• 电池和储热罐的 SOC 全程落在上下界内，且终值与初值偏差在允许范围；
• 5 个随机种子下的最优成本极差不超过 1%；
• 结果符合基本经济逻辑：谷时段电池充电、峰时段放电，气价低于边际电价时微燃机多发。

跑完这五项，再把负荷数组换成另一天的预测，重复一遍流程，PSO 调度代码才算真正接进你的运行环境。

本文还有配套的精品资源，点击获取