2026/9/9 20:55:23
FEATURED · 精选文章 
最近在搭一个计及需求响应和电能交互的多主体综合能源系统调度算例时又被主从博弈折腾了一轮。这套东西看着高深说白了就是解决一个现实矛盾综合能源系统运营商想让收益最大用户想让用能成本最低两边都不想被对方完全拿捏。集中式优化在这里行不通因为用户没有义务服从你制定的用能计划。所以要用主从博弈让运营商作为领导者先定能源价格用户作为跟随者再根据价格调整负荷最终收敛到一个双方都能接受的均衡点。这篇文章我打算用Matlab代码把整个思路落地把建模、KKT转化、Yalmip实现和调参踩坑完整梳理一遍。正在做综合能源系统优化调度、微电网/园区能源管理相关课题或者想用博弈论做论文创新点的同学可以参考。1. 多主体综合能源系统调度为何需要主从博弈1.1 从集中式优化到主从博弈一个实际场景先说一个我实际遇到的场景。一个园区综合能源系统里面有风电、光伏、燃气轮机、电锅炉和储能下面带着商业楼、工厂、居民区好几个用户主体。最开始我按传统思路做集中式优化采集所有用户负荷数据全局求解一个最小化总成本或最大化总收益的问题然后给每个用户下达“最优用电指令”。模型跑得很顺结果也漂亮但导师一句话问住我用户凭什么都听你的现实中的用户有独立利益商业楼希望电费最低但舒适度不能降工厂希望生产计划不被破坏居民用户则更在意电价波动和用能习惯。强行统一调度即使全局最优也落不了地。这就是多主体综合能源系统调度和传统单主体调度的本质区别信息私有、目标冲突、决策分散。主从博弈Stackelberg game恰好是刻画这种层级决策关系的工具。它天然有一个“领导者”和多个“跟随者”领导者先行动公布策略跟随者观察到领导者的策略后再做出最优反应领导者又会根据跟随者的反应调整策略直到没人愿意单方面改变。在综合能源系统里最自然的角色分配就是综合能源服务商Integrated Energy Service Operator简称IESO当领导者制定售电和售热价格多个负荷聚合商Load Aggregator简称LA当跟随者根据价格决定用能计划和需求响应量。最后得到的均衡解就是双方在利益博弈下都能接受的结果。1.2 需求响应与电能交互在博弈中的角色有了博弈框架需求响应和电能交互这两个词就不再是孤立的技术名词而是博弈中真正的“耦合变量”。需求响应本质上是给跟随者用户松绑让用户对价格信号产生弹性。如果用户负荷完全刚性运营商定价定多少用户都得买那博弈就没有意义。可一旦用户具备可转移负荷、可削减负荷、蓄热式电采暖等柔性资源价格升高时用户会主动削减或转移用电价格降低时用户会增加用电。运营商想高价卖电用户会“用脚投票”这就是供需双方的相互制约。电能交互体现在两个层面。一是运营商与上级电网之间的购售电交互运营商既可以从上级电网买电也可以在本地富余时向上级电网卖电这个交互功率本身就是运营商的重要决策变量。二是运营商与多个用户主体之间的电能交互用户向运营商购电在不同价格信号下改变购电量本质上是一种“主从之间”的功率交互。在我的模型里所有用户主体的购电需求汇聚到一起再由运营商统一与上级电网交互形成一个电能交换的枢纽。多主体的博弈结果最终都会反映到这些交互功率曲线里。1.3 本文采用的系统结构与博弈关系我搭建的系统结构比较典型也方便用Matlab实现。考虑一个IESO作为上层领导者服务N个负荷聚合商作为下层跟随者。IESO拥有分布式电源风电、光伏、燃气轮机、电锅炉、储能装置向下层用户出售电和热IESO同时与上级电网进行电能交互电价按分时电价执行。每个LA内部包含柔性负荷可转移、可削减和蓄热式电采暖负荷LA可以自主决定从IESO购买多少电、多少热以及参与需求响应的功率调整量。博弈流程如下IESO首先公布每个时段的售电价格和售热价格LA们看到价格后各自以用能成本最小化为目标优化购电计划、购热计划和需求响应量IESO根据LA的购能计划核算自己的收益和运行约束再调整价格重复这个过程直到价格和用能计划都不再变化。这个均衡点就是Stackelberg均衡。如果读者想把这个模型扩展到多个IESO之间通过联络线交互电能可以在上层增加相邻系统交互功率变量但求解复杂度会明显上升。下面我讲的求解方法本质上仍可以扩展到多领导者场景只是收敛性验证会更麻烦建议大家先把单领导者多跟随者跑通再说。2. 上层运营商优化模型定价、购售电与运行约束2.1 目标函数收益最大化上层IESO的目标是最大化整个调度周期内的净收益。这里有个容易混淆的地方很多教材喜欢把上层写成“成本最小化”但成本最小化看不到定价带来的收益弹性。在主从博弈里领导者要决定价格所以目标函数写成收益最大化更自然。上层目标函数包含这几部分售电收入所有LA在t时段的购电功率之和乘以IESO制定的售电价格售热收入所有LA在t时段的购热功率之和乘以售热价格向上级电网售电收入如果IESO本地发电量超过负荷需求可以把多余的电卖给上级电网向上级电网购电成本本地发电不足时从上级电网买电的费用燃气轮机的燃料成本、储能和新能源设备的运维成本对用户需求响应的补偿成本用户削减或转移负荷需要给予一定补偿否则用户没有动力参与。为了简洁上层目标可以写成下面的形式max Σ_t [ p_t^e * (Σ_i L_i,t^e) p_t^h * (Σ_i H_i,t^h) p_t^grid_out * P_t^out - p_t^grid_in * P_t^in - C_gt(P_gt,t) - C_om(...) - C_dr(...) ]其中 L_i,t^e 和 H_i,t^h 不是上层直接控制的变量而是下层LA对上层公布的 p_t^e、p_t^h 的最优响应。也就是说上层目标里的售电收入项 p_t^e * L_i,t^e 同时包含上层变量价格和下层变量购电量这是一个双线性项也是后面KKT转化时要重点处理的难点。2.2 约束条件功率平衡、交互容量与设备运行范围上层模型需要满足一系列运行约束这些约束在Matlab里用Yalmip写起来很直接但容易被忽略。我列一下我在代码里实际用到的约束电功率平衡约束燃气轮机出力 风电出力 光伏出力 储能放电功率 上级电网购电功率 LA总购电功率 上级电网售电功率 储能充电功率 电锅炉耗电功率。这里每个时段的等式平衡都是强约束求解不可行时第一个检查它。热功率平衡约束燃气轮机余热回收 电锅炉产热 蓄热装置放热 LA总购热功率 蓄热装置充热。热网动态可以简化成静态平衡如果考虑管道延时模型会变成复杂的时空耦合建议初学者先忽略。上级电网交互约束购电和售电不能同时为真需要引入0-1变量约束 P_t^in * P_t^out 0或者用一个大M约束限制。实际计算中如果不加互斥约束求解器经常给出既买又卖的荒唐结果。燃气轮机约束出力上下限、爬坡约束。注意爬坡约束要按分钟或小时折算。储能约束荷电状态SOC动态、充放电功率上下限、以及始末SOC相等。最后一个约束很多论文不提但实际上为了保证调度方案的可持续性必须加。电锅炉约束输入电功率与输出热功率之间有转换效率还有容量上下限。价格约束售电/售热价格有上下限这既来自监管也避免领导者随意定价导致下层问题无界。这些约束都相对成熟不再展开推导。重点是上层问题的可行域里价格变量 p_t^e 和 p_t^h 是决策变量而不是固定参数这会让整个问题变成带均衡约束的数学规划也就是MPEC。MPEC在Matlab里直接写给求解器是解不动的必须处理下层反应函数。2.3 为什么说“定价购售电”是上层博弈的核心抓手我在调试过程中最大的感受是上层模型里真正体现“博弈能力”的决策变量就是价格和与电网交互的功率。机组出力、储能充放电这些说到底是在响应价格信号和负荷需求属于执行层。如果只把储能和机组出力当作上层变量、把价格固定那模型就退化成了单主体的运行优化主从博弈就名存实亡。所以大家在复现论文时一定要确认清楚上层有没有把价格作为决策变量如果没有那只是“考虑需求响应的多主体优化”而不是严格意义上的主从博弈。审稿人或者答辩老师其实很看重这个逻辑。3. 下层用户需求响应模型柔性负荷与用能成本最小化3.1 可转移负荷、可削减负荷与蓄热式电采暖下层用户不是铁板一块我按负荷性质把每个LA内部的用电需求拆成三部分刚性负荷、可转移负荷、可削减负荷。刚性负荷没有调节空间就是必须消费的可转移负荷允许在时间轴上平移比如工厂的流水线、洗衣房的洗涤任务可削减负荷允许在高峰时段削减一部分比如空调温度上调一度、照明亮度下调一档代价是要给用户补偿。可转移负荷的建模要点转移前后总用电量保持守恒。比如某个用户一天需要10 MWh的转移负荷不能转移到最后总电量变成12 MWh。这个守恒约束用数学表达很简单Σ_t P_shift,i,t E_shift,i同时每个时段的转移功率要满足上下限并且只能转移到允许的时间窗内。很多初学者容易漏掉“同一个转移任务不能被拆散”的约束如果忽略离散档位纯连续模型也是常用简化。可削减负荷的建模要点削减量有上限通常表示成基线负荷的一定比例同时要限制削减次数或连续削减时长否则求解器会让用户在某个时段削减到零这在现实中是 unacceptable 的。我习惯加一个0-1变量表示用户在某时段是否参与削减然后令削减功率小于等于该变量乘以上限值。蓄热式电采暖是北方园区常见的柔性热负荷。它本质是一个小型储热装置SOC状态方程和电储能类似但输入是热功率并且热负荷有舒适度区间约束允许室温在某个范围内波动。如果你觉得同时建模三类负荷太复杂最低限度也要包含可转移和可削减两类这样需求响应才算真正“计及”了。3.2 用户用能成本最小化目标每个LA的决策目标是自己在调度周期内的净用能成本最小。净用能成本购电成本购热成本-需求响应补偿收益。注意需求响应补偿是IESO付给用户的所以对用户来说是负成本。min Σ_t [ p_t^e * (P_base,i,t^e P_shift,i,t P_cut,i,t) p_t^h * H_i,t λ_shift * P_shift,i,t^abs λ_cut * P_cut,i,t - ρ_t^dr * (P_shift,i,t P_cut,i,t) ]这里有个关键点当IESO公布价格 p_t^e、p_t^h 和补偿价格 ρ_t^dr 后用户的目标函数是线性的约束也是线性的所以下层是一个线性规划问题。线性规划非常重要因为它保证了后续KKT转化的正确性。如果下层是非线性凸规划KKT条件也能用但需要求解器支持二次约束如果是非凸问题KKT条件只是必要条件Stackelberg均衡的存在性就无法保证。所以我在构建下层模型时会刻意避免引入非线性项比如避免用“目标函数里含二次成本”这类导致下层非线性的设计。3.3 需求响应约束的凸性分析与KKT可行性为什么我反复强调下层要是线性规划因为主从博弈求解有一种非常漂亮的方式用下层问题的KKT条件替代下层优化问题。KKT条件成立的前提是下层问题满足约束规格比如Slater条件并且是凸的。线性规划天然满足。如果下层模型里出现了双线性项比如用户购电量乘以用户自己的某个决策变量那就把凸性破坏了后面所有KKT推导都失去依据。在Matlab里验证下层问题凸性很简单把所有约束打出来看一遍如果有变量相乘、绝对值、分段函数等就需要先线性化。我遇到过一个案例用户侧电动汽车的充电费用包含了“充电量乘以充电效率”的变量相乘必须把效率归到系数里避免出现变量相乘。这种细节写论文时很多人不注意但是到了编程阶段就一定会暴露出来。KKT可行性的另一个注意点是下层变量的非负约束也要写进KKT的互补松弛条件里。Yalmip里定义变量时可以声明为nonnegative这样在手动写KKT时要记得为每个非负变量生成对应的互补条件。漏掉一个就会导致求解出的“均衡”实际上只是单方面最优下层并没有真正达到最优。4. 主从博弈求解KKT转化与线性化处理4.1 用KKT条件将双层问题化为单层主从博弈直接求解很难问题在于上层目标里包含下层的最优反应函数。好在当下层是连续线性规划时我们可以用它的KKT条件把“最优反应”完整地刻画成一组约束和互补条件。这样双层问题就变成单层MPEC问题。具体做法对每个LA的优化模型写出拉格朗日函数对每个变量求偏导得到驻点条件列原始可行性约束列对偶可行性约束再叠加互补松弛条件。这些条件一起构成了“用户最优反应”的数学描述。把这组KKT条件替代到上层问题中就得到一个单层的约束优化问题。这一步说说容易手推真的很费时间。以4个LA、24个时段为例每个LA有约10个决策变量、15个不等式约束KKT条件会多出来上百条互补约束。全部手推容易出错我曾经因为漏了一个对偶变量结果算出来的电价高得离谱排查了整整两天。4.2 互补松弛条件的线性化KKT条件里最麻烦的是互补松弛条件形式是 λ * g(x) 0其中 λ≥0 是对偶变量g(x)≥0 是原始不等式约束。非线性乘积没法直接交给MILP求解器。标准做法是用大M法引入二进制变量把互补条件拆成一堆线性逻辑约束。以条件 λ * g(x) 0 为例等价于当 λ 0 时g(x) 0当 g(x) 0 时λ 0。用大M法写g(x) ≤ M * z λ ≤ M * (1 - z) λ ≥ 0, z ∈ {0, 1}这里z是二进制变量M是一个足够大的正数。这样就把乘积关系变成了线性约束。需要注意的是一个不等式约束对应一组互补约束每个互补约束需要一个二进制变量。如果下层约束多二进制变量数量会涨得很厉害所以实际求解时M的值和初始点对求解速度影响巨大。4.3 强对偶定理消除双线性项上层目标函数里的 p_t^e * L_t^e 仍然是个双线性项。因为L_t^e是下层变量它和上层价格相乘属于两个不同层级变量相乘的项MILP求解器不接受。我的处理方法是利用强对偶定理当下层LP达到最优时它的目标函数值等于其对偶问题的目标函数值。于是可以把下层目标函数中用价格乘以购电量产生的线性项替换成对偶问题里对应的对偶变量乘以常数系数的项。这样原始目标中的双线性项就消失了换来的是对偶变量和常数相乘的线性项。具体推导需要根据下层模型的每个约束逐一配平网上很多论文里叫“强对偶转化”或“dual rewrite”。这一步是整个算法里最容易出错的环节。我推荐一个检查技巧在完成强对偶替换后把求解结果代入分别计算下层原问题和对偶问题的目标值看是否相等。如果差了很多说明替换公式有误。我在调试时加了一个校准函数专门做这个检查。经过上述三步——KKT条件、互补松弛线性化、强对偶消双线性——原双层问题就转化成了一个混合整数线性规划问题可以直接用Gurobi或Cplex求解。5. Matlab代码实现框架与核心片段5.1 主程序流程与数据准备Matlab环境建议R2020b以上配Yalmip工具箱和一个商业求解器。我个人用Gurobi 10Cplex也可以但Gurobi在MILP求解上通常速度更快。主程序流程我分成五步写入基础数据24个时段、机组参数、储能参数、上级电网分时电价、用户基线负荷、可转移负荷量、可削减比例、价格上下限、需求响应补偿价格。用Yalmip定义所有变量区分上层变量、下层变量、对偶变量、二进制变量。添加上层约束、下层KKT约束、互补松弛线性化约束和强对偶目标替换。调用求解器求解MILP。提取结果绘制电价曲线、负荷响应曲线、储能SOC和收益柱状图。数据准备部分我直接写脚本生成没有用Excel读取。理由很简单调参时改数据太频繁每次改Excel还要重新读表不如在Matlab脚本里用结构体统一管理。比如% 基础时段 T 24; % 上级电网分时购电价单位元/kWh grid_buy_price [0.35*ones(1,8), 0.75*ones(1,7), 1.2*ones(1,5), 0.75*ones(1,4)]; % 用户个数 N 4; % 可转移负荷总量kWh每个用户一个值 E_shift_total [1200, 2400, 800, 1500]; % 可削减比例上限 cut_ratio_max [0.15, 0.20, 0.10, 0.25];5.2 Yalmip建模核心代码定义变量阶段最需要注意的是价格变量是上层决策变量但在下层KKT条件中又作为固定参数出现。实际处理时价格变量仍然用sdpvar声明因为在KKT中它要参与拉格朗日函数的偏导。用Matlab的大致写法如下% 上层变量 pe sdpvar(1, T); % 售电价 ph sdpvar(1, T); % 售热价 Pgrid_in sdpvar(1, T); % 上级电网购电 Pgrid_out sdpvar(1, T); % 向上级电网售电 Pgt sdpvar(1, T); % 燃气轮机出力 Pes_c sdpvar(1, T); % 储能充电 Pes_d sdpvar(1, T); % 储能放电 SOC sdpvar(1, T1); % 储能SOC Heb sdpvar(1, T); % 电锅炉产热 % 下层变量 Pshift sdpvar(N, T); % 可转移负荷 Pcut sdpvar(N, T); % 可削减负荷 Hbuy sdpvar(N, T); % 用户购热 Heat_store sdpvar(N, T); % 蓄热式电采暖充放热 % 对偶变量 lambda_eq sdpvar(N, T); % 用于下层等式约束 mu_ineq sdpvar(N, T); % 用于下层不等式约束 ...写约束和KKT时可以封装成函数把对偶变量和原始变量传进去返回约束数组。这样主函数清爽排查问题时也方便。5.3 求解器配置与结果输出求解前用Yalmip的ops设置求解器ops sdpsettings(solver, gurobi, verbose, 2); ops.gurobi.MIPGap 0.001; ops.gurobi.TimeLimit 600; ops.gurobi.NumericFocus 3;NumericFocus 3这个选项非常重要因为经过KKT和大M转化后模型数值尺度差异很大不同数量级变量混在一起默认参数容易让求解器在单纯形迭代里出现数值困难。我一开始没设置这个Gurobi经常报“Numerical trouble”警告设成3之后明显稳定。结果输出部分我一般会把决策结果保存到结构体里然后用plot画四张图电价与负荷响应曲线、各类电源出力柱状堆叠图、储能SOC和交互功率、各用户成本占比。如果有对比场景就画柱状图对比收益和成本方便写报告。6. 算例验证与结果分析6.1 参数设置与场景设计为了让结果能说明问题我设计了三个对比场景场景A无需求响应、无电能交互用户负荷固定IESO只能向上级电网购电满足刚性负荷。场景B仅考虑需求响应用户具备可转移/可削减负荷但不允许IESO向上级电网售电。场景C同时考虑需求响应和电能交互即本文的完整主从博弈模型。基础参数如下24个时段4个LA上级电网购电价采用峰谷平三段峰段1.2元/kWh、平段0.75元/kWh、谷段0.35元/kWhIESO售电价上下限为0.3~1.5元/kWh售热价上下限为0.2~0.8元/kWh燃气轮机容量2 MW效率0.4天然气热值按9.7 kWh/m³计算储能容量1 MWh充放电效率0.95最大充放电功率0.25 MW电锅炉效率0.9容量1.5 MW用户可转移负荷比例约15%~25%可削减比例约10%~25%。6.2 博弈均衡结果分析以场景C为例求解结束后得到的售电价曲线很有意思高峰时段电价明显抬升但并不是直接顶到上限1.5元/kWh而是停留在1.32元/kWh附近。原因是如果电价定得太高用户的可转移负荷会大规模转移出高峰时段虽然单价高但销售量下降总收益反而变小如果定得太低高峰时段上级电网购电成本高IESO利润空间不够。均衡电价正好是两者博弈后的折中。低谷时段电价下探到0.38元/kWh吸引用户把可转移负荷转移到夜间既提高了用户侧的用能经济性也提高了IESO低谷时段风电的消纳能力。用户侧的需求响应行为也比较直观高峰时段4个LA的可削减负荷平均削减了基线负荷的18%左右可转移负荷中约70%转移到了22:00到次日6:00的谷段。蓄热式电采暖在谷段自动充热早高峰时段放热热负荷的峰谷差显著减小。IESO与上级电网的电能交互上购电功率主要出现在16:00到21:00的峰段而凌晨2:00到6:00之间由于燃气轮机和风电出力充裕IESO会选择向上级电网售电。这里的“交互”不是固定的而是根据价格信号动态调整的充分体现了“电能交互”在博弈中的作用。6.3 需求响应和电能交互带来的效益对比三个场景的优化结果如下表指标场景A无需求响应/无交互场景B仅需求响应场景C需求响应电能交互IESO总收益元256002970034150用户总用能成本元201001760015700系统峰时购电功率MW2.411.931.67系统总购电成本元12600114009800从表里能看出场景C相比场景AIESO收益提高了33%用户成本降低了22%峰时购电功率下降了30%。需求响应让用户把负荷从高峰挪到低谷电能交互让IESO能利用地区间的价差灵活买卖电两个机制叠加实现了双赢。有一点必须提醒这些数字是特定参数下的结果不代表所有场景下都有这么大的提升。如果把上级电网的峰谷价差调小电能交互的收益就会变小如果把用户可转移负荷比例调低需求响应的作用也会降低。所以论文里做敏感性分析是有必要的这能验证模型在不同环境下的适应性。7. 实际代码调参中的踩坑记录7.1 迭代法 vs KKT单层法为什么我放弃迭代一开始我图省事用“交替迭代”实现主从博弈上层先给定初始价格求解下层得到负荷响应再把负荷响应代回上层用启发式规则更新价格重复直到收敛。结果发现价格经常在上下限之间来回振荡怎么调步长都压不住而且更新规则没有理论保证审稿人一问“怎么证明收敛”我就没法回答。后来老老实实改用KKT单层转化。虽然推导和编码量大了不少但求解过程是确定性的Gurobi给出的是全局最优解在给定离散变量取值下不需要担心振荡问题。如果你的目标只是快速验证想法迭代法可以跑着看趋势但要做严谨的科研或者实际项目KKT转化是更稳妥的路线。7.2 大M值选取和数值病态互补松弛线性化里的M值选择是影响求解稳定性的最直接因素。M取得太小可能把本该可行的解错误截断M取得太大会让约束矩阵条件数变差求解器陷入数值困难。我踩过的最夸张一次M取了1e8结果Gurobi直接警告“Matrix is ill-conditioned”求解出的电价出现很多肉眼可见的不合理尖峰。我的经验是M值根据对应约束中变量的数量级来取通常取该约束可能达到的最大值的2~3倍。比如功率平衡约束右侧最大负荷接近3 MWM取10就够对偶变量的M值可以按对偶变量理论上限取比如1000。不要图省事全模型统一用一个很大的M。必要的话可以先解一个RELAXATION版本统计原始约束的最大mismatch再据此定M。7.3 不可行问题的系统排查经验加了KKT和互补条件后模型规模激增不可行几乎是必然遇到的情况。我排查不可行时有一套固定流程先去掉所有互补松弛约束求解一个松弛问题看看原始约束和目标函数是否可行。如果这时候还不可行说明问题出在功率平衡、设备容量这些基础约束上检查和上层、下层约束定义是否冲突。如果松弛问题可行再逐一加入互补约束定位是哪一组互补约束导致不可行。通常是大M值太小或者是某个变量忘记定义非负。另外Yalmip默认的求解信息里会给出“infeasible”提示但不会直接告诉你哪条约束坏了。我习惯用yalmip(clear)之后把约束数组拆分成小块逐一验证配合check(Constraint)查看残差。这个方法比较笨但确实能救命。7.4 关于Yalmip的kkt命令Yalmip自带一个kkt函数可以自动生成大规模优化问题的KKT条件。听起来很方便但实测有几个限制它要求问题必须是凸的而且对变量类型和约束形式有比较严格的要求遇到含二进制变量的上层问题时直接把kkt用在目标里会报错生成的KKT条件会比较冗余导致求解性能下降。另外kkt生成的互补条件也是非线性仍需自己处理。所以我的建议是模型简单时可以用kkt快速验证思路但正式项目还是手动推导更可控。手动推导的过程虽然痛苦但能让你对模型结构理解更深后面扩展模型时也不至于一脸懵。以上这些经验都是我从最初能跑通简单算例到最终实现完整主从博弈模型的过程中一步步攒下来的。代码不是一次写成的模型也不是一步到位的。如果你现在正在折腾双层优化、多主体博弈建议先把一个最小的两主体、两时段模型跑通再慢慢扩展到24时段、多用户。小模型上的KKT推导正确了大模型只是“复制粘贴改索引”的事但前提是你真的理解了每一行KKT条件是从哪条约束来的。
主从博弈下综合能源系统调度：需求响应与电能交互的Matlab实现
最近在搭一个计及需求响应和电能交互的多主体综合能源系统调度算例时又被主从博弈折腾了一轮。这套东西看着高深说白了就是解决一个现实矛盾综合能源系统运营商想让收益最大用户想让用能成本最低两边都不想被对方完全拿捏。集中式优化在这里行不通因为用户没有义务服从你制定的用能计划。所以要用主从博弈让运营商作为领导者先定能源价格用户作为跟随者再根据价格调整负荷最终收敛到一个双方都能接受的均衡点。这篇文章我打算用Matlab代码把整个思路落地把建模、KKT转化、Yalmip实现和调参踩坑完整梳理一遍。正在做综合能源系统优化调度、微电网/园区能源管理相关课题或者想用博弈论做论文创新点的同学可以参考。1. 多主体综合能源系统调度为何需要主从博弈1.1 从集中式优化到主从博弈一个实际场景先说一个我实际遇到的场景。一个园区综合能源系统里面有风电、光伏、燃气轮机、电锅炉和储能下面带着商业楼、工厂、居民区好几个用户主体。最开始我按传统思路做集中式优化采集所有用户负荷数据全局求解一个最小化总成本或最大化总收益的问题然后给每个用户下达“最优用电指令”。模型跑得很顺结果也漂亮但导师一句话问住我用户凭什么都听你的现实中的用户有独立利益商业楼希望电费最低但舒适度不能降工厂希望生产计划不被破坏居民用户则更在意电价波动和用能习惯。强行统一调度即使全局最优也落不了地。这就是多主体综合能源系统调度和传统单主体调度的本质区别信息私有、目标冲突、决策分散。主从博弈Stackelberg game恰好是刻画这种层级决策关系的工具。它天然有一个“领导者”和多个“跟随者”领导者先行动公布策略跟随者观察到领导者的策略后再做出最优反应领导者又会根据跟随者的反应调整策略直到没人愿意单方面改变。在综合能源系统里最自然的角色分配就是综合能源服务商Integrated Energy Service Operator简称IESO当领导者制定售电和售热价格多个负荷聚合商Load Aggregator简称LA当跟随者根据价格决定用能计划和需求响应量。最后得到的均衡解就是双方在利益博弈下都能接受的结果。1.2 需求响应与电能交互在博弈中的角色有了博弈框架需求响应和电能交互这两个词就不再是孤立的技术名词而是博弈中真正的“耦合变量”。需求响应本质上是给跟随者用户松绑让用户对价格信号产生弹性。如果用户负荷完全刚性运营商定价定多少用户都得买那博弈就没有意义。可一旦用户具备可转移负荷、可削减负荷、蓄热式电采暖等柔性资源价格升高时用户会主动削减或转移用电价格降低时用户会增加用电。运营商想高价卖电用户会“用脚投票”这就是供需双方的相互制约。电能交互体现在两个层面。一是运营商与上级电网之间的购售电交互运营商既可以从上级电网买电也可以在本地富余时向上级电网卖电这个交互功率本身就是运营商的重要决策变量。二是运营商与多个用户主体之间的电能交互用户向运营商购电在不同价格信号下改变购电量本质上是一种“主从之间”的功率交互。在我的模型里所有用户主体的购电需求汇聚到一起再由运营商统一与上级电网交互形成一个电能交换的枢纽。多主体的博弈结果最终都会反映到这些交互功率曲线里。1.3 本文采用的系统结构与博弈关系我搭建的系统结构比较典型也方便用Matlab实现。考虑一个IESO作为上层领导者服务N个负荷聚合商作为下层跟随者。IESO拥有分布式电源风电、光伏、燃气轮机、电锅炉、储能装置向下层用户出售电和热IESO同时与上级电网进行电能交互电价按分时电价执行。每个LA内部包含柔性负荷可转移、可削减和蓄热式电采暖负荷LA可以自主决定从IESO购买多少电、多少热以及参与需求响应的功率调整量。博弈流程如下IESO首先公布每个时段的售电价格和售热价格LA们看到价格后各自以用能成本最小化为目标优化购电计划、购热计划和需求响应量IESO根据LA的购能计划核算自己的收益和运行约束再调整价格重复这个过程直到价格和用能计划都不再变化。这个均衡点就是Stackelberg均衡。如果读者想把这个模型扩展到多个IESO之间通过联络线交互电能可以在上层增加相邻系统交互功率变量但求解复杂度会明显上升。下面我讲的求解方法本质上仍可以扩展到多领导者场景只是收敛性验证会更麻烦建议大家先把单领导者多跟随者跑通再说。2. 上层运营商优化模型定价、购售电与运行约束2.1 目标函数收益最大化上层IESO的目标是最大化整个调度周期内的净收益。这里有个容易混淆的地方很多教材喜欢把上层写成“成本最小化”但成本最小化看不到定价带来的收益弹性。在主从博弈里领导者要决定价格所以目标函数写成收益最大化更自然。上层目标函数包含这几部分售电收入所有LA在t时段的购电功率之和乘以IESO制定的售电价格售热收入所有LA在t时段的购热功率之和乘以售热价格向上级电网售电收入如果IESO本地发电量超过负荷需求可以把多余的电卖给上级电网向上级电网购电成本本地发电不足时从上级电网买电的费用燃气轮机的燃料成本、储能和新能源设备的运维成本对用户需求响应的补偿成本用户削减或转移负荷需要给予一定补偿否则用户没有动力参与。为了简洁上层目标可以写成下面的形式max Σ_t [ p_t^e * (Σ_i L_i,t^e) p_t^h * (Σ_i H_i,t^h) p_t^grid_out * P_t^out - p_t^grid_in * P_t^in - C_gt(P_gt,t) - C_om(...) - C_dr(...) ]其中 L_i,t^e 和 H_i,t^h 不是上层直接控制的变量而是下层LA对上层公布的 p_t^e、p_t^h 的最优响应。也就是说上层目标里的售电收入项 p_t^e * L_i,t^e 同时包含上层变量价格和下层变量购电量这是一个双线性项也是后面KKT转化时要重点处理的难点。2.2 约束条件功率平衡、交互容量与设备运行范围上层模型需要满足一系列运行约束这些约束在Matlab里用Yalmip写起来很直接但容易被忽略。我列一下我在代码里实际用到的约束电功率平衡约束燃气轮机出力 风电出力 光伏出力 储能放电功率 上级电网购电功率 LA总购电功率 上级电网售电功率 储能充电功率 电锅炉耗电功率。这里每个时段的等式平衡都是强约束求解不可行时第一个检查它。热功率平衡约束燃气轮机余热回收 电锅炉产热 蓄热装置放热 LA总购热功率 蓄热装置充热。热网动态可以简化成静态平衡如果考虑管道延时模型会变成复杂的时空耦合建议初学者先忽略。上级电网交互约束购电和售电不能同时为真需要引入0-1变量约束 P_t^in * P_t^out 0或者用一个大M约束限制。实际计算中如果不加互斥约束求解器经常给出既买又卖的荒唐结果。燃气轮机约束出力上下限、爬坡约束。注意爬坡约束要按分钟或小时折算。储能约束荷电状态SOC动态、充放电功率上下限、以及始末SOC相等。最后一个约束很多论文不提但实际上为了保证调度方案的可持续性必须加。电锅炉约束输入电功率与输出热功率之间有转换效率还有容量上下限。价格约束售电/售热价格有上下限这既来自监管也避免领导者随意定价导致下层问题无界。这些约束都相对成熟不再展开推导。重点是上层问题的可行域里价格变量 p_t^e 和 p_t^h 是决策变量而不是固定参数这会让整个问题变成带均衡约束的数学规划也就是MPEC。MPEC在Matlab里直接写给求解器是解不动的必须处理下层反应函数。2.3 为什么说“定价购售电”是上层博弈的核心抓手我在调试过程中最大的感受是上层模型里真正体现“博弈能力”的决策变量就是价格和与电网交互的功率。机组出力、储能充放电这些说到底是在响应价格信号和负荷需求属于执行层。如果只把储能和机组出力当作上层变量、把价格固定那模型就退化成了单主体的运行优化主从博弈就名存实亡。所以大家在复现论文时一定要确认清楚上层有没有把价格作为决策变量如果没有那只是“考虑需求响应的多主体优化”而不是严格意义上的主从博弈。审稿人或者答辩老师其实很看重这个逻辑。3. 下层用户需求响应模型柔性负荷与用能成本最小化3.1 可转移负荷、可削减负荷与蓄热式电采暖下层用户不是铁板一块我按负荷性质把每个LA内部的用电需求拆成三部分刚性负荷、可转移负荷、可削减负荷。刚性负荷没有调节空间就是必须消费的可转移负荷允许在时间轴上平移比如工厂的流水线、洗衣房的洗涤任务可削减负荷允许在高峰时段削减一部分比如空调温度上调一度、照明亮度下调一档代价是要给用户补偿。可转移负荷的建模要点转移前后总用电量保持守恒。比如某个用户一天需要10 MWh的转移负荷不能转移到最后总电量变成12 MWh。这个守恒约束用数学表达很简单Σ_t P_shift,i,t E_shift,i同时每个时段的转移功率要满足上下限并且只能转移到允许的时间窗内。很多初学者容易漏掉“同一个转移任务不能被拆散”的约束如果忽略离散档位纯连续模型也是常用简化。可削减负荷的建模要点削减量有上限通常表示成基线负荷的一定比例同时要限制削减次数或连续削减时长否则求解器会让用户在某个时段削减到零这在现实中是 unacceptable 的。我习惯加一个0-1变量表示用户在某时段是否参与削减然后令削减功率小于等于该变量乘以上限值。蓄热式电采暖是北方园区常见的柔性热负荷。它本质是一个小型储热装置SOC状态方程和电储能类似但输入是热功率并且热负荷有舒适度区间约束允许室温在某个范围内波动。如果你觉得同时建模三类负荷太复杂最低限度也要包含可转移和可削减两类这样需求响应才算真正“计及”了。3.2 用户用能成本最小化目标每个LA的决策目标是自己在调度周期内的净用能成本最小。净用能成本购电成本购热成本-需求响应补偿收益。注意需求响应补偿是IESO付给用户的所以对用户来说是负成本。min Σ_t [ p_t^e * (P_base,i,t^e P_shift,i,t P_cut,i,t) p_t^h * H_i,t λ_shift * P_shift,i,t^abs λ_cut * P_cut,i,t - ρ_t^dr * (P_shift,i,t P_cut,i,t) ]这里有个关键点当IESO公布价格 p_t^e、p_t^h 和补偿价格 ρ_t^dr 后用户的目标函数是线性的约束也是线性的所以下层是一个线性规划问题。线性规划非常重要因为它保证了后续KKT转化的正确性。如果下层是非线性凸规划KKT条件也能用但需要求解器支持二次约束如果是非凸问题KKT条件只是必要条件Stackelberg均衡的存在性就无法保证。所以我在构建下层模型时会刻意避免引入非线性项比如避免用“目标函数里含二次成本”这类导致下层非线性的设计。3.3 需求响应约束的凸性分析与KKT可行性为什么我反复强调下层要是线性规划因为主从博弈求解有一种非常漂亮的方式用下层问题的KKT条件替代下层优化问题。KKT条件成立的前提是下层问题满足约束规格比如Slater条件并且是凸的。线性规划天然满足。如果下层模型里出现了双线性项比如用户购电量乘以用户自己的某个决策变量那就把凸性破坏了后面所有KKT推导都失去依据。在Matlab里验证下层问题凸性很简单把所有约束打出来看一遍如果有变量相乘、绝对值、分段函数等就需要先线性化。我遇到过一个案例用户侧电动汽车的充电费用包含了“充电量乘以充电效率”的变量相乘必须把效率归到系数里避免出现变量相乘。这种细节写论文时很多人不注意但是到了编程阶段就一定会暴露出来。KKT可行性的另一个注意点是下层变量的非负约束也要写进KKT的互补松弛条件里。Yalmip里定义变量时可以声明为nonnegative这样在手动写KKT时要记得为每个非负变量生成对应的互补条件。漏掉一个就会导致求解出的“均衡”实际上只是单方面最优下层并没有真正达到最优。4. 主从博弈求解KKT转化与线性化处理4.1 用KKT条件将双层问题化为单层主从博弈直接求解很难问题在于上层目标里包含下层的最优反应函数。好在当下层是连续线性规划时我们可以用它的KKT条件把“最优反应”完整地刻画成一组约束和互补条件。这样双层问题就变成单层MPEC问题。具体做法对每个LA的优化模型写出拉格朗日函数对每个变量求偏导得到驻点条件列原始可行性约束列对偶可行性约束再叠加互补松弛条件。这些条件一起构成了“用户最优反应”的数学描述。把这组KKT条件替代到上层问题中就得到一个单层的约束优化问题。这一步说说容易手推真的很费时间。以4个LA、24个时段为例每个LA有约10个决策变量、15个不等式约束KKT条件会多出来上百条互补约束。全部手推容易出错我曾经因为漏了一个对偶变量结果算出来的电价高得离谱排查了整整两天。4.2 互补松弛条件的线性化KKT条件里最麻烦的是互补松弛条件形式是 λ * g(x) 0其中 λ≥0 是对偶变量g(x)≥0 是原始不等式约束。非线性乘积没法直接交给MILP求解器。标准做法是用大M法引入二进制变量把互补条件拆成一堆线性逻辑约束。以条件 λ * g(x) 0 为例等价于当 λ 0 时g(x) 0当 g(x) 0 时λ 0。用大M法写g(x) ≤ M * z λ ≤ M * (1 - z) λ ≥ 0, z ∈ {0, 1}这里z是二进制变量M是一个足够大的正数。这样就把乘积关系变成了线性约束。需要注意的是一个不等式约束对应一组互补约束每个互补约束需要一个二进制变量。如果下层约束多二进制变量数量会涨得很厉害所以实际求解时M的值和初始点对求解速度影响巨大。4.3 强对偶定理消除双线性项上层目标函数里的 p_t^e * L_t^e 仍然是个双线性项。因为L_t^e是下层变量它和上层价格相乘属于两个不同层级变量相乘的项MILP求解器不接受。我的处理方法是利用强对偶定理当下层LP达到最优时它的目标函数值等于其对偶问题的目标函数值。于是可以把下层目标函数中用价格乘以购电量产生的线性项替换成对偶问题里对应的对偶变量乘以常数系数的项。这样原始目标中的双线性项就消失了换来的是对偶变量和常数相乘的线性项。具体推导需要根据下层模型的每个约束逐一配平网上很多论文里叫“强对偶转化”或“dual rewrite”。这一步是整个算法里最容易出错的环节。我推荐一个检查技巧在完成强对偶替换后把求解结果代入分别计算下层原问题和对偶问题的目标值看是否相等。如果差了很多说明替换公式有误。我在调试时加了一个校准函数专门做这个检查。经过上述三步——KKT条件、互补松弛线性化、强对偶消双线性——原双层问题就转化成了一个混合整数线性规划问题可以直接用Gurobi或Cplex求解。5. Matlab代码实现框架与核心片段5.1 主程序流程与数据准备Matlab环境建议R2020b以上配Yalmip工具箱和一个商业求解器。我个人用Gurobi 10Cplex也可以但Gurobi在MILP求解上通常速度更快。主程序流程我分成五步写入基础数据24个时段、机组参数、储能参数、上级电网分时电价、用户基线负荷、可转移负荷量、可削减比例、价格上下限、需求响应补偿价格。用Yalmip定义所有变量区分上层变量、下层变量、对偶变量、二进制变量。添加上层约束、下层KKT约束、互补松弛线性化约束和强对偶目标替换。调用求解器求解MILP。提取结果绘制电价曲线、负荷响应曲线、储能SOC和收益柱状图。数据准备部分我直接写脚本生成没有用Excel读取。理由很简单调参时改数据太频繁每次改Excel还要重新读表不如在Matlab脚本里用结构体统一管理。比如% 基础时段 T 24; % 上级电网分时购电价单位元/kWh grid_buy_price [0.35*ones(1,8), 0.75*ones(1,7), 1.2*ones(1,5), 0.75*ones(1,4)]; % 用户个数 N 4; % 可转移负荷总量kWh每个用户一个值 E_shift_total [1200, 2400, 800, 1500]; % 可削减比例上限 cut_ratio_max [0.15, 0.20, 0.10, 0.25];5.2 Yalmip建模核心代码定义变量阶段最需要注意的是价格变量是上层决策变量但在下层KKT条件中又作为固定参数出现。实际处理时价格变量仍然用sdpvar声明因为在KKT中它要参与拉格朗日函数的偏导。用Matlab的大致写法如下% 上层变量 pe sdpvar(1, T); % 售电价 ph sdpvar(1, T); % 售热价 Pgrid_in sdpvar(1, T); % 上级电网购电 Pgrid_out sdpvar(1, T); % 向上级电网售电 Pgt sdpvar(1, T); % 燃气轮机出力 Pes_c sdpvar(1, T); % 储能充电 Pes_d sdpvar(1, T); % 储能放电 SOC sdpvar(1, T1); % 储能SOC Heb sdpvar(1, T); % 电锅炉产热 % 下层变量 Pshift sdpvar(N, T); % 可转移负荷 Pcut sdpvar(N, T); % 可削减负荷 Hbuy sdpvar(N, T); % 用户购热 Heat_store sdpvar(N, T); % 蓄热式电采暖充放热 % 对偶变量 lambda_eq sdpvar(N, T); % 用于下层等式约束 mu_ineq sdpvar(N, T); % 用于下层不等式约束 ...写约束和KKT时可以封装成函数把对偶变量和原始变量传进去返回约束数组。这样主函数清爽排查问题时也方便。5.3 求解器配置与结果输出求解前用Yalmip的ops设置求解器ops sdpsettings(solver, gurobi, verbose, 2); ops.gurobi.MIPGap 0.001; ops.gurobi.TimeLimit 600; ops.gurobi.NumericFocus 3;NumericFocus 3这个选项非常重要因为经过KKT和大M转化后模型数值尺度差异很大不同数量级变量混在一起默认参数容易让求解器在单纯形迭代里出现数值困难。我一开始没设置这个Gurobi经常报“Numerical trouble”警告设成3之后明显稳定。结果输出部分我一般会把决策结果保存到结构体里然后用plot画四张图电价与负荷响应曲线、各类电源出力柱状堆叠图、储能SOC和交互功率、各用户成本占比。如果有对比场景就画柱状图对比收益和成本方便写报告。6. 算例验证与结果分析6.1 参数设置与场景设计为了让结果能说明问题我设计了三个对比场景场景A无需求响应、无电能交互用户负荷固定IESO只能向上级电网购电满足刚性负荷。场景B仅考虑需求响应用户具备可转移/可削减负荷但不允许IESO向上级电网售电。场景C同时考虑需求响应和电能交互即本文的完整主从博弈模型。基础参数如下24个时段4个LA上级电网购电价采用峰谷平三段峰段1.2元/kWh、平段0.75元/kWh、谷段0.35元/kWhIESO售电价上下限为0.3~1.5元/kWh售热价上下限为0.2~0.8元/kWh燃气轮机容量2 MW效率0.4天然气热值按9.7 kWh/m³计算储能容量1 MWh充放电效率0.95最大充放电功率0.25 MW电锅炉效率0.9容量1.5 MW用户可转移负荷比例约15%~25%可削减比例约10%~25%。6.2 博弈均衡结果分析以场景C为例求解结束后得到的售电价曲线很有意思高峰时段电价明显抬升但并不是直接顶到上限1.5元/kWh而是停留在1.32元/kWh附近。原因是如果电价定得太高用户的可转移负荷会大规模转移出高峰时段虽然单价高但销售量下降总收益反而变小如果定得太低高峰时段上级电网购电成本高IESO利润空间不够。均衡电价正好是两者博弈后的折中。低谷时段电价下探到0.38元/kWh吸引用户把可转移负荷转移到夜间既提高了用户侧的用能经济性也提高了IESO低谷时段风电的消纳能力。用户侧的需求响应行为也比较直观高峰时段4个LA的可削减负荷平均削减了基线负荷的18%左右可转移负荷中约70%转移到了22:00到次日6:00的谷段。蓄热式电采暖在谷段自动充热早高峰时段放热热负荷的峰谷差显著减小。IESO与上级电网的电能交互上购电功率主要出现在16:00到21:00的峰段而凌晨2:00到6:00之间由于燃气轮机和风电出力充裕IESO会选择向上级电网售电。这里的“交互”不是固定的而是根据价格信号动态调整的充分体现了“电能交互”在博弈中的作用。6.3 需求响应和电能交互带来的效益对比三个场景的优化结果如下表指标场景A无需求响应/无交互场景B仅需求响应场景C需求响应电能交互IESO总收益元256002970034150用户总用能成本元201001760015700系统峰时购电功率MW2.411.931.67系统总购电成本元12600114009800从表里能看出场景C相比场景AIESO收益提高了33%用户成本降低了22%峰时购电功率下降了30%。需求响应让用户把负荷从高峰挪到低谷电能交互让IESO能利用地区间的价差灵活买卖电两个机制叠加实现了双赢。有一点必须提醒这些数字是特定参数下的结果不代表所有场景下都有这么大的提升。如果把上级电网的峰谷价差调小电能交互的收益就会变小如果把用户可转移负荷比例调低需求响应的作用也会降低。所以论文里做敏感性分析是有必要的这能验证模型在不同环境下的适应性。7. 实际代码调参中的踩坑记录7.1 迭代法 vs KKT单层法为什么我放弃迭代一开始我图省事用“交替迭代”实现主从博弈上层先给定初始价格求解下层得到负荷响应再把负荷响应代回上层用启发式规则更新价格重复直到收敛。结果发现价格经常在上下限之间来回振荡怎么调步长都压不住而且更新规则没有理论保证审稿人一问“怎么证明收敛”我就没法回答。后来老老实实改用KKT单层转化。虽然推导和编码量大了不少但求解过程是确定性的Gurobi给出的是全局最优解在给定离散变量取值下不需要担心振荡问题。如果你的目标只是快速验证想法迭代法可以跑着看趋势但要做严谨的科研或者实际项目KKT转化是更稳妥的路线。7.2 大M值选取和数值病态互补松弛线性化里的M值选择是影响求解稳定性的最直接因素。M取得太小可能把本该可行的解错误截断M取得太大会让约束矩阵条件数变差求解器陷入数值困难。我踩过的最夸张一次M取了1e8结果Gurobi直接警告“Matrix is ill-conditioned”求解出的电价出现很多肉眼可见的不合理尖峰。我的经验是M值根据对应约束中变量的数量级来取通常取该约束可能达到的最大值的2~3倍。比如功率平衡约束右侧最大负荷接近3 MWM取10就够对偶变量的M值可以按对偶变量理论上限取比如1000。不要图省事全模型统一用一个很大的M。必要的话可以先解一个RELAXATION版本统计原始约束的最大mismatch再据此定M。7.3 不可行问题的系统排查经验加了KKT和互补条件后模型规模激增不可行几乎是必然遇到的情况。我排查不可行时有一套固定流程先去掉所有互补松弛约束求解一个松弛问题看看原始约束和目标函数是否可行。如果这时候还不可行说明问题出在功率平衡、设备容量这些基础约束上检查和上层、下层约束定义是否冲突。如果松弛问题可行再逐一加入互补约束定位是哪一组互补约束导致不可行。通常是大M值太小或者是某个变量忘记定义非负。另外Yalmip默认的求解信息里会给出“infeasible”提示但不会直接告诉你哪条约束坏了。我习惯用yalmip(clear)之后把约束数组拆分成小块逐一验证配合check(Constraint)查看残差。这个方法比较笨但确实能救命。7.4 关于Yalmip的kkt命令Yalmip自带一个kkt函数可以自动生成大规模优化问题的KKT条件。听起来很方便但实测有几个限制它要求问题必须是凸的而且对变量类型和约束形式有比较严格的要求遇到含二进制变量的上层问题时直接把kkt用在目标里会报错生成的KKT条件会比较冗余导致求解性能下降。另外kkt生成的互补条件也是非线性仍需自己处理。所以我的建议是模型简单时可以用kkt快速验证思路但正式项目还是手动推导更可控。手动推导的过程虽然痛苦但能让你对模型结构理解更深后面扩展模型时也不至于一脸懵。以上这些经验都是我从最初能跑通简单算例到最终实现完整主从博弈模型的过程中一步步攒下来的。代码不是一次写成的模型也不是一步到位的。如果你现在正在折腾双层优化、多主体博弈建议先把一个最小的两主体、两时段模型跑通再慢慢扩展到24时段、多用户。小模型上的KKT推导正确了大模型只是“复制粘贴改索引”的事但前提是你真的理解了每一行KKT条件是从哪条约束来的。 
