最近在复现一篇配电网韧性方向的SCI一区论文做的是应急移动电源Mobile Power Source, MPS的预配置和动态调度。这个课题很有意思也很有工程价值台风、冰灾这类极端灾害来临前怎么提前把移动电源放到最合适的位置、配多大的容量灾害发生后再怎么动态调度它们去恢复供电。整个研究分上下两部分这篇先聊上篇——MPS预配置的建模和Matlab实现。我尽量把思路、公式、代码细节和踩过的坑都摊开来讲希望能帮到正在做韧性电网、移动储能、应急电源优化配置方向的同学。先说清楚这篇博文适合谁如果你已经在做配电网优化、韧性评估或者刚接触移动储能调度手头有Matlab和基本的优化求解器YALMIP、Gurobi或CPLEX那这篇文章可以直接当复现笔记用如果你是刚入门的小白也不用慌我会把数学模型里的每个约束都解释一遍再把代码结构拆开讲你跟着搭一遍基本能跑通。1. 项目概述与核心问题1.1 配电网韧性与应急移动电源先交代一下背景。配电网韧性Resilience指电网对极端事件的预防、吸收、适应和快速恢复能力。传统可靠性分析更多考虑随机故障比如线路老化、设备失效但极端天气下的多重故障、大规模停电是另一套逻辑。台风可能同时吹断几十条馈线变电站也可能受灾这时候光靠静态的网架重构是不够的必须有额外的移动资源进场。应急移动电源就是这么个角色。它可以是一台移动柴油发电机、一辆储能车也可以是一台可移动的分布式电源能在灾前部署到关键节点灾后给重要负荷供电。它的核心优势是“灵活”——配电网的网络拓扑是固定的但移动电源可以提前放到预测的故障区域附近等故障隔离后迅速接入恢复供电。这里有个关键点预配置和动态调度是两个时间尺度的问题。预配置发生在灾前决策变量是“每个候选节点放多少台MPS、容量多大”动态调度发生在灾中及灾后决策变量是“什么时候移动到哪个节点、接入哪个负荷、输出多少功率”。上篇只做预配置表示的是“把资源提前布下去”这个动作下篇再考虑动态调度这里先不展开。1.2 为什么需要“预先配置”而不是事后调配可能有同学会问既然MPS是移动的灾害发生后根据实际故障情况再调配不行吗为什么非得在灾前预配置这个问题我在复现时也问过自己细想之后发现逻辑很严密。第一交通网络在灾害中同样可能损毁。台风过后道路受阻、桥梁中断移动电源未必能及时从仓库运到现场。如果灾前提前部署到易受灾区域附近就能规避交通不确定性。第二预配置可以结合气象预测。灾害路径和强度预测通常提前24到72小时发布这段时间足够把移动电源从驻地移动到候选节点。第三从优化角度来说预配置决策是一个“随机规划”问题需要在不确定故障情景下做资源分配目标是最大化工后恢复效果。静态事后调度只能针对已发生的故障做优化但资源数量、位置已经固定可能错过最佳恢复机会。所以论文里通常构建一个两阶段随机优化模型第一阶段是预配置决策第二阶段是故障情景确定后的负荷恢复优化。上篇实现的就是第一阶段。1.3 复现目标与整体思路我复现的目标很明确给定一个配电网拓扑、负荷数据、候选MPS节点集合、灾害故障情景集合求解“在哪些节点配置多少台MPS”这个0-1整数规划问题目标是最小化极端事件下的失负荷期望值同时兼顾经济成本。输出结果包括各候选节点的MPS配置数量、预期失负荷量、关键负荷恢复率以及敏感性分析表。整体思路是先建立物理模型把配电网线性化DistFlow模型再把随机故障情景融入优化然后用Matlab YALMIP建模调用商用求解器求精确解。如果问题规模大可以改用启发式算法但上篇先不追求大规模重点是把数学模型和求解链路跑通。2. MPS预配置数学模型拆解2.1 目标函数经济性还是可靠性优先预配置模型的核心目标不是单纯最小化投资成本也不是单纯最大化恢复负荷而是一个权衡。我这里参考了原论文的思路目标函数包含两部分预配置投资成本和灾后失负荷惩罚。目标函数可以写成[ \min \left( \sum_{i \in \Omega_{MPS}} c_{inv} \cdot x_i \sum_{s \in \Omega_S} \pi_s \cdot \sum_{t} \sum_{j \in \Omega_L} w_j \cdot (P_{j,t}^0 - P_{j,t}^s) \right) ]其中( x_i ) 是节点 ( i ) 处配置的MPS数量非负整数( c_{inv} ) 是单台MPS的日折算投资/租赁成本( \pi_s ) 是故障情景 ( s ) 的发生概率( P_{j,t}^0 ) 是负荷 ( j ) 在时段 ( t ) 的原始有功需求( P_{j,t}^s ) 是情景 ( s ) 下该负荷实际获得的有功功率( w_j ) 是负荷权重反映负荷的重要程度。注意这里没有把MPS运行成本放进第一阶段目标因为运行成本属于第二阶段的调度问题如果在预配置阶段就考虑会让模型变得笨重。从工程角度看预配置阶段的核心矛盾是“投资多少资源”和“预期能减少多少失负荷”所以目标函数里两类成本的量纲要统一。如果投资成本单位是万元/MWh失负荷惩罚单位也必须折算成万元/kWh我当时就是在这里被单位坑了后面细说。2.2 约束条件网络、资源与负荷恢复预配置模型不直接调度MPS出力但需要保证“配置方案在任意故障情景下都有可行的调度空间”。所以约束条件要涵盖配电网运行约束、资源数量约束、负荷恢复约束三块。第一类是配电网网络约束。这里采用线性化的DistFlow方程对于每条支路 ( (i,j) ) 有[ P_{ij} \sum_{k \in N(j)} P_{jk} P_j^L ][ Q_{ij} \sum_{k \in N(j)} Q_{jk} Q_j^L ][ V_j V_i - \frac{r_{ij} P_{ij} x_{ij} Q_{ij}}{V_0} ]其中 ( P_j^L, Q_j^L ) 是节点净注入( V_0 ) 是基准电压。这个线性化忽略了网损和电压降的二次项在配电网辐射状结构下精度是够用的而且能大幅降低求解复杂度。实际中如果节点电压偏差超过5%建议用二阶锥松弛版本不过上篇用线性版本就够了。第二类是MPS资源约束。每个候选节点能配置的MPS数量有上限[ 0 \le x_i \le \bar{x}i, \quad x_i \in \mathbb{Z}{\ge 0} ]还有一个总资源约束比如总共最多可调用的MPS台数或总容量[ \sum_{i \in \Omega_{MPS}} x_i \le \bar{N} ]第三类是负荷恢复约束。在情景 ( s ) 下如果节点 ( j ) 失电可以由MPS供电也可以由电网其他带电区域通过联络开关转供。但预配置阶段不详细建模时序一般只要求[ 0 \le P_{j,t}^s \le P_{j,t}^0 ][ \sum_{j: j \in \Omega_{MPS}} P_{j,t}^s \le \sum_{i: i \in \Omega_{MPS}} x_i \cdot P_{MPS}^{rated} ]这个约束的意思是MPS总容量不能超过预配置资源的总和负荷恢复量不能超过原始需求。有些模型还会加入“一个节点只有配置了MPS才能恢复”的逻辑约束比如[ P_{j,t}^s \le M \cdot \sum_{i \in \Gamma(j)} x_i ]其中 ( \Gamma(j) ) 表示能为节点 ( j ) 供电的MPS候选节点集合( M ) 是一个足够大的数。这个约束非常关键它能防止模型在不配置MPS的情况下“凭空”恢复负荷。2.3 不确定性建模场景生成与概率灾害故障场景怎么来论文通常用蒙特卡洛模拟生成。基本思路是根据历史灾害数据或气象预测给出每条线路的故障概率 ( p_e )然后随机抽样得到一组故障线路集合对应一个故障场景。重复生成 ( N_s ) 个场景就构成场景集合。场景生成这一步影响非常大但很容易被忽略。如果场景数太少预配置方案对真实故障的适应性差场景数太多求解时间爆炸。我在复现时用了100个场景YALMIP建模后求解时间大约在2到5分钟还能接受。如果场景数到500个CPLEX可能半小时都算不完这时就要考虑场景削减技术比如用概率距离快速前推法Fast Forward Selection把场景从500个缩减到50个损失很小但求解效率能提升近10倍。原论文里对不确定性还有更精细的处理比如考虑负荷波动和MPS容量衰退。不过上篇复现时我先假设负荷是确定性曲线故障场景是唯一的随机源这样模型更干净也更容易验证代码正确性。3. Matlab实现思路与代码架构3.1 参数定义与数据准备我的代码结构分成四个部分参数定义、场景生成、模型构建、求解与后处理。参数定义放在一个config.m脚本里方便统一修改。包括网络参数、负荷参数、MPS参数和场景参数。网络参数我直接用IEEE 33节点配电网测试系统这是最常用的算例。节点坐标、支路阻抗、负荷数据网上都能找到就不贴完整数据表了只说明格式bus是节点编号branch是支路起止节点、电阻、电抗load是节点有功/无功需求。MPS参数简化为单台容量500 kW最大配置节点数5个总配置台数上限8台单台日成本根据文献取800元。这里要注意容量和成本都要折算成同一时段的量纲。如果负荷数据是24小时曲线投资成本也要折算成“日成本”否则目标函数里两个量纲不一致优化结果会偏向某一项。场景生成我用了randsample函数。假设每条线路故障概率为0.05生成一个0/1向量表示线路状态然后从故障线路集合中随机选若干条作为同时故障线路。考虑到灾害场景通常伴随多条线路同时故障我设定故障线路数量服从泊松分布期望值5再随机抽取对应数量的线路。代码大致是num_lines size(branch, 1); lambda 5; Ns 100; scenarios zeros(Ns, num_lines); for s 1:Ns n_fault min(poissrnd(lambda), num_lines); fault_lines randperm(num_lines, n_fault); scenarios(s, fault_lines) 1; end注意randperm要求第二个参数不超过总数量所以要用min截断。这个细节我一开始没注意运行时报错“n must be less than or equal to N”查了半天才发现是泊松抽样可能大于线路总数。3.2 优化求解器选择YALMIP CPLEX/GurobiMatlab下建优化模型我强烈推荐 YALMIP。它是建模工具箱可以无缝对接CPLEX、Gurobi、MOSEK等求解器。这个问题的决策变量有整数变量 ( x_i ) 和连续变量 ( P_{j,t}^s )属于混合整数线性规划MILP首选求解器是Gurobi或CPLEX。学术用户可以用Gurobi的免费License性能很给力。YALMIP里用binvar或intvar声明变量。注意MPS数量理论上可以是2台、3台所以要用整数变量intvar而不是二值变量。不过为了简化我也可以把每台MPS作为一个独立的二值变量比如x_{i,k}表示节点 ( i ) 的第 ( k ) 台候选电源是否配置这样模型会自动满足整数性而且后续添加逻辑约束更方便。两种写法本质上等价但后者在YALMIP里更容易处理“一个节点是否配置了任意电源”这类约束。我最终选用了整数变量写法配合value()函数取值。代码核心部分大概是x intvar(length(candidate_nodes), 1); P_rec sdpvar(n_load, Ns, full); cons []; % 资源约束 cons [cons, 0 x max_per_node]; cons [cons, sum(x) total_mps]; % 负荷恢复约束 for s 1:Ns cons [cons, 0 P_rec(:,s) load_demand(:,s)]; cons [cons, sum(P_rec(:,s)) sum(x) * P_mps_rated]; % 关联约束没配置MPS时不能恢复负荷 for j 1:n_load % MPS coverage matrix if cover(j) 0 cons [cons, P_rec(j,s) 0]; end end end注意上面的load_demand我用了确定值但实际上应该根据故障情景判断节点是否失电。如果某个节点本身没故障且与上级电网连通那么即使不配置MPS也能从电网取电。所以正确的建模需要区分“正常供电节点”和“由MPS恢复的节点”否则模型会误伤正常负荷。这个问题我在3.4节详细讲。3.3 网络连通性与失负荷判断这里有个建模难点要确定在故障情景 ( s ) 下哪些负荷是失电的。最直观的方法是对每个场景做连通性分析断开故障线路后从上级变电站松弛节点出发用深度优先搜索或graph对象的conncomp函数找出连通区域。只有与松弛节点相连的节点才能从主网获得电能其余节点只能靠MPS供电。我当时用Matlab自带图函数实现for s 1:Ns G graph(branch(:,1), branch(:,2)); % 删除故障线路 fault_idx find(scenarios(s,:) 1); G rmedge(G, branch(fault_idx,1), branch(fault_idx,2)); bins conncomp(G); % 松弛节点编号 slack bin_slack bins(slack_bus); is_connected (bins bin_slack); % 生成节点失电标识 outage_nodes{s} find(~is_connected); end这个思路是对的但有一个坑graph的节点编号必须连续而且rmedge需要指定端点编号如果你的 branch 数据里节点编号不是从1开始连续排列要先重映射。IEEE 33节点系统编号正好连续所以没问题。连通性分析得到失电节点后就可以设置恢复约束了。只有满足下面条件的节点才能在场景 ( s ) 中恢复负荷节点本身属于MPS覆盖范围候选节点附近一定距离内或者通过与故障区域相邻的联络开关从其他馈线转供。后者涉及网络重构建模复杂我这版先只考虑MPS恢复联络转供留到下篇的动态调度里做。因此对于失电节点如果它不在任何MPS候选节点的覆盖半径内则对应恢复功率强制为0如果在覆盖半径内则可以由配置在该候选节点的MPS供电。3.4 核心模型求解与结果输出把目标函数和约束都放进YALMIP后调用Gurobi求解ops sdpsettings(solver, gurobi, verbose, 2); optimize(cons, objective, ops); x_opt value(x); P_rec_opt value(P_rec);求解时间取决于场景数和整数变量个数。我是33节点系统、5个候选节点、100个场景整数变量5个连续变量3300个约束7000多个Gurobi跑下来大概2分钟。如果场景数再翻倍时间指数上升就需要削减场景。求解完建议马上做合法性检查x是否整数、是否满足总台数限制、总恢复功率是否超过总MPS容量。还要检查每个场景的失负荷是否非负。我当时跑完第一版发现某些场景恢复功率高于该场景失电负荷查了半天发现是目标函数里只惩罚了失负荷没有限制恢复功率不能超过原始需求导致模型“多恢复”以降低惩罚。加上约束P_rec demand后问题就没了。现在总结起来很简单但实际排查花了几个小时。结果输出我建议存成三张表配置方案表候选节点编号、配置台数、总容量场景指标表每个场景的恢复负荷量、失负荷量、恢复率综合评价表期望失负荷、重要负荷恢复率、目标函数值。可视化方面我画了系统拓扑图用颜色标记MPS配置位置用柱状图展示各节点恢复率。Matlab的plot加scatter就够了不用额外工具包。图的价值在于方便写论文汇报但自检时更重要还是看数据表。4. 关键细节与避坑指南4.1 大M约束的正确姿势前面提到关联约束 ( P_{j,t}^s \le M \cdot y_j ) 需要大M但这个M不能随便选。如果M太小可能把可行域切掉如果M太大会导致求解器数值问题出现奇怪的割平面和长时间收敛。经验做法是M取该节点最大负荷的1.1或者2倍只要能保证当 y1 时约束不紧y0 时 P 被迫为0就行。宁小勿大。我一开始图省事设M10000Gurobi求解时间是M100时的3倍多换适当M后快多了。另外如果要表达“某个节点只有配置了MPS才能恢复”大M要乘以该候选节点是否存在配置的变量。如果 ( x_i ) 是连续整数变量直接写 ( P_j \le M \cdot x_i ) 会引入非线性因为 ( x_i ) 是变量。稳妥做法是引入指示变量 ( y_i 1 ) 当且仅当 ( x_i 0 )y binvar(n_candidate, 1); for i 1:n_candidate cons [cons, x(i) max_per_node * y(i)]; cons [cons, x(i) y(i)]; end这样就建立了整数变量和二值变量的等价关系。再用 y 去约束 P 就安全了。这个技巧非常实用很多初学YALMIP的人会卡在这里。4.2 场景削减与求解速度的平衡在复现过程中我对比了三组场景配置50个场景、100个场景、200个场景。50个场景求解约30秒100个场景约2分钟200个场景约8分钟。但结果差异很小因为故障场景重合度高冗余场景太多。想要既保留精度又提速可以用语义丰富的场景削减。YALMIP本身不做削减但可以用Matlab写一个简单的快速前推算法% 简化的场景削减示意 distance pdist2(scenarios, scenarios, squaredeuclidean); selected 1:Ns; while length(selected) N_target % 找与其他场景距离最近的场景 min_dist min(distance(selected, selected), [], 2); [~, idx] min(min_dist); selected(idx) []; end这个算法把聚类中心附近的冗余场景剔除保留代表性场景。注意削减后的概率要重新归一化。削减前后目标函数值差异控制在2%以内完全可接受。4.3 线性潮流模型的选择原论文可能用了二阶锥潮流但我在做预配置复现时用了更简单的线性DistFlow因为预配置阶段不涉及详细的电压无功调整主要关注有功平衡。如果你要对比电压分布建议改用二阶锥松弛YALMIP可以用optimize直接处理SOCP。不过在33节点系统下线性模型的电压误差很小因为馈线负载率不高、电压降落不明显。这里有个个人体会复现论文第一步不要一上来就重建整个复杂模型。先把“预配置 失负荷最小化”这个主干实现跑通后再逐步增加约束比如电压约束、无功约束、MPS移动时间窗。这样每一步都有清晰的验证节点出问题也好定位。我就是一开始想一步到位加上所有细节结果模型写了几百行一跑就是错误最后删掉重来才顺利。4.4 参数量纲与单位的一致性这是一个看似不起眼但杀伤力极大的坑。原论文里投资成本可能是$/台负荷功率是MW失负荷惩罚可能是$/MWh。你如果不做统一优化结果要么全配置MPS要么一台都不配完全失真。我采用的方式是把MPS成本折算为单台套日成本单位元/台·日负荷按日电量kWh统计失负荷价值按元/kWh折算。以某重要负荷为例失负荷惩罚取10元/kWh而一台500kW MPS日成本800元意味着这台电源如果一天能恢复200 kWh以上重要负荷就是划算的。这样把量纲统一后目标函数才有可比性。建议代码里增加一行注释标清楚所有单位否则隔两周自己再看都容易懵。5. 常见问题与排查技巧实录5.1 求解器报错“License”问题YALMIP调用Gurobi时最常见的是License expired或者No license found。学校如果有校园网License还好办个人用户去Gurobi官网申请学术License需要学校域名邮箱。CPLEX也有类似限制。如果实在搞不到商业求解器可以先用开源的HiGHSYALMIP也支持性能比Gurobi差一些但小规模算例完全够用。安装方式是下载HiGHS二进制文件然后设置sdpsettings(solver,highs)。我试过33节点100场景HiGHS也能在三分钟内解出来够复现用。5.2 出现“Infeasible”问题怎么找原因MILP不可行是新手最头疼的事。我的排查顺序是先去掉目标函数只求解约束看是否可行用check(cons)检查每条约束的残差逐步注释掉约束组二分定位不可行约束。最常见的原因有二一是总MPS容量小于某些关键场景的最小恢复需求导致约束矛盾二是关联约束写错导致模型认为“不配置电源也能恢复负荷”或“配置了电源却无法供出功率”。我实际遇到过一种有趣情况某个失电节点既不在MPS覆盖范围又没有联络线路但我把它的恢复功率上限设成了原始负荷值。模型发现无论如何都无法恢复该负荷但又被目标函数引导去恢复它于是只能报不可行。解决办法是把该节点的恢复功率直接强制为0不参与优化。也就是说先做连通性分析把不可达负荷识别出来再用约束固定P0模型立刻可行。5.3 结果与直觉不符为什么配置了MPS但恢复率不高有次跑完结果发现某候选节点配置了三台MPS但该节点覆盖的负荷恢复率只有40%。我检查后发现问题出在故障场景里该节点的联络线路也断了MPS接入后孤岛内只有部分负荷处于覆盖范围其余负荷电气距离太远线路上无法传输足够功率。这说明预配置不能只看节点位置还要考虑局部网架的供电路径。如果候选节点的下游只有一条馈线且该馈线在故障中受损那这台MPS的效用就打折扣。所以后来我加入了一个“有效覆盖范围”预筛选先把故障概率高的线路和重要负荷节点挑出来再逆向选取MPS候选节点确保候选节点的下游网架相对可靠或至少有多条供电路径。这个思路在论文里可能就一句话但实际工程中非常关键。5.4 如何验证模型正确性我在复现时做了几个自检实验大家也可以参考。第一把投资成本设为零此时目标函数退化为最小化失负荷应该会在所有候选节点配置满上限台数。第二把所有故障概率设为零只有一个正常场景那么最优配置应该是0台因为没有灾害就不需要预置电源。第三把总MPS台数设成1而且所有候选节点离重要负荷都非常远此时结果应该是不配置任何MPS因为配置了也无法恢复负荷只会增加成本。这几种极端情况能快速发现模型约束写错或参数设置异常。我用这套自检流程抓到过两个隐蔽bug一个是场景概率没归一化导致期望失负荷偏大另一个是P_rec的维度写反导致负荷恢复张冠李戴。6. 上篇实现中的个人体会到这里MPS预配置的核心模型和Matlab实现就讲得差不多了。整个过程走下来我最大的体会是复现一篇论文真正花时间的不是抄代码而是理解每个约束为什么这么写。比如大M约束、场景生成、连通性判断这些细节如果照搬论文而不理解遇到问题完全无从下手。另外一个很实用的经验是先跑通一个小规模算例再逐步扩大。我一开始直接上了IEEE 123节点系统结果模型几十万个变量求解器跑了一晚上都没收敛差点放弃。后来回到33节点系统把逻辑理顺再重新做大系统就顺手多了。如果你也在复现强烈建议从小算例开始。关于下篇的动态调度我已经在构思了。动态调度比预配置更复杂因为要引入时间维度和MPS的移动路径约束目标也不再是简单的期望失负荷最小而是多时段的切负荷和恢复顺序优化。等我把下篇复现完会继续写一篇关于MPS动态调度的Matlab实现笔记包括移动时间窗、交通约束、联络开关重构等内容。这篇先到这儿如果你对某一部分有疑问欢迎在评论区留言我尽量针对实际问题回复。