电力市场自调度中基于分布鲁棒优化与CVaR的建模和MATLAB实现

做电力市场优化的同行应该都有这种体验：你辛辛苦苦把机组约束、网络约束、投标策略都建好模，最后发现最不可控的变量是明天的电价。它不给你面子，负荷预测偏了它涨，新能源大发它跌，某条通道检修它直接飙升。最近我把一组分布鲁棒优化（DRO）模型在IEEE 6、30、118节点系统上完整跑了一遍，模型用基于矩的模糊集去刻画电价不确定性，并在目标函数里揉入了条件风险价值（CVaR），用来处理发电商自调度问题。这篇博文把整个建模思路、对偶推导、MATLAB实现细节和踩坑记录一次性讲清楚，适合正在做电力市场优化、风险评估，或者刚入行想搞懂DRO到底怎么落地的人。

1. 为什么自调度问题绕不开价格不确定性这道坎

1.1 自调度在电力市场里到底解决什么问题

在现货市场环境下，发电商需要在日前市场申报第二天的发电曲线，系统出清后按市场价结算。所谓自调度，就是发电商在给定自身机组物理特性的前提下，决定各时段启停和出力，使期望收益最大。注意这和经济调度有本质区别：调度中心关心的是全系统社会福利，发电商关心的只是自己的利润。

如果电价是确定常数，这就是个混合整数线性规划，任何一个商业求解器都能秒解。但现实是电价本身就是随机变量，而且它的分布你并不清楚。尤其在新能源渗透率高的系统里，价格尖峰、负电价、爬坡段飙价已经成为常态。你提前一天的报价决策，可能在一个极端价格场景下就把整个月的利润亏回去。

1.2 确定性、随机规划和纯鲁棒优化为什么都不够

确定性方法最简单，用点预测电价去优化，算完再做个敏感分析。问题也很明显：点预测误差是必然存在的，而自调度问题的收益函数对价格非常敏感，一旦实际价格偏离预测，要么少赚，要么亏损，某些极端场景下连启动成本都收不回来。

随机规划把电价建模成已知概率的场景树，然后求期望利润最大化。思路没问题，但有一个致命前提——你必须知道真实分布。现实里你只有历史样本，没有那个“真实分布”。如果你假设正态分布而真实电价是重尾分布，解出来的方案在真实环境下反而很差，这就是分布误差带来的“不合意最优”。

纯鲁棒优化走向另一个极端：把所有可能的价格都装进一个不确定集合，在集合内最坏情形下做决策。好处是很安全，代价是过度保守。为了应对1%概率的极端尖峰，模型可能让机组在平时也压低出力，导致期望收益被砍掉一大截。

DRO站在中间：我只知道电价分布的均值、协方差这些矩信息，真实分布在满足这些矩条件的集合里，我就在集合中找最坏情况求最优。它比随机规划抗分布误差，又不像鲁棒优化那样把极端场景和常规场景一视同仁。这就是“分布鲁棒”四个字的来历。

2. 基于矩的模糊集：用均值和协方差把真实分布圈起来

2.1 模糊集的数学定义和直觉

DRO的核心是模糊集（ambiguity set）。基于矩的模糊集长这样：

[ \mathcal{D} = \left{ P \in \mathcal{P}(\Xi):\ \mathbb{E}_P[\xi]=\mu,\ \mathbb{E}_P[(\xi-\mu)(\xi-\mu)^\top] \preceq \Sigma,\ P(\xi \in \Xi)=1 \right} ]

这里的 (\xi) 是整个决策周期T个时段的价格向量，(\mu) 和 (\Sigma) 从历史样本估计，(\Xi) 是价格支持集。协方差用“(\preceq)（矩阵半正定意义下的小于）”而不是“=”，意思是只要真实分布的协方差不超过这个矩阵上界就算合格，这样做可以把矩估计误差和分布不确定同时包进去。信息量稍微松一点，换来的是模型对“分布到底是什么样”不再那么敏感。

直觉上可以这么理解：你站在市场外面，看不清明天电价到底怎么分布，但你能统计出过去几百天的均价和波动范围。真实分布虽然未知，但总不能偏离这些统计量太远。模糊集就是把所有“统计上和历史数据一致”的分布都圈进来，再求一个最坏情况下的方案。只要真实分布确实在这个圈里，你的决策就不会翻车。

2.2 为什么矩信息是工程上最划算的假设

矩信息是所有分布信息里最容易从历史数据里获得的。电价的均值、协方差有明确的工程意义——均价决定基准收益，方差刻画波动程度。而高阶矩（偏度、峰度）估计对样本量非常敏感，看着一天的8760个小时数据很多，但分时段、分季节之后，每个场景的样本数根本不够，偏度估计出来噪声很大，加了反而坏事。

和Wasserstein距离模糊集相比，基于矩的集合对偶后结构清晰，可以直接提炼成半定约束。Wasserstein模糊集在小样本数据驱动场景下确实更灵活，但对偶里会出现带运输项的结构，和机组组合这种混合整数模型纠缠在一起时，商业求解器很难直接啃下来。实操中，我不建议一上来就上Wasserstein，先把基于矩的模型跑通、跑稳，再根据数据量决定要不要换更精细的模糊集。

2.3 矩估计的三个注意点

第一，用滚动窗口估计，不要用全部历史平均。电力市场有明显的季节性和周内结构，全历史均值会把夏天和冬天的价差糊在一起，均值估计没有代表性。我习惯用过去30到60天、分时段的滑动窗口，再叠加一个简单的季节调整项。

第二，分时段估计。市场电价下午高峰段的方差和凌晨低谷段完全不是一个量级。如果用一个全局协方差矩阵去描述，高峰段的风险会被低估，低谷段又被高估。分时段处理后，协方差矩阵会出现明显的“分段对角”结构，这正好和机组运行特性对应。

第三，协方差矩阵维度和决策时段数相同（(T \times T)），当T=24甚至168时，样本协方差很容易病态。后面求解部分我会详细讲正则化处理，这里先提个醒，别拿到样本就直接算协方差往模型里塞。

3. CVaR在DRO框架里的两种嵌法：惩罚目标还是风险约束

3.1 CVaR和VaR的本质区别

很多资料讲CVaR喜欢直接上公式：

[ \mathrm{CVaR}\alpha(L) = \min{\beta} \left{ \beta + \frac{1}{1-\alpha} \mathbb{E}[(L-\beta)_+] \right} ]

其中 ((L-\beta)_+ = \max{L-\beta, 0})。这个公式被学术界用烂了，但工程上真正重要的是它和VaR的区别。VaR说的是“有(\alpha)的概率损失不超过某个值”，超过之后怎么分布它完全不关心；CVaR说的是“最差那 (1-\alpha) 概率情景下的平均损失”。前者管不住尾部，后者直接把尾部期望管住了。

对发电商来说，CVaR的工程含义非常直观：你做了一百次市场仿真，把最差的五次结果拿出来平均，看平均亏多少钱。把这个值限制住，就相当于管住了“黑天鹅”情景下的平均表现。而且CVaR满足次可加性，是凸的，做优化时能保持问题的可解性——这一点对DRO模型尤其重要，因为凸性意味着你可以放心地做对偶变换。

3.2 目标函数嵌入和约束嵌入

CVaR进入自调度模型有两种常见方式。

方式一是均值-CVaR目标：

[ \max_{x \in X} \ \mathbb{E}[\mathrm{Profit}] - \lambda \cdot \mathrm{CVaR}_\alpha(-\mathrm{Profit}) ]

(\lambda) 是风险权重，越大越保守；(\lambda=0) 时退化为纯期望利润最大化。

方式二是CVaR硬约束：

[ \max \ \mathbb{E}[\mathrm{Profit}] \quad \text{s.t.} \quad \mathrm{CVaR}\alpha(\mathrm{Loss}) \le V{\text{th}} ]

方式一容易调节风险偏好，风险权重连续变化，可以直接画出“利润—风险前沿”；方式二适合企业有明确风险限额的场合。两种方式在DRO下都会形成“最坏分布下的CVaR”，但数学结构不同，实现难度也不同。我在项目里优先用方式一，因为风险前沿曲线对管理层解释模型最直观。

3.3 DRO与CVaR叠加出的双层结构

在DRO里，你已经对真实分布取了外层的最坏情况 (\sup_{P \in \mathcal{D}})，CVaR内部又有一个期望算子。所以DRO-CVaR模型实际上是一个双层最坏结构：

[ \sup_{P \in \mathcal{D}} \ \mathbb{E}P \left[ \mathrm{CVaR}\alpha(-\mathrm{Profit}(x,\xi)) \right] ]

这里不能简单地把CVaR样本平均再求最坏，因为内层CVaR的最优点 (\beta) 也依赖于分布P。需要在模糊集内部同时优化β和期望，这就涉及下一节的minimax交换和对偶化。很多人卡在这一步，以为用YALMIP定义一个CVaR就能和DRO直接拼接，实际上要先把内层LP结构写出来，再做整体对偶。

4. 完整模型怎么落地：从双层最坏情形到半定规划

4.1 机组组合约束不能只写出力上下限

自调度模型的决策变量包括机组启停状态 (u_{i,t})、出力 (p_{i,t})，以及启动/停机状态变量。约束除了出力上下限，还有最小启停时间、爬坡约束、合同电量平衡等。这部分是标准MILP，细节不展开，但提醒几个容易错的地方：启动成本要区分热启动和冷启动，否则检修后的机组启动成本会被严重低估；爬坡约束要同时写“上一时段到本时段”和“本时段到下一时段”两个方向；最小启停时间约束如果写成一连串递推式会引入大量二进制变量，建议用文献里经典的紧凑不等式，避免重复建模。

如果你在118节点上还打算加网络约束，那就不是简单的自调度了，而是考虑节点电价的博弈问题，模型复杂度会指数级上升。这篇项目里默认机组报价和售电量在单节点聚合，网络约束留作扩展方向。

4.2 目标函数：确定性收益减去最坏CVaR

把收益拆成确定性部分和随机部分。确定性部分包括燃料成本、启停成本、固定合同收入；随机部分主要是现货市场收益，它和电价向量 (\xi) 线性相关。目标函数写成：

[ \max_{x \in X, \beta, z} \quad \text{确定项}(x) - \kappa \cdot V(x) ]

其中：

[ V(x) = \sup_{P \in \mathcal{D}} \mathbb{E}P \left[ \min{\beta} \left( \beta + \frac{1}{1-\alpha} \mathbb{E}P[(-\mathrm{Profit}(x,\xi) - \beta)+] \right) \right] ]

(\kappa) 就是前面的风险权重(\lambda)，这里用(\kappa)是想强调它和CVaR内部(\beta)无关。

4.3 对偶化之后为什么变成半定规划

把CVaR的LP形式展开，外层 (\sup) 对模糊集 (\mathcal{D}) 取拉格朗日对偶，minimax交换条件在Slater条件下成立，最后会得到一个关于不确定变量(\xi)的二次齐次表达式。这个表达式要求对任意(\xi \in \Xi)都成立，等价于一个线性矩阵不等式（LMI）。

最终整个模型变成MISDP（混合整数半定规划），或者在某些简化场景下退化为MISOCP / MIQP。对偶后的目标函数会出现类似这样的项：

[ \tau + \mu^\top \theta + \langle \Sigma, Y \rangle ]

其中(\tau, \theta, Y)是对偶变量，(Y)是(T \times T)半定矩阵变量，(\langle \Sigma, Y \rangle = \mathrm{trace}(\Sigma Y))。这个结构是“基于矩”这个选择带来的最大红利：换成Wasserstein模糊集，这里不会出现整齐的半定项，求解器拿它没办法。

4.4 求解策略：整体求解还是Benders

MISDP直接丢给商业求解器是不现实的。Gurobi能解混合整数线性/二次，但SDP约束需要MOSEK这类锥规划求解器，两者耦合起来很难处理。我实际的做法分两档：

6节点系统规模小，直接用YALMIP把对偶后的SDP约束和机组组合一起建模，用MOSEK解根节点，再用割平面枚举处理二进制变量，整体求解没问题。30节点开始，我会固定机组组合后把DRO-CVaR子问题交给MOSEK，外层用枚举或者简单的启发式遍历启停组合，虽然不保证全局最优，但工程上足够。

118节点必须上Benders分解或者外逼近：主问题是机组组合MILP，子问题是给定启停方案下的DRO-CVaR评估（SDP），子问题返回最坏情形的收益和次梯度作为割，不断迭代。这条路我跑通之后，118节点168时段的场景下，单次迭代SDP大约几十秒，总求解时间在几十分钟量级，算可以接受。别指望一行YALMIP代码把MISDP整体丢给求解器——内存会先爆炸。

5. MATLAB实现路线：6节点调通、30节点调参、118节点验规模

5.1 三个测试系统的分工

IEEE 6节点系统机组少、节点少，适合把每个环节跑通，和确定性模型逐项对账，验证DRO-CVaR代码有没有bug。30节点是标准测试系统，机组数和约束复杂度适中，做灵敏度分析最方便，论文里的“利润—风险前沿”多数是在这个规模上画的。118节点是公认的中大型测试系统，机组几十台，约束上千，用来验证算法在真实规模下的可扩展性。建议严格按这个顺序来，不要一上来就跑118，SDP和机组组合耦合在一起，调试难度会叠加到让人怀疑人生。

5.2 环境、求解器和代码结构

我的环境是MATLAB R2022b + YALMIP + Gurobi + MOSEK。YALMIP负责建模，Gurobi解混合整数部分，MOSEK解SDP。只是6节点的话，也可以用YALMIP的CUTSDP把SDP割成线性约束交给Gurobi，但30节点以上别这么玩，割平面层数太多，计算量反而失控。

项目文件建议按模块分目录：

• data/放IEEE系统参数、价格样本、机组参数
• models/放四个模型：确定性、SP、RO、DRO-CVaR
• results/放输出表格和调度曲线
• utils/放矩估计、协方差正则化、结果后处理函数

这样换一个IEEE系统时，只需要替换data目录下的数据文件，模型代码一行不用改。

5.3 核心代码骨架

矩估计和协方差正则化：

% price_samples: T x M 矩阵，T个时段，M个历史样本 mu = mean(price_samples, 2); Sigma = cov(price_samples'); % 正则化，防止协方差矩阵病态 % 1e-4 是经验值，数据噪声大时放到 1e-2 Sigma = Sigma + 1e-4 * diag(diag(Sigma));

DRO-CVaR对偶后的YALMIP骨架：

p = sdpvar(T,1); % 机组出力 u = binvar(T,1); % 启停状态 beta = sdpvar(1,1); % CVaR辅助变量 % 对偶变量 tau = sdpvar(1,1); theta = sdpvar(T,1); Y = sdpvar(T,T,'symmetric'); % 半定约束和线性约束 Constraints = [Y >= 0, ...]; % 补充线性矩阵不等式的具体项 Constraints = [Constraints, ...]; % 机组约束 % 最坏情形CVaR上界（手推对偶后得到） worst_cvar = tau + mu'*theta + trace(Sigma*Y); % 目标：确定性利润 - kappa * 最坏CVaR Objective = base_profit - kappa * worst_cvar;

注意上面的theta和Y不是随便定义的，它们来自对偶推导，每一项都必须和原问题的约束一一对应。我建议不要依赖YALMIP的自动对偶，手推一遍对偶表达式再写代码，否则出错了根本不知道错在哪。

求解器配置：

ops = sdpsettings('solver','gurobi','gurobi.MIPGap',0.01,... 'verbose',2); ops_sdp = sdpsettings('solver','mosek',... 'mosek.SDP.MSK_DPAR_INTPNT_TOL_PFEAS',1e-6);

118节点用Benders分解时，主问题用Gurobi，子问题用MOSEK，两个求解器之间通过割传递信息，每次迭代都重新构造YALMIP模型，不要尝试在同一个模型里同时调两个求解器。

5.4 没有真实电价数据时怎么生成样本

没有现成现货价格数据时，可以用时间序列模型生成样本。我的简单做法：ARIMA拟合日内价格趋势，GARCH描述波动率聚集，再叠加一个尖峰分布做蒙特卡洛抽样。核心目的是让样本具备电价该有的自相关和重尾特征，而不是做精确预测。

如果样本量太少，比如只有20个，μ和Σ估计方差很大，DRO结果会偏保守。一个处理办法是给矩约束加松弛变量，让模糊集自动放大；另一个更省事的办法是干脆采集更长时间窗口，把样本拉到200个以上。样本充足时，基于矩的DRO几乎不会出现“保守到没用的局面”。

6. 结果对比与灵敏度：风险偏好如何改写机组组合

6.1 四种模型在三个节点系统上的结果对比

以一台燃煤机组加一台燃气机组为自调度主体，基于模拟价格样本，对比四个模型的核心指标。数值是相对值，但趋势在三个节点系统上完全一致：

DRO-CVaR牺牲了约9%的期望利润，换来的是和鲁棒优化接近的最坏情形保障，代价远低于纯鲁棒优化28%的利润牺牲。这个“中间地带”正是DRO的核心价值。在118节点系统上，求解时间差异更明显，确定性模型秒出，DRO-CVaR需要Benders迭代，但结果回报是值得的。

6.2 风险权重和置信水平的灵敏度

把风险权重 (\kappa) 从0慢慢拉到1，期望利润单调下降，CVaR单调改善。有趣的是，调度计划会发生结构性变化：(\kappa) 较小时，燃气机组在价格高峰时段满发，追逐尖峰收益；