王克文，冶夢雨，劉艷紅

(鄭州大學 電氣工程學院，河南 鄭州 450001)

0 引言

特征值分析法以狀態(tài)空間模型為基礎，是電力系統(tǒng)小干擾穩(wěn)定性分析的常用方法之一[1]。隨著現(xiàn)代電力系統(tǒng)規(guī)模不斷擴大[2]，求解狀態(tài)空間矩陣所需計算量急劇增大，國內外學者開始尋找快速形成狀態(tài)矩陣的方法。

目前，國外通常采用擬合的方法求解狀態(tài)矩陣，比如Saitoh等[3]用廣域監(jiān)測系統(tǒng)(WAMS)搜集各個離散時間點的值，并利用最小二乘法擬合出狀態(tài)矩陣。國內有兩種主流方法:一種是用解析的方法直接得出狀態(tài)矩陣，該方法計算速度很快，但是預處理較困難；另一種則是利用插入式模擬技術(PMT)[4]，先將所有系統(tǒng)元件轉化為兩種基本傳輸模塊，再利用傳輸模塊兩端的節(jié)點編號組成關聯(lián)矩陣，最終得到狀態(tài)空間方程。第2種方法較為靈活，適用于對元件內部模塊間變量的監(jiān)控和靈敏度計算，但是矩陣運算較多，漸漸陷入了矩陣“維數(shù)災難”。因此，一些學者在該方法的基礎上進行了改進，比如羅丹等[5]將狀態(tài)矩陣形成過程中的相關算式進行QR優(yōu)化和重組，通過矩陣降階求逆的方法提高計算效率；胡崴[6]則是利用Hager算法縮減矩陣的外形，再引入并行計算提高計算效率，但是Hager縮減過程中沒有與并行計算相結合，此過程本身也會耗費一定時間。除此之外，為了縮短特征值的計算時間，有些學者利用改進的Rayleigh商逆迭代法來獲得機電模式的特征值:首先用相關因子將機電模式分組，再參照模型降階的方法，計算出大規(guī)模電力系統(tǒng)機電模式特征值[7]。

在PMT的基礎上充分利用相關矩陣的稀疏特性及Open MP技術的并行計算功能：首先將求解過程中的關鍵矩陣進行ILUTP預處理，合理舍棄掉一些非對角非零元素，使系數(shù)矩陣譜分布更集中，再采用Krylov子空間迭代法中的BICGSTAB算法，對得到的預處理模塊進行快速迭代計算。ILUTP技術與BICGSTAB算法均能很好地與Open MP結合使用，進一步加快了迭代計算速度。稀疏矩陣存儲采用行壓縮矩陣存儲的方法降低存儲量，減小計算量，從而提高整體計算效率。

1 BICGSTAB算法及ILUTP技術

1.1 BICGSTAB算法

Krylov子空間迭代法[8]適用于求解大型線性方程組。給定線性方程組：

Ax=b。

(1)

式中：A∈Rn×n；x，b∈Rn。

設Km為m維子空間，一般投影方法是從m維仿射子空間x0+Km中尋找上式的近似解xm，x0為迭代初值，使相應的殘差滿足Petrov-Galerkin條件:

rm=b-Axm。

(2)

其中，要求rm與預設的約束空間Lm正交[9]。

記Km=Km(A,r0)為Krylov子空間，定義為：

Km(A,r0)=span{r0,Ar0,A2r0,…,Am-1r0}。

(3)

約束空間Lm的選擇對迭代過程具有重要影響，本文中矩陣具有非對稱性，故選用雙正交化方法，取Lm=Km(AT,r0)，其中BICG法是雙正交化方法中最典型的算法。

BICGSTAB算法[10]是在BICG的基礎上發(fā)展起來的，其避免了BICG中對AT的計算，收斂速度比BICG更快。

(4)

(5)

BICGSTAB算法能夠在保證精度和穩(wěn)定性的同時，達到快速收斂的效果。但是其收斂速度非常依賴于系數(shù)矩陣特征值的分布，所以在一些情況中，直接使用BICGSTAB算法會導致收斂速度很慢，甚至根本不能收斂[11]。可使用預條件技術將原線性方程組轉化為更利于迭代收斂的線性方程組[11]，保證迭代過程不中斷，提高收斂速度。

1.2 ILUTP預處理

ILUTP方法[15]的優(yōu)勢在于其可以降低分解因子中非零元素的個數(shù)，從而縮減不完全分解時間和預條件處理的迭代時間；同時，當矩陣性質較差時，可以通過調整參數(shù)τ與p的值來提高分解因子的質量，增強迭代有效性。但是在特定情況下，需要不斷嘗試才能預先選擇出合適的參數(shù)τ和p，否則精度將無法得到保證。因此，最優(yōu)值應當結合實際問題和矩陣特點來合理選定。

引入一個名為droptol(閾值)的參數(shù)以確定是否改變變量值，使得矩陣的系數(shù)譜變得更加集中，從而將矩陣轉換為對角占優(yōu)形矩陣，可提高計算效率。但是由于舍棄掉了一些非零元素，計算準確度無法保證，因此又引入一個名為tol(絕對誤差限)的參數(shù)以保證計算精度。通過不斷嘗試，在速度和精度之間找出最佳平衡點。

2 狀態(tài)空間方程的建立及特點分析

2.1 插入式建模技術

在插入式建模技術中，將電力系統(tǒng)描述為傳輸模塊集合：非狀態(tài)變量傳輸塊(零階)、狀態(tài)變量傳輸塊(一階)和5類基本參數(shù)k、ka、kb、Ta、Tb，再根據(jù)各個系統(tǒng)元件和電力網絡中傳輸模塊兩端的節(jié)點編號列出關聯(lián)矩陣L:

(6)

式中：Xi和Mi分別為狀態(tài)變量和非狀態(tài)變量模塊的輸入;X和M則為狀態(tài)變量和非狀態(tài)變量模塊的輸出;R、Y對應系統(tǒng)輸入、輸出的列向量。

第n個零階和一階的傳輸塊方程可表示為：

(7)

那么，所有零階和一階的傳輸塊方程可表示為：

(8)

式(7)、式(8)中，kn為第n個零階傳輸塊的傳輸增益；K為由所有的kn組成的對角矩陣；Ka、Kb和Kt分別為所有一階模塊參數(shù)kan/Tbn、kbn/Tbn和Tan/Tbn共同組成的對角矩陣。

由式(6)～(8)消去Xi、Mi和M，得到電力系統(tǒng)狀態(tài)空間方程如下:

(9)

其中，各系數(shù)矩陣的具體表達可詳見文獻[4]。

狀態(tài)矩陣A的表達式為：

A=(I-Kt(L1-L3L9-1))-1×

(Ka(L1-L3L9-1L7)-Kb)。

(10)

2.2 矩陣分析

基于插入式建模中節(jié)點電壓和電流是非狀態(tài)變量，節(jié)點導納矩陣可以直接被插入到分塊子矩陣L9中，得到下式：

(11)

式中：G、B分別為導納矩陣的實部和虛部；KVR、KVJ、KIR、KIJ和KO為非狀態(tài)模塊的參數(shù)矩陣；I為單位陣。

形成狀態(tài)空間方程時，內存的需求主要取決于L9子矩陣，求解L9矩陣的逆矩陣是求解A矩陣的關鍵，也是計算過程中耗時最長的環(huán)節(jié)，因此加速求解L9逆矩陣可加快狀態(tài)空間方程的形成。

L9系數(shù)矩陣為大型稀疏矩陣，其稀疏度隨著數(shù)據(jù)增多而變高，所以可充分利用其稀疏特性優(yōu)化求逆過程。對于大型非對稱的稀疏矩陣，首先考慮用BICGSTAB 算法迭代。圖1為L9矩陣的數(shù)據(jù)分布示意圖，“×”表示非零元素，其余為零元素。觀察圖1可以看出L9矩陣條件數(shù)很大且為非對角占優(yōu)矩陣，直接用BICGSTAB算法無法收斂，因此本文最終選用ILUTP結合BICGSTAB求解方程，迭代計算則采用Open MP技術的并行計算功能。

圖1 L9矩陣Figure 1 The L9 matrix

2.3 算法分析

并行處理是將原有單線程執(zhí)行的程序變?yōu)槎嗑€程執(zhí)行，增加線程并不會影響算法的時間復雜度，但是會提高執(zhí)行效率。優(yōu)化算法與原算法相比，與Open MP的并行功能結合度更高，因此執(zhí)行效率更高，計算時間會有明顯縮短。

3 加速求解狀態(tài)矩陣的程序實現(xiàn)

本文求解狀態(tài)矩陣A的步驟如下：

Step1輸入稀疏矩陣A，矩陣采用行壓縮存儲；

Step2ILU分解，得到L、U矩陣;

Step3BICGSTAB，求解線性方程組Ax=b;

Step4求解方程組采用Open MP并行;

Step5輸出結果，將L9的逆矩陣代入原程序，求A矩陣及特征值。

3.1 行壓縮存儲

鑒于稀疏矩陣A的特點，矩陣采用行壓縮存儲(CRS)[16]，該格式只顯示保留每行第一個非零元素的位置，具體在圖2所示例子中可以看到：假設有稀疏矩陣B，創(chuàng)建一個浮點型數(shù)組Val和兩個整型數(shù)組Colind和Rowp。Val數(shù)組的大小為矩陣B中非零元素的個數(shù)，保存矩陣B的非零元素(按從上往下，從左往右的行遍歷方式訪問元素)；Colind數(shù)組和Val數(shù)組一樣，大小為矩陣B的非零元素的個數(shù)，保存Val數(shù)組中元素的列索引；Rowp數(shù)組大小為矩陣B的行數(shù)，保存矩陣B的每行第一個非零元素在Val中的索引。該方法可以節(jié)省很多空間，只需要2nnz+n+1個存儲單元，而不是n2個單元，其中nnz為稀疏矩陣中非零元素的個數(shù)。

圖2 稀疏矩陣B的行壓縮示意圖Figure 2 Schematic of sparse matrix B using CRS

3.2 預處理參數(shù)設定

integer::ILU_mode=1

real*8::droptol=0.001

integer::lfil=10

real*8::permtol=0.5

integer::mbloc=0

通過閾值來舍棄數(shù)值小的元素，從而影響生成的ILU矩陣大小。這個值設定得太小會導致無法收斂，太大則無法保證精度。本文經過多次修改對比，最終選用droptol_def=0.001，tol=0.000 01達到理想的平衡點。下降公差lfil用來控制L和U中每行限定的最大非對角元素個數(shù)，填充參數(shù)設定得較小是為了保證計算精度。列交換容差permtol合理的值通常介于0.5和0.01，本文矩陣維數(shù)很大，因此將列交換容差設為0.5。圖3為預處理之后得到的L9矩陣的數(shù)據(jù)分布示意圖，經過ILUTP預處理之后，矩陣轉換為對角占優(yōu)型矩陣，一方面更有利于BICGSTAB算法的求逆迭代，另一方面更易與Open MP的并行功能相結合。

圖3 轉換后的L9矩陣Figure 3 The transformed L9 matrix

3.3 BICGSTAB算法實現(xiàn)

程序中BICGSTAB算法流程如下：

Step4xi=xi-1+αpi+ωis；

Step5若xi達到精度要求，則轉Step 7；

Step6ri=s-ωit，i=i+1，若i

Step7輸出xi，算法結束。

3.4 Open MP并行實現(xiàn)

Open MP作為并行提速計算的一種工具，適合處理獨立循環(huán)或者分開的子任務，程序編寫簡單，具有強擴展性，可支持Fortran等多種編程語言[8]。

在狀態(tài)矩陣形成的程序中有很多DO循環(huán)語句，如對矩陣L1、L3、L7、L9賦值、對矩陣L9求逆過程中的迭代等，這些任務在程序執(zhí)行過程中可以獨立執(zhí)行，因此結合Open MP進行并行處理。其中BICGSTAB算法中求解方程組時的并行語句如下：

!S|OMP Parallel

!S|OMP do Private(i,x0,x,bi)

doi=1,N

bi=0.d0;bi(i)=1.d0;x0=0.d0

call Precond(bi,x0)

x0=x0(Iperm(1:N))

call Bicgstab(bi，x0，x，max_it，tol)

a(1:N,i)=x(Iperm(1:N))

enddo

!S|OMP enddo

!S|OMP end Parallel

4 算例分析

某省電力網絡的兩個實際算例分別包含23臺發(fā)電機和98臺發(fā)電機，發(fā)電機均采用六階發(fā)電機模型，勵磁調節(jié)模塊與原動機調速模塊均為系統(tǒng)的實際參數(shù)。

程序采用Fortran語言編寫，編寫環(huán)境為Visual Studio 2010平臺、Intel Fortran 12，在2.14 GHz 主頻、32 GB 內存、8核計算機上運行實現(xiàn)。存儲方式為行壓縮存儲形式。

通過對兩種求解狀態(tài)矩陣方法的所用時間進行比較來驗證可行性。兩種方法分別為原方法和現(xiàn)方法，均利用插入式建模技術及Open MP的并行功能。其中,原方法為文獻[4]中基礎的PMT方法，矩陣求逆采用直接LU分解法，矩陣存儲采用三元組技術存儲；現(xiàn)方法為預處理的不完全LU分解法，矩陣求逆采用BIGSTAB算法迭代求逆，矩陣存儲采用行壓縮存儲技術。分別將原方法及現(xiàn)方法無并行所用時間、原方法及現(xiàn)方法加入并行所用時間進行比對，驗證現(xiàn)方法的加速效果。

算例1為23臺發(fā)電機，系統(tǒng)中包含90個節(jié)點，210條支路，狀態(tài)矩陣的階數(shù)為257階。算例1中矩陣L1、L3、L9及狀態(tài)矩陣A的階數(shù)、非零元素個數(shù)和對應的稀疏程度如表1所示。本算例為了對比算法方面的改進，設為單線程執(zhí)行，各方法運行時間的對比如表2所示。

表1 算例1中各矩陣的信息Table 1 Information of each matrix in example 1

表2 算例1中各方法的運行時間Table 2 Run time of each method in example 1

由表2得出，在矩陣階數(shù)較少時，BICGSTAB算法及預處理技術的運用對計算效率并沒有非常明顯的提高。

算例2為98臺發(fā)電機，系統(tǒng)中包含2 394個節(jié)點，5 512條支路，狀態(tài)矩陣的階數(shù)為1 736階。其中矩陣L1、L3、L9及狀態(tài)矩陣A的矩陣規(guī)模、非零元素數(shù)和對應的稀疏程度如表3所示。算例2中加入了并行計算功能，為八線程執(zhí)行，各方法運行時間的對比如表4所示。

表3 算例2中各矩陣的信息Table 3 Information of each matrix in example 2

表4 算例2中各方法的運行時間Table 4 Run time of each method in example 2

綜合各表中數(shù)據(jù)可得：當矩陣規(guī)模較小時，ILU預處理及BICGSTAB算法對運行時間有一定改善，但效果并不明顯；當系統(tǒng)規(guī)模較大時，L9矩陣稀疏度增高，越來越接近于1，稀疏技術的利用率變高。

通過原方法無并行和現(xiàn)方法無并行的運算時間對比，說明ILUTP及BICGSTAB算法可以提高迭代速度；引入并行計算后，加速效果更加明顯，并行加速比接近于3.

5 結論

在插入式模擬技術的基礎上，對L9矩陣的求逆過程進行優(yōu)化，利用矩陣的稀疏特性引入ILUTP預處理技術，并與BICGSTAB算法相結合，提高了迭代速度，使矩陣轉化為更有利于并行計算的形式；然后通過行壓縮存儲技術和Open MP中的并行功能，加速了整個求逆過程；最后分別對兩個不同規(guī)模的算例進行分析計算，驗證了本文方法能夠縮短狀態(tài)矩陣形成的時間，有利于電力系統(tǒng)小干擾穩(wěn)定性的計算與分析。

在之后的研究中，一方面可以嘗試不同的算法轉換矩陣形式，比如將L9矩陣變換成對角加邊的形式，再進行并行化處理；另一方面也可嘗試使用 Open MP+MPI 混合編程，利用不同的并行處理方式加快程序執(zhí)行。