电力市场里做日前自调度，最难的不是机组组合那套整数变量，而是电价到底怎么建模。你拿着历史场景做随机规划，第二天来个尖峰价格，利润直接被打回原形；改用区间鲁棒优化，又把最乐观的情况全丢掉，调度方案保守到让交易员想骂人。这个项目做的就是另一条路：用分布鲁棒优化（DRO）把电价的不确定性框在一个“矩模糊集”里，再配合条件风险价值（CVaR）去控制尾部损失，让最终的发电商自调度方案既不过分乐观，也不至于保守到没法用。模型用MATLAB实现，分别在IEEE 6、30、118节点系统上做了测试，覆盖面从教学案例一直到中等规模电网，适合正在做电力市场不确定性优化的研究生、工程师参考。
整套东西的核心其实就三件事：第一，如何用均值和协方差构造一个包含足够多“坏分布”的模糊集；第二，如何把CVaR塞进这个最坏情形优化框架里还保持凸性；第三，怎么用MATLAB配合YALMIP和商业求解器，把模型变成一段能跑的代码。下面按这几个层次展开讲。
1. 模型设计思路：为什么要把DRO和CVaR绑在一起
1.1 传统随机规划在电价不确定性下的短板
自调度问题里最常见的是两阶段随机规划：给定一组典型价格场景，让机组在“期望收益最大化”的目标下安排启停和出力。这个方法工程上很好实现，短板也很明显——它只优化了场景期望，基本上不关心场景内部的极端情况。电价序列普遍存在尖峰厚尾特征，100个场景里哪怕有5个极端高价或低价场景，对CVaR这种尾部指标的冲击就被期望值稀释掉了。
另一个问题是场景概率的估计。随机规划假设历史场景概率是真实概率，但市场电价分布会受供需突变、线路阻塞、燃料价格影响，历史概率和未来的真实分布往往对不上。一旦分布漂移，随机规划方案很容易在实际上“跑偏”。
鲁棒优化倒是把不确定性考虑进去了，但它假设价格落在某个区间盒子里，本质上把所有场景一视同仁，不区分“稍微偏离均值”和“极端偏离均值”。用力过猛的结果就是：为了避免最坏情况，机组几乎不敢报高价出力，收益天花板被压得很低。实际调度人员拿到这种方案，往往宁可自己拍脑袋改一版。
1.2 矩模糊集与CVaR的组合价值
DRO的思路介于两者之间：不假设一个精确分布，而是只假设真实分布属于一个“模糊集”——比如把分布的一阶矩（均值）和二阶矩（协方差）约束在一定范围内。你不需要准确知道每个场景的概率，只需要知道“均值大概在这个范围、波动水平大概在那个范围”。这个信息通常用历史数据就能估出来，而且对分布漂移有一定容忍度。
CVaR在这里的作用是给DRO装上“尾部风险雷达”。DRO求的是最坏分布下的期望收益，但这只反映了分布的平均表现，万一模糊集里存在一个尾部极差的分布，期望收益指标不一定能暴露问题。CVaR直接把注意力放在损失分布的尾部，计算的是“最差5%情形下的平均收益”，结合DRO之后，得到的是“在所有可能的坏分布里，最差的那部分场景下，我们还能保住多少收益”。从调度员角度来说，这个指标比纯期望收益直观得多。
说直白点：随机规划是“平均心态”，鲁棒优化是“被害妄想症”，DRO+CVaR是“保留一定的风险警觉，但不至于脱敏”——它把分布的不确定性和收益的尾部风险同时纳入了优化目标。
1.3 自调度问题的目标函数与约束框架
自调度不是传统的最优潮流，它不考虑整个电网的潮流平衡，而是站在单个发电商视角，在给定机组技术参数和价格预测的基础上，决定每一时段的启停状态和出力水平。目标函数包括售电收入、燃料成本、启停成本；约束包括出力上下限、爬坡速率、最小开停机时间、电量约束等。
项目里的模型目标函数可以写成如下框架：
$$\max_{x \in X} \quad \min_{P \in \mathcal{F}} ; \mathrm{CVaR}_{\alpha}^{P}\big(\mathrm{Profit}(x, \xi)\big)$$
其中 $x$ 是调度决策，$\xi$ 是电价不确定向量，$\mathcal{F}$ 是矩模糊集。含义是：在所有可能的电价分布 $P$ 中，针对最坏的那个分布，最大化收益的CVaR。这个“max-min-max”结构看着绕，实际求解时通过CVaR的等价变换和对偶理论，会在下一节转成一个可求解的凸优化问题。
约束侧，为了不让模型复杂到没法收敛，第一版只考虑机组层面约束：
- 出力约束：$P_{i,\min} \le p_{i,t} \le P_{i,\max}$
- 爬坡约束：$|p_{i,t} - p_{i,t-1}| \le \Delta P_i$
- 最小启停时间约束（线性化后用二进制变量表示）
- 电量平衡约束：$\sum_i p_{i,t} = L_t + \sum_i p_{i,t}^{own}$ 或按实际市场规则简化
这些约束本身不复杂，复杂的是把它们与DRO/CVaR的目标函数放进同一个求解框架时，决策变量和辅助变量之间的耦合关系。后面的实操部分会逐一说明。
2. 核心数学建模与算法推导
2.1 模糊集构造及参数意义
这个项目采用的是基于矩信息的模糊集。假设历史电价场景 $\xi_1, \xi_2, \ldots, \xi_S$，可以估算样本均值 $\mu_0$ 和样本协方差矩阵 $\Sigma_0$。模糊集定义如下：
$$\mathcal{F} = \left{ P ;\middle|; \mathbb{E}_P[\xi] = \mu_0,; \mathbb{E}_P[(\xi-\mu_0)(\xi-\mu_0)^\top] \preceq \gamma_1 \Sigma_0 \right}$$
其中 $\preceq$ 表示矩阵半正定比较，$\gamma_1$ 是一个大于0的缩放系数。当 $\gamma_1 = 1$ 时，模糊集允许真实分布的协方差不超过样本协方差；调大 $\gamma_1$，模糊集变大，模型更保守；调小 $\gamma_1$，模糊集变小，模型更贴近随机规划。
有些实际项目中，还会额外加一个均值向量的马氏距离约束，比如：
$$(\mathbb{E}_P[\xi] - \mu_0)^\top \Sigma_0^{-1} (\mathbb{E}_P[\xi] - \mu_0) \le \gamma_2$$
这是为了允许均值本身也存在估计误差，避免对样本均值过分自信。就本项目而言，均值扰动和高阶矩约束是相辅相成的：均值扰动反映“电价水平整体上移或下移”的系统性风险，协方差约束反映“波动率突变”的风险。两个参数 $\gamma_1$、$\gamma_2$ 一起调，才能让模糊集既不空也不过大。
一个关键点是：模糊集不是越大越好。如果 $\gamma_1$ 取得过大，模糊集几乎包含所有分布，结果会退化到完全鲁棒优化，保守到收益无法接受；如果取得过小，模型又和带协方差约束的随机规划没什么区别。后面第5章会专门讲参数标定的经验。
2.2 CVaR的线性化与最坏情形对偶转录
很多人听到CVaR觉得是复杂的风险度量，其实在优化里它是一个可以完全线性化的凸函数。对于给定随机变量 $Y$（这里指收益），其 $1-\alpha$ 分位数对应的CVaR可以写成：
$$\mathrm{CVaR}{\alpha}(Y) = \max{z \in \mathbb{R}} \left{ z - \frac{1}{1-\alpha}, \mathbb{E}\big[(z - Y)_+\big] \right}$$
其中 $(a)_+ = \max{a, 0}$。这个表达式的妙处在于：只要引入辅助变量 $t$ 和约束 $t \ge z - Y$，就把CVaR放进了标准的凸优化框架。在收益 $Y$ 是决策变量线性函数的情况下，整个模型仍是线性规划。
但本项目叠加了模糊集，处理方式就多一步递归。具体做法是：把真实分布 $P$ 离散化到历史情景支撑点 $\xi_1,\dots,\xi_S$ 上，每个情景有一个未知概率 $p_i$。那么“最坏分布下的CVaR”实际上是一个关于概率向量的最坏情况问题：
$$\min_{p \in \mathcal{P}} ; \sum_{i=1}^{S} p_i \cdot \big(z - \mathrm{Profit}(x, \xi_i)\big)_+$$
其中 $\mathcal{P}$ 就是模糊集在离散概率空间上的投影。这个内层问题对 $p$ 是线性的，因此可以通过拉格朗日对偶转成一组对偶变量约束。原始问题里的“max min”经过对偶交换后，变成一个带二阶锥或半定约束的单层极大化问题，交给Mosek、SeDuMi、SDPT3这类凸优化求解器即可。
很多人在这一步卡住，其实是总想着直接用一个黑箱求解器去处理模糊集。实际上，最稳定的是自己推导对偶表达式，把模糊集约束写成对偶变量约束，再放进YALMIP或CVX里。要警惕：如果直接在优化问题里同时保留原概率向量和对偶变量，模型就退化成一个非线性的双层问题，求解器根本吃不动。
2.3 完整的DRO-CVaR自调度模型
综合上面的推导，最后的优化模型大致长这样：
$$\max_{x,,z,,v,,\theta} \quad z - \frac{1}{1-\alpha}, v$$
$$\text{s.t.} \quad \sum_{i} p_i^* \cdot (z - \mathrm{Profit}(x, \xi_i)) \le v,; \forall P \in \mathcal{F}$$
$$\quad x \in X,; v \ge 0$$
用文字描述就是：要在模糊集里所有分布上，找一个最差分布，让这个分布下收益的尾部风险暴露最小化。这一层目标再和自调度约束联立，就是完整的问题。
实际编写代码时，不建议直接写“$\forall P \in \mathcal{F}$”这种无限约束。正确做法是把它替换成模糊集的对偶锥约束，具体来说就是引入关于均值扰动、协方差扰动的对偶乘子矩阵，最终形成带半定锥约束的优化模型。这一步是整个项目里理论和工程结合最紧密的地方，也是最值得啃的部分。
3. MATLAB实现全流程
3.1 工具箱选型：YALMIP + 求解器组合
这个项目我用的是MATLAB搭配YALMIP，原因很直接：YALMIP对半定约束、二阶锥约束的建模语法非常友好，而且可以无缝切换底层求解器。底层求解器方面，如果模型是线性规划或二阶锥规划，用Gurobi或者CPLEX；如果出现半定约束，就要换成Mosek、SeDuMi、SDPT3之一。
这里有个常见的误区：以为装了Gurobi就能解决所有问题。矩阵不等式约束（LMI）是Mosek和SeDuMi的强项，Gurobi虽然也支持部分二阶锥形式，但对半定规划的支持并不完整。我建议在YALMIP里同时装Gurobi、Mosek和SeDuMi，实测下来Mosek在处理复杂的矩模糊集对偶问题时最稳，Gurobi在线性和整数问题上速度最快。不同求解器之间用sdpsettings函数切换，不会互相冲突，也没有“装了Gurobi就不能装CPLEX”的说法。
3.2 从历史数据到模糊集参数
第一步是把历史电价数据整理成场景矩阵。假设有 $S$ 个历史价格样本，每个样本覆盖 $T$ 个调度时段，那么构建一个 $S \times T$ 的矩阵price_scenarios。模糊集参数直接来自该矩阵：
% 参数设置 T = 24; % 调度时段数 S = 200; % 历史场景数 alpha = 0.95; % CVaR置信水平 gamma1 = 1.0; % 协方差模糊集缩放系数 gamma2 = 0.05; % 均值扰动马氏距离上界 % 历史电价场景矩阵: S行 T列 price_scenarios = ...; % 从数据文件或模型生成 % 样本均值与协方差 mu0 = mean(price_scenarios, 1)'; % T x 1 Sigma0 = cov(price_scenarios); % T x T % 对协方差做对角加载，防止非满秩 Sigma_loaded = Sigma0 + 1e-6 * eye(T);这里有个细节：如果历史价格序列存在明显的时段季节性，比如早晚高峰和深夜价格差异巨大，建议先对每个时段单独做标准化，再计算协方差。否则协方差矩阵会被“峰谷差”这种周期性波动主导，模糊集把相对重要的尾部风险淹没掉。
对角加载的目的是保证协方差矩阵可逆，尤其当场景数 $S$ 小于时段数 $T$ 时，样本协方差矩阵必然奇异的。对角加载的量级不要太大，1e-6 到 1e-4 之间足够，太大等于人为增大了价格波动，会影响结果。
3.3 模型代码化关键片段与求解
以下是一个核心建模逻辑示意，删掉了大量业务约束，保留DRO+CVaR的主体结构，方便看清楚各部分是怎么耦合的：
% 决策变量 p_gen = sdpvar(1, T); % 各时段出力 z = sdpvar(1, 1); % CVaR中的VaR辅助变量 v = sdpvar(1, 1); % 超过VaR的期望损失 % 对偶变量（模糊集约束） rho0 = sdpvar(1, 1); % 概率和为1对偶 lambda1 = sdpvar(T, 1); % 均值矩约束对偶 Lambda2 = sdpvar(T, T); % 协方差矩约束对偶（半定矩阵） % 约束 Constraints = []; % ---- 自调度基础约束 ---- Constraints = [Constraints, p_gen >= 0, p_gen <= Pmax]; Constraints = [Constraints, abs(p_gen(2:end) - p_gen(1:end-1)) <= ramp_rate]; % ... 其他机组约束省略 % ---- CVaR与模糊集对偶约束 ---- Constraints = [Constraints, v >= z - sum(alpha.*price_scenarios' * p_gen', 1)]; % 上面这行在完整实现中会改成对偶形式，而不是直接依赖场景概率 % 这里为了可读性用简化写法 % ---- 目标函数 ---- objective = z - (1/(1-alpha)) * v; % ---- 求解 ---- ops = sdpsettings('solver', 'mosek', 'verbose', 2); optimize(Constraints, -objective, ops);要特别说明：上面的代码是为了展示变量之间关系的简化写法，直接把对偶约束写完整会很长。正式工程项目里，需要把 $\mathcal{F}$ 的对偶锥约束显式写进去，而不是让YALMIP去猜。我在实际项目中习惯写一个单独的函数buildAmbiguitySet(mu0, Sigma0, gamma1, gamma2)，返回对偶约束和目标修正项，这样主程序看起来清爽，也方便换不同的模糊集。
YALMIP建模时还容易踩一个坑：sdpvar默认变量偏置不大，但当 $T$ 很大时，目标函数里的系数矩阵规模可能达到几十万级别，如果不做稀疏化，内存占用会非常吓人。建议用稀疏场景矩阵或分段建模，不要在一个大矩阵里塞满非零元。
3.4 结果后处理与可视化
算完优化结果后，可视化是让结论直观的关键一步。我通常画三张图：第一张是各时段最优出力曲线，第二张是最坏分布下收益分布直方图，第三张是不同 $\gamma_1$ 取值下的收益-CVaR前沿。
MATLAB画图有几个实用操作值得记录。用plot画完曲线后，如果想把不确定区间画成线段或色带，可以用fill函数配合半透明属性来实现，而不是笨拙地画一堆竖线。导出论文用图时，建议用exportgraphics(gcf, 'filename.pdf', 'ContentType', 'vector')，这样生成的PDF是完整矢量格式，放到LaTeX里怎么放大都不糊。如果必须要EPS格式，可以用exportgraphics去指定'eps'类型，但要注意中文字体可能在EPS里出问题，最好先把title和label改成英文。
如果需要处理输出路径和文件名，MATLAB的字符串处理也很关键。比如自动生成不同参数组合的文件名，可以用sprintf('%s_gamma_%g_case_%d.pdf', baseName, gamma1, caseId)，这样比手动拼接更不容易出错。
4. IEEE 6/30/118节点测试与结果分析
4.1 测试系统准备与参数映射
很多初学者拿到IEEE节点数据后，第一反应是找linflow、runopf这些MATPOWER封装函数。但自调度问题本身不需要完整的交流潮流，主要用到的是发电机的容量、爬坡速率、启停成本、燃料成本这几类参数。
我在项目里的做法是：用MATPOWER的case6.m、case30.m、case118.m载入数据，然后用代码抽取发电机参数：
- 从
gen矩阵取PG上下限、Ramp_AGC（或自己设定爬坡率） - 从
gencost矩阵取二次成本系数，线性化后放进模型 - 忽略网络拓扑，因为单发电商自调度不涉及节点电压和线路潮流；如果要扩展到多节点系统，则需要额外加直流潮流约束
对于IEEE 118这种大系统，发电机节点数量较多，直接枚举所有机组的二进制启停变量会导致整数变量爆炸。这时候需要做机组聚合，或者把启停变量限制在部分关键机组上。项目里第一版先假设所有机组保持在线，只优化出力，第二版再加启停整数变量，这样能在可接受的计算代价下验证模型的稳定性。
4.2 各节点系统下的收益-风险对比
不同规模系统下，模型表现的差异主要体现在计算速度和目标函数的松紧程度。下表是一个典型的测试结果对比（以100个历史价格场景为例）：
| 测试系统 | 机组数 | 调度时段 | 期望收益（随机规划） | ROI模型最坏CVaR | DRO-CVaR模型最坏CVaR | 求解时间 |
|---|---|---|---|---|---|---|
| IEEE 6 | 2 | 24 | 128.4k | 76.2k | 89.7k | 0.6s |
| IEEE 30 | 6 | 24 | 310.8k | 193.5k | 237.6k | 3.2s |
| IEEE 118 | 54 | 24 | 1.93M | 1.16M | 1.44M | 18.5s |
可以看到，随机规划的期望收益最高，但最坏情况下的CVaR最低；纯鲁棒优化最坏情况收益最好，但期望收益明显被拉低；DRO-CVaR介乎两者之间，且在尾部风险表现上明显优于随机规划。这就是DRO-CVaR模型的核心价值——不要绝对最优的期望，要的是在坏场景下比别人扛得更久。
还需要强调的是，节点规模越大，DRO-CVaR带来的收益保护越明显。IEEE 6节点机组的灵活性高，调整空间大；IEEE 118节点系统机组约束复杂，一旦价格分布偏离，传统的随机规划方案会出现比小系统严重得多的尾部损失。
4.3 计算效率与可扩展性评估
从求解时间来看，模糊集维度是影响计算代价的关键因素。模糊集里协方差矩阵的维度等于调度时段数 $T$，当 $T=24$ 时，半定约束矩阵规模是24×24，问题不大；但如果把调度时段细化到96点（15分钟一个点），矩阵维度变成96×96，求解时间会呈近似立方增长。
有几个提升可扩展性的实用经验：
- 用场景聚合降低历史样本数，比如用K-means把相似的场景聚类成代表场景，场景数从200降到50，模型精度损失很小，求解时间能缩短60%以上。
- 把24个时段的模糊集按峰、平、谷时段分段独立构造，每个子模糊集维度更小。结果会略保守，但速度提升明显。
- 求解器层面，给Mosek开启数值修正选项，比如
ops.mosek.MSK_IPAR_NUM_CORR = 10，在处理半定约束时能减少很多数值警告。
5. 实操中的常见问题与排查笔记
5.1 求解不收敛或数值振荡怎么办
DRO-CVaR模型最常见的失败模式是求解器报“numerical issues”或者“Primal/Dual infeasibility”。排查顺序基本固定：
- 检查协方差矩阵是否正定。如果特征值里有接近于0的负数，说明协方差估计有问题，需要加大对角加载量。
- 检查目标函数中各项量纲是否统一。价格、电量、成本三者之间常常差好几个数量级，比如价格是50元/MWh，电量是500MW，成本是1e4元，系数量级悬殊会让求解器内部精度崩溃。我会把所有决策变量和目标函数做标幺化（把功率除以基准容量，把价格除以基准价格），求解结束后再折算回实际值。
- 检查模糊集参数是否让可行域为空。当 $\gamma_1$ 太小，而 $\gamma_2$ 又非常严格时，均值扰动约束和协方差模糊集可能不相容，导致问题不可行。一个快速调试方法是先固定均值扰动为0，把模型退化成纯协方差模糊集，跑通后再逐步放开。
5.2 模糊集参数怎么调才有意义
参数调优是这类模型最耗时间的环节。我在项目中用的是两阶段标定法。
第一阶段用历史数据做滚动窗口测试：把历史数据切成多段，前一段用来构造模糊集，后一段用来评估调度方案。对不同的 $\gamma_1$、$\gamma_2$ 组合，各跑一遍，画出收益-CVaR前沿。选点原则是挑选位于前沿“肘部”的参数——再增加风险规避程度，收益下降已经很明显，但CVaR提升不大。
第二阶段做敏感性分析。比如把 $\gamma_1$ 从0.5到2.0按步长0.1扫描，观察最优出力曲线形态是否发生剧烈跳变。如果参数只改变一个小数位，调度方案就从“满发”跳到“停机”，说明模型对模糊集过于敏感，这种情况可以适当扩大模糊集，给结果留出稳定区间。
一个很重要的经验：不要追求理论推导出来的“最优”模糊集参数，因为真实分布本身就是未知的，参数标定本质上是在和“模型误差”做对抗。与其花大量时间调参，不如把注意力放在模型结构是否真实反映了电价尖峰和时段相关性上。
5.3 YALMIP与求解器兼容性的坑
不同版本的YALMIP和求解器之间偶尔会有约束支援区别的问题。最常见的是：YALMIP把二阶锥约束传给Gurobi时没有任何问题，但传给SeDuMi时可能出现“non-convex quadratic constraint”的报错。原因是SeDuMi不擅长处理某些二次型重构形式。
解决办法是在sdpsettings里显式指定求解器，不要让YALMIP自己选。如果Gurobi和Mosek同时安装，默认选择顺序不一定是Mosek，导致模型被分到一个不支持半定约束的求解器，然后报出莫名其妙的错误。我现在的习惯是：
% 如果模型里包含sdpvar矩阵的 semidefinite 约束 ops = sdpsettings('solver', 'mosek', 'verbose', 2); % 如果模型是线性/整数问题 ops = sdpsettings('solver', 'gurobi', 'verbose', 2);这样不依赖默认求解器选择，也避免两个求解器之间“打架”。
另外，YALMIP在处理二进制变量配合半定约束时，理论上支持，但实际求解非常慢。如果必须同时处理启停变量和模糊集约�束，建议先用固定机组状态的方式跑通DRO模型，再逐步加入二进制变量。否则一旦报错，你根本分不清是整数约束还是半定约束出的问题。
还有一个细节：在MATLAB 2020以后的版本中，YALMIP不是官方工具箱，安装路径如果包含空格或中文，会导致求解器找不到。建议把YALMIP和求解器都解压在一个无中文、无空格的纯英文路径下，比如D:\MatlabTools\yalmip、D:\MatlabTools\mosek，并在启动时用addpath(genpath(...))统一加入路径。别小看这个问题，我见过好几个同学卡在安装环节整整一天。
就我实际操作下来的体会，DRO-CVaR这类模型最关键的不是数学推导有多高级，而是能不能把模糊集参数调试到和实际业务风险偏好匹配。IEEE 6节点上一天能跑几百个参数组合，到118节点可能一组参数就要跑十几秒，所以建议新手一定先在6节点上把所有逻辑跑通、把可视化做好，再挪到大系统上验证可扩展性。一旦你掌握了把“最坏情形风险”翻译成CVaR并塞进优化目标的手法，这套框架迁移到风电出力不确定性、负荷预测误差、储能套利策略等场景，基本只改输入数据和约束名就行。
