做电力系统韧性研究的同行应该都绕不开应急移动电源MPS这个方向。特别是最近几年极端灾害下的配电网韧性提升成了热门而MPS凭借“可移动、可调度、即插即用”的特点成了灾前预配置和灾中应急响应的核心手段。我这次复现的是一篇SCI一区论文的下半部分——MPS动态调度对应的Matlab代码已经跑通。这篇文章我会把动态调度的建模思路、代码实现细节和调试过程中踩过的坑全部摊开来讲给正在复现或者打算做类似方向的朋友一个参考。坦白说动态调度比上篇的预配置要复杂一个量级。预配置好歹是“一次性决策”把MPS数量、容量、初始位置定下来就完事动态调度则要在故障发生后持续滚动优化每隔一段时间重新计算一次MPS该往哪走、在哪接入、出多少力。这个“滚动”的过程牵涉到时间、空间、能量三个维度的耦合再加上配电网潮流和拓扑约束整个模型写出来是个不折不扣的混合整数规划。我复现时前后花了三周时间中间踩了不少坑今天把完整的思路和可用的Matlab代码框架分享出来。1. 动态调度问题建模从“预配置”到“实时响应”的关键逻辑1.1 为什么动态调度比预配置更麻烦上篇预配置解决的是灾前“把MPS放在哪里、配多少容量”的问题本质是一个随机规划故障场景固定决策只做一次。而动态调度是灾中实时的MPS动作规划——故障发生后线路状态不断变化MPS需要从初始位置出发沿着路网移动在合适的时间接入节点并调整输出功率。这里面多出了“时间维”和“空间维”的耦合决策变量从静态的“在哪”变成了动态的“何时在哪、何时移动、何时出力”。如果只做预配置MPS灾后的行驶路径和接入时刻都是事先假设的遇到实际故障信息更新就容易失效。动态调度每过一个时段就要重新优化把最新的故障状态和系统潮流信息反馈进来相当于在滚动时域的框架里反复求解同一个混合整数优化问题。这也是为什么很多论文里MPS动态调度要用模型预测控制MPC的原因——不需要等整个故障全过程已知每个时段只对未来有限时域做决策滚动推进。实际运行时故障信息是一个小时接一个小时更新的MPC天然适配这种“走一步看一步”的应急场景。1.2 目标函数与韧性指标怎么量化复现论文时最重要的一件事是先搞清它定义的韧性指标是什么。我复现的这篇论文用的是“负荷恢复率”的面积积分——把每个时段的恢复负荷功率累计起来除以理想状态下应恢复的总负荷能量得到一个介于0到1之间的韧性指标。动态调度的目标就是在每个滚动时域内最大化这个累计恢复负荷的加权和权重可以取该时段对用户的重要程度。在Matlab里目标函数写成objective -sum(sum(W .* (1 - L_shed))); % 最小化失负荷L_shed是切负荷比例变量我这里用的是最小化切负荷的等价形式。权重矩阵W如果是常数就相当于追求恢复总能量最大如果按一天中不同时段给不同权重可以体现负荷优先级。真实场景下医院、通信基站这样的重要负荷权重高普通居民负荷权重低这个权重设置会直接影响MPS的调度路径。我的调试经验是先跑一组所有权重都等于1的算例确认恢复曲线形状合理再换加权权重对比这样能更快定位目标函数的问题。2. MPS动态调度的核心约束与变量设计2.1 时间-空间-能量三维耦合动态调度最难处理的不是潮流方程而是MPS的移动-接入-出力三件事之间的时序耦合。MPS在t时刻的位置决定了它能否在t1时刻给某个节点供电从节点i移动到节点j需要时间ΔT在移动过程中MPS不能出力接入节点后它的出力受剩余电量限制。这三层关系必须完整建模少一个约束求解器就会给出荒谬的结果——比如前一刻还在节点A下一时刻同时出现在节点B出力。我把位置变量设成三维数组X(N_node, N_time, N_mps)取值0或1表示第m个MPS在t时刻是否位于节点i。“同一时刻只能在一个节点”的约束用YALMIP写就是% 每个MPS每个时段只能在一个节点 for t 1:N_time for m 1:N_mps Constraints [Constraints, sum(X(:, t, m)) 1]; end end移动约束需要额外引入一个“移动状态”变量M_move(m, t)表示MPS在t时段是否正在移动。如果M_move为1则该时段MPS不能在任何节点出力。这里有一个常见的建模方法定义移动开始变量Start(m,i,j,t)为1表示第m个MPS在t时段从节点i出发去节点j同时该MPS需要经过ceil(τij/dt)个时段才能到达。为了代码简洁可以把目标节点的位置变量固定在经过时间之后才允许为1。更工程化的做法是引入“中途状态”节点让MPS在移动期间始终处于一个虚拟的移动节点中这样位置约束统一写为sum(X)等于1或0即可不用单独处理移动标志。能量约束要同时考虑MPS的电池容量和充放电功率限制。MPS输出功率P_mps不能超过额定功率P_rated且SOC要保持在安全区间比如[0.1, 0.9]。SOC的递推关系是SOC(t1) SOC(t) - P_mps(t) * dt / E_capacity注意如果MPS支援负荷是放电P_mps为正。这个递推在滚动时域里必须带入下一时段的初始SOC不然每次优化都从满电开始结果会过度乐观。我早期复现时就没带上一次SOC导致调度方案在第一个时段把电全部放光后面的时段全部失负荷算出来的韧性指标虚高改过来之后才正常。2.2 配电网重构与MPS移动路径的协同很多复现动态调度的论文除了MPS调度还会同时考虑配电网的拓扑重构——即利用联络开关改变馈线结构把失电区域转接到其他供电路径上。MPS和重构是协同关系重构能恢复一部分负荷MPS能补充剩余缺口。但两者的决策耦合在一起会让模型变成非常难解的混合整数规划。我在代码里把线路状态变量L_sw(N_branch, N_time)引入和MPS位置变量同步优化。辐射状拓扑约束用经典的“虚拟流法”或者生成树约束需要保证每个时段网络都保持辐射状。这个约束在Matlab中实现起来比较繁琐常见做法是引入辅助父节点变量要求每个供电节点只有一个父节点且没有环。如果只做MPS动态调度不重构可以跳过这部分但如果想完整复现SCI论文重构约束基本绕不开。在IEEE 33节点系统上我用的辐射状约束写法是定义潮流方向变量F_b(i,j,t)表示虚拟功率流根节点变电站注入虚拟功率其余节点流出负荷需求同时每条闭合线路上的虚拟潮流不能超过一个大M值。这样生成的树约束能保证每个电气岛都是辐射状结构。加了重构约束之后模型规模会翻倍但恢复效果提升通常在5%到10%之间所以论文里一般都会包含。3. Matlab代码实现从数学公式到可运行程序的关键步骤3.1 代码总体架构与数据结构设计我建议把整个动态调度工程拆成五部分数据读取、场景生成、参数设置、滚动时域主循环、结果分析。不要把所有代码塞进一个脚本里否则后期调试定位问题会非常痛苦。我在复现时按这个思路做的每个部分一个函数主脚本只负责调用。数据结构上用struct存系统参数比如System.baseMVA、System.branch、System.bus用cell数组存多个MPS的属性包括初始位置、容量、额定功率、移动速度。滚动时域里每一轮优化都会更新当前故障信息所以把故障场景也单独存成struct。这样改参数、加场景都很方便不会牵连主循环。一个典型的主脚本结构如下%% 初始化 mpc loadcase(case33); System extract_system(mpc); % 提取节点、支路、负荷数据 MPS define_mps(); % 定义MPS参数 Fault generate_scenario(System); % 生成故障场景 %% 滚动时域主循环 for k 1:N_time-1 [Constraints, obj, init_soc] build_model(k, H, System, MPS, Fault); ops sdpsettings(solver, gurobi, verbose, 1, mipgap, 0.001); optimize(Constraints, obj, ops); % 执行第一个时段的决策 record_action(k, value(X), value(P_mps)); end3.2 核心代码片段讲解动态调度主循环与求解器配置滚动时域主循环的核心逻辑是在每个时刻k读取当前故障状态构建从k到kH的预测时域优化模型求解然后只执行第一时段的控制动作其他时段的解丢弃。这叫滚动时域控制也叫MPC。主循环代码大致是for k 1:N_time-1 % 构建当前时刻的预测模型 [Constraints, obj, init_soc] build_model(k, H, System, MPS, Fault(k)); % 求解 ops sdpsettings(solver, gurobi, verbose, 2, mipgap, 0.001); optimize(Constraints, obj, ops); % 提取第一个时段的决策结果 X_opt value(X); P_opt value(P_mps); % 记录结果更新时间 end这里的build_model函数是你自己写的用来定义YALMIP变量、目标、约束并返回初始SOC状态。要注意的是YALMIP变量在循环里重复定义时必须用assign或者重新创建变量不要让上一轮的变量定义残留到下一轮否则求解规模会不断膨胀。build_model函数内部需要定义三类基本变量位置变量X二进制、出力变量P_mps连续、切负荷比例L_shed连续。我贴一个简化版的变量定义和约束构建过程function [Constraints, obj, init_soc] build_model(k, H, System, MPS, Fault) N_node System.N_node; N_branch System.N_branch; N_mps length(MPS); t_horizon k:kH-1; % 定义变量 X binvar(N_node, H, N_mps, full); P_mps sdpvar(N_node, H, N_mps, full); L_shed sdpvar(N_node, H, full); % 1表示完全切负荷0表示恢复 Constraints []; % 位置唯一性约束 for m 1:N_mps for t 1:H Constraints [Constraints, sum(X(:, t, m)) 1]; end end % SOC递推约束 for m 1:N_mps soc MPS(m).SOC_init; for t 1:H Constraints [Constraints, ... MPS(m).SOC_min soc MPS(m).SOC_max]; soc soc - sum(P_mps(:, t, m)) * System.dt / MPS(m).E_capacity; % 出力上下限 Constraints [Constraints, ... sum(P_mps(:, t, m)) MPS(m).P_rated * sum(X(:, t, m))]; Constraints [Constraints, ... sum(P_mps(:, t, m)) 0]; end end % 功率平衡与切负荷约束、潮流约束、重构约束等略 obj -sum(sum(System.Weight .* (1 - L_shed))); init_soc soc; % 返回当前时域结束时的SOC作为下一轮的初始值 end注意上面代码里“约束(略)”的部分实际工程中包括DistFlow潮流方程、线路容量约束、电压约束以及重构约束。完整代码里这些约束是用YALMIP的sdpvar和binvar组合构建的篇幅很长但核心逻辑就是按照论文公式逐条翻译。3.3 参数设置与场景生成如何复现论文中的算例论文算例通常给的是IEEE 33节点负荷数据和线路参数在Matpower里都有现成的。我从Matpower的case33中提取数据再根据故障场景随机断开若干条线路。场景生成我建议用两种方式一是论文里写死的几个典型极端场景二是随机生成100个场景求期望值。前者用来对齐论文结果后者用来做敏感性分析。场景中的故障信息包括故障线路、故障时段、预计修复时间。动态调度里故障修复后系统可以恢复正常所以故障时段是有限的。我把Fault(k)设计成这样一个structFault(k).outage_branchFault(k).end_time。MPS的初始位置来自上篇预配置的结果如果没有预配置代码可以手动指定几个关键节点作为MPS初始位置。一个比较实用的参数设置技巧是把时间步长dt设为1小时预测时域H设为6小时这样可以平衡求解速度和调度效果。如果时间步长太长MPS移动时间难以准确离散太短则优化变量暴增求解器很容易卡死。我实测下来33节点系统加3个MPS加24小时总时段用Gurobi求解单次滚动优化大约需要20秒到一分钟可以接受。4. 常见问题与调试心得让复现少走三个月弯路4.1 求解器安装与YALMIP配置细节YALMIP和Gurobi的配置是我见过复现场最常见卡壳的地方。首先YALMIP要从官网或GitHub下载最新版加入MATLAB路径。Gurobi要装对应版本并获取license。装完后在MATLAB里运行yalmiptest看到gurobi那一行显示successful才算配置好。如果yalmiptest里gurobi显示not installed多半是Gurobi的MATLAB接口路径没加对。Gurobi安装目录下的matlab文件夹要加入MATLAB路径同时确保Gurobi的版本跟MATLAB版本兼容。另外Gurobi是命名冲突的大户——如果你的当前工作区或者路径里有其他变量叫gurobi也会干扰调用。我建议在脚本开头加一段路径配置addpath(C:\gurobi1000\win64\matlab); addpath(C:\YALMIP-master);然后运行yalmiptest验证。如果求解器状态是empty或者failed先检查license是否正常再检查路径。这一步卡两周的人都不少见所以单独拿出来提醒。4.2 移动时间离散化与路径约束的坑这部分是我踩得最深的一个坑。一开始我把MPS移动简化为“一步到位”——只要在t时刻决定从i移动到j就假设t1时刻已经在j。结果模型倒是能解出来但结果完全不可用因为实际移动时间可能超过一个时段。后来我把移动时间按离散时段取整引入“中途节点”来模拟移动过程中的位置。更精确的做法是用时间扩展图把每个节点按时间维度复制成多个状态路径约束变成时间上的流约束。效果很好但变量数目会显著增加。实际操作中如果MPS从一个节点到另一个节点需要1.5个时段我建议向上取整为2个时段在这2个时段内MPS只能处于“移动中”状态不能出力和接入。这样虽然有点保守但保证了解可行性和实际物理一致性。我在代码里具体实现是增加了一个虚拟节点编号为N_node1表示“移动中”。位置约束改为每个MPS每个时段要么在某个实际节点要么在虚拟节点。当MPS从i到j需要d个时段时就强制它在虚拟节点停留d-1个时段然后再允许出现在j% 示例设置MPS从节点i到节点j的移动时间d for t k:kH-1 % 如果t时刻在i且计划去j那么td-1之前不能在j Constraints [Constraints, ... X(j, t d - 1, m) 1 - sum(X(i, t, m))]; end这种写法有一点保守但胜在简单可靠。4.3 韧性指标的收敛性与结果分析技巧韧性指标计算出来是个标量但只看数字很容易出问题。我在复现时发现不同时段的恢复曲线形状比单一指标更重要。论文里一般给的是负荷恢复曲线图横轴时间、纵轴恢复比例。你的代码里要输出这个曲线和论文的曲线对比趋势是否一致如果差很多说明目标函数权重或者约束有问题。还有一个技巧在滚动时域中把每个时段的“已实现负荷恢复量”累积起来同时把预测时域里的“计划恢复量”也画出来两条线之间的gap能反映调度保守程度。gap越大说明环境不确定性影响越大。论文中往往会对不同MPS容量、不同预配置方案做对比曲线复现时保留这些对比逻辑可以快速验证你的代码是否复现了论文的核心结论。我自己的习惯是每次跑完结果后先不急着分析指标而是画出三个关键图第一张是MPS位置随时间的变化热力图第二张是各时段MPS出力柱状图第三张是负荷恢复率曲线。位置热力图能一眼看出调度路径是否有跳变出力柱状图能看出MPS是否频繁充放电恢复率曲线能看出整体趋势。这三张图一旦符合物理直觉再去看数字指标基本就不会出大错。5. 从MATLAB结果到论文绘图复现工作最后一步5.1 把优化结果整理成可复用的表格格式论文里通常有MPS调度结果表包含每个时段的接入节点、功率、SOC值。我建议在MATLAB里把优化结果存成一个结构体Result然后用writetable输出成csv方便后续生成表格。代码片段如下T table(); T.Time (1:N_time); T.MPS1_Node squeeze(X_opt(?, :, 1)); % 根据你的变量维度调整 T.MPS1_Power squeeze(P_opt(?, :, 1)); writetable(T, mps_schedule.csv);实际上如果N_node33X_opt是一个3维数组提取每个MPS所在节点需要先找位置变量的最大索引用find函数。我写了一个小函数function node extract_node(X_matrix, t, m) % X_matrix是N_node x N_time x N_mps的数值矩阵 idx find(X_matrix(:, t, m) 0.5); if isempty(idx) node NaN; % 移动中 else node idx(1); end end注意值取出来是0/1浮点数有时会落在0.49或0.51这种阈值附近所以判断临界要大于0.5不要用等于1。这也是从Gurobi返回的精度问题导致的一个小坑。5.2 一个提高调试效率的小技巧先跑确定性场景我复现时最大的教训是不要一上来就跑完整随机场景否则模型不可行时你根本不知道是哪个约束引起的。我建议先固定一个简单故障场景比如只有一条支路发生故障MPS数量设成1预测时域设成3然后逐步放大。任何一项参数改动导致结果突变时都能通过回溯定位问题。另一个技巧是先去掉非关键约束比如电压幅值约束可以先放宽等MPS调度逻辑跑通后再把电压约束加回来。如果加了之后模型不可行说明约束参数设置有问题通常是大M值不够大或者初始SOC设置不合理。先把大M值设成1000再根据实际潮流结果缩小到合适范围。最后提醒一下YALMIP的变量定义顺序会影响求解效率。尽量把二进制变量集中定义连续变量放在后面这样Gurobi预处理时能识别出结构。如果二进制变量分散在大量连续变量之间求解时间可能增加好几倍。我实测同一个模型调整变量定义顺序后求解时间从90秒降到了25秒效果非常明显。复现SCI论文确实是个磨人的过程尤其是动态调度这种混合整数规划公式推导看着简单落成代码到处都是细节。我上面分享的这些经验都是我踩过坑之后沉淀出来的希望能帮你少走一些弯路。如果你也在复现类似论文卡在某个环节可以先对着这几条自查一遍——八成是移动时间离散、SOC递推或者求解器配置这几个地方出了问题。