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

PAGEPAGE10《现代控制理论》MATLAB实践指导书刘红军 编写《现代控制理论》MATLAB实践指导书概述LABoratory的缩写，早期主要用于现代控制中复杂的矩阵、向量的各种运算。由于MTLAB工具包（toolbox，如控制系统工具包（controlsystemstoolbox；系统辨识工具包（systemidentificationtoolbo；信号处理工具包（signalprocessingtoolbo；鲁棒控制工具包（robustcontroltoolbo；最优化工具包(optimization等等。由于功能的不断扩展，所以现在的已不仅仅局限与现代控制系统分析和综合应用，它已是一种包罗众多学科的功能强大的“技术计算语言（TheLanguageofTechnicalComputin。MathWorks公司于1992年推出了具有划时代意义的4.0互式模型输入与仿真系统，它使得控制系统的仿真与CAD应用更加方便、快捷，用户可以方便地在计算机上建模和仿真实验。1997MathWorks推出的5.0版允许了更多的数据结构，1999的5.3版在很多方面又进一步改进了语言的功能。2000年底推出的6.0。MATLAB以矩阵作为基本编程单元，它提供了各种矩阵的运算与操作，并有较强的绘图功能。集科学计算、图像处理、声音处理于一身，是一个高度的集成系统，有良好的用户界面，并有良好的帮助功能。工程、语音处理、图像处理、信号分析、计算机技术等各行各业中都有极广泛的应用。如何获得MATLAB帮助在MATLAB主窗口中键入help，即可获得第一层帮助：help％加重型内容为用户键入的内容，其它为执行后显示的内容。HELPtopics:matlab\general Generalpurposecommands.matlab\ops Operatorsandspecialcharacters.matlab\lang Programminglanguageconstructs.matlab\elmat Elementarymatricesandmatrixmanipulation.matlab\elfun Elementarymathfunctions.matlab\specfun Specializedmathfunctions.matlab\matfun Matrixfunctions-numericallinearalgebra.simulink\simulink Simulinksimulink\blocks Simulinkblocklibrary.simulink\simdemos Simulink3demonstrationsandsamples.simulink\dee DifferentialEquationEditor(Notableofcontentsfile)toolbox\local Preferences.如果用户对MATLAB的语言结构lang感兴趣，想进一步了解，则键入：helplangProgramminglanguageconstructs.Controlflow.if Conditionallyexecutestatements.else IFstatementcondition.elseif IFstatementcondition.end TerminatescopeofFOR,WHILE,SWITCH,TRYandIFstatements.for Repeatstatementsaspecificnumberoftimes.while Repeatstatementsanindefinitenumberoftimes.如果想进一步了解for语句，则键入：helpforFORRepeatstatementsaspecificnumberoftimes.ThegeneralformofaFORstatementis:FORvariable=expr,statement, ,statementENDThecolumnsoftheexpressionarestoredoneatatimeinthevariableandthenthefollowingstatements,uptotheEND,areexecuted.……Someexamples(assumeNhasalreadybeenassignedavalue).FORI=1:N,FORJ=1:N,A(I,J)=1/(I+J-1);ENDEND同样，如果想了解MATLAB中有关矩阵的操作运算函数，可以键入：helpmatfunMatrixfunctions-numericallinearalgebra.Matrixanalysis.norm Matrixorvectornorm.normest Estimatethematrix2-norm.rank Matrixrank.det Determinant.trace Sumofdiagonalelements.null Nullspace.orth Orthogonalization.rref Reducedrowechelonform.subspace--Anglebetweentwosubspaces.Eigenvaluesandsingularvalues.eig Eigenvaluesandeigenvectors.svd Singularvaluedecomposition.gsvd Generalizedingularvaluedecomposition.eigs Afeweigenvalues.svds Afewsingularvalues.poly Characteristicpolynomial.polyeig Polynomialeigenvalueproblem.condeig Conditionnumberwithrespecttoeigenvalues.hess Hessenbergform.qz QZfactorizationforgeneralizedeigenvalues.schur Schurdecomposition.Matrixfunctions.expm Matrixexponential.logm Matrixlogarithm.sqrtm Matrixsquareroot.funm Evaluategeneralmatrixfunction.上面所列的都是有关矩阵的操作函数。如eig(A)可求出A的特征根及其特征向量，具体执行方法为：A=[01;-6-5]％输入AA0 1-6 -5E=eig(A)%求出方阵A的特征根EE=-2-3AVAV=0.4472-0.8944D=-20

-0.31620.94870-3基本功能我们下面给出一些及其众多TOOLBOX系统中用命令查阅其它功能。MATLAB的主要线性代数运算helpmatfun常用线性代数函数矩阵转置C=A+B矩阵相加C=A*B矩阵相乘C=A^k矩阵幂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)两多项式相乘常用的控制系统处理函数TF2SS将传递函数转换到状态空间表达式[A,B,C,D]=TF2SS(NUM,DEN)将系统：G(s)

NUM(s)

b smb sm1 bsbm m1 1 0DEN(s转换成：AX YCX DU其中:

sna

n

sn1 a1

sa0NUM=[bm,bm-1,…,b1,b0]，DEN=[1,an-1,an-2,…,a1,a0]a

n1

a an2 1

a 0

1 1A 00

0 01 0 0 1

0 0 0

0B 00Cbn1

bn

b

Dbn0ZP2SS将零极点型传递函数转换到状态空间表达式0[A,B,C,D]=ZP2SS(Z,P,K)K(sz)(sz )(sz )除了Gs)

(sp1

1)(sp

2 m) (sp 2 n

以外，其它与TF2SS相同。SS2TF将状态空间表达式转换到传递函数[NUM,DEN]=SS2TF(A,B,C,D,iu)即求第iu个输入信号对输出y(t)的传递函数，即：G(s)

Y(s)

NUM(s) C(sIA)1B

NUM(s) U (s)iu

DEN(s)

sna

n1

sn1 a1

sa0SS2TF的调用返回值为G(s)的分子多项式的系数矩阵NUM和分母多项式的系数向量DEN。SS2ZP将状态空间表达式到零极点形式的传递函数的转换[Z,P,K]=SS2ZP(A,B,C,D,iu)（5）TF2ZP一般传递函数转换到零极点型传递函数[Z,P,K]=TF2ZP（NUM,DEN）ZP2TF零极点型传递函数转换到一般传递函数[NUM,DEN]=ZP2TF(Z,P,K)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))CANON求状态空间表达式的对角标准型[As,Bs,Cs,Ds,Ts]=canon(A,B,C,D,'mod')其中Ts为变换矩阵，注意变换方程为：Xs=TsX。CTRBM=ctrb(A,B)OBSVN=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连续系统状态方程转换为离散状态方程,T为采样周期[G,H]=c2d(A,B,T)相关的函数还有D2C,D2D。LYAPAPPATQATPPAQ代入，即：P=Lya（’，。PLACE极点配置F=PLACE（A，B，P）PLOT绘图函数plot(X,Y,’str’)可以用不同颜色、不同符号绘制曲线，其中‘str’可以是下列三组选项的任意组合。y yellowm magentac cyanr redg greenb bluew whitek black4

.pointo circlex x-mark+ plus* stars squared diamondv triangle(down)^ triangle(up)< triangle(left)> triangle(right)p pentagramh hexagram

-solid:dotted-.dashdot--dashed[例1]给定某控制系统的状态空间描述为：0

1 0 x66y

11110

6x5

0u1试求对角规范型和变换矩阵W，并根据其对角规范型绘制系统在初值x10 15u0解：closeall;%状态方程模型A=[01-1;-6-116;-6-115]; B=[0;0;1]; C=[100];%求取对角规范型[W,lamda]=eig(A); L=inv(W)*A*W; b=inv(W)*B; c=C*W;%显示结果disp('TheDiagonalCanonicalFromofSystemis:');Lbcdisp('TransformationMatrix W%仿真数据初始化t=0:0.01:3; x0=[5;10;15];xx0=inv(W)*x0; x=zeros(3,n); xx=zeros(3,n);%求解状态变量fori=1:nxx(:,i)=expm(L*t(i))*xx0;x(:,i)=W*xx(:,i);end%绘图plot(t,x)title('SystemTimeResponsewithu=0');xlabel('Time-sec'); ylabel('Response-value');text(0.8,3.7,'x(1)');text(0.8,1.5,'x(3)');text(0.8,-0.4,'x(2)');>>TheDiagonalCanonicalFromofSystemis:L=-1.0000-0.0000-0.00000.0000-2.0000-0.0000-0.00000.0000-3.0000b=-2.8284-13.747710.8628c=0.7071-0.2182-0.0921TransformationMatrix is:W=0.7071 -0.2182 -0.09210.0000-0.4364-0.55230.7071-0.8729-0.8285的标量的指数即可。在本例具体的求解的过程中，需要注意的是求得对角规范型的状态xxWx，然后画图。从系统的响应曲线来看，各状态变量最终都趋于零。这一方面是因为输入量u=0，另一方面是因为初值不为零，但矩阵A的特征值均为零，系统的零输入响应呈衰减趋势。但状态变量x（3）下降幅度和速度最大，x（1）下降幅度和速度最小。这是因为状态变量x（3）对应的特征值是-3，而x（1）对应的特征值是-1，其衰减速率不同。[例2]给定线性定常系统的状态空间描述为：2 0 0 4 x

1 13 0x u0 0

4 1

1y

0 0 2 21 0 0试判断系统的能控性。如果系统状态完全能控，试求其能控标准型和变换矩阵。解：%系统状态方程模型A=[2001;0413;0041;000B=[1;0;1;2];C=[1100];%系统阶次n=length(A);%求解系统能控性矩阵Q=ctrb(A,B);%系统能控性矩阵的秩m=rank(Q);%判断系统是否状态完全能控，并求解能控标准型ifm==nAc1=inv(Q)*A*Q;elseend

Bc1=inv(Q)*B;Cc1=C*Q;disp('SystemisControllable.');disp('SystemFirstControllableCanonnicalFormAc1Bc1Cc1disp('TheTransformationMatrixis:');Qdisp('SystemStateVariablecannotbetotallyControlled.');disp('TherankofSystemControllableMatrixis:');m运行结果：SystemisControllable.SystemFirstControllableCanonnicalFormis:Ac1=0 -0.0000 -0.0000 -64.00001.0000 -0.0000 -0.0000 96.00000 1.0000 -0.0000 -52.00000.00000.00001.000012.0000Bc1=1000Cc1=11158268TheTransformationMatrixis:Q= 1412320746236162812024816例3]求其能观标准型和变换矩阵。解：%系统状态方程模型A=[2001;0413;0041;000B=[1;0;1;2];C=[1100];%系统阶次n=length(A);P=obsv(A,C);%系统能观性矩阵的秩m=rank(P);ifm==nAo1=P*A*inv(P);Bo1=P*B;elseend

Co1=C*inv(P);disp('SystemisObservable.');disp('SystemFirst ObservableCanonnicalFormis:');Ao1Bo1Co1disp('TheTransformationMatrixis:');inv(P)disp('SystemStateVariablecannotbetotallyObserved.');disp('TherankofSystemObservableMatrixis:');m运行结果：ystemisObservable.SystemFirst ObservableCanonnicalFormis:Ao1=010000100001-64 96 -52 12Bo1=11158268Co1=1 0 0 0TheTransformationMatrixis:ans=-14.0000 16.0000 -5.3750 0.562515.0000-16.00005.3750-0.562501.0000-0.75000.1250-8.00008.0000-2.50000.2500[例4]某控制系统的状态方程描述如下。试判断其稳定性并绘制其时间响应来验证上述判断。3 8 2 4 11 0 x

0 0x u0 1 0 0

0 0 0 0y0 1 解：z=-1p=-1.4737+2.2638i-1.4737-2.2638i-0.0263+0.7399i-0.0263-0.7399ik=1Systemisstable-0.0263，非常靠近虚轴，在实际过程中，例5]含积分环节的伺服系统设计，设对象为G(s) 1 s(s1)(s2)设计控制器uKx解：%Poleplacementcloseall;a=[010;001;0-2-3];b=[0;0;1];c=[100];d=[0];

r，使闭环系统具有极点2j2 3,10。1disp('Therankofcontrollabilitymatrix')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=b*k(1);c1=c;d1=d;figure(1)step(a1,b1,c1,d1)title('StepResponeseofDesignedServoSystem')figure(2)[y,x,t]=step(a1,b1,c1,d1);plot(t,x)title('StepResponeseCurvesforx1,x2,x3')gridonTherankofcontrollabilitymatrixrc=3k= 160.0000 54.0000 11.0000 [例6]开环系统的状态方程为0 1 0 00 0 1x

0u 6y

110

6

1设计全维状态观测器，使观测器的闭环极点为2j2 3,5。解：状态观测器的结构见教材图5-5，为求出状态观测器的反馈增益阵K，先为原系统构造一对偶系统apcqbp然后采用极点配置方法对对偶系统进行闭环极点培植，得到反馈增益阵kkek。程序如下：%Designofafull-orderstateobservera=[010;001;-6-11-6];b=[0;0;1];c=[100];%Checkobservabilitydisp('Therankofobservabilitymatrix')r0=rank(obsv(a,c))%Constructthedaulsystem(a1,b1,c1)a1=a';b1=c';c1=b';%sloveKusingpolesplacementp=[-2+2*sqrt(3)*j-2-2*sqrt(3)*j-5];k=acker(a1,b1,p);ke=k'运行结果：Therankofobservabilityr0=3ke=3.00007.0000-1.00005上机习题1假设系统由两个1

s)和G2s)串联连接而成，已知G(s) s1 且

s3s5 )1 s23s4 2

s44s33s22s1若想求出总系统的状态方程模型，请在MATLAB下比较下面两种方法有何不同结果：将两个传递函数模型进行串联连接，然后求出整个系统的状态方程模型。求出两个模型的状态方程表示，然后求出整个系统的状态方程模型。某位置控制系统方块图如下图所示，求其状态空间表达式。5-3试将5-2题所得动态方程转换成Jordan2=54语言求解。5-4k=2，k1=－0.05，k2=0.16，k3=.245-2所示系统的单位阶跃响应。5-5判断