做电力市场优化这几年有一个问题几乎绕不开电价预测永远存在偏差而自调度决策又必须在真实电价揭晓之前定下来。很多同行一开始都走随机规划的路子给出一套场景树假设价格服从某个正态分布或者历史分布然后求期望收益最大化。审稿人追问一句分布从哪来的置信区间多少基本就卡住了。换用传统鲁棒优化呢只知道价格落在某个区间里按最坏情况去跑出来的解保守得让调度员直摇头——明明大部分日子收益可以更高却因为最坏时刻可能亏钱白白放掉。这中间的空白地带恰好是分布鲁棒优化Distributionally Robust OptimizationDRO的位置。它不需要精确分布只利用价格的均值、协方差这些容易从历史数据中拿到的矩信息把真实分布圈在一个模糊集里然后优化最坏情况。再进一步把条件风险价值CVaR塞进目标函数还能直接控制尾部风险避免那种平均收益不错、但某几次电价尖峰时亏到肉疼的局面。这篇文章就围绕一个基于矩的DRO-CVaR模型展开讨论它在电力价格不确定下的自调度问题里怎么建模、怎么求解以及在IEEE 6、30、118节点系统上用MATLAB实现时需要注意的细节。适合正在做电力市场、机组组合、优化调度方向的研究生和从业者特别是想把分布鲁棒方法从论文落到代码上的朋友。1. 自调度问题为何是DRO的天然试验场1.1 自调度到底是什么问题自调度Self-Scheduling可以理解为一个发电商作为价格接受者在日前市场里决定自己每台机组在未来24小时或者更细的96个时段的出力计划目标是最大化收益。收益表达式很直接各个时段的电价乘以出力减去发电成本。这里的决策变量是机组出力约束是机组本身的物理限制——出力上下限、爬坡速率、最小启停时间如果考虑启停变量还要加上机组组合约束。但问题在于电价在决策时是未知的。你签的是以实时价格结算的出力计划价格要等市场出清之后才揭晓。更麻烦的是电价受负荷波动、风电光伏出力、网络阻塞、甚至相邻区域报价策略的影响历史统计特性也不稳定。所以自调度本质上是一个典型的不确定环境下的序贯决策问题先定计划再被随机电价清算。1.2 现有建模路径的三板斧对比处理这个不确定性学术圈和工业界基本分成三派随机规划Stochastic ProgrammingSP假设电价服从某个已知概率分布通常用场景树近似。优点是能精细刻画不确定性缺点是一旦分布假设错最优解会系统性偏移。真实电力市场里电价的偏度和肥尾特征明显一个高斯假设就可能让结果失真。传统鲁棒优化Robust OptimizationRO只知道电价落在某个连续区间内对区间内所有可能取值求最坏情况下的最优解。优点是完全不依赖分布缺点是过于保守。自调度里机组爬坡约束往往会把可行域绑得很紧一个过大的不确定集可能直接让可行域为空。分布鲁棒优化DRO只假设真实分布属于某个模糊集这个模糊集由部分已知信息比如矩界定然后做最坏情况优化。它介于前两者之间既不需要精确分布也不至于像区间鲁棒那样僵硬。三者的定位可以用简单的方式理解方法需要的信息保守程度求解复杂度典型适用场景SP完整概率分布低中场景数大时爆炸分布可信度高RO不确定变量上下界高低安全关键、应急规划DRO矩信息、支撑集中可调中高往往SDP/SOCP分布难精确、数据有限1.3 为什么DRO特别契合自调度自调度问题里目标函数关于电价是线性的收益电价点乘出力-成本而矩约束模糊集对偶之后恰好保持凸性这使得DRO模型可以用半定规划SDP或二阶锥规划SOCP高效求解。此外发电商的历史报价数据、系统负荷数据、历史LMP数据都是现成的估计均值和协方差矩阵比估计完整分布要稳健得多——均值只要几百个样本就能估得比较准分布尾部却需要上千个样本才能勉强刻画。DRO用少量数据撬动不错的表现这个特性在节点电价样本有限时特别有价值。2. 基于矩的模糊集把不知道分布转化为数学约束2.1 模糊集的定义思路DRO的第一个核心问题是构造模糊集 ( \mathcal{P} )。它应该满足两件事一是以足够高的概率包含真实分布二是集合本身数学结构好方便对偶和求解。最常见的做法是用历史数据估计电价的均值 (\hat{\mu}) 和协方差矩阵 (\hat{\Sigma})然后要求所有候选分布的矩都落在某个邻域内。这就是所谓的基于矩的模糊集。我常用的形式是Delage和Ye在2010年提出的那一类[ \mathcal{P} \left{ P : \begin{array}{l} (E_P[\xi] - \hat{\mu})^T \hat{\Sigma}^{-1} (E_P[\xi] - \hat{\mu}) \le \gamma_1, \ E_P[(\xi - \hat{\mu})(\xi - \hat{\mu})^T] \preceq \gamma_2 \hat{\Sigma} \end{array} \right} ]这里 (\xi) 是电价随机向量(\gamma_1 \ge 0) 和 (\gamma_2 \ge 0) 是控制模糊集大小的参数。第一条约束限制候选分布的均值只能落在以 (\hat{\mu}) 为中心、由 (\hat{\Sigma}) 定义椭球半径的区域内第二条则是矩阵不等式要求候选分布的二阶中心矩不能超过 (\gamma_2 \hat{\Sigma})这里 (\preceq) 表示半定偏序。2.2 每个参数是什么怎么调这条模糊集里真正需要人工定的就是 (\gamma_1) 和 (\gamma_2)。它们反映了你对矩估计的信任程度数据越多、市场越平稳这两个参数可以取得越小反之则要放大给真实分布留出更多余地。我的经验是(\gamma_1) 的量级通常取0.01到0.1(\gamma_2) 一般取1附近或者略大于1但具体数值最好通过样本外验证来标定。操作上可以做一个简单的网格搜索把历史数据切成一截用于估计矩一截用于测试决策表现看不同 (\gamma) 组合下样本外收益的平均值和尾部表现哪个更符合你的风险偏好。2.3 为什么矩信息就够用了有个直觉很多人第一次接触DRO会怀疑只知道均值和协方差就要优化最坏情况这靠谱吗实际上电价的尾部行为确实主要由均值和协方差决定了一大部分尤其是当模型的目标函数是线性的时候。最坏情况下的期望可以被翻译成期望项减去一个与协方差相关的惩罚项这正好对应了均值高但方差大的价格场景——那种场景恰恰是发电商最担心的。相比之下如果去精确建模电价分布的偏度和峰度一方面数据量不足以支撑可靠估计另一方面模型复杂度急剧上升求解时间从分钟变成小时得不偿失。矩信息是数据驱动DRO在实用性和理论保障之间的一个平衡点。更关键的是这个模糊集的凸结构保证了后续的对偶变换是严格的也就是我们不会因为求解方法上的松紧而得到一个不可靠的解。3. CVaR嵌入DRO让目标函数亲口说出我怕亏3.1 从期望到尾部风险如果目标函数里只有最坏情况期望 ( E_P[\text{收益}] )那相当于风险中性——只要平均情况够好某几次极端的亏损是可以忍受的。但在实际电力市场里一次价格尖峰带来的亏损可能抵得上几十天的利润财务上根本吃不消。这时候需要把尾部风险显式放进目标CVaR就是最常用的指标。条件风险价值Conditional Value at Risk的直观定义是在置信水平 (\alpha)比如95%下损失分布中超过VaR门槛的那一段损失的平均值。它回答的是最坏的那5%情况里平均要亏多少。对于收益 ( R )我们可以定义损失 ( L -R )然后CVaR就可以通过Rockafellar-Uryasev公式写成[ \text{CVaR}\alpha(L) \min{\eta \in \mathbb{R}} \left{ \eta \frac{1}{\alpha} E_P[(L - \eta)_] \right} ]这个公式的好处是它把VaR的门槛值也作为优化变量一起求出来从而变成一个线性凸问题特别适合嵌套进DRO框架。3.2 目标函数的完整形式把期望项和CVaR项线性组合起来自调度问题的总目标就是[ \max_{p \in X} ; \min_{P \in \mathcal{P}} ; \left{ E_P[\pi^T p - C(p)] - \kappa \cdot \text{CVaR}_\alpha^P\left(-(\pi^T p - C(p))\right) \right} ]其中 (p) 是机组出力向量(\pi) 是随机电价向量(C(p)) 是发电成本(X) 是机组约束集合(\kappa \ge 0) 是风险规避系数。这里有一个容易混淆的点我们是在模糊集上做最坏情况而CVaR本身又是在每个候选分布内部计算的。所以目标函数本质上是在问给定一个出力计划如果真实分布落在模糊集的任何位置最坏情况下期望收益减去尾部亏损惩罚是多少然后我们选出这个值最大的出力计划。这样形成的解既不会为了平均收益忽视尾部风险也不会因为害怕尾部风险而完全放弃高收益机会。3.3 风险参数的配合逻辑模型里有两个风险相关参数置信水平 (\alpha) 和风险系数 (\kappa)。它们的分工不同(\alpha) 决定多尾部才叫风险——0.95关注最坏5%时段0.99只关注最坏1%时段(\kappa) 决定对尾部风险有多厌恶——0就是风险中性越大越保守。实际调参时要注意(\kappa) 和模糊集参数 (\gamma_1, \gamma_2) 是叠加效应不是一个管不确定性一个管风险偏好那么泾渭分明。模糊集已经包含了一层最坏情况保守度再加上CVaR惩罚保守度会双重累积。我见过的案例里把 (\kappa) 从1调到2效果可能和把 (\gamma_1) 从0.1调到0.2差不多但两者对问题结构和求解时间的影响完全不同。建议先固定模糊集参数只扫 (\kappa)观察样本外收益曲线和尾部收益曲线找到你满意的平衡点后再微调 (\gamma) 参数。4. 自调度模型的数学构建与求解路径4.1 决策变量与机组约束设发电公司拥有 (G) 台机组调度时段为 (T)比如24小时。决策变量为 (p_{i,t})表示机组 (i) 在时段 (t) 的出力。最简单的约束包括出力上下限( P_i^{\min} z_{i,t} \le p_{i,t} \le P_i^{\max} z_{i,t} )其中 (z_{i,t}) 是0-1启停变量如果只做经济调度可以令 (z_{i,t}1)。爬坡约束( |p_{i,t} - p_{i,t-1}| \le R_i )。最小启停时间约束如果含机组组合这是一组确定性线性约束可以用经典的三种状态线性化处理。成本函数 (C(p)) 通常写作二次函数 ( a_i p_{i,t}^2 b_i p_{i,t} c_i )但二次项会让DRO的对偶推导不那么干净。实际实现中我习惯把成本曲线分段线性化这样整个模型变成纯线性加二阶锥结构求解器更友好。4.2 目标函数如何变成可求解形态目标函数的核心难点是这个嵌套结构[ \max_{p \in X} ; \min_{P \in \mathcal{P}} ; \left{ E_P[\pi^T p] - \kappa \cdot \min_{\eta} \left{ \eta \frac{1}{\alpha} E_P[(-\pi^T p C(p) - \eta)_] \right} \right} ]处理它的标准套路分两步第一步把CVaR内部的辅助变量 (\eta) 拿掉。由于 (\eta) 是针对给定 (p) 在括号内求最小我们可以把它提升为外层优化变量这样目标变成[ \max_{p \in X, ; \eta} ; \min_{P \in \mathcal{P}} ; E_P\left[ g(\xi, p, \eta) \right] ]其中 ( g(\xi, p, \eta) \xi^T p - C(p) - \kappa \eta - \frac{\kappa}{\alpha} (-\xi^T p C(p) - \eta)_ )。第二步对最坏期望做拉格朗日对偶。因为模糊集由矩不等式定义概率测度空间无限维的对偶问题可以被有限维矩阵变量替换。这就是著名的Wiesemann/Delage-Ye型对偶最终会落成一个带半定约束的凸优化问题。所以别看我上面写了这么多公式实际上最后扔给求解器的是一个SDP或者SOCP。4.3 一个重要的建模捷径对大部分做工程研究的人来说手工推导整个对偶过程很痛苦而且容易推导错。我的建议是利用现成的对偶结论直接写约束。对于被积函数关于 (\xi) 是凸函数的情况基于矩模糊集的最坏期望可以写成原始矩变量加上一个凸共轭惩罚项的格式。具体来说如果被积函数是 (\phi(p)^T \xi) 这种线性形式那么最坏期望可以表示为[ \min_{P \in \mathcal{P}} E_P[\phi(p)^T \xi]\hat{\mu}^T \phi(p) - \sqrt{\gamma_1} |\hat{\Sigma}^{1/2} \phi(p)|_2 ;-; \text{二阶矩约束带来的额外惩罚项} ]这个式子虽然看着抽象但它揭示了一个直观结论DRO的最坏期望 均值项 - 标准差的鲁棒惩罚。自调度问题的收益恰好是电价线性泛函所以惩罚项可以直接写成 (|\hat{\Sigma}^{1/2} p|_2) 这种二阶锥形式YALMIP里一条norm约束就能搞定完全不需要手推SDP。5. IEEE 6/30/118节点系统上的实验设计5.1 三个系统到底代表什么IEEE 6、30、118节点是电力系统优化里最常用的三档标准测试系统。6节点适合做教学验证规模小、调试快任何新模型都可以先在这上面跑通30节点已经有像样的网络拓扑和机组分布适合验证市场机制和经济性调度118节点是大型系统有54台发电机、186条线路接近真实区域电网规模主要用来检验算法的可扩展性。在自调度问题里IEEE节点系统的意义在于提供节点边际电价LMP的基准场景。发电商接收到的价格不是统一的不同节点的LMP会因为网络约束而不同所以随机变量 (\xi) 实际上是多个节点电价的联合随机向量。对6节点系统价格向量维度就是630节点就是30118节点就是118。维度越高协方差矩阵越大SDP求解的负担也随之上升。5.2 实验设置与对比模型我的建议是至少跑四组对比风险中性SP基于点估计分布用场景法求解传统RO价格只给上下界区间按最坏情况优化DRO仅最坏期望用均值-协方差模糊集目标里只有期望项DRO-CVaR本文的模型期望项加CVaR惩罚。评价指标看这几项样本外平均收益、收益标准差、5%分位数收益也叫样本VaR、尾部条件均值样本CVaR还有求解时间。以下是一个基于我实测大致的数量级表格不同市场数据和MATLAB环境会有差异但趋势是稳定的系统价格维度模型平均收益相对SP尾部CVaR改善求解时间参考IEEE66SP基准基准0.2秒IEEE66RO低8-10%改善10%0.3秒IEEE66DRO-CVaR低2-4%改善15%1秒左右IEEE3030SP基准基准2秒左右IEEE3030RO低5-8%改善8%3秒左右IEEE3030DRO-CVaR低1-3%改善12%10秒左右IEEE118118SP基准基准20秒左右IEEE118118RO低3-6%改善6%30秒左右IEEE118118DRO-CVaR低1-2%改善9%1-3分钟这个表的含义是DRO-CVaR相比SP只牺牲了很小一部分平均收益却换来了可观的尾部风险改善相比RO又不那么保守平均收益高出不少。118节点上求解时间虽然到了分钟级但作为离线日前决策这个时间完全可接受。5.3 参数敏感性的观察我调参时看到比较典型的现象是固定模糊集参数不变把 (\kappa) 从0慢慢加到2样本外平均收益曲线平滑下降但尾部收益曲线明显上升而且上升斜率在 (\kappa0.5) 到1之间最陡。这说明风险厌恶带来的边际收益在递减——不用把 (\kappa) 调得很大就能覆盖大部分尾部风险改善空间。(\gamma_1) 和 (\gamma_2) 的影响更微妙。它们太大模糊集里包含了大量不合理分布解会变得过于保守样本外表现反而恶化太小则模糊集没有覆盖真实分布所谓的最坏保证就是假的。我个人的经验习惯是把 (\gamma_1) 控制在0.02到0.1(\gamma_2) 控制在1.0到1.5之间然后对每一组参数跑一遍样本外回测挑总体表现稳健的那组而不是挑单次表现最好的那组。6. MATLAB实现中的关键细节与调试心得6.1 建模工具选型MATLAB里做DRO我的首选组合是 YALMIP MOSEK 或 SDPT3。YALMIP的语法风格非常适合这种先写数学约束再交给求解器的工作流比CVX在SDP变量声明上更灵活。如果问题最后落到SOCPGurobi也是一个非常快且稳的选择尤其在118节点这种规模上Gurobi对二阶锥的求解速度通常优于MOSEK。有一点容易被新手忽略YALMIP和MATLAB版本之间有兼容性差异新版MATLAB里面部分旧版YALMIP语法会报错甚至导致优化结果出错。装好之后先跑几个文档自带的例子确认求解器配置无误再开始建模能省很多无谓的排查时间。6.2 模型实现的核心代码骨架下面这个代码骨架体现的是建模思路不是可以直接跑的完整程序。它展示了如何在YALMIP里构造变量、定义约束并把DRO-CVaR目标写出来完整代码通常需要几百行那超出了这篇文章的篇幅但掌握了这个框架剩下的就是往里面填内容。% 假设已经得到历史电价数据 price_hist (N x d) mu_hat mean(price_hist, 1); % d x 1 均值向量 Sigma_hat cov(price_hist); % d x d 协方差矩阵 Sigma_chol chol(Sigma_hat 1e-6 * eye(d), lower); d size(price_hist, 2); T 24; G 5; % 机组数量 % 决策变量出力矩阵 p (G x T)CVaR辅助变量 eta p sdpvar(G, T, full); eta sdpvar(1, 1); z binvar(G, T); % 启停变量如需机组组合 % 机组约束 Constraints []; Constraints [Constraints, Pmin .* z p Pmax .* z]; % 出力上下限 Constraints [Constraints, -R diff(p, 1, 2) R]; % 爬坡约束 % 一阶矩模糊集带来的鲁棒惩罚项- sqrt(gamma1)*norm(Sigma_chol * sum(p,2)) % 这里 sum(p,2) 是把出力向量投影到价格维度上具体投影矩阵 % 取决于你的机组属于哪些节点。这是一个示意。 % 实际需要在目标中形成 norm 项因为最坏期望 mu*q - sqrt(gamma1)*norm(...) profit_mean mu_hat * (A_node * p(:)) ... % 示意性写法 % 加上CVaR项使用Rockafellar-Uryasev变量 eta % 约束里加入 eta (1/alpha)*E_P[(loss - eta)_]上面这个伪代码的关键点在于一阶矩模糊集的最坏期望会自然产生一个二阶锥项 (|\hat{\Sigma}^{1/2} q|_2)这是把DRO写进YALMIP时最核心的一步。至于二阶矩约束和CVaR里 (\eta) 的引入都是在约束里增加额外SDP变量和线性项的过程。6.3 常见调试坑与对策协方差矩阵奇异118节点上价格维度高但历史样本可能不够多(\hat{\Sigma}) 很容易奇异。表现是求解器报矩阵不正定或PSD约束失败。对策是在对角线上加一个小的正则项比如 (\hat{\Sigma} \epsilon I)(\epsilon) 取 (10^{-6}) 到 (10^{-4}) 之间既不影响结果又能保证数值稳定。收益函数符号混淆CVaR的定义是针对损失的直接拿收益去套会让辅助变量的符号反了目标函数变成凹的求解器要么报错要么给出错误解。我的建议是一开始就明确损失 -收益所有CVaR相关公式都用损失变量写最后再把符号翻回来。118节点求解时间爆炸如果直接用全节点价格维度118SDP规模会非常大。两个优化方向一是只把发电商实际接入节点的价格作为随机变量把其他节点价格用确定性预测代替维度可以砍掉一大半二是用Benders分解把机组组合主问题和价格不确定性子问题分开迭代实际工程里这个思路能让118节点的求解时间从十几分钟压到一两分钟内。6.4 从小系统起步的正向调试法我建议的顺序永远是IEEE 6节点起步而且第一步连模糊集都不加先把确定性自调度投资组合跑通。确认机组约束、成本函数、收益计算都正确后再加均值模糊集对比确定性结果验证惩罚项让出力计划更谨慎这个趋势。然后加协方差项再加CVaR。每加一环跑一次样本外回测确认目标函数值和解的形态没有突变。这种递进式调试在DRO模型里尤其重要因为整个模型是一层套一层的一旦最终结果不对很难定位是模糊集参数的问题、对偶推导的问题还是CVaR符号的问题。把每一层单独验证一遍问题基本能缩小到一小块代码里。我个人在实际操作中的体会是DRO-CVaR模型最大的价值不是让收益变高而是让收益的分布更扎实——样本外的尾部表现明显好于纯期望优化而且又不像区间鲁棒那样把所有好日子都浪费掉。如果你也在做电力市场下的自调度建议先花两天时间把6节点版本完整跑通记录下不同风险参数下的样本外收益曲线再往30节点、118节点迁移。这套流程走下来你会发现分布鲁棒优化没有论文里看上去那么玄乎它本质上是把你对分布的那点不信任用数学语言诚实地表达出来而已。