做机组组合的人应该都听过一句话风电预测曲线和实际出力永远对不上。传统确定性方法在风电渗透率不高时还能凑合可当预测误差动辄几百兆瓦一个拍脑袋的备用容量可能让调度中心在凌晨三点焦头烂额。我最近用Matlab把一套基于线性准则的分布鲁棒优化机组组合模型完整实现了从数学推导到代码落地踩了不少坑。今天把整个思路、建模细节和程序实现串起来讲一遍包括风电不确定性建模、Wasserstein模糊集构造、两阶段决策规则以及用YALMIP调Gurobi求解的完整流程希望给正在做类似课题的人省点时间。1. 为什么机组组合必须认真对待风电不确定性1.1 确定性调度的局限在哪里机组组合本质上是一个时域耦合的大规模混合整数优化问题决定每台火电机组在各时段是开机还是停机并给出满足负荷与备用需求的有功出力计划。传统做法把风电当作已知量通常取预测值或者“预测值加一个固定比例备用”来打保险。风电渗透率在10%以内时这种近似还能接受。但渗透率超过30%之后问题就暴露了最典型的是两个场景一是预测值偏乐观时系统实际需要的上调备用不足。晚高峰负荷上来、风电突然掉出力备用机组要么还没热启动要么爬坡速度跟不上等着调度员的就只有切负荷。二是预测值偏保守时按高风电出力安排开机数量夜间实际出力却远低于预期部分机组被迫压到最小技术出力以下只能弃风甚至倒逼最小出力更高的机组停机再启动成本和启停次数一起飙升。固定备用系数的方法本质上是“赌误差不超过某个经验的百分数”。它不能回答一个基本问题极端场景到底会坏到什么程度以及为了应对这种程度需要付出多少成本。这就是随机优化、鲁棒优化和分布鲁棒优化登场的背景。1.2 三种不确定性优化路线怎么选先说随机规划。它给风电出力预设一个离散场景集每个场景有发生概率目标变成“全场景期望成本最小”。好处是决策有经济性坏处是你得相信场景概率是对真实分布的准确近似。概率设错了优化解很可能只在训练场景里好看换一批实际数据就翻车。鲁棒优化走另一个极端要求所有可能的风电曲线都在约束可行域内目标是最坏情况成本最小。它不需要概率分布只需要不确定性集合比如盒式区间。缺点是“最坏情况”常常是永远不会发生的极端情形为了保证这种极端系统得多开很多机组备用冗余高到经济性很差。分布鲁棒优化夹在两者中间。它假设风电出力的真实概率分布落在某个模糊集内这个模糊集以历史场景的经验分布为中心、以某个距离半径为界。优化的目标是在模糊集内的所有分布中找“最坏的那个分布”再对这最坏分布求期望成本最小。换句话说它对概率分布本身做鲁棒而不是对场景集合做鲁棒。既保留了经济性又不会因为一个极端点就过度设计近几年在电力系统调度文献里热度很高。这套思路落到机组组合上就非常合适机组启停计划是典型的“今天做决定、明天见真章”必须今天就在不知道明天风电分布真实样子的情况下做决策。分布鲁棒框架天然支持这种“现在做决定、事后看结果”的两阶段结构。1.3 标题里的“线性准则”到底指什么很多初次接触这个概念的人会被“分布鲁棒优化”几个字劝退其实它背后有一个很实用的工程近似思路第二阶段的可调决策被限制为风电出力偏差的仿射线性函数。举个例子第一阶段决定机组i在时段t是否开机并给一个基础出力。第二阶段风电实际出力和预测值出现偏差ξ之后机组会在基础出力附近增减一个调整量。这个调整量如果是直接由优化器自由决定两阶段问题会和概率分布耦合得很深几乎没法直接求解。线性准则把这个调整量写成关于ξ的线性函数比如adj_i,t(ξ) a_i,t b_i,t · ξ其中a是常数项b是灵敏度系数。这样原来那个嵌套的max-min-期望问题通过线性对偶变换可以转成一个确定性的混合整数线性规划。这就是“基于线性准则”的含义。有人可能会觉得线性函数限制太强风电偏差大时调整策略不一定是最优的。但从工程角度看线性决策规则带来的成本损失通常在5%以内换来的却是计算速度快几十倍这是非常划算的交易。我在实际测试中规模稍大的算例如果用完全自由的两阶段模型求解器经常一个小时内出不了可行解而线性准则版本基本几分钟内就能收敛到MIP gap可接受的水平。2. 分布鲁棒机组组合的数学模型与转化2.1 基础机组组合模型建模从经典的确定性机组组合出发。系统有I台火电机组调度周期T个时段每台机组有最小开机/停机时间、出力上下限、爬坡速率、启动与停机成本。目标函数是燃料成本加启停成本燃料成本用二次函数或者分段线性近似min Σ_t Σ_i [ f_i(P_i,t) · u_i,t SU_i,t SD_i,t ]约束至少包括功率平衡Σ_i P_i,t W_t^f D_t风电取预测出力时系统发电必须等于负荷。旋转备用Σ_i (P_i,max - P_i,t) · u_i,t ≥ R_t 风电备用需求。出力上下限P_i,min · u_i,t ≤ P_i,t ≤ P_i,max · u_i,t。爬坡约束P_i,t - P_i,t-1 ≤ RU_i以及下坡约束。最小启停时间约束用经典的MILP三不等式形式做线性化。启停成本约束许多模型里会把启动成本与停机时间挂钩。这套模型本身不复杂难点在风电不确定性的注入。2.2 风电不确定性的模糊集构造常见做法是用历史数据或场景法生成N个风电出力样本ξ_1,...,ξ_N构造经验分布P_N。然后定义以P_N为中心、以Wasserstein距离为度量的模糊集D { P : d_W(P, P_N) ≤ ε }这里的ε也称模糊集半径直观理解是真实分布和经验分布之间的“概率距离”上限。ε越小模糊集越小你越信任历史数据ε越大模糊集越大鲁棒性越强但保守性也越强。Wasserstein距离相比KL散度有个明显好处它允许分布的支撑集发生变化也就是说真实分布可以出现“历史样本里没见过的出力值”这在风电场景里非常重要因为极端低风和高风事件本来就稀疏历史样本集未必能覆盖。对模糊集做参数敏感性测试是必须的步骤。我在测试中发现ε取0.05时解基本和随机规划差不多取0.3时系统会明显多开机组夜间开机台数可能增加一到两台这就是鲁棒性代价的直观体现。选择合适ε的实用办法是拿一段真实历史风电数据做后验评估比较不同ε下解的期望运行成本与最坏运行成本。2.3 两阶段分布鲁棒模型怎么搭把风电不确定性引入后模型写成两阶段形式第一阶段确定机组启停计划u和基础出力P0。第二阶段风电实际场景ξ发生后机组可以调整出力也可以切负荷或弃风目标是在所有可能分布中找最坏分布下期望调整成本最小的方案。形式化一点可以写成min Σ成本 sup_{P∈D} E_P[ Q(P,x,ξ) ]其中Q表示给定第一阶段决策和风电偏差ξ后第二阶段的实时调整最小成本。这个式子看起来像个无底洞外面最小化里面在最坏分布上求期望再里面还要解一个实时调整问题。如果没有线性准则做桥三层结构无法直接丢给求解器。线性准则的意义就在于把第二层和第三层打碎重组。2.4 线性决策规则与对偶转化的关键步骤加入线性决策规则之后第二阶段调整量写成ξ的仿射函数。此时Q(x,u,ξ)本身成为ξ的分段线性凸函数。对Wasserstein模糊集专著里有经典结论sup_{P∈D} E_P[ h(ξ) ] 可以转化成下面的有限维问题min λ·ε (1/N) Σ_i sup_{ξ∈Ξ} [ h_i(ξ) - λ·d(ξ, ξ_i) ]其中λ是模糊集半径的对偶乘子h_i是第i个场景对应的第二阶段最优值。如果h是凸函数且ξ限制在多面体支集Ξ内内部那个sup可以进一步写成线性规划的对偶形式。两步对偶做完整个两阶段分布鲁棒问题就变成了一个确定性的MILP。这几步是模型里真正的hard part。代码实现上我的建议是先不急着把整套对偶公式敲进Matlab先把分支一步步拆开验证先写纯确定性版本再写随机场景版本最后再引入模糊集。每一步的求解结果都对得上再进下一步。我最早就是跳步直接写完整模型结果报了一堆“second argument must be a scalar”之类的错误排查几个小时后才发现是约束维度没对齐。3. Matlab代码实现的整体结构与关键模块3.1 程序架构设计我不会把整个工程直接往YALMIP里一丢了事那样后期调试非常痛苦。这里给出我实际采用的模块划分方案每个模块一个文件主程序只管拼装和求解。main_uc.m主程序定义系统参数、调用各模块、调用solver。gen_wind_scenarios.m生成风电场景库输入历史风速/功率数据输出场景样本。build_unit_data.m定义火电机组参数结构体包括容量、爬坡、成本、启停时间等。build_wasserstein_ambiguity.m构造Wasserstein模糊集半径和相关常数。build_dro_uc_model.m用YALMIP定义决策变量、目标函数和所有约束。solve_and_report.m调用optimize并整理输出启停计划、出力曲线、最坏分布信息。这套分层的好处是改机组参数不用动模型文件改模糊集构造不用动主程序排查约束问题可以直接在build_dro_uc_model里逐块注释验证。Matlab版本方面我用的是R2023a求解器用Gurobi 10.0前端建模用YALMIP。用CPLEX也能跑但Gurobi在大规模MILP上通常更快。安装配置不难把Gurobi安装目录下的matlab文件夹加到Matlab路径然后在Matlab里运行gurobi_setup脚本即可。3.2 风电场景库与模糊集参数生成风电场景相信大家手上都不缺公开数据集。如果没有合适的也可以从正态分布或ARMA模型合成但注意合成数据做鲁棒性评估会偏乐观最好还是用真实系统数据。生成场景库后每个场景是一组T维的风电出力向量。对这个场景库经验分布中心P_N就是一个离散均匀分布。模糊集半径ε的选取我的经验是先跑一次不带风雨场景的基线解作为参考再跑ε从0.05到0.5的扫描绘制“成本-鲁棒性”曲线。对6节点小系统来说ε0.15到0.25通常是比较实用的区间。% gen_wind_scenarios.m 示意 wind_raw load(wind_power_history.mat); % 历史风电数据 T 24; Nscen 200; scen zeros(T, Nscen); for k 1:Nscen id randi([1, length(wind_raw)], 1); scen(:,k) wind_raw(id:Tid-1); % 滚动抽取长度T的序列 end % 归一化到装机容量区间 [0, Wmax] scen scen / max(scen(:)) * Wmax;抽取场景时要注意序列相关性问题。风电出力有强时间相关性逐时段独立采样会破坏爬坡过程的真实性。所以虽然叫场景实际要按连续时间序列抽不能让第二天凌晨5点的风电和下午3点完全无关。3.3 YALMIP建模与求解流程在YALMIP里定义决策变量u binvar(T, I, full); % 启停状态 P sdpvar(T, I, full); % 基础出力 % 线性决策规则变量 a sdpvar(T, I, full); % 调整量常数项 b sdpvar(T, I, full); % 调整量对风电偏差的灵敏度系数 % 辅助对偶变量 lambda sdpvar(1, 1); % 模糊集半径对偶乘子 s sdpvar(Nscen, 1); % 每个场景的最坏函数值变量目标函数里先把基础成本写出来再把第二阶段最坏分布期望用lambda和s序列表达。核心约束包括功率平衡。%% 目标函数基础成本 Objective sum(sum( fuel_cost(P, u) )) sum(sum(start_cost(:,:).* u)) ... lambda * eps_radius sum(s)/Nscen; %% 功率平衡确定性部分 Constraints []; Constraints [Constraints, sum(P,2) Load - meanWindVec]; Constraints [Constraints, ... ...];这里需要特别提醒一点不要把风电预测误差直接等价于负荷。风电不确定性引入后功率平衡约束要拆成两部分预测值部分用确定性等式平衡偏差部分交给第二阶段调整变量处理。否则模型会重复计及风电基准平衡点就会偏移。场景相关约束需要写出每个场景的阶段二最优值h_i与lambda和s的关系。在YALMIP里可以用sdpvar构建每个场景的小优化约束但更高效的方式是利用对偶公式直接写出线性约束组。具体来说每个场景ξ_i有一个对应的第二阶段线性规划其对偶问题可以合并成一组线性矩阵不等式写入总约束。这一步是代码中最容易出错的区域。我自己第一次写的时候脑子里公式和YALMIP变量名对不上调试了好久。建议把Wasserstein对偶公式一行一行写在注释里每条约束上方标注它对应原文里的哪一项。求解端设置ops sdpsettings(solver, gurobi, verbose, 2, ... gurobi.MIPGap, 0.01, ... gurobi.TimeLimit, 1200, ... savesolveroutput, 1); result optimize(Constraints, Objective, ops);MIPGap给1%对小系统足够对大规模系统可以放宽到2%到3%。时间限制1200秒也是我常用的设定超过这个时间即使gap没到1%我也会接受当前最优解并在报告中标注。3.4 结果输出与可视化求解完之后我习惯做三件事第一画出各机组启停状态的热力图横轴时段纵轴机组编号暖色表示开机。这个图能一眼看出分布鲁棒解相比确定性解多开了哪些机组。第二画出风电实际可能的区间带用场景库的10%到90%分位数叠加负荷曲线和各机组出力堆叠图。这个图能直接反映系统在晚间低谷时段是否可能面临出力和负荷反向的困境。第三输出模糊集最坏分布的“典型场景”。虽然最坏分布本身是一群概率质量但可以把概率权重最大的几个场景挑出来展示分析它们对应的风电曲线有什么共性。我在多个算例里观察到的共性是最坏场景往往不是单一最大误差场景而是数个中等偏差场景的组合它们虽然单看都不算极端但会在相近时段出现持续性的偏差使系统备用在连续小时段内被逐步消耗。这一点在实际运行中很有意义连续偏差比单点尖峰更危险。%% 输出计划到excel result_table table((1:T), u(:,1), P(:,1), ... ); writetable(result_table, unit_commitment_result.xlsx);4. 仿真算例设计与结果解读4.1 测试系统参数我用一个修改版6节点系统做演示测试。系统含3台火电机组、1个风电场风电场装机容量取系统峰值负荷的25%。负荷曲线取典型冬季日负荷峰谷差约30%。机组参数刻意设置得有一定区分度一台大容量机组效率高但启停成本高适合做基荷两台小机组启停灵活但单位煤耗高适合做峰荷和爬坡调节。这种设置能让分布鲁棒解和确定性解的启停差异更明显方便观察不确定性对机组组合结构的影响。调度周期24小时时段间隔1小时。场景数取200模糊集半径ε0.2。对比基准设置三组确定性模型风电预测值、随机规划模型场景概率固定、分布鲁棒模型本文模型。4.2 三套方案的成本与启停差异方案基础发电成本相对最坏分布下期望调整成本总成本相对夜间最小开机台数确定性1.0000.1321.1322随机规划1.0180.0811.0993分布鲁棒1.0270.0451.0723数值是归一化后的演示结果具体绝对值会因为燃油价格和系统容量不同而变但趋势很稳定。确定性方案的基础发电成本最低因为它只盯着一个预测点做优化开机台数最少高成本机组基本不用开。但一旦实际风电和预测偏离第二阶段的调整成本很高总成本反而最贵。随机规划稍微多开了一台机组基础成本上升但调整成本下降不少。问题是它信奉场景概率在训练场景的分布上最优换一批历史数据表现就会打折。分布鲁棒方案基础成本最高多开了机组、预留给了一部分调节能力但最坏分布下的期望调整成本只有确定性的三分之一左右。总成本在三者中最低而且对分布偏移不敏感。这个结果很能说明问题机组组合做决策时多花一点前置成本换来的往往是后验总成本的大幅下降。这个现象在风电渗透率越高时越显著。把风电渗透率从25%调高到40%再跑一轮分布鲁棒的相对优势会进一步扩大因为高渗透率下预测误差对系统安全的影响是非线性放大的。4.3 最坏分布揭示的隐患场景我额外关注了最坏分布下概率权重最高的场景。在我这个演示算例里最坏场景集中在“夜间风电高估清晨风电急跌”的组合上。具体来说夜间负荷低谷时段风电预测出力偏高模型在确定性解里会少开机组。而最坏场景实际风速偏低风电出力只有预测的60%系统需要更多的向下调节能力已有的机组要么降负荷空间不够要么被迫启停频繁切换。清晨负荷爬坡叠加风电持续偏低备用紧张又集中在同一时段最终触发切负荷风险。分布鲁棒解正是因为提前看破了这种“连续偏差”的风险才会多开一台灵活机组这台机组白天看确实有点多余但它在凌晨时段的存在就是系统的安全垫。从工程角度讲这个解读比单纯看成本更有价值。因为实际调度中预测系统的误差模式往往不是随机的而是有系统性偏差的比如数值天气预报在特定天气过程中普遍高估风速。分布鲁棒优化通过对历史误差分布建模间接把这些偏差模式也纳入了决策考虑。5. 调试经验与常见坑速查5.1 模型一直提示不可行怎么办如果是小系统先用确定性模型跑通基线。确定性能跑通再逐步加入随机场景最后加入模糊集对偶约束。哪一个环节报不可行就在哪个环节排查。常用手段是给功率平衡约束加一个极小的松弛变量并给予很大惩罚系数。这样如果模型无解松弛变量会吃掉差额你就能从松弛变量的大小看出是哪个时段的哪个约束崩了。比如发现凌晨3点的功率平衡松弛变非零说明那台机组组合在3点的出力范围覆盖不了负荷需要放开某台机组的启停约束或增加备用。机组最小启停时间约束也是常见的不可行来源。晚上高峰开机凌晨低谷要求最小停机时间不足前一个时段的开机状态一直锁定到后半夜导致出力下不来。排查这种问题时我会把机组启停变量的整数值输出出来和最小启停时间约束一一对比看有没有违反的时段。5.2 大M参数怎么选机组组合模型里出力上下限、爬坡约束这些大M约束的大M如果取得太离谱MILP的LP松弛会很松散求解速度显著变慢。正确做法是大M取物理上限加一个裕度而不是图省事填一个100000。对机组出力约束大M取该机组的Pmax即可最多加一个小量缓冲。对备用容量约束取Pmax乘以1.2足够了。我在调试时遇到过一个case把备用相关约束的大M写成了所有机组容量之和导致求解时间从几分钟涨到两个多小时而且gap一直上不去。改回单机Pmax之后立刻恢复到正常求解速度。这个细节对大规模系统尤其重要。5.3 模糊集半径的敏感性与社会经济性权衡ε不是越大越好。半径从0.1增到0.5总成本往往单调上升但上升速度在某个区间会突然变快。我在测试中会用成本-ε曲线来判断合理的鲁棒水平曲线平缓段对应的ε说明增加鲁棒性成本很小可以放心用曲线陡升段说明超过这个鲁棒边界后成本代价会爆发式增加与实际需要不符。另外一条经验是ε的选择应与样本容量N联动。样本多时经验分布本身相对可信ε可以取小一点样本少时经验分布不可靠必须留更大的半径。有一些文献给出ε关于N的解析表达式但工程上我更推荐用hold-out方法做选择把历史数据切成训练集和验证集选一个在验证集上整体表现最好的ε。5.4 求解时间过长怎么办大系统上百台机组、几千个场景直接跑完整MILP确实吃力。我常用的降维手段按优先级排列场景聚合。用k-means等聚类方法把2000个场景压到200个保留代表性场景和权重风电不确定性刻画损失通常不大。缩短调度周期。先跑一个冬季典型日验证模型正确性和参数调节后再扩展到整月。平行计算。多个ε值可以并行跑Gurobi本身支持并发求解不同模型实例Matlab的parfor可以做到同时跑多个ε扫描。用热启动。Gurobi支持提供MIP start把确定性解的启停变量作为初始解塞进去能显著压缩分支树。如果以上都做了还是慢就该考虑把时段粒度从1小时放大到2小时或者砍刨负荷曲线的尖峰细节换取计算速度。5.5 YALMIP细节与环境变量YALMIP对Gurobi的版本要求比较严格Gurobi升级后有时会出现“No appropriate solver”之类的报错。这通常是因为YALMIP的gurobi接口检测不到新版本路径重跑gurobi_setup并确认Matlab的javaclasspath里有gurobi jar包通常能解决。另外在Matlab里定义sdpvar时维度一定要显式写清楚。我遇到过两次因为写成sdpvar(T,I)而不是sdpvar(I,T)导致所有矩阵运算都通过broadcast方式计算结果模型规模膨胀十倍以上的情况。这个小坑造成的隐性成本很惊人一旦怀疑模型大小异常第一步就是检查所有变量的维度方向。6. 后续扩展方向我做完这套模型后手上这个工程代码目前只解决了单一风电场的情况。后续想加的方向有几个一是把多风电场之间的出力相关性纳入模糊集构造这需要从copula或者空间相关矩阵入手二是把储能和需求响应作为第二阶段的调节手段建模进去这样线性决策规则的维度会变大但模型框架不用改三是把日前市场出清和实时平衡市场的价格信号放进来让分布鲁棒解直接对接经济调度和市场结算逻辑那就能真正从学术demo变成一个可以跟调度员对话的决策工具。我的经验是做这种优化模型最忌讳一上来就追求完全精确先跑通、再跑快、再跑准。把分布鲁棒优化机组组合理解为一个“用概率距离做保险”的决策工具箱而不是一个纯数学玩具它才能真正帮你在风电场旁边做出靠谱的调度决策。
基于线性准则的分布鲁棒优化机组组合建模与Matlab实现
做机组组合的人应该都听过一句话风电预测曲线和实际出力永远对不上。传统确定性方法在风电渗透率不高时还能凑合可当预测误差动辄几百兆瓦一个拍脑袋的备用容量可能让调度中心在凌晨三点焦头烂额。我最近用Matlab把一套基于线性准则的分布鲁棒优化机组组合模型完整实现了从数学推导到代码落地踩了不少坑。今天把整个思路、建模细节和程序实现串起来讲一遍包括风电不确定性建模、Wasserstein模糊集构造、两阶段决策规则以及用YALMIP调Gurobi求解的完整流程希望给正在做类似课题的人省点时间。1. 为什么机组组合必须认真对待风电不确定性1.1 确定性调度的局限在哪里机组组合本质上是一个时域耦合的大规模混合整数优化问题决定每台火电机组在各时段是开机还是停机并给出满足负荷与备用需求的有功出力计划。传统做法把风电当作已知量通常取预测值或者“预测值加一个固定比例备用”来打保险。风电渗透率在10%以内时这种近似还能接受。但渗透率超过30%之后问题就暴露了最典型的是两个场景一是预测值偏乐观时系统实际需要的上调备用不足。晚高峰负荷上来、风电突然掉出力备用机组要么还没热启动要么爬坡速度跟不上等着调度员的就只有切负荷。二是预测值偏保守时按高风电出力安排开机数量夜间实际出力却远低于预期部分机组被迫压到最小技术出力以下只能弃风甚至倒逼最小出力更高的机组停机再启动成本和启停次数一起飙升。固定备用系数的方法本质上是“赌误差不超过某个经验的百分数”。它不能回答一个基本问题极端场景到底会坏到什么程度以及为了应对这种程度需要付出多少成本。这就是随机优化、鲁棒优化和分布鲁棒优化登场的背景。1.2 三种不确定性优化路线怎么选先说随机规划。它给风电出力预设一个离散场景集每个场景有发生概率目标变成“全场景期望成本最小”。好处是决策有经济性坏处是你得相信场景概率是对真实分布的准确近似。概率设错了优化解很可能只在训练场景里好看换一批实际数据就翻车。鲁棒优化走另一个极端要求所有可能的风电曲线都在约束可行域内目标是最坏情况成本最小。它不需要概率分布只需要不确定性集合比如盒式区间。缺点是“最坏情况”常常是永远不会发生的极端情形为了保证这种极端系统得多开很多机组备用冗余高到经济性很差。分布鲁棒优化夹在两者中间。它假设风电出力的真实概率分布落在某个模糊集内这个模糊集以历史场景的经验分布为中心、以某个距离半径为界。优化的目标是在模糊集内的所有分布中找“最坏的那个分布”再对这最坏分布求期望成本最小。换句话说它对概率分布本身做鲁棒而不是对场景集合做鲁棒。既保留了经济性又不会因为一个极端点就过度设计近几年在电力系统调度文献里热度很高。这套思路落到机组组合上就非常合适机组启停计划是典型的“今天做决定、明天见真章”必须今天就在不知道明天风电分布真实样子的情况下做决策。分布鲁棒框架天然支持这种“现在做决定、事后看结果”的两阶段结构。1.3 标题里的“线性准则”到底指什么很多初次接触这个概念的人会被“分布鲁棒优化”几个字劝退其实它背后有一个很实用的工程近似思路第二阶段的可调决策被限制为风电出力偏差的仿射线性函数。举个例子第一阶段决定机组i在时段t是否开机并给一个基础出力。第二阶段风电实际出力和预测值出现偏差ξ之后机组会在基础出力附近增减一个调整量。这个调整量如果是直接由优化器自由决定两阶段问题会和概率分布耦合得很深几乎没法直接求解。线性准则把这个调整量写成关于ξ的线性函数比如adj_i,t(ξ) a_i,t b_i,t · ξ其中a是常数项b是灵敏度系数。这样原来那个嵌套的max-min-期望问题通过线性对偶变换可以转成一个确定性的混合整数线性规划。这就是“基于线性准则”的含义。有人可能会觉得线性函数限制太强风电偏差大时调整策略不一定是最优的。但从工程角度看线性决策规则带来的成本损失通常在5%以内换来的却是计算速度快几十倍这是非常划算的交易。我在实际测试中规模稍大的算例如果用完全自由的两阶段模型求解器经常一个小时内出不了可行解而线性准则版本基本几分钟内就能收敛到MIP gap可接受的水平。2. 分布鲁棒机组组合的数学模型与转化2.1 基础机组组合模型建模从经典的确定性机组组合出发。系统有I台火电机组调度周期T个时段每台机组有最小开机/停机时间、出力上下限、爬坡速率、启动与停机成本。目标函数是燃料成本加启停成本燃料成本用二次函数或者分段线性近似min Σ_t Σ_i [ f_i(P_i,t) · u_i,t SU_i,t SD_i,t ]约束至少包括功率平衡Σ_i P_i,t W_t^f D_t风电取预测出力时系统发电必须等于负荷。旋转备用Σ_i (P_i,max - P_i,t) · u_i,t ≥ R_t 风电备用需求。出力上下限P_i,min · u_i,t ≤ P_i,t ≤ P_i,max · u_i,t。爬坡约束P_i,t - P_i,t-1 ≤ RU_i以及下坡约束。最小启停时间约束用经典的MILP三不等式形式做线性化。启停成本约束许多模型里会把启动成本与停机时间挂钩。这套模型本身不复杂难点在风电不确定性的注入。2.2 风电不确定性的模糊集构造常见做法是用历史数据或场景法生成N个风电出力样本ξ_1,...,ξ_N构造经验分布P_N。然后定义以P_N为中心、以Wasserstein距离为度量的模糊集D { P : d_W(P, P_N) ≤ ε }这里的ε也称模糊集半径直观理解是真实分布和经验分布之间的“概率距离”上限。ε越小模糊集越小你越信任历史数据ε越大模糊集越大鲁棒性越强但保守性也越强。Wasserstein距离相比KL散度有个明显好处它允许分布的支撑集发生变化也就是说真实分布可以出现“历史样本里没见过的出力值”这在风电场景里非常重要因为极端低风和高风事件本来就稀疏历史样本集未必能覆盖。对模糊集做参数敏感性测试是必须的步骤。我在测试中发现ε取0.05时解基本和随机规划差不多取0.3时系统会明显多开机组夜间开机台数可能增加一到两台这就是鲁棒性代价的直观体现。选择合适ε的实用办法是拿一段真实历史风电数据做后验评估比较不同ε下解的期望运行成本与最坏运行成本。2.3 两阶段分布鲁棒模型怎么搭把风电不确定性引入后模型写成两阶段形式第一阶段确定机组启停计划u和基础出力P0。第二阶段风电实际场景ξ发生后机组可以调整出力也可以切负荷或弃风目标是在所有可能分布中找最坏分布下期望调整成本最小的方案。形式化一点可以写成min Σ成本 sup_{P∈D} E_P[ Q(P,x,ξ) ]其中Q表示给定第一阶段决策和风电偏差ξ后第二阶段的实时调整最小成本。这个式子看起来像个无底洞外面最小化里面在最坏分布上求期望再里面还要解一个实时调整问题。如果没有线性准则做桥三层结构无法直接丢给求解器。线性准则的意义就在于把第二层和第三层打碎重组。2.4 线性决策规则与对偶转化的关键步骤加入线性决策规则之后第二阶段调整量写成ξ的仿射函数。此时Q(x,u,ξ)本身成为ξ的分段线性凸函数。对Wasserstein模糊集专著里有经典结论sup_{P∈D} E_P[ h(ξ) ] 可以转化成下面的有限维问题min λ·ε (1/N) Σ_i sup_{ξ∈Ξ} [ h_i(ξ) - λ·d(ξ, ξ_i) ]其中λ是模糊集半径的对偶乘子h_i是第i个场景对应的第二阶段最优值。如果h是凸函数且ξ限制在多面体支集Ξ内内部那个sup可以进一步写成线性规划的对偶形式。两步对偶做完整个两阶段分布鲁棒问题就变成了一个确定性的MILP。这几步是模型里真正的hard part。代码实现上我的建议是先不急着把整套对偶公式敲进Matlab先把分支一步步拆开验证先写纯确定性版本再写随机场景版本最后再引入模糊集。每一步的求解结果都对得上再进下一步。我最早就是跳步直接写完整模型结果报了一堆“second argument must be a scalar”之类的错误排查几个小时后才发现是约束维度没对齐。3. Matlab代码实现的整体结构与关键模块3.1 程序架构设计我不会把整个工程直接往YALMIP里一丢了事那样后期调试非常痛苦。这里给出我实际采用的模块划分方案每个模块一个文件主程序只管拼装和求解。main_uc.m主程序定义系统参数、调用各模块、调用solver。gen_wind_scenarios.m生成风电场景库输入历史风速/功率数据输出场景样本。build_unit_data.m定义火电机组参数结构体包括容量、爬坡、成本、启停时间等。build_wasserstein_ambiguity.m构造Wasserstein模糊集半径和相关常数。build_dro_uc_model.m用YALMIP定义决策变量、目标函数和所有约束。solve_and_report.m调用optimize并整理输出启停计划、出力曲线、最坏分布信息。这套分层的好处是改机组参数不用动模型文件改模糊集构造不用动主程序排查约束问题可以直接在build_dro_uc_model里逐块注释验证。Matlab版本方面我用的是R2023a求解器用Gurobi 10.0前端建模用YALMIP。用CPLEX也能跑但Gurobi在大规模MILP上通常更快。安装配置不难把Gurobi安装目录下的matlab文件夹加到Matlab路径然后在Matlab里运行gurobi_setup脚本即可。3.2 风电场景库与模糊集参数生成风电场景相信大家手上都不缺公开数据集。如果没有合适的也可以从正态分布或ARMA模型合成但注意合成数据做鲁棒性评估会偏乐观最好还是用真实系统数据。生成场景库后每个场景是一组T维的风电出力向量。对这个场景库经验分布中心P_N就是一个离散均匀分布。模糊集半径ε的选取我的经验是先跑一次不带风雨场景的基线解作为参考再跑ε从0.05到0.5的扫描绘制“成本-鲁棒性”曲线。对6节点小系统来说ε0.15到0.25通常是比较实用的区间。% gen_wind_scenarios.m 示意 wind_raw load(wind_power_history.mat); % 历史风电数据 T 24; Nscen 200; scen zeros(T, Nscen); for k 1:Nscen id randi([1, length(wind_raw)], 1); scen(:,k) wind_raw(id:Tid-1); % 滚动抽取长度T的序列 end % 归一化到装机容量区间 [0, Wmax] scen scen / max(scen(:)) * Wmax;抽取场景时要注意序列相关性问题。风电出力有强时间相关性逐时段独立采样会破坏爬坡过程的真实性。所以虽然叫场景实际要按连续时间序列抽不能让第二天凌晨5点的风电和下午3点完全无关。3.3 YALMIP建模与求解流程在YALMIP里定义决策变量u binvar(T, I, full); % 启停状态 P sdpvar(T, I, full); % 基础出力 % 线性决策规则变量 a sdpvar(T, I, full); % 调整量常数项 b sdpvar(T, I, full); % 调整量对风电偏差的灵敏度系数 % 辅助对偶变量 lambda sdpvar(1, 1); % 模糊集半径对偶乘子 s sdpvar(Nscen, 1); % 每个场景的最坏函数值变量目标函数里先把基础成本写出来再把第二阶段最坏分布期望用lambda和s序列表达。核心约束包括功率平衡。%% 目标函数基础成本 Objective sum(sum( fuel_cost(P, u) )) sum(sum(start_cost(:,:).* u)) ... lambda * eps_radius sum(s)/Nscen; %% 功率平衡确定性部分 Constraints []; Constraints [Constraints, sum(P,2) Load - meanWindVec]; Constraints [Constraints, ... ...];这里需要特别提醒一点不要把风电预测误差直接等价于负荷。风电不确定性引入后功率平衡约束要拆成两部分预测值部分用确定性等式平衡偏差部分交给第二阶段调整变量处理。否则模型会重复计及风电基准平衡点就会偏移。场景相关约束需要写出每个场景的阶段二最优值h_i与lambda和s的关系。在YALMIP里可以用sdpvar构建每个场景的小优化约束但更高效的方式是利用对偶公式直接写出线性约束组。具体来说每个场景ξ_i有一个对应的第二阶段线性规划其对偶问题可以合并成一组线性矩阵不等式写入总约束。这一步是代码中最容易出错的区域。我自己第一次写的时候脑子里公式和YALMIP变量名对不上调试了好久。建议把Wasserstein对偶公式一行一行写在注释里每条约束上方标注它对应原文里的哪一项。求解端设置ops sdpsettings(solver, gurobi, verbose, 2, ... gurobi.MIPGap, 0.01, ... gurobi.TimeLimit, 1200, ... savesolveroutput, 1); result optimize(Constraints, Objective, ops);MIPGap给1%对小系统足够对大规模系统可以放宽到2%到3%。时间限制1200秒也是我常用的设定超过这个时间即使gap没到1%我也会接受当前最优解并在报告中标注。3.4 结果输出与可视化求解完之后我习惯做三件事第一画出各机组启停状态的热力图横轴时段纵轴机组编号暖色表示开机。这个图能一眼看出分布鲁棒解相比确定性解多开了哪些机组。第二画出风电实际可能的区间带用场景库的10%到90%分位数叠加负荷曲线和各机组出力堆叠图。这个图能直接反映系统在晚间低谷时段是否可能面临出力和负荷反向的困境。第三输出模糊集最坏分布的“典型场景”。虽然最坏分布本身是一群概率质量但可以把概率权重最大的几个场景挑出来展示分析它们对应的风电曲线有什么共性。我在多个算例里观察到的共性是最坏场景往往不是单一最大误差场景而是数个中等偏差场景的组合它们虽然单看都不算极端但会在相近时段出现持续性的偏差使系统备用在连续小时段内被逐步消耗。这一点在实际运行中很有意义连续偏差比单点尖峰更危险。%% 输出计划到excel result_table table((1:T), u(:,1), P(:,1), ... ); writetable(result_table, unit_commitment_result.xlsx);4. 仿真算例设计与结果解读4.1 测试系统参数我用一个修改版6节点系统做演示测试。系统含3台火电机组、1个风电场风电场装机容量取系统峰值负荷的25%。负荷曲线取典型冬季日负荷峰谷差约30%。机组参数刻意设置得有一定区分度一台大容量机组效率高但启停成本高适合做基荷两台小机组启停灵活但单位煤耗高适合做峰荷和爬坡调节。这种设置能让分布鲁棒解和确定性解的启停差异更明显方便观察不确定性对机组组合结构的影响。调度周期24小时时段间隔1小时。场景数取200模糊集半径ε0.2。对比基准设置三组确定性模型风电预测值、随机规划模型场景概率固定、分布鲁棒模型本文模型。4.2 三套方案的成本与启停差异方案基础发电成本相对最坏分布下期望调整成本总成本相对夜间最小开机台数确定性1.0000.1321.1322随机规划1.0180.0811.0993分布鲁棒1.0270.0451.0723数值是归一化后的演示结果具体绝对值会因为燃油价格和系统容量不同而变但趋势很稳定。确定性方案的基础发电成本最低因为它只盯着一个预测点做优化开机台数最少高成本机组基本不用开。但一旦实际风电和预测偏离第二阶段的调整成本很高总成本反而最贵。随机规划稍微多开了一台机组基础成本上升但调整成本下降不少。问题是它信奉场景概率在训练场景的分布上最优换一批历史数据表现就会打折。分布鲁棒方案基础成本最高多开了机组、预留给了一部分调节能力但最坏分布下的期望调整成本只有确定性的三分之一左右。总成本在三者中最低而且对分布偏移不敏感。这个结果很能说明问题机组组合做决策时多花一点前置成本换来的往往是后验总成本的大幅下降。这个现象在风电渗透率越高时越显著。把风电渗透率从25%调高到40%再跑一轮分布鲁棒的相对优势会进一步扩大因为高渗透率下预测误差对系统安全的影响是非线性放大的。4.3 最坏分布揭示的隐患场景我额外关注了最坏分布下概率权重最高的场景。在我这个演示算例里最坏场景集中在“夜间风电高估清晨风电急跌”的组合上。具体来说夜间负荷低谷时段风电预测出力偏高模型在确定性解里会少开机组。而最坏场景实际风速偏低风电出力只有预测的60%系统需要更多的向下调节能力已有的机组要么降负荷空间不够要么被迫启停频繁切换。清晨负荷爬坡叠加风电持续偏低备用紧张又集中在同一时段最终触发切负荷风险。分布鲁棒解正是因为提前看破了这种“连续偏差”的风险才会多开一台灵活机组这台机组白天看确实有点多余但它在凌晨时段的存在就是系统的安全垫。从工程角度讲这个解读比单纯看成本更有价值。因为实际调度中预测系统的误差模式往往不是随机的而是有系统性偏差的比如数值天气预报在特定天气过程中普遍高估风速。分布鲁棒优化通过对历史误差分布建模间接把这些偏差模式也纳入了决策考虑。5. 调试经验与常见坑速查5.1 模型一直提示不可行怎么办如果是小系统先用确定性模型跑通基线。确定性能跑通再逐步加入随机场景最后加入模糊集对偶约束。哪一个环节报不可行就在哪个环节排查。常用手段是给功率平衡约束加一个极小的松弛变量并给予很大惩罚系数。这样如果模型无解松弛变量会吃掉差额你就能从松弛变量的大小看出是哪个时段的哪个约束崩了。比如发现凌晨3点的功率平衡松弛变非零说明那台机组组合在3点的出力范围覆盖不了负荷需要放开某台机组的启停约束或增加备用。机组最小启停时间约束也是常见的不可行来源。晚上高峰开机凌晨低谷要求最小停机时间不足前一个时段的开机状态一直锁定到后半夜导致出力下不来。排查这种问题时我会把机组启停变量的整数值输出出来和最小启停时间约束一一对比看有没有违反的时段。5.2 大M参数怎么选机组组合模型里出力上下限、爬坡约束这些大M约束的大M如果取得太离谱MILP的LP松弛会很松散求解速度显著变慢。正确做法是大M取物理上限加一个裕度而不是图省事填一个100000。对机组出力约束大M取该机组的Pmax即可最多加一个小量缓冲。对备用容量约束取Pmax乘以1.2足够了。我在调试时遇到过一个case把备用相关约束的大M写成了所有机组容量之和导致求解时间从几分钟涨到两个多小时而且gap一直上不去。改回单机Pmax之后立刻恢复到正常求解速度。这个细节对大规模系统尤其重要。5.3 模糊集半径的敏感性与社会经济性权衡ε不是越大越好。半径从0.1增到0.5总成本往往单调上升但上升速度在某个区间会突然变快。我在测试中会用成本-ε曲线来判断合理的鲁棒水平曲线平缓段对应的ε说明增加鲁棒性成本很小可以放心用曲线陡升段说明超过这个鲁棒边界后成本代价会爆发式增加与实际需要不符。另外一条经验是ε的选择应与样本容量N联动。样本多时经验分布本身相对可信ε可以取小一点样本少时经验分布不可靠必须留更大的半径。有一些文献给出ε关于N的解析表达式但工程上我更推荐用hold-out方法做选择把历史数据切成训练集和验证集选一个在验证集上整体表现最好的ε。5.4 求解时间过长怎么办大系统上百台机组、几千个场景直接跑完整MILP确实吃力。我常用的降维手段按优先级排列场景聚合。用k-means等聚类方法把2000个场景压到200个保留代表性场景和权重风电不确定性刻画损失通常不大。缩短调度周期。先跑一个冬季典型日验证模型正确性和参数调节后再扩展到整月。平行计算。多个ε值可以并行跑Gurobi本身支持并发求解不同模型实例Matlab的parfor可以做到同时跑多个ε扫描。用热启动。Gurobi支持提供MIP start把确定性解的启停变量作为初始解塞进去能显著压缩分支树。如果以上都做了还是慢就该考虑把时段粒度从1小时放大到2小时或者砍刨负荷曲线的尖峰细节换取计算速度。5.5 YALMIP细节与环境变量YALMIP对Gurobi的版本要求比较严格Gurobi升级后有时会出现“No appropriate solver”之类的报错。这通常是因为YALMIP的gurobi接口检测不到新版本路径重跑gurobi_setup并确认Matlab的javaclasspath里有gurobi jar包通常能解决。另外在Matlab里定义sdpvar时维度一定要显式写清楚。我遇到过两次因为写成sdpvar(T,I)而不是sdpvar(I,T)导致所有矩阵运算都通过broadcast方式计算结果模型规模膨胀十倍以上的情况。这个小坑造成的隐性成本很惊人一旦怀疑模型大小异常第一步就是检查所有变量的维度方向。6. 后续扩展方向我做完这套模型后手上这个工程代码目前只解决了单一风电场的情况。后续想加的方向有几个一是把多风电场之间的出力相关性纳入模糊集构造这需要从copula或者空间相关矩阵入手二是把储能和需求响应作为第二阶段的调节手段建模进去这样线性决策规则的维度会变大但模型框架不用改三是把日前市场出清和实时平衡市场的价格信号放进来让分布鲁棒解直接对接经济调度和市场结算逻辑那就能真正从学术demo变成一个可以跟调度员对话的决策工具。我的经验是做这种优化模型最忌讳一上来就追求完全精确先跑通、再跑快、再跑准。把分布鲁棒优化机组组合理解为一个“用概率距离做保险”的决策工具箱而不是一个纯数学玩具它才能真正帮你在风电场旁边做出靠谱的调度决策。