WO2015118736A1 - システム同定装置 - Google Patents

システム同定装置 Download PDF

Info

Publication number
WO2015118736A1
WO2015118736A1 PCT/JP2014/079257 JP2014079257W WO2015118736A1 WO 2015118736 A1 WO2015118736 A1 WO 2015118736A1 JP 2014079257 W JP2014079257 W JP 2014079257W WO 2015118736 A1 WO2015118736 A1 WO 2015118736A1
Authority
WO
WIPO (PCT)
Prior art keywords
dimension
output
input
matrix
dynamic
Prior art date
Application number
PCT/JP2014/079257
Other languages
English (en)
French (fr)
Inventor
光伯 齊藤
百合夏 金井
Original Assignee
三菱電機株式会社
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 三菱電機株式会社 filed Critical 三菱電機株式会社
Priority to DE112014006135.2T priority Critical patent/DE112014006135T5/de
Priority to CN201480074855.7A priority patent/CN105960614B/zh
Priority to US15/114,729 priority patent/US20160342731A1/en
Priority to JP2015561158A priority patent/JP6009105B2/ja
Publication of WO2015118736A1 publication Critical patent/WO2015118736A1/ja

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/10—Geometric CAD
    • G06F30/17—Mechanical parametric or variational design
    • G—PHYSICS
    • G05—CONTROLLING; REGULATING
    • G05B—CONTROL OR REGULATING SYSTEMS IN GENERAL; FUNCTIONAL ELEMENTS OF SUCH SYSTEMS; MONITORING OR TESTING ARRANGEMENTS FOR SUCH SYSTEMS OR ELEMENTS
    • G05B13/00—Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion
    • G05B13/02—Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion electric
    • G05B13/04—Adaptive 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/042—Adaptive 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/044—Adaptive 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

Definitions

  • the present invention relates to a system identification apparatus that constructs a mathematical model of a system based on system input / output when a pseudo-random input is applied to a target dynamic system.
  • Non-Patent Document 1 As a conventional system identification device using pseudo-random input, for example, there is one based on the N4SID method described in Non-Patent Document 1.
  • this N4SID method a block Hankel related to system input based on system input / output when a pseudo-random input is applied to a dynamic system described by a linear discrete time system (A d , B d , C d, D d ).
  • a matrix (U p , U f ) and a block Hankel matrix (Y p , Y f ) relating to the system output are generated.
  • the data matrix obtained by combining the block Hankel matrices is subjected to LQ decomposition, and a parallel projection ⁇ is generated from the partial matrix obtained by the LQ decomposition and the block Hankel matrices U p and Y p .
  • the parallel projection ⁇ is decomposed into singular values, and the number of singular values with significant values is determined as the system dimension. From the result of the singular value decomposition and the determined system dimensions, the state vector ( ⁇ X K , ⁇ X K + 1 ) is calculated.
  • Patent Document 1 As another example of a conventional system identification apparatus using pseudo-random input, for example, there are an exposure apparatus and a vibration isolation apparatus, a system identification apparatus, and a method thereof described in Patent Document 1.
  • the dynamics are determined by a subspace method represented by the N4SID method. Identify the equation of state of a dynamic system. At this time, by matching the system dimension of the state equation to be identified with the system dimension determined from the equation of motion of the dynamic system, it is possible to compare the characteristic equation based on the equation of motion with the characteristic equation based on the identified state equation. An unknown physical parameter included in the equation is identified.
  • the system dimension of the target dynamic system is determined from the number of singular values having significant values or the system dimension determined from the dynamic system equation of motion.
  • the present invention has been made in view of the above, and the singular value of the parallel projection ⁇ calculated from the actual system input / output is gently monotonously decreased. Therefore, the singular value having a significant value and the negligible minute value can be ignored. It is an object of the present invention to obtain a system identification apparatus that can eliminate the trial and error from the determination of the system dimension and determine the optimum system dimension even when the boundary with the singular value as the value becomes unclear.
  • an object of the present invention is to obtain a system identification device that can be identified by limiting to a stable system when it is clear that an actual dynamic system is stable.
  • the system identification apparatus is configured to provide system input / output and a specified system dimension when a pseudo random input is applied to a dynamic system to be identified.
  • a system identification device having a search range as input, a system input / output extraction unit for extracting identification input / output data to be applied to identification from system input / output of the dynamic system, and a block based on the identification input / output data
  • a block Hankel matrix generator for generating a Hankel matrix, an input / output vector generator for generating an input vector and an output vector of the dynamic system based on the block Hankel matrix, and a data matrix by combining the block Hankel matrix
  • An LQ decomposition unit that generates a partial matrix obtained by LQ decomposition of the data matrix, and the unit
  • a parallel projection generating unit that generates a parallel projection based on the matrix and the block Hankel matrix; and a first orthogonal matrix having the singular vector of the parallel projection as a column vector by singular value de
  • the singular value of the parallel projection calculated from the actual system input / output is gently monotonously decreased, and therefore, the singular value having a significant value and the negligible small value Even when the boundary with the singular value becomes unclear, trial and error is eliminated from the determination of the system dimension, and the optimal system dimension can always be determined and the linear discrete-time system describing the dynamic system can be identified It becomes.
  • FIG. 1 is a block diagram showing the overall configuration of the system identification apparatus according to the first and second embodiments.
  • FIG. 2 is a schematic diagram showing time waveforms of system input / output in the system identification apparatus of the first embodiment.
  • FIG. 3 is a schematic diagram illustrating a relationship between a singular value and a dimension of parallel projection in the system identification apparatus according to the first and second embodiments.
  • FIG. 4 is a block diagram showing an internal configuration of a system dimension determining unit in the system identification apparatus of the first embodiment.
  • FIG. 5 is a schematic diagram showing the relationship between the norm of error square sum and dimension in the time domain or frequency domain of the identified linear discrete-time system in the system identification apparatus according to the first and second embodiments. .
  • FIG. 1 is a block diagram showing the overall configuration of the system identification apparatus according to the first and second embodiments.
  • FIG. 2 is a schematic diagram showing time waveforms of system input / output in the system identification apparatus of the first embodiment.
  • FIG. 3 is a schematic diagram
  • FIG. 6 is a schematic diagram showing system input / output time waveforms when the dynamic system in the system identification apparatus of the second embodiment is subjected to M-sequence excitation.
  • FIG. 7 is a block diagram showing an internal configuration of a system dimension determining unit in the system identification apparatus of the second embodiment.
  • FIG. 8 is a block diagram showing an overall configuration according to the third embodiment.
  • FIG. 1 is a block diagram showing the overall configuration of the system identification apparatus according to the first embodiment
  • FIG. 2 is a schematic diagram showing system input / output time waveforms in the system identification apparatus of the first embodiment.
  • the input / output vector generation unit 3 generates an input vector ⁇ U K
  • the LQ decomposition unit 4 generates a data matrix obtained by combining the block Hankel matrices U p , U f , Y p , Y f , and generates and outputs partial matrices L 22 , L 32 obtained by performing LQ decomposition on the data matrix.
  • the parallel projection generation unit 5 is based on the partial matrices L 22 and L 32 output from the LQ decomposition unit 4 and the block Hankel matrices U p and Y p output from the block Hankel matrix generation unit 2. Generate a parallel projection ⁇ .
  • the singular value decomposition unit 6 performs singular value decomposition on the parallel projection ⁇ output from the parallel projection generation unit 5, and a first orthogonal matrix U having the left singular vector of the parallel projection ⁇ as a column vector, and the right singularity of the parallel projection ⁇ .
  • the system dimension n is determined by comparison with (described in FIG. 1 as system characteristics of the dynamic system).
  • the system matrix identification unit 9 includes the dynamic system input vector ⁇ U K
  • FIG. FIG. 5 is a block diagram showing an internal configuration of the system dimension determination unit 7 in the system identification apparatus 10 of FIG. 1, and FIG. 5 shows the system output of the linear discrete-time system identified in the system identification apparatus 10 of Embodiment 1 and the dynamics norm of the error sum of squares in the time domain and the actual system output of the system
  • a dimension n i (i 1,2, ... , a) is a schematic diagram showing the relationship between.
  • the relationship shown in the singular value distribution 21 is obtained.
  • the boundary between the singular value having a significant value and the singular value that becomes a negligible minute value becomes unclear, and the optimum system dimension n is not always determined. Therefore, there arises a problem that trial and error are required for determining the system dimension n.
  • the system dimension determination unit 7 executes the process shown in FIG. Specifically, it is as follows.
  • the system dimension determination unit 7 includes a recursive system matrix estimation unit 31, a system characteristic estimation unit 32, and a system dimension estimation unit 33.
  • n i (n 1 , n 2 ,..., N a )
  • n 1 ⁇ n 2 ⁇ ... ⁇ N a the system matrix A d, n i , B d, n i , C d, n i , D d, n i output from the recursive system matrix estimation unit 31
  • the system input u (jT S ) to the dynamic system is composed of pseudo random inputs
  • the system input u (jT S ) and the system output y (jT S ) corresponding to the above [Equation 1] are For example, time waveforms such as the system input 11 and the system output 12 shown in FIG.
  • system input / output extraction unit 1 extracts the system input 11 and the system output 12 after the pseudo random input application time j min T S using the following equations.
  • system input / output extraction unit 1 uses the values extracted using the above [Equation 3] as identification input data u id (jT S ) and identification output data y id (jT S ), Remove system static time domain data before applying pseudo-random input from the system input / output of the target dynamic system.
  • the input / output vector generator 3 inputs the dynamic system input vector ⁇ U K
  • the LQ decomposition unit 4 generates a data matrix given by the following expression obtained by combining the block Hankel matrices U p , U f , Y p , Y f .
  • the LQ decomposition unit 4 performs LQ decomposition on the data matrix as in the following equation, and outputs partial matrices L 22 and L 32 from the elements of the LQ decomposed matrix.
  • the parallel projection generation unit 5 is defined by the following equation based on the partial matrices L 22 and L 32 output from the LQ decomposition unit 4 and the block Hankel matrices U p and Y p output from the block Hankel matrix generation unit 2. Generate a parallel projection ⁇ of the dynamic system to be performed.
  • the singular value decomposition unit 6 performs a singular value decomposition on the parallel projection ⁇ represented by the above equation, whereby the first orthogonal matrix U and parallel using the left singular vector u j of the parallel projection ⁇ given by the following equation as a column vector.
  • n of the target dynamic system is based on the following relationship in which n have a significant value in the singular value of the parallel projection ⁇ , and n + 1 and beyond are sufficiently smaller than those. Can be determined.
  • the relationship shown in the distribution 21 is obtained.
  • the boundary ⁇ n >> ⁇ n + 1 between the singular value having a significant value and the singular value that becomes a negligible minute value becomes unclear. Therefore, the conventional method does not always determine the optimum system dimension n, and there is a problem that trial and error are required to determine the optimum system dimension n.
  • the system dimension determination unit 7 determines the optimum system dimension n under the premise that “the system is most suitable for actual system input / output in the time domain”.
  • K output from the unit 3, and the system dimension search range n i (n 1 , n 2 , ..., n specified by the operator a ) (where n 1 ⁇ n 2 ⁇ ...
  • the recursive system matrix estimation unit 31 performs a system dimension search range n i (n 1 , n 2 ,..., N a ) (provided that n 1 is specified). ⁇ n 2 ⁇ ...
  • of the sum of squared errors expressed by the above equation is the system dimension n that “best fits actual system input / output in the time domain”.
  • decreases monotonously with the increase of dimension n i regardless of the noise level, and exceeds a certain dimension as shown in FIG. The value is almost constant. Therefore, here, in order to avoid that the estimated value of the system dimension n becomes higher than necessary, an error square sum norm threshold 42 given by the following equation is defined.
  • the system dimension estimating unit 33 determines and outputs the smallest dimension among the dimensions in which the distribution 41 of the norm of error square sum
  • is equal to or less than the error square sum norm threshold value 42 as the system dimension n. (In the example of FIG. 5, system dimension n n 6 ).
  • system matrix identification unit 9 outputs the dynamic system input vector ⁇ U K
  • the presence of the recursive system matrix estimation unit 31 makes it possible to reduce the amount of calculation for determining the system dimension n having a high degree of coincidence with an actual dynamic system.
  • the system output when the actual identification input data is applied to the linear discrete time system is calculated as a system characteristic, and the system output and the actual dynamic system are calculated.
  • the system dimension n is determined as the system dimension n.
  • the system characteristic of the linear discrete-time system is calculated as a frequency response, and the frequency response and the actual frequency response obtained from the input / output data for identification of the dynamic system are calculated.
  • the system dimension n may be determined based on the error sum of squares in the frequency domain.
  • a weighting function is further determined based on the actual frequency response of the dynamic system, and the weight of the error square value in the frequency domain between the frequency response of the linear discrete-time system and the actual frequency response of the dynamic system is determined.
  • the system dimension n may be determined based on the added value obtained by multiplying the function.
  • FIG. 2 a system identification device according to Embodiment 2 will be described.
  • and dimension n i (i 1,2, ..., a) between the frequency response of the linear discrete-time system and the actual frequency response of the dynamic system
  • the schematic diagram showing the relationship is the same as FIG. 1, FIG. 3, and FIG. 5 used in the description of the first embodiment.
  • FIG. 6 is a schematic diagram showing system input / output time waveforms when the dynamic system in the system identification apparatus of the second embodiment is subjected to M-sequence excitation.
  • FIG. 7 is a block diagram showing an internal configuration of the system dimension determining unit 7 in the system identification apparatus of the second embodiment. 7, components denoted by the same reference numerals as those in FIG. 4 are the same or equivalent components as those in the first embodiment, and a system stability evaluation unit 34 is added.
  • the system dimension search range n i (n 1 , n 2 ,..., N a ) (where n 1 ⁇ For each dimension n i belonging to n 2 ⁇ ... ⁇ n a ), the system matrices A d, n i , B d, n i , C d, n i , D identified by the recursive system matrix estimation unit 31 are used. Based on d and n i , the system stability evaluation unit 34 evaluates the stability of the linear discrete-time system.
  • the system characteristic estimator 32 applies the system matrix A d, n i , B d, n i output from the recursive system matrix estimator 31 to the dimension determined to be a stable system by the system stability evaluator 34. , C d, n i , D d, n i , calculate the frequency response for the identified linear discrete-time system.
  • the system dimension estimation unit 33 determines a weighting function based on the actual frequency response (described as system characteristics of the dynamic system in FIG. 7) obtained from the system input / output of the dynamic system, and outputs the weight function from the system characteristic estimation unit 32.
  • Added value en i (n i : stable system) obtained by multiplying the squared error value in the frequency domain between the frequency response of the linear discrete-time system and the actual frequency response of the dynamic system 5), and the distribution 41 of the norm
  • the system input / output extraction unit 1 includes a preset ratio threshold and the maximum value of the system input 11 (u (jT S )).
  • [Equation 2] is multiplied by the system input threshold 13, and the minimum value of the time when the absolute value of the system input 11 is equal to or greater than the system input threshold 13 is specified as the M-sequence signal application time j min T s (FIG. 6).
  • j min T s 2T s ).
  • system input / output extraction unit 1 extracts the system input 11 and the system output 12 after the M-sequence signal application time j min T s by [Equation 3], and inputs each of the extracted input data u id (jT S ) And identification output data y id (jT S ), the system static time domain data before application of the M-sequence signal is removed from the system input / output of the target dynamic system.
  • the block Hankel matrix generation unit 2 generates block Hankel matrices U p , U f , Y p , Y f given by [Equation 4], and the input / output vector generation unit 3 [
  • K given by Equation 5] are generated, and the LQ decomposition unit 4 combines the block Hankel matrices U p , U f , Y p , Y f
  • the parallel projection generation unit 5 generates a parallel projection ⁇ of the dynamic system defined by [Equation 8], and the singular value decomposition unit 6 performs singular value decomposition on the generated parallel projection ⁇ , ],
  • the system dimension estimation unit 33 determines the optimum system dimension n under the premise that “the most suitable frequency response in the frequency domain is most suitable”. Specifically, it is as follows.
  • is the stable system dimension n that “best fits the actual frequency response in the frequency domain according to the weight function”.
  • the system matrix identification unit 9 inputs the dynamic system input vector ⁇ U K
  • the system dimension n having a high degree of coincidence can be determined for the system in accordance with the weight function in the frequency domain, and a linear discrete-time system describing a dynamic system can be identified.
  • the presence of the recursive system matrix estimation unit 31 makes it possible to reduce the amount of calculation for determining the system dimension n having a high degree of coincidence with an actual dynamic system.
  • system stability evaluation unit 34 makes it possible to identify a linear discrete-time system limited to a stable system when it is clear that the actual dynamic system is a stable system.
  • the system characteristic of the linear discrete time system is calculated as the frequency response, and the frequency response and the actual frequency response obtained from the input / output data for identification of the dynamic system are calculated.
  • the norm distribution 41 of the sum of squared errors in the frequency domain is equal to or less than a preset threshold value 42
  • the smallest dimension is determined as the system dimension n.
  • the present invention is not limited to this, and the system output when the actual identification input data is applied to the linear discrete-time system is calculated as a system characteristic.
  • the system dimension n may be determined based on the sum of squared errors in the time domain with the identification output data.
  • FIG. 8 is a block diagram showing an overall configuration according to the third embodiment.
  • system identification device 10 shown in FIG. 8 has the same or equivalent configuration as system identification device 10 according to the first embodiment shown in FIG.
  • the system identification device 10 receives the system input / output and the system dimension search range as inputs, and identifies a linear discrete time system describing the DC servo motor 51.
  • the system identification device 10 servos the linear discrete time system in order to determine a system dimension having a high degree of coincidence with an actual dynamic system and to identify a linear discrete time system describing the dynamic system. It can be used for parameter design in a motor control system, filter parameter design, and the like.

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Geometry (AREA)
  • Evolutionary Computation (AREA)
  • Health & Medical Sciences (AREA)
  • Automation & Control Theory (AREA)
  • Software Systems (AREA)
  • Medical Informatics (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Artificial Intelligence (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Optimization (AREA)
  • General Engineering & Computer Science (AREA)
  • Computer Hardware Design (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Computational Mathematics (AREA)
  • Complex Calculations (AREA)
  • Testing And Monitoring For Control Systems (AREA)
  • Feedback Control In General (AREA)

Abstract

　システム同定装置は、指定されたシステム次元の探索範囲に属するそれぞれの次元に対して再帰的手法によって線形離散時間システムを同定し、各次元に対応する線形離散時間システムに対して、実際の同定用入力データを適用した場合のシステム出力をシステム特性として算出し、当該システム出力と、動的システムの実際の同定用出力データとの時間領域における誤差２乗和のノルム分布（４１）が閾値（４２）以下となる次元のうち、最小の次元をシステム次元ｎとして決定すると共に、動的システムの入出力ベクトルおよび、決定されたシステム次元を用いて生成された状態ベクトルに基づいて線形離散時間システムのシステム行列を同定する。

Description

システム同定装置
　この発明は、対象とする動的システムに対して擬似ランダム入力を印加した場合のシステム入出力に基づいて、当該システムの数学モデルを構築するシステム同定装置に関する。
　従来の擬似ランダム入力によるシステム同定装置として、例えば非特許文献１に記載されたＮ４ＳＩＤ法に基づくものがある。このＮ４ＳＩＤ法では、線形離散時間システム（Ａd,Ｂd,Ｃd,Ｄd）で記述される動的システムに擬似ランダム入力を印加した場合のシステム入出力に基づいて、システム入力に関するブロックハンケル行列（Ｕp,Ｕf）、およびシステム出力に関するブロックハンケル行列（Ｙp,Ｙf）を生成し、さらにブロックハンケル行列（Ｕf,Ｙf）に基づいて、入出力ベクトル（~ＵK|K,~ＹK|K）を生成する。なお、「~」の表記は、本来であれば「Ｕ」の文字の上部に横線（オーバーバー）を引くべきところであるが、その表記ができない。このため、本明細書では、イメージで挿入する数式部分を除き、横線（オーバーバー）を「~」で代替する。
　次に、上記ブロックハンケル行列を結合したデータ行列をＬＱ分解し、ＬＱ分解により得られた部分行列と、ブロックハンケル行列Ｕp,Ｙpから平行射影Θを生成する。この平行射影Θを特異値分解し、有意な値を持つ特異値の個数をシステム次元として決定して、特異値分解の結果と決定したシステム次元から、動的システムの状態ベクトル（~ＸK,~ＸK+1）を算出する。最後に、入出力ベクトル（~ＵK|K,~ＹK|K）および状態ベクトル（~ＸK,~ＸK+1）に対して最小２乗法を適用することで、動的システムを記述する線形離散時間システム（Ａd,Ｂd,Ｃd,Ｄd）を同定している。
　また、従来の擬似ランダム入力によるシステム同定装置の他の例として、例えば特許文献１に記載された露光装置及び除振装置、システム同定装置及びその方法がある。
　この露光装置及び除振装置、システム同定装置及びその方法では、対象とする動的システムに擬似ランダム入力を印加した場合のシステム入出力に基づいて、Ｎ４ＳＩＤ法に代表される部分空間法により、動的システムの状態方程式を同定する。このとき、同定する状態方程式のシステム次元を、動的システムの運動方程式から定まるシステム次元に一致させることで、運動方程式に基づく特性方程式と、同定した状態方程式に基づく特性方程式との比較により、運動方程式に含まれる未知の物理パラメータを同定している。
特開２０００－８２６６２号公報
システム同定－部分空間法からのアプローチ、朝倉書店　ｐｐ．１１７－１２０
　このような擬似ランダム入力によるシステム同定装置では、有意な値を持つ特異値の個数、または動的システムの運動方程式から定まるシステム次元から、対象とする動的システムのシステム次元を決定している。
　しかしながら、現実のシステム入出力から算出した平行射影Θの特異値は、なだらかな単調減少となる場合が多数存在し、この場合は有意な値を持つ特異値と、無視できる微小な値となる特異値との境界が不明確となる。したがって、非特許文献１に記載された従来のシステム同定装置では、システム次元の決定が作業者の判断に依存し、常に最適なシステム次元を決定しているとは限らないか、もしくはシステム次元の決定に関して試行錯誤が必要になるという問題があった。
　また、動的システムのモデリングによって得られる運動方程式では、動的システムの実際の動特性を全て記述することは困難であり、一般に、“運動方程式から定まるシステム次元＜動的システムの実際のシステム次元”となる。したがって、特許文献１に記載された従来のシステム同定装置では、そもそも動的システムを記述する最適なシステム次元を決定することができないという問題があった。
　加えて、従来の擬似ランダム入力によるシステム同定装置では、同定結果として得られる線形離散時間システム（Ａd,Ｂd,Ｃd,Ｄd）の安定性は全く考慮されていないため、現実の動的システムが安定であるにも関わらず、不安定システムとして同定される場合があるという問題もあった。
　本発明は、上記に鑑みてなされたものであって、現実のシステム入出力から算出した平行射影Θの特異値がなだらかな単調減少となり、したがって有意な値を持つ特異値と、無視できる微小な値となる特異値との境界が不明確となる場合においても、システム次元の決定から試行錯誤を排除し、最適なシステム次元を決定することができるシステム同定装置を得ることを目的とする。
　また、本発明は、現実の動的システムが安定であることが明確である場合に、安定システムに限定して同定することができるシステム同定装置を得ることを目的とする。
　上述した課題を解決し、目的を達成するため、本発明に係るシステム同定装置は、同定対象とする動的システムに対して擬似ランダム入力を印加した場合のシステム入出力および指定されたシステム次元の探索範囲を入力とするシステム同定装置において、前記動的システムのシステム入出力から同定に適用する同定用入出力データを抽出するシステム入出力抽出部と、前記同定用入出力データに基づいて、ブロックハンケル行列を生成するブロックハンケル行列生成部と、前記ブロックハンケル行列に基づいて、前記動的システムの入力ベクトルおよび出力ベクトルを生成する入出力ベクトル生成部と、前記ブロックハンケル行列を結合してデータ行列を生成し、当該データ行列をＬＱ分解した部分行列を出力するＬＱ分解部と、前記部分行列と前記ブロックハンケル行列に基づいて、平行射影を生成する平行射影生成部と、前記平行射影の特異値分解により、前記平行射影の特異ベクトルを列ベクトルとする第１の直交行列、当該平行射影の右特異ベクトルを列ベクトルとする第２の直交行列および当該平行射影の特異値を出力する特異値分解部と、前記第２の直交行列および前記特異値と、前記動的システムの入力ベクトルおよび出力ベクトルと、前記探索範囲に基づいて、当該探索範囲に属するそれぞれの次元に対して動的システムを記述する線形離散時間システムのシステム行列を同定し、さらに当該システム行列に基づいて算出した線形離散時間システムのシステム特性と、動的システムの実際のシステム特性との比較から、システム次元を決定するシステム次元決定部と、前記第２の直交行列および特異値と、前記決定されたシステム次元に基づいて、前記動的システムの状態ベクトルを生成する状態ベクトル生成部と、前記動的システムの入力ベクトルおよび出力ベクトルならびに前記動的システムの状態ベクトルに基づいて、当該動的システムを記述する線形離散時間システムのシステム行列を同定するシステム行列同定部と、を備え、前記同定されたシステム行列を、前記動的システムを記述する線形離散時間システムとして出力することを特徴とする。
　この発明によれば、同定対象とする動的システムにおいて、現実のシステム入出力から算出した平行射影の特異値がなだらかな単調減少となり、したがって有意な値を持つ特異値と、無視できる微小な値となる特異値との境界が不明確となる場合においても、システム次元の決定から試行錯誤を排除し、常に最適なシステム次元の決定と、動的システムを記述する線形離散時間システムの同定が可能となる。
図１は、実施の形態１および実施の形態２に係るシステム同定装置の全体構成を示すブロック線図である。 図２は、実施の形態１のシステム同定装置におけるシステム入出力の時間波形を示す概略図である。 図３は、実施の形態１および実施の形態２に係るシステム同定装置における平行射影の特異値と次元との関係を示す概略図である。 図４は、実施の形態１のシステム同定装置におけるシステム次元決定部の内部構成を示すブロック線図である。 図５は、実施の形態１および実施の形態２に係るシステム同定装置において、同定した線形離散時間システムの時間領域または周波数領域における誤差２乗和のノルムと次元との関係を示す概略図である。 図６は、実施の形態２のシステム同定装置における動的システムをＭ系列加振した場合のシステム入出力の時間波形を示す概略図である。 図７は、実施の形態２のシステム同定装置におけるシステム次元決定部の内部構成を示すブロック線図である。 図８は、実施の形態３に係る全体構成を示すブロック線図である。
　以下、添付図面を参照し、本発明の実施の形態に係るシステム同定装置について説明する。なお、以下に示す実施の形態により本発明が限定されるものではない。
実施の形態１．
　図１は、実施の形態１に係るシステム同定装置の全体構成を示すブロック線図であり、図２は、実施の形態１のシステム同定装置におけるシステム入出力の時間波形を示す概略図である。
　実施の形態１に係るシステム同定装置１０では、図１および図２に示すように、同定対象とする動的システムに対して擬似ランダム入力を印加した場合のシステム入力１１（ｕ(jＴS)(ｊ=0,1,2,…)）および、システム出力１２（ｙ(jＴS)(ｊ=0,1,2,…)）を入力とする。
　システム入出力抽出部１は、予め設定された比率閾値とシステム入力１１の最大値とを乗じた値で決定されるシステム入力閾値１３に対して、システム入力１１の絶対値がシステム入力閾値１３以上となる時刻の最小値を擬似ランダム入力印加時刻（図２では３ＴS)として、擬似ランダム入力印加時刻以降のシステム入力１１およびシステム出力１２を、それぞれ同定用入力データ（ｕid(jＴS)(ｊ=0,1,2,…)）および、同定用出力データ（ｙid(jＴS)(ｊ=0,1,2,…)）として抽出して出力する。
　ブロックハンケル行列生成部２は、システム入出力抽出部１から出力される同定用入力データｕid(jＴS)(ｊ=0,1,2,…)および、同定用出力データｙid(jＴS)(ｊ=0,1,2,…）に基づいて、ブロックハンケル行列Ｕp,ＵfおよびＹp,Ｙfを生成する。
　入出力ベクトル生成部３は、ブロックハンケル行列Ｕp,Ｕf,Ｙp,Ｙfに基づいて、動的システムの入力ベクトル~ＵK|Kおよび出力ベクトル~ＹK|Kを生成する。
　ＬＱ分解部４は、ブロックハンケル行列Ｕp,Ｕf,Ｙp,Ｙfを結合したデータ行列を生成し、当該データ行列をＬＱ分解した部分行列Ｌ22,Ｌ32を生成して出力する。
　平行射影生成部５は、ＬＱ分解部４から出力される部分行列Ｌ22,Ｌ32と、ブロックハンケル行列生成部２から出力されるブロックハンケル行列Ｕp,Ｙpに基づいて、動的システムの平行射影Θを生成する。
　特異値分解部６は、平行射影生成部５から出力される平行射影Θを特異値分解し、平行射影Θの左特異ベクトルを列ベクトルとする第１の直交行列Ｕ、平行射影Θの右特異ベクトルを列ベクトルとする第２の直交行列Ｖおよび、平行射影Θの特異値σi(i=1,2,3…)を出力する。
　システム次元決定部７は、特異値分解部６から出力される第２の直交行列Ｖおよび特異値σi(i=1,2,3…)と、入出力ベクトル生成部３から出力される動的システムの入力ベクトル~ＵK|Kおよび出力ベクトル~ＹK|Kならびに、作業者が指定したシステム次元の探索範囲ｎi=(ｎ1,ｎ2,…,ｎa)（ただし、ｎ1<ｎ2<…<ｎa）に基づいて、当該探索範囲に属する次元ｎi(i=1,2,…,a)のそれぞれに対して動的システムを記述する線形離散時間システムのシステム行列を同定する。さらに当該システム行列に基づいて、当該探索範囲に属する次元ｎi(i=1,2,…,a)のそれぞれに対応する線形離散時間システムに対して実際の同定用入力データｕid(jＴS)(ｊ=0,1,2,…)を適用した場合のシステム出力を算出し、動的システムの実際の同定用出力データｙid(jＴS)(ｊ=0,1,2,…）（図１では動的システムのシステム特性と記載）との比較から、システム次元ｎを決定する。
　状態ベクトル生成部８では、特異値分解部６から出力される第２の直交行列Ｖおよび特異値σi(i=1,2,3…)および、システム次元決定部７から出力されるシステム次元ｎに基づいて、動的システムの状態ベクトル~ＸK+1,~ＸKを生成する。
　システム行列同定部９は、入出力ベクトル生成部３から出力される動的システムの入力ベクトル~ＵK|Kおよび、出力ベクトル~ＹK|Kと、状態ベクトル生成部８から出力される動的システムの状態ベクトル~ＸK+1,~ＸKに基づいて、動的システムを記述する線形離散時間システムのシステム行列Ａd,Ｂd,Ｃd,Ｄdを同定して出力する。
　図３は、実施の形態１のシステム同定装置１０における平行射影Θの特異値σiと次元(i=1,2,3…)の関係を示す概略図であり、図４は、実施の形態１のシステム同定装置１０におけるシステム次元決定部７の内部構成を示すブロック線図であり、図５は、実施の形態１のシステム同定装置１０における同定した線形離散時間システムのシステム出力と、動的システムの実際のシステム出力との時間領域における誤差２乗和のノルム||ｅn||と次元ｎi(i=1,2,…,a)との関係を示す概略図である。
　図３に示すように、動的システムのシステム入出力から算出した平行射影Θの特異値σi(i=1,2,3…)は、理想的には次元(i=1,2,3…)に対して、例えば特異値分布２１に示す関係となる。この場合、有意な値を持つ特異値の個数を明確に規定することができ、当該個数が動的システムのシステム次元ｎに対応する（図３の場合はシステム次元ｎ＝４）。
　一方、観測ノイズ等の影響を受ける現実のシステム入出力に基づいて算出した特異値σiは、次元(i=1,2,3…)に対して例えば特異値分布２２に示す関係となるため、有意な値を持つ特異値と、無視できる微小な値となる特異値との境界が不明確となり、常に最適なシステム次元ｎを決定しているとは限らない。よって、システム次元ｎの決定に関して試行錯誤が必要になるという問題が発生する。
　そこで、実施の形態１に係るシステム同定装置１０では、システム次元決定部７により、図４に示す処理を実行する。具体的には、以下の通りである。
　システム次元決定部７には、再帰的システム行列推定部３１、システム特性推定部３２およびシステム次元推定部３３が設けられる。
　再帰的システム行列推定部３１は、作業者が予め指定したシステム次元の探索範囲ｎi(ｎ1,ｎ2,…,ｎa)（ただし、ｎ1<ｎ2<…<ｎa）に属する第１の次元ｎiに対応するシステム行列の同定に関して、第１の次元ｎiよりも１段階低い第２の次元ｎi-1に対応したシステム行列Ａd,ni-1,Ｂd,ni-1,Ｃd,ni-1,Ｄd,ni-1の同定結果と、特異値分解部６から出力される第２の直交行列Ｖおよび特異値σi(i=1,2,3…)のうち、それぞれ第２の次元ｎi-1よりも大きく、第１の次元ｎi以下となる右特異ベクトルｖjおよび特異値σj(j=ｎi-1+1,ｎi-1+2,…,ｎi)と、入出力ベクトル生成部３から出力される動的システムの入力ベクトル~ＵK|K、および出力ベクトル~ＹK|Kとを用いて、再帰的手法によって第１の次元ｎiに対応するシステム行列Ａd,ni,Ｂd,ni,Ｃd,ni,Ｄd,niを同定する。
　次にシステム特性推定部３２は、システム次元の探索範囲ｎi=(ｎ1,ｎ2,…,ｎa)（ただし、ｎ1<ｎ2<…<ｎa）に属するそれぞれの次元に対して、再帰的システム行列推定部３１から出力されるシステム行列Ａd,ni,Ｂd,ni,Ｃd,ni,Ｄd,niに基づいて、同定した線形離散時間システムに対して、実際の同定用入力データｕid(jＴS)(ｊ=0,1,2,…)を適用した場合のシステム出力を算出する。
　なお、再帰的システム行列推定部３１およびシステム特性推定部３２の処理は、ｉをインクリメントし、ｉ＝ａとなるまで実行される。
　システム次元推定部３３は、システム特性推定部３２から出力される線形離散時間システムのシステム出力と、動的システムの実際の同定用出力データｙid(jＴS)(ｊ=0,1,2,…)（図４では動的システムのシステム特性と記載）との時間領域における誤差２乗和ｅni(i=1,2,…,ａ)を算出し、図５に示したように誤差２乗和のノルム||ｅni||の分布４１が、予め設定された誤差２乗和ノルム閾値４２以下となる次元のうち、最小の次元をシステム次元ｎとして決定して出力する構成としている（図５の場合はシステム次元ｎ＝ｎ6）。
　次に、実施の形態１に係るシステム同定装置の動作について説明する。
　同定対象とする動的システムが、次式のように、１入力Ｐ出力のｎ次元線形離散時間システムで記述できるとする。
Figure JPOXMLDOC01-appb-M000001
　上記の動的システムへのシステム入力ｕ(jＴS)を擬似ランダム入力で構成すると、当該システム入力ｕ(jＴS)および、上記［数１］式に対応するシステム出力ｙ(jＴS)は、例えば図２に示したシステム入力１１、およびシステム出力１２のような時間波形となる。
　ここで、図１および図２の説明で述べたように、予め設定された比率閾値とシステム入力１１（ｕ(jＴS)）の最大値とを乗じた次式をシステム入力閾値１３として使用する。
Figure JPOXMLDOC01-appb-M000002
　システム入出力抽出部１は、システム入力１１の絶対値がシステム入力閾値１３以上となる時刻の最小値を擬似ランダム入力印加時刻ｊminＴSとして特定する（図２の場合はｊminＴS＝３ＴS）。
　また、システム入出力抽出部１は、擬似ランダム入力印加時刻ｊminＴS以降のシステム入力１１およびシステム出力１２を次式を用いて抽出する。
Figure JPOXMLDOC01-appb-M000003
　さらに、システム入出力抽出部１は、上記［数３］を用いて抽出したそれぞれの値を同定用入力データｕid(ｊＴS)および同定用出力データｙid(ｊＴS)とすることで、対象とする動的システムのシステム入出力から擬似ランダム入力印加前のシステム静止時間領域データを除去する。
　ブロックハンケル行列生成部２は、システム入出力抽出部１から出力される同定用入力データｕid(jＴS)(ｊ=0,1,2,…)および同定用出力データｙid(jＴS)(ｊ=0,1,2,…)に基づいて、次式で与えられるブロックハンケル行列Ｕp,Ｕf,Ｙp,Ｙfを生成する。
Figure JPOXMLDOC01-appb-M000004
　入出力ベクトル生成部３は、このブロックハンケル行列Ｕp,Ｕf,Ｙp,Ｙfに基づいて、次式で与えられる動的システムの入力ベクトル~ＵK|Kおよび出力ベクトル~ＹK|Kを生成する。
Figure JPOXMLDOC01-appb-M000005
　ＬＱ分解部４は、ブロックハンケル行列Ｕp,Ｕf,Ｙp,Ｙfを結合した次式で与えられるデータ行列を生成する。
Figure JPOXMLDOC01-appb-M000006
　また、ＬＱ分解部４は、上記データ行列を次式のようにＬＱ分解すると共に、ＬＱ分解した行列の要素から、部分行列Ｌ22,Ｌ32を出力する。
Figure JPOXMLDOC01-appb-M000007
　平行射影生成部５は、ＬＱ分解部４から出力される部分行列Ｌ22,Ｌ32と、ブロックハンケル行列生成部２から出力されるブロックハンケル行列Ｕp,Ｙpに基づいて、次式で定義される動的システムの平行射影Θを生成する。
Figure JPOXMLDOC01-appb-M000008
　特異値分解部６は、上式で示される平行射影Θを特異値分解することで、次式で与えられる平行射影Θの左特異ベクトルｕjを列ベクトルとする第１の直交行列Ｕ、平行射影Θの右特異ベクトルｖjを列ベクトルとする第２の直交行列Ｖおよび、平行射影Θの特異値σi(i=1,2,3…)を出力する。
Figure JPOXMLDOC01-appb-M000009
　対象とする動的システムのシステム次元ｎは、平行射影Θの特異値において有意な値を持つものはｎ個であり、それらに比してｎ＋１個目以降は十分小さいとする以下の関係に基づいて決定することができる。
Figure JPOXMLDOC01-appb-M000010
　図３に示したように、動的システムのシステム入出力から算出した平行射影Θの特異値σiは、理想的には次元(i=1,2,3…)に対して、例えば特異値分布２１に示す関係となる。この場合、有意な値を持つ特異値の個数を明確に規定することができ、当該個数から動的システムのシステム次元ｎを決定することができる（図３の場合はシステム次元ｎ＝４）。一方、観測ノイズ等の影響を受ける現実のシステム入出力に基づいて算出した特異値σiは、次元(i=1,2,3…)に対して、例えば特異値分布２２に示す関係となるため、有意な値を持つ特異値と、無視できる微小な値となる特異値との境界σn>>σn+1が不明確となる。したがって、従来の手法では、常に最適なシステム次元ｎを決定しているとは限らず、最適なシステム次元ｎを決定するためには、試行錯誤が必要になるという課題があった。
　そこで、実施の形態１に係るシステム同定装置１０では、“時間領域において実際のシステム入出力に最も適合する”という前提の下、システム次元決定部７において最適なシステム次元ｎを決定する。システム次元決定部７では、図１に示したように特異値分解部６から出力される第２の直交行列Ｖおよび特異値σi(i=1,2,3…)と、入出力ベクトル生成部３から出力される動的システムの入力ベクトル~ＵK|Kおよび出力ベクトル~ＹK|Kならびに、作業者が指定したシステム次元の探索範囲ｎi=(ｎ1,ｎ2,…,ｎa)（ただし、ｎ1<ｎ2<…<ｎa）に基づいて、当該探索範囲に属する次元ｎi(i=1,2,…,a)のそれぞれに対して動的システムを記述する線形離散時間システムのシステム行列を同定する。さらに当該システム行列に基づいて、当該探索範囲に属する次元ｎi(i=1,2,…,a)のそれぞれに対応する線形離散時間システムに対して、実際の同定用入力データｕid(jＴS)(ｊ=0,1,2,…)を適用した場合のシステム出力を算出し、動的システムの実際の同定用出力データｙid(jＴS)(ｊ=0,1,2,…）（図１では動的システムのシステム特性と記載）との比較から、システム次元ｎを決定する。
　具体的には図４に示したように、再帰的システム行列推定部３１は、作業者が指定したシステム次元の探索範囲ｎi(ｎ1,ｎ2,…,ｎa)（ただし、ｎ1<ｎ2<…<ｎa）に属する第１の次元ｎiに対応するシステム行列の同定に関して、第１の次元ｎiよりも１段階低い第２の次元ｎi-1に対応したシステム行列Ａd,ni-1,Ｂd,ni-1,Ｃd,ni-1,Ｄd,ni-1の同定結果と、特異値分解部６から出力される第２の直交行列Ｖおよび特異値σi(i=1,2,3…)のうち、それぞれ第２の次元ｎi-1よりも大きく、第１の次元ｎi以下となる右特異ベクトルｖjおよび特異値σj(j=ｎi-1+1,　ｎi-1+2,…,ｎi)と、入出力ベクトル生成部３から出力される動的システムの入力ベクトル~ＵK|Kおよび出力ベクトル~ＹK|Kとを用いて、次式に示す再帰的手法によって第１の次元ｎiに対応するシステム行列Ａd,ni,Ｂd,ni,Ｃd,ni,Ｄd,niを同定する。
Figure JPOXMLDOC01-appb-M000011
　システム特性推定部３２は、システム次元の探索範囲ｎi=(ｎ1,ｎ2,…,ｎa)（ただし、ｎ1<ｎ2<…<ｎa）に属する次元ｎiのそれぞれに対して、再帰的システム行列推定部３１から出力されるシステム行列Ａd,ni,Ｂd,ni,Ｃd,ni,Ｄd,niに基づいて、同定した線形離散時間システムに対して実際の同定用入力データｕid(jＴS)(ｊ=0,1,2,…)（［数３］参照）を適用した場合のシステム出力^ｙid,ni(jＴS)(ｊ=0,1,2,…）を算出する。なお、「^ｙ」の表記は、「ｙ」という文字の上部に「^」の記号を付したことを意味する代替表記である。
　また、システム次元推定部３３は、システム特性推定部３２から出力される線形離散時間システムのシステム出力^ｙid,ni(jＴS)(ｊ=0,1,2,…）と、動的システムの実際の同定用出力データｙid(jＴS)(ｊ=0,1,2,…）（図４では動的システムのシステム特性と記載）との時間領域における誤差２乗和ｅni(i=1,2,…,a）を次式を用いて算出する。
Figure JPOXMLDOC01-appb-M000012
　上式で示される誤差２乗和のノルム||ｅni||を最小とする次元ｎiが、“時間領域において実際のシステム入出力に最も適合する”システム次元ｎとなる。一方、実際のノルム||ｅni||は、観測ノイズが白色性であれば、そのノイズレベルによらず次元ｎiの増加に伴って単調減少し、図５に示すようにある次元以上でほぼ一定値となる。そこでここでは、システム次元ｎの推定値が必要以上に高次元となることを避けるため、次式で与えられる誤差２乗和ノルム閾値４２を規定する。
Figure JPOXMLDOC01-appb-M000013
　システム次元推定部３３は、誤差２乗和のノルム||ｅni||の分布４１が上記誤差２乗和ノルム閾値４２以下となる次元のうち、最小の次元をシステム次元ｎとして決定して出力する（図５の例では、システム次元ｎ＝ｎ6）。
　状態ベクトル生成部８は、特異値分解部６から出力される第２の直交行列Ｖおよび特異値σi(i=1,2,3…)と、システム次元決定部７から出力されるシステム次元ｎに基づいて、動的システムの状態ベクトル~ＸK,~ＸK+1を次式に従って生成する。
Figure JPOXMLDOC01-appb-M000014
　最後に、システム行列同定部９は、入出力ベクトル生成部３から出力される動的システムの入力ベクトル~ＵK|Kおよび出力ベクトル~ＹK|Kならびに、状態ベクトル生成部８から出力される動的システムの状態ベクトル~ＸK,~ＸK+1に基づいて、動的システムを記述する線形離散時間システムのシステム行列Ａd,Ｂd,Ｃd,Ｄdを次式を用いて同定して出力する。
Figure JPOXMLDOC01-appb-M000015
　このように、実施の形態１に係るシステム同定装置１０によれば、現実のシステム入出力から算出した平行射影Θの特異値σi(i=1,2,3…)がなだらかな単調減少となり、したがって有意な値を持つ特異値と、同定において無視できる微小な値となる特異値との境界が不明確となる場合においても、システム次元ｎの決定から試行錯誤を排除し、現実の動的システムに対して時間領域で一致度の高いシステム次元ｎを決定することができ、動的システムを記述する線形離散時間システムの同定が可能となる。
　また、動的システムの現実のシステム入出力から、擬似ランダム入力印加前のシステム静止時間領域データを除去することで、同定精度を向上させることも可能となる。
　さらに、再帰的システム行列推定部３１の存在により、現実の動的システムに対して一致度の高いシステム次元ｎを決定するための演算量を低減することが可能となる。
　なお、実施の形態１のシステム同定装置１０では、線形離散時間システムに対して実際の同定用入力データを適用した場合のシステム出力をシステム特性として算出し、当該システム出力と、動的システムの実際の同定用出力データとの時間領域における誤差２乗和のノルム分布４１が閾値４２以下となる次元のうち、最小の次元をシステム次元ｎとして決定している。しかしながら、本発明はこれに限定するものではなく、線形離散時間システムのシステム特性を周波数応答として算出し、当該周波数応答と、動的システムの同定用入出力データから得られる実際の周波数応答との周波数領域における誤差２乗和に基づいて、システム次元ｎを決定してもよい。この場合、さらに動的システムの実際の周波数応答に基づいて重み関数を決定し、線形離散時間システムの周波数応答と、動的システムの実際の周波数応答との周波数領域における誤差２乗値に当該重み関数を乗じた値の加算値に基づいて、システム次元ｎを決定してもよい。
実施の形態２．
　つぎに、実施の形態２に係るシステム同定装置について説明する。なお、実施の形態２に係るシステム同定装置の全体構成を示すブロック線図、平行射影Θの特異値σiと次元(i=1,2,3…)の関係を示す概略図および、同定した線形離散時間システムの周波数応答と動的システムの実際の周波数応答との周波数領域における誤差２乗和のノルム||ｅni||と次元ｎi(i=1,2,…,a)との関係を示す概略図は、それぞれ実施の形態１の説明において使用した図１、図３および図５と同一である。
　図６は、実施の形態２のシステム同定装置における動的システムをＭ系列加振した場合のシステム入出力の時間波形を示す概略図である。
　図６に示したように、実施の形態２に係るシステム同定装置１０では、同定対象とする動的システムに対してＭ系列信号を印加した場合のシステム入力１１（ｕ(jＴS)(ｊ=0,1,2,…)）およびシステム出力１２（ｙ(jＴS)(ｊ=0,1,2,…））に基づいて、当該システムを記述する線形離散時間システムを同定する。
　図７は、実施の形態２のシステム同定装置におけるシステム次元決定部７の内部構成を示すブロック線図である。図７において、図４と同一符号を付したものは、実施の形態１と同一または同等の構成要素であり、システム安定性評価部３４が追加されている。
　実施の形態２によるシステム次元決定部７では、図７に示すように、作業者が指定したシステム次元の探索範囲ｎi=(ｎ1,ｎ2,…,ｎa)（ただし、ｎ1<ｎ2<…<ｎa）に属する次元ｎiのそれぞれに対して、再帰的システム行列推定部３１で同定したシステム行列Ａd,ni,Ｂd,ni,Ｃd,ni,Ｄd,niに基づき、システム安定性評価部３４にて線形離散時間システムの安定性を評価する。
　システム特性推定部３２は、システム安定性評価部３４にて安定システムと判断された次元に対して、再帰的システム行列推定部３１から出力されるシステム行列Ａd,ni,Ｂd,ni,Ｃd,ni,Ｄd,niに基づいて、同定した線形離散時間システムに関する周波数応答を算出する。
　システム次元推定部３３は、動的システムのシステム入出力から得られる実際の周波数応答（図７では動的システムのシステム特性と記載）に基づいて重み関数を決定し、システム特性推定部３２から出力される線形離散時間システムの周波数応答と、動的システムの実際の周波数応答との周波数領域における誤差２乗値に対して、当該重み関数を乗じた値の加算値ｅni（ｎi：安定システムとなる次元）を算出し、図５に示したように当該加算値のノルム||ｅni||の分布４１が予め設定された誤差２乗和ノルム閾値４２以下となる次元のうち、最小の次元をシステム次元ｎとして決定して出力する（図５の場合はシステム次元ｎ＝ｎ6）。
　次に、実施の形態２に係るシステム同定装置の動作について説明する。
　同定対象とする動的システムが、１入力Ｐ出力のｎ次元線形離散時間システムとして［数１］で記述できるとする。この動的システムへのシステム入力ｕ(jＴS)をＭ系列信号で構成すると、当該システム入力ｕ(jＴS)および［数１］に対応するシステム出力ｙ(jＴS)は、例えば図６に示したシステム入力１１およびシステム出力１２のような時間波形となる。
　実施の形態２に係るシステム同定装置１０では、図１および図６に示すように、システム入出力抽出部１は、予め設定された比率閾値とシステム入力１１（ｕ(jＴS)）の最大値とを乗じた［数２］をシステム入力閾値１３とし、システム入力１１の絶対値がシステム入力閾値１３以上となる時刻の最小値をＭ系列信号印加時刻ｊminＴsとして特定する（図６の例はｊminＴs＝２Ｔs）。
　また、システム入出力抽出部１は、Ｍ系列信号印加時刻ｊminＴs以降のシステム入力１１およびシステム出力１２を［数３］で抽出し、抽出したそれぞれを同定用入力データｕid(jＴS)および同定用出力データｙid(jＴS)とすることで、対象とする動的システムのシステム入出力からＭ系列信号印加前のシステム静止時間領域データを除去する。
　次に実施の形態１と同様に、ブロックハンケル行列生成部２は［数４］で与えられるブロックハンケル行列Ｕp,Ｕf,Ｙp,Ｙfを生成し、入出力ベクトル生成部３は［数５］で与えられる動的システムの入力ベクトル~ＵK|Kおよび出力ベクトル~ＹK|Kを生成し、ＬＱ分解部４はブロックハンケル行列Ｕp,Ｕf,Ｙp,Ｙfを結合したデータ行列（［数６］）をＬＱ分解し（［数７］）、部分行列Ｌ22,Ｌ32を出力する。
　平行射影生成部５は、［数８］で定義される動的システムの平行射影Θを生成し、特異値分解部６は、生成された平行射影Θを特異値分解することで、［数９］で与えられる第１の直交行列Ｕ、第２の直交行列Ｖおよび、特異値σi(i=1,2,3…)を出力する。
　システム次元決定部７では、図７に示す処理が実行される。まず、再帰的システム行列推定部３１は、作業者が指定したシステム次元の探索範囲ｎi=(ｎ1,ｎ2,…,ｎa)（ただし、ｎ1<ｎ2<…<ｎa）に属する次元ｎiのそれぞれに関して、［数１１］に示される再帰的手法によって、対応するシステム行列Ａd,ni,Ｂd,ni,Ｃd,ni,Ｄd,niを同定する。
　次にシステム安定性評価部３４は、作業者が指定したシステム次元の探索範囲ｎi=(ｎ1,ｎ2,…,ｎa)（ただし、ｎ1<ｎ2<…<ｎa）に属する次元ｎiのそれぞれに対して、再帰的システム行列推定部３１が同定したシステム行列Ａd,niに基づいて、線形離散時間システムの安定性を、以下の内容にて評価する。
Figure JPOXMLDOC01-appb-M000016
　システム特性推定部３２は、システム安定性評価部３４で安定システムと判断された次元に対して、再帰的システム行列推定部３１が生成したシステム行列Ａd,ni,Ｂd,ni,Ｃd,ni,Ｄd,niに基づいて、同定された線形離散時間システムの周波数応答^Ｈni(kΔf)(k=0,1,2,…,N/2-1）を算出する。
　実施の形態２のシステム同定装置１０では、“周波数領域において実際の周波数応答に最も適合する”という前提の下、システム次元推定部３３にて最適なシステム次元ｎを決定する。具体的には、以下の通りである。
　まず、次式で与えられる同定用入出力データｕid(jＴS)，ｙid(jＴS)の有限離散化フーリエ変換Ｕid(kΔf),Ｙid(kΔf)(k=0,1,2,…,N/2-1）から、次々式で得られる動的システムの実際の周波数応答Ｈ(kΔf)(k=0,1,2,…,N/2-1）（図７では動的システムのシステム特性と記載）を算出する。
Figure JPOXMLDOC01-appb-M000017
Figure JPOXMLDOC01-appb-M000018
　次に周波数応答Ｈ(kΔf)(k=0,1,2,…,N/2-1）に基づいて、例えば高ゲインかつ低周波領域に重みを持たせた次式で示されるような重み関数Ｗ(kΔf)(k=0,1,2,…,N/2-1）を決定する。
Figure JPOXMLDOC01-appb-M000019
　そして、システム特性推定部３２から出力される線形離散時間システムの周波数応答^Ｈni(kΔf)と、動的システムの実際の周波数応答Ｈ(kΔf)との周波数領域における誤差２乗値に対して、重み関数Ｗ(kΔf)を乗じた値の加算値ｅni（ｎi：安定システムとなる次元）を次式にて算出する。
Figure JPOXMLDOC01-appb-M000020
　この重み付き誤差２乗和のノルム||ｅni||を最小とする次元ｎiが、“重み関数に応じて、周波数領域において実際の周波数応答に最も適合する”安定なシステム次元ｎとなる。ここでは、図５に示すように重み付き誤差２乗和のノルム||ｅni||の分布４１が、［数１３］で与えられる誤差２乗和ノルム閾値４２以下となる次元のうち、最小の次元をシステム次元ｎとして決定して出力する（図５の例は、システム次元ｎ＝ｎ6）。
　状態ベクトル生成部８は、特異値分解部６から出力される第２の直交行列Ｖおよび特異値σi(i=1,2,3…)と、システム次元決定部７から出力されるシステム次元ｎに基づいて、動的システムの状態ベクトル~ＸK,~ＸK+1を［数１４］を用いて生成する。
　最後にシステム行列同定部９は、入出力ベクトル生成部３から出力される動的システムの入力ベクトル~ＵK|Kおよび出力ベクトル~ＹK|Kならびに、状態ベクトル生成部８から出力される動的システムの状態ベクトル~ＸK,~ＸK+1に基づいて、動的システムを記述する線形離散時間システムのシステム行列Ａd,Ｂd,Ｃd,Ｄdを［数１５］を用いて同定して出力する。
　このように、実施の形態２に係るシステム同定装置１０によれば、現実のシステム入出力から算出した平行射影Θの特異値σi(i=1,2,3…)がなだらかな単調減少となり、したがって有意な値を持つ特異値と、同定において無視できる微小な値となる特異値との境界が不明確となる場合においても、システム次元ｎの決定から試行錯誤を排除し、現実の動的システムに対して、周波数領域での重み関数に応じて一致度の高いシステム次元ｎを決定することができ、動的システムを記述する線形離散時間システムの同定が可能となる。
　また、動的システムの現実のシステム入出力から、Ｍ系列信号印加前のシステム静止時間領域データを除去することで、同定精度を向上させることも可能となる。
　さらに、再帰的システム行列推定部３１の存在により、現実の動的システムに対して一致度の高いシステム次元ｎを決定するための演算量を低減することが可能となる。
　加えて、システム安定性評価部３４の存在により、現実の動的システムが安定システムであることが明確である場合に、安定システムに限定した線形離散時間システムの同定が可能となる。
　なお、実施の形態２のシステム同定装置１０では、線形離散時間システムのシステム特性を周波数応答として算出し、当該周波数応答と、動的システムの同定用入出力データから得られる実際の周波数応答との周波数領域における誤差２乗和のノルム分布４１が予め設定された閾値４２以下となる次元のうち、最小の次元をシステム次元ｎとして決定している。しかしながら、本発明はこれに限定するものではなく、線形離散時間システムに対して実際の同定用入力データを適用した場合のシステム出力をシステム特性として算出し、当該システム出力と、動的システムの実際の同定用出力データとの時間領域における誤差２乗和に基づいて、システム次元ｎを決定してもよい。
実施の形態３．
　実施の形態３では、同定対象とする動的システムがＤＣサーボモータである場合について説明する。図８は、実施の形態３に係る全体構成を示すブロック線図である。本実施の形態において、図８に示すシステム同定装置１０は、図１に示した実施の形態１に係るシステム同定装置１０と同一または同等の構成をとる。本実施の形態では、ＤＣサーボモータ５１の入力電流[Ａｒｍｓ]として、例えばＭ系列信号などの疑似ランダム信号を入力し、当該疑似ランダム信号を動的システムのシステム入力１１（ｕ(jＴS)(ｊ=0,1,2,…)）とする。また、角速度[rad/s]を動的システムのシステム出力１２（ｙ(jＴS)(ｊ=0,1,2,…)）として取得する。システム同定装置１０は、これらのシステム入出力および、システム次元の探索範囲を入力とし、ＤＣサーボモータ５１を記述する線形離散時間システムを同定する。このとき、システム次元の探索範囲は、ｎi＝(1,2,…,50)のように、予測されるシステム次元に対して十分な幅を取って設定すればよい。システム同定装置１０では、現実の動的システムに対して一致度の高いシステム次元の決定と、動的システムを記述する線形離散時間システムの同定とを可能とするため、当該線形離散時間システムをサーボモータの制御系におけるパラメータ設計、フィルタのパラメータ設計などに利用することができる。
　１　システム入出力抽出部、２　ブロックハンケル行列生成部、３　入出力ベクトル生成部、４　ＬＱ分解部、５　平行射影生成部、６　特異値分解部、７　システム次元決定部、８　状態ベクトル生成部、９　システム行列同定部、１０　システム同定装置、１１　システム入力、１２　システム出力、１３　システム入力閾値、２１　（理想的なシステム入出力における平行射影の）特異値分布、２２　（現実のシステム入出力における平行射影の）特異値分布、３１　再帰的システム行列推定部、３２　システム特性推定部、３３　システム次元推定部、３４　システム安定性評価部、４１　（時間領域または周波数領域における）誤差２乗和のノルム分布、４２　（時間領域または周波数領域における）誤差２乗和ノルム閾値、５１　ＤＣサーボモータ。

Claims (7)

  1. 　同定対象とする動的システムに対して擬似ランダム入力を印加した場合のシステム入出力および指定されたシステム次元の探索範囲を入力とするシステム同定装置において、
    　前記動的システムのシステム入出力から同定に適用する同定用入出力データを抽出するシステム入出力抽出部と、
    　前記同定用入出力データに基づいて、ブロックハンケル行列を生成するブロックハンケル行列生成部と、
    　前記ブロックハンケル行列に基づいて、前記動的システムの入力ベクトルおよび出力ベクトルを生成する入出力ベクトル生成部と、
    　前記ブロックハンケル行列を結合してデータ行列を生成し、当該データ行列をＬＱ分解した部分行列を出力するＬＱ分解部と、
    　前記部分行列と前記ブロックハンケル行列に基づいて、平行射影を生成する平行射影生成部と、
    　前記平行射影の特異値分解により、前記平行射影の左特異ベクトルを列ベクトルとする第１の直交行列、当該平行射影の右特異ベクトルを列ベクトルとする第２の直交行列および当該平行射影の特異値を出力する特異値分解部と、
    　前記第２の直交行列および前記特異値と、前記動的システムの入力ベクトルおよび出力ベクトルと、前記探索範囲に基づいて、当該探索範囲に属するそれぞれの次元に対して動的システムを記述する線形離散時間システムのシステム行列を同定し、さらに当該システム行列に基づいて算出した線形離散時間システムのシステム特性と、動的システムの実際のシステム特性との比較から、システム次元を決定するシステム次元決定部と、
    　前記第２の直交行列および特異値と、前記決定されたシステム次元に基づいて、前記動的システムの状態ベクトルを生成する状態ベクトル生成部と、
    　前記動的システムの入力ベクトルおよび出力ベクトルならびに前記動的システムの状態ベクトルに基づいて、当該動的システムを記述する線形離散時間システムのシステム行列を同定するシステム行列同定部と、を備え、
    　前記同定されたシステム行列を、前記動的システムを記述する線形離散時間システムとして出力することを特徴とするシステム同定装置。
  2. 　前記システム次元決定部は、
    　前記探索範囲に属するそれぞれの次元に対して、
    　同定した線形離散時間システムに関して実際の同定用入力データを適用した場合のシステム出力を算出し、当該システム出力を線形離散時間システムのシステム特性として出力するシステム特性推定部と、
    　前記線形離散時間システムのシステム出力と、前記動的システムの実際の同定用出力データとの時間領域における誤差２乗和のノルムが設定閾値以下となる次元のうち、最小の次元をシステム次元として決定して出力するシステム次元推定部と
    　を備えたことを特徴とする請求項１に記載のシステム同定装置。
  3. 　前記システム次元決定部が、
    　前記探索範囲に属するそれぞれの次元に対して、
    　同定した線形離散時間システムの周波数応答を算出し、当該周波数応答を線形離散時間システムのシステム特性として出力するシステム特性推定部と、
    　前記線形離散時間システムの周波数応答と、前記動的システムのシステム入出力から得られる実際の周波数応答との周波数領域における誤差２乗和のノルムが設定閾値以下となる次元のうち、最小の次元をシステム次元として決定して出力するシステム次元推定部と、
    　を備えたことを特徴とする請求項１に記載のシステム同定装置。
  4. 　前記システム次元推定部は、
    　前記動的システムのシステム入出力から得られる実際の周波数応答に基づいて重み関数を決定し、
    　前記線形離散時間システムの周波数応答と、前記動的システムの実際の周波数応答との周波数領域における誤差２乗値に対して、当該重み関数を乗じた値の加算値を算出し、
    　当該加算値のノルムが設定閾値以下となる次元のうち、最小の次元をシステム次元として決定して出力する
    　ことを特徴とする請求項３に記載のシステム同定装置。
  5. 　前記システム次元決定部は、
    　前記探索範囲に属する第１の次元に対応するシステム行列の同定に関して、
    　当該探索範囲の中で第１の次元よりも１段階低い第２の次元に対応したシステム行列の同定結果と、
    　前記第２の直交行列および特異値のうち、それぞれ前記第２の次元よりも大きく、前記第１の次元以下となる右特異ベクトルおよび特異値、ならびに、前記動的システムの入力ベクトルおよび出力ベクトルを用いて、再帰的手法により前記第１の次元に対応するシステム行列を同定する再帰的システム行列推定部と、
    　を備えたことを特徴とする請求項１に記載のシステム同定装置。
  6. 　前記システム入出力抽出部は、
    　設定された比率閾値とシステム入力の最大値とを乗じた値をシステム入力閾値とし、システム入力の絶対値がシステム入力閾値以上となる時刻の最小値を擬似ランダム入力印加時刻として、擬似ランダム入力印加時刻以降のシステム入力およびシステム出力を、それぞれ同定用入力データおよび同定用出力データとして抽出することを特徴とする請求項１に記載のシステム同定装置。
  7. 　前記システム次元決定部は、
    　前記探索範囲に属するそれぞれの次元に対して線形離散時間システムの安定性を評価するシステム安定性評価部を備え、
    　安定システムとなる次元に対応する線形離散時間システムのシステム特性から、システム次元を決定することを特徴とする請求項１に記載のシステム同定装置。
PCT/JP2014/079257 2014-02-07 2014-11-04 システム同定装置 WO2015118736A1 (ja)

Priority Applications (4)

Application Number Priority Date Filing Date Title
DE112014006135.2T DE112014006135T5 (de) 2014-02-07 2014-11-04 Systemidentifikationsvorrichtung
CN201480074855.7A CN105960614B (zh) 2014-02-07 2014-11-04 系统辨识装置
US15/114,729 US20160342731A1 (en) 2014-02-07 2014-11-04 System identification device
JP2015561158A JP6009105B2 (ja) 2014-02-07 2014-11-04 システム同定装置

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
JP2014022814 2014-02-07
JP2014-022814 2014-02-07

Publications (1)

Publication Number Publication Date
WO2015118736A1 true WO2015118736A1 (ja) 2015-08-13

Family

ID=53777555

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/JP2014/079257 WO2015118736A1 (ja) 2014-02-07 2014-11-04 システム同定装置

Country Status (5)

Country Link
US (1) US20160342731A1 (ja)
JP (1) JP6009105B2 (ja)
CN (1) CN105960614B (ja)
DE (1) DE112014006135T5 (ja)
WO (1) WO2015118736A1 (ja)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2020121527A1 (ja) * 2018-12-14 2020-06-18 三菱電機株式会社 学習識別装置、学習識別方法、及び、学習識別プログラム
WO2020188776A1 (en) * 2019-03-19 2020-09-24 Nec Corporation System identification device, non-transitory computer readable medium, and system identification method

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPS61267102A (ja) * 1985-05-22 1986-11-26 Toshiba Corp プラント・モデリング装置
JPH0433102A (ja) * 1990-05-30 1992-02-04 Toshiba Corp モデル予測制御装置
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 モデル同定装置およびモデル同定プログラム

Family Cites Families (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP5196653B2 (ja) * 2006-04-14 2013-05-15 国立大学法人岩手大学 システム同定方法及びプログラム及び記憶媒体、システム同定装置
DE102009038011B9 (de) * 2009-08-20 2018-04-12 Schenck Rotec Gmbh Verfahren zur automatischen Erfassung und Erkennung von Fehlern an einer Auswuchtmaschine
CN102183699B (zh) * 2011-01-30 2012-12-26 浙江大学 化工过程多变量预测控制系统模型失配检测与定位方法
US9558300B2 (en) * 2011-11-11 2017-01-31 Carnegie Mellon University Stochastic computational model parameter synthesis system
CN103501150B (zh) * 2013-10-12 2017-01-25 上海联孚新能源科技集团有限公司 一种内嵌式永磁同步电机参数辨识装置及方法
CN104699894B (zh) * 2015-01-26 2017-07-28 江南大学 基于实时学习的高斯过程回归多模型融合建模方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPS61267102A (ja) * 1985-05-22 1986-11-26 Toshiba Corp プラント・モデリング装置
JPH0433102A (ja) * 1990-05-30 1992-02-04 Toshiba Corp モデル予測制御装置
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 モデル同定装置およびモデル同定プログラム

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
GABRIE ELKAIM: "System Identification of a Farm Vehicle Using Carrier-Phase Differential GPS", PROCEEDINGS OF THE 9TH INTERNATIONAL TECHNICAL MEETING OF THE SATELLITE DIVISION OF THE INSTITUDE OF NAVIGATION, 20 September 1996 (1996-09-20), pages 485 - 494, XP055219386 *
HIROSHI OKU: "Identification Experiment and Control Design of an Inverted Pendulum System via Closed-Loop Subspace Model Identification", JOURNAL OF THE SOCIETY OF INSTRUMENT AND CONTROL ENGINEERS, vol. 49, no. 7, 10 July 2010 (2010-07-10), pages 457 - 462 *
TAKASHI KIDA: "Attitude Control Technology of Large Spacecraft", JOURNAL OF THE SOCIETY OF INSTRUMENT AND CONTROL ENGINEERS, vol. 35, no. 11, 10 November 1996 (1996-11-10), pages 833 - 836, XP055219421, DOI: doi:10.11499/sicejl1962.35.11_833 *
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, 19 May 2004 (2004-05-19), pages 165 - 166 *

Also Published As

Publication number Publication date
US20160342731A1 (en) 2016-11-24
DE112014006135T5 (de) 2016-09-29
JPWO2015118736A1 (ja) 2017-03-23
JP6009105B2 (ja) 2016-10-19
CN105960614B (zh) 2020-11-27
CN105960614A (zh) 2016-09-21

Similar Documents

Publication Publication Date Title
Wang et al. Recursive least squares estimation algorithm applied to a class of linear-in-parameters output error moving average systems
Rajagopal et al. Dynamical analysis and FPGA implementation of a novel hyperchaotic system and its synchronization using adaptive sliding mode control and genetically optimized PID control
Hušková et al. Bootstrapping sequential change-point tests for linear regression
EP2870450A1 (en) Method for determining modal parameters
Akçay Subspace-based spectrum estimation in frequency-domain by regularized nuclear norm minimization
CN112861066A (zh) 基于机器学习和fft的盲源分离信源数目并行估计方法
JP6009105B2 (ja) システム同定装置
JP6009106B2 (ja) システム同定装置
Xia et al. A fast algorithm for AR parameter estimation using a novel noise-constrained least-squares method
WO2004090782A1 (en) Accurate linear parameter estimation with noisy inputs
Kay A new nonstationarity detector
Chen et al. Maximum likelihood based recursive parameter estimation for controlled autoregressive ARMA systems using the data filtering technique
Timoshenkova et al. On the possibility of the forecast correction for inaccurate observations based on data assimilation
Ding et al. Parameter and state estimator for state space models
Rahman et al. A unified analysis of proposed wavelet transform domain LMS-algorithm for ARMA process
Ogunfunmi et al. An alternative kernel adaptive filtering algorithm for quaternion-valued data
Soverini et al. Frequency domain identification of autoregressive models in the presence of additive noise
Nikolov et al. Discrete adaptive real-time state observer development using least-squares method
Ishihara Robust Kalman filtering for nonlinear systems with parameter uncertainties
Xia et al. Application of Kalman filter in microseismic data denoising based on identified signal model
JP2009098203A (ja) 信号推定装置、その方法、そのプログラム、その記録媒体
Ding Recursive least squares algorithm for parameter identification of multi-input output-error systems using the data filtering
JP6504381B2 (ja) システム同定装置及びそのプログラム
Costa et al. New insights into the weight behaviour of the Affine Projection algorithm
Ghanbari et al. Identification Of Flexible Robot Arm System Using Extended Volterra Series By Kautz Orthogonal Functions

Legal Events

Date Code Title Description
121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 14881990

Country of ref document: EP

Kind code of ref document: A1

ENP Entry into the national phase

Ref document number: 2015561158

Country of ref document: JP

Kind code of ref document: A

WWE Wipo information: entry into national phase

Ref document number: 15114729

Country of ref document: US

WWE Wipo information: entry into national phase

Ref document number: 112014006135

Country of ref document: DE

122 Ep: pct application non-entry in european phase

Ref document number: 14881990

Country of ref document: EP

Kind code of ref document: A1