这几年做配电网方向的课题韧性这个词几乎绕不开。所谓配电网韧性通俗说就是电网被台风、冰灾这类极端事件打坏以后能不能扛得住、恢复得快。扛得住靠网架恢复快靠抢修资源和应急电源。这篇SCI一区论文复现项目的核心是解决恢复阶段里最现实的一个问题手头有几台应急移动电源MPSMobile Power Source灾害还没发生之前应该把它们停在哪几个节点、预留多大容量才能让后面真正出故障时失负荷最小。项目标题里的“上篇”只做MPS预配置动态调度留给下篇。这篇文章我就用Matlab把模型拆开讲从数学建模、场景生成到求解器调用按实际复现的顺序走一遍既有代码思路也有踩坑记录。1. MPS预配置在配电网韧性研究里到底解决什么问题1.1 先搞清楚时间线灾前、灾中、灾后在做这个复现之前建议先把配电网韧性提升的时间线理清楚否则很容易把预配置和调度两个概念混在一起。通常来说电网应对极端事件分三个阶段。灾前阶段叫预防与准备主要是加固线路、修剪树木、调整运行方式也包括把移动电源提前布置到关键节点。灾中阶段叫抵御与响应故障刚开始发生开关动作、隔离故障、维持剩余网架供电。灾后阶段叫恢复这时候才轮到抢修队伍、修复线路、转移负荷MPS也会在这个阶段被派到故障点附近。预配置解决的问题非常具体灾害来临前资源应该摆在哪里。以台风为例预报中心给出路径和风速分布后调度人员需要决定移动储能车是停在变电站、重要医院门口还是临时放到某个容易孤岛运行的分叉节点。这个决策看起来是“车停哪儿”本质上是一个带有不确定性的资源配置问题。因为直到故障发生后你才知道哪些线路断了、哪些负荷失电但车必须提前停好。这和连锁超市在台风前把发电机、饮用水提前配送到门店是一个道理事后调度再快也弥补不了事前布局的错误。MPS之所以比固定储能更值得研究就是因为“可移动”三个字。固定储能只能服务于固定节点而MPS可以在灾前灵活预置灾后再移动调度实现跨区域的支援。但可移动也带来新的成本需要运输时间、需要接入条件、需要考虑道路能不能通行。所以论文里把MPS分成“预配置”和“动态调度”两个阶段前者决定位置和容量后者决定实际接入节点和出力曲线。1.2 为什么选“两阶段随机规划”而不是“确定性问题”很多人第一次看韧性优化会有一个疑惑直接把故障场景设为“最严重情况”然后优化MPS布局不行吗没那么简单。最严重故障场景往往把所有关键线路同时断开需要配置的MPS数量会非常夸张投资成本完全不可接受。而实际场景是故障线路服从一定的概率分布有的线路可能断了有的没断。我们要做的是在“所有可能故障场景”的期望意义下找一个成本与可靠性之间的平衡点。这正是两阶段随机规划框架的用武之地。第一阶段决策叫做“现在就要定的决策”对应预配置方案哪些候选点安装MPS、配置多大容量。第二阶段决策叫做“看到场景后才做的决策”对应故障发生后每个场景下的失负荷量、MPS出力、潮流分布。目标函数是第一阶段的投资成本加上第二阶段失负荷惩罚的期望值。用数学语言描述就是min 投资成本 运维成本 Σ p_s × Q(x, ξ_s)其中p_s是场景s的概率Q(x, ξ_s)是在预配置决策x给定的条件下场景ξ_s发生时最优的失负荷成本。注意上篇虽然叫“预配置”求解时仍要把每个故障场景下的失负荷结果算出来因为目标函数里包含期望值。这就是很多初学者跑代码时最大的认知误区以为预配置只做灾前优化不需要管灾后潮流。实际上不建模灾后恢复过程预配置就没有评价标准优化出来的方案也是盲目的。从数学性质看这个问题是一个混合整数线性规划MILP或者混合整数二阶锥规划。MIP的好处是能求得全局最优解坏处是大规模场景下求解时间长。所以在复现时场景削减、约束简化、线性化处理都是绕不开的关键动作。1.3 复现之前先想清楚选哪条技术路线我复现这套代码时选的是Matlab YALMIP Gurobi的组合。YALMIP负责把优化问题建模成标准形式Gurobi负责解MIP。如果你没有GurobiMatlab自带的intlinprog也能跑只是规模大了会慢很多建议先把算例缩小再跑。为什么不直接手写单纯形法或者粒子群因为MPS预配置问题本质是离散决策加连续潮流的耦合粒子群这类启发式算法虽然也能解但每次结果不稳定而且很难证明接近最优。既然论文是SCI一区审稿人默认要求全局最优或者有理论保证的求解方法所以复现时最好也用商业求解器不要在算法上自行降级。除了求解器还需要Matlab统计工具箱里的随机数生成、矩阵操作等基础功能。数据算例一般用IEEE 33节点配电网系统这个系统结构简单、数据公开、故障场景容易可视化非常适合做预配置研究。如果后面想发论文或者做毕设再扩展到IEEE 123节点也不迟。2. 预配置优化模型的数学结构与关键约束2.1 目标函数投资成本、运维成本和失负荷惩罚怎么组合复现论文时目标函数是最容易看半天却不知道参数怎么取的部分。我建议先不要纠结具体系数而是把结构理解透。总成本通常写成三项相加min f C_inv C_op C_shed第一项是投资成本C_inv。MPS预配置不是“从仓库调车”而是指在候选节点上预先放置并配置一定容量的移动电源因此要计算设备购置或租赁成本。若候选点i配置了MPS投资成本可以写成单位容量成本C_unit乘以配置容量E_i再乘一个表示是否配置的0-1变量y_i。论文里常见的处理是C_inv Σ c_inv × E_i × y_i也有按“台数×单价”计费的复现时要看原文定义。第二项是运维成本C_op。MPS在故障过程中实际出力会产生燃料、电池损耗等运行成本。它和输出功率P_MPS以及持续运行时间T有关可以写成c_op × Σ P_MPS × Δt。对纯电池储能来说这部分相对小但对柴油发电车来说燃料成本不可忽略。第三项是失负荷惩罚C_shed也是最关键的韧性指标。每个故障场景s下如果部分节点供电不足会产生失负荷功率P_shed配以惩罚系数c_shed再乘以停电时间。最终失负荷惩罚的期望为Σ p_s × c_shed × P_shed × Δt。从数值量级上我采用的一组参考参数是MPS投资成本按1000~1500元/kWh折算运行维护成本取0.1元/kWh失负荷惩罚系数取5000元/kWh左右。惩罚系数设置很关键定太低优化结果倾向于“不装MPS直接停电”定太高则倾向于“到处都装严重超预算”。一般通过试探确定先跑一版结果观察失负荷率是否在可接受范围内再微调惩罚系数。2.2 基本约束安装状态、容量上限与潮流可行性约束条件可以分为三大类预配置决策约束、MPS运行约束、电网运行约束。预配置决策约束最常见的是数量上限和容量耦合。全配电网可选安装点是有限的通常设一个最大可配置数量N_max即Σ y_i ≤ N_max。容量上每个候选点配置的MPS容量不能超过该站点允许的上限E_max同时必须和安装状态耦合E_i ≤ E_max × y_i。这个约束保证了“没安装就没容量安装了才有容量”。MPS运行约束主要限制故障时的出力。例如接入某节点的MPS输出功率不能超过其额定功率P_max且不能超过该节点可消纳的功率上限。如果是储能型MPS还要考虑荷电状态SOC的时序耦合SOC(t) SOC(t-1) η_ch × P_ch × Δt - P_dis × Δt / η_dis以及SOC上下限约束。上篇如果只做静态预配置时间维度可能很简单但为了给下篇动态调度留接口我建议代码里把时间层一并加上。电网运行约束比较麻烦。最基本的节点功率平衡对每个故障场景、每个时段各节点的注入功率、MPS注入功率、负荷消耗功率必须满足潮流方程。配电网通常呈辐射状可以用DistFlow方程描述P_j P_{ij} - r_{ij} × I_{ij}^2 - Σ P_{jk} U_j^2 U_i^2 - 2(r_{ij}P_{ij} x_{ij}Q_{ij}) (r_{ij}^2x_{ij}^2)I_{ij}^2这个方程是非线性的直接放进MILP会爆炸。复现中最常用的办法是线性化DistFlow忽略网损项并把节点电压近似为1.0 p.u.附近的小扰动得到线性化潮流。这样处理在配电网韧性研究中是够用的因为网损占比小真正决定失负荷的是拓扑连通性和有功平衡。还有一类约束是故障场景与拓扑的耦合。场景s中某些线路故障断开网络变成若干孤岛MPS只能在未断开且可接入的节点出力。这个耦合在代码里通常通过预先生成每个场景的供电拓扑矩阵来实现而不是在优化模型里显式列写所有连通性约束。2.3 非线性项怎么线性化才能交给求解器复现过程中最大的“深坑”不是建模本身而是把公式翻译成求解器能接受的形式。Gurobi这类求解器很挑食它只能处理线性约束、二次约束、整数变量等标准形式不能直接接收“y_i乘以E_i”这种双线性项。第一个典型双线性项是投资成本中y_i和E_i相乘。虽然E_i的连续变量与y_i的0-1变量相乘看上去很简单但不能直接写。标准做法是大M法引入辅助连续变量z_i y_i × E_i用三个约束限制z_iz_i ≤ E_max × y_iz_i ≤ E_iz_i ≥ E_i - E_max × (1 - y_i)这样z_i就等价于“如果y_i1则z_iE_i如果y_i0则z_i0”。在YALMIP里可以直接用implies或者手动写约束我习惯手动写因为在不同求解器之间迁移更透明。第二个麻烦是潮流方程中电压平方项U_i^2。线性化时直接设U_i ≈ 1把U_i^2作为变量v_i那么DistFlow方程中v_j v_i - 2(rP xQ)。这样v_i是连续变量约束全部线性。如果论文模型用的是二阶锥松弛则需要写成v_j ≤ v_i - 2(rPxQ) ...的形式并配合锥约束此时求解器必须支持二阶锥Gurobi是支持的intlinprog不支持。第三个容易忽略的地方是失负荷量P_shed不能简单设成一个自由变量。它应该满足0 ≤ P_shed ≤ P_load否则求解器会让 “有功超发” 来平衡目标函数。另外MPS接入的节点只有在拓扑通电时才能出力这一条如果不加优化结果可能会出现“MPS在孤岛外给孤岛负荷供电”的荒谬场景。处理方式是提前生成每个场景的可供电孤岛列表把MPS出力约束到可达节点集合上。3. Matlab代码实现从数据到结果的完整链路3.1 算例准备IEEE 33节点系统的数据组织我复现时用的算例是经典的IEEE 33节点配电网基准电压12.66 kV基准容量100 MVA首端节点1通过变电站供电。这套系统有32条支路、5个联络开关负荷数据到处都能找到但要注意不同文献把负荷标幺值定义得不太一样复现论文前先确认基准值。Matlab代码里我习惯用结构体组织数据不用一堆散落变量。比如定义系统参数表节点编号1~33有功负荷P_load单位kW无功负荷Q_load单位kvar支路起止节点支路电阻r、电抗x单位Ω候选MPS安装点例如节点2、8、15、25、30为什么要保留无功和电抗因为虽然线性化潮流可能忽略网损但无功平衡仍然影响电压约束。有些简化版复现把Q全部忽略导致结果偏差巨大。建议至少把Q_load放进功率平衡方程。候选点选择也有讲究。不能随便选33个点原论文一般会设定MPS只能放在变电站、重要负荷节点或馈线分段开关附近。我用的候选点以“覆盖重要负荷”和“具备停车场与接入条件”为标准选5个左右。如果候选点太多0-1变量数量会剧增求解时间显著变长。故障场景体现在支路状态上。我用一个矩阵branch_status(n_line, n_scenario)来标记每条支路是否故障断开。这样后续潮流约束可以通过矩阵运算批量生成不需要for循环套着写约束速度会快很多。3.2 故障场景抽样与场景削减场景生成这一步直接决定结果质量。设定线路故障概率后最简单的方式是蒙特卡洛抽样。具体做法是对每一条线路生成一个[0,1]均匀随机数如果小于故障概率p_fault就认为该线路断开。为了避免每次运行结果波动要在Matlab代码开头固定随机种子rng(2024);然后生成故障矩阵fault_matrix rand(n_line, n_scenario) p_fault;这里有几个细节。第一线路故障概率可以全局统一取0.05~0.15也可以根据历史台风风速按线路暴露程度差异化设置差异化更贴合实际。第二抽样数建议至少1000个但不能直接拿1000个场景去跑MIP求解器会卡死。所以需要场景削减。场景削减的目标是用少量典型场景逼近原始场景集合的统计特性同时保留每个原始场景的概率权重。业界常用同步回代削减法每次合并两个概率距离最近的场景把其中一个的场景概率加到另一个上直到场景数达到目标值。在Matlab里可以自己写也可以用一些开源函数。我不建议直接调用kmeans因为kmeans会改变场景向量本身而且没有保留概率语义用概率距离去削减更符合随机规划的理论。场景削减函数大概长这样[reduced_scen, reduced_prob] scenario_reduction(fault_matrix, ones(n_scenario,1)/n_scenario, K);削减数量K一般取5~20个。太少丢失极端场景优化结果太乐观太多求解时间直线上升。我复现时用10个场景作为折中点既能覆盖主要故障模式又能在几分钟内完成求解。3.3 YALMIP建模与求解器调用建模阶段我习惯把变量分为三类二进制变量y、连续变量E、场景相关的运行变量P_shed和P_mps。用YALMIP声明变量时让第一维是候选点或节点第二维是场景这样目标函数里可以直接用矩阵运算避免循环。一个简化的YALMIP建模骨架如下y binvar(n_candidate, 1, full); E sdpvar(n_candidate, 1); P_shed sdpvar(n_bus, n_scenario, full); P_mps sdpvar(n_candidate, n_scenario, full); Constraints []; % 最大安装数量 Constraints [Constraints, sum(y) N_max]; % 容量耦合 E_min E_max * 0.1; Constraints [Constraints, E E_max .* y]; Constraints [Constraints, E E_min .* y]; % 失负荷约束 for s 1:n_scenario Constraints [Constraints, 0 P_shed(:,s) P_load]; % 潮流平衡、电压约束等... end Objective sum(c_inv .* E) sum(scenario_prob .* (C_shed * sum(P_shed) * delta_t)); ops sdpsettings(solver, gurobi, mipgap, 0.01, verbose, 2); sol optimize(Constraints, Objective, ops); y_opt value(y); E_opt value(E);有几个关键点需要强调。一是求解器设置里的mipgap论文复现时我可以接受1%的次优性这样求解时间能大幅缩短。二是YALMIP默认可能把目标函数里的E和y相乘当成双线性项必须提前用大M法手动线性化否则会报“bilinear problem”或求解极慢。三是如果只用intlinprog需要把所有sdpvar变量合并成一个决策向量并手动拼接约束矩阵虽然麻烦但能脱离YALMIP依赖。我测下来在5个候选点、10个场景、33节点的算例中Gurobi一般10秒到5分钟能收敛到1%间隙。如果跑半小时还没出结果大概率是约束里出现了我没发现的双线性项或者是场景数太多。3.4 结果可视化与方案输出复现论文不能只给一个“最优成本”数字还要把MPS预配置方案可视化出来。我常用的三张图是第一张是网络拓扑图在节点坐标上用红色五角星标出MPS预配置位置用加粗线标出高故障概率线路。这个图能直观看出MPS是否覆盖了易损区域。节点坐标可以手动写也可以用论文里的相对位置。第二张是成本构成图用柱状图或堆叠图展示投资成本、运维成本、失负荷惩罚在三类场景下的比例。如果失负荷惩罚占比过大说明韧性提升不足如果投资成本占比过大说明MPS配置冗余。第三张是针对每个候选点的容量分布图直观展示“哪个节点装了多少容量”。绘图代码不复杂plot和scatter就能搞定。关键是地图坐标要合理否则画出来像是蜘蛛网。我会先定义节点相对坐标再在图上叠加支路最后用scatter标记MPS站点这样一个图就能撑起论文里一整个结果小节。4. 复现中的坑和调试经验4.1 场景随机性带来的重复性问题头几次跑代码我发现连续两次运行结果不一样而且差异不小。排查后发现问题出在随机数种子没有固定。Matlab的rand每次默认种子不一定相同蒙特卡洛抽样就会变场景削减结果跟着变最后预配置方案自然不一样。解决方法是每次运行前固定rng并把这个种子值写进日志。固定种子之后每次结果可复现调试时就不会怀疑“代码有bug”还是“随机性导致”。另外如果做参数扫描比如改变故障概率从0.05扫到0.15建议不同参数用不同种子否则所有方案都是同一批随机场景的产物缺乏统计意义。4.2 求解时间失控怎么办这个坑几乎每个人都会遇到。原论文用的可能是一个大规模多时段模型复现时如果直接照搬五个候选点加十个场景看似不大但一旦把每个节点每个时间段的变量都建出来整数变量数量会让人崩溃。我的经验是先做“松弛测试”把所有0-1变量设置成[0,1]之间的连续变量看求解时间是多少。如果松弛后的LP都跑得很慢说明约束矩阵本身有问题比如出现了冗余的大M导致的病态数值如果松弛很快、MIP很慢那问题出在整数变量的分支界定上。对付MIP慢的手段包括降低mipgap到0.05或0.1限制最大求解时间优先尝试“固定部分候选点”来缩小整数空间或者减少场景数到5个做敏感性分析。SCI一区论文能跑到0.1%间隙复现代码能跑到1%间隙就已经很有工程价值了不需要追求一模一样的数值。4.3 结果明显不合理时的排查清单我遇到过最无语的结果是优化出来的方案把MPS集中配置在同一个节点其他节点全都不装成本还巨高。排查下来发现是约束条件漏写或者惩罚系数设置不合理。遇到不合理结果时我的排查顺序是先看失负荷惩罚系数是否比投资成本低一个量级如果太低失负荷根本进不了目标函数再看容量耦合约束会不会因为E_max设置太大导致一个点的MPS容量就覆盖了所有失负荷需求然后看潮流约束是不是没有限制MPS只能在通电节点出力导致模型“隔空送电”最后检查场景概率是否归一化如果p_s总和不为1目标函数期望值计算是错误的。这里可以整理成一张速查表方便后面复现的人快速定位异常现象最常见原因排查手段所有MPS集中在一点容量上限太大或惩罚系数失衡检查容量耦合约束限制单点最大容量失负荷量很小但成本巨大惩罚系数太大过度投资降低c_shed观察投资/失负荷的权衡曲线求解器报双线性错误变量乘积未线性化搜索代码中的“.*”乘号改用大M法结果不可复现随机种子未固定在脚本开头加rng固定种子求解时间超长场景数多、整数变量多场景削减、增大mipgap、先跑松弛解4.4 数值尺度的经验Matlab和求解器对数值尺度很敏感。MPS容量单位是kWh动辄几千投资成本单位是万元失负荷成本也很大。变量数值相差六个量级时求解器内部的数值容差会出问题可能出现“所有约束都能过但目标值莫名其妙”的情况。我的做法是统一把功率和容量变成标幺值或相对值。比如以100 kW为基准把负荷、MPS容量、投资成本全部除以基准值让变量集中在0~100之间。目标函数如果用“万元”表示再除以一个系数让整体数值在几十到几百之间。这样Gurobi几乎不需要额外调参速度和稳定性都有明显改善。5. 从预配置到动态调度扩展方向与实用建议5.1 预配置结果如何衔接灾后动态调度上篇的预配置结果不是一个“最终答案”而是下篇动态调度模型的输入边界。灾害发生之后故障场景已经确定系统需要决定每个MPS什么时候从预置位置出发、去哪个故障区域支援、每个时段输出多少功率。这个阶段比预配置更复杂因为加入了MPS的移动时间、道路阻断、抢修队协同等因素。预配置模型输出的每个候选点容量E_i会成为动态调度模型中该点可用MPS的容量上限。如果预配置阶段没在某个区域放MPS灾后再想调过去运输时间可能让失负荷损失进一步扩大。所以在复现时我建议把预配置结果单独保存为MAT文件下篇直接加载避免重新优化一遍。动态调度还要考虑多时段的耦合一个MPS在t时段位于节点At1时段如果转移到节点B需要移动时间约束。这部分典型建模是旅行时间约束可以用弧变量来表示。预配置阶段的候选点数量不多动态调度通常可以把MPS视为“移动储能车队”逐台编号而不是像预配置那样按节点聚合。5.2 工程落地时还需要考虑哪些因素论文算例毕竟是理想化的真要把这套代码用于工程评估有几个因素必须补上。第一个是负荷的时间波动。配电网负荷不是恒定的居民区、商业区、工业区的负荷曲线差异很大。预配置阶段如果只按峰值负荷设计方案会过于保守如果能拿到典型日负荷曲线把预配置评估放到24小时维度会更贴近真实韧性水平。第二个是MPS的交通可达性。移动电源要进入配网节点必须考虑道路是否因灾中断、桥梁是否损坏。更精细的模型会引入图论中的最短路径约束把MPS从一个节点到另一个节点的可行路径和通行时间作为约束条件。这一点在下篇动态调度中最重要。第三个是电池寿命衰减和循环次数。现代移动储能车配备锂电池频繁深度放电会加速老化。如果论文模型不考虑退化成本运维成本会被低估。复现时可以加一个简单的线性退化模型每次充放电积累“吞吐量”超过阈值后运维成本上升。第三个扩展方向是韧性指标的输出。单独看“失负荷期望最小”不够直观建议计算系统韧性提升指标例如配置MPS前后失负荷电量减少多少、关键负荷恢复时间缩短多少、单位投资带来的韧性提升是多少。这些指标对论文和工程项目汇报都非常有价值。5.3 给新手的复现路线图如果之前没接触过这类模型我建议不要一上来就瞄着SCI论文公式抄而是按三个阶段走。第一阶段先搭建一个确定性的单场景模型假设某一条特定线路故障优化MPS预配置目标是在给定预算下让失负荷最小。这个阶段不涉及随机变量代码结构清晰能快速跑通潮流约束和MPS容量耦合。第二阶段加入随机故障场景固定随机种子生成10个场景跑两阶段期望模型。这个阶段重点学场景削减和期望成本的计算方式。第三阶段引入多时段和动态调度把上篇预配置和下篇路径调度串起来最终形成一个完整的韧性提升方案。到这一步论文中的绝大多数图都能复现出来。我给新手的最重要建议是先跑通小算例再追求大算例先接受1%间隙再研究收紧求解器参数先用YALMIP快速建模再考虑是否需要手写intlinprog。这一个流程走下来你对配电网韧性优化、场景法、MILP建模的理解会非常扎实。补充一句个人的实在体会这套预配置模型看起来复杂拆开之后的核心其实是两个变量、三组约束、一次期望求解。“两个变量”是安装状态和容量“三组约束”是投资约束、潮流约束、场景耦合约束。我自己在复现过程中真正耗时最多的反而是场景削减和数值调参而不是模型本身。工具选型上Matlab加YALMIP的优势在于建模快、可视化好、调试方便只要把大M线性化和场景抽样这两个点盯住整体跑通并不难。后面等我整理完下篇再把动态调度部分怎么把MPS路径决策履历出来包括移动时间约束、多周期调度和抢修协同那部分一起聊。