1. 为什么预配置算完了还不够MPS动态调度的必要性1.1 配电网韧性评估从“抗灾”到“快速恢复”的核心指标先说个背景。配电网韧性这个概念前几年国内研究热度还没那么高但这几年极端天气频发台风、冰灾、暴雨导致的大规模停电事件大家都有切身感受。传统配电网规划主要解决的是“正常态”下的供电可靠性N-1通过率、供电可用率这些指标但真正面对极端灾害时这些指标是不够用的。韧性Resilience关注的是系统在遭受小概率-高影响事件时能不能扛得住、恢复得快不快量化指标通常包括灾害过程中的失负荷量、恢复时间、系统性能曲线的凹陷面积等。我复现的这篇文章说的是“应急移动电源”也就是MPSMobile Power Source移动储能电源在配电网韧性提升里的作用。这里面有个很关键的区分预配置和动态调度是两个不同阶段的问题。预配置阶段是在灾害预警信息发布之后、灾害来临之前在不知道故障具体位置的情况下决定把MPS预先部署在哪些节点附近待命动态调度则是灾害发生过程中随着故障信息逐步明确、网络拓扑状态发生变化实时决策MPS往哪走、在哪个节点接入、以什么功率充放电。很多刚开始做这个方向的同学会困惑既然预配置阶段已经安排了MPS的位置为什么还要做动态调度这不是重复工作吗其实不是。预配置阶段手里的信息是“预测的灾害强度、可能影响的线路区间”本质是随机优化或者鲁棒优化而动态调度阶段手里的信息是“已经发生的故障、实时的网络连通状态、最新的负荷需求”本质是跟踪实际灾情修正之前的部署计划。两者配合才能把移动储能这种灵活性资源的价值真正释放出来。1.2 预配置与动态调度的分工一阶段“蹲点”二阶段“跑位”如果拿一个类比预配置阶段的MPS就像是消防队根据火警预报提前把消防车开到几个可能发生火灾的小区门口蹲点而动态调度阶段就是在火灾真正发生后指挥中心根据哪个楼栋火势大、哪个路段封路了实时调度消防车从蹲点位置赶赴具体着火点并且决定每辆消防车喷多少水、喷多久。放在配电网场景中预配置阶段一般建模为两阶段鲁棒优化的第一阶段或随机规划的前期决策决策内容是MPS的初始停靠节点位置。动态调度阶段则建模为灾害过程中的时空调度决策变量包括MPS在各时段所处的节点位置、是否移动、在接入节点注入的有功/无功功率、以及荷电状态SOC的演变。两个阶段之间通过MPS的初始位置约束衔接起来动态调度开始时MPS的起点就是预配置阶段给出的待命节点。这种“先蹲点、再跑位”的思想其实抓住了移动储能最核心的优势相比于固定储能MPS可以跨节点转移跑到故障区域附近的变电站或关键负荷节点去供电相比于应急发电车MPS响应速度更快、噪音排放更小、接入更灵活。但与此同时MPS的调度也受交通网络约束限制从一个节点移动到另一个节点需要时间这个时间成本必须在模型里明确刻画否则算出来的调度方案在现实中根本没法定时抵达。1.3 动态调度和普通网络重构、负荷转供的区别还要说清楚一个问题配电网灾后恢复本身已经有大量研究比如网络重构、孤岛划分、分布式电源DG出力调整、需求响应等。MPS动态调度和这些措施的差异在于MPS引入了一类特殊的时空耦合约束——MPS的可用性取决于它的物理位置而位置的变化需要消耗时间。换句话说MPS在t时段能不能在节点i提供功率取决于它在t时刻是否已经抵达或本来就停靠在节点i这一条件就把时间维度耦合进空间决策里了。这一点直接导致模型的求解难度上升了一个档次。普通网络重构问题虽然也是混合整数规划但网络拓扑变量一般是静态的或者准静态的而MPS调度多了位置状态变量、移动时间变量每个时段的状态转移存在时序关系。所以动态调度模型比单纯的故障恢复模型多考虑一整套移动资源约束这也是这部分工作能发在SCI一区的原因之一——问题建模够复杂、求解方法有讲究、应用场景有实际价值。我在复现时的体会是理解清楚动态调度在整个韧性提升框架中的定位比急着打开Matlab敲代码要重要得多。模型定位清楚了代码写起来才不容易乱。2. MPS动态调度的数学建模目标函数与关键约束2.1 目标函数失负荷价值损失怎么量化文献原文的目标函数形式可能有差异但核心逻辑是一致的在灾害影响时段内最小化系统总的失负荷量有时候会带上负荷权重系数反映不同等级负荷比如医院、通信基站、居民负荷的重要性差异。我复现的工作里目标函数如下MatlabYalmip表达objective 0; for t 1:T for i 1:N objective objective w(i) * P_load_demand(i,t) * (1 - x_serve(i,t)) * dt; end end其中w(i)是节点i的负荷权重x_serve(i,t)是0-1变量或者连续变量表示节点i在t时段是否被服务/服务比例P_load_demand(i,t)是该时段的负荷需求。注意这里是“削负荷”还是“恢复供电”的两种建模视角有些论文用“切除负荷量最小化”有些用“恢复负荷量最大化”本质等价但代码实现时变量的含义不同容易出错。更精细的目标函数还会加入MPS移动成本的惩罚项。比如MPS从节点a移动到节点b会产生交通成本虽然对于韧性提升来说这个成本不是最关键的但加上这一项可以有效避免目标值相同的情况下调度的“无意义乱跑”。不过这一项的权重要调得小一些否则模型会牺牲供电恢复去节省移动成本得不偿失。2.2 移动储能路径与功率协同约束MPS动态调度模型最核心的部分是移动储能的时空状态约束。我用整数变量u(i,t)表示MPS在t时段是否连接在节点im(i,j,t)表示t时段MPS正在从i向j移动。这两个变量的关系要通过“每个MPS在同一时刻只能处于一种状态”来约束% 每个时段的单点约束 for k 1:M for t 1:T sum(u(:,t,k)) sum(m(:,:,t,k), all) 1; end end这里的sum(u(:,t,k))是MPS k在t时段接入的节点数之和只能接入0或1个节点sum(m(:,:,t,k), all)是正在移动中的状态标识。还有状态转移约束如果MPS在t时段接入节点i那么t1时段要么继续接入节点i要么开始向某个节点j移动不能瞬移。移动时间的建模也很有讲究。从节点i到节点j需要travel_time(i,j)个时段这个时间矩阵一般根据交通网络距离除以MPS运输速度得到。在代码里我直接用了一个预计算的矩阵travel_time(i,j)然后用“延迟到达”方式来约束for t 1:T-1 for i 1:N for j 1:N if travel_time(i,j) T-t % 如果t时段开始从i向j移动那么到达时刻不能早于ttravel_time(i,j) m(i,j,t) 1 - u(j, ttravel_time(i,j)); end end end end这只是其中一种实现的方式实际代码需要配合状态变量一起把逻辑闭合。MPS在节点接入后充放电功率要满足上下限约束同时SOC状态转移SOC(k,t1) SOC(k,t) (eta_c * P_ch(k,t) - P_dis(k,t) / eta_d) * dt / Cap(k);其中eta_c、eta_d是充放电效率Cap(k)是储能容量。注意充放电功率不能同时非零需要引入互补约束或者0-1变量辅助。2.3 配电网运行约束潮流方程与安全边界动态调度不是孤立的储能优化而是储能和电网的联合优化。配电网潮流约束采用DistFlow形式比较常见。对每个节点i和时段t% 有功平衡 P_inj(i,t) P_g(i,t) P_dis_sum(i,t) - P_ch_sum(i,t) - P_load_served(i,t); % 其中P_dis_sum和P_ch_sum是所有MPS在该节点注入/充电的聚合功率支路潮流约束线性化DistFlowP_branch(l,t) sum(P_inj(in_nodes_before_l, t)); % 简化表达 P_branch(l,t) P_line_max(l);节点电压约束V(i,t) V_min; V(i,t) V_max;想要处理快、收敛稳电压约束通常用线性化的DistFlow忽略网损项或者用线性网损近似。SCI论文里一般会用二阶锥松弛SOCP把DistFlow精确化Matlab里用Yalmip可以很方便地写锥约束但求解时间会显著增加。如果只是复现一个中等规模算例比如33节点系统SOCP压力不大。2.4 动态调度与预配置的时序衔接动态调度模型的初始条件就是预配置阶段给出的MPS初始位置。如果只做动态调度模块这一部分通常简化为“给定初始位置”参数。我在代码里用一个矩阵initial_position(k, i)来存储预配置结果然后在约束中把t1时刻的接入位置固定for k 1:M for i 1:N if initial_position(k,i) 1 u(i,1,k) 1; % t1时MPS停在预配置节点 end end end如果你要复现整篇论文预配置动态调度联合仿真那就是两层优化迭代外层生成预配置方案内层跑动态调度用动态调度的结果去评估外层方案的韧性指标。这样理解就会清晰很多。3. 模型怎么才能算得动解耦策略与求解器选型3.1 为什么不能直接丢给求解器一个“大杂烩”模型动态调度模型里有0-1变量接入位置、移动状态、充放电状态还有连续变量功率、SOC、电压是个典型的混合整数规划MIP。如果系统规模不大比如IEEE 33节点、2台MPS、24个时段直接用YalmipGurobi/Cplex求解是可行的。但如果规模放大到100节点、4~6台MPS、96个时段那MIP的求解时间会指数级增长直接求解很容易卡死。SCI论文里通常会做一些数学变换来加速求解。最常见的做法是线性化——把双线性项比如0-1变量乘以连续变量用大M法线性化或者用McCormick包络。这一部分代码写得好不好直接影响求解速度。3.2 大M法处理逻辑约束以充放电互斥为例MPS在同一节点同一时段不能同时充放电这个逻辑约束可以直接写为P_ch(i,t) M * z_ch(i,t); P_dis(i,t) M * z_dis(i,t); z_ch(i,t) z_dis(i,t) 1;z_ch和z_dis是0-1变量M取一个足够大的正数一般取储能额定功率的1.05~1.2倍就够了不要取10^6这种“大得吓人”的数。我见过很多同学在写大M约束时随便取M1e6结果导致数值条件变差求解器迭代半天出不来还以为模型本身有问题。M的取值原则是“比物理最大值稍大即可”。3.3 求解器参数与性能表现Gurobi实测数据我在复现时对比了几种求解器的表现用的算例是IEEE 33节点系统3台MPS24个时段场景是某一线路发生永久性故障对比结果如下表求解器求解时间目标值总失负荷量/kWh说明Gurobi 11.08.4s1276.5默认参数MIP Gap 0.01%Cplex 12.1015.2s1276.5略慢结果一致GLPK时间过长未收敛—开源求解器处理MIP明显吃力SCIP46.7s1276.5能算但速度一般Gurobi在解决这类配电网MIP问题上是目前表现最好的强烈建议用MatlabYalmipGurobi的组合。如果实验室没有Gurobi授权可以先申请学术版免费但需要提交申请或者用Cplex。开源方案里SCIP是勉强能用的但体验差距明显。求解器参数层面我习惯设置outputflag 0关闭求解日志如果调试阶段则打开mipgap 1e-3或更小。Yalmip里设置方式如下ops sdpsettings(solver, gurobi, verbose, 1, gurobi.mipgap, 1e-3);3.4 滚动时域控制退化大规模场景的另一个选择如果调度时域很长比如96时段一次性求解整段MIP很可能内存爆掉。这时候可以用滚动时域控制的思路把调度问题按时间窗口切为多段子问题每次只优化未来H个时段执行第一个时段的决策然后滚动推进。这种方式虽然不是全局最优但抗不确定性更强——因为未来的故障信息本来就是逐步明确的。滚动时域和“全局优化”的取舍SCI论文里一般会做对比。我复现时也发现全局优化在某些场景下目标值确实好一些但如果把故障信息的不确定性引进来滚动时域的解反而更接近实际可执行方案。如果你的目标不是复现答辩而是实际科研出图建议两种都跑一下用对比曲线来支撑结论。4. Matlab代码实现从矩阵搭建到结果可视化4.1 数据结构设计节点、线路、MPS三类对象的组织Matlab写优化模型最忌讳的就是变量满天飞没有结构。我强烈建议用结构体来组织数据。下面是我复现时采用的框架%% 基础数据 mpc loadcase(case33bw); % IEEE 33节点配电网数据 N 33; % 节点数 L 37; % 支路数 T 24; % 调度时段1h间隔 %% 负荷数据从原始数据提取 load_demand mpc.bus(:,3) / 1000; % 单位转为MW P_load repmat(load_demand, 1, T) .* load_profile; % load_profile是24h标幺曲线 %% MPS参数 M 3; % MPS数量 mp_cap [500; 500; 1000]; % 容量/kWh mp_pmax [100; 100; 200]; % 最大充放电功率/kW mp_s0 [0.3; 0.5; 0.8]; % 初始SOC节点和线路数据直接用Matpower格式读取即可但要注意电压等级和功率基准值的换算。如果基准值是100kVA、10kV那么负荷数据就要统一换算成标幺值或者同一单位体系否则目标函数和约束里的数量级会出问题。这个单位换算问题我踩过坑一开始直接拿原始kW数据写约束结果潮流约束里支路功率和节点注入功率差了三个数量级求解器疯狂报数值错误。4.2 核心代码模块拆解目标函数、约束生成、求解调用我整个代码分成四个文件逻辑非常清楚main_schedule.m主文件初始化数据、调用建模、求解、输出结果build_network_matrix.m生成节点关联矩阵、支路首末端矩阵、DistFlow系数矩阵build_mps_constraints.m生成MPS移动、容量、SOC相关约束plot_results.m可视化韧性恢复曲线建模的主体用Yalmip表达。下面的代码展示了核心约束的生成方式%% 决策变量 u binvar(N, T, M, full); % MPS k在t时段是否在节点i m binvar(N, N, T, M, full); % MPS k在t时段是否正在从i向j移动 P_ch sdpvar(N, T, M, full); % 充电功率 P_dis sdpvar(N, T, M, full); % 放电功率 SOC sdpvar(M, T, full); % 荷电状态 z_ch binvar(N, T, M, full); % 充电标志 z_dis binvar(N, T, M, full); % 放电标志 Constraints []; %% 唯一接入约束 for t 1:T for k 1:M Constraints [Constraints, sum(u(:,t,k)) 1]; Constraints [Constraints, sum(m(:,:,t,k), all) 1 - sum(u(:,t,k))]; end end %% 移动时间约束简化为固定移动时间的状态转移 for t 1:T-1 for k 1:M for i 1:N for j 1:N if travel_time(i,j) T-t Constraints [Constraints, ... m(i,j,t,k) 1 - u(j, ttravel_time(i,j), k)]; end end end end end这里要注意binvar(N, T, M, full)创建的是一个N×T×M的三维二进制变量数组在Yalmip里的索引方式要小心别把维度搞混了。我刚开始写的时候经常把u(i,t,k)和u(t,i,k)搞混调试时浪费了不少时间。建议在代码开头先加一句注释或输出验证维度和物理含义。潮流约束用DistFlow线性化表达对每个时段和每条支路定义注入功率与支路功率的关系%% 潮流约束线性DistFlow for t 1:T for i 1:N P_inj(i,t) P_g(i,t) sum(P_dis(i,t,:)) - sum(P_ch(i,t,:)) - P_serve(i,t); end for l 1:L from branch_from(l); to branch_to(l); Constraints [Constraints, P_branch(l,t) P_inj(from,t) ...]; Constraints [Constraints, -P_line_max(l) P_branch(l,t) P_line_max(l)]; end end从节点i注入功率到支路潮流的对应关系需要根据网络拓扑的关联矩阵来换算。这一部分可以直接用Matpower的makeBdc或者自己写一个incidence matrix来生成转换矩阵。我一般是自己写因为SCI复现时常常要修改网络拓扑结构断线、重构用现成函数反而麻烦。4.3 求解与结果输出Yalmip调用Gurobi的完整流程求解部分代码很短但有几个坑值得注意%% 求解 ops sdpsettings(solver,gurobi,verbose,1,gurobi.mipgap,1e-3,gurobi.timelimit,300); result optimize(Constraints, objective, ops); if result.problem 0 disp(求解成功); else disp([求解失败: , result.info]); return; end如果求解失败最常见的错误提示是“Infeasible problem”。这时候别慌按顺序排查先去掉目标函数改成可行性问题feasibility problem看约束是否本身就有冲突再检查MPS初始SOC设置是否合理——如果初始SOC才10%然后要求它输出额定功率6小时那SOC下限约束基本必冲突最后检查MPS移动时间约束里时间索引是否越界。4.4 调度结果可视化韧性曲线怎么画才好看且有说服力可视化是论文复现里占用时间比较长、但最容易被忽视价值的一环。好的图能直接把“动态调度优于无调度”“预配置动态调度优于只做预配置”这些结论一目了然地展示出来。我复现时画了四张图第一张是系统韧性曲线横轴为时间纵轴为“当前可用负荷/总负荷”的百分比。图上画三条曲线无MPS调度、只有预配置MPS原地不动、预配置动态调度。从曲线上你能直观看到动态调度能显著缩短失负荷状态的持续时间并把最低负荷水平抬高。figure; plot(t_hours, available_ratio_no_mps, r--, LineWidth, 1.5); hold on; plot(t_hours, available_ratio_pre, b-., LineWidth, 1.5); plot(t_hours, available_ratio_dynamic, g-, LineWidth, 2); legend({无MPS, 仅预配置, 预配置动态调度}, Location, SouthEast); xlabel(时间/h); ylabel(可用负荷比例); grid on; axis([0 24 0.6 1.05]);第二张是MPS的时空轨迹图用Matlab的imagesc画MPS在各时段接入的节点编号颜色深浅表示接入状态。这张图能直观展示MPS从预配置位置向故障区域移动的过程是审稿人非常喜欢的一种呈现方式。第三张是节点负荷恢复时序图选择几个关键负荷节点用堆叠面积图展示各时段的恢复情况。第四张是SOC变化曲线配合MPS在不同节点的充放电功率柱状图验证调度结果的物理合理性。% MPS轨迹图示例 figure; imagesc(t_hours, 1:N, squeeze(u_solution(:, :, 1))); colorbar; xlabel(时间/h); ylabel(节点编号); title(MPS-1 时空轨迹);5. 复现过程中的坑与调参总结5.1 移动时间矩阵最容易忽略的物理约束很多复现者把移动时间设为固定的常数比如任意两个节点之间都是1小时这在简单算例里还能跑通但如果换到IEEE 123节点这类规模大的系统就会出现严重问题——两个相距几十公里的节点移动时间怎么可能和相邻节点一样物理上完全不合理审稿人一眼就能看出来。正确的做法是根据配电网的地理信息或者你自己构造的节点坐标计算交通距离再除以MPS的运输速度。举个简单的例子假设节点坐标xy(i,:) [x_i, y_i]MPS运输车的平均速度v 40 km/h节点间直线距离乘以绕行系数1.5实际道路不可能是直线那么for i 1:N for j 1:N dist_km norm(xy(i,:) - xy(j,:)); travel_time(i,j) ceil(dist_km * 1.5 / v); % 向上取整单位h end endceil向上取整是必要的因为真实系统不会允许“0.3小时”这种碎片化的移动时间通常以小时或半小时为单位做离散调度。5.2 预配置和动态调度迭代时结果不自洽的问题如果你要跑整篇论文的联合仿真预配置外层动态调度内层最容易遇到的问题就是外层给出一组预配置位置内层动态调度跑出来的目标值和外层评估时用的目标值不一致。这个问题我在复现过程中也遇到了排查了半天发现是外层鲁棒优化用的故障场景集合和内层动态调度用的实际故障场景不一致导致的。预配置阶段面对的是不确定的故障场景会采用鲁棒优化或者场景法来构建一个对所有可能故障都“不至于太差”的方案而动态调度阶段面对的是“实际发生”的故障场景。这本来就是两个信息层级不一致是正常的。但如果你发现连“已知实际故障场景”下的预配置位置都能让动态调度结果变差那就要检查是不是MPS初始SOC、位置约束写错了。5.3 Yalmip和Gurobi的版本兼容性Matlab的版本更新太快Yalmip和Gurobi的接口有时会抽风。我用的配置是Matlab R2022a Yalmip 2023版 Gurobi 11.0。如果装的是R2023b及以上注意Yalmip是否更新到支持该版本的版本否则会出现无法识别求解器的问题。另外Gurobi在Linux服务器上跑通常比Windows上稳定内存管理也好一些。如果算例规模大建议部署在服务器上跑。服务器上装Gurobi有一步容易忽略要手动把Gurobi的matlab接口路径加入Matlab路径否则Yalmip找不到求解器。% 在Matlab中手动添加Gurobi接口路径每次启动可能需要 addpath(/opt/gurobi1101/linux64/matlab); savepath; % 保存到默认路径避免下次重复添加如果提示找不到gurobi_mex检查一下Gurobi的许可文件是否配置好了用gurobi_setup或gurobi_test做环境检测。5.4 参数敏感性分析SOC上下限和移动速度的影响我复现过程中对比了几组参数的影响分享几个结论供你参考参数变化对结果的影响SOC下限从0.1提高到0.3有效放电容量减少系统失负荷量增加约8%~15%MPS数量从2台增至4台失负荷量下降40%以上但求解时间从5秒涨到30秒移动速度从30km/h提到60km/h失负荷量下降约5%因为MPS能更快到达故障区域负荷权重医院负荷×5居民×1调度结果会优先恢复高权重节点总失负荷量可能稍增但韧性价值更高这些敏感性分析对于论文讨论部分特别有用。建议在复现代码时就把这些参数设计成外部可调参数而不是硬编码在代码里——否则每次改参数都要翻代码找位置效率极低。我自己在实际复现过程中还有几个小技巧调试阶段别一上来就跑24时段、3台MPS的完整算例。先跑一个4时段、1台MPS的“玩具模型”确认逻辑正确后再逐步放大规模。这种方法帮我至少节省了两天调试时间。每次修改模型或参数先把运行结果自动保存为mat文件这样后面画图和分析可以直接加载不需要重复求解。用tic/toc把求解时间记录到日志里方便做求解性能对比。说到最后MPS动态调度这个方向最大的魅力在于它把电网调度和交通网络、储能运行耦合在一起建模空间很大做出来的结果也直观——一条韧性曲线预配置和动态调度的优势一目了然。复现过程中踩的坑也是一种成长如果你在复现时卡在某一步不妨照着上面的思路一步步排查大部分问题都会迎刃而解。
