经典辨识方法之Hankel矩阵法在辨识过程中的增益倍数问题

最近在学习系统辨识这门课的时候产生了一个很大的疑问,就是经典辨识方法中的Hankel矩阵法辨识的结果与真实的传递函数之间总是存在一个与采样时间T0有关的倍数。
一开始我以为这个倍数是因为信号在采样和恢复过程中出现的衰减(1/T0倍),但后来发现直接使用离散的传递函数来进行辨识,仍然会出现1/T0倍的问题,这就说明不是采样的问题。老师解释说是因为c2d(matlab)函数的关系。。。,这个解答显然非常的不搭边,因为我们辨识的离散传递函数就已经和真实的离散传递函数有倍数之差了,上例子:
真实离散传递函数,采样时间为0.5S:
G(z)=0.3041z2+0.4988z+0.04133z30.8945z2+0.7885z0.04979G(z)=\frac{0.3041z^2+0.4988z+0.04133}{z^3-0.8945z^2+0.7885z-0.04979}
使用Hankel矩阵法辨识出来的离散传递函数:
G(z)=0.6083z2+0.9975z+0.08266z30.8945z2+0.7885z0.04979G(z)=\frac{0.6083z^2+0.9975z+0.08266}{z^3-0.8945z^2+0.7885z-0.04979}
可以发现二者分母多项式系数相同,分子多项式系数为1/T0倍。
发现问题,就一定要解决问题!!!
出现这种情况的原因还需要从Hankel矩阵的原理来讲(下面可能会有一些枯燥):
现在有一个n阶的脉冲传递函数需要我们辨识,可以将其写为:
G(z1)=b1z1+b2z2+...+bnzn1+a1z1+a2z2+...+anznG(z^{-1})=\frac{b_1z^{-1}+b_2z^{-2}+...+b_nz^{-n}}{1+a_1z^{-1}+a_2z^{-2}+...+a_nz^{-n}}
然后将其写成状态方程 的形式:
{x(k+1)=A(T)x(k)+B(T)u(k)y(k)=cx(k) \left\{ \begin{array}{c} x(k+1)=A(T)x(k)+B(T)u(k)\\ y(k)=cx(k) \ \end{array}\right.
上面的A(T)和B(T)表示这个矩阵是和采样时间有关的,故脉冲传递函数:
G(z1)=c(zIA)1bG(z^{-1})=c(zI-A)^{-1}b
也一定是和采样时间有关的。
接下来是一个巧妙的步骤,将G(z) 用脉冲响应表达:
(IAz)i=0(Az)i=i=0(Az)ii=0(Az)i+1=(Az)0+i=0(Az)i+1i=0(Az)i+1=I(I-\frac{A}{z})\sum_{i=0}^{\infty}(\frac{A}{z})^i= \sum_{i=0}^{\infty}(\frac{A}{z})^i-\sum_{i=0}^{\infty}(\frac{A}{z})^{i+1}\\=(\frac{A}{z})^0+ \sum_{i=0}^{\infty}(\frac{A}{z})^{i+1}-\sum_{i=0}^{\infty}(\frac{A}{z})^{i+1}=I
则:
(IAz)i=0(Az)i=Ii=0(Az)i=z(zIA)1(zIA)1=z1i=0(Az)i(I-\frac{A}{z})\sum_{i=0}^{\infty}(\frac{A}{z})^i=I\Longrightarrow \sum_{i=0}^{\infty}(\frac{A}{z})^i=z(zI-A)^{-1}\\\Longrightarrow(zI-A)^{-1}=z^{-1}\sum_{i=0}^{\infty}(\frac{A}{z})^i
即:
(zIA)1=A0z1+A1z2+A2z3+...(zI-A)^{-1}=A^0z^{-1}+A^1z^{-2}+A^2z^{-3}+...
在根据脉冲传递函数:
G(z1)=c(zIA)1b=g(1)z1+g(2)z2+g(3)z3+...G(z^{-1})=c(zI-A)^{-1}b=g(1)z^{-1}+g(2)z^{-2}+g(3)z^{-3}+...
这样就可以得到an,bn和g(n)的关系了:
G(z1)=b1z1+b2z2+...+bnzn1+a1z1+a2z2+...+anznb1z1+b2z2+...+bnzn=G(z1)(1+a1z1+a2z2+...+anzn)(b1z1+b2z2+...+bnzn)=g(1)z1+[g(2)+a1g(1)]z2+...+[g(n)+i=1n1(g(i)ani]zn+[g(n+1)+i=1n1(g(i)an+1i]zn1+...+[g(2n)+i=12n1(g(i)a2ni]z2nG(z^{-1})=\frac{b_1z^{-1}+b_2z^{-2}+...+b_nz^{-n}}{1+a_1z^{-1}+a_2z^{-2}+...+a_nz^{-n}}\\\Longrightarrow\\b_1z^{-1}+b_2z^{-2}+...+b_nz^{-n}=G(z^{-1})*{(1+a_1z^{-1}\\+a_2z^{-2}+...+a_nz^{-n}})\\\Longrightarrow\\(b_1z^{-1}+b_2z^{-2}+...+b_nz^{-n})=g(1)z^{-1}+[g(2)+a_1g(1)]z^{-2}+...\\+[g(n)+\sum_{i=1}^{n-1}(g(i)a_{n-i}]z^{-n}+[g(n+1)+\sum_{i=1}^{n-1}(g(i)a_{n+1-i}]z^{-n-1}+...\\+[g(2n)+\sum_{i=1}^{2n-1}(g(i)a_{2n-i}]z^{-2n}
比较幂级数,构造Hankel矩阵:
经典辨识方法之Hankel矩阵法在辨识过程中的增益倍数问题
这样就可以求出来分子分母的多项式系数了,感觉是不是很完美…(实际上辨识出来的函数确实很接近!!!)
但这个地方会有问题:
1.上面我们求出来的G(z)=c(zIA)1BG(z)=c(zI-A)^{-1}B是带有带有采样时间T的,而且应该在分母上(具体的没有进行推,但可以用个例子试试,不过要注意是离散的阵),然后使用长除法可以将其变化为g(n)zng(n)z^{-n}的形式,故我们采样得到的数据是隐含有采样时间T的。在求分母多项式时,两边将T的影响消除,故没有影响,但是在求分子多项式式时,就出现了1/T的问题,所以我们求出来以后,要*T进行补偿。
2.只有在g(n)是无穷多项的时候才能完全等于G(z),而我们只用了n项,这就产生了截断误差,故是一个误差来源。

上面就是我对这个问题的一点看法,希望大家有不同见解的及时交流!!!
(其实,我认为Hankel矩阵在判别系统的阶次时是一个非常不错的方法,下次我会简单介绍一下,大家关注一下我哈!!!手动狗头)