自動控制原理的MATLAB仿真與實(shí)踐第3章-線性系統(tǒng)的時域分析課件

11第3章主要內(nèi)容：

3.1時域分析的拉普拉斯變換法3.1.1連續(xù)時間函數(shù)的拉普拉斯變換3.1.2時域函數(shù)的拉氏反變換法3.1.3時域函數(shù)的部分分時展開法3.2時域分析的函數(shù)命令方法3.3作業(yè)與實(shí)驗(yàn)2第3章主要內(nèi)容：23.1時域分析的拉普拉斯變換法3

利用拉普拉斯變換法對線性系統(tǒng)進(jìn)行時域分析的一般步驟：（1）對系統(tǒng)的傳遞函數(shù)模型進(jìn)行部分分式展開，將其變?yōu)楹唵蝹鬟f函數(shù)之和；（2）利用拉普拉斯反變換，得到系統(tǒng)的輸出時間響應(yīng)函數(shù)；（3）繪制系統(tǒng)的響應(yīng)曲線;（4）通過改變系統(tǒng)的參數(shù)，觀察系統(tǒng)輸出響應(yīng)的變化情況，對系統(tǒng)的時域特性（如性能指標(biāo)等）進(jìn)行分析。本節(jié)介紹上述過程的MATLAB實(shí)現(xiàn)方法。3.1時域分析的拉普拉斯變換法3利用拉普拉斯變2.1.1生成模型的常用函數(shù)命令格式4連續(xù)時間信號f(t)的單邊拉普拉斯變換F(s)的定義：記為：拉普拉斯反變換的定義：記為：2.1.1生成模型的常用函數(shù)命令格式4連續(xù)時間信號f(t)s=tf(‘s’)：生成以s為變量的遞函數(shù)s。此時，s既是傳遞函數(shù)也是指定變量。

sys=tf(num,den)：生成傳遞函數(shù)模型，num,den分別為模型的分子和分母多項(xiàng)式系數(shù)向量。sys=tf(num,den,‘td’,v)：生成延遲時間td=v的傳遞函數(shù)模型。sys=zpk(z,p,k)：生成零極點(diǎn)增益模型，z,p,k分別為零點(diǎn)、極點(diǎn)和增益向量。[num,den]=rmodel(n,P)：隨機(jī)生成一個n階連續(xù)的傳遞函數(shù)模型，該系統(tǒng)具有P個輸出。[num,den]=ord2(wn,z)：生成固有頻率為wn，阻尼系數(shù)為z的連續(xù)二階系統(tǒng)模型系統(tǒng)。[num,dent]=pade(L,n)：返回延遲環(huán)節(jié)G(s)=e-Ls近似為n階多項(xiàng)式傳遞函數(shù)的num和den。5s=tf(‘s’)：生成以s為變量的遞函數(shù)s。此時，s既是傳2.1.2有理多項(xiàng)式分式傳遞函數(shù)模型的建立6

控制理論中，系統(tǒng)傳遞函數(shù)的有理多項(xiàng)式分式模型通常表示為：

在MATLAB中，多項(xiàng)式是用其降冪系數(shù)的排列構(gòu)成的數(shù)組來表示的，稱為多項(xiàng)式系數(shù)向量。

系統(tǒng)傳遞函數(shù)的分子和分母的多項(xiàng)式系數(shù)向量分別表示為：num=[bm,bm-1,…,b1,b0]den=[an,an-1,…,a1,a0]2.1.2有理多項(xiàng)式分式傳遞函數(shù)模型的建立67

在MATLAB語言中有專門對信號進(jìn)行正、反拉氏變換的函數(shù)命令。常用拉氏變換的函數(shù)命令格式及其說明見表3-1。F=laplace(f)：f(t)的拉氏變換，結(jié)果為F(s)，默認(rèn)變量為sF=laplace(f,v)：f(t)拉氏變換，并用v代替s，結(jié)果為F(v)F=laplace(f,v,u)：f(v)拉氏變換，并用u代替v，結(jié)果為F(u)f=ilaplace(F)

：F(s)的拉氏反變換，結(jié)果為f(t)，變量為tf=ilaplace(F,u)：F(s)拉氏反變換，用u代替s，結(jié)果為f(u)f=ilaplace(F,v,u)：

F(v)拉氏反變換，用u代v，結(jié)果為f(u)說明：默認(rèn)情況下，MTLAB給出的是單邊拉氏變換；使用函數(shù)laplace()及ilaplace()之前，需要用syms命令對所有要用到的變量進(jìn)行符號變量的定義。7在MATLAB語言中有專門【例3-1】試求函數(shù)

的拉氏變換式，并用拉氏反變換觀察變換結(jié)果。解：MATLAB程序如下：

>>clear;%清除所有變量symstAwbs

%定義符號變量t,A,w,b,sft=A*sin(w*t+b);%定義f(t)的符號函數(shù)ft的表達(dá)式Fs=laplace(ft)

%求ft的拉氏變換式Fs,即F(s)運(yùn)行結(jié)果：Fs=A*(cos(b)*w/(s^2+w^2)+sin(b)*s/(s^2+w^2))8【例3-1】試求函數(shù)當(dāng)取參數(shù)A=1、b=0、w=1時，帶入Fs>>A=1;b=0;w=1;%設(shè)定參數(shù)值Fs=A*(cos(b)*w/(s^2+w^2)+sin(b)*s/(s^2+w^2))

%計算給定參數(shù)下的F(s)運(yùn)行結(jié)果：Fs=1/(s^2+1)即得到

F(s)=L[f(t)]=L[sin(t)]=1/(s2+1)9當(dāng)取參數(shù)A=1、b=0、w=1時，帶入Fs9

可利用拉氏反變換對上述結(jié)果進(jìn)行檢驗(yàn)：>>ft=ilaplace(Fs)%求Fs的拉氏反變換式ft運(yùn)行結(jié)果：ft=sin(t)即

f(t)=L-1[F(s)]=L-1[1/(s2+1)]=sin(t)說明:MATLAB中的拉氏變換的函數(shù)命令laplace()，在默認(rèn)情況下是指拉氏右變換，其運(yùn)行結(jié)果是單邊函數(shù)。如L[f(t)]=laplace(sin(t))，實(shí)際上是指對ft=sin(t)u(t)的雙邊拉氏變換，

其中u(t)為單位階躍函數(shù)。10可利用拉氏反變換對上述結(jié)果進(jìn)行檢驗(yàn)：10

利用MATLAB拉氏反變換可以求出傳遞函數(shù)的時間函數(shù)，或系統(tǒng)的時域輸出響應(yīng)?！纠?-3】試求傳遞函數(shù)

的時間函數(shù)g(t)。解：求解G(s)的拉氏反變換程序如下：113.1.2時域函數(shù)的拉氏反變換法利用MATLAB拉氏反變換可以求出傳遞函數(shù)1>>clear;symss;Gs=(s^3+4*s^2+8*s+5)/(s^2+3*s+4);gt=ilaplace(Gs)運(yùn)行結(jié)果：gt=Dirac(t)+Dirac(1,t)-1/7*exp(-3/2*t)*7^(1/2)*sin(1/2*7^(1/2)*t)+exp(-3/2*t)*cos(1/2*7^(1/2)*t)其中，即得：12>>clear;symss;12【例3-4】某質(zhì)量-彈簧-阻尼(m-k-c)系統(tǒng)如圖3-1所示。其中，f(t)為外加力，x(t)為m的位移。試求：1）系統(tǒng)的MATLAB傳遞函數(shù)；2）輸入為單位階躍（f(t)=u(t)）時的時間響應(yīng)函數(shù)x(t)；3）繪制m=3;k=2;c=1時的時域響應(yīng)曲線。13【例3-4】某質(zhì)量-彈簧-阻尼(m-k-c)系統(tǒng)如圖3-1所解：1）由牛頓第二定律可得到該系統(tǒng)模型的數(shù)學(xué)表示：上式的拉氏變換后的傳遞函數(shù)模型為：生成上述模型的MATLAB模型：>>symsmkcsGs=1/(m*s^2+c*s+k)%生成G(s)的符號傳遞函數(shù)Gs運(yùn)行結(jié)果：Gs=1/(m*s^2+c*s+k)14解：1）由牛頓第二定律可得到該系統(tǒng)模型的數(shù)學(xué)表示：14152）f(t)=u(t)的拉氏變換為F(s)=1/s。所以輸出傳函為求上式的MATLAB拉氏反變換：>>Xs=Gs*(1/s)%建立輸出傳函X(s)xt=ilaplace(Xs)%求解輸出X(s)的時域響應(yīng)

152）f(t)=u(t)的拉氏變換為F(s)=1/s。所以16運(yùn)行結(jié)果：xt=-4*exp(-1/2*c/m*t)/(4*k*m-c^2)*m*cos(1/2*((4*k*m-c^2)/m^2)^(1/2)*t)+1/k*exp(-1/2*c/m*t)/(4*k*m-c^2)*c^2*cos(1/2*((4*k*m-c^2)/m^2)^(1/2)*t)-1/k*exp(-1/2*c/m*t)/(4*k*m-c^2)*m*((4*k*m-c^2)/m^2)^(1/2)*c*sin(1/2*((4*k*m-c^2)/m^2)^(1/2)*t)+1/k3）繪制m=3;k=2;c=1時的響應(yīng)曲線。程序如下：>>symstm=3;k=2;c=1;%參數(shù)賦值

16運(yùn)行結(jié)果：xt=-4*exp(-1/2*c/m*t)/(4*k*m-c^2)*m*cos(1/2*((4*k*m-c^2)/m^2)^(1/2)*t)+1/k*exp(-1/2*c/m*t)/(4*k*m-c^2)*c^2*cos(1/2*((4*k*m-c^2)/m^2)^(1/2)*t)-1/k*exp(-1/2*c/m*t)/(4*k*m-c^2)*m*((4*k*m-c^2)/m^2)^(1/2)*c*sin(1/2*((4*k*m-c^2)/m^2)^(1/2)*t)+1/k運(yùn)行結(jié)果：xt=-1/2*exp(-1/6*t)*cos(1/6*23^(1/2)*t)-1/46*exp(-1/6*t)*23^(1/2)*sin(1/6*23^(1/2)*t)+1/217xt=-4*exp(-1/2*c/m*t)/(4*k*m-c>>t=[0:0.1:10]*pi;%取繪圖區(qū)間采樣點(diǎn)xt=-1/2*exp(-1/6*t).*cos(1/6*23^(1/2)*t)-...1/46*exp(-6*t).*23^(1/2).*sin(1/6*23^(1/2)*t)+1/2;%計算所有采樣點(diǎn)上的函數(shù)值plot(t,xt);grid%繪圖，加?xùn)鸥窠Y(jié)果如圖3-2所示。18>>t=[0:0.1:10]*pi;1919

對于復(fù)雜的分式多項(xiàng)式傳函模型，首先，采用部分分式法將其分解成如下部分分式展開的形式：其中，pi為極點(diǎn)，

ri為有理函數(shù)Y(s)的留數(shù)，即k為余項(xiàng)；其次，對其求拉氏反變換；最后，由拉氏變換簡表得到時間響應(yīng)函數(shù)。203.1.3時域函數(shù)的部分分時展開法對于復(fù)雜的分式多項(xiàng)式傳函模型，首先，采用部20321通常有下列幾中情況及其相應(yīng)的處理方法：1）單極點(diǎn)時，即pi為不相同極點(diǎn)。則Y(s)的拉普拉斯逆變換為：2）有m重極點(diǎn)時。假設(shè)有m個重極點(diǎn)

，Y(s)的拉普拉斯逆變換為：3）有共軛極點(diǎn)時，逆變換較為復(fù)雜，需具體分析處理（參考教材），

21通常有下列幾中情況及其相應(yīng)的處理方法：，22

在MATLAB中，對于上述三種情況均可利用如下函數(shù)命令格式[r,p,k]=residue(num,den)對分式多項(xiàng)式進(jìn)行部分分式展開，然后再根據(jù)極點(diǎn)的不同情況進(jìn)行相應(yīng)的處理。其中，num,den為Y(s)的分子和分母多項(xiàng)式系數(shù)向量，[r,p,k]為返回的留數(shù)，極點(diǎn)，

和余項(xiàng)。22在MATLAB中，對于上述三種情況均可利用如下【例3-5】試求輸出傳遞函數(shù)的時間響應(yīng)。解：MATLAB程序如下：>>clear;num=[1-2];den=[13310];%系數(shù)向量[r,p,k]=residue(num,den)%求部分分式參數(shù)運(yùn)行結(jié)果：r=2.00002.00003.0000-2.000023【例3-5】試求輸出傳遞函數(shù)23p=-1.0000-1.0000-1.00000k=[]可見，系統(tǒng)有3個重根p1=p2=p3=-1和p4=0，則由式（3-3）和式（3-4）得到由式（3-4）知，m=3,對應(yīng)的逆變換為y(t)=L-1[Y(s)]?？梢灾苯永美献儞Q表得出結(jié)果。也可以通過MATLAB求解，程序如下：24p=24>>symssft1=ilaplace(3/(s+1)^3),ft2=ilaplace(2/(s+1)^2),ft3=ilaplace(2/(s+1)),ft4=ilaplace(-2/s)%求各分式逆變換得結(jié)果：ft1=3/2*t^2*exp(-t)ft2=2*t*exp(-t)ft3=2*exp(-t)ft4=-2>>ft=ft1+ft2+ft3+ft4%求f(t)的結(jié)果：ft=3/2*t^2*exp(-t)+2*t*exp(-t)+2*exp(-t)-2即25>>symss25【例3-7】設(shè)一個RLC網(wǎng)絡(luò)結(jié)構(gòu)如圖所示。其中ui為輸入信號，uc為輸出信號，網(wǎng)絡(luò)各變量的初始狀態(tài)均為零。試求：1）建立RLC網(wǎng)絡(luò)的傳遞函數(shù)

；2）當(dāng)ui=Eu(t)，用拉氏逆變換函數(shù)命令求uc；3）繪制不同參數(shù)時的單位階躍響應(yīng)uc(t)。解：1）建立RLC網(wǎng)絡(luò)的傳遞函數(shù)

(1)建立時域模型。根據(jù)R、L、C電壓定律及回路電壓方程

，可得RLC網(wǎng)絡(luò)的微分方程模型：其中，

26【例3-7】設(shè)一個RLC網(wǎng)絡(luò)結(jié)構(gòu)如圖所示。26（2）建立傳函模型

對上述微分模型的兩邊取拉氏變換，經(jīng)整理得到RLC網(wǎng)絡(luò)拉氏變換的標(biāo)準(zhǔn)二階系統(tǒng)傳遞函數(shù)：2）當(dāng)ui=Eu(t)時，由于L[ui]=L[E]=E/s，所以輸出的拉氏變換為：利用MATLAB拉氏逆變換函數(shù)命令求uc。程序如下：27（2）建立傳函模型27>>clear

;symstsEwnetwd%定義變量t,s,E,ωn,ζUc=E*wn^2/((s^2+2*et*wn*s+wn^2)*s);%建立函數(shù)uc=ilaplace(Uc,s,t)%拉氏逆變換求輸出uc(t)simple(uc)%簡化結(jié)果結(jié)果：ans=E*wn^2*(1/wn^2+1/(4*et^2*wn^2-4*wn^2)^(1/2)*(1/(-et*wn+1/2*(4*et^2*wn^2-4*wn^2)^(1/2))*exp((-et*wn+1/2*(4*et^2*wn^2-4*wn^2)^(1/2))*t)-1/(-et*wn-1/2*(4*et^2*wn^2-4*wn^2)^(1/2))*exp((-et*wn-1/2*(4*et^2*wn^2-4*wn^2)^(1/2))*t)))28>>clear

;283）繪制不同參數(shù)時的單位階躍響應(yīng)uc(t)。>>clear

;E=1;wn=3;t=0:0.02:6;%設(shè)置參數(shù)E和wn,時間采樣點(diǎn)tet=[0.1,0.2,0.5,0.95];%設(shè)置4個不同的阻尼向量etfori=1:4;%對et循環(huán)計算uc=E*wn^2*(1/wn^2+1/(4*et(i)^2*wn^2-4*wn^2)^(1/2)*(1/(-et(i)*wn+1/2*(4*et(i)^2*wn^2-4*wn^2)^(1/2))*exp((-et(i)*wn+1/2*(4*et(i)^2*wn^2-4*wn^2)^(1/2)).*t)-1/(-et(i)*wn-1/2*(4*et(i)^2*wn^2-4*wn^2)^(1/2)).*exp((-et(i)*wn-1/2*(4*et(i)^2*wn^2-4*wn^2)^(1/2)).*t)));plot(t,uc,'k-');holdon;%繪制響應(yīng)曲線并保持圖形end;293）繪制不同參數(shù)時的單位階躍響應(yīng)uc(t)。29title('uc(t)的二階欠阻尼響應(yīng)曲線')gtext('et=0.1')%在圖中對應(yīng)的響應(yīng)曲線的合適位置%點(diǎn)擊鼠標(biāo)左鍵添加阻尼參數(shù)標(biāo)示'et=0.1'gtext('et=0.2')gtext('et=0.5')gtext('et=0.95')xlabel('t');ylabel('uc(t)');%添加x和y軸標(biāo)示't'和'uc(t)'grid;%添加網(wǎng)格運(yùn)行得到響應(yīng)曲線見圖3-4。30title('uc(t)的二階欠阻尼響應(yīng)曲線')3031313.2時域分析的函數(shù)命令方法

上述介紹的MATLAB拉氏變換法，不但可以用來對通常情況下的微分方程、分式多項(xiàng)式、傳遞函數(shù)及系統(tǒng)結(jié)構(gòu)等函數(shù)進(jìn)行數(shù)學(xué)上的處理與分析，還可用于控制系統(tǒng)的時域分析。MATLAB還提供了特殊的函數(shù)命令，專用于對系統(tǒng)的分析。本節(jié)介紹常用MATLAB線性系統(tǒng)時域分析的函數(shù)命令格式及其使用方法。表3-2列出了常用時域響應(yīng)函數(shù)命令格式及說明。323.2時域分析的函數(shù)命令方法上述介紹的MA【例3-8】已知系統(tǒng)傳遞函數(shù)為

試求其單位脈沖響應(yīng)。解：使用專用函數(shù)命令來實(shí)現(xiàn)，程序如下：>>clear;n=[3];d=[1,1,3];%傳遞函數(shù)的系數(shù)向量Gs=tf(n,d);%建立模型G(s)impulse(Gs)%G(s)的單位脈沖響應(yīng)運(yùn)行結(jié)果如圖3-4所示。33【例3-8】已知系統(tǒng)傳遞函數(shù)為333434【例3-9】已知系統(tǒng)傳遞函數(shù)為求其單位階躍響應(yīng)。解：程序如下：>>clear;n=[1,7];d=[1,1.5,8,7];Gs=tf(n,d);step(Gs)運(yùn)行結(jié)果如圖3-5所示。35【例3-9】已知系統(tǒng)傳遞函數(shù)為3536【例3-14】已知某二階系統(tǒng)的超調(diào)量Mp=16.3%。試根據(jù)下列條件確定二階系統(tǒng)的開環(huán)傳遞函數(shù)G(s)：1）上升時間tr=0.1s；2）峰值時間tp=0.114s；3）調(diào)節(jié)時間ts=2s。

解：由二階系統(tǒng)模型結(jié)構(gòu)可知，其傳遞函數(shù)由兩個參數(shù)決定：阻尼比ζ和自然振蕩頻率ωn?？筛鶕?jù)二階系統(tǒng)的性能指標(biāo)中的任意兩個進(jìn)行求解即可確定模型參數(shù)。所以可以利用MATLAB的方程求解函數(shù)命令實(shí)現(xiàn)。具體格式：

X=solve('eq1','eq2',...,'eqn',val1,val2,...,valn)其中，eqi和vali，i=1,2,…,n,分別為方程和待求變量。36【例3-14】已知某二階系統(tǒng)的超調(diào)量Mp=16.3%。371）求解同時滿足超調(diào)量Mp=16.3%和上升時間tr=0.1s時的阻尼比ζ和自然振蕩頻率ωn。即求解下列方程組中的兩個變量ζ和ωn。設(shè)ζ=x，ωn=y，則上述方程改寫為371）求解同時滿足超調(diào)量Mp=16.3%和上升時間tr=038計算程序如下：>>clear;symsxyMptrx=solve('Mp-exp(-x*pi/sqrt(1-x^2))=0',x);%求x=ζx=subs(x,Mp,0.163)y=solve('tr-(pi-atan(sqrt(1-x^2)/x))/(y*sqrt(1-x^2))=0',y);%求y=ωny=subs(y,tr,0.1);y=subs(y,x)得到的結(jié)果為：x=0.5000；y=24.1852。即得二階系統(tǒng)模型為：38計算程序如下：392）求解同時滿足超調(diào)量Mp=16.3%和峰值時間tp=0.114s時的阻尼比ζ和自然振蕩頻率ωn。即求解下列方程組計算程序如下：>>clear;symsxyMptpx=solve('Mp-exp(-x*pi/sqrt(1-x^2))=0',x);%求阻尼比x=subs(x,Mp,0.163)y=solve('tp-pi/(y*sqrt(1-x^2))=0',y);%求y=ωny=subs(y,tp,0.114);y=subs(y,x)所得結(jié)果：x=ζ=0.5000；y=ωn=31.8219。將參數(shù)ζ和ωn代入G(s)即可。392）求解同時滿足超調(diào)量Mp=16.3%和峰值時間tp=040【例2-12】已知傳遞函已知傳遞函數(shù)G1和G2分別為試以G1作為前向通道，以G2作為反饋通道時的負(fù)反饋等效閉環(huán)模型G。解：實(shí)現(xiàn)反饋的等效傳函程序如下：方法1：利用多項(xiàng)式系數(shù)求解。>>clear;num1=[1,11];den1=[1,3,21];num2=[2,1];den2=[3,5,