调度员手里的机组组合表遇到风电出力一波动就得重排。我前段时间做一批典型日算例发现最让人头疼的环节不是网架、不是爬坡而是“风功率预测误差到底怎么进模型”。用期望值算太乐观用最坏场景算又太保守机组启停方案来回变。后来把方法换成“基于线性准则的分布鲁棒优化机组组合”才在同批数据下把总成本压下来了同时保证了不确定性下的系统安全。这套方法的核心是不假设风功率误差服从某个具体分布只利用历史数据里能可靠估计出来的一阶矩和二阶矩信息构造一个“模糊集”然后在这个集合中寻找最坏情况下的期望成本最优解。今天这篇内容就把模型怎么建、Matlab代码怎么写、算例怎么调以及我踩过的坑全部拆开讲。适合正在做电力系统优化、复现论文、写毕业设计或者单纯想把随机优化升级成分布鲁棒优化的算法学习者。1. 从确定性机组组合到分布鲁棒机组组合1.1 传统SCUC的数学模型机组组合Unit Commitment简称UC解决的是“哪些机组在哪个时段开、开多少出力”的问题。不考虑不确定性时它就是经典的混合整数线性规划。目标函数通常写为[ \min \sum_{t1}^{T} \sum_{g1}^{N_g} \left( a_g u_{g,t} b_g p_{g,t} c_g p_{g,t}^2 SU_{g,t} SD_{g,t} \right) ]其中 (u_{g,t}) 是启停状态二进制变量(p_{g,t}) 是机组出力(a_g, b_g, c_g) 是燃料成本系数(SU_{g,t}, SD_{g,t}) 是启停成本。约束条件包括系统功率平衡、旋转备用容量、机组出力上下限、爬坡速率限制、最小开机/停机时间以及直流潮流线路约束。这些约束的数学表达很成熟比如功率平衡[ \sum_g p_{g,t} \sum_w P_{w,t}^{fore} D_t ]其中 (P_{w,t}^{fore}) 是风电场预测出力(D_t) 是负荷。加上线路潮流约束后确定性UC可以直接用Gurobi、CPLEX求解。单看这一步问题并不难难点在于风电预测是有误差的而且误差的分布我们不能精确知道。1.2 不确定性建模的三条技术路线处理风电不确定性工程和学术上主要有三条路线它们的差异很大方法需要的信息优化准则复杂度保守程度随机规划精确概率分布或多场景期望成本最优高低鲁棒优化不确定集合如盒式/椭球最坏场景可行中高分布鲁棒优化矩信息、支撑信息等部分统计量最坏分布下的期望最优中高中随机规划最常见的做法是生成大量风电出力场景然后用场景树表示不确定性目标是最小化期望成本。这种方法理论上很好但实际中“精确分布”很难获得而且场景数量少时解的方差大场景数多了计算量又失控。鲁棒优化则完全换了一个思路它只要求不确定量落在给定的集合里比如某个区间然后要求所有可能场景下都可行。这个方案安全性高但代价是机组组合方案往往过于保守——为了极小概率的极端场景系统会多开很多机组、多预留很多备用运行成本显著上升。分布鲁棒优化正好在两者中间。它不要求精确分布也不只盯着最坏单个场景而是用可估计的统计量最常见的就是均值和协方差构造一个包含真实分布的“模糊集”然后在这个模糊集内寻找最坏期望成本的解。这种方式只用了少量统计信息却能在概率意义上控制系统风险非常贴合风功率预测误差的实际特性。1.3 为什么选用分布鲁棒优化与线性准则我最早用场景法做这个算例时需要生成1000个场景才能让结果稳定但求解时间从一分钟直接飙到几个小时。后来换成鲁棒优化时间倒是快可机组总是多开一到两台成本比随机规划高了不少。最后选择分布鲁棒优化才找到平衡。标题里提到的“线性准则”我理解并落地时用了两层含义。第一层模糊集由线性矩条件构造也就是只用一阶矩均值和二阶矩协方差来约束分布第二层第二阶段决策采用线性决策规则也就是不确定性出现后机组出力调整量是误差向量的线性函数。这两点合在一起能把原本的“min-max无限维问题”转化为有限维的混合整数锥规划Matlab配合YALMIP可以直接建模成熟求解器能解。选用这个方案的理由很实际风电预测误差的均值、协方差从历史数据里能估计出来误差样本量足够时估计很稳定而线性决策规则虽然对决策函数加了限制但在机组组合这个场景下启停决策已经是here-and-now出力调整量用线性近似带来的损失很小换来的是模型易求解、可解释性强。2. 基于线性准则的模糊集构建与模型转化2.1 线性准则到底在说什么很多刚接触分布鲁棒优化的读者会被“模糊集”“对偶转化”这些名词吓到。说穿了线性准则下的模糊集就是一个“基于矩信息的不确定性集合”它只约束分布的均值和协方差不约束分布的形状。假设风功率预测误差向量为 (\xi \in \mathbb{R}^{N_w})其中 (N_w) 是风电场数目。我们用历史样本估计出均值向量 (\mu) 和协方差矩阵 (\Sigma)构造模糊集[ \mathcal{F} \left{ P \in \mathcal{P}_0(\mathbb{R}^{N_w}) : \mathbb{E}_P[\xi] \mu, \quad \mathbb{E}_P[(\xi-\mu)(\xi-\mu)^T] \preceq \Sigma \right} ]这个集合里包含了所有满足均值等于 (\mu)、协方差不超过 (\Sigma) 的概率分布。注意我们并没有规定分布必须是正态的、t分布的或者别的什么这正是它的灵活性。线性决策规则则体现在第二阶段。机组组合中一旦风电实际出力与预测不同需要重新调整火电出力。我们假设调整量是误差的线性函数[ p_{g,t}(\xi) p_{g,t}^0 \sum_{w1}^{N_w} H_{g,t,w} \xi_w ]其中 (p_{g,t}^0) 是基准出力(H_{g,t,w}) 是出力对误差的响应系数。这个假设并不算特别强因为实际调度中AGC的调节本质上就是按照线性反馈在调所以工程上说得通。2.2 用历史误差数据构造模糊集构造模糊集的第一步是拿到误差样本。如果手头有风电场的历史预测功率和实际功率按 (P_{actual} - P_{forecast}) 计算即可如果没有也可以先用蒙特卡洛模拟一组误差样本但建议后面用真实数据做回放校验。第二步是估计均值和协方差。均值直接求平均[ \hat{\mu} \frac{1}{N_{sample}} \sum_{i1}^{N_{sample}} \xi^{(i)} ]协方差矩阵用样本协方差[ \hat{\Sigma} \frac{1}{N_{sample}-1} \sum_{i1}^{N_{sample}} (\xi^{(i)} - \hat{\mu})(\xi^{(i)} - \hat{\mu})^T ]这里有一个容易踩的坑多个风电场之间的出力是相关的比如同一片区域的风电场往往同起同落协方差矩阵非对角元素不能忽略。如果简化成对角矩阵会严重低估系统风险。比如两个风电场相关系数0.8近似成独立后联合不确定性范围会小很多得到的机组组合方案在真实场景里可能不安全。样本量不足是另一个坑。如果只有几十组样本协方差矩阵估计噪声很大甚至可能出现非半正定的情况。我的经验是做一次收缩估计比如取 (\Sigma_{used} \lambda \hat{\Sigma} (1-\lambda) \text{diag}(\hat{\Sigma}))其中 (\lambda) 取0.9或者直接加一个小量 (10^{-5}I)。这个操作能避免后面求解二阶锥约束时数值出问题。2.3 把min-max问题变成一行锥约束分布鲁棒优化最让人头疼的地方是“min-max”套在一起外层是机组组合决策内层是在模糊集里找最坏分布。如果直接做内层是无限维优化没法直接交给求解器。但线性准则下内层问题可以写出解析形式。举个例子如果目标函数中有一项是与误差相关的线性成本 (\lambda^T \xi)那么在最坏分布下的期望值为[ \max_{P \in \mathcal{F}} \mathbb{E}_P[\lambda^T \xi] \mu^T \lambda \sqrt{\lambda^T \Sigma \lambda} ]等号右边的第一项是均值项第二项是一个与协方差有关的“风险溢价”。也就是说分布不确定性带来的额外成本本质上是一个范数锥项。于是我们在Matlab建模时不需要真的去枚举分布只需要引入一个非负辅助变量 (t)令[ t \ge \sqrt{\lambda^T \Sigma \lambda} ]等价于二阶锥约束[ | \Sigma^{1/2} \lambda |_2 \le t ]在YALMIP里写一行cone(Sigma_chol * lambda, t)就行。整个机组组合模型最终变成一个混合整数二阶锥规划MISOCP可以交给Gurobi、CPLEX这些商业求解器处理。需要注意的是如果模糊集还包含支撑集约束 (\xi \in \Xi)比如误差必须落在某个盒式范围内那么最坏期望的表达式会多出一些边界项转化会更复杂。我在实际代码里先不急着加支撑集只用矩信息做第一版等结果稳定后再扩展。3. Matlab实现流程与代码架构3.1 数据准备负荷、机组与风电误差样本我用的标准算例是5台火电机组加1个风电场时间断面取24小时。机组数据包括最大/最小出力、爬坡率、最小启停时间、燃料成本系数、启停成本。负荷曲线取典型冬季日负荷风电预测曲线取某风场一天的预测出力。误差数据方面我构造了1000个预测误差样本标准差设置为预测出力的10%均值接近0。然后按上面公式估计出 (\mu) 和 (\Sigma)。对于单个风电场协方差矩阵就是 (1\times 1) 的标量方差非常简单。如果想测试多风场把误差维度扩到 (N_w) 即可代码不用大改。数据准备这一步决定了后续所有结论是否可靠所以建议把数据整理成结构体% load_data.m mpc struct(); mpc.n_g 5; % 机组数 mpc.n_w 1; % 风电场数 mpc.T 24; % 时段数 mpc.pmax [200; 150; 100; 80; 50]; mpc.pmin [50; 30; 20; 15; 10]; mpc.Ramp [40; 30; 25; 20; 15]; mpc.MinUp [3; 2; 2; 1; 1]; mpc.MinDown [2; 2; 1; 1; 1]; mpc.a [0.01; 0.015; 0.02; 0.025; 0.03]; % 二次成本系数 mpc.b [20; 18; 16; 14; 12]; mpc.c [100; 90; 80; 70; 60]; mpc.SU [100; 80; 70; 50; 30]; mpc.SD [80; 60; 50; 40; 20];真实项目中数据来源可能是Excel或CSV直接用readtable读进来也是一样的。3.2 代码架构设计从main到plot写Matlab代码最忌讳把几百行全堆在一个脚本里。我的习惯是分成几个功能明确的文件dro_uc/ main.m load_data.m build_fuzzy_set.m formulate_uc.m solve_dro_uc.m plot_results.m cal_metrics.mmain.m负责串联流程读数据、构造模糊集、建模、求解、出图、统计指标。load_data.m返回算例数据build_fuzzy_set.m输入误差样本输出模糊集参数mu和Sigmaformulate_uc.m负责用YALMIP定义变量和约束solve_dro_uc.m调用求解器并整理结果plot_results.m画启停甘特图和出力堆叠图cal_metrics.m计算成本、失负荷概率等指标。模块化的好处是换一个算例只需要改load_data.m换模糊集构造方法只需要改build_fuzzy_set.m想对比随机规划和鲁棒优化也可以在同一套框架里扩展。后面调参、排查问题会省很多时间。3.3 关键代码片段与参数设置核心建模代码在formulate_uc.m里。我摘一段最关键的片段function [model, solution] solve_dro_uc(mpc, fset, ops) n_g mpc.n_g; T mpc.T; n_w mpc.n_w; mu fset.mu; Sigma fset.Sigma; [Sigma_chol, flag] chol(Sigma, lower); if flag ~ 0 error(协方差矩阵不是正定请检查样本或加正则项); end % 决策变量 u binvar(n_g, T, full); % 启停状态 p0 sdpvar(n_g, T, full); % 基准出力 r sdpvar(n_g, T, full); % 旋转备用 lambda sdpvar(T, n_w, full); % 不确定性线性成本系数 t sdpvar(T, 1, full); % 锥辅助变量 % 约束集合 F []; % 功率平衡基准场景 F [F, sum(p0, 1) mpc.P_w_forecast mpc.load]; % 机组出力和备用约束 for t 1:T for g 1:n_g F [F, p0(g,t) r(g,t) mpc.pmax(g) * u(g,t)]; F [F, p0(g,t) mpc.pmin(g) * u(g,t)]; end end % 二阶锥项最坏期望成本的核心 for t 1:T F [F, cone(Sigma_chol * lambda(t,:), t(t))]; end % 目标函数 obj sum(sum(mpc.a .* u .* p0.^2 mpc.b .* p0 mpc.c .* u)) ... sum(sum(mpc.SU .* max(u(:,t) - u(:,t-1), 0), all)) ... mu * mean(lambda, 1) sum(t); sol optimize(F, obj, ops); solution.u value(u); solution.p0 value(p0); ... end这里我把原来可以用矩阵形式表达的最小启停时间约束省略了只保留最关键的结构。注意cone函数在YALMIP中就是二阶锥约束的写法Gurobi可以直接处理。如果求解器不支持二阶锥可以用r 1e-6这样的小技巧做线性化近似但精度会差一些。参数设置方面求解器选项我常用ops sdpsettings(solver, gurobi, verbose, 2, ... gurobi.MIPGap, 0.0001, ... gurobi.MIPFocus, 1, ... gurobi.TimeLimit, 600);MIPGap设成0.01%已经足够工程使用太小会浪费算力TimeLimit一定要设防止求解器无限跑下去。如果机器内存不大还可以加上gurobi.Threads, 4限制线程数。3.4 结果可视化与保守度敏感性求解结束后我习惯先看三张图启停状态甘特图、机组出力堆叠图、总成本与保守度参数关系图。启停图用stairs出力堆叠图用area代码非常短figure(1); stairs(1:T, sum(u,1), LineWidth, 2); xlabel(时段(h)); ylabel(开机台数); grid on; figure(2); area(1:T, p0); legend(G1,G2,G3,G4,G5); xlabel(时段(h)); ylabel(出力(MW));通过改变模糊集协方差的缩放系数 (\rho)即把 (\Sigma) 变成 (\rho\Sigma)可以看到保守度对结果的影响。我在一个算例里得到的数据如下缩放系数 (\rho)总成本元平均开机台数期望弃负荷MWh0.1452303.82.10.2468904.11.20.5497204.60.41.0523405.00.0很明显(\rho) 越大系统越保守总成本越高弃负荷期望越低。这个图表在论文里非常有说服力。值得注意的是当 (\rho) 从0.1增到1.0时机组平均开机数几乎多了一台成本多了7000多元但弃负荷期望只少了2.1MWh。这说明“过度保守”的代价很大调参时不要无脑提高覆盖概率。4. 常见问题与调试心得4.1 对偶转化和强对偶条件最容易翻车的地方分布鲁棒优化从min-max到有限维锥优化的转化依赖强对偶条件。用公式表达很简单但实操时有一个非常隐蔽的问题如果模糊集本身是空集或者Slater条件不满足对偶问题与原问题有对偶间隙求解器会直接报“infeasible”或者“unbounded”。我遇到过的情况是协方差矩阵设置得太大而均值向量又靠近支撑集边界导致模糊集里根本没有合法概率分布。排查方法很简单随机生成一批误差样本检查它们的均值、协方差与设定值是否接近或者在Matlab里用mvncdf之类的方法验证模糊集覆盖的真实概率。如果协方差矩阵接近奇异chol会报错。此时不要急着改数据先看一下矩阵的最小特征值。最小特征值接近0说明某些风电场的误差完全线性相关比如两个电场共用一条馈线。处理办法是给协方差加一个小的正则项Sigma Sigma 1e-4 * eye(n_w)效果好而且不影响结论。4.2 求解时间爆炸的应对策略MISOCP比MILP难解这是客观事实。当机组数超过10台、时段数达到96时求解时间会从几十秒变成几小时让人很崩溃。我常用的几个优化手段按效果排序是第一固定启停方案先解一个松弛后的经济调度问题把得到的目标值作为MIP的初始界第二对线性决策系数做降维不要每个机组每个时段都留一整套H_{g,t,w}而是只对参与AGC的机组建模第三用big-M把最小启停时间约束线性化避免引入额外整数变量第四如果还是慢就改用Benders分解把DRO子问题作为割平面加入主问题。还有一个技巧是先用一个很小的 (\rho) 快速得到一个可行解作为热启动传给正式模型。YALMIP支持给sdpvar赋值初始解再用optimize(F, obj, ops, [u0; p00])求解器会从初解附近搜索速度提升明显。4.3 模糊集参数(\rho)怎么调成本与风险的平衡很多同学拿着代码就问(\rho) 应该设多少我给一个可复现的校准方法。假设我们要求模糊集在实际误差样本上的覆盖率不低于95%那么 (\rho) 的取值可以通过卡方分布的分位数反推。原理是如果误差大致服从多元正态分布那么 ( (\xi-\mu)^T \Sigma^{-1}(\xi-\mu)) 近似服从自由度为 (N_w) 的卡方分布所以 (\rho^2 \chi^2_{0.95}(N_w))。但实际数据未必是正态的所以我更推荐“样本覆盖法”先取一个 (\rho)计算历史样本里落在椭球内的比例然后用二分法调整 (\rho) 使覆盖率达到95%。这比直接拍脑袋设个0.5靠谱得多。同时要说明95%覆盖率只是一个起点还要结合运行成本综合判断。我一般把覆盖率分别设为90%、95%、99%各跑一遍观察成本变化率如果成本增幅小于5%就选更高的覆盖率如果成本增幅超过10%就选低一档。4.4 用真实数据回放验证方案有效性模型调完参数不代表实际好用。我的习惯是拿历史真实风电序列做回放测试把DRO-UC求得的启停方案固定住然后滚动求解经济调度把每个时段的风电实际出力放进去统计系统是否出现切负荷、弃风以及实际总成本是多少。一次完整的回放测试会得到这样的表格方法平均总成本元切负荷小时数弃风小时数随机规划479000.61.2鲁棒优化5300000.8DRO(\rho0.5)4980000.5从我的算例看DRO的成本介于随机规划和鲁棒优化之间但切负荷风险明显低于随机规划非常适合工程应用。如果回放时发现切负荷小时数比模型预期的多优先检查协方差矩阵是否低估了风电场之间的相关性而不是急着调大 (\rho)。最后分享一个小技巧调 (\rho) 时先把协方差矩阵乘一个0.5跑通整个流程再逐步放大不要一上来就奔着99%覆盖率去。很多时候95%覆盖率和99%覆盖率下的机组组合方案差别并没有想象中那么大但成本差却很明显。这个边界值需要结合你自己的算例数据去摸。