控制系统仿真设计报告-加热炉出口温度仿真设计

控制系统仿真设计报告-加热炉出口温度仿真设计 本文关键词：仿真，加热炉，设计，控制系统，温度

控制系统仿真设计报告-加热炉出口温度仿真设计 本文简介：内蒙古科技大学控制系统仿真课程设计报告题目：加热炉出口温度仿真设计学生姓名：学号：专业：班级：指导教师：目录第1章概述1第2章总体方案设计22.1设计思路22.1.1单回路控制22.1.2串级控制系统2第3章Simulink建模53.1单回路控制整定53.2串级控制参数整定7第4章Simulink仿

控制系统仿真设计报告-加热炉出口温度仿真设计 本文内容：

内蒙古科技大学

控制系统仿真课程设计报告

题

目：加热炉出口温度仿真设计

学生姓名：

学

号：

专

业：

班

级：

指导教师：

目

录

第1章

概述1

第2章

总体方案设计2

2.1

设计思路2

2.1.1

单回路控制2

2.1.2

串级控制系统2

第3章

Simulink建模5

3.1

单回路控制整定5

3.2

串级控制参数整定7

第4章

Simulink仿真与优化设计10

第5章

总结13

参考文献14

第1章

概述

图1-1所示为某工业生产中的加热炉，其任务是将被加热物料加热到一定温度,然后送到下道工序进行加工。加热炉工艺过程为：加热物料流过排列炉膛四周的管道后,加热到炉出口工艺所要求的温度。在加热用的燃料油管道上装有一个调节阀,用以控制燃料油流量,以达到控制出口温度的目的。

图

1-1

加热炉出口温度系统

但是，由于加热炉时间常数大，而且扰动的因素多，单回路反馈控制系统不能满足工艺对加热炉出口温度的要求。为了提高控制质量，采用串级控制系统，运用副回路的快速作用，有效地提高控制质量，满足生产要求。

第2章

总体方案设计

2.1

设计思路

2.1.1

单回路控制

在这里先可以尝试采用使用单闭环实现。加热炉出口温度单回路反馈控制系统结构框图如图2-1。

图2-1

单回路控制系统框图

控制器使用PI控制规则，使用衰减曲线法对控制器参数进行整定。

2.1.2

串级控制系统

由于主对象时间常数较大，直接采用串级控制。

在这里选择之后较小的炉膛温度为副变量，炉出口温度为主变量。则该系统的结构图如图2-2。

图2-2

串级控制系统图

因此可以得到系统的方框图，如图2-3。

图2-3

串级控制系统框图

主控制器采用PI控制规则，副控制器直接采用P控制器规则。使用一步法直接确定副控制器的放大倍数，再使用单回路控制系统的整定方法来整定主控制器的参数。

第3章

Simulink建模

3.1

单回路控制整定

使用衰减曲线法来整定PI调节器，将积分时间T设为最大值，比例增益K设为一个较小值。给系统输入阶跃信号，调整K使系统的响应出现4:1的衰减过程。并记录下此时的比例环节增益K和震荡周期T。

Simulink仿真环节如图3-1：

图3-1

单回路P控制Simulink图

经过尝试，当K=1.254时，此时T=16.6s正好是4：1衰减。如图3-2所示：

图3-2

单回路P控制4：1响应曲线

由图3-2可以明显看出系统存在一定的余差，需要在控制器中加入积分来消除系统的稳态偏差。按表3-1的经验公式进行调节。

表3-1

衰减曲线法整定计算公式

由得到PI调节器的正定参数：K=1.254/1.2=1.045，T=16.6×0.5=8.30s。加入PI缓解后系统Simulink如图3-3所示：

图3-3

单闭环PI控制Simulink仿真图

经过略微调整K

和T使系统满足4：1的衰减变化。最终使K=0.909，T=13.65s时于是得到如图3-4的仿真曲线。

图3-4

单闭环PI控制4：1响应曲线

曲线较好的呈现4：1的衰减比，经过两个震荡周期便进入了稳态。由于积分的作用，系统已消除了静差。

3.2

串级控制参数整定

由一步法参照副控制器参数经验设置值表来整定副控制器。如表3-2所示。

表3-2

副控制器参数经验设置值

这里副变量是炉膛的温度，则由上表直接整定副控制器参数为K=3，K=5

按单回路控制系统整定方法直接整定主控制器参数。将内环直接看作主环内的一个环节G，因为此时副回路可以看作一个理想的随动系统。使用衰减曲线法4：1，这时候同样按照单闭环的方法来整定外环，Simulink如图3-5示：

图3-5

串级控制主P控制Simulink仿真图

当K=19.4时外环出现4：1，如图3-6所示：

图3-6

串级控制主P控制器4：1响应曲线

参考表3-1中衰减曲线法的经验公式，则有K=19.4/1.2=16.17，T=16.6×0.5=8.30s。加入PI环节后系统Simulink如图3-7所示。

图3-7

串级主PI控制Simulink仿真图

经过多次尝试调整K和T使系统满足4：1的衰减变化。最终使K=13，T=20s时于是得到如图3-8的仿真曲线：

图3-8

串级主PI控制4：1响应曲线

曲线较好的呈现4：1的衰减比，经过两个震荡周期便进入了稳态。由于积分的作用，系统已消除了静差。

第4章

Simulink仿真与优化设计

单闭环控制器和串级控制器参数均已调好。如表4-1所示，串级控制的优势十分明显。这里分别从衰减率、静态偏差、系统的工作频率、对干扰的抑制作用和副对象干扰的抑制作用这几个方面进行计算和分析。

由表中所见，在同样为4:1衰减比的情况，两者均因为PI控制器的作用，使得系统无静态偏差。但串级系统由于内环的存在，提高了工作频率，由0.0499提高到0.0873，几乎两倍。副回路存在使得系统响应变化，更早的达到峰值。在加入干扰后，串级控制的自适应能力便脱颖而出，由图4-1的响应曲线所见，单回路虽然能克服干扰，但最大偏差明显大于串级控制。当副回路中引入干扰时，串级控制使单回路的动态偏差由0.112降到0.0032；即使是进入主回路的干扰(t)=0.1，最大动态偏差也由单回路的0.023降到0.0102。可见，串级控制系统对扰动有一定的抑制作用，尤其是进入到副回路内的扰动。

表4-1

串级控制效果与单回路控制比较的结果

串级控制和单回路控制的仿真对比结果以图形反映如下，明显看出两者差别。

给定单位阶跃变化时的响应曲线

作用时的响应曲线

(t)=0.1作用时的响应曲线

副对象放大系数变为2时响应曲线

最后一图显示当副对象由时，串级控制的响应曲线基本不变，而单回路出现次数较多的震荡效果明显变坏，可见串级控制系统有一定的自适应能力。

第5章

总结

这次课程设计主要是学习整定串级系统的控制参数。这次课程设计在参数整定上遇到不少问题，最终通过耐心地调整以及同学帮助，终于得到了较好的控制效果。通过本次课程设计，可以加深对过程控制理论知识的理解，同时更加熟练的使用matlab以及Simulink进行控制系统的仿真实验。将所学的理论与方法应用于实际之中，提高了自己理论联系实际的能力，也丰富了自己的阅历。

参考文献

[1]

方康玲

过程控制与集散系统(第一版)[M]，电子工业出版社，2009.1

[2]

吴怀宇

自动控制原理(第一版)[M]，华中科技大学出版社，2007.9

[3]

郑阿奇

MATLAB实用教程(第二版)[M]，电子工业出版社，2007.8

篇2：微波技术与天线仿真实验报告

微波技术与天线仿真实验报告 本文关键词：天线，微波，仿真，实验，报告

微波技术与天线仿真实验报告 本文简介：微波技术与天线仿真实验报告姓名：柳立志学号：110250208ANSOFTHFSS软件的使用与魔T的仿真一、实验内容1、下载并且安装ANSOFTHFSS软件10.0版本；2、学习使用该软件；3、仿真魔T；4、写出仿真使用后的报告。二、软件介绍与使用AnsoftHFSS软件设置步骤如下：1.打开Ans

微波技术与天线仿真实验报告 本文内容：

微波技术与天线

仿真实验报告

姓名：柳立志

学号：110250208

ANSOFT

HFSS软件的使用与魔T的仿真

一、实验内容

1、下载并且安装ANSOFT

HFSS软件10.0版本；

2、学习使用该软件；

3、仿真魔T；

4、写出仿真使用后的报告。

二、软件介绍与使用

Ansoft

HFSS软件设置步骤如下：

1.打开Ansoft

HFSS

V10.0，点击TOLLS栏进行软件仿真设置；

2.点击Options中的HFSS

Options，在General选项卡中将“Use

wizards

for

data

input

when

creating

new

boundary“和“Duplicate

boundaries

with

geometry“前的复选框打钩，点击OK；

3.点击Options中的3D

Modeler

Options，在Operation选项卡中将“Automatically

cover

closed

ployline“前的复选框打钩，在Drawing选项卡中将“Edit

property

new

primitives“前的复选框打钩，点击OK；

三、仿真魔T

1、创建三维模型

1.设置模型单位

点击3D

Modeler，选择Units，设置单位为“mm“，点击OK；

2.设置默认材料

在3D

Modeler

Materials

Tollbar中选择材料类型为vacuum；

3.设置第一个模块参数

点击Draw中Box，设置起始坐标为X:-25.0，Y:-10.0，Z:0.0，按下回车，接着输入相对坐标dX：50.0，dY：20.0，dZ：75.0，按下回车，按快捷键Ctrl+D，使模块以适当大小显示在屏幕内；

4.设置波端口激励

点击Edit栏Select中Faces，选择模块的顶面（Z=75.0），点击HFSS栏Excitations=>Assign中的Wave

Port，打开Wave

Port对话框，输入名称为p1，Next=>Next=>Finish,完成对波端口激励的设置；

5.更改选择目标

点击Edit栏Select中的Object，完成设置；

6.创建第二个模块

按下Ctrl+A选择全部可视，点击Edit栏Duplicate中Around

Axis，使第一个模块按轴线旋转复制，选择Axis为X轴，设置角度为90度，点击OK；

7.创建第三、四个模块

点击Edit栏Select中Select

Name，选择第二个模块同样绕Z轴旋转90度复制，再以第三个模块绕绕Z轴旋转90度复制获得最终模型；

8.组合所有模块

按下Ctrl+A选择全部模块，点击3D

Modeler栏Boolean中Unites，将所有模块合并，Ctrl+D使模块以适当大小显示在屏幕内，模型如下图所示；

9.边界显示设置

保存工程，点击HFSS栏中Boundary

Display(Solve

view)，选择要显示的边界，outer为背景显示（图中蓝色边界），黑、红、绿、黄分别分边界p1、p2、p3、p4的边界显示；

点击View栏Active

View

Visibility

.，将所有实体隐藏，只显示边界，如下图所示：

2、分析设置

1.创建一个分析设置

点击HFSS栏Analysis

Setup中Add

Solution

Setup.，在General选项卡中设置频率为4GHz，最大量程为5，最大增量为0.02，点击OK；

2.添加频率扫描

点击HFSS栏Analysis

Setup中Add

Sweep.，选择要设置项Setup1，点击OK进入Edit

Sweep对话框，选择扫描类型为Fast，设置频率类型为Linear

Setup，开始频率为3.4GHz，停止频率为4.0GHz，计数1001次（即设置Size约为600KHz），在Save

Fields前的复选框打钩，点击OK；

3.保存工程

点击File栏Save

As，更改工程名为hfss_magic_T，点击Save；

3、分析

4.模型验证

点击HFSS栏Validation

Check进行检查验证，没有错误时点击Close；

5.分析

点击HFSS栏Analysis

All，进行分析；

4、生成报告

通过仿真得到如下两种不同结果：

5、动态仿真

最终进行动态仿真，由软件观察魔T的传输使用过程，静态图如下：

软件使用体会：我使用的是ANSOFT

HFSS14.0版本，软件操作语言为英文，它是三维电磁仿真软件。HFSS提供了一简洁直观的用户设计界面，可以仿真三维电磁场，可直接得到特征阻抗、传播常数、S参数及电磁场、辐射场、天线方向图等结果。

实验结论：波导魔T是微波、毫米波电路中的重要器件,能进行功率的分配与合成,由于其具有功率容量高、端口性能好等优点,被广泛应用于微波集成电路、电子对抗设备以及制导系统等。

篇3：系统仿真实验报告

系统仿真实验报告 本文关键词：仿真，实验，报告，系统

系统仿真实验报告 本文简介：系统仿真实验报告班级：电气工程及其自动化1301班学号：姓名：指导老师：完成时间：2015年4月19日目录实验一MATLAB中矩阵与多项式的基本运算-1-实验二MATLAB绘图命令8实验三MATLAB程序设计10实验四MATLAB的符号计算与SIMULINK的使用11实验五MATLAB在控制系统分析

系统仿真实验报告 本文内容：

系统仿真实验报告

班

级：电气工程及其自动化1301班

学

号：

姓

名：

指导老师：

完成时间：2015年4月19日

目

录

实验一

MATLAB中矩阵与多项式的基本运算-

实验二

MATLAB绘图命令8

实验三

MATLAB程序设计10

实验四

MATLAB的符号计算与SIMULINK的使用11

实验五

MATLAB在控制系统分析中的应用13

实验六

连续系统数字仿真的基本算法22

实验一

MATLAB中矩阵与多项式的基本运算

实验任务

1．了解MATLAB命令窗口和程序文件的调用。

2．熟悉如下MATLAB的基本运算：

①

矩阵的产生、数据的输入、相关元素的显示；

②

矩阵的加法、乘法、左除、右除；

③

特殊矩阵：单位矩阵、“1”矩阵、“0”矩阵、对角阵、随机矩阵的产生和运算；

④

多项式的运算：多项式求根、多项式之间的乘除。

【命令】与【运行结果】

eye(3)

ans

ones(3)

ans

zeros(3,4)

ans

rand(3,4)

ans

0.8147

0.9134

0.2785

0.9649

0.9058

0.6324

0.5469

0.1576

0.1270

0.0975

0.9575

0.9706

diag(3)

ans

A=rand(3),diag(A)

0.9572

0.1419

0.7922

0.4854

0.4218

0.9595

0.8003

0.9157

0.6557

ans

0.9572

0.4218

0.6557

A=[1

3;4

6;7

9];.

B=[1

5;7

2;4

8];.

A＼B,B/A,inv(A)*B,B*inv(A),inv(B)*A

Warning:

Matrix

singular

badly

scaled.

Results

may

inaccurate.

RCOND

1.541976e-018.

ans

1.0e+016

4.0532

-4.0532

-8.1065

8.1065

4.0532

-4.0532

Warning:

Matrix

singular

working

precision.

ans

NaN

-Inf

Inf

NaN

-Inf

Inf

NaN

-Inf

Inf

Warning:

Matrix

singular

badly

scaled.

Results

may

inaccurate.

RCOND

1.541976e-018.

ans

1.0e+016

4.0532

-4.0532

-8.1065

8.1065

4.0532

-4.0532

Warning:

Matrix

singular

badly

scaled.

Results

may

inaccurate.

RCOND

1.541976e-018.

ans

1.0e+016

0.0000

4.0532

-8.1065

4.0532

ans

3.5000

3.0000

2.5000

-2.5000

-2.0000

-1.5000

1.0000

【结果分析】

题目中的矩阵A是没有逆的，但由于计算机精度的问题，使得上面的结果比较奇怪：

首先，B/A应该和B*inv（A）都出错显示结果为无定义，但是事实上却是前者结果无定义，而后者却又结果。另外A/B和inv(A)*B

应该也是没有定义的，但还是有结果，因此，我们有必要检查运算结果，鉴于软件在精度方面不可避免的误差。

>>syms

f=3*x^5+4*x^3-5*x^2-7*x+5;

p=[3,0,4,-5,-7,5];

X=roots(p)

poly(X)

-0.3074

1.6176i

-0.3074

1.6176i

-1.0000

1.0000

0.6147

ans

1.0000

0.0000

1.3333

-1.6667

-2.3333

1.6667

A=fix(10*rand(1,3)),B=fix(20*rand(1,4))

conv(A,B)，[Q,r]=deconv(B,A)

ans

123

1.5000

2.9167

1.5833

8.0833

A=fix(10*rand(3)),B=fix(20*rand(3))

A*B，A.*B

ans

110

155

273

342

149

151

ans

135

【结果分析】：A*B是正常的规矩乘法，而A.*B是一种扩展功能，其结果Cij=Aij*Bij。

who

Your

variables

are:

ans

whos

Name

Size

Bytes

Class

Attributes

3x3

double

3x3

double

1x2

double

5x1

double

complex

ans

3x3

double

1x1

170

sym

1x6

double

1x4

double

1x1

126

sym

A=fix(10*rand(1,3)),B=fix(20*rand(3)),C=fix(10*rand(3,2)),D=[3]

diag(A),diag(B),diag(C),diag(D)

ans

disp(

lalalalaalalllallla

);

A=fix(10*rand(3,2)),B=zeros(3),C=ones(3,6)

size(A),length(A),size(B),length(B),size(C),length(C)

lalalalaalalllallla

ans

【结果分析】：disp不受”;分号”影响，size输出矩阵的行数和列数，length输出行列数较大值

实验二

MATLAB绘图命令

实验任务

熟悉MATLAB基本绘图命令，掌握如下绘图方法：

1．坐标系的选择、图形的绘制；

2．图形注解（题目、标号、说明、分格线）的加入；

3．图形线型、符号、颜色的选取。

【程序】

x=linspace(0,100,pi);

y1=x.*x;

y2=x.*x.*x;

subplot(2,2,1);plot(x,y1,-.b,x,y2,r

),legend(

y=x^{2},y=x^{3}

);

grid

off;

subplot(2,2,2);loglog(x,y1,rp-,x,y2,:b

),title(

说啥好呢

),grid

on;

subplot(2,2,3);semilogx(x,y1,b*-,x,y2,r+-

),xlabel(

),ylabel(

);

subplot(2,2,4);semilogy(x,y1,g-+,x,y2,m

),text(50,8000,y=x^{2}

),text(40,100000,y=x^{3}

)

【结果截图】

【程序】

t=0:0.01:2*pi;

r=sin(2*t).*cos(2*t);

subplot(2,2,1);polar(t,r)

subplot(2,2,2);bar(t,r,r

),axis([0,1,-0.5,0.5]),title(

bar

)

subplot(2,1,2);stairs(t,r,c

),axis([0,1,-0.5,0.5]),title(

stairs

)

【结果截图】

【程序】

[x,y,z]=peaks(100);

subplot(1,2,1),mesh(x,y,z)

subplot(1,2,2),contour(x,y,z)

【结果截图】

实验三

MATLAB程序设计

实验任务

用一个MATLAB语言编写一个程序：输入一个自然数，判断它是否是素数，如果是，输出“It

one

prime”，如果不是，输出“It

not

one

prime.”。要求通过调用子函数实现。最好能具有如下功能：①设计较好的人机对话界面，程序中含有提示性的输入输出语句。②能实现循环操作，由操作者输入相关命令来控制是否继续进行素数的判断。如果操作者希望停止这种判断，则可以退出程序。③如果所输入的自然数是一个合数，除了给出其不是素数的结论外，还应给出至少一种其因数分解形式。例：输入

6，

因为6不是素数。则程序中除了有“It

not

one

prime”的结论外，还应有：“6=2*3”的说明。

【程序】：

a=input(

请输入一个整数(end

0):

);

while(a~=0)

a==1

disp(

not

prime

composite

number.

)

elseif

isprime(a)==1

fprintf(

prime./n,a)

elseif

isprime(a)==0

f=factor(a);m=length(f);

fprintf(

not

prime,%d=%d,a,a,f(1))

for

i=2:m

fprintf(%d,f(i));

end

fprintf(

);

end

a=input(

请输入一个整数(end

0):

);

end

【运行结果】：

请输入一个整数(end

0):1

not

prime

composite

number.

请输入一个整数(end

0):2

prime.

请输入一个整数(end

0):35

not

prime,35=5*7

请输入一个整数(end

0):0

实验四

MATLAB的符号计算与SIMULINK的使用

实验任务

1．掌握MATLAB符号计算的特点和常用基本命令；

2．掌握SIMULINK的使用。

【程序1】

syms

A=[a

b;c

d];B=[1,4;2,9];

det(A),det(B),eig(A),eig(B)

【运行结果1】

ans

a*d-b*c

ans

1/2*a+1/2*d+1/2*(a^2-2*a*d+d^2+4*b*c)^(1/2)

1/2*a+1/2*d-1/2*(a^2-2*a*d+d^2+4*b*c)^(1/2)

ans

0.1010

9.8990

【程序2】

r=solve(

x^3+5*x^2-9*x+17

)

r4=vpa(r,4)

r10=vpa(r,10)

【运行结果2】

-1/3*(557+3*18849^(1/2))^(1/3)-52/3/(557+3*18849^(1/2))^(1/3)-5/3

1/6*(557+3*18849^(1/2))^(1/3)+26/3/(557+3*18849^(1/2))^(1/3)-5/3+1/2*i*3^(1/2)*(-1/3*(557+3*18849^(1/2))^(1/3)+52/3/(557+3*18849^(1/2))^(1/3))

1/6*(557+3*18849^(1/2))^(1/3)+26/3/(557+3*18849^(1/2))^(1/3)-5/3-1/2*i*3^(1/2)*(-1/3*(557+3*18849^(1/2))^(1/3)+52/3/(557+3*18849^(1/2))^(1/3))

-6.717

.858-1.339*i

.858+1.339*i

r10

-6.716750613

.858375306-1.339469138*i

.858375306+1.339469138*i

【simulink仿真】

【仿真结果】

实验五

MATLAB在控制系统分析中的应用

一、实验任务

1．掌握MATLAB在控制系统时间响应分析中的应用；

2．掌握MATLAB在系统根轨迹分析中的应用；

掌握MATLAB控制系统频率分析中的应用；

掌握MATLAB在控制系统稳定性分析中的应用

求下面系统的单位阶跃响应

【程序】

num=[4]

;

den=[1,1,4]

;step(num,den)

[y,x,t]=step(num,den)

;

tp=spline(y,t,max(y))

%计算峰值时间

ymaxmax(y)

%计算峰值

【结果】

=1.5708

ymax

=1.4419

求下面系统的单位脉冲响应：

【程序】

num=[4]

;

den=[1,1,4]

;

impulse(num,den)

【结果如右图】

3.已知二阶系统的状态方程为：

求系统的零输入响应和脉冲响应。

【程序如下】：

【结果如下图】

a=[0,1

;

-10,-2]

;

b=[0

;

c=[1,0]

;

d=[0]

;

x0=[1,0];

subplot(1,2,1)

;

initial(a,b,c,d,x0)

subplot(1,2,2)

;

impulse(a,b,c,d)

4.系统传递函数为：

输入正弦信号时，观察输出信号的相位差。

【程序如下】：

【结果如下图】：

num=[1]

;

den=[1,1]

;

t=0

0.01

;

u=sin(2*t)

;

hold

plot(t,u,‘r’)

lsim(num,den,u,t)

5.有一二阶系统，求出周期为4秒的方波的输出响应

【程序如下】：

【结果如下图】：

num=[2

1];

den=[1

3];

t=(0:.1:10);

period=4;

u=(rem(t,period)>=period./2);

lsim(num,den,u,t);

6.已知开环系统传递函数，绘制系统的根轨迹，并分析其稳定性

【程序如下】：

num=[1

2];den1=[1

3];den=conv(den1,den1);figure(1)

rlocus(num,den)

[k,p]=

rlocfind(num,den)

figure(2)

k=55;num1=k*[1

2];den=[1

3];den1=conv(den,den);[num,den]=cloop(num1,den1,-1);

impulse(num,den)

title(

impulse

response

(k=55)

)

figure(3)

k=56;num1=k*[1

2];den=[1

3];den1=conv(den,den);[num,den]=cloop(num1,den1,-1);

impulse(num,den)

title(

impulse

response(k=56)

)

【结果】：

K=55时系统稳定，k=56时系统发散。

作如下系统的bode图

【程序如下】：

【结果】：

n=[1,1]

;

d=[1,4,11,7]

;

bode(n,d)

8.系统传函如下

求有理传函的频率响应，然后在同一张图上绘出以四阶伯德近似表示的系统频率响应

【程序如下】：

num=[1];den=conv([1

2],conv([1

2],[1

2]));

w=logspace(-1,2);

t=0.5;

[m1,p1]=bode(num,den,2);

p1=p1-t*w180/pi;

[n2,d2]=pade(t,4);

numt=conv(n2,num);dent=(conv(den,d2));

[m2,p2]=bode(numt,dent,w);

subplot(2,1,1);

semilogx(w,20*log10(m1),w,20*log10(m2),k--

);

title(

bode

plot

);xlabel(

frequency

);ylabel(

gain

);

subplot(2,1,2);

semilogx(w,p1,w,p2,k--

);

xlabel(

frequency

);ylabel(

phase

);

【结果如下图】：

9.已知系统模型为

求它的幅值裕度和相角裕度

【程序如下】：

n=[3.5];

d=[1

2];

[Gm,Pm,Wcg,Wcp]=margin(n,d)

【结果】：

=1.1433

=7.1688

Wcg

=1.7323

Wcp

=1.6541

10.二阶系统为：

令wn=1,分别作出ξ=2,1,0.707,0.5时的nyquist曲线。

【程序如下】：

n=[1]

;

d1=[1,4,1]

;

d2=[1,2,1]

;

d3=[1,1.414,1];

d4=[1,1,1];

nyquist(n,d1)

;

hold

nyquist(n,d2)

;

nyquist(n,d3)

;

nyquist(n,d4)

;

【结果如下图】：

11.一多环系统，其结构图如下，使用Nyquist频率曲线判断系统的稳定性。

【程序如下】：

k1=16.7/0.0125;z1=[0];

p1=[-1.25

-4

-16];

[num1,den1]=zp2tf(z1,p1,k1);

[num,den]=cloop(num1,den1);

[z,p,k]=tf2zp(num,den);

figure(1)

nyquist(num,den)

figure(2)

[num2,den2]=cloop(num,den);

impulse(num2,den2);

【结果如下图】：

12.

已知系统为：

作该系统的nichols曲线。

【程序如下】：

n=[1]

;

d=[1,1,0]

;

ngrid(‘new’)

;

nichols(n,d)

;

【结果如下图】：

实验六

连续系统数字仿真的基本算法

实验任务

1．理解欧拉法和龙格-库塔法的基本思想；

2．理解数值积分算法的计算精度、速度、稳定性与步长的关系；

取h=0.2，试分别用欧拉法、RK2法和RK4法求解微分方程的数值解，并比较计算精度。

注：解析解：

【程序】

clear

t(1)=0;

置自变量初值

y(1)=1;

y_euler(1)=1;

y_rk2(1)=1;

y_rk4(1)=1;

置解析解和数值解的初值

h=0.2;

步长

求解析解

for

k=1:5

t(k+1)=t(k)+h;

y(k+1)=sqrt(1+2*t(k+1));

end

利用欧拉法求解

for

k=1:5

y_euler(k+1)=y_euler(k)+h*(y_euler(k)-2*t(k)/y_euler(k));

end

利用RK2法求解

for

k=1:5

k1=y_rk2(k)-2*t(k)/y_rk2(k);

k2=(y_rk2(k)+h*k1)-2*(t(k)+h)/(y_rk2(k)+h*k1);

y_rk2(k+1)=y_rk2(k)+h*(k1+k2)/2;

end

利用RK4法求解

for

k=1:5

k1=y_rk4(k)-2*t(k)/y_rk4(k);

k2=(y_rk4(k)+h*k1/2)-2*(t(k)+h/2)/(y_rk4(k)+h*k1/2);

k3=(y_rk4(k)+h*k2/2)-2*(t(k)+h/2)/(y_rk4(k)+h*k2/2);

k4=(y_rk4(k)+h*k3)-2*(t(k)+h)/(y_rk4(k)+h*k3);

y_rk4(k+1)=y_rk4(k)+h*(k1+2*k2+2*k3+k4)/6;

end

disp(

时间

解析解

欧拉法

RK2法

RK4法

)

yt=[t,y,y_euler,y_rk2,y_rk4

];

disp(yt)

【运行结果】

时间

解析解

欧拉法

RK2法

RK4法

1.0000

0.2000

1.1832

1.2000

1.1867

1.1832

0.4000

1.3416

1.3733

1.3483

1.3417

0.6000

1.4832

1.5315

1.4937

1.4833

0.8000

1.6125

1.6811

1.6279

1.6125

1.0000

1.7321

1.8269

1.7542

1.7321

%H=0.05

时间

解析解

欧拉法

RK2法

RK4法

1.0000

0.0500

1.0488

1.0500

1.0489

1.0488

0.1000

1.0954

1.0977

1.0956

1.0954

0.1500

1.1402

1.1435

1.1403

1.1402

0.2000

1.1832

1.1876

1.1834

1.1832

0.2500

1.2247

1.2301

1.2250

1.2247

H=0.001

时间

解析解

欧拉法

RK2法

RK4法

1.0000

0.0010

1.0010

0.0020

1.0020

0.0030

1.0030

0.0040

1.0040

0.0050

1.0050

%欧拉法需要步长足够小时才逼近解析解，RK4逼近得最快，其次是RK2

考虑如下二阶系统：

在上的数字仿真解（已知：，），

并将不同步长下的仿真结果与解析解进行精度比较。

说明：

已知该微分方程的解析解分别为：

采用RK4法进行计算，选择状态变量：

则有如下状态空间模型及初值条件

采用RK4法进行计算。程序如下：

clear

h=input(

请输入步长h=

);

输入步长

M=round(10/h);

置总计算步数

t(1)=0;

置自变量初值

y_0(1)=100;

y_05(1)=100;

置解析解的初始值(y_0和y_05分别对应于为R=0和R=0.5)

x1(1)=100;

x2(1)=0;

置状态向量初值

y_rk4_0(1)=x1(1);

y_rk4_05(1)=x1(1);

置数值解的初值

求解析解

for

k=1:M

t(k+1)=t(k)+h;

y_0(k+1)=100*cos(t(k+1));

y_05(k+1)=100*exp(-t(k+1)/2).*cos(sqrt(3)/2*t(k+1))+100*sqrt(3)/3*exp(-t(k+1)/2).*sin(sqrt(3)/2*t(k+1));

end

利用RK4法求解

R=0

for

k=1:M

k11=x2(k);

k12=-x1(k);

k21=x2(k)+h*k12/2;

k22=-(x1(k)+h*k11/2);

k31=x2(k)+h*k22/2;

k32=-(x1(k)+h*k21/2);

k41=x2(k)+h*k32;

k42=-(x1(k)+h*k31);

x1(k+1)=x1(k)+h*(k11+2*k21+2*k31+k41)/6;

x2(k+1)=x2(k)+h*(k12+2*k22+2*k32+k42)/6;

y_rk4_0(k+1)=x1(k+1);

end

R=0.5

for

k=1:M

k11=x2(k);

k12=-x1(k)-x2(k);

k21=x2(k)+h*k12/2;

k22=-(x1(k)+h*k11/2)-(x2(k)+h*k12/2);

k31=x2(k)+h*k22/2;

k32=-(x1(k)+h*k21/2)-(x2(k)+h*k22/2);

k41=x2(k)+h*k32;

k42=-(x1(k)+h*k31)-(x2(k)+h*k32);

x1(k+1)=x1(k)+h*(k11+2*k21+2*k31+k41)/6;

x2(k+1)=x2(k)+h*(k12+2*k22+2*k32+k42)/6;

y_rk4_05(k+1)=x1(k+1);

end

求出误差最大值

err_0=max(abs(y_0-y_rk4_0));

err_05=max(abs(y_05-y_rk4_05));

输出结果

disp(

最大误差（R=0)

最大误差（R=0.5)

)

err_max=[err_0,err_05];

disp(err_max)

步长h取0.0001

到0.5之间变化，并且将数值解直接与解析解比较，求出最大误差数据列入下表中。

步长h

0.0001

0.0005

0.001

0.005

0.01

0.05

0.1

0.5

R=0

最大误差

5.4330×10-10

1.6969×10-10

1.0574×10-10

4.1107×10-9

6.6029×10-8

4.1439×10-5

6.6602×10-4

4.2988×10-1

R=0.5

最大误差

2.7649×10-11

6.8123×10-12

5.3753×10-12

4.0902×10-10

6.5425×10-9

4.1365×10-6

6.7152×10-5

4.5976×10-2

从上表中可以看出，当步长h=0.001时，总误差最小；当步长h小于0.001时，由于舍入误差变大而使总误差增加；当步长h大于0.001时，则由于截断误差的增加也使得总误差加大。另外，当系统的解变化激烈时（如R=0），误差对步长的变化较为敏感；当系统的解变化平稳时，步长的变化对误差的影响就要缓和得多。数值积分算法确定以后，在选择步长时，需要综合考虑。