常微分方程的求解及系统数学模型的转换实验报告

《计算机仿真》实验报告

实验一常微分方程的求解及系统数学模型的转换一．实验目的通过实验熟悉计算机仿真中常用到的Matlab指令的使用方法，掌握常微分方程求解指令和模型表示及转换指令，为进一步从事有关仿真设计和研究工作打下基础。二.实验设备个人计算机，Matlab软件。三.实验准备预习本实验有关内容（如教材第2、3、5章中的相应指令说明和例题），编写本次仿真练习题的相应程序。四.实验内容1.Matlab中常微分方程求解指令的使用题目一：请用MATLAB的ODE45算法分别求解下列二个方程。要求:1.编写出Matlab仿真程序；2.画出方程解的图形并对图形进行简要分析；3.分析下列二个方程的关系。2． 解：程序：新建f1：functionf1=f1(t,x)f1=-x^2;end新建f2：functionf2=f2(t,x)f2=x^2;end运行：[t,x]=ode45(@f1,[0,40],[1])figure(1);plot(t,x);[t,x]=ode45(@f2,[0,40],[-1])figure(2);plot(t,x)运行截图：结论：方程一和方程二的图形关于横轴对称题目二：下面方程组用在人口动力学中，可以表达为单一化的捕食者-被捕食者模式（例如，狐狸和兔子）。其中表示被捕食者，表示捕食者。如果被捕食者有无限的食物，并且不会出现捕食者。于是有，则这个式子是以指数形式增长的。大量的被捕食者将会使捕食者的数量增长；同样，越来越少的捕食者会使被捕食者的数量增长。而且，人口数量也会增长。请分别调用ODE45、ODE23算法求解下面方程组。要求编写出Matlab仿真程序、画出方程组解的图形并对图形进行分析和比较。解：OED45：程序：新建f3：functionf3=f3(t,x)f3=[x(1)-0.1*x(1)*x(2)+0.01*t;-x(2)+0.02*x(1)*x(2)+0.04*t];end运行：[t,x]=ode45('f3',[0,20],[30,20]);plot(t,x);title('ode45×÷ͼ');xlabel('ºìÏßΪ²¶Ê³Õߣ¬À¶ÏßΪ±»²¶Ê³Õß');grid运行截图OED23：程序;新建f4：functionf4=f4(t,x)f4=[x(1)-0.1*x(1)*x(2)+0.01*t;-x(2)+0.02*x(1)*x(2)+0.04*t];end运行：[t,x]=ode23('f4',[0,20],[30,20]);plot(t,x);title('ode23×÷ͼ');xlabel('ºìÏßΪ²¶Ê³Õߣ¬À¶ÏßΪ±»²¶Ê³Õß');grid运行截图：2.Matlab中模型表示及模型转换指令的使用题目三：若给定系统的的传递函数为请用MATLAB编程求解其系统的极零点模型。解：程序：num=[612610];den=[12311];[z,p,k]=tf2zp(num,den);Gzpk=zpk(z,p,k)[zeropolesk]=zpkdata(Gzpk,'v')运行结果：Gzpk=6(s+1.929)(s^2+0.07058s+0.8638)-----------------------------------------------(s^2+0.08663s+0.413)(s^2+1.913s+2.421)Continuous-timezero/pole/gainmodel.zero=-1.9294+0.0000i-0.0353+0.9287i-0.0353-0.9287ipoles=-0.9567+1.2272i-0.9567-1.2272i-0.0433+0.6412i-0.0433-0.6412ik=6题目四：例题2.9，求其伴随标准型。解：程序为：num=[41716];den=[171612];sys=tf(num,den)canon(sys,'companion')运行结果为：timusisys=4s^2+17s+16-----------------------s^3+7s^2+16s+12Continuous-timetransferfunction.ans=A=x1x2x3x100-12x210-16x301-7B=u1x11x20x30C=x1x2x3y14-1129D=u1y10Continuous-timestate-spacemodel.题目五：习题5.8，采样周期改变为0.05s，求系统的数学仿真模型。解：程序：dnum=[10];dden=[1-0.30.02];sys=tf(dnum,dden,0.1)D=d2d(sys,0.05)运行结果:timuwusys=z------------------z^2-0.3z+0.02Sampletime:0.1secondsDiscrete-timetransferfunction.D=0.6222z-0.09725-----------------------z^2-0.7634z+0.1414Sampletime:0.05secondsDiscrete-timetransferfunction.五．总结与体会通过本次实验，学会了使用matlab新建和运行函数，和使用ODE23和ODE45函数，以及如何利用matlab画图，学会了matlab中模型表示及模型转换的使用；学会了如何查看和使用图上的信息；学会了求解系统的极零点模型。学会了函数伴随标准型的产生以及系统的数学仿真模型，为以后的学习和实验打下了坚实的基础。

实验二Matlab优化工具箱的使用一．实验目的通过上机操作熟悉Matlab优化工具箱的主要功能及其使用方法，掌握优化工具箱中常用函数的功能和语法，并利用其进行极值运算、求解线性和非线性问题等，为进一步的仿真设计和研究打下基础。二.实验设备个人计算机，Matlab软件。三.实验准备预习本实验有关内容（如教材第6章中的相应指令说明和例题），编写本次仿真练习题的相应程序。四.实验内容1.应用Matlab优化工具箱求解优化问题例题6.6原函数：新建myfun211functionF=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(@myfun211,x0,options)结果为：NormofFirst-orderTrust-regionIterationFunc-countf(x)stepoptimalityradius0347071.22.29e+0411612003.415.75e+031293147.0211.47e+031312854.45213881415239.5271107151867.0412130.8162116.704219.0517242.4278812.2618270.0326580.7595110.2062.59307.03149e-060.1119270.002942.510333.29525e-130.001691326.36e-072.5Equationsolved.fsolvecompletedbecausethevectoroffunctionvaluesisnearzeroasmeasuredbythedefaultvalueofthefunctiontolerance,andtheproblemappearsregularasmeasuredbythegradient.<stoppingcriteriadetails>x=0.56710.5671fval=1.0e-06*-0.4059-0.4059更改参数后：新建myfun2111:functionF=myfun(x)F=[3*x(1)-1*x(2)-exp(-x(1));-2*x(1)+4*x(2)-exp(-x(2))];调用优化函数：x0=[-5;-5];options=optimset('Display','iter');[x,fval]=fsolve(@myfun2111,x0,options)运行结果为：yunxing2111NormofFirst-orderTrust-regionIterationFunc-countf(x)stepoptimalityradius0350189.52.4e+041161337016.17e+031293740.911.67e+0313121108.1414751415343.11911451518104.483147162126.215115.417242.6223414.0318270.009749870.5367630.2152.59301.25882e-070.03695340.0007372.510332.06617e-170.0001337829.08e-092.5Equationsolved.fsolvecompletedbecausethevectoroffunctionvaluesisnearzeroasmeasuredbythedefaultvalueofthefunctiontolerance,andtheproblemappearsregularasmeasuredbythegradient.<stoppingcriteriadetails>x=0.35170.3517fval=1.0e-08*-0.3800-0.2495结论：更改参数前后，运行过程中使用的方法是一样的，运行结果结构也类似。例题6.8:原函数：新建fun212:functionf=fun212(x)f=exp(x(1))*(4*x(1)^2+2*x(2)^2+4*x(1)*x(2)+2*x(2)+1);运行函数：x0=[-11];options=[];[x,options]=fminunc('fun212',x0,options)结果为：>>yunxing212Localminimumfound.Optimizationcompletedbecausethesizeofthegradientislessthanthedefaultvalueoftheoptimalitytolerance.<stoppingcriteriadetails>x=0.5000-1.0000options=3.6609e-16更改后：新建fun2121：functionf=fun212(x)f=exp(x(1))*(5*x(1)^2+6*x(2)^2+7*x(1)*x(2)+8*x(2)+1);运行：x0=[-11];options=[];[x,options]=fminunc('fun2121',x0,options)结果为：>>yunxing2121Localminimumfound.Optimizationcompletedbecausethesizeofthegradientislessthanthedefaultvalueoftheoptimalitytolerance.<stoppingcriteriadetails>x=1.2671-1.4058options=-10.04832.应用Matlab优化工具箱求解极值问题已知函数f(x)=10*exp(-x)*cos(x)，求函数的极值。（1）x∈[2，5]时，求函数的最小值，并画出函数的曲线。解：程序为：x=fminbnd('10*exp(-x)*cos(x)',2,5)f='10*exp(-x)*cos(x)';fplot(f,[2,5])y=10*exp(-x)*cos(x)运行结果为：x=2.3562y=-0.6702（2）x∈[3，9]时，求函数的最大值，并画出函数的曲线。解：程序为：x=fminbnd('-(10*exp(-x)*cos(x))',3,9)f='10*exp(-x)*cos(x)';fplot(f,[3,9])y=10*exp(-x)*cos(x)结果为：x=5.4978y=0.0290五．总结与体会通过本实验熟悉了Matlab优化工具箱的主要功能及其使用方法，掌握了优化工具箱中常用函数的功能和语法，能利用其进行极值运算、求解线性和非线性问题，为接下来的学习和实验进一步打下了基础。

实验三利用Matlab和Simulink进行系统仿真设计一．实验目的通过实验对一个汽车运动控制系统进行实际设计与仿真，掌握控制系统性能的分析和仿真处理过程，熟悉用Matlab和Simulink进行系统仿真的基本方法。二.实验设备个人计算机，Matlab软件。三.实验准备预习本实验相关说明，复习PID控制器的原理和作用，明确汽车运动控制系统问题的描述及其模型表示，编写本次仿真练习的相应程序。四.实验说明本实验是对一个汽车运动控制系统进行实际设计与仿真，其方法是先对汽车运动控制系统进行建摸，然后对其进行PID控制器的设计，建立了汽车运动控制系统的模型后，可采用Matlab和Simulink对控制系统进行仿真设计。注意：设计系统的控制器之前要观察该系统的开环阶跃响应，采用阶跃响应函数step()来实现，如果系统不能满足所要求达到的设计性能指标，需要加上合适的控制器。然后再按照仿真结果进行PID控制器参数的调整，使控制器能够满足系统设计所要求达到的性能指标。五.实验内容1.问题的描述如下图所示的汽车运动控制系统，设该系统中汽车车轮的转动惯量可以忽略不计，并且假定汽车受到的摩擦阻力大小与汽车的运动速度成正比，摩擦阻力的方向与汽车运动的方向相反，这样，该汽车运动控制系统可简化为一个简单的质量阻尼系统。根据牛顿运动定律，质量阻尼系统的动态数学模型可表示为：系统的参数设定为：汽车质量m=1000kg，比例系数b=50N·s/m，汽车的驱动力u=500N。根据控制系统的设计要求，当汽车的驱动力为500N时，汽车将在5秒内达到10m/s的最大速度。由于该系统为简单的运动控制系统，因此将系统设计成10％的最大超调量和2％的稳态误差。这样，该汽车运动控制系统的性能指标可以设定为：上升时间：tr<5s；最大超调量：σ％<10％；稳态误差：essp<2％。2.系统的模型表示假定系统的初始条件为零，则该系统的Laplace变换式为：即则该系统的传递函数为：如果用Matlab语言表示该系统的传递函数模型，相应的程序代码如下：num=1;den=[100050];G=tf(num,den)结果为：G=1-----------1000s+50同时，系统的数学模型也可写成如下的状态方程形式：如果用Matlab语言表示该系统状态空间模型，相应的程序代码如下：A=-1/20;B=1/1000;C=1;sys=ss(A,B,C,0)sys=A=x1x1-0.05B=u1x10.001C=x1y11D=u1y103.系统的仿真设计利用Matlab进行仿真设计=1\*ROMANI．求系统的开环阶跃响应在Matlab命令窗口输入相应的程序代码，得出该系统的模型后，接着输入下面的指令：step(u*sys)可得到该系统的开环阶跃响应曲线，如下图所示：从图上可看出该系统不能满足系统设计所要求达到的性能指标，需要加上合适的控制器。=2\*ROMANII．PID控制器的设计PID控制器的传递函数为：在PID控制中，比例（P）、积分（I）、微分（D）这三种控制所起的作用是不同的（请注意在实验总结中进行归纳）。下面分别讨论其设计过程。（1）比例（P）控制器的设计增加比例控制器之后闭环系统的传递函数为：由于比例控制器可以改变系统的上升时间，现在假定Kp＝100，观察一下系统的阶跃响应。在MATLAB命令窗口输入指令：num=100;den=[1000150];sys=tf(num,den)step(sys)可得系统阶跃响应如下：由此仿真结果，分析系统的稳态值是否满足设计要求，系统的稳态误差和上升时间能不能满足设计要求？由上图知系统的稳态误差为σ%=(0.667-0.6)/0.6*100%=11.17%上升时间为t,=14.6s均不满足设计要求若减小汽车的驱动力为10N，重新进行仿真，仿真结果为：如果所设计的比例控制器仍不能满足系统的稳态误差和上升时间的设计要求，则可以通过提高控制器的比例增益系数来改善系统的输出。例如把比例增益系数Kp从100提高到10000重新计算该系统的阶跃响应，结果为：num=10000;den=[100010050];sys=tf(num,den)step(sys)此时系统的稳态误差接近为零，系统上升时间也降到了0.5s以下。这样做虽然满足了系统性能要求，但实际上该控制过程在现实中难以实现。因此，引入比例积分（PI）控制器来对系统进行调节。（2）比例积分（PI）控制器的设计采用比例积分控制的系统闭环传递函数可表示为：增加积分环节的目的是减小系统的稳态误差，假设比例系数Kp=600，积分系数KI=1，编写相应的MATLAB程序代码如下：num=[6001];den=[10006501];sys=tf(num,den)step(sys)运行上述程序后可得系统阶跃响应曲线为：可以调节控制器的比例和积分系数来满足系统的性能要求。例如选择比例系数KP＝800，积分系数KI=40时，可得系统阶跃响应曲线为：可见，此时的控制系统已经能够满足系统要求达到的性能指标设计要求。但此控制器无微分项，而对于有些实际控制系统往往需要设计完整的PID控制器，以便同时满足系统的动态和稳态性能要求。（3）比例积分微分（PID）控制器的设计采用PID控制的系统闭环传递函数为：假设该控制器的比例系数KP=1，积分系数KI=1，微分系数KD=1，编写MATLAB程序代码如下：num=[111];den=[10015