JP6981526B2 - システム同定装置、システム同定方法及びコンピュータプログラム - Google Patents
システム同定装置、システム同定方法及びコンピュータプログラム Download PDFInfo
- Publication number
- JP6981526B2 JP6981526B2 JP2020501746A JP2020501746A JP6981526B2 JP 6981526 B2 JP6981526 B2 JP 6981526B2 JP 2020501746 A JP2020501746 A JP 2020501746A JP 2020501746 A JP2020501746 A JP 2020501746A JP 6981526 B2 JP6981526 B2 JP 6981526B2
- Authority
- JP
- Japan
- Prior art keywords
- unit
- response function
- system identification
- analysis
- analysis target
- 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.)
- Active
Links
- 238000000034 method Methods 0.000 title claims description 32
- 238000004590 computer program Methods 0.000 title claims 2
- 238000004458 analytical method Methods 0.000 claims description 79
- 238000005316 response function Methods 0.000 claims description 76
- 238000005259 measurement Methods 0.000 claims description 28
- 238000009434 installation Methods 0.000 claims description 15
- 238000013016 damping Methods 0.000 claims description 7
- 238000006073 displacement reaction Methods 0.000 claims description 7
- 230000001133 acceleration Effects 0.000 claims description 6
- 238000012545 processing Methods 0.000 description 26
- 238000010586 diagram Methods 0.000 description 12
- 230000006870 function Effects 0.000 description 8
- 238000002474 experimental method Methods 0.000 description 7
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 5
- 238000012360 testing method Methods 0.000 description 4
- 238000005516 engineering process Methods 0.000 description 3
- 230000005284 excitation Effects 0.000 description 3
- 239000011159 matrix material Substances 0.000 description 3
- 238000005070 sampling Methods 0.000 description 3
- 238000001228 spectrum Methods 0.000 description 3
- 230000004044 response Effects 0.000 description 2
- 238000010079 rubber tapping Methods 0.000 description 2
- 239000013598 vector Substances 0.000 description 2
- 241001123248 Arma Species 0.000 description 1
- 229910001018 Cast iron Inorganic materials 0.000 description 1
- 229910001141 Ductile iron Inorganic materials 0.000 description 1
- 238000004364 calculation method Methods 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 239000000284 extract Substances 0.000 description 1
- 238000012544 monitoring process Methods 0.000 description 1
- 230000003595 spectral effect Effects 0.000 description 1
- 230000001131 transforming effect Effects 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
- G01N29/00—Investigating or analysing materials by the use of ultrasonic, sonic or infrasonic waves; Visualisation of the interior of objects by transmitting ultrasonic or sonic waves through the object
- G01N29/44—Processing the detected response signal, e.g. electronic circuits specially adapted therefor
- G01N29/4472—Mathematical theories or simulation
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01H—MEASUREMENT OF MECHANICAL VIBRATIONS OR ULTRASONIC, SONIC OR INFRASONIC WAVES
- G01H13/00—Measuring resonant frequency
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01M—TESTING STATIC OR DYNAMIC BALANCE OF MACHINES OR STRUCTURES; TESTING OF STRUCTURES OR APPARATUS, NOT OTHERWISE PROVIDED FOR
- G01M7/00—Vibration-testing of structures; Shock-testing of structures
- G01M7/02—Vibration-testing by means of a shake table
- G01M7/025—Measuring arrangements
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01M—TESTING STATIC OR DYNAMIC BALANCE OF MACHINES OR STRUCTURES; TESTING OF STRUCTURES OR APPARATUS, NOT OTHERWISE PROVIDED FOR
- G01M7/00—Vibration-testing of structures; Shock-testing of structures
- G01M7/08—Shock-testing
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
- G01N29/00—Investigating or analysing materials by the use of ultrasonic, sonic or infrasonic waves; Visualisation of the interior of objects by transmitting ultrasonic or sonic waves through the object
- G01N29/04—Analysing solids
- G01N29/12—Analysing solids by measuring frequency or resonance of acoustic waves
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
- G01N29/00—Investigating or analysing materials by the use of ultrasonic, sonic or infrasonic waves; Visualisation of the interior of objects by transmitting ultrasonic or sonic waves through the object
- G01N29/44—Processing the detected response signal, e.g. electronic circuits specially adapted therefor
- G01N29/4409—Processing the detected response signal, e.g. electronic circuits specially adapted therefor by comparison
- G01N29/4418—Processing the detected response signal, e.g. electronic circuits specially adapted therefor by comparison with a model, e.g. best-fit, regression analysis
-
- 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
- G05B23/00—Testing or monitoring of control systems or parts thereof
- G05B23/02—Electric testing or monitoring
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
- G01N2291/00—Indexing codes associated with group G01N29/00
- G01N2291/26—Scanned objects
- G01N2291/262—Linear objects
- G01N2291/2626—Wires, bars, rods
Landscapes
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Pathology (AREA)
- Life Sciences & Earth Sciences (AREA)
- Chemical & Material Sciences (AREA)
- Analytical Chemistry (AREA)
- Biochemistry (AREA)
- General Health & Medical Sciences (AREA)
- Immunology (AREA)
- Health & Medical Sciences (AREA)
- Engineering & Computer Science (AREA)
- Signal Processing (AREA)
- Algebra (AREA)
- Mathematical Optimization (AREA)
- Mathematical Physics (AREA)
- Pure & Applied Mathematics (AREA)
- Mathematical Analysis (AREA)
- Acoustics & Sound (AREA)
- Automation & Control Theory (AREA)
- Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)
- Testing And Monitoring For Control Systems (AREA)
Description
本発明は、システム同定装置、システム同定方法及び記録媒体に関する。
石油、ガス、水道などのプラントシステムや、産業用ロボット等の物理システムを対象として、ＩｏＴ（Internet of Things）による監視や制御を行う際には、対象とするシステムのモデリング技術が重要である。このモデリング技術として、観測データに基づいて物理システムの数理的モデリングを行うものがある（例えば、特許文献１〜４、非特許文献１参照）。
中溝 高好、「信号解析とシステム同定」、コロナ社、ｐ．２２−２４、４９−５３，１２１−１２７、１９８８年
しかし、モデリングの対象システムが、うなり現象や共鳴が発生するような、近接固有値を有する系である場合、特許文献１〜４および非特許文献１に記載の技術では、システム同定が困難な場合があった。なお、近接固有値を有する系は、固有値が縮退した重根を有する系も含むものとする。
本発明は、近接固有値を有する系のシステム同定を行うことができるシステム同定装置等を提供することを目的としている。
本発明は、近接固有値を有する系のシステム同定を行うことができるシステム同定装置等を提供することを目的としている。
本発明の第１の態様によれば、システム同定装置は、分析対象を励振させた位置において計測された入力信号及び出力信号に基づいて自己周波数応答関数を算出し、算出された前記自己周波数応答関数から得られるインパルス応答関数と、前記分析対象がモデル化された仮想二自由度モデルのインパルス応答関数とを用いて前記分析対象のシステム同定を行う分析部と、を備える。
本発明の第２の態様によれば、システム同定方法は、分析対象を励振させる加振ステップと、前記加振ステップにおいて前記分析対象を励振させた位置における入力信号及び出力信号を計測する計測ステップと、前記計測ステップにおいて計測された前記入力信号及び前記出力信号に基づいて自己周波数応答関数を算出し、算出された前記自己周波数応答関数から得られるインパルス応答関数と、前記分析対象がモデル化された仮想二自由度モデルのインパルス応答関数とを用いて前記分析対象のシステム同定を行う分析ステップと、を有する。
本発明の第３の態様によれば、記録媒体は、コンピュータに、分析対象を励振させた位置において計測された入力信号及び出力信号に基づいて自己周波数応答関数を算出し、算出された前記自己周波数応答関数から得られるインパルス応答関数と、前記分析対象がモデル化された仮想二自由度モデルのインパルス応答関数とを用いて前記分析対象のシステム同定を行う分析ステップ、を実行させるプログラムを記録する。
本発明によれば、近接固有値を有する系のシステム同定を行うことができる。
以下、図面を参照して本発明の一実施の形態を説明する。
観測データに基づく物理システムの数理的なモデリング手法は、「システム同定問題」と呼ばれている。本問題では大別して、（１）システムの入力および出力の信号が既知の場合と、（２）入力が未知の場合とに分類される。さらには、時間領域信号又は周波数領域信号を用いた技術も知られている。
観測データに基づく物理システムの数理的なモデリング手法は、「システム同定問題」と呼ばれている。本問題では大別して、（１）システムの入力および出力の信号が既知の場合と、（２）入力が未知の場合とに分類される。さらには、時間領域信号又は周波数領域信号を用いた技術も知られている。
時間領域情報を用いた技術では、ＡＲモデル(Autoregressive model：自己回帰モデル)やＭＡモデル(Moving Average model：移動平均モデル)、ＡＲＭＡモデル(Autoregressive moving average model：自己回帰移動平均モデル)、ＡＲＸモデル（Auto-Regressive eXogeneous model）などの多項式モデルが用いられる。多項式モデルでは、周波数領域が平坦化され、異なる固有値が近接している近接固有値を有する系への適用が困難であった。一方で、周波数領域情報を用いた場合には、固有値双方のピーク位置が不明瞭となる場合が多く、曲線適合法を適用することができない。特に、フーリエ変換による周波数領域のサンプル数によっては、ピーク位置を視認することすらできない。
上記のことから、近接固有値を有する系では固有値が近接するため、周波数領域同定法および時間領域同定法では、システム同定が困難となる。このように、近接固有値を有する系のシステム同定問題は未だ十分に解決されていない問題である。本実施形態のシステム同定装置等は、このような問題を解決し、近接固有値を有する系（固有値が重なる重根を有する系も含む）のシステム同定を行う。
［第１の実施形態］
図１は、第１の実施形態によるシステム同定装置１の構成を示すブロック図である。システム同定装置１は、設置位置決め部１０１、加振部１０２、計測部１０３、信号収集部１０４及び分析部１０５を有する。対象物理システム１０６は、システム同定装置１による同定対象である。設置位置決め部１０１は、対象物理システム１０６に、加振部１０２及び計測部１０３を設置する。加振部１０２は、設置位置決め部１０１を介して対象物理システム１０６を励振させる。計測部１０３は、加振部１０２が設置位置決め部１０１を介して対象物理システム１０６を励振させたときの、対象物理システム１０６への入力信号と出力信号を検知する。信号収集部１０４は、計測部１０３が検知した入力信号及び出力信号をデータ化する。分析部１０５は、信号収集部１０４により得られたデータを分析し、対象物理システム１０６のシステム同定を行う。
図１は、第１の実施形態によるシステム同定装置１の構成を示すブロック図である。システム同定装置１は、設置位置決め部１０１、加振部１０２、計測部１０３、信号収集部１０４及び分析部１０５を有する。対象物理システム１０６は、システム同定装置１による同定対象である。設置位置決め部１０１は、対象物理システム１０６に、加振部１０２及び計測部１０３を設置する。加振部１０２は、設置位置決め部１０１を介して対象物理システム１０６を励振させる。計測部１０３は、加振部１０２が設置位置決め部１０１を介して対象物理システム１０６を励振させたときの、対象物理システム１０６への入力信号と出力信号を検知する。信号収集部１０４は、計測部１０３が検知した入力信号及び出力信号をデータ化する。分析部１０５は、信号収集部１０４により得られたデータを分析し、対象物理システム１０６のシステム同定を行う。
図２は、システム同定装置１の処理を示すフローチャートである。測定者は、対象物理システム１０６に対し、設置位置決め部１０１を介して加振部１０２及び計測部１０３を設置する（ステップＳ１１０）。このとき、入力位置と出力位置とを一致させるようにする。入力位置とは、加振部１０２が設置される、すなわち、振動の原因が入力される位置である。出力位置とは、計測部１０３が設定される、すなわち、対象物理システム１０６の振動を計測する位置である。
自己周波数応答関数を計測するため、加振部１０２により、設置位置決め部１０１を介して対象物理システム１０６を励振させる。計測部１０３は、対象物理システム１０６への振動の入力信号と対象物理システム１０６からの振動の出力信号を検知する。信号収集部１０４は、計測部１０３が検知した入力信号及び出力信号をデータ化し、分析部１０５に出力する。分析部１０５は、得られたデータを分析する。具体的には、分析部１０５は、最初に、入力信号と出力信号をそれぞれ高速フーリエ変換（ＦＦＴ：fast fourier transform）する。分析部１０５は、周波数領域で出力信号を入力信号により除算することにより、自己周波数応答関数を求める（ステップＳ１２０）。
次いで、分析部１０５は、自己周波数応答関数において、対象とする近接固有値が存在する周波数帯域のみズーミングを行う（ステップＳ１３０）。分析部１０５は、ズーミングした自己周波数応答関数に逆フーリエ変換（Inverse fourier transform）を施すことにより、自己周波数応答関数のインパルス応答関数を求める（ステップＳ１４０）。
分析部１０５は、次のステップＳ１６０において用いられる初期値及びステップ幅の入力を受ける（ステップＳ１５０）。分析部１０５は、ステップＳ１４０で得られたインパルス応答関数に対し、後述する仮想二自由度系のインパルス応答関数を用いた多変数ニュートン法を適用する（ステップＳ１６０）。分析部１０５は、多変数ニュートン法の実行時に、ステップＳ１５０において入力された初期値及びステップ幅を用いる。
分析部１０５は、解が収束しないと判定した場合（ステップＳ１７０：ＮＯ）、ステップＳ１５０に戻って、新たに使用する初期値及びステップ幅の入力を受け、ステップＳ１６０の処理を行う。分析部１０５は、解が収束したと判定した場合（ステップＳ１７０：ＹＥＳ）、その収束したときの仮想二自由度系のインパルス応答関数に基づいて対象物理システム１０６の系の質量、剛性定数及び減衰係数を求め、システム同定を行う（ステップＳ１８０）。
本実施形態によれば、近接固有値を有する系のシステム同定が可能である。
本実施形態によれば、近接固有値を有する系のシステム同定が可能である。
［第２の実施形態］
本実施形態は、第１の実施形態を管路のシステム同定に適用する。
図３は、第２の実施形態によるシステム同定装置３の構成を示すブロック図である。システム同定装置３は、設置位置決め部３０１、加振部３０２、計測部３０３、信号収集部３０４、分析部３０５、記憶部３０６及び初期値設定部３０７を有する。対象物理システム３０８は、管路であり、システム同定装置３による同定対象である。
本実施形態は、第１の実施形態を管路のシステム同定に適用する。
図３は、第２の実施形態によるシステム同定装置３の構成を示すブロック図である。システム同定装置３は、設置位置決め部３０１、加振部３０２、計測部３０３、信号収集部３０４、分析部３０５、記憶部３０６及び初期値設定部３０７を有する。対象物理システム３０８は、管路であり、システム同定装置３による同定対象である。
設置位置決め部３０１、加振部３０２、計測部３０３及び信号収集部３０４はそれぞれ、第１の実施形態の設置位置決め部１０１、加振部１０２、計測部１０３及び信号収集部１０４と同様の機能を有する。記憶部３０６は、管路管理台帳データを記憶する。管路管理台帳データは、対象物理システム３０８の物理データである口径、材種、管厚の情報を含む。初期値設定部３０７は、記憶部３０６から読み出した口径、材種、管厚のデータに基づいて、分析部３０５において用いられる多変数ニュートン法のパラメータの初期値を算出する。分析部３０５は、初期値設定部３０７が算出した初期値を多変数ニュートン法において用いる点以外は、第１の実施形態の分析部１０５と同様の機能を有する。
図４は、システム同定装置３の処理を示すフローチャートである。測定者は、対象物理システム３０８に対し、設置位置決め部３０１を介して加振部３０２及び計測部３０３を設置する（ステップＳ３１０）。加振部３０２により、設置位置決め部３０１を介して対象物理システム３０８を励振させる。計測部３０３は、加振位置における対象物理システム３０８への入力信号と対象物理システム３０８からの出力信号を検知する。信号収集部３０４は、計測部３０３が検知した入力信号及び出力信号をデータ化する。分析部３０５は、入力信号と出力信号を用いて自己周波数応答関数を求める（ステップＳ３２０）。分析部３０５は、自己周波数応答関数において、対象とする近接固有値が存在する周波数帯域のみズーミングを行う（ステップＳ３３０）。分析部３０５は、ズーミング部分を逆フーリエ変換し、自己周波数応答関数のインパルス応答関数を求める（ステップＳ３４０）。
初期値設定部３０７は、記憶部３０６から対象物理システム３０８の口径、材種、管厚のデータを読み出す（ステップＳ３５０）。初期値設定部３０７は、読み出したデータに基づいて多変数ニュートン法のパラメータの初期値を算出する（ステップＳ３６０）。分析部３０５は、ステップ幅の入力を受ける（ステップＳ３７０）。分析部３０５は、ステップＳ３４０において得られたインパルス応答関数に対し、後述する仮想二自由度系のインパルス応答関数を用いた多変数ニュートン法を適用する（ステップＳ３８０）。多変数ニュートン法においては、ステップＳ３６０において算出された初期値と、ステップＳ３７０において入力されたステップ幅とを用いる。
分析部３０５は、解が収束しないと判定した場合（ステップＳ３９０：ＮＯ）、ステップＳ３７０に戻って、新たに使用するステップ幅の入力を受け、ステップＳ３８０の処理を行う。分析部３０５は、解が収束したと判定した場合（ステップＳ３９０：ＹＥＳ）、その収束したときの仮想二自由度系のインパルス応答関数に基づいて対象物理システム３０８の系の質量、剛性定数及び減衰係数を求め、システム同定を行う（ステップＳ４００）。
本実施形態によれば、管路のシステム同定が可能である。
本実施形態によれば、管路のシステム同定が可能である。
［第３の実施形態］
ここでは、第１の実施形態を水道管路のシステム同定に適用する。
図５は、本実施形態のシステム同定装置５の構成を示すブロック図である。システム同定装置５は、消火栓カプラ５０１、ハンマー５０２、センサ５０３、データロガー５０４及び同定処理部５０５を有する。管路５０６は、システム同定装置５によるシステム同定対象の水道管路である。ハンマー５０２として、力センサを内蔵したインパルスハンマー、市販ハンマーに加速度ピックアップを設置したもの、電磁式加振器などが例示できる。センサ５０３として、加速度ピックアップ、レーザードップラー型速度計、レーザー変位計、接触式変位計などが例示できる。同定処理部５０５は、例えば、プロセッサ、メモリ及びＨＤＤ（Hard disk drive）により実現される。プロセッサは、処理をコンピュータに実行させるための同定処理プログラムをＨＤＤから読み出して実行することにより同定処理部５０５として動作する。
ここでは、第１の実施形態を水道管路のシステム同定に適用する。
図５は、本実施形態のシステム同定装置５の構成を示すブロック図である。システム同定装置５は、消火栓カプラ５０１、ハンマー５０２、センサ５０３、データロガー５０４及び同定処理部５０５を有する。管路５０６は、システム同定装置５によるシステム同定対象の水道管路である。ハンマー５０２として、力センサを内蔵したインパルスハンマー、市販ハンマーに加速度ピックアップを設置したもの、電磁式加振器などが例示できる。センサ５０３として、加速度ピックアップ、レーザードップラー型速度計、レーザー変位計、接触式変位計などが例示できる。同定処理部５０５は、例えば、プロセッサ、メモリ及びＨＤＤ（Hard disk drive）により実現される。プロセッサは、処理をコンピュータに実行させるための同定処理プログラムをＨＤＤから読み出して実行することにより同定処理部５０５として動作する。
図６は、管路５０６に設置された消火栓カプラ５０１の状況を示す図である。計測者は、管路５０６に消火栓カプラ５０１及びセンサ５０３を設置する（図２のステップＳ１１０）。計測者は、ハンマー５０２により図６に示す消火栓カプラ５０１を打撃し、管路５０６を励振する。センサ５０３は、打撃位置（Tapping Point）と同じ測定位置（Measurements point for vibration response）における励振後の出力信号を検出する。データロガー５０４は、ハンマー５０２の入力信号及びセンサ５０３の出力信号を収集する。データロガー５０４は、収集した入力信号及び出力信号をデータ化して同定処理部５０５に出力する。
同定処理部５０５は、入力信号と出力信号それぞれに高速フーリエ変換（ＦＦＴ）処理を行う。ＦＦＴにより得られた入力信号のスペクトルをＸ（ω）、出力信号のスペクトルをＹ（ω）とそれぞれ表す。ωは周波数である。同定処理部５０５は、各周波数領域のスペクトルを除し、自己周波数応答関数Ｌ（ω）＝Ｙ（ω）／Ｘ（ω）を求める（図２のステップＳ１２０）。ここで、自己周波数応答関数の算出にあたり、Ｌ（ω）＝（Ｙ（ω）・Ｘ＊（ω））／（Ｘ（ω）・Ｘ＊（ω））をＨ１推定、Ｌ（ω）＝（Ｙ（ω）・Ｙ＊（ω））／（Ｘ（ω）・Ｙ＊（ω））をＨ２推定と呼ぶが、どちらを用いても良い。なお、Ｘ＊（ω）はＸ（ω）の複素共役、Ｙ＊（ω）はＹ（ω）の複素共役である。
同定処理部５０５は、得られた自己周波数応答関数Ｌ（ω）において、着目する近接固有値が存在する周波数帯域のみをズーミングする（図２のステップＳ１３０）。着目する近接固有値が存在する周波数帯域にはピークが現れる。そこで、同定処理部５０５は、自己周波数応答関数Ｌ（ω）におけるピークを検出することにより、着目する近接固有値が存在する周波数帯域を特定する。あるいは、同定処理部５０５は、自己周波数応答関数Ｌ（ω）をシステム同定装置５が備える表示装置に表示し、表示を確認したユーザが、ピークが現れる周波数帯域を入力してもよい。同定処理部５０５は、特定された周波数帯域の自己周波数応答関数Ｌ（ω）を抽出し、閾値よりも小さい値については０に置き換えるズーミングを行う。同定処理部５０５は、ズーミングした結果を逆ＦＦＴすることにより、インパルス応答関数ｇｅ（ｔ）を得る。なお、ｔは時間を表す（図２のステップＳ１４０）。
ここで、システム同定用モデルとして仮想二自由度モデルを仮想する。なお、仮想二自由度モデルは、対称型二自由度ばね質量系とも呼ばれる。
図７は、仮想二自由度モデルを示す図である。仮想二自由度モデルは、等しい質量、ばね定数、及び、減衰係数からなる一自由度ばね質量系が、ばね及びダッシュポットによって接続された系である。Ｍは質量、Ｋはばね定数、Ｃは減衰係数、Ｆは外力ベクトル、ｘ１、ｘ２は変位ベクトルである。また、ΔＫはばね定数の変化、ΔＣは減衰係数の変化を示す。仮想二自由度モデルは、運動方程式を表す質量行列、剛性行列、減衰行列が対称行列となるシステムであり、固有ベクトルが対称となる性質を有しており、近接固有値を有する系のシステム同定に好適となる。
図７は、仮想二自由度モデルを示す図である。仮想二自由度モデルは、等しい質量、ばね定数、及び、減衰係数からなる一自由度ばね質量系が、ばね及びダッシュポットによって接続された系である。Ｍは質量、Ｋはばね定数、Ｃは減衰係数、Ｆは外力ベクトル、ｘ１、ｘ２は変位ベクトルである。また、ΔＫはばね定数の変化、ΔＣは減衰係数の変化を示す。仮想二自由度モデルは、運動方程式を表す質量行列、剛性行列、減衰行列が対称行列となるシステムであり、固有ベクトルが対称となる性質を有しており、近接固有値を有する系のシステム同定に好適となる。
図７に示す仮想二自由度モデルの自己周波数応答関数Ｌ１１を求めると、以下の式（１）となる。ｓは複素数を表す。ｓはｓ＝ｉ・ωで表される（ただし、ωは周波数、ｉは虚数単位）。
式（１）をラブラス逆変換し、仮想二自由度モデルのインパルス応答関数ｇ１１を求めると、式（２）となる。なお、ｓ＝ｉ・ωを式（１）に代入した場合、これを逆フーリエ変換すると、式（２）となる。
δ（ｔ）はディラックのデルタ関数である。第１項の（１／Ｍ）・δ（ｔ）は、無視できる程度に他の項よりも小さな値である。
なお、ここで、各パラメータωｄ１、ωｄ２、α、β、γを式（３）とした。
なお、ここで、各パラメータωｄ１、ωｄ２、α、β、γを式（３）とした。
これより、未知パラメータは計５つとなる。実験値のインパルス応答関数ｇｅ（ｔ）と式（２）の差の二乗和Ｊを目的関数とした多変数ニュートン法によって、パラメータ推定を行う更新式を求める。目的関数Ｊを式（４）に示し、更新式を式（５）に示す。ｉはサンプリングの番号、ｔｉはｉ番目のサンプリングを行った時間を表す。
ここで、λはステップ調整用パラメータであり、パラメータ推定アルゴリズムが発散する場合に０．００１〜０．１の間で調整すると、収束性が改善される。特に、管路系のシステム同定では、λを０．０１程度に設定することが好ましい。ｇ＾１１は、現在のパラメータ値を用いて式（２）を算出したものである。図２のステップＳ１５０においては、各パラメータ（ωｄ１、ωｄ２、α、β、γ）の初期値、及び、ステップ幅となるλの値が入力される。なお、同定処理部５０５は、初期値の算出に用いられる情報の入力を受け、入力された情報に基づいて初期値を算出してもよい。
同定処理部５０５は、ステップＳ１４０により得られたインパルス応答関数ｇｅ（ｔｉ）と、現在の各パラメータの値を用いて式（２）により算出したｇ＾１１（ｔｉ）とに基づいて、式（４）により二乗和Ｊを算出する。同定処理部５０５は、二乗和Ｊが閾値以下となるように、式（５）により各パラメータの値を更新する（図２のステップＳ１６０）。
同定処理部５０５は、パラメータの更新を繰り返し、Ｊの値やその値の変化によって収束したか否かを判断する。同定処理部５０５は、収束しないと判断した場合（図２のステップＳ１７０：ＮＯ）、新たな初期値及びλの入力を受ける（ステップＳ１５０）。同定処理部５０５は、収束したと判断した場合、そのときのパラメータ値を用い、式（５）の関係に基づいて、質量Ｍ、ばね定数Ｋ及び減衰係数Ｃを算出する（ステップＳ１８０）。
数値実験により、本実施形態によるシステム同定法の動作検証をおこなった。口径１００ｍｍのダクタイル鋳鉄管で構成された管路を模擬するため、真値がＭ＝１３．３０７０ｋｇ、Ｋ＝２．８９９４×１０９Ｎ／ｍ、Ｃ＝１０００Ｎｓ／ｍ、ΔＫ＝５．７９８９×１０７Ｎ／ｍ、ΔＣ＝６６６．６６６７Ｎｓ／ｍとした。本条件よりサンプリング周波数５０ｋＨｚで２０ｍｓのインパルス応答関数を生成し、さらに平均０、分散１の正規性白色雑音をインパルス応答関数に加算し、試験用のテストデータとした。初期値は各真値に０．９５を乗じた値を用い、ステップ調整用パラメータは０．０１を用いた。
図８は、上記の条件によるシステム同定結果を示す図である。同図において横軸は周波数を表し、縦軸はアクセレランス（Accelerance）を表す。同図には、同定したパラメータを用いて計算した自己周波数応答関数（Identified）と、真値の自己周波数応答関数（Experiment）をそれぞれ示している。同図によれば、真値と、同定結果の両者の良い一致が確認できる。なお、推定値は、Ｍ＝１３．３０７３ｋｇ、Ｋ＝２．８９８０×１０９Ｎ／ｍ、Ｃ＝９９９．９３１２Ｎｓ／ｍ、ΔＫ＝６．０６５２×１０７Ｎ／ｍ、ΔＣ＝６６６．６７３０Ｎｓ／ｍであった。以上により、本実施形態による同定アルゴリズムの動作が確認できる。
なお、第２の実施形態のように管路のシステム同定を行う場合、パラメータωｄ１、ωｄ２、α、β、γの初期値は以下のように算出される。これらの初期値を決める主たるパラメータである質量Ｍ及びばね定数Ｋは、以下の式（６）を用いて算出される。
ここで、Ｒは半径、Ａは断面積である。管長をＬ、管厚をｈとした場合、Ａ＝ｈＬの関係がある。半径Ｒ、管厚ｈは管路管理台帳データから読み出した口径、管厚から求められる。管長Ｌは、管路管理台帳データから読み出してもよく、ユーザが入力してもよい。Ｅは管の弾性係数、Ｉは断面二次モーメントであり、Ｉ＝Ｌ・ｈ３／１２である。ρは、管の密度である。弾性係数Ｅ及び管の密度ρは、管路管理台帳データから読み出した材種に応じた値としてもよく、予め決められた値としてもよい。ΔＫはばね定数Ｋの１／１００、ΔＣは減衰係数Ｃの１／３程度が望ましい。初期値の算出に用いる減衰係数Ｃ、減衰係数の変化ΔＣは、管路管理台帳データから読み出した口径、管厚、材種のうち一以上に応じた値としてもよく、予め決められた値としてもよい。初期値設定部３０７は、これらの値を用いて式（３）により各パラメータωｄ１、ωｄ２、α、β、γの初期値を算出する。
以上が初期値設定部の具体的な表式である。
以上が初期値設定部の具体的な表式である。
［実験］
供用後の水道管を試験用管路に設置し、通水環境下でシステム同定実験を実施した。供試管路は口径１００ｍｍ、肉厚１０ｍｍの普通鋳鉄管で供用済みの管を用いた。管路の上部には消火栓を設置し、カプラ上に加速度センサを設置した。
供用後の水道管を試験用管路に設置し、通水環境下でシステム同定実験を実施した。供試管路は口径１００ｍｍ、肉厚１０ｍｍの普通鋳鉄管で供用済みの管を用いた。管路の上部には消火栓を設置し、カプラ上に加速度センサを設置した。
図９は、システム同定実験の結果を示す図である。同図において横軸は周波数を表し、縦軸はアクセレランス（Accelerance）を表す。また図中の丸のプロットは自己周波数応答関数の実験値（Experiment）を表し、実線は同定結果（Identified）を表している。同図に示すように、実験値と推定値のピーク位置やスペクトル形状について、良好な一致が確認された。なお、得られた推定値はＭ＝１４．８４４６ｋｇ、Ｋ＝３．４３４６×１０９Ｎ／ｍ、Ｃ＝９０３．５４９１Ｎｓ／ｍ、ΔＫ＝１．２８７７×１０８Ｎ／ｍ、ΔＣ＝９０３．５４９１Ｎｓ／ｍであった。
関連技術と比較した優位性を示すため、特許文献１で用いられているＡＲモデルを用いた同定を行った。そして、ＡＲモデルを用いた同定方法と本実施形態による同定方法とを比較した。
その結果、ＡＲ法による同定結果は、実験値との大きな乖離が確認され、同定困難であることがわかった。これは、インパルス応答関数がうなり波形をともなっており、時間領域同定法で用いる多項式モデルでは、本波形の特徴を記述するのに限界が生じていることを示している。以上の検討により、本実施形態が優れた効果を奏することがわかった。
図１０は、本発明の実施形態によるシステム同定装置の最小構成を示すブロック図である。同図に示す最小構成のシステム同定装置１ａは、上述した第１の実施形態の分析部１０５を少なくとも備えればよい。分析部１０５は、分析対象を励振させた位置において計測された入力信号及び出力信号に基づいて自己周波数応答関数を算出する。そして、分析部１０５は、算出された自己周波数応答関数から得られるインパルス応答関数と、分析対象がモデル化された仮想二自由度モデルのインパルス応答関数とを用いて、分析対象のシステム同定を行う。
本実施形態によれば、仮想二自由度モデルの自己周波数応答関数に対応するインパルス応答関数を用いて、近接固有値を有する系のシステム同定を行うことができる。
上述した実施形態におけるシステム同定装置１、１ａ、３、５は、バスで接続されたＣＰＵ（Central Processing Unit）やメモリや補助記憶装置などを備え、システム同定プログラムを実行することによって上述した実施形態におけるシステム同定装置１、１ａ、３、５の一部の機能を実現する。なお、システム同定装置１、１ａ、３、５の一部の機能は、ＡＳＩＣ（Application Specific Integrated Circuit）やＰＬＤ（Programmable Logic Device）やＦＰＧＡ（Field Programmable Gate Array）等のハードウェアを用いて実現されても良い。システム同定プログラムは、コンピュータ読み取り可能な記録媒体に記録されても良い。コンピュータ読み取り可能な記録媒体とは、例えばフレキシブルディスク、光磁気ディスク、ＲＯＭ（Read Only Memory）、ＣＤ−ＲＯＭ（Compact Disc Read only memory）等の可搬媒体、コンピュータシステムに内蔵されるハードディスク等の記憶装置である。システム同定プログラムは、電気通信回線を介して送信されても良い。
上記の実施形態の一部又は全部は、以下の付記のようにも記載されうるが、以下には限られない。
（付記１）分析対象を励振させた位置において計測された入力信号及び出力信号に基づいて自己周波数応答関数を算出し、算出された前記自己周波数応答関数から得られるインパルス応答関数と、前記分析対象がモデル化された仮想二自由度モデルのインパルス応答関数とを用いて前記分析対象のシステム同定を行う分析部、を備えるシステム同定装置。
（付記２）前記分析部は、多変数ニュートン法により仮想二自由度モデルの前記インパルス応答関数を推定し、推定により得られた前記インパルス応答関数に基づいて前記分析対象のシステム同定を行う、付記１に記載のシステム同定装置。
（付記３）前記分析対象の物理データに基づいて多変数ニュートン法に用いられる初期値を算出する初期値設定部をさらに備える、付記２に記載のシステム同定装置。
（付記４）前記分析対象は、管路である、付記１に記載のシステム同定装置。
（付記５）前記分析対象を励振させる加振部と、前記加振部により前記分析対象を励振させた位置における入力信号及び出力信号を計測する計測部と、前記加振部が前記分析対象を励振させる位置と前記計測部が前記分析対象を計測する位置とを揃えるように前記加振部及び前記計測部を設置する設置位置決め部とをさらに備える、付記１に記載のシステム同定装置。
（付記６）前記加振部は、力センサを内蔵したインパルスハンマー又は電磁式加振器である、付記５に記載のシステム同定装置。
（付記７）前記計測部は、加速度ピックアップ、レーザー変位計、レーザードップラー型速度計又は接触式変位計である、付記５に記載のシステム同定装置。
（付記８）前記設置位置決め部は、消火栓カプラである、付記５に記載のシステム同定装置。
（付記９）前記分析部は、前記システム同定により前記分析対象の系の質量、剛性定数及び減衰係数を得る、付記１に記載のシステム同定装置。
（付記１０）分析対象を励振させる加振ステップと、前記加振ステップにおいて前記分析対象を励振させた位置における入力信号及び出力信号を計測する計測ステップと、前記計測ステップにおいて計測された前記入力信号及び前記出力信号に基づいて自己周波数応答関数を算出し、算出された前記自己周波数応答関数から得られるインパルス応答関数と、前記分析対象がモデル化された仮想二自由度モデルのインパルス応答関数とを用いて前記分析対象のシステム同定を行う分析ステップと、を有するシステム同定方法。
（付記１１）コンピュータに、分析対象を励振させた位置において計測された入力信号及び出力信号に基づいて自己周波数応答関数を算出し、算出された前記自己周波数応答関数から得られるインパルス応答関数と、前記分析対象がモデル化された仮想二自由度モデルのインパルス応答関数とを用いて前記分析対象のシステム同定を行う分析ステップ、を実行させるプログラム。
以上、実施形態（及び実施例）を参照して本願発明を説明したが、本願発明は上記実施形態（及び実施例）に限定されるものではない。本願発明の構成や詳細には、本願発明のスコープ内で当業者が理解し得る様々な変更をすることができる。
この出願は、２０１８年２月２１日に出願された日本出願特願２０１８−０２９２１８を基礎とする優先権を主張し、その開示の全てをここに取り込む。
以上、実施形態（及び実施例）を参照して本願発明を説明したが、本願発明は上記実施形態（及び実施例）に限定されるものではない。本願発明の構成や詳細には、本願発明のスコープ内で当業者が理解し得る様々な変更をすることができる。
この出願は、２０１８年２月２１日に出願された日本出願特願２０１８−０２９２１８を基礎とする優先権を主張し、その開示の全てをここに取り込む。
近接固有値を有する系のシステム同定に適用可能である。よってその工業的価値は高い。
１、１ａ、３、５…システム同定装置
１０１、３０１…設置位置決め部
１０２、３０２…加振部
１０３、３０３…計測部
１０４、３０４…信号収集部
１０５、３０５…分析部
１０６、３０８…対象物理システム
３０６…記憶部
３０７…初期値設定部
５０１…消火栓カプラ
５０２…ハンマー
５０３…センサ
５０４…データロガー
５０５…同定処理部
５０６…管路
１０１、３０１…設置位置決め部
１０２、３０２…加振部
１０３、３０３…計測部
１０４、３０４…信号収集部
１０５、３０５…分析部
１０６、３０８…対象物理システム
３０６…記憶部
３０７…初期値設定部
５０１…消火栓カプラ
５０２…ハンマー
５０３…センサ
５０４…データロガー
５０５…同定処理部
５０６…管路
Claims (10)
- 分析対象を励振させた位置において計測された入力信号及び出力信号に基づいて自己周波数応答関数を算出し、算出された前記自己周波数応答関数から得られるインパルス応答関数と、前記分析対象がモデル化された仮想二自由度モデルのインパルス応答関数とを用いて前記分析対象のシステム同定を行う分析部、
を備えるシステム同定装置。 - 前記分析部は、多変数ニュートン法により仮想二自由度モデルの前記インパルス応答関数を推定し、推定により得られた前記インパルス応答関数に基づいて前記分析対象のシステム同定を行う、
請求項１に記載のシステム同定装置。 - 前記分析対象の物理データに基づいて多変数ニュートン法に用いられる初期値を算出する初期値設定部をさらに備える、
請求項２に記載のシステム同定装置。 - 前記分析対象は、管路である、
請求項１に記載のシステム同定装置。 - 前記分析対象を励振させる加振部と、
前記加振部により前記分析対象を励振させた位置における入力信号及び出力信号を計測する計測部と、
前記加振部が前記分析対象を励振させる位置と前記計測部が前記分析対象を計測する位置とを揃えるように前記加振部及び前記計測部を設置する設置位置決め部とをさらに備える、
請求項１に記載のシステム同定装置。 - 前記加振部は、力センサを内蔵したインパルスハンマー又は電磁式加振器である、
請求項５に記載のシステム同定装置。 - 前記計測部は、加速度ピックアップ、レーザー変位計、レーザードップラー型速度計又は接触式変位計である、
請求項５に記載のシステム同定装置。 - 前記分析部は、前記システム同定により前記分析対象の系の質量、剛性定数及び減衰係数を得る、請求項１に記載のシステム同定装置。
- 分析対象を励振させる加振ステップと、
前記加振ステップにおいて前記分析対象を励振させた位置における入力信号及び出力信号を計測する計測ステップと、
前記計測ステップにおいて計測された前記入力信号及び前記出力信号に基づいて自己周波数応答関数を算出し、算出された前記自己周波数応答関数から得られるインパルス応答関数と、前記分析対象がモデル化された仮想二自由度モデルのインパルス応答関数とを用いて前記分析対象のシステム同定を行う分析ステップと、
を有するシステム同定方法。 - コンピュータに、
分析対象を励振させた位置において計測された入力信号及び出力信号に基づいて自己周波数応答関数を算出し、算出された前記自己周波数応答関数から得られるインパルス応答関数と、前記分析対象がモデル化された仮想二自由度モデルのインパルス応答関数とを用いて前記分析対象のシステム同定を行う分析ステップ、
を実行させるコンピュータプログラム。
Applications Claiming Priority (3)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP2018029218 | 2018-02-21 | ||
| JP2018029218 | 2018-02-21 | ||
| PCT/JP2019/005805 WO2019163701A1 (ja) | 2018-02-21 | 2019-02-18 | システム同定装置、システム同定方法及び記録媒体 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| JPWO2019163701A1 JPWO2019163701A1 (ja) | 2021-02-04 |
| JP6981526B2 true JP6981526B2 (ja) | 2021-12-15 |
Family
ID=67688228
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| JP2020501746A Active JP6981526B2 (ja) | 2018-02-21 | 2019-02-18 | システム同定装置、システム同定方法及びコンピュータプログラム |
Country Status (4)
| Country | Link |
|---|---|
| US (1) | US20210010980A1 (ja) |
| EP (1) | EP3757704A4 (ja) |
| JP (1) | JP6981526B2 (ja) |
| WO (1) | WO2019163701A1 (ja) |
Families Citing this family (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111830062A (zh) * | 2020-07-22 | 2020-10-27 | 深圳市赛龙自动化科技有限公司 | 一种陶瓷瓶瑕疵检测设备及其检测方法 |
| CN112765863B (zh) * | 2021-02-04 | 2022-06-10 | 上海交通大学 | 基于位姿依赖特性和交叉耦合项的机器人刀尖频响预测方法和系统 |
| CN116878801B (zh) * | 2023-06-01 | 2024-08-30 | 哈尔滨工业大学(深圳)(哈尔滨工业大学深圳科技创新研究院) | 随机激励作用下的结构状态识别方法、装置、终端 |
Family Cites Families (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JPH03217901A (ja) * | 1990-01-24 | 1991-09-25 | Toshiba Corp | システムの同定装置 |
| JPH0477798A (ja) * | 1990-07-19 | 1992-03-11 | Oki Electric Ind Co Ltd | 周波数包絡線成分の特徴量抽出方法 |
| US20050072234A1 (en) * | 2003-05-20 | 2005-04-07 | Weidong Zhu | System and method for detecting structural damage |
| EP2682729A1 (en) * | 2012-07-05 | 2014-01-08 | Vrije Universiteit Brussel | Method for determining modal parameters |
| WO2015059956A1 (ja) * | 2013-10-23 | 2015-04-30 | 日本電気株式会社 | 構造物診断装置、構造物診断方法、及びプログラム |
| US10387116B2 (en) * | 2014-02-07 | 2019-08-20 | Mitsubishi Electric Corporation | System identification device |
| JPWO2016114136A1 (ja) * | 2015-01-14 | 2017-11-02 | 日本電気株式会社 | 配管検査システム、配管検査装置、配管検査方法、及び、プログラム |
| JP6725782B2 (ja) | 2016-08-15 | 2020-07-22 | エイチ・シー・ネットワークス株式会社 | アンテナ装置 |
-
2019
- 2019-02-18 JP JP2020501746A patent/JP6981526B2/ja active Active
- 2019-02-18 US US16/970,754 patent/US20210010980A1/en not_active Abandoned
- 2019-02-18 WO PCT/JP2019/005805 patent/WO2019163701A1/ja unknown
- 2019-02-18 EP EP19756575.7A patent/EP3757704A4/en active Pending
Also Published As
| Publication number | Publication date |
|---|---|
| US20210010980A1 (en) | 2021-01-14 |
| EP3757704A1 (en) | 2020-12-30 |
| JPWO2019163701A1 (ja) | 2021-02-04 |
| WO2019163701A1 (ja) | 2019-08-29 |
| EP3757704A4 (en) | 2021-04-07 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| JP6981526B2 (ja) | システム同定装置、システム同定方法及びコンピュータプログラム | |
| EP2904368B1 (en) | Turbine blade fatigue life analysis using non-contact measurement and dynamical response reconstruction techniques | |
| US7779690B2 (en) | Vibrating wire sensor using spectral analysis | |
| JP4994448B2 (ja) | モーダルパラメータの推定方法および装置 | |
| Larsen et al. | Modal analysis of wind turbine blades | |
| Law et al. | Crack identification in beam from dynamic responses | |
| KR101747116B1 (ko) | 주파수 응답 기반 교량 내하력 평가 방법 | |
| Zhou et al. | Efficient hysteresis loop analysis-based damage identification of a reinforced concrete frame structure over multiple events | |
| JP3837099B2 (ja) | 構造物の損傷推定システムおよびプログラム | |
| Catbas et al. | Modal analysis of multi-reference impact test data for steel stringer bridges | |
| Debut et al. | Identification of the nonlinear excitation force acting on a bowed string using the dynamical responses at remote locations | |
| Smail et al. | ARMA models for modal analysis: effect of model orders and sampling frequency | |
| JP6364742B2 (ja) | 構造物診断装置、構造物診断方法、及びプログラム | |
| JP2005083975A (ja) | 構造性能指標推定装置、及び構造物の構造性能リアルタイムモニタリング方法 | |
| Ganguly et al. | Performance assessment of time-domain damage indicators based on output-only measurement and Poincaré map: A comparative review on nonlinear structures | |
| Kasinos et al. | Using the vibration envelope as a damage-sensitive feature in composite beam structures | |
| Smail et al. | Assessment of optimal ARMA model orders for modal analysis | |
| Rossing | Modal analysis | |
| JPWO2019093294A1 (ja) | 推定装置、推定方法及びコンピュータプログラム | |
| Wu et al. | Input force identification: application to soil–pile interaction | |
| US20220397500A1 (en) | Apparatus and method for detecting microcrack using orthogonality analysis of mode shape vector and principal plane in resonance point | |
| WO2020179241A1 (ja) | 構造物診断装置、構造物診断方法、及びコンピュータ読み取り可能な記録媒体 | |
| KR101557270B1 (ko) | 스마트 구조물의 유지 관리를 위한 est 기반 단일 계측 시스템 | |
| JP2018200217A (ja) | 打音検査装置及び打音検査方法 | |
| Shariati et al. | Oversampling in virtual visual sensors as a means to recover higher modes of vibration |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| A521 | Request for written amendment filed |
Free format text: JAPANESE INTERMEDIATE CODE: A523 Effective date: 20200813 |
|
| A621 | Written request for application examination |
Free format text: JAPANESE INTERMEDIATE CODE: A621 Effective date: 20200813 |
|
| A521 | Request for written amendment filed |
Free format text: JAPANESE INTERMEDIATE CODE: A523 Effective date: 20200903 |
|
| TRDD | Decision of grant or rejection written | ||
| A01 | Written decision to grant a patent or to grant a registration (utility model) |
Free format text: JAPANESE INTERMEDIATE CODE: A01 Effective date: 20211019 |
|
| A61 | First payment of annual fees (during grant procedure) |
Free format text: JAPANESE INTERMEDIATE CODE: A61 Effective date: 20211101 |
|
| R150 | Certificate of patent or registration of utility model |
Ref document number: 6981526 Country of ref document: JP Free format text: JAPANESE INTERMEDIATE CODE: R150 |