做电力系统优化的朋友应该都碰过单元承诺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.2429.8412.6最坏分布下成本万元529.3451.2435.7求解时间秒18.645.352.1切负荷期望MW85.423.615.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那个甜点区会让你对不确定性建模有全新的感觉。
