现代控制理论MATLAB上机实验指导

1、现代控制理论matlab实践指导书刘红军 编写 现代控制理论matlab实践指导书1 matlab概述matlab是matrix laboratory的缩写，早期主要用于现代控制中复杂的矩阵、向量的各种运算。由于 matlab提供了强大的矩阵处理和绘图功能，很多专家因此在自己擅长的领域用它编写了许多专门的matlab工具包（toolbox），如控制系统工具包（control systems toolbox）；系统辨识工具包（system identification toolbox）；信号处理工具包（signal processing toolbox）；鲁棒控制工具包（robust contr

2、ol toolbox）；最优化工具包(optimization toolbox）等等。由于matlab功能的不断扩展，所以现在的matlab已不仅仅局限与现代控制系统分析和综合应用，它已是一种包罗众多学科的功能强大的“技术计算语言（the language of technical computing）”。mathworks 公司于 1992 年推出了具有划时代意义的 matlab 4.0 版本，并推出了交互式模型输入与仿真系统 simulink，它使得控制系统的仿真与 cad应用更加方便、快捷，用户可以方便地在计算机上建模和仿真实验。1997 年mathworks推出的 matlab 5.0

3、 版允许了更多的数据结构， 1999 年初推出的 matlab 5.3 版在很多方面又进一步改进了 matlab 语言的功能。2000 年底推出的 matlab 6.0。最新版本是matlab7.0。matlab以矩阵作为基本编程单元，它提供了各种矩阵的运算与操作，并有较强的绘图功能。matlab集科学计算、图像处理、声音处理于一身，是一个高度的集成系统，有良好的用户界面，并有良好的帮助功能。matlab不仅流行于控制界，在机械工程、生物工程、语音处理、图像处理、信号分析、计算机技术等各行各业中都有极广泛的应用。2 如何获得matlab帮助在matlab主窗口中键入help，即可获得第一层帮助

4、：help 加重型内容为用户键入的内容，其它为执行后显示的内容。help topics:matlabgeneral - general purpose commands.matlabops - operators and special characters.matlablang - programming language constructs.matlabelmat - elementary matrices and matrix manipulation.matlabelfun - elementary math functions.matlabspecfun - specialized

5、 math functions.matlabmatfun - matrix functions - numerical linear algebra.simulinksimulink - simulink simulinkblocks - simulink block library.simulinksimdemos - simulink 3 demonstrations and samples.simulinkdee - differential equation editormatlabr11work - (no table of contents file)toolboxlocal -

6、preferences.如果用户对matlab的语言结构lang感兴趣，想进一步了解，则键入：help langprogramming language constructs.control flow.if - conditionally execute statements.else - if statement condition.elseif - if statement condition.end - terminate scope of for, while, switch, try and if statements.for - repeat statements a specif

7、ic number of times.while - repeat statements an indefinite number of times. 如果想进一步了解for语句，则键入：help forfor repeat statements a specific number of times.the general form of a for statement is:for variable = expr, statement, ., statement endthe columns of the expression are stored one at a time in the

8、variable and then the following statements, up to the end, are executed. some examples (assume n has already been assigned a value).for i = 1:n,for j = 1:n,a(i,j) = 1/(i+j-1);endend同样，如果想了解matlab中有关矩阵的操作运算函数，可以键入：help matfunmatrix functions - numerical linear algebra.matrix analysis.norm - matrix or

9、 vector norm.normest - estimate the matrix 2-norm.rank - matrix rank.det - determinant. trace - sum of diagonal elements.null -null space.orth - orthogonalization.rref - reduced row echelon form.subspace - angle between two subspaces.eigenvalues and singular values.eig - eigenvalues and eigenvectors

10、.svd - singular value decomposition.gsvd - generalized ingular value decomposition.eigs - a few eigenvalues.svds - a few singular values.poly - characteristic polynomial.polyeig - polynomial eigenvalue problem.condeig - condition number with respect to eigenvalues.hess - hessenberg form.qz - qz fact

11、orization for generalized eigenvalues.schur - schur decomposition.matrix functions.expm - matrix exponential.logm - matrix logarithm.sqrtm - matrix square root.funm - evaluate general matrix function.上面所列的都是有关矩阵的操作函数。如eig(a)可求出a的特征根及其特征向量，具体执行方法为：a=0 1;-6 -5 输入a矩阵a = 01-6-5e=eig(a) %求出方阵a的特征根ee =-2-

12、3v,d=eig(a) %求出方阵a的特征向量v及其a的对角型dv =0.4472-0.3162-0.89440.9487d =-200-33 matlab基本功能我们下面给出一些matlab的常用的功能，不过这只是matlab及其众多toolbox中的极少极少部分。用户可以参阅有关matlab的手册，或直接在matlab系统中用help命令查阅其它功能。31 matlab的主要线性代数运算 如表a1所示为常用的矩阵和线性代数运算函数，用户可以用help matfun获得更多内容。表a1常用线性代数函数b=a矩阵转置c=a+b矩阵相加c=a*b矩阵相乘c=ak矩阵幂c=a.*b矩阵点乘，即两维

13、数相同的矩阵各对应元素相乘expm(a)指数矩阵，也就是inv(a)矩阵的逆矩阵det(a)矩阵的行列式的值rank(a)计算矩阵的秩eig(a)矩阵的特征值x,d=eig(a)矩阵的特征向量x和以特征值为元素的对角阵dp=poly(a)矩阵的特征多项式r=roots(p)特征多项式方程的根conv(p1,p2)两多项式相乘32常用的控制系统处理函数（1） tf2ss 将传递函数转换到状态空间表达式a,b,c,d = tf2ss(num,den) 将系统：转换成：其中:num=bm,bm-1,b1,b0，den=1,an-1,an-2,a1,a0 （2）zp2ss将零极点型传递函数转换到状态空

14、间表达式a,b,c,d = zp2ss(z,p,k) 除了以外，其它与tf2ss相同。（3） ss2tf 将状态空间表达式转换到传递函数num,den = ss2tf(a,b,c,d,iu) 即求第iu个输入信号对输出y(t)的传递函数，即：ss2tf的调用返回值为g(s)的分子多项式的系数矩阵num和分母多项式的系数向量den。（4） ss2zp 将状态空间表达式到零极点形式的传递函数的转换z,p,k = ss2zp(a,b,c,d,iu) （5 ）tf2zp 一般传递函数转换到零极点型传递函数z,p,k = tf2zp（num,den） （6） zp2tf 零极点型传递函数转换到一般传递函

15、数num,den = zp2tf(z,p,k) （7） ss2ss 状态空间表达式的线性变换a1,b1,c1,d1=ss2ss(a,b,c,d,t)其中t为变换矩阵。注意变换方程为：x1=tx，而不是常见的x=tx1。所以要与用户习惯的变换方程一致，我们必须用t的逆代入上式，即：a1,b1,c1,d1=ss2ss(a,b,c,d,inv(t)（8）canon 求状态空间表达式的对角标准型as,bs,cs,ds,ts=canon(a,b,c,d,mod)其中ts为变换矩阵，注意变换方程为：xs=tsx。 （9） ctrb 计算系统的可控判别矩阵mm=ctrb(a,b)（10） obsv 计算系统

16、的可观判别矩阵nn=obsv(a,c)（11） impulse 求系统的单位脉冲响应y,x=impulse(a,b,c,d,in,t)y,x=impulse(num,den,t)（12） step 求系统的单位阶跃响应y,x=step(a,b,c,d,in,t)y,x=step(num,den,t)（13） lsim 求系统对任意输入函数u(t)的响应y,x=lsim(a,b,c,d,u,t)y,x=lsim(num,den,u,t)（14） c2d 连续系统状态方程转换为离散状态方程,t为采样周期g,h=c2d(a,b,t)相关的函数还有d2c,d2d。（15） lyap 求解如下形式的李雅

17、普诺夫方程：所以求解用户习惯的李雅普诺夫方程，我们必须用a的转置a代入，即：p=lyap（a，q）。（16）place 极点配置f=place（a，b，p）（17） plot 绘图函数plot(x,y,str)可以用不同颜色、不同符号绘制曲线，其中str可以是下列三组选项的任意组合。y -yellow . -point - -solidm -magentao-circle: -dottedc -cyanx-x-mark-. -dashdot r -red + -plus- -dashed g -green* -starb -blues -squarew -white d -diamondk -

18、black v -triangle (down) -triangle (up) -triangle (right)p -pentagramh -hexagram4例题分析例1 给定某控制系统的状态空间描述为：试求对角规范型和变换矩阵w，并根据其对角规范型绘制系统在初值，时的响应曲线。解：close all;%状态方程模型a=0 1 -1; -6 -11 6;-6 -11 5; b=0;0;1; c=1 0 0;%求取对角规范型w,lamda=eig(a); l=inv(w)*a*w; b=inv(w)*b; c=c*w;%显示结果disp(the diagonal canonical from

19、 of system is:); l b c disp(transformation matrix is:);w%仿真数据初始化t=0:0.01:3; x0=5;10;15; xx0=inv(w)*x0; n=length(t);x=zeros(3,n); xx=zeros(3,n);%求解状态变量for i=1:n xx(:,i)=expm(l*t(i)*xx0; x(:,i)=w*xx(:,i);end%绘图plot(t,x)title(system time response with u=0);xlabel(time-sec); ylabel(response-value);text(

20、0.8,3.7,x(1); text(0.8,1.5,x(3); text(0.8,-0.4,x(2); the diagonal canonical from of system is:l = -1.0000 -0.0000 -0.0000 0.0000 -2.0000 -0.0000 -0.0000 0.0000 -3.0000b = -2.8284 -13.7477 10.8628c = 0.7071 -0.2182 -0.0921transformation matrix is:w = 0.7071 -0.2182 -0.0921 0.0000 -0.4364 -0.5523 0.70

21、71 -0.8729 -0.8285分析：通过对角规范型求借系统状态方程的优点是，不必计算矩阵指数，只需计算单个的标量的指数即可。在本例具体的求解的过程中，需要注意的是求得对角规范型的状态xx的时间响应之后，还要用变换矩阵w将其转化为原状态变量x，然后画图。从系统的响应曲线来看，各状态变量最终都趋于零。这一方面是因为输入量u=0，另一方面是因为初值不为零，但矩阵a的特征值均为零，系统的零输入响应呈衰减趋势。但状态变量x（3）下降幅度和速度最大，x（1）下降幅度和速度最小。这是因为状态变量x（3）对应的特征值是-3，而x（1）对应的特征值是-1，其衰减速率不同。 例2 给定线性定常系统的状态空间

22、描述为：试判断系统的能控性。如果系统状态完全能控，试求其能控标准型和变换矩阵。解：%系统状态方程模型a=2 0 0 1;0 4 1 3;0 0 4 1;0 0 0 2;b=1;0;1;2;c=1 1 0 0;%系统阶次n=length(a);%求解系统能控性矩阵q=ctrb(a,b);%系统能控性矩阵的秩m=rank(q);%判断系统是否状态完全能控，并求解能控标准型if m=n ac1=inv(q)*a*q; bc1=inv(q)*b; cc1=c*q; disp(system is controllable.); disp(system first controllable canonni

23、cal form is:); ac1 bc1 cc1 disp(the transformation matrix is:); qelse disp(system state variable cannot be totally controlled.); disp(the rank of system controllable matrix is:); mend运行结果：system is controllable.system first controllable canonnical form is:ac1 = 0 -0.0000 -0.0000 -64.0000 1.0000 -0.0

24、000 -0.0000 96.0000 0 1.0000 -0.0000 -52.0000 0.0000 0.0000 1.0000 12.0000bc1 =1 0 0 0cc1 = 1 11 58 268the transformation matrix is:q = 1 4 12 32 0 7 46 236 1 6 28 120 2 4 8 16例3 对于上例的线性定常系统，试判断系统的能观性，如果系统状态完全能观的，试求其能观标准型和变换矩阵。解：%系统状态方程模型a=2 0 0 1;0 4 1 3;0 0 4 1;0 0 0 2;b=1;0;1;2;c=1 1 0 0;%系统阶次n=l

25、ength(a);p=obsv(a,c);%系统能观性矩阵的秩m=rank(p);if m=n ao1=p*a*inv(p); bo1=p*b; co1=c*inv(p); disp(system is observable.); disp(system first observable canonnical form is:); ao1 bo1 co1 disp(the transformation matrix is:); inv(p)else disp(system state variable cannot be totally observed.); disp(the rank of

26、 system observable matrix is:); mend运行结果：ystem is observable.system first observable canonnical form is:ao1 = 0 1 0 0 0 0 1 0 0 0 0 1 -64 96 -52 12bo1 =1 11 58 268co1 = 1 0 0 0the transformation matrix is:ans = -14.0000 16.0000 -5.3750 0.5625 15.0000 -16.0000 5.3750 -0.5625 0 1.0000 -0.7500 0.1250 -

27、8.0000 8.0000 -2.5000 0.2500例4某控制系统的状态方程描述如下。试判断其稳定性并绘制其时间响应来验证上述判断。解：z =-1p =-1.4737 + 2.2638i -1.4737 - 2.2638i -0.0263 + 0.7399i -0.0263 - 0.7399ik =1system is stable分析：系统的极点均在复平面的左半平面，因而系统是稳定的。从系统的时间响应曲线上看也是如此，系统的输出振荡衰减，但可以看出，系统衰减响应曲线的振荡频率很高，衰减速率很慢。这是因为系统有两个极点的实部为-0.0263，非常靠近虚轴，在实际过程中，这种系统属于性能非常

28、差的一类。稳定裕量很小，稍有干扰系统的极点就可能越过虚轴跑到复平面的右半面去。从而造成系统的不稳定。对于这类系统，一般都要加校正手段，使其极点远里离虚轴。例5含积分环节的伺服系统设计，设对象为设计控制器，使闭环系统具有极点。解：%pole placementclose all;a=0 1 0;0 0 1;0 -2 -3;b=0;0;1;c=1 0 0;d=0;disp(the rank of controllability matrix)rc=rank(ctrb(a,b)p=-2+2*sqrt(3)*j -2-2*sqrt(3)*j -10;k=place(a,b,p)a1=a-b*k;b1=

29、b*k(1);c1=c;d1=d;figure(1)step(a1,b1,c1,d1)title(step responese of designed servo system)figure(2)y,x,t=step(a1,b1,c1,d1);plot(t,x)title(step responese curves for x1,x2,x3)grid onthe rank of controllability matrixrc = 3k = 160.0000 54.0000 11.0000例6开环系统的状态方程为设计全维状态观测器，使观测器的闭环极点为。解：状态观测器的结构见教材图5-5，为求出状态观测器的反馈增益阵k，先为原系统构造一对偶系统然后采用极点配置方法对对偶系统进行闭环极点培植，得到反馈增益阵k，从而可由对偶原理得到原系统的状态观测器增益阵。程序如下：%design of a full-order state observera=0 1 0;0 0 1;-6 -11 -6;b=0;0;1;c=1 0 0;%check observabilitydisp(the rank of observability matrix)r0=rank(obsv(a,c)%construct the daul system (