基于预处理共轭梯度迭代法的电力系统状态估计算法

李建斌，王鹏程，傅 侃，方 睿，董树锋

（1. 国网浙江杭州市萧山区供电有限公司，浙江省杭州市 311200；2. 杭州电力设备制造有限公司萧山欣美成套电气制造分公司，浙江省杭州市 311200；3. 浙江大学电气工程学院，浙江省杭州市 310027）

0 引言

电力系统状态估计是现代能量管理系统（energy management system，EMS）的基础，为EMS中的高级应用提供底层支撑。随着中国电力系统省地一体化［1］和输配一体化［2］的发展，电力系统计算维度越来越高，计算时间急剧增加，传统状态估计算法难以满足快速增长的计算需求。因此，亟须研究更加快速准确的电力系统状态估计算法。

目前，电力系统中最广泛使用的状态估计算法是加权最小二乘（weighted least squares，WLS）法［3］。这种方法假设量测量服从正态分布，数学模型较为简单、迭代次数少、计算速度快，在测点集中、无不良数据时估计效果较好。WLS 状态估计方法以牛顿-拉夫逊法为基础，通过逐次线性化、多次求解线性方程组的方式，反复迭代得到状态变量的最终解。在WLS 状态估计方法的求解过程中，需要耗费大量的时间求解高维稀疏矩阵乘法和高维稀疏线性方程组，占据了绝大部分的计算时间。同时，由于WLS 状态估计方法抗差性较差，有学者也提出了一些抗差状态估计方法［4-5］，能够减少不良数据对状态估计结果的影响。

求解大规模稀疏线性方程组主要可以分为直接法和迭代法。直接法是利用矩阵分解和变换技术对原线性方程组进行消去，代表方法有高斯消去法［6］、LU 分解法［7］等。这类方法的特点是对主元的选取较为苛刻，对于病态的稀疏矩阵需要通过置换行列的方式得到合适的对角主元。同时，直接法占用的内存较多，难以并行计算，在方程规模达到一定数量级后求解效率较低，不适合计算大型稀疏线性方程组。文献［8-10］使用直接法求解电力系统潮流和状态估计问题；文献［11-12］通过图形处理器（graphics processing unit，GPU）对状态估计算法进行并行加速，在稀疏线性方程组求解方面采用并行LU 分解算法，但并行的LU 分解法加速效果不明显，不适合运用在大规模系统中。

相比于直接法，迭代法对于大型稀疏线性方程组的计算具有较大的优势。迭代法每次计算时的内存占用低，且非常适合并行，当矩阵规模增大时迭代法的计算效率不会受到影响。但其主要问题在于收敛性，受矩阵的条件数影响较大。目前，针对这一问题，主要采用预处理技术降低条件数，代表方法有不完全LU 分解法［13］、不完全Cholesky 分解法等［14］。近年来，由于电力系统规模的不断增大，已有不少研究将迭代法运用在电力系统分析和计算领域。文献［15-17］对潮流计算中的雅可比矩阵进行预处理，并利用迭代法并行计算线性方程组，解决了大规模潮流计算的效率问题。文献［18］提出了一种基于不完全LU 分解预处理迭代法的潮流计算法，但该方法主要针对的是潮流计算问题而没有能够扩展到状态估计问题。

本文基于预处理共轭梯度迭代法，针对WLS 状态估计的特点，提出一种快速的电力系统状态估计算法。该算法在预处理后矩阵条件数得到明显下降，同时基于GPU 并行计算技术实现矩阵乘法和迭代法的高效并行。算例分析验证了本文算法的快速性、收敛性和低内存占用等特点。

1 状态估计基本原理

1.1 WLS 状态估计

电力系统的非线性量测状态估计方程可以表示为：

式中：z为系统量测向量；x为系统状态变量；h（x）为在状态x下的量测函数；v为量测误差。

定义残差矢量为：

求解状态估计问题变为求解潮流问题的扩展，可以建立优化模型：

式中：R为量测误差方差矩阵，表示各量测量的准确程度。

为了使式（3）的目标函数值最小，则有：

采用牛顿法求解，可以得到迭代修正量［19］为：

1.2 并行性分析

GPU 的统一计算设备架构（compute unified device architecture，CUDA）［20］具备并行计算功能，非常适合处理矩阵和向量的并行运算。

从式（5）可以看出，WLS 状态估计中耗时部分主要是矩阵乘法与线性方程组求解这2 个步骤。在牛顿法的每一次迭代过程中都需要求解矩阵乘法和线性方程组，其中涉及大量运算，占据了耗时的绝大部分。

矩阵乘法步骤包括矩阵A和向量b的运算，其计算公式如下。

由于矩阵H(x̂k)为高维稀疏矩阵，因此经过矩阵乘法后矩阵A仍是高维稀疏矩阵。在向量b中，由于向量z−h(x̂k)为稠密向量，因此b也是高维稠密向量。 此步骤采用CUDA 中的运算库cuSPARSE［21］即可并行高效计算出矩阵A和向量b的结果，根据CUDA 的描述，计算速度比纯CPU 替代产品快2 至5 倍。

矩阵线性方程组计算部分最为耗时。从式（9）可以看出，矩阵A是矩阵转置、误差权重矩阵和矩阵本身之间的乘积，是对称正定矩阵。当系统规模较小时，矩阵规模较小，适合使用LU 分解等直接法对线性方程组进行计算。常用的高性能线性方程组求解库SuperLU［22］经过了高度优化，属于直接法，效率较高。但对于大规模稀疏线性方程组，矩阵条件数较高，直接法的计算效率难以满足要求，适合在GPU 上使用迭代法进行并行计算。

2 迭代法基本原理

2.1 Krylov 子空间法

Krylov 子空间法是求解大规模线性方程组的首选方法。求解一般线性方程组：

投影方法的基本思想是从一个维数更小的子空间Km内寻找近似的解。这个子空间Km被称为搜索空间，其维度为m。

此时，设置m个约束，要求残差矢量r满足m个正交条件，即Petrov-Galerkin 条件：

式中：Lm为另一个m维的子空间，被称为约束空间。当Lm=Km时，该方法为正交投影法，否则为斜投影法。

给定迭代初值x0，采用仿射空间x0+Km可以得到：

式中：x͂为迭代值，初始残差r0=b−Ax0。

在Krylov 子空间方法中，搜索空间Km就是Krylov 子空间，定义为：

式中：r可选为初始残差r0，Krylov 子空间方法就是在Krylov 子空间中寻找近似解。

当选用不同的约束空间Lm时，对于迭代过程会有比较大的影响。主要可以分为如下方法。

1）正交投影法：Lm=Km(A，r0)，代表方法有共轭梯度法。该方法计算步骤较为简单，但要求矩阵A对称正定。

2）正交化方法：Lm=AKm(A，r0)，代表方法有广义极小残差方法。该方法数值稳定性高，但计算量会随迭代的步数线性增长。

3）双正交化方法：Lm=AKm(AT，r0)，代表方法有稳定的双正交共轭梯度法。该方法适用性最广且稳定性较高，但计算步骤较为复杂，计算量比前2 种方法要大。

搜索空间Km的选取不会对最终状态估计结果的精度有影响，其选取依据主要是状态估计中矩阵A的特性，例如是否对称、是否正定等。在状态估计的系数矩阵A对称正定情况下，可以使用正交投影法确定约束空间Lm，此方法计算量较小、步骤较简单。本文采用共轭梯度法进行求解。从式（9）可以看出，矩阵A为对称正定矩阵，满足共轭梯度法的要求，能够充分发挥其计算步骤简单、效率高的优势。

2.2 预处理方法

如果直接使用上述Krylov 子空间方法进行迭代，原始矩阵条件数过高，可能会出现收敛性差、迭代次数多等问题。采用合适的预处理方法可以降低矩阵条件数。预处理方法的原理是通过预处理子M将原始线性方程组转化为另一个易于求解的线性方程组。

常见的预处理方法可分为如下几类［23］。

1）左预处理：对M−1Ax=M−1b使用迭代法。

2）右预处理：对AM−1u=b使用迭代法，其中x=M−1u。

3）分裂预处理：对M1−1u=M1−1b使用迭代法，预处理子M=M1M2。

对于大规模电力系统，其矩阵条件数较高，利用预处理技术可以减少迭代次数，易于问题求解。不完全LU 分解预处理方法适用范围广，且在文献［18］中已经被证明在潮流计算中有较好的预处理效果。因此，本文选用不完全LU 预处理方法。

2.3 预处理共轭梯度迭代法

考虑对称正定矩阵A，并假设有一个预处理子M可供使用，且M是对称正定的。为了保持对称性，注意到M−1A对M内积(x，y)M≡(Mx，y)=(x，My)是自伴的，因为

式中：x和y为任意向量。

所以，一种选择是用M内积代替共轭梯度迭代法中的欧氏内积，如果按这种新内积重写共轭梯度迭代法，用rj=b−Axj表示原残差，而用zj=M−1rj表示预处理方程组的残差，可以得到本文所使用的预处理共轭梯度迭代法，且不需要显式地计算M内积。预处理共轭梯度迭代法的具体步骤见3.1 节。

3 GPU 并行加速状态估计算法

3.1 算法步骤

本文提出一种基于预处理共轭梯度迭代法的电力系统WLS 状态估计算法，其计算步骤如下。

步骤1：初始化，形成节点导纳矩阵，给状态变量赋初值形成x̂0。

步骤2：设置迭代变量k=0 和最大迭代次数kmax。

步骤3：根据当前状态变量，计算雅可比矩阵H()，计算公式如文献［19］附录中所示。

步骤4：在GPU 上利用cuSparse 库，根据式（9）计算矩阵A和向量b。

步骤5：在GPU 上求解线性方程组Ax=b，利用矩阵A对称正定的特点，使用预处理共轭梯度法进行迭代求解［24］，具体做法如下。

①对矩阵A进行ILU（0）分解，ILU（0）分解是不完全LU 分解的一种形式，形成矩阵A的预处理子：

式中：L和U分别为ILU（0）分解出的上三角矩阵和下三角矩阵。设置迭代次数i=0 和最大迭代次数imax，设x的初始猜测为x0，计算初始残差r0及其2 范数‖r0‖。

②根据L和U求解方程组Mz=ri，其中ri为迭代次数为i时的残差，得到z后计算ρi=(ri，z)。

③如果i=0，则令pi=z，否则令β=ρi/ρi−1，并计算pi=z+βpi−1。

步骤6：根据步骤5 中解出的x即可得到牛顿法迭代的修正量Δx̂k，根据式（7）得到状态变量x̂k+1。

步骤7：令k=k+1，若不满足式（8）且k没有达到kmax则转至步骤3，否则退出状态估计过程。

3.2 算法特点

1）算法采用ILU（0）预处理方法，保证了残差矩阵满足ILU（0）的分解条件。这种预处理方法产生的预处理子不会注入非零元。经过预处理后，能够保证预处理子的稀疏性。在大型稀疏矩阵中，保证预处理子的稀疏性可以在矩阵运算中节省计算量和内存，同时利用预处理子与原矩阵相同的稀疏性能够更快地加速迭代法的计算。从迭代步骤可以看出，计算只涉及稀疏矩阵向量线性乘法，且由于没有更多的非零元注入，理论上不会产生内存泄漏。

2）算法采用了共轭梯度法，WLS 状态估计中线性方程组系数矩阵A对称正定的特性使共轭梯度法苛刻的适用条件得到满足。在大型稀疏线性方程组求解过程中，迭代法的计算效率更高，同时共轭梯度法作为迭代法中计算步骤最简单的方法，具有最少的计算量和最高的计算效率。

3）算法采用了GPU 并行计算架构，充分利用CUDA 的高性能矩阵向量计算技术和cuSparse 库的乘法运算，加速矩阵A和向量b的形成。同时，在共轭梯度法迭代过程中，使用GPU 快速计算各个中间变量，保证了迭代法的快速计算。

4 算例分析

为了分析本文所提算法的有效性，选取部分算例进行测试。其中，CASE300 是IEEE 标准测试算例，分别将10、20、50、100 个CASE300 算例以平衡节点为中心进行拼接，可以得到拼接算例CASE2991、 CASE5981、 CASE14951和CASE29901。本文的量测数据是在潮流计算的基础上加入2%的高斯噪声所形成的。算例测试环境为64 bit 的Windows 10 操作系统，CPU 型号为Intel Core i7-5820K，GPU 型号为NVIDIA GeForce GTX1080。其中，牛顿法的迭代精度设置为10−3，线性方程组求解的迭代精度设置为10−10。

4.1 计算耗时对比

为了分析本文算法的并行加速效果，将本文算法与直接法进行比较。其中，直接法使用高性能SuperLU 库，能够较快地使用列选主元高斯消去法求解高性能稀疏线性方程组，计算结果如表1所示。

表1 本文算法与直接法计算效果对比Table 1 Comparison of calculation results between method of this paper and direct method

从表1 可以看出，本文算法随着算例规模的增大，计算加速效果越来越明显。当算例规模较小时，如CASE300 算例中本文算法计算速度比直接法更慢，这是由于直接法在线性方程组维数较小时LU分解速度较快。但当系统规模增大时，本文算法受问题规模的影响较小，而此时直接法的计算时间会急剧增大。上述算例证明在大规模电力系统中，本文算法能够大幅提高计算速度，与直接法相比优势明显。

4.2 计算精度对比

为了分析本文算法的估计精度误差，将本文算法与直接法进行比较。通过计算本文算法和直接法估计结果与真值之间的平均相对误差来评估计算精度，计算结果如表2 所示。

表2 本文算法与直接法计算精度对比Table 1 Comparison of calculation precision between method of this paper and direct method

从表2 可以看出，使用本文算法状态估计结果的平均相对误差在2%以内，与直接法精度相比几乎没有差别。从计算步骤上来看，本文方法与直接法的不同之处在于求解线性方程的步骤，而由于线性方程组的迭代精度设置较为严苛，设为10−10，求得方程的解与直接法相差无几。算例结果表明本文算法不会对状态估计精度造成过大的影响。

4.3 预处理技术有效性分析

良好的预处理有助于减少迭代次数，改善迭代法的稳定性。为了验证本文所采用的ILU（0）预处理方法的有效性，对预处理前后矩阵的条件数和迭代次数进行分析，计算结果如表3 所示。其中，预处理条件数和迭代次数取各次迭代的平均值。

表3 ILU(0)预处理方法效果分析Table 3 Effect analysis of ILU(0)preprocessing method

从表3 可以看出，经过ILU（0）预处理后，预处理子的条件数大幅下降。在迭代法中，每次求解都需要计算预处理子所形成的线性方程组，条件数下降意味着更容易对该方程组进行求解。因此，预处理后共轭梯度法的迭代次数也随之下降，甚至可以改善原本不收敛的线性方程组的收敛性。

减少迭代次数能够缩短计算时间，且该预处理方法不注入非零元，不会因稀疏性被破坏而大幅增加计算时间。算例结果表明，本文所采用的预处理方法能够有效减少共轭梯度法的迭代次数，提高本文算法的计算效率。

4.4 内存占用测试

在大型电力系统计算中，内存占用是一个重要的性能指标。为了验证本文算法的内存占用情况，对内存占用和显存占用进行测试，计算结果如表4所示。

表4 本文算法内存占用测试Table 4 Memory footprint test of method in this paper

本文算法为了节约内存，矩阵存储方式为压缩矩阵形式。从表4 可以看出，即使对于节点数约3 万的大规模电力系统，本文算法内存占用不超过700 Mbit，显存占用不超过600 Mbit。整体而言，本文所用算法内存和显存占用峰值都较低，在普通的计算机上即可进行计算。

5 结语

针对大型电力系统状态估计速度较慢的问题，本文提出了一种基于预处理共轭梯度迭代法的电力系统状态估计算法。本文算法是在WLS 状态估计法的基础上，针对牛顿法求解过程中矩阵乘法和大规模稀疏线性方程组求解效率较低的特点，在GPU上进行并行加速。由于WLS 状态估计中线性方程组稀疏矩阵为对称正定矩阵，满足共轭梯度法的使用条件，本文提出采用不完全LU 分解预处理方法和共轭梯度迭代法进行线性方程组的求解。算例分析验证了本文算法计算速度快、内存占用低的特点，并验证了预处理方法的有效性。综上所述，本文所提算法能够满足大规模电力系统状态估计的实时性要求，未来可考虑多机多卡下的分布式计算，进一步提高状态估计的计算效率。