基于线性决策规则的分布鲁棒优化机组组合Matlab实现

搞电力调度的人最怕的就是风电出力预测不准带来的连锁反应。你按照预报值把明天96个时段的机组开机计划排好，后半夜风突然小了，火电还顶着最小技术出力，系统频率往下掉，备用被吃掉一大块；后半夜风又突然大了，火电没办法那么快压下来，只能眼睁睁看着弃风。这种场景在新能源渗透率高的电网里几乎是常态。传统确定性机组组合对这种随机性基本没有招架之力，于是大家开始把不确定性处理直接塞进调度模型。这篇文章要分享的，就是目前学术和工程界都很有分量的做法：基于线性准则（线性决策规则，LDR）的分布鲁棒优化机组组合，并用Matlab完整落地实现。它解决的核心问题是"风电出力到底怎么波动、误差分布到底长什么样"这个不确定性源头下的最优决策，适合电力系统优化方向的研究生、做调度算法开发的工程师、以及对两阶段鲁棒优化感兴趣的程序员参考。

1. 从确定性机组组合到不确定性环境下的建模逻辑

1.1 经典机组组合的本质和痛点

机组组合（Unit Commitment，简称UC）是电力日前调度里最核心的优化问题，本质上是问两件事：明天哪些机组要开机，以及开着的机组每个时段发多少电。决策变量分两类，一类是0/1整数变量，代表机组启停状态；另一类是连续变量，代表每台机组在某个时段的出力。目标函数通常由启动成本、空载成本和燃料成本三部分组成，约束条件包括功率平衡、旋转备用、出力上下限、最小开停机时间、爬坡速率等等。标准模型是混合整数线性规划或混合整数二次规划，规模一放大到IEEE 118节点这种量级，变量数量就奔着十万以上去了，对建模质量和求解器性能都非常有考验。

我的体会是，经典UC最常见的坑反而不是模型本身难建，而是参数设置不严谨。比如爬坡约束的单位问题，冷启动和热启动成本的区分，最小开停机时间的递推初始化，这些细节稍微写错，求解器就会给出一个看似合理、实则无法落地的调度计划。所以在做不确定性问题之前，我强烈建议先把确定性UC在几个标准算例上跑通，再往上叠加不确定性建模，否则出了问题都没法定位，到底是基础模型写错了，还是鲁棒化过程引入的bug。

1.2 风电不确定性让调度变得多难

风电出力的不确定性来自天气过程的随机性。预测模型给一条期望出力曲线，但实际出力会在预测值附近波动，波动幅度大时能到额定容量的百分之二三十。当风电场数量多、分布广时，预测误差之间还带有空间相关性，整个系统的总误差不再是单机误差的简单叠加。把这些不确定性写进机组组合模型，意味着排计划时必须给系统留出足够余量，既要能应对风电出力不足时的缺电风险，也要能应对出力盈余时的调峰压力。

确定性模型处理这种问题的传统办法是加固定比例备用，比如让开机机组的可用容量比负荷预测高一个固定的百分比。这招简单有效，但在风电占比高的系统里，固定比例备用要么太保守，要么太激进。更关键的是，备用约束只卡住容量这一个维度，根本不关心风电误差的分布形态，也就没法回答"最有可能发生的坏事是什么、概率有多大"这类问题。要做更精细的调度，必须把风电的随机特性真正放进优化模型，而不是简单加一个备用系数了事。

1.3 三种主流处理不确定性的策略对比

学术界处理调度不确定性的主流框架有三类：随机优化、传统鲁棒优化、分布鲁棒优化。我用一个表格先把它们放在一起对比，后面再做详细解释。

用一个生活化的类比帮助理解。随机优化就像出门前看天气预报说下雨概率70%，你带把伞；传统鲁棒优化则是不管天气预报怎么报，都穿全套雨衣雨靴，因为"最坏情况就是下大暴雨"；分布鲁棒优化介于两者之间——参考天气预报，但同时又考虑"万一这个天气预报本身不太准，真实天气状况可能会偏到哪个方向"，然后针对这种不确定性再做一个更保守的决策。三种方法没有绝对优劣，关键看场景、数据条件和保守度接受程度。

1.4 引入线性决策规则的动机

在风电不确定环境下的机组组合，天然是一个两阶段决策问题。第一阶段在日前提前安排机组启停；第二阶段等实际风电出力揭晓后，再对火电机组出力做实时调整。严格来说，第二阶段决策应该是"实际风电出力向量"到"调整量"的任意函数关系，这个函数空间是无限维的，求解器根本没有办法直接处理。

线性决策规则（在标题里叫"线性准则"）就是把这个任意函数限制成仿射函数。举个例子：火电机组i的实时调整量写成dg_i = a_i + b_i乘以风电预测误差，其中a_i和b_i是待求的常数系数。这样一来，第二阶段策略就只由有限个系数完全刻画，无限维的函数空间自然塌缩成一个有限维的系数空间。从电力系统运行的角度看，这也不是一个脱离实际的假设，调度中心安排机组参与AGC调节时，本来就是按偏差比例分配调节量，本质上就是一种线性规则，所以线性决策规则这类方法不仅仅是数学上的权宜之计，它也有合理的工程背景。

2. 分布鲁棒机组组合的模型与线性决策规则

2.1 模糊集：给概率分布"上保险"

分布鲁棒优化要解决的第一个问题是：我该怎么描述"真实分布到底在哪里"。答案就是构造一个集合，把可能的概率分布都圈在集合里，这个集合就叫模糊集。工程上最常看到的是两类模糊集：矩约束模糊集和Wasserstein距离模糊集。

矩约束模糊集限制分布的均值、协方差等矩信息在一定范围内。它的好处是表达式简洁，对偶转换相对直接；缺点是只刻画低阶矩，很多形态完全不同的分布在同一个模糊集里是难以区分的。Wasserstein距离模糊集则是以历史样本构造的经验分布为中心，以某个半径ε为半径的Wasserstein球。从几何意义上讲，Wasserstein距离度量的是把一个分布"搬移"成另一个分布的最小运输成本，所以这个模糊集可以直观控制"真实分布和样本分布差多远"。

实际项目里我更倾向于用Wasserstein模糊集，因为它直接吃历史样本，不需要假设分布形态，而且半径的大小和鲁棒程度是一一对应的，调参方向非常清楚。统计学里还有一个有用的结论：样本量越大，真实分布越可能落在一个半径跟样本量倒数平方根同阶的Wasserstein球里。这可以作为初始半径的上界参考，后面再用交叉验证修正。

2.2 两阶段模型与min-max-min结构

加入模糊集之后，分布鲁棒机组组合的完整结构是一个三层的min-max-min问题。最外层min对应第一阶段的机组启停决策，这部分在日前阶段必须提前敲定。最内层min对应第二阶段实时调整，即风电出力揭晓后，在已知基点出力的前提下微调各机组。夹在中间的max是模糊集内最坏分布的选取，我们希望在"最不利的现实分布"下，期望运行成本仍然尽可能低。三层套在一起，直观的理解就是：启停方案可以先定下来，但必须保证无论现实按哪种合理分布发展，通过最优实时调整都能把成本控制在可接受范围。

单独看这个min-max-min结构，直接求解是不可能的。求解器的能力边界在有限维确定性规划，所以必须想办法把无穷维的分布变量和函数变量消掉。模糊集的对偶变换对应了max这一层，线性决策规则对应了内层min这一半，两个工具搭在一起，才能把整个问题压成可解的确定性规划。这个点理解透了，后面看代码和推导都会顺畅很多。

2.3 线性决策规则如何化简模型

用线性决策规则替换第二阶段决策后，实时调整量变成不确定性变量ξ的仿射函数：d = D0 + D1乘以ξ，其中D0是常数向量，D1是系数矩阵，都是待优化变量。把这种表达式代入功率平衡约束、爬坡约束、出力上下限约束以后，每条约束都变成含随机变量ξ的线性不等式。

对等式约束来说，因为要求"对任意ξ都成立"，所以ξ的常数项和一次项系数必须分别满足等式约束，这个处理比较直接。对不等式约束来说，则需要做鲁棒化处理，把"对任意ξ都满足不等式"的要求转换成一组确定的约束。在分布鲁棒框架下，这个转换通常会引入与模糊集半径相关的惩罚项，惩罚项最终会和目标函数的期望成本联合在一起，形成新的确定性表达式。如果你不用分布鲁棒而用盒式支撑集，那其实就是退化成了传统鲁棒优化，所有点都要求满足，表达式更简单但结果更保守。

2.4 对偶转换后得到的确定性MILP

经过模糊集对偶和线性决策规则替换这两个步骤，原来的min-max-min三层结构被压成一个确定性的混合整数线性规划。目标函数变成：第一阶段启动与空载成本，加上第二阶段期望成本的上界；约束条件包含原本的机组运行约束、LDR系数相关约束，以及对偶变换产生的辅助变量约束。这个确定性MILP可以直接交给Gurobi或CPLEX求解。

这里要说明一点，对偶转换的数学推导步骤比较多，实际建模时可以利用YALMIP这类建模工具减少手算负担，但前提是你要清楚模型里哪些变量是决策变量、哪些是不确定参数、哪些是辅助对偶变量。如果只是照抄别人的代码，一旦遇到无解或者结果不对，很难定位问题。我自己做这个项目时，最开始就是手推一遍小规模的6节点算例，然后才在Matlab里扩展，后续调试才没那么痛苦。

3. Matlab实现：从数据准备到求解器调用

3.1 代码架构与文件模块划分

这类优化代码最忌讳所有东西堆在一个脚本里。我建议把项目拆成几个模块：main主程序管整体流程，case_data负责定义机组参数和负荷风电数据，build_deterministic_uc负责构建确定性UC模型，build_dro_constraints负责叠加模糊集、LDR系数和对偶约束，solve_milp负责调用求解器并检查状态，plot_results负责出图。文件划分清晰了，后面改参数、换数据、对比不同模糊集半径时能省下大量时间。

模块化还有一个隐藏的好处：调试时可以先用deterministic版本验证基础模型完全正确，再打开DRO开关。出问题时定位范围一下就缩小了——是先确认基础模型能跑通，再去研究新增的分布鲁棒约束哪里写得不对。

3.2 数据准备：风电样本与机组参数

数据准备第一件事是统一单位。风电出力、负荷、机组容量尽量都换算成同一个基准值。我一般用100 MVA做标幺化，这样目标函数数值不至于因为数量级差异在求解器里出现数值病态。我踩过一次很深的坑：负荷数据用MW，风电数据也用了MW，但忘记了把某个风电场的历史数据减去基准值，导致功率平衡约束在数值上总有几百的残差，求解结果画出来曲线看不出问题，一校验才发现特定时段严重不满足平衡。

风电不确定性样本的构造方式取决于你手上的数据。如果有历史实测出力，可以直接把预测误差作为样本；如果只有预测和实测两条曲线，可以按小时构造误差样本，并考虑时序相关性。具体到线性决策规则，稳妥的做法是把每个时段的风电误差都作为多维随机变量ξ的一个分量，样本集就是历史误差向量的集合。用经验分布作为模糊集的中心，再结合交叉验证选半径，要比随手指定一个半径可靠得多。如果原始数据里有明显的坏数据或者通信中断导致的零值，记得先做数据清洗，否则模糊集中心会被垃圾样本带偏。

3.3 YALMIP建模与求解器配置

Matlab里最省事的优化建模路线是YALMIP加Gurobi或CPLEX。YALMIP负责把数学表达式翻译成求解器输入，Gurobi负责解混合整数线性规划。安装和配置不展开说了，重点讲代码逻辑和容易出错的地方。

%% 主流程核心代码片段（YALMIP + Gurobi） nUnits = 3; % 机组数 nPeriods = 24; % 时段数 nSamples = 500; % 风电误差样本数 u = binvar(nUnits, nPeriods, 'full'); % 机组启停 p0 = sdpvar(nUnits, nPeriods, 'full'); % 基点出力 a0 = sdpvar(nUnits, nPeriods, 'full'); % LDR常数项系数 A1 = sdpvar(nUnits*nPeriods, nSamples, 'full'); % LDR对风电误差的线性系数，实际建模需按系统维度展平 %% 目标函数：启动成本 + 空载成本 + 燃料成本 + 分布鲁棒惩罚项 objective = startupCost(u) + noLoadCost(u) + fuelCost(p0) + droPenalty(a0, A1, radius, samples); %% 约束：功率平衡、机组限值、爬坡、最小启停时间、分布鲁棒约束 constraints = [powerBalance, unitLimit, rampLimit, minUpDownTime, droConstraints]; ops = sdpsettings('solver', 'gurobi', 'verbose', 2, ... 'gurobi.MIPGap', 0.001, 'gurobi.TimeLimit', 600); optimize(constraints, objective, ops);

代码的核心是变量层级要规划清楚。u和p0是第一层变量，直接出现在目标函数和确定性约束里；a0和A1是线性决策规则层的系数，几乎每个与风电误差相关的约束里都要用到。写代码时要特别注意A1的索引顺序，别把时段、机组和样本的顺序搞混。我建议每定义一个变量都写清注释，否则这类三维以上的数据展开在调试时真的会让人怀疑人生。分布式鲁棒惩罚项在代码里通常不是一个简单的函数，而是通过对偶辅助变量构成的线性项，放在独立函数里更利于测试。

3.4 结果分析与后处理要点

模型求解完之后别急着出图，按顺序做四件事。第一，检查求解器返回状态，确认MIPgap是否在可接受范围。第二，检查机组启停状态是否满足最小开停时间约束，有时候初始状态设置不当会出现机组连续启停的不合理序列。第三，把最优解代回原始的两阶段问题做蒙特卡洛回验，在随机风电样本下统计系统实际的失负荷量和弃风量。这个回验非常重要，因为它才能说明这个分布鲁棒调度方案在模拟真实运行中到底靠不靠谱。第四，做不同Wasserstein半径的敏感性分析，画出"半径—期望成本/最坏成本"曲线，看结果的鲁棒性如何随保守程度变化。

后处理是论文里最有说服力的部分。我的习惯是把确定性UC、传统鲁棒优化、DRO不同半径的三组方案放在同一张图上，横轴是风电出力的波动幅度，纵轴是系统运行成本或者失负荷量。这样能直观展示DRO相对传统方法的优势和成本代价，也能在审稿人面前展现出对模型特性的深入理解。

4. 调试经验与常见问题速查

4.1 求解慢的排查方案

求解慢是这类项目最常见的难题，但"慢"的原因往往不在模型，而在于建模质量。我按经验列一个排查顺序。第一，看MIP下界是不是长时间不动。如果下界一点改善都没有，多半是约束里大M常数太多或者太大，尝试把所有大M压缩到"恰好够用"的量级，求解速度通常立刻提升。第二，看整数变量规模是不是过大。机组数量多时，可以尝试对称性破缺，把同类型机组合并成机组簇，或者添加对称性破缺约束，能有效压缩搜索空间。第三，调整求解器参数。Gurobi的MIPFocus可以设为1或2，Threads适当增加，对中小规模算例一般有正向收益。

分布式鲁棒模型求解慢还有一个来源：模糊集半径设置过大，导致对偶辅助变量数量膨胀，线性规划松弛质量变差。如果你只是做原理验证，先在3机6节点系统上测试通，再逐步放大到IEEE 30节点和118节点，别一上来就拿大系统调参，那样很容易被漫长的求解时间磨掉耐心。

4.2 Wasserstein半径的调节技巧

半径的选择本质上是风险偏好的量化。太小，模型几乎完全信任样本分布，相当于退化成随机优化；太大，模型倾向于为极端概率分布做准备，成本一路走高。我的调节思路是两步走。先用较粗的网格扫描半径，比如从0.01到0.5，观察成本曲线变化，找到曲线斜率变化最明显的拐点；然后在拐点附近加密扫描。拐点处的半径通常是性价比最高的选择——再增大半径，成本大幅上升但鲁棒性收益有限；再缩小半径，鲁棒性明显下降但成本节省不多。

如果手头数据是真实风电历史出力，半径还可以参考样本量做上界估计。样本量越大，经验分布越接近真实分布，半径应取得越小。这种基于统计直觉的调参方式写进论文里，比"人为设定半径"要显得扎实得多。另外提醒一句，每次改完半径，都要重新跑一遍蒙特卡洛回验，因为纸上谈兵的"最坏分布成本"和实际随机样本下的表现不一定完全对应。

4.3 无可行解问题分析

遇到无可行解，我的习惯不是盯着报错信息看，而是先做"松弛定位"。具体做法是逐个去掉约束组，看看哪一组约束拿掉之后问题重新可行。根据经验，常见原因有这么几个。最小开停机时间和爬坡约束在小系统中特别容易打架，尤其当负荷曲线峰谷差大、机组最小启停时间又长时；备用约束和机组最大出力之间的配合不好，导致某个时段总可用容量低于负荷加备用；还有一个容易被忽略的情况，模糊集半径设得太大、样本里存在极端点，会让对偶辅助变量的可行域被冲散。

另外，线性决策规则系数A1如果没有设置对应的物理上下限，优化器可能会利用变量自由度产生技术可行但实际反物理的调度方案。做可行性测试的时候，把机组出力值打印出来和上下限肉眼比对一下，很多问题一眼就能看出来。如果发现某台机组出力长期贴在边界上，同时启停状态还在频繁切换，大概率是初始状态或者爬坡约束出了问题。

4.4 隐藏在实现细节中的坑

整理几个我实测下来最容易踩的坑。第一，YALMIP中二元变量和连续变量相乘时，需要用big-M做线性化，这个M的值必须基于物理边界设置得尽量紧。如果M太大，线性松弛后的可行域过宽，求解过程会异常缓慢；M太小，又会切掉真实可行解，导致结果不对。第二，Gurobi许可证和Matlab版本之间的兼容性，装不上就换Mosek或者用YALMIP自带求解器，但默认求解器对大规模MILP非常吃力，不建议用在大算例上。第三，做蒙特卡洛后验时，抽样样本数量一定要大于模糊集中心的样本数量，不然你验证的"稳健性"可能只是对少量样本的重放，统计意义不足。第四，所有风电相关数据进入模型前务必统一为标幺值。我因为这个吃过一次大亏，某节点的功率平衡残差显示500多MW，排查了半天，最后发现只是一个风电场的数据漏了基准值换算。

这些坑单个看都不大，但叠加起来足以让一个好端端的模型变成"看起来很对、跑起来就废"的鸡肋。做项目时，最好把数据校验和结果校验写成独立脚本，每次改动模型后自动跑一遍，在早期就能拦住绝大部分问题。

4.5 可以继续扩展的方向

这个模型还有不少值得扩展的空间。实际系统往往同时存在负荷不确定性、风电不确定性，甚至电价和检修计划也是随机的，把多源不确定性都纳入模糊集是更贴近现实的建模方式。另一个有价值的方向是把线性决策规则升级为分段线性决策规则，模型会变大，但保守性会明显降低。还有些团队把分布鲁棒机组组合和储能规划、需求响应联合优化，在计算复杂度可控的前提下有效降低系统运行成本。如果你打算在这个方向发论文，这几个角度比单纯换算例改参数更有含金量。

我个人做完这个模型的体会是，数学推导再漂亮，最终还是要落到代码细节上才见真章。分布鲁棒、线性决策规则这些名词听起来很前沿，但实现过程中的很多问题其实是基础性的：单位换算了没有，变量索引对没对，大M松不松，MIPgap设得合不合理。把这些基础层面的东西打磨扎实，模型自然就能跑得动、也跑得稳。希望这篇文章能帮你少走一点弯路，把风电不确定性调度这个方向的研究和工程落地做得更顺手。