现代控制理论实验指导书

1、现代控制理论实验指导书现代控制理论实验指导书 信息与电子工程学院2008年1月一、 MATLAB简介1980年美国学者MOLER博士在研究线性代数时，为了解决矩阵运算非常烦琐这一问题，构思并开发了一种用于矩阵运算的矩阵实验室（MATrix LABoratory）软件。由于MATLAB提供了功能非常强大的矩阵处理和绘图功能，吸引了很多控制界的名家在自己擅长的领域编写了一些具有特殊意义的MATLAB工具箱，从而空前扩大了MATLAB的功能，使其成为了国际上最流行的控制系统计算机辅助设计的软件工具。下面简单介绍自动控制理论实验涉及到的MATLAB知识（一）SIMULINK（模块图仿真）的使用1、 进

2、入WINDOWS：2、 进入MATLAB：点击WIN窗口中MATLAB图标。3、 进入SIMULINK：在MATLAB窗口中键入SIMULINK命令，即可弹出SIMULINK模块库。如图1-1所示：图1-1 SIMULINK模块库从左到右分别为信号源、输出、离散、线性、非线性、其他七个子模块库，每个子模块库分别包含若干个不同的模块。4、建立自己的实验文件：（1）FILENEW。（2）FILESAVE AS。（键入自己的文件名。注意：文件名不能用纯数字）5、从SIMULINK模块库中挑选所需的模块：（1）把SIMULIMK窗和自己的实验窗并列平铺。（2）在SIMULINK窗中打开SOURCES模

3、块库。（3）从SOURCES库(信号源模块库)中拖出STEP INPUT（单位阶跃输入） 模块至自己的实验窗口。（4）类似（2）（3）从SINKS库中拖出AUTO-SCALE GRAPH( 自动图形输出) 模块。（5）从LINEAR、UNLINEAR库中分别选出INTEGRATOR（积分）、GAIN （比例）、TRANFER-Fun（传递函数）等模块。（6）通过TRANFER-FORNT模块构造积分、惯性、振荡环节。 双击TRANFER-FORNT模块，弹出如图1-2对话框图1-2 Transfer Fun 模块对话框其中Numerator为分子多项式的系数 ， Denominator为分母多

4、项式的系数，通过修改多项式的系数实现不同环节的转变。举例如下：例1：要实现传递函数，则令Numerator为1，2；Denominator为2，3，4；例2：要将传递函数变为积分环节：，则令Numerator为1；Denominator为1，0；例3：要将传递函数变为振荡环节：，则令Numerator为1；Denominator为1，1，1；（在此传函中阻尼系数为0.5）例4：要将传递函数变为实际微分环节：，则令Numerator为1，0；Denominator为1，1；实际微分环节的传递函数为:分子分母同除以Td,则为因此,上式中分子s的系数即为Kd值,分母中常数项为Td的倒数（2）设仿真过

5、程参数：SIMULATION菜单PARAMETERS菜单项。如图1-3：图1-3 仿真时间的设置其中： Start Time为仿真开始时间，计算机默认的值为0.0。 Stop time为仿真终止时间，计算机默认的值为999999。 Mix Step Size为仿真最小步长。 Max Step Size为仿真最大步长。， Tolerance为仿真精度。仿真开始前应对Stop Time 进行修改，如改为10秒，50秒或200秒，再根据实际情况进行调整。注：在此使用的仿真算法为五阶的龙格库塔法（RUNGE-KUTTA5）（3）进行仿真：SIMULATION菜单中START菜单项。 注：要使曲线光滑可

6、以通过提高仿真精度或增加采样点来实现。（二）频域和时域分析法1 根轨迹图的绘制（1） 乘积多项式的系数可利用多项式相乘函数CONV（X，Y），该函数直接在 MATLAB COMMAND WINDOWS 窗口中使用。例：A(z)=(z-1)B(z)=(z2+1.15z-1)求A(z)*B(z)的多项式系数，可按如下步骤计算：X=1，1；Y=1，1.15，1；C=CONV（X，Y）上面式子等效于C=CONV(1，-1，1，1.15，1)注：无常数项时，常数项系数记为0（2） 绘制根轨迹图例：绘制如下传递函数的根轨迹图G(s)=可在MATLAB COMMAND WINDOW窗口中输入下列命令：NUM

7、=1，1；DEN=2，0，3；RLOCUS(NUM，DEN)（3） 绘制Z平面的根轨迹图在MATLAB COMMAND WINDOWS窗口中输入下列命令：ZGRID（NEW）NUM=A，B，C.；DEN=D，E，F.；RLOCUS（NUM，DEN）（4） 在根轨迹图上求系统的闭环主导极点在MATLAB COMMAND WINDOW窗口中输入下列语句：R，K=RLOCUS(NUM，DEN)即可得出不同的K值对应的不同根值R，从中找出尽量接近要求的主导极点及其对应的K值。在MATLAB COMMAND WINDOW窗口中输入下列语句： R=RLOCUS(NUM，DEN，K)在得出的K值基础上不断调

8、整K值, 得出满足Y/X=tgq的复根，则此时的K值即为开环放大系数。（已知闭环主导极点的阻尼比z为已知值，则以原点为起点作角度为q=COS-1z的直线，与根轨迹相交得S1点，S1与实轴相对称的点即为S2，若S1，2=XYj则tgq=Y/X，即求主导极点的问题转化为求一复根，使其满足Y/X=tg q的要求。）2伯德图的绘制（1）绘制伯德图 在MATLAB COMMAND WINDOW窗口中输入下列语句： NUM=A，B； DEN=C，D，E； BODE(NUM，DEN)（2）绘制离散控制系统BODE图 绘制离散控制系统BODE图时需将Z域的函数表达式变换为W域的，使用 BODE（NUM，DEN

9、）绘制BODE图。（3）求系统的增益裕量(Gm)、相位裕量(Pm)和穿越频率Wcp 使用如下语句： Gm，Pm，Wcg，Wcp=MARGIN(NUM，DEN)（4）求新的增益穿越频率Wc 设在该穿越频率Wc处，G0(jW)的相角为q，使用如下语句： MAR，PHA，W=BODE(NUM，DEN) 即可得出不同的W值对应相角(PHA)和增益(MAR)，找出近似满足PHA=q 时对应的W值，即为Wc。（三）状态空间分析法1、 矩阵的基本运算（1）求矩阵的秩设矩阵A=0，2；0，-2，则它的秩可直接在MATLAB WINDOWS COMMOND 窗口中调用函数R=RANK（A）得到。（2）矩阵的乘法

10、设有A，B两个矩阵，若求其相乘的结果可在MATLAB WINDOWS COMMOND窗口中输入如下命令：A=1，2；3，4B=3，5；6，8C=A*B则C即为矩阵A和B 相乘得到的新矩阵。（3）矩阵的转置设A矩阵为1 ，2；3，4 ，则它的转置矩阵为C=A（4）矩阵的逆 C=inv(A)2、由状态空间模型转换成系统的传递函数（阵）SS2TF在MATLAB WINDOWS COMMAND窗口中直接调用如下函数：NUM，DEN=SS2TF（A，B，C，D，IU）即可求出传递函数分子和分母多项式的系数，由此可直接写出传递函数的表达式。其中A，B，C，D矩阵表示系统的状态空间模型，而IU为输入的代号，

11、对于单输入系统来说，IU=1。注意：矩阵的输入是按行进行的，先输入第一行，再输入第二行，依次类推。行与行之间用分号隔开，每行中不同元素用空格或逗号隔开。3、由系统的传递函数（阵）转换成状态空间模型TF2SSA，B，C，D=TF2SS（NUM，DEN）为能控型（I型）实现（A为友矩阵）。4、实现状态空间模型的线性变换，转化为任意型SS2SSA1，B1，C1，D1=SS2SS（A，B，C，D，T）原系统 通过Z=TX（T为任意非奇异矩阵），即X=T-1Z 的线性变换，得 5、实现状态空间模型的线性变换，转化为对角型CANONA1，B1，C1，D1，T=CANON（A，B，C，D，mod）原系统 通

12、过Z=TX（T为非奇异矩阵），即X=T-1Z 的线性变换，得 为对角型 注意：仅适用于系统有互异特征根的情况，当有相同特征根时，运行结果有错。6、用ACKERMAN算法计算状态反馈矩阵K及状态观测器输出误差反馈矩阵G在MATLAB WINDOWS COMMOND 窗口中调用函数K=ACKER（A，B，P），其中A，B为状态方程系数矩阵，P为要求配置到的极点。例如，A=0，0，4；-3，-3，-2；0，1，2；B=2；1；3；若想引入状态反馈矩阵K，使得闭环系统的极点位置为-3，-1+j, -1-j ,则可以使用下面的MATLAB命令来完成A=0，0，4；-3，-3，-2；0，1，1；B=0；1

13、；3；P=-3；-1+j；-1-j；K=ACKER（A，B，P）即可求出状态反馈矩阵K。同理使用函数G=ACKER（A,C,P）可求出状态观测器输出误差反馈矩阵G。二、实验指导实验1 基于MATLAB的系统数学模型的转换1.0 实验设备 PC计算机1台（要求P4-1.8G以上），MATLAB6.X或MATLAB7.X软件1套。1.1实验目的学习多变量系统状态空间表达式的建立方法、了解系统状态空间表达式与传递函数相互转换的方法；通过编程、上机调试，掌握多变量系统状态空间表达式与传递函数相互转换方法。1.2实验原理说明设系统的模型如式(1.1) 所示 (1.1)其中A为维系数矩阵、B为维输入矩阵，

14、C为维输出矩阵，D为传递阵，一般情况下为0，只有n和m维数相同时， D为维。系统的传递函数阵和状态空间表达式之间的关系如式(1.2)示。 (1.2)式中，表示传递函数矩阵的分子的系数矩阵，其维数是；表示传递函数阵的最小共倍式的系数向量，并按s降幂排列的分母。若在SISO系统中，传递函数阵退化为传递函数，表示传递函数的分子多项式的系数向量，表示传递函数的分母多项式的系数向量。1.3实验步骤 根据所给系统的传递函数或（A、B、C阵），依据系统的传递函数阵和状态空间表达式之间的关系如式(1.2)，采用MATLAB的file.m编程。注意：ss2tf和tf2ss是互为逆转换的指令； 在MATLAB界面

15、下调试程序，并检查是否运行正确。例1.1 已知SISO系统的状态空间表达式为(1.3)，求系统的传递函数。 (1.3)程序:%首先给A、B、C阵赋值；A=0 1 0;0 0 1;-4 -3 -2;B=1;3;-6;C=1 0 0;D=0;%状态空间表达式转换成传递函数阵的格式为num,den=ss2tf(A,B,C,D,IU)num,den=ss2tf(A,B,C,D,1) 程序运行结果:num = 0 1.0000 5.0000 3.0000den = 1.0000 2.0000 3.0000 4.0000从程序运行结果得到:系统的传递函数为: (1.4)例1.2 从系统的传递函数(1.4)

16、式求状态空间表达式。程序:num =0 1 5 3; %在给num赋值时，在系数前补0，使num和den赋值的个数相同。den =1 2 3 4;A,B,C,D=tf2ss(num,den)程序运行结果:A = -2 -3 -4 1 0 0 0 1 0B = 1 0 0C = 1 5 3D =0由于一个系统的状态空间表达式并不唯一, 例1.2程序运行结果虽然不等于式(1.3)中的A、B、C阵，但该结果与式(1.3)是等效的。不妨对上述结果进行验证。例1.3 对上述结果进行验证编程%将例1.2上述结果赋值给A、B、C、D阵A =-2 -3 -4;1 0 0; 0 1 0;B =1;0;0;C =

17、1 5 3;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.0000程序运行结果与例1.1完全相同，说明从状态空间模型到传递函数之间的相互转换完全正确。例1.1所给的系统状态空间模型和例1.2的运行结果虽然不同，但两者是等效的。例1.4 已知SISO系统的状态空间表达式为(1.3)，用SS2SS，CANON函数进行线性变换。（同学自己编程）1.4 实验要求 在运行以上例程序的基础上，应用MATLAB对(1.5)系统仿照例1.2编程，求系统的A、B、C、D阵；然

18、后再仿照例1.3进行验证。并写出实验报告。 (1.5)提示：num =0 0 1 2；0 1 5 3；在MIMO系统的状态空间模型转换成传递函数阵时，D与C阵必须维数相同。实验2 利用MATLAB分析状态空间模型2.0 实验设备 PC计算机1台（要求P4-1.8G以上），MATLAB6.X或MATLAB7.X软件1套。2.1实验目的学习多变量系统状态转移矩阵和状态空间表达式的求解方法；通过编程、上机调试，掌握多变量系统状态空间表达式的计算机仿真分析。2.2实验原理说明设系统的模型如式(2.1) 所示 (2.1)初始状态为状态响应为：非齐次状态的解由两部分组成：第一部分是在初始状态作用下的自由运

19、动，第二部分是在系统输入作用下的强迫运动。2.3实验步骤1状态转移矩阵的求解例：已知：，求。源程序：% ex1.m% 状态转移矩阵的指数矩阵计算法A=0 1；-2 -3；syms t； % 定义变量eat1=expm(A*t) % 求源程序：% ex2.m% 状态转移矩阵的拉氏反变换计算法A=0 1；-2 -3；syms t s; G=inv(s*eye(size(A)-A) % 求的逆eat2=ilaplace(G) % 求的逆的拉氏反变换2.4 实验要求已知系统方程为， 当单位阶跃函数时，用MATLAB编程求系统的状态响应和输出响应。（提示：分别求出零输入响应和零状态响应，再求和。；）实验

20、3 系统能控性和能观性分析3.0 实验设备PC计算机1台（要求P4-1.8G以上），MATLAB6.X或MATLAB7.X软件1套。3.1 系统的能控性分析实验3.1.1 实验目的学习多变量系统状态能控性的定义；通过用MATLAB编程、上机调试，掌握多变量系统能控性判别方法。3.1.2 实验原理说明设系统的状态空间表达式 （3.1）系统的能控分析是多变量系统设计的基础，包括能控性的定义和能控性的判别。系统状态能控性的定义的核心是：对于线性连续定常系统（3.1）,若存在一个分段连续的输入函数U（t），在有限的时间t0，t1内，能把任一给定的初态x（t0）转移至预期的终端x（t1），则称此状态是能

21、控的。若系统所有的状态都是能控的，则称该系统是状态完全能控的。系统输出能控性是指输入函数U（t）加入到系统，在有限的时间t0，t1内，能把任一给定的初态x（t0）转移至预期的终态输出y（t1）。能控性判别分为状态能控性判别和输出能控性判别。状态能控性分为一般判别和直接判别法，后者是针对系统的系数阵A是对角标准形或约当标准形的系统，状态能控性判别时不用计算，应用公式直接判断，是一种直接简易法；前者状态能控性分为一般判别是应用最广泛的一种判别法。输出能控性判别式为： （3.2）状态能控性判别式为： （3.3）3.1.3 实验步骤 根据所给系统的统阵A和输入阵B，依据能控性判别式，采用MATLAB的

22、file.m编程； 在MATLAB界面下调试程序，并检查是否运行正确。例3.1 已知系统矩阵A和输入阵B分别如下，判断系统的状态能控性。， 程序：A =6.6667,-10.6667,-0.3333;1.0000,0,1;0,1.0000,2; B=0; 1; 1; q1=B; q2=A*B; %将AB的结果放在q2中 q3=A2*B; %将A2B的结果放在q3中 Qc=q1 q2 q3 %将能控矩阵Qc显示在MATLAB的窗口Q=rank(Qc) %能控矩阵Qc的秩放在Q程序运行结果:Qc = 0 -11.0000 -85.0003 1.0000 1.0000 -8.0000 1.0000

23、3.0000 7.0000Q = 3从程序运行结果可知，能控矩阵Qc的秩为3=n，所以系统是状态能控性的。3.1.4 实验要求在运行以上例程序的基础上，编程（M 文件）判别下面系统的能控性。 提示：从B阵看，输人维数m=2，Qc的维数为n(mn)=48,而Q=rank(Qc)语句要求Qc是方阵，所以先令，然后Q=rank(R)。对MATLAB新的版本，Qc可为任意阵。3.2 系统的能观性分析实验3.2.1 实验目的学习多变量系统状态能观性的定义；通过编程、上机调试，掌握多变量系统能观性判别方法。3.2.2 实验原理说明设系统的状态空间表达式如式（3.1）所示。系统的能观性分析是多变量系统设计的

24、基础，包括能观性的定义和能观性的判别。系统状态能观性的定义：对于线性连续定常系统（3.1）,如果对t0时刻存在ta，t0ta，根据t0，ta上的y(t)的测量值，能够唯一地确定S系统在t0时刻的任意初始状态x0，则称系统S在t0时刻是状态完全能观测的，或简称系统在t0，ta区间上能观测。状态能观性分为一般判别和直接判别法，后者是针对系统的系数阵A是对角标准形或约当标准形的系统，状态能观性判别时不用计算，应用公式直接判断，是一种直接简易法；前者状态能观性分为一般判别是应用最广泛的一种判别法。状态能控性判别式为： （3.4）3.2.3 实验步骤根据所给系统的系统阵A和输入阵B，采用MATLAB的f

25、ile.m编程； 在MATLAB界面下调试程序，并检查是否运行正确。例3.2 已知系统阵A和输入阵C分别如下，判断系统的状态能观性。， 程序：A =6.6667,-10.6667,-0.3333;1.0000,0,1;0,1.0000,2; C=1 0 2; q1=C; q2=C*A; %将CA的结果放在q2中 q3=C*A2; %将CA2的结果放在q3中 Qo=q1; q2; q3 %将能观矩阵Qo显示在MATLAB的窗口Q=rank(Qo) %能观矩阵Qo的秩放在Q程序运行结果:Qo = 1.0000 0 2.0000 6.6667 -8.6667 3.6667 35.7782 -67.4

26、450 -3.5553Q =3从程序运行结果可知，能控矩阵Qo的秩为3=n，由式（2.4）可知，系统是状态完全能观性的。3.2.4 实验要求在运行以上例程序的基础上，编程判别3.1.4所给系统的能观性。3.3 系统的结构分解按能控性结构分解: A1,B1,C1,T1,K1=ctrbf(A,B,C)按能观性结构分解: A2,B2,C2,T2,K2=obsvf(A,B,C)T1、 T2为变换阵； K1、K2为能控、能观状态变量个数。试编程，对下列系统进行能控、能观性结构分解， 实验4 多变量系统的稳定性分析4.0 实验设备 PC计算机1台（要求P4-1.8G以上），MATLAB6.X或MATLAB

27、7.X软件1套。4.1实验目的 学习多变量系统稳定性分析方法； 通过用MATLAB编程、上机调试，掌握多变量系统稳定性分析方法4.2实验原理说明设系统的模型如式(4.1) 所示 (4.1)只要系统的A的特征根实部为负，系统就是状态稳定的，可写成： (4.2)当状态方程是系统的最小实现时，系统的状态渐近稳定与系统的BIBO（有界输入有界输出）稳定等价；当时，若系统状态渐近稳定则系统一定是BIBO稳定的。4.3 实验步骤根据所给系统的A、B、和C阵，依据3.2实验原理，采用MATLAB的file.m编程；在MATLAB界面下调试程序，并检查是否运行正确。例4.1 已知系数阵A、B、和C阵分别如下,

28、分析系统的状态稳定性。 (4.3)编程:A=0 1 0;0 0 1;-4 -3 -2;B=1;3;-6;C=1 0 0;D=0;z,p,k=ss2zp(A,B,C,D,1)程序运行结果:z = -4.3028 -0.6972p = -1.6506 -0.1747 + 1.5469i -0.1747 - 1.5469ik = 1 由于系统的零、极点均具有负的实部，则系统是BIBO稳定的；又因为状态方程是系统的最小实现，系统的状态渐近稳定与系统的BIBO稳定等价，所以系统是状态渐近稳定的。4.4 实验要求分析系统式(4.4)的状态稳定性，用MATLAB编程，并写出实验报告。 (4.4)提示：编程时

29、，先从传递函数G(s)求出A、B、C阵，然后仿照例4.1编程。实验5 状态观测器的设计5.0 实验设备 PC计算机1台（要求P4-1.8G以上）、MATLAB6.X或MATLAB7.X软件1套。5.1 实验目的学习用MATLAB软件编程设计状态观测器的方法，并在SIMULINK仿真界面下搭建设计全阶状态观测器，用示波器观察的变化情况。学习掌握现代控制理论的实际工程应用。5.2 实验原理说明设系统的模型如式(5.1)所示 (5.1)系统状态观测器包括全阶观测器和降阶观测器。设计全阶状态观测器的条件是系统状态完全能观。全阶状态观测器的方程为： (5.2)5.3 实验步骤根据所给系统的A、B、和C阵

30、，依据（5.2）式，采用MATLAB的file.m编程；在MATLAB界面下调试程序，并检查是否运行正确。例5.1 已知系数阵A、B、和C阵分别如式(5.3)示,设计全阶状态观测器，要求状态观测器的极点为-3 3。， ， (5.3)首先验证系统是状态完全能观的。通过该系统的对偶系统极点配置求出状态反馈阵K，则状态观测器的增益阵为G=-KT。虽然所给系统是能控标准型，但其对偶系统不是能控标准型，设状态反馈阵，先按能控标准型进行极点配置求出，再还原到非能控标准型状态空间中，因此，Tc为将对偶系统化成能控标准型的变换阵。 编程:% 对该系统的对偶系统极点配置求出状态反馈阵K，则状态观测器的输出误差反

31、馈阵为G=-K。%设，Tc为将对偶系统化成能控标准型的变换阵。A=0 1;-2 -1; % 所给系统的对偶系统的状态矩阵A=AC=0;1; % 对偶系统的输出矩阵C=BB=1 0; %对偶系统的输入矩阵B=CD=0;num,den=ss2tf(A,B,C,D,1); %求出对偶系统特征多项式denf=1 6 9; %期望极点的特征多项式的系数阵k2=den(:,2)-denf(:,2); %计算k2=d1-a1，特征方程S的系数与期望极点的特征多项式S的系数相减k1=den(:,3)-denf(:,3); %计算k1=d2-a2，常数项相减 %显示阵Qc=B A*B; %计算对偶系统的能控阵p

32、1=0 1*inv(Qc);Tc=p1; p1*A %计算将对偶系统化成能控标准型的变换阵Tc并显示变换阵 %计算对偶系统状态反馈阵K并显示G=-K %计算状态观测器输出误差反馈阵G并显示或直接使用函数G=ACKER（A,C,P），可求出状态观测器输出误差反馈矩阵G。程序运行结果:G = 5 2 所以，状态观测器输出误差反馈阵G=5 2T。则状态观测器的方程为 在SIMULINK仿真界面下画出带有状态观测器的状态结构图，并进行仿真。观察e=y-z的变化。仿真结果如下：从仿真结果看，所给系统的输出y(即x1)和状态观测器输出z(即z1)的误差e在10-16内变化。因此，设计的状态观测器完全可以替代原系统。为实际中危险、高温或无法接近的那些系统的状态观测提供了一种有效途径。5.4 实验要求 (5.4) 对系统式(5.3)所示系统，采用例5.1的思路，用MATLAB编程求状态观测器的增益阵G=k1 k2 k3T ； 在SIMULINK仿真