MATLAB仿真课程设计报告

北华大学

《MATLAB仿真》课程设计

姓名：

班级学号：

实习日期：

辅导教师：

前言

科学技术的发展使的各种系统的建模与仿真变得日益复杂起来。如何快速有效的构建系统并进行系统仿真，已经成为各领域学者急需解决的核心问题。特别是近几十年来随着计算机技术的迅猛发展，数字仿真技术在各个领域都得到了广泛的应用与发展。而MATLAB作为当前国际控制界最流行的面向工程和科学计算的高级语言，能够设计出功能强大、界面优美、稳定可靠的高质量程序，而且编程效率和计算效率极高。MATLAB环境下的Simulink是当前众多仿真软件中功能最强大、最优秀、最容易使用的一个系统建模、仿真和分析的动态仿真环境集成工具箱，并且在各个领域都得到了广泛的应用。

本次课程设计主要是对磁盘驱动读取系统校正部分的设计，运用自动控制理论中的分析方法，利用MATLAB对未校正的系统进行时域和频域的分析，分析各项指标是否符合设计目标,若有不符合的，根据自动控制理论中的校正方法，对系统进行校正，直到校正后系统满足设计目标为止。我组课程设计题目磁盘驱动读取系统的开环传递函数为是设计一个校正装置，使校正后系统的动态过程超调量δ%≤7%，调节时间ts≤1s。

电锅炉的温度控制系统由于存在非线性、滞后性以及时变性等特点,常规的PID控制器很难达到较好的控制效果。考虑到模糊控制能对复杂的非线性、时变系统进行很好的控制, 但无法消除静态误差的特点, 本设计将模糊控制和常规的PI D控制相结合, 提出一种模糊自适应PID控制器的新方法。并对电锅炉温度控制系统进行了抗扰动的仿真试验, 结果表明, 和常规的PI D控制器及模糊PI D复合控制器相比,模糊自适应PI D控制改善了系统的动态性能和鲁棒性, 达到了较好的控制效果。

目录

1未校正前系统的分析 (3)

1.1时域分析 (3)

1.2根轨迹分析 (3)

1.3频域分析 (5)

2串联校正及校正后系统分析 (6)

2.1校正原理 (6)

2.2校正过程 (6)

2.2.1选择增益系数 (6)

2.2.2校正环节 (6)

2.3校正后系统验证系统性能 (7)

2.3.1时域分析 (7)

2.3.2频域分析 (8)

3校正装置的实现 (9)

3.1无源校正装置 (9)

3.2有源校正装置 (10)

4对校正后的系统利用连续系统按环节离散化的数字进行仿真 (10)

4.1连续系统按环节离散化仿真原理 (10)

4.2Matlab实现仿真过程 (11)

5实习总结 (14)

6附录 (15)

7参考文献 (19)

1未校正前系统的分析

衡量一个系统的好坏主要是通过性能指标，而其中最主要的分析方法是时域分析、根轨迹和频域分析。

1.1 时域分析

时域分析法是根据系统的微分方程， 以拉氏变换为工具，直接解出控制系统的时间响应，根据响应表达式及响应曲线来分析系统的稳定性、快速性、准确性等。

我们的题目是磁盘驱动读取系统的开环传递函数为Gk(s)=1/s(s+20)(s+10)，对这个函数进行展开得到Gk(s)=1/s^3+30s^2+200s 。运用MATLAB 的step()函数对此系统进行仿真，得到系统单位阶跃响应曲线如1-1所示。

未校正系统阶跃响应

Time (sec)

A m p l i t u d e

0500

10001500

0.2

0.4

0.6

0.8

1

System: sys

Settling Time (sec): 782

System: sys

Time (sec): 1.27e+003Amplitude: 0.998

图1-1 单位阶跃响应曲线

M 文件如下:

num0=1;

den0=conv([1,0],conv([1,20],[1,10])); [numb,denb]=cloop(num0,den0); step(numb,denb);

title('未校正系统阶跃响应');

由以上分析，得到系统的时域指标δ%=0%,ts=782s 。从调节时间上看，远远地不符合设计要求。

1.2 根轨迹分析

根轨迹法是根据反馈控制系统开环和闭环传递函数之间的关系提出的一种由开环传递函数求闭环特征根的简便方法，它是分析和设计线性定常控制系统的图解方

法，使用十分简便，在工程上获得了广泛的应用。

运用根轨迹分析法，可以对系统的稳定性进行分析而这又为系统的校正提供依据。

在调用rlocus()函数之后，调用根轨迹增益函数rlocfind()，可得出系统的根轨迹曲线如1-2所示。

未校正前系统根轨迹

Real Axis

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

-80

-60

-40

-20

20

40

-50

-40-30-20

-100102030

40500.2

0.42

0.58

0.72

0.83

0.91

0.960.990.2

0.420.580.72

0.830.91

0.960.99

20

406080

图1-2 根轨迹曲线

MATLAB 的程序为； num0=1;

den0=[1,30,200,0]; axis equal;

rlocus(num0,den0);

[K,poles]=rlocfind(num0,den0); title('未校正前系统根轨迹'); 执行以上程序，并移动鼠标到根轨迹与虚轴的交点处单击鼠标左键后可得到如下结果：

selected_point = 0 +14.0727i

K =5.9413e+003 poles =-29.9464

-0.0268 +14.0853i

-0.0268 -14.0853i

由此可见，根轨迹与虚轴交点处的增益K=5941.3这说明，K<1.6时系统稳定；

当K>941.6时系统不稳定。利用rlocfind()函数可找出根轨迹从实轴上的分离点处的增益K=384，说明当0

1.3 频域分析

频域分析法是应用频域特性分析线性系统的方法。它是以传递函数为基础的一种图解分析法，对与高阶系统的分析非常实用，它同时也适用于系统的设计。而此次的系统为三阶系统，所以对系统的频域分析就显得尤为重要。

系统的bode 图可用函数bode(num,den)生成，而频域指标幅值裕量、相位裕量、相位穿越频率和幅值穿越频率通过[Gm,Pm,Wcg,Wcp]=margin(num,den)获得。

所编写的MATLAB 程序为： num0=1;

den0=[1,30,200,0]; g = tf(num0,den0); bode(g);

[Gm,Pm,Wcg,Wcp]=margin(g); title('未校正前系统伯德图');

执行后得如下数据及如图1-3所示的bode 图

Gm =

6.0000e+003 Pm =

89.9570 Wcg =

14.1421 Wcp =

0.0050

-200-150-100-500

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

)10

-3

10

-2

10

-1

10

10

1

10

2

10

3

-270

-225-180-135

-90P h a s e (d e g )

未校正前系统伯德图

Frequency (rad/sec)

图1-3 系统bode图

从图中可以看出相位裕量γ=89.9570度，即裕量为正，开环系统稳定，还有较好的动态性能指标。

2串联校正及校正后系统分析

自动控制系统是由被控对象和控制器两大部分组成的，当被控对象被确定后，对系统的设计实际上就是对控制器的设计，这就称为对控制系统的校正。

由于控制系统不满足控制质量的要求，需根据系统预先给定的性能指标重新设计一个满足性能要求的控制系统，具体任务是选择校正方式，确定校正装置的类型以及计算出具体参数等。

2.1校正原理

校正装置是为了改善系统控制性能而人为的引入的控制部分根据校正装置在控制系统中的位置不同，可分为串联校正和并联校正两类基本形式。串联校正分为超前校正、滞后校正、滞后-超前校正三种结合目前的分析方法对串联校正的选择方法又可分为bode图法、根轨迹法和频率特性法。根据他们各自的特点和本次设计的要求最后选择了频率特性法来进行校正装置的设计。

结合前面对系统性能指标的分析，我们得出系统的调节时间过长，主要是增益过小导致的，所以进行校正前的首要任务是确定系统的增益系数。确定了增益系数之后就可以借助于频域特性法来设计校正装置。

2.2校正过程

2.2.1选择增益系数

根据增益系数对系统稳定性影响的特点以及前面对系统根轨迹分析得出的系统稳定时增益的范围（0

2.2.2校正环节

选择了K=300之后，观察到系统的动态响应的快速性仍有待提高，所以就选择了串联超前校正。

接下来我们根据串联校正编写了MATLAB程序:

num0=300;

den0=conv([1,0],conv([1,20],[1,10]));

g=tf(num0,den0);

[mag1,phase1]=bode(g);

[Gm,Pm,Wcg,Wcp]=margin(g);

r=62;w=logspace(-3,1);%给定期望系统的相位域量

for epsilon=5:15

r0=(-180+r+epsilon);

[i1,ii]=min(abs(phase1-r0));

wc=w(ii);

alpha=mag1(ii);

T=5/wc;

numc=[T,1];denc=[alpha*T,1];

[num,den]=series(num0,den0,numc,denc);

[Gm,Pm,Wcg,Wcp]=margin(g);

if(Pm>=r);break;end;

end

printsys(numc,denc);printsys(num,den);

g=tf(numc,denc);

figure(1);

bode(g) %画出校正环节伯德图

title('校正环节伯德图');

%[mag2,phase2]=bode(g);

g=tf(num0,den0);

figure(3);

bode(g); %画出有增益未校正系统伯德图

title('有增益未校正系统bode图');

g=tf(num,den);

figure(2);

bode(g); %画出校正后环节伯德图

title('有增益已校正系统bode图');

%[mag,phase]=bode(g);

%未校正前

figure(3);

subplot(2,1,1);

[numb,denb]=cloop(num0,den0);

step(numb,denb); %有增益未校正阶跃响应曲线title('未校正加入增益前系统单位阶跃响应曲线');

%校正后

subplot(2,1,2);

[numb,denb]=cloop(num,den);

step(numb,denb); %有增益已校正阶跃响应曲线title('校正后系统单位阶跃响应曲线');

经过多次调试得到校正环节的传递函数Gc=（359.8428 s + 1）/ （201.3615 s + 1）。从而得到Gk（s）=（107952.851 s + 300）/（201.3615 s^4 + 6041.8461 s^3 + 40302.307 s^2 + 200 s）。

2.3校正后系统验证系统性能

2.3.1时域分析

00.51 1.52 2.53

0.20.40.60.81未校正加入增益前系统单位阶跃响应曲线

Time (sec)

A m p l i t u d e

校正后系统单位阶跃响应曲线

Time (sec)

A m p l i t u d e

0.5

1 1.5

00.20.40.60.81

1.21.4System: sys P eak amplitude: 1.01

Overshoot (%): 1.22At time (sec): 1.14

System: sys

Settling Time (sec): 0.83

图2-3-1 校正后-校正前系统动态响应

所以由仿真图得到超调量δ=1.22%，调节时间ts=0.83s ，符合期望的系统性能。

2.3.2频域分析

-150-100

-50

50

100M a g n i t u d e (d

B )10

-4

10

-3

10

-2

10

-1

10

10

1

10

2

10

3

-270

-225-180-135-90

-45P h a s e (d e g )

有增益已校正系统bode 图

Frequency (rad/sec)

图2-3-2 校正后的系统bode 图

所以由仿真结果得到Gm = 11.1953，Pm =68.2751，Wcg =14.1445，Wcp =2.5746，所以得到相位裕量γ=68.2751度，即裕量为正，开环系统稳定，还有很好的动态性能指标。

3校正装置的实现

校正装置按其使用的动力源和信号性质的不同，可分为电气型、气动型、液压型等，根据此次系统设计的要求采用电气型校正装置。超前校正装置，既可采用RC 无源网络组成，又可用由运算放大器加入适当电路的有源网络组成，下面分别介绍无源超前网络和有源超前网络在本次校正环节中的具体应用。

3.1无源校正装置

由MATLAB 仿真程序得出校正环节传递函数： Gc=（359.8428 s + 1）/ （201.3615 s + 1）(1) 根据无源超前网络传递函数Gc(s)=α（Ts+1）/(αTs+1)，α=R2/(R1+R2);T=R1C(2) 有（1）(2)式对比计算得出，R1=100K,R2=150K,C=3.6Uf

网络衰减α=3/5，放大器的放大系数就的增大1/α倍，既α’=5/3, 所以K=k*α’=300*5/3=500.

所以画出图3-1-1无源超前校正网络。

图3-1无源超前校正网络

3.2有源校正装置

由MATLAB仿真程序得出校正环节传递函数：Gc=（359.8428 s + 1）/ （201.3615 s + 1）,可得到如下曲线图3-2：

图3-2有源超前校正网络

4对校正后的系统利用连续系统按环节离散化的数字进行仿真

4.1连续系统按环节离散化仿真原理

连续系统数学模型的离散化是通过数值积分法实现的，尽管面向结构图的仿真方法是按环节给定参数，但是在仿真计算时还是按整个系统进行离散化这就不便于引进非线性环节系统的仿真所以就要进行连续系统按环节离散化的仿真过程。

首先是典型环节的离散系数及其差分方程的确定。一般情况下，典型环节的仿真模型归纳为一个统一的公式；

z[(k+1)T]=Ez(kT)+Fu(kT)+Gu[(k-1)T]

x[(k+1)T]=Hz[(k+1)T]+Lu(kT)+Qu[(k-1)T]

式中，E，F，G，H，L，Q是差分方程的系数，它们的数值根据典型环节系数a,b的

不同情况可由表4-1确定。

a ≠0,b=0 a=0,

b ≠0 a

≠0,b ≠0

E F G H L Q

0 0 0 0

(c+d/T)/a -d/(aT)

1

cT/(2b) F 1 d/b 0

Exp[-(a/b)T]

(d/b-c/a)[(1-E)*b/(aT)-1] (d/b-c/a)[1+(E-1)(1+b/(aT))] 1 d/b 0

表4—1 E , F, G , H, L, Q 的系数

4.2Matlab 实现仿真过程

图4-2-1 系统结构图

根据系统结构图可得，

u1

0 0 0 -1

x1 1 x1 u2 = 1 0 0 0 x2 + 0 , y=[0 0 0 1] x2 u3 0 1 0 0 x3 0 x3 u4

0 0 1 0

x4 0 x4

0 0 0 -1

1

则有 W= 1 0 0 0 ， W0= 0 , Wc =[0 0 0 1]

0 1 0 0 0

0 0 1 0 0

a1 b1 c1 d1 FZ1 s1 1 201.3615 1 359.8428 0 0 而P= a2 b2 c2 d2 FZ2 s2 = 0 1 1 0 0 0 a3 b3 c3 d3 FZ3 s3 10 1 300 0 0 0 a4 b4 c4 d4 FZ4 s4 20 1 1 0 0 0

仿真程序如下： R=1;

P=[ 1 201.3615 1 359.8428 0 0; %A 0 1 1 0 0 0; %B 10 1 300 0 0 0; %C

20 1 1 0 0 0]; %D

W=[ 0 0 0 -1;

1 0 0 0;

0 1 0 0;

0 0 1 0];

W0=[1;0;0;0];

Wc=[0 0 0 1];

Tf=input('仿真时间Tf>>');

T=input(' 步长时间T>>');

A=P(:,1);

B=P(:,2);

C=P(:,3);

D=P(:,4);

FZ=P(:,5);

S=P(:,6);

n=length(A);

for i=1:n

if(A(i)~=0)

if(B(i)==0)

E(i)=0;F(i)=0;G(i)=0;H(i)=0;

L(i)=(C(i)+D(i)/T)/A(i);Q(i)=-D(i)/(A(i)*T);

else

E(i)=exp(-A(i)*T/B(i));

F(i)=(D(i)/B(i)-C(i)/A(i))*((1-E(i))*B(i)/(A(i)*T)-1);

G(i)=(D(i)/B(i)-C(i)/A(i))*(1+(E(i)-1)*(1+B(i)/(A(i)*T)));

H(i)=1;L(i)=D(i)/B(i);Q(i)=0;

end

else

if(B(i)~=0)

E(i)=1;F(i)=0.5*C(i)*T/B(i);G(i)=F(i);

E(i)=1;F(i)=0.5*C(i)*T/B(i);G(i)=F(i);

H(i)=1;L(i)=D(i)/B(i);Q(i)=0;

else

disp('A(i)=B(i)=0');

end

end

end

d

end

x=zeros(length(A),1);x0=x;z=x;

u=zeros(length(A),1);u0=u;

y=zeros(length(Wc(:,1)),1);t=0;

for j=1:Tf/T

u1=u;u=W*x+W0*R;

x1=x; for i=1:n

z(i)=E(i)*z(i)+F(i)*u(i)+G(i)*u1(i); x(i)=H(i)*z(i)+L(i)*u(i)+Q(i)*u1(i); end

y=[y,Wc*x];t=[t,t(j)+T];

if((abs(y-R)/R)<=0.02) data_adjust=j; end end

plot(t,y);

[pos,tr,ts,tp]=stepchar(t,y,0.02);

disp(['最大超调量M=' num2str(pos) '%']); disp([' 调节时间ts=' num2str(ts) 's']); grid on; hold on;

取仿真时间：Tf=5 计算步长：T=0.001

在matlab 环境下执行以上程序可得到如图4-3所示的仿真曲线和如下结果。 最大超调量δ%=1.5471%，调节时间ts=0.813s 。

0123456

0.2

0.4

0.6

0.8

1

1.2

1.4

数字仿真曲线

图4-2-2所示的仿真曲线5实习总结

6附录

Code1:未校正前系统时域响应

num0=1;

den0=conv([1,0],conv([1,20],[1,10]));

[numb,denb]=cloop(num0,den0);

step(numb,denb);

title('未校正系统阶跃响应');

Code2:未校正前系统根轨迹

num0=1;

den0=[1,30,200,0];

axis equal;

rlocus(num0,den0);

[K,poles]=rlocfind(num0,den0);

title('未校正前系统根轨迹');

Code3:未校正前系统频域分析

num0=1;

den0=[1,30,200,0];

g = tf(num0,den0);

bode(g);

[Gm,Pm,Wcg,Wcp]=margin(g);

title('未校正前系统伯德图');

Code4:校正环节及其校正后系统相应指标分析num0=300;

den0=conv([1,0],conv([1,20],[1,10]));

g=tf(num0,den0);

[mag1,phase1]=bode(g);

[Gm,Pm,Wcg,Wcp]=margin(g);

r=62;w=logspace(-3,1);%给定期望系统的相位域量for epsilon=5:15

r0=(-180+r+epsilon);

[i1,ii]=min(abs(phase1-r0));

wc=w(ii);

alpha=mag1(ii);

T=5/wc;

numc=[T,1];denc=[alpha*T,1];

[num,den]=series(num0,den0,numc,denc);

[Gm,Pm,Wcg,Wcp]=margin(g);

if(Pm>=r);break;end;

end

printsys(numc,denc);printsys(num,den);

g=tf(numc,denc);

figure(1);

bode(g) %画出校正环节伯德图

title('校正环节伯德图');

%[mag2,phase2]=bode(g);

g=tf(num0,den0);

figure(3);

bode(g); %画出有增益未校正系统伯德图title('有增益未校正系统bode图');

g=tf(num,den);

figure(2);

bode(g); %画出校正后环节伯德图[Gm,Pm,Wcg,Wcp]=margin(g);

title('有增益已校正系统bode图');

%[mag,phase]=bode(g);

%未校正前

figure(3);

subplot(2,1,1);

[numb,denb]=cloop(num0,den0);

step(numb,denb); %有增益未校正阶跃响应曲线title('未校正加入增益前系统单位阶跃响应曲线');

%校正后

subplot(2,1,2);

[numb,denb]=cloop(num,den);

step(numb,denb); %有增益已校正阶跃响应曲线title('校正后系统单位阶跃响应曲线');

Code5:对校正后系统进行离散化数字仿真

R=1;

P=[ 1 201.3615 1 359.8428 0 0; %A

0 1 1 0 0 0; %B

10 1 300 0 0 0; %C

20 1 1 0 0 0]; %D

W=[ 0 0 0 -1;

1 0 0 0;

0 1 0 0;

0 0 1 0];

W0=[1;0;0;0];

Wc=[0 0 0 1];

Tf=input('仿真时间Tf>>');

T=input(' 步长时间T>>');

A=P(:,1);

B=P(:,2);

C=P(:,3);

D=P(:,4);

FZ=P(:,5);

S=P(:,6);

n=length(A);

for i=1:n

if(A(i)~=0)

if(B(i)==0)

E(i)=0;F(i)=0;G(i)=0;H(i)=0;

L(i)=(C(i)+D(i)/T)/A(i);Q(i)=-D(i)/(A(i)*T);

else

E(i)=exp(-A(i)*T/B(i));

F(i)=(D(i)/B(i)-C(i)/A(i))*((1-E(i))*B(i)/(A(i)*T)-1);

G(i)=(D(i)/B(i)-C(i)/A(i))*(1+(E(i)-1)*(1+B(i)/(A(i)*T)));

H(i)=1;L(i)=D(i)/B(i);Q(i)=0;

end

else

if(B(i)~=0)

E(i)=1;F(i)=0.5*C(i)*T/B(i);G(i)=F(i);

H(i)=1;L(i)=D(i)/B(i);Q(i)=0;

else

disp('A(i)=B(i)=0');

end

end

end

x=zeros(length(A),1);x0=x;z=x;

u=zeros(length(A),1);u0=u;

y=zeros(length(Wc(:,1)),1);t=0;

for j=1:Tf/T

u1=u;u=W*x+W0*R;

x1=x;

for i=1:n

z(i)=E(i)*z(i)+F(i)*u(i)+G(i)*u1(i);

x(i)=H(i)*z(i)+L(i)*u(i)+Q(i)*u1(i);

end

y=[y,Wc*x];t=[t,t(j)+T];

if((abs(y-R)/R)<=0.02)

data_adjust=j;

end

end

plot(t,y);

[pos,tr,ts,tp]=stepchar(t,y,0.02);

disp(['最大超调量M=' num2str(pos) '%']); disp([' 调节时间ts=' num2str(ts) 's']); grid on;

7参考文献

1 李国勇，谢克明，杨丽娟。计算机仿真技术与CAD——基于MATLAB的控制系统（第二版）

2高国燊，余文烋，彭康拥，陈好来。自动控制原理（第二版）。华南理工大学出版社2005

3 翁思义. 自动控制系统计算机仿真与辅助设计. 西安：西安交通大学出版社，1987 4刘豹主编. 现代控制理论（第二版. 北京：机械工业出版社，1988）

5张培强.MATLAB语言.合肥：中国科技大学出版社，1995