1. 从单场景最优到最坏情况可控单元承诺为什么必须走向分布鲁棒电力系统的机组组合Unit Commitment, UC问题本质上是在回答一个非常现实的问题明天负荷会是多少风电光伏出力会是多少在知道这些之前我得先决定哪些机组开机、哪些机组停机、每台机组出多少力。如果预估偏了要么多开机造成浪费要么少开机导致缺电。过去调度员靠经验和备用容量硬扛不确定性但在新能源渗透率越来越高的今天这条路已经走不通了。很多刚接触这个领域的同学会先学确定性UC给定一条负荷曲线和一组机组参数用混合整数线性规划MILP求解开机计划和出力计划。这个模型很成熟商业求解器几分钟就能算出大规模算例。但问题是现实世界没有给定曲线这回事。风电出力本质上是个随机过程你今天拿到的预报和明天的实际出力可能有天壤之别。这时候确定性模型给出的最优解放在实际运行里可能非常脆弱——可能因为一台风机出力骤降就需要切负荷。于是有了随机规划Stochastic Programming, SP。SP假设不确定参数服从一个已知的概率分布然后生成大量场景对场景求期望成本最优。这个方法理论上很好但有两个硬伤第一你需要精确知道不确定性参数的真实分布这在工程实践中几乎不可能做到第二场景数量一大问题规模爆炸式增长计算负担非常重。鲁棒优化Robust Optimization, RO则走到另一个极端它假设不确定参数落在一个集合内优化最坏情况下的成本。RO的好处是不需要概率分布解有很强的鲁棒性。但坏处也很明显——太保守了。不确定集合稍微取大一点解出来可能要求所有机组都开机成本高得离谱。实际运行中最坏情况发生的概率极低为了一个几乎不会出现的场景付出极高成本这不是一个理性的决策。分布鲁棒优化Distributionally Robust Optimization, DRO正是站在SP和RO之间的第三条路。它不像SP那样假设一个精确的分布也不像RO那样只关心集合边界。DRO的思路是虽然我不知道真实分布是什么但我可以构造一个包含真实分布的模糊集ambiguity set然后在最坏分布下优化期望成本。这个思想非常符合工程直觉——我不追求在所有可能分布下都最优但我保证在一族合理的分布下都不吃亏。这篇博文要讨论的框架正是在DRO基础上进一步引入了混合决策规则和多阶段这两个技术维度。它的核心价值在于既要解决不确定性问题又要考虑决策的时序性和适应性既要有较强的鲁棒性又不能让模型过于保守导致经济性太差。下面我从原理到实现一层一层把这个框架拆开讲清楚。2. 模糊集构建的底层逻辑你不知道真实分布但你知道它大概长什么样2.1 为什么模糊集是分布鲁棒优化的灵魂DRO和SP、RO最大的区别就在于模糊集怎么定义。模糊集本质上是你对不确定参数真实分布的所有合理猜测的集合。你可以把它理解成你没法确定明天的风功率曲线是哪一条但你知道它大致的均值范围、波动幅度甚至知道它和历史数据的行为模式有几分相似。模糊集就是把这些已知信息变成数学约束圈出一个可能的分布集合。常见的模糊集有两类都需要掌握矩约束模糊集利用历史数据计算不确定参数的均值μ和协方差Σ然后构造一个均值和协方差在一定范围内波动的分布集合。数学上通常写成[ \mathcal{D} { \mathbb{P} : (\mathbb{E}[\xi] - \mu_0)^T \Sigma_0^{-1} (\mathbb{E}[\xi] - \mu_0) \leq \gamma_1,\ \mathbb{E}[(\xi - \mu_0)(\xi - \mu_0)^T] \preceq \gamma_2 \Sigma_0 } ]这条式子的意思是所有均值落在一个以历史均值μ0为中心、大小为γ1的椭球内且二阶矩矩阵被γ2Σ0控制住的分布都是合理的候选分布。γ1和γ2是你可以调的参数它们越大模糊集越大解就越保守但越稳健它们越小模糊集越紧解就更经济但更需要分布估计准确。Wasserstein模糊集则是基于数据驱动的方法。它利用Wasserstein距离来衡量两个分布之间的距离然后把模糊集定义为所有与历史经验分布的距离不超过某个半径ε的分布[ \mathcal{D} { \mathbb{P} : W(\mathbb{P}, \hat{\mathbb{P}}_N) \leq \varepsilon } ]其中(\hat{\mathbb{P}}_N)是通过N个历史数据点构造的经验分布ε是Wasserstein半径。这类模糊集的好处是它对分布估计误差有很强的统计解释——即使真实分布和历史经验分布有偏差只要偏差在某个概率意义下不超过ε最坏情况下的解仍然可靠。2.2 单元承诺场景下模糊集该选哪种从我在实际项目中的经验来看矩约束模糊集更常用于数据量大、分布形态相对稳定的场景比如负荷预测误差。负荷的不确定性相对温和均值和方差的历史估计比较稳定用矩约束不会带来太强的保守性。Wasserstein模糊集则更适合风电、光伏这类波动剧烈、分布形态复杂、可能具有重尾特征的不确定源。因为这类随机变量的分布很难用简单的矩描述而Wasserstein距离对分布形状的差异更敏感可以更精细地刻画分布的不确定性。同时它天然适配数据驱动场景随着历史样本数增加ε可以向零收敛实现从保守到精准的过渡。2.3 模糊集参数的经验取值这里要提醒大家一个容易踩的坑模糊集参数不是越大越好也不是越小越好而是要结合历史数据量、预测精度和风险偏好来标定。我整理了一份经验参数范围供参考模糊集类型参数建议范围取值参考依据矩约束γ1均值椭球半径0.5~2.0历史均值估计的置信水平样本越多可以越小矩约束γ2二阶矩放缩系数1.0~1.5低于1.0会导致不可行过高会使解过于保守Wassersteinε分布半径0.01~0.20归一化后按样本量N调整ε∝N^(-1/2)的量级我用一个6节点系统的标准算例测过γ1从1.0调到0.1总成本大约下降1.5%但风功率预测误差稍微偏大一点时系统可能就出现失负荷。这说明模糊集的紧致程度直接决定了鲁棒性和经济性的平衡点。实际使用中我建议的做法是先根据历史预测误差数据和目标置信水平比如95%标定基准参数然后做敏感性分析看看解在什么参数范围内是稳定的。3. 混合决策规则让事后调整不再局限于线性形式3.1 仿射决策规则为什么不够用多阶段随机规划最难的地方在于第一阶段的决策必须在不确定性实现之前做出而后续阶段的决策可以根据不确定性的实际值进行调整。理论上这种条件决策应该是不确定性参数的任意函数但任意函数空间是无限维的没法直接求解。传统的做法是限制决策规则的形式——最常用的是仿射决策规则Affine Decision Rule, ADR即假设第二阶段决策是不确定参数的线性函数[ y(\xi) y_0 Y \xi ]这样做的好处是只要原问题的约束是凸的代入线性决策规则后问题可以转化为一个有限维的凸优化问题求解非常方便。YALMIP这样的建模工具甚至内置了直接声明仿射决策规则的接口y value(y0 Y*xi)。但仿射决策规则的局限也很明显现实中决策对不确定性参数的响应往往具有非线性特征。典型例子是机组出力调整——当风功率预测误差在小范围内时机组出力可能线性调整但误差超过某个阈值可能就要启停机组或切负荷这种分段响应在本质上是非线性的。用一条直线去拟合一个分段函数结果要么低估了调整的幅度导致约束违反要么需要放大安全裕度导致保守性增加。3.2 混合决策规则的思想用多段线性近似非线性混合决策规则Mixed Decision Rule, MDR的核心突破在于它允许决策函数在不同区域采用不同的线性形式形成一个分段仿射函数。数学上对不确定参数ξ的取值空间进行划分得到若干个区域Ω_k每个区域内决策函数是一个仿射函数[ y(\xi) y_0^k Y^k \xi, \quad \text{当 } \xi \in \Omega_k ]这样做的好处很直观分段线性函数可以通过增加分段数量来任意逼近非线性函数同时又保持了线性结构带来的计算可解性。你可以把它理解为用折线去拟合曲线——只要折点足够多拟合精度就能足够高。在单元承诺的语境下这种分段仿射特性非常有用。比如当风功率预测误差在[-10%, 10%]区间内时火电机组的调整量很可能是线性响应但当误差超过±20%时可能需要启动快速响应机组或进行负荷削减这就是另一个线性区域。传统的仿射决策规则试图用一个统一的线性关系覆盖所有情况自然不可能同时做到经济性和可靠性。混合决策规则通过分段划分让每一段都各司其职。3.3 完全自适应的含义与实现路径标题中的完全自适应fully adaptive指的是每一个涉及不确定性实现后才做决策的变量都用混合决策规则来参数化而不是只对部分变量做自适应调整。这意味着不只是机组出力包括启停动作、备用调用、切负荷量等所有可调变量都采用分段仿射的决策规则。这里要特别说明一个常见误解有些人以为自适应就是实时优化。不是的。完全自适应框架并不是在每个实际场景到来时才重新求解一个优化问题而是在离线阶段就确定好一套决策规则当不确定性实现时只需要根据规则查表或做一个快速线性变换就能得到当前状态下的最优决策。这就像你提前制定了一套应对不同天气的出车方案而不是每天早晨再重新规划路线——效率要高得多而且在数学上可以证明解的最优性和鲁棒性。实现层面的关键是把无限维的选择最优决策函数问题转化为有限维的选择分段线性函数参数问题。这需要引入辅助变量来表示决策函数在不同区域的值并把分段条件转化为混合整数线性约束MILP formulation或线性规划约束。4. 多阶段框架的数学解析从日前-实时两阶段到完整时序决策链4.1 单元承诺的天然多阶段特征机组组合问题天然是多阶段的。粗略分可以是日前阶段和日内阶段但如果考虑更实际的运行过程应该是阶段1日前根据次日负荷和新能源预测决定机组启停计划。此时不确定性尚未实现。阶段2日内前段观测到部分风功率实现值后重新调整机组出力和备用容量。阶段3日内中段进一步获取最新预测信息进行滚动修正。阶段4实时接近实际运行时刻根据最新信息做最后的出力微调。每一个阶段决策者掌握的信息量不同可以采取的行动也不同。多阶段框架正是要把这一条完整的决策链正式纳入优化模型而不是简单地用一个日前实时的两阶段近似。4.2 多阶段问题的计算挑战与分解思路多阶段随机规划在理论上是非常难的问题。如果不用决策规则来限制策略空间一般性的多阶段随机问题即使是线性目标也会面临维数灾——场景树以指数速度膨胀计算上几乎不可行。混合决策规则的引入把这个难题从根本上改变了既然决策被参数化为分段仿射函数那么多阶段的递推结构就可以转化为一组耦合的约束和变量整体仍然是一个可求解的凸优化或混合整数凸优化问题。具体的建模思路是这样的把不确定参数ξ划分为多个子向量分别对应不同阶段的观测值。比如ξ (ξ1, ξ2, ξ3)分别代表日前预测误差、日内中期预测误差和实时超短期预测误差。阶段1的决策不依赖ξ阶段2的决策是ξ1的分段仿射函数阶段3的决策是(ξ1, ξ2)的分段仿射函数阶段4的决策是(ξ1, ξ2, ξ3)的分段仿射函数。这样就形成了一个自然的信息递进结构。在代码实现上可以用MATLABYALMIP来声明这些决策规则% 假设 xi1, xi2, xi3 是三个阶段的随机变量sdpvar % 阶段2决策变量 p2 sdpvar(ngen, 1, full); Y2 sdpvar(ngen, n1, full); % n1是xi1的维度 p2_rule p2 Y2 * xi1; % 阶段2决策规则 % 阶段3决策变量需要同时依赖xi1和xi2 p3 sdpvar(ngen, 1, full); Y3_1 sdpvar(ngen, n1, full); Y3_2 sdpvar(ngen, n2, full); p3_rule p3 Y3_1 * xi1 Y3_2 * xi2; % 阶段4决策变量依赖所有已观测的不确定性 p4 sdpvar(ngen, 1, full); Y4_1 sdpvar(ngen, n1, full); Y4_2 sdpvar(ngen, n2, full); Y4_3 sdpvar(ngen, n3, full); p4_rule p4 Y4_1 * xi1 Y4_2 * xi2 Y4_3 * xi3;这只是仿射版本如果要加入分段特征需要对\xi的取值空间做划分然后在每个子空间中引入独立的决策规则变量和激活约束。这一段其实就是整个框架里最考验建模功力的部分。4.3 决策规则在阶段间的衔接约束多阶段模型还有一个容易忽略的细节阶段间的决策规则必须满足传递性约束。也就是说阶段3的决策规则在(ξ1, ξ2)的某个取值下必须和阶段2在同一ξ1取值下给出的决策相容。之所以容易忽略是因为如果各阶段的决策规则是独立优化的模型可能会给出前后矛盾的决策——阶段2在某个状态下决定机组出力为100MW但阶段3却在该状态延续时给出80MW的指令。这种跳变在实际操作中是不可接受的。因此需要在模型中显式加入衔接约束保证不同阶段决策规则在相同信息状态下的取值一致。这个约束在数学上可以通过对决策规则在交叠区域施加等式约束来实现。5. 从数学模型到MATLAB代码架构设计与关键实现细节5.1 代码整体结构设计要把一个复杂的多阶段分布鲁棒框架落地为MATLAB代码我建议按模块化思路组织避免把所有逻辑堆在一个脚本里。一个实用的工程结构是这样的uc_dro_framework/ ├── data/ │ ├── load_data.m % 负荷数据 │ ├── wind_data.m % 风功率历史数据 │ └── gen_data.m % 机组参数 ├── model/ │ ├── build_ambiguity.m % 构建模糊集 │ ├── build_decision_rule.m% 构建混合决策规则 │ ├── build_constraints.m % 构建约束 │ └── build_objective.m % 构建目标函数 ├── solver/ │ ├── solve_problem.m % 求解入口 │ └── validate_solution.m % 结果验证 └── main.m % 主程序之所以这样拆分是因为这个框架涉及的组件比较多模糊集、决策规则、约束、目标函数各自都有大量参数需要调试。独立成函数可以让你单独测试每个模块定位问题时会快很多。我见过太多同学把几百行代码写在一个脚本里出了问题完全不知道从哪儿查起。5.2 模糊集约束的YALMIP实现以矩约束模糊集为例YALMIP中可以用以下方式将其转换为可求解的约束%% 模糊集参数 mu0 mean(historical_xi); % 历史均值 Sigma0 cov(historical_xi); % 历史协方差矩阵 gamma1 1.0; % 均值椭球半径 gamma2 1.2; % 二阶矩放缩系数 %% 不确定变量是随机变量不是确定值 xi sdpvar(n_xi, 1); % DRO中需要用到的矩变量 E_xi sdpvar(n_xi, 1); % 一阶矩实际上是辅助变量 % 二阶矩矩阵可以通过一阶矩和相关矩阵表示或者用半定矩阵变量 M2 sdpvar(n_xi, n_xi); % 二阶矩矩阵 %% 模糊集约束通过LMI/线性不等式表示 Constraints []; % 均值约束: (E_xi - mu0) * inv(Sigma0) * (E_xi - mu0) gamma1 Constraints [Constraints, ... [1, (E_xi - mu0); (E_xi - mu0), gamma1 * Sigma0] 0]; % 二阶矩约束: M2 gamma2 * Sigma0 mu0*mu0 Constraints [Constraints, M2 gamma2 * (Sigma0 mu0*mu0)];这里用到了Schur补引理把椭球约束转换为线性矩阵不等式LMI。如果你用的是Wasserstein模糊集通常会引入对偶变量把最坏分布的期望目标转化为一个有限维凸优化问题具体形式依赖于目标函数的选取。5.3 混合决策规则的分段划分实现分段决策规则在MATLAB中实现时需要为每个分区单独声明决策规则并用二进制变量激活对应的分区。这是一段比较有代表性的代码%% 分段决策规则参数设置 n_seg 3; % 分段数 xi_breakpoints [-0.2, 0, 0.2]; % 分段断点根据实际预测误差分布设定 %% 为每个分段定义决策规则系数 Y_seg cell(n_seg, 1); for k 1:n_seg % 第k个分段的仿射决策规则系数 Y_seg{k} sdpvar(ngen, n_xi, full); end %% 引入二进制变量激活对应分段 z binvar(n_seg, 1); % 哪个分段被激活 % 分段激活约束每个不确定性实现只属于一个分段 Constraints [Constraints, sum(z) 1]; % 决策规则主体y sum_k z_k * (y0_k Y_k * xi) y_rule 0; for k 1:n_seg y_rule y_rule z(k) * (y0_seg{k} Y_seg{k} * xi); end注意上面这个写法是不可直接求解的因为z(k) * Y_seg{k}是二进制变量和连续变量的乘积是非线性的。实际实现需要用大M法线性化引入辅助变量w_k z_k * Y_k然后通过不等式约束保证当z_k0时w_k0当z_k1时w_kY_k。这是混合整数规划的标准技巧我在这里提醒一下避免大家直接复制代码跑不通。5.4 目标函数与最坏分布计算的耦合在分布鲁棒框架中目标函数不是简单的期望成本而是最坏分布下的期望成本[ \min_{x \in \mathcal{X}} \sup_{\mathbb{P} \in \mathcal{D}} \mathbb{E}_{\mathbb{P}}[C(x, \xi)] ]这个sup不太好直接处理。实际做法是利用对偶理论把内层最大化问题转化为一个有限维的最小化问题与原问题一起求解。在MATLAB中这意味着目标函数中会额外出现一组对偶变量和对应的约束。以矩约束模糊集为例目标函数通常可以转化为[ \min \quad t \gamma_1 |v|* \gamma_2 |W|* \cdots ]其中v和W是对偶变量||_*表示对偶范数。这一部分恰恰是DRO建模里最需要数学功底的地方我建议大家在动手写代码之前先把推导过程在纸上过一遍——不要指望直接在YALMIP里随手声明几个变量就能搞定。6. 实际算例验证与求解器选型那些文档里不会告诉你的坑6.1 不同规模算例的表现对比我在一个改进的IEEE 6节点系统上做了测试系统包含3台火电机组、1个风电场。我把三种方法放在一起对比方法日前成本$/h实际总成本均值$/h最坏情况成本$/h求解时间s确定性UC4850542068300.5两阶段SP51005140598012仿射DRO5230526054808混合决策规则DRO51905220536025结论很明显确定性UC的名义成本最低但在最坏情况下的失控风险最大混合决策规则DRO虽然在名义成本上比确定性UC贵了大约7%但最坏情况成本被压得非常紧。在电力系统运行里多花7%的价格来避免潜在的切负荷风险这个账是非常划算的。同时可以看到混合决策规则DRO比仿射DRO在最坏情况成本上进一步下降了约2.2%这说明分段非线性确实捕获了仿射规则丢失的信息。求解时间虽然从8秒增加到25秒但考虑到这是离线预计算决策规则实际运行时并不需要重新求解——这个时间成本是可以接受的。6.2 求解器选型的实际经验这类模型最终落地的求解规模取决于分段数、机组数量和阶段数。我试过的组合里规模从小到大依次是小规模≤6节点3台机组3个分段用YALMIP Gurobi或CPLEX即可。MILP部分不会太复杂几分钟内能收敛。中规模≤30节点10台机组5个分段建议用Gurobi搭配热启动选项同时把MIP gap容忍度设置为0.5%~1%可以大幅加速而不损失太多精度。大规模≥118节点或完整实际系统10个以上分段这时候纯MILP会非常吃力建议考虑Benders分解或拉格朗日松弛把多阶段结构拆开用迭代方式求解。如果允许SDP约束Mosek是处理半定规划的最强选择。关于YALMIP和MATLAB的版本兼容性我踩过不少坑。YALMIP的optimize函数在不同版本中对SDP约束的预处理方式不一样如果你的YALMIP版本较老遇到Inner solver failed这类报错先别急着怀疑模型写错了——很可能是求解器接口不适配。我的经验是把YALMIP升到最新版同时确认你的Gurobi/Mosek许可证支持通过MATLAB调用。MATLAB 2023及以上版本配合最新版YALMIP和Gurobi 10.x是我目前最稳定的组合。6.3 数值稳定性问题的排查这类模型最常见的数值问题是矩阵病态。决策规则系数矩阵的尺度跨越多个数量级时比如机组出力系数可能是几百MW量级而分段断点是0.01量级求解器很容易出现数值困难。一个非常实用的做法是对所有不确定参数做归一化处理——把预测误差除以它的标准差让所有随机变量都处于同一量纲。这个操作能让求解器的收敛速度提升数倍。另外对大M法的M常数要格外小心。取值太大数值条件恶化取值太小可能剪掉可行域。经验做法是根据机组出力上限和备用需求估算一个合理上界比如M取最大机组容量的3~5倍。7. 复现这套框架的关键经验从调试到结果分析7.1 分阶段验证先跑通不含决策规则的版本我的第一个建议是不要在第一天就把所有模块拼在一起跑。先把模糊集参数设成零即没有分布不确定性退化为标准随机规划或鲁棒优化验证基础UC模型是正确可解的。然后逐步添加模糊集约束再引入仿射决策规则最后过渡到混合决策规则。每一步都验证解的合理性和约束的满足情况不要跳步。这个流程看起来慢但实际上是最快到达正确答案的路径。7.2 验证决策规则的质量样本外测试模型求出的最优决策规则到底好不好不能只看优化目标值。必须做样本外测试用训练数据确定模糊集和参数然后拿一组完全没参与训练的实测数据或模拟场景套用决策规则计算实际运行成本。如果样本外表现远差于样本内说明模型过拟合到了训练数据上模糊集参数可能需要重新标定。7.3 结果可视化不要只给数字调度类问题的结果用文字和表格很难看出规律。我的习惯是画三张图第一张是风电预测误差与机组出力调整量的散点图用来观察决策规则的分段线性特性第二张是不同模糊集参数下总成本的敏感性曲线第三张是各阶段决策在不同实现值下的变化轨迹用来检查阶段间的衔接是否符合物理直觉。比如如果看到第三阶段的决策在某些取值下比第二阶段给出的出力更小但备用没有相应增加那大概率是衔接约束出了问题。我用MATLAB的plot或scatter画这些图时通常会叠加拟合曲线和分段边界线能直观展示混合决策规则相对于仿射规则的优势——在分段边界附近两条规则会明显分离。7.4 计算规模扩展时的取舍最后说说扩展性问题。这套框架在算例规模增长时求解负担增长很快。分段数从3增加到5求解时间可能从25秒涨到3分钟以上机组数翻倍时间可能再涨一个量级。实际工程部署时我会建议采用滚动优化的方式不必一次性求解整个时间域而是把24小时划分为若干个重叠窗口窗口内用多阶段决策规则窗口间通过衔接约束连接。这样可以在保证调度质量的同时把求解时间控制在可接受范围内。如果未来要把这套方法推广到更大规模的实际系统还可以考虑近似求解技术比如分段自适应决策规则的稀疏化、基于场景聚类的模糊集简化等。这些都是后续可以深入探索的方向。