李建斌，王鵬程，傅 侃，方 睿，董樹鋒

（1. 國網(wǎng)浙江杭州市蕭山區(qū)供電有限公司，浙江省杭州市 311200；2. 杭州電力設(shè)備制造有限公司蕭山欣美成套電氣制造分公司，浙江省杭州市 311200；3. 浙江大學(xué)電氣工程學(xué)院，浙江省杭州市 310027）

0 引言

電力系統(tǒng)狀態(tài)估計是現(xiàn)代能量管理系統(tǒng)（energy management system，EMS）的基礎(chǔ)，為EMS中的高級應(yīng)用提供底層支撐。隨著中國電力系統(tǒng)省地一體化［1］和輸配一體化［2］的發(fā)展，電力系統(tǒng)計算維度越來越高，計算時間急劇增加，傳統(tǒng)狀態(tài)估計算法難以滿足快速增長的計算需求。因此，亟須研究更加快速準(zhǔn)確的電力系統(tǒng)狀態(tài)估計算法。

目前，電力系統(tǒng)中最廣泛使用的狀態(tài)估計算法是加權(quán)最小二乘（weighted least squares，WLS）法［3］。這種方法假設(shè)量測量服從正態(tài)分布，數(shù)學(xué)模型較為簡單、迭代次數(shù)少、計算速度快，在測點集中、無不良數(shù)據(jù)時估計效果較好。WLS 狀態(tài)估計方法以牛頓-拉夫遜法為基礎(chǔ)，通過逐次線性化、多次求解線性方程組的方式，反復(fù)迭代得到狀態(tài)變量的最終解。在WLS 狀態(tài)估計方法的求解過程中，需要耗費大量的時間求解高維稀疏矩陣乘法和高維稀疏線性方程組，占據(jù)了絕大部分的計算時間。同時，由于WLS 狀態(tài)估計方法抗差性較差，有學(xué)者也提出了一些抗差狀態(tài)估計方法［4-5］，能夠減少不良數(shù)據(jù)對狀態(tài)估計結(jié)果的影響。

求解大規(guī)模稀疏線性方程組主要可以分為直接法和迭代法。直接法是利用矩陣分解和變換技術(shù)對原線性方程組進行消去，代表方法有高斯消去法［6］、LU 分解法［7］等。這類方法的特點是對主元的選取較為苛刻，對于病態(tài)的稀疏矩陣需要通過置換行列的方式得到合適的對角主元。同時，直接法占用的內(nèi)存較多，難以并行計算，在方程規(guī)模達到一定數(shù)量級后求解效率較低，不適合計算大型稀疏線性方程組。文獻［8-10］使用直接法求解電力系統(tǒng)潮流和狀態(tài)估計問題；文獻［11-12］通過圖形處理器（graphics processing unit，GPU）對狀態(tài)估計算法進行并行加速，在稀疏線性方程組求解方面采用并行LU 分解算法，但并行的LU 分解法加速效果不明顯，不適合運用在大規(guī)模系統(tǒng)中。

相比于直接法，迭代法對于大型稀疏線性方程組的計算具有較大的優(yōu)勢。迭代法每次計算時的內(nèi)存占用低，且非常適合并行，當(dāng)矩陣規(guī)模增大時迭代法的計算效率不會受到影響。但其主要問題在于收斂性，受矩陣的條件數(shù)影響較大。目前，針對這一問題，主要采用預(yù)處理技術(shù)降低條件數(shù)，代表方法有不完全LU 分解法［13］、不完全Cholesky 分解法等［14］。近年來，由于電力系統(tǒng)規(guī)模的不斷增大，已有不少研究將迭代法運用在電力系統(tǒng)分析和計算領(lǐng)域。文獻［15-17］對潮流計算中的雅可比矩陣進行預(yù)處理，并利用迭代法并行計算線性方程組，解決了大規(guī)模潮流計算的效率問題。文獻［18］提出了一種基于不完全LU 分解預(yù)處理迭代法的潮流計算法，但該方法主要針對的是潮流計算問題而沒有能夠擴展到狀態(tài)估計問題。

本文基于預(yù)處理共軛梯度迭代法，針對WLS 狀態(tài)估計的特點，提出一種快速的電力系統(tǒng)狀態(tài)估計算法。該算法在預(yù)處理后矩陣條件數(shù)得到明顯下降，同時基于GPU 并行計算技術(shù)實現(xiàn)矩陣乘法和迭代法的高效并行。算例分析驗證了本文算法的快速性、收斂性和低內(nèi)存占用等特點。

1 狀態(tài)估計基本原理

1.1 WLS 狀態(tài)估計

電力系統(tǒng)的非線性量測狀態(tài)估計方程可以表示為：

式中：z為系統(tǒng)量測向量；x為系統(tǒng)狀態(tài)變量；h（x）為在狀態(tài)x下的量測函數(shù)；v為量測誤差。

定義殘差矢量為：

求解狀態(tài)估計問題變?yōu)榍蠼獬绷鲉栴}的擴展，可以建立優(yōu)化模型：

式中：R為量測誤差方差矩陣，表示各量測量的準(zhǔn)確程度。

為了使式（3）的目標(biāo)函數(shù)值最小，則有：

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

1.2 并行性分析

GPU 的統(tǒng)一計算設(shè)備架構(gòu)（compute unified device architecture，CUDA）［20］具備并行計算功能，非常適合處理矩陣和向量的并行運算。

從式（5）可以看出，WLS 狀態(tài)估計中耗時部分主要是矩陣乘法與線性方程組求解這2 個步驟。在牛頓法的每一次迭代過程中都需要求解矩陣乘法和線性方程組，其中涉及大量運算，占據(jù)了耗時的絕大部分。

矩陣乘法步驟包括矩陣A和向量b的運算，其計算公式如下。

由于矩陣H(x?k)為高維稀疏矩陣，因此經(jīng)過矩陣乘法后矩陣A仍是高維稀疏矩陣。在向量b中，由于向量z?h(x?k)為稠密向量，因此b也是高維稠密向量。 此步驟采用CUDA 中的運算庫cuSPARSE［21］即可并行高效計算出矩陣A和向量b的結(jié)果，根據(jù)CUDA 的描述，計算速度比純CPU 替代產(chǎn)品快2 至5 倍。

矩陣線性方程組計算部分最為耗時。從式（9）可以看出，矩陣A是矩陣轉(zhuǎn)置、誤差權(quán)重矩陣和矩陣本身之間的乘積，是對稱正定矩陣。當(dāng)系統(tǒng)規(guī)模較小時，矩陣規(guī)模較小，適合使用LU 分解等直接法對線性方程組進行計算。常用的高性能線性方程組求解庫SuperLU［22］經(jīng)過了高度優(yōu)化，屬于直接法，效率較高。但對于大規(guī)模稀疏線性方程組，矩陣條件數(shù)較高，直接法的計算效率難以滿足要求，適合在GPU 上使用迭代法進行并行計算。

2 迭代法基本原理

2.1 Krylov 子空間法

Krylov 子空間法是求解大規(guī)模線性方程組的首選方法。求解一般線性方程組：

投影方法的基本思想是從一個維數(shù)更小的子空間Km內(nèi)尋找近似的解。這個子空間Km被稱為搜索空間，其維度為m。

此時，設(shè)置m個約束，要求殘差矢量r滿足m個正交條件，即Petrov-Galerkin 條件：

式中：Lm為另一個m維的子空間，被稱為約束空間。當(dāng)Lm=Km時，該方法為正交投影法，否則為斜投影法。

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

式中：x?為迭代值，初始?xì)埐顁0=b?Ax0。

在Krylov 子空間方法中，搜索空間Km就是Krylov 子空間，定義為：

式中：r可選為初始?xì)埐顁0，Krylov 子空間方法就是在Krylov 子空間中尋找近似解。

當(dāng)選用不同的約束空間Lm時，對于迭代過程會有比較大的影響。主要可以分為如下方法。

1）正交投影法：Lm=Km(A，r0)，代表方法有共軛梯度法。該方法計算步驟較為簡單，但要求矩陣A對稱正定。

2）正交化方法：Lm=AKm(A，r0)，代表方法有廣義極小殘差方法。該方法數(shù)值穩(wěn)定性高，但計算量會隨迭代的步數(shù)線性增長。

3）雙正交化方法：Lm=AKm(AT，r0)，代表方法有穩(wěn)定的雙正交共軛梯度法。該方法適用性最廣且穩(wěn)定性較高，但計算步驟較為復(fù)雜，計算量比前2 種方法要大。

搜索空間Km的選取不會對最終狀態(tài)估計結(jié)果的精度有影響，其選取依據(jù)主要是狀態(tài)估計中矩陣A的特性，例如是否對稱、是否正定等。在狀態(tài)估計的系數(shù)矩陣A對稱正定情況下，可以使用正交投影法確定約束空間Lm，此方法計算量較小、步驟較簡單。本文采用共軛梯度法進行求解。從式（9）可以看出，矩陣A為對稱正定矩陣，滿足共軛梯度法的要求，能夠充分發(fā)揮其計算步驟簡單、效率高的優(yōu)勢。

2.2 預(yù)處理方法

如果直接使用上述Krylov 子空間方法進行迭代，原始矩陣條件數(shù)過高，可能會出現(xiàn)收斂性差、迭代次數(shù)多等問題。采用合適的預(yù)處理方法可以降低矩陣條件數(shù)。預(yù)處理方法的原理是通過預(yù)處理子M將原始線性方程組轉(zhuǎn)化為另一個易于求解的線性方程組。

常見的預(yù)處理方法可分為如下幾類［23］。

1）左預(yù)處理：對M?1Ax=M?1b使用迭代法。

2）右預(yù)處理：對AM?1u=b使用迭代法，其中x=M?1u。

3）分裂預(yù)處理：對M1?1u=M1?1b使用迭代法，預(yù)處理子M=M1M2。

對于大規(guī)模電力系統(tǒng)，其矩陣條件數(shù)較高，利用預(yù)處理技術(shù)可以減少迭代次數(shù)，易于問題求解。不完全LU 分解預(yù)處理方法適用范圍廣，且在文獻［18］中已經(jīng)被證明在潮流計算中有較好的預(yù)處理效果。因此，本文選用不完全LU 預(yù)處理方法。

2.3 預(yù)處理共軛梯度迭代法

考慮對稱正定矩陣A，并假設(shè)有一個預(yù)處理子M可供使用，且M是對稱正定的。為了保持對稱性，注意到M?1A對M內(nèi)積(x，y)M≡(Mx，y)=(x，My)是自伴的，因為

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

所以，一種選擇是用M內(nèi)積代替共軛梯度迭代法中的歐氏內(nèi)積，如果按這種新內(nèi)積重寫共軛梯度迭代法，用rj=b?Axj表示原殘差，而用zj=M?1rj表示預(yù)處理方程組的殘差，可以得到本文所使用的預(yù)處理共軛梯度迭代法，且不需要顯式地計算M內(nèi)積。預(yù)處理共軛梯度迭代法的具體步驟見3.1 節(jié)。

3 GPU 并行加速狀態(tài)估計算法

3.1 算法步驟

本文提出一種基于預(yù)處理共軛梯度迭代法的電力系統(tǒng)WLS 狀態(tài)估計算法，其計算步驟如下。

步驟1：初始化，形成節(jié)點導(dǎo)納矩陣，給狀態(tài)變量賦初值形成x?0。

步驟2：設(shè)置迭代變量k=0 和最大迭代次數(shù)kmax。

步驟3：根據(jù)當(dāng)前狀態(tài)變量，計算雅可比矩陣H()，計算公式如文獻［19］附錄中所示。

步驟4：在GPU 上利用cuSparse 庫，根據(jù)式（9）計算矩陣A和向量b。

步驟5：在GPU 上求解線性方程組Ax=b，利用矩陣A對稱正定的特點，使用預(yù)處理共軛梯度法進行迭代求解［24］，具體做法如下。

①對矩陣A進行ILU（0）分解，ILU（0）分解是不完全LU 分解的一種形式，形成矩陣A的預(yù)處理子：

式中：L和U分別為ILU（0）分解出的上三角矩陣和下三角矩陣。設(shè)置迭代次數(shù)i=0 和最大迭代次數(shù)imax，設(shè)x的初始猜測為x0，計算初始?xì)埐顁0及其2 范數(shù)‖r0‖。

②根據(jù)L和U求解方程組Mz=ri，其中ri為迭代次數(shù)為i時的殘差，得到z后計算ρi=(ri，z)。

③如果i=0，則令pi=z，否則令β=ρi/ρi?1，并計算pi=z+βpi?1。

步驟6：根據(jù)步驟5 中解出的x即可得到牛頓法迭代的修正量Δx?k，根據(jù)式（7）得到狀態(tài)變量x?k+1。

步驟7：令k=k+1，若不滿足式（8）且k沒有達到kmax則轉(zhuǎn)至步驟3，否則退出狀態(tài)估計過程。

3.2 算法特點

1）算法采用ILU（0）預(yù)處理方法，保證了殘差矩陣滿足ILU（0）的分解條件。這種預(yù)處理方法產(chǎn)生的預(yù)處理子不會注入非零元。經(jīng)過預(yù)處理后，能夠保證預(yù)處理子的稀疏性。在大型稀疏矩陣中，保證預(yù)處理子的稀疏性可以在矩陣運算中節(jié)省計算量和內(nèi)存，同時利用預(yù)處理子與原矩陣相同的稀疏性能夠更快地加速迭代法的計算。從迭代步驟可以看出，計算只涉及稀疏矩陣向量線性乘法，且由于沒有更多的非零元注入，理論上不會產(chǎn)生內(nèi)存泄漏。

2）算法采用了共軛梯度法，WLS 狀態(tài)估計中線性方程組系數(shù)矩陣A對稱正定的特性使共軛梯度法苛刻的適用條件得到滿足。在大型稀疏線性方程組求解過程中，迭代法的計算效率更高，同時共軛梯度法作為迭代法中計算步驟最簡單的方法，具有最少的計算量和最高的計算效率。

3）算法采用了GPU 并行計算架構(gòu)，充分利用CUDA 的高性能矩陣向量計算技術(shù)和cuSparse 庫的乘法運算，加速矩陣A和向量b的形成。同時，在共軛梯度法迭代過程中，使用GPU 快速計算各個中間變量，保證了迭代法的快速計算。

4 算例分析

為了分析本文所提算法的有效性，選取部分算例進行測試。其中，CASE300 是IEEE 標(biāo)準(zhǔn)測試算例，分別將10、20、50、100 個CASE300 算例以平衡節(jié)點為中心進行拼接，可以得到拼接算例CASE2991、 CASE5981、 CASE14951和CASE29901。本文的量測數(shù)據(jù)是在潮流計算的基礎(chǔ)上加入2%的高斯噪聲所形成的。算例測試環(huán)境為64 bit 的Windows 10 操作系統(tǒng)，CPU 型號為Intel Core i7-5820K，GPU 型號為NVIDIA GeForce GTX1080。其中，牛頓法的迭代精度設(shè)置為10?3，線性方程組求解的迭代精度設(shè)置為10?10。

4.1 計算耗時對比

為了分析本文算法的并行加速效果，將本文算法與直接法進行比較。其中，直接法使用高性能SuperLU 庫，能夠較快地使用列選主元高斯消去法求解高性能稀疏線性方程組，計算結(jié)果如表1所示。

表1 本文算法與直接法計算效果對比Table 1 Comparison of calculation results between method of this paper and direct method

從表1 可以看出，本文算法隨著算例規(guī)模的增大，計算加速效果越來越明顯。當(dāng)算例規(guī)模較小時，如CASE300 算例中本文算法計算速度比直接法更慢，這是由于直接法在線性方程組維數(shù)較小時LU分解速度較快。但當(dāng)系統(tǒng)規(guī)模增大時，本文算法受問題規(guī)模的影響較小，而此時直接法的計算時間會急劇增大。上述算例證明在大規(guī)模電力系統(tǒng)中，本文算法能夠大幅提高計算速度，與直接法相比優(yōu)勢明顯。

4.2 計算精度對比

為了分析本文算法的估計精度誤差，將本文算法與直接法進行比較。通過計算本文算法和直接法估計結(jié)果與真值之間的平均相對誤差來評估計算精度，計算結(jié)果如表2 所示。

表2 本文算法與直接法計算精度對比Table 1 Comparison of calculation precision between method of this paper and direct method

從表2 可以看出，使用本文算法狀態(tài)估計結(jié)果的平均相對誤差在2%以內(nèi)，與直接法精度相比幾乎沒有差別。從計算步驟上來看，本文方法與直接法的不同之處在于求解線性方程的步驟，而由于線性方程組的迭代精度設(shè)置較為嚴(yán)苛，設(shè)為10?10，求得方程的解與直接法相差無幾。算例結(jié)果表明本文算法不會對狀態(tài)估計精度造成過大的影響。

4.3 預(yù)處理技術(shù)有效性分析

良好的預(yù)處理有助于減少迭代次數(shù)，改善迭代法的穩(wěn)定性。為了驗證本文所采用的ILU（0）預(yù)處理方法的有效性，對預(yù)處理前后矩陣的條件數(shù)和迭代次數(shù)進行分析，計算結(jié)果如表3 所示。其中，預(yù)處理條件數(shù)和迭代次數(shù)取各次迭代的平均值。

表3 ILU(0)預(yù)處理方法效果分析Table 3 Effect analysis of ILU(0)preprocessing method

從表3 可以看出，經(jīng)過ILU（0）預(yù)處理后，預(yù)處理子的條件數(shù)大幅下降。在迭代法中，每次求解都需要計算預(yù)處理子所形成的線性方程組，條件數(shù)下降意味著更容易對該方程組進行求解。因此，預(yù)處理后共軛梯度法的迭代次數(shù)也隨之下降，甚至可以改善原本不收斂的線性方程組的收斂性。

減少迭代次數(shù)能夠縮短計算時間，且該預(yù)處理方法不注入非零元，不會因稀疏性被破壞而大幅增加計算時間。算例結(jié)果表明，本文所采用的預(yù)處理方法能夠有效減少共軛梯度法的迭代次數(shù)，提高本文算法的計算效率。

4.4 內(nèi)存占用測試

在大型電力系統(tǒng)計算中，內(nèi)存占用是一個重要的性能指標(biāo)。為了驗證本文算法的內(nèi)存占用情況，對內(nèi)存占用和顯存占用進行測試，計算結(jié)果如表4所示。

表4 本文算法內(nèi)存占用測試Table 4 Memory footprint test of method in this paper

本文算法為了節(jié)約內(nèi)存，矩陣存儲方式為壓縮矩陣形式。從表4 可以看出，即使對于節(jié)點數(shù)約3 萬的大規(guī)模電力系統(tǒng)，本文算法內(nèi)存占用不超過700 Mbit，顯存占用不超過600 Mbit。整體而言，本文所用算法內(nèi)存和顯存占用峰值都較低，在普通的計算機上即可進行計算。

5 結(jié)語

針對大型電力系統(tǒng)狀態(tài)估計速度較慢的問題，本文提出了一種基于預(yù)處理共軛梯度迭代法的電力系統(tǒng)狀態(tài)估計算法。本文算法是在WLS 狀態(tài)估計法的基礎(chǔ)上，針對牛頓法求解過程中矩陣乘法和大規(guī)模稀疏線性方程組求解效率較低的特點，在GPU上進行并行加速。由于WLS 狀態(tài)估計中線性方程組稀疏矩陣為對稱正定矩陣，滿足共軛梯度法的使用條件，本文提出采用不完全LU 分解預(yù)處理方法和共軛梯度迭代法進行線性方程組的求解。算例分析驗證了本文算法計算速度快、內(nèi)存占用低的特點，并驗證了預(yù)處理方法的有效性。綜上所述，本文所提算法能夠滿足大規(guī)模電力系統(tǒng)狀態(tài)估計的實時性要求，未來可考慮多機多卡下的分布式計算，進一步提高狀態(tài)估計的計算效率。