1. 别急着写代码：先把“韧性提升+MPS预配置”这件事想清楚
1.1 配电网韧性和移动电源到底解决什么问题
先说一个很常见的场景：台风过境，或者冰灾压垮线路，配电网最容易出现的情况是“一条主馈线断掉，后面一串负荷全黑”。传统配电网可靠性研究关注的是随机故障，处理起来大多是故障后转供、抢修恢复；但极端天气下往往是多重故障同时发生，转供线路也可能被波及，这时候常规手段就不够用了。所以这十多年，“韧性（Resilience）”这个概念才在电力系统里火起来——它关心的不是故障发生的概率，而是故障发生之后，系统能不能扛住、能不能快速恢复。
移动电源在韧性提升里扮演的是一个“应急替补队员”的角色。不管是移动储能车、应急发电车，还是那种可移动的电池集装箱，它们的共同特点是可以跨区域调度，哪里缺电就开到哪里去。问题是，移动电源数量有限、容量有限，灾害到来前它们停在哪里，直接决定了故障发生后能不能在最短时间内覆盖到最需要的地方。这就是“预配置”的意义：不是等停电了再满世界找车，而是提前把力量布置到关键节点上。
我在复现的那篇SCI一区论文里，核心就是把这种“灾害前决策”建成了一个数学优化模型：第一阶段决定移动电源的预配置位置和容量，第二阶段再根据具体故障场景，去做动态调度。这篇博客先讲第一阶段的预配置怎么做，用Matlab把模型跑通，顺便把我在实现过程中踩过的坑都写出来。如果你在做韧性分析、孤岛划分、移动储能调度或者类似的课题，这篇内容应该能直接帮你省下几天试错时间。
1.2 SCI原论文的思路与复现目标
这类论文的套路其实非常清晰，基本是“两阶段随机优化”或者“两阶段鲁棒优化”的框架。第一阶段是“这里、现在”的决策：我要在哪些候选节点预先放置移动电源，每个点放多大的容量。这个决策必须在灾害信息完全确定之前做出来，所以叫“预配置”。
第二阶段则是“那里、事后”的决策：灾害发生后，具体哪几条线路断了、哪些负荷失电，这时候移动电源怎么移动、接到哪个节点、输出多少功率。但因为第二阶段还没真正发生，第一阶段做决策时必须把第二阶段所有可能的应对成本都算进去，于是就有了“期望成本”或“最坏情况成本”。
我的复现目标是把第一阶段模型写成一个可解的混合整数线性规划（MILP），用Matlab + Yalmip + Gurobi这套组合求解。具体来说，给定一个配电网算例，一组负荷曲线和一组故障场景，模型要输出：在哪些节点配置几台移动电源。至于故障发生后的动态调度，我打算放到下一篇博客里单独讲，先把预配置这个地基打好。
1.3 复现时的工作范围与输入输出
我在动手前先把输入输出列了个清单，这能避免写着写着就跑偏。输入方面，需要准备：
- 配电网拓扑：节点编号、支路编号、线路电阻电抗、连接关系。我用了IEEE 33节点系统，这个算例资料多、规模适中，很适合做算法验证。
- 负荷数据：每个节点的有功、无功负荷，以及峰谷时段的变化曲线。一般用24小时或几个典型时段。
- 移动电源参数：单台MPS的容量（kWh）、额定功率（kW）、数量上限、每节点最多配置几台。
- 故障场景：灾害可能破坏哪些线路，我用随机抽样生成若干场景，每个场景对应一组断线集合。
- 候选预配置点：不是所有节点都适合停移动电源，所以模型里我预先筛掉一部分出线少、路况差的节点，只保留一部分作为候选。
输出则很明确：一份预配置方案，包括每个候选点配置的MPS数量。同时还会输出对应的期望负荷削减量、配置成本等中间结果，用来验证模型行为是否合理。
2. 数学模型拆解：预配置问题的核心是“把移动电源放在哪、放多少”
2.1 目标函数：最小化什么，最大化什么
刚开始复现时，最容易犯的错是看到“韧性提升”就觉得目标函数应该是“最大化恢复负荷”。但实际建模往往不是这么单独设目标的，因为移动电源有成本、数量有限，你不能无限往网络上塞。所以论文里的常见做法是：最小化总成本，包括预配置成本和灾害后的损失期望。
我这里采用的目标函数是这个样子：
[ \min \quad \sum_{i \in N_{cand}} c_{pre} \cdot x_i + \sum_{s \in S} w_s \cdot \sum_{t \in T} \sum_{j \in N_{load}} p_{s,t,j}^{curt} ]
其中：
- (x_i) 是在候选节点 (i) 预配置的MPS台数，整数变量。
- (c_{pre}) 是单台MPS的预配置成本系数，包括租赁费、场地费等。
- (S) 是故障场景集合，每个场景有概率权重 (w_s)（通常等概率，也可以按灾害强度设置）。
- (p_{s,t,j}^{curt}) 是场景 (s) 在第 (t) 个时段在节点 (j) 的负荷削减量，这个变量由第二阶段决定。
预配置成本决定了“放多少”，而场景期望削减量决定了“放哪里更有价值”。这两个量需要设置合适的量纲和权重。如果配置成本太高，模型会干脆不配任何MPS；如果太低，模型就会把所有候选点都塞满。这不算模型错了，而是参数调校的问题，我在第四章会专门讲。
2.2 关键约束：负荷、容量、时序、连通性
模型里不能只有目标函数，约束才是让结果符合物理规律的关键。我复现时主要考虑了四类约束。
第一类，预配置数量约束。候选节点上配置的MPS台数不能超过该节点最大允许数，且总台数不能超过可用MPS总数：
[ 0 \le x_i \le \bar{x}_i, \quad \sum_i x_i \le M^{total} ]
第二类，MPS接入能力约束。第二阶段里，一台MPS只能在某个时段接入到某个节点，不能同一时间出现在两个地方。这个约束对预配置本身不直接起作用，但它是第二阶段调度模型的基础。因为本篇复现只聚焦预配置，我为了验证方案可行性，把每个故障场景简化成“故障后系统进入孤岛运行，MPS可以直接向相邻负荷供电”，然后通过一组连续潮流近似来检验。
第三类，功率平衡与负荷削减约束。在场景 (s)、时段 (t) 中，节点 (j) 的恢复供电量加上负荷削减量要等于原始负荷值：
[ p_{s,t,j}^{rec} + p_{s,t,j}^{curt} = P_{t,j}^{load} ]
且每个节点的恢复供电量不能超过该节点可用的MPS出力总和：
[ p_{s,t,j}^{rec} \le \sum_{k \in \Phi(j)} P^{MPS} \cdot z_{s,t,k} ]
其中 (z_{s,t,k}) 是MPS在场景 (s) 时段 (t) 是否接入节点 (k) 的0-1变量，(\Phi(j)) 表示能覆盖节点 (j) 的MPS接入点集合。这里隐含了“移动电源接入某个节点后，只对该节点或者其下游孤岛供电”的简化，对预配置阶段的评估足够了。
第四类，连通性约束。故障后恢复的区域必须满足辐射状网络的要求，不能出现环网，也不能形成不可行的孤岛。严格意义上是二阶锥约束或者角度约束，但为了保持MILP可解，我采用“单商品流”方法，把连通性约束线性化。
2.3 参数设定与场景生成
参数设定直接影响结果质量，不能拍脑袋。我用IEEE 33节点系统，基准电压12.66kV，总负荷约3.7MW。移动电源单台参数参考市面上常见的移动储能车：额定功率100kW，容量400kWh，单台预配置成本设为8000元（这个值只是为了验证流程，实际项目可以按租赁合同修改）。候选节点选了8个：5、8、13、18、22、25、29、32，这些节点分散在各分支末端和关键连接处，覆盖效果比较好。
故障场景生成我一般这么做：设定一个灾害强度，比如“同时故障线路数为3条”，每次随机从全部33条支路中抽3条断开，生成50个场景。抽样时注意不要让同一条线路反复出现在同一个场景里，也不要出现两个场景完全一样。场景数量太少会让结果随机性很大，太多则MILP规模爆炸。我实测50个场景对于33节点系统已经能把预配置位置稳定下来了。
2.4 求解器与线性化处理
预配置本身是整数决策，再加上每个场景里的连续变量和0-1变量，这个模型是一个典型的MILP。求解器我用的是Gurobi，专门对付这种大规模整数规划。Matlab里用Yalmip做建模语言，把优化问题写成人类可读的形式，再交给Gurobi求解。
为什么不用现成的Matlab工具箱？因为intlinprog虽然能解MILP，但遇到非线性或复杂约束时需要手动构造稀疏矩阵，非常痛苦。Yalmip允许直接用x = binvar(n,1)定义0-1变量，用Constraints = [...]一条条叠加，代码可读性高很多，改起模型来也快。
但要注意，Yalmip本身不是为了求解，它只是个“翻译器”。你要先把模型中有平方项、绝对值项、逻辑或关系等非线性全部线性化，才能交给Gurobi。比如目标里如果出现“负荷削减量的期望”，那是线性项，没问题；但如果有“MPS接入节点后最多同时通过一个子区域供电”这类逻辑，就需要引入辅助0-1变量和big-M惩罚。
3. Matlab代码实现：从零搭一个MPS预配置求解框架
3.1 数据准备：节点系统、负荷曲线与移动电源参数
我用的是Matpower的case33bw算例，不过在Matpower里加载后还要自己提取节点和支路数据。第一步先把数据转换成求解模型用的结构体。
mpc = loadcase('case33bw'); baseMVA = mpc.baseMVA; bus = mpc.bus; branch = mpc.branch; % 节点导纳、拓扑 nBus = size(bus, 1); nBranch = size(branch, 1); fromBus = branch(:, 1); toBus = branch(:, 2); % 负荷：直接把节点有功、无功取出来 Pload = bus(:, 3) / baseMVA; % 转化为标幺值 % 更推荐保留有名值，方便后续读结果 Pload_kW = bus(:, 3);这里有个小坑：Matpower里节点序号不一定是连续整数，有的算例节点编号会跳。所以我在构建所有矩阵时都用”序号数组“去映射，而不是直接用节点编号做索引，否则后面定义Yalmip变量时会越界。
移动电源参数我用一个结构体存：
MPS.comm = 6; % 可用移动电源总数 MPS.P = 100; % 单台额定功率 kW MPS.E = 400; % 单台容量 kWh MPS.costPre = 8000; % 单台预配置成本 MPS.maxAtNode = 2; % 单个候选点最多放置台数候选节点直接写在数组里：
candNodes = [5 8 13 18 22 25 29 32];负荷曲线这一块，我直接对每时段做了一个缩放系数：白天1.0，晚上0.6，凌晨0.3。如果做更精细的研究，可以填入真实的24小时负荷曲线。但在预配置阶段，用几个典型时段的峰值负荷就够了，因为我们要保证移动电源在最坏负荷时段也能顶上。
3.2 决策变量与约束构建的写法
预配置变量用整数变量，每个候选节点一个：
x = intvar(1, length(candNodes), 'full'); Constraints = []; for i = 1:length(candNodes) Constraints = [Constraints, 0 <= x(i) <= MPS.maxAtNode]; end Constraints = [Constraints, sum(x) <= MPS.comm];接下来是每个故障场景的变量。由于故障场景里要评估每个时段每个节点的供电恢复情况，我定义了一个Pcurt的三维变量。Yalmip支持x = sdpvar(nBus, nTime)，场景多了就用循环内部定义。
Pcurt = sdpvar(nBus, nTime, 'full'); % 恢复功率也定义 Prec = sdpvar(nBus, nTime, 'full');然后在场景循环里加入功率平衡约束：
for t = 1:nTime for j = 1:nBus Constraints = [Constraints, Prec(j,t) + Pcurt(j,t) == Pload_curve(j,t)]; Constraints = [Constraints, Prec(j,t) >= 0, Pcurt(j,t) >= 0]; % 恢复功率上限：该节点所在孤岛可用的MPS总功率 Constraints = [Constraints, Prec(j,t) <= MPS.P * sum(x(find(candNodes == ...)))]; end end注意上面这个写法和严谨模型还有距离，尤其是”故障后哪些负荷能由哪些MPS供电“，需要根据拓扑连通性去计算。我在简化版本里先设定：每个候选MPS接入点只对其“下游侧”负荷有供电能力，这个关系可以从网络拓扑推出来。实际写代码时，我会预计算一个reachableMatrix(nCand, nBus)，标识候选接入点与负荷节点的连通关系，然后再乘上x。这样就避免了在Yalmip约束里做动态索引。
大体约束框架就是如此。因为每个场景都包含Prec和Pcurt，50个场景下来变量数会到几千个，好在Yalmip处理这个规模很轻松。
3.3 求解流程与结果输出
约束装完之后，直接调optimize：
ops = sdpsettings('solver', 'gurobi', 'verbose', 2, 'showprogress', 1); optimize(Constraints, Objective, ops);目标函数的拼装我放在场景循环里：
Objective = sum(MPS.costPre * x); for s = 1:nScenario Objective = Objective + w_s * sum(sum(Pcurt_s{s})); end求解结束后，提取x的办法是：
x_sol = value(x); fprintf('预配置方案：\n'); for i = 1:length(candNodes) if x_sol(i) > 0 fprintf('节点 %d 配置 %d 台MPS\n', candNodes(i), x_sol(i)); end end整个求解过程在50场景、3个时段、33节点系统上，Gurobi大约10到30秒就能出结果。如果场景数加到200，时间会明显上升，这时候就需要考虑场景削减策略了。
3.4 结果可视化：地图上看出预配置位置
先把网络拓扑画出来，再把预配置位置标出来，这是很直观的验证方式。我用Matlab自带的plot和line画节点和支路，再用红色方块标出MPS：
figure; hold on; % 画支路 for k = 1:nBranch x1 = bus( fromBus(k), 6 ); y1 = bus( fromBus(k), 7 ); x2 = bus( toBus(k) , 6 ); y2 = bus( toBus(k) , 7 ); plot([x1 x2], [y1 y2], 'k-', 'LineWidth', 1); end % 画节点 plot(bus(:,6), bus(:,7), 'bo', 'MarkerSize', 4); % 画MPS位置 idx = find(x_sol > 0); plot(bus(candNodes(idx), 6), bus(candNodes(idx), 7), 'rs', 'MarkerSize', 10, 'LineWidth', 2);如果结果科学，红色方块应该出现在网络的关键分岔点或负荷较重、恢复路径容易切断的地方，而不是集中在变电站附近。万一看到MPS全堆在电源侧，那大概率是约束没加对，或者场景选择太温和，第四章会聊这个问题。
4. 踩坑实录与复现心得
4.1 求解器报错与调试技巧
我最常遇到的报错是Yalmip提示NaN in constraints或者Unable to evaluate。这种问题九成是变量维度不匹配导致的。比如我早期用Pcurt(j,t) == Pload_curve(j,t)时，Pload_curve是33×3，但循环里j是从1到33，t是从1到3，看着没问题，可一旦Pload_curve中有节点编号错位，就会把空的索引传给Yalmip。
调试办法很简单：在约束构建完后，用check(Constraints)看每行约束的状态，Yalmip会返回约束是否可行。然后逐步注释掉可疑约束，定位到出问题的那一行。还有一个实用习惯是定义变量前先把输入数据的维度打印一遍，尤其是bus、branch来自Matpower这种外部工具时，千万不要直接信任它的编号顺序。
Gurobi这边也常出现Q matrix is not positive definite之类的报错，这通常说明我把某两个0-1变量相乘了，产生了二次项。预配置模型必须是线性的，遇到变量相乘，必须用辅助变量和big-M法线性化。
4.2 结果不合理？检查这几个地方
如果模型跑完，发现所有节点都不配置MPS，第一反应是“程序bug”，但多数时候其实是目标函数里预配置成本相对停电损失太高了。这时把MPS.costPre调小，或者把负荷削减的惩罚单价调高，比如把单位停电损失设为每千瓦时20元，很快就会出现MPS配置。
另一种极端情况是几乎每个候选点都塞满了MPS。这说明惩罚成本太夸张，或者MPS总数约束没写进去。我建议先给目标函数里的两类成本定一个数量级参考：单台MPS一天租赁费如果是8000元，而它一天最多能减少的电量约400kWh，那停电损失单价如果超过20元/kWh，模型就会倾向把所有MPS都用上。所以调权重不是瞎调，得按经济指标来。
还有一种隐蔽问题：MPS设置了，但第二阶段的恢复功率上不去。原因往往是我只约束了“MPS总功率覆盖负荷”，但没有考虑移动电源自身的电量约束。预配置阶段如果不看电量，模型就可能把移动电源放在一个重要节点，但故障发生后根本没有足够时间对负荷持续供电。后来我在预配置模型里加入了单位时段最大放电量约束，才让结果更可信。
4.3 如何把预配置结果扩展到动态调度
虽然这篇只讲预配置，但我还是要提一句动态调度，因为它和预配置是强耦合的。预配置解决的是“在哪里放”，而动态调度解决的是“故障发生后，第t小时的MPS移到哪里去”。通常动态调度模型会引入0-1变量z_{s,t,k}，表示在场景s、时段t，MPS是否接入节点k。这个变量和预配置决策x_i之间的关系是：任何时段，MPS接入点必须来自预配置集合，且一台MPS在时段t只能在一个节点：
[ \sum_k z_{s,t,k} \le \sum_i x_i, \quad \forall s,t ]
[ \sum_k z_{s,t,k} \le M \cdot \sum_{i \in N_{pre}} x_i ]
通过调整z变量在不同时段的变化，就能模拟MPS从A点到B点的移动路径，同时还要考虑移动时间约束和MPS电量时序转移。这部分我准备在下一篇博客里把完整代码放出来，包括怎么把故障场景下的时间序列负荷和分时段调度写进去。这里先记住一个点：预配置模型里一定要把MPS的总数和候选位置写成参数和变量，这样动态调度模型才能直接调用，不然两阶段的接口对不上，还得返工。
4.4 个人经验：算例规模与运行时间控制
我跑过IEEE 33节点和IEEE 123节点两种系统。33节点配50个场景，Gurobi很快；123节点配100个场景，MILP变量数量会到几万个，求解时间可能从几十秒变成几十分钟。这种时候有几个非常有效的瘦身方法。
第一，场景削减。与其随机生成100个相似场景，不如用k-means对故障线路集合聚类，选出20到30个有代表性的场景，再给每个场景分配概率。这样可以保持精度的同时大幅减少变量数。
第二，候选节点预筛选。不是每个节点都适合放MPS。我会先用一个启发式规则：计算每个节点的“失电负荷节点数”，如果节点处于网络末端且下游负荷很小，把它从候选集合里去掉。候选点从33个减到8个以后，整数变量一下子少了四分之三，求解速度提升明显。
第三，时段聚合。24小时调度问题如果按小时建模，变量数会涨24倍。但对预配置评估来说，白天、夜间、凌晨三个典型时段已经足够，把时间维度从24个小时聚合成3个时段，模型变得非常清爽。动态调度阶段当然需要更细的时间分辨率，那时候再单独细化。
我在实际操作中发现，这三个方法组合起来，123节点系统也能在5分钟内求出较优解，而直接暴力建模可能需要两小时还解不到gap 1%。所以遇到大规模算例，先别急着加服务器内存，先把模型“减脂”再求解。
复现这篇SCI论文上下两篇对我来说收获真的很大。尤其预配置这一阶段，看起来只是“放几台车”的小问题，但真正建模时会牵扯到故障场景、时序负荷、MILP求解效率等一大堆东西。这篇把预配置部分代码思路和调试过程都讲了一遍，希望能给你省点时间。如果你的算例规模更大，或者想换成鲁棒优化框架，关键改动其实只在于把期望项换成max-min结构，但Yalmip和Gurobi这套求解路径是一样的。下一篇动态调度，我会把z_{s,t,k}的建模细节、MPS移动路径约束和分时功率分配代码补上，到时候可以在本文的预配置结果基础上直接接着跑。
