Matlab与控制系统仿真部分习题答案

【4.2】 程序：

num=[5,0];den=conv([1,1],conv([1,2],[1,3]));

[numc,denc]=cloop(num,den);

[z,p,k]=tf2zp(numc,denc);

[A,B,C,D]=tf2ss(numc,denc);

g_zp=zpk(z,p,k)

g_tf=tf(numc,denc)

g_ss=ss(A,B,C,D)

运行结果：

Zero/pole/gain:

5 s

----------------------------------

(s+0.4432) (s^2 + 5.557s + 13.)

Transfer function:

5 s

----------------------

s^3 + 6 s^2 + 16 s + 6

a =

x1 x2 x3

x1 -6 -16 -6

x2 1 0 0

x3 0 1 0

b =

u1

x1 1

x2 0

x3 0

c =

x1 x2 x3

y1 0 5 0

d =

u1

y1 0

【4.3】程序：

A=[0 0 0 -1;1 0 0 -2;0 1 0 -3;0 0 1 -4];

B=[0;0;0;1];

C=[1 0 0 0];

g_ss=ss(A,B,C,D)

[num,den]=ss2tf(A,B,C,D);

g_tf=tf(num,den)

[z,p,k]=ss2zp(A,B,C,D);

g_zpk=zpk(z,p,k)

运行结果：

a =

x1 x2 x3 x4

x1 0 0 0 -1

x2 1 0 0 -2

x3 0 1 0 -3

x4 0 0 1 -4

b =

u1

x1 0

x2 0

x3 0

x4 1

c =

x1 x2 x3 x4

y1 1 0 0 0

d =

u1

y1 0

Continuous-time model.

Transfer function:

-3.109e-015 s^3 - s^2 - 3.331e-015 s - 4.441e-016

-------------------------------------------------

s^4 + 4 s^3 + 3 s^2 + 2 s + 1

Zero/pole/gain:

- s^2

----------------------------------------------

(s+0.6724) (s+3.234) (s^2 + 0.0936s + 0.4599)

【5.1】

（1） 程序

num=[0,10];

den=conv([1,0],[1,7,17]);

[numc,denc]=cloop(num,den,-1);

G=tf(numc,denc)

[y,t]=step(G);

plot(t,y,'b-')

C=dcgain(G);

n=1;

while y(n)<0.1*C

n=n+1;

end

m=1;

while y(m)<0.9*C

m=m+1;

end

risetime=t(m)-t(n)

[Y,k]=max(y);

percentovershoot=100*(Y-C)/C

i=length(t);

while(y(i)>0.98*C)&(y(i)<1.02*C)

i=i-1;

end

settlingtime=t(i)

运行结果：

Transfer function:

10

-----------------------

s^3 + 7 s^2 + 17 s + 10

risetime =

2.7312

percentovershoot =

-0.4399

settlingtime =

5.1372

10.90.80.70.60.50.40.30.20.1001234567图：

（2） 程序

k=[10,100,1000];

t=linspace(1,20,200);

num=1;

den=conv([1,0],[1,7,17]);

for j=1:3;

s1=tf(num*k(j),den);

sys=feedback(s1,1)

y(:,j)=step(sys,t);

end

plot(t,y(:,1),'r',t,y(:,2),'b',t,y(:,3),'g')

gtext('k=10');gtext('k=100');gtext('k=1000')

运行结果：

Transfer function:

10

-----------------------

s^3 + 7 s^2 + 17 s + 10

Transfer function:

100

------------------------

s^3 + 7 s^2 + 17 s + 100

Transfer function:

1000

-------------------------

s^3 + 7 s^2 + 17 s + 1000

图：

1.81.61.41.210.80.60.40.202468101214161820

图：

4x 10223210-1-2-302468101214161820

【6.1】程序：

（1）

num1=[1,1];

den1=conv([1,0,0],conv([1,2],[1,4]));

sys1=tf(num1,den1)

rlocus(sys1)

运行结果：

Root Locus82Imaginary Axis0-2-4-6-8-12-10-8-6-4Real Axis-2024

（2）

num2=[1,1];

den2=conv([1,0],conv([1,-1],[1,4,16]));

sys2=tf(num2,den2)

rlocus(sys2)

运行结果：

Root Locus82Imaginary Axis0-2-4-6-8-10-8-6-4Real Axis-2024

（3）

num3=[1,8];

den3=conv([1,0,0],conv([1,3],conv([1,5],conv([1,7],[1,15]))));

sys3=tf(num3,den3)

rlocus(sys3)

运行结果：

Root Locus2015105Imaginary Axis0-5-10-15-20-30-25-20-15-10-5051015Real Axis

【6.3】

程序：

num=[1,2];

den=conv([1,0],conv([1,4],conv([1,8],[1,2,5])));

sys=tf(num,den)

rlocus(sys)

[k,poles]=rlocfind(sys)

运行结果：

Transfer function:

s + 2

---------------------------------------

s^5 + 14 s^4 + 61 s^3 + 124 s^2 + 160 s

Select a point in the graphics window

selected_point =

0.0296 + 2.2826i

k =

135.8815

poles =

-7.3248

-5.4104

0.0145 + 2.3021i

0.0145 - 2.3021i

-1.2939

图：

Root Locus15105Imaginary Axis0-5-10-15-20-15-10-5Real Axis051015

【7.3】

程序（1）画波特图

num=[50];

den=conv([1,0],conv([1,10],[3,1]));

sys=tf(num,den)

sys1=feedback(sys,1)

bode(sys)

grid

图（1）

Bode Diagram10050Magnitude (dB)Phase (deg)0-50-100-150-90-135-180-225-27010-210-1100101102103Frequency (rad/sec)

程序（2）

画奈奎斯特图

num=[50];

den=conv([1,0],conv([1,10],[3,1]));

sys=tf(num,den)

sys1=feedback(sys,1)

nyquist(sys)

grid

图（2）

Nyquist Diagram3000 dB200100Imaginary Axis0-100-200-300-16-14-12-10-8Real Axis-6-4-20

程序（3）画零极点图

num=[50];

den=conv([1,0],conv([1,10],[3,1]));

sys=tf(num,den)

sys1=feedback(sys,1)

pzmap(sys1)

gird

图（3）

Pole-Zero Map1.510.5Imaginary Axis0-0.5-1-1.5-12-10-8-6Real Axis-4-20

程序（4）

计算相角裕量和幅值裕量

num=[50];

den=conv([1,0],conv([1,10],[3,1]));

sys=tf(num,den)

sys1=feedback(sys,1)

[gm,pm,wcg,wcp]=margin(sys)

运行结果

Transfer function:

50

---------------------

3 s^3 + 31 s^2 + 10 s

Transfer function:

50

--------------------------

3 s^3 + 31 s^2 + 10 s + 50

gm =

2.0667

pm =

7.5615

wcg =

1.8257

wcp =

1.25

程序（5）

绘制阶跃响应曲线

num=[50];

den=conv([1,0],conv([1,10],[3,1]));

sys=tf(num,den)

sys1=feedback(sys,1)

step(sys1)

图（5）

Step Response21.81.61.41.2Amplitude10.80.60.40.200102030Time (sec)40506070

【7.4】

程序如下：

num=[300];

den=conv([1,0,0],conv([0.2,1],[0.02,1]));

sys=tf(num,den)

margin(sys)

grid

波特图如下：

Bode DiagramGm = Inf , Pm = -78 deg (at 11 rad/sec)10050Magnitude (dB)Phase (deg)0-50-100-150-180-225-270-315-36010-1100101102103Frequency (rad/sec)

【9.3】程序：

A=[-2 2 -1;0 -2 0;1 -4 0];

B=[0;0;1];

C=[1,0,0];

D=0;

M=ctrb(A,B)

m=rank(M)

if m==3;

disp('系统可控')

else

disp('系统不可控')

end

N=obsv(A,C)

n=rank(N)

if n==3;

disp('系统可观')

else

disp('系统不可观')

end

sys=ss(A,B,C,D)

[num,den]=ss2tf(A,B,C,D)

sys1=tf(num,den)

[z,p,k]=ss2zp(A,B,C,D)

运行结果：

M =

0 -1 2

0 0 0

1 0 -1

m =

2

系统不可控

N =

1 0 0

-2 2 -1

3 -4 2

n =

2

系统不可观

a =

x1 x2 x3

x1 -2 2 -1

x2 0 -2 0

x3 1 -4 0

b =

u1

x1 0

x2 0

x3 1

c =

x1 x2 x3

y1 1 0 0

d =

u1

y1 0

Continuous-time model.

num =

0 0 -1 -2

den =

1 4 5 2

Transfer function:

-s - 2

---------------------

s^3 + 4 s^2 + 5 s + 2

z =

-2

p =

-1

-1

-2

k = -1

图形：

【10.1】

（1）程序：

A=[0,1,0,0;0,5,0,0;0,0,-7,0;0,0,0,-8];

B=[0;1;0;1];

C=[1,2,3,4];

D=zeros(1,1);

G_ss=ss(A,B,C,D)

运行结果：

a =

x1 x2 x3 x4

x1 0 1 0 0

x2 0 x3 0 x4 0 b =

u1

x1 0

x2 1

x3 0

x4 1

c =

5 0 0 -7 0 0 0

0

-8

x1 x2 x3 x4

y1 1 2 3 4

d =

u1

y1 0

（2）：程序：

[num1,den1]=ss2tf(A,B,C,D);

p=roots(den1)

i=0;

for k=1:1:length(p)

if real(p(k))>0

i=i+1;

end

end

if i>0

disp('系统不稳定');

else

disp('系统稳定');

end

运行结果：

p =

0

5.0000

-8.0000

-7.0000

系统不稳定

（3）（4）程序：

AA=[0,1,0;0,5,0;0,0,-8];

BB=[0;1;1];

P=[-1,-2,-8];

K=acker(AA,BB,P);

i=4;

K(4)=0;

K

pp=eig(A-B*K)

sys1=tf(num1,den1);

[y1,t]=step(sys1);

plot(t,y1)

hold on

A_feedback=A-B*K;

[num2,den2]=ss2tf(A_feedback,B,C,D);

sys2=tf(num2,den2);

[y2,t]=step(sys2);

plot(t,y2,'r')

grid

gtext('反馈前')

gtext('反馈后')

运行结果：

K =

2 8 0 pp =

-8

0

-2

-1

-7

图形：

6反馈前321反馈后00123456

【13.1】程序：

A=[0,1;0,0];

B=[0;1];

C=[1,0];

D=zeros(1,1);

G_ss=ss(A,B,C,D)

M=ctrb(A,B);

if rank(M)==2

disp('系统完全能控');

else

disp('系统不完全能控');

end

S=[1,0];

N=obsv(A,S);

if rank(N)==2

disp('(A,S)可观测');

else

disp('(A,S)不可观测');

end

R=1;

Q=[1,0;0,0];

[K,P,E]=Lqr(A,B,Q,R)

A_new=A-B*K;

G_new=ss(A_new,B,C,D);

t=linspace(0,5,100)';

y1=step(G_ss,t);

y2=step(G_new,t);

plot(t,y1,'r:',t,y2,'b-')

grid

gtext('反馈前')

gtext('反馈后')

运行结果：

a =

x1 x2

x1 0 1

x2 0 0

b =

u1

x1 0

x2 1

c =

x1 x2

y1 1 0

d =

u1

y1 0

Continuous-time model.

系统完全能控

(A,S)可观测

K =

1.0000 1.4142

P =

1.4142 1.0000

1.0000 1.4142

E =

-0.7071 + 0.7071i

-0.7071 - 0.7071i

图形：

14121086反馈前42反馈后000.511.522.533.4.55