最近在复现一篇SCI一区期刊上关于配电网韧性提升的文章重点把其中应急移动电源Mobile Power SourceMPS的动态调度部分在Matlab里完整跑通了。这个方向现在确实很热——台风、暴雨、覆冰这些极端天气一上来配电网最容易发生大面积停电靠传统抢修恢复太慢所以很多研究都把重心转向“灾前预配置灾中动态调度”的组合拳。上篇预配置解决的是“移动电源提前停在哪里、配多大容量”下篇动态调度解决的是“灾害发生后故障点已经明确这些移动电源怎么陆续赶到灾区、优先给哪些负荷供电”。这篇博文就是要把动态调度这部分原理解透并把Matlab复现的过程记录下来。想搞配电网韧性、移动储能、应急电源优化调度的研究生或者刚从传统电力系统优化转过来的同学看完应该能少走不少弯路。复现之前我先把整个问题拆了一遍。最深的感受是这类文章表面上在写“优化算法”实际上真正的难点在于约束建模尤其是移动电源的空间转移和时间过程怎么统一到同一个模型里。这一篇我会把动态调度模型、Yalmip建模技巧、求解器调参、常见坑一次讲清楚。1. 先把问题边界划清楚动态调度到底在调度什么1.1 预配置和动态调度如何衔接从题目就能看出来“预配置”和“动态调度”是同一个研究链条里的两个环节。预配置阶段一般在灾前完成输入是灾害预测信息或者历史典型场景输出是MPS的初始部署位置、数量和容量。到了动态调度阶段假设灾害已经发生调度中心拿到故障线路、停电负荷、可用MPS位置这些信息需要决定每个MPS后续怎么移动、何时接入、给谁供电。复现的时候这两个阶段在代码层面其实是分开的。上篇的结果就是一个静态数组比如mps(:,1)表示每个MPS的初始节点编号mps_cap(m)表示容量。到了下篇动态调度这些值直接作为固定参数写进约束里。这里有一个容易忽视的衔接问题如果上篇的预配置方案是随机场景优化的那么下篇在具体故障场景下做调度时初始位置未必是“最优”的但模型必须接受这个设定。这也是很多论文里说的“两阶段决策”结构——第一阶段做决策时不依赖第二阶段的具体实现第二阶段则基于第一阶段的决策做适应性调整。我在代码里用了一个固定的初始化结构mps_init_pos [14, 25, 30]; % 预配置位置来自上篇结果 mps_capacity [500, 300, 300]; % 单位kWh应急电源容量如果你的复现目标里没有上篇代码完全可以手动设置一个合理初始位置不影响下篇调度逻辑的验证。1.2 动态调度的三个关键决策动态调度模型的核心可以拆成三个“要决策什么”第一分配给谁。每个MPS容量有限、位置不同而每个故障节点的负荷重要程度和恢复价值不同所以要把MPS和节点做匹配。第二何时到达。MPS从初始位置到目标节点需要通行时间这个时间取决于节点之间的道路距离和移动速度。早到一个小时关键负荷就少停一个小时所以在时间轴上做优化比单纯空间分配更有价值。第三供多少功率。到了节点之后MPS在同一时刻能输出的功率受容量限制而且总放电能量也受储能容量限制不是想供多少就供多少。这三层决策耦合在一个时间轴上本质上是一个时空网络流问题。我用一个生活化类比辅助理解几个水电工师傅从各自家里出发要处理多个不同优先级的故障点到达每个点需要不同时间处理每个点消耗不同工时和物料怎么排最合理。这个类比虽然不完全等价但能很快抓住MPS调度的本质——空间转移有代价时间分配有先后负荷恢复有优先级。1.3 目标函数怎么定配电网韧性提升的量化指标有很多种动态调度里最常用的是“加权负荷恢复量最大化”。这里的权重直接对应负荷等级比如医院、通信基站、供水设施的权重最高普通居民负荷权重低一些。极端灾害下无法恢复全部负荷所以优先保障重要负荷是模型的核心逻辑。有的论文会在目标函数里加一个惩罚项比如MPS移动距离的惩罚或者移动次数的惩罚目的是防止模型为了极小收益让MPS频繁移动。复现时可以先不加惩罚项跑通之后再对比加惩罚项的结果差异这本身也是一种很好的模型验证手段。目标函数的标准形式是objective -sum(omega(n) * p_load(n,t) * dt) lambda * sum(dep(...));这里omega是负荷权重dt是时间步长dep是移动弧变量lambda是移动成本系数。注意Yalmip默认是最小化所以目标函数前面加负号或者把Yalmip的sense设成minimize之后按上式写。我习惯直接写Objective -sum(sum(omega * pL)) lambda * sum(dep(:));其中pL是恢复负荷矩阵omega是列向量。2. 动态调度的核心数学模型从时间扩展图入手2.1 时间扩展图建模思路动态调度里最容易卡住的地方是MPS在空间上移动需要时间而模型中所有决策都发生在离散时段上怎么把“移动过程”表达清楚。我强烈建议用时间扩展图time-expanded network来建模。时间扩展图的思想很简单把每个时段t的配电网节点i看作一个时空节点(i,t)每个MPS在时空节点之间移动。节点i在t时段停留就相当于占用时空节点(i,t)从一个节点i移动到节点j需要T_ij个时段就相当于在时空图上从(i,t)走到(j,tT_ij)的一条弧。这样整个调度问题就变成了在时空网络上为每台MPS找一条从初始节点出发的路径路径经过的节点可以在对应时段为负荷供电。这种建模方式好处很明显一是逻辑清晰移动和供电两个状态天然分开二是约束容易写成线性形式配合二进制变量直接交给求解器三是调试直观把时空路径画出来就是MPS的完整轨迹。2.2 位置、移动与供电三类变量怎么定义复现中我定义了三组核心0-1变量和一组连续变量pos(m,n,t)MPS m在时段t是否位于节点n这是“停留可用”状态。dep(m,i,j,t)MPS m在时段t是否处于从i到j的移动过程中表示移动占用状态。s(m,n,t)MPS m在时段t是否在节点n处接入电网供电它必须约束在pos1的前提下。p(m,n,t)MPS m在时段t向节点n注入的有功功率连续非负变量。变量多了之后矩阵维度容易搞混我的习惯是先固定维度顺序为(M,N,T)也就是MPS数×节点数×时段数。下面的代码是定义变量的示例pos binvar(Nmps, Nnode, T, full); dep binvar(Nmps, Nnode, Nnode, T, full); s binvar(Nmps, Nnode, T, full); p sdpvar(Nmps, Nnode, T, full); pL sdpvar(Nnode, T, full);有个细节p(m,n,t)的变量在MPS没有到达节点n时没有意义完全可以用固定值0代替。变量越多求解越慢所以最好做一次变量裁剪。一般做法是只对MPS初始位置、故障节点以及它们附近的潜在接入节点保留p变量其他位置直接置0能显著降低模型规模。2.3 核心约束状态互斥与移动时延动态调度最核心的约束是“一个MPS在任意时刻只能处于一种状态”。具体来说每个MPS在时段t要么停留在某个节点要么在某个移动过程中不能同时出现在两个地方也不能边移动边供电。写成约束就是for m 1:Nmps for t 1:T % 状态互斥停留 移动占用 恰好一种 [con, con] ... [con, sum(pos(m,:,t), 2) sum(sum(dep(m,:,:,t)))]; Constraints [Constraints, sum(pos(m,:,t),2) sum(sum(dep(m,:,:,t))) 1]; % 只有停留可用时才能接入供电 Constraints [Constraints, s(m,:,t) pos(m,:,t)]; % 供电功率受接入状态和容量限制 Constraints [Constraints, p(m,:,t) MPS_Pmax(m) * s(m,:,t)]; end end这里有一个在复现时容易踩的坑如果MPS从节点i到节点j的移动时间T_ij大于1个时段那么dep变量表示的是“正在移动中”这个状态而不是“开始移动”的事件。这样设计的好处是状态互斥约束可以统一写成“停留移动1”而不需要关心移动是从哪个时刻开始的。但代价是需要额外加一条“移动结束时才能到达节点j”的约束% 从i出发到达j需要T_ij个时段只有完成移动后才能在j点出现 % 即 dep(m,i,j,t) 1 时pos(m,j,tT_ij) 1 % 实际写成线性表达式 for m 1:Nmps for t 1:T for i 1:Nnode for j 1:Nnode if TravelTime(i,j) 0 t TravelTime(i,j) T Constraints [Constraints, ... pos(m,j,tTravelTime(i,j)) dep(m,i,j,t)]; end end end end end同时为了防止提前到达需要在中间时段把pos锁死为0。不过因为状态互斥已经把pos和dep的关系绑定了当dep1时所有pos都必须为0这个约束自然就保证了移动期间不会出现在除终点之外的其他节点。如果你的模型里dep含义是“开始移动”而不是“正在移动”那中间时段的pos锁定就必须额外加这点一定要看清原文的定义方式。2.4 配电网潮流约束与二阶锥松弛MPS动态调度不是简单的路径规划接入配电网之后功率注入会改变潮流分布所以还要考虑潮流约束。对于辐射状配电网最常用的是DistFlow模型。以IEEE 33节点系统为例每条支路(i,j)的DistFlow方程可以写成P_ij(t) - sum(P_jk(t)) - r_ij * l_ij(t) p_mps(j,t) p_load(j,t)这里P_ij是流入支路的有功l_ij是支路电流的平方p_mps是MPS注入功率。电压约束用U_i - U_j ≥ 2(r_ij P_ij x_ij Q_ij) - (r_ij^2 x_ij^2) l_ij 来近似其中U_i是节点i电压幅值的平方。DistFlow本身包含非线性项P_ij^2 Q_ij^2直接求解很麻烦。好在大量研究已经证明在目标函数单调的前提下可以把二次等式松弛为二阶锥不等式即2P_ij^2 2Q_ij^2 (l_ij - U_i)^2 ≤ (l_ij U_i)^2用Yalmip表达这个二阶锥约束很简洁Constraints [Constraints, ... norm([2*Pij(b,t); 2*Qij(b,t); lij(b,t)-U(i,t)], 2) lij(b,t)U(i,t)];有人会问松弛之后还是原问题最优解吗这就是所谓“精确凸松弛”问题。大多数配电网辐射状且目标函数是恢复负荷最大化的场景下松弛是紧的也就是松弛解就是原问题的全局最优解。复现时如果你发现结果里某条支路的二阶锥约束不紧要检查是不是目标函数里加了不恰当的移动惩罚项或者负荷权重设置导致目标函数对潮流不敏感。判断方法很简单看每个支路对应约束左右两边的差值如果差很小比如小于1e-4说明松弛紧结果可信。3. MatlabYalmip复现的完整流程3.1 环境准备与数据组织先明确一下环境我用的Matlab R2022b加上Yalmip求解器用的Gurobi。如果没有Gurobi用Cplex或者Mosek也行如果是小规模算例甚至Cbc这种开源求解器也能跑只是速度慢一些。Yalmip的安装和配置这里不展开网上资料很多。代码组织上我分成下面几个文件main_mps_dispatch.m % 主程序 data_ieee33.m % 配电网参数与负荷数据 gen_travel_time.m % 生成节点间通行时间矩阵 build_model.m % 构建优化模型并求解 plot_results.m % 可视化结果数据准备是关键的一步。IEEE 33节点系统的支路参数、负荷数据都有公开版本但不同论文使用的基准容量、电压等级可能不同复现前一定要和原文核对清楚尤其是功率基准值和时间步长dt。我复现时用1小时作为时间步长总调度周期取24小时相当于覆盖一个典型灾后恢复日。通行时间矩阵的生成也需要提前处理。最简化方式是用节点之间的线路长度除以MPS移动速度。更精细的方式是用道路网络距离但一般论文里不会给那么详细的数据所以用线路长度代替就行。关键是TravelTime(i,j)矩阵要在建模之前就生成好并且保证对角线为0非对角线为正整数。% 生成通行时间矩阵 TravelTime zeros(Nnode, Nnode); for i 1:Nnode for j 1:Nnode if i ~ j dist norm(bus_coord(i,:) - bus_coord(j,:)); TravelTime(i,j) max(1, round(dist / MPS_speed)); end end end这里有个经验移动时间必须取整到时间步长的整数倍。如果步长是1小时而两节点之间行车需要2.4小时向上取整成3小时会让MPS“迟到”一些向下取整会导致移动时间不合理。优先向上取整保证物理可实现。3.2 模型构建核心代码build_model.m是整个复现的核心。我先把所有变量定义好然后按四类约束添加MPS运行约束、配电网潮流约束、负荷恢复约束、目标函数。MPS运行约束这一块除了前面讲的状态互斥和移动时延还有两个容易被忽略的点。第一MPS总放电能量不能超过容量。写成% MPS储能容量约束 for m 1:Nmps Constraints [Constraints, ... sum(sum(p(m,:,:))) * dt mps_capacity(m)]; end第二MPS在初始时段必须位于初始位置这是一个初始条件约束for m 1:Nmps Constraints [Constraints, pos(m, mps_init_pos(m), 1) 1]; end负荷恢复约束的关键在于恢复量不能超过原负荷需求。我一般设定一个连续恢复比例变量alpha(n,t)取值范围0~1然后pL(n,t)alpha(n,t)*load_demand(n,t)。因为无功负荷也要同步恢复所以无功恢复量用功率因数绑定。潮流部分我以DistFlow为基础按3.3节的二阶锥约束逐一添加。注意每条支路在t时段都要加一套约束循环嵌套要注意效率。Matlab里用for循环直接加约束在规模不大33节点×24时段时可接受如果系统扩大到几百节点就要考虑用矩阵方式一次性添加约束否则建模时间会暴涨。3.3 求解与结果校验求解设置方面我一般配置MIPGap为1%TimeLimit为3600秒。这里有个容易忽略的点MISOCP问题如果Gap太紧比如0.01%以下求解时间会指数级上升而对复现验证来说1%的精度已经足够判断模型逻辑是否正确。ops sdpsettings(solver,gurobi,... verbose,2,... gurobi.MIPGap,0.01,... gurobi.TimeLimit,3600,... gurobi.Threads,8); result optimize(Constraints, Objective, ops);求解完成后第一件事不是看结果而是检查求解状态。如果result.problem不为0说明模型有问题要回看约束。如果求解成功把三个关键结果画出来恢复负荷曲线、MPS轨迹、电压分布。画MPS轨迹我用的方式是对pos变量取最大值的索引figure; hold on; for m 1:Nmps [~, loc] max(pos(m,:,:), [], 2); loc squeeze(loc); stairs(1:T, loc, LineWidth, 1.8); end xlabel(时段/h); ylabel(MPS所在节点编号);这张图能很直观地看出每台MPS什么时候在哪如果曲线出现跳变比如从节点14直接跳到节点25而中间时段没有经过节点那一定是移动时延约束或状态互斥约束写错了。4. 求解过程中最容易被坑的几个地方4.1 时间步长与移动耗时矩阵的统一我最早跑模型的时候TravelTime矩阵用的是浮点数结果求解器疯狂报数值错误。后来才发现因为步长是整数小时移动时间必须取整才能匹配到整数时段。这里建议在生成TravelTime矩阵之后立刻检查assert(all(all(TravelTime round(TravelTime))), 移动时间必须为整数);另外如果某两个节点之间的移动时间大于总调度周期T那这条移动弧其实是没有意义的可以直接禁用避免增加无用的二进制变量。4.2 big-M不要取得过大在MPS接入功率约束里我用了p(m,n,t) MPS_Pmax(m) * s(m,n,t)这种形式这里MPS_Pmax就是天然的上界不需要额外再设大M。但有些约束必须用继电器形式处理比如移动时延约束里pos(m,j,tTravelTime) dep(m,i,j,t)这种逻辑约束本质是“上升沿检测”不需要大M。真正需要大M的地方是潮流约束或负荷恢复约束里某段非线性逻辑建议大M取该变量的物理上限乘以1.1不要拍脑袋填1e6。大M太大不仅数值条件差还会让MIP问题更难解。4.3 目标函数量纲与负荷权重目标函数里负荷权重omega和MPS移动惩罚lambda的量纲如果不一致可能导致优化结果完全偏向某一边。比如omega以元/kWh为单位lambda以元/次移动为单位那么lambda设为1e3可能就太小了模型会频繁移动MPS。合理做法是跑一组灵敏度分析固定omega把lambda从0开始逐渐增大观察MPS移动次数和恢复负荷的变化曲线选一个折中值。这个步骤也是论文里常说的“参数敏感性分析”复现时顺手做出来还能当额外贡献。4.4 Gurobi求解MISOCP的调参技巧Gurobi求解二阶锥MIP问题时有时候会卡在某个整数节点上半天不动。我的经验是把NumericFocus调到1或2开启Presolve的强约束检测然后MIPGap设到0.01。如果问题规模太大还可以先把潮流约束的SOCP松弛放宽一点求一个近似可行解然后再逐步收紧相当于热启动策略。具体代码如下ops.gurobi.NumericFocus 2; ops.gurobi.Presolve 2; ops.gurobi.MIPFocus 2;再补充一个实用技巧先用一个小规模算例比如IEEE 13节点、6个时段验证模型逻辑跑通之后再放大到33节点、24时段。这样排查约束错误的时间能节省一大半。5. 常见问题与排查速查表5.1 典型求解报错与处理现象可能原因解决思路Yalmip报“Solver not found”求解器路径没有配置好运行yalmiptest检查求解器是否可用确认Gurobi/Cplex已安装且Yalmip能找到模型显示“Infeasible problem”移动时延约束或状态互斥约束过紧先打开约束松弛模式逐步放宽大M或增加MPS容量定位第几类约束导致无解Gurobi警告“Numerical trouble”变量尺度差异过大常见是电压标幺值和小数约束混用全部使用标幺值U用电压标幺值的平方功率用基准功率归一化求解时间很长且Gap不下降MIP变量太多MISOCP分支搜索困难减小时间步长数量、裁剪无效移动弧、设置MIPGap为0.02或0.03先求参考解5.2 结果异常的逻辑检查求解成功但结果不合理的现象往往是建模逻辑漏洞比报错更难发现。我把检查步骤总结成一张清单第一打印每个时段的pos变量看MPS是否出现“瞬移”。如果某一台MPS在t时刻还在节点14t1时刻突然出现在节点25那说明移动时延约束没起作用。第二检查“移动中供电”是否发生。把s(m,n,t)和dep(m,i,j,t)叠加在同一条时间轴上如果同一时段既有移动又有接入供电说明状态互斥约束写错了。第三校验能量平衡。统计每台MPS所有时段注入的总能量是否严格小于等于容量。如果相等或者超出大概率是p变量的时段统计口径有问题。第四对比负荷恢复率。单时段最大恢复负荷不能超过该节点原负荷需求也不能超过MPS最大功率。超了说明pL约束或恢复比例约束写松了。第五画电压曲线。如果某个节点电压长时间越限要检查二阶锥约束的系数方向是否正确DistFlow公式里的r和x是否搞混。这里特别想强调第一点和第二点。我复现这个题目时最耗时的排查就是发现MPS会在移动过程中供电物理上完全说不通但目标函数很喜欢这种“免费供电”所以模型会钻这个空子。这类逻辑漏洞必须靠状态变量可视化来抓不能只看目标函数数值是不是合理。5.3 与预配置结果的衔接问题还有一个常见问题出现在上篇预配置和下篇动态调度的数据传递上。如果预配置输出的MPS初始位置和容量超出了动态调度允许的接入范围比如某个MPS初始位置在一个不能接入的节点那动态调度就会直接无解。我在代码里加了一个断言assert(all(ismember(mps_init_pos, candidate_nodes)), 初始位置不在可用节点集合中);同时预配置节点如果是变电站节点或者没有负荷的联络节点也要确认配电网模型里有对应的潮流变量。否则MPS到了那个节点却无处可接模型约束就会很怪。6. 写在最后复现这个题目的一点体会复现SCI一区论文最重要的是不要上来就动键盘。我最早拿到这个题目的时候以为难点在目标函数和求解器调参结果真正花时间的几乎都在约束建模尤其是MPS的时空转移逻辑。如果你也打算复现类似文章我强烈建议先把变量定义、约束类型、时间步长、节点编号全部用注释写在代码头部然后按“MPS层—负荷层—潮流层”的顺序一层层加约束每加一组约束就跑一遍看可行性这样能快速定位问题。另一个体会是这个模型后续可以扩展的方向真的很多。比如把MPS换成柴油发电车和移动储能车的混合车队考虑它们的成本差异和启动时间或者把移动时间做成随机变量用鲁棒优化或者机会约束表达交通不确定性再进一步还可以把MPS调度和抢修队伍调度统一建模实现“先复电、后修复”的协同恢复策略。这套动态调度的底层框架只要建好了扩展起来都是自然的。至少我现在好几个后续想法都是在这套代码基础上改出来的。 配电网韧性提升中的移动电源动态调度：Matlab+Yalmip建模与求解实践
最近在复现一篇SCI一区期刊上关于配电网韧性提升的文章重点把其中应急移动电源Mobile Power SourceMPS的动态调度部分在Matlab里完整跑通了。这个方向现在确实很热——台风、暴雨、覆冰这些极端天气一上来配电网最容易发生大面积停电靠传统抢修恢复太慢所以很多研究都把重心转向“灾前预配置灾中动态调度”的组合拳。上篇预配置解决的是“移动电源提前停在哪里、配多大容量”下篇动态调度解决的是“灾害发生后故障点已经明确这些移动电源怎么陆续赶到灾区、优先给哪些负荷供电”。这篇博文就是要把动态调度这部分原理解透并把Matlab复现的过程记录下来。想搞配电网韧性、移动储能、应急电源优化调度的研究生或者刚从传统电力系统优化转过来的同学看完应该能少走不少弯路。复现之前我先把整个问题拆了一遍。最深的感受是这类文章表面上在写“优化算法”实际上真正的难点在于约束建模尤其是移动电源的空间转移和时间过程怎么统一到同一个模型里。这一篇我会把动态调度模型、Yalmip建模技巧、求解器调参、常见坑一次讲清楚。1. 先把问题边界划清楚动态调度到底在调度什么1.1 预配置和动态调度如何衔接从题目就能看出来“预配置”和“动态调度”是同一个研究链条里的两个环节。预配置阶段一般在灾前完成输入是灾害预测信息或者历史典型场景输出是MPS的初始部署位置、数量和容量。到了动态调度阶段假设灾害已经发生调度中心拿到故障线路、停电负荷、可用MPS位置这些信息需要决定每个MPS后续怎么移动、何时接入、给谁供电。复现的时候这两个阶段在代码层面其实是分开的。上篇的结果就是一个静态数组比如mps(:,1)表示每个MPS的初始节点编号mps_cap(m)表示容量。到了下篇动态调度这些值直接作为固定参数写进约束里。这里有一个容易忽视的衔接问题如果上篇的预配置方案是随机场景优化的那么下篇在具体故障场景下做调度时初始位置未必是“最优”的但模型必须接受这个设定。这也是很多论文里说的“两阶段决策”结构——第一阶段做决策时不依赖第二阶段的具体实现第二阶段则基于第一阶段的决策做适应性调整。我在代码里用了一个固定的初始化结构mps_init_pos [14, 25, 30]; % 预配置位置来自上篇结果 mps_capacity [500, 300, 300]; % 单位kWh应急电源容量如果你的复现目标里没有上篇代码完全可以手动设置一个合理初始位置不影响下篇调度逻辑的验证。1.2 动态调度的三个关键决策动态调度模型的核心可以拆成三个“要决策什么”第一分配给谁。每个MPS容量有限、位置不同而每个故障节点的负荷重要程度和恢复价值不同所以要把MPS和节点做匹配。第二何时到达。MPS从初始位置到目标节点需要通行时间这个时间取决于节点之间的道路距离和移动速度。早到一个小时关键负荷就少停一个小时所以在时间轴上做优化比单纯空间分配更有价值。第三供多少功率。到了节点之后MPS在同一时刻能输出的功率受容量限制而且总放电能量也受储能容量限制不是想供多少就供多少。这三层决策耦合在一个时间轴上本质上是一个时空网络流问题。我用一个生活化类比辅助理解几个水电工师傅从各自家里出发要处理多个不同优先级的故障点到达每个点需要不同时间处理每个点消耗不同工时和物料怎么排最合理。这个类比虽然不完全等价但能很快抓住MPS调度的本质——空间转移有代价时间分配有先后负荷恢复有优先级。1.3 目标函数怎么定配电网韧性提升的量化指标有很多种动态调度里最常用的是“加权负荷恢复量最大化”。这里的权重直接对应负荷等级比如医院、通信基站、供水设施的权重最高普通居民负荷权重低一些。极端灾害下无法恢复全部负荷所以优先保障重要负荷是模型的核心逻辑。有的论文会在目标函数里加一个惩罚项比如MPS移动距离的惩罚或者移动次数的惩罚目的是防止模型为了极小收益让MPS频繁移动。复现时可以先不加惩罚项跑通之后再对比加惩罚项的结果差异这本身也是一种很好的模型验证手段。目标函数的标准形式是objective -sum(omega(n) * p_load(n,t) * dt) lambda * sum(dep(...));这里omega是负荷权重dt是时间步长dep是移动弧变量lambda是移动成本系数。注意Yalmip默认是最小化所以目标函数前面加负号或者把Yalmip的sense设成minimize之后按上式写。我习惯直接写Objective -sum(sum(omega * pL)) lambda * sum(dep(:));其中pL是恢复负荷矩阵omega是列向量。2. 动态调度的核心数学模型从时间扩展图入手2.1 时间扩展图建模思路动态调度里最容易卡住的地方是MPS在空间上移动需要时间而模型中所有决策都发生在离散时段上怎么把“移动过程”表达清楚。我强烈建议用时间扩展图time-expanded network来建模。时间扩展图的思想很简单把每个时段t的配电网节点i看作一个时空节点(i,t)每个MPS在时空节点之间移动。节点i在t时段停留就相当于占用时空节点(i,t)从一个节点i移动到节点j需要T_ij个时段就相当于在时空图上从(i,t)走到(j,tT_ij)的一条弧。这样整个调度问题就变成了在时空网络上为每台MPS找一条从初始节点出发的路径路径经过的节点可以在对应时段为负荷供电。这种建模方式好处很明显一是逻辑清晰移动和供电两个状态天然分开二是约束容易写成线性形式配合二进制变量直接交给求解器三是调试直观把时空路径画出来就是MPS的完整轨迹。2.2 位置、移动与供电三类变量怎么定义复现中我定义了三组核心0-1变量和一组连续变量pos(m,n,t)MPS m在时段t是否位于节点n这是“停留可用”状态。dep(m,i,j,t)MPS m在时段t是否处于从i到j的移动过程中表示移动占用状态。s(m,n,t)MPS m在时段t是否在节点n处接入电网供电它必须约束在pos1的前提下。p(m,n,t)MPS m在时段t向节点n注入的有功功率连续非负变量。变量多了之后矩阵维度容易搞混我的习惯是先固定维度顺序为(M,N,T)也就是MPS数×节点数×时段数。下面的代码是定义变量的示例pos binvar(Nmps, Nnode, T, full); dep binvar(Nmps, Nnode, Nnode, T, full); s binvar(Nmps, Nnode, T, full); p sdpvar(Nmps, Nnode, T, full); pL sdpvar(Nnode, T, full);有个细节p(m,n,t)的变量在MPS没有到达节点n时没有意义完全可以用固定值0代替。变量越多求解越慢所以最好做一次变量裁剪。一般做法是只对MPS初始位置、故障节点以及它们附近的潜在接入节点保留p变量其他位置直接置0能显著降低模型规模。2.3 核心约束状态互斥与移动时延动态调度最核心的约束是“一个MPS在任意时刻只能处于一种状态”。具体来说每个MPS在时段t要么停留在某个节点要么在某个移动过程中不能同时出现在两个地方也不能边移动边供电。写成约束就是for m 1:Nmps for t 1:T % 状态互斥停留 移动占用 恰好一种 [con, con] ... [con, sum(pos(m,:,t), 2) sum(sum(dep(m,:,:,t)))]; Constraints [Constraints, sum(pos(m,:,t),2) sum(sum(dep(m,:,:,t))) 1]; % 只有停留可用时才能接入供电 Constraints [Constraints, s(m,:,t) pos(m,:,t)]; % 供电功率受接入状态和容量限制 Constraints [Constraints, p(m,:,t) MPS_Pmax(m) * s(m,:,t)]; end end这里有一个在复现时容易踩的坑如果MPS从节点i到节点j的移动时间T_ij大于1个时段那么dep变量表示的是“正在移动中”这个状态而不是“开始移动”的事件。这样设计的好处是状态互斥约束可以统一写成“停留移动1”而不需要关心移动是从哪个时刻开始的。但代价是需要额外加一条“移动结束时才能到达节点j”的约束% 从i出发到达j需要T_ij个时段只有完成移动后才能在j点出现 % 即 dep(m,i,j,t) 1 时pos(m,j,tT_ij) 1 % 实际写成线性表达式 for m 1:Nmps for t 1:T for i 1:Nnode for j 1:Nnode if TravelTime(i,j) 0 t TravelTime(i,j) T Constraints [Constraints, ... pos(m,j,tTravelTime(i,j)) dep(m,i,j,t)]; end end end end end同时为了防止提前到达需要在中间时段把pos锁死为0。不过因为状态互斥已经把pos和dep的关系绑定了当dep1时所有pos都必须为0这个约束自然就保证了移动期间不会出现在除终点之外的其他节点。如果你的模型里dep含义是“开始移动”而不是“正在移动”那中间时段的pos锁定就必须额外加这点一定要看清原文的定义方式。2.4 配电网潮流约束与二阶锥松弛MPS动态调度不是简单的路径规划接入配电网之后功率注入会改变潮流分布所以还要考虑潮流约束。对于辐射状配电网最常用的是DistFlow模型。以IEEE 33节点系统为例每条支路(i,j)的DistFlow方程可以写成P_ij(t) - sum(P_jk(t)) - r_ij * l_ij(t) p_mps(j,t) p_load(j,t)这里P_ij是流入支路的有功l_ij是支路电流的平方p_mps是MPS注入功率。电压约束用U_i - U_j ≥ 2(r_ij P_ij x_ij Q_ij) - (r_ij^2 x_ij^2) l_ij 来近似其中U_i是节点i电压幅值的平方。DistFlow本身包含非线性项P_ij^2 Q_ij^2直接求解很麻烦。好在大量研究已经证明在目标函数单调的前提下可以把二次等式松弛为二阶锥不等式即2P_ij^2 2Q_ij^2 (l_ij - U_i)^2 ≤ (l_ij U_i)^2用Yalmip表达这个二阶锥约束很简洁Constraints [Constraints, ... norm([2*Pij(b,t); 2*Qij(b,t); lij(b,t)-U(i,t)], 2) lij(b,t)U(i,t)];有人会问松弛之后还是原问题最优解吗这就是所谓“精确凸松弛”问题。大多数配电网辐射状且目标函数是恢复负荷最大化的场景下松弛是紧的也就是松弛解就是原问题的全局最优解。复现时如果你发现结果里某条支路的二阶锥约束不紧要检查是不是目标函数里加了不恰当的移动惩罚项或者负荷权重设置导致目标函数对潮流不敏感。判断方法很简单看每个支路对应约束左右两边的差值如果差很小比如小于1e-4说明松弛紧结果可信。3. MatlabYalmip复现的完整流程3.1 环境准备与数据组织先明确一下环境我用的Matlab R2022b加上Yalmip求解器用的Gurobi。如果没有Gurobi用Cplex或者Mosek也行如果是小规模算例甚至Cbc这种开源求解器也能跑只是速度慢一些。Yalmip的安装和配置这里不展开网上资料很多。代码组织上我分成下面几个文件main_mps_dispatch.m % 主程序 data_ieee33.m % 配电网参数与负荷数据 gen_travel_time.m % 生成节点间通行时间矩阵 build_model.m % 构建优化模型并求解 plot_results.m % 可视化结果数据准备是关键的一步。IEEE 33节点系统的支路参数、负荷数据都有公开版本但不同论文使用的基准容量、电压等级可能不同复现前一定要和原文核对清楚尤其是功率基准值和时间步长dt。我复现时用1小时作为时间步长总调度周期取24小时相当于覆盖一个典型灾后恢复日。通行时间矩阵的生成也需要提前处理。最简化方式是用节点之间的线路长度除以MPS移动速度。更精细的方式是用道路网络距离但一般论文里不会给那么详细的数据所以用线路长度代替就行。关键是TravelTime(i,j)矩阵要在建模之前就生成好并且保证对角线为0非对角线为正整数。% 生成通行时间矩阵 TravelTime zeros(Nnode, Nnode); for i 1:Nnode for j 1:Nnode if i ~ j dist norm(bus_coord(i,:) - bus_coord(j,:)); TravelTime(i,j) max(1, round(dist / MPS_speed)); end end end这里有个经验移动时间必须取整到时间步长的整数倍。如果步长是1小时而两节点之间行车需要2.4小时向上取整成3小时会让MPS“迟到”一些向下取整会导致移动时间不合理。优先向上取整保证物理可实现。3.2 模型构建核心代码build_model.m是整个复现的核心。我先把所有变量定义好然后按四类约束添加MPS运行约束、配电网潮流约束、负荷恢复约束、目标函数。MPS运行约束这一块除了前面讲的状态互斥和移动时延还有两个容易被忽略的点。第一MPS总放电能量不能超过容量。写成% MPS储能容量约束 for m 1:Nmps Constraints [Constraints, ... sum(sum(p(m,:,:))) * dt mps_capacity(m)]; end第二MPS在初始时段必须位于初始位置这是一个初始条件约束for m 1:Nmps Constraints [Constraints, pos(m, mps_init_pos(m), 1) 1]; end负荷恢复约束的关键在于恢复量不能超过原负荷需求。我一般设定一个连续恢复比例变量alpha(n,t)取值范围0~1然后pL(n,t)alpha(n,t)*load_demand(n,t)。因为无功负荷也要同步恢复所以无功恢复量用功率因数绑定。潮流部分我以DistFlow为基础按3.3节的二阶锥约束逐一添加。注意每条支路在t时段都要加一套约束循环嵌套要注意效率。Matlab里用for循环直接加约束在规模不大33节点×24时段时可接受如果系统扩大到几百节点就要考虑用矩阵方式一次性添加约束否则建模时间会暴涨。3.3 求解与结果校验求解设置方面我一般配置MIPGap为1%TimeLimit为3600秒。这里有个容易忽略的点MISOCP问题如果Gap太紧比如0.01%以下求解时间会指数级上升而对复现验证来说1%的精度已经足够判断模型逻辑是否正确。ops sdpsettings(solver,gurobi,... verbose,2,... gurobi.MIPGap,0.01,... gurobi.TimeLimit,3600,... gurobi.Threads,8); result optimize(Constraints, Objective, ops);求解完成后第一件事不是看结果而是检查求解状态。如果result.problem不为0说明模型有问题要回看约束。如果求解成功把三个关键结果画出来恢复负荷曲线、MPS轨迹、电压分布。画MPS轨迹我用的方式是对pos变量取最大值的索引figure; hold on; for m 1:Nmps [~, loc] max(pos(m,:,:), [], 2); loc squeeze(loc); stairs(1:T, loc, LineWidth, 1.8); end xlabel(时段/h); ylabel(MPS所在节点编号);这张图能很直观地看出每台MPS什么时候在哪如果曲线出现跳变比如从节点14直接跳到节点25而中间时段没有经过节点那一定是移动时延约束或状态互斥约束写错了。4. 求解过程中最容易被坑的几个地方4.1 时间步长与移动耗时矩阵的统一我最早跑模型的时候TravelTime矩阵用的是浮点数结果求解器疯狂报数值错误。后来才发现因为步长是整数小时移动时间必须取整才能匹配到整数时段。这里建议在生成TravelTime矩阵之后立刻检查assert(all(all(TravelTime round(TravelTime))), 移动时间必须为整数);另外如果某两个节点之间的移动时间大于总调度周期T那这条移动弧其实是没有意义的可以直接禁用避免增加无用的二进制变量。4.2 big-M不要取得过大在MPS接入功率约束里我用了p(m,n,t) MPS_Pmax(m) * s(m,n,t)这种形式这里MPS_Pmax就是天然的上界不需要额外再设大M。但有些约束必须用继电器形式处理比如移动时延约束里pos(m,j,tTravelTime) dep(m,i,j,t)这种逻辑约束本质是“上升沿检测”不需要大M。真正需要大M的地方是潮流约束或负荷恢复约束里某段非线性逻辑建议大M取该变量的物理上限乘以1.1不要拍脑袋填1e6。大M太大不仅数值条件差还会让MIP问题更难解。4.3 目标函数量纲与负荷权重目标函数里负荷权重omega和MPS移动惩罚lambda的量纲如果不一致可能导致优化结果完全偏向某一边。比如omega以元/kWh为单位lambda以元/次移动为单位那么lambda设为1e3可能就太小了模型会频繁移动MPS。合理做法是跑一组灵敏度分析固定omega把lambda从0开始逐渐增大观察MPS移动次数和恢复负荷的变化曲线选一个折中值。这个步骤也是论文里常说的“参数敏感性分析”复现时顺手做出来还能当额外贡献。4.4 Gurobi求解MISOCP的调参技巧Gurobi求解二阶锥MIP问题时有时候会卡在某个整数节点上半天不动。我的经验是把NumericFocus调到1或2开启Presolve的强约束检测然后MIPGap设到0.01。如果问题规模太大还可以先把潮流约束的SOCP松弛放宽一点求一个近似可行解然后再逐步收紧相当于热启动策略。具体代码如下ops.gurobi.NumericFocus 2; ops.gurobi.Presolve 2; ops.gurobi.MIPFocus 2;再补充一个实用技巧先用一个小规模算例比如IEEE 13节点、6个时段验证模型逻辑跑通之后再放大到33节点、24时段。这样排查约束错误的时间能节省一大半。5. 常见问题与排查速查表5.1 典型求解报错与处理现象可能原因解决思路Yalmip报“Solver not found”求解器路径没有配置好运行yalmiptest检查求解器是否可用确认Gurobi/Cplex已安装且Yalmip能找到模型显示“Infeasible problem”移动时延约束或状态互斥约束过紧先打开约束松弛模式逐步放宽大M或增加MPS容量定位第几类约束导致无解Gurobi警告“Numerical trouble”变量尺度差异过大常见是电压标幺值和小数约束混用全部使用标幺值U用电压标幺值的平方功率用基准功率归一化求解时间很长且Gap不下降MIP变量太多MISOCP分支搜索困难减小时间步长数量、裁剪无效移动弧、设置MIPGap为0.02或0.03先求参考解5.2 结果异常的逻辑检查求解成功但结果不合理的现象往往是建模逻辑漏洞比报错更难发现。我把检查步骤总结成一张清单第一打印每个时段的pos变量看MPS是否出现“瞬移”。如果某一台MPS在t时刻还在节点14t1时刻突然出现在节点25那说明移动时延约束没起作用。第二检查“移动中供电”是否发生。把s(m,n,t)和dep(m,i,j,t)叠加在同一条时间轴上如果同一时段既有移动又有接入供电说明状态互斥约束写错了。第三校验能量平衡。统计每台MPS所有时段注入的总能量是否严格小于等于容量。如果相等或者超出大概率是p变量的时段统计口径有问题。第四对比负荷恢复率。单时段最大恢复负荷不能超过该节点原负荷需求也不能超过MPS最大功率。超了说明pL约束或恢复比例约束写松了。第五画电压曲线。如果某个节点电压长时间越限要检查二阶锥约束的系数方向是否正确DistFlow公式里的r和x是否搞混。这里特别想强调第一点和第二点。我复现这个题目时最耗时的排查就是发现MPS会在移动过程中供电物理上完全说不通但目标函数很喜欢这种“免费供电”所以模型会钻这个空子。这类逻辑漏洞必须靠状态变量可视化来抓不能只看目标函数数值是不是合理。5.3 与预配置结果的衔接问题还有一个常见问题出现在上篇预配置和下篇动态调度的数据传递上。如果预配置输出的MPS初始位置和容量超出了动态调度允许的接入范围比如某个MPS初始位置在一个不能接入的节点那动态调度就会直接无解。我在代码里加了一个断言assert(all(ismember(mps_init_pos, candidate_nodes)), 初始位置不在可用节点集合中);同时预配置节点如果是变电站节点或者没有负荷的联络节点也要确认配电网模型里有对应的潮流变量。否则MPS到了那个节点却无处可接模型约束就会很怪。6. 写在最后复现这个题目的一点体会复现SCI一区论文最重要的是不要上来就动键盘。我最早拿到这个题目的时候以为难点在目标函数和求解器调参结果真正花时间的几乎都在约束建模尤其是MPS的时空转移逻辑。如果你也打算复现类似文章我强烈建议先把变量定义、约束类型、时间步长、节点编号全部用注释写在代码头部然后按“MPS层—负荷层—潮流层”的顺序一层层加约束每加一组约束就跑一遍看可行性这样能快速定位问题。另一个体会是这个模型后续可以扩展的方向真的很多。比如把MPS换成柴油发电车和移动储能车的混合车队考虑它们的成本差异和启动时间或者把移动时间做成随机变量用鲁棒优化或者机会约束表达交通不确定性再进一步还可以把MPS调度和抢修队伍调度统一建模实现“先复电、后修复”的协同恢复策略。这套动态调度的底层框架只要建好了扩展起来都是自然的。至少我现在好几个后续想法都是在这套代码基础上改出来的。