分布式鲁棒电力系统状态估计PSSE复现：WLAV与ADMM的关键工程细节

最近在整理一套状态估计的复现代码时，偶然翻到早前做的一个项目——“分布式鲁棒电力系统状态估计源代码：高水平复现的PSSE方法及其实验验证”。这套代码当时是给一个区域电网的量测质量分析任务做的，前后折腾了大概三周。回头复盘，最值钱的不是“跑通了”，而是把鲁棒损失函数、区域分解、ADMM协调这三件事粘在一起的工程细节。这篇博文就围绕这套复现项目，把PSSE（Power System State Estimation，电力系统状态估计）是什么、为什么需要鲁棒、分布式怎么拆、源代码怎么组织、实验怎么验证，完整捋一遍。

1. 常规状态估计算法在真实量测面前为什么不够用

1.1 WLS方法的基本原理与数学框架

电力系统状态估计的核心任务，是利用SCADA/PMU上送的冗余量测，估计出系统各处节点的电压幅值和相角。这套机制是能量管理系统（EMS）的基石，后续的潮流计算、安全分析、经济调度都建立在它输出的“最可信系统断面”之上。

传统实现以加权最小二乘（WLS）为主。设系统状态向量为
x = [V_1, ..., V_n, θ_1, ..., θ_n]
其中V是节点电压幅值，θ是节点相角（参考节点θ=0）。量测向量z包含注入有功/无功、支路有功/无功、电压幅值等，量测方程为：

z = h(x) + e

e是服从均值为0、协方差矩阵R的高斯噪声。WLS的目标函数是：

min J(x) = [z - h(x)]^T W [z - h(x)]

其中W = R^(-1)，是量测权重矩阵。由于h(x)是非线性潮流方程，标准解法是Gauss-Newton迭代：在每步线性化得到雅可比矩阵H，然后解正规方程

Δx = (H^T W H)^(-1) H^T W [z - h(x)]

这个框架的理论性质非常漂亮：当量测噪声严格服从高斯分布且没有异常值时，WLS估计是最大似然估计，统计效率极高，方差最小。我在复核代码时用IEEE 118节点系统做基准测试，理想噪声环境下电压幅值平均估计误差能到0.0009 p.u.左右，这个精度对调度侧足够用了。

1.2 坏数据注入后WLS估计结果为何会系统性偏移

问题出在“没有异常值”这个前提上。真实量测系统里，坏数据不是小概率事件：CT/PT断线、通信通道误码、变电站在遥信变位瞬间的数据冻结、量测终端软件初始化阶段吐出的半个断面……这些都会产生远离真实值的局部极大偏差。

WLS对坏数据的反应是灾难性的。原因在于L2范数对残差取平方，残差被放大后，优化器会拼命去压低最大的那几个残差，导致整个状态向量向坏数据方向倾斜。坏数据的“污染”还会通过雅可比矩阵的耦合行传播到相邻节点——一个间隔的有功量测错了，附近几个节点的电压相角估计全部跟着偏移。

更隐蔽的是标准化残差检测（LNR检验）自身的缺陷。LNR检验基于估计残差判断坏数据，但多个坏数据同时存在时会发生masking效应：两个方向一致的坏数据互相掩饰，真实大的残差在标准化之后反而不超阈值；而好的量测可能被错误标记为坏数据（swamping效应），这就是“剔除坏数据”流程有时候越剔越乱的原因。

用一句直白的话说：真实调度环境里，你很难指望纯WLS在坏数据比例超过3%~5%时还能给出可信断面。这片土壤上，鲁棒状态估计的适用空间非常大。

2. 分布式鲁棒PSSE的整体方案与算法选型

2.1 为什么要做分布式：区域自治、数据隔离与故障隔离

既然要做鲁棒化，为什么还要叠加“分布式”这一层？我在项目里实际面对的需求是：电网按调度管辖范围分为多个区域，各区域只愿共享边界联络线的协调信息，绝不愿意把本地全部量测数据上送到一个中心节点。

这个约束在跨省跨网场景下非常真实。即便技术上中心式估计算力完全够用，但数据归属、隐私和调度权限决定了你必须把问题拆开。分布式还有一层意外收获：坏数据的传播范围被物理隔离在所属区域内部，一个区域的量测大面积异常时，相邻区域受到的冲击大幅减弱，这正好和鲁棒估计形成互补。

区域分解的做法很简单：按电气联系把全网节点划分成若干子区域，每个区域保留内部节点和边界节点。相邻区域共享边界节点电压，通过一致性约束连接——这保证了全网解的拓扑一致性。

2.2 算法选型比较：为什么是ADMM而不是分散式梯度或者辅助问题原理

分布式优化算法候选方案不少：

• 分散式梯度下降（distributed gradient descent）：实现简单，但对步长敏感，收敛慢，且需要每轮交换完整梯度向量。
• 辅助问题原理（APP）：收敛较慢，工程参数多。
• 交替方向乘子法（ADMM）：把目标函数拆分到各子区域独立求解，协调层只需要处理边界变量；有经典收敛性理论支撑；惩罚参数调好后收敛速度和稳定性都不错。

我在代码里最终选ADMM作为协调框架。ADMM有一个很实在的优势：每个区域内部子问题就像一个“小状态估计器”，你可以自由选择WLS、WLAV或者Huber等任意损失函数作为子问题的目标函数，外层算法结构不需要改。换句话说，鲁棒化改造被限制在子问题内部，协调层完全无感知，这为后续替换算法留下了极大灵活性。

2.3 鲁棒损失函数的选择：WLAV（加权最小绝对值）路线

鲁棒估计的常用路线有两种：一是对WLS迭代过程中的残差做统计检验再剔除坏数据（残差过滤法），二是直接改用对重尾噪声不敏感的损失函数（鲁棒估计法）。

残差过滤法实现简单，但在坏数据比例高、存在masking效应时效果不稳定。我采用的是WLAV（Weighted Least Absolute Value），把目标函数改为：

min Σ w_i |z_i - h_i(x)|

从统计学角度，L1回归对响应变量中的异常值有天然的抵抗能力；从计算角度，WLAV问题可以转换成线性规划（LP），每一轮迭代都有成熟的求解器可以接。它在IEEE标准算例下对10%~20%坏数据的抑制效果非常直观：坏数据量测对应的残差不参与主导目标，状态估计值不会被拉偏。

3. 核心数学模型与迭代公式推导

3.1 分布式鲁棒状态估计的统一数学形式

设全网分为 R 个区域。对区域 r，其状态向量 x_r 由内部节点和边界节点两部分组成，量测方程为：

z_r = h_r(x_r) + e_r

相邻区域之间对共享边界节点的状态一致性约束为：

x_r^bd - x_s^bd = 0

把 WLAV 和一致性约束写进同一个优化问题，得到：

min Σ_r ||W_r (z_r - h_r(x_r))||_1

s.t. x_r^bd - x_s^bd = 0, ∀(r,s) ∈ 邻接关系

这个形式的特点是：目标函数按区域可分，约束只作用在边界变量上，非常适合用 ADMM 做分解。

3.2 WLAV子问题的线性规划转换

对单个区域的 WLAV 子问题，核心在于把绝对值项线性化。引入非负松弛变量 u_i, v_i，令：

z_i - h_i(x) = u_i - v_i

则 |z_i - h_i(x)| = u_i + v_i。子问题变成：

min Σ_i w_i (u_i + v_i)

s.t. z_i - h_i(x) - u_i + v_i = 0

u_i ≥ 0, v_i ≥ 0

由于 h_i(x) 非线性，外循环采用逐次线性化：当前迭代点 x^k 附近把 h_i(x) 展开为 H_i Δx + h_i(x^k)，则每次迭代求解一个 LP 或 QP 获得增量 Δx，再更新状态。实测情况下，这种序列线性规划方法在状态估计问题上收敛性很好，一般 5~8 轮外循环即可达到 10^-5 的精度。

3.3 ADMM外层迭代的三个核心步骤

引入全局辅助变量 y 表示边界节点的“协调值”，把一致性约束等价写成 x_r^bd - y = 0。对区域 r，增广拉格朗日函数为：

L_ρ = ||W_r(z_r - h_r(x_r))||_1 + λ_r^T(x_r^bd - y) + (ρ/2) ||x_r^bd - y||_2^2

ADMM 在每次迭代内执行三步：

第一步（本地子问题求解）：固定协调变量 y^k 和对偶变量 λ^k，各区域并行求解：

x_r^{k+1} = argmin_{x_r} ||W_r(z_r - h_r(x_r))||_1 + (ρ/2) ||x_r^bd - y^k + u_r^k||_2^2

其中 u_r^k = λ_r^k / ρ。这一步每个区域独立执行，互不干扰，是天然的可并行部分。

第二步（协调变量更新）：收集所有区域的边界状态，取各区域边界状态与对偶项的平均值：

y^{k+1} = (1/N) Σ_r (x_r^{bd,k+1} + u_r^k)

这里 N 是区域内边界节点的数量。平均值规则产生的 y 会作为下一次迭代各区域子问题的基准。

第三步（对偶变量更新）：更新每个区域的对偶项：

u_r^{k+1} = u_r^k + ρ (x_r^{bd,k+1} - y^{k+1})

对偶更新的物理含义是：如果某个区域边界状态与全局协调值存在偏差，就对它施加矫正性惩罚，促使最终收敛时边界变量全网一致。

3.4 收敛判据与关键参数

ADMM 迭代是否该停下来，我同时监控两种残差：原始残差 s^k = Σ_r ||x_r^{bd,k} - y^k||_2，以及对偶残差 t^k = ρ ||y^k - y^{k-1}||_2。当两者同时小于阈值（我通常设 10^-4）时判定收敛。

惩罚参数 ρ 是唯一最需要调节的参数。ρ 太小，协调变量更新导致的状态变化幅度大，收敛很慢；ρ 太大，子问题里二次惩罚项权重过高，目标函数中 WLAV 部分的影响力被稀释，鲁棒性反而变差。在 IEEE 118 节点系统上，我取 ρ=10 时收敛速度和鲁棒表现都处于良好区间。

4. 源代码工程化实现的几个关键细节

4.1 代码仓库结构：关注点分离

整套复现代码我按模块组织，核心思路是“数据、算法、实验独立演进”：

• data/：存放测试系统（IEEE CDF/MATPOWER格式）和量测数据生成脚本
• core/：网络拓扑解析、量测映射、雅可比矩阵构造
• algorithms/：wls_core.py（集中式参考实现）、wlav_admm.py（分布式鲁棒估计主模块）
• experiments/：run_experiments.py、结果对比与指标计算
• utils/：收敛判据、误差指标、稀疏矩阵工具

algorithms/wlav_admm.py 里最核心的伪代码结构如下：

python复制def run_admm_wlav(network, measurements, regions, rho=10.0, max_iter=50):
    # 初始化：各区域状态变量取平启动
    x = {r: flat_start(network, regions[r]) for r in regions}
    # 边界节点全局协调值
    y = init_boundary_coordination(network, regions)
    # 各区域对偶变量
    u = {r: np.zeros_like(y[r]) for r in regions}

    for k in range(max_iter):
        # 第一步：各区域并行求解 WLAV 子问题
        for r in regions:
            x[r] = solve_wlav_subproblem(
                network, measurements[r], regions[r],
                y, u[r], rho
            )
        # 第二步：更新协调变量
        y_new = update_coordination(x, u)
        # 第三步：更新对偶变量
        for r in regions:
            u[r] = u[r] + rho * (x[r]["bd"] - y_new)
        # 收敛检查
        primal_res = compute_primal_residual(x, y)
        dual_res = rho * np.linalg.norm(y_new - y)
        if primal_res < 1e-4 and dual_res < 1e-4:
            break
        y = y_new
    return assemble_global_state(x)

这段代码最大的好处是：区域循环是完全解耦的，换成多进程或者分布式消息队列只需要在这个循环上做并行化改造，算法逻辑一行不用动。

4.2 雅可比矩阵构造：稀疏模式与缓存

状态估计里最容易被性能卡住的就是雅可比矩阵 H。在 WLS 迭代中，每次状态更新都需要重新计算 H，而 H 的结构由量测方程决定：注入功率量测只关联本节点及相邻节点状态，支路功率量测只关联两端节点状态。因此 H 的稀疏模式与节点导纳矩阵几乎一致。

我在代码里用 scipy.sparse 的 CSR 格式存储 H，并且在一开始就把稀疏模式固定下来，迭代中只更新非零元数值，不重新分配内存。这样在 118 节点系统、约 400 个量测的情况下，单次雅可比计算和正规方程求解耗时在 30ms 以内，整轮 WLS 5 次迭代总耗时不到 0.2 秒。

4.3 求解器的选择：从CVXPY原型到定制求解

做算法原型阶段，我直接用 CVXPY 表达 WLAV 子问题，一行建模、底层自动调求解器，非常省事。但它的缺点是每次子问题求解都有建模和编译开销，在 ADMM 需要迭代 20~30 轮的应用场景下显得笨重。

后来我把子问题改成基于内点法思想的定制求解器：利用子问题稀疏的约束结构，对 KKT 系统做稀疏 Cholesky 分解，迭代过程中复用分解模式。这一步重构之后，单轮 ADMM 整体耗时从 2 秒左右降到 0.4 秒。如果你复现时对性能没有极端要求，CVXPY 版本完全够用，它更易读、更易改。

4.4 澄清：这里的“分布式”是算法分解，并非大数据框架

搜索相关关键词时，不少开发者把“分布式状态估计”和“Hadoop 分布式计算”混为一谈。需要明确：PSSE 分布式化的核心是把网络按拓扑切块，让各区域独立解子问题再协调边界，它并不需要 HDFS，也不需要 Zookeeper 或 Spark 集群。代码可以在单机上跑，区域之间的“通信”通过函数调用或者本地消息传递完成。

不过，如果你要处理的是数千节点、量测上万的超大规模系统，ADMM 的区域循环天然适合并行化——用 multiprocessing、Ray 或者其他并发框架把区域子问题分发到多核甚至多机，完全能无缝迁移。这也是我选择 ADMM 的一个重要考量。

5. 实验验证：从IEEE标准算例到坏数据场景

5.1 测试系统与实验配置

复现实验选用 IEEE 118 节点系统，这是状态估计领域最常用的中等规模标准算例。系统被切分为 3 个区域，区域间通过联络线连接。量测冗余度控制在 3.0 左右，量测噪声设置为：功率量测标准差 0.02 p.u.，电压幅值量测标准差 0.005 p.u.。

实验分为两组场景：第一组是无坏数据的纯净量测，第二组在量测数据中加入坏数据。坏数据生成方式：在每个区域内随机抽取 20% 的功率量测，注入 3~8 倍标准差的偏移，模拟传感器偏差或通信误码场景。

5.2 正常工况下的对比结果

先看无坏数据场景下的基准测试结果，以电压幅值平均绝对误差（MAE）、相角 MAE 和迭代次数为指标：

可以看到，在理想噪声环境下，鲁棒方法付出了少量精度代价：电压幅值误差略高约 0.0002 p.u.，相角误差略高约 0.0002 rad。这个代价是 L1 损失函数在高斯噪声下统计效率略低于 L2 的必然结果，工程上完全可以接受。

5.3 坏数据场景下的有效性验证

加入 20% 坏数据后，结果分化非常明显：

集中式 WLS 在坏数据场景下电压幅值误差放大到了原来的 5 倍以上，相角误差甚至超过 0.006 rad，整个估计断面已经偏离真实运行工况。加上坏数据剔除流程后虽然有所改善，但受 masking 效应影响，仍无法恢复到理想噪声水平。分布式 WLAV+ADMM 则保持与纯净场景几乎一致的误差水平，证明 20% 坏数据比例下鲁棒状态估计能有效抵抗污染。

5.4 分布式求解效率与收敛性分析

从收敛曲线来看，ADMM 原始残差在前 5 轮迭代内快速下降约两个数量级，之后进入线性收敛区间。实际迭代 15~18 轮即可满足 10^-4 的收敛阈值，总耗时约 0.4 秒。相比集中式 WLS 多出约 0.2 秒，但这在 EMS 的断面刷新周期（通常不小于 1 秒）内完全可接受。

真正值得注意的是：ADMM 各区域子问题在收敛过程中呈现出一定程度的“独立波动”，一个区域内部坏数据导致该区域状态在小范围内抖动，但通过边界协调变量，这种抖动不会传导到其他区域。这正好验证了分布式架构对坏数据传播的天然隔离效果。

6. 复现过程中踩过的坑与调参心得

6.1 节点编号陷阱：从1-indexed到0-indexed的转换

IEEE 标准测试系统的节点编号通常从 1 开始，但 Python 环境以及 scipy 稀疏矩阵使用 0-indexed。这个看似简单的转换，实际调试中花费了大量时间：雅可比矩阵的行列索引错位不会立刻报错，而是在迭代两三轮之后产生异常跳变，非常难定位。

建议在读取测试系统文件时，统一维护一个 bus_id_mapping 字典，显式把原始编号映射到连续整数索引，后续所有计算都基于映射后的索引，避免在代码各处零零散散做加减一的转换。

6.2 坏数据比例的临界点：鲁棒方法也会失效

你可能会以为 WLAV 是“万能盾牌”，能抵抗任意比例的坏数据。实验表明，当坏数据比例超过 35%~40% 时，L1 目标函数的鲁棒性急剧恶化。原因是坏数据占比过高后，目标函数中“正常量测”的总权重不再占优，优化器会倾向于拟合坏数据簇，导致状态估计整体偏离。

这一点在真实系统里意义重大：如果发现某个区域超过三分之一量测同时异常，不要指望鲁棒估计器帮你扛，此时正确的做法是直接隔离该区域，先查量测系统故障源。

6.3 惩罚参数 ρ 的调节经验

ADMM 对 ρ 的敏感性是复现中最大的调参痛。根据经验，ρ 对收敛行为的影响有规律：

• ρ 过小（例如 0.1）：原始残差下降极慢，迭代 50 轮后仍未收敛，且各区域边界状态长期不一致。
• ρ 过大（例如 1000）：子问题中二次惩罚项占据主导，WLAV 丢失鲁棒特性，坏数据场景下误差指标反弹明显。
• ρ 在 5~20 区间表现稳定。

我最后的做法是：先用纯净量测场景跑一个 ρ 扫描，画出迭代轮数-误差曲线，选一个曲线平坦区域的中间值，然后固定 ρ 再做坏数据实验。实际项目中，ρ 不需要每个场景都重新调，一套测试系统选定后可以长期复用。

6.4 初始化选择对非线性迭代的影响

WLAV 子问题采用逐次线性化求解，对初始点有要求。我实验中发现，如果取经典“平启动”（所有电压幅值=1.0，所有相角=0.0），在坏数据场景下个别区域会出现小幅振荡，收敛轮数增加。后来改成先用 WLS 做一个一步预估计，把 WLS 的解作为 ADMM 的初始状态，振荡问题基本消失，收敛速度也有提升。

这个技巧成本极低，收益却非常明显：预估计只需要一次全局或区域级 WLS 迭代，却能显著改善后续 L1 优化的线性化起点质量。

6.5 对偶更新的缩放：保持数值尺度一致性

电力系统状态变量有两个不同的量纲：电压幅值在 1.0 p.u. 附近，相角在零点零几弧度量级。如果对电压和相角直接施加相同的 ADMM 更新步长，收敛过程中边界相角的一致性约束会在数值尺度上“碾压”电压幅值变量，导致收敛失衡。

我的处理方式是在一致性约束中给相角乘上一个尺度因子 α（例如 10），让两个状态分量在 ADMM 的罚项中具有相近的数值阶数。这个操作完全等价于改变罚项单位的度量选择，但能让收敛行为明显更平滑。

复现这套代码后的一些体会

把 PSSE 从集中式 WLS 升级到分布式鲁棒版本，并不是一个单纯的算法替换过程。最花精力的其实是三件事：一是把鲁棒损失函数正确地嵌套进 ADMM 子问题结构，二是让雅可比矩阵和稀疏求解器在迭代循环中发挥出真实性能，三是理解坏数据在不同场景下的“攻防博弈”——鲁棒估计防御的是常规量测污染，数据质量大面积崩塌时更应该做的是设备级隔离而不是算法修正。

这套代码最终在项目里沉淀为两条使用路径：对正常调度断面，集中式 WLS 快速出数；对事故后或者数据质量存疑场景，切到分布式鲁棒估计器做交叉验证。两条路径共用一套拓扑解析和量测映射模块，维护成本可控，实验对比也方便。希望这篇复盘对你复现类似系统时少踩几个坑。