做机组组合调度的人大概都体会过这种纠结日前排计划的时候风电出力的预测曲线明明已经做得很细了可一到实际运行预测误差照样会“打脸”。偏差小的时候还行偏差一旦到了风电装机容量的百分之十几火电机组的爬坡和启停方案就完全不够看了。这个时候再用传统的确定性机组组合模型调度结果往往不是导致备用不足就是导致经济性太差两头吃亏。我近几年一直在做新能源不确定性下的机组组合问题试过随机规划也试过区间鲁棒优化各有各的短板。后来把方向转向了分布鲁棒优化Distributionally Robust OptimizationDRO并且用线性决策规则Linear Decision RuleLDR来处理风电不确定性整体效果让我比较满意。这篇就当是项目复盘把建模思路、数学推导、Matlab实现要点和踩坑记录一起写出来。内容偏技术适合正在做电力系统调度、研究鲁棒优化或者准备复现同类论文的同学参考。1. 为什么机组组合里的风电不确定性不能继续用传统手段处理1.1 从点预测到区间风电预测误差的真实形态风电预测误差不是简单的“均值附近对称分布”。实际统计下来它往往带有明显的厚尾和偏态特征。尤其在大风天气过程切换、季节交替、阵风时段误差分布形态会显著变化甚至出现双峰。拿过去某天某风电场的数据来看预测误差有时候能摸到风电装机容量的20%以上而且这种极端样本不是小概率事件一年里能碰上好几回。在机组组合模型中如果只用一个确定性的点预测值比如预测风电有功出力为300 MW那么实际上等于默认“预测即实际”没有任何安全裕度这在风电占比不高的时候还能靠系统热备用硬扛起来。但风电占比一旦提升点预测方案在极端场景下往往直接导致切负荷。1.2 随机规划与鲁棒优化的两端局限随机规划的经典做法是假设风电预测误差服从某个已知概率分布比如正态分布然后生成大量场景做随机机组组合。问题在于真实场景分布我们根本不知道正态分布假设只是数学上的便利。样本有限、场景生成偏差、概率密度估计误差都会让“最优解”变成“次优解”。而传统鲁棒优化则走另一个极端它把不确定性描述成一个确定性的有界集合比如让风电出力在每个时段都落在 ([P_{min}, P_{max}]) 的区间里然后在这个区间内找最坏情况下的可行解。听起来很稳妥但其保守性非常强。尤其是当需要同时考虑96个时段的风电出力不确定性时区间鲁棒模型会让每个时段都按照最极端的情况去预留调整能力最后调度出来的成本可能比随机规划高出一大截。1.3 分布鲁棒的现实定位我只知道一阶和二阶统计量分布鲁棒优化其实站在了两者中间。它的思路是我不需要知道风电预测误差的完整概率分布也不需要强迫自己相信“误差服从正态分布”这种过于理想的假设。我只需要从历史数据中估计出一阶矩均值和二阶矩协方差矩阵然后构造一个包含真实分布可能位置的模糊集ambiguity set在这个模糊集内做最坏期望优化。用大白话说就是既然我搞不清楚风电预测误差的真实分布那我把所有均值、协方差都符合历史数据的分布都看成“可能的对手”调度方案必须对这群分布里最坏的那个分布仍然可行、仍然经济。这才是标题里“分布鲁棒”四个字的真正含义。下表是我在模型选型时做的简单对比方便你理解几种方法的差异方法对分布信息的依赖保守程度求解难度适用场景确定性UC完全不考虑最低低风电占比低、预测准随机规划需要精确分布较低高场景多有大量可靠历史场景区间鲁棒只需上下界很高中重视安全但成本较敏感分布鲁棒只需均值、协方差适中、可调较高风电占比高的日前调度2. 两阶段DRO-UC模型风电场不确定集和决策变量的分工2.1 第一阶段与第二阶段各自承担什么责任机组组合本身天然适合两阶段建模。第一阶段是“日前决策”在还没看到实际风电出力之前就要决定哪些机组开机、哪些机组停机、各机组日前计划出力是多少。第二阶段是“实时调整”当风电实际出力确定下来之后调度员会在日前计划基础上对各机组进行再调度用爬坡能力去弥补风电偏差同时决定是否弃风、是否存在切负荷。在本文的分布鲁棒框架下第一阶段的启停变量是二元的这个没法含糊。第二阶段再调度量是连续变量而且它必须随预测误差呈线性变化。之所以引入这种依赖关系是因为实时调度不可能等不确定性实现之后再去重新求解一个大规模优化问题调度策略必须预先给定而且要足够简单到可以快速执行。2.2 用预测误差作为不确定性向量设风电预测误差向量为 (\xi)它的维数等于参与调度的风电场个数乘以调度时段数。每个时段的风电预测误差用 (\xi_t) 表示实际风电可用出力就可以写成[ P_{wind,t}^{actual} P_{wind,t}^{forecast} \xi_t ]这里的 (\xi_t) 不是区间鲁棒里那个简单的上下界变量而是带有统计约束的随机向量。我们构造的支持域 (\Xi) 用来描述 (\xi_t) 的物理边界比如[ \Xi {\xi_t : \xi_t^{min} \le \xi_t \le \xi_t^{max}} ]同时我们还构造一个模糊集[ \mathcal{D} \left{ \mathbb{P} : \mathbb{E}[\xi] \mu,\ \mathbb{E}[\xi \xi^\top] \preceq \Sigma,\ \mathbb{P}(\xi \in \Xi) 1 \right} ]其中 (\mu) 是历史预测误差样本均值(\Sigma) 是样本协方差矩阵。注意二阶矩约束这里用的是半定小于等于关系含义是真实分布的“能量”不能超过历史样本给出的二阶统计量这样构造出的模糊集既包含分布的位置信息又保留了尾部风险。2.3 两阶段决策变量设计为了后面写代码方便我们把关键决策变量列成一张表。第一阶段变量是常规的机组组合变量第二阶段变量全部写成预测误差 (\xi) 的仿射函数。阶段变量类型含义第一阶段(u_{i,t} \in {0,1})机组i在时段t的启停状态第一阶段(p_{i,t}^{sch})机组i在时段t的日前计划出力第二阶段(\Delta p_{i,t}(\xi))机组i在时段t的再调度出力第二阶段(\Delta w_{k,t}^{cur}(\xi))风电场k在时段t的弃风量第二阶段(l_{t}(\xi))时段t的失负荷量3. 线性决策规则如何把“看风调电”变成一行仿射函数3.1 线性准则的核心定义“线性准则”在本文里指的就是线性决策规则。它的思想非常直接把第二阶段变量写成不确定性向量 (\xi) 的仿射函数。以机组再调度出力为例[ \Delta p_{i,t}(\xi) \alpha_{i,t} \beta_{i,t}^\top \xi ]其中 (\alpha_{i,t}) 是常数项(\beta_{i,t}) 是线性系数向量二者都是待优化变量。类似地弃风量和失负荷量也可以表示成[ \Delta w_{k,t}^{cur}(\xi) \eta_{k,t} \theta_{k,t}^\top \xi ][ l_t(\xi) \gamma_t \phi_t^\top \xi ]这样处理之后实时调度策略就完全被有限个系数所刻画不再需要“等风来再算优化问题”只需要在实际误差 (\xi) 实现之后代入这几个公式即可。3.2 为什么线性规则能在两阶段问题里打开局面两阶段随机鲁棒优化的最大难点在于第二阶段的可行域和最优值都依赖于不确定性实现直接求解是一个无穷维优化问题。如果限制第二阶段变量与 (\xi) 成仿射关系那么对于任意给定的第一阶段决策变量和线性规则系数第二阶段所有约束都可以写成关于 (\xi) 的仿射不等式。接下来只需要保证这些仿射不等式在整个支持域 (\Xi) 上都成立即可。于是“看风调电”的无穷维问题就被压缩成了一个有限维的线性规划/半定规划问题。这个压缩是有代价的因为仿射策略并不等于全局最优策略它相当于在“策略空间”里限定了一个子空间。不过工程实践表明在机组组合这类问题里线性决策规则的性能损失很小却能换来巨大的计算优势。3.3 平衡约束与容量约束如何变成等式和不等式组举个例子系统功率平衡约束可以写成[ \sum_{i} \left( p_{i,t}^{sch} \Delta p_{i,t}(\xi) \right)\sum_k \left( P_{k,t}^{forecast} \xi_{k,t} - \Delta w_{k,t}^{cur}(\xi) \right)l_t(\xi) D_t ]把 (\Delta p)、(\Delta w)、(l_t) 全部代入仿射表达式后上式会化简成[ C_0 C_1^\top \xi 0, \quad \forall \xi \in \Xi ]由于 (\xi) 在支持域内是可以变化的要让这个等式恒成立唯一的办法就是常数项 (C_0) 和线性系数向量 (C_1) 都等于零。这个化简是很关键的它直接让随机约束变成有限个确定性约束。机组容量约束和爬坡约束同理。例如容量约束[ p_i^{min} u_{i,t} \le p_{i,t}^{sch} \alpha_{i,t} \beta_{i,t}^\top \xi \le p_i^{max} u_{i,t}, \quad \forall \xi \in \Xi ]由于左边和右边都是 (\xi) 的仿射函数只要在支持域的所有端点处验证不等式即可。当 (\Xi) 是一个盒式集合时端点个数是有限的因此可以枚举端点也可以利用对偶理论把它转化为更紧的半定约束。实际代码里我用的是YALMIP内置的鲁棒约束处理再把端点枚举作为交叉验证手段。4. 从模型到MatlabYALMIP实现的目标函数、模糊集与LDR约束4.1 目标函数里为什么一定要保留二次项很多初学者第一次看到DRO模型会问如果我目标函数里全是 (\xi) 的线性项那期望最坏值不就直接等于均值吗分布鲁棒唯一性体现在哪里这是一个非常关键的技术点。答案是只有当第二阶段成本或者惩罚项中含有 (\xi) 的二次项时模糊集中的二阶矩约束才真正起作用。在实际机组组合中火电成本往往用二次函数来拟合比如[ Cost(p) a p^2 b p c ]把 (p p^{sch} \alpha \beta^\top \xi) 代进去就会出现 ((\beta^\top \xi)^2) 这样的项。对这样的项取期望得到的是 (\beta^\top \mathbb{E}[\xi \xi^\top] \beta)它依赖概率分布的二阶矩。只有在此时模糊集 (\mathcal{D}) 的 (\mathbb{E}[\xi \xi^\top] \preceq \Sigma) 才真正约束了优化问题的方向。所以我在建模时目标函数写为[ \min_{u,p^{sch},\alpha,\beta,\cdots} \left[ \sum_{i,t} \left( Su_{i,t} a_i (p_{i,t}^{sch})^2 b_i p_{i,t}^{sch} \right)\sup_{\mathbb{P} \in \mathcal{D}} \mathbb{E}{\mathbb{P}} \left[ \sum{i,t} a_i \left( \alpha_{i,t} \beta_{i,t}^\top \xi \right)^2V_{l} l_t(\xi) V_{w} \Delta w_{k,t}^{cur}(\xi) \right] \right] ]其中 (V_l) 是失负荷惩罚价格(V_w) 是弃风惩罚价格。最外层 (\sup) 的含义是在模糊集里所有可能的分布中取最坏的那个来计算期望成本调度方案必须为这个最坏期望做优化。4.2 主脚本框架数据、变量、约束、求解Matlab实现的整体流程如下% 1. 载入系统参数机组、负荷、风电历史误差 load(uc_data.mat); % 2. 定义YALMIP变量 u binvar(nGen, T, full); % 机组启停状态 p_sch sdpvar(nGen, T, full); % 日前计划出力 alpha sdpvar(nGen, T, full); % 再调度常数项 beta sdpvar(nGen, nWind, T, full);% 再调度线性系数 % 3. 定义不确定性变量 xi sdpvar(nWind, T, full); uncertain(xi); % 4. 支持域和模糊集 W [xi xi_min, xi xi_max]; % 物理边界 Fmu mean(xi) mu; % 一阶矩约束 Fsig xi*xi Sigma; % 二阶矩约束的SDP表达这里特别注意uncertain(xi)声明的作用。在YALMIP中一旦某个变量被声明为不确定变量那么其他确定性变量与它相乘或参与约束时YALMIP会自动进行鲁棒对等转换。对于线性决策规则我会把再调度量显式写成p_sch alpha beta*xi_vec的形式而不是直接定义一个带不确定性的变量。完整的主脚本可以分为七个模块% 模块1数据准备 % 模块2变量定义 % 模块3目标函数 objective sum(sum(start_cost)) sum(sum(fuel_cost)) worst_case_cost; % 模块4约束条件 Constraints []; % 最小启停时间、爬坡率、功率平衡、容量约束 % 模块5调用求解器 options sdpsettings(solver, mosek, verbose, 2); optimize(Constraints, objective, options); % 模块6结果提取与后处理 % 模块7蒙特卡洛验证4.3 模糊集内部的最坏期望转成半定约束对于包含 (\sup_{\mathbb{P} \in \mathcal{D}} \mathbb{E}[\xi^\top Q \xi q^\top \xi c]) 的子问题我利用矩模糊集对偶理论把它转化为一个半定约束。具体来说引入对偶变量 (S)、(\lambda_0)、(\lambda)可以得到一个形如[ t \ge \lambda_0 \mu^\top \lambda \langle \Sigma, S \rangle ]的辅助约束同时要求 (S \succeq 0)并且对所有 (\xi \in \Xi) 满足[ \lambda_0 \lambda^\top \xi \xi^\top S \xi \ge \xi^\top Q \xi q^\top \xi c ]这一组约束可以通过支持域 (\Xi) 的端点枚举来检验也可以用半无限规划的对偶再处理一次。实践中我倾向于把 (\Xi) 取成各时段独立的盒式集合这样端点枚举的代价是可控的。5. 数值实验结果调度成本、失负荷风险与迭代求解的观察5.1 测试系统与风电场参数设置我在复现时用了一个典型的小型测试系统6台火电机组24个调度时段负荷峰值为560 MW接入两座风电场总装机容量为200 MW。风电预测误差的历史样本按以下规则生成均值取0标准差取风电出力预测值的8%误差上下界取装机容量的±20%。为了验证模型对不同风况的适应性我分别做了三组实验小风电渗透率、中风电渗透率、高风电渗透率。渗透率通过改变负荷曲线和风电场出力基准来实现。核对了机组参数和负荷数据之后模型规模大概是二元变量144个连续变量超过4000个加上SDP约束之后整体是一个混合整数半定规划问题。这个规模在普通笔记本电脑上用MOSEK求解大约需要10到30分钟。5.2 三种对比方案的差异我把方案分为三组方案A确定性机组组合忽略所有风电预测误差。方案B区间鲁棒机组组合要求风电在每个时段的任意偏差下都必须可行。方案C本文的分布鲁棒机组组合使用均值、协方差模糊集加线性决策规则。求解完成后把三种方案的日前计划放到1000组蒙特卡洛场景下回测统计实际运行成本和切负荷概率。蒙特卡洛回测的结果大致如下方案日前计划成本万元平均再调度成本万元切负荷概率确定性118.46.28.7%区间鲁棒129.63.10%分布鲁棒124.83.80.4%5.3 从结果反推模型合理性确定性方案虽然日前成本最低但在实际大风偏差下再调度成本最高而且切负荷概率接近9%这基本是不可接受的。区间鲁棒方案完全没有切负荷风险但代价是计划成本最高而且所有机组都为了最极端场景预留了非常大的调整空间很多时段的开机组合明显比实际需要的更“厚”。分布鲁棒方案处于折中位置。它把切负荷概率压到了0.4%以下同时计划成本比区间鲁棒低了不少。这个0.4%并不是因为模型做了妥协而是蒙特卡洛样本里存在个别极端场景超出了模糊集支持域的覆盖范围。如果想进一步压低可以把支持域尾部再扩宽一点或者把协方差矩阵乘以一个大于1的系数模型会自动变得更保守。这类结果其实是DRO最典型的特征通过调节模糊集的大小可以连续控制方案的保守程度。不像随机规划需要重新生成分布假设也不像区间鲁棒只能全有或全无。6. 复现这类DRO-UC代码时绕不开的细节坑6.1 协方差矩阵病态导致求解器报错第一次跑通模型之前我遇到的最大问题是协方差矩阵 (\Sigma) 半正定且接近奇异。原因是风电预测误差的历史样本之间相关性很强尤其相同时刻的不同风电场误差几乎高度相关。YALMIP在构建SDP约束时经常出现数值警告MOSEK甚至直接返回“problem is infeasible”。解决办法是加一个小量正则项把 (\Sigma) 替换为[ \Sigma_{reg} \Sigma \epsilon I ]其中 (\epsilon) 取 (10^{-6}) 到 (10^{-4}) 之间。正则项对最终结果影响很小但能极大改善求解稳定性。6.2 支持域过大导致线性规则过于保守线性决策规则的可行性约束要求对所有 (\xi \in \Xi) 成立。如果把支持域取成上下界整个盒式区域那么所有不等式都会被推到盒子的极端顶点处这会让再调度系数变得很大、很保守。我建议把支持域分成两部分理解物理边界用盒式集合统计约束用模糊集。物理边界可以设置成“实际可能出现的极端情况”但没必要让优化模型在每个极端点都强迫自己满足非常严格的容量约束。比如可以在某些极端方向上允许少量失负荷只要失负荷惩罚够大模型会自动权衡而不是粗暴地让所有约束都满足。6.3 相邻时段爬坡约束的耦合会让端点枚举爆炸如果不做线性决策规则爬坡约束的鲁棒转换很容易因为相邻时段不确定性而变得非常复杂。(\Delta p_{i,t}(\xi)) 可以依赖所有时段的 (\xi)这意味着爬坡约束会同时涉及 (t) 和 (t1) 两个时段的误差变量端点枚举数量随时段数指数增长。我的处理方式是限制线性决策规则的稀疏结构让 (\Delta p_{i,t}) 只依赖当前时段和前后一个时段的预测误差也就是给 (\beta_{i,t}) 加上结构化稀疏约束。这样既能保留线性规则的核心能力又能减少端点枚举和半定约束的规模。牺牲一点点最优性换来的是模型能算得动。6.4 求解器选型MILP还是MISDP需要提前想清楚的是这个模型本质上不是纯MILP而是带有SDP约束的混合整数问题。YALMIP里可以设置多种求解器组合比如用分支定界处理二元变量内点法处理SDP。如果直接用Gurobi去求解带有矩阵锥约束的模型Gurobi是不支持的。我个人的经验是先用MOSEK测试纯连续变量的SDP子问题确认对偶转换没问题之后再加入二元变量。如果机器性能一般也可以把SDP约束进一步松弛为二阶锥约束虽然会损失部分精度但求解速度会明显加快。6.5 蒙特卡洛验证是必须做的一步模型跑完不等于工作结束。无论你用的是矩模糊集还是别的模糊集最终都要用蒙特卡洛场景去回测。我会从历史数据中随机抽取大量场景模拟真实运行中可能出现的各种风电误差组合然后检查每个场景下日前计划是否可行、严重程度如何。这个方法能帮你发现很多建模时看不到的问题。有一次我发现模型在某个罕见的极端组合下出现了较大的弃风量排查了很久才发现是支持域边界设窄了把那个场景挡在了模糊集外面。调整支持域边界并重新训练协方差正则系数之后问题就消失了。我觉得做这类DRO-UC项目最值得花时间的不是把公式写得漂漂亮亮而是把不确定性建模的每一个细节都落地到代码里。线性决策规则的价值在于它让“最坏期望”变得可计算但只有当你对支持域、模糊集、对偶转换和求解器特性都足够熟悉时这套方法才能真正帮你做出既经济又安全的调度方案。如果你也在复现类似模型可以从一个小型6机系统开始先把线性规则和矩模糊集跑通再去扩展到大系统这样排查问题会轻松很多。
