多速率系统的高速频率响应识别法和装置

摘要

本发明提供多速率系统的高速频率响应识别法和装置，通过FIR模型表示控制对象，基于M系列信号的周期M

说明书

技术领域

本发明的实施方式涉及多速率系统的高速频率响应识别法和高速频率 响应识别装置，涉及例如在控制输出的采样周期为控制输入的采样周期的 偶数倍的多重输入多速率系统中，高速地识别控制对象在控制输入侧的采 样周期下的频率响应的算法。

背景技术

在磁盘装置的头定位控制系统等所用的多重输入多速率系统中，有时 需要在控制系统设计时识别控制对象在控制输入侧的采样周期下的频率响 应(传递函数)。在多重输入多速率系统中，由于控制输出(观测输出) 被减少(降采样)，所以一般的频率响应的识别方法不能够直接应用。因 而，需要考虑了多重输入多速率系统的特征的识别算法。

发明内容

作为本发明的一方式的多速率系统的高速频率响应识别法，识别控制 输出的采样周期为控制输入的采样周期的偶数P倍的多重输入多速率系统 的、控制对象在控制输入侧的采样周期下的频率响应，包括：M系列信号 生成步骤、脉冲响应推定步骤、频率响应识别步骤。

前述M系列信号生成步骤，基于作为模拟白色信号的一种的M系列 信号的周期Mp，生成数据取得长度Mp×P-1量的M系列信号。

前述脉冲响应推定步骤，基于前述生成的M系列信号和输入前述生成 的M系列信号而从前述控制对象得到的输出数据，计算前述控制对象的脉 冲响应推定值。

前述频率响应识别步骤，通过对前述脉冲响应推定值进行离散傅立叶 变换，识别前述控制对象的频率响应。

附图说明

图1是表示本发明的实施方式的识别算法中的数据取得长度L的图。

图2是表示本发明的实施方式的执行步骤的流程图。

图3是表示计算机模拟中的输入输出关系的框线图。

图4是表示本发明的实施方式的满足了数学式(15)的情况下的频率 响应识别结果的图。

图5是表示本发明的实施方式的不满足数学式(15)的情况下的频率 响应识别结果的图。

图6是本发明的实施方式所涉及的装置的功能框图。

具体实施方式

本发明的实施方式所涉及的识别算法，作为识别输入使用作为模拟白 色信号的一种的M系统信号，作为识别模型使用FIR(Finite Impulse Response，有限脉冲响应)模型。此外，以多速率比为偶数P倍、即控制 输出的采样周期为控制输入的采样周期的2、4、6...倍为条件。

在此，已知多重输入多速率系统的传递函数由以下数学式(1)提供。

$Y\left(z\right)=\frac{1}{P}\underset{n=0}{\overset{P-1}{\Sigma }}{G}_0\left(z{\phi }^n\right)·U\left(z{\phi }^n\right),\phi =\exp \left(j\frac{2\pi }{P}\right)···\left(1\right)$

在数学式(1)中，Y(z)为输出的传递函数，U(z)为输入的传递函数， G₀(z)为控制对象的在输入侧的采样周期下的传递函数，P为多速率比。 在本实施方式中应该识别的控制对象的传递函数为数学式(1)中的G₀(z)。 数学式(1)，例如在多速率比P＝2时，成为下式。

$Y\left(z\right)=\frac{1}{2}{G}_0\left(z\right)U\left(z\right)+\frac{1}{2}{G}_0(-z)U(-z)···\left(2\right)$

以下，假定使用多速率比P＝2，说明本实施方式的识别算法。

首先，在本实施方式中，使用未知参数将G₀(z)的识别 模型设定为以下数学式(3)的FIR模型。

${\hat{G}}_0\left(z\right)={\hat{b}}_0+{\hat{b}}_1{z}^{-1}={\hat{b}}_2{z}^{-2}+···+{\hat{b}}_n{z}^{-n}···\left(3\right)$

此外，如果将输入的时间序列设定为u[0]、u[1]、...、u[m]，则数学 式(2)的U(z)与数学式(3)同样地，能够由以下的数学式(4)的FIR 模型表示。

U(z)＝u[0]+u[1]z^-1+u[2]z^-2+…+u[m]z^-m

…………………………(4)

进而，成为下式。

${\hat{G}}_0(-z)={\hat{b}}_0-{\hat{b}}_1{z}^{-1}+{\hat{b}}_2{z}^{-2}-···+{(-1)}^n{\hat{b}}_n{z}^{-n}$

U(-z)＝u[0]-u[1]z^-1+u[2]z^-2-…+(-1)^mu[m]z^-m

…………………(5)

在此，将数学式(2)的G₀(z)置换为的情况下的输出、即模型 输出的各时刻k＝0、1、2、...的时间序列能够利用数学式(3)-(5) 如以下那样通过单纯的卷积来计算。

(在k＝0时)

$\hat{y}\left[0\right]=\frac{1}{2}\left(\begin{array}{ccccc}{\hat{b}}_0& {\hat{b}}_1& {\hat{b}}_2& ···& {\hat{b}}_n\end{array}\right)·{\left(\begin{array}{cccccccc}u\left[0\right]& 0& 0& 0& 0& 0& ···& 0\end{array}\right)}^T$

$+\frac{1}{2}\left(\begin{array}{ccccc}{\hat{b}}_0& -{\hat{b}}_1& {\hat{b}}_2& ···& {(-1)}^n{\hat{b}}_n\end{array}\right)·{\left(\begin{array}{cccccccc}u\left[0\right]& 0& 0& 0& 0& 0& ···& 0\end{array}\right)}^T$

$={\hat{b}}_{0}u\left[0\right]$

(在k＝1时)

$\hat{y}\left[1\right]=\frac{1}{2}\left(\begin{array}{ccccc}{\hat{b}}_0& {\hat{b}}_1& {\hat{b}}_2& ···& {\hat{b}}_n\end{array}\right)·{\left(\begin{array}{cccccccc}u\left[1\right]& u\left[0\right]& 0& 0& 0& 0& ···& 0\end{array}\right)}^T$

$+\frac{1}{2}\left(\begin{array}{ccccc}{\hat{b}}_0& -{\hat{b}}_1& {\hat{b}}_2& ···& {(-1)}^n{\hat{b}}_n\end{array}\right)·{\left(\begin{array}{cccccccc}-u\left[1\right]& u\left[0\right]& 0& 0& 0& 0& ···& 0\end{array}\right)}^T$

$=0$

(在k＝2时)

$\hat{y}\left[2\right]=\frac{1}{2}\left(\begin{array}{ccccc}{\hat{b}}_0& {\hat{b}}_1& {\hat{b}}_2& ···& {\hat{b}}_n\end{array}\right)·{\left(\begin{array}{cccccccc}u\left[2\right]& u\left[1\right]& u\left[0\right]& 0& 0& 0& ···& 0\end{array}\right)}^T$

$+\frac{1}{2}\left(\begin{array}{ccccc}{\hat{b}}_0& -{\hat{b}}_1& {\hat{b}}_2& ···& {(-1)}^n{\hat{b}}_n\end{array}\right)·{\left(\begin{array}{cccccccc}u\left[2\right]& u\left[1\right]& u\left[0\right]& 0& 0& 0& ···& 0\end{array}\right)}^T$

$={\hat{b}}_{0}u\left[2\right]+{\hat{b}}_1\left[1\right]+{\hat{b}}_{2}u\left[0\right]$

(在k＝3时)

$\hat{y}\left[3\right]=\frac{1}{2}\left(\begin{array}{ccccc}{\hat{b}}_0& {\hat{b}}_1& {\hat{b}}_2& ···& {\hat{b}}_n\end{array}\right){·\left(\begin{array}{cccccccc}u\left[3\right]& u\left[2\right]& u\left[1\right]& u\left[0\right]& 0& 0& ...& 0\end{array}\right)}^T$

$+\frac{1}{2}\left(\begin{array}{ccccc}{\hat{b}}_0& -{\hat{b}}_1& {\hat{b}}_2& ···& {(-1)}^n{\hat{b}}_n\end{array}\right)·{\left(\begin{array}{cccccccc}-u\left[3\right]& u\left[2\right]& -u\left[1\right]& u\left[0\right]& 0& 0& ···& 0\end{array}\right)}^T$

$=0$

(在k＝4时)

$\hat{y}\left[4\right]=\frac{1}{2}\left(\begin{array}{ccccc}{\hat{b}}_0& {\hat{b}}_1& {\hat{b}}_2& ···& {\hat{b}}_n\end{array}\right)·{\left(\begin{array}{cccccccc}u\left[4\right]& u\left[3\right]& u\left[2\right]& u\left[1\right]& u\left[0\right]& 0& ···& 0\end{array}\right)}^T$

$+\frac{1}{2}\left(\begin{array}{ccccc}{\hat{b}}_0& -{\hat{b}}_1& {\hat{b}}_2& ···& {(-1)}^n{\hat{b}}_n\end{array}\right)·{\left(\begin{array}{cccccccc}u\left[4\right]& -u\left[3\right]& u\left[2\right]& -u\left[1\right]& u\left[0\right]& 0& ···& 0\end{array}\right)}^T$

$={\hat{b}}_{0}u\left[4\right]+{\hat{b}}_{1}u\left[3\right]+{\hat{b}}_{2}u\left[2\right]+{\hat{b}}_{3}u\left[1\right]+{\hat{b}}_{4}u\left[0\right]$

(以下同样)。

根据以上，可以理解，模型输出的时间序列在k为奇数时成为输 出减少(降采样)的响应。

为了推定识别模型的未知参数只要解决将输入信号 u[0]、u[1]、...、u[m]施加于实际的硬件(例如磁盘装置的头定位控制系统) 的情况下的输出的时间序列与模型输出的时间序列的误差e[k]的平方 和最小化的问题(最小二乘法)即可。在多速率比P＝2的多重输入多速率 系统中，在k为奇数时y[k]、都减小、成为0，所以e[k](k＝0、1、2、...) 成为以下数学式(6)。

$e\left[0\right]=y\left[0\right]-\hat{y}\left[0\right]=y\left[0\right]-{\hat{b}}_{0}u\left[0\right]$

$e\left[1\right]=y\left[1\right]-\hat{y}\left[1\right]=0$

$e\left[2\right]=y\left[2\right]-\hat{y}\left[2\right]=y\left[2\right]-\left({\hat{b}}_{0}u\right[2]+{\hat{b}}_{1}u[1]+{\hat{b}}_{0}u[0\left]\right)$

$e\left[3\right]=y\left[3\right]-\hat{y}\left[3\right]=0$

$e\left[4\right]=y\left[4\right]-\hat{y}\left[4\right]=y\left[4\right]-\left({\hat{b}}_{0}u\right[4]+{\hat{b}}_{1}u[3]+{\hat{b}}_{2}u[2]+{\hat{b}}_{3}u[1]+{\hat{b}}_{4}u[0\left]\right)$

……………(6)

从而，应该解开的最小二乘法成为下式。

$\underset{l=0}{\overset{N}{\Sigma }}e{\left[2l\right]}^2\to \min ···\left(7\right)$

数学式(7)的解，使用包含未知参数和输入u[0]、u[1]、u[2]、... 的以下数学式(8)的下述向量：

$\hat{\theta }={\left(\begin{array}{ccccc}{\hat{b}}_0& {\hat{b}}_1& {\hat{b}}_2& ···& {\hat{b}}_n\end{array}\right)}^T···\left(8\right)$

u[k]＝[u[k] u[k-1] u[k-2]…u[k-n]]^T

，由以下数学式(9)提供。

$\hat{\theta }={\left(\underset{l=0}{\overset{N}{\Sigma }}u{\left[2l\right]u\left[2l\right]}^T\right)}^{-1}\left(\underset{l=0}{\overset{N}{\Sigma }}u\right[2l\left]y\right[2l\left]\right)···\left(9\right)$

接着，考虑实际计算数学式(9)。在本实施方式中，对于u[k]使用作 为模拟白色信号的一种的M系列，所以若通过M系列的周期M_p将数学 式(9)的总和次数N选择为：

N＝M_p-1……………………………………………(10)

，则数学式(9)右边的第一个括弧内

$\Lambda =\underset{l=0}{\overset{N}{\Sigma }}u{\left[2l\right]u\left[2l\right]}^T···\left(11\right)$

成为下式。

Λ＝u[0]u[0]^T+u[2]u[2]^T+…+u[2(M_p-2)]u[2(M_p-2)]^T+u[2(M_p-1)]u[2(M_p-1)]^T

                                          …………………(12)

在此，若考虑M系列的周期性，则数学式(12)能够变形为：

$\Lambda =u{\left[0\right]u\left[0\right]}^T+{u\left[2\right]u\left[2\right]}^T+···+u{[M_p-1]u[M_p-1]}^T$

$+u{\left[1\right]u\left[1\right]}^T+u{\left[3\right]u\left[3\right]}^T+···+u{[M_p-2]u[M_p-2]}^T$

$=\underset{l=0}{\overset{M_p-1}{\Sigma }}u{\left[l\right]u\left[l\right]}^T···\left(13\right)$

，可以看出成为M系列1周期量的自相关矩阵。因而，根据白色信号的自 相关的性质，数学式(13)的非对角项能够如以下数学式(14)那样近似 地设定为0。

数学式(14)的σ_u²是u[k]的方差。根据以上，数学式(9)右边的逆矩阵 运算能够置换为单纯的除法，能够大幅地降低计算负荷。也就是说，不需 要进行计算负荷大的逆矩阵运算。使数学式(14)成立所需最小限的识别 实验时的取得数据长度L成为：

L＝M_pP-1    ………………………………………………(15)

。将识别实验时的取得数据长度的概念示于图1。在图1中，示出多速率 比P＝2、M系列的周期Mp＝7的情况下的数据取得长度。

根据以上，多速率比P＝2的情况下的未知参数推定算法成为：

$\hat{\theta }=\frac{\left(\underset{l=0}{\overset{M_p-1}{\Sigma }}u\right[2l\left]y\right[2l\left]\right)}{{\sigma }_u^2}···\left(16\right)$

。一般的偶数倍多速率比(P＝2、4、6、...)的情况成为下式。

$\hat{\theta }=\frac{\left(\underset{l=0}{\overset{M_p-1}{\Sigma }}u\right[\mathrm{Pl}\left]y\right[\mathrm{Pl}\left]\right)}{{\sigma }_u^2}···\left(17\right)$

以上的数学式(15)、(17)是本实施方式的基本的算法。此外，根 据数学式(17)计算的是G₀(z)的脉冲响应推定值(FIR模型的分子多 项式系数)，所以如果对进行离散傅立叶变换，则能够求出G₀(z)的频 率响应。离散傅立叶变换能够使用FFT(Fast Fourier Transform，快速傅 立叶变换)容易地执行。

以上所示的本实施方式的算法，由于将识别输入信号规定为M系列信 号，所以一般与正弦波扫频法比较能够缩短识别实验所需的时间。此外， 由于数学式(17)由单纯的加法及乘除法构成，所以能够以简单的描述安 装于应用本实施方式的硬件所搭载的微处理器和/或外部的计算机。

表示本实施方式的执行步骤的流程图为图2。按照该执行步骤的处理， 由硬件所搭载的微处理器和/或外部的计算机执行。

在步骤S1，判定取得的数据长度L是否达到了Mp×P-1。

在未达到Mp×P-1时，生成下一时刻的M系列信号(S2)，将M系 列信号输入于识别对象(例如磁盘装置)(S3)。另外，M系列信号的生 成算法预先提供。

观测控制对象(识别对象)的输出(S4)，将M系列和输出数据保存 于内部存储器或外部存储装置等存储部(S5)。反复以上的步骤S1～S4， 直至数据长度L达到Mp×P-1为止。

若数据长度L达到了Mp×P-1，即若取得了数据长度L量的M系列信 号及输出数据，则计算控制对象的脉冲响应推定值。具体地，根据数学式 (17)计算控制对象的识别模型、即FIR模型的脉冲响应推定值并且，对进行离散傅立叶变换(S7)。离散傅立叶变换的结果，为应该 求取的识别对象的频率响应。

在图6中示出本实施方式所涉及的装置的功能框图。

M系列信号生成部11按数据取得长度L＝Mp×P-1的量，生成作为模 拟白色信号的一种的M系列信号。

将所生成的M系列信号的各时刻的信号分别输入于识别对象12，并 从控制对象(识别对象)12输出输出数据。

存储部15存储由生成部11生成的M系列信号和与各时刻的信号对应 地从识别对象12得到的输出数据。

脉冲响应推定部13基于存储于存储部15的M系列信号及输出数据， 计算控制对象12的脉冲响应推定值。该计算根据数学式(A)进行。

脉冲响应推定值$\hat{\theta }=\frac{\left(\underset{l=0}{\overset{M_p-1}{\Sigma }}u\right[\mathrm{Pl}\left]y\right[\mathrm{Pl}\left]\right)}{{\sigma }_u^2}···\left(A\right)$

数学式(A)与前述的数学式(17)相同。如前所述，u[k]是在时刻k (第k个)取得的M系列信号的分量，y[k]是在将u[k]设定为识别对象12 的输入时得到的输出数据，σ_u²是u[k]的方差。

傅立叶变换部14通过对进行离散傅立叶变换，识别控制对象的频率 响应。

以上，根据本实施方式，在多重输入多速率系统中，能够通过高速且 简单的运算识别控制对象在控制输入侧的采样周期下的频率响应。

以下，通过计算机模拟表示本发明的实施例。

首先，将假定了具有多个谐振模式的机械系统的、以下数学式(18) 的离散时间序列传递函数设定为计算机模拟中的识别对象。

$P\left(z\right)=\frac{\left(\begin{array}{c}-0.009539z^9+0.08884z^8+0.1452z^7+0.3556z^6+0.4401z^5\\ +0.4295z^4+0.3348z^3+0.1316z^2+0.07566z-0.005744\end{array}\right)}{\left(\begin{array}{c}{z}^{10}-0.4028z^9+1.818z^8-2.063z^7+1.624z^6-2.951z^5\\ +1.453z^4-1.915z^3+1.352z^2-0.3971z+0.6743\end{array}\right)}···\left(18\right)$

采样周期设定为50[μs]。对该识别对象以采样周期50[μs]输入M系列信 号，通过计算机模拟得到以多速率比P＝2降采样后的输出(采样周期100[μ s])。在图3中示出表示此时的识别对象的输入输出关系的框线图。M系 列的移位寄存器阶数设定为13阶(周期M_p＝8191)。此外，由数学式(17) 推定的脉冲响应的长度设定为1024。这相当于在数学式(8)中设定为 n＝1023，即：$\hat{\theta }={\left(\begin{array}{ccccc}{\hat{b}}_0& {\hat{b}}_1& {\hat{b}}_2& ···& {\hat{b}}_{1023}\end{array}\right)}^T$

u[k]＝[u[k] u[k-1] u[k-2]…u[k-1023]]^T

…………………(19)

。最后，对执行数学式(17)而得到的脉冲响应推定值执行1024点的FFT， 得到识别对象的频率响应。在此，在根据本实施方式的数学式(15)将输 入输出的取得数据长度设定为

L＝M_pP-1

＝8191×2-1………………………………………………(20)

＝16381

的情况下和在不根据数学式(15)而将输入输出的取得数据长度设定为任 意设定的值

L＝13000   …………………………………………………(21)

的情况下进行数学式(17)的计算，检查取得数据长度对识别精度产生的 影响。

将在以上的设定下执行计算机模拟的结果(识别对象的频率响应)示 于图4及图5。在图4及图5中，实线是所识别的频率响应，点线是识别 对象的真实的频率响应。此外，图中在5000[Hz]处示出的虚线表示采样周 期100[Hz]时的奈奎斯特频率。与本实施方式的识别方法不同，在不正面 考虑多速率系统的特性的识别方法中，难以识别奈奎斯特频率以上(虚线 右侧)的频率响应。

从图4可以看出，在设定了满足本实施方式的数学式(15)的取得数 据长度L的情况下，所识别的频率响应与真实的频率响应高度一致，在奈 奎斯特频率以上的频带也能够正确地识别。另一方面，从图5可以看出， 在设定了不满足数学式(15)的任意的取得数据长度L的情况下，虽然可 以看出频率响应的大致的一致，但是与满足数学式(15)的情况相比，与 真实的频率响应的误差较大。通过以上的计算机模拟所实现的实施例，验 证了本实施方式(15)、(17)的效果。

另外，对于如磁盘装置的头定位控制系统那样在奈奎斯特频率以上具 有谐振模式的控制对象，通过本实施方式的算法得到的控制输入侧的采样 周期下的频率响应能够用于为了将谐振模式稳定化而使用的多速率陷波滤 波器的设计。