如果你做过含风电的机组组合一定遇见过这种尴尬凌晨预测风大实际却一片寂静火电机组按计划启动后只能压着最低技术出力运行煤白白烧着反过来预测没风实际却狂风大作机组来不及降下来只能眼睁睁看着风机限功率。问题不在预测算法而在于模型把“预测值”当成了“真实值”。这几年业内逐渐把目光从随机规划和鲁棒优化转向了分布鲁棒优化DRO而“线性准则”就是让DRO从理论走向可计算的关键一步。这篇文章我以机组组合为背景把“基于线性准则的分布鲁棒优化”这条路线彻底拆开从数学原理到Matlab实现完整过一遍顺便说说那些论文里不会写的坑。1. 风电不确定性对机组组合的影响先看一个两机系统1.1 确定性机组组合为什么在风电接入后失效传统机组组合模型里风电出力通常被当作确定性参数即给定一个预测值然后在这个预测值下优化启停和出力。这样做在风电占比很低时没问题因为预测误差相对于系统负荷很小机组爬坡和备用可以自然吸收。但当风电渗透率超过20%甚至30%时预测误差的量级可能达到数百兆瓦相当于好几台大型火电机组的容量。此时如果仍然按确定性模型调度实际运行中可能出现两种情况风电实际出力低于预测系统功率不足需要紧急启动快速机组或切负荷代价极高。风电实际出力高于预测火电下调能力不足不得不弃风造成清洁能源浪费。两种情况的本质是机组组合的启停决策在信息尚未揭晓之前就必须做出而确定性模型没有为“预测错误”预留任何决策余地。1.2 一个两机算例直观展示风险假设系统只有一台燃煤机组和一座风电场负荷固定为300 MW。煤电机组最小出力100 MW最大出力500 MW风电预测出力为150 MW。确定性模型给出方案煤机出力150 MW风电出力150 MW。但如果实际风电只有30 MW煤机虽然可以向上爬坡到270 MW但受限于每小时60 MW的爬坡速率在15分钟调度间隔内只能增加15 MW系统依然缺电需要切负荷。这个例子说明机组组合关心的是“如果风电不像预测那样系统能不能应对”。于是就有了不确定性下的机组组合问题在第一阶段决定启停第二阶段根据风电实际值调整出力。1.3 三种主流建模范式对比处理不确定性的经典方法有三类方法假设优势劣势随机规划已知概率分布经济性最优可处理风险分布估计误差会影响结果鲁棒优化不确定量落在区间内绝对安全不用分布过于保守成本过高分布鲁棒优化真实分布属于一个模糊集兼具鲁棒性与经济性模型复杂度高求解难度大随机规划需要假设风电出力服从某个精确分布但实际中这个分布本身是未知的我们只有历史数据。鲁棒优化则假设风电出力在某个区间内任意取值但它不考虑各个取值发生的可能性导致方案过于保守。分布鲁棒优化的思路是我们不需要知道精确分布只需要知道分布应该满足的一些统计特征比如均值、协方差然后优化最坏情况下的期望成本。这样既利用了分布信息又对“分布估计错误”免疫。2. 分布鲁棒优化的核心模糊集与线性决策规则2.1 模糊集如何描述“分布不确定”分布鲁棒优化里不确定量 ξ 代表风电出力向量真实概率分布 P 未知但知道它属于一个模糊集Ambiguity Set D。最常见的模糊集是矩模糊集D { P : E_P[ξ] μ, E_P[(ξ-μ)(ξ-μ)^T] Σ }也就是所有均值为 μ、协方差为 Σ 的分布构成的集合。μ 和 Σ 可以从历史数据估计。更精细的模糊集还可以加入支撑集约束比如 ξ 必须落在某个区间内或概率距离约束如 Wasserstein 距离。决策变量分为两阶段第一阶段变量 x机组启停状态二进制变量和基准出力连续变量必须在风电不确定性实现前确定。第二阶段变量 y(ξ)根据实际风电出力所做的调整量比如机组出力增量、切负荷量、弃风量是 ξ 的函数。目标函数是min_x { c^T x sup_{P∈D} E_P[ Q(x, ξ) ] }其中 Q(x, ξ) 是给定 x 和 ξ 后的最小调整成本。如果 Q 是一个凸优化问题那么整体就是一个两阶段分布鲁棒优化问题。2.2 线性决策规则把函数优化变成参数优化第二阶段决策 y(ξ) 一般是一个任意函数这导致问题在函数空间中优化无法直接求解。线性决策规则Linear Decision Rule, LDR做了一个关键近似假设 y(ξ) 是 ξ 的仿射函数即y(ξ) y0 Y ξ其中 y0 是常数列向量Y 是系数矩阵。这个近似在工程上非常常用因为很多最优策略在不确定性范围不太大时已经接近线性。使用 LDR 后第二阶段问题被转化为对有限个变量y0 和 Y 的元素的优化从而把无穷维问题降为有限维。2.3 从两阶段DRO到可计算的单阶段问题将 LDR 代入第二阶段目标函数和约束后目标函数中的 sup_{P∈D} E_P[·] 可以解析地写成关于均值、协方差的线性或二次函数。约束条件如果涉及 P则需要保证对所有 P∈D 都成立。对于矩模糊集这种“对任意分布成立”的约束可以转化为半定规划SDP约束使用 YALMIP MOSEK 或 CVX SDPT3 即可求解。在实际机组组合中第二阶段约束通常包括功率平衡、机组出力上下限、爬坡约束等。这些约束如果不确定量 ξ 是线性的那么“对所有 P∈D 都成立”往往等价于一个线性矩阵不等式LMI。这正是标题中“线性准则”的含义通过线性化处理将复杂的分布鲁棒机组组合问题变成一个可计算的凸优化问题。3. 机组组合的DRO-LDR完整建模过程3.1 问题设置与符号定义以一个简化但完整的机组组合为例。设调度时段数为 T发电机组集合为 G。风电场出力随机变量为 ξ ∈ R^T每个时段一个值其均值 μ 和协方差 Σ 由历史预测误差数据得到。备用容量需求来自系统安全规程。决策变量分为两部分第一阶段变量u_{g,t} ∈ {0,1}机组 g 在时段 t 的启停状态p_{g,t}^0机组 g 在时段 t 的基准出力第二阶段变量通过 LDR 近似Δp_{g,t}(ξ) a_{g,t} b_{g,t}^T (ξ - μ)机组出力的调整量满足仿射结构r_{g,t}(ξ) c_{g,t} d_{g,t}^T (ξ - μ)弃风量或切负荷量视场景而定为简化我们假设存在一个系统级调整变量 Δp_t(ξ) 表示在时段 t 的系统功率不平衡调整量或者更精细地建模每台机组。3.2 目标函数基准成本 最坏情况期望调整成本目标函数包含三部分基准场景下的燃料成本和启停成本 Σ_{g,t} [ a_g p_{g,t}^0 b_g u_{g,t} S_g (1-u_{g,t-1}) u_{g,t} ] 其中 a_g、b_g 是燃料成本系数S_g 是启动成本。第二阶段调整成本的最坏情况期望 sup_{P∈D} E_P[ Σ_{g,t} c_g Δp_{g,t}^ c_g^r r_{g,t} ] 这里 Δp^ 表示向上调整量r 表示弃风或切负荷惩罚。因为 LDR 假设 Δp 是 ξ 的仿射函数期望中只涉及 ξ 的均值和协方差因此 sup 可以解析计算。实际上由于 D 是矩模糊集且目标函数关于 ξ 二次最坏情况下的期望等于在均值处取值加上一个与协方差和系数矩阵相关的项。3.3 约束条件第一阶段必须满足基准功率平衡Σ_g p_{g,t}^0 μ_t D_t其中 D_t 是负荷。除基准平衡外更关键的是对任意风电实现 ξ系统必须满足调整后的功率平衡Σ_g [ p_{g,t}^0 Δp_{g,t}(ξ) ] ξ_t r_{g,t}(ξ) D_t由于 Δp 是仿射函数这个等式对所有 ξ 成立的条件是常数项和系数项分别相等。于是可以写成两组线性约束。机组出力上下限约束对每个可能的 ξp_g^min u_{g,t} ≤ p_{g,t}^0 Δp_{g,t}(ξ) ≤ p_g^max u_{g,t}这个约束要求对模糊集内所有分布都成立等价于求一个线性函数在模糊集支撑上的最大/最小值。如果模糊集还包含支撑信息比如 ξ ∈ [ξ_min, ξ_max]则该约束可转化为对支撑端点的检查如果只有矩约束则需要使用鲁棒约束转化技巧最后变成半定约束。爬坡约束考虑时段间调整量p_{g,t}^0 Δp_{g,t}(ξ) - p_{g,t-1}^0 - Δp_{g,t-1}(ξ) ≤ R_g^up同样需要鲁棒化。注意爬坡约束是相邻时段之差这会给 LDR 系数矩阵带来耦合增加求解复杂度。最小启停时间约束只涉及第一阶段变量 u保持不变Σ_{τmax(1,t-T_on1)}^{t} u_{g,τ} ≥ T_on (u_{g,t} - u_{g,t-1})以及对应的停机约束。备用约束在不确定风电发生后系统需要足够的向上备用和向下备用。比如向上备用约束可以写成Σ_g (p_g^max - p_{g,t}^0 - Δp_{g,t}(ξ)) r_t^down_reserve ≥ Reserve_t这是描述当风电低于期望时系统还能提供多少额外的功率。这类约束同样是仿射函数的不等式鲁棒约束。3.4 鲁棒约束的转化方法与最终模型形式对形如 w^T ξ v ≤ 0 的约束其中 w 和 v 是待定决策变量的仿射函数要求它对所有 P∈D 都成立的充分条件是 w^T μ v κ ||Σ^{1/2} w||_2 ≤ 0其中 κ 是与置信水平相关的常数通常取 1.645 或 3取决于希望覆盖概率。这个转化来自于统计中的矩不等式——如果只知道均值和协方差那么对任意分布随机变量落在某个椭球内的概率有下界。通过选择 κ可以控制保守程度。因此最终模型是一个混合整数半定规划MISDP。其中 0-1 变量来自启停状态SDP 约束来自鲁棒化后的随机约束。Matlab 中推荐使用 YALMIP 作为建模层求解器用 MOSEK商业或 SDPT3开源。整数部分则由 YALMIP 内置分支定界处理或者配合外部求解器如 Gurobi 处理外层的整数需要将 SDP 约束作为回调。4. Matlab代码实现从数据构造到求解器配置4.1 代码总体结构我建议把代码分成四个模块数据生成模块生成负荷、机组参数、风电预测误差样本。模糊集构造模块从数据估计均值、协方差并定义模糊集。模型建模模块用 YALMIP 定义决策变量、目标函数和约束。求解与后处理模块调用求解器输出启停计划和调度结果。这种分层结构便于调参和扩展。下面给出核心代码思路。4.2 风电场景生成与模糊集参数估计第一个关键的参数是风电预测误差的均值和协方差。假设历史数据中有多组“预测值 vs 实际值”的样本误差向量 e 实际值 - 预测值。我们假设预测值已知因此 ξ 可以直接用实际值建模或者用误差建模。为方便我们令 ξ 为各时段的风电出力误差向量。代码示例% 假设 errors 是一个 N x T 矩阵N 为历史样本数T 为时段数 mu mean(errors, 1); % T x 1 均值向量 Sigma cov(errors); % T x T 协方差矩阵注意如果 T 比较大而 N 较小协方差矩阵可能病态。我通常会在协方差矩阵上加一个小的正则项Sigma Sigma 1e-4 * eye(T);避免数值问题。模糊集就定义为所有均值为 mu、协方差为 Sigma 的分布支撑集取 [err_min, err_max]每个时段有上下界。支撑集的信息也很重要因为它能显著降低保守度。4.3 YALMIP建模在 YALMIP 中定义二进制变量和连续变量u binvar(G, T, full); % 启停状态 p0 sdpvar(G, T, full); % 基准出力 % LDR 系数 A_adj sdpvar(G, T, full); % 常数项 Y_adj sdpvar(G, T, T); % 线性系数Y_adj(g,t,j) % 第二阶段变量调整量 A_adj(g,t) sum_j Y_adj(g,t,j)*(xi(j)-mu(j))这里的 Y_adj 是一个三维变量实际上 YALMIP 不支持三维 sdpvar 直接做矩阵运算需要用元胞数组或者把 j 维度拆开。更高效的方式是用二维矩阵Y_adj {t}代表每个时段的系数。不过如果 T24直接按二维定义Y_adj sdpvar(G, T, T)在 YALMIP 中会解析为三维但 YALMIP 允许sdpvar(g,t,k)只是操作时要小心。我有时会将三维变量变形为二维块矩阵。为降低复杂度可以省略 Y_adj 中大量交叉项只保留对角项即每个时段的调整量只依赖当前时段的 ξ不考虑时段间相关性。这样会损失一些最优性但求解速度提升巨大而且在时段间相关性较弱的场景下效果很好。实际工程中我一般先用对角结构跑通再视情况加入交叉项。4.4 目标函数和约束的YALMIP表达目标函数可以直接写% 第一阶段成本 fuel_cost sum(sum(repmat(a,1,T) .* p0)) sum(sum(repmat(b,1,T) .* u)); start_cost sum(sum(repmat(S,1,T) .* max(0, diff(u,1,2)))); % 第二阶段成本的最坏情况期望简化版 % 假设向上调整成本系数为 c_up向下调整成本系数为 c_down % 利用二次函数在矩模糊集下的期望公式 % sup E_P[ (A sum_j Y_j (xi_j - mu_j))^2 ] A^2 sum_j Y_j^2 * Sigma(j,j) 2 * ... % 具体需要展开这里略去细节只是示意。实际中如果第二阶段成本是线性的那么 sup E_P[...] 等于在均值处的值加上一个与协方差和 Y 相关的鲁棒项。比较复杂的是约束条件。对于功率平衡约束要求对任何 ξ ∈ D 成立。写成% 对每个时段 t需要满足 % sum_g (p0(g,t) A_adj(g,t) sum_j Y_adj(g,t,j) * (xi_j - mu_j)) xi_t D_t % 该约束对所有 xi 成立逐系数列出常数项约束sum(p0(:,t)) sum(A_adj(:,t)) mu(t) D_t;线性系数项约束对每个 j满足sum(Y_adj(:,t,j)) (tj) 0;这个 (tj) 表示 ξ_t 的系数为 1。类似地机组出力上下限鲁棒约束% p0(g,t) A_adj(g,t) sum_j Y_adj(g,t,j) * (xi_j - mu_j) p_max * u(g,t) % 对所有 xi 在支撑内成立 % 如果支撑为盒式 [xi_min, xi_max]则可以转化为对端点校验 % 但由于 Y 矩阵较大使用端点校验会枚举 2^T 个端点不可行。 % 因此通常利用矩模糊集转化为 SDP 约束对应约束是p0(g,t) A_adj(g,t) w*(xi-mu) ≤ p_max*u(g,t)对所有 P∈D 成立其中 w 是 Y_adj(g,t,:) 的列向量。利用鲁棒约束转化等价于p0 A mu?? 不需要仔细推导。实际上鲁棒约束w^T ξ v ≤ 0对所有均值为 μ 协方差为 Σ 的分布成立最坏情况是 ξ 在某个方向上的极端值。但如果没有支撑集约束ξ 可能取任意大的值约束不可能对所有分布成立。因此必须加入支撑集或者只在概率意义下约束机会约束。在 DRO 中常见的是将硬约束放松为机会约束即要求概率至少 1-ε。这也是为什么“线性准则”之外还需要引入“风险”概念。为避免过于复杂我们可以假设支撑集是一个椭圆或盒式然后使用线性矩阵不等式处理。在 YALMIP 中我们不需要手动推导所有约束可以直接利用 YALMIP 的鲁棒建模能力。YALMIP 支持使用uncertain和expectation等命令但对于自定义模糊集通常需要结合solvesdp或optimizer与手动 SDP 约束。我个人更倾向于手动写出 SDP 约束因为这样可控性高且更直观。手动构造 SDP 约束示例假设使用 YALMIP 的sdpvar% 对每个 g,t考虑约束 p0(g,t)A_adj(g,t)w*(xi-mu) p_max*u(g,t) % 将其改写为鲁棒线性约束等价于存在标量 s 0 使得 % [s*Sigma w*w w*(p0A-p_max*u w*mu); % w(p0A-p_max*u w*mu) (p0A-p_max*u w*mu)^2 - s] 0 % 这个形式需要根据具体推导调整。这里的推导确实容易出错我在实际编码中会用一个小型例子验证约束是否正确。4.5 求解器配置与性能优化一旦模型构建完毕可以用下面的命令求解options sdpsettings(solver,mosek,verbose,1,debug,1); result optimize(constraints, objective, options);如果是 MISOYALMIP 会自动使用分支定界。MOSEK 支持整数变量 二阶锥/半定约束。但 MOSEK 对整数 SDP 支持相对较弱遇到大规模问题可能非常慢。可行的改进方案是使用 Benders 分解把二进制变量作为主问题LDR 和 SDP 约束作为子问题。对于需要发论文或解决大规模场景的朋友我建议先用小系统测试数学正确性再考虑分解算法。5. 求解中的坑与工程化经验5.1 协方差矩阵估计中的数值病态问题当 T 达到 24 甚至 96 时协方差矩阵是 96x96如果历史样本只有几百条矩阵的秩会不足半定规划求解器会报错或收敛缓慢。我常用的处理手段有三个增加正则项Sigma Sigma 1e-3 * eye(T)。使用样本协方差的特征值截断保留前 K 个主成分重构协方差矩阵。使用收缩估计比如 Ledoit-Wolf 方法直接在 Matlab 中实现。实际测试中Ledoit-Wolf 收缩对模糊集的合理性影响很大。收缩后的协方差矩阵不仅能保证可逆还能避免分布鲁棒模型过度信任样本协方差。5.2 LDR近似的适用边界线性决策规则假设调整量与风电误差是线性关系。当风电误差异常大或机组爬坡约束很强时真正的最优调整策略可能是非线性甚至非凸的比如触发启停调整。LDR不一定能捕捉到这种非线性导致模型给出的可行域比实际小或者在真实运行中调整策略不可行。一个验证 LDR 适用性的经验法则比较“仅基准工况”和“最坏工况”下系统功率差额如果差额不超过机组总容量的30%LDR误差通常可以接受。如果差额很大考虑引入分段线性决策规则把 ξ 的取值范围分成若干段每段用一个仿射函数。虽然分段变量会增加计算量但效果显著。5.3 模糊集保守度的调节矩模糊集使用均值 μ 和协方差 Σ但如果我们对这两个统计量的估计不准确DRO模型会过于保守。实际中可以将 μ 和 Σ 周围再套一个球或其他集合形成“模糊集的模糊集”但那样计算更复杂。更常用的方法是引入一个缩放参数 ρ使模糊集定义变为D {P : E_P[ξ]μ, E_P[(ξ-μ)(ξ-μ)^T] ρ Σ}ρ 越大模糊集越大结果越保守。调节 ρ 可以权衡经济性与鲁棒性。在项目交付时我一般会画一条“ρ-总成本”曲线由决策者根据风险偏好选取 ρ。这个曲线通常呈凹函数形态边际成本会随着 ρ 增加而快速上升因此找一个“拐点”比较有意义。5.4 二阶段变量太多时的降维技巧LDR 系数矩阵的大小直接决定求解速度。假设 G10, T24Y_adj 的维度是 10x24x24 5760 个变量。再加上 0-1 变量MISDP 会非常庞大。我的经验是采取以下降维策略将第二阶段调整变量按机组聚合比如所有火电机组等效为一台聚合机组仅保留启停细节在第二阶段之外。只对风电出力影响较大的关键时段使用完整 LDR其余时段使用简化决策规则例如只依赖本时段 ξ。使用滚动时域每次只求未来 6 小时滚动更新运行计划。这些方法在工程上非常有效但要注意别破坏原问题的结构。6. 扩展从线性规则到更高级的决策规则LDR 解决了“可计算性”但代价是最优性损失。理论研究表明当随机参数维数较低时LDR 在最优值上界的误差通常可以控制在 5% 以内。但对高维、强非线性的场景LDR 可能会产生比较大的性能损失。替代方向有两个分段线性决策规则将 ξ 的支撑集分成有限个区域在每个区域上定义不同的 LDR。分段策略能逼近任意最优策略但代价是引入整型变量且整数变量代表“区域索引”会让问题变成混合整数问题。无分布假设下的仿射决策规则直接假设最优调整策略是仿射函数但不对分布做任何假设仅求最坏情况下的负效用。这种思路本质上变成了一个鲁棒线性规划也可以结合 LDR 求解。未来机组组合还会引入储能、需求响应等更多灵活性资源它们的决策变量包含能量状态这使得多阶段问题更加复杂。LDR 在多阶段问题中仍然有效因为能量状态与出力之间的关系是线性的可以将储能电量状态也写作 ξ 的仿射函数。我在实际项目中当系统规模较大时会优先使用 LDR 做初算得到一套启停计划然后采用蒙特卡洛模拟验证该启停计划在不同分布场景下的表现。如果模拟发现某些场景下存在切负荷或弃风再把对应场景放入模糊集加强约束。这种“先求解再对抗验证”的流程比一开始就保守建模要高效得多。如果你正在复现这篇思路建议先用一个简化算例比如 3 台机组、6 个时段把 DRO-LDR 模型跑通再逐步扩大规模。Matlab 配合 YALMIP 和 MOSEK 是目前比较顺手的组合。代码不要追求一步到位我踩过的最大教训是一上来就构建完整 24 时段、10 台机组的模型结果根本不知道是模型错还是数值错。先小后大把每一段约束单独验证正确再拼装最终你会拿到一套稳定可靠的分布鲁棒机组组合代码。