冷热电联供系统经济优化MATLAB建模与求解全解析

做CCHP冷热电联供系统经济优化运行，最怕的其实不是模型复杂，而是你辛辛苦苦用MATLAB算出一个“全局最优解”，拿给运行工程师一看，对方轻飘飘一句“这跟我手调的结果差不多，程序还慢半拍”，瞬间就把你噎住。我在综合能源方向折腾了几年，从单台设备建模做到多能源系统优化，最大的体会是：CCHP的程序难点根本不在求解器，而在怎么把冷、热、电三种能量流在设备层耦合清楚，再把这个耦合关系翻译成优化模型。这篇就围绕我自己维护的一套冷热电联供经济优化MATLAB程序，把模型怎么搭、约束怎么设、代码怎么组织、算例结果怎么验证，一次说透。

这套程序面向的是典型的园区综合能源站场景，做的是日前经济调度：给定未来24小时的冷负荷、热负荷、电负荷预测曲线，给定分时电价和天然气价格，优化燃气轮机、余热回收、吸收式制冷、电制冷、燃气锅炉以及储能的逐时出力，目标是让系统全天运行成本最低。程序支持纯线性规划LP和混合整数线性规划MILP两种模式，模型文件、数据文件、结果绘图分开，改参数不需要动主程序。无论你是做毕业设计、科研复现，还是实际项目的前期可研，都可以直接拿这个框架改。

1. 冷热电联供经济优化到底在优化什么

1.1 为什么“多能互补”经常算错账

很多人一听“冷热电联供”，第一反应是燃气轮机发电、余热制冷制热，效率高还环保，那优化不就是让设备尽量满发吗？这是最典型的误区。CCHP的“经济优化”不是让单台设备最经济，而是让整个系统在一天24小时里，每一时刻的电、热、冷都平衡，同时总成本最低。这里的成本包括买电的钱、买气的钱，还要扣除可能的售电收益。设备之间是强耦合的：燃气轮机发一度电，同时会带出一股余热；这股余热可以拿去制热，也可以驱动吸收式制冷，但热和冷在时序上未必刚好匹配。比如夏天中午冷负荷高、热负荷低，余热多出来怎么办？是开到吸收式制冷去，还是直接排掉，甚至为了平衡热而降低发电出力？这些取舍才是优化的价值所在。

另外一个容易算错账的地方是分时电价。燃气轮机发一度电的燃料成本通常比电网购电价低，但低多少取决于气价、发电效率、当前时段电价。低谷电价可能只要0.35元/度，而天然气发一度电折算下来可能要0.55元/度，这时候让燃气轮机满发反而亏钱。所以正确做法是：低谷多买电、少发电，高峰多发、少买电。这个“削峰填谷”的节奏，必须由优化模型自己通过电价信号算出来，不是拍脑袋定一个“高峰满发、低谷停机”策略就能覆盖所有场景的。

1.2 程序定位：日前经济调度的典型日优化

这套程序解决的是“给定未来24小时预测数据，预先制定每个设备逐时出力计划”的问题，也就是日前调度。时间步长取1小时，优化周期24个点，决策变量包括每小时的购电功率、燃气轮机发电功率、燃气锅炉产热功率、电制冷功率、吸收式制冷功率、储能充电/放电功率、储能蓄热/放热功率，还有可能的启停状态0-1变量。目标函数是全天总运行费用，包括购电费、购气费，售电收益作为负项进入目标。

程序定位为“典型日优化”，意味着它适合做方案比选、经济性分析和调度策略研究。如果要投入实际运行，还需要在日前计划的基础上加日内滚动修正，用更短的预测窗口校正负荷误差，这个我在第六章里展开讲。但哪怕只做典型日，模型的约束设置、求解器配置和结果校验方法都是通用的，先把这套基本功打扎实，后面加滚动优化只是多包一层循环的事。

2. 设备级建模：先把能量流拆干净

2.1 燃气轮机与余热回收的耦合关系

燃气轮机（或燃气内燃机）是CCHP的“心脏”，它具有发电和产热双重特性。建模时最常用的简化模型是：

• 发电效率η_gt，一般取0.28~0.35；
• 余热热电比β，表示发1 kW电同时能回收多少kW热，小型机一般取1.0~1.5；
• 发电功率P_gt(t)，出力范围通常在额定功率的30%~100%之间；
• 天然气耗量V_gt(t) = P_gt(t) / (η_gt × LHV)，其中LHV是天然气低位热值，取9.7 kWh/Nm³左右；
• 回收余热Q_hr(t) = β × P_gt(t)。

注意，这里的余热回收量是“回收”后的热，不是烟气里的全部热量。实际系统中烟气还要经过换热器、余热锅炉等设备，换热效率和介质温度都会影响可回收热量。做方案阶段用β简化没问题，但如果你已经有厂家设备样本，最好按样本上的“余热回收功率-发电功率”曲线分段线性化，这样结果更可信。

这里有一个非常关键的约束：燃气轮机的发电量和余热产出是同比例的，不能“只发电不要热”。如果某个时段热负荷很低，余热用不完，模型必须允许“排热”变量（Q_waste≥0），否则会出现为了满足热平衡被迫降发电出力、从而多买高价电的反直觉结果。程序里我专门加了Q_waste，并且在结果里绘制余热利用率的曲线，方便检查是否存在大量排热。

2.2 制冷侧的“两条腿”：吸收式与电制冷如何取舍

CCHP系统的冷负荷由两条通道提供：一是烟气/缸套余热驱动的吸收式制冷机，二是电网或燃气轮机发电驱动的电压缩式制冷机（电制冷）。

吸收式制冷机的建模很简单：Q_cool_ac(t) = COP_ac × Q_hr_in(t)，COP_ac一般取1.1~1.4。注意这个概念和电制冷COP的差别，电制冷COP动辄4~5，吸收式只有1.3左右，看起来效率很低。但吸收式用的是“本来可能被排掉的余热”，这部分热量的成本几乎为零，所以经济上仍然划算。电制冷COP_ec取3.5~5.0，效率高，但它消耗的是高位电能，运行成本取决于购电价或燃气轮机发电成本。

制冷侧的优化逻辑其实是一场“热量如何分配到热负荷与冷负荷”的博弈。当一天里冷负荷大、热负荷小时，尽量让吸收式制冷多消纳余热；当热负荷本身就很大时，余热优先供热，冷负荷由电制冷补足。这种分配关系用人工经验很难逐时排准，必须在模型里把热平衡、冷平衡耦合起来让求解器去算。

2.3 储能环节：建与不建的边界

储能是CCHP系统里最容易被高估收益的设备。我的程序里支持储电和储热两类，建模采用通用形式：

• SOC(t) = SOC(t-1) + η_ch × P_ch(t) × Δt − P_dis(t) × Δt / η_dis
• 0 ≤ SOC(t) ≤ E_cap
• 0 ≤ P_ch(t) ≤ P_ch_max，0 ≤ P_dis(t) ≤ P_dis_max
• SOC(0) = SOC(24) = 初始设定值

储电用于“低谷充电、高峰放电”的电价套利，储热用于把余热或低谷电产生的热转移到热负荷高峰。但储能本身有效率损耗，每一轮“充放循环”要付出至少5%~10%的能量代价。只有当峰谷电价差足够大、且充放电功率足够覆盖负荷波动时，储能才有经济价值。我在程序里默认把储电容量设成1000 kWh、最大充放功率300 kW，储热容量2000 kWh、最大充放功率400 kW，但算完常常发现储电一天下来几乎不动。这不是程序bug，而是参数和电价结构决定的，后面第五章我会专门讲这个现象。

3. 经济优化模型的核心结构与约束设计

3.1 目标函数：购电、购气、售电三项怎么权衡

经济优化目标函数我写成：

min C = Σ_t [ c_ele(t) × P_buy(t) − c_sell(t) × P_sell(t) + c_gas × V_gt(t) + c_gas × V_boiler(t) ]

其中t=1..24，c_ele(t)是分时电价，c_sell(t)是上网电价，c_gas是天然气价格，V_gt(t)是燃气轮机耗气量，V_boiler(t)是燃气锅炉耗气量。如果没有余电上网条件，就强制P_sell(t)=0。

这里有个工程细节：如果系统里有光伏或风电，购电项要改成P_buy(t)，售电项要加上光伏盈余上网。程序在数据文件里预留了P_pv(t)和P_wind(t)两个序列，没有新能源时设成全0，模型自动退化为传统联供系统。

目标函数里还要注意天然气价格单位换算。如果天然气按元/立方米计，燃气轮机耗气量就要按立方米/小时算，公式V_gt = P_gt / (η_gt × LHV)，LHV取9.7 kW·h/Nm³。很多初学者直接拿“每度电耗燃气多少立方米”的经验系数，而不考虑设备效率随负荷率变化，会造成高峰时段燃气轮机出力安排失真。

3.2 能量平衡与出力边界：最容易出bug的两类约束

程序把约束分成四类，每一类作用不同：

1. 电平衡约束：P_gt(t) + P_buy(t) + P_pv(t) + P_dis(t) = P_load_e(t) + P_ec(t) + P_ch(t) + P_sell(t)。意思是发电、购电、新能源、储电放电之和，等于电负荷、电制冷耗电、储电充电、售电之和。电制冷功率P_ec(t)在等式右侧，因为它消耗电能产生冷量。

2. 热平衡约束：Q_hr(t) + Q_boiler(t) + Q_dis_h(t) = Q_load_h(t) + Q_h_ch(t) + Q_waste(t)。这里Q_hr是燃气轮机余热，Q_boiler是燃气锅炉产热，Q_dis_h是储热放热，右侧是热负荷、储热充电、排热。

3. 冷平衡约束：Q_ac(t) + Q_ec_cool(t) = Q_load_c(t)，其中Q_ac(t) = COP_ac × Q_hr_ac(t)，Q_ec_cool(t) = COP_ec × P_ec(t)。还要加一个分配约束Q_hr(t) = Q_hr_heat(t) + Q_hr_ac(t)，表示余热要么供热、要么给吸收式制冷、要么排热。

4. 设备出力边界：每台设备有上下限约束，发电机组配备最小出力限制（如1000 kW机组最小300 kW），制冷机、锅炉同理。

实际运行时最容易出的bug是“变量维度不一致”。YALMIP里定义sdpvar时如果写成sdpvar(24,1)，后面所有约束向量都必须是24×1；如果某个负荷数据是从Excel读进来是行向量1×24，直接相加会报维度错误。我习惯在数据读入后统一加一句P_load = P_load(:)，把所有向量转成列向量，这个习惯帮我省了大量查错时间。

3.3 启停和爬坡：什么时候需要0-1变量

如果只做典型日经济优化、且机组最小出力约束存在，可以把燃气轮机视为全天在线，不引入0-1变量，这样模型是一个纯线性规划LP，求解速度极快。但很多工程场景要求机组可以启停，例如夜间低谷电价极低、燃气轮机发电成本明显高于购电成本，最优策略是夜间停机。这时需要引入启停状态u_gt(t)：

• P_gt(t) ≤ P_gt_max × u_gt(t)
• P_gt(t) ≥ P_gt_min × u_gt(t)
• P_gt(t) − P_gt(t−1) ≤ R_up × Δt + M×(1−u_gt(t))
• P_gt(t−1) − P_gt(t) ≤ R_dn × Δt + M×(1−u_gt(t−1))

第一个和第二个约束让u_gt(t)=0时P_gt必须为0，u_gt(t)=1时出力在上下限之间。爬坡约束用大M法松弛，M取P_gt_max就行，不要取成1e6，否则数值求解器很容易出现病态。

加上0-1变量后，模型从LP变成MILP。24小时只有24个机组启停变量时，Gurobi求解毫秒级；但如果再加入储能充放状态、吸收式制冷启停、售电状态等，变量会膨胀。我的经验是：能不加0-1就不加，先跑纯LP看结果趋势，再决定要不要细化。

4. MATLAB程序落地：从数据表到Gurobi求解

4.1 输入数据怎么组织：负荷曲线、电价、气价

程序目录分三块：data放输入数据，src放主程序和建模函数，result放输出结果和图片。数据文件用一个load_case_data.m脚本组织，全部以结构体形式返回，这样保证主程序代码不需要因为换算例而改动。

%% load_case_data.m caseData.T = 24; caseData.dt = 1; % 时间步长, h caseData.P_load_e = [350; 320; ...]; % 24x1 电负荷, kW caseData.Q_load_h = [480; 450; ...]; % 24x1 热负荷, kW caseData.Q_load_c = [200; 180; ...]; % 24x1 冷负荷, kW caseData.P_pv = [0; 0; ...]; % 24x1 光伏预测, kW caseData.price_buy = [0.35; 0.35; ...]; % 24x1 购电价, 元/kWh caseData.price_sell = 0.40; % 上网电价, 元/kWh caseData.gas_price = 2.8; % 天然气价格, 元/Nm3 caseData.LHV = 9.7; % 天然气低位热值, kWh/Nm3

设备参数我建议单独放一个device_params.m函数，用结构体逐台填写。比如燃气轮机额定功率1000 kW、发电效率0.30、余热热电比1.2、最小出力系数0.3；燃气锅炉额定功率800 kW、效率0.90；电制冷COP 4.5、额定功率500 kW；吸收式COP 1.3、最大输入热功率1200 kW；储电容量1000 kWh、最大充放功率300 kW、效率0.92；储热容量2000 kWh、最大充放功率400 kW、效率0.90。

4.2 用YALMIP搭模型的核心代码结构

主程序用YALMIP建模，求解器调Gurobi。之所以选YALMIP而不是直接调Gurobi的MATLAB接口，是因为YALMIP的约束写法直观，而且可以一行命令切换求解器，方便前期用linprog验证、后期切换到商业求解器。

%% main_cchp_optimization.m addpath(genpath('src')); caseData = load_case_data(); device = device_params(); T = caseData.T; % 决策变量 P_gt = sdpvar(T,1); % 燃气轮机发电功率 P_buy = sdpvar(T,1); % 购电功率 P_sell = sdpvar(T,1); % 售电功率 P_ec = sdpvar(T,1); % 电制冷耗电功率 Q_boiler = sdpvar(T,1); % 燃气锅炉产热 Q_hr_heat = sdpvar(T,1); % 余热供热量 Q_hr_ac = sdpvar(T,1); % 余热给吸收式制冷量 Q_waste = sdpvar(T,1); % 排热量 P_ch = sdpvar(T,1); % 储电充电 P_dis = sdpvar(T,1); % 储电放电 Q_h_ch = sdpvar(T,1); % 储热蓄热 Q_h_dis = sdpvar(T,1); % 储热放热 SOC_e = sdpvar(T,1); % 储电状态 SOC_h = sdpvar(T,1); % 储热状态 Constraints = []; % 燃气轮机出力边界 Constraints = [Constraints, device.P_gt_min <= P_gt <= device.P_gt_max]; % 燃气轮机余热关系 Q_hr = device.beta * P_gt; Constraints = [Constraints, Q_hr_heat + Q_hr_ac + Q_waste == Q_hr]; % 电平衡 Constraints = [Constraints, ... P_gt + P_buy + caseData.P_pv + P_dis == ... caseData.P_load_e + P_ec + P_ch + P_sell]; % 热平衡 Constraints = [Constraints, ... Q_hr_heat + Q_boiler + Q_h_dis == ... caseData.Q_load_h + Q_h_ch]; % 冷平衡 Constraints = [Constraints, ... device.COP_ac * Q_hr_ac + device.COP_ec * P_ec == caseData.Q_load_c]; % 储能约束 Constraints = [Constraints, ... SOC_e(2:T) == SOC_e(1:T-1) + device.eta_e_ch*P_ch(2:T)*caseData.dt ... - P_dis(2:T)*caseData.dt/device.eta_e_dis]; Constraints = [Constraints, SOC_e(1) == 0.5*device.E_e_cap]; Constraints = [Constraints, SOC_e(T) == 0.5*device.E_e_cap]; % 储热同理 Objective = sum(caseData.price_buy.*P_buy) - sum(caseData.price_sell*P_sell) ... + sum(caseData.gas_price * (P_gt/(device.eta_gt*caseData.LHV))) ... + sum(caseData.gas_price * (Q_boiler/(device.eta_boiler*caseData.LHV))); ops = sdpsettings('solver','gurobi','verbose',2); result = optimize(Constraints, Objective, ops);

注意Q_hr = device.beta * P_gt这一行不是约束，而是直接定义了一个中间表达式，YALMIP里可以这样用。如果改成写成约束Q_hr == device.beta * P_gt，同时又单独定义Q_hr变量，也可以，但没必要多一个辅助变量。

4.3 求解器选择和调试参数

我测试过三套求解方案：纯MATLAB的linprog和intlinprog、YALMIP+Gurobi、YALMIP+Cplex。对这套不含0-1变量的LP模型，linprog完全够用，适合学生和白嫖党；但一旦加上启停0-1变量，intlinprog在小规模24时段模型还能跑，扩展到96时段多机组时就很吃力，强烈建议用Gurobi或Cplex。我电脑上实测24时段、4台设备、2个储能的MILP，Gurobi在0.3秒内收敛到gap小于0.01%，intlinprog通常要5~20秒。差距主要在大规模整数规划的分支定界效率上。

调试时遇到Infeasible problem，不要急着怀疑求解器。先做三件事：第一，检查所有等式约束的维度是否一致，可以用size(Constraints) == [1 1]?来判断；第二，把储能SOC约束改回SOC_e(T) >= 0而不是等式，排除末端约束过紧导致的无解；第三，单独注释掉冷平衡约束，看模型能否求解，如果能求，就说明冷侧参数有冲突，比如吸收式COP乘余热最大值还是低于某时刻冷负荷。

5. 算例结果复盘与调优经验

5.1 一套典型日参数下的优化结果

我随手构造了一个夏季典型日算例：燃气轮机额定1000 kW，发电效率0.30，余热热电比1.2；电制冷COP 4.5；吸收式COP 1.3；燃气锅炉800 kW；分时电价峰1.1、平0.68、谷0.35元/kWh；气价2.8元/Nm³。冷负荷峰值约1800 kW，出现在15:00；电负荷峰值约1200 kW；热负荷较小，主要是生活热水，峰值约400 kW。

优化结果符合工程直觉：夜间23:00到次日6:00低谷电价时段，燃气轮机尽量压低到最小出力300 kW，不足电量主要从电网买，储电几乎不充电；上午电价进入平段后，燃气轮机出力抬升到600~700 kW；下午高峰时段13:00~17:00燃气轮机满发1000 kW，余热优先供应吸收式制冷，吸收式制冷机几乎满负荷，电制冷补足剩余冷负荷，购电只买到刚够电负荷的水平。全天总运行成本大约1.86万元。同一个算例如果改成“燃气轮机全天满发”的傻瓜策略，总成本约2.15万元，优化带来了约13.5%的成本下降。

储能结果很有意思：储热在12:00~15:00余热较充裕时蓄热，在18:00~20:00热负荷略有上升且电价较高时段放热；但储电全天充放电量很小，因为峰谷价差0.75元/度，按储电效率0.92折算后单次循环套利空间约0.69元/度，再扣除电池充放损耗后，套利利润勉强覆盖成本，模型自然趋向于“少动”。如果储能单位投资成本或运维成本计入目标，它甚至可能完全不出力。这个结论说明一件事：储能不是CCHP系统的“标配”，经济性要算过才知道。

5.2 储能在经济模型里“失灵”的三个原因

第一是峰谷价差不够大。很多地方分时电价峰谷差在0.6~0.8元/kWh，看似可观，但储电充放一次要损失约8%电量，实际套利空间被压缩。第二是负荷曲线和光伏曲线导致“无谷可充”或“无峰可放”，比如夜间电负荷本身很低，储电充满后白天放电又遇上光伏大发，放电价值不高。第三是模型缺少储能参与需量电费或容量电费的收益项。实际项目里储能收益大头往往来自降低变压器最大需量、减少容量电价，而不是峰谷套利。典型日成本模型如果只算电费，储能很容易被“优化”成闲置资产。所以在结果里看到储能量极少，别急着调大效率参数，先想一想收益项是不是漏了。

5.3 无解、收敛差、结果震荡的排查方法

我在迭代程序过程中遇到过几次典型的求解异常，也总结出了对应的排查顺序。

结果震荡（相邻两个时段的出力大幅跳变）往往来自爬坡约束缺失或大M取值不当。如果模型没有加爬坡约束，优化问题在每个时段独立约束下求解，结果可能出现300 kW和1000 kW来回跳。加爬坡约束时，大M不要取太大，我通常取设备出力的最大值，也就是1000，这样既保证了启停时约束松弛，又不会引入数值误差。另一个震荡原因是负荷数据本身有高频尖峰，给模型传导成了机组出力尖峰。这时先平滑负荷曲线，不要急着改模型。

无解问题，前面提过先用变量维度检查，再检查储能SOC末值约束。实际项目里最阴间的无解原因是“冷平衡约束里COP乘错位置”：吸收式制冷消耗的是热功率，输出冷量是COP_ac × Q_hr_ac；电制冷消耗的是电功率，输出冷量是COP_ec × P_ec。如果误写成Q_hr_ac == COP_ac × 冷负荷，冷负荷大时直接把余热需求推到超过燃气轮机最大余热，必然无解。要养成一个习惯：先把约束的物理单位写清楚，再翻译成代码，不要直接从论文公式抄。

收敛差的最后一个常见原因是目标函数里的成本量级差异太大。购电费一天几万元，储能运维成本一天几毛钱，两者相差1e4以上，求解器可能忽略小项或产生数值噪声。建议把成本统一折算成“元”或“万元”，不要让同一目标函数里出现1e4和1e-2同时存在的情况。

6. 从“能跑”到“可信”：后续还能改什么

6.1 结果合理性校验：手算+极端工况测试

程序跑通只是第一步，结果可信才是关键。我每次改完模型都会做三组校验。第一组是“退化测试”：把燃气轮机效率设成一个很低的值、购电价设成0，模型应该完全不发电只购电；再把设备全设为0，模型应该无解并报出明确的冲突约束。第二组是“人工核验”：挑出负荷最平稳的凌晨时段，手动算一遍电平衡、热平衡和成本，看和优化结果是否吻合。第三组是“能量守恒检查”：单独把储电、储热初末SOC设为相同值后，统计全天的充电量与放电量、充放电损耗，检查SOC递推式残差是否在1e-6以下。

这一套校验做完，程序才谈得上“可用于分析”。否则你拿一个变量名写错、单位写错的模型去跟别人讲“我们优化后成本降低了15%”，那是自己骗自己。我在实际项目里见过不止一次，所谓的最优方案比现场人工调度还差，追查下来都是效率参数填错或者余热关系耦合错误。

6.2 程序的可扩展方向

这套框架后续可以继续加：

1. 目标函数增加碳交易成本，变成经济+低碳双目标，或者通过权重系数把碳排放折算进总费用。
2. 增加需求响应机制，把可平移负荷、可中断负荷建模成弹性负荷变量。
3. 把日前调度改成两层滚动优化：外层做日前计划，内层每15分钟滚动修正，应对光伏和负荷预测误差。
4. 设备效率曲线从常数改成部分负荷率的分段线性函数，更贴近真实机组。分段线性化用YALMIP的binvar和implies命令可以实现，但在模型里引入大量整数变量，求解时间会明显增加。
5. 加入电转热、电转冷等更多能量转换环节，变成真正的多能源系统优化。

我个人在实际操作中的体会是，CCHP经济优化程序的功力不体现在能不能调通一个求解器，而体现在你能不能把每个约束的物理含义讲清楚。刚接触这套程序时，我经常往模型里堆约束，觉得约束越多越精确，后来发现大量冗余约束只会拖慢求解速度，并不会改变最优解。现在我的习惯是：每个约束加进模型前都先问自己一句，如果去掉它，最优解会变成什么样？如果答案是不确定，再把它加回来。这套“先做减法、再做加法”的思路，比任何求解器调参技巧都管用。