系统建模过程辨识脉冲响应

实验二：基于三阶系统的脉冲响应求传递函数
一、目的
所谓辨识，就是在输入和输出数据的基础上，从一组给定的模型类中，确定一个与所测系统等价的模型。

此次作业任务是基于脉冲响应（既包含了输入数据也能得到输出数据）对一个三阶模型进行过程辨识，整个辨识路线如图1所示。

图1基于脉冲响应的过程辨识路线
二、辨识方法
1）差分方程法
假设待辨识的过程传递函数为：
G(s)=b m s m+b m−1s m−1+⋯+b1s+b0
a n s+a n−1s+⋯+a1s+1
,n≥m（1）如果特征方程a n s n+a n−1s n−1+⋯+a1s+1=0具有n个单根s1、s2、…s n，那么传递函数可以写成：
G(s)=c1
s−s1+c2
s−s2
+⋯+c n
s−s n
（2）
对应的脉冲响应:
g t=c1e s1t+c2e s1t+⋯+c n e s n t（3）
如果特征方程a n s n+a n−1s n−1+⋯+a1s+1=0具有重根(其中s1、s2…s n−r为单根，s0为r阶重根)时，那么传递函数可以写成：
G(s)=c1
s−s1+c2
s−s2
+⋯+c n−r
s−s n−r
+c n−r+1
s−s0
+c n−r+2
(s−s0)2
+⋯+c n
(s−s0)r
（4）
对应的脉冲响应:
g t=c1e s1t+c2e s1t+⋯+c n−r e s n−r t+⋯+c n−r+1e s0t
+c n−r+2te s0t+⋯+c n t r−1e s0t（5）可以根据待定过程的脉冲响应估计值g(t)来确定c i和s i。

从所获得的g(t)按照采样周期T0，选取前(n+1)个坐标点，各坐标点上的脉冲响应分别记为g(k)、g(k+1)、…、g(k+n)。

组成一个AR模型：
g k +α1g k +1 +⋯+αn g k +n =0
（6）
其中，α1、α2、…、αn 为待定系数。

如果特征方程
1+α1x +α2x 2+⋯+αn x n =0
（7）
有一个单根x i ，则x i k 必是AR 模型的解，它们的线性组合也是AR 模型的解。

当特征方程式（7）具有n 个单根时，AR 模型的解可以表示为:
g k =β1x 1k +β2x 2k +⋯+βn x n k
（8）
当特征方程式（7）具有重根时，AR 模型的解可以表示为:
g k =β1x 1k +β2x 2k +⋯+βn−r x n−r k
+βn−r +1x 0k
+βn−r +2kx 0k +⋯+βn k r−1x 0k
（9）
其中，其中x 1、x 2…x n−r 为单根，x 0为r 阶重根。

将式子（3）和式子（5）分别与式子（8）和式子（9）进行比较，我们可以得到：
c i =βi s i =
1T 0
logx i
（10）
这样就可以求出待辨识过程的传递函数。

2） Hankel 矩阵法
考虑一个n 阶的脉冲传递函数
G z
−1
=
b 1z −1+b 2z −2+⋯+b n z −n 1+a 1z +a 2z +⋯+a n z （11）
将传递函数转化称状态方程后，经过进一步推导，可以得到
G z −1 =g 1 z −1+g 2 z −2+g 3 z −3+⋯
（12）
由式（11）和（12）可得：
b 1z −1+b 2z −2+⋯+b n z −n
=[g 1 z −1+g 2 z −2+g 3 z −3+⋯][1+a 1z −1+a 2z −2+⋯+a n z −n ]
=g 1 z −1+⋯+[a 2+a 1g 1 ]z −2 +⋯+[g n + g i n−1i =1a n−i ] z
−n
+⋯ （13） 比较式（13）同幂次项系数，有：
g 1 g 2 …g n g 2 g 3
…g n +1 ⋮g n ⋮g n +1 ⋮⋮
…g 2n −1
a n a n−1⋮a 1 = - g n +1
g n +2 ⋮g 2n
（14）
以及
b 1b
2⋮
b n
= 1
…00a 1
1…00⋮a n−1⋮
a n−2
⋮⋮⋮…
a 10 g 1 g 2
⋮g n
（15）
根据Hankel 矩阵的定义
H l,k=
g k g k+1⋯g k+l−1
g k+1g k+2⋯g k+l
⋮
g k+l−1
⋮
g k+l
⋮
⋯
⋮
g k+2l−2
（16）
式（14）左边的矩阵是一种待定的Hankel矩阵，记作H（n,1），并且它是可逆的。

因此只要将所获得的脉冲响应估计值g(k)填入式（14）和式（15）的相应位置，便可以求得式（11）的脉冲传递函数估计值。

三、实例
1）考虑传递函数为G(s)=1
s+6s+11s+6
的过程；
2）对上述过程基于Matlab进行仿真，得到脉冲响应曲线如图2所示.
3）设置采样周期为0.5s,采样脉冲响应数据，记录在表1。

图2给定过程的脉冲响应曲线
表1脉冲响应数据（采样周期T0=0.5s）
4）确定待定系数αi
由于待辨识的系统是三阶的，根据式（6），取k=0、1、2，可得如下方程组：
0+0.0469α1+0.0735α2+0.0672α3=0
0.0469+0.0735α1+0.0672α2+0.0507α3=0
0.0735+0.0672α1+0.0507α2+0.0346α3=0
（11）
解得：
α1=−16.5874
α2=49.8207
α3=−42.9166
（12）5）求特征方程的根
根据式（7），对应的特征方程为：
1-16.5874x+49.8207x2-42.9166x3=0 （13）解得
x1= 0.0769
x2.3= 0.5420± 0.0971j（14）求传递函数的极点
s1=−3.7455
s2.3=−1.1934± 0.3545j
根据（3）式（如果特征方程具有重根，则根据（4）式），取K=0,1,2，得如下方程组：
β1+β2+β2=0
0.0769β1+0.5420+ 0.0971jβ2+0.5420− 0.0971jβ3=0.0469
(0.0769)2β1+0.5420+ 0.0971j2β2+(0.5420− 0.0971j)2β3=0.0735
（15）解得
β1=0.1004
β2.3=−0.0502 ±0.4819j（16）6）令c i=βi
7）确定传递函数
根据式（2），得出待辨识过程的传递函数为G(s)=−0.5979s−1.573
s+6.132s+10.49s+5.805。

其实际传递函
数为G(s)=1
s3+6s2+11s+6
，可见，利用差分方程法可以获得较好的近似传递函数。

8）构造Hankel矩阵H(n,1)
H3,1=g(1)g(2)g(3)
g(2)g(3)g(4)
g(3)g(4)g(5)
=
0.04690.07350.0672
0.07350.06720.0507
0.06720.05070.0346
（17）
9）根据式（14）和式（15），分别求得：
a1
a2 a3=
−1.9298
1.2616
−0.2930
b1
b2
b3
=
0.0469
−0.0170
−0.0155
（18）
根据式（11），过程脉冲传递函数的估计值为
G z−1=（19）
10）经过双线性变换后，连续传递函数的估计值为
G(s)=0.02714s2+0.4456s+0.2105
（20）
s3+2.455s2+2.483s+0.5671
对比式（20）和实际传递函数，发现辨识结果偏离真实结果较远，猜想可能是采样周期设置过大，加之数据太小，于是出现了较大的偏差。