电力市场里做日前自调度最难的不是机组组合那套整数变量而是电价到底怎么建模。你拿着历史场景做随机规划第二天来个尖峰价格利润直接被打回原形改用区间鲁棒优化又把最乐观的情况全丢掉调度方案保守到让交易员想骂人。这个项目做的就是另一条路用分布鲁棒优化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之后得到的是“在所有可能的坏分布里最差的那部分场景下我们还能保住多少收益”。从调度员角度来说这个指标比纯期望收益直观得多。说直白点随机规划是“平均心态”鲁棒优化是“被害妄想症”DROCVaR是“保留一定的风险警觉但不至于脱敏”——它把分布的不确定性和收益的尾部风险同时纳入了优化目标。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_{i1}^{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 模型代码化关键片段与求解以下是一个核心建模逻辑示意删掉了大量业务约束保留DROCVaR的主体结构方便看清楚各部分是怎么耦合的% 决策变量 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模型最坏CVaRDRO-CVaR模型最坏CVaR求解时间IEEE 6224128.4k76.2k89.7k0.6sIEEE 30624310.8k193.5k237.6k3.2sIEEE 11854241.93M1.16M1.44M18.5s可以看到随机规划的期望收益最高但最坏情况下的CVaR最低纯鲁棒优化最坏情况收益最好但期望收益明显被拉低DRO-CVaR介乎两者之间且在尾部风险表现上明显优于随机规划。这就是DRO-CVaR模型的核心价值——不要绝对最优的期望要的是在坏场景下比别人扛得更久。还需要强调的是节点规模越大DRO-CVaR带来的收益保护越明显。IEEE 6节点机组的灵活性高调整空间大IEEE 118节点系统机组约束复杂一旦价格分布偏离传统的随机规划方案会出现比小系统严重得多的尾部损失。4.3 计算效率与可扩展性评估从求解时间来看模糊集维度是影响计算代价的关键因素。模糊集里协方差矩阵的维度等于调度时段数 $T$当 $T24$ 时半定约束矩阵规模是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并塞进优化目标的手法这套框架迁移到风电出力不确定性、负荷预测误差、储能套利策略等场景基本只改输入数据和约束名就行。