应急移动电源动态调度：配电网韧性提升的Matlab实现与优化

配电网韧性和应急移动电源调度这个方向这几年在电力系统顶刊里出现频率很高尤其极端天气频发之后台风、冰灾一过配网大面积停电的场景大家都不陌生。传统抢修是人力巡检加定点修复恢复速度慢停电期间医院、通信基站、应急指挥中心这类重要用户很难保障供电。这时候应急移动电源MPSMobile Power Source的价值就体现出来了——本质上它就是一台可以“跑起来”的储能电源车在故障发生后移动到关键节点附近临时给重要负荷供电。但问题来了MPS数量有限、移动速度有限、配电网拓扑又会随故障变化怎么在正确的时间把电源派到正确的位置就是动态调度要解决的核心问题。我最近在复现一篇SCI一区论文里的MPS动态调度部分这篇文章分了两个阶段第一是灾前的预配置第二是灾后的动态调度。上篇已经整理过预配置阶段这篇重点讲下篇——MPS动态调度的Matlab实现包括数学建模思路、代码架构、求解器选型以及我复现过程中踩过的一些坑。适合正在做配电网韧性研究、需要复现调度类论文的研究生也适合想用Matlab做移动应急资源优化调度的工程师参考。1. 问题拆解MPS动态调度到底在调度什么1.1 配电网韧性提升的三个层次韧性Resilience这个概念通俗理解就是配电网对极端事件的“扛打击”和“恢复”能力它跟可靠性Reliability不是一回事。可靠性处理的是发生概率较高的常规故障比如设备随机失效、线路雷击跳闸故障概率可以用历史统计数据描述。韧性处理的则是概率低、但后果严重的大规模极端事件比如台风导致多条馈线同时断线、大面积负荷失电这种场景用期望值建模意义不大更需要关注最坏情况下的系统应对能力。学术上通常把韧性研究拆成三个层次第一是系统级韧性评估通过构造韧性指标负荷损失率、恢复时间、韧性曲线下面积来量化系统的应对能力第二是网络加固与资源预配置在灾前对线路设备进行升级改造或者在关键位置预先部署应急资源第三是运行调度优化在灾害发生过程中和发生后通过调度手段线路重构、分布式电源出力、MPS移动最大化恢复供电。MPS动态调度属于第三层但它跟前两层的耦合很强。你在论文里会看到动态调度模型里通常要嵌入配电网潮流约束这就涉及系统建模能力而MPS的初始位置又来自第一阶段的预配置结果这就把第二层和第三层串了起来。复现的时候如果只盯着第三层的调度模型忽略它和前一层的关系代码里很容易出现初始条件对不上、结果物理上不成立的问题。1.2 预配置与动态调度两阶段问题的逻辑关系这种两阶段结构在调度类文献里非常常见本质上是一个“先决策、后观察、再决策”的框架。第一阶段是灾前预配置极端事件还没发生电网只能根据气象预报估计灾害影响范围不确定性很强。此时要决定MPS的初始部署位置和数量分配目标是让MPS处在一个“进可攻退可守”的位置——既能覆盖大概率受灾区域又保留一定的机动能力。第二阶段就是本文的核心灾后动态调度。此时故障信息逐步确认哪条线路断了、哪些负荷失电调度员需要动态调整MPS的移动路径和接入节点让有限的MPS在合适的时间到达合适的位置最大化恢复供电。这里要特别强调一点预配置和动态调度并不是两个独立的优化问题而是时序耦合的。预配置阶段输出的MPS初始位置必须作为动态调度阶段的初始状态约束。我在复现时发现如果代码里漏掉这层约束优化模型会“自由地”让MPS出现在任意节点求解结果虽然好看但物理上完全不成立。所以写代码之前建议先在纸上把两阶段的数据流图画清楚预配置输出哪些变量哪些变量作为动态调度的输入两个文件之间如何传递参数一清二楚之后再动手写代码。2. 数学建模把调度问题写成可求解的优化模型2.1 目标函数怎么设计动态调度问题的目标函数论文里通常会在几个候选里选一个复现时别急着全部实现先从最基础的目标入手。最常用的目标函数是恢复供电量最大化。把调度周期离散成时段后每台MPS对某个负荷节点供电该时段内就有一定的恢复电量把全部MPS、全部时段、全部恢复负荷的电量累加就是目标值。公式形式上很简单但它的权重设计值得留意重要负荷医院、应急指挥中心和不重要负荷的权重往往不同论文里会给每个负荷一个权重系数复现时要用原文的权重否则恢复策略的重点会偏离。第二个常见目标是韧性指标优化。有些论文直接用韧性曲线的积分面积作为目标横轴是时间纵轴是系统可供电负荷比例曲线下面积越大说明系统恢复越快、断电影响越小。这个目标跟第一个目标在很多时候是等价的但前者的表达更贴近“韧性提升”这个研究主题。第三类目标加上了调度成本在恢复收益里扣掉MPS的移动成本、运行维护成本变成一个净收益最大化问题。这类目标建模更复杂因为要刻画移动成本的函数形式通常是线性或阶梯函数但物理意义更强。我复现的时候建议先做“恢复电量最大化”最简单直观跟算例结果也好对照。代码跑通之后再往目标函数里加成本项一步步扩展不要一开始就追求跟论文里的复杂模型完全一致。2.2 关键约束条件逐条拆解约束条件才是MPS动态调度的难点。我列几个核心约束每一条都是代码里容易出错的地方。配电网潮流约束。MPS接入配电节点后系统潮流分布会改变严格的做法是使用DistFlow支路潮流模型对每个节点建立有功、无功平衡方程。这里有两个关键点一是辐射状配电网的DistFlow方程可以用Big-M法线性化把线路开关状态、MPS接入状态通过二进制变量耦合进去二是线路开断之后对应的支路潮流量必须强制为0否则优化模型会“利用”断开的线路虚拟送电解出来是一个假的最优值。我在复现时就吃过这个亏检查了很久才发现是断线支路的潮流上限没有归零。MPS移动约束。每台MPS在同一时刻只能位于一个节点从节点i移动到节点j需要时间t_ij这个约束跟经典的车辆路径问题很相似。移动时间矩阵通常预先算好假设MPS沿道路网以固定速度行驶用Floyd或Dijkstra算法求所有节点对之间的最短路径时间。注意移动时间矩阵需要预先算好因为把它嵌入MILP里会增加大量非线性约束求解器的负担会迅速变大。功率约束。MPS有额定容量、最大输出功率、充放电效率等参数限制这些约束必须精确表达。尤其是“MPS接入节点后提供的功率不超过额定容量同时不超过该节点允许的最大注入功率”这一条很多初次建模的人会漏掉导致MPS在单个节点注入的功率超过线路容量结果不满足物理实际。状态转移约束。MPS的位置变量在不同时段之间要满足转移逻辑上一时段在节点A下一时段要么还在A要么移动到A的邻近可达节点。这个约束要和移动时间矩阵配合起来写否则模型会出现MPS“瞬移”的不合理结果。2.3 时间离散化与场景处理动态调度是一个时序决策问题必须把连续时间轴离散成时段。常用做法是把灾害响应周期切成等长时段比如以1小时为步长共24个时段每个时段内认为系统状态保持不变。时段长度是权衡结果太短决策变量会爆炸每个MPS每个时段都要定义位置变量、接入变量、功率变量问题规模指数增长太长MPS的移动过程和被忽略的负荷变化会导致调度策略失真。场景处理也是动态调度建模里绕不开的点。论文里通常会构造多个故障场景比如不同位置的线路断线组合对应不同的失电范围。复现时可以先从单一场景开始跑通之后再扩展成多场景鲁棒优化或场景树。多场景会增加约束数量代码结构不变但求解难度明显上升这个要心里有数。3. Matlab实现代码架构与求解器选型3.1 整体代码架构怎么组织复现SCI论文里的调度模型代码绝不能是一坨全写在一个脚本里。我建议按模块化组织每个模块只干一件事main.m主程序负责数据读取、模型构建、求解、结果输出data/存放系统参数、负荷数据、故障场景数据model/构建优化模型的函数目标函数和约束在这个目录里定义solver/封装求解器调用的函数utils/工具函数比如最短路径计算、韧性指标计算output/结果输出与可视化脚本这种组织方式对调试特别重要。学术复现的代码很少能一次写对模块化之后你可以单独验证路径计算模块单独验证潮流约束模块出了问题不用在几千行的代码里反复翻找。我自己的习惯是每写完一个模块就先用一个最小用例测一下确认输出合理再继续下一个。等全部模块拼起来的时候绝大部分低级错误已经提前排掉了。3.2 MPS移动路径与最短路径矩阵计算MPS的移动时间矩阵是整个调度模型的核心输入之一。很多论文为简化直接假设MPS在节点之间沿直线移动移动时间等于欧氏距离除以平均速度。如果原文给了移动时间矩阵直接用原文数据如果没给就用节点坐标自己算。深度优先遍历或Floyd算法都可以实现最短路径时间计算。我习惯用Floyd代码短、思路清晰33节点系统上运行耗时基本可以忽略。算出所有节点之间的最短路径时间后存成一个N×N的矩阵N是节点总数。这个矩阵在MILP建模时直接以参数形式传入不需要参与优化变量的构建。用YALMIP建模时MPS移动约束可以写成如下形式% 假设 n_mps 台MPSn_t 个时段n_b 个节点 % 决策变量 x_mps(m, i, t) 表示第m台MPS在时段t是否位于节点i二元变量 x_mps binvar(n_mps, n_b, n_t, full); % 每一台MPS在每个时段只能位于一个节点 for m 1:n_mps for t 1:n_t F [F, sum(x_mps(m, :, t)) 1]; end end % 状态转移约束位置变化必须匹配移动时间矩阵 travel_time(i,j) % 实际实现时需要根据移动时间与时段长度的大小关系构造可达性矩阵这里最容易出错的是移动时间和时段长度的匹配。假设设定的时段长度是2小时但MPS从节点5到节点8需要3小时那么在这2小时内的状态切换逻辑就不能简单用相邻时段的位置变量来约束。一种处理方式是引入“移动中”状态变量另一种是把时段长度缩小到小于最短移动时间。后者实现简单但会增加变量数量。我建议先按论文实际设定来如果论文的时段长度和移动时间存在冲突那就优先修改移动速度参数让它落在合理的范围内。3.3 求解器选型intlinprog还是YALMIP加外部求解器MPS动态调度本质上是混合整数线性规划问题。Matlab环境里做MILP主要有两条路。第一条路是直接用Matlab自带的intlinprog优点是不需要额外安装求解器开箱即用也没有许可证问题。缺点非常明显对大中型MILP问题求解效率偏低整数变量多了之后经常要跑几十分钟甚至几小时。我实测过IEEE 33节点、24时段、3台MPS这个规模intlinprog的求解时间通常超过30分钟而且不一定能证明全局最优。第二条路是用YALMIP建模底层调用Gurobi或CPLEX求解。这是学术复现的主流配置求解速度比intlinprog快一个数量级。同样规模的问题Gurobi常常几分钟就能收敛到很紧的MIP Gap以内。代价是许可证不过Gurobi和CPLEX都有学术license学校邮箱申请很方便。用YALMIP建模的好处除了速度还有建模的直观性。矩阵形式的约束intlinprog需要A·x≤b在约束条目多时非常容易拼错而YALMIP里直接写约束表达式再用[]拼接可读性和可维护性高很多。建模完成后调用optimize求解再通过value()取出变量值进行分析。% YALMIP建模示例示意 x binvar(n_mps, n_b, n_t, full); % MPS位置 p sdpvar(n_b, n_t, full); % MPS注入功率 F []; for t 1:n_t for i 1:n_b F [F, 0 p(i, t) p_max]; % 功率上限 F [F, p(i, t) sum(x(:, i, t)) * p_max]; % 位置耦合 end end optimize(F, -objective, options); x_opt value(x); p_opt value(p);options里可以设置求解时间上限、MIP Gap等参数。Gurobi的TimeLimit和MIPGap控制在复现阶段非常实用先用一个宽松的时间上限把模型跑出可行解再收紧参数提高精度。4. 复现实战数据准备、代码调试与结果可视化4.1 数据准备与算例选择复现调度的第一步是确定算例系统。绝大多数配电网韧性论文用的是IEEE 33节点系统因为节点少、参数公开、拓扑结构清晰适合做算法验证。也有一些论文用IEEE 123节点或实际馈线系统没有原文系统参数时先用IEEE 33节点跑通最稳妥。数据准备工作量不小我一般把整块数据拆成几类处理线路参数电阻、电抗、容量、负荷参数每个节点的有功、无功需求以及重要负荷权重、故障场景数据线路断线位置、失电负荷集合、MPS参数数量、容量、最大功率、移动速度。把这些数据统一放在data/目录下用脚本读取这样换算例、改参数都不用动代码主体。这里有个实用建议把负荷权重单独做成一个向量与目标函数里的权重系数直接对应。很多论文里的目标函数有一个权重向量W我在第一次复现时把权重和节点混在了一起导致约束和目标里的参数搞反结果解出来的恢复策略把所有电源都集中到权重最高的节点上其他节点全不管。单独管理权重向量之后这种问题就很好排查。4.2 代码调试过程中的典型难点调试阶段会遇到几个反复出现的坑我先挑典型的说一下。索引和编号方向容易错位。Matlab数组索引从1开始论文里的节点编号通常也从1开始但支路数据的起始节点和终止节点方向容易在构造邻接矩阵时搞混。我建议统一用“起点-终点”边列表记录线路避免邻接矩阵行列含义混淆。同时要把节点编号和线路编号分开管理掐指一算就能对应上。Big-M参数取值很关键。线性化二进制变量与连续变量的耦合时M值太大会导致数值稳定问题太小会截断可行域。经验值是取该支路或该节点功率上限的2到3倍。比如线路容量是5MW取M10既不会截断可行域数值上也不会出现极端的尺度差异。还有我前面反复提到的时间步长与MPS移动时间匹配问题这是最隐蔽的坑之一。模型看似能跑出最优解但仔细检查MPS移动轨迹你会发现它在一个时段内从节点A跳到了距离很远的节点B这就是移动约束没有生效的典型表现。排查方法是在结果分析阶段逐时段打印MPS位置对照移动时间矩阵检查是否合理。4.3 结果可视化与韧性指标计算复现完成之后用贴近论文风格的图形展示结果既是给自己检查也是后续组会汇报的素材。至少应该输出四类图。第一张是配电网拓扑图在图上标注故障线路、失电负荷区域和MPS的最终接入位置。第二张是各时段恢复负荷比例的阶梯图或柱状图直观展示恢复进程。第三张是MPS移动轨迹图把每台MPS在不同时段的位置在拓扑图上连成路径用来验证移动约束的合理性。第四张是韧性曲线横轴时间纵轴系统可供电负荷比例或恢复电量计算曲线下面积作为韧性指标。韧性曲线下面积用Matlab的trapz函数直接算这是最简单的梯形积分几行代码就够% 恢复比例向量 load_recovery_ratio每个时段一个值 % 时段长度向量 time_step小时 resilience_index trapz(time_steps, load_recovery_ratio);这个指标可以直接和基准场景无MSP调度对比量化MPS动态调度带来的韧性提升。如果你的复现结果里韧性指标比论文低了不少优先检查两件事一是故障场景是否一致二是MPS的移动约束是否被模型简化掉了。这两个原因导致的指标偏差最常见。5. 常见问题与排查技巧实录5.1 求解时间过长怎么办MILP求解时间爆炸是最常见的问题。变量数量约等于时段数乘以节点数乘以MPS数量再乘以状态变量类型数规模稍大就容易卡住。应对思路按以下顺序排查。第一步压缩问题规模。考虑减少时段数量比如从24时段改成12时段或者把负荷节点适当聚合。第二步给求解器设置时间上限和MIP Gap。Gurobi的TimeLimit参数intlinprog用MaxTime选项先设一个600秒的上限如果MIP Gap能在5%以内结果基本可接受如果超过5%再收紧。第三步用启发式解作为热启动。先在一个简化模型上跑一个粗略解比如不考虑MPS移动时间作为MILP初始可行解传入可以有效减少分支定界树探索量。具体到我个人经验Gurobi的MIPFocus参数值得提一下。默认值偏平衡如果模型明显受限于可行解的搜索速度可以尝试MIPFocus1偏向快速找可行解如果上界已经不错但下界提升慢用MIPFocus2。这个小参数在复现时经常发挥奇效。5.2 模型报不可行怎么排查求解器返回infeasible是另一个高频问题。我的排查套路是逐层放松约束做二分定位。先去掉MPS移动约束看模型是否可行。如果去掉后可行说明问题出在移动约束和功率约束的耦合上。再检查负荷恢复变量与节点状态的关系失电节点的负荷恢复变量必须强制为0。有些模型忘了加这条“激活”逻辑或者加的方式不对数学上会产生矛盾约束导致全局不可行。YALMIP有个很实用的命令是查看约束冗余和可行性虽然不可行时的报错信息有时比较模糊但可以用check(F)逐条检查约束满足情况它会告诉你哪些约束有较大残差。结合缩小规模到7节点系统排查效果很好。7个节点、3个时段的小系统可以很快定位问题在哪条约束然后针对性地修正大模型里的对应部分。5.3 结果不合理先查约束残差如果模型有解但结果不合理比如MPS“瞬移”、恢复负荷超过该节点最大负荷、线路潮流超过容量先在Matlab里把变量值代入约束表达式一列一列检查残差。具体操作是求解完成后用value()取回所有变量再重新计算每条约束表达式的取值跟上下界比对。比如检查恢复负荷约束就把恢复功率变量求和看看是否超过该节点最大有功需求。检查潮流约束就把支路潮流量跟线路容量上限比对。这种方法虽然笨但对定位约束写错、索引错位、Big-M取值不当这些问题非常有效。另外一个值得注意的细节是结果不合理往往并不是某一处写错而是数值问题引起的微小违规。比如节点电压幅值约束的松弛量不够或Big-M给的边界刚好卡在可行域边上。此时看残差大小如果只是0.001量级的偏差基本不影响宏观决策不必追求数值上的绝对完美。最后分享一点个人体会我复现这类调度论文的最大感受是可复现性不高的原因往往不在优化算法本身而在模型与数据的对应关系。论文里的公式链条很长每个符号对应哪一行代码、每条约束对应哪一条表达式必须先理清逻辑否则调试时根本无从下手。MPS动态调度的核心卖点是“移动电源在正确的时间出现在正确的位置”所以移动过程建模的质量直接决定结论的可信度。复现时不要为了省事简化移动约束一定要把移动时间矩阵和状态转移约束写扎实。最后分享一个小技巧动手写代码之前先手工画一张时间-位置状态转移表把MPS在典型故障场景下的期望动作推演一遍。比如故障发生后预期第一台MPS应该在哪个时段到达哪个节点第二台需要跨几个时段移动。然后用这张预期表去对照程序输出一旦程序结果和预期不一致大概率能快速定位问题出在哪个模块。这个小习惯帮我省了大量调试时间也推荐给正在复现这类代码的同学。