系統(tǒng)仿真技術(shù)第6章 病態(tài)系統(tǒng)仿真.ppt

系統(tǒng)仿真技術(shù) 第6章 病態(tài)系統(tǒng)仿真,剡昌鋒 劉軍 蘭州理工大學(xué)機(jī)電工程學(xué)院,6.1病態(tài)系統(tǒng)的定義,系統(tǒng)中各環(huán)節(jié)的時(shí)間常數(shù)差異巨大。為保證仿真計(jì)算的穩(wěn)定性，由于仿真步長(zhǎng)必須限制在最小時(shí)間常數(shù)的數(shù)量級(jí)而選得很小，然而仿真結(jié)束的時(shí)間則決定于系統(tǒng)中的最大時(shí)間常數(shù)，若按滿足穩(wěn)定性要求所選擇的步長(zhǎng)進(jìn)行仿真，則不僅整個(gè)仿真所花費(fèi)的時(shí)間非常長(zhǎng)，甚至由于計(jì)算的舍入誤差而導(dǎo)致整個(gè)仿真的失敗。這就是所謂“病態(tài)系統(tǒng)仿真”問題。,病態(tài)系統(tǒng)的定義（續(xù)）,病態(tài)系統(tǒng)定義: （1） 令 （2） J稱為系統(tǒng)的雅可比矩陣。 若J的特征值全部具有負(fù)實(shí)部，且有： ，則該系統(tǒng)稱為病態(tài)系統(tǒng)，在某些文獻(xiàn)中也叫做剛性系統(tǒng)(Stiff)，而: 稱為病態(tài)比，一般在50以上。,6.2 線性病態(tài)系統(tǒng)仿真,對(duì)線性定常系統(tǒng)，我們可用如下狀態(tài)方程進(jìn)行一般性描述： (3) 1. 增廣矩陣法 將u(t)作為系統(tǒng)的增廣狀態(tài),線性病態(tài)系統(tǒng)仿真(續(xù)）,其中: 由于病態(tài)系統(tǒng)特征值相差倍數(shù)很大，必須用加速收斂的方法計(jì)算該狀態(tài)轉(zhuǎn)移矩陣的值。 2.蛙跳算法 基本思想： （1）考慮作用函數(shù)為階躍函數(shù)，則增廣狀態(tài)十分簡(jiǎn)單。 選擇q，使得：,線性病態(tài)系統(tǒng)仿真(續(xù)）,(3) 仿真計(jì)算時(shí)采用如下“蛙跳”方式：,線性病態(tài)系統(tǒng)仿真(續(xù)）,即在qh以前采用加倍跳躍式計(jì)算，而在qh以后每隔qh計(jì)算一次。 優(yōu)點(diǎn)是： h可以取得很小（可按最小時(shí)間常數(shù)考慮），從而保證初始階段的精度而計(jì)算量卻不大，而到qh以后，小時(shí)間常數(shù)的作用完成，則加大步長(zhǎng)計(jì)算，從而加快仿真計(jì)算速度。 從x(h)到x(qh)都是以x(0)為基礎(chǔ)進(jìn)行計(jì)算，所以誤差傳播比較均勻（僅僅是狀態(tài)轉(zhuǎn)移矩陣的誤差）。,6.3 非線性病態(tài)系統(tǒng)仿真,一般非線性系統(tǒng)的仿真大多采用數(shù)值積分法。而數(shù)值積分法一般又只具有有限的穩(wěn)定域，典型的如龍格庫塔法，仿真步長(zhǎng)限定在系統(tǒng)最小時(shí)間常數(shù)的數(shù)量級(jí)，才能保證計(jì)算的穩(wěn)定性，而系統(tǒng)的過渡過程時(shí)間卻決定于最大時(shí)間常數(shù)，因而對(duì)病態(tài)系統(tǒng)來說計(jì)算量極大，加上存在誤差傳播，仿真的精度甚至穩(wěn)定性也會(huì)受到影響。,6.3.1 吉爾(Gear)法,6.3.1.1.Stiff穩(wěn)定域 Gear研究后發(fā)現(xiàn)，并不要求一定采用恒穩(wěn)方法，而只要具有所謂Stiff穩(wěn)定域就可以了。 Stiff穩(wěn)定域定義：對(duì)實(shí)際的物理系統(tǒng)，時(shí)間常數(shù)一般小于零。選擇仿真步長(zhǎng)h若滿足： 與 可保證仿真的穩(wěn)定性，稱該算法具有stiff穩(wěn)定域。,Stiff穩(wěn)定域（續(xù)）,實(shí)際上具有Stiff穩(wěn)定域的方法與恒穩(wěn)方法只在近虛軸處有一點(diǎn)差別，即如果系統(tǒng)中的極點(diǎn)全部為實(shí)極點(diǎn)，那么無論選擇多大的步長(zhǎng)，計(jì)算是恒穩(wěn)的。如果系統(tǒng)中有復(fù)極點(diǎn)（實(shí)部仍為負(fù)數(shù)），只要步長(zhǎng)的選擇滿足上述條件，也能保證算法穩(wěn)定。,Stiff穩(wěn)定域（續(xù)）,Stiff域中與的確定： 按病態(tài)系統(tǒng)的大特征值來選擇步長(zhǎng)： 該特征值所對(duì)應(yīng)的模態(tài)大約要經(jīng)過4倍左右時(shí)間常數(shù)的時(shí)間才能有效地衰減掉，即 ，也就是 這樣，此時(shí)即使加大步長(zhǎng)h，也能保持計(jì)算的穩(wěn)定性?；谶@一考慮，可設(shè)- 4。,Stiff穩(wěn)定域（續(xù)）,另一方面，考慮到系統(tǒng)特征值為復(fù)數(shù)，它所對(duì)應(yīng)的瞬態(tài)響應(yīng)呈振蕩型。一個(gè)振蕩周期內(nèi)至少計(jì)算N個(gè)點(diǎn)。最小振蕩周期為： 其中h為計(jì)算步長(zhǎng)，若選擇N8，則有： ，因此可選/4。 綜上所述，如果選擇某一種方法，其穩(wěn)定域/4，且|4，則從使用的角度來看，圖6.1所示的穩(wěn)定域與恒穩(wěn)域沒有差別，從而完全可以用于病態(tài)系統(tǒng)的仿真。,6.3.1.2吉爾（Gear）法的基本原理,設(shè)系統(tǒng)： 滿足Stiff穩(wěn)定域的多步法：Gear提出的用于病態(tài)系統(tǒng)仿真的計(jì)算公式是： (1),用于病態(tài)系統(tǒng)仿真的Gear公式的系數(shù)表,吉爾（Gear）法的基本原理(續(xù)）,穩(wěn)定域如圖6.2所示。 從圖上可以看出，該方法在5階以下（包括5階）的穩(wěn)定域滿足Stiff穩(wěn)定域的條件（/4,|4，而且還可能穿過負(fù)實(shí)軸。,吉爾（Gear）法的基本原理(續(xù)）,在用Gear法仿真非線性病態(tài)系統(tǒng)時(shí)，有以下三個(gè)基本問題需要解決： 1）啟動(dòng)問題 上述Gear法本質(zhì)上是隱式多步法。 對(duì)于初值問題，困難：隱式方法一般用顯式方法啟動(dòng)，即先進(jìn)行預(yù)報(bào)，然后通過迭代進(jìn)行校正。如果迭代方法的收斂性不好，可能引起計(jì)算發(fā)散或計(jì)算量加大。 即使選擇的迭代方法收斂性滿足要求，顯式多步法預(yù)報(bào)，仍然難以啟動(dòng)，必須采用單步法啟動(dòng)，由于單步法不具有Stiff穩(wěn)定域，因而很難保證計(jì)算的穩(wěn)定性。 2）變步長(zhǎng)策略 非線性病態(tài)系統(tǒng)仿真往往采用變步長(zhǎng)策略，如何適時(shí)地將步長(zhǎng)調(diào)整到合適長(zhǎng)度，以同時(shí)滿足仿真精度和速度的要求。 3）加速迭代 為了提高計(jì)算效能，加速迭代也是非線性病態(tài)系統(tǒng)仿真中重要問題。,6.3.1.3 單步多值法,以三階為例，采用顯式多步法進(jìn)行預(yù)報(bào)，然后用隱式法校正。 其顯式預(yù)報(bào)的公式是： （2） 三階隱式Gear公式校正： （3）,單步多值法（續(xù)）,其中等式右邊第4項(xiàng)為導(dǎo)函數(shù)項(xiàng)，它是通過將第i次迭代所得到的y的預(yù)報(bào)值代入導(dǎo)函數(shù)后計(jì)算得到的。為便于程序?qū)崿F(xiàn)，由（2）式，并令：,單步多值法（續(xù)）,迭代的校正公式可表示成：,單步多值法（續(xù)）,對(duì)一般情形，令：,單步多值法（續(xù)）,其中p表示Gear法的階次。對(duì)三階預(yù)報(bào)迭代校正算法： 用矩陣的方法加以表示： （4） （5） 其中B，C分別是相應(yīng)的系數(shù)矩陣。,單步多值法（續(xù)）,然而，對(duì)初值問題，上式是不能自啟動(dòng)的。 解決方法-單步多值法：用高階導(dǎo)數(shù)值來取代前幾步的y及f的值。 先定義一個(gè)向量，稱之為Nordsieck向量： （6） （7） 需要確定Z向量與Y向量之間的關(guān)系。,單步多值法（續(xù)）,采用多項(xiàng)式逼近，以三階為例： （8） 在 處，有 （9）,單步多值法（續(xù)）,同樣，也可以得到： 由（8）式及（9）式消去 ，整理后可得到：,單步多值法（續(xù)）,寫成矩陣形式，就是： 簡(jiǎn)記為： （10） Q陣就是Y向量與Z向量之間的變換陣。對(duì)p階Gear法，Q陣為（p+1）階非奇異方陣。,單步多值法（續(xù)）,將（10）式代入（4）及（5）式，可得到單步多值法的計(jì)算公式： (11) 說明：（11）式中計(jì)算 時(shí)以 為自變量，但由于 是一個(gè)標(biāo)量，所以實(shí)際上只是用到它們的第一個(gè)分量 。,單步多值法（續(xù)）,若記 （12） 則（11）式可簡(jiǎn)寫為： （13） 在三階的情況下，P及L的值如下：,6.3.1.4 誤差估計(jì)與控制,當(dāng)用單值多步的Gear法對(duì)非線性病態(tài)系統(tǒng)進(jìn)行仿真時(shí)，它要求從初值開始，必須依靠顯式法來啟動(dòng)，即先從 及 開始，按一階公式計(jì)算，然后逐次升階，計(jì)算出 在這種升階過程中必須滿足誤差要求，對(duì)k階多步法，其截?cái)嗾`差為: (14) 其中 是t=點(diǎn)上y(t)的（p+1）階導(dǎo)數(shù)值，而是所討論區(qū)間中的某一個(gè)點(diǎn)。,誤差估計(jì)（續(xù)）,例如，對(duì)三階Gear法，就是 區(qū)間上的某一個(gè)點(diǎn)。若設(shè) ，則第k步的截?cái)嗾`差為： (15) 為了估計(jì) ，首先要估計(jì) ，已知： 現(xiàn)在要用 的差分來近似估計(jì) ，即：,誤差估計(jì)（續(xù)）,兩邊同乘以 則可得： （16） 將（16）式代入(15)式可得： （17） 仿真中一般采用的是相對(duì)誤差，若要求每一步的相對(duì)誤差不大于 ，即 （18） 其中 為到目前為止已出現(xiàn)過的y的最大絕對(duì)值（注：若y的初值為0，則應(yīng)取 1為宜。,誤差估計(jì)（續(xù)）,若系統(tǒng)為微分方程組，狀態(tài)變量y為N個(gè)（y(1)，y(2)，y(N)），相對(duì)誤差可定義為： （19） 若 ，則本步計(jì)算結(jié)果有效，進(jìn)入下一步；如果不滿足，則需要減小步長(zhǎng)或者采取其它措施。,誤差控制技術(shù),誤差控制：一是改變步長(zhǎng)，其二是改變階次。 無論是 或 ，均需要考慮變階或（與）變步長(zhǎng)，即通過改變p及h以使 。 變階或（與）變步長(zhǎng)的原則： 首先考慮僅改變步長(zhǎng)h（階次不變），設(shè)新的步長(zhǎng)為 ，且 。那么， 應(yīng)取多大為宜呢？,誤差控制技術(shù)（續(xù)）,為簡(jiǎn)便起見，我們根據(jù)單變量表達(dá)式（15）來分析，即: 即 根據(jù)（16）式， 故有：,誤差控制技術(shù)（續(xù)）,可得到： 考慮到誤差僅僅是估計(jì)值，考慮經(jīng)驗(yàn)系數(shù)1/1.2，這樣， 表達(dá)式如（20）式： （20） 若改變步長(zhǎng)的效果不理想，則考慮要提升仿真方法的階次。,誤差控制技術(shù)（續(xù)）,設(shè)當(dāng)前為p階，考慮用p1階計(jì)算。假設(shè)此時(shí)步長(zhǎng) 能使相對(duì)誤差接近規(guī)定的要求，即： 采用二階差分來近似 ，即：,誤差控制技術(shù)（續(xù)）,考慮到： 則不難得到： 也考慮經(jīng)驗(yàn)系數(shù)（這里取其為1/1.4，），則可得到升階時(shí) 的表達(dá)式如（21）式： （21） 如果仿真步長(zhǎng)縮短到相當(dāng)小而其誤差仍然達(dá)不到要求，則可能是因階次過高而穩(wěn)定域達(dá)不到要求的緣故。,誤差控制技術(shù)（續(xù)）,若將當(dāng)前仿真的階次p降低一階，仿真步長(zhǎng)h變?yōu)?，此時(shí)要使仿真誤差接近規(guī)定要求，則： 考慮經(jīng)驗(yàn)系數(shù)為1/1.3，則可得到： (22) 式中 ，它是向量 的最后一個(gè)分量。,誤差控制技術(shù)（續(xù)）,步長(zhǎng)發(fā)生變化時(shí)，如何由當(dāng)前的 產(chǎn)生在 下的 已知 只變步長(zhǎng)：,誤差控制技術(shù)（續(xù)）,降階變步長(zhǎng)：從p階降為p-1階，步長(zhǎng)由h變?yōu)?，這時(shí)， 、 用 表示：,誤差控制技術(shù)（續(xù)）,升階變步長(zhǎng)：當(dāng)由p階升到p+1階，步長(zhǎng)由h變?yōu)?，則這時(shí)， 、 用 表示： 然而， 是未知的，為此，必須由已知的數(shù)據(jù)來估計(jì)。若采用差分法，即：,誤差控制技術(shù)（續(xù)）,則得到 的表達(dá)式如下：,6.3.1.5 加速收斂問題,Gear法用于病態(tài)系統(tǒng)進(jìn)行仿真時(shí)，必須先用顯式公式啟動(dòng)，然后用隱式法進(jìn)行校正，經(jīng)過多次迭代，達(dá)到適當(dāng)精度后，再往前推進(jìn)。顯然，迭代的收斂速度極大地影響著仿真效能。 三階Gear法為例 (1) 其中 就是 的簡(jiǎn)寫，它一般是 的非線性函數(shù)，（1）式是一非線性方程。,加速收斂問題(續(xù)）,迭代計(jì)算時(shí)，我們采用的迭代公式為： （2） 其中 是向量L的第一個(gè)分量，這稱為Picard迭代。對(duì)該迭代過程： （3）,加速收斂問題(續(xù)）,為保證收斂，則要求： 即： 一般 大約為13之間，可見，h不能太大，否則迭代將不收斂或收斂速度十分緩慢。這大大限制了Gear法的有效性。 牛頓迭代法： （1）式可改寫為： （4）,加速收斂問題(續(xù)）,簡(jiǎn)記為： （5） 令 （6） 顯然，（5）式與g(y)=0同解。由牛頓迭代法，不難得到： （7）,加速收斂問題(續(xù)）,從(6)式求得 的表達(dá)式，并將 代入 及（7）式，可得： (8) 將 代入(5)式，所得結(jié)果在代入(8)式，可得：,加速收斂問題(續(xù)）,(9),加速收斂問題(續(xù)）,多了一個(gè)“加速”因子： 對(duì)微分方程組來說，y是N維向量，(9)式成為如下形式： (10) 其中Y，F(xiàn)，G為N維向量，I為N階單位陣，而 為NN的方陣，亦稱為雅可比矩陣：,加速收斂問題(續(xù)）,(11) 采用牛頓法： (12) 其中Z為（K1）1的向量，而G為標(biāo)量。,加速收斂問題(續(xù)）,對(duì)于y為N維的情況，則有： (13) 式中的Z為（K1）N的矩陣，而G為N1的向量。 牛頓法能加快迭代過程的收斂，這是該方法的優(yōu)點(diǎn)，然而，因每次迭代時(shí)都必須計(jì)算雅可比陣，還要計(jì)算矩陣求逆，計(jì)算量將增加較多。,6.3.1.6 病態(tài)性探測(cè),系統(tǒng)并不是每一步均處于病態(tài)。從數(shù)值計(jì)算的角度進(jìn)行定義“病態(tài)性”和“非病態(tài)性”： 如果用非Stiff法仿真，其步長(zhǎng)限制是由于計(jì)算精度原因引起的而不是由于穩(wěn)定性引起的，則稱系統(tǒng)呈現(xiàn)“非病態(tài)性”，反之，若步長(zhǎng)是因穩(wěn)定性而受到限制時(shí)，則系統(tǒng)呈現(xiàn)“病態(tài)性”。 當(dāng)系統(tǒng)呈現(xiàn)病態(tài)性時(shí)應(yīng)采用Stiff法，而當(dāng)系統(tǒng)呈現(xiàn)非病態(tài)性時(shí)應(yīng)采用非Stiff法。,病態(tài)性探測(cè) (續(xù)）,病態(tài)性探測(cè)的三種方法： 1）穩(wěn)定半徑法 非Stiff法的穩(wěn)定域的形狀接近半圓，其穩(wěn)定域與實(shí)軸的交點(diǎn)稱為穩(wěn)定半徑（ ）。若記為系統(tǒng)的最大特征值，為保證計(jì)算穩(wěn)定，則要求： (1) 各種方法的穩(wěn)定半徑已知，而的值可由雅可比陣的范數(shù) 來估計(jì)，即： 因此，在用stiff方法對(duì)病態(tài)系統(tǒng)仿真過程中，雅可比矩陣的值可以得到，則在下一步計(jì)算時(shí)，先用（1）式進(jìn)行判斷病態(tài)性，再?zèng)Q定是否仍采用stiff方法仿真。,病態(tài)性探測(cè) (續(xù)）,2）嵌入低階大穩(wěn)定域法 如果當(dāng)前采用的是非Stiff法，而計(jì)算誤差不能滿足要求，為判斷其穩(wěn)定性，可降低階次，如果此時(shí)誤差減少，則說