两阶段分布鲁棒优化机组组合：线性决策规则与Matlab实现

先说结论：这套代码的核心并不复杂，一句话可以讲清楚——把风电出力的随机性用一个分布模糊集装起来，第二阶段的机组调整量写成不确定量的仿射函数，也就是所谓的线性决策规则，然后在一个分布鲁棒优化框架里同时优化一阶段启停和二阶段调整系数。直接用Matlab + YALMIP + Gurobi能跑通，适合电气工程专业做机组组合、微网调度、电力市场方向的研究生，也适合刚接触分布鲁棒优化的工程师拿来当模板改自己的模型。

我更想强调的是，很多人在“分布鲁棒优化”这个名字上先被吓住了。实际上，你只要把期望、方差、支撑集这几个概念理清楚，把两阶段决策变量降维成线性形式，剩下的工作就是让求解器去算。这篇文章会把模型、数学推导、代码实现、踩坑记录全部串起来，尽量按照我自己调试这个项目时的顺序来写。

1. 项目整体设计与思路拆解

1.1 机组组合：一个看起来简单、算起来头疼的调度问题

机组组合问题本质上是决定“哪些机组在哪些时段开着”，然后分配每台机组的出力。传统模型里，负荷和风电出力都是确定数值，目标函数是燃料成本加上启停成本，约束包括功率平衡、机组出力上下限、爬坡约束、最小启停时间、旋转备用等等。这是一个混合整数线性规划或者混合整数二次规划问题，量级不大的时候Gurobi/CPLEX可以秒解。

但一旦把风电出力当成随机变量，问题就完全变味了。风电出力的预测误差会传导到功率平衡约束里，如果继续用确定性模型，可能在实际运行时出现电量不足或者弃风过多的情况。更麻烦的是，启停决策需要提前一天定下来，而实际风电出力只能在当天滚动揭晓，这就形成了一个天然的两阶段决策结构：第一阶段定启停，第二阶段等风电出力知道了再调整各机组出力。

我一直觉得，理解这个两阶段结构是读懂整个代码的钥匙。第一阶段变量里藏着0/1整数变量，第二阶段变量是连续量，它们在时间上是先后关系，在模型里却是同时求出来的。解决这种问题，不能用简单的一次性优化，必须引入不确定性优化框架。

1.2 不确定性建模的三条路：场景法、鲁棒优化、分布鲁棒优化

处理风电不确定性，行业内用得最多的是三条技术路线。

第一条是随机规划，尤其是场景法。思路很简单：用历史数据生成几百上千个风电出力场景，把期望成本最小化。这种方法的优点是直观，但缺点是场景数量直接决定求解规模，而且场景生成的质量严重影响结果，更重要的是，你本质上在假设那些场景的概率分布就是真实分布，这在实际中是很难保证的。

第二条是鲁棒优化。它不关心概率，只关心不确定量在某个集合内取值的所有可能情况，然后让约束在最坏情况下也满足。鲁棒优化的优点是稳健性极强，缺点同样是稳健性极强——可能为了一个几乎不可能出现的极端场景，把调度方案变得非常保守，经济性差。

第三条就是分布鲁棒优化。它走的是中间路线：不假设风电出力服从某个精确分布，但也不像鲁棒优化那样完全无视概率信息。你构建一个包含所有“看起来合理”的概率分布的集合，比如均值落在某个范围、协方差矩阵有界，然后优化在这个模糊集合中最坏情况下的期望成本。这样，结果对分布误差有抵抗力，又不像纯鲁棒那样过度保守。

下面这张表是我个人常用的对比方式：

这个项目选择分布鲁棒优化，主要原因是风电预测误差的真实分布很难准确刻画，但均值、协方差这些统计量相对容易从历史数据里估计出来。用矩信息构造模糊集，既不过度依赖分布假设，又能保留概率信息，电网调度场景里很合适。

1.3 线性准则（线性决策规则）在这里扮演什么角色

分布鲁棒模型建好之后，马上会遇到一个尴尬的问题：第二阶段决策变量是“不确定量的函数”，这个函数可以是任意函数，也就是说它是一个无限维变量，计算机没法直接处理。

这时候就需要用线性准则，或者更准确地叫线性决策规则。它的思想极其朴素：假设第二阶段变量和风电出力之间是线性关系。比如风电场w的出力随机变量是ξ_w，那么机组i在t时段的调整量不再是任意函数，而是写成：

r_i,t(ξ) = r_i,t^0 + Σ_w R_i,t,w * ξ_w

这个假设看起来简单，但作用非常大：它把一个无限维的函数优化问题，转化成有限维的系数优化问题。你只需要优化常数项r^0和线性系数R，优化变量个数从“无穷多”直接变成“有限个”。这个转化代价是有可能牺牲一部分最优性，但实际算下来，在机组组合这种调整量对不确定性响应相对平滑的问题里，线性规则往往已经相当接近最优解。

为了让这个概念更好懂，我经常拿它和“用直线拟合曲线”类比。真实最优的二阶段策略可能是一条弯曲的复杂曲线，线性决策规则就是在这条曲线上画一条直线，直线离曲线在一些点上有偏差，但整体能用。如果你觉得一条直线不够，还可以做分段线性决策规则，那就是用折线去逼近，精度更高，代价是变量更多。这个项目名字里明确写着“线性准则”，所以核心就是用一次仿射函数做近似。

2. 核心建模细节与数学推导

2.1 两阶段分布鲁棒机组组合模型

我把模型写成下面这个抽象形式，实际代码里的约束会更多，但骨架是它：

第一阶段：

min Σ_t Σ_i [ c_i * P_i,t + SU_i,t + SD_i,t ] + sup_{P∈D} E_P[ Q(ξ) ]

这里P_i,t是机组i在时段t的计划出力，SU和SD是启停成本，Q(ξ)是在风电出力ξ实现后的第二阶段最优调整成本，D是分布模糊集。

第二阶段问题Q(ξ)的定义是：

min Σ_t Σ_i d_i * (r_i,t^+) + d_i * (r_i,t^-) + pen * (L_shed,t + W_spill,t)

s.t. 对于所有t，功率平衡方程在加入调整量、切负荷、弃风后仍然成立；各机组调整后的出力不超过运行上限；爬坡约束也包含调整量。

第二阶段变量是r(ξ)和切负荷/弃风量。如果没有线性规则，这个函数优化问题直接无解。用了线性规则之后，r(ξ)就替换成r^0 + R ξ，切负荷、弃风量也写成仿射形式。然后问题就变成一个带期望项和半无限约束的有限维优化。

这里有个重要细节：为什么切负荷和弃风要加惩罚系数pen？因为如果风电不足或者过剩，模型需要有一个“泄压阀”，否则功率平衡约束可能无解。工程上我们会把切负荷惩罚设得非常高，比如几千美元每兆瓦时，弃风惩罚相对低一些，这样模型不会轻易用切负荷来优化成本。

2.2 模糊集怎么构造：矩约束为主

模糊集是分布鲁棒优化最核心的设计自由度。这个项目用的是矩信息模糊集，一般写成：

D = { P ∈ P(Ξ) : E_P[ξ] = μ0， E_P[(ξ - μ0)(ξ - μ0)^T] ⪯ Σ0， P(ξ ∈ Ξ) = 1 }

也就是所有均值固定为μ0、协方差矩阵被Σ0按半定序约束住、并且支撑在集合Ξ内的概率分布。这里μ0和Σ0可以从风电功率预测误差的历史样本中估计。支撑集Ξ通常是盒式约束，例如ξ∈[ξ_min, ξ_max]，对应风电出力的物理上下限。

为什么用均值和协方差约束？因为这两个量估计起来稳定，而且对偶推导后能得到一个半定规划或者二阶锥规划，YALMIP可以直接建模。如果你只有一批历史场景，也可以用经验均值、经验协方差，再乘一个保守系数作为Σ0，这就是实践中常用的做法。

还有一点值得提醒：支撑集不能取得太保守，否则效果会接近纯鲁棒；也不能取得太松，否则模糊集里会混入大量不切实际的分布，导致结果变得很奇怪。我一般把ξ_min和ξ_max取为历史观测的最小值和最大值，再稍微外扩5%左右。

2.3 线性决策规则下的约束转换

采用仿射规则后，几乎所有含ξ的约束都变成“对任意ξ∈Ξ都成立”的半无限约束。以功率平衡为例，假设原来某条约束是：

A(x, r^0) + B(x, R) ξ + C ≥ 0， ∀ξ ∈ Ξ

如果Ξ是一个多面体，可以用鲁棒优化的标准对偶去做转换。比如Ξ是盒式集合，这个半无限约束就等价于对R和A的线性约束转换，核心是利用每个ξ分量的端点值。具体推导时要注意，R的每一项都可能带符号，不能简单地把ξ都换成最大值，而要分别判断系数正负，或者用对偶变量统一处理。

代码实现里，我推荐直接用YALMIP的uncertain变量和robust命令，它会把这类半无限约束自动重写。不过实际使用中YALMIP对非线性的不确定性描述支持有限，对盒式集合表现不错。如果我们自己手推对偶，就要额外小心正负号，这也是后面第五部分要讲的常见bug来源。

2.4 最坏期望成本的对偶处理

目标函数里的sup_{P∈D} E_P[Q(ξ)]是最难啃的骨头。好在矩约束模糊集有非常漂亮的对偶形式。如果Q(ξ)是凸函数，那么在一个只有均值固定、协方差有界的模糊集上求期望的最大值，可以转化为一个有限维的对偶问题。

稍微写一下思路。对于包含均值和协方差约束的模糊集，对偶问题会引入两组变量：对应均值等式的向量v，对应协方差不等式半定约束的矩阵W。最后最坏期望可以写成：

inf_{W⪰0, v} v^T μ0 + tr(W Σ0) + something over support

这个“something”在支撑集是盒式时，会变成一个在端点集合上的极大值项，需要用额外的变量去做线性化。如果Q(ξ)本身是仿射函数，那它会特别简化。正因为我们用了线性决策规则，Q(ξ)对所有ξ其实是一个二次函数（因为r是仿射的，再乘上惩罚成本，本质是分段线性），这会带来一些二次项处理，最终模型会变成混合整数二阶锥规划（MISOCP）。Gurobi 9.0以上可以直接解MISOCP。

我在实际建模时，为了减少调试难度，会把第二阶段目标函数的二次项做一次处理，切成一个个辅助变量和线性约束。这样虽然约束数量变多，但数值稳定性更好，YALMIP和Gurobi配合起来也更快。追求完美精确解的话，也可以保留二次项，让求解器内部处理。

3. Matlab代码实现与实操过程

3.1 运行环境准备：Matlab + YALMIP + Gurobi

这套代码的硬性依赖有三个：Matlab、YALMIP工具箱、Gurobi求解器。如果你想用CPLEX或MOSEK也可以，但下面的讲解以Gurobi为主。

Matlab版本建议R2020a以上，太老的版本对YALMIP某些函数支持不好。YALMIP直接去GitHub下载并添加路径即可，不需要额外安装。Gurobi需要申请一个学术许可证，下载后其实就是配置环境变量。安装完之后，在Matlab里执行：

yalmip('clear') which sdpvar which gurobi

如果都能找到路径，说明环境配置成功。我见过大量同学卡在“找不到gurobi”上，主要原因是安装Gurobi之后没有重启Matlab，或者没有运行Gurobi自带的Matlab接口设置命令：

gurobi_setup

这个小命令会帮你在Matlab的pathdef.m里加入Gurobi的Matlab接口目录，重启后依然有效。

YALMIP和Gurobi都建议只用正版学术授权，学校IP申请很方便，别去碰不明来源的工具包。

3.2 代码整体架构和关键变量定义

以典型的6节点3机系统为例，我习惯把代码分成下面几个文件：

主程序的核心变量不要一股脑全用sdpvar，我建议分三层定义。第一层是确定性变量，包括负荷、机组参数、预测值；第二层是一阶段决策量，包括启停变量U（binvar）和计划出力P（sdpvar）；第三层是二阶段线性决策系数，包括常数项R0和系数矩阵R。

一个常见的坑是变量维度。比如机组有3台，时段24个，风电场1个，那么R0的维度是(3, 24)，R的维度是(3, 24, 1)，但YALMIP的sdpvar不支持三维数组，需要把第二维拍平。我通常写成R = sdpvar(3, 24 * nW, 'full')，然后再reshape。第一次写代码的人很容易在这里被矩阵维度错误折磨半天。

3.3 关键代码片段：YALMIP建模

下面给出一个经过简化的核心片段，主要展示线性决策规则如何落到YALMIP里：

% 系统参数 nG = 3; % 机组数 T = 24; % 时段数 nW = 1; % 风电场数 % 一阶段变量 U = binvar(nG, T, 'full'); % 启停状态 P = sdpvar(nG, T, 'full'); % 计划出力 % 二阶段线性决策规则变量 R0 = sdpvar(nG, T, 'full'); % 调整量常数项 Rxi = sdpvar(nG, T * nW, 'full'); % 调整量对风功率的线性系数 % 风电不确定变量（用于声明） xi = sdpvar(nW, 1); % 第二阶段调整量：r = R0 * ones + Rxi * xi r_adj = R0 + reshape(Rxi, nG, T, nW) * xi; % 这里省略功率平衡约束的展开，核心是： % constraints = [constraints, P + r_adj >= Pmin .* U]; % constraints = [constraints, P + r_adj <= Pmax .* U]; % 目标函数：一阶段成本 + DRO最坏期望成本 objective = sum(sum(C * P)) + sum(sum(SU_cost .* U)) + dro_worst_expect_cost; % 求解 ops = sdpsettings('solver', 'gurobi', 'verbose', 2, 'debug', 1); optimize(constraints, objective, ops);

注意这里有个关键点：YALMIP本身并不直接运行“对任意xi都满足”的约束，你需要用uncertain(xi)和with等命令来声明不确定变量。如果你自己把对偶约束已经写出来了，那就不需要YALMIP的robust命令，直接把转换后的约束加进去，求解更快。

完整的模型里，你还需要把风电预测值写成xi = w_forecast + delta_xi，其中delta_xi才是零均值随机量。这样线性规则中的常数项可以拆成“针对预测值的调整项”和“针对误差的响应项”，物理含义更清晰。

3.4 案例测试：6节点3机系统

我用一个标准6节点系统做了测试。3台火电机组，额定容量分别取200MW、150MW、100MW；负荷曲线取典型冬季日负荷，峰值280MW；风电场额定容量120MW，预测出力采用某风场的历史数据，预测误差标准差约为额定容量的10%。

模糊集构造如下：均值μ0设为预测误差的样本均值（基本接近0），协方差Σ0设成样本协方差矩阵的1.2倍，支撑集设为预测值上下50%范围。Gurobi求解时间在同配置电脑上大约40秒到2分钟，取决于模糊集大小。

结果出来后，我通常绘制三张图：第一张是各机组出力曲线和风电出力曲线堆叠图；第二张是启停状态图；第三张是不同模糊集大小下的总成本对比图。可视化在Matlab里很简单，用stairs画机组出力，用area画堆叠负荷，用bar画启停状态即可。

我还试过把第二步的实测风电数据回放进去，计算实际运行成本，发现分布鲁棒方案比确定性方案的平均成本更低，尤其在风电误差较大的时段，这是因为分布鲁棒方案会预留更多爬坡空间，减少了切负荷惩罚。

4. 与场景法和鲁棒优化的对比分析

4.1 分布鲁棒结果到底比场景法好在哪

很多人会问，场景法也可以考虑风电不确定性，为什么非得用分布鲁棒？我的实测感受是，如果场景数足够多，比如2000个以上，场景法的结果和分布鲁棒结果很接近，但场景法的求解时间会成倍增加。分布鲁棒相当于把场景信息压缩成矩信息，在不需要很多场景的情况下得到稳健解。

下面的表是我在某次测试中记录的对比结果：

当然，这个表不是通用结论，不同系统的差异会很大。但它至少反映出分布鲁棒能在求解时间和保守性之间做一个较好的平衡。

4.2 线性决策规则带来的近似损失

线性决策规则虽然方便，但毕竟限制了解空间的形状。如果第二阶段调整策略是非线性的，比如某些机组存在启停状态与出力联合调整的强耦合，仿射近似可能误差变大。

我自己的经验是，当风电渗透率不高（低于30%）时，线性规则的结果非常接近真实最优；当渗透率超过50%后，误差会明显增加。此时可以考虑两个改进方向：一是把风电不确定量按区间分段，每段用不同的线性规则；二是把启停变量本身也做成不确定量的仿射函数，但这样会引入0/1仿射变量，求解难度陡增，一般不推荐。

所以，线性准则适合作为第一版实现，先把框架跑通，再根据实际精度需求决定是否升级。

5. 常见问题与排查技巧实录

5.1 YALMIP/Gurobi求解时报不可行

这是最多人问的问题。第一个要排查的是功率平衡约束的系数是否写反，尤其是线性决策规则展开后，所有包含ξ的项必须在“对任意ξ成立”的约束下逐项匹配。最简单的方法是先固定住Rxi=0，只优化R0，跑一个纯确定性模型，如果能通过，说明模型骨架没问题。然后把Rxi放开，如果此时不可行，再检查模糊集和支撑集是否冲突。

我建议在调试阶段把敏感约束打印出来：

assign(U, initialGuess); check(constraints)

YALMIP的check函数会显示每个约束的残差，残差为负的约束就是违规约束。

5.2 模糊集参数选择不当导致结果异常

均值μ0如果偏离实际太大，会导致模型过于乐观；协方差Σ0如果太小，模糊集没有覆盖真实分布，结果会不可靠；如果太大，结果又趋于保守。我的调参方式是做一个敏感性扫描：把Σ0从0.5倍样本协方差扫描到2倍样本协方差，观察总成本变化。如果成本上升过快，说明模型对模糊集很敏感，此时需要增加支撑集约束来抑制极端分布进入模糊集。

另外，样本量太少时，直接用样本协方差会低估不确定性，建议用统计方法做一点保守化处理，比如加上一个与样本量成反比的修正项。

5.3 计算时间爆炸怎么办

计算时间主要来自三个方面：整数变量数量、二阶段决策规则系数数量、以及半定约束的规模。如果一个系统规模较大，比如上百台机组，建议先用聚合机组方式降低机组数量，再或者采用时段聚合，把24小时聚合成典型时段集合。

还可以考虑把模糊集简化成一个基于1范数或无穷范数的矩约束，这样最坏期望成本的对偶会变成线性规划而不是半定规划，Gurobi求解速度快很多。虽然精度略降，但工程上完全够用。

5.4 代码运行中的其他小坑

第一个坑是binvar与full参数。在MATLAB较新版本中，binvar(nG, T, 'full')没问题，但如果你直接写binvar(nG, T)，默认可能是稀疏结构，后续参与计算时容易出现维度不匹配。我一般统一写成'full'。

第二个坑是求解器的容差设置。默认的MIP gap如果是1e-4，对于这个模型通常已经够用，但如果你发现目标值在重复运行中有细微波动，可以把sdpsettings('gurobi.MIPGap', 1e-3)设得松一点，换取速度。

第三个坑是Gurobi在Windows和Linux上的接口路径不同。在Linux下，gurobi_setup有时不会自动生效，建议在startup.m里手动添加路径，否则每次重启Matlab都要重新设置。

6. 可以怎么继续扩展

这个项目做到能跑、能出图之后，可以往很多方向延伸。如果你手头有风电历史数据，可以先做一步风功率预测，然后用预测误差的历史样本构造模糊集。这时候可以接一个BP神经网络做误差修正，再用修正后的误差序列重新估计协方差，能进一步降低模糊集的保守性。我自己试过用长短时记忆网络做预测，再接入这个DRO框架，效果也不错，但要注意预测模块和调度模块的时间尺度要匹配。

如果系统里还有储能，可以把储能当成特殊机组，增加SOC（荷电状态）的时序约束，这时第二阶段调整量不只有出力，还有充放电功率，线性决策规则依然适用，只是变量维度增加。如果要考虑多时段耦合的爬坡约束，二阶段线性规则的展开会更复杂，因为ξ不只是一个时刻的向量，而是整个调度周期内的随机过程，系数矩阵会变成块三对角结构。这个时候务必采用稀疏矩阵建模，否则内存容易爆。

另外，线性准则升级为分段线性决策规则也是一个很好的研究方向，实现方式是把支撑集划分成几个子区域，每个区域用一套线性规则。求解时会引入一些整数变量用来选择区域，模型就从MISOCP变成MIQP或者混合整数线性规划，Gurobi都能处理，但调试难度会高不少。

我在实际使用这个代码的过程中最深的体会是：分布鲁棒优化的核心并不在求解器，也不在代码本身，而在于你对“不确定性到底怎么描述”这个问题想得有多清楚。模糊集选得好，方法效果立竿见影；模糊集选得随意，再高级的算法也是空中楼阁。最后再分享一个小技巧——跑完代码后，把风电误差回放进去，算一下实际成本，并且把不同方法的结果放在同一张图上比较，这比任何理论分析都更能说服审稿人或者领导。