電力系統(tǒng)暫態(tài)穩(wěn)定緊急控制并行序貫優(yōu)化算法

甘國曉，耿光超，仲悟之，江全元，徐振華，蘇 毅，黃 霆

(1. 浙江大學 電氣工程學院，浙江 杭州 310027；2. 中國電力科學研究院有限公司，北京 100192；3. 國網(wǎng)福建省電力有限公司電力科學研究院，福建 福州 350007)

0 引言

電力系統(tǒng)遭受嚴重擾動后，必須盡快采取控制措施保證系統(tǒng)的安全穩(wěn)定運行。切機切負荷作為一種有效的暫態(tài)穩(wěn)定緊急控制手段已得到廣泛認可。然而，隨著電力系統(tǒng)以及電力市場化的發(fā)展，電力系統(tǒng)緊急控制問題變得更加復雜。

目前，對于暫態(tài)穩(wěn)定緊急控制問題的研究方法主要分為兩大類：基于能量函數(shù)的直接法和基于最優(yōu)控制理論的優(yōu)化算法。直接法主要有基于穩(wěn)定域邊界的主導不穩(wěn)定平衡點法[1]、暫態(tài)能量函數(shù)法[2]和擴展等面積法[3-4]等。它們都是構(gòu)造一個穩(wěn)定裕度，根據(jù)穩(wěn)定裕度相對控制變量的靈敏度計算滿足要求的控制策略。此類方法由于模型適應性較差，存在計算結(jié)果不準確的問題。基于最優(yōu)控制理論的方法將緊急控制問題描述為包含微分代數(shù)方程(DAE)的最優(yōu)控制問題[5]，采用動態(tài)優(yōu)化算法求解。此類方法具有廣泛的模型適應性，但存在求解效率不高的缺點。

求解動態(tài)優(yōu)化問題的常用方法有聯(lián)立法[6-8]和序貫類[9-11]的方法。聯(lián)立法通過將DAE離散為非線性代數(shù)方程組，從而使動態(tài)優(yōu)化問題轉(zhuǎn)化為大規(guī)模非線性規(guī)劃(NLP)問題，使用非線性優(yōu)化方法直接求解。文獻[7]和文獻[8]將緊急控制問題建模為大規(guī)模NLP問題，分別基于并行簡約空間內(nèi)點法和偽普法進行求解，在一定程度上提升了求解效率。但隨著系統(tǒng)規(guī)模的增大，“維數(shù)災”的問題不可避免，計算速度仍然較慢。序貫類的方法則將DAE與優(yōu)化問題分開求解。文獻[9-10]基于約束轉(zhuǎn)換，將最優(yōu)切負荷問題分解為靈敏度計算和NLP問題的交替求解。文獻[11]采用序貫法框架將動態(tài)優(yōu)化問題分為模擬層和優(yōu)化層，在模擬層基于正交配置求解DAE，但是靈敏度分析的計算量大，耗時較長。序貫法將DAE剝離原優(yōu)化問題進行求解，雖然形成了較小規(guī)模的NLP問題，但同時也帶來了DAE求解和靈敏度分析耗時過大的缺點，這在大規(guī)模系統(tǒng)的動態(tài)優(yōu)化中尤為突出。

本文基于序貫法的框架將原問題分為仿真層和優(yōu)化層進行交替求解。為提高仿真層的計算效率，基于有限元正交配置，并行求解DAE與靈敏度。在優(yōu)化層采用預測校正內(nèi)點法求解較小規(guī)模的NLP問題。不同方法的計算結(jié)果表明，有限元正交配置能以較少的離散點數(shù)獲得較精確的結(jié)果，減少了計算量；并行計算技術(shù)加速了仿真層的計算過程；另外，由于不需要直接求解大規(guī)模NLP問題，本文方法比聯(lián)立法的計算速度更快。

1 電力系統(tǒng)緊急控制數(shù)學模型

暫態(tài)穩(wěn)定緊急控制問題可描述成一個含有DAE約束的動態(tài)優(yōu)化問題，其一般形式為：

(1)

1.1 目標函數(shù)

為保持暫態(tài)穩(wěn)定性，以切除最少的發(fā)電機和負荷為目標：

(2)

其中，u=[uGi,uL i]，uGi、uL i分別為各個發(fā)電機節(jié)點切機比例和負荷節(jié)點切負荷比例；rG、rL分別為可切發(fā)電機和可切負荷總數(shù)；aGi、aL i分別為切機、切負荷控制代價的權(quán)值；PGi、PL i分別為各可切發(fā)電機節(jié)點的容量和各可切負荷節(jié)點的總負荷量。

1.2 等式約束

a. 動態(tài)方程。

電力系統(tǒng)動態(tài)方程包含發(fā)電機轉(zhuǎn)子運動方程、勵磁系統(tǒng)與調(diào)速系統(tǒng)動態(tài)方程等，可以描述為：

(3)

其中，x為包括發(fā)電機功角、角速度等在內(nèi)的微分狀態(tài)變量；V為包括節(jié)點電壓等在內(nèi)的代數(shù)變量。

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

網(wǎng)絡方程可以表示成如下形式：

(4)

其中，G、B分別為節(jié)點導納矩陣的電導部分和電納部分；Vx、Vy分別為節(jié)點電壓實部和虛部；Ix、Iy分別為節(jié)點注入電流的實部和虛部。

1.3 不等式約束

不等式約束包括暫態(tài)穩(wěn)定約束和切機切負荷比例的上下限約束。本文采用各發(fā)電機功角與慣性中心角的偏差不超過指定角度作為暫態(tài)穩(wěn)定約束：

(5)

其中，δi為第i臺發(fā)電機功角；NG為發(fā)電機總數(shù)；δCOI為發(fā)電機功角的慣性中心，可由式(6)計算。

(6)

其中，TJi為第i臺發(fā)電機的慣性時間常數(shù)。

切機切負荷上下限約束：

(7)

2 并行序貫優(yōu)化算法

2.1 序貫優(yōu)化算法

序貫法將原動態(tài)優(yōu)化問題分為仿真和優(yōu)化2個子問題交替迭代，構(gòu)成一個序貫的求解過程，算法框架如圖1所示。仿真子問題求解式(1)中的DAE初值問題：

圖1 本文算法框架Fig.1 Framework of proposed method

(8)

給定狀態(tài)變量的初始值x0，在仿真層通過數(shù)值方法求解動態(tài)優(yōu)化問題式(1)中的DAE約束Ψ，可以得到狀態(tài)變量x在整個時域上的值。求解DAE的常用方法有改進歐拉法、隱式梯形積分等方法，在基于聯(lián)立法[7]或者序貫法[12]的暫態(tài)穩(wěn)定控制優(yōu)化算法中，常用隱式梯形積分對微分方程進行離散化。隱式梯形積分求解DAE時采用較小的步長，計算量大。而有限元正交配置[11,13]具有離散點數(shù)少、計算精度高、穩(wěn)定性好等優(yōu)點，為了提高仿真層的計算效率，本文用其求解DAE。

在仿真層采用數(shù)值方法求出狀態(tài)變量后，可以得到各個狀態(tài)變量在各個離散時間點上的數(shù)值解，再將狀態(tài)變量的不等式約束加于所有離散點上，可使過程約束H近似為狀態(tài)變量在各個離散點處的約束。此時DAE等式約束已經(jīng)被消除，狀態(tài)變量x可以視為控制變量u的隱函數(shù)，仍需計算狀態(tài)變量對控制變量的靈敏度信息。則原動態(tài)優(yōu)化問題可以轉(zhuǎn)化為一個較小規(guī)模的NLP問題：

(9)

2.2 基于有限元正交配置的DAE求解

如圖2所示，在時域t0～T上配置NL個有限元，每個有限元有NC+1個配置點。為保持狀態(tài)變量的連續(xù)性，本文采用Radau 配置點法，即前一有限元的最后一個配置點等于后一有限元的起始配置點，以NC=3為例，如圖2所示。式(3)中的微分變量x可近似為一組拉格朗日多項式的線性組合：

i=1,2，…,NC；n=1,2，…,NL

(10)

其中，xn(tn,i)為狀態(tài)變量在第n個有限元的多項式近似；xn,j為狀態(tài)變量在第n個有限元第j個配置點的值；lj(tn,i)為拉格朗日多項式；tn,i為第n個有限元的第i個配置點。

圖2 Radau正交配置點Fig.2 Radau orthogonal collocation points

對式(10)的拉格朗日多項式求一階導可得：

i=1,2，…,NC；n=1,2，…,NL

(11)

通過正交配置法，式(3)的微分方程可以轉(zhuǎn)化為一系列的非線性代數(shù)方程：

i=1,2，…,NC；n=1,2，…,NL

(12)

求解上述非線性代數(shù)方程組即可得到狀態(tài)變量在各個離散點處的值。本文采用非誠實牛頓法(VDHN)[14]求解正交配置離散后的DAE。

2.3 并行靈敏度計算

在第n個有限元，離散后的DAE可以表示為：

(13)

其中，F(xiàn)n和gn分別為第n個有限元處離散后的動態(tài)方程和網(wǎng)絡方程；xn為第n個有限元中NC個配置點上的狀態(tài)變量；xn,0為第n個有限元中初始配置點的值；Vn為第n個有限元中NC個配置點上的代數(shù)變量。

為了求出各個離散點處狀態(tài)變量對控制量的靈敏度，對式(13)采用鏈式求導法則可得：

(14)

定義Jn如下：

(15)

第n個有限元的靈敏度方程可以表示為：

(16)

由于采用Radau 配置點法，即前一有限元的最后一個配置點等于后一有限元的起始配置點，則有：

(17)

通過式(16)、(17)，狀態(tài)變量對控制變量的靈敏度可以從一個有限元到下一個有限元順序進行求解，這與DAE的求解順序是一致的，因此可以將它們并行執(zhí)行。值得注意的是，在求解第n個有限元處的DAE時，第k次牛頓迭代需要求解線性方程組：

(18)

圖3 仿真層框架Fig.3 Framework of simulation layer

另外，靈敏度方程式(16)是一個含多右邊項的線性方程組，利用求解DAE時Jn的因式分解結(jié)果，只需進行前代回代過程即可得到第n個有限元處的靈敏度。由于多個右邊項的前代回代過程相互獨立，為了進一步提升效率，可以對不同右邊項的回代過程采用多個計算核心并行執(zhí)行，如圖4所示。

圖4 并行靈敏度計算Fig.4 Parallel sensitivity calculation

計算出狀態(tài)變量x對于控制變量u的靈敏度即可求出不等式約束對控制量的一階梯度。不僅一階梯度，二階梯度對基于梯度的優(yōu)化算法也至關(guān)重要。由于對控制變量的二階導數(shù)難以用解析的方法獲得，BFGS(Broyden-Fletcher-Goldfarb-Shanno)[15]方法被用于近似海森矩陣：

(19)

(20)

其中，un為第n次迭代的控制變量值；yn和ζn如式(21)所示。

(21)

其中，▽uΦn為第n次迭代時目標函數(shù)對控制量的一階梯度；▽uHn為第n次迭代時不等式約束對控制量的一階梯度；zH和wH為內(nèi)點法中不等式約束對應的拉格朗日乘子。

3 算例分析

本文算法采用C++ 編程實現(xiàn)，測試環(huán)境為Intel Xeon E5-2650 2 GHz CPU，采用OpenMP[16]實現(xiàn)多線程并行。其中，仿真過程采用了KLU[17]解法器進行稀疏線性方程組的求解，優(yōu)化過程調(diào)用了英特爾提供的高性能多線程線性代數(shù)庫(MKL)進行稠密矩陣的運算。測試算例參數(shù)如表1所示，其中nB、nL、nG、nLoad分別為系統(tǒng)節(jié)點數(shù)、線路數(shù)、發(fā)電機數(shù)和負荷數(shù)。發(fā)電機考慮二階經(jīng)典模型，負荷采用恒定阻抗模型，對所有算例優(yōu)化1.8 s，t=0時刻發(fā)生三相接地故障，t=0.2 s切除故障線路，t=0.3 s實施緊急控制措施。有限元個數(shù)NL=8，NC=3。暫態(tài)穩(wěn)定約束中相對功角的上下限取為±110°。

表1 測試算例Table 1 Test cases

3.1 計算結(jié)果比較

采用不同的方法對算例進行計算，結(jié)果對比如表2所示。對比方法為基于隱式梯形積分的并行序貫法，即在本文算法框架下，在求解DAE時，將有限元正交配置法換成隱式梯形積分，離散步長為0.01 s。從表2可以看出，2種方法的優(yōu)化結(jié)果、迭代次數(shù)相差不大，但在同樣采用并行計算的情況下本文方法耗時更少。與基于隱式梯形積分的并行序貫法對比，本文采用有限元正交配置法離散DAE，能夠以較少的離散點數(shù)達到較好的計算精度，計算效率更高。

表2 計算結(jié)果比較Table 2 Comparison of calculation results

3.2 計算效率分析

表3給出了本文方法的各部分計算耗時，其中tCPU、tS、tO和PS分別為計算總耗時、仿真層耗時、優(yōu)化層耗時和仿真層占總耗時的比值。從表3可以看出，仿真層的計算時間占了總時間的絕大部分，這是由于原動態(tài)優(yōu)化問題的DAE約束已在仿真層計算，轉(zhuǎn)化得到較小規(guī)模的NLP問題，所以優(yōu)化層的求解耗時很少。另外，本文方法只需幾秒鐘的耗時即可得出解決方案，這得益于仿真與優(yōu)化的分離，避免了直接求解大規(guī)模NLP問題，而且并行計算加速了仿真層的求解。

表3 計算時間Table 3 Calculation time

為了表征算法的并行效率，定義了靈敏度計算加速比和仿真層計算加速比。其中，靈敏度計算加速比為單線程靈敏度計算時間除以多線程靈敏度計算時間，仿真層計算加速比為單線程仿真層計算時間除以多線程仿真層計算時間。圖5給出了在不同計算線程數(shù)下靈敏度計算過程的加速比。由圖5可以看出，靈敏度計算的過程幾乎可以達到線性加速比，加速效果明顯。圖6則給出了整個仿真層的加速比，隨著線程數(shù)的增加，加速比增長放緩，8線程時最大加速比達到4.8。通常，測試算例的規(guī)模越大，加速效果越好。因此，本文算法對大規(guī)模系統(tǒng)暫態(tài)穩(wěn)定緊急控制問題的求解有較好的潛力。

圖5 靈敏度計算過程加速比Fig.5 Speedup ratio in process of sensitivity calculation

圖6 仿真層加速比Fig.6 Speedup ratio in simulation layer

圖7 算例仿真結(jié)果Fig.7 Simulative results of test cases

3.3 仿真驗證

對算例進行時域仿真驗證本文算法的有效性，各個算例的發(fā)電機功角曲線如圖7所示。其中，虛線表示系統(tǒng)發(fā)生故障后沒有緊急控制的功角曲線；實線表示執(zhí)行切機切負荷操作的功角曲線。由圖7可以看出，如果不執(zhí)行緊急控制措施，電力系統(tǒng)發(fā)生故障后將失去穩(wěn)定；采用本文算法的計算結(jié)果進行切機切負荷操作，系統(tǒng)的暫態(tài)穩(wěn)定性得到明顯改善。

4 結(jié)論

本文提出一種基于有限元正交配置的并行序貫法求解電力系統(tǒng)緊急控制問題。計算結(jié)果表明，相比隱式梯形積分，有限元正交配置能以較少的離散點數(shù)獲得較精確的結(jié)果，減少了計算量；靈敏度計算過程重用DAE求解時的雅可比因式分解，并將這2個過程并行化，加速了仿真層的計算，提升了序貫優(yōu)化算法的效率。另外，由于不需要直接求解大規(guī)模NLP問題，本文方法計算速度快，適合快速求解大規(guī)模系統(tǒng)的緊急控制問題。