CN105960613A - 系统辨识装置 - Google Patents

系统辨识装置 Download PDF

Info

Publication number
CN105960613A
CN105960613A CN201480074854.2A CN201480074854A CN105960613A CN 105960613 A CN105960613 A CN 105960613A CN 201480074854 A CN201480074854 A CN 201480074854A CN 105960613 A CN105960613 A CN 105960613A
Authority
CN
China
Prior art keywords
dimension
matrix
identification
input
identifying device
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Pending
Application number
CN201480074854.2A
Other languages
English (en)
Inventor
齐藤光伯
金井百合夏
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Mitsubishi Electric Corp
Original Assignee
Mitsubishi Electric Corp
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Mitsubishi Electric Corp filed Critical Mitsubishi Electric Corp
Publication of CN105960613A publication Critical patent/CN105960613A/zh
Pending legal-status Critical Current

Links

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F7/00Methods or arrangements for processing data by operating upon the order or content of the data handled
    • G06F7/38Methods or arrangements for performing computations using exclusively denominational number representation, e.g. using binary, ternary, decimal representation
    • GPHYSICS
    • G05CONTROLLING; REGULATING
    • G05BCONTROL OR REGULATING SYSTEMS IN GENERAL; FUNCTIONAL ELEMENTS OF SUCH SYSTEMS; MONITORING OR TESTING ARRANGEMENTS FOR SUCH SYSTEMS OR ELEMENTS
    • G05B13/00Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion
    • G05B13/02Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion electric
    • G05B13/04Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion electric involving the use of models or simulators
    • G05B13/042Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion electric involving the use of models or simulators in which a parameter or coefficient is automatically adjusted to optimise the performance
    • G05B13/044Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion electric involving the use of models or simulators in which a parameter or coefficient is automatically adjusted to optimise the performance not using a perturbation signal
    • HELECTRICITY
    • H03ELECTRONIC CIRCUITRY
    • H03MCODING; DECODING; CODE CONVERSION IN GENERAL
    • H03M13/00Coding, decoding or code conversion, for error detection or error correction; Coding theory basic assumptions; Coding bounds; Error probability evaluation methods; Channel models; Simulation or testing of codes
    • H03M13/29Coding, decoding or code conversion, for error detection or error correction; Coding theory basic assumptions; Coding bounds; Error probability evaluation methods; Channel models; Simulation or testing of codes combining two or more codes or code structures, e.g. product codes, generalised product codes, concatenated codes, inner and outer codes
    • HELECTRICITY
    • H03ELECTRONIC CIRCUITRY
    • H03MCODING; DECODING; CODE CONVERSION IN GENERAL
    • H03M13/00Coding, decoding or code conversion, for error detection or error correction; Coding theory basic assumptions; Coding bounds; Error probability evaluation methods; Channel models; Simulation or testing of codes
    • H03M13/61Aspects and characteristics of methods and arrangements for error correction or error detection, not provided for otherwise
    • H03M13/611Specific encoding aspects, e.g. encoding by means of decoding

Abstract

为了能够从系统维数的决定排除反复试验,始终决定出最佳的系统维数,并辨识记述动态系统的线性离散时间系统,系统辨识装置将作为辨识对象的动态系统的脉冲响应、和系统维数的搜索范围作为输入,具备:直通项辨识部;块汉克尔矩阵生成部;奇异值分解部,通过块汉克尔矩阵的奇异值分解,输出第1正交矩阵、第2正交矩阵以及奇异值;系统维数决定部,根据第1正交矩阵、第2正交矩阵、奇异值以及搜索范围,针对属于搜索范围的各维数,辨识直通项以外的系统矩阵,根据基于系统矩阵以及直通项计算出的系统特性和实际的系统特性的比较,决定系统维数;和系统矩阵辨识部,根据第1正交矩阵、第2正交矩阵、奇异值以及系统维数,辨识直通项以外的系统矩阵。

Description

系统辨识装置
技术领域
本发明涉及系统辨识装置。
背景技术
作为以往的基于脉冲响应的系统辨识装置,在例如非专利文献1中公开了基于Ho-Kalman的方法的系统辨识装置。在该方法中,根据用线性离散时间系统(Ad,Bd,Cd,Dd)记述的动态系统的脉冲响应(G0,G1,G2,…),依据G0决定系统的直通项Dd,并且依据G1、G2、…生成块汉克尔矩阵Hkl。接下来,对块汉克尔矩阵Hkl进行奇异值分解,将具有有效值的奇异值的个数决定为系统维数,根据奇异值分解的结果和所决定的系统维数,计算扩展可观测矩阵Ok以及扩展可到达矩阵Cl。最后,根据扩展可观测矩阵Ok以及扩展可到达矩阵Cl,计算系统矩阵Ad、Bd、Cd,从而辨识记述动态系统的线性离散时间系统(Ad,Bd,Cd,Dd)。
另外,作为以往的基于脉冲响应的系统辨识装置的其他例子,在例如专利文献1中公开了对象建模(plant modeling)装置。在该对象建模装置中,对动态系统的脉冲响应(G0,G1,G2,…)应用上述Ho-Kalman的方法,而选择性地应用接下来的2种方法作为系统维数的决定方法。第1系统维数的决定方法是指如下方法:在图形终端中按照对数尺度显示奇异值和与该奇异值对应的次数的关系,作业者决定具有有效值的奇异值的个数、即系统维数。第2系统维数的决定方法是指如下方法:应用基于奇异值的变化率和观测噪声的评价函数,将使该评价函数成为最小的维数自动地决定为系统维数。
在这样的基于脉冲响应的系统辨识装置中,根据具有有效值的奇异值的个数来决定系统维数。
专利文献1:日本特开昭61-267102号公报
非专利文献1:片山徹著,“システム同定-部分空間法からのアプローチー(系统辨识-根据子空间法的方法)”,朝仓书店,2004年2月,p.102-107
发明内容
但是,根据上述以往的技术,根据现实的脉冲响应计算出的块汉克尔矩阵的奇异值为平滑的单调减少的情况较多,在该情况下,具有有效值的奇异值和成为可忽略的微小的值的奇异值的边界变得不明确。因此,在基于第1系统维数的决定方法的以往的系统辨识装置中,存在如下问题:系统维数的决定依赖于作业者的判断,未必始终决定出最佳的系统维数、或者在系统维数的决定中需要反复试验。
另外,作为应对这样的问题的方法,在专利文献1中,应用了第2系统维数的决定方法,但即使在该方法中,根据评价函数的提供方式而所决定的系统维数也发生变化,所以系统维数的决定依赖于评价函数的提供方式,即使在该方法中也未必始终决定出最佳的系统维数、或者关于系统维数的决定需要反复试验,尚未解决问题。
本发明是鉴于上述而完成的,其目的在于得到一种系统辨识装置,即使在根据现实的脉冲响应计算出的块汉克尔矩阵的奇异值为平滑的单调减少,具有有效值的奇异值和成为可忽略的微小的值的奇异值的边界变得不明确的情况下,也能够从系统维数的决定排除反复试验,始终决定出最佳的系统维数,进行记述动态系统的线性离散时间系统的辨识。
为了解决上述课题并实现目的,本发明的系统辨识装置将作为辨识对象的动态系统的脉冲响应、和被指定的系统维数的搜索范围作为输入,其特征在于,具备:直通项辨识部,根据所述脉冲响应,辨识并输出记述所述动态系统的线性离散时间系统的直通项;块汉克尔矩阵生成部,根据所述脉冲响应,生成并输出块汉克尔矩阵;奇异值分解部,通过从所述块汉克尔矩阵生成部输出的所述块汉克尔矩阵的奇异值分解,输出将述块汉克尔矩阵的左奇异矢量作为列矢量的第1正交矩阵、将所述块汉克尔矩阵的右奇异矢量作为列矢量的第2正交矩阵以及所述块汉克尔矩阵的奇异值;系统维数决定部,根据所述第1正交矩阵、所述第2正交矩阵、所述奇异值以及所述搜索范围,辨识针对属于该搜索范围的各维数的所述线性离散时间系统的系统矩阵中的所述直通项以外的系统矩阵,根据基于该系统矩阵以及所述直通项计算出的所述线性离散时间系统的系统特性和所述动态系统的实际的系统特性的比较,决定并输出系统维数;以及系统矩阵辨识部,根据所述第1正交矩阵、所述第2正交矩阵、所述奇异值以及从所述系统维数决定部输出的所述系统维数,辨识所述线性离散时间系统的系统矩阵中的所述直通项以外的系统矩阵,所述系统辨识装置将由所述直通项辨识部辨识出的所述直通项以及由所述系统矩阵辨识部辨识出的所述系统矩阵作为所述线性离散时间系统输出。
根据本发明,起到如下效果:在作为辨识对象的动态系统中,即使在根据现实的脉冲响应计算出的块汉克尔矩阵的奇异值为平滑的单调减少,具有有效值的奇异值和成为可忽略的微小的值的奇异值的边界变得不明确的情况下,也能够从系统维数的决定排除反复试验,始终进行最佳的系统维数的决定以及记述动态系统的线性离散时间系统的辨识。
附图说明
图1是示出实施方式1的系统辨识装置的整体结构的框图。
图2是示出在实施方式1的系统辨识装置中,系统输入输出的时间波形的概略图。
图3是示出在实施方式1的系统辨识装置中,块汉克尔矩阵的奇异值与维数的关系的概略图。
图4是示出在实施方式1的系统辨识装置中,系统维数决定部的内部结构和动作的框图。
图5是示出在实施方式1的系统辨识装置中,辨识出的线性离散时间系统的推测脉冲响应和系统的实际的脉冲响应的误差平方和范数与维数的关系的概略图。
图6是示出在实施方式2的系统辨识装置中,对动态系统进行了冲击振动的情况下的系统输入输出的时间波形的概略图。
图7是示出在实施方式2的系统辨识装置中,将对动态系统进行了冲击振动的情况下的系统输入输出变换为该系统的脉冲响应的步骤的框图。
图8是示出在实施方式2的系统辨识装置中,系统维数决定部的内部结构和动作的框图。
图9是示出在实施方式3的系统辨识装置中,将动态系统的频率响应变换为该系统的脉冲响应的步骤的框图。
图10是示出实施方式4的系统辨识装置的整体结构的框图。
(符号说明)
1:直通项辨识部;2:块汉克尔矩阵生成部;3:奇异值分解部;4、4a:系统维数决定部;5:系统矩阵辨识部;10、70:系统辨识装置;11:系统输入;12:系统输出;13:系统输入阈值;21:(理想的脉冲响应中的)奇异值分布;22:(现实的脉冲响应中的)奇异值分布;31:递归系统矩阵推测部;32:系统特性推测部;33:系统维数推测部;34:系统稳定性评价部;41:误差平方和范数分布;42:误差平方和范数阈值;50:基于冲击振动的系统辨识装置;51:系统输入施加时刻确定部;52:脉冲响应变换部;60:基于频率响应的系统辨识装置;61:逆傅立叶变换部;71:DC伺服马达;72:频率响应变换部;73:控制系统或者滤波器参数设计部。
具体实施方式
以下,根据附图,详细说明本发明的系统辨识装置的实施方式。另外,本发明不限于该实施方式。
实施方式1.
图1是示出本发明的系统辨识装置的实施方式1的整体结构的框图。图2是示出在本实施方式的系统辨识装置中,系统输入输出的时间波形的概略图。在图1所示的系统辨识装置10中,将作为辨识对象的动态系统的脉冲响应g(jTs)(j=0,1,2,…)作为输入。如图2所示,该脉冲响应g(jTs)(j=0,1,2,…)通过用理想的脉冲输入构成了向动态系统的系统输入11的情况下的时刻0以后的系统输出12给出,该理想的脉冲输入在时刻0振幅为1且在该时刻以外振幅为0。
在系统辨识装置10中,根据这样得到的动态系统的脉冲响应g(jTs)(j=0,1,2,…),由直通项辨识部1辨识记述动态系统的线性离散时间系统的直通项Dd,进而由块汉克尔矩阵生成部2生成块汉克尔矩阵Hkl
在奇异值分解部3中,对从块汉克尔矩阵生成部2输出的块汉克尔矩阵Hkl进行奇异值分解,输出将块汉克尔矩阵Hkl的左奇异矢量作为列矢量的第1正交矩阵U、将块汉克尔矩阵Hkl的右奇异矢量作为列矢量的第2正交矩阵V以及块汉克尔矩阵Hkl的奇异值σi(i=1,2,3,…)。
在系统维数决定部4中,根据从奇异值分解部3输出的第1正交矩阵U、第2正交矩阵V及奇异值σi(i=1,2,3,…)、和作业者所指定的系统维数的搜索范围ni=(n1,n2,…,na)(其中n1<n2<…<na),针对属于该搜索范围的各个维数ni(i=1,2,…,a),辨识记述动态系统的线性离散时间系统的除了直通项Dd以外的系统矩阵。进而,根据该系统矩阵以及从直通项辨识部1输出的直通项Dd,针对属于该搜索范围的各个维数ni(i=1,2,…,a),计算线性离散时间系统的推测脉冲响应,根据与动态系统的实际的脉冲响应(图1中的动态系统的系统特性)的比较,决定系统维数n。
在系统矩阵辨识部5中,根据从奇异值分解部3输出的第1正交矩阵U、第2正交矩阵V及奇异值σi(i=1,2,3,…)、和从系统维数决定部4输出的系统维数n,辨识线性离散时间系统的除了直通项Dd以外的系统矩阵Ad、Bd、Cd
在系统辨识装置10中,将由直通项辨识部1辨识出的直通项Dd以及由系统矩阵辨识部5辨识出的系统矩阵Ad、Bd、Cd作为记述动态系统的线性离散时间系统而最终地输出。
图3是示出在本实施方式的系统辨识装置10中,块汉克尔矩阵Hkl的奇异值σi和维数(i=1,2,3,…)的关系的概略图。图4是示出在本实施方式的系统辨识装置10中,系统维数决定部4的内部结构和动作的框图。图5是示出在本实施方式的系统辨识装置10中,辨识出的线性离散时间系统的推测脉冲响应和动态系统的实际的脉冲响应在时域中的误差平方和范数||eni||与维数ni(i=1,2,…,a)的关系的概略图。另外,在图4中,“i=1”的处理、“i++”的处理以及“i≤a”的判断既可以由系统特性推测部32或者系统维数推测部33承担,也可以由未图示的其他结构承担。
如图3所示,根据动态系统的脉冲响应g(jTs)(j=0,1,2,…)而生成的块汉克尔矩阵Hkl的奇异值σi在理想的脉冲响应中,相对维数(i=1,2,3,…)成为例如奇异值分布21所示的关系。
在该情况下,能够明确地规定具有有效值的奇异值的个数,该个数与动态系统的系统维数n对应(在图3的情况下系统维数n=4)。另一方面,根据受到观测噪声等的影响的现实的脉冲响应计算出的奇异值σi相对维数(i=1,2,3,…)成为例如奇异值分布22所示的关系,所以具有有效值的奇异值和成为可忽略的微小的值的奇异值的边界变得不明确,存在未必始终决定出最佳的系统维数n、或者关于系统维数n的决定需要反复试验这样的问题。
因此,如图4所示,系统维数决定部4通过递归系统矩阵推测部31,关于属于作业者所指定的系统维数的搜索范围ni=(n1,n2,…,na)(其中n1<n2<…<na)的第1维数ni,使用与比第1维数ni低1个等级的第2维数ni-1对应的系统矩阵Ad,ni-1、Bd,ni-1、Cd,ni-1的辨识结果、和从奇异值分解部3输出的第1正交矩阵U、第2正交矩阵V及奇异值σi(i=1,2,3,…)中分别比第2维数ni-1大且在第1维数ni以下的左奇异矢量uj、右奇异矢量vj及奇异值σj(j=ni-1+1,ni-1+2,…,ni),通过递归方法来辨识与第1维数ni对应的系统矩阵Ad,ni、Bd,ni、Cd,ni
在系统辨识装置10中,系统维数决定部4具备递归系统矩阵推测部31,该递归系统矩阵推测部31使用与在搜索范围ni=(n1,n2,…,na)(其中n1<n2<…<na)中比第1维数ni低1个等级的第2维数ni-1对应的系统矩阵Ad,ni-1、Bd,ni-1、Cd,ni-1的辨识结果、和从第1正交矩阵U、第2正交矩阵V及奇异值σi(i=1,2,3,…)中分别比第2维数ni-1大且在第1维数ni以下的左奇异矢量uj、右奇异矢量vj及奇异值σj(j=ni-1+1,ni-1+2,…,ni),通过递归方法来辨识与第1维数ni对应的系统矩阵Ad,ni、Bd,ni、Cd,ni,该系统辨识装置10能够降低用于决定相对现实的动态系统一致度高的系统维数的运算量。
接下来,系统特性推测部32针对属于系统维数的搜索范围ni=(n1,n2,…,na)(其中n1<n2<…<na)的各个维数,根据从递归系统矩阵推测部31输出的系统矩阵Ad,ni、Bd,ni、Cd,ni以及从直通项辨识部1输出的直通项Dd,计算辨识出的线性离散时间系统的推测脉冲响应。然后,对i加上1,在i≤a的情况下,由递归系统矩阵推测部31辨识系统矩阵Ad,ni、Bd,ni、Cd,ni,在i>a的情况下进入到系统维数推测部33的处理。
系统维数推测部33计算从系统特性推测部32输出的线性离散时间系统的推测脉冲响应和动态系统的实际的脉冲响应(图4中的动态系统的系统特性)在时域中的误差平方和eni(i=1,2,…,a),将如图5所示误差平方和范数分布41在误差平方和范数阈值42以下的维数中的最小的维数决定为系统维数n而输出(在图5的情况下系统维数n=n6)。
即,在系统辨识装置10中,系统维数决定部4具备系统特性推测部32和系统维数推测部33,该系统特性推测部32针对属于搜索范围ni=(n1,n2,…,na)的各维数,计算线性离散时间系统的系统特性作为推测脉冲响应并输出,该系统维数推测部33将从系统特性推测部32输出的线性离散时间系统的推测脉冲响应和动态系统的实际的脉冲响应在时域中的误差平方和范数在阈值以下的维数中的最小的维数决定为系统维数并输出,在该系统辨识装置10中,即使在根据现实的脉冲响应计算出的块汉克尔矩阵的奇异值为平滑的单调减少,具有有效值的奇异值和成为可忽略的微小的值的奇异值的边界变得不明确的情况下,也能够从系统维数的决定排除反复试验,进行相对现实的动态系统在时域中一致度高的系统维数的决定、和记述动态系统的线性离散时间系统的辨识。
接下来,说明动作。设为能够通过下述的式(1)所示的1输入P输出的n维线性离散时间系统来记述作为辨识对象的动态系统。如果通过下述的式(2)构成向上述动态系统的系统输入u(jTs),以成为图2所示的系统输入11、即理想的脉冲输入,则与下述的式(1)对应的系统输出y(jTs)、即图2所示的系统输出12成为下述的式(3),通过时刻0以后的系统输出12给出动态系统的脉冲响应g(jTs)(j=0,1,2,…)。
[式1]
x ( ( j + 1 ) T s ) = A d x ( jT s ) + B d u ( jT s ) y ( jT s ) = C d x ( jT s ) + D d u ( jT s ) ... ( 1 )
其中,状态矢量x∈Rn,系统输入u∈R,系统输出y∈RP,系统矩阵Ad∈Rn×n、Bd∈Rn、Cd∈RP×n、Dd∈RP,脉冲响应列g(jTs)(j=0,1,2,…)。
[式2]
u ( jT s ) = 1 ( j = 0 ) 0 ( j &NotEqual; 0 ) ... ( 2 )
[式3]
y ( 0 ) = g ( 0 ) = D d y ( T s ) = g ( T s ) = C d B d y ( 2 T s ) = g ( 2 T s ) = C d A d B d y ( 3 T s ) = g ( 3 T s ) = C d A d 2 B d . . . ... ( 3 )
如图1所示,本实施方式的系统辨识装置10将以上得到的动态系统的脉冲响应g(jTs)(j=0,1,2,…)作为输入,在直通项辨识部1中用下述的式(4)辨识用上述式(1)记述动态系统的线性离散时间系统的直通项Dd,进而通过块汉克尔矩阵生成部2生成用下述的式(5)给出的块汉克尔矩阵Hkl
[式4]
Dd=g(0) …(4)
[式5]
接下来,通过奇异值分解部3,对从块汉克尔矩阵生成部2输出的块汉克尔矩阵Hkl进行奇异值分解,输出将用下述的式(6)给出的块汉克尔矩阵Hkl的左奇异矢量uj作为列矢量的第1正交矩阵U、将块汉克尔矩阵Hkl的右奇异矢量vj作为列矢量的第2正交矩阵V以及块汉克尔矩阵Hkl的奇异值σi(i=1,2,3,…)。
[式6]
Hkl=U∑VT …(6)
其中,第1正交矩阵U=[u1u2…ukP]∈RkP×kP,第2正交矩阵V=[v1v2…ul]∈Rl×l,块汉克尔矩阵Hkl的奇异值σ1≥σ2≥…≥σn≥σn+1≥σn+2≥…,Σ用下述的式(7)表示。
[式7]
能够根据下述的式(8)的关系决定作为对象的动态系统的系统维数n,在下述的式(8)中,在块汉克尔矩阵Hkl的奇异值中具有有效值的奇异值是n个,与这些奇异值相比在第n+1个以后充分小,而成为可忽略的微小的值。
[式8]
σ1≥σ2≥…≥σn>>σn+1≥σn+2≥… …(8)
如图3所示,根据动态系统的脉冲响应g(jTs)(j=0,1,2,…)生成的块汉克尔矩阵Hkl的奇异值σi在理想的脉冲响应中相对维数(i=1,2,3,…)成为例如奇异值分布21所示的关系。在该情况下,能够明确地规定具有有效值的奇异值的个数,能够根据该个数决定动态系统的系统维数n(在图3的情况下系统维数n=4)。
另一方面,根据受到观测噪声等的影响的现实的脉冲响应计算出的奇异值σi相对维数(i=1,2,3,…)成为例如奇异值分布22所示的关系,所以具有有效值的奇异值和成为可忽略的微小的值的奇异值的边界变得不明确。因此,在以往的方法中,存在未必始终决定出最佳的系统维数n、或者关于系统维数n的决定需要反复试验这样的课题。因此,在本实施方式的系统辨识装置10中,当“在时域中最适合于实际的脉冲响应”这样的前提下,在系统维数决定部4中决定最佳的系统维数n。
系统维数决定部4如图1所示,根据从奇异值分解部3输出的第1正交矩阵U、第2正交矩阵V及奇异值σi(i=1,2,3,…)、和作业者所指定的系统维数的搜索范围ni=(n1,n2,…,na)(其中n1<n2<…<na),针对属于该搜索范围的各个维数ni(i=1,2,…,a),辨识记述动态系统的线性离散时间系统的除了直通项Dd以外的系统矩阵。
进而,根据该系统矩阵以及从直通项辨识部1输出的直通项Dd,计算针对属于该搜索范围的各个维数ni(i=1,2,…,a)辨识出的线性离散时间系统的推测脉冲响应,根据与动态系统的实际的脉冲响应(图1中的动态系统的系统特性)的比较,决定系统维数n。具体而言,如图4所示,通过递归系统矩阵推测部31,关于属于作业者所指定的系统维数的搜索范围ni=(n1,n2,…,na)(其中n1<n2<…<na)的第1维数ni,使用与比第1维数ni低1个等级的第2维数ni-1对应的系统矩阵Ad,ni-1、Bd,ni-1、Cd,ni-1的辨识结果、和从奇异值分解部3输出的第1正交矩阵U、第2正交矩阵V及奇异值σi(i=1,2,3,…)中分别比第2维数ni-1大且在第1维数ni以下的左奇异矢量uj、右奇异矢量vj以及奇异值σj(j=ni-1+1,ni-1+2,…,ni),通过下式所示的递归方法来辨识与第1维数ni对应的系统矩阵Ad,ni、Bd,ni、Cd,ni。首先,通过下述的式(9)表示与第1维数ni对应的扩展可观测矩阵。
[式9]
O k , n i = U n i &Sigma; n i 1 / 2 = O k , n i - 1 &sigma; n i - 1 + 1 u n i - 1 + 1 ... &sigma; n i u n i &Element; R k P &times; n i = U ( : , 1 : n i ) &Sigma; ( 1 : n i , 1 : n i ) 1 / 2 ( n i = n 1 ) O k , n i - 1 U ( : , n i - 1 + 1 : n i ) &Sigma; ( n i - 1 + 1 : n i , n i - 1 + 1 : n i ) 1 / 2 ( n i > n 1 ) ... ( 9 )
通过下述的式(10)表示与第1维数ni对应的扩展可到达矩阵。
[式10]
C l , n i = &Sigma; n i 1 / 2 V n i T = C l , n i - 1 &sigma; n i - 1 + 1 v n i - 1 + 1 T . . . &sigma; n i v n i T &Element; R n i &times; l = &Sigma; ( 1 : n i , 1 : n i ) 1 / 2 V ( : , 1 : n i ) T ( n i = n 1 ) C l , n i - 1 &Sigma; ( n i - 1 + 1 : n i , n i - 1 + 1 : n i ) 1 / 2 V ( : , n i - 1 + 1 : n i ) T ( n i > n 1 ) ... ( 10 )
然后,通过下述的式(11)表示与第1维数ni对应的系统矩阵Ad,ni、Bd,ni、Cd,ni
[式11]
接下来,在系统特性推测部32中,针对属于系统维数的搜索范围ni=(n1,n2,…,na)(其中n1<n2<…<na)的各个维数ni,根据从递归系统矩阵推测部31输出的系统矩阵Ad,ni、Bd,ni、Cd,ni以及从直通项辨识部1输出的直通项Dd,应用上述式(2)作为针对上述式(1)的系统输入u(jTs),从而计算辨识出的线性离散时间系统的推测脉冲响应g^ni(jTs)(j=0,1,2,…)。
另外,在本说明书中,“g^”是在“g”之上配置了“^(抑扬符(circumflex))”的字符的代替标记。
系统维数推测部33用下述的式(12)计算从系统特性推测部32输出的线性离散时间系统的推测脉冲响应g^ni(jTs)(j=0,1,2,…)和动态系统的实际的脉冲响应g(jTs)(j=0,1,2,…)(图4中的动态系统的系统特性)在时域中的误差平方和eni(i=1,2,…,a)。
[式12]
e n i = &Sigma; j = 0 k + l - 1 ( g ( jT s ) - g ^ n i ( jT s ) ) 2 &Element; R P ... ( 12 )
使该误差平方和范数||eni||成为最小的维数ni为“在时域中最适合于实际的脉冲响应的”系统维数n。另一方面,如果观测噪声是白噪声,则无论其噪声等级如何,伴随维数ni的增加,实际的范数||eni||都单调减少,如图5所示,在预定的维数(在图5的情况下为n6)以上成为大致恒定值。因此,此处,为了避免系统维数n的推测值超出必要地成为高维数,规定用下述的式(13)给出的误差平方和范数阈值42,将误差平方和范数分布41在误差平方和范数阈值42以下的维数中的最小的维数决定为系统维数n而输出(在图5的情况下系统维数n=n6)。
[式13]
最后,在系统矩阵辨识部5中,根据从奇异值分解部3输出的第1正交矩阵U、第2正交矩阵V及奇异值σi(i=1,2,3,…)、和从系统维数决定部4输出的系统维数n,通过下述的式(14)来辨识线性离散时间系统的除了直通项Dd以外的系统矩阵Ad、Bd、Cd
[式14]
系统辨识装置10将由系统矩阵辨识部5辨识出的系统矩阵Ad、Bd、Cd以及由直通项辨识部1辨识出的直通项Dd作为记述动态系统的线性离散时间系统最终地输出。
这样,通过本实施方式的系统辨识装置10,即使在根据现实的脉冲响应计算出的块汉克尔矩阵Hkl的奇异值σi(i=1,2,3,…)为平滑的单调减少,具有有效值的奇异值和成为在辨识中可忽略的微小的值的奇异值的边界变得不明确的情况下,也能够从系统维数n的决定排除反复试验,进行相对现实的动态系统在时域中一致度高的系统维数n的决定以及记述动态系统的线性离散时间系统的辨识。
另外,通过递归系统矩阵推测部31,能够降低用于决定相对现实的动态系统一致度高的系统维数n的运算量。
另外,在本实施方式的系统辨识装置10中,计算线性离散时间系统的系统特性而作为推测脉冲响应,将该脉冲响应和动态系统的实际的脉冲响应在时域中的误差平方和范数分布41在误差平方和范数阈值42以下的维数中的最小的维数决定为系统维数n。但是,本发明不限于此,也可以计算线性离散时间系统的系统特性而作为推测频率响应,根据该频率响应和对动态系统的脉冲响应进行傅立叶变换而得到的实际的频率响应在频域中的误差平方和,决定系统维数n。在该情况下,也可以进而根据动态系统的实际的频率响应来决定权重函数,根据对线性离散时间系统的推测频率响应和动态系统的实际的频率响应在频域中的误差平方值乘以该权重函数而得到的值的相加值,决定系统维数n。
实施方式2.
在本实施方式中,说明还能够应用通过动态系统的冲击振动得到的现实的系统输入输出,并且在现实的动态系统是稳定的情形明确的情况下能够限定于稳定系统而进行辨识的系统辨识装置。
作为用于得到动态系统的脉冲响应的方法,最一般的方法是通过脉冲锤对作为对象的动态系统进行冲击振动并测量此时的系统输出的方法。但是,通过冲击振动给出的现实的系统输入为正弦半波状,所以与真正的脉冲输入不同。因此,通过该输入得到的系统输出也与动态系统的真正的脉冲响应不同,所以存在无法将该系统输出直接应用于Ho-Kalman的方法这样的问题。
另外,在以往的基于脉冲响应的系统辨识装置中,作为辨识结果得到的线性离散时间系统(Ad,Bd,Cd,Dd)的稳定性完全未被考虑,所以还存在尽管现实的动态系统是稳定的,有时还是被辨识为不稳定系统这样的问题。
在本实施方式中,鉴于上述问题,其目的在于得到一种系统辨识装置,该系统辨识装置还能够应用通过动态系统的冲击振动得到的现实的系统输入输出,并且在现实的动态系统是稳定的情形明确的情况下,能够限定于稳定系统而进行辨识。
在本实施方式中,示出系统辨识装置的整体结构的框图与实施方式1的图1相同,示出块汉克尔矩阵Hkl的奇异值σi和维数(i=1,2,3,…)的关系的概略图与实施方式1的图3相同,示出辨识出的线性离散时间系统的推测系统输出和动态系统的实际的系统输出在时域中的误差平方和范数||eni||与维数ni(i=1,2,…,a)的关系的概略图与实施方式1的图5相同。
图6是示出在本实施方式的系统辨识装置中,对动态系统进行了冲击振动的情况下的系统输入输出的时间波形的概略图。图7是示出在本实施方式的系统辨识装置中,将对动态系统进行了冲击振动的情况下的系统输入输出变换为该系统的脉冲响应的步骤的框图。此时,作为图7的脉冲响应变换部52的输出的动态系统的脉冲响应成为图1的直通项辨识部1以及块汉克尔矩阵生成部2的输入。将对作为辨识对象的动态系统进行了冲击振动的情况下的系统输入11以及系统输出12作为输入,将图7所示的系统输入施加时刻确定部51及脉冲响应变换部52以及图1的系统辨识装置10整体一并地作为基于冲击振动的系统辨识装置50。
如图6所示,在本实施方式的基于冲击振动的系统辨识装置50中,根据对作为辨识对象的动态系统进行了冲击振动的情况下的系统输入11以及系统输出12,辨识记述该系统的线性离散时间系统。关于该系统输入11以及系统输出12,如图7所示,通过系统输入施加时刻确定部51来确定系统输入11具有有效值的多个时刻。接下来,通过脉冲响应变换部52,将与从系统输入施加时刻确定部51输出的多个时刻对应的系统输入的相加值作为脉冲输入振幅,将从系统输入施加时刻确定部51输出的多个时刻的最大值作为脉冲输入施加时刻,将脉冲输入施加时刻以后的系统输出除以脉冲输入振幅而得到的信号作为动态系统的脉冲响应而输出。
与脉冲响应对应的信息是对动态系统进行了冲击振动时的系统输入11以及系统输出12,基于冲击振动的系统辨识装置50具备系统输入施加时刻确定部51以及脉冲响应变换部52,该系统输入施加时刻确定部51确定系统输入11具有有效值的多个时刻,该脉冲响应变换部52将与从系统输入施加时刻确定部51输出的多个时刻对应的系统输入的相加值作为脉冲输入振幅,将从系统输入施加时刻确定部51输出的多个时刻的最大值作为脉冲输入施加时刻,将脉冲输入施加时刻以后的系统输出除以所述脉冲输入振幅而得到的信号作为所述脉冲响应输出,该基于冲击振动的系统辨识装置50将对作为辨识对象的动态系统进行了冲击振动的情况下的系统输入11及系统输出12以及搜索范围ni=(n1,n2,…,na)(其中n1<n2<…<na)作为输入来辨识线性离散时间系统,能够根据对动态系统进行冲击振动而得到的现实的系统输入输出来辨识记述该系统的线性离散时间系统。
图8是示出在本实施方式的系统辨识装置中,系统维数决定部4a的内部结构和动作的框图。在图8中,附加与图4相同的符号的构成要素是与实施方式1相同或者等同的构成要素。如图8所示,在本实施方式的系统维数决定部4a中,针对属于作业者所指定的系统维数的搜索范围ni=(n1,n2,…,na)(其中n1<n2<…<na)的各个维数ni,根据由递归系统矩阵推测部31辨识出的系统矩阵Ad,ni,系统稳定性评价部34评价线性离散时间系统的稳定性。另外,在图8中,“i=1”的处理、“i++”的处理以及“i≤a”的判断既可以由系统特性推测部32或者系统维数推测部33承担,也可以由未图示的其他结构承担。
系统特性推测部32针对在系统稳定性评价部34中判断为稳定系统的维数,根据从递归系统矩阵推测部31输出的系统矩阵Ad,ni、Bdni、Cd,ni以及从直通项辨识部1输出的直通项Dd,关于辨识出的线性离散时间系统,计算与现实的系统输入11对应的推测系统输出。
在系统辨识装置10中,系统维数决定部4a具备针对属于搜索范围ni=(n1,n2,…,na)(其中n1<n2<…<na)的各维数评价线性离散时间系统的稳定性的系统稳定性评价部34,系统辨识装置10根据与成为稳定系统的维数对应的线性离散时间系统的系统特性来决定系统维数,在现实的动态系统是稳定系统的情形明确的情况下,能够限定于稳定系统而辨识记述系统的线性离散时间系统。
系统维数推测部33计算从系统特性推测部32输出的线性离散时间系统的推测系统输出和动态系统的实际的系统输出12(图8中的动态系统的系统特性)在时域中的误差平方和eni(ni:成为稳定系统的维数),将如图5所示误差平方和范数分布41在误差平方和范数阈值42以下的维数中的最小的维数决定为系统维数n而输出(在图5的情况下系统维数n=n6)。
接下来,说明动作。设为能够通过上述式(1)将作为辨识对象的动态系统记述为1输入P输出的n维数线性离散时间系统。对该动态系统进行了冲击振动的情况下的系统输入u(jTs)成为例如图6所示的系统输入11那样的时间波形,用下述的式(15)给出该情况下的系统输入u(jTs)。
[式15]
u ( jT s ) = { a j ( - 2 &le; j &le; 0 ) 0 ( j < - 2 , 0 < j ) ... ( 15 )
此时,与上述式(1)对应的系统输出y(jTs)、即图6所示的系统输出12成为下述的式(16),所以难以根据系统输出y(jTs)直接地求出脉冲响应g(jTs)(j=0,1,2,…)。
[式16]
y ( 0 ) = g ( 0 ) a 0 + g ( T s ) a - 1 + g ( 2 T s ) a - 2 = D d a 0 + C d B d a - 1 + C d A d B d a - 2 y ( T s ) = g ( T s ) a 0 + g ( 2 T s ) a - 1 + g ( 3 T s ) a - 2 = C d B d a 0 + C d A d B d a - 1 + C d A d 2 B d a - 2 y ( 2 T s ) = g ( 2 T s ) a 0 + g ( 3 T s ) a - 1 + g ( 4 T s ) a - 2 = C d A d B d a 0 + C d A d 2 B d a - 1 + C d A d 2 B d a - 2 y ( 3 T s ) = g ( 3 T s ) a 0 + g ( 4 T s ) a - 1 + g ( 5 T s ) a - 2 = C d A d 2 B d a 0 + C d A d 3 B d a - 1 + C d A d 4 B d a - 2 . . . ... ( 16 )
因此,在本实施方式的基于冲击振动的系统辨识装置50中,在冲击振动中的加振时间很短这样的前提下,通过图2那样的理想的脉冲输入来近似图6所示的现实的系统输入11。具体而言,通过图7所示的系统输入施加时刻确定部51,将根据叠加于系统输入u(jTs)的噪声等级而规定的比例阈值和系统输入u(jTs)的最大值相乘而得到的下述的式(17)作为系统输入阈值13,将系统输入u(jTs)在系统输入阈值13以上的时刻确定为系统输入u(jTs)具有有效值的多个时刻j’Ts(在图6的情况下j’=-2,-1,0)。
[式17]
系统输入比例阈值·max(u(jTs)) …(17)
在基于冲击振动的系统辨识装置50中,系统输入施加时刻确定部51将比例阈值和系统输入的最大值相乘而得到的值作为系统输入阈值13,将系统输入成为系统输入阈值13以上的时刻确定为系统输入具有有效值的多个时刻,该基于冲击振动的系统辨识装置50能够从对动态系统进行冲击振动而得到的现实的系统输入,正确地提取系统输入具有有效值的加振时刻。
接下来,通过脉冲响应变换部52,与从系统输入施加时刻确定部51输出的多个时刻j’Ts对应地,通过下述的式(18)将系统输入u(jTs)的相加值作为脉冲输入振幅a导出,通过下述的式(19)将从系统输入施加时刻确定部51输出的多个时刻j’Ts的最大值作为脉冲输入施加时刻j”Ts(在图6的情况下j”=0)导出,通过下述的式(20)计算动态系统的脉冲响应g(jTs)(j=0,1,2,…),而作为将脉冲输入施加时刻j”Ts以后的系统输出y(jTs)除以脉冲输入振幅a而得到的信号。
另外,在本说明书中,“a”是在“a”之上配置了“(波浪符)”的字符的代替标记。
[式18]
a ~ = &Sigma; j = j &prime; u ( jT s ) ... ( 18 )
[式19]
j″Ts=max(j′Ts) …(19)
[式20]
g ( jT s ) = 1 a ~ y ( ( j &prime; &prime; + j ) T s ) , ( j = 0 , 1 , 2 , ... ) ... ( 20 )
在本实施方式中,根据对动态系统进行了冲击振动的情况下的系统输入输出,将用上述式(20)得到的动态系统的脉冲响应g(jTs)(j=0,1,2,…)作为向图1所示的系统辨识装置10的输入。在本实施方式中,如图1所示,在系统辨识装置10中,将用上述式(20)得到的动态系统的脉冲响应g(jTs)(j=0,1,2,…)作为输入,在直通项辨识部1中用上述式(4)辨识记述动态系统的上述式(1)的线性离散时间系统的直通项Dd,并且通过块汉克尔矩阵生成部2生成用下述的式(21)给出的块汉克尔矩阵Hkl。接下来,通过奇异值分解部3,对从块汉克尔矩阵生成部2输出的块汉克尔矩阵Hkl进行奇异值分解,输出用上述式(6)给出的块汉克尔矩阵Hkl的第1正交矩阵U、第2正交矩阵V以及奇异值σi(i=1,2,3,…)。
[式21]
系统维数决定部4a如图8所示通过递归系统矩阵推测部31,关于属于作业者所指定的系统维数的搜索范围ni=(n1,n2,…,na)(其中n1<n2<…<na)的各个维数ni,通过上述式(9)~(11)所示的递归方法来辨识对应的系统矩阵Ad,ni、Bd,ni、Cd,ni
系统稳定性评价部34针对属于作业者所指定的系统维数的搜索范围ni=(n1,n2,…,na)(其中n1<n2<…<na)的各个维数ni,根据由递归系统矩阵推测部31辨识出的系统矩阵Ad,ni,通过下述的式(22)来评价线性离散时间系统的稳定性。
[式22]
接下来,在系统特性推测部32中,针对在系统稳定性评价部34中判断为稳定系统的维数,根据从递归系统矩阵推测部31输出的系统矩阵Ad,ni、Bd,ni、Cd,ni以及从直通项辨识部1输出的直通项Dd,计算对上述式(1)应用了现实的系统输入u(jTs)(例如式(15))的情况下的推测系统输出y^ni(jTs)。
另外,在本说明书中,“y^”是在“y”之上配置了“^(抑扬符)”的字符的代替标记。
系统维数推测部33用下述的式(23)计算从系统特性推测部32输出的线性离散时间系统的推测系统输出y^ni(jTs)和动态系统的实际的系统输出y(jTs)(图8中的动态系统的系统特性)在时域中的误差平方和eni(ni:成为稳定系统的维数)。
[式23]
e n i = &Sigma; j - min ( j &prime; ) min ( j &prime; ) + k + l - 1 ( y ( jT s ) - y ^ n i ( jT s ) ) 2 &Element; R P ... ( 23 )
使该误差平方和范数||eni||成为最小的维数ni为“在时域中最适合于实际的系统输出的”稳定的系统维数n。此处,将如图5所示误差平方和范数分布41在用上述式(13)给出的误差平方和范数阈值42以下的维数中的最小的维数决定为系统维数n而输出(在图5的情况下系统维数n=n6)。
最后,在系统矩阵辨识部5中,根据从奇异值分解部3输出的第1正交矩阵U、第2正交矩阵V及奇异值σi(i=1,2,3,…)、和从系统维数决定部4a输出的系统维数n,用上述式(14)辨识线性离散时间系统的除了直通项Dd以外的系统矩阵Ad、Bd、Cd,将由直通项辨识部1辨识出的直通项Dd以及由系统矩阵辨识部5辨识出的系统矩阵Ad、Bd、Cd作为记述动态系统的线性离散时间系统输出。
这样,通过本实施方式的基于冲击振动的系统辨识装置50,将现实的系统输入输出变换为脉冲响应,即使在根据该脉冲响应计算出的块汉克尔矩阵Hkl的奇异值σi(i=1,2,3,…)为平滑的单调减少,具有有效值的奇异值和成为在辨识中可忽略的微小的值的奇异值的边界变得不明确的情况下,也能够从系统维数n的决定排除反复试验,进行相对现实的动态系统在时域中一致度高的系统维数n的决定以及记述动态系统的线性离散时间系统的辨识。
另外,通过递归系统矩阵推测部31,能够降低用于决定相对现实的动态系统一致度高的系统维数n的运算量。
进而,从对动态系统进行冲击振动而得到的现实的系统输入输出正确地提取系统输入具有有效值的加振时刻,并据此变换为动态系统的脉冲响应,从而能够从对动态系统进行了冲击振动的情况下的现实的系统输入输出,正确地辨识记述该系统的线性离散时间系统。
另外,通过系统稳定性评价部34,在现实的动态系统是稳定系统的情形明确的情况下,能够进行限定于稳定系统的线性离散时间系统的辨识。
另外,在本实施方式的基于冲击振动的系统辨识装置50中,计算线性离散时间系统的系统特性而作为推测系统输出,将该系统输出和动态系统的实际的系统输出在时域中的误差平方和范数分布41在误差平方和范数阈值42以下的维数中的最小的维数决定为系统维数n。但是,本发明不限于此,也可以计算线性离散时间系统的系统特性而作为推测频率响应,根据该频率响应和实际的频率响应在频域中的误差平方和,决定系统维数n,该实际的频率响应是针对变换动态系统的系统输入输出而得到的脉冲响应对该脉冲响应进行傅立叶变换而得到的。
在该情况下,也可以进而根据动态系统的实际的频率响应来决定权重函数,根据对线性离散时间系统的推测频率响应和动态系统的实际的频率响应在频域中的误差平方值乘以该权重函数而得到的值的相加值,决定系统维数n。
根据本实施方式,能够得到一种系统辨识装置,该系统辨识装置还能够应用通过动态系统的冲击振动得到的现实的系统输入输出,并且在现实的动态系统是稳定的情形明确的情况下,能够限定于稳定系统而进行辨识。
实施方式3.
在本实施方式中,说明与脉冲响应对应的信息是动态系统的频率响应的情况。在本实施方式中,示出系统辨识装置的整体结构的框图与实施方式1的图1相同,示出块汉克尔矩阵Hkl的奇异值σi和维数(i=1,2,3,…)的关系的概略图与实施方式1的图3相同,示出辨识出的线性离散时间系统的推测频率响应和动态系统的实际的频率响应在频域中的误差平方和范数||eni||与维数ni(i=1,2,…,a)的关系的概略图与实施方式1的图5相同。另外,在本实施方式中,系统维数决定部与实施方式2的图8所示的系统维数决定部4a相同。
图9是示出在本实施方式的系统辨识装置中,将动态系统的频率响应变换为该系统的脉冲响应的步骤的框图。此时,作为图9的逆傅立叶变换部61的输出的动态系统的脉冲响应成为图1的直通项辨识部1以及块汉克尔矩阵生成部2的输入。将动态系统的频率响应作为输入,将包括图9所示的逆傅立叶变换部61以及图1的系统辨识装置10的系统辨识装置作为基于频率响应的系统辨识装置60。在本实施方式的基于频率响应的系统辨识装置60中,根据作为辨识对象的动态系统的频率响应,辨识记述该系统的线性离散时间系统。如图9所示,通过逆傅立叶变换部61对该动态系统的频率响应进行逆傅立叶变换,并作为动态系统的脉冲响应输出。
与脉冲响应对应的信息是动态系统的频率响应,该基于频率响应的系统辨识装置60具备通过频率响应的逆傅立叶变换输出对应的脉冲响应的逆傅立叶变换部61,将作为辨识对象的动态系统的频率响应以及搜索范围ni=(n1,n2,…,na)(其中n1<n2<…<na)作为输入,辨识线性离散时间系统,该基于频率响应的系统辨识装置60能够根据动态系统的频率响应,辨识记述该系统的线性离散时间系统。
接下来,说明动作。设为能够用上述式(1)将作为辨识对象的动态系统记述为1输入P输出的n维数线性离散时间系统。在给出了该动态系统的(复数)频率响应H(kΔf)(k=0,1,2,…,N-1)的情况下,通过图9所示的逆傅立叶变换部61,用下述的式(24)进行逆傅立叶变换,得到动态系统的脉冲响应g(jTs)(j=0,1,2,…)。
[式24]
g ( jT s ) = 1 N &Sigma; k = 0 N - 1 H ( k &Delta; f ) e x p &lsqb; i &CenterDot; 2 &pi; j k N &rsqb; , ( j = 0 , 1 , 2 , ... , N - 1 ) ... ( 24 )
其中,采样周期Ts=T/N,采样频率fs=1/Ts=N/T,频率分辨率Δf=1/T,时刻t=jTs=jT/N,频率f=kΔf=k/T。
在本实施方式中,根据动态系统的频率响应H(kΔf)(k=0,1,2,…,N-1),将用上述式(24)得到的动态系统的脉冲响应g(jTs)(j=0,1,2,…)作为向系统辨识装置10的输入。在本实施方式中,如图1所示,在系统辨识装置10中,将用上述式(24)得到的动态系统的脉冲响应g(jTs)(j=0,1,2,…)作为输入,在直通项辨识部1中用式(4)辨识记述动态系统的上述式(1)的线性离散时间系统的直通项Dd,并且通过块汉克尔矩阵生成部2生成用下述的式(25)给出的块汉克尔矩阵Hkl
[式25]
接下来,通过奇异值分解部3,对从块汉克尔矩阵生成部2输出的块汉克尔矩阵Hkl进行奇异值分解,输出用上述式(6)给出的块汉克尔矩阵Hkl的第1正交矩阵U、第2正交矩阵V以及奇异值σi(i=1,2,3,…)。
系统维数决定部4a如图8所示,通过递归系统矩阵推测部31,关于属于作业者所指定的系统维数的搜索范围ni=(n1,n2,…,na)(其中n1<n2<…<na)的各个维数ni,通过上述式(9)~(11)所示的递归方法来辨识对应的系统矩阵Ad,ni、Bd,ni、Cd,ni
系统稳定性评价部34针对属于作业者所指定的系统维数的搜索范围ni=(n1,n2,…,na)(其中n1<n2<…<na)的各个维数ni,根据由递归系统矩阵推测部31辨识出的系统矩阵Ad,ni,用上述式(22)评价线性离散时间系统的稳定性。
接下来,在系统特性推测部32中,针对在系统稳定性评价部34中判断为稳定系统的维数,根据从递归系统矩阵推测部31输出的系统矩阵Ad,ni、Bd,ni、Cd,ni以及从直通项辨识部1输出的直通项Dd,计算线性离散时间系统的推测频率响应H^ni(kΔf)(k=0,1,2,…,N-1)。
另外,在本说明书中,“H^ni”是在“Hni”之上配置了“^(抑扬符)”的字符的代替标记。
在本实施方式的基于频率响应的系统辨识装置60中,在“在频域中最适合于实际的频率响应”这样的前提下,在系统维数推测部33中决定最佳的系统维数n。具体而言,根据动态系统的实际的频率响应H(kΔf)(k=0,1,2,…,N-1),决定例如对高增益且低频区域附加了权重的下述的式(26)那样的权重函数W(kΔf)(k=0,1,2,…,N-1),用下述的式(27)计算对从系统特性推测部32输出的线性离散时间系统的推测频率响应H^ni(kΔf)和动态系统的实际的频率响应H(kΔf)(图8中的动态系统的系统特性)在频域中的误差平方值乘以权重函数W(kΔf)而得到的值的相加值eni(ni:成为稳定系统的维数)。
[式26]
W ( k &Delta; f ) = | H ( k &Delta; f ) | k &Delta; f , ( k = 0 , 1 , 2 , ... , N - 1 ) ... ( 26 )
[式27]
e n i = &Sigma; k = 0 N - 1 ( H ( k &Delta; f ) - H ^ n i ( k &Delta; f ) ) 2 &CenterDot; W ( k &Delta; f ) ... ( 27 )
使该加权误差平方和范数||eni||成为最小的维数ni为“根据权重函数,在频域中最适合于实际的频率响应的”稳定的系统维数n。此处,将如图5所示误差平方和范数分布41在用上述式(13)给出的误差平方和范数阈值42以下的维数中的最小的维数决定为系统维数n而输出(在图5的情况下系统维数n=n6)。
最后,在系统矩阵辨识部5中,根据从奇异值分解部3输出的第1正交矩阵U、第2正交矩阵V及奇异值σi(i=1,2,3,…)、和从系统维数决定部4a输出的系统维数n,用上述式(14)辨识线性离散时间系统的除了直通项Dd以外的系统矩阵Ad、Bd、Cd,将由直通项辨识部1辨识出的直通项Dd以及由系统矩阵辨识部5辨识出的系统矩阵Ad、Bd、Cd作为记述动态系统的线性离散时间系统输出。
这样,通过本实施方式的基于频率响应的系统辨识装置60,将现实的频率响应变换为脉冲响应,即使在根据该脉冲响应计算出的块汉克尔矩阵Hkl的奇异值σi(i=1,2,3,…)为平滑的单调减少,具有有效值的奇异值和成为在辨识中可忽略的微小的值的奇异值的边界变得不明确的情况下,也能够从系统维数n的决定排除反复试验,决定在频域中根据权重函数相对现实的动态系统一致度高的系统维数n以及辨识记述动态系统的线性离散时间系统。
另外,通过递归系统矩阵推测部31,能够降低用于决定相对现实的动态系统一致度高的系统维数n的运算量。
进而,根据逆傅立叶变换,将动态系统的频率响应变换为脉冲响应,从而能够根据动态系统的现实的频率响应,正确地辨识记述该系统的线性离散时间系统。
另外,通过系统稳定性评价部34,在现实的动态系统是稳定系统的情形明确的情况下,能够进行限定于稳定系统的线性离散时间系统的辨识。
另外,在本实施方式的基于频率响应的系统辨识装置60中,计算线性离散时间系统的系统特性而作为推测频率响应,将该频率响应和动态系统的实际的频率响应在频域中的加权误差平方和范数分布41在误差平方和范数阈值42以下的维数中的最小的维数决定为系统维数n。但是,本发明不限于此,也可以计算线性离散时间系统的系统特性而作为推测脉冲响应,根据该脉冲响应和对动态系统的频率响应进行逆傅立叶变换而得到的实际的脉冲响应在时域中的误差平方和,决定系统维数n。
在系统辨识装置中,系统维数决定部4a具备系统特性推测部32和系统维数推测部33,该系统特性推测部32针对属于搜索范围ni=(n1,n2,…,na)(其中n1<n2<…<na)的各维数计算线性离散时间系统的系统特性而作为推测频率响应并输出,该系统维数推测部33将从系统特性推测部32输出的线性离散时间系统的推测频率响应和对脉冲响应进行傅立叶变换而得到的实际的频率响应在频域中的误差平方和范数在阈值42以下的维数中的最小的维数决定为系统维数而输出,在该系统辨识装置中,即使在根据现实的脉冲响应计算出的块汉克尔矩阵的奇异值为平滑的单调减少,具有有效值的奇异值和成为可忽略的微小的值的奇异值的边界变得不明确的情况下,也能够从系统维数的决定排除反复试验,进行相对现实的动态系统在频域中一致度高的系统维数的决定、和记述动态系统的线性离散时间系统的辨识。
在系统辨识装置中,系统维数推测部33根据进行傅立叶变换而得到的实际的频率响应决定权重函数,计算对从系统特性推测部32输出的线性离散时间系统的推测频率响应和动态系统的实际的频率响应在频域中的误差平方值乘以权重函数而得到的值的相加值,将相加值的范数在阈值以下的维数中的最小的维数决定为系统维数而输出,在该系统辨识装置中,即使在根据现实的脉冲响应计算出的块汉克尔矩阵的奇异值为平滑的单调减少,具有有效值的奇异值和成为可忽略的微小的值的奇异值的边界变得不明确的情况下,也能够从系统维数的决定排除反复试验,决定在频域中根据权重函数相对现实的动态系统一致度高的系统维数以及辨识记述动态系统的线性离散时间系统。
实施方式4.
在本实施方式中,说明作为辨识对象的动态系统是DC伺服马达的情况。
图10是示出本发明的系统辨识装置的实施方式4的整体结构的框图。本实施方式的系统辨识装置70具备DC伺服马达71、频率响应变换部72、基于频率响应的系统辨识装置60以及控制系统或者滤波器参数设计部73。基于频率响应的系统辨识装置60与在实施方式3中使用图1以及图9说明的装置相同。
在本实施方式中,作为DC伺服马达71的输入电流[Arms],输入例如正弦扫描,将其作为系统输入。另外,取得角速度[rad/s]而作为系统输出。通过频率响应变换部72,将系统输入输出变换为频率响应H(kΔf)。频率响应变换部72通过例如FFT(Fast FourierTransform:快速傅立叶变换)等,根据下述的式(28)输出频率响应。
[式28]
基于频率响应的系统辨识装置60将该动态系统的频率响应以及系统维数的搜索范围作为输入,辨识记述DC伺服马达71的输入电流[Arms]至角速度[rad/s]的动态系统的线性离散时间系统,而输出到控制系统或者滤波器参数设计部73。此时,如ni=(1,2,…,50)那样,相对预测的系统维数具有充分的宽度地设定系统维数的搜索范围即可。
在基于频率响应的系统辨识装置60中,能够进行相对现实的动态系统一致度高的系统维数的决定、和记述动态系统的线性离散时间系统的辨识,所以能够将该线性离散时间系统用于伺服马达的控制系统或者滤波器的参数设计。
产业上的可利用性
如以上那样,本发明的系统辨识装置作为基于脉冲响应的系统辨识装置是有用的。

Claims (10)

1.一种系统辨识装置,将作为辨识对象的动态系统的脉冲响应、和被指定的系统维数的搜索范围作为输入,其特征在于,具备:
直通项辨识部,根据所述脉冲响应,辨识并输出记述所述动态系统的线性离散时间系统的直通项;
块汉克尔矩阵生成部,根据所述脉冲响应,生成并输出块汉克尔矩阵;
奇异值分解部,通过从所述块汉克尔矩阵生成部输出的所述块汉克尔矩阵的奇异值分解,输出将所述块汉克尔矩阵的左奇异矢量作为列矢量的第1正交矩阵、将所述块汉克尔矩阵的右奇异矢量作为列矢量的第2正交矩阵以及所述块汉克尔矩阵的奇异值;
系统维数决定部,根据所述第1正交矩阵、所述第2正交矩阵、所述奇异值以及所述搜索范围,辨识针对属于该搜索范围的各维数的所述线性离散时间系统的系统矩阵中所述直通项以外的系统矩阵,根据基于该系统矩阵以及所述直通项计算出的所述线性离散时间系统的系统特性和所述动态系统的实际的系统特性的比较,决定并输出系统维数;以及
系统矩阵辨识部,根据所述第1正交矩阵、所述第2正交矩阵、所述奇异值以及从所述系统维数决定部输出的所述系统维数,辨识所述线性离散时间系统的系统矩阵中所述直通项以外的系统矩阵,
所述系统辨识装置将由所述直通项辨识部辨识出的所述直通项以及由所述系统矩阵辨识部辨识出的所述系统矩阵作为所述线性离散时间系统而输出。
2.根据权利要求1所述的系统辨识装置,其特征在于,
所述系统维数决定部具备:
系统特性推测部,针对属于所述搜索范围的各维数,计算所述线性离散时间系统的系统特性作为推测脉冲响应而输出;以及
系统维数推测部,将从所述系统特性推测部输出的所述线性离散时间系统的所述推测脉冲响应和所述动态系统的实际的脉冲响应在时域中的误差平方和范数在阈值以下的维数中的最小的维数决定为系统维数而输出。
3.根据权利要求1所述的系统辨识装置,其特征在于,
所述系统维数决定部具备:
系统特性推测部,针对属于所述搜索范围的各维数,计算所述线性离散时间系统的系统特性作为推测频率响应而输出;以及
系统维数推测部,将从所述系统特性推测部输出的所述线性离散时间系统的推测频率响应和对所述脉冲响应进行傅立叶变换而得到的实际的频率响应在频域中的误差平方和范数在阈值以下的维数中的最小的维数决定为系统维数而输出。
4.根据权利要求3所述的系统辨识装置,其特征在于,
所述系统维数推测部根据进行所述傅立叶变换而得到的所述实际的频率响应来决定权重函数,
所述系统维数推测部计算对从所述系统特性推测部输出的所述线性离散时间系统的推测频率响应和所述动态系统的所述实际的频率响应在频域中的误差平方值乘以该权重函数而得到的值的相加值,将该相加值的范数在阈值以下的维数中的最小的维数决定为系统维数而输出。
5.根据权利要求2至4中的任意一项所述的系统辨识装置,其特征在于,
所述阈值是通过误差平方和容许值与误差平方和范数的最小值之积规定的。
6.根据权利要求1所述的系统辨识装置,其特征在于,
所述系统维数决定部具备递归系统矩阵推测部,该递归系统矩阵推测部使用与在所述搜索范围中比第1维数低1个等级的第2维数对应的系统矩阵的辨识结果、和从所述第1正交矩阵、所述第2正交矩阵以及所述奇异值中分别比所述第2维数大且在所述第1维数以下的左奇异矢量、右奇异矢量及奇异值,通过递归方法来辨识与所述第1维数对应的系统矩阵。
7.根据权利要求1所述的系统辨识装置,其特征在于,
与所述脉冲响应对应的信息是对所述动态系统进行了冲击振动时的系统输入以及系统输出,
所述系统辨识装置具备:
系统输入施加时刻确定部,确定所述系统输入具有有效值的多个时刻;以及
脉冲响应变换部,将与从所述系统输入施加时刻确定部输出的多个时刻对应的系统输入的相加值作为脉冲输入振幅,将从所述系统输入施加时刻确定部输出的多个时刻的最大值作为脉冲输入施加时刻,将该脉冲输入施加时刻以后的系统输出除以所述脉冲输入振幅而得到的信号作为所述脉冲响应而输出,
所述系统辨识装置将对作为所述辨识对象的所述动态系统进行了冲击振动的情况下的系统输入及系统输出以及所述搜索范围作为输入,辨识所述线性离散时间系统。
8.根据权利要求7所述的系统辨识装置,其特征在于,
所述系统输入施加时刻确定部将比例阈值和系统输入的最大值相乘而得到的值作为系统输入阈值,将系统输入在系统输入阈值以上的时刻确定为系统输入具有有效值的多个时刻。
9.根据权利要求1所述的系统辨识装置,其特征在于,
与所述脉冲响应对应的信息是所述动态系统的频率响应,
所述系统辨识装置具备逆傅立叶变换部,所述逆傅立叶变换部通过所述频率响应的逆傅立叶变换,输出对应的脉冲响应,
所述系统辨识装置将作为所述辨识对象的所述动态系统的频率响应以及所述搜索范围作为输入,辨识所述线性离散时间系统。
10.根据权利要求1所述的系统辨识装置,其特征在于,
所述系统维数决定部具备系统稳定性评价部,所述系统稳定性评价部针对属于所述搜索范围的各维数,评价所述线性离散时间系统的稳定性,
所述系统维数决定部根据与成为稳定系统的维数对应的所述线性离散时间系统的系统特性,决定系统维数。
CN201480074854.2A 2014-02-07 2014-11-05 系统辨识装置 Pending CN105960613A (zh)

Applications Claiming Priority (3)

Application Number Priority Date Filing Date Title
JP2014-022813 2014-02-07
JP2014022813 2014-02-07
PCT/JP2014/079339 WO2015118737A1 (ja) 2014-02-07 2014-11-05 システム同定装置

Publications (1)

Publication Number Publication Date
CN105960613A true CN105960613A (zh) 2016-09-21

Family

ID=53777556

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201480074854.2A Pending CN105960613A (zh) 2014-02-07 2014-11-05 系统辨识装置

Country Status (5)

Country Link
US (1) US10387116B2 (zh)
JP (1) JP6009106B2 (zh)
CN (1) CN105960613A (zh)
DE (1) DE112014006345T5 (zh)
WO (1) WO2015118737A1 (zh)

Families Citing this family (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107632522B (zh) * 2017-08-31 2020-06-19 南京理工大学 一种质子交换膜燃料电池非线性状态空间模型辨识方法
EP3757704A4 (en) * 2018-02-21 2021-04-07 NEC Corporation SYSTEMS IDENTIFICATION DEVICE, SYSTEM IDENTIFICATION PROCESS AND RECORDING MEDIA
CN109165432B (zh) * 2018-08-09 2022-12-13 厦门理工学院 一种基于部分奇异值和的磁共振波谱重建方法
CN110826017A (zh) * 2019-09-25 2020-02-21 中国地质大学(武汉) 一种基于参数优化Hankel矩阵和奇异值分解的信号去噪方法

Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2005078559A (ja) * 2003-09-03 2005-03-24 Fuji Electric Holdings Co Ltd 特性不明システムの同定装置
JP2006195543A (ja) * 2005-01-11 2006-07-27 Fuji Electric Holdings Co Ltd モデル同定装置およびモデル同定プログラム
CN101421923A (zh) * 2006-04-14 2009-04-29 国立大学法人岩手大学 系统辨识方法及程式、存储媒体和系统辨识装置
US20100198897A1 (en) * 2009-01-30 2010-08-05 Can Evren Yarman Deriving a function that represents data points
US20110125684A1 (en) * 2009-11-24 2011-05-26 Hussain Al-Duwaish Method for identifying multi-input multi-output hammerstein models
US20110161059A1 (en) * 2009-12-30 2011-06-30 Ankur Jain Method for Constructing a Gray-Box Model of a System Using Subspace System Identification
CN103501150A (zh) * 2013-10-12 2014-01-08 上海联孚新能源科技有限公司 一种内嵌式永磁同步电机参数辨识装置及方法

Family Cites Families (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPH0677211B2 (ja) 1985-05-22 1994-09-28 株式会社東芝 プラント・モデリング装置
JP2954660B2 (ja) 1990-05-30 1999-09-27 株式会社東芝 モデル予測制御装置
JP2005050057A (ja) 2003-07-31 2005-02-24 Fuji Electric Holdings Co Ltd 特性不明システムの伝達関数推定装置
JP2007260504A (ja) 2006-03-27 2007-10-11 Mitsubishi Chemical Engineering Corp システムモデル構築方法
JP4815391B2 (ja) 2007-05-15 2011-11-16 株式会社神戸製鋼所 モデルパラメータ推定演算装置及び方法、モデルパラメータ推定演算処理プログラム並びにそれを記録した記録媒体

Patent Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2005078559A (ja) * 2003-09-03 2005-03-24 Fuji Electric Holdings Co Ltd 特性不明システムの同定装置
JP2006195543A (ja) * 2005-01-11 2006-07-27 Fuji Electric Holdings Co Ltd モデル同定装置およびモデル同定プログラム
CN101421923A (zh) * 2006-04-14 2009-04-29 国立大学法人岩手大学 系统辨识方法及程式、存储媒体和系统辨识装置
US20100198897A1 (en) * 2009-01-30 2010-08-05 Can Evren Yarman Deriving a function that represents data points
US20110125684A1 (en) * 2009-11-24 2011-05-26 Hussain Al-Duwaish Method for identifying multi-input multi-output hammerstein models
US20110161059A1 (en) * 2009-12-30 2011-06-30 Ankur Jain Method for Constructing a Gray-Box Model of a System Using Subspace System Identification
CN103501150A (zh) * 2013-10-12 2014-01-08 上海联孚新能源科技有限公司 一种内嵌式永磁同步电机参数辨识装置及方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
HIROSHI OKU 等: "Identification experiment and conrtol design of an inverted pendulum system via closed-loop subspace model identification", 《JOURNAL OF THE SOCIETY OF INSTRUMENT AND CONTROL ENGINEERS》 *
YOSHITO HIRAI 等: "Model reduction for linear time-invariant discrete-time systems using matrix inequalities", 《PROCEEDINGS OF THE 48TH ANNUAL CONFERENCE OF THE INSTITUTE OF SYSTEMS,CONTROL AND INFORMATION ENGINEERS》 *

Also Published As

Publication number Publication date
US10387116B2 (en) 2019-08-20
DE112014006345T5 (de) 2016-10-20
US20170010861A1 (en) 2017-01-12
JP6009106B2 (ja) 2016-10-19
WO2015118737A1 (ja) 2015-08-13
JPWO2015118737A1 (ja) 2017-03-23

Similar Documents

Publication Publication Date Title
Subramani et al. Stochastic time-optimal path-planning in uncertain, strong, and dynamic flows
Staszewski et al. Wavelet-based frequency response function for time-variant systems—An exploratory study
Christensen et al. Simulating weather regimes: Impact of stochastic and perturbed parameter schemes in a simple atmospheric model
CN107247687B (zh) 一种基于特征正交分解的非平稳随机过程快速模拟方法
CN101916241B (zh) 一种基于时频分布图的时变结构模态频率辨识方法
Hu et al. Model order determination and noise removal for modal parameter estimation
CN107290589A (zh) 基于短时分数阶傅里叶变换的非线性信号时频分析方法
CN104793253A (zh) 基于数学形态学的航空电磁数据去噪方法
CN105960613A (zh) 系统辨识装置
Gunn et al. On validating numerical hydrodynamic models of complex tidal flow
Liu et al. Weak-mode identification and time-series reconstruction from high-level noisy measured data of offshore structures
CN105205461A (zh) 一种用于模态参数识别的信号降噪方法
Yue et al. A Bayesian wavelet packet denoising criterion for mechanical signal with non-Gaussian characteristic
Hertwig et al. LES validation of urban flow, part II: eddy statistics and flow structures
CN105991492A (zh) 跳频信号识别方法
CN105844217A (zh) 一种基于量测驱动新生目标强度估计的phd多目标跟踪方法
Dash et al. A hybrid unscented filtering and particle swarm optimization technique for harmonic analysis of nonstationary signals
CN110782041B (zh) 一种基于机器学习的结构模态参数识别方法
CN107329115A (zh) 基于gcrbf网络的lfm信号参数估计方法
Garnier et al. The contsid toolbox: A software support for data-based continuous-time modelling
Uyanık et al. Parametric identification of hybrid linear-time-periodic systems
Farzampour et al. Unsupervised identification of arbitrarily-damped structures using time-scale independent component analysis: Part I
Marano Non-stationary stochastic modulation function definition based on process energy release
Olsen et al. Condensation of the correlation functions in modal testing
CN103698757B (zh) 低频段雷达目标微动特性估计方法

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
RJ01 Rejection of invention patent application after publication

Application publication date: 20160921

RJ01 Rejection of invention patent application after publication