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和Pcurt50个场景下来变量数会到几千个好在Yalmip处理这个规模很轻松。3.3 求解流程与结果输出约束装完之后直接调optimizeops 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画节点和支路再用红色方块标出MPSfigure; 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到33t是从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、时段tMPS是否接入节点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移动路径约束和分时功率分配代码补上到时候可以在本文的预配置结果基础上直接接着跑。