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

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

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

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

5、進行辨識、修正模型結(jié)構(gòu)和數(shù)據(jù),3. 火電廠熱工過程建模與仿真,火電廠熱工過程特點： 系統(tǒng)復(fù)雜、多變量、大遲延、非線性、時變性 建立模型時要有針對性： (1) 要明確建模的目的及模型的用途； (2) 確定建模的方法 (3) 對系統(tǒng)進行合理的簡化 (4) 確定模型的輸入、輸出 (5) 建立模型 (6) 分析、驗證、修正模型 (7) 模型應(yīng)用,火電廠設(shè)備,計算機接口,實際DCS操作,數(shù)學(xué)模型,計算機接口、網(wǎng)絡(luò),仿真DCS操作,實際電廠,仿真機,機理法建模實例：彈性阻尼系統(tǒng)建模實例 問題：研究行使在不平路面上的汽車的顛簸情況。 方法：物理模型研究 過程：（1）采用類比方法得到相應(yīng)的概念模型 （2）建立

6、數(shù)學(xué)模型 建立彈性阻尼系統(tǒng)模型過程： （i）分析 （ii）假定 K、B線性 （iii）建模：牛頓第二定律 得到彈性阻尼系統(tǒng)模型：,四、數(shù)字仿真的基本概念 數(shù)字仿真：求解數(shù)學(xué)模型數(shù)值解的過程 仿真分類：按模型類型分；按系統(tǒng)的連續(xù)性分；按時間尺度分； 實時仿真、非實時仿真,彈性阻尼系統(tǒng)建模實例,（3）模型求解,（i）模型變換 轉(zhuǎn)換為狀態(tài)方程：,初始條件：x1(0)=0, x2(0)=0;,（ii）選擇數(shù)值計算方法得到仿真模型 采用歐拉法可寫出迭代公式如下：,iii）選擇編程語言，編程并進行仿真研究 例如：Matlab、C、Basic語言等 （4）對仿真結(jié)果進行分析并驗證模型 一般需要試驗數(shù)據(jù)的驗證

7、，或采用實驗數(shù)據(jù)修正模型 （5）模型應(yīng)用，在確認模型正確后，可應(yīng)用該模型對彈性阻尼系統(tǒng)進行研究,物理原型概念模型（分析、類比得到） 數(shù)學(xué)模型（假設(shè)，采用物理定律得到） 仿真模型（采用一定的數(shù)值計算方法） 數(shù)字模型（程序代碼、通過編程得到） 實用模型（多次的模型分析、檢驗、修 正） 模型應(yīng)用（通過模型研究物理系統(tǒng)的特性） 建模與仿真需要明確的問題： 明確建模目的 模型的輸入、輸出 模型的描述形式 建模對象機理及特性分析 編程語言的選擇 模型的驗證方法：靜態(tài)工況、動態(tài)過程,動態(tài)過程建模與仿真研究的一般步驟,仿真的應(yīng)用領(lǐng)域，發(fā)展史，火電站仿真 機的發(fā)展及現(xiàn)狀 航空航天 核電、火電 交通、運輸 通訊、

8、電子 經(jīng)濟、人類 我?；痣姺抡婕夹g(shù)的發(fā)展 1111111111111,五、 計算機輔助分析 、數(shù)學(xué)模型之間的相互轉(zhuǎn)化 、控制系統(tǒng)穩(wěn)定性分析 頻域法：由W(S)分析系統(tǒng)的穩(wěn)定性， 勞斯判據(jù)、波特圖、奈奎斯特圖、 根軌跡等 時域法：狀態(tài)空間法。,六、 計算機輔助設(shè)計 頻域法： 基于傳遞函數(shù)的控制器設(shè)計方法； 例如：串聯(lián)校正器的概念及設(shè)計方法； 極點配置設(shè)計方法； 時域法： 基于狀態(tài)空間模型的控制器設(shè)計方法； 線性二次型最優(yōu)調(diào)節(jié)器的設(shè)計方法； 觀測器的基本設(shè)計方法； 過程控制系統(tǒng)：PID控制器設(shè)計及參數(shù)整定方法；,機理法建模： 實例4、椅子著地問題： 4條腿長度相等的椅子，放在起伏不平的地面上，4條

9、腿是否能同時著地？,假定地面為連續(xù)平面，則轉(zhuǎn)動椅 子在90內(nèi)，總能找到一點使四 條腿同時著地。 讓椅子原地旋轉(zhuǎn)，x為旋轉(zhuǎn)角度。 設(shè)f(x)為A、C兩腿距地面距離 之和，g(x)為B、D兩腿距地面距 離之和。,由于任意三腿總在一個平面上，所以任意地點三腿可同時著地?？傆衒(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

10、(x)=0,f(x)0 可得：g(90)=f(0)0, f(90)=g(0)=0 即：h(90)= f(90)-g(90)狀態(tài)方程 方法（1）：傳遞函數(shù) 微分方程狀態(tài)方程,例3：求下列傳遞函數(shù)對象的狀態(tài)方程描述,解：先將傳遞函數(shù)進行拉氏反變換得到微分方程,再將微分方程轉(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)方程 方法（3）：直接轉(zhuǎn)化 并聯(lián)法,等效框圖：,并聯(lián)法的特點：方法通用，A為對角陣,同理：,狀態(tài)方程：,輸出方程：,初值：x1

11、(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.05x3+x1,把e,e1代入上面各式，整理可得狀態(tài)方程描述:,狀態(tài)方程：,輸出方程： y=x3,初值：x1(0)

12、=0; x2(0)=0; x3(0)=0;,本章復(fù)習(xí)題,2. 將微分方程轉(zhuǎn)為狀態(tài)方程描述，假定初值為0,3. 將傳遞函數(shù)轉(zhuǎn)為狀態(tài)方程描述，要求A為對角陣。,簡述熱工過程機理建模的基本假設(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的值 即：給定一個時間計算步距DT，求y(DT), y(2DT), y(3DT), y(4DT),.y(kDT).y(nDT),例：彈性阻尼系

13、統(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ù)值積分仿真算法原理,對于： X=AX+Bu 設(shè) F(t)=A X + Bu，則：X=F(t),兩邊積分,得X計算通式：,取t=kDT,

14、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)會造成較大的誤差。,歐拉公式幾何解釋：,精確解：,歐拉公式： 用矩形面積近似代替積分面積,DT*F(K),誤差,分析：歐拉公式簡單實用，但當DT較大時，會造成較大的誤差。要想提高精度，必須選擇較小的計算

15、那步距DT。DT越小，精度越高。,DT的提高受計算機精度、速度的限制，在DT一定的情況下，如何提高計算精度？,問題：,（2）梯形公式：,用梯形面積近似代替積分面積,?。?由此，得到梯形公式仿真算法：,梯形公式存在的問題：,計算F(k+1)時需要用到X(k+1),無法使用！,怎么辦？,放棄,改進,如何改進？,把歐拉公式引入梯形公式，先采用歐拉公式預(yù)測F(k+1),再采用梯形公式計算。,改進方法：,改進后的梯形公式：,預(yù)報校正公式,改進的梯形公式幾何解釋,在DT一定的情況下，能否進一步提高計算精度？,問題：,歐拉公式：只采用k點的信息，計算F*,算法分析：,梯形公式：采用了k點和k+1點的信息，計

16、算F*,是否可以采用k點，k+1,k+1/2點的信息計算F*,進一步提高精度？,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ù)測k+1/2點的F值F(k+1/2), 記做F3,5、計算F*:,4. 根據(jù)F3預(yù)測k+1點的F值F(k+1),

17、記做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,反之，若|1+DTa|1 則數(shù)值計算收斂,顯然，數(shù)值計算發(fā)散,精度與穩(wěn)定性,需要注意： 1.精度與穩(wěn)定性是兩個不同的

18、概念,兩者之間沒有 直接的聯(lián)系； 2.上述每種數(shù)值積分算法都有自己的收斂區(qū)間； 3.收斂區(qū)間與算法有關(guān)，也與對象特性有關(guān); 4.仿真步距（計算步距）DT，越小，算法穩(wěn)定性越好；,5. 其它仿真算法,（1）三階龍格-庫塔法,以上數(shù)值計算方法為單步法，由 K 時刻的值推導(dǎo)K+1 時刻的值。,（2） 阿達姆斯法,多步法：,由 K，K-1，k-2.時刻的值推導(dǎo)K+1 時刻的值。,誤差 :O(T4),單步法:,阿達姆斯法：,（3）隱式歐拉公式,若取 F*=F(k+1) , 得：,對于高階系統(tǒng)需要求逆矩陣，采用高級語言較為復(fù)雜，但對于單個微分方程：,把F(k+1) 代入得：,隱式歐拉公式如何求解？,隱式歐

19、拉公式的優(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,F2, ST=20,DT=0.1; main() NP = (int)(ST / DT);,C語言仿真程序

20、設(shè)計,for(i=0;ihelp 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),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); tit

21、le(彈性阻尼系統(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)定狀態(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

22、, 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=K1*u1; W2=K2*sqrt(L); Y(i)=L; L=L+DT*(W1-W2)/(RO*F); end plot(1:NP)*DT,Y); title(單容水箱仿真);,水箱水位仿真研究,水箱水位仿真曲線,DT仿真研究-顯式歐拉公式,

23、水箱水位仿真研究-隱式歐拉公式,D. 仿真模型：隱式歐拉公式,非線性方程，無法直接求得L(k+1)的表達式,需要把上式中的非線性部分進行線性化后，才能得到L(k+1)的表達式。,常用處理方法：,即：把原來非線性項分為兩部分：系數(shù)部分和線性部分。系數(shù)部分用K時刻的值近似計算，得到B2(k),整理得：,for i=1:NP B2=K2/sqrt(L); L=( 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)

24、過程應(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 （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)過程,不同仿真步距是輸出曲線：顯式歐拉,隱式歐拉公式仿真模型如下：,按照同樣的

25、線性化方法可得：,整理得仿真模型：,不同仿真步距是輸出曲線：隱式歐拉,如圖所示：流體網(wǎng)絡(luò)中一個壓力節(jié)點,建模時假定： (1) P點周圍的管道容積都集中在P處 (2) 忽略流體密度變化,根據(jù)質(zhì)量守衡可得：,二、 質(zhì)量守恒應(yīng)用之二：流體網(wǎng)絡(luò)壓力節(jié)點模型,對于熱工過程而言，一般情況下，溫度的變化速度要比壓力的變化速度慢得多，忽略溫度的變化,K1的大小反映了流體的可壓縮性，K1值越大，說明流體可壓性越強，反之，可壓性則越小，對于不可壓流體K10。 需要注意的是，只有理想流體才可能是不可壓流體，實際流體都是可壓的。,忽略流體密度變化，則各路流量與壓差關(guān)系如下：,若采用顯式歐拉公式：,若采用隱式歐拉公式：

26、,按前述方法對上式進行線性化，設(shè)：,壓力節(jié)點仿真模型：,考慮流體網(wǎng)絡(luò)中各節(jié)點的計算順序，在計算P(k+1)時，P1(k+1)、P2(k+1)、P3(k+1)可能還是未知，為使模型通用，所有相臨節(jié)點壓力都取為K時刻的值，即：,如圖， 對于小流量，為簡單起見，直接引用流量,若采用隱式歐拉公式：,對于任意多支路的壓力節(jié)點，假定周圍有i個壓力節(jié)點，j路微小流量：,建模時假設(shè)： （1）外壁絕熱； （2）采用集總參數(shù)法，假定金屬管道壁溫均勻一致； （3）忽略工質(zhì)定壓比熱Cp隨溫度的變化,三、 能量守恒應(yīng)用之一：絕熱蒸汽管道模型,根據(jù)質(zhì)量守衡：,根據(jù)能量守衡,對于工質(zhì)有：,對于金屬有：,工質(zhì)和金屬之間的換熱

27、量：,對流換熱系數(shù)：,1. 基于溫度的集總參數(shù)模型,由于實際過程中壓力、流量變化較快，而溫度、焓值變化較慢，建立模型時一般把流動過程與傳熱過程分開考慮。質(zhì)量守衡動態(tài)過程在壓力節(jié)點模型或水箱水位模型中考慮，能量守衡方程在溫度、焓值模型中考慮。,取蒸汽側(cè)集總參數(shù)T=Tout ，采用顯式歐拉公式得：,若采用隱式的歐拉公式，注意到Q與T、Tm有關(guān)，為了增加模型的穩(wěn)定性，需要將Q帶入,考慮到計算順序，前面公式中用到的Tm(k+1)要用Tm(k)代替才能計算,由于水蒸汽的定壓比熱隨壓力、溫度的變化很大，定值假定會造成較大誤差，在建模時常采用焓值作為熱平衡計算的基礎(chǔ)，建立基于焓值的集總參數(shù)模型:,H為集總參

28、數(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，模型中蒸汽側(cè)不變，只需在金屬側(cè)減去管道對環(huán)境散熱量 ：,計算Tm的遞推公式也應(yīng)進行相應(yīng)的改動。,建模時假定：,（1）整個換熱器用一段逆流管代替，沿管長方向吸熱均勻，金屬溫度均勻一致，為Tm ； （2）忽略煙氣定壓比熱隨溫度的變化； （3）忽略流動的動

29、態(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, Ms 輸出：hs2, Ts2, Tg2, Tm,整理后模型的一般形式,整理好的仿真模型,整理好的仿真模型如下,混合聯(lián)箱：聯(lián)箱內(nèi)有一電加熱器，熱功率為Q,建模時假定： （1）忽略流體密度變化，水箱等截面，面積為F； （2）集總參數(shù)法，整個聯(lián)箱內(nèi)工

30、質(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，采用隱式歐拉公式：,仍按前面的線性化處理方法可得隱式歐拉公式模型:,模型輸入: W1，t1, W2, t2, Q,模型輸出: L，t3,已知：,第四章、離散相似法,一、離散相似法原理,狀態(tài)方程：,輸出方程：,初值：,求：Y的數(shù)值解,方法一：前述的數(shù)值

31、積分法,方法二：離散相似法,通解：,K時刻：,K+1：,如何計算eAT及積分是該算法的關(guān)鍵！,遞推關(guān)系：,與數(shù)值積分法比較,數(shù)值積分法：近似處理F,離散相似法：近似處理,取t1=(K+1)T-t,進行積分變換：,零階保持器離散相似算法如下：,在區(qū)間kT,(k+1)T, 近似取U(t)=U(kT),零階保持器,（1）、計算(T) ：,一般指數(shù) ：,矩陣指數(shù) ：,近似取前幾項,（2）、計算m(T) ：,A、當系統(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,近似取,代入得：,整理得：,其中：,

32、三角階保持器,三角保持器離散相似算法如下：,方法三、通過拉氏變化求解：,求拉氏變換： 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ù)卷積定理：,對于：,設(shè)：,由此可得：,與離散相似法比較：,可得：,即：可以通過拉氏反變換計算,同理：x2=-1/5x2+1/5x1 e=u-x2 ,整理得：,解: (1) 先求系統(tǒng)的狀態(tài)方程 設(shè)定狀態(tài)變量x1,x2,例1：對如下系統(tǒng)進行仿真,x1=-1/4x1+

33、1/4e,求(T),m(T), n(T),由于不是對角陣，求逆矩陣比較困難!,若A是對角陣，例如,可見，離散相似法過于繁瑣!,解決方法1：,系統(tǒng)為一階系統(tǒng)時，A=a, eAT=eaT,對每個一階系統(tǒng)分別進行離散相似,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)+ (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/5x

34、2(K)+ (1-e-T/5)x1(K) y(K+1)= x2(K+1),仿真模型如下：,（1）原來是對u進行近似處理，這次分別對兩個子系統(tǒng)的輸入e,x1進行了近似處理，誤差增大。 （2）當u是階躍輸入時，對u進行近似處理不會引入誤差，而第二種方法對e,x1近似處理，引入了兩次誤差。,注意：,分析發(fā)現(xiàn)，計算x1(k+1)用到e(K+1) e(K+1)=u(K+1)-x2(K+1) 可見，該模型是無法使用的。,為了減小誤差，采用三角保持器：,即：,A在區(qū)間kT,(k+1)T,整理得：,處理方法,方法1：對系統(tǒng)1（與反饋環(huán)節(jié)相鄰）使用零階保持器,方法2：改進三角保持器,一階保持器,B在區(qū)間kT,(

35、k+1)T,取,這樣近似處理后，算法實用化。,滯后一拍三角保持器,整理得：,二、保持器的概念,b. 一階保持器,c. 三角保持器,d. 滯后一拍的三角保持器,采樣系統(tǒng)：采樣開關(guān)、保持器 采樣處理離散，信號恢復(fù)（保持）相似,2保持器種類及頻率特性,1. 離散相似概念的來源：,a. 零階保持器,3. 各種離散相似處理后信號的變化,（1）.多次離散后造成信號有較大的相位滯后,（2）.多次離散后造成信號有一定的幅值衰減,為了提高精度，對離散相似信號補償，包括對幅值及相位的補償,為了簡單,一般取,4. 補償器,補償器：,即:輸入信號超前一拍,例如對于零階保持器,加入補償環(huán)節(jié)后:,小結(jié)：,1. 對u離散相

36、似，各種算法都能使用，但需要求得整個系統(tǒng)的狀態(tài)方程； 2對于負反饋系統(tǒng)，對e離散相似，三角保持器不能使用，對x1離散相似，各種算法都能使用； 3. 對e離散相似, 不能使用補償器； 4. 零階保持器可以無失真的再現(xiàn)階躍信號； 5. 其它保持器可以無失真的再現(xiàn)斜波信號； 6. 對各個子系統(tǒng)分別進行離散相似可以簡化仿真算法，但會引入新的誤差；,離散相似法程序設(shè)計,采用matlab 函數(shù)進行彈性阻尼系統(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

37、=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);,對于單容水箱系統(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;

38、 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-FA); Y(i)=L; L=FA*L+FAM*u1; end plot(1:NP)*DT,Y,r); title(采用離散相似法進行單容水箱仿真);,仿真結(jié)果分析： 離散相似法與 隱式歐拉公式對比,三、Z變換法,原理：脈沖傳遞函數(shù)直接對應(yīng)差

39、分方程,1. 脈沖傳遞函數(shù)的求?。狠斎胼敵龅膠變換之比,拉普拉斯變換定義:,：拉普拉斯算子，,z變換定義：,已知X(t)的采樣信號x*(t),對于kQBQC, 即QA最差，QB為除QA之外的最差點。 則新點A1: x10*=x11+x12-x10; x20*=x21+x22-x20; if (QA-QC)/QAEPSQ, 結(jié)束。,(3) 單純形的壓縮： A、若QA1QB, 則新三角形中A1最差，若按前面方法求反射點，則反射點為原來的z1。此時，說明新點沿反射方向前進太多，應(yīng)適當減少前進距離。 即取A2: x10*-x1c=k(x1c-x10) (-1k1 ) x10*=(x11+x12)*(k+1)/2-kx10; 同理：x20*=(x21+x22)*(k+1)/2-kx20;,4) 在(3)B 中取k=0.5若QA2QB,取k= - 0.5 求得A3點的QA3重新判斷。,(5) 若QA3QB,取k=0 求得A4點, 無法形成三角形 可構(gòu)造新點C1，x1c1=x12+delta1*lx1, x2c1=x22+delta1*lx2, 由C1,B,C三點構(gòu)成新三角形，轉(zhuǎn)入（2）。,編程改進： 1。為保證三角形的性狀，將dta,Ti尋優(yōu)區(qū)間等量化，例如：dta=0.1,10,ti=1,10