第七章系统仿真的MATLAB实现由于计算机技术的高速发展，我们可以借助计算机完成系统的数字仿真。
综前所述，数字仿真实质上是根据被研究的真实系统的模型，利用计算机进行实验研究的一种方法。
仿真的主要过程是：建立模型、仿真运行和分析研究仿真结果。
仿真运行就是借助一定的算法，获得系统的有关信息。
MATLAB是一种面向科学与工程计算的高级语言，它集科学计算、自动控制、信号处理、神经网络和图像处理等学科的处理功能于一体，具有极高的编程效率。
MATLAB是一个高度集成的系统，MATLAB提供的Simulink是一个用来对动态系统进行建模、仿真和分析的软件包，它支持线性和非线性系统，能够在连续时间域、离散时间域或者两者的混合时间域里进行建模，它同样支持具有多种采样速率的系统。
在过去几年里，Simulink已经成为数学和工业应用中对动态系统进行建模时使用得最为广泛的软件包。
MATLAB仿真有两种途径：（1）MATLAB可以在SIMULINK窗口上进行面向系统结构方框图的系统仿真；（2）用户可以在MATLAB的COMMAND窗口下，用运行m文件，调用指令和各种用于系统仿真的函数，进行系统仿真。
这两种方式可解决任意复杂系统的动态仿真问题，前者编辑灵活，而后者直观性强，实现可视化编辑。
下面介绍在MATLAB上实现几类基本仿真。
7.1 计算机仿真的步骤在学习计算机仿真以前，让我们先总结一下计算机仿真的步骤。
计算机仿真，概括地说是一个“建模—实验—分析”的过程，即仿真不单纯是对模型的实验，还包括从建模到实验再到分析的全过程。
因此进行一次完整的计算机仿真应包括以下步骤：（1）列举并列项目每一项研究都应从说明问题开始，问题由决策者提供或由熟悉问题的分析者提供。
（2）设置目标及完整的项目计划目标表示仿真要回答的问题、系统方案的说明。
项目计划包括人数、研究费用以及每一阶段工作所需时间。
（3）建立模型和收集数据模型和实际系统没有必要一一对应，模型只需描述实际系统的本质或者描述系统中所研究部分的本质。
因此，最好从简单的模型开始，然后进一步建立更复杂的模型。
（4）编制程序和验证利用数学公式、逻辑公式和算法等来表示实际系统的内部状态和输入/输出的关系。
建模者必须决定是采用通用语言如MATLAB、FORTRAN、C还是专用仿真语言来编制程序。
在本教材中，我们选择的是MATLAB和其动态仿真工具Simulink。
（5）确认确认指确定模型是否精确地代表实际系统。
它不是一次完成，而是比较模型和实际系统特性的差异，不断对模型进行校正的迭代过程。
（6）实验设计确定仿真的方案、初始化周期的长度、仿真运行的长度以及每次运行的重复次数。
（7）生产性运行和分析通常用于估计被仿真系统设计的性能量度。
利用理论定性分析、经验定性分析或系统历史数据定量分析来检验模型的正确性，利用灵敏度分析等手段来检验模型的稳定性。
（8）文件清单和报表结果（9）实现图7.1是计算机仿真的程序图。
图7.1 计算机仿真程序流图7.2 基于数值积分法的连续系统仿真7.2.1 数值积分法的MATLAB实现MATLAB的工具箱提供了各种数值积分方法函数，这些函数是ODE23、ODE45、ODE113和ODE15s。
这些函数均是m文件，还有一个函数是ode1.C，是直接用C语言编写的。
函数ode23( )是用Runge-Kutta法求解微分方程。
它是一种采用三阶积分算法、二阶误差估计、变积分步长的低阶积分算法，调用格式为[T, Y] = ode23 ( 'F', TSPAN, YO, OPTIONS )其中，F为系统模型文件名，模型为y' = f( t, y )形式；TSPAN = [ T o TFINAL] 为积分计算时间，初值为T o ，终值为TFINAL ； YO 为系统输出初始值；OPTIONS 选项积分计算相对允差 'RelTol' 和绝对允差 'AbsTol'，当缺省时，Reltol ＝1e-3, AbsTol ＝1e-6T 为计算点时间向量，Y 为微分方程的解。
函数ode45( )也是用Runge-Kutta 法求解微分方程，它是变步长的一种中等阶次积分算法，调用格式为[T, Y] = ode45 ( 'F' , TSPAN, YO, OPTIONS )各项含义同上。
函数ode113( )是变阶的Adams-Bashforth-Moulton ，用变阶方法解微分方程，采用多步法，调用格式为[T, Y] = ode113 ( 'F', TSPAN, YO, OPTIONS )各项含义同上。
函数odel5s( )采用改进的Gear 法解微分方程，调用格式为 [T, Y] = odel5s ( 'F', TSPAN, YO, OPTIONS ) 各项含义同上。
MATLAB 还提供了解微分方程函数ode23s 。
函数ode1.C 是用Euler 法求解微分方程。
在SIMULINK 中，系统仿真有变步距和固定步距两种工作方式。
在变步距仿真中，可选用ode45、ode23、ode113、ose15s 和ode23s 。
在固定步距仿真中，可选用材ode5、ode4、ode3、ode2和ode1。
ode5是5阶Runge-Kutta 法，ode4是4阶Runge-Kutta 法，odel 是Euler 法。
【例7.1】求微分方程1 x 10,t 0 5,x x0=≤≤+= 。
用MATLAB 编写程序： 建立一个m 函数文件dfun.mfuncion y ＝dhfn ( t ，x ) y ＝sqrt ( 3*x )＋5 ;在MATLAB COMMAND 窗口下运行 % MATLAB PROGRAM 7-1 [ t, x ]＝ode23 ( 'dfun'，[0 l0]，1) ; Plot ( t ，x ) ;7.2.2 基于数值积分法的连续系统的数字仿真对于一个n 阶微分方程表示的连续系统，也可以用具有n 个状态变量的状态空间表达式来描述：DU CX YBU AX X+=+=状态方程BU AX X += 实际上是n 个一阶微分方程组成的方程组。
如果应用四阶Runge-Kutta 法进行仿真计算。
则该式可以改写为h)BU(t )hK A(X K ))2h BU(t K 2h A(X K )2h BU(t )K 2h A(X K BU AX K )K 2K 2K (K 6hX X k 3k 4k 2k 3k 1k 2k k 14321k 1k +++=+++=+++=+=++++=+式中，K 1，K 2，K 3，K 4均为n 维列向量。
系统输出为y k ＝CX k ＋DU( t k )类似地，可将各种数值积分方法用于连续系统的仿真中。
【例7.2】用MATLAB编写图7.2所示液压控制系统的仿真程序。
图7.2 系统方框图用MATLAB编写仿真程序，采用四阶Runge-Kutta法：% MATLAB PROGRAM 7-2% ****** Input system data ***** %调入数据文件hynat;% Input system function; %调入系统模型ypfun1 = ' valve ';ypfun = ' hysys ';% Initialization %初始化yref = 5000; % Referent value of system outputx0 = [0 0]; % Initial value of servo valvey0 = [0 0 0]; % Initial value of hydraulic cylinderu0 = 0;t0 = 0; % Start time of simulationtfinal = 1; % End time of simulationtsamp = 0.01; % Sample periodh = 0.001; % Simulation step sizeY1imit = 512;Ilimit = 40;max_epoch = fix ( tfinal / h )－1;t = t0;u1 = u0;x = x0';y = y0;tout = zeros ( max_epoch, 1 );uout = zeros ( max_epoch, 1);yd = zeros ( max_epoch, 1 );yout = zeros ( max_epoch, legth ( y ) );i = 1;tout ( i ) = t;uout ( i ) = 1;yd ( i ) = Kf * y ( l );yout ( i, : ) = y';% The main loopfor i = 1 : max_epoch% Compute output of valvesvl = feval ( ypfun1, t, u, x, av, bv );sv2 = fcval (ypfun1, t+h/2, u, x+h*sv1/2, av, hv ); %模块valve仿真sv3 = feval (ypfun1, t+h/2, u, x+h*sv2/2, av, hv );sv4 = feval (ypfun1, t+h, u, x+h*sv3, av, hv);x = x+h*(svl+2*sv2+2*sv3+sv4)/6;vo = cv*x;% Compute output of cylinderxl = feval (ypfun, t, vo, y, a, b);s2 = feval (ypfun, t+h/2, vo, y+h*s1/2, a, h); %模块hysya 仿真s3 = feval (ypfun, t+h/2, vo, y+h*s2/2, a, b);s4 = feval (ypfun, t+h u, vo, y+h*s3, a, b);y = y+h* (sl+2*s2+2*s3+s4)/6;i=i+lt=t+h;tout(i) =t;uout(i) =u;yd(i)=Kf*y(l);yout(i, :) = y';% Discrete control process %离散控制量计算if abs (round (t/tsamp)-t/tsamp)<le－9ye = yref－y(l)*Kf;ul=Kd*ye;% Saturation blockif ul>Ylimitxl=ylimit ;else if ul<－ylimitxl =－Ylimit;else xl = ul;endend% D/A conversionx2 = Kda*xl;% Amplifier gainx3=Ka*x2;% V/I conversionu4=Kvi*x3;% limit of currentif u4>Ilimitx4=Ilimit;else if u4<－Ilimitx4=－Ilimit;else x4 =u4;endendu = Kq*x4;end % for discrete section end % for main loop% save data to file %存储仿真数据hout = [tout uout yout];save hout.dat hout -ascii;plot (tout, yd, ' y');gird;m 函数文件valve .mfunction xd = valve (t, u, x, av, bv)xd = av*x+hv*u;m 函数文件hysys. Mfunction yd=hysys(t, u, y, a, b)yd = a* y+h*u上例的仿真结果如图7.3所示。