仿真课件 4 经典的连续系统仿真建模方法学2.ppt

Chapter 4,经典的连续系统仿真建模方法学,主要内容,3.1 连续系统离散化原理 数值积分原理，欧拉法，梯形法 3.2 龙格－库塔法 基本原理，步长控制等 3.3 稳定性分析 欧拉法，龙格－库塔法步长限制 3.4 数值积分仿真编程 Matlab程序结构，卫星发射轨道仿真,稳定性分析,,3.3,数值积分算法误差来源： 1）截断误差：基于泰勒公式展开后，舍入高阶项产生的误差截断误差与算法阶次和步长有关，如欧拉法截断误差为O(h2)，RK4法截断误差为O(h5) 2）舍入误差：由于计算机的有限字长引起的每一步计算引起的误差舍入误差会积累，步长越小，舍入误差越大稳定性问题：指的是误差的积累是否受到控制的问题 误差（计算值－实际值）逐步减小，则计算稳定； 否则，计算就不稳定稳定性分析,,3.3,1. 稳定性测试方程,对高阶系统的数值计算稳定性难以分析，采用一个一阶方程来分析显然，若 ，则系统的解析解是稳定的，,(3-23),该系统的解析解为,(3-24),稳定性分析,,3.3,2. 欧拉法计算稳定性分析,欧拉法：,设xk步的误差为εk ，则有,与前式相减，得,很显然，要使误差是逐渐减小得，就要求,(3-25),如果λ0，系统的解析解就不稳定，欧拉法计算总是不稳定。

即,稳定性分析,,3.3,3. RK-4法计算稳定性分析【不讲】,,稳定性分析,,3.3,所以，RK-4法计算稳定的条件是,(3-26),可得前后两步误差之间的关系为,稳定性分析,,3.3,【例3.3.1】针对下面的系统，用解析法、欧拉法和RK4分别求解，计算欧拉法最大允许步长，将步长从0.1逐渐增大，比较三种解的效果解：1）该系统是稳定的，解析解为,2）用欧拉法计算本例时，其步长应该满足(3-25)式，其中λ=-4,3）RK4步长条件式(3-26)是一个高阶不等式，无法直接求解，只能用试探法确定RK4的步长上限3-27),稳定性分析,,3.3,4）下面是步长从0.1不断增大时的仿真结果,1】RK4法的曲线与解析解重合，说明RK4法非常准确 2】欧拉法有较大误差h=0.1时的仿真曲线,稳定性分析,,3.3,1】RK4法还是比较准确 2】欧拉法误差很大，出现衰减震荡 3】两种数值算法的步长条件仍然满足，算法仍然稳定h=0.4时的仿真曲线稳定性分析,,3.3,欧拉法的步长条件已经不满足，仿真解发散,h=0.6时的仿真曲线,稳定性分析,,3.3,RK4法的步长条件仍然满足但误差增大,h=0.6时的仿真曲线,稳定性分析,,3.3,RK4法的步长条件已经不满足，仿真结果发散。

h=0.7时的仿真曲线,数值积分仿真程序,,3.4,3.4.1 编程思路,目的：针对下面的系统编写通用的数值积分法仿真程序,(3-28),对这样的系统进行仿真，实际上涉及到控制的计算、状态的数值积分计算和输出的计算3个函数g,f,h确定后，就可以完整地描述一个系统给定初值：u0，x0，向量都采用列向量表示数值积分仿真程序,,3.4,function [t,y]=w_DigiInteSimu(tstart,tstop,h,x0,u0,cnty,InteMethod,StateModel, OutputFile, ControlFile) % 函数功能：用数值积分法仿真 % 输入参数：tstart, tstop,h 分别是起始时间、结束时间和仿真步长，是标量 % x0,u0是状态、输入的初值，都是列向量 % cnty是输出变量的个数 % InteMethod是数值积分方法，可以是'EUL' ,'RK2','RK4' % StateModel是状态模型的文件名 % ControlFile是控制作用的文件名 % OutputFile是系统输出的文件名 % 输出参数：t是仿真结果的时间序列 % y是仿真结果系统的输出序列,(1)数值积分法仿真框架函数,数值积分仿真程序,,3.4,function [t,y]=w_DigiInteSimu(tstart,tstop,h,x0,u0,cnty,InteMethod,StateModel,OutputFile, ControlFile) t=[tstart:h:tstop];%t数一个行序列 cntt=size(t,2);%返回列数 y=zeros(cnty,cntt);%构造一个空矩阵，用来存储结果 y0=eval([OutputFile,'(tstart,x0,u0)']);%计算初始输出 y(:,1)=y0;%将cury作为输出的第1列 curx=x0; %当前一步的x curu=u0; %当前一步的u cury=y0; %当前一步的y for i=1:1:cntt-1 curu=eval([ControlFile,'(t(i),h,curx,curu,cury)']);%计算控制时传递的参数：当前时间，步长，当前状态和输出 curx=w_StepIntegral(t(i),h,curx,curu,InteMethod,StateModel);%单步积分运算 cury=eval([OutputFile,'(t(i),curx,curu)']);%计算输出 y(:,i+1)=cury;%将输出加入到输出序列里 end,数值积分仿真程序,,3.4,function NewX=w_StepIntegral(curt,h,curx,curu,InteMethod,StateModel) %函数功能：单步积分运算,与模型方程有关 % 输入参数：curt是当前时间，h是数值积分步长 % curx,curu分别是当前的状态和控制向量 % InteMethod是积分算法，字符串类型，可以取'EUL','RK2','RK4' % StateModel是状态模型文件名称，字符串类型 % 输出参数：NewX是这一步计算的新的状态向量,(2)单步数值积分法函数,单步数值积分函数只是对微分方程组StateModel进行一步的计算，计算法由InteMethod参数指定，可以上欧拉法，RK2或RK4。

数值积分仿真程序,,3.4,function NewX=w_StepIntegral(curt,h,curx,curu,InteMethod,StateModel); if InteMethod == 'RK4' k1=eval([StateModel,'(curt,curx,curu)']); k2=eval([StateModel,'(curt+0.5*h,curx+0.5*h*k1,curu)']); k3=eval([StateModel,'(curt+0.5*h,curx+0.5*h*k2,curu)']); k4=eval([StateModel,'(curt+h,curx+h*k3,curu)']); NewX=curx+h*(k1+2*k2+2*k3+k4)/6; elseif InteMethod == 'RK2' k1=eval([StateModel,'(curt,curx,curu)']); k2=eval([StateModel,'(curt+h,curx+h*k1,curu)']); NewX=curx+0.5*h*(k1+k2); else %欧拉法EUL Qk=eval([StateModel,'(curt,curx,curu)']);%eval用来执行一个函数，传递函数名和函数的输入参数 NewX=curx+h*Qk; end,数值积分仿真程序,,3.4,函数w_DigiInteSimu和w_StepIntegral构造了一个数值积分法仿真的框架，并不涉及具体的系统。

具体的系统由StateModel,ControlFile,OutputFile参个参数决定，实际上就是三个函数文件名，这三个函数输入输出参数必须遵循特定的格式，在准备好由这3个函数描述的系统后，调用w_DigiInteSimu即可进行仿真 还需要一个调用w_DigiInteSimu进行仿真的程序，它指定模型文件，指定初始参数，并且对仿真结果绘图数值积分仿真程序,,3.4,3.4.2 例3.3.1的编程实现,在对例3.3.1进行编程时，我们使用3.4.1介绍的两个基本函数，同时要针对系统模型编写3个函数来描述该系统，最后编写一个实现仿真初值设置、仿真求解和仿真结果显示的函数function derX=w_SysStateEqu(curt,curx,curu) % function derX=w_stateEqu(curt,curx,curu) % 函数功能：在此函数里写系统的状态方程 % 输入参数：curt当前时间，curx和curu是当前状态和控制 % 输出参数：derX是状态的X的导数值 derX(1)=-4*curx(1);,（1）状态模型函数w_SysStateEqu 在该函数里描述(3-27)式的状态微分方程,数值积分仿真程序,,3.4,function OutputY=w_SysOutput(curt,curx,curu) % function OutputY=w_SysOutput(curt,curx,curu) % 函数功能：计算系统的输出向量 % 输入参数：curt是当前时间，curx,curu分别是当前的状态和控制向量 % 输出参数：OutputY是计算出来的输出向量 OutputY(1)=curx(1);%根据具体的模型写输出方程,（2）系统输出函数w_SysOutput,function ControlU=w_SysControl(curt,h,curx,curu,cury) % ControlU=w_sysControl(curt,h,curx,curu,cury) % 函数功能： 计算系统的控制，如u=PID(t,curx,cury) % 输入参数：curt是当前时间，h是数值积分步长 % curx,curu,cury分别是当前的状态、控制和输出向量 % 输出参数：ControlU是计算出来的控制向量 ControlU=curu; %根据实际系统写控制,（3）系统控制计算函数w_SysControl,数值积分仿真程序,,3.4,（4）仿真组织函数,function simu_Stability % 函数功能：稳定性分析仿真 tstart=0;tstop=7;h=0.7; StateModel='w_SysStateEqu'; ControlFile='w_SysControl'; OutputFile='w_SysOutput'; x0=[1];u0=[0];cnty=1; [t,y1]=w_DigiInteSimu(tstart,tstop,h,x0,u0,cnty,'EUL',StateModel,OutputFile,ControlFile); %用欧拉法求解 [t,y4]=w_DigiInteSimu(tstart,tstop,h,x0,u0,cnty,'RK4',StateModel,OutputFile,ControlFile); %用RK4求解 y5=exp(-4*t);% 计算解析解 plot(t,y1,'--m',t,y4,'-.r',t,y5,'-k');%绘制叠加曲线 xlabel('时间'); ylabel('y'); title('算法稳定性分析仿真');,该函数的主要任务是： 1) 设置仿真起始时间、结束时间、仿真步长、初始的x和u 2) 设置模型的输出方程、状态方程和控制方程所在函数文件名 3) 调用w_DigiInteSimu进行仿真计算 4) 将计算结果绘图进行比较。

数值积分仿真程序,,3.4,3.4.3 卫星发射仿真,卫星发射运动方程,G为重力系数,系统是高阶微分方程，不能直接用于仿真计算，需要转化为状态方程定义状态。