1、1计算机仿真技术济南大学控制科学与工程学院授课教师：李实cse_1教1007室23.2 刚性系统的特点及算法宇航飞行器自动控制系统一般包含两个相互作用但效应速度相差十分悬殊的子系统，一个是控制飞行器质心运动的系统，当飞行器速度较大时，质心运动惯性较大，因而相对来说变化缓慢；另一个是控制飞行器运动姿态的系统，由于惯性小，相对来说变化很快，因而整个系统就是一个刚性系统。导弹、鱼雷等航行器的惯性导航也是如此，质点运动速度较慢，偏航与俯仰运动较快，导致两个变化速度相差十分大的系统.这些现象，主要是由于系统模型中的一些小参数，如小时间常数、小质量等存在而引起的.3刚性系统在用微分方程描述的一个变化过程中
2、，若往往又包含着多个相互作用但变化速度相差十分悬殊的子过程，这样一类过程就认为具有“刚性”。描述这类过程的微分方程初值问题称为“刚性问题”，该微分方程成为“刚性方程”，该系统成为“刚性系统”.4刚性系统例3.3以一个二阶微分方程组为例：212119999991998998yyyy满足初始条件：Ty010 经计算，该方程的解析解为：tttteetyeety1000110001)(2)(方程矩阵A的特征值为：1000,121可以看到，方程的解由两个特征值决定，因为特征值2=-1000非常小，对应系统的最小时间常数，当时间t增大，该项在方程的解中很快减小，并可以忽略.5刚性系统该表给出例3.3微分方
3、程的两个解的具体数值，其中y*代表将2忽略，只保留1的方程解的值，进行比较：t0120-10.0011.631.998-0.63112-0.990.011.980051.98010-0.99000-0.990050.11.809681.80968-0.90484-0.90484tteety100012)(tety 2)(*1tteety10002)(tety)(*2在0.01s时，精确解与简化解已经几乎一致，即在0.01s以后，可以将2忽略.而在t=0到0.01之间时，可以看到精确解与简化解的差别是很大的，这时不能将2忽略.这就是系统的边界层效应，即奇异摄动系统边界层效应，即奇异摄动系统.6刚
4、性系统例3.4代数方程的刚性系统，方程：7.91.46.67.98.21.421xx经计算，该方程解为：0,121xx如果将b摄动为：Tbb7.911.4这两个特征值的差别很大，该方程为刚性方程，说明了b摄动后引起了解的较大的变动.这时，方程的解变为：97.0,34.021xx如果我们计算矩阵A的特征值，可以得到：709.10,00934.0217刚性系统对于刚性系统来说，仿真算法的步长要由最大的特征仿真算法的步长要由最大的特征值来决定值来决定，也就是步长取很小的值，但会使求解时间加长，而随着时间增大，实际上可以忽略掉最大的特征值作用，不必要用很小的步长浪费计算机资源求解.一般使用Gear方法
5、方法来求解刚性系统，Gear首先进入刚性稳定性的概念，可以满足稳定性，同时又降低对步长的要求.对于刚性系统的认识，我们会在后面MATLAB建模与仿真中举例说明.83.3 实时半实物仿真实时数字仿真实时数字仿真通常是指把一个数字仿真过程嵌入到一个具有实物模型的实际系统或仿真系统的运行过程中.这种系统必须按照实际系统运行的时序要求来完成数字仿真过程.半实物仿真半实物仿真是指在仿真实验系统的仿真回路中接入部分实物的实时仿真.准确定义为“Hardware In the Loop Simulation(HILS)”，即回路中含有实物的仿真.实时性是半实物仿真的必要前提.实时半实物仿真实时半实物仿真具有更
6、高真实度，是仿真技术中置信度最高的一种仿真方法.因为允许在系统中接入部分实物，意味着把部分实物放在系统中考察，可以提高系统设计的可靠性和研制质量.随着科学研究和大型工程设计的需要，计算机技术的迅猛发展，实时数字仿真已经在航空、航天、核工业、电子、电力工业等领域中得到广泛的应用.9实时半实物仿真假设一个实际的动力学系统由实物系统过程A和B组成，在实时仿真过程中，可以将系统的一部分，比如A用数用数学模型来代替学模型来代替，在实际实现时，用计算机的数字处理过程来代替，如图所示.实物系统过程B实物系统过程A实物系统过程B计算机处理过程A采样A/Dz保持D/Ayzmym10计算机数字处理过程A必须在实物
7、系统同步的条件下获得动态输入信号zm，并实时的产生动态输出响应ym.zm和ym都是具有固定采样周期的数值序列.ym对输入zm的延迟称为响应时间.实时仿真要求响应时间满足系统的时间限制.实时半实物仿真实物系统过程B计算机处理过程A采样A/D保持D/Azyzmym11实时仿真算法由于动力学系统的实时仿真有实物系统介入仿真模型，所以要求仿真模型的时间比例尺完全等于原始模型的时间比例尺仿真模型的时间比例尺完全等于原始模型的时间比例尺.进行实时仿真时必须采用相应的实时仿真算法.要求实时仿真算法具有快速性、可行性、鲁棒性和相容性的特点实时仿真算法的快速性快速性实物系统过程B计算机处理过程A采样A/D保持D
8、/Azyzmym假设采样周期为T，zm=z(mT)，ym=y(mT)，由zm和ym计算出T+1时刻的ym+1值，因此计算时间需要小于采计算时间需要小于采样时间样时间T，否则实物系统得不到及时的数据.对于快速性的要求，希望采用单遍算单遍算法法(只需计算一次右函数值，比如显式方法)和大步长大步长积分算法.但对于右函数复杂，刚性问题或高频振荡的，则需要小步长，可以采用变化速率方法或构造并行算法.12实时仿真算法实时仿真算法的可行性可行性算法所用的输入信息应该是已经从实物系统或其它过程获取的，在计算时不允许使用还没有获得的信息.即kT时刻，只知道x(kT),x(k-1)T),但x(k+1)T)未知不能
9、使用.例如，二阶龙格-库塔法：1121211,5.05.0nnnnnnuKxhfKuxhfKKKxx(3.17)取步长h=T，即步长等于采样时间.在计算x(n+1)T)时，需要计算两次函数f(t,x)的值，在计算k2时需要知道u在n+1时刻的值，该值并未获得，所以该算法不能使用.13实时仿真算法改变后的二阶龙格-库塔法：2112121,21,nnnnnnuKxhfKuxhfKKxx(3.17)取u的采样周期为T=0.5h，将两次函数求值分为两步，每步在0.5h时间内求值.在时间区间nT nT+T上，计算函数f求得K1在时间区间nT+T nT+2T上求K2，这时已知u(nT+T)的采样值.这样，
10、仿真模型的周期为h=2T.这类算法称为实时二阶龙格实时二阶龙格-库塔算库塔算法法.14实时仿真算法实时仿真算法的鲁棒性鲁棒性(强壮性强壮性)算法的鲁棒性指在不同的复杂计算环境下，都能给出合理的计算结果，即适应不同环境和干扰的能力叫做鲁棒性，适应性越强，鲁棒性越好，算法和程序对不同的环境和异常干扰均具有良好的运行能力.由于这种要求，实时仿真算法应该具有处理异常因素的能力，必要时能够对计算流程进行重组、切换或使算法具有容错能力，可靠性高.如果算法具有迭代过程，应该能够保证在规定次数内结束该迭代过程.绝对不允许迭代不收敛和计算时间大于规定的时间.因此，为了避免超时现象发生，要有措施对数字处理过程测试
11、其计算过程的随机尖峰负载时间.15实时仿真算法实时仿真算法的相容性相容性相容性指当一个系统中的某个子系统由数字处理过程替代时，这个数字处理过程所用的数字仿真算法能保证在替代后的系统具有与原系统相同的动态特性.如用数字系统A与实物系统B构成的替代系统，与原系统有几个差别：实物系统过程B计算机处理过程A采样A/D保持D/Azyzmym实物过程A由数字过程替代，会引起模型误差模型误差和离散化误差和离散化误差.即数序模型不匹配或参数失真的误差，和离散化数值方法引起的误差实物A与B之间的信息传输由数字系统A与实物系统B之间的信息传输取代，采样、保持会引起量化误差量化误差原连续系统A与B是同步的，但数字系
12、统A与实物系统B是异步异步的，实物系统B给出的输出值需经过一个步长h时间后才能得到数字系统A给出的响应输出值，这里有一个延迟延迟时间h.16实时仿真算法：亚当斯-贝希霍斯显式算法(Adams-Bashforth)(AB)AB 1 阶法=向前欧拉法：nnnhfxx1AB 2 阶法：11321nnnnffhxx21151623121nnnnnfffhxxAB 3 阶法：AB 4 阶法：3211937595524nnnnnnffffhxx这些算法都是单遍算法，在采样时间间隔T=h时，在一个积分步长中只需要采样输入信息u一次，并对右函数f求值一次，计算工作量较小.17实时仿真算法：亚当斯-莫尔顿隐式算
13、法(Adams-Moulton)(AM)AM 1 阶法=向后欧拉法：11nnnhfxxAM 2 阶法：nnnnffhxx112111185121nnnnnfffhxxAM 3 阶法：AM 4 阶法：211151924nnnnnnffffhxx隐式的稳定区域要比显式算法大得多，但是隐式方法每积分一个步长需要求解一个非线性方程组，在复杂的函数情况下无法求得，在实时仿真中一般采用预测-校正法：将显式作为预测值，将显式法得到的预测值代入隐式法，最终得到校正值.18实时龙格-库塔法实时二阶龙格-库塔法：2112121,21,nnnnnnuKxhfKuxhfKKxx对于三阶算法，取u的采样周期为T=h/3
14、，分为三次函数求值：在时间区间nT nT+T上，计算函数f求得K1，采样获得u(nT+T)时刻值，在时间区间nT+T nT+2T上求K2，采样获得u(nT+2T)时刻值，在时间区间nT+2T nT+3T上求K3.实时三阶龙格-库塔法：322331121311,32,31,341nnnnnnnnuKxhfKuKxhfKuxhfKKKxx193.4 采样控制系统的仿真方法数字控制系统，或计算机控制系统，其控制器是由数字计算机组成.它的输入变量和控制变量只是在采样点（时刻）取值的间断的脉冲序列信号，描述控制器的数学模型是离散的控制器的数学模型是离散的差分方程差分方程或者离散状态方程或者离散状态方程.
15、而被控对象是时间连续的，其数学模型是连被控对象是时间连续的，其数学模型是连续时间模型续时间模型，它们之间用采样器和保持器相连接，整个系统可看作一个连续-离散混合系统.主要包括连续的控制对象、离散的控制器、采样器（A/D转换器）、保持器（D/A转换器）等几个环节组成，如图所示的典型采样系统.控制对象保持器（D/A）数字控制器r(t)e(t)e(kT)u(kT)u(t)y(t)20kT时刻的误差信号e(kT),经计算机计算后得到u(kT)，对于零阶保持器来说，u(kT)保持T时间，直到下一轮新的计算机输出结果u(kT+T)到达才改变一次数值.采样器，数字控制器，保持器并不是同步并行的工作方式，而是
16、串行流水的工作方式：但通常认为采样、控制器、保持器三者的计算处理时间是瞬时完成，即三者总和时间可以忽略不计，若将三者完成任务的时间加入，等于系统中增加一个纯滞后环节.控制对象保持器（D/A）数字控制器r(t)e(t)e(kT)u(kT)u(t)y(t)21对连续系统离散化，等于对被控对象加了一个虚拟的采样器，一般假设虚拟采样时间和仿真步长是一致的.而该步长与实际的采样周期可能相同，也可能不同，如何来确定连续系统离散化的仿真步长？也就是说，如何解决实际的离散部分仿真与虚拟的连续部分仿真的接口问题.可以考虑三种情况实际采样周期T与连续系统仿真步长h相等采样周期T大于仿真步长h采样周期T小于仿真步长
17、h控制对象保持器（D/A）数字控制器r(t)e(t)e(kT)u(kT)u(t)y(t)22采样周期T等于仿真步长h两者相等，则等于实际的采样开关与虚拟的采样开关是同步工作的，因此，这种仿真与连续系统仿真完全相同，是最为理想的情况，可以简化仿真模型，缩短仿真程序，提高仿真速度.当采样周期T较小时，可以取仿真步长h=T，满足仿真精度和稳定性.采样周期T大于仿真步长h这是实际系统仿真中最为常见的情况.采样周期T由系统频带宽度采样开关硬件性能和数字控制器执行时间所决定，一般大于仿真步长h.一般取T=kh，即取采样时间为仿真步长的k整数倍.离散部分按采样周期T进行仿真，将输入按保持器要求保持，连续部分
18、仿真模型按步长h计算k次，将最后第k次结果输出作为离散部分下一采样周期的输入.采样周期T小于仿真步长h有时为了减小计算量，加快仿真速度，采用大的仿真步长h来仿真，会出现h大于T的情况.这时需要对仿真模型作出修改.控制对象保持器（D/A）数字控制器r(s)e(s)e(z)u(z)u(s)y(s)23采样系统仿真方法将连续系统部分进行离散化处理，得到离散差分方程或脉冲传递函数.数字部分本身已经给出了离散差分方程或脉冲传递函数.就可以将两个系统相连，求出整体闭环系统的输出.例3.5 系统如图所示，采用计算机控制，控制器的脉冲传递函数和被控对象的传递函数已经给出，采用零阶保持器，T为采样周期.系统的初
19、始状态为：e(0)=0,u(0)=0,y(0)=0.求输出量y的表达式？TT保持器（D/A）r(t)e(t)e(n)u(n)u(t)y(t)11111zzkseTs1100sTk zGu(n)y(n)24y(n)根据离散化方法求连续传递函数的脉冲传递函数：0000/00/0000000001111)/1(11)1(111)(TTTTTTTTTsezeTkezzezzzTkTssZzzTksTskZzzsTkseZzG控制器的脉冲传递函数为：1111)1(zzkzCTT保持器（D/A）r(t)e(t)e(n)u(n)u(t)y(t)11111zzkseTs1100sTk令 ，得到脉冲传递函数为：
20、0/TTe100Tk111)(zzzzG25y(n)得到开环的脉冲传递函数为：y对设定值r的闭环脉冲传递函数为：TT保持器（D/A）r(t)e(t)e(n)u(n)u(t)y(t)11111zzkseTs1100sTk2121111111)1(11)1(1)()()(zzzzkzzkzzzCzGzGBO2111211)()1(1)(1)()(zkzkzzkzGzGzGBOBOBF得到系统的差分方程为：)2()1(2)(1)1(111nrnrknyknykny26y(n)如果连续系统的仿真步长h小于采样时刻T，取T=Nh，即要将连续系统仿真N个步长才能将输入y用于数字控制器.那么零阶保持器在Nh
21、内输入u(n)不变，此时连续系统的差分方程变为：TT保持器（D/A）r(t)e(t)e(n)u(n)u(t)y(t)11111zzkseTs1100sTk nuTknynuTkhNhnyNhnynuTkhnyhnynuTknyhnyNN111210000000027换成采样时间，T=Nh，得到：nuTknynyNN1100与T=h不同的是，这里仅把替换成为N，N=T/h，即可以得到总的差分方程：变换为脉冲传递函数，T=Nh，得到：hThTzTkzG/001)1()()2()1(2)(1)1(1/11/nrnrknyknyknyhThT0/TTeNTk100283.5 分布参数系统的数字仿真以上
22、章节介绍的都是常微分方程(ODE:ordinary differential equation)的数字仿真及模型，它们又叫集中参数模型.实际上，在工程应用中，更为广泛的是偏微分方程(PDE:partial differential equation)，又叫分布参数模型，即变量除了随时间轴变化，也随空间轴变化.比如热传递、流体、波的传递、电磁PDE的形式更为复杂，比ODE问题复杂得多.实际上，ODE往往是PDE问题的一种简化形式.比如，我们可以将3维空间(x,y,z三个空间坐标轴)里面的热传递简化为二维空间(x,y两个方向)热传递，进而可以简化为一维空间，即x方向，最后将空间维度忽略，即空间内各
23、处温度相等，只考虑时间变化，即变为ODE问题.T(t)293.5 分布参数系统的数字仿真利用计算机求解PDE问题时，可以将空间内的物体分解成无限小的小块，对每个小块进行ODE问题求解，实际上是无限个ODE的集合，并具有一定的边界约束.那么相应的数值运算工作量变大.比如，一维空间的热传导问题：因为存在着温度的不同，使得热量从温度高的地方传递到温度低的地方，叫热的传导(conduction).微观上温度可以看作分子活跃状态，即热能从温度高的地方传递到温度低的地方.取 u为x轴上某点在时间t时刻的热量(J)，可以写出热传导的方程式：22xubtux=0 x=lu(t,x)Fourier law:th
24、e time rate of heat transfer through a material is proportional to the negative gradient in the temperature and to the area.3022xubtux=0 x=lu(t,x)Fourier law:the time rate of heat transfer through a material is proportional to the negative gradient in the temperature and to the area.还要满足一定的边界条件：在x=0
25、和x=l时，热量u(t)对时间t的函数：)()(210tuutuulxx xut0在稳态情况下，或静态情况下，热量u(x)对坐标轴x的函数：Ttlx0,0最后，考虑t和x的边界约束条件：31PDE问题的数值求解方法有限差分法(Finite difference method)基本思想是把连续的区域用有限个离散点构成的网格来代替，这些离散点称作网格的节点；把连续区域上的连续变量的函数用在网格上定义的离散变量函数来近似；把原方程和定解条件中的微商用差商来近似，积分用积分和来近似，于是原微分方程和定解条件就近似地代之以代数方程组，即有限差分方程组，解此方程组就可以得到原问题在离散点上的近似解。txx
26、nxn+1Xn-1tntn+132PDE问题的数值求解方法有限体积法(Finite volume method)将计算区域划分为一系列不重复的控制体积，并使每个网格点周围有一个控制体积；将待解的微分方程对每一个控制体积积分，便得出一组离散方程。有限元法(Finite element method)将连续的求解域离散为一组单元的组合体，用在每个单元内假设的近似函数来分片的表示求解域上待求的未知场函数，近似函数通常由未知场函数及其导数在单元各节点的数值插值函数来表达。从而使一个连续的无限自由度问题变成离散的有限自由度问题。有限元法在工程应用中是最为常用的一种高效能算法，广泛的应用在各种物理场(力场，磁场，流体，传热，波)的描述中.33二维有限元法分析静态磁场34有限元法分析汽车撞击后的变形
