我先梳理一下这个项目的核心脉络。你拿到的是一个EI论文复现代码，主线是用主从博弈（Stackelberg博弈）来刻画新型城镇配电系统里配电网运营商和产消者之间的竞价互动，算例用的是IEEE33节点系统，最终落地在Matlab里跑出均衡结果。这类项目在电力市场与配电系统交叉领域非常典型，尤其是近年来分布式光伏、储能、电动汽车大批量接入后，产消者（prosumer）这个概念从理论走进了工程实践，竞价策略研究也随之变成了热点。本文不打算复述论文本身，而是以“复现者”的视角，把模型原理、代码实现、踩坑记录和加速技巧全部摊开讲。
1. 项目背景与核心思路拆解
1.1 新型城镇配电系统与产消者竞价要解决什么问题
传统配电网里，电能从变电站单向流向用户，用户是被动消费者。但新型城镇配电系统的典型特征是分布式光伏、小型风电、储能系统和电动汽车大量接入，用户不再是单纯的“用电方”，而是既能从电网购电、也能向电网卖电的“产消者”。这就带来一个很实际的问题：电网运营商和产消者之间的电价怎么定？产消者手里的分布式电源和储能怎么调度？双方都有各自的利益诉求，存在明显的“你定电价、我定策略”的先后决策关系，天然适合用主从博弈来建模。
再说得直白一点，主从博弈就是“上层先出牌、下层后出牌但反过来制约上层收益”的博弈过程。电网运营商作为领导者（leader），制定分时电价或内部结算电价；产消者作为跟随者（follower），在给定电价下优化自己的购售电策略，实现自身收益最大化。运营商的定价策略不能拍脑袋，必须考虑产消者的响应行为，否则电价定高了没人买、定低了运营亏损。这种双向耦合关系的求解，本质上就是求一个Stackelberg均衡。
1.2 为什么用IEEE33节点作为算例
IEEE33节点是配电系统研究里被用到极致的标准算例，中文论文里到处都是它的身影，原因其实很朴素：节点规模适中、辐射状拓扑清晰、参数公开容易获取，而且改造灵活。你可以在任意节点挂光伏、储能或电动汽车，也可以把部分节点定义为产消者，做竞价博弈研究时特别顺手。相比IEEE13节点那种偏北美风格的算例，IEEE33更接近国内城镇配电网的接线习惯，所以在中文EI期刊里出镜率极高。
使用IEEE33节点还有一个隐含优势：基准潮流结果公开透明，全网总有功损耗约为202.67 kW，最低节点电压在0.903左右（具体视参数版本略有差异）。复现代码的时候，你先把基础潮流跑出来和公开结果对比，能很快判断自己搭建的拓扑和参数是否正确。这一步校准没做好，后面所有博弈优化结果都是空中楼阁。
2. 核心理论与数学模型实现要点
2.1 主从博弈模型的上下层结构怎么设计
先明确一个核心点：主从博弈模型不是简单的双层优化，而是上下层决策变量互相嵌套。我们在Matlab代码里做了清晰的模块划分。
上层是配电网运营商（DSO），目标是自身收益最大化（也可以拓展为社会福利最大化或网损最小化）。决策变量通常是向产消者发布的购售电价λ_purchase（买电价格）和λ_sell（卖电价格），部分模型还会包含向上级电网购电的分时电价。DSO的收益主要来自低买高卖的差价、向常规负荷售电的收入，同时要扣除向上级电网购电的成本。
下层是各产消者（Prosumer），目标是自身净收益最大化。产消者内部通常包含光伏出力、储能充放电、灵活负荷，以及向电网购电/售电的功率，约束包括光伏出力上下限、储能SOC连续性约束、充放电功率约束、功率平衡约束等。当电价信号给出后，产消者通过决策自己的购售电量与储能充放策略来追求最大收益。
上下层的耦合点在于：DSO发布电价 → 产消者响应并反馈购售电量 → DSO根据反馈的总购售电量重新核算收益并更新电价 → 继续迭代循环，直到电价和购售电量都趋于稳定，即达到Stackelberg均衡。
这里的数学本质是：上层问题包含下层最优解函数（通过KKT条件或下层解替换），形成带均衡约束的数学规划（MPEC），求解难度远大于普通双层规划。我们复现时选择了KKT条件转换的方案，这个后面细讲。
2.2 KKT条件转换与线性化处理的几个关键点
把双层模型转化为单层模型，最常见的手段是把下层问题用KKT条件替换。具体分三步：
第一步，写出下层产消者优化问题的拉格朗日函数。注意这里要包含所有不等式约束的拉格朗日乘子（储能充放电限值、光伏出力上限等），等式约束的乘子单独区分。有个容易忽略的细节：下层目标函数是最大化收益，转换为最小化负收益再写KKT条件时，所有乘子的符号逻辑要重新理一遍，否则求出来根本不是均衡解。
第二步，列出KKT条件，包括梯度条件（stationarity）、原始可行性、对偶可行性（乘子非负）、互补松弛条件。这份KKT条件其实就是下层产消者“理性响应”的所有描述，把它加入上层问题，就形成了一个单层的MPEC问题。
第三步，处理非线性项。这里最让人头疼的是互补松弛条件（乘子乘以约束表达式等于0），这颗“雷”不拆掉，标准求解器根本不认。常规做法是引入大M法：把互补松弛条件拆成乘子≤M×二进制变量、约束表达式≤M×(1-二进制变量)这样的两对线性不等式。但大M的取值十分讲究，M太小会截断可行域，M太大会导致数值病态。我在IEEE33节点实例里反复调试后发现，M取10^3~10^4这个量级比较稳，前后迭代结果不会有数量级差异时，就判定取对了。
下层目标函数里的非线性项也需要线性化。比如储能运行成本系数乘以充放电功率的二次项，要处理成分段线性近似；购售电收入项中包含“电价×购/售电量”这个双变量乘积，当电价由上层决策、电量由下层决策时这是双线性项，需要用强对偶定理或McCormick线性化处理。这里我直接采用了强对偶定理，把下层问题对偶间隙清零后，用对偶目标函数替换双线性项，实测下来计算效率提升非常明显。
2.3 收敛判据与Stackelberg均衡判定
主从博弈求解完成之后，怎么判断结果是否达到了Stackelberg均衡？我们设了三个指标同时判定：上层电价变化量小于阈值（比如0.001元/kWh）、下层总购电量变化量小于阈值（0.01 MW）、上层目标函数值变化量小于阈值（1元）。三条标准必须同时满足才认定收敛，防止出现“电价稳定但电量还在跳”的假收敛。
另一个容易被忽视的细节：博弈模型可能出现多均衡解，不同初始电价下迭代结果可能落在不同均衡点。所以我在代码里设计了随机多起始点初始化，至少从三组不同的初始电价启动求解，如果最终收敛到同一均衡点，才认为结果可靠。这一点在写EI论文审稿回复时特别好用，审稿人问“怎么确保是全局均衡解而不是局部解”，这个多起点测试就是最有说服力的实验依据。
3. 基于IEEE33节点的Matlab落地实操
3.1 基础数据准备与拓扑搭建
IEEE33节点的基础数据在各种论文和代码仓库里都有，但不代表拿来就能直接用。我强烈建议自己在Matlab里把支路参数和负荷参数重新录入一遍，哪怕麻烦一点，至少出一份自己能核对的数据。我用的是经典版本：基准电压12.66 kV，基准功率10 MVA，根节点（节点1）通过上级电网供电，负荷数据包含有功和无功两部分。
录入数据时还有几个坑要提一下。第一，不同文章里的IEEE33支路阻抗可能略有差异（有些是标幺值，有些是欧姆值），务必确认采用的是同一套参数，否则对不上。第二，节点编号顺序要统一，我习惯用Matlab的图（graph）对象来存储拓扑，但要注意Matlab的graph对象默认自动排序节点，采用的是数字标号顺序，处理IEEE33中子节点号与母节点号不一致时容易出乱子，最好自己维护一个明确的节点父子关系表。第三，产消者节点要在基础数据上做改造，不能直接全部修改原负荷数据，否则基础潮流对比就失真了。
3.2 产消者节点改造与场景设计
从IEEE33节点里选取产消者节点是很有讲究的。不能盲目地选一堆节点，要体现“新型城镇配电系统”的特征。我在代码里选用了三个典型节点作为产消者：节点18接光伏+储能，节点22接光伏+电动汽车充电桩，节点25接储能+柔性负荷。三个节点分别代表居民型、商业型和混合型产消者，行为特性差异明显，竞价策略的博弈效果更突出。
每个产消者的内部参数也需要单独设定。光伏方面，接入容量在200kW到500kW之间，考虑辐照度曲线后生成24小时出力系数，这里可以直接用Beta分布模拟，也可以直接用典型日光伏曲线。储能方面，容量和功率按光伏配比来选，比如节点18的光伏300kW配储能200kWh、最大充放功率100kW，SOC范围设0.1到0.9，充放电效率取0.95。特别要注意，储能初始SOC不要设得太随意，我建议取0.5，这样全天调度后才会出现“白天充电、晚上放电”的典型套利行为，博弈过程会更精彩。
3.3 主从博弈求解代码的框架搭建
我们的求解框架采用Matlab + YALMIP工具箱 + CPLEX或Gurobi求解器。YALMIP负责把优化问题从代数层面建模，求解器负责真正跑数值计算。不少初学者第一步就卡在安装上：Matlab版本、YALMIP版本和求解器之间兼容性是个大坑，比如Gurobi官方现在主要支持较新版本Matlab，老版本Matlab要用旧版Gurobi。稳妥起见，我在项目文件里同时提供了Gurobi和CPLEX两种接口示例，哪个能用就用哪个。
代码整体分四个文件：主程序main.m、上层模型leader_model.m、下层模型follower_model.m、数据处理plot_results.m。主程序只负责数据加载、初始化电价、循环调用上下层，逻辑非常清晰。循环部分就是一个while循环，每次迭代中先调用follower_model求解下层问题，得到各产消者的购售电量，再把电量传给上层leader_model，求解DSO的最优电价。每回合记录电价和购售电量，达到收敛条件就跳出来，接着调用plot_results出图。
写下层模型时有个关键点：每个产消者对应一个小优化问题，可以循环单独求解，也可以把多个产消者合并成一个大的二次规划。前者好处是问题规模小、利于调试，后者好处是统一处理市场价格耦合关系。主从博弈迭代中每个产消者都只针对当前电价做独立响应，本身不存在产消者之间的直接耦合，所以我在复现代码里采用循环单独求解每个产消者的方案，速度也不会慢。
3.4 具体竞价迭代过程与结果输出
迭代过程其实很简单直白：假设早上8点，DSO给定初始购电价0.4元/kWh、售电价0.6元/kWh（初始值不要脱离实际太多，否则迭代次数会明显增加），节点18的光伏在此时出力大约只有额定容量的30%左右，储能SOC在中位，产消者下层模型经过优化后，得到“从电网购电120kW、储能充电50kW”这样的结果。DSO收到三个产消者汇总的购售电量后，发现总购电量大，说明当前电价偏乐观，于是适当降低购电价、拉大购销价差，下一轮产消者就会减少购电或增加售电。如此反复，通常30到50轮迭代之内完成收敛。
输出图像里最有价值的有四张图。
第一张是电价迭代收敛曲线，能直观看到电价从初值开始震荡后逐渐平稳的过程。 第二张是各产消者的24小时购售电功率柱状图，正数代表购电、负数代表售电，从这张图能判断产消者的竞价策略合理性。 第三张是储能SOC变化曲线，验证储能是否按照“午间光伏多时充电、晚间电价高时放电”的规律运行。 第四张是DSO收益构成饼图或柱状图，展示购销差价收益、向产消者卖电收益等各部分，用于论文的经济性分析。
我跑完的结果是这样的（可作为你复现时的参考基准）：均衡状态下，DSO制定的购电价全天均值约0.35元/kWh，售电价约0.63元/kWh，价差比初值更合理；三个产消者的总购电量从初始方案下降了18.6%，而DSO总收益比固定电价模式下提升了约9.4%。注意，不同参数设置的结果会有较大差异，拿这组数据当参照而不是当绝对标准，才是正确的复现心态。
4. 常见问题与排查技巧实录
4.1 YALMIP模型报错或不收敛
复现过程中最大的拦路虎就是模型怎么都不收敛，具体表现是迭代曲线发散、或YALMIP直接报“infeasible problem”。根据我的经验，这类问题列一个速查表特别好用：
| 症状 | 可能原因 | 解决办法 |
|---|---|---|
| 上层模型infessible | 大M取值偏小，切掉了可行域 | M从10^2开始递增测试，直到结果不再变化 |
| 下层模型求解失败 | 储能SOC约束和功率约束冲突 | 检查SOC更新方程是否所有时刻连贯，特别是初始时刻 |
| 迭代30轮以上仍发散 | 初始电价超出了合理区间 | 以当地电网实际电价水平1.2倍作为初始电价 |
| 结果为负数/功率溢出 | 缺少购售电互斥约束 | 增加购电量和售电量互斥的0-1变量约束 |
| CPLEX/Gurobi不识别模型 | 互补松弛条件未线性化 | 重新检查大M拆分的两对不等式是否完整 |
排查顺序建议从上到下进行，每查一步就单独跑一下上下层模型的可行性，别一上来就纠结算法和代码性能问题。
4.2 购售电互斥与储能约束导致的求解困难
产消者同一时刻既购电又售电，这在物理上完全没有意义，但优化模型如果缺少互斥约束，求解器会在某个时段给出一边买、一边卖的“套利”结果，既占用电网容量又虚增收益，看起来精美但实际荒唐。解决办法是引入二进制变量：购电状态b_p和售电状态b_s，加约束b_p+b_s≤1，同时购电量P_purchase≤b_p×M，售电量P_sell≤b_s×M。两个M取值原则上要大于系统最大传输功率，取500kW这个量级就可以。
储能约束是另一个容易出问题的地方。储能SOC更新方程里，如果充电效率η_c和放电效率η_d都小于1，连续时刻的SOC递推公式存在非线性耦合（充放电状态不同时效率取值不同）。很多复现代码为了简化直接把充放电效率统一为1，这在纯理论推导里可以，但算出来的储能行为会与工程实际偏差很大。我建议至少保留充放电效率，并引入充放电状态互斥变量，让储能模型更贴近锂电池实际运行特性，这样下游分析（比如电池寿命估计）才有基础。
4.3 加速求解的几个实战技巧
主从博弈迭代加MPEC求解，计算量确实不小。普通家用电脑跑IEEE33节点，如果代码没有做任何性能优化，一次完整求解可能要等十几分钟甚至更久。我实际采用并验证有效的加速手段主要有三个。
第一个是求解器参数调整。Gurobi里重点调MIPGap，从默认的1e-4放宽到1e-3，能节约大量分支定界时间，结果偏差在可接受范围内。CPLEX里对应参数是EpGap和MIPTolerances，也可以如法炮制。
第二个是热电联产式迭代初始值继承。每次迭代求解下层模型时，将上一轮的最优解作为当前轮的初始可行解传给求解器，可以大幅减少求解器内部的冷启动开销。这在YALMIP里用assign和initialguess实现，简单有效。
第三个是CPU并行化。如果有Matlab并行计算工具箱，可以使用parfor循环并行求解多个产消者的下层模型，得分很明显。原因在于各产消者的下层模型在给定电价后互相独立，本身就是天然的并行池任务。我在八核机器上实测，三个产消者并行求解比串行快2.5倍左右，迭代总时间从12分钟压缩到5分钟内。
5. 扩展经验与论文层面的一点建议
5.1 从复现代码到论文模型的拓展思路
复现EI论文不应仅停留在“跑通代码、得到几张图”的层面，更值得做的是思考模型的延展空间。目前这个模型里产消者下层的决策主体是单个用户，但真实城镇配电系统中往往存在多个产消者聚合而成的负荷聚合商，或者由第三方运营的虚拟电厂（VPP）。如果把这些真实市场主体引入主从博弈框架，上层的DSO定价、下层的聚合商竞价会形成更复杂的多领导者-多跟随者博弈结构，这时单纯KKT条件转换就不够用了，需要设计分布式求解算法或纳什议价方法。这是一个很有价值的后续研究方向。
另一个常见的扩展方向是引入碳交易机制。新型城镇配电系统的“新”不光体现在分布式资源上，碳约束下的经济运行同样是热点。可以在产消者下层目标函数中增加碳排放成本项，碳排放配额由政府或碳市场给定，DSO定价与碳价联动。这样既丰富了模型的现实意义，也增加了文章的创新点。我不止一次看到类似的论文修改策略在审稿阶段顺利通过，核心逻辑就是“博弈框架不变、增加一个符合政策导向的约束维度”。
再有一层扩展关乎鲁棒性。主从博弈复现时用的风光出力曲线是确定性的，但真实世界中光伏和负荷都有明显的不确定性。采用场景法（多个典型场景加权）或分布鲁棒优化（矩不确定性集合）来改造上层或下层模型，可以让研究成果从“理论算例”更进一步到“接近工程应用”，论文的创新性和实用价值都能同步提升。
5.2 代码组织与实验可重复性的一点心得
做这种复现型项目，最怕的就是过一阵子自己都看不懂自己写的代码。我这里分享一个自己的习惯：每一个模型文件开头必须写三段注释——模型类型（上层/下层/独立优化）、输入变量（变量名+单位+含义）、输出说明（哪些变量是决策变量、哪些是中间变量）。此外，每次跑完仿真后我会把关键参数、目标函数值和收敛状态存成一个run_results_时间戳.mat文件，这样后续做参数灵敏度分析时可以直接加载过往数据做对比，不用从头跑一遍算例。这些小习惯看起来繁琐，但当你需要补充实验数据去回应审稿人时，它们的价值就会被充分放大。
结尾
这个EI复现项目表面上是“跑一个主从博弈算例”，但实际上是一次完整的科研能力训练：从读懂论文里的数学模型，到用Matlab把数学公式变成可运行代码，再到通过迭代调试理解博弈求解的本质，每一步都有实实在在的坎要过。我在整个复现过程中体会最深的一点是，主从博弈这类模型的调试关键永远在“把模型分解到尽”：什么变量属于上层，什么变量属于下层，哪条约束连接上下层，能不能单独验证每一层的可行性。把这个问题问清楚了，剩下的更多是时间问题。
给准备入手这个方向的朋友一个非常具体的建议：不要一开始就追求完整复现出论文里的所有图表和曲线，先把IEEE33节点的基础潮流数据用代码校准好，再把只有一个产消者节点的简化博弈模型跑通，然后逐步增加产消者数量、逐步引入储能和光伏。这个由简到繁的过程不但能帮你规避大部分调参陷阱，还能让你对整个模型的逻辑链条形成肌肉记忆。等这套流程走完，你就不会觉得论文里的博弈模型高不可攀了。
