基于OpenMP并行簡(jiǎn)約空間內(nèi)點(diǎn)法的暫態(tài)穩(wěn)定緊急控制

王 云，江全元

（浙江大學(xué) 電氣工程學(xué)院，浙江 杭州 310027）

0 引言

暫態(tài)穩(wěn)定緊急控制是維持電力系統(tǒng)安全穩(wěn)定運(yùn)行的重要手段。暫態(tài)穩(wěn)定緊急控制問(wèn)題通常建模為考慮暫態(tài)穩(wěn)定約束和系統(tǒng)經(jīng)濟(jì)運(yùn)行的最優(yōu)化問(wèn)題。

目前，求解暫態(tài)穩(wěn)定緊急控制問(wèn)題的算法主要有試湊法或啟發(fā)式算法[1-3]、基于暫態(tài)能量函數(shù)的直接法和基于最優(yōu)控制理論的優(yōu)化算法。暫態(tài)穩(wěn)定緊急控制問(wèn)題的難點(diǎn)之一在于暫態(tài)穩(wěn)定約束的描述，基于能量函數(shù)的方法在暫態(tài)穩(wěn)定預(yù)防控制和暫態(tài)穩(wěn)定緊急控制方面已有廣泛的研究[4-11]。時(shí)域仿真法是解決暫態(tài)穩(wěn)定緊急控制的另一種方法，基于最優(yōu)控制理論，各種緊急決策計(jì)算均描述為以控制代價(jià)最小為目標(biāo)、以維持系統(tǒng)安全穩(wěn)定性為約束條件的最優(yōu)控制問(wèn)題，這類最優(yōu)控制問(wèn)題實(shí)質(zhì)上就是最優(yōu)參數(shù)選取問(wèn)題。文獻(xiàn)[12-13]將切負(fù)荷問(wèn)題描述為最優(yōu)控制問(wèn)題，提出了一種系統(tǒng)化地決定最小切負(fù)荷量及其在各切負(fù)荷點(diǎn)分配的算法。文獻(xiàn)[14]基于時(shí)域仿真得到的系統(tǒng)受擾軌跡給出暫態(tài)穩(wěn)定緊急控制的非線性模型，然后采用近似規(guī)劃法和擬貪婪法求解。

簡(jiǎn)約空間技術(shù)主要用來(lái)求解自由度較小的非線性規(guī)劃問(wèn)題，已在化工領(lǐng)域獲得了廣泛應(yīng)用[15]。文獻(xiàn)[16-17]用簡(jiǎn)約空間內(nèi)點(diǎn)算法求解暫態(tài)穩(wěn)定最優(yōu)潮流問(wèn)題，計(jì)算效率獲得極大的提升。文獻(xiàn)[18]應(yīng)用簡(jiǎn)約空間二次規(guī)劃法求解暫態(tài)電壓穩(wěn)定緊急控制問(wèn)題，顯著提高了計(jì)算效率。由于暫態(tài)穩(wěn)定緊急控制問(wèn)題中可切負(fù)荷/發(fā)電機(jī)節(jié)點(diǎn)相對(duì)于問(wèn)題規(guī)模非常小，即該類問(wèn)題自由度很低，本文運(yùn)用簡(jiǎn)約空間內(nèi)點(diǎn)算法進(jìn)行求解。此外，算法采用C++編程實(shí)現(xiàn)，對(duì)于算法的關(guān)鍵耗時(shí)環(huán)節(jié)，充分利用多線程并行技術(shù)。多個(gè)算例測(cè)試結(jié)果表明，本文提出的基于OpenMP并行簡(jiǎn)約空間內(nèi)點(diǎn)法在求解緊急控制問(wèn)題中，計(jì)算時(shí)間明顯縮短，占用內(nèi)存減小，能夠求解大規(guī)模系統(tǒng)的緊急控制問(wèn)題。

1 暫態(tài)穩(wěn)定緊急控制數(shù)學(xué)模型

電力系統(tǒng)最優(yōu)切機(jī)切負(fù)荷控制問(wèn)題是一個(gè)大規(guī)模動(dòng)態(tài)優(yōu)化問(wèn)題，其一般形式為：

其中，x為優(yōu)化變量；F為優(yōu)化問(wèn)題目標(biāo)函數(shù)；H為含微分方程等式約束；G為含微分方程不等式約束；、分別為不等式約束上、下限。

采用隱式梯形積分方法對(duì)上述動(dòng)態(tài)優(yōu)化問(wèn)題差分化，動(dòng)態(tài)優(yōu)化問(wèn)題（1）可轉(zhuǎn)化為如下非線性規(guī)劃問(wèn)題：

1.1 目標(biāo)函數(shù)

暫態(tài)穩(wěn)定緊急控制的目的是切除最少的發(fā)電機(jī)和負(fù)荷來(lái)保持系統(tǒng)暫態(tài)穩(wěn)定，因此，本文將優(yōu)化目標(biāo)函數(shù)設(shè)為切機(jī)和切負(fù)荷量加權(quán)最?。?/p>

其中，u=[uG，uL]，為優(yōu)化控制變量，uG、uL分別為各可切發(fā)電機(jī)節(jié)點(diǎn)切機(jī)比例組成的向量和可切負(fù)荷節(jié)點(diǎn)切負(fù)荷比例組成的向量；cG、cL分別為切機(jī)、切負(fù)荷控制代價(jià)的權(quán)值，用來(lái)表征可切發(fā)電機(jī)和負(fù)荷的重要性，本文中將權(quán)值均取為1，即所有可切發(fā)電機(jī)、可切負(fù)荷重要性相同；PG、PL分別為各可切發(fā)電機(jī)節(jié)點(diǎn)裝機(jī)容量和各可切負(fù)荷節(jié)點(diǎn)總負(fù)荷。

1.2 等式約束

a.發(fā)電機(jī)動(dòng)態(tài)方程。

電力系統(tǒng)經(jīng)典模型對(duì)于系統(tǒng)第一擺穩(wěn)定分析非常有效實(shí)用[19]。本文中發(fā)電機(jī)模型采用經(jīng)典模型，負(fù)荷采用恒阻抗模型，因此，各發(fā)電機(jī)可由如下兩階微分方程描述：

其中，i=1，2，…，NG（NG為發(fā)電機(jī)節(jié)點(diǎn)數(shù)）；為發(fā)電機(jī)暫態(tài)電勢(shì)，經(jīng)典模型下為恒定值，E′i、x′di分別為發(fā)電機(jī)節(jié)點(diǎn) i暫態(tài)電勢(shì)、直軸暫態(tài)電抗，分別為發(fā)電機(jī)節(jié)點(diǎn)i在t時(shí)刻的功角、節(jié)點(diǎn)電壓、節(jié)點(diǎn)電壓實(shí)部、節(jié)點(diǎn)電壓虛部；TJi、Pmi分別為發(fā)電機(jī)節(jié)點(diǎn)i慣性時(shí)間常數(shù)和有功出力；ωti為發(fā)電機(jī)節(jié)點(diǎn)i在t時(shí)刻的角速度；ωs為同步角速度有名值。

b.電網(wǎng)方程。

將發(fā)電機(jī)表示成電流源形式，并將發(fā)電機(jī)等值導(dǎo)納和負(fù)荷等值導(dǎo)納并入網(wǎng)絡(luò)，電網(wǎng)方程可寫成：

其中，Ixti、Ityi分別為t時(shí)刻發(fā)電機(jī)節(jié)點(diǎn)i的注入電流實(shí)部和虛部；Gij、Bij分別為節(jié)點(diǎn)i和節(jié)點(diǎn)j節(jié)點(diǎn)導(dǎo)納實(shí)部和虛部；nB為系統(tǒng)節(jié)點(diǎn)數(shù)。

對(duì)于可切機(jī)的發(fā)電機(jī)節(jié)點(diǎn)，有：

對(duì)于不可切機(jī)的發(fā)電機(jī)節(jié)點(diǎn)，有：

對(duì)于可切負(fù)荷節(jié)點(diǎn)，有：

對(duì)于不可切負(fù)荷節(jié)點(diǎn)，有：

其中，G′ii、B′ii分別為發(fā)電機(jī)等值導(dǎo)納和負(fù)荷等值導(dǎo)納并入網(wǎng)絡(luò)前節(jié)點(diǎn)自導(dǎo)納實(shí)部和虛部；GLi、BLi分別為負(fù)荷等值導(dǎo)納的實(shí)部和虛部；uGi、uLi分別為切機(jī)點(diǎn)切機(jī)比例和切負(fù)荷點(diǎn)切負(fù)荷比例。

1.3 不等式約束

a.控制變量上下界約束：b.暫態(tài)穩(wěn)定約束。

暫態(tài)穩(wěn)定判據(jù)采用中性慣量形式表示：

2 基于OpenMP并行簡(jiǎn)約空間內(nèi)點(diǎn)法

2.1 簡(jiǎn)約空間內(nèi)點(diǎn)算法

采用預(yù)測(cè)－校正內(nèi)點(diǎn)法求解非線性規(guī)劃問(wèn)題（2），其一階最優(yōu)必需條件（KKT條件）求解包括預(yù)測(cè)和校正2步，其主要計(jì)算量在求解以下對(duì)稱線性方程組[20]：

簡(jiǎn)約空間技術(shù)適用于自由度較低的非線性規(guī)劃問(wèn)題，用來(lái)求解對(duì)稱線性方程組（12）。將變量x分成程空間DY和零空間DZ兩部分，DY由m個(gè)非獨(dú)立變量組成，DZ由n-m個(gè)獨(dú)立變量組成。變量x的解可寫成：

其中，pY?DY；pZ?DZ；Y?Rn×m、Z?Rn×（n-m）分別稱為程空間和零空間基矩陣，[Y Z]非奇異。零空間基矩陣Z滿足：

則 pY、pZ、Δy 分別如下求得：

詳細(xì)的簡(jiǎn)約空間內(nèi)點(diǎn)算法原理及實(shí)施細(xì)則可參考文獻(xiàn)[16]。采用簡(jiǎn)約空間內(nèi)點(diǎn)算法求解緊急控制問(wèn)題時(shí)，方程組（15）—（17）的求解通過(guò)調(diào)用 KLU[21]計(jì)算，其對(duì)矩陣C的分解可復(fù)用，每次迭代只需做1次LU分解。算例分析表明，簡(jiǎn)約空間下緊急控制算法的主要計(jì)算時(shí)間集中在矩陣B和Zm的計(jì)算上，本文通過(guò)多線程技術(shù)提高這部分的計(jì)算效率。

2.2 OpenMP并行算法

2.2.1 OpenMP編程模型

為充分利用多核處理器計(jì)算平臺(tái)的計(jì)算性能，針對(duì)算法可解耦部分利用多線程技術(shù)實(shí)現(xiàn)并行計(jì)算從而加快計(jì)算速度。OpenMP[22]提供了一套共享存儲(chǔ)體系結(jié)構(gòu)下的多線程編程模型，從而實(shí)現(xiàn)本文算法多線程并行技術(shù)。

如圖1所示，OpenMP采用Fork-Join的并行執(zhí)行模型：程序首先以單線程（主線程）開始，類似一個(gè)串行程序，當(dāng)遇到并行結(jié)構(gòu)，主線程會(huì)產(chǎn)生一組線程（Fork動(dòng)作），每個(gè)線程都執(zhí)行并行動(dòng)態(tài)擴(kuò)展域中的代碼；并行結(jié)構(gòu)執(zhí)行完后，只有主線程繼續(xù)執(zhí)行，其他所有的線程結(jié)束執(zhí)行（Join動(dòng)作）。本文程序并行化部分，系統(tǒng)給每個(gè)線程分配一個(gè)CPU核心。

圖1 Fork-Join模型Fig.1 Fork-Join model

2.2.2 并行簡(jiǎn)約空間內(nèi)點(diǎn)算法

在簡(jiǎn)約空間內(nèi)點(diǎn)算法框架下，多線程技術(shù)可通過(guò)下面2種方式實(shí)現(xiàn)。

a.矩陣Zm的多線程求解。

矩陣N是由n-m個(gè)m維列向量構(gòu)成，在回代過(guò)程，矩陣N各個(gè)列向量回代過(guò)程相互獨(dú)立。將N分解成p（p為線程數(shù)）個(gè)m維列向量組成的矩陣，即N=[N1N2… Np]?？蓪⑦@p個(gè)矩陣分配給各計(jì)算線程實(shí)現(xiàn)并行回代，從而加快計(jì)算速度，如圖2所示。KLU解法器是一個(gè)串行解法器，不支持并行計(jì)算。本文對(duì)KLU代碼進(jìn)行了修改，使其按圖2所示支持多線程回代。

圖2 Zm多線程計(jì)算Fig.2 Multi-thread calculation of Zm

b.矩陣B的多線程計(jì)算。

本文算法采用C++編程實(shí)現(xiàn)，算法實(shí)施過(guò)程中涉及到較多矩陣操作，尤其是B的計(jì)算占據(jù)較多時(shí)間。目前有許多高性能線程代數(shù)庫(kù)可供調(diào)用，很多都支持多處理器平臺(tái)下的多線程并行。在本文中，基礎(chǔ)線程代數(shù)庫(kù)（BLAS）被用來(lái)作一些矩陣與矩陣相乘運(yùn)算和矩陣與向量相乘運(yùn)算，在計(jì)算矩陣B時(shí)主要調(diào)用DCSRMM和DGEMM函數(shù)進(jìn)行求解，這2個(gè)函數(shù)均支持OpenMP下多線程計(jì)算，程序?qū)崿F(xiàn)時(shí)通過(guò)調(diào)用MKL_NUM_THREADS函數(shù)設(shè)置線程數(shù)，開啟多線程從而加快計(jì)算速度。此外，線性代數(shù)解法器（LAPACK）被用來(lái)求解稠密線性方程組（18）。

2.3 算法實(shí)施流程

采用本文提出的算法求解暫態(tài)穩(wěn)定緊急控制問(wèn)題的主要步驟如下：

a.讀入算例數(shù)據(jù)，包括故障信息；

b.對(duì)故障后系統(tǒng)進(jìn)行暫態(tài)穩(wěn)定仿真計(jì)算，若系統(tǒng)保持穩(wěn)定，不需要切機(jī)切負(fù)荷，跳到步驟g；

c.按照第1節(jié)內(nèi)容對(duì)暫態(tài)穩(wěn)定緊急控制問(wèn)題數(shù)學(xué)建模；

d.對(duì)系統(tǒng)變量及收斂條件初始化，k=0，算法優(yōu)化開始；

制定智能化建筑用電節(jié)能設(shè)計(jì)時(shí)，一定要整體把握用電設(shè)施的布局分配、負(fù)荷容量及功率大小等信息，以此為基礎(chǔ)，選擇恰當(dāng)、合理的節(jié)能供配電設(shè)備，保證用電設(shè)備正常工作的同時(shí)，最大限度地降低能耗的損失。

e.按照第2節(jié)提出的并行簡(jiǎn)約空間內(nèi)點(diǎn)算法計(jì)算變量x更新步長(zhǎng)；

f.判斷是否達(dá)到收斂條件，若未收斂，k=k+1，跳到步驟e；

g.算法退出，輸出結(jié)果。

3 算例分析

3.1 測(cè)試算例

本文采用基于OpenMP并行簡(jiǎn)約空間預(yù)測(cè)校正內(nèi)點(diǎn)法對(duì)4個(gè)算例進(jìn)行了計(jì)算，這4個(gè)算例是在IEEE標(biāo)準(zhǔn)算例的基礎(chǔ)上修改得到，根據(jù)切機(jī)量和切負(fù)荷量相對(duì)功角穩(wěn)定的靈敏度大小來(lái)確定可切發(fā)電機(jī)和負(fù)荷。表1給出4個(gè)測(cè)試算例的參數(shù)，其中nB、nL、NG、nG、Nload、nload分別表示系統(tǒng)節(jié)點(diǎn)數(shù)、系統(tǒng)線路數(shù)、系統(tǒng)發(fā)電機(jī)數(shù)、系統(tǒng)可切發(fā)電機(jī)數(shù)、系統(tǒng)負(fù)荷數(shù)、系統(tǒng)可切負(fù)荷數(shù)。對(duì)所有測(cè)試系統(tǒng)，在t=0時(shí)刻發(fā)生三相短路故障，t=0.2 s切除故障線路，t=0.3 s執(zhí)行切機(jī)切負(fù)荷操作。tsim為算法仿真時(shí)間，即為切負(fù)荷時(shí)刻至末端時(shí)刻的時(shí)間差。采用隱式梯形積分法進(jìn)行暫態(tài)穩(wěn)定分析時(shí)，積分步長(zhǎng)在本文中取為Δt。程序中各發(fā)電機(jī)相對(duì)于慣性中心的偏差不超過(guò)±110°作為暫態(tài)穩(wěn)定判據(jù)。

表1 測(cè)試算例參數(shù)Table 1 Parameters for case tests

表2給出了4個(gè)測(cè)試算例采用隱式梯形積分差分化方法離散微分方程和網(wǎng)絡(luò)方程后問(wèn)題規(guī)模，其中n表示問(wèn)題狀態(tài)變量數(shù)，m表示問(wèn)題等式約束數(shù)，δDOF表示問(wèn)題自由度即為n與m之差，rDIM、NNZ分別表示原對(duì)偶系統(tǒng)式（12）的維數(shù)和非零元數(shù)量?？梢钥闯?，差分化后優(yōu)化問(wèn)題自由度很小，即滿足，從而該問(wèn)題非常適合簡(jiǎn)約空間下求解。

表2 測(cè)試算例問(wèn)題規(guī)模Table 2 Scales of case tests

本文算法采用C++編程實(shí)現(xiàn)，測(cè)試環(huán)境為Intel Xeon 2 quad-core CPU，OpenMP被用來(lái)實(shí)現(xiàn)多線程并行。通過(guò)調(diào)用英特爾提供的高性能多線程線性代數(shù)庫(kù)（BLAS和LAPACK）實(shí)現(xiàn)算法中的一些矩陣運(yùn)算。KLU解法器被用來(lái)對(duì)算法中稀疏線性方程組進(jìn)行求解。

3.2 測(cè)試結(jié)果

本文算法是基于OpenMP并行簡(jiǎn)約空間內(nèi)點(diǎn)法。為表述簡(jiǎn)潔，本文稱傳統(tǒng)內(nèi)點(diǎn)算法為全空間內(nèi)點(diǎn)算法。為了比較本文算法與常規(guī)算法的性能優(yōu)劣，對(duì)上面4個(gè)算例用隱式梯形差分化方法離散微分方程和網(wǎng)絡(luò)方程，并在全空間、簡(jiǎn)約空間以及并行簡(jiǎn)約空間下求解，測(cè)試結(jié)果如表3所示。表3中“×”表示因時(shí)間或者內(nèi)存不足導(dǎo)致問(wèn)題不可解。

表3 測(cè)試結(jié)果Table 3 Results of case tests

從表3中可以看出，采用隱式梯形積分差分微分方程和網(wǎng)絡(luò)方程，簡(jiǎn)約空間下計(jì)算速度明顯比全空間下計(jì)算速度要快，并且可以求解更大規(guī)模系統(tǒng)緊急控制問(wèn)題。因簡(jiǎn)約空間技術(shù)只是在求取式（12）時(shí)采取的策略，其并不改變迭代過(guò)程，簡(jiǎn)約空間和全空間下問(wèn)題迭代進(jìn)程是完全一樣的，得出的優(yōu)化結(jié)果也一樣。因此，本文將把全空間和簡(jiǎn)約空間得出問(wèn)題的迭代次數(shù)統(tǒng)一列出。

基于OpenMP并行簡(jiǎn)約空間算法是在簡(jiǎn)約空間算法的框架下對(duì)一些計(jì)算步驟利用多線程技術(shù)并行求解，它不改變算法的計(jì)算進(jìn)程，因而兩者有相同的迭代進(jìn)程和相同的優(yōu)化結(jié)果。從表3可以看出，在8核處理器計(jì)算平臺(tái)下，并行簡(jiǎn)約空間內(nèi)點(diǎn)算法計(jì)算效率進(jìn)一步提高。

本文選取切機(jī)切負(fù)荷作為緊急控制措施來(lái)保證暫態(tài)穩(wěn)定性。以CASE300為例，列出優(yōu)化后切機(jī)、切負(fù)荷量如表4所示。

表4 切機(jī)切負(fù)荷優(yōu)化結(jié)果Table 4 Results of generator trip and load shedding optimization

3.3 并行性能分析

首先對(duì)串行簡(jiǎn)約空間下算法各部分計(jì)算時(shí)間進(jìn)行分析如表5所示。從表5可以看出，簡(jiǎn)約空間內(nèi)點(diǎn)算法下暫態(tài)穩(wěn)定緊急控制問(wèn)題的計(jì)算時(shí)間主要集中在B計(jì)算時(shí)間和Zm回代時(shí)間。因而本文對(duì)這部分進(jìn)行多線程并行求解，提高計(jì)算效率。

表5 串行簡(jiǎn)約空間內(nèi)點(diǎn)算法計(jì)算效率分析Table 5 Calculation efficiency analysis of serial reduced-space interior point method s

圖3 Zm回代過(guò)程加速比Fig.3 Acceleration ratio of Zmback substitution

圖4 B計(jì)算過(guò)程加速比Fig.4 Acceleration ratio of B calculation

為了分析B計(jì)算和Zm回代的并行效率，線程數(shù)p取1到8分別進(jìn)行測(cè)試，測(cè)試結(jié)果如圖3、圖4所示。從圖3、圖4可以看出，當(dāng)線程數(shù)較少時(shí)，本文算法幾乎可以達(dá)到線性加速比，隨著線程數(shù)的增加，加速比增長(zhǎng)放緩，在8核處理器下，加速比最大可達(dá)4.5。一般而言，測(cè)試算例規(guī)模越大，本文算法可獲得越好的加速比。因此，本文算法對(duì)大規(guī)模系統(tǒng)暫態(tài)穩(wěn)定緊急控制問(wèn)題的求解有較好的潛力。

3.4 時(shí)域仿真驗(yàn)證

通過(guò)時(shí)域仿真來(lái)驗(yàn)證上面4個(gè)算例的第一擺穩(wěn)定性，圖5給出了各個(gè)測(cè)試算例發(fā)電機(jī)搖擺曲線，圖中實(shí)線為系統(tǒng)發(fā)生故障后不執(zhí)行緊急控制操作的功角曲線；虛線為執(zhí)行切機(jī)切負(fù)荷操作的功角曲線。因全空間內(nèi)點(diǎn)算法、簡(jiǎn)約空間內(nèi)點(diǎn)算法和并行簡(jiǎn)約空間內(nèi)點(diǎn)算法三者有相同的迭代進(jìn)程，因此3種算法得出的優(yōu)化曲線一致，本文只給出簡(jiǎn)約空間下發(fā)電機(jī)搖擺曲線。由圖5可以看出，電力系統(tǒng)發(fā)生故障后不執(zhí)行切機(jī)切負(fù)荷操作將失去穩(wěn)定。采用本文提出的算法執(zhí)行切機(jī)切負(fù)荷操作，系統(tǒng)能夠保持暫態(tài)穩(wěn)定。

圖5 測(cè)試算例時(shí)域仿真結(jié)果Fig.5 Time-domain simulative results of case tests

4 結(jié)語(yǔ)

電力系統(tǒng)暫態(tài)穩(wěn)定緊急控制問(wèn)題是一類大規(guī)模非線性動(dòng)態(tài)優(yōu)化問(wèn)題。本文提出了一種基于OpenMP的并行簡(jiǎn)約空間內(nèi)點(diǎn)法求解該類問(wèn)題。測(cè)試結(jié)果表明，相對(duì)傳統(tǒng)內(nèi)點(diǎn)算法，本文提出的并行簡(jiǎn)約空間內(nèi)點(diǎn)算法計(jì)算效率明顯提高，能計(jì)算更大規(guī)模系統(tǒng)，在多核處理器平臺(tái)下，多線程技術(shù)進(jìn)一步提高了算法計(jì)算速度。隨著多核處理器的普及，本文算法實(shí)際效益更為明顯。