馮軍基，耿光超，王 云，江全元

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

0 引言

隨著電力市場化改革的深入和電網(wǎng)規(guī)模的擴大，電力系統(tǒng)在遭受大擾動或突發(fā)故障下的安全穩(wěn)定問題日益突出[1]。為保證電力系統(tǒng)在故障切除后仍能穩(wěn)定運行，需要及時對系統(tǒng)進行緊急控制，切機、切負(fù)荷是應(yīng)用最廣的2種措施。

暫態(tài)穩(wěn)定約束緊急控制TSCEC（Transient Stability Constrained Emergency Control）問題的解決方法主要有3類：第一類是基于暫態(tài)能量函數(shù)的直接法[2-4]；第二類是試湊法和啟發(fā)式算法[5-6]；第三類方法將原問題描述為包含微分代數(shù)方程組（DAE）約束的最優(yōu)控制問題OCP（Optimal Control Problem）進行求解[7-9]。

文獻[2]采用基于最大勢能搜索的勢能邊界面法（PEBS）計算臨界勢能，但計算結(jié)果存在可靠性問題。文獻[3]采用主導(dǎo)不穩(wěn)定平衡點法（BCU）制定緊急控制策略，有可能因為失穩(wěn)模式判斷錯誤引起誤差，計算結(jié)果偏于保守。文獻[4]基于擴展等面積法則（EEAC）提出了一種動態(tài)安全域方法，但是受等效模型限制，不適用于非線性時變因素較強的大系統(tǒng)。

文獻[5]和文獻[6]分別使用神經(jīng)網(wǎng)絡(luò)算法和蟻群算法進行切機、切負(fù)荷計算，啟發(fā)式算法直觀簡單，易于獲得可行解，但花費的計算時間較多。

文獻[7]使用約束轉(zhuǎn)化方法和控制向量參數(shù)化方法求解最優(yōu)控制問題，將TSCEC分解為靈敏度計算和非線性規(guī)劃（NLP）問題交替求解。文獻[8]和文獻[9]分別使用隱式梯形積分算法和有限元正交配置 OCFE（Orthogonal Collocation on Finite Elements）算法將暫態(tài)穩(wěn)定約束離散化后，對生成的NLP問題進行求解。

基于最優(yōu)控制理論的TSCEC求解算法首先對暫態(tài)穩(wěn)定約束進行處理，然后利用非線性優(yōu)化算法進行迭代計算。OCP的求解策略包括間接法（變分法）和直接法，后者具有較好的收斂性，能方便地處理不等式約束，成為最優(yōu)控制的主要研究方向[10]。根據(jù)對DAE的處理方式，直接法又分為序貫法和聯(lián)立法。序貫法包括控制向量參數(shù)化和多重打靶法；聯(lián)立法包括隱式梯形積分、有限元正交配置和偽譜法PS（Pseudo-Spectral method）等方法。近年來，偽譜法作為一種高階方法，得到了迅速的發(fā)展，并在航天、化工等領(lǐng)域廣泛應(yīng)用[11-14]。

在使用低階的隱式梯形積分等聯(lián)立法求解TSCEC問題時，會由大量的離散點生成高維系統(tǒng)，導(dǎo)致計算困難。本文將高階的偽譜法引入TSCEC問題，以偽譜法配置離散點，采用內(nèi)點法進行求解，并在不同規(guī)模系統(tǒng)上與有限元正交配置方法進行對比，該方法具有以下特點：數(shù)值解收斂半徑大，對初值猜測不敏感，魯棒性強；計算精度高，可以用較少的離散點數(shù)量獲得較精確的解，降低NLP系統(tǒng)規(guī)模；通過細(xì)粒度并行可以大幅提高計算性能。

1 暫態(tài)穩(wěn)定緊急控制問題的建模

暫態(tài)穩(wěn)定約束的緊急控制問題可以描述成一個DAE約束的最優(yōu)控制問題：以切機與切負(fù)荷量最小為目標(biāo)，使系統(tǒng)運行滿足暫態(tài)穩(wěn)定約束。

其中，u=[uG，uL]為切機、切負(fù)荷比例的控制向量，其為非時變量；cG、cL與PG、PL分別為可切負(fù)荷與可切發(fā)電機節(jié)點的權(quán)值與功率向量，cGi、cLi、PGi、PLi分別為對應(yīng)向量的第i個元素；x為系統(tǒng)的狀態(tài)變量，y為系統(tǒng)的代數(shù)變量，x、y均為時變量。

1.1 等式約束

TSCEC的等式約束是一組DAE方程組。微分方程組 x˙=g（x，y，u）包含發(fā)電機轉(zhuǎn)子運動方程、勵磁系統(tǒng)與調(diào)整系統(tǒng)動態(tài)方程；代數(shù)方程組 0=φ（x，y，u）主要為電力網(wǎng)絡(luò)方程。

第i臺發(fā)電機的轉(zhuǎn)子運動方程可以描述為：

其中，NG為發(fā)電機數(shù)目；δi為發(fā)電機功角；ωi為發(fā)電機角速度；PGi為發(fā)電機有功功率；Pei為發(fā)電機電磁功率；TJi為發(fā)電機慣性時間常數(shù)。

電力網(wǎng)絡(luò)方程YU=I可寫為：

其中，Y為節(jié)點導(dǎo)納矩陣；Yij=Gij+jBij為節(jié)點i和j的互導(dǎo)納；Yii=Gii+jBii為節(jié)點 i的自導(dǎo)納；Uxi、Uyi分別為節(jié)點i電壓的實部、虛部；Ixi、Iyi分別為注入電流的實部、虛部，非發(fā)電機節(jié)點注入電流為0。

對于有可切發(fā)電機的節(jié)點i，假設(shè)發(fā)電機等值導(dǎo)納為YGi=GGi+jBGi，則切除后節(jié)點i的自導(dǎo)納為G′ii=Gii-uGiGGi、B′ii=Bii-uGiBGi，注入電流為 I′xi= （1-uGi）Ixi、I′yi= （1-uGi）Iyi。

對于有可切負(fù)荷的節(jié)點i，假設(shè)負(fù)荷等值導(dǎo)納為YLi=GLi+jBLi，則切除后節(jié)點i的自導(dǎo)納為G′ii=GiiuLiGLi、B′ii=Bii-uLiBLi。

1.2 不等式約束

不等式約束包括暫態(tài)穩(wěn)定約束和控制向量的上下限約束。本文用各發(fā)電機功角與慣性中心角的偏差不超過指定角度作為暫態(tài)穩(wěn)定判據(jù)。不等式可以用式（4）表示。

2 偽譜法原理

偽譜法是一種全局正交配置方法，它通過對狀態(tài)變量和代數(shù)變量進行離散，將含有時變變量的約束轉(zhuǎn)化為非時變的代數(shù)等式/不等式約束，使含DAE約束的最優(yōu)控制問題轉(zhuǎn)化為普通的NLP問題。

2.1 全局配置法

全局法在整個問題區(qū)間上通過一系列多項式逼近狀態(tài)變量。偽譜法通常選取三角多項式、勒讓德多項式等正交多項式為基函數(shù)，在正交配置點上離散變量，用有限序列之和逼近真實解。

假設(shè)在優(yōu)化區(qū)間[t0，tf]上配置有 R 個點[t1，t2，…，tR]，狀態(tài)變量 x（t）和代數(shù)變量 y（t）在 ti點的離散值分別為 Xi、Yi，為求取 x˙（t）在離散點的數(shù)值，將 x（t）用拉格朗日多項式表示如下：

對式（5）進行微分得到：

其中，矩陣 D=｛Dki｝為微分矩陣。

2.2 配置點選取

如果用等距的配置點進行多項式插值，一般可以通過配置更多的點提高逼近程度，但在一些情況下，點數(shù)增多會導(dǎo)致插值區(qū)間邊緣產(chǎn)生無窮大誤差，稱為龍格（Runge）現(xiàn)象[15]。選用正交多項式的根作為配置點是消除龍格現(xiàn)象的一種有效方法。

勒讓德正交多項式滿足遞推關(guān)系[16]和正交關(guān)系分別如式（7）、（8）所示。

其中，δkj為 Kronecker函數(shù)。

常用如下3類配置點。

a.Legendre-Gauss點：選用 LN+1（v）=0 的根 vj（0≤j≤N）。

b.Legendre-Gauss-Radau 點：選用 LN（v）+LN+1（v）=0 的根 vj（0≤j≤N）。

c.Legendre-Gauss-Lobatto 點：v0=-1，vN=1，｛vj｝N-1j=1為 L′N（v）=0 的根。

在[-1，1]區(qū)間上用3類方法進行5節(jié)點配置的點分布如圖1所示，Legendre-Gauss法不包含優(yōu)化區(qū)間端點，Radau法僅包含1個端點，Lobatto法包含2個端點。

圖1 正交配置點分布Fig.1 Distribution of orthogonal collocation points

理論上，Legendre-Gauss-Radau點在處理剛性問題時具有良好的穩(wěn)定性[17]，最優(yōu)化軟件包GPOPS-Ⅱ就以該方法為基礎(chǔ)實現(xiàn)離散點配置[14]。本文也選用Radau方法進行最優(yōu)控制。

3 算法實現(xiàn)

3.1 原問題到NLP的轉(zhuǎn)化

設(shè)區(qū)間[-1，1]上的 Legendre-Gauss-Radau 配置點為[τ0，τ1，…，τM]，其對應(yīng)的標(biāo)準(zhǔn)化微分矩陣為：

其中，k=0，1，…，M-1；i=0，1，…，M。

將變量τ∈[-1，1]變換到優(yōu)化區(qū)間 t∈[t0，tf]：

區(qū)間[t0，tf]對應(yīng)的離散點為[t0，t1，…，tN]，微分矩陣為：

在 t=tk處等式約束 x˙=g（x，y，u）的離散化形式為：

其中，Yk是代數(shù)狀態(tài)變量在t=tk處的值。

DAE約束的緊急控制問題可以轉(zhuǎn)化為下面的NLP問題：

其中，Δk?Xk、ΔCOIk?Yk（k=0，1，…，M）。 優(yōu)化變量為Xk、Yk、u=[uG，uL]。

3.2 NLP求解

通過離散生成的NLP問題是一個不含微分方程的大規(guī)模稀疏問題。常用的NLP求解算法有序列二次規(guī)劃 SQP（Sequential Quadratic Programming）法、內(nèi)點法 IPM（Interior Point Method）等。

在求解時，需要提供初始點、雅可比矩陣和海森矩陣信息。

a.初始點：在控制時刻前進行暫態(tài)仿真，選取u（0）=｛u0G，u0L｝進行全時段仿真，將離散點處的仿真值作為相應(yīng)優(yōu)化變量的初值。

b.雅可比矩陣：將各約束分為線性部分和非線性部分，線性部分雅可比矩陣只需求解一次，由當(dāng)前優(yōu)化變量和約束解析表達式求出非線性部分雅可比矩陣后疊加[18-19]。離散生成的雅可比矩陣結(jié)構(gòu)如圖2所示，矩陣非零元個數(shù)為55916。

圖2 雅可比矩陣結(jié)構(gòu)Fig.2 Structure of Jacobi matrix

c.海森矩陣：采用擬牛頓法近似海森矩陣。求取算法采用基于有限內(nèi)存的L-BFGS（Limited-memory BFGS）算法[20]。

4 算例測試

所用測試算例見表1，模擬故障為三相接地，仿真時長為1.8s，計算平臺為Intel Xeon X5650 2.67GHzPC機。其中Case 118和Case 300分別由IEEE118和IEEE300標(biāo)準(zhǔn)測試系統(tǒng)修改得到；Case 678、Case 2 052和Case 2 383分別由2007年意大利、2005年日本和2000年波蘭的實際電力系統(tǒng)簡化得到。

表1 測試算例Table 1 Cases for test

在計算時，有限元正交配置法配置3段共14個點，偽譜法使用14個點。為達到理想的精度，隱式梯形積分使用40個差分點，步長Δt=0.0375 s。

IPOPT是一個開源的優(yōu)化軟件包，使用C++編寫，通過內(nèi)點法求解大規(guī)模非線性優(yōu)化問題[21]。

本文用IPOPT求解產(chǎn)生的NLP問題，為提高性能，利用Intel MKL庫重新編譯IPOPT。稀疏線性解法器使用 PARDISO[22]，通過OpenMP并行化框架在多處理器主機上實現(xiàn)細(xì)粒度并行計算[23]。

4.1 有效性驗證

在單線程下，分別使用隱式梯形積分、有限元正交配置和偽譜法對各算例進行計算，結(jié)果如表2所示。表中niter表示迭代次數(shù)，tCPU表示總計算時間，POBJ表示切除的總功率，“Fail”表示算例不收斂。

表2 不同離散方法下的優(yōu)化結(jié)果Table 2 Optimization results for different discrete methods

優(yōu)化完成后，根據(jù)計算出的切機切負(fù)荷量，使用4階龍格-庫塔算法（步長Δt=0.01 s）進行時域仿真，驗證優(yōu)化結(jié)果，均滿足暫態(tài)穩(wěn)定條件。圖3是對118節(jié)點算例進行的仿真驗證（點劃線表示未施加控制時暫態(tài)不穩(wěn)定的情況，實線表示優(yōu)化后切機切負(fù)荷滿足暫態(tài)穩(wěn)定約束）。

圖3 Case118節(jié)點算例仿真驗證Fig.3 Simulative verification by Case 118-bus system

在Case 118算例測試中，通過隨機函數(shù)在約束范圍內(nèi)生成隨機數(shù)作為控制參數(shù) u（0）= [u0G，u0L]的初值，偽譜法求解均能收斂。表3給出了選取9組不同初值時，3種算法的迭代次數(shù)?？梢钥闯?，偽譜法對初值的要求不高，隱式梯形積分方法使用表中4組初值計算時無法收斂。在仿真驗證中，偽譜法對初值的要求最低，具有很好的魯棒性和收斂性。

表3 不同初值下的迭代次數(shù)Table 3 Iteration times for different initial values

4.2 偽譜法與局部方法的對比

利用優(yōu)化計算出的切機切負(fù)荷量和4階龍格-庫塔算法進行時域仿真生成的狀態(tài)變量記為x～（t），優(yōu)化后各離散點對應(yīng)的狀態(tài)變量記為X*k，定義狀態(tài)變量誤差矩陣E，其第k列表示為：

其中，tk為使用的正交配置點。取E各行的最大值（即各狀態(tài)變量在配置點的最大誤差）為誤差向量vE=｛vi｝iN=s1（Ns為狀態(tài)變量總數(shù)）。 記狀態(tài)變量誤差：

采用40個差分節(jié)點時，隱式梯形積分方法的狀態(tài)變量誤差見表4。

表2中的正交配置方法均使用14個離散點，與隱式梯形積分方法相比，離散生成的NLP問題維數(shù)大幅降低，在計算時間上有明顯的優(yōu)勢。此時，雖然有限元正交配置法和偽譜法使用了更少的離散點，但其誤差精度仍高于隱式梯形積分方法（參考表4中配置點數(shù)為14的誤差結(jié)果）。2種正交配置方法均遠(yuǎn)優(yōu)于隱式梯形積分方法。

分別在配置點數(shù)為10、12、14、16下計算有限元正交配置法與偽譜法的誤差，結(jié)果見表4（有限元正交配置法使用3個分段）。

在相同的配置點數(shù)下，有限元正交配置法和偽譜法的迭代次數(shù)和計算時間都是相當(dāng)?shù)模捎趥巫V法的微分矩陣為稠密矩陣，單次迭代花費的時間略大于有限元正交配置方法。

從誤差上看，偽譜法的精度略優(yōu)于有限元正交配置法，這是因為偽譜法的離散化方程包含了整個區(qū)間上的配置點信息，而有限元正交配置法僅包含1個分段內(nèi)的信息。

表4 配置點與誤差Table 4 Collocation points and errors

4.3 并行加速

利用內(nèi)點法求解稀疏大規(guī)模NLP問題時，大部分時間在求解原對偶系統(tǒng)線性方程組，對線性方程組求解實現(xiàn)細(xì)粒度并行可以大幅縮短計算時間。

PARDISO線性解法器可在共享內(nèi)存SMP平臺實現(xiàn)OpenMP并行加速。分別在線程數(shù)為1、2、4、8下統(tǒng)計偽譜法的優(yōu)化運行時間，結(jié)果如表5所示。

表5 不同線程下計算時間Table 5 Computation time for different thread numbers

并發(fā)線程數(shù)的改變主要影響線性方程求解效率，對算法其他部分影響很小。圖4表示了線程數(shù)與線性方程求解時間占總時間百分比的關(guān)系，單線程下線性方程求解占用了大部分求解時間。隨著線程數(shù)加倍，線程方程組求解時間所占的比例下降，雖然可以通過并行縮短計算時間，但通過增加線程獲得的總體效率提升空間有限。線程數(shù)與線性方程求解加速比的關(guān)系見圖5，當(dāng)算例規(guī)模擴大時，加速比未表現(xiàn)出單調(diào)性的改變，加速比與算例規(guī)模大小的關(guān)系不明顯。

圖4 線程數(shù)與線性方程求解占總時間百分比的關(guān)系Fig.4 Relationship between thread number and time percentage for solving linear equations

圖5 線程數(shù)與線性方程求解加速比的關(guān)系Fig.5 Relationship between thread number and speedup ratio for solving linear equations

5 結(jié)語

本文使用求解最優(yōu)控制問題的偽譜法計算電力系統(tǒng)緊急控制問題，并與局部聯(lián)立法進行比較。不同規(guī)模電力系統(tǒng)的測試說明該算法具有以下特點。

a.在聯(lián)立法中，偽譜法的計算性能和精度顯著高于隱式梯形積分方法；在相同的配置點數(shù)下，偽譜法與有限元正交配置方法計算時間相當(dāng)，而誤差略小于有限元正交配置方法。偽譜法的應(yīng)用為該問題求解帶來了精度上的提升。

b.對于大規(guī)模系統(tǒng)，線性方程組求解成為計算性能瓶頸。利用OpenMP在共享內(nèi)存平臺上實現(xiàn)細(xì)粒度并行，可以大幅提高問題的求解效率。