时域离散相似法

第三章 时域离散相似法用数字计算机对一个连续系统进行仿真时，必须将这个系统看作一个时间离散系统。

也就是说，我们只能计算到各状态量在各计算步距点上的数值，它们是一些时间离散点的数值。

在第二章中主要是从数值积分法的角度来讨论数字仿真问题，没有显式地涉及到“离散”这个概念。

史密斯从控制和工程的概念出发提出离散相似问题[1]，并导出离散相似法。

所谓“离散相似法”就是将一个连续系统进行离散化处理，然后求得与它等价的离散模型。

由于连续系统的模型可以用传递函数来表示，也可以用状态空间模型来表示，因此，与连续系统等价的离散模型可以通过两个途径获得，其一是对传递函数作离散化处理得离散传递函数(或脉冲传递函数)，称为频域离散相似模型。

其二是基于状态方程离散化，得到时域离散相似模型。

本章介绍时域离散相似法，第四章介绍频域离散相似法。

3.1 时域离散相似法基本原理3.1.1基本方法假设有一个连续系统，它由以下状态方程描述：xAx B =+u (3.1) 对于(3.1)式描述的连续系统进行离散化处理，如图3.1所示。

在系统的输入端加上虚拟采样开关和虚拟信号重构器，输出端加一个虚拟采样开关。

虚拟采样周期为T 且同步。

其中，u (t )是系统输入；u (k )是加虚拟采样开关后，在kT 时刻系统输入；x (k )是加虚拟采样开关后在kT 时刻系统输出；~()ut 、~()x t 是等价的连续信号。

只要~()ut 能足够精确地表示u (t )，那么~()x t 也就能足够精确地表示x t ()，这样，就能获得与连续系统等价的时域离散相似模型。

对该连续系统进行离散化处理后可以得到系统离散相似模型如式3.2所示：x [(k +1)T )]=)(T Φ x (kT )+)(T m Φ u (kT )+ Φm (T) ()u kT (3.2) 其中：T 是采样时间间隔(或称采样周期)；u (k )、x (k )为系统kT 时刻的输入及状态量；)(ˆ)()(T T T mm ΦΦΦ、、为离散化后与系统模型有关的系数。

下面将讨论怎样从(3.1)式经离散化处理获得(3.2)式所示的时域离散相似模型。

将方程两边取拉氏变换，经整理得到：x (s )=(s I －A )-1x (0)+(S I －A )-1x (0)+(s I －A )-1B u (s ) (3.3) 设 Φ(t)=L -1[(s I －A )-1 ] (3.4) 则Φ(t)=e At 称作状态转移矩阵。

将(3.4)式代入(3.3)，可得到：x (s )=L [Φ(t )]x (0)+L [Φ(t)]B u (s ) (3.5)由卷积公式：L f g t d F s G s t[()()]()()τττ0⎰-=称作f 和g 的卷积，可以写作f *g ，其中，Lf (t) =F (s)，Lg (t) = G (s )．(3.5)式取反拉氏变换，运用卷积公式得到：x t ex e AtA t t()()()=+-⎰00τB ~()ud ττ (3.6)图3.1连续系统的离散化处理对离散化处理后的系统，设kT 及(k +1)T 为两个依次相连的采样瞬时，则有：x (kT ) = eA kTx (0)+e A kT kT()-⎰τ0B ~()ud ττ (3.7) x [ (k +1)T )] = eA (k+1)Tx (0)+e A k T k T[()]()+-+⎰101τB ~()ud ττ (3.8) 将(3.8)式－(3.7)式*e A T,可得：x [ (k +1)T )] = e A Tx (kT )+e A k T kTk T[()]()+-+⎰11τB ~()ud ττ (3.9) 由于(3.9)式右端的积分与k 无关，故可令k =0，若信号重构器使kT 与(k+1)T 之间的~()u t 不变，即积分式中的~()uτ保持常数~()()u u kT τ=,那么，(3.9)式可改写为： x [ (k +1)T )] = e A Tx (kT )+e A T T()-⎰τ0B d u kT τ()=Φ(T )x (kT )+T⎰Φ(T -τ)Bd u kT τ() (3.10)若令T⎰Φ(T －τ)B d τ=Φm (T )，则有：x [ (k +1)T )] =Φ(T )x (kT )+Φm (T )u kT () (3.11)这就是一个离散状态方程。

(3.11)式是在假定信号重构器使输入量u (t )在两个采样时刻之间保持不变这样的前提下推导出来的，如图3.2中的矩形近似。

实际上，输入量在两个采样时刻之间是变化的，这样就会引起误差。

为了减小误差，可以假定在两个采样时刻之间，信号重构器使~()u τ为一斜坡函数，即用图3.2中的梯形近似。

此时，在kT 与(k +1)T 之间，相对矩形近似有一个∆u k ()τ存在∆u u k T u kT Tk ()[()][]ττ=+-1 ≅u kT ()τ (3.12)对应∆u k ()τ，对 x [ (k +1)T )]引起的变化量为：∆x [(k +1)T )]=⎰0TeA (T-τ)B ∆u k ()τd τ≅T T0⎰eA (T-τ)B d τ ()ukT (3.13) 若令ττ⎰e A (T-τ)B d τ = Φm (T)， 则整个离散状态方程为： x [(k +1)T )]= Φ(T )x (kT )+ Φm (T )u (kT )+ Φm (T) ()u kT (3.14) 至此，我们得到式(3.2)中与系统模型有关的系数为：Φ(T )＝e A T(状态转移矩阵)u (t 图3.2矩形近似与梯形近似Φm (T )＝0T⎰Φ (T －τ)B dτ (输入信号采用零阶重构器引入的系数矩阵)Φm (T)＝0T ⎰τ e A (T-τ)B d τ (输入信号采用三角形信号重构器后叠加的系数矩阵) 下面举一个例子来说明如何用离散相似法来进行数字仿真。

假设有一系统，如图3.3所示。

由图可知，该系统的开环部分是由一个非线性环节(饱和环节)及一个线性环节Ks s ()+1所组成。

今将线性部分化成时域离散相似模型。

已知线性部分状态方程为：x=Ax +B uy=Cx(3.15)其中：A ＝0011-⎡⎣⎢⎤⎦⎥，B =K 0⎡⎣⎢⎤⎦⎥，C =[]01，，根据前述，Φ( T )=e A T =L -1[(s I －A )-1] 因为 s I －A =ss 011-+⎡⎣⎢⎤⎦⎥， 所以有： (s I －A )-1=101111s s s s ()++⎡⎣⎢⎢⎢⎤⎦⎥⎥⎥ Φ( T )= L -1[(s I －A )-1] =101-⎡⎣⎢⎤⎦⎥--ee TT Φm ( T )＝T⎰Φ (T －τ)B d τ＝10100-⎡⎣⎢⎤⎦⎥⎡⎣⎢⎤⎦⎥---⎰e e K d T T ()τττ =KK e d T T (()10-⎡⎣⎢⎤⎦⎥--⎰ττ=KT K T e T ()+-⎡⎣⎢⎤⎦⎥-1则由(3.14)式可得线性动态环节的离散状态方程x k x k e e x k x k KT K T e u k T T T 1212111011()()()()()()++⎡⎣⎢⎤⎦⎥=-⎡⎣⎢⎤⎦⎥⎡⎣⎢⎤⎦⎥++-⎡⎣⎢⎤⎦⎥--- (3.16)式中(k +1)即表示[(k +1)T ]，(k)即表示（kT ），以后为简单起见均如此表示。

利用(3.16)式即可进行数字仿真，仿真程序的框图如图3.4所示。

用此程序可研究放大倍数K 对系统输出响应之影响。

例如，设输入信号r (t )为阶跃函数，幅度为25，且c 1=20，图3.3 被仿真的系统结构图并选采样周期T＝0.05，则可得K＝1及K＝2表3.1仿真结果t y(K=1) y(K=2)0 0 01 7.95649 15.17392 20.3824 30.59993 27.4777 30.1134 28.6091 24.57835 26.9278 23.08746 25.2142 24.51337 24.4764 25.548 24.5136 25.31119 24.8072 24.904110 25.0235 24.862711 25.0875 24.992112 25.06 25.045613 25.0155 25.016414 24.9906 24.98915 24.9871 24.9909图3.4用离散相似法进行数字仿真与第二章介绍的数值积分方法比较，方程系数ΦΦΦ、、可以一次求出，()() ()T T Tm m这样每做一步积分只要计算一次右端函数。

而四阶龙格－库塔数值积分方法，每做一步积分要计算4次右端函数。

显然，用离散相似方法做仿真运算比用数值积分法作仿真运算运算量小，对于同一系统作仿真时，离散相似方法运算速度快。

由本小节所介绍的离散相似法原理可知，离散相似法的基本问题是计算、、，而它们的计算主要归结为关于e A T的计算。

下面我们将进一步讨论T T TΦΦΦ()() ()m m这一问题。

3.1.2 状态转移矩阵的计算离散相似方法中的一个问题是状态转移矩阵e A T的计算。

一些常用的一阶、二阶环节，它们的ΦΦΦ、、可以很方便地通过直接解析的方法来求解。

对于高阶系统，若T T T()() ()m m要求列出系统的离散状态方程，那么e A T必须利用数值计算才能得到。

1.泰勒级数展开法e A T的算法很多[2]，但是用得最多的还是泰勒级数展开方法，即将e A T展成泰勒级数形式，然后截取其中的前若干项。

这个方法首先是由（Lion）提出的。

状态转移矩阵e A T可以用无穷级数来表示，如(3.17)式所示。

e A T＝A T i i ii !=∞∑0， A 0＝I (3.17) 这个级数在有限的时间间隔内是均匀收敛的，因此可以在一定精度内计算e A T 。

若级数在i =L 处截断，则我们可以写成e A T＝A T i i ii L!=∑+0A T i i ii L !=+∞∑1＝M ＋R (3.18) 式(3.18)中的第一项为级数的近似值，而第二项为余数项。

要求： r E m ij ij ≤ E =10-d ，d 为正整数或 min max m E r ≤ (3.19) 其中r ij 和m ij 对应为R 与M 的元素， m ax r 为r ij 中最大元素。

式中m min 是容易求出的，但 r max 却无法求出，因为R 仍是一个无穷项的和。

下面来估计m ax r 。

令R 为矩阵R 的范数，根据矩阵范数的定义，R ∑==nj i ijr1,，有≤m ax r R ，由 ∑∑∞+=∞+=≤=11!!L i iiL i ii i T A i TA R))2)(3(21()!1(22!1+++++++=++L L T A L T A L T AL L))2(21()!1(22211++++++⋅≤++L T A L T A L T AL L令ε＝＋2L T A (3.20)则有 )1()!1(211++++⋅≤++εεL T AR L L如果1<ε，则 )11()!1(11ε-+⋅≤++L T AR L L因此，若E ≤-⋅+⋅++ε11111)!(L T A L L m min (3.21) 则满足（3.19）式，e A T 可以按照以下迭代过程来计算：(1) 选择初始L ； (2) 计算矩阵M 及m min ；(3) 用式(3.20)求ε；(4) 用式(3.21)判别是否满足精度的要求。