上一期我们把预配置部分讲完了这一期接着来聊MPS动态调度。这篇论文来自配电网韧性方向的一个一区工作整体框架是两阶段鲁棒优化灾害前做应急移动电源的预配置决策灾害后做移动电源的动态调度。上篇解决“电源放哪”下篇解决“电源怎么动、怎么出电”。如果你是刚接触韧性这个方向建议先看一下上篇里的概念铺垫这篇我会把建模细节、求解框架、Matlab实现和调试过程中遇到的那些坑一次性讲透。1. 问题背景与研究思路拆解1.1 配电网韧性问题到底在解决什么先说清楚配电网韧性resilience和传统可靠性的区别。可靠性关注的是日常场景下的小概率短时故障比如某条馈线被施工挖断了一般几分钟到几小时就能恢复。而韧性关注的是极端事件——台风、暴雨、冰灾这类大面积自然灾害往往导致多条线路同时故障、变电站停运停电范围大、持续时间长常规的备用手段根本不够用。这时候应急移动电源Mobile Power SourceMPS就成了关键手段。你可以把它理解成一种“大号充电宝”通常是移动储能车或者应急发电车平时停在某个固定位置灾后由运输车辆把它们运到关键负荷附近临时建立孤岛供电把重要用户先保下来。核心思路就一句话在电网被“打散”成若干个孤岛之后用移动电源去填补那些孤立无援的负荷区。1.2 预配置和动态调度在模型里的分工这篇论文把应急电源调度拆成了两个阶段放在一个两阶段鲁棒优化框架里。预配置是第一阶段的here-and-now决策。在灾害发生之前我们不知道具体哪些线路会断、哪些负荷会失电只有基于气象预报的预测信息。这个时候要先把MPS放到一些预先选好的位置上这些位置必须提前确定因为灾后道路可能被破坏MPS不一定能顺利到达所有地方。预配置做得好不好直接影响后续恢复的灵活度。动态调度是第二阶段的wait-and-see决策。灾害发生后看到实际的网络损坏情况和负荷需求了再决定每个时段里MPS从哪里挪到哪里、接在哪个节点、输出多少功率。这个阶段是滚动式的时间维度从灾后一直延伸到恢复结束。两者不是割裂的。预配置的结果会作为动态调度的初始状态动态调度的可行性又反过来约束预配置的选择。常见的做法是用两阶段鲁棒优化同时建模这两个决策第一阶段变量是预配置位置第二阶段变量是动态调度全过程中间通过不确定集把“最恶劣的灾害场景”嵌进去。1.3 为什么动态调度才是真正难啃的部分预配置本质上是选址问题整数变量集中在“MPS放在哪个节点”决策空间不算大。动态调度则完全不同它把时间维度和空间维度耦合在一起每个时段MPS的位置、充放电功率、接入节点再加上每个时段的网络拓扑状态变量数量指数增长。更难的是移动电源的调度带有强烈的时空耦合约束。MPS从节点A移动到节点B需要时间移动过程中不能供电到了新位置接入电网后又要满足配电网的潮流约束和电压约束。电量也是有限资源白天放了多少电晚上就没得用必须在整个恢复周期内统筹分配。这些约束叠在一起模型复杂度一下子就上来了。我复现的时候最深的一个体会是动态调度这一步的求解效率基本上决定了整篇论文能不能复现成功。模型规模一大求解器直接卡死是常有的事。2. 动态调度模型的核心难点与数学建模2.1 目标函数损失什么才算最小动态调度部分的优化目标论文里最常用的写法是“最小化整个恢复周期内的加权负荷削减量”。权重怎么定不同负荷的重要程度不一样。医院、通信基站、供水系统这些关键负荷权重最高商业负荷次之普通居民负荷权重最低。目标函数的形式大致是min Σ(t) Σ(i) w_i * (P_load(i,t) - P_rec(i,t))其中w_i是节点i的负荷权重P_load(i,t)是t时段节点i的需求P_rec(i,t)是恢复的功率两者之差就是削减量。相当于把每个时段、每个节点的失电损失按权重累加起来。这个目标看起来简单实际上隐含了一个先验假设权重系数怎么设直接决定了调度策略的偏向。有些论文会用持续时间惩罚有些会考虑二次型来惩罚大范围同时削减还有的是直接最大化最小恢复率。复现的时候一定要看清原文目标函数的确切形式不同形式对结果影响很大。2.2 关键约束移动时间、电量、网络潮流动态调度模型里最核心的约束可以分成四组缺一组都会出问题。第一组是MPS的位置与移动约束。用二进制变量代表某个MPS在某个时段是否接入某个节点同时要保证每个MPS每个时段最多只能在一个节点。移动约束是关键中的关键如果MPS在t1时刻从节点A出发经过移动时间travel_time(A,B)之后到达节点B那么在t1到t1travel_time之间的时段MPS不能出现在任何节点。这个约束漏了MPS就会在结果里“闪现”恢复曲线漂亮得一塌糊涂实际根本没法落地。第二组是MPS的电量约束。每辆MPS有电池容量每个时段的输出功率不能超过额定功率电量变化等于放电量如果有分布式电源补充的话还包含充电部分。需要注意的是MPS的剩余电量在每个时段结束时要保持一致不能算着算着电量变成负数或者凭空增加。第三组是配电网潮流约束。这一块通常采用DistFlow线性化模型即把非线性潮流方程在正常工作点附近线性化得到有功、无功平衡方程以及电压降落方程。线性化之后才能把整个模型写成MILP否则非线性项会让求解器无从下手。电压约束一般写成节点电压幅值上下限线路功率约束限制每条支路的传输容量。第四组是孤岛运行约束。MPS接入配电网后恢复区域内的负荷和MPS构成一个孤岛必须满足孤岛内有功、无功功率平衡。如果恢复区域内没有MPS接入那么这个区域无法恢复供电。2.3 不确定集合最恶劣场景怎么刻画两阶段鲁棒模型的核心是“先决策后观察”。预配置在灾害前确定动态调度在灾害后实现。两者之间用一个不确定集合来连接。论文里最常见的是盒式不确定集加预算约束。不确定参数通常是每个节点的负荷需求和线路的故障状态。负荷预测有误差不确定集就用区间刻画比如真实需求在预测值的±20%范围内波动。线路故障状态更直接用0-1变量表示“断”或“不断”不确定集里约束断线数量不能超过某个上限。预算约束的作用是控制保守度。如果不加预算约束意味着所有负荷同时出现最坏情况、所有线路同时断这种极端场景几乎不会发生结果会过度保守MPS会被配置得冗余到离谱。加了预算约束之后不确定集限定了“最多有多少个参数同时偏离最坏值”这才符合实际。复现时这个预算参数是个调参关键。设太小模型过于乐观紧急情况下MPS不够用设太大求解时间爆炸结果保守得没法看。一般做法是从小到大扫描几个值做敏感性分析。3. 复现过程中的关键技术决策3.1 两阶段鲁棒优化求解框架选型处理两阶段鲁棒优化主流方法是CCG列与约束生成也叫Benders分解的加速变体。核心逻辑是把问题拆成主问题和子问题迭代求解。主问题MPMaster Problem包含第一阶段预配置变量和一个辅助变量θ每一轮迭代会添加一组对应特定场景的约束所以主问题规模会越来越大。主问题的最优值给鲁棒问题的下界LB。子问题SPSubproblem给定主问题解出的第一阶段变量在最恶劣的不确定参数下求解第二阶段的动态调度问题。注意这里要最大化第二阶段的目标值这就形成了max-min结构。子问题的最优值给上界UB。迭代过程就是不断“主问题给方案、子问题找茬”直到UB和LB的gap小于设定阈值。这个算法在论文里写起来很漂亮复现的时候麻烦事一堆后面我会专门讲调试的问题。3.2 CCG算法的实现路径子问题的max-min结构是这个算法最费神的地方。内层min是一个动态调度问题里面有MPS位置、接入状态这些整数变量直接做对偶不成立。常见的处理方式有两种。一种是利用强对偶把内层min问题对偶化成max问题这样max-min就变成了max-max可以合成一个单层的max问题但由于内层存在整数变量对偶并不直接可行需要借助KKT条件或者Big-M线性化来处理互补松弛条件。另一种更省事的方式是如果问题规模可控直接把子问题作为一个MILP交给求解器处理让求解器内部的算法来解这个max-min问题不加手动变换。第二种方式对于复现阶段来说更稳妥不容易因为对偶推导出错而卡住。我复现时用的是第二种方案子问题直接写成MILP丢给求解器虽然理论上不如手动对偶推导那么符合论文的“优雅形式”但工程上省事且可靠。如果你的场景规模特别大再考虑手动做KKT变换也不迟。每次调整不确定集参数后记得清理求解器缓存否则后续迭代会带着上一次的结果算容易出错。3.3 参数设置与工具链选择工具链我选的是Matlab YALMIP Gurobi这套组合在电力系统优化复现里最主流。Matlab负责数据处理和结果可视化YALMIP负责建模Gurobi负责求解MILP。如果你手头没有Gurobi许可证用Cplex也行但Gurobi的MILP求解速度在同类问题上实测更稳。关键参数设置上有几个值得注意恢复周期T论文里常见的是24小时或72小时求解时间随T增加几乎是线性放大我建议先用T12小时做调试模型跑通了再扩大到24小时。MPS数量一般设2到4辆。数量越多求解难度指数上升因为MPS之间有协同调度的问题整数变量的组合空间急剧膨胀。时间粒度1小时间隔是常见选择。粒度太细比如15分钟会让移动时间约束变得极其碎片化求解器压力很大粒度太粗比如4小时又会让MPS的调度失去灵活性。Big-M的值这个参数让我吃了不少苦头。M值应该取一个略大于对应量的合理上界的数不要贪图“足够大”太大的M值会引发数值病态。4. Matlab代码实现的核心环节4.1 数据准备与测试系统搭建我选用的测试系统是常见的IEEE 33节点配电网系统。这个系统有33个节点、32条支路带一个主变电站。做韧性分析时会把主变电站视作出口极端场景下失去主网供电系统只能靠MPS孤岛运行。数据准备阶段需要准备的东西整理如下数据项内容来源/说明节点负荷数据每个节点的有功、无功需求标准测试系统参数线路参数电阻、电抗、容量上限标准测试系统参数MPS参数容量、最大输出功率、移动速度论文中给出或自行设定节点间移动时间每对节点间的道路通行时间根据网络拓扑估算负荷权重每个节点的重要等级按负荷类型设定不确定参数负荷波动范围、断线预算基于预测信息设定MPS参数我设置成这样每辆MPS容量为500kWh最大输出功率为100kW平均移动速度取10km/h考虑灾后道路状况打折在33节点系统里相邻节点间移动时间大约0.2到0.5小时所以把移动时间取整到1小时粒度调度上会更清晰。节点负荷数据直接读入Matlab后需要按时间序列扩展。论文通常会给一个日负荷曲线把每个节点的峰值负荷乘以归一化曲线系数得到每个时段的负荷值。这里注意不确定集要在这个基础负荷上叠加扰动扰动幅度我取的是±20%。4.2 YALMIP建模关键代码梳理变量定义这部分是最考验细心的。MPS的位置变量、接入状态变量、出力变量、电量变量还有配电网各节点的电压相角、线路功率每个变量都要小心定义维度。核心变量的定义代码如下% 参数定义 nb 33; % 节点数 nl 32; % 支路数 T 24; % 恢复周期小时 nMPS 2; % MPS数量 MPS_cap 500; % MPS容量 kWh MPS_pmax 100; % MPS最大功率 kW % MPS位置变量 x(m,i,t)MPS m 在 t 时段接入节点 i x binvar(nMPS, nb, T, full); % MPS出力变量 p(m,t)MPS m 在 t 时段的输出功率 p sdpvar(nMPS, T, full); % MPS电量变量 e(m,t)MPS m 在 t 时段结束时的剩余电量 e sdpvar(nMPS, T1, full); % 节点恢复功率变量 prec(i,t) prec sdpvar(nb, T, full); % 线路功率变量 pl(l,t) pl sdpvar(nl, T, full);接下来是约束的添加。位置约束要注意“每个MPS每个时段最多接入一个节点”这个约束用求和形式写% 每个MPS每时段最多接入一个节点 for t 1:T for m 1:nMPS Constraints [Constraints, sum(x(m,:,t)) 1]; end end移动时间约束是这个模型里最容易出错的。假设travel_time(i,j)表示从节点i到节点j的通行时间单位小时如果MPS在t时段在节点i那么它在t1到ttravel_time(i,j)-1之间不能出现在节点j。写成逻辑约束% 移动时间约束从节点i到节点j需要时间tau期间不能出现在j for t 1:T for m 1:nMPS for i 1:nb for j 1:nb if travel_time(i,j) 1 i ~ j tau travel_time(i,j); % 如果t时段在i那么t1到min(T, ttau)时段不能在j for tt t1:min(T, ttau) Constraints [Constraints, x(m,i,t) x(m,j,tt) 1]; end end end end end end电量约束相对直接每个时段结束时的电量等于上一时段电量减去该时段放电量同时电量始终在0到容量之间% 电量平衡约束 for m 1:nMPS Constraints [Constraints, e(m,1) MPS_cap]; % 初始满电 for t 1:T Constraints [Constraints, e(m,t1) e(m,t) - p(m,t)]; Constraints [Constraints, 0 e(m,t1) MPS_cap]; Constraints [Constraints, 0 p(m,t) MPS_pmax]; end end潮流约束我用DistFlow的线性化版本有功和无功都考虑进去。这部分代码量比较大建议做成独立的函数文件方便不同场景复用别把所有代码堆在一个脚本里否则后面调试起来找错会非常痛苦。4.3 结果分析与可视化模型求解完成后输出结果的整理和可视化是验证复现是否成功的重要环节。我常用的输出维度和可视化方式有这么几类第一类是MPS的轨迹图。把每个时段的MPS位置画在配电网拓扑图上用箭头标注移动路径。这能直观看到MPS是否在预配置位置出发、中途如何迁移、最终停在哪里。如果轨迹里有“闪现”立刻能发现问题。第二类是负荷恢复率曲线。横轴是时间纵轴是加权恢复率即当前恢复负荷占权重总负荷的比例。好的结果应该呈现阶梯式上升也就是随着MPS逐步到位恢复率逐步提高而且关键负荷优先恢复。第三类是MPS的出力功率图和SOC曲线。出力功率图能看出MPS是否满负荷运行还是留有裕度SOC曲线则反映电量消耗的平滑程度。如果SOC曲线出现突然跳变多半是电量约束加错了。复现论文最大的一个挑战是结果对不上原文图表。我遇到的典型情况是恢复率曲线趋势一致但绝对数值差了几个百分点。后来发现是负荷权重的设定和原文不一致有些论文把所有负荷权重设成1有些按重要等级分档。这个细节光看文字描述容易忽略一定要仔细对照原文的案例设置。5. 常见问题与排查技巧实录5.1 求解时间爆炸怎么办动态调度模型规模一大Gurobi求解MILP经常跑几小时都出不来。这种时候不要盲目等先检查是不是模型本身有太多整数变量。我遇到一次求解时间爆炸是在T24、MPS3的设置下整数变量数量达到了上万Gurobi报出根松弛就要花很长时间。后来做了三件事解决了第一把MPS位置变量从二维的(m,i,t)压缩成带可行域限制的稀疏变量。有些节点根本不具备MPS接入条件比如没有联络开关或没有可用场地这些节点直接从索引集合里剔除变量数量立刻降下来。第二根据移动时间矩阵提前剪枝。如果节点i到节点j的移动时间大于剩余恢复时间那么这条路径上根本不可能有调度动作对应的约束就不用添加。第三给求解器设置合理的MIP Gap容忍度。复现阶段Gap设到2%完全够用没必要追求0.01%的精确解。Gurobi里直接用GurobiOptions.MIPGap参数控制实测能省掉一大半求解时间。5.2 MPS轨迹紊乱的根因前面提到的“闪现”问题我记得很清楚有一次结果图画出来MPS在t3时在节点12t4就跑到了节点25这两个节点之间移动时间至少要2小时结果模型给出的方案完全违反了物理约束。排查下来发现原因出在移动时间约束的写法上。我最初用的是“x(m,i,t) x(m,j,t1) 1”这只能限制相邻时段的直接跳转没有约束中间时段。正确的写法是把移动时间内的所有时段都加进约束也就是上面代码里那种双重循环。这个错误特别隐蔽因为目标函数值可能看起来没毛病恢复率也正常但MPS的物理轨迹完全不可行。另外还有一个坑移动时间的取整逻辑。如果移动时间是1.5小时而调度时间粒度是1小时那么到底算需要2个时段还是1个时段这里的约定会影响约束的严格性。稳妥起见统一向上取整保证不会低估移动时间。5.3 对偶gap不收敛的处理CCG迭代的核心指标是gap值。我遇到过一种情况gap从15%降到5%之后怎么迭代都降不下去一直卡在5%附近震荡。这种不收敛往往来自子问题的解的稳定性问题。具体来说就是子问题在每次迭代中找到的最恶劣场景对应的最优解变化很小但目标值波动明显导致UB不稳定。几个排查方向第一检查不确定集的预算约束是否合理。如果预算设得太大子问题会在场景组合里反复横跳每次都是不同的极端组合UB自然稳定不下来。把预算缩小一些试试。第二给求解器设置DeterministicTimeLimit确保子问题每次求解都在同样的时间预算内完成避免因为求解器随机分支策略导致每次结果有细微差异。第三检查主问题添加的CCG割是否有重复。有些版本的YALMIP在处理索引变量时会产生重复约束累积到后期会让主问题退化。加一个查重的预处理逻辑能改善收敛速度。5.4 几个值得保留的小经验复现过程中沉淀下来几条经验写在这里供参考。Big-M参数的选择要“刚刚好”。我一开始图省事把M值设成了10^6结果数值病态严重求解器报了一堆数值警告。后来改成按实际物理上界设置有功功率的M取线路容量乘2电压的M取0.1pu量级合适之后求解稳定性明显改善。初始可行解能极大加速求解。用启发式方法先生成一个粗糙的可行调度方案作为MIP start喂给求解器Gurobi可以在几分钟内找到好解而冷启动可能要半小时以上。生成初始可行解的办法是先把MPS固定在预配置位置不移动只做出力优化得到一个保守但可行的方案再放开移动变量做完整求解。不确定性参数的场景生成别贪多。理论上不确定参数组合数是指数级的但实际迭代中只要10到20轮CCG就能收敛到稳定结果。每轮迭代生成的场景可以保存下来做可视化观察哪些场景最恶劣这本身也是对结果的一种合理性验证。还有一个容易被忽略的点求解结束后一定要做可行性校验。把最优解带回去逐个检查原始约束是否满足特别是移动时间约束和电量约束。有时候求解器报告的“最优解”因为数值问题在原始约束下是轻微违规的这一步校验能避免把错误结果写进报告里。最后说一个通用建议无论复现哪篇论文一定不要只盯代码看要把论文中的每个公式编号和代码里的约束一一对应起来。这张“公式-代码映射表”在调试阶段价值极大哪块结果不对顺着映射表定位到对应的约束很快就能找出问题所在。这个项目做完之后我最大的体会是两阶段鲁棒优化的难点不在于理解算法流程而在于所有细节叠加在一起时如何系统地定位和解决那些隐形的建模错误。动态调度部分的MPS轨迹验证就是我排查过程中最有效的一把钥匙。
