从读者视角来看这种SCI论文复现类的项目往往有两个痛点一是论文里公式连篇不知道从哪下手建模二是就算看懂了模型也不知道怎么用Matlab给它落地。这篇博文我打算用一次完整的项目复盘来带大家把MPS预配置这个环节彻底走通。1. 项目背景与核心问题拆解1.1 极端灾害下的配电网韧性到底是什么先聊一个基本概念。配电网韧性英文叫resilience和传统的可靠性reliability不是一回事。可靠性关注的是日常小故障比如某条馈线被施工挖断了、某个变压器过负荷烧了这些属于“常态事件”发生概率相对高、影响范围相对小。而韧性关注的是极端低概率-高影响事件比如台风把整排电杆刮倒、冰灾把导线压断、洪水泡了变电站。这类事件的特点是一来就是大面积停电而且恢复周期可能以天甚至周为单位。在做韧性评估时行业里常用的量化指标包括系统负荷供应能力LSCLoad Supply Capability、系统运行状态转移曲线就是所谓的“韧性三角形”或“韧性梯形”以及各类基于失负荷量的期望值指标。这篇论文的落脚点是通过提前部署应急移动电源Mobile Power SourceMPS来抬高韧性曲线在“灾后响应恢复阶段”的下限说白了就是让关键负荷在电网修复期间尽量不停电或者少停电。MPS这个概念可能有些人不太熟悉。它指的是可以移动的发电/储能单元通常装在拖车或卡车上常见形式包括移动柴油发电机、移动储能车、以及带有并网接口的移动充电宝站。在配电网韧性研究里MPS的优势在于机动性强可以从非受灾区域调往受灾区域灵活接入不同节点形成临时的微电网或孤岛供电。1.2 预配置与动态调度怎么分工论文标题里有“预配置”和“动态调度”两个关键词这两个阶段在时间维度上是错开的。预配置pre-positioning发生在极端事件来临之前属于“灾前准备”。这时候气象部门已经发出预警但配电网还没有实际受损我们需要基于预测信息决定MPS应该提前放到哪些节点每个节点放多少台每台配置多大容量动态调度dynamic dispatch则发生在灾害发生之后属于“灾后响应”。此时的线路断线情况已经逐步暴露我们需要在预配置的基础上实时决定MPS从初始位置往哪里移动、接入哪个节点、给哪些负荷供电以及随着修复进度不断调整供电方案。这两个阶段是典型的序贯决策问题。预配置做得好能给动态调度提供更多的操作空间和更优的初始条件预配置做得差比如把MPS全放在了一个线路大概率会断掉的区域那灾后想调也调不出来。这篇博文聚焦的是“上篇”——MPS预配置也就是把第一阶段的问题建好、解好、代码实现好。动态调度的内容我后面会单独写一篇下篇来展开。1.3 解决这个问题到底难在哪预配置问题看着简单似乎就是一个“把移动电源放到哪里”的选址问题但真正做起来难点有三个。第一个难点是灾前信息的高度不确定性。极端事件发生时哪些线路会断、断几条、断多久都是随机事件我们只能基于气象预测、历史灾害数据做概率推断。这就导致预配置模型天然是一个不确定性优化问题不能当成确定性选址问题来解。第二个难点是配电网本身的运行约束。MPS接入后要考虑节点电压约束、支路电流约束、功率平衡约束还要考虑配电网的辐射状结构约束。配电网和输电网不一样输电网一般是环网结构而配电网通常按辐射状运行这就要求我们在建模时处理好拓扑约束。第三个难点是计算规模。一旦把不确定性场景离散采样出来每个场景都对应一套灾后运行变量模型规模会急剧膨胀尤其是和拓扑重构相关的0-1变量组合爆炸问题非常明显。写代码的时候如果不注意求解效率光等求解器跑完就能把人等崩溃。2. 预配置数学模型构建全解析2.1 目标函数成本还是失负荷还是两个都要我复现这篇论文时最先仔细研究的其实是目标函数。因为目标函数直接决定了优化的导向。常见的MPS预配置论文里目标函数有三种写法第一种是纯成本最小化把MPS的购置/租赁成本、运输成本、运维成本全部折算到目标函数里本质上是一个经济调度问题。第二种是纯韧性指标最大化比如最大化最小失负荷率、最大化灾后小时级负荷供应量这种写法能直接体现韧性提升但在实际工程里缺乏经济性约束项目评审时容易被质疑成本不可行。第三种也是论文采用的写法——两阶段随机优化的期望总成本最小化把预配置的固定成本和灾后各个场景下的失负荷惩罚成本、MPS运行成本统一进一个目标函数。我个人比较认可第三种写法。它的好处是通过失负荷惩罚成本系数来体现“供电可靠性价值”这个系数可以根据不同用户的单位停电损失例如居民、商业、医疗、通信等不同负荷类型差异很大来设定更加符合实际。而且从数学形式上看两阶段随机规划的目标函数天然适配“预配置决策先做、运行调度后做”的信息结构。在具体公式上目标函数可以写成第一阶段总成本 MPS的调用/准备成本 预布置相关的运输成本 E[各个随机场景下的灾后运行成本]其中灾后运行成本包括未满足负荷的惩罚成本、MPS在孤岛供电中的燃料或电池损耗成本甚至还可以考虑线路修复前的临时供电成本。这里有个关键点需要提醒在预配置阶段我们并不知道哪些线路实际会断所以目标函数里带了一个期望算子E这个期望是针对不确定场景求的。代码实现时这个期望往往会离散化为对多个场景的加权平均求和。2.2 约束条件哪些约束是决定模型质量的命门约束条件是配电网韧性优化里最容易写错也最考验功力的部分。我梳理一下论文模型里最重要的几类约束并逐个说明它在代码中怎么体现。第一类是MPS能力约束。每台MPS有额定功率和额定容量预配置到某个节点后它在该节点的可用功率不能超过额定值同时总能量输出不能超过容量。如果MPS是移动储能车还要考虑SOC荷电状态的动态平衡甚至在预配置阶段可以预先充电。如果是移动柴油发电机则更多是燃料约束。论文里的MPS我按移动储能/发电一体机来处理写约束时同时考虑功率上界和能量上界。第二类是潮流约束。配电网运行必须满足基尔霍夫定律。为了规避非凸的交流潮流AC OPF求解困难论文采用了DistFlow分支潮流方程并做了线性化处理。线性化之后电压幅值可以表示为父节点电压减去支路有功/无功与电阻/电抗的线性组合非线性项在正常运行范围内误差很小工程上完全够用。这个处理很关键因为MILP求解器对角色的要求非常高如果保留非凸约束那就沦为一个巨大的NLP求解稳定性会非常差。第三类是配电网辐射状拓扑约束。对于灾后可能形成的微电网孤岛每个孤岛内部必须是辐射状结构。经典做法是利用单商品流single-commodity flow或生成树约束假设虚拟源节点向每个实际节点供应单位虚拟流量每条闭合支路最多贡献一个单位的流动所有节点都得到供应。这组约束在代码里是纯粹的0-1变量线性约束求解器友好度很高。第四类是预配置决策和灾后调度之间的耦合约束。比如某节点在灾后能够由MPS供电前提是在预配置阶段该节点或该区域的接入点被选中且MPS容量足够。这种“先选位置、再定使用”的关系在数学上通过大M法和指示变量来处理。2.3 不确定性建模场景法还是分布鲁棒不确定性建模是这篇论文的另一个技术核心。论文里采用的是基于场景的随机优化也就是预先采样出若干个线路故障场景每个场景给出哪些线路断线、故障修复时间多长然后在这个场景集上做期望优化。场景生成的具体方式我是这样处理的根据台风风圈和线路历史故障概率对不同线路给定不同的断裂概率然后使用蒙特卡洛抽样生成一个包含数百个场景的初始集合。接着用场景缩减技术SCENRED或fast forward selection把场景数量缩减到20~30个这样既能保留不确定性分布特征又能让模型规模控制在可求解范围内。这里有一个我踩过的坑如果直接把几百个场景全塞进MILP模型内存占用会非常夸张求解时间会不可接受。Gurobi面对这类大规模MILP虽然很强势但母问题规模一大Branch and Cut的迭代次数会急剧上升。说到底场景缩减不是可选项而是必选项。除此之外论文提高模型鲁棒性的一个细节是引入“最坏情况”调节因子。通俗地说就是在期望优化基础上对尾部风险极端高风险场景进行惩罚避免优化结果在某个特别恶劣的场景下彻底崩溃。这个思路有点类似于金融里的条件风险价值CVaR优化。如果读者想进一步扩展可以考虑把模型改造成分布鲁棒优化DRO但计算复杂度会显著上升。3. 代码实现MatlabYalmipGurobi落地MPS预配置3.1 环境选型与数据准备我复现这个项目时使用的环境是Matlab R2023a Yalmip工具箱 Gurobi求解器。选择Yalmip而不是直接用Matlab内置的intlinprog最大的原因是Yalmip可以用非常接近数学语言的语法来建模优化变量的声明、约束的批量写入、求解器的参数配置都极其方便。对于这种需要反复调整模型细节的复现工作Yalmip的调试效率比纯脚本高很多。Gurobi之所以不用CPLEX或MOSEK是因为在这个测试系统中Gurobi的MILP求解速度通常最快而且学术许可获取方便。对于中小规模配网测试系统而言即使只用intlinprog也能解但随着节点数和场景数上升还是建议上Gurobi。数据准备方面我采用的是IEEE 33节点配电网测试系统作为算例。这是一个非常经典的辐射状配网有32条支路、5个联络开关基准电压12.66kV总负荷大概3.7MW左右。系统数据在Matpower或文献里都能找到。为了适配韧性研究我对它做了一些改造给每条线路设定一个故障概率模拟台风场景下的故障形态把负荷分成一级负荷医疗、通信权重高、二级负荷商业权重中等、三级负荷居民权重较低三类并赋予不同的失负荷惩罚单价。3.2 变量定义与核心代码逐段拆解在Yalmip中最关键的一步是定义优化变量。对于MPS预配置需要定义三类变量第一类是预配置决策变量这是一个二维的0-1变量表示第r台MPS是否被预配置到节点i。论文里还允许一个节点配置多台MPS所以变量维数是MPS数量乘以节点数量。实际上论文用的方法是连续变量和离散变量混合但为了求解速度和直观理解我们这里先按0-1整数来处理。第二类是功率相关变量包括从变电站上网的功率、各个MPS在预配置节点上的注入功率、负荷是否被切除的0-1变量。第三类是潮流辅助变量包括节点电压幅值的平方、支路电流的平方、支路有功和无功功率。核心代码用Yalmip的写法大致长这样%% 定义变量 X binvar(n_mps, n_bus, full); % MPS预配置决策变量 P_sub sdpvar(n_bus, n_scenario, full); % 变电站注入有功 Q_sub sdpvar(n_bus, n_scenario, full); % 变电站注入无功 P_mps sdpvar(n_mps, n_bus, n_scenario, full); % MPS注入功率 z_line binvar(n_line, n_scenario, full); % 线路状态1闭合0断开 P_flow sdpvar(n_line, n_scenario, full); % 支路有功 Q_flow sdpvar(n_line, n_scenario, full); % 支路无功 U_sq sdpvar(n_bus, n_scenario, full); % 节点电压平方 I_sq sdpvar(n_line, n_scenario, full); % 支路电流平方注意上面用binvar声明的0-1变量有三个分别是MPS预配置位置、线路连通状态。线路状态在预配置阶段是随机参数但在每个场景里是已知的输入。真正需要求解的0-1变量主要是X和z_line以及负荷切除变量。这样声明有一个好处就是约束可以批量写在维度上而不是靠for循环一个个加执行效率高很多。3.3 约束写入与求解器调参约束写入是建模的核心。Yalmip里用很多个约束表达式来拼接约束集。为了便于后期排查问题我会按类别给每个约束命名并且严格区分“所有场景通用约束”和“单场景约束”。对于潮流约束我们的常见写法是Constraints []; % 节点有功平衡 for k 1:n_scenario for j 1:n_bus Constraints [Constraints, P_sub(j,k) sum(P_mps(:,j,k)) - P_load(j,k) P_cut(j,k) sum(P_flow(find(headj),k))]; end end这里head是支路首端节点向量。要注意负荷切除变量P_cut要设成0到P_load之间的连续变量但如果你想让切除行为更符合实际要么全供要么全断就把它设成0-1变量乘负荷量但这会显著增加MIP的难度。经过权衡我在预配置阶段用连续变量松弛因为预配置更关注整体站位精确到单条负荷的供电与否可以留给动态调度阶段以更细粒度启动灵活。电压约束方面DistFlow线性化后的核心关系是U_sq(j,k) U_sq(i,k) - 2*(R_ij*P_flow(l,k) X_ij*Q_flow(l,k));R_ij和X_ij分别是支路电阻和电抗。这个等式在MIQP里被直接写为线性约束。节点电压幅值的上下限则用U_sq的取值范围来表示比如0.95^2到1.05^2之间。MPS接入约束是用来耦合X决策变量和功率变量的关键约束。对于节点j在该节点预配置的MPS台数必须不小于该节点接入的MPS功率数量。实际代码里我用大M法for k 1:n_scenario for j 1:n_bus Constraints [Constraints, sum(P_mps(:,j,k)) sum(X(:,j)) * M_mps]; end end这里M_mps是一个足够大的正数比如取MPS最大功率乘以预配置台数上限。这样当节点j没有被选中预配置任何MPS时右边是0此时P_mps必须为0一旦选中了X容量上限也有了物理意义。求解阶段的核心代码特别简单ops sdpsettings(solver, gurobi, verbose, 2, gurobi.MIPGap, 0.01); optimize(Constraints, Objective, ops);我设置了MIPGap为1%来加速求解。对于单次案例仿真求解时间大约在几十秒到几分钟不等具体要看场景数和MPS数量。如果你是第一次接触这种问题建议先从MPS数量2、场景数5开始试跑确认模型没问题后再把规模加大到论文水平逐步逼近。3.4 结果输出预配置位置分析求解完成后核心的预配置结果可以从X变量里提取。例如假设我们有3台MPS求解得到的X矩阵是一个3乘33的0-1矩阵。统计每列的和就能知道每个节点配置了几台MPS。在IEEE 33节点的经典算例下论文方法的典型结果往往是MPS不会全部放在同一个变电站节点而是分散布置在远离变电站、负载较重、线路故障概率较高的中后段区域。这个结果在直觉上也合理——灾后如果变电站到末端节点的供电线路断了末端微电网可以用MPS组建孤岛供电关键负荷就能保住。预配置位置如果全部放在变电站附近那灾后线路一旦断开MPS根本过不去等于白配。4. 仿真设计、参数设置与结果分析4.1 算例系统说明与参数表在这个部分我把复现时用到的关键参数梳理成了一份可参考的清单。测试系统是IEEE 33节点基准电压12.66kV、基准功率10MVA。线路故障概率的设置不能拍脑袋我采用的是与线路长度正相关的函数同时叠加台风空间分布影响因子。MPS单台额定功率设为150kW储能容量600kWh4小时满功率放电每台MPS的调用成本、布置成本和单位失负荷惩罚成本等参数如下表所示参数项数值备注MPS台数3台可根据算例调整单台额定功率150 kW移动储能车典型参数单台储能容量600 kWh4小时满功率持续供电单位失负荷惩罚系数一级负荷1000 元/kWh医疗/通信等高敏感负荷单位失负荷惩罚系数二级负荷300 元/kWh商业负荷单位失负荷惩罚系数三级负荷50 元/kWh居民负荷MPS预配置固定成本1200 元/台含安装接线、准备费用平均修复时间4小时用于动态调度场景预配置阶段可不涉及MIPGap阈值1%用于Gurobi求解4.2 场景生成与缩减过程场景生成本身就是一个值得单列的项目。我最初直接按300个场景全量建模结果模型里光线路状态0-1变量就有300×329600个再加上其他变量Gurobi的预求解阶段就要花掉好几分钟。后来我改用场景缩减采用了快速前向选择法把300个场景缩减到20个。缩减后的20个场景每个场景携带一个归一化的概率权重。缩减算法过程其实不复杂先定义一个场景间的距离度量论文里常用的是坎托-鲁宾斯坦距离即场景变量和概率分布的加权距离然后迭代地合并距离最近的两个场景直到达到目标场景数。Matlab里可以直接写这个算法也可以用GAMS里的SCENRED但在Matlab里自己写一个百来行的版本并不难。场景缩减完预配置阶段目标函数就写成了总成本 第一阶段固定成本 所有场景概率加权平均的灾后运行成本。由于每个场景的概率权重归一化求和为1期望成本的计算在Yalmip里只需用一个加权求和表达式就能实现。4.3 结果对比预配置与不预配置为了检验预配置策略的有效性我做了三组对照组第一组是“无MPS”基准组即不配置任何移动电源台风灾害后只依赖变电站恢复和线路维修。第二组是“随机放置”组MPS随机分散布置在三个节点但不经过优化模拟没有科学指导的运维决策。第三组是“优化预配置”组使用论文模型求解得到最优位置。三组都评估在20个故障场景下的期望失负荷量。结果非常直观随机放置相比无MPS有一定的失负荷下降但效果有限优化预配置相比随机放置期望失负荷还能再降低大约35%~45%具体数值随故障概率设定不同而波动。这个结果上的差距来源就是优化模型把MPS放在了那些“故障概率高、负荷重要度高、孤岛供电受益面大”的关键节点上。值得一提的是在某个极端场景三条主要馈线同时断线下随机放置组甚至出现了MPS所在地完好无损、而故障区完全无电源可用的窘境。优化预配置组则在同样的极端场景下通过提前把MPS布置在断线高发区域附近成功组建了一个包含核心负荷的微网孤岛一级负荷供应率超过了90%。这一个对比案例非常有说服力建议复现时一定要在论文里配图展示。4.4 灵敏度分析MPS数量与故障概率的影响文章里还做了一组灵敏度分析考察不同MPS数量和不同线路故障概率下的成本变化趋势。整体趋势是MPS数量从1台增加到4台时失负荷惩罚带来的效益改善是边际递减的——第1台MPS带来的失负荷下降幅度最大从优化结果看能降低约20%的期望失负荷惩罚后续每增加一台边际改善逐步减小到第4台时单位成本的收益已经显著变差。这说明在实际工程部署中MPS数量并非越多越好而要结合采购预算和负荷重要性做一个综合权衡。故障概率的灵敏度测试更容易理解随着线路故障概率整体提高预配置MPS的价值会越来越凸显因为更多场景下会出现“变电站供电不可用”的情况MPS在孤岛中的角色权重更大。这个趋势在论文输出里表现为目标函数中对失负荷惩罚项的占比不断上升。5. 常见问题与排查技巧实录5.1 求解时间爆炸怎么定位瓶颈这是复现这类论文时最高频的问题。模型写完一跑五分钟、十分钟没反应先别急着换求解器参数第一步是判断瓶颈在约束规模还是变量整数性。我常用的做法是用Yalmip的size命令查看模型规模model export(Constraints, Objective, ops); size(model.A)如果是行数和列数都到了百万级别那就是场景数和节点数太大导致的“维度灾难”。优先考虑加大场景缩减程度或者把MPS变量由“每台单独0-1”改为“每个节点配置台数整数变量”即把台数上限设为5用整数变量表示而不是用5个0-1变量。后一种做法能把相关0-1变量的数量从n_mps×n_bus降到n_bus级别模型求解难度呈指数级下降。如果是整数变量并不多但耗时依然很长那多半是约束里的大M值设置不合理导致LP松弛过于宽松分支定界效率低下。此时可以把大M值收紧到紧约束范围或者用indicator约束Yalmip的implies函数替代大M法让求解器在预处理阶段就能识别逻辑关系。5.2 模型求解时报不可行怎么排查不可行infeasible是另一大坑。如果solve之后返回infeasible千万别直接去翻代码先用Yalmip的conflict分析功能sol optimize(Constraints, Objective, ops); if sol.problem 1 close find(conflict(Constraints)); end这个功能会高亮与无解冲突相关的约束子集。我的经验里最常见的不可行原因有三类一是负荷需求和变电站供电上限、MPS容量上限之间存在硬冲突尤其当所有MPS都集中在前半段节点而末端负荷无法由任何电源覆盖时功率平衡约束和线路状态约束会打架二是电压等式约束在某个场景下出现迭代结构错误我曾在手写DistFlow方程时把首末端节点搞反了导致等式方向反转结果模型永远无解三是场景数据里有明显异常比如负荷量忘记乘基准值导致单位不统一数值上出现上千倍的偏差。排错的核心思路就是把大模型拆成小模块先只验证潮流约束子模型再逐步加入MPS预配置约束。5.3 结果反直觉优化后MPS反而放在变电站旁边偶尔会有读者问为什么我的模型跑出来MPS预配置结果全放在变电站出口节点这看起来不符合“孤岛供电”的逻辑。这种结果十有八九是目标函数里失负荷惩罚成本设置得太低或者MPS的布置成本太高导致优化器认为“把MPS放在任何地方都收益不大”干脆选择放在成本最低或约束最宽松的节点上。也可能是场景集里线路断线概率设置不合理比如全部线路的故障概率都设得特别小导致最优预配置就是不放或者少放MPS纯靠变电站供电。此时建议把一级负荷的失负荷惩罚系数拉高一个数量级或者加大线路故障概率再观察结果是否改变。另外还有一种情况节点电压约束上限太紧导致MPS放在末端时因线路压降问题无法向关键负荷输送足够的功率模型宁可把MPS放在变电站附近直接调节全网电压。这种情况不是模型错误而是系统物理约束的自然结果。遇到这种情况先检查电压限幅和线路阻抗参数是否符合IEEE 33节点标准再考虑是否需要对末端线路进行加固。5.4 MPS位置结果的工程可行性检查学术模型的解在工程上未必可实施尤其是当你把某个节点选为MPS接入点时需要确认该节点是否具备场地条件、是否有并网接口。我在仿真输出后增加了一个后处理脚本遍历所有选中的预配置节点检查其到道路节点的距离是否小于设定阈值以及该节点是否为变电站或开关站站点。如果结果选在了一个拓扑上极难到达的“孤岛末端”我会在报告中额外提示在第二阶段的动态调度里可能需要额外的运输时间修正。写在最后一点个人复现体会这套MPS预配置模型的复现前前后后花了大概两周时间其中一半时间花在模型调试和数据整理上真正写代码的时间反而没那么长。我的最大体会是论文里的公式和代码实现之间存在大量没有写出来的中间假设。比如DistFlow的线性化到底保留了哪些项、场景缩减的样本数取多少、MPS的容量是否按能量约束处理这些在论文的方法部分往往一带而过但在代码里每一个都对应一个具体的变量或约束。另外想分享一个实用技巧每当拿到一篇新的论文先别急着写代码把论文里的算例参数表格先摘录出来做一个带单位、带含义的参数清单再对照Matpower系统数据逐一核对。很多复现失败的根源就是单位错乱或遗漏参数。一次把参数梳理干净后面的代码实现会顺畅很多。这篇内容只覆盖了MPS的灾前预配置。动态调度部分设计的变量和约束会更复杂会涉及MPS在灾后的路径重规划、多时段耦合的负荷恢复顺序以及和配电网拓扑重构的联合优化。我已经在整理相关内容后续会作为“下篇”单独输出。配电网韧性这个方向学术研究和工程落地之间还有很长一段路要走但每复现一篇扎实的论文都能让你往前走一步。