电热综合能源系统的动态定价说穿了就是“谁来定价、怎么定价、定了之后用户听不听”的问题。我最早做这个项目的时候用传统的集中式优化把整个园区的电、热一起调度模型收敛得很快结果拿到实际场景里却没人买账——园区用户不是整体每个楼宇、每条产线都有自己的利益诉求运营商定的价格如果只让运营商赚钱用户就会用脚投票。后来我换成主从博弈框架把运营商设为上层领导者、用户设为下层跟随者问题才真正立起来。这篇文章就把我基于Matlab实现的一套“主从博弈电热综合能源系统动态定价能量管理”的完整思路和踩坑记录分享出来适合正在做综合能源优化调度、或者想用博弈论方法做定价研究的同学参考。1. 为什么电热综合能源定价会变成一个博弈问题1.1 运营商和用户的利益不一致才是定价问题的本质很多人在第一次接触“综合能源系统优化”时第一反应是建一个大目标函数把购电成本、购气成本、设备运维成本、用户满意度全部加权求和然后一次性求最优。这种做法的隐含假设是“整个系统是一个利益共同体”——但现实完全不是这样。综合能源系统中至少存在两类主体一类是综合能源服务商或者叫能源运营商它拥有CHP机组、电锅炉、储能设备负责从上级电网买电、从气网购气再把这些能源加工成电和热出售给用户另一类是终端用户它们关心的是用能成本最小化会随着电价、热价的变化调整自己的用电用热行为。这两类主体的目标函数天然冲突运营商想让售能收入最大化用户想让购能成本最小化。如果运营商盲目抬高价格用户会削减用能运营商的收入反而可能下降如果运营商压低价格用户用能增加但单位利润变薄。这是一个典型的利益博弈关系不是单个优化问题能刻画清楚的。1.2 主从博弈为什么天然适配“先定价、后响应”的模式主从博弈Stackelberg Game的核心思想是博弈双方的地位不对等一方先行动领导者另一方观察到领导者的行动后再做决策跟随者。领导者做决策时必须把跟随者的最优反应函数纳入自己的优化模型中。放到电热综合能源系统里这个结构非常自然领导者上层综合能源服务商。先制定下一个调度周期内比如24小时的电价、热价它的目标是净收益最大化扣除购能成本和设备运行成本。跟随者下层用户。给定电价和热价后用户调整自己的购电、购热计划目标是用能成本最小化同时满足自身的用能舒适度约束。这种“先定价、后响应”的时序关系和现实中综合能源服务商发布分时电价、用户根据价格调整空调温度设定或电动车充电时段的场景完全一致。1.3 和集中式优化比主从博弈模型带来了什么新信息集中式优化得到的结果是一个系统层面的“全局最优”但它回答不了“这个最优是怎么通过市场机制实现的”这个问题。主从博弈模型则给出了一个价格信号形成的机制上层先定价格下层做出最优响应上层再根据响应调整价格直到达到均衡。最终得到的电价、热价曲线是博弈均衡解而不是单纯的最优解。举个例子集中式优化可能得到“在某个时段大量用电”的调度方案但它不会告诉你电价应该是每千瓦时多少——电价是输入给定的。主从博弈模型却把“价格”本身作为决策变量这正好对应了电力市场改革后售电公司自主定价的现实需求。当然主从博弈也有代价模型是非凸的双层优化问题数学上很难直接求解需要做KKT条件转换、线性化等一系列处理。这也是这篇文章后面重点展开的部分。2. 电热耦合系统建模先把物理图景搭起来要写代码先得把数学模型写清楚。主从博弈不是凭空套壳下层用户能优化出什么东西取决于上层运营商的系统里有哪些设备、怎么运行。我先给出一套标准的电热综合能源系统拓扑和对应的数学模型。2.1 一个典型的电热综合能源系统长什么样我做项目时一般假设这样一个园区级系统调度周期为24小时时间尺度为1小时能源输入上级电网购电、天然气网购气。能量转换设备CHP机组热电联产同时产电产热、燃气锅炉或电锅炉。储能设备蓄电池、蓄热罐。负荷侧电负荷、热负荷由多个用户楼宇聚合而成。系统的能量流可以简化描述为天然气进入CHP机组后同时输出电功率和热功率燃气锅炉补充供热电网购电和CHP发电一起供应电负荷多余的可以充入蓄电池热负荷由CHP余热、锅炉供热和蓄热罐放热共同满足。热网部分可以做一个合理简化不细究管网拓扑和流体动态把热网看作一个传输环节主要考虑热功率平衡和蓄热罐的储放热约束。对于定价研究来说这种精度的模型已经够用细节热网模型会让博弈问题复杂到很难求解。2.2 设备模型CHP机组、电锅炉、储能、新能源的数学描述下面给出我常用的设备模型这些方程会直接进入上层优化问题。CHP机组CHP机组的热电耦合特性可以用“热电比”来刻画。最简单也最常用的模型是固定热电比模型电功率输出P_CHP(t)热功率输出H_CHP(t) r_CHP × P_CHP(t)其中r_CHP是热电比常数。更精细的做法是用一个可行域来描述CHP的运行范围这时H_CHP(t)和P_CHP(t)之间存在线性不等式约束。我在项目中用的是后者因为固定热电比的模型太理想无法反映机组在部分负荷下的调节能力。约束如下P_CHP(t) 在最小技术出力P_CHP_min和最大出力P_CHP_max之间。热出力与电出力满足上、下爬坡约束限制机组在相邻时段的出力变化量。燃料成本可以用二次函数近似但为了后续线性化处理会分段线性化。电锅炉/燃气锅炉电锅炉H_EB(t) η_EB × P_EB(t)输入电功率P_EB(t)输出热功率H_EB(t)。效率η_EB一般取0.9左右。燃气锅炉H_GB(t) η_GB × F_GB(t)输入天然气功率F_GB(t)输出热功率H_GB(t)。效率η_GB一般取0.85~0.95。蓄电池用荷电状态SOC来描述SOC(t1) SOC(t) η_ch × P_ch(t) × Δt / E_cap - P_dis(t) × Δt / (η_dis × E_cap)其中P_ch、P_dis是充放电功率η_ch、η_dis是充放电效率E_cap是额定容量。充放电功率都有上下限SOC也在[0.1, 0.9]之间且调度周期结束时SOC要回到初始值这是为了保证储能运行的可持续性。蓄热罐蓄热罐的模型和蓄电池非常像用储热量S(t)来描述S(t1) S(t) H_ch_st(t) × Δt - H_dis_st(t) × Δt蓄热罐的提升在于“热惰性”——热负荷高峰时可以通过放热支撑热负荷低谷时吸收多余热量。在动态定价问题中蓄热罐的配置会显著影响运营商的定价策略。2.3 用户侧建模价格弹性、热负荷特性、需求响应约束下层用户不是简单的“固定负荷曲线”而是会对价格做出响应。建立下层模型时我需要描述用户如何根据电价和热价调整用能。我采用一个常见的需求响应模型用户的总用能效用函数最大化。比如用户用电和用热的效用可以用二次效用函数描述U_e(P_L(t)) a_e × P_L(t) - 0.5 × b_e × P_L(t)^2U_h(H_L(t)) a_h × H_L(t) - 0.5 × b_h × H_L(t)^2那么用户的决策问题就是给定电价ρ_e(t)和热价ρ_h(t)选择用电功率P_L(t)、用热功率H_L(t)最大化总效用减去购能成本max Σ_t [U_e(P_L(t)) - ρ_e(t) × P_L(t) U_h(H_L(t)) - ρ_h(t) × H_L(t)]s.t. P_L_min ≤ P_L(t) ≤ P_L_maxH_L_min ≤ H_L(t) ≤ H_L_max这里的上下限就体现了用户的刚性需求哪怕电价再高基本负荷也得用哪怕热价再低热负荷也不能超过散热器或地暖系统允许的上限。这个二次效用函数有一个非常好的性质用户的最优响应是价格的线性函数——电价升高用电量线性下降。这一点是后面通过KKT条件把下层问题变成上层约束的关键基础。2.4 上层目标函数和价格约束上层运营商的优化目标是调度周期内净收益最大化max Σ_t [ρ_e(t) × P_L(t) ρ_h(t) × H_L(t) - C_buy(t) - C_fuel(t) - C_OM(t)]其中ρ_e(t) × P_L(t) 是售电收入ρ_h(t) × H_L(t) 是售热收入C_buy(t) 是向上级电网购电的成本按分时电价结算C_fuel(t) 是天然气购气成本C_OM(t) 是设备运维成本价格不能乱定要加约束ρ_e_min ≤ ρ_e(t) ≤ ρ_e_maxρ_h_min ≤ ρ_h(t) ≤ ρ_h_max这个约束代表的是监管要求或市场规则运营商不能无限抬价。做上层模型时如果不加价格上限优化算法很容易给出一个极高的价格你一看结果就知道这个模型跑飞了——下层用户的负荷会被压到下限整个系统失去意义。3. 双层模型的数学转换KKT条件与强对偶怎么用模型建好之后摆在面前的是一个双层优化问题上层运营商定价格最大化收益下层用户根据价格最小化用能成本等价于最大化效用减去成本这个双层问题不能直接丢给求解器。求解器能处理的要么是线性规划、要么是混合整数线性规划顶多是非线性规划但双层结构它不认识。标准的处理方法是把下层问题用KKT条件替换把双层问题转换成单层问题。3.1 先把双层问题写成标准形式为了推导方便我把下层用户问题简化表示为一个凸优化问题max_{x} f(x, ρ)s.t. g(x) ≤ 0其中ρ是上层给定的价格向量x是用户的用能向量f是用户的效用函数g是不等式约束集合。用户问题是凸的二次效用函数的Hessian矩阵是负定的约束是线性的所以KKT条件是用户问题的充要条件。这意味着把KKT条件加上就等价于把用户问题“嵌入”了上层模型。3.2 下层问题的KKT条件推导以下层用户的功率平衡约束为例。状态变量是P_L(t)和H_L(t)约束是上下限。我对每个不等式约束引入拉格朗日乘子λ_min、λ_max写出拉格朗日函数再对决策变量求偏导。用户问题的KKT条件包括三部分平稳性条件拉格朗日函数对每个决策变量的偏导数为0。对电功率P_L(t)∂L/∂P_L(t) a_e - b_e × P_L(t) - ρ_e(t) λ_e_min(t) - λ_e_max(t) 0对热功率H_L(t)同理∂L/∂H_L(t) a_h - b_h × H_L(t) - ρ_h(t) λ_h_min(t) - λ_h_max(t) 0原始可行性原问题的约束必须成立也就是P_L_min ≤ P_L(t) ≤ P_L_max、H_L_min ≤ H_L(t) ≤ H_L_max。互补松弛条件每个不等式约束和它对应的乘子不能同时“松弛”。写成数学形式就是λ_e_min(t) × (P_L_min - P_L(t)) 0λ_e_max(t) × (P_L(t) - P_L_max) 0对热负荷约束同理。乘子和约束的松弛量都是非负的。3.3 互补松弛的大M线性化互补松弛条件是非线性的两个变量相乘等于0而且包含整数逻辑求解器没法直接处理。行业里的标准做法是引入二进制变量用大M法线性化。比如对λ_e_min(t) × (P_L_min - P_L(t)) 0我引入二进制变量z1(t) ∈ {0,1}写成两个约束λ_e_min(t) ≤ M × z1(t)P_L_min - P_L(t) ≤ M × (1 - z1(t))这里M是一个足够大的正数大M。逻辑是如果z1(t)1则第一个约束给乘子上限第二个约束强制P_L(t)P_L_min即约束紧满足互补松弛如果z1(t)0则强制λ_e_min(t)0第二个约束被松弛到恒成立。每个互补松弛条件需要引入一个二进制变量。对于24小时的调度问题引入的二进制变量数量是可控的每类约束24个总共不过几十个到上百个二进制变量Cplex/Gurobi处理这个规模的混合整数规划毫无压力。3.4 非线性目标项的强对偶消去把KKT条件替换进上层模型后上层目标函数里仍然有ρ_e(t) × P_L(t)和ρ_h(t) × H_L(t)这类乘积项——价格变量和负荷变量相乘这不是线性的。怎么处理用强对偶定理。因为下层用户问题是凸优化强对偶成立下层问题的原目标最优值等于对偶问题的最优值。具体来说用户问题的目标函数是Σ_t [a_e × P_L(t) - 0.5 × b_e × P_L(t)^2 - ρ_e(t) × P_L(t) a_h × H_L(t) - 0.5 × b_h × H_L(t)^2 - ρ_h(t) × H_L(t)]它的对偶函数可以通过拉格朗日函数对决策变量求极大推导出来。利用强对偶我可以把上层目标函数中非线性项ρ_e(t) × P_L(t)替换成用户效用函数加上对偶函数的形式从而实现对目标函数的线性化。具体推导过程比较长这里直接给结论在引入KKT条件后上层目标函数中每个时段的收入和用户负荷的乘积项可以用用户效用项与拉格朗日乘子的内积等价替换。这样整个模型变成一个混合整数线性规划MILP可以直接用MatlabYalmipCplex求解。这一步是整个转换中最容易出错的地方也是很多论文里“推导略”的部分。强烈建议自己动手推一遍不要直接抄结果。4. MatlabYalmip代码实现从数学模型到可运行程序模型转换完之后编码实现就是体力活了。我用的工具链是Matlab Yalmip Cplex。下面把程序框架和关键实现细节列出来这套代码跑通之后你就可以在它的基础上改设备参数、改价格约束、换用户效用函数。4.1 为什么选YalmipCplex而不是手写求解器如果你刚接触优化建模我强烈建议先别自己写求解算法。主从博弈问题经过KKT转换后是一个MILP手写单纯形法或分支定界法的工作量完全不可接受。Yalmip是Matlab下面的一个建模工具箱它把优化建模抽象成“声明变量、写约束、调求解器”三步极大降低出bug的概率。求解器的选择上Cplex和Gurobi都是MILP领域的标配。我习惯用Cplex主要原因是它在学术界的license申请方便、和Yalmip的接口稳定。如果你有Gurobi的license完全可以直接换成gurobi代码不需要改太多。4.2 程序整体框架和模块划分整个程序我分成四个模块各司其职参数初始化模块定义24小时的电网分时购电价、天然气价格、CHP机组参数、储能参数、用户效用函数参数、价格上下限。变量声明模块定义上层的价格变量、设备出力变量、储能SOC变量以及下层用户的负荷变量和经KKT转换引入的拉格朗日乘子变量、二进制变量。约束组装模块把设备运行约束、用户KKT条件约束、互补松弛线性化约束逐一写入Yalmip的约束集合。模型求解与结果输出模块调用Cplex求解MILP解出后把电价、热价、用户负荷、各设备出力曲线保存成结构体或表格方便后续画图和对比。下面是一个关键代码片段用来创建上层运营商的决策变量% 时间尺度 N 24; % 上层变量 rho_e sdpvar(1, N, full); % 电价决策变量 rho_h sdpvar(1, N, full); % 热价决策变量 P_CHP sdpvar(1, N, full); % CHP电出力 H_CHP sdpvar(1, N, full); % CHP热出力 P_EB sdpvar(1, N, full); % 电锅炉电功率 SOC sdpvar(1, N1, full); % 蓄电池SOC含初始状态 % 下层变量跟随者的负荷决策 P_L sdpvar(1, N, full); % 电负荷 H_L sdpvar(1, N, full); % 热负荷 % KKT拉格朗日乘子 lambda_e_min sdpvar(1, N, full); % 电价相关下限乘子 lambda_e_max sdpvar(1, N, full); % 电价相关上限乘子 lambda_h_min sdpvar(1, N, full); lambda_h_max sdpvar(1, N, full); % 互补松弛线性化需要的二进制变量 z_e_min binvar(1, N, full); z_e_max binvar(1, N, full); z_h_min binvar(1, N, full); z_h_max binvar(1, N, full);4.3 设备约束和KKT约束怎么往Yalmip里写设备约束直接按前文的数学公式写注意SOC的索引错位问题。Yalmip支持用循环写约束也可以用矩阵形式批量写建议把每一天的约束写在一个for循环里逻辑清晰一点。下面是我写约束的示例这里只列出关键部分Constraints []; % 价格上下限 Constraints [Constraints, rho_e_min rho_e rho_e_max]; Constraints [Constraints, rho_h_min rho_h rho_h_max]; % CHP出力约束 for t 1:N Constraints [Constraints, P_CHP_min P_CHP(t) P_CHP_max]; Constraints [Constraints, H_CHP(t) r_CHP * P_CHP(t)]; end % 蓄电池SOC递推 Constraints [Constraints, SOC(1) SOC_0]; for t 1:N Constraints [Constraints, SOC(t1) SOC(t) ... - P_dis(t)*dt/eta_dis/E_cap eta_ch*P_ch(t)*dt/E_cap]; Constraints [Constraints, SOC_min SOC(t1) SOC_max]; end % 功率平衡约束 for t 1:N Constraints [Constraints, P_grid(t) P_CHP(t) P_dis(t) - P_ch(t) P_L(t)]; Constraints [Constraints, H_CHP(t) H_EB(t) H_dis_st(t) - H_ch_st(t) H_L(t)]; endKKT条件部分先把平稳性条件写进去% KKT平稳性条件 for t 1:N Constraints [Constraints, ... a_e - b_e * P_L(t) - rho_e(t) lambda_e_min(t) - lambda_e_max(t) 0]; Constraints [Constraints, ... a_h - b_h * H_L(t) - rho_h(t) lambda_h_min(t) - lambda_h_max(t) 0]; end然后是互补松弛条件的大M线性化M 1000; % 大M需要根据问题规模评估 for t 1:N % 电负荷下限互补条件 Constraints [Constraints, lambda_e_min(t) M * z_e_min(t)]; Constraints [Constraints, P_L_min - P_L(t) M * (1 - z_e_min(t))]; % 电负荷上限互补条件 Constraints [Constraints, lambda_e_max(t) M * z_e_max(t)]; Constraints [Constraints, P_L(t) - P_L_max M * (1 - z_e_max(t))]; % 热负荷同理 Constraints [Constraints, lambda_h_min(t) M * z_h_min(t)]; Constraints [Constraints, H_L_min - H_L(t) M * (1 - z_h_min(t))]; Constraints [Constraints, lambda_h_max(t) M * z_h_max(t)]; Constraints [Constraints, H_L(t) - H_L_max M * (1 - z_h_max(t))]; end大M的取值要特别小心。M太小会把最优解卡掉M太大又会导致数值病态、求解精度下降。我的做法是把用户负荷的上下限差值作为基准再乘上价格上限估算出一个合理的M一般取M 2 × max(价格上限) × max(负荷范围)。在Cplex里还可以设置数值精度的相关参数来缓解大M带来的问题。4.4 两种求解思路KKT合并一次性求解 vs 迭代优化上面给的是“把KKT条件写进上层模型、一次性求解MILP”的方案这是最常用的做法。它的优点是全局最优性有保证因为是等价的单层MILP缺点是用户问题必须足够“良构”——凸、可导、强对偶成立一旦用户模型加入非凸约束就行不通。还有一种思路是迭代求解先猜一组电价热价把下层用户问题作为QP求解得到用户负荷再把用户负荷代入上层模型更新价格重复直到价格不再变化。这种思路更像数值迭代逼近博弈均衡能处理更复杂的下层模型但收敛性没有保证、速度也慢我一般只在KKT方案处理不了的时候才用。代码实现层面我建议先跑KKT方案确认模型无误后再尝试迭代方案做对比验证。5. 算例调试与避坑实录模型和数据都准备好了求解器一跑问题就来了。我把自己调试过程中遇到的最多的三类坑整理出来希望对你有帮助。5.1 坑一求解器报Infeasible但你的模型看起来“完全合理”这是新手最崩溃的报错。我遇到这个问题的根源大多数时候是互补松弛条件大M线性化写错了方向导致约束集合里出现了一个根本不可能同时满足的组合二进制变量z取0时lambda被强制为0但原始约束又必须被满足z取1时lambda被放开但原始约束被收紧。一旦逻辑写反模型会瞬间infessible。排查步骤可以这样先注释掉KKT条件和互补松弛约束只保留上层设备约束确认上层模型本身可行。再单独求解下层用户问题确认给定合理价格时用户问题的解存在。最后把KKT条件逐块加回上层模型每加一块求解一次定位到具体哪个约束导致冲突。我用这个三步骤法基本都能找到问题所在。另外Yalmip强烈建议加上check(Constraints)检查约束的一致性能快速发现变量维度不匹配等问题。5.2 坑二大M取值不当结果“看起来对但经不起推敲”大M的选取是个经验活。我见过M取100000、结果求解出的价格曲线高频振荡的情况也见过M取100、结果和真实最优解偏差很大的情况。M的核心逻辑是它必须大到足以“压过”乘子的任何合理取值但又不能大到让求解器的数值精度失效。Cplex默认的可行性容差是1e-6如果M和模型中的其他系数相差十几个数量级求解器在判断约束是否满足时可能出现误判。我总结的经验公式是M的取值参考乘子的理论最大值即目标函数中价格的上限乘以对偶变量的规模。价格上限已经约束在[ρ_min, ρ_max]区间时乘子的数量级通常不会超过价格上限本身。按这个量级取M再乘一个1.5~2倍的放大系数一般比较稳妥。调试时还可以这样做求解完后检查互补松弛条件是否真的满足——计算每个lambda和对应松弛量的乘积看它是否接近0。如果不满足说明M太大或太小需要调整。5.3 坑三热负荷模型简化过度导致定价结果偏离实际热力系统和电力系统的一大差异是“惯性”。我给热负荷建模时如果只做一个简单的热功率平衡不考虑建筑热惯性和热网管道的储热效应那么热价信号会非常敏感地跟随热负荷波动实际工程中根本没法执行。更合理的做法是给热负荷加入一个热惯性约束比如室内温度保持在舒适区间的约束。这样用户的热负荷就具备一定的可平移性热价高的时候可以稍微降低供热热价低的时候预供热模型更有物理含义。在代码里这体现为用户热负荷约束从简单的上下限变成温度动态约束T_in(t1) T_in(t) k1 × (H_L(t) - H_need(t)) - k2 × (T_in(t) - T_out(t))T_min ≤ T_in(t) ≤ T_max这个约束是非线性的温度乘以时间常数但可以线性化处理会增加一定的变量数量换来的是热负荷模型更真实。动态定价研究里我强烈建议采用带热惯性的用户模型否则热价的博弈结果没有说服力。5.4 结果分析主从博弈均衡解长什么样调通之后用一组典型参数跑出来的结果很有意思。以冬季典型日为例电价和热价会在早晚两个高峰时段明显抬升用户电负荷和热负荷在高峰时段被压低、在低谷时段回升蓄热罐在低价时段储热、高价时段放热CHP机组在电价高峰时加大电出力多余的余热存入蓄热罐而不是直接供给热负荷。这些现象本身就能验证模型的正确性如果跑出来的价格曲线完全不反映供需变化说明模型或数据有问题。我建议大家在结果分析时做一个对比实验场景一用固定电价热价比如按用户成本加成定价场景二用主从博弈动态定价对比两个场景下运营商的收益和用户的用能效用。你会看到动态定价在保障用户效用的同时运营商的收益有明显的提升空间——这就是博弈定价的价值所在。6. 做完这个项目后的几点想法这套主从博弈的框架不是万能的用的时候有几个边界条件值得注意。第一下层用户必须是价格接受者也就是单个用户不能影响价格如果园区里有一个能级很大的工业用户、能单独影响价格博弈就不是Stackelberg而是再复杂的寡头模型第二用户的效用函数参数怎么标定是个大问题实际场景里没有谁能告诉你a_e、b_e是多少通常要靠历史数据和价格弹性分析来拟合第三模型求解速度对算例规模敏感如果用户数量成百上千、又要考虑多时段热网动态MILP规模会爆炸式增长那时候就得考虑分解算法或者元启发式算法。从我的实际经验看主从博弈动态定价最合适的应用场景是区域级的综合能源系统比如一个园区、一个社区、一组商业楼宇用户数量从几个到几十个的规模。在这个规模下KKT转换加MILP求解的方法论是成熟、稳健、可复现的。如果往大了做就需要分布式求解的理论储备了。最后分享一个做这类研究的心得代码跑通只是第一步真正花时间的是“验证你的解确实是博弈均衡”——把求解出的价格代回下层用户问题看用户是否愿意维持原来的用能决策或者用迭代算法从不同的初始价格出发看最终是否收敛到同一个均衡点。这个验证步骤论文里写起来只是一句话但实际做起来能帮你发现模型里的很多隐藏问题。 基于主从博弈的电热综合能源系统动态定价与Matlab实现
电热综合能源系统的动态定价说穿了就是“谁来定价、怎么定价、定了之后用户听不听”的问题。我最早做这个项目的时候用传统的集中式优化把整个园区的电、热一起调度模型收敛得很快结果拿到实际场景里却没人买账——园区用户不是整体每个楼宇、每条产线都有自己的利益诉求运营商定的价格如果只让运营商赚钱用户就会用脚投票。后来我换成主从博弈框架把运营商设为上层领导者、用户设为下层跟随者问题才真正立起来。这篇文章就把我基于Matlab实现的一套“主从博弈电热综合能源系统动态定价能量管理”的完整思路和踩坑记录分享出来适合正在做综合能源优化调度、或者想用博弈论方法做定价研究的同学参考。1. 为什么电热综合能源定价会变成一个博弈问题1.1 运营商和用户的利益不一致才是定价问题的本质很多人在第一次接触“综合能源系统优化”时第一反应是建一个大目标函数把购电成本、购气成本、设备运维成本、用户满意度全部加权求和然后一次性求最优。这种做法的隐含假设是“整个系统是一个利益共同体”——但现实完全不是这样。综合能源系统中至少存在两类主体一类是综合能源服务商或者叫能源运营商它拥有CHP机组、电锅炉、储能设备负责从上级电网买电、从气网购气再把这些能源加工成电和热出售给用户另一类是终端用户它们关心的是用能成本最小化会随着电价、热价的变化调整自己的用电用热行为。这两类主体的目标函数天然冲突运营商想让售能收入最大化用户想让购能成本最小化。如果运营商盲目抬高价格用户会削减用能运营商的收入反而可能下降如果运营商压低价格用户用能增加但单位利润变薄。这是一个典型的利益博弈关系不是单个优化问题能刻画清楚的。1.2 主从博弈为什么天然适配“先定价、后响应”的模式主从博弈Stackelberg Game的核心思想是博弈双方的地位不对等一方先行动领导者另一方观察到领导者的行动后再做决策跟随者。领导者做决策时必须把跟随者的最优反应函数纳入自己的优化模型中。放到电热综合能源系统里这个结构非常自然领导者上层综合能源服务商。先制定下一个调度周期内比如24小时的电价、热价它的目标是净收益最大化扣除购能成本和设备运行成本。跟随者下层用户。给定电价和热价后用户调整自己的购电、购热计划目标是用能成本最小化同时满足自身的用能舒适度约束。这种“先定价、后响应”的时序关系和现实中综合能源服务商发布分时电价、用户根据价格调整空调温度设定或电动车充电时段的场景完全一致。1.3 和集中式优化比主从博弈模型带来了什么新信息集中式优化得到的结果是一个系统层面的“全局最优”但它回答不了“这个最优是怎么通过市场机制实现的”这个问题。主从博弈模型则给出了一个价格信号形成的机制上层先定价格下层做出最优响应上层再根据响应调整价格直到达到均衡。最终得到的电价、热价曲线是博弈均衡解而不是单纯的最优解。举个例子集中式优化可能得到“在某个时段大量用电”的调度方案但它不会告诉你电价应该是每千瓦时多少——电价是输入给定的。主从博弈模型却把“价格”本身作为决策变量这正好对应了电力市场改革后售电公司自主定价的现实需求。当然主从博弈也有代价模型是非凸的双层优化问题数学上很难直接求解需要做KKT条件转换、线性化等一系列处理。这也是这篇文章后面重点展开的部分。2. 电热耦合系统建模先把物理图景搭起来要写代码先得把数学模型写清楚。主从博弈不是凭空套壳下层用户能优化出什么东西取决于上层运营商的系统里有哪些设备、怎么运行。我先给出一套标准的电热综合能源系统拓扑和对应的数学模型。2.1 一个典型的电热综合能源系统长什么样我做项目时一般假设这样一个园区级系统调度周期为24小时时间尺度为1小时能源输入上级电网购电、天然气网购气。能量转换设备CHP机组热电联产同时产电产热、燃气锅炉或电锅炉。储能设备蓄电池、蓄热罐。负荷侧电负荷、热负荷由多个用户楼宇聚合而成。系统的能量流可以简化描述为天然气进入CHP机组后同时输出电功率和热功率燃气锅炉补充供热电网购电和CHP发电一起供应电负荷多余的可以充入蓄电池热负荷由CHP余热、锅炉供热和蓄热罐放热共同满足。热网部分可以做一个合理简化不细究管网拓扑和流体动态把热网看作一个传输环节主要考虑热功率平衡和蓄热罐的储放热约束。对于定价研究来说这种精度的模型已经够用细节热网模型会让博弈问题复杂到很难求解。2.2 设备模型CHP机组、电锅炉、储能、新能源的数学描述下面给出我常用的设备模型这些方程会直接进入上层优化问题。CHP机组CHP机组的热电耦合特性可以用“热电比”来刻画。最简单也最常用的模型是固定热电比模型电功率输出P_CHP(t)热功率输出H_CHP(t) r_CHP × P_CHP(t)其中r_CHP是热电比常数。更精细的做法是用一个可行域来描述CHP的运行范围这时H_CHP(t)和P_CHP(t)之间存在线性不等式约束。我在项目中用的是后者因为固定热电比的模型太理想无法反映机组在部分负荷下的调节能力。约束如下P_CHP(t) 在最小技术出力P_CHP_min和最大出力P_CHP_max之间。热出力与电出力满足上、下爬坡约束限制机组在相邻时段的出力变化量。燃料成本可以用二次函数近似但为了后续线性化处理会分段线性化。电锅炉/燃气锅炉电锅炉H_EB(t) η_EB × P_EB(t)输入电功率P_EB(t)输出热功率H_EB(t)。效率η_EB一般取0.9左右。燃气锅炉H_GB(t) η_GB × F_GB(t)输入天然气功率F_GB(t)输出热功率H_GB(t)。效率η_GB一般取0.85~0.95。蓄电池用荷电状态SOC来描述SOC(t1) SOC(t) η_ch × P_ch(t) × Δt / E_cap - P_dis(t) × Δt / (η_dis × E_cap)其中P_ch、P_dis是充放电功率η_ch、η_dis是充放电效率E_cap是额定容量。充放电功率都有上下限SOC也在[0.1, 0.9]之间且调度周期结束时SOC要回到初始值这是为了保证储能运行的可持续性。蓄热罐蓄热罐的模型和蓄电池非常像用储热量S(t)来描述S(t1) S(t) H_ch_st(t) × Δt - H_dis_st(t) × Δt蓄热罐的提升在于“热惰性”——热负荷高峰时可以通过放热支撑热负荷低谷时吸收多余热量。在动态定价问题中蓄热罐的配置会显著影响运营商的定价策略。2.3 用户侧建模价格弹性、热负荷特性、需求响应约束下层用户不是简单的“固定负荷曲线”而是会对价格做出响应。建立下层模型时我需要描述用户如何根据电价和热价调整用能。我采用一个常见的需求响应模型用户的总用能效用函数最大化。比如用户用电和用热的效用可以用二次效用函数描述U_e(P_L(t)) a_e × P_L(t) - 0.5 × b_e × P_L(t)^2U_h(H_L(t)) a_h × H_L(t) - 0.5 × b_h × H_L(t)^2那么用户的决策问题就是给定电价ρ_e(t)和热价ρ_h(t)选择用电功率P_L(t)、用热功率H_L(t)最大化总效用减去购能成本max Σ_t [U_e(P_L(t)) - ρ_e(t) × P_L(t) U_h(H_L(t)) - ρ_h(t) × H_L(t)]s.t. P_L_min ≤ P_L(t) ≤ P_L_maxH_L_min ≤ H_L(t) ≤ H_L_max这里的上下限就体现了用户的刚性需求哪怕电价再高基本负荷也得用哪怕热价再低热负荷也不能超过散热器或地暖系统允许的上限。这个二次效用函数有一个非常好的性质用户的最优响应是价格的线性函数——电价升高用电量线性下降。这一点是后面通过KKT条件把下层问题变成上层约束的关键基础。2.4 上层目标函数和价格约束上层运营商的优化目标是调度周期内净收益最大化max Σ_t [ρ_e(t) × P_L(t) ρ_h(t) × H_L(t) - C_buy(t) - C_fuel(t) - C_OM(t)]其中ρ_e(t) × P_L(t) 是售电收入ρ_h(t) × H_L(t) 是售热收入C_buy(t) 是向上级电网购电的成本按分时电价结算C_fuel(t) 是天然气购气成本C_OM(t) 是设备运维成本价格不能乱定要加约束ρ_e_min ≤ ρ_e(t) ≤ ρ_e_maxρ_h_min ≤ ρ_h(t) ≤ ρ_h_max这个约束代表的是监管要求或市场规则运营商不能无限抬价。做上层模型时如果不加价格上限优化算法很容易给出一个极高的价格你一看结果就知道这个模型跑飞了——下层用户的负荷会被压到下限整个系统失去意义。3. 双层模型的数学转换KKT条件与强对偶怎么用模型建好之后摆在面前的是一个双层优化问题上层运营商定价格最大化收益下层用户根据价格最小化用能成本等价于最大化效用减去成本这个双层问题不能直接丢给求解器。求解器能处理的要么是线性规划、要么是混合整数线性规划顶多是非线性规划但双层结构它不认识。标准的处理方法是把下层问题用KKT条件替换把双层问题转换成单层问题。3.1 先把双层问题写成标准形式为了推导方便我把下层用户问题简化表示为一个凸优化问题max_{x} f(x, ρ)s.t. g(x) ≤ 0其中ρ是上层给定的价格向量x是用户的用能向量f是用户的效用函数g是不等式约束集合。用户问题是凸的二次效用函数的Hessian矩阵是负定的约束是线性的所以KKT条件是用户问题的充要条件。这意味着把KKT条件加上就等价于把用户问题“嵌入”了上层模型。3.2 下层问题的KKT条件推导以下层用户的功率平衡约束为例。状态变量是P_L(t)和H_L(t)约束是上下限。我对每个不等式约束引入拉格朗日乘子λ_min、λ_max写出拉格朗日函数再对决策变量求偏导。用户问题的KKT条件包括三部分平稳性条件拉格朗日函数对每个决策变量的偏导数为0。对电功率P_L(t)∂L/∂P_L(t) a_e - b_e × P_L(t) - ρ_e(t) λ_e_min(t) - λ_e_max(t) 0对热功率H_L(t)同理∂L/∂H_L(t) a_h - b_h × H_L(t) - ρ_h(t) λ_h_min(t) - λ_h_max(t) 0原始可行性原问题的约束必须成立也就是P_L_min ≤ P_L(t) ≤ P_L_max、H_L_min ≤ H_L(t) ≤ H_L_max。互补松弛条件每个不等式约束和它对应的乘子不能同时“松弛”。写成数学形式就是λ_e_min(t) × (P_L_min - P_L(t)) 0λ_e_max(t) × (P_L(t) - P_L_max) 0对热负荷约束同理。乘子和约束的松弛量都是非负的。3.3 互补松弛的大M线性化互补松弛条件是非线性的两个变量相乘等于0而且包含整数逻辑求解器没法直接处理。行业里的标准做法是引入二进制变量用大M法线性化。比如对λ_e_min(t) × (P_L_min - P_L(t)) 0我引入二进制变量z1(t) ∈ {0,1}写成两个约束λ_e_min(t) ≤ M × z1(t)P_L_min - P_L(t) ≤ M × (1 - z1(t))这里M是一个足够大的正数大M。逻辑是如果z1(t)1则第一个约束给乘子上限第二个约束强制P_L(t)P_L_min即约束紧满足互补松弛如果z1(t)0则强制λ_e_min(t)0第二个约束被松弛到恒成立。每个互补松弛条件需要引入一个二进制变量。对于24小时的调度问题引入的二进制变量数量是可控的每类约束24个总共不过几十个到上百个二进制变量Cplex/Gurobi处理这个规模的混合整数规划毫无压力。3.4 非线性目标项的强对偶消去把KKT条件替换进上层模型后上层目标函数里仍然有ρ_e(t) × P_L(t)和ρ_h(t) × H_L(t)这类乘积项——价格变量和负荷变量相乘这不是线性的。怎么处理用强对偶定理。因为下层用户问题是凸优化强对偶成立下层问题的原目标最优值等于对偶问题的最优值。具体来说用户问题的目标函数是Σ_t [a_e × P_L(t) - 0.5 × b_e × P_L(t)^2 - ρ_e(t) × P_L(t) a_h × H_L(t) - 0.5 × b_h × H_L(t)^2 - ρ_h(t) × H_L(t)]它的对偶函数可以通过拉格朗日函数对决策变量求极大推导出来。利用强对偶我可以把上层目标函数中非线性项ρ_e(t) × P_L(t)替换成用户效用函数加上对偶函数的形式从而实现对目标函数的线性化。具体推导过程比较长这里直接给结论在引入KKT条件后上层目标函数中每个时段的收入和用户负荷的乘积项可以用用户效用项与拉格朗日乘子的内积等价替换。这样整个模型变成一个混合整数线性规划MILP可以直接用MatlabYalmipCplex求解。这一步是整个转换中最容易出错的地方也是很多论文里“推导略”的部分。强烈建议自己动手推一遍不要直接抄结果。4. MatlabYalmip代码实现从数学模型到可运行程序模型转换完之后编码实现就是体力活了。我用的工具链是Matlab Yalmip Cplex。下面把程序框架和关键实现细节列出来这套代码跑通之后你就可以在它的基础上改设备参数、改价格约束、换用户效用函数。4.1 为什么选YalmipCplex而不是手写求解器如果你刚接触优化建模我强烈建议先别自己写求解算法。主从博弈问题经过KKT转换后是一个MILP手写单纯形法或分支定界法的工作量完全不可接受。Yalmip是Matlab下面的一个建模工具箱它把优化建模抽象成“声明变量、写约束、调求解器”三步极大降低出bug的概率。求解器的选择上Cplex和Gurobi都是MILP领域的标配。我习惯用Cplex主要原因是它在学术界的license申请方便、和Yalmip的接口稳定。如果你有Gurobi的license完全可以直接换成gurobi代码不需要改太多。4.2 程序整体框架和模块划分整个程序我分成四个模块各司其职参数初始化模块定义24小时的电网分时购电价、天然气价格、CHP机组参数、储能参数、用户效用函数参数、价格上下限。变量声明模块定义上层的价格变量、设备出力变量、储能SOC变量以及下层用户的负荷变量和经KKT转换引入的拉格朗日乘子变量、二进制变量。约束组装模块把设备运行约束、用户KKT条件约束、互补松弛线性化约束逐一写入Yalmip的约束集合。模型求解与结果输出模块调用Cplex求解MILP解出后把电价、热价、用户负荷、各设备出力曲线保存成结构体或表格方便后续画图和对比。下面是一个关键代码片段用来创建上层运营商的决策变量% 时间尺度 N 24; % 上层变量 rho_e sdpvar(1, N, full); % 电价决策变量 rho_h sdpvar(1, N, full); % 热价决策变量 P_CHP sdpvar(1, N, full); % CHP电出力 H_CHP sdpvar(1, N, full); % CHP热出力 P_EB sdpvar(1, N, full); % 电锅炉电功率 SOC sdpvar(1, N1, full); % 蓄电池SOC含初始状态 % 下层变量跟随者的负荷决策 P_L sdpvar(1, N, full); % 电负荷 H_L sdpvar(1, N, full); % 热负荷 % KKT拉格朗日乘子 lambda_e_min sdpvar(1, N, full); % 电价相关下限乘子 lambda_e_max sdpvar(1, N, full); % 电价相关上限乘子 lambda_h_min sdpvar(1, N, full); lambda_h_max sdpvar(1, N, full); % 互补松弛线性化需要的二进制变量 z_e_min binvar(1, N, full); z_e_max binvar(1, N, full); z_h_min binvar(1, N, full); z_h_max binvar(1, N, full);4.3 设备约束和KKT约束怎么往Yalmip里写设备约束直接按前文的数学公式写注意SOC的索引错位问题。Yalmip支持用循环写约束也可以用矩阵形式批量写建议把每一天的约束写在一个for循环里逻辑清晰一点。下面是我写约束的示例这里只列出关键部分Constraints []; % 价格上下限 Constraints [Constraints, rho_e_min rho_e rho_e_max]; Constraints [Constraints, rho_h_min rho_h rho_h_max]; % CHP出力约束 for t 1:N Constraints [Constraints, P_CHP_min P_CHP(t) P_CHP_max]; Constraints [Constraints, H_CHP(t) r_CHP * P_CHP(t)]; end % 蓄电池SOC递推 Constraints [Constraints, SOC(1) SOC_0]; for t 1:N Constraints [Constraints, SOC(t1) SOC(t) ... - P_dis(t)*dt/eta_dis/E_cap eta_ch*P_ch(t)*dt/E_cap]; Constraints [Constraints, SOC_min SOC(t1) SOC_max]; end % 功率平衡约束 for t 1:N Constraints [Constraints, P_grid(t) P_CHP(t) P_dis(t) - P_ch(t) P_L(t)]; Constraints [Constraints, H_CHP(t) H_EB(t) H_dis_st(t) - H_ch_st(t) H_L(t)]; endKKT条件部分先把平稳性条件写进去% KKT平稳性条件 for t 1:N Constraints [Constraints, ... a_e - b_e * P_L(t) - rho_e(t) lambda_e_min(t) - lambda_e_max(t) 0]; Constraints [Constraints, ... a_h - b_h * H_L(t) - rho_h(t) lambda_h_min(t) - lambda_h_max(t) 0]; end然后是互补松弛条件的大M线性化M 1000; % 大M需要根据问题规模评估 for t 1:N % 电负荷下限互补条件 Constraints [Constraints, lambda_e_min(t) M * z_e_min(t)]; Constraints [Constraints, P_L_min - P_L(t) M * (1 - z_e_min(t))]; % 电负荷上限互补条件 Constraints [Constraints, lambda_e_max(t) M * z_e_max(t)]; Constraints [Constraints, P_L(t) - P_L_max M * (1 - z_e_max(t))]; % 热负荷同理 Constraints [Constraints, lambda_h_min(t) M * z_h_min(t)]; Constraints [Constraints, H_L_min - H_L(t) M * (1 - z_h_min(t))]; Constraints [Constraints, lambda_h_max(t) M * z_h_max(t)]; Constraints [Constraints, H_L(t) - H_L_max M * (1 - z_h_max(t))]; end大M的取值要特别小心。M太小会把最优解卡掉M太大又会导致数值病态、求解精度下降。我的做法是把用户负荷的上下限差值作为基准再乘上价格上限估算出一个合理的M一般取M 2 × max(价格上限) × max(负荷范围)。在Cplex里还可以设置数值精度的相关参数来缓解大M带来的问题。4.4 两种求解思路KKT合并一次性求解 vs 迭代优化上面给的是“把KKT条件写进上层模型、一次性求解MILP”的方案这是最常用的做法。它的优点是全局最优性有保证因为是等价的单层MILP缺点是用户问题必须足够“良构”——凸、可导、强对偶成立一旦用户模型加入非凸约束就行不通。还有一种思路是迭代求解先猜一组电价热价把下层用户问题作为QP求解得到用户负荷再把用户负荷代入上层模型更新价格重复直到价格不再变化。这种思路更像数值迭代逼近博弈均衡能处理更复杂的下层模型但收敛性没有保证、速度也慢我一般只在KKT方案处理不了的时候才用。代码实现层面我建议先跑KKT方案确认模型无误后再尝试迭代方案做对比验证。5. 算例调试与避坑实录模型和数据都准备好了求解器一跑问题就来了。我把自己调试过程中遇到的最多的三类坑整理出来希望对你有帮助。5.1 坑一求解器报Infeasible但你的模型看起来“完全合理”这是新手最崩溃的报错。我遇到这个问题的根源大多数时候是互补松弛条件大M线性化写错了方向导致约束集合里出现了一个根本不可能同时满足的组合二进制变量z取0时lambda被强制为0但原始约束又必须被满足z取1时lambda被放开但原始约束被收紧。一旦逻辑写反模型会瞬间infessible。排查步骤可以这样先注释掉KKT条件和互补松弛约束只保留上层设备约束确认上层模型本身可行。再单独求解下层用户问题确认给定合理价格时用户问题的解存在。最后把KKT条件逐块加回上层模型每加一块求解一次定位到具体哪个约束导致冲突。我用这个三步骤法基本都能找到问题所在。另外Yalmip强烈建议加上check(Constraints)检查约束的一致性能快速发现变量维度不匹配等问题。5.2 坑二大M取值不当结果“看起来对但经不起推敲”大M的选取是个经验活。我见过M取100000、结果求解出的价格曲线高频振荡的情况也见过M取100、结果和真实最优解偏差很大的情况。M的核心逻辑是它必须大到足以“压过”乘子的任何合理取值但又不能大到让求解器的数值精度失效。Cplex默认的可行性容差是1e-6如果M和模型中的其他系数相差十几个数量级求解器在判断约束是否满足时可能出现误判。我总结的经验公式是M的取值参考乘子的理论最大值即目标函数中价格的上限乘以对偶变量的规模。价格上限已经约束在[ρ_min, ρ_max]区间时乘子的数量级通常不会超过价格上限本身。按这个量级取M再乘一个1.5~2倍的放大系数一般比较稳妥。调试时还可以这样做求解完后检查互补松弛条件是否真的满足——计算每个lambda和对应松弛量的乘积看它是否接近0。如果不满足说明M太大或太小需要调整。5.3 坑三热负荷模型简化过度导致定价结果偏离实际热力系统和电力系统的一大差异是“惯性”。我给热负荷建模时如果只做一个简单的热功率平衡不考虑建筑热惯性和热网管道的储热效应那么热价信号会非常敏感地跟随热负荷波动实际工程中根本没法执行。更合理的做法是给热负荷加入一个热惯性约束比如室内温度保持在舒适区间的约束。这样用户的热负荷就具备一定的可平移性热价高的时候可以稍微降低供热热价低的时候预供热模型更有物理含义。在代码里这体现为用户热负荷约束从简单的上下限变成温度动态约束T_in(t1) T_in(t) k1 × (H_L(t) - H_need(t)) - k2 × (T_in(t) - T_out(t))T_min ≤ T_in(t) ≤ T_max这个约束是非线性的温度乘以时间常数但可以线性化处理会增加一定的变量数量换来的是热负荷模型更真实。动态定价研究里我强烈建议采用带热惯性的用户模型否则热价的博弈结果没有说服力。5.4 结果分析主从博弈均衡解长什么样调通之后用一组典型参数跑出来的结果很有意思。以冬季典型日为例电价和热价会在早晚两个高峰时段明显抬升用户电负荷和热负荷在高峰时段被压低、在低谷时段回升蓄热罐在低价时段储热、高价时段放热CHP机组在电价高峰时加大电出力多余的余热存入蓄热罐而不是直接供给热负荷。这些现象本身就能验证模型的正确性如果跑出来的价格曲线完全不反映供需变化说明模型或数据有问题。我建议大家在结果分析时做一个对比实验场景一用固定电价热价比如按用户成本加成定价场景二用主从博弈动态定价对比两个场景下运营商的收益和用户的用能效用。你会看到动态定价在保障用户效用的同时运营商的收益有明显的提升空间——这就是博弈定价的价值所在。6. 做完这个项目后的几点想法这套主从博弈的框架不是万能的用的时候有几个边界条件值得注意。第一下层用户必须是价格接受者也就是单个用户不能影响价格如果园区里有一个能级很大的工业用户、能单独影响价格博弈就不是Stackelberg而是再复杂的寡头模型第二用户的效用函数参数怎么标定是个大问题实际场景里没有谁能告诉你a_e、b_e是多少通常要靠历史数据和价格弹性分析来拟合第三模型求解速度对算例规模敏感如果用户数量成百上千、又要考虑多时段热网动态MILP规模会爆炸式增长那时候就得考虑分解算法或者元启发式算法。从我的实际经验看主从博弈动态定价最合适的应用场景是区域级的综合能源系统比如一个园区、一个社区、一组商业楼宇用户数量从几个到几十个的规模。在这个规模下KKT转换加MILP求解的方法论是成熟、稳健、可复现的。如果往大了做就需要分布式求解的理论储备了。最后分享一个做这类研究的心得代码跑通只是第一步真正花时间的是“验证你的解确实是博弈均衡”——把求解出的价格代回下层用户问题看用户是否愿意维持原来的用能决策或者用迭代算法从不同的初始价格出发看最终是否收敛到同一个均衡点。这个验证步骤论文里写起来只是一句话但实际做起来能帮你发现模型里的很多隐藏问题。