《 计 算 机 仿 真》
上 机 实 验 报 告
姓名：
学号：
专业：
班级：
实验一 常微分方程的求解及系统数学模型的转换
一．实验目的
通过实验熟悉计算机仿真中常用到的Matlab指令的使用方法，掌握常微分方程求解指令和模型表示及转换指令，为进一步从事有关仿真设计和研究工作打下基础。
二. 实验设备
个人计算机，Matlab软件。
三. 实验准备
预习本实验有关内容（如教材第2、3、5章中的相应指令说明和例题），编写本次仿真练习题的相应程序。
四. 实验内容
1. Matlab中常微分方程求解指令的使用
题目一：请用MATLAB的ODE45算法分别求解下列二个方程。要求:1.编写出Matlab仿真程序；2.画出方程解的图形并对图形进行简要分析；3.分析下列二个方程的关系。
1． 2．
1.仿真程序
xp1=@(t,x)-x^2;
[t,x]=ode45(xp1,[0,20],[1]);
plot(t,x);title('xp1')
grid
xp2=@(t,x)x^2;
[t,x]=ode45(xp2,[0,20],[-1]);
figure (2);plot(t,x);xlabel('t');
ylabel('x')
grid
2.方程解的图形图形进行简要分析
3.二个方程的关系
题目二：下面方程组用在人口动力学中，可以表达为单一化的捕食者-被捕食者模式（例如，狐狸和兔子）。其中
表示被捕食者， 
表示捕食者。如果被捕食者有无限的食物，并且不会出现捕食者。于是有
，则这个式子是以指数形式增长的。大量的被捕食者将会使捕食者的数量增长；同样，越来越少的捕食者会使被捕食者的数量增长。而且，人口数量也会增长。请分别调用ODE45、ODE23算法求解下面方程组。要求编写出Matlab仿真程序、画出方程组解的图形并对图形进行分析和比较。
function xp3=xp3(t,x)
xp3=[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('xp3',[0,20],[30;20]);
plot(t,x);grid
[t,x]=ode23('xp3',[0,20],[30;20]);
plot(t,x);grid
ODE45
ODE23
2. Matlab中模型表示及模型转换指令的使用
题目三：若给定系统的的传递函数为
请用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=conv([1 0],[1 5 6]);
sys=tf(num,den);
canon(sys,'modal')
题目五：习题5.8
采样周期 Ts=0.02s 先在t=0.1s仿真
dnum=[1,0 ];
dden=[1 -0.3 0.02];
sys1=tf(dnum,dden,0.02)
sys2=d2d(sys1,0.1)
五．总结与体会
实验二 Matlab优化工具箱的使用
一．实验目的
通过上机操作熟悉Matlab优化工具箱的主要功能及其使用方法，掌握优化工具箱中常用函数的功能和语法，并利用其进行极值运算、求解线性和非线性问题等，为进一步的仿真设计和研究打下基础。
二. 实验设备
个人计算机，Matlab软件。
三. 实验准备
预习本实验有关内容（如教材第6章中的相应指令说明和例题），编写本次仿真练习题的相应程序。
四. 实验内容
1. 应用Matlab优化工具箱求解优化问题
例题6.6~6.10，选做2题，要求自行修改方程系数，并比较运行结果。
例6.6 例如我们希望求解非线性方程：
初始值为x0 = [-5 -5]
function F = myfun(x)
F = [2*x(1) - x(2) - exp(-x(1)); -x(1) + 2*x(2) - exp(-x(2))];
x0 = [-5; -5];
options=optimset('Display','iter');
[x,fval] = fsolve(@myfun,x0,options)
例6.7 利用Matlab语言求解下列线性规划问题
f=[-2 -1.3 -5]';
A=[1 2 4 -1; 2 3 -1 1; 1 0 1 1; -1 0 0 0; 0 -1 0 0; 0 0 -1 0; 0 0 0 -1];
b=[6 12 4 0 0 0 0]';
X=linprog(f,A,b)
例6.8：利用Matlab命令求解下面的无约束非线性规划问题。
function f=fun(x,y) f=exp(x(1))*(4*x(1)^2+2*x(2)^2+4*x(1)*x(2)+2*x(2)+1)
function f=fun(x,y) f=exp(x)*(4*x^2+2*y^2+4*x*y+2*y+1);
options=[];
[x,options]=fminunc(fun,[-1 1],options)
2. 应用Matlab优化工具箱求解极值问题
已知函数f(x)=10*exp(-x)*cos(x) ，求函数的极值。
（1）x∈[2，5]时，求函数的最小值，并画出函数的曲线。
（2）x∈[3，9]时，求函数的最大值，并画出函数的曲线。
(1)[x,f]=fminbnd('10*exp(-x)*cos(x)',2,5),fplot('10*exp(-x)*cos(x)',[2 5])
(2)f='-10*exp(-x)*cos(x)';[x,f]=fminbnd(f,3,9),fplot('10*exp(-x)*cos(x)',[3 9])
a=inline('-10*exp(-x)*cos(x)','x');[x,a]=fminbnd(a,3,9)
f=-10*exp(-x)*cos(x);a=inline('f');[x,a]=fminbnd(a,3,9)
五．总结 与体会
实验三利用Matlab和Simulink进行系统仿真设计
一．实验目的
通过实验对一个汽车运动控制系统进行实际设计与仿真，掌握控制系统性能的分析和仿真处理过程，熟悉用Matlab和Simulink进行系统仿真的基本方法。
二. 实验设备
个人计算机，Matlab软件。
三. 实验准备
预习本实验相关说明，复习PID控制器的原理和作用，明确汽车运动控制系统问题的描述及其模型表示，编写本次仿真练习的相应程序。
四. 实验说明
本实验是对一个汽车运动控制系统进行实际设计与仿真，其方法是先对汽车运动控制系统进行建摸，然后对其进行PID控制器的设计，建立了汽车运动控制系统的模型后，可采用Matlab和Simulink对控制系统进行仿真设计。
注意：设计系统的控制器之前要观察该系统的开环阶跃响应，采用阶跃响应函数step( )来实现，如果系统不能满足所要求达到的设计性能指标，需要加上合适的控制器。然后再按照仿真结果进行PID控制器参数的调整，使控制器能够满足系统设计所要求达到的性能指标。
五. 实验内容
1. 问题的描述
如下图所示的汽车运动控制系统，设该系统中汽车车轮的转动惯量可以忽略不计，并且假定汽车受到的摩擦阻力大小与汽车的运动速度成正比，摩擦阻力的方向与汽车运动的方向相反，这样，该汽车运动控制系统可简化为一个简单的质量阻尼系统。
根据牛顿运动定律，质量阻尼系统的动态数学模型可表示为：
系统的参数设定为：汽车质量m=1000kg，
比例系数b=50 N·s/m，
汽车的驱动力u=500 N。
根据控制系统的设计要求，当汽车的驱动力为500N时，汽车将在5秒内达到10m/s的最大速度。由于该系统为简单的运动控制系统，因此将系统设计成10％的最大超调量和2％的稳态误差。这样，该汽车运动控制系统的性能指标可以设定为：
上升时间：tr<5s；
最大超调量：σ％<10％；
稳态误差：essp<2％。
2.系统的模型表示
假定系统的初始条件为零，则该系统的Laplace变换式为：
即
则该系统的传递函数为：
如果用Matlab语言表示该系统的传递函数模型，相应的程序代码如下：
g=tf(1,[1000,50]),
同时，系统的数学模型也可写成如下的状态方程形式：
如果用Matlab语言表示该系统状态空间模型，相应的程序代码如下：
sys=ss(-0.05,0.001,1,0),
3. 系统的仿真设计
● 利用Matlab进行仿真设计
I．求系统的开环阶跃响应
在Matlab命令窗口输入相应的程序代码，得出该系统的模型后，接着输入下面的指令：
step(500*g),step(500*sys),
可得到该系统的开环阶跃响应曲线，如下图所
从图上可看出该系统不能满足系统设计所要求达到的性能指标，需要加上合适的控制器。
II．PID控制器的设计
PID控制器的传递函数为：
在PID控制中，比例（P）、积分（I）、微分（D）这三种控制所起的作用是不同的（请注意在实验总结中进行归纳）。下面分别讨论其设计过程。
（1）比例（P）控制器的设计
增加比例控制器之后闭环系统的传递函数为：
由于比例控制器可以改变系统的上升时间，现在假定Kp＝100，观察一下系统的阶跃响应。在MATLAB命令窗口输入指令：
g=tf(100,[1000,150]),step(500*g)
可得系统阶跃响应如下：
由此仿真结果，分析系统的稳态值是否满足设计要求，系统的稳态误差和上升时间能不能满足设计要求？
若减小汽车的驱动力为10N，重新进行仿真，仿真结果为：
g=tf(100,[1000,150]),step(10*g)
如果所设计的比例控制器仍不能满足系统的稳态误差和上升时间的设计要求，则可以通过提高控制器的比例增益系数来改善系统的输出。例如把比例增益系数Kp从100提高到10000重新计算该系统的阶跃响应，结果为：
g=tf(10000,[1000,10050]),step(500*g)
此时系统的稳态误差接近为零，系统上升时间也降到了0.5s以下。这样做虽然满足了系统性能要求，但实际上该控制过程在现实中难以实现。因此，引入比例积分（PI）控制器来对系统进行调节。
（2）比例积分（PI）控制器的设计
采用比例积分控制的系统闭环传递函数可表示为：
增加积分环节的目的是减小系统的稳态误差，假设比例系数Kp=600，积分系数KI=1，编写相应的MATLAB程序代码如下：
g=tf([600,1],[1000,650,1]),step(500*g)
运行上述程序后可得系统阶跃响应曲线为：
可以调节控制器的比例和积分系数来满足系统的性能要求。例如选择比例系数KP＝800，积分系数KI=40时，可得系统阶跃响应曲线为：
g=tf([800,40],[1000,850,40]),step(500*g)
可见，此时的控制系统已经能够满足系统要求达到的性能指标设计要求。但此控制器无微分项，而对于有些实际控制系统往往需要设计完整的PID控制器，以便同时满足系统的动态和稳态性能要求。
（3）比例积分微分（PID）控制器的设计
采用PID控制的系统闭环传递函数为：
假设该控制器的比例系数KP=1，积分系数KI =1，微分系数KD=1，编写MATLAB程序代码如下：
g=tf([1,1,1],[1001,51,1]),step(500*g),grid
运行上述程序，并且调整PID控制器的控制参数，直到控制器满足系统设计的性能指标要求为止。
最后，选择KP= 800 ，KI = 60 ，KD= 1 ，此时系统的阶跃响应曲线如下：
从图中可以看出该系统能够满足设计的总体性能要求。
● 利用Simulink进行仿真设计
I．求系统的开环阶跃响应
利用Simulink建立系统阶跃响应模型，如下图所示。双击Step模块，设置模块属性：跳变时间为0；初始值为0；终止值为10；采样时间为0。
单击◢按钮开始仿真，双击Scope模块，可得系统阶跃响应曲线。
分析此时系统的性能指标是否满足设计要求？若不能则应加上合适的控制器。
II．PID控制器的设计
在Simulink的模型窗口建立一个包含PID控制器的闭环系统阶跃响应模型，如图所示：
分别双击Kp、Ki、Kd模块设定比例、积分、微分系数，点击◢按钮开始仿真，双击Scope模块，观察系统的阶跃响应曲线，直到满足要求为止。
最终选取KP= ，KI= ，KD= ，此时控制器能满足系统设计所要求达到的性能指标，Simulink仿真的汽车运动PID控制系统的阶跃响应曲线如下：
六．总结与体会
实验总结：
（总结实验中遇到的问题及解决办法，分析PID控制器三种控制所起的作用等。）
第二篇：计算机仿真实验报告-实验三
一、实验内容：
实验三 利用欧拉法、梯形数法和二阶显式Adams法对RLC串联电路的仿真
1前向欧拉法
状态方程：
然后根据前向欧拉法（其中h为步长）
即可得到系统的差分方程
2后向欧拉法
根据前向欧拉法得到的系统状态方程，结合后向欧拉法（其中h为步长）
即可得到系统的差分方程
3梯形法
由前面的系统状态方程，结合梯形法
即可得到系统的差分方程
4二阶显式Adams方法
由前面的状态方程，结合二阶显式Adams方法
即可得到系统的差分方程
但是二阶显式Adams法不能自起步，要使方程起步，需要知道开始的三个值，但是我们只知道第一个值。经过分析后，二阶显式Adams方法精度是二阶的，而梯形法精度也是二阶的，因此我们可以先借助梯形法得到输出的前三个值，以达到起步的目的，然后借助上面得到的差分方程对其进行求解。
二、实验波形：
下图为前向欧拉法、后向欧拉法、梯形法、二阶显式Adams方法的系统差分方程得到相应的输出波形：
图1 h=
时四种方法的输出波形
图2 h=
时四种方法的输出波形
图3 h=
时四种方法的输出波形
图4 h=
时四种方法的输出波形
三、实验分析：
由输出波形可以看到各种方法的特点（在图中蓝色线均表示连续系统模型的实际输出波形，红色线表示在对应方法下系统的输出波形。）：
1前向欧拉法和二阶显式Adams方法对步长的要求很强。步长 太大，最后的到的结果不是绝对收敛，而是发散。在小步长 下才显得收敛，这也从另一方面验证，步长越小，截断误差 越小；
2步长不能太小，太小的步长相应的舍入误差和累积误差也会 增大；
3前向欧拉法也可称为显式欧拉法，后向欧拉法也可称为隐式 欧拉法，可以看到，后向欧拉法的稳定域要比前向欧拉法大， 计算精度也要高一些。这也从另一方面证明同一算法隐式算法 要比显式算法稳定域大，计算精度高；
4以h=
为例来说明各种方法的另外的特点。可以明显看 到，前向欧拉法并没有绝对收敛，到后来发散了，纵向比较来 看，其对步长的依赖性比较大，稳定域比较小。后向欧拉法收 敛速度块，但是跟随性相对来说差了一些，因此计算精度也就 相对差了一些，但是纵向比较来看其对步长的依赖性比较小， 稳定域还是比较大。梯形法跟随性看起来比较好，但是并不是 和实际值完全一致，因此计算精度稍微好一些，当然其很好地 反映了实际的进行过程，纵向比较来看，其对步长的依赖性比 较小，稳定域相对比较大。二阶显式Adams方法收敛速度比 后向欧拉法慢了一些，跟随性比梯形法差了一些，所以计算精 度上比梯形法好了一些，反映实际的过程进行也要比梯形法好 一些，但是在纵向来看，其对步长的依赖比较大，而且稳定域 小一些，还有为了使其启动，还得先借助其他方法得到最初的 几个初值，这客观上增加了实现算法的难度。
四、实验结论：
综合来讲，根据问题的实际需要，选择恰当地处理方法。在此基础上合理地选择步长，以及选择合适的处理机，以及其他的一些方面，借助计算机仿真从而得到合理的结果。