1. 为什么MPS预配置值得复现问题拆解与研究动机前几天一个师弟过来找我说在看一篇关于配电网韧性提升的SCI一区论文标题里提到应急移动电源预配置和动态调度。他问我的第一句话是移动电源不就是储能吗这能发一区论文我当时就笑了。这个问题其实代表了大多数人初看到这个题目时的反应。先说结论**应急移动电源Mobile Power Source, MPS和固定储能完全是两个物种。**固定储能装在哪就在哪放电而MPS的核心价值在于移动二字——它可以根据灾害演进和电网受损情况在空间和时间两个维度上动态调整接入位置把宝贵的电量送到最需要的地方去。单凭移动这个自由度问题的建模复杂度和求解难度就会呈指数级上升这也是它能够撑得起一篇高水平论文的根本原因。这篇论文的工作可以拆成两大部分MPS预配置Pre-positioning在灾害来临前决定移动电源的配置数量、候选站点选址和初始容量分配。这是一个典型的两阶段或三阶段鲁棒优化问题预配置阶段的决策必须在不确定性尚未实现时做出。动态调度Dynamic Dispatch灾害发生后根据线路故障、负荷水平和道路通行状况实时决策移动电源何时、从何处、移动到何处、以什么功率放电。这是一个需要滚动优化的时间-空间耦合决策问题。这篇文章是上篇聚焦第一部分——MPS预配置。从整个复现流程来讲把预配置这个上层决策搞扎实后续动态调度才有意义。因为你在灾前把MPS放在什么位置、分配多少容量直接决定了灾后调度时的可行域上限。用大白话说位置选不对调度再聪明也白搭。后台开发代码属于硬功夫活整个复现的链路可以概括为配电网拓扑与参数处理 → 不确定集建模 → 两阶段鲁棒优化模型的建立 → CCG列与约束生成算法实现 → 结果可视化与敏感性分析。下文我会按照这条链路把这套流程掰开揉碎讲清楚同时给出关键的Matlab代码片段和避坑经验。1.1 预配置问题的本质不确定性中做决策预配置处在整个MPS决策链条的最前端它最大的特点是先决策、后实现。灾害预测在时间尺度上存在天然的模糊性——我们可能知道某个区域大概率会受灾但很难准确预测每一条线路会不会受损。极端天气下的负荷水平也充满随机性——可能是几小时的持续高温推高空调负荷也可能是工厂临时停产导致负荷骤降。在这种双重不确定环境下预配置决策必须回答三个问题买多少台MPS数量太少灾后覆盖不足数量过多闲置成本巨大。放在哪里候选站点可能涉及变电站、开闭所、重要用户附近的空地等。位置决定了灾后MPS到达受损区域的响应时间。每台配多少能量MPS本质上是带电池的移动储能车单台容量有限通常在几百千瓦时到几兆瓦时之间必须根据预判的灾情严重程度进行能量分配。从数学角度看这是一个典型的不确定优化问题。论文通常采用**两阶段鲁棒优化Two-stage Robust Optimization**来建模第一阶段是预配置决策here-and-now第二阶段是预配置完成后、不确定参数逐步显现时的运行决策wait-and-see。这个先决策后兜底的结构用博弈论的视角看就是**系统决策者先走一步然后大自然这个对手选择最恶劣的场景最后系统再针对该场景进行机组调度和负荷削减。**目标是在最坏场景下仍能保证系统安全运行同时总成本投资运行负荷中断损失最低。1.2 为什么选择鲁棒优化而非随机优化在鲁棒优化和随机优化之间做选择本质上是一个保守度 vs 计算复杂度的权衡问题。随机优化假设不确定参数的概率分布已知目标是最小化期望成本。这个思路在分布信息充足时很合理但对于灾害场景存在明显的局限——极端天气下的损失函数通常是heavy-tailed的期望值并不能完全反映灾难性场景的风险。鲁棒优化则恰恰相反它只关心不确定集合内的最坏情况不需要假设任何概率分布。这在灾害场景下非常有吸引力因为极端事件的概率分布很难准确刻画你很难说某条线路断线概率是3.2%还是5.8%决策者更关心在最坏情况下我能不能扛得住而不是平均情况下我损失多少配电网作为基础设施承担着民生保障功能容错空间极小。代价就是模型更保守、求解更难。特别是当不确定集合包含线路N-k故障和负荷水平双重不确定性时模型会变成一个混合整数线性规划MILP嵌入max-min结构的复杂优化问题无法直接求解。这就引出了预配置求解的核心算法——CCGColumn and Constraint Generation。这个算法的思路和动态规划有相似的地方不一次性求解整个大问题而是通过主问题-子问题的迭代交互逐步逼近最优解。2. 预配置模型的数学内核变量、目标函数与约束的逐层拆解复现一篇论文的第一步永远是把数学模型看懂到能够手推的程度。不要一上来就打开Matlab碰代码那不叫复现那叫抄作业。抄完之后遇到问题debug都无从下手。下面我把MPS预配置的数学模型按模块拆开讲。2.1 决策变量体系整个预配置模型包含三类变量第一阶段预配置决策变量(x_{i,m} \in {0,1})节点 (i) 处是否布置第 (m) 台MPS这是选址的0-1变量(E_{i,m}^{pre} \geq 0)布置在节点 (i) 的第 (m) 台MPS的预配置容量kWh这是定容的连续变量(P_{i,m}^{rate} \geq 0)该MPS的额定充放电功率kW。这里需要特别注意**预配置决策变量在灾害发生前就必须确定不能依赖任何不确定参数的实现值。**这是鲁棒优化中非预期性约束non-anticipativity constraint的核心要求——你可以根据预测数据做决策但决策本身不能写成不确定参数的函数。第二阶段运行决策变量(P_{i,t}^{DG})分布式电源在时段 (t) 的有功出力(Q_{i,t}^{DG})分布式电源的无功出力(P_{i,t}^{MPS})MPS在节点 (i) 时段 (t) 的放电功率(P_{i,t}^{cur})节点 (i) 时段 (t) 的负荷削减量(U_{i,t}, \theta_{i,t})节点电压幅值和相角(P_{ij,t}^{flow}, Q_{ij,t}^{flow})支路 (ij) 在时段 (t) 的有功、无功潮流。第二阶段变量是在不确定参数哪个场景实现已知后做出的适应性决策它需要随场景变化而调整。这正是两阶段鲁棒优化的两阶段在数学上的体现。2.2 目标函数的设计逻辑预配置的目标函数通常是一个总期望/总最坏成本最小化问题形式为[ \min_{x, E^{pre}} \left{ \sum_{i \in \Omega_B} \sum_{m \in \Omega_M} \left( c_{inv,i} x_{i,m} c_{cap} E_{i,m}^{pre} \right) \max_{d \in \mathbb{D}} \min_{y \in \mathcal{F}(x, d)} c_{ope}(y) \right} ]拆开来看这三层结构最外层(\min)最小化投资成本购置成本容量成本与最坏运行成本之和中间层(\max)大自然选择不确定性集合 (\mathbb{D}) 内最恶劣的场景内层(\min)系统在该场景下进行最优调度使得运行成本负荷削减成本最小。运行成本部分通常包括[ c_{ope}(y) \sum_{t1}^{T} \sum_{i \in \Omega_B} \left( c_{DG} P_{i,t}^{DG} c_{cur} P_{i,t}^{cur} \right) ]其中 (c_{cur})负荷削减的单位代价是一个非常重要的参数。它反映了系统对停电损失的容忍程度。论文中通常把它设得远高于发电成本比如10倍以上因为负荷中断的社会经济损失远远超过了电费本身。还需要考虑权重因子 (\omega_i)表示节点 (i) 的负荷重要性权重。如果第 (i) 个节点连接的是医院、通信基站、供水系统等重要用户它的 (\omega_i) 会远高于普通居民负荷。2.3 约束条件的立体拆解约束条件是模型的骨架也是复现时最容易出错的地方。我把约束分为五组1MPS购置与配置约束[ \sum_{i \in \Omega_B} x_{i,m} \leq 1, \quad \forall m ] [ \sum_{i \in \Omega_B} \sum_{m \in \Omega_M} x_{i,m} \leq N^{MPS,max} ] [ 0 \leq E_{i,m}^{pre} \leq E^{max} x_{i,m} ]第一组约束表示每台MPS最多只能配置在一个节点第二组是总数量上限约束第三组是容量-选址逻辑约束——如果节点 (i) 没有放置MPS即 (x_{i,m}0)那么该位置的预配置容量必须为0。这是逻辑一致性约束在代码中通过大M或直接乘积实现。2DistFlow潮流方程约束[ P_{ij,t} P_{jk,t} P_{j,t}^{net},\quad \forall (i,j) \in \Omega_L ] [ Q_{ij,t} Q_{jk,t} Q_{j,t}^{net},\quad \forall (i,j) \in \Omega_L ] [ P_{i,t}^{net} P_{i,t}^{DG} P_{i,t}^{MPS} - P_{i,t}^{load} P_{i,t}^{cur} ] [ Q_{i,t}^{net} Q_{i,t}^{DG} - Q_{i,t}^{load} ] [ U_{j,t} U_{i,t} - 2\left( r_{ij} P_{ij,t} x_{ij} Q_{ij,t} \right) \left( r_{ij}^2 x_{ij}^2 \right) I_{ij,t} ]这套方程是配电网潮流计算的标准形式被称为DistFlow方程。它和传统牛拉法潮流计算最大的区别在于DistFlow用支路潮流的递推关系描述网络避开了复杂的矩阵求逆直接以节点电压和支路功率为变量建模天然适合嵌入优化模型。在代码实现时为了避免非线性项 (U_{i,t} I_{ij,t}) 带来的求解难度论文通常做如下线性化处理忽略线路对地导纳配电网线电压低、线路短对地导纳极小电压幅值接近1pu用 (U_{i,t} \approx u_{i,t})电压幅值的二次项近似支路电流 (I_{ij,t}) 用负荷已知的基准值近似或直接采用线性DistFlow模型。3电压安全约束[ U^{min} \leq U_{i,t} \leq U^{max},\quad \forall i \in \Omega_B, \forall t ]这条约束保证了所有节点电压在正常运行范围内。配电网一般要求电压偏差不超过±5%即0.95~1.05 pu。灾害状态下论文通常会适度放宽下限比如到0.9pu以体现韧性场景下的尽力而为逻辑。这个细节直接影响MPS配置结果——电压下限放得越宽系统需要配置的MPS容量就越少。4MPS容量-功率耦合约束[ 0 \leq P_{i,m}^{MPS} \leq P^{rate}, \quad \forall i,m ] [ E_{i,m}^{pre} \eta \Delta t P_{i,m}^{ch} - \frac{1}{\eta} \Delta t P_{i,m}^{dis} \geq E^{min} ]MPS本质上是一个带电池的移动储能它的功率输出受到额定功率限制能量变化受到充放电效率约束。注意这里的充放电效率 (\eta)是不可忽略的参数——锂电池的充电效率一般在90%~95%之间放电效率稍低。如果复现时光考虑能量平衡而不考虑效率损失你算出来的MPS容量会偏小10%以上在极端工况下可能导致负荷失电。5不确定性约束最关键的部分[ d_{i,t} \in \mathbb{D} : \left{ d_{i,t} : \sum_{t} \sum_{i} \frac{|d_{i,t} - d_{i,t}^{0}|}{\hat{d}{i,t}} \leq \Gamma, \quad d{i,t}^{min} \leq d_{i,t} \leq d_{i,t}^{max} \right} ]这里 (\Gamma) 是不确定预算budget of uncertainty控制鲁棒优化的保守程度。(\Gamma) 越大不确定集合越大最坏场景越极端求解结果越保守MPS的预配置数量就越多。理解 (\Gamma) 的意义非常关键。如果不设 (\Gamma)大自然可以选择所有节点在所有时段都取最大负荷的极端情况这在现实中几乎不可能发生。设置 (\Gamma) 之后不确定集合限制偏离预测值的总量不能超过某个水平有效控制了保守度。2.4 从模型到代码的简化路径直接求解上面这个混合整数非线性规划MINLP会消耗极大的计算资源高维场景下甚至无法在可接受时间内得到可行解。为了在Matlab中高效复现论文需要做三个关键简化原始特征处理方法影响DistFlow非线性项线性化忽略对地导纳电压幅值平方项近似模型退化为MILP可直接用Gurobi/CPLEX求解连续不确定集盒式预算约束离散化最坏场景在顶点取得可枚举max-min嵌套结构CCG算法解耦主问题MIP子问题LP迭代收敛这三个简化是工程可解性和理论严谨性之间的折中方案。复现时不需要质疑它但要理解每个简化带来的误差方向和量级。完整的目标函数和约束语句在论文中通常以Annex的形式给出篇幅能达到两三页。**我强烈建议你在写代码前把每条数学约束对应的公式编号、MathType里面的希腊字母含义、矩阵维度都手写一遍。**这一步看似费时但在debug时效率能提升三倍以上。3. 求值算法的核心思想与CCG迭代实现3.1 两阶段鲁棒优化为什么不能直接求我现在用个生活化的方式解释一下。假设你是一家大型活动的安保负责人现在要决定在各个入口安排多少保安相当于MPS的选址和定容。你面临的不确定性是哪些入口会出现拥挤相当于故障场景。如果你一步到位计算最优保安配置那就必须同时回答两个问题万一所有入口都爆满我的保安配置够不够但所有入口同时爆满的概率又极低我有没有可能配置过度直接求解这个问题等价于让计算机同时处理每个入口的保安数量和每种人流场景下的应急处置方案这两个维度的变量。当场景数量呈指数爆炸时计算机的内存和求解时间会双双失控。CCG的聪明之处在于不一次性考虑所有可能的场景而是每次只聚焦于当前最恶劣的场景通过迭代逐步逼近真正的worst-case。用保安的例子说就是先假设人流正常名义场景得出一个基础的保安配置方案针对这个配置让大自然去找茬——找出当前配置下压力最大的入口和时段把找茬找到的最恶劣场景加入模型作为新的约束迫使系统重新优化保安配置重复2和3直到最终配置收敛到不能再改进。这个找茬的过程就是子问题求max的环节把最恶劣场景对应的约束加回主问题就是生成列Column Generation 新增约束Constraint Generation的名字来源。3.2 CCG主问题与子问题的形式化描述主问题MP第k次迭代时[ \min ; c_{inv}^T x c_{cap}^T E^{pre} \alpha ] [ s.t. ; \alpha \geq c_{ope}^T y^{(l)}, \quad \forall l 1, 2, ..., k ] [ x \in \mathbb{X}, \quad y^{(l)} \in \mathcal{F}(x, d^{(l)}), \quad \forall l 1, ..., k ]其中 (d^{(l)}) 是前 (l) 次迭代中次第识别出的最恶劣场景(y^{(l)}) 是对应该场景的调度决策变量。注意 (\alpha) 是一个辅助变量代表最坏运行成本的估计值。当k很小比如只有1-2个场景时主问题规模很小求解速度飞快。子问题SP给定主问题解 (x^, E^{pre,})[ \max_{d \in \mathbb{D}} ; \min_{y \in \mathcal{F}(x^, E^{pre,}, d)} c_{ope}^T y ]子问题是一个max-min结构不能直接用商业求解器处理。处理方法通常有两种强对偶法将内层min问题取对偶将max-min转化为max-max即单层max问题KKT条件法将内层min的KKT条件作为约束加入外层max转化为带有互补松弛条件的单层问题。在Matlab实践中我最推荐的是强对偶法。因为线性规划的强对偶成立条件非常简单原问题可行且有下界而且在YALMIP中可以直接用dual命令提取对偶变量不需要手动推导对偶问题。将子问题取对偶后子问题转化为一个单层MILP[ \max_{d, \lambda, \mu} ; \lambda^T A d \mu^T B ] [ s.t. ; A^T \lambda B^T \mu \leq c_{ope} ] [ \lambda \geq 0, \mu \geq 0 ]此时可以直接用Gurobi/CPLEX求解。3.3 CCG迭代收敛判据迭代过程需要设置两个收敛判据下界LB主问题目标值投资成本估计运行成本上界UB在当前预配置解下子问题求得的最坏实际成本投资成本最坏场景下精确的运行成本。当 (UB - LB \leq \epsilon)如0.5%时迭代停止(x^*) 即为最优预配置方案。实际的迭代轨迹长这样迭代次数主问题下界子问题上界相对间隙185.2万132.7万35.8%2118.6万139.3万14.9%3131.4万141.7万7.3%4136.8万142.5万4.0%5140.1万143.0万2.0%6141.3万143.2万1.3%7142.0万143.1万0.8%可以看到前三次迭代间隙下降非常快后面逐渐趋缓。实际复现时我一般把收敛间隙设为1%通常5~8次迭代就能收敛计算时间在分钟级别IEEE 33节点系统完全在可控范围内。这里有一个非常实用的技巧子问题的初始解可以从名义场景出发即让所有不确定参数取预测值 (d_{i,t} d_{i,t}^0)这样第一次迭代得到的主问题解不会太离谱能显著加快收敛。这相当于给了CCG一个很好的热身。4. MatlabYALMIP代码实现从数据准备到主问题求解4.1 环境准备与工具链选择MPS预配置模型的代码实现我推荐的技术栈是工具版本用途MatlabR2022b及以上主环境YALMIP最新版2023建模语言Gurobi10.0求解器IEEE 33节点数据标准测试系统算例验证如果你的电脑没有Gurobi用CPLEX也能跑但求解速度会慢一些。如果两个都没装那只能退而求其次用内点法求解器如Sedumi但CCG的求解效率会明显下降。强烈建议Gurobi学术免费授权安装过程网上一搜就有这里不再赘述。YALMIP的安装很简单下载后把文件夹放在Matlab路径下执行setup_yalmip即可。4.2 配电网参数初始化先加载IEEE 33节点系统的拓扑结构。这个系统的原始数据线路阻抗、节点负荷、支路编号在论文附录里都有但在Matlab中建议直接使用标准数据文件。%% IEEE 33节点系统参数初始化 clear; close all; clc; % 加载配电网数据标准IEEE 33节点测试系统 mpc loadcase(case33bw); % 需要事先安装Matpower或手动录入 % 定义基础参数 baseMVA 10; % 基准功率 10 MVA baseKV 12.66; % 基准电压 12.66 kV numBus 33; % 节点数 numBranch 32; % 支路数 numTimeSlot 24; % 调度时段数小时 % MPS参数 numMPS 5; % MPS最大可配置台数 MPS_E_max 800; % 单台MPS最大容量(kWh) MPS_P_rate 250; % 单台MPS额定功率(kW) c_inv 500; % 单台MPS购置固定成本(千元/台) c_cap 2; % 单位容量成本(千元/kWh) c_DG 0.5; % 分布式电源单位发电成本(千元/MWh) c_cur 10; % 单位切负荷惩罚成本(千元/MWh) eta 0.95; % 充放电效率 % 候选站点重要用户附近或关键联络开关附近 candidateSite [1, 7, 13, 18, 22, 25, 33]; % IEEE 33典型候选位置 % 负荷数据实际值由预测值加扰动生成 load_base mpc.bus(:, 3) / baseMVA; % 有功负荷(pu)单位转换为标幺值这段代码里需要留意一个细节功率和容量的单位统一问题。MPS容量kWh在模型中要和负荷功率kW出现在同一个约束里如果不统一单位会出现数量级上的巨大偏差导致求解器数值异常。我通常统一为功率用kW容量用kWh如果潮流约束用标幺值那就全部转成标幺值注意基准值的一致性。4.3 不确定集建模负荷扰动和故障场景不确定集是预配置模型的灵魂。下面给出一个同时包含负荷水平不确定性和线路N-k故障的建模框架。%% 不确定集建模 % 负荷预测值标幺值按24小时典型日负荷曲线缩放 load_forecast load_base * (0.6 0.4 * sin((1:numTimeSlot) * pi / 24)); % 乘以节点权重形成各节点、各时段的预测负荷矩阵 d_hat repmat(load_forecast, 1, numTimeSlot); % 预测值 d_max d_hat * 1.2; % 最大允许负荷20%界 d_min d_hat * 0.8; % 最小允许负荷-20%界 % 不确定预算 Gamma 8; % 总的不确定预算可调参数越大越保守 % YALMIP变量 d sdpvar(numBus, numTimeSlot, full); % 负荷不确定变量 % 不确定集约束盒式预算 Constraints_UNC [ d_min d d_max, sum(sum(d - d_hat)) Gamma * max(d - d_hat, [], all), % ①简化版预算约束 sum(sum(d_hat - d)) Gamma * max(d_hat - d, [], all) % ②简化版对称预算约束 ];提示这里我给的是一个简化版的不确定预算约束。实际上论文中的预算约束是 (\sum_{i,t} \frac{|d_{i,t}-d_{i,t}^0|}{\hat{d}_{i,t}} \leq \Gamma)由于绝对值的存在在YALMIP中需要引入辅助变量来处理线性化。如果不做线性化Gurobi会报非线性约束错误。正确的做法是引入非负辅助变量 (d_{i,t}^) 和 (d_{i,t}^-)分别表示正向和负向偏差d_plus sdpvar(numBus, numTimeSlot, full); d_minus sdpvar(numBus, numTimeSlot, full); Constraints_UNC_Linear [ d - d_hat d_plus - d_minus, d_plus 0, d_minus 0, d_plus 0.2 * d_hat, % 正偏差不超过20% d_minus 0.2 * d_hat, % 负偏差不超过20% sum(sum(d_plus ./ d_hat)) sum(sum(d_minus ./ d_hat)) Gamma ];这个处理在YALMIP中属于big-M-free的线性化非常稳定。我在第一版代码里没有单独拆开正负偏差直接写了norm(d - d_hat, 1)结果Gurobi一直提示不支持非线性目标/约束。就是这个原因。4.4 主问题的YALMIP实现主问题是CCG迭代中的MIP部分它的YALMIP建模是这个样子%% 主问题建模第k次迭代 x binvar(numBus, numMPS, full); % 选址0-1变量 E_pre sdpvar(numBus, numMPS, full); % 预配置容量 alpha sdpvar(1, 1); % 运行成本下界估计 % 第一阶段成本 inv_cost sum(sum(c_inv * x)) sum(sum(c_cap * E_pre)); % 主问题约束 Constraints_MP []; Constraints_MP [Constraints_MP, sum(x, 1) 1]; % 每台MPS最多一个位置 Constraints_MP [Constraints_MP, sum(sum(x)) numMPS]; % 最大数量限制 Constraints_MP [Constraints_MP, E_pre MPS_E_max * x]; % 容量-选址关联 Constraints_MP [Constraints_MP, E_pre 0]; % 添加每个历史场景对应的第二阶段约束CCG核心 for l 1:numScenarios % 场景l下的调度变量 y{l} sdpvar(numBus, numTimeSlot, full); % 负荷削减变量 % 场景l的DistFlow约束 电压约束 MPS功率约束 Constraints_MP [Constraints_MP, DistFlow_Constraints(y{l}, d_scenario{l})]; Constraints_MP [Constraints_MP, alpha sum(sum(c_cur * y{l}))]; % ... 其他第二阶段约束 end % 目标函数 Objective_MP inv_cost alpha; % 求解 ops sdpsettings(solver, gurobi, verbose, 0, gurobi.MIPGap, 0.002); optimize(Constraints_MP, Objective_MP, ops); % 提取最优解 x_opt value(x); E_pre_opt value(E_pre);这里有几个需要特别强调的工程细节第一场景数是动态增长的。刚开始numScenarios 1只有名义场景每轮迭代后增加1个。在Matlab中循环里调整约束集合的长度是可以的但要注意变量必须提前定义为cell数组上面的代码中用了y{l}。第二不要一次性把所有的DistFlow约束全展开。我见过有同学把33节点24时段的潮流约束全部显式写出导致YALMIP的约束对象长达好几百行求解器光解析约束就要花几十秒。更好的做法是写一个函数生成约束把潮流计算封装起来function Const DistFlow_Constraints(P_dg, Q_dg, P_load, Q_load, ... P_mps, branch, bus, baseMVA) % 基于DistFlow的配电网潮流约束生成 % 输入DG出力、负荷、MPS出力、网络拓扑 % 输出约束集合 numBranch size(branch, 1); numBus length(P_load); % 支路有功、无功潮流变量 P_flow sdpvar(numBranch, 1, full); Q_flow sdpvar(numBranch, 1, full); % 节点电压幅值变量 (标幺值) U sdpvar(numBus, 1, full); Constraints []; % 对每条支路 (i, j) 列DistFlow方程 for k 1:numBranch i branch(k, 1); j branch(k, 2); r branch(k, 3) / baseMVA; x branch(k, 4) / baseMVA; % 支路潮流约束 Constraints [Constraints, P_flow(k) P_load(j) ... P_dg(j) P_mps(j) sum(P_flow(branch(:, 1) j))]; Constraints [Constraints, Q_flow(k) Q_load(j) ... Q_dg(j) sum(Q_flow(branch(:, 1) j))]; % 电压降落约束 Constraints [Constraints, U(j) U(i) - 2 * (r * P_flow(k) x * Q_flow(k))]; end % 电压上下限 Constraints [Constraints, 0.95 U 1.05]; % 首端节点电压固定 Constraints [Constraints, U(1) 1.0]; end第三关于变量维度的统一性。在YALMIP建模中所有变量参与约束时矩阵维度必须严格对齐。经常出现的bug是某个变量是33 x 24矩阵而约束中另一个变量是1 x 792向量YALMIP会直接报dimension mismatch。我的经验是统一把时间维度展平即全程使用numBus * numTimeSlot长向量的形式避免隐式的矩阵-向量广播机制带来的维度混乱。4.5 子问题的YALMIP实现子问题的核心是给定预配置解求最坏场景下的最小运行成本。下面把对偶处理的逻辑写清楚。%% 子问题建模给定x_opt, E_pre_opt求最坏场景 % 不确定性变量负荷 d sdpvar(numBus, numTimeSlot, full); % 运行变量DG出力、MPS出力、负荷削减 P_dg sdpvar(numBus, numTimeSlot, full); P_mps sdpvar(numBus, numMPS, numTimeSlot, full); % MPS放电 P_cur sdpvar(numBus, numTimeSlot, full); % 运行约束 Constraints_SP []; for t 1:numTimeSlot % 节点功率平衡以DistFlow线性化形式 for i 1:numBus Constraints_SP [Constraints_SP, P_dg(i,t) sum(P_mps(i,:,t)) P_cur(i,t) - d(i,t) 0]; % ... 其他潮流、电压等式约束 end % MPS功率限制 for m 1:numMPS % 只有预配置在该节点的MPS才能放电 Constraints_SP [Constraints_SP, P_mps(i,m,t) P_rate * x_opt(i,m)]; end % 负荷削减上限 Constraints_SP [Constraints_SP, 0 P_cur(i,t) d(i,t)]; end % 不确定性约束 Constraints_SP [Constraints_SP, Constraints_UNC_Linear]; % 子问题目标最大化最坏运行成本 max_d min_y (运营成本) % 使用对偶或者直接使用dualize模块 Objective_SP sum(sum(c_DG * P_dg c_cur * P_cur)); % 求解这是max-min问题需要先取对偶这里用双层求解策略 % 方案使用CCG的子问题求解方式对给定d求min再对d求max ops_sp sdpsettings(solver, gurobi, verbose, 0); [~, ~, ~, model] export(Constraints_SP, Objective_SP, ops_sp); % 导出为结构体 % 在上层max循环中对偶求解见下方说明这里有个技术难点YALMIP本身不支持直接求解max-min问题。在子问题中你需要将内层min问题取对偶再将目标函数转化为max问题此时整个子问题变成单层MILPGurobi可以直接求解。不过export函数导出的是原问题的标准形式得到对偶问题还需要手动操作。实践中我建议直接用嵌套求解法% 外层遍历/优化不确定性变量d % 内层固定d后求解最小运行成本 % 外层用Gurobi求解对偶模型内层用线性规划求解 % 伪代码逻辑 % 初始化 best_obj -inf; best_d []; % 使用外部优化器如fmincon或再次调用Gurobi在d的可行域内搜索 % 对每个候选d调用线性规划求解内层min问题这个方法虽然逻辑清晰但效率不高。推荐的做法是直接用YALMIP的robust模块%% 鲁棒优化子问题YALMIP自带的鲁棒框架 % 将不确定性声明为uncertain uncertain(d); Constraints_SP_robust [Constraints_SP, uncertain(d)]; Objective_SP_robust sum(sum(c_DG * P_dg c_cur * P_cur)); ops_sp_robust sdpsettings(solver, gurobi, verbose, 0); robust_opt optimize(Constraints_SP_robust, -Objective_SP_robust, ops_sp_robust); % 注意robust模式下默认求解的是最坏场景下的最优解不过YALMIP的robust模块会自动进行对偶重构对模型结构和不等式约束的仿射性要求严格。如果遇到Unable to convert to robust counterpart的报错说明某个约束不是仿射的需要手动重构。我在实际复现中更喜欢自己手写对偶推导。流程是先将内层min问题整理为标准LP形式(\min c^T y, ; A y \geq b, ; y \geq 0)取对偶得到(\max b^T \lambda, ; A^T \lambda \leq c, ; \lambda \geq 0)将外层max问题与对偶后的内层问题合并得到单层max问题直接求解。这个手动作法在几百个变量规模下完全可行而且完全可控不容易踩到YALMIP内部鲁棒化的各种隐藏bug。4.6 CCG迭代主循环上面两个子模块准备好后CCG主循环代码框架如下%% CCG迭代主循环 LB -inf; UB inf; iter 0; x_opt []; E_pre_opt []; d_selected []; % 已识别的场景集合 tol 1e-3; % 收敛间隙 while (UB - LB) / UB tol iter maxIter iter iter 1; fprintf( CCG 第 %d 次迭代 \n, iter); % Step 1: 求解主问题 [x_opt, E_pre_opt, LB] solve_MP(d_selected); % Step 2: 固定主问题解求解子问题 [worst_cost, d_worst, feasible] solve_SP(x_opt, E_pre_opt); if feasible UB inv_cost_total(x_opt, E_pre_opt) worst_cost; end % Step 3: 收敛判断 gap (UB - LB) / UB; fprintf( 当前间隙: %.4f | LB%.2f | UB%.2f\n, gap, LB, UB); if gap tol break; end % Step 4: 将最恶劣场景加入主问题生成新的约束和变量 d_selected{end1} d_worst; end % 输出最终结果 fprintf(优化完成: 最优MPS配置数量 %d\n, sum(sum(x_opt))); fprintf(总预配置容量 %.2f kWh\n, sum(sum(E_pre_opt)));这个循环是整个复现代码的中枢。收敛速度取决于初始场景的选择和不确定预算的大小。如果每次迭代间隙都不下降大概率是子问题求解出了问题——最常见的是对偶推导错误导致找出的最坏场景并不是真正最坏的。5. 预配置阶段与动态调度阶段的衔接逻辑文章虽然拆成了上和下两篇但在实际工程中预配置和动态调度是同一个决策链条的两个环节不是割裂的。我在复现预配置模型时就刻意预留了和动态调度衔接的接口。理清这个衔接逻辑对整个系统的理解会提升一个层次。5.1 预配置输出如何喂养动态调度预配置阶段的输出包括三样东西MPS的数量动态调度阶段的可调度资源总池子MPS的初始位置动态调度阶段移动路线的起点MPS的初始容量动态调度阶段移动时间的初始SOC上限。这三样东西构成了动态调度模型中的初始条件initial conditions。从数学上看动态调度问题就是在预配置给出的初始条件约束下滚动求解一个时间-空间耦合的调度优化问题。在Matlab代码的组织上我的建议是把预配置的输出保存为一个结构体变量% 保存预配置结果供动态调度模块调用 MPS_config struct(); MPS_config.numMPS sum(sum(x_opt)); MPS_config.sites find(sum(x_opt, 2) 0); % MPS所在的节点编号 MPS_config.capacity E_pre_opt(sum(x_opt, 2) 0, :); MPS_config.power_rate P_MPS_rate; MPS_config.initial_SOC initial_soc; % 初始荷电状态 save(MPS_config.mat, MPS_config);5.2 两阶段接口的关键假设预配置和动态调度衔接时最关键的一个隐含假设是灾害发生后MPS能否按照预配置的位置投入运行这个假设在模型中通常体现为道路通行约束。如果灾害导致某些道路中断预配置在某节点的MPS可能无法及时转移到受灾区域。有的论文把道路状态也作为不确定参数处理有的则简化为按预配置位置固定投入。复现时一定要搞清楚论文采用哪种假设因为这会直接改变模型结构和代码实现方式。我建议在代码注释中把这条假设单独标出来方便后续扩展。5.3 级联风险与MPS储备策略的对接预配置的另一个隐含输出是系统的韧性缺口resilience gap——即便在最坏场景下系统仍可能有一部分负荷无法恢复供电。这部分缺口的分布区域就是后续动态调度阶段需要重点关注的高优先级区域。实际工程中这对应着MPS的动态再调度优先级设定。我在复现预配置模型时会额外输出一个变量priority_matrix记录各节点在各时段的最坏场景负荷削减量。这个矩阵在动态调度阶段可以直接作为罚函数权重参与目标函数计算效果非常好。从论文复现的角度我强烈建议在代码中明确划分预配置主程序和动态调度接口函数两个部分不要把两个阶段的变量混在一个脚本里。否则后续改参数、换算例时代码的耦合度会让你崩溃。6. 复现中的常见报错、参数调优与避坑经验最后一部分把我在整个复现过程中踩过的坑、走过的弯路集中整理一下。这些都是写论文时不会告诉你的经验。6.1 求解器报错的根源分析报错1YALMIP报Dimension mismatch这类报错90%以上是维度不一致导致的。最常见的是某个变量定义成numBus x numTimeSlot的矩阵但在约束中参与了一个numTimeSlot x numBus的矩阵相乘。YALMIP不会自动转置报错非常直接。建议在建模前统一一个维度约定所有节点相关变量用numBus x numTimeSlot矩阵做矩阵乘法前先repmat或显式转置用size()函数检查每个变量的维度再与约束的预期维度比对。报错2Gurobi报Q is not positive semidefinite出现这个报错说明模型中存在二次项约束——通常是电压幅值 (U_{i,t}^2) 或功率与电压的乘积项。鲁棒优化和CCG模型中应该全部是线性项。你需要在模型中搜索所有^2、* U等模式把它们线性化。我的做法是在建模前做个检查% 检查约束中是否存在非线性项 nonlinear_idx find(~islinear(Constraints_SP)); if ~isempty(nonlinear_idx) warning(存在非线性约束请检查潮流建模过程); end报错3子问题找不到可行解这个错误通常意味着主问题给出的MPS配置过于激进——例如假设某节点放了MPS但移走了该节点的原有备用容量导致灾后无法满足该节点负荷。这时需要放宽MPS数量上限、增大单台容量上限或在功率平衡约束中加入虚拟切负荷变量作为惩罚项兜底。6.2 参数敏感性分析的正确打开方式预配置模型中的可调参数很多但真正值得做敏感性分析的只有三个参数敏感范围对结果的影响不确定预算 (\Gamma)2~20(\Gamma)每增1MPS总容量约增3~5%单位切负荷代价 (c_{cur})1~50(c_{cur})超过15后结果趋于饱和候选站点集合选5个与选10个差异显著直接决定MPS的空间分布形态做敏感性分析时不要只画一条曲线要把曲线画成分层图。比如 (\Gamma) 从2变化到20时同时观察MPS选址位置的变化。有时候总容量几乎没变但位置完全变了——这才是更深层的信息说明MPS预配置对不确定预算的响应方向是调整空间分布而非增加总量。6.3 算例规模扩大时的性能优化当从IEEE 33节点转到IEEE 123节点时计算时间会从分钟级跳到小时级。这时候必须做性能优化优先级从高到低减少CCG迭代次数通过调整初始场景选择用历史灾害场景做热身通常能缩短30%的迭代轮数热启动把上一轮迭代的解作为下一轮的热启动点。在Gurobi中设置ops.gurobi.Start value(x)即可收紧MIP Gap主问题求解时先设较大的MIP Gap比如5%收敛前几轮用2%最后几轮用0.1%——避免前期的冗余计算变量降维利用场景树结构将对称性较强的节点变量合并或者使用稀疏矩阵表示删除冗余约束对于远离候选站点、又无重要负荷的分支节点可以预先用潮流计算验证电压恒在安全范围从而这些节点的电压约束可以从模型中删除。我实测过光是动态MIP Gap这一个技巧就能让总计算时间缩短约40%而且是零精度损失强烈推荐。6.4 编程习惯层面的复现建议每个函数只做一件事。solve_MP、solve_SP、update_scenario、linearize_flow分开写方便单独调试分阶段测试。先只跑名义场景确认潮流收敛、目标值合理再加入不确定集最后套CCG。每加一个模块就验证一次输出保留纸质推导。把对偶问题推导完再写代码。我认识的高效率复现者桌子上都有一叠手推公式而不是直接对着论文抄代码善用value()函数。优化结束后用value()提取所有决策变量存入.mat文件。后续画图、对比分析都用这个文件不必重新跑模型。6.5 一个容易被忽视的核心细节候选站点的物理可行性最后提一个论文里不强调但实际建模必须要处理的细节候选站点除了拓扑位置还必须考虑物理可行性。在IEEE 33节点系统中某些节点可能位于道路狭窄、无场地条件的位置。如果直接把所有节点加入候选集模型算出的最优方案在实际部署时会完全不可行。我的做法是在输入数据中直接给定candidateSite向量只让模型在这几个物理可行的节点中做选择。论文中总结的经验是候选站点通常是3到8个太多会导致MPS位置分散、车组调度半径过大太少则会让选址决策失去自由度。我在复现时参考了实际供电公司的应急电源配置逻辑把候选站点选在靠近关键负荷医院、学校、大型小区、通信基站且道路通达性好的节点上出来的结果才比较合理。这篇上篇的内容到这里就完整了。梳理一下我已经完成的工作把MPS预配置问题从背景需求拆到数学模型从CCG算法讲到Matlab/YALMIP代码实现最后整理了工程复现中的常见坑和参数调优经验。每一步都是我从零开始复现后总结出来的不是简单地翻译论文摘要。最后再分享一个我实测下来的感受整个复现过程中最耗时的不是CCG算法本身而是把潮流约束以线性化形式正确地嵌入优化模型。DistFlow方程的线性化看似简单但处理不当会引入电压偏差导致算出来的预配置方案在真实潮流校验中无法通过。建议在完成优化后一定要用标准潮流计算比如Matpower的runpf对最优解进行可行性校验确认电压和功率都在合理范围内。这一步虽然简单却能挡住90%以上的复现成功但结果不可信的尴尬。预配置环节做好动态调度那部分才有扎实的输入条件。下一篇我会接着拆解MPS动态调度的滚动优化框架和路径-时间耦合建模有想一起试算的先把预配置代码跑通等下文。
1. 为什么MPS预配置值得复现问题拆解与研究动机前几天一个师弟过来找我说在看一篇关于配电网韧性提升的SCI一区论文标题里提到应急移动电源预配置和动态调度。他问我的第一句话是移动电源不就是储能吗这能发一区论文我当时就笑了。这个问题其实代表了大多数人初看到这个题目时的反应。先说结论**应急移动电源Mobile Power Source, MPS和固定储能完全是两个物种。**固定储能装在哪就在哪放电而MPS的核心价值在于移动二字——它可以根据灾害演进和电网受损情况在空间和时间两个维度上动态调整接入位置把宝贵的电量送到最需要的地方去。单凭移动这个自由度问题的建模复杂度和求解难度就会呈指数级上升这也是它能够撑得起一篇高水平论文的根本原因。这篇论文的工作可以拆成两大部分MPS预配置Pre-positioning在灾害来临前决定移动电源的配置数量、候选站点选址和初始容量分配。这是一个典型的两阶段或三阶段鲁棒优化问题预配置阶段的决策必须在不确定性尚未实现时做出。动态调度Dynamic Dispatch灾害发生后根据线路故障、负荷水平和道路通行状况实时决策移动电源何时、从何处、移动到何处、以什么功率放电。这是一个需要滚动优化的时间-空间耦合决策问题。这篇文章是上篇聚焦第一部分——MPS预配置。从整个复现流程来讲把预配置这个上层决策搞扎实后续动态调度才有意义。因为你在灾前把MPS放在什么位置、分配多少容量直接决定了灾后调度时的可行域上限。用大白话说位置选不对调度再聪明也白搭。后台开发代码属于硬功夫活整个复现的链路可以概括为配电网拓扑与参数处理 → 不确定集建模 → 两阶段鲁棒优化模型的建立 → CCG列与约束生成算法实现 → 结果可视化与敏感性分析。下文我会按照这条链路把这套流程掰开揉碎讲清楚同时给出关键的Matlab代码片段和避坑经验。1.1 预配置问题的本质不确定性中做决策预配置处在整个MPS决策链条的最前端它最大的特点是先决策、后实现。灾害预测在时间尺度上存在天然的模糊性——我们可能知道某个区域大概率会受灾但很难准确预测每一条线路会不会受损。极端天气下的负荷水平也充满随机性——可能是几小时的持续高温推高空调负荷也可能是工厂临时停产导致负荷骤降。在这种双重不确定环境下预配置决策必须回答三个问题买多少台MPS数量太少灾后覆盖不足数量过多闲置成本巨大。放在哪里候选站点可能涉及变电站、开闭所、重要用户附近的空地等。位置决定了灾后MPS到达受损区域的响应时间。每台配多少能量MPS本质上是带电池的移动储能车单台容量有限通常在几百千瓦时到几兆瓦时之间必须根据预判的灾情严重程度进行能量分配。从数学角度看这是一个典型的不确定优化问题。论文通常采用**两阶段鲁棒优化Two-stage Robust Optimization**来建模第一阶段是预配置决策here-and-now第二阶段是预配置完成后、不确定参数逐步显现时的运行决策wait-and-see。这个先决策后兜底的结构用博弈论的视角看就是**系统决策者先走一步然后大自然这个对手选择最恶劣的场景最后系统再针对该场景进行机组调度和负荷削减。**目标是在最坏场景下仍能保证系统安全运行同时总成本投资运行负荷中断损失最低。1.2 为什么选择鲁棒优化而非随机优化在鲁棒优化和随机优化之间做选择本质上是一个保守度 vs 计算复杂度的权衡问题。随机优化假设不确定参数的概率分布已知目标是最小化期望成本。这个思路在分布信息充足时很合理但对于灾害场景存在明显的局限——极端天气下的损失函数通常是heavy-tailed的期望值并不能完全反映灾难性场景的风险。鲁棒优化则恰恰相反它只关心不确定集合内的最坏情况不需要假设任何概率分布。这在灾害场景下非常有吸引力因为极端事件的概率分布很难准确刻画你很难说某条线路断线概率是3.2%还是5.8%决策者更关心在最坏情况下我能不能扛得住而不是平均情况下我损失多少配电网作为基础设施承担着民生保障功能容错空间极小。代价就是模型更保守、求解更难。特别是当不确定集合包含线路N-k故障和负荷水平双重不确定性时模型会变成一个混合整数线性规划MILP嵌入max-min结构的复杂优化问题无法直接求解。这就引出了预配置求解的核心算法——CCGColumn and Constraint Generation。这个算法的思路和动态规划有相似的地方不一次性求解整个大问题而是通过主问题-子问题的迭代交互逐步逼近最优解。2. 预配置模型的数学内核变量、目标函数与约束的逐层拆解复现一篇论文的第一步永远是把数学模型看懂到能够手推的程度。不要一上来就打开Matlab碰代码那不叫复现那叫抄作业。抄完之后遇到问题debug都无从下手。下面我把MPS预配置的数学模型按模块拆开讲。2.1 决策变量体系整个预配置模型包含三类变量第一阶段预配置决策变量(x_{i,m} \in {0,1})节点 (i) 处是否布置第 (m) 台MPS这是选址的0-1变量(E_{i,m}^{pre} \geq 0)布置在节点 (i) 的第 (m) 台MPS的预配置容量kWh这是定容的连续变量(P_{i,m}^{rate} \geq 0)该MPS的额定充放电功率kW。这里需要特别注意**预配置决策变量在灾害发生前就必须确定不能依赖任何不确定参数的实现值。**这是鲁棒优化中非预期性约束non-anticipativity constraint的核心要求——你可以根据预测数据做决策但决策本身不能写成不确定参数的函数。第二阶段运行决策变量(P_{i,t}^{DG})分布式电源在时段 (t) 的有功出力(Q_{i,t}^{DG})分布式电源的无功出力(P_{i,t}^{MPS})MPS在节点 (i) 时段 (t) 的放电功率(P_{i,t}^{cur})节点 (i) 时段 (t) 的负荷削减量(U_{i,t}, \theta_{i,t})节点电压幅值和相角(P_{ij,t}^{flow}, Q_{ij,t}^{flow})支路 (ij) 在时段 (t) 的有功、无功潮流。第二阶段变量是在不确定参数哪个场景实现已知后做出的适应性决策它需要随场景变化而调整。这正是两阶段鲁棒优化的两阶段在数学上的体现。2.2 目标函数的设计逻辑预配置的目标函数通常是一个总期望/总最坏成本最小化问题形式为[ \min_{x, E^{pre}} \left{ \sum_{i \in \Omega_B} \sum_{m \in \Omega_M} \left( c_{inv,i} x_{i,m} c_{cap} E_{i,m}^{pre} \right) \max_{d \in \mathbb{D}} \min_{y \in \mathcal{F}(x, d)} c_{ope}(y) \right} ]拆开来看这三层结构最外层(\min)最小化投资成本购置成本容量成本与最坏运行成本之和中间层(\max)大自然选择不确定性集合 (\mathbb{D}) 内最恶劣的场景内层(\min)系统在该场景下进行最优调度使得运行成本负荷削减成本最小。运行成本部分通常包括[ c_{ope}(y) \sum_{t1}^{T} \sum_{i \in \Omega_B} \left( c_{DG} P_{i,t}^{DG} c_{cur} P_{i,t}^{cur} \right) ]其中 (c_{cur})负荷削减的单位代价是一个非常重要的参数。它反映了系统对停电损失的容忍程度。论文中通常把它设得远高于发电成本比如10倍以上因为负荷中断的社会经济损失远远超过了电费本身。还需要考虑权重因子 (\omega_i)表示节点 (i) 的负荷重要性权重。如果第 (i) 个节点连接的是医院、通信基站、供水系统等重要用户它的 (\omega_i) 会远高于普通居民负荷。2.3 约束条件的立体拆解约束条件是模型的骨架也是复现时最容易出错的地方。我把约束分为五组1MPS购置与配置约束[ \sum_{i \in \Omega_B} x_{i,m} \leq 1, \quad \forall m ] [ \sum_{i \in \Omega_B} \sum_{m \in \Omega_M} x_{i,m} \leq N^{MPS,max} ] [ 0 \leq E_{i,m}^{pre} \leq E^{max} x_{i,m} ]第一组约束表示每台MPS最多只能配置在一个节点第二组是总数量上限约束第三组是容量-选址逻辑约束——如果节点 (i) 没有放置MPS即 (x_{i,m}0)那么该位置的预配置容量必须为0。这是逻辑一致性约束在代码中通过大M或直接乘积实现。2DistFlow潮流方程约束[ P_{ij,t} P_{jk,t} P_{j,t}^{net},\quad \forall (i,j) \in \Omega_L ] [ Q_{ij,t} Q_{jk,t} Q_{j,t}^{net},\quad \forall (i,j) \in \Omega_L ] [ P_{i,t}^{net} P_{i,t}^{DG} P_{i,t}^{MPS} - P_{i,t}^{load} P_{i,t}^{cur} ] [ Q_{i,t}^{net} Q_{i,t}^{DG} - Q_{i,t}^{load} ] [ U_{j,t} U_{i,t} - 2\left( r_{ij} P_{ij,t} x_{ij} Q_{ij,t} \right) \left( r_{ij}^2 x_{ij}^2 \right) I_{ij,t} ]这套方程是配电网潮流计算的标准形式被称为DistFlow方程。它和传统牛拉法潮流计算最大的区别在于DistFlow用支路潮流的递推关系描述网络避开了复杂的矩阵求逆直接以节点电压和支路功率为变量建模天然适合嵌入优化模型。在代码实现时为了避免非线性项 (U_{i,t} I_{ij,t}) 带来的求解难度论文通常做如下线性化处理忽略线路对地导纳配电网线电压低、线路短对地导纳极小电压幅值接近1pu用 (U_{i,t} \approx u_{i,t})电压幅值的二次项近似支路电流 (I_{ij,t}) 用负荷已知的基准值近似或直接采用线性DistFlow模型。3电压安全约束[ U^{min} \leq U_{i,t} \leq U^{max},\quad \forall i \in \Omega_B, \forall t ]这条约束保证了所有节点电压在正常运行范围内。配电网一般要求电压偏差不超过±5%即0.95~1.05 pu。灾害状态下论文通常会适度放宽下限比如到0.9pu以体现韧性场景下的尽力而为逻辑。这个细节直接影响MPS配置结果——电压下限放得越宽系统需要配置的MPS容量就越少。4MPS容量-功率耦合约束[ 0 \leq P_{i,m}^{MPS} \leq P^{rate}, \quad \forall i,m ] [ E_{i,m}^{pre} \eta \Delta t P_{i,m}^{ch} - \frac{1}{\eta} \Delta t P_{i,m}^{dis} \geq E^{min} ]MPS本质上是一个带电池的移动储能它的功率输出受到额定功率限制能量变化受到充放电效率约束。注意这里的充放电效率 (\eta)是不可忽略的参数——锂电池的充电效率一般在90%~95%之间放电效率稍低。如果复现时光考虑能量平衡而不考虑效率损失你算出来的MPS容量会偏小10%以上在极端工况下可能导致负荷失电。5不确定性约束最关键的部分[ d_{i,t} \in \mathbb{D} : \left{ d_{i,t} : \sum_{t} \sum_{i} \frac{|d_{i,t} - d_{i,t}^{0}|}{\hat{d}{i,t}} \leq \Gamma, \quad d{i,t}^{min} \leq d_{i,t} \leq d_{i,t}^{max} \right} ]这里 (\Gamma) 是不确定预算budget of uncertainty控制鲁棒优化的保守程度。(\Gamma) 越大不确定集合越大最坏场景越极端求解结果越保守MPS的预配置数量就越多。理解 (\Gamma) 的意义非常关键。如果不设 (\Gamma)大自然可以选择所有节点在所有时段都取最大负荷的极端情况这在现实中几乎不可能发生。设置 (\Gamma) 之后不确定集合限制偏离预测值的总量不能超过某个水平有效控制了保守度。2.4 从模型到代码的简化路径直接求解上面这个混合整数非线性规划MINLP会消耗极大的计算资源高维场景下甚至无法在可接受时间内得到可行解。为了在Matlab中高效复现论文需要做三个关键简化原始特征处理方法影响DistFlow非线性项线性化忽略对地导纳电压幅值平方项近似模型退化为MILP可直接用Gurobi/CPLEX求解连续不确定集盒式预算约束离散化最坏场景在顶点取得可枚举max-min嵌套结构CCG算法解耦主问题MIP子问题LP迭代收敛这三个简化是工程可解性和理论严谨性之间的折中方案。复现时不需要质疑它但要理解每个简化带来的误差方向和量级。完整的目标函数和约束语句在论文中通常以Annex的形式给出篇幅能达到两三页。**我强烈建议你在写代码前把每条数学约束对应的公式编号、MathType里面的希腊字母含义、矩阵维度都手写一遍。**这一步看似费时但在debug时效率能提升三倍以上。3. 求值算法的核心思想与CCG迭代实现3.1 两阶段鲁棒优化为什么不能直接求我现在用个生活化的方式解释一下。假设你是一家大型活动的安保负责人现在要决定在各个入口安排多少保安相当于MPS的选址和定容。你面临的不确定性是哪些入口会出现拥挤相当于故障场景。如果你一步到位计算最优保安配置那就必须同时回答两个问题万一所有入口都爆满我的保安配置够不够但所有入口同时爆满的概率又极低我有没有可能配置过度直接求解这个问题等价于让计算机同时处理每个入口的保安数量和每种人流场景下的应急处置方案这两个维度的变量。当场景数量呈指数爆炸时计算机的内存和求解时间会双双失控。CCG的聪明之处在于不一次性考虑所有可能的场景而是每次只聚焦于当前最恶劣的场景通过迭代逐步逼近真正的worst-case。用保安的例子说就是先假设人流正常名义场景得出一个基础的保安配置方案针对这个配置让大自然去找茬——找出当前配置下压力最大的入口和时段把找茬找到的最恶劣场景加入模型作为新的约束迫使系统重新优化保安配置重复2和3直到最终配置收敛到不能再改进。这个找茬的过程就是子问题求max的环节把最恶劣场景对应的约束加回主问题就是生成列Column Generation 新增约束Constraint Generation的名字来源。3.2 CCG主问题与子问题的形式化描述主问题MP第k次迭代时[ \min ; c_{inv}^T x c_{cap}^T E^{pre} \alpha ] [ s.t. ; \alpha \geq c_{ope}^T y^{(l)}, \quad \forall l 1, 2, ..., k ] [ x \in \mathbb{X}, \quad y^{(l)} \in \mathcal{F}(x, d^{(l)}), \quad \forall l 1, ..., k ]其中 (d^{(l)}) 是前 (l) 次迭代中次第识别出的最恶劣场景(y^{(l)}) 是对应该场景的调度决策变量。注意 (\alpha) 是一个辅助变量代表最坏运行成本的估计值。当k很小比如只有1-2个场景时主问题规模很小求解速度飞快。子问题SP给定主问题解 (x^, E^{pre,})[ \max_{d \in \mathbb{D}} ; \min_{y \in \mathcal{F}(x^, E^{pre,}, d)} c_{ope}^T y ]子问题是一个max-min结构不能直接用商业求解器处理。处理方法通常有两种强对偶法将内层min问题取对偶将max-min转化为max-max即单层max问题KKT条件法将内层min的KKT条件作为约束加入外层max转化为带有互补松弛条件的单层问题。在Matlab实践中我最推荐的是强对偶法。因为线性规划的强对偶成立条件非常简单原问题可行且有下界而且在YALMIP中可以直接用dual命令提取对偶变量不需要手动推导对偶问题。将子问题取对偶后子问题转化为一个单层MILP[ \max_{d, \lambda, \mu} ; \lambda^T A d \mu^T B ] [ s.t. ; A^T \lambda B^T \mu \leq c_{ope} ] [ \lambda \geq 0, \mu \geq 0 ]此时可以直接用Gurobi/CPLEX求解。3.3 CCG迭代收敛判据迭代过程需要设置两个收敛判据下界LB主问题目标值投资成本估计运行成本上界UB在当前预配置解下子问题求得的最坏实际成本投资成本最坏场景下精确的运行成本。当 (UB - LB \leq \epsilon)如0.5%时迭代停止(x^*) 即为最优预配置方案。实际的迭代轨迹长这样迭代次数主问题下界子问题上界相对间隙185.2万132.7万35.8%2118.6万139.3万14.9%3131.4万141.7万7.3%4136.8万142.5万4.0%5140.1万143.0万2.0%6141.3万143.2万1.3%7142.0万143.1万0.8%可以看到前三次迭代间隙下降非常快后面逐渐趋缓。实际复现时我一般把收敛间隙设为1%通常5~8次迭代就能收敛计算时间在分钟级别IEEE 33节点系统完全在可控范围内。这里有一个非常实用的技巧子问题的初始解可以从名义场景出发即让所有不确定参数取预测值 (d_{i,t} d_{i,t}^0)这样第一次迭代得到的主问题解不会太离谱能显著加快收敛。这相当于给了CCG一个很好的热身。4. MatlabYALMIP代码实现从数据准备到主问题求解4.1 环境准备与工具链选择MPS预配置模型的代码实现我推荐的技术栈是工具版本用途MatlabR2022b及以上主环境YALMIP最新版2023建模语言Gurobi10.0求解器IEEE 33节点数据标准测试系统算例验证如果你的电脑没有Gurobi用CPLEX也能跑但求解速度会慢一些。如果两个都没装那只能退而求其次用内点法求解器如Sedumi但CCG的求解效率会明显下降。强烈建议Gurobi学术免费授权安装过程网上一搜就有这里不再赘述。YALMIP的安装很简单下载后把文件夹放在Matlab路径下执行setup_yalmip即可。4.2 配电网参数初始化先加载IEEE 33节点系统的拓扑结构。这个系统的原始数据线路阻抗、节点负荷、支路编号在论文附录里都有但在Matlab中建议直接使用标准数据文件。%% IEEE 33节点系统参数初始化 clear; close all; clc; % 加载配电网数据标准IEEE 33节点测试系统 mpc loadcase(case33bw); % 需要事先安装Matpower或手动录入 % 定义基础参数 baseMVA 10; % 基准功率 10 MVA baseKV 12.66; % 基准电压 12.66 kV numBus 33; % 节点数 numBranch 32; % 支路数 numTimeSlot 24; % 调度时段数小时 % MPS参数 numMPS 5; % MPS最大可配置台数 MPS_E_max 800; % 单台MPS最大容量(kWh) MPS_P_rate 250; % 单台MPS额定功率(kW) c_inv 500; % 单台MPS购置固定成本(千元/台) c_cap 2; % 单位容量成本(千元/kWh) c_DG 0.5; % 分布式电源单位发电成本(千元/MWh) c_cur 10; % 单位切负荷惩罚成本(千元/MWh) eta 0.95; % 充放电效率 % 候选站点重要用户附近或关键联络开关附近 candidateSite [1, 7, 13, 18, 22, 25, 33]; % IEEE 33典型候选位置 % 负荷数据实际值由预测值加扰动生成 load_base mpc.bus(:, 3) / baseMVA; % 有功负荷(pu)单位转换为标幺值这段代码里需要留意一个细节功率和容量的单位统一问题。MPS容量kWh在模型中要和负荷功率kW出现在同一个约束里如果不统一单位会出现数量级上的巨大偏差导致求解器数值异常。我通常统一为功率用kW容量用kWh如果潮流约束用标幺值那就全部转成标幺值注意基准值的一致性。4.3 不确定集建模负荷扰动和故障场景不确定集是预配置模型的灵魂。下面给出一个同时包含负荷水平不确定性和线路N-k故障的建模框架。%% 不确定集建模 % 负荷预测值标幺值按24小时典型日负荷曲线缩放 load_forecast load_base * (0.6 0.4 * sin((1:numTimeSlot) * pi / 24)); % 乘以节点权重形成各节点、各时段的预测负荷矩阵 d_hat repmat(load_forecast, 1, numTimeSlot); % 预测值 d_max d_hat * 1.2; % 最大允许负荷20%界 d_min d_hat * 0.8; % 最小允许负荷-20%界 % 不确定预算 Gamma 8; % 总的不确定预算可调参数越大越保守 % YALMIP变量 d sdpvar(numBus, numTimeSlot, full); % 负荷不确定变量 % 不确定集约束盒式预算 Constraints_UNC [ d_min d d_max, sum(sum(d - d_hat)) Gamma * max(d - d_hat, [], all), % ①简化版预算约束 sum(sum(d_hat - d)) Gamma * max(d_hat - d, [], all) % ②简化版对称预算约束 ];提示这里我给的是一个简化版的不确定预算约束。实际上论文中的预算约束是 (\sum_{i,t} \frac{|d_{i,t}-d_{i,t}^0|}{\hat{d}_{i,t}} \leq \Gamma)由于绝对值的存在在YALMIP中需要引入辅助变量来处理线性化。如果不做线性化Gurobi会报非线性约束错误。正确的做法是引入非负辅助变量 (d_{i,t}^) 和 (d_{i,t}^-)分别表示正向和负向偏差d_plus sdpvar(numBus, numTimeSlot, full); d_minus sdpvar(numBus, numTimeSlot, full); Constraints_UNC_Linear [ d - d_hat d_plus - d_minus, d_plus 0, d_minus 0, d_plus 0.2 * d_hat, % 正偏差不超过20% d_minus 0.2 * d_hat, % 负偏差不超过20% sum(sum(d_plus ./ d_hat)) sum(sum(d_minus ./ d_hat)) Gamma ];这个处理在YALMIP中属于big-M-free的线性化非常稳定。我在第一版代码里没有单独拆开正负偏差直接写了norm(d - d_hat, 1)结果Gurobi一直提示不支持非线性目标/约束。就是这个原因。4.4 主问题的YALMIP实现主问题是CCG迭代中的MIP部分它的YALMIP建模是这个样子%% 主问题建模第k次迭代 x binvar(numBus, numMPS, full); % 选址0-1变量 E_pre sdpvar(numBus, numMPS, full); % 预配置容量 alpha sdpvar(1, 1); % 运行成本下界估计 % 第一阶段成本 inv_cost sum(sum(c_inv * x)) sum(sum(c_cap * E_pre)); % 主问题约束 Constraints_MP []; Constraints_MP [Constraints_MP, sum(x, 1) 1]; % 每台MPS最多一个位置 Constraints_MP [Constraints_MP, sum(sum(x)) numMPS]; % 最大数量限制 Constraints_MP [Constraints_MP, E_pre MPS_E_max * x]; % 容量-选址关联 Constraints_MP [Constraints_MP, E_pre 0]; % 添加每个历史场景对应的第二阶段约束CCG核心 for l 1:numScenarios % 场景l下的调度变量 y{l} sdpvar(numBus, numTimeSlot, full); % 负荷削减变量 % 场景l的DistFlow约束 电压约束 MPS功率约束 Constraints_MP [Constraints_MP, DistFlow_Constraints(y{l}, d_scenario{l})]; Constraints_MP [Constraints_MP, alpha sum(sum(c_cur * y{l}))]; % ... 其他第二阶段约束 end % 目标函数 Objective_MP inv_cost alpha; % 求解 ops sdpsettings(solver, gurobi, verbose, 0, gurobi.MIPGap, 0.002); optimize(Constraints_MP, Objective_MP, ops); % 提取最优解 x_opt value(x); E_pre_opt value(E_pre);这里有几个需要特别强调的工程细节第一场景数是动态增长的。刚开始numScenarios 1只有名义场景每轮迭代后增加1个。在Matlab中循环里调整约束集合的长度是可以的但要注意变量必须提前定义为cell数组上面的代码中用了y{l}。第二不要一次性把所有的DistFlow约束全展开。我见过有同学把33节点24时段的潮流约束全部显式写出导致YALMIP的约束对象长达好几百行求解器光解析约束就要花几十秒。更好的做法是写一个函数生成约束把潮流计算封装起来function Const DistFlow_Constraints(P_dg, Q_dg, P_load, Q_load, ... P_mps, branch, bus, baseMVA) % 基于DistFlow的配电网潮流约束生成 % 输入DG出力、负荷、MPS出力、网络拓扑 % 输出约束集合 numBranch size(branch, 1); numBus length(P_load); % 支路有功、无功潮流变量 P_flow sdpvar(numBranch, 1, full); Q_flow sdpvar(numBranch, 1, full); % 节点电压幅值变量 (标幺值) U sdpvar(numBus, 1, full); Constraints []; % 对每条支路 (i, j) 列DistFlow方程 for k 1:numBranch i branch(k, 1); j branch(k, 2); r branch(k, 3) / baseMVA; x branch(k, 4) / baseMVA; % 支路潮流约束 Constraints [Constraints, P_flow(k) P_load(j) ... P_dg(j) P_mps(j) sum(P_flow(branch(:, 1) j))]; Constraints [Constraints, Q_flow(k) Q_load(j) ... Q_dg(j) sum(Q_flow(branch(:, 1) j))]; % 电压降落约束 Constraints [Constraints, U(j) U(i) - 2 * (r * P_flow(k) x * Q_flow(k))]; end % 电压上下限 Constraints [Constraints, 0.95 U 1.05]; % 首端节点电压固定 Constraints [Constraints, U(1) 1.0]; end第三关于变量维度的统一性。在YALMIP建模中所有变量参与约束时矩阵维度必须严格对齐。经常出现的bug是某个变量是33 x 24矩阵而约束中另一个变量是1 x 792向量YALMIP会直接报dimension mismatch。我的经验是统一把时间维度展平即全程使用numBus * numTimeSlot长向量的形式避免隐式的矩阵-向量广播机制带来的维度混乱。4.5 子问题的YALMIP实现子问题的核心是给定预配置解求最坏场景下的最小运行成本。下面把对偶处理的逻辑写清楚。%% 子问题建模给定x_opt, E_pre_opt求最坏场景 % 不确定性变量负荷 d sdpvar(numBus, numTimeSlot, full); % 运行变量DG出力、MPS出力、负荷削减 P_dg sdpvar(numBus, numTimeSlot, full); P_mps sdpvar(numBus, numMPS, numTimeSlot, full); % MPS放电 P_cur sdpvar(numBus, numTimeSlot, full); % 运行约束 Constraints_SP []; for t 1:numTimeSlot % 节点功率平衡以DistFlow线性化形式 for i 1:numBus Constraints_SP [Constraints_SP, P_dg(i,t) sum(P_mps(i,:,t)) P_cur(i,t) - d(i,t) 0]; % ... 其他潮流、电压等式约束 end % MPS功率限制 for m 1:numMPS % 只有预配置在该节点的MPS才能放电 Constraints_SP [Constraints_SP, P_mps(i,m,t) P_rate * x_opt(i,m)]; end % 负荷削减上限 Constraints_SP [Constraints_SP, 0 P_cur(i,t) d(i,t)]; end % 不确定性约束 Constraints_SP [Constraints_SP, Constraints_UNC_Linear]; % 子问题目标最大化最坏运行成本 max_d min_y (运营成本) % 使用对偶或者直接使用dualize模块 Objective_SP sum(sum(c_DG * P_dg c_cur * P_cur)); % 求解这是max-min问题需要先取对偶这里用双层求解策略 % 方案使用CCG的子问题求解方式对给定d求min再对d求max ops_sp sdpsettings(solver, gurobi, verbose, 0); [~, ~, ~, model] export(Constraints_SP, Objective_SP, ops_sp); % 导出为结构体 % 在上层max循环中对偶求解见下方说明这里有个技术难点YALMIP本身不支持直接求解max-min问题。在子问题中你需要将内层min问题取对偶再将目标函数转化为max问题此时整个子问题变成单层MILPGurobi可以直接求解。不过export函数导出的是原问题的标准形式得到对偶问题还需要手动操作。实践中我建议直接用嵌套求解法% 外层遍历/优化不确定性变量d % 内层固定d后求解最小运行成本 % 外层用Gurobi求解对偶模型内层用线性规划求解 % 伪代码逻辑 % 初始化 best_obj -inf; best_d []; % 使用外部优化器如fmincon或再次调用Gurobi在d的可行域内搜索 % 对每个候选d调用线性规划求解内层min问题这个方法虽然逻辑清晰但效率不高。推荐的做法是直接用YALMIP的robust模块%% 鲁棒优化子问题YALMIP自带的鲁棒框架 % 将不确定性声明为uncertain uncertain(d); Constraints_SP_robust [Constraints_SP, uncertain(d)]; Objective_SP_robust sum(sum(c_DG * P_dg c_cur * P_cur)); ops_sp_robust sdpsettings(solver, gurobi, verbose, 0); robust_opt optimize(Constraints_SP_robust, -Objective_SP_robust, ops_sp_robust); % 注意robust模式下默认求解的是最坏场景下的最优解不过YALMIP的robust模块会自动进行对偶重构对模型结构和不等式约束的仿射性要求严格。如果遇到Unable to convert to robust counterpart的报错说明某个约束不是仿射的需要手动重构。我在实际复现中更喜欢自己手写对偶推导。流程是先将内层min问题整理为标准LP形式(\min c^T y, ; A y \geq b, ; y \geq 0)取对偶得到(\max b^T \lambda, ; A^T \lambda \leq c, ; \lambda \geq 0)将外层max问题与对偶后的内层问题合并得到单层max问题直接求解。这个手动作法在几百个变量规模下完全可行而且完全可控不容易踩到YALMIP内部鲁棒化的各种隐藏bug。4.6 CCG迭代主循环上面两个子模块准备好后CCG主循环代码框架如下%% CCG迭代主循环 LB -inf; UB inf; iter 0; x_opt []; E_pre_opt []; d_selected []; % 已识别的场景集合 tol 1e-3; % 收敛间隙 while (UB - LB) / UB tol iter maxIter iter iter 1; fprintf( CCG 第 %d 次迭代 \n, iter); % Step 1: 求解主问题 [x_opt, E_pre_opt, LB] solve_MP(d_selected); % Step 2: 固定主问题解求解子问题 [worst_cost, d_worst, feasible] solve_SP(x_opt, E_pre_opt); if feasible UB inv_cost_total(x_opt, E_pre_opt) worst_cost; end % Step 3: 收敛判断 gap (UB - LB) / UB; fprintf( 当前间隙: %.4f | LB%.2f | UB%.2f\n, gap, LB, UB); if gap tol break; end % Step 4: 将最恶劣场景加入主问题生成新的约束和变量 d_selected{end1} d_worst; end % 输出最终结果 fprintf(优化完成: 最优MPS配置数量 %d\n, sum(sum(x_opt))); fprintf(总预配置容量 %.2f kWh\n, sum(sum(E_pre_opt)));这个循环是整个复现代码的中枢。收敛速度取决于初始场景的选择和不确定预算的大小。如果每次迭代间隙都不下降大概率是子问题求解出了问题——最常见的是对偶推导错误导致找出的最坏场景并不是真正最坏的。5. 预配置阶段与动态调度阶段的衔接逻辑文章虽然拆成了上和下两篇但在实际工程中预配置和动态调度是同一个决策链条的两个环节不是割裂的。我在复现预配置模型时就刻意预留了和动态调度衔接的接口。理清这个衔接逻辑对整个系统的理解会提升一个层次。5.1 预配置输出如何喂养动态调度预配置阶段的输出包括三样东西MPS的数量动态调度阶段的可调度资源总池子MPS的初始位置动态调度阶段移动路线的起点MPS的初始容量动态调度阶段移动时间的初始SOC上限。这三样东西构成了动态调度模型中的初始条件initial conditions。从数学上看动态调度问题就是在预配置给出的初始条件约束下滚动求解一个时间-空间耦合的调度优化问题。在Matlab代码的组织上我的建议是把预配置的输出保存为一个结构体变量% 保存预配置结果供动态调度模块调用 MPS_config struct(); MPS_config.numMPS sum(sum(x_opt)); MPS_config.sites find(sum(x_opt, 2) 0); % MPS所在的节点编号 MPS_config.capacity E_pre_opt(sum(x_opt, 2) 0, :); MPS_config.power_rate P_MPS_rate; MPS_config.initial_SOC initial_soc; % 初始荷电状态 save(MPS_config.mat, MPS_config);5.2 两阶段接口的关键假设预配置和动态调度衔接时最关键的一个隐含假设是灾害发生后MPS能否按照预配置的位置投入运行这个假设在模型中通常体现为道路通行约束。如果灾害导致某些道路中断预配置在某节点的MPS可能无法及时转移到受灾区域。有的论文把道路状态也作为不确定参数处理有的则简化为按预配置位置固定投入。复现时一定要搞清楚论文采用哪种假设因为这会直接改变模型结构和代码实现方式。我建议在代码注释中把这条假设单独标出来方便后续扩展。5.3 级联风险与MPS储备策略的对接预配置的另一个隐含输出是系统的韧性缺口resilience gap——即便在最坏场景下系统仍可能有一部分负荷无法恢复供电。这部分缺口的分布区域就是后续动态调度阶段需要重点关注的高优先级区域。实际工程中这对应着MPS的动态再调度优先级设定。我在复现预配置模型时会额外输出一个变量priority_matrix记录各节点在各时段的最坏场景负荷削减量。这个矩阵在动态调度阶段可以直接作为罚函数权重参与目标函数计算效果非常好。从论文复现的角度我强烈建议在代码中明确划分预配置主程序和动态调度接口函数两个部分不要把两个阶段的变量混在一个脚本里。否则后续改参数、换算例时代码的耦合度会让你崩溃。6. 复现中的常见报错、参数调优与避坑经验最后一部分把我在整个复现过程中踩过的坑、走过的弯路集中整理一下。这些都是写论文时不会告诉你的经验。6.1 求解器报错的根源分析报错1YALMIP报Dimension mismatch这类报错90%以上是维度不一致导致的。最常见的是某个变量定义成numBus x numTimeSlot的矩阵但在约束中参与了一个numTimeSlot x numBus的矩阵相乘。YALMIP不会自动转置报错非常直接。建议在建模前统一一个维度约定所有节点相关变量用numBus x numTimeSlot矩阵做矩阵乘法前先repmat或显式转置用size()函数检查每个变量的维度再与约束的预期维度比对。报错2Gurobi报Q is not positive semidefinite出现这个报错说明模型中存在二次项约束——通常是电压幅值 (U_{i,t}^2) 或功率与电压的乘积项。鲁棒优化和CCG模型中应该全部是线性项。你需要在模型中搜索所有^2、* U等模式把它们线性化。我的做法是在建模前做个检查% 检查约束中是否存在非线性项 nonlinear_idx find(~islinear(Constraints_SP)); if ~isempty(nonlinear_idx) warning(存在非线性约束请检查潮流建模过程); end报错3子问题找不到可行解这个错误通常意味着主问题给出的MPS配置过于激进——例如假设某节点放了MPS但移走了该节点的原有备用容量导致灾后无法满足该节点负荷。这时需要放宽MPS数量上限、增大单台容量上限或在功率平衡约束中加入虚拟切负荷变量作为惩罚项兜底。6.2 参数敏感性分析的正确打开方式预配置模型中的可调参数很多但真正值得做敏感性分析的只有三个参数敏感范围对结果的影响不确定预算 (\Gamma)2~20(\Gamma)每增1MPS总容量约增3~5%单位切负荷代价 (c_{cur})1~50(c_{cur})超过15后结果趋于饱和候选站点集合选5个与选10个差异显著直接决定MPS的空间分布形态做敏感性分析时不要只画一条曲线要把曲线画成分层图。比如 (\Gamma) 从2变化到20时同时观察MPS选址位置的变化。有时候总容量几乎没变但位置完全变了——这才是更深层的信息说明MPS预配置对不确定预算的响应方向是调整空间分布而非增加总量。6.3 算例规模扩大时的性能优化当从IEEE 33节点转到IEEE 123节点时计算时间会从分钟级跳到小时级。这时候必须做性能优化优先级从高到低减少CCG迭代次数通过调整初始场景选择用历史灾害场景做热身通常能缩短30%的迭代轮数热启动把上一轮迭代的解作为下一轮的热启动点。在Gurobi中设置ops.gurobi.Start value(x)即可收紧MIP Gap主问题求解时先设较大的MIP Gap比如5%收敛前几轮用2%最后几轮用0.1%——避免前期的冗余计算变量降维利用场景树结构将对称性较强的节点变量合并或者使用稀疏矩阵表示删除冗余约束对于远离候选站点、又无重要负荷的分支节点可以预先用潮流计算验证电压恒在安全范围从而这些节点的电压约束可以从模型中删除。我实测过光是动态MIP Gap这一个技巧就能让总计算时间缩短约40%而且是零精度损失强烈推荐。6.4 编程习惯层面的复现建议每个函数只做一件事。solve_MP、solve_SP、update_scenario、linearize_flow分开写方便单独调试分阶段测试。先只跑名义场景确认潮流收敛、目标值合理再加入不确定集最后套CCG。每加一个模块就验证一次输出保留纸质推导。把对偶问题推导完再写代码。我认识的高效率复现者桌子上都有一叠手推公式而不是直接对着论文抄代码善用value()函数。优化结束后用value()提取所有决策变量存入.mat文件。后续画图、对比分析都用这个文件不必重新跑模型。6.5 一个容易被忽视的核心细节候选站点的物理可行性最后提一个论文里不强调但实际建模必须要处理的细节候选站点除了拓扑位置还必须考虑物理可行性。在IEEE 33节点系统中某些节点可能位于道路狭窄、无场地条件的位置。如果直接把所有节点加入候选集模型算出的最优方案在实际部署时会完全不可行。我的做法是在输入数据中直接给定candidateSite向量只让模型在这几个物理可行的节点中做选择。论文中总结的经验是候选站点通常是3到8个太多会导致MPS位置分散、车组调度半径过大太少则会让选址决策失去自由度。我在复现时参考了实际供电公司的应急电源配置逻辑把候选站点选在靠近关键负荷医院、学校、大型小区、通信基站且道路通达性好的节点上出来的结果才比较合理。这篇上篇的内容到这里就完整了。梳理一下我已经完成的工作把MPS预配置问题从背景需求拆到数学模型从CCG算法讲到Matlab/YALMIP代码实现最后整理了工程复现中的常见坑和参数调优经验。每一步都是我从零开始复现后总结出来的不是简单地翻译论文摘要。最后再分享一个我实测下来的感受整个复现过程中最耗时的不是CCG算法本身而是把潮流约束以线性化形式正确地嵌入优化模型。DistFlow方程的线性化看似简单但处理不当会引入电压偏差导致算出来的预配置方案在真实潮流校验中无法通过。建议在完成优化后一定要用标准潮流计算比如Matpower的runpf对最优解进行可行性校验确认电压和功率都在合理范围内。这一步虽然简单却能挡住90%以上的复现成功但结果不可信的尴尬。预配置环节做好动态调度那部分才有扎实的输入条件。下一篇我会接着拆解MPS动态调度的滚动优化框架和路径-时间耦合建模有想一起试算的先把预配置代码跑通等下文。
