风电出力一上来机组组合这事就再也不是查表填数那么简单了。以前调度员排计划煤机、气机、水电机组的启停和出力基本是“确定题”算一次就完了。现在风电场动不动几十上百台风机出力受风速影响忽高忽低预报曲线和实际曲线之间总有偏差而这偏差直接决定了你排出来的开停机方案靠不靠谱——排太满风电一超发就废风排太保守风电一掉下来就缺电。这篇要聊的就是我在Matlab里实际跑通的一条技术路线用“线性准则”把考虑风电不确定性的分布鲁棒优化机组组合问题变成可求解、可复现的优化模型并给出整套代码实现思路。这个方向适合三类人一是电力系统优化方向的研究生论文需要方法对比和算例支撑二是调度运行相关工程师想了解分布鲁棒优化怎么从论文落到工程三是刚接触优化建模、想用YALMIP加求解器解决实际问题的程序员。我会把数学推导、建模细节、代码结构、踩过的坑全部拆开讲不是教科书式的概念堆砌而是按我实际调通这套流程的先后顺序来。1. 机组组合遇到风电不确定性到底难在哪1.1 一个调度问题的日常困境机组组合Unit CommitmentUC本质上是这样一件事在满足负荷需求、备用要求、机组出力上下限、最小启停时间、爬坡率等一系列约束的前提下决定未来24小时哪些机组开机、哪些停机、每台机组发多少电目标是让总成本最小。确定性版本很好理解把风电当成固定的负负荷代入约束求解即可。但真实运行不是这样风电预测不是铁板钉钉预测误差会随时间、季节、天气变化。如果你把预测值当成真值用当天实际风电低于预测时就得让别的机组紧急爬坡甚至启停调整成本飙升和调度风险都来了。我举一个实际算例中的场景。一个6机系统某时段负荷是450MW风电预测是120MW那么净负荷是330MW。如果实际风电只有80MW净负荷变370MW多出来的40MW缺口就需要其他机组补上。但上一时段你可能刚好把一台100MW的煤机停掉了此时再启动它需要额外费用和时间还可能违反最小停机时间约束——整个优化方案直接失效。这就是不确定性给UC带来的核心麻烦开停机决策必须在不确定尚未揭晓前定下来而出力调整可以在不确定发生后做。1.2 传统处理方式的三个天花板处理风电不确定性的传统方法大概分三路。第一路是蒙特卡洛模拟配合随机规划。你生成大量风电场景每个场景对应一个概率然后求期望成本最小。这个思路直观而且精度高但前提是你能拿到准确的概率分布。实际中风电预测误差分布随时间、季节、气象条件漂移你根本说不清“真实分布”长什么样。更致命的是两阶段随机规划求解规模随场景数爆炸式增长场景一多MILP直接卡死。第二路是鲁棒优化。它不依赖分布只要求保证所有不确定参数落在给定集合比如风电出力区间内时方案都可行。好处是绝对安全坏处是过于保守。你把风电出力设成最坏情况几乎为零来做机组组合意味着几乎全天都要把所有火电拉满经济性惨不忍睹。实际中“最坏情况”发生的概率极低为它牺牲太多正常工况的成本不划算。第三路是机会约束规划。它允许约束在某个概率下失效比纯鲁棒灵活但求解时需要知道分布信息而且机会约束本身不好处理一般要转成近似的确定性等价形式。这个近似的质量决定了结果是否可靠而近似过程中你又回到了“分布到底长什么样”的原始问题。这三路方法各有短板本质上是因为它们只用了“分布已知”或“分布完全未知”两个极端信息。而现实中你手里通常有点数据历史风电误差样本、预测系统的偏差统计。这些数据能告诉你均值、方差、甚至更高阶矩的大致范围但不足以确定准确分布。于是就有了分布鲁棒优化的登场时机。1.3 分布鲁棒组合拳模糊集 最坏期望分布鲁棒优化Distributionally Robust OptimizationDRO的思路介于随机规划和鲁棒优化之间你不需要知道真实分布但你知道真实分布应该属于一个“模糊集”也就是一组满足特定统计特征的分布的集合。优化目标变成在所有可能的分布中找那个让期望成本最大的“最坏分布”再优化你的决策来对抗这个最坏分布。这样做的好处很清楚只用分布的部分信息矩、支撑集范围、距离约束比纯随机规划更“抗分布偏差”又比纯鲁棒优化多了统计信息约束不会无限保守。模糊集的选择直接决定了方法的保守程度和求解难度。我实际使用的是基于矩的模糊集形式如下D { P ∈ P(R^N) : E_P[ξ] μ, E_P[(ξ - μ)(ξ - μ)^T] ⪯ Σ, P(ξ ∈ Ξ) 1 }翻译成人话就是真实分布P的均值等于历史样本均值μ协方差不超过历史样本协方差Σ这里用半定矩阵不等式而且所有可能的ξ都落在支撑集Ξ内。μ描述“平均偏差点”Σ描述“波动幅度”Ξ描述“极端边界”。这三样东西历史数据都能给出相对靠谱的估计。有了模糊集机组组合的目标函数就从“固定场景下的期望成本”变成了“最坏分布下的期望成本”。后者通常写成一个min-max问题直接解非常困难因为内层max是在无穷维概率分布空间上做的。这正是线性准则要解决的问题。2. 线性准则把无穷维决策拉回有限维2.1 两阶段机组组合的数学模样先把问题形式化。我把机组组合建模成两阶段优化第一阶段在风电实际出力不确定之前决定机组启停状态u和基准出力p。这一阶段决策必须对任何可能的风电场景都可行。第二阶段在风电不确定参数ξ揭晓之后对机组出力做调整r(ξ)。这个调整量依赖于实际的ξ是一个“函数”不是固定数值。整体模型写下来长这样min Σ_t Σ_i [ c_i^u u_it c_i^p p_it ] E_P[ Q(u, p, ξ) ]其中Q(u, p, ξ)是给定第一阶段决策后第二阶段的最小再调度成本通常包含机组爬坡调整费用、负荷平衡惩罚等。约束包括系统功率平衡、机组出力上下限、爬坡约束、最小启停时间等。关键点来了第二阶段决策r(ξ)是定义在连续空间上的函数。如果允许它是任意函数这个优化问题在数学上属于无限维问题即使不考虑分布鲁棒也已经足够棘手。加上模糊集后内层要对分布求期望再最大化问题进一步复杂化。常规求解器只能处理有限维决策所以必须对r(ξ)的函数形式做限制。2.2 仿射决策规则推导线性准则也叫线性决策规则、仿射决策规则英文是Linear Decision Rule / Affine Decision Rule就是做这样一个限制第二阶段调整量是风电不确定参数的仿射线性加常数函数。r(ξ) r^0 R ξ这里r^0是常数向量R是系数矩阵。换句话说不管风电实际偏差是多少你只按照“初始调整量 线性修正系数 × 偏差值”的方式响应。比如某机组对风电偏差的响应系数是0.3风电偏差10MW它就多调3MW偏差-20MW它就少调6MW。就这么直接。为什么这个限制合理从工程角度看机组对净负荷变化的响应本来就近似线性AGC系统的调节特性、机组爬坡的线性化处理都让“按偏差比例调整出力”成为实际运行中的常用策略。从数学角度看把r(ξ)代入约束和目标函数后期望部分只需要用到ξ的一阶矩μ和二阶矩Σ——这两个值从历史数据很容易算。而最坏期望部分在矩模糊集下也有闭合表达式。整个问题就从无限维退化成有限维凸优化问题。代入后目标函数的期望部分可以这样处理。第二阶段成本中线性项q^T r(ξ) q^T r^0 (R^T q)^T ξ它的期望是q^T r^0 (R^T q)^T μ。但别忘了我们是在分布鲁棒框架下内层要对P做最大化。在模糊集D下最坏期望有显式形式sup_{P∈D} E_P[q^T r(ξ)] q^T r^0 (R^T q)^T μ γ ||Σ^{1/2} R^T q||_2其中γ是一个由模糊集参数和支撑集共同决定的常数。这个公式是整套方法的核心灵魂第一项是常数部分第二项是均值效应第三项是协方差放大的“风险溢价”。它告诉我们最坏期望成本不仅取决于平均偏差还取决于调整策略在协方差方向上的“暴露程度”。如果你让系数矩阵R很大响应很剧烈最坏情况下的期望成本也会被抬高——这是模糊集里协方差约束的自然结果。约束部分的处理同样关键。第二阶段约束要求A (r^0 R ξ) B u C p ≤ d D ξ对一切ξ∈Ξ成立。这是一个含不确定参数的线性不等式需要转成对所有ξ都成立。如果支撑集Ξ是箱式每维有上下界或椭球式这种“鲁棒约束”可以用对偶理论转化为有限个线性矩阵不等式LMI直接交给求解器处理。我在代码中就是用YALMIP的鲁棒建模功能或者手工对偶转换来做这一步的。2.3 线性准则的适用边界与精度损失讨论线性准则显然是次优的——真实的最优第二阶段策略不一定是线性的可能是分段线性、甚至非线性。但它的优势无可替代可求解、可解释、可工程部署。需要说明的是在线性准则下得到的最优解给出的是原问题的最优值上界因为可行域缩小了。如果你论文里需要下界可以用拉格朗日对偶或者场景松弛去构造。我在对比实验中发现对6机24时段的系统线性准则给出的成本通常比完全场景树下的两阶段随机规划高2%到5%左右但求解时间能快一到两个数量级而且不需要事先指定准确分布。这个精度损失换来了鲁棒性和计算效率对工程实践是划算的。另外需要注意线性准则的效果和不确定参数的维度强相关。风电节点少、误差模式简单时线性准则精度很高当系统有大量风电场、且它们之间的相关性复杂时线性准则可能会过度简化响应模式。此时可以考虑分段仿射决策规则或二次决策规则的扩展方向但建模复杂度大幅上升我暂时没有在Matlab里跑通那部分后续有时间再单独写一篇。3. Matlab实现从数学公式到可运行代码3.1 代码总框架与数据流我先交代一下整个程序的结构。完整实现分四个模块数据模块读取负荷数据、机组参数、风电预测与误差样本生成模糊集参数μ、Σ和支撑集Ξ。建模模块用YALMIP定义第一阶段决策变量启停、基准出力和第二阶段的线性决策系数。求解模块设置目标函数与约束调用外部求解器求解。结果分析模块对比确定性、传统鲁棒、分布鲁棒三种方案的成本、机组组合结果和运行风险指标。数据流大致是原始数据 → 数据处理脚本 → 生成模糊集 → YALMIP建模 → 求解器 → 结果结构化输出。我在实践中把模糊集生成单独写成一个函数文件build_ambiguity_set.m输入是风电误差样本矩阵输出是μ、Σ、支撑集上下界。这样将来换数据或者换模糊集类型只改这个文件就行。3.2 风电误差不确定集构建风电误差数据是整个方法的数据基础。我从历史预测数据和实际出力数据中计算偏差生成每个时段的风电预测误差样本。% 输入: wind_forecast (N_T, N_W, N_S) 预测值; wind_actual (N_T, N_W, N_S) 实际值 % 输出: mu (N_T*N_W, 1), Sigma (N_T*N_W, N_T*N_W) T size(wind_forecast, 1); W size(wind_forecast, 2); S size(wind_forecast, 3); % 误差场景矩阵: 每个样本是一个 T*W 维向量 xi zeros(T*W, S); for s 1:S err wind_actual(:, :, s) - wind_forecast(:, :, s); xi(:, s) err(:); end % 一阶矩 mu mean(xi, 2); % 二阶矩/协方差注意对奇异矩阵做正则化修正 Sigma cov(xi); Sigma 0.5 * (Sigma Sigma); % 强制对称 % 关键: 加正则项防止SDP求解器因数值问题崩溃 reg 1e-4 * eye(T*W); Sigma Sigma reg;这里有个实际经验原始协方差矩阵经常接近奇异原因是风电误差场景之间高度相关——比如某时段所有风电场同时受同一阵风影响。不加正则项的话YALMIP的SDP求解器经常会报“矩阵非正定”的错误。加一个很小的单位矩阵倍数1e-4到1e-6量级不会改变解的本质但能让数值稳定性大幅提升。3.3 YALMIP核心建模与求解器取舍Matlab里做分布鲁棒优化我优先推荐YALMIP其次CVX。YALMIP的优势在于对半定规划SDP和混合整数规划MILP的支持都很成熟能在一个框架内同时处理0-1变量、连续变量、线性矩阵不等式。核心建模代码骨架如下yalmip(clear); % 机组数、时段数 N 6; T 24; % 第一阶段变量 u binvar(N, T, full); % 启停状态 p sdpvar(N, T, full); % 基准出力 % 第二阶段线性决策系数: 每台机组、每个时段对每个风电节点的响应系数 % 这里为清晰起见简化为块结构实际可按需求展开 r0 sdpvar(N, T, full); % 仿射常数项 R sdpvar(N*T, T*W, full); % 线性响应系数矩阵 % 模糊集参数 (由build_ambiguity_set计算得到) % mu, Sigma, xi_lb, xi_ub ... % 目标函数: 确定性成本 最坏期望调整成本 cost_det sum(sum(c_start .* u c_prod .* p)); % 最坏期望部分见公式: q*r0 (R*q)*mu gamma * norm(Sigma^(1/2)*R*q) cost_adj q_lin * r0(:) (R*q_lin) * mu gamma * norm(sqrtm(Sigma) * R * q_lin); objective cost_det cost_adj; % 约束: 功率平衡、出力上下限、爬坡、最小启停时间等 Constraints []; % ... 添加常规UC约束 ... % 求解 ops sdpsettings(solver, mosek, verbose, 2, debug, 1); optimize(Constraints, objective, ops);这一段代码体现了整个建模的关键YALMIP可以混合处理binvar和sdpvar、norm和sqrtm这种非线性表达式。目标函数里那个norm(Sigma^(1/2) * R * q_lin)YALMIP会自动把它转成一个二阶锥约束的形式交给MOSEK处理。求解器方面我的建议是SDP部分用MOSEKMILP部分用Gurobi。MOSEK对半定规划的支持很稳Gurobi的混合整数求解速度业界顶尖。你可以在solve前用sdpsettings分别指定或者干脆在sdpsettings里写solver,mosek遇到整数变量时YALMIP会自动调用MOSEK的MISDP求解能力。不过说实话纯MOSEK处理大规模MISDP性能一般我的做法是对规模较小的测试系统比如6机24时段直接用MOSEK一把梭没问题。对更大规模系统用Benders分解主问题是启停决策的MILPGurobi子问题是给定启停下的SDPMOSEK割平面迭代求解。这个方案我后续有机会单独写。4. 实操记录一个6机系统的完整跑通流程4.1 算例设置与三种对比方案我用的是一个经典的6机24时段测试系统机组参数参照IEEE RTS-24系统的典型数据做了简化。负荷曲线用典型日负荷峰谷差约40%。风电数据用的是某个风电场的历史出力数据预测值用持续法persistence method生成误差样本直接取历史误差。三种方案对比确定性方案风电取预测值直接求解UC。传统鲁棒方案风电取区间最坏值预测值减一个偏差带用鲁棒优化求解。分布鲁棒方案本文方法模糊集取矩约束 箱式支撑集。目标函数统一用总运行成本启停成本燃料成本备用约束统一设5%负荷。4.2 关键代码走读功率平衡约束与线性决策的展开功率平衡约束不能简单写成一个等式。确定性部分要求基准出力加上调整量的期望等于净负荷Σ_i p_it μ_t D_t - W_f_t这里W_f_t是风电预测值μ_t是该时段风电误差的均值。由于误差均值通常很小这个等式基本和确定性模型差别不大。真正关键的是第二阶段约束。对每个可能的风电偏差ξ需要满足Σ_i [r^0_it Σ_w R_it,w ξ_w] (W_f_t ξ_t) ≥ 0 不弃风即功率平衡下界 Σ_i [r^0_it Σ_w R_it,w ξ_w] (W_f_t ξ_t) ≤ L_imb 超过可接受弃风这类约束必须对支撑集Ξ内所有ξ成立。我用YALMIP的robust方法实现% 支撑集: 箱式区间 xi_unc uncertain(xi_sdp); Constraints [Constraints, ... sum(r0(:, t)) sum(R(:, t)) * xi_unc (W_f(t) xi_unc(t)) 0 : balance_lb]; Constraints [Constraints, ... sum(r0(:, t)) sum(R(:, t)) * xi_unc (W_f(t) xi_unc(t)) L_imb : balance_ub];YALMIP会自动对这种线性不确定约束做鲁棒对偶转换生成等价的确定性约束。需要注意的是uncertain变量必须和sdpvar区分开而且YALMIP只支持线性或二次约束下的自动对偶。如果支撑集太复杂建议手动对偶可控性更强。目标函数中的最坏期望部分% 先确定第二阶段单位调整成本向量 c_adj q_lin repmat(c_adj(:), 1, 1); % 维度 N*T × 1 % 最坏期望调整成本 常数 均值效应 协方差风险项 term_mean q_lin * r0(:) (R * q_lin) * mu; term_risk gamma * norm(sqrtm(Sigma) * R * q_lin); objective cost_det term_mean term_risk;这里gamma是根据模糊集和支撑集算出来的系数。如果模糊集只约束均值和协方差且协方差约束取“⪯”形式最坏期望风险项前系数γ通常为1即直接是协方差范数的对偶形式当也有支撑集约束时γ会随支撑集大小变化。我实际测试时γ的范围大概在0.8到1.2之间具体数值通过求解一个小的对偶问题得到。代码里我先固定γ做初步求解再通过数值检查调整效果稳定。4.3 结果分析的四个看点跑完结果后我建议从四个维度分析第一总成本对比。分布鲁棒方案的成本通常比确定性方案高3%到6%比纯鲁棒方案低8%到15%。这个区间很有规律如果你的模糊集构造合理分布鲁棒的成本会落在两者之间且更靠近确定性方案。第二开机组合的稳健性。确定性方案在风电偏差较大时会频繁启停机组来补差额而分布鲁棒方案的开机组合在大多数场景下都能保持稳定。我统计过确定性方案在500个随机测试场景中有41个场景发生失负荷或弃风分布鲁棒方案只有3个。这就是鲁棒性的直接表现。第三LDR响应系数R的分布。看R矩阵的热力图可以发现偏差响应主要集中在爬坡速度快、调节成本低的机组上煤机等慢速机组的响应系数很小。这说明模型自动做了“最优响应分工”是线性准则带来的可解释性红利。第四模糊集参数灵敏度。把Σ放大到原来的1.5倍总成本上升约4%把均值μ偏移一个标准差成本上升约2%。这个灵敏度分析能告诉调度员最值得提升的是预测精度减小μ还是出力波动控制减小Σ对实际运行有指导价值。5. 常见的坑与排查技巧实录5.1 求解器报错与不可行问题我遇到最多的错误是“Infeasible problem”。大部分情况下不是模型真的无解而是数值问题协方差矩阵非正定加正则项前面已经说过。支撑集设置过窄如果ξ的箱式范围比数据中实际出现的误差小鲁棒约束可能无解。解决办法是先看误差数据的分位数用90%或95%分位数作为箱式边界别拍脑袋设。目标函数中sqrtm(Sigma)数值病态改用chol(Sigma)但注意YALMIP里对chol的支持一般不如直接用sqrtm加正则。YALMIP的debug开关很有用。在sdpSettings(debug, 1)下YALMIP会给出更详细的不可行原因提示。配合check(Constraints)逐条检查约束残差能快速定位是哪条约束出了问题。5.2 SDPMILP的性能瓶颈6机24时段规模不大MOSEK一般几分钟内能出结果。但如果你想扩展到时序更长或者机组更多性能会急剧恶化。原因是SDP的尺度随不确定参数维度三次方增长——风电节点数从3加到10SDP矩阵规模翻好几倍。我实践中的应对方案降维如果多个风电场相关性很强用PCA把误差场景降维到3到5个主成分再用主成分基做模糊集。精度损失很小求解速度提升明显。支撑集简化把多维箱式支撑集改成一个椭球约束可以显著减少鲁棒约束的对偶变量数量。热启动用确定性方案的解作为MILP部分的初始解能减少约30%的求解时间。方法是先算一次确定性UC把启停变量赋值给u的初始值。5.3 保守性与经济性的天平分布鲁棒优化不是万能的它的保守程度完全由模糊集决定。模糊集越大结果越保守越小越接近随机规划。这里的关键是让模糊集“贴着”真实数据的可信范围走。我踩过的坑是一开始把模糊集构造得过大协方差乘了3倍结果成本比纯鲁棒还高方法优势全没了。后来改用交叉验证来确定模糊集参数——把历史数据分成训练集和验证集在训练集上构造模糊集在验证集上看方案的失负荷率。最终选定能保证验证集失负荷率小于2%的最小模糊集。这个方法比拍脑袋调参靠谱得多。另外风电误差的均值μ如果估计不准对结果影响很大。我建议用带漂移项的动态均值模型替代简单历史平均尤其在风电出力有明显季节性波动的场景下。动态均值的实现不复杂按月份或天气类型分组分别计算每个组内的均值μ_k然后模糊集取所有组的均值包络。这算是一个低成本高收益的改进方向。6. 个人体会与扩展方向跑完这套分布鲁棒机组组合我最直观的感受是鲁棒性和经济性的平衡不是靠拍脑袋调参数而是靠把不确定性的“信息量”准确传进模型里。线性准则的智慧在于它承认了我们无法预知风电的每一次波动但利用了“偏差和调整量之间大致呈线性关系”这个工程常识把复杂的无限维问题压缩到可以求解的尺度。相比之下纯鲁棒像是一刀切的“最坏情况预案”随机规划像是依赖天价场景的“精算师”而分布鲁棒优化更像是“拿着统计报告做预案的老调度”——知道大概范围、知道平均趋势、知道波动幅度然后据此排计划。代码层面的几条经验再总结一下模糊集的正则化处理是数值稳定的关键YALMIP的uncertain变量极大简化了鲁棒约束的建模求解器的选择上MOSEKGurobi组合最稳结果分析一定要做灵敏度测试否则你的模糊集参数在审稿人或者领导面前经不起追问。后续我打算在这个框架上做三个扩展一是加入多风电场联合模糊集用Copula结构建模相关性的高阶信息二是把时间耦合约束储能的充放电决策纳入线性准则框架做多时段联合调度三是将这套方法迁移到配电网的有功-无功联合优化里。等代码跑通了再来更新。