復现一篇 SCI 一区论文的配电网韧性研究最难的不是看懂模型而是把论文里轻描淡写的约束还原成能跑出结果的 Matlab 代码。这篇专门讲应急移动电源动态调度MPS Dynamic Dispatch的实现接上篇的预配置问题。上篇我们解决了“灾前把移动电源放哪”下篇解决的是“灾中每台移动电源该去哪、该发多少电”。这个问题和静态预配置完全不同预配置是在灾害之前基于预测场景做一次性选址本质是一个随机优化或鲁棒优化问题动态调度则要在灾害发生后面对不断更新的线路故障信息和负荷状态滚动地给每台移动电源规划接入节点和出力既要照顾电网潮流约束又要兼顾车辆移动时间、电池容量甚至交通条件。如果这两块叠在一起模型复杂度立刻上一个台阶。我这次复现用的是 Matlab YALMIP CPLEX 的组合算例是经典的 33 节点辐射状配电网。下面从数学模型开始一步步说到代码实现中间穿插我在调试时踩过的六个坑。如果你也在复现韧性调度类论文或者在做应急电源配置课题这篇大概率能帮你省下至少两周的调代码时间。1. MPS动态调度到底在解决什么从“放哪”到“怎么跑”的衔接1.1 韧性研究的四个阶段与MPS的定位配电网韧性研究一般把极端事件过程拆成四个阶段灾前预防、灾中抵御、灾后响应、恢复重建。传统可靠性关心的是日常故障的平均停电时长而韧性关心的是“台风、冰灾这种低概率高影响事件下系统能不能扛住、能不能快速恢复”。MPS应急移动电源可以是移动储能车也可以是应急发电车在其中的角色很有意思。固定储能装在变电站一块极端事件把线路打断末端重要负荷照样带不上。MPS的价值恰恰在于“机动”故障点在哪电源就能跟到哪通过计划孤岛方式直接接入断线区附近的关键负荷。在抢修班组还没到达现场的时候MPS就是恢复供电最快的手段。1.2 预配置和动态调度的边界划分很多论文把预配置和动态调度放在一个两阶段框架里看。预配置是“灾前决策”输入是台风的路径预报、历史负荷数据输出是每辆 MPS 的初始停靠位置动态调度是“灾中决策”输入是已经暴露的故障线路集合、当前各节点负荷状态、预配置给出的 MPS 初始位置输出是未来若干时段内每辆 MPS 的完整行动计划。要特别提示的是动态调度并不是一次求解就完事。灾害过程中的故障信息是滚动披露的比如台风路径变了、某条线路提前抢通、某条新线路断开了这些信息很难在建模型时全部提前知道所以主流方案是用滚动时域的方式做模型预测控制每个时段求解一次未来时间窗内的优化问题但只执行当前时段的决策下一时段再带着最新状态重新求解。1.3 一个可复现的输入输出框架动手写代码前最好先把输入输出理清楚。我这次复现时给模型定义了这样一套边界输入具体内容网络拓扑节点数、支路数、线路电阻电抗、线路传输容量上限负荷数据各节点 24 小时有功/无功负荷曲线、负荷权重故障场景断线支路集合、断线发生的时段、预计抢通时段MPS 参数电池容量、最大注入功率、移动速度、初始位置预配置结果上篇算出的每辆 MPS 灾前停靠节点输出具体内容时段调度结果每辆 MPS 在每时段接入的节点编号出力计划各时段 MPS 的有功、无功注入量储能状态每辆 MPS 的 SOC 变化曲线负荷恢复结果各节点每时段的切负荷比例、恢复比例韧性指标系统总失负荷能量 ENS、重要负荷恢复时长等这套输入输出框架先定下来后面无论改成 123 节点还是 118 节点系统代码结构都不用重写只需要改数据表。2. 动态调度模型把“MPS该去哪”写成可求解的数学问题2.1 时段离散化与决策变量定义把调度周期分成 T 个时段每个时段时长 Δt。以 24 小时调度为例T24Δt1 小时。节点集合用 N 表示支路集合用 E 表示MPS 的集合用 M 表示。决策变量我分成四组接入状态变量 z(m,i,t) ∈ {0,1}第 m 辆 MPS 在时段 t 是否接入节点 i出力变量 p(m,i,t)、q(m,i,t)第 m 辆 MPS 在时段 t 接入节点 i 时注入的有功、无功功率储能状态变量 soc(m,t)第 m 辆 MPS 在时段 t 末的荷电状态切负荷变量 u(i,t)节点 i 在时段 t 被切除的负荷比例。有些论文把 u 写成 0-1 二进制表示“接上/没接上”但第一版建议先用连续变量 [0,1]跑通后再改二进制因为二进制会让求解时间成倍增加。此外还有潮流变量支路功率 P(ij,t)、Q(ij,t)节点电压幅值平方 V(i,t)。在 DistFlow 模型里V 直接表示电压幅值的平方这样做的好处是很多约束变成线性或二阶锥求解器友好。2.2 潮流约束DistFlow 与二阶锥松弛配电网是辐射状结构经典做法是用 DistFlow 方程描述潮流对于每条支路 (i,j)有功和无功的流动关系写成P(ij,t) - Σ P(jk,t) p_j_inj(t)Q(ij,t) - Σ Q(jk,t) q_j_inj(t)V(j,t) V(i,t) - 2(r·P(ij,t) x·Q(ij,t)) 电压损耗项这里的 r、x 是支路电阻和电抗。电压损耗项包含支路电流平方直接写会导致非凸约束所以标准做法是用二阶锥松弛|| [2·P(ij,t) 2·Q(ij,t) V(i,t)-V(j,t)] ||_2 ≤ V(i,t)V(j,t)这个锥约束把非凸问题松弛成凸问题对辐射状配电网在大多数工况下是精确的误差在可接受范围内。在 YALMIP 里写法非常直接for t 1:T for k 1:nl con [con, cone([2*P(k,t); 2*Q(k,t); V(bus_i(k),t)-V(bus_j(k),t)], ... V(bus_i(k),t)V(bus_j(k),t))]; end end节点注入功率由本地负荷和 MPS 出力共同决定p_j_inj(t) p_j_load(t)·(1 - u(j,t)) - Σ_m p(m,j,t)这个式子的意思是节点注入功率等于“没有被切除的负荷”减去“MPS 注入的功率”。如果 MPS 从节点往外供负荷对网络来说相当于一个负的负荷。2.3 MPS时空耦合约束防止车辆“瞬移”这部分是我复现时最容易出问题的环节也是论文里往往一笔带过的地方。先加最基本的约束每辆 MPS 每个时段最多只能接入一个节点Σ_i z(m,i,t) ≤ 1再加“接入才能出力”的约束p(m,j,t) ≤ P_max · z(m,j,t) q(m,j,t) ≤ Q_max · z(m,j,t)这两条是线性大M约束的紧凑形式本质上就是告诉求解器车不在节点就不能往网络里送功率。然后是移动时间约束。假设车辆从节点 i 移动到节点 j 需要 K(i,j) 个时段那么从 i 出发后K(i,j) 个时段之内是不可能出现在 j 的。翻译成约束就是对任意 m、任意 i≠j、任意 t、任意 k1,…,K(i,j)-1成立z(m,i,t) z(m,j,tk) ≤ 1这个不等式的含义很直白如果车辆在 t 时段停在 i那么在 t1 到 tK-1 这些时段里都不能出现在 j。换个角度看这就是“车辆不会瞬移”的数学表达。它不会限制车辆一直停在原地也不会限制它在满足移动时长之后正常出发只管住“不该出现的时候别出现”。这个时间窗写法的好处是不需要额外引入大量“移动中”状态变量约束数量虽然因为循环嵌套会多一些但对 33 节点系统完全可控。如果你的论文用了更严格的时空网络流模型原理是一样的只是变量更多、求解更慢。2.4 储能SOC约束MPS 作为移动储能车时SOC 的计算不复杂soc(m,t) soc(m,t-1) - [Σ_i p(m,i,t)] · Δt / C_m其中 C_m 是第 m 辆车的电池容量。同时加上容量边界soc_min ≤ soc(m,t) ≤ soc_max如果论文里考虑了 MPS 应急发电车的情况原理也一样——只不过容量上限换成燃料储量SOC 换成油量百分比。多数韧性论文里 MPS 只对外放电不反向充电所以这里不需要充放电的二进制切换变量。2.5 目标函数把韧性变成可最小化的数值韧性提升最终要落成一个可以优化的指标。最常用的目标是最小化加权失负荷能量 ENSmin Σ_t Σ_i w(i) · p_i_load(t) · u(i,t) · Δt其中 w(i) 是节点 i 的负荷权重。医院、通信基站、供水这些重要负荷权重取大比如 100普通居民负荷权重取 1。这样求解器在资源有限时会优先恢复重要负荷。目标函数的符号方向特别容易搞反。有些论文写的是“最大化恢复电量”也就是 Σ w(i)·p_i_load(t)·(1-u(i,t))·Δt和最小化失负荷只差一个常数项总负荷固定时两者等价。但如果总负荷随时间变化就不完全等价。复现时我建议先按论文原文的目标写算出来的数值再和论文表格对比方向和量级都对齐之后再做优化调整。3. Matlab代码实现从模型到可运行程序3.1 数据准备33节点算例与MPS参数我这次用的算例是经典的 33 节点辐射状配电网基准电压 12.66kV总负荷约 3.7MW。数据表格里需要节点编号、节点有功/无功负荷、支路首端/末端节点、支路电阻电抗、传输容量上限。MPS 参数表我设计成这样的形式MPS 编号容量/kWh最大有功/kW最大无功/kVar移动单位时间/(节点/时段)初始位置MPS1150050030018MPS21000350200125MPS3800250150130故障场景我设为时段 8 起线路 5-6、9-10、14-15 断开预计抢通时间分别为第 16、18、14 小时。这里要提醒一句不同论文的故障场景差异很大复现时务必逐条核对论文表格里的断线集合和时段错一个数字结果就对不上。3.2 YALMIP建模变量声明与约束组装YALMIP 的建模思路是声明变量、写约束、设目标、调求解器。核心代码结构如下%% 决策变量声明 z binvar(M, Nb, T, full); % 接入状态 pM sdpvar(M, Nb, T, full); % MPS有功出力 qM sdpvar(M, Nb, T, full); % MPS无功出力 soc sdpvar(M, T, full); % 荷电状态 u sdpvar(Nb, T, full); % 切负荷比例 0~1 P sdpvar(Nl, T, full); % 支路有功 Q sdpvar(Nl, T, full); % 支路无功 V sdpvar(Nb, T, full); % 节点电压幅值平方 %% 约束组装 con []; % 每个时段最多接入一个节点 for t 1:T for m 1:M con [con, sum(z(m,:,t)) 1]; end end % MPS出力限制接入才能出力 for t 1:T for m 1:M for i 1:Nb con [con, pM(m,i,t) P_max(m) * z(m,i,t)]; con [con, qM(m,i,t) Q_max(m) * z(m,i,t)]; end end end % SOC更新 for t 2:T for m 1:M con [con, soc(m,t) soc(m,t-1) - sum(pM(m,:,t))*dt / C(m)]; end end % 二阶锥潮流约束 for t 1:T for k 1:Nl i branch_i(k); j branch_j(k); con [con, cone([2*P(k,t); 2*Q(k,t); V(i,t)-V(j,t)], ... V(i,t)V(j,t))]; end end写完约束后还需要在支路首末端功率、节点功率平衡之间建立关系式。这部分代码比较机械核心是遍历每条支路和每个节点把 DistFlow 的线性等式写进去。我建议先把节点注入功率表达式抽出来做一个中间变量。比如用 inj_p(node,t) 表示节点净注入再用它去写潮流平衡索引出问题时会容易定位。3.3 求解器配置与求解求解配置我建议写成这样ops sdpsettings(solver,cplex, ... verbose,2, ... cplex.mip.tolerances.mipgap,0.01, ... cplex.timelimit,600); sol optimize(con, objective, ops); if sol.problem 0 disp(求解成功); else % 这里一定要看很多隐性错误会通过problem flag暴露 sol.info yalmiperror(sol.problem) end如果机器上没有 CPLEX license可以把 solver 换成 Gurobi或者用开源的 SCIP。但提醒一下带二阶锥约束的混合整数优化规模一大开源求解器速度会明显差一截33 节点 24 时段内还可以接受再大的算例就要认真考虑求解器选型。3.4 滚动时域调度的循环结构动态调度的核心不是单次求解而是多次滚动。我的实现思路是这样% 初始化MPS状态 soc_now soc_init; z_now z_init; results {}; for t0 1:T % 1. 根据t0时刻已知信息生成预测窗口t0 到 t0H-1 % H是预测窗口长度比如6越往后信息越全可以适当延长 % 2. 建立窗口内的优化问题变量、约束、目标 % 窗口内的起始状态直接用当前实际状态 soc_now, z_now % 3. 调用求解器求解窗口问题 % 4. 只提取 t0 时段的决策 % z_execute value(z(:,:,1)); % p_execute value(pM(:,:,1)); % soc_execute value(soc(:,1)); % 5. 更新当前状态 soc_now、z_now把执行结果存入 results end这里最容易踩的坑是“把整个窗口的求解结果都当成执行结果”。记住MPC 原则永远是求解未来 H 个时段、只执行当前 1 个时段。下一时刻重新求解时窗口里的第一个时段就是新的当前时段。索引上稍微偷懒整个结果就会错位。4. 复现中踩过的六个坑以及完整排查链路4.1 坑一MPS瞬移我的第一版模型没写时间窗约束只写了“每时段最多接入一个节点”结果求解出来的轨迹完全离谱MPS1 时段 3 还在节点 8时段 4 直接出现在节点 30中间隔着半个配电网明显不合理。排查方法很简单把每辆车的 z 矩阵抽出来横轴画时段、纵轴画节点编号车辆的运动轨迹一眼就能看出有没有瞬移。我当时画完图就明白了然后补上第 2.3 节的时间窗约束轨迹立刻变得物理可行。这里给个建议无论模型多复杂先把 z 出行轨迹图再做其他分析。4.2 坑二非线性潮流数值病态第一版我直接用有名值写潮流参数电压单位是 V功率单位是 W电阻是 Ω结果二阶锥约束数值跨度从 1e-4 到 1e4CPLEX 求解时反复出现数值警告有些锥约束干脆不闭合。解决办法是把整个系统转到标幺制电压基准取 12.66kV功率基准取 1MVA。改完之后数值量级基本在 0.1~2 之间求解器收敛很快稳定性明显提升。规则很简单所有跟电压、功率、阻抗相关的参数和变量统一换算到标幺值不要混用单位。4.3 坑三大M参数选得太大接入与出力的关联约束里Z 取 0 时要强制 p0。常规写法是加一个足够大的 M于是很多新手会写 M1e6。结果 M 太大LP 松弛后的可行域过松整数变量分支剪枝效率很差求解时间翻了几倍。我最终把 M 取成对应 MPS 最大出力的 1.2 倍比如最大出力 500kWM 就取 600kW。约束依然能有效钳制变量但数值特性好很多。经验法则大M取物理量上限的 1.2~1.5 倍就够不需要无穷大。4.4 坑四求解时间爆炸33 节点、24 时段、3 辆 MPSz 变量就有 3×33×24≈2376 个二进制加上切负荷二进制变量规模不小。第一次全二进制建模CPLEX 跑了半小时还没出最优解。我的压缩手段按顺序进行了三种缩小候选接入节点把断线区域之外、且不承担重要负荷的纯末端节点从候选集里删掉z 变量减少约 40%u 变量从 0-1 改成连续 [0,1]二进制变量数量大幅下降设置 mipgap0.01允许 1% 最优性间隙。改完后单次求解稳定在 2~4 分钟滚动 24 次总时长约 80 分钟。如果你的课题更强调“足够快”还可以在求解前用贪心策略给 CPLEX 喂一个初始可行解效果立竿见影。4.5 坑五结果和论文对不上当我终于跑出曲线后和论文里的数值对比发现 ENS 差了不少。当时第一反应是模型错了后来逐项核对才发现论文用的是“切负荷 0-1 状态”而不是“连续比例”而且故障时段比我设置的晚了一个时段。调整之后数值基本吻合了。我建议的排查顺序是故障线路表 → 故障时段 → MPS 初始位置 → 负荷权重 → 目标函数符号 → 移动时间定义。还有一个隐蔽点论文里“移动一个节点”的时间到底是 1 个时段还是 0.5 个时段代码里 K(i,j) 矩阵的定义会影响整个调度结果必须核对清楚。4.6 坑六滚动时域首时段执行逻辑这个问题非常隐蔽。我第一次实现滚动时域时求解完一个窗口后直接把窗口内所有时段的 z 都记录进了 results导致后面的车辆轨迹跟“上帝视角”一样提前知道了未来故障结果优于论文当时还挺高兴——后来一分析才发现逻辑错了。正确做法是每次只取窗口第一个时段的决策更新当前 MPS 状态和 SOC然后进入下一轮。判断自己的滚动是否正确的简单方法故意在第二个时段增加一个突发故障如果车辆轨迹在故障发生前就有了反应说明你记录错了执行时段。5. 结果分析与扩展方向从复现到真正做研究5.1 结果怎么呈现能读出什么复现完代码之后最需要的是把结果可视化出来方便和论文对比也方便自己理解动态调度的价值。我个人比较推荐的三个图韧性曲线横轴时间、纵轴累计失负荷能量。动态调度、固定接入、无 MPS 三种方案画在同一张图里曲线之间的空隙就是论文的核心贡献。MPS 轨迹图横轴时段、纵轴节点每辆车一条折线配合断线节点标注。能直观看到 MPS 如何从预配置位置出发、移动到关键负荷区域、最后撤离。SOC 曲线每辆车一个子图看储能是否在关键时段保持足够余量。如果某辆车总是处于满电或空电的极端状态说明容量配置或者调度策略有问题。绘图代码本身不复杂figure; plot(1:T, ens_dynamic, LineWidth, 2); hold on; plot(1:T, ens_fixed, --, LineWidth, 2); plot(1:T, ens_nomeps, :, LineWidth, 2); legend(动态调度,固定接入,无MPS); xlabel(时段); ylabel(累计失负荷能量(kWh));5.2 从复现代码出发的三个扩展方向代码跑通之后如果想把这套逻辑延伸到自己的研究里我建议从三个方向切入。第一多类型移动资源协同。把 MPS 从“只能放电”扩展为“可充可放的移动储能”甚至加入应急发电车、移动制氢车等对象车辆之间可以交接供电任务。模型只需要增加充放电二进制变量和交接时序约束整体框架不用动。第二MPS 与配电网其他灵活资源联合调度。配电网里往往还有分布式光伏、固定储能、软开关 SOP。把它们纳入同一个调度模型MPS 可以只负责“补缺口”其他资源负责削峰填谷恢复效果和经济性都可能更好但模型量级会上去计算复杂度也会明显增加。第三交通网络与配电网耦合。现实中的 MPS 移动时间不是欧氏距离决定的而是受路网可达性、道路拥堵甚至洪水淹没道路影响。可以预先用路网最短路径算法生成节点间移动时间矩阵 K(i,j)替换掉原来的直线距离矩阵本质上不改变模型结构但结果会贴近实际。5.3 一点经验收尾这次复现最花的部分不是数学模型推导不是求解器调参而是把约束条件一条条翻译成 YALMIP 时保持索引一致。变量一多任何一个维度错位都会导致求解器报不可行而报错信息又往往不直接指向问题根源。我后来养成一个习惯写完约束先打印每个 sdpvar 的 dimension再用一个极小的 3 节点 3 时段算例跑通最后才上 33 节点全规模。多时段模型出了问题90% 出在时间索引和 MPS 移动约束上这个判断在之后好几个项目里都成立。如果你刚跑通预配置的静态模型我建议先别急着上 MPC先把单时段动态调度模型跑稳、确认目标值和论文对得上再往多时段滚动扩展。方法再先进还是要从能复现的第一行代码开始。