系统时间响应及其仿真概述

第七章系统时间响应及其仿真本章主要内容以下：7.1仿真算法

7.2系统仿真MATLAB函数系统时间响应是指系统在输入信号或初始状态作用下，系统输出随时间改变情况。系统时间响应反应了系统特征和性能，如快速性、稳定性等。对系统时间响应分析是我们设计、校正系统基础。

系统时间响应及其仿真概述第1页第七章系统时间响应及其仿真7.1仿真算法对系统时间响应进行动态仿真，采取什么样仿真算法是一个至关主要问题。对连续时间系统进行数字动态仿真，主要是两种方法：基于数值积分仿真方法;基于离散相同法仿真方法。因为后者包括到离散控制系统理论，所以本节重点介绍基于数字积分连续系统仿真方法。系统时间响应及其仿真概述第2页7.1仿真算法7.1.1数值求解基本概念设微分方程为则求解方程中函数y(t)问题，就是已知初值常微分方程求解问题。所谓数值求解就是要在时间区间[a,b]中取若干离散点

求出微分方程在这些时刻近似值常微分方程数值求解基本方法是数值积分法。系统时间响应及其仿真概述第3页7.1仿真算法7.1.2数值积分基本原理积分区间划分将区间[a,b]分成N个小区间,时间间隔h（)也称为积分步长,在第k个间隔t=[tk,tk+1]内积分：则可用yk（k=0,1,…N）作为解y(t)近似值，如图所表示。abtky0ykyt数值积分图解tk+1系统时间响应及其仿真概述第4页7.1仿真算法7.1.2数值积分基本原理数值积分展开式为防止(2)式中积分项，将y在tk

，以h为增量展开成Taylor级数：式(3)是一个递推公式。积分值与实际微分方程解误差取决于步长h和计算所用阶数，它是数值积分基础。系统时间响应及其仿真概述第5页7.1仿真算法7.1.2数值积分基本原理相关概念单步法和多步法单步法指计算yk+1值只需利用tk时刻信息，也称为自开启算法；多步法在计算yk+1值时，则需利用tk，tk-1，…时刻信息。显示法和隐式法显示法在计算yk+1时所需数据均已算出；隐式法在计算yk+1时需用到tk+1时刻数据，该算法必须借助予估公式。定步长和变步长定步长为积分步长在仿真运行过程中一直不变；变步长指在仿真运行过程中自动修改步长。系统时间响应及其仿真概述第6页7.1仿真算法7.1.3数值积分几个算法欧拉算法在(3)式中取前两项：可得欧拉算法：t0t1t2t3hy(t)y0y1yt欧拉近似解欧拉法【说明】欧拉法是用一条过各点切线取代曲线来迫近准确解。该算法简单，计算量小，但精度较低。系统时间响应及其仿真概述第7页7.1仿真算法7.1.3数值积分几个算法梯度法梯度法是欧拉法改进。与欧拉法相比，梯度法是用两个点(tm，ym)、(tm+1,ym+1)斜率平均值来确定下一点y值。因为上式计算时需要用到ym+1值，而ym+1不能预先知道，故梯度法需要和欧拉法结合使用,即用欧拉法对ym+1

进行予估，再由梯度法计算ym+1系统时间响应及其仿真概述第8页7.1仿真算法7.1.3数值积分几个算法龙格-库塔法龙格-库塔法基本思想欧拉算法精度较低，主要是其微分方程解y

Taylor展开式所取项数太少。显然为了提升计算精度，应该取泰勒公式（3）更高阶项。即使增加高阶项可提升计算精度，但也同时带来了需要计算高阶导数困难。龙格-库塔法关键是利用低阶导数组成曲线去拟合含有高阶导数曲线，从而防止了计算高阶导数问题。系统时间响应及其仿真概述第9页7.1仿真算法7.1.3数值积分几个算法龙格-库塔法二阶龙格-库塔（RK）法取(3)式前三项，则有设原微分方程(1)式解含有以下形式:式中，a1，a2，b1，b2为待定系数。系统时间响应及其仿真概述第10页将(8)式中K2按二元函数展开成泰勒级数，并取前三项将K1，K2代入(8)式：比较(6-10)、(6-7)式：系统时间响应及其仿真概述第11页显然由(11)式并不能唯一确定a1,a2,b1,b2，因为只有三个方程。所以对于同一个算法能够有不一样表现形式。【说明】因为该算法只取到泰勒展开式二阶导数项，所以称为二阶龙格-库塔法。但由(8)~(12)式可知，算法并没有用y二阶导数。若设a1=a2，则即二阶RK法公式为系统时间响应及其仿真概述第12页7.1仿真算法7.1.3数值积分几个算法龙格-库塔法龙格-库塔（RK）法普通形式

式中，

i为待定权系数，ai,bij为待定系数，r为使用Ki个数(即级数)，Ki为所取各点导数f值。

Ki个数与yk+1泰勒展开式所取项数相关(即RK算法阶数)，同时还与计算区间内所取导数值点数相关。系统时间响应及其仿真概述第13页7.1仿真算法7.1.3数值积分几个算法龙格-库塔法四阶RK公式四阶RK公式用到了y泰勒展开式四阶导数。在RK算法普通公式（13）中，取r=4可得：因为(14)式在同级RK算法中，计算精度较高，计算量较少，而在系统仿真数值积分中应用十分广泛。称之为四阶四级RK公式。系统时间响应及其仿真概述第14页7.1仿真算法7.1.3数值积分几个算法Gear算法“病态”常微分方程(刚性方程)系数矩阵A特征值含有以下特征：则称为“病态”方程。系统时间响应及其仿真概述第15页7.1仿真算法7.1.3数值积分几个算法Gear算法控制系统仿真中“病态”问题病态系统中绝对值最大特征值对应于系统动态性能解中瞬态分量衰减最快部分，它反应了系统动态响应和系统反应灵敏度。普通与系统中含有最小时间常数Tmin步骤相关，要求计算步长h取得很小。病态系统中绝对值最小特征值对应于系统动态性能解中瞬态分量衰减最慢部分,它决定了整个系统动态过渡过程时间长短。普通与系统中含有最大时间常数Tmax步骤相关，要求计算步长h取得很大。对于病态问题仿真需要寻求愈加合理算法，以处理病态系统带来选取计算步长与计算精度、计算时间之间矛盾。系统时间响应及其仿真概述第16页7.1仿真算法7.1.3数值积分几个算法Gear算法Gear算法Gear算法适合用于病态系统仿真，该算法类似于四阶RK算法系统时间响应及其仿真概述第17页7.1仿真算法7.1.4数值积分方法选择在选择积分方法时应考虑以下几个问题。计算精度数值积分方法所得到离散数值解只是准确解近似，其误差来自两个方面，即舍入误差和局部截断误差。舍入误差：由计算机字长有限而造成计算时舍入误差，它随计算次数增加而增加。所以舍入误差与计算步长

h

成反比。局部截断误差：由积分方法和阶次限制而引发误差。这种误差与h成正比。截断误差舍入误差总误差eh

误差与积分步长

显然选择一个适当积分步长可使总误差到达最小。系统时间响应及其仿真概述第18页7.1.4数值积分方法选择积分步长选择和控制积分步长选择标准在确保数值积分稳定性和精度前提下，尽可能选则较大积分步长，以降低仿真计算次数和仿真时间。固定步长与变步长固定步长：在整个仿真计算过程中，积分步长h一直不变。其算法简单，但极难确保步长最优。另外，h还应与模型信号响应情况相关，比如在稳态时，可取较大步长，见上图。变步长：在仿真计算过程中依据计算误差大小来改变步长。其目标是在确保一定计算精度前提下，尽可能选择较大步长。系统时间响应及其仿真概述第19页7.2系统仿真MATLAB函数7.2.1数值积分方法MATLAB函数对于用数值方法求解常系数微分方程(OrdinaryDifferentialEquation,简写为ODE)或微分方程组，MATLAB提供了七种解函数，最惯用是ODE45（四阶RK算法，单步、变步长，用五阶RK算法估算局部截断误差），其调用格式为：[T,Y]=ode45(‘f’,tspan,y0)【说明1】‘f’为常微分方程（组）或系统模型文件名；

tspan=[t0,tfinal]即积分时间初值和终值；

y0是积分初值；

T为计算时间点时间向量；

Y为对应微分方程解数据向量或矩阵。系统时间响应及其仿真概述第20页7.2系统仿真MATLAB函数7.2.1数值积分方法MATLAB函数【说明2】对于刚性微分方程（特征值数值相差较大），可用ode15s，其调用格式与ode45相同。ode函数只能用于求解一阶微分方程或一阶微分方程组。若系统数学模型为高阶微分方程，则应将高阶微分方程转化成一阶微分方程组。所以在用MATLABode函数求解微分方程时，应首先建立描述系统模型一阶微分方程（组）函数‘f’。系统时间响应及其仿真概述第21页【例2】已知二阶微分方程求时间区间t=[0，20]微分方程解。解：（1）将微分方程表示为一阶微分方程组【说明】这种描述系统微分方程函数与ODE函数配套使用，其格式是固定。dy为2*1数组，其维数等于微分方程阶数。（2）建立描述系统微分方程m-函数文件vdp.mfunctiondy=vdp(t,y)dy=zeros(2,1);%生成2行1列零阵dy(1)=y(2);%dy(2)=(1-y(1)^2)*y(2)-y(1);%系统时间响应及其仿真概述第22页（3）编写MATLAB主程序[T,Y]=ode45('VDPd',[020],[0,1]);%调用ode45产生离散点时间向量和解向量plot(T,Y(:,1),'r-',T,Y(:,2),'b:')title('Solution')xlabel('times'),ylabel('PositionY')legend('y1','y2')运行结果如右图所表示。其中y1(红线)为微分方程解。系统时间响应及其仿真概述第23页7.2系统仿真MATLAB函数7.2.2时间响应仿真MATLAB函数对于线性时不变系统，MATLAB直接提供了在各种输入作用下时间响应函数。阶跃响应仿真函数（STEP）基本调用格式对于LTI连续(或离散)时间系统，以下调用格式可用于绘制系统单位阶跃响应曲线。step(sys)step(sys,Tfinal)其中，sys为系统模型（传递函数）；Tfinal为仿真终止时间，若省略则由系统默认。系统时间响应及其仿真概述第24页【例4】已知系统模型，求其单位阶跃响应。sys=tf([1,-1],[1,1,5])subplot(1,2,1),step(sys,20)subplot(1,2,2),step(sys)建立系统模型指定阶跃响应时间不指定阶跃响应时间系统时间响应及其仿真概述第25页7.2.2时间响应仿真MATLAB函数阶跃响应仿真函数（STEP）多系统阶跃响应调用格式在同一幅图中绘制多个系统单位阶跃响应曲线，可用以下调用格式：这种调用格式，还可定义每个系统响应曲线颜色、线型和标志，比如返回仿真输出调用格式其中，Y为输出响应，T为仿真时间向量。这种调用格式不绘制仿真曲线图。step(sys1,sys2,…)step(sys1,'r',sys2,'y--',sys3,'gx')

[Y,T]=step(sys)系统时间响应及其仿真概述第26页7.2系统仿真MATLAB函数7.2.2时间响应仿真MATLAB函数脉冲响应仿真函数（IMPULSE）IMPULSE函数用来计算LTI系统单位脉冲响应。其调用格式与STEP函数相同。impulse(sys)impulse(sys,Tfinal)impulse(sys1,sys2,…)[Y,T]=impulse(sys)系统时间响应及其仿真概述第27页7.2.2时间响应仿真MATLAB函数信号发生器和任意输入响应函数MATLAB也可计算LTI系统在任意输入作用下时间响应。信号发生器函数GENSIGGENSIG可为系统时间响应产生周期输入信号，其调用格式为：其中,Type为信号类型：‘sin’—正弦波

‘square’—方波

'pulse'—周期脉冲波Tau为信号周期；U为信号值向量;T为与U对应时间向量;Tf为信号时间区间;Ts为采样周期。[U,T]=gensig(Type,Tau)[U,T]=gensig(Type,Tau,Tf,Ts)系统时间响应及