基于IEEE33节点的移动储能两阶段预布局与动态调度策略复现

这两年做配电网韧性方向的仿真一个结论越来越清晰移动储能的价值能否兑现很大程度上不取决于电池本身而取决于“什么时候知道该去哪”。同样是两辆2 MW/4 MWh的移动储能车预布局做得好的场景故障后关键负荷几乎不掉电布局没跟上故障发生6小时后储能车还堵在错误的节点上干等。这篇文章要拆解的就是这一类论文最核心的两阶段决策框架——面向配电网韧性提升的移动储能预布局与动态调度策略并基于IEEE33节点测试系统用Matlab完整复现。这个复现方案覆盖了从灾前选址、故障场景生成、两阶段优化建模到求解器调用的完整链路适合正在做电力系统韧性、主动配电网恢复、储能优化调度方向的研究生和工程师。你不需要自己从头搭模型只需要理解每一段代码背后的逻辑就能把它移植到自己的算例里。部分参数设置来自该方向文献中的常见配置复现时我按照最典型的方案补齐了细节后面会逐条说明哪里是通用的、哪里需要根据你自己的故障场景调整。1. 为什么移动储能能改变配电系统的韧性天花板1.1 从极端灾害后的负荷恢复场景说起配电系统韧性研究通常从一张灾后恢复曲线讲起极端事件导致多条馈线同时断线部分负荷脱离主网供电抢修队伍需要数小时甚至数十小时才能恢复送电这段时间内关键负荷只能靠本地电源支撑。实际的痛点在于传统配电自动化开关能切除故障区域但无法给孤岛内的负荷提供能量柴油发电机机动性强却存在噪声、燃料补给和碳排放问题。移动储能车的定位刚好卡在这个空隙里——它本质是一个可移动的电源缓冲器可以在故障前到达预测位置故障后根据实际断线情况再移动到最需要的地方临时支撑孤岛运行。我在复现许多论文时注意到一个共性事实这类模型很少追求复杂的灾后重构而是把重心放在“储能车的位置决策”上。因为配电网故障后的在线重构受到开关动作次数、保护配合等多重限制而移动储能的位置变化却能直接决定哪些负荷能被支撑。位置选对了整个恢复方案的可行域都会变大。1.2 移动储能与固定储能差别不在电池在“位置”固定储能和移动储能之间最本质的差异不是电池容量、不是功率等级而是位置是否随时间和故障场景变化。固定储能一旦安装其供电范围就受限于所在馈线的连通性如果故障点恰好把储能和关键负荷隔开储能再大也无能为力。对比维度固定储能移动储能位置可变性固定不可变灾前预部署、灾中动态转移对故障场景的适应性依赖安装位置可根据预测提前靠近风险区域利用率低平时可能闲置灾前可参与调峰灾后作为应急电源调度建模复杂度单层充放电决策位置路径充放电联合决策移动储能需要一体两面的“预布局”和“动态调度”来支撑预布局解决的是“灾前把车放在哪”的问题动态调度解决的是“灾后车往哪跑、在哪个节点出力”的问题。两件事分开看都不难难在它们要放在同一个优化框架里协同还要考虑路网通行时间、电池SOC变化、配电网潮流约束。这也正是本篇复现的核心价值所在。如果你之前只做过固定储能参与配电网调度的模型那么你已有的DistFlow约束、SOC递推约束都可以直接复用。真正需要从头写的是位置变量和移动约束这部分我会在数学模型章节重点拆解。1.3 这个复现方案解决什么问题具体到代码层面这个项目解决的是这样一类优化问题在极端事件预警发布后已知台风路径预测、线路故障概率和若干典型故障场景如何在IEEE33节点配电网中为若干台移动储能车确定灾前预布局阶段的最优初始停靠节点灾后动态调度阶段每个时段的移动路径和停靠节点每个时段移动储能在接入节点处的充放电功率系统整体的负荷削减量和韧性提升效果。这套代码跑完之后你能得到的东西很直观移动储能车的时空轨迹、节点电压分布、关键负荷恢复曲线、系统韧性指标改善程度等。它面向的读者既有刚入门的硕士生也有需要在算例中快速验证自己想法的科研人员。代码本身不需要太大改动就可以换成其他节点系统只要你能提供对应的导纳矩阵和拓扑信息。2. 两阶段决策框架拆解预布局和动态调度如何分工2.1 第一阶段预布局在不确定性中选一个“不后悔”的初始位置预布局发生在灾害尚未发生的预警阶段。这个阶段最大的麻烦是你根本不知道哪些线路会真实断线只能根据气象预报和历史统计知道大致概率。如果只按单一预测场景去布局一旦实际故障点偏移储能车可能离真正的失电区域很远。所以第一阶段的本质是随机规划里的“here and now”决策在故障场景不确定的情况下决定储能车的初始位置使得所有可能发生的故障场景下后续恢复效果的期望值最优。代码里通常用一组离散场景代表可能的故障情况每个场景包括断线支路集合和故障持续时间。预布局模型的目标函数通常写成$$\min_{x^{pre}} \sum_{s \in S} p_s \cdot Q(x^{pre}, \xi_s)$$其中 (x^{pre}) 是预布局决策(Q) 是第二阶段在具体故障场景 (\xi_s) 下的最优恢复代价(p_s) 是场景概率。也就是说第一阶段的“好位置”必须在所有可能故障场景的期望意义下表现不差。这是一种保守但工程上常用的妥协不确定性真实存在时追求“全场景平均最优”比“单一场景最优”务实得多。在IEEE33节点系统上MES初始位置通常被限定在一些关键候选节点集合中比如负荷中心的邻近节点、联络开关附近的节点。全节点放开也不是不行但会显著增加二进制变量数量求解时间会呈指数增长。我复现时采用的做法是先把33个节点按负荷重要度和支路故障概率排序再选出8到10个候选节点只在这些候选节点上给预布局变量。2.2 第二阶段动态调度故障信息明确后的路径-功率联合优化故障发生后断线集合、故障时段都成了已知信息。第二阶段要做的是在一个确定性环境下把移动储能车的“路径”和“出力”一起优化。第二阶段模型里有两个关键要素拓扑连通性和移动代价。由于故障线路被切断部分支路不可用储能车从一个节点到达另一个节点只能通过还在运行的线路或道路网。为了简化IEEE33节点的代码里一般只考虑配电线路拓扑作为移动通道默认相邻节点之间移动需要一个时段。第二阶段的决策包括每个时段储能车停靠在哪一个节点何时从当前节点转移到目标节点接入节点后充放电功率如何分配各时段哪些负荷被削减。目标函数一般是最小化加权失负荷量负荷权重根据重要度设定比如医院、通信基站等关键负荷权重高普通居民负荷权重低。权重设计直接影响优化结果这一点在做敏感性分析时要格外注意。2.3 两阶段衔接逻辑与随机规划视角从数学上看这是一个典型的两阶段随机规划。预布局变量是第一阶段的原始决策动态调度变量和恢复代价是第二阶段的响应。代码实现时最常见的做法有两种。第一种是把所有场景的两阶段模型写成一个大规模MILP一次性求解。这种方式理论最优性有保障代价是二进制变量非常多33节点系统的人规模灾害场景可能都跑不动。第二种是先独立求解第一阶段得到预布局位置后再逐场景求解第二阶段动态调度。这种方式计算压力小很多代价是一阶段没有显式考虑二阶段的响应容易得到保守方案。我在代码里采用的是折中思路第一阶段用场景期望成本作为目标但实际求解时只对每个场景做二阶段的预评估再用迭代方式逼近两阶段联合最优。如果你只是复现出基本曲线用一次性MILP就够了如果你的场景数超过10个建议用我这种迭代方法否则求解时间会很难看。3. IEEE33节点系统构建和故障场景设置3.1 为什么用IEEE33而不是更大或更小的系统IEEE33节点系统几乎是配电网优化研究的事实标准。它由33个节点、32条支路组成正常运行呈辐射状结构基准电压12.66 kV基准功率10 MVA总负荷约3715 kW加2300 kvar。规模不大但足够承载故障重构、储能调度、分布式电源接入等多种典型场景。选择它有三个现实原因。其一是参数公开几乎所有论文附录都能找到完整支路阻抗和节点负荷数据复现时不用花时间在数据清洗上。其二是节点规模适中二阶锥松弛后的MILP模型能在几秒到几分钟内求解调试效率远高于动辄上百节点的实际馈线模型。其三是IEEE33本身带有联络开关便于在后续扩展中对比“仅储能调度”和“储能网络重构”两种策略的效果差异。如果你的研究方向是韧性提升策略的机理分析而非工程落地IEEE33完全够用如果偏实际工程验证后续可以换成IEEE123节点代码的主体结构只需要改拓扑数据文件。3.2 从原始参数表到Matlab数据结构在Matlab中搭建IEEE33系统核心是把网络拓扑和电气参数组织成可供YALMIP建模使用的矩阵形式。我复现时常用的数据结构是这样的% 支路数据: [起始节点, 终止节点, 电阻(pu), 电抗(pu), 容量(pu)] branch [ 1 2 0.0922 0.0470 0.6; 2 3 0.4930 0.2511 0.6; 3 4 0.3660 0.1864 0.5; % ... 其余支路省略 ]; % 节点数据: [节点编号, 有功负荷(kW), 无功负荷(kvar), 重要度权重] bus [ 1 100 60 0.4; 2 90 40 0.5; 3 120 80 1.0; % ... 其余节点省略 ];节点编号从1开始节点1是松弛节点代表上级电网。支路电阻和电抗要换算成标幺值后再写入约束否则数值会出现巨大差异。换算公式是$$Z_{pu} \frac{Z_{ohm}}{Z_{base}}, \quad Z_{base} \frac{V_{base}^2}{S_{base}}$$12.66 kV、10 MVA时(Z_{base})约为16.03欧姆。之前有人直接用欧姆值代入DistFlow方程得到的结果电压全部偏离实际范围就是没做标幺化导致的。系统参数我会在代码注释里完整给出核心表如下参数数值说明基准电压12.66 kV线电压有效值基准功率10 MVA三相功率基准节点数33含根节点支路数32不含联络开关总负荷3715 kW 2300 kvar峰值工况3.3 极端事件故障场景的典型建模方式极端事件故障场景的建模在代码里落地的形式就是一组“断线支路集合持续时间”。以台风灾害为例通常的做法是从风险分析中得到每条支路的故障概率再通过蒙特卡洛采样生成典型场景或者直接人工挑选高风险支路构造代表场景。我在复现中采用三组典型场景覆盖不同严重程度场景断线支路故障时段场景概率S17-8第2至第7时段0.5S214-15, 24-25第3至第10时段0.3S37-8, 14-15, 24-25第2至第11时段0.2选择7-8、14-15、24-25是因为它们分处不同的馈线分支故障后形成的孤岛在空间上错开能充分测试储能车在不同方向之间转移的能力。这些支路断线后节点8至18、节点24至25一带会脱离主网形成两个互不相连的孤岛移动储能车只能通过还在连通的弧段间接靠近。故障场景的构造直接决定了预布局的成果。如果你想让储能车全部集中在故障概率最高的支路附近可以缩小场景集合如果你想测试脚本的鲁棒性就增加更多随机故障组合。推荐先跑通三场景再逐步扩展到更多场景观察计算时间随场景数的增长趋势。4. 数学模型怎么搭目标、约束和线性化处理4.1 韧性指标与目标函数韧性提升效果的量化常见指标是失负荷量或可供电量。考虑到不同节点负荷重要度不同我会在目标函数里给各节点设置权重 (\omega_i)目标是整个调度周期内加权失负荷总量最小$$\min \sum_{t \in \mathcal{T}} \sum_{i \in \mathcal{B}} \omega_i \cdot P^{shed}_{i,t} \cdot \Delta t$$其中 (\Delta t) 是调度时段长度通常取1小时。权重高的节点会被优先保障这符合实际应急响应的优先级逻辑。另外如果灾害期间上级电网还能对部分区域供电那么根节点注入功率的爬坡限制也要加入约束避免模型通过无限增大根节点出力来“作弊”似的消除全部失负荷。这样设置的意义是当配电网内部断开时上级电网能供给的能量有限移动储能的补充才有意义。4.2 移动储能的时空-电量约束移动储能建模需要同时处理位置、路径、充放电功率和SOC四个维度。最容易漏掉的是路径与位置之间的逻辑关系我复现时用下面的变量表示pos(i,t,k) % 第k辆储能在t时段是否停靠节点i二进制变量 P_dis(i,t,k) % 第k辆储能在t时段节点i处的放电功率 P_ch(i,t,k) % 第k辆储能在t时段节点i处的充电功率 soc(k,t) % 第k辆储能在t时段的荷电状态约束分四组位置唯一约束同一时刻每辆车只能停靠一个节点移动距离约束相邻时段之间停靠节点必须在网络连通弧段内防止“瞬移”充放电功率上下限和SOC递推约束接入节点才允许放电的逻辑约束。用数学形式写前三条如下$$\sum_i pos(i,t,k) 1, \quad \forall t,k$$$$pos(i,t,k) pos(j,t1,k) \le 1, \quad \forall (i,j) \notin \mathcal{E}$$$$SOC(k,t1) SOC(k,t) \left(\eta_{ch} P^{ch}{k,t} - \frac{P^{dis}{k,t}}{\eta_{dis}}\right)\Delta t / E_k$$第三条里的 (E_k) 是储能容量标幺值(\eta_{ch}) 和 (\eta_{dis}) 分别是充放电效率典型值取0.95和0.95。这些参数直接影响储能车能支撑多长时间一定不能让SOC超过[0.1, 0.9]的合理运行区间。4.3 DistFlow潮流方程与二阶锥松弛配电网的辐射状结构使DistFlow方程成为业界标准选择。对任意支路 (i \rightarrow j)完整方程包含有功、无功和电压三个递推关系在标幺值下可以写成$$P_{ij,t} P_{j,t}^{load} - P_{j,t}^{dis} P_{j,t}^{ch} \sum_{m \in child(j)} P_{jm,t}$$$$Q_{ij,t} Q_{j,t}^{load} \sum_{m \in child(j)} Q_{jm,t}$$$$v_{j,t} v_{i,t} - 2(r_{ij}P_{ij,t} x_{ij}Q_{ij,t})$$严格来说DistFlow还包含一个非线性项 ((r_{ij}^2x_{ij}^2)(P_{ij,t}^2Q_{ij,t}^2)/v_{i,t})忽略它需要足够小的支路阻抗。在IEEE33节点系统上忽略后的误差通常小于1%因此绝大多数论文都采用上述线性化形式。需要更精确的表达时可以用二阶锥松弛$$\left| \begin{bmatrix} 2P_{ij,t} \ 2Q_{ij,t} \ v_{i,t} - l_{ij,t} \end{bmatrix} \right|2 \le v{i,t} l_{ij,t}$$加上这个约束后模型变成二阶锥规划SOCP求解器如Gurobi可以高效处理。如果你的模型里出现了潮流约束的乘积项导致无法收敛优先考虑是不是漏了锥约束。4.4 配电网辐射状拓扑约束故障发生后网络可能因为部分支路被切开而分裂成若干个孤岛。如果允许网络重构就必须显式约束每个孤岛内节点之间的连通性比较简单的做法是单商品流法。在代码里我引入一个虚拟的“流”变量 (flow_{ij,t})在每个时段确保每个非根节点都能从某个父节点获得单位虚拟流且虚拟流只能从根节点沿线注入。这样能防止优化结果中出现环网或脱离电网的孤立带电孤岛。不过如果你一开始就不打算开放联络开关只需要保证储能接入节点和需要恢复的负荷节点处于同一个连通区域拓扑约束可以大幅简化。我复现时优先保证了预布局位置和动态接入位置均在正常连通区域内所以这部分约束写得比较轻。一旦你后续尝试“孤岛划分储能支撑”的联合优化这里就需要认真补齐。5. Matlab代码实现从建模到求解的关键链路5.1 开发环境YALMIP Gurobi/CplexMatlab平台上做优化建模我最推荐的环境组合是YALMIP加Gurobi。YALMIP是建模语言层负责把约束和目标函数转成标准优化模型Gurobi是底层求解器负责实际求解MILP或MISOCP。环境配置本身不复杂但有几个注意点YALMIP需要放到Matlab路径下并确保当前版本与Matlab版本匹配Gurobi需要安装并配置许可证激活之后做一次gurobi_setup在YALMIP中通过solvesdp或optimize时指定solver选项为gurobi不要让它自动选求解器。首次运行可以用check命令逐条查看约束是否满足这个习惯能帮你在早期发现建模错误比等模型解完再处理高效得多。5.2 变量定义与约束拼装的核心代码核心变量定义如下nbus 33; % 节点数 NT 24; % 调度时段数 nes 2; % 移动储能车数量 pos binvar(nbus, NT, nes, full); % 位置指示 P_ch sdpvar(nbus, NT, nes, full); % 充电功率 P_dis sdpvar(nbus, NT, nes, full); % 放电功率 soc sdpvar(nes, NT, full); % 荷电状态 P_shed sdpvar(nbus, NT, full); % 失负荷功率约束拼装时建议把同类约束单独放到一个函数里不要在总脚本里堆千行代码。我通常会把约束分成constraints_mes.m、constraints_distflow.m、constraints_network.m三个模块出问题时可以定位到具体模块。移动储能相关的约束在代码中长这样Constraints []; % 每个时段每辆车只能停在唯一节点 for k 1:nes for t 1:NT Constraints [Constraints, sum(pos(:,t,k)) 1]; end end % SOC递推与充放电逻辑 for k 1:nes for t 1:NT-1 Constraints [Constraints, ... soc(k,t1) soc(k,t) ... (eta_ch * sum(P_ch(:,t,k)) - sum(P_dis(:,t,k))/eta_dis) * dt / E_cap]; end end接入节点约束也很重要如果储能车没有停在某个节点那么它在这个节点的充放电功率必须为零。用一个大M约束实现M 5; % 大于储能最大充放电功率的上限 for k 1:nes for t 1:NT Constraints [Constraints, ... P_dis(:,t,k) M * pos(:,t,k)]; Constraints [Constraints, ... P_ch(:,t,k) M * pos(:,t,k)]; end end如果不加这两个约束模型会“找出”一个漏洞储能车通过不在某个节点却在该节点放电来满足负荷需求这在物理上不可能。这个坑我见过不止一次踩进去后排查了很长时间。5.3 参数调试与求解配置求解器的参数设置直接影响收敛速度和结果质量。我在复现中用的是options sdpsettings(solver,gurobi, ... verbose, 2, ... gurobi.MIPGap, 0.01, ... gurobi.TimeLimit, 600);MIPGap设为0.01表示允许1%的最优性间隙对于这种规模的两阶段模型来说是合理折中既能保证结果不差又能避免求解器在最优解附近反复验证花费过长时间。如果你需要更快的调试体验可以先把调度周期缩短到6到8个时段验证模型逻辑无误后再恢复到24时段。这个方法看起来笨实际却非常管用尤其是处理二进制变量相关约束时一次全时段求解失败往往很难判断是哪个约束写错了。5.4 结果后处理与图形输出调试正确后输出结果我一般用三类图形第一类是电源出力与失负荷曲线。横轴是时段纵轴是功率叠加显示储能放电功率和系统总失负荷功率能够直观看到移动储能在故障时段顶住了多少负荷。第二类是储能车位置热力图。按节点号按时段展示储能车停靠位置检查是否符合直觉。如果某辆车在故障后第3个时段出现在距离故障点很远的节点上且中间没有经过任何路径大概率是移动距离约束漏写了。第三类是节点电压分布曲线用来检查潮流结果是否越限。故障期间部分节点电压可能偏低但不应低于0.90 pu如果图中出现低于0.85的节点要先排查DistFlow方程的方向和符号。6. 复现过程中避不开的坑和排查思路6.1 求解器直接报Infeasible先查三件事第一次跑通全模型时遇到Infeasible别慌按下面的顺序排查第一检查量纲。支路阻抗、负荷功率、储能容量是否全部统一到了标幺值或同一单位体系。混用kW和MW、欧姆和标幺值是初学者最容易犯的错。第二用check(Constraints)逐条查看约束是否满足。YALMIP会输出每条约束的最大残差残差很大的那条通常就是问题所在。第三检查故障场景下的网络连通性。如果某个场景中断线支路把一部分节点完全孤立而储能车又无法进入该区域时模型天然不可行。这时需要给失负荷变量设置一个足够大的上限或允许该区域完全断电。症状可能原因排查方法全模型Infeasible量纲不一致统一标幺值特定场景Infeasible网络断线导致孤岛不可达检查断线组合结果退化大M值过小调大M并验证极限情况6.2 储能车在结果里“瞬移”问题出在哪结果图上储能车位置在相邻时段从节点5跳到了节点25中间跨越了故障断开的线路这种情况就是“瞬移”。根因通常是漏了移动距离约束或者约束写得太宽松只限制了两个时段间不能完全任意但没有限制在拓扑连通距离内。IEEE33规模的网络比较稳妥的写法是维护一个节点邻接矩阵adj若两个节点之间有运行支路则adj(i,j)1否则为无穷大。然后约束% 禁止储能车在不可达节点间瞬移 for k 1:nes for t 1:NT-1 for i 1:nbus for j 1:nbus if adj_full(i,j) 1 Constraints [Constraints, ... pos(i,t,k) pos(j,t1,k) 1]; end end end end end这种方式虽然循环层数多但约束很直观。想要更高效可以只对adj值大于1的节点对生成约束避免空转到全节点组合。6.3 二阶锥松弛不收敛的调试顺序如果你在模型中保留了完整的DistFlow非线性项并使用二阶锥松弛偶尔会遇到求解不收敛或对偶间隙无法下降的情况。我建议按以下顺序排查先检查锥约束是否被YALMIP正确识别。YALMIP中用cone()函数显式定义锥约束比写成不等式更稳定。其次检查松驰的紧致性。求解后计算每个时段每条支路的锥松弛间隙如果间隙过大说明松弛不够紧模型结果离物理可行解有距离。可以尝试在目标函数中加入一个很小的锥松弛惩罚项促使求解器收敛到紧致解。最后检查求解器选项。BarQCPConvTol这个参数对MISOCP求解影响很大默认值在部分模型上表现不好可以试着从1e-6调整到1e-8代价是求解时间上升。一般情况下IEEE33节点系统不需要调整到这个精度就能得到满意的结果。6.4 求解时间爆炸用什么思路压缩当场景数量超过5个二进制变量会迅速膨胀Gurobi有时要跑十几分钟才能找到可行解。常见的压缩思路有三种。第一种是合并时段。把1小时粒度改为2小时粒度时段数从24降到12二进制变量直接减半精度损失在可接受范围内。第二种是候选节点筛选。预布局阶段只开放一小部分候选节点不要全节点招投标。第三种是场景削减。对原始气象场景聚类用少数代表性场景替代完整场景集这是随机规划的常用技巧。我实际测试下来在IEEE33节点上采用8个候选节点、24时段、3个典型场景、2辆储能车Gurobi约在30到60秒内能稳定求出MIPGap在1%以内的解。如果你的求解器跑到几分钟还没动静检查是不是预布局候选节点没有限制住。7. 复盘这套代码跑通之后还可以往哪个方向延伸整个复现做完我比较大的感受是两阶段模型的价值不在代码本身而在你能拿它快速验证什么。所有做韧性研究的人最终都会面临同一类问题——如果极端事件持续时间更长如果储能车数量不同如果负荷权重改变系统韧性改善的边际效应是多大用这套代码你可以在几分钟内通过改参数回答上述问题。比如调整储能车数量从2辆到4辆观察总失负荷量下降的幅度就能粗略估算当前系统下移动储能的最优配置数量调整故障持续时间从6小时到10小时就能看到预布局位置权重是否需要向更关键的节点偏移。这些敏感性分析在实际项目中往往是拼论文或者写报告最有价值的部分。我个人踩过几次坑后也比较确定地认为现阶段做移动储能调度还值得往两个方向扩展一是时变路网和交通约束的融合让储能车移动时间不再是固定值而是与道路拥堵情况相关二是与抢修队联合调度让储能车和维修人员共享决策。如果你已经有这套两阶段预布局动态调度代码往这两个方向加约束都不会太困难核心变量和框架完全兼容。