复现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.95,SOC初始值0.5,SOC范围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(T+1,1); C = []; C = [C, SOC(1) == param.SOC0 * param.E]; for t = 1:T C = [C, SOC(t+1) == 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,键是价格向量的哈希,值是对应的真实响应。这样一来,同一组价格在迭代中再次出现时,直接查表返回结果,整个复现时间能省掉一半以上。这种工程细节,论文里永远不会写,但对实际跑通项目的帮助是实打实的。