☰
区域综合能源系统低碳经济调度中的主从博弈分层优化与Matlab实现
2026/9/30 4:01:16 网站建设 项目流程

如果你做过区域综合能源系统的调度优化,大概率经历过这种尴尬:模型建得很漂亮,结果一跑,总是有人不满意。运营商觉得利润太低,用能用户觉得成本太高,两边博弈出的调度方案怎么也落不下去。原因挺简单——传统集中式优化默认“所有设备属于同一个大脑”,但实际区域里每一方都有自己的利益诉求。所以这几年我越来越倾向于用多主体主从博弈来构建区域综合能源系统低碳经济优化调度模型,并在Matlab里把它落地成了一套可复现的分层模型求解代码。这篇文章就从模型结构、数学推导、代码实现和调试踩坑几个角度,把完整思路给你捋一遍。

1. 为什么要用多主体主从博弈,而不是集中式优化

1.1 集中式调度的两难

区域综合能源系统里通常有电、气、热、冷多种能源耦合,设备侧有热电联产机组(CHP)、燃气锅炉、电制冷机、储能、P2G等,负荷侧有居民、商业、工业等多种用户。早期做这类调度,最常见做法是把整个区域当成一个整体,用一个目标函数(比如总成本最小)一次性求解所有设备的出力计划。这种思路在数学上很干净,但工程上经常被挑战:区域能源运营商和用户根本不是同一个利益主体。

举个最直观的例子:运营商为了降低自身购电成本,倾向于在低谷时段大幅压低售电价格引导用户多用电,但用户不一定愿意配合,因为改变自己的用能习惯可能带来生产不便。反过来,用户希望在尖峰时段自己多发电减少购电,但运营商售电收入就会减少,双方利益直接冲突。集中式模型把这种冲突“压平”了,求出来的解是“上帝视角”下的事故调度,执行时合作方根本不接受。

1.2 领导者-跟随者结构:运营商的定价权与用户的响应权

主从博弈(Stackelberg博弈)正好适合描述这种“先决策、后响应”的关系。在区域综合能源系统里,能源运营商天然处于领导者地位,因为它掌握电网、气网的接入权,能制定售电、售热价格以及需求响应补贴策略。用户作为跟随者,在看到价格信号后,根据自己的设备条件调整购能计划和自发电计划。

两者之间是一个典型的双层决策过程:

  • 上层(领导者):运营商决策各时段售电、售热价格、购能计划、碳交易策略,目标是自身利润最大化;
  • 下层(跟随者):各类用户决策从运营商购买的电、热功率,以及自有设备(光伏、储能、小型CHP等)的出力,目标是自身综合用能成本最小化。

这种结构天然就是分层模型:上层优化结果通过价格传递给下层,下层优化结果通过负荷需求反馈到上层,反复迭代直到双方都没有单方面改变策略的动机。这时候得到的解才是真正可落地的纳什均衡解。

1.3 分层模型在工程上到底意味着什么

很多刚接触的人会把“分层”理解成“先算上层再算下层”的简单串联,这其实是误区。分层模型的核心在于上下层变量互相耦合,不能分开独立求解。上层改变价格,下层需求就变;下层需求变了,上层利润就变。两者之间是一个闭环。

我做这个项目时最深刻的体会是:分层模型不是算法的选择,而是现实关系在数学结构上的映射。区域里产权边界清晰、运营主体多元,用主从博弈描述就是最自然的建模方式。后面搭建的低碳经济调度目标、碳排放约束、设备运行约束,全都依附在这个“领导者-跟随者”骨架之上。

2. 低碳经济调度的目标与约束怎么形式化

2.1 上层决策变量与经济目标

运营商的决策变量分两类:价格类变量和计划类变量。

价格类变量包括每个时段的售电价格、售热价格,以及给参与需求响应用户的补偿单价。计划类变量包括从上级电网购电功率、购买天然气量,以及碳交易量。目标函数可以写成:

利润 = 售电收入 + 售热收入 - 购电成本 - 购气成本 - 设备运维成本 - 碳排放成本

如果用数学公式表达,大概是这样:

max Σ(rP_e,t·P_load_e,t + rH_h,t·H_load_h,t) - Σ(rGrid_t·P_buy,t) - Σ(rGas·V_gas,t) - C_om,t - C_CO2

这里有个细节容易忽略:用户在低谷用不用电、用多少电,取决于价格,所以上层利润函数里收入项rP_e,t·P_load_e,t包含了上下层耦合变量,两者都是决策变量而且相乘,这是一个双线性项。后期处理时需要特别注意。

2.2 下层用能主体的运行模型

下层模型我会按“典型主体”来建,比如工业用户、商业园区、居民小区各有各的设备构成和用能特征。拿一个有代表性的综合用户来说,它的可调资源包括:

  • 热电联产机组(CHP),电热联供;
  • 电锅炉或热泵,用于补充供热;
  • 储能电池,可充可放;
  • 一部分可平移负荷(比如工业流水线、空调等)。

下层目标函数是用户总成本最小:

min Σ(购电费用 + 购热费用 + 设备运行成本 - 需求响应补贴) - 用能满意度收益

约束条件需要包含能量平衡约束:

P_load_e,t = P_buy_user,t + P_chp_e,t + P_battery_discharge,t - P_battery_charge,t

H_load_h,t = H_buy_user,t + Q_chp_h,t + Q_boiler,t

以及设备出力上下限约束、爬坡约束、储能SOC循环约束等。如果用户侧需求响应要建模为可转移负荷,还需要引入时间顺序约束和时段迁移的整数变量:

sum(Δt·P_shiftable,(start,end)) = E_shiftable_total

这类整数变量在下层模型中出现,会让问题变成一个混合整数线性规划(MILP),后续用KKT条件转化时也会带来对应的对偶变量。

2.3 阶梯碳交易与碳排放约束的实现

低碳经济调度不是说“碳排越少越好”,而是要在经济性和低碳性之间找平衡。我做这个模型时没有把碳排放当成硬约束,而是用阶梯碳交易机制,让碳成本随排放量增长而加速,这样模型会自发地权衡减碳投资和经济成本。

阶梯碳交易成本函数一般写成:

C_CO2 = λ_c·E_carbon,且当E_carbon处于不同区间时,λ_c采用不同档位。

比如:

  • 第一档:碳排放在配额Q以内时,按基准碳价λ1计算;
  • 第二档:超过配额但小于1.2倍配额时,按λ2 = 1.2λ1计算;
  • 第三档:超过1.2倍配额时,按λ3 = 1.5λ1计算。

这里碳排放量E_carbon怎么算也很关键。我采用的是从上级电网购电对应的间接碳排放加天然气直接燃烧排放,再扣减P2G设备或碳捕集装置的固碳量:

E_carbon = Σ(α_grid·P_buy,t + β_gas·V_gas,t) - E_capture

阶梯碳价是非线性分段函数,在Matlab里用整数变量和Heger不等式转化成线性约束。这个转化在求解器眼里是最基本的操作,但写错的人非常多,后面我会专门讲坑。

3. 双层优化问题的求解路径:从KKT到MILP

3.1 把下层问题用KKT条件“吸”到上层

上下层问题不能分开独立求解,就必须用数学手段把下层问题“嵌入”到上层问题里。最经典的做法就是KKT条件替换法——把下层优化问题的最优性条件作为约束加到上层问题中。

下层用户的目标是线性或凸的优化问题,当它满足凸性和约束规范(比如Slater条件)时,KKT条件是原问题的充分必要条件。下层KKT条件包含四组:

  • 拉格朗日函数对下层决策变量求导为零(驻点条件);
  • 原问题约束可行(原始可行性);
  • 对偶变量不小于零(对偶可行性);
  • 对偶变量与对应松弛变量乘积为零(互补松弛条件)。

把这些条件全部加到上层模型里,双层优化就变成了带互补约束的单层数学模型。这时模型在数学上等价于原来的双层问题,但可以直接交给优化求解器处理。

3.2 强对偶松弛:避坑利器

KKT方法数学上干净,但实际求解时会带来非常多整数变量,问题规模大了以后求解速度非常不理想。尤其当用户侧含储能、可平移负荷这些带整数变量的设备时,下层本身就是个MILP,再对MILP取KKT条件就困难了——整数变量的对偶理论不像连续问题那么直接。

我实际项目里更偏好用强对偶松弛来替代其中的一部分。思路是这样:下层问题是线性规划时,原问题的最优值等于对偶问题的最优值。利用这个性质,可以把下层目标函数值用对偶变量表达,从而消掉上层目标函数和约束里的双线性项,也就是价格变量和电量变量的乘积。

具体来说,下层目标函数中所有“价格×电量”的乘积项,对偶后会变成“对偶变量×常量系数”的形式。如果下层满足强对偶条件,这些项就可以完全改写。这个技巧在实际计算中非常有效,因为普通KKT方法处理双线性项只能通过引入大量辅助变量,而强对偶松弛直接绕过了一部分双线性项。

3.3 互补松弛与双线性项的线性化

无论用KKT还是强对偶,最后一定会碰到两个需要线性化的地方:

第一个是互补松弛条件。对于类似λ·(g_max - g) = 0这种乘积为零的形式,标准做法是引入二进制变量z和足够大的常数M:

λ ≤ M·z

g_max - g ≤ M·(1 - z)

这里的M取值非常讲究,如果太小会错误切割可行域,导致求解器直接报“Infeasible”;如果太大又会让松弛空间过大,增加求解时间。我一般先用连续松弛版本算一遍,看对偶变量的大致量级,然后取M为这个量级上限的10倍以上,再微调。

第二个是目标函数中的双线性项,比如上层收入项rP_e,t·P_load_e,t。当用户侧模型被KKT或对偶改写后,这些项往往可以消掉。如果消不掉,就需要做分段线性化或McCormick包络近似。我的经验是:尽量通过强对偶重构把双线性项的结构性消掉,实在消不掉再考虑近似方法,因为任何近似都会破坏均衡解的严格性。

4. Matlab代码架构与核心实现要点

4.1 整体求解流程与代码目录

当你把双层问题转化成单层MILP之后,Matlab里的实现就回到了“建模-求解-后处理”的标准流程。我习惯用Yalmip做建模层,求解器用Gurobi或Cplex,效率比直接用linprog或fmincon高一个量级。整个代码目录大概是这样的:

RIES_Stackelberg/ │ main.m % 主程序:参数初始化、建模、求解、结果输出 │ case_data.m % 算例参数:负荷曲线、能源价格、设备参数 │ build_upper_model.m % 上层目标函数与约束 │ build_lower_kkt.m % 下层问题KKT条件生成 │ linearize_complement.m % 互补松弛线性化 │ plot_results.m % 结果可视化

main.m的推进逻辑很固定:先加载数据,然后生成上层变量和下层的原始变量、再调用build_lower_kkt将下层问题转为KKT系统、线性化、组合成单层模型、设置求解器选项,最后用optimize求解并取出结果。

4.2 Yalmip建模的关键片段

上层变量定义示范:

T = 24; nb = 3; % 用户数量 % 上层变量:运营商价格 rho_e = sdpvar(1, T, 'full'); % 售电价格 rho_h = sdpvar(1, T, 'full'); % 售热价格 P_buy = sdpvar(1, T, 'full'); % 从上级电网购电 V_gas = sdpvar(1, T, 'full'); % 购气量 % 下层变量:每个用户的购电量和购热量 P_load_user = sdpvar(nb, T, 'full'); H_load_user = sdpvar(nb, T, 'full'); % 下层自设备变量(行列对应用户数和时段数) P_chp = sdpvar(nb, T, 'full'); E_sto = sdpvar(nb, T, 'full'); % 储能充放电净功率

下层KKT条件如果手动写,代码会比较长。一个简化做法是直接用Yalmip内置的kkt函数:

% 定义下层问题 x_user = [P_load_user(:); H_load_user(:); P_chp(:); E_sto(:)]; % 注意:这里要求下层目标是线性/二次凸,约束为线性 Constraints_user = [ ... ]; Objective_user = 购电成本 + 购热成本 + 设备成本 - 需求响应补贴; [KKT_sys, details] = kkt(Constraints_user, Objective_user, x_user);

kkt函数会自动生成驻点条件、互补松弛条件和整数变量标记,省去手动编写大量拉格朗日求导的环节。但要注意,kkt函数只适用于凸问题和线性目标,使用前务必确认下层模型是凸的。我早期的设备模型里带了一个凹的供能效率函数,结果KKT生成后求解器老是报错,排查了很久才意识到是凸性条件没满足。

4.3 求解器配置与收敛性控制

单层MILP模型规模通常比较大,变量数几千、约束数接近上万很正常。求解器选项设置直接决定能不能在合理时间拿到解。我习惯的配置如下:

options = sdpsettings('solver', 'gurobi', ... 'gurobi.MIPGap', 1e-3, ... % 相对MIP间隙,1e-3足够工程用 'gurobi.TimeLimit', 3600, ... 'gurobi.NumericFocus', 2, ... 'savesolveroutput', 1);

NumericFocus是我特别关注的一项。因为互补松弛线性化引入了大M,模型数值条件普遍不好,把NumericFocus开到2或3能让求解器多花时间做预处理和数值稳定处理,很多“求解器报数值错误”的问题都能这样解决。

还有一个容易被忽略的点:初始化。MILP里二进制变量数量庞大,给一个靠谱的初始解能大幅缩短求解时间。我的做法是先用固定典型价格传给下层模型求出用户响应,再用这个响应结果作为MILP的热启动点。Gurobi会通过start属性接收部分变量的初始值。

5. 典型算例结果:调度曲线与低碳经济性

5.1 算例场景设置

我给自己搭的典型算例设了1个区域能源运营商和3类用户,24小时调度周期。上级电网采用分时购电价,天然气价格固定。用户1是工业用户,带小型CHP和可平移负荷;用户2是商业园区,带储能和电锅炉;用户3是居民小区,以纯电负荷为主,热负荷用燃气锅炉满足。

碳排放配额按历史负荷乘以配额系数确定,基准碳价设0.25元/kg,阶梯系数按1.0、1.2、1.5倍递进。整个模型求解规模大概是连续变量3200个,二进制变量480个,Gurobi在MIPGap=1e-3下运行约280秒收敛。

5.2 静态策略对比与结果解读

作为对照,我还跑了传统分时电价下的集中式调度模型。两种模式的核心结果对比如下:

指标集中式分时电价主从博弈分层模型
运营商总利润(元)84209135
用户总用能成本(元)2068019540
系统碳排放(kgCO2)1486013620
尖峰负荷(kW)12501140
谷时负荷(kW)620710

从结果能看出几个有意思的现象,也可以说是主从博弈模型最典型的行为特征:

第一,运营商利润和用户成本不是零和博弈。运营商调低峰时段电价、拉高尖峰时段电价,引导用户把可平移负荷从尖峰挪到谷段,尖峰购电成本大幅下降,这部分利润增量足以弥补峰时少量售电量下降带来的收入损失。用户因为整体购电结构优化,总成本也降了。这就是主从博弈和集中式模型的本质区别——它找到的不是“总成本最小”,而是“双方都不吃亏”的均衡。

第二,低碳指标明显改善。碳排下降接近8.4%,主要来自两个渠道:一是CHP更倾向于在气价相对购电价有利的时段满发,减少从电网买高碳火电;二是用户侧储能充电时机会更贴近谷时低价电,整体电量结构变“绿”。值得注意的是,这个改善不是靠硬约束逼出来的,而是阶梯碳价让运营商主动把碳成本纳入价格策略的结果。

第三,负荷曲线被“削峰填谷”了。尖峰负荷从1250kW降到1140kW,谷时负荷从620kW升到710kW,日负荷率明显提升。对于配电网来说这个效果很有价值,因为它意味着上级电网扩容压力变小,从区域综合能源系统整体来看也是一种低碳贡献。

6. 我踩过的坑与后续扩展建议

6.1 大M值选取的教训

第一次把下层KKT互补松弛线性化后,我直接拍脑袋选了M=1e5,结果模型要么Infesible要么解出来的价格在几个时段跳得极其不合理。排查后才反应过来:M太大时二进制变量的约束几乎不起作用,对偶变量可以在大范围里自由漂移,数值解实际上被“松弛”坏了。后来我用连续松弛先求解一遍,把对偶变量量级摸清楚——一般在10左右——然后取M=50,模型稳定性立刻上去了,求解时间也降了一个数量级。

6.2 交替迭代求解的稳定性问题

项目初期我还试过另一种流行思路:上层用粒子群或遗传算法,下层用线性规划,两层交替迭代直到收敛。这种方法代码写起来直观,但实际表现很不稳定,经常在几个方案之间来回震荡,很难判定是否真的到达均衡。后来换成KKT/强对偶的单层MILP转化,问题一次性求解,均衡解的定义就严格多了。

我的建议是:除非模型规模大到单层转化无法求解,否则尽量用单层化方法。如果实在要用交替迭代,至少要加解一致性约束和动态惯性权重,否则就是在碰运气。

6.3 扩展方向:多领导者、不确定性和碳捕集耦合

这套分层模型框架的可扩展性很好,我目前正在做两个方向:第一是把单个领导者扩展成“多领导者-多跟随者”结构,比如同时存在多个综合能源运营商时,上层需要通过更复杂的均衡条件耦合;第二是引入风光出力的不确定性和分布式鲁棒优化,让价格策略在风光波动时依然可靠。

碳捕集与P2G设备的耦合也是我很看好的方向。区域综合能源系统的低碳价值不仅体现在用电结构上,还要考虑二氧化碳捕集后用于制天然气循环利用。把这个环节纳入阶梯碳交易体系后,上层决策变量会多出捕集率、储碳量等维度,模型闭环会更完整——只是求解规模又上一个台阶,这也是下一步最头疼也最值得投入的地方。

最后分享一个我个人的实操建议:不要一开始就追求模型面面俱到。先把一个最简单的主从博弈结构跑通——上层只做分时电价,下层只有一种用户一台CHP,理解透KKT和强对偶的转化关系后再逐步加储能、加阶梯碳价、加多类用户。这样每加一层,你都能明确知道是哪些约束在影响均衡结果,而不是面对一个几千行代码的黑盒模型束手无策。我的这套Matlab框架就是从这个小版本一步步长起来的,最大的教训就是:分层模型难的不是模型本身,而是每一层之间的因果链条必须清清楚楚。

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

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

立即咨询