连续系统仿真方法学PPT学习教案

1、会计学1 连续系统仿真方法学连续系统仿真方法学 第1页/共72页 第2页/共72页 第3页/共72页 模型模型 第4页/共72页 第5页/共72页 11 110110 11 nnmm nmm nnmm d ydydyd udu aaa ybbbub u dtdtdtdtdt 11 110110 11 nnmm nmm nnmm d ydydyd udu aaa ybbbub u dtdtdtdtdt 第6页/共72页 1 110 1 110 ( ) ( ) ( ) mm mm nn n b sbsb sbY s W s U ssasa sa 12 ( )( ) ( )( )( )( ) l Y

2、 sW s U sW sW sW s 1 12 ( ) ( )( )( )( ) l y tLY sy ty ty t 第7页/共72页 权函数权函数 g(t)指初始条件为指初始条件为0时系统在理想脉冲时系统在理想脉冲 函数函数(t)作用下的响应，又称脉冲过渡函数作用下的响应，又称脉冲过渡函数 系统对任意输入的响应可由卷积积分公式求出系统对任意输入的响应可由卷积积分公式求出 0 ( )( ) () t y tug td 权函数与传递函数有如下关系：权函数与传递函数有如下关系： ( )( )L g tG s 第8页/共72页 为为n维状态向量，维状态向量， u为为r维输入向量，维输入向量，y为为

3、m维输出向量维输出向量 动态系统的状态是指能完全描述系统行为的最小一组变量动态系统的状态是指能完全描述系统行为的最小一组变量 若知道若知道t=t0时刻的初始状态向量时刻的初始状态向量x0及及tt0时的输入时的输入u就能完全就能完全 确定系统在确定系统在tt0时刻的行为时刻的行为 状态空间表达式由状态方程和输出方程组成状态空间表达式由状态方程和输出方程组成 xAxBu yCxDu 1 ,T n xxx n n AA 为参数矩阵（或称动态矩阵），为参数矩阵（或称动态矩阵）， n r BB 为输入矩阵，为输入矩阵， m n CC 为输出矩阵为输出矩阵 m r DD 运用矩阵计算方法且借助于计算机很容

4、易运用矩阵计算方法且借助于计算机很容易 对状态空间方程求解对状态空间方程求解 为关联矩阵（输入和输出直接关联）为关联矩阵（输入和输出直接关联） 第9页/共72页 第10页/共72页 110 10 () (1) (1) () () (1) ()() n mm y kn Tay knTa y kTa y kT c u km Tcu kmTc u kTmn 第11页/共72页 1 110 1 110 ( ) ( ) ( ) mm mm nn n b zbzb zbY z W z U zzczc zc 第12页/共72页 0 ( )( ) () k i y ku i h ki ( )( )Z h kW

5、 z 第13页/共72页 (1)( )( ) ( )( ) x kFx kGu k y kx k 第14页/共72页 第15页/共72页 设有微分方程设有微分方程 先考虑右边仅含先考虑右边仅含u的形式，令的形式，令 = 则有则有 11 110110 11 nnmm nmm nnmm d ydydyd udu aaa ybbbub u dtdtdtdtdt 12 , T n xx xx (1) , , nT y yy 12 23 1 01121 nn nnn xx xx xx xa xa xaxu 第16页/共72页 写成矩阵形式：写成矩阵形式： 输出矩阵写为输出矩阵写为 其中其中 xAxBu

6、1 011 0 0 n n I A aaa 0 0 1 B 100yCxC 第17页/共72页 当右式包含导数项时，状态方程形式为当右式包含导数项时，状态方程形式为 A,C与前相同，与前相同， 其中其中 xAxBu yCxDu 10 , T n BD 1 ,(0,1,) i in injij j bain 现代控制理论中还介绍了其它形式的转换方程，现代控制理论中还介绍了其它形式的转换方程， 如能控标准型、能观标准型等如能控标准型、能观标准型等 第18页/共72页 连续状态方程对应的离散状态表达式为连续状态方程对应的离散状态表达式为 T为采样周期或者计算步长为采样周期或者计算步长 为确定为确定(

7、T)和和H(T)，可利用连续状态方程解，可利用连续状态方程解 其中其中 At e为系统的矩阵指数或状态转移矩阵，为系统的矩阵指数或状态转移矩阵， x(0)为初始状态向量为初始状态向量 当采用零阶保持器时，当采用零阶保持器时， 即认为即认为u(t)在每个采样周期内保持常值在每个采样周期内保持常值 u(t)=u(kT)，( kTt(k+1)T ) (1) ( ) ()( ) ()(0,1,)x kTT x kTH T u kTk 0 ( )(0)( ) t AtAtA x te xeeBud 第19页/共72页 (1) (1)(1) 0 0 (1) (1) 0 (1) (0)( )(1) (0)(

8、 )(2) (1)(2) (1) ()( ) ()()( () () kT A kTA kTA kT AkTAkTA AT kT ATA kTA kT T ATATAt ATA x kTexeeBud x kTexeeBud e x kTex kTeeBud ex kTeeBu kT dtu kT ex kTe 与 得： 不变 0 () ( ) ()( ) () T Bu kT dt T x kTH T u kT 则有则有 其中其中和和H与与T有关，当有关，当T确定后，确定后，和和H为常值矩阵为常值矩阵 第20页/共72页 离散化公式的核心在于计算矩阵指数及其积分，离散化公式的核心在于计算矩阵

9、指数及其积分， 常用级数展开的算法，即常用级数展开的算法，即 其余离散化的表示形式与连续形式之间的转换其余离散化的表示形式与连续形式之间的转换 在计算机控制中介绍在计算机控制中介绍 22 11 2! ATkk eIATA TA T k 第21页/共72页 第22页/共72页 10cxm 第23页/共72页 第24页/共72页 若运算过程中计算误差不断增长，称算法为数值不稳定的若运算过程中计算误差不断增长，称算法为数值不稳定的 反之则为稳定的反之则为稳定的 例：计算例：计算 由分部积分得递推公式：由分部积分得递推公式： 用用Taylor展开计算：展开计算： 若取若取k=7，并保留，并保留4位小数

10、，可得位小数，可得 截断误差：截断误差： 1 1 0 (0,1,) nx n Iex e dxn 1 11 10 0 1,1 x nn InIIee dxe 2 1 ( 1)( 1) 1( 1) 2! k e k 1 0.3679e 14 7 11 0.367910 8!4 Re 第25页/共72页 只考虑初值误差，对只考虑初值误差，对I00.6321，递推计算，结果如表中第一行，递推计算，结果如表中第一行 该积分不可能为负值，显然算法有问题该积分不可能为负值，显然算法有问题 分析计算误差：分析计算误差： 满足关系：满足关系： 误差增长迅速误差增长迅速 n0123456789 0.63210.

11、36790.26420.20740.17040.14800.11200.2160-0.7287.552 0.63210.36790.26430.20730.17080.14550.12680.11210.10350.0684 nnn EII 1nn EnE 0 ( 1) ( !) n n En E 第26页/共72页 如果换一种算法，可以减小误差，考虑积分估计值：如果换一种算法，可以减小误差，考虑积分估计值： 逆向计算：逆向计算： 取取n=9时，时， ，结果如表中第二行，结果如表中第二行 分析计算误差，满足关系：分析计算误差，满足关系： 显然误差一直减小显然误差一直减小 1 11 11 000

12、101 1 (min)(max) 11 xnxn n xx e eex dxIeex dx nn * 1 1 (1) nn II n 1 9 11 ()0.0684 2 1010 e I 1 1 nn EE n 第27页/共72页 如线性方程组如线性方程组 精确解：精确解： 12 12 3.0004.12715.41(1) 1.0001.3745.147(2) xx xx 1 2 13.6658 6.2 x x 如果用如果用4位有效数字进行运算：位有效数字进行运算：(2)-(1)/3.000， 逐步计算后，可得逐步计算后，可得x2=-5.000 ，误差很大，误差很大 当方程特征根相差太大时出现

13、病态问题，也称刚性当方程特征根相差太大时出现病态问题，也称刚性(Stiff)问题问题 需要设计有效的算法需要设计有效的算法 第28页/共72页 第29页/共72页 设有模型方程设有模型方程 ，初始条件，初始条件 欧拉法用欧拉法用tk点切线近似该点附近的曲线点切线近似该点附近的曲线f(t,y)，则有则有 其中其中是曲线上的点，是曲线上的点，是切线上的点是切线上的点 称为第称为第k 步的计算步长步的计算步长 此类方法称为此类方法称为“微分方程初值问题的数值计算法微分方程初值问题的数值计算法”， 也称也称“数值积分法数值积分法” 优点：简单易行优点：简单易行 缺点：缺点：h取得大时，计算速度快，单步

14、误差大；取得大时，计算速度快，单步误差大； h取得小，计算速度慢，且累计误差大取得小，计算速度慢，且累计误差大 1 () k y t 1k y ( , ) dy f t y dt 00 ( )y ty 111 ()( ,)() kkkkkkk y tyyf tytt 1kkk htt 第30页/共72页 二阶形式：二阶形式： 其中其中 迭代公式由迭代公式由Taylor展开并保留展开并保留h2项获得项获得 注意：注意：R-K方法实质是用均差代替导数，方法实质是用均差代替导数， 其中其中k项的加权系数可任选项的加权系数可任选 1112 ()() 2 kkk h y tyykk 1 21 ( ,)

15、(,) kk kk kf ty kf th yk h 第31页/共72页 四阶形式（固定步长）：四阶形式（固定步长）： 其中其中 四阶形式在精度和复杂度方面都有较好的表现，最常用四阶形式在精度和复杂度方面都有较好的表现，最常用 Euler法与法与R-K法计算时仅用到前一步的结果，法计算时仅用到前一步的结果， 称单步法，已知初值后可自启动称单步法，已知初值后可自启动 111234 ()(22) 6 kkk h y tyykkkk 1 21 32 43 ( ,) (/2,/2) (/2,/2) (,) kk kk kk kk kf ty kf thyk h kf thyk h kf th yk h

16、 第32页/共72页 Euler法是用矩形公式（面积）近似定积分，在曲线法是用矩形公式（面积）近似定积分，在曲线 与矩形的边之间有较大误差；与矩形的边之间有较大误差； Adams法考虑用梯形公式（面积）代替矩形公式法考虑用梯形公式（面积）代替矩形公式 称二阶隐式称二阶隐式Adams法公式法公式 因公式右端包含未知项，不能直接求解，可以用迭因公式右端包含未知项，不能直接求解，可以用迭 代法求解，设有初值代法求解，设有初值 迭代公式：迭代公式： 直至规定的精度直至规定的精度 1 111 ( , ) (,)( ,)() 22 k k t kkkkkk t hh f t y dtf tyf tyff

17、11 () 2 kkkk h yyff (0) 1k y ( )(1) 111 (,)( ,) 2 nn kkkkkk h yyf tyf ty 第33页/共72页 隐式迭代计算步数太多，为此可降低精度，设计显式隐式迭代计算步数太多，为此可降低精度，设计显式 Adams法，公式：法，公式： Adams法的统一形式为：法的统一形式为： 算法特点：多步法，不能自启动，隐式方法还需迭代求解算法特点：多步法，不能自启动，隐式方法还需迭代求解 实际应用中，先用显式法计算初值，再用隐式法校正一次实际应用中，先用显式法计算初值，再用隐式法校正一次 ，称预报，称预报-校正法校正法 与与R-K法比较，同样阶次和

18、精度下法比较，同样阶次和精度下Adams法计算次数较少法计算次数较少 1111 ()3 ( ,)(,) 2 kkkkkkk h y tyyf tyf ty 111011 kkkkNkN yyh B fB fBf 第34页/共72页 试验方程：试验方程： 若数值积分公式为若数值积分公式为 其中其中 是一个高阶多项式函数是一个高阶多项式函数 仅当仅当 时算法才稳定时算法才稳定 如如Euler法稳定条件为法稳定条件为 隐式一阶、二阶隐式一阶、二阶Adams法恒稳定，更高阶条件稳定法恒稳定，更高阶条件稳定 R-K法正好相反，阶次越高稳定域越大法正好相反，阶次越高稳定域越大 ,Re0 dy yi dt

19、1 () kk yp hy ()p h ()1p h 11h 第35页/共72页 两个原则：保证稳定性，要求一定的计算精度两个原则：保证稳定性，要求一定的计算精度 受稳定性限制，受稳定性限制，h应在系统中最小时间常数量级应在系统中最小时间常数量级 如如R-K4，要求步长小于系统中最小时间常数的，要求步长小于系统中最小时间常数的2.78倍倍 实际应用中对大量的仿真计算可以考虑采用变步长法实际应用中对大量的仿真计算可以考虑采用变步长法 自动改变步长自动改变步长 Matlab中中Simulink仿真时缺省的数值积分法为变步长仿真时缺省的数值积分法为变步长 法法 第36页/共72页 首先进行误差估计：

20、首先进行误差估计： 分别找一个三阶和一个四阶分别找一个三阶和一个四阶R-K公式：公式： 11345 1145 (39126) 30 (4) 6 kk kk h yykkkk h yykkk 其中其中 1 21 312 413 5134 ( ,) (/3,/3) (/3,/6/6) (/2,/83/8) (,/23/22) kk kk kk kk kk kf ty kf thyk h kf thyk hk h kf thyk hk h kf th yk hk hk h 第37页/共72页 则误差为则误差为 1111345 (298) 30 kkk h Eyykkkk 变步长策略：变步长策略： 设

21、定一个最小误差限设定一个最小误差限 ，一个最大误差限，一个最大误差限 每一步的局部误差取为每一步的局部误差取为 ，第，第k+1步有效，下一步用步有效，下一步用2h积分；积分； ，保持，保持h不变；不变； ，第，第k+1步无效，步长变为步无效，步长变为h/2 min max 111 /(| 1) kkk eEy 1mink e min1maxk e 1maxk e 第38页/共72页 形式：形式： 传递函数传递函数Z传递函数：传递函数：Z域离散相似模型，域离散相似模型， 状态空间模型状态空间模型离散状态方程：时域离散相似模型，离散状态方程：时域离散相似模型， 右边第三项表示使用一阶保持器增加的项

22、右边第三项表示使用一阶保持器增加的项 ( )( ) ( ) h G zZ G s G s (1) ( ) ()( ) ()( ) () mm x kTT x kTT u kTT u kT xAxBu 第39页/共72页 由采样定理，为使离散相似模型中重构的信号能由采样定理，为使离散相似模型中重构的信号能 精确表示原信号，应有采样时间小于系统最小时精确表示原信号，应有采样时间小于系统最小时 间常数的一半，或者采样频率是最大信号频率的间常数的一半，或者采样频率是最大信号频率的 两倍两倍 离散相似法的优点：不易受模型方程特性的影响离散相似法的优点：不易受模型方程特性的影响 ，尤其对特征根相差较大的系

23、统十分有效，计算，尤其对特征根相差较大的系统十分有效，计算 速度也更快；但有的模型不易离散化速度也更快；但有的模型不易离散化 第40页/共72页 的计算误差；的计算误差；u(t)在采样间隔中的近似处理在采样间隔中的近似处理 对对后者后者，当输入是典型函数时，可通过增广矩阵，当输入是典型函数时，可通过增广矩阵 法消除误差法消除误差 考虑时域离散相似法，误差来源：考虑时域离散相似法，误差来源： AT e 状态方程状态方程 的解为：的解为： 对典型函数，考虑将输入对典型函数，考虑将输入u(t)增广到状态向量增广到状态向量x(t)中，中， 得得 齐次解，避免积分近似出现的误差齐次解，避免积分近似出现的

24、误差 0 ( )(0)( ) t AtAtA x te xeeBud xAxBu ( )(0) At x te x 第41页/共72页 假设系统为假设系统为n阶，模型阶，模型 xAxBu yCx 0 (0)xx 对阶跃输入对阶跃输入 0 ( )1( )u tUt 定义定义 1( ) ( ) n xtu t ，则，则1( ) 0 n xt 增广后的状态方程及输出方程为增广后的状态方程及输出方程为 11 00 nn xxAB xx 1 0 n x yC x 初始条件初始条件 0 10 (0) (0) n xx xU 第42页/共72页 系统模型同上系统模型同上，对斜坡输入，对斜坡输入 0 ( )u

25、 tU t 定义：定义： 1( ) ( ) n xtu t 210 ( )( )( ) nn xtxtu tU 则则 2( ) 0 n xt ，增广后的状态方程及输出方程为，增广后的状态方程及输出方程为 11 22 0 001 000 nn nn xxAB xx xx 1 2 00 n n x yCx x 初始条件初始条件 0 1 02 (0) (0)0 (0) n n xx x Ux 第43页/共72页 第44页/共72页 其中信号比较环节可在控制器内、外进行其中信号比较环节可在控制器内、外进行 第45页/共72页 第46页/共72页 采样控制系统方块图采样控制系统方块图 其中其中G(s)为

26、被控对象传递函数，为被控对象传递函数，H(s)为保持器为保持器 传递函数，传递函数，D(z)为数字控制器的为数字控制器的z传递函数，传递函数， Ts为实际采样周期为实际采样周期 仿真步距仿真步距T的选择必须根据被控对象结构、采的选择必须根据被控对象结构、采 样周期大小、保持器类型及仿真精度和仿真速样周期大小、保持器类型及仿真精度和仿真速 度的要求综合考虑度的要求综合考虑 第47页/共72页 第48页/共72页 第49页/共72页 第50页/共72页 第51页/共72页 0.98 ( )2.62 0.64 z D z z Ts=0.04s，取Ts*=0.1s仿真，求D(z) 0 0 0 0 0.

27、64 ln()/11.16 ln()/5.505 0.98 p pps zzs z z szT szT z 0.3277 0.9508 s p s z T s p T s z ze ze 0.9508 ( ) 0.3277 z z D zk z 第52页/共72页 再根据稳态值相等原则确定再根据稳态值相等原则确定kz，注意与输入信号有关，注意与输入信号有关 本例中若输入单位阶跃信号，由终值定理本例中若输入单位阶跃信号，由终值定理 1 1 ( )lim( )0.14556 1 z zz uD z zz 同样同样 1 0.9508 lim0.145561.989 0.3277 zz z z kk

28、z 第53页/共72页 第54页/共72页 第55页/共72页 =0 (y,t) y dy L y=0 T(L,t) 第56页/共72页 考虑厚度为考虑厚度为dy的一段扭杆，其力矩平衡方程为：的一段扭杆，其力矩平衡方程为： 22222 222 2 PPPP R dyR I GI GI GdyIdy yyytt 上段力矩下段力矩 其中其中 4 2 P IR为圆截面的极惯性矩为圆截面的极惯性矩 2 2 2 2 t G y 其中其中是杆的线密度，是杆的线密度，G G是材料的剪切弹性模量。是材料的剪切弹性模量。 这是一个物理学经典方程，直接求解非常繁琐。这是一个物理学经典方程，直接求解非常繁琐。 得一

29、维波动方程：得一维波动方程： 第57页/共72页 采用频域分析，先对采用频域分析，先对t进行拉氏变换，有：进行拉氏变换，有： 22 22 ()() tt LL yGt 2 2 2 ( , )( , )( , )( , )y ssy ssy sy s yG 取初始条件为取初始条件为0，则有，则有 22 2 0 ds dyG 解的形式为解的形式为 12 ( , )expexp(1)y scsGycsGy 第58页/共72页 由于杆的固定端有由于杆的固定端有(0,t)=0，则，则(0,s)=0， 可得可得 c1= - c2 (2) 在在y=L处，处，( , )( , ) ( , ), ( , )(3

30、) PP L tL s T L tI GT L sI G yy (1)式对式对y求偏导：求偏导： 12 expexp(4)scsyscsy yGGGG 令令y=L，由，由(2),(3),(4)可求出可求出c1，c2（略），最终得（略），最终得 expexp( , ) ( , )(expexp) sinh()1 cosh() P P sGysGyy s T L sIGssGysGy sGy IGssGy 第59页/共72页 令令s=jw可得可得y=L处角运动的频率响应处角运动的频率响应 sin()( ,)1 (5) ( ,)cos() P LGL j T L jIGLG 显然其中有无限多显然其中

31、有无限多 个固有频率：个固有频率： ,(1,3,5,) 2 k kG k L 若结构以其中某一固有频率振动，则此时的动态若结构以其中某一固有频率振动，则此时的动态 扭转曲线称为振型扭转曲线称为振型 (5)式可用来求振型式可用来求振型 ( ,)( ,)( ,) sin() ( ,)( ,)( ,)2 kkk kkk y jy jT L jky L jT L jL jL 第60页/共72页 4 21 22/ i KKKr G L 扭振系统可用物理集总参数法离散形式（集总块）扭振系统可用物理集总参数法离散形式（集总块） 近似，集总块数目可通过经验或试验的方法确定。近似，集总块数目可通过经验或试验的方

32、法确定。 设为设为2个集总块个集总块 选择每个集总块质心作为集总选择每个集总块质心作为集总 惯量惯量J 所在位置，并以此确定所在位置，并以此确定 K（质心之间杆长的弹性），（质心之间杆长的弹性）， 有有 4 12 /4JJLr 第61页/共72页 原系统可以表示为一原系统可以表示为一 个两段集总模型：个两段集总模型： 2 4 2 4 21 4 1 4 21 4 4 2 )( 4 )( Lr L Gr L Gr Lr L Gr Ti 可求得两个固有频率可求得两个固有频率 1 2 1.53 3.70 G L G L 而准确值为而准确值为1.57和和4.71 当集总数量增大时，可预测出更多的固有频率，当集总数量增大时，可预测出更多的固有频率， 数值也更精确数值也更精确 第62页/共72页 对应物理离散法有数学离散法，如有限差分法对应物理离散法有数学离散法，如有限差分法 ，其中常用的为中心差分法，其中常用的为中心差分法 设设y=f(x,t)，当，当t为常值时，中心差分法即用切点为常值时，中心差分法即用切点 Pn上的中心差分近似函数曲线的斜率上的中心差分近似函数曲线的斜率 11 2 nn n yyy xh 二阶导数类似二阶导数类似 11 11 22 , nnnn nn yyyyyy xhxh 即一阶导数为即一阶导数为 第63页/共72页