配电网韧性提升：MPS动态调度复现全解析与Matlab+Yalmip实现

聊到我最近在做的这个复现项目之前，先讲个有点丢人的事：静态预配置部分我一周就跑通了，MPS动态调度却拖了将近一个月。那段时间几乎每天都能碰到新问题——数学模型理解偏差、Yalmip里二阶锥约束写法不对、Cplex数值警告刷屏、韧性曲线跟论文差距大。后来把所有约束拆开逐个排查，终于把整套代码稳定了下来。这中间踩过的坑，我觉得比论文本身还值得写下来。

这篇博文是复现工作的下半部分，对应的是“应急移动电源预配置和动态调度”中的MPS动态调度模块，全程基于Matlab代码实现。如果你正在复现类似的配电网韧性提升论文，或者手里已经有一套预配置结果但不知道怎么把移动电源“跑起来”，这篇应该能帮你省掉不少试错时间。我会把模型怎么建、Yalmip代码怎么搭、滚动时域怎么循环、以及哪些细节会导致复现结果跟论文对不上，一次讲清楚。

1. 动态调度与预配置到底在做什么

1.1 预配置给了什么，动态调度又补了什么

先理清两个阶段的关系。预配置解决的是“放哪”的问题：在灾害来临前，根据预测的故障场景，把有限数量的应急移动电源布放到若干关键节点。动态调度解决的是“现在该去哪”的问题：灾害发生后，故障信息逐步披露，配电网拓扑随时在变，原来预置的节点可能已经因为线路断开而失去作用，移动电源需要重新移动到其他支撑节点。

有个容易混淆的点：动态调度不是修改预配置结果，而是在预配置给出的初始位置上做时序滚动优化。我一开始就理解错了，总想着把“初始位置”和“移动路径”写进同一个模型里一起优化，结果模型直接变成一个混合整数非线性规划，Cplex根本解不动。后来把“预配置”和“动态调度”拆成两级，把预配置解出的位置当作动态调度的初值，问题规模才降到合理区间。

1.2 滚动时域控制像一个“边走边看”的计划

动态调度通常采用滚动时域控制（RHC），也叫模型预测控制的一种离散实现。说直白点：像开车导航。预配置相当于出发前根据天气预报选定了大致路线，动态调度则是导航每5分钟根据实时路况重新规划路径。

具体到代码实现，就是for t = 1:T的大循环，循环体里求解一个有限时域的优化问题，然后只执行当前时段的动作，更新系统状态，再进入下一轮求解。听上去很简单，但循环里的状态更新、移动电源位置维护、故障场景刷新，任何一个环节写错都会导致结果发散，或者收敛到一个完全没有物理意义的解。我调试过程中见过最离谱的结果：移动电源在移动途中持续给母线供电，算出来的韧性提升效果比理想调度还夸张，一看约束，移动互斥条件漏写了。

1.3 为什么许多复现者卡在这一步

我身边不少同学卡在SCI一区论文的复现上，往往不是数学公式看不懂，而是不知道“上下两篇”到底是什么关系。如果直接从动态调度开始建模型，会发现缺了一堆边界条件：MPS初始位置从哪来？可用移动电源数量多少？故障隔离时间多长？这些没定好，模型怎么调都跟论文对不上。

我的建议是：先把上篇的预配置结果跑出来，至少生成一张初始布点表，再把它作为动态调度的输入文件。别在这上面省事，后面所有对比实验都依赖这个接口的一致性。

2. MPS动态调度模型怎么建

2.1 目标函数：以“恢复负荷”还是“减少负荷削减”为目标

动态调度的目标函数一般有两种写法：一种是最大化调度周期内系统恢复的加权负荷总量，另一种是最小化整个周期内的加权负荷削减量。从数学上说两种形式等价，但我更推荐用后者，因为和韧性三角形里的“损失面积”指标能直接对应，画图时比较好解释。

具体形式可以写成：

min ∑_{t∈T} ∑_{i∈N} ω_i · (1 − α_{i,t}) · P_i^load · Δt

其中ω_i是节点i的负荷权重，α_{i,t}是t时段节点i的供电恢复比例，P_i^load是节点的峰值负荷，Δt是时间步长。负荷削减变量一般取连续型，它在0到1之间，表示部分负荷被切除的比率。

这里有个经验：权重向量不要自己随便拍脑袋定，论文里通常会给出节点负荷等级的说明，或者用表格列出每类负荷的权重。如果原文没有明确说明，我建议用负荷功率本身作为权重，至少这样结果在物理上更直观，审稿人也挑不出毛病。

2.2 MPS移动的时空转移约束

移动电源的状态必须用二进制变量刻画。常见做法是定义变量x_{m,i,j,t}，表示第m台移动电源在t时段从节点i移动到节点j。这个变量一旦置1，意味着移动电源在移动途中，不能向任何节点供电。

核心约束有三个：

• 每台移动电源每个时段最多只能处于一种状态：要么在某个节点待命并充放电，要么在移动途中。
• 移动需要时间，从i到j需要至少travel_time(i,j)个时段，不能瞬间到位。
• 如果移动电源在t时段决定从i移动到j，那么它在t+travel_time(i,j)之前不能出现在任何节点出力。

我调试时经常发现，把“移动途中不能供电”漏掉之后，优化器会让移动电源一边“瞬移”一边供电，目标函数值异常好看，但毫无实际意义。这个约束在Yalmip里通常用一组互斥不等式来表达，代码不复杂，但逻辑上非常容易遗漏。

2.3 潮流约束用DistFlow还是线性化

配电网一般是辐射状结构，经典做法是用DistFlow潮流方程。但DistFlow包含非凸二次项，直接放进动态调度模型，再加上时间维度的整数变量，整个问题会变成混合整数非线性规划，现代求解器基本无能为力。

我用的策略是：把电压幅值的平方设为新变量，把支路功率二次项用数学恒等式转成二阶锥约束，整体变成混合整数二阶锥规划，也就是MISOCP。Yalmip里写二阶锥约束很自然，Cplex可以直接处理。如果你不想自己推导松弛过程，可以用Yalmip的implies搭配norm约束来建模，但算起来会慢不少。

有的论文会直接用线性化DistFlow，也就是忽略网损、假设电压近似恒定。这个做法在动态调度阶段也很常见，因为调度本身是粗粒度决策，不需要精确到每条支路的电压降落。如果你的论文复现只要求趋势一致，线性化是性价比最高的选择；如果要求最终数值接近原文，建议用SOCP松弛。

2.4 储能约束和节点互斥

移动电源本质上是移动储能装置，所以必须考虑电池容量、充放电功率限制和荷电状态变化。常用约束是：

• SOC_{m,t+1} = SOC_{m,t} + η_charge · P_{m,t}^charge · Δt − P_{m,t}^discharge · Δt / η_discharge
• SOC有上下限，通常在10%到90%之间。
• 充放电功率不能同时为正。

真正容易出错的是“位置变量”和“充放电变量”之间的互斥关系：移动电源只有在t时段位于节点i时，才能在该节点注入有功功率。如果漏掉这个互斥约束，优化器会让同一个移动电源同时在多个节点同时供电，功率越算越多，结果完全失真。

3. Matlab代码实现骨架：Yalmip是核心利器

3.1 工具链推荐：Yalmip + Cplex

Matlab里快速搭建这类优化模型，我强烈推荐Yalmip做建模层，底层求解器用Cplex或Gurobi。Yalmip最大的好处是：你不需要手动把二阶锥约束转成求解器的标准形式，它会在内部自动识别锥结构并完成转换。底层求解器的选择上，Gurobi对MISOCP的求解速度通常优于Cplex，但Cplex在很多高校有现成授权，而且对“收敛但gap大”这类问题的数值稳定性更好。

环境配置有两点提醒：第一，Yalmip定期更新，旧版本对某些非线性约束的锥识别有bug；第二，求解器许可证和Matlab路径都要配好，最容易出的问题是你装了Gurobi但没装可用许可证，导致Yalmip直接报“No suitable solver found”。

3.2 主循环里调用滚动时域求解函数

核心结构是这样的：一个run_dispatch_rhc.m作为主脚本，内部做参数初始化、故障场景加载、初始状态赋值，然后进入滚动时域循环。循环里调用solve_rhc_step.m，这个函数负责构建Yalmip模型、调用求解器、返回第一步的控制决策。

一个精简版的结构如下：

%% 主循环：滚动时域调度 T_sim = 24; % 总时段数 H_pred = 8; % 预测时域长度 MPS_ini = pre_config; % 预配置结果，初始位置 for t = 1:T_sim % 状态刷新：加载当前时段的故障信息 fault_state = load_fault_state(t); % 构建并求解当前时刻的MPC问题 [u_dispatch, sol_info] = solve_rhc_step(t, MPS_ini, fault_state, H_pred); % 只执行第一步控制指令 act_on_decision(u_dispatch); % 更新状态：MPS位置、SOC、线路状态 MPS_ini = update_MPS_state(MPS_ini, u_dispatch, t); end

在solve_rhc_step内部，用Yalmip定义优化变量和约束，示例代码如下：

function [u, sol] = solve_rhc_step(t_now, MPS_state, fault_state, H) % 定义预测时域内的决策变量 P_dc = sdpvar(n_MPS, n_node, H); % 移动电源出力 SOC = sdpvar(n_MPS, H + 1); % 荷电状态 x_pos = binvar(n_MPS, n_node, H); % 位置0-1变量 alpha = sdpvar(n_node, H); % 供电恢复比例 obj = 0; cons = []; for k = 1:H obj = obj + sum(alpha(:,k) .* weight_vec .* P_load); % 互斥约束：同一MPS只能在一个位置 cons = [cons, sum(x_pos(:,:,k), 2) <= 1]; % 潮流和社会约束... end assign(SOC, SOC_ini); optimize(cons, obj, sdpsettings('solver', 'cplex', 'verbose', 1)); u = value([P_dc(:,:,1), x_pos(:,:,1)]); end

代码里有个小技巧：把SOC定义成H+1个时段，是为了匹配相邻时段的递推约束。我一开始用H个时段，结果最后一步SOC无处存放，Yalmip一直报维度不匹配。

3.3 如何检查求解结果是合理的

很多新手拿到sol.info为“成功”就完事了，这是一个坑。求解器返回成功只能说明“求解过程收敛”，不等于“模型建得正确”。我每次跑完都会检查三个东西：

• 是否有单台移动电源同时出现在两个位置；
• 是否有移动电源在移动途中还在出力；
• SOC是否落在上下限之间。

这三个检查各需要一个sum命令和一句assert，加起来不到10行。我在调试阶段靠这三个assert抓出了七八个建模错误，强烈建议你保留这段自检代码。

4. 复现过程中最容易卡住的四个细节

4.1 时间粒度与MPS移动时间的匹配

论文里移动时间通常给的是“分钟”，而调度时段可能按“小时”划分。如果直接拿分钟数除以60再四舍五入，会引入很大的误差。比如移动时间是30分钟，时段长度是1小时，那么30分钟到底算0个时段还是1个时段？算0个时段等于允许移动电源瞬移，算1个时段又会过度约束。

我采用的稳妥方案是：把移动时间矩阵向上取整到小时，并写进参数表统一调用。这样做虽然牺牲了一点灵活性，但至少保证了可行性。如果你需要更高精度，可以把时段长度从1小时改成15分钟，代价是求解时间会明显上升。

4.2 二阶锥约束的Yalmip写法

DistFlow松弛后，支路约束会出现形如P^2 + Q^2 ≤ S_max^2的形式。在Yalmip里写成norm([P; Q]) <= S_ub_abs，然后加上双向不等式约束。有个细节：Yalmip老版本对norm锥的识别不如新版本，如果你发现Yalmip把它当作非线性约束交给fmincon处理，速度会慢到让你怀疑人生。赶紧升级Yalmip，并确认输出中显示“Convex solver”而不是“Nonlinear solver”。

4.3 数值缩放：单位统一是门玄学

Cplex虽然号称数值稳定，但对量级相差悬殊的模型一样会崩溃。我遇到过的情况是：有功功率用MW，SOC用0到1，储能容量用kWh，结果三个变量的系数差了好几个数量级，Cplex一直在报数值警告，gap怎么都降不下去。后来把所有功率统一到kW、储能容量统一到kWh、SOC保持标幺值，再用最大功率做归一化，求解速度立刻快了几倍。详细的坑和对应方案我整理在下面：

4.4 求解时间爆炸的实用对策

动态调度最让人头疼的就是求解时间。我实测下来，IEEE 33节点系统加3台移动电源、预测时域24、总时段24的问题，直接跑MISOCP，Cplex需要几个小时才能给出满意gap。解决方案有三个，按优先级推荐：

• 降低预测时域到8到10个时段，滚动更新本来就不需要太长前瞻；
• 先解连续松弛，把性能上限算出来，再决定是否需要整数解；
• 对移动电源的总移动次数做限制，比如一个调度周期内最多移动两次，这相当于加一条割平面，效果立竿见影。

我最后用H=10 + 每台移动电源最多移动3次，整个动态调度24时段只花了不到15分钟，和论文图表基本吻合。

5. 韧性曲线怎么验证：复现结果到底对不对

5.1 韧性三角形和自己画R曲线

韧性评估最直观的图表就是“系统供电能力R(t)随时间变化曲线”。灾害时刻t0发生后，系统供电能力从100%快速下降；随着故障隔离和移动电源接入，供电能力逐渐回升；最终恢复到100%或接近100%。整个曲线围出的凹陷面积，就是韧性损失，面积越小，代表系统韧性越强。

动态调度的价值，就在于让R(t)曲线“掉得更浅、回得更快”。具体做法是：把每个时段所有节点的恢复负荷累计起来除以总负荷，得到R(t)序列，再用cumsum算面积。

5.2 对比实验怎么设置才可信

如果只是为了验证复现结果，单独一条动态调度曲线是不够的。我通常设置四个方案：

做完这个对比之后，论文里的核心结论基本都能复现出来：方案B比方案A好，方案C比方案B更好，方案D给出理论上限。如果你的曲线在方案C不能稳定优于B，说明动态调度的优势没跑出来，大概率是约束写松了或者时间窗口设置不合理。

5.3 结果对不上的时候先查什么

我复现时发现自己的方案C比论文差了一块，反复调模型都没用。最后对比发现，论文用了更细的15分钟时段，而我用1小时时段，导致移动时间的量化误差累计放大。把时段精度提上去之后，数值就基本贴合了。

所以当你发现结果对不上时，排查顺序建议是：先查时段长度和移动时间矩阵，再查故障场景的持续时间，最后才考虑是不是模型约束写错。这三个顺序我是按实际概率排的，前两个出问题的可能性远大于最后一个。

最后的几点实操提醒

多说一句数据接口的事：整个复现项目里，最消耗时间的不是建模，也不是调求解器，而是模块之间的数据传递。预配置输出、故障场景加载、MPS状态更新、韧性指标统计，多个脚本之间如果统一一套数据标准，后面会轻松很多。

我的做法是定义了一个公共参数结构体，节点数、时段数、移动时间矩阵、MPS数量、权重向量全部塞进去，所有脚本统一调用这个结构体。这样改一个参数，整套流程自动同步，再也不用担心预配置的结果在动态调度里读出来对不上号。

另外调试阶段一个小技巧：先用随机故障场景和固定随机种子跑通全套代码，再去加载论文的原始场景。固定随机种子可以保证每次调试看到的差异都来自代码改动，而不是随机噪声，这一点在穷举式调参时特别管用。

等这一整套流程跑通，你再去回头看那篇SCI一区论文，会发现模型本身并没有想象中那么高不可攀，真正的门槛全藏在那些一行就能写完、但漏掉就万劫不复的约束条件里。希望这篇能帮你少踩几个我没躲开的坑。