第九章-模型阶次的确定

一、Hankel矩阵法
可控可观的SISO系统，利用相关矩阵 法——脉冲响应序列——参数模型 时，需要确定阶次。 方法：利用Hankel矩阵的秩确定模型 阶次。
{g(k)}, k = 1,", L为脉冲响应序列，其构 成的Hankel矩阵为
⎡ g(k) g(k +1) " g(k + l −1) ⎤
⎢
⎢
H
(l,
k
)
=
⎢ ⎢
g(k +1) #
g(k + 2) "
#
#
g(k + l)
⎥ ⎥
#⎥
⎡c⎤
⎢ ⎣
g(k
+
l
−
1)
g(k + l)
"
g
(k
+
2l
−
2)
⎥ ⎦
l×l
⎢ H (l, k) = ⎢
⎢
cA #
⎥ ⎥ Ak−1[b, Ab,", Al−1b] ⎥
g(k ) = cAk−1b
⎢⎣cAl
−1
⎥ ⎦
(
L
−
n
)
⎥ ⎦
⎪ ⎪ ⎪ ⎪En ⎪ ⎪⎩
=
⎡e(0)
⎢ ⎢
e(1)
⎢
⎢ ⎣
e(
L
− 1)
e(−1) e(0)
⋅⋅⋅⋅⋅⋅ e(L − 2)
⋅ ⋅ ⋅ e(1 − n) ⎤
⋅⋅⋅e(2 − n)
⎥ ⎥
⎥
⋅⋅⋅
e(L
−
n)
⎥ ⎦
其中，L为数据长度，则（13.4.1）式可写成
zn = Hnθn + vn
⋅ ⋅ ⋅ z(2 − n) # u (1) #
⋅ ⋅ ⋅ z( L − n)# u ( L − 1)
u(0) ⋅⋅⋅⋅⋅⋅ u(L − 2)
⋅ ⋅ ⋅ u (1 − n) ⎤
⋅⋅⋅u(2
−
n)
⎥ ⎥
⎥
⋅
⋅
⋅
u
(
L
−
n)
⎥ ⎦
⎪ ⎪⎪
Δ ⎡⎣
Zn
#
Un
]
⎨ ⎪
Δ ⎡⎣
Yn + En
#
Un
]
⎪ ⎡ y(0)
τ n
⋅
⋅
⎤ ⎥ ⎥
⎢⎣
U
τ n
⎥⎦
v n0
⎫ ⎪ ⎬ ⎪ ⎭
考虑到v(k)是零均值的白噪声，当L → ∞时，上式右边第三、四项为零，则
[ ] P
lim
L→∞
V1
(nˆ)
=
P
lim(
L→∞
1 L
~x τ
~x )
+
P
lim(
L→∞
1 L
vτ n0
vn0
)
因参数估计值 θˆnˆ有下列性质
P
lim
L→∞
θˆnˆ
⎧⎪ A(z−1
⎨ ⎪⎩
B(
z
−1
) )
= =
1+ a1z−1 + ⋅⋅⋅ + an b1z−1 + b2 z−2 + ⋅ ⋅⋅
z−n + bn
z−n
⎪
⎡ z(0) z(−1) ⋅ ⋅ ⋅ z (1 − n) # u (0) u (−1)
⎪
⎪ ⎪
H
n
=
⎢ ⎢
z
(1)
⎢
⎪ ⎪
⎢ ⎣
z
(
L
−
1)
z(0) ⋅⋅⋅⋅⋅⋅ z(L − 2)
1
即⎧
⎡c⎤
⎪ ⎪
⎢ ⎢
cA
⎥ ⎥
⎪
⎢#⎥
⎪ ⎨
rank
⎢⎢cAn0
−1
⎥ ⎥
= n0
⎪ ⎪ ⎪
⎢#⎥
⎢ ⎣
cAl
⎥ ⎦
⎪⎩rank[b, Ab,", An0−1b,", Alb] = n0
将脉冲响应与状态方程 ( A, b, c)的关系
g(k) =cAk−1b
代入Hankel矩阵，有
有
⎡ g(k) g(k +1) " g(k + l −1) ⎤
H
(l,
k
)
=
⎢ ⎢
g(k +1) #
g(k + 2) "
#
#
g(k + l)
⎥ ⎥
#⎥
⎢ ⎣
g
(k
+
l
−1)
g(k + l)
"
g
(k
+
2l
−
2)
⎥ ⎦l×l
l 决定Hankel矩阵的维数；k 是 1 ≤ k ≤ (L − 2l + 2) 的
任意整数，决定由哪些脉冲响应序列组成Hankel
矩阵。 l ≥ n0 （系统的真实阶次）时，Hankel矩
第九章 模型阶次的确定
前面各章讨论差分方程参数的辨识 方法时，都假定差分方程的阶是已 知的。在一些实际问题中，模型的 阶可以按理论推导获得，而在另一 些实际问题中，模型的阶次却无法 用理论推导的方法确定，需要对模 型的阶进行辨识。
几种常用的模型阶的确定方法
Hankel矩阵法 利用行列式比法定阶 利用残差的方差估计模型的阶次
在最小二乘意义下，参数θ
的估计值为
n
^
θn
=
(
H
τ nˆ
H
n
)−1
H
τ nˆ
zn0
其中，nˆ表示模型的阶次估计值，n0为过程的真实阶次。于是模型的输出残差可写成
式中
znˆ = zn0 − Hnˆθˆnˆ = x + vn0 ~x = θ H n 0 n 0 − H nˆθˆnˆ
θn0 表示过程模型的真实参数，那么残差的方差为
rank(H (l, k)) = n0 ,l ≥ n0 , ∀k
此式表明当 l 大于系统的真实阶次 n0 时，
Hankel矩阵的秩仍然等于n0 。 判定方法：通过对不同的 l 值，判断
Hankel矩阵的奇异性来确定系统的阶次。
1、无噪声情况
脉冲序列 g(1), g(2),", g(L) 不含噪声的情况。
2、弱噪声情况
当脉冲序列 g(1), g(2),", g(L) 含有噪声，这样即 使 l 已增加到了n0 = nˆ ，但对所有的k ， Hankel矩阵的行列式不会都绝对为零。这 样难以按照无噪声的情况来确定模型的阶 次。
如果脉冲响应所含的噪声比较小，则可引
进 Hankel 矩 阵 行 列 式 的 平 均 比 值 Dl 值 来观察Hankel矩阵是否已由非奇异变成奇
v
n0
)
+
2θ
τ n0
P
⎧
lim
⎪⎪ 1 ⎨
L→ ∞ ⎪ L
⎪⎩
⎡ ⎢
Y
τ n0
+
Eτ n0
⎢⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅
⎢ ⎣
Uτ n0
⎤ ⎥ ⎥ ⎥ ⎦
v n0
⎫ ⎪⎪ ⎬ ⎪ ⎪⎭
−
2P
lim
L→ ∞
(
θˆ
τ nˆ
)P
⎧ ⎪1 lim ⎨ L→ ∞ ⎪ L ⎩
⎡ ⎢ ⎢
Y
τ n
⋅⋅⋅
⋅
+ ⋅⋅
E ⋅⋅
现在要讨论的定阶方法则需要在获得模 型参数估计值之后，求得模型残差序列，并 借以统计假设检验方法对残差的方差进行显 著性检验来确定模型的阶次，它和参数辨识 方法是密切相关的。
1、 残差方差分析
考虑用如下模型描述的SISO过程
A(z−1)z(k ) = B(z−1)u(k ) + v(k )
其中，u(k )和v(k )是过程的输入输出变量；v(k )是均值为零，
=
⎢ ⎢
y(n +1) #
y(n) " y(2) u(n +1)
#
%#
#
u(n)
"
u(2)
⎥ ⎥
#
% #⎥
⎢ ⎢⎣ y(n + L −1)
y(n + L − 2)
"
y(L) u(n + L −1)
u(n + L − 2)
"
⎥ u(L)⎥⎦
= [Yn |Un ]
det(H (nˆ)) = 0
则应取 nˆ −1 作为系统的阶次。
方差为σ
2的不相关随机噪声；且
v
⎪⎧ A( ⎨ ⎪⎩B(
z z
−1 −1
) )
= =
1+ b1z
a1z−1 + −1 + b2z
⋅⋅⋅ + an −2 + ⋅⋅⋅
z−n + bn
z
−n
z(k )
=
B( z −1) A( z −1 )
u(k)
+
1 A( z −1 )
Δ
v(k) =
y(k)
+
e(k )
置
⎡ ρ(k) ρ(k +1) " ρ(k + l −1) ⎤
⎢
H
(l,
k)
=
⎢ ⎢
ρ(k + 1) #
ρ(k + 2) " #%
ρ(k + l)
⎥ ⎥
⎥
⎢ ⎣
ρ
(k
+
l
− 1)
ρ(k + l)
"
ρ(k + 2l − 2)⎥⎦
ρ(k) = Rg (k) Rg (0)
∑ Rg (k)
=
L
1 −
k
L−k i =1
g (i ) g (i
+
k)
2
例 系统的脉冲响应序列表
方法一： 矩阵H（2，k）行列式的平均值＝0.00087872 矩阵H（3，k）行列式的平均值＝-0.00029311 矩阵H（4，k）行列式的平均值＝-3.214×10-7 矩阵H（5，k）行列式的平均值＝-5.709×10-9 D2＝2.998，D3＝913.1，D4＝64.2，因此可确定 系统的阶数为3。 方法二、求出脉冲响应序列的相关系数值为
3
受计算误差等因素的影响，真正让 det(H (nˆ)) = 0
是比较困难的。为了提高阶次估计值的准 确性，利用行列式比来确定模型的阶次。 定义 DR(nˆ) = det(H (nˆ))
det(H (nˆ +1))
当 nˆ 从1开始逐一增加时，若DR(nˆ)较 有 DR(nˆ +1)
显著的增加，则此时的 nˆ 可认为比较接近
真实阶次，即应取 n0 = nˆ 。
特点：不利用参数的估计值，而是采用系 统输入和输出的量测值，因而在开始估计 参数之前就能确定系统的阶次。
三、 利用残差的方差估计模型的阶次
上面讨论的Hankel矩阵判秩定阶法或行 列式比定阶法在获得模型参数估计值之前就 可预先确定模型的阶次，它基本上与参数识 别方法无关。
⎪ ⎪⎪Yn
=
⎢ ⎢
y
(1)
⎢
y ( − 1) y(0) ⋅⋅⋅⋅⋅⋅
⋅ ⋅ ⋅ y(1 − n) ⎤
⋅⋅⋅ y(2 − n)
⎥ ⎥
⎥
z(k )
=
B( z −1) A( z −1 )
u(k) +
1 A( z −1 )
Δ
v(k ) =
y(k)
+
e(k )
⎪ ⎪
⎢ ⎣
y
(
L
−
1)
y(L − 2)
⋅
⋅
⋅
y
首先，构造Hankel矩阵 H (l, k)
然后，l 逐一增加，计算 det[H (l, k)]，若 当 l = 1,", nˆ ，对所有 k ， 有 det[H (l, k)] ≠ 0；而当 l ≥ nˆ +1 时，对所 有 k，有det[H (l, k)] = 0，这说明Hankel矩阵 在 l = nˆ +1 处由非奇异阵变成奇异阵，据 此可以判断系统模型的阶次是 n0 = nˆ 。
二、利用行列式比法定阶
考虑无观测噪声时的情况，系统模型为
y(k) = −a1y(k −1) −" − an y(k − n) + b1u(k −1) +"+ bnu(k − n)
利用输入输出数据确定模型阶次 n 。
⎡ y(n)
y(n −1) " y(1) u(n)
u(n −1) " u(1) ⎤
⎢
置H n
vτ n0
vn0
)
=
σ
2 v
,
nˆ
≥
n0
V1 (nˆ)
异。
∑ 1
L−2l+2
det[H (l, k)]
∑ Dl
=
L − 2l 1
+ 2 k =1
L−2l
det[H (l +1, k)]
L − 2l k=1
当 Dl 从1逐渐增加时，不断计算 Dl 值，可 取 Dl 达到最大值的 l 作为模型的阶次，见 图。
Dl
l n0
3、强噪声情况
Hankle矩阵构成不能直接用脉冲响应序 列，需用脉冲响应序列的相关系数。除此 以外，其它计算与前面相同。
阵的秩等于 n0 ，即
。
rank(H (l, k)) = n0 ,l ≥ n0 ,∀k
证明：
考虑可观可控SISO系统
x(k +1) = Ax(k) + bu(k)
z(k) = cx(k)
其可控性矩阵和可观测性矩阵的秩均为n0 。
这意味着可观测性矩阵(或可控性矩阵)中 的 n0 行或 n0列是线性独立的。若在可观 测性矩阵（或可控性矩阵）中增加若干行 向量 cA（l 或若干列 Alb ），其中 l > n0 ，则 可观测性矩阵（或可控性矩阵）的秩不会 改变。
det H (2,0) = 0.014937371
det H (3,0) = −0.00001282
det H (4,0) = −0.000000058
由行列式的值可知，系统模型的阶次可以 定为3阶，也可定为2阶，因为detH(3,0)已 经很小。由行列式比值可知，系统的阶次 可定为2阶。
说明：由以上两种定阶方法的计算结果来 看，这两种方法在定阶时是存在一定差异 的。
=
⎧ ⎪ ⎨ ⎪⎩=
≠ θn0 = θn0 [θ T n0 ,0T
nˆ nˆ ]
< n0 = n0 nˆ >
n0
P
lim(
L→∞
1 L
~x τ
~x
)
=
⎧= ⎩⎨>
0 0
nˆ ≥ n0 nˆ < n0
可见 nˆ > n0 时，输出残差趋于固定值
[ ] P
lim
L→∞
பைடு நூலகம்
V1
(nˆ)
=
P
lim(
L→∞
1 L
V1 (nˆ)
=
1 L
znτˆ znˆ
=
1 L
(
xτ
x
+
2θ
τ n0
H
τ n0
vn0
−
2θ
τ n
H
τ n
vn0
+
vτ n0
vn0
)
对上式运用定理 D.3，可得
[ ] P
lim
L→ ∞
V 1 ( nˆ )
=
P
lim ( 1 L→ ∞ L
~x τ ~x ) +
P
lim ( 1 L→ ∞ L
v
τ n0
ρ0 = 1, ρ1 = 0.88052126, ρ2 = 0.79025506 ρ3 = 0.72231277, ρ4 = 0.67060564, ρ5 = 0.62999127 ρ6 = 0.60107303, ρ7 = 0.5769552,"
以相关系数值为元素构造Hankel矩阵并计 算Hankel矩阵的行列式，得
A(z−1)z(k) = B(z−1)u(k) + v(k)
置 ⎧θ
⎪
n
⎪ zn
= =
[− a1, − a2 , ⋅ ⋅ ⋅, − an , b1, b2 , ⋅ ⋅ ⋅, bn ]τ [ z(1), z(2), ⋅ ⋅ ⋅, z( L)]τ
⎪ ⎪
vn
=
[v(1), v(2), ⋅ ⋅ ⋅, v( L)]τ