隱式Taylor級(jí)數(shù)暫態(tài)穩(wěn)定計(jì)算變步長方法初探

鄭煥坤，孫耀芹，常鮮戎

（1.華北電力大學(xué) 電氣與電子工程學(xué)院，保定071003；2.冀北保定電力職業(yè)技術(shù)學(xué)院，保定071000）

暫態(tài)穩(wěn)定算法主要有時(shí)域仿真法和直接法。直接法計(jì)算速度快，但其模型簡(jiǎn)單結(jié)果偏保守。時(shí)域仿真法利用各電氣設(shè)備的數(shù)學(xué)模型，通過求解微分方程組給出變量隨時(shí)間變化的曲線，計(jì)算結(jié)果準(zhǔn)確，但其計(jì)算速度較慢。這兩種方法相輔相成，在電力系統(tǒng)離線分析和在線安全分析中都得到廣泛應(yīng)用。高階Taylor級(jí)數(shù)法是一種優(yōu)秀的時(shí)域仿真法，文獻(xiàn)［1］第一次將Taylor級(jí)數(shù)法運(yùn)用到暫態(tài)穩(wěn)定計(jì)算中。文獻(xiàn)［2～4］提出快速高階Taylor級(jí)數(shù)暫態(tài)穩(wěn)定算法，并做了深入研究。文獻(xiàn)［5，6］提出隱式Taylor級(jí)數(shù)方法，并通過調(diào)諧參數(shù)的合理設(shè)置，進(jìn)一步得出了具有A穩(wěn)定性的高精度隱式Taylor級(jí)數(shù)暫態(tài)穩(wěn)定算法的計(jì)算格式。其中文獻(xiàn)［6］給出具有A穩(wěn)定性的高精度隱式調(diào)諧Taylor級(jí)數(shù)法。

近年來，針對(duì)時(shí)域仿真法計(jì)算速度慢的特點(diǎn)，提出很多改進(jìn)措施。文獻(xiàn)［8］對(duì)電力系統(tǒng)暫態(tài)穩(wěn)定仿真和中長期穩(wěn)定性仿真中的變階、變步長技術(shù)進(jìn)行了深入研究。文獻(xiàn)［9］研究了暫態(tài)穩(wěn)定仿真中的并行計(jì)算方法，以適應(yīng)大規(guī)模電力系統(tǒng)在線實(shí)時(shí)仿真要求。文獻(xiàn)［10］利用高階Taylor級(jí)數(shù)法的特點(diǎn)，研究了快速高階Taylor級(jí)數(shù)暫態(tài)穩(wěn)定仿真方法中步長和階數(shù)的動(dòng)態(tài)控制，提高了算法的計(jì)算速度。

本文研究了隱式Taylor級(jí)數(shù)暫態(tài)穩(wěn)定計(jì)算中步長的動(dòng)態(tài)控制。首先對(duì)算法的計(jì)算量進(jìn)行了分析，指出變步長在節(jié)省計(jì)算時(shí)間上的重要性，然后給出了誤差估計(jì)方法，并在此基礎(chǔ)上實(shí)現(xiàn)了變步長功能。仿真算例的結(jié)果表明：算法在采用定步長時(shí)，為取得較高的計(jì)算精度只能采用小步長，步長的減小直接導(dǎo)致了仿真速度的減慢；而算法在變步長時(shí)，步長的大小根據(jù)計(jì)算精度進(jìn)行自動(dòng)調(diào)整，既保證了計(jì)算精度又提高了計(jì)算速度。

1 算法的計(jì)算量分析

隱式Taylor級(jí)數(shù)暫態(tài)穩(wěn)定算法的計(jì)算量主要集中在以下兩個(gè)方面：①發(fā)電機(jī)各狀態(tài)變量的各階導(dǎo)數(shù)的遞推求?。虎谛纬梢蜃颖聿⒗盟鼘?duì)各階導(dǎo)數(shù)進(jìn)行網(wǎng)絡(luò)方程的解算。接下來針對(duì)這兩個(gè)方面對(duì)算法的計(jì)算量進(jìn)行分析。

在消去負(fù)荷節(jié)點(diǎn)和聯(lián)絡(luò)節(jié)點(diǎn)后，Taylor級(jí)數(shù)法中網(wǎng)絡(luò)方程的導(dǎo)納矩陣收縮至發(fā)電機(jī)節(jié)點(diǎn)，可將其看成是滿陣。為了確切表示每一步求解過程的計(jì)算量，對(duì)于各狀態(tài)量的m階導(dǎo)數(shù)的求取均以需用的乘法次數(shù)來表示。設(shè)一個(gè)系統(tǒng)中包含ng個(gè)發(fā)電機(jī)節(jié)點(diǎn)，nf個(gè)故障節(jié)點(diǎn)，則由文獻(xiàn)［11］可知形成因子表所需的計(jì)算量為

進(jìn)行一次前推、回代求解網(wǎng)絡(luò)方程的計(jì)算量為

發(fā)電機(jī)各狀態(tài)變量、坐標(biāo)轉(zhuǎn)換矩陣及諾頓電流的m階導(dǎo)數(shù)的迭代求取的計(jì)算量為

對(duì)于任意m階導(dǎo)數(shù)，發(fā)電機(jī)各狀態(tài)變量的求解運(yùn)算量為

設(shè)Taylor級(jí)數(shù)法的最高導(dǎo)數(shù)階數(shù)為p＋2，由文獻(xiàn)［2］可知實(shí)際參加求解過程為p階，因此從0到p的各階導(dǎo)數(shù)總計(jì)算量為

則對(duì)于每一步的計(jì)算量為

設(shè)定仿真時(shí)間T內(nèi)的仿真總步數(shù)為k，每一步收斂時(shí)迭代的次數(shù)為dk，則仿真時(shí)間內(nèi)總的計(jì)算量為：

由式（7）可見，仿真時(shí)間內(nèi)總的計(jì)算量是總仿真步數(shù)k和每步迭代次數(shù)dk的函數(shù)。在滿足精度要求的前提下，與定步長相比，變步長可使總的仿真步數(shù)減少并保證每一步的迭代次數(shù)最少，從而可減少算法的計(jì)算量，進(jìn)而提高其計(jì)算速度。

2 算法的變步長控制

解微分方程的數(shù)值解法都采用“步進(jìn)式”。步長對(duì)計(jì)算量和計(jì)算精度有較大影響，隨著步長取值的減小，每一步截?cái)嗾`差會(huì)隨之變小，但步數(shù)相應(yīng)增加。步數(shù)的增加，不但引起計(jì)算量的增大，而且可能導(dǎo)致舍入誤差的嚴(yán)重積累。變步長方法能根據(jù)計(jì)算情況對(duì)步長的大小進(jìn)行自動(dòng)調(diào)整，一般可起到既保證計(jì)算精度又節(jié)省計(jì)算量和計(jì)算時(shí)間的作用。

步長的自動(dòng)調(diào)整需要考慮兩個(gè)問題：①如何估計(jì)計(jì)算結(jié)果的精度，即給出該步的估計(jì)誤差；②如何根據(jù)所估計(jì)的誤差去調(diào)整步長。下面就從這兩個(gè)方面介紹如何在保證計(jì)算精度的前提下，對(duì)隱式Taylor級(jí)數(shù)暫態(tài)穩(wěn)定算法進(jìn)行變步長控制。

2.1 誤差估計(jì)

步長要實(shí)現(xiàn)自動(dòng)調(diào)整，需要依據(jù)該步的估計(jì)誤差和給定誤差ε的比較來進(jìn)行，而其中的誤差估計(jì)是關(guān)鍵。在實(shí)際計(jì)算中，常采用外推法估計(jì)誤差，即用同一個(gè)公式，從節(jié)點(diǎn)xn出發(fā)，分別用步長h計(jì)算一步，用h／2計(jì)算兩步到同一節(jié)點(diǎn)xn＋1。然后將二者結(jié)果進(jìn)行比較近似的確定誤差，但是這個(gè)方法有花費(fèi)時(shí)間多的弊病。

本文針對(duì)隱式Taylor級(jí)數(shù)暫態(tài)穩(wěn)定算法的預(yù)估 －校正的迭代格式，提出一種新的誤差估計(jì)方法，此方法與外推法相比，可節(jié)省很大的計(jì)算量。

如圖1所示，其中（xk，ηk）為解曲線上的初值，F(xiàn)（xk，ηk）和D（xk，ηk）分別表示某數(shù)值解方向和精確解方向。以步長h計(jì)算一步產(chǎn)生的局部離散誤差如式（8）所示。

該數(shù)值方法的整體離散誤差如式（9）所示。

一般情況下，常微分方程的解析解很難得到。因此，精確解方向D（xk，ηk）也很難精確確定。在進(jìn)行步長動(dòng)態(tài)控制作用一般采用比較精確的數(shù)值解近似替代精確解方法。進(jìn)而可以采用公式（9）進(jìn)行誤差控制。

隱式調(diào)諧Taylor級(jí)數(shù)暫態(tài)穩(wěn)定算法的迭代格式及變量定義參照文獻(xiàn)［6］。此算法是一個(gè)具有高精度及A穩(wěn)定性的優(yōu)秀算法，當(dāng)預(yù)估式即顯式Taylor級(jí)數(shù)的最高導(dǎo)數(shù)階數(shù)為p階時(shí)，相應(yīng)的校正式即隱式Taylor級(jí)數(shù)可取得2p階精度。利用隱式Taylor級(jí)數(shù)法的預(yù)估校正迭代格式及A穩(wěn)定性和高精度特點(diǎn)，推導(dǎo)下述理論。

對(duì)于隱式Taylor級(jí)數(shù)暫態(tài)穩(wěn)定算法的預(yù)估式，設(shè)預(yù)估式的最高導(dǎo)數(shù)階數(shù)為p階，在xk＋1處得到的數(shù)值解為y1，則有下式成立：

其中，d1（xk）與函數(shù)的高階導(dǎo)數(shù)有關(guān)。

當(dāng)預(yù)估式的最高導(dǎo)數(shù)階數(shù)為p階時(shí)，校正式的精度為2p階，設(shè)校正式在xk＋1處得到的數(shù)值解為y2，則有下式成立：

式（10）減去式（11）得

將式（12）帶入式（10）得

圖1 局部離散誤差示意圖Fig.1 Schematic diagram of the local discretization error

由式（13）可見，利用預(yù)估式和校正式在xk＋1處得到的數(shù)值解的差值，可估計(jì)精確解與數(shù)值解的局部誤差。與傳統(tǒng)外推法相比，此方法充分利用了預(yù)估校正的迭代格式，不需每步都用h／2計(jì)算兩步到同一節(jié)點(diǎn)來估計(jì)誤差，從而節(jié)省了大量計(jì)算時(shí)間。

2.2 步長的動(dòng)態(tài)控制

在上一節(jié)誤差估計(jì)δ的基礎(chǔ)上，依據(jù)給定的精度要求ε1和ε2即可實(shí)現(xiàn)步長的自動(dòng)調(diào)整。

步長的調(diào)整方法通常采用加倍和折半法。但是當(dāng)初始步長已經(jīng)接近于合適值時(shí)，這種變步長方法將很不合理。比如，當(dāng)估計(jì)誤差δ比精度要求ε略小時(shí)，將步長加倍后極有可能出現(xiàn)δ?ε的情況，這時(shí)不得不再將步長折半一次，回到原來的步長；同理，當(dāng)估計(jì)誤差δ比精度要求ε略大時(shí)，也會(huì)出現(xiàn)類似情況而增加后繼的計(jì)算量。由于連貫性，相鄰各步的合適步長的變化通常是比較平緩的，而加倍和折半處理使得步長變化過快，在整個(gè)計(jì)算過程中，這種折半和加倍的變步長方法將導(dǎo)致大量回退性的計(jì)算，增加計(jì)算量。因此，采用較小的步長變化量進(jìn)行連續(xù)變化更加合理。當(dāng)需要增加和減少步長時(shí)，本文取0.01s作為每次步長的變化量，并且配合相應(yīng)的修正策略一起進(jìn)行步長控制。

步長的動(dòng)態(tài)控制判據(jù)如下：

1）如果δ＜ε2，則表示當(dāng)前步長時(shí)該步的估計(jì)誤差比精度要求小。將步長增大（可取h＝h＋0.01）進(jìn)入下一步計(jì)算。

2）如果ε2＜δ＜ε1，則表示當(dāng)前步長時(shí)該步的誤差滿足精度要求。步長保持不變，進(jìn)入下一步計(jì)算。

3）如果δ＞ε1，則表示當(dāng)前步長時(shí)該步的估計(jì)誤差比精度要求大。將步長減小（可取h＝h－0.01），進(jìn)入下一步計(jì)算。

2.3 變步長仿真中的修正策略

在程序變步長仿真過程中，為避免步長頻繁切換，以加快收斂速度，本文采取如下修正策略：

（1）在仿真到網(wǎng)絡(luò)故障時(shí)刻時(shí)，設(shè)置一個(gè)恒定步長h（可取0.01s），該步長維持到故障切除時(shí)刻，再進(jìn)行可變步長的判別修改（因?yàn)榇藭r(shí)網(wǎng)絡(luò)的變化可能會(huì)導(dǎo)致故障開始時(shí)期和故障恢復(fù)初期的狀態(tài)變量變化較大，從而影響到仿真速度和精度）。

（2）設(shè)置最小步長hmin（取0.005s）和最大步長hmax（取0.08s）。雖然算法具有A穩(wěn)定性，但是當(dāng)步長過大時(shí)可能會(huì)出現(xiàn)收斂困難的情況；而步長過小會(huì)使仿真步數(shù)太多，計(jì)算量增加。因此，在變步長仿真中設(shè)置最大和最小的步長限制。

3 算例分析

本文采用了忽略發(fā)電機(jī)次暫態(tài)過程的實(shí)心轉(zhuǎn)子電機(jī)模型，并采用IEEEI型勵(lì)磁系統(tǒng)。具體公式和參數(shù)詳見文獻(xiàn)［3］。

測(cè)試系統(tǒng)采用某省網(wǎng)72臺(tái)發(fā)電機(jī)、432節(jié)點(diǎn)系統(tǒng)。計(jì)算機(jī)采用Intel（R）Core（TM）2Duo 2.20GHz處理器，操作系統(tǒng)采用Windows XP sp3。開發(fā)語言采用C＋＋，開發(fā)工具采用Visual C＋＋6.0。

仿真時(shí)間20s，設(shè)定故障為482號(hào)支路50%處（距482號(hào)支路電氣距離最近的是14號(hào)發(fā)電機(jī)），5 s時(shí)發(fā)生三相短路故障，5.12s時(shí)故障消除。取仿真迭代收斂精度為10－5，各狀態(tài)變量的單步誤差精度要求ε1取10－5、ε2取10－6，最高導(dǎo)數(shù)階數(shù)為5階。

3.1 仿真精度對(duì)比

在定步長的情況下，程序采用不同的仿真步長時(shí)仿真精度是不一樣的。本節(jié)給出了定步長時(shí)步長為0.01s和0.02s時(shí)的仿真結(jié)果，最后給出了程序在變步長時(shí)的仿真結(jié)果。

步長為0.01s時(shí)，14號(hào)發(fā)電機(jī)的轉(zhuǎn)速ω仿真曲線如圖2所示。

圖2 定步長0.01s時(shí)，14號(hào)發(fā)電機(jī)的轉(zhuǎn)速ω仿真曲線Fig.2 Simulation curve for the rotation rateωof No.14 generator at the fixed－step of 0.01s

步長為0.02s時(shí)，14號(hào)發(fā)電機(jī)的轉(zhuǎn)速ω故障曲線如圖3所示。

圖3 定步長0.02s時(shí)，14號(hào)發(fā)電機(jī)的轉(zhuǎn)速ω仿真曲線Fig.3 Simulation curve for the rotation rateωof No.14 generator at the fixed－step of 0.02s

變步長時(shí)，14號(hào)發(fā)電機(jī)的轉(zhuǎn)速ω故障曲線如圖4所示。

圖4 變步長時(shí)，14號(hào)發(fā)電機(jī)的轉(zhuǎn)速ω仿真曲線Fig.4 Simulation curve for the rotation rateωof No.14 generator at variable steps

由圖2和圖3可見，在定步長0.01s的仿真曲線中，最大和最小峰值分別為1.0042和0.9955，而定步長0.02s的仿真曲線的最大和最小峰值卻只有1.003和0.997，并且曲線在隨后的震蕩過程中簡(jiǎn)化太多。這充分說明定步長0.02s的仿真精度下降。然而，由圖4可見，變步長的仿真曲線和定步長0.01s時(shí)的曲線幾乎完全一致，變步長可以取得類似定步長0.01s時(shí)的計(jì)算精度。在變步長仿真過程中，最小步長0.005s出現(xiàn)在故障后曲線震蕩較劇烈的時(shí)段；最大步長0.08s出現(xiàn)在故障前和故障后的曲線平穩(wěn)時(shí)期。

3.2 仿真速度對(duì)比

由上節(jié)可知，定步長0.02s較0.01s的仿真精度下降。而變步長可以取得定步長0.01s時(shí)較高的計(jì)算精度，并且，變步長在計(jì)算速度上也提高了很多。三種情況下的計(jì)算速度對(duì)比見表1。

表1 定步長和變步長的總仿真步數(shù)和迭代次數(shù)對(duì)比Tab.1 Contrast of the simulation step number and the iteration number between fixed－step and variable step

從表1的仿真結(jié)果可看出，定步長0.02s較定步長0.01s在仿真速度上提高了將近一倍，但是由上節(jié)的仿真曲線對(duì)比可知，其仿真精度下降。這充分說明：算法在定步長時(shí)為了取得較高的計(jì)算精度，只能采用小步長仿真，從而不能兼顧計(jì)算速度。仿真20s，定步長0.01s需要仿真2000步，總共迭代2031次；而變步長只需要仿真631步，總共迭代670次。與定步長0.01s相比，變步長的計(jì)算速度提高了67%，顯著地提高了計(jì)算效益。

4 結(jié)語

本文首先對(duì)算法的計(jì)算量進(jìn)行了分析，得出步長變化與計(jì)算時(shí)間的關(guān)系，然后結(jié)合算法的預(yù)估校正迭代格式提出了一種新的誤差估計(jì)方法，實(shí)現(xiàn)了隱式Taylor級(jí)數(shù)暫態(tài)穩(wěn)定算法的變步長控制，并給出了步長修正策略。結(jié)合仿真算例，從計(jì)算精度和計(jì)算速度兩個(gè)方面給出了定步長和變步長的對(duì)比結(jié)果。仿真結(jié)果表明：算法在定步長時(shí)，為了提高仿真速度而采用大步長仿真，最終導(dǎo)致仿真精度下降；而在變步長時(shí)，不僅保持了算法的高精度特點(diǎn)，同時(shí)其計(jì)算速度提高了67%，有效地提高了暫態(tài)穩(wěn)定計(jì)算的效率。

［1］ Daozhi Xia，Heydt G T.On－line transient stability evaluation by system decomposition－aggregation and high order derivatives［J］.IEEE Trans on Power Apparatus and Systems，1983，102（7）：2038－2046.

［2］ 白 雪 峰，郭 志 忠，王 永 剛 （Bai Xuefeng，Guo Zhizhong，Wang Yonggang）.一種簡(jiǎn)捷快速的暫態(tài)穩(wěn)定仿真算法（A new fast transient stability simulation method）［J］.繼電器（Relay），2000，28（1）：1－2，45.

［3］ 郭志忠，柳焯（Guo Zhizhong，Liu Zhuo）.快速高階Taylor級(jí)數(shù)法暫態(tài)穩(wěn)定計(jì)算（Fast transient stability simulation by higher order Taylor series expansions）［J］.中國電機(jī)工程學(xué)報(bào)（Proceedings of the CSEE），1991，11（3）：8－16.

［4］ 于繼來，郭志忠，柳焯（Yu Jilai，Guo Zhizhong，Liu Zhuo）.基于能量函數(shù)高階Taylor級(jí)數(shù)展開技術(shù)的直接法（Direct method with high order Taylor series expansions of energy function）［J］.電網(wǎng)技 術(shù)（Power System Technology），1995，19（2）：18－20，24.

［5］ 王宇賓，常鮮戎，羅艷，等（Wang Yubin，Chang Xianrong，Luo Yan，et al）.基于隱式Taylor級(jí)數(shù)法的電力系統(tǒng)暫態(tài)穩(wěn)定計(jì)算（An implicit Taylor series method for simulation of power system transient）［J］.華北電力大學(xué)學(xué)報(bào)（Journal of North China Electric Power University），2005，32（2）：1－6.

［6］ 鄭煥坤，常鮮戎，王正輝（Zheng Huankun，Chang Xianrong，Wang Zhenghui）.高精度A穩(wěn)定隱式調(diào)諧Taylor級(jí)數(shù)法在電力系統(tǒng)中的應(yīng)用（Application of high precision and A－stability implicit tuned Taylor series method in power system）［J］.電工技術(shù)學(xué)報(bào)（Transactions of China Electrotechnical Society），2012，27（1）：217－223，230.

［7］ 鄧陽，吳政球，容文光，等（Deng Yang，Wu Zhengqiu，Rong Wenguang，et al）.基于泰勒級(jí)數(shù)的N－1網(wǎng)絡(luò)快速靈敏度修正計(jì)算（Sensitivity of N－1system fast correction calculation based on Taylor series）［J］.繼電器（Relay），2007，35（16）：23－26.

［8］ 宋新立，湯涌，劉文卓，等（Song Xinli，Tang Yong，Liu Wenzhuo，et al）.電力系統(tǒng)全過程動(dòng)態(tài)仿真的組合數(shù)值積分算法研究（Mixed numerical integral algorithm for full dynamic simulation of the power system）［J］.中國電機(jī)工程學(xué)報(bào) （Proceedings of the CSEE），2009，29（28）：23－29.

［9］ 徐箭，陳允平（Xu Jian，Chen Yunping）.基于改進(jìn)通信算法的暫態(tài)穩(wěn)定并行仿真（Parallel simulation for transient stability based on improved communication algorithm in cluster environment）［J］.中國電機(jī)工程學(xué)報(bào)（Proceedings of the CSEE），2006，26（15）：12－18.

［10］ 白雪峰，郭志忠（Bai Xuefeng，Guo Zhizhong）.Taylor級(jí)數(shù)法暫態(tài)穩(wěn)定計(jì)算中階數(shù)的動(dòng)態(tài)控制（The dynamic control of order selection in fast transient stability simulation by higher order Taylor series expansions）［J］.電力系統(tǒng)自動(dòng)化（Automation of Electric Power Systems），1999，23（22）：5－7.

［11］周全仁，張清益.電網(wǎng)計(jì)算與程序設(shè)計(jì)［M］.長沙：湖南技術(shù)出版社，1984.

［12］黃舒予，牟龍華，石林（Huang Shuyu，Mu Longhua，Shi Lin）.自適應(yīng)變步長MPPT算法（Adaptive variable step size MPPT algorithm）［J］.電 力 系 統(tǒng) 及其 自 動(dòng) 化 學(xué) 報(bào)（Proceedings of the CSU－EPSA），2011，23（5）：26－30.

［13］ 王晨煒，靳希（Wang Chenwei，Jin Xi）.支持向量機(jī)動(dòng)態(tài)訓(xùn)練算法電力系統(tǒng)暫態(tài)穩(wěn)定評(píng)估（New algorithm in support vectors machine dynamic training for transient stability assessment）［J］.電 力 系 統(tǒng)及 其 自 動(dòng) 化 學(xué) 報(bào)（Proceedings of the CSU－EPSA），2009，21（2）：31－34.

［14］ 李傳棟，楊金剛，陳家榮（Li Chuandong，Yang Jingang，K.W.CHAN）.暫態(tài)穩(wěn)定約束下的最優(yōu)潮流（Optimal power flow with transient stability constraints）［J］.電 力 系 統(tǒng) 及 其 自 動(dòng) 化 學(xué) 報(bào)（Proceedings of the CSU－EPSA），2009，21（6）：45－50.

［15］汪芳宗，何仰贊，陳德樹（Wang Fangzong，He Yangzan，Chen Deshu）.空間時(shí)間混合變步長暫態(tài)穩(wěn)定仿真新算法（Power system transient stability simulation with space－and－time varying step method）［J］.電力系統(tǒng)及其自動(dòng)化學(xué)報(bào)（Proceedings of the CSU－EPSA），1991，3 （1）：24－30，23.