尧图网站开发 
简介面向能源动力与电气工程相关专业学生及研究人员围绕微型燃气轮机冷热电联供系统的日前优化调度问题展开。程序选取夏季典型日仅考虑电负荷与冷负荷两类需求由燃气轮机出力、吸收式制冷与电制冷协同供能并可向电网购电以全天运行费用最小为目标函数采用粒子群算法PSO求解最优调度方案。资源包共8个文件含4个m脚本、2个jpg结果图及2个asv自动保存文件所有代码逐行注释、结构简练便于复现与二次开发压缩包仅61KB轻量易用。目前已有253人学习浏览特别适合希望掌握PSO算法在综合能源系统优化中应用方法的入门及进阶学习者。通过运行源码与查看功率曲线图可直观理解冷热电联供系统的建模过程、负荷平衡约束处理及粒子群寻优细节是学习智慧能源优化调度的实用参考资料。1. 微型燃气轮机冷热电联供系统优化的关键矛盾微型燃气轮机冷热电联供系统的运行优化难在电、热、冷三股能量互相耦合调任何一个出力都会牵动另外两条平衡。很多团队初期按“以热定电”或“以电定热”的规则排程一旦叠加上峰谷电价、蓄电池和蓄热罐人工排出来的方案距离经济最优往往差出一大截。PSO算法把一天的设备出力编码成粒子位置靠粒子之间的信息共享迭代搜索不要求目标函数可导、不要求模型凸性几百行代码就能形成一套可用的求解框架。这篇按建模、编码、实现、调参、验证的顺序写适合调度工程师和刚转入综合能源方向的算法工程师直接上手。2. CCHP系统建模与PSO算法的对接方式2.1 能量流与设备组成微型燃气轮机联供系统的骨架是一条能量梯级利用链天然气在微燃机里燃烧发电排出的高温烟气和缸套水通过余热回收产生热水或蒸汽一部分直接供暖一部分驱动溴化锂吸收式机组制冷。在这个基础上系统还会配电制冷机、燃气锅炉、蓄电池和蓄热罐用来补冷差、补热差、做跨时段平移。整个模型里微燃机每发 1 kWh 电同时产生约 1.6 kWh 的燃料余热按 η_e0.30、η_r0.48 折算这个“发电和产热没法拆开”的特性是所有优化矛盾的根源。提示余热回收效率取的是相对燃料热值的系数别和微燃机排烟温度简单等同。不同工况下效率有偏差第一版模型用常数没问题后续精修就按制造商部分负荷曲线插值。下面是一组常见的小型系统参数后面的示例代码直接按这套数值写设备/参数数值说明微燃机额定功率100 kW发电上限下限按 20 kW发电效率 η_e0.30相对燃料低位热值余热回收效率 η_r0.48烟气缸套水回收吸收式制冷 COP1.3热驱动溴化锂机组电制冷 COP3.0压缩式制冷燃气锅炉效率0.90补热备用蓄电池容量100 kWhSOC 0.1~0.9蓄热罐容量60 kWh热水蓄能这些参数是工程里比较典型的取值实际项目要以设备铭牌和热力学核算为准。尤其微燃机的部分负荷效率在 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_costP_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。不等式约束可以分成设备上下限、爬坡和储能三类约束类型典型表达式设备出力区间不等式P_min ≤ P_mgt ≤ P_max0 ≤ Q_abs ≤ Q_abs_maxQ_ec ≤ Q_ec_max爬坡不等式|P_mgt[t] - P_mgt[t-1]| ≤ ramp_max电池 SOC不等式SOC_min ≤ SOC[t] ≤ SOC_maxSOC 连续性差方程SOC[t] SOC[t-1] - P_bt[t]·Δt / η_dch放电时等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_size40, max_iter200, w_max0.9, w_min0.4, c12.0, c22.0, seed0): 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_size40, max_iter300)跑完后 best_cost 会比规则调度低收敛曲线的形态通常是前 30 代快速下降、100 代以后趋于平缓。把 best_x 按四个 T 段拆分可以看到符合经济直觉的排程时段P_mgt(kW)P_grid(kW)Q_abs(kW)电池 P_bt(kW)储热 P_ht(kW)0:00 低谷21.328.60.0-14.2充电-6.8蓄热12:00 午段88.26.531.421.0放电12.6放热20:00 晚峰96.00.024.818.3放电9.2放热深夜电价低电池选择充电而不是让微燃机多发电午后冷负荷上来吸收式制冷承担主要冷量电制冷只补差额晚高峰微燃机接近满发电池和储热同时放电/放热来压低购电价。这些行为不是人工写进去的而是目标函数和价格信号共同逼出来的。4. PSO参数整定与收敛诊断的实战细节4.1 惯性权重w的线性衰减w 控制上一代速度有多少保留到下一代。w 接近 1 时粒子主要沿原方向飞探索范围大w 接近 0 时粒子迅速向 pbest/gbest 靠拢收敛快但容易早熟。工程中最常见的是线性衰减w w_max - (w_max - w_min) * it / max_iter对应代码里 w_max0.9、w_min0.4。对 CCHP 这种 96 维变量的问题衰减速度可以再慢一点比如把 0.4 改成 0.5代价是迭代次数要增加 50~100 代。如果发现收敛曲线后期还在大幅跳动优先把 w_min 调高而不是调低跳动能帮助逃离局部最优。4.2 学习因子、种群规模与迭代次数的搭配c1 让粒子倾向自己的历史最优c2 让粒子倾向群体最优。c1 过大时粒子各飞各的收敛慢c2 过大时群体最优像磁铁容易把粒子吸进局部极值。CCHP 调度目标里惩罚项占了一部分权重适应度面比较崎岖c1 和 c2 的差距不宜拉太大。给一张工程上常用的取值范围表参数推荐范围对结果的影响调试优先级pop_size40~80太小早熟太大单轮评估耗时成倍高max_iter200~500不足则收敛不完整中w0.9 → 0.4~0.5全局搜索和局部收敛的平衡中c11.5~2.0个人记忆权重低c21.8~2.5群体引导权重低v_max0.1~0.2×(ub-lb)过大震荡过小早熟中种群规模对单次 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_size40, max_iter150) 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-3P_grid 序列无负值Q_boiler 无负值电池和储热罐的 SOC 全程落在上下界内且终值与初值偏差在允许范围5 个随机种子下的最优成本极差不超过 1%结果符合基本经济逻辑谷时段电池充电、峰时段放电气价低于边际电价时微燃机多发。跑完这五项再把负荷数组换成另一天的预测重复一遍流程PSO 调度代码才算真正接进你的运行环境。本文还有配套的精品资源点击获取 
简介面向能源动力与电气工程相关专业学生及研究人员围绕微型燃气轮机冷热电联供系统的日前优化调度问题展开。程序选取夏季典型日仅考虑电负荷与冷负荷两类需求由燃气轮机出力、吸收式制冷与电制冷协同供能并可向电网购电以全天运行费用最小为目标函数采用粒子群算法PSO求解最优调度方案。资源包共8个文件含4个m脚本、2个jpg结果图及2个asv自动保存文件所有代码逐行注释、结构简练便于复现与二次开发压缩包仅61KB轻量易用。目前已有253人学习浏览特别适合希望掌握PSO算法在综合能源系统优化中应用方法的入门及进阶学习者。通过运行源码与查看功率曲线图可直观理解冷热电联供系统的建模过程、负荷平衡约束处理及粒子群寻优细节是学习智慧能源优化调度的实用参考资料。1. 微型燃气轮机冷热电联供系统优化的关键矛盾微型燃气轮机冷热电联供系统的运行优化难在电、热、冷三股能量互相耦合调任何一个出力都会牵动另外两条平衡。很多团队初期按“以热定电”或“以电定热”的规则排程一旦叠加上峰谷电价、蓄电池和蓄热罐人工排出来的方案距离经济最优往往差出一大截。PSO算法把一天的设备出力编码成粒子位置靠粒子之间的信息共享迭代搜索不要求目标函数可导、不要求模型凸性几百行代码就能形成一套可用的求解框架。这篇按建模、编码、实现、调参、验证的顺序写适合调度工程师和刚转入综合能源方向的算法工程师直接上手。2. CCHP系统建模与PSO算法的对接方式2.1 能量流与设备组成微型燃气轮机联供系统的骨架是一条能量梯级利用链天然气在微燃机里燃烧发电排出的高温烟气和缸套水通过余热回收产生热水或蒸汽一部分直接供暖一部分驱动溴化锂吸收式机组制冷。在这个基础上系统还会配电制冷机、燃气锅炉、蓄电池和蓄热罐用来补冷差、补热差、做跨时段平移。整个模型里微燃机每发 1 kWh 电同时产生约 1.6 kWh 的燃料余热按 η_e0.30、η_r0.48 折算这个“发电和产热没法拆开”的特性是所有优化矛盾的根源。提示余热回收效率取的是相对燃料热值的系数别和微燃机排烟温度简单等同。不同工况下效率有偏差第一版模型用常数没问题后续精修就按制造商部分负荷曲线插值。下面是一组常见的小型系统参数后面的示例代码直接按这套数值写设备/参数数值说明微燃机额定功率100 kW发电上限下限按 20 kW发电效率 η_e0.30相对燃料低位热值余热回收效率 η_r0.48烟气缸套水回收吸收式制冷 COP1.3热驱动溴化锂机组电制冷 COP3.0压缩式制冷燃气锅炉效率0.90补热备用蓄电池容量100 kWhSOC 0.1~0.9蓄热罐容量60 kWh热水蓄能这些参数是工程里比较典型的取值实际项目要以设备铭牌和热力学核算为准。尤其微燃机的部分负荷效率在 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_costP_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。不等式约束可以分成设备上下限、爬坡和储能三类约束类型典型表达式设备出力区间不等式P_min ≤ P_mgt ≤ P_max0 ≤ Q_abs ≤ Q_abs_maxQ_ec ≤ Q_ec_max爬坡不等式|P_mgt[t] - P_mgt[t-1]| ≤ ramp_max电池 SOC不等式SOC_min ≤ SOC[t] ≤ SOC_maxSOC 连续性差方程SOC[t] SOC[t-1] - P_bt[t]·Δt / η_dch放电时等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_size40, max_iter200, w_max0.9, w_min0.4, c12.0, c22.0, seed0): 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_size40, max_iter300)跑完后 best_cost 会比规则调度低收敛曲线的形态通常是前 30 代快速下降、100 代以后趋于平缓。把 best_x 按四个 T 段拆分可以看到符合经济直觉的排程时段P_mgt(kW)P_grid(kW)Q_abs(kW)电池 P_bt(kW)储热 P_ht(kW)0:00 低谷21.328.60.0-14.2充电-6.8蓄热12:00 午段88.26.531.421.0放电12.6放热20:00 晚峰96.00.024.818.3放电9.2放热深夜电价低电池选择充电而不是让微燃机多发电午后冷负荷上来吸收式制冷承担主要冷量电制冷只补差额晚高峰微燃机接近满发电池和储热同时放电/放热来压低购电价。这些行为不是人工写进去的而是目标函数和价格信号共同逼出来的。4. PSO参数整定与收敛诊断的实战细节4.1 惯性权重w的线性衰减w 控制上一代速度有多少保留到下一代。w 接近 1 时粒子主要沿原方向飞探索范围大w 接近 0 时粒子迅速向 pbest/gbest 靠拢收敛快但容易早熟。工程中最常见的是线性衰减w w_max - (w_max - w_min) * it / max_iter对应代码里 w_max0.9、w_min0.4。对 CCHP 这种 96 维变量的问题衰减速度可以再慢一点比如把 0.4 改成 0.5代价是迭代次数要增加 50~100 代。如果发现收敛曲线后期还在大幅跳动优先把 w_min 调高而不是调低跳动能帮助逃离局部最优。4.2 学习因子、种群规模与迭代次数的搭配c1 让粒子倾向自己的历史最优c2 让粒子倾向群体最优。c1 过大时粒子各飞各的收敛慢c2 过大时群体最优像磁铁容易把粒子吸进局部极值。CCHP 调度目标里惩罚项占了一部分权重适应度面比较崎岖c1 和 c2 的差距不宜拉太大。给一张工程上常用的取值范围表参数推荐范围对结果的影响调试优先级pop_size40~80太小早熟太大单轮评估耗时成倍高max_iter200~500不足则收敛不完整中w0.9 → 0.4~0.5全局搜索和局部收敛的平衡中c11.5~2.0个人记忆权重低c21.8~2.5群体引导权重低v_max0.1~0.2×(ub-lb)过大震荡过小早熟中种群规模对单次 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_size40, max_iter150) 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-3P_grid 序列无负值Q_boiler 无负值电池和储热罐的 SOC 全程落在上下界内且终值与初值偏差在允许范围5 个随机种子下的最优成本极差不超过 1%结果符合基本经济逻辑谷时段电池充电、峰时段放电气价低于边际电价时微燃机多发。跑完这五项再把负荷数组换成另一天的预测重复一遍流程PSO 调度代码才算真正接进你的运行环境。本文还有配套的精品资源点击获取