自动控制原理Matlab时域分析实验.

《自动控制原理》课程实验报告

实验名称基于MATLAB仿真的系统时域分析专业班级 11级过程自动化

学

号

姓名

指导教师李离

学院名称电气信息学院

2012 年 11 月 20 日

基于MATLAB仿真的系统时域分析

一、实验目的

（1）学习如何利用MATLAB 分析控制系统的时域性能和比较系统的近似模型和实际模型；

（2）巩固系统绝对稳定性和相对稳定性的概念；

（3）掌握利用MATLAB 进行Routh-Hurwitz稳定性检验的方法; （4）学习利用MATLAB进行系统参数设计的方法。

二、实验设备

（1）硬件：个人计算机；

（2）软件：MATLAB 仿真软件（版本6.5 或以上）。

三、实验内容和步骤

本实验借助MATLAB 仿真，分析控制系统关于给定输入信号的瞬态响应和稳态跟踪误差，观察系统所实现的性能指标水平；同时，观察系统简化带来的性能变化情况。

验2-1实标准二阶系统的阶跃响应及性能分析

考虑图2.1 所示的标准二阶系统。假设1

w，观察当ζ=0.1、

n

0.2、0.4、0.7、1.0、2.0 时的系统单位阶跃响应, 估计各自对应的性能水平，并将其与理论值进行比较。

图2.1

源程序代码：

t=[0:0.1:12]; num=[1];

zeta1=0.1; den1=[1 2*zeta1 1]; sys1=tf(num,den1);

zeta2=0.2; den2=[1 2*zeta2 1]; sys2=tf(num,den2);

zeta3=0.4; den3=[1 2*zeta3 1]; sys3=tf(num,den3);

zeta4=0.7; den4=[1 2*zeta4 1]; sys4=tf(num,den4);

zeta5=1.0; den5=[1 2*zeta5 1]; sys5=tf(num,den5);

zeta6=2.0; den6=[1 2*zeta6 1]; sys6=tf(num,den6);

%

[y1,T1]=step(sys1,t); [y2,T2]=step(sys2,t);

[y3,T3]=step(sys3,t); [y4,T4]=step(sys4,t);

[y5,T5]=step(sys5,t); [y6,T6]=step(sys6,t);

plot(T1,y1,T2,y2,T3,y3,T4,y4,T5,y5,T6,y6)

%

xlabel('\omega_n t'), ylabel('y(t)')

title('\zeta = 0.1, 0.2, 0.4, 0.7, 1.0, 2.0'), grid

运行结果：

24

681012

00.20.40.60.81

1.21.41.6

1.8ωn t

y (t )

ζ = 0.1, 0.2, 0.4, 0.7, 1.0, 2.0

结果分析：

可以看出图中从上到下分别是1=n w ，ζ=0.1、 0.2、0.4、0.7、1.0、2.0 的图线。从图中我们可以得出：

系统逐渐从欠阻尼系统过渡到临界阻尼系统再到过阻尼系统，系统随阻尼比的增大，上升时间变长，调整时间变短，超调量减小，系统越稳定。 由理论公式：

上升时间21?

π-Θ

-=n r w t 调整时间n s w t ?4

= 超调量2

1??π

σ--

=e p 可知

结果与理论分析是一致的。

实验2-2标准二阶系统的脉冲响应

仍然考虑图 2.1 所示系统和假设1=n w 观察当ζ = 0.1、 0.25、0.5、1.0时的系统单位脉冲响应。 源程序代码：

t=[0:0.1:10]; num=[1];

zeta1=0.1; den1=[1 2*zeta1 1]; sys1=tf(num,den1); zeta2=0.25; den2=[1 2*zeta2 1]; sys2=tf(num,den2); zeta3=0.5; den3=[1 2*zeta3 1]; sys3=tf(num,den3); zeta4=1.0; den4=[1 2*zeta4 1]; sys4=tf(num,den4); %

[y1,T1]=impulse(sys1,t); [y2,T2]=impulse(sys2,t); [y3,T3]=impulse(sys3,t); [y4,T4]=impulse(sys4,t); %

plot(t,y1,t,y2,t,y3,t,y4)

xlabel('\omega_n t'), ylabel('y(t)/\omega_n') title('\zeta = 0.1, 0.25, 0.5,1.0'), grid

运行结果：

12345678910

-0.8-0.6

-0.4

-0.2

0.2

0.4

0.6

0.8

1

ωn

t

y (t )/ωn

ζ = 0.1, 0.25, 0.5,1.0

y1 v s. t

y2 v s. t y3 v s. t

y4 v s. t

结果分析：

从图中可知从上到下依次是1=n w 时ζ = 0.1、 0.25、0.5、1.0的脉冲响应。从图中我们可以得知：

?越小，超调量越大震荡的越厉害，频率越大响应的时间越快。

由二阶系统的脉冲响应()()

θ??

?+--=-t w e w t y n t w n n 22

1sin 1可知实验的结

果与理论的分析是一致的。

实验2-3移动机器人驾驶控制系统关于三角波输入的响应 移动机器人驾驶控制系统如图 2.3 所示。其中()s K K s G 211+= 运行程序 观察当系统输入如图 2.4 所示时的系统响应, 估计其稳态误差，并将其与理论值进行比较。利用函数 lsim 可对闭环系统关于斜坡输入的响应进行仿真。

图2.3 移动机器人驾驶控制系统

图2.4 移动机器人控制系统

源程序代码：

numg=[10 20]; deng=[1 10 0]; sysg=tf(numg,deng);

sys=feedback(sysg,[1]);

t1=[0:0.1:2]';

t2=[2.1:0.1:6]';

t3=[6.1:0.1:8]';

t=[0:0.1:8];

v1=[0:0.1:2]';v2=[1.9:-0.1:-2]';v3=[-1.9:0.1:0]';

% v1=[0:0.1:2]';v2=[2:-0.1:-2]';v3=[-2:0.1:0]';

u=[v1;v2;v3];

[y,T]=lsim(sys,u,t);

plot(t1,v1,'b--',t2,v2,'b--',t3,v3,'b--',T,y,'k-'),

%figure(1);

%plot(t1,v1)

%hold on; ishold;

%plot(t2,v2)

%hold on; ishold;

%plot(t3,v3)

%hold on; ishold;

xlabel('Time (seconds)'), ylabel('\theta (radians)'), grid %hold off; ishold;

xlabel('Time (seconds)'), ylabel('\theta (radians)'), grid %hold off; ishold

运行结果：

0123

45

678

-2

-1.5-1-0.500.5

11.5

2Time (seconds)

θ (r a d i a n s )

结果分析：

从运行的结果可以估计稳定误差5.0=ss e 。 从理论分析：定义系统的开环传递函数为

())

11.0()

15.0(2102010++=++s s s s s s ，系

统为Ⅰ型系统，开环增益K=2，当输入斜波函数时，系统的稳定误差

5.01

==

K

e ss . 实验2-4高阶模型的低阶近似

三阶系统 ()6

1166

2

3+++=

s s s s H （1） 的二阶近似模型为 ()6

.1584.26

.1L 2++=s s s （2）

运行程序 ，观察系统（1）和（2）的单位阶跃响应, 并就其各个性能指标水平进行比较。 源程序代码

num1=[6]; den1=[1 6 11 6]; sys1=tf(num1,den1);

num2=[1.6]; den2=[1 2.584 1.6]; sys2=tf(num2,den2); t=[0:0.1:8];

[y1,T1]=step(sys1,t); [y2,T2]=step(sys2,t);

plot(T1,y1,T2,y2,'--'), grid

xlabel('Time (seconds)'), ylabel('Step Response')

运行结果：

1

2

3

45

6

7

8

00.10.20.30.40.50.60.7

0.80.91Time (seconds)

S t e p R e s p o n s e

y1y2

结果分析：

从图中可知原三阶系统为下面的y1，化简后的二阶系统为上面的y2。

从图中可以看出两个系统都没有超调量，三阶系统的上升时间r

t 约等于（3.4-0.7）=2.7s,调节时间s t （%5±）约等于4s ；二阶系统的上升时间r t 约等于（3.1-0.4）=2.7s ，调节时间s t （%5±）约等于3.9s 。可见三阶系统化简成二阶系统对阶跃的响应相似。

实验2-5Routh-Hurwitz 稳定性检验

Routh-Hurwitz 稳定判据是一个关于系统稳定性的充要判据。如果系统特征方程的系数均已确定，则其在左半 s 平面上、右半 s 平面上以及 s 平面虚轴上根的数目可由 Routh-Hurwitz 稳定判据来确定。

调用 Matlab 函数 pole 和 roots, 可通过直接求解系统特征方程的根（即闭环传递函数的极点）来验证利用 Routh-Hurwitz 稳定判据得到的结果。

本项实验内容为：首先对下述系统或系统的特征方程运用 Routh-Hurwitz 稳定判据判断其特征根在s 平面上的分布情况，然后编写 Matlab 仿真程序加以验证。

⑴010*********=+++++s s s s s ⑴01222345=+++++s s s s s ⑶某闭环系统如图 2.5 所示；其中()()

50

351024

23+++=

s s s s s G

图2.5

⑴Routh-Hurwitz 判稳

10

10-72-2410

12-4610

42

11210

2

12

3

45s

s s s s s ε

εεεεε

其中ε为很小但大于0的数。

劳斯表中第一列元素符号的变化的次数为两次，说明特征方程有两个根在s 平面的右边，所以系统是不稳定的。 源程序代码

sys=[1 2 2 4 11 10]; r=roots(sys)

运行的结果

r =

0.8950 + 1.4561i 0.8950 - 1.4561i -1.2407 + 1.0375i -1.2407 - 1.0375i -1.3087 + 0.0000i

结果分析

运行得到有两个根0.8950 + 1.4561i 0.8950 - 1.4561i 在s 平面的右边，与Routh-Hurwitz 判稳的结果相同。 ⑵ Routh-Hurwitz 判稳

1

1

11

211210

1234

5s s s s s s ε

ε

ε ()01122

224=+=++s s s 其中ε为很小但大于0的数。

由劳斯表中4s 行都为0，可知特征方程在s 平面的虚轴上有两对重根。 源程序代码：

sys=[1 1 2 2 1 1]; r=roots(sys)

运行结果：

r =

-1.0000 + 0.0000i 0.0000 + 1.0000i 0.0000 - 1.0000i -0.0000 + 1.0000i -0.0000 - 1.0000i

结果分析：

从运行的结果中可以看出在虚轴上有两对重根i ±，与Routh-Hurwitz 判稳的结果相同。 ⑶ Routh-Hurwitz 判稳

系统的传递函数为()24

50351024

112

34++++=+=

s s s s s G sys ，得到特征方程024*********=++++s s s s 。

24

42243050

1024

351

1234s s s s s

劳斯表的第一列都是正数，所以系统稳定，特征方程的根都在s 平面的左边。 源程序代码：

numg=[24];deng=[1 10 35 50 0];sysg=tf(numg,deng); sys=feedback(sysg,1); pole(sys)

运行结果

-4.0000 -3.0000 -2.0000 -1.0000

结果分析：

运行结果特征方程的根都在s 平面的左边，系统稳定，与劳斯判定一样的结果。

实验2-6系统相对稳定性分析

Routh-Hurwitz 稳定判据还可用来确定使系统稳定的参数取值范围。例如，对于特征方程为

()04223=+++=K s s s s q

的系统，由 Routh-Hurwitz 稳定判据可以确定，要使之稳定，必须满足 0 < K <8 。这个结果可通过Matlab 以图解的方式加以验证。 本项实验内容包括：

1）对下列系统的特征方程运用 Routh-Hurwitz 稳定判据确定系统稳定时参数 K 的取值范围，并编写 Matlab 仿真程序加以验证： ()0234=++++=K s s s s s q

2）对图 2.6所示系统运用 Routh-Hurwitz 稳定判据决定系统稳定时参数 K 、p 和z 的取值范围，假设 p 和 z 均大于零。试运行 Matlab 仿真程序段加以验证，并在程序运行后生成的图形中标注说明 K 、p 和 z 的值域。

图2.6

⑴ Routh-Hurwitz 判稳：

K

s K s K

s s 0

1

238241-

劳斯表中的第一列都要大于0，系统稳定，即： ??

?>>0

-8K K 所以 0

源程序代码：

K=[0:0.5:20];

for i=1:length(K) q=[1 2 4 K(i)]; p(:,i)=roots(q); end

plot(real(p),imag(p),'x'),grid

xlabel('Real axis'),ylabel('imaginary axis')

运行结果：

-3

-2.5-2

-1.5-1-0.500.5

-3-2

-1

1

2

3

Real axis

i m a g i n a r y a x i s

结果分析：

从图中看出过了极点j 2±系统就不问定了，将s=2j 带入

()04223=+++=K s s s s q 中得k=8.

Routh-Hurwitz 判稳

K

s K s K

s s K

s 0

123411

11ε

εε-ε为为很小大于0的数

系统稳定有劳斯表的第一列都大于0，若K>0,那么首列变号两次;若K<0,首列变号一次；若K=0，系统临界稳定。 源程序代码：

K=[0:0.5:20];

for i=1:length(K) q=[1 1 1 1 K(i)]; p(:,i)=roots(q); end

plot(real(p),imag(p),'x'),grid

xlabel('Real axis'),ylabel('imaginary axis')

程序运行结果:

-2

-1.5-1

-0.500.51 1.5

-2-1.5-1-0.500.511.5

2Real axis

i m a g i n a r y a x i s

结果分析:

从运行的结果中看出无论K 取什么值系统都不稳定，与上面劳斯表分析的一样。 ⑵劳斯判稳： 系统的特征方程:

()()0123=+-+-+kz s p k s p s 系统是一三阶系统则有：

()()????

??

?>-->>->-kz

p k p kz p k p 10001

????

?

????

>-->>>p p z k z p k )11(0

1 源程序代码：

[p,z]=meshgrid(1.2:0.2:10,0.1:.2:10);

k=p.*(p-1)./(p-1-z);

rc=size(k);r=rc(1);c=rc(2); for i=1:r for j=1:c

if abs(z(i,j)-p(i,j)+1)<1.0e-0.3 k(i,j)=0; end

if k(i,j)<0; k(i,j)=0; end end end

M=[-50 25]; mesh(k,M)

运行结果：

20

40

60

10

20

30

40

50

020*********

结果分析：

运行的结果与理论的分析一致。

思考题

（1）关于标准二阶系统的脉冲响应，如果ζ值不变，而ωn 发生变化，响应曲线会有何变化?请用实验数据或图表加以说明。

源程序代码：

t=[0:0.1:12]; zeta=[1.2];

wn1=0.5;num1=[wn1*wn1]; den1=[1 2*zeta*wn1 wn1*wn1]; sys1=tf(num1,den1);

wn2=0.75;num2=[wn2*wn2]; den2=[1 2*zeta*wn2 wn2*wn2]; sys2=tf(num2,den2);

wn3=1;num3=[wn3*wn3]; den3=[1 2*zeta*wn3 wn3*wn3]; sys3=tf(num3,den3);

wn4=1.25;num4=[wn4*wn4]; den4=[1 2*zeta*wn4 wn4*wn4]; sys4=tf(num4,den4);

wn5=1.5;num5=[wn5*wn5]; den5=[1 2*zeta*wn5 wn5*wn5]; sys5=tf(num5,den5);

%

[y1,T1]=impulse(sys1,t); [y2,T2]=impulse(sys2,t);

[y3,T3]=impulse(sys3,t); [y4,T4]=impulse(sys4,t);

[y5,T5]=impulse(sys5,t);

%

plot(T1,y1,T2,y2,T3,y3,T4,y4,T5,y5)

xlabel('\omega_n t'), ylabel('y(t)')

title('zata=0.707 wn= 0.5, 0.75, 1, 1.25, 1.5'), grid

运行结果：

结果分析：

当 不变时，

w增大，标准二阶系统脉冲响应的上升时间减小，

n

反应变快，峰值变大，调节时间变短。

（2）某闭环控制系统如图 2.7 所示。试确定参数 a 和 K 的值使该系统稳定，并且系统对于斜坡输入的稳态误差不大于输入幅值的25%。（要求编制并运行 Matlab 程序加以说明）

图2.7闭环系统

解：

系统的开环传递函数为：

()()()()()?

?

?

??+??? ??++?

?? ??+=++++1511211101521s s s s a s aK s s s s a s K 可得开环传递函数的系数为

10aK ，那么稳定误差为aK

10

,设斜波输入为(0-1),则有: aK

10

<0.25 所以aK>40.

系统的特征方程为：

()010178234=+++++aK s K s s s 劳斯表：

()aK

s K a K s aK

k

s K s aK

s 0

2

12341260

641168

12610817

1

+-+--+

由劳斯表的首列为正有：

()??

?

??>+-+->>-01260641160

01262K a K aK K 编程序：

syms k

a1=-k/64+116/64+1260/64/k; ezplot(a1) hold on a2=40/k; ezplot(a2)

结果：