风电高比例接入之后调度部门最头疼的往往不是“风没有了怎么办”而是“风出力到底按多少算”。机组组合Unit CommitmentUC这个经典问题一旦叠加上风电的不确定性就从确定性优化变成了带随机参数的复杂决策。这个项目的标题很直白用线性准则下的分布鲁棒优化来解机组组合问题并给出Matlab实现。简单说就是在一个“不知道风电真实概率分布、只知道大概统计信息”的前提下让系统既能守住安全底线又不会因为过度保守而烧掉太多煤。这套思路对写论文的研究生、做电网调度算法开发的工程师、以及刚接触分布鲁棒优化DRO的算法爱好者都有参考价值。这篇文章我会把这套方法的建模思路、关键数学处理、Matlab代码框架以及我调试过程中踩过的坑全部摊开来讲。1. 项目概述与问题定位1.1 这个问题到底在解什么机组组合问题的本质是决定“哪些机组开、哪些机组关、每台机组发多少电”让整个系统在满足负荷和各类运行约束的前提下总成本最低。它是个标准的混合整数规划问题机组启停是0/1整数变量出力是连续变量再加上功率平衡、爬坡约束、最小启停时间约束组合爆炸的规模让求解从来都不轻松。风电一进来问题就变了味道。风电出力不是一个固定值而是围绕预测值上下波动的随机量。你今天定好的机组启停方案到了实际运行时刻风一调皮功率不平衡就出现了。传统的确定性UC只能针对一个预测场景做优化遇到风电出力大偏差就得靠旋转备用硬扛。于是就有了各种不确定性处理方法——随机优化、鲁棒优化、分布鲁棒优化。这个项目的核心就是在机组组合的框架里把风电的不确定性用“分布鲁棒”的方式刻画进去。它不假定风电误差服从某个精确的概率分布而是认为真实分布落在某个“模糊集”内然后基于最坏情况分布来优化决策。换句话说我不知道风到底怎么变但我给所有可能出现的概率分布画了一个圈站在圈里最恶劣的那个分布上做决策保证方案依然可行且经济性可接受。1.2 随机优化、鲁棒优化、分布鲁棒三者的取舍很多人刚接触这个概念时觉得DRO高深其实它是对随机规划和鲁棒优化的一次折中把两者的优缺点看明白你就知道为什么选择它。方法不确定性建模方式决策尺度问题随机优化假设误差服从已知分布通常用场景集近似目标为期望成本分布假设一旦错了方案就偏鲁棒优化误差落在确定的集合内但集合内不再谈概率目标为最坏情况成本过度保守成本虚高分布鲁棒优化真实分布属于一个模糊集只利用矩或距离信息目标为最坏分布下的期望成本分布信息少但决策更稳健我个人的体会是随机优化在风电建模里最大的软肋是“你根本拿不到真实分布”。风电预测误差的分布跟季节、天气过程、地形都有关拿历史数据拟合出来的分布换一个气象过程就失真。鲁棒优化则是另一个极端它只关心风电出力落在哪个区间完全不看概率结果就是为了一年可能只出现几次的极端场景把系统调度得非常保守成本高得离谱。分布鲁棒优化站在中间你不需要给出精确分布只需要给出分布应该满足的约束条件。比如均值大约是多少、协方差大概是什么范围或者真实分布离历史经验分布的距离不超过某个阈值。它把决策建立在“最坏但仍可能发生”的分布上既不像随机优化那样轻信假设又不像鲁棒优化那样完全无视概率。1.3 “线性准则”为什么是核心这个项目的副标题是“线性准则”很多人一眼看过去不知道在说什么。这里的线性准则指的就是两阶段分布鲁棒优化里的线性决策规则Linear Decision RuleLDR。两阶段模型的意思是第一天先决定机组启停和基准出力第一阶段第二天风电的实际出力偏差出现了调度员再根据偏差量做实时调整第二阶段。问题是第二阶段的调整量是风电偏差的函数这个函数是无穷维的——理论上任何映射关系都有可能。不处理这个无穷维问题模型根本没法求解。线性决策规则做了一件非常“工程化”的事情限制这个调整函数必须是偏差量的线性函数。也就是说实时调整量等于一个常数项加上一个线性系数矩阵乘以风电偏差向量。这样一来原本要在一堆函数空间里做优化的问题退化成了一堆实数系数的优化问题数学上立刻从“几乎不可解”变成“可以交给标准求解器解决”。更妙的是电力系统里这种做法其实有物理背景。调度员在实际运行中做AGC自动发电控制分配时用的本来就是参与因子法——按机组的调节能力把系统功率偏差按比例分摊给各台机组。这本质上就是线性决策规则。所以LDR不是纯数学上的妥协它在电力系统里有天然的工程合理性这也是这套方法不显得“空中楼阁”的重要原因。2. 分布鲁棒机组组合模型从思路到数学表达2.1 确定性机组组合的数学框架在加入风电不确定性之前先把确定性UC的底子打好。标准的UC模型通常长这样目标函数是燃料成本加上启停成本。燃料成本用二次函数近似经常要分段线性化这是后面用Matlab建模时第一个要处理的细节。启停成本则按启动方式和停机时长分类。功率平衡约束算是整套模型的灵魂所有开机机组的出力之和必须等于系统负荷不含风电时。这个约束在确定性模型里是等式在含风电模型里则要小心处理后面会细说。除此之外还有出力上下限约束、爬坡约束、最小启停时间约束、旋转备用约束。这些约束里面最小启停时间约束实现起来最麻烦需要引入额外的辅助变量或者用特殊的大M写法。2.2 风电不确定性与再调度机制引入风电之后确定性UC里的等式功率平衡不能直接用了因为风电出力在实时阶段是一个不确定量。处理的办法是引入再调度re-dispatch机制。假设风电预测出力是已知的实际出力是预测值加上误差ξ。在第一阶段调度员先制定机组启停计划和基准出力。到了第二阶段风电实际值揭晓如果系统出现了功率偏差调度员就要调整各台机组的出力来重新满足平衡。这里有一个很微妙的点启停变量在第一阶段就固定了第二阶段只能调整出力幅度不能再改变机组开停状态。这在物理上是合理的实际调度中不可能因为风电抖动就让机组频繁启停。所以第二阶段的调整量存在上下界变相限制了系统应对不确定性的能力。2.3 模糊集给概率分布“画圈”分布鲁棒优化的核心构件是模糊集Ambiguity Set。说白了就是把你认为“可能出现的概率分布”都装进一个集合里然后在这个集合里找最坏的那个分布来评估成本。目前主流有两种模糊集。第一种是基于矩信息的模糊集它约束真实分布的均值落在某个区间、协方差落在某个半定锥内。这类模糊集的优点是结构清楚对偶转化有现成的数学工具但缺点是它对分布形态的限制比较粗糙而且如果矩估计本身偏差大模型也会跟着偏。第二种是基于Wasserstein距离的模糊集它以历史数据的经验分布为中心画一个Wasserstein半径真实分布只要离经验分布不超过这个半径就算“合法”。这类模糊集的优势在于半径的选择与数据量有直接关系样本越多半径可以设得越小模型也越精确。近几年的电力系统DRO论文里Wasserstein模糊集几乎成了默认选项。模糊集半径怎么选是实操中的关键问题。理论上Wasserstein半径应当随样本量按某个速率衰减但实际中并没有一个统一公式能把半径跟经济性直接挂钩。我的经验是先按理论公式算一个基准值然后以5%为步长做几次灵敏度扫描画一条“成本-半径”曲线选在曲线斜率突变之前的那个半径值。这个做法工程上很实用比死记公式靠谱。2.4 两阶段分布鲁棒模型与线性决策规则的嵌入把前几节拼起来的完整模型本质上是一个min-max-min结构。最外层最小化总成本开机成本基准燃料成本实时调整成本中间层最大化也就是找最坏情况分布最内层最小化实时调整成本。三层结构在数学上是不能直接求解的必须层层处理。这里是整套算法最挑人的地方也是很多论文实现卡壳的根源。第一步处理最内层的最小化问题。对给定的一阶段决策和给定的风电偏差值最内层是一个线性规划因为再调度量对ξ的依赖是线性的线性决策规则风电偏差ξ本身是确定的。线性规划的极小值问题可以用对偶理论转成极大值问题或者直接跟外层的最大化问题合并。第二步引入线性决策规则。把第二阶段再调度量写成r(ξ) r0 R·ξ其中r0是常数向量R是系数矩阵。这里R矩阵的第(i, j)个元素表示第j台风电机组出力偏差对第i台火电机组再调度量的影响强度。把这个规则代入约束条件原本关于函数r(ξ)的约束就退化成了关于r0和R的矩阵不等式约束。第三步处理中间层的最大化问题。对矩模糊集通常用对偶变换和S-lemma消去概率测度变量对Wasserstein模糊集可以利用其强对偶性把最坏分布下的期望成本转化为一个包含半径惩罚项的确定性表达式。很多论文里一大段对偶推导实际落到代码里就是几个矩阵形式的约束替换。做完这三步原来的三层问题变成了一层一个包含决策变量x、r0、R的凸优化问题。如果原本的成本函数和约束都是线性的最终模型就是LP或者MILP用Gurobi就能跑如果涉及到半定约束那就要升级到Mosek这类SDP求解器。这一点直接决定Matlab代码里用哪个工具箱后面会细说。3. Matlab代码实现全流程3.1 环境准备与工具箱选型Matlab环境下做这类优化建模YALMIP是绕不开的工具。YALMIP是一个建模层你不用手写各种标准形式的矩阵而是用接近自然语言的语法定义变量和目标它会自动把模型翻译成求解器需要的标准格式。求解器方面要分两种情况。如果最终模型是MILP整数变量来自机组启停我推荐Gurobi或CPLEX性能在电力系统这种大规模MILP上明显优于MATLAB自带的intlinprog。如果最终模型含半定约束Mosek是SDP求解里的可靠选择。版本上提个醒YALMIP对Matlab新版本的支持偶尔会滞后我在R2023b上遇到过YALMIP的sdpvar在特定用法下报错的情况后来升级了YALMIP的GitHub版本就好了。遇到诡异报错优先检查YALMIP是不是最新版别一上来就怀疑模型写错了。3.2 测试系统与数据准备我去掉具体的算例细节用经典的6节点系统做测试是最省心的选择3台火电机组1个风电场3条负荷母线。机组参数按标准算例设置包括出力上下限、爬坡率、燃料成本系数以及单位启动成本。风电数据用预测出力和一组预测误差的历史样本。误差样本的生成要注意不能直接用正态分布随机数那等于自己给自己预设了分布假设违背了分布鲁棒的初衷。正确做法是用某个季节的历史实测误差曲线做抽样或者基于NWP数值天气预报的历史误差统计生成。3.3 YALMIP建模关键代码下面是整个模型最核心的代码逻辑简化版。第一步定义变量。% 机组数、风电场数、时段数 nG 3; nW 1; T 24; % 第一阶段变量 u binvar(nG, T, full); % 启停状态 p0 sdpvar(nG, T, full); % 基准出力 % 线性决策规则变量 r0 sdpvar(nG, T, full); % LDR常数项 R sdpvar(nG, T, full); % LDR线性系数简化版 % 这里R的下标含义第i台机组在第t时段对风电偏差的响应系数第二步构建目标函数。注意目标里有一部分是“最坏分布下的实时调整期望成本”在Wasserstein模糊集下经过对偶转换后这部分会变成一个包含半径惩罚项的表达式。代码里我会用一个函数封装这个惩罚项避免主程序里一长串公式影响可读性。% 目标函数开机成本 基准燃料成本 DRO调整成本 objective sum(sum(c_start .* max(0, diff(u, 1, 2)), 2)) ... sum(sum(fuel_cost(p0))) ... dro_adjustment_cost(r0, R, W_error_samples, radius);第三步加约束。功率平衡约束要特别注意在LDR框架下系统功率平衡必须对任意风电偏差ξ都成立所以它是一个“对任意ξ”的约束最终要转换成关于R的等式约束。constraints []; % 每个时段的功率平衡任意风电偏差下 for t 1:T % 火电出力 风电实际出力 负荷 % 风电实际出力 风电预测 偏差 % 火电出力 p0 r0 R * xi % 因此必须满足p0 r0 R*xi (Wforecast xi)*ones(nW,1) Load % 对任意xi成立需要常数项和系数项分别恒等 constraints [constraints, sum(p0(:,t)) sum(r0(:,t)) Wforecast(t) Load(t)]; constraints [constraints, sum(R(:,t)) 1 0]; % xi的系数项恒等于0 end % 出力上下限约束同样要保证任意偏差下都不越界 % 这类约束最终转化为 R 矩阵的线性不等式约束 常数项约束 for i 1:nG for t 1:T constraints [constraints, p0(i,t) r0(i,t) Pmin(i) * u(i,t)]; constraints [constraints, p0(i,t) r0(i,t) Pmax(i) * u(i,t)]; % 偏差上下的边界由 R 和模糊集半径共同决定 end end这里要解释一个新手很容易踩坑的点约束“对任意ξ成立”的处理方式。以功率平衡为例等式左边出现了一个含ξ的项要做的是把常数项和ξ的系数分开然后分别令它们满足等式。这样写出来的约束才是完整的否则你只是让某个特定ξ场景下平衡成立那整个分布鲁棒模型就名存实亡了。最后调用求解器ops sdpsettings(solver, gurobi, verbose, 2, showprogress, 1); optimize(constraints, objective, ops);求解完成后通过value()函数取出u、p0、r0、R的数值再算各项成本的明细。3.4 结果输出与对比跑完模型以后不要只看一个总成本数字就完事一定要做三组对比确定性UC、传统鲁棒UC、分布鲁棒UC。这三组对比的数字很有说服力。确定性UC的成本最低但只要你把真实风电偏差样本代进去做事后校验大概率会出现出力越界或者功率不平衡鲁棒UC永远安全但成本比确定性UC高出很多分布鲁棒UC处于两者之间——它不会出现越界只要模糊集画得合理成本也只比确定性UC高一点点。这个“经济性-鲁棒性”的平衡曲线就是你论文里最值得展示的核心结果。可视化方面我建议画两类图。第一类是24小时机组启停状态甘特图能直观看出DRO方案相比确定性问题哪些机组被提前开了、哪些时段留了更多备用。第二类是风电偏差-再调度量散点图把你解出来的R矩阵画成一条直线对比实际的最优再调度曲线能直观展示线性决策规则的逼近效果。4. 常见问题与排查技巧汇总4.1 求解时间爆炸怎么办机组组合本身是MILP加入分布鲁棒部分后变量数目进一步膨胀求解时间很容易从几分钟跳到几小时。我的排查路径是第一步看模型里有没有不必要的整数变量。比如机组启动变量能用3台机组的状态差分推导出来就没必要再单独建一个辅助整数变量。第二步检查是否可以用Benders分解把混合整数主问题只有u和p0和分布鲁棒子问题只有r0和R分开迭代求解。第三步如果坚持用单层模型务必打开求解器的MIP emphasis选项并给一个合理的MIP gap比如1%别试图求到绝对最优解工程上没那个必要。4.2 模型无解或数值异常YALMIP返回infeasible时第一步不是改代码而是用assign和check命令逐条检查约束的残差。我在实际项目中遇到过的情况是约束本身没有问题但变量量纲差别太大——出力是百兆瓦级成本系数是$/MWh两者相乘之后数值上了1e6而偏差量的数值只有0.1量级导致矩阵条件数极差求解器直接判定不可行。解决办法是把所有变量归一化到同一量纲或对目标函数整体除以一个基准值再优化。还有一个隐蔽问题线性决策规则里的系数矩阵R如果某一行的所有元素都是0说明这台机组完全不参与再调度。这在数学上合法但物理上意味着这台机组只是个“纯基荷”如果它的出力上下限比较宽实际运行中可能反而造成调节能力浪费。这种结果通常不是模型bug而是机组成本参数设置不合理导致的需要回到数据层面调整。4.3 模糊集半径怎么调Wasserstein半径的选择直接影响模型的经济性。半径设得太大DRO就退化成鲁棒优化成本虚高半径设得太小DRO跟随机优化差不多抵御风险的能力不足。一个可操作的经验做法是用历史数据做留一法交叉验证。把样本切成训练集和测试集训练集用来构造经验分布和模糊集测试集用来仿真决策方案的越界频率。不同半径下的越界频率画出来通常是一条单调下降的曲线选越界频率低于某个阈值比如0.1%且半径尽量小的点即可。这个方法工作量不大但有说服力。4.4 对偶转换中的强对偶条件分布鲁棒模型里最脏最累的活是对偶而这个环节恰恰最容易出错。Wasserstein模糊集的对偶化依赖强对偶条件成立而强对偶条件成立需要内点条件等前提。很多论文里轻描淡写一句“by strong duality”实际代码里你就会发现转换后的目标函数符号反了、约束方向反了、惩罚项系数差了倍数。我的做法是先用一个只包含单个机组、单时段退化到手工可算的小例子在Matlab里逐个数验证对偶公式的常数项、一阶项、惩罚项分别对不对。这个小例子花半天时间但能省下后面一星期的Debug时间。5. 个人实操体会机组组合问题做了这么多年最深的感受是模型能跑通不等于方案能落地。分布鲁棒优化给了我们一个很优雅的数学框架但它的实用性很大程度上取决于你对模糊集的选择和对数据的理解。线性决策规则看起来是个“近似”但实际上跟电力系统的AGC参与因子分配机制高度契合所以它不只是数值上的技巧更是工程直觉的数学化表达。如果你打算把这套代码扩展到更大规模的系统我建议优先考虑三个方面一是把风电的多时段相关性加进去目前很多模型都假设各时段偏差独立这跟实际情况出入很大二是把储能系统纳入再调度变量储能是天然的线性决策规则执行器两者结合会非常自然三是用真实的风电历史数据替换掉模型里的正态假设数据质量往往比模型复杂度更影响最终方案的可信度。最后分享一个小技巧在所有结果展示里加一张“确定性方案在真实偏差下的平衡越界图”这个图不需要任何额外数据就是把你求出的确定性UC方案拿去做事后仿真把每个时段的平衡偏差画出来。往往就是这张图最直观地说明了“为什么我们非要用分布鲁棒不可”——理论讲得再好看也不如一张波动曲线有说服力。