双碳目标落到工程层面,核心就一件事:在保证冷热电可靠供应的前提下,把碳排放从“事后核算”变成“事中约束”,让调度优化程序自己算清楚未来24小时每台设备的启停和出力。我这一年多一直在用Matlab做综合能源系统低碳运行优化调度的仿真验证,从模型搭建到求解落地完整跑了一遍,踩了不少坑,也攒了很多可以直接“抄作业”的写法。这篇文章就把整套思路拆开讲清楚:模型怎么建、碳成本怎么进目标函数、约束怎么写,以及用Matlab+YALMIP实现混合整数线性规划调度模型的完整步骤和避坑实录。适合正在做综合能源调度方向毕业论文、或者刚到设计院/研究院接触碳交易和能量管理的朋友,有一定电力系统或运筹优化基础但还没动手写过调度代码的最好。
1. 从“为什么要做”说起:双碳约束下的调度难题
1.1 综合能源系统低碳调度的现实需求
综合能源系统跟传统电力系统最大的区别在于:电、气、热、冷四种能量流耦合在一起,同一种负荷需求往往有多条能量转换路径可选,而每条路径的碳排放强度可能差出好几倍。
举个例子,同样是满足冬季园区热负荷,你可以开燃气锅炉直接烧气产热,也可以让燃气轮机多发电、再把烟气余热回收利用;同样是满足电负荷,你可以从电网买电,也可以让燃气轮机自发电。如果只看经济成本,分时电价低的时候买电可能很划算;如果叠加上碳成本,这种决策就可能完全反过来——因为电网购电的间接碳排放因子通常比天然气机组还要高。这种情况在“双碳”目标提出来之前基本没人认真考虑,因为碳排放没有价格,自然就不会进入优化目标。但当碳配额和碳交易机制逐步落地,碳排放变成一项有真金白银成本的约束之后,调度优化就必须从“最小化运行成本”升级为“在碳排放代价与运行成本之间找平衡点”。
我把这个问题拆成三层来看:
- 设备层:燃气轮机、电锅炉、吸收式制冷机、储电、储热,每一台设备的出力区间、爬坡能力、效率特性都不一样。
- 系统层:电、热、冷三条母线的功率平衡需要同时满足,储能跨时段耦合带来记忆性约束。
- 政策层:碳排放配额是多少、碳价怎么定、超排怎么罚,这些参数直接改变目标函数的形状。
1.2 为什么选Matlab做这件事
很多人问我:现在Python这么火,深度学习也能挂在优化外面做预测,你干嘛还用Matlab?我的回答一直是:Matlab在这个特定场景下不是因为它“新”,而是因为它“稳”。
综合能源系统调度在数学上是一个典型的多时段混合整数线性规划问题,里面有0-1启停变量、连续功率变量、跨时间的储能递推约束。YALMIP这个建模工具箱跟Matlab的结合非常成熟,写变量约束就像写数学公式一样直接。虽然Python有Pyomo,但我个人实测下来,YALMIP在语法简洁度、报错提示的清晰度、以及跟求解器的衔接上,对做调度优化的人来说还是最顺手的一个。另外Matlab内置的Optimization Toolbox、Sensor数据清洗、绘图工具配合起来也很省事,尤其是做敏感性分析的时候,一套代码循环跑几十个碳价场景,结果自动出图,比反复切换Python的matplotlib和pandas要直观得多。
1.3 建模语言选用背后的代价考量
用YALMIP并不意味着所有问题都迎刃而解,它只是一个建模层的语法糖。真正决定模型能不能快速求解的还是数学形式。我在实际项目里坚持三个原则:
- 能用线性约束绝不用非线性约束。哪怕需要牺牲一点点精度,只要结果对工程判断没有实质影响,就坚持线性化。
- 整数变量的数量要严格控制。每加一个0-1变量,求解时间都可能是非线性增长的事。
- 参数全部集中放在一个data结构中管理,不要散落到脚本各处。后期换数据、跑敏感性分析,直接改结构体字段就行。
2. 模型怎么搭:冷热电联供系统的数学框架
2.1 设备模型与能量流拓扑
一个典型的综合能源园区,能量流拓扑是这样的:电网入口、天然气入口、光伏和风电作为不可控电源进入电力母线;燃气轮机发电上电力母线,余热进入余热回收装置以补热力母线;电锅炉直接电转热;吸收式制冷机从热力母线取热产生冷量进入冷力母线;蓄电和蓄热分别跨时段转移电量和热量。冷负荷可以来自吸收式制冷机,也可以来自压缩式电制冷机(如果园区有)。
数学上,三母线功率平衡可以写成:
电力平衡: P_gt(t) + P_pv(t) + P_wt(t) + P_es_d(t) + P_grid(t) = P_load(t) + P_eb(t) + P_ec(t) + P_es_c(t)
热力平衡: H_gt(t) + H_eb(t) + H_hs_d(t) = H_load(t) + H_ar(t) + H_hs_c(t)
冷力平衡: C_ar(t) + C_ec(t) = C_load(t)
其中P_es_c和P_es_d是蓄电的充/放功率,同一时刻只能有一个方向打开;蓄热类似。这组平衡约束是每天24小时都成立的,每个时段都要写一组等式进去。刚开始建模的人最容易漏掉的是电制冷机这条支路——如果你园区里同时有吸收式制冷和电制冷,冷母线两条来源都要写进去,否则求出来的结果可能在物理上不可行。
2.2 目标函数里的“碳账本”
低碳调度跟传统经济调度的根本区别,就在目标函数里多了一项碳排放相关的成本项。我常用的是“碳交易成本”框架,公式写出来是:
min F = Σ_t [ C_fuel(t) + C_grid(t) + C_start(t) + C_stop(t) ] + C_carbon_total
其中C_fuel是天然气燃料成本(与燃气轮机、燃气锅炉耗气量有关),C_grid是购电费用(分时电价×网购电量),C_start和C_stop是机组启停成本,最后一项C_carbon_total是整个调度周期的碳交易结算金额。
碳交易成本又拆成两块:
实际碳排放量。电网购电的间接排放 = Σ 购电量 × 电网排放因子;天然气燃烧的直接排放 = Σ 耗气量 × 天然气排放因子。这里的排放因子我建议用当地碳市场公布的最新值,例如某省电网排放因子约0.55 kgCO2/kWh,天然气机组约0.20 kgCO2/kWh,不同区域差异很大,直接抄网上的数很容易算错。
碳配额差额。系统会拿到一个免费的碳排放配额E_limit,如果实际排放E_actual高于配额,就要按碳价λ购买差额;如果低于配额,多余部分可以卖出获得收益。所以:
C_carbon_total = λ × (E_actual - E_limit)
λ也就是碳市场的交易价格,是模型里最敏感的一个参数。后面会在案例里专门做敏感性分析。
2.3 约束条件的物理意义
约束条件我分成四类,每一类都有对应的物理含义,不能乱写,否则模型会“数学上可解、物理上胡说”。第一类是设备出力上下限,比如燃气轮机出力在30%-100%额定功率之间,风电光伏的出力上限就是该时段的预测值。第二类是爬坡约束,燃气轮机和电锅炉的出力变化速率必须落在规定范围内,否则就是纯理论上的“瞬间拉满”,工程中不可能实现。第三类是储能约束,包括充放功率上限、容量上限、充放效率,以及状态量SOC的递推关系,每天最后一时刻SOC还要等于初始值,否则就相当于从系统外白拿了一块能量。第四类是我个人特别强调的碳排放硬约束,在碳交易机制之外还可以叠加一个调度周期碳排放总量上限,这样即使碳价波动,系统也不会突破某个绝对排放值。
这些约束写出来之后,整个模型就是一个标准的MILP。为什么非要是MILP?因为储能充放方向必须互斥,燃气轮机启停是0-1状态,这些天然就是整数变量;如果不处理成MILP而做成连续线性规划,求出来的解虽然快,但不符合实际操控逻辑。
3. Matlab实现的关键环节
3.1 数据准备与基础参数设置
建模之前数据准备这步最容易被忽视,但往往就是结果离谱的根源。我的数据清单包括:
- 24小时电、热、冷负荷曲线(室内温度、生产班次等共同决定,通常从DEST或EnergyPlus仿真导出)。
- 光伏、风电24小时出力预测曲线(典型日可由历史数据聚类得到)。
- 分时电价曲线(峰谷平时段单价差异明显)。
- 天然气价格、电网排放因子、天然气排放因子。
- 碳配额总量和碳交易价格。
- 设备参数表:额定容量、出力上下限、爬坡速率、效率、启停成本。设备参数齐了,就可以把这些数据统一整理成一个结构体,类似这样:
par.T = 24; par.grid.pmax = 30; % 网购电上限 MW par.grid.price = [0.42*ones(1,8), 0.82*ones(1,8), 0.57*ones(1,8)]; % 分时电价 par.gt.pmax = 40; par.gt.pmin = 12; par.gt.eta_e = 0.35; % 发电效率 par.gt.eta_h = 0.45; % 余热回收效率 par.gt.ramp = 8; % 爬坡 MW/h par.carbon.lambda = 100; % 碳价 元/t par.carbon.quota = 400; % 配额 t par.eb.pmax = 15; par.eb.eta = 0.95; par.es.pmax = 5; % 储电最大充放功率 MW par.es.capacity = 20; % 储电容量 MWh par.es.init_soc = 10; % 初始SOC MWh为什么要用结构体而不是一堆散变量?因为后期跑不同场景、不同参数组合时,只需修改这个结构体的字段,脚本主体不用动。如果把参数散写在脚本里,改一次参数就得全局翻一遍,非常容易改漏。
3.2 用YALMIP定义决策变量与约束
yalmip跟Matlab的矩阵思维天然契合,可以用sdpvar定义一个行向量来表示24小时的决策量,所有约束一次性向量化写入,避免写24个for循环。基本骨架我一般这样写:
T = par.T; % 连续变量:各设备24小时出力 Pgt = sdpvar(1, T, 'full'); Peb = sdpvar(1, T, 'full'); Pgrid = sdpvar(1, T, 'full'); Ppv = sdpvar(1, T, 'full'); Pwt = sdpvar(1, T, 'full'); Soc_es = sdpvar(1, T, 'full'); Pes_c = sdpvar(1, T, 'full'); % 充电功率 Pes_d = sdpvar(1, T, 'full'); % 放电功率 Har = sdpvar(1, T, 'full'); % 吸收式制冷机耗热量 % 0-1变量:燃气轮机启停、储能充放方向 u_gt = binvar(1, T); u_es_c = binvar(1, T); u_es_d = binvar(1, T);变量定义完之后开始写约束。这里有个经验:约束先写物理边界(上下限、爬坡),再写耦合约束(母线平衡),最后写储能跨时段耦合约束。注意储能充放电的互斥约束,一定要用同一个时间段的binvar把它们关联起来:
Constraints = []; % 燃气轮机约束 Constraints = [Constraints, par.gt.pmin*u_gt <= Pgt <= par.gt.pmax*u_gt]; Constraints = [Constraints, -par.gt.ramp <= Pgt(2:T) - Pgt(1:T-1) <= par.gt.ramp]; % 电平衡 Constraints = [Constraints, Pgt + Ppv + Pwt + Pes_d + Pgrid == ... Pload + Peb + Pec + Pes_c]; % 热平衡 Constraints = [Constraints, par.gt.eta_h/par.gt.eta_e*Pgt + Peb*par.eb.eta ... + Phs_d == Hload + Har + Phs_c]; % 储能约束 Constraints = [Constraints, Soc_es(1) == par.es.init_soc]; Constraints = [Constraints, Soc_es(2:T) == Soc_es(1:T-1) + ... (Pes_c(1:T-1)*par.es.eta_c - Pes_d(1:T-1)/par.es.eta_d)]; Constraints = [Constraints, 0 <= Soc_es <= par.es.capacity]; Constraints = [Constraints, 0 <= Pes_c <= par.es.pmax*u_es_c]; Constraints = [Constraints, 0 <= Pes_d <= par.es.pmax*u_es_d]; Constraints = [Constraints, u_es_c + u_es_d <= 1]; Constraints = [Constraints, Soc_es(T) == par.es.init_soc];这里有几个坑必须强调。第一个坑是SOC递推公式里的效率位置——充电效率乘在充电功率上,放电效率除在放电功率上,位置写反会导致储能系统“凭空造能量”,第二天开局SOC对不上。第二个坑是电平衡等式右边别忘了电制冷机Pec和电锅炉Peb,很多初学者把Peb当成热负荷的一部分写在热平衡里,结果电力平衡写少了,求解器为了凑等式会凭空增加购电。第三个坑是0-1变量的互斥约束u_es_c + u_es_d <= 1,不加这个约束的话,同一个储能设备可能被优化成同一时刻既充电又放电,两边功率都在跑,系统还能倒赚效率差的能量。
3.3 求解器配置与目标函数实现
目标函数分成购电成本、燃料成本、启停成本和碳交易成本四部分。燃料成本的天然气的量,通过燃气轮机的发电量除以发电效率换算,表达式写起来要仔细:
% 购电成本 C_grid = sum(par.grid.price .* Pgrid); % 燃料成本(燃气轮机耗气量×天然气单价) C_fuel = sum(Pgt / par.gt.eta_e * par.gas_price / 3600 * 10^3); % 启停成本 C_switch = sum(par.gt.start_cost * max(u_gt(2:T) - u_gt(1:T-1), 0) + ... par.gt.stop_cost * max(u_gt(1:T-1) - u_gt(2:T), 0)); % 碳排放量 E_grid = sum(Pgrid) * par.carbon.ef_grid; % 购电间接碳排放 E_gas = sum(Pgt / par.gt.eta_e) * par.carbon.ef_gas; % 天然气燃烧排放 E_actual = E_grid + E_gas; % 碳交易成本 C_carbon = par.carbon.lambda * (E_actual - par.carbon.quota); % 总目标 Objective = C_grid + C_fuel + C_switch + C_carbon;启停成本那里用了max函数,这在YALMIP里会引入额外的二进制变量,实现上稍微慢一点。如果场景中启停次数很少,可以简化成只计启动成本或者忽略。目标函数写完,调用求解器:
ops = sdpsettings('solver', 'gurobi', 'verbose', 1); ops.gurobi.MIPGap = 0.01; ops.gurobi.TimeLimit = 120; sol = optimize(Constraints, Objective, ops); if sol.problem == 0 Pgt_opt = value(Pgt); ... else disp('求解失败,请检查约束有无冲突'); check(Constraints); endcheck(Constraints)这一步非常重要,它会逐一报告每条约束的残差,快速定位出错的是功率平衡还是储能约束。我强烈建议新手养成求解后立即跑check的习惯,而不是直接去看结果。
3.4 结果可视化与敏感性分析
求解完成之后,通常需要画两组图。第一组是各设备24小时出力堆叠图,用bar函数画电功率和热功率的堆叠柱状图,一眼能看出哪个时段用了什么能量来源。第二组是储能SOC曲线和碳排放累计曲线,看储能是否起到了削峰填谷、以及碳排放是否真的降下来了。
我习惯把结果同时导出成一个结构体,方便后续对比不同碳价下的调度方案:
result.Pgt = value(Pgt); result.Pgrid = value(Pgrid); result.Pes_c = value(Pes_c); result.Pes_d = value(Pes_d); result.Carbon = value(E_actual); result.Cost = value(Objective);之后做敏感性分析,就是把碳价lambda作为外层循环变量,里面重复调用optimize。这个流程我一般写在独立脚本里,跑一次大概几分钟,最后把所有结果横向对比。
4. 我在实际调参中踩过的坑
4.1 非线性项处理:设备效率怎么线性化
燃气轮机的发电效率并不是固定值,低负载的时候效率明显下降。如果对精度要求高,最简单的线性化手段是把出力区间分成三段,每段用不同的线性效率表达式。比如40MW机组,12-20MW为一档,20-32MW为一档,32-40MW为一档,每档对应一个常数效率。YALMIP里可以用implies语句或者Big-M方法来实现分段函数的建模,但这需要引入额外的0-1变量。
我在工程中往往直接采用恒定效率,尤其是初步可行性研究阶段。因为燃气轮机的变工况效率数据很多时候厂家都不一定给得全,而且分段线性化带来的求解时间增长和数值稳定性问题可能得不偿失。如果你做的是学术研究,那分段线性化是加分项,但必须在论文里说明分段依据和误差来源。
4.2 碳配额与碳价设置的坑
碳配额E_limit不是瞎拍的。我见过很多人直接把配额设成一个特别小的数,结果模型直接无解——因为所有设备全开满都达不到排放要求。正确做法是先跑一次不带碳排放项的纯经济调度,记录下此时的碳排放量E_base,然后把配额设成E_base的某个比例,比如0.8×E_base,这样模型一定有可行解,同时又构成了实质性的碳约束。
碳价参数的设置同样要谨慎。碳价太低,模型会完全忽略碳约束,结果跟经济调度几乎一样;碳价太高,模型可能会把电网购电全砍掉,让燃气轮机超负荷运行,这种结果经济上不可行、工程上也不合理。所以我每次做案例都先扫一遍λ的敏感性曲线,大概了解“拐点”在哪——碳价超过某个值之后碳排放量下降开始停滞,说明系统已经到了减排极限。
4.3 求解器选型与收敛性判断
同一套YALMIP模型,换不同的求解器,结果和速度可能差异很大。我的经验是:中小规模算例(24时段、几十个节点)用Gurobi或CPLEX都很快,几秒到几十秒能收敛到1%的MIPGap。如果只有Windows机器没有商业求解器授权,用GLPK或者SCIP也凑合能跑,但速度会慢不少,大一点的问题可能跑到几百秒。
还有一类问题值得警惕:YALMIP默认的求解器如果检测不到,它会自动用一个内置的纯线性求解器,那个对MILP基本无能为力,可能跑很久不收敛,也可能给出奇怪结果。所以调用前务必检查sdpsettings里面solver字段是否指定正确。另外我建议在求解设置里加上TimeLimit,避免出现优化几小时出不来结果的情况。工程调度是滚动运行的,单次求解超过10分钟基本就失去实际意义了。
5. 案例复现:一个典型园区的低碳调度结果
5.1 算例场景与基础数据
为了直观展示低碳调度跟传统经济调度的差异,我构造了一个简化的园区算例。园区装机:燃气轮机40MW、光伏30MW、风电20MW、电锅炉15MW、吸收式制冷机15MW、蓄电10MW/20MWh、蓄热10MW/20MWh。冬季典型日,电负荷峰值约68MW,热负荷峰值约25MW,冷负荷很小这里先忽略。分时电价峰段1.0元/kWh、平段0.6元/kWh、谷段0.3元/kWh,天然气价格2.5元/m3,电网排放因子0.55 kgCO2/kWh,天然气机组排放因子按0.20 kgCO2/kWh煤气计算,碳配额取纯经济调度排放量的75%,碳价先设为100元/吨。
5.2 两种调度结果详细对比
跑完模型后我把纯经济调度(目标函数去掉碳交易项)和低碳调度(完整目标函数)的结果放在同一张表里对比,这里给出一组典型日结果:
| 指标 | 纯经济调度 | 低碳调度 | 变化 |
|---|---|---|---|
| 全天购电量(MWh) | 320 | 215 | -32.8% |
| 燃气轮机发电量(MWh) | 520 | 640 | +23.1% |
| 光伏消纳量(MWh) | 288 | 288 | 0% |
| 风电消纳量(MWh) | 195 | 195 | 0% |
| 碳排放总量(t) | 452 | 378 | -16.4% |
| 总运行成本(万元) | 38.6 | 41.2 | +6.7% |
| 碳交易结算成本(万元) | 0 | 3.35 | — |
这个表很有信息量。购电量降下来之后,碳排放确实少了,但总成本却上升了,因为天然气比高峰电价时段买电要贵。碳价100元/吨的时候,系统每年“少排”的碳折成钱还不足以覆盖燃料成本增加,所以这其实是一个权衡空间很大的问题。
5.3 结果解读与调度策略分析
再看设备出力曲线,低碳调度下燃气轮机的运行时段明显拉长,以前在后夜低谷时段可能关掉靠买电,现在保持低负载运行,为的就是减少高峰时段电网购电的间接排放。蓄电SOC曲线也有意思:纯经济调度是典型的“谷充峰放”,低碳调度则倾向于在光伏出力高的中午时段充电——因为那时候燃气轮机已经满载,如果电力母线还多出光伏电,与其卖给电网不如存起来,相当于间接降低了系统整体碳排放。
用数据说话,这个模型输出的结论对实际决策很有价值:如果未来碳价从100元/吨涨到200元/吨,购电量还会继续下降,但下降斜率明显趋缓。这说明园区在现有技术条件下“减排潜力”是有边界的,再要继续减排就得靠上更大容量的储能或者增加高效设备,光靠调度算法优化已经无法进一步挖掘空间了。
6. 常见问题速查表与心得
6.1 高频问题排查速查表
| 现象 | 可能原因 | 解决办法 |
|---|---|---|
| 求解器报错Infeasible | 碳排放配额设置过严,或某个功率平衡等式漏写了负荷项 | 先去掉碳项跑一遍,确认基础可行域;再用check(Constraints)逐个检查约束残差 |
| 储能既充电又放电 | 缺少u_es_c + u_es_d <= 1互斥约束 | 补上储能方向互斥条件 |
| SOC到第二天对不上 | SOC递推公式效率位置写反,或首末平衡约束缺失 | 检查SOC递推里充电效率乘在哪一侧,确保Soc_es(T)==Soc_es(0) |
| 结果中光伏、风电弃电却还购电 | 机组爬坡约束太紧,或者最小出力限制了消纳空间 | 放宽爬坡限制,或者在电力平衡里加入弃风弃光变量,让模型来决定是否弃电 |
| Gurobi跑很久不出解 | 整数变量过多或MIPGap设置过严 | 设置MIPGap=0.01或TimeLimit=120,接受次优解 |
| 碳价很高但碳排放没降 | 系统结构调整空间已耗尽 | 检查敏感性曲线拐点,考虑改变装机方案而不是继续调碳价 |
6.2 参数可解释性检查
很多同学拿到优化结果先看成本,再看出力曲线,却不检查结果是否符合物理直觉。我一般会让程序额外输出三个“体检指标”:一是系统弃光弃风率,正常情况下应该很低,如果超过5%,说明储能不能匹配光伏波动;二是燃气轮机启停次数,一天超过四次就要检查是不是碳价导致的频繁切换;三是储能SOC曲线斜率的方向,充电时段应该对应低成本或高光伏时段,如果反过来了,多半是约束写错。
这三个指标不直接出现在目标函数里,但能帮我们判断整套模型的逻辑是否正常,属于“用经验去验证程序”。
6.3 我个人实操下来的几点体会
跑了大半年调度模型,最深的感触是:模型百分之七八十的精力都花在约束表达和参数合理性上,真正调求解器参数的时间反而很少。我一开始也是拿到代码就想立刻优化,结果被无效约束坑了两周。后来学乖了,每加一类新设备,就单独标定它的子系统,确保单设备平衡约束没问题再整合到整体模型里。这套渐进式建模习惯,帮我省下的时间比任何求解器优化都多。
另外,千万别只看单一场景的结果。我在案例里做碳价敏感性扫描时发现,100元/吨和300元/吨的结论差异非常大,甚至可能反过来影响“要不要新建储能”的投资决策。做调度优化,本质上做的是marginal analysis,不是求一个解就结束。建议每个人在做完项目后,都去跑一遍关键参数从低到高的扫描曲线,你对你系统的理解会明显比只看一组结果清晰很多。