1
一、实验目的：
1.熟悉线性系统的数学模型、模型转换。
2.了解MATLAB 中相应的函数。
二、实验内容及步骤：
1.给定系统的传递函数为
18s+36G(s)=3s+40.3s2+391s+150 要求（1）将其用Matlab 表达；
（2）生成状态空间模型。
2.在Matlab 中建立如下离散系统的传递函数模型
y(k+2)+5y(k+1)+6y(k)=u(k+2)+2u(k+1)+u(k) 3.在Matlab 中建立如下传递函数阵的Matlab 模型
⎡s2+2s+1s+5⎤⎢⎥2s+2⎥G(s)=⎢s+5s+62s+36⎥⎢⎢⎣s3+6s2+11s+62s+7⎥⎦ 4.给定系统的模型为
18(s+2)G(s)=(s+15)(s+25)(s+0.4) 求（1）将其用Matlab 表达；
（2）生成状态空间模型。
5.给定系统的状态方程系数矩阵如下：
⎡-40.4-138-160⎤⎡1⎤⎢⎥⎢⎥A=⎢100⎥,B=⎢0⎥,C=[018360],D=0
⎢0⎢0⎥10⎥⎣⎦⎣⎦ 用Matlab 将其以状态空间模型表示出来。
6.输入零极点函数模型，零点z=1,-2;极点p=-1,2,-3 增益k=1；求相应的传递函
数模型、状态空间模型。
三、实验结果及分析：
1.运行程序为：
num = [18 36];
den = [1 40.3 391 150];
tf(num,den)
[A,B,C,D]=tf2ss(num,den)
ss(A,B,C,D)
运行结果为：
Transfer function:
18 s + 36
----------------------------
s^3 + 40.3 s^2 + 391 s + 150
A =
-40.3000 -391.0000 -150.0000
1.0000 0 0
0 1.0000 0
B =
1
C =
0 18 36
D =
a =
x1 x2 x3
x1 -40.3 -391 -150
x2 1 0 0
x3 0 1 0
b =
u1
x1 1
x2 0
x3 0
c =
x1 x2 x3
y1 0 18 36
d =
u1
y1 0
Continuous-time model.
2.运行程序为：
clear all
clc
a1 = 5; a0 = 6; b2 = 1;
b1 = 2;b0 = 1;
c0 = b2;
c1 = b1 - a1*c0;
c2 = b0 - a0*c0 - a1*c1;
G = [0 1;-a0 -a1]; H = [c1;c2];
C = [1 0]; D = [c0];
syms z;
GG = simple(C*inv(z*eye(2)-G)*H)
%求传递函数
运行结果为：
GG =
1/(z + 2) - 4/(z + 3)
3.运行程序为：
num = {[1 2 1];[1 5];[2 3];[6]};
den = {[1 5 6];[1 2];[1 6 11 6];[2 7]};
tf(num,den)
运行结果为：
Transfer function from input to output...
s^2 + 2 s + 1
#1: -------------
s^2 + 5 s + 6
s + 5
#2: -----
s + 2
2 s + 3
#3: ----------------------
s^3 + 6 s^2 + 11 s + 6
6
#4: -------
2 s + 7
4.运行程序为：
num = conv(18,[1 2]);
den = conv([1 15],conv([1 25],[1 0.4]));
tf(num,den)
[A,B,C,D] = tf2ss(num,den)
ss(A,B,C,D)
运行结果：
Transfer function:
18 s + 36
----------------------------
s^3 + 40.4 s^2 + 391 s + 150
A =
-40.4000 -391.0000 -150.0000
1.0000 0 0
0 1.0000 0
B =
1
C =
0 18 36
D =
a =
x1 x2 x3
x1 -40.4 -391 -150
x2 1 0 0
x3 0 1 0
b =
u1
x1 1
x2 0
x3 0
c =
x1 x2 x3
y1 0 18 36
d =
u1
y1 0
Continuous-time model.
5.运行程序为：
A = [-40.4 -138 -160;1 0 0;0 1 0];
B = [1;0;0];
C = [0 18 360];
D = 0;
ss(A,B,C,D)
运行结果为：
a =
x1 x2 x3
x1 -40.4 -138 -160
x2 1 0 0
x3 0 1 0
b =
u1
x1 1
x2 0
x3 0
c =
x1 x2 x3
y1 0 18 360
d =
u1
y1 0
Continuous-time model.
6.运行程序为：
clear all
clc
z = [1;2];p = [-1;2;-3];k = 1;
[num,den] = zp2tf(z,p,k);
[A,B,C,D] = zp2ss(z,p,k);
sys1=tf(num,den)
sys2=ss(A,B,C,D)
运行结果为：
sys1 =
s^2 - 3 s + 2
---------------------
s^3 + 2 s^2 - 5 s - 6
Continuous-time transfer function.
sys2 =
a =
x1 x2 x3
x1 2 0 0
x2 1 -4 -1.732
x3 0 1.732 0
b =
u1
x1 1
x2 0
x3 0
c =
x1 x2 x3
y1 1 -7 -0.5774
d =
u1
y1 0
Co
ntinuous-time state-space model
四、实验总结：
通过对实验一的练习，掌握了系统模型之间的转换，并且了解了Matlab中的一些函数及应用，从而更加深刻的了解了状态空间表达式的涵义。
实验二 状态空间标准形与控制系统的运动分析
一、实验目的：
1.掌握线性系统的对角线标准形、约当标准形、能控标准形和能观测标准型的表
示及相应变换阵的求解。
深入理解状态空间模型的相关理论。
2. 掌握利用Matlab进行矩阵指数函数的数值计算和符号计算方法；对定常连续系
统和定常离散系统的状态空间模型进行求解，分析其运动规律；对连续系统进行离散化。
二、实验原理：
1.规范形转换函数canon() ，可将传递函数模型转换，得到状态空间的能控标准
型和对角规范形，其调用格式为
连续系统：[con_ss,T]=canon (con_tf, ' type ')
离散系统：[dis_ss,T]=canon (dis_tf, ' type ')
说明：con_tf为状态空间模型或传递函数模型；′type′指定变换标准型的类型，有
两个选项：′modal′和′companion′，分别对应对角标准型和能控标准形，T为返回的相似变换矩阵，可缺省。
2.约当标准型
1）特征值、特征向量计算函数eig()，调用格式为：
d = eig(A)
[V,D] = eig(A)
其中，第1种格式只计算所有特征值，输出格式将所有特征值排成向量；第2
种格式可同时得到所有特征向量和特征值，输出格式为所有特征值为对角线元素的对角线矩阵D，所有特征向量为列向量并排成矩阵V。
2）广义特征向量计算函数jordan()，调用格式为：
J = jordan(A)
[V,J] = jordan(A)
其中，第1种调用格式只计算A矩阵对应的约旦矩阵J；第2格式可同时得到所有广
义特征向量和约旦矩阵J，广义特征向量为列向量，并排成矩阵V。
3）一般状态空间模型到约旦规范形的变换
Matlab没有将一般状态空间模型变换成约旦规范形(对角线规范形为其一个特例)
的函数，可利用符号计算工具箱的计算约旦矩阵和广义特征向量的函数jordan()，求解广义特征向量，进而构造变换矩阵求解约旦规范形。
3.非奇异变换
在matlab中用ss2ss来实现动态方程之间的变换
sys_out=ss2ss(sys_in,T)，其中T为非奇异变换阵，其变换关系为z = Tx
z = [TAT-1] z + [TB]u
y = [CT-1] z + Du
4.矩阵eAt的数值计算函数expm()，采用帕德(Pade)逼近法计算矩阵指数eAt,精度高,
数值稳定性好。其调用格式为：
Y = expm(X)
其中,X为输入的需计算矩阵指数的矩阵,Y为计算的结果。
5. Matlab中，线性定常连续系统的状态空间模型求解(即系统运动轨迹的计算)
的主要函数有初始状态响应函数initial()、阶跃响应函数step()以及可计算任意输入的系统响应函数lsim()。
三、实验内容：
1.将实验一的第2题用对角标准型实现。
2.系统的动态方程如下：
⎡0⎡1⎤10⎤∙⎢⎥⎢⎥x=⎢001⎥x+⎢0⎥u,y=[001]x
⎢-6-11-6⎥⎢0⎥⎣⎦⎣⎦
1）求对角标准型实现，并写出实现变换的非奇异阵和变换关系。
2）求可控准型实现，并写出实现变换的非奇异阵和变换关系。
3.计算如下矩阵的特征值和广义特征向量。
⎡-4-3-6⎤⎢⎥A=⎢102⎥
⎢111⎥⎣⎦
4.将如下状态空间模型变换为约旦规范形。
⎡0⎡1⎤10⎤∙⎢⎥⎢⎥x=⎢001⎥x+⎢0⎥u
⎢-6-11-6⎥⎢0⎥⎣⎦⎣⎦
y=[100]x
5.将下列状态方程化为约当型，并写出实现变换的非奇异阵和变换关系。
⎡3-11⎡1100⎤0⎤⎢⎥⎢⎥11-1-100-11⎢⎥⎢⎥⎢00⎢22011⎥1⎥A=⎢,B=⎥⎢⎥ 02-1-1⎥⎢00⎢0-1⎥⎢00⎢00011⎥2⎥⎢⎥⎢⎥⎢⎢1⎥0011⎥0⎦⎣00⎦⎣
6.在Matlab 中计算矩阵A 在t=0.3 时的矩阵指数eAt的值。
⎡01⎤A=⎢⎥
⎣-2-3⎦
7.计算系统在[0,10s]内， T=3s 的单位方波输入的状态响应。
∙⎡0⎡0⎤⎡1⎤1⎤x=⎢x+u,x=0⎥⎢⎥⎢⎥
⎣-2-3⎦⎣1⎦⎣2⎦
四、思考题：
将模型实现为某种标准型（可控标准型，可观测标准型和Jondan标准型）的条件
是什么？
答：将模型实现为某种标准型时，一定要判断模型是否符合某种标准型的要求。 比如可控标准型的要求是 rankQc=B[ABA2B...An-1B=n ； ]
⎡C⎤⎢⎥CA⎥ 可观测标准型的要求是 rankQ0=⎢=n； ⎢⋯⎥⎢n-1⎥⎣CA⎦
Jondan标准型的要求是 A矩阵有相同的特征值。
五、实验报告要求：
1.明确状态空间实现的三种标准型的原理；实验的目的及步骤。
2.扼要记录实验数据，分析实验数据结果。
1.运行程序为：
clear all
clc
a1 = 5; a0 = 6;
b2 = 1; b1 = 2; b0 = 1;
c0 = b2;
c1 = b1 - a1*c0;
c2 = b0 - a0*c0 - a1*c1;
G = [0 1;-a0 -a1];
H = [c1;c2];
C = [1 0];
D = [c0];
[Q,d] = eig(G); %求特征向量和特征值
P = inv(Q); %求变换矩阵
Gb = P*G*Q %求变换之后的各矩阵
Hb = P*H
Cb = C*Q
Db = D
运行结果为：
Gb =
-2.0000 0.0000
-0.0000 -3.0000
Hb =
2.2361
12.6491
Cb =
0.4472 -0.3162
Db =
1
2.(1)运行程序：
clear all
clc
A = [0 1 0;0 0 1;-6 -11 -6];
B = [1;0;0];
C = [1 1 0];
D = 0;
con_ss = canon(A,B,C,D,'modal') %对角标准型的实现
[Q,d] = eig(A);
P = inv(Q); %求变换矩阵P
A1 = P*A*Q %求变换之后的各矩阵
B1 = P*B
C1 = C*Q
D1 = D
运行结果为：
con_ss =
-1.0000 0 0
0 -2.0000 0
0 0 -3.0000
A1 =
-1.0000 -0.0000 0.0000
0.0000 -2.0000 0.0000
0.0000 -0.0000 -3.0000
B1 =
-5.1962
-13.7477
-9.5394
C1 =
0.0000 -0.2182 0.2097
D1 =
2.(2)运行程序为：
clear all
A = [0 1 0;0 0 1;-6 -11 -6];
B = [1;0;0];
C = [1 1 0];
D = 0;
%判断系统是否可控
Qc = ctrb(A,B)
syms s;
sys=det(s*eye(3)-A)
if (rank(Qc)==3)
disp('The system is controllable')
else
disp('The system is uncontrollable')
end
运行结果为：
Qc =
1 0 0
0 0 -6
0 -6 36
sys =
s^3+6*s^2+11*s+6
The system is controllable
% a0 = 6; a1 = 11; a2 = 6;
% 标准能控型的实现，其程序为：
Q = Qc*[11 6 1;6 1 0;1 0 0]
P = inv(Q)
Ab = P*A*Q
Bb = P*B
Cb = C*Q
Db = D
运行结果为：
Q =
11 6 1
-6 0 0
0 -6 0
P =
-0.0000 -0.1667 -0.0000
0 0 -0.1667
1.0000 1.8333 1.0000
Ab =
0.0000 1.0000 0.0000
0.0000 0 1.0000
-6.0000 -11.0000 -6.0000
Bb =
-0.0000
1.0000
Cb =
5 6 1
3. 运行程序为：
clear all
clc
A = [-4 -3 -6;1 0 2;1 1 1];
[V,D] = eig(A) %求特征值及广义特征向量
运行结果为：
V =
-0.9045 -0.9045 -0.8046
0.3015 0.3015 -0.2615
0.3015 0.3015 0.5331
D =
-1.0000 0 0
0 -1.0000 0
0 0 -1.0000
4.运行程序为：
clear all
clc
A = [0 1 0;0 0 1;-6 -11 -6];
B = [1;0;0];
C = [1 0 0];
[Q,J] = jordan(A); %求约旦阵J
P = inv(Q); %求变换矩阵P
Ab = J %求转换之后的各矩阵
Bb = P*B
Cb = C*P
运行结果为：
Ab =
-3 0 0
0 -2 0
0 0 -1
Bb =
9.0000
-12.0000
3.0000
Cb =
9.0000 13.5000 4.5000
5.运行程序为：
clear all
clc
A = [3 -1 1 1 0 0;1 1 -1 -1 0 0;0 0 2 0 1 1;0 0 0 2 -1 -1;0 0 0 0 1 1;0 0 0 0 1 1];
B = [1 0;-1 1;2 1;0 -1;0 2;1 0];
[Q,J] = jordan(A) %求约旦阵J
P = inv(Q) %求变换矩阵P
Ab = J %求转换之后的各矩阵
Bb = P*B
运行结果为：
Q =
2.0000 2.0000 1.0000 0 0 0 2.0000 0 0 0 0 0 0 0 1.0000 0 1.0000 0 0 0 0 0 -1.0000 0 0 0 0 0.5000 0 0.5000 0 0 0 -0.5000 0 0.5000 J =
2 1 0 0 0 0 0 2 1 0 0 0 0 0 2 0 0 0 0 0 0 0 0 0 0 0 0 0 2 1 0 0 0 0 0 2 P =
0 0.5000 0 0 0 0.5000 -0.5000 -0.5000 -0.5000 0 0 0 1.0000 1.0000 0 0 0 0 0 1.0000 0 0 0 -1.0000 0 0 0 0 0 1.0000 Ab =
2 1 0 0 0 0 0 2 1 0 0 0 0 0 2 0 0 0 0 0 0 0 0 0 0 0 0 0 2 1 0 0 0 0 0 2 Bb =
-0.5000 0.5000 0 -0.5000 2.0000 0 -1.0000 2.0000 0 1.0000 1.0000 2.0000 6.运行程序为： clear all clc
A = [0 1;-2 -3]; syms s t;
F = s*eye(length(A))-A; Q = inv(F); f = ilaplace(Q); x = subs(f,'t',0.3) 运行结果为： x =
[ 2*exp(-3/10) - exp(-3/5), exp(-3/10) - exp(-3/5)] [ 2*exp(-3/5) - 2*exp(-3/10), 2*exp(-3/5) - exp(-3/10)]
简便的方法：
0 0 0 -1.0000 0 1.0000
clear all clc
A = [0 1;-2 -3]; t = 0.3;
y = expm(A*t) 运行结果为： y =
0.9328 0.1920 -0.3840 0.3568
7.运行程序为： clear all clc
A = [0 1;-2 -3]; B = [0;1]; C = []; D = []; x0 = [1;2];
sys = ss(A,B,C,D);
[u,t] = gensig('square',3,10,0.1); %产生周期为3s，时间为10s，采样周期为0.1s的方波信号
[y,t,x] = lsim(sys,u,t,x0); plot(t,u,t,x) 运行结果为：
实验三系统的能控性和能观性
一、实验目的:
系统的能控性和能观性关系到系统的极点配置法设计，和最优控制。通过本实验，掌握判断系统能控性能观性的条件和方法。 二、实验原理：
Matlab 用于状态能控性、能观性判定方法如下：
能控性矩阵函数ctrb()、能观性矩阵函数obsv()和能控性/能观性格拉姆矩阵函数gram()，通过对这些函数计算所得的矩阵求秩就可以很方便地判定系统的状态能控性、能观性。
1)检验系统的能控性和能观性。
⎡11⎤⎡1-1⎤⎡10⎤A=⎢,B=,C=⎥⎢1-1⎥⎢01⎥,D=0 4-2⎣⎦⎣⎦⎣⎦2)判定离散系统的状态能控性。
⎡100⎤⎡1⎤
⎥x(k)+⎢2⎥u(k) x(k+1)=⎢02-2⎢⎥⎢⎥
⎢⎢⎣-110⎥⎦⎣1⎥⎦
3)用格拉姆矩阵判据判断下面系统的可控性和可观测性。 ⎡-1-2-3⎤⎡2⎤∙
⎥x+⎢0⎥ux=⎢0-11⎢⎥⎢⎥
⎢⎢⎣10-1⎥⎦⎣1⎥⎦
y=[110]x
4)将下面系统化为可控标准型和可观测标准型。
⎡120⎤⎡2⎤
⎥x+⎢1⎥ux=⎢3-11⎢⎥⎢⎥
⎢⎢⎣020⎥⎦⎣1⎥⎦y=[001]x结果如下： 1.
A = [1 1;4 -2]; B = [1 -1;1 -1]; C = [1 0;0 1]; D = 0;
Qc = ctrb(A,B) Qo = obsv(A,C)
if rank(Qc) == length(A)
disp('The system is controllable') else
disp('The system is uncotrollable') end
if rank(Qo) == length(A)
disp('The system is observable') else
disp('The system is not observable') end
运行结果为： Qc =
1 -1 2 -2 1 -1 2 -2 Qo =
1 0 0 1 1 1 4 -2
The system is uncotrollable The system is observable 2.
G = [1 0 0;0 2 -2; -1 1 0]; H = [1;2;1]; Qc = ctrb(G,H)
if rank(Qc) == length(A)
∙
disp('The system is controllable') else
disp('The system is uncotrollable') end
运行结果为： Qc =
1 1 1 2 2 2 1 1 1
The system is uncontrollable
3. A = [-1 -2 -3;0 -1 1;1 0 -1]; B = [2;0;1]; C = [1 1 0]; syms t s;
F = inv(s*eye(length(A))-A) f = ilaplace(F)
4. %求能控标准型 clear all clc
A = [1 2 0;3 -1 1;0 2 0]; B = [2;1;1]; C = [0 0 1]; D = 0;
Qc = ctrb(A,B) %判断其系统是否可控 syms s;
sys=det(s*eye(3)-A); if (rank(Qc)==3)
disp('The system is controllable') else
disp('The system is uncontrollable') end
运行结果为：
sys = s^3-9*s+2
The system is controllable
%由上面可知：a0 = 2; a1 = -9; a2 = 0; a3 = 1; Q = Qc*[-9 0 1;0 1 0;1 0 0]; %能控标准型的实现 P = inv(Q) Ab = P*A*Q Bb = P*B Cb = C*Q Db = D
P =
-0.1250 0 0.2500 -0.1250 0.2500 -0.0000 0.6250 -0.5000 0.2500 Ab =
0 1.0000 0 0.0000 -0.0000 1.0000 -2.0000 9.0000 0.0000 Bb =
0 -0.0000 1.0000 Cb =
3 2 1 Db =
%求能观标准型 clear all clc
A = [1 2 0;3 -1 1;0 2 0]; B = [2;1;1]; C = [0 0 1]; D = 0;
Qo = obsv(A,C) %判断其系统是否可观 syms s;
sys=det(s*eye(3)-A) if (rank(Qo)==3)
disp('The system is observable') else
disp('The system is not observable') end
sys = s^3-9*s+2
The system is observable
%由上面可知：a0 = 2; a1 = -9; a2 = 0; a3 = 1; P = [-9 0 1;0 1 0;1 0 0]*Qo %能观标准型的实现 Q = inv(P); Ab = P*A*Q Bb = P*B Cb = C*Q Db = D
运行结果为： P =
6 -2 -7 0 2 0 0 0 1 Ab =
0 0 -2.0000 1.0000 -0.0000 9.0000 0 1.0000 0 Bb =
3 2 1 Cb =
0 0 1 Db =
四、思考题：
控制系统的能控性和能观性的物理意义是什么？ 能控性是描述状态变量与输入量的关系，而能观测性是描述输出量与状态变量之间 的关系。
五、实验报告要求：
1.掌握判断系统能控性和能观性的原理和方法；实验的目的及步骤。 2.精确记录实验数据，分析实验数据结果。
实验四李雅普诺夫(Lyapunov)判据实验
一、实验目的：
现代控制系统分析和设计中，很多问题归结于李雅普诺夫 (Lyapunov)方程求解，如稳定性的研究等。通过实验掌握求解方法；判断系统的稳定性。 二、实验原理： 1）在Matlab 控制箱中，函数LYAP 和DLYAP 用于求解Lyapunov 方程。其中LYAP 用于求解连续系统的Lyapunov 方程，DLYAP 用于求解离散系统的Lyapuno 方程。调用格式：
p=lyap(A,Q) 该函数可以求解如ATP+PA=-Q的Lyapunov 方程。其中Q 和A（系统矩阵）是具有相同维数的方阵，如Q 对称，则P 也是对称的。
2）线性系统的稳定性的充分必要条件：微分方程的全部根都是负实数或实部为负的复数，亦即全部根都位于左半复平面。 三、实验内容：
1.在Matlab 中判定如下系统的李雅普诺夫稳定性，并通过各阶主子式来判断定号性。 ⎡∙⎤⎡01⎤⎡x1⎤1⎥⎢x=⎢∙⎥⎢⎥ ⎢x⎥⎣-1-1⎦⎣x2⎦⎣2⎦
⎡000⎤
⎥，求Lyapunov 方程0002.已知系统的方框图如下， 选择正半定对称矩阵Q=⎢⎢⎥
⎢⎣001⎥⎦的解，试在Matlab 中判定系统的稳定性。
⎡100⎤
⎥，求稳定性判别矩阵P，判断0103.系统状态空间模型如下，选取正定矩阵Q=⎢⎢⎥
⎢⎣001⎥⎦
该系统的稳定性，并求该系统的传递函数矩阵。
⎡-3-8-2-4⎤⎡1⎤⎢1⎥⎢0⎥000⎥,B=⎢⎥,C=[0011],D=0 A=⎢⎢0⎢0⎥100⎥⎢⎥⎢⎥
010⎦⎣0⎣0⎦
实验结果如下：
1.
clear all clc
A = [0 1;-1 -1]; Q = eye(length(A)); P = lyap(A,Q);
if ((P(1,1)>0)&(det(P)>0)) disp('系统稳定') else
disp('无法判断系统的稳定性') end
运行结果为： 系统稳定 2.
clear all clc
num1 = 5;
den1 = conv([1 0],conv([1 1],[1 2])); [num,den] = cloop(num1,den1,-1); [A,B,C,D] = tf2ss(num,den); Q = [0 0 0;0 0 0;0 0 1]; P = lyap(A,Q);
if((P(1,1)>0)&(P(1,1)*P(2,2)-P(1,2)*P(2,1)>0)&(det(P)>0)) disp('系统稳定') else
disp('无法判断系统的稳定性') end
运行结果为： 系统稳定
3. clear all clc
num1 = 5;
den1 = conv([1 0],conv([1 1],[1 2])); [num,den] = cloop(num1,den1,-1); [A,B,C,D] = tf2ss(num,den); Q = [0 0 0;0 0 0;0 0 1]; P = lyap(A,Q); for i=1:length(P)
H = det(P(1:i,1:i)) end
四、思考题：
一个系统如果没有找到Lyapunov 函数，是否能对系统的稳定性做出否定的结论？ 答 ：不能，Lyapunov 函数和系统的稳定性之间不是充分必要条件，Lyapunov 函数是系统稳定的充分条件. 五、实验报告要求：
1.深刻理解系统稳定性的充要条件；Lyapunov 第二方法判断系统稳定性的条件；明
确实验的目的及步骤。
2.如实记录实验数据，分析实验数据结果。
21
1
一、实验目的：
1.熟悉线性系统的数学模型、模型转换。
2.了解MATLAB 中相应的函数。
二、实验内容及步骤：
1.给定系统的传递函数为
18s+36G(s)=3s+40.3s2+391s+150 要求（1）将其用Matlab 表达；
（2）生成状态空间模型。
2.在Matlab 中建立如下离散系统的传递函数模型
y(k+2)+5y(k+1)+6y(k)=u(k+2)+2u(k+1)+u(k) 3.在Matlab 中建立如下传递函数阵的Matlab 模型
⎡s2+2s+1s+5⎤⎢⎥2s+2⎥G(s)=⎢s+5s+62s+36⎥⎢⎢⎣s3+6s2+11s+62s+7⎥⎦ 4.给定系统的模型为
18(s+2)G(s)=(s+15)(s+25)(s+0.4) 求（1）将其用Matlab 表达；
（2）生成状态空间模型。
5.给定系统的状态方程系数矩阵如下：
⎡-40.4-138-160⎤⎡1⎤⎢⎥⎢⎥A=⎢100⎥,B=⎢0⎥,C=[018360],D=0
⎢0⎢0⎥10⎥⎣⎦⎣⎦ 用Matlab 将其以状态空间模型表示出来。
6.输入零极点函数模型，零点z=1,-2;极点p=-1,2,-3 增益k=1；求相应的传递函
数模型、状态空间模型。
三、实验结果及分析：
1.运行程序为：
num = [18 36];
den = [1 40.3 391 150];
tf(num,den)
[A,B,C,D]=tf2ss(num,den)
ss(A,B,C,D)
运行结果为：
Transfer function:
18 s + 36
----------------------------
s^3 + 40.3 s^2 + 391 s + 150
A =
-40.3000 -391.0000 -150.0000
1.0000 0 0
0 1.0000 0
B =
1
C =
0 18 36
D =
a =
x1 x2 x3
x1 -40.3 -391 -150
x2 1 0 0
x3 0 1 0
b =
u1
x1 1
x2 0
x3 0
c =
x1 x2 x3
y1 0 18 36
d =
u1
y1 0
Continuous-time model.
2.运行程序为：
clear all
clc
a1 = 5; a0 = 6; b2 = 1;
b1 = 2;b0 = 1;
c0 = b2;
c1 = b1 - a1*c0;
c2 = b0 - a0*c0 - a1*c1;
G = [0 1;-a0 -a1]; H = [c1;c2];
C = [1 0]; D = [c0];
syms z;
GG = simple(C*inv(z*eye(2)-G)*H)
%求传递函数
运行结果为：
GG =
1/(z + 2) - 4/(z + 3)
3.运行程序为：
num = {[1 2 1];[1 5];[2 3];[6]};
den = {[1 5 6];[1 2];[1 6 11 6];[2 7]};
tf(num,den)
运行结果为：
Transfer function from input to output...
s^2 + 2 s + 1
#1: -------------
s^2 + 5 s + 6
s + 5
#2: -----
s + 2
2 s + 3
#3: ----------------------
s^3 + 6 s^2 + 11 s + 6
6
#4: -------
2 s + 7
4.运行程序为：
num = conv(18,[1 2]);
den = conv([1 15],conv([1 25],[1 0.4]));
tf(num,den)
[A,B,C,D] = tf2ss(num,den)
ss(A,B,C,D)
运行结果：
Transfer function:
18 s + 36
----------------------------
s^3 + 40.4 s^2 + 391 s + 150
A =
-40.4000 -391.0000 -150.0000
1.0000 0 0
0 1.0000 0
B =
1
C =
0 18 36
D =
a =
x1 x2 x3
x1 -40.4 -391 -150
x2 1 0 0
x3 0 1 0
b =
u1
x1 1
x2 0
x3 0
c =
x1 x2 x3
y1 0 18 36
d =
u1
y1 0
Continuous-time model.
5.运行程序为：
A = [-40.4 -138 -160;1 0 0;0 1 0];
B = [1;0;0];
C = [0 18 360];
D = 0;
ss(A,B,C,D)
运行结果为：
a =
x1 x2 x3
x1 -40.4 -138 -160
x2 1 0 0
x3 0 1 0
b =
u1
x1 1
x2 0
x3 0
c =
x1 x2 x3
y1 0 18 360
d =
u1
y1 0
Continuous-time model.
6.运行程序为：
clear all
clc
z = [1;2];p = [-1;2;-3];k = 1;
[num,den] = zp2tf(z,p,k);
[A,B,C,D] = zp2ss(z,p,k);
sys1=tf(num,den)
sys2=ss(A,B,C,D)
运行结果为：
sys1 =
s^2 - 3 s + 2
---------------------
s^3 + 2 s^2 - 5 s - 6
Continuous-time transfer function.
sys2 =
a =
x1 x2 x3
x1 2 0 0
x2 1 -4 -1.732
x3 0 1.732 0
b =
u1
x1 1
x2 0
x3 0
c =
x1 x2 x3
y1 1 -7 -0.5774
d =
u1
y1 0
Co
ntinuous-time state-space model
四、实验总结：
通过对实验一的练习，掌握了系统模型之间的转换，并且了解了Matlab中的一些函数及应用，从而更加深刻的了解了状态空间表达式的涵义。
实验二 状态空间标准形与控制系统的运动分析
一、实验目的：
1.掌握线性系统的对角线标准形、约当标准形、能控标准形和能观测标准型的表
示及相应变换阵的求解。
深入理解状态空间模型的相关理论。
2. 掌握利用Matlab进行矩阵指数函数的数值计算和符号计算方法；对定常连续系
统和定常离散系统的状态空间模型进行求解，分析其运动规律；对连续系统进行离散化。
二、实验原理：
1.规范形转换函数canon() ，可将传递函数模型转换，得到状态空间的能控标准
型和对角规范形，其调用格式为
连续系统：[con_ss,T]=canon (con_tf, ' type ')
离散系统：[dis_ss,T]=canon (dis_tf, ' type ')
说明：con_tf为状态空间模型或传递函数模型；′type′指定变换标准型的类型，有
两个选项：′modal′和′companion′，分别对应对角标准型和能控标准形，T为返回的相似变换矩阵，可缺省。
2.约当标准型
1）特征值、特征向量计算函数eig()，调用格式为：
d = eig(A)
[V,D] = eig(A)
其中，第1种格式只计算所有特征值，输出格式将所有特征值排成向量；第2
种格式可同时得到所有特征向量和特征值，输出格式为所有特征值为对角线元素的对角线矩阵D，所有特征向量为列向量并排成矩阵V。
2）广义特征向量计算函数jordan()，调用格式为：
J = jordan(A)
[V,J] = jordan(A)
其中，第1种调用格式只计算A矩阵对应的约旦矩阵J；第2格式可同时得到所有广
义特征向量和约旦矩阵J，广义特征向量为列向量，并排成矩阵V。
3）一般状态空间模型到约旦规范形的变换
Matlab没有将一般状态空间模型变换成约旦规范形(对角线规范形为其一个特例)
的函数，可利用符号计算工具箱的计算约旦矩阵和广义特征向量的函数jordan()，求解广义特征向量，进而构造变换矩阵求解约旦规范形。
3.非奇异变换
在matlab中用ss2ss来实现动态方程之间的变换
sys_out=ss2ss(sys_in,T)，其中T为非奇异变换阵，其变换关系为z = Tx
z = [TAT-1] z + [TB]u
y = [CT-1] z + Du
4.矩阵eAt的数值计算函数expm()，采用帕德(Pade)逼近法计算矩阵指数eAt,精度高,
数值稳定性好。其调用格式为：
Y = expm(X)
其中,X为输入的需计算矩阵指数的矩阵,Y为计算的结果。
5. Matlab中，线性定常连续系统的状态空间模型求解(即系统运动轨迹的计算)
的主要函数有初始状态响应函数initial()、阶跃响应函数step()以及可计算任意输入的系统响应函数lsim()。
三、实验内容：
1.将实验一的第2题用对角标准型实现。
2.系统的动态方程如下：
⎡0⎡1⎤10⎤∙⎢⎥⎢⎥x=⎢001⎥x+⎢0⎥u,y=[001]x
⎢-6-11-6⎥⎢0⎥⎣⎦⎣⎦
1）求对角标准型实现，并写出实现变换的非奇异阵和变换关系。
2）求可控准型实现，并写出实现变换的非奇异阵和变换关系。
3.计算如下矩阵的特征值和广义特征向量。
⎡-4-3-6⎤⎢⎥A=⎢102⎥
⎢111⎥⎣⎦
4.将如下状态空间模型变换为约旦规范形。
⎡0⎡1⎤10⎤∙⎢⎥⎢⎥x=⎢001⎥x+⎢0⎥u
⎢-6-11-6⎥⎢0⎥⎣⎦⎣⎦
y=[100]x
5.将下列状态方程化为约当型，并写出实现变换的非奇异阵和变换关系。
⎡3-11⎡1100⎤0⎤⎢⎥⎢⎥11-1-100-11⎢⎥⎢⎥⎢00⎢22011⎥1⎥A=⎢,B=⎥⎢⎥ 02-1-1⎥⎢00⎢0-1⎥⎢00⎢00011⎥2⎥⎢⎥⎢⎥⎢⎢1⎥0011⎥0⎦⎣00⎦⎣
6.在Matlab 中计算矩阵A 在t=0.3 时的矩阵指数eAt的值。
⎡01⎤A=⎢⎥
⎣-2-3⎦
7.计算系统在[0,10s]内， T=3s 的单位方波输入的状态响应。
∙⎡0⎡0⎤⎡1⎤1⎤x=⎢x+u,x=0⎥⎢⎥⎢⎥
⎣-2-3⎦⎣1⎦⎣2⎦
四、思考题：
将模型实现为某种标准型（可控标准型，可观测标准型和Jondan标准型）的条件
是什么？
答：将模型实现为某种标准型时，一定要判断模型是否符合某种标准型的要求。 比如可控标准型的要求是 rankQc=B[ABA2B...An-1B=n ； ]
⎡C⎤⎢⎥CA⎥ 可观测标准型的要求是 rankQ0=⎢=n； ⎢⋯⎥⎢n-1⎥⎣CA⎦
Jondan标准型的要求是 A矩阵有相同的特征值。
五、实验报告要求：
1.明确状态空间实现的三种标准型的原理；实验的目的及步骤。
2.扼要记录实验数据，分析实验数据结果。
1.运行程序为：
clear all
clc
a1 = 5; a0 = 6;
b2 = 1; b1 = 2; b0 = 1;
c0 = b2;
c1 = b1 - a1*c0;
c2 = b0 - a0*c0 - a1*c1;
G = [0 1;-a0 -a1];
H = [c1;c2];
C = [1 0];
D = [c0];
[Q,d] = eig(G); %求特征向量和特征值
P = inv(Q); %求变换矩阵
Gb = P*G*Q %求变换之后的各矩阵
Hb = P*H
Cb = C*Q
Db = D
运行结果为：
Gb =
-2.0000 0.0000
-0.0000 -3.0000
Hb =
2.2361
12.6491
Cb =
0.4472 -0.3162
Db =
1
2.(1)运行程序：
clear all
clc
A = [0 1 0;0 0 1;-6 -11 -6];
B = [1;0;0];
C = [1 1 0];
D = 0;
con_ss = canon(A,B,C,D,'modal') %对角标准型的实现
[Q,d] = eig(A);
P = inv(Q); %求变换矩阵P
A1 = P*A*Q %求变换之后的各矩阵
B1 = P*B
C1 = C*Q
D1 = D
运行结果为：
con_ss =
-1.0000 0 0
0 -2.0000 0
0 0 -3.0000
A1 =
-1.0000 -0.0000 0.0000
0.0000 -2.0000 0.0000
0.0000 -0.0000 -3.0000
B1 =
-5.1962
-13.7477
-9.5394
C1 =
0.0000 -0.2182 0.2097
D1 =
2.(2)运行程序为：
clear all
A = [0 1 0;0 0 1;-6 -11 -6];
B = [1;0;0];
C = [1 1 0];
D = 0;
%判断系统是否可控
Qc = ctrb(A,B)
syms s;
sys=det(s*eye(3)-A)
if (rank(Qc)==3)
disp('The system is controllable')
else
disp('The system is uncontrollable')
end
运行结果为：
Qc =
1 0 0
0 0 -6
0 -6 36
sys =
s^3+6*s^2+11*s+6
The system is controllable
% a0 = 6; a1 = 11; a2 = 6;
% 标准能控型的实现，其程序为：
Q = Qc*[11 6 1;6 1 0;1 0 0]
P = inv(Q)
Ab = P*A*Q
Bb = P*B
Cb = C*Q
Db = D
运行结果为：
Q =
11 6 1
-6 0 0
0 -6 0
P =
-0.0000 -0.1667 -0.0000
0 0 -0.1667
1.0000 1.8333 1.0000
Ab =
0.0000 1.0000 0.0000
0.0000 0 1.0000
-6.0000 -11.0000 -6.0000
Bb =
-0.0000
1.0000
Cb =
5 6 1
3. 运行程序为：
clear all
clc
A = [-4 -3 -6;1 0 2;1 1 1];
[V,D] = eig(A) %求特征值及广义特征向量
运行结果为：
V =
-0.9045 -0.9045 -0.8046
0.3015 0.3015 -0.2615
0.3015 0.3015 0.5331
D =
-1.0000 0 0
0 -1.0000 0
0 0 -1.0000
4.运行程序为：
clear all
clc
A = [0 1 0;0 0 1;-6 -11 -6];
B = [1;0;0];
C = [1 0 0];
[Q,J] = jordan(A); %求约旦阵J
P = inv(Q); %求变换矩阵P
Ab = J %求转换之后的各矩阵
Bb = P*B
Cb = C*P
运行结果为：
Ab =
-3 0 0
0 -2 0
0 0 -1
Bb =
9.0000
-12.0000
3.0000
Cb =
9.0000 13.5000 4.5000
5.运行程序为：
clear all
clc
A = [3 -1 1 1 0 0;1 1 -1 -1 0 0;0 0 2 0 1 1;0 0 0 2 -1 -1;0 0 0 0 1 1;0 0 0 0 1 1];
B = [1 0;-1 1;2 1;0 -1;0 2;1 0];
[Q,J] = jordan(A) %求约旦阵J
P = inv(Q) %求变换矩阵P
Ab = J %求转换之后的各矩阵
Bb = P*B
运行结果为：
Q =
2.0000 2.0000 1.0000 0 0 0 2.0000 0 0 0 0 0 0 0 1.0000 0 1.0000 0 0 0 0 0 -1.0000 0 0 0 0 0.5000 0 0.5000 0 0 0 -0.5000 0 0.5000 J =
2 1 0 0 0 0 0 2 1 0 0 0 0 0 2 0 0 0 0 0 0 0 0 0 0 0 0 0 2 1 0 0 0 0 0 2 P =
0 0.5000 0 0 0 0.5000 -0.5000 -0.5000 -0.5000 0 0 0 1.0000 1.0000 0 0 0 0 0 1.0000 0 0 0 -1.0000 0 0 0 0 0 1.0000 Ab =
2 1 0 0 0 0 0 2 1 0 0 0 0 0 2 0 0 0 0 0 0 0 0 0 0 0 0 0 2 1 0 0 0 0 0 2 Bb =
-0.5000 0.5000 0 -0.5000 2.0000 0 -1.0000 2.0000 0 1.0000 1.0000 2.0000 6.运行程序为： clear all clc
A = [0 1;-2 -3]; syms s t;
F = s*eye(length(A))-A; Q = inv(F); f = ilaplace(Q); x = subs(f,'t',0.3) 运行结果为： x =
[ 2*exp(-3/10) - exp(-3/5), exp(-3/10) - exp(-3/5)] [ 2*exp(-3/5) - 2*exp(-3/10), 2*exp(-3/5) - exp(-3/10)]
简便的方法：
0 0 0 -1.0000 0 1.0000
clear all clc
A = [0 1;-2 -3]; t = 0.3;
y = expm(A*t) 运行结果为： y =
0.9328 0.1920 -0.3840 0.3568
7.运行程序为： clear all clc
A = [0 1;-2 -3]; B = [0;1]; C = []; D = []; x0 = [1;2];
sys = ss(A,B,C,D);
[u,t] = gensig('square',3,10,0.1); %产生周期为3s，时间为10s，采样周期为0.1s的方波信号
[y,t,x] = lsim(sys,u,t,x0); plot(t,u,t,x) 运行结果为：
实验三系统的能控性和能观性
一、实验目的:
系统的能控性和能观性关系到系统的极点配置法设计，和最优控制。通过本实验，掌握判断系统能控性能观性的条件和方法。 二、实验原理：
Matlab 用于状态能控性、能观性判定方法如下：
能控性矩阵函数ctrb()、能观性矩阵函数obsv()和能控性/能观性格拉姆矩阵函数gram()，通过对这些函数计算所得的矩阵求秩就可以很方便地判定系统的状态能控性、能观性。
1)检验系统的能控性和能观性。
⎡11⎤⎡1-1⎤⎡10⎤A=⎢,B=,C=⎥⎢1-1⎥⎢01⎥,D=0 4-2⎣⎦⎣⎦⎣⎦2)判定离散系统的状态能控性。
⎡100⎤⎡1⎤
⎥x(k)+⎢2⎥u(k) x(k+1)=⎢02-2⎢⎥⎢⎥
⎢⎢⎣-110⎥⎦⎣1⎥⎦
3)用格拉姆矩阵判据判断下面系统的可控性和可观测性。 ⎡-1-2-3⎤⎡2⎤∙
⎥x+⎢0⎥ux=⎢0-11⎢⎥⎢⎥
⎢⎢⎣10-1⎥⎦⎣1⎥⎦
y=[110]x
4)将下面系统化为可控标准型和可观测标准型。
⎡120⎤⎡2⎤
⎥x+⎢1⎥ux=⎢3-11⎢⎥⎢⎥
⎢⎢⎣020⎥⎦⎣1⎥⎦y=[001]x结果如下： 1.
A = [1 1;4 -2]; B = [1 -1;1 -1]; C = [1 0;0 1]; D = 0;
Qc = ctrb(A,B) Qo = obsv(A,C)
if rank(Qc) == length(A)
disp('The system is controllable') else
disp('The system is uncotrollable') end
if rank(Qo) == length(A)
disp('The system is observable') else
disp('The system is not observable') end
运行结果为： Qc =
1 -1 2 -2 1 -1 2 -2 Qo =
1 0 0 1 1 1 4 -2
The system is uncotrollable The system is observable 2.
G = [1 0 0;0 2 -2; -1 1 0]; H = [1;2;1]; Qc = ctrb(G,H)
if rank(Qc) == length(A)
∙
disp('The system is controllable') else
disp('The system is uncotrollable') end
运行结果为： Qc =
1 1 1 2 2 2 1 1 1
The system is uncontrollable
3. A = [-1 -2 -3;0 -1 1;1 0 -1]; B = [2;0;1]; C = [1 1 0]; syms t s;
F = inv(s*eye(length(A))-A) f = ilaplace(F)
4. %求能控标准型 clear all clc
A = [1 2 0;3 -1 1;0 2 0]; B = [2;1;1]; C = [0 0 1]; D = 0;
Qc = ctrb(A,B) %判断其系统是否可控 syms s;
sys=det(s*eye(3)-A); if (rank(Qc)==3)
disp('The system is controllable') else
disp('The system is uncontrollable') end
运行结果为：
sys = s^3-9*s+2
The system is controllable
%由上面可知：a0 = 2; a1 = -9; a2 = 0; a3 = 1; Q = Qc*[-9 0 1;0 1 0;1 0 0]; %能控标准型的实现 P = inv(Q) Ab = P*A*Q Bb = P*B Cb = C*Q Db = D
P =
-0.1250 0 0.2500 -0.1250 0.2500 -0.0000 0.6250 -0.5000 0.2500 Ab =
0 1.0000 0 0.0000 -0.0000 1.0000 -2.0000 9.0000 0.0000 Bb =
0 -0.0000 1.0000 Cb =
3 2 1 Db =
%求能观标准型 clear all clc
A = [1 2 0;3 -1 1;0 2 0]; B = [2;1;1]; C = [0 0 1]; D = 0;
Qo = obsv(A,C) %判断其系统是否可观 syms s;
sys=det(s*eye(3)-A) if (rank(Qo)==3)
disp('The system is observable') else
disp('The system is not observable') end
sys = s^3-9*s+2
The system is observable
%由上面可知：a0 = 2; a1 = -9; a2 = 0; a3 = 1; P = [-9 0 1;0 1 0;1 0 0]*Qo %能观标准型的实现 Q = inv(P); Ab = P*A*Q Bb = P*B Cb = C*Q Db = D
运行结果为： P =
6 -2 -7 0 2 0 0 0 1 Ab =
0 0 -2.0000 1.0000 -0.0000 9.0000 0 1.0000 0 Bb =
3 2 1 Cb =
0 0 1 Db =
四、思考题：
控制系统的能控性和能观性的物理意义是什么？ 能控性是描述状态变量与输入量的关系，而能观测性是描述输出量与状态变量之间 的关系。
五、实验报告要求：
1.掌握判断系统能控性和能观性的原理和方法；实验的目的及步骤。 2.精确记录实验数据，分析实验数据结果。
实验四李雅普诺夫(Lyapunov)判据实验
一、实验目的：
现代控制系统分析和设计中，很多问题归结于李雅普诺夫 (Lyapunov)方程求解，如稳定性的研究等。通过实验掌握求解方法；判断系统的稳定性。 二、实验原理： 1）在Matlab 控制箱中，函数LYAP 和DLYAP 用于求解Lyapunov 方程。其中LYAP 用于求解连续系统的Lyapunov 方程，DLYAP 用于求解离散系统的Lyapuno 方程。调用格式：
p=lyap(A,Q) 该函数可以求解如ATP+PA=-Q的Lyapunov 方程。其中Q 和A（系统矩阵）是具有相同维数的方阵，如Q 对称，则P 也是对称的。
2）线性系统的稳定性的充分必要条件：微分方程的全部根都是负实数或实部为负的复数，亦即全部根都位于左半复平面。 三、实验内容：
1.在Matlab 中判定如下系统的李雅普诺夫稳定性，并通过各阶主子式来判断定号性。 ⎡∙⎤⎡01⎤⎡x1⎤1⎥⎢x=⎢∙⎥⎢⎥ ⎢x⎥⎣-1-1⎦⎣x2⎦⎣2⎦
⎡000⎤
⎥，求Lyapunov 方程0002.已知系统的方框图如下， 选择正半定对称矩阵Q=⎢⎢⎥
⎢⎣001⎥⎦的解，试在Matlab 中判定系统的稳定性。
⎡100⎤
⎥，求稳定性判别矩阵P，判断0103.系统状态空间模型如下，选取正定矩阵Q=⎢⎢⎥
⎢⎣001⎥⎦
该系统的稳定性，并求该系统的传递函数矩阵。
⎡-3-8-2-4⎤⎡1⎤⎢1⎥⎢0⎥000⎥,B=⎢⎥,C=[0011],D=0 A=⎢⎢0⎢0⎥100⎥⎢⎥⎢⎥
010⎦⎣0⎣0⎦
实验结果如下：
1.
clear all clc
A = [0 1;-1 -1]; Q = eye(length(A)); P = lyap(A,Q);
if ((P(1,1)>0)&(det(P)>0)) disp('系统稳定') else
disp('无法判断系统的稳定性') end
运行结果为： 系统稳定 2.
clear all clc
num1 = 5;
den1 = conv([1 0],conv([1 1],[1 2])); [num,den] = cloop(num1,den1,-1); [A,B,C,D] = tf2ss(num,den); Q = [0 0 0;0 0 0;0 0 1]; P = lyap(A,Q);
if((P(1,1)>0)&(P(1,1)*P(2,2)-P(1,2)*P(2,1)>0)&(det(P)>0)) disp('系统稳定') else
disp('无法判断系统的稳定性') end
运行结果为： 系统稳定
3. clear all clc
num1 = 5;
den1 = conv([1 0],conv([1 1],[1 2])); [num,den] = cloop(num1,den1,-1); [A,B,C,D] = tf2ss(num,den); Q = [0 0 0;0 0 0;0 0 1]; P = lyap(A,Q); for i=1:length(P)
H = det(P(1:i,1:i)) end
四、思考题：
一个系统如果没有找到Lyapunov 函数，是否能对系统的稳定性做出否定的结论？ 答 ：不能，Lyapunov 函数和系统的稳定性之间不是充分必要条件，Lyapunov 函数是系统稳定的充分条件. 五、实验报告要求：
1.深刻理解系统稳定性的充要条件；Lyapunov 第二方法判断系统稳定性的条件；明
确实验的目的及步骤。
2.如实记录实验数据，分析实验数据结果。
21