最近刚把一个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_{t1} 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 竞价机制的衔接规则主从博弈里的竞价机制具体指向的是内部市场怎么结算。复现时最常用的是统一出清价格机制CPSClearing Price Settlement所有成交的购电和售电都按同一时段统一的结算价格\lambda(t)来结。这就比PABPay-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 feasibilityAx \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数值不稳定过小会截断可行域。我最初用M10000CPLEX一直报数值警告部分互补松弛条件趋近于不满足后来把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效率高得多。另外储能充放电互斥的二进制变量尽量不要省省掉的话模型在低价时段会同时充电放电空转白白消耗能量结果出来你都不知道哪里错了。