配电网韧性提升：MPS预配置两阶段随机规划及Matlab实现

这个标题一出来做配电网方向的人应该都能get到重点极端事件下移动电源到底该提前放在哪、放多少才能在灾后快速恢复关键负荷同时又不至于白花钱。我当初复现这类SCI一区论文时最大的卡点反而不是那些密密麻麻的公式而是“预配置”这个概念如何落到Matlab代码里——它不像常规配电网优化那样开头就能给出明确的目标函数和约束而是必须先把“灾前决策”和“灾后运行”这两层逻辑在脑子里理顺否则写出来的模型要么是纯投资规划要么就是纯调度问题两头都不靠根本对不起“韧性提升”这四个字。这篇文章我分上下两篇来写。上篇聚焦第一阶段的MPS预配置为什么需要预配置、数学上是怎么建模的、故障场景如何生成与削减、Matlab/Yalmip实现的核心框架以及复现过程中最容易被论文省略、但实际跑代码时必踩的五个大坑。下篇再展开灾后阶段的动态调度包括MPS路径决策、功率时变输出和多个移动电源协调的完整代码实现。这里先说明一下适用人群如果你是刚开始接触配电网韧性的研究生需要有电力系统潮流与运筹学的基本概念如果你已经在做相关方向但卡在了代码细节上这篇文章能帮你把“论文符号”翻译成“能跑的约束”。我会按一条从物理问题到数学模型、再到代码落地的完整链路来讲尽量把每个“为什么这样建模”都交代清楚。1. 从灾区抢修场景理解为什么MPS要先“预配置”1.1 韧性研究的本质是“灾前准备”极端天气——台风、冰冻、暴雨——对配电网的破坏不是慢慢发生的而是一瞬间大面积故障。输配电网的韧性resilience研究关注的正是系统在极端扰动下维持关键负荷供电、快速恢复的能力。传统配电网规划关注N-1可靠性但韧性研究的视角已经拉到N-k甚至整片网络损毁的程度。应急移动电源Mobile Power Source, MPS在这里扮演的角色就像灾害救援中的流动医院平时在仓库里待命一旦主干线路断了、变电站出不了电它可以第一时间移动到关键负荷附近提供临时电力支撑。很多SCI论文都用MPS作为灾后恢复的主要手段因为它灵活、部署周期短、对固定基础设施依赖低而且退役电池梯次利用、移动储能车等工程实践都在推动它的成本持续下降。但MPS有一个天然的约束——数量有限移动需要时间。灾害发生之后你不可能指望从远方临时调集大量电源因为道路可能被毁、通信可能中断、抢修力量分布在很多故障点。这就逼出一个核心问题灾害发生前这些宝贵的MPS应该配置在哪些位置这本质上是一个“不确定环境下的设施选址问题”也就是预配置。1.2 预配置阶段到底在做什么决策预配置阶段需要回答的问题可以拆成三个放在哪位置、放几台数量、每台多大容量。这几项决策在数学上通常全是0-1整数变量因为位置是离散候选点数量是整数台容量也是从几个标准型号里选。这里特别容易混淆的一个点是预配置不是简单的“把电源放在负荷最多的地方”。它需要考虑三件事候选位置本身是否可能在灾害中受损。比如放在地下停车场暴雨一来可能直接被淹放在变电站附近如果变电站是故障源MPS也跟着遭殃。从配置点到故障区域的移动时间受路网、地形、道路损毁影响。台风过后倒树和积水经常阻路移动时间不能只按直线距离估。不同故障场景下同一个预配置方案表现差异巨大。某个位置在场景A里能覆盖医院但在场景B里可能被困住出不来。所以预配置方案必须通过所有可能故障场景的期望表现来评估。这也解释了为什么这个问题的数学模型天然是一个两阶段随机规划而不是一个简单的选址模型。1.3 直接堆在变电站为什么不行很多第一次接触这个问题的人会问为什么不把MPS全放在变电站变电站是配电网的电源点放在那里用起来最方便。这个问题实际跑一次案例就能看得很清楚极端灾害下变电站本身就可能失效——高压侧进线断线、变压器被水淹、保护装置误动MPS全堆在变电站里根本无法输出电力。另外配电网的故障点往往分布在中低压线路沿线的各个区段如果所有MPS都堆在变电站它们灾后要大老远跑向故障现场关键负荷的恢复时间要白白延长好几个小时。预配置的核心权衡是在配置成本买几台、多大容量、占用场地与灾后失负荷损失之间找平衡。论文里的目标函数最常见的写法就是“配置成本 失负荷惩罚”权重设得不同优化结果就指向完全不同的决策风格惩罚系数大方案偏向多买电源、广覆盖惩罚系数小方案偏向少花钱、接受一定程度的失负荷。2. 论文数学模型拆解第一阶段的变量、目标与约束边界2.1 两阶段优化框架随机规划的标准姿势复现这类SCI论文第一步就是把两阶段模型的结构看懂。我按这个方向多篇文献的通行建模框架来梳理下面会注明哪些是常见处理方式。第一阶段预配置阶段在只有灾害预测信息、具体故障未知的情况下决定MPS在哪些候选节点配置、配置多少。第二阶段动态调度阶段在某个故障场景ξ实现之后给定第一阶段的配置方案决定MPS向哪里移动、何时到达、给哪些负荷供多少电目标是最小化该场景下的失负荷量。数学上可以写成min c^T x E_ξ[ Q(x, ξ) ]其中x是第一阶段的预配置决策变量c是配置成本ξ是随机故障场景Q(x, ξ)是第二阶段的最优运行损失一般是失负荷惩罚。E_ξ表示对所有可能场景取期望实际操作中用离散场景加权求和来近似。写代码之前一定要分清哪些变量是first-stage只和位置、数量有关哪些是second-stage和场景、时段有关。如果混在一起写Gurobi会非常痛苦你自己看结果也看不懂。最常见的错误是只建了第一阶段模型、用固定参数代替第二阶段那算出来的预配置方案就是“闭眼拍脑袋”的方案。2.2 目标函数的三个成本块权重怎么配预配置阶段的目标函数一般包含三部分每一部分在Matlab里对应不同的变量成本类型对应变量典型表示MPS预配置投资成本x_place0-1配置变量每台MPS的购置/租赁费 场地改造成本灾后运行成本p_outMPS输出功率变量燃料/电池损耗费 移动成本失负荷惩罚成本p_curt切负荷变量单位失负荷量 × 节点权重 × VOLL这三个部分不是简单相加关键在于权重协调。论文里经常把负荷节点分成一级、二级、三级——医院、水厂、通信基站是一级大型商业是二级普通居民负荷是三级——用不同权重体现“保电优先级”。复现时有个特别容易被忽略的细节失负荷惩罚系数如果设得过大比如大到比配置成本高几个数量级优化器会把所有注意力放在“任何情况下都尽量不失负荷”上结果就是每个候选点都配满MPS方案看似完美实际上已经失去了优化的意义。合理做法是让两类成本在数值上处于同一量级——配置一台MPS的年度化成本大约等于它所能减少的期望失负荷惩罚时才是一个有区分度的优化问题。2.3 约束条件的五层拼图预配置模型的约束大致分五层每一层在代码里都要有对应的方程网络潮流约束用DistFlow配电网辐射状潮流模型近似包含节点有功/无功平衡、电压降落方程、支路电流或容量限制。网络拓扑约束配电网要求辐射状运行故障后允许形成孤岛但孤岛内部也必须保持树状结构不能出现环网。MPS运行约束每台MPS同一时刻只能去一个节点未到达前不能供电输出功率不超过额定容量移动期间不能供电。时间序列约束将灾后恢复过程离散成多个时段MPS移动时间、充电时间、出力上下限都按时序约束。需求约束切负荷量不能超过该节点负荷需求MPS供电量不能超过负荷需求也不能给故障节点供电。这五层约束不是并列关系而是像齿轮一样咬合在一起。特别是潮流约束和MPS运行约束之间的耦合——MPS接入某个节点会改变该节点的注入功率进而影响整条馈线的电压分布这个耦合关系正是“韧性提升”效果的主要来源也是模型最难解的部分。3. 把不确定性装进模型故障场景生成与削减的实操做法3.1 故障场景怎么来的蒙特卡洛抽样是最朴素的做法预配置问题最核心的不确定性是“灾害发生后哪条线路断了”。论文里最常见的做法分三步为每条线路设定一个故障概率根据灾害强度、线路历史故障率、沿线植被/地形条件修正。用蒙特卡洛抽样生成大量故障场景。每个场景是一组“断线集合”。对场景做削减保留少量代表性场景交给优化器求解。具体到Matlab抽样代码其实很简短n_line size(branch, 1); n_scenarios_raw 2000; fault_prob branch.failure_prob; % 每条线路的故障概率 scenario_raw rand(n_scenarios_raw, n_line) fault_prob; % 每个场景是0/1向量1表示该线路故障抽样数量我通常取1000~2000个。太少覆盖不了低概率高损失的极端场景太多后续场景削减耗费时间而且最终保留的代表场景数并不会因此变多。为了让结果更贴近实际部分论文还会把故障概率和灾害强度挂钩台风强度越大每条线路的故障概率越高。一天内的降雨量还会导致多个时段故障概率不同。如果论文研究的是台风场景通常会把“一个灾害过程”作为一个整体抽象成一组故障概率而不是逐小时变化这样计算量更可控。3.2 场景削减别让计算量爆炸直接把2000个场景全部代入两阶段优化Gurobi只能望洋兴叹——二进制变量乘以场景数问题规模立刻爆炸。场景削减就成了必须做的一步。两种最常用的方法K-means聚类法把每个场景看作一个高维向量维数线路数用Matlab自带的kmeans聚成K类每一类用其质心或距离质心最近的场景作为代表场景权重等于该类的样本占比。这种方法简单快速但对极端场景的保留效果一般因为聚类倾向于把包含多条断线的“坏场景”合并掉。快速前向/后向回代消除法这是随机规划文献里的经典算法。从N个场景开始每一步删除一个对概率分布影响最小的场景并把它的概率累加到距离它最近的剩余场景上直到剩下K个。这个算法在论文里通常叫fast forward selection或者backward reduction实现起来就是不断算场景之间的欧氏距离然后找离得最近的场景合并。实操中我的经验是K取5~10左右往往就能在精度和计算量之间取得很好的平衡。你可以先固定一个小规模算例从K3、5、8、12各跑一遍观察目标函数值的变化幅度——如果K8和K12只差不到2%那就没必要用12。削减完的场景要顺手做一步归一化把所有权重加起来归一化为1否则目标函数里的期望计算会偏掉。3.3 模型简化论文不写但你必须知道的取舍复现论文不等于照搬所有细节。在很多SCI论文里MPS移动时间都被简化成只和距离、速度有关忽略路网实际通达性。这样处理在数学上很方便只需要在预处理阶段算一个“候选点到故障节点的最短时间矩阵”。但如果你做的是工程实际项目我建议至少把“道路损毁概率”作为一个调整因子乘上去——台风过后大概率有倒树阻路10公里的路程跑2小时不是段子。另一个常见简化是把故障修复时间设置成定值。实际中修复时间是随机变量而且往往和故障点位置、抢修队配置有关。但第一阶段的预配置决策本质上取决于“哪些区域可能失电、失电多久”所以很多论文用期望修复时间或者最坏情况时长作为近似。这两种简化在场景削减框架下都可以接受你只需要在论文复现报告里注明假设即可不用强行把随机修复时间加进模型——那样会引入大量额外变量对预配置阶段的决策改善非常有限。4. Matlab实现从数据结构到Yalmip建模的核心步骤4.1 建模工具选型为什么是YalmipGurobi这类两阶段MILP问题的标准做法是Matlab Yalmip Gurobi或CPLEX。Yalmip负责建模Gurobi负责求解。核心原因是这类问题最终落在MILP上而Yalmip是目前Matlab生态里把优化模型表达得最接近数学语言的工具代码可读性和可维护性都远高于直接用Gurobi的C接口。如果课题组没有商用求解器的License可以用SCIP或者Mosek兜底。性能会差一些但中小规模算例33节点、123节点一般都能在可接受时间内跑完。Yalmip的核心调用模式非常统一% 定义变量 x binvar(n_candidate, n_type, full); y sdpvar(n_bus, n_period, n_scenario, full); % 定义目标 obj cost_config * x sum(w_scenario .* sum(VOLL .* p_curt, [1,2])); % 添加约束 Constraints []; Constraints [Constraints, sum(x, 2) 1]; % 每个候选点最多放一台 % 求解 ops sdpsettings(solver, gurobi, verbose, 2); optimize(Constraints, obj, ops);4.2 数据结构设计矩阵优先别到处撒变量写这个复现程序数据组织比建模技巧更重要。我第一次复现时把所有参数都定义为零散变量结果代码写到一半自己都忘了mps_vel和mps_velocity_limit哪个才是移动速度。后来老老实实改成结构体问题解决。我常用的数据组织方式% 基础网络数据 net.nb 33; % 节点数 net.nl 32; % 支路数 net.bus bus_table; % 节点表编号、有功/无功负荷、权重 net.branch branch_table; % 支路表首端、末端、电阻、电抗、容量、故障概率 % MPS参数 mps.num 3; % 可用MPS台数 mps.capacity [500 800 1000]; % 各台容量(kWh) mps.speed 30; % 移动速度(km/h) mps.config_cost [10 15 20] * 1e4; % 各台配置成本(元) % 场景数据 sce.prob prob_vector; % 削减后各场景概率 sce.fault fault_matrix; % 场景×线路 的0/1矩阵 % 候选配置点 cand.node [6 11 18 25 33]; % 允许配置MPS的节点 cand.max_per_node 1; % 每个节点最多放几台数据结构以矩阵/表格为主写约束时用索引提取代码既清晰又不容易出错。一个额外的建议把所有参数的单位统一成标幺值或同一量纲不然你会在调试时被101和9.8这样量级迥异的数字搞晕。4.3 第一阶段核心约束的Yalmip写法下面给一个比较完整的预配置模型Yalmip实现骨架以IEEE 33节点系统为例论文里最常用的测试系统。这个代码不是某篇论文的完整原代码而是从复现视角整理出的“满足预配置两阶段建模要求的可运行框架”。变量定义% 预配置决策候选点×MPS类型 的0-1矩阵 x_place binvar(n_cand, n_type, full); % 各场景、各时段MPS出力 p_out sdpvar(n_bus, n_period, n_sce, full); % 各场景、各时段切负荷量 p_curt sdpvar(n_bus, n_period, n_sce, full); % 各场景、各时段MPS是否在某节点供电 z_supply binvar(n_bus, n_period, n_sce, full);目标函数config_cost sum(sum(x_place .* mps.config_cost)); % 预配置投资成本 load_loss 0; for s 1:n_sce load_loss load_loss sce.prob(s) * sum(sum(w_load .* p_curt(:,:,s))); end objective config_cost load_loss;节点功率平衡约束线性化形式Constraints []; for s 1:n_sce for t 1:n_period % 节点注入 常规电源出力 MPS出力 - 负荷 切负荷 Constraints [Constraints, ... P_inj(:,t) p_out(:,t,s) - P_load(:,t) p_curt(:,t,s) 0]; end end这里的P_inj来自DistFlow模型的上游注入功率实际复现时要把DistFlow的支路潮流递推方程一并加进来。为了控制篇幅这里只展示主逻辑。MPS时空可达性约束% 若节点i在场景s时段t有MPS供电则必须存在一个已配置的候选点j % 使得从j到i的最短时间 t * dt for s 1:n_sce for t 1:n_period for i 1:n_bus % reachable(i,:) 是候选点到各节点的可达性0/1向量 Constraints [Constraints, ... z_supply(i,t,s) sum(x_place .* reachable(i,:), all)]; end end end这一条是整个预配置模型里最关键的约束——它把第一阶段的“放哪”和第二阶段的“能否在t时刻到达”联系起来了。如果候选点离故障节点太远超过允许的移动时间z_supply只能为0MPS在该场景下对该节点没有供电能力。这就是时空耦合约束的核心。MPS出力上限约束% MPS在某节点供电时输出功率不超过其容量且不供电时输出为0 for s 1:n_sce for t 1:n_period bigM max(mps.capacity) * 1.1; Constraints [Constraints, ... p_out(:,t,s) bigM * z_supply(:,t,s)]; Constraints [Constraints, ... p_out(:,t,s) 0]; end end求解ops sdpsettings(solver, gurobi, verbose, 2, gurobi.TimeLimit, 3600); sol optimize(Constraints, objective, ops); if sol.problem 0 x_result value(x_place); fprintf(目标函数值: %.2f\n, value(objective)); % 输出配置位置与数量 else disp(求解失败检查约束和参数); end4.4 求解规模控制变量爆炸时的保底操作预配置模型规模瓶颈主要来自二进制变量数量特别是z_supply节点×时段×场景。新手最容易踩的坑就是场景数一多、时段一长变量数直接上万Gurobi直接内存报警。常规破解思路有几种增大时间步长把1小时改成2小时或3小时时段数减少一半变量数直线下降。代价是精度损失但预配置阶段对时间的敏感度通常没那么高可以接受。减少候选配置点不要把所有节点都作为候选。根据负荷重要度和故障概率先筛出前20~30个节点作为候选点问题的整数维度大幅降低。给Gurobi设时间限制和MIP gap比如gurobi.MIPGap, 0.02允许2%的次优解换求解速度。学术复现足够用工程上2%误差也微不足道。这些手段是工程上的“保底操作”论文里不会写但复现的人必须会否则一个问题挂机半天不出解会很影响研究节奏。5. 复现过程中最容易翻车的五个细节5.1 故障场景和辐射状约束的冲突很多配电网测试系统的原始拓扑是辐射状的故障断线后网络会分成多个孤岛。如果某个孤岛既没有变电站供电MPS也没有过去那这个孤岛的负荷就全部丢失——这是模型自然得到的结果没问题。但复现时经常出现这样的错误故障后仍然允许了环网结构。原因是在约束中使用了“线路投切变量”去决定某些支路是否运行而故障线路的断开状态没有被强制写入约束。结果是优化器为了减少失负荷虚拟出一条实际上不存在的连接路径把两个孤岛“缝合”起来算出来的恢复效果虚高。避免方法在场景数据预处理阶段直接把故障线路的连通性置0模型里不要再去“选择”是否断开故障线路。故障是既定事实不是决策变量。5.2 Big-M参数的取值模型里一定有类似P_ij - M*z 0的约束用于在线路停运时强制潮流为0。M如果取太小比如100kW但实际线路容量是5MW约束会错误地限制正常运行时的潮流如果取太大比如1e9MILP的线性松弛边界太松Gurobi分支定界效率极差解算时间暴增数倍甚至数十倍。合理做法是M取支路容量上限的2~3倍。比如支路容量是5MWM取10~15MW。如果想更严谨可以根据电压和线路导纳参数先计算一个理论上限但大多数算例里直接用容量倍数就够了。5.3 时间离散化的一致性预配置模型里时段划分直接影响MPS移动时间的表达。如果MPS从候选点到故障现场需要1.5小时而你用的是2小时一个时段那它最早只能在下下个时段到达——这会错误地扩大失负荷时间预配置方案会偏向于多配MPS来补偿这个“移动时间误差”。解决办法有两个方向一是把时段长度加密到0.5小时或1小时这会同时增加决策变量数量二是干脆把“到达时刻”建模成一个整数变量用到达时间约束替代时间段硬编码。第二种在写代码时复杂一些但求解效率显著更好尤其当MPS数量和节点数都不小的时候。5.4 失负荷惩罚系数的量纲前面提到VOLL不要设得过大这里用一个具体例子说明。假设配置一台MPS的成本是20万元单位失负荷惩罚是5000元/MWh。如果一个方案能在某个场景下减少40MWh失负荷那减少的惩罚成本是20万元正好和一台MPS的成本打平。这个边界条件让优化器真正去权衡“买”还是“不买”结果才合理。如果你的惩罚系数不小心设成了50000元/MWh那减少40MWh失负荷的价值是200万元相当于10台MPS优化器会倾向于“买最多的MPS、覆盖所有节点”算出来的结果没有区分度。论文里的灵敏度分析往往就是在讨论这个系数的变化如何影响配置方案但它对第一次复现的干扰很大——建议先用同一量级跑通再去做系数扫描。5.5 测试系统数据的“版本陷阱”复现常用IEEE 33节点、IEEE 123节点系统但网上的版本极多。有的节点编号从0开始有的从1开始有的支路参数带单位kV、MVA有的已经是标幺值甚至有的版本里支路电阻/电抗的顺序都反了。如果不对数据做归一化潮流约束算出负电压都不奇怪。我的做法是一开始就写一个数据自检脚本三个检查跑完再建模检查支路首末端节点编号是否都在节点表范围内。检查正常状态下全网功率平衡总供电总负荷总损耗标幺值模式下特别明显。检查支路参数是否在合理数量级电阻0.1~10Ω电抗0.1~5Ω容量0.1~10MVA。三个检查都通过再往下写模型。这能帮你省掉后面90%的排查时间别问我怎么知道的。6. 从算例结果反推优化逻辑结果怎么读6.1 典型配置结果长什么样跑通之后第一次看到预配置结果很多人会嘀咕“就这”——MPS往往被配置在少数几个关键节点上而这些节点通常有共同特征负荷量大且等级高医院、数据中心、水厂。处于线路末端离变电站远故障后恢复时间很长。邻近区域线路故障概率较高。把“位置-负荷-故障概率”三张图叠在一起看就能明白优化器的逻辑它本质上是在找一个“最坏情况下损失最小”的位置组合也就是风险最小化的思路。MPS不是放在变电站旁边而是放在“变电站救不到、负荷又重要、故障概率又高”的区域附近。这种“读结果”的能力比跑通代码更值钱。6.2 单场景与多场景方案的巨大差异如果只用某个特定故障场景去优化预配置方案会非常“偏科”它会集中在对应故障线路周边的候选点。一旦实际故障发生在另一个区域这个方案的韧性表现会非常差——可能连一台能用的MPS都调不过去。多场景优化得到的方案则“顾全大局”它不会在任何一个特定场景里做到最优但在所有可能故障场景的平均表现最好。这也是预配置论文一定要用随机规划框架的原因。复现时你可以做一个对比实验把单场景优化得到的配置方案放回全场景集合里评估期望失负荷再和多场景方案做比较。两者差距就是随机规划带来的“鲁棒性收益”也是论文里最常出现的“韧性提升”曲线。6.3 灵敏度分析怎么做预配置结果的工程解释力一般通过几组参数扰动来体现参数扰动的目的MPS数量上限看期望失负荷下降曲线找边际收益递减的饱和点VOLL惩罚系数看配置方案对风险偏好的敏感度反映决策者风格场景削减数量K验证K取值是否合理目标函数是否收敛MPS移动速度看恢复时间对MPS机动性的敏感度为设备选型提供依据每组参数跑完画一条曲线或者整理成表格在论文/报告里就是一组非常有说服力的结果。而且这一步对理解“预配置到底在优化什么”帮助巨大——当你看到MPS数量从1台到2台时期望失负荷大幅下降、从3台到4台时曲线变平缓你就真正理解了为什么这个模型的最优解不是“马配越多越好”。7. 下篇预告与个人体会这篇文章只拆解了MPS预配置的第一阶段。动态调度部分——MPS灾后的路径规划、各时段供电功率决策、多台MPS之间的协调比预配置维度更高、变量更多对求解算法的要求也更高。最后说一点个人体会这类论文的数学符号看再多遍不亲手在Matlab里跑通一遍永远理解不了“预配置”和“调度”的耦合关系到底有多微妙。我复现过程中最大的顿悟时刻是理解到预配置看似只定位置实际上它在预测调度的需求——哪个时刻、哪个节点最需要帮助。时空耦合的思维一旦建立起来后续做动态调度就是水到渠成的事。下篇我会把第二阶段的动态调度模型、路径约束的两种常见建模方式时间扩展图 vs 整数到达时间、以及MPS接力供电的完整案例拉通讲一遍并给出完整的Matlab/Yalmip代码框架方便大家直接对照修改和复现。