控制系統(tǒng)計算機(jī)輔助設(shè)計與仿真

1、華北電力大學(xué)自動化系 劉長良 電話Email: changliang_,控制系統(tǒng)計算機(jī)輔助設(shè)計與仿真,教材,1、韓璞，自動控制系統(tǒng)數(shù)字仿真，中國電力出版社，1996 韓璞，控制系統(tǒng)數(shù)字仿真技術(shù)，中國電力出版社，2007 2、薛定宇，控制系統(tǒng)仿真與計算機(jī)輔助設(shè)計，機(jī)械工業(yè)出版社，2005,參考教材,1、薛定宇，控制系統(tǒng)計算機(jī)輔助設(shè)計，清華大學(xué)出版社，1996 2、孫增圻，控制系統(tǒng)計算機(jī)輔助設(shè)計與仿真，清華大學(xué)出版社，1988 3、劉金琨，先進(jìn)PID控制及其MATLAB仿真，電子工業(yè)出版社，2003 4、張葛祥，MATLAB 仿真技術(shù)與應(yīng)用，清華大學(xué)出版社，2003,一、

2、本課程的特點及要求,本課程是一門綜合應(yīng)用課程。,目的：在自動控制原理、過程控制的基礎(chǔ)上，進(jìn)一步掌握控制系統(tǒng)計算機(jī)輔助設(shè)計與仿真的有關(guān)理論與方法以及應(yīng)用這些理論與方法進(jìn)行控制系統(tǒng)研究的實用技能。 掌握：控制系統(tǒng)數(shù)字仿真、控制系統(tǒng)計算機(jī)輔助設(shè)計、以MATLAB語言為基礎(chǔ)的控制系統(tǒng)計算機(jī)輔助設(shè)計與仿真的程序設(shè)計方法以及這些方法在工程實際中的應(yīng)用。 培養(yǎng)：運用相關(guān)理論和計算機(jī)輔助手段解決實際問題的能力。,第一章 緒論,計算機(jī)輔助建模 計算機(jī)仿真 控制系統(tǒng)計算機(jī)輔助分析 控制系統(tǒng)計算機(jī)輔助設(shè)計,二、主要內(nèi)容,三、動態(tài)過程數(shù)學(xué)模型概述,1. 數(shù)學(xué)模型: 是對于現(xiàn)實世界的一個特定對象，根據(jù)其內(nèi)在規(guī)律，做出一

3、些必要的假設(shè)，運用適當(dāng)?shù)臄?shù)學(xué)工具，得到一個數(shù)學(xué)結(jié)構(gòu)。 數(shù)學(xué)模型又可分為靜態(tài)模型與動態(tài)模型 靜態(tài)模型：用來描述系統(tǒng)在穩(wěn)定狀態(tài)或平衡狀態(tài)下各種輸入變量與輸出變量之間的關(guān)系。 例如：當(dāng)機(jī)組運行在穩(wěn)定狀態(tài)時，輸入的物質(zhì)及能量保持不變，機(jī)組各系統(tǒng)的參數(shù)也將保持穩(wěn)定，這些穩(wěn)定工況下各參數(shù)之間的關(guān)系便可用靜態(tài)模型描述。靜態(tài)模型主要用于機(jī)組的設(shè)計計算及校核計算，一般要求具有較高的精度。,動態(tài)模型： 用來描述系統(tǒng)在過渡過程中各種變量隨時間變化的關(guān)系。 當(dāng)系統(tǒng)從一個穩(wěn)定狀態(tài)變化到另一穩(wěn)定狀態(tài)時，哪些參數(shù)會發(fā)生變化，其變化的速度及歷程如何，這些都屬于動態(tài)模型要研究的問題。 例如，當(dāng)燃料量變化時，機(jī)組原來的平衡狀態(tài)就

4、會受到破壞，電功率等參數(shù)都將發(fā)生變化，經(jīng)過一段時間運行，機(jī)組又將達(dá)到新的平衡狀態(tài)。這個動態(tài)過程中電功率的變化規(guī)律需要用動態(tài)模型描述。,2、 動態(tài)過程建模的一般方法 (1). 機(jī)理法白盒法 依據(jù)基本的物理定律 優(yōu)點：范圍廣、物理意義明確 缺點：精度差，與簡化程度有關(guān)，與過程的不確定性有關(guān)。 (2). 試驗建模黑盒法 系統(tǒng)辨識 時域法： 飛升曲線法 頻域法： 正弦信號、離散信號 相關(guān)分析法：隨機(jī)信號 參數(shù)估計法：最小二乘法、極大似然法等 優(yōu)點：與實際設(shè)備特性一致，模型精度高，包含了各種干擾模型 缺點：需要大量的人力物力，模型的應(yīng)用應(yīng)限制在試驗范圍內(nèi)。 (3). 灰盒法 ： 先建立機(jī)理模型，再進(jìn)行辨

5、識、修正模型結(jié)構(gòu)和數(shù)據(jù),3. 火電廠熱工過程建模與仿真,火電廠熱工過程特點：系統(tǒng)復(fù)雜、多變量、大遲延、非線性、時變性建立模型時要有針對性：(1) 要明確建模的目的及模型的用途；(2) 確定建模的方法(3) 對系統(tǒng)進(jìn)行合理的簡化(4) 確定模型的輸入、輸出(5) 建立模型(6) 分析、驗證、修正模型(7) 模型應(yīng)用,國電電力大連莊河電廠2X600MW機(jī)組,國華盤山電廠2X500MW機(jī)組,大唐國際盤山2X600MW機(jī)組,大唐國際潮州2X600MW機(jī)組,華潤滄州2X325MW機(jī)組,國投集團(tuán)300MW火電機(jī)組集控運行值班員技能大賽,秦山300MW核電機(jī)組全范圍仿真機(jī),大慶石油培訓(xùn)中心變電仿真,500M

6、W火電機(jī)組全仿真機(jī),200MW火電機(jī)組的全仿真機(jī),大唐國際寧德電廠2X600MW超臨界機(jī)組仿真機(jī),神六航天模擬器,火電廠設(shè)備,計算機(jī)接口,實際DCS操作,數(shù)學(xué)模型,計算機(jī)接口、網(wǎng)絡(luò),仿真DCS操作,實際電廠,仿真機(jī),機(jī)理法建模實例：彈性阻尼系統(tǒng)建模實例 問題：研究行使在不平路面上的汽車的顛簸情況。 方法：物理模型研究 過程：（1）采用類比方法得到相應(yīng)的概念模型 （2）建立數(shù)學(xué)模型 建立彈性阻尼系統(tǒng)模型過程： （i）分析 （ii）假定 K、B線性 （iii）建模：牛頓第二定律 得到彈性阻尼系統(tǒng)模型：,四、數(shù)字仿真的基本概念 數(shù)字仿真：求解數(shù)學(xué)模型數(shù)值解的過程 仿真分類：按模型類型分；按系統(tǒng)的連續(xù)

7、性分；按時間尺度分； 實時仿真、非實時仿真,彈性阻尼系統(tǒng)建模實例,（3）模型求解,（i）模型變換 轉(zhuǎn)換為狀態(tài)方程：,初始條件：x1(0)=0, x2(0)=0;,（ii）選擇數(shù)值計算方法得到仿真模型 采用歐拉法可寫出迭代公式如下：,iii）選擇編程語言，編程并進(jìn)行仿真研究 例如：Matlab、C、Basic語言等 （4）對仿真結(jié)果進(jìn)行分析并驗證模型 一般需要試驗數(shù)據(jù)的驗證，或采用實驗數(shù)據(jù)修正模型 （5）模型應(yīng)用，在確認(rèn)模型正確后，可應(yīng)用該模型對彈性阻尼系統(tǒng)進(jìn)行研究,物理原型概念模型（分析、類比得到） 數(shù)學(xué)模型（假設(shè)，采用物理定律得到） 仿真模型（采用一定的數(shù)值計算方法） 數(shù)字模型（程序代碼、通

8、過編程得到） 實用模型（多次的模型分析、檢驗、修 正） 模型應(yīng)用（通過模型研究物理系統(tǒng)的特性） 建模與仿真需要明確的問題： 明確建模目的 模型的輸入、輸出 模型的描述形式 建模對象機(jī)理及特性分析 編程語言的選擇 模型的驗證方法：靜態(tài)工況、動態(tài)過程,動態(tài)過程建模與仿真研究的一般步驟,仿真的應(yīng)用領(lǐng)域，發(fā)展史，火電站仿真 機(jī)的發(fā)展及現(xiàn)狀 航空航天 核電、火電 交通、運輸 通訊、電子 經(jīng)濟(jì)、人類 我?；痣姺抡婕夹g(shù)的發(fā)展,五、 計算機(jī)輔助分析 、數(shù)學(xué)模型之間的相互轉(zhuǎn)化 、控制系統(tǒng)穩(wěn)定性分析 頻域法：由W(S)分析系統(tǒng)的穩(wěn)定性， 勞斯判據(jù)、波特圖、奈奎斯特圖、 根軌跡等 時域法：狀態(tài)空間法。,六、 計算機(jī)

9、輔助設(shè)計 頻域法： 基于傳遞函數(shù)的控制器設(shè)計方法； 例如：串聯(lián)校正器的概念及設(shè)計方法； 極點配置設(shè)計方法； 時域法： 基于狀態(tài)空間模型的控制器設(shè)計方法； 線性二次型最優(yōu)調(diào)節(jié)器的設(shè)計方法； 觀測器的基本設(shè)計方法； 過程控制系統(tǒng)：PID控制器設(shè)計及參數(shù)整定方法；,機(jī)理法建模：實例4、椅子著地問題：4條腿長度相等的椅子，放在起伏不平的地面上，4條腿是否能同時著地？,假定地面為連續(xù)平面，則轉(zhuǎn)動椅 子在90內(nèi)，總能找到一點使四 條腿同時著地。 讓椅子原地旋轉(zhuǎn)，x為旋轉(zhuǎn)角度。 設(shè)f(x)為A、C兩腿距地面距離 之和，g(x)為B、D兩腿距地面距 離之和。,由于任意三腿總在一個平面上，所以任意地點三腿可同時

10、著地?？傆衒(x) , g(x)之一為0，即f(x)g(x)=0 數(shù)學(xué)描述： 已知：f(x) , g(x)為x的連續(xù)函數(shù)，已知g(0)=0,f(0)0 且 f(x)g(x)=0 求證：存在x0 使得g(x0)=0 證明：令h(x)=f(x) - g(x), 則h(0)=f(0) - g(0)=f(0)0 將椅子轉(zhuǎn)動90度，對角互換，則：由g(x)=0,f(x)0 可得：g(90)=f(0)0, f(90)=g(0)=0 即：h(90)= f(90)-g(90)0, 則由中值定理可知，必定存在x0使得h(x0)=0, 由f(x0)g(x0)=0, 可知f(x0)=g(x0)=0 問題得證。,第二

11、章、動態(tài)過程數(shù)學(xué)模型描述方法,、動態(tài)過程建模的一般方法 機(jī)理法、試驗法、兩者結(jié)合,二、熱工過程機(jī)理建模的基本依據(jù),基本假設(shè)：,采用集總參數(shù)法，忽略系統(tǒng)參數(shù)沿空間的分布情況，只考慮時間導(dǎo)數(shù)項 ； 假定煙氣、空氣為理想氣體，滿足理想氣體狀態(tài)定律； 各系統(tǒng)滿足基本的物理及熱力學(xué)定律，如質(zhì)量守恒、能量守恒、動量守恒、傳熱方程、熱力學(xué)狀態(tài)參數(shù)方程等。,1質(zhì)量守恒：,其中，V 系統(tǒng)容積，工質(zhì)密度，Win 入口質(zhì)量流量，Wout 出口質(zhì)量流量。,A. 建摸時假定： (1) 進(jìn)水管上游壓力為定值，流量W1和閥門開度u1為線性關(guān)系; (2) 忽略流體密度變化，假定水箱等截面，面積為F； (3) 水箱與大氣相通，

12、流量W2只與水箱液位有關(guān)；,C. 動態(tài)模型：,質(zhì)量守恒應(yīng)用1：求水箱水位,B. 模型輸入：u1，輸出：L,質(zhì)量守恒應(yīng)用2：求流體網(wǎng)絡(luò)的節(jié)點壓力,A. 假定： (1) 節(jié)點周圍的管道容積都集中在節(jié)點處，容積為V (2) 忽略節(jié)點容積V的變化,B. 模型輸入：P1，P2, P3 輸出：P,C. 動態(tài)模型：,壓力節(jié)點動態(tài)模型：,2能量守恒：,應(yīng)用：過熱器、省煤器等單項介質(zhì)換熱器,（1）以焓值作為集總參數(shù)：,集總參數(shù)法：,（2）以溫度作為集總參數(shù)：,3. 伯努力方程,忽略密度變化：,應(yīng)用：管路中流量與壓差的關(guān)系,V1 、V2 ：平均流速；Z1 、Z2 ：高度； P1 、P2 壓力；hw 管道上總的能量

13、損失；,忽略高度，速度變化：,4. 換熱方程:,煙氣-金屬:,工質(zhì)-金屬:,斯蒂芬波爾茲曼定律：,5熱力學(xué)參數(shù)及熱力學(xué)方程：P、T、H、S、V,水和水蒸汽模型 : HH(P、T), SS(P、T), VV(P、T), TT(P、H),計算方法：,（1）. 查表計算 （2）. 國際公式化委員會(IFC)于1967年提出的IFC-67公式; 國際水和水蒸氣協(xié)會（IAPWS）于1997年采納了新型的水和水蒸氣性質(zhì)計算模型IAPWS-IF97; （3）. 其他擬合公式,1. 高階微分方程 2. 響應(yīng)函數(shù)： 單位脈沖響應(yīng)，單位單位階躍響應(yīng)， 方波響應(yīng) 3. 傳遞函數(shù) 4. 脈沖傳遞函數(shù)（差分方程） 5.

14、 狀態(tài)方程描述,三、動態(tài)過程數(shù)學(xué)模型描述方法,狀態(tài)方程描述，矩陣形式， 包括三部分：狀態(tài)方程、輸出方程、初值 對于MIMO系統(tǒng)：,四、動態(tài)過程數(shù)學(xué)模型描述方法轉(zhuǎn)換 1. 高階微分方程狀態(tài)方程,設(shè)：x1=x, x2=x 可得：,輸出方程： x=x2 初值：x1(0)=0, x2(0)=0;,例1：,設(shè)：x1=y, x2=y, x3=y 則：x1=x2, x2=x3, x3=-4x1-3x2-2x3+5u,寫成矩陣形式如下,狀態(tài)方程：,輸出方程： y=x1,狀態(tài)變量初值： x1(0)=0,x2(0)=0,x3(0)=0,例2：,整理得：,方法2: 逐步積分法 兩邊積分：,設(shè)x1=5u-4y , 再

15、積分:,設(shè)x2=x1-3y+6u , 再積分：,設(shè)x3=x2-2y 得: y=x3,狀態(tài)變量初值： x1(0)=0,x2(0)=0;,對于通式的高級微分方程：,利用方法二可得:,輸出方程：y=xn，初值 x (i)=0;,2. 傳遞函數(shù)狀態(tài)方程 方法（1）：傳遞函數(shù) 微分方程狀態(tài)方程,例3：求下列傳遞函數(shù)對象的狀態(tài)方程描述,解：先將傳遞函數(shù)進(jìn)行拉氏反變換得到微分方程,再將微分方程轉(zhuǎn)換成狀態(tài)方程,例4：,2. 傳遞函數(shù)狀態(tài)方程 方法（2）：直接轉(zhuǎn)化 串聯(lián)法,等效框圖：,解：,同理：,狀態(tài)方程：,輸出方程：,初值：x1(0)=0; x2(0)=0; x3(0)=0;,2. 傳遞函數(shù)狀態(tài)方程 方法（

16、3）：直接轉(zhuǎn)化 并聯(lián)法,等效框圖：,并聯(lián)法的特點：方法通用，A為對角陣,同理：,狀態(tài)方程：,輸出方程：,初值：x1(0)=0; x2(0)=0; x3(0)=0;,注意： 1. 狀態(tài)方程描述不是唯一的，不同變化方法得到的狀態(tài)方程及輸出方程是不一樣的 ，但階次不變。 2. 要想得到A為對角陣，必須采用并聯(lián)法！,例如：前述的例題也可以采用并聯(lián)法求解,3. 傳遞函數(shù)框圖形式的描述狀態(tài)方程,例5：求如下控制系統(tǒng)的狀態(tài)方程描述,等效框圖：,按照前述方法，可得： x1=0.04e， x2=-0.1x2+0.1e1 x3=-0.05x3+0.05x2 e=u-x3， e1=0.05e+x1=0.05u-0.

17、05x3+x1,把e,e1代入上面各式，整理可得狀態(tài)方程描述:,狀態(tài)方程：,輸出方程： y=x3,初值：x1(0)=0; x2(0)=0; x3(0)=0;,本章復(fù)習(xí)題,2. 將微分方程轉(zhuǎn)為狀態(tài)方程描述，假定初值為0,3. 將傳遞函數(shù)轉(zhuǎn)為狀態(tài)方程描述，要求A為對角陣。,簡述熱工過程機(jī)理建模的基本假設(shè)，并給出常用的公式,4. 求以下控制系統(tǒng)狀態(tài)方程描述，輸入為u,輸出為y,假定初值為0.,第三章：數(shù)字仿真算法-數(shù)值積分法,問題：已知狀態(tài)方程，如何求狀態(tài)變量的數(shù)值解？,已知： 狀態(tài)方程： X=AX+Bu 輸出方程： y=Cx+Du 初值： X(0) 求：各個時刻輸出y的值 即：給定一個時間計算步距

18、DT，求y(DT), y(2DT), y(3DT), y(4DT),.y(kDT).y(nDT),例：彈性阻尼系統(tǒng)數(shù)值計算,設(shè)狀態(tài)變量：x1=x, x2=x 得：,X1=x2, x2=-k/m x1 b/m X2 +1/m u,如何計算x1,x2在各個時刻的數(shù)值？,對于：,在某個時刻t0, 近似有：,由此可得：,即：對于微小的DT，可以由t0時刻的x1,x2, 近似計算t0+DT時刻的x1,由前面，已知x1(0) ,x2(0), 這樣有如下遞推公式：,.,對于任意一點， 有如下遞推公式：,簡記為：,按此遞推公式即可求得各個時刻x1的數(shù)值,一、數(shù)值積分法仿真算法,1. 數(shù)值積分仿真算法原理,對于

19、： X=AX+Bu 設(shè) F(t)=A X + Bu，則：X=F(t),兩邊積分,得X計算通式：,取t=kDT, t=(k+1)DT,由高數(shù)可知： 在kT,(k+1)DT區(qū)間，存在一點，F(xiàn)*,可得：,如何計算積分項是數(shù)值積分法的關(guān)鍵問題！,如何計算F*是關(guān)鍵問題！,2. 各種數(shù)值積分仿真算法,（1）歐拉公式：,取 F*=F(kDT) , 得：,簡記：,注意到：,F(t)=A X + Bu,得：,由此，得到歐拉公式仿真算法：,問題： 在區(qū)間KDT，(k+1)DT內(nèi)F(t)是變化的，近似取F(t)=F(KDT)會造成較大的誤差。,歐拉公式幾何解釋：,精確解：,歐拉公式： 用矩形面積近似代替積分面積,

20、DT*F(K),誤差,分析：歐拉公式簡單實用，但當(dāng)DT較大時，會造成較大的誤差。要想提高精度，必須選擇較小的計算那步距DT。DT越小，精度越高。,DT的提高受計算機(jī)精度、速度的限制，在DT一定的情況下，如何提高計算精度？,問題：,（2）梯形公式：,用梯形面積近似代替積分面積,?。?由此，得到梯形公式仿真算法：,梯形公式存在的問題：,計算F(k+1)時需要用到X(k+1),無法使用！,怎么辦？,放棄,改進(jìn),如何改進(jìn)？,把歐拉公式引入梯形公式，先采用歐拉公式預(yù)測F(k+1),再采用梯形公式計算。,改進(jìn)方法：,改進(jìn)后的梯形公式：,預(yù)報校正公式,改進(jìn)的梯形公式幾何解釋,在DT一定的情況下，能否進(jìn)一步提

21、高計算精度？,問題：,歐拉公式：只采用k點的信息，計算F*,算法分析：,梯形公式：采用了k點和k+1點的信息，計算F*,是否可以采用k點，k+1,k+1/2點的信息計算F*,進(jìn)一步提高精度？,B、采用歐拉公式預(yù)測k+1點的F值F(k+1),C、計算F*:,D、采用歐拉公式計算，X(K+1)=X(k)+DT F*,梯形公式計算步驟：,A、先計算k點的F值F(k),（3）四階龍格庫塔法：,2. 根據(jù)F1預(yù)測k+1/2點的F值F(k+1/2), 記做F2,6. 根據(jù)F*, 采用歐拉公式計算，X(K+1)=X(k)+DT F*,計算步驟：,1、先計算k點的F值F(k)，記做F1,3. 根據(jù)F2重新預(yù)測

22、k+1/2點的F值F(k+1/2), 記做F3,5、計算F*:,4. 根據(jù)F3預(yù)測k+1點的F值F(k+1), 記做F4,計算步驟：,1.,2.,3.,4.,5.,6.,總體截斷誤差： 1、歐拉公式 O(DT) 2、梯形公式 O(DT2) 3、四階龍格-庫塔法 O(DT4),3. 仿真算法精度分析：,算法精度：,總體截斷誤差,局部截斷誤差,4. 仿真算法穩(wěn)定性分析,考慮較為簡單的情況，對于一階系統(tǒng),若采用歐拉公式計算，公式如下：,由k時刻誤差引起的k+1時刻的誤差計算如下：,若|1+DTa|1，注意一般a0，例如：1+DTa = -2，則：,由-1-2, a0, 得收斂區(qū)間：DT-2/a,反之

23、，若|1+DTa|1 則數(shù)值計算收斂,顯然，數(shù)值計算發(fā)散,精度與穩(wěn)定性,需要注意：1.精度與穩(wěn)定性是兩個不同的概念,兩者之間沒有 直接的聯(lián)系；2.上述每種數(shù)值積分算法都有自己的收斂區(qū)間；3.收斂區(qū)間與算法有關(guān)，也與對象特性有關(guān);4.仿真步距（計算步距）DT，越小，算法穩(wěn)定性越好；,5. 其它仿真算法,（1）三階龍格-庫塔法,以上數(shù)值計算方法為單步法，由 K 時刻的值推導(dǎo)K+1 時刻的值。,（2） 阿達(dá)姆斯法,多步法：,由 K，K-1，k-2.時刻的值推導(dǎo)K+1 時刻的值。,誤差 :O(T4),單步法:,阿達(dá)姆斯法：,（3）隱式歐拉公式,若取 F*=F(k+1) , 得：,對于高階系統(tǒng)需要求逆矩

24、陣，采用高級語言較為復(fù)雜，但對于單個微分方程：,把F(k+1) 代入得：,隱式歐拉公式如何求解？,隱式歐拉公式的優(yōu)點：,穩(wěn)定性好,對仿真步距DT沒有要求，不管DT取多大，計算不會發(fā)散,注意：雖然不會發(fā)散，但誤差可能很大！,說明：精度與穩(wěn)定性是兩個不同的概念,二、仿真程序設(shè)計,例：對于前述彈性阻尼系統(tǒng)模型,采用高級語言編程仿真，例如Matlab，C，Basic等,采用歐拉公式如下：,#include ； #include ； #include ； float K=0.7,B=0.8; X1 = 0, X2 = 0, M = 1, u = 1; int i,NP; float y1000,F1,F

25、2, ST=20,DT=0.1; main() NP = (int)(ST / DT);,C語言仿真程序設(shè)計,for(i=0;iNP;i+) F1 = X2; F2 = -K / M * X1 - B / M * X2 + u / M; X1 = X1 + DT * F1; X2 = X2 + DT * F2; yi = X1; printf(“i=%d, y(i)=%fn, i, Xi); plot(X,NP,ST); ,printf(Input Simulation Time ST=(0-100); scanf(%f,ST=20,DT=0.1; 也可以運行時輸入,彈性阻尼系統(tǒng)MATLAB程

26、序如下： clear; X1= 0; X2= 0; M= 1; K= 1; B=1;u =1.0; ST=20; DT=0.1; NP = ST / DT; for i=1:NP F1=X2; F2=-K/M*X1-B/M*X2+u/M; Y(i)=X1; X1 = X1 + DT * F1; X2 = X2 + DT * F2; end plot(1:NP)*DT,Y); title(彈性阻尼系統(tǒng));,clear; %清空內(nèi)存變量 A=0 1;-1 -1; B=0;1;C=1 0；%A B C X Y初值 X =0;0; u=1; % X 初值，u的大小 ST=20; DT=0.01; NP

27、=ST/DT； %仿真時間步距點數(shù) for i=1:NP; F=A*X+B*u; X=X+DT*F; Y(i)=C*X; end plot(1:NP)*DT,Y); title(彈性阻尼系統(tǒng));,矩陣運算形式MATLAB程序,歐拉公式:,彈性阻尼系統(tǒng)仿真曲線,仿真程序結(jié)構(gòu)：三個部分,1. 初始化部分,變量定義，狀態(tài)變量初值X1(0),X2(0)，系統(tǒng)參數(shù)M,K,B，仿真時間ST，仿真步距DT,2. 運算部分,先計算仿真點數(shù)NP，再循環(huán)計算X(K),3. 輸出部分,plot(1:NP)*DT,Y); title(彈性阻尼系統(tǒng));,（1）繪圖 （2）數(shù)據(jù)文件,仿真時間ST和仿真步距DT的選取,1.

28、 ST：根據(jù)經(jīng)驗估計完成動態(tài)過程所需要的時間，由此確定仿真時間；若對該過程沒有足夠認(rèn)識，可先任取一個仿真時間，例如，取ST=100 S，然后根據(jù)曲線的形狀確定ST是否需要調(diào)整。 2. DT：先取一個較小的數(shù)值，保證算法不發(fā)散，待ST確定之后，一般取DT=ST/100ST/1000即可。,程序設(shè)計中需要注意的問題,1. 變量、數(shù)組的定義 ，例如y(1000) 狀態(tài)變量的初值：X1= 0; X2= 0;,2. 變量的遞推關(guān)系如何實現(xiàn)？,X1 = X1 + DT * F1; X2 = X2 + DT * F2;,(1) k+1和k時刻的X采用同一個變量,(2) k+1和k時刻的X采用不同變量，另外給

29、出遞推關(guān)系,X1 = X10 + DT * F1; X2 = X20 + DT * F2; X10 = X1； X20 = X2；,梯形公式程序,for i=1:NP F10=X20; F20=-K/M*X10-B/M*X20+u/M; X11 = X10 + DT * F10; X21 = X20 + DT * F20; F11=X21; F21=-K/M*X11-B/M*X21+u/M; X1 = X10 + DT * (F10+F11)/2; X2 = X20 + DT * (F20+F21)/2; X10 = X1; X20 = X2; Y(i)=X1; end,MATLAB程序如下,

30、clear; %清空內(nèi)存變量 A=0 1;-1 -1; B=0;1;C=1 0；%A B C X Y初值 X =0;0; u=1; % X 初值，u的大小 ST=20; DT=0.01; NP=ST/DT； %仿真時間步距點數(shù) for i=1:NP; F0=A*X+B*u; X1=X+DT*F0; F1=A*X1+B*u; X=X+DT*(F0+F1)/2; Y(i)=C*X; end plot(1:NP)*DT,Y); title(彈性阻尼系統(tǒng));,矩陣運算形式MATLAB程序,梯形公式:,隱式歐拉公式程序,1.對于簡單的系統(tǒng)可直接推導(dǎo)模型,F1(K+1)= X2(K+1); F2(K+1)

31、= -K/M*X1(K+1) - B/M*X2 (k+1)+ u/M; X1(k+1)= X1(k) + DT*F1(k+1); X2(k+1)= X2(k) + DT*F2(k+1); 代入得： X1(k+1)= X1(k) + DT*X2(K+1); X2(k+1)=X2(k)+DT*(-K/M*X1(K+1)- B/M*X2(k+1)+u/M),由以上兩式可推得彈性阻尼系統(tǒng)隱式歐拉公式的計算模型如下：,編程是注意變量的設(shè)置，分清k時刻和k+1時刻的值,彈性阻尼系統(tǒng)MATLAB程序如下：隱式歐拉公式 clear; x10= 0; x20= 0; M= 1; K= 1; B=1;u =1.0

32、; ST=20; DT=0.1; NP = ST / DT; for i=1:NP x2=(x20-K*DT/M*x10+DT/M*u/(1+ K*DT*DT/M+B/M*DT); x1 = x10 + DT * x2 ; Y(i)=x1; x10=x1;x20=x2; end plot(1:NP)*DT,Y); title(彈性阻尼系統(tǒng));,clear; A=0 1;-1 -1; B=0;1;C=1 0；X =0;0; u=1; ST=20; DT=0.01; NP=ST/DT; I=eye(2); %單位陣 G= inv(I-DT*A); %求逆矩陣 for i=1:NP; X=G*X+D

33、T*G*B*u; Y(i)=C*X; end plot(1:NP)*DT,Y); title(彈性阻尼系統(tǒng));,矩陣運算形式MATLAB程序,隱式歐拉公式:,上機(jī)作業(yè)：,已知一個彈性阻尼系統(tǒng)，取M1，K1，B0.5 ，初始狀態(tài)為零，輸入為階躍擾動。,分別利用1、歐拉公式，2、梯形公式3、四階龍格庫塔法4、隱式歐拉公式，編程對該系統(tǒng)進(jìn)行仿真，改變步距T的大小，觀察各種仿真方法對T的敏感程度。,通過實例對上述各種仿真算法的穩(wěn)定性進(jìn)行分析,三、MATLAB仿真程序設(shè)計,1、MATLAB的M函數(shù)（7.1版）,（1）直接在Editor輸入程序，形成M文件，在Editor下直接運行。也可以存為M文件，在命

34、令行下運行。,MATLAB的命令窗口,把剛才的文件存為ABC.M文件，在命令行下輸入ABC運行。,Malab是解釋型語言，所有的語句在Matlab的命令窗口可直接運行。有些語句的調(diào)試可直接在命令窗口進(jìn)行。,函數(shù)格式如下： function 返回變量列表=函數(shù)名(輸入變量列表) 注釋說明段 函數(shù)體語句,（2）帶有函數(shù)名的M文件,可以與文件名稱一致，也可以不同。不同時以文件名為準(zhǔn)。,function Q=Simulation(ST,DT) % 彈性阻尼系統(tǒng)MATLAB程序：隱式歐拉公式 % 函數(shù)調(diào)用格式為ABC(ST,DT) x10= 0; x20= 0; M= 1; K= 1; B=1;u =1

35、.0; NP = ST / DT; for i=1:NP F1=X2; F2=-K/M*X1-B/M*X2+u/M; Y(i)=X1; X1 = X1 + DT * F1; X2 = X2 + DT * F2 end Q=x1; plot(1:NP)*DT,Y); title(彈性阻尼系統(tǒng));,在命令窗口輸入AB（20，0.1）可得到函數(shù)的返回值及曲線輸出。AB為文件名字AB.m，函數(shù)名最好與文件名一致。,帶有函數(shù)名稱的M函數(shù)，可以作為子函數(shù)被調(diào)用。,例如，可以在另一個程序中調(diào)用剛才的函數(shù)：,Untitled2.m ST=20; DT=0.1; p=AB(ST,DT);,運行上面程序，可得到同

36、樣的結(jié)果。,Matlab語言概述,1. 變量基本類型為矩陣型,A=1,2,3; 4,5,6; 7,8,9,特殊矩陣的實現(xiàn)：,a.零矩陣：A=zeros(m,m) mxn的零矩陣 b.單位陣：A=eye(n) nxn單位陣 c.隨機(jī)元素矩陣：A=rand(n) nxn隨機(jī)陣,例如： A=1,2; 3,4, B=5,6;7,8,A+B=6,8; 10,12,A-B=-4,-4; -4,-4,A*B=19,22; 43,50,A/B=3,-2; 2,-1,B-1=inv(B)=-4,3; 3.5,-2.5,A/B=A*B-1=3,-2; 2,-1,矩陣運算 +， ， X， /,注意：A，B的階次有一

37、定的要求,2. 變量名 區(qū)分大小寫,3. 常用的語句,（1）語句以, 或; 結(jié)束. 用,顯示運算結(jié)果，用;不顯示,（2）注釋語句 %,（3）循環(huán) for,例如：y=1+2+3+4+5+6y=0;for i=1:1:6 y=y+i;end,1 可缺省,包括6,結(jié)果：y=21,y=0; i=1; while( i=6) y=y+i; i=i+1; end,循環(huán) while,（4）if (條件） 語句塊 end,y=0; for i=1:1:6 if(i3), break; end y=y+i; end,結(jié)果：y=6,結(jié)果：y=21,（5）clear 清除工作空間的變量（刪去）,例如：clear a

38、 c 清除變量a,c,clear 清除所有變量,（6）輸入輸出語句,例如：ST=input(Enter ST=);,DT=input(Enter DT=);,先顯示 Enter ST=，然后等待輸入ST，再將輸入的數(shù)值賦值給ST,（7）用help在命令窗口可以得到幫助,例如：help AB 顯示： 彈性阻尼系統(tǒng)仿真程序 函數(shù)格式調(diào)用為ABC(ST,DT),（8）繪圖函數(shù),x=1,2,3,4,5,6; y=110,120,150,120,145,133;plot(x,y),注意：x, y數(shù)據(jù)個數(shù)必須相同,Plot( (1:NP)*DT,y) 等價于,T= (1:NP)*DT;Plot( T,y)

39、,for i=1:NP F1=X2; F2=-K/M*X1-B/M*X2+u/M; T(i)=i*DT, Y(i)=X1; X1 = X1 + DT * F1; X2 = X2 + DT * F2 end Q=x1; plot(T,Y); title(彈性阻尼系統(tǒng));,程序,（9）其它繪圖函數(shù),grid 網(wǎng)格,title(圖形的名稱);,hold on 鎖住屏幕,xlabel (時間 t/s),ylabel (位移 y/cm),四、熱工過程仿真,質(zhì)量守恒應(yīng)用1：求水箱水位,A. 假定：,B. 模型輸入：u1，輸出：L,C. 動態(tài)模型：,D. 仿真模型：,顯式歐拉公式,研究重點：系統(tǒng)從一個穩(wěn)定狀

40、態(tài)變化到另一個穩(wěn)定狀態(tài)時所經(jīng)歷的歷程。,E. 仿真研究,（1）確定系統(tǒng)初始穩(wěn)定狀態(tài)； （2）改變系統(tǒng)的輸入，觀察輸出的響應(yīng)過程,假定初始穩(wěn)態(tài)如下：L0=100 cm, F=900 cm2 , W10=W20=900 g/s, ro=1 g/cm3, u10=0.5 。,由此可計算得到動態(tài)模型中的系數(shù)： K1=1800, K2=90。,仿真時假定閥門開度由50%變到60%，相應(yīng)的matlab仿真程序如下：,clear; L=100; F=900;K1=1800; K2=90; RO=1; u1=0.6; ST = 2000; DT = 10;NP = ST / DT; for i=1:NP W1

41、=K1*u1; W2=K2*sqrt(L); Y(i)=L; L=L+DT*(W1-W2)/(RO*F); end plot(1:NP)*DT,Y); title(單容水箱仿真);,水箱水位仿真研究,水箱水位仿真曲線,DT仿真研究-顯式歐拉公式,水箱水位仿真研究-隱式歐拉公式,D. 仿真模型：隱式歐拉公式,非線性方程，無法直接求得L(k+1)的表達(dá)式,需要把上式中的非線性部分進(jìn)行線性化后，才能得到L(k+1)的表達(dá)式。,常用處理方法：,即：把原來非線性項分為兩部分：系數(shù)部分和線性部分。系數(shù)部分用K時刻的值近似計算，得到B2(k),整理得：,for i=1:NP B2=K2/sqrt(L); L

42、=( RO*F/DT*L+K1*u1)/(RO*F/DT+B2); Y(i)=L; end,不同仿真步距是輸出曲線： DT=10s, 500s,1000s,2000s,5000s,隱式歐拉公式仿真曲線,該系統(tǒng)由兩個水箱組成，是二階系統(tǒng)，所以其動態(tài)過程應(yīng)包括兩個動態(tài)方程。,對于水箱1、2分別列出其動態(tài)模型如下：,多容水箱系統(tǒng),例：如圖水箱系統(tǒng)：,顯式歐拉公式：,仿真模型：,仿真過程如下：,（1）確定初始平衡狀態(tài)； 假定L10=100 cm, F1=900 cm2 , W10=W20=W30=900 g/s, Ro=1 g/cm3, u10=0.5, L20=64 cm, F2=600 cm2 （

43、2）計算動態(tài)模型中的系統(tǒng)常數(shù)： 根據(jù)上面的初始條件，計算可得：K1=1800, K2=150, K3=112.5。 （3）假定調(diào)節(jié)閥開度u1由0.5變成0.6，研究動態(tài)過程,不同仿真步距是輸出曲線：顯式歐拉,隱式歐拉公式仿真模型如下：,按照同樣的線性化方法可得：,整理得仿真模型：,不同仿真步距是輸出曲線：隱式歐拉,如圖所示：流體網(wǎng)絡(luò)中一個壓力節(jié)點,建模時假定： (1) P點周圍的管道容積都集中在P處 (2) 忽略流體密度變化,根據(jù)質(zhì)量守衡可得：,二、 質(zhì)量守恒應(yīng)用之二：流體網(wǎng)絡(luò)壓力節(jié)點模型,對于熱工過程而言，一般情況下，溫度的變化速度要比壓力的變化速度慢得多，忽略溫度的變化,K1的大小反映了流

44、體的可壓縮性，K1值越大，說明流體可壓性越強(qiáng)，反之，可壓性則越小，對于不可壓流體K10。 需要注意的是，只有理想流體才可能是不可壓流體，實際流體都是可壓的。,忽略流體密度變化，則各路流量與壓差關(guān)系如下：,若采用顯式歐拉公式：,若采用隱式歐拉公式：,按前述方法對上式進(jìn)行線性化，設(shè)：,壓力節(jié)點仿真模型：,考慮流體網(wǎng)絡(luò)中各節(jié)點的計算順序，在計算P(k+1)時，P1(k+1)、P2(k+1)、P3(k+1)可能還是未知，為使模型通用，所有相臨節(jié)點壓力都取為K時刻的值，即：,如圖， 對于小流量，為簡單起見，直接引用流量,若采用隱式歐拉公式：,對于任意多支路的壓力節(jié)點，假定周圍有i個壓力節(jié)點，j路微小流量

45、：,建模時假設(shè)： （1）外壁絕熱； （2）采用集總參數(shù)法，假定金屬管道壁溫均勻一致； （3）忽略工質(zhì)定壓比熱Cp隨溫度的變化,三、 能量守恒應(yīng)用之一：絕熱蒸汽管道模型,根據(jù)質(zhì)量守衡：,根據(jù)能量守衡,對于工質(zhì)有：,對于金屬有：,工質(zhì)和金屬之間的換熱量：,對流換熱系數(shù)：,1. 基于溫度的集總參數(shù)模型,由于實際過程中壓力、流量變化較快，而溫度、焓值變化較慢，建立模型時一般把流動過程與傳熱過程分開考慮。質(zhì)量守衡動態(tài)過程在壓力節(jié)點模型或水箱水位模型中考慮，能量守衡方程在溫度、焓值模型中考慮。,取蒸汽側(cè)集總參數(shù)T=Tout ，采用顯式歐拉公式得：,若采用隱式的歐拉公式，注意到Q與T、Tm有關(guān)，為了增加模型

46、的穩(wěn)定性，需要將Q帶入,考慮到計算順序，前面公式中用到的Tm(k+1)要用Tm(k)代替才能計算,由于水蒸汽的定壓比熱隨壓力、溫度的變化很大，定值假定會造成較大誤差，在建模時常采用焓值作為熱平衡計算的基礎(chǔ)，建立基于焓值的集總參數(shù)模型:,H為集總參數(shù)，取H=hout。由于計算Q時要用到蒸汽溫度，還需要根據(jù)焓值計算溫度，然后再計算傳熱量Q。,假定蒸汽壓力為P,則：,2. 基于焓值的集總參數(shù)模型,由于Q與H沒有直接的函數(shù)關(guān)系，所以無論用顯式還是隱式歐拉公式，只能采用Q(k),Tm的遞推公式同前：,在該模型中，若考慮由于外壁與環(huán)境的自然對流換熱，則可假設(shè)總的散熱系數(shù)為Kamb，環(huán)境溫度為Tamb，模型

47、中蒸汽側(cè)不變，只需在金屬側(cè)減去管道對環(huán)境散熱量 ：,計算Tm的遞推公式也應(yīng)進(jìn)行相應(yīng)的改動。,建模時假定：,（1）整個換熱器用一段逆流管代替，沿管長方向吸熱均勻，金屬溫度均勻一致，為Tm ； （2）忽略煙氣定壓比熱隨溫度的變化； （3）忽略流動的動態(tài)過程，取Ws2=Ws1=Ws, Wg2=Wg1=Wg ，工質(zhì)壓力為Ps；,四. 能量守恒應(yīng)用之二：單項介質(zhì)換熱器模型,煙氣側(cè)：,動態(tài)模型如下：,工質(zhì)側(cè)：,金屬側(cè)：,換熱計算中蒸汽、煙氣取平均溫度，采用隱式歐拉公式：,整理得：,輸入：Ws, hs1, Ps, Wg, Tg1 系數(shù)：Kg1, Kg2, Ks1, Ks2, Mg ,Cpg, Mm ,Cm,

48、 Ms輸出：hs2, Ts2, Tg2, Tm,整理后模型的一般形式,整理好的仿真模型,整理好的仿真模型如下,混合聯(lián)箱：聯(lián)箱內(nèi)有一電加熱器，熱功率為Q,建模時假定： （1）忽略流體密度變化，水箱等截面，面積為F； （2）集總參數(shù)法，整個聯(lián)箱內(nèi)工質(zhì)溫度均勻一致； （3）忽略工質(zhì)定壓比熱Cp隨溫度的變化；,五、 質(zhì)量守恒-能量守恒同時應(yīng)用：混合聯(lián)箱模型,質(zhì)量守恒：,能量守恒：,雖然上式中V與L有關(guān)，但建立動態(tài)模型時一般并不將上式按微分式子展開，而是處理（2）式先將L按常數(shù)處理，得到仿真模型，再和（1）一起迭代運算。,動態(tài)模型：,微分項的近似處理方法：,簡化如下：,集總參數(shù)T取為出口溫度t3，采用隱

49、式歐拉公式：,仍按前面的線性化處理方法可得隱式歐拉公式模型:,模型輸入: W1，t1, W2, t2, Q,模型輸出: L，t3,已知：,第四章、離散相似法,一、離散相似法原理,狀態(tài)方程：,輸出方程：,初值：,求：Y的數(shù)值解,方法一：前述的數(shù)值積分法,方法二：離散相似法,通解：,K時刻：,K+1：,如何計算eAT及積分是該算法的關(guān)鍵！,遞推關(guān)系：,與數(shù)值積分法比較,數(shù)值積分法：近似處理F,離散相似法：近似處理,取t1=(K+1)T-t,進(jìn)行積分變換：,零階保持器離散相似算法如下：,在區(qū)間kT,(k+1)T, 近似取U(t)=U(kT),零階保持器,（1）、計算(T) ：,一般指數(shù) ：,矩陣指數(shù)

50、 ：,近似取前幾項,（2）、計算m(T) ：,A、當(dāng)系統(tǒng)為一階系統(tǒng)時，A變?yōu)槌?shù)，求(T)、 m(T)變得較簡單； B、 MATLAB有相應(yīng)的函數(shù),(T)、 m(T)的計算很復(fù)雜！,在區(qū)間kT,(k+1)T,近似取,代入得：,整理得：,其中：,三角階保持器,三角保持器離散相似算法如下：,方法三、通過拉氏變化求解：,求拉氏變換： sX(s)-X(0)=AX(s)+BU(s) (sI-A)X(s)=X(0)+BU(s) X(s)=(sI-A)-1X(0)+ (sI-A)-1BU(s) 求拉氏反變換 X(t)=L-1(sI-A)-1 X(0)+L-1(sI-A)-1BU(s),根據(jù)卷積定理：,對于：

51、,設(shè)：,由此可得：,與離散相似法比較：,可得：,即：可以通過拉氏反變換計算,同理：x2=-1/5x2+1/5x1 e=u-x2 ,整理得：,解: (1) 先求系統(tǒng)的狀態(tài)方程 設(shè)定狀態(tài)變量x1,x2,例1：對如下系統(tǒng)進(jìn)行仿真,x1=-1/4x1+1/4e,求(T),m(T), n(T),由于不是對角陣，求逆矩陣比較困難!,若A是對角陣，例如,可見，離散相似法過于繁瑣!,解決方法1：,系統(tǒng)為一階系統(tǒng)時，A=a, eAT=eaT,對每個一階系統(tǒng)分別進(jìn)行離散相似,x1=-1/4x1+1/4e,A=-1/4, B=1/4, (t)=e-t/4, (T)=e-T/4,x1(K+1)= e-T/4x1(K)

52、+ (1-e-T/4)e(K) x2(K+1)= e-T/5x2(K)+ (1-e-T/5)x1(K),同理：,e(K)=u(K)-x2(K) x1(K+1)= e-T/4x1(K)+ (1-e-T/4)e(K) x2(K+1)= e-T/5x2(K)+ (1-e-T/5)x1(K) y(K+1)= x2(K+1),仿真模型如下：,（1）原來是對u進(jìn)行近似處理，這次分別對兩個子系統(tǒng)的輸入e,x1進(jìn)行了近似處理，誤差增大。 （2）當(dāng)u是階躍輸入時，對u進(jìn)行近似處理不會引入誤差，而第二種方法對e,x1近似處理，引入了兩次誤差。,注意：,分析發(fā)現(xiàn)，計算x1(k+1)用到e(K+1) e(K+1)=u

53、(K+1)-x2(K+1) 可見，該模型是無法使用的。,為了減小誤差，采用三角保持器：,即：,A在區(qū)間kT,(k+1)T,整理得：,處理方法,方法1：對系統(tǒng)1（與反饋環(huán)節(jié)相鄰）使用零階保持器,方法2：改進(jìn)三角保持器,一階保持器,B在區(qū)間kT,(k+1)T,取,這樣近似處理后，算法實用化。,滯后一拍三角保持器,整理得：,二、保持器的概念,b. 一階保持器,c. 三角保持器,d. 滯后一拍的三角保持器,采樣系統(tǒng)：采樣開關(guān)、保持器 采樣處理離散，信號恢復(fù)（保持）相似,2保持器種類及頻率特性,1. 離散相似概念的來源：,a. 零階保持器,3. 各種離散相似處理后信號的變化,（1）.多次離散后造成信號有

54、較大的相位滯后,（2）.多次離散后造成信號有一定的幅值衰減,為了提高精度，對離散相似信號補償，包括對幅值及相位的補償,為了簡單,一般取,4. 補償器,補償器：,即:輸入信號超前一拍,例如對于零階保持器,加入補償環(huán)節(jié)后:,小結(jié)：,1. 對u離散相似，各種算法都能使用，但需要求得整個系統(tǒng)的狀態(tài)方程； 2對于負(fù)反饋系統(tǒng)，對e離散相似，三角保持器不能使用，對x1離散相似，各種算法都能使用； 3. 對e離散相似, 不能使用補償器； 4. 零階保持器可以無失真的再現(xiàn)階躍信號； 5. 其它保持器可以無失真的再現(xiàn)斜波信號； 6. 對各個子系統(tǒng)分別進(jìn)行離散相似可以簡化仿真算法，但會引入新的誤差；,離散相似法程序

55、設(shè)計,采用matlab 函數(shù)進(jìn)行彈性阻尼系統(tǒng)仿真的 程序如下(離散相似法零階保持器) clear all；clc； m=1;b=0.5;k=1;u=1;X=0 0;T=0.01;ST=30; A=0 1;-k/m -b/m;B=0 1/m;NP=round(ST/T); Q=expm(A*T); E=inv(symsub(symmul(s,eye(2,2),A); F=ilaplace(E) G=symmul(F,B),Qm=numeric(int(G,t,0,0.01) for i=1:NP y(i)=X(1); X=Q*X+Qm*u; end plot(0:NP-1)*T,y,y);,對于

56、單容水箱系統(tǒng)：,需要注意，上述算法是對于線性系統(tǒng)推導(dǎo)得出的，只能用于線性系統(tǒng)。,A=-B2/(roF)， B=K1/(roF) (T)=eAT,m(T)=K1/B2(1-eAT) L(K+1)= (T)L(K)+ m(T)u(K),clear; L = 100; F=900;K1=1800; K2=90; RO=1;u1=0.6; ST = 2000; DT =100; NP = round(ST / DT); for i=1:NP W1=K1*u1; B2=K2/sqrt(L); AA=-B2/(RO*F); BB=K1/(RO*F); FA=exp(AA*DT); FAM=K1/B2*(1

57、-FA); Y(i)=L; L=FA*L+FAM*u1; end plot(1:NP)*DT,Y,r); title(采用離散相似法進(jìn)行單容水箱仿真);,仿真結(jié)果分析： 離散相似法與 隱式歐拉公式對比,三、Z變換法,原理：脈沖傳遞函數(shù)直接對應(yīng)差分方程,1. 脈沖傳遞函數(shù)的求?。狠斎胼敵龅膠變換之比,拉普拉斯變換定義:,：拉普拉斯算子，,z變換定義：,已知X(t)的采樣信號x*(t),對于k0有x(kT)=0，,對上式進(jìn)行拉普拉斯變換可得：,定義：,得到：,即：z變換是離散拉普拉斯變換。,例1：求x(t)=1(t)的z變換。,例2求x(t)=t的z變換,變換的性質(zhì) ： 線性性質(zhì)：Zx1(t)+x

58、2(t)=x1(z)+x2(z)， Zmx(t)=mx(z) (m 為實數(shù)) (2) 位移定理：Zx(k-n)=z-nx(z),() 初值定理,(4) 終值定理,(5) 非一一對應(yīng)性： z變換與原函數(shù)不是一一對應(yīng)的關(guān)系。,脈沖傳遞函數(shù)的求?。?（1）對連續(xù)系統(tǒng)離散化（需要采樣開關(guān)、保持器） （2）求輸出輸入的z變換之比 由此可知：z 變換不但與原系統(tǒng)傳遞函數(shù)有關(guān)， 還與采樣開關(guān)、保持器性質(zhì)有關(guān)。 一般：G(z)G(s),典型傳遞函數(shù)的z變換 ：,對于w(s),若采用零階保持器:,若采用其它保持器，也可得到類似的公式。,分別對e、x1進(jìn)行離散相似，采用零階保持器：,x1(K+1)= e-T/4x1(K)+ (1-e-T/4)e(K),x2(K+1)= e-T/5x2(K)+ (1-e-T/5)x1(K) 仿真模型：e(K)=u(K)-x2(K) x1(K+1)= e-T/4x1(K)+ (1-e-T/4)e(K) x2(K+1)= e-T/5x2(K)+ (1-e-T/5)x1(K) y(K+1)= x2(K+1) 與前述（-16頁）離散相似法得到的模型是一樣的,四、 數(shù)字控制系統(tǒng)的數(shù)字仿真,1、系統(tǒng)構(gòu)成：數(shù)字控制器、實際對象,數(shù)字PID調(diào)節(jié):,例如：取Kp=1， Ki=0.3， Kd=1,數(shù)字控制器與常規(guī)PID對比,采樣