离散相似法

第五章 连续系统离散相似法数字仿真
上一章讲的数值积分法，其结果得到的是递推公式。

其出发点不是建立一个与被仿真的连续系统等价（某种误差意义下）的离散模型，而是建立在数值计算的基础上，即对积分进行数值计算（某种求和计算）。

或从另一角度说，对微分算子s（以一阶微分方程为对象） 找更好的替代方法。

所得的公式，只是这种思路的结果。

也就是说，人们为了计算方便，将其设法整理成了递推公式，递推公式在固定步长的情况下，就变成了通常所用的常系数的差分方程，可用z变换与频域联系起来（也正因为这点，使其结果从形式上与离散相似累积法类似）。

对这种方法精度的分析，是从截断、舍入误差累积角度进行的。

其改进思路，是减小这两种误差。

时域看失真。

可为线性系统，也可为非线性系统。

离散相似法的出发点与数值积分法不同。

它的出发点在于：
采样系统可用差分方程 和z变换来描述，而用计算机求解差分方程式很方便的（差分方程可看作多步法递推公式）。

所以，设法找到一个采样系统，使其输入—输出值与被仿真的连续系统在采样时刻的值近似相等。

这个采样系统称为连续系统的离散等价模型（二次模型、仿真模型）。

从此模型，可得出原系统的一系列值。

“等价”程度的高低，就是此法的精度问题。

所以精度分析从保持器引入的幅度与相位失真进行。

频域看失真。

只有线性时，才可列出整个系统的差分方程。

上述二法从结果（即所得递推公式）看有某种相似性。

所以精度稳定性等的解释可互换借鉴。

改进思路的不同，可看作看问题的角度不同（或出发点不同）而已。

因角度等价手段不同，各自能达到的“等价”程度也就有差别。

小结：
1.数值积分法：
出发点： 对积分进行数值计算，对微分算子s找更好的替代方法。

结果： 递推公式。

定步长时，为差分方程。

精度分析： 从截断、舍入误差累积角度进行。

时域看失真。

2.离散相似法：
出发点：出发点与数值积分法不同。

用计算机解差分很方便，找一个采样系统，使其I-O值与被仿真连续系统（在采样时刻）近似相等。

此采样系统称为连续系统的离散等价模型，用此采样系统的差分方程作为仿真模型。

结果：差分方程。

精度分析：从采样——保持器引入的失真进行。

频域看失真。

3.上述二法从结果（即递推公式）上看相似。

改进思路不同，可看作看问题的角度不同（或出发点不同）。

因为角度手段不同，各自能达到的“等价”程度也就有差别。

§1. 连续系统的离散化模型
可由下述5步获得（离散化过程）： ()()u t y t →→连续系统个G(s)
1. 对输入信号()u t 采样（加采样开关）：()u t ⎯⎯
→⎯采样
离散信号()u KT 。

2. 模型等价→信号等价：加上一个信号重构器（保持器），使()u KT ⎯→⎯()h
u t 。

3. 在系统的输出端也加一同步采样开关，()()y t y KT →。

4. 对()u KT 及()y KT 分别取z 变换，可得()U z 、()Y z ，从而得到：
[]()
()()()()
h Y z G z Z G s G s U z =
=， ()G z 为与原系统等价（在采样时刻）的离散模型。

5. 对()G z 取z 反变换，得到差分方程⎯→⎯
可放到计算机上求解：此差分为仿真模型。

为使()G z 能较准确地代表()G s ，要注意：
1.选择合适的采样周期T：
()u t ⎯⎯→⎯采样
离散信号()u KT ，为使()u KT 保持()u t 的全部信息，从抽样定理知，要求：
12m f T ≥，其中m f 是连续信号的最高频率。

则1
2m
T f ≤。

2.选择合适的信号重构器（保持器）：为从()u KT 完全恢复()u t ，即()()h u t u t =，要求: 理想的信号重构器（理想低通滤波器）+理想的低通信号()u t 。

这是不可实现的。

相频特性为水平O 线。

工程上可实现的重构器的特征，要尽量（取决于目的、精度）接近于理想低通滤波器。

也就是说，要选择合适的信号重构器。

小结： 当精度要求高时，要求()h G s 的幅、相频率特性尽量接近理想低通滤波器的，同时，信号本身频谱有限。

§2 信号重构器的特性及传递函数
工程上常用的重构器有：零阶、一阶、三角形保持器。

为有助于选择合适的重构器，下面讨论它的特性。

一. 零阶保持器：
对其他更为一般的信号，以上二保持器都将引起较大失真。

三角形保持器是一种较理想的保持器（∵h（t）为三角形，∴称为三角形保持器） 。

但有时不能实现，于是产生了滞后一拍的三角形保持器：
四．G（z）的求法：
∵Z[G（s）]常常可查表得到 ∴对上表进行简化（近似），得出二次近似的仿真模型。

例： 采用零阶保持器，求G(s)/s的z变换二次简化法
二次近似的应用还有：
例：
a. 一次近似离散化
b. 二次近似离散化
b精度＜a精度，但b的G(z)计算方便一些且可仿真非线性系统。

五．常用环节的离散化举例
2.惯性环节：
G s。

如果被仿真系统为齐次方程（即输入
离散相似法的误差仅仅来源于采样开关及()
h
=0），则不产生任何误差。

如上表，令u=0，则y K+1=e-aT y K，是精确解。

而Ealer，梯形法当u=0时，仍不为精确解。

∴T可取得大一些。

讨论：
①当原系统稳定时，由一次近似离散化法得到的差分方程本身是稳定的（由上表，令u=0，可知；或从差分方程的极点）。

∴一次离散相似法不论T取多大，由其得到的差分方程总是稳定的（但精度未必高）。

∴T可取的很大
G s会带来误差，这些误差表现为幅度与相位失真。

T↑则失真↑；当T很
②∵开关+()
h
大时，产生的延迟过大→整个系统不稳定（闭环时）。

∴T仍需限制。

①相当于从一个环节着眼看稳定性→T可无限大
②相当于从整个系统（有很多个采样保持器及反馈）着眼看稳定性→T必须小于
某值。

③比较数值积分法与离散相似法：
数值积分法：仿真模型从数值计算法得到，有较明确的几何意义，存在计算的不稳定性问题（即原系统稳定，而h大于一定值后，差方不稳定）。

离散相似法：仿真模型从离散等价模型角度得到，有较明确的物理意义，不存在数值计算的不稳定性问题（即原系统稳定，则一次离散化的差方稳定）。

由于上二特点，许多快速数字仿真法从离散法推出。

§3 线性连续系统状态方程的离散化（时域离散相似法，状态转换法，状态相似法）
一． 离散化：
对①两边取拉氏变换，得：
sX(s)—X(0) = AX(s)+BU(s)
或 （sI-A）X(s)=X(0)+BU（s）。

以（sI-A）-1左乘上式两边，得：
X(s)=（sI-A）-1X(0)+（sI-A）-1BU(s) ② 令 L-1[(SI-A)-1]= Φ(t), Φ(t)为状态转移矩阵。

则：
X(s)=L[Φ(t)]X(0)+L[Φ(t)]BU(s)
对上式两边进行反变换，得：
X(t)= Φ(t)X(0)+∫0tφ(t-τ)Bu(τ)dτ ③
可以证明，有：
Φ(t)=e At， e At=I+At+(At)/2！+···
所以，③亦可写成：
X(t)=e At X(0)+∫0t e A(t- τ)Bu(τ)dτ ④ 上式为连续系统状态方程的解。

令t=KT、(K+1)T,有：
X(kT)=e AKT X(0)+∫0KT e A(KT-τ)Bu(τ)di ⑤
X[(K+1)T]=e A(K+1)T X(0)+∫0(K+1)T e A[(k+1)T-τ]Bu(τ)dτ ⑥
⑥—e AT·⑤：
X[(K+1)T]= e AT X（KT）+∫KT(K+1)T e A[(k+1)T-τ]Bu(τ)dτ ⑦ ⑦为精确解.
1) 若U（t）在k→k+1间保持不变(相当于加了一个零阶保持器)，则：
X[(K+1)T]= e AT X（KT）+[∫0T e A(T-τ)Bdτ]u(KT) ⑧
讨论：当u（t）为阶跃信号时，⑧式为精确解（要求U(0)=U(0+)）。

即⑧的离散化对阶跃信号不引入误差。

（广义函数在积分意义下存在）。

已知e AT =φ（T），定义∫0T e A(T-τ)Bdτ=∫0Tφ(T-τ)Bdτ=φm(T)，则有：
X[(k+1)T]= φ(T)X(kT)+ φm(T)u(kT) ⑨ 上式即为连续系统离散化后的解。

该解是零阶保持器意义下的，对阶跃输入是精确解。

若已知A、B，则φ(T)，φm(T)阵可事先求出来。

这样，利用差分方程⑨可十分容易地求出不同采样时刻的状态变量值。

若u（t）为阶跃，则该解还是精确解。

2) 假定u(t)在两相邻采样点间斜线变化（相当于加了一个三角形保持器）：
则⑦为：
X(k+1)=Φ(T)X(k)+Φm(T)u(k)+
ˆ()
m T
φ()
u k•
(11)
上法较精确，用其构成的仿真模型，误差来自于： c输入采样+保持
dΦ、Φm、
ˆ()
m T
φ
计算。

对(11)而言：若原系统稳定，则（11）必稳定∴恒稳计算式。

此法与数值积分法的差别在于：
数值积分法： X(k+1)=X(k)+
(1)
[] k T
kT
AX BU dt +
+
∫
状态相似法： X(k+1)=AT
e X(k)+(1)()
()k T
A T kT
Bu d e
τττ
+−∫
上面两个公式均需要对积分进行近似计算，但前者对右函数进行近似积分，后者仅对输入进行近似积分。

所以，时域离散相似法是恒稳的（为何？），数值积分法得到的差分方程不是恒稳的。

若整个系统不用状态方程表示，而用典型环节表示。

对各典型环节分别离散化，则整个系统的差分方程本身可能不稳定。

原因为：从物理上说，采样-保持引入滞后；从数学上说，各环节输入信号u（t）为x(t)的组合，对u（t）近似积分相当于对x(t) 近似积分。

二. Φ(T)、Φm(T)、ˆ
()m T φ阵的求法： 1. 解析法：
优：精确解，所得差分方程恒稳。

缺：麻烦，不适合计算机求解。

2.级数展开法：
优点：计算方便，缺点：所得差分方程不一定恒稳。

§4.常用典型环节的离散化 肖P.66
§5.连续系统（包括非线性原件）的离散化数字仿真程序 肖P.68
第六章 连续系统快速数字仿真方法
思路：①加大步距，减少步数。

------>从恒稳算式处得到 ②加快每步计算速度。

一般，速度的提高以精度的降低为代价。

§1.时域矩阵法
数字信号处理中：输出＝输入*脉冲响应--->
∑
或矩阵乘法。

时域中用无穷矩阵来分析和设计的方法，特别适合于采样系统。

但通过引入
也可以对连续系统进行研究。

好处：若能获得系统的脉冲响应g(t)，则采用时域矩阵法对各种不同输入信号下的系统仿真时，计算量小，且与系统的阶次无关。

缺点：g(t)的要花费时间。

一、时域矩阵：初始条件为
0。

设 g(t)⎯→
←L
G(s)，则有： C(K)= r(n)
n)-(K 0∑=K
n g （1）
其中：g(K)=KT
t =g(t)。

将（1）展开，得：
⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡)(┋)2()1()0(k C C C C = ⎥⎥⎥
⎥
⎥⎥⎦
⎤⎢
⎢⎢⎢⎢⎢⎣⎡−−)0(...)
2()1(g(k)┋...┋┋┋0...)0()1()2(0...0)
0()1(0 0
0)0(g k g k g g g g g g g ⎥⎥⎥⎥⎥⎥⎦⎤
⎢⎢⎢⎢⎢⎢⎣⎡r(k)┋)2()1()0(r r r （2）
C = G * R
若响应在所有采样点上均要求出，则G 阵阶为无穷大，称G 为无穷矩阵。

在此称C、G、R 为时域矩阵，以与状态空间法中的状态阵A、B 等区别。

从（2）知，若已知G、R，则可求出输出C。

关键是G，G 阵称之为时域转移矩阵。

二、时域转移矩阵G 的求取：
g(t)=L 1−[])(s G
g(k)=g(t)
kT
t = 或g(k)⎯→
←Z
G(Z)=L 1−[])(s G
三、求闭环系统的响应：
由（2）知： ⎩⎨
⎧==C -R E E
*G C 。

所以，有：
⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡)(┋)2()1()0(k C C C C = ⎥⎥⎥
⎥
⎥⎥⎦⎤⎢
⎢⎢⎢⎢⎢⎣⎡−−)0(...)
2()1(g(k)┋...┋┋┋0...)0()1()2(0...0)0()1(0 0
0)0(g k g k g g g g g g g e(k)
┋)2()1()
0(c(k)-r(k)┋)2()2()1()1()0()0(e e e c r c r c r ⎥⎥⎥
⎥⎥
⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡−−− 又 []
⎪⎩
⎪
⎨⎧+=∴−=)0(1)
0()0()0()0()0()0()0(g r g c c r g c ∵ [][][]{})1()0()0()0()1()
0(11
)1()1()1()0()0()0()1()1(r g c r g g c c r g c r g c +−+=
⇒−+−=∴ 如此下去可求出c(1)、c(2)....c(k)，较麻烦。

对一般的物理系统，输出不能瞬时建立起来，总有一小段延迟。

故g(0)=0 → e(0)=r(0),c(0)=0。

例1：
1
1
+s
：
此时，有：
(0)(1)(2)()C C C C k ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦ = 000
...0(1)00...0(2)(1)0...0...g(k)(1)(2)...
0g g g g k g k ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢
⎥−−⎣⎦ (0)(1)(1)(2)(2)r(k)-c(k)r r c r c ⎡⎤⎢⎥−⎢⎥⎢⎥−⎢⎥
⎢⎥⎢⎥⎣⎦
对此阵，当求c(i)时，已知c(i-1)、c(i-2)...c(1)、c(0)，所以：
[][]()()(0)(1)(1)(1)...(1)(1)(1)c i g i r g i r c g r i c i =+−−++−−−
这可容易地算出来。

∴仿真过程是一个行向量与一个列向量相乘的问题，计算很快。

讨论：
①G(s)较复杂时，求g(t)费力；
②K 越大，求C(k)的计算量越大；
③已知g(t)，对不同r(t)求c(t)，很方便； ④与阶次无关，仿真高阶的快；
⑤优点：计算量小，精度较高；缺点：高阶系统g(t)难求。

精度： 取决于采样+保持。

时域矩阵法程序图（N：计算总步数+1； T：采样间隔）
作业：
其中：G 0(s)=)1(10
+s s ，
D(z)=)1)(1()1)(1(1
111−−−−+−−−cz z bz az k 。

若直接以T 为采样周期，应用时域矩阵法，则只能求得C(KT)，K=0，1，2，...N。

现
还想知道C(t)在采样点之间的几个值，则可对图1做如下变换：
则可求出C(KT 1)，T 1=
l
T
，1l ≥且为正整数。

G(s)=)1(1011+−−s s s e S T =10(1-e S
T 1−)(21s -s 1+1
1+s )→g(t)→g(KT 1)
要求：T=1，T 1=
l
T ，计算C(KT 1)。

（1）时域矩阵法程序框图。

（2）对T,N=6；T1,N=1l ≥。

编出Matlab 程序。

§2 增广矩阵法
一、基本思想
设连续系统状态方程为： X
=AX+Bu （1） 其中u 为系统输入信号。

则：
X(t)=e
At
X(0)+
τττ)d (0)
(Bu e
t
t A −∫ t≥0 （2）
当u=0时，有： X(t)=e
At
X(0)。

已知 e
At
=I+At+()!
22
At +…+
N At N )(!
1
+…。

若取前五项，则精度与四阶龙格-库塔法相同。

就是说，对X
=AX 有X(t)=e At
X(0)。

选定仿真步长为T，且只取e
AT
前五项，则有：
X ()1n T +⎡⎤⎣⎦=⎥⎦
⎤⎢⎣⎡+++!4AT I 4
4T A …X(nT) （3）
其截断误差为 0(T 5
)。

∵A、A 2、A 3、A 4可以事先算出。

∴（3）右端的系数项可事先求出→递推计算变式
简单的乘法。

现要求解的方程是（1），他的解（2）除一个自由项外，还有一个强制项:
dt t Bu e t
t A )(0)(τ−∫。

为了能利用（3）法，就要求将u(t)增广到状态变量中去，使（1）
成为齐项方程。

当u(t)是一些典型函数时，上述增广的思想很容易实现。

二、典型输入函数时的增广矩阵
见教材P150-152。

今假设u 为标量，即系统为单输入。

优点：计算量小，精度高；
缺点：对特定输入项须特定处理，只适用于典型输入x。

三、病态系统仿真
系统有的时常数很大，有的很小。

设τ
min =0.05，max τ=10。

则病态系数s=
05
.010
=200,为保证计算稳定及一定精度，计算步长T≤0.05。

若用等间隔计算，系统过渡时间约为50s，则总共需要计算1000步以上。

→时间↑，舍入误差↑
解决的一个方法如下： 跳跃式计算
)0()(X e T X AT =
)0()0()()2(2X e X e e T X e T X AT AT AT AT === )0()0()2()4(4222X e X e e T X e T X AT AT AT AT ===
)0()0()2(2
22
11
X e X e
e T X AT
AT
AT K K
K K ==−−
上面AT
K
e 2
不必重新算，只需让4
2）
（）、（AT AT AT …各乘以系数相加即可。

在t=T K 2后，采用等距计算，计算步距T T K
21=。

则 )2()()2(2
111
T X e T X e T X K AT
AT K
==
)2()3(111
T X e T X AT = )3()4(111
T X e T X AT =
这种方法特别适合于S 很大的系统。

开始时系统主要由min τ起作用，此时步长较小。

后来主要由max τ起作用，此时步长1T 很大。

中间，变步长。

仍以上例为例，选跳跃式。

T=0.03125＜0.05，226
1==T T 秒。

跳跃式计算，当1T =2秒后，改为等步计算（到50秒需要25步），则总共计算32步即可把整个过渡过程计算出来。

§3 替换法
上述二法共同优点：①计算量小 ②精度较高 →较适合快速数字仿真。

但二法在使用上却有其局限性。

时法求高阶系统的g(t)较难，增法只是用于一些典型输入函数的情况。

所以，要求进一步探索更为简单有效的方法。

快速数字仿真
⎯⎯⎯→⎯往往要求快速性↑，精度要求较低。

对高阶G(s)，如能方便地从它的传递函数G(s)直接推出与其相匹配的，T 可较大的G(z)，对快速数字仿真十分有利。

“相匹配”指： ①G(s)稳定，则G(z)亦稳定。

②对同一输入，由G(z)、G(s)各自求出的输出有相同特征，比如终值相等，等等。

从G(s)
⎯⎯→⎯直接G(z)的方法有两种： ① 替换法 ② 根匹配法
替换法：找到s 与z 的一个对应公式，然后将G(s)中的s 换成z →G(z)
根匹配法：找到一个G(z)，它与G(s)有相同的零、极点，系数由中值相等确定。

一、欧拉公式中找：
已知：s、z 的关系为： sT z e =
但用其导出的G(z)难以化成有理函数→难以得到差方。

设),()(t y f t y
= ，则由欧拉公式得： ),(1n n n n t y Tf y y +=+n n y
T y += ∴y
T y z =−)1( ∴11−==z T
y y s ∴T
z s 1−= 上式虽简单但不实用，∵在T 较大时这样得到的G(z)不稳定。

从s、z 平面的对应关系
来看，若G(z)的极点均在z 平面的单位圆内，则其对应的系统稳定。

G(s)的四个极点如上图所示。

当T 较大时，用T
z s 1
−=转到z 平面，将有两个极点在单位圆外。

∴要T ↓→非快速性。

∴G(s)稳定⎯→⎯↑
T G(z)不稳定。

∴这种替换不能用。

二、图士丁（Tutin）替换法：双线性替换
已知梯形积分替换公式 ：
)(2
)(2111
+++++=++=n n n n n n n y y T
y f f T y y ∴ y z T y z )1(2)1(+=
−
y z z T y 112+−•= ∵
sy y
= ∴112+−=
z z T s s T s T
z 2
121−+=
映射关系：
∴恒稳
梯形法为恒稳算式（隐式稳定性高的原因）。

但因f k+1不能准确计算→迭代数步↑
误
差↑
不恒稳。

但用在替换法中恒稳。

已知： sT
z e = ，所以 231211111ln ()()13151z z z s z T T z z z −−−⎡⎤
=
=+++⎢⎥+++⎣⎦
∴相当于取级数的第一相近似→图示丁近似
例. 见教材P153 或本讲P40
∴此法速度高，精度相当于梯形法（或二阶龙格—库法），但恒稳；使用场合：线性时不变；数模形式：传递函数。

§4．根匹配法
1． 零极点一一对应，K 由终值决定。

2． 零点的对应：s 有无穷远零点时，z 的零点定在：
① Z=0处
② Z=-1处（由Tustin 法得到：11
2+−=z z T s ）
③ (0，-1)之间，具体是什么位置，用转折频率处相位误差最小决定；同时，K 由转折频率处幅度误差最小决定。