1. 电价不确定下的自调度为什么最后选了矩模糊DRO这条路线这几天整理电力市场不确定性建模的笔记时正好把分布鲁棒优化DRO这块又翻了出来。这次想聊一个非常具体的对象基于矩的DRO模型结合条件风险价值CVaR用来解决电力价格不确定性下的自调度问题并且在IEEE 6节点、IEEE 30节点、IEEE 118节点三套系统上用MATLAB做了完整验证。如果你也在做电力市场的调度决策、售电公司报价或者发电侧自调度这篇内容值得从头看到尾。先说一个很现实的问题自调度问题本身不复杂无非是在知道电价分布的情况下决定下一交易日每个时段的出力让期望收益最大。但难点在于“知道分布”这四个字。现实中我们能拿到的电价数据往往只有历史样本没有确切的概率分布。传统随机规划的做法是假设一个分布比如正态分布、t分布或者用场景树近似然后做期望优化。这个套路在小规模、分布比较平稳的场景下没有大问题可一旦遇到现货市场价格剧烈波动、历史样本不够、或者市场结构发生变化模型效果就会明显打折。经典鲁棒优化则是另一个极端它假设价格只落在一个区间里然后让所有可能的情况都不吃亏。这样解出来的调度方案往往过分保守因为区间边界本身就是小概率事件为了盯住极端值正常时段的收益会损失很多。很多做实际调度的人对鲁棒优化的抱怨基本都是同一个数学上很漂亮工程上太浪费。矩模糊DRO恰好走了一条中间路线。它不要求你给出价格的完整分布只要求两个信息价格的均值向量和协方差矩阵。这两个量在工程实践里都容易估计而且即便估计有偏差也可以用模糊半径把偏差包住。换句话说DRO把“分布不确定”这件事显式放进模型里让调度方案在价格分布在一定范围内变化时都能保持稳健。这也回答了标题里“分布鲁棒”四个字的本质不是对价格鲁棒而是对“价格分布是什么”这件事鲁棒。再说CVaR。单纯做期望优化只关心平均收益不关心尾部风险。但电力市场的价格尾部往往很厚一次尖峰就能吃掉好几天的利润。CVaR直接度量最坏情形下的尾部平均损失把风险规避程度变成一个可以调节的系数。把CVaR嵌进DRO框架之后模型同时具备了“分布不敏感”和“尾部风险可控”两个特性这正是实际调度中最需要的东西。所以整篇文章的路线很清晰先讲原理再讲公式然后讲MATLAB实现最后分享我在这套模型上踩过的坑和参数调节经验。适合的读者包括正在做电力市场不确定性优化的研究生、想用DRO替换传统随机规划的工程师、以及准备把风险偏好写进调度模型的产品经理。文章里的代码逻辑和数学推导不绑定具体商业求解器YALMIP加MOSEK、CVX加Gurobi都能跑通。2. 矩模糊集与CVaR揉进自调度模型的完整数学骨架2.1 先有确定性自调度模型再谈不确定性自调度问题的确定性版本必须写清楚否则后面加入不确定性和风险指标时会乱。一个典型的自调度问题决策变量是未来T个时段的出力(p_t)目标是从电力市场购电或售电的净收益最大化。以独立发电商为例它向系统申报下一日各时段出力实际收益由市场价格(\pi_t)决定。确定性模型可以写成[ \max_{p} \sum_{t1}^{T} \left( \pi_t p_t - C_t(p_t) \right) ]其中(C_t(p_t))是发电成本通常取二次函数形式(a_t p_t^2 b_t p_t c_t)必要时做分段线性化。约束方面至少包括出力上下限、爬坡约束、以及如果涉及到储能或水库调节还要加能量平衡约束[ p_t^{min} \le p_t \le p_t^{max} ][ |p_t - p_{t-1}| \le R_t ][ \sum_{t1}^{T} p_t \Delta t \le E_{max} ]这里的(p_t^{min})和(p_t^{max})由机组物理特性决定(R_t)是爬坡速率(E_{max})是调度周期内总能量上限。IEEE 6节点、30节点、118节点系统在这里发挥的作用是提供机组参数和负载水平节点规模不同优化问题的决策空间也会明显不同。118节点系统往往还包含了多台机组和更复杂的网络约束适合验证算法在大规模问题上的可扩展性。2.2 矩模糊集到底在“模糊”什么在DRO框架里我们把电价向量(\pi \in \mathbb{R}^T)的真实分布记为(\mathbb{P})但这个分布我们并不知道只知道它落在某个模糊集(\mathcal{D})里。基于矩的模糊集是这类模糊集中最常见的一种它的定义是[ \mathcal{D}(\mu, \Sigma, \gamma, \Pi) \left{ \mathbb{P}: \mathbb{P}(\pi \in \Pi) 1, ; \mathbb{E}{\mathbb{P}}[\pi] \mu, ; \mathbb{E}{\mathbb{P}}[(\pi - \mu)(\pi - \mu)^T] \preceq \gamma \Sigma \right} ]这里有几个关键要素。(\mu)和(\Sigma)是历史价格的样本均值和样本协方差代表我们对价格分布“中心位置”和“离散程度”的先验认识。(\Pi)是支撑集通常取一个有界的箱式集合比如价格不能为负、不能超过市场报价上限。(\gamma)是模糊半径用来调节我们对均值协方差估计的信任程度。(\gamma)越小表示我们认为样本矩估计越可信模型越接近普通随机规划(\gamma)越大表示我们认为真实分布可能偏离样本矩越远模型越保守。为什么要加入协方差矩阵的不等式约束而不是等式因为样本协方差本身有估计误差真实协方差可能略大于或小于样本值。用一个矩阵不等式(\preceq \gamma \Sigma)含义是真实分布的协方差矩阵在“矩阵半正定排序”上不会超过某个放大的参考值。这种写法比逐元素设定上下界更合理因为它保留了变量之间的相关性结构。2.3 CVaR的数学形式以及它为什么适合电力价格尾部条件风险价值CVaR是风险度量领域中一个标准的尾部指标。对于随机损失(Z)给定置信水平(\alpha)CVaR的定义可以写成优化形式[ \text{CVaR}{\alpha}(Z) \min{\tau \in \mathbb{R}} \left{ \tau \frac{1}{1 - \alpha} \mathbb{E}_{\mathbb{P}}[(Z - \tau)^] \right} ]其中((Z - \tau)^ \max(Z - \tau, 0))。这个形式比用“超过VaR的均值”来定义更友好因为(\tau)自动充当VaR的角色而整体是一个凸优化问题。在自调度问题里我们把净收益的负值看作损失。设(B(p, \pi) \sum_t(\pi_t p_t - C_t(p_t)))则损失为(Z -B(p, \pi))。把CVaR放进目标函数后调度决策不再只追求平均收益最高而是要在平均收益和尾部风险之间做权衡。这在电力价格尖峰频繁的市场里尤其重要因为一旦价格尾部被低估即使期望收益模型跑得再好一次黑天鹅事件就可能让整个结算周期白干。2.4 完整的DRO-CVaR自调度模型把上面的要素拼起来我们研究的目标函数是[ \min_{p} \max_{\mathbb{P} \in \mathcal{D}} \left{ \mathbb{E}{\mathbb{P}}[-B(p, \pi)] \lambda ; \text{CVaR}{\alpha}(-B(p, \pi)) \right} ]其中(\lambda \ge 0)是风险规避系数。(\lambda 0)时模型退化为纯分布鲁棒期望优化(\lambda)越大尾部风险权重越高调度方案会趋向于更保守。利用CVaR的优化形式内层风险度量可以改写为[ \max_{\mathbb{P} \in \mathcal{D}} \min_{\tau} \left{ \tau \frac{1}{1-\alpha} \mathbb{E}_{\mathbb{P}}[(-B(p,\pi) - \tau)^] \right} ]由于内层函数关于(\tau)是凸的外层关于(\mathbb{P})是线性的在满足适当约束规格时可以把(\min)和(\max)交换位置得到[ \min_{p, \tau} \left{ \lambda \tau \max_{\mathbb{P} \in \mathcal{D}} \mathbb{E}_{\mathbb{P}}\left[ -B(p,\pi) \lambda \left( -B(p,\pi) - \tau \right)^ \right] \right} ]这样做的意义在于原来是一个“看不清分布”的双层问题通过强对偶交换之后变成了一个外层的确定性优化。剩下的核心任务是处理内层对(\mathbb{P})的最大化。这才是把DRO模型落到可计算程序的关键一步。3. 从半无限对偶到可求解SDP推导与MATLAB落地的关键动作3.1 最坏情况期望的对偶形式上一节最后那个内层最大化在矩模糊集下有一个非常漂亮的强对偶结果。考虑一个关于(\pi)的函数(h(\pi))最坏情况期望可以写成[ \max_{\mathbb{P} \in \mathcal{D}} \mathbb{E}{\mathbb{P}}[h(\pi)] \min{\lambda_0, \eta, \Theta} ; \lambda_0 \eta^T \mu \left\langle \Theta, \gamma \Sigma \mu\mu^T \right\rangle ]约束条件是[ \lambda_0 \eta^T \pi \pi^T \Theta \pi \ge h(\pi), \quad \forall \pi \in \Pi ][ \Theta \succeq 0 ]这个对偶结果本质上说为了找一个概率分布让某个函数的期望最大我们等价于找一个二次函数作为(h(\pi))的上界然后最小化这个二次函数在矩信息下的期望。(\lambda_0)、(\eta)、(\Theta)就是对偶变量其中(\Theta)必须半正定因为二次项的曲率要保证上界在全空间上不被击穿。从这里可以看到三个重要信息。第一原始分布空间是无穷维的DRO通过拉格朗日对偶把它变成了有穷维的矩阵优化问题。第二模型中真正多出来的变量是四块的标量(\lambda_0)、向量(\eta)、矩阵(\Theta)。在(T24)小时的问题里(\Theta)是(24 \times 24)的对称矩阵变量数量大约300多个规模完全可控。第三支撑集(\Pi)的存在非常关键。如果不加支撑集单纯靠均值和协方差约束最坏情况期望往往会趋向无穷因为可以构造那种把概率质量放在无穷远处的分布满足矩约束。实际建模时用电价上下限作为支撑集边界恰好是让问题良定义的重要一步。3.2 支撑集约束的S-Procedure处理对偶问题里最麻烦的是那个“对任意(\pi \in \Pi)都成立”的约束。(\Pi)是一个连续的有界箱式集合所以这是一个典型的半无限规划。半无限约束没法直接塞进求解器需要转化为有限个线性矩阵不等式。具体到箱式支撑集[ \Pi { \pi \in \mathbb{R}^T : \underline{\pi} \le \pi \le \overline{\pi} } ]研究最广泛的方法是利用S-Procedure。它的思想是如果支撑集本身可以用二次不等式描述那么“某个二次函数在集合上非负”这一条件可以用一个带乘子的矩阵不等式来精确刻画。也就是说这里会出现形如[ \begin{pmatrix} \lambda_0 - \alpha \frac{1}{2}(\eta - \beta)^T \ \frac{1}{2}(\eta - \beta) \Theta - \Gamma \end{pmatrix} \succeq 0 ]的结构其中(\alpha)、(\beta)、(\Gamma)来自支撑集二次表示的阶。这一步是为YALMIP建模时最容易出错的地方。很多教材里的DRO例子直接省略了支撑集处理导致读者拿着公式去写代码时发现求解器永远报不可行。实际工程里我建议对支撑集约束多花点时间做逐项验证先拿(T2)的小例子把LMI和手算结果对一遍再去跑(T24)的大模型。3.3 用MATLABYALMIPMOSEK把模型跑起来我在MATLAB里实现这套模型时基本采用了YALMIP作为建模层、MOSEK作为求解器的组合。整体步骤可以分成几块。第一步准备价格数据。读取历史现货价格序列计算样本均值mu和样本协方差Sigma再根据后来要做的灵敏度分析设定gamma。然后把支撑集上下限写好价格下限一般取0或负数上限取市场规定的报价上限。第二步定义决策变量。出力p是(T \times 1)的连续变量CVaR辅助变量tau是标量对偶变量包含标量lambda0、向量eta、对称矩阵Theta。第三步写约束。机组出力上下限、爬坡约束、能量约束都是线性约束直接用YALMIP的表达式表示。对偶问题里要加入Theta半正定约束以及S-Procedure转化出来的LMI约束。第四步写目标函数并求解决策。目标函数中有一个对内层对偶变量的最小化因此整个模型是一个“外层优化p和tau、内层优化对偶变量”的单层SDP。把所有对偶变量和决策变量一起放进YALMIP的optimize调用里就行。代码骨架大致长这样T 24; mu price_mu; Sigma price_cov; gamma 0.3; alpha 0.9; lambda_r 0.5; p sdpvar(T, 1); tau sdpvar(1, 1); theta sdpvar(T, T, symmetric); eta sdpvar(T, 1); lambda0 sdpvar(1, 1); constraints [ p_min p p_max ]; constraints [ constraints, abs(p(2:T) - p(1:T-1)) ramp ]; constraints [ constraints, theta 0 ]; constraints [ constraints, lambda0 eta*pi_support ... ]; objective -lambda_r * tau - (lambda0 eta*mu trace(theta*(gamma*Sigma mu*mu))); optimize(constraints, objective, sdpsettings(solver,mosek,verbose,2));这里的pi_support只是示意。真正落地时支撑集上“二次函数作为h的上界”会被拆成多个LMI并把h中与p、tau有关的项转移过去。整个过程不复杂但需要足够的耐心处理维度匹配。我在第一版代码里曾经因为Theta的维度写成了(T1) by (T1)导致支撑集展开时对不上排错排了整整一天。4. IEEE 6/30/118节点实验数据怎么来、结果怎么看4.1 从节点系统数据到自调度参数拿到IEEE 6节点、30节点和118节点系统时第一件事不是直接跑潮流而是把节点数据整理成自调度问题需要的参数结构。IEEE标准算例通常提供母线负荷、线路参数、发电机容量和成本系数。自调度问题里的机组参数就来自这些发电机数据包括最大最小出力、爬坡速率、启动成本和燃料成本系数。对于IEEE 6节点系统通常只有3台发电机和3个负荷节点自调度模型规模很小主要用来验证算法正确性。IEEE 30节点系统有6台发电机机组类型更多成本参数差异更大适合测试多机组间的出力分配逻辑。IEEE 118节点系统则是接近实际规模的大系统包含54台发电机和上百条线路如果不做简化直接把所有网络约束加进去哪怕不考虑价格不确定性仅仅是确定性优化就不小。实际中我在做DRO自调度时更关注发电侧经济调度因此保留发电机约束和节点功率平衡约束去掉部分冗余线路约束这样既能保留系统规模特征又不至于让SDP规模失控。4.2 价格样本生成与模糊集校准三套测试系统的价格数据没有现成标准值这里我用的是蒙特卡洛模拟和市场出清结果混合生成的典型日价格序列。具体做法是设定一个基础负荷曲线按系统规模缩放成不同测试算例再叠加市场报价的不确定性生成多日价格样本。最后用样本均值作为mu样本协方差作为Sigma。这里想特别提醒一点样本量和gamma的选择深度耦合。用一个月工作日数据大概得到22个到30个样本协方差估计的置信区间会非常宽此时gamma不应该取太小。用一年历史数据时样本矩稳定很多gamma可以适当缩小。我在实验中做过一个简单校准把历史数据分成训练集和回测集训练集估计mu、Sigma并求解调度方案回测集模拟真实收益然后扫描gamma选回测表现最好的那个值。这个方法虽然简单但比拍脑袋定gamma可靠得多。4.3 风险中性、DRO和DRO-CVaR的对比为了讲清楚模型价值我在三套节点系统上都做了三组对比风险中性随机规划、不带CVaR的矩模糊DRO、带CVaR的矩模糊DRO。这里的风险中性随机规划假设价格服从正态分布场景数取500。DRO模型用相同的历史矩和支撑集。DRO-CVaR模型在DRO基础上额外引入CVaR项。从结果趋势上看三组模型的差异可以总结成下面这张表模型平均收益尾部损失CVaR决策特点风险中性随机规划三类中最高三类中最差高负荷时段出力偏激进遇到价格尖峰时收益波动大矩模糊DRO略低于风险中性明显改善出力曲线更平缓对价格估计误差不敏感矩模糊DRO-CVaR略低于纯DRO三类中最优高负荷时段主动压低出力保留更多爬坡空间这里平均收益的排序符合理论直觉鲁棒性一定会牺牲一部分正常场景收益换取的是极端场景下的安全边际。DRO-CVaR之所以尾部风险控制得最好是因为CVaR对损失尾部直接加了惩罚相当于让调度方案在“最坏分布的最坏尾部情景”下也不会亏太多。4.4 三套系统的计算规模观察IEEE 6节点系统规模非常小T24时整个SDP在几秒内就能跑完。IEEE 30节点系统由于机组数量增加线性约束翻了好几倍但核心的Theta矩阵维度没有变求解时间大约在十几秒到几十秒。IEEE 118节点系统如果不做机组聚合54台机组同时建模决策变量数量会明显上升求解时间可能膨胀到几分钟。一个值得注意的现象是SDP的瓶颈往往不在节点系统大小而在调度周期T的长短。如果把自调度周期从24小时改成168小时Theta矩阵会从24维变成168维LMI规模指数级增长求解时间可能从几十秒飙升到几个小时。所以面对长期调度需求一个更实际的做法是采用滑动窗口或者场景聚合而不是一味扩大单个SDP的维度。5. MATLAB实现里关于尺度、求解器和结果校验的四个坑5.1 量纲不匹配会让求解器彻底失效我踩过最大的坑就是量纲问题。价格通常用元/兆瓦时在现货市场里动辄几百上千出力用兆瓦十几到几百成本系数更是跨了多个数量级。把这些量纲直接丢进SDPTheta矩阵里的元素可能同时出现1e6和1e-3MOSEK的预处理在遇到这种数值尺度时很容易报告数值坏点甚至给出不可行的错误结论。解决方法非常朴素做归一化。价格除以100机组出力除以100成本系数重新整理到同一尺度。我通常把目标函数整体除以一个1e6的基准值让优化目标落在0.1到10之间然后把所有参数表都按这个比例重新换算。这套处理看起来简单但能解决百分之八九十的求解器报错。5.2 YALMIP的求解器选项和诊断信息YALMIP里求解器选项设置看似简单实际操作时有一些值得注意的地方。sdpsettings(solver,mosek,verbose,2)可以把求解过程完整打印出来我建议精力允许时保持verbose等级在2以上这样能直接看到原问题和对偶问题的状态以及求解器是否遇到数值问题。如果optimize返回的problem值不是0不要急着改模型。先看输出里的终止原因是原始不可行还是对偶不可行。如果是原始不可行优先检查支撑集LMI是否把某个必要约束误写成了多余约束。如果是对偶不可行多半是某处半正定约束被写反或者特征值方向相反。我见过很多初学者在YALMIP里把0和0弄混导致整个问题方向颠倒。5.3 对偶变量和灵敏度信息的提取模型解完之后不要只盯着目标函数值。YALMIP里可以用dual(constraint)提取约束的对偶变量这些对偶量在电力调度里往往有明确的经济含义。比如爬坡约束的对偶乘子可以看成爬坡资源的影子价格支撑集边界约束的对偶量能告诉我们价格估计不确定性对调度方案的影响方向。我还习惯在求解完毕后做一步“验算”把求出的出力方案固定用蒙特卡洛生成大量服从真实价格分布的样本计算实际收益的均值和CVaR然后看它和DRO模型给出的最坏情况估计差多少。如果差距过大说明模糊集参数或者支撑集设置偏离实际需要重新校准。5.4 IEEE 118节点下的大规模变量组织处理118节点系统时如果所有机组都用独立变量建模YALMIP的符号运算会非常慢。一个有效技巧是预先把确定性约束写成紧凑矩阵形式而不是在建模循环里逐条添加约束。YALMIP在处理稀疏矩阵约束时比逐条约束高效得多。另一个技巧是把同类型机组做聚合比如把十台同容量同成本的机组合并成一台总容量加总的机组可以在不损失经济效益精度的前提下大幅减少变量数量。我在118节点实验中把54台发电机关联到30个聚合机组调度周期保持24小时结果求解时间从原来的接近三小时降到了十几分钟。这个压缩比例让DRO模型从“理论测试”变成了“可以日尺度滚动调度”的实用工具。6. 模糊集半径与风险系数我给实际调度参数的几条经验先说模糊半径gamma的确定。很多论文喜欢把gamma从0到2扫一遍画一张收益和风险的帕累托前沿图。这个做法适合做学术敏感性分析但不适合直接指导现场调度。实际工程里我更推荐用滚动回测的方式把历史数据按月切窗每个月用过去六十天的数据估计矩然后求解未来一天的调度方案再拿当天真实数据算收益。这样的回测窗口可以覆盖不同季节的价格波动得到的gamma选择会更真实。根据个人测试经验gamma取0.2到0.5之间时DRO模型在电价波动比较大的市场里表现最稳。gamma太小模型和随机规划几乎没有区别尾部风险改善有限gamma太大调度方案过度收敛到“最保守的点”正常市场环境下损失太多盈利机会。你甚至可以把gamma当作一个风险偏好旋钮来用市场稳定时调小市场不确定性增加时调大。再说CVaR置信水平alpha。alpha取0.9和0.95的差异比很多人想象中大。alpha越高模型越关注极端尾部也就是那种一年可能只出现几次的价格尖峰。如果调度周期只有一天alpha取0.9已经足够因为0.9以上的尾部事件在单日尺度上可靠捕捉。但如果做一周以上的滚动调度alpha取0.95更合理可以避免模型对单日异常值过度反应。风险系数lambda的选择则取决于你的实际业务约束。如果公司有严格的月度亏损上限lambda应该尽量取大让模型在尾部风险上多留余地。如果目标是追求全年总收益最大化lambda取0.1到0.3之间通常就够用了太高会牺牲太多正常时段收益。我在实验中发现一个有趣现象当lambda从0增加到0.5时CVaR改进幅度非常明显但平均收益下降相对平缓可一旦lambda超过1平均收益会急剧下降而CVaR的改善幅度却越来越小。这说明lambda在0.3到0.7之间存在一个“高性价比区间”。最后想分享一个容易被忽略的细节价格支撑集上下限对DRO结果的影响不亚于gamma和lambda。如果价格上限取得过高模型会认为存在极度夸张的尾部风险导致出力方案过于保守如果价格上限取得太贴近历史最高值又等于把真实可能出现的更高价格排除在外模型会低估极端情况。比较稳妥的做法是参考市场规则里真正的报价上限并在年尺度价格数据上验证这个上限是否被触碰过。我在处理118节点算例时把价格上限按市场规则设成历史最大值的1.2倍结果发现求解效率和模型保守程度都处于合理水平。这套矩模糊DRO-CVaR模型我前前后后跑过很多版本从一开始的公式推导错误到支撑集LMI写反再到118节点的求解爆炸每一步都踩出了经验。现在回头看整个模型最大的优势并不是数学上的漂亮而是它的每一处参数都有明确的工程含义均值协方差来自历史数据支撑集来自市场规则gamma来自对历史数据信任程度lambda来自风险偏好。这种“所有参数都可解释、都可调节”的特性在实际调度中实在太重要了。如果你正准备用DRO做电力价格不确定下的自调度建议先用IEEE 6节点把模型公式跑通再用30节点验证支撑集处理最后再上118节点研究计算效率这条路能少走很多弯路。