電子科大-工程控制基礎(chǔ)上機(jī)(高端版)

基于Matlab的控制系統(tǒng)

的計(jì)算機(jī)輔助設(shè)計(jì)、分析及仿真1、Matlab簡介2、控制系統(tǒng)的數(shù)學(xué)描述與建模3、連續(xù)系統(tǒng)模型的離散化4、控制系統(tǒng)分析方法（穩(wěn)定性、時(shí)域和頻域分析）5、控制系統(tǒng)的設(shè)計(jì)與仿真6、Simulink仿真1MATLAB簡介MATLAB語言的簡潔高效性MATLAB語言的科學(xué)運(yùn)算功能MATLAB語言的繪圖功能MATLAB龐大的工具箱與模塊集MATLAB強(qiáng)大的動(dòng)態(tài)系統(tǒng)仿真功能

MATLAB語言是當(dāng)前國際上控制領(lǐng)域的首選計(jì)算機(jī)語言，也是很多理工科專業(yè)最適合的計(jì)算機(jī)數(shù)學(xué)語言。通過學(xué)習(xí)可更深入理解和掌握數(shù)學(xué)問題的求解思想，提高求解數(shù)學(xué)問題的能力，為今后其他專業(yè)課程的學(xué)習(xí)提供幫助。直接賦值語句賦值變量＝賦值表達(dá)式

例：>>a=pi^2a=9.8696

例：表示矩陣>>B=[1+9i,2+8i,3+7j;4+6j5+5i,6+4i;7+3i,8+2j1i]B=1.0000+9.0000i2.0000+8.0000i3.0000+7.0000i4.0000+6.0000i5.0000+5.0000i6.0000+4.0000i7.0000+3.0000i8.0000+2.0000i0+1.0000iMATLAB的基本語句函數(shù)調(diào)用語句[返回變量列表］＝函數(shù)名（輸入變量列表）例：［a,b,c]=my_fun(d,e,f,c)冒號表達(dá)式

v=s1:s2:s3

該函數(shù)生成一個(gè)行向量v，其中s1是起始值，s2是步長（若省略步長為1），s3是最大值。

例：用不同的步距生成(0,p)間向量。>>v1=0:0.2:piv1=Columns1through900.20000.40000.60000.80001.00001.20001.40001.6000Columns10through161.80002.00002.20002.40002.60002.80003.0000循環(huán)結(jié)構(gòu)for結(jié)構(gòu)while結(jié)構(gòu)MATLAB語言流程控制二維圖形繪制二維圖形繪制基本語句構(gòu)造向量:2、控制系統(tǒng)的數(shù)學(xué)描述與建模

控制系統(tǒng)的數(shù)學(xué)建模函數(shù)命令tf

()，用法：sys=tf(num,den)orsysd=tf(num,den,T)

——建立控制系統(tǒng)的傳遞函數(shù)模型，或者將零極點(diǎn)模型或者狀態(tài)空間模型轉(zhuǎn)換為傳遞函數(shù)模型；函數(shù)命令zpk()，用法：sys=zpk(z,p,k)orsysd=zpk(z,p,k,T)

——建立控制系統(tǒng)的零極點(diǎn)模型;函數(shù)命令ss()，用法：sys=ss(A,B,C,D)orsysd=ss(A,B,C,D,T)

——建立控制系統(tǒng)的狀態(tài)空間模型;函數(shù)命令tf2ss(),tf2zp(),zp2ss(),zp2tf(),ss2tf(),ss2zp()——數(shù)學(xué)模型的轉(zhuǎn)化;基于Simulink的動(dòng)態(tài)結(jié)構(gòu)模型.傳遞函數(shù)

單輸入單輸出線性連續(xù)系統(tǒng)的傳遞函數(shù)為其中m≤n。G(s)的分子多項(xiàng)式的根稱為系統(tǒng)的零點(diǎn),分母多項(xiàng)式的根稱為系統(tǒng)的極點(diǎn)。令分母多項(xiàng)式等于零,得系統(tǒng)的特征方程:D(s)=a0sn+a1sn－1+……+an－1s+an=0

因傳遞函數(shù)為多項(xiàng)式之比,所以我們先研究MATLAB是如何處理多項(xiàng)式的。MATLAB中多項(xiàng)式用行向量表示,行向量元素依次為降冪排列的多項(xiàng)式各項(xiàng)的系數(shù),例如多項(xiàng)式P(s)=s3+2s+4,其輸入為

>>P=［1024］

注意盡管s2項(xiàng)系數(shù)為0,但輸入P(s)時(shí)不可缺省0。

MATLAB下多項(xiàng)式乘法處理函數(shù)調(diào)用格式為C=conv(A,B)

例如給定兩個(gè)多項(xiàng)式A(s)=s+3和B(s)=10s2+20s+3,求C(s)=A(s)B(s),則應(yīng)先構(gòu)造多項(xiàng)式A(s)和B(s),然后再調(diào)用conv()函數(shù)來求C(s)>>A=［1,3］;B=［10,20,3］；>>C=conv(A,B)

C=1050639即得出的C(s)多項(xiàng)式為10s3

+50s2+63s+9

MATLAB提供的conv()函數(shù)的調(diào)用允許多級嵌套,例如

G(s)=4(s+2)(s+3)(s+4)可由下列的語句來輸入

>>G=4*conv(［1,2］,conv(［1,3］,［1,4］))對于其它復(fù)雜的表達(dá)式,如可由下列語句來輸入

>>num=conv(［1,1］,conv(［1,2,6］,［1,2,6］));>>den=conv(［1,0,0］,conv(［1,3］,［1,2,3,4］));>>G=tf(num,den)Transferfunction:

控制系統(tǒng)的典型連接串聯(lián)環(huán)節(jié)，用法：

sys=series(sys1,sys2)orsys=sys1*sys2*sys3*…*sysn

并聯(lián)環(huán)節(jié)，用法：

sys=parallel(sys1,sys2)orsys=sys1+sys2+sys3+…+sysn

反饋環(huán)節(jié)，用法：等效傳遞函數(shù)——G=feedback(G1,G2,sign)解編程如下：

num1=[0.1,1];den1=[0.4,0];sys1=tf(num1,den1);

num2=[15];den2=[0.054,1];sys2=tf(num2,den2);

num3=[1.5];den3=[0.12,1];sys3=tf(num3,den3);

sys123=sys1*sys2*sys3

運(yùn)行結(jié)果所以系統(tǒng)等效傳遞函數(shù)為

Transferfunction:2.25s+22.5---------------------------------0.002592s^3+0.0696s^2+0.4s已知某系統(tǒng)前向通道三個(gè)模塊的傳遞函數(shù)分別為。試求串聯(lián)連接的等效傳遞函數(shù)解編程如下

num1=[3];den1=[1,1];sys1=tf(num1,den1);

num2=[6,10];den2=[1,2,1];sys2=tf(num2,den2);

sys=sys1+sys2;num=sys.num{1}

den=sys.den{1}

運(yùn)行結(jié)果所以系統(tǒng)的等效傳遞函數(shù)為

num=092213den=1331已知兩子系統(tǒng)傳遞函數(shù)分別為：試求兩系統(tǒng)并聯(lián)連接的等效傳遞函數(shù)的num與den向量。3、連續(xù)系統(tǒng)模型的離散化

函數(shù)命令c2d()，用法:sysd=c2d(sys,T)orsysd=c2d(sys,T,method)

輸入?yún)⒘縮ys為連續(xù)時(shí)間模型對象,T,采樣周期;

Sysd,帶采樣時(shí)間t的離散事件模型;Method,用來指定離散化采用的方法:

′zoh′——采用零階保持器法;′foh′——采用一階保持器法;′tustin′——采用雙線形變換法;′prewarp′——采用頻率預(yù)曲折的雙線性變化

法;′matched′——采用零極點(diǎn)匹配法.解編程如下：

k=10;z=[-5];p=[-1,-3,-8];sys=zpk(z,p,k);t=0.1;sys1=tf(sys);

disp('Discretesystem--usingc2dwithzoh')s1=c2d(sys1,t,'zoh')

disp('Discretesystem--usingc2dwithmatched')s2=c2d(sys1,t,'matched')運(yùn)行結(jié)果Discretesystem--usingc2dwithzoh

Transferfunction:0.0401z^2+0.007485z-0.01929---------------------------------z^3-2.095z^2+1.41z-0.3012

Samplingtime:0.1Discretesystem--usingc2dwithmatched

Transferfunction:0.03596z^2+0.01415z-0.02181---------------------------------z^3-2.095z^2+1.41z-0.3012

Samplingtime:0.1已知連續(xù)系統(tǒng)零極點(diǎn)增益模型為試采用零階保持器與零極點(diǎn)匹配法求其離散傳遞函數(shù).設(shè)采樣周期T=0.1s.4控制系統(tǒng)的分析方法

控制系統(tǒng)的穩(wěn)定性分析特征方程的根，即閉環(huán)極點(diǎn)應(yīng)具有負(fù)實(shí)部或分布在左半s平面上?！苯优袛嗬щy。勞斯（Routh）穩(wěn)定性判據(jù)：由特征方程的系數(shù)來判斷。根軌跡法頻率響應(yīng)特性 根據(jù)S平面和Z平面之間的關(guān)系，離散系統(tǒng)的穩(wěn)定性可以由特征方程穩(wěn)定區(qū)域的根在Z平面中的位置來確定——必須位于Z平面中單位園的內(nèi)部。如果有一個(gè)根恰好位于單位園上則系統(tǒng)處于臨界穩(wěn)定，臨界穩(wěn)定在實(shí)踐中屬于不穩(wěn)定。穩(wěn)定區(qū)域不穩(wěn)定區(qū)域臨界穩(wěn)定Z平面的穩(wěn)定性條件運(yùn)行結(jié)果：已知系統(tǒng)的閉環(huán)傳遞函數(shù)為：

試判斷該閉環(huán)系統(tǒng)的穩(wěn)定性。解：根據(jù)題意，運(yùn)行下列MATLAB程序：num=[5410.630.5];den=[100000];[z,p]=tf2zp(num,den)ii=find(real(p)>0);n1=length(ii);if(n1>0),disp('SystemisUnstable');else,disp('SystemisStable');endz=-0.7822+0.5660i-0.7822-0.5660i0.4681+0.6367i0.4681-0.6367i-0.1718p=00000SystemisStable解：根據(jù)題意，運(yùn)行下列MATLAB程序：num=[1];den=[31-1-21];[z,p]=tf2zp(num,den)ii=find(abs(p)>1);n1=length(ii);if(n1>0),disp('SystemisUnstable');else,disp('SystemisStable');endz=

Emptymatrix:0-by-1

p=

-0.7357+0.6859i-0.7357-0.6859i0.5690+0.0753i0.5690-0.0753iSystemisUnstable運(yùn)行結(jié)果：已知系統(tǒng)的特征方程為如下所示,試判斷該系統(tǒng)穩(wěn)定性:控制系統(tǒng)的時(shí)域分析方法直接在時(shí)域內(nèi)對系統(tǒng)進(jìn)行分析；從輸出量的時(shí)域表達(dá)式分析系統(tǒng)的穩(wěn)定性、動(dòng)態(tài)性能和穩(wěn)態(tài)誤差；常用命令函數(shù)，step(),impulse(),initial(),dstep,dimpulse,dinitial,nec.任意函數(shù)作用下系統(tǒng)的響應(yīng)

用線性仿真函數(shù)lsim來實(shí)現(xiàn)，其調(diào)用格式為

［y，x］=lsim(num，den，u，t)

式中；y(t)為系統(tǒng)輸出響應(yīng)；x(t)為系統(tǒng)狀態(tài)響應(yīng)；u為系統(tǒng)輸入信號；t為仿真時(shí)間。

求的單位階躍響應(yīng)。解：編制Matlab程序如下：num=[25];den=[1,4,25];g=tf(num,den);step(g)或num=[25];den=[1,4,25];step(num,den)求當(dāng)時(shí)的單位沖擊響應(yīng)。解：

編制Matlab程序如下：wn=5;zeta=0.6;num=wn.^2;den=[1,2*zeta*wn,wn.^2];impulse(num,den)設(shè)線性離散系統(tǒng)如圖所示，且a=1，K=1，T=1s，輸入為單位階躍序列。試分析系統(tǒng)的過渡過程。解：被控對象Z傳遞函數(shù)編制Matlab程序如下：num=[0.3680.264];den=[1-10.632];t=0:1:25;dstep(num,den,t)或[y,x]=dstep(num,den,t)程序運(yùn)行結(jié)果：y=00.36801.00001.39941.39941.14700.89460.80170.86830.99371.07691.08091.03230.98120.96080.97270.99751.01471.01631.00700.99670.99230.99430.99921.00281.0033x=[]調(diào)整時(shí)間ts=12s（12個(gè)采樣周期），超調(diào)量p=40％，峰值時(shí)間tp=3s

，振蕩次數(shù)N=1.5次，衰減比=5:1

，穩(wěn)態(tài)誤差ess=0。例：已知一單位負(fù)反饋系統(tǒng)，其開環(huán)零極點(diǎn)增益模型為試分析在單位階躍及斜坡信號作用下系統(tǒng)的穩(wěn)態(tài)誤差。利用Matlab進(jìn)行穩(wěn)態(tài)誤差分析單位階躍信號t=0.1;k=6;z=-0.5;p=[-210];[n1,d1]=zp2tf(z,p,k);s=tf(n1,d1);sys=feedback(s,1);roots(sys.den{1})ans=-0.1084+1.9541i-0.1084-1.9541i-0.7832系統(tǒng)穩(wěn)定（a）系統(tǒng)穩(wěn)定性判斷sysd=c2d(sys,t)step(sys);t1=[0:t:300]';y=step(sysd,t1);subplot(121),plot(t1,y),grid;subplot(122),ess=1-y;plot(t1,ess),gridess(length(ess))階躍響應(yīng)及穩(wěn)態(tài)誤差ans=-8.5931e-014系統(tǒng)為一型系統(tǒng)，其Kp=，即ess=0。t=0.1;k=6;z=-0.5;p=[-210];[n1,d1]=zp2tf(z,p,k);s=tf(n1,d1);sys=feedback(s,1);num=sys.num{1};den=[sys.den{1},0];sys=tf(num,den);sysd=c2d(sys,t)step(sys);t1=[0:t:50]';y=step(sysd,t1);subplot(121),plot(t1,[t1,y]),grid;subplot(122),ess=t1-y;plot(t1,ess),gridess(length(ess))斜坡響應(yīng)及穩(wěn)態(tài)誤差系統(tǒng)為一型系統(tǒng)，其Kv=C，即ess=C。ans=-0.6678控制系統(tǒng)的頻域分析方法以系統(tǒng)的頻率特性來研究系統(tǒng)的性能；對數(shù)頻率特性圖（Bode圖），

[mag,phase,w]=dbode(sys,w)極坐標(biāo)圖（Nyquist圖）

[re,im,w]=dnyquist(sys,w)求系統(tǒng)增益裕度和相位裕度的函數(shù)，margin();[Gm,Pm,Wcp,wcg]=margin(sys)根軌跡圖，[p,z]=pzmap(sys),[p,z]=rlocus(sys).對數(shù)頻率特性Matlab提供的相關(guān)函數(shù)為：函數(shù)dbode()繪制線性離散系統(tǒng)對數(shù)頻率特性圖。函數(shù)dbode()的使用格式如下：[mag,phase,w]=dbode(num,den,Ts,w)當(dāng)包含左方變量時(shí)，屏幕上不顯示伯德圖，返回角頻率點(diǎn)w處的幅值和相位角（單位為度）。函數(shù)[Gm,Pm,ωcg,ωcp]=margin(sys)除繪制線性離散系統(tǒng)對數(shù)頻率特性圖外，還可以求出系統(tǒng)的幅值裕度Gm、相位裕度Pm以及所在點(diǎn)對應(yīng)的頻率值ωcg,ωcp試求所示離散系統(tǒng)，在T=1s，K=1時(shí)的開環(huán)對數(shù)頻率特性（伯德圖）。在T=1s，K=1時(shí)，編制Matlab程序如下：num=[0.368,0.264];den=[1,-1.368,0.386];sys=tf(num,den,1)sysw=d2c(sys,'tustin')syszpk=zpk(sysw)margin(syszpk)title('離散系統(tǒng)伯德圖');grid解：可知，系統(tǒng)的開環(huán)傳遞函數(shù)為直接用dbode()

函數(shù)：num=[0.368,0.264];den=[1,-1.368,0.368];dbode(num,den,1)title('離散系統(tǒng)伯德圖');grid繪制程序2系統(tǒng)的相位裕度幅值裕度系統(tǒng)性能預(yù)測：穩(wěn)定，有一定超調(diào)。結(jié)論：極坐標(biāo)頻率特性Matlab提供的相關(guān)函數(shù)為：

Dnyquist(num,den,t)試求上例所示離散系統(tǒng)，在T=1s，K=1時(shí)的極坐標(biāo)頻率特性曲線。解：編制程序如下：num=[0.368,0.264];den=[1,-1.368,0.386];dnyquist(num,den,1)title('離散系統(tǒng)極坐標(biāo)圖');已知系統(tǒng)傳遞函數(shù)為繪制零極點(diǎn)圖和根軌跡圖.MATLAB命令為:num=[0,0,2,8,12,8,2];den=[1,5,10,10,5,1,0];sys=tf(num,den);subplot(1,2,1)pzmap(sys)[p,z]=pzmap(sys)subplot(1,2,2)rlocus(sys)運(yùn)行結(jié)果：p=0-1.0010+0.0007i-1.0010-0.0007i-0.9996+0.0012i-0.9996-0.0012i-0.9988

z=-1.0002-1.0000+0.0002i-1.0000-0.0002i-0.9998

解：系統(tǒng)的開環(huán)Ｚ傳遞函數(shù)試求所示離散系統(tǒng)，在T=1s和2s時(shí)的根軌跡圖，并求出相應(yīng)的臨界增益。例：編制Matlab程序如下：t1=1;t2=2;k=1;n=2;fori=1:nifi<=1t=t1;elseifi<=2t=t2;endb1=t-1+exp(-t);b2=1-exp(-t)-t*exp(-t);a0=1;a1=1+exp(-t);a2=exp(-t);num=[k*b1,k*b2];den=[a0,-a1,a2];figure(i)

rlocus(num,den)title('根軌跡圖')holdp=0:0.07:2*pi;x=sin(p);y=cos(p);plot(x,y,'.')holdv=[-4,2,-2,2];axis(v)grid[k,poles]=rlocfind(num,den)endT=1秒T=2秒零界增益K=2.52零界增益K=0.58結(jié)論采樣周期T與系統(tǒng)的臨界增益K密切相關(guān)。T增長，K降低。系統(tǒng)的零極點(diǎn)改變時(shí)，根軌跡的形狀也發(fā)生變化，系統(tǒng)的性能也跟著變化。為了滿足性能指標(biāo)的要求，可以采用零極點(diǎn)對消的辦法，或加入新的零極點(diǎn)，來改善系統(tǒng)的性能。5控制系統(tǒng)的設(shè)計(jì)和仿真頻域設(shè)計(jì)方法，根據(jù)系統(tǒng)的傳遞函數(shù)進(jìn)行設(shè)計(jì)。時(shí)域設(shè)計(jì)方法，基于系統(tǒng)的狀態(tài)空間模型進(jìn)行設(shè)計(jì)。頻域設(shè)計(jì)方法超前校正的設(shè)計(jì)步驟如下：

(1)根據(jù)穩(wěn)態(tài)性能要求，確定系統(tǒng)開環(huán)增益K。

(2)根據(jù)外環(huán)增益K，畫出校正前系統(tǒng)的Bode圖，并計(jì)算未校正系統(tǒng)的相位裕度等性能指標(biāo)．以檢測是否滿足要求.(3)確定需要增加的最大相位超前角：式中,γ0為期望相位裕度,γ1為原來開環(huán)系統(tǒng)相位裕度.(4)由確定α值及最大相位超前角所對應(yīng)的頻率ωm,并取新的幅值穿越頻率ωcg=

ωm。

(5)確定超前校正器傳遞函數(shù)，并繪制校正后的Bode圖和進(jìn)行驗(yàn)證。經(jīng)典控制理論設(shè)計(jì)已知單位正反饋系統(tǒng)被控對象的傳遞函數(shù)為試用Bode圖設(shè)計(jì)方法對系統(tǒng)進(jìn)行超前校正設(shè)計(jì)。性能指標(biāo)要求如下：(1)在斜坡信號r(t)=v0t,ess≤0.0001v0;(2)在系統(tǒng)校正后,相位裕量48o<γ<53o.解:(1)求K0速度誤差系數(shù)Kv=K=K0,(K是系統(tǒng)的開環(huán)增益)

Ess=v0/Kv=v0/K0≤0.0001v0所以,Kv=K0≥1000,取K0=1000(2)作