王 云，江全元

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

0 引言

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

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

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

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

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

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

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

1.1 目標函數(shù)

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

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

1.2 等式約束

a.發(fā)電機動態(tài)方程。

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

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

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

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

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

對于可切機的發(fā)電機節(jié)點，有：

對于不可切機的發(fā)電機節(jié)點，有：

對于可切負荷節(jié)點，有：

對于不可切負荷節(jié)點，有：

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

1.3 不等式約束

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

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

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

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

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

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

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

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

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

2.2 OpenMP并行算法

2.2.1 OpenMP編程模型

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

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

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

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

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

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

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

圖2 Zm多線程計算Fig.2 Multi-thread calculation of Zm

b.矩陣B的多線程計算。

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

2.3 算法實施流程

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

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

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

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

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

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

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

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

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

3 算例分析

3.1 測試算例

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

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

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

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

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

3.2 測試結(jié)果

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

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

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

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

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

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

3.3 并行性能分析

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

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

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

圖4 B計算過程加速比Fig.4 Acceleration ratio of B calculation

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

3.4 時域仿真驗證

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

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

4 結(jié)語

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