企业官网建设流程全解析

台风还有48小时登陆，你手里有3台移动储能车，该把车提前停到哪几个节点？这道题看起来像"选个车位"，实际上是把配电网韧性、随机优化、时序耦合全部串起来的硬骨头。我最近在IEEE33节点系统上完整复现了一套"移动储能预布局+动态调度"策略，跑通的那天凌晨一点半，Gurobi的求解日志停在最优解上，我看着失负荷曲线一点一点降下来，确实挺有成就感。这篇就把整个复现过程、数学模型、代码逻辑和踩过的坑全部摊开讲，给同方向的研究生和工程师做个参考。

这篇博文适合三类人：一是正在做配电网韧性方向课题的研究生，想找一套能跑通的两阶段优化写法；二是做配网调度系统、应急保电方案设计的工程师，想评估移动储能实际能发挥多大作用；三是想了解IEEE33节点系统、YALMIP建模和线性化处理技巧的Matlab玩家。全文不绕弯子，只说从实际运行调试中得到的结论和经验。

1. 预布局和动态调度到底在解决什么——先把这个策略的本质逻辑捋清楚

1.1 配电网韧性和常规可靠性是两套评估逻辑

很多人刚接触"韧性"这个词时，容易和配电网可靠性混在一起。可靠性关注的是日常大概率故障下的供电能力，比如N-1校验、转供成功率、平均停电时间这类指标，这一套在配电自动化里已经非常成熟。韧性关注的是极端小概率、大影响事件下的表现，比如台风刮倒一排电杆、冰灾压断多条支路，此时配电网已经不是一个"失去单个元件"的状态，而是出现多个故障点、甚至局部孤岛的情况。韧性评估用的是"事件发生前—事件冲击中—恢复过程中"三阶段曲线，看的是系统能扛住多少冲击、损失多大、恢复多快。

这就带来一个本质区别：常规可靠性优化是在"正常拓扑+少数故障"的假设下做调度方案，而韧性提升必须在"最坏情况下怎么保住关键负荷"的思维下做决策。移动储能在这种场景下的价值就出来了——台风来临前把电池运到重负荷、孤岛后能独立供电的关键节点，故障发生后再根据实际断线情况调整接入点，这就是标题里"预布局"和"动态调度"的核心分工。

1.2 移动储能的核心优势不是"能储电"，而是"能移动"

固定储能的位置一旦安装就锁死，而移动储能(MESS)装备在车辆上，在时间尺度上引入了空间自由度。理解这一点需要拿消防队做类比：一栋楼装了固定灭火器当然有用，但消防车可以提前停在风险高的片区、火灾发生后快速机动到最需要的地方，这种"可重定位"能力在事件演变不确定时价值极大。

回到配电网的数学建模上，这个"可移动"特性给优化模型引入了两个麻烦：一是储能接入位置的决策变量从连续变成了离散（它到底停在哪，是0-1变量）；二是储能从一个节点转移到另一个节点需要时间和路径约束（移动期间既不能充电也不能放电），这就让模型从"纯容量/功率约束"上升到"时空耦合约束"。我在复现过程中最大的体会是：固定储能模型写起来半小时，移动储能模型光约束推导就占了整整两页纸，后面每一个求解报错都跟这些时空约束脱不开干系。

1.3 这套策略能回答的关键工程问题

复现这套策略，最终是要回答三个工程问题：灾前气象信息还不完全确定时，移动储能的最佳预置点在哪几个节点；灾中实际故障状态慢慢暴露后，哪些储能车应该继续待在原位、哪些应该转移去支援新的孤岛；以及牵扯到多台储能车时，容量分配、充放电时序和移动路径之间怎么协同才能保证关键负荷损失最小。这三个问题恰好对应两阶段优化模型的三个核心决策块，也是整篇代码里最值得花时间理解的部分。

2. IEEE33节点系统与移动储能数学建模——复现前必须打好的地基

2.1 IEEE33节点：标准测试系统的参数从哪来、怎么看

IEEE33节点配电网是配电网研究中出现频率最高的测试系统之一，33个节点、32条支路的辐射状结构，基准电压12.66kV，额定总有功负荷大约3.7MW左右、无功2.3Mvar。它不算大，但节点和支路数量足够体现优化模型的组合复杂度，所以非常适合做移动储能预布局这种带0-1变量的规划问题。

具体参数获取有两个常用渠道：一是Matpower自带的case33bw数据文件，装好Matpower后用loadcase('case33bw')直接导入，里面包含了完整的支路电阻、电抗、对地电纳和节点负荷参数，是我这次复现主要用的数据源；二是很多论文附录会给出完整参数表，逐条录入也行，但手输容易出错，我建议能用数据文件就尽量用数据文件。拿到数据后第一步要做的是把支路首末端节点号、阻抗值、各节点有功无功负荷整理成四个标准矩阵，后面建模、写约束、可视化都靠它们。

2.2 移动储能约束：时空耦合是建模的核心难点

移动储能的数学模型可以拆成三个层次来看。

第一层是储能本体的运行约束，和固定储能一致。设第(b)台MESS在时段(t)的荷电状态为(S_{b,t})，充放电功率分别为(P^{ch}{b,t})和(P^{dis}{b,t})，电池容量为(C_b)，那么SOC递推公式就是：

[ S_{b,t+1} = S_{b,t} + \left( \eta^{ch} P^{ch}{b,t} - \frac{P^{dis}{b,t}}{\eta^{dis}} \right) \frac{\Delta t}{C_b} ]

充电功率和放电功率不能同时为正，并且受额定功率上限(P^{max}b)约束：(0 \le P^{ch}{b,t} \le u^{ch}{b,t} P^{max}b)，(0 \le P^{dis}{b,t} \le u^{dis}{b,t} P^{max}b)，其中(u^{ch}{b,t} + u^{dis}_{b,t} \le 1)。

第二层是位置约束。定义0-1变量(x_{b,n,t})表示MESS (b)在时段(t)是否接入节点(n)，那么任一时刻每台MESS最多接入一个节点：(\sum_{n \in \Omega_N} x_{b,n,t} = 1)。这个约束看似简单，却是整个模型里把"可移动"落到数学空间的关键。

第三层是转移约束，这是和固定储能最大的区别。从节点(i)移动到节点(j)需要时间(T^{move}_{ij})，在移动期间储能处于不可用状态，因此需要保证：

[ x_{b,i,t} = x_{b,j,t + T^{move}_{ij}} = 1 \quad \text{且移动中间时段不接入任何节点} ]

实际建模中更常用的写法是引入0-1移动变量(y_{b,ij,t})表示MESS在时段(t)开始从(i)向(j)转移，然后建立位置变量和移动变量之间的耦合关系：

[ x_{b,i,t} - x_{b,i,t+1} \le \sum_{j \in \Omega_N, j \neq i} y_{b,ij,t} + M(1 - z_{b,t}) ]

在代码里，这一步我走了不少弯路——一开始直接用"时刻t在A、时刻t+1在B"来约束转移，结果模型忽略了移动过程的中间时段，导致算出来的方案里电池一边移动一边还在供电，完全不物理。后来改成引入辅助变量、对移动过程做时段占用约束才算解决。

2.3 目标函数和整体约束的全貌

这套模型的目标函数不是简单的运行成本最小，而是直接瞄准韧性——最小化所有故障场景下的加权失负荷量：

[ \min \sum_{k \in \Omega_S} \pi_k \sum_{t \in T} \sum_{i \in \Omega_N} w_i \cdot P^{shed}_{i,t,k} ]

其中(\pi_k)是场景(k)的概率，(w_i)是节点负荷的重要度权重，(P^{shed}_{i,t,k})是场景(k)下时段(t)节点(i)的切负荷量。用失负荷期望作为目标函数，好处是结果可以直观地和韧性曲线对应起来，坏处是数值量级往往很小，求解器收敛判据如果不小心设置，容易被提前终止并得到次优解，这点后文还会提。

整体约束除了储能本身的约束，还包括配电网潮流约束。这里用DistFlow支路潮流方程，适用于辐射状配电网，节点(i)的有功、无功平衡表示为：

[ P_{i,t,k} = \sum_{j: (j,i) \in E} P_{ji,t,k} - \sum_{l: (i,l) \in E} (P_{il,t,k} - r_{il} I^2_{il,t,k}) + P^{MESS}{i,t,k} + P^{shed}{i,t,k} ]

节点电压满足(V^2_{i,t,k} = V^2_{j,t,k} - 2(r_{ij}P_{ij,t,k} + x_{ij}Q_{ij,t,k}) + (r^2_{ij} + x^2_{ij})I^2_{ij,t,k})，联合起来还需要线路载流量约束和节点电压上下限约束。后面的线性化处理中，这个非线性方程是首先要处理的对象。

3. 两阶段优化框架：预布局与动态调度是怎样咬合在一起的

3.1 为什么必须用两阶段，而不是一阶段全规划好

最初我想过一个"一步到位"的建模方案：直接把所有时段的(x_{b,n,t})都设成决策变量，一次性求出每台MESS在每个时段的位置和功率。这个思路理论上没错，但实际求解时问题很大。原因在于信息结构和决策时序——预布局阶段是在灾害发生前做的决策，基于的是台风的预测路径；而动态调度是灾中实时决策，基于的是已经暴露的实际故障信息。两个阶段的决策时序和信息基础不同，塞进同一个优化问题里，等于让模型提前"看到"了未来故障的实现情况，得到的预布局方案在现实中根本不可行。

两阶段框架的正确逻辑是：第一阶段在灾害前确定MESS的初始布局位置（这里，"位置"是不确定场景的"前瞻"决策），随后第二阶段在各个故障场景下，基于已知的预布局位置去优化MESS的转移路径、充放电功率和切负荷量。翻译成数学语言就是"在这里决策"—"在那里看到结果"—"再在这里调整"的关系。

3.2 主问题与子问题的分工协作

两阶段模型落地到代码，用的是主问题-子问题迭代结构。主问题只求解第一阶段变量（MESS预布局位置矩阵），不考虑具体故障场景的潮流和切负荷细节；子问题根据主问题给的位置，循环求解每个故障场景下的动态调度子问题，并把切负荷目标值反馈回主问题，形成割约束，循环迭代直至收敛。这个过程和列与约束生成(CCG)算法很像，但因为我们算例就33节点，场景数量几十个，直接采用"枚举场景+全场景耦合求解"的方式也能在可接受时间内得到最优解，不需要太复杂的分解算法。

顺带一提，我之前在写这段逻辑时经常搞混"先给位置后求调度"和"位置调度一起求"的边界，导致每次主问题更新位置后，子问题里预布局位置的变量没有同步锁定，结果子问题偷偷把位置也改了，算出来的目标值偏乐观。解决办法很简单——子问题里把第一阶段变量全部设成常数，一个都不能漏。

3.3 故障场景生成：蒙特卡洛抽样和典型场景选择

场景是驱动两阶段模型运转的"燃料"。最常见的方法是对每个线路设一个故障概率，假设台风期间线路独立故障，用蒙特卡洛抽样大量生成场景，再从中选择概率较高或者失负荷较为严重的场景子集进入优化。但直接抽样得到的场景集合往往大而冗余——有些场景拓扑结构几乎相同，计算代价却成倍增加。所以更实际的做法是先用1000-5000次蒙特卡洛抽样得到初步的场景池，统计各线路故障频次和负荷损失，再用K-means或者K-medoids聚类手段提取10-20个典型场景，把聚类后场景的发生概率按簇内样本比例折算。

我在这个环节看到过不少复现代码直接用均匀随机概率抽故障线路组合，从不做场景缩减，导致50个场景求解时子问题计算量巨大，到后面根本算不动。实际调试下来，10-15个典型场景已经能稳定体现韧性差异，再往上加场景对结果的影响非常小，但求解时间基本呈线性增长。

4. Matlab+YALMIP代码实现：核心逻辑逐段拆解与求解器配置

4.1 代码整体架构和模块划分

我复现时的代码结构是四个模块划分：参数数据脚本、场景生成脚本、两阶段优化主脚本、结果可视化脚本。所有数据输入集中在最前面的参数模块里，后续脚本只需调用，改算例时不用一寸寸翻代码。

参数模块的核心是IEEE33节点的线路参数和负荷数据。我从Matpower的case33bw里抽出数据，组织成Branch和Load两个结构体数组，顺便定义了储能参数：

%% IEEE33节点基础参数读取 mpc = loadcase('case33bw'); Branch = [mpc.branch(:, 1:2), mpc.branch(:, 3:4)]; % [首端, 末端, R(pu), X(pu)] Load = [mpc.bus(:, 1), mpc.bus(:, 3), mpc.bus(:, 4)]; % [节点, P(pu), Q(pu)] %% 移动储能参数 MESS = struct(); MESS.num = 3; % 移动储能台数 MESS.Cap = 1.0; % 储能容量 1.0 p.u.（基准功率100MVA下对应10MWh） MESS.Pmax = 0.25; % 最大充放电功率 0.25 p.u. MESS.eta = 0.95; % 充放电效率 MESS.SOC0 = 0.5; % 初始荷电状态 MESS.SOCmin = 0.1; MESS.SOCmax = 0.9; MESS.Tmove = 1; % 单位节点间移动所需时间步长（可扩展为节点距离矩阵）

这里在数据组织上有个值得注意的地方：原始case33bw里的单位是有名值，我全部转成标幺值，这样潮流方程数值比较干净，对求解器的数值稳定性大有帮助。

4.2 第一阶段预布局模型写法

第一阶段模型用YALMIP定义，核心是位置0-1变量和最少接入节点数量的约束：

%% 第一阶段：预布局决策 X_pre = binvar(33, MESS.num, 'full'); % X_pre(n,b)=1表示储能b预布局在节点n Constraints_pre = []; % 每台MESS只能预布局在一个节点 for b = 1:MESS.num Constraints_pre = [Constraints_pre, sum(X_pre(:, b)) == 1]; end % 预布局节点不能是联络开关节点（可选约束） % 预布局任意两台MESS不能在同一节点（可选，实际风灾场景可放宽） for b1 = 1:MESS.num for b2 = b1+1:MESS.num Constraints_pre = [Constraints_pre, ... sum(X_pre(:, b1) .* X_pre(:, b2)) == 0]; % 线性化后处理 end end Objective_pre = % 传递给子问题的辅助目标，通常设为0或者最小化接入距离

真实场景下，"两台MESS不能在同一节点"是非线性约束，我采用了线性化展开的方式转化：定义辅助0-1变量(z_{n,b1,b2})约束(z \ge x_{n,b1}+x_{n,b2}-1)，然后求(\sum_{n} z_{n,b1,b2} = 0)。很多代码里为了省事直接用乘积形式，在小规模算例中也能跑出来，但YALMIP遇到非线性0-1乘积时会引入额外二元变量，求解效率明显下降。我全部统一改成线性化写法后，MIP的求解速度提高了快一倍。

4.3 第二阶段动态调度模型写法

第二阶段是最容易让人头大的部分，因为要处理整个时间轴上的MESS状态变化和潮流约束。核心也是用sdpvar定义连续变量、用binvar定义位置和移动变量，然后一步步堆约束。

%% 第二阶段：动态调度（针对给定故障场景） P_dch = sdpvar(33, T, MESS.num, 'full'); % 放电功率 P_ch = sdpvar(33, T, MESS.num, 'full'); % 充电功率 P_shed = sdpvar(33, T, 'full'); % 切负荷量 SOC = sdpvar(MESS.num, T+1, 'full'); % 荷电状态 X_dyn = binvar(33, T, MESS.num, 'full'); % 动态位置变量 Constraints = []; % 初始位置锁定为第一阶段结果 for n = 1:33 for b = 1:MESS.num Constraints = [Constraints, X_dyn(n, 1, b) == X_pre(n, b)]; end end % 每时段每台MESS只在一个节点 for t = 1:T for b = 1:MESS.num Constraints = [Constraints, sum(X_dyn(:, t, b)) == 1]; end end % 充放电功率上限与位置耦合 for t = 1:T for b = 1:MESS.num for n = 1:33 Constraints = [Constraints, ... P_dch(n, t, b) <= MESS.Pmax * X_dyn(n, t, b)]; Constraints = [Constraints, ... P_ch(n, t, b) <= MESS.Pmax * X_dyn(n, t, b)]; Constraints = [Constraints, ... P_dch(n, t, b) >= 0, P_ch(n, t, b) >= 0]; end end end

这里我踩过一个很严重的坑：直接把充放电变量定义成(节点数 × 时段数 × 储能数)的三维变量，YALMIP里三维sdpvar旋转很麻烦，后来改成了把变量压平为二维，即用一个大索引n*T*T替代，代码才清爽了许多。所以如果你也打算复现，建议第一步就把变量维度规划好，三维变量实在不方便就别硬用。

4.4 DistFlow线性化和二阶锥松弛处理

DistFlow潮流方程里有(V_i^2)和(I_{ij}^2)两个二次项，处理起来最常规的做法是用电压平方(U_i = V_i^2)和电流平方(L_{ij} = I_{ij}^2)替换，然后利用二阶锥松弛把等式松弛成不等式。对一个支路((i,j))，二阶锥约束写作：

[ \left| \begin{matrix} 2P_{ij} \ 2Q_{ij} \ L_{ij} - U_i \end{matrix} \right|2 \le L{ij} + U_i ]

在YALMIP里，这个约束用cone函数一行就能完成：

Constraints = [Constraints, cone([2*P_ij; 2*Q_ij; L_ij - U_i], L_ij + U_i)];

二阶锥松弛的精确性对配电网来说通常是有保障的（满足一定的网络条件），但别掉以轻心——求解完后要检查支路电流松弛间隙，如果某些支路松弛间隙过大，计算结果虽然"最优"但不物理。我复现时发现末端节点附近的松弛误差有时会偏大，解决办法是对相关支路加大惩罚系数或者改内点法参数重新求解。

4.5 求解器配置：Gurobi参数调优实战

这个模型是混合整数二阶锥规划(MISOCP)，我用的求解器是Gurobi，通过YALMIP调用。求解器配置里最关键的有三个参数：

options = sdpsettings('solver', 'gurobi', ... 'gurobi.MIPGap', 0.01, ... 'gurobi.TimeLimit', 1800, ... 'gurobi.NumericFocus', 1, ... 'gurobi.MIPFocus', 2, ... 'verbose', 1);

MIPGap设为0.01，即1%的优化间隙就够了——论文复现不是生产系统，追求0.01%的绝对最优纯属浪费计算资源，1%的间隙完全不影响结论判断。NumericFocus设为1用来增强数值稳定性，IEEE33的标幺值数据本身还算规整，但这个设置能预防一些奇怪的数值警告。MIPFocus设为2让求解器偏向于证明最优性而非快速找可行解，因为两阶段框架里我们更需要确切的上下界来判断循环是否收敛。

有时会遇到Gurobi许可证到期的问题，这时候可以换开源的SCIP求解器。但说句实话，SCIP在MISOCP上的性能跟Gurobi差距明显，小算例勉强可以，场景一多就难顶。有条件还是用Gurobi或者Cplex，学生可以用学术许可证。

4.6 完整的两阶段迭代求解流程

整个主脚本的循环逻辑是：

%% 两阶段迭代求解主循环 LB = -inf; UB = inf; iter = 0; X_pre_current = initial_guess(); % 初始可行解：等概率随机布点 while (UB - LB) / max(1e-6, abs(UB)) > 0.01 && iter < 10 iter = iter + 1; % 固定X_pre，求解所有场景的动态调度子问题 total_obj = 0; for s = 1:S [obj_s, P_shed_s, MESS_result_s] = solve_subproblem(X_pre_current, scenario(s)); total_obj = total_obj + scenario(s).prob * obj_s; end UB = total_obj; % 添加Benders割约束到主问题（或直接采用CCG框架，在此简略） % 重新求解主问题，得到新的X_pre X_pre_current = solve_masterproblem(added_cuts); % 主问题目标值作为新的LB（简化处理，严格意义上有改进） LB = get_master_objective(); end

如果场景数不多且目标是复现文献结果，甚至可以绕开这个迭代框架：把所有场景全部写进一个大优化问题里一次性求解，第一阶段变量是公共的，第二阶段的动态调度变量在各场景内独立。这种方式在10个场景以内时其实更快——省去了迭代通信的反复优化开销，但代价是约束规模大了不少，Gurobi在内存占用和预求解上压力也更大。两者各有取舍，我建议：场景少于15个用全场景耦合，场景更多用迭代分解。

5. 复现过程中绕不过去的坑：调试记录与解决思路

5.1 场景概率归一化问题

一个非常不起眼但后果很严重的坑是场景概率没有归一化。我一开始从蒙特卡洛抽样得到5000个场景，聚类成10个典型场景之后，直接拿每个簇的样本数除以5000作为概率，理论上这没问题。但后来我在抽样时加入了"至少三条线路故障"的筛选条件，导致能进入场景池的场景概率总和小于1，而优化模型里各场景概率之和必须严格等于1，否则最终目标函数会被系统性低估。解决方式很简单：进入优化器之前强制做一次归一化(p_k' = p_k / \sum_k p_k)。

这个坑的危害在于它不报错，所有约束都满足，目标函数数值也对，但算出来的预布局方案在概率意义上是错的。我是在对比不同故障筛选条件之间的结果差异时发现对不上的，追查了半天，最后打印出概率和才暴露。建议你的调试脚本里在求解前至少加一行disp(sum(prob_vector))，保证它恒等于1。

5.2 不可行解的根因定位：锁定阶段一拍错，全盘皆输

两阶段模型一个极其常见的报错是第二阶段求解返回"Infeasible problem"。我遇到过的根因有两类。

第一类是第一阶段的预布局位置在某些故障场景下无法满足任何约束。比如台风导致某片区域的联络线全断了，而MESS恰好被要求预布局在那片区域的孤岛上，但MESS的额定功率又不足以支撑孤岛电压——此时该场景下子问题无解。排查办法是逐场景求解子问题，打印出Io矛盾行，看是哪条电压约束还是功率平衡约束被违反。

第二类是移动时间约束和调度周期不匹配。我模拟的时间尺度是15分钟一个时段，总调度时长8小时，共32个时段。如果从节点10移到节点27需要跨5个时段，而调度周期末尾才安排移动，移动完成已经超出T的边界，模型就会无解。处理办法有两种：一是允许最后一个时段位置继续保持不移动（即终止条件放宽），二是把调度周期延长至少超过最大移动时间加最大充电需求时长。我实测定的是前者，因为延长调度周期会让SOC变量维度变大，求解时间增加不少。

5.3 二阶锥松弛间隙异常偏大的处理经验

正常情况下IEEE33节点的DistFlow二阶锥松弛误差在1e-5量级，完全不用管。但某次我把故障线路的切换状态做成了变量，某场景下部分线路断开、系统分成两个孤岛，这时候某些支路的电流平方变量突然变得异常，松弛间隙飙到0.3以上。仔细分析后发现是断开线路的潮流方程里，(P_{ij})和(Q_{ij})被约束为0，但(L_{ij})没有被强制为0，导致锥约束被"软"化了。解决办法是在线路断开场景下，把断开支路的(P, Q, L)全部固定为0。这个细节如果不处理，得到的负荷削减量会偏低近10%，结论方向不会错，但数值差异会被人审稿人质疑。

5.4 结果可视化：画出一条真正能讲故事的韧性曲线

复现完成后，最值得做的可视化有两张图。第一张是预布局方案图：用plot画出IEEE33的拓扑，用不同颜色标出MESS的预布局节点，配合负荷大小做节点大小映射，一眼就能看出布局逻辑是否符合工程直觉。第二张是韧性曲线：横轴是调度时段，纵轴是系统总负荷满足率（1 - 总失负荷/总负荷），对比"有MESS预布局优化"和"无MESS"两条曲线。优化后的曲线在故障初期的跌落深度明显更小，恢复阶段更陡，这两条曲线之间的面积差正好对应韧性的提升量。

我在画这个图时用了阴影面积填充，在Matlab里用area函数叠两层数据，视觉上非常直观。论文里展示时，再配上各时刻MESS位置变化的甘特图（转移路径），整个方案的时空动态就清清楚楚了。

5.5 求解时间控制：工程可用性和学术探索的平衡

IEEE33节点、10个场景、32个时段、3台MESS，全场景耦合求解在Gurobi下大约要7-15分钟，具体取决于MIPGap的设置和初始可行解的质量。如果只是验证算法"能不能算对"，这个耗时可以接受；但如果想做参数敏感性分析，比如扫MESS数量从1到5、容量从0.5到2.0，那就得在求解时间上下功夫了。我的实践建议是：先用少量场景（比如5个）做粗扫描，锁定了有趋势的参数区间之后，再在关键点用全场景精细求解。这样可以节省至少一半的总计算时间。

6. 从复现到改进：这套模型还能往哪些方向延伸

跑通标准IEEE33节点的复现只是第一步，模型框架本身的可扩展性才是这个方向真正值钱的地方。

最自然的扩展是多时段恢复调度。刚才两阶段模型里的T是单一日内调度，但实际台风恢复过程往往持续数天，MESS需要在多个日调度周期之间持续贡献电量。把单日模型扩展为多日模型，要额外考虑MESS如何回充电站补能、SOC跨日传递怎么建模，以及和抢修队伍排程的协同。这一块我还没完全跑通，但已经开始在做了，预计模型规模会比现在大一个量级。

另一个非常实际的方向是多类型应急资源联合调度。移动储能不是孤零零在战斗——应急发电车、抢修队伍、无人机巡检都是灾后恢复资源。把这些资源统一纳入两阶段框架，让MESS预布局和抢修排程互相配合，数学上更复杂，工程上也更有说服力。这个方向论文发了不少，但真正能跑的代码公开的很少。

7. 写在最后：复现这件事，真正的门槛在模型理解

回到开头那个问题——3台移动储能车，48小时后台风登陆，停在哪儿？跑完这套复现之后，我给自己的答案是：第一台提前布局在馈线中段的重负荷节点，第二台靠近分段开关，第三台作为机动预备队放在转移时间最短的中心节点。这个答案不是拍脑袋拍出来的，而是从上千个场景里优化出来的，每条约束都有出处，每个数字都能追溯。

复现这套策略的收获，最大不在于学会了YALMIP的语法或是Gurobi的参数调试，而是真正理解了"信息时序决定决策结构"这个优化思想——什么时候该做不可逆的决策，什么时候可以等待更多信息再做调整，这套方法论放在配电网韧性里是移动储能预布局，放在供应链里就是安全库存预置，放在灾害应急里就是物资预调配。建模技术是通用的，问题场景只是它的一种投影方式。

如果你也在复现类似的模型，卡在某个约束写不对或者求解器报错，可以沿着我上面总结的排查链路走一趟：先打印概率归一化，再检查阶段变量锁定，再看二阶锥松弛间隙。大概率能帮你少绕几个弯。