1. 先搞清楚风力发电不确定性到底难在哪里1.1 风电出力的三副面孔做电力系统优化的人嘴上天天挂着不确定性但真正落到机组组合Unit CommitmentUC模型里风电不确定性其实有三副面孔理解不到位后面全白搭。第一副面孔是随机性。风电出力不是某个确定值而是服从某种分布的随机变量。今天10点钟风大明天同一时刻可能风小这种波动本质上是随机的。第二副面孔是间歇性。风机启动、切出都有风速阈值风速低了不发风速高了限发出力曲线不是连续的而是带着大量毛刺和断层。第三副面孔是预测误差的时空相关性。风电预测不是只预测一个点的出力而是预测一条出力曲线。误差在时间上相关——前一小时预测偏大后一小时往往也偏大在空间上也相关——同一风电场的风机受同一股风影响误差方向常常一致。这种相关性如果建模时不考虑调度方案就会过于乐观实际运行中容易出现断面越限。在做机组组合时传统做法是把风电预测值当作确定性的净负荷输入然后留一定备用容量。这种做法不是不行但问题在于备用留少了风电预测偏差一大就可能切负荷备用留多了火电机组压出力运行煤耗和碳排放都上去了经济性很差。所以业界和学界这些年一直在探索怎么在建模阶段就把风电的不确定性内化到优化决策里而不是靠事后加备用。分布鲁棒优化Distributionally Robust OptimizationDRO就是这条路线上很有代表性的一个方向。1.2 从随机优化到分布鲁棒优化一次保守度的权衡要理解分布鲁棒优化得先看它的两个亲戚。**随机优化Stochastic Optimization**假设风电出力的概率分布完全已知比如假设它服从某个正态分布或者用一堆历史场景来近似。思路清晰但痛点也很明显真实分布永远不可能精确已知你假设正态实际可能是偏态、多峰的。模型算出来的最优解只有在你假定的分布下是最优的换到真实分布下可能就差很多。这就是所谓的分布错配问题。**传统鲁棒优化Robust Optimization**走向另一个极端只考虑风电出力落在某个不确定集合比如区间盒式集合内然后优化最坏情况下的成本。这样做的好处是解在任何可能出力下都可行但代价就是过于保守。因为鲁棒优化保证的是所有情况都能扛住可实际运行中很多极端场景出现的概率几乎为零为了这些几乎不可能发生的场景去抬高运行成本经济上划不来。分布鲁棒优化正好卡在两者中间我不假设真实分布是某个具体分布而是假设它属于一个模糊集Ambiguity Set——这个集合里包含了所有满足已知统计信息的分布。比如我知道历史数据里风电出力的均值是某个向量、协方差矩阵是某个矩阵那我就把所有均值和协方差都匹配这个统计信息的分布都收进模糊集里。然后我优化的目标是在模糊集里最坏的那个分布下期望运行成本最小。用大白话说就是随机优化赌分布是对的鲁棒优化不怕任何情况分布鲁棒优化则在两者之间找一个最坏情况下的期望最优。提示如果你对数学形式比较敏感分布鲁棒优化的目标函数通常写成 min max E_P[f(x, ξ)] 的形式其中内层 max 是在模糊集中找最坏分布外层 min 是决策变量。这个min-max结构是整个模型的核心。这个思路好是好但有个现实问题两阶段机组组合模型一旦引入模糊集内层求期望、外层做决策双层结构直接暴力求解计算量惊人。这时候就需要线性准则也叫线性决策规则Linear Decision RuleLDR出场了。2. 核心建模线性准则如何把难题变简单2.1 什么是线性决策规则LDR先回顾一下两阶段机组组合的基本结构。第一阶段是日前决策在风电出力还没完全确定时决定哪些机组开机、哪些机组停机、每台机组的基准出力是多少。这些决策用 0-1 变量和部分连续变量表示是看到不确定性之前就要定下来的。第二阶段是实时调整风电出力揭晓后发现实际出力和预测有偏差这时候机组通过上调、下调出力来平衡系统。这个阶段的决策是看到不确定性之后才做的。问题在于第二阶段决策本质上是风电出力 ξ 的函数也就是说不同的风电出力会导致不同的调整策略。理论上这个函数可以是任意形式可以是高度非线性的。但如果我们允许这个函数是任意函数问题就变成了无限维优化根本没法直接求解。线性决策规则的核心思想很朴素把第二阶段的调整策略限制为风电出力 ξ 的仿射函数。数学表达就是y(ξ) y₀ W ξ其中 y₀ 是基准调整量不依赖 ξ 的部分W 是线性系数矩阵反映对风电波动的响应灵敏度。这个限制把无限维的决策空间压缩成了有限维只要确定 y₀ 和 W第二阶段在任何 ξ 下的调整量都能直接算出来。模型从难以求解的泛函优化变成了能上求解器的有限维优化。可能有人会问这样限制之后解还会是最优的吗说实话LDR 一般求不到全局最优它得到的是原问题的上界对最小化问题而言。但它有两个实实在在的好处第一可解性。引入 LDR 之后很多分布鲁棒机组组合模型可以转化为线性规划LP或二阶锥规划SOCP配合 YALMIP 调用 Gurobi 或 CPLEX求解效率非常可观。第二工程合理性。在实际电力系统运行中自动发电控制AGC对机组出力的调整本来就是近似线性的——频率偏差信号进来机组按照有功-频率静态特性增减出力。所以调整量随不确定性线性变化这个假设和实际运行机制是吻合的。注意LDR 近似在不确定性波动范围较小时精度很好。波动越大线性假设带来的近似误差越大。如果风电渗透率特别高、出力波动幅度达到装机容量的 50% 以上建议考虑分段线性决策规则或加辅助随机变量来提升精度这个后面实测部分再说。2.2 模糊集的构造矩信息驱动的分布簇LDR 解决的是第二阶段决策长什么样的问题模糊集解决的是最坏分布在哪里找的问题。学术文献里模糊集的构造方式五花八门但从工程可操作性来看矩模糊集是最常用也最容易落地的一种。它的思路是利用风电预测误差的历史数据统计出均值向量和协方差矩阵然后定义模糊集为所有满足这两个统计约束的分布集合F { P : E_P[ξ] μ₀, E_P[(ξ-μ₀)(ξ-μ₀)ᵀ] Σ₀, P ∈ P(Ξ) }其中 μ₀ 是历史预测误差的样本均值Σ₀ 是样本协方差矩阵P(Ξ) 表示支撑集 Ξ 上的所有概率分布。支撑集 Ξ 一般取成盒式约束也就是每个时间段的预测误差落在 [ξ_min, ξ_max] 区间内。这个区间可以根据历史数据的最大值和最小值来定也可以按 3σ 原则取。构造模糊集时有三个关键细节直接影响解的保守度均值 μ₀ 的估计直接用样本均值没问题但如果历史数据量很少比如只有几十天的数据建议用缩尾均值或加正则项避免对离群点过度敏感。协方差 Σ₀ 的估计样本协方差矩阵在高维情况下比如 24 时段往往病态特征值差异极大。实操中我一般做两个处理一是加一个小的正则项 Σ₀ εI保证正定性二是如果数据量不够用对角协方差矩阵代替全矩阵牺牲部分相关性信息换计算稳定性。支撑集 Ξ 的选取盒式支撑集如果取得太宽最坏分布可能把所有概率质量都堆到盒子的角点上导致结果过于保守。我的经验是支撑集按历史数据的 90% 分位数来截断而不是简单取 min/max这样能排除极端离群点的影响。模糊集确定之后原始的分布鲁棒模型里的内层问题——即在模糊集中寻找最坏分布使期望成本最大——可以通过对偶理论转化为一个有限维凸优化问题。这一步转化是数学上最重的一环但好消息是文献中已经有非常成熟的结论当模糊集由均值、协方差和支撑集定义时对偶问题可以写成半定规划SDP或 SOCP 形式直接用现成求解器处理即可。2.3 完整的两阶段机组组合模型框架把 LDR 和模糊集拼到一起整个模型的骨架就出来了。目标函数分两部分第一阶段成本机组的启停成本 基准出力对应的燃料成本。这部分不依赖 ξ是确定性的。第二阶段成本针对模糊集中最坏分布下的期望调整成本。包括上调备用成本、下调备用成本、弃风惩罚等。由于 LDR 把调整量写成 y₀ Wξ这部分期望成本可以解析地表示为 y₀、W 和模糊集统计量μ₀、Σ₀的显式函数不需要枚举场景。约束条件也分两部分第一阶段约束机组启停逻辑约束、最小启停时间约束、基准出力上下限约束、系统功率平衡约束用风电预测值、线路潮流约束用预测值。第二阶段约束对任意 ξ ∈ Ξ调整后的机组出力仍要在上下限内调整后的系统功率平衡要成立线路潮流要不过载。第二阶段约束是对任意 ξ 都成立的鲁棒约束这个对任意 ξ让问题带上了半无限约束的性质。但好消息是由于调整量 y(ξ) y₀ Wξ 是 ξ 的仿射函数约束可以写成 A ξ b ≤ 0 的形式只要对支撑集 Ξ 的所有顶点检验即可。如果 Ξ 是盒式约束那只需要检验 2^T 个顶点T 是时段数。24 时段就是 2^24 ≈ 1600 万个顶点太多了但好在这些约束的结构是线性且单调的可以转化为对每个时段独立检验最坏情况计算复杂度直接降下来。模型整体写出来是一个混合整数规划MIP——第一阶段有 0-1 变量转换后的第二阶段约束是线性或二阶锥约束。这类模型在 Gurobi 或 CPLEX 上中等规模系统比如 10 台机组、24 时段通常能在几分钟内求到良好可行解。3. Matlab 实现全流程拆解3.1 代码整体架构与文件组织说实话网上关于分布鲁棒优化的代码不少但大多是为了发论文写的可读性差到离谱。自己动手实现的话我建议按下面的文件结构组织跑通一个算例再逐步扩展DRUC/ ├── main.m % 主程序入口 ├── data/ │ ├── gen_data.m % 机组参数定义 │ ├── wind_data.m % 风电出力与误差数据 │ └── load_data.m % 负荷数据 ├── model/ │ ├── build_first_stage.m % 第一阶段约束与目标 │ ├── build_second_stage.m % 第二阶段约束与目标 │ └── build_ambiguity_set.m % 构造模糊集参数 ├── solver/ │ ├── solve_druc.m % 调用求解器求解 │ └── post_process.m % 结果分析与可视化 └── utils/ ├── get_scenarios.m % 场景生成用于对比验证 └── plot_results.m % 绘图工具这个架构的出发点很简单不确定性的数据准备、模型构建、求解、后处理四块完全解耦。后期如果想换模糊集类型比如从矩模糊集换成 Wasserstein 模糊集只需要改build_ambiguity_set.m和build_second_stage.m其他部分不用动。建模工具我强烈推荐YALMIP原因有三个一是语法和数学表达式一一对应调试的时候不容易晕二是内置了对 LDR 的支持yalmip(definerules)可以自动把仿射决策规则展开三是后端求解器切换方便Gurobi、CPLEX、MOSEK、SDPT3 随便换。如果不想用 YALMIP直接用 CVX 或者纯 Gurobi 回调建模也行但代码量至少多一倍而且 LDR 展开容易出错。3.2 关键代码模糊集参数与场景生成构造模糊集的代码不长但每一步都有讲究。function amb build_ambiguity_set(history_err, T, alpha, lambda) % history_err: T x N 矩阵N 天历史预测误差数据 % alpha: 置信水平参数用于支撑集截断 % lambda: 协方差正则化系数 % 1. 样本均值向量 mu0 mean(history_err, 2); % 2. 样本协方差矩阵带正则化 Sigma0 cov(history_err); Sigma0 Sigma0 lambda * eye(T); % 保证正定 % 3. 支撑集按分位数截断 xi_min quantile(history_err, (1-alpha)/2, 2); xi_max quantile(history_err, 1 - (1-alpha)/2, 2); amb.mu0 mu0; amb.Sigma0 Sigma0; amb.xi_min xi_min; amb.xi_max xi_max; amb.type moment; end这里有三个实操注意点cov(history_err)返回的是 T×T 的协方差矩阵这个方向别搞反了。history_err的行是时段、列是天数所以转置后求协方差才能得到时段之间的相关结构。正则化系数lambda我一般取1e-3 * trace(Sigma0) / T这个量纲与协方差矩阵对角元素相近能起到稳定作用又不至于让正则项喧宾夺主。支撑集的分位数截断参数(1-alpha)/2alpha 取 0.1 意味着保留中间 90% 的历史观察到的误差范围排除两头的极端离群点。这个策略比直接取 min/max 稳得多。场景生成是另一码事。虽然 DRO 模型本身不需要场景但验证模型时要和随机优化SO做对比SO 需要场景。场景生成用简单的拉丁超立方采样加 Cholesky 分解即可function scenarios get_scenarios(mu0, Sigma0, N_sce) % 从正态分布采样但截断到支撑集范围内 T length(mu0); L chol(Sigma0, lower); raw randn(T, N_sce); scenarios mu0 L * raw; % 截断到支撑集 scenarios max(scenarios, repmat(amb.xi_min, 1, N_sce)); scenarios min(scenarios, repmat(amb.xi_max, 1, N_sce)); end提示如果历史数据能给出更准确的边缘分布可以用 copula 采样替代纯正态假设。但实测下来只要样本量不是太小正态假设加支撑集截断的效果已经够用不必过度设计。3.3 关键代码LDR 转化后的主问题建模这一步是整个代码的核心。为了说清楚我简化一下场景假设系统只有一台火电机组和一座风电场T 24 时段。决策变量分三层% 第一阶段决策 u binvar(1, T); % 启停状态 p0 sdpvar(1, T); % 基准出力 % 第二阶段 LDR 参数 w0 sdpvar(1, T); % 调整量基准项 W sdpvar(T, T); % 调整量对不确定性的线性响应矩阵这里 W 是 T×T 的矩阵W(i,j) 表示第 i 时段的调整量对第 j 时段风电预测误差的响应。如果认为不同时段之间的误差完全解耦可以把 W 限制为对角矩阵变量数从 T² 降到 T计算速度大幅提升但代价是忽略了误差的时间相关性。我建议先跑对角版本验证程序正确性再放开为全矩阵看相关性的影响。第二阶段出力表达式为xi sdpvar(T, 1); % 不确定性变量预测误差 p_adj w0 W * xi; % LDR 定义的调整量 p_total p0 p_adj; % 实际出力第二阶段期望成本的建模是利用模糊集的统计信息解析展开的。以调整成本为例假设上下调成本系数分别为 c_up、c_down这里简化为只写上调用成本% 期望上调成本 % E[(p_adj)_] 的精确计算需要知道分布但通过 LDR 和矩信息 % 可以化为凸逼近或引入辅助变量。 % 常见的处理引入辅助变量 t满足 t 0, t p_adj % 然后 E[t] E[p_adj max(0, -p_adj)] 需要进一步凸松弛。 % 实际工程中的简化处理期望调整成本 c_up * abs_sum(w0 W*mu0) % 加上协方差修正项 expected_up_cost c_up * (sum(w0) sum(W * amb.mu0) ... sqrt(sum(sum(W.^2 .* amb.Sigma0, 1), 2)));这个式子里的平方根项来自协方差矩阵——如果调整量是 ξ 的线性函数且 ξ 的协方差为 Σ₀则调整量的标准差正比于 sqrt(W Σ₀ Wᵀ)。用均值加一个标准差的启发式来估计期望调整成本在工程上是常用的近似做法。如果论文要求严格应该把期望成本写成对偶后的半定约束形式这个后面在提升精度小节展开。第二阶段鲁棒约束的处理也一样。比如出力上下限约束% 对任意 xi 属于支撑集p_total 都要在 [Pmin, Pmax] 内 % p0 w0 W*xi Pmin对任意 xi ∈ [xi_min, xi_max] % 等价于 p0 w0 W*xi_min Pmin 最坏情况取 xi_min 处 % 以及 p0 w0 W*xi_max Pmax 最坏情况取 xi_max 处 constraints [constraints, ... p0 w0 W * amb.xi_min Pmin, ... p0 w0 W * amb.xi_max Pmax];因为 LDR 是线性的不等式在盒式支撑集上取到极值的点一定在顶点所以只需要检验端点即可。这是整个 LDR 方法能高效求解的核心原因。把目标函数和约束拼起来调用求解器optimize(constraints, objective, sdpsettings(solver, gurobi, ... verbose, 2, mipgap, 0.01));3.4 求解器选择与求解参数调优求解器的选择取决于模型转换后的具体形式如果第二阶段期望成本用的是前述启发式均值 标准差整个模型是混合整数二阶锥规划MISOCPGurobi 或 CPLEX 都能直接解。如果期望成本通过半定规划对偶严格建模模型变成混合整数半定规划MISDP这个 Gurobi 不支持得用 MOSEK 的混合整数锥优化模块或者用 Benders 分解把 SDP 子问题剥离出来。对于中等规模测试系统10 台机组、24 时段我在 Gurobi 11 上的实测表现是对角 W 模型 30 秒内收敛到 1% 间隙全矩阵 W 模型大约 5-10 分钟视机组数量而定。求解参数调整上有两个小技巧第一0-1 变量优先。机组组合问题的难点在启停变量不是连续变量。用sdpsettings(gurobi.MIPFocus, 1)让求解器先集中精力找可行解再提升目标精度比默认参数快不少。第二边界推导。给启停变量加一个基于负荷预测的可行性剪枝约束比如负荷峰谷差超过所有机组最大可调范围时限制最小开机台数。这种启发式裁剪可以显著减少分支定界的节点数。4. 实测中的坑与排查技巧4.1 模糊集参数到底怎么设这是所有人上手 DRO 时问得最多的一个问题。模糊集参数直接决定了模型的保守程度但文献里很少给出工程上的参考值。均值 μ₀这个没太多调参空间就是历史数据的样本均值。但如果历史数据覆盖的季节跨度很大比如风电的冬夏特性差异明显建议按季节分别建模而不是混在一起。我吃过这个亏全年数据混着估计 μ₀结果春季的预测误差特性完全被平均掉了调度方案在春季的适用性很差。协方差 Σ₀除了加正则项之外一个实用的做法是给协方差乘一个缩放系数 κΣ₀ κ · Σ₀κ 1 放大不确定度模型更保守κ 1 缩小不确定度模型更激进。这个 κ 可以直接作为保守度的调节旋钮。实际使用中从 κ 1 起步每次增加 0.2对比解的成本和实际仿真中的失负荷率找到既不切负荷又不浪费备用的临界值。支撑集 Ξ前面提过用 90% 分位数截断。但注意支撑集和协方差是联动的——如果把支撑集压得太窄最坏分布几乎不可能出现极端出力模型就过于乐观如果把协方差调大但支撑集不变模型又可能过保守。我习惯的做法是先固定协方差扫描支撑集分位数90%、95%、99%看目标函数值的变化趋势如果目标变化超过 5%说明支撑集对结果影响过大需要重新审视数据的离群情况。4.2 线性决策规则近似精度不够怎么办LDR 的近似误差在以下三类场景中会明显放大风电出力波动范围大ξ 的支撑集很宽机组调节能力弱出力上下限范围窄不确定性对系统平衡的影响是非线性的比如网络阻塞导致潮流与出力关系非线性。如果实测发现 LDR 解对应的真实成本用 Monte Carlo 场景返还测试和模型预测成本差距超过 15%就说明 LDR 近似太粗糙了。我的处理策略是按优先级逐一尝试**第一优先加辅助随机变量。**把部分控制变量从用 ξ 直接表示改成用 ξ 的一个线性变换表示让 W 矩阵有更多自由度。具体做法是引入辅助变量 z 作为 ξ 的扩展比如 z [ξ; ξ²] 或者 z 是 ξ 的滞后项LDR 变成 y y₀ W z。这种方法在文献中称为增强型 LDR实现简单效果立竿见影。**第二优先分段线性决策规则。**把支撑集分成若干个子区域在每个子区域内各用一条 LDR。这相当于用分段仿射函数逼近真实决策规则精度显著提升但代价是模型规模乘以子区域数求解时间相应增加。**第三优先换成二次决策规则。**直接把 y 写成 ξ 的二次函数。二次决策规则在理论上比 LDR 更一般但目标函数和约束不再是凸的求解难度指数级上升一般做法是引入 McCormick 包络做松弛对松弛误差需要额外验证。从工程项目的角度我强烈建议先走第一优先路线改动最小且 90% 的精度问题都能解决。4.3 求解性能上不去的几个真因如果你的模型不大10 台机组以内但求解器跑几个小时还收不了官大概率问题不在求解器而在建模方式。**原因一W 矩阵太密。**全矩阵 W 虽然有理论上的最优性保证但实际场景相邻时段的误差之间相关程度有限远端的 W(i,j) 往往算出来接近零。可以在建模阶段就给 W 加稀疏模式只允许 |i - j| ≤ 2 的时段之间有交叉响应。这样 W 的变量数从 T² 降到大约 5TMIP 的求解难度指数级下降。**原因二目标函数的绝对值项没做线性化。**带绝对值的最优性条件让 MIP 求解器非常痛苦。手写转换时遇到 |x| 一定要引入辅助变量拆成 x⁺ x⁻ 的形式而不是让求解器自己处理非光滑项。**原因三大 M 约束的 M 取值过大。**机组组合里的形如出力上限 × 启停状态的约束通常会引入大 MM 取 1e6 这种值看着安全但会让 LP 松弛的边界变得极度松散分支定界效率骤降。正确做法是 M Pmax - Pmin 备用容量上限能用物理含义推导出来的尽量推导。4.4 对比实验怎么设计才有说服力做 DRO 项目免不了要和随机优化SO、传统鲁棒优化RO做对比。对比实验设计得好不好直接决定结论的可信度。我建议至少做三个维度的对比**维度一不同不确定性处理方式在相同数据下的表现。**SO、DRO、RO 三种模型用同一套历史数据训练、同一套场景做测试。重点比较三个指标总运行成本、失负荷概率、弃风率。预期结果是 SO 名义成本最低但失负荷概率最高RO 不失负荷但成本最高DRO 居中。**维度二模糊集参数敏感性。**固定其他条件扫描 κ协方差缩放系数从 0.5 到 1.5看 DRO 解的成本和失负荷率怎么变。如果 κ 稍微一变结果就大幅波动说明模型对模糊集选择高度敏感需要在报告中明确说明。**维度三样本外测试。**用训练期之后的一个季度真实风电数据做模拟运行验证 DRO 调度方案的样本外表现。这一步必须有否则你只是在证明模型在训练数据上表现好而 DRO 的初衷恰恰是要应对分布未知的情况。实操中我还有一个小技巧在论文或报告里画成本-稳定性前沿曲线横轴是失负荷概率纵轴是运行成本把 SO、DRO、RO 三个点标在同一条曲线上。这条曲线直观地展示了三种方法在经济性-鲁棒性之间的权衡位置是审稿人和项目验收专家最容易理解的形式。回到代码项目本身我的一个体会是分布鲁棒机组组合的实现难点其实不在数学推导而在模型转换后的每一步都要反复验证——模糊集参数有没有被错误传递、LDR 展开后的维度对不对、鲁棒约束的极值点检验方向对不对这些细节任何一个出错结果都会静默地变得不合理只有从简单到复杂逐步加约束才能有效定位问题。做这类项目慢就是快先把 3 台机 24 时段的模型跑透再扩展到大系统能省下大量排查时间。