复现EI论文的活干多了你会发现最难的不是把摘要里的公式看懂而是把论文里没写的工程细节一点点补出来。这两天刚把一个基于元模型优化算法的多虚拟电厂主从博弈动态定价与能量管理算例在Matlab里完整跑通代码是论文附带的底子里面很多环节需要我们做研究的人自己拿主意。这篇笔记把整个项目的建模思路、求解框架、代码结构、调试过程都盘一遍给准备碰虚拟电厂、双层规划、博弈优化方向的朋友做一个参照。这个项目很适合在读研究生复现参考也适合工程师快速理解动态定价机制在VPP集群协调里到底怎么落地。先说清楚这里面的几个关键词免得后面对不上号元模型优化算法国内论文里也叫代理模型、替代模型英文对应metamodel或者surrogate model本质是用一个便宜的函数去替代昂贵的真实评估过程主从博弈就是Stackelberg博弈一个领导者先出价多个跟随者再响应多虚拟电厂就是多个VPP各自聚合分布式电源、储能和可调负荷动态定价和能量管理则是上层定价、下层调度的完整闭环。这四个东西串在一起就是一个典型的多主体协调优化问题。1. 项目在解决什么问题多虚拟电厂定价与能量管理的整体设计思路1.1 多VPP协同场景里定价和能量管理为什么绑在一起虚拟电厂不是一台具体设备它是一个聚合体。光伏、风电、储能、可调负荷这些分散的资源被聚合成一个整体以类似电厂的形态参与市场或接受调度。问题在于多个VPP往往归属不同运营主体各自的资源结构、成本函数、风险偏好都不一样谁也没有权力直接命令别的VPP怎么调度集中式优化在这类场景里是走不通的。价格信号是最自然的协调手段。上层运营方定一个动态电价各个VPP看到这个价格后会基于自己的内部约束做能量管理决定多买电、少买电还是干脆把储能放出来卖电。这个逻辑很像我们在电商平台看到的定价机制平台改价格消费者改购买量平台再根据购买量调整价格最后收敛到一个双方都能接受的状态。只不过在电力系统里这个过程有物理约束比电商补货复杂得多。所以动态定价和能量管理天然就是一对主从关系。上层管价格下层管电量响应两者必须放在同一个模型里算单独看哪一层都不完整。做这个复现项目的时候我也是先把这层业务逻辑理顺了才开始写代码否则连上层的目标函数定义成什么都很容易搞错。1.2 主从博弈框架的建模选择主从博弈在博弈论里叫Stackelberg博弈。它描述的是一个先后决策的场景领导者先出一个策略跟随者看到之后做出最优反应领导者预测到跟随者的反应之后再选择让自己利益最大化的策略。在这个项目里上层是虚拟电厂聚合商或者配电系统运营商负责制定各时段的购售电价格下层是多个VPP根据价格做内部能量管理。用Stackelberg而不是集中式优化不是数学上的偏好问题而是现实约束决定的。集中式优化要求所有VPP把内部信息完整上报储能SOC、柔性负荷的舒适度约束、分布式电源的检修计划这些信息涉及运营隐私和商业利益VPP不一定愿意交出来。主从博弈只需要上层给价格下层返回购电量和售电量信息交互最小更符合实际市场机制。我在复现时选的是典型的三层结构上层定价中层市场交互下层VPP内部优化。也可以理解为二层结构因为中层只是买卖电量的结算环节。不同论文的层级划分略有区别但核心思想是一致的上层做价格决策下层做能量管理决策通过迭代或者等价变换求均衡解。1.3 元模型优化算法在这里的真实作用如果不用元模型这类问题怎么解最常见的是迭代法上层先给一组初始电价调用下层优化求解所有VPP的响应然后用梯度法或启发式算法更新电价再重新调用下层。问题是下层每个VPP都是带储能的优化模型这一步本身就需要求解一个混合整数规划三个VPP就要解三次。上层要搜索整个价格空间往往需要评估上百组甚至几百组候选价格算下来就是几百次下层优化求解计算量瞬间就上来了。元模型解决的就是这个痛点。它的思路非常朴素:先用一批采样价格求解真实的下层优化拿到一组“价格-响应”或“价格-上层收益”的样本然后训练一个便宜的替代模型比如高斯过程回归或者RBF神经网络。之后上层的大部分搜索评估都在这个代理模型上完成只有少数有潜力的候选点才拿去做真实的下层求解再把这些真实样本加进训练集反复迭代让代理模型越来越准。这个思路在我复现的工程里效果很明显。没有代理模型之前跑一次动态定价加能量管理的完整算例可能要二十分钟加了元模型之后几分钟内就能收敛到比较稳定的解。这也是这类EI论文为什么喜欢把元模型和博弈框架放在一起的原因思路新颖又能实实在在解决嵌套优化的计算瓶颈。2. 数学模型与关键公式从博弈论到可计算的形式2.1 上层动态定价模型上层模型的核心决策变量是各时段的动态电价。我这里按24个时段建模每个时段可以单独定价但实际运行中如果价格维度太高采样和代理模型拟合都会遇到困难后面我会讲怎么降维。上层优化目标常见有两种定义方式一种是最大化聚合商收益即从VPP售电收入减去从上级电网购电的成本另一种是最大化系统社会福利把VPP的成本也纳入目标函数。我复现阶段采用第一种因为它在商业逻辑上更直观也更容易解释均衡结果。如果不加约束上层可能会把电价推到上限因为VPP无论如何都得用电价格越高聚合商收入越高。所以动态定价模型必须带上限约束通常还会设置一个合理的峰谷价差区间防止价格信号过度扭曲。约束可以写成各时段电价在0.3到1.2元/kWh之间同时峰时段与谷时段价差不超过某个阈值。这些参数在不同论文里取值不同但作用都是把价格控制在一个可接受的范围内。2.2 下层VPP能量管理模型下层每个VPP在收到电价后解决的是一个典型的日前调度问题。决策变量包括从上级电网的净购电量、储能充放电功率、光伏出力的削减量以及柔性负荷的调整量。目标是最小化运行成本包括购电成本、储能损耗成本和可能的弃光惩罚。约束条件是这个模型的重头戏。第一是功率平衡约束发电加购电加储能放电要等于负荷加储能充电。第二是储能约束包括SOC递推方程、充放电功率上下限、充放电互斥约束以及SOC的上下限。第三是联络线约束VPP与上级电网交互的功率不能超过线路允许值。如果考虑负荷可调还需要加上可调范围和能量守恒约束。储能约束里最容易出错的是SOC递推的单位。我复现时搞混过一次SOC的单位是MWh但算例里储能容量给的也是MWh如果充放电功率用的是MW乘上1小时就能对得上如果调度步长不是1小时一定要乘上时间步长否则SOC曲线会凭空多出一块能量。2.3 均衡求解路径与元模型加速求解主从博弈均衡工程上主要有三条路线。第一条是把下层模型用KKT条件替换把双层问题转成单层数学规划再用现成的商业求解器求解。这条路线数学严谨但要求下层模型是凸的而且多个VPP的KKT条件堆在一起变量和约束会膨胀得很厉害YALMIP建模时特别容易出错。第二条是分布式迭代。给定电价用优化器解每个VPP的子问题返回响应后更新电价循环直到收敛。这条路线实现最简单但计算量大而且电价更新步长设不好会震荡。第三条就是我复现用的元模型加速方案。大体流程是先用拉丁超立方采样在价格可行域内生成一批初始候选点求解下层真实优化拿到响应数据和上层目标值训练高斯过程回归模型然后在代理模型上寻找最有潜力的新候选点对候选点做真实结算加进样本库重新训练。这里有两个细节值得注意第一我选择直接对上层目标函数建代理因为它是一个标量拟合稳定如果你选择对VPP响应函数建代理则需要在多维输出上做处理难度会高不少。第二加点准则我用的是期望改进EI它能在探索和利用之间做平衡有效避免代理模型只盯着已经找到的最优点。3. Matlab代码实现与关键模块解读3.1 工程文件结构与整体流程复现这种带博弈和代理优化的项目代码组织特别重要。我实际跑的时候把整个工程拆成了这样几个文件main.m主程序负责初始化、循环、结果输出data_loader.m读取负荷曲线、光伏曲线等基础数据sim_setup.m设置VPP数量、储能参数、价格边界等upper_objective.m计算上层聚合商的目标函数值lower_vpp.m求解单个VPP的能量管理子问题surrogate_train.m训练和更新高斯过程代理模型acquisition.m计算加点准则选择新的候选价格plot_results.m画动态定价曲线、VPP响应曲线和SOC曲线文件之间的调用关系很清晰main调用sim_setup初始化参数调用data_loader读取曲线然后进入外层迭代循环。在循环里surrogate_train负责拟合代理acquisition负责选点lower_vpp负责真实结算upper_objective负责评估目标。把不同功能拆开的好处是后面哪一步出了问题直接单独调试那个函数就行不用在几百行的主程序里翻来翻去。3.2 参数设置与数据准备数据准备这一步很多人不重视其实最影响复现结果。我这里把默认参数列出来你可以直接抄作业三个VPP调度周期24小时电价下限0.3元/kWh上限1.2元/kWh储能容量分别取10、15、12MWh储能充放功率上限分别取1、2、1.5MW充放电效率0.95SOC初始值0.5SOC范围0.2到0.9光伏装机分别取20、30、25MW负荷峰值分别取30、40、35MW。负荷曲线和光伏出力曲线我放在CSV文件里用readtable读入。这里推荐直接用readtable不要再用老旧的csvread或者importdata前者对表头、数据类型的处理要省心得多data readtable(vpp_profile.csv); load_curve [data.Load_VPP1, data.Load_VPP2, data.Load_VPP3]; pv_curve [data.PV_VPP1, data.PV_VPP2, data.PV_VPP3];读完之后第一件事不是写优化模型而是画一下曲线看看数据有没有明显异常。我接过不少复现项目很多人卡在最后结果对不上排查半天发现是光伏曲线单位和负荷曲线差了一千倍。这种低级错误画图一眼就能看出来。3.3 上层定价循环与代理模型拟合上层优化的第一步是在价格可行域内生成初始样本。我推荐用lhsdesign做拉丁超立方采样它的好处是能保证样本在价格空间里分布得比较均匀不会像纯随机采样那样扎堆。样本数量上如果价格维度较低一般取40到60个初始点就够了如果价格维度高到接近24维那建议你先把价格归并成峰平谷几段再采样否则样本量不够代理模型很难拟合准。代理模型训练用fitrgp这是MATLAB统计和机器学习工具箱里的高斯过程回归函数。核心的调用如下gprMdl fitrgp(price_samples, response, ... KernelFunction, ardsquaredexponential, ... Standardize, true);核函数选ardsquaredexponential也就是ARD形式的平方指数核。ARD的意思是每个输入维度自动学一个长度尺度参数对价格维度重要性差异比较大的场景特别有用比如某个时段的电价对VPP响应影响小自动会被分配一个较长的长度尺度。Standardize设为true把训练数据标准化能提升拟合的数值稳定性。在代理模型上选新的候选点我用的是EI加点。EI的核心思想是某个候选点的价值不仅要看它预测的目标值有多好还要看它的不确定性有多大两者综合起来算一个期望改进。下面的代码演示了EI的一种简化写法[mu, sigma] predict(gprMdl, candidate_prices); improve (mu - best_so_far); term1 improve .* normcdf(improve ./ max(sigma, 1e-10)); term2 sigma .* normpdf(improve ./ max(sigma, 1e-10)); EI term1 term2; [~, idx] max(EI); new_price candidate_prices(idx, :);拿到new_price后调用真实的下层优化去求解得到真实的响应和上层目标值然后把它加入样本集重新训练。整个循环就完成了“采样-拟合-选点-加点-再拟合”的闭环。3.4 下层VPP决策模块下层VPP模型是整个项目里最容易被写崩的地方尤其是储能充放电互斥约束。我这里用YALMIP建模逻辑会更直观。YALMIP不是MATLAB自带工具箱需要自己去GitHub下载安装但它能大幅降低建模复杂度强烈建议装一个。单个VPP的核心模型可以这样写function [cost, q] lower_vpp(price, pv, load, param) T param.T; p_ch sdpvar(T,1); % 储能充电功率 p_dch sdpvar(T,1); % 储能放电功率 x_ch binvar(T,1); % 充电状态标志 x_dch binvar(T,1); % 放电状态标志 q sdpvar(T,1); % 净购电量负值表示售电 SOC sdpvar(T1,1); C []; C [C, SOC(1) param.SOC0 * param.E]; for t 1:T C [C, SOC(t1) SOC(t) ... (param.eta_ch * p_ch(t) - p_dch(t) / param.eta_dch) * param.dt]; C [C, p_ch(t) 0, p_ch(t) param.Pch_max * x_ch(t)]; C [C, p_dch(t) 0, p_dch(t) param.Pdch_max * x_dch(t)]; C [C, x_ch(t) x_dch(t) 1]; C [C, SOC(t) param.SOC_min * param.E]; C [C, SOC(t) param.SOC_max * param.E]; C [C, q(t) pv(t) p_dch(t) load(t) p_ch(t)]; end objective sum(price .* q) param.c_om * sum(p_ch p_dch); ops sdpsettings(solver, gurobi, verbose, 0); optimize(C, objective, ops); cost value(objective); q value(q); end这段代码里x_ch和x_dch是二进制变量配合大M约束限制储能不能同时充电和放电。如果你不想装YALMIP也可以用MATLAB自带的intlinprog把二进制变量显式拼进决策变量向量只不过约束矩阵要自己手工拼调错一个索引就会出非常隐蔽的问题整体开发效率低很多。功率平衡约束里我把q定义为净购电量正数是从电网买电负数表示向电网售电。这样一来q可以直接和上层收益对接上层收入等于电价乘所有VPP的净购电量和。注意如果q为负上层收入会减少这正好反映了聚合商在低谷时段从VPP买电再转售给其他用户或上级电网的商业逻辑。3.5 主程序循环与收敛判据主程序的循环逻辑其实不复杂我把它概括成下面的步骤先初始化参数和样本接着做若干轮迭代每轮训练代理、选点、真实结算、判断收敛。收敛判据我用了两个满足任意一个就退出一是价格变化量小于阈值二是EI改进值小于阈值。价格更新这里有个大坑。如果直接拿EI选出来的最优价格替换当前价格经常会出现震荡也就是这轮选高下轮选低永远稳定不下来。我的做法是加一个阻尼系数新价格等于上一次价格加一个较小的比例乘以候选价格和上一次价格的差。阻尼系数从0.3开始调如果还震荡就降到0.1。这个技巧虽然简单但能省掉大量调试时间。外层迭代次数一般设50到100次就够用了没必要设太大。每轮迭代最耗时的就是一次真实下层求解如果用Gurobi一次三个VPP的求解大概在几秒量级总的跑下来也就几分钟。如果用的是MATLAB内置求解器时间会稍微长一点但也属于可接受范围。4. 复现结果、调试实录与避开常见坑4.1 如何判断复现结果是对的这个问题看起来基础却是所有复现项目里最考验经验的环节。结果不是跑出来就能交差的必须能解释得通。我一般会先看三条曲线动态定价曲线、各VPP的购电响应曲线、储能SOC曲线。动态定价曲线应该呈现明显的峰谷形态晚高峰时段电价走高深夜低谷时段电价走低。各VPP的购电响应应该和电价反向电价高的时候VPP购电量下降甚至变为负值也就是开始卖电。储能SOC曲线则应该表现为电价低时充电SOC爬升电价高时放电SOC下降。如果这些行为对不上比如价格高的时候储能反而在充电那基本可以断定模型里某个约束写反了或者符号定义错了。还有一个判断方法是看目标函数随迭代的变化。正常情况下上层收益应该是先快速上升然后慢慢平稳形成一条类似学习曲线的形状。如果目标函数值反复跳动说明收敛判据或者阻尼系数有问题需要回头检查。4.2 现场遇过的典型问题与处理建议复现过程中我踩了不少坑整理成表格方便你对照排查。现象可能原因解决方案迭代价格震荡不收敛价格更新步长太大加阻尼系数p_new p_old 0.3*(p_candidate - p_old)震荡则继续降到0.1代理模型预测严重偏离真实值初始样本太少或价格维度太高把24小时价格归并成峰平谷3到5段初始样本N0提高到80以上下层优化返回无可行解功率平衡约束配平不了允许q为负数即允许VPP向电网售电或放宽联络线功率上限SOC曲线出现突变跳变SOC递推单位没乘时间步长检查dt是否等于1如果不是1小时必须乘上dt结果曲线锯齿明显收敛阈值太松把上层价格变化阈值从1e-2降到1e-4同时增大内层样本量GPR训练报错或者预测值全是常数核函数参数不合适或者数据未标准化fitrgp里设Standardize为true替换核函数为ardsquaredexponential这些坑里最隐蔽的是第一个和第四个。价格震荡通常是阻尼系数的问题但如果你模型里同时有多层循环还要确认一下是不是内外层判据写反了。SOC跳变如果不是单位问题那就要看储能初始SOC和最终SOC的边界条件很多论文要求调度周期末SOC回到初始值附近这个约束不加的话储能会在最后一个时段把能量全部放光曲线看起来就像断崖一样。4.3 MATLAB版本、工具箱与求解器选型建议MATLAB版本对复现结果影响不大但工具箱必须装全。这个项目依赖Optimization Toolbox、Statistics and Machine Learning Toolbox如果自定义启发式算法还想用Global Optimization Toolbox。平时我跑这种项目用的是R2023b其实R2021b往上的版本差别不大核心还是工具箱是否完整。有一个验证工具箱的小技巧。在MATLAB命令行里敲ver命令它会列出所有已安装的工具箱清单。如果运行fitrgp时报“未定义函数或变量”十有八九是统计和机器学习工具箱没装直接看ver输出就能确认不用去代码里瞎猜。求解器方面YALMIP加Gurobi是当前学术复现场景里的黄金组合求解速度比内置的intlinprog快不少。如果只是先跑通逻辑用内置intlinprog也完全可以毕竟下层模型规模并不大三五个VPP、24小时变量数量级在几百个内置求解器足够应付。等确认模型没问题、需要跑大规模算例的时候再切换到Gurobi不迟。另外关于AI辅助写Matlab代码最近问的人很多。简单函数和数据处理AI写得又快又稳但像这种主从博弈加代理优化的闭环逻辑AI生成的代码容易在迭代顺序和样本更新上出错还是得自己把控整体结构。4.4 复现论文结果时的对数字技巧如果你拿到一篇EI论文作者没有公开完整代码只有算例参数和结果图复现时该怎么对数字我自己的经验是分三步走。第一步先对量级。论文里如果写目标函数是几百万那你的目标函数不应该是几十也不应该是几亿先把量级对齐才能继续往下比。量级不对八成是单位问题或者目标函数里漏了某一项成本。第二步对曲线的形状特征。峰谷出现的时间点是否一致储能开始放电的时段是否一致这些不需要精确匹配但趋势必须吻合。比如论文里的储能是在18点开始放电你复现出来在14点就开始放那大概率是负荷曲线或者电价边界设置有出入。第三步对具体数值。数值对不上不一定是代码错也可能是论文参数没写全。这时候可以结合曲线反推参数比如从SOC曲线斜率反推储能容量从购电曲线反推负荷基线。这一步比较费时间但一旦反推成功论文里的底层数据基本就摸清了后面的调参也就有了方向。做这种EI复现项目我的一个体会是不要迷信论文里写的每一个参数有些论文算例参数是从其他文献挪过来的有笔误很正常。你需要在合理范围内自己调整只要最终结果曲线在趋势和量级上和论文一致复现就算成功。最后再分享一个小技巧。代理模型加点的时候候选点经常重复踩在历史最优附近每次都要重新求解下层优化造成大量重复计算。我后来把所有已经解过的“价格-响应”配对存进一个containers.Map键是价格向量的哈希值是对应的真实响应。这样一来同一组价格在迭代中再次出现时直接查表返回结果整个复现时间能省掉一半以上。这种工程细节论文里永远不会写但对实际跑通项目的帮助是实打实的。
