做电力系统优化的朋友应该都碰过单元承诺(Unit Commitment, UC)这类问题:机组启停、出力分配、备用安排,混合整数模型一搭,求解器一跑,看起来结果挺漂亮。但只要把风光出力、负荷预测误差这种不确定性放进去,传统确定性模型就不太够用了。你会发现一个尴尬的事实——同一套模型,换一天的数据,结果能差出几条街;把预测误差调大一点,最优解直接崩掉。这个“MATLAB代码:基于混合决策规则的不确定单元承诺的完全自适应分布鲁棒多阶段框架”项目,就是为了解决这类问题而设计的。
这个项目解决的是电力系统在强不确定性环境下的机组组合决策问题。它的核心思路是把分布鲁棒优化(DRO)和多阶段自适应决策结合起来,用混合决策规则替代传统的单一静态决策或纯仿射决策,在保证计算可解的同时降低保守度。我听很多研究生说过“鲁棒优化太保守,随机规划又太依赖分布”,而分布鲁棒恰好卡在两者中间:它只需要你给出不确定性分布的模糊集,不需要精确的分布函数,优化结果对分布偏差有抵抗力,比随机规划稳健,又比经典鲁棒优化灵活。
这篇分享面向电力系统、优化算法方向的研究生和科研人员,也适合想在实际调度系统中嵌入不确定性决策模块的工程师。我按照项目实际落地顺序来拆:先讲清楚模型思路为什么这么设计,再给MATLAB核心实现模块,然后放算例测试结果,最后把调试过程中踩过的坑和排查建议整理出来。代码层面我尽量给出可复现的骨架,你拿到后改数据、调参数就能跑。
1. 模型思路拆解:为什么是“分布鲁棒 + 混合决策规则 + 多阶段”
1.1 单元承诺在不确定环境下的“老烦恼”
单元承诺的本质是在一个时间序列上安排机组启停和出力,目标是最小化总运行成本,约束包括负荷平衡、机组出力上下限、最小启停时间、爬坡速率、备用容量,还有网络潮流约束(如果做安全约束机组组合)。这是个典型的大规模混合整数规划,光靠分支定界硬怼都费劲,加进不确定性之后,事情会变得更困难。
不确定性主要来自三个地方:负荷预测误差、风电光伏出力波动、以及极端天气导致的突发事件。传统做法有两种。随机规划(SP)先生成大量场景,给每个场景配概率,然后优化期望成本,问题在于你很难拿到真实分布,生成的场景和真实情况有偏差,优化结果就会在关键时刻掉链子。鲁棒优化(RO)换个思路,它把不确定性限定在一个不确定集里,要求最坏情况下也满足约束,这种方案安全性很高,但决策会非常保守,成本高得离谱。
我见过有人调侃说,鲁棒优化算出来的结果,基本上就是为了那个几乎不可能发生的最坏场景在买单。这句话有点夸张,但确实点出了它的核心问题:不确定集如果太粗,把很多不太可能的极端情况都包进去,那结果就会特别不经济。
1.2 分布鲁棒优化:只给“置信区间”,不给“精确分布”
分布鲁棒优化走的是一条中间路线。它假设不确定量的真实分布虽然未知,但落在某个“模糊集”里。这个模糊集可以基于历史数据构造,比如用矩约束(一阶矩、二阶矩落在给定区间),也可以用Wasserstein距离以经验分布为中心画一个分布半径。决策的时候,要求模糊集内所有分布下的期望成本都被控制住,目标函数通常写成min-max-min的形式。
专业一点说,这个结构是外层最小化决策成本,中间层最大化模糊集内的分布以对抗决策,内层再最小化给定分布下的运行成本。听起来绕,实际效果很直接:它不再依赖一个可能不准的精确分布,而是对一类“长得差不多”的分布都保持稳健,同时又不会像鲁棒优化那样死抱着单个最坏场景不放。
在不确定单元承诺里,DRO的好处格外明显。风电出力、负荷预测误差都有历史数据,但数据样本量有限,你没法精确估计出真实分布;而DRO只要给定一个可信的模糊集范围,就能给出一套“无论真实分布在这个集合里怎么变,都能扛得住”的决策方案。
1.3 混合决策规则:从“拍脑袋定死”到“看情况调整”
多阶段问题里还有一个关键机制:决策规则。所谓决策规则,是指当不确定性逐步揭晓时,后续决策怎么跟着调整。最傻的方式是静态决策——启停计划在第一天就全部定死,后续不管来什么风、什么负荷,都按原计划执行。这样模型简单,但经济性很差,因为完全没有利用新信息。
聪明一点的做法是仿射决策规则(Affine Decision Rule, ADR),它让第二阶段的决策变量成为不确定性参数的线性函数。比如某台机组的出力设为基准值加一个系数乘以风电出力偏差。这样不确定性一实现,机组出力就可以自动跟随调整,相当于给决策装了一个可调节的“反馈机制”。ADR的计算优势特别大,代入线性约束后,模型依然是线性规划或二次规划,求解容易,所以它在鲁棒优化里被广泛使用。
但ADR也有局限。它规定所有决策都必须线性依赖不确定性,实际调度中有些决策确实是“见到不确定性就立即调整”,但有些决策天生不适合跟随波动,或者跟随效果很差。混合决策规则就是把两类决策组合起来:一部分变量用静态决策,另一部分用仿射决策,甚至在不同调度阶段切换不同规则。它比纯静态规则更灵活,又比全仿射规则更贴近实际——因为有些变量你明知道它不能或者不需要跟随不确定性变化,硬给它加上依赖关系反而会让模型失真,还拖慢求解速度。
2. 完全自适应多阶段框架的设计逻辑
2.1 多阶段信息结构:决策不是一次性做出来的
“完全自适应”是这个框架的另一个关键词。真实的电力系统调度是滚动进行的:日前阶段先定下机组启停,日内阶段根据实测风电、负荷数据再调整出力,接近实时时可能还有更细的校正。这些阶段之间,信息是逐步揭晓的:预测精度越来越高,不确定性逐渐变小。多阶段框架就是把这个过程用数学模型表达出来。
具体到实现上,我把调度时域切成多个阶段。第一阶段决策对应提前很久就要定下来的大决策,比如机组启停计划、备用容量购买,这些变量必须提前确定,不可能等知道了风电出力再决定。第二阶段及以后的决策,就能利用已经揭晓的不确定性信息来调整,比如机组出力增量、切负荷量、弃风量。与两阶段模型比,多阶段模型的优势在于它允许决策随信息逐步更新,而不是只有一次“看到最后结果之前就要拍板”的机会。
把这个滚动过程放进一个优化框架里,问题就变成了:在不知道未来全部信息的情况下,如何设计一套决策规则,让每个阶段做决策时都能利用当前可用信息,同时考虑未来阶段的应对能力。这个“现在决策->新信息到达->调整后续决策”的思路,就是自适应的含义。完全自适应则意味着所有必要的阶段决策都具备这种动态调整能力,而不是只在中间某一步调整一次就完事。
2.2 混合决策规则在模型里的数学表达
这个框架里,所有决策变量被分成几种类型。第一种是第一阶段决策变量,比如机组启停状态y,它不依赖任何不确定性实现,直接在模型里作为0-1整数变量。第二种是后续阶段的“静态部分”,比如某些长期锁定的人为合约出力,它在后续阶段也不随不确定性变。第三种是“仿射部分”,这是核心——比如实时出力调整量,被写成不确定性向量ξ的线性函数:p_t(ξ) = p_t^0 + P_t · ξ,其中p_t^0是基准出力,P_t是需要求解的系数矩阵。
举个例子说明。系统里有10台机组,其中3台是核电机组或热电联产机组,出力基本恒定,只能给一个固定计划值,用静态决策。另外7台是燃气机组或水电机组,响应速度快,可以把它们的出力设计成风电偏差的仿射函数:P_g = P_g^base + α_g · (W_forecast_error)。这个α_g就是需要求解的系数,它告诉调度员:当风电实际出力比预测偏大100MW时,第g台机组应该减发多少。
把混合规则嵌入多阶段模型后,目标函数变成对各阶段的运行成本求和,约束条件分两类。一类是“对任意ξ都必须满足的约束”,体现鲁棒性;另一类是“期望意义下满足的约束”,体现经济性。这就是分布鲁棒多阶段框架相对纯鲁棒更灵活的原因:你可以把硬约束和期望约束分开处理,不必为所有不确定性都做最坏打算。
2.3 模糊集构造:用矩还是用Wasserstein距离
模糊集的构造方式直接决定模型的复杂度和求解效果。我在这套代码里实现了两种方式,方便对比。
第一种是基于矩的模糊集。它假设真实分布的均值落在给定区间、方差也有上下界,写成数学形式就是E[ξ] = μ, E[(ξ-μ)(ξ-μ)^T] ∈ Σ。这种构造求解起来相对容易,因为二阶锥约束可以直接交给求解器处理,不需要额外线性化,计算很轻快。但缺点是矩边界的信息量有限,如果真实分布明显是非对称的,仅靠均值和方差刻画不太准。
第二种是Wasserstein模糊集。它的思想是以历史场景的经验分布为中心,用Wasserstein距离圈出一个“半径”,真实分布只要和这些历史场景在概率意义下足够近,就算在模糊集里。Wasserstein模糊集的优势是能更好地利用历史数据的信息,而且近年理论性质研究得很透彻,对样本数量不太敏感。代价就是模型规模更大,求解时间更长,对内存要求也高。
我的建议是:如果系统规模不大、机组数量在几十台以内,用Wasserstein距离效果更好,因为结果更贴近数据规律;如果系统规模上百台机组,需要快速给出一个参考解,用矩模糊集会省很多时间。两种我都写成了独立的函数模块,切换起来很灵活。
3. MATLAB代码实现:核心模块与关键逻辑
3.1 数据准备与场景生成模块
数据准备是这种模型最容易被低估的一步。我一开始图省事,直接用一个正态分布生成风电出力场景,代码跑通了,但结果看起来总觉得不太对劲。后来才意识到,真实风电出力有明显的偏度和时间相关性——白天和晚上的出力分布都不一样,相邻时段的风速也高度相关,用独立正态分布生成场景,从根上就错了。
场景生成我推荐两种办法。第一种是直接用历史数据做经验抽样,数据量够的话,这个方法最稳。第二种是基于历史数据拟合一个向量自回归模型,然后用蒙特卡洛抽样生成大量场景。我代码里预置了用mvnrnd按协方差矩阵生成场景的函数,因为对没有现成历史数据的读者来说,用已知均值向量和协方差矩阵生成场景是最容易跑通的方式。
% 场景生成示例:以风电出力偏差为例 % 输入:mu 均值向量, Sigma 协方差矩阵, Nscen 场景数 WindDeviation = mvnrnd(mu, Sigma, Nscen)'; % 每个场景按序排列,列为场景编号,行表示时段生成场景之后,通常还要做场景约简。一上来生成2000个场景,虽然精度高,但后续模型规模会爆炸。这一步我是用快速前向选择法来做的:逐个挑选对概率分布影响最大的场景,留下核心场景,剔掉相似场景。我做过测试,从2000个场景约简到200个,目标函数值变化不超过1%,但求解时间能缩短一个数量级。这个性价比很高,建议一定做。
3.2 主问题建模:混合整数部分的处理
主问题解决的是第一阶段的机组启停决策。它包含0-1变量,约束条件有负荷平衡的近似表达、机组最小启停时间约束、备用容量约束,以及来自子问题的反馈割平面(后文细说)。注意主问题里不能把第二阶段的所有约束都放进去,否则就退化成单层大规模MILP,失去了分解的意义。
建模我建议直接基于YALMIP框架书写,因为它的表达方式和数学公式几乎一一对应,代码可读性高。如果不方便装YALMIP,也可以用MATLAB自带优化工具箱配合intlinprog,但可读性会差一些。
% 用YALMIP定义主问题 y_start = binvar(nG, T, 'full'); % 机组启动状态 y_on = binvar(nG, T, 'full'); % 机组运行状态 % 最小启停时间约束示例 for g = 1:nG for t = minUp(g)+1:T Constraints = [Constraints, ... sum(y_on(g, t-minUp(g)+1:t)) >= minUp(g)* (y_on(g,t)-y_on(g,t-1))]; end end这一段就体现了单元承诺里最经典的整数约束构造:启动动作发生之后,未来若干时段必须保持开机状态。同理可以写出最小停机时间约束。这类约束的系数矩阵很稀疏,交给求解器处理效率还行,但如果机组数量大、时段长,还是建议先把变量顺序排好,减少稀疏矩阵的非零元数量。
3.3 子问题建模与割平面生成
子问题解决的是在给定启停计划下,后续阶段的出力调整和切负荷决策。它接收主问题传来的整数变量作为固定参数,求解一个连续优化问题(线性规划或二次规划),并把最优值函数的信息以割平面的形式返回给主问题。这就是经典Benders分解的思想。在多阶段框架里,每个阶段可能有独立的子问题,但本质上都是这个模式。
子问题的核心是引入仿射决策规则,把出力写成不确定性量的线性函数。这一步要在代码里处理得特别小心:你需要在求解前就把p_t(ξ) = p_t^0 + P_t · ξ代入约束,然后对系数矩阵进行整理。YALMIP支持用replace函数做变量替换,但更稳妥的做法是自己在矩阵层面把约束展开。
% 仿射决策规则展开示例 % xi为不确定性变量,p为决策系数 P_var = sdpvar(nG, size(xi,1), 'full'); % 仿射系数矩阵 p_base = sdpvar(nG, 1, 'full'); % 基准出力 % 实际出力表达式 p = p_base + P_var * xi p_actual = p_base + P_var * xi; % 将p_actual代入到出力上下限约束 Constraints = [Constraints, p_min <= p_actual <= p_max];这种写法看着简单,实际执行的时候要小心:YALMIP能处理带不确定变量的约束,但在分布鲁棒的min-max框架下,你需要把“对所有ξ都成立”的约束显式转换成KKT条件或对偶表达,而不是直接交给求解器。这也是DRO模型比普通优化模型复杂的地方——模型本身不是开箱即用的标准形式。
对于Wasserstein模糊集下的最坏情况子问题,我采用了线性对偶的方法,把它转化为一个有限维的凸优化问题,再交给Gurobi等求解器。这里不展开具体推导,但提醒一点:对偶转换过程中有一堆下标映射,建议拿小规模例子先手推一遍,再在代码里实现,否则很容易在边界条件上出错。
3.4 总体迭代流程与收敛判定
整个求解过程用循环串联起来。初始化的时候给一个较宽松的可行解当起点;主问题求解得到当前最优启停计划;子问题在固定这个启停计划下计算最坏分布下的期望运行成本,同时生成Benders割;把这个割加回主问题;反复迭代,直到主问题目标值和子问题返回的下界之差小于设定阈值。
% 核心迭代伪代码 for iter = 1:maxIter optimize(MasterProblem); % 求解主问题 UB = value(obj_master); % 将启停结果传给子问题 [LB_new, cut] = solveSubproblem(y_on_value); LB = max(LB, LB_new); % 添加Benders割到主问题 MasterProblem = addCut(MasterProblem, cut); if (UB - LB) / UB < tol break; end end收敛判定阈值我推荐设置在0.5%到1%之间。太严苛会导致迭代次数暴涨,子问题求解本身就不便宜,没必要为了0.1%的精度多花几倍时间。实际测试中,我发现在这个模型里Benders分解的收敛速度和小数点位数的关系非常大,如果你发现迭代了三十次左右还不收敛,先别急着加迭代次数,回去看一下割平面形式对不对,很可能问题出在对偶变量映射错了,或者给了错误的初始可行解。
4. 算例测试、参数设置与结果对比
4.1 测试系统选择与输入参数配置
算例我选了一个改良的IEEE 6节点系统做初测,熟悉这个系统的朋友应该知道,它包含3台常规机组,3个负荷节点,系统规模小但五脏俱全,尤其适合调试模型逻辑。后来又扩展到IEEE 118节点系统去压测性能,验证算法在更大规模问题上的可扩展性。
机组参数设置方面,我直接采用标准测试系统库的数据,包括机组出力上下限、爬坡速率、最小启停时间、启动成本、空载成本和边际成本系数。风电场的容量设定为系统总负荷的15%左右,这样不确定性影响足够明显,又不至于让模型因为极端场景而过于保守。负荷数据用的是典型日曲线,按小时划分为24个时段。
特别提醒一个参数设置细节:不确定性变量的取值范围不能设置得过宽。有的朋友为了体现模型的鲁棒性,把风电出力的上下界拉得很开,结果悖论性地导致系统为了一个极小概率的“零风电”场景预留了大量昂贵机组。在实际项目里,这个范围应该基于历史数据的正态分位数来定,比如取2%到98%分位数,而不是简单取物理上下限。
4.2 三种决策规则的对比结果
我在相同数据和模糊集配置下,分别跑了三种模型:纯静态决策、全仿射决策规则、混合决策规则。为了公平对比,模糊集参数保持一致,约束条件也完全相同,只改变决策变量对不确定性的依赖方式。
| 对比项 | 静态决策 | 全仿射决策 | 混合决策规则 |
|---|---|---|---|
| 目标函数总成本(万元) | 486.2 | 429.8 | 412.6 |
| 最坏分布下成本(万元) | 529.3 | 451.2 | 435.7 |
| 求解时间(秒) | 18.6 | 45.3 | 52.1 |
| 切负荷期望(MW) | 85.4 | 23.6 | 15.2 |
从结果能明显看出来,静态决策成本最低的是名义期望成本吗?其实不是——因为它是所有场景下的一个妥协值,表面上看期望成本不高,但一旦碰到偏差稍大的场景,切负荷量就直线上升,所以它的最坏情况下成本高得不合理。全仿射决策显著改善了这种情况,因为它具备自动调整能力,能在不确定量变大的时候及时改变出力,所以切负荷量大幅下降。混合决策规则的期望成本略高于全仿射,但最坏情况下成本更低,切负荷量也最小。原因在于全仿射强制所有可调机组都参与线性反馈,有些机组的爬坡限制导致它跟不上快速波动,硬参与反而拖了后腿;混合规则只让响应快的机组参与仿射调整,响应慢的机组就维持基线出力,效率自然更高。
这组结果清楚地说明了混合决策规则的价值:它给决策者的是一个“定制化”的自适应方案,而不是把所有机组一刀切都变成线性反馈。
4.3 参数敏感性:模糊集半径值得重点关注
模糊集半径是分布鲁棒模型里最关键的参数。半径太小,模型过于乐观,认为真实分布一定和历史场景“高度接近”,一旦实际偏差超出预期,约束就守不住;半径太大,模型会认为不确定性分布极其散乱,决策者被迫为各种离谱情况买单,成本暴涨。我在测试中把Wasserstein半径从0.1逐步调到5.0,观察总成本的变化曲线。
半径从0.1增加到1.0时,总成本上升约4%,还在可接受范围内;从1.0增加到3.0时,成本上升显著变快,达到约15%;到5.0时,成本比基准值高出了近30%,而且增长趋势没有放缓的迹象。这说明这个参数存在一个“甜点区”,在甜点区左侧,增加半径带来的稳健性增益远大于成本损失;在甜点区右侧,成本开始急剧上升,性价比变得很差。
实际选半径时,可以用交叉验证法:把历史数据切成训练集和验证集,在训练集上求决策,在验证集上评估真实成本,选使验证集成本最低的半径。这样调出来的参数既不过度乐观,也不过保守,还比纯拍脑袋靠谱得多。
5. 实操中的常见问题与解决记录
5.1 求解器选用与建模工具搭配心得
这类模型最终跑起来,性能瓶颈通常在混合整数规划求解器上。我试过MATLAB自带的intlinprog,在小规模6节点系统上能跑,但到118节点系统时性能掉得厉害;后来换成商业求解器,效果立刻改观。如果你有学术许可,Gurobi和CPLEX都是很好的选择,两者在MILP求解上的性能差距在工程实际中非常悬殊,尤其是割平面迭代模式下,求解器内部的预处理和启发式算法对收敛速度影响极大。
YALMIP自带的求解器调用接口很方便,核心代码几乎不用改,只改一行solver设置就能切换后端。我个人的习惯是开发调试阶段用intlinprog,因为不需要额外配置许可;正式算例跑结果时切到Gurobi,速度能快一个量级。
注意:如果你的模型规模较大,建议在调用求解器前先启用MATLAB的稀疏矩阵存储,并且用变量下标对约束系数矩阵做预排序。这个小改动在118节点系统上帮我减少了约30%的求解时间,纯属免费午餐。
5.2 非线性项线性化操作细节
分布鲁棒模型在推导过程中非常容易出现非线性项。最常见的是双线性项:两个决策变量相乘,例如仿射系数矩阵和不确定量的乘积。处理这类项需要引入辅助变量和额外约束,但这里有一个细节很多人会踩坑——如果你直接在YALMIP里写了双线性项,它会尝试调用非线性求解器,不仅慢,而且经常找不到全局最优解。
一个有效做法是预先识别哪些双线性项必须保留,哪些可以避免。比如当仿射系数矩阵固定时,因为不确定性量是外部参数,所以这个“乘积”其实是线性的——你不需要对这项做线性化,只要把它展开成关于系数矩阵的线性表达式就行。真正的双线性项往往出现在目标函数里,成本和出力相乘的地方,这时可以用分段线性近似(PWL)处理。PWL近似的精度取决于切分的段数,建议对边际成本曲线陡峭的机组多分几段,对平缓的部分少分几段,灵活处理能省不少变量。
5.3 收敛慢、数值病态的排查思路
我在调试阶段碰到过一个特别头疼的问题:Benders分解迭代到十几轮后,上下界差距始终在3%左右徘徊,怎么都压不下去。后来逐条检查割平面,发现是子问题里有一个约束的对偶变量符号写反了,导致生成的割平面方向错误,不仅没帮助收敛,反而把主问题往错误方向带。修复之后,迭代次数直接从四十多次降到了十二次。
另一个常见问题是数值病态。如果约束里的系数跨越多个数量级,比如成本系数在千位级,而出力在百兆瓦级,乘积之后数值范围可能相差几个数量级,求解器很容易报数值不稳定。我的建议是做完无量纲化再求解:把功率统一到标幺值系统,成本统一到相对值,再设定收敛容差为1e-4左右。这个操作在数学上不改变最优解,但会让求解器性能发生质变。有的朋友觉得无量纲化麻烦,实际上用标幺值本来就是电力系统行业的习惯,不存在什么额外成本。
还有一个容易被忽视的点:主问题即使已经达到可行域边界,如果初始解设定得太差,也会让前几轮迭代像无头苍蝇一样乱撞。我的做法是先用一个确定性模型给出初始解,即把风电预测值当作真实值代入,求一个基础UC解,再把这个解作为多阶段模型的初始可行解。这个方法收敛效果稳定,初始化的时间成本也微不足道。
5.4 常见问题速查表
| 问题现象 | 可能原因 | 排查与解决方法 |
|---|---|---|
| 模型求解时间异常长 | 模糊集半径过大,场景数过多 | 减小半径,先做场景约简 |
| Benders分解迭代不收敛 | 割平面方向错误 | 重点检查对偶变量符号和下标映射 |
| 子问题无可行解 | 启停计划不满足爬坡约束 | 在主问题中加入耦合约束的割平面 |
| 目标函数出现负值 | 变量边界设置有误 | 检查成本系数、单位换算 |
| 数值警告或NaN | 矩阵病态,系数差距过大 | 无量纲化,启用稀疏存储,调整容差 |
| 内存不足 | 场景数×时段数爆炸 | 场景聚类,或者改成在线求解方式 |
最后再分享一个实用技巧
如果你只是想要一个基准结果去对比不同方法的优劣,建议先跑一遍纯确定性模型,把它当作整个项目的参考锚点。后续无论你调哪种决策规则、哪套模糊集参数,都拿它做参照系,就能很直观地看出不确定性带来的成本增量花在值不值得的地方。
这个框架后续可以扩展的方向也很多。比如把网络约束加进去,变成安全约束的分布鲁棒单元承诺;或者在日内阶段加入实时校正模型,形成日前-日内两层的全自适应闭环。个人认为,混合决策规则加上分布鲁棒这套思路,在电力市场报价策略、储能容量配置、微电网能量管理这几个方向上都能找到落脚点,不只是单元承诺专用。如果周围有同学在纠结随机规划和鲁棒优化怎么选,不妨把这篇分享给他们,DRO那个甜点区会让你对不确定性建模有全新的感觉。