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

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