第一講微分方程

1、第一節(jié) 微分方程的數(shù)值計(jì)算1-1-1 引言電力系統(tǒng)的動(dòng)態(tài)過(guò)程通常用微分方程來(lái)描述。在微分方程中，自變量是時(shí)間，因變量是系統(tǒng)中各物理變量。對(duì)于階常微分方程，可以采用引入新變量的方法把它轉(zhuǎn)化為個(gè)一階常微分方程來(lái)求解。例如，對(duì)于二階常微分方程（1-1）可以表示成：（1-2）一般來(lái)說(shuō)，依據(jù)電工學(xué)和動(dòng)力學(xué)建立起來(lái)的電力系統(tǒng)數(shù)學(xué)模型是可以用來(lái)數(shù)字仿真的。但是，在復(fù)雜電力系統(tǒng)中，由于數(shù)學(xué)模型很復(fù)雜，仿真的計(jì)算量很大，影響了仿真效率。有些環(huán)節(jié)對(duì)于所研究的問(wèn)題影響很小，在建立數(shù)學(xué)模型時(shí)可以把它們忽略。這種簡(jiǎn)化嚴(yán)格來(lái)說(shuō)應(yīng)該用實(shí)際實(shí)驗(yàn)，或者用精確的物理模擬實(shí)驗(yàn)來(lái)檢驗(yàn)，以確定簡(jiǎn)化的合理性。數(shù)學(xué)模型中的參數(shù)是數(shù)字仿真的

2、基礎(chǔ)，一般可以通過(guò)各種測(cè)量和實(shí)驗(yàn)得到，但這是一件十分煩瑣的工作。目前，有些參數(shù)只能憑經(jīng)驗(yàn)給出，其合理性有待實(shí)驗(yàn)的嚴(yán)格檢驗(yàn)。電力系統(tǒng)仿真的另一個(gè)重要問(wèn)題是數(shù)學(xué)模型的求解問(wèn)題。對(duì)于微分方程，除少數(shù)可以得到解析解以外，大多數(shù)只能采用數(shù)值解法。早在18世紀(jì)末，就有很多人提出過(guò)求解微分方程精確而有效的數(shù)值積分方法。微分方程的積分是一簇曲線。通常，在初值確定之后，數(shù)值解才能夠確定。得到微分方程的數(shù)值近似解有兩種基本的方法。一種方法是把近似解表示成有限個(gè)獨(dú)立函數(shù)之和，另一種方法是差分法。差分方法是尋求在一系列離散點(diǎn)上的近似值，這些離散點(diǎn)稱為結(jié)點(diǎn)。在多數(shù)情況下，這些結(jié)點(diǎn)是等距的，即（1-3）稱為步長(zhǎng)。差分方法

3、是一種遞推算法，它使求解過(guò)程能順著結(jié)點(diǎn)的順序一步一步向前推進(jìn)，即可用前一個(gè)結(jié)點(diǎn)上的值(單步法)或者前面幾個(gè)結(jié)點(diǎn)上的值(多步法)來(lái)計(jì)算當(dāng)前步上的近似值。1-1-2 單步法所謂單步法就是用前一結(jié)點(diǎn)上的值來(lái)計(jì)算當(dāng)前結(jié)點(diǎn)上的近似值的方法。1歐拉法歐拉法是一種最簡(jiǎn)單的單步法。對(duì)一階微分方程（1-3）假定已給定，可計(jì)算出。把一階微分方程寫成數(shù)值積分法的計(jì)算形式（1-4）如果十分靠近，則有（1-5）當(dāng)然，也可以把一階微分方程寫成（1-6）即用時(shí)的導(dǎo)數(shù)表示在區(qū)間，內(nèi)的平均增量。結(jié)點(diǎn)和結(jié)點(diǎn)之間的關(guān)系可表示為（1-7）盡管這種方法精度不高，但卻提出了一種設(shè)想，可以通過(guò)遞推計(jì)算方法得出微分方程的近似解，這種方法稱為

4、歐拉法。因?yàn)槭怯谜劬€來(lái)近似表示曲線，故又可稱為折線法。2隱式梯形法從歐拉法可以看出，為提高精度，可以用和的平均值作為在區(qū)間，內(nèi)的平均增量值。（1-8）（1-9）在等式的兩邊均含有未知量，故稱為隱式梯形法。隱式梯形公式不能用遞推的方式直接做數(shù)值計(jì)算。如果用歐拉法求得的值作為預(yù)測(cè)值，則有：（1-10）這種方法稱為改進(jìn)歐拉法。利用歐拉法和改進(jìn)歐拉法求微分方程的數(shù)值解，如果選取步長(zhǎng)相同，改進(jìn)歐拉法的計(jì)算量大，但改進(jìn)歐拉法的精度高。反過(guò)來(lái)說(shuō)，如果兩種方法要求精度相同，則改進(jìn)歐拉法可以選取較大步長(zhǎng)，總的計(jì)算量可以節(jié)省，舍入誤差也可較小。3龍格一庫(kù)塔法隱式梯形法是把和的平均值作為在區(qū)間，內(nèi)的平均增量值。同理

5、可以設(shè)想在，區(qū)間內(nèi)，取不同的后再求加權(quán)平均值，用它作為該區(qū)間的平均增量，即（1-11）其中這種計(jì)算式稱為自啟動(dòng)的N階龍格一庫(kù)塔法。當(dāng)N=l時(shí)為歐拉法，當(dāng)N=2時(shí)為改進(jìn)歐拉法。常采用的各階龍格一庫(kù)塔法計(jì)算公式的系數(shù)如表l所示。表1 各階龍格一庫(kù)塔法的系數(shù)表通常采用的名稱10歐拉法2改進(jìn)歐拉法3三階龍格一庫(kù)塔法4四階龍格一庫(kù)塔法四階龍格一庫(kù)塔法（1-12）其中1-1-3多步法多步法的一般計(jì)算公式為（1-13）其中和如果為顯式公式，否則是隱式的。這種方法不僅用到前面一個(gè)結(jié)點(diǎn)上的值，而且還用到前面幾個(gè)結(jié)點(diǎn)上的值。常用的方法是，其它的，系數(shù)如表2所示 通常采用的名稱顯式1 歐拉法2梯形法3三點(diǎn)亞當(dāng)姆斯4

6、四點(diǎn)亞當(dāng)姆斯隱式1歐拉法2梯形法3三點(diǎn)亞當(dāng)姆斯4四點(diǎn)亞當(dāng)姆斯在電力系統(tǒng)動(dòng)態(tài)仿真計(jì)算中，常采用的多步計(jì)算方法是用四點(diǎn)亞當(dāng)姆斯顯式公式計(jì)算預(yù)測(cè)值，用四點(diǎn)亞當(dāng)姆斯隱式公式計(jì)算校正值。為了提高精度，在誤差分析的基礎(chǔ)上，分別用亞當(dāng)姆斯四點(diǎn)顯式公式計(jì)算預(yù)測(cè)值及其修正后的預(yù)測(cè)值，用亞當(dāng)姆斯四點(diǎn)隱式公式計(jì)算校正值及其修正后的終值。計(jì)算過(guò)程如下：1) 計(jì)算2) 用顯式公式計(jì)算預(yù)測(cè)值 3) 計(jì)算誤差修正后的預(yù)測(cè)值 ，其中是在第步計(jì)算中得到的校正值。4)計(jì)算5) 用隱式公式計(jì)算校正值 6) 計(jì)算誤差修正后的終值 1-1-4數(shù)值積分法的分類按照數(shù)值積分時(shí)被積函數(shù)的近似表達(dá)式的不同，分為單步法和多步法。若采用線性插值函

7、數(shù)代替被積函數(shù)，則可得出單步法計(jì)算公式：若采用高階插值函數(shù)近似代替被積函敷，可推導(dǎo)出多步法計(jì)算公式。從理論上講，多步法比單步法更為有效，可取用更多結(jié)點(diǎn)的數(shù)據(jù)來(lái)計(jì)算，精度較高。但是在突變點(diǎn)，因?yàn)樽兞客蛔?，前面幾個(gè)步長(zhǎng)的數(shù)據(jù)不能用來(lái)估算本步長(zhǎng)的值。多步法一般需要用同樣精度的單步法啟動(dòng)。按選取的插值結(jié)點(diǎn)是否包含結(jié)點(diǎn)來(lái)劃分，分為顯式法和隱式法。如果不包含結(jié)點(diǎn)，這種算法稱為顯式法。因?yàn)轱@式公式等號(hào)右邊不包含未知量，可以通過(guò)遞推方法連續(xù)計(jì)算下去。如果包含結(jié)點(diǎn)，所導(dǎo)出的數(shù)值積分公式稱為隱式公式，公式等號(hào)兩邊都含有未知量，不能直接通過(guò)遞推過(guò)程連續(xù)計(jì)算，一般先由顯式公式求得預(yù)測(cè)值，然后代入隱式公式求出校正值。第

8、二節(jié) 數(shù)值積分法的誤差、數(shù)值穩(wěn)定性與剛性在使用數(shù)值積分法中應(yīng)該注意分析各種不同算法的誤差。首先要了解各種計(jì)算方法中誤差的來(lái)源，其次要了解在相繼的各步長(zhǎng)計(jì)算中誤差的傳遞和影響，即數(shù)值穩(wěn)定性。這個(gè)問(wèn)題不僅與算法選取有關(guān)，而且與微分方程及步長(zhǎng)的選取有密切關(guān)系。第三是要了解動(dòng)態(tài)系統(tǒng)各環(huán)節(jié)時(shí)間常數(shù)的較大差異與取同一步長(zhǎng)計(jì)算的協(xié)調(diào)問(wèn)題，如果不協(xié)調(diào)將出現(xiàn)剛性問(wèn)題，得不到真實(shí)的解。這些問(wèn)題相互之間有一定聯(lián)系和影響，在電力系統(tǒng)動(dòng)態(tài)仿真研究中應(yīng)給予足夠重視。1-2-1解的誤差在電力系統(tǒng)仿真計(jì)算中誤差來(lái)自許多方面，主要有以下幾方面。1)數(shù)值計(jì)算中，例如數(shù)值微分，數(shù)值積分，無(wú)窮級(jí)數(shù)計(jì)算等都是通過(guò)有限次的近似計(jì)算得到近

9、似解。近似解與真解的誤差稱為截?cái)嗾`差。2)由于計(jì)算機(jī)表示數(shù)的位數(shù)有限，因此在運(yùn)算中出現(xiàn)了舍入誤差，這種誤差有時(shí)候會(huì)隨著運(yùn)算量的增加而增大。例如有些病態(tài)線性方程組在求解時(shí)，舍入誤差可以在計(jì)算過(guò)程中繁殖起來(lái)產(chǎn)生錯(cuò)誤的解。3)當(dāng)采用分割法求解微分方程組和代數(shù)方程組時(shí)，可能產(chǎn)生交接誤差。4)由于對(duì)仿真系統(tǒng)變量特性作了某些簡(jiǎn)化而造成的誤差，或者在給定步長(zhǎng)內(nèi)方程組線性化造成的誤差。這種誤差稱為近似化誤差。5)在動(dòng)態(tài)仿真中，有些變量有限值約束，如果限值約束不是恰好在求解時(shí)間間隔的結(jié)點(diǎn)上，這時(shí)不能算出也不能使用變量的限值和對(duì)限值的補(bǔ)償，由此引起的誤差稱為限值誤差。當(dāng)采用適當(dāng)小的步長(zhǎng)時(shí)，截?cái)嗾`差，交接誤差、近似

10、化誤差和限值誤差都可減小到容許范圍內(nèi)。但是，由于計(jì)算量增加會(huì)引起舍入誤差，對(duì)仿真的整體效果不一定有利，因此對(duì)不同的系統(tǒng)和不同仿真問(wèn)題，應(yīng)選取合理的步長(zhǎng)。1-2-2數(shù)值積分法的穩(wěn)定性1數(shù)值穩(wěn)定性的定義在數(shù)值積分法計(jì)算過(guò)程中，一個(gè)步長(zhǎng)終點(diǎn)上的誤差是上一步長(zhǎng)計(jì)算所引起的誤差和本步長(zhǎng)計(jì)算所留下誤差之和。這些相繼步長(zhǎng)上的誤差可能相互抵消而衰減下去，也可能不斷增殖而使結(jié)果無(wú)法使用。數(shù)值穩(wěn)定性是對(duì)方法而言的，如果用某種方法求解方程式的數(shù)值解，初值的微小變化對(duì)近似解只產(chǎn)生有界的變化，那么我們稱這種情況為穩(wěn)定的。數(shù)值穩(wěn)定性對(duì)于高階微分方程式是程復(fù)雜的，為了解釋這個(gè)問(wèn)題，用一個(gè)簡(jiǎn)單的一階方程式來(lái)說(shuō)明。考慮一階微分

11、方程（1-14）用顯式歐拉法求解時(shí)，計(jì)算公式為（1-15）計(jì)算時(shí)，實(shí)際上只能得到近似值（1-16）誤差為（1-17）只有當(dāng)時(shí)才能夠使誤差在計(jì)算過(guò)程中衰減，誤差不會(huì)增殖。也就是說(shuō)，是數(shù)值穩(wěn)定的。因此采用顯式歐拉法時(shí)，步長(zhǎng)不能取得太大。用隱式梯形法計(jì)算時(shí)，計(jì)算公式為隱式公式（1-18）（1-19）（1-20）顯然， 。由此可見(jiàn)，在計(jì)算過(guò)程中誤差是衰減的，能保證數(shù)值穩(wěn)定性。這也意味著隱式梯形法可以采取較大步長(zhǎng)。2數(shù)值穩(wěn)定域考慮一階微分方程 是復(fù)數(shù)（1-21）用歐拉法求解時(shí)，為保證數(shù)值穩(wěn)定性，要求。此不等式在復(fù)平面上表示為以=-1為中心，以1為半徑的圓，該區(qū)域稱為絕對(duì)穩(wěn)定區(qū)域。絕對(duì)穩(wěn)定區(qū)域越大，這種數(shù)

12、值方法的適應(yīng)性越強(qiáng)。如果穩(wěn)定區(qū)域是的整個(gè)左半復(fù)平面，我們稱這種數(shù)值方法是A穩(wěn)定的。因此，歐拉法不是A穩(wěn)定的，僅僅是以=-1為中心，以1為半徑的圓。后退歐拉法的計(jì)算公式為 （1-22）后退歐拉法是隱式方法，計(jì)算時(shí)先用顯式歐拉法計(jì)算初值，再用隱式公式作迭代，計(jì)算公式為 （1-23） （1-24）用迭代法求解時(shí)應(yīng)考慮到迭代過(guò)程的收斂條件 （1-25）其中為李氏常數(shù)，當(dāng)時(shí)迭代過(guò)程收斂。為說(shuō)明后退歐拉法的數(shù)值穩(wěn)定域，有 （1-26）誤差方程關(guān)系是 （1-27）只要，后退歐拉方法是A穩(wěn)定的。由于迭代時(shí)要求，在比較大時(shí)，仍要受到限制。隱式梯形法也是A穩(wěn)定的，但迭代過(guò)程收斂的條件是。1-2-3剛性方程問(wèn)題1剛

13、性方程問(wèn)題的定義在許多仿真計(jì)算中，特別是網(wǎng)絡(luò)分析和仿真中，常出現(xiàn)剛性方程式。微分方程式的剛性問(wèn)題在性質(zhì)上與代數(shù)方程式的病態(tài)問(wèn)題相類似，它同動(dòng)態(tài)系統(tǒng)中各環(huán)節(jié)時(shí)間常數(shù)的差別有很大關(guān)系。更確切說(shuō)，剛性問(wèn)題是用線性化方程的最大特征值和最小特征值的比率來(lái)衡量的。對(duì)于某個(gè)系統(tǒng)，如果有一些分量反應(yīng)很慢，另一些分量反應(yīng)很快，為了防止計(jì)算中截?cái)嗾`差增大必須根據(jù)反應(yīng)很快的分量選取步長(zhǎng)。例如，考慮如下系統(tǒng) （1-28）它的解是 從穩(wěn)定性考慮，對(duì)于，用歐拉法時(shí)，要滿足，即要求。而計(jì)算幾個(gè)步長(zhǎng)以后，的值已很快衰減到與相比可忽略不計(jì)的程度，但是由此點(diǎn)往后，由于穩(wěn)定性的要求，仍只好采用非常小的步長(zhǎng)，這就是剛性方程問(wèn)題。一般

14、來(lái)說(shuō)，對(duì)于一般的一階常微分方程組 （1-29）如果的特征值滿足 則稱該方程組為剛性方程組。把稱為剛性比。2解決剛性問(wèn)題的辦法解決剛性問(wèn)題的有效辦法是選取小步長(zhǎng)，使計(jì)算結(jié)果準(zhǔn)確反映系統(tǒng)中那些迅速變化的分量，使截?cái)嗾`差和其它誤差保持在安全的低水平上。另外，數(shù)值穩(wěn)定性好的數(shù)值積分法因?yàn)樵谟?jì)算中誤差不易增殖，容許每步計(jì)算有較大誤差。即，在總誤差相同的條件下，穩(wěn)定性好的數(shù)值積分法可以采用較大步長(zhǎng)。由此看來(lái)，解決剛性問(wèn)題的有效方法是選用對(duì)步長(zhǎng)不加限制的數(shù)值積分法。顯然，前面所討論的A穩(wěn)定的方法能解決剛性問(wèn)題，例如后退歐拉法和隱式梯形法。C.w.Gear提出的stiff穩(wěn)定方法也是解決剛性問(wèn)題的一種重要方法

15、。在電力系統(tǒng)數(shù)字仿真中，如果剛性問(wèn)題不嚴(yán)重，那么穩(wěn)定性好的計(jì)算方法優(yōu)越性也不明顯。電力系統(tǒng)暫態(tài)穩(wěn)定研究中的經(jīng)典模型(只考慮發(fā)電機(jī)轉(zhuǎn)子運(yùn)動(dòng)的模型)不存在剛性問(wèn)題，只有采用詳細(xì)模型才可能出現(xiàn)剛性問(wèn)題。有一些剛性問(wèn)題的影響常在大擾動(dòng)之后的很短時(shí)間內(nèi)才表現(xiàn)出來(lái)，隨后很快衰減趨向平穩(wěn)。解決這種剛性問(wèn)題采用穩(wěn)定性好的方法固然很好，但是有時(shí)除非步長(zhǎng)很小很小，否則這類數(shù)值積分也不能夠精確地計(jì)算出這些響應(yīng)的過(guò)程。因此，應(yīng)該根據(jù)有關(guān)變量對(duì)整個(gè)電力系統(tǒng)動(dòng)態(tài)過(guò)程的影響來(lái)采取適當(dāng)?shù)姆椒āＳ袝r(shí)候不能僅限于在數(shù)值算法上想辦法，實(shí)際上還可采用其它特殊方法，例如以下幾種辦法有時(shí)很有效。1)改造數(shù)學(xué)模型。根據(jù)經(jīng)驗(yàn)對(duì)系統(tǒng)各環(huán)節(jié)時(shí)間

16、常數(shù)進(jìn)行檢驗(yàn)，可以省略掉某些時(shí)間常數(shù)值很小，衰減很快的變量，使模型簡(jiǎn)化，或者人為適當(dāng)增大時(shí)間常數(shù)。更重要的是在仿真時(shí)對(duì)工程中所考慮的模型和求解方法做巧妙的配合。2)混合步長(zhǎng)。對(duì)于不同的系統(tǒng)采用不同的步長(zhǎng)。其做法是先給出標(biāo)準(zhǔn)步長(zhǎng)，在計(jì)算中按照系統(tǒng)各環(huán)節(jié)中最小的時(shí)間常數(shù)選擇最小步長(zhǎng)。在一個(gè)的間隔內(nèi)求解一次代數(shù)方程，其它時(shí)間采用外推法求出代數(shù)方程。當(dāng)然，這種方法只適用于分割法。3)混合積分算法。對(duì)于剛性方程適宜采用穩(wěn)定性好的數(shù)值積分方法，而對(duì)于非剛性方程適宜采用精度高的數(shù)值積分方法?；旌戏e分法就是對(duì)方程組中的剛性方程和非剛性方程采用不同算法。1-2-4精度、速度與方法選擇在電力系統(tǒng)動(dòng)態(tài)仿真中，評(píng)價(jià)一個(gè)程序性能的好壞除了其適用性，方便性以外，計(jì)算時(shí)間少，解的可靠性與精確性也是很重要的指標(biāo)。但是，這些要求很難同時(shí)滿足，而且它們之間往往是相互矛盾的。因?yàn)榉抡嬗?jì)算中每一步計(jì)算誤差是前面各步誤差的函數(shù)，所以要想使總的誤差限制在可以接受的范圍內(nèi)，唯一實(shí)用的辦法是對(duì)不同計(jì)算題目和程序采取實(shí)驗(yàn)的辦法來(lái)了解整個(gè)計(jì)算過(guò)程對(duì)于局部誤差會(huì)有多大的影響。多數(shù)短期動(dòng)態(tài)穩(wěn)定程序采用定步長(zhǎng)，在中期和長(zhǎng)期動(dòng)態(tài)仿真計(jì)算中，可以采用變步長(zhǎng)或自動(dòng)變步長(zhǎng)的方法，以保證計(jì)算的精度和機(jī)時(shí)消耗達(dá)到理想的程度。另外，代數(shù)方程的矩陣分解法、結(jié)點(diǎn)編號(hào)順序優(yōu)化和壓縮存貯等技術(shù)，對(duì)于不同的仿真計(jì)算題目也會(huì)有不