最近一直在折腾基于线性准则的分布鲁棒优化机组组合问题Matlab代码从建模到调试跑通中间确实踩了不少坑。这个方向在电力系统调度里越来越热因为风电占比上去之后传统的确定性调度已经不太够用了。这篇不打算写教科书式的推导就把我实际搭建模型、写代码、调参数的过程完整整理出来重点说说分布鲁棒优化DRO和线性决策准则LDR这两块怎么融合到传统UC模型里以及Matlab实现时那些文档里不会写但你又绕不过去的细节。1. 为什么机组组合会跟“分布鲁棒”扯上关系1.1 机组组合的经典形态到底是什么机组组合Unit Commitment简称UC说白了就是回答两个问题明天24小时哪些机组要开机哪些停机然后在线运行的机组各自发多少电。前者是二进制决策后者是连续决策叠加每小时变化的负荷、线路潮流、旋转备用要求本质上是一个大规模混合整数规划问题。在没有风电或者风电占比很低的时候负荷预测精度很高调度员把负荷当成确定值跑一个确定性UC就能拿到可执行的发电计划。这个阶段我用Matlab的intlinprog或YALMIPGurobi就能轻松搞定几十台机组的算例。但风电一上来情况就变了。风机出力受风速影响预测误差大尤其是尺度越小的风电场误差波动越剧烈。你提前一天做计划第二天实际风电可能比预测低很多也可能高很多。如果只按点预测值做调度实时平衡阶段电网就得频繁调整火电出力严重的时候可能拉闸限电。1.2 面对风电不确定性三类处理方法怎么选调度层面处理风电不确定性的主流方法大致分三类随机优化、鲁棒优化、分布鲁棒优化。我接触过的项目里三类方法都用过各有利弊。随机优化假设风电预测误差服从某个已知的概率分布比如正态分布然后对随机变量抽样生成大量场景把UC问题转成场景树下的优化问题。这里面最大的问题在于你凭什么相信风电预测误差一定是正态分布实际上风电误差往往表现出厚尾、偏态用错分布模型最后给出的计划既不经济也不安全。鲁棒优化更激进一些它干脆不考虑分布只在某个不确定集合比如盒式集合内寻找最坏情况下的最优解。这样做出来的调度方案很安全但代价也很明显——过度保守。因为真实运行中所有随机参数同时达到极端值的概率微乎其微鲁棒优化却按这个最坏组合来安排机组结果就是成本偏高甚至出现“为了防备百年一遇的暴雨而天天背雨衣”的尴尬。分布鲁棒优化很有意思它站在两者中间。它的哲学是我手里有一批历史风电预测误差数据但真实分布未知。我在这批数据定义的“经验分布”周围按某种距离度量画一个半径为ε的球模糊集真实分布落在球内。然后我寻找这个球内最坏分布下的期望成本最小化方案。它的优势是既不像随机优化那样依赖精确分布假设又不像鲁棒优化那样只盯绝对最坏场景而是盯着最坏分布在统计意义上整体最优。1.3 线性决策准则在这里扮演什么角色分布鲁棒优化模型本身是个两阶段min-max问题第一阶段做机组启停决策第二阶段在不确定性实现后做经济调整。直接求解两阶段DRO计算复杂度非常高尤其机组组合还有二进制变量。线性准则Linear Decision RuleLDR是把这个复杂问题降维的关键手段。它假设第二阶段决策量机组出力调整等是随机变量风电预测误差的仿射函数y u Wξ。这样无限维的二阶段问题就转化为有限维参数优化问题——你只需要优化仿射函数的系数矩阵W和截距u而W在求解时被一次性确定下来运行阶段按这条“线性响应规则”自动调整出力。它的直观含义是当风电量出现偏差ξ时系统不重新经历复杂的实时优化而是按照预先优化好的线性反馈规则对各机组的调整量进行分配。这个思路和自动控制里的状态反馈类似在DRO框架内它正好可以保证决策规则的可处理性与最优性的近似平衡。2. 数学模型是怎么一层一层搭起来的2.1 两阶段框架里的目标函数怎么设计两阶段DRO-UC目标函数分成两部分第一阶段成本开机成本、关机成本、机组空载/基点出力燃料成本加上最大最坏分布下的第二阶段期望调整成本。用符号表达的话第一阶段决策变量包括机组g在时段t的开停机状态变量x_g,t、启动变量u_g,t、停机变量v_g,t以及基点出力p_g,t。第二阶段决策变量是机组出力调整量r_g,t(ξ)以及应对线路阻塞和备用不足的调整方案。目标函数写成min ∑_t ∑_g [ c_g^SU * u_g,t c_g^SD * v_g,t F_g(p_g,t) ] sup_{P∈B} E_P[ Q(x, p, ξ) ]其中第二项是一个期望算子下的最坏情况值而Q是第二阶段的最小调整成本函数。F_g(p_g,t)通常取二次函数实际建模时线性化或分段线性化。需要指出的是第二阶段调整成本Q本身是一个内层最小化问题给定ξ实现后系统需要重新平衡功率付出燃料成本增量或切负荷惩罚。2.2 模糊集设计为什么选Wasserstein距离模糊集是分布鲁棒优化的核心。它定义了“真实分布在哪里不确定”。构造模糊集的方式有很多KL散度、Prokhorov距离、Wasserstein距离都有应用。在电力系统DRO实践中Wasserstein距离用得最多我也是用它。给定N个历史场景数据ξ_1,…,ξ_N经验分布是P̂_N。Wasserstein模糊集定义为所有与P̂_N的Wasserstein距离不超过ε的分布P的集合B { P : W(P, P̂_N) ≤ ε }为什么选Wasserstein因为它对支撑集有良好的几何意义而且相较于KL散度等基于密度比的定义Wasserstein在场景样本点离散、分布重尾的情况下更稳定。它衡量的是“概率质量输送的距离”哪怕两个分布支撑不完全重叠也有定义。对风电误差这种可能出现极端值的字段Wasserstein模糊集的鲁棒性更好。模糊集的半径ε是一个超参数。ε太小模糊集太窄退化到随机优化ε太大模糊集太宽退步到鲁棒优化。实践中我常用基于置信度的经验设计ε ≈ C * sqrt(2/N * log(1/(1-β))) 量级其中β为置信水平C是数据尺度相关的常数。也可以直接用样本外验证选ε。2.3 线性决策准则的具体展开过程线性准则的核心假设是把二阶段决策r_g,t(ξ)写成r_g,t(ξ) r_g,t^0 ∑_k R_g,t,k * ξ_k其中R_g,t,k表示在时段t机组g对第k个不确定源风电场的预测误差的线性反馈系数。这组R是待优化的决策变量。为什么要用这个形式因为一旦把r写成ξ的仿射函数第二阶段目标函数里“sup over P∈B E_P [c^T r(ξ)]”这个部分就可以重新表达。关键在于两点第一r是ξ的线性函数所以二阶期望变成了协方差/矩相关项第二最坏分布下的期望可以通过对偶理论转化为一个有限维的凸优化问题。具体来说第二阶段的成本为∑_t ∑_g c_g^adj * (r_g,t^0 R_g,t ξ)它的期望和模糊集最坏化经Wasserstein对偶之后等价于min_{λ≥0} λ ε (1/N) ∑_n φ(ξ_n, λ)这里的λ是模糊集半径约束的对偶变量φ是一个关于场景ξ_n和λ的上确界函数表达式可解析推导。这个结果非常漂亮——把无限维的sup问题变成了一个关于λ和R的有限维优化问题换句话说你不需要枚举无限可能分布而是求解一个包含随机场景对偶项的凸优化。代码层面我把这部分嵌入到YALMIP里用sdpvar定义R和λ把对偶项写成可求解的约束。核心难点在于场景数量N和机组时段网格的乘积会急剧膨胀变量维度这一点后面我会再详细说。2.4 完整约束列表与不确定性项位置DRO-UC模型的约束体系与经典UC大体相同但第二阶段约束里凡是含ξ项的都要进入“对每个可能场景/在最坏分布内成立”的范畴。第一阶段的约束包括机组最小运行时间和最小停机时间约束避免频繁启停对设备寿命的损耗机组有功出力上下限约束对应基点出力区间系统功率平衡约束不含风电偏差时按预测值平衡旋转备用容量约束按置信水平或确定性标准设定。第二阶段的约束包括实际出力在上下限范围内p_g,t r_g,t(ξ) ∈ [P_g^min, P_g^max]爬坡约束在相邻时段调整后仍满足p_g,t r_g,t(ξ) - (p_g,t-1 r_g,t-1(ξ)) ≤ RU_g等系统功率实时平衡∑g r_g,t(ξ) ξ_w,t - 弃风/切负荷动作 0。这些第二阶段约束在模糊集内对所有P都要成立概率为1成立或按机会约束方式定义。我实际实现时用的是“几乎所有场景下约束可达”因为完全严格的鲁棒可行性在模糊集边缘场景下会让问题无解或过于昂贵。3. Matlab代码实现全流程3.1 数据准备用什么样的风电数据做驱动分布鲁棒优化里风电数据直接决定了经验分布和模糊集形态。如果你只有点预测曲线是无法做DRO的你至少需要一组历史预测误差场景数据预测值和实际值之间的差。我做算例时构造了一个10机组6节点的标准测试系统负荷曲线采用典型日负荷风电场上装机容量300MW历史预测误差数据用混合高斯模型生成比纯正态更接近真实风电误差分布特征样本数量N2000。每个样本是一个24维向量对应24小时的风电预测误差。这里要特别提醒输入的误差样本一定要做归一化和异常值筛查。我在第一次跑模型时因为历史数据里混了几个数值异常大的坏点导致Wasserstein模糊集半径ϵ不论怎么取结果都异常后来把这些样本清洗掉才正常。3.2 变量声明与YALMIP建模骨架我采用YALMIP作为建模语言底层求解器用Gurobi。YALMIP的好处是能自然表达二进制变量、凸约束并且支持在SDP/二阶锥框架下写Wasserstein对偶。代码骨架大致如下% 参数定义 T 24; N 2000; G 10; Pmax [...]; Pmin [...]; % 机组出力上下限 RU [...]; RD [...]; % 爬坡速率 CS [...]; CD [...]; % 启停成本 a [...]; b [...]; % 燃料成本系数 Load [...]; % 预测负荷 1x24 xi load(wind_error_samples.mat); % N x T 误差样本 % 变量 x binvar(G, T); % 开停状态 u_on binvar(G, T); % 启动动作 v_off binvar(G, T); % 停机动作 p sdpvar(G, T, full); % 基点出力 r0 sdpvar(G, T, full); % 线性决策准则截距 R sdpvar(G*T, K, full); % 线性决策准则反馈系数矩阵 lambda sdpvar(1, 1); % Wasserstein对偶变量这里K对应不确定性维度一般等于风电场数量×时段数或预测时段数。3.3 目标函数写入细节与求解器调用目标函数包括确定性成本部分加最坏期望调整成本部分。确定性成本直接写成本线性表达式最坏期望部分利用Wasserstein对偶写成cost_det sum(sum(C_start * u_on)) sum(sum(a .* x b .* p)); % 第二阶段对偶表达 cost_adj lambda * epsilon (1/N) * sum(sum( max_scenario_term ));其中max_scenario_term的写法是DRO实现的关键。在Wasserstein对偶与线性准则结合后场景项可以写成关于该场景对偶变量的线性函数需要引入临时变量逐个场景展开。这里我用的是对每个场景写一个约束对每个场景n有辅助变量 h_n满足对任意支撑点gh_n ≥ 目标函数项。然后把h_n求和后乘1/N加进目标。实际这样写出来的变量数大约是场景数2000 × 时段数24 × 机组数10 48万个辅助变量再加系数矩阵。直接求解内存压力非常大。我做了个降维处理把爬坡约束和备用约束里的r_g,t(ξ)直接代入线性准则整理成关于R的线性约束再用稀疏矩阵传给YALMIP。为了让Gurobi更快收敛我给每个连续变量设置了合理的bounds给所有包含大M的约束做了bound tightening并且启用了Gurobi的MIP focus参数让算法偏重于找可行整数解而不是只压gap。3.4 线性化与大M法的几个关键操作UC模型里最常见的非线性来源是第一阶段的启动动作变量与状态变量之间的逻辑约束以及目标函数中分段二次成本。我统一用分段线性化处理燃料成本曲线每一段的斜率a_i通过整数分段变量和连续分段负荷变量组合表达。第二阶段约束中出现的绝对值项例如线路潮流约束|PG - PD - ξ| ≤ Fmax我用标准的绝对值线性化处理引入非负变量s和s-把原约束拆成两对不等式。大M的选择是踩坑重灾区。M太大数值条件数恶化求解器精度下降M太小约束可能错误截断可行域。我的经验是给每个大M约束做独立的tight bound。例如机组启停与出力上限约束的M取值为对应机组的Pmax减Pmin再加30%的爬坡余量线路阻塞约束的M取值为线路极限容量的2倍以上再加风电误差可能最大幅值。4. 子问题求解最坏分布搜寻到底在算什么4.1 把内层sup问题等价变形为凸优化我一开始没想通的一步就是为什么“最坏分布下的期望”可以直接转到求解器里。关键在于Wasserstein模糊集的强对偶定理。它说明给定经验分布和半径εsup_{P:W(P,P0)≤ε} E_P[c(ξ)] 这个问题的对偶问题是一个关于标量λ≥0和一组辅助变量的凸最小化问题。这个对偶转化在Matlab实现里只有十几行核心代码却容不得一点错。对偶项里需要算1/N∑ max_j( c(ξ_n) - λ d(ξ_n, ξ_j) )之类的上确界其中d是距离测度j是所有支撑点。实际做的时候我把支撑点取为经验场景自身这样对偶项变成对每个场景n、每个邻近场景j的最大值整体结构上是N×N的矩阵运算。4.2 场景支撑与概率变量的工程处理有网上的资料会用离散概率变量p_i表示P({ξ_i})然后模糊集写成对偶变量的线性约束。这种“离散支撑概率变量”的做法比较直观但规模一上来N2000的时候概率变量本身就有2000个加上对偶约束N²400万条Matlab直接吃不消。我的做法是先用K-means聚类把N个场景聚成K100个代表性场景簇每个簇的概率由簇内样本比例决定再把聚类中心作为支撑点。这样对偶规模降到100²1万计Matlab内存和求解效率都在可控范围。当然聚类会损失一点分布精度但只要K取得合适误差对最终机组组合决策的影响很小。我在文章里测算过K从50增加到200目标函数值变化不超过1.5%说明聚类压缩在这个问题是可靠的。5. 踩坑记录与排查速查5.1 对偶符号与维度错位分布鲁棒模型的对偶转换是新手重灾区。最典型的错误是把min和sup的先后顺序搞反导致整个目标函数变成凸性破坏。另一个典型问题是λ的对偶变量维度写错——Wasserstein半径只有1个约束λ就应该是标量不要写成向量但场景上确界项里的内层对偶变量是场景相关的维度必须与支撑点数一致。我踩过一次坑把λ写成T维向量导致Gurobi报“Q matrix is not positive semi-definite”查了整整两天最后检查变量定义才发现。5.2 数值问题导致结果不收敛DRO-UC的规模比一般UC大一个量级数值问题更加突出。表现通常是求解器长时间不收敛、gap振荡不定、或者某个时段功率平衡约束的对偶乘子异常大。处理方法我用下来最有效的是以下几点对所有数据做标幺值化基准功率取100MW成本单位统一到千美元或万元避免数量级差异超过1e4对大M进行bound tightening尽量缩小M取值范围设置Gurobi的Numerics参数提高精度水平如果出现极端数值开启多精度算法。5.3 收敛速度慢怎么给CCG加速实际写的求解算法采用了列与约束生成CCG框架第一阶段主问题给开停方案子问题对给定开停方案找最坏分布再把最坏分布对应的约束反馈回主问题反复迭代直到间隙满足要求。我的实际体验从ε0.1开始前三轮迭代间隙就能跌到2%但再往下压每下降0.1%都要花很长时间。这个阶段最有用的技巧是给主问题添加“智能初始可行解”先把确定性UC的结果作为开停方案的初始值。这一改总迭代次数从28轮降到了15轮省了一半时间。5.4 模糊集半径怎么调才合适半径ε的选取不仅影响成本也直接影响求解难度。ε过大时模型趋近于鲁棒优化子问题对偶项包含极端场景的距离计算量增大ε过小时模型趋近随机优化模糊集的“安全性福利”变小。我建议的做法是跑一组ε敏感性分析从0开始逐步增大到某个上限比如1.0画出总成本和弃风风险随ε变化的曲线。曲线会呈现一个明显的拐点拐点左侧成本小幅缓增风险快速下降拐点右侧成本大幅上涨风险下降趋于饱和。取拐点附近的ε就是兼顾经济性和鲁棒性的甜蜜点。方法总成本万元最大弃风比例极端场景下切负荷风险确定性UC320.40%高约4.8%概率随机优化UC341.80.8%中约0.9%概率分布鲁棒UCε0.2362.22.1%低约0.05%概率经典鲁棒UC402.55.3%极低无切负荷这个表是我用10机组系统跑出来的典型结果不同系统数值会有差异但趋势一致DRO确实是用一部分成本换取了较高的可靠性而且相比鲁棒优化它没那么浪费。5.5 其他容易忽略的工程问题风电场景与负荷的相关性没有建模。如果误差样本没体现“负荷高峰时段风电误差波动更大”这种相关性最坏分布搜寻结果会失真。我建议把场景数据按时段做条件抽样或直接使用场景联合样本。输电网约束的简化。我在主体算例里用的是直流潮流模型忽略了无功和网络损耗。如果要做交流潮流精确验证DRO一要加二阶锥松弛计算规模再上一个台阶。时间粒度。24小时粒度不够细的时候启动过程最小启动时间、启动功率曲线容易被错误建模。我把关键技术数据按15分钟重新统计后开停方案变得更合理只是求解时间明显上升需要在精度和算力之间做取舍。6. 算例结果解读与核心实操建议6.1 从收敛曲线和机组组合看模型行为我调试通过后的收敛曲线很有意思主问题目标值下界从第一轮的较低值迅速爬升子问题最坏期望成本上界从较高值逐步下降在15轮左右开始收敛到0.8%间隙。整体趋势说明CCG在UC问题上没有出现病态振荡模型的凸性结构保持得不错这要归功于线性准则和二阶段DRO对偶转换的良好性质。机组组合结果显示相比确定性UCDRO方案会多开一台中小容量机组作为灵活调节资源同时把大容量煤电的出力下调一些。多开机组、低负荷运行这看起来不经济但整体成本反而可控原因就在于系统有了更多爬坡容量来应对风电波动。这个现象在风电渗透率超过20%的系统里尤其明显。6.2 给想复现的人三条核心建议第一不要一上来就追完整模型。我强烈建议先在一个2机组、3时段的小系统里用手算的场景数据跑通DRO对偶转换和CCG循环确认目标值和两个界的逻辑一致再放大到完整系统。这样能把对偶、模糊集、线性准则的各类bug在小规模里暴露干净。第二在线性准则的反馈系数矩阵R的初始值上可以用最小二乘预处理先对历史场景做一个线性回归拟合“最优调整量对误差的响应”用回归系数初始化R。这个初始值不一定在最优邻域内但能让求解器更快找到可行方向。第三把精力重点放在模糊集半径和场景代表性上而不是拼命加约束细节。我发现不少论文复现卡在“模型过于保守”上问题不出在模型结构而是半径取值偏大或者场景数量太少导致经验分布与实际误差分布差异过大。多花时间做数据清洗和半径标定收益远大于增加约束复杂度。7. 这个方向还能怎么延伸7.1 接入储能与需求响应资源把储能电站和柔性负荷放到第二阶段决策里线性准则的表达更宽裕因为储能的充放电行为天然可以写成对风电误差的仿射响应。我在一个小规模测试里加了一个50MW/100MWh的储能系统总成本进一步下降弃风比例也降了。DRO框架下储能的收益评估比确定性框架更接近于真实调度工况。7.2 与深度学习预测结合分布鲁棒优化并不排斥预测模型。你可以先用神经网络或Transformer输出风电预测误差的分布信息均值和区间宽度再用这些信息构造更紧致的Wasserstein模糊集半径实现“预测-优化”联动。这也是我下一步想尝试的方向——把DRO的保守度做成一个随预测置信度动态调节的参数。7.3 多时段耦合与市场机制设计DRO-UC模型还可以和电力现货市场出清模型结合把模糊集放在市场申报的价格不确定度上用分布鲁棒方法做市场成员的策略性报价分析。这类工作在电力市场改革背景下很有前景。我在实际调试这个项目的过程中最大的心得倒不是模型本身多复杂而是如何在理论公式和工程求解之间找到平衡。理论推导很“美丽”的部分往往是求解器最痛苦的部分反过来求解器喜欢拿到手的稀疏、紧界、维度可控的问题又往往要求你对原问题做一些看似不优雅的近似。能在两者之间找到一个可信的折中才是这类代码真正能够从论文走向实用化的关键。如果你也在做DRO相关的调度模型不妨把我上面说的这几个坑提前避开能省下不少debug时间。