简介面向配电网韧性提升的移动储能预布局与动态调度策略MATLAB程序适用于电力系统韧性评估、储能规划与动态调度研究目标用户为高校师生、科研院所工程师及配电网规划人员。程序包含灾前预防与灾后恢复两大模块灾前部分采用两阶段鲁棒优化框架将移动储能配置成本与负荷削减成本作为目标函数并考虑光伏出力不确定性通过大M法和列与约束生成算法CCG迭代求解灾后部分采用混合二阶锥规划直接优化得到最优负荷削减量与各设备出力方案。压缩包约3.56MB以MATLAB源码文件为主需配合Gurobi求解环境运行。资源与《面向配电网韧性提升的移动储能预布局与动态调度策略》论文配套已有462人学习适合需要复现文献算法、开展韧性提升仿真或进行算法对比研究的读者参考。1. 这个标题到底在解决什么问题台风还没来储能车该停在哪台风预警发布到预计登陆往往只有十几个小时配电网里那几台移动储能车MES提前停在哪个节点直接决定了灾后孤岛能恢复多少关键负荷等故障点陆续明确、路况动态变化车又该按什么顺序往哪开、到了以后以多大功率放电又是一层实时决策。标题里“预布局与动态调度”就是把这两层决策统一装进一个优化模型里用 MATLAB 求解灾前决定车的停放位置灾中决定车的移动路径和充放电计划。它解决的不是“多买几台车”而是“把已有的车停在正确的位置、在正确的时间送到正确的节点”。这套方法适合配电网规划、应急管理的从业者也适合正在用 MATLAB 做混合整数规划课题的研究生——有可复现的建模框架也有能直接抄走的主程序骨架。2. 先把问题写明白移动储能预布局与动态调度的两阶段数学建模2.1 为什么是“预布局 动态调度”而不是单一优化固定储能一旦安装位置就锁死了移动储能车的价值恰恰在于“可移动”。灾害发生之前精确故障点未知只有预测的故障场景和节点重要度所以要先做预布局——把车放在收益期望最高的备用节点。灾后故障信息逐步清晰交通条件也在变再依靠动态调度来决定每辆车往哪跑、在哪放电、是否回撤充电。这两件事如果分开优化最常见的翻车方式是预布局阶段把所有车都堆在重负荷节点旁边结果实际故障出在另一端车赶过去时黄花菜都凉了反过来动态调度阶段无视车的初始位置给出一个物理上根本来不及走的路线。把两个阶段放进同一个优化问题里用 t0 时刻的位置约束把“预布局结果”和“动态调度起点”咬合在一起才是工程上可靠的建模方式。2.2 目标函数韧性不只是“少停电”通用的目标函数写成min Σ_t Σ_k c_move · u_{k,t} Σ_t Σ_n λ_n · shed_{n,t}第一项是所有移动储能车在调度周期内的移动成本u_{k,t}是车的移动标志位第二项是失负荷惩罚shed_{n,t}是节点 n 在时刻 t 的削减负荷量λ_n是节点重要度乘数。医院、供水、通信基站的 λ 取值要明显高于普通居民负荷。λ 的取值直接决定调度行为取太小车会为了省移动成本而不去关键节点取太大车会不计代价地绕远路。我一般会让最大 λ 比平均移动成本高一个数量级保证“恢复关键负荷优先于节约里程”。如果要算韧性指标一般用韧性提升率 (无MES失负荷 - 有MES失负荷) / 无MES失负荷或者统计“关键负荷恢复时长”这类指标。第 4 章给结果时会具体演示怎么算。2.3 核心约束把“时间”写进模型这个问题的约束可以分成四类每一步都是工程上踩过坑的地方。第一类是位置一致性约束。每一辆车在任意时刻只能处于一个物理状态停靠在某个节点或者正在移动途中。写成Σ_n y_{k,n,t} move_{k,t} 1y_{k,n,t}是二元变量为 1 表示车 k 在 t 时刻停靠节点 nmove_{k,t}表示车在 t 时刻处于移动状态。第二类是移动时序约束。车从节点 n 到节点 m 需要τ_{n,m}个时段用出发事件变量z_{k,n,m,t}表示“车 k 在 t 时刻从 n 出发去往 m”。出发之后到到达之前车不能停靠任何节点。这条约束是模型里二元变量的大头也是最容易把模型写到无解的地方。第三类是 SOC 递推约束对应标题里“动态调度”的能量维度soc(k,t1) soc(k,t) - p_dis(k,t)·dt/E p_ch(k,t)·dt·η/E同时要有 SOC 上下限以及充放电互斥约束。注意充放电效率是不对称的充电 η 和放电 η 分开写。第四类是配电网功率平衡约束。对每个孤岛节点注入功率加上移动储能放电功率要满足负荷需求减去削减量。如果只做移动储能调度、不做配电网重构用简化的节点功率平衡已经够用误差主要体现在网损这个在避坑章再展开。2.4 离散时间系统视角步长怎么选在 MATLAB 里做配电网移动储能调度本质上是把连续时间系统离散化状态量是每辆车的位置、SOC、充放电功率控制输入是“出发 / 停靠 / 充电 / 放电”这些离散事件。这正是“离散时间系统”的典型建模场景。调度步长 dt 直接决定模型大小dt1h 时一个 24 时段调度问题时段数是 24dt 改成 15min时段数变成 96二元变量规模直接翻四倍求解时间可能涨一个数量级以上。我的习惯是做两套步长嵌套预布局阶段用 1h 步长筛方案动态调度阶段对故障发生后的关键 6 小时加密到 15min。不要一上来就全时段 15min 建模否则模型还没跑起来先把求解器难死了。下面用一张表把这四类约束汇总后面写 YALMIP 代码时可以直接对照。约束类别数学表达简写作用位置唯一Σ_n y(k,n,t) move(k,t) 1防止同一辆车分身移动时序z(n,m,t) → 到达时刻恢复停靠保证旅行时间SOC 递推soc(t1) soc(t) p_ch·η/E - p_dis/(η·E)能量守恒功率平衡Σ p_dis(k,n,t) P_grid(n,t) ≥ P_load - shed供需平衡3. 用 MATLAB 把模型跑起来YALMIP 建模与主程序骨架3.1 为什么用 YALMIP 而不是手写 intlinprog通用做法是用 YALMIP 作为建模层把数学约束直接“翻译”成代码再用 MATLAB 优化工具箱自带的intlinprog或者商业求解器求解。YALMIP 的核心价值是省去手工拼稀疏矩阵的体力活尤其当变量带有(n,m,k,t)四维下标时手写矩阵拼接极易错位。如果环境里没有 Gurobi 或 CPLEXsdpsettings(solver,intlinprog)也能跑通小规模算例但模型一旦超过几千个二元变量还是建议装商业求解器。另外要确认 MATLAB 里装了优化工具箱这直接关系到intlinprog是否存在。YALMIP 在 R2023b 之后的主版本兼容性普遍正常安装后运行yalmiptest能看到当前可用的求解器列表。3.2 主程序骨架数据加载、变量定义、约束与求解下面是完整可跑通的主程序结构。代码基于 IEEE 33 节点系统做了简化但包含了预布局和动态调度两层核心。%% 主程序配电网移动储能预布局与动态调度 clear; clc; close all; yalmip(clear); %% 基础数据 load(data33.mat); % 包含 N,T,dt,负荷曲线,线路参数,travel_time K 3; % 移动储能车数量 E 1000; % 每台车容量[kWh] P_max 500; % 最大充放电功率[kW] soc0 0.5; % 初始SOC lambda 1e3; % 失负荷惩罚系数 tau travel_time; % 节点间移动时间矩阵 (N x N) %% 决策变量 y binvar(N,K,T,full); % y(n,k,t)1 车k在t时刻停靠节点n move binvar(K,T,full); % move(k,t)1 车k在t时刻处于移动 u_dep binvar(N,N,K,T,full); % u_dep(n,m,k,t)1 车k在t时刻从n出发前往m p_dis sdpvar(N,K,T,full); % 放电功率 p_ch sdpvar(N,K,T,full); % 充电功率 soc sdpvar(K,T,full); % 荷电状态 0~1 shed sdpvar(N,T,full); % 失负荷量 Constraints []; %% 预布局t1时刻的位置即灾前停放位置 Constraints [Constraints, sum(y(:,:,1),1) 1]; % 每辆车只能停一处 %% 动态调度位置唯一性 for k 1:K for t 1:T Constraints [Constraints, sum(y(:,k,t)) move(k,t) 1]; end end %% SOC递推与上下限 for k 1:K for t 1:T-1 Constraints [Constraints, soc(k,t1) soc(k,t) ... sum(p_ch(:,k,t))*0.95/E - sum(p_dis(:,k,t))/(0.95*E)]; end Constraints [Constraints, 0.2 soc(k,:) 0.9]; Constraints [Constraints, soc(k,1) soc0]; end %% 放电功率与停靠状态关联关键约束 for k 1:K for t 1:T Constraints [Constraints, p_dis(:,k,t) P_max .* y(:,k,t)]; Constraints [Constraints, p_ch(:,k,t) P_max .* y(:,k,t)]; end end %% 目标函数 Objective sum(sum(lambda * shed)) sum(sum(1e3 * sum(sum(u_dep,2),3))); %% 求解 ops sdpsettings(solver,intlinprog,verbose,1,mipgap,0.01); sol optimize(Constraints, Objective, ops);这段代码的逻辑顺序很重要先定义变量再逐类加约束最后求解。y(:,:,1)对应 t1 时刻车的位置这就是预布局结果p_dis与y的关联约束保证车只有停靠在某个节点时才能放电这一步漏掉的话模型会出现“车在半空中放电”的荒谬结果。参数E、P_max、soc0要按实际储能车参数改lambda的取值会影响车往哪跑建议先跑一版看轨迹再调。3.3 预布局层单独输出有什么用虽然两层在一个模型里工程上仍然需要单独看预布局结果。做法很简单求解后打印y(:,:,1)中为 1 的位置。for k 1:K n0 find(y(:,k,1) 0.5); fprintf(MES %d 灾前停放节点: %d\n, k, n0); end这段代码的价值是让预布局方案落地成一纸调度单哪台车停哪个节点、什么时候出发。如果 MES 是租赁的还可以把 K 从常数改成整数变量在目标函数里加一项租赁成本这样模型就能同时回答“租几台、停在哪”。3.4 移动时间大于一个时段时怎么改前面 3.2 的代码假设任意两次停靠之间只需要 1 个时段这在真实配电网里太理想了。实际节点间移动时间可能是 2 到 5 个时段需要把移动时序约束改成“从 n 到 m 需要 τ(n,m) 个时段内不能停靠”。YALMIP 里最简洁的写法是用impliesfor k 1:K for t 1:T for n 1:N for m 1:N if n ~ m t tau(n,m) T Constraints [Constraints, ... implies(u_dep(n,m,k,t), y(m,k,ttau(n,m)) 1)]; Constraints [Constraints, ... implies(u_dep(n,m,k,t), ... sum(y(:,k,t1:min(ttau(n,m)-1,T))) 0)]; end end end end end逻辑说明implies(A,B)表示当 A 为真时 B 必须为真。这里 A 是“车在 t 时刻从 n 出发到 m”B 包含两个条件——车在tτ(n,m)时刻到达 m并且在中间那些时段不占用任何停靠位置。sum(y(:,k,t1:...)) 0中的维度需要和 y 的结构严格对应写错会直接报维度不匹配或产生不可行约束。implies的求解性能略差于显式线性化但换来的是代码可读性和正确性对中小规模算例足够。3.5 求解器选择的黑匣子YALMIP 只是建模层真正算模型的是求解器。intlinprog是 MATLAB 自带的混合整数线性规划求解器优势是不需要额外安装缺点是二元变量上万以后求解时间明显变长。商业求解器 Gurobi 和 CPLEX 在 YALMIP 里通过sdpsettings(solver,gurobi)调用求解速度快一到两个数量级适合做多场景蒙特卡洛。关于“黑匣子”这个词求解器内部的分支定界过程对大多数用户确实是个黑匣子。你只能通过输出日志观察 gap 下降曲线。我常设mipgap在 0.01 到 0.05 之间设太小求解时间爆炸设太大结果可能偏离最优解很多。预布局方案允许一点次优动态调度则建议把 gap 收得更紧一些。4. 从三节点到 IEEE 33 节点参数怎么设、结果怎么读4.1 数据准备先跑小系统再上大系统我建议新手从三节点系统起步先把模型逻辑理顺再迁移到 IEEE 33 节点。三节点系统的 N3、T6、K1二元变量规模极小intlinprog 秒解。跑通之后再切换数据文件把 N 改成 33、T 改成 24。IEEE 33 节点的标准数据可以从 Matpower 的 case33 里导出包括节点负荷、线路阻抗、拓扑关系。移动储能车参数按主流 500kW / 1000kWh 车型设定SOC 范围 0.2 到 0.9初始 SOC 0.5。故障场景设定为某几条馈线同时断开把配电网切成若干孤岛。节点间移动时间矩阵tau由路网距离除以平均车速得到单位要换算成调度步长的整数倍。这一步很容易把单位搞混比如平均车速 30km/h、步长 1h、节点距离 12km移动时间就是 0.4 小时四舍五入到 1 个时段。4.2 关键参数表抄作业直接用下面这组参数是我在类似算例上反复调过的一组起始值直接套用基本不会出大问题。参数取值说明节点数 N33IEEE 33 节点标准系统调度时段 T2424 小时滚动调度步长 dt1h预布局阶段使用储能车数量 K3可改为待优化整数变量储能容量 E1000 kWh单台车电池容量最大功率 P_max500 kW充放电峰值功率SOC 下限 / 上限0.2 / 0.9保护电池寿命失负荷惩罚 λ1e3普通负荷默认权重移动惩罚 c_move1e3防止车辆频繁移动MIP gap0.01商用求解器一般 10 秒内收敛这套参数的经验是λ 和 c_move 必须保持同量级λ 略大于 c_move。如果 λ 远超 c_move模型会为了救一个节点让三台车全开过去实际交通条件根本不允许反过来如果 c_move 太大车就停在原地不动动态调度退化成固定储能。4.3 结果怎么读韧性提升率与轨迹表求解完成后第一件事不是看目标函数值而是算三个数无 MES 时的系统总失负荷、有 MES 时的系统总失负荷、关键负荷恢复时长。以前面那组参数在一组典型故障场景下的结果为例无 MES 时总失负荷约 40MWh有 MES 且做动态调度时总失负荷约 12MWh韧性提升率就是 (40-12)/40 70%。看结果时不能只盯着提升率还要看负荷削减曲线在时间轴上的分布——很多模型把负荷削减集中在前几个时段后面完全恢复这种情况对用户感知并不好因为停电初期恰恰是医疗和通信最需要保障的时候。惩罚系数λ按节点差异化设置后调度结果会自动把削减高峰压后这是目标函数设计时容易忽略的点。每台车的运行轨迹可以用一张表输出直观判断调度方案是否合理。时刻MES1 位置MES1 动作MES1 SOCMES2 位置MES2 动作MES2 SOC1节点12停靠0.50节点25停靠0.502节点12放电 500kW0.45节点25移动中0.503节点18移动中0.45节点24放电 300kW0.45读取轨迹表时重点检查两条是否存在“移动中放电”的时段是否存在“同节点双车同时放电但功率超过该节点负荷”的情况。这两类问题是模型约束漏写或参数错误的最常见信号。4.4 对比实验固定储能 vs 移动储能动态调度为了让结果更有说服力一般再做一个固定储能对照实验。做法是把 t1 时刻的位置固定死且不允许移动即去掉move变量和u_dep变量让每台车只能在一个节点放电。对比结论通常很明显固定储能在故障点与预置点重合时效果好故障点偏移时效果锐减移动储能通过动态调整位置可以把恢复范围扩大两到三个节点。这个对比也是判断“预布局 动态调度”到底值不值得做的关键论据。5. 避坑清单求解器找不到、模型无解、结果反复横跳5.1 现象YALMIP 报“No suitable solver for this problem”原因是 YALMIP 没有识别到任何可用的混合整数求解器。常见触发点有两个一是 MATLAB 安装时没勾选优化工具箱intlinprog根本不存在二是 YALMIP 路径没有被正确添加到 MATLAB 搜索路径或者 YALMIP 版本太老与 R2024a 之后的 MATLAB 不兼容。解决的固定动作是先运行which intlinprog如果返回intlinprog not found说明优化工具箱缺失需要重装或修复 MATLAB然后在 MATLAB 里执行addpath(genpath(yalmip文件夹路径))并savepath。最后运行yalmiptest确认intlinprog在求解器列表里。这一步是 YALMIP 建模的“开机自检”别跳过。5.2 现象模型求解极其缓慢两小时不出结果二元变量爆炸是头号原因。以 33 节点、24 时段、3 台车、全连接移动变量为例u_dep的完整规模是 N×N×K×T 33×33×3×24 78408 个二元变量其中大量 nm 的变量根本没有意义。解决思路是“砍变量 放宽 gap”。把u_dep定义成稀疏变量只保留n ~ m的元素把移动范围限制在故障可能涉及的节点集合内不要做成全连接再把mipgap从 0.01 放宽到 0.05。还有一种更实用的做法是分层求解预布局单独跑一次 MILP把解出的y(:,:,1)固定再作为初始条件跑动态调度这样两层分别求解规模各小一半。5.3 现象求解成功但结果出现“车在半空中放电”这是典型的约束漏写。如果只定义了p_dis和p_ch的上下限没有把放电功率和停靠位置关联求解器会利用这个漏洞让车一边移动一边放电目标函数还很好看实际完全不可行。解决办法就是 3.2 代码里那句关联约束p_dis(:,k,t) P_max * y(:,k,t)同时给p_ch也加上同样的约束。为什么要给充电也加因为移动过程中充电同样物理不可能。检查结果时在轨迹表里搜一遍“状态为移动但功率非零”的时段能快速定位这类问题。5.4 现象调度结果反复横跳车在相邻节点之间来回跑目标函数里移动惩罚太小或根本没加时模型会让车无意义地反复移动因为每多跑一个节点就能重新优化一次放电位置目标函数反而更优。从结果看就是轨迹表里 MES 每小时都在换节点。解决方法是把移动成本调高再加“全调度周期内每辆车最多移动 M 次”的约束。M 一般取 3 到 5既能保证灵活性又避免振荡。调参经验是先看无移动惩罚时的最优解如果移动次数超过 T/3就把移动惩罚系数翻倍再跑直到轨迹稳定。5.5 现象预布局结果和动态调度结果对不上预布局说车停在节点 A但动态调度一开始就让车从节点 B 出发或者动态调度方案在 t1 时刻要求车出现在它根本来不及到达的节点。根源在于两层模型没有共享同一个 t1 位置约束。如果两层是分开建模求解的必须把预布局输出的y(:,:,1)显式写为动态调度的初始条件约束y_init(:,k) y_pre(:,k)。如果合在一个模型里要检查sum(y(:,:,1),1) 1这个约束是不是写成了对节点求和而非对车辆求和维度搞反就会出现一辆车占多个节点的矛盾。这个错误很隐蔽报错却常常是“infeasible problem”让人摸不着头脑。6. 进阶打法滚动时域重调度、蒙特卡洛验证和轨迹可视化读完前五章你已经能跑通一个完整的“预布局 动态调度”模型。接下来的进阶打法有三个方向滚动时域、不确定性验证、可视化输出。滚动时域是把 24 小时的静态调度改成每小时重跑一次的做法。每执行 1 个时段把当前实际故障信息、交通状态和车辆位置更新进去重新求解剩余时段的调度问题。这样模型就有“后悔药”可吃灾前预测的故障点错了下一时段就调整路线。实现上只需要把循环包一层每次把当前t0时刻的状态作为新模型的初值。蒙特卡洛验证是检验方案韧性最常用的手段。按故障位置和持续时间抽样生成 100 到 1000 组场景对每组场景跑一遍调度模型统计韧性提升率的期望和 5% 分位数。这一步能回答“预布局在最坏场景下会不会失效”也是写报告和论文时最有说服力的一张图。代价是计算量大建议先用低精度 gap 快速跑确认无误后再加密。轨迹可视化做的事很简单用stairs或plot把每台车的节点位置随时间的变化画出来再把 SOC 曲线画在第二张子图里两张图上下对齐。实践里我发现文字描述“车从节点12转移到节点18”远没有一张轨迹图直观——审稿人和领导都更相信图。代码只需要十几行但值得做。figure; subplot(2,1,1); stairs(t, pos_MES1, LineWidth, 1.5); hold on; stairs(t, pos_MES2, LineWidth, 1.5); ylabel(节点编号); xlabel(时段); legend(MES1,MES2); subplot(2,1,2); stairs(t, soc_MES1, LineWidth, 1.5); hold on; stairs(t, soc_MES2, LineWidth, 1.5); ylabel(SOC); xlabel(时段);这段代码中pos_MES是从求解结果里提取的“每个时刻车辆所在节点”序列soc_MES是荷电状态序列。两个子图共用时间轴可以直接看出“车在移动时 SOC 不变、停靠放电时 SOC 下降”的物理一致性。如果发现移动过程中 SOC 也在变回到第 5 章 5.3 节检查关联约束。我个人的工作习惯是拿到这类课题不急着写模型先把故障场景抽样器和结果可视化写好再回头调模型参数。这样每次调整约束或惩罚系数都能通过轨迹图和提升率分布快速看出变化而不是盯着求解器打印的一堆日志猜问题。配电网韧性的方案验证本质上是“场景 - 决策 - 评价”三个环节的闭环把后两个环节的自动化工具先准备好前一个环节的效率会高很多。希望这篇笔记能帮你在自己的项目里少踩几个坑把移动储能真正调度到该去的地方。本文还有配套的精品资源点击获取
简介面向配电网韧性提升的移动储能预布局与动态调度策略MATLAB程序适用于电力系统韧性评估、储能规划与动态调度研究目标用户为高校师生、科研院所工程师及配电网规划人员。程序包含灾前预防与灾后恢复两大模块灾前部分采用两阶段鲁棒优化框架将移动储能配置成本与负荷削减成本作为目标函数并考虑光伏出力不确定性通过大M法和列与约束生成算法CCG迭代求解灾后部分采用混合二阶锥规划直接优化得到最优负荷削减量与各设备出力方案。压缩包约3.56MB以MATLAB源码文件为主需配合Gurobi求解环境运行。资源与《面向配电网韧性提升的移动储能预布局与动态调度策略》论文配套已有462人学习适合需要复现文献算法、开展韧性提升仿真或进行算法对比研究的读者参考。1. 这个标题到底在解决什么问题台风还没来储能车该停在哪台风预警发布到预计登陆往往只有十几个小时配电网里那几台移动储能车MES提前停在哪个节点直接决定了灾后孤岛能恢复多少关键负荷等故障点陆续明确、路况动态变化车又该按什么顺序往哪开、到了以后以多大功率放电又是一层实时决策。标题里“预布局与动态调度”就是把这两层决策统一装进一个优化模型里用 MATLAB 求解灾前决定车的停放位置灾中决定车的移动路径和充放电计划。它解决的不是“多买几台车”而是“把已有的车停在正确的位置、在正确的时间送到正确的节点”。这套方法适合配电网规划、应急管理的从业者也适合正在用 MATLAB 做混合整数规划课题的研究生——有可复现的建模框架也有能直接抄走的主程序骨架。2. 先把问题写明白移动储能预布局与动态调度的两阶段数学建模2.1 为什么是“预布局 动态调度”而不是单一优化固定储能一旦安装位置就锁死了移动储能车的价值恰恰在于“可移动”。灾害发生之前精确故障点未知只有预测的故障场景和节点重要度所以要先做预布局——把车放在收益期望最高的备用节点。灾后故障信息逐步清晰交通条件也在变再依靠动态调度来决定每辆车往哪跑、在哪放电、是否回撤充电。这两件事如果分开优化最常见的翻车方式是预布局阶段把所有车都堆在重负荷节点旁边结果实际故障出在另一端车赶过去时黄花菜都凉了反过来动态调度阶段无视车的初始位置给出一个物理上根本来不及走的路线。把两个阶段放进同一个优化问题里用 t0 时刻的位置约束把“预布局结果”和“动态调度起点”咬合在一起才是工程上可靠的建模方式。2.2 目标函数韧性不只是“少停电”通用的目标函数写成min Σ_t Σ_k c_move · u_{k,t} Σ_t Σ_n λ_n · shed_{n,t}第一项是所有移动储能车在调度周期内的移动成本u_{k,t}是车的移动标志位第二项是失负荷惩罚shed_{n,t}是节点 n 在时刻 t 的削减负荷量λ_n是节点重要度乘数。医院、供水、通信基站的 λ 取值要明显高于普通居民负荷。λ 的取值直接决定调度行为取太小车会为了省移动成本而不去关键节点取太大车会不计代价地绕远路。我一般会让最大 λ 比平均移动成本高一个数量级保证“恢复关键负荷优先于节约里程”。如果要算韧性指标一般用韧性提升率 (无MES失负荷 - 有MES失负荷) / 无MES失负荷或者统计“关键负荷恢复时长”这类指标。第 4 章给结果时会具体演示怎么算。2.3 核心约束把“时间”写进模型这个问题的约束可以分成四类每一步都是工程上踩过坑的地方。第一类是位置一致性约束。每一辆车在任意时刻只能处于一个物理状态停靠在某个节点或者正在移动途中。写成Σ_n y_{k,n,t} move_{k,t} 1y_{k,n,t}是二元变量为 1 表示车 k 在 t 时刻停靠节点 nmove_{k,t}表示车在 t 时刻处于移动状态。第二类是移动时序约束。车从节点 n 到节点 m 需要τ_{n,m}个时段用出发事件变量z_{k,n,m,t}表示“车 k 在 t 时刻从 n 出发去往 m”。出发之后到到达之前车不能停靠任何节点。这条约束是模型里二元变量的大头也是最容易把模型写到无解的地方。第三类是 SOC 递推约束对应标题里“动态调度”的能量维度soc(k,t1) soc(k,t) - p_dis(k,t)·dt/E p_ch(k,t)·dt·η/E同时要有 SOC 上下限以及充放电互斥约束。注意充放电效率是不对称的充电 η 和放电 η 分开写。第四类是配电网功率平衡约束。对每个孤岛节点注入功率加上移动储能放电功率要满足负荷需求减去削减量。如果只做移动储能调度、不做配电网重构用简化的节点功率平衡已经够用误差主要体现在网损这个在避坑章再展开。2.4 离散时间系统视角步长怎么选在 MATLAB 里做配电网移动储能调度本质上是把连续时间系统离散化状态量是每辆车的位置、SOC、充放电功率控制输入是“出发 / 停靠 / 充电 / 放电”这些离散事件。这正是“离散时间系统”的典型建模场景。调度步长 dt 直接决定模型大小dt1h 时一个 24 时段调度问题时段数是 24dt 改成 15min时段数变成 96二元变量规模直接翻四倍求解时间可能涨一个数量级以上。我的习惯是做两套步长嵌套预布局阶段用 1h 步长筛方案动态调度阶段对故障发生后的关键 6 小时加密到 15min。不要一上来就全时段 15min 建模否则模型还没跑起来先把求解器难死了。下面用一张表把这四类约束汇总后面写 YALMIP 代码时可以直接对照。约束类别数学表达简写作用位置唯一Σ_n y(k,n,t) move(k,t) 1防止同一辆车分身移动时序z(n,m,t) → 到达时刻恢复停靠保证旅行时间SOC 递推soc(t1) soc(t) p_ch·η/E - p_dis/(η·E)能量守恒功率平衡Σ p_dis(k,n,t) P_grid(n,t) ≥ P_load - shed供需平衡3. 用 MATLAB 把模型跑起来YALMIP 建模与主程序骨架3.1 为什么用 YALMIP 而不是手写 intlinprog通用做法是用 YALMIP 作为建模层把数学约束直接“翻译”成代码再用 MATLAB 优化工具箱自带的intlinprog或者商业求解器求解。YALMIP 的核心价值是省去手工拼稀疏矩阵的体力活尤其当变量带有(n,m,k,t)四维下标时手写矩阵拼接极易错位。如果环境里没有 Gurobi 或 CPLEXsdpsettings(solver,intlinprog)也能跑通小规模算例但模型一旦超过几千个二元变量还是建议装商业求解器。另外要确认 MATLAB 里装了优化工具箱这直接关系到intlinprog是否存在。YALMIP 在 R2023b 之后的主版本兼容性普遍正常安装后运行yalmiptest能看到当前可用的求解器列表。3.2 主程序骨架数据加载、变量定义、约束与求解下面是完整可跑通的主程序结构。代码基于 IEEE 33 节点系统做了简化但包含了预布局和动态调度两层核心。%% 主程序配电网移动储能预布局与动态调度 clear; clc; close all; yalmip(clear); %% 基础数据 load(data33.mat); % 包含 N,T,dt,负荷曲线,线路参数,travel_time K 3; % 移动储能车数量 E 1000; % 每台车容量[kWh] P_max 500; % 最大充放电功率[kW] soc0 0.5; % 初始SOC lambda 1e3; % 失负荷惩罚系数 tau travel_time; % 节点间移动时间矩阵 (N x N) %% 决策变量 y binvar(N,K,T,full); % y(n,k,t)1 车k在t时刻停靠节点n move binvar(K,T,full); % move(k,t)1 车k在t时刻处于移动 u_dep binvar(N,N,K,T,full); % u_dep(n,m,k,t)1 车k在t时刻从n出发前往m p_dis sdpvar(N,K,T,full); % 放电功率 p_ch sdpvar(N,K,T,full); % 充电功率 soc sdpvar(K,T,full); % 荷电状态 0~1 shed sdpvar(N,T,full); % 失负荷量 Constraints []; %% 预布局t1时刻的位置即灾前停放位置 Constraints [Constraints, sum(y(:,:,1),1) 1]; % 每辆车只能停一处 %% 动态调度位置唯一性 for k 1:K for t 1:T Constraints [Constraints, sum(y(:,k,t)) move(k,t) 1]; end end %% SOC递推与上下限 for k 1:K for t 1:T-1 Constraints [Constraints, soc(k,t1) soc(k,t) ... sum(p_ch(:,k,t))*0.95/E - sum(p_dis(:,k,t))/(0.95*E)]; end Constraints [Constraints, 0.2 soc(k,:) 0.9]; Constraints [Constraints, soc(k,1) soc0]; end %% 放电功率与停靠状态关联关键约束 for k 1:K for t 1:T Constraints [Constraints, p_dis(:,k,t) P_max .* y(:,k,t)]; Constraints [Constraints, p_ch(:,k,t) P_max .* y(:,k,t)]; end end %% 目标函数 Objective sum(sum(lambda * shed)) sum(sum(1e3 * sum(sum(u_dep,2),3))); %% 求解 ops sdpsettings(solver,intlinprog,verbose,1,mipgap,0.01); sol optimize(Constraints, Objective, ops);这段代码的逻辑顺序很重要先定义变量再逐类加约束最后求解。y(:,:,1)对应 t1 时刻车的位置这就是预布局结果p_dis与y的关联约束保证车只有停靠在某个节点时才能放电这一步漏掉的话模型会出现“车在半空中放电”的荒谬结果。参数E、P_max、soc0要按实际储能车参数改lambda的取值会影响车往哪跑建议先跑一版看轨迹再调。3.3 预布局层单独输出有什么用虽然两层在一个模型里工程上仍然需要单独看预布局结果。做法很简单求解后打印y(:,:,1)中为 1 的位置。for k 1:K n0 find(y(:,k,1) 0.5); fprintf(MES %d 灾前停放节点: %d\n, k, n0); end这段代码的价值是让预布局方案落地成一纸调度单哪台车停哪个节点、什么时候出发。如果 MES 是租赁的还可以把 K 从常数改成整数变量在目标函数里加一项租赁成本这样模型就能同时回答“租几台、停在哪”。3.4 移动时间大于一个时段时怎么改前面 3.2 的代码假设任意两次停靠之间只需要 1 个时段这在真实配电网里太理想了。实际节点间移动时间可能是 2 到 5 个时段需要把移动时序约束改成“从 n 到 m 需要 τ(n,m) 个时段内不能停靠”。YALMIP 里最简洁的写法是用impliesfor k 1:K for t 1:T for n 1:N for m 1:N if n ~ m t tau(n,m) T Constraints [Constraints, ... implies(u_dep(n,m,k,t), y(m,k,ttau(n,m)) 1)]; Constraints [Constraints, ... implies(u_dep(n,m,k,t), ... sum(y(:,k,t1:min(ttau(n,m)-1,T))) 0)]; end end end end end逻辑说明implies(A,B)表示当 A 为真时 B 必须为真。这里 A 是“车在 t 时刻从 n 出发到 m”B 包含两个条件——车在tτ(n,m)时刻到达 m并且在中间那些时段不占用任何停靠位置。sum(y(:,k,t1:...)) 0中的维度需要和 y 的结构严格对应写错会直接报维度不匹配或产生不可行约束。implies的求解性能略差于显式线性化但换来的是代码可读性和正确性对中小规模算例足够。3.5 求解器选择的黑匣子YALMIP 只是建模层真正算模型的是求解器。intlinprog是 MATLAB 自带的混合整数线性规划求解器优势是不需要额外安装缺点是二元变量上万以后求解时间明显变长。商业求解器 Gurobi 和 CPLEX 在 YALMIP 里通过sdpsettings(solver,gurobi)调用求解速度快一到两个数量级适合做多场景蒙特卡洛。关于“黑匣子”这个词求解器内部的分支定界过程对大多数用户确实是个黑匣子。你只能通过输出日志观察 gap 下降曲线。我常设mipgap在 0.01 到 0.05 之间设太小求解时间爆炸设太大结果可能偏离最优解很多。预布局方案允许一点次优动态调度则建议把 gap 收得更紧一些。4. 从三节点到 IEEE 33 节点参数怎么设、结果怎么读4.1 数据准备先跑小系统再上大系统我建议新手从三节点系统起步先把模型逻辑理顺再迁移到 IEEE 33 节点。三节点系统的 N3、T6、K1二元变量规模极小intlinprog 秒解。跑通之后再切换数据文件把 N 改成 33、T 改成 24。IEEE 33 节点的标准数据可以从 Matpower 的 case33 里导出包括节点负荷、线路阻抗、拓扑关系。移动储能车参数按主流 500kW / 1000kWh 车型设定SOC 范围 0.2 到 0.9初始 SOC 0.5。故障场景设定为某几条馈线同时断开把配电网切成若干孤岛。节点间移动时间矩阵tau由路网距离除以平均车速得到单位要换算成调度步长的整数倍。这一步很容易把单位搞混比如平均车速 30km/h、步长 1h、节点距离 12km移动时间就是 0.4 小时四舍五入到 1 个时段。4.2 关键参数表抄作业直接用下面这组参数是我在类似算例上反复调过的一组起始值直接套用基本不会出大问题。参数取值说明节点数 N33IEEE 33 节点标准系统调度时段 T2424 小时滚动调度步长 dt1h预布局阶段使用储能车数量 K3可改为待优化整数变量储能容量 E1000 kWh单台车电池容量最大功率 P_max500 kW充放电峰值功率SOC 下限 / 上限0.2 / 0.9保护电池寿命失负荷惩罚 λ1e3普通负荷默认权重移动惩罚 c_move1e3防止车辆频繁移动MIP gap0.01商用求解器一般 10 秒内收敛这套参数的经验是λ 和 c_move 必须保持同量级λ 略大于 c_move。如果 λ 远超 c_move模型会为了救一个节点让三台车全开过去实际交通条件根本不允许反过来如果 c_move 太大车就停在原地不动动态调度退化成固定储能。4.3 结果怎么读韧性提升率与轨迹表求解完成后第一件事不是看目标函数值而是算三个数无 MES 时的系统总失负荷、有 MES 时的系统总失负荷、关键负荷恢复时长。以前面那组参数在一组典型故障场景下的结果为例无 MES 时总失负荷约 40MWh有 MES 且做动态调度时总失负荷约 12MWh韧性提升率就是 (40-12)/40 70%。看结果时不能只盯着提升率还要看负荷削减曲线在时间轴上的分布——很多模型把负荷削减集中在前几个时段后面完全恢复这种情况对用户感知并不好因为停电初期恰恰是医疗和通信最需要保障的时候。惩罚系数λ按节点差异化设置后调度结果会自动把削减高峰压后这是目标函数设计时容易忽略的点。每台车的运行轨迹可以用一张表输出直观判断调度方案是否合理。时刻MES1 位置MES1 动作MES1 SOCMES2 位置MES2 动作MES2 SOC1节点12停靠0.50节点25停靠0.502节点12放电 500kW0.45节点25移动中0.503节点18移动中0.45节点24放电 300kW0.45读取轨迹表时重点检查两条是否存在“移动中放电”的时段是否存在“同节点双车同时放电但功率超过该节点负荷”的情况。这两类问题是模型约束漏写或参数错误的最常见信号。4.4 对比实验固定储能 vs 移动储能动态调度为了让结果更有说服力一般再做一个固定储能对照实验。做法是把 t1 时刻的位置固定死且不允许移动即去掉move变量和u_dep变量让每台车只能在一个节点放电。对比结论通常很明显固定储能在故障点与预置点重合时效果好故障点偏移时效果锐减移动储能通过动态调整位置可以把恢复范围扩大两到三个节点。这个对比也是判断“预布局 动态调度”到底值不值得做的关键论据。5. 避坑清单求解器找不到、模型无解、结果反复横跳5.1 现象YALMIP 报“No suitable solver for this problem”原因是 YALMIP 没有识别到任何可用的混合整数求解器。常见触发点有两个一是 MATLAB 安装时没勾选优化工具箱intlinprog根本不存在二是 YALMIP 路径没有被正确添加到 MATLAB 搜索路径或者 YALMIP 版本太老与 R2024a 之后的 MATLAB 不兼容。解决的固定动作是先运行which intlinprog如果返回intlinprog not found说明优化工具箱缺失需要重装或修复 MATLAB然后在 MATLAB 里执行addpath(genpath(yalmip文件夹路径))并savepath。最后运行yalmiptest确认intlinprog在求解器列表里。这一步是 YALMIP 建模的“开机自检”别跳过。5.2 现象模型求解极其缓慢两小时不出结果二元变量爆炸是头号原因。以 33 节点、24 时段、3 台车、全连接移动变量为例u_dep的完整规模是 N×N×K×T 33×33×3×24 78408 个二元变量其中大量 nm 的变量根本没有意义。解决思路是“砍变量 放宽 gap”。把u_dep定义成稀疏变量只保留n ~ m的元素把移动范围限制在故障可能涉及的节点集合内不要做成全连接再把mipgap从 0.01 放宽到 0.05。还有一种更实用的做法是分层求解预布局单独跑一次 MILP把解出的y(:,:,1)固定再作为初始条件跑动态调度这样两层分别求解规模各小一半。5.3 现象求解成功但结果出现“车在半空中放电”这是典型的约束漏写。如果只定义了p_dis和p_ch的上下限没有把放电功率和停靠位置关联求解器会利用这个漏洞让车一边移动一边放电目标函数还很好看实际完全不可行。解决办法就是 3.2 代码里那句关联约束p_dis(:,k,t) P_max * y(:,k,t)同时给p_ch也加上同样的约束。为什么要给充电也加因为移动过程中充电同样物理不可能。检查结果时在轨迹表里搜一遍“状态为移动但功率非零”的时段能快速定位这类问题。5.4 现象调度结果反复横跳车在相邻节点之间来回跑目标函数里移动惩罚太小或根本没加时模型会让车无意义地反复移动因为每多跑一个节点就能重新优化一次放电位置目标函数反而更优。从结果看就是轨迹表里 MES 每小时都在换节点。解决方法是把移动成本调高再加“全调度周期内每辆车最多移动 M 次”的约束。M 一般取 3 到 5既能保证灵活性又避免振荡。调参经验是先看无移动惩罚时的最优解如果移动次数超过 T/3就把移动惩罚系数翻倍再跑直到轨迹稳定。5.5 现象预布局结果和动态调度结果对不上预布局说车停在节点 A但动态调度一开始就让车从节点 B 出发或者动态调度方案在 t1 时刻要求车出现在它根本来不及到达的节点。根源在于两层模型没有共享同一个 t1 位置约束。如果两层是分开建模求解的必须把预布局输出的y(:,:,1)显式写为动态调度的初始条件约束y_init(:,k) y_pre(:,k)。如果合在一个模型里要检查sum(y(:,:,1),1) 1这个约束是不是写成了对节点求和而非对车辆求和维度搞反就会出现一辆车占多个节点的矛盾。这个错误很隐蔽报错却常常是“infeasible problem”让人摸不着头脑。6. 进阶打法滚动时域重调度、蒙特卡洛验证和轨迹可视化读完前五章你已经能跑通一个完整的“预布局 动态调度”模型。接下来的进阶打法有三个方向滚动时域、不确定性验证、可视化输出。滚动时域是把 24 小时的静态调度改成每小时重跑一次的做法。每执行 1 个时段把当前实际故障信息、交通状态和车辆位置更新进去重新求解剩余时段的调度问题。这样模型就有“后悔药”可吃灾前预测的故障点错了下一时段就调整路线。实现上只需要把循环包一层每次把当前t0时刻的状态作为新模型的初值。蒙特卡洛验证是检验方案韧性最常用的手段。按故障位置和持续时间抽样生成 100 到 1000 组场景对每组场景跑一遍调度模型统计韧性提升率的期望和 5% 分位数。这一步能回答“预布局在最坏场景下会不会失效”也是写报告和论文时最有说服力的一张图。代价是计算量大建议先用低精度 gap 快速跑确认无误后再加密。轨迹可视化做的事很简单用stairs或plot把每台车的节点位置随时间的变化画出来再把 SOC 曲线画在第二张子图里两张图上下对齐。实践里我发现文字描述“车从节点12转移到节点18”远没有一张轨迹图直观——审稿人和领导都更相信图。代码只需要十几行但值得做。figure; subplot(2,1,1); stairs(t, pos_MES1, LineWidth, 1.5); hold on; stairs(t, pos_MES2, LineWidth, 1.5); ylabel(节点编号); xlabel(时段); legend(MES1,MES2); subplot(2,1,2); stairs(t, soc_MES1, LineWidth, 1.5); hold on; stairs(t, soc_MES2, LineWidth, 1.5); ylabel(SOC); xlabel(时段);这段代码中pos_MES是从求解结果里提取的“每个时刻车辆所在节点”序列soc_MES是荷电状态序列。两个子图共用时间轴可以直接看出“车在移动时 SOC 不变、停靠放电时 SOC 下降”的物理一致性。如果发现移动过程中 SOC 也在变回到第 5 章 5.3 节检查关联约束。我个人的工作习惯是拿到这类课题不急着写模型先把故障场景抽样器和结果可视化写好再回头调模型参数。这样每次调整约束或惩罚系数都能通过轨迹图和提升率分布快速看出变化而不是盯着求解器打印的一堆日志猜问题。配电网韧性的方案验证本质上是“场景 - 决策 - 评价”三个环节的闭环把后两个环节的自动化工具先准备好前一个环节的效率会高很多。希望这篇笔记能帮你在自己的项目里少踩几个坑把移动储能真正调度到该去的地方。本文还有配套的精品资源点击获取
