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

鄭煥坤，孫耀芹，常鮮戎

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

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

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

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

1 算法的計算量分析

隱式Taylor級數(shù)暫態(tài)穩(wěn)定算法的計算量主要集中在以下兩個方面：①發(fā)電機(jī)各狀態(tài)變量的各階導(dǎo)數(shù)的遞推求??；②形成因子表并利用它對各階導(dǎo)數(shù)進(jìn)行網(wǎng)絡(luò)方程的解算。接下來針對這兩個方面對算法的計算量進(jìn)行分析。

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

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

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

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

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

則對于每一步的計算量為

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

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

2 算法的變步長控制

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

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

2.1 誤差估計

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

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

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

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

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

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

對于隱式Taylor級數(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階時，校正式的精度為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ù)值解的差值，可估計精確解與數(shù)值解的局部誤差。與傳統(tǒng)外推法相比，此方法充分利用了預(yù)估校正的迭代格式，不需每步都用h／2計算兩步到同一節(jié)點來估計誤差，從而節(jié)省了大量計算時間。

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

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

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

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

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

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

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

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

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

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

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

3 算例分析

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

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

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

3.1 仿真精度對比

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

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

圖2 定步長0.01s時，14號發(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時，14號發(fā)電機(jī)的轉(zhuǎn)速ω故障曲線如圖3所示。

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

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

圖4 變步長時，14號發(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，并且曲線在隨后的震蕩過程中簡化太多。這充分說明定步長0.02s的仿真精度下降。然而，由圖4可見，變步長的仿真曲線和定步長0.01s時的曲線幾乎完全一致，變步長可以取得類似定步長0.01s時的計算精度。在變步長仿真過程中，最小步長0.005s出現(xiàn)在故障后曲線震蕩較劇烈的時段；最大步長0.08s出現(xiàn)在故障前和故障后的曲線平穩(wěn)時期。

3.2 仿真速度對比

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

表1 定步長和變步長的總仿真步數(shù)和迭代次數(shù)對比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é)的仿真曲線對比可知，其仿真精度下降。這充分說明：算法在定步長時為了取得較高的計算精度，只能采用小步長仿真，從而不能兼顧計算速度。仿真20s，定步長0.01s需要仿真2000步，總共迭代2031次；而變步長只需要仿真631步，總共迭代670次。與定步長0.01s相比，變步長的計算速度提高了67%，顯著地提高了計算效益。

4 結(jié)語

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

［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）.一種簡捷快速的暫態(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級數(shù)法暫態(tài)穩(wěn)定計算（Fast transient stability simulation by higher order Taylor series expansions）［J］.中國電機(jī)工程學(xué)報（Proceedings of the CSEE），1991，11（3）：8－16.

［4］ 于繼來，郭志忠，柳焯（Yu Jilai，Guo Zhizhong，Liu Zhuo）.基于能量函數(shù)高階Taylor級數(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級數(shù)法的電力系統(tǒng)暫態(tài)穩(wěn)定計算（An implicit Taylor series method for simulation of power system transient）［J］.華北電力大學(xué)學(xué)報（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級數(shù)法在電力系統(tǒng)中的應(yīng)用（Application of high precision and A－stability implicit tuned Taylor series method in power system）［J］.電工技術(shù)學(xué)報（Transactions of China Electrotechnical Society），2012，27（1）：217－223，230.

［7］ 鄧陽，吳政球，容文光，等（Deng Yang，Wu Zhengqiu，Rong Wenguang，et al）.基于泰勒級數(shù)的N－1網(wǎng)絡(luò)快速靈敏度修正計算（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)全過程動態(tài)仿真的組合數(shù)值積分算法研究（Mixed numerical integral algorithm for full dynamic simulation of the power system）［J］.中國電機(jī)工程學(xué)報 （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é)報（Proceedings of the CSEE），2006，26（15）：12－18.

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

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

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

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

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