李建斌，王鵬程，傅 侃，方 睿，董樹(shù)鋒

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

0 引言

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

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

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

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

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

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

1.1 WLS 狀態(tài)估計(jì)

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

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

定義殘差矢量為：

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

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

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

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

1.2 并行性分析

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

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

矩陣乘法步驟包括矩陣A和向量b的運(yùn)算，其計(jì)算公式如下。

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

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

2 迭代法基本原理

2.1 Krylov 子空間法

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

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

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

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

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

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

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

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

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

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

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

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

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

2.2 預(yù)處理方法

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

常見(jiàn)的預(yù)處理方法可分為如下幾類(lèi)［23］。

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

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

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

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

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

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

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

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

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

3.1 算法步驟

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

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

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

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

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

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

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

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

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

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

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

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

3.2 算法特點(diǎn)

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

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

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

4 算例分析

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

4.1 計(jì)算耗時(shí)對(duì)比

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

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

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

4.2 計(jì)算精度對(duì)比

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

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

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

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

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

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

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

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

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

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

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

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

5 結(jié)語(yǔ)

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