第四章 Simulink与微分方程组数解

第五章Simulink与微分方程组数解•Simulink简介控制系统的运动规律用微分方程表示，其动态响应可用计算机数值求解，并进而求出系统的超调量，恢复时间等性能参数。Matlab任意微分方程均可用函数ODE45数值求解，同时它还提供了另外一种更直接的图示法对控制系统进行仿真计算，即利用Simulink求线性控制系统的动态响应。•Simulink操作步骤（1）执行File-New-Model菜单MATLAB7.1（2）点击库窗口SIMULINK—CONTNUOUS，选择建立模型方式：传递函数形式（TRANSFERFCN），或状态方程模式（STATE-SPACE），一般可选传递函数，并拉入模型窗口。(3)双击模型，修改参数，按OK或APPLY。可能的话可以重新命名新名字，或改变显示方向。MTALAB6.5分子多项式系数分母多项式系数显示方向设置（FORMAT—FLIP（ROTATE）：MATLAB7.1分子多项式定义分母多项式定义（4）用和、乘、分支、建立复杂模型，和，积用MATH里的函数，分支直接从输入箭头拉直线到输出线，使两条线相交。信号的+，-，比例放大在MathOperation项内信号加：可将+改成-，或增加+，表示多个信号加（5）定义输入、输出，在源（SOURCE）选择输入信号源类型，在SINK选择输出类型，一般设置成SCOPE（示波器观察）。（6）设置计算参数，按SIMULATION—PARAMETER。计算开始时间计算结束时间计算方法MATLAB6.5（7）按SIMULATION—START开始仿真计算。（8）点击SCOPE观察计算结果，或者查看WORKSPACE，观看变量值，根据需要保存绘制计算值。SCOPE:•例212diLURiCdtdJCiMdt已知直流发动机的控制系统的微分方程为：1221216,0.5,0.03,1.3,1.6,107,220ULRCCMLmHRJkgmCCMNmUV为直流电机输入电压，为直流电机的电感，为直流电机的电枢电阻，为电势常数，为力矩常数，为外部负载。其中方法1：Laplace变换简化112()LsiURiCLsRiUCJsCiM根据第1个方程，可得：1C电流iU从Mathoperation选放大器，加法器，把机械运动方程加进去。比例放大，数据可以是具体数值，也可以是变量转速外部负载积分定义输出：从sink选输出(out)示波器看输出输出电流输出转速定义输出计算选项设置：设置微分方程求解方法，计算时间，输出项数。输出变量名，一般不改变保存输出变量个数计算结果矩阵下标序号每隔1输出按Simulation----Start开始计算。查看计算结果，或用plot(tout,yout(:,1))等命令绘制计算结果012345678910050100150绘制转速曲线012345678910-100-50050100150200250300350方法2：根据微分方程直接利用积分，微分，求和绘制系统模型。12.......didLURiCJCiMdtdt积分电流转速第1个方程输入电压U第2个方程：设置输出点：输出转速2dJCiMdt积分电流负载输出点2输出点1012345678910050100150200250电流计算结果转速计算结果：012345678910-50050100150200250方法3：方法2需要对每一项求和，求积，操作复杂，可直接利用微分方程也即状态模型来计算动态响应。状态模型：.....(1)........(2)dXAXBuXUdtYCXDuYXANNBNMUMYrCrNDrM此为状态方程，为状态，为输入此为输出方程，为输出,其中为维数为N的向量，为的矩阵，为的矩阵为维数为的向量，为维数为的向量,为的矩阵，为的矩阵.12.......didLURiCJCiMdtdt如果输出电流和转速，则输出方程为：例：直流电机模型12,010,10CRiLLXACJUdXLBuAXBuMdtJ状态：，1000,,,0100iCYDYCXDu（1）根据上述方程定义矩阵A，B，C，D。（2）设置矩阵，双击选项，出现设置菜单。将默认值1修改成矩阵常数或已经定义好的矩阵变量（3）设置输入，输出设置输入值，可以是向量常数或是已经定义的向量变量（4）设置计算选项，开始计算，并选取计算结果012345678910050100150200250电流计算结果012345678910-50050100150200250转速计算结果：•子系统设计可将模型保存为子系统，求解其他问题时如果要引用该模型，可直接打开复制。输入1输入2子系统模型建立步骤：（1）定义模型（2）定义输入输出（3）全选模型，并选Edit---CreateSubsystem(4)保存子系统到一文件，以后打开该文件显示如下输出输入点子系统名字，可修改（5）引用子系统方法（A）打开保存子系统的文件，出现：（B）修改输入输出，构建新计算模型（C）计算•数值积分计算(1)定义函数2200120sin()cos()|1cos(2)11sin()xdxxdxax菜单:File-new-M_file进入编辑器写程序:functionf=f1(t)f=sin(t);函数编写格式：Function返回变量=函数名(参数)…函数的matlab命令(命令以‘;’结束)返回变量=函数值(最后一条命令必须将函数值赋给返回变量)写好程序后保存到一文件,文件名一般与函数名相同(2)调用函数quad(‘name’,tmin,tmax)求数值积分Name:函数名,tmin,tmax,积分下限,上限例:quad(‘f1’,0,pi/2)Functiony=tuoyuan(t)X=sin(t);Y=1./(1-0.5.*x.*x).^5;2110.5*sin()x例：编写函数NOTE：函数内*，/一般采用.*,./面积分(二重积分)函数dblquad()用法:dblquad(‘name’,xmin,xmax,ymin,ymax)Name:函数名xmin,xmax,ymin,ymax积分区域例:functionout=f2(x,y)out=y*sin(x)+x*cos(y);v=dblquad(‘f2’,0,1,0,5);•微分方程数值求解1.多元一次微分方程组数值解1112221212(,,,...,)(,,,...,)...(,,,...,)nnnnndyftyyydtdyftyyydtdyftyyydt计算步骤如下:(1)定义函数例如:functiondf=f1(t,y)df=[y(1)-y(2);y(1)+y(2)](2)调用函数[t,y]=ode45(‘name’,[tmin,tmax],[y0;y1;…])t为保存自变量的计算点向量变量,y保存计算点处各变量值,为矩阵.y(:,1)表示第1个变量在t各计算点的值,y(1,2)表示第2个变量在t的第1个计算点的值第1各参数为函数名第2个参数为积分范围第3个参数为变量初值(3)绘制曲线Plot(t,y(:,1))第1个参数为自变量第2个参数为向量变量返回值为列向量[T,Y]=solver(odefun,tspan,y0)[T,Y]=solver(odefun,tspan,y0,options)[T,Y,TE,YE,IE]=solver(odefun,tspan,y0,options)sol=solver(odefun,[t0tf],y0...)Solver:ode45,ode113,ode15s,ode23s,ode23t,ode23tb非刚性矩阵：ode45,ode113刚性矩阵：ode15s,ode23s,ode23t,ode23tb刚性矩阵（条件数很大）的矩阵12121212,,(0)0,(0)0dydyyyyydtdtyy例：求微分方程组的数值解已知Functiondy=f2(t,y)dy=[y(1)+y(2);y(1)-y(2)];2.一元高次微分方程数值求解()(1)(1)1(1)(1)(,,','',...)...(,,','',...)(,,','',...)nnnnnnyatyyyyyatyyyyybtyyyy11221...nndyydtdyydtdyydt1111122(,,,...,)...nnnnndybtyyyayayaydt0yy令001111110100000...00,00...1............nnnnyyyydYYdtaaayybFunctiondf=f1(t,y)df=[y(1);-y(1)-(1-y(1)^2)*y(2)];[ts]=ode45(‘f1’,[0,10],[1;1]);例:求微分方程yy2)y+y=0数值解,初值y=1,y=11211[(1)]dyydtdyyyydt