第一講微分方程

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

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

3、是一種遞推算法，它使求解過程能順著結(jié)點(diǎn)的順序一步一步向前推進(jìn)，即可用前一個(gè)結(jié)點(diǎn)上的值(單步法)或者前面幾個(gè)結(jié)點(diǎn)上的值(多步法)來計(jì)算當(dāng)前步上的近似值。1-1-2 單步法所謂單步法就是用前一結(jié)點(diǎn)上的值來計(jì)算當(dāng)前結(jié)點(diǎn)上的近似值的方法。1歐拉法歐拉法是一種最簡單的單步法。對一階微分方程（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è)想，可以通過遞推計(jì)算方法得出微分方程的近似解，這種方法稱為

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

5、可以設(shè)想在，區(qū)間內(nèi)，取不同的后再求加權(quán)平均值，用它作為該區(qū)間的平均增量，即（1-11）其中這種計(jì)算式稱為自啟動的N階龍格一庫塔法。當(dāng)N=l時(shí)為歐拉法，當(dāng)N=2時(shí)為改進(jìn)歐拉法。常采用的各階龍格一庫塔法計(jì)算公式的系數(shù)如表l所示。表1 各階龍格一庫塔法的系數(shù)表通常采用的名稱10歐拉法2改進(jìn)歐拉法3三階龍格一庫塔法4四階龍格一庫塔法四階龍格一庫塔法（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)動態(tài)仿真計(jì)算中，常采用的多步計(jì)算方法是用四點(diǎn)亞當(dāng)姆斯顯式公式計(jì)算預(yù)測值，用四點(diǎn)亞當(dāng)姆斯隱式公式計(jì)算校正值。為了提高精度，在誤差分析的基礎(chǔ)上，分別用亞當(dāng)姆斯四點(diǎn)顯式公式計(jì)算預(yù)測值及其修正后的預(yù)測值，用亞當(dāng)姆斯四點(diǎn)隱式公式計(jì)算校正值及其修正后的終值。計(jì)算過程如下：1) 計(jì)算2) 用顯式公式計(jì)算預(yù)測值 3) 計(jì)算誤差修正后的預(yù)測值 ，其中是在第步計(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ù)來計(jì)算，精度較高。但是在突變點(diǎn)，因?yàn)樽兞客蛔儯懊鎺讉€(gè)步長的數(shù)據(jù)不能用來估算本步長的值。多步法一般需要用同樣精度的單步法啟動。按選取的插值結(jié)點(diǎn)是否包含結(jié)點(diǎn)來劃分，分為顯式法和隱式法。如果不包含結(jié)點(diǎn)，這種算法稱為顯式法。因?yàn)轱@式公式等號右邊不包含未知量，可以通過遞推方法連續(xù)計(jì)算下去。如果包含結(jié)點(diǎn)，所導(dǎo)出的數(shù)值積分公式稱為隱式公式，公式等號兩邊都含有未知量，不能直接通過遞推過程連續(xù)計(jì)算，一般先由顯式公式求得預(yù)測值，然后代入隱式公式求出校正值。第

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

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

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

11、方程（1-14）用顯式歐拉法求解時(shí)，計(jì)算公式為（1-15）計(jì)算時(shí)，實(shí)際上只能得到近似值（1-16）誤差為（1-17）只有當(dāng)時(shí)才能夠使誤差在計(jì)算過程中衰減，誤差不會增殖。也就是說，是數(shù)值穩(wěn)定的。因此采用顯式歐拉法時(shí)，步長不能取得太大。用隱式梯形法計(jì)算時(shí)，計(jì)算公式為隱式公式（1-18）（1-19）（1-20）顯然， 。由此可見，在計(jì)算過程中誤差是衰減的，能保證數(shù)值穩(wěn)定性。這也意味著隱式梯形法可以采取較大步長。2數(shù)值穩(wěn)定域考慮一階微分方程 是復(fù)數(shù)（1-21）用歐拉法求解時(shí)，為保證數(shù)值穩(wěn)定性，要求。此不等式在復(fù)平面上表示為以=-1為中心，以1為半徑的圓，該區(qū)域稱為絕對穩(wěn)定區(qū)域。絕對穩(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)考慮到迭代過程的收斂條件 （1-25）其中為李氏常數(shù)，當(dāng)時(shí)迭代過程收斂。為說明后退歐拉法的數(shù)值穩(wěn)定域，有 （1-26）誤差方程關(guān)系是 （1-27）只要，后退歐拉方法是A穩(wěn)定的。由于迭代時(shí)要求，在比較大時(shí)，仍要受到限制。隱式梯形法也是A穩(wěn)定的，但迭代過程收斂的條件是。1-2-3剛性方程問題1剛

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

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

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

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