☰
移动储能预布局提升配电网韧性:MATLAB+YALMIP建模与求解实战
2026/10/6 14:24:59 网站建设 项目流程

前几天帮实验室复现移动储能预布局方向的论文,卡了好几个晚上才把代码跑通。说实话,这类文章学术价值很高,但作者往往把精力花在理论上,程序实现里那些坑——约束怎么排、求解器为什么报不可行、结果怎么验证——论文里一个都不会写。今天把整个过程整理出来,包括我自己踩过的坑和调试思路,给准备做类似方向的同学一份能直接上手的参考。

先说清楚这篇东西覆盖什么:移动储能(Mobile Energy Storage,MES)预布局提升配电网韧性问题的建模、求解框架和 MATLAB 实现路径。参考的那篇《面向配电网韧性提升的移动储能预布局...》里面,核心讨论的是极端灾害来临前怎么提前安排移动储能车的位置,灾害发生后怎么调度它们去恢复关键负荷。说白了就是一个“先布点、再调度”的两阶段决策问题。如果你正在做配电网韧性、防灾应急、储能调度相关的研究,或者导师丢给你一个带“移动储能”关键词的项目,这篇文章应该能帮你省下不少试错时间。

1. 跨过第一道坎:搞清楚“预布局”比“事后调度”难在哪

1.1 韧性、可靠性、安全性:三个概念别混着用

做这个方向之前,我一直把配电网的可靠性和韧性混为一谈,直到被导师纠了一次才彻底分清。可靠性(Reliability)关注的是正常状态下系统持续供电的能力,比如 N-1 准则、平均停电时间(SAIDI)、平均停电频率(SAIFI)这些老牌指标,它们描述的是“电网平时稳不稳”。韧性(Resilience)描述的是另一种能力:面对极端小概率事件时,系统能不能扛住冲击、在故障后能不能快速恢复,它关注的不是“常态”,而是“非常态”。

这么一区分,问题就清楚了。台风来了,线路断了一堆,你不可能靠日常的可靠性措施完全预防;你只能让关键负荷(医院、通信基站、水厂)在馈线失电后依然有电,或者尽可能快速恢复供电。韧性提升的核心是“灾前布局 + 灾后响应”的组合拳,移动储能恰好在这两个环节都能派上用场:灾前给它定好位置,灾中它可以脱离电网独立供电,灾后它又能在极端条件下作为黑启动电源或者临时电源。

1.2 固定储能与移动储能的本质差异:决策空间多出一个“位置”维度

固定储能(BESS)的特点是容量大、响应快、安装位置固定,它的决策变量只有“何时充电、何时放电、发多少功率”。移动储能(MES)通常安装在车辆或者可搬运的储能单元上,除了充放电状态,你还得决定“往哪走、什么时候走、到了之后接入哪个节点”。

别小看这个“移动”维度。固定储能相当于你把资源放在那,赌灾害发生的位置正好在你的覆盖范围内;移动储能则允许你基于灾前预报信息调整资源位置,不确定性被显式处理进决策里。

这也解释了为什么移动储能的建模复杂性远高于固定储能。因为充放电变量还是连续变量,但“位置选择”和“移动路径”是离散变量,最后模型会变成一个混合整数线性规划(MILP),求解难度随节点数和场景数成倍增长。很多论文里常用的 IEEE 33 节点、123 节点配电网系统,运行一次到最优解往往需要几分钟,大系统可能要几十分钟甚至更久。后面我会详细讲怎么解决这个效率问题。

1.3 预布局问题的两阶段性质:灾前不知道、灾后才知晓

“预布局”之所以难,本质上是信息不对称造成的。灾害发生前,我们只能根据气象预报或历史灾损模型估计哪些区域可能受损,但具体到某条线路断还是不断、某个节点是不是失电,都是不确定的。你必须在不确定条件下做位置决策,然后在不确定性实现后再做调度决策。

在数学上,这个结构对应的是两阶段优化:第一阶段决策(预布局位置)必须在不确定性揭晓之前做出;第二阶段决策(实际充放电功率、部分移动资源重新调动)则是在灾害场景明确之后做的适应性调整。经典的处理办法有两种:

  • 随机规划(Stochastic Programming):给每个场景设定一个概率,目标是最小化所有场景下的期望代价。
  • 鲁棒优化(Robust Optimization):设定一个不确定集合,目标是最小化“最坏情况下的代价”。

参考的那篇文献主要走的是鲁棒优化路线,原因是极端灾害场景样本稀缺、概率难以准确估计,用区间不确定集合描述“哪些节点可能失电”更贴合实际。我在复现时也优先采用了鲁棒框架。

2. 核心数学模型:从韧性指标到混合整数规划

2.1 目标函数:最小化失负荷代价加上移动储能使用成本

目标函数怎么定,直接决定了模型能不能反映“韧性提升”这个诉求。我见过不少论文直接以“负荷削减总量最小”为目标,但实际写代码时你会发现这样设定太粗糙——它区分不了一个医院和一个普通居民负荷谁更重要,也不体现恢复速度。

我更推荐的建模方式是分层量化:

  1. 给负荷分类:一类负荷(重要负荷)、二类负荷(较重要)、三类负荷(一般负荷),分别给不同的权重(比如 100、10、1),这样求解器会优先恢复一类负荷。
  2. 目标函数包含所有时段内的负荷削减量乘以权重,再累加。
  3. 对于“韧性提升”,引入恢复时间相关指标——同样是削减同样多的电量,早恢复和晚恢复在韧性曲线上体现出来的价值完全不同。

这里可以用一个简化的表达式描述目标函数:

$$ \min \sum_{t \in T} \sum_{i \in N} w_i \cdot P^{cut}{i,t} \cdot \Delta t + \sum{m \in M} c_m \cdot \text{moveCost}_m $$

其中第一项是加权失负荷代价,第二项是移动储能的移动成本。实际代码里,我们还经常把“灾后恢复的负荷曲线下面积(AUC)”作为韧性指标来验证模型效果,后面会详细讲。

2.2 配电网潮流与节点失负荷约束:DistFlow 是合适的选择

配电网是典型的辐射状网络,直接解交流潮流(AC Power Flow)在 MILP 里不现实,绝大多数相关研究都采用线性化的 DistFlow 分支潮流模型。

DistFlow 的核心思想是把支路潮流写成功率流动逐步递推的形式。简化之后,每条支路的有功、无功、电压关系可以写成以下几组约束(线性化版本):

  • 节点功率平衡:流入节点的功率 + 储能发出的功率 = 负荷消耗功率 + 节点注入功率 - 失负荷。
  • 支路电压关系:节点 j 的电压幅值平方 ≈ 节点 i 的电压幅值平方 - 2(RP+XQ),其中 R 和 X 是支路阻抗。
  • 支路容量约束:流过的视在功率不超过上限。

在 MATLAB + YALMIP 框架下,这些约束直接写成矩阵乘法形式的等式和不等式即可。我建议把网络参数提前存成邻接矩阵形式,这样写约束循环时思路非常清晰。那些直接在目标函数里写线性潮流的人十有八九最后会吃电压越限的亏——母线电压压降超过允许范围,结果曲线很漂亮但物理上完全不合理。

2.3 移动储能的时空耦合约束:最容易出错的一段

这一节是整个模型的重头戏。移动储能的约束大致分三类:运行状态约束、容量约束、空间移动约束。分开说。

运行状态约束定义储能在一个时段内要么充电要么放电或者闲置,不能同时充放电。一般情况下:

$$ 0 \le P^{ch}{m,t} \le B^{ch}{m,t} \cdot P^{rate}, \quad 0 \le P^{dis}{m,t} \le B^{dis}{m,t} \cdot P^{rate} $$

其中 B 是二进制变量,且 B_ch + B_dis ≤ 1。这个约束很多初学者会漏掉,结果求解器算出一个“又充又放”的诡异解,还白白消耗储能寿命。

容量约束描述的是荷电状态(SOC)随时间的变化:

$$ SOC_{m,t} = SOC_{m,t-1} + \eta_{ch} P^{ch}{m,t} \Delta t - \frac{P^{dis}{m,t}}{\eta_{dis}} \Delta t - E^{move}_{m,t} $$

注意最后一项:移动本身要消耗一部分电量。很多论文会忽略移动消耗,代码跑出来恢复效果当然“更好”,但实际车辆移动是有能量损耗的,建议保留这一项,模型会更可信。

空间移动约束是移动储能最有特色、也最容易出错的地方。移动储能车在 t 时段处于节点 i 的二进制变量为 1 时,它才能向该节点放电或从该节点充电;移动时间需要显式建模:

$$ \sum_{m} \text{loc}_{m,i,t} = 1 \quad \forall i,t? \quad \text{(上式仅在限定时段适用)} $$

实际写代码时,我们通常用一个小的时间滞后矩阵来表示“从节点 i 转移到节点 j 需要几个时段”。我最初就是因为没考虑移动时间,让储能车“瞬移”,结果求解器给出的调度计划在物理上根本不可能实现。具体怎么修改约束,在第 4 节我会专门讲这个坑。

2.4 韧性指标量化:光有目标函数还不够

目标函数解决“怎么优化”,韧性指标解决“怎么评价、怎么对比”。做实验的时候,如果只有目标函数值,审稿人或导师会问:你的方案比固定储能方案好在哪里?这就需要画韧性曲线。

通常的韧性曲线是横轴为时间(灾前、灾中、灾后恢复),纵轴为系统性能水平(可理解为负荷供给率、供电能力等)。一次灾害事件过后,曲线先快速跌落,然后随着抢修和储能支援逐渐恢复。韧性指标可以用曲线下面积(Area Under Curve)衡量:

$$ R = \int_{t_0}^{t_f} \text{performance}(t) , dt $$

值越大,说明灾害过程中整体损失越小、恢复越好。我在实验里通常把预布局策略、固定储能策略、无储能策略三条韧性曲线放在同一张图里对比,直观说明移动储能的优势。

3. 求解框架与 MATLAB 代码设计:让问题落地

3.1 场景集合与鲁棒模型的处理:从“不确定”到“可求解”

移动储能预布局问题最大的难点是求解复杂度。如果直接枚举所有故障场景,组合数爆炸,任何商业求解器都扛不住。参考论文里使用的是场景集合与鲁棒优化结合的方式,具体做法是:

  • 不确定集合定义为“可能受到灾害影响的线路/节点集合”。比如台风路径预测里可能受损的 10 条线路就是不确定集合的成员。
  • 外层问题(主问题)决策预布局位置;内层问题(子问题)在给定布局下求解最坏场景的灾后调度。通过逐次生成割平面(C&CG 算法)反复迭代直到收敛。

这种方法比单纯的随机期望优化稳健,也比全枚举高效。实现时不需要自己写底层的 Benders/C&CG 代码,YALMIP 可以配合求解器直接求解紧凑重构形式,或者用 bait 风格手动实现主问题-子问题迭代。

3.2 YALMIP 建模的代码结构实践

我复现时采用的 MATLAB 代码结构大致如下:

  1. load_case.m:定义配电网拓扑、负荷曲线、线路参数。我用的 IEEE 33 节点系统,节点负荷数据来自论文附录。
  2. scenario_generator.m:生成灾害场景集。做法是随机抽取不同线路故障组合,输出受影响节点集合。
  3. build_model.m:调用 YALMIP 定义变量、写约束、设定目标函数。
  4. solve_case.m:调用 Gurobi/CPLEX 求解 MILP。
  5. plot_results.m:距离曲线、韧性曲线、SOC 变化图、储能位置图。

以 YALMIP 代码片段为例,核心建模思路是这样:

% 定义变量 P_dis = sdpvar(n_mes, n_node, n_time, 'full'); % 放电功率 P_ch = sdpvar(n_mes, n_node, n_time, 'full'); % 充电功率 SOC = sdpvar(n_mes, n_node, n_time, 'full'); % 荷电状态 B_dis = binvar(n_mes, n_node, n_time); % 放电状态 B_ch = binvar(n_mes, n_node, n_time); % 充电状态 loc = binvar(n_mes, n_node, n_time); % 位置变量 % 约束:一个移动储能同一时刻只能在一个节点 for t = 1:n_time Constraints = [Constraints, sum(loc(:,:,t), 2) == 1]; end % 约束:充放电不能同时进行 Constraints = [Constraints, B_dis + B_ch <= 1]; Constraints = [Constraints, P_dis <= M .* B_dis]; Constraints = [Constraints, P_ch <= M .* B_ch]; % 目标函数 Objective = sum(w_cut .* P_cut(:)) + sum(alpha .* move_dist(:));

这里注意:大 M 的取值不是随意写的,它会直接影响求解器的数值稳定性,这一块也是我踩过的重灾区,具体放第 4 节。

3.3 求解器选型与参数配置

对于 MILP 问题,我试过三种搭配:

求解器优势劣势适用场景
Gurobi对 MILP 支持最好,并行性能强,默认参数就能跑得不错商业授权,学术版免费但有时限推荐首选
CPLEX老牌求解器,稳定性好在新版本中部分 MIP 策略不如 Gurobi 激进备选
SCIP开源,可自由使用速度明显落后于商业求解器小规模测试、课程作业

我用 Gurobi 跑 IEEE 33 节点、24 时段、5 辆移动储能车的模型,默认参数下大约 30 秒能收敛到 1% 的 MIP gap;如果模型规模翻倍,时间会指数增长。建议调用的求解器参数如下:

ops = sdpsettings('solver', 'gurobi', 'verbose', 2); ops.gurobi.MIPGap = 0.01; % 设置 1% 的精度门槛,够用了 ops.gurobi.TimeLimit = 3600; % 防止长时间卡死 ops.gurobi.Threads = 8; % 多核并行

提示:不要一开始就追求 MIPGap = 0,那会让求解时间暴涨,而且工程上 1%-5% 的 gap 完全足够支撑结论。

3.4 初始解的生成与热启动

还有一个容易被忽略的操作:热启动。MILP 求解器如果从一个较好的可行解开始,分支定界树的剪枝效率会大幅提升。

我采用的做法是:先忽略移动时序约束(允许储能车“瞬移”),快速求一个松弛解作为初始可行解,再传给 Gurobi 做 warm start。实测在部分场景下能把总求解时间缩短 20%-30%,代价只是预处理阶段多花十来秒。

% 热启动示例 assign(P_dis, P_dis_initial); assign(B_dis, B_dis_initial); ops.gurobi.MIPStart = 1;

这种思路对于大型系统非常实用,研究过程中没必要每次都从零开始暴力搜索。

4. 复现中的关键细节与踩坑记录:从报错到可重复实验

4.1 大 M 值不是拍脑袋定的:数值病态与不可行解的根源

这是我第一次跑通代码时最扎心的教训。在移动储能建模里,大 M 通常用来表示“如果不在这个节点,就不能充放电”之类的逻辑约束。起初我图省事,M 直接设为 1e6,结果求解器给出的结果要么不可行,要么迭代半天不收敛。

原因是过大的 M 值会让线性松弛后的可行域过于宽松,分支定界效率大幅下降,同时可能导致数值精度问题。正确的做法是:M 取该储能额定功率或容量上限的 1.05-1.1 倍即可。比如储能额定功率 500 kW,M 就取 550,而不是 1000000。

% 推荐:用储能参数推导 M,而不是写死 M_power = 1.1 * max(P_rate(:));

4.2 移动时间约束缺失导致“瞬移”问题

前面提到的空间移动约束,是我调试时间最长的地方。忽略移动耗时的模型跑出来的最优策略会在同一个时段内出现在两个不同节点,显然不现实。

解决方法是引入移动状态矩阵:储能车从节点 i 转移到节点 j 需要 d_ij 个时段(通常取两者最短路径长度除以车速取整)。在这 d 个时段内,储能车既不能充电也不能放电。我在构建模型时增加了一组“转场禁止”约束:

% loc(i,t)=1 表示储能车在 t 时刻位于节点 i % move(i,j,t)=1 表示 t 时刻开始从 i 转移到 j for t = 1:n_time if t + travel_time(i,j) <= n_time Constraints = [Constraints, ... loc(:, t+travel_time(i,j)-1) >= move(i,j,:,t)]; end end

类似地,移动期间不充放电的约束就是强制该时段功率为 0。

4.3 求解时间爆炸:场景削减与分解迭代

如果是期末展示或论文复现,你可以一次性跑一两个小时等结果,但如果你要调参数、比较多个方案,求解时间直接决定你的工作效率。我在实验中尝试过以下几种缓解策略,效果排序如下:

  1. 场景削减(Scenario Reduction):先用 K-means 或快速前向选择法把原始数百个灾害场景聚合成 10-20 个代表性场景,目标函数改成期望损失,求解速度会快一个数量级。
  2. C&CG 外层主问题-子问题迭代:不需要一次性列全部场景,每轮只加入当前最坏场景,逐步逼近鲁棒最优解。
  3. 固定部分整数变量:先固定位置变量求连续功率,再用启发式方法微调位置,属于次优但实用的工程解法。

这三种方法我按顺序全部实现了。场景削减写起来最快,适合做初版实验;C&CG 结构最完整,适合作为论文的核心算法;固定整数变量启发式法适合做对照实验或快速估算。

4.4 结果验证:怎么证明你的代码算对了

这个坑我必须要讲:很多同学代码跑出结果就以为大功告成,但其实不同设置下跑出的目标函数值完全不可比。

我常用的验证方法有三步:

  • 设定一个“无储能”的基线场景,跑一遍得到韧性曲线,确认失负荷量大于所有有储能场景,如果比有储能场景还小,说明代码有 bug。
  • 把移动储能的行驶速度参数设成无穷大(或移动耗时设为零),看结果是否退化为“移动储能全知全知”的理想上界。如果退化后的结果和论文中列出的理想值不一致,说明约束有遗漏。
  • 将每个场景下的最优解代入原始潮流方程做一次物理校验,逐节点检查电压、功率是否越限。由于我们用的是线性化潮流,误差通常在 3%-5% 以内,如果偏差过大就要怀疑确定性模型参数是否出错。

注意:线性化潮流的误差会随负载率上升而放大,重载场景下最好引入电压约束的凸包络修正,或者直接改用二阶锥松弛(SOCP)形式,YALMIP 里只需几行改动就能切换。

4.5 参数敏感性分析:别让结论建立在运气上

等到模型跑通,我建议不要急着下结论,先把关键参数的敏感性分析做一遍。至少这几个维度值得测:

  • 移动储能数量从 1 到 10 递增,韧性指标的变化斜率——斜率放缓说明资源饱和。
  • 储能容量从 200 kWh 到 2000 kWh 递增,看是“数量”还是“单机容量”更影响恢复效果。
  • 不同灾害场景集合的鲁棒优化和随机规划结果对比——确认你的结论在不同不确定建模方法下不完全相反。

我在做敏感性分析时比较意外的一个发现是:在中小规模系统中,移动储能的“灵活性收益”主要来自布局阶段的策略性位置选择,而调度阶段的路径优化对韧性指标的边际贡献反而没那么大。这提醒我,花太多精力去优化灾后路径细节,可能不如仔细研究灾前场景生成更有效。

5. 从复现到改造:代码怎么扩展成自己的论文方案

5.1 从单目标到多目标:加入经济性维度

参考论文的核心优化目标是韧性提升,但你投稿时如果完全沿用一个目标,评审很可能说创新性不足。一个比较容易加入的扩展方向是把“韧性提升”和“全寿命周期成本”放在同一个框架里权衡。

具体实现时,把移动储能购置/租赁成本、运维成本、移动成本全部折算成年值,然后在目标函数里加一个权重系数 λ:

$$ \min \left( \sum_{t}\sum_{i} w_i P^{cut}{i,t} \Delta t \right) + \lambda \cdot \text{Cost}{total} $$

用参数扫描画出 “代价-韧性”帕累托前沿,这类结果在电网公司实际规划中说服力非常强。

5.2 从理想化到工程化:柴油发电机与移动储能的联合调度

现实中,移动储能车往往不是现场唯一可用资源,应急柴油发电机、可调度负荷、分布式光伏都可能参与恢复。如果你想做更有工程价值的扩展,把移动储能和应急柴油发电机放在同一个调度框架里是非常自然的切入点。

这个扩展的建模改动不大:柴油机引入燃料约束、碳排放约束和启动成本约束,目标函数相应增加燃料和碳排放惩罚。移动储能负责快速响应和零排放,柴油机负责支撑大功率长时间负荷,两者互补。

5.3 从离线到在线:滚动时域控制(RHC)

另一个更贴近实际的方向是滚动时域控制。预布局之后,灾害过程中信息不断更新(某条线路实际没断、某节点负荷实际比预测小),固定一次求解的结果就过时了。我在扩展实验里用多阶段滚动时域框架,每个时段重新求解未来 T 个时段的子问题,实测恢复效果比单次优化好很多。

for t = 1:n_time-T_win+1 window_data = load_data(t : t+T_win-1); x_opt = solve_window_model(window_data, current_SOC); apply_commands(x_opt(:, 1)); % 只执行第一个时段的指令 end

这种在线决策模式虽然概念上是把多阶段优化拆成多个单次优化,但文章的故事性会强很多,非常适合作为论文创新点。

6. 最后的实操建议:给正在跑仿真的人

如果你现在正卡在“移动储能 + 配电网韧性”仿真上下不来,我有几点心得可以直接抄:

  • 先跑通小算例再放大。不要一上来就上 123 节点、几百个场景,先用 6 节点或 33 节点、3 个场景、1 辆储能车把模型逻辑调通,确认韧性曲线形态合理再逐步扩大。否则连不可行原因都查不清。
  • 数据尽量用公开算例。配电网常用 IEEE 33/123 节点系统,负荷数据、线路阻抗都能直接找到,避免因为编造数据浪费大量时间。
  • 把调试信息打开。YALMIP 的yalmiptime和求解器的 log 日志要养成看的习惯。如果约束数量和解算时间明显异常,问题往往出在循环写约束时维度匹配出错。
  • 版本控制别偷懒。我改用 Git 管理 MATLAB 脚本后,调试效率提升明显,因为你随时可以回退到“上一个能跑通的版本”,而不是在乱改中越陷越深。

代码方面,YALMIP 本身有部分示例,但移动储能预布局的完整开源代码确实少,建议以参考那篇文献的思路为骨架,按我上面说的模块自己搭建。Gurobi 学术授权可以免费申请,MATLAB 2023a 以上版本对 YALMIP 的兼容性都不错。求解过程中如果遇到数值警告,优先检查 M 值和二进制变量的上下界,这能解决 80% 的模型病态问题。

最后再补充一个我自己的小技巧:把每个实验的随机种子固定下来——MATLAB 里用rng(seed)控制场景生成。否则你每跑一次结果都在变,很难判断是自己算法改进带来的收益还是随机性带来的波动,这一点在写论文实验部分时尤其重要。固定种子后,任何对比实验的差异都能准确归因,实验结果才能服人。

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询