火電廠仿真(課件)

1、火電廠仿真太原理工大學電氣與動力工程學院熱能系太原理工大學電氣與動力工程學院熱能系3.提高機組的安全性、效率和使用壽命 由于培訓和研究開發(fā)試驗不在現(xiàn)場進行，可以避免現(xiàn)場的人為失誤、減少停機次數(shù)、提高設(shè)備的安全性、提高機組使用效率及壽命。 按真實系統(tǒng)的物理性質(zhì)構(gòu)造系統(tǒng)的物理模型，按真實系統(tǒng)的物理性質(zhì)構(gòu)造系統(tǒng)的物理模型，再現(xiàn)系統(tǒng)的一些特性。再現(xiàn)系統(tǒng)的一些特性。例如例如:新型鍋爐生產(chǎn)前，各部件按一定比例縮小，新型鍋爐生產(chǎn)前，各部件按一定比例縮小，構(gòu)成一個幾何相似的物理模型。構(gòu)成一個幾何相似的物理模型。 對此分析研究可發(fā)現(xiàn)該鍋爐在結(jié)構(gòu)和形態(tài)上對此分析研究可發(fā)現(xiàn)該鍋爐在結(jié)構(gòu)和形態(tài)上是否合理，安裝、操作和

2、檢修是否便利。是否合理，安裝、操作和檢修是否便利。 這種模型僅反映了其結(jié)構(gòu)特性，而未能反映這種模型僅反映了其結(jié)構(gòu)特性，而未能反映鍋爐內(nèi)部傳熱學、熱力學、流體力學的特性。鍋爐內(nèi)部傳熱學、熱力學、流體力學的特性。一、計算機仿真的概念第二節(jié)計算機仿真技術(shù)簡介第二節(jié)計算機仿真技術(shù)簡介 按真實系統(tǒng)的數(shù)學關(guān)系構(gòu)造系統(tǒng)的數(shù)按真實系統(tǒng)的數(shù)學關(guān)系構(gòu)造系統(tǒng)的數(shù)學模型，即將實際系統(tǒng)的運動規(guī)律用數(shù)學學模型，即將實際系統(tǒng)的運動規(guī)律用數(shù)學形式表達出來，并在數(shù)學模型上進行試驗形式表達出來，并在數(shù)學模型上進行試驗，再現(xiàn)系統(tǒng)的某些特性。，再現(xiàn)系統(tǒng)的某些特性。數(shù)學模型能精確反映系統(tǒng)內(nèi)部的各種靜數(shù)學模型能精確反映系統(tǒng)內(nèi)部的各種靜態(tài)

3、和動態(tài)特性，如鍋爐運行中的燃料化學反態(tài)和動態(tài)特性，如鍋爐運行中的燃料化學反應(yīng)、傳熱過程、能量儲存與釋放、工質(zhì)循環(huán)應(yīng)、傳熱過程、能量儲存與釋放、工質(zhì)循環(huán)流動的特性。流動的特性。數(shù)學仿真一般是利用計算機對系統(tǒng)數(shù)學數(shù)學仿真一般是利用計算機對系統(tǒng)數(shù)學模型進行運算和試驗的。因此，利用計算機模型進行運算和試驗的。因此，利用計算機實現(xiàn)的系統(tǒng)數(shù)學仿真也稱為：實現(xiàn)的系統(tǒng)數(shù)學仿真也稱為：計算機仿真計算機仿真。 將真實系統(tǒng)的一部分用數(shù)學模型將真實系統(tǒng)的一部分用數(shù)學模型描述，并放到計算機上運行，而另一描述，并放到計算機上運行，而另一部分則構(gòu)造其物理模型或直接采用實部分則構(gòu)造其物理模型或直接采用實物，然后將它們連接成系

4、統(tǒng)，并在此物，然后將它們連接成系統(tǒng)，并在此系統(tǒng)上進行試驗，再現(xiàn)系統(tǒng)的某些特系統(tǒng)上進行試驗，再現(xiàn)系統(tǒng)的某些特性。性。這種仿真又稱：這種仿真又稱：半物理仿真。半物理仿真。4. 計算機仿真計算機仿真定義：定義：將一個能近似描述實際系統(tǒng)的數(shù)學模型進將一個能近似描述實際系統(tǒng)的數(shù)學模型進行第二次模化，轉(zhuǎn)變?yōu)橐粋€仿真模型，并將其放行第二次?；D(zhuǎn)變?yōu)橐粋€仿真模型，并將其放到計算機上進行運行，以再現(xiàn)系統(tǒng)某些特性的過到計算機上進行運行，以再現(xiàn)系統(tǒng)某些特性的過程。程。第一次?；貉芯繉嶋H系統(tǒng)與數(shù)學模型之間的關(guān)系；第一次?；貉芯繉嶋H系統(tǒng)與數(shù)學模型之間的關(guān)系；第二次模化：研究數(shù)學模型與計算機之間的關(guān)系。第二次模化：

5、研究數(shù)學模型與計算機之間的關(guān)系。計算機仿真是本課程討論的主題。計算機仿真是本課程討論的主題。5.計算機仿真系統(tǒng)計算機仿真系統(tǒng)用于系統(tǒng)仿真的一套軟、硬設(shè)備。用于系統(tǒng)仿真的一套軟、硬設(shè)備。其主體是計算機。其主體是計算機。按計算機的類型分，有以下仿真系統(tǒng)：按計算機的類型分，有以下仿真系統(tǒng)：1 1、模擬仿真系統(tǒng)、模擬仿真系統(tǒng)2 2、數(shù)字仿真系統(tǒng)、數(shù)字仿真系統(tǒng)3 3、數(shù)模混合仿真系統(tǒng)、數(shù)?；旌戏抡嫦到y(tǒng)用模擬計算機組成的仿真系統(tǒng)用模擬計算機組成的仿真系統(tǒng)。優(yōu)點：優(yōu)點：運算是運算是 “ “并行并行”的、的、 “ “連續(xù)連續(xù)”的的。因此，運算速度快，且更接近連。因此，運算速度快，且更接近連續(xù)系統(tǒng)。續(xù)系統(tǒng)。缺點

6、：缺點：運算精度低、線路復(fù)雜、對采樣和邏運算精度低、線路復(fù)雜、對采樣和邏輯系統(tǒng)的仿真比較困難、仿真的自動輯系統(tǒng)的仿真比較困難、仿真的自動化程度低（依靠排題板接線）?；潭鹊停ㄒ揽颗蓬}板接線）。用數(shù)字計算機組成的仿真系統(tǒng)用數(shù)字計算機組成的仿真系統(tǒng)。優(yōu)點：優(yōu)點：被仿真的系統(tǒng)包含在一組程序中，仿真自被仿真的系統(tǒng)包含在一組程序中，仿真自動化程度高，使用方便；運算精度高；易實現(xiàn)邏動化程度高，使用方便；運算精度高；易實現(xiàn)邏輯處理和非線性環(huán)節(jié)；程序和參數(shù)修改容易。輯處理和非線性環(huán)節(jié)；程序和參數(shù)修改容易。缺點：缺點：運算過程是運算過程是“串行串行”的，運算速度相對較的，運算速度相對較低、實時仿真和尋優(yōu)計算等不

7、如模擬仿真系統(tǒng)快低、實時仿真和尋優(yōu)計算等不如模擬仿真系統(tǒng)快。 用微型數(shù)字計算機陣列組成的仿真系統(tǒng)，用微型數(shù)字計算機陣列組成的仿真系統(tǒng)， 稱為稱為 用模擬計算機和數(shù)字計算機組成用模擬計算機和數(shù)字計算機組成的仿真系統(tǒng)。的仿真系統(tǒng)。由于模擬計算機和數(shù)字計由于模擬計算機和數(shù)字計算機的優(yōu)缺點是互補的。因此，該系統(tǒng)達算機的優(yōu)缺點是互補的。因此，該系統(tǒng)達到了揚長避短的目的。到了揚長避短的目的。該系統(tǒng)適用于：該系統(tǒng)適用于：（1 1）要求與實物連接，又有許多復(fù)雜函數(shù)需計算的實時）要求與實物連接，又有許多復(fù)雜函數(shù)需計算的實時仿真；仿真；（2 2）需要進行反復(fù)迭代計算（如統(tǒng)計分析、參數(shù)尋優(yōu)）需要進行反復(fù)迭代計算（如

8、統(tǒng)計分析、參數(shù)尋優(yōu)）的仿真；的仿真；（3 3）對計算機控制系統(tǒng)的仿真；）對計算機控制系統(tǒng)的仿真；三、火電廠仿真技術(shù)的應(yīng)用三、火電廠仿真技術(shù)的應(yīng)用 火電廠計算機仿真技術(shù)是一門綜火電廠計算機仿真技術(shù)是一門綜合技術(shù)。它是以數(shù)學、物理、化學合技術(shù)。它是以數(shù)學、物理、化學、傳熱學、熱力學、流體力學、控、傳熱學、熱力學、流體力學、控制理論、計算機技術(shù)、熱能動力、制理論、計算機技術(shù)、熱能動力、電工學、熱工儀表及電氣儀表等多電工學、熱工儀表及電氣儀表等多學科專業(yè)的理論為基礎(chǔ)，以計算機學科專業(yè)的理論為基礎(chǔ)，以計算機和各種物理效應(yīng)設(shè)備（它再現(xiàn)真實和各種物理效應(yīng)設(shè)備（它再現(xiàn)真實環(huán)境）為工具，對真實火電廠發(fā)電環(huán)境）為

9、工具，對真實火電廠發(fā)電機組及其運行進行仿真的技術(shù)。機組及其運行進行仿真的技術(shù)。 火力發(fā)電機組的設(shè)備龐大，系統(tǒng)復(fù)雜?；鹆Πl(fā)電機組的設(shè)備龐大，系統(tǒng)復(fù)雜。因此：因此：對設(shè)備的研究、對系統(tǒng)的試驗、對參對設(shè)備的研究、對系統(tǒng)的試驗、對參數(shù)的校正、對運行安全和經(jīng)濟的分析、對運數(shù)的校正、對運行安全和經(jīng)濟的分析、對運行值班員的培訓等，行值班員的培訓等，在實際發(fā)電機組上直接在實際發(fā)電機組上直接進行，或是很困難或根本不可能。利用仿真進行，或是很困難或根本不可能。利用仿真技術(shù)的特點、在模型上進行試驗研究，在仿技術(shù)的特點、在模型上進行試驗研究，在仿真機上培訓運行值班員，則上述問題可迎刃真機上培訓運行值班員，則上述問題可

10、迎刃而解。而解。 等之分。等之分。數(shù)學模型：數(shù)學模型：是系統(tǒng)的數(shù)學描述，是系是系統(tǒng)的數(shù)學描述，是系統(tǒng)研究的基礎(chǔ)，是計算機仿真的依據(jù)。統(tǒng)研究的基礎(chǔ)，是計算機仿真的依據(jù)。實物模型：實物模型：根據(jù)相似性建立的形象模型，具有實體性屬物理模擬試驗技術(shù)范疇。（不介紹）數(shù)學模型：數(shù)學模型：用符號和數(shù)學方程式表示的系統(tǒng)模型，具有抽象性。（本課程內(nèi)容）系系 統(tǒng)統(tǒng) 模模 型型按描述的狀態(tài)分按描述的方式分按系統(tǒng)的性質(zhì)分按求解的方法分按獲取的方法分靜態(tài)模型動態(tài)模型連續(xù)模型離散模型線性模型非線性模型分析求解模型數(shù)字求解模型理論模型黑箱模型反映變量間變量間的相互函數(shù)關(guān)系反映變量與時間變量與時間的函數(shù)關(guān)系是時間的連續(xù)函數(shù)

11、是時間的離散函數(shù)可以用線性微分方程描述不能用線性微分方程描述利于解析法求解 利于計算機求解 理論推導所得的模型試驗研究所得的模型 注：（1）一般，仿真所要建立的是系統(tǒng)的動態(tài)模型，靜態(tài)模型包含在動態(tài)模型之中。（2）電廠的系統(tǒng)仿真，往往綜合采用連續(xù)離散線性非線性理論黑箱等模型，以保證模型精度。線性特性汽機功率調(diào)節(jié)級壓力抽氣量汽機功率(凝汽式機組) 非線性特性爐膛傳熱溫度流量壓差轉(zhuǎn)速油壓 可以應(yīng)用疊加原理不能應(yīng)用疊加原理必須有效地協(xié)調(diào)線性與非線性之間的相互關(guān)系必要時,非線性特性線性化 工質(zhì)狀態(tài)工質(zhì)狀態(tài) （溫度） 工況工況環(huán)境環(huán)境條件條件 （過熱器管長） 時間時間建模時多采建模時多采用用分區(qū)集總分區(qū)集

12、總方法方法，即將，即將三維空間三維空間的分布參數(shù)簡化的分布參數(shù)簡化為一維空間。否為一維空間。否則無法求解。則無法求解。汽機甩負荷轉(zhuǎn)速煙溫主汽溫時間常數(shù)小響應(yīng)快燃料汽壓減溫水量主汽溫時間常數(shù)大響應(yīng)慢建模時應(yīng)處理好模型的簡潔與精確之間的關(guān)系對響應(yīng)速度差異很大的對象，一般分別建對響應(yīng)速度差異很大的對象，一般分別建立動態(tài)數(shù)學模型，可采用不同的時間步長立動態(tài)數(shù)學模型，可采用不同的時間步長進行仿真計算。進行仿真計算。分解集合把實際系統(tǒng)分成若干子系統(tǒng)，分別建立其模型，然后綜合。把若干子系統(tǒng)視為一個整體，統(tǒng)一建立其模型。模型精細復(fù)雜模型簡潔粗糙因此，在建立復(fù)雜系統(tǒng)的動態(tài)模型時，應(yīng)善于對系統(tǒng)進行合理的分解，即按

13、不同的化學、物理過程（如燃燒、傳熱、流動、做功等），把系統(tǒng)分成許多子系統(tǒng)，這些子系統(tǒng)并不都是簡單地按設(shè)備的類型來劃分的，而根據(jù)系統(tǒng)的屬性來劃分更為合理。在分解系統(tǒng)時，正確處理各子系統(tǒng)間的關(guān)系是十分重要的，因此要正確地選擇系統(tǒng)中的狀態(tài)變量，充分保證信息的傳遞如實際過程那樣進行。為滿足實時運行要求，仿真機應(yīng)采用速度快、為滿足實時運行要求，仿真機應(yīng)采用速度快、容量大的計算機。容量大的計算機。 （7 7）模型計算涉及的數(shù)據(jù)量大）模型計算涉及的數(shù)據(jù)量大數(shù)學模型運算涉及的數(shù)學模型運算涉及的I/OI/O變量、中間變量、常變量、中間變量、常數(shù)等數(shù)量巨大。例如：數(shù)等數(shù)量巨大。例如：一臺采用常規(guī)儀表盤一臺采用常規(guī)

14、儀表盤/ /臺的臺的300300MWMW燃煤發(fā)電機組燃煤發(fā)電機組的全范圍仿真機有：的全范圍仿真機有：50005000多個多個I/OI/O變量；變量； 10000 10000多多個中間變量；個中間變量； 建模是綜合知識的體現(xiàn)，它涉及的知識廣泛。建模是綜合知識的體現(xiàn)，它涉及的知識廣泛。如：如：理論知識理論知識物理學化學數(shù)學熱力學傳熱學物理學化學數(shù)學熱力學傳熱學專業(yè)知識專業(yè)知識鍋爐原理汽輪機原理電機原理機械原理、控鍋爐原理汽輪機原理電機原理機械原理、控制理論制理論實踐知識實踐知識結(jié)構(gòu)知識運行知識試驗技術(shù)結(jié)構(gòu)知識運行知識試驗技術(shù)系統(tǒng)數(shù)學模型建立在每一個方框基礎(chǔ)上系統(tǒng)數(shù)學模型建立在每一個方框基礎(chǔ)上分解分

15、解模型。模型。準確的模型是取得正確仿真結(jié)果的基礎(chǔ) 系統(tǒng)簡化 與 的正確與否密切相關(guān) 假設(shè)條件 物理規(guī)律依賴于對系統(tǒng)內(nèi)部 化學規(guī)律 的熟悉程度和概括能力 運作機理 是指與熱平衡相關(guān)的系統(tǒng)。例如主燃料系統(tǒng)、風煙系統(tǒng)、爐膛、主蒸汽系統(tǒng)、汽輪機和抽汽系統(tǒng)、發(fā)電機、凝汽器、凝結(jié)水系統(tǒng)、給水系統(tǒng)、循環(huán)水系統(tǒng)、冷卻塔和和加熱器疏水系統(tǒng)等。ucdtudcdtudcyadtdyadtydadtydnnnnnnnnnnn12211101111自由項強迫項(2-1)式中：y=y（t）系統(tǒng)的輸出量 u=u（t）系統(tǒng)的輸入量 ai， cI系數(shù)dtdp pny+a1pn-1y+an-1py+any=c0pn-1u+c1p

16、n-2u+cn-1unjnjjjnjjnupcypa0101或：或：（其中a0=1）njjjnnjjjnpapcuy0101(2-2)()()()()()()()()()()()(2211222221212112121111tuytaytaytadtdytuytaytaytadtdytuytaytaytadtdynnnnnnnnnnnaij隨時間而變時變系統(tǒng)aij不隨時間變定常系統(tǒng)Ui(t)全部恒等于零齊次線性微分方程組Ui(t)有一個不為零非齊次線性微分方程組式中：yi 局部輸出量； ui 局部輸入量； aij, 系數(shù)（2-3）)(,1122tudtyddtyddtdyytfdtydnnnn

17、令 y1 y2 y3 yn （2-4）一種表達式為則一階微分方程組的另)(,211132221tuyyytfdtdyydtdyydtdyydtdynnnnnn式中：y1，y2，yn為自變量t的n個未知函數(shù)。求解式（2-4）相當于求解式（2-5）(2-5)()()()()()()(12110111sUcsUscsUscsYassYasYsasYsnnnnnnn（2-6）即： sn 替換 dn/dtn s為拉普拉斯算子 Y(s)替換y(t) Y(s)為輸出量y(t)的拉普拉斯變換 U(s)替換u(t) U(s)為輸入量u(t)的拉普拉斯變換 此式與式（2-2）比較可知，在初值為零的情況下，用算子P

18、表示的式子與用傳遞函數(shù)G(s)表示的式子在形式完全相同。nnnnnnnasasascscscsUsYsG11112110)()()(（2-7）將上式整理可得系統(tǒng)傳遞函數(shù)G（s）的規(guī)范式為：nnnnnnnnajajajcjcjcjcjUjYjG)()()()()()()()()(111122110（2-8）上述模型只描述了系統(tǒng)輸入與輸出之間的關(guān)系未描述系統(tǒng)內(nèi)部的情況稱為：外部模型仿 真在計算機上對系統(tǒng)的數(shù)學模型進行試驗必須在計算機上復(fù)現(xiàn)（實現(xiàn)）系統(tǒng)即復(fù)現(xiàn)系統(tǒng)輸入量復(fù)現(xiàn)系統(tǒng)輸出量復(fù)現(xiàn)系統(tǒng)的內(nèi)部變量又稱：狀態(tài)變量仿真要求將系統(tǒng)辨識或其它方法建立的外部模型轉(zhuǎn)化為：內(nèi)部模型內(nèi)部模型狀態(tài)空間模型狀態(tài)空間模

19、型狀態(tài)方程狀態(tài)方程控制理論中稱之為：實現(xiàn)問題實現(xiàn)問題傳遞函數(shù)G(s)U1（s）U2（s）UL（s）Y1（s）Y2（s）Ym（s）當引入狀態(tài)變量后，框圖變成了：狀態(tài)變量xi UYlnmU、Y、X均為時間的函數(shù)。luuuU21（2-9）myyyY21（2-10） m個輸出變量yI，輸出變量的集合可用輸出向量Y來表示；nxxxX21（2-11） n個狀態(tài)變量xI，狀態(tài)變量的集合可用狀態(tài)向量X來表示：),(),()(00ttUtXftX),(),()(00ttUtXgtY（2-12）（2-13）系統(tǒng)的狀態(tài)方程系 統(tǒng) 的 輸 出 方程式中：X（t0）初始狀態(tài)向量 U（t0,t）初始時刻t0到t的系統(tǒng)輸入

20、向量 f 表示一個單值函數(shù) g 也表示一個單值函數(shù)若用微分方程表示，則其表示式為：)(),()(tUtXFtX)(),()(tUtXGtY（2-14）（2-15）若用線性微分方程表示，則其表示式如下：)()()()()(tUtBtXtAtX)()()()()(tUtDtXtCtY（2-16）（2-17）系統(tǒng)數(shù)學模型的狀態(tài)方程規(guī)范式在狀態(tài)空間中： A（t）(nn階)系統(tǒng)矩陣（系數(shù)矩陣、狀態(tài)矩陣）B（t）(nl階)控制矩陣（驅(qū)動矩陣）C（t）(mn階)輸出矩陣D（t）(ml階)傳遞矩陣（傳輸矩陣）當A(t)、B(t)、C(t)、D(t)隨時間而變化時，系統(tǒng)稱為變系數(shù)系統(tǒng)（時變系統(tǒng)）；當A(t)、B

21、(t)、C(t)、D(t)不隨時間而變時，可用A、B、C、D來表示，系統(tǒng)則稱為常系數(shù)系統(tǒng)（定常系統(tǒng)）。D(t)U(t)項表示系統(tǒng)的輸入和輸出通過D(t)直接聯(lián)系的部分。在普通的控制系統(tǒng)中，通常D(t)=0，則式（2-17）變成：（2-18）)()()(tXtCtYa、微分方程的強迫項無無導數(shù)項時的狀態(tài)方程（規(guī)范式）)(11111tucyadtdyadtdadtydnnnnnnn在該情況下，式（2-1）中的微分方程為：(2-19）其中：nnnnnnndtydxdtydxxdtydxxdtdyxxyx/1112223121今引進n個狀態(tài)變量一階微分方程組nnndtydx )()(112112111

22、222111tucxaxaxaxatucyadtdyadtydadtydannnnnnnnnnnn把上述一階微分方程組寫成矩陣形式，可得：ucxxxaaaaxxxXnnnnnn1211212100101000010（2-20）XxxxYn0, 0, 0, 1 0, 0, 0, 1 21(2-21)121101000010aaaaAnnn100ncB0 , 0 , 0 , 1CCXYBuAXX（222）（223）顯然，狀態(tài)變量的初值為：)0()0()0()0()0()0(121nnyyyxxx.uydtdydtyddtyd54322233且：且：3)0(,2)0(, 1)0(yyy 試求該系統(tǒng)的

23、試求該系統(tǒng)的狀狀態(tài)方程，并將給出的初態(tài)方程，并將給出的初始條件用狀態(tài)變量的初值表示。始條件用狀態(tài)變量的初值表示。解：令22321dtydyxdtdyyxyx 則有uxxxxxxxx523432133221將上式寫成矩陣形式，uydtdydtyddtyd54322233 321001xxxy狀態(tài)變量的初值為：3)0()0(2)0()0(1)0()0(321yxyxyx 即可得狀態(tài)方程：uxxxxxx500234100010321321規(guī)范式（規(guī)范式（）b、微分方程強迫項有有導數(shù)項時的狀態(tài)方程（規(guī)范式、）ucdtudcdtudcyadtdyadtydadtydnnnnnnnnnnn12211101

24、111dtdp pny+a1pn-1y+an-1py+any=c0pn-1u+c1pn-2u+cn-1unjnjjjnjjnupcypa0101或：（其中a0=1）njjjnnjjjnpapcuy0101njjjnxpau0引進n個狀態(tài)變量：xpdtxdxxxpdtxdxxpxdtdxxxxxnnnnn111122223121即： pjx=xj+1 (j=0, 1, 2, , n-1)（2-24）則式（2-24）可寫為：xpaxpapxaxaxpaunnnnnjjjn022101001njnjjnxpaxau101njjjnnnuxaxxp（*）將式（2-24）代入式（2-2）可得：njnjn

25、jjjnjjnjjnxpapcypa01001（*）由式（*）和（*）可得：（2-25）BuAXX1010111njnjjjnjjnCXxcxpcy規(guī)范式（）121101000010aaaaAnnn其中：100B,0, 2, 1cccCnn)0()0()0()0()0()0(121nnxxxxxx規(guī)范式（）的應(yīng)用：是系統(tǒng)的數(shù)學模型可直接寫成狀態(tài)方程形式，且初始條件也是按各狀態(tài)變量的初值給出時，計算將十分方便。假定給出的微分方程組為：njnjjjnjjnupcypa00（2-30）ucupcupcyapyaypaypannnnnnn1101110即u和y的導數(shù)階次相等若令：nnnjjjnnnpx

26、ucyapucyapucyapucyap)()()()(1111100則有：ucyapxnnn若將式（2-30）兩邊同時取不定積分一次并令：122112001)()()(nnnnnpxucyapucyapucyap則有：ucyaxpxnnnn111并令：233113002)()()(nnnnnpxucyapucyapucyap則有：ucyaxpxnnnn2212同理有ucyaxpxjjjj1ucyapxxjjjj1或：則有：ucyapxxjjjj111)()()()(1122112001ucyaucyapucyapucyapxjjjjjjj即：（j=0，1，2，n）（2-31）（#）又因為在j

27、=0時，由式（#）有：00010ucyaxx 則在a0=1時，有：求初值用求初值用ucxy01uaccxaucucxaxucyaxxuaccxxaucucxaxucyaxxuaccxxaucucxaxucyaxxuaccxxaucucxaxucyaxxnnnnnnnnnnnnnnnnnnnn)()()()()()()()(0120111101111011111202312201232232101211101121121此時由式（#）可得： 由上述討論，可得式（2-30）在a0=1時的狀態(tài)方程和輸出方程如下：BuAXuaccaccaccaccxxxxaaaaxxxxXnnnnnnnnnn0101

28、202101121121121000100010001DuCXucxxxxynn01210, 0, 0, 1 （2-32）（2-33）規(guī)范式（規(guī)范式（）uBAXuccccxxxxaaaaxxxxXnnnnnnnn121121121121000100010001CXxxxxynn1210,0,0, 1 （2-32）（2-33）規(guī)范式（規(guī)范式（）jiijiijijucyax1)(1)(1)0()0()0(其中： ( j=1，2，n )udtduydtdydtyddtyd654322233式中：u為單位階躍函數(shù)，且初始條件為：0)0(u3)0(,2)0(, 1)0(yyy 試求該系統(tǒng)的狀態(tài)方程和輸出

29、方程，并將給出的初始條件用狀態(tài)變量的初值表示。解：根據(jù)題意，有：n=3， a0=1， a1=2， a2=3， a3=4，c0=0， c1=0， c2=5， c3=6可由規(guī)范式（），得系統(tǒng)的狀態(tài)方程和輸出方程：uxxxxxxX6500041030123213213210, 0, 0, 1 xxxy其中，各狀態(tài)變量由式（2-31）可求得：uydtdydtyducyaucyapucyapxydtdyucyaucyapxyucyax532)()()(2)()(222211002311002001jiijiijijucyax1)(1)(1)0()0()0(可求得：100343)0()0()0()0()0

30、()0()0()0()0()0()0(422)0()0()0()0()0()0()0(1)0()0()0()0(2212211003111002001ucyayayucyaucyaucyaxyayucyaucyaxyucyax 應(yīng)該指出：一個系統(tǒng)可以引入不同組合的狀態(tài)變量。那么，所設(shè)的狀態(tài)變量不同，獲得的狀態(tài)方程也會不同。即一個外部模型，可能有很多不同的內(nèi)部模型（狀態(tài)方程）一個系統(tǒng)的狀態(tài)方程形式并非是唯一的。)()()()(1122112001ucyaucyapucyapucyapxjjjjjjj（j=0，1，2，n）（2-31）由式（由式（2-31）可知，所設(shè)的各狀態(tài)變量與可知，所設(shè)的各狀態(tài)

31、變量與u的導數(shù)項有關(guān)，因此而得的狀態(tài)方程的導數(shù)項有關(guān)，因此而得的狀態(tài)方程:BuAXuaccaccaccaccxxxxaaaaxxxxXnnnnnnnnnn0101202101121121121000100010001DuCXucxxxxynn01210,0,0, 1（2-32）（2-33）規(guī)范式（規(guī)范式（）其解其解可能可能不是唯一的。不是唯一的。ucupcupcyapyaypaypnnnnnnn110111（2-34）如取以下如取以下n個變量作為一組狀態(tài)變量：個變量作為一組狀態(tài)變量：uxuuuuyxuxuuuyxuxuuyxuyxnnnnnnnn1112211012221031110201 （

32、3-35）且式中，且式中， 0， 1， 2， n分別為：分別為：a0=1011221103122133021122011100nnnnnnaaaacaaacaacacc（2-36）則可保證狀態(tài)方程解的唯一性。根據(jù)上面一組狀態(tài)變量，對于式（2-34）所描述的系統(tǒng)，可得出如下規(guī)范式：DuCXYBuAXX（2-37）（2-38）式中：00121121121121,0, 0, 0, 1 ,100001000010,cDCBxxxxXaaaaAxxxxXnnnnnnnnn狀態(tài)方程規(guī)規(guī)范范式式（）1對于以上各規(guī)范式：若A A，B B，C C，D D隨時間變化稱為：變系數(shù)系統(tǒng)變系數(shù)系統(tǒng)或時變系統(tǒng)時變系統(tǒng)若A

33、 A，B B，C C，D D不隨時間變稱為：常系數(shù)系統(tǒng)常系數(shù)系統(tǒng)或定常系統(tǒng)定常系統(tǒng)2一個系統(tǒng)的數(shù)學模型并非唯一的。表示一個系統(tǒng)的模型，可用不同的數(shù)學方法。最常用的方法為傳遞函數(shù)，其次是微分方程，狀態(tài)方程。3采用微分方程微分方程進行連續(xù)系統(tǒng)的仿真時：微分方程一階微分方程組代數(shù)方程采用變量分離法一般把轉(zhuǎn)換成轉(zhuǎn)換成 （a）應(yīng)根據(jù)初始值是否為零，選擇相應(yīng)的規(guī)范式。初始值為零時可選用規(guī)范式（）或規(guī)范式（）。 （b）若高階微分方程的初始條件是輸入u，輸出y及其導數(shù)的初值，應(yīng)變換成狀態(tài)方程的初值。 （c）若微分方程的強迫項含有導數(shù)項，可采用規(guī)范式（）或規(guī)范式（），但此時可能沒有唯一的解，若要獲得唯一的解，則

34、應(yīng)采用規(guī)范式（）。5、對已確定的系統(tǒng)，仿真所采用的數(shù)學方法可以不同，但研究結(jié)果應(yīng)與系統(tǒng)的實際情況一致。數(shù)值積分是數(shù)值分析的一個基本問題。數(shù)值積分是數(shù)值分析的一個基本問題。也是復(fù)雜計算問題中的一個基本組成部分。也是復(fù)雜計算問題中的一個基本組成部分。數(shù)值積分往往用極簡單的方法就能較好地得數(shù)值積分往往用極簡單的方法就能較好地得出對所求解的具體數(shù)值問題的解答。出對所求解的具體數(shù)值問題的解答。但數(shù)值積分的難點在于計算時間有時會過長但數(shù)值積分的難點在于計算時間有時會過長，有時會出現(xiàn)數(shù)值不穩(wěn)定現(xiàn)象。，有時會出現(xiàn)數(shù)值不穩(wěn)定現(xiàn)象。另外，數(shù)值積分的理論性較強。其理論和方另外，數(shù)值積分的理論性較強。其理論和方法都已

35、經(jīng)比較成熟，計算精度也比較高。法都已經(jīng)比較成熟，計算精度也比較高。3.1 仿真中研究數(shù)值積分法的意義仿真中研究數(shù)值積分法的意義數(shù)值解的一種近似方法。對于連續(xù)系統(tǒng)的高階微分方程，可化為若干個一階微分方程組成的方程組。數(shù)值積分法是求解微分方程：00)(),(ytyytfdtdyy 例如：下式所示的狀態(tài)方程BuAxx可以化為一個一階微分方程組ubxaxaxaxubxaxaxaxubxaxaxaxnnnnnnnnnnn2211222221212112121111所以，連續(xù)系統(tǒng)的仿真就是從給定的初始條件出發(fā)，對描述系統(tǒng)動態(tài)特性的常微分方程或常微分方程組進行求解，從而得到系統(tǒng)在一定輸入作用下的變化過程。在

36、求解這些微分方程時，最常用、也是最有效的一種方法就是數(shù)值積分法。3.2 數(shù)值積分法仿真的基本原理數(shù)值積分法仿真的基本原理對微分方程（3-1）兩端同時取積分，可得dyftytytt0),()()(0dyftytytt0),()()(0當 時，上式變?yōu)?：1nttdtytftytynttn10),()()(01將積分項拆成兩項 dtytfdtytftytynnnttttn10),(),()()(01)(nty故上式可寫為： dtytftytynnttnn1),()()(1此式是方程（3-1）在tn+1時刻的精確解。 在數(shù)值解法中，希望用近似解 ：nnnQyy1代替準確解，其中 ： 1ny為)(1n

37、tyny為)(ntynQ為1),(nnttdtytf的近似值 httnn1令：稱為計算步長或步距 是從常微分方程（是從常微分方程（3-1）出發(fā)建）出發(fā)建立的離散數(shù)學模型立的離散數(shù)學模型差分方程。差分方程。 dtytftytynnttnn1),()()(1nnnQyy11ny為)(1ntyny為)(ntynQ為1),(nnttdtytfnnnQyy11ny為)(1ntyny為)(ntynQ為1),(nnttdtytf由此可見，數(shù)值積分法就是在已知微分方程初值的情況下，求解該方程在一系列離散點 處的近似值，其特點是步進式根據(jù)初始值逐步遞推地計算出以后各時刻的值。從式（3-8）可知，數(shù)值積分法的主要

38、問題歸結(jié)為如何對f（t，y）進行數(shù)值積分求出f（t，y）在區(qū)間tn， tn+1上定積分的近似值Qn。采用不同的方法求Qn，就出現(xiàn)了各種各樣的數(shù)值積分方法。不同的數(shù)值積分將對求解的精度、速度和數(shù)值穩(wěn)定性會產(chǎn)生不同的影響，這將在下述內(nèi)容中具體介紹。 數(shù)值積分法種類繁多，在此從實用角度介紹幾種基本的方法 3.3 歐拉歐拉(Euler)法法3.3.1 簡單歐拉法簡單歐拉法 歐拉法是一種最簡單的數(shù)值積分法，對于方程：00)(),(ytyytfdtdyy 在區(qū)間tn， tn+1上求積分，得到： dtytftytynnttnn1),()()(1若區(qū)間tn， tn+1足夠小，則tn， tn+1上的f（t，y）

39、可近似地看成常數(shù)f（tn，yn） 。故可用矩形面積近似代替),(1ytfnntt即：tntn+1f（tn，yn）于是有：11),()(nnnnnyhytfyty將此式寫成差分方程為： , 3 , 2 , 1 , 0),(1nhytfyynnnn著名的歐拉公式著名的歐拉公式 3.3.2 改進的歐拉法改進的歐拉法 如果用梯形面積而不是矩形面積來代替每一個小區(qū)間上的曲線積分，就可以提高計算精度，梯形法的計算公式為： ),(),(2111nnnnnnytfytfhyy（3-11）式中的右端含有待求量yn+1，因而它是隱函數(shù)形式。這種方法不能自行啟動運算，需要依賴其它算法的幫助。 每次計算都用歐拉法算出

40、y（t n+1 ）的近似值 ，以此計算近似值 ，然后利用梯形公式（3-11）求出修正后的 。即有： 幫助方法：幫助方法：Pny1),(111PnnPnytfyCny1),(1nnnpnythfyy),(),(2111PnnPnnnCnytfytfhyy預(yù)估式 校正式 （3-12）改進的歐拉公式改進的歐拉公式 3.3.3 幾個基本概念幾個基本概念 簡單的歐拉法是用前一時刻tn的yn求出后一時刻的yn+1，這種方法稱為單步法，它是一種自行啟動的算法。如果求yn+1時需要tn , tn-1 , tn-2 時刻yn , yn-1 , yn-2 的值，這種方法為多步法（改進的歐拉法為兩步法），它是一種不

41、能自行啟動的算法。1、單步法與多步法、單步法與多步法 簡單的歐拉法表達式的右端，計算 所用的數(shù)據(jù)均已求出，這種公式稱為顯式公式。 改進的歐拉法表達式的右端，有待求量 ，這種公式稱為隱式公式。 隱式公式不能自行啟動，需要用預(yù)估-校正法。 單步法和顯式在計算上比多步法和隱式方便，但有時為了滿足精度、穩(wěn)定性等方面的要求，需要采用隱式算法。 2、顯式與隱式、顯式與隱式3、截斷誤差、截斷誤差 這里用泰勒級數(shù)為工具來分析數(shù)值積分公式的精度假定yn是精確的，用泰勒級數(shù)表示 處的精確解，即： 1nt)(0)(!)(! 2)()()(1)()2(2)1(1rnrrnnnnhtyrhtyhthytyty顯然，簡單

42、的歐拉法是從以上精確解中取前兩項之和來近似計算，每一步由這種方法引入的誤差稱為局部截斷誤差，簡稱截斷誤差。簡單的歐拉法的截斷誤差為： )(0)(211hytynn不同的數(shù)值方法有不同的截斷誤差。一般若截斷誤差為 ，則方法為r階的。所以方法的階數(shù)可以作為衡量方法精確度的一個重要標志。)(01rh4、舍入誤差、舍入誤差 由于計算機的字長有限，數(shù)字不能表示得完全精確，在計算過程中不可避免地會產(chǎn)生舍入誤差。舍入誤差與計算步長成反比。如果計算步長小，計算次數(shù)多，則舍入誤差就大。舍入誤差除了與計算機字長有關(guān)以外，還與計算機所使用的數(shù)字系統(tǒng)、數(shù)的運算次序以及計算所用的子程序的精度等因素有關(guān)。 采用數(shù)值積分法

43、求解穩(wěn)定的常微分方程，應(yīng)該保持原系統(tǒng)的穩(wěn)定特征。但是：（1）在計算機逐步計算時，初始數(shù)據(jù)的誤差及計算過程的舍入誤差對后面的計算結(jié)果將產(chǎn)生影響。（2）如果計算步長取的不合格，有可能使仿真出現(xiàn)不穩(wěn)定的結(jié)果。下面我們簡單討論一下這個問題。 差分方程的解與微分方程的解類似，可分為特解和通解兩部分。與穩(wěn)定性有關(guān)的是方程的通解，它取決于差分方程的特征根是否滿足穩(wěn)定性條件。例如，為了考查歐拉法的穩(wěn)定性，我們研究檢驗方程（Test Equation）：yy其中， 為方程的特征根 nnnnyhhyyy)1 (111 h即：2h表明：為使數(shù)字仿真穩(wěn)定，對計算步長應(yīng)有所限制表明：為使數(shù)字仿真穩(wěn)定，對計算步長應(yīng)有所限

44、制。 另外，穩(wěn)定性還與系統(tǒng)的特性以及數(shù)值積分方法有關(guān)。上述分析歐拉法穩(wěn)定性的思想，同樣適用于其它數(shù)值積分方法。 （3-14）（3-15）3.4 龍格龍格-庫塔庫塔（Runge-Kutta）法法由前面的分析可知，將泰勒展開式多取幾項以后截斷，就能提高截斷誤差的階數(shù)和計算精度。然而，直接采用泰勒展開方法要計算函數(shù)的高階導數(shù)，運用起來不便。龍格-庫塔方法的基本思想是：用幾個點上的函數(shù)值的線性組合代替函數(shù)的各階導數(shù)，然后按泰勒級數(shù)展開確定其中的系數(shù)，這樣既可避免計算高階導數(shù)，又可提高積分的精度。龍格-庫塔法有多種形式，以下從實用角度直接給出公式的形式和相應(yīng)的精度。3.4 .1 龍格龍格-庫塔方法的基本

45、思想庫塔方法的基本思想3.4.2 二階龍格二階龍格-庫塔方法庫塔方法2階龍格-庫塔方法的公式為：),(),()(2/(121211hkyhtfkytfkkkhyynnnnnn（3-16）上式表示的數(shù)值解是用的泰勒級數(shù)在2階導數(shù)以后截斷所求得的，因此稱為2階方法。故2階龍格-庫塔法與式改進的歐拉法相比，實質(zhì)完全相同。所以改進的歐拉法實質(zhì)上是2階龍格-庫塔法。)(03h截斷誤差為：實時仿真實時仿真的的2階龍格階龍格-庫塔方法庫塔方法)2,2(),(12121khyhtfkytfkhkyynnnnnn（3-17）)(03h其截斷誤差為：3.4.3 四階龍格四階龍格-庫塔方法庫塔方法4階龍格-庫塔法是

46、一種最常用的方法。其經(jīng)典表達式為：),()2,2()2,2(),()22(6342312143211hkyhtfkkhyhtfkkhyhtfkytfkkkkkhyynnnnnnnnnn（3-18）其截斷誤差為：)(05h顯然， 4階龍格-庫塔法的計算量較大，但計算精度較高，在比較不同算法的計算精度時，常以它的計算結(jié)果作為標準。實時仿真實時仿真的的4階龍格階龍格-庫塔方法庫塔方法)2103,54()52,53()52,52()5,5(),()105515(2441521413121543211khkhyhtfkhkkhyhtfkkhyhtfkkhyhtfkytfkkkkkkhyynnnnnnnn

47、nnnn（3-19）以上公式都是標量形式，如果要換成向量形式，只要把式中的標量y，f，k換成向量Y，F(xiàn)，K即可。從理論上講，可以構(gòu)造任意階數(shù)的龍格-庫塔方法，但是，精度的階數(shù)與計算函數(shù)值f 的次數(shù)之間并非等量增加的關(guān)系，見下表所列：表 3.1f 的計算次數(shù)與精度階數(shù)的關(guān)系每一步計算 f的次數(shù)234567N8精度階數(shù)234456n-2 由此可見，4階經(jīng)典龍格-庫塔方法有其優(yōu)越性，而4階以上的龍格-庫塔方法計算f所需的次數(shù)比階數(shù)多，增加了計算量，從而限制了更高龍格-庫塔方法的應(yīng)用。對于大量的實際問題，4階方法已可滿足精度要求，所以得到了廣泛的應(yīng)用。我們?nèi)圆捎脵z驗方程 進行討論，對它利用泰勒級數(shù)展開

48、得：3.4.4 龍格龍格-庫塔公式的穩(wěn)定區(qū)域庫塔公式的穩(wěn)定區(qū)域yy)(0!1)(11rinriinnhyrhyy對于 ，有 ，將它代入上式得：yyyyii)()(0)(!1)(! 211 121rnrnhyhrhhy（3-20）（3-21）令： 將其代入上式得到該式的穩(wěn)定條件穩(wěn)定條件為：hh1!1!21121rhrhh（ 3-22 ）由此穩(wěn)定條件，下表給出了各階龍格-庫塔公式的穩(wěn)定區(qū)域。表表3.2 龍格龍格-庫塔公式的穩(wěn)定區(qū)域庫塔公式的穩(wěn)定區(qū)域r1h所在的實穩(wěn)定區(qū)域1h1（2，0）22/12hh （2，0）36/2/132hhh（2.51，0）424/6/2/1432hhhh（2.78，0）在

49、使用龍格-庫塔公式時，選取的步長應(yīng)使落在穩(wěn)定區(qū)域內(nèi)。否則，在計算時會產(chǎn)生很大的誤差，從而得不到穩(wěn)定的數(shù)值解。例如：用4階龍格-庫塔方法解：1)0(,20yyy （1）用解析法求解：本例是穩(wěn)定的；（2）用數(shù)值法求解： 當h=0.1 時，數(shù)值解是穩(wěn)定的; 當h=0.2時，數(shù)值解就不穩(wěn)定了，這時因為：4)20(2 . 0hh此數(shù)值在穩(wěn)定區(qū)間以外，所以數(shù)值解不能收斂。這種對步長有限制的數(shù)值積分法稱為條件 穩(wěn)定積分法穩(wěn)定積分法。從4階龍格-庫塔方法的穩(wěn)定條件：中可以看出，系統(tǒng)的特征根越大，需要的積分步長就越小，這一點可作為選擇步長的依據(jù)。步長的大小，除了與數(shù)值積分方法的階數(shù)有關(guān)外，還與方程本身的性質(zhì)有關(guān)

50、。078. 2h除以上介紹的歐拉法、龍格-庫塔法外，數(shù)值積分方法還有許多種，如亞當姆斯方法、吉爾方法等等。此處不一一介紹。3.5 計算方法的選擇計算方法的選擇數(shù)值積分方法很多，在實際使用時存在一個選擇問題。對于一個具體問題，如何選擇具有一定的難度，至今尚無一種確定的方法。一般來說，數(shù)值積分方法的選擇應(yīng)考慮的因素有：1、精度要求、精度要求數(shù)值積分法的精度受截斷誤差舍入誤差積累誤差的影響，而這些誤差與積分方法、步長、計算時間、計算機精度等有關(guān)。一般：（1）積分方法的階數(shù)越高，截斷誤差越小，精度越高；（2）步長越小，精度越高；（3）多步法的精度高于單步法；（4）隱式算法的精度高于顯式算法。 因此，當

51、需要高精度時，可采用高階、多步、較小步長、隱式算法。 但是，步長的減小往往會增加迭代次數(shù)并增大舍入誤差和積累誤差。與此相反，在精度要求不高時，最好使用低階方法。計算速度取決于計算步數(shù)、每步積分所需的時間。而每步的計算時間又與積分方法有關(guān)，它主要取決于計算導數(shù)的次數(shù)（4階龍格-庫塔方法每步要計算4次導數(shù)），導數(shù)的計算是最費時的；為了加快計算速度，在積分方法已定的條件下，應(yīng)要保證精度的前提下，盡量選擇較大步長，以縮短積分時間。2、計算速度、計算速度數(shù)值解的穩(wěn)定性必須保證，否則，計算結(jié)果將失去真實意義。從穩(wěn)定性來看，不同的數(shù)值積分算法有不同的穩(wěn)定性。應(yīng)用時應(yīng)控制步長h，使數(shù)值積分算法在穩(wěn)定域內(nèi)。3、

52、數(shù)值穩(wěn)定性、數(shù)值穩(wěn)定性步長的選擇很重要，步長過大會增大截斷誤差，甚至出現(xiàn)數(shù)值不穩(wěn)定現(xiàn)象，過小了又因增加了步數(shù)，而使舍入誤差增大。所以，仿真的總誤差與步長的關(guān)系不是單調(diào)函數(shù)關(guān)系，而是一個具有極值的函數(shù)。如圖所示。3.6 計算步長的選擇計算步長的選擇誤差總誤差截斷誤差舍入誤差步長由圖可知，存在一個最佳步長。當積分方法確定以后，在選取步長時，需要考慮的一個重要因素就是仿真系統(tǒng)的動態(tài)響應(yīng)特性。如果系統(tǒng)的動態(tài)響應(yīng)快，導數(shù)變化激烈，則應(yīng)取高階方法和小步長進行計算。為了保證數(shù)值穩(wěn)定性，步長應(yīng)限制在最小時間常數(shù)（相當于最大特征值的倒數(shù)）的數(shù)量級上。用經(jīng)驗方法經(jīng)驗方法選取步長的兩種方法：（1）由系統(tǒng)方程中最小時

53、間常數(shù)Tmin來決定：min)05. 02 . 0(Th （2）由系統(tǒng)開環(huán)頻率特性的剪切頻率來決定：1)205(ch3.7 面向微分方程的仿真程序構(gòu)成面向微分方程的仿真程序構(gòu)成一般的仿真程序的組成可用下圖表示：圖中各方塊的功能為：主程序：主程序：實現(xiàn)對整個仿真計算的邏輯控制。輸入或預(yù)置參數(shù)模塊：輸入或預(yù)置參數(shù)模塊：輸入系統(tǒng)參數(shù)初值，計算步長、時間等參數(shù)。運行管理模塊：運行管理模塊：這是數(shù)字仿真程序的核心，對仿真計算進行時間控制，以保證計算機按要求進行計算及輸出。輸出及顯示模塊：輸出及顯示模塊：將仿真結(jié)果以數(shù)據(jù)或圖表的形式輸出給用戶。4階龍格-庫塔法仿真程序設(shè)計 目前，應(yīng)用于系統(tǒng)仿真的商用軟件包

54、諸多，但無論是開發(fā)者還是應(yīng)用者，了解程序的設(shè)計思想、分析程序的結(jié)構(gòu)與功能，完善和編寫仿真程序都是必要的。因此，此處通過對經(jīng)典龍格-庫塔方法（即式（7-22）仿真程序的介紹，使讀者了解和掌握編寫仿真程序的基本技術(shù)。 為了便于說明程序的實質(zhì)，我們用具體的例子來討論。 假設(shè)有一個2階系統(tǒng)： )( 1, 0)0(, 0)0(,44828. 222tuyyuydtdydtyddtdyyyy21,令：則：122214828. 24,yyudtdyydtdy以下是用C語言編寫的原理性程序：程序中變量說明：i，n: 分別為循環(huán)控制和積分器個數(shù)變量。h，t，tmax： 分別為步長、時間和仿真總時間變量。yn10

55、: 存放tn時刻y值的臨時變量數(shù)組。y10: 存放y值的變量數(shù)組。doty10: 存放 的數(shù)組。k110，k210，k310，k410:分別為存放K1，K2，K3，K4的數(shù)組。)( fy 程序清單 ：（“/*”和“* /”之間的文字是非執(zhí)行語句，作為注釋使用。 ）#includevoid diffeq();int i,n;float h, tmax, t, yn10,y10,doty10,k110,k210,k310,k410; main()/* 主函數(shù) */ n=2;/*設(shè)置積分器個數(shù)*/ y1=0;y2=0 /*設(shè)置積分器初值*/ h=0.01;tmax=20;t=0 /*設(shè)置步長、仿真總

56、時間和初始時間*/ whilt(t=tmax) /*判斷是否仿真時間到*/ printf(“t=%f”,t)；/*打印時間*/ printf(“y1=%fy2=%frn”,y1,y2); /*打印積分器輸出*/ for (i=1,i=n,i+) yni=yi; /*保留tn時刻的初值*/ dif(); /*計算 */ for(i=1;i=n;i+) k1i=dotyi; /*計算K1*/ yi=yni+(h/2)*k1i; )( fy t=t+h/2dif();/*計算 */for(i=1;i=n;i+) k2i=dotyi;/*計算K2*/ yi=yni+(h/2)*k2i; dif();

57、/*計算 */for(i=1;i=n;i+) k3i=dotyi;/*計算K3*/ yi=yni+(h/2)*k3i; t=t+h/2 dif();/*計算 */ )( fy )( fy )( fy for(i=1;i=n;i+) k4i=dotyi;/*計算K4*/ yi=yni+(h/6.0)*(k1i+2*k2i +2*k3i+k4 i ); /*計算tn+1時刻的y值*/ void dif() /*計算變量導數(shù)的函數(shù)*/float u; u=1;/*計算u=1(t)*/ doty1=y2; doty2=4*u-2.828*y2-4*y1;下面是程序的主函數(shù)main()的處理流程圖 四階

58、龍格四階龍格- -庫塔方法的程序處理流圖庫塔方法的程序處理流圖4.2 狀態(tài)方程的離散化 4.1 傳遞函數(shù)的離散化 本章介紹本章介紹:4.3 典型環(huán)節(jié)的離散化模型 4.4 常用非線性環(huán)節(jié)的仿真 4.5 離散相似法仿真的特點 )(sG保持器TTuy圖圖4-1 連續(xù)系統(tǒng)離散化結(jié)構(gòu)圖連續(xù)系統(tǒng)離散化結(jié)構(gòu)圖 smaxmax2s的條件時，可由采樣后的信號唯一地確定原始信號。 （2）把采樣后的離散信號通過一個低通濾波器，即可實現(xiàn)信號的重構(gòu)。常用的濾波器是零階保持器和一階保持器等。 因此，研究系統(tǒng)G(s) ，在一定的條件下可能等效地研究圖4-1所示的系統(tǒng)。 注意：圖4-1所示系統(tǒng)的采樣開關(guān)和保持器實際上是不存在

59、的，而是為了將(4-1)式離散化而虛構(gòu)的。)()()()()(sGsGZzUzYzGh（4-1）表表4-1 不同保持器的不同保持器的G(z) seTs1211)(seTTsTs221TseeTsTs)(22221111TssGZzzTsTssGZzzssGZzz)()()()()(ssG)(21TsTssG)(ssG)()(sGTTs1圖圖4-24-2若在積分環(huán)節(jié)之前再加一個采樣器和保持器，如下圖: )(sGTTs1圖圖4-34-3TseTs1)()()(sGZzTsseZsGZssGZTs111 )()(sGZzTzG（4-2）（4-3） 從上述可看出，虛擬采樣器和保持器可以不止一次地使用。

60、在必要的地方加入，可為求脈沖傳遞函數(shù)帶來方便。 但須注意，由于保持器的頻譜特性并非理想矩形，所以每加一次采樣器和保持器都會帶來誤差。因此(4-3)式較之表4-1中的式子誤差要大些，所以要盡量減少虛擬采樣器和保持器的使用。 asKsG)()()()()(aTaTezaeKassKZzzzG11所以得差分方程為 :111naTnaTnueaKyey)(也可根據(jù)(4-3)式來求G(z)： aTaTezKTezKzzTsGZzTzG)()(其對應(yīng)差分方程為 :11nnaTnKTuyey（4-4）（4-5）BuAxx 若人為地在系統(tǒng)的輸入端及輸出端加上采樣開關(guān)，同時為了使輸入信號復(fù)原為原來的信號，在輸入