配电网韧性下移动电源动态调度的Matlab复现详解

先把结论放这儿：这套复现代码跑通了，动态调度部分比预配置复杂一个量级，但效果也确实明显。下面是我从建模、代码到结果完整走一遍的经验记录，希望对正在啃韧性相关论文的朋友有点帮助。

之所以把动态调度单独拎出来写，是因为“上篇”的MPS预配置解决的是“灾害来之前，移动电源放在哪里”的问题，而“下篇”的MPS动态调度解决的是“灾害发生之后，这些移动电源往哪儿跑、怎么出力”的问题。这两个问题在数学建模上是完全不同的：预配置是一个带场景集的随机优化问题，动态调度则是一个带时间维度的序贯决策问题。你可以把预配置理解为“赛前布局”，把动态调度理解为“赛后临场指挥”——前者定的是站位，后者定的是动作，两者配合好了才能把配电网韧性真正提起来。

从复现角度讲，动态调度这块更值得花时间。原因也很简单：预配置的结果可以直接用混合整数线性规划（MILP）解出来，而动态调度引入了“移动”这个动作，模型的变量数和约束数会成倍上涨，稍不注意就会出现求解器内存爆炸、约束写漏、结果不符合物理规律这些坑。

1. 内容整体设计与思路拆解

1.1 从预配置到动态调度：两个阶段到底在配合什么

很多刚开始接触配电网韧性的朋友会有一个困惑：既然MPS预配置已经给出了初始部署位置，为什么还要再做一次动态调度？

我用自己的话解释一下。预配置阶段的本质是“基于预测的场景去做决策”。台风过境之前，我们知道的是哪几条线路大概率会断、哪些负荷节点最脆弱，但实际的故障情况往往和预测有偏差。预配置把所有可能发生的故障场景都考虑进去，给出一个不把鸡蛋放在同一个篮子里的部署方案，它的目标是让移动电源的初始位置在期望意义上最优。

但灾害发生之后，真实故障状态已经确定了，比如实际断了两条线路、某个变电站失压了。这时候如果还按预配置的方案一动不动，那MPS的位置未必是最优的。动态调度要做的就是在故障清晰之后，重新优化每个时段MPS的位置、接入节点和输出功率，让它从预配置的位置出发，移动到真正需要它的地方去。

从数学角度讲，预配置解出来的MPS初始位置，刚好可以作为动态调度的初值条件。这就是两阶段方法最自然的衔接方式：预配置给动态调度提供“起点”，动态调度把“起点”变成“终点”。

1.2 为什么选择Matlab作为复现工具

说实话，配电网韧性调度这类模型，用Python也能做，但我在复现论文代码时坚持用Matlab，主要有三个原因。

第一，Matlab的YALMIP工具箱极大降低了建模成本。论文里的目标函数和约束条件是用数学公式写的，YALMIP的建模语法几乎和数学表达式一一对应，你可以直接照着论文公式把约束一条条写进去，不用像Python的PuLP那样刻意调整数据结构和语法。

第二，Matlab自带的绘图功能非常顺手。动态调度结果要输出MPS移动路径、负荷恢复曲线、电压分布这些图，Matlab的plot、bar、gantt图都很好用，尤其是论文复现阶段需要对照原论文出图，这个优势非常明显。

第三，Matlab的面向对象（OOP）架构很适合做这种带物理实体的仿真。我把配电网节点、线路、MPS都抽象成了类，每个类里封装了各自的属性（比如MPS的容量、SOC、额定功率）和行为（移动、充电、放电），代码结构比一长串脚本清晰得多，也方便后续改参数、换算例。

1.3 动态调度模型的总体思路

开始写代码之前，我先把动态调度模型的总体框架画在纸上（建议你也这么做）。整个框架包含四层：

• 数据层：配电网拓扑数据（节点、线路、负荷）、故障场景数据（断线位置）、MPS参数（容量、功率、速度）
• 决策层：每个时段MPS在哪个节点、是否接入、输出多少功率、是否处于移动状态
• 约束层：潮流约束（DistFlow模型）、电压约束、线路容量约束、MPS时空约束、SOC动态约束
• 目标层：最小化负荷削减惩罚和MPS运行成本

这四层之间的关系是：数据层给决策层提供输入，决策层的变量代入约束层形成可行域，最后在可行域内优化目标层。代码的结构也是按这个逻辑组织的，后面会详细讲。

2. 核心模型细节与实操要点

2.1 动态调度的目标函数：不只是“少切负荷”

动态调度的目标函数，不同论文写法略有差异，但核心思路一致。我复现的这个版本用的是加权多目标：

第1项是负荷削减惩罚。每个负荷节点都有重要程度，用权重w_i表示。灾害恢复场景下，一级负荷权重最高，二级次之，三级最低。第2项是MPS运行成本，包括移动成本（烧油/耗电）和输出功率成本。目标函数把这二者统一到一个量纲下（通常都折算成钱），然后整体最小化。

一个关键点是：这里的负荷削减是一个连续变量P_{curt,i,t}，取值范围是0到该节点的实际负荷。连续变量比0-1变量好解很多，但要注意，如果想要“优先保住一级负荷”这种硬约束，就得引入0-1变量来控制切负荷顺序，这会增加求解难度，需要在论文方案和可解性之间做权衡。

2.2 MPS时空状态：动态调度的灵魂

动态调度最难的约束就是MPS的时空状态描述。原因很简单：MPS不是静止的，它在某个时段可能接入节点A供电，下一个时段可能正在去节点B的路上，再下一个时段才接入节点B。这三态切换在数学上要用一组变量和约束来精确刻画。

我用两个核心变量来表示MPS的时空状态：

• s_{m,t}：表示第m个MPS在t时段所处的状态（0=移动中，1=接入节点i运行）
• b_{m,i,t}：表示第m个MPS在t时段是否接入节点i（0-1变量）

这样，“移动中”和“接入某个节点”就是互斥的，体现在约束上就是t时段要么在移动，要么在唯一的节点接入，不能同时出现两个状态。

另外一个容易踩坑的地方是移动时间的建模。MPS从一个节点到另一个节点需要时间，这个时间取决于地理距离和移动速度，折算成时段数后，约束是：如果MPS在t时段从节点i出发去节点j，那么它在t、t+1、...、t+time(i,j)-1这期间都得处于移动状态，time(i,j)表示两节点间的行程时段数。这个“移动期内不允许接入”的约束，我一开始漏写了，导致求解结果里MPS一边移动一边给负荷供电，物理上就不可能。

2.3 潮流约束：DistFlow模型及其线性化

配电网潮流分析最常用的是DistFlow模型，它描述的是放射状配电网中功率从根节点流向末端节点的过程，递推式是：

P_{ij,t} - P_{jk,t} - r_{jk} * I_{jk,t}^2 = P_{load,k,t} - P_{gen,k,t}（有功平衡） Q_{ij,t} - Q_{jk,t} - x_{jk} * I_{jk,t}^2 = Q_{load,k,t} - Q_{gen,k,t}（无功平衡） V_{j,t}^2 = V_{i,t}^2 - 2(r_{jk}P_{jk,t} + x_{jk}Q_{jk,t}) + (r_{jk}^2 + x_{jk}^2)I_{jk,t}^2（电压降落）

这个模型是非线性的，因为I^2和V^2之间存在非线性的乘积关系。直接求解非线性规划非常慢，论文复现时一般要做处理。常见的做法有两种：

第一种是完全线性化，忽略线路损耗项I^2*r，把潮流约束转成线性方程。这个方案在配电网负荷不太重的情况下精度足够，求解速度最快。

第二种是二阶锥松弛SOCP，把I^2和V^2换元成新变量，将等式约束松弛成锥约束，精度高但求解时间更长。

我复现时先用的线性化版本，把线路损耗项按常数近似处理，求解速度非常快，结果和原论文数值对比误差在可接受范围内。如果你对精度要求高，再考虑SOCP版本，后面会讲怎么切换。

2.4 MPS的SOC与功率约束

MPS的本质是一个移动储能电池，所以SOC（荷电状态）的动态变化必须建模。约束逻辑是：

SOC_{m,t+1} = SOC_{m,t} - P_{out,m,t} * Δt / Cap_m + P_{in,m,t} * Δt * η / Cap_m

其中P_{out}是放电功率（为正），P_{in}是充电功率（为负），η是充放电效率。放电时SOC下降，充电时SOC上升。

还要限制变压器/MPS的充放电功率不能超过额定值，SOC不能超过上下限（比如20%-90%）。这些约束都要结合MPS接入节点变量b_{m,i,t}来写：只有当b_{m,i,t}=1时，MPS的功率才有非零值；b=0时功率必须为0。这种“条件性约束”用big-M法实现，M值取MPS额定功率即可。

实操中我建议把SOC约束单独写成一个函数，因为论文里经常会改动这个部分，单独封装好调试。

3. 实操过程与核心环节实现

3.1 Matlab代码结构：用OOP把每个实体变成类

先说结论：用Matlab OOP写这个项目的代码，最大的收获是“好改”。论文复现是一个反复迭代的过程，你可能今天改MPS数量，明天改故障场景，后天换配电网算例。如果用纯脚本写，这些改动会牵一发动全身；用类封装之后，大部分改动只是改参数或加方法，代码不会乱。

我采用的目录结构如下：

MPS_Scheduling/ ├── main_dynamic_scheduling.m # 主程序入口 ├── data/ │ ├── network_properties_xlsx.xlsx # 配电网拓扑与线路参数 │ ├── load_curve.mat # 负荷曲线数据 │ └── fault_scenarios.mat # 故障场景数据 ├── classes/ │ ├── Network.m # 配电网类：节点、线路、潮流计算 │ ├── MobileStorage.m # MPS类：容量、功率、SOC、位置 │ └── SchedulingModel.m # 调度模型类：约束构建、求解 ├── utils/ │ ├── load_data.m # 数据加载与预处理 │ ├── build_travel_time_matrix.m # 计算MPS移动时间矩阵 │ ├── check_feasibility.m # 求解结果物理合理性检查 │ └── plot_results.m # 结果可视化 └── results/ └── case_ieee33/

这里我特别想强调的是build_travel_time_matrix这个函数。动态调度里MPS的移动时间矩阵是核心输入，它是用节点间的实际道路距离（不是直线距离）除以MPS平均移动速度算出来的，再除以单个时段的时长向上取整。不同算例的节点地理坐标格式可能不一样，这个函数要做成通用的。

3.2 主程序运行流程

主程序main_dynamic_scheduling.m的执行流程是这样的：

第一步，加载基础数据。读入配电网拓扑（节点编号、线路阻抗）、负荷曲线、MPS参数。这里有一个细节：负荷曲线建议按标幺值存储，实际负荷=基准负荷*标幺值曲线，这样后续切换算例只需要改基准参数，曲线不用重新造。

第二步，生成故障场景。在配电网韧性研究中，故障一般指线路断开。数据里fault_scenarios.mat存储了一个结构体数组，每个元素描述一个故障场景：断哪条线路、发生在哪个时段、预计恢复时段。注意，这里的“预计恢复时段”不是固定的，动态调度需要在这个窗口内决策MPS动作。

第三步，构建优化模型。实例化SchedulingModel类，把网络类、MPS类、故障场景传进去，在类的内部方法里完成YALMIP变量创建和约束搭建。

第四步，调用求解器求解。我用的是YALMIP+Cplex/Gurobi的组合。

第五步，结果校验与可视化。先调用check_feasibility检查物理合理性（电压越限、SOC越界等），再调用plot_results画图。

3.3 关键代码段讲解：构建YALMIP变量

下面这段是我在SchedulingModel类里创建决策变量的核心代码，我把它贴出来讲一些容易理解偏的地方：

function model = build_decision_variables(model, network, mps, time_horizon, num_mps) % 创建变量 % P_line: 线路有功功率 [n_line, T] % Q_line: 线路无功功率 [n_line, T] % V_square: 节点电压平方 [n_bus, T] % P_curt: 切负荷量 [n_load, T] % P_mps_out: MPS放电功率 [num_mps, T] % P_mps_in: MPS充电功率 [num_mps, T] % SOC: MPS荷电状态 [num_mps, T] % zeta: MPS接入节点指示变量 [num_mps, n_bus, T] 0-1变量 % y_move: MPS移动状态指示变量 [num_mps, T] 0-1变量 n_line = network.num_line; n_bus = network.num_bus; n_load = network.num_load; T = time_horizon; % 连续变量 P_line = sdpvar(n_line, T, 'full'); Q_line = sdpvar(n_line, T, 'full'); V_square = sdpvar(n_bus, T, 'full'); P_curt = sdpvar(n_load, T, 'full'); P_mps_out = sdpvar(num_mps, T, 'full'); P_mps_in = sdpvar(num_mps, T, 'full'); SOC = sdpvar(num_mps, T, 'full'); % 0-1变量 zeta = binvar(num_mps, n_bus, T, 'full'); y_move = binvar(num_mps, T, 'full'); model.P_line = P_line; model.Q_line = Q_line; model.V_square = V_square; model.P_curt = P_curt; model.P_mps_out = P_mps_out; model.P_mps_in = P_mps_in; model.SOC = SOC; model.zeta = zeta; model.y_move = y_move; end

很多第一次用YALMIP的朋友会困惑sdpvar和binvar的区别，这里解释一下：sdpvar声明的是连续决策变量，binvar声明的是0-1决策变量。配电网里的切负荷量、功率、SOC是连续变量（可以是小数），接入位置、移动状态是逻辑判断，必须用binvar。如果该用binvar的地方错用sdpvar求解，模型性质就会从MILP变成LP，结果必然不对。

还要注意变量维度。我的习惯是统一用(n_variable, T)的矩阵表示每个时刻的变量值，第t列就是该变量在第t时段的值。这样约束循环的时候索引很直观。

3.4 约束构建实战：把这些式子敲进YALMIP

约束构建是整个代码里最需要细心的地方，我列举几个关键约束和它们对应的YALMIP写法。

一是潮流约束（线性化版本）：

% 线路潮流平衡（DistFlow线性化，忽略线损） % P_inj = 节点注入功率 for line_idx = 1:n_line bus_from = network.line_from(line_idx); bus_to = network.line_to(line_idx); constraints = [constraints, ... P_line(line_idx, :) == ... P_line(find(network.line_from == bus_to), :) + P_curt_load(bus_to, :) + ... % 下游负荷 ]; end

注意这里“下游”的确定是用拓扑数据预先算好的，写代码之前先把网络的邻接矩阵和上下游关系算清楚，否则循环逻辑容易出错。

二是电压约束：

% 节点电压上下限 V_min = (0.95)^2; % 下限0.95pu，转换成平方 V_max = (1.05)^2; % 上限1.05pu constraints = [constraints, V_min <= V_square <= V_max];

电压用的是平方形式，这是DistFlow模型的常见做法，直接约束V平方的上下限就行。

三是MPS的单节点接入约束：

% 每时刻MPS只能在一个节点接入 for m = 1:num_mps for t = 1:T constraints = [constraints, sum(zeta(m, :, t)) == 1 - y_move(m, t)]; end end

这个约束的逻辑是：如果y_move(m,t)=1，表示MPS在移动，此时sum(zeta)=0，即MPS没有接入任何节点；如果y_move(m,t)=0，表示MPS没在移动，sum(zeta)=1，即必须接入某个节点。这就天然包含了“移动状态和接入状态互斥”的含义。这个写法比我一开始用两个独立约束的方式简洁得多，推荐。

需要展示具体代码的话，约束部分还会包括SOC更新、充放电上限、切负荷上下限、移动路程时序约束、潮流平衡约束等。下面把SOC约束的写法也贴出来：

% SOC动态约束 for m = 1:num_mps for t = 1:T-1 constraints = [constraints, ... SOC(m, t+1) == SOC(m, t) ... - P_mps_out(m, t) * dt / cap_mps(m) / eta_discharge ... + P_mps_in(m, t) * dt * eta_charge / cap_mps(m)]; end end

充电功率取正数还是负数一定要看论文定义，最好统一成“充电为正、放电为负”，或者反过来，但一致性必须保证。我一开始没注意这个符号问题，SOC一直算不对，折腾了很久。

3.5 求解器设置与参数调优

模型建好后，求解器调用如下：

options = sdpsettings('solver', 'gurobi', 'verbose', 2, ... 'gurobi.TimeLimit', 3600, 'gurobi.MIPGap', 0.01); sol = optimize(constraints, objective, options);

有一点要留意：Gurobi/Cplex的许可证配置。YALMIP调用第三方求解器，要求求解器本身安装并添加过系统路径。如果你的机器上没装Gurobi，直接用内置求解器（默认的linprog等）解MILP非常慢，建议装一个学术版，速度快非常多。

参数调优方面，我实测下来的经验是：

• TimeLimit设成3600秒（1小时），防止个别场景陷进去出不来
• MIPGap设成1%（0.01），这是学术复现比较合适的gap，能平衡最优性和求解时间
• 如果算例规模确实很大，可以先设gap=5%跑一遍快速看结构，确认无误后再缩小gap

3.6 结果可视化：让论文图片活起来

动态调度的结果需要可视化展示，我做了一套输出图：

第一张图是负荷恢复曲线。横轴是时段，纵轴是各节点负荷恢复比例，一级负荷会有单独的高亮曲线，一眼就能看出MPS动态调度后关键负荷的供电情况。

第二张图是MPS移动路径图。在配电网拓扑图上标出MPS的移动轨迹，用颜色深浅表示时间先后，这比单纯看数据表格直观得多。实现方式是把MPS每个时段的所在节点坐标取出来，用plot和箭头连接，箭头颜色用colormap映射到时间。

画图之前，要先把每个节点的XY坐标准备好。IEEE 33节点系统有标准的坐标数据，网上很容易找到，但如果换算例，就要自己去坐标化处理。

第三张图是MPS的SOC曲线。把每台MPS的充放电状态和SOC变化画在一起，可以验证模型是否合理。比如MPS如果从预配置点出发去往故障区，那出场时SOC应该是满的，途中放电SOC下降，如果之后休整了，SOC又要在允许的范围内；SOC曲线异常往往意味着约束写错了。

4. 常见问题与排查技巧实录

4.1 问题一：求解器报错“Infeasible problem”

这是我在调试过程中最常遇到的一个错误，基本所有配电网优化新手都会踩。除非约束逻辑本身有漏洞，否则“不可行”十有八九是因为变量上下限设置矛盾，或者时段衔接约束漏了变量。

我的排查套路是：

1. 先把目标函数去掉，变成纯可行性问题（feasibility problem），只用constraints调用optimize，看能不能找到可行解
2. 如果可行域为空，逐条注释掉约束，二分法定位哪一条约束导致不可行
3. 检查每个变量的上下限是否物理一致，比如V_square下限是0.9025（0.95平方），如果写成了0.05，那约束必然不可行
4. 检查MPS移动约束和接入约束的衔接：如果MPS某个时段既不在移动也不在接入，约束就漏了

这招“先跑可行性，再跑优化”非常重要，能省下大量的排查时间。我自己调试时经常一上来就带目标函数跑，结果运气好一次过了，运气不好就不知道是约束还是目标的问题。

4.2 问题二：求解时间特别长，十几分钟不见底

动态调度模型比预配置多了时间维度和移动决策变量，求解时间倍增是很正常的，但耗时过长就需要优化。

我实测下来的经验有这几条：

1. 减少时段数。如果原论文是24时段、每时段1小时，可以先压缩成12时段、每时段2小时跑通逻辑，确认代码没问题再恢复全时段。KGT这里的经验是，减少一半时段，求解时间能降到原来的20%以下。

2. 固定一些无关紧要的变量。比如某些负荷节点不重要而且没有MPS经过，可以把切负荷量设为已知值，不用当变量优化。

3. 调整MIPGap。做灵敏度分析时用5%的gap，最终确定论文结果时再跑1%的gap。

4. 给求解器一个热启动初值。YALMIP支持通过assign给变量赋初值，然后optimize时加'gurobi.MIPStart', 1，虽然设置起来繁琐，但对大规模算例，热启动有时能把求解时间从一小时缩短到几分钟。

4.3 问题二：电压曲线出现异常的尖峰或凹陷

如果画出来的电压曲线在某个时段突然掉得很低或者升得很高，大概率不是数学问题，而是MPS在某时段接入/断开时功率阶跃导致的电压突变，也可能是潮流方程里节点注入功率表达式写错了。

排查方法：

1. 检查MPS接入或退出时段的功率值：如果这个时段MPS同时放电又充电（不满足专门约束），那么注入功率就会计算错。
2. 检查潮流方程里的负荷方向：负荷是消耗功率的，发电是注入功率的，符号反了电压就会越走越高。
3. 用evaluate函数手动算几个时段的潮流结果，跟求解器的输出比较，确认无误后再找潮流模型的问题。

4.4 问题四：SOC计算结果跟手算不一致

这个问题是我在复现过程中花时间最久的。排查下来发现，是因为我在充电效率η和放电效率η的处理上搞混了。放电时SOC下降量是放电功率乘ΔT再除以放电效率，充电时SOC上升量是充电功率乘ΔT再乘充电效率。有些论文把两个效率都取0.95，复现起来省心；但有的论文用了不对称效率（充电效率0.95、放电效率0.9），如果你统一取0.95，SOC曲线就会明显偏离论文结果。

4.5 避坑速查表

我把上面所有经验顺手整理成一个速查表，方便调试时对照：

5. 复现过程中一点体会

这篇SCI一区论文的动态调度部分，我前后折腾了两周，最大的感触是：论文里的数学公式写得再漂亮，落到代码里都会遇到无数个小坑。特别是MPS的时空状态约束，看起来就是一个简单的0-1变量，但实际建模时要考虑的边界条件非常多，“从节点A到节点B的行程中不允许接入任何负荷”这种约束，论文里可能只是一句话，代码里却要拆成好几个矩阵运算才能精确表达。

对我个人来说，这套代码跑通之后，收获的不仅仅是一个可复现的结果，更重要的是理解了“动态调度”这类序贯决策问题在配电网韧性场景里是怎么建模的。后续想扩展方向也很方便，比如把MPS换成移动储能车与固定储能结合、在目标函数里加碳排放约束，或者把单目标扩展成多目标，代码框架都能直接复用，只需要改约束和目标函数模块。

如果你也在复现类似的论文，我的建议是：不要急着抄代码，先把论文里的每个公式对应的变量维度写清楚，再动手。变量维度清晰了，代码写起来会顺畅很多，调试时遇到矩阵维度不匹配的概率也会大大降低。希望这篇记录能让你少踩几个坑。