合作型Stackelberg博弈微电网调度：Matlab+Cplex建模与求解

这些年一直在做微电网调度相关的仿真，经常碰到一圈讨论博弈论优化的人，嘴上说着纳什均衡、帕累托最优，手里拿的却还是单层确定性优化代码。真正要把一套“分布式主体互动”的模型落到可复现的Matlab+Cplex工程里，中间隔着一大段没人细讲的坑。这篇文章我就拿自己近期完整跑通的一套基于合作型Stackelberg博弈的微网运行策略来拆解，从为什么这么建模、数学模型长什么样，到Cplex求解前的每层变换、代码逻辑和调试经验，全部按复现顺序说一遍，给正在做微网优化、需求响应、储能调度的朋友一个能直接借鉴的工程参考。

这套策略解决的典型问题是：一个微网里，聚合商、光伏用户和储能之间利益不一致，聚合商要定内部电价，用户要决定用多少电、充多少电，双方互相影响。经典调度把用户负荷当成固定值，显然不合理；简单博弈模型又不考虑合作分配，算出来的方案没人愿意执行。合作型Stackelberg博弈恰好把“价格制定+负荷响应+收益再分配”闭环起来，而Matlab+Cplex则是把这种多层优化转成可解数学规划的最成熟组合。

1. 为什么微网运行策略需要Stackelberg博弈框架

1.1 一个网络，多个利益主体

微网的最大特征就是“角色多”。大电网侧有购电电价，微网内部有聚合商或运营商，再往下有带屋顶光伏的居民用户、有公共或私有充电桩、有储能系统、有柔性负荷。传统优化调度往往默认这一切都听一个“调度中心”指挥，目标函数写成一个总成本最小化，约束里把负荷曲线拍死在某个仿真日。这在物理上没问题，在经济学上却不成立——用户凭什么按照调度中心给的曲线用电？如果没有内部价格机制的激励，所谓“优化结果”在真实运行中根本执行不下去。

我在实际项目中做过一个对比测试：同样一套24小时数据，用集中式优化算出“最优”储能充放电计划和用户用电计划之后，再把计划里的电价拿出来让用户重新计算他自己的最优负荷，结果用户的真实响应和原先预想的方案差了将近15%。这个偏差在配网容量紧张的场景里足以导致过载。所以只要有多个决策主体，就不能回避“博弈”这个视角。

1.2 单层优化解决不了“价格—响应”闭环

单层模型里，内部电价通常是决策变量，同时用户负荷也是决策变量，二者被放进同一个目标函数和约束集。看起来自由度很大，但实际上违反了主从逻辑：电价和负荷不是并联关系，而是串联关系——聚合商先报电价，用户根据电价调整行为。

Stackelberg博弈的框架正好抓住这个先后的“承诺-响应”结构：领导者（聚合商）先做决策，宣布内部购售电价或激励价格；跟随者（用户、产消者）观察到价格后再做自身用能优化。求解出来的均衡点满足一个基本性质：给定领导者的价格策略，跟随者选的是最优响应；给定跟随者的响应规律，领导者也在约束范围内实现了自身收益最大化。

用生活类比来理解就是：商场先定折扣规则，顾客再决定买不买，商场不能既定折扣又替顾客决定购物清单。微网里就是这个逻辑。

1.3 合作型跟非合作型的本质差异

经典的Stackelberg模型偏“非合作”，领导者只在乎自己利润最大，跟随者只在乎自己成本最小，最后均衡点可能让某一方受益明显，另一方利益受损。对于配电网内的聚合商和居民用户来说，这种结果很难在现实中推行，因为用户会有抵触，甚至干脆不参与价格响应。

合作型Stackelberg在这里做了三个关键扩展：

• 把联盟整体收益也放进目标函数，领导者不会为了自己利润压榨用户到零。
• 引入合作剩余（cost savings / benefit surplus）的概念，即“大家一起协作”和“各干各的”之间的差额。
• 用Shapley值或其他合作分配规则把合作的额外收益在参与者之间分配，保证每个参与者分到的收益至少不比单独行动时差。

这相当于在非合作博弈的基础上加了一层“合作契约”。所以它在学术论文里看着高级，在工程上也很实用：能直接回答“用户为什么要参与需求响应”这个问题——因为参与之后分到的钱比不参与多。

2. 合作型Stackelberg博弈的数学模型搭建

2.1 上层领导者模型：微网聚合商的决策变量

上层角色我假设为微网聚合商，它负责从大网购电或向大网售电，经营储能，并制定内部电价信号。典型目标函数写成：

[ \max \sum_{t} \left( \text{售电收入}_t - \text{购电成本}_t - \text{储能运维成本}_t \right) ]

其中售电收入来自内部用户购电量和向大网售电量，购电成本来自大网电价和从内部光伏用户购电的电价。决策变量包括每个时段内部购/售电价 (\lambda_{t})、大网交换功率 (P_{grid,t})、储能充放电功率 (P_{ch,t}, P_{dis,t})。

约束要覆盖几类：

• 功率平衡：(P_{grid,t} + P_{pv,t} + P_{dis,t} = P_{load,t} + P_{ch,t})，写成允许弃光或允许甩负荷的松弛形式会更好。
• 储能SOC递推：(SOC_{t+1} = SOC_t + \eta_{ch} P_{ch,t} \Delta t - (1/\eta_{dis})P_{dis,t} \Delta t)，注意充放电效率别写成同一个符号。
• 容量约束：储能SOC上下限、充电功率上限、放电功率上限，以及“不能同时充放电”这类逻辑约束，如果变量连续可以用二者之和不超过上限来近似。
• 电价约束：内部电价一般限制在合理区间，比如大网电价的某倍范围，防止模型解出极端价格。

这里有一个我特别提醒的点：储能运维成本如果按充放电电量线性计费，做的时候很容易跟能量损耗混在一起，导致目标函数里重复扣费。最好只选其一，要么用折旧系数乘以绝对充放电量，要么用效率折损来体现损耗。

2.2 下层跟随者模型：用户的响应逻辑

下层跟随者可以建模为一个或多个用户群，每个用户群内部假设行为一致。用户的目标函数通常是最小化购电成本，加上对可转移负荷或EV充电的效用损失惩罚：

[ \min \sum_t \lambda_{t} P_{load,t}^{\text{总}} + \sum_t \omega (P_{shift,t} - P_{base,t})^2 ]

其中 (P_{shift}) 是用户的柔性负荷功率，(\omega) 是用户调整负荷的不舒适成本系数。用户约束包括：

• 负荷上下限：刚性负荷+柔性负荷的合理范围。
• 日电量约束：洗衣机、洗碗机这类可转移负荷一天总用电量不变，只是时段挪移。
• EV电池约束：到达时可充容量、离开时需要的最低SOC、最大充电功率。

为什么要写成一个二次惩罚项而不直接用分段线性？因为二次项让下层问题的Hessian矩阵正定，KKT条件简洁，也便于上层做迭代求解时保证下层最优解唯一。若下层问题线性但退化（多个最优解），KKT条件会捕获你不想要的那个解，后面调试会非常痛苦。这个经验值值得记住。

2.3 合作收益量化与Shapley值分配

合作博弈部分需要定义特征函数：任意一个参与者联盟 (S)，计算出该联盟单独运行时取得的总收益 (v(S))。比如：

• 用户单独行动：按大网零售电价购电，无内部优惠。
• 聚合商单独行动：无用户参与需求响应，只能按固定负荷购电。
• 全体联盟：采用Stackelberg均衡策略，实现系统总收益最大化。

Shapley值的计算原理是参与者 (i) 在所有可能联盟中的边际贡献期望：

[ \phi_i = \sum_{S \subseteq N \setminus {i}} \frac{|S|! (n-|S|-1)!}{n!} \left[ v(S \cup {i}) - v(S) \right] ]

这个公式对三个主体以内还能手算，对更多主体就组合爆炸了。工程上常用的做法是先跑主从优化，得到总收益增量，再按各主体边际贡献比例或按负荷参与量比例做近似分配。严格Shapley值适合结果后处理，不适合放进每个时段的优化迭代里，否则计算量受不了。

2.4 双层问题为什么难解，以及转单层的基本思路

上面的上层和下层合在一起是一个典型双层优化（bilevel optimization）。直接求解非常困难，因为下层问题是嵌套在约束里的优化问题，不是一组普通不等式。主流解法有两类：

• 第一类是直接用KKT条件替换下层问题，把双层转化为带互补约束的单层数学规划（MPEC），再通过大M法把互补约束线性化，交给MILP求解器。
• 第二类是用启发式迭代：上层固定电价，下层求响应；上层根据响应修正电价，循环到收敛。

我这次用的方案是第一类，因为Matlab+Cplex的组合对这种MPEC形式的MILP处理最直接。第二类在算力上省，但收敛到均衡的条件比较脆弱，适合做大规模场景下的快速近似。

3. Matlab + Cplex 工具链的前期准备

3.1 为什么选择Matlab+Cplex而不是Gurobi或自带求解器

Matlab自带优化工具箱的linprog、intlinprog对小规模教学问题完全够用，但双层转化后的MILP往往带有几百个0-1变量和上千条约束，自带求解器在速度和数值稳定性上明显不足。而且intlinprog对部分参数暴露不够细，想调分支策略、MIP gap、线程数都费劲。

Cplex作为商业求解器里最老的牌面之一，在电力系统领域积累很深。它对大型MILP、二次约束规划（MIQP）的支持非常成熟，尤其擅长处理带大M的0-1互补约束线性化问题。Cplex还有一个专门接口可以直接从Matlab调用，不需要写文件来回导入导出，对反复改参数做仿真实验很友好。

Gurobi确实也很强，某些稀疏问题上甚至比Cplex快，但如果你所在的课题组或公司已经有Cplex授权，或者你用的是Cplex社区版，那完全没必要同时维护两套接口。工程项目的真实约束往往不是“哪个求解器最强”，而是“哪个环境最容易跟团队现有代码跑通”。

3.2 环境配置与许可注意事项

我用的是Matlab R2023b配合Cplex 12.10的官方Matlab接口。装好之后验证环境是否正常，随手跑这几行：

try cplex = Cplex('test'); cplex.Model.sense = 'minimize'; cplex.Model.obj = [1; 2]; cplex.Model.A = sparse([1, 1]); cplex.Model.lhs = -inf; cplex.Model.rhs = 5; cplex.Model.lb = [0; 0]; cplex.solve(); disp(cplex.Solution.status); catch ME disp(ME.message); end

如果你看到1或者101这类状态码，说明求解器能正常启动。常见的坑有三个：

• Cplex官方文档里的addRows、addCols老接口仍然可用，但社区版经常提示许可证类型限制，求解超过一定规模变量数会报错。你需要在官网申请合适的学术或者社区版本许可，而不是用网上流传的零散文件。
• Cplex的Matlab接口对Matlab版本有兼容范围，不是你装最新的Cplex就一定支持最新Matlab。建议装之前查一下官方支持矩阵，否则Cplex类压根加载不出来。
• 很多人在旧版Matlab（比如R2016a）下用旧版Cplex还坚持用YALMIP，这没问题，但如果直接用Cplex对象会碰到数据结构不兼容的问题。我的建议是：能用纯Cplex接口就尽量不引入额外建模层，除非你的模型用YALMIP已经写好了不想重写。

3.3 建模前的数据准备清单

开工之前，把所有输入参数整理成标准结构体比什么都重要。我一般建一个inputData结构体，字段包括：

• timeHorizon = 24;仿真时段数
• loadBase = [...];每个时段的基线负荷
• pvOutput = [...];每个时段的光伏出力上限
• priceBuy = [...];大网购电电价
• priceSell = [...];大网售电电价
• essCapacity、essPowerMax、essSocMin、essMax、etaCh、etaDis
• flexRatio、evNum等用户侧参数

把所有参数集中在一个结构体里看起来多此一举，但在后面调参时有非常大的价值：你不需要在代码里找散落的数字，只需要改inputData的几个字段就能做敏感性分析。单位统一也很关键，我见过好几次因为MW和kW混用，优化器给出的“最优解”是物理上完全荒谬的数值却还显示在可行域内。

4. 核心实现过程与关键代码逻辑

4.1 整体求解框架：KKT变换加Cplex MILP

我的求解框架是“单层MPEC”路线。上层保持原样，下层问题拆出KKT条件塞进约束，叠加互补松弛线性化，最后形成一个单层MILP。流程上可以画成以下几个阶段：

1. 初始化输入数据，把每个时段、每个用户群的索引建好。
2. 构造上层目标函数，目标系数向量f，线性约束矩阵A。
3. 构造下层问题的拉格朗日函数，写KKT条件。
4. 把互补松弛条件线性化，需要引入辅助0-1变量。
5. 处理目标函数或约束里的双线性项。
6. 调用Cplex求解，读取结果，核算Shapley值。

这个框架最大的好处是求解一次得到的是一个数学上严格的均衡解，而不是通过迭代碰运气。缺点就是问题规模变大，变量数量大概增加三倍，但以Cplex的处理能力，微网这种几十上百个节点的规模完全在可接受范围内。

4.2 下层问题的KKT变换细节

下层如果是一个二次规划，形式为：

[ \min_x \frac{1}{2} x^T H x + c^T x \quad \text{s.t. } A x \le b ]

KKT条件包括：

• 平稳性：(H x + c + A^T u = 0)
• 原问题可行：(A x \le b)
• 对偶可行：(u \ge 0)
• 互补松弛：(u_i (b_i - A_{(i,:)} x) = 0)

其中只有互补松弛是非线性的，也不是凸的，需要用大M法转成混合整数线性约束：

[ \begin{array}{l} u_i \le M_u \cdot z_i \ b_i - A_{(i,:)} x \le M_s \cdot (1 - z_i) \end{array} ]

这里z_i是0-1变量。M_u和M_s这两组大M常数需要分别取合理上界：M_u可以从对偶变量的物理意义上估计，M_s则用对应松弛量的最大可能值来定。比如某个约束是储能不超过容量上限，那么松弛量最大就是容量上限本身，M_s取容量上界乘个1.5安全余量就够，千万不要统一取1e9。

大M选太大，数值求解会出问题；选太小，可能把实际可行解切掉。这是我调试中花时间最多的一个点，后面专门说。

4.3 上层目标函数与约束的Cplex建模

在Matlab里直接用Cplex对象建模，很多人第一次接触会觉得麻烦，但其实套路很固定。我在代码里一般是这样的逻辑：

nTotal = nCont + nBin; % 连续变量数 + 0-1变量数 f = zeros(nTotal, 1); % 目标系数 % 设置连续变量下界上界 lb = [-inf * ones(nCont,1); zeros(nBin,1)]; ub = [inf * ones(nCont,1); ones(nBin,1)]; % 不等式约束：Aineq * x <= bineq Aineq = []; bineq = []; % 等式约束：Aeq * x == beq Aeq = []; beq = []; % 变量类型：'C' 表示连续，'B' 表示0-1 ctype = [repmat('C', nCont, 1); repmat('B', nBin, 1)];

然后统一交给Cplex：

cplex = Cplex('stackelberg_mpec'); cplex.Model.sense = 'minimize'; cplex.Model.obj = f; cplex.Model.lb = lb; cplex.Model.ub = ub; cplex.Model.A = [Aineq; Aeq]; cplex.Model.lhs = [-inf * ones(size(bineq,1),1); beq]; cplex.Model.rhs = [bineq; beq]; cplex.Model.ctype = ctype; % 设置MIP gap和求解时间上限 cplex.Param.mip.tolerances.mipgap.Cur = 1e-4; cplex.Param.timelimit.Cur = 600; cplex.solve(); if cplex.Solution.status == 101 || cplex.Solution.status == 102 x_opt = cplex.Solution.x; else error('模型求解失败，状态码%d', cplex.Solution.status); end

这个模板可以直接套用。lhs和rhs的组合是Cplex的一种风格，既支持等式也支持不等式，比分开添加更简洁。注意目标函数里如果有常数项，直接加在cplex.Model.obj对应的连续变量系数里，而不是额外加一个常数变量。

4.4 双线性项线性化的处理

上层目标函数里，聚合商的内部售电收入是内部电价乘以用户负荷，也就是 (\lambda_t \times P_{load,t})，这天然是双线性的。如果直接把这个放进MILP，Cplex会拒绝，因为这是非凸二次项。

处理思路分两种：

• 如果 (\lambda_t) 和 (P_{load,t}) 其中一个被KKT条件替换后变成对偶变量的函数，有时可以通过强对偶定理消掉乘积项。具体是利用下层问题的最优值等于其对偶最优值，把下层目标函数替换成只含对偶变量的表达式，从而把上层目标中的某些双线性项变成线性项。
• 如果消不掉，就要用McCormick包络做四边形松弛，或者按价格离散成多个区间做分段线性近似。对每个时段拆成 (K) 个价格区间，(\lambda_t) 写成 (\sum_k \lambda_{t,k} \cdot \delta_{t,k})，其中 (\delta_{t,k}) 是0-1变量，再引入辅助变量表示 (\lambda_t \cdot P_{load,t}) 的线性上下界。

我实际测试下来，能用强对偶消项就优先消项，因为McCormick松弛会增加很多变量，还会带来松弛误差。只有在某个约束里实在消不掉，才做分段线性化，而且分段数控制在4到6段比较合适，太多段求解时间成倍上升。

4.5 SOC递推与双向功率约束的建模细节

储能建模是实现中出错率最高的地方。SOC递推本身是线性等式，但充放电同时开启的问题需要处理。如果用两套非负变量建模，必须加“不能同时充放电”约束。最直接的方式：

[ P_{ch,t} \le P_{ch}^{\max} \cdot z_{ch,t},\quad P_{dis,t} \le P_{dis}^{\max} \cdot z_{dis,t},\quad z_{ch,t} + z_{dis,t} \le 1 ]

其中 (z_{ch,t}, z_{dis,t}) 是0-1变量。如果问题规模大，不想因此引入太多0-1变量，可以用一个可变“净功率”变量 (P_{ess,t}) 表示正为充电、负为放电，SOC递推里同时用效率系数分段表达：

[ SOC_{t+1} = SOC_t + P_{ess,t}^+ \eta_{ch} \Delta t + P_{ess,t}^- \frac{1}{\eta_{dis}} \Delta t ]

再把 (P_{ess,t}) 分解成正负两部分，这样能减少一半的0-1变量，缺点是效率带来的非线性仍然要用分段或者整型变量表达。我的经验是：如果是24小时单阶段求解，直接用两套变量加0-1约束更清晰，代码可读性好，调试也方便；如果你要滚动优化跑几个月，才考虑压缩变量数量。

4.6 求解结果的后处理与收敛质量检查

Cplex求解完之后，不要只看Solution.status等于101就直接收数据。我一般做三层检查：

• 第一层，检查cplex.Solution.bestobjective和gap，确认MIP gap真的到了设定阈值，而不是靠时间限制强制截止。
• 第二层，把解出来的负荷、电价、光伏出力代回原始约束里逐一验证，至少跑一个校验脚本输出“最大功率不平衡量”和“SOC越限量”。
• 第三层，把各主体的收益分别计算一遍，和合作前做对比，确认没有参与者收益下降。如果某个用户算出来收益比单独行动还差，说明Shapley值分配环节权重设置有问题。

我遇到过不少次“求解器说最优但物理上不对”的情况，最后都是出在校验脚本没跑或者校验逻辑写错上。所以强烈建议把校验脚本当成和主程序同等重要的代码来写。

5. 常见问题与排查经验实录

5.1 “求解不可行”的7个高频原因

不可行问题在MPEC中最常见，我列一个排查顺序表，按照这个表去查基本能定位80%的问题：

调试时我常常用到一个小技巧：先固定某些变量，比如把0-1变量全部固定到一组试解，再把互补约束去掉，看看剩余LP是否可行。如果LP可行，说明问题出在整数因素或互补松弛上；如果LP都不可行，那就是原问题约束冲突，先去查功率平衡和边界条件。

5.2 Cplex求解速度慢怎么办

单层MPEC转化后的MILP规模可能是原问题的3到5倍，一旦出现求解慢，不要急着压缩模型，先按下面顺序做优化：

• 打开Cplex的日志，看是下界上不去，还是上界下不来。如果下界上升缓慢，多半是MILP的LP松弛太弱，需要加有效不等式或者加强大M；如果上界下降慢，问题在找到的好整数解太少，这时候要设cplex.Param.mip.strategy.nodeselect.Cur = 2，让搜索更偏向深度优先。
• 减少非必要的0-1变量。储能“不能同时充放电”这种约束，如果能用连续变量加一个小的惩罚项达到近似效果，就尽量别引入整数变量。
• 给求解器一个热启动初始解。我先用上一轮迭代中得到的Stackelberg解赋值给cplex.Model.start.x，这样Cplex可以在MIP搜索早期就切掉大量分支。实测某些算例能提速两三倍。

cplex = Cplex('mipHotStart'); % 设置模型... cplex.start.x = x_start; cplex.solve();

• 收紧MIP gap。如果工程上不需要1e-6级最优，把mipgap设在5e-4到1e-3，求解时间能差出一个数量级。

5.3 数值稳定性：为什么同样的模型换个数据就出错

另一个让我印象极深的坑是数值稳定性。同样的代码，用一组平滑负荷数据跑得好好的，换一组尖峰特别明显的负荷数据，突然报“infeasible”或者“unbounded”。查来查去，往往问题出在约束矩阵的条件数上。

大M变量一旦取值偏大，矩阵数值跨度会到1e6甚至1e9，Cplex内部预求解器很可能直接把这些行判成病态。处理办法：

• 尽可能把单位统一成“标幺值”或相同数量级。比如功率用100kW作为基准的标幺值，电价用元/kWh的10倍作为标幺。
• 大M值逐个计算最小可行值，而不是统一一个很大的数。
• 如果还不行，用cplex.Param.emphasis.numerical.Cur = 1开启数值稳定性侧重模式，但这会在速度上有些代价。

6. 结果解读、实验对比与几个扩展方向

6.1 如何衡量这套策略真的比传统方案好

代码跑通之后，我会做三组横向对比：

第一组是“集中式优化”：把所有主体合并成一个利益体，不存在博弈，只有一个总成本最小。这组结果作为理论上限。

第二组是“非合作Stackelberg”：聚合商只优化自身利润，用户只响应价格，但不做收益共享和合作分配。这组结果代表没有合作契约的均衡状态。

第三组是“合作型Stackelberg”：完整模型，包含联盟收益最大化、KKT下层响应和Shapley值分配。

实际数据跑下来，集中式优化的系统总成本最低，但用户侧满意度往往差；非合作Stackelberg下聚合商利润最高，但用户成本反而可能比不参与更高；合作型Stackelberg的总成本介于前两者之间，但系统总剩余（即用户效用加聚合商利润）最优，而且按Shapley值分配后各方都比单独行动时好。这个结果其实很有说服力——它说明博弈优化不是为了追求一个极端最优值，而是为了找到一个能落地的稳定方案。

在评估指标上，除了系统总成本和用户成本，建议再统计几个工程指标：

• 储能循环次数，看合作机制是否让储能利用率更合理。
• 负荷峰谷差，如果电价信号有效，峰时负荷应该被抑制。
• 弃光率，合作型策略下光伏用户应该更愿意在光伏大发时段调整用电。

6.2 后续可以怎么扩展这套代码

这套模型最大的扩展价值在于它的框架是通用的。我已经在好几个方向上试过延伸，改动都不需要伤筋动骨：

• 滚动时域：把24小时模型改成分时段滚动，每次只执行前1个小时的策略，然后重新优化，可以应对光伏预测误差。
• 多层博弈：在聚合商之上再加一个配网运营商的电价信号，形成三层Stackelberg结构，数学处理上再多一层KKT替换。
• 鲁棒优化：把光伏出力和负荷预测误差建模成不确定集，在下层问题的参数上做区间扰动，Cplex依然能处理线性化的鲁棒对等模型。
• EV充放电协调：把用户层细化成电动汽车集群，增加可移动负荷约束，就可以直接当作车网互动调度策略使用。

我个人在实际操作中的体会是，博弈类微网模型最难的不是数学公式，也不是Cplex调用，而是你能不能把现实中的利益关系准确翻译成目标函数和约束。翻译对了，后面的KKT变换和线性化全是体力活；翻译错了，再花哨的算法也是空转。建议刚上手的读者先从两主体简单场景做起，确认均衡逻辑通顺之后，再逐步加储能、加多用户、加合作分配，这套路踩坑最少。