要说这次复现经历最让我头疼也最让我兴奋的就是这篇配电网韧性提升论文中MPS动态调度部分。上篇已经把应急移动电源的预配置问题讲透了也就是在灾害来临前怎么把移动电源提前放到最有利的位置而下篇的MPS动态调度才是真正考验模型和算法功力的地方——如何在灾害发生后的连续时段内让这些移动电源在配电网里跑起来最大化恢复失电负荷。这两者合在一起才构成了完整的预配置动态调度闭环。而这篇博文我主要聚焦在下篇的MPS动态调度Matlab实现上把自己踩过的坑、摸索出的调试技巧、以及最终能跑通的代码逻辑全部摊开来聊。先说一句公道话这个问题的难点不在于动态调度四个字本身而在于MPS应急移动电源有着双重属性它既是电源又是可移动的储能。你得同时建模它的空间位置变化、荷电状态变化、以及接入电网后的功率注入行为。三种时间尺度耦合在一起再加上配电网的辐射状拓扑约束模型复杂度直接拉满。我在复现过程中最大的体会就是模型不是越复杂越好而是要把关键约束用最稳健的方式表达出来否则求解器分分钟给你脸色看。1. 整体设计与问题建模为什么MPS动态调度这么难1.1 配电网韧性提升的核心逻辑配电网韧性简单说就是配电网遭遇极端事件台风、冰灾、地震后能够感知、抵御、吸收并快速恢复的能力。传统的配电网规划主要考虑正常工况和N-1故障但极端灾害往往是多个设备同时故障、甚至大范围停电这时候常规备用手段就不够用了。应急移动电源Mobile Power SourceMPS作为临时电源可以灵活接入配电网的关键节点为失电负荷供电是提升韧性的重要手段。我复现的这篇SCI一区论文其亮点在于将MPS的预配置和动态调度纳入统一优化框架。预配置解决的是灾害前MPS放在哪的问题动态调度解决的是灾害后MPS怎么走、怎么用的问题。上篇文章已经介绍了预配置部分本文着重讲动态调度——也就是在灾害发生后的恢复期内如何通过动态调整MPS的位置和出力使得整体恢复效果最优。1.2 动态调度模型的目标与约束拆解MPS动态调度的目标通常是最小化失电负荷总量或者最大化恢复供电的负荷价值也可以加权重表示关键负荷优先恢复。在我复现的这个模型里目标函数包含两部分MPS在节点接入后的供电收益恢复负荷的权重减去MPS移动和运行的相关成本。约束条件则复杂得多主要分成四类MPS荷电状态约束MPS本身储能有限放电会减少SOCState of Charge可能还可以接入充电点补充电量所以要追踪每个时段的SOC变化。MPS移动约束MPS从一个节点移动到另一个节点需要时间且同一时刻只能在一个节点接入移动路径要符合配电网拓扑比如只能沿线路走。配电网运行约束包括节点功率平衡、线路潮流约束、电压上下限等通常用DistFlow或线性化潮流方程表达。时间耦合约束动态调度是一个多时段决策问题MPS在某时段的位置、出力和SOC会影响到后续时段所以需要跨时段的状态耦合。如果把MPS的数量记为M网络节点数为N时段数为T那么决策变量就包含了MNT个二进制接入位置变量加上连续的出力变量和SOC变量问题规模随着T的增长会变得相当可敬。因此几乎所有的这类论文都会采用混合整数线性规划MILP来建模再用成熟的商业求解器求解。这也就要求所有非线性项都必须线性化。1.3 为什么用MILP而不是启发式算法有一些人会问这种动态调度问题用遗传算法、粒子群这类智能优化算法是不是更简单我的经验是如果只是做一个简化版的仿真智能算法确实能很快出结果但它有两个致命弱点一是无法保证收敛到全局最优二是每次运行结果可能不一样这对于想复现论文数据的人来说就是灾难。MILP模型配合商业求解器如GUROBI、CPLEX则能给出确定性的全局最优解只要建模正确结果可复现性极强。当然前提是你得把模型约束写对不然求解器也会给出最优但其实违反物理规律的解。后面我会专门讲几个我踩过的坑。2. Matlab实现核心细节工具选择与代码架构2.1 测试系统与参数设置我采用的是配电网领域最常用的IEEE 33节点系统基准电压12.66kV基准功率1MVA总负荷约3.7MW。在Matlab中可以用matpower读取标准数据也可以直接写成结构体。为了方便做动态调度我把一天分成24个时段也可以分96个时段但求解时间会急剧上升假设灾害在第1时段发生造成若干线路断开形成孤岛。MPS参数设置如下表所示参数数值说明MPS数量2台每台容量500kWh额定功率200kW放电最大功率SOC初始值0.9灾害发生时的荷电状态SOC下限/上限0.1 / 0.95保护电池移动速度5km/h考虑实际道路距离单位移动成本0.5 /km主要算燃油/损耗需要注意的是移动速度这个参数很关键。我一开始设得太快30km/h导致求解器总是让MPS满地图乱跑结果看似恢复了更多负荷但实际上忽略了移动路径上可能存在的道路损坏问题。后来我把速度调低并增加了同一节点恢复后MPS不能随意离开的约束结果才合理。2.2 求解工具链YALMIP GUROBI我在Matlab中的建模首选是YALMIP工具箱它能把MILP建模过程变得非常优雅。求解器我用的是GUROBI因为它在处理大规模MILP时速度明显优于CPLEX至少在我的问题上如此。如果你没有GUROBI的许可证也可以用Matlab自带的intlinprog但求解速度会慢一个数量级而且在大规模问题上容易内存溢出。YALMIP的安装很简单去GitHub下载并添加到路径即可。GUROBI则需要注册学术许可证安装后记得在Matlab中运行gurobi_setup然后可以用yalmiptest验证是否配置成功。我的经验是首选GUROBI其次CPLEXintlinprog只作为兜底方案。2.3 代码结构设计我习惯按照数据-参数-模型-求解-绘图五个模块来组织代码。具体文件结构如下MPS_dynamic_scheduling/ ├── main.m // 主程序入口 ├── case33.m // 电网数据定义 ├── config.m // 参数配置 ├── build_model.m // 用YALMIP构建优化模型 ├── solve_model.m // 调用求解器求解 ├── plot_results.m // 结果可视化 ├── utils/ │ ├── distflow.m // DistFlow约束构建 │ ├── mps_mobility.m // MPS移动约束 │ └── soc_update.m // SOC动态约束主程序的大致流程是%% 主程序 clear; clc; close all; addpath(utils); run config.m; % 加载参数 mpc case33; % 加载电网系统 model build_model(mpc, params); % 构建模型 result solve_model(model); % 求解 plot_results(mpc, params, result); % 绘图这种模块化设计的好处是当你想修改目标函数中的权重系数或者加入新的约束时不需要牵一发动全身。比如我在后期加入了关键负荷优先恢复的权重参数只需要在config.m里定义权重向量然后在build_model.m中修改一行代码即可。3. MPS动态调度实操关键约束的代码实现3.1 MPS时空网络状态变量设计MPS动态调度的核心是把位置和电量两个状态关联起来。我定义两个二进制变量x(m,n,t)MPS m在时段t接入节点n取1或不在该节点取0。move(m,n1,n2,t)MPS m在时段t从节点n1移动到节点n2。同时定义连续变量p_dch(m,n,t)MPS m在时段t于节点n的放电功率。soc(m,t)MPS m在时段t结束时的荷电状态。YALMIP定义如下x binvar(M, N, T); % 位置指示 move binvar(M, N, N, T); % 移动指示可能维度很大需注意内存 p_dch sdpvar(M, N, T); % 放电功率 soc sdpvar(M, T1); % SOCt0为初始值这里要提醒一点move变量是四维的当N33时单个MPS就有3333T个二进制变量内存和求解负担非常重。我的优化方式是把移动变量压缩成move(m,path_idx,t)其中path_idx只列出实际可达的节点对剔除同一节点和不可达的组合。这一步能将问题规模直接砍掉一半以上。3.2 位置唯一性与移动约束每个MPS在任一时刻必须且只能存在于一个地方——要么接入某个节点要么在移动的路上。这个约束用数学表达就是% 每个MPS每个时段最多接入一个节点 for m 1:M for t 1:T Model [Model, sum(x(m,:,t)) 1]; end end移动约束是建模中最容易出错的部分。我们首先要保证如果MPS在t时刻接入节点n1且t1时刻接入节点n2n1≠n2那么它必须在t时段内完成从n1到n2的移动且移动时间不能超过一个时段长度。简化起见可以设定移动时间不超过1个时段这样只需要约束% 移动起止约束 for m 1:M for t 1:T-1 for n1 1:N for n2 1:N if n1 ~ n2 % 从n1到n2的移动需要x(m,n1,t)1且x(m,n2,t1)1 Model [Model, move(m,n1,n2,t) x(m,n1,t) x(m,n2,t1) - 1]; % 反向移动则禁止防止来回跑 Model [Model, move(m,n2,n1,t) 1 - x(m,n1,t) x(m,n2,t1)]; end end end end end上面的线性化技巧是把与关系转化为不等式这是MILP建模的经典手法。如果你是新手可能觉得这些约束有些绕但记住一个原则所有逻辑关系都要表达成线性不等式不能直接写x(1) x(2)这样的逻辑运算否则YALMIP会把模型变成非凸问题求解器根本处理不了。3.3 SOC动态约束与功率限制SOC的变化相对直观放电时SOC减少如果允许接入充电桩则可以增加。我这里假设MPS不能从电网取电只作为临时电源所以SOC更新方程是% SOC动态约束 (t1:T) for m 1:M Model [Model, soc(m,t1) soc(m,t) - sum(p_dch(m,:,t)*delta_t) / cap_m]; end这里cap_m是MPS容量delta_t是时段时长小时。注意离散化时如果时段较长功率乘以时间才是能量。很多人在这一步直接把功率加到SOC里导致单位错乱结果完全不可用。我建议在建模时统一单位功率用kW时间用h容量用kWh这样SOC变化量就是纯数值。功率限制必须与接入位置绑定% 如果x(m,n,t)0则p_dch(m,n,t)必须为0 for m 1:M for t 1:T Model [Model, p_dch(m,:,t) repmat(P_max, 1, N) .* x(m,:,t)]; end end这个约束用大M法实际上这里P_max就是M强制放电功率与接入位置一致。同样SOC上下限也要约束Model [Model, soc_min soc(m,:) soc_max];3.4 配电网潮流约束DistFlow模型配电网潮流计算有多种方法在调度优化中最常用的是DistFlow分支潮流模型。对于辐射状网络DistFlow可以写成节点电压降V_j^2 V_i^2 - 2(R_ij P_ij X_ij Q_ij) (R_ij^2 X_ij^2) * (P_ij^2 Q_ij^2) / V_i^2这个方程是非线性的。在MILP建模中我们通常忽略第二项并使用线性化近似的DistFlow或者用二阶锥规划SOCP松弛。但既然我们选择了MILP那就得用线性DisfFlow。我的做法是假设电压偏差不大忽略网络损耗得到% 线性化DistFlow电压降约束 for ij 1:L f branch(ij); t branch(ij); Model [Model, V(f) - V(t) 2*(R(ij)*Pij(ij,t) X(ij)*Qij(ij,t))]; end这里V是电压平方变量Pij/Qij是线路有功/无功功率。此外节点功率平衡方程需要把MPS注入功率加入进去% 节点功率平衡有功 for n 1:N Model [Model, sum(Pij(in_edges, t)) - sum(Pij(out_edges, t)) ... Pd(n,t) - sum(p_dch(m,n,t))]; end注意这里Pd是负荷有功如果MPS接入则相当于减少从电网侧汲取的有功。如果负荷可以切负荷那么还要加入zeta变量表示未恢复的负荷比例。我这篇复现里做了简化假设负荷只能全部恢复或全部不恢复用二进制recover(n,t)表示是否恢复节点n在时段t的负荷。这样目标函数就可以写成最大化加权恢复负荷。3.5 求解设置与结果展示构建完All Model后设置求解器参数很重要。我的经验是对于中等规模T24N33M2GUROBI的MIP Gap默认是1e-4如果直接求解可能要几分钟到几十分钟视计算机而定。为了提高速度可以设置ops sdpsettings(solver,gurobi,verbose,2); ops.gurobi.MIPGap 0.01; % 允许1%的gap ops.gurobi.TimeLimit 300; % 最多求解5分钟 ops.gurobi.MIPFocus 1; % 偏向快速找到可行解这样设置后通常几十秒内就能得到一个较优解。绘图方面我画出三个图1每个时段的恢复负荷比例曲线2MPS在空间上的移动轨迹用地理坐标画线3MPS SOC随时间的变化曲线。结果图上能很直观地看到MPS是如何先救援重要负荷再逐步向边缘移动SOC曲线则呈阶梯状下降恢复负荷比例逐渐上升。4. 复现中的常见问题与调试经验4.1 约束被错误松弛小心隐形的移动耗电我第一次运行模型时发现结果中MPS一个时段内能从电网的一端瞬移到另一端而且SOC不降。后来排查发现是因为我忘了在SOC更新方程中加入移动耗电项。在真实场景中MPS移动是需要消耗电量的通常可以设为每公里消耗的SOC百分比。加了这个约束后移动行为变得谨慎多了这也更符合论文的假设。建议在SOC更新方程中增加一项move_consumption(m,t)Model [Model, soc(m,t1) soc(m,t) - sum(p_dch(m,:,t))*delta_t/cap_m - move_energy(m,t)];其中move_energy是对所有移动距离的求和乘系数% 移动耗能约束 for m 1:M for t 1:T-1 Model [Model, move_energy(m,t) sum(sum(move(m,:,:,t .* dist(n1,n2)))) * k_consume / cap_m]; end end注意这里的dist(n1,n2)是节点间的距离矩阵k_consume是单位距离耗能比例。4.2 求解时间爆炸如何有效降维MPS动态调度最大的痛点就是规模。特别是移动变量move是四维的当M2, N33, T24时二进制变量数量约为23332*23≈4.8万加上位置和恢复变量总二进制变量轻松超过5万。如果初始模型不加任何削减GUROBI可能几个小时都找不到最优解。我的解决办法有几点限制最大移动距离实际中MPS不可能在短时间内跨越半个电网所以只允许满足dist(n1,n2) v_max * delta_t的移动候选其他变量直接不创建。时段聚合灾害初期的抢修阶段可以先用较粗的时段比如T12跑通后再细化。一开始我直接用T96模型爆炸到求解器直接Out of Memory改成T24后问题迎刃而解。对称性消除如果两台MPS是同质的求解器会把它们视为可互换对象导致大量对称分支。我通过给MPS增加初始接入节点编号小的优先的约束来破坏对称性。4.3 结果出现节点功率倒流问题某次运行时我发现恢复负荷比例虽然有上升但某些线路上的功率方向变成从MPS接入点向外送而MPS电量和容量明明不足以支撑。查了一圈问题出在潮流的功率平衡约束上我用的简化DistFlow没有考虑线路上两个方向都可能有功率的情况导致在支路功率变量上允许了正负同时存在即变量无界优化器就利用这个漏洞凭空创造功率。解决方法是给支路功率变量加上合理的上下界或者把支路功率拆分成正向和反向两个非负变量。4.4 常见问题速查表问题现象可能原因解决方案MPS瞬移移动时间约束未生效检查move变量约束是否需要x(t)x(t1)-1SOC变化单位不对功率*时间与容量单位不一致统一为kW、h、kWh求解器计算缓慢二进制变量过多削减移动候选、增大MIPGap、破坏对称性负荷恢复比例异常节点功率平衡中负荷方向错误检查Pd的正负号约定结果对Gap敏感目标函数数值尺度差异大对目标各项系数做归一化这些坑我想只要复现过类似问题的人大概率都遇到过。如果你准备上手建议先跑一个N5节点的简化网络把所有约束调试正确再替换成33节点。这样能极大减少排错地狱的时间。5. 复现心得与扩展建议5.1 如何让复现结果更贴近论文论文里的结果通常是在特定算例下得到的你想完全复现一模一样的数据几乎不可能除非作者公开了全部代码和参数。所以我的目标不是复现出同样的数字而是复现出同样的规律和曲线形态。如果你发现论文里恢复负荷比例是98%你的结果是95%先别慌看看是不是MPS容量、速度、负荷权重的参数设定有差异。特别是负荷权重论文里可能用关键负荷等级加权而你只用总计恢复电量这会导致调度策略有明显区别。我的建议是先按照论文的假设重新推导一遍模型再把每个参数的物理意义对照清楚最后才能放心修改成自己的场景。另外很多论文的潮流约束使用的是二阶锥规划SOCP我用的线性DistFlow其实是简化版所以结果有5%以内的偏差是正常的。如果你追求精度可以改用YALMIP的optimize(..., ops)直接求解SOCP但那样就不属于严格的MILP了求解时间也会变长。5.2 扩展多类型移动电源与协同调度这个项目基础打牢后你可以很自然地扩展到以下方向多类型MPS比如同时有大型储能车和小型移动发电机它们的容量、速度、成本不同模型里只需增加MPS类型索引即可。与固定储能协同把MPS和固定储能ESS放在同一框架下共享节点功率平衡约束但ESS不能移动约束要少很多。考虑交通路网约束MPS移动时间不是简单距离除以速度还要考虑道路拥堵和损坏。可以建立一个交通网络模型用最短路径时间矩阵代替简单距离矩阵。我自己正在尝试的扩展是把动态调度与网络重构开关操作联合优化。传统的配电网重构是控制联络开关和分段开关来改变拓扑从而 transfer 负荷。如果把MPS移动和开关重构放在一起问题的可行域会变得非常复杂但恢复能力也会大幅提升。不过这个模型的求解难度系数呈指数增长我已经做好和求解器长期斗争的准备。5.3 给初学者的三条建议第一不要直接啃大代码。我见过很多同学下来一个几百行的复现包打开后一脸蒙。正确的打开方式是自己先写一个小模型比如两节点、一台MPS、三个时段把每个约束都亲手写出来理解它们是如何把物理规则翻译成数学不等式的。这个过程虽然慢但能让你在调试大模型时迅速定位问题。第二学会用YALMIP的debug工具。当模型不可行时YALMIP会提示infeasible problem但不会告诉你是哪条约束出了问题。这时可以借助optimize返回的diagnostics信息结合check命令检查每个约束的残差逐步注释可疑约束二分法定位问题源。我至少有一半的建模错误是靠这个方式找出来的。第三注意Matlab版本和工具箱兼容性。YALMIP目前对R2022b以上的版本支持很好但GUROBI插件有时因为Java版本不匹配而加载失败。建议在开始前仔细阅读官方setup文档并测试一个简单LP问题比如min x约束x1确保求解器调用链路通畅。最后分享一个小技巧在跑大规模算例时可以先用较小的MIPGap比如0.1快速得到一个可行解把这个解作为初始解传给后续细化求解通过assign和optimize的x0参数能够明显加快收敛。这是我多次实验后发现的作弊通道在复现大论文时特别管用。希望这些踩坑经验和代码思路能让你在MPS动态调度这条路上少走几步弯路。
