1、第七章 系统时间响应及其仿真仿真技术第七章 系统时间响应及其仿真本章主要内容如下：7.1 仿真算法7.2 系统仿真的MATLAB函数 系统的时间响应是指系统在输入信号或初始状态作用下，系统输出随时间变化的情况。系统的时间响应反映了系统的特征和性能，如快速性、稳定性等。对系统时间响应的分析是我们设计、校正系统的基础。第七章 系统时间响应及其仿真仿真技术第七章 系统时间响应及其仿真7.1 仿真算法n对系统的时间响应进行动态仿真，采用什么样的仿真算法是一个至关重要的问题。对连续时间系统进行数字动态仿真，主要是两种方法：基于数值积分的仿真方法;基于离散相似法的仿真方法。由于后者涉及到离散控制系统理论，
2、因此本节重点介绍基于数字积分的连续系统仿真方法。第七章 系统时间响应及其仿真仿真技术7.1 仿真算法7.1.1 数值求解的基本概念n设微分方程为 则求解方程中函数y(t)问题，就是已知初值的常微分方程求解问题。n所谓数值求解就是要在时间区间a,b中取若干离散点 求出微分方程在这些时刻的近似值n常微分方程数值求解的基本方法是数值积分法。00(,)()dyf t ydty ty01Natttb012Ny y yy()0,1,2kkyy tkN第七章 系统时间响应及其仿真仿真技术7.1 仿真算法7.1.2 数值积分的基本原理n积分区间的划分)1()(),(,(0ytybtatytfdtdyo将区间a
3、,b分成 N个小区间,时间间隔h（)也称为积分步长,在第k个间隔 t=tk,tk+1内积分：Nabh)2(),(11kkttkkdtytfyy则可用yk（k=0,1,N）作为解y(t)的近似值，如图所示。abtky0ykyt数值积分图解tk+1第七章 系统时间响应及其仿真仿真技术7.1 仿真算法7.1.2 数值积分的基本原理n数值积分的展开式 为避免(2)式中的积分项，将y在tk，以h为增量展开成Taylor级数：式(3)是一个递推公式。积分值与实际微分方程解的误差取决于步长h和计算所用的阶数，它是数值积分的基础。),2,1,0()3(3)3(!312)2(!211Nkhyhyhyyykkkk
4、k)2(),(11kkttkkdtytfyy第七章 系统时间响应及其仿真仿真技术7.1 仿真算法7.1.2 数值积分的基本原理n有关概念n单步法和多步法 单步法指计算yk+1值只需利用tk时刻的信息，也称为自启动算法；多步法在计算yk+1值时，则需利用tk，tk-1，时刻的信息。n显示法和隐式法 显示法在计算yk+1时所需数据均已算出；隐式法在计算yk+1时需用到tk+1时刻的数据，该算法必须借助予估公式。n定步长和变步长 定步长为积分步长在仿真运行过程中始终不变；变步长指在仿真运行过程中自动修改步长。第七章 系统时间响应及其仿真仿真技术7.1 仿真算法7.1.3 数值积分的几个算法n欧拉算法
5、 在(3)式中取前两项：可得欧拉算法：t0t1t2t3hy(t)y0y1yt欧拉近似解欧拉法n【说明】欧拉法是用一条过各点的切线取代曲线来逼近精确解。该算法简单，计算量小，但精度较低。)4(1mmmhyyy)5(),(),(),(111120000100hytfyyhytfyyhytfyhyyymmmmty)3(3)3(!312)2(!211kkkkkhyhyhyyy第七章 系统时间响应及其仿真仿真技术7.1 仿真算法7.1.3 数值积分的几个算法n梯度法 梯度法是欧拉法的改进。n与欧拉法相比，梯度法是用两个点(tm，ym)、(tm+1,ym+1)的斜率的平均值来确定下一点的y值。n由于上式计
6、算时需要用到ym+1的值，而ym+1不能预先知道，故梯度法需要和欧拉法结合使用,即用欧拉法对ym+1 进行予估，再由梯度法计算ym+1)3(3)3(!312)2(!211kkkkkhyhyhyyy),(),(11211mmmmmmytfytfhyy)6(),(),()(),(112111（校正）予估mmmmmmmmmmytfytfhyyhytfyy第七章 系统时间响应及其仿真仿真技术7.1 仿真算法7.1.3 数值积分的几个算法n龙格龙格-库塔法库塔法n龙格龙格-库塔法的基本思想库塔法的基本思想 欧拉算法的精度较低，主要是其微分方程解 y 的 Taylor 展开式所取的项数太少。显然为了提高计
7、算精度，应当取泰勒公式（3）更高阶项。虽然增加高阶项可提高计算精度，但也同时带来了需要计算高阶导数的困难。龙格龙格-库塔法库塔法的关键是利用低阶导数构成的曲线去拟合含有高阶导数的曲线，从而避免了计算高阶导数的问题。)3(3)3(!312)2(!211kkkkkhyhyhyyy第七章 系统时间响应及其仿真仿真技术7.1 仿真算法7.1.3 数值积分的几个算法n龙格龙格-库塔法库塔法n二阶龙格二阶龙格-库塔（库塔（RK）法）法 取(3)式的前三项，则有)3(3)3(!312)2(!211kkkkkhyhyhyyy)7(),(),(2212)2(211hytfythfyhyhyyymyfmmmtfm
8、mmmmmm设原微分方程(1)式解具有以下形式:)8()(),(),(2211112121KaKahyyhKbyhbtfKytfKmmmmmm式中，a1，a2，b1，b2为待定系数。),(),()2(ytfdtdyytfdtdy第七章 系统时间响应及其仿真仿真技术)7(),(),(2211hytfythfyymyfmmmtfmmmm)8()(),(),(2211112121KaKahyyhKbyhbtfKytfKmmmmmm 将(8)式中K2按二元函数展开成泰勒级数，并取前三项)9(),(),(1211212myfmtfmmmyfmtfmmKbbhytfhKbhbytfK将K1，K2代入(8)
9、式：)10(),(),()(),(),(),(21222121211myfmmmtfmmmmyfmmmtfmmmmmmytfbbhaythfaayythfbhbytfhaythfayy比较(6-10)、(6-7)式：)11(2/12/11221221babaaa第七章 系统时间响应及其仿真仿真技术 显然由(11)式并不能唯一确定a1,a2,b1,b2，因为只有三个方程。因此对于同一种算法可以有不同的表现形式。n【说明】由于该算法只取到泰勒展开式的二阶导数项，所以称为二阶龙格-库塔法。但由(8)(12)式可知，算法并没有用y的二阶导数。若设a1=a2，则12/12121bbaa即二阶RK法公式为
10、)12()(2),(),(211121KKhyyhKyhtfKytfKmmmmmm)11(2/12/11221221babaaa第七章 系统时间响应及其仿真仿真技术7.1 仿真算法7.1.3 数值积分的几个算法n龙格龙格-库塔法库塔法n龙格龙格-库塔（库塔（RK）法的一般形式）法的一般形式 )13(,2,1),(1111riKbyhatfKKhyyijjijmimiiriimm式中，i为待定权系数，ai,bij为待定系数，r为使用Ki的个数(即级数)，Ki为所取各点导数f的值。Ki的个数与yk+1泰勒展开式所取的项数有关(即RK算法的阶数)，同时还与计算区间内所取导数值的点数有关。第七章 系统
11、时间响应及其仿真仿真技术7.1 仿真算法7.1.3 数值积分的几个算法n龙格龙格-库塔法库塔法n四阶四阶RK公式公式 四阶RK公式用到了y的泰勒展开式的四阶导数。在RK算法的一般公式（13）中，取r=4可得：)14(),(),(),(),(223422123121221432161hKyhtfKhKytfKhKytfKytfKKKKKyymmmhmmhmmmhmm 由于(14)式在同级的RK算法中，计算精度较高，计算量较少，而在系统仿真的数值积分中应用十分广泛。称之为四阶四级RK公式。第七章 系统时间响应及其仿真仿真技术7.1 仿真算法7.1.3 数值积分的几个算法nGear算法n“病态”常微
12、分方程(刚性方程)的系数矩阵A的特征值具有如下特征：00()()(),()X tAX tBU tX tX则称为“病态”方程。)17(50ReminRemax),2,1(0Reiiiini第七章 系统时间响应及其仿真仿真技术7.1 仿真算法7.1.3 数值积分的几个算法nGear算法n控制系统仿真中的“病态”问题a)病态系统中绝对值最大的特征值对应于系统动态性能解中瞬态分量衰减最快的部分，它反映了系统的动态响应和系统的反应灵敏度。一般与系统中具有最小时间常数Tmin的环节有关，要求计算步长h取得很小。b)病态系统中绝对值最小的特征值对应于系统动态性能解中瞬态分量衰减最慢的部分,它决定了整个系统的
13、动态过渡过程时间的长短。一般与系统中具有最大时间常数Tmax的环节有关，要求计算步长h取得很大。c)对于病态问题的仿真需要寻求更加合理的算法，以解决病态系统带来的选取计算步长与计算精度、计算时间之间的矛盾。第七章 系统时间响应及其仿真仿真技术7.1 仿真算法7.1.3 数值积分的几个算法nGear算法nGear算法Gear算法适用于病态系统的仿真，该算法类似于四阶RK算法)18(21121,2112121,22,2,324213121KKyhthfKKKyhthfKKyhthfKythfKmmmmmmmm)19(6121131211316143211KKKKyymm第七章 系统时间响应及其仿真
14、仿真技术7.1 仿真算法7.1.4 数值积分方法的选择 在选择积分方法时应考虑以下几个问题。a)计算精度 数值积分方法所得到的离散数值解只是精确解的近似，其误差来自两个方面，即舍入误差和局部截断误差。n舍入误差舍入误差：由计算机字长有限而造成的计算时的舍入误差，它随计算次数的增加而增加。因此舍入误差与计算步长 h 成反比。n局部截断误差局部截断误差：由积分方法和阶次的限制而引起的误差。这种误差与h成正比。截断误差舍入误差总误差eh 误差与积分步长 显然选择一个合适的积分步长可使总误差达到最小。第七章 系统时间响应及其仿真仿真技术7.1.4 数值积分方法的选择b)积分步长的选择和控制n积分步长的
15、选择原则 在保证数值积分稳定性和精度的前提下，尽可能选则较大的积分步长，以减少仿真计算次数和仿真时间。n固定步长与变步长 固定步长固定步长：在整个仿真计算过程中，积分步长h始终不变。其算法简单，但很难保证步长最优。n此外，h还应与模型的信号响应情况有关，例如在稳态时，可取较大的步长，见上图。6.283212.566418.84960.9511.051.2843 y=1-e-tcost=0.3=0.7 镇定时间 ts=9.6133 t y 变步长：变步长：在仿真计算过程中根据计算误差的大小来改变步长。其目的是在保证一定计算精度的前提下，尽可能选择较大步长。第七章 系统时间响应及其仿真仿真技术7.
16、2 系统仿真的MATLAB函数7.2.1 数值积分方法的MATLAB函数n对于用数值方法求解常系数微分方程(Ordinary Differential Equation,简写为ODE)或微分方程组，MATLAB提供了七种解函数，最常用的是ODE45（四阶RK算法，单步、变步长，用五阶RK算法估算局部截断误差），其调用格式为：T,Y=ode45(f,tspan,y0)n【说明1】f 为常微分方程（组）或系统模型的文件名；tspan=t0,tfinal 即积分时间初值和终值；y0是积分初值；T为计算时间点的时间向量；Y为相应的微分方程解数据向量或矩阵。第七章 系统时间响应及其仿真仿真技术7.2 系
17、统仿真的MATLAB函数7.2.1 数值积分方法的MATLAB函数n【说明2】对于刚性微分方程（特征值数值相差较大），可用ode15s，其调用格式与ode45相同。ode函数只能用于求解一阶微分方程或一阶微分方程组。若系统的数学模型为高阶微分方程，则应将高阶微分方程转化成一阶微分方程组。因此在用MATLAB的ode函数求解微分方程时，应首先建立描述系统模型的一阶微分方程（组）函数f。第七章 系统时间响应及其仿真仿真技术【例2】已知二阶 微分方程 求时间区间t=0，20微分方程的解。1)0(,0)0(0)1(2yyyyyy 解：解：（1）将微分方程表示为一阶微分方程组12212211)1(yyy
18、yyyyyn【说明】这种描述系统微分方程的函数与ODE函数配套使用，其格式是固定的。dy为2*1数组，其维数等于微分方程的阶数。（2）建立描述系统微分方程的m-函数文件vdp.mfunction dy=vdp(t,y)dy=zeros(2,1);%生成2行1列的零阵dy(1)=y(2);%dy(2)=(1-y(1)2)*y(2)-y(1);%21yy 12212)1(yyyy第七章 系统时间响应及其仿真仿真技术（3）编写MATLAB主程序T,Y=ode45(VDPd,0 20,0,1);%调用ode45产生离散点时间向量和解向量plot(T,Y(:,1),r-,T,Y(:,2),b:)titl
19、e(Solution)xlabel(time s),ylabel(Position Y)legend(y1,y2)运行结果如右图所示。其中y1(红线)为微分方程的解。第七章 系统时间响应及其仿真仿真技术7.2 系统仿真的MATLAB函数7.2.2 时间响应仿真的MATLAB函数 对于线性时不变系统，MATLAB直接提供了在各种输入作用下的时间响应函数。n阶跃响应仿真函数（STEP）n基本调用格式 对于LTI连续(或离散)时间系统，以下调用格式可用于绘制系统单位阶跃响应曲线。step(sys)step(sys,Tfinal)其中，sys为系统模型（传递函数）；Tfinal为仿真终止时间，若省略则
20、由系统默认。第七章 系统时间响应及其仿真仿真技术【例4】已知系统模型 ，求其单位阶跃响应。51)(2ssssGsys=tf(1,-1,1,1,5)subplot(1,2,1),step(sys,20)subplot(1,2,2),step(sys)建立系统模型指定阶跃响应时间不指定阶跃响应时间第七章 系统时间响应及其仿真仿真技术7.2.2 时间响应仿真的MATLAB函数n阶跃响应仿真函数（STEP）n多系统阶跃响应调用格式 在同一幅图中绘制多个系统的单位阶跃响应曲线，可用以下调用格式：这种调用格式，还可定义每个系统响应曲线的颜色、线型和标志，例如n返回仿真输出的调用格式 其中，Y为输出响应，T
21、为仿真时间向量。这种调用格式不绘制仿真曲线图。step(sys1,sys2,)step(sys1,r,sys2,y-,sys3,gx)Y,T=step(sys)第七章 系统时间响应及其仿真仿真技术7.2 系统仿真的MATLAB函数7.2.2 时间响应仿真的MATLAB函数n脉冲响应仿真函数（IMPULSE）IMPULSE函数用来计算LTI系统的单位脉冲响应。其调用格式与STEP函数相同。impulse(sys)impulse(sys,Tfinal)impulse(sys1,sys2,)Y,T=impulse(sys)第七章 系统时间响应及其仿真仿真技术7.2.2 时间响应仿真的MATLAB函数
22、n信号发生器和任意输入响应函数MATLAB也可计算LTI系统在任意输入作用下的时间响应。n信号发生器函数 GENSIG GENSIG可为系统时间响应产生周期输入信号，其调用格式为：其中,Type为信号类型：sin正弦波 square方波 pulse周期脉冲波 Tau为信号周期；U为信号值向量;T为与U对应的时间向量;Tf为信号的时间区间;Ts为采样周期。U,T=gensig(Type,Tau)U,T=gensig(Type,Tau,Tf,Ts)第七章 系统时间响应及其仿真仿真技术7.2.2 时间响应仿真的MATLAB函数n信号发生器和任意输入响应函数n 任意输入响应函数 LSIM LSIM用来
23、仿真系统对任意输入的时间响应，并绘制响应曲线。其调用格式为：其中，sys为系统模型；U为输入信号向量；T为和输入对应的时间向量;Ys为响应值向量;Ts为与Ys相对应的时间向量。lsim(sys,U,T)%基本调用格式lsim(sys1,sys2,U,T)%绘制多个系统对同一个任意输入响应曲线Ys,Ts=lsim(sys,U,T)%该格式不绘制响应曲线第七章 系统时间响应及其仿真仿真技术【例7】已知系统模型 ，计算系统在周期为5s的方波信号作用下的响应。sys=tf(3,100,1,10,40,100);u,t=gensig(square,5,10);%产生方波信号数据lsim(sys,r,u,t),hold on%产生方波响应并绘曲线plot(t,u,-.)%在同一坐标系绘方波波形hold offtext(1.3,0.8,输入rightarrow)text(5.4,0.8,leftarrow输出)10040101003)(23sssssG第七章 系统时间响应及其仿真仿真技术练习教材126-127页：1、2、3题
