台风刚过境的那个夜晚配电网上游主网失电多条馈线同时跳闸城市边缘的几块区域瞬间成了孤岛。这时候成本最低、见效最快的补救手段就是应急移动电源MPS——说白了就是一辆能跑能供电的储能车。问题来了车平时放在哪、灾害发生后开到哪、什么时候接进配网、以多大功率放电这些动作叠在一起就是一个典型的预配置动态调度两阶段优化问题。上一篇文章我们把预配置阶段讲透了也就是“提前把车部署到哪些候选点”这篇专门讲下篇MPS动态调度怎么做以及对应Matlab代码怎么落地。我会从数学模型、代码架构、测试场景到复现中踩到的坑逐层拆开来讲。这篇文章适合两类人一类是在做配电网韧性、应急电源优化调度方向的研究生想把论文里的算法快速复现成可跑的实验代码另一类是电力应急装备规划的业务人员想搞清楚移动储能车在灾害场景里到底怎么分配才算最优。我的实现思路用一个词概括整数线性规划。把时间、位置、能量三个维度全离散化最后形成一个混合整数线性规划MILP模型用Matlab优化建模工具箱声明变量、写好约束交给求解器去算。代码不算复杂但细节不少尤其是移动时间建模和潮流约束这两块最容易出错。1. 上篇回顾预配置解决“电源放哪”动态调度解决“怎么走、怎么用”1.1 两阶段模型在韧性评估里各管哪一段配电网韧性这个概念通常被拆成四个阶段预防、抵御、响应、恢复。移动电源车在当中能干的活主要集中在响应和恢复阶段。灾害前我们对灾后的故障场景只能做概率猜测所以只能先把有限的MPS部署到几个可能用得上的候选节点这是预配置阶段对应论文里的第一阶段决策。一旦灾害真的发生故障线路集合、负荷损失情况基本明朗了这时候车往哪里开、在哪个节点接入、每个小时放多少电就是动态调度阶段也就是第二阶段决策。你可以把这两个阶段理解成消防系统的逻辑消防站建在哪是第一步规划火灾发生后消防车走什么路线、在哪灭火是第二步调度。两者强耦合——预配置阶段给出的初始位置直接决定了动态调度阶段每台车从哪个节点出发。所以在复现代码时第一阶段输出要作为第二阶段的固定参数输入不能各算各的。1.2 动态调度模型的输入输出我复现的动态调度模型输入包括五块网络拓扑数据节点编号、支路阻抗、线路容量、联络开关位置负荷数据各节点24小时有功和无功曲线MPS参数数量、额定容量、最大充放电功率、充放电效率、初始SOC、移动速度预配置结果每台车在灾害开始时的初始位置故障场景哪些线路在哪些时段断开。模型输出的是一个完整的调度方案包括每一台MPS在每一个时段的接入节点、移动状态、充放电功率、剩余电量以及配电网里每个节点的负荷恢复比例、电压和潮流分布。这里有一个很重要的前提假设故障场景被视为已知或可以预测。实际场景里线路故障往往是逐步发生、不断更新的所以更工程化的做法是滚动优化——每1到2个小时重新求解一次用最新信息修正方案。这个机制我放在后面详细说你先记住模型本身是面向单个确定性场景的但代码架构要预留滚动重算的接口。2. 动态调度的数学模型目标、变量、约束怎么一步步搭起来2.1 目标函数先想清楚你要优化什么先把目标函数写出来因为我发现很多初学者第一步就卡在这到底是最小化失负荷量还是最大化恢复量其实两者等价但权重不同结果差异很大。我采用的目标是最小化加权失负荷总量表达式如下[ \min \quad \sum_{t1}^{T}\sum_{i1}^{N} \lambda_i (1-\alpha_{i,t}) P_{i,t}^{D}\Delta t\gamma \sum_{k1}^{K}\sum_{t1}^{T}\sum_{m} y_{k,m,t}^{move} ]其中 (\alpha_{i,t}) 是节点 (i) 在时段 (t) 的负荷恢复比例取值0到1(P_{i,t}^{D}) 是该节点原始有功负荷(\lambda_i) 是节点重要度权重一级负荷权重取10二级取3三级取1第二项是MPS移动次数的惩罚(\gamma) 设得很小比如0.001目的是让求解器在多组可行解里优先选动作更少的方案。为什么要设权重因为灾害恢复不是“能送多少电就送多少电”而是要让有限电力优先供应医院、通信基站、供水设施这类关键负荷。如果不加权求解器会把电力平均分给所有节点这在工程上是不合理的。2.2 决策变量时间、空间、能量三个维度跑动态调度模型本质上是把MPS的运行状态拆成离散变量组合。我的模型里主要有三类决策变量第一类是0-1整数变量包括(x_{k,i,t})MPS (k) 在时段 (t) 是否接入节点 (i)(y^{move}_{k,m,t})MPS (k) 在时段 (t) 是否开始沿路径 (m) 移动(z_{ij,t})时段 (t) 线路 (ij) 是否闭合且为投运状态。第二类是连续变量包括(\alpha_{i,t})节点 (i) 的负荷恢复比例(p^{dis}{k,t})、(p^{ch}{k,t})MPS放电、充电功率(E_{k,t})MPS剩余电量潮流变量 (P_{ij,t})、(Q_{ij,t})、(U_{i,t})、(l_{ij,t})支路有功、无功、节点电压幅值平方、支路电流幅值平方。第三类是辅助连续变量比如大M线性化用的中间变量这在后面的约束里会用到。变量规模很好算假设33个节点、2台MPS、24个时段那么 (x) 变量的数量就是 (33 \times 2 \times 24 1584) 个二进制变量加上移动变量、线路投运变量总共两三千个整数变量。这个规模对现代求解器来说完全不是问题但如果你直接把节点数放大到100以上变量数会迅速膨胀后面我会讲怎么削减。2.3 约束条件MPS自身约束和配电网约束MPS自身的约束第一个是能量递推[ E_{k,t1} E_{k,t} \eta^{ch} p^{ch}{k,t}\Delta t - \frac{p^{dis}{k,t}}{\eta^{dis}}\Delta t ][ 0 \le E_{k,t} \le E_k^{max} ]第二是放电功率必须要和接入状态绑定只有接入某个节点才能对外放电。[ 0 \le p^{dis}{k,t} \le P_k^{max} \sum{i \in \mathcal{N}} x_{k,i,t} ]第三是接入和移动互斥同一时刻一台MPS要么停在某个节点放电要么正在路上移动不能两边同时发生。[ \sum_{i \in \mathcal{N}} x_{k,i,t} \sum_{m \in \mathcal{M}} y^{move}_{k,m,t} \le 1,\quad \forall k,t ]移动时间约束是这里最容易写错的地方。如果两点之间的通行时间超过一个时段你不能只用一个0-1变量表示“正在移动”因为那样等于允许MPS瞬移。我的处理方式是引入“在途剩余时间”变量 (r_{k,t})当某台车开始移动时(r_{k,t}) 被置为所需的出行时段数之后每个时段减1在 (r_{k,t}0) 期间该车不能接入任何节点。具体代码实现我会在下一节展示。配电网侧的约束核心是潮流方程。我用的是一阶线性化的DistFlow模型这也是很多韧性研究论文的常规选择因为它能保持整个模型线性求解稳定又高效。对每条支路 (ij)三个方程[ P_{ij,t}-\sum_{m \in \delta(j)} P_{jm,t} - r_{ij}l_{ij,t} P_{j,t}^{net} ][ Q_{ij,t}-\sum_{m \in \delta(j)} Q_{jm,t} - x_{ij}l_{ij,t} Q_{j,t}^{net} ][ U_{j,t} U_{i,t} - 2(r_{ij}P_{ij,t}x_{ij}Q_{ij,t}) (r_{ij}^2x_{ij}^2)l_{ij,t} ]其中 (P_{j,t}^{net}) 是节点 (j) 的净注入有功等于该节点的负荷消耗减去MPS注入。如果节点 (j) 被MPS接入这个净注入就会变小甚至变成负的相当于电源向网络送电。还有辐射状拓扑约束。配电网正常运行时不允许合环所以每个非根节点在同一时刻只能有一个父节点用线路投运变量来表达[ \sum_{i \in \pi(j)} z_{ij,t} \le 1,\quad \forall j \ne r ]故障线路对应的 (z_{ij,t}) 直接固定为0。MPS接入节点时可以把该节点临时当成一个孤岛电压源也就是说它可以在不依赖上游线路的情况下维持局部供电。这个逻辑必须写进约束否则求解器会让没有通电的节点通过MPS恢复供电导致潮流方程不自洽。2.4 从离线全局到滚动优化上面建立的是全时段联合优化的模型——假设求解器知道未来24小时的所有故障信息和负荷曲线一次性算出全局最优。这在复现论文时是对的论文给的结果大多就是这个“上帝视角”解。但真实灾害中故障可能随着灾情演化不断变化负荷也在实时波动你不可能提前知道完整曲线。所以工程落地时我更推荐改成滚动优化每个时段开始时用当前状态作为初始条件只对未来 (T_w) 个时段求解一次执行前 (T_c) 个时段的方案然后更新时间窗口再求解。代码只要把场景数据和初始状态参数化就能在同一套建模函数里反复调用增加的代码量很小。这也是为什么很多论文会在实验里对比“离线全局优化”和“滚动优化”前者给出性能上限后者给出工程可行性。你复现的时候两组实验都做出来结果表会更有说服力。3. Matlab代码架构拆解从数据结构到求解循环3.1 模块划分与数据组织先把代码目录搭好我习惯分成四个目录data放网络参数和场景配置model放建模函数solver放求解入口utils放后处理和画图工具。主脚本只负责调度流程。核心数据结构我用三个structnetwork、mps、scenario。network.nb 33; % 节点数 network.lines [...]; % 每条线路 [首端 末端 r(pu) x(pu)] network.load load_24h; % Nx24 有功负荷(pu)单位MW network.qload ...; % Nx24 无功负荷(MVar) network.sbase 1.0; % 功率基准 1MVA network.vmin 0.95^2; network.vmax 1.05^2; % 电压幅值平方上下限 network.linecap ...; % 支路容量上限(pu) mps.K 2; mps.cap [0.5; 0.5]; % 单位MWh mps.pmax [0.1; 0.1]; % 单位MW mps.soc0 [0.9; 0.9]; % 初始SOC mps.initNode [12; 26]; % 来自预配置阶段结果 mps.eta [0.95; 0.95]; mps.speed 2; % 每时段可经过的线路条数 scenario.faultLines [5 6; 8 9; 12 13; 28 29]; scenario.faultStart 9; scenario.faultEnd 18; scenario.T 24; scenario.dt 1;主脚本大概长这样clear; close all; clc; addpath(./data,./model,./solver,./utils); setup_case33; model build_dispatch_model(network, mps, scenario); ops sdpsettings(solver,*,verbose,0,debug,1); sol optimize(model.cons, model.obj, ops); if sol.problem ~ 0 warning(求解失败: %s, sol.info); return; end x_val value(model.x); alpha_val value(model.alpha); E_val value(model.E); post_process(x_val, alpha_val, E_val, network, mps, scenario);那个solver,*的意思是让Matlab优化建模工具箱自己检测本机可用的MILP求解器只要你装了任何一个商用级求解器都能跑起来。3.2 关键代码段变量的声明与约束组装建模函数build_dispatch_model是我认为整段代码中最值得反复检查的部分。别急着一次性把所有约束写完先声明变量再把约束分组每写完一组就加一条注释。% ---------- 决策变量 ---------- x binvar(N, mps.K, T, full); % 接入状态 y binvar(mps.K, T, 3, full); % 移动方向状态预留3条可能路径 z binvar(N, N, T, full); % 线路投运状态 alpha sdpvar(N, T, full); % 负荷恢复比例 pdis sdpvar(mps.K, T, full); % 放电功率 pch sdpvar(mps.K, T, full); % 充电功率 E sdpvar(mps.K, T1, full); % 剩余电量 r_ontrip sdpvar(mps.K, T1, full); % 在途剩余时段 Pij sdpvar(N, N, T, full); % 支路有功 Qij sdpvar(N, N, T, full); U sdpvar(N, T, full); % 电压幅值平方 l sdpvar(N, N, T, full);特意说明两点。第一alpha是我在全时段统一引入的负荷恢复比例不拍脑袋设0-1整数而是连续变量因为一个节点可以部分恢复供电第二z虽然是在全网节点对之间定义但实际只有存在支路的节点对才会加约束其他地方会自动被限制为0。接着判断是否该写约束了。先把变量组合进约束里cons []; % ---------- MPS能量约束 ---------- for k 1:mps.K for t 1:T cons [cons, E(k,t1) E(k,t) mps.eta(k)*pch(k,t)*scenario.dt ... - pdis(k,t)*scenario.dt/mps.eta(k)]; cons [cons, 0 E(k,t1) mps.cap(k)]; cons [cons, 0 pdis(k,t) mps.pmax(k)*sum(x(:,k,t))]; cons [cons, 0 pch(k,t) 0]; % 灾害孤岛阶段通常不充电先留接口 end end灾害恢复阶段MPS接入孤岛时主要任务是放电但如果某个孤岛上有新能源出力MPS也可以吸收多余电量充当储能。代码里我先把这个通道关掉等你有新能源模型时再把上限改成mps.pmax*sum(x(:,k,t))。移动约束这块我写一个简化但正确的方案假设每台车同一时段只能选择一个方向移动移动所需的时段数通过travelTime矩阵查表得到。% travelTime(i,j): 节点i到j需要经过的线路数按速度折算 cons [cons, sum(sum(x,1),2) sum(y(k,:,:),all) 1]; % 在途状态开始移动时置为对应旅行时间 % 这里简化每时段最多移动1条边多段行程用连续多个时段的y实现如果你不想太过简化可以引入旅行时间表。我在完整代码里用的是“节点到节点直达表”的方式预计算travelTime矩阵然后当y触发时用约束把r_ontrip设为travelTime再随时间递减。这样虽然要多写几行约束但模型不会出现瞬移。3.3 潮流约束与辐射状拓扑的组装潮流约束是建模里最绕的部分但套路固定。先初始化全部线路断开cons [cons, z network.topology]; % 无线路的地方永远不能闭合 for idx 1:size(scenario.faultLines,1) i scenario.faultLines(idx,1); j scenario.faultLines(idx,2); cons [cons, z(i,j,:) 0]; % 故障时段强制断开 end这里我没把“故障时段”和“全天断开”分开实际代码里你可以针对faultStart:faultEnd区间单独约束。然后把DistFlow方程逐条写进循环for t 1:T for i 1:N % 从父节点流到i的功率 节点净负荷 流向子节点的功率 % 这里用支路功率变量组装注意方向性 cons [cons, Pnet(i,t) network.load(i,t)*(1-alpha(i,t)) ... - sum(pdis(:,t).*x(i,:,t)) ]; end end这一段的最终形态会很长但我强烈建议你分步验证先不管MPS单独跑通潮流约束看看电压是否在范围、支路功率是否平衡再加入MPS变量。全线一起建模出了bug非常难找。求解器调用上除了设verbose,0还要设mipgap。我一般设成1e-3也就是最多容忍千分之一的目标值差距求解速度能快好几倍复现论文时误差完全可接受。4. 测试场景设计与结果验证33节点系统上的三组对比实验4.1 测试系统与故障场景设置我用的是经典的33节点配电网测试系统在这个复现里简称为“某33节点系统”。系统基准电压12.66kV功率基准1MVA32条支路总峰值负荷约3.8MW。我在原始拓扑上做了两处改动给节点12、25、30这三个位置预留了移动电源接入接口把MPS初始位置设为预配置阶段给出的节点12和节点26。故障场景设置成灾害持续10小时线路5-6、8-9、12-13、28-29断开导致网络分裂成三个孤岛。每个孤岛里至少有一个节点具备MPS接入条件这正是预配置阶段想达到的效果。MPS参数如下表参数数值说明额定容量500 kWh相当于中等规模储能车最大放电功率100 kW出口功率上限初始SOC90%预配置阶段已充满移动速度2节点/时段折算为每小时能经过的节点数充放电效率95%单程损耗5%负荷权重10 / 3 / 1关键节点/重要节点/一般节点4.2 三组对比实验与结果解读我设计了三个对照场景核心问题是MPS到底该不该动态调度以及动态调度比静态部署强在哪。Case A无MPS参与。所有孤岛只能依靠网络重构恢复失负荷最严重。Case B静态接入。MPS在预配置节点上全程固定放电不做任何移动。Case C动态调度。MPS按照优化结果在每小时决策是否移动、移动到哪里、放多少电。一个典型场景下的结果对比如下场景累计失负荷量(MWh)负荷恢复率(%)MPS移动次数Case A 无MPS87.4120Case B 静态接入55.2440Case C 动态调度34.1656从这个表能很直观地看出动态调度相比静态接入失负荷量下降了38%来源就是MPS在负荷高峰时段主动转移到负荷最重的孤岛。静态方案里停在节点12的车只能在北区放电哪怕南区负荷是北区两倍也过不去动态方案里车会从节点12先给北区顶住前两小时然后转移到南区支援实现“哪缺电去哪”。4.3 滚动优化与全局优化的差距我还加了Case D用滚动优化替代离线全局优化窗口宽度6小时每2小时重算一次。结果如下场景累计失负荷量(MWh)单次求解时间(s)Case C 全局离线34.1约260Case D 滚动优化35.4约12可以看到滚动优化只比全局最优多损失约4%的负荷但单次求解时间从260秒降到12秒量级完全能满足实时调度的要求。这个结论对工程实践很有参考意义论文边的全局最优解告诉你方法的上限而实际部署时滚动优化的版本才是能落地用的版本。复现的时候两组都保留图里把滚动优化的边界画出来审稿人和导师都会觉得你想得比普通复现者深一步。5. 复现过程中最容易踩的五个坑5.1 索引错位预配置输出和动态调度输入的节点编号不一致这个坑我见过不止一次。预配置阶段如果用0起始编号的矩阵运算函数处理数据输出初始节点时得到的是“0到N-1”的编号到了动态调度建模阶段如果沿用1起始编号的节点习惯两者的编号体系就对不上。现象很隐蔽求解结果里每台MPS都能正常放电但位置节点始终对不上网络拓扑画图时MPS出现在完全没接入的节点上。解决办法很简单全项目统一用1起始编号在预配置阶段的输出端加一行initNode initNode 1或读入端检查范围并写一个assert( all(initNode1) all(initNodenetwork.nb) )。这种低级错误在大型复现项目里浪费的时间比建模本身还多。5.2 移动时间离散化不足导致MPS“瞬移”如果你的移动约束只是“每时段最多选一条边移动”但网络里两个节点之间实际要经过3条边、耗时3小时那么求解器会钻空子它让MPS第一小时从A点消失第三小时在B点出现中间的2小时凭空调度。这种解看着目标函数很漂亮但实际上完全不可行。根因是移动状态和旅行时间脱节。我建议引入“在途剩余时间”变量把旅行状态显式建模开始移动时设置剩余时间每个时段结束减少1只要剩余时间大于0该车的接入变量强制为0。代码量增加不大但模型的物理合理性立刻就有了。这也是论文方不方便用文字写清楚、但复现时绕不过去的细节。5.3 量纲不统一引起的数值病态配电网里电阻、电抗的数值都很小如果用有名值线路阻抗可能是0.01量级节点负荷却是几千千瓦量级两者相差七八个数量级商用求解器的内部数值算法会被这种系数拖垮。结果表现为求解特别慢甚至报数值警告解出来的电压、功率乱作一团。解决方法是全部采用标幺值。功率基准1MVA电压基准12.66kV把负荷、阻抗、功率全部折算到0到几的范围内。这样做不仅求解快而且约束之间的尺度接近求解器数值稳定性会好一个档次。输出结果时再做一次逆变换不影响论文要的物理量。5.4 二进制变量爆炸与求解时间失控33节点、2台车、24时段这个规模求解很快。但如果你想复现更大规模比如123节点、10台车、96个时段二进制变量数量直接超过十万单纯砸求解器性能不是好办法。我常用的两个加速手段一是候选节点削减。每台车在同一时段能到达的节点是有限的因为移动速度和一个时段的时长限制了它的活动半径。在建模前把所有不可达节点的接入变量直接剔除变量能减少一半以上。二是加对称破缺约束。如果两台MPS型号相同、初始位置对称模型里存在大量对称解求解器会在等价解之间反复横跳。这时候可以按“编号更小的车初始位置更靠前”之类的规则加约束切断对称性。5.5 大M取值与约束松弛度之间的矛盾在把逻辑约束线性化的过程中大M几乎是绕不开的工具。比如“线路断开时支路功率必须为0”这类约束需要写成 (P_{ij,t} \le M z_{ij,t})。M如果取太大比如10000数值问题立刻显现如果取太小比如刚好等于支路最大值某些截断又可能把可行解切掉。我的经验法则是 (M) 取该支路允许的最大传输功率的1.05倍最多不超过1.1倍。如果你用了MILP求解器优先考虑用指示器约束或凸松弛替代大M线性化能省掉这一整类调参痛苦。DistFlow里的电压-功率乘积约束也可以用二阶锥凸松弛处理规避大M。最后再分享一个心得这类动态调度模型我第一次跑通整整花了一周后来发现最耗时的是调试约束之间的逻辑矛盾。后来改成从3节点、2台MPS、2时段的极小样例开始验证用Matlab把每一组约束单独打出来检查全部确认无误后再换到33节点系统。建议你也从这个小样例开始它会把“建模逻辑正确性”和“数据规模”这两个维度分离排查问题的速度快很多。这套框架其实还有不少扩展空间目前我默认MPS在道路网上的移动是简化成节点跳变的如果你们场景对交通拥堵敏感可以把道路通行时间建模成时段相关的参数也可以把多类型移动储能、移动发电机和固定储能同时纳入调度。后面有机会我会把这几条扩展路线逐一写出来尤其是“MPS与配电网重构联合优化”那部分坑更多但收益也更大。