建立電力系統(tǒng)狀態(tài)空間方程的并行方法

王克文，冶夢(mèng)雨，劉艷紅

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

0 引言

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

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

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

1 BICGSTAB算法及ILUTP技術(shù)

1.1 BICGSTAB算法

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

Ax=b。

(1)

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

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

rm=b-Axm。

(2)

其中，要求rm與預(yù)設(shè)的約束空間Lm正交[9]。

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

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

(3)

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

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

(4)

(5)

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

1.2 ILUTP預(yù)處理

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

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

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

2.1 插入式建模技術(shù)

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

(6)

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

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

(7)

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

(8)

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

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

(9)

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

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

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

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

(10)

2.2 矩陣分析

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

(11)

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

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

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

圖1 L9矩陣Figure 1 The L9 matrix

2.3 算法分析

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

3 加速求解狀態(tài)矩陣的程序?qū)崿F(xiàn)

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

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

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

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

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

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

3.1 行壓縮存儲(chǔ)

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

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

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

integer::ILU_mode=1

real*8::droptol=0.001

integer::lfil=10

real*8::permtol=0.5

integer::mbloc=0

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

圖3 轉(zhuǎn)換后的L9矩陣Figure 3 The transformed L9 matrix

3.3 BICGSTAB算法實(shí)現(xiàn)

程序中BICGSTAB算法流程如下：

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

Step5若xi達(dá)到精度要求，則轉(zhuǎn)Step 7；

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

Step7輸出xi，算法結(jié)束。

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

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

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

!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 算例分析

某省電力網(wǎng)絡(luò)的兩個(gè)實(shí)際算例分別包含23臺(tái)發(fā)電機(jī)和98臺(tái)發(fā)電機(jī)，發(fā)電機(jī)均采用六階發(fā)電機(jī)模型，勵(lì)磁調(diào)節(jié)模塊與原動(dòng)機(jī)調(diào)速模塊均為系統(tǒng)的實(shí)際參數(shù)。

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

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

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

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

表2 算例1中各方法的運(yùn)行時(shí)間Table 2 Run time of each method in example 1

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

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

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

表4 算例2中各方法的運(yùn)行時(shí)間Table 4 Run time of each method in example 2

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

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

5 結(jié)論

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

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