系统建模与仿真习题2及答案

系统建模与仿真习题二及答案
1. 考虑如图所示的典型反馈控制系统框图
(1)假设各个子传递函数模型为
66.031.05
.02)(232++-+=s s s s s G ，s s s G c 610)(+=，2
1)(+=s s H 分别用feedback （）函数以及G*Gc/(1+G*Gc*H)（要最小实现）方法求该系统的传递函数模型。

(2) 假设系统的受控对象模型为s e s s s G 23
)1(12
)(-+=，控制器模型为 s
s s G c 32)(+=，并假设系统是单位负反馈，分别用feedback （）函数以及G*Gc/(1+G*Gc*H)（要最小实现）方法能求出该系统的传递函数模型？如果不能，请近似该模型。

解：
（1）
clc;clear;
G=tf([2 0 0.5],[1 -0.1 3 0.66]);
Gc=tf([10 6],[1 0]);
H=tf(1,[1 2]);
G1=feedback(G*Gc,H)
G2=G*Gc/(1+G*Gc*H)
Gmin=minreal(G2)
结果：
Transfer function:
20 s^4 + 52 s^3 + 29 s^2 + 13 s + 6
s^5 + 1.9 s^4 + 22.8 s^3 + 18.66 s^2 + 6.32 s + 3
Transfer function:
20 s^8 + 50 s^7 + 83.8 s^6 + 179.3 s^5 + 126 s^4 + 57.54 s^3 + 26.58 s^2 + 3.96 s
s^9 + 1.8 s^8 + 25.61 s^7 + 22.74 s^6 + 74.11 s^5 + 73.4 s^4 + 30.98 s^3+ 13.17 s^2 + 1.98 s Transfer function:
20 s^4 + 52 s^3 + 29 s^2 + 13 s + 6
s^5 + 1.9 s^4 + 22.8 s^3 + 18.66 s^2 + 6.32 s + 3
(2)
由于
s c e s s s s G s G 23
2)1(36
24)(*)(-++= 方法1：将s e 2-转换为近似多项式。

clc;clear;
s=tf('s');
G=(24*s+36)/(s^2*(s+1)^3);
[num,den]=pade(2,2);
G1=feedback(tf(num,den)*G,1)
结果：
Transfer function:
24 s^3 - 36 s^2 - 36 s + 108
------------------------------------------------------------
s^7 + 6 s^6 + 15 s^5 + 19 s^4 + 36 s^3 - 33 s^2 - 36 s + 108
方法2：将G*Gc/(1+G*Gc*H)中的分母中的s e 2-转换为近似多项式。

clc;clear;
s=tf('s');
G=(24*s+36)/(s^2*(s+1)^3);
[num,den]=pade(2,2);
G1=feedback(G,tf(num,den));
G1.iodelay=2
Transfer function:
24 s^3 + 108 s^2 + 180 s + 108
exp(-2*s) * ------------------------------------------------------------------------------- s^7 + 6 s^6 + 15 s^5 + 19 s^4 + 36 s^3 - 33 s^2 - 36 s + 108
2. 假定系统为：
)(0001)(111000100001024269)(t u t x t x ⎥⎥⎥⎥⎦
⎤⎢⎢⎢⎢⎣⎡+⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡----= [])(2110)(t x t y =
请检查该系统是否为最小实现，如果不是最小实现，请从传递函数的角度解释该模型为何不是最小实现，并求其最小实现。

解：
clc;clear;
A=[-9 -26 -24 0;1 0 0 0;0 1 0 0;0 1 1 -1];
B=[1;0;0;0];
C=[0 1 1 2];
D=0;
G=ss(A,B,C,D);
sys=tf(G)
Gmin=minreal(sys)
G1=ss(Gmin)
sys=zpk(G)
结果：
传递函数表示：
Transfer function:
s^2 + 4 s + 3
---------------------------------------------
s^4 + 10 s^3 + 35 s^2 + 50 s + 24
最小实现后的传递函数模型：
Transfer function:
1
---------------------
s^2 + 6 s + 8
最小实现后的状态方程：
a =
x1 x2
x1 -6 -2
x2 4 0
b =
u1
x1 0.5
x2 0
c =
x1 x2
y1 0 0.5
d =
u1
y1 0
Continuous-time model.
不是最小实现的原因是分子分母没有对消
Zero/pole/gain:
(s+1) (s+3)
-----------------------
(s+1) (s+2) (s+3) (s+4)
3. 双输入双输出系统的状态方程：
)(20201000)()(20224264)(75.025.075.125
.1125.15.025.025.025.125.425.25.025.1525.2)(t x t y t u t x t x ⎥⎦⎤⎢⎣⎡=⎥⎥⎥⎥⎦
⎤⎢⎢⎢⎢⎣⎡+⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡------------= （1）试将该模型输入到MATLAB 空间，并求出该模型相应的传递函数矩阵。

（2）将该状态空间模型转化为零极点增益模型，确定该系统是否为最小实现模型。

如果不是，请将该模型的传递函数实现最小实现。

（3）若选择采样周期为s T 1.0=，求出离散后的状态方程模型和传递函数模型。

（4）对离散的状态空间模型进行连续变化，测试一下能否变回到原来的系统。

解：
（1）
clc;clear;
A=[2.25 -5 -1.25 -0.5;2.25 -4.25 -1.25 -0.25;...
0.25 -0.5 -1.25 -1;1.25 -1.75 -0.25 -0.75];
B=[4 6;2 4;2 2;0 2];
C=[0 0 0 1;0 2 0 2];
D=[0 0; 0 0];
G=ss(A,B,C,D)
Gtf=tf(G)
结果：
a =
x1 x2 x3 x4
x1 2.25 -5 -1.25 -0.5
x2 2.25 -4.25 -1.25 -0.25
x3 0.25 -0.5 -1.25 -1
x4 1.25 -1.75 -0.25 -0.75
b =
u1 u2
x1 4 6
x2 2 4
x3 2 2
x4 0 2
c =
x1 x2 x3 x4
y1 0 0 0 1
y2 0 2 0 2
d =
u1 u2
y1 0 0
y2 0 0
Continuous-time model.
Transfer function from input 1 to output...
s^2 + 3 s + 2.25
#1: --------------------------------------------------------- s^4 + 4 s^3 + 6.25 s^2 + 5.25 s + 2.25
4 s^3 + 14 s^2 + 22 s + 15
#2: ---------------------------------------------- s^4 + 4 s^3 + 6.25 s^2 + 5.25 s + 2.25
Transfer function from input 2 to output...
2 s^
3 + 6.5 s^2 + 7.75 s + 3.75
#1: --------------------------------------
s^4 + 4 s^3 + 6.25 s^2 + 5.25 s + 2.25
12 s^3 + 32 s^2 + 37 s + 17
#2: --------------------------------------
s^4 + 4 s^3 + 6.25 s^2 + 5.25 s + 2.25
（2）
clc;clear;
A=[2.25 -5 -1.25 -0.5;2.25 -4.25 -1.25 -0.25;...
0.25 -0.5 -1.25 -1;1.25 -1.75 -0.25 -0.75];
B=[4 6;2 4;2 2;0 2];
C=[0 0 0 1;0 2 0 2];
D=[0 0; 0 0];
G=ss(A,B,C,D);
Gtf=zpk(G)
sys=minreal(Gtf)
结果：
Zero/pole/gain from input 1 to output...
(s+1.5)^2
#1: ------------------------
(s+1.5)^2 (s^2 + s + 1)
4 (s+1.5) (s^2 + 2s + 2.5)
#2: ---------------------------
(s+1.5)^2 (s^2 + s + 1)
Zero/pole/gain from input 2 to output...
2 (s+1.5) (s^2 + 1.75s + 1.25)
#1: -------------------------------
(s+1.5)^2 (s^2 + s + 1)
12 (s+1) (s^2 + 1.667s + 1.417)
#2: --------------------------------
(s+1.5)^2 (s^2 + s + 1)
说明不是最小实现模型，最小实现模型如下：
Zero/pole/gain from input 1 to output...
1
#1: --------------
(s^2 + s + 1)
4 (s^2 + 2s + 2.5)
#2: ----------------------
(s+1.5) (s^2 + s + 1)
Zero/pole/gain from input 2 to output...
2 (s^2 + 1.75s + 1.25)
#1: -----------------------
(s+1.5) (s^2 + s + 1)
12 (s+1) (s^2 + 1.667s + 1.417)
#2: --------------------------------
(s+1.5)^2 (s^2 + s + 1)
（3）
clc;clear;
A=[2.25 -5 -1.25 -0.5;2.25 -4.25 -1.25 -0.25;...
0.25 -0.5 -1.25 -1;1.25 -1.75 -0.25 -0.75];
B=[4 6;2 4;2 2;0 2];
C=[0 0 0 1;0 2 0 2];
D=[0 0; 0 0];
G1=ss(A,B,C,D);
G2=tf(G1);
s1=c2d(G1,0.1)
s2=c2d(G2,0.1)
结果：
a =
x1 x2 x3 x4 x1 1.192 -0.4455 -0.1013 -0.04215
x2 0.2008 0.6124 -0.1058 -0.01884
x3 0.01526 -0.03499 0.8849 -0.09054
x4 0.1147 -0.1622 -0.01973 0.9279
b =
u1 u2
x1 0.3833 0.5527
x2 0.1906 0.3694
x3 0.1879 0.1764
x4 0.004833 0.1927
c =
x1 x2 x3 x4
y1 0 0 0 1
y2 0 2 0 2
d =
u1 u2
y1 0 0
y2 0 0
Sampling time: 0.1
Discrete-time model.
Transfer function from input 1 to output...
0.004833 z^3 - 0.003645 z^2 - 0.004467 z + 0.003463 #1: ---------------------------------------------------
z^4 - 3.617 z^3 + 4.908 z^2 - 2.962 z + 0.6703
0.3908 z^3 - 1.038 z^2 + 0.9242 z - 0.2754
#2: ----------------------------------------------
z^4 - 3.617 z^3 + 4.908 z^2 - 2.962 z + 0.6703 Transfer function from input 2 to output...
0.1927 z^3 - 0.5181 z^2 + 0.465 z - 0.1392
#1: ----------------------------------------------
z^4 - 3.617 z^3 + 4.908 z^2 - 2.962 z + 0.6703
1.124 z^3 - 3.078 z^2 +
2.817 z - 0.8611
#2: ----------------------------------------------
z^4 - 3.617 z^3 + 4.908 z^2 - 2.962 z + 0.6703 Sampling time: 0.1
（4）
clc;clear;
A=[2.25 -5 -1.25 -0.5;2.25 -4.25 -1.25 -0.25;...
0.25 -0.5 -1.25 -1;1.25 -1.75 -0.25 -0.75];
B=[4 6;2 4;2 2;0 2];
C=[0 0 0 1;0 2 0 2];
D=[0 0; 0 0];
G=ss(A,B,C,D);
s1=c2d(G,0.1);
sys1=d2c(s1)
结果：
a =
x1 x2 x3 x4
x1 2.25 -5 -1.25 -0.5
x2 2.25 -4.25 -1.25 -0.25
x3 0.25 -0.5 -1.25 -1
x4 1.25 -1.75 -0.25 -0.75
b =
u1 u2
x1 4 6
x2 2 4
x3 2 2
x4 1.387e-015 2
c =
x1 x2 x3 x4
y1 0 0 0 1
y2 0 2 0 2
d =
u1 u2
y1 0 0
y2 0 0
Continuous-time model.
4. 假设系统的传递函数模型为：
222
)(2+++=s s s s G
系统状态的初始值为⎥⎦
⎤⎢⎣⎡-21，假设系统的输入为t e t u 2)(-=。

（1）将该传递函数模型转化为状态空间模型。

（2）利用公式 ⎰--+=t
t t A t t A d Bu e t x e
t x 00)()()()(0)(τττ求解],0[t 的状态以及系统输出的解析解。

（3）根据上述的解析解作出s ]10,0[时间区间的状态以及系统输出曲线。

（4）采用lsim 函数方法直接作出s ]10,0[时间区间的状态以及系统输出曲线，并与（3）的结果作比较。

解：
（1）
clc;clear;
num=[1 2];
den=[1 2 2];
G=tf(num,den);
sys=ss(G)
结果：
a =
x1 x2
x1 -2 -2
x2 1 0
b =
u1
x1 2
x2 0
c =
x1 x2
y1 0.5 1
d =
u1
y1 0
Continuous-time model.
于是系统的状态空间方程为：
[])
(15.0)()(02)(0122)(t x t y t u t x t x =⎥⎦
⎤⎢⎣⎡+⎥⎦⎤⎢⎣⎡--=
（2）
clc;clear;
syms t t1 x
A=[-2 -2;1 0];
B=[2;0];
C=[0.5 1];
D=0;
x0=[1;-2];
C1=expm(A*t);
C2=int((expm(A*(t-t1)))*B*exp(-2*t1),'t1',0,t);
x=C1*x0+C2;
x1=x(1,:)
x2=x(2,:)
y=C*x
结果：
x1 =
exp(-t)*cos(t)+3*exp(-t)*sin(t)-2*(cos(t)^2+sin(t)^2-cos(t)*exp(t))/exp(t)^2
x2 =
-exp(-t)*sin(t)-2*exp(-t)*cos(t)+(cos(t)^2+sin(t)^2+exp(t)*sin(t)-cos(t)*exp(t))/exp(t)^2
y =
-3/2*exp(-t)*cos(t)+1/2*exp(-t)*sin(t)-(cos(t)^2+sin(t)^2-cos(t)*exp(t))/exp(t)^2+(cos (t)^2+sin(t)^2+exp(t)*sin(t)-cos(t)*exp(t))/exp(t)^2
(3)
clc;clear;
syms t t1 x
A=[-2 -2;1 0];
B=[2;0];
C=[0.5 1];
D=0;
x0=[1;-2];
C1=expm(A*t);
C2=int((expm(A*(t-t1)))*B*exp(-2*t1),'t1',0,t); x=C1*x0+C2;
y=C*x;
m=0:0.1:10;
x=subs(x,{t},{m});%给syms定义的变量t赋值y=subs(y,{t},{m});
x1=x(1,:);
x2=x(2,:);
plot(m,x1,m,x2,m,y);
legend('x1','x2','y')
grid on
结果：
（4）
clc;clear;
t=0:0.1:10;
A=[-2 -2;1 0];
B=[2;0];
C=[0.5 1];
D=0;
x0=[1;-2];
u=exp(-2*t);
[y,x]=lsim(A,B,C,D,u,t,x0);
plot(t,x(:,1),t,x(:,2),t,y);
legend('x1','x2','y')
5. 已知矩阵
⎥⎥⎥⎦
⎤⎢⎢⎢⎣⎡----=212332110A （1）取1:1.0:0=t ，利用expm(At)函数绘制求A 的状态转移矩阵，看运行的速度如何？
（2）采用以下程序绘制A 的状态转移矩阵的曲线，看运行的速度如何？ clc;clear;
A=[0 1 -1;-2 -3 3;2 1 -2];
t=0:0.1:2;
Nt=length(t);
for k=1:Nt
F(:,:,k)=expm(A*t(k));
end
z=reshape(F,[9,Nt]);
plot(t,z)
grid
title('系统的状态转移矩阵')
解：
（1）
syms t
A=[0 1 -1;-2 3 3;2 1 -2];
Q=expm(A*t)
计算结果复杂，速度慢！
（2）
clc;clear;
A=[0 1 -1;-2 3 3;2 1 -2];
t=0:0.1:2;
Nt=length(t);
for k=1:Nt
F(:,:,k)=expm(A*t(k));
end
z=reshape(F,[9,Nt]);
plot(t,z)
grid
title('系统的状态转移矩阵')。