现代控制理论实验报告

实验报告 ( 20１6-2017 年度第二学期) 名

称:《现代控制理论基础》

题

目:状态空间模型分析 院

系:控制科学与工程学院

班

级:＿__

学

号:＿_

学生姓名:__＿___

指导教师:______＿

成

绩:

日期: 20１7 年 4 月 1５日

线控实验报告

一 、实验目得: :

l。加强对现代控制理论相关知识得理解; 2、掌握用 matlaｂ 进行系统李雅普诺夫稳定性分析、能控能观性分析; 二、实验内容

第一题:已知某系统得传递函数为

求解下列问题: (1)用 mａtlab 表示系统传递函数

ｎum＝［1];

dｅn＝［1 3 2];

sys=ｔf(num,dｅn);

ｓys1=zpｋ([],[-１ －2］,1); 结果:

syｓ =

1

—－---－-－——-－—

s^2 ＋ 3 s ＋ ２

sys1 ＝

1

-－——－——--——

(s＋１) (ｓ+2) (２)求该系统状态空间表达式: [A1,Ｂ1,C1,D1]=tf2ss(nuｍ,ｄen); Ａ ＝

-3

—2

1

0 B =

1

０ Ｃ =

0

1 第二题:已知某系统得状态空间表达式为::求解下列问题: （1）求该系统得传递函数矩阵: （2）该系统得能观性与能空性: （3）求该系统得对角标准型: （4）求该系统能控标准型: （5）求该系统能观标准型:

（6）求该系统得单位阶跃状态响应以及零输入响应: 解题过程: 程序:Ａ=[—３ －２;1 0］;B=［1 ０]＇;C=[０ 1];D=0; ［ｎuｍ,den］=ss２ｔｆ(A,B,Ｃ,Ｄ); ｃo=ctｒb(A,Ｂ); t１＝ｒanｋ(co); ob=ｏbｓv(Ａ,C); t2＝raｎk(ob); [Aｔ,Bt,Ct,Ｄt,T］=canｏn(A,B,Ｃ,Ｄ,＇ｍodal’); ［Ａｃ,Bc,Cc,Dｃ,Ｔｃ]=cａｎoｎ(Ａ,B,C,D,"cｏmｐａnｉon＇); Ao=Ａc"; Ｂo＝Cc"; Ｃo=Bc＇; 结果: (1)nｕm ＝

０

0

１ den =

1

3

２ (2)能控判别矩阵为: cｏ =

１

—3

0

1 能控判别矩阵得秩为: ｔ１ ＝

2

故系统能控。

(3)能观判别矩阵为: ｏb ＝

0

1

1

0 能观判别矩阵得秩为: t2 =

2 故该系统能观、(4)该系统对角标准型为: At =

－2

0

0

-1 Bt =

-1、414２

－1、１1８０ Ct =

0。70７1

-0.８944 (5)该系统能观标准型为:

Ao ＝

0

-2

1

-3 Bｏ ＝

1

0 Cｏ =

0

１ (6)该系统能控标准型为: Ａc =

０

1 -2

－３ Bc =

0

1 Ｃc =

1

０ (７)系统单位阶跃状态响应; G＝sｓ(A1,Ｂ１,C１,Ｄ１); ［y,t,x］=sｔｅp(G); fiｇure(1) plot(t,x);

(8)零输入响应: x0＝［0 1];

［y,t,ｘ］=inｉｔial(Ｇ,x０); fiｇuｒe(2) plot(t,x)

第三题:已知某系统得状态空间模型各矩阵为: ,求下列问题: （1）按能空性进行结构分解: （2）按能观性进行结构分解: ｃlear

Ａ=[0 0 -１;１ 0 —3;0 １ -３］; Ｂ=[1 1 0］"; C=[０ 1 -2］; tｃ＝rａnk(ctｒb(A,Ｂ)); to=raｎk(oｂｓｖ(A,C)); ［A１,B1,C1,ｔ1,k1］=ctrｂf(A,B,C); ［Ａ２,B２,C2,t２,k２]=ctrbf(A,Ｂ,C); 结果: 能控判别矩阵秩为: tc =

2 可见,能空性矩阵不满秩,系统不完全能控。

A1 =

－1、0０00

－０、0000

—0.0000

２。１213

-2。5000

0、86６０

1.２247

—2。5981

０、5０００

B１ =

0。0000

0.0000 1。4142 Ｃ1 =

1、7３２１

-1.２247

０。７０71 t1 ＝

-0、５7７４

0、5774

—０、5774

-０、40８2

0、4０８2

0、８16５ 0.70７１

0、707１

0 k1 =

1

1

0 能观性判别矩阵秩为: to =

2 可见,能观性判别矩阵不满秩,故系统不完全能观。

A2 =

－1、0000

1、3４16

3、8341

0.000０

—０。400０

—0。7348 0。00０0

0。48９9

－1、6０００ B2 =

1。2２47

0。5477 0。4４72 C2 =

0

－０。０00０

2。236１ t2 ＝

0、4082

0.8165

0、４０82

0、９1２9

-０.3651

-０.18２6

0

0、4472

-０、8944 k2 =

1

1

0 第四题:已知系统得状态方程为:

希望极点为—2,－3,-4.试设计状态反馈矩阵Ｋ,并比较状态反馈前后输出响应。

A=[1 2 ３;4 5 6;7 8 9］; B=［０ 0 1]＇; C=[0 1 0］; D=0; tc=rａnk(ctｒb(A,B)); ｐ=［—２ -３ -4］; K=ｐｌace(A,B,p); t=0:0.01:5; U=0。025*ｏneｓ(siｚｅ(t));

[Y1,Ｘ１]＝lsiｍ(A,B,C,D,U,t); ［Y2,X２］=lsｉm(A－B＊K,B,C,D,U,t); ｆigure(1) ｐｌｏt(t,Y1); grｉd oｎ title(’反馈前"); fｉguｒe(2) plot(t,Y２) title(’反馈后") 结果: ｔc =

3 可见,能观判别矩阵满秩,故系统能进行任意极点配置。

反馈矩阵为: Ｋ =

1５。３333

23、6667

２4.000０ 反馈前后系统输出对比:

第五题。已知某线性定常系统得系统矩阵为:,判断该系统稳定性。

clｅａr

clc Ａ=［－１ １;２ -３］; Ａ=Ａ’; Q=eye(2); P＝lyａp(Ａ,Q); det(P); 结果: 求得得 P 矩阵为: Ｐ =

1、７5０0

0、6２5０ ０.６250

0。３75０ 且Ｐ阵得行列式为: 〉＞ ｄeｔ(Ｐ) ａns = 0。2656 可见,P 矩阵各阶主子行列式均大于 0,故 P 阵正定,故该系统稳定、

实验一 线性定常系统模型

一 实验目的

1.掌握线性定常系统的状态空间表达式。学会在MATLAB中建立状态空间模型的方法。 2.掌握传递函数与状态空间表达式之间相互转换的方法。学会用MATLAB实现不同模型之间的相互转换。

3.熟悉系统的连接。学会用MATLAB确定整个系统的状态空间表达式和传递函数。

4.掌握状态空间表达式的相似变换。掌握将状态空间表达式转换为对角标准型、约当标准型、能控标准型和能观测标准型的方法。学会用MATLAB进行线性变换。

二 实验原理

1.线性定常系统的数学模型

在MATLAB中，线性定常（linear time invariant, 简称为 LTI）系统可以用4种数学模型描述，即传递函数(TF)模型、零极点增益(ZPK)模型和状态空间(SS)模型以及SIMULINK结构图。前三种数学模型是用数学表达式表示的，且均有连续和离散两种类型，通常把它们统称为LTI模型。

1) 传递函数模型（TF 模型）

令单输入单输出线性定常连续和离散系统的传递函数分别为

Y(s)bmsmbmsmb1sb0

(1-1) G(s)nU(s)san1sn1a1sa0和

Y(z)bmzmbmzmb1zb0。

(1-2) G(z)nn1U(z)zan1za1za0在MATLAB中，连续系统和离散系统的传递函数都用分子/分母多项式系数构成的两个行向量num和den表示，即

numbmb1b0，den1an1a0

系统的传递函数模型用MATLAB提供的函数tf( )建立。函数tf ( )不仅能用于建立系统传递函数模型，也能用于将系统的零极点增益模型和状态空间模型转换为传递函数模型。该函数的调用格式如下：

,de)n 返回连续系统的传递函数模型G。

Gtf(num

Gtf(num,den,Ts) 返回离散系统的传递函数模型G。Ts为采样周期，当Ts=-1或者Ts=[]时，系统的采样周期未定义。

Gtftf(G) 可将任意的LTI模型G转换为传递函数模型Gtf。

2) 零极点增益模型（ZPK模型）

系统的零极点增益模型是传递函数模型的一种特殊形式。令线性定常连续和离散系统的零极点形式的传递函数分别为

G(s)

(sz1)(sz2)(szm)Y(s) (1-3) KU(s)(sp1)(sp2)(spn)

和

G(z)(zz1)(zz2)(zzm)Y(z) (1-4) KU(z)(zp1)(zp2)(zpn)在MATLAB中，连续和离散系统的零点和极点都用行向量z和p表示，即

zz1z2zm，pp1p2pn。

系统的零极点增益模型用MATLAB提供的函数zpk ( )建立。函数zpk( )不仅能用来建立系统零极点增益模型，也能用于将系统的传递函数模型和状态空间模型转换为零极点增益模型。该函数的调用格式如下：

Gzpk(z,p,k) 返回连续系统的零极点增益模型G。

Gzpk(z,p,k,Ts) 返回离散系统的零极点增益模型G。Ts为采样周期，当Ts=-1或者Ts=[]时，系统的采样周期未定义。

Gzpkzpk(G) 可将任意的LTI模型G转换为零极点增益模型Gzpk。 3) 状态空间模型(SS模型) 令多输入多输出线性定常连续和离散系统的状态空间表达式分别为

(t)Ax(t)Bu(t) xy(t)Cx(t)Du(t) （1-5）

和

x(k1)Ax(k)Bu(k)

y(k)Cx(k)Du(k) （1-6）

在MATLAB中，连续系统和离散系统的状态空间模型都用MATLAB提供的函数 ( )建立。函数 ( )不仅能用于建立系统的状态空间模型，也能用于将系统的传递函数模型和零极点增益模型转换为状态空间模型。该函数的调用格式如下：

G(A,B,C,D) 返回连续系统的状态空间模型G。

G(A,B,C,D,Ts) 返回离散系统的状态空间模型G。Ts为采样周期，当Ts=1或者Ts=[]时，系统的采样周期未定义。

G(G) 可将任意的LTI模型G转换为状态空间模型G。

2．模型转换

上述三种LTI模型之间可以通过函数tf( )，zpk( )和( )相互转换。线性定常系统的传递函数模型和零极点增益模型是唯一的，但系统的状态空间模型是不唯一的。函数( )只能将传递函数模型和零极点增益模型转换为一种指定形式的状态空间模型。

三 实验内容

1.已知系统的传递函数

s26s84(a) G(s) (b) G(s)2 2s4s3s(s1)(s3)（1）建立系统的TF或ZPK模型。

（2）将给定传递函数用函数( )转换为状态空间表达式。再将得到的状态空间表达式用函数tf( )转换为传递函数，并与原传递函数进行比较。 解：（a）G(s)4

s(s1)2(s3)（1）TF模型

在命令窗中运行下列命令

>> num=4;den=[1 5 7 3];G=tf(num,den)

Transfer function:

4 --------------------- s^3 + 5 s^2 + 7 s + 3

ZPK模型

在命令窗中运行下列命令

>> z=[];p=[0 -1 -1 -3];k=4;G=zpk(z,p,k)

Zero/pole/gain:

4 --------------- s (s+1)^2 (s+3)

（2）在命令窗中运行下列命令

>> num=4;den=[1 5 7 3];Gtf=tf(num,den); >> G=(Gtf)

a =

x1

x2

x3

x1

-5

-0.875 -0.09375

x2

0

0

x3

0

0

b =

u1

x1 0.25

x2

0

x3

0

c =

x1

x2

x3

y1

0

0 0.5

d =

u1

y1

0

Continuous-time model.>> Gtf1=tf(G)

Transfer function:

4 --------------------- s^3 + 5 s^2 + 7 s + 3

s26s8（b）G(s)2

s4s3（1）TF模型

在命令窗中运行下列命令

>> num=[1 6 8];den=[1 4 3];G=tf(num,den)

Transfer function: s^2 + 6 s + 8 ------------- s^2 + 4 s + 3

ZPK模型

在命令窗中运行下列命令

>> z=[-2 -4];p=[-1 -3];k=1;G=zpk(z,p,k)

Zero/pole/gain: (s+2) (s+4) ----------- (s+1) (s+3)

(2) 在命令窗中运行下列命令

>> num=[1 6 8];den=[1 4 3];Gtf=tf(num,den); >> G=(Gtf)

a =

x1

x2

x1

-4 -0.75

x2

0

b =

u1

x1

x2

0

c =

x1

x2

y1

1 0.625

d =

u1

y1

1

Continuous-time model.>> Gtf1=tf(G)

Transfer function: s^2 + 6 s + 8 ------------- s^2 + 4 s + 3

2.已知系统的状态空间表达式

100(a) xx1u y11x

56

1002x1u y111x 302 (b) x12767（1）建立给定系统的状态空间模型。用函数eig( ) 求出系统特征值。用函数tf( ) 和zpk( )将这些状态空间表达式转换为传递函数，记录得到的传递函数和它的零极点。比较系统的特征值和极点是否一致，为什么? （2）用函数canon( )将给定状态空间表达式转换为对角标准型。用函数eig( )求出系统特征值。比较这些特征值和（1）中的特征值是否一致，为什么? 再用函数tf( )和zpk( )将

对角标准型或约当标准型转换为传递函数。比较这些传递函数和（1）中的传递函数是否一致，为什么? 解： (a) x100xu

y11x

561（1）在命令窗中运行下列命令

>> A=[0 1;-5 -6];B=[0;1];C=[1 1];D=0;G=(A,B,C,D)

a =

x1 x2

x1

0

x2 -5 -6

b =

u1

x1

0

x2

c =

x1 x2

y1

d =

u1

y1

0

Continuous-time model.>> Geig=eig(G)

Geig =

-3

-1

>> Gtf=tf(G)

Transfer function: s^2 + 6 s + 8 ------------- s^2 + 4 s + 3

>> Gzpk=zpk(G)

Zero/pole/gain: (s+4) (s+2) ----------- (s+3) (s+1)

分析：z=-4,-2;p=-3,-1 系统的特征值和极点一致。

（2）在命令窗中运行下列命令

>> A=[0 1;-5 -6];B=[0;1];C=[1 1];D=0;G=(A,B,C,D);GJ=canon(G,\\"model\\")

a =

x1 x2

x1 -1

0

x2

0 -5

b =

u1

x1 0.3536

x2

1.275

c =

x1

x2

y1

0 0.7845

d =

u1

y1

0

Continuous-time model.>> Geig=eig(GJcanon) ??? Undefined function or variable \\"GJcanon\\".

>> A=[0 1;-5 -6];B=[0;1];C=[1 1];D=0;G=(A,B,C,D); >> Gcanon=canon(G)

a =

x1 x2

x1 -3

0

x2

0 -1

b =

u1

x1

-5

x2 -4.123

c =

x1

x2

y1

-0.1 -0.3638

d =

u1

y1

1

Continuous-time model.>> Geig=eig(Gcanon)

Geig =

-3

-1

>> Gtf=tf(Gcanon)

Transfer function: s^2 + 6 s + 8 ------------- s^2 + 4 s + 3

>> Gzpk=zpk(Gcanon)

Zero/pole/gain: (s+4) (s+2) ----------- (s+3) (s+1) 分析：这些特征值和（1）中的特征值一致；这些传递函数和（1）中的传递函数一致。

1002x1u y111x 302(b) x12767

（1）在命令窗中运行下列命令

>> A=[0 1 0;3 0 2;-12 -7 -6];B=[2;1;7];C=[1 1 1];D=0;G=(A,B,C,D)

a =

x1

x2

x3

x1

0

0

x2

0

x3 -12

-7

-6

b =

u1

x1

x2

x3

c =

x1 x2 x3

y1

d =

u1

y1

0

Continuous-time model.>> Geig=eig(G)

Geig =

-3

-1

>> Gtf=tf(G)

Transfer function: s^2 + 6 s + 8 ------------- s^2 + 4 s + 3

>> Gzpk=zpk(G)

Zero/pole/gain:

(s+4) (s+2) ----------- (s+3) (s+1)

(2) >> A=[0 1 0;3 0 2;-12 -7 -6];B=[2;1;7];C=[1 1 1];D=0;G=(A,B,C,D)

a =

x1

x2

x3

x1

0

0

x2

0

x3 -12

-7

-6

b =

u1

x1

x2

x3

c =

x1 x2 x3

y1

d =

u1

y1

0

Continuous-time model.

>> Geig=eig(Gcanon)

Geig =

-3

-1

>> Gtf=tf(Gcanon)

Transfer function: s^2 + 6 s + 8 ------------- s^2 + 4 s + 3

>> Gzpk=zpk(Gcanon)

Zero/pole/gain: (s+4) (s+2) ----------- (s+3) (s+1)

>> A=[0 1 0;3 0 2;-12 -7 -6];B=[2;1;7];C=[1 1 1];D=0;G=(A,B,C,D)

a =

x1

x2

x3

x1

0

0

x2

0

x3 -12

-7

-6

b =

u1

x1

x2

x3

c =

x1 x2 x3

y1

d =

u1

y1

0

Continuous-time model.>> Geig=eig(G)

Geig =

-3

-1

>> Gtf=tf(G)

Transfer function:

s^2 + 6 s + 8 ------------- s^2 + 4 s + 3

>> Gzpk=zpk(G)

Zero/pole/gain: (s+4) (s+2) ----------- (s+3) (s+1)

>> Geig=eig(Gcanon)

Geig =

-3

-1

>> Gtf=tf(Gcanon)

Transfer function: s^2 + 6 s + 8 ------------- s^2 + 4 s + 3

>> Gzpk=zpk(Gcanon)

Zero/pole/gain: (s+4) (s+2) ----------- (s+3) (s+1)

四．实验总结

1.通过实验，掌握了线性定常系统的状态空间表达式、传递函数与状态空间表达式之间相互转换的方法、状态空间表达式的相似变换、将状态空间表达式转换为对角标准型、约当标准型。

2.学会在MATLAB中建立状态空间模型的方法、实现不同模型之间的相互转换、进行线性变换。

实验二 线性定常系统状态方程的解

一、实验目的

1.掌握状态转移矩阵的概念。学会用MATLAB求解状态转移矩阵。 2.掌握线性系统状态方程解的结构。学会用MATLAB求解线性定常系统的状态响应和输出响应，并绘制相应曲线。

二 实验原理

1、线性定常连续系统状态转移矩阵的计算

线性定常连续系统的状态转移矩阵为(t)eAtL1[(sIA)1]。 （3-2-1） 在MATLAB中, 状态转移矩阵可直接用指数矩阵法和拉氏反变换法计算。 2.线性定常连续系统的状态方程求解

如果线性定常连续系统的状态空间表达式为

AxBu xyCxDu

且初始状态为x(0)，那么状态方程解的拉氏变换式为

x(s)(sIA)1x(0)(sIA)1Bu(s)

（3-2-2）

其解为

tx(t)ex(0)eA(t)Bu()d

（3-2-3） At0其中零输入响应为

ex(0)或L{(sIA)}x(0)

（3-2-4） 零状态响应为

At11t0eA(t)Bu()d或L1{(sIA)1Bu(s)}

（3-2-5）

111系统的输出响应为

L{C(sIA)x(0)C(sIA)Bu(s)}Du(t)

（3-2-6）

三、实验内容

1.求下列系统矩阵A对应的状态转移矩阵

0100001001 （c）A00 （a）A （b）A4025400解：（a）A01 40指数矩阵法：

在命令窗中运行下列命令

>> A=[0 -1;4 0];syms t;phet=expm(A*t)

phet =

[

cos(2*t), -1/2*sin(2*t)] [

2*sin(2*t),

cos(2*t)]

拉氏反变换法：

在命令窗中运行下列命令

>> A=[0 -1;4 0];syms s;G=inv(s*eye(size(A))-A)

G =

[ s/(s^2+4), -1/(s^2+4)] [ 4/(s^2+4), s/(s^2+4)] 即(sIA)1。再对其进行拉氏逆变换，即在命令窗中输入语句 >> phet=ilaplace(G)

phet =

[

cos(4^(1/2)*t), -1/4*4^(1/2)*sin(4^(1/2)*t)] [

4^(1/2)*sin(4^(1/2)*t),

cos(4^(1/2)*t)]

010（b）A001 254指数矩阵法：

在命令窗中运行下列命令

>> A=[0 1 0;0 0 1;2 -5 4];syms t;phet=expm(A*t)

phet =

[

-2*t*exp(t)+exp(2*t), exp(2*t)-exp(t)-t*exp(t)] [ 2*exp(2*t)-2*exp(t)-2*t*exp(t), 2*exp(2*t)-2*exp(t)-t*exp(t)] [ -2*t*exp(t)+4*exp(2*t)-4*exp(t), -3*exp(t)+4*exp(2*t)-t*exp(t)]

拉氏反变换法：

在命令窗中运行下列命令

>> A=[0 1 0;0 0 1;2 -5 4];syms s;G=inv(s*eye(size(A))-A)

-2*exp(2*t)+2*exp(t)+3*t*exp(t),

5*exp(t)+3*t*exp(t)-4*exp(2*t),

-8*exp(2*t)+8*exp(t)+3*t*exp(t),

G =

[ (s^2-4*s+5)/(s^3-4*s^2+5*s-2),

(s-4)/(s^3-4*s^2+5*s-2),

1/(s^3-4*s^2+5*s-2)] [

2/(s^3-4*s^2+5*s-2),

s*(s-4)/(s^3-4*s^2+5*s-2),

s/(s^3-4*s^2+5*s-2)] [

2*s/(s^3-4*s^2+5*s-2),

-(5*s-2)/(s^3-4*s^2+5*s-2),

s^2/(s^3-4*s^2+5*s-2)]

即(sIA)1。再对其进行拉氏逆变换，即在命令窗中输入语句 >> phet=ilaplace(G)

phet =

[

-2*t*exp(t)+exp(2*t), exp(2*t)-exp(t)-t*exp(t)] [ 2*exp(2*t)-2*exp(t)-2*t*exp(t), 2*exp(2*t)-2*exp(t)-t*exp(t)] [ -2*t*exp(t)+4*exp(2*t)-4*exp(t), 4*exp(2*t)-t*exp(t)-3*exp(t)]

-2*exp(2*t)+2*exp(t)+3*t*exp(t),

-4*exp(2*t)+3*t*exp(t)+5*exp(t),

-8*exp(2*t)+3*t*exp(t)+8*exp(t),

00（c）A00

00指数矩阵法：

在命令窗中运行下列命令

>> A=[3 0 0;0 3 0;0 0 3];syms t;phet=expm(A*t)

phet =

[ exp(3*t),

0,

0] [

0, exp(3*t),

0] [

0,

0, exp(3*t)]

拉氏反变换法：

在命令窗中运行下列命令

>> A=[3 0 0;0 3 0;0 0 3];syms s;G=inv(s*eye(size(A))-A)

G =

[ 1/(s-3),

0,

0] [

0, 1/(s-3),

0]

[

0,

0, 1/(s-3)]

即(sIA)1。再对其进行拉氏逆变换，即在命令窗中输入语句 >> phet=ilaplace(G)

phet =

[ exp(3*t),

0,

0] [

0, exp(3*t),

0] [

0,

0, exp(3*t)]

2.已知系统

100xxu y10x 651（1）令初始状态为x(0)，输入为零。

a) 用MATLAB求状态方程的解析解。选择时间向量t，绘制系统的状态响应曲线。观察并记录这些曲线。

b) 用函数initial( )计算系统在初始状态作用下状态响应和输出响应的数值解, 并用函数plot( ) 绘制系统的状态响应曲线和输出响应曲线。观察并记录这些响应曲线，然后将这一状态响应曲线与a)中状态响应曲线进行比较。 (2) 令初始状态为零，输入为u(t)1(t)。

a) 用MATLAB求状态方程的解析解。选择时间向量t，绘制系统的状态响应曲线。观察并记录这些曲线。

b) 用函数initial( )计算系统在初始状态作用下状态响应和输出响应的数值解, 并用函数plot( ) 绘制系统的状态响应曲线和输出响应曲线。观察并记录这些响应曲线，然后将这一状态响应曲线与a).中状态响应曲线进行比较。

101（3）令初始状态为x(0)，输入为u(t)1(t)。求系统状态响应和输出响应的数值

1解，绘制系统的状态响应曲线、输出响应曲线和状态轨迹。观察和分析这些响应曲线和状态轨迹是否是（1）和（2）中的响应曲线和状态轨迹的叠加。

解：x100x1u y10x

6510（1）令初始状态为x(0)，输入为零

（a）编制程序%ex22求输入为零时状态方程的解。该程序如下：

>> A=[0 1;-6 -5];syms s;G=inv(s*eye(size(A))-A);phet=ilaplace(G);X0=[1 0]\\";Xt1=phet*X0

Xt1 =

[ -2*exp(-3*t)+3*exp(-2*t)] [ -6*exp(-2*t)+6*exp(-3*t)]

>> B=[0 1]\\";Xt2=ilaplace(G*B*1)

Xt2 =

[

exp(-2*t)-exp(-3*t)] [ 3*exp(-3*t)-2*exp(-2*t)] 其中xt1为零输入响应，xt2为零状态响应。 上述得到的是状态方程的解析解。

状态响应曲线：

（b）在命令窗中运行下列命令，建立状态空间模型，计算系统在初始状态作用下的状态响应和输出响应，并绘制相应的响应曲线。

>> A=[0 1;-6 -5];B=[0;1];C=[1 0];D=0;G=(A,B,C,D); >> t=0:0.5:10;x0=[1;0]; >> [yo,t,xo]=initial(G,x0,t);plot(t,xo,\\":\\",t,yo,\\"-\\") 返回图1

图1状态响应

在命令窗中继续运行下列命令，计算系统在输入作用下的状态响应和输出响应，并绘制相应的响应曲线。

>> figure(\\"pos\\",[50 50 200 150],\\"color\\",\\"w\\"); u=ones(size(t)); [yu,t,xu]=lsim(G,u,t); plot(t,xu,\\":\\",t,yu,\\"-\\") 返回图2。

图2 输出响应

再继续运行下列命令求系统总的状态响应和输出响应，并绘制相应的响应曲线。 >>y=yo+yu; x=xo+xu; plot(t,x,\\":\\",t,y,\\"-\\") 返回图3。

图3

（2）令初始状态为零，输入为u(t)1(t)。

编制程序%ex22求输入为u(t)1(t)时状态方程的解。该程序如下：

>> A=[0 1;-6 -5];syms s;G=inv(s*eye(size(A))-A);phet=ilaplace(G);X0=0;Xt1=phet*X0

Xt1 =

[ 0, 0] [ 0, 0]

>> B=[0 1]\\";Xt2=ilaplace(G*B*(1/s))

Xt2 =

[ 1/6-1/2*exp(-2*t)+1/3*exp(-3*t)] [

exp(-2*t)-exp(-3*t)]

在命令窗中运行下列命令，建立状态空间模型，计算系统在初始状态作用下的状态响应和输出响应，并绘制相应的响应曲线。

>> A=[0 1;-6 -5];B=[0;1];C=[1 0];D=0;G=(A,B,C,D); t=0:0.5:10;x0=[1;0];

[yo,t,xo]=initial(G,x0,t);plot(t,xo,\\":\\",t,yo,\\"-\\") 返回图4。

图4 状态响应

在命令窗中继续运行下列命令，计算系统在输入作用下的状态响应和输出响应，并绘制相应的响应曲线。

>> figure(‘pos’,[50 50 200 150],’color’,’w’); u=ones(size(t)); [yu,t,xu]=lsim(G,u,t); plot(t,xu,’:’,t,yu,’-‘) 返回图5。

图5 输出响应

再继续运行下列命令求系统总的状态响应和输出响应，并绘制相应的响应曲线。

>> y=yo+yu; x=xo+xu; plot(t,x,\\":\\",t,y,\\"-\\") 返回图6。

图6 (3) 在命令窗中运行下列命令

>> A=[0 1;-6 -5];B=[0;1];C=[1 0];D=0;G=(A,B,C,D);t=0:0.5:20;u=exp(-t); [y,t,x]=lsim(G,u,t);plot(t,x,\\":k\\",t,y,\\"-k\\") 可得状态响应和输出响应的数值解以及相应的曲线，如图7。

图7 也可编制如下程序%ex24，先求状态方程的解析解再求数值解，然后绘制曲线。

>> figure(\\"pos\\",[50 50 200 150],\\"color\\",\\"w\\");A=[0 1;-6 -5];B=[0;1];C=[1 0]; syms s;G=inv(s*eye(size(A))-A);phet=ilaplace(G);u=1/s;x=ilaplace(G*B*u);y=C*x; for i=1:61 tt=0.1*(i-1);xt(:,i)=subs(x(:),\\"t\\",tt);yt(i)=subs(y,\\"t\\",tt);

end >> plot(0:60,xt,\\":k\\",0:60,yt,\\"-k\\") >> gtext(\\"y\\",\\"FontSize\\",8) >> gtext(\\"x\\",\\"FontSize\\",8)

在命令窗中运行该程序得到状态和输出响应解析解和数值解，以及相应的曲线如图8。

图8

四．实验总结

1.通过实验，掌握了状态转移矩阵的概念、线性系统状态方程解的结构。

2.学会用MATLAB求解状态转移矩阵、求解线性定常系统的状态响应和输出响应，并绘制相应曲线。

实验三 线性定常系统的能控性和能观测性

一、实验目的

1.掌握能控性和能观测性的概念。学会用MATLAB判断能控性和能观测性。 2.掌握系统的结构分解。学会用MATLAB进行结构分解。 3.掌握最小实现的概念。学会用MATLAB求最小实现。

二 实验原理 1.能控性

1）线性定常系统状态能控性的判断

n阶线性定常连续或离散系统(A,B)状态完全能控的充分必要条件是：能控性矩阵

UcBABA2BAn1B的秩为n。

能控性矩阵可用MATLAB提供的函数ctrb( )自动产生，其调用格式为: Ucctrb(A,B)

其中A,B分别为系统矩阵和输入矩阵，Uc为能控性矩阵。

能控性矩阵的秩即rank(Uc)称为能控性指数，表示系统能控状态变量的数目，可由MATLAB提供的函数rank( )求出。 2）线性定常系统输出能控性的判断

m(n1)r矩阵线性定常连续或离散系统(A,B,C,D)输出能控的充分必要条件是：UyCBCABCA2BCAn1BD的秩为m，其中r为系统的输入个数，m为输出个数。

矩阵Uy可以通过能控性矩阵Uc得到，即UyC*Uc2.能观测性

n阶线性定常连续或离散系统(A,C)状态完全能观测的充分必要条件是：能观测性矩D

CCA2阵VoCA的秩为n。

n1CA能观测性矩阵可以用MATLAB提供的函数obsv( )自动产生，其调用格式为: Voobsv(A,C)

其中A, C分别为系统矩阵和输出矩阵，Vo为能观测性矩阵。

能观测性矩阵的秩即rank(Vo)称为能观测性指数，表示系统能观测状态变量的数目。可由MATLAB提供的函数rank( )求出。 3.最小实现

MATLAB 提供的函数minreal( )可直接得出系统的最小实现，其调用格式为

Gmminreal(G)

其中G为系统的LTI对象，Gm为系统的一个最小实现。

三 实验内容 1.已知系统

344xx1u y11x

10（1）判断系统状态的能控性和能观测性，以及系统输出的能控性。说明状态能

控性和输出能控性之间有无联系。

（3） 将给定的状态空间表达式变换为对角标准型，判断系统的能控性和能观测性，与（1）的结果是否一致？为何？ 解：（1）在命令窗中运行下列命令

>> A=[-3 -4;-1 0]; B=[4;1]; Uc=ctrb(A,B); rank(Uc)

ans =

1 因为rank(Uc)=1n=2，所以系统的状态不完全能控.>> A=[-3 -4;-1 0]; C=[-1 -1]; Vo=obsv(A,C); rank(Vo)

ans =

1 因为rank(Vo)=1n=2，故系统状态不完全能观测

>> A=[-3 -4;-1 0];B=[4;1];C=[-1 -1];D=0;Uc=ctrb(A,B);Uy=[C*Uc D];rank(Uy)

ans =

1 因为rank(Uy)=1=m，故系统是输出能控的。 状态能控性和输出能控性之间没有任何联系。

（3）在命令窗中运行下列命令

>> A=[-3 -4;-1 0];B=[4 1]\\";C=[-1 -1];G=(A,B,C,0); G1=canon(G)

a =

x1 x2

x1 -4

0

x2

0

b =

u1

x1 -4.123

x2

0

c =

x1

x2

y1 1.213

0

d =

u1

y1

0

Continuous-time model.>> A=[-4 0;0 1]; B=[-4.123;0]; Uc=ctrb(A,B); rank(Uc)

ans =

1 因为rank(Uc)=1n=2，所以系统的状态不完全能控.>> A=[-4 0;0 1]; C=[1.213 0]; Vo=obsv(A,C); rank(Vo)

ans =

1 因为rank(Vo)=1n=2，故系统状态不完全能观测。

变换为对角标准型系统的能控性和能观测性与（1）的结果一致，因为变换为对角标准型系统状态矩阵之间秩没变。

3.已知系统

（b）G(s)s1

(s1)(s2)(s3)用函数minreal( )求最小实现。判断所得系统的能控性和能观测性，验证其是否最小实现。 解：在命令窗中运行下列命令

>> z=[-1];p=[-1,-2,-3];k=1;Gzpk=zpk(z,p,k);G=(Gzpk);Gm=minreal(G)

1 state removed.

a =

x1 x2

x1 -2

x2

0 -3

b =

u1

x1

0

x2 0.5

c =

x1

x2

y1 0.5

0

d =

u1

y1

0

Continuous-time model.

>> >> A=[-2 4;0 -3]; B=[0;0.5]; Uc=ctrb(A,B); rank(Uc)

ans =

2 因为rank(Uc)=2= n=2，所以系统的状态完全能控。 >> A=[-2 4;0 -3]; C=[0.5 0]; Vo=obsv(A,C); rank(Vo)

ans =

2 因为rank(Vo)=2=n=2，故系统状态完全能观测。

由于系统既能控又能观，所以系统的实现是最小实现。

四．实验总结

1.通过实验，掌握了能控性和能观测性的概念和最小实现的概念。

2.学会用MATLAB判断能控性和能观测性、用函数minreal( )求最小实现。

字字精辟！

好好感受，好好思考，不要浪费了作者的良苦用心！

推荐访问:理论 实验 控制