最近刚把一个EI论文里的"基于主从博弈的新型城镇配电系统产消者竞价策略"在IEEE33节点系统上完整复现了一遍，Matlab代码从双层模型搭建到KKT条件转化，再到CPLEX求解和结果验证，前前后后折腾了不少时间。这个方向确实是当下的热点——分布式光伏、储能、充电桩大规模涌入配电系统，传统"用户"变成了"产消者"，既买电又卖电，配电网的运行逻辑整个变了。但我也发现一个问题：论文很多，真正能把主从博弈从公式落地成可运行代码、并讲清楚每个环节为什么这么做的资料却很少。这篇文章就把我复现的全过程、模型推导和踩过的坑一次说透，适合正在做配电网优化、电力市场、需求响应方向的研究生和工程师参考。
1. 为什么新型城镇配电系统的竞价问题非得用主从博弈
1.1 从"被动负荷"到"主动产消者"，优化问题的属性变了
传统的配电系统优化，不管是最优潮流还是经济调度，本质上都是一个"集中式决策"问题：调度中心掌握全部负荷和电源信息，直接下发指令，所有用户按指令执行。这个模式成立的前提是用户是"被动负荷"——只能按固定曲线用电，没有任何自主调节能力。
但新型城镇配电系统里出现了大量分布式光伏、储能、电动汽车充电桩，情况完全不同了。一个装了屋顶光伏和家庭储能的小区用户，电价高的时候可以少买电甚至卖电，电价低的时候可以把储能充满留到高峰用。这类用户具备"生产+消费"双重属性，圈内叫prosumer（产消者）。他们不再是执行指令的被动节点，而是会基于价格信号最大化自身利益的独立决策主体。
这时候再用"调度中心一杆子管到底"的思路就不成立了。你发指令让某个用户晚上必须充多少电，人家不一定乐意，因为影响的是他自己的电费账单。本质上，这个系统的决策权已经分散到了产消者手里，优化问题从单层变成了多层。
1.2 双层博弈恰好匹配配电市场的"先后决策"结构
既然决策权分散了，那怎么描述DSO（配电系统运营商）和产消者之间的互动关系？这就回到博弈论。
主从博弈，学术上叫Stackelberg博弈，刻画的是一个"Leader-Follower"结构：上层决策者先行动，公布一个策略（比如电价）；下层决策者观察到上层策略后，做出自己的最优响应；上层再根据下层的响应调整策略，最终达到一个双方都不愿单方面改变的均衡状态。
你细品一下，这个结构跟实际的配电市场交易顺序几乎一模一样：配电公司先公布内部结算电价或者交易规则，产消者看到价格之后决定多买还是多卖、储能充还是放。所以用主从博弈来描述"配电运营商定价-产消者响应"的互动关系，非常自然。
相比之下，其他博弈模型就不太合适。Nash博弈假设所有参与者地位平等同时决策，但实际上DSO拥有网络信息和定价权，产消者没有；合作博弈需要设计利润分配机制，而且默认大家会遵守联盟约定，这在市场环境下不太现实。主从博弈其实反映了电力市场中天然存在的"时序差异"和"信息不对称"，这是它成为该问题主流建模方法的根本原因。
1.3 竞价策略的落点：内部出清电价怎么定
论文标题里的"竞价策略"直白讲就是：DSO制定内部市场的结算电价，产消者根据电价各自优化购售电决策。这里面有个关键点——DSO不是简单配电网公司卖电给用户，而是在一个内部电能交易市场里扮演"做市商"角色，既要保障系统安全运行，又要激励产消者主动参与调节。
所以整个模型追求的目标是双重的：产消者实现自身效益最大化，DSO通过电价信号引导所有主体的行为，最终使配电系统运行在安全且高效的运行点上。
2. IEEE33节点测试系统建模：网络参数与产消者接入方案
2.1 33节点系统的底子参数
IEEE33节点系统是配电网研究里最经典的测试算例，没有之一。它是一个12.66 kV的辐射状配电网，33个节点、32条支路，总负荷大概有功3715 kW、无功2300 kVAr，末端还带5条联络开关。做复现之前，先要把这个系统的参数表拉出来看清楚，尤其是支路阻抗和节点负荷这两块，大部分复现结果对不上都是因为这两组数据抄错了。
常用基准参数如下表：
| 参数项 | 数值 |
|---|---|
| 基准电压 | 12.66 kV |
| 基准功率 | 1 MVA（或10 MVA，看论文定义） |
| 节点数 | 33 |
| 支路数 | 32 |
| 总有功负荷 | 3715 kW |
| 总无功负荷 | 2300 kVAr |
| 联络开关 | 5条 |
这里提醒一句：论文里所有标幺值换算都要先确认基准功率选的是1 MVA还是10 MVA，不然潮流结果整体差一个数量级，后面KKT条件、约束越限全都对不上。
2.2 潮流方程用DistFlow而不用牛拉法
复现这种双层优化模型，潮流约束必须嵌入到优化问题里一起求解。如果用传统牛顿-拉夫逊潮流，那是一套非线性方程组，没法直接丢给优化求解器。所以主流论文都用DistFlow方程描述辐射状配电网潮流：
P_{j} = P_{i} - r_{ij}\frac{P_{i}^{2}+Q_{i}^{2}}{V_{i}^{2}} - P_{load,j}
Q_{j} = Q_{i} - x_{ij}\frac{P_{i}^{2}+Q_{i}^{2}}{V_{i}^{2}} - Q_{load,j}
V_{j}^{2} = V_{i}^{2} - 2(r_{ij}P_{i} + x_{ij}Q_{i})
其中r和x是线路阻抗，P和Q是支路有功无功功率。这个方程的好处是结构清晰，而且可以通过忽略网损项（即\frac{P_{i}^{2}+Q_{i}^{2}}{V_{i}^{2}}近似为0）得到线性DistFlow，直接把潮流方程变成线性约束，求解稳定性和速度都有很大提升。
在实际复现中我建议先用线性DistFlow跑通整体框架，再根据论文需要决定要不要加回网损项。多数EI论文里的算例采用线性化处理，因为主从博弈转换后本身已经很复杂，没有必要在潮流精度上给自己找麻烦。
2.3 产消者节点怎么选：不是随便挑的
IEEE33节点系统里哪些节点可以作为产消者的接入点？这不是随便选几个就行，要从物理背景出发。一般来说，选末端重负荷节点或者有分布式光伏接入条件的节点更合理。节点18、25、33这些靠近馈线末端的节点电压偏低，接入光伏和储能在调压和降损上的效果最明显，产消者参与竞价的潜力也最大。
我在复现时设了3个产消者集群，分别挂在节点18、25、33，每个集群内包含光伏、储能和可调负荷。这样设置的好处是：节点分布在不同的馈线分支上，产消者之间的距离足够远，不会因为电气距离太近导致电价信号互相干扰太严重，博弈结果更容易体现出策略差异。
2.4 产消者侧设备建模的三个关键块
每个产消者内部的设备模型分三块：
光伏模型。光伏出力按照典型日光照曲线给定，实际上就是一条时间序列P_{pv}(t)。复现时注意光照曲线是不是和电网负荷曲线同一天，很多论文用夏季典型日，你换冬季数据结果肯定不一样。
储能模型。储能约束是最容易出问题的地方，包括SOC（荷电状态）转移方程、充放电功率上限、以及最关键的"不能同时充放电"约束。SOC方程是：
SOC_{t+1} = SOC_{t} + (\eta_{ch}P_{ch,t} - \frac{P_{dis,t}}{\eta_{dis}})\Delta t
充放电互斥可以用二进制变量处理，这一步在后面线性化时要非常小心。
可调负荷。可调负荷的功率在一定范围内连续可调，或者按需求响应方式把部分负荷转移到低电价时段。对产消者来说，可调负荷是应对电价波动、优化自身成本的重要抓手。
产消者的效用函数通常用二次函数表示，这保证了后面KKT转化时下层问题的凸性，这个性质非常重要，后面第4节细讲。
3. 主从博弈的双层数学模型：逐层拆解目标函数和约束
3.1 上层：配电系统运营商D怎么决策
上层的决策者是配电系统运营商（用D表示），它的核心决策变量是内部结算电价\lambda(t)，即每个时段向产消者购电和售电的统一价格。
上层的目标函数需要和论文核对清楚，常见的有两种：一种是追求DSO自身收益最大化（低买高卖赚价差），另一种是社会福利最大化或者系统运行成本最小化。不同目标函数直接决定了博弈均衡点完全不同，复现前先看论文第几节把目标函数写清楚了。
除了目标函数，上层还背着网络运行约束：节点电压上下限（一般0.95~1.05 p.u.）、支路容量约束、以及DSO从上级电网购电的功率限制。这些约束加进去，电价就不能想定多少定多少——如果某个电价会导致电压越限，这个解就得被砍掉。
3.2 下层：产消者的收益最大化问题
每个产消者i在给定电价策略\lambda(t)后，求解自己的优化问题。决策变量包括：从电网购电功率、向电网售电功率、储能充放电功率、可调负荷的用电量。目标函数通常是最大化自身净收益：
效用(用电) - 购电成本 + 售电收入
这里有个细节：产消者的售电价格和购电价格可能不同，有的论文直接用统一出清价格\lambda，也有的论文设置购电价\lambda_{buy}和售电价\lambda_{sell}有差值。这个差值是DSO盈利空间的来源，也是储能套利的驱动因素。如果做成完全相同的价格，储能在经济性上就没有套利动机了。
3.3 竞价机制的衔接规则
主从博弈里的"竞价机制"具体指向的是内部市场怎么结算。复现时最常用的是统一出清价格机制（CPS，Clearing Price Settlement）：所有成交的购电和售电都按同一时段统一的结算价格\lambda(t)来结。这就比PAB（Pay-As-Bid，按报价支付）要简单得多，模型写起来也清爽。
产消者的功率平衡约束是：
P_{pv,i}(t) + P_{dis,i}(t) + P_{buy,i}(t) = P_{load,i}(t) + P_{ch,i(t)} + P_{sell,i}(t)
也就是说光伏出力加储能放电加购电，要等于负荷加储能充电加卖电。这个等式约束在下层优化中会对应一个很重要的对偶变量，后面KKT条件里用得到。
3.4 博弈均衡的数学含义
这个主从博弈的均衡解叫做Stackelberg均衡：在均衡状态下，DSO给定的电价策略是在考虑产消者最优响应后的最优定价；同时产消者在给定电价下确实选择了自己的最优策略，没有任何一方可以通过单方面改变策略而获益。
实现均衡的前提是下层问题对每个产消者都是凸优化问题（或者能转化成凸问题），这样才能保证KKT条件的必要性，进而通过KKT条件把下层问题"编码"进上层问题的约束里。这也是为什么前面说要选二次效用函数——一旦选了非凸的目标函数，整个博弈均衡的存在性都会出问题。
4. KKT条件转化与大M线性化：把双层博弈变成可计算的形式
4.1 为什么不能直接写个循环迭代求解
很多人第一次看到双层模型，脑子里冒出的想法是：外层定电价，内层解优化，循环迭代不就行了吗？这个思路没问题，但要注意论文发表时审稿人一般不会接受纯启发式迭代——因为迭代过程不保证收敛到均衡，也不保证解的最优性。EI论文里的主从博弈求解，主流做法是：利用下层问题的KKT条件把下层优化"嵌入"上层优化，把双层问题转化为单层数学规划，再交给商业求解器一次求解。
4.2 KKT条件的四大件
把下层优化问题写成标准形式：
min_{x} f(x), s.t. Ax \leq b
其KKT条件包括四部分：
- 稳定性条件（Stationarity）：拉格朗日函数对各决策变量的偏导为零
- 原始可行性（Primal feasibility）：Ax \leq b
- 对偶可行性（Dual feasibility）：\lambda \geq 0
- 互补松弛条件（Complementary slackness）：\lambda_{i}(Ax - b)_{i} = 0，对每个约束i成立
对于凸优化问题，KKT条件是充分必要条件，所以下层优化可以被这组KKT条件完全等价替换。
实操中这里有个最常踩的坑：遗漏对偶变量。公式推导时手写的拉格朗日函数少写了一项，或者约束漏了一个非负限制，都会导致KKT条件不完整。我建议在Matlab里建模时，把下层每个约束都对应一个对偶变量，并在代码注释里逐一标注映射关系，这样后期查错能省大量时间。
4.3 互补松弛条件的线性化：大M法
互补松弛条件\lambda_{i}(Ax - b)_{i} = 0是一个双线性等式（两个非负量的乘积为零），这是非线性的，求解器处理不了。标准处理手法是大M法：
引入二进制变量z_{i}，将互补松弛条件转为两组线性约束：
\lambda_{i} \leq M z_{i}
Ax - b \leq M (1 - z_{i})
逻辑是：如果\lambda_{i} > 0则z_{i}=1从而(Ax - b){i}必须为0；反过来如果(Ax - b){i} < 0则z_{i}=0从而\lambda_{i}必须为0。这样就实现了两个量至少一个为0的效果。
M值怎么取是个学问。M太小会砍掉可行域，导致求出的解根本不在真实可行域里；M太大会让松弛问题产生数值病态，尤其CPLEX处理大M值大整数变量时经常报数值警告。我的经验是先看等式/不等式各变量的理论上下界，然后取上界值的10~100倍，再用一组已知可行解回测，不断调整。
4.4 上层目标里的双线性项也得处理
KKT转化嵌入下层后，上层目标函数里往往还残留着电价乘以功率的双线性项——比如DSO的收入是\lambda(t)乘以总成交电量，这里的\lambda和电量都是变量，两个变量相乘又是个非线性项。
处理思路是利用强对偶定理。如果下层问题满足强对偶（凸+Slater条件），那下层原问题的最优目标值等于对偶问题的最优目标值。利用这个等量关系，可以把下层目标函数里的双线性项替换成对偶表达式，从而把上层目标函数转化为线性形式。这一步是很多复现代码里最绕的地方，推完之后整个模型就变成一个混合整数线性规划（MILP），交给CPLEX/Gurobi就非常好解了。
4.5 最终单层MILP的问题结构
做完KKT转化加线性化，原本的双层博弈变成一个单层MILP：
- 变量：上层电价变量+下层决策变量+下层对偶变量+大M二进制变量+储能充放电互斥二进制变量
- 目标：DSO的目标函数（已经线性化）
- 约束：潮流约束、设备约束、KKT条件（转化成线性约束）
这里要留意问题的规模：IEEE33节点配上十几个时段，再加上二进制变量，问题规模不小。如果CPLEX默认参数跑起来太慢，建议把MIP gap设为0.5%~1%，结果差别很小但速度能快不少。
有些论文用"数学规划+均衡约束"（MPEC）来描述这个单层问题，本质上是一样的东西，只是叫法不同。
5. Matlab代码实现：从数据初始化到结果输出的完整流程
5.1 环境配置与工具箱选择
Matlab代码实现主从博弈，离不开YALMIP做建模层，求解器用CPLEX或Gurobi。推荐环境是Matlab R2020b以上版本，YALMIP用Github上的最新发布版，CPLEX 12.10或Gurobi 9.5以上。这几个版本之间兼容性比较成熟，不容易出奇怪的报错。
建模之前先确认求解器装好了，YALMIP里用solvesdp或者optimize的时候能识别到求解器。命令很方便：
% 检查求解器可用性 yalmiptest如果CPLEX没识别出来，八成是YALMIP找不到CPLEX的路径，需要手动addpath到CPLEX的matlab目录。
5.2 代码结构全景
我的复现代码分五个文件块，思路是按"数据-建模-求解-后处理"分离：
main.m % 主程序，串联整个流程 case33_data.m % IEEE33节点参数、产消者设备参数 build_upper_problem.m % 上层DSO模型 build_lower_kkt.m % 下层产消者模型+KKT条件转化 plot_results.m % 结果可视化主程序里的核心逻辑是：先加载数据，再构建上层目标函数和约束，然后构建下层KKT约束（等价替换下层问题），加上大M线性化辅助变量，最后yalmip的optimize求解，输出电价和产消者策略。
5.3 核心建模代码片段
用YALMIP定义变量和约束的代码长这样：
% 加载33节点系统数据 mpc = case33_data(); % 定义决策变量 lambda = sdpvar(n_prosumer, T); % 内部结算电价 p_buy = sdpvar(n_prosumer, T); % 产消者购电 p_sell = sdpvar(n_prosumer, T); % 产消者售电 p_ch = sdpvar(n_prosumer, T); % 储能充电 p_dis = sdpvar(n_prosumer, T); % 储能放电 soc = sdpvar(n_prosumer, T); % 储能SOC u_ch = binvar(n_prosumer, T); % 充电互斥二进制 u_dis = binvar(n_prosumer, T); % 放电互斥二进制 % 再加上下层对偶变量... lambda_dual = sdpvar(size(A_lower,1), T); % 下层约束对偶变量定义完变量后，就是构建约束集合。下层KKT条件的稳定性条件用YALMIP写大概是这样：
% 构造拉格朗日函数（简化示意） L = sum(f_lower) + sum(lambda_dual .* (A_lower * x_lower - b_lower)); % 稳定性条件：L对x的梯度为0 Constraints = [Constraints, jacobian(L, x_lower) == 0]; % 对偶可行性 Constraints = [Constraints, lambda_dual >= 0]; % 互补松弛条件用大M法 M_val = 100; % 根据实际数据量级调整 z = binvar(size(A_lower,1), T); Constraints = [Constraints, lambda_dual <= M_val*z]; Constraints = [Constraints, A_lower*x_lower - b_lower <= M_val*(1-z)];这里用jacobian函数对向量求偏导，在YALMIP里可直接写成等式约束，不需要手动展开每一项梯度，代码会清爽很多。当然梯度项也不复杂，只是手写容易出错。
5.4 求解器参数与调用
模型构建完成后直接调用求解函数：
% 配置求解器参数 options = sdpsettings('solver','cplex',... 'verbose',2,... 'cplex.mip.tolerances.mipgap',0.005,... 'cplex.mip.tolerances.integrality',1e-6,... 'cpu',5); % 求解 sol = optimize(Constraints, Objective, options);跑完之后先判断sol.problem是否为零，非零说明模型有问题，我用一个统一检查函数把YALMIP的problem code映射成可读错误信息，排错效率高很多。
5.5 迭代式求解的替代路线（临时方案）
如果你只是想快速看到趋势效果、或者KKT转化一时推不出来，可以先用启发式迭代法跑通流程：初始化电价，固定电价求各产消者最优策略，再把产消者总响应带回上层重新定价，反复迭代直到电价变化小于阈值。
% 主从博弈数值迭代（简化示意） lambda = 0.4 * ones(n_prosumer, T); % 初始电价 damping = 0.5; % 阻尼系数，防震荡 while norm(delta_lambda, 'fro') > 1e-4 && iter < 50 % 下层：给定lambda，各产消者自优化 for i = 1:n_prosumer [p_buy(i,:), p_sell(i,:), soc(i,:)] = solve_prosumer(i, lambda(i,:)); end % 上层：根据总响应重新定价 lambda_new = solve_dso_pricing(p_buy, p_sell); % 阻尼更新，避免来回震荡 delta_lambda = lambda_new - lambda; lambda = lambda + damping * delta_lambda; iter = iter + 1; end这个方案的问题在于：不保证收敛到Stackelberg均衡，结果不能作为正式复现结论。我的建议是把它当作调试辅助工具，理解博弈互动的直观动态可以，最终以KKT转化+MILP的结果为准。
5.6 结果输出与对标检查
求解完成后，把关键结果输出成图表：
- 内部结算电价曲线（24时段）
- 各产消者的购售电策略和净收益
- 储能SOC变化曲线
- 节点电压分布（重点检查是否在0.95~1.05 p.u.范围内）
- 配电网网损，和基准场景（无竞价机制）做对比
对标的重点是：论文算例表格里的DSO收益、产消者收益、网损下降比例、电压改善幅度这几个数字，如果差得远，大概率是参数或曲线数据不一致，需要回到第2节排查。
6. 复现中的典型坑：不收敛、奇异解与对偶变量丢失
6.1 对偶变量遗漏导致的"模型无解"
这是KKT转化路线里最隐蔽的坑。拉格朗日函数漏写一个约束项，稳定性条件就不完整，但求解器不一定报错——它可能求出一个满足所有约束但根本不是原问题最优解的解，或者直接报infeasible。如果是infeasible还好查，最怕的是"看起来有解但结果诡异"。
排查思路：把下层原问题和KKT转化后的单层问题分别求解，对比下层目标函数值。如果跟原问题的精确最优解对不上，说明KKT条件有遗漏。我在代码里做了一个自动测试模块：随机抽取一百组电价，分别用"直接求解下层原问题"和"用KKT约束求下层"，比对最优目标值和决策变量，误差大于阈值就报错定位。
6.2 大M值的量级陷阱
大M法处理互补松弛条件，M取值过大会导致CPLEX数值不稳定，过小会截断可行域。我最初用M=10000，CPLEX一直报数值警告，部分互补松弛条件趋近于不满足；后来把M缩小到100~200量级（根据电量上界和电价上界估算），数值稳定很多。
实操建议是：先根据约束表达式的物理上界算出一个理论M值，再上下调整各一个数量级做对比实验。如果结果不随M变化，说明M取值在这个范围内是安全的；如果结果明显依赖M，说明M的取值有问题，需要重估。
6.3 迭代法不收敛的两个常见原因
如果走数值迭代路线，最常见的不收敛原因是：初始电价离均衡点太远，产消者响应在迭代中来回震荡；以及没有阻尼系数，迭代轨迹发散。给电价更新加个0.3~0.6的阻尼系数能有效缓解震荡；
另外注意收敛判据不要只看电价差，还要看产消者策略是否稳定。我见过电价已经收敛但储能SOC还在缓慢漂移的情况，这其实是收敛判据太宽了。建议同时检查电价和储能策略两个量的相对变化。
6.4 跟论文结果对不上时的排查顺序
复现结果和论文数值对不上，不要急着改参数硬凑。我总结的排查顺序是：
- 基准值检查——基准功率是不是1MVA，转为标幺值的基准是不是一致
- 负荷曲线检查——论文用的负荷曲线峰值、形状和你的数据是否一致
- 光照曲线检查——光伏出力峰值时刻、峰值功率是否一致
- 储能的初始SOC——很多论文初值设0.5，有的设0.2，终点SOC有没有约束，直接影响储能收益
- 结算机制——CPS还是PAB，购售电价有没有差值
- 对偶问题检查——强对偶替换用的对偶问题有没有写对
每检查一项就重新跑一遍看结果偏移方向，不要多个因素同时改，否则你根本不知道是哪个因素导致的对不上。
6.5 结果可信度的自检清单
复现完成后的最后一步，一定要做结果合理性自检：
| 检查项 | 合格标准 |
|---|---|
| 节点电压 | 所有节点电压在0.95~1.05 p.u. |
| 支路潮流 | 不过载 |
| 储能SOC | 全过程在0.1~0.9范围内，满足转移方程 |
| 充放电互斥 | 同一时段不会同时充电和放电 |
| 市场出清条件 | 产消者总购电量与总售电量之差等于上级电网购电量 |
| 互补松弛余量 | 所有互补松弛条件乘积接近0（小于小阈值） |
如果这些检查项全部通过，基本可以确认复现的代码逻辑是正确的。
我在实际复现里还有一个体会：主从博弈这类模型，公式推导花的时间只占三成，剩下七成都在跟"非线性项线性化""求解器数值问题""对偶变量的调试"搏斗。如果你也是第一次复现，建议先从单时段简化版跑通整个代码框架，再加多时段和完整约束——一步一步来，比一口气写完再回头debug效率高得多。另外储能充放电互斥的二进制变量尽量不要省，省掉的话模型在低价时段会同时充电放电"空转"，白白消耗能量，结果出来你都不知道哪里错了。
