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

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

2、ontrol toolbox )；最优化工具包(optimization toolbox )等 等。由于MATLAB功能的不断扩展，所以现在的MATLAB已不仅仅局限与现代控制系统分析和综合应用，它已是一种包罗众多学科的功能强大的“技术计算语言(The Language ofTechnical Computing)”。MathWorks公司于1992年推出了具有划时代意义的MATLAB版本，并推出了交互式模型输入与仿真系统SIMULINK,它使得控制系统的仿真与CAD应用更加方便、快捷，用户可以方便地在计算机上建模和仿真实验。1997年MathWorks推出的 MATLAB 版允许了更多的数据

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

4、容。HELP topics:matlabgeneralGeneral purpose commands.matlabopsOperators and special characters.matlablangProgramming language constructs.matlabelmatElementary matrices and matrix manipulation.matlabelfunElementary math functions.matlabspecfunSpecialized math functions.matlabmatfunMatrix functions - n

5、umerical linear algebra.:csimulinksimulinkSimulinksimulinkblocksSimulink block library.simulinksimdemos - Simulink 3 demonstrations and samples. simulinkdeeDifferential Equation EditorMATLABR11work(No table of contents file)toolboxlocalPreferences.如果用户对 MATLAB的语言结构lang感兴趣，想进一步了解，则键入:help langProgram

6、ming language constructs.Control flow.ifConditionally execute statements.elseIF statement condition.elseif - IF statement condition.endTerminate scope of FOR, WHILE, SWITCH, TRY and IF statements.forRepeat statements a specific number of times.while - Repeat statements an indefinite number of times.

7、如果想进一步了解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 variable and then the following statements, up to the END, are executed. Som

8、e examples (assume N has already beenassigned 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.normMatrix or vector norm.normest - Estimate the matrix 2-norm.rankMatrix rank.detDeterminant.

9、traceSum of diagonal elements.nullNull space.orthOrthogonalization.rrefReduced rowechelon cnsubspace - Angle between two subspaces.Eigenvalues and singular values.svdeigEigenvalues and eigenvectors.svdSingular value decomposition.gsvdGeneralized ingular value decomposition.eigsA few eigenvalues.svds

10、A few singular values.polyCharacteristic polynomial.polyeig Polynomial eigenvalue problem.condeig - Condition number with respect to eigenvalues.hessHessenberg form.qzQZ factorization for generalized eigenvalues.schurSchur decomposition.Matrix functions.expmMatrix exponential.logmMatrix logarithm.sq

11、rtmMatrix square root.funmEvaluate general matrix function.上面所列的都是有关矩阵的操作函数。如eig(A)可求出A的特征根及其特征向量，具体执行方法为：A=0 1;-6 -5 %输入 A 矩阵A =01-6-5E=eig(A) %出方阵A的特征根EE =-2 -3V,D=eig(A) %求出方阵A的特征向量 V及其A的对角型D-200-33 MATLAB基本功能我们下面给出一些 MATLAB的常用的功能，不过这只是 MATLAB及其众多TOOLBOX中 的极少极少部分。用户可以参阅有关 MATLAB的手册，或直接在 MATLAB系统中

12、用HELP命 令查阅其它功能。3. 1 MATLAB的主要线性代数运算如表A-1所示为常用的矩阵和线性代数运算函数，用户可以用help matfun获得更多内容。表A- 1常用线性代数函数B=A矩阵转置C=A+B矩阵相加C=A*B矩阵相乘C=AAk矩阵哥C=A.*B矩阵点乘，即两维数相同的矩阵各对应元素相乘expm(A)指数矩阵，也就是eAinv(A)矩阵的逆矩阵det(A)矩阵的行列式的值rank(A)计算矩阵的秩eig(A)矩阵的特征值X,D=eig(A)矩阵的特征向量X和以特征值为元素的对角阵Dp=poly(A)矩阵的特征多项式r=roots(p)特征多项式方程的根conv(p1,p2)

13、两多项式相乘2常用的控制系统处理函数TF2SS将传递函数转换到状态空间表达式A,B,C,D = TF2SS(NUM,DEN)将系统:G(s) SDEN(s)b_ m mSG(s) SDEN(s)b_ m mSbm 1m 1Sn 1an 1Sbis b0as a。转换成:X AX BUY CX DU其中:an 1an 2a1a0110000A1B000000100Cbn 1bn 2b0DbnNUM=bm,bm- 1,，b1,b0DEN=1,an-1,an-2,，a1,a0ZP2SS将零极点型传递函数转换到状态空间表达式A,B,C,D = ZP2SS(Z,R)除了 G(s) K(s Z1)(s Z

14、2) (s Zm)以外，其它与 TF2SS1 同。 (s P1)(s P2) (s Pn)SS2TF将状态空间表达式转换到传递函数NUM,DEN = SS2TF(A,B,C,D,iu)即求第iu个输入信号对输出y(t)的传递函数，即:Y(s) NUM (s)1G(s) 山 C(sI A) 1B DU iu (s) DEN (s)ns anNUM (s)n 11saisa0SS2TF的调用返回值为 G(s)的分子多项式的系数矩阵NUM和分母多项式的系数向量DEN。SS2ZP将状态空间表达式到零极点形式的传递函数的转换Z,PK = SS2ZP(A,B,C,D,iu)(5 ) TF2ZP 一般传递函

15、数转换到零极点型传递函数Z,PK = TF2ZP ( NUM,DEN)ZP2TF零极点型传递函数转换到一般传递函数NUM,DEN = ZP2TF(Z,fK)SS2SS状态空间表达式的线性变换A1,B1,C1,D1=ss2ss(A,B,C,D,T)X=TX1o所以要与用户习其中T为变换矩阵。注意变换方程为：X1=TX,X=TX1o所以要与用户习惯的变换方程一致，我们必须用T的逆代入上式，即:A1,B1,C1,D1=ss2ss(A,B,C,D,inv(T)CANON求状态空间表达式的对角标准型As,Bs,Cs,Ds,Ts=canon(A,B,C,D,mod)其中Ts为变换矩阵，注意变换方程为：Xs

16、=TsXCTRB计算系统的可控判别矩阵M=ctrb(A,B)OBSV计算系统的可观判别矩阵N=obsv(A,C)IMPULSE求系统的单位脉冲响应y,x=impulse(A,B,C,D,in,t)y,x=impulse(num,den,t)STEP求系统的单位阶跃响应y,x=step(A,B,C,D,in,t)y,x=step(num,den,t)LSIM求系统对任意输入函数u(t)的响应y,x=lsim(A,B,C,D,u,t)y,x=lsim(num,den,u,t)C2D连续系统状态方程转换为离散状态方程1为采样周期G,H=c2d(A,B,T)相关的函数还有D2C,D2D。(15) LY

17、AP求解如下形式的李雅普诺夫方程:AP PAt Q所以求解用户习惯的 ATP PAQ李雅普诺夫方程，我们必须用 A的转置A代入,即：P=Lyap (A, Q)。PLACE极点配置F=PLACE(A, B, P)PLOT绘图函数10t(X,Y stt以用不同颜色、不同符号绘制曲线，其中str可以是下列三组选项的任意组合。y yellow.pointsolidmocircle:dottedmagentac cyanxx-mark-. dashdotrred+-plus dashedg green*starb -blues-squarew -whited-diamondk -blackv-trian

18、gle(down)A triangle (up) -triangle (right)p pentagram h -hexagram4例题分析例1给定某控制系统的状态空间描述为：0110 x 611 6 x 0 u61151y 1 0 0 x试求对角规范型和变换矩阵 W,并根据其对角规范型绘制系统在初值 _一一一 x 510 15 , u 0时的响应曲线。解：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;

19、%显7K结果 disp(The Diagonal Canonical From of System is:); L b c disp(Transformation Matrix is:); W %仿真数据初始化 t=0:3; x0=5;10;15; xx0=inv(W)*x0; n=length(t);x=zeros(3,n); xx=zeros(3,n); %求解状态变量 for i=1:nxx(:,i)=expm(L*t(i)*xx0;end图plot(t,x)title(System Time Response with u=0);xlabel(Time-sec); ylabel(Res

20、ponse-value);text,x(1); text,x(3); text,x(2); The Diagonal Canonical From of System is:L =b =c =Transformation Matrix is:W =1G14121uSQ42n 1G14121uSQ42n 交由三。二(LISJzldzlJJQZ0,611.S22.63Time-sec分析：通过对角规范型求借系统状态方程的优点是,不必计算矩阵指数，只需计算单个分析：通过对角规范型求借系统状态方程的优点是,不必计算矩阵指数，只需计算单个 TOC o 1-5 h z 的标量的指数即可。在本例具体的求解的

21、过程中，需要注意的是求得对角规范型的状态xx的时间响应之后，还要用变换矩阵W将其转化为原状态变量x,然后画图。从系统的响应曲线来看，各状态变量最终都趋于零。这一方面是因为输入量u=0,另一方面是因为初值不为零，但矩阵 A的特征值均为零，系统的零输入响应呈衰减趋势。但状 态变量x (3)下降幅度和速度最大，x (1)下降幅度和速度最小。这是因为状态变量x (3)对应的特征值是-3,而x (1)对应的特征值是-1,其衰减速率不同。例2给定线性定常系统的状态空间描述为： HYPERLINK l bookmark138 o Current Document 20011 HYPERLINK l book

22、mark121 o Current Document 04130 HYPERLINK l bookmark119 o Current Document xx u HYPERLINK l bookmark140 o Current Document 00411 HYPERLINK l bookmark136 o Current Document 00022y 1 1 0 0 x试判断系统的能控性。如果系统状态完全能控，试求其能控标准型和变换矩阵。解：%系统状态方程模型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=leng

23、th(A);%求解系统能控性矩阵Q=ctrb(A,B);%系统能控性矩阵的秩m=rank(Q);处U断系统是否状态完全能控，并求解能控标准型if m=nAc1=inv(Q)*A*Q;Bc1=inv(Q)*B;Cc1=C*Q;disp(System is Controllable.);disp(System First Controllable Canonnical Form is:);Ac1Bc1Cc1disp(The Transformation Matrix is:);Qelsedisp(System State Variable cannot be totally Controlled.

24、);disp(The rank of System Controllable Matrix is:);mend运行结果：System is Controllable.System First Controllable Canonnical Form is:Ac1 = 00Bc1 =1000Cc1 = 1115826812The Transformation Matrix is:1232 TOC o 1-5 h z 46236162812024816例3对于上例的线性定常系统，试判断系统的能观性，如果系统状态完全能观的，试 求其能观标准型和变换矩阵。解：系统状态方程模型A=2 0 0 1;0 4

25、 1 3;0 0 4 1;0 0 0 2;B=1；0；1；2;C=1 1 0 0;%系统阶次n=length(A);P=obsv(A,C);% 系统能观性矩阵的秩m=rank(P);if m=nAo1=P*A*inv(P);Bo1=P*B;Co1=C*inv(P);disp(System is Observable.);disp(System First Observable Canonnical Form is:);Ao1Bo1Co1disp(The Transformation Matrix is:);inv(P)elsedisp(System State Variable cannot

26、be totally Observed.);disp(The rank of System Observable Matrix is:);end运行结果:ystem is Observable.System First Observable Canonnical Form is: TOC o 1-5 h z Ao1 =0100 HYPERLINK l bookmark144 o Current Document 00100001-6496-5212Bo1 =11158268Co1 = 1000The Transformation Matrix is:ans =例4某控制系统的状态方程描述如下。

27、试判断其稳定性并绘制其时间响应来验证上述判断。10000 x 01000001y 0 0 1 1x解：z =-1+-k =1System is stable6 5 4 3 2 10 1 6 5 4 3 2 10 1 口D,D.口 口 口 a马口口沈0|0:分析：系统的极点均在复平面的左半平面，因而系统是稳定的。 从系统的时间响应曲线上看也是如此，系统的输出振荡衰减，但可以看出，系统衰减响应曲线的振荡频率很高，衰减速率很慢。这是因为系统有两个极点的实部为,非常靠近虚轴，在实际过程中，这种系统减速率很慢。这是因为系统有两个极点的实部为,非常靠近虚轴，在实际过程中，这种系统属于性能非常差的一类。稳定

28、裕量很小，稍有干扰系统的极点就可能越过虚轴跑到复平面的右半面去。从而造成系统的不稳定。 对于这类系统，一般都要加校正手段，使其极点远里离 虚轴。s(s 1)( s 2)例s(s 1)( s 2)设计控制器u Kx k1r，使闭环系统具有极点2 j2j3, 10。解:%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

29、(a,b,p)a1=a-b*k;b1=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 = 3例6开环系统的状态方程为例6开环系统的状态方程为 TOC o 1-5 h z 0100 x001x0u61161y 1 0

30、0 x设计全维状态观测器，使观测器的闭环极点为2 j2j3, 5。解：状态观测器的结构见教材图5-5,为求出状态观测器的反馈增益阵K,先为原系统构造一对偶系统p a p c vq bp然后采用极点配置方法对对偶系统进行闭环极点培植，得到反馈增益阵k,从而可由对偶原理得到原系统的状态观测器增益阵ke 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

31、(obsv(a,c)%Construct the daul system (a1,b1,c1) a1=a;b1=c;c1=b;%slove K using poles placementp=-2+2*sqrt(3)*j -2-2*sqrt(3)*j -5;k=acker(a1,b1,p);ke=k运行结果： The rank of observability matrixr0 = 3 ke =5上机习题5-1假设系统由两个5-1假设系统由两个Gi(s)和G2(s)串联连接而成，已知Gi(s)s 1-s 3s 5Gi(s)且 G2(s)-32s2 3s 4s4 4s3s2 2s 1若想求出总系统的状态方程模型，请在MATLAB下比较下面两种方法有何不同结果:将两个传递函数模型进行串联连接，然后求出整个系统的状态方程模型。求出两个模型的状态方程表示，然后求出整个系统的状态方程模型。5-2某位置控制系统方块图如下图所示，求其状态空间表