基于CVaR的微网动态定价与调度策略及Matlab实现
做微网优化的人应该都遇到过这种纠结光伏出力飘忽不定批发市场电价上蹿下跳靠期望值做出来的调度方案看着利润挺高但一遇到极端场景就“翻车”不是购电成本爆表就是被考核罚款。这两年“基于条件风险价值CVaR的微网动态定价与调度策略”成了热门研究方向就是因为CVaR能把这种“尾部风险”量化出来让决策者在利润和风险之间找一个自己舒服的位置。这篇文章我想把这套方法从原理到Matlab实现完整梳理一遍包括CVaR怎么建模、动态定价和调度怎么联动、案例代码怎么搭以及我实际调参踩过的坑。适合正在做微网优化调度、电力市场策略、需求响应方向的研究生和工程师参考。我要讲清楚的不只是公式更是“为什么这么做”以及“跑代码时什么地方最容易卡住”。1. 微网动态定价与调度策略的整体思路1.1 为什么要做动态定价与协同调度微网本质上是一个小型发配用电系统里面有光伏、储能、柴油机或燃气轮机还带着一堆用户负荷。传统做法是把用户负荷当成固定值微网只管“发够电、保平衡”这是单边决策。但现实中用户会对电价做出反应——电价高的时候少用点电价低的时候多安排用电。既然用户有弹性那定价就不该是静态的应该跟着供需关系走这就是动态定价的价值。更进一步定价和调度不是两件事。你把电价定高了用户负荷降下来那发电侧就要相应调整出力反过来储能什么时候充、什么时候放直接影响你向主网购电的曲线也影响你给用户报什么价。所以必须把定价和调度放进同一个优化框架里一起算否则两边各自为政算出的结果在工程上根本没法用。这也是“动态定价与调度协同优化”这几年被反复研究的原因。1.2 CVaR到底解决了什么问题CVaRConditional Value at Risk条件风险价值最早是金融领域用来衡量投资组合尾部风险的指标。什么叫尾部风险就是那些概率不高、但一旦发生就损失很大的情景。在微网场景里光伏出力的预测误差和批发市场电价的波动就是典型的不确定源。我用一个例子说明假设某个调度方案在100个随机场景里有99个场景能赚5000块但有1个极端场景会亏3万。如果只看期望值这个方案利润不错但那个3万的亏损可能就是运营方承受不起的。VaR风险价值告诉你在95%置信水平下最大损失是多少但它不关心超过这个阈值之后损失有多大。CVaR则更进一步它计算的是“一旦损失超过VaR平均会损失多少”。所以CVaR衡量的是“最坏情况下的平均损失”这正好是微网运营者做决策时真正关心的事。在优化模型里我们通常把CVaR作为目标函数里的一项或者把它作为约束条件加进去。这样得到的最优解不是“期望利润最高”的方案而是“在风险可控的前提下利润最高”的方案。说白了CVaR给决策者提供了一个风险偏好旋钮——你越厌恶风险就把权重调大系统就会主动选择更保守的调度策略。1.3 方案框架从不确定场景到优化决策这套策略的完整实现链路分四步不确定性建模对光伏出力和市场电价分别建随机模型用蒙特卡洛抽样生成大量场景代表未来的各种可能。场景处理场景数量直接影响求解规模需要用场景缩减方法快速前向选择法把几千个场景缩到几十个代表性场景同时保持概率分布特征。优化建模把动态定价、机组调度、储能充放电和CVaR风险约束写进一个混合整数线性规划模型。求解分析用Matlab调用求解器Gurobi、CPLEX等求解再做不同风险权重下的对比分析。这个框架的好处是模块化。你替换掉场景生成的方式或者换一种需求响应模型其他模块基本不用大改。我后面讲Matlab实现时也会按这个模块来拆分代码。2. 核心数学模型CVaR的建模原理与线性化2.1 从VaR到CVaR一次说清两个指标的区别先给出数学定义这个搞清楚后面建模才不迷糊。假设微网在某个调度方案下每个随机场景得到净收益或运行成本我们把损失定义为负收益。设损失随机变量为(L)置信水平为(\beta)常见取0.9、0.95那么VaR损失分布 (\beta) 分位数即 (VaR_\beta(L)\inf{l\mid P(L\le l)\ge\beta})。它回答“最坏情况下有(\beta)的概率损失不超过多少”。CVaR损失超过VaR部分的期望即 (CVaR_\beta(L) E[L \mid L \ge VaR_\beta(L)])。它回答“如果灾难真的发生平均损失是多少”。两者关系可以用一句话概括VaR是门槛CVaR是跨过门槛之后的平均水深。VaR有个致命缺点是不满足次可加性也就是说组合后的风险可能看起来比各部分风险之和还大这在优化里会造成“把鸡蛋分到多个篮子反而被认为风险更高”的悖论。CVaR是一致性风险度量满足次可加性更适合作为优化目标或约束。2.2 基于情景法的CVaR线性化表达直接算CVaR需要知道损失分布的分位数这没法直接塞进线性优化模型。但Rockafellar和Uryasev在2000年给出一个著名的变换公式把CVaR的计算转化为一个线性规划问题[ \min_{\alpha} \left{ \alpha \frac{1}{(1-\beta)} \sum_{s1}^{S} p_s \cdot \max(L_s - \alpha, 0) \right} ]其中 (S) 是场景总数(p_s) 是第(s)个场景的概率等概率抽样时 (p_s1/S)(L_s) 是场景(s)下的损失(\alpha) 就是VaR值优化出来的目标函数值就是CVaR。这个公式里有一个 (\max) 函数它不是线性的但可以引入辅助变量 (z_s) 把它线性化[ z_s \ge L_s - \alpha, \quad z_s \ge 0 ]然后把目标函数写成[ \min \left{ \alpha \frac{1}{(1-\beta)S} \sum_{s1}^{S} z_s \right} ]这样就完全变成了线性表达式可以直接放进LP或MILP里求解。很多新手在这卡住不清楚 (\alpha) 和 (z_s) 是什么角色。记住(\alpha) 是一个自由变量最优时等于VaR(z_s) 是每个场景的辅助损失变量。注意我只列出最基础的单层CVaR形式。实际微网调度里经常把CVaR放进目标函数作为惩罚项即目标 期望收益 - λ × CVaR损失或者目标 期望成本 λ × CVaR。λ就是风险厌恶系数需要做灵敏度分析来选。2.3 目标函数怎么设计风险和收益的权衡在这个模型里目标函数有三种常见设计方式期望收益最大化完全不考虑风险(\lambda0)等价于随机期望规划。CVaR最大化或损失最小化只考虑最坏情况相当于极度保守。均值-CVaR加权最大化 (\mu - \lambda \cdot CVaR)其中 (\mu) 是期望收益(\lambda) 越大越保守。第三种最实用因为可以通过调整 (\lambda) 画出一条“有效前沿”曲线展示利润和风险之间的权衡关系决策者可以按自己的风险偏好选取运行点。我在实际项目中一般会把 (\lambda) 从0开始逐步增大到某个上限每隔0.5取一个点计算得到一系列调度方案再从中选择。值得提醒的是CVaR在目标函数里是以损失形式出现的所以目标函数写法从“最大化收益”变成“最大化收益减CVaR惩罚”需要把收益写成每个场景下的表达式这个我在下面的Matlab实现部分会给出具体写法。另外目标函数中的各项数量级要匹配否则CVaR惩罚项可能过小或过大导致结果失控。我一般会把收益归一化或者用权重折算。3. 动态定价与调度的联动机制3.1 微网内部分布式资源建模微网内的可控资源我通常建模为以下几类每类的数学表达差异很大光伏不可控每个场景下给定出力曲线作为参数输入。注意光伏出力场景之间是相关的不能独立抽样否则会高估光伏的平滑效应。用Copula或直接对总出力做时序抽样更合理。储能可控但耦合时间储能是跨时段耦合的需要引入荷电状态SOC变量 (E_t)满足[ E_{t1} E_t \eta_c P_t^c \Delta t - \frac{P_t^d}{\eta_d} \Delta t ]其中 (P_t^c, P_t^d) 是充放电功率(\eta_c, \eta_d) 是充放电效率。储能建模有几个常见的坑充放电不能同时进行需要二进制变量构成MILPSOC有上下限循环寿命损耗要不要计入目标函数取决于研究精度。如果全部用线性约束储能模型还是一个LP求解速度会快很多但“不能同时充放”必须引入二进制变量所以整个问题基本是MILP。燃气轮机/柴油机可控但非线性发电成本通常是出力的二次函数 (aP^2bPc)。在MILP里我习惯用分段线性化处理分成3-5段即可满足工程精度。还有最小启停时间约束、爬坡约束这些都会增加整数变量数量。3.2 动态定价与需求响应耦合动态定价的核心是价格弹性系数。简单说用户在当前时刻的用电量不仅受当前电价影响还受其他时段电价影响交叉弹性。最常用的是自弹性系数[ q_t q_t^0 \cdot \left(1 \varepsilon_t \cdot \frac{\rho_t - \rho_t^0}{\rho_t^0}\right) ]其中 (q_t^0) 是基准负荷(\rho_t^0) 是基准电价(\rho_t) 是优化出的实时电价(\varepsilon_t) 是自弹性系数通常为负电价越高负荷越低。如果考虑用户跨时段转移负荷还要加交叉弹性矩阵模型会更复杂。这里有个很重要的细节动态定价不是随随便便给电价加个上下限就完事。电价过高会导致用户负荷大幅下降影响微网收入电价过低则无法引导削峰填谷。我一般会设置一个“价格上限不超过基准电价的1.2倍、下限不低于0.8倍”的约束同时要求微网整体收益不低于某个基准值。这样才能保证定价结果既有激励效果又不会让用户流失。3.3 完整约束体系与变量设计一个完整的微网调度模型变量总结如下变量名含义类型(P_t^g, P_t^{pv})向主网购电/售电功率连续(P_t^{mt})燃气轮机出力连续(P_t^c, P_t^d)储能充/放电功率连续(u_t^{c}, u_t^{d})储能充/放状态二进制(E_t)储能SOC连续(P_t^{e})可调负荷削减量连续(\rho_t)零售电价连续(q_t)用户需求响应后负荷连续约束体系包括功率平衡(P_t^{pv} P_t^{mt} P_t^d P_t^g P_t^{e} q_t P_t^c P_t^{s})含网损修正储能约束SOC更新、充放电上下限、不能同时充放机组约束出力上下限、爬坡约束、最小启停时间如需要定价约束电价上下限、负荷响应表达式CVaR约束每个场景下的损失与辅助变量关系这里每个约束在Matlab里都是用YALMIP的约束对象来定义的后面我会给代码片段。注意功率平衡约束里每个场景都要满足这是随机规划里“非预期性”的基本要求有些新手只对期望场景做平衡约束那就不是严格意义上的随机优化了。4. Matlab实现代码架构与核心模块4.1 工具箱选择与环境配置做这套模型我强烈推荐YALMIP Gurobi或CPLEX组合。YALMIP是Matlab下的建模语言你能用接近数学公式的方式直接写约束非常省心。求解器方面Gurobi对MILP的求解速度明显优于MATLAB自带的intlinprog尤其是场景数量上来以后差距会很可观。需要说明的是YALMIP只是建模接口真正的MILP求解靠后端求解器。在安装时把Gurobi的Matlab接口路径添加到Matlab路径中即可。有些国内用户可能不方便直接下载完整版求解器可以用MATLAB自带的intlinprog先跑通小规模模型再迁移到Gurobi。这个迁移在YALMIP框架下基本就是一行代码的事。4.2 场景生成模块怎么做场景生成有几种做法蒙特卡洛直接抽样假设光伏服从Beta分布电价服从正态分布直接抽样生成场景。好处是简单坏处是场景数多求解慢。场景缩减先抽3000-5000个场景再用快速前向选择法缩减到20-50个。这一步能极大降低计算负担效果几乎不损失。拉丁超立方抽样LHS比纯蒙特卡洛效率更高能更好地覆盖尾部区域。我在代码里提供一个简化实现用randn生成相关随机数再加均值扰动。核心代码如下% 场景生成光伏出力 预测值 误差扰动 % 电价场景 预测值 随机扰动 T 24; % 调度时段数 S 500; % 场景数 pv_scen zeros(T, S); price_scen zeros(T, S); mu_pv [0, 0.2, 0.5, 0.8, 1.0, 0.9, ...]; % 各时段预测出力 sigma_pv 0.15 * mu_pv 0.02; % 标准差随出力变化 for s 1:S pv_scen(:, s) mu_pv sigma_pv .* randn(T, 1); price_scen(:, s) mu_price sigma_price .* randn(T, 1); end % 截断处理避免负值 pv_scen(pv_scen 0) 0; price_scen(price_scen 0) min(price_scen(price_scen 0));注意一个细节直接抽样产生的场景往往带有不合理的极端值需要做截断或者拒绝采样。另外如果需要更精细的相关性建模建议用copularnd生成相关性场景而不是独立的randn。4.3 YALMIP建模实现CVaR调度模型的代码框架下面是整个模型的核心代码框架我用YALMIP语法写。为了便于阅读我先给出简化版暂不考虑动态定价的交叉弹性。% 参数设置 T 24; S 50; beta 0.95; lambda 2.0; % 决策变量 P_g sdpvar(T, 1, full); % 与主网交换功率正为购电 P_mt sdpvar(T, 1, full); % 微型燃气轮机出力 P_c sdpvar(T, 1, full); % 储能充电功率 P_d sdpvar(T, 1, full); % 储能放电功率 u_c binvar(T, 1); % 充电状态 u_d binvar(T, 1); % 放电状态 E sdpvar(T1, 1, full); % SOC第1个元素为初始SOC P_cut sdpvar(T, 1, full); % 可调负荷削减量 rho sdpvar(T, 1, full); % 零售电价 % 场景相关变量 P_g_s sdpvar(T, S, full); % 各场景购电功率 P_cut_s sdpvar(T, S, full); % 各场景削减负荷 loss_s sdpvar(1, S, full); % 各场景损失 z_s sdpvar(1, S, full); % CVaR辅助变量 alpha sdpvar(1, 1); % VaR变量 % 约束 Constraints []; % 储能约束 Constraints [Constraints, E(1) 0.2 * E_max]; for t 1:T Constraints [Constraints, ... E(t1) E(t) eta_c * P_c(t) - P_d(t) / eta_d, ... 0 P_c(t) P_c_max * u_c(t), ... 0 P_d(t) P_d_max * u_d(t), ... u_c(t) u_d(t) 1, ... E_min E(t1) E_max]; end Constraints [Constraints, E(T1) 0.2 * E_max]; % 调度周期末SOC恢复 % 功率平衡约束每个场景下都需要满足 for s 1:S Constraints [Constraints, ... P_g_s(:, s) P_pv(:, s) P_mt P_d - P_c ... P_load * (1 - 0.2 * (rho - rho_base) / rho_base) - P_cut_s(:, s)]; Constraints [Constraints, ... 0 P_g_s(:, s), P_g_s(:, s) P_grid_max, ... 0 P_cut_s(:, s), P_cut_s(:, s) 0.1 * P_load]; end % 定价约束 Constraints [Constraints, 0.8 * rho_base rho 1.2 * rho_base]; % 损失与CVaR约束 for s 1:S % 场景s的净收益售电收入 - 购电成本 - 燃料成本 - 削减补偿 revenue_s sum(rho .* (P_load * (1 - 0.2 * (rho - rho_base) / rho_base) - P_cut_s(:, s))) ... sum(c_sell(s, :) .* max(0, -P_g_s(:, s))); cost_s sum(c_buy(s, :) .* max(0, P_g_s(:, s))) ... sum(c_gas * P_mt) sum(c_cut * P_cut_s(:, s)); loss_s(s) - (revenue_s - cost_s) / 1000; % 损失 -净利润 Constraints [Constraints, z_s(s) loss_s(s) - alpha, z_s(s) 0]; end % 目标函数 % 期望收益 - lambda * CVaR Obj - mean(loss_s) - lambda * (alpha 1 / ((1 - beta) * S) * sum(z_s)); % 求解 ops sdpsettings(solver, gurobi, verbose, 1); result optimize(Constraints, -Obj, ops);这段代码比较简略但整体结构是完整的。在实际项目中我还会加上燃气轮机的最小运行时间约束、爬坡约束、以及与主网交换功率的购售电状态变量不能同时买和卖。这些细节代码量不少但换汤不换药关键是理解每个约束的物理含义。4.4 求解与结果输出求解完成后最重要的是结果分析和可视化。我会把结果整理成以下输出24小时调度曲线图光伏出力、负荷、储能充放电、购电功率、燃气轮机出力画在同一张图上直观看到削峰填谷效果。电价与负荷联动图零售电价曲线和用户响应后的负荷曲线放在同一坐标系看定价机制是否有效引导了负荷转移。CVaR有效前沿图横轴是CVaR值风险纵轴是期望收益不同λ点拟合出一条曲线。这里特别提醒MILP求解器返回的double()必须用在所有决策变量上否则你会发现自己后续画图全是NaN。这是新手最常犯的错误在YALMIP里写value(P_g)而不是double(P_g)有时会出问题一定要检查变量类型。5. 典型仿真结果分析与参数调试经验5.1 一个可复现的测试场景我用一个典型微网算例做说明参数如下微网峰值负荷800kW光伏装机600kW储能容量300kWh最大充放电功率100kW燃气轮机额定功率400kW。零售基准电价0.8元/kWh向主网购电分时电价在0.5-1.2元/kWh之间波动。光伏出力在中午时段达到高峰与负荷高峰错位。场景缩减后保留30个代表性场景置信水平β0.95风险厌恶系数λ从0开始以0.5为步长增加到4。5.2 不同风险偏好下的调度结果对比λ期望收益(元)CVaR(元)储能日均充放次数燃气轮机出力占比06832-32501.222%16515-18701.526%26288-11201.831%45805-6402.235%看这个结果能读出不少信息λ增大后期望收益下降但CVaR最坏情况损失显著减小说明系统从“追求高利润”转向“保证下限”。同时储能充放次数和燃气轮机出力占比都在上升因为系统更倾向于用可控电源保证供电可靠性而不是依赖不确定的光伏和波动的市场电价。这个规律在大多数算例下都成立可以作为验证模型正确性的参照。5.3 三段式调试法从易到难排查问题调试这类优化模型我总结了一套三段式方法能少走很多弯路第一段小规模验证。先把T缩到4小时、S缩到5个场景关闭CVaR约束λ0看模型能不能稳定求解。这一阶段主要排查约束写没写错、变量维度是否匹配、求解器是否报infeasible。第二段逐步加复杂度。先加CVaR约束λ从0.5开始再加动态定价模型最后加场景缩减。每加一个模块都跑一次确认结果合理后再往下走。这样出了问题能立刻定位是哪个模块的锅。第三段灵敏度分析。固定其他参数扫描λ、β、储能容量这几个关键参数看结果变化趋势是否符合物理直觉。如果λ增大反而CVaR增大那基本可以断定代码有bug。6. 常见问题与避坑指南6.1 求解器报infeasible怎么办这是最让人头疼的问题。我的排查顺序是先检查功率平衡约束负荷充电是否可能大于所有电源的最大出力之和再检查储能SOC约束初始SOC 20%最终要求20%但容量太小、充放电效率太低时可能达不到。检查购电功率上限电网容量给得太小极端场景下无法满足负荷。检查需求响应后的负荷是否为负电价降低导致负荷增加但可能超过系统承载上限。建议把不可行约束排查工具打开Gurobi有IIS不可行子集计算功能YALMIP里用check(Constraints)也能定位不可行约束。我发现大多数infeasible问题都出在储能SOC的“最终值约束”上因为充放电效率和容量约束叠加后可能确实无法在24小时内把SOC恢复到初始水平。这时候要么放宽终值约束要么增大储能容量。6.2 场景数量与求解时间的平衡场景数量是本模型复杂度的重要影响因素。S30时Gurobi在几十秒内求解S100时可能要几分钟到十几分钟S1000时基本没法实际使用。我的经验是用场景缩减把原始5000个场景缩到30-50个既能保留分布特征又不会让求解时间爆炸。缩得太少比如5个会导致CVaR估计偏差非常大缩得太多又让模型没法实用30个左右是一个经验上的甜点位。6.3 定价结果震荡与储能频繁充放有时候优化出来的电价曲线会出现“锯齿状”波动相邻时段电价忽高忽低。这不是正常的动态定价结果而是因为目标函数里缺少对电价平顺性的惩罚或者价格弹性设置过小导致用户响应不足。解决办法是加电价变化幅度的罚项[ -\theta \cdot \sum_{t2}^{T} (\rho_t - \rho_{t-1})^2 ]储能频繁充放也有类似原因——目标函数没有对充放次数或SOC变化率做惩罚。可以加SOC变化惩罚项或者限制一天内充放电切换次数这在MILP里可以用额外的二进制变量实现。6.4 常见问题速查表现象可能原因排查方向求解不可行储能终值SOC无法恢复放宽E(T1)约束求解结果光伏弃光率异常光伏场景生成边界处理不当检查pv_scen截断逻辑峰谷电价效果不明显价格弹性系数太小增大ε绝对值极端场景下损失巨大λ太小或β太小增大λ或β求解时间过长场景数过多或二进制变量过多场景缩减、减少分段数7. 从代码到论文进一步提升研究方向这套模型框架跑通之后可以往以下几个方向延伸一是多微网协同。单个微网的CVaR调度只考虑自身利益但多个微网之间可以通过能量共享降低整体风险。这一般写成两阶段鲁棒优化或者纳什议价模型求解复杂度上了一个台阶。二是考虑电碳耦合。当前政策环境下碳排放约束不能忽视可以把碳配额和碳交易成本纳入目标函数观察碳价对调度策略和动态定价的影响。三是用强化学习替代场景优化。场景优化的瓶颈在于需要事先知道不确定参数的概率分布。用深度强化学习可以做到无模型决策直接在训练中学习风险敏感策略这是目前很火的方向。但我要强调一点不要为了创新而创新。CVaR微网动态定价这个方向已经比较成熟想要出成果关键在“结合点”上下功夫。比如你是做通信的可以研究5G基站参与微网调度的风险建模你是做交通的可以考虑电动汽车充放电的时空不确定性对微网CVaR的影响。把自己的领域优势和微网调度结合起来比单纯卷CVaR建模算法更容易出成果。8. 我对这套方案的个人体会做了几个月的CVaR微网调度项目我最大的感受是模型复杂程度和实际工程价值之间的平衡比想象中难把握。理论文章喜欢把模型写得越精细越好但实际运行中很多参数比如价格弹性根本测量不到那么精确即便算出一个“最优电价”执行时也可能因为用户行为偏差而达不到预期。所以我现在做项目优先保证模型可解释性再做精细化。代码架构上我强烈建议从一开始就按模块拆分文件generate_scenarios.m、build_model.m、solve_model.m、plot_results.m。不要在一个脚本里从头写到尾否则每次调参数都痛苦不堪。场景缩减、CVaR线性化、需求响应模型都做成独立函数后面换模型或换数据都方便。另外一个小技巧在求解MILP时设置一个合理的MIPGap比如1%或0.5%能显著缩短求解时间而且结果误差完全可以接受。Gurobi里用ops.gurobi.MIPGap 0.01CPLEX里用ops.cplex.mip.tolerances.mipgap 0.01。很多论文里的“最优解”其实也就是在一定gap下的近似最优解这不影响结论的可靠性。如果你准备自己也动手搭这套模型我的建议是先不要急着上动态定价和需求响应先把“光伏储能燃气轮机CVaR约束”这个子集跑通再加价格联动模型。这个顺序能帮你少踩至少一半的坑。这里说的每一条都是我自己踩过之后才总结出来的照着走可以让你的调试过程顺很多。