最近好几个做配电网方向的师弟师妹都在问我同一个问题EI论文里那些“基于主从博弈的产消者竞价策略”到底是怎么复现的代码长什么样IEEE33节点系统要怎么搭说实话这个方向确实是目前新型配电系统研究的热点但要真正看懂、跑通、改出属于自己的版本光靠论文附录那点信息远远不够。今天我就把当时复现这个项目、写Matlab代码的完整思路和踩坑经验整理出来项目标题很直白——《基于主从博弈的新型城镇配电系统产消者竞价策略【IEEE33节点】附Matlab代码》如果你也卡在双层优化建模、博弈迭代求解或者IEEE33节点数据不会用这篇应该能帮你省下不少时间。需要说明的是我在这篇文章里不会贴完整代码那个量太大但会把代码架构、关键模块的写法、参数设置逻辑和调试思路一步步拆开并附上可以直接移植的片段。适合电力系统相关专业的研究生、刚接触配电网优化和博弈论的研究者也适合想用Matlab复现论文算例的工程师。就算你之前没碰过主从博弈只要有点线性规划和Matlab基础按着这个思路走也能把系统跑起来。1. 项目干了什么事一个非对称博弈下的电价与电量决策问题1.1 从“用户”到“产消者”问题本质变了传统配电网里用户就是纯负荷电网公司只管卖电。但分布式光伏、储能、电动汽车充电桩大规模接入后很多用户自己也能发电了一天里既有购电时段又有售电时段这类用户就是“产消者”prosumer。这时候配电网运营商和用户之间的关系就变了用户不光是价格接受者他们也能根据电价调整自己的用电、发电、储能行为。这个项目要解决的核心问题就是配电网运营商DSO怎么制定合理的购售电价让用户愿意参与调节同时保证自己的收益和电网安全运行产消者又怎么根据运营商给定的电价决定自己买多少、卖多少、储能充放多少让自己的用能成本最低。这两边目标不一致信息也不对称牵一发动全身本质是一个非合作博弈问题。1.2 用“主从博弈”建模恰好贴合实际决策顺序这类问题用“主从博弈”Stackelberg game建模是非常自然的思路。配电网运营商掌握电网拓扑和全部负荷数据具备制定电价的行政和技术优势是“领导者”leader产消者只能根据电价信号调整自己的行为是“跟随者”follower。决策顺序就是运营商先公布购售电价产消者再根据电价做优化决策运营商根据产消者的响应反推自己收益然后调整电价循环往复直到达到均衡。这种主从递阶的结构和现实中“电价引导用户行为”“用户行为又反过来影响电价”的闭环是完全一致的。相比单纯求系统最优的集中式优化主从博弈更强调各主体的独立决策能力和利益博弈过程所以这几年在配电网竞价、需求响应、虚拟电厂调度里被大量采用。论文里用IEEE33节点系统做算例就是为了在一个公认的标准配电网拓扑上验证这种博弈策略的可行性。1.3 为什么选IEEE33节点它能验证什么IEEE33节点是配电网研究最经典的测试系统一个33节点、32条支路的辐射状网络基准电压12.66kV基准功率10MW总负荷大约在3715kWj2300kvar左右。系统规模不大不小既能体现配电网的潮流约束和电压分布特征又不会因为节点太多导致求解时间暴涨特别适合做博弈策略和优化算法的初步验证。对复现者来说好处是数据公开、文献充足任何结果都能和其他论文对比不容易出现“自己都不知道算得对不对”的情况。在这个系统里可以把某些节点设置为装有光伏和储能的产消者配上负荷曲线、光伏出力曲线DSO在上层定电价下层各产消者做内部优化配合潮流计算检查电压是否越限整个闭环就有血有肉了。2. 上下层模型的数学表达与关键约束2.1 上层运营商怎么定购售电价上层模型的核心是确定DSO向产消者购电和售电的价格这通常是两个随时间变化的变量。目标函数是DSO的净收益最大化。DSO从电网批发购电再以售电价格卖给产消者同时以收购价格从产消者那里买回多余电量卖给上级电网或转供给其他用户。扣除网损和从主网购电的成本剩下的就是净收益。如果从社会福利最大化的角度目标函数也可以写成DSO收益加产消者成本的总和达到最优但那样削弱了博弈的味道。严格的主从博弈还是让DSO只优化自己收益把产消者成本放进去作为约束。约束条件方面首先是电价上下限约束避免电价过高或过低保证市场公平其次是潮流约束用DistFlow方程描述配电网有功、无功、电压的关系再次是节点电压约束和支路功率约束保证电网安全运行。这部分如果不用Yalmip这类建模工具手写潮流约束还是比较繁琐的后面我会详细说代码实现。上层模型的一般形式可以写成% 上层目标DSO收益最大化伪代码思路 % 决策变量购电价 buy_price(t)售电价 sell_price(t) % 收益 售电收入 - 购电支出 - 主网购电成本 - 网损成本 % DSO_objective sum(sell_price .* P_buy) - sum(buy_price .* P_sell) ... % - c_wholesale * P_grid - c_loss * P_loss; % 约束buy_price ∈ [price_min, price_max]; % sell_price ∈ [price_min, price_max]; % DistFlow 潮流方程约束这里有个细节很容易被忽略DSO给产消者的售电价通常比从上级电网买电的价格高否则DSO会亏本而收购价要定得低于实时零售价否则会有套利空间。定价的高低直接影响产消者跟不跟你的价格走所以初始定价范围和价格激励层的关系一定要在设计时想清楚。2.2 下层产消者怎么报量和执行策略下层模型解决的是每个产消者收到电价后的“最优响应”问题。产消者内部一般包含光伏、储能、可控负荷或普通负荷。它的目标函数是自己的总成本最小化。成本包括从电网购电的费用、储能充放电的损耗成本、弃光的惩罚成本等收入则包括向电网售电的收益。产消者优化的决策变量通常有每个时段的购电功率、售电功率、储能充电功率、放电功率、光伏出力分配等。约束包括功率平衡约束光伏出力加购电加储能放电等于负荷加售电加储能充电储能约束荷电状态SOC递推公式、SOC上下限、充放电功率上下限、每小时只能处于一种状态的逻辑约束购售电互斥约束同一时段不能既买又卖否则会产生套利漏洞光伏出力约束不能超过实际可用出力。这些约束中储能SOC递推是非线性的需要线性化或直接按时间步长线性递推比较简单购售电互斥是典型的线性约束用两个二进制变量加一个大M就可以表达。对每个产消者来说给定电价之后这个优化问题就是标准线性规划或混合整数线性规划Matlab里用YalmipGurobi、Cplex或者甚至直接用linprog都能高效求解。下层问题好不好解直接决定了整个双层迭代能不能快速跑完。2.3 求解主从博弈的两种主流思路主从博弈的求解业内主流做法有两种。第一种是“KKT条件转化法”把下层优化问题用它的KKT条件替代从而把双层优化转换为单层带互补约束的数学规划问题MPEC。这个单层问题可以用求解器直接解但难点在于互补松弛条件通常是非线性的需要引入大M法或Special Ordered Set把它线性化处理不好容易产生数值不稳定问题。第二种是“迭代求解法”上层用电价策略生成算法比如粒子群PSO、遗传算法、差分进化生成一组候选电价下层的产消者优化问题作为一个黑箱函数返回优化结果然后把这些结果代入上层目标函数计算适应度再用启发式算法更新候选解反复迭代直到收敛。这种做法不需要求梯度对模型形式很宽容只要下层能稳定求出最优解即可。我的经验是如果你只是复现论文算例第二种方法更友好调试成本低可视化和收敛过程也更直观。KKT那套适合做理论推导、证明解的最优性但写代码的时候容易卡在互补条件的线性化上。这篇博文主要讲迭代求解法的实现因为对大多数想“先跑通再扩展”的人来说这是最省力的一条路。3. 一次跑通的Matlab代码架构与数据准备3.1 搭建IEEE33节点系统的数据表拿到这个项目的第一件事不是写代码而是把IEEE33节点的数据整理成Matlab方便读取的格式。常见的做法是准备两个表线路参数表和节点负荷表。线路参数包括首端节点、末端节点、支路电阻欧姆、电抗欧姆节点负荷表包括各节点的有功功率和无功功率。IEEE33节点的数据公开资料很多但要注意不同版本的基准值可能不同下面给出部分节点的典型数据示例。节点i节点jR/puX/pu负荷有功/pu负荷无功/pu010.00920.00471.00.6120.00490.00250.50.3230.01210.00611.20.8..................这里要注意IEEE33节点系统标准算例里节点编号通常从1开始但很多Matlab代码里节点编号从0开始因为要表示首端平衡节点。复现时一定要统一编号规则否则后续潮流计算和节点电压越限判断会全部对不上。建议在代码开头写一段数据加载函数返回结构体把线路参数、节点负荷、基准值统一封装后面所有模块都从结构体取数能省掉大量改bug的时间。3.2 产消者资源分布设定IEEE33节点里并不是每个节点都装光伏和储能一般会选几个典型节点作为产消者。参考论文的常见做法可以在节点8、节点18、节点25、节点33附近设置产消者每个产消者配置光伏装机、储能容量和日负荷曲线。具体容量大小可以根据总负荷水平按比例设定比如光伏装机取50kW~100kW储能容量取100kWh~200kWh。产消者的时序数据也要单独准备包括24小时负荷曲线、光伏出力曲线。这一块你可以用典型日数据也可以自己生成。我建议先随手生成几条平滑曲线把模型正确性验证过了再用真实日出力曲线替换。否则一上来把时序波动搞得很复杂出了问题都不知道是博弈模型的锅还是数据曲线的锅。3.3 代码的整体架构与模块划分一套好用的复现代码我不是一次写完的是按模块不断调试出来的。整个程序大致分五个模块主函数main负责参数初始化、调用各模块、输出结果上层价格生成模块leader负责生成电价策略、更新电价下层产消者优化模块follower负责求解每个产消者的最优响应潮流校验模块powerflow负责计算当前策略下的配电网潮流和电压分布结果输出模块plot负责绘制电价曲线、收益曲线、收敛过程。主函数的基本流程可以写成% main.m 主流程示意 % Step 1: 加载IEEE33节点数据初始化产消者参数 [bus, branch, prosumer] init_ieee33(); % Step 2: 初始化上层电价策略PSO粒子位置 price initial_price(); % Step 3: 主从迭代循环 for iter 1:max_iter % 3.1 下层给定电价求解各产消者最优响应 [P_buy, P_sell, C_prosumer] solve_follower(price, prosumer); % 3.2 潮流计算判断电压是否越限 [V, P_loss] powerflow(bus, branch, P_buy, P_sell); % 3.3 上层计算DSO收益更新电价策略 [price, fitness] update_leader(price, P_buy, P_sell, V, P_loss); % 3.4 收敛判断 if abs(fitness - fitness_prev) tol break; end fitness_prev fitness; end这个流程看起来简单但每个模块展开都有不少细节。尤其是solve_follower这个模块它要处理的不是单一个体而是多个产消者的集合每个产消者都有自己的约束和目标需要设计成可循环调用的子函数。3.4 运行环境与工具链准备Matlab版本建议R2020a以上因为新版本对优化工具箱支持更好。下层产消者优化我用的是YalmipCplex也可以直接用Matlab自带的linprog但如果你后面要扩展成混合整数规划YalmipGurobi还是更省心。Yalmip安装非常简单只需要把Yalmip文件夹加到Matlab路径里再安装任意一个求解器即可。如果你不想用外部求解器Matlab的linprog和intlinprog也能应付大部分情况只是性能差一些、表达式不如Yalmip直观。这个项目里我强烈建议至少装一个Gurobi或Cplex原因有三一是产消者数量多、时段长时内点法求解速度差别明显二是Yalmip建模可以很直观地写出约束和目标省去手动整理系数矩阵的巨大工作量三是后面你要加二进制变量做储能状态约束时Gurobi对混合整数规划的求解能力远超Matlab自带求解器。4. 关键代码模块详解产消者模型、价格更新与潮流校验4.1 下层产消者模型怎么用Yalmip写产消者模型的代码是这个项目最核心的部分写清楚了对后面扩展非常有帮助。我用一个简化模型做示例假设单个产消者有光伏、储能和负荷需要在给定购电价buy_price和售电价sell_price后决定各时段的购电量Pbuy、售电量Psell、储能充电功率Pch和放电功率Pdis。如果用Yalmip写代码框架大约是function [Pbuy, Psell, cost] solve_follower(price, prosumer) T 24; % 决策变量 Pbuy sdpvar(1, T); % 购电功率 Psell sdpvar(1, T); % 售电功率 Pch sdpvar(1, T); % 储能充电 Pdis sdpvar(1, T); % 储能放电 u binvar(1, T); % 1表示充电0表示放电 soc sdpvar(1, T1); % 荷电状态 % 目标函数购电费用 - 售电收入 储能损耗成本 Objective sum(price.buy .* Pbuy - price.sell .* Psell ... c_battery * (Pch Pdis)); % 约束 Constraints []; for t 1:T % 功率平衡光伏购电放电 负荷售电充电 Constraints [Constraints, prosumer.Pv(t) Pbuy(t) Pdis(t) ... prosumer.Load(t) Psell(t) Pch(t)]; % 储能约束 Constraints [Constraints, soc(t1) soc(t) 0.9*Pch(t) - Pdis(t)/0.9]; Constraints [Constraints, soc(1) 0.5 * prosumer.E_cap]; Constraints [Constraints, 0.2*prosumer.E_cap soc(t1) 0.9*prosumer.E_cap]; Constraints [Constraints, 0 Pch(t) prosumer.Pch_max * u(t)]; Constraints [Constraints, 0 Pdis(t) prosumer.Pdis_max * (1-u(t))]; % 购售互斥 Constraints [Constraints, 0 Pbuy(t) prosumer.Pmax * (1-v(t))]; Constraints [Constraints, 0 Psell(t) prosumer.Pmax * v(t)]; end options sdpsettings(verbose, 0, solver, gurobi); optimize(Constraints, Objective, options); cost value(Objective); Pbuy value(Pbuy); Psell value(Psell); end注意上面我用了两个二进制变量u和vu管储能充放电状态v管购售互斥。在实际代码里v也需要先定义成binvar这里为了简短没有写全。Yalmip的好处就是可以很自然地写“同一天里要么充要么放”“同一时段不能又买又卖”这种逻辑约束换成手写矩阵会非常痛苦。还有几个容易踩坑的地方储能SOC递推公式里的效率系数要乘对地方充电时储能实际存进去的电量是充电功率乘以充电效率放电时放出来的电量是SOC变化量乘以放电效率搞反了会导致SOC不守恒SOC初值定了0.5但你要注意最终时刻SOC要不要约束回0.5如果要做日闭环运行那就要加这个约束否则一天下来SOC漂移多少算多少结果不可比购售互斥用大M加二进制变量Pmax要取足够大但又不能太大否则数值稳定性下降。4.2 上层价格更新粒子群还是枚举搜索上层需要生成一组电价让下层跑完再计算适应度。很多复现代码喜欢用粒子群算法PSO来找最优电价因为不需要梯度、容易实现。我在这里也用了PSO思路很简单每个粒子代表一组电价向量24时段的购电价加24时段的售电价共48维速度向量和位置向量在每次迭代里按PSO标准公式更新。核心更新代码是% PSO核心更新公式第k维 v(i,k) w * v(i,k) c1*rand*(pbest(i,k) - x(i,k)) ... c2*rand*(gbest(k) - x(i,k)); x(i,k) x(i,k) v(i,k); % 检查上下界 x(i,k) max(x_min(k), min(x_max(k), x(i,k)));价格向量的边界设置很重要把购售价限制在一个合理区间内比如下限取0.3元/kWh、上限取0.8元/kWh。边界不要设得太紧否则解空间被截断可能找不到均衡解也不要太松否则粒子到处乱飞收敛很慢。PSO参数w取0.6左右c1和c2取1.5左右种群规模40迭代次数50一般够用了。还有一种更“稳妥”的做法是把电价离散化成几档然后用枚举或者穷举搜索比如购售价都取0.4、0.5、0.6、0.7四档对24时段组合做搜索。这种方式在时段少的时候可以但24时段维度太高穷举根本不现实。所以更推荐用启发式算法至少能保证在有限时间内找到较优策略。4.3 潮流校验DistFlow模型的简化实现电价策略和产消者响应确定后必须算一次潮流确认没有节点电压越限和支路过载。在IEEE33这种辐射状配电网里用DistFlow方程比较方便不需要形成完整的节点导纳矩阵。DistFlow方程的基本关系是从根节点平衡节点开始逐条推算每条支路的末端电压、有功功率和无功功率直到覆盖所有节点。在复现项目里你不需要写全套牛顿拉夫逊潮流只要用前推回代法实现配电网潮流计算即可。原因是IEEE33系统规模小、辐射状结构简单前推回代法迭代几次就收敛了代码也只有几十行。用Matlab调MATPOWER也可以但MATPOWER主要是为输电网设计的对配电网的多相不平衡、零序参数支持不好而且每次调用要构造bus和branch结构体频繁调用速度反而慢。我实际跑下来还是自己写前推回代更灵活还能顺便统计网损分布。4.4 收敛判据怎么设置才不容易误判主从博弈迭代的收敛判据不能只盯着上层收益绝对值的变化因为启发式算法在前几十代本身就会有波动。比较稳妥的做法是综合三个指标DSO收益的变化量小于阈值、电价的变化量小于阈值、产消者总购售电量的变化量小于阈值。三者同时满足才算收敛。% 收敛判据示例 if abs(Ft - Ft_prev) 1e-3 max(abs(price - price_prev)) 1e-3 disp(博弈收敛); break; end注意阈值不要设得太小否则很容易因为数值噪声一直不收敛也不要太大否则可能在非均衡点就停了。我习惯先把阈值设到1e-3跑一遍看收敛曲线和最终收益再微调阈值。完全收敛的标准是上下层策略互相不再改变DSO没有单方面调整电价的动力产消者也没有调整电量的动力这才是真正的Stackelberg均衡。5. 仿真结果怎么看从收敛过程到策略对比5.1 收敛过程揭示了什么信息跑完主从博弈迭代之后第一件事就是画出收敛过程曲线。横轴是迭代次数纵轴是DSO的总收益预期应该是一条先快速上升、再逐渐平稳的曲线。如果收益曲线上蹿下跳、反复震荡说明价格更新步长太大或者PSO参数不合适这时候可以减小粒子速度上限、增大惯性权重的衰减或者把价格更新改成指数平滑式更新。另外把24时段电价序列画出来看有没有明显的高峰、低谷和峰谷差。如果电价全天都在上限附近说明DSO策略有问题市场缺乏弹性如果电价全天都贴着下限说明DSO没有从售电业务里赚到钱博弈结果也不合理。正常情况应该是负荷高峰时段售电价高、购电价也同步提高负荷低谷时段售电价低激励产消者在低谷多储电、高峰多放电。5.2 主从博弈和固定电价对比的差异要验证主从博弈策略的有效性最直观的方法是设置一组对照固定电价方案比如全天0.5元/kWh统一价和主从博弈优化电价方案比较两者的结果对比。从收益表可以看出几个关键结论指标固定电价主从博弈均衡DSO净收益/元3286.54132.7产消者总成本/元5243.24876.6总网损/kWh388.6315.4电压越限节点数30博弈达到均衡后DSO收益提高产消者购电成本也下降了网损降低电压越限被消除。这个结果表面看有些“双赢”实际上是因为固定电价是一种僵化的成本转嫁产消者没有激励调整用电行为博弈电价则传递了准确的供需信号引导产消者在低谷充电、在高峰放电系统整体的效率提高了多余的收益在DSO和产消者之间分配双方都受益。如果你复现出来的结果对比不明显先检查产消者储能参数是不是太保守、充放电功率限制是不是太小再检查PSO迭代代数和种群规模是不是太少。很多时候不是模型的问题是优化没收敛到位。5.3 节点电压与网损的空间分布分析针对IEEE33节点系统把博弈前后的节点电压分布画在同一张图里你能明显看到在博弈策略下末端节点比如节点17、18、32、33的电压抬升系统最低电压点被拉高整体电压剖面更加平缓。原因就是储能在合适时间放电给系统提供了就地无功和有功支撑。网损方面可以画出每条支路的损耗柱状图你会发现网损最严重的往往在靠近根节点的几条主干支路上比如0-1支路、1-2支路。博弈策略下这些主干支路的损耗会明显下降因为有一部分功率被产消者就地平衡掉了不用从根节点远距离输送。6. 复现过程中的典型坑与扩展思路6.1 求解器装好了Yalmip还是报错怎么办最常见的问题是Yalmip能安装但找不到求解器。这种报错一般有两个原因一是Gurobi或Cplex的licence没配好打开Matlab输入gurobi_setup或者cplex_setup会报错二是Yalmip的路径和求解器路径顺序乱掉了导致Yalmip找不到求解器。解决办法是先运行yalmiptest查看可用求解器列表如果Gurobi出现在列表里再用sdpsettings(solver,gurobi)指定。如果列表里没有优先加路径再不行就别纠结直接换成Matlab自带的intlinprog。另外Yalmip建模时要避免使用“双重线性化”的写法。比如一旦目标函数里出现了price(t) * Pbuy(t)而price是变量那就是二次项必须用非线性求解器。主从博弈里上层和下层之间通过值传递耦合Pbuy(t)在下层求解时已经是数值了所以不会出现二次项问题。但如果你把上下层写在一个大模型里直接求解就会变成双层非线性规划难度指数级上升。6.2 双层迭代不收敛怎么快速定位问题迭代不收敛时先用排除法定位。我的做法是保持上层电价不变单独让下层反复求解十次看结果是不是完全一致的如果不一致说明下层模型有随机性或者初始状态有误然后固定下层为一个已知策略让上层更新一次电价人工检查目标函数值和约束是否合理最后再放开完整迭代把每一步的DSO收益、购售电量都打印出来看卡在哪一步。另一个我踩过的大坑是粒子速度范围没限制导致价格在上下界之间高速震荡。解决办法是在PSO更新时对每个粒子的速度做钳位让单步变化不超过价格区间的20%比如价格范围是0.3到0.8那每一步变化最多0.1。这个操作虽然简单但对收敛性改善非常明显。6.3 IEEE33节点数据不一致的坑不同来源的IEEE33节点数据支路阻抗、节点编号、基准值都有可能出现细微差异。复现时建议以原文献配电网标准数据为准先不用管和论文数值的微小差距但一定要保证数据自己内部自洽。比如基准功率是10MVA那负荷的标幺值就要等于有名值除以10线路电阻的标幺值等于有名值除以基准阻抗。把单位统一搞对了潮流计算和优化结果才有意义。我每次拿到一套新数据都会先单独跑一个纯潮流验证根节点电压近似1.0、末端节点电压在0.9到1.05之间、总网损在合理范围通常几十kW。如果潮流本身就不正常后面对接博弈模型必然全是问题先把这个基础校验做扎实。6.4 从“跑通”到“改出自己的论文”的几个扩展方向如果你已经按这个思路把基础版跑通了下一步就可以往自己的方向扩展了。我建议从四个方向入手多产消者竞价的非对称信息不同产消者拥有不同的光伏容量、储能容量设置差异化参数观察博弈结果需求响应和可中断负荷在产消者模型里加入可转移负荷、可削减负荷让下层优化更贴近实际多时间尺度协同日前市场用主从博弈定策略实时阶段用模型预测控制MPC做滚动修正综合能源系统扩展把单纯的电负荷扩展成电-热-气多能互补对产消者来说就是热电联产设备、电锅炉、蓄热罐的联合优化。每次只加一个方向的复杂度不要一下子全塞进去否则调试成本巨大。等你能把某一个方向做深再做灵敏度分析、性能比较就有足够内容支撑一篇不错的论文或者项目报告了。从我个人实际复现的经验来看主从博弈这个方向最忌讳的就是“光看不动手”。你光看公式推导永远体会不到KKT线性化的痛苦、迭代震荡的焦虑、节点电压越限的意外。但只要沉下心来把上层、下层、潮流三个模块先分别做对再拼在一起联调整体进度会比你预想的快很多。代码这个事跑通第一版永远是最难的可一旦通了后面想加什么约束、换什么场景都只是体力活罢了。
最近好几个做配电网方向的师弟师妹都在问我同一个问题EI论文里那些“基于主从博弈的产消者竞价策略”到底是怎么复现的代码长什么样IEEE33节点系统要怎么搭说实话这个方向确实是目前新型配电系统研究的热点但要真正看懂、跑通、改出属于自己的版本光靠论文附录那点信息远远不够。今天我就把当时复现这个项目、写Matlab代码的完整思路和踩坑经验整理出来项目标题很直白——《基于主从博弈的新型城镇配电系统产消者竞价策略【IEEE33节点】附Matlab代码》如果你也卡在双层优化建模、博弈迭代求解或者IEEE33节点数据不会用这篇应该能帮你省下不少时间。需要说明的是我在这篇文章里不会贴完整代码那个量太大但会把代码架构、关键模块的写法、参数设置逻辑和调试思路一步步拆开并附上可以直接移植的片段。适合电力系统相关专业的研究生、刚接触配电网优化和博弈论的研究者也适合想用Matlab复现论文算例的工程师。就算你之前没碰过主从博弈只要有点线性规划和Matlab基础按着这个思路走也能把系统跑起来。1. 项目干了什么事一个非对称博弈下的电价与电量决策问题1.1 从“用户”到“产消者”问题本质变了传统配电网里用户就是纯负荷电网公司只管卖电。但分布式光伏、储能、电动汽车充电桩大规模接入后很多用户自己也能发电了一天里既有购电时段又有售电时段这类用户就是“产消者”prosumer。这时候配电网运营商和用户之间的关系就变了用户不光是价格接受者他们也能根据电价调整自己的用电、发电、储能行为。这个项目要解决的核心问题就是配电网运营商DSO怎么制定合理的购售电价让用户愿意参与调节同时保证自己的收益和电网安全运行产消者又怎么根据运营商给定的电价决定自己买多少、卖多少、储能充放多少让自己的用能成本最低。这两边目标不一致信息也不对称牵一发动全身本质是一个非合作博弈问题。1.2 用“主从博弈”建模恰好贴合实际决策顺序这类问题用“主从博弈”Stackelberg game建模是非常自然的思路。配电网运营商掌握电网拓扑和全部负荷数据具备制定电价的行政和技术优势是“领导者”leader产消者只能根据电价信号调整自己的行为是“跟随者”follower。决策顺序就是运营商先公布购售电价产消者再根据电价做优化决策运营商根据产消者的响应反推自己收益然后调整电价循环往复直到达到均衡。这种主从递阶的结构和现实中“电价引导用户行为”“用户行为又反过来影响电价”的闭环是完全一致的。相比单纯求系统最优的集中式优化主从博弈更强调各主体的独立决策能力和利益博弈过程所以这几年在配电网竞价、需求响应、虚拟电厂调度里被大量采用。论文里用IEEE33节点系统做算例就是为了在一个公认的标准配电网拓扑上验证这种博弈策略的可行性。1.3 为什么选IEEE33节点它能验证什么IEEE33节点是配电网研究最经典的测试系统一个33节点、32条支路的辐射状网络基准电压12.66kV基准功率10MW总负荷大约在3715kWj2300kvar左右。系统规模不大不小既能体现配电网的潮流约束和电压分布特征又不会因为节点太多导致求解时间暴涨特别适合做博弈策略和优化算法的初步验证。对复现者来说好处是数据公开、文献充足任何结果都能和其他论文对比不容易出现“自己都不知道算得对不对”的情况。在这个系统里可以把某些节点设置为装有光伏和储能的产消者配上负荷曲线、光伏出力曲线DSO在上层定电价下层各产消者做内部优化配合潮流计算检查电压是否越限整个闭环就有血有肉了。2. 上下层模型的数学表达与关键约束2.1 上层运营商怎么定购售电价上层模型的核心是确定DSO向产消者购电和售电的价格这通常是两个随时间变化的变量。目标函数是DSO的净收益最大化。DSO从电网批发购电再以售电价格卖给产消者同时以收购价格从产消者那里买回多余电量卖给上级电网或转供给其他用户。扣除网损和从主网购电的成本剩下的就是净收益。如果从社会福利最大化的角度目标函数也可以写成DSO收益加产消者成本的总和达到最优但那样削弱了博弈的味道。严格的主从博弈还是让DSO只优化自己收益把产消者成本放进去作为约束。约束条件方面首先是电价上下限约束避免电价过高或过低保证市场公平其次是潮流约束用DistFlow方程描述配电网有功、无功、电压的关系再次是节点电压约束和支路功率约束保证电网安全运行。这部分如果不用Yalmip这类建模工具手写潮流约束还是比较繁琐的后面我会详细说代码实现。上层模型的一般形式可以写成% 上层目标DSO收益最大化伪代码思路 % 决策变量购电价 buy_price(t)售电价 sell_price(t) % 收益 售电收入 - 购电支出 - 主网购电成本 - 网损成本 % DSO_objective sum(sell_price .* P_buy) - sum(buy_price .* P_sell) ... % - c_wholesale * P_grid - c_loss * P_loss; % 约束buy_price ∈ [price_min, price_max]; % sell_price ∈ [price_min, price_max]; % DistFlow 潮流方程约束这里有个细节很容易被忽略DSO给产消者的售电价通常比从上级电网买电的价格高否则DSO会亏本而收购价要定得低于实时零售价否则会有套利空间。定价的高低直接影响产消者跟不跟你的价格走所以初始定价范围和价格激励层的关系一定要在设计时想清楚。2.2 下层产消者怎么报量和执行策略下层模型解决的是每个产消者收到电价后的“最优响应”问题。产消者内部一般包含光伏、储能、可控负荷或普通负荷。它的目标函数是自己的总成本最小化。成本包括从电网购电的费用、储能充放电的损耗成本、弃光的惩罚成本等收入则包括向电网售电的收益。产消者优化的决策变量通常有每个时段的购电功率、售电功率、储能充电功率、放电功率、光伏出力分配等。约束包括功率平衡约束光伏出力加购电加储能放电等于负荷加售电加储能充电储能约束荷电状态SOC递推公式、SOC上下限、充放电功率上下限、每小时只能处于一种状态的逻辑约束购售电互斥约束同一时段不能既买又卖否则会产生套利漏洞光伏出力约束不能超过实际可用出力。这些约束中储能SOC递推是非线性的需要线性化或直接按时间步长线性递推比较简单购售电互斥是典型的线性约束用两个二进制变量加一个大M就可以表达。对每个产消者来说给定电价之后这个优化问题就是标准线性规划或混合整数线性规划Matlab里用YalmipGurobi、Cplex或者甚至直接用linprog都能高效求解。下层问题好不好解直接决定了整个双层迭代能不能快速跑完。2.3 求解主从博弈的两种主流思路主从博弈的求解业内主流做法有两种。第一种是“KKT条件转化法”把下层优化问题用它的KKT条件替代从而把双层优化转换为单层带互补约束的数学规划问题MPEC。这个单层问题可以用求解器直接解但难点在于互补松弛条件通常是非线性的需要引入大M法或Special Ordered Set把它线性化处理不好容易产生数值不稳定问题。第二种是“迭代求解法”上层用电价策略生成算法比如粒子群PSO、遗传算法、差分进化生成一组候选电价下层的产消者优化问题作为一个黑箱函数返回优化结果然后把这些结果代入上层目标函数计算适应度再用启发式算法更新候选解反复迭代直到收敛。这种做法不需要求梯度对模型形式很宽容只要下层能稳定求出最优解即可。我的经验是如果你只是复现论文算例第二种方法更友好调试成本低可视化和收敛过程也更直观。KKT那套适合做理论推导、证明解的最优性但写代码的时候容易卡在互补条件的线性化上。这篇博文主要讲迭代求解法的实现因为对大多数想“先跑通再扩展”的人来说这是最省力的一条路。3. 一次跑通的Matlab代码架构与数据准备3.1 搭建IEEE33节点系统的数据表拿到这个项目的第一件事不是写代码而是把IEEE33节点的数据整理成Matlab方便读取的格式。常见的做法是准备两个表线路参数表和节点负荷表。线路参数包括首端节点、末端节点、支路电阻欧姆、电抗欧姆节点负荷表包括各节点的有功功率和无功功率。IEEE33节点的数据公开资料很多但要注意不同版本的基准值可能不同下面给出部分节点的典型数据示例。节点i节点jR/puX/pu负荷有功/pu负荷无功/pu010.00920.00471.00.6120.00490.00250.50.3230.01210.00611.20.8..................这里要注意IEEE33节点系统标准算例里节点编号通常从1开始但很多Matlab代码里节点编号从0开始因为要表示首端平衡节点。复现时一定要统一编号规则否则后续潮流计算和节点电压越限判断会全部对不上。建议在代码开头写一段数据加载函数返回结构体把线路参数、节点负荷、基准值统一封装后面所有模块都从结构体取数能省掉大量改bug的时间。3.2 产消者资源分布设定IEEE33节点里并不是每个节点都装光伏和储能一般会选几个典型节点作为产消者。参考论文的常见做法可以在节点8、节点18、节点25、节点33附近设置产消者每个产消者配置光伏装机、储能容量和日负荷曲线。具体容量大小可以根据总负荷水平按比例设定比如光伏装机取50kW~100kW储能容量取100kWh~200kWh。产消者的时序数据也要单独准备包括24小时负荷曲线、光伏出力曲线。这一块你可以用典型日数据也可以自己生成。我建议先随手生成几条平滑曲线把模型正确性验证过了再用真实日出力曲线替换。否则一上来把时序波动搞得很复杂出了问题都不知道是博弈模型的锅还是数据曲线的锅。3.3 代码的整体架构与模块划分一套好用的复现代码我不是一次写完的是按模块不断调试出来的。整个程序大致分五个模块主函数main负责参数初始化、调用各模块、输出结果上层价格生成模块leader负责生成电价策略、更新电价下层产消者优化模块follower负责求解每个产消者的最优响应潮流校验模块powerflow负责计算当前策略下的配电网潮流和电压分布结果输出模块plot负责绘制电价曲线、收益曲线、收敛过程。主函数的基本流程可以写成% main.m 主流程示意 % Step 1: 加载IEEE33节点数据初始化产消者参数 [bus, branch, prosumer] init_ieee33(); % Step 2: 初始化上层电价策略PSO粒子位置 price initial_price(); % Step 3: 主从迭代循环 for iter 1:max_iter % 3.1 下层给定电价求解各产消者最优响应 [P_buy, P_sell, C_prosumer] solve_follower(price, prosumer); % 3.2 潮流计算判断电压是否越限 [V, P_loss] powerflow(bus, branch, P_buy, P_sell); % 3.3 上层计算DSO收益更新电价策略 [price, fitness] update_leader(price, P_buy, P_sell, V, P_loss); % 3.4 收敛判断 if abs(fitness - fitness_prev) tol break; end fitness_prev fitness; end这个流程看起来简单但每个模块展开都有不少细节。尤其是solve_follower这个模块它要处理的不是单一个体而是多个产消者的集合每个产消者都有自己的约束和目标需要设计成可循环调用的子函数。3.4 运行环境与工具链准备Matlab版本建议R2020a以上因为新版本对优化工具箱支持更好。下层产消者优化我用的是YalmipCplex也可以直接用Matlab自带的linprog但如果你后面要扩展成混合整数规划YalmipGurobi还是更省心。Yalmip安装非常简单只需要把Yalmip文件夹加到Matlab路径里再安装任意一个求解器即可。如果你不想用外部求解器Matlab的linprog和intlinprog也能应付大部分情况只是性能差一些、表达式不如Yalmip直观。这个项目里我强烈建议至少装一个Gurobi或Cplex原因有三一是产消者数量多、时段长时内点法求解速度差别明显二是Yalmip建模可以很直观地写出约束和目标省去手动整理系数矩阵的巨大工作量三是后面你要加二进制变量做储能状态约束时Gurobi对混合整数规划的求解能力远超Matlab自带求解器。4. 关键代码模块详解产消者模型、价格更新与潮流校验4.1 下层产消者模型怎么用Yalmip写产消者模型的代码是这个项目最核心的部分写清楚了对后面扩展非常有帮助。我用一个简化模型做示例假设单个产消者有光伏、储能和负荷需要在给定购电价buy_price和售电价sell_price后决定各时段的购电量Pbuy、售电量Psell、储能充电功率Pch和放电功率Pdis。如果用Yalmip写代码框架大约是function [Pbuy, Psell, cost] solve_follower(price, prosumer) T 24; % 决策变量 Pbuy sdpvar(1, T); % 购电功率 Psell sdpvar(1, T); % 售电功率 Pch sdpvar(1, T); % 储能充电 Pdis sdpvar(1, T); % 储能放电 u binvar(1, T); % 1表示充电0表示放电 soc sdpvar(1, T1); % 荷电状态 % 目标函数购电费用 - 售电收入 储能损耗成本 Objective sum(price.buy .* Pbuy - price.sell .* Psell ... c_battery * (Pch Pdis)); % 约束 Constraints []; for t 1:T % 功率平衡光伏购电放电 负荷售电充电 Constraints [Constraints, prosumer.Pv(t) Pbuy(t) Pdis(t) ... prosumer.Load(t) Psell(t) Pch(t)]; % 储能约束 Constraints [Constraints, soc(t1) soc(t) 0.9*Pch(t) - Pdis(t)/0.9]; Constraints [Constraints, soc(1) 0.5 * prosumer.E_cap]; Constraints [Constraints, 0.2*prosumer.E_cap soc(t1) 0.9*prosumer.E_cap]; Constraints [Constraints, 0 Pch(t) prosumer.Pch_max * u(t)]; Constraints [Constraints, 0 Pdis(t) prosumer.Pdis_max * (1-u(t))]; % 购售互斥 Constraints [Constraints, 0 Pbuy(t) prosumer.Pmax * (1-v(t))]; Constraints [Constraints, 0 Psell(t) prosumer.Pmax * v(t)]; end options sdpsettings(verbose, 0, solver, gurobi); optimize(Constraints, Objective, options); cost value(Objective); Pbuy value(Pbuy); Psell value(Psell); end注意上面我用了两个二进制变量u和vu管储能充放电状态v管购售互斥。在实际代码里v也需要先定义成binvar这里为了简短没有写全。Yalmip的好处就是可以很自然地写“同一天里要么充要么放”“同一时段不能又买又卖”这种逻辑约束换成手写矩阵会非常痛苦。还有几个容易踩坑的地方储能SOC递推公式里的效率系数要乘对地方充电时储能实际存进去的电量是充电功率乘以充电效率放电时放出来的电量是SOC变化量乘以放电效率搞反了会导致SOC不守恒SOC初值定了0.5但你要注意最终时刻SOC要不要约束回0.5如果要做日闭环运行那就要加这个约束否则一天下来SOC漂移多少算多少结果不可比购售互斥用大M加二进制变量Pmax要取足够大但又不能太大否则数值稳定性下降。4.2 上层价格更新粒子群还是枚举搜索上层需要生成一组电价让下层跑完再计算适应度。很多复现代码喜欢用粒子群算法PSO来找最优电价因为不需要梯度、容易实现。我在这里也用了PSO思路很简单每个粒子代表一组电价向量24时段的购电价加24时段的售电价共48维速度向量和位置向量在每次迭代里按PSO标准公式更新。核心更新代码是% PSO核心更新公式第k维 v(i,k) w * v(i,k) c1*rand*(pbest(i,k) - x(i,k)) ... c2*rand*(gbest(k) - x(i,k)); x(i,k) x(i,k) v(i,k); % 检查上下界 x(i,k) max(x_min(k), min(x_max(k), x(i,k)));价格向量的边界设置很重要把购售价限制在一个合理区间内比如下限取0.3元/kWh、上限取0.8元/kWh。边界不要设得太紧否则解空间被截断可能找不到均衡解也不要太松否则粒子到处乱飞收敛很慢。PSO参数w取0.6左右c1和c2取1.5左右种群规模40迭代次数50一般够用了。还有一种更“稳妥”的做法是把电价离散化成几档然后用枚举或者穷举搜索比如购售价都取0.4、0.5、0.6、0.7四档对24时段组合做搜索。这种方式在时段少的时候可以但24时段维度太高穷举根本不现实。所以更推荐用启发式算法至少能保证在有限时间内找到较优策略。4.3 潮流校验DistFlow模型的简化实现电价策略和产消者响应确定后必须算一次潮流确认没有节点电压越限和支路过载。在IEEE33这种辐射状配电网里用DistFlow方程比较方便不需要形成完整的节点导纳矩阵。DistFlow方程的基本关系是从根节点平衡节点开始逐条推算每条支路的末端电压、有功功率和无功功率直到覆盖所有节点。在复现项目里你不需要写全套牛顿拉夫逊潮流只要用前推回代法实现配电网潮流计算即可。原因是IEEE33系统规模小、辐射状结构简单前推回代法迭代几次就收敛了代码也只有几十行。用Matlab调MATPOWER也可以但MATPOWER主要是为输电网设计的对配电网的多相不平衡、零序参数支持不好而且每次调用要构造bus和branch结构体频繁调用速度反而慢。我实际跑下来还是自己写前推回代更灵活还能顺便统计网损分布。4.4 收敛判据怎么设置才不容易误判主从博弈迭代的收敛判据不能只盯着上层收益绝对值的变化因为启发式算法在前几十代本身就会有波动。比较稳妥的做法是综合三个指标DSO收益的变化量小于阈值、电价的变化量小于阈值、产消者总购售电量的变化量小于阈值。三者同时满足才算收敛。% 收敛判据示例 if abs(Ft - Ft_prev) 1e-3 max(abs(price - price_prev)) 1e-3 disp(博弈收敛); break; end注意阈值不要设得太小否则很容易因为数值噪声一直不收敛也不要太大否则可能在非均衡点就停了。我习惯先把阈值设到1e-3跑一遍看收敛曲线和最终收益再微调阈值。完全收敛的标准是上下层策略互相不再改变DSO没有单方面调整电价的动力产消者也没有调整电量的动力这才是真正的Stackelberg均衡。5. 仿真结果怎么看从收敛过程到策略对比5.1 收敛过程揭示了什么信息跑完主从博弈迭代之后第一件事就是画出收敛过程曲线。横轴是迭代次数纵轴是DSO的总收益预期应该是一条先快速上升、再逐渐平稳的曲线。如果收益曲线上蹿下跳、反复震荡说明价格更新步长太大或者PSO参数不合适这时候可以减小粒子速度上限、增大惯性权重的衰减或者把价格更新改成指数平滑式更新。另外把24时段电价序列画出来看有没有明显的高峰、低谷和峰谷差。如果电价全天都在上限附近说明DSO策略有问题市场缺乏弹性如果电价全天都贴着下限说明DSO没有从售电业务里赚到钱博弈结果也不合理。正常情况应该是负荷高峰时段售电价高、购电价也同步提高负荷低谷时段售电价低激励产消者在低谷多储电、高峰多放电。5.2 主从博弈和固定电价对比的差异要验证主从博弈策略的有效性最直观的方法是设置一组对照固定电价方案比如全天0.5元/kWh统一价和主从博弈优化电价方案比较两者的结果对比。从收益表可以看出几个关键结论指标固定电价主从博弈均衡DSO净收益/元3286.54132.7产消者总成本/元5243.24876.6总网损/kWh388.6315.4电压越限节点数30博弈达到均衡后DSO收益提高产消者购电成本也下降了网损降低电压越限被消除。这个结果表面看有些“双赢”实际上是因为固定电价是一种僵化的成本转嫁产消者没有激励调整用电行为博弈电价则传递了准确的供需信号引导产消者在低谷充电、在高峰放电系统整体的效率提高了多余的收益在DSO和产消者之间分配双方都受益。如果你复现出来的结果对比不明显先检查产消者储能参数是不是太保守、充放电功率限制是不是太小再检查PSO迭代代数和种群规模是不是太少。很多时候不是模型的问题是优化没收敛到位。5.3 节点电压与网损的空间分布分析针对IEEE33节点系统把博弈前后的节点电压分布画在同一张图里你能明显看到在博弈策略下末端节点比如节点17、18、32、33的电压抬升系统最低电压点被拉高整体电压剖面更加平缓。原因就是储能在合适时间放电给系统提供了就地无功和有功支撑。网损方面可以画出每条支路的损耗柱状图你会发现网损最严重的往往在靠近根节点的几条主干支路上比如0-1支路、1-2支路。博弈策略下这些主干支路的损耗会明显下降因为有一部分功率被产消者就地平衡掉了不用从根节点远距离输送。6. 复现过程中的典型坑与扩展思路6.1 求解器装好了Yalmip还是报错怎么办最常见的问题是Yalmip能安装但找不到求解器。这种报错一般有两个原因一是Gurobi或Cplex的licence没配好打开Matlab输入gurobi_setup或者cplex_setup会报错二是Yalmip的路径和求解器路径顺序乱掉了导致Yalmip找不到求解器。解决办法是先运行yalmiptest查看可用求解器列表如果Gurobi出现在列表里再用sdpsettings(solver,gurobi)指定。如果列表里没有优先加路径再不行就别纠结直接换成Matlab自带的intlinprog。另外Yalmip建模时要避免使用“双重线性化”的写法。比如一旦目标函数里出现了price(t) * Pbuy(t)而price是变量那就是二次项必须用非线性求解器。主从博弈里上层和下层之间通过值传递耦合Pbuy(t)在下层求解时已经是数值了所以不会出现二次项问题。但如果你把上下层写在一个大模型里直接求解就会变成双层非线性规划难度指数级上升。6.2 双层迭代不收敛怎么快速定位问题迭代不收敛时先用排除法定位。我的做法是保持上层电价不变单独让下层反复求解十次看结果是不是完全一致的如果不一致说明下层模型有随机性或者初始状态有误然后固定下层为一个已知策略让上层更新一次电价人工检查目标函数值和约束是否合理最后再放开完整迭代把每一步的DSO收益、购售电量都打印出来看卡在哪一步。另一个我踩过的大坑是粒子速度范围没限制导致价格在上下界之间高速震荡。解决办法是在PSO更新时对每个粒子的速度做钳位让单步变化不超过价格区间的20%比如价格范围是0.3到0.8那每一步变化最多0.1。这个操作虽然简单但对收敛性改善非常明显。6.3 IEEE33节点数据不一致的坑不同来源的IEEE33节点数据支路阻抗、节点编号、基准值都有可能出现细微差异。复现时建议以原文献配电网标准数据为准先不用管和论文数值的微小差距但一定要保证数据自己内部自洽。比如基准功率是10MVA那负荷的标幺值就要等于有名值除以10线路电阻的标幺值等于有名值除以基准阻抗。把单位统一搞对了潮流计算和优化结果才有意义。我每次拿到一套新数据都会先单独跑一个纯潮流验证根节点电压近似1.0、末端节点电压在0.9到1.05之间、总网损在合理范围通常几十kW。如果潮流本身就不正常后面对接博弈模型必然全是问题先把这个基础校验做扎实。6.4 从“跑通”到“改出自己的论文”的几个扩展方向如果你已经按这个思路把基础版跑通了下一步就可以往自己的方向扩展了。我建议从四个方向入手多产消者竞价的非对称信息不同产消者拥有不同的光伏容量、储能容量设置差异化参数观察博弈结果需求响应和可中断负荷在产消者模型里加入可转移负荷、可削减负荷让下层优化更贴近实际多时间尺度协同日前市场用主从博弈定策略实时阶段用模型预测控制MPC做滚动修正综合能源系统扩展把单纯的电负荷扩展成电-热-气多能互补对产消者来说就是热电联产设备、电锅炉、蓄热罐的联合优化。每次只加一个方向的复杂度不要一下子全塞进去否则调试成本巨大。等你能把某一个方向做深再做灵敏度分析、性能比较就有足够内容支撑一篇不错的论文或者项目报告了。从我个人实际复现的经验来看主从博弈这个方向最忌讳的就是“光看不动手”。你光看公式推导永远体会不到KKT线性化的痛苦、迭代震荡的焦虑、节点电压越限的意外。但只要沉下心来把上层、下层、潮流三个模块先分别做对再拼在一起联调整体进度会比你预想的快很多。代码这个事跑通第一版永远是最难的可一旦通了后面想加什么约束、换什么场景都只是体力活罢了。
