系统建模与仿真习题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('系统的状态转移矩阵')