☰
综合能源系统优化调度:碳交易与需求响应的Matlab建模实战
2026/10/10 22:26:03 网站建设 项目流程

做综合能源系统调度的时候,很多人会把"碳交易"和"需求响应"当成两个孤立的功能模块:碳交易嘛,就是在目标函数里加个碳价乘排放量;需求响应嘛,就是让负荷曲线削个峰填个谷。但实际上,这两个机制一旦同时进入优化模型,它们会通过设备出力、购电策略、储能充放电互相影响,甚至可能改变整个系统的运行方式。这篇文章就来拆解我是怎么在Matlab里把碳交易和需求响应同时放进综合能源系统优化调度模型,并且真正跑出可解释的结果。

这篇内容适合几类人:正在做综合能源、微电网、园区多能互补相关毕设或课题的同学,需要复现"碳交易+需求响应"优化调度代码的工程师,以及想了解从数学建模到Matlab求解全流程的入门者。我会把碳配额的核算方式、阶梯碳价的线性化处理、价格型与激励型需求响应的建模思路,以及Yalmip+Cplex/Gurobi这套求解组合的关键代码都过一遍。文末还会聊聊我实际调试中踩过的一些坑——这些东西光看论文是学不到的。

1. 为什么要把碳交易和需求响应放进同一个优化模型

1.1 传统经济调度只算"用能成本",不算"碳排放成本"

过去做综合能源系统调度,目标函数通常只有三类成本:购能成本(买电、买气)、设备运行维护成本,以及可能的弃风弃光惩罚。这个框架本身没问题,但在双碳背景下它有一个明显的盲区——碳排放变成了一个"外部性",系统不会主动去约束它。

举个例子,如果天然气价格很低、而电网购电价格相对较高,传统调度模型会让燃气轮机(CHP)满发,因为燃气发电的单位成本更低。但如果把碳交易成本算进来,天然气燃烧的直接排放和电网购电的间接排放都要承担碳成本,这时候CHP的"低价优势"可能就被碳价抵消了。更复杂的是,碳配额的计算方式会影响这个判断——如果系统获得了免费配额,只要实际排放低于配额,甚至还能通过出售富余配额赚钱。这就不是简单加个碳价的问题了。

1.2 碳交易与需求响应如何互相影响调度结果

需求响应进入模型后,系统端的负荷曲线不再是固定的"刚性需求",而是可调节的。这里有个容易忽略的联动关系:

  • 价格型DR通过分时电价引导用户把峰时负荷挪到谷时,系统侧的购电曲线随之变化。谷时负荷增加,意味着原本在谷时被压低的CHP出力可以适当上调,或者储能可以在谷时充入更多电能。
  • 碳交易机制会改变设备的边际成本排序。碳价高的时候,高排放设备(燃气锅炉、没装CCUS的CHP)的出力会被压降,低排放设备(电锅炉、热泵)的出力会上升。如果DR把负荷转移到了这些低排放设备的可用时段,碳排放和运行成本就能同时降低。
  • 储能的作用也会被重新定义。传统调度里储能主要是套利——低充高放;引入碳交易后,储能的充放电策略还会影响购电量和自发电量的配比,进而影响间接排放和碳配额盈亏。

所以两件事必须放到同一个优化模型里联立求解,分开做等价于"先定负荷曲线再定设备出力",而实际上负荷和设备出力是同时被决策出来的。

2. 碳交易机制建模:配额核算、差额结算与阶梯碳价线性化

2.1 免费配额怎么算:基准线法下的排放差额

碳交易建模的第一步是确定配额。国内目前电力行业用基准线法比较多,综合能源系统里常见做法是对供电量和供热量分别给定单位配额系数:

E_quota = δ_e × ΣP_e,t + δ_h × ΣH_t

其中δ_e是单位电量的免费配额(单位tCO2/MWh),δ_h是单位热量的免费配额(单位tCO2/GJ)。这个系数的取值通常参考行业基准,比如供电配额0.45~0.7 tCO2/MWh,供热配额0.2~0.3 tCO2/GJ。如果是做算例设计,也可以自己设定,但要注意合理性——配额太紧会导致碳成本占比过高,太松则碳交易形同虚设。

实际排放量这边,主要包含两个来源:一个是天然气消耗的直接排放,一个是外购电力的间接排放:

E_emission = Σ F_gas,t × EF_gas + Σ P_buy,t × EF_grid

天然气排放因子EF_gas按低位热值折算大约是0.2 kgCO2/kWh(实际项目里用0.185~0.21都正常),电网排放因子EF_grid国内不同区域差异很大,从0.3到0.7 kgCO2/kWh都有,算例里取0.4~0.6比较常见。

然后定义排放差额:

ΔE = E_emission - E_quota

ΔE>0说明配额不足,需要在碳市场购买;ΔE<0说明有富余配额,可以出售。要注意符号方向——很多初次写代码的人在这里搞反,最后目标函数里碳成本变成了"负数越大越好",整个调度逻辑就崩了。

2.2 阶梯碳价的线性化:0-1变量与大M约束

实际碳市场往往采用阶梯价格,目的是让排放量越高的主体边际惩罚越大。常见设定是:

  • 若 ΔE ≤ 0,则碳收益为 p_c × ΔE(此时为负成本,也就是收益),其中p_c是基准碳价;
  • 若 0 < ΔE ≤ λ1,碳成本为 p_c × ΔE;
  • 若 λ1 < ΔE ≤ λ2,碳成本为 p_c × λ1 + 1.1p_c × (ΔE - λ1);
  • 若 ΔE > λ2,则超出部分按1.2p_c计,即 p_c × λ1 + 1.1p_c × (λ2 - λ1) + 1.2p_c × (ΔE - λ2)。

这个函数是分段线性的,不能在Matlab里直接用if-else写进目标函数——因为优化求解器需要的是显式的数学表达式,而不是程序逻辑。分段线性函数的标准做法是引入0-1变量把定义域切分成几个区间,再用大M约束把每个区间的成本和排放差额对应起来。

具体实现时,可以定义三个连续非负变量E1、E2、E3,分别表示落在三个梯度区间内的排放差额部分,以及对应的0-1变量z1、z2、z3:

  • E1代表第一档区间 [0, λ1] 内的排放量,约束 E1 ≤ λ1 × z1;
  • E2代表第二档区间 (λ1, λ2] 内的排放量,约束 λ1 × z2 ≤ E2 ≤ λ2 × z2;
  • E3代表第三档区间 (λ2, ∞) 内的排放量,约束 E3 ≥ 0,且用到时 z3=1。

同时令 ΔE_pos = E1 + E2 + E3,且 z1 + z2 + z3 = 1(只有当ΔE_pos>0时才需要这些变量,ΔE≤0的部分另算收益)。目标函数里碳成本部分就可以写成:

C_carbon = p_c × E1 + 1.1 × p_c × E2 + 1.2 × p_c × E3 - p_c × E_sell

其中E_sell是富余配额出售量,约束 E_sell ≥ -ΔE(ΔE为负时成立)。

这套东西看起来繁琐,但它是把"不可导的分段函数"变成"线性约束+MILP可解形式"的关键。我在Yalmip里通常用binvar定义z变量,再用大M系数(M取远大于物理排放量上限的数,比如1e4)来写互斥条件。

2.3 配额制度下的一个特殊边界:负差额是否允许出售

有些算例会假设富余配额可以自由出售,有些则不允许(设一个最低持有量)。这个设定会显著影响系统行为。如果允许出售,系统可能因为电锅炉供热替代燃气锅炉而获得额外的碳收益;如果不允许出售,排放低于配额的部分就只是"不亏",系统的减排动力会弱一些。

我的建议是算例里把两种情况都跑一遍,对比结果,这样论文或报告里能多一个有意义的敏感性分析。实现方式也很简单:出售变量E_sell的增加一个上限约束(比如不超过0)就是在代码里一行的事。

3. 需求响应建模:价格型弹性矩阵与激励型可中断负荷

3.1 价格型DR:用弹性系数矩阵描述负荷转移

价格型需求响应的经典模型是弹性矩阵。它的思路是:用户对电价的反应可以用自弹性和交叉弹性来描述。自弹性是当前时段价格变化对当前时段负荷的影响(通常为负,即电价越高负荷越低),交叉弹性是其他时段价格变化对本时段负荷的影响(通常为正,即其他时段涨价会把负荷挤到本时段)。

用数学式表达:

L_dr,t = L_base,t + Σ_t' ε(t,t') × (π_t' - π_0,t') / π_0,t' × L_base,t

其中π_t'是优化后(或者分时电价政策)的价格,π_0,t'是基准价格。这个公式里负荷转移量和价差、基准负荷、弹性系数三者直接挂钩,非常直观。

不过在实际算例中,弹性矩阵法有一个麻烦:价格变量π本身往往也是调度模型的决策变量(尤其在引入需求响应后,有些模型把电价设计成实时电价),这时候公式里出现π_t' × L_base,t这种双变量乘积,会让问题变成非线性。

我的处理办法是避开来解决——把需求响应当作已知的分时电价方案,比如峰谷平时段各设几档固定价格,弹性矩阵计算出来的负荷变化率直接在优化前算好,再代入调度模型。这样既保留了DR对负荷曲线的重塑效果,又不引入非线性。如果一定要做电价内生化,那就要用二元乘积线性化的手段,复杂度会上升一个档次。

3.2 激励型DR:可中断负荷的0-1变量处理

激励型DR相对更好建模。它的核心逻辑是:系统在高峰时段可以呼叫用户中断一部分负荷,作为交换,用户获得补偿,像签订一个"可中断负荷合同"一样。

数学上,每个时段定义可中断负荷量ΔL_DR,t,配有0-1状态变量u_DR,t:

0 ≤ ΔL_DR,t ≤ ΔL_DR,max × u_DR,t

∑_t u_DR,t ≤ N_max(全天最多中断N次)

∑ ΔL_DR,t ≤ E_DR,max(全天累计中断电量上限)

目标函数里加上补偿成本:C_DR = Σ c_DR × ΔL_DR,t

关键是要把"中断电量"从该时段的电负荷里扣掉——它等于系统少供的那部分电。这个负向负荷会直接影响功率平衡约束:P_supply,t = L_base,t - ΔL_DR,t + P_storage_net + ...。我在初版代码里就吃过一次亏:DR变量定义好了,目标函数也加了成本,但忘了更新功率平衡方程,结果DR形同虚设。

3.3 两类DR联合使用的协调约束

价格型DR的变化量受弹性矩阵约束,是"连续但总量守恒"的——用户转移负荷总量不减少,只是挪到别的时段。激励型DR则是"真减少",中断的负荷就是少用了。

当两者同时存在时,要注意:价格型DR把峰荷转移去谷时,会抬高谷时负荷;激励型DR在峰时中断部分负荷,会进一步压低谷时需要充入的功率。如果储能策略跟不上,可能出现"谷时负荷依然很高、储能没空间充电、光伏被弃"的尴尬结果。所以调度模型里这两类DR的变量都要参与功率平衡,不能先算好价格型再叠加激励型,否则会重复计算负荷削减量。

4. Matlab建模与求解:Yalmip与求解器选型

4.1 为什么不用Matlab内置优化工具箱

很多人一上来用fmincon或者linprog做调度优化,遇到0-1变量就只能用intlinprog。小规模算例(比如单设备、24时段)intlinprog勉强能跑,但综合能源系统一旦加上CHP热电耦合、储能SOC时序约束、碳交易阶梯区间变量和DR中断变量,整个模型的决策变量数量和约束矩阵规模会迅速膨胀,intlinprog的求解效率会变得很难看,而且数值稳定性一般。

我更推荐Yalmip + Cplex/Gurobi这套组合。Yalmip是一个Matlab下的建模层,它把优化问题用人类友好的方式写出来(sdpvar、binvar、constraints、objective),底层调用商业求解器。Gurobi在MILP上几乎是当前最快的,学术license申请也很方便。Cplex如果拿不到新license,老版本配合旧版Matlab也能用。我自己现在主力是Gurobi,备胎是Cplex。

4.2 核心代码骨架:从变量定义到optimize调用

下面给一个能跑通基本框架的Matlab+Yalmip片段,覆盖了我们前面讨论的关键要素:

%% 基础数据 T = 24; dt = 1; % 时段长度1h load_base = [...]; % 基础电负荷,1x24 heat_load = [...]; % 热负荷,1x24 price_buy = [...]; % 分时购电价 price_DR = [...]; % 价格型DR套用的电价 %% 决策变量 P_chp = sdpvar(1, T); H_chp = sdpvar(1, T); P_gb = sdpvar(1, T); % 燃气锅炉 P_eb = sdpvar(1, T); % 电锅炉 P_pv = sdpvar(1, T); % 光伏实际出力(<=预测值) P_wt = sdpvar(1, T); P_buy = sdpvar(1, T); P_sell = sdpvar(1, T); % 向电网售电 P_es_c = sdpvar(1, T); % 电储能充电 P_es_d = sdpvar(1, T); % 电储能放电 SOC = sdpvar(1, T+1); % 电储能SOC u_es = binvar(1, T); % 储能充放互斥 u_chp = binvar(1, T); % CHP启停 u_dr = binvar(1, T); % 激励型DR中断状态 L_dr_inc = sdpvar(1, T); % 激励型DR中断量 L_dr_price = sdpvar(1, T); % 价格型DR变化量(可为负,表示负荷转移到该时段)

然后是约束的骨架:

Constraints = []; % 电功率平衡 Constraints = [Constraints, P_buy - P_sell + P_pv + P_wt + P_chp ... - P_eb - P_es_c + P_es_d == load_base + L_dr_price - L_dr_inc]; % 热功率平衡 Constraints = [Constraints, H_chp + P_gb + P_eb == heat_load]; % CHP热电耦合(简化定热电比) k_chp = 0.9; Constraints = [Constraints, H_chp == k_chp * P_chp]; Constraints = [Constraints, P_chp >= 0.3 * P_chp_max * u_chp]; Constraints = [Constraints, P_chp <= P_chp_max * u_chp];

目标函数就比较直白了:

C_fuel_chp = sum(P_chp / eta_chp * price_gas); % CHP耗气成本 C_fuel_gb = sum(P_gb / eta_gb * price_gas); C_grid = sum(price_buy .* P_buy) - sum(price_sell .* P_sell); C_om = sum(om_chp * P_chp) + sum(om_gb * P_gb) + ...; C_carbon = ...; % 按2.2节分段线性 C_dr = sum(c_dr * L_dr_inc) + ...; % 如果有DR补偿 Objective = C_fuel_chp + C_fuel_gb + C_grid + C_om + C_carbon + C_dr; ops = sdpsettings('solver', 'gurobi', 'verbose', 2); Diagnostics = optimize(Constraints, Objective, ops); if Diagnostics.problem ~= 0 disp('求解失败'); end

需要提醒的是,L_dr_price是通过弹性矩阵在优化前计算好的"固定值"还是"决策变量",要分清楚。如果作为决策变量,就必须额外加上弹性约束矩阵;如果提前算好,直接当作已知量代入平衡方程即可。两种都有人用,但混用会出错。

4.3 目标函数里非线性项的线性化处理

这个模型里最容易出现非线性项的地方有三处:

一是CHP热电耦合如果是运行域多边形,可能引入"出力点是否在多边形内部"这类约束,处理方式是顶点组合法或线性不等式组近似,不建议用二次约束。

二是1.1节提到的碳交易阶梯价格,用分段线性化。

三是储能充放电的功率损耗建模。如果写成P_es_c × η_es这样的形式还好,一旦要表达"充电时的损耗和放电时的损耗不对称",就必须用两个0-1变量互斥,再加辅助变量,避免出现P_es_c × P_es_d这种产品项。

凡是出现两个sdpvar相乘的地方,都要警惕。Yalmip会尝试自动处理部分情况,但代价是模型变成非线性半定规划(SDP),Gurobi就不认了。

5. 算例设计与结果解读:怎么让数字说话

5.1 测试系统构成与参数设置

我不建议一上来就搭一个巨型系统,先把一个中型园区模型跑通更有价值。我常用的算例配置如下:

设备容量/参数说明
CHP机组2 MW电出力,热电比0.9,电效率0.35主要电源
燃气锅炉4 MW,热效率0.9调峰热源
电锅炉1 MW,热效率0.95低碳热源
光伏+风电3 MW + 2 MW可再生出力曲线给出
电储能1 MW / 2 MWh,效率0.95SOC范围0.1~0.9
热储能2 MWh热力时移

电价采用峰平谷三段,比如峰时1.2元/kWh、平时0.8元/kWh、谷时0.4元/kWh,天然气价格取2.8元/m³并折算成单位热值价格。碳价基准设为50元/tCO2,排放因子按0.6 kgCO2/kWh(电网)和0.2 kgCO2/kWh(天然气)来设。

这些参数的选取没有放诸四海皆准的标准,关键是保持内部一致性。我见过不少论文把天然气价和碳价取到某个比例之后,结果里燃气设备从头到尾都不开机,这明显是参数失配。

5.2 三个场景怎么设置才有对比价值

我会固定所有设备参数和负荷曲线,只改动三个变量:

  • 场景A:不引入碳交易、不考虑DR,即传统经济调度;
  • 场景B:引入碳交易机制,但没有DR;
  • 场景C:碳交易和两类DR同时启用。

这个三场景对比的好处是能拆分两个机制的贡献。B和A对比看碳交易的减排效果和成本影响,C和B对比看DR在碳约束下的进一步优化能力。如果想把需求响应的效果单独拆出来,还可以加一个场景D:只有DR没有碳交易。四个场景一起跑,能展现的结论链条就非常完整了。

5.3 结果里最值得关注的三个指标

跑完优化后,我习惯先看三样东西,而不是直接贴出一堆曲线:

第一是总成本和成本构成变化。碳交易引入后总成本不一定上升——如果系统通过电锅炉替代燃气锅炉获得了碳收益,反而可能下降。DR引入后购电成本会降,但补偿成本会增加,要看净效果。

第二是碳排放总量的变化。这是碳交易建模的"验收指标",如果考虑碳交易后排放不降反升,那多半是配额系数给得太松,或者排放因子设置有问题。

第三是负荷曲线形态。把场景A和场景C的L_base+L_dr_price-L_dr_inc画出来,对比峰谷差缩小比例。这个指标能直观反映DR有没有起作用,也是报告里最有说服力的图之一。

我还会把每个时段的CHP出力、储能SOC、购电量画在一张图上,检查是否存在不合理的跳变——比如SOC突然从0.9掉到0.1再冲回0.9,那大概率是储能约束写错了。

6. 实际调试中的坑与排查思路

6.1 求解器报infeasible时的排查顺序

Yalmip/Gurobi返回不可行时,第一反应不要去看数学公式,先做三件事:

一是检查功率平衡约束的"方向"对不对。等号约束最容易因为符号习惯出错,尤其是把售电、储能放电、DR削减这些带方向性的量写反一个符号,模型立刻不可行。

二是检查储能SOC约束的初始条件。如果SOC(1)没赋初值,或者SOC(T)被限定为0.5,而T时刻前后出现了奇怪的强制关系,也很容易不可行。我习惯先松开末端SOC约束跑一遍,确认模型是"真不可解"还是"约束太紧"。

三是检查0-1变量的互斥约束。比如储能充放互斥,如果两个0-1变量之和≤1约束写成了≥1,系统要求每个时段必须充或放,在某个零边界时段就会冲突。

另一个定位技巧是逐组注释约束:先把全部约束注释掉,只跑目标函数,然后逐步放开约束组,看哪一组放开之后才变得不可行,问题就出在哪里。这个方法笨但极有效。

6.2 阶梯碳价0-1变量导致的非线性陷阱

碳交易阶梯成本如果实现不当,很容易引入sdpvar的乘积项。举个例子,有人在Yalmip里直接写:

C_carbon = 0.5 * p_c * lambda1 * z1 + 1.1 * p_c * (delta_E - lambda1) * z2;

这里z2是binvar,delta_E是sdpvar,两者相乘就是一个双线性项。Gurobi不能直接处理,Yalmip要么把它转为非凸二次规划,要么报错。解决方案就是我2.2节说的:把每个梯度区间的排放量拆成独立的非负连续变量E1、E2、E3,让0-1变量只和这些"切片变量"的上限绑定,而不是直接参与乘法。

此外,大M量的选择也容易出问题。M取太小会导致可行域被错误截断,M取太大会引起数值振荡。我通常按"该变量的物理上限×10"来定,E的物理上限就是最大排放量,算一下单位数量级再给M,不要随手写个1e6。

6.3 DR参数不合理导致的"伪优化"结果

需求响应参数如果设得太激进,会出现一种看起来漂亮但不真实的调度结果:系统通过DR把大量负荷挪到谷时,然后谷时电价和碳排放同时很低,设备全开,总成本大幅下降,峰谷差几乎抹平。现实中这不可能,因为用户的负荷转移意愿有限。

我的经验是价格型DR的弹性系数绝对值不要超过0.5,激励型DR的中断容量不要超过峰值负荷的15%,中断次数限制在2~3次以内。另外,DR削减的负荷总量要监控:削减时段负荷下降,转移时段负荷上升,但如果转移后的负荷超过了设备总出力上限,又会出现不可行。这时候要回头检查弹性矩阵配置是否合理,而不是盲目加设备容量。

最后分享一个我自己的调试习惯:每跑完一组场景,我会把SOC曲线、DR削减量曲线和碳成本曲线同时画出来看。这三个量的行为能反映模型里90%的逻辑错误。调度优化这东西,模型能跑通只是万里长征第一步,结果能不能解释得通、物理上是否自洽,才是最花时间的地方。希望这篇分享能帮你少走一段弯路。

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询