配电网全过程韧性规划：Wasserstein模糊集与CVaR分布鲁棒优化实战

简介：这份资源面向电力系统研究人员、配电网规划工程师及韧性电网方向学者，针对当前韧性提升策略单阶段关注、资源类型有限、线路故障不确定性刻画不足等痛点，复现了基于Wasserstein距离与CVaR理论的分布鲁棒机会约束规划方法。资源包内含1个PDF文件，约790KB，集中呈现论文复现的完整代码与逐段解释，涵盖不确定性模糊集构建、预防-应急-抢修双层三阶段模型、混合整数二阶锥规划转化及结果可视化等环节。读者可据此理解多类型韧性资源协同配置的建模思路，掌握分布鲁棒优化处理极端事件尾部风险的关键技巧，并借助可运行代码完成案例复现与参数调整。目前已有147人学习，适合希望深入智能电网韧性规划的中高级技术读者参考。

1. 台风季里的配电网规划：为什么“韧性”不能只算期望值

每年台风季，沿海城市的配电网调度员最怕的不是“平均停电时间”这个指标，而是那种十年一遇的极端天气把几条主干线同时打断、抢修队被积水堵在路上、备用电源撑不到天亮的场景。传统规划模型算的是期望成本，可极端灾害恰恰落在概率分布的尾部——期望值根本管不住它。这篇要拆解的，就是一套把“配电网全过程韧性”拆成灾前、灾中、灾后三个阶段，用 Wasserstein 距离构造模糊集、用 CVaR 度量尾部风险、再套上机会约束的多类型资源分布鲁棒规划方法。它解决的核心问题是：在台风路径和故障概率都不确定的情况下，储能、联络开关、移动发电车、可中断负荷这几类资源到底该配多少、配在哪，才能让最坏情况下的韧性成本可控。适合做配电网规划、韧性评估、鲁棒优化方向的工程师和研究生，尤其是需要复现论文又卡在“模糊集怎么建、对偶怎么推”这一步的人。

2. 全过程韧性建模：三个阶段的目标函数怎么拆

2.1 灾前、灾中、灾后各自管什么

配电网韧性（resilience）和可靠性（reliability）最大的区别在于时间尺度和事件强度。可靠性管的是高频小扰动，韧性管的是低频高损事件。把全过程拆开，是因为三个阶段的可控变量完全不同。

灾前阶段（pre-event）是规划视角，决策变量是资源容量和位置：储能装多大、移动发电车预置在哪个节点、哪些线路加装自动化开关。这个阶段的成本是投资成本，目标是让系统在灾害来临前就具备足够的“缓冲”。

灾中阶段（during-event）是运行视角，台风已经登陆，线路开始故障。此时决策变量是拓扑重构、孤岛划分、储能放电策略。这个阶段的核心约束是功率平衡和辐射状拓扑，目标是最小化失负荷量。

灾后阶段（post-event）是恢复视角，故障线路逐步修复，决策变量是抢修顺序和负荷恢复顺序。目标是最小化恢复时间和累计失负荷。

三个阶段耦合在一起，灾前的资源配置决定了灾中能撑多久，灾中的运行策略决定了灾后恢复的起点。常见做法是用场景法把三个阶段串起来，但场景法需要知道故障概率分布——而台风场景下，这个分布恰恰是最不确定的。

2.2 为什么用 Wasserstein 距离而不是经验分布

如果直接用历史台风数据生成的经验分布做优化，会有一个致命问题：样本外表现极差。历史数据里没出现过的故障组合，模型完全没准备。这就是“过拟合”在规划问题里的翻车现场。

Wasserstein 距离的作用是构造一个以经验分布为中心的“模糊集”（ambiguity set），把所有与经验分布距离不超过 ε 的分布都纳入考虑。数学上，1-Wasserstein 距离定义为：

W(P, P̂) = inf_{π∈Π(P,P̂)} ∫ d(x,y) dπ(x,y)

其中 P̂ 是经验分布，P 是真实分布，d(x,y) 是场景之间的距离度量。模糊集就是 {P : W(P, P̂) ≤ ε}。ε 越大，模型越保守；ε 越小，越接近经验分布。

选 Wasserstein 而不是 KL 散度或 φ-散度的理由很实际：Wasserstein 距离对支撑集不重叠的分布仍然有定义，而 KL 散度在支撑集不重叠时会发散。台风场景下，极端故障组合的支撑集和训练集很可能不重叠，KL 直接失效。

2.3 机会约束和 CVaR 怎么配合

机会约束的形式是 P(失负荷 ≤ L_max) ≥ 1-α，意思是失负荷不超过阈值的概率至少是 1-α。但机会约束本身是非凸的，直接求解很困难。CVaR（条件风险价值）提供了一个凸近似。

CVaR 的定义是：CVaR_α(X) = E[X | X ≥ VaR_α(X)]，即在最坏的 α 比例场景下的平均损失。机会约束 P(X ≤ L_max) ≥ 1-α 可以转化为 CVaR 约束：CVaR_α(X) ≤ L_max。这个转化是保守的，但保证了凸性，可以用线性规划求解。

在分布鲁棒框架下，CVaR 要在模糊集内取最坏情况：

sup_{P∈模糊集} CVaR_α^P(X) ≤ L_max

这个 sup 问题可以通过对偶理论转化成有限维的线性规划，这是整篇论文最核心的推导。

2.4 多类型资源的建模差异

储能、联络开关、移动发电车、可中断负荷这四类资源在模型里的处理方式完全不同：

储能的核心约束是 SOC 连续性：E(t+1) = E(t) + η_ch * P_ch(t) - P_dis(t)/η_dis。联络开关的核心约束是拓扑辐射状，通常用单商品流或生成树约束来保证。移动发电车的难点在于调度路径和到达时间的耦合，常见做法是把路径简化为“预置点到故障点的最短距离”。可中断负荷需要建模用户响应意愿，通常用价格弹性或中断成本函数。

3. 分布鲁棒对偶推导：从 sup 到线性规划

3.1 模糊集的定义和 ε 的选取

模糊集的形式是：

P = {P : W_1(P, P̂_N) ≤ ε}

其中 P̂_N = (1/N) Σ δ_{ξ_i} 是 N 个历史场景的经验分布，δ 是狄拉克测度。ε 的选取直接决定保守程度。常见做法有两种：一是用理论公式 ε = C * sqrt(1/N) * (1 + log(1/β))，其中 β 是置信水平；二是用交叉验证，在验证集上调整 ε 使目标函数最优。

我一般会先用理论公式给一个初值，再在 [0.5ε_0, 2ε_0] 范围内做敏感性分析。如果 ε 太小，模型退化成随机规划，样本外表现差；如果 ε 太大，投资成本飙升，方案没有经济性。

3.2 对偶转化的关键步骤

原始问题是：

min_x sup_{P∈P} E_P[CVaR_α(f(x, ξ))]

其中 f(x, ξ) 是失负荷函数，x 是规划决策，ξ 是随机场景。CVaR 可以写成：

CVaR_α(f) = min_{t} { t + (1/α) E_P[(f - t)^+] }

代入后，sup 和 min 交换顺序（强对偶成立），得到：

min_{x, t} { t + (1/α) sup_{P∈P} E_P[(f - t)^+] }

关键一步：sup_{P∈P} E_P[g(ξ)] 在 Wasserstein 模糊集下有对偶形式：

sup_{P: W(P, P̂) ≤ ε} E_P[g] = inf_{λ≥0} { λε + (1/N) Σ_i sup_ξ [g(ξ) - λ d(ξ, ξ_i)] }

这个对偶公式把无限维的分布优化转化成了有限维的凸优化。λ 是对偶变量，可以理解为“模糊集半径的影子价格”。

3.3 代码实现：用 Python 构造模糊集和对偶问题

import numpy as np from scipy.optimize import linprog def wasserstein_dual(scenarios, g_values, epsilon, distance_matrix): """ 求解 Wasserstein 模糊集下的对偶问题 scenarios: 历史场景 (N, T) g_values: 每个场景下的损失值 (N,) epsilon: 模糊集半径 distance_matrix: 场景间距离 (N, N) """ N = len(scenarios) # 对偶变量: lambda (1个) + s_i (N个) # 目标: min lambda * epsilon + (1/N) * sum(s_i) c = np.concatenate([[epsilon], np.ones(N) / N]) # 约束: s_i >= g_j - lambda * d(i,j) 对所有 i,j # 即: s_i + lambda * d(i,j) >= g_j A_ub = [] b_ub = [] for i in range(N): for j in range(N): row = np.zeros(N + 1) row[0] = distance_matrix[i, j] # lambda 系数 row[i + 1] = 1.0 # s_i 系数 A_ub.append(-row) # 转为 <= 形式 b_ub.append(-g_values[j]) # 变量下界: lambda >= 0, s_i 无约束 bounds = [(0, None)] + [(None, None)] * N result = linprog(c, A_ub=np.array(A_ub), b_ub=np.array(b_ub), bounds=bounds, method='highs') if result.success: lambda_opt = result.x[0] s_opt = result.x[1:] return lambda_opt, s_opt, result.fun else: raise ValueError(f"对偶问题求解失败: {result.message}") # 示例: 10个场景, 距离用欧氏距离 np.random.seed(42) N = 10 scenarios = np.random.randn(N, 3) # 3维场景 g_values = np.random.rand(N) * 100 # 损失值 # 计算距离矩阵 distance_matrix = np.zeros((N, N)) for i in range(N): for j in range(N): distance_matrix[i, j] = np.linalg.norm(scenarios[i] - scenarios[j]) epsilon = 0.5 lambda_opt, s_opt, obj = wasserstein_dual(scenarios, g_values, epsilon, distance_matrix) print(f"最优 lambda: {lambda_opt:.4f}") print(f"对偶目标值: {obj:.4f}")

这段代码的核心逻辑是：把 sup_{P∈P} E_P[g] 的对偶形式写成线性规划。变量是 λ 和 s_i，约束是 s_i + λ d(i,j) ≥ g_j。求解后得到的 λ 就是模糊集的“影子价格”，s_i 是每个场景的调整后损失。

参数说明：epsilon 控制保守程度，一般取 0.1 到 1.0 之间；distance_matrix 的度量方式很关键，常见做法是用场景向量的欧氏距离，但如果场景是故障组合，可以用汉明距离或加权距离。g_values 是每个场景下的损失，需要先固定规划决策 x 才能计算。

3.4 机会约束的 CVaR 近似误差分析

CVaR 近似机会约束是保守的，保守程度取决于 α 的选取。理论上，如果 CVaR_α(X) ≤ L_max，那么 P(X ≤ L_max) ≥ 1-α。但反过来不成立：P(X ≤ L_max) ≥ 1-α 不能推出 CVaR_α(X) ≤ L_max。

实际中，α 一般取 0.05 到 0.1，对应 90% 到 95% 的置信水平。如果 α 太小，CVaR 只关注极端尾部，模型会过度保守；如果 α 太大，CVaR 接近期望值，失去尾部控制能力。我一般会做 α 的敏感性分析，看目标函数和失负荷概率的变化曲线，选一个拐点。

4. 求解流程：从数据到规划方案的完整链路

4.1 场景生成和距离矩阵构造

台风场景的生成需要三个维度的数据：台风路径、风速衰减、线路故障概率。常见做法是用蒙特卡洛模拟，但纯随机采样效率低。我一般会用拉丁超立方采样（LHS）先生成均匀覆盖的场景，再用重要性采样把极端场景的权重调高。

from scipy.stats import qmc import numpy as np def generate_typhoon_scenarios(n_scenarios, n_dims, seed=42): """ 用拉丁超立方采样生成台风场景 n_scenarios: 场景数 n_dims: 维度 (路径x, 路径y, 风速, 持续时间...) """ sampler = qmc.LatinHypercube(d=n_dims, seed=seed) samples = sampler.random(n=n_scenarios) # 假设前两维是路径坐标, 第三维是风速, 第四维是持续时间 # 路径坐标范围 [-100, 100] km samples[:, 0] = samples[:, 0] * 200 - 100 samples[:, 1] = samples[:, 1] * 200 - 100 # 风速范围 [10, 60] m/s samples[:, 2] = samples[:, 2] * 50 + 10 # 持续时间 [1, 12] 小时 samples[:, 3] = samples[:, 3] * 11 + 1 return samples def compute_distance_matrix(scenarios, weights=None): """ 计算场景间的加权欧氏距离 """ n = len(scenarios) dist = np.zeros((n, n)) for i in range(n): for j in range(n): diff = scenarios[i] - scenarios[j] dist[i, j] = np.sqrt(np.sum(diff**2)) return dist scenarios = generate_typhoon_scenarios(50, 4) dist_matrix = compute_distance_matrix(scenarios) print(f"场景形状: {scenarios.shape}") print(f"距离矩阵范围: [{dist_matrix.min():.2f}, {dist_matrix.max():.2f}]")

LHS 的优势是用较少的场景数覆盖整个参数空间。50 个场景在 4 维空间里已经能给出不错的覆盖，比纯随机采样效率高 3 到 5 倍。距离矩阵用欧氏距离是最简单的做法，但如果各维度量纲差异大，需要先归一化。

4.2 主问题-子问题迭代框架

分布鲁棒优化通常用 Benders 分解或列与约束生成（C&CG）算法求解。主问题是规划决策，子问题是给定规划决策下的最坏情况评估。

主问题： min c^T x + θ s.t. θ ≥ 对偶子问题的目标值（逐步添加割平面）

子问题： max_{P∈P} E_P[CVaR_α(f(x, ξ))]

迭代流程：先给一个初始规划方案 x_0，求解子问题得到最坏情况损失，把对应的割平面加回主问题，再求解主问题得到新的 x，直到上下界收敛。

def ccg_algorithm(max_iter=20, tol=1e-4): """ C&CG 主循环框架 """ # 初始化 x = np.zeros(10) # 规划决策: 储能容量等 upper_bound = np.inf lower_bound = -np.inf for k in range(max_iter): # 步骤1: 求解子问题, 得到最坏情况损失 worst_loss, worst_scenario = solve_subproblem(x) # 更新上界 upper_bound = min(upper_bound, worst_loss) # 步骤2: 添加割平面到主问题 # 割平面: theta >= worst_loss - lambda * (x - x_current) add_cut_to_master(worst_loss, worst_scenario, x) # 步骤3: 求解主问题 x_new, theta_new = solve_master() lower_bound = max(lower_bound, theta_new) # 收敛判断 gap = (upper_bound - lower_bound) / abs(upper_bound) print(f"迭代 {k+1}: 上界={upper_bound:.4f}, 下界={lower_bound:.4f}, gap={gap:.6f}") if gap < tol: print(f"收敛于第 {k+1} 次迭代") break x = x_new return x, upper_bound def solve_subproblem(x): """求解给定规划决策下的最坏情况""" # 这里需要根据具体模型实现 # 返回最坏情况损失和对应场景 pass def solve_master(): """求解主问题""" # 这里需要根据具体模型实现 pass def add_cut_to_master(loss, scenario, x): """添加割平面""" pass

C&CG 的收敛速度取决于割平面的质量。我一般会在前几轮迭代时用“最坏场景”加割，后面用“多个次坏场景”加割，加速收敛。如果迭代 20 次还不收敛，通常是子问题求解不精确或割平面形式有误。

4.3 参数敏感性分析怎么做

ε（模糊集半径）、α（CVaR 置信水平）、N（场景数）是三个最关键的参数。敏感性分析的做法是固定两个、变化一个，画目标函数和失负荷概率的曲线。

ε 的拐点通常在 0.5 到 1.0 之间，超过 1.0 后投资成本急剧上升但失负荷改善有限。α 一般取 0.05 到 0.1，太小会导致模型过度保守。N 取 50 到 100 就足够稳定，再多边际收益很低。

5. 避坑与排查：复现时最容易翻车的五个地方

5.1 对偶推导符号搞反导致无界

现象：求解对偶问题时，linprog 返回“问题无界”或目标值为负无穷。

原因：Wasserstein 对偶的约束方向写反了。原始约束是 s_i ≥ g_j - λ d(i,j)，转成 linprog 的 ≤ 形式时，系数和常数都要取负。如果只取了一边，约束方向就错了。

解决：检查 A_ub 和 b_ub 的符号。建议先用一个小规模问题（N=3）手动验证对偶目标值和原始问题一致，再扩展到大规模。

5.2 距离矩阵量纲不统一导致 ε 失效

现象：模糊集半径 ε 调到很大，但模型仍然不保守，失负荷概率很高。

原因：距离矩阵各维度的量纲差异太大。比如路径坐标范围是 [-100, 100]，风速范围是 [10, 60]，欧氏距离被路径坐标主导，风速的变化被淹没。

解决：在计算距离矩阵前，先对场景做归一化。常见做法是 z-score 标准化或 min-max 归一化。归一化后，各维度对距离的贡献均衡，ε 的物理意义更清晰。

5.3 CVaR 近似机会约束的保守性被忽略

现象：模型求解成功，但实际失负荷概率远高于 1-α。

原因：CVaR 近似是保守的，但保守方向是“CVaR ≤ L_max 推出 P(X ≤ L_max) ≥ 1-α”，不是反过来。如果 L_max 设得太大，CVaR 约束虽然满足，但实际概率可能不达标。

解决：求解后要做样本外测试。用独立的测试场景集评估实际失负荷概率，如果低于 1-α，说明 L_max 需要收紧。我一般会把 L_max 设为理论值的 0.8 倍，留出保守裕度。

5.4 场景数太少导致模糊集退化

现象：N 小于 20 时，对偶问题的解不稳定，每次运行结果差异很大。

原因：Wasserstein 模糊集的经验分布 P̂_N 在 N 很小时，与真实分布偏差大。对偶问题的最优 λ 对场景数敏感，N 越小 λ 波动越大。

解决：N 至少取 50，最好 100。如果计算资源有限，可以用场景削减技术，先用 LHS 生成 200 个场景，再用 k-means 或快速前向选择削减到 50 个代表性场景。

5.5 拓扑约束和功率平衡的耦合导致不可行

现象：灾中阶段的优化问题频繁返回“不可行”。

原因：拓扑辐射状约束和功率平衡约束耦合太紧。如果某个孤岛内的负荷大于分布式电源出力，功率平衡无法满足，但拓扑约束又要求必须辐射状供电。

解决：在灾中阶段引入“切负荷”变量作为松弛，允许部分负荷被切除。切负荷的成本要设得足够高，让模型优先满足功率平衡，但在必要时可以切负荷。这样既保证了可行性，又不会让切负荷成为常态。

6. 进阶技巧：用样本外测试验证鲁棒性

复现论文最容易被忽略的一步是样本外测试。训练集上的表现好不代表模型真的鲁棒。我一般会留出 20% 的历史场景作为测试集，或者用独立的蒙特卡洛模拟生成测试场景。

具体做法是：固定规划决策 x*，在测试集上评估三个指标：平均失负荷、CVaR_α 失负荷、最大失负荷。如果测试集上的 CVaR 和训练集上的 CVaR 差距超过 15%，说明模糊集半径 ε 设得太小，模型过拟合了训练集。

def out_of_sample_test(x_star, test_scenarios, alpha=0.05): """ 样本外测试: 评估规划方案的鲁棒性 x_star: 规划决策 test_scenarios: 测试场景集 """ losses = [] for scenario in test_scenarios: loss = evaluate_loss(x_star, scenario) losses.append(loss) losses = np.array(losses) mean_loss = np.mean(losses) # CVaR 计算 sorted_losses = np.sort(losses)[::-1] n_tail = max(1, int(len(losses) * alpha)) cvar_loss = np.mean(sorted_losses[:n_tail]) max_loss = np.max(losses) print(f"平均失负荷: {mean_loss:.4f}") print(f"CVaR_{alpha} 失负荷: {cvar_loss:.4f}") print(f"最大失负荷: {max_loss:.4f}") return mean_loss, cvar_loss, max_loss # 假设已有规划决策和测试场景 # x_star = ccg_algorithm()[0] # test_scenarios = generate_typhoon_scenarios(100, 4, seed=123) # out_of_sample_test(x_star, test_scenarios)

另一个技巧是“压力测试”：人为构造比训练集更极端的场景，比如台风路径直接穿过变电站、风速超过历史最大值 20%。如果模型在这些场景下仍然能保持失负荷低于阈值，说明鲁棒性足够。如果不行，需要增大 ε 或增加资源容量。

我自己的习惯是：每次调完参数，先跑样本外测试，再看压力测试，最后才看训练集上的目标函数值。训练集上的数字好看没用，样本外不翻车才是真的稳。这套流程帮我省了很多后悔药——有几次训练集上目标函数降了 10%，但样本外 CVaR 涨了 30%，果断回退参数。

希望帮到你。

本文还有配套的精品资源，点击获取