简介这份资源面向电力系统研究人员、配电网规划工程师及韧性电网方向学者围绕台风等极端灾害下配电网全过程韧性提升问题复现了基于Wasserstein距离与CVaR理论的多类型韧性资源分布鲁棒机会约束规划方法。内容涵盖不确定性模糊集构建、预防-应急-抢修双层三阶段模型、混合整数二阶锥规划转化与求解并给出完整Python代码及逐段解释便于读者理解建模思路并动手复现。资源包共1个PDF文件约790KB集中呈现论文复现说明、数学模型推导与代码实现结构紧凑适合按章节顺序研读。目前已有147人学习。读者可从中获得从理论框架到代码落地的完整链路包括模糊集与CVaR处理尾部风险的技巧、三阶段决策变量的定义方式、求解与可视化流程以及不同策略效果对比的案例分析对开展配电网韧性规划、分布式能源与储能协同配置等研究具有直接参考价值。1. 台风季配电网规划为什么“韧性”不能只算设备加固费每年台风季沿海城市配电网的抢修队伍几乎连轴转。杆塔倒伏、线路跳闸、开关站进水这些场景对一线运维来说并不陌生。但真正让规划岗头疼的不是灾后修得多快而是灾前这笔钱怎么花才不冤——全换成地下电缆当然韧性高可投资回收期长到让财务皱眉只做常规N-1校验一场十七级风圈过境多个变电站同时失压N-1准则直接失效。这就是“全过程韧性”要解决的问题不是单看某个设备扛不扛得住而是把灾害演进拆成预防、抵抗、恢复三个阶段分别量化资源投入的效果。而多类型资源指的是储能、可中断负荷、移动发电车、联络开关这些手段各有各的响应时间和空间约束不能简单加总。分布鲁棒机会约束规划则是在台风路径和强度不确定的前提下既不假设已知概率分布也不做最坏情况的绝对保守决策而是用Wasserstein距离构造一个围绕经验分布的模糊集把机会约束嵌进去再用CVaR把尾部风险折算成可优化的线性项。这套组合拳在论文里逻辑自洽但复现时参数一改结果可能天差地别。下面按我实际跑通这套模型的顺序把每个环节拆开讲。2. 从台风路径到模糊集Wasserstein距离和CVaR怎么接进配电网模型2.1 为什么不用鲁棒优化或随机规划偏要选分布鲁棒常规鲁棒优化把台风强度设成区间上限所有资源按最坏情况配置结果往往是储能容量翻倍、联络线截面加大投资高到评审会直接打回。随机规划则要求知道台风路径和强度的联合概率分布实际中样本量少则几十条历史台风记录多则几百条拟合出来的分布尾部误差极大机会约束的置信水平根本兜不住。分布鲁棒机会约束的折中逻辑是我不信你给的经验分布但我也不跑到无穷保守而是以经验分布为中心用Wasserstein距离画一个半径ε的球所有落在球内的分布都算“可能”。机会约束要求在这个球内最坏的分布下系统失负荷概率仍低于阈值。这样既保留了数据驱动又给尾部风险留了缓冲。Wasserstein距离的好处是它对样本扰动的敏感度有理论界且和CVaR有对偶关系能把无限维的分布优化转成有限维的线性规划。2.2 机会约束转CVaR把概率约束变成可求解的线性项机会约束原始形式是 ( \mathbb{P}(g(x,\xi) \leq 0) \geq 1-\alpha )其中ξ是台风随机变量x是资源规划决策。直接处理概率不可导常规做法是用CVaR近似当CVaR_{α}(g(x,ξ)) ≤ 0时机会约束成立。更精确地说对于任意分布P( \mathbb{P}(g \leq 0) \geq 1-\alpha ) 等价于 ( \text{CVaR}_{\alpha}(g) \leq 0 )。而分布鲁棒版本进一步要求模糊集内最坏分布的CVaR也≤0。Wasserstein模糊集下的最坏CVaR有对偶形式最终可以写成[ \inf_{\lambda \geq 0} \left{ \lambda \varepsilon \frac{1}{N}\sum_{i1}^{N} \max(0, g(x,\xi_i) - \lambda \cdot d(\xi_i, \xi_0)) \right} ]其中ε是Wasserstein半径N是样本数d是样本到参考点的距离。这个形式是凸的可以用现成求解器处理。实际写代码时我一般把max项用辅助变量线性化避免直接写max导致求解器报错。2.3 配电网全过程韧性的三阶段指标怎么定义预防阶段看的是灾前资源预布局是否合理指标用“关键负荷供电裕度”即台风登陆前24小时重要负荷医院、水厂、通信基站的备用电源可支撑时长。抵抗阶段看灾害发生时的网络连通性用“失负荷期望”除以总负荷但要注意这里的时间粒度要细到小时级因为台风过境通常6到12小时粗粒度会平滑掉峰值。恢复阶段看灾后恢复速度用“负荷恢复曲线下面积”与理想恢复曲线的比值越接近1越好。这三个指标不能各自独立优化否则预防阶段多投储能抵抗阶段可能因为储能放电深度限制反而降低支撑时长。我的做法是把三阶段指标加权成一个综合韧性指数权重根据当地台风频率和负荷重要性调整。沿海高发区抵抗阶段权重给到0.5内陆低发区预防阶段权重可以到0.4。2.4 多类型资源的建模差异储能、可中断负荷、移动发电车储能的核心参数是额定功率、容量、充放电效率、SOC上下限。在韧性模型里还要加一个“灾前预留SOC”变量因为台风预警后储能可能被调度到满充状态但满充会加速老化需要在目标函数里加退化成本项。可中断负荷不是简单切负荷而是分等级一级负荷不可中断二级负荷可中断但需提前通知三级负荷可即时中断。建模时用0-1变量表示中断状态但要注意通知时间约束提前通知时间越长可中断容量越大。移动发电车的关键是位置和路径。灾前预布局时发电车停在哪个变电站决定了灾后能多快接入关键负荷。我一般把候选位置设为配电网节点用运输时间矩阵约束接入时间运输时间由路网距离和台风风速共同决定风速超过25m/s时运输时间乘1.5倍。2.5 最小可复现算例IEEE 33节点加台风场景生成复现时不要一上来就上真实城市配电网先用IEEE 33节点系统跑通。台风场景生成用蒙特卡洛抽样路径用随机游走模拟强度用对数正态分布每次抽样得到一条台风轨迹和对应的风速-时间序列。然后根据风速计算线路故障概率用状态抽样得到故障场景集。import numpy as np from scipy.stats import lognorm def generate_typhoon_scenarios(n_scenarios, n_nodes, base_wind30): 生成台风场景每条场景包含各节点风速和线路故障状态 n_scenarios: 场景数建议先设50跑通 n_nodes: 配电网节点数IEEE33就是33 base_wind: 基准风速 m/s scenarios [] for _ in range(n_scenarios): # 台风路径随机游走影响各节点风速 path_offset np.cumsum(np.random.randn(n_nodes) * 0.5) wind_speed base_wind * np.exp(-0.1 * np.abs(path_offset)) # 风速超过35m/s线路故障概率显著上升 fault_prob 1 / (1 np.exp(-(wind_speed - 35) / 3)) line_status (np.random.rand(n_nodes) fault_prob).astype(int) scenarios.append({ wind: wind_speed, line_status: line_status }) return scenarios # 跑50个场景看分布 scenarios generate_typhoon_scenarios(50, 33) wind_all np.array([s[wind] for s in scenarios]) print(f风速均值 {wind_all.mean():.1f} m/s, 最大 {wind_all.max():.1f} m/s)这段代码的逻辑是先模拟台风路径偏移再折算成各节点风速最后用Sigmoid函数把风速映射成故障概率。参数base_wind根据当地历史台风调整沿海强台风区可以设到35内陆设25。n_scenarios先跑50看计算时间如果求解器能在5分钟内收敛再增加到200。注意故障概率函数不是唯一选择也可以用威布尔分布拟合历史故障数据但Sigmoid在风速阈值附近更平滑数值稳定性好。3. 分布鲁棒机会约束的代码实现从对偶形式到求解器调用3.1 模糊集半径ε的选取不能拍脑袋也不能全信理论界Wasserstein半径ε决定了模糊集的大小。ε0退化成随机规划ε→∞变成最坏情况鲁棒优化。理论上有有限样本下的置信界公式但那个界通常偏大实际用会过于保守。我的经验做法是先用理论界算一个上界然后在这个上界和0之间取0.3到0.5倍作为初始值跑一遍看失负荷概率是否接近目标置信水平。如果实际失负荷概率远低于目标说明ε偏大可以调小如果高于目标ε偏小。具体操作设目标机会约束置信水平1-α0.95即允许5%的失负荷概率。跑完模型后用测试场景集不参与训练的场景验证实际失负荷频率。如果实际频率是2%说明模型过于保守ε可以降到原来的0.7倍如果实际频率是8%ε需要增大到1.3倍。迭代两三轮基本能收敛。3.2 对偶转化后的线性化技巧辅助变量和big-M对偶形式里的max项直接写进求解器会报错因为大多数求解器不支持max函数。标准做法是引入辅助变量t_i ≥ g(x,ξ_i) - λ·d(ξ_i,ξ_0)且t_i ≥ 0目标函数里最小化λ·ε (1/N)Σt_i。这样就把max转成了线性约束。但要注意big-M的选取。如果g(x,ξ)的量级是兆瓦级big-M设太小会导致约束失效设太大会让松弛问题病态。我一般先跑一次不带big-M的模型看g(x,ξ)的最大值然后big-M取这个最大值的1.5倍。如果求解器报“numerical trouble”优先检查big-M是不是过大。import cvxpy as cp import numpy as np def solve_dro_cvar(n_scenarios, n_resources, g_values, distances, epsilon): 求解分布鲁棒CVaR近似 g_values: 每个场景下的g(x,ξ)值形状(n_scenarios,) distances: 每个场景到参考分布的距离形状(n_scenarios,) epsilon: Wasserstein半径 N n_scenarios lam cp.Variable(nonnegTrue) t cp.Variable(N, nonnegTrue) # 目标lambda*epsilon 平均t objective cp.Minimize(lam * epsilon cp.sum(t) / N) constraints [] for i in range(N): # t_i g_i - lambda * d_i constraints.append(t[i] g_values[i] - lam * distances[i]) prob cp.Problem(objective, constraints) prob.solve(solvercp.ECOS) return lam.value, t.value, prob.value # 模拟数据测试 np.random.seed(42) g_vals np.random.randn(50) * 10 20 # 模拟g值 dists np.random.rand(50) * 0.5 lam_val, t_val, obj_val solve_dro_cvar(50, 4, g_vals, dists, epsilon0.1) print(flambda{lam_val:.3f}, 目标值{obj_val:.3f})这段代码用cvxpy搭建对偶问题ECOS求解器适合小规模凸问题。参数epsilon先设0.1跑通再根据验证结果调整。g_values在实际模型中不是预先算好的而是决策变量x的函数所以需要把这段嵌到主优化问题里用x表示g。如果求解规模大ECOS可能慢可以换MOSEK或Gurobi但要注意许可证。3.3 机会约束置信水平的验证用测试集算实际失负荷频率训练完模型后必须用独立测试集验证。测试集场景数至少是训练集的1/3且要包含极端场景比如风速超过历史最大值10%。验证指标是实际失负荷频率即测试集中失负荷场景数除以总场景数。如果实际频率低于目标置信水平对应的α说明模型保守但安全如果高于α说明模糊集没覆盖住真实分布需要增大ε或增加训练样本。我一般会画一张图横轴是ε纵轴是实际失负荷频率和投资成本。随着ε增大失负荷频率下降但投资成本上升拐点通常在ε0.15到0.3之间。选拐点附近的ε性价比最高。3.4 求解器选择与收敛性排查ECOS、MOSEK、Gurobi的实测差异小规模算例IEEE 33节点50场景用ECOS足够求解时间在30秒以内。但场景数增加到200、节点数增加到123时ECOS会明显变慢甚至不收敛。这时候换MOSEK内点法对凸问题的收敛性更稳但MOSEK对模型形式要求严格max项必须线性化干净否则报“model not convex”。Gurobi适合混合整数问题因为配电网里可中断负荷和移动发电车有0-1变量Gurobi的分支定界比ECOS强。但Gurobi的学术许可证申请需要时间商用许可证贵。如果只是复现论文ECOS加MOSEK组合基本够用。收敛性排查顺序先检查big-M是否过大再检查Wasserstein距离计算是否有数值误差样本距离矩阵如果出现NaN后面全崩最后检查CVaR的α是否设得太小α0.01时尾部样本太少对偶问题不稳定。4. 复现避坑Wasserstein距离计算、场景生成和参数标定里的翻车点4.1 现象Wasserstein距离矩阵出现NaN求解器直接报错原因样本距离计算时用了欧氏距离但台风场景里风速和故障状态量纲差几个数量级风速是几十故障状态是0或1直接算距离会导致数值下溢。另外如果两个样本完全相同距离为0后续除以距离的步骤会出Inf。解决先对每个维度做归一化风速除以最大风速故障状态保持0-1。距离用加权欧氏距离权重按维度重要性设风速权重0.7故障状态权重0.3。如果还有0距离加一个极小值1e-6。4.2 现象机会约束置信水平设0.95但实际失负荷频率高达15%原因训练场景太少经验分布根本没覆盖尾部。50个场景里可能只有2个极端场景Wasserstein球半径再大也兜不住。另外α设0.05意味着允许5%失负荷但实际台风频率可能高于5%模型没考虑台风发生概率本身。解决训练场景增加到至少200且要分层抽样极端场景占20%。α根据当地台风年发生概率调整如果每年平均1.5次台风α可以设0.1留更多裕度。4.3 现象储能容量优化结果比常规规划还小韧性反而下降原因目标函数里储能退化成本权重设太大优化器为了省退化成本宁愿少投储能。或者灾前预留SOC约束没加储能被调度到深度放电灾时没电可用。解决退化成本权重先设0跑一遍看储能容量上限再逐步增加权重到容量下降10%左右。灾前预留SOC约束必须硬约束SOC_min在台风预警后提高到0.8。4.4 现象移动发电车预布局结果全挤在同一个变电站原因运输时间矩阵没考虑路网容量所有发电车都选最近的变电站但实际路网可能拥堵。另外目标函数里发电车接入时间权重太低优化器觉得放哪都一样。解决运输时间矩阵加拥堵系数台风期间主干道拥堵系数1.8。目标函数里接入时间权重提高到和失负荷成本同一量级。还可以加地理分散约束同一变电站最多停2辆。4.5 现象CVaR对偶问题求解时间随场景数指数增长原因对偶形式里每个场景一个辅助变量t_i场景数N500时变量数500约束数500内点法复杂度O(N^3)。如果还嵌在主问题里整体规模更大。解决场景削减。用快速前向选择法把500场景削减到100保留尾部极端场景。削减后再跑求解时间从小时级降到分钟级。削减时注意保留至少10%的极端场景否则尾部风险被削掉。5. 进阶技巧用场景削减和对偶校准把复现时间压到可接受范围场景削减是复现这套模型最实用的加速手段。快速前向选择法的逻辑是先选一个代表性场景然后每次选一个使Wasserstein距离增量最大的场景直到达到目标场景数。这样保留的场景在分布上最分散尾部场景不会被均匀采样稀释。def fast_forward_selection(scenarios, n_keep): 快速前向选择场景削减 scenarios: 原始场景列表每个场景是字典 n_keep: 保留场景数 n len(scenarios) # 提取特征向量 features np.array([np.concatenate([s[wind], s[line_status]]) for s in scenarios]) selected [0] # 先选第一个 candidates list(range(1, n)) while len(selected) n_keep: best_idx None best_dist -1 for c in candidates: # 计算候选场景到已选场景集的最小距离 min_dist min(np.linalg.norm(features[c] - features[s]) for s in selected) if min_dist best_dist: best_dist min_dist best_idx c selected.append(best_idx) candidates.remove(best_idx) return [scenarios[i] for i in selected] # 从200削减到80 reduced fast_forward_selection(scenarios, 80) print(f削减后场景数: {len(reduced)})这段代码的核心是贪心选择距离已选集合最远的场景保证保留的场景在特征空间里尽量分散。n_keep一般设原始场景数的30%到40%200场景保留80个求解时间能降60%以上。注意特征向量要先归一化否则风速会主导距离计算。对偶校准是另一个技巧。Wasserstein半径ε和对偶变量λ在最优解处满足互补松弛条件如果λ0说明模糊集约束不起作用ε可以调小如果λ很大说明模糊集约束紧ε可能偏小。我一般跑完模型后检查λ值如果λ在0.1到1之间说明ε选得合理如果λ10ε需要增大如果λ0.01ε可以减小。最后说一个我踩过的坑不要一上来就调参先把模型在50场景、ε0.1下跑通确认目标函数值和约束满足再逐步增加场景数和调整ε。很多复现失败不是模型错是数值问题在早期没暴露到大规模时集中爆发。希望帮到你。本文还有配套的精品资源点击获取 台风季配电网韧性规划：分布鲁棒机会约束与多资源协同复现
简介这份资源面向电力系统研究人员、配电网规划工程师及韧性电网方向学者围绕台风等极端灾害下配电网全过程韧性提升问题复现了基于Wasserstein距离与CVaR理论的多类型韧性资源分布鲁棒机会约束规划方法。内容涵盖不确定性模糊集构建、预防-应急-抢修双层三阶段模型、混合整数二阶锥规划转化与求解并给出完整Python代码及逐段解释便于读者理解建模思路并动手复现。资源包共1个PDF文件约790KB集中呈现论文复现说明、数学模型推导与代码实现结构紧凑适合按章节顺序研读。目前已有147人学习。读者可从中获得从理论框架到代码落地的完整链路包括模糊集与CVaR处理尾部风险的技巧、三阶段决策变量的定义方式、求解与可视化流程以及不同策略效果对比的案例分析对开展配电网韧性规划、分布式能源与储能协同配置等研究具有直接参考价值。1. 台风季配电网规划为什么“韧性”不能只算设备加固费每年台风季沿海城市配电网的抢修队伍几乎连轴转。杆塔倒伏、线路跳闸、开关站进水这些场景对一线运维来说并不陌生。但真正让规划岗头疼的不是灾后修得多快而是灾前这笔钱怎么花才不冤——全换成地下电缆当然韧性高可投资回收期长到让财务皱眉只做常规N-1校验一场十七级风圈过境多个变电站同时失压N-1准则直接失效。这就是“全过程韧性”要解决的问题不是单看某个设备扛不扛得住而是把灾害演进拆成预防、抵抗、恢复三个阶段分别量化资源投入的效果。而多类型资源指的是储能、可中断负荷、移动发电车、联络开关这些手段各有各的响应时间和空间约束不能简单加总。分布鲁棒机会约束规划则是在台风路径和强度不确定的前提下既不假设已知概率分布也不做最坏情况的绝对保守决策而是用Wasserstein距离构造一个围绕经验分布的模糊集把机会约束嵌进去再用CVaR把尾部风险折算成可优化的线性项。这套组合拳在论文里逻辑自洽但复现时参数一改结果可能天差地别。下面按我实际跑通这套模型的顺序把每个环节拆开讲。2. 从台风路径到模糊集Wasserstein距离和CVaR怎么接进配电网模型2.1 为什么不用鲁棒优化或随机规划偏要选分布鲁棒常规鲁棒优化把台风强度设成区间上限所有资源按最坏情况配置结果往往是储能容量翻倍、联络线截面加大投资高到评审会直接打回。随机规划则要求知道台风路径和强度的联合概率分布实际中样本量少则几十条历史台风记录多则几百条拟合出来的分布尾部误差极大机会约束的置信水平根本兜不住。分布鲁棒机会约束的折中逻辑是我不信你给的经验分布但我也不跑到无穷保守而是以经验分布为中心用Wasserstein距离画一个半径ε的球所有落在球内的分布都算“可能”。机会约束要求在这个球内最坏的分布下系统失负荷概率仍低于阈值。这样既保留了数据驱动又给尾部风险留了缓冲。Wasserstein距离的好处是它对样本扰动的敏感度有理论界且和CVaR有对偶关系能把无限维的分布优化转成有限维的线性规划。2.2 机会约束转CVaR把概率约束变成可求解的线性项机会约束原始形式是 ( \mathbb{P}(g(x,\xi) \leq 0) \geq 1-\alpha )其中ξ是台风随机变量x是资源规划决策。直接处理概率不可导常规做法是用CVaR近似当CVaR_{α}(g(x,ξ)) ≤ 0时机会约束成立。更精确地说对于任意分布P( \mathbb{P}(g \leq 0) \geq 1-\alpha ) 等价于 ( \text{CVaR}_{\alpha}(g) \leq 0 )。而分布鲁棒版本进一步要求模糊集内最坏分布的CVaR也≤0。Wasserstein模糊集下的最坏CVaR有对偶形式最终可以写成[ \inf_{\lambda \geq 0} \left{ \lambda \varepsilon \frac{1}{N}\sum_{i1}^{N} \max(0, g(x,\xi_i) - \lambda \cdot d(\xi_i, \xi_0)) \right} ]其中ε是Wasserstein半径N是样本数d是样本到参考点的距离。这个形式是凸的可以用现成求解器处理。实际写代码时我一般把max项用辅助变量线性化避免直接写max导致求解器报错。2.3 配电网全过程韧性的三阶段指标怎么定义预防阶段看的是灾前资源预布局是否合理指标用“关键负荷供电裕度”即台风登陆前24小时重要负荷医院、水厂、通信基站的备用电源可支撑时长。抵抗阶段看灾害发生时的网络连通性用“失负荷期望”除以总负荷但要注意这里的时间粒度要细到小时级因为台风过境通常6到12小时粗粒度会平滑掉峰值。恢复阶段看灾后恢复速度用“负荷恢复曲线下面积”与理想恢复曲线的比值越接近1越好。这三个指标不能各自独立优化否则预防阶段多投储能抵抗阶段可能因为储能放电深度限制反而降低支撑时长。我的做法是把三阶段指标加权成一个综合韧性指数权重根据当地台风频率和负荷重要性调整。沿海高发区抵抗阶段权重给到0.5内陆低发区预防阶段权重可以到0.4。2.4 多类型资源的建模差异储能、可中断负荷、移动发电车储能的核心参数是额定功率、容量、充放电效率、SOC上下限。在韧性模型里还要加一个“灾前预留SOC”变量因为台风预警后储能可能被调度到满充状态但满充会加速老化需要在目标函数里加退化成本项。可中断负荷不是简单切负荷而是分等级一级负荷不可中断二级负荷可中断但需提前通知三级负荷可即时中断。建模时用0-1变量表示中断状态但要注意通知时间约束提前通知时间越长可中断容量越大。移动发电车的关键是位置和路径。灾前预布局时发电车停在哪个变电站决定了灾后能多快接入关键负荷。我一般把候选位置设为配电网节点用运输时间矩阵约束接入时间运输时间由路网距离和台风风速共同决定风速超过25m/s时运输时间乘1.5倍。2.5 最小可复现算例IEEE 33节点加台风场景生成复现时不要一上来就上真实城市配电网先用IEEE 33节点系统跑通。台风场景生成用蒙特卡洛抽样路径用随机游走模拟强度用对数正态分布每次抽样得到一条台风轨迹和对应的风速-时间序列。然后根据风速计算线路故障概率用状态抽样得到故障场景集。import numpy as np from scipy.stats import lognorm def generate_typhoon_scenarios(n_scenarios, n_nodes, base_wind30): 生成台风场景每条场景包含各节点风速和线路故障状态 n_scenarios: 场景数建议先设50跑通 n_nodes: 配电网节点数IEEE33就是33 base_wind: 基准风速 m/s scenarios [] for _ in range(n_scenarios): # 台风路径随机游走影响各节点风速 path_offset np.cumsum(np.random.randn(n_nodes) * 0.5) wind_speed base_wind * np.exp(-0.1 * np.abs(path_offset)) # 风速超过35m/s线路故障概率显著上升 fault_prob 1 / (1 np.exp(-(wind_speed - 35) / 3)) line_status (np.random.rand(n_nodes) fault_prob).astype(int) scenarios.append({ wind: wind_speed, line_status: line_status }) return scenarios # 跑50个场景看分布 scenarios generate_typhoon_scenarios(50, 33) wind_all np.array([s[wind] for s in scenarios]) print(f风速均值 {wind_all.mean():.1f} m/s, 最大 {wind_all.max():.1f} m/s)这段代码的逻辑是先模拟台风路径偏移再折算成各节点风速最后用Sigmoid函数把风速映射成故障概率。参数base_wind根据当地历史台风调整沿海强台风区可以设到35内陆设25。n_scenarios先跑50看计算时间如果求解器能在5分钟内收敛再增加到200。注意故障概率函数不是唯一选择也可以用威布尔分布拟合历史故障数据但Sigmoid在风速阈值附近更平滑数值稳定性好。3. 分布鲁棒机会约束的代码实现从对偶形式到求解器调用3.1 模糊集半径ε的选取不能拍脑袋也不能全信理论界Wasserstein半径ε决定了模糊集的大小。ε0退化成随机规划ε→∞变成最坏情况鲁棒优化。理论上有有限样本下的置信界公式但那个界通常偏大实际用会过于保守。我的经验做法是先用理论界算一个上界然后在这个上界和0之间取0.3到0.5倍作为初始值跑一遍看失负荷概率是否接近目标置信水平。如果实际失负荷概率远低于目标说明ε偏大可以调小如果高于目标ε偏小。具体操作设目标机会约束置信水平1-α0.95即允许5%的失负荷概率。跑完模型后用测试场景集不参与训练的场景验证实际失负荷频率。如果实际频率是2%说明模型过于保守ε可以降到原来的0.7倍如果实际频率是8%ε需要增大到1.3倍。迭代两三轮基本能收敛。3.2 对偶转化后的线性化技巧辅助变量和big-M对偶形式里的max项直接写进求解器会报错因为大多数求解器不支持max函数。标准做法是引入辅助变量t_i ≥ g(x,ξ_i) - λ·d(ξ_i,ξ_0)且t_i ≥ 0目标函数里最小化λ·ε (1/N)Σt_i。这样就把max转成了线性约束。但要注意big-M的选取。如果g(x,ξ)的量级是兆瓦级big-M设太小会导致约束失效设太大会让松弛问题病态。我一般先跑一次不带big-M的模型看g(x,ξ)的最大值然后big-M取这个最大值的1.5倍。如果求解器报“numerical trouble”优先检查big-M是不是过大。import cvxpy as cp import numpy as np def solve_dro_cvar(n_scenarios, n_resources, g_values, distances, epsilon): 求解分布鲁棒CVaR近似 g_values: 每个场景下的g(x,ξ)值形状(n_scenarios,) distances: 每个场景到参考分布的距离形状(n_scenarios,) epsilon: Wasserstein半径 N n_scenarios lam cp.Variable(nonnegTrue) t cp.Variable(N, nonnegTrue) # 目标lambda*epsilon 平均t objective cp.Minimize(lam * epsilon cp.sum(t) / N) constraints [] for i in range(N): # t_i g_i - lambda * d_i constraints.append(t[i] g_values[i] - lam * distances[i]) prob cp.Problem(objective, constraints) prob.solve(solvercp.ECOS) return lam.value, t.value, prob.value # 模拟数据测试 np.random.seed(42) g_vals np.random.randn(50) * 10 20 # 模拟g值 dists np.random.rand(50) * 0.5 lam_val, t_val, obj_val solve_dro_cvar(50, 4, g_vals, dists, epsilon0.1) print(flambda{lam_val:.3f}, 目标值{obj_val:.3f})这段代码用cvxpy搭建对偶问题ECOS求解器适合小规模凸问题。参数epsilon先设0.1跑通再根据验证结果调整。g_values在实际模型中不是预先算好的而是决策变量x的函数所以需要把这段嵌到主优化问题里用x表示g。如果求解规模大ECOS可能慢可以换MOSEK或Gurobi但要注意许可证。3.3 机会约束置信水平的验证用测试集算实际失负荷频率训练完模型后必须用独立测试集验证。测试集场景数至少是训练集的1/3且要包含极端场景比如风速超过历史最大值10%。验证指标是实际失负荷频率即测试集中失负荷场景数除以总场景数。如果实际频率低于目标置信水平对应的α说明模型保守但安全如果高于α说明模糊集没覆盖住真实分布需要增大ε或增加训练样本。我一般会画一张图横轴是ε纵轴是实际失负荷频率和投资成本。随着ε增大失负荷频率下降但投资成本上升拐点通常在ε0.15到0.3之间。选拐点附近的ε性价比最高。3.4 求解器选择与收敛性排查ECOS、MOSEK、Gurobi的实测差异小规模算例IEEE 33节点50场景用ECOS足够求解时间在30秒以内。但场景数增加到200、节点数增加到123时ECOS会明显变慢甚至不收敛。这时候换MOSEK内点法对凸问题的收敛性更稳但MOSEK对模型形式要求严格max项必须线性化干净否则报“model not convex”。Gurobi适合混合整数问题因为配电网里可中断负荷和移动发电车有0-1变量Gurobi的分支定界比ECOS强。但Gurobi的学术许可证申请需要时间商用许可证贵。如果只是复现论文ECOS加MOSEK组合基本够用。收敛性排查顺序先检查big-M是否过大再检查Wasserstein距离计算是否有数值误差样本距离矩阵如果出现NaN后面全崩最后检查CVaR的α是否设得太小α0.01时尾部样本太少对偶问题不稳定。4. 复现避坑Wasserstein距离计算、场景生成和参数标定里的翻车点4.1 现象Wasserstein距离矩阵出现NaN求解器直接报错原因样本距离计算时用了欧氏距离但台风场景里风速和故障状态量纲差几个数量级风速是几十故障状态是0或1直接算距离会导致数值下溢。另外如果两个样本完全相同距离为0后续除以距离的步骤会出Inf。解决先对每个维度做归一化风速除以最大风速故障状态保持0-1。距离用加权欧氏距离权重按维度重要性设风速权重0.7故障状态权重0.3。如果还有0距离加一个极小值1e-6。4.2 现象机会约束置信水平设0.95但实际失负荷频率高达15%原因训练场景太少经验分布根本没覆盖尾部。50个场景里可能只有2个极端场景Wasserstein球半径再大也兜不住。另外α设0.05意味着允许5%失负荷但实际台风频率可能高于5%模型没考虑台风发生概率本身。解决训练场景增加到至少200且要分层抽样极端场景占20%。α根据当地台风年发生概率调整如果每年平均1.5次台风α可以设0.1留更多裕度。4.3 现象储能容量优化结果比常规规划还小韧性反而下降原因目标函数里储能退化成本权重设太大优化器为了省退化成本宁愿少投储能。或者灾前预留SOC约束没加储能被调度到深度放电灾时没电可用。解决退化成本权重先设0跑一遍看储能容量上限再逐步增加权重到容量下降10%左右。灾前预留SOC约束必须硬约束SOC_min在台风预警后提高到0.8。4.4 现象移动发电车预布局结果全挤在同一个变电站原因运输时间矩阵没考虑路网容量所有发电车都选最近的变电站但实际路网可能拥堵。另外目标函数里发电车接入时间权重太低优化器觉得放哪都一样。解决运输时间矩阵加拥堵系数台风期间主干道拥堵系数1.8。目标函数里接入时间权重提高到和失负荷成本同一量级。还可以加地理分散约束同一变电站最多停2辆。4.5 现象CVaR对偶问题求解时间随场景数指数增长原因对偶形式里每个场景一个辅助变量t_i场景数N500时变量数500约束数500内点法复杂度O(N^3)。如果还嵌在主问题里整体规模更大。解决场景削减。用快速前向选择法把500场景削减到100保留尾部极端场景。削减后再跑求解时间从小时级降到分钟级。削减时注意保留至少10%的极端场景否则尾部风险被削掉。5. 进阶技巧用场景削减和对偶校准把复现时间压到可接受范围场景削减是复现这套模型最实用的加速手段。快速前向选择法的逻辑是先选一个代表性场景然后每次选一个使Wasserstein距离增量最大的场景直到达到目标场景数。这样保留的场景在分布上最分散尾部场景不会被均匀采样稀释。def fast_forward_selection(scenarios, n_keep): 快速前向选择场景削减 scenarios: 原始场景列表每个场景是字典 n_keep: 保留场景数 n len(scenarios) # 提取特征向量 features np.array([np.concatenate([s[wind], s[line_status]]) for s in scenarios]) selected [0] # 先选第一个 candidates list(range(1, n)) while len(selected) n_keep: best_idx None best_dist -1 for c in candidates: # 计算候选场景到已选场景集的最小距离 min_dist min(np.linalg.norm(features[c] - features[s]) for s in selected) if min_dist best_dist: best_dist min_dist best_idx c selected.append(best_idx) candidates.remove(best_idx) return [scenarios[i] for i in selected] # 从200削减到80 reduced fast_forward_selection(scenarios, 80) print(f削减后场景数: {len(reduced)})这段代码的核心是贪心选择距离已选集合最远的场景保证保留的场景在特征空间里尽量分散。n_keep一般设原始场景数的30%到40%200场景保留80个求解时间能降60%以上。注意特征向量要先归一化否则风速会主导距离计算。对偶校准是另一个技巧。Wasserstein半径ε和对偶变量λ在最优解处满足互补松弛条件如果λ0说明模糊集约束不起作用ε可以调小如果λ很大说明模糊集约束紧ε可能偏小。我一般跑完模型后检查λ值如果λ在0.1到1之间说明ε选得合理如果λ10ε需要增大如果λ0.01ε可以减小。最后说一个我踩过的坑不要一上来就调参先把模型在50场景、ε0.1下跑通确认目标函数值和约束满足再逐步增加场景数和调整ε。很多复现失败不是模型错是数值问题在早期没暴露到大规模时集中爆发。希望帮到你。本文还有配套的精品资源点击获取