中国矿业大学《控制系统计算机仿真》实验试题答案

实验一 1、37 例2.1）已知一阶微分方程：

试编写程序，用欧拉法求 ，即

的值；并画出图形。 解：程序如下： clc clear

Tf=input('仿真时间Tf='); %输入仿真时间 h=input('计算步长h='); %输入仿真步长 x0=1/3; t=0; x=x0; %输入初始值 for i=1:Tf/h K1=-30*x0; x1=x0+h*K1;

x=[x;x1]; %变量以向量形式保存

t=[t;t(i)+h]; %对应时刻以向量形式保存 x0=x1; end

[t,x] %以数据形式输出 plot(t,x) %以曲线形式输出 所绘图形如下：

00.20.40.60.81 1.2 1.4 1.6

-0.2

-0.1

0.1

0.2

0.3

0.4

0.5

X: 1.5

Y: 3.104e-010

2、2. （教材P79习题2.5）已知系统的状 态方程和输出方程为： 式中 u(t)=1(t)。初始条件为： x1(0)=x2(0)=x3(0)=0。取h=0.05， 试用RK4法求t=0.5时的y 值。

解：程序如下：

05

.0,3/1)0(,

30==-=?

h y y y )5.1(=t y 30y ????=??

????

????+????????????????????---=??????????21321321]0 0 1[1040 0 0 121 0 190 1

8x x y u x x x x x x

clc clear u=1;

x0=0; t=0;y=0; x0=[0;0;0];

A=[-8 1 0;-19 0 1;-12 0 0]; B=[0;4;10]; C=[1 0 0];

Tf=input('仿真时间Tf='); h=input('计算步长h='); x=x0;

for i=1:Tf/h K1=A*x+B*u;

K2=A*(x+h*K1/2)+B*u; K3=A*(x+h*K2/2)+B*u; K4=A*(x+h*K3)+B*u;

x=x+h*(K1+2*K2+2*K3+K4)/6; y=[y;C*x]; t=[t;t(i)+h]; end y

plot(t,y)

所绘图形如下：

05101520253035

0.10.20.30.4

0.50.6

0.70.80.900.050.10.150.20.250.30.350.40.450.5

00.05

0.1

0.150.2

0.25

X: 0.5Y: 0.2201

3、已知系统的传递函数， 试求其状态方程及输出方程。 （教材P26页 习题1.7 ）

解：程序如下： den=[1 3];

num=[1 3 3 1];

[A,B,C,D]=tf2ss(den,num) 运行结果如下： A =

-3 -3 -1 1 0 0 0 1 0 B =

1 0 0 C =

0 1 3 D =

4、已知线性定常系统的状态方程

及输出方程，试求系统的传递函数。 （教材P26页 习题1.8 ） 解：程序如下：

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

B=[0;0;1];

C=[1 0 1];

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

运行结果如下：

1

333

)(2

3

++++=s s s s s G ????

??????=??

????????+????????????????????--=??????????321321321]1 0 1[100 2- 0 0 0 1 0 0 1 1x x x y u x x x x x x

den =

0 1 2 1

num =

1 4 5 2

5、（教材P48页例2.4 ）如图所示的直流电机拖动系统，试在SIMULINK环境下建立系统仿真模型：

1）采用欧拉法与四阶龙格库塔法进行仿真，并确定其步长的稳定区间，取不同步长比较仿真结果；

2）采用RK4研究外环PI控制器参数的改变(0.17 0.5 1 1.5)，观其阶跃响应对系统动态参数的影响。

（1）、欧拉法

clc

clear

Tf=10; %输入仿真时间

h=input('计算步长h='); %输入仿真步长

P=[1 0.01 0.1 0;

0 0.85 1 0.17;

1 0.01 1 0;

0 0.051 1 0.15;

1 0.0067 70 0;

1 0.15 0.21 0;

0 1 130 0;

1 0.01 0.1 0;

1 0.01 0.0044 0]; %输入各环节参数

W0=zeros(9,1);W0(1,1)=1; %输入外部链接矩阵

W=zeros(9,9);W(2,1)=1;W(2,9)=-1;W(3,2)=1;W(4,3)=1;

W(4,8)=-1;W(5,4)=1;W(6,5)=1;W(6,7)=-0.212;W(7,6)=1;

W(8,6)=1;W(9,7)=1; %输入系统连接矩阵

A=diag(P(:, 1));B=diag(P(:, 2));

C=diag(P(:, 3));D=diag(P(:, 4)); %生成系数矩阵 Q=B-D*W;Q1=inv(Q); R=C*W-A; V=C*W0;

A1=Q1*R; B1=Q1*V; %生成闭环系数矩阵 C1=[0 0 0 0 0 0 1 0 0];

u=1; %阶跃输入幅值

x0=zeros(9,1); t=0;y=0; %设置初值 x=x0;

for i=1:Tf/h

K1=A1*x+B1*u; x=x+h*K1; y=[y;C1*x]; t=[t;t(i)+h]; end

plot(t,y)

计算步长h=0.011 计算步长h=0.012

0123456789100

51015

20

25

303540

012345678910

-2

-1.5-1-0.5

0.511.522.5x 10

50

四阶龙格库塔法： clc clear

Tf=input('仿真时间Tf='); %输入仿真时间 h=input('计算步长h='); %输入仿真步长 P=[1 0.01 0.1 0; 0 0.85 1 0.17; 1 0.01 1 0;

0 0.051 1 0.15; 1 0.0067 70 0; 1 0.15 0.21 0; 0 1 130 0;

1 0.01 0.1 0;

1 0.01 0.0044 0]; %输入各环节参数 W0=zeros(9,1);W0(1,1)=1; %输入外部链接矩阵

W=zeros(9,9);W(2,1)=1;W(2,9)=-1;W(3,2)=1;W(4,3)=1; W(4,8)=-1;W(5,4)=1;W(6,5)=1;W(6,7)=-0.212;W(7,6)=1; W(8,6)=1;W(9,7)=1; %输入系统连接矩阵 A=diag(P(:, 1));B=diag(P(:, 2));

C=diag(P(:, 3));D=diag(P(:, 4)); %生成系数矩阵 Q=B-D*W;Q1=inv(Q); R=C*W-A; V=C*W0;

A1=Q1*R; B1=Q1*V; %生成闭环系数矩阵 C1=[0 0 0 0 0 0 1 0 0];

u=1; %阶跃输入幅值

x0=zeros(9,1); t=0;y=0; %设置初值 x=x0;

for i=1:Tf/h

K1=A1*x+B1*u;

K2=A1*(x+h*K1/2)+B1*u; K3=A1*(x+h*K2/2)+B1*u; K4=A1*(x+h*K3)+B1*u;

x=x+h*(K1+2*K2+2*K3+K4)/6; y=[y;C1*x]; t=[t;t(i)+h]; end

plot(t,y)

仿真时间Tf=20 仿真时间Tf=5 仿真时间Tf=100 计算步长h=0.015 计算步长h=0.015 计算步长h=0.016

02468101214161820

510152025303540

00.51 1.52 2.53 3.54 4.55

510

1520253035

40020406080100120

-18

-16

-14-12

-10-8-6-4-202x 10

302

（2）、

clc

clear

Tf=input('仿真时间Tf='); %输入仿真时间

h=input('计算步长h='); %输入仿真步长

P=[1 0.01 0.1 0;

0.17 0.5 1 1.5;

1 0.01 1 0;

0 0.051 1 0.15;

1 0.0067 70 0;

1 0.15 0.21 0;

0 1 130 0;

1 0.01 0.1 0;

1 0.01 0.0044 0]; %输入各环节参数

W0=zeros(9,1);W0(1,1)=1; %输入外部链接矩阵

W=zeros(9,9);W(2,1)=1;W(2,9)=-1;W(3,2)=1;W(4,3)=1; W(4,8)=-1;W(5,4)=1;W(6,5)=1;W(6,7)=-0.212;W(7,6)=1; W(8,6)=1;W(9,7)=1; %输入系统连接矩阵

A=diag(P(:, 1));B=diag(P(:, 2));

C=diag(P(:, 3));D=diag(P(:, 4)); %生成系数矩阵

Q=B-D*W;Q1=inv(Q);

R=C*W-A; V=C*W0;

A1=Q1*R; B1=Q1*V; %生成闭环系数矩阵

C1=[0 0 0 0 0 0 1 0 0];

u=1; %阶跃输入幅值

x0=zeros(9,1); t=0;y=0; %设置初值

x=x0;

for i=1:Tf/h

K1=A1*x+B1*u;

K2=A1*(x+h*K1/2)+B1*u;

K3=A1*(x+h*K2/2)+B1*u;

K4=A1*(x+h*K3)+B1*u;

x=x+h*(K1+2*K2+2*K3+K4)/6;

y=[y;C1*x];

t=[t;t(i)+h];

end

plot(t,y)

仿真时间Tf=20 计算步长h=0.015

02468101214161820

5

10

15

20

25

30

仿真时间Tf=5 计算步长h=0.015

00.51 1.52 2.53 3.54 4.55

5

10

15

20

25

30

SIMULINK 仿真

0.0044

0.01s+1

0.10.01s+1

130s

0.210.15s+1

700.00167s+1

0.15s+10.051s

10.01s+1

0.17s+10.85s

0.10.01s+1

y1To Workspace

Saturation1

0.212

0.0044

0.01s+1

0.10.01s+1

130s

0.210.15s+1

700.00167s+1

0.15s+10.051s

10.01s+1

0.17s+10.85s

0.10.01s+1

y1To Workspace

0.212

实验二

1、已知某四阶非线性系统结构框图如图1所示。当系统输入幅值为10的阶跃信号时，试仿真分析系统输出响应：

1）在SIMULINK 环境下建立仿真模型，输出仿真结果；

2）采用MATLAB 语言编程实现仿真程序，与1)进行结果比较；

3）在第一环节前加饱和非线性环节，分析限幅C1=6、2时系统的输出响应，并与没限幅前比较其超调、峰值时间及调节时间。

解：（1）SIMULINK 图形如下

1s

s+0.785s+0.16

10.1s+1

10.5s+1

1s

s+0.785s+0.16

10.1s+1

10.5s+1

10.1s+1

10.5s+1

1s

s+0.785s+0.16

Saturation1

Saturation

1

0.1s+1

10.5s+11s s+0.785s+0.16

（2）

子程序：satura.m

function uc=satura(ur,s1)

%ur=1;s1=5;

if (ur>=s1)

uc=s1;

elseif (ur<=-1*s1)

uc=-1*s1;

else uc=ur;

end

子程序：dead.m

function uc=dead(ur,s1)

%ur=-1.5; s1=0.5;

if (abs(ur)>=s1)

if (ur>0)

uc=ur-s1;

else uc=ur+s1;

end

else uc=0;

End

主程序：

clc

clear

Tf=20; %输入仿真时间

h=0.01; %输入仿真步长

z=input('非线性标志向量z='); %输入非线性标志向量

c1=input('饱和非线性环节限幅c1=');%输入饱和非线性环节参数

a1=6;%输入死区非线性环节参数

P=[0.16 1 0.785 1;

0 1 1 0;

1 0.5 1 0;

1 0.1 1 0]; %输入各环节参数

A=P(:,1);B=P(:,2);C=P(:,3);D=P(:,4);

W0=zeros(4,1);W0(1,1)=1; %输入外部链接矩阵

W=zeros(4,4);W(1,4)=-1;W(2,1)=1;W(3,2)=1;W(4,3)=1;

for i=1:4

if(A(i)==0);

F(i)=1;

G(i)=h*C(i)/B(i);

H(i)=0.5*h*h*C(i)/B(i);

if (D(i)==0)

C1(i)=1;D1(i)=0;

else

C1(i)=1;D1(i)=D(i)/B(i); %求积分、比例积分环节离散系数end

else

F(i)=exp(-h*A(i)/B(i));

if (D(i)==0)

G(i)=(1- F(i))*C(i)/A(i);

H(i)=h*C(i)/A(i)-G(i)*B(i)/A(i);

C1(i)=1;D1(i)=0;

else

G(i)=(1- F(i))*D(i)/A(i);

H(i)=h*D(i)/A(i)-G(i)*B(i)/A(i);

C1(i)=C(i)/D(i)-A(i)/B(i);

D1(i)=D(i)/B(i); %求惯性、比例惯性环节离散系数end

end

end

%F1=diag(F(:, 1));G1=diag(G(:, 1));H1=diag(H(:, 1));

%C2=diag(C1(:, 1));D2=diag(D1(:, 1));

%求各环节输入

Y=zeros(4,1);X=Y;y=zeros(1,4);r=10;Uk=zeros(4,1);

t=0:h:Tf; N=length(t); %求采样点个数

for k=1:N-1

Ub=Uk; %保存前一次输入

Uk=W*Y+W0*r; %求当前个环节输入v

for i=1:4;

if(z(i)~=0)

if(z(i)==1)

Uk(i,1)=satura(Uk(i,1),c1);

end

if (z(i)==2)

Uk(i,1)=dead(Uk(i,1),a1); end end end

Udot=(Uk-Ub)/h; %求当前各环节输入导数 Uf=2*Uk-Ub; %求下一拍输入 %求各环节状态变量及输出

X=F'.*X+G'.*Uk+H'.*Udot; Y=C1'.*X+D1'.*Uf; for i=1:4;

if (z(i)~=0) if (z(i)==3)

Y(i,1)=satura(Y(i,1),c1); end

if (z(i)==4)

Y(i,1)=dead(Y(i,1),a1); end end end

y=[y;Y']; end

y4=y(:,4); plot(t,y4)

非线性标志向量z=[0 0 0 0] 饱和非线性环节限幅c1=6

02468101214161820

246810121416X: 2.74Y: 15.17

（3）

非线性标志向量z=[1 0 0 0] 饱和非线性环节限幅c1=6

02468101214161820

5

10

15

X: 3.12Y: 14.7

非线性标志向量z=[1 0 0 0] 饱和非线性环节限幅c1=2

2468101214161820

02

4

6

8

10

12

14

X: 4.91Y: 13.18

2、 某控制系统如图所示，选择增益K 的值，使系统阶跃响应的超调整量小于20%，且调节时间小于5 s 。（教材P80习题2.11）

K(s+1)(s+1)(s+8)(s+20)Zero-Pole

1s +3.2s +3.56s 32Transfer Fcn

Step

Scope

解：SIMULINK 图形如下：

1550(s+1)(s+1)

(s+8)(s+20)

Zero-Pole

1

s +3.2s +3.56s 32Transfer Fcn

Step

Scope

670(s

+1)(s +1)(s +8)(s +20)

Zero-Pole

1

s +3.2s +3.56s 32Trans fer Fcn

Step

Scope

3、某非线性系统如图所示，试求 时系统的动态特性。 (P52

页例题2.7)

解：SIMULINK 图形如下：

2s+1s +2s +s 32Transfer Fcn

Step

Scope

Saturation

1

Gain

)(12)(t t r ?=2s +1s +2s +s

32Trans fer Fcn

Step Scope

Saturation

5Gain

2s +1s +2s +s 32Trans fer Fcn

Step Scope

MATLAB

Function MATLAB Fcn

5Gain

2s +1s +2s +s 32T rans fer Fcn1

2s +1s +2s +s 32T rans

fer Fcn

Step1

Step

Scope

Saturation

5

Gain1

5

Gain

u y

fcn Embedded

MAT LAB Function

其中：

function y = fcn(u) %#eml if u>=1 y=1;

elseif u<=-1 y=-1; else y=u; end

4、考虑如图所示的单位反馈系统，试用双线性变换和根匹配法对其进行仿真（r(t)=1(t),T=0.05)：

（1）通过“c2d ”指令将系统各环节进行离散化；

（2）通过“tools —Control Design – Model Discretize ”将连续系统离散化； （3）试绘制闭环频率特性曲线。

解：（1） t=0.05;

n1=[1 3.33];d1=[1 24];[nd1,dd1]=c2dm(n1,d1,t,'tustin') n2=[10];d2=[1 0];[nd2,dd2]=c2dm(n2,d2,t,'tustin') n3=[1];d3=[1 10 16];[nd3,dd3]=c2dm(n3,d3,t,'tustin') 运行结果： nd1 =

0.6770 -0.5730

dd1 =

1.0000 -0.2500

nd2 =

0.2500 0.2500

dd2 =

1 -1

nd3 =

1.0e-003 *

0.4960 0.9921 0.4960

dd3 =

1.0000 -1.5714 0.6032

10s Transfer Fcn4s+3.33s+24Transfer Fcn31s +10s+16

2Transfer Fcn1Step

57.66

（2）

(

s +3.33)(s +24)t u st i n Zer o -P ol e

1s +10s+162t u st i n Tr a nsf e r F cn1

10s t u st i n Tr a nsf e r F cn

St e p

Scope

57.66

Gai n