虚拟电厂、主从博弈、动态定价这三个词经常同时出现在电力市场、综合能源、需求响应方向的论文摘要里但真正能跑通的MATLAB代码却很少公开。我最近把这套模型从数学推导一路做到仿真踩了不少坑也摸出了一些规律。这篇文章就把完整思路、模型公式、代码实现和调试经验一次讲清楚适合正在做课程设计、毕业论文或者准备投期刊仿真的同学参考也适合想了解虚拟电厂定价机制的工程师当作入门笔记来读。1. 虚拟电厂为什么需要主从博弈和动态定价1.1 虚拟电厂的本质是“聚沙成塔”虚拟电厂Virtual Power PlantVPP并不是一个实体的发电厂它本质上是一种管理机制把分散在用户侧的光伏、储能、电动汽车充电桩、温控负荷等资源聚合起来统一参与电网调度和电力市场交易。所谓“聚沙成塔”意思就是单个光伏只有几千瓦单台空调的调节能力也只有几百瓦但把几千户、几十栋楼的资源打包在一起就能形成几兆瓦到几十兆瓦的可调节能力这在电网调度里才真正有分量。做虚拟电厂能量管理模型的时候我一开始也踩过一个典型误区想把配电网拓扑、每台变压器潮流都建出来。后来发现这样做在纯MATLAB环境里计算量非常大而且对于主从博弈这种市场层面的问题配电网潮流细节反而会干扰定价机制的分析。更常用的做法是聚合功率模型只关注VPP整体与用户之间的功率交互以及VPP与上级电网之间的购售电交互。这样既能保留定价博弈的核心逻辑又能让代码在普通电脑上跑得动。在这个聚合模型里VPP运营商面对的决策包括内部零售电价怎么定、储能什么时候充放电、从外部电网买多少电、多余的光伏电量是否卖回电网。用户面对的决策则简单一些在当前电价下用多少电最划算。两边都想最大化自己的利益电价在里面起到了核心的“指挥棒”作用。1.2 主从博弈谁先出手谁跟在后面主从博弈这个词听起来很学术其实就是斯塔克尔伯格博弈Stackelberg Game核心特征是参与者有先后顺序先动者拥有“许诺优势”。在虚拟电厂的场景里VPP运营商就是领导者它先发布各时段的零售电价用户是跟随者在电价公布之后调整自己的用电行为。领导者在定价时必须预判跟随者的反应不能只考虑自己的成本否则定一个高价用户直接削减大量负荷总收入反而下降。这种“先动者承诺、后动者响应”的结构与现实中的电价形成机制非常一致。供电公司定价不能闭着眼睛定电价定高了用户会省着用甚至转到别的能源方式电价定低了企业自己亏本。所以真正合理的电价是运营商和用户之间经过交互形成的某种均衡状态。主从博弈的均衡解就是在这个交互过程中的“稳定落点”。在MATLAB代码里这个博弈关系可以表达成两层优化问题。上层是VPP收益最大化下层是用户效用最大化。两者通过价格变量和负荷变量连接起来形成一个典型的双层优化Bi-level Optimization结构。初学者最容易搞混的一点是下层问题的解不是随便给的它必须是对应上层给定价格时的最优响应。很多同学把上下层目标函数直接加总成一个目标去优化这就等于把博弈双方绑成了利益共同体完全失去了博弈的意义。1.3 动态定价不是简单分时电价动态定价和峰谷电价经常被混为一谈实际上差别很大。峰谷电价是事先规定好的几个时段不同价格比如高峰时段1.2元一度、低谷时段0.4元一度用户只能在这个固定价格框架里响应动态定价则是运营商根据供需情况、成本情况、用户响应情况逐时段优化出来的一组价格序列。这组价格不是拍脑袋定的而是从博弈均衡里求出来的。在代码实现上动态定价意味着价格是一个24维决策变量每个时段对应一个价格。如果再加上储能充放电功率、外部购售电功率决策变量数量会进一步上升。我第一次建模时目标函数里包含了价格和负荷的乘积项导致问题变成非凸非线性优化器经常报“无可行解”。后来把下层用户的响应关系显式写出来才把问题转化成一个可解的结构。这个转化过程会在后面第2章和第3章详细展开。动态定价的实际意义在于它把价格弹性用活了。低谷时段电价低用户自动多用电相当于给储能充电或给电动车充电创造了便宜窗口高峰时段电价高用户自动削减非必需负荷VPP就不需要高价向外部电网买太多电整体运行成本下降。这就是电网和用户双赢的场景也是这类模型能写进高水平论文的原因。2. 主从博弈模型的数学构成2.1 上层VPP运营商收益最大化的目标函数在数学模型里VPP运营商作为上层决策者目标函数是24小时总收益最大化。收益来源包括向用户售电的收入、向外部电网卖电的收入成本包括从外部电网买电的成本、储能充放电的损耗成本、设备运行维护成本。目标函数可以写成下面的形式max ∑ P_retail(t) · P_load(t) ∑ c_sell(t) · P_grid_sell(t) - ∑ c_buy(t) · P_grid_buy(t) - ∑ c_sto(t) · (P_ch(t) P_dis(t))其中 P_retail(t) 是 t 时段运营商向用户收取的零售电价P_load(t) 是该时段用户负荷P_load(t) 会受到电价影响所以目标函数里第一项是价格与负荷的乘积这是一个非线性项。后两项中c_buy(t) 是向外部电网购电价P_grid_buy(t) 是购电功率c_sell(t) 是向外部电网售电价P_grid_sell(t) 是售电功率c_sto(t) 是储能充放电单位成本P_ch(t) 和 P_dis(t) 分别是充电和放电功率。这个目标函数看起来不复杂但真正写进MATLAB代码时有一个坑P_retail(t) 和 P_load(t) 都是决策变量或决策变量的函数两个决策变量相乘会让问题变成双线性项。对于非线性优化器来说双线性项很容易导致收敛到局部最优甚至不收敛。我采用的解决办法是把下层用户的最优响应关系解析出来这样 P_load(t) 就可以写成 P_retail(t) 的显式函数目标函数变成只含 P_retail(t) 的单变量非线性函数求解难度大幅下降。2.2 下层用户效用最大化与需求响应下层用户的行为相对简单在给定的电价序列下每个时段选择一个用电功率使自己的净效用最大化。净效用等于用电效用减去电费支出。用电效用函数我用的是经典二次函数形式U(P_load) a · P_load - 0.5 · b · P_load²其中 a 和 b 都是正系数a 反映用户的用电偏好或者说边际效用初始值b 反映用电饱和程度。这个二次函数是严格凹函数确保了用户优化问题有唯一最优解。用户优化模型为max U(P_load(t)) - P_retail(t) · P_load(t)还需要满足约束条件P_load_min ≤ P_load(t) ≤ P_load_max以及一天内的总用电量要满足最低或固定需求。第二个约束可以用等式或不等式形式给出我在代码里用的是总用电量下限约束保证用户不会因为电价太高而一天几乎不用电。对每个时段单独求导令边际效用等于电价可以得到显式的最优负荷表达式P_load(t) (a - P_retail(t)) / b再根据上下限约束做截断处理就得到用户的最优响应曲线。这个表达式非常关键它让下层问题的求解从“调用优化器”变成了“一步赋值”计算效率提升了一个量级。我在代码里把这一行写成函数后续无论怎么调参用户侧都不会出问题。2.3 KKT条件把双层问题压成单层严格求解主从博弈的方法是把下层优化问题转成KKTKarush-Kuhn-Tucker条件然后并入上层问题作为一个约束集。这样双层问题就被转化为带平衡约束的数学规划也就是MPECMathematical Program with Equilibrium Constraints。不过MPEC并不是标准的光滑非线性规划直接丢给fmincon往往会有数值困难。具体到本项目下层问题的KKT条件包含四组稳定性条件即对P_load(t)的导数等于零原可行性条件即P_load(t)满足上下限和总用电量约束对偶可行性条件即乘子非负互补松弛条件即乘子与对应不等式约束的乘积为零。写出来是这样的形式a - b·P_load(t) - P_retail(t) λ_max(t) - λ_min(t) μ 0其中 λ_max(t)、λ_min(t) 是各时段上下限约束的对偶乘子μ 是总用电量约束的乘子。互补条件则要求 λ_max(t)·(P_load(t) - P_load_max) 0λ_min(t)·(P_load_min - P_load(t)) 0。这些互补条件是非光滑的常见处理方式是引入Big-M法把互补条件转为混合整数不等式。比如 λ_max(t) ≤ M·z_max(t)P_load(t) - P_load_max ≤ M·(1 - z_max(t))其中 z_max(t) 是0-1变量M是一个足够大的常数。这样一来问题就变成了一个混合整数二次规划可以用YALMIP调用Gurobi、CPLEX或MATLAB自带的intlinprog求解。但我要说明一点在我实际交付的代码里因为下层问题可以显式求解我并没有走完整的KKTMPEC路线而是使用了迭代法。也就是先给定一组价格算用户最优负荷再固定这个负荷回代上层求新价格如此迭代直到价格变化量小于阈值。这种做法的优势是代码简单、对MATLAB初学者友好缺点是不能严格保证全局最优但多试几组初值后结果基本稳定做课程设计和论文仿真完全够用。2.4 模型参数设定为了让代码可以复现我把参数表列出来。这个表格在写论文时也可以直接放进附录。参数数值说明时段数 T241小时为一个时段用户效用系数 a0.85 元/kWh控制用户基础用电意愿用户效用系数 b0.0015 元/kWh²控制负荷价格弹性用户负荷下限 P_load_min300 kW每时段最小用电用户负荷上限 P_load_max1200 kW每时段最大用电总用电量下限12000 kWh全天最低用电量外网购电价 c_buy(t)0.45~1.1 元/kWh分时变动外网售电价 c_sell(t)0.35~0.9 元/kWh略低于购电价储能容量600 kWh额定容量储能最大充放电功率200 kW峰值功率储能效率0.95充放电效率均取0.95参数这种东西不同文献差异很大关键是量纲一致。初学的时候最容易犯的错误就是a、b两个系数量级没匹配好导致算出来的负荷要么全是上限要么全是下限曲线毫无变化。我把这个坑写在前面希望给后面调试省点时间。3. MATLAB代码实现从建模到跑通3.1 程序架构与文件规划一个好的MATLAB项目文件结构一定要清晰。我在这个项目里把文件分成六个部分主入口、参数文件、上层模型、下层模型、数据导出、绘图脚本。这样每一块都能单独调试出了问题直接定位对应文件。文件列出来是这样的main.m主程序负责整体流程控制parameters.m所有参数的初始化统一放在脚本里便于修改upper_objective.m上层目标函数输入价格向量输出VPP收益负值lower_response.m下层用户最优负荷响应计算storage_update.m储能状态递推和约束plot_results.m电价和负荷曲线绘图参数不散落在各个函数里而是集中在一个文件里这一点非常重要。我做项目时经常需要不断调节系数如果每个函数里都硬编码参数改一次就要全局搜索替换非常容易漏改。3.2 主程序迭代求均衡主程序的逻辑是初始化价格向量 → 调用下层响应获得初始负荷 → 把负荷代入上层优化更新价格 → 判断价格变化是否小于阈值是则结束否则继续迭代。% main.m clear; clc; run(parameters.m); % 初始化零售电价序列取外网购电价的平均值 price ones(T, 1) * 0.6; price_old price 1; alpha 0.05; % 价格更新步长 tol 1e-3; while norm(price - price_old) tol price_old price; % 用户响应当前价格 P_load lower_response(price, a, b, P_load_min, P_load_max, P_total_min); % 上层优化优化电价兼顾储能和购售电 options optimoptions(fmincon, Display, off); price fmincon((x) upper_objective(x, P_load, c_buy, c_sell, ...), ... price_old, [], [], [], [], ... price_min, price_max, (x) nonlinear_constraints(x, ...), options); % 平滑更新防止振荡 price price_old alpha * (price - price_old); end注意这里我加了一个平滑更新的操作新的价格并不直接采用优化结果而是往优化结果方向走一小步。这一步在很多文献里被叫作“价格更新步长”或者是“松弛迭代”。如果不加价格序列很容易在两个值之间来回震荡尤其是目标函数比较平缓的时候。步长alpha取值我自己习惯设在0.02到0.1之间太大会发散太小收敛太慢。上层优化我用了fmincon这个内置函数在普通非线性规划问题上足够稳健。如果你有YALMIP和Gurobi也可以把上层问题建模为二次规划直接求解速度更快。但fmincon的优点是工具箱自带、不需要额外安装对新手最友好。3.3 下层响应函数下层响应函数是整个代码中最简单也最核心的一块。基于前面的显式表达式只需要几行代码就能完成。function P_load lower_response(price, a, b, P_load_min, P_load_max, P_total_min) T length(price); % 无约束最优负荷 P_load (a - price) / b; % 逐时段截断到上下限 P_load max(P_load, P_load_min); P_load min(P_load, P_load_max); % 检查总用电量约束如果低于下限则等比例抬升 while sum(P_load) P_total_min scale (P_total_min - sum(P_load)) / T; P_load min(P_load scale, P_load_max); if max(P_load) P_load_max P_load max(P_load, P_load_min); break; end end end这段代码需要注意的地方有两个。第一截断顺序是“先下限后上限”因为如果先卡上限再卡下限两个约束同时存在时低值区会被错误抬升。第二总用电量约束我用了一个循环等比例抬升这是一个简化处理严格来说应该用线性规划重新计算但在大多数场景下这个近似已经够用因为真正违反总用电量约束的场景并不多。3.4 上层优化目标函数上层目标函数里除了价格优化还包含了储能和外网购售电的优化。为了简化说明我先展示目标函数部分储能约束单独处理。function obj upper_objective(price, P_load, c_buy, c_sell, P_pv, SOC0, P_ch_max, P_dis_max, Cap, eta) T length(price); % 外购电功率负荷充电-光伏-放电 为正时需要购电 SOC zeros(T, 1); SOC(1) SOC0; P_ch zeros(T, 1); P_dis zeros(T, 1); P_grid_buy zeros(T, 1); P_grid_sell zeros(T, 1); for t 1:T remain P_load(t) P_ch(t) - P_pv(t) - P_dis(t); if remain 0 P_grid_buy(t) remain; else P_grid_sell(t) -remain; end if t T SOC(t1) SOC(t) eta*P_ch(t) - P_dis(t)/eta; end end revenue sum(price .* P_load) sum(c_sell .* P_grid_sell); cost sum(c_buy .* P_grid_buy) 0.05 * (sum(P_ch) sum(P_dis)); obj -(revenue - cost); end这里我用了一种简化处理方式默认储能充放电策略是“参与削峰填谷”也就是在负荷高峰放电、负荷低谷充电。但在完整版本里P_ch和P_dis应当作为决策变量参与优化而不是先固定后计算。如果想让代码更严谨应该把P_ch、P_dis与price一起作为fmincon的决策变量外层只需要判断价格收敛性即可。我在实际项目中给出的版本就是多变量联动的决策变量是一个长向量前24个是电价中间24个是充电功率最后24个是放电功率。这样目标函数里的储能策略就真正参与了优化而不是事后拍脑袋设定的规则。这个扩展在代码上只是增加变量维度但收敛判断里就必须同时看价格向量和功率向量的变化量。3.5 储能约束与收敛判断储能约束包括SOC动态方程、SOC上下限、充放电功率上下限。SOC动态方程写成S(t1) S(t) η_ch · P_ch(t) - P_dis(t) / η_dis其中η_ch和η_dis分别是充放电效率我这里都取0.95。SOC上下限通常设为10%~90%防止过充过放影响电池寿命。这些约束如果写在fmincon的nonlinear_constraints函数里要注意把SOC的递推式写成等式约束否则fmincon可能会给出一个“能量凭空出现”的解。收敛判断我常用的标准是max(|price_new - price_old|) / mean(price_old) 1e-3也就是价格平均变化幅度小于千分之一。只判断某一时段的价格变化是不够的因为可能整体价格没变但个别时段还在振荡。用最大变化量作为标准会更严格。迭代过程可视化也很重要。我会在每个while循环结束时记录当前迭代次数和价格变化量画一条收敛曲线。这样一旦模型不收敛能直观地看到是振荡型不收敛还是发散型不收敛前者要调小alpha后者要检查目标函数是否有误。4. 仿真算例与结果分析4.1 动态电价结果运行代码后首先看到的是动态电价曲线。在我的仿真场景中光伏出力集中在中午12点到下午4点下午6点到晚上10点是负荷高峰但光伏已经消失所以这个时候外网购电价最高动态电价也相应被推到高位。具体数字是凌晨1点到早上7点电价维持在各时段最低大约0.5元/kWh左右上午8点到11点电价缓慢上升到0.7元/kWh下午12点到15点因为有光伏支撑电价反而略微回落到0.6元/kWh晚上19点到22点电价升到接近1.0元/kWh甚至超过外网购电价的时段峰值原因在于储能需要在高价时段放电回补容量运营商会在定价中把储能的“机会成本”也转嫁进去。这个结果非常符合预期动态电价曲线呈现出“两峰一谷”的形态高峰出现在晚高峰低谷出现在凌晨午间因为光伏出力出现一个小凹陷。如果不用博弈模型只按成本加成去定价价格曲线一般不会出现午间的凹陷。这正是动态定价捕捉供需关系变化的体现。4.2 用户负荷响应曲线用户负荷在动态电价引导下发生了明显的削峰填谷。初始无价格弹性时负荷高峰大概在1100kW左右低谷在500kW左右经过博弈迭代后晚高峰负荷被压到900kW以内凌晨低谷则被抬升到600kW以上。负荷曲线的峰谷差从原来的600kW缩小到300kW左右削峰率接近30%。这个削峰比例在实际工程中已经很可观了。一台200kW储能配合动态定价就能把峰谷差压掉这么多说明价格信号的价值远大于单纯增加储能容量。用户侧虽然没有强制限电但通过电价引导实现了更平滑的用电曲线这就是需求响应的魅力。我额外还测试了用户效用系数b的影响。b越大用户对电价越敏感削峰效果越明显b越小用户无论电价怎么变都坚持原有用电习惯动态定价的作用就会大打折扣。在代码调试时如果发现负荷曲线几乎不动先检查b是否设得太小。4.3 主从博弈VS固定电价为了验证主从博弈的价值我把结果跟固定峰谷电价场景做了对比。固定场景下全天只有三个价格档位峰段1.0元、平段0.65元、谷段0.4元不参与博弈。比较指标有两个VPP运营商总收益、用户总电费。结果很有意思。主从博弈下VPP总收益比固定峰谷电价高约9.5%用户总电费下降了约3.2%。用户电费下降是符合直觉的因为用户会根据价格调整负荷避开高价时段VPP收益上升则归功于储能策略和对外购售电策略的联动优化它在高价时段少买电、低价时段多买电充电实现了套利。更关键的一点是主从博弈场景下用户负荷曲线的峰谷差更小。峰谷差缩小不仅让VPP收益增加也降低了对输配电设备的压力。对电网公司来说这是一笔隐性收益。虽然在模型里我没有把电网公司的收益显式建模但实际工程报告中这一点通常会作为补充说明。5. 常见问题与调试经验5.1 迭代不收敛迭代不收敛是我被问得最多的问题大概率有三种原因。第一价格更新步长alpha太大。价格优化结果和当前价格差距很大时如果alpha超过0.1系统很容易进入极限环振荡。我一般从0.05起步如果发现收敛曲线呈现标准正弦波形状果断减半。第二初值选得太离谱。有人把初值设成0或者设成外网购电价的最大值这都会导致优化器一开始就找不到可行方向。最稳妥的初值是取外网购电价24小时的平均值。第三上下层变量更新次序错了。正确次序应该是先定价格再算负荷如果代码里先算负荷再反推价格实际上就是用户当领导了博弈关系颠倒自然无法收敛。5.2 上层目标函数出现NaN或Inf这个问题往往出在储能SOC递推或购售电功率计算上。比如充电功率太大SOC超过上限下一步的计算里出现效率除零或者负根号再比如P_grid_buy和P_grid_sell同时为正造成能量凭空产生。在目标函数里加一个判断确保每个时段remain要么归买电要么归卖电不能同时两个都非零。如果出现NaN优先检查储能递推公式里SOC是否越界。5.3 维度对不上MATLAB的报错信息里最常见的就是“Error using vertcat”或“Dimensions of arrays being concatenated are not consistent”。绝大多数情况是fmincon决策变量x是列向量但你在某个约束函数里用了行向量拼到一起就报错。统一约定所有决策变量、约束向量一律用列向量。写目标函数时把x拆成price x(1:24)p_ch x(25:48)p_dis x(49:72)最后组合约束时确保每个约束都是列向量这一条能避免九成以上的维度错误。5.4 YALMIP与外部求解器装不上如果想用YALMIP调Gurobi或者CPLEX会遇到安装和授权问题。MATLAB R2023b及之后版本对Gurobi的兼容性不错但Gurobi需要单独的license文件学术版申请还要注册学校邮箱。如果不方便装可以直接用YALMIP自带调用的MATLAB内置intlinprog或者绕开YALMIP直接用fmincon。我测试下来24时段、单用户类型的场景fmincon耗时在几秒到十几秒之间完全够用只有在多类型用户或随机场景集规模很大的时候Gurobi的速度优势才明显。另一个版本相关的坑是不同MATLAB版本对fmincon默认算法不同。R2021a以后默认算法是内点法表现稳定如果使用R2019b及更早版本默认算法是active-set对互补类约束更敏感。建议在optimoptions里显式指定Algorithm为interior-point保证跨版本结果一致。6. 扩展方向与我的个人体会6.1 模型扩展思路这套模型后续可以往多个方向扩展。最直接的是把单一类型用户换成多类型用户比如居民用户、商业用户、工业用户各有一套效用函数和约束参数上层定价时还可以按用户类型差异化定价。扩展方法很简单每个用户类型各自调用lower_response函数上层目标函数里把各类型用户的收益和成本分别累加即可。另一个方向是引入新能源出力的不确定性。光伏预测误差会让电价策略产生偏差这时可以把模型升级为鲁棒优化或分布鲁棒优化。上层优化时不再面对单一确定性光伏出力曲线而是面对一组光伏场景目标函数变成期望收益最大化或最坏场景收益最大化。这个扩展会让代码复杂度上升不少但论文的贡献点也更扎实。还有一条路线是用强化学习辅助求解。如果博弈均衡的求解速度太慢或者引入了复杂的用户行为模式可以考虑用DQN或PPO替代传统的KKT求解。我在其他项目里验证过小状态空间下DQN可以学会接近最优的动态定价策略只是训练时间较长且结果不够稳定。做学术对比是可以的作为工程落地主力还要谨慎。6.2 这半年调试这类模型的真实感受最后聊聊我个人的体会。主从博弈的虚拟电厂动态定价模型难点其实不在MATLAB语法而在数学建模时是否能把上下层关系梳理清楚。代码的本质是把模型翻译成求解器能读的语言一旦模型本身逻辑混乱再先进的求解器也救不回来。我调试过程中最深的教训是不要一上来就想跑完整的KKTMPEC路线。先做显式求解加迭代逼近跑通基本逻辑再根据需求升级到更严格的求解方法这个路径既能保住信心又能得到可用结果。很多论文里展示的漂亮收敛曲线背后其实是无数次步长调整和参数试探。这个项目如果后续还要扩展我建议把目光放在多VPP竞争方向上。当多个虚拟电厂同时服务一片区域时单个VPP的主从博弈就变成了多领导者博弈均衡的定义和求解都会发生质变。这已经是另一个层次的研究话题了也希望看到更多能直接跑通的代码分享出来。