☰
基于线性准则的分布鲁棒优化机组组合建模与Matlab实现
2026/10/6 5:15:49 网站建设 项目流程

做机组组合的人应该都听过一句话:风电预测曲线和实际出力永远对不上。传统确定性方法在风电渗透率不高时还能凑合,可当预测误差动辄几百兆瓦,一个拍脑袋的备用容量可能让调度中心在凌晨三点焦头烂额。我最近用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 ]

约束至少包括:

  1. 功率平衡:Σ_i P_i,t + W_t^f = D_t,风电取预测出力时系统发电必须等于负荷。
  2. 旋转备用:Σ_i (P_i,max - P_i,t) · u_i,t ≥ R_t + 风电备用需求。
  3. 出力上下限:P_i,min · u_i,t ≤ P_i,t ≤ P_i,max · u_i,t。
  4. 爬坡约束:P_i,t - P_i,t-1 ≤ RU_i,以及下坡约束。
  5. 最小启停时间约束:用经典的MILP三不等式形式做线性化。
  6. 启停成本约束,许多模型里会把启动成本与停机时间挂钩。

这套模型本身不复杂,难点在风电不确定性的注入。

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:T+id-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确实吃力。我常用的降维手段按优先级排列:

  1. 场景聚合。用k-means等聚类方法把2000个场景压到200个,保留代表性场景和权重,风电不确定性刻画损失通常不大。
  2. 缩短调度周期。先跑一个冬季典型日,验证模型正确性和参数调节后再扩展到整月。
  3. 平行计算。多个ε值可以并行跑,Gurobi本身支持并发求解不同模型实例,Matlab的parfor可以做到同时跑多个ε扫描。
  4. 用热启动。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变成一个可以跟调度员对话的决策工具。

我的经验是:做这种优化模型最忌讳一上来就追求完全精确,先跑通、再跑快、再跑准。把分布鲁棒优化机组组合理解为一个“用概率距离做保险”的决策工具箱,而不是一个纯数学玩具,它才能真正帮你在风电场旁边做出靠谱的调度决策。

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询