这两年我一直在做配电网分布式交易和产消者竞价方向的仿真经常被同行问同一个问题你那个基于主从博弈的Matlab代码跑通了没有问的人多了我发现这个方向的代码资料其实不少但很多给出来也看不明白——要么只有模型推公式不给实现细节要么给了代码跑起来到处报错。这篇文章就把我整理的1029号模型也就是基于主从博弈的新型城镇配电系统产消者竞价模型从问题定义、数学建模、Matlab代码实现到算例复现完整讲一遍。参考资料以《基于主从博弈的新型城镇配电系统产消者竞价》为主我把读论文时容易卡住的点以及代码里对应的实现方式放到一起说。想用主从博弈做配电网多主体优化、光储协同、电动汽车有序充电这些方向的同学可以直接拿这份笔记去对照自己的代码。1. 主从博弈建模思路先把“谁主导、谁跟随”想清楚1.1 为什么城镇配电网里的产消者竞价天然适合主从博弈新型城镇配电网和传统配电网最大的区别就是用户侧不再只是被动消费负荷。分布式光伏、小型储能、电动汽车充电桩大量接入之后居民和商业用户可以在某个时段从电网买电也可以在另一个时段向电网卖电这就是“产消者”的由来。用户在博弈中不再只是价格接受者而是会成为价格的响应者但又不是完全的价格制定者。实际交易场景中配电系统运营商掌握电网拓扑、变压器容量、线路载流量这些公共信息也掌握外部电网购电成本和内部负荷预测因此它有资格制定内部交易价格而产消者数量众多、分散在不同节点每个产消者都只能根据运营商给出的购电价和售电价去优化自己的用电、发电和储能策略。这个决策顺序非常清晰上层先出价下层再响应下层响应反过来影响上层收益。这就是经典的Stackelberg主从博弈结构。主从博弈在这个场景里的优势我自己的体会是有三点。第一它建模逻辑贴合实际城镇配电网的运营商确实拥有主导定价权产消者确实只能在价格信号下做局部优化这种非对称决策关系用主从博弈来描述很自然。第二它能很好保护产消者的隐私产消者不需要把内部负荷、储能SOC等敏感数据全部交给运营商只需要在给定价格下上报购售电量适合多主体分布式交易的仿真需求。第三它比完全集中式优化更符合市场规律集中式建模要求所有主体把目标和约束全部放到一个优化问题里这在现实里根本做不到。和其他竞争模型相比差异也很明显。古诺模型里多个发电商同时决策产量大家都在同一层级没有先后之分伯川德模型里大家同时决策价格同样没有上下层关系而主从博弈强调领导者和跟随者的先后顺序正好匹配城镇配电系统运营商和产消者之间的真实关系。我常用下面这个表来向学生解释博弈模型决策主体决策顺序适用场景古诺模型多个同层主体同时决策产量发电商批发市场竞争伯川德模型多个同层主体同时决策价格售电公司价格竞争主从博弈领导者跟随者先定价再响应配电系统运营商与产消者竞价1.2 代码里的博弈层级与信息流设计拿到一个主从博弈题目最容易犯的错就是一上来就写代码。我建议先把手画一张信息流图哪怕只是在纸上草草画一下。1029这套代码里的信息流是这样的配电系统运营商在负荷预测和光伏预测的基础上先给出内部购电价格和售电价格各产消者收到价格后各自求解一个局部优化问题得到最优购电量、售电量以及储能充放电策略产消者把购售电量结果上报给运营商运营商汇总所有产消者的电量后检查全网功率平衡和线路潮流约束并根据电量偏差调整价格然后开始下一轮直到价格和电量都不再发生明显变化就认为达到了主从博弈均衡。代码的组织方式和这个信息流是一一对应的。主程序负责初始化参数、加载负荷和光伏数据、驱动迭代循环上层模型封装成upper_dso()函数输入是产消者上报的电量输出是新的价格下层模型封装成lower_prosumer()函数输入是价格输出是该产消者的购售电量和储能出力迭代结束后再调用一个绘图函数把价格曲线、电量曲线、储能SOC曲线画出来。这样模块化设计的好处是以后想改上层目标函数或者想给某个产消者增加约束只需要动对应函数不需要把整个代码推倒重来。2. 数学模型搭建成本、收益、约束一个都不能少2.1 上层DSO的定价优化模型在1029模型里上层配电系统运营商的目标函数不只是“卖电赚差价”这么简单。它要向外部电网购电来满足本地负荷需求也要向产消者售电赚取收入同时还要从产消者那里收购多余的光伏电量。目标函数可以写成在调度周期内最大化自身净收益max F_DSO Σ_t [ Σ_i (λ_s,i(t) * P_b,i(t) - λ_b,i(t) * P_s,i(t)) - c_grid(t) * P_grid(t) - c_loss * P_loss(t) ]其中λ_s,i(t)表示DSO向产消者i的售电价P_b,i(t)是产消者i在t时段的购电量λ_b,i(t)是DSO从产消者i的购电价P_s,i(t)是产消者i在t时段的售电量c_grid(t)是外部电网购电价格P_grid(t)是从上级电网购买的有功功率。注意这里的P_b,i和P_s,i并不是DSO直接决定的变量而是下层产消者优化后的响应结果这就把上下层耦合起来了。上层约束一般有三个层次第一是配电系统内部的功率平衡约束所有产消者购电量之和加上网损必须等于所有产消者售电量之和加上外部电网注入功率第二是线路容量约束每条馈线的功率潮流不能超过最大允许值第三是DSO报价约束价格不能超过政府允许的上限否则可能被监管机构叫停。我在代码里还额外加了变压器容量约束因为城镇配电网的变电站容量通常就是限制光伏接入的主要瓶颈。有一点需要特别提醒很多初学者会把上层目标函数写成“最大化全社会福利”这在集中式优化里没问题但在主从博弈里容易失去博弈意义。上层DSO既然处在领导者位置它就应该以自身利益最大化为目标而不是一上来就想着社会福利最大化。只有在论文想要对比“市场效率”和“系统安全”两个场景时才需要额外设计社会福利最大化的对照模型。2.2 下层产消者的响应模型每个产消者的角色可以看成是一个“小型电力公司”。它自己有光伏、储能和负荷既要从电网买电也可能向电网卖电。下层目标函数是最小化它在调度周期内的净成本min C_i Σ_t [ λ_s,i(t) * P_b,i(t) - λ_b,i(t) * P_s,i(t) c_pv,curtail * P_curtail,i(t) c_bat * (P_ch,i(t) P_dis,i(t)) ]这个式子的逻辑很直接从电网购电是一笔支出向电网售电是一笔收入光伏弃光会产生惩罚成本储能充放电会带来电池损耗成本。产消者需要在电价低的时候适当多买电给储能充电在电价高的时候释放储能甚至向电网卖电从而实现自己的利益最大化。实际代码里还要为每个产消者加上自己的功率平衡约束P_b,i(t) - P_s,i(t) P_load,i(t) P_ch,i(t) - P_dis,i(t) - P_pv,i(t) P_curtail,i(t)这个约束的含义是产消者从电网的净购电量等于自身负荷加上储能充电功率减去储能放电功率和光伏出力再加上弃光功率。如果右边是负的说明光伏出力大于本地需求产消者就会向电网卖电。储能约束是模型里最容易写错的部分。储能SOC递推方程是时域耦合的如果不加周期始末约束模型很容易把所有电量都放到最后一小时放完导致结果失真。我在代码里的写法是SOC_i(t1) SOC_i(t) η_ch * P_ch,i(t) * Δt / E_i - (P_dis,i(t) * Δt) / (η_dis * E_i) SOC_i(1) SOC_i(T1)其中SOC_i(t)是电池荷电状态η_ch是充电效率η_dis是放电效率E_i是储能容量。最后一个等式强制储能在一个调度周期结束后回到初始SOC这个约束对于日内调度模型非常关键否则优化结果在物理上不可执行。如果做的是多日滚动调度这个约束也可以改成每天SOC差不超限但总的来说SOC约束一定要仔细检查。产消者之间是相互独立的它们不会共享自己的负荷曲线和储能参数所以在代码里每个产消者对应一个独立的结构体里面存自己的参数。这也是主从博弈下“分布式决策”思想在编程上的体现。2.3 主从博弈均衡与KKT转化思路主从博弈的均衡点是上层DSO无法通过单方面改变价格来获得更多收益、下层产消者无法通过单方面改变交易策略来降低成本的状态。数学上如果上层问题是一个凹优化问题下层每个产消者问题是凸优化问题那么理论上Stackelberg均衡是存在的。实际代码里我们一般不去证明存在性而是通过两层迭代不停逼近看到价格和电量收敛到一个稳定点就认为达到了均衡。不过在Matlab里直接写两个独立优化问题来回迭代虽然直观但收敛速度慢而且有时候会在两个解之间振荡很难判断是否真正收敛。更规范的做法是把下层每个产消者的优化问题用KKT条件写成一组等价约束把这组约束放进上层问题里整个主从博弈就变成一个单层带均衡约束的优化问题也就是MPEC问题。这个思路在参考文献《基于主从博弈的新型城镇配电系统产消者竞价》中有明确的推导代码实现也是按这个逻辑来的。下一节我会给出具体怎么在Matlab中用Yalmip做这件事。3. Matlab代码实现从双层模型到可运行的仿真程序3.1 整体结构与数据对象我先说一下代码文件组织方式。很多科研代码最大的问题是所有内容都堆在主脚本里变量命名混乱后期想改参数只能全局搜索。1029代码我按功能拆成了下面几个文件1029_main.m % 主程序 data_params.m % 系统参数初始化 load_data.m % 负荷、光伏、电价数据加载 upper_dso.m % 上层DSO优化函数 lower_prosumer.m % 下层产消者优化函数 kkt_reform.m % 下层问题KKT重构 update_price.m % 迭代定价更新模块 plot_result.m % 结果可视化模块主程序就是一个循环加函数调用逻辑非常清楚。所有参数我统一放在一个结构体params里包括params.N产消者数量、params.T调度时段数、params.dt时间步长、params.PV每个产消者的光伏容量、params.Battery储能参数等。每个产消者单独用一个结构体数组prosumer(k)保存。用结构体而不是一堆散落的全局变量后续调参真的会轻松很多这是我从项目里踩出来的经验。初始化代码大概长这样%% 基础参数 params.N 8; % 产消者数量 params.T 24; % 调度时段数 params.dt 1; % 时间步长小时 %% 能耗数据 for k 1:params.N prosumer(k).load load_data(:, k) * 1e3; % 单位kW prosumer(k).pv pv_profile(:, k) * 1e3; prosumer(k).E 500; % 储能容量kWh prosumer(k).Pch_max 100; % 最大充电功率kW prosumer(k).Pdis_max 100; % 最大放电功率kW prosumer(k).soc0 0.5; % 初始SOC prosumer(k).eta_ch 0.95; prosumer(k).eta_dis 0.9; end如果需要用IEEE 33节点配网结构我会把网络数据也放到params.net结构体里用DistFlow方程做线性化潮流约束。3.2 用Yalmip定义上层和下层变量在Matlab里做优化建模我个人强烈建议用Yalmip工具箱而不是自己写求解器接口。Yalmip的语法非常接近数学表达式写起来不容易出错可以无缝切换使用Gurobi、Cplex、Mosek等商业求解器。对于双层模型Yalmip还提供了生成KKT条件的内置函数这是很多论文复现代码选择它的原因。下层产消者k的优化变量定义如下P_b sdpvar(1, T); % 购电功率 P_s sdpvar(1, T); % 售电功率 P_ch sdpvar(1, T); % 充电功率 P_dis sdpvar(1, T); % 放电功率 P_cur sdpvar(1, T); % 弃光功率 SOC sdpvar(1, T1); % 荷电状态然后写约束cons []; cons [cons, 0 P_b P_b_max]; cons [cons, 0 P_s P_s_max]; cons [cons, 0 P_ch Pch_max]; cons [cons, 0 P_dis Pdis_max]; cons [cons, 0 P_cur P_pv]; cons [cons, P_b - P_s P_load P_ch - P_dis - P_pv P_cur]; cons [cons, SOC(2:T1) SOC(1:T) eta_ch * P_ch * params.dt / E - P_dis * params.dt / (eta_dis * E)]; cons [cons, SOC_min SOC SOC_max, SOC(1) soc0, SOC(T1) soc0];这组约束看起来不复杂但已经把产消者的物理模型全部覆盖了。目标函数需要用到当前迭代轮次的价格lambda_s和lambda_b这两个价格在上层模型中其实是待求变量但在KKT重构中它们会同时出现在上层目标函数和下层KKT约束里形成一个大的优化问题。3.3 两种求解路线KKT单层化还是迭代定价在代码实现时我同时保留了两套求解逻辑。第一套是严谨的KKT单层化路线把下层每个产消者的问题写出来然后用拉格朗日函数对变量求导得到Stationarity条件、Primal feasibility、Dual feasibility以及Complementary Slackness条件再把这些条件全部塞进上层模型。Yalmip里可以直接用kkt()这个命令生成一个优化问题的KKT系统代码特别简洁% 先定义下层问题 objective_i sum(lambda_s .* P_b - lambda_b .* P_s c_bat * (P_ch P_dis) c_pv * P_cur); % 生成KKT约束 [kkt_cons, kkt_details] kkt(optimize(cons_i, objective_i), [], sdpsettings(solver,cplex));生成的kkt_cons就是下层问题的KKT条件集合。把这些kkt_cons叠加到上层问题里原来两层耦合的问题就变成了一个单层的、带互补约束的数学规划问题。互补约束是非线性的需要用大M法把它线性化成混合整数约束这部分我在3.4节展开。这种做法的优点是精度高可直接用商业求解器一次求解论文里写“Stackelberg均衡”也有底气。第二套是工程上更常见也更适合调试验证的迭代最优响应路线。思路很简单先给定一组初始价格依次求解所有产消者的下层问题拿到总购售电量上层根据这些电量重新更新价格然后反复迭代直到价格变化小于阈值。这个逻辑在代码里就是lambda_s 0.5 * ones(1, T); lambda_b 0.3 * ones(1, T); alpha 0.4; % 阻尼系数 for iter 1:200 P_b_sum zeros(1, T); P_s_sum zeros(1, T); for k 1:params.N [P_b, P_s] lower_prosumer(prosumer(k), lambda_s, lambda_b); P_b_sum P_b_sum P_b; P_s_sum P_s_sum P_s; end lambda_new update_price(lambda_s, lambda_b, P_b_sum, P_s_sum, params); if norm(lambda_new - [lambda_s; lambda_b]) 1e-4 break; end % 加阻尼防止振荡 lambda_s alpha * lambda_new(1,:) (1 - alpha) * lambda_s; lambda_b alpha * lambda_new(2,:) (1 - alpha) * lambda_b; end两种路线怎么选我建议如果你是验证模型可行性、还在调参数阶段先用迭代法跑通如果你准备写论文需要严格求解均衡或者想对比不同场景下的博弈结果就花时间把KKT重构路线跑通。我实际跑下来KKT路线在产消者数量少比如8个以内时特别快但如果产消者数量增加到几十个变量和约束数量爆炸MPEC求解时间会迅速拉长迭代法反而更稳健。3.4 大M线性化与互补松弛条件的处理互补松弛条件是KKT系统里最麻烦的部分。对于一个约束g(x) 0它的对偶变量是μ 0互补条件写成μ * g(x) 0。这个等式是双线性约束求解器不能直接处理。标准做法是引入一个二进制变量z和一个足够大的正数M把互补条件拆成下面两组约束μ M * z g(x) M * (1 - z)当z1时μ被压到0g(x)可以为正当z0时g(x)被压到0μ可以为正。这样就把非线性互补约束变成了混合整数线性约束。在Matlab里的写法是M 1000; % 大M值需要根据实际量级调整 z binvar(1, length(g)); cons [cons, mu M * z]; cons [cons, g M * (1 - z)];这个M的取值非常关键。取小了可能会剪掉真实的可行域导致求出来的解不是真正的均衡取大了会让求解器数值稳定性变差尤其在使用Cplex或Gurobi时常出现“数值病态”警告。我一般是这样处理的先不带互补约束用一个普通优化算一遍把对偶变量和约束函数值的大致范围打出来再根据量级设M通常会在100到10000之间。这个细节看着不起眼但能让求解成功率提高一大截。4. 实操复现算例参数与结果解读4.1 我采用的算例系统与参数复现1029模型时我没有直接套用完整IEEE 33节点原始数据而是采用了一个改进的城镇配电系统拓扑保留主要辐射状结构在8个节点接入不同类型的产消者。这样的好处是既能体现网络约束的影响又不会让MPEC问题的规模大到没法求解。如果直接用IEEE 33节点每个节点都设置产消者变量数量会非常大新手跑起来基本都会卡死。基础参数我列在下面方便你直接抄作业参数数值调度周期24小时步长1小时产消者数量8个含居民型、商业型、工业型光伏渗透率每个产消者5kW~200kW不等储能容量100kWh~1000kWh储充/放功率上限0.2C~0.5C荷电状态范围0.1~0.9购电价范围0.3~0.9元/kWh售电价范围0.2~0.6元/kWh外部电网购电价峰时0.85平时0.55谷时0.25元/kWh收敛精度价格变化小于1e-4负荷曲线和光伏曲线我尽量用真实典型日数据而不是拍脑袋填随机数组。居民负荷典型特征是早晚两个高峰商业负荷集中在白天工业负荷比较平稳。光伏出力用夏季典型日S型曲线中午12点到14点最高。把这些曲线画出来后我再把它们叠加到算例系统里。这样得到的仿真结果在提交论文或汇报时更有说服力。4.2 典型结果价格收敛与储能“低充高放”跑完代码后最先看的是价格收敛曲线。我的实测结果里初始价设定为售电价0.5元/kWh、购电价0.3元/kWh经过大概18次迭代后电价稳定在售电价0.58元/kWh左右购电价0.27元/kWh左右。价格不再大幅波动说明主从博弈已经到达均衡点。然后是产消者储能策略。结果非常清晰凌晨0点到5点外部购电价低DSO给出的内部购电价也低产消者倾向于从电网买电给储能充电午间光伏大发时居民型产消者的光伏出力超过了本地负荷它们开始向电网售电晚上18点到21点负荷高峰期间外部购电价最高DSO的售电价也贵储能开始放电满足本地负荷减少从电网购电。整个过程就是照着“低谷充电、高峰放电”的规律来走但具体充放电时刻会根据每个产消者的负荷曲线和光伏曲线错峰进行。我还做了固定电价模式的对照把内部交易电价锁定为常数这时候产消者没有主动调节的激励大量光伏电量在午间返送电网造成线路倒送压力到了傍晚负荷高峰又要从上级电网大量购电。而在主从博弈模式下DSO通过动态调整价格引导储能和负荷转移系统净负荷峰值下降了大约8%左右光伏消纳率也更高。这些数字虽然随算例参数不同会有变化但整体趋势是稳定的。4.3 参考资料和代码之间的对应关系如果你手上有那篇参考文献会发现论文里的公式很多初看很劝退。但对照代码来读会好很多论文里每个下标i和t在代码里基本都是矩阵的维度论文里的KKT条件推导对应代码里kkt()命令生成的约束集合论文里主从博弈迭代流程对应代码主程序的for iter 1:max_iter循环。我建议复现时先按“数据准备→下层模型→上层模型→KKT推导→迭代求解”的顺序自己画一个映射表把论文公式编号和代码文件位置对应起来后面调试会非常高效。参考资料的价值不是给你背诵公式而是告诉你为什么这么建模代码只是把“为什么”落地成“怎么算”而已。5. 常见问题与调试经验这些坑我都踩过5.1 迭代不收敛或者价格振荡怎么办这是主从博弈代码最常遇到的问题。现象是价格在某个区间来回震荡不收敛到稳定值。我试过几种解决手段最有效的是在价格更新时加入阻尼系数也就是把新价格和旧价格做加权平均lambda_new alpha * lambda_candidate (1 - alpha) * lambda_old;alpha通常取0.3到0.6之间。如果还振荡可以把alpha调小到0.1或0.2。另一个办法是检查价格更新的步长不要一次变太多可以在update_price()函数里限制单次价格变化不超过0.02元/kWh这样能显著提升稳定性。还有一种情况是模型本身存在多个均衡这往往出现在产消者类型差别过大、储能容量差异悬殊的时候。这时可以尝试更换初始价格多跑几组初始点看是否收敛到同一个均衡如果收敛到不同均衡说明模型有多均衡论文里要如实讨论这也是很有价值的研究点。5.2 KKT线性化后求解器报错或结果异常用KKT路线求解时如果Cplex或Gurobi报“infeasible”或“unbounded”我首先怀疑的是大M参数设置。M太小会强行约束变量域导致原本可行的均衡被切掉M太大又会导致数值病态。解决办法是分两步走先用一个小规模松弛模型得到对偶变量的数量级再回到原模型里把M设成该数量级的10到100倍。如果还不行检查一下互补约束拆分的公式是不是把不等号方向写反了。另一个常见错误是漏了KKT里的dual feasibility条件即所有对偶变量必须非负忘记加这个约束会让求解器跑出一个看似合理、实际不符合KKT条件的解。5.3 Yalmip版本和求解器兼容性Yalmip版本更新很快不同版本对kkt()的支持程度不同。我最早用老版本跑kkt()命令结果直接报错“kkt only works with ...”后来升级到较新的Yalmip版本才解决。建议你直接用Matlab 2021a以上版本配合最新Yalmip求解器优先选Gurobi或Cplex因为它们对混合整数规划支持最好而且线性化后的MPEC问题本质是MILP或MIQP。如果没有商业求解器许可证可以在学术许可下申请Gurobi免费license比Cplex更容易拿到。实在不行用Yalmip自带的solve加默认求解器也能跑通小规模算例但收敛速度会慢不少不推荐做论文主算例。5.4 DistFlow潮流约束过紧导致无解如果加了网络潮流约束后模型突然无解我一般先检查线路容量设置是否合理。城镇配电网的线路载流量有限高光伏渗透时会反向潮流某些支路可能接近或超过上限。这时不要急着把容量上限放大而是应该调整储能的安装位置和容量配置或者增加弃光机制。代码里可以用DistFlow的二阶锥形式也可以用文献里常见的线性化DistFlow。线性化形式简单、计算快但会忽略网损对电压的影响在长线路或末端节点误差较大。如果只是复现论文线性化DistFlow一般够用如果做工程评估建议用二阶锥松弛Yalmip加Cplex处理起来也不算难。6. 扩展方向与实用建议6.1 在主从博弈框架里增加不确定性城镇配电系统的负荷和光伏出力都有很强的随机性确定性主从博弈模型虽然结构清晰但工程实用性有限。扩展方向之一是把下层产消者的优化改成两阶段随机规划或者引入分布鲁棒优化让产消者根据最坏情况下的成本做决策。上层DSO也可以考虑用条件风险价值来约束价格波动风险。代码层面不需要推翻重来只需要把下层目标函数里的成本项替换成风险度量形式KKT推导会复杂一些但整体框架还是能沿用。6.2 加入电动汽车有序充电新型城镇配电系统里电动汽车的占比越来越高产消者类型也可以扩展成“光储充”一体的智慧能源站。电动汽车充电负荷和普通家用负荷不一样它具有很强的时空转移特性。把电动汽车接入产消者模型时需要在功率平衡约束里加入充电桩功率同时新增“车辆接入时段”和“离开时SOC需求”约束。这部分我建议不要一开始就在每个产消者里都加而是先单独建一个“含EV产消者”模块跑通后再接入主循环避免修改一处影响全局。6.3 多主多从博弈和分布式求解如果未来要考虑多个配电台区之间的竞争或者多虚拟电厂之间的博弈单领导者的主从博弈就不够用了需要扩展成多领导者Stackelberg博弈。这类问题求解难度会指数级上升一种可行方案是把上层多领导者之间的竞争用广义纳什均衡来描述再用分布式算法分解求解。目前Matlab代码可以实现小规模场景但大规模场景还是建议转到Python或Julia环境配合商业求解器。这不是说Matlab不行而是模型规模到了一定程度语言生态和并行计算能力会有明显差别。6.4 最后分享一个实用小技巧调试这类模型我强烈建议从小算例开始。先用“2个产消者6个时段”的简化模型跑通代码确认博弈过程能收敛、KKT条件没写错再逐步扩大产消者数量和时段数。千万不要一开始就上“33节点24小时8个产消者”的完整算例一旦结果不对你很难判断是模型参数的问题还是代码逻辑的问题。还有一个很多人会忽略的点在迭代过程中把每一轮的价格、电量和目标函数值都打印出来哪怕只是用disp()打印到命令行也能帮你快速定位到底是上层更新出了问题还是下层优化不响应。我调试时几乎都会盯着这些中间量比最后只看一堆曲线图高效得多。
