做电力市场优化的人大概都有这种经历预测模型把明天的结算电价估得头头是道真到日前出清时一个高温天气、一条线路检修价格就能走出三倍方差。我最初做的自调度程序也是这样把点预测当成真值喂给LP目标函数和约束都对但结果到了结算日总在利润表的边缘来回试探。问题不在于优化器而在于我把“价格不确定”压缩成了“价格等于期望值”。后来我把模型换成一种结合CVaR的基于矩的分布鲁棒优化DRO在IEEE6、IEEE30和IEEE118节点算例上重跑了一遍才算真正把尾部风险装进了调度决策。这篇文章就记录这个模型的建模思路、MATLAB实现和三种规模算例下的实测现象适合正在做自调度、随机优化和风险决策的电力方向研究生和工程师参考。1. 自调度问题比想象中更需要处理价格尾部风险1.1 自调度的目标函数和传统约束自调度英文叫self-scheduling是发电商、储能系统或聚合商在日前市场上自主决定机组启停和出力计划的过程。它和传统集中调度的最大区别在于集中调度是系统调度机构统一优化全部机组而自调度完全站在市场主体自己的角度目标很直接要么在满足自身技术约束的前提下让总收益最大要么让总成本最小。收益或者成本里最大的不可控项就是市场结算边际电价。如果把电价当成已知参数这个优化问题其实不复杂。典型模型可以写成min_{q,u} ∑_{i,t}(c_i q_{i,t} S_i u_{i,t}) ∑_{t1}^{T} π_t q_buy,ts.t. q_{i,t}^{min} u_{i,t} ≤ q_{i,t} ≤ q_{i,t}^{max} u_{i,t}爬坡约束|q_{i,t} - q_{i,t-1}| ≤ R_i电量平衡∑_i q_{i,t} q_buy,t D_t最小启停时间约束、机组状态耦合约束等。这里q是机组出力u是启停0-1变量π_t是时段电价q_buy是从市场买入的电量。电价π_t如果已知这是一个典型的混合整数线性规划问题MATLAB里用intlinprog都能跑规模再大一点用Gurobi也很快。但现实世界里π_t要等到日前市场出清之后才知道它是一个随机向量不是你随意指定的一个数。于是问题就从“确定性MILP”变成了“随机优化问题”。很多初做自调度的人包括我自己第一版代码都默认电价等于预测值。预测模型给一个期望电价就直接丢进优化器。这样做在普通日子里问题不大因为电价在均值附近波动出清结果和优化结果相差有限。可一旦出现极端天气、网络阻塞或者某台大机组临时检修电价会直接冲到尾部分位数。此时前期基于点预测做的购电计划、机组启停方案全部失真轻则利润缩水重则实时平衡市场倒贴罚款。这正是自调度问题区别于其他确定性调度问题的地方损失函数不是对称的。电价暴涨时如果你是购电商一次尾部事件造成的损失可能覆盖几十个正常交易日的盈利如果你是发电商电价暴跌时收益同样可能一夜归零。自调度的风险本质上是价格分布尾部的风险而尾巴恰恰是最难预测的部分。1.2 期望值最优与CVaR最优的差距常见的随机自调度模型会用“期望成本最小”或者“情景平均成本最小”作为目标函数。构造几百上千个电价场景然后求平均成本最低的调度方案。这个方法比点预测前进了一大步因为它在一定程度上包含了电价波动信息。但它有一个根本问题期望值把好结果和坏结果平均掉了。设想你明天有95%的概率电价稳定在50美元/兆瓦时5%的概率电价暴涨到300美元/兆瓦时。期望电价是62.5美元/兆瓦时一个以期望成本为目标的模型会觉得“电价还好嘛”于是继续按着原来的方式大量购电。然而那5%的尾部事件一旦发生单次结算就足以让你把前面几十天的利润全部吐回去。要刻画这种尾部风险最有用的指标就是CVaR条件风险价值。CVaR_α定义的是损失分布中超过α分位点的那部分损失的期望值。以α0.95为例CVaR_95是所有结果里最坏那5%情况下的平均损失。它不是把尾部平均到全体情景里而是只看尾部本身。这个指标的好处是它天然只关心极端风险并且可以通过一个简单的辅助变量η线性化到优化模型中。在自调度问题里加入CVaR项模型会自动在电价极易出现极端上涨的时段减少购电或者在电价极易暴跌的时段减少发电而不是无动于衷地站在原地。举一个我实际算过的简单例子。某机组出力上限100MW发电成本30美元/兆瓦时日前电价有90%概率落在40美元/兆瓦时附近10%概率跌到15美元/兆瓦时。用期望电价优化机组会保持满发因为平均价格确实高于成本。可一旦遇到那10%的低电价事件每发一度电都在亏钱。加入CVaR_95之后优化器会主动把部分时段出力压低甚至停机因为停机成本比硬着头皮满发的尾部损失要小。这件事靠期望目标模型是做不到的哪怕你把场景数从500加到5000也做不到因为平均下来尾部损失总会被稀释。1.3 CVaR与DRO为什么放在一起CVaR虽然好但它要求我们知道电价的完整概率分布至少要有足够可靠的情景数据。现实情况下我们往往只有历史电价的时间序列或者来自点预测模型的一阶矩和二阶矩信息。真实分布到底是正态、t分布、还是带有明显偏峰厚尾的混合分布我们很难准确判断。这时候硬套一个正态分布做CVaR估计出来的尾部风险本身就有很大误差。DRO就是专门对付这种“分布不确定”的框架。它不试图估计出一个精确分布而是构造一个分布模糊集把所有与已知统计量保持一致的分布都放进这个集合里然后求最坏情况下的目标函数值。CVaR解决的是“风险怎么度量”DRO解决的是“分布不知道怎么办”两者放在一起是天然互补的关系。在本文模型中我采用的DRO目标是min_{q,u} max_{P∈D} E_P[L(q,u,π)] λ · max_{P∈D} CVaR_α,P[L(q,u,π)]其中L是总成本D是基于一阶矩和二阶矩信息构造的模糊集。这里max算子表示我们宁可面对一组最坏的合法电价分布也要保证调度方案还能接受。λ是风险偏好系数λ越大模型越保守越愿意牺牲期望成本来换取尾部风险降低。2. 一个能落地的基于矩DROCVaR模型2.1 有限支撑上的矩模糊集很多同学一听到基于矩的DRO、模糊集、概率测度空间这些词就觉得头大。其实落到数值实现上有一个非常实用的做法先用历史数据或预测误差生成一批价格场景π_s然后用概率权重p_s来描述所有可能的离散分布。在这个有限支撑集上矩模糊集可以写成D { p ≥ 0, ∑_s p_s 1 :(∑_s p_s π_s - μ)^T Σ^{-1} (∑_s p_s π_s - μ) ≤ ε_1,∑_s p_s (π_s - μ)(π_s - μ)^T ⪯ ε_2 Σ }这里μ和Σ是参考均值向量和参考协方差矩阵通常由历史数据或预测模型给出。ε_1和ε_2控制模糊集的大小。第一个约束限制权重分布的一阶矩不能偏离参考均值太远第二个约束限制二阶中心矩不能超过参考协方差的某个倍数。用矩阵不等式的写法看起来复杂但本质上就是给概率权重加了一组二次约束和半定约束。当ε_1 ε_2 0时模型退化为均值-协方差完全锁定的分布集当ε值逐渐增大模型对分布误设的容忍度也增大调度结果自然更保守。我选择基于矩而不是基于Wasserstein距离的模糊集原因很实际一阶矩和二阶矩在电力市场中容易获得调度部门手上通常有点预测结果和预测误差协方差矩阵但很难给出一个精确的Wasserstein半径。而且矩模糊集和CVaR本身的矩特性很搭风险指标的波动范围能被比较自然地控制住。2.2 CVaR目标在有限支撑下的计算Rockafellar和Uryasev在2000年给出了一个经典的CVaR等价表达式CVaR_α(L) min_η { η 1/(1-α) E_P[ max(0, L - η) ] }在有限支撑集上E_P[max(0, L - η)]展开为∑_s p_s max(0, L_s - η)。由于p_s是模糊集里的未知权重目标函数里出现了p_s与max项相乘的双线性结构。这个问题不能直接丢进通用优化器需要先做变量拆分。我的做法是分两层求解。外层变量是调度计划q和启停状态u内层在给定计划的情况下求解最坏概率权重p以及CVaR辅助变量η和z_s。内层是一个线性规划max_{p, η, z} ∑_s p_s L_s λ { η (1/(1-α)) ∑_s p_s z_s }s.t. z_s ≥ L_s - η, z_s ≥ 0,(∑_s p_s π_s - μ)^T Σ^{-1} (∑_s p_s π_s - μ) ≤ ε_1,∑_s p_s (π_s - μ)(π_s - μ)^T ⪯ ε_2 Σ,∑_s p_s 1, p ≥ 0.给定外层调度计划后L_s是已知常数所有变量都是线性的或者只有二次/半定约束这样的内层问题用Gurobi或Mosek都能稳定求解。外层再用MATLAB的fmincon或者遗传算法等无梯度优化器搜索调度计划。这个内外层结构听起来没有直接写一个大型MILP那么优雅但实际跑起来非常稳。它能避开p_s和z_s乘积带来的大规模双线性项让我在IEEE6和IEEE30上都能快速迭代。2.3 完整模型的紧凑写法把机组约束和模糊集目标拼在一起完整模型可以用如下紧凑形式描述min_{q,u} F(q,u)其中F(q,u) max_{p∈D} { E_p[L(q,u,π)] λ · CVaR_α,p[L(q,u,π)] }s.t. 机组出力上下限 q_{i,t}^{min} u_{i,t} ≤ q_{i,t} ≤ q_{i,t}^{max} u_{i,t}爬坡约束 |q_{i,t} - q_{i,t-1}| ≤ R_i最小启停时间约束根据具体机组类型设置on/off持续时间电能量平衡 ∑_i q_{i,t} q_buy,t D_tq_buy,t ≥ 0总购电量与售电收益之间的结算由L(q,u,π)线性表示。注意T个时段的价格向量π不再是一个固定参数而是在模糊集D中变化。由于目标是最大化最坏情况成本模型会自动把调度计划往“在所有合理分布下都不差”的方向推。这个模型对电力市场的几个核心特征做了简化未显式建模网络潮流只用了节点总平衡这是自调度问题的常见处理方式未考虑备用市场收益只关注电能量市场开关机成本和燃料成本采用线性近似避免非线性项干扰DRO求解。在实际项目中这些简化可以根据需要逐步加入。加入网络约束时可以在矩模糊集之外再叠加一个潮流约束用DCOPF或者交流潮流进行安全校验加入备用收益时只需要在目标函数中增加一个与价格弱相关的二次项即可。3. MATLAB实现和三种IEEE节点数据接入3.1 用MATPOWER准备机组参数做IEEE节点系统仿真我强烈建议直接用MATPOWER工具箱不要自己手工录入机组数据。MATPOWER自带case6ww、case30、case118等经典测试系统只需要三行代码就能读出所有参数mpc loadcase(case6ww); % 换成 case30 或 case118 即可切换系统 gen mpc.gen; bus mpc.bus; branch mpc.branch;MATPOWER的gen矩阵里包含发电机所在母线编号、有功出力上限、下限、爬坡率、启停成本、燃料成本系数等字段。bus矩阵里包含负荷有功、无功和母线电压基准。case6ww有3台机组适合快速调试case30有6台机组是中等规模测试的最佳选择case118有54台机组已经逼近大规模实际系统的计算复杂度。我拿到这些数据后第一步不是直接建优化模型而是先画一张机组参数分布图看看系统的边际成本大概在什么范围发电机组的启停成本各是多少。因为电价场景的均值、方差要跟系统边际成本大致匹配否则DRO优化结果会明显失真。比如case6ww里机组容量很小如果我把电价均值设成120美元/兆瓦时那所有机组都会满发鲁棒性优化就没有意义了。3.2 电价场景生成与矩参数设置电价场景质量直接决定DRO模型效果。我采用两条腿走路先用多元正态分布作为基准生成器这一步只是为了拿到一个具有合理时空相关性的初始场景集随后再用矩匹配方法对场景进行修正让场景集的一阶矩和二阶矩严格匹配目标μ和Σ。一个简单的生成代码如下T 24; mu 50 10*sin(2*pi*(0:T-1)/24); Sigma 15^2 * exp(-abs((0:T-1)-(0:T-1))/3); rng(2024); Pi mvnrnd(mu, Sigma, S); Pi max(Pi, 5); % 电价不能是负的简单截断这组参数的含义是日均电价在40-60美元/兆瓦时之间波动相邻时段电价相关性随间隔衰减时间常数3小时。实际项目中μ和Σ应该来自预测系统的误差统计或者从历史日度电价序列估计。截断到5美元/兆瓦时是为了避免生成不合理负电价虽然部分市场确实会出现负电价但那会让DRO模型的尾部更加极端首次调试不建议混入这种特殊工况。生成场景后我会用样本矩修正工具把场景均值拉回μ把协方差修正到Σ附近。这一步很重要因为直接从mvnrnd生成的有限样本会有抽样误差内层模糊集约束里用的μ和Σ如果不匹配内层LP在ε较小时甚至会出现无解。3.3 求解器与YALMIP核心代码在我这个模型中内层LP是每次外层迭代都要重复求解的所以求解速度非常重要。我用YALMIP做建模层后端接Gurobi求解LP和MILP。YALMIP的好处是二次约束和半定约束可以直接写成约束表达式不需要自己手动推导对偶形式。内层最坏分布求解的核心代码结构如下p sdpvar(S, 1, full); eta sdpvar(1, 1); z sdpvar(S, 1, full); constraints [sum(p) 1, p 0, z 0]; for s 1:S constraints [constraints, z(s) Ls(s) - eta]; end % 一阶矩约束 dev Pi * p - mu; constraints [constraints, dev * inv(Sigma) * dev eps1]; % 二阶矩约束半定形式 M zeros(T, T); for s 1:S d Pi(s,:) - mu; M M p(s) * (d * d); end constraints [constraints, M eps2 * Sigma]; objective sum(p .* Ls) lam * (eta sum(p .* z) / (1 - alpha)); ops sdpsettings(solver, gurobi, verbose, 2); optimize(constraints, -objective, ops);外层调度计划我用YALMIP的sdpvar定义机组出力、启停状态和买入电量然后把内层最优值作为目标函数的一部分返回给fmincon。整体时间在一个IEEE6算例上大约能控制在几十秒到几分钟取决于场景数和外层迭代次数。如果对YALMIP不熟也可以直接用MATLAB的quadprog和mosek接口写但那样代码量会翻倍。YALMIP在处理半定约束时特别省事尤其是二阶矩矩阵不等式的写法手写要维护一大堆矩阵索引YALMIP直接声明矩阵变量就行。3.4 从IEEE6到IEEE118的数据扩展IEEE6节点系统只有3台机组、24个时段外层变量约72个连续变量和72个整数变量代码跑起来非常快适合验证模型逻辑。IEEE30节点系统有6台机组变量数翻倍但依然在Gurobi的可控范围内。真正让我头疼的是IEEE11854台机组、24个时段如果全部建模机组组合变量规模超过2600个加上爬坡和最小启停时间约束外层迭代一次就要好几分钟。我的处理方式是分三步先用一个简化机组组合模型跑一次固定少数关键机组的启停状态再把固定后的启停状态带入完整DRO模型只优化连续出力最后进行一轮启停邻域搜索看是否有必要改变个别机组的启停状态。这样牺牲了一部分全局最优性但能获得一个工程上可接受的近似最优解。在电力市场实际应用中市场出清时间窗口有限这种“先粗后精”的两阶段策略远比硬啃大规模MILP实用。4. 三种节点系统下的调度结果和风险指标4.1 IEEE6基准模型的保守性肉眼可见在IEEE6节点系统上我以α0.95、λ0.5、ε10.1、ε21.2为例做了对比实验。表格如下模型期望成本最坏情景成本CVaR_95CPU时间确定性点预测102.8k128.4k119.1k0.8s随机期望值103.2k126.9k116.5k15.3sDROCVaR104.1k119.8k108.6k47.6s注意这里成本单位是“千美元”只是针对我构造的模拟电价场景不是真实市场数据。关键结论是DROCVaR让期望成本从102.8k上升到104.1k只增加了1.3%但CVaR_95从119.1k降到了108.6k下降了8.8%。最坏情景成本也从128.4k降到119.8k。这意味着模型用很小的期望成本牺牲换来了尾部风险的大幅压缩。在调度曲线上这个变化非常直观。确定性模型在电价最低时段依然安排了机组运行因为它没有预见到极端低价会长期持续DRO模型则在尾部风险高的时段直接压低甚至关停机组即使这些时段在期望情景里看起来还有利润。4.2 IEEE30风险系数λ怎么改变调度策略IEEE30节点系统比IEEE6复杂不少机组数量和负荷曲线都更接近实际。我固定α0.95、ε10.1、ε21.2把λ从0逐渐增加到5观察调度策略变化。λ取值期望成本CVaR_95高价时段购电占比0208.6k242.3k65%0.2209.8k235.1k59%0.5211.2k222.7k51%1.0214.9k213.4k44%2.0219.7k207.6k38%5.0228.3k201.2k33%可以看到λ从0增加到5时期望成本从208.6k上升到228.3k增幅接近9.4%但CVaR_95从242.3k下降到201.2k降幅达到17%。高价时段购电比例从65%降到33%说明模型在风险厌恶程度提升时会主动放弃价格尖峰时段的购电机会把电量转移到本地机组或者更便宜的时段。这条曲线对于实际操作特别有意义。调度员不需要机械地选一个λ而是可以通过λ-风险曲线在期望成本和尾部风险之间找平衡点。如果企业风险承受能力低就选曲线右端如果更看重平均利润就选左端。4.3 IEEE118计算规模与折中IEEE118是我这次实验里最接近实际规模的一个系统。54台机组、24时段、500个价格场景直接用完整MILP跑会非常吃力。我采用固定启停连续出力两阶段策略后单次外层迭代从无法收敛降到几百秒内。从结果看DROCVaR模型在IEEE
电力市场自调度风险优化：基于矩DRO与CVaR模型及MATLAB实现
做电力市场优化的人大概都有这种经历预测模型把明天的结算电价估得头头是道真到日前出清时一个高温天气、一条线路检修价格就能走出三倍方差。我最初做的自调度程序也是这样把点预测当成真值喂给LP目标函数和约束都对但结果到了结算日总在利润表的边缘来回试探。问题不在于优化器而在于我把“价格不确定”压缩成了“价格等于期望值”。后来我把模型换成一种结合CVaR的基于矩的分布鲁棒优化DRO在IEEE6、IEEE30和IEEE118节点算例上重跑了一遍才算真正把尾部风险装进了调度决策。这篇文章就记录这个模型的建模思路、MATLAB实现和三种规模算例下的实测现象适合正在做自调度、随机优化和风险决策的电力方向研究生和工程师参考。1. 自调度问题比想象中更需要处理价格尾部风险1.1 自调度的目标函数和传统约束自调度英文叫self-scheduling是发电商、储能系统或聚合商在日前市场上自主决定机组启停和出力计划的过程。它和传统集中调度的最大区别在于集中调度是系统调度机构统一优化全部机组而自调度完全站在市场主体自己的角度目标很直接要么在满足自身技术约束的前提下让总收益最大要么让总成本最小。收益或者成本里最大的不可控项就是市场结算边际电价。如果把电价当成已知参数这个优化问题其实不复杂。典型模型可以写成min_{q,u} ∑_{i,t}(c_i q_{i,t} S_i u_{i,t}) ∑_{t1}^{T} π_t q_buy,ts.t. q_{i,t}^{min} u_{i,t} ≤ q_{i,t} ≤ q_{i,t}^{max} u_{i,t}爬坡约束|q_{i,t} - q_{i,t-1}| ≤ R_i电量平衡∑_i q_{i,t} q_buy,t D_t最小启停时间约束、机组状态耦合约束等。这里q是机组出力u是启停0-1变量π_t是时段电价q_buy是从市场买入的电量。电价π_t如果已知这是一个典型的混合整数线性规划问题MATLAB里用intlinprog都能跑规模再大一点用Gurobi也很快。但现实世界里π_t要等到日前市场出清之后才知道它是一个随机向量不是你随意指定的一个数。于是问题就从“确定性MILP”变成了“随机优化问题”。很多初做自调度的人包括我自己第一版代码都默认电价等于预测值。预测模型给一个期望电价就直接丢进优化器。这样做在普通日子里问题不大因为电价在均值附近波动出清结果和优化结果相差有限。可一旦出现极端天气、网络阻塞或者某台大机组临时检修电价会直接冲到尾部分位数。此时前期基于点预测做的购电计划、机组启停方案全部失真轻则利润缩水重则实时平衡市场倒贴罚款。这正是自调度问题区别于其他确定性调度问题的地方损失函数不是对称的。电价暴涨时如果你是购电商一次尾部事件造成的损失可能覆盖几十个正常交易日的盈利如果你是发电商电价暴跌时收益同样可能一夜归零。自调度的风险本质上是价格分布尾部的风险而尾巴恰恰是最难预测的部分。1.2 期望值最优与CVaR最优的差距常见的随机自调度模型会用“期望成本最小”或者“情景平均成本最小”作为目标函数。构造几百上千个电价场景然后求平均成本最低的调度方案。这个方法比点预测前进了一大步因为它在一定程度上包含了电价波动信息。但它有一个根本问题期望值把好结果和坏结果平均掉了。设想你明天有95%的概率电价稳定在50美元/兆瓦时5%的概率电价暴涨到300美元/兆瓦时。期望电价是62.5美元/兆瓦时一个以期望成本为目标的模型会觉得“电价还好嘛”于是继续按着原来的方式大量购电。然而那5%的尾部事件一旦发生单次结算就足以让你把前面几十天的利润全部吐回去。要刻画这种尾部风险最有用的指标就是CVaR条件风险价值。CVaR_α定义的是损失分布中超过α分位点的那部分损失的期望值。以α0.95为例CVaR_95是所有结果里最坏那5%情况下的平均损失。它不是把尾部平均到全体情景里而是只看尾部本身。这个指标的好处是它天然只关心极端风险并且可以通过一个简单的辅助变量η线性化到优化模型中。在自调度问题里加入CVaR项模型会自动在电价极易出现极端上涨的时段减少购电或者在电价极易暴跌的时段减少发电而不是无动于衷地站在原地。举一个我实际算过的简单例子。某机组出力上限100MW发电成本30美元/兆瓦时日前电价有90%概率落在40美元/兆瓦时附近10%概率跌到15美元/兆瓦时。用期望电价优化机组会保持满发因为平均价格确实高于成本。可一旦遇到那10%的低电价事件每发一度电都在亏钱。加入CVaR_95之后优化器会主动把部分时段出力压低甚至停机因为停机成本比硬着头皮满发的尾部损失要小。这件事靠期望目标模型是做不到的哪怕你把场景数从500加到5000也做不到因为平均下来尾部损失总会被稀释。1.3 CVaR与DRO为什么放在一起CVaR虽然好但它要求我们知道电价的完整概率分布至少要有足够可靠的情景数据。现实情况下我们往往只有历史电价的时间序列或者来自点预测模型的一阶矩和二阶矩信息。真实分布到底是正态、t分布、还是带有明显偏峰厚尾的混合分布我们很难准确判断。这时候硬套一个正态分布做CVaR估计出来的尾部风险本身就有很大误差。DRO就是专门对付这种“分布不确定”的框架。它不试图估计出一个精确分布而是构造一个分布模糊集把所有与已知统计量保持一致的分布都放进这个集合里然后求最坏情况下的目标函数值。CVaR解决的是“风险怎么度量”DRO解决的是“分布不知道怎么办”两者放在一起是天然互补的关系。在本文模型中我采用的DRO目标是min_{q,u} max_{P∈D} E_P[L(q,u,π)] λ · max_{P∈D} CVaR_α,P[L(q,u,π)]其中L是总成本D是基于一阶矩和二阶矩信息构造的模糊集。这里max算子表示我们宁可面对一组最坏的合法电价分布也要保证调度方案还能接受。λ是风险偏好系数λ越大模型越保守越愿意牺牲期望成本来换取尾部风险降低。2. 一个能落地的基于矩DROCVaR模型2.1 有限支撑上的矩模糊集很多同学一听到基于矩的DRO、模糊集、概率测度空间这些词就觉得头大。其实落到数值实现上有一个非常实用的做法先用历史数据或预测误差生成一批价格场景π_s然后用概率权重p_s来描述所有可能的离散分布。在这个有限支撑集上矩模糊集可以写成D { p ≥ 0, ∑_s p_s 1 :(∑_s p_s π_s - μ)^T Σ^{-1} (∑_s p_s π_s - μ) ≤ ε_1,∑_s p_s (π_s - μ)(π_s - μ)^T ⪯ ε_2 Σ }这里μ和Σ是参考均值向量和参考协方差矩阵通常由历史数据或预测模型给出。ε_1和ε_2控制模糊集的大小。第一个约束限制权重分布的一阶矩不能偏离参考均值太远第二个约束限制二阶中心矩不能超过参考协方差的某个倍数。用矩阵不等式的写法看起来复杂但本质上就是给概率权重加了一组二次约束和半定约束。当ε_1 ε_2 0时模型退化为均值-协方差完全锁定的分布集当ε值逐渐增大模型对分布误设的容忍度也增大调度结果自然更保守。我选择基于矩而不是基于Wasserstein距离的模糊集原因很实际一阶矩和二阶矩在电力市场中容易获得调度部门手上通常有点预测结果和预测误差协方差矩阵但很难给出一个精确的Wasserstein半径。而且矩模糊集和CVaR本身的矩特性很搭风险指标的波动范围能被比较自然地控制住。2.2 CVaR目标在有限支撑下的计算Rockafellar和Uryasev在2000年给出了一个经典的CVaR等价表达式CVaR_α(L) min_η { η 1/(1-α) E_P[ max(0, L - η) ] }在有限支撑集上E_P[max(0, L - η)]展开为∑_s p_s max(0, L_s - η)。由于p_s是模糊集里的未知权重目标函数里出现了p_s与max项相乘的双线性结构。这个问题不能直接丢进通用优化器需要先做变量拆分。我的做法是分两层求解。外层变量是调度计划q和启停状态u内层在给定计划的情况下求解最坏概率权重p以及CVaR辅助变量η和z_s。内层是一个线性规划max_{p, η, z} ∑_s p_s L_s λ { η (1/(1-α)) ∑_s p_s z_s }s.t. z_s ≥ L_s - η, z_s ≥ 0,(∑_s p_s π_s - μ)^T Σ^{-1} (∑_s p_s π_s - μ) ≤ ε_1,∑_s p_s (π_s - μ)(π_s - μ)^T ⪯ ε_2 Σ,∑_s p_s 1, p ≥ 0.给定外层调度计划后L_s是已知常数所有变量都是线性的或者只有二次/半定约束这样的内层问题用Gurobi或Mosek都能稳定求解。外层再用MATLAB的fmincon或者遗传算法等无梯度优化器搜索调度计划。这个内外层结构听起来没有直接写一个大型MILP那么优雅但实际跑起来非常稳。它能避开p_s和z_s乘积带来的大规模双线性项让我在IEEE6和IEEE30上都能快速迭代。2.3 完整模型的紧凑写法把机组约束和模糊集目标拼在一起完整模型可以用如下紧凑形式描述min_{q,u} F(q,u)其中F(q,u) max_{p∈D} { E_p[L(q,u,π)] λ · CVaR_α,p[L(q,u,π)] }s.t. 机组出力上下限 q_{i,t}^{min} u_{i,t} ≤ q_{i,t} ≤ q_{i,t}^{max} u_{i,t}爬坡约束 |q_{i,t} - q_{i,t-1}| ≤ R_i最小启停时间约束根据具体机组类型设置on/off持续时间电能量平衡 ∑_i q_{i,t} q_buy,t D_tq_buy,t ≥ 0总购电量与售电收益之间的结算由L(q,u,π)线性表示。注意T个时段的价格向量π不再是一个固定参数而是在模糊集D中变化。由于目标是最大化最坏情况成本模型会自动把调度计划往“在所有合理分布下都不差”的方向推。这个模型对电力市场的几个核心特征做了简化未显式建模网络潮流只用了节点总平衡这是自调度问题的常见处理方式未考虑备用市场收益只关注电能量市场开关机成本和燃料成本采用线性近似避免非线性项干扰DRO求解。在实际项目中这些简化可以根据需要逐步加入。加入网络约束时可以在矩模糊集之外再叠加一个潮流约束用DCOPF或者交流潮流进行安全校验加入备用收益时只需要在目标函数中增加一个与价格弱相关的二次项即可。3. MATLAB实现和三种IEEE节点数据接入3.1 用MATPOWER准备机组参数做IEEE节点系统仿真我强烈建议直接用MATPOWER工具箱不要自己手工录入机组数据。MATPOWER自带case6ww、case30、case118等经典测试系统只需要三行代码就能读出所有参数mpc loadcase(case6ww); % 换成 case30 或 case118 即可切换系统 gen mpc.gen; bus mpc.bus; branch mpc.branch;MATPOWER的gen矩阵里包含发电机所在母线编号、有功出力上限、下限、爬坡率、启停成本、燃料成本系数等字段。bus矩阵里包含负荷有功、无功和母线电压基准。case6ww有3台机组适合快速调试case30有6台机组是中等规模测试的最佳选择case118有54台机组已经逼近大规模实际系统的计算复杂度。我拿到这些数据后第一步不是直接建优化模型而是先画一张机组参数分布图看看系统的边际成本大概在什么范围发电机组的启停成本各是多少。因为电价场景的均值、方差要跟系统边际成本大致匹配否则DRO优化结果会明显失真。比如case6ww里机组容量很小如果我把电价均值设成120美元/兆瓦时那所有机组都会满发鲁棒性优化就没有意义了。3.2 电价场景生成与矩参数设置电价场景质量直接决定DRO模型效果。我采用两条腿走路先用多元正态分布作为基准生成器这一步只是为了拿到一个具有合理时空相关性的初始场景集随后再用矩匹配方法对场景进行修正让场景集的一阶矩和二阶矩严格匹配目标μ和Σ。一个简单的生成代码如下T 24; mu 50 10*sin(2*pi*(0:T-1)/24); Sigma 15^2 * exp(-abs((0:T-1)-(0:T-1))/3); rng(2024); Pi mvnrnd(mu, Sigma, S); Pi max(Pi, 5); % 电价不能是负的简单截断这组参数的含义是日均电价在40-60美元/兆瓦时之间波动相邻时段电价相关性随间隔衰减时间常数3小时。实际项目中μ和Σ应该来自预测系统的误差统计或者从历史日度电价序列估计。截断到5美元/兆瓦时是为了避免生成不合理负电价虽然部分市场确实会出现负电价但那会让DRO模型的尾部更加极端首次调试不建议混入这种特殊工况。生成场景后我会用样本矩修正工具把场景均值拉回μ把协方差修正到Σ附近。这一步很重要因为直接从mvnrnd生成的有限样本会有抽样误差内层模糊集约束里用的μ和Σ如果不匹配内层LP在ε较小时甚至会出现无解。3.3 求解器与YALMIP核心代码在我这个模型中内层LP是每次外层迭代都要重复求解的所以求解速度非常重要。我用YALMIP做建模层后端接Gurobi求解LP和MILP。YALMIP的好处是二次约束和半定约束可以直接写成约束表达式不需要自己手动推导对偶形式。内层最坏分布求解的核心代码结构如下p sdpvar(S, 1, full); eta sdpvar(1, 1); z sdpvar(S, 1, full); constraints [sum(p) 1, p 0, z 0]; for s 1:S constraints [constraints, z(s) Ls(s) - eta]; end % 一阶矩约束 dev Pi * p - mu; constraints [constraints, dev * inv(Sigma) * dev eps1]; % 二阶矩约束半定形式 M zeros(T, T); for s 1:S d Pi(s,:) - mu; M M p(s) * (d * d); end constraints [constraints, M eps2 * Sigma]; objective sum(p .* Ls) lam * (eta sum(p .* z) / (1 - alpha)); ops sdpsettings(solver, gurobi, verbose, 2); optimize(constraints, -objective, ops);外层调度计划我用YALMIP的sdpvar定义机组出力、启停状态和买入电量然后把内层最优值作为目标函数的一部分返回给fmincon。整体时间在一个IEEE6算例上大约能控制在几十秒到几分钟取决于场景数和外层迭代次数。如果对YALMIP不熟也可以直接用MATLAB的quadprog和mosek接口写但那样代码量会翻倍。YALMIP在处理半定约束时特别省事尤其是二阶矩矩阵不等式的写法手写要维护一大堆矩阵索引YALMIP直接声明矩阵变量就行。3.4 从IEEE6到IEEE118的数据扩展IEEE6节点系统只有3台机组、24个时段外层变量约72个连续变量和72个整数变量代码跑起来非常快适合验证模型逻辑。IEEE30节点系统有6台机组变量数翻倍但依然在Gurobi的可控范围内。真正让我头疼的是IEEE11854台机组、24个时段如果全部建模机组组合变量规模超过2600个加上爬坡和最小启停时间约束外层迭代一次就要好几分钟。我的处理方式是分三步先用一个简化机组组合模型跑一次固定少数关键机组的启停状态再把固定后的启停状态带入完整DRO模型只优化连续出力最后进行一轮启停邻域搜索看是否有必要改变个别机组的启停状态。这样牺牲了一部分全局最优性但能获得一个工程上可接受的近似最优解。在电力市场实际应用中市场出清时间窗口有限这种“先粗后精”的两阶段策略远比硬啃大规模MILP实用。4. 三种节点系统下的调度结果和风险指标4.1 IEEE6基准模型的保守性肉眼可见在IEEE6节点系统上我以α0.95、λ0.5、ε10.1、ε21.2为例做了对比实验。表格如下模型期望成本最坏情景成本CVaR_95CPU时间确定性点预测102.8k128.4k119.1k0.8s随机期望值103.2k126.9k116.5k15.3sDROCVaR104.1k119.8k108.6k47.6s注意这里成本单位是“千美元”只是针对我构造的模拟电价场景不是真实市场数据。关键结论是DROCVaR让期望成本从102.8k上升到104.1k只增加了1.3%但CVaR_95从119.1k降到了108.6k下降了8.8%。最坏情景成本也从128.4k降到119.8k。这意味着模型用很小的期望成本牺牲换来了尾部风险的大幅压缩。在调度曲线上这个变化非常直观。确定性模型在电价最低时段依然安排了机组运行因为它没有预见到极端低价会长期持续DRO模型则在尾部风险高的时段直接压低甚至关停机组即使这些时段在期望情景里看起来还有利润。4.2 IEEE30风险系数λ怎么改变调度策略IEEE30节点系统比IEEE6复杂不少机组数量和负荷曲线都更接近实际。我固定α0.95、ε10.1、ε21.2把λ从0逐渐增加到5观察调度策略变化。λ取值期望成本CVaR_95高价时段购电占比0208.6k242.3k65%0.2209.8k235.1k59%0.5211.2k222.7k51%1.0214.9k213.4k44%2.0219.7k207.6k38%5.0228.3k201.2k33%可以看到λ从0增加到5时期望成本从208.6k上升到228.3k增幅接近9.4%但CVaR_95从242.3k下降到201.2k降幅达到17%。高价时段购电比例从65%降到33%说明模型在风险厌恶程度提升时会主动放弃价格尖峰时段的购电机会把电量转移到本地机组或者更便宜的时段。这条曲线对于实际操作特别有意义。调度员不需要机械地选一个λ而是可以通过λ-风险曲线在期望成本和尾部风险之间找平衡点。如果企业风险承受能力低就选曲线右端如果更看重平均利润就选左端。4.3 IEEE118计算规模与折中IEEE118是我这次实验里最接近实际规模的一个系统。54台机组、24时段、500个价格场景直接用完整MILP跑会非常吃力。我采用固定启停连续出力两阶段策略后单次外层迭代从无法收敛降到几百秒内。从结果看DROCVaR模型在IEEE
)
