这段时间一直在折腾一篇电力系统方向SCI一区论文的复现工作——基于配电网韧性提升的应急移动电源预配置与动态调度。论文的核心思路并不复杂但真正落地到Matlab代码时很多细节远超我的预期。尤其在处理MPSMobile Power Source移动电源预配置这个子问题时我踩了不少坑也积累了不少经验。趁着休息把上篇整理出来重点讲清楚预配置阶段的模型构建、代码实现和避坑手法给同样在复现这类文章的朋友做一个参考。先说下这篇论文解决什么问题台风、冰灾这类极端天气来袭前配电网需要提前把应急移动电源布置到关键位置灾后再根据实际故障情况调整这些电源去恢复失电负荷。预配置就是在灾害发生前回答应急电源应该提前放在哪、放几台、容量多大这三个问题。这篇博文是我整个复现过程的上半部分核心聚焦在预配置阶段的Matlab实现动态调度的下半部分后面会单独整理。无论你是在做韧性评估、孤岛划分还是应急电源优化调度这套思路和代码框架都值得你花点时间看完尤其是那些刚接触YALMIP配Cplex、正在纠结怎么把论文公式变成可运行代码的朋友。1. 项目概述为什么应急电源要提前部署而不是事后调动1.1 配电网韧性和应急移动电源的基本概念配电网韧性Resilience这个词在2012年桑迪飓风之后被电力行业重点讨论。学术界一般把它定义成配电网对极端事件的抵御、吸收、适应和恢复能力核心关注的是大扰动、低概率、高影响的事件。和传统的可靠性Reliability不同可靠性研究的是N-1这种常规故障关注停电频率和时长韧性面对的是极端灾害下可能同时出现的大量故障关注的是系统能不能快速恢复。应急移动电源就是用来提升韧性的重要手段主要包括移动储能车、应急发电车、移动UPS这类可以灵活转场的电源装备。它们的共同特点是平时在指定站点待命灾害来临时可以移动到需要供电的节点和被动的固定应急电源相比最大优势是灵活。但这份灵活带来一个关键难题——灾害发生时道路可能受损移动时间受限如果等故障出现了再临时调度往往来不及。所以论文提出要在灾害预警阶段就把移动电源预配置到风险较高的区域附近这就是预配置pre-positioning的动机。提示MPS这类设备在英文文献里也叫MESSMobile Energy Storage System或者MEGMobile Emergency Generator搜索相关论文时这几个关键词都可以试不同文章的侧重点不一样MESS偏储能约束MEG偏出力约束。1.2 预配置阶段决策什么预配置阶段发生在灾害预警期还不确定具体哪个节点会故障只能基于气象预报和历史数据估计每个节点的故障概率。这个阶段要做的决策包括三块确定预配置站点从候选节点里选出哪些节点作为移动电源的初始接入点确定预配置数量在选出的站点投入多少台MPS确定预配置容量每台MPS带多少电能/燃料带多了浪费带少了恢复效果不够。这三块决策的背后其实是成本约束下的风险对冲问题。论文用数学语言表达为在预算限制下选择一个预配置方案使得后续所有可能故障场景下的负荷损失期望值最小。仅仅看到这里你应该已经意识到这个问题的本质是一个两阶段随机优化第一阶段在做决策时不知道灾害后果第二阶段等故障场景清楚了再调整移动电源的位置和出力。MATLAB的强大之处在于它提供了优化工具箱和YALMIP这样的高级建模框架让我可以用接近论文公式的形式写出数学模型而不是自己去手写求解算法。1.3 我的复现路线图整个复现路线我拆成了四步每步对应一个核心模块数据准备模块构建IEEE 33节点配电网参数生成候选场景集模型构建模块根据论文公式写出目标函数和约束条件采用YALMIP建模求解配置模块调用Cplex/Gurobi求解MILP提取最优预配置方案结果可视化模块绘制系统拓扑图和预配置效果对比图直观展示方案有效性。这四个模块中最耗时间的不是写代码而是理解论文的模型假设和变量定义。很多公式在附录里是缩写的如果不搞清楚上标、下标的确切含义代码写到一半就会卡壳。2. 论文模型解构预配置和动态调度是怎么打包进一个模型的2.1 两阶段随机规划框架复现的第一步是把论文框架看懂。我选的这篇论文走的是两阶段随机规划路线这也是目前电力系统不确定性决策用的最多的框架。第一阶段是预配置决策在灾害预报阶段完成第二阶段是动态调度决策在故障信息明确后执行。两个阶段通过预配置的位置这个变量耦合在一起。用公式框架表达就是min 第一阶段的预配置成本 E[第二阶段的负荷损失] s.t. 第一阶段预算约束 第二阶段所有场景下的运行约束这里E表示对所有故障场景取期望。用YALMIP写这种模型时核心思路就是把随机变量场景显式写出来把第二阶段问题对每个场景分别建模再用期望算子汇总到目标函数。要注意的是很多论文里第一阶段和第二阶段共享同一个MPS空间位置变量——预配置决定了移动电源的初始位置动态调度时可以从这个初始位置出发也可以部分改变位置。这意味着第一阶段决策对第二阶段的影响不是硬性的必须留在原地而是从某个初始点出发。代码实现时需要把移动电源的可移动性建模为初始接入点固定后续节点可通过路径变量切换的形式。2.2 目标函数不只是最小化负荷损失我复现的论文中目标函数由三部分构成预配置成本包括移动电源租赁成本、部署到站点的运输成本、充电/加油成本预期负荷损失所有场景下失负荷量的加权期望权重代表不同等级负荷的重要程度惩罚项对预配置方案违反约束的罚函数在部分论文中会通过权重系数并入目标。这里最关键的设计是负荷权重的设定。如果所有负荷一视同仁优化结果往往会优先恢复距离MPS初始位置近的小负荷而忽略重要的医院、数据中心等关键负荷。所以模型中每个节点都有重要度系数ω_i目标函数里用的是ω_i乘以失负荷量而不是单纯的失电量。这一点在复现时容易被忽略但恰恰是论文能发SCI一区的关键贡献之一。2.3 关键约束条件解读约束条件中处理起来最棘手的是这几条移动时间约束MPS从当前位置移动到目标接入点的时间不能超过允许的最长移动时间。这部分在论文里通常用分段线性函数或者交通模型来表达但复现时一般简化为移动时间与距离成正比即t_ij d_ij / v。还需要加一个0-1变量表示是否从i移动到j。潮流约束配电网要满足DistFlow方程即节点功率平衡、支路潮流和电压降落关系。这是非线性方程需要线性化处理。常见的做法是采用线性DistFlow模型LDF假设电压幅值近似为1pu忽略高阶小项。实际复现中如果系统是33节点这种规模LDF的精度完全够用而且能大幅降低求解难度。MPS容量约束每台MPS的容量有限提供的功率不能超过额定功率累计放电量不能超过额定容量。如果MPS是储能车而不是发电车还需要考虑充放电效率以及每个时段的SOC状态转移方程。辐射状拓扑约束配电网运行必须满足辐射状拓扑不能形成环网。在预配置阶段主要考虑的是MPS接入后的孤岛划分。如果MPS接入到某个节点该节点下游的负荷可以被孤岛供电。实现时可以用单商品流模型或生成树约束来保证拓扑是辐射状但这类约束往往会让求解时间急剧上升。2.4 不确定性建模场景生成与缩减两阶段随机规划的前提是场景集。我复现时用了蒙特卡洛采样生成故障场景每个场景记录了哪些线路故障、哪些节点失电。直接蒙特卡洛采几百个场景扔进模型求解时间会爆炸所以必须做场景缩减。目前最常用的手段是K-means聚类K均值中心选取% 场景缩减K-means聚类 % scenarios为n_scen*n_line的0-1矩阵表示每条线路是否故障 [idx, C] kmeans(scenarios, K, Distance, hamming); reduced_scenarios C; % K个聚类中心作为代表性场景这里一个比较容易被忽视的细节是聚类中心虽然是0-1矩阵但K-means得到的中心可能是小数。所以场景缩减后需要把中心重新二值化按0.5阈值取整否则在第二阶段模型里会出现线路半故障这种物理上不存在的场景。注意场景缩减不是随便删几个场景就行。要检查缩减前后各节点的故障概率是否基本一致否则优化结果会偏移。每次缩减后我都对比原始场景和缩减场景的节点失电概率分布误差超过5%就要重新选K值或换距离度量。3. Matlab代码实现从论文公式到可运行代码3.1 代码整体架构用OOP理清复杂关系复现初期我把所有代码塞进一个脚本结果改一个参数就要跑半天调试极其痛苦。后来下决心重构了代码结构参考了铜陵学院李光耀那套基于Matlab OOP架构的多算法融合系统设计思路把整个项目拆成了三个层次数据层存放电网参数、场景数据、MPS参数的类主要负责数据读写和预处理模型层构建优化模型的类负责定义变量、目标函数和约束求解层调用求解器处理求解结果并输出。这套OOP架构的好处是更换测试系统时只需要改数据层的对象更换论文目标函数时只需要改模型层的目标函数方法更换求解器时只需要改求解层的配置。对需要反复实验的复现工作来说这种模块化改造的投入产出比非常高。我在代码里定义了以下几个核心类classdef MPS properties id % 移动电源ID type % 类型储能/发电 rated_power % 额定功率 (kW) rated_energy% 额定容量 (kWh) initial_node% 预配置初始位置 speed % 平均移动速度 (km/h) cost % 租赁/使用成本 end end3.2 YALMIP建模代码精讲这是整个复现最核心的部分。要明确的是直接调用优化工具箱里的intlinprog是可行的但如果有多阶段、多周期的约束嵌套代码量会非常庞大而且可读性差。所以我的选择是YALMIP——Matlab下的免费优化建模工具箱它最大的价值是把数学表达式和求解器解耦让你只关心模型本身。先看变量定义部分% 决策变量 % y(i) 二进制变量节点i是否被选为预配置站点 % n_mps(i)整数变量节点i预配置的MPS数量 % p_out(s,i,t)连续变量场景s下节点i在t时刻的实际出力 % l_shed(s,i,t)连续变量场景s下节点i在t时刻的失负荷量 % move(s,i,j)二进制变量场景s下MPS是否从i移动到j y binvar(n_node, 1); n_mps intvar(n_node, 1); p_out sdpvar(n_scen, n_node, n_time); l_shed sdpvar(n_scen, n_node, n_time); move binvar(n_scen, n_node, n_node);这里一个关键点预配置决策y和n_mps是不带场景索引的因为它们在第一阶段就已经确定不随场景变化。而p_out、l_shed、move这些第二阶段变量必须带场景索引。不少复现者在这里犯错误——把y也定义了场景维度结果模型变的巨大无比求解速度极慢而且逻辑上也不对第一阶段决策不能依赖未来场景信息。目标函数写成objective sum(y .* c_pre) sum(n_mps .* c_rent) ... (1/n_scen) * sum(sum(sum(w .* l_shed, 3), 2), 1);对应的约束条件以容量约束和移动约束为例% 容量约束MPS在任意时刻出力不能超过额定功率 for s 1:n_scen for i 1:n_node Constraints [Constraints, ... p_out(s,i,:) n_mps(i) * P_rated * y(i)]; end end % 移动距离约束简化版 for s 1:n_scen for i 1:n_node for j 1:n_node Constraints [Constraints, ... move(s,i,j) y(i) * available_path(i,j)]; end end end3.3 求解器选择与配置YALMIP只是建模语言真正求解还要靠底层的MILP求解器。我测试了三款常见求解器的表现求解器许可证适用规模备注Gurobi商业有学术许可中大规模求解速度最快能处理二次约束我最终选它Cplex商业有学术许可中大规模和Gurobi不相上下部分版本在纯MILP上更稳intlinprogMatlab自带中小规模免费无需额外安装但30节点以上多场景时速度明显下降我的测试配置是33节点、20个缩减场景、24个时段。intlinprog跑了将近40分钟还没收敛换成Gurobi之后大约4分钟就出结果了。如果你没有商业求解器的学术许可建议先用intlinprog验证小规模模型正确性再考虑申请Gurobi的学术授权基本秒批。调用求解器的代码很简单ops sdpsettings(solver, gurobi, verbose, 2, ... debug, 1, showprogress, 1); sol optimize(Constraints, objective, ops);这里debug1参数值得单独说一下。YALMIP的debug模式能在模型不可行时给出具体是哪条约束出了问题复现阶段80%的时间是在排错中度过的这参数简直是救命的。4. 复现实操全流程我就是这么一步步走通的4.1 IEEE 33节点系统数据准备绝大多数配电网韧性论文都在IEEE 33节点系统上验证节点和线路参数在Matpower里可以直接调出来也可以网上搜标准数据。需要注意的是论文里可能不会明确给出所有参数部分数据需要参考原文的算例部分自己推演。我的做法是把每个节点的有功/无功负荷、线路阻抗、开关状态整理成结构体统一存到电网数据类里。bus_data [ 1, 0, 0, 0, 0; 2, 100, 60, 0, 0; 3, 90, 40, 0, 0; ... ]; % 列节点编号有功负荷(kW)无功负荷(kvar)注入有功注入无功这里一个容易坑人的细节是IEEE标准系统的负荷是恒功率负荷但在极端事件下的负荷模型可能需要考虑电压依赖性即实际消耗功率随电压降低而减少。大部分论文简化为恒功率模型但如果你复现的论文用了ZIP负荷模型那在代码里要把电压项引入潮流方程复杂度会高一个量级。我建议优先复现恒功率版本验证了框架后再扩展。4.2 场景生成蒙特卡洛采样故障场景按论文假设生成每条线路有独立的故障概率p_fault故障事件之间相互独立。这个假设在真实灾害中不完全成立台风往往沿着路径损坏线路但在论文复现时按原文假设即可。n_scen_raw 5000; % 原始采样场景数 n_scen_red 20; % 缩减后场景数 scen_raw rand(n_scen_raw, n_line) p_fault; % 过滤掉全0场景无故障场景 scen_raw scen_raw(sum(scen_raw, 2) 0, :); % K-means聚类缩减 [~, Centers] kmeans(scen_raw, n_scen_red, Distance, hamming); scen_reduced double(Centers 0.5);场景数选择上有一个权衡太少则随机性刻画不足方案过于保守或过于激进太多则模型规模指数爆炸。我实测下来33节点系统用20~30个场景是性价比最高的区间。如果你不想用K-means也可以用快速前向选择算法Fast Forward Selection效果类似但实现起来复杂一些。4.3 求解与结果输出一个完整的预配置方案跑通模型后程序输出预配置结果。一个典型的输出长这样 MPS预配置方案IEEE 33节点系统 预配置站点 MPS数量 MPS类型 总容量(kWh) 节点 12 2 储能车 1200 节点 25 1 发电车 800 节点 30 1 储能车 600 期望负荷损失优化前856.3 kW 期望负荷损失优化后412.7 kW 韧性提升率51.8% 这个方案在20个场景下的平均表现是失负荷量下降超过一半。当然方案是否最优取决于预算上限和MPS总量约束需要做敏感性分析来验证总预算提高20%会不会带来等比例的韧性提升这个分析对论文审稿人来说很关键对复现者来说也能检验模型逻辑是否正确。数据可视化上Matlab自带的绘图功能完全够用。我用plot函数配合自定义的节点坐标画出33节点系统拓扑图用红色圆点标记预配置站点用绿色区域标记孤岛恢复范围再用柱状图对比优化前后各场景的失负荷量。这套图做出来之后对理解模型行为帮助巨大。有几次我以为是代码bug结果看图才发现是场景数据本身有偏可视化在调试中能帮你快速定位问题别偷懒省掉这一步。4.4 灵敏度分析验证模型的合理性复现完主模型后强烈建议做一组参数灵敏度实验。我自己做的几组实验中最有价值的是MPS总量变化从3台增加到10台观察韧性提升曲线的变化。这条曲线会呈现边际递减趋势说明在某个数量之后增加MPS对韧性的提升不再显著这个拐点对实际应急资源配置有指导意义。移动速度变化把MPS平均移动速度从20km/h提高到60km/h观察预配置方案是否发生变化。直觉上速度越快预配置站点应该越分散因为MPS有更强的机动性去覆盖更多区域。预算约束变化收紧预算时模型应该自动放弃负荷重要性低的节点优先保重要负荷。如果这些定性结论在代码里体现不出来就要检查模型是不是有逻辑问题了。我第一个版本跑出来的结果是增加MPS数量韧性反而下降最后查出来是预配置的固定成本在目标函数里权重设置不合理导致模型宁愿少配MPS也不愿意承担高额的预配置成本——这类问题光靠看代码很难发现必须靠灵敏度实验来暴露。5. 常见问题与排查技巧实录5.1 YALMIP报错No suitable solver found这个错误90%是因为YALMIP没有找到可用的MILP求解器。检查两步第一步看Matlab路径里有没有添加YALMIP和求解器的路径第二步用yalmiptest命令测试哪些求解器可用。yalmiptest % 检查solver列表里是否有gurobi或cplex如果yalmiptest里gurobi显示Fail大概率是Gurobi的安装路径没有正确添加到Matlab搜索路径或者是Gurobi版本和你的Matlab版本比如2026b不兼容。这时候可以安装写addpath命令把Gurobi的matlab接口目录加进来。5.2 模型不可行Infeasible problem这是复现阶段最让人崩溃的问题。YALMIP的debug1参数会告诉你是哪条约束导致了不可行。我遇到的情况分成三类移动约束过强可移动路径矩阵把某些合法路径排除了导致MPS无法从预配置点到达目标节点。检查available_path矩阵的非零元素是否连通。容量约束冲突n_mps(i)定义为整数且下限为1时如果预配置站点y(i)0约束p_out n_mps(i) * P_rated * y(i)自动变成p_out 0也就是说MPS容量被错误清零了。正确逻辑应该是y(i)0时n_mps(i)也必须为0这需要利用Big-M约束绑定两个变量。功率平衡约束不满足场景生成时如果多个相邻线路同时故障有些节点变成孤岛且没有MPS接入功率平衡就无法满足。这需要在模型中允许失负荷变量l_shed非零而不是强制完全平衡。5.3 求解时间过长求解时间爆炸一般有两个原因原因一Big-M值过大。Big-M取的越大MILP的松弛越松求解器需要探索的分支越多。我刚复现时给移动时间约束取了M10000结果求解器卡了几个小时。后来把M缩小到最大移动时间的1.5倍求解时间从几小时降到几分钟。这不是什么花哨的技巧但非常有效。原因二对称性问题。两台完全相同的MPS模型分不清谁是谁导致求解器在等价分支之间浪费时间。解决办法是给MPS加编号约束比如要求MPS的预配置节点按编号升序排列n_mps(i)和n_mps(j)之间增加排序约束。这种方法能显著减少求解空间。5.4 不可忽视的数值稳定问题Matlab默认双精度足够支撑大部分优化问题但在MILP中二进制变量和连续变量刚量级差异巨大时约束矩阵的条件数会很大求解器可能出现数值困难。我的建议是所有参数统一量纲。有功功率用kW而不是MW距离用km而不是m别混用。对成本参数做归一化处理。把预配置成本除以总预算让目标函数的各项系数保持在0.1~10这个量级能明显改善求解稳定性。这个经验是我被一个hidden bug折磨了三天后总结出来的原来那篇论文的原文里电费损失的单位是MWh而我换算时多乘了一个1000导致目标函数中负荷损失项的权重异常大优化结果完全偏向某几个节点方案失去参考意义。6. 下篇预告从预配置到动态调度的衔接预配置做完下半场就是MPS的动态调度了。动态调度和预配置最大的不同在于变量从静态选址变成了动态路径规划每个时段都要决策MPS往哪走、给谁供电。这会引入时间维度上的状态转移约束模型规模比预配置大一个量级。目前我整理的思路是这样的动态调度阶段会把故障场景固定成确定性事件即假设灾后信息已知然后建立多时段混合整数线性规划模型目标函数是最小化整个恢复周期内的综合损失。MPS的路径用时空网络来描述每个时段MPS要么在某个节点供电要么在移动。这里和预配置阶段最重要衔接点是动态调度的初始位置就是预配置方案的输出。所以我的代码架构里预配置类的结果直接作为动态调度类的输入参数两个阶段通过MPS对象的initial_node属性天然衔接。如果你也打算复现这类论文一个建议是不要急着写动态调度的代码先把预配置阶段的每个约束吃透特别是把潮流线性化方法和Big-M约束的选取逻辑搞明白。这些基础打牢了动态调度阶段虽然变量多但代码风格和建模模式是一样的无非是多套一层时间循环。目前这套代码在两个项目里被复用过了一个是沿海城市的台风灾害应急预演另一个是配电网灾后孤岛划分策略验证。后续我还会在现有模型基础上加入交通网络约束和移动充电时间窗把预配置和动态调度的耦合做得更细一些等代码重新梳理完后再来更新。复现SCI论文不是一个轻松活尤其从零开始的时候。但这个过程对理解和掌握优化建模技术的帮助是巨大的——你不仅学会了YALMIP怎么写约束更学会了怎么把一个复杂的物理问题转换成数学问题再转换成代码。希望这篇上篇能帮你少踩几个坑下篇我们接着聊MPS怎么在路上跑起来。