微网里的能量博弈可比菜市场砍价刺激多了——菜市场砍价最多是一买一卖两头拉锯,微网里却是零售商、产消者、光伏、储能、不确定的天气和看不见的市场风险挤在同一张结算表里角力。前阵子我啃了一套基于CVaR风险管理的动态定价MATLAB代码,折腾了差不多两个星期才把逻辑理顺。今天不聊虚的,直接从利益格局、风险度量、数学模型到代码实现,把整套思路拆开讲清楚,面向想用MATLAB在能源交易场景里做仿真决策的人,尤其是那些在研究微网优化调度、需求响应定价方向卡在建模和编程之间的同行。
1. 先搞懂利益格局:零售商和产消者在争什么
1.1 这盘棋里每个玩家想要什么
微网里有两类核心角色,一个叫零售商,一个叫产消者。零售商从电网或者上级市场批发电能,再卖给微网里的用户;产消者更特殊,家里装了光伏、储能甚至小风机,既能用电也能发电,多出来的电可以卖给零售商。听起来很简单,但一旦引入动态电价,双方就进入了持续的博弈状态。
零售商想的是:我定的电价要能让产消者在用电高峰少用电、在光伏大发时段多用电,这样能减少我从外部市场高价买电的冲击。同时我还得保证自己不亏本,毕竟外部市场价格波动、光伏出力不确定性都可能导致成本失控。产消者想的是:电价便宜我就多买电、甚至给储能充电,电价贵我就少买电或者把光伏多余的电卖回给电网,偶尔还能靠储能套个利。两边都有自己的优化目标,而电价就是连接两边决策的桥梁。
这套代码里,零售商和产消者不是同时做决策的,而是按主从博弈来走:零售商先公布一个电价策略,产消者再根据这个电价调整用能计划,零售商看到产消者的反应后重新调整电价,直到双方都不愿意再改变自己的策略。用术语说就是Stackelberg博弈,零售商是领导者,产消者是跟随者。这种结构很贴近现实,因为电价往往是运营商先出,用户后响应。
1.2 为什么动态定价不是拍脑袋
动态定价的核心不是“今天比昨天贵一点”这种傻瓜式调整,而是要让价格能反映电力在不同时间段的真实稀缺程度。你在用电高峰把电价抬高,产消者的负荷曲线自然会被压低,这就相当于用价格杠杆做了需求响应。同时,光伏出力的间歇性会让净负荷变得难以预测,如果电价一直不变,产消者就没有动力在正午光伏过剩时把负荷挪过来,零售商也没法用价格信号去平衡供需。
所以动态定价需要解决一个双层优化问题:上层是零售商的定价模型,目标是自己收益最大,但必须考虑产消者会如何响应;下层是产消者的用能优化模型,目标是自己用电成本最小或收益最大,约束条件包括储能充放电、可平移负荷、光伏出力上限等。上下两层互相嵌套,这就是为什么这套MATLAB代码不是简单单层优化,而是要在迭代中反复调用下层模型。
我最开始犯过一个错误,就是试图把两层模型合并成一个大规模优化问题一次求解,结果约束条件和变量耦合到一起,求解器直接报内存不足。后来改成主从迭代,每轮先固定电价,解一出产消者的最优响应,再根据响应更新零售商的边际收益,这才跑通。很多刚接触微网建模的人也会踩这个坑,建议一开始就按主从结构拆开写,别贪图一步到位。
2. CVaR风险管理:为什么比“平均损失”靠谱
2.1 VaR和CVaR到底差在哪
微网优化不能只看期望收益,因为不确定性意味着你有可能踩中极端场景。比如某天光伏出力大幅低于预测,外部市场价格又突然飙高,零售商当天的利润可能直接被打穿。如果决策目标里只写“期望收益最大化”,那做的其实是一个风险中性决策,对极端情况视而不见。
业界常用的风险度量有两种,VaR和CVaR。VaR是“在给定置信水平下,最坏情况下的最大损失”,比如95%置信水平下VaR等于10万元,意思是只有5%的可能性损失会超过10万元。但是VaR有个明显缺点,它只告诉你分位点在哪,完全不关心超过这个分位点之后损失有多惨。也就是说,VaR忽视了尾部风险。
CVaR把这个问题补上了:它衡量的是“损失超过VaR那一部分的条件期望”。打个比方,如果你开车出门,VaR告诉你“有5%的概率撞车,损失超过修车费”,CVaR则告诉你“如果真撞了,平均下来要花多少钱修”。显然CVaR更符合风险管理的直觉。在微网环境里,光伏骤减、负荷突升、市场电价暴涨这些极端情况都属于尾部分布,你当然希望决策时把这些场景也纳入约束。
2.2 把CVaR翻译成数学公式
假设零售商在某个决策策略下面临一组不确定场景,每个场景s对应一个净收益损失值f(s),场景概率为π(s)。在置信水平β下,CVaR的定义是:
CVaR_β = min { η + (1/(1-β)) * Σ_s π(s) * max(f(s) - η, 0) }
这里η其实是一个辅助变量,它等于VaR的最优近似值。CDaR、CVaR这类公式在金融工程里很常见,搬到微网里也适用,只是把“投资组合损失”换成“零售商日运行成本”或者“负收益”。
在实际MATLAB代码里,我会把目标函数改造成一个带辅助变量η和u_s的最小化问题。u_s是场景s下的超额损失,约束条件要加一个u_s大于等于f(s)减去η,并且u_s大于等于0。这样一来,原来难以处理的分位数问题就被变成了线性规划或者凸优化问题。这也是为什么很多论文里说CVaR可以线性化。
不要被公式吓到。在MATLAB中实现时,你只需要增加一个维度为场景数量的变量u,以及每个场景对应的约束,优化器会自动处理这些线性等式。我写代码时习惯把CVaR部分单独封装成一个函数,输入是各场景下的损失向量和置信水平,输出是CVaR值和对应的辅助变量,这样主程序里就能很干净地把它叠加到目标函数上。
3. 数学建模与MATLAB代码架构
3.1 主从博弈模型:零售商领跑还是跟随?
我在这套代码里用的是零售商作为领导者,产消者作为跟随者的主从博弈。为什么不让产消者当领导者?因为现实里零售商掌握着市场接口和定价权,产消者只能被动接受电价,哪怕嘴里说“再贵我就不买了”,最后还是只能在给定电价下做局部优化。从计算角度看,领导者问题包含跟随者的最优反应函数,通常用KKT条件把下层优化转化为上层约束,然后用求解器直接解MPEC(数学规划均衡约束)问题。但MPEC问题往往非凸,路径依赖严重,直接用现成求解器容易卡壳。
所以代码里我采用了交替迭代法:外层循环更新电价,内层固定电价求解产消者的最优负荷曲线,然后用产消者的需求曲线更新零售商的边际收益函数,再重新计算电价。每一步都只解一个相对简单的优化问题,稳定性和可调试性都更好。数据流大概是这样的:首先是基础参数初始化,包括24小时预测负荷、光伏出力曲线概率场景、储能参数、批发市场价格场景;然后进入迭代;每一轮里,零售商先根据当前价格p和风险偏好给出定价,产消者求解自己的用能优化得到购电功率P_buy;零售商再用P_buy更新自己的收益,检查前后两轮电价差距,如果小于阈值就收敛。
3.2 代码模块划分和数据流
整个项目我按功能拆成五个脚本模块,这样不用每次跑全流程都重新复制一大坨代码。首先是parameters.m,定义所有常量,比如24时段、储能的容量和最大充放电功率、产消者负荷范围、光伏容量、置信水平β、风险厌恶系数λ、场景数N。然后是scenario_generation.m,用来生成光伏和负荷的不确定性场景,我会用拉丁超立方采样先剪出一组初步样本,再通过同步回代削减法把场景数砍到合适的数量,比如从500个降到20个,这样既能保留尾部极端情况又不会让求解慢到离谱。
接下来是retailer_profit.m,它负责计算零售商在给定电价和产消者响应下的收入、购电成本、风险损失。这个函数要能够接受一个电价向量和负荷响应向量,返回期望收益和CVaR损失值。第四个模块是prosumer_response.m,这是下层优化,输入是电价向量,输出是产消者最优购电功率、储能充放电计划和总支出。最后一个模块是main_iteration.m,它把上面几个模块串起来,实行双层迭代,输出最终电价曲线和各时段功率曲线。
实际跑起来你会发现,光有模块还不够,变量结构组织不好照样乱套。我习惯用结构体struct保存所有中间结果,比如sceanario.pv保存光伏场景矩阵,result.price保存迭代过程中的电价历史,这样后面做图和调试都能快速定位。
3.3 核心函数:目标函数+约束条件怎么写
在MATLAB里,我们通常用fmincon或quadprog来解这些优化问题。零售商的目标函数是求最小化以下综合指标:
G = (1-λ) * ExpectedCost - λ * CVaR
其中ExpectedCost是零售商的期望净成本,λ是风险厌恶系数,0到1之间。λ=0就是完全风险中性,只看平均成本;λ=1就只关心最坏情况下的尾部平均损失,会变得极度保守。这个写法把风险偏好直接变成了一个标量权重,特别适合做敏感性分析。
细分到具体代码,零售商决策变量是卖给产消者的电价p(t),可能还要加上从主网购电的功率P_grid(t)。约束条件包括电价上下限、功率平衡方程、零售商的购电成本构成等。产消者的下层优化里,决策变量包括储能充电功率P_ch(t)、放电功率P_dis(t)、可平移负荷的状态,以及向电网购电功率P_buy(t)。约束条件会包含储能SOC的状态转移、功率约束、负荷平衡约束等。
这里给你一段产消者响应函数的伪结构代码,我用MATLAB的optimvar和prob2struct来做变量声明,这样后续换成其他求解器也方便:
function [res] = prosumer_response(p, param) N = param.N; % 时段数 prob = optimproblem('ObjectiveSense','min'); P_buy = optimvar('P_buy', N, 1, 'LowerBound', 0, 'UpperBound', param.PbuyMax); P_ch = optimvar('P_ch', N, 1, 'LowerBound', 0, 'UpperBound', param.PchMax); P_dis = optimvar('P_dis', N, 1, 'LowerBound', 0, 'UpperBound', param.PdisMax); SOC = optimvar('SOC', N+1, 1, 'LowerBound', param.SOCmin, 'UpperBound', param.SOCmax); % 目标:购电成本最小 prob.Objective = sum(p .* P_buy) * param.dt; % 约束:负荷平衡、储能SOC更新 prob.Constraints.energy = P_buy + P_dis - P_ch == param.Load - param.PV; prob.Constraints.soc_update = SOC(2:N+1) == SOC(1:N) + param.dt * (param.eta_ch * P_ch - P_dis / param.eta_dis); prob.Constraints.soc_boundary = SOC(1) == param.SOC0; prob.Constraints.soc_final = SOC(N+1) == param.SOC0; % 求解 opts = optimoptions('quadprog','Display','off'); [sol, fval, exitflag] = solve(prob, 'Options', opts); if exitflag <= 0 error('下层优化求解失败'); end res.P_buy = sol.P_buy; res.P_dis = sol.P_dis; res.P_ch = sol.P_ch; res.SOC = sol.SOC; res.Cost = fval; end注意这里我把约束写成等式,是因为简化模型假定微网内功率平衡时所有缺额都需要从零售商处购买,多出来的电也只能卖给零售商。如果后续要加入向主网卖电,那就需要额外变量和对应价格,模型会更复杂,但基本骨架不变。
4. 动态定价求解实操:从标称场景到风险场景
4.1 场景生成:光伏出力和负荷的不确定性怎么模拟
前面说了CVaR需要一组场景,场景的质量直接决定风险度量的可信度。光伏出力我一般会用一个基础预测曲线加上随机扰动来生成,扰动来源包括云层波动、温度影响、组件衰减等。更讲究一点的做法是用Beta分布来拟合光伏出力的归一化值,因为Beta分布定义在0到1之间,和辐照度归一化后的分布非常契合。负荷方面,则会对每个时段的基准负荷乘以一个正态扰动项,扰动项均值是1,标准差根据历史数据设定,比如早晚高峰时段标准差大,夜间低谷时段标准差小。
场景生成完之后还要削减,否则N=500个场景会让每个时期的求解都慢得让人抓狂。同步回代削减法的核心思想是:计算每对场景之间的欧氏距离,把距离最近的场景合并,并更新概率权重,直到场景数降到目标值。我做完之后通常会做一个验证,检查削减后的场景集合下,期望负荷和原始场景的期望负荷差距控制在2%以内,否则就增大保留场景数。
4.2 迭代求解与参数更新
主从迭代的收敛判据我取的是电价向量的二范数相对变化小于1e-4。每轮迭代里,零售商会根据当前产消者响应后的净负荷,重新计算边际成本和风险项,然后更新电价。为了避免电价在迭代初期来回震荡,我加了一个阻尼更新策略,新电价等于旧电价加上阻尼系数乘以优化方向,阻尼系数通常取0.3到0.5。千万不要在没加阻尼的情况下一路用最新解覆盖,我试过,电价会在两个极端之间反复跳,根本收敛不了。
另一个关键点是零售商优化内部的CVaR约束。要注意,CVaR只是目标函数里的一项,你还得把场景逻辑传到下层优化里。有些实现会把零售商的购电成本建成一个场景相关的函数,但产消者在此轮看到的电价依然是单一确定性价格,也就是说产消者不需要关心场景,它只对当前公布的电价做响应。这样做的好处是两层的输入输出维度很清晰,外层负责处理所有不确定性,内层只需要解确定性优化,避免了场景维度的嵌套爆炸。
每完成一轮迭代,我会保存本轮的电价曲线、产消者购电功率、储能SOC、零售商的期望成本和CVaR值,存到一个日志表里。之后可以把迭代过程中电价的演化画成三维图或热力图,观察价格如何从初始值逐渐稳定到博弈均衡解。我见过很多代码跑完只输出一个终值,中间过程全丢,这样出问题很难排查,建议无论如何都要把迭代历史存下来。
4.3 结果输出和图谱分析
最终结果按照惯例会生成几张核心图,第一张是24小时电价曲线,第二张是产消者购电功率与光伏出力、负荷曲线的对比,第三张是储能SOC变化,第四张是不同风险厌恶系数下的高效前沿图。高效前沿图特别重要,它能直观展示风险和收益之间的权衡关系,横轴是CVaR损失,纵轴是期望收益。你可以画一组不同λ下的解,得到一条向右下方倾斜的曲线,那就是风险收益帕累托前沿。这比单纯列一个数字有说服力得多。
我跑出来的结果里有一个典型现象:当λ从0升到1,零售商的电价波动会变小,更倾向于用稳定价格换取可预测的负荷响应,而产消者的购电成本波动也随之下降。这就说明风险厌恶系数不只是数学里的一个权重,它真的会改变市场参与者的行为。如果你在汇报里只讲均价和总收益,听众可能无感,但把不同λ下的电价曲线画在同一个坐标系里,所有人立刻能看出风险管理带来的“削峰”效果。
5. 踩坑实录与调参心得
5.1 常见报错与解决方案
这个项目里最常见的报错出现在求解器层面。新手最常遇到的是fmincon给出的exitflag = -2,这意味着模型没有可行解。大部分原因是约束冲突,比如储能SOC终值必须等于初值,同时某个时段又要同时充电放电和满足过高的负荷,结果功率上限不够。排查方法很简单,先把SOC终值约束放松,看看是否还有其他冲突,或者把最大放电功率临时调大,验证是不是功率边界卡住了。
还有一种报错是prob2struct生成的A矩阵条件数过大,导致求解器数值不稳定。出现这个问题的场景通常是变量单位差异太大,比如电价在小数点后两位的级别,而储能功率可能到几百千瓦,目标函数里两者乘积的量级悬殊,最终让Hessian矩阵病态。解决办法是对变量做归一化,比如把功率统一除以基准功率,电价统一除以基准电价,然后再求解。
5.2 风险厌恶系数怎么选
λ的选取没有绝对标准,更多看你的业务诉求。零售商的资金实力越弱,越需要规避尾部风险,λ可以取0.6到0.9;如果零售商背后有强大资本支撑,能扛住偶发亏损,那λ可以取0.1到0.3。我个人的习惯是跑一组λ从0到1的灵敏度分析,步长0.1,画出收益和CVaR的前沿,再去决策层挑一个“膝盖点”,也就是曲线曲率最大处的λ。这个点和经济学里的边际效用递减有点像,选在那个位置通常性价比最高。
实际执行时,当你调整λ,CVaR部分的权重变化可能会让整体目标函数数值范围发生跳动。为了避免求解器因为量级变化而迷失方向,我会把期望成本项除以基准成本,把CVaR项除以一个基准风险值,再做加权,这样无论λ怎么变,目标函数的数值范围都不会剧烈波动。
5.3 性能优化建议
这套代码如果场景数很大,迭代次数又多,跑起来可能让人等到怀疑人生。我曾经试过直接上200个场景加100轮迭代,结果一次实验要跑将近一小时。后来我做了三板斧优化。第一板斧是场景削减,把200个场景削减到20个,CVaR值变化不到5%,时间缩短到原来的六分之一。第二板斧是热启动,每一轮迭代里,产消者下层优化的初始值用上一轮的最优解。
第三板斧是并行计算,因为多个场景的损失计算彼此独立,我用parfor来并行评估场景损失,特别是外层零售商评估CVaR时要遍历每个场景,这对并行特别友好。我提醒一句,用parfor之前注意变量分类,避免不同worker之间频繁广播大数据矩阵,否则通信开销比计算还大。
提示:场景削减不是只追求数量少,要看削减后风险指标是否稳定。我每次削减后都会对比一下完整场景和削减场景下的CVaR数值,偏差控制在可接受范围才继续往下跑。
我在实际使用中发现,这套模型虽然名字叫“微网能量博弈”,但它的思想完全可以迁移到虚拟电厂日前申报、社区级充电桩聚合商定价、园区型多能互补系统里。核心思想都是“价格信号引导灵活资源”加“CVaR控住尾部风险”。如果你要在自己的项目里复现,不用急着一次性把代码写完,先把主从迭代的框架搭出来,再用最简单的确定性场景跑通一条曲线,最后再引入CVaR和不确定性场景,每一步都确认结果合理再向上叠,比一口气写完全部逻辑再调bug要高效得多。我自己就因为在刚开始叠了太多功能,导致后面调试时根本分不清是博弈不收敛、场景问题还是CVaR约束写错。一点一点来,反而快。