做综合能源调度的人十有八九都被“冷热电怎么同时平衡”这个问题折磨过。燃气轮机发出来的电好安排但排出的余热用不完就是浪费追着热负荷跑电价尖峰时又得从电网高价买电。我去年在某个综合能源示范项目上做运行优化用Matlab把一整套冷热电联供CCHP系统模型搭起来再用粒子群算法PSO做24小时运行调度才算是把“电、热、冷”这三条能源总线理顺了。这篇东西就是这套方法的一个完整复盘从问题出发到数学模型、算法设计、代码实现再到踩坑记录适合正在做综合能源系统优化、或者想用智能算法做能量调度的朋友参考。1. 为什么冷热电联供系统需要运行优化1.1 系统构成与痛点冷热电联供系统的核心思路很简单让燃料先走一遍发电设备把发出来的电送给电负荷同时把发电过程产生的余热回收起来一部分直接供热另一部分送进吸收式制冷机变成冷量。听起来很完美但实际跑起来全是“互相撞车”的矛盾。以我做的那个园区项目为例系统里主要设备包括燃气轮机、余热锅炉或者叫余热回收装置、燃气锅炉、溴化锂吸收式制冷机组、电制冷机。电力母线旁边还配置了一套蓄能装置可以在谷电时段充电在峰电时段放电。这套系统在夏季典型日的最大问题是中午冷负荷飙到峰值但热负荷几乎为零燃气轮机若是按电负荷运行余热回收的热量根本驱动不了满负荷的吸收式制冷机如果反过来让燃气轮机去追热负荷电负荷缺口又得靠电网高价补。这就是典型的“跟随电负荷”和“跟随热负荷”策略都不好用的场景。那时候我需要解决的就是给出一套每天96点或者24点的调度方案每一台设备每个时段出力多少电网买多少电蓄能装置什么时候充、什么时候放让全天的运行成本最低、碳排放尽量少同时还得满足所有能量平衡和设备出力限制。1.2 运行优化的目标与约束运行优化本质上是一个“有约束的非线性优化问题”。目标函数可以根据项目需求选通常有三种写法经济性目标购电费用加燃气费用再减去可能的售电收入环保性目标全周期碳排放量最低包括购电折算碳排放和燃气直接碳排放综合目标把经济指标和环境指标按权重合并成一个单目标或者干脆做多目标帕累托优化。约束条件则更麻烦。首先要保证电、热、冷三条母线的能量守恒其次每个设备都有出力上下限和爬坡限制蓄能装置有容量上下限、充放电速率限制还要满足调度周期开始时和结束时的蓄能量相等。这些约束不是随便写写就完事的稍微漏掉一个等式约束优化结果就可能在运行层面完全不可行。1.3 为什么选择粒子群算法在选定算法之前我其实也对比过传统线性规划和遗传算法。CCHP系统里燃气轮机效率随部分负荷率变化、吸收式制冷机COP随输入热参数变化这些非线性环节如果硬用线性规划建模阶段就得做大量线性化近似模型复杂度上去了精度还不一定保得住。遗传算法能用但算子多、参数多调起来反而费劲。粒子群算法吸引我的点在于实现代码量小只需要维护速度和位置两个向量没有复杂的交叉变异算子对非线性、非凸的目标函数适应能力强参数就那么几个惯性权重、学习因子、种群规模工程上很容易摸到经验范围。用粒子群做系统运行优化一个重要价值是“快速试错”负荷曲线变了、电价方案变了、设备容量变了改几个参数再跑一轮就能得到参考调度方案这对前期方案比选非常有帮助。2. 系统数学模型搭建——优化的第一步2.1 设备建模很多初学者上来就把粒子群代码写得飞起结果目标函数都是拍脑袋糊的跑出来的结果自然没法看。建模这一步看着简单实际上最耗时间。燃气轮机模型我用了工程上常见的简化公式燃料消耗与发电出力的关系F_gt a * P_gt b * P_gt_rated其中P_gt是当前发电功率P_gt_rated是额定功率a和b是我根据燃气轮机效率曲线拟合出来的系数。如果只要粗略估算也可以用F_gt P_gt / η_gtη_gt取额定电效率0.35左右但要注意部分负荷下效率会掉。余热回收热功率H_gt F_gt * η_recη_rec取0.45~0.5表示燃料余热中能被回收利用的比例。溴化锂吸收式机组核心参数是制冷性能系数COP一般取1.2~1.4。它消耗的不是电而是高温热水或者蒸汽。如果吸收式机组输出冷量Q_abs那它需要从余热/锅炉侧取走的热量大约是Q_abs / COP_abs。这个量要写进热平衡式子里很多新手容易漏。电制冷机就简单得多输入电功率P_ec输出冷量Q_ec COP_ec * P_ecCOP电制冷机取3.5~4.5。燃气锅炉作为辅助热源热效率取0.9左右H_boiler η_boiler * F_fuel_boiler。蓄能装置如果是蓄电或者蓄冷蓄热需要单独写状态转移方程SOC(t1) SOC(t) η_ch * P_ch(t) - P_disch(t) / η_dis。这个方程会把调度周期内各时段的状态串起来让问题的求解难度提升一个档次。2.2 目标函数经济性目标函数看起来长其实逻辑很直白min Cost Σ [电价(t) * 购电功率(t) * Δt 天然气价 * 燃气功率(t) * Δt]如果系统允许并网售电那购电功率为负时就可以记为售电收益。为了让成本模型贴近真实分时电价一定要做进去。我用的典型峰平谷电价是峰时1.1元/kWh、平时0.7元/kWh、谷时0.4元/kWh。碳排放目标可以写成E Σ [电网排放因子(t) * 购电功率(t) 天然气排放因子 * 燃气功率(t)]如果要单目标优化我习惯把经济性成本和碳排放成本加权C_total C λ * Eλ就是环境代价折算系数。λ取0的时候是纯经济调度λ调大之后算法会主动减少高排放时段的外购电。2.3 约束条件约束条件必须写全否则优化结果看起来成本很低实际根本跑不出那样的运行状态。我总结下来至少要包含这几类电功率平衡购电 燃气轮机发电 电负荷 电制冷耗电 蓄能充电功率放电时蓄能功率取负热功率平衡燃气轮机余热 锅炉供热 蓄热放热 热负荷 吸收式制冷耗热冷功率平衡吸收式制冷 电制冷 蓄冷放冷 冷负荷设备出力上下限燃气轮机、锅炉、吸收式、电制冷都不能超过容量爬坡约束燃气轮机相邻时段输出变化不能超过某个限值蓄能设备SOC边界和首尾相等约束。等式约束在粒子群算法里不好直接硬处理。我采用的方案是罚函数法把平衡等式两侧的差值平方后乘以一个罚系数加进目标函数里。这个罚系数不能太小太小的话最后的调度方案会存在能量对不上的点也不能太大太大会让粒子群的搜索方向被罚函数主导。实战里我一般从1000起步跑几十次观察结果。2.4 典型算例参数设置为了让这个复盘有可复现性我把当时用的参数大致列一下读者可以根据自己的系统改设备参数数值燃气轮机额定电功率1000 kW燃气轮机额定电效率0.35余热回收余热回收效率0.45燃气锅炉额定热功率800 kW燃气锅炉热效率0.90溴化锂机组制冷容量500 kW溴化锂机组COP1.2电制冷机制冷容量800 kW电制冷机COP4.0蓄能装置容量300 kWh电网购电峰/平/谷电价1.1 / 0.7 / 0.4 元/kWh天然气单位热值价格0.35 元/kWh典型日负荷我采用的是某园区夏季数据电负荷最高大约1200 kW热负荷最高大约600 kW冷负荷最高大约900 kW。三条负荷曲线明显错峰这样才能看出优化调度的价值。3. 粒子群算法设计与改进3.1 标准粒子群原理粒子群算法的生物学背景是老生常谈鸟群找食物时每只鸟既记住自己找到的最好位置也参考整个鸟群发现的最好位置。对应到优化问题里每个粒子就是一组候选解粒子有自己的位置向量X和速度向量V。标准更新公式为V_new w * V c1 * r1 * (pbest - X) c2 * r2 * (gbest - X)X_new X V_new其中w是惯性权重控制粒子继承上一时刻速度的程度c1、c2是学习因子r1、r2是0到1之间的随机数。如果w过大粒子会飞过头难以收敛w过小容易早熟收敛到局部最优。工程上常用的做法是让w从0.9线性递减到0.4。3.2 编码方式我一开始犯了个错误把每个设备每个时段的出力都塞进粒子向量结果维度暴涨。后来做了简化决策变量只取购电功率、燃气轮机发电功率、燃气锅炉供热功率、溴化锂机组输出冷量、电制冷机耗电功率、蓄能装置充放电功率这6类变量。调度周期是24小时所以每个粒子就是一个6×24144维的向量。每个维度的上下界直接由设备容量决定。比如说燃气轮机发电功率设置在0~1000 kW电制冷耗电功率设置在0~200 kW因为800 kW制冷容量对应COP4蓄能装置充放电则用正负号表示充/放数值限制在对应速率范围内。这种编码方式很直接但有个隐患粒子在边界附近更新时很容易超出设备允许范围。所以我每次更新完位置后都会做一次边界修正把超出上限的拉回上限低于下限的拉回下限。速度也要做钳制防止粒子横冲直撞。3.3 约束处理方式罚函数法是粒子群处理约束的首选因为代码好写而且适用于各种复杂约束。具体做法是构造一个适应度函数fitness 目标函数值 ρ * 惩罚项对等式约束比如电功率不平衡量eq_e(t)惩罚项加上ρ * Σ eq_e(t)^2对不等式约束比如设备出力越界g 0惩罚项加上ρ * Σ max(0, g)^2。我个人的经验是使用动态罚系数迭代前期ρ取小一点让算法有足够空间探索可行域迭代后期ρ逐渐加大让粒子被迫回到可行解附近。动态系数可以写成ρ ρ_0 * (0.1 iter / maxIter)或者更简单一点每迭代10次把ρ乘以1.5。这样做的好处是既不容易一开始就陷入不可行区域又能在后期保证约束严格满足。3.4 改进策略标准粒子群在低维问题上很好用但144维这种中高维问题容易有两个毛病收敛慢、容易掉进局部最优。我在实现时做了三处小改进效果比较明显。第一惯性权重线性递减。w从0.9降到0.4前前期重探索后阶段重开发。第二增加随机重置机制。每次迭代后随机选两个粒子把它们的位置重新初始化到搜索空间内。这样等于每隔一会儿就往群体里注入新鲜血液有效防止种群停滞。第三对全局最优解加微小扰动。用gbest 高斯扰动生成一个新候选解和目标函数比较如果更好就替换gbest。这个其实就是一种局部搜索让算法在收敛末段还有机会跳出小坑。4. Matlab实现与核心代码4.1 代码整体框架Matlab代码按功能分成四块参数初始化、目标函数、粒子群迭代主循环、结果绘图。如果只是快速验证方法不用搞复杂的文件结构全部写在一个m文件里也没问题。但为了后续扩展我倾向于把参数和目标函数独立成子函数。这里给出最核心的框架读者可以直接抄下来改参数。4.2 参数初始化代码%% 参数初始化 clc; clear; close all; rng(42); % 调度参数 T 24; % 调度时段数 dt 1; % 单位时段小时 nVar 6 * T; % 决策变量维度6类设备 * 24小时 % 粒子群参数 nP 60; % 粒子数 maxIter 200; % 最大迭代次数 wStart 0.9; % 初始惯性权重 wEnd 0.4; % 结束惯性权重 c1 2; % 个体学习因子 c2 2; % 全局学习因子 % 决策变量上下限购电、燃气轮机、锅炉、吸收式、电制冷、蓄能 lb zeros(1, nVar); ub zeros(1, nVar); lb(1:T) 0; ub(1:T) 1200; % 购电功率 lb(T1:2*T) 0; ub(T1:2*T) 1000; % 燃气轮机电功率 lb(2*T1:3*T) 0; ub(2*T1:3*T) 800; % 锅炉热功率 lb(3*T1:4*T) 0; ub(3*T1:4*T) 500; % 溴化锂机组输出冷量 lb(4*T1:5*T) 0; ub(4*T1:5*T) 200; % 电制冷耗电功率 lb(5*T1:6*T) -200; ub(5*T1:6*T) 200; % 蓄能功率放电为正 % 初始化粒子位置与速度 X repmat(lb, nP, 1) rand(nP, nVar) .* (repmat(ub-lb, nP, 1)); V zeros(nP, nVar); pbest X; pbest_fit inf(nP, 1); gbest zeros(1, nVar); gbest_fit inf;注意这里的购电上限我设成了1200 kW只是演示用实际要按园区变压器容量和电网合同设定。电制冷耗电功率上限200 kW是“制冷容量800 kW / COP 4”倒推出来的这样代入目标函数时才能对应到800 kW冷量。4.3 目标函数与罚函数实现目标函数文件是整套代码的核心我把它单独写成一个子函数。为了简化这里假设蓄能装置是电蓄能热蓄能和冷蓄能暂时不建模读者需要根据自己系统扩展。function cost objFun(x, params) T params.T; dt params.dt; % 拆出决策变量 P_grid x(1:T); P_gt x(T1:2*T); H_boiler x(2*T1:3*T); Q_abs x(3*T1:4*T); P_ec x(4*T1:5*T); P_es x(5*T1:6*T); % 正为放电负为充电 % 燃气轮机运行特性 eta_gt 0.35; eta_rec 0.45; F_gt P_gt / eta_gt; % 输入燃料热功率 H_gt F_gt * eta_rec; % 回收余热 % 电制冷输出冷量 Q_ec 4.0 * P_ec; % 基础运行成本 fuelCost params.gasPrice * sum(F_gt H_boiler / 0.9) * dt; elecCost sum(params.ePrice .* P_grid) * dt; baseCost fuelCost elecCost; % 等式约束——电平衡 eqE P_grid P_gt - params.eLoad - P_ec - P_es; % 等式约束——热平衡吸收式制冷消耗热量 Q_abs/COP_abs eqH H_gt H_boiler - params.hLoad - Q_abs / 1.2; % 等式约束——冷平衡 eqC Q_abs Q_ec - params.cLoad; % 不等式约束蓄能SOC首尾相等可以用惩罚项 % 这里简单用蓄能累计量近似完整模型需要SOC状态转移 eqSOC sum(P_es); % 理想情况下一天累计净放电为0 % 罚函数 penalty 1e4 * sum(eqE.^2) 1e4 * sum(eqH.^2) 1e4 * sum(eqC.^2) 1e3 * (eqSOC)^2; cost baseCost penalty; end这个目标函数写得比较直接。等号约束的罚系数统一取1e4SOC首尾约束取1e3。实际应用时罚系数需要反复试我见过有人用1e6甚至更高但如果罚系数过大粒子群在探索阶段会非常僵硬所有粒子都挤在可行域边缘收敛质量反而下降。4.4 粒子群主循环主循环逻辑不复杂关键是把速度和位置更新写对并且每次更新后重新评估适应度。%% 粒子群主循环 fitness zeros(nP, 1); for iter 1:maxIter w wStart - (wStart - wEnd) * iter / maxIter; for i 1:nP % 越界处理把位置限制在边界内 X(i, :) max(X(i, :), lb); X(i, :) min(X(i, :), ub); % 计算适应度 fitness(i) objFun(X(i, :), params); % 更新个体最优 if fitness(i) pbest_fit(i) pbest(i, :) X(i, :); pbest_fit(i) fitness(i); end % 更新全局最优 if fitness(i) gbest_fit gbest_fit fitness(i); gbest X(i, :); end end % 速度与位置更新 r1 rand(nP, nVar); r2 rand(nP, nVar); V w * V ... c1 * r1 .* (pbest - X) ... c2 * r2 .* (repmat(gbest, nP, 1) - X); % 速度钳制 V max(V, -200); V min(V, 200); X X V; % 随机重置两个粒子 resetIdx randperm(nP, 2); X(resetIdx, :) repmat(lb, 2, 1) rand(2, nVar) .* repmat(ub-lb, 2, 1); fprintf(Iter %d, best cost %.2f\n, iter, gbest_fit); end速度钳制上限我随意写了200实际要根据变量量纲去调整。如果你的决策变量数量级在0~1000速度上限设为变量范围的10%~20%比较合理。4.5 结果分析与绘图迭代结束后把gbest拆回各设备出力就能画出24小时调度图。我一般会画三张图第一张是电力平衡图把购电、燃气轮机发电、电制冷耗电、蓄能充放电、电负荷画成堆叠图看电功率在每一小时到底平不平衡。第二张是冷功率平衡图重点看吸收式机组和电制冷机的出力占比。价格高的时候算法通常会让电制冷少出力尽量用吸收式制冷。第三张是收敛曲线图横轴是迭代次数纵轴是目标函数值。我看这条曲线的主要目的是判断迭代次数是否够了。如果200代之后曲线还在明显下降说明设定得太短要加迭代数。5. 仿真结果与讨论5.1 不同运行策略的对比为了验证优化的价值我把三种策略做了对比跟随电负荷、跟随热负荷、粒子群优化调度。表格里的数据是当时那个算例的结果不同系统参数下数字会不同但趋势基本一致。策略日运行成本元日碳排放kg说明跟随电负荷2435015200余热过剩锅炉几乎不用冷负荷靠电制冷硬扛跟随热负荷2612016900余热不足时段要锅炉补燃且电负荷缺口大PSO优化调度2178013840电、热、冷同步平衡低谷充电高峰放电以这个算例来看优化调度比单纯跟随电负荷节省了约10.6%的日运行成本碳排放降低了约9.0%。成本节省主要来自两个方向一是让燃气轮机在高电价时段多发电同时把余热送去制冷二是在谷电时段给蓄能充电峰电时段放电赚取峰谷价差。5.2 改进粒子群与标准粒子群的收敛性标准粒子群用同样参数跑同一算例大概要60代左右才稳定而且中间会出现明显的平台期像是在局部最优附近转圈。加了惯性权重线性递减、随机重置和gbest扰动之后收敛速度明显变快大约30代就进入稳定区域。这个对比说明了一个问题对于144维的优化问题粒子群本身的“基础版本”并不是万能的必须针对维度高、约束多的特点做改进。如果什么都不改就硬跑结果大概率不够理想。5.3 参数敏感性分析我还做了一组简单的敏感性测试调整目标函数里的环境权重λ。当λ从0逐渐增大时系统外购电量明显减少燃气轮机发电时长增加即使某些时段燃气发电的经济性不如直接买电算法也会为了降低碳排放而选择燃气。这说明在综合能源系统优化里目标函数权重本质上决定了系统的运行偏好。另外天然气的价格变化也非常敏感。气价每上涨10%优化结果里的锅炉补燃量会明显下降燃气轮机会更倾向在效率最高点附近运行。做项目汇报时这几组敏感性曲线常常比最终调度方案更有说服力因为它解释的是“为什么这套方案经济性好”。6. 常见问题与调参心得6.1 结果不收敛或收敛很慢这个问题我遇到得太多了。最常见的原因是粒子数太少或者迭代次数不够。144维的问题用20个粒子根本不够种群至少要放到50以上。其次惯性权重的起始值不要设得太小0.9起步比较稳。如果100代之后曲线还在斜坡式下降大概率是粒子数不足不是迭代次数不够。还有一种隐蔽的情况目标函数计算里有bug比如变量拆错位置导致相邻维度的数据互相串位。遇到结果完全没规律时先把决策变量维度打印出来人工核对一遍。6.2 等式约束最后仍然对不上罚函数法做等式约束最怕罚系数不够大导致结果里电平衡差值还有几十千瓦。我当时排查的步骤是把所有时段的平衡残差单独打印出来哪几个小时差得多就去看那几小时的负荷和出力数据找出是哪个变量被边界卡住了。如果要求严格等式平衡可以在粒子群迭代完成后加一步“可行化修正”保持设备出力不变把购电功率作为松弛变量直接由电平衡反算出来。这个修正能让最终结果在工程上是严格等价的。6.3 陷入局部最优很多情况下粒子群最后给出的方案比简单跟随策略好但和更精细的混合整数规划结果比还有差距。这未必是算法不行可能只是没加跳出机制。我后来在迭代中加入了变异操作每轮随机重置两个粒子效果立竿见影。如果还不行可以考虑用多个初始种群分别跑再从中选最优。6.4 运行时间太长144维、60个粒子、200次迭代在普通办公电脑上用Matlab跑大约一到两分钟。如果模型更复杂运行时间会显著增加。我常用的办法是先用24小时96个时段验证逻辑最后再用1小时分辨率跑正式结果。粒子数的选择不要盲目加大60到80在大多数场景下够用。再往上就是维度本身的问题可以改用分时段滚动优化。6.5 调参速查表下面这张表是我从多次调试里总结出来的经验范围不一定精确但能帮新人少走很多弯路。参数建议范围我的调试心得粒子数50~100维度高就取100维度低取50足够迭代次数150~300看收敛曲线稳定后再加50代确认惯性权重0.9递减至0.4比固定值效果好学习因子1.5~2.5c1、c2相等即可差异太大易震荡速度上限变量范围的10%~20%防止粒子飞散罚系数1e4~1e6先小后大采用动态罚函数这门技术我后来又扩展了很多方向比如加上蓄冷罐、改成多目标帕累托、用混合算法做冷热电联供和地源热泵的协同调度。回头看整个项目最值钱的其实不是粒子群算法本身而是把设备模型和约束条件一条条捋清楚的过程。有一次跑出结果蓄能装置自动在谷电时段充电、峰电时段放电曲线和分时电价严丝合缝那一刻是真的有成就感。如果非要说经验我会建议先在24小时两设备的小系统上把代码跑通再逐步加复杂度千万不要一上来就搞超大维度。希望这篇复盘对你能有点帮助。