哈工大机电系统控制基础大作业一

《机电系统控制基础》大作业一

基于MATLAB的机电控制系统响应分析

哈尔滨工业大学

2013年11月8日

1

作业题目

1. 用MATLAB 绘制系统2

()25

()()425

C s s R s s s Φ=

=++的单位阶跃响应曲线、单位斜坡响应曲线。 2. 用MATLAB 求系统2

()25

()()425

C s s R s s s Φ==++的单位阶跃响应性能指标：上升时间、峰值时间、调节时间和超调量。

3. 数控直线运动工作平台位置控制示意图如下：

X i

伺服电机原理图如下：

L R

(1)假定电动机转子轴上的转动惯量为J 1，减速器输出轴上的转动惯量为J 2，减速器减速比为i ，滚珠丝杠的螺距为P ，试计算折算到电机主轴上的总的转动惯量J ；

(2)假定工作台质量m ，给定环节的传递函数为K a ，放大环节的传递函数为K b ，包括检测装置在内的反馈环节传递函数为K c ，电动机的反电势常数为K d ，电动机的电磁力矩常数为K m ，试建立该数控直线工作平台的数学模型，画出其控制系统框图；

(3)忽略电感L 时，令参数K a =K c =K d =R=J=1，K m =10，P/i =4π，利用MATLAB 分析kb 的取值对于系统的性能的影响。

2

题目1 （1）操作过程

在Matlab 中进行编程，然后进行仿真,得到曲线。

（2）仿真结果及分析

00.5

1

1.5

22.533.5t(sec)

x (t )

00.2

0.4

0.6

0.8

1

1.2

1.4

t(sec)

x (t )

左图为单位阶跃响应曲线

右图为单位斜坡响应曲线

（3）源代码

t = [0:0.01:3.5];u = t nG = [25]; dG = [1 4 25]; G = tf(nG,dG); [y,T] = step(G,t); z = lsim(G,u,t);

subplot(121),plot(t,z,'-')

grid;xlabel('t(sec)'),ylabel('x(t)'); subplot(122),plot(T,y,'-') xlabel('t(sec)'),ylabel('x(t)'); grid on ;

3 题目2

（1）操作过程

按照要求编写程序，将写好的程序输入matlab中得到所要求的结果。

(2)仿真结果

（3）源代码

t=[0:0.001:10];

yss=1,dta=0.02;

num=25;

den=[1,4,25];

G1=tf(num,den);

y1=step(G1,t);

r=1;while y1(r)

tr1=(r-1)*0.001;

[ymax,tp]=max(y1);tp1=(tp-1)*0.001;

mp1=(ymax-yss)/yss;

s=5001;while y1(s)>1-dta & y1(s)<1+dta;s=s-1;end

ts1=(s-1)*0.001;

[tr1 tp1 ts1 mp1]

4 题目3

(一)公式推导

（1）机械部分

由公式：

各级轴折算转动惯量： 222311123224

(

)()()n n n

J J J J n n n ∑=+++

工作台折算转动惯量： 2()2L J m π

= 得：

2

122211()2P J J J m i i π

=+

+ （2）电路部分

①电枢回路电压平衡方程：

()

()()a a a di t u t L

Ri t E dt

=++ （1）

式中()a E V 是电枢反电势，其大小与激磁磁通及转速成正，d ()(//)a E k w t V rad s =。

②电磁转矩方程：

()()m m a M t k i t =

（2）

③电动机上的转矩平衡方程：

()

()()m c dw t J

M t M t dt

=- （3）

J 是电动机和负载折合到电动机轴上的转动惯量。

由式（2-1）、式（2-2）和式（2-3）中消去中间变量()a i t 、a E 及()m M t 便可得到以()w t 为输出量，以()u t 为输入量的微分方程：

()()

()()2222m d m d w t dw t LJ RJ k k w t k u t d t dt

++= （4）

对式（2-4）进行拉普拉斯变换，得：

22()2()()()m d m LJs W s RJsW s k k W s k U s ++=

所以

02()

()()22m m d

k W s G s U s LJs RJs k k =

=

++ 而由：

()()W s s s θ=

故，

12

1

()22m m d k G s LJs RJs k k si

=

++ 式中i 为减速比。

又:

2()2p G s π

=

控制系统框图：

（二） （1）操作过程

由已知条件电感L 时，令参数K a =K c =K d =R=J=1，K m =10，P/i =4π，所以

系统的开环传递函数为：

310()(5)

b

k G s s s =

+

系统的闭环传递函数为：

2

10()510b

b

k G s s s k =

++ （2）仿真结果及分析

Ⅰ.系统的时域分析:

n ω=

，?=

①过阻尼:

当1?<,即00.625b k <<时,见图1

②临界阻尼:

当0.625?=,即0.625b k =时,见图2

③欠阻尼:

当01?<<,即0.625b k >时,见图3

④零阻尼:

当0?=,即b k =+∞时,见图4

1

23

4

00.20.40.6

0.8

G

1(t )

1

234

00.5

1

G 2(t )

1

23

4

00.51

1.5

t(sec)

G 3(t )

1

234

00.511.5

2t(sec)

G 4(t )

由图4知，与书上结果基本符合。

再取欠阻尼时，取0.625,1,5,10b b b b k k k k ====时，得到单位阶跃响应如图5所示：

图3 图4

t(sec)

X 0(t )

图5

由图5可知，当b k 增大时，系统的最大超调量p M 和震荡次数都在增加，系统稳定性下降。但上升时间和调整时间都有减小，系统的快速性有所提高。 通过Bode 图分析系统的稳定性，如图6所示：

M a g n i t u d e (d B )

10

10

10

10

10

10

P h a s e (d e g )Bode Diagram Frequency (rad/s)

-50050M a g n i t u d e (d B )

10

10

10

10

10

10

-90-450P h a s e (d e g )Bode Diagram

Frequency (rad/s)

M a g n i t u d e (d B )

10

10

10

10

10

10

P h a s e (d e g )Bode Diagram

Frequency (rad/s)

-50050M a g n i t u d e (d B )

10

10

10

10

10

10

-90-450P h a s e (d e g )Bode Diagram

Frequency (rad/s)

图6

由Bode 图可知，幅频特性曲线无限趋近于2

π

-

线，但是无交点。由此可知，系统的幅值裕量g K =+∞，

相角裕量0γ>，由几何判据可知，系统稳定。

综上，当0.625b k <时，由Bode 图可知，系统恒稳定。但是，随着b k 的增加，最大超调量和振荡次数均增加，系统的相对稳定性降低，但上升时间和调整时间都有减小，系统的快速性有所提高。故而b k 的值只有维持在适当的范围，才能使系统获得良好的稳定性。 （3）源代码

图1-4代码： t = [0:0.01:4];

k = 0.1;nG1 = [10*k];dG1 = [1 5 10*k];G1=tf(nG1,dG1); k = 0.625;nG2 = [10*k];dG2 = [1 5 10*k];G2=tf(nG2,dG2); k = 1;nG3 = [10*k];dG3 = [1 5 10*k];G3=tf(nG3,dG3);

k = 9999999999999999999;nG4 = [10*k];dG4 = [1 5 10*k];G4=tf(nG4,dG4); [y1,T] = step(G1,t); [y2,T] = step(G2,t); [y3,T] = step(G3,t); [y4,T] = step(G4,t); subplot(221),plot(T,y1,'-') legend('kb=0.1')

xlabel('t(sec)'),ylabel('G1(t)');grid on ; subplot(222),plot(T,y2,'-') legend('kb=0.625')

xlabel('t(sec)'),ylabel('G2(t)');grid on ;

subplot(223),plot(T,y3,'-')

legend('kb=1')

xlabel('t(sec)'),ylabel('G3(t)');grid on;

subplot(224),plot(T,y4,'-')

legend('kb=9999')

xlabel('t(sec)'),ylabel('G4(t)');grid on;

图5代码：

t = [0:0.01:4];

k = 0.625;nG1 = [10*k];dG1 = [1 5 10*k];G1=tf(nG1,dG1); k = 1;nG2 = [10*k];dG2 = [1 5 10*k];G2=tf(nG2,dG2);

k = 10;nG3 = [10*k];dG3 = [1 5 10*k];G3=tf(nG3,dG3);

k = 50;nG4 = [10*k];dG4 = [1 5 10*k];G4=tf(nG4,dG4); [y1,T] = step(G1,t);

[y2,T] = step(G2,t);

[y3,T] = step(G3,t);

[y4,T] = step(G4,t);

plot(T,y1,'-',T,y2,'--',T,y3,'-.',T,y4,'r-')

legend('k=0.625','k=1','k=10','k=50')

xlabel('t(sec)'),ylabel('X0(t)');grid on;

图6代码：

k = 0.625;nG1 = [10*k];dG1 = [1 5];G1=tf(nG1,dG1);

k = 1;nG2 = [10*k];dG2 = [1 5];G2=tf(nG2,dG2);

k = 10;nG3 = [10*k];dG3 = [1 5];G3=tf(nG3,dG3);

k = 50; nG4 = [10*k];dG4 = [1 5];G4=tf(nG4,dG4);

w=logspace(-2,3,100);

subplot(221);bode(nG1,dG1,w);legend('k=0.625');grid on; subplot(222);bode(nG2,dG2,w);legend('k=1');grid on; subplot(223);bode(nG3,dG3,w);legend('k=10');grid on; subplot(224);bode(nG4,dG4,w);legend('k=50');grid on;