计算机仿真实验报告

计算机仿真实验报告

班级：

姓名：

学号：

指导教师：

昆明理工大学信息工程与自动化学院

2012年12月

实验一常微分方程的求解及系统数学模型的转换一．实验目的

通过实验熟悉计算机仿真中常用到的Matlab指令的使用方法，掌握常微分方程求解指令和模型表示及转换指令，为进一步从事有关仿真设计和研究工作打下基础。

二. 实验设备

个人计算机，Matlab软件。

三. 实验准备

预习本实验有关内容（如教材第2、3、5章中的相应指令说明和例题），编写本次仿真练习题的相应程序。

四. 实验内容

1. Matlab中常微分方程求解指令的使用

题目一：请用MATLAB的ODE45算法分别求解下列二个方程。要求:1.编写出Matlab 仿真程序；2.画出方程解的图形并对图形进行简要分析；3.分析下列二个方程的关系。

1．2．

方程1：

M文件：

function xprim=xprim(t,x)

xprim=-x^2

程序：

[t,x]=ode45('xprim',[0,1],1);

plot(t,x,'o',t,x,'-');

xlabel('time t0=0,tt=1');

ylabel('X values X(0)=1');

grid

图像：

方程2：

M文件：

function xprim=xprim2(t,x) xprim=x^2

程序：

[t,x]=ode45('xprim2',[0,1],-1); plot(t,x,'o',t,x,'-');

xlabel('time t0=0,tt=1'); ylabel('X values X(0)=-1'); grid

图像：

题目二：下面方程组用在人口动力学中，可以表达为单一化的捕食者-被捕食者模式（例如，狐狸和兔子）。其中1x 表示被捕食者， 2x 表示捕食者。如果被捕食者有无限的食物，并且不会出现捕食者。于是有1'1x x ，则这个式子是以指数形式增长的。大量的被捕食者将会使捕食者的数量增长；同样，越来越少的捕食者会使被捕食者的数量增长。而且，人口数量也会增长。请分别调用ODE45、ODE23算法求解下面方程组。要求编写出Matlab 仿真程序、画出方程组解的图形并对图形进行分析和比较。

题目二程序及其运行结果： ODE45算法： M 文件：

function fun=fun(t,x)

fun=[x(1)-0.1*x(1)*x(2)+0.01*t;-x(2)+0.02*x(1)*x(2)+0.04*t]; 程序：

[t,x]=ode45('fun',[0,20],[30;20]);

plot(t,x);

xlabel('time t0=0,tt=20');

ylabel('x values x1(0)=30,x2(0)=20');

ODE23解法：

M文件：

function fun=fun(t,x)

fun=[x(1)-0.1*x(1)*x(2)+0.01*t;-x(2)+0.02*x(1)*x(2)+0.04*t]; 程序：

[t,x]=ode23('fun',[0,20],[30;20]);

plot(t,x);

xlabel('time t0=0,tt=20');

ylabel('x values x1(0)=30,x2(0)=20');

2. Matlab 中模型表示及模型转换指令的使用

题目三：若给定系统的的传递函数为

1

13210

6126)(2

3423+++++++=s s s s s s s s G 请用MATLAB 编程求解其系统的极零点模型。 程序：

num=[6,12,6,10];den=[1,2,3,1,1];

[z,p,k]=tf2zp(num,den);sys=zpk(z,p,k) 运行结果：

题目四：习题2.4 num=[1,4,5]; den=[1,6,11,6]; sys=tf(num,den); canon(sys,'modal')

题目五：习题5.8

程序：

D=zpk(0,[0.1,0.2],1,0.02)

D2=d2d(D,0.1)

程序运行结果：

>> D=zpk(0,[0.1,0.2],1,0.02)

Zero/pole/gain:

z

---------------

(z-0.1) (z-0.2)

Sampling time: 0.02

>> D2=d2d(D,0.1)

Zero/pole/gain:

1.3881 (z+0.0002381)

----------------------

(z-1e-005) (z-0.00032) Sampling time: 0.1 五．总结与体会

通过本次实验了解了计算机仿真中常用到的Matlab 指令的使用方法，掌握了常微分方程求解指令和模型表示及转换指令，为进一步学习仿真设计打下了基础。

实验二 Matlab 优化工具箱的使用

一．实验目的

通过上机操作熟悉Matlab 优化工具箱的主要功能及其使用方法，掌握优化工具箱中常用函数的功能和语法，并利用其进行极值运算、求解线性和非线性问题等，为进一步的仿真设计和研究打下基础。

二. 实验设备

个人计算机，Matlab 软件。

三. 实验准备

预习本实验有关内容（如教材第6章中的相应指令说明和例题），编写本次仿真练习题的相应程序。

四. 实验内容

1. 应用Matlab 优化工具箱求解优化问题

例题6.6~6.10，选做2题，要求自行修改方程系数，并比较运行结果。 例题6.7（原）：min 432142x x x x +++ s.t. 8434321≤-++x x x x 6324321≤+-+x x x x 94431≤++x x x 0,,,4321≥x x x x 解：f=[2 4 1 1];

A=[1 3 4 -1;2 1 -3 1;1 0 4 1]; b=[8 6 9];

LBnd=[0 0 0 0];

[x,fval]=linprog(f,A,b,[],[],LBnd,[]) Optimization terminated.

x =

1.0e-012* 0.0078 0.0259 0.2098 0.0860

fval =

4.1508e-013

6.7（修改）：min 432122x x x x +++ s.t. 62424321≤-++x x x x 8324321≤+-+x x x x 622431≤++x x x 0,,,4321≥x x x x 解：f=[1 2 2 1];

A=[1 2 4 -2;2 1 -3 1;1 0 2 2]; b=[6 8 6];

LBnd=[0 0 0 0];

[x,fval]=linprog(f,A,b,[],[],LBnd,[]) Optimization terminated. x =

1.0e-009 *

0.2255 0.0452 0.0324 0.0467

fval =

4.2736e-010

6.9 （原）求解下面的约束非线性规划问题

s.t. 122

12122121212min ()(42421)

1.50100

x f x e x x x x x x x x x x x =+++++--≤--≤

function [f,g]=fun(x)

f=exp(x(1))*(4*x(1)^2+2*x(2)^2+4*x(1)*x(2)+2*x(2)+1); g(1)=1.5+x(1)*x(2)-x(1)-x(2); g(2)=-x(1)*x(2)-10; 调用：

x0=[-1 1];options=[];[x,options]=constr('fun',x0,options) x =

-9.5474 1.0474

options =

Columns 1 through 15

0 0.0001 0.0001 0.0000 0 0 0 0.0236 0 29.0000 10.0000 10.0000 0 200.0000 0

Columns 16 through 18

0.0000 0.1000 1.0000

输出最优解为options(8)=0.0236

6.9（修改）求解下面的约束非线性规划问题

s.t.

1

22

12122

1212

12

min()(221)

1.50

100

x

f x e x x x x x

x x x x

x x

=++++

+--≤

--≤

function [f,g]=fun(x)

f=exp(x(1))*(2*x(1)^2+x(2)^2+2*x(1)*x(2)+x(2)+1);

g(1)=1.5+x(1)*x(2)-x(1)-x(2);

g(2)=-x(1)*x(2)-10;

调用：

x0=[-1 1];options=[];[x,options]=constr('fun',x0,options)

x =

-9.5474 1.0474

options =

Columns 1 through 15

0 0.0001 0.0001 0.0000 0 0

0 0.0118 0 25.0000 9.0000 9.0000 0 200.0000 0

Columns 16 through 18

0.0000 0.1000 1.0000

输出最优解为options(8)=0.0118

对修改前后分析可以看出，前七列的值相同，后六列的值也相同，同时最优解均出现在第八列上。

2. 应用Matlab优化工具箱求解极值问题

已知函数f(x)=10*exp(-x)*cos(x) ，求函数的极值。

（1）x∈[2，5]时，求函数的最小值，并画出函数的曲线。

（2）x∈[3，9]时，求函数的最大值，并画出函数的曲线。

(1)

f='10*exp(-x)*cos(x)';

fplot(f,[2,5])

X=fmin('10*exp(-x)*cos(x)',2,5)

X =

2.3562

Y=10*exp(-2.3562)*cos(2.3562)

Y =

-0.6702

(2)

f='10*exp(-x)*cos(x)';

fplot(f,[3,9]);

X=fmin('-10*exp(-x)*cos(x)',3,9)

X =

5.4978

Y=10*exp(-5.4978)*cos(5.4978)

Y =

0.0290

五．总结与体会

通过这次实验让我对matlab优化工具箱有了一定的了解，对一些函数的使用基本比较熟悉了，像线性规划、最大最小值函数、约束条件下的优化函数都能使用它们来解决一些问题，特别是习题中的问题，对它们的功能和语法都掌握了。

实验三 利用Matlab 和Simulink 进行系统仿真设计

一．实验目的

通过实验对一个汽车运动控制系统进行实际设计与仿真，掌握控制系统性能的分析和仿真处理过程，熟悉用Matlab 和Simulink 进行系统仿真的基本方法。

二. 实验设备

个人计算机，Matlab 软件。

三. 实验准备

预习本实验相关说明，复习PID 控制器的原理和作用，明确汽车运动控制系统问题的描述及其模型表示，编写本次仿真练习的相应程序。

四. 实验说明

本实验是对一个汽车运动控制系统进行实际设计与仿真，其方法是先对汽车运动控制系统进行建摸，然后对其进行PID 控制器的设计，建立了汽车运动控制系统的模型后，可采用Matlab 和Simulink 对控制系统进行仿真设计。

注意：设计系统的控制器之前要观察该系统的开环阶跃响应，采用阶跃响应函数step( )来实现，如果系统不能满足所要求达到的设计性能指标，需要加上合适的控制器。然后再按照仿真结果进行PID 控制器参数的调整，使控制器能够满足系统设计所要求达到的性能指标。

五. 实验内容 1. 问题的描述

如下图所示的汽车运动控制系统，设该系统中汽车车轮的转动惯量可以忽略不计，并且假定汽车受到的摩擦阻力大小与汽车的运动速度成正比，摩擦阻力的方向与汽车运动的方向相反，这样，该汽车运动控制系统可简化为一个简单的质量阻尼系统。

根据牛顿运动定律，质量阻尼系统的动态数学模型可表示为：

?

?

?==+v y u bv v

m 系统的参数设定为：汽车质量m =1000kg ，

比例系数b =50 N ·s/m ， 汽车的驱动力u =500 N 。

根据控制系统的设计要求，当汽车的驱动力为500N 时，汽车将在5秒内达到10m/s 的最大速度。由于该系统为简单的运动控制系统，因此将系统设计成10％的最大超调量和2％的稳态误差。这样，该汽车运动控制系统的性能指标可以设定为：

上升时间：t r <5s ；

最大超调量：σ％<10％； 稳态误差：e ssp <2％。

2.系统的模型表示

假定系统的初始条件为零，则该系统的Laplace 变换式为：

?

?

?==+)()()

()()(s V s Y s U s bV s msV 即 )()()(s U s bY s msY =+ 则该系统的传递函数为：

b

ms s U s Y +=1

)()( 如果用Matlab 语言表示该系统的传递函数模型，相应的程序代码如下： num=1;den=[1000,50];sys=tf(num,den)

同时，系统的数学模型也可写成如下的状态方程形式：

?????=+-=v

y u m v m b v 1 如果用Matlab 语言表示该系统状态空间模型，相应的程序代码如下：

A=-50/1000;B=1/1000;C=1;D=0;sys=ss(A,B,C,D)

3. 系统的仿真设计

● 利用Matlab 进行仿真设计

I ．求系统的开环阶跃响应

在Matlab 命令窗口输入相应的程序代码，得出该系统的模型后，接着输入下面的指令：

step(u*sys)

可得到该系统的开环阶跃响应曲线，如下图所示：

从图上可看出该系统不能满足系统设计所要求达到的性能指标，需要加上合适的控制器。

II ．PID 控制器的设计

PID 控制器的传递函数为：

s

K s K s K s K s K K I

P D D I P ++=++2在PID 控制中，比例（P ）、积分（I ）、

微分（D ）这三种控制所起的作用是不同的（请注意在实验总结中进行归纳）。

下面分别讨论其设计过程。

（1）比例（P ）控制器的设计

增加比例控制器之后闭环系统的传递函数为：

)

()()(P P

K b ms K s U s Y ++=

由于比例控制器可以改变系统的上升时间，现在假定Kp ＝100，观察一下系统的阶跃响应。在MATLAB 命令窗口输入指令：

num=100;den=[1000,150];sys=tf(num,den);step(500*sys)

由此仿真结果，分析系统的稳态值是否满足设计要求，系统的稳态误差和上升时间能不能满足设计要求？

上升时间为40s远远大于5s，不能满足设计要求，稳态误差为（500-333）/500远远大于2%，因此系统不满足设计要求

若减小汽车的驱动力为10N，重新进行仿真，仿真结果为：

num=100;den=[1000,150];sys=tf(num,den);step(10*sys)

如果所设计的比例控制器仍不能满足系统的稳态误差和上升时间的设计要求，则可以通过提高控制器的比例增益系数来改善系统的输出。例如把比例增益系数

Kp 从100提高到10000重新计算该系统的阶跃响应，结果为：

程序：num=10000;den=[1000,10050];sys=tf(num,den);step(10*sys)

此时系统的稳态误差接近为零，系统上升时间也降到了0.5s 以下。这样做虽然满足了系统性能要求，但实际上该控制过程在现实中难以实现。因此，引入比例积分（PI ）控制器来对系统进行调节。 （2）比例积分（PI ）控制器的设计

采用比例积分控制的系统闭环传递函数可表示为：

I

P I P K s K b ms K s K s U s Y ++++=

)()()

(2 增加积分环节的目的是减小系统的稳态误差，假设比例系数Kp =600，积分系数K I =1，编写相应的MATLAB 程序代码如下：

num=[600,1];den=[1000,650,1];sys=(num,den)

可以调节控制器的比例和积分系数来满足系统的性能要求。例如选择比例系数K P ＝800，积分系数K I =40时，可得系统阶跃响应曲线为：

num=[800,40];den=[1000,850,40];sys=(num,den);step(500*sys)

可见，此时的控制系统已经能够满足系统要求达到的性能指标设计要求。但此控制器无微分项，而对于有些实际控制系统往往需要设计完整的PID 控制器，以便同时满足系统的动态和稳态性能要求。 （3）比例积分微分（PID ）控制器的设计

采用PID 控制的系统闭环传递函数为：

I

P D I P d K s K b s K m K s K s K s U s Y ++++++=

)()()()()

(22 假设该控制器的比例系数K P =1，积分系数K I =1，微分系数K D =1，编写MATLAB 程序代码如下：num=[1,1,1];den=[1001,51,1];sys=(num,den)

运行上述程序，并且调整PID控制器的控制参数，直到控制器满足系统设计的性能指标要求为止。

num=[10,650,50];den=[1010,700,50];sys=tf(num,den);step(500*sys)

最后，选择K P=650，K I =50，K D=10，此时系统的阶跃响应曲线如下：

从图中可以看出该系统能够满足设计的总体性能要求。

●利用Simulink进行仿真设计

I．求系统的开环阶跃响应

利用Simulink建立系统阶跃响应模型，如下图所示。双击Step模块，设置模块属性：跳变时间为0；初始值为0；终止值为10；采样时间为0。

单击◢按钮开始仿真，双击Scope模块，可得系统阶跃响应曲线。