MATLAB控制工具箱的应用及线性系统的运动分析

信控学院上机实验信控学院上机实验第页共页第页共页实验报告课程 线性系统理论基础 实验日期 专业班姓名 学号 同组 实验名称MATLAB控制工具箱的应用及线性系统的运动分析评批阅教师签 一、实验目的1、学习掌握MATLAB控制工具箱中的基本命令的操作方法；2、掌握线性系统的运动分析方法。二、实验内容自选控制对象模型，应用以下命令，并写出结果。1）step,damp,pzmap,rlocus,rlocfind,bode,margin,2）tf2ss,ss2tf,tf2zp,zp2ss；3）ss2ss,jordan,canon,eig。掌握线性系统的运动分析方法 1)已知A0 1 2 3

，求e

At。(用三种方法求解) MATLAB2.8曲线，求解时域性能指标。(加图标题、坐标轴标注及图标)利用MATLAB求解书上例2.12（c1 0并画出状态应和输出响应曲线。(加图标题、坐标轴标注及图标)4)P361.4（2）1.5（3）；P562.3（3）三、实验环境1、计算机120台；2、MATLAB6.X软件1套。四、实验原理（或程序框图）及步骤1MATLAB设系统的模型如式(1-1)所示：xAxBu(1-1)

yCxDu

xRn uRm yRp其中A为n×n维系数矩阵；B为n×m维输入矩阵；C为p×n维输出矩阵；D为p×m维传递矩阵，一般情况下为0。系统的传递函数阵和状态空间表达式之间的关系如式(1-2)所示：(1-2)

G(s)

num((s)den(s)

C(sIA)1BD式(1-2)num(sm；den(s表示传递函数阵的分母多项式，按s降幂排列的后，各项系数用向量表示。五、程序源代码自选控制对象模型，应用以下命令，并写出结果。step,damp,pzmap,rlocus,rlocfind,bode,margin,Step（系统单位阶跃响应）程序：num=[1,4];den=[1,2,8];step(num,den)Damp（计算系统模型的固有频率）程序：num=[1];den=[1,10,10];sys=tf(num,den);[wzeta]=damp(sys)w=3.87303.8730zeta=0.64550.6455Pzmap:（绘制连续系统的零极点图）程序：H=tf([1-13],[21-1]);pzmap(H)Rlocus:（求系统根轨迹）程序：num=[0001];den=conv([10],[14+sqrt(-1)]);den=conv(den,[14-sqrt(-1)]);rlocus(num,den);v=[-41-22];axis(v);Rlocfind:确定闭环根位置对应增益值K程序：G=tf([1,5,6],[1,8,3,25]);rlocus(G);[k,r]=rlocfind(G)G_c=feedback(G,1);Bode:程序：num=[1,4];den=[1,2,8];bode(num,den)Margin:从频率响应数据中计算出幅值裕度、相角裕度以及对应的频率num=[4];den=[12margin(num,den),grid;w=0.1:1000;[mag,phase,w]=bode(num,den);Nyquist：绘制系统奈氏图n=[0011];m=[1200];sys=tf(n,m);nyquist(sys);tf2ss,ss2tf,tf2zp,tf2ss:传递函数阵转换为状态空间模型num=[0153];den=[1234];[A,B,C,D]=tf2ss(num,den)A=-2 -3 -41 0 00 1 0B=100C=1 5 3D=0Ss2tf：状态空间模型转换为状态空间A=[010;001;-4-3-2];B=[1;3;-6];C=[100];D=0;[num,den]=ss2tf(A,B,C,D,1)num=0 1.0000 5.0000 3.0000den=1.0000 2.0000 3.0000 4.0000Tf2zp：将系统的传递函数形式转换为零点，极点，增益形式num=[0153];den=[1234];[z1,p1,k1]=tf2zp(num,den)z1=-4.3028-0.6972p1=-1.6506+0.0000i-0.1747+1.5469i-0.1747-1.5469ik1=1Zp2ss：将系统零极点形式转换为传递函数形式z=[-2-3];p=[0;-1;-3;-5];k=3;[A,B,C,D]=zp2ss(z,p,k)A=-8.0000 -3.8730 0 03.8730 0 0 0-3.0000 -2.3238 -1.0000 00 0 1.0000 0B=1010C=0 0 0 3D=0ss2ss,jordan,canon,Ss2ss:线性非奇异变换A=[010;001;-4-3-2];B=[1;3;-6];C=[100];D=0;P=[111;-1-2-3;149];sys_in=ss(A,B,C,D);sys_out=ss2ss(sys_in,inv(P))sys_out=a=x1x2x3x1-217x22-4-14x3-114b=u1x1 x2 -9x3 2.5c=x1 x2 x3y1 1 1 d=u1y1 Continuous-timestate-spacemodel.Jordan：将广义特征向量计算矩阵对应的约旦向量A=[010;001;-4-3-2];J=jordan(A)J=-1.6506+0.0000i0.0000+0.0000i0.0000+0.0000i0.0000+0.0000i-0.1747-1.5469i0.0000+0.0000i0.0000+0.0000i0.0000+0.0000i-0.1747+1.5469iCanon:规范形转换函数num_1=[121];den_1=[156];sys_1=tf(num_1,den_1,-1)sys_can_ss=canon(sys_1,'model')sys_1=z^2+2z+1z^2+5z+6Sampletime:unspecifiedDiscrete-timetransferfunction.sys_can_ss=a=x1x2x1-30x20-2b=u1x1-7.211x2-5.657c=x1x2y10.5547-0.1768d=u1y11Sampletime:unspecifiedDiscrete-timestate-spacemodel%即为所求的模型的对角线规范形实现Eig：求矩阵的全部特征值A=[010;001;-4-3-2];D=eig(A)D=-1.6506+0.0000i-0.1747+1.5469i-0.1747-1.5469i掌握线性系统的运动分析方法 1)已知A0 1 2 3

，求e

At。(用三种方法求解) 法一：状态转移矩阵的指数矩阵计算法a=[01;-2-3];symst;eat1=expm(a*t)法二：拉氏反变换计算法a=[01;-2-3];symsst;G=inv(s*eye(size(a))-a)eat2=ilaplace(G)G=[(s+3)/(s^2+3*s+2), 1/(s^2+3*s+2)][ -2/(s^2+3*s+2), 法三：非奇异变换法a=[01;-2-3];symst;[P,D]=eig(a);Q=inv(P);eat3=P*expm(D*t)*QMATLAB2.8曲线，求解时域性能指标。(加图标题、坐标轴标注及图标)a=[-10;0-2;];b=[1;1]c=[1.50.5];d=0;G=ss(a,b,c,d);x0=[2;3]symsst;G0=inv(s*eye(size(a))-a);x1=ilaplace(G0)*x0G1=inv(s*eye(size(a))-a)*bx2=ilaplace(G1/s)x=x1+x2y=c*xforI=1:61tt=0.1*(I-1)xt(:,I)=subs(x(:),'t',tt)yt(I)=subs(y,'t',tt)endplot(0:60,[xt;yt])利用MATLAB求解书上例2.12（c1 0并画出状态应和输出响应曲线。(加图标题、坐标轴标注及图标)G=[01;-0.16-1]h=[1;1]x0=[1;-1]symsznkthta=inv(z*eye(size(G))-G)*zthtak=iztrans(thta,k)uz=z/(z-1)xk=iztrans(thta*x0+thta/z*h*uz)4)P361.4（2）1.5（3）；P562.3（3）1.4(2)A=[214;020;001]B=[10;34;21]C=[351]D=[00][num,den]=ss2tf(A,B,C,D,1)验证：num=[0 20.0000 55.0000 -4.0000]den=[1 2 -4 -8][A,B,C,D]=tf2ss(num,den)num=02055-4den=12-4-81．5（3）num=[1422;0311]den=[1232][A,B,C,D]=tf2ss(num,den)2．3a=[01;-6b=[1;0]c=[1-1]d=0G=ss(a,b,c,d)x0=[1;1]symsstG0=inv(s*eye(size(a))-a)x1=ilaplace(G0)*x0G1=inv(s*eye(size(a))-a)*bx2=ilaplace(G1/s)x=x1+x2y=c*x六、实验数据、结果分析（2）掌握线性系统的运动分析方法 1)已知A0 1 2 3

，求e

At。(用三种方法求解) eat1=[ -exp(-2*t)+2*exp(-t), exp(-t)-exp(-2*t)][-2*exp(-t)+2*exp(-2*t), eat2=[ -exp(-2*t)+2*exp(-t), exp(-t)-exp(-2*t)][-2*exp(-t)+2*exp(-2*t), eat3=[ -exp(-2*t)+2*exp(-t), exp(-t)-exp(-2*t)][-2*exp(-t)+2*exp(-2*t), 2*exp(-2*t)-exp(-t)]MATLAB2.8曲线，求解时域性能指标。(加图标题、坐标轴标注及图标)4.51.1系统的状态响应和输出响应43.532.52输出响应1.5状态响应*210.501020状态响应t(s)30405060利用MATLAB4.51.1系统的状态响应和输出响应43.532.52输出响应1.5状态响应*210.501020状态响应t(s)30405060G= 0 1.0000-0.1600 -1.0000h= 11x0= 1-1thta=[25*z*(z+1)/(25*z^2+25*z+4), 25*z/(25*z^2+25*z+4)][ -4*z/(25*z^2+25*z+4), thtak=[4/3*(-1/5)^k-1/3*(-4/5)^k,5/3*(-1/5)^k-5/3*(-4/5)^k][-4/15*(-1/5)^k+4/15*(-4/5)^k,-1/3*(-1/5)^k+4/3*(-4/5)^k]uz=z/(z-1)xk=[-17/6*(-1/5)^n+22/9*(-4/5)^n+25/18]图1.3图1.3系统的状态响应和输出相应32.52输出响应1.51状态响应*10.50-0.5-1状态相应*2-1.50102030t(s)4050604)P361.4（2）num=020.000055.0000-4.0000den=12-4-8验证：A=-2 4 81 0 00 1 0B= 100C=20 55 -4D=01.5（3）；num=14220311den=1232A= -2-3-21000B= 10100C= 2-103D= 1011P562.3（3）a=01-6-5b=10c=d=01-1a= x1 x1 0 1x2 -6 -5b= u1x1 1x2 0c= x1 y1 1 -1d= u1y1 0Continuous-timemodel.x0=11G0=[(s+5)/(s^2+5*s+6), 1/(s^2+5*s+6)][