刚看到这个题目时，我以为是“应急电源车选址”的简单变种，真正把MPS预配置的模型读进去、再用Matlab逐行复现出来，才发现这里面的门道比想象中深很多。它并不回答“某条线路坏了怎么救”，而是在极端灾害还没发生之前，就把移动电源放到最合理的节点、配好容量，让灾后的动态调度不至于“巧妇难为无米之炊”。这篇复现文章前后陆陆续续花了我三周，中间踩了不少坑，也把Yalmip建模、二阶锥松弛、求解器调试的细节摸了一遍。下面就把从模型拆解到代码落地的完整过程整理出来，给正在做配电网韧性、应急电源优化配置的同行一点参考。
这个模型非常适合正在复现电力系统优化类论文的研究生、刚接触配电网韧性计算的分析人员，也适合想用Matlab把混合整数二阶锥规划落到实际工程场景的开发工程师。无论你是想把论文公式变成可运行代码，还是想弄明白“预配置”和“动态调度”在数学上到底怎么衔接，这篇梳理都能帮你少走弯路。
1. 动手前我先想清楚的事情：这个预配置到底在配置什么
1.1 为什么预配置看起来简单，跑起来却很烧脑
应急移动电源系统（Mobile Power System，简称MPS）在配电网韧性研究里，本质上是“移动式的备用电源”，可以是移动储能车，也可以是车载柴油发电机。预配置解决的问题是：灾害还没发生、破坏情况未知的时候，手头有限数量的MPS应该放到哪些节点、每个节点配置多大容量，才能让灾害发生后的负荷恢复效果最好。
乍一听这就是一个选址加容量的组合优化问题，但它真正麻烦的点在于“不知道故障在哪”。极端天气下的配电网故障通常有成百上千种可能场景，比如台风路径不同、倒塔位置不同、线路断线组合不同。如果等到故障发生后再临时调度MPS，从出发点到目标节点的移动时间可能长达几个小时，关键负荷早就停电了。预配置的价值，正是把“资源提前摆到最有价值的位置”，用事前准备对冲事后不确定性。
建模的时候，这个“不确定性”会直接抬高问题的复杂度。如果和灾后动态调度一起建模求解，决策变量会同时包含预配置选址、灾后MPS移动路径、输出功率、配网潮流状态，再加上多场景多时段，问题规模会爆炸。这篇工作的做法是把整个决策链条拆成两阶段：预配置阶段只管“放哪、放多大”，动态调度阶段再管“怎么移动、怎么用”。两阶段之间通过资源约束衔接，既保证了问题可解，也符合应急管理的实际操作顺序。
动手复现之前，我花了一个晚上把原文的目标函数、约束条件和假设梳理成表格。强烈建议你也这样做：把每个符号代表什么、哪些是第一阶段决策、哪些是第二阶段决策、哪些是公共参数列出来，再开始写代码。跳过这一步直接抄公式，后期调试时大概率会一头雾水。
1.2 我选的复现技术路线：MISOCP加Yalmip
原问题严格来说是一个混合整数非线性规划，配电网潮流方程里面电压幅值和功率乘积天生就是二次项，直接扔给求解器很难处理。我采取的方案是把非线性潮流约束做二阶锥松弛，转换成混合整数二阶锥规划（MISOCP），然后交给Gurobi求解。
这里补充一个背景知识：MISOCP是混合整数线性规划的升级版，允许变量之间的二次关系以锥形式出现，适合处理带二次约束但又有离散决策变量的优化问题。配电网的DistFlow方程经过变量替换后，可以把非凸二次等式松弛成凸的二阶锥不等式，这个操作在理论上有严格条件保证最优解一致，在工程复现中也足够稳定。
Matlab这边用的是Yalmip工具箱。Yalmip的核心价值是让建模语法贴近数学表达式：sdpvar声明变量、约束写在方括号里、optimize一行调用外部求解器。相比手动把所有约束写成矩阵系数形式，Yalmip至少能省掉一半建模时间，而且修改目标函数、增加约束都很方便。
Gurobi做MISOCP的优势在于自带二阶锥割平面、对称性探测、场景预处理等技术，收敛速度和数值稳定性比Matlab自带求解器好太多。尤其是故障场景数量超过30个之后，开源求解器那种“跑一晚上还不收敛”的情况，在Gurobi上通常几十分钟内就能给出高质量解。
我还额外做了一件事：先写一个单时段、单一故障场景的简化版模型验证逻辑，确认变量定义、约束格式都没有问题，再往多时段、多场景扩展。这个小技巧帮我省了至少两天排查时间，后面第五部分会详细展开。
2. MPS预配置模型的数学拆解：目标函数、决策变量和约束
2.1 目标函数到底在优化什么
预配置的目标函数通常写成最小化“加权失负荷量”，这里的核心是“加权”两个字。配电网里不同节点的重要程度差别很大，医院、通信基站、应急指挥中心这些一级负荷的停电损失远高于普通居民负荷，所以目标函数要按节点重要程度分配权重系数。
数学模型可以写成：
[ \min \sum_{s\in S} p_s \sum_{t\in T} \sum_{i\in N} w_i \cdot L_i \cdot (1 - z_{i,t,s}) ]
其中S是故障场景集合，p_s是场景s的发生概率，T是调度时段集合，N是配电网节点集合，w_i是节点i的重要度权重，L_i是节点i的负荷值，z_{i,t,s}是0-1变量，表示场景s下时段t节点i是否恢复供电。这个目标函数的值越小，说明整个系统的预期负荷损失越小，韧性越高。
复现时有一个容易忽略的点：很多论文会在目标函数里加上MPS预配置的固定成本项，比如移动电源的部署费用、运输准备费用。加成本项会改变最优解的偏好，让模型倾向少配置电源。我建议先把目标函数里的“物理量”跑通，也就是只看失负荷最小化，确认结果符合直觉后再加入经济项做灵敏度分析。否则模型不收敛时，你很难判断是目标函数权重配比问题，还是约束本身有问题。
2.2 决策变量和约束体系
MPS预配置阶段的决策变量主要有三个维度：在哪里配置、配置多大容量、是否作为应急电源接入节点。具体来说：
- 预配置位置变量x_i：0-1变量，表示节点i是否作为MPS配置节点。
- 预配置容量变量C_i：连续变量，表示在节点i配置的MPS容量，受单台MPS容量上限约束。
- MPS接入状态变量u_{i,t,s}：0-1变量，表示故障场景s下时段t节点i是否有MPS可接入供电。
约束体系可以从四个层面理解。第一是节点功率平衡约束，这是配电网潮流的基础，MPS输出功率、配电网注入功率、负荷消耗功率必须在每个节点每个时段保持平衡。第二是线路容量约束，每条线路的传输功率不能超过热稳定极限，极端工况下允许部分过载但需要明确定义。第三是MPS资源约束，所有节点配置的容量总和不能超过可用MPS总量，单个节点容量不能超过额定值。第四是运行逻辑约束，比如某节点没有预配置MPS就不能在故障后接入输出，这组约束把第一阶段决策和第二阶段状态联系起来。
这四类约束缺一不可。我最早复现时漏掉了第四类关联约束，结果模型给出的方案是在一个没有预配置电源的节点上疯狂输出功率，目标函数值低得离谱，但物理上一眼就能看出不可能。这就是“阶段衔接约束”的意义。
2.3 非线性项的处理：二阶锥松弛怎么落下来
配电网潮流的精确约束是非线性的，直接建模会让问题变成难以精确求解的混合整数非线性规划。工程上最常见的做法是使用DistFlow分支潮流方程，通过变量代换把二次项转成线性项和二阶锥项。
变量代换分三步：把节点电压幅值的平方定义成新变量，等于每条支路潮流；电流幅值的平方也定义成新变量，等于功率与电压之比；节点功率平衡方程里的二次项会被这几组新变量替换，得到一组线性等式约束加上一个二阶锥不等式约束：
[ \left| \begin{array}{c} 2P_{ij} \ 2Q_{ij} \ V_i^2 - I_{ij}^2 \end{array} \right|2 \leq V_i^2 + I{ij}^2 ]
这个不等式在Yalmip里可以用cone函数直接表达，不需要手动展开，非常方便。把原来的非凸等式松弛成凸不等式后，求解难度大幅下降。理论上有论文证明在一定条件下松弛是精确的，也就是说松弛后的最优解等价于原问题的最优解，但在复现时仍要留意个别场景下松弛间隙过大的情况。我碰到过两次松弛误差达到百分之几的场景，排查下来是故障场景切负荷比例过高导致辐射状网络边界条件恶化，后面第5部分会说处理方法。
3. 用Matlab实现MPS预配置：数据准备、建模和求解
3.1 运行环境与工具箱版本
我使用的环境是Matlab R2021a配合Yalmip R20230630版本和Gurobi 9.5.2。这里特别提醒一点：Yalmip和Gurobi的兼容性并不是越新越好，早期版本Yalmip在某些情况下对Gurobi 10的接口支持不完整，建议用相对成熟的组合。Matlab自带的优化工具箱也能处理MISOCP，但性能上不如专用商业求解器，尤其当场景数超过20个时差距非常明显。
安装Yalmip后第一步先运行yalmiptest命令，确认所有求解器接口正常识别。Gurobi安装时要确保licence文件可以被当前用户读取，很多看似Yalmip报错的问题，实际是许可证路径没有配置好。
另外一个容易被忽略的环境问题是路径冲突。如果你电脑上装了多个版本的Yalmip或者多个求解器目录同时出现在Matlab路径里，调用optimize时很可能出现了不同求解器同时候选的情况，导致结果异常。干净的做法是新建一个专用工作区，把Yalmip、Gurobi对应目录按唯一优先级加入路径。别嫌麻烦，我曾在两个旧项目目录里残留了旧版Yalmip，结果同一个模型跑出来的最优值和论文对不上，花了两天才排查出是路径污染。
3.2 数据准备：以IEEE 33节点系统为例
复现配电类论文，最常用的标准算例是IEEE 33节点配电网。这个系统有33个节点、32条支路，基准电压12.66千伏，总负荷约3715千瓦加2300千乏，网络结构是典型辐射状。数据一般以两个列表形式给出：节点数据（编号、有功负荷、无功负荷）和支路数据（起端节点、末端节点、电阻、电抗、容量上限）。
读取数据时我建议直接在Matlab脚本里以矩阵形式硬编码，或者用readtable读Excel文件。硬编码的好处是算例数据一眼就能看到，改起来方便，适合复现阶段；如果后续要做批量测试，再改用统一的输入文件。对IEEE 33节点系统，每个节点的负荷和每条支路的阻抗参数必须仔细核对，网上流传的版本偶尔会把电阻和电抗单位搞混，导致潮流计算出的网损完全对不上。
另外，MPS的数量和容量参数也是关键输入。常用的假设是：可用MPS总数为3台，单台容量在200千瓦到500千瓦之间，移动时间忽略不计或按支路连通性简化处理。这些参数在原文中都有设定，复现时不要拍脑袋改参数，先按默认值跑通，再在灵敏度分析中调整。
3.3 Yalmip建模核心代码片段
下面这段代码展示了预配置模型最核心的部分：决策变量声明、目标函数和关键约束。这个片段不是完整代码，但已经把建模骨架讲清楚了。
%% 参数初始化 nb = 33; % 节点数 nl = 32; % 支路数 nt = 24; % 时段数 ns = 50; % 故障场景数 %% 决策变量 x = binvar(nb, 1, 'full'); % 预配置位置 e = sdpvar(nb, 1, 'full'); % 预配置容量 P_mps = sdpvar(nb, nt, ns, 'full'); % MPS有功输出 Q_mps = sdpvar(nb, nt, ns, 'full'); % MPS无功输出 z = binvar(nb, nt, ns, 'full'); % 负荷是否恢复 %% 目标函数 objective = 0; for s = 1:ns for t = 1:nt for i = 1:nb objective = objective + p(s) * w(i) * P_load(i) * (1 - z(i,t,s)); end end end在看到变量定义时候，有两点值得解释清楚。第一，P_mps声明为三维变量，维度是节点数乘时段数乘场景数。如果你的Matlab内存比较紧张，这种全排列方式在场景数超过100时会有一定负担，可以考虑把场景维度拆到循环里面做加和，但会让代码复杂不少。我优先用了全排列，逻辑直观、后续加约束方便。第二，z是负荷恢复变量，它和目标函数直接挂钩，建模时不要把P_load乘到z外面，否则约束和目标之间的因果链会变得非常难调试。
关键约束在Yalmip里这样写：
Constraints = []; %% 节点功率平衡 for s = 1:ns for t = 1:nt Constraints = [Constraints, ... P_mps(:,t,s) + P_sub(:,t,s) - P_load .* z(:,t,s) == A_bus * P_branch(:,t,s)]; end end %% MPS容量与位置绑定 for i = 1:nb Constraints = [Constraints, ... e(i) <= M * x(i)]; % 容量与位置绑定 Constraints = [Constraints, ... e(i) >= C_min * x(i)]; % 配置位置必须有最小容量 end %% MPS功率输出上限 Constraints = [Constraints, ... P_mps(:,:,s) <= e .* ones(1, nt)];这里面M是足够大的常数，作用是把“没配置的节点容量强制设为0”。这里要注意Big-M的取值不能太大，太大会导致数值病态，建议取最大单台MPS容量的1.1到1.2倍即可。这个细节不处理好的话，Gurobi求解时会出现很多数值警告，最优解精度也会变差。
3.4 求解器调用与参数调优
Yalmip调用Gurobi的语法很简单，核心就是一行：
ops = sdpsettings('verbose', 2, 'solver', 'gurobi', ... 'gurobi.MIPGap', 0.01, 'gurobi.TimeLimit', 3600); optimize(Constraints, objective, ops);这里MIPGap设成0.01，意思是求解器在目标值与最优下界差距小于1%时提前停止。复现论文时别贪心设成0或者1e-4，因为预配置模型在上百个场景下要跑很久，1%的精度已经足够工程使用。TimeLimit设成3600秒是防止某些场景出现完全陷死的情况。
还有一个容易被忽略的参数是线性化可行性容差。Gurobi的MISOCP问题默认对锥约束的可行性要求比较严格，如果约束数值尺度差距太大（比如电压平方约10000伏特平方、功率约几千瓦），建议对模型做“标幺化”处理。把功率和电压都折算到标幺值后，数值尺度会集中在0.1到10之间，求解器内部计算会更稳定，收敛速度提升非常明显。我实测同一个模型，标幺化前后求解时间能差出3倍以上。
4. 实操过程全纪录：从简化版到完整模型
4.1 第一步：单时段单场景先把模型跑通
我复现的第一步是在IEEE 33节点系统上跑一个“单时段、单场景、固定故障”的简化模型。这个模型的规模比完整模型小几十倍，目的纯粹是验证变量定义、约束格式和求解器接口是否正常。
具体做法是：假设某条支路断开，配电网形成一开一合的状态，MPS可以选择在任意节点预配置，目标是让加权失负荷最小。简化版的结果大概在几秒内就能算出来，而且解非常直观：MPS一定会跑到一个既能避开故障线路、又能靠近重要负荷的节点上。看到这个符合物理常识的结果，我才放心往完整模型扩展。
这一步有一个调试技巧：把求解器verbose设置为2，看Gurobi日志输出。如果出现“numerical issues”或者“constraint violation”相关警告，先不要往下扩展，回头检查Big-M取值、变量标幺化、约束是否写错。日志中持续出现No feasible solution也不要慌，先用linprog或者cplex跑一个无整数约束的松弛版本，判断是模型不可行还是求解器没找到可行解。
4.2 第二步：多时段扩展，把时序关系捋清楚
确认简化模型没问题之后，第二步把时间维度加进来。多时段的MPS预配置模型中，MPS的输出功率在每个时段可以独立变化，但它的总能量是有限的，所以需要增加一个“剩余电量”约束：每个时段的累计放电量不能超过预配置容量。
时序约束在Yalmip里实现也不复杂，核心就是循环约束：
Constraints = [Constraints, ... SOC(:,t+1,s) == SOC(:,t,s) - P_mps(:,t,s) * delta_t, ... SOC(:,1,s) == e, ... 0 <= SOC(:,t,s) <= e];这段约束物理意义很清楚：MPS的电量是一个状态变量，每个时段的电量等于上一时段电量减去放电量。delta_t是时段时长，如果你用的是小时级时段，delta_t就等于1。复现时最容易忘记的是把SOC变量的维度方向写对：SOC和e都是nb乘1的向量，不要写成标量。
多时段模型跑通后，我建议先固定预配置方案不变，随机抽几个场景去检查MPS功率输出曲线是否合理，比如是否出现了“配置了MPS的节点在故障后用不上”或者“同一时段两台电源同时给一个节点送电”这类逻辑错误。各类论文里都会加“节点同时接入电源数量不超过1台”这类整数互斥约束，千万不要漏掉。
4.3 第三步：完整多场景模型和结果校验
完整模型把场景数扩展到50个甚至更多。每个场景对应一种故障组合，比如两到三条线路同时断开。场景概率p_s需要提前算好，常见的做法是根据线路故障率和蒙特卡洛抽样得到，或者直接用原文给定的场景集。
这一阶段最大的变化是约束数量成倍增长。50个场景乘24个时段，再加上33节点潮流约束，模型变量数会上到几十万，约束数接近百万量级。这就是为什么我在第一、二步坚持用简化模型验证逻辑——如果直接在完整模型上调试，修改一次约束格式可能就得等半小时才能看到结果。
结果校验我采用三个指标：目标函数值（加权失负荷期望）、MPS配置位置分布、平均空闲容量。目标函数值用来和论文数值做对比，差距在2%以内可以接受；配置位置分布用来做物理合理性检查，比如是否集中在关键负荷节点附近；平均空闲容量用来判断配置是否过度保守，如果所有场景下MPS的利用率都很低，说明预配置方案大概率偏离了最优。
4.4 跑完模型后怎么读结果
我在调试时习惯把每个场景的失负荷结果画成热力图，节点编号是纵轴、时段是横轴、颜色深浅代表失负荷率。这样能一眼看出来模型把负荷损失集中在哪些时段、哪些节点，也能快速判断预配置方案是否把MPS摆在了“每一条可能故障路径的交汇位置”。这个可视化技巧对写论文汇报和答辩特别有用，因为评审看到一张清晰的热力图，比看一行行数字直观得多。
另外，把目标函数中各项拆开看也很重要：是哪个场景主导了整个目标值？如果某一个极端场景的概率不大却贡献了超过一半的失负荷，说明预配置方案过于保守，可能需要调整场景概率设置或者增加对不同场景的鲁棒性要求。
5. 复现中最常踩坑的三类问题
5.1 求解器不收敛或长时间无解
这个问题在MISOCP里很常见，原因大概率不是模型本身不可行，而是数值条件不好。我遇到过两种典型情况：一是Big-M值取得太大，Gurobi在松弛过程中碰到很大数值的系数，数值溢出或精度退化；二是变量单位不统一，功率用千瓦、电压用伏特、阻抗用欧姆，所有量纲混在一起，约束矩阵条件数恶劣。
排查方法分三步。第一步对所有连续变量做标幺化处理，基准功率取系统总容量、基准电压取12.66千伏；第二步检查约束中所有常数项的数值范围，尽量控制在1e-3到1e3之间；第三步开启Gurobi的NumericalFocus参数设为2或3，让求解器主动加强数值处理。
5.2 MPS移动路径约束和预配置逻辑冲突
动态调度阶段需要描述MPS从一个节点移动到另一个节点的过程，而在预配置阶段，往往已经确定了MPS的初始位置。两者衔接时最常见的冲突是：预配置阶段认为某节点已放置MPS，但动态调度阶段又给这个MPS指派了“从远方移动过来”的任务，逻辑上就矛盾了。
解决办法是增加一个特殊的“初始位置固定约束”：在调度阶段，MPS的0时刻位置必须等于预配置阶段给出的方案。这组约束在模型里是纯粹的线性约束，实现没有难度，但需要在建模前就把两个阶段变量关系理清楚。建议把我的建议落在代码注释里，把“which variables come from stage1”和“which variables are stage2”写清楚，避免后期自己看代码都理不清逻辑。
5.3 二阶锥松弛间隙过大的问题
松弛间隙是指在松弛后的凸锥模型求出的目标值，和真实非凸模型的最优解之间存在的偏差。处理这个问题需要区分两种情况。第一种情况是个别场景下松弛精确度依然很高，只有个别场景出现偏差，这种直接忽略或稍微收紧锥约束；第二种情况是系统性偏差比较大，那就需要在目标函数中加入对潮流越限的惩罚项，或者先固定整数变量后重新求解原始非线性潮流方程做校验与修正。
我在复现中还额外测试过一个办法：在预配置模型求解完成后，把得到的整数解固定下来，再用非线性潮流工具（比如Matpower）做一次确定性潮流计算，验证每个场景下电压、功率是否越限。这套“优化加仿真闭环校验”的方法虽然会多花一些计算时间，但能让复现结果的可信度大幅提升，答辩或写报告时也更有底气。
6. 复现完成后我能留下的三条经验
MPS预配置这个模型做完之后，我最大的感受不是“多会了一个公式”，而是更清楚论文复现该怎么推进。第一条经验是：复现电力系统优化类论文，千万别一上来就追求完整复现，一定先做最小化原型。哪怕只是单时段、单节点、单场景的“幼稚模型”，跑通后你也能立刻发现变量维度、约束格式、求解器接口这些基础问题。
第二条经验是：Matlab加Yalmip的组合非常适合做这类模型的快速验证，但工程落地时一定要关注数值标定。Yalmip让你表达模型很容易，可它不会替你去处理量纲不统一、Big-M过大这些“物理层面的脏活”。标幺化、参数归一化、容差调整这些看起来不起眼的操作，往往是复现成败的分水岭。
第三条经验是：复现工作一定要留下“可对答案的基准”。我在跑通完整模型后，专门把论文中的关键参数和对应结果做成了一个参照表，每次修改决策假设、调整场景概率，都把新结果和参照表对比。这样既能保证复现的准确性，后续做改进实验时也有清晰的基线。这种“论文数据对照表”的方式，在学术写作和实验管理里非常值得推广。
最后再分享一个小技巧：画热力图时可以把MPS配置位置用黑色方框叠加在图上，一眼就能看出预配置方案是不是和失负荷热点区域对齐。这个可视化操作成本极低，但对理解模型行为和向别人解释方案逻辑，帮助远超一大堆密密麻麻的表格数字。希望这套复盘对你的复现工作有点启发，动手试起来吧。
