特許 6981526 システム同定装置、システム同定方法及びコンピュータプログラム

(19)【発行国】日本国特許庁(JP)

(12)【公報種別】特許公報(B2)

(11)【特許番号】6981526

(24)【登録日】2021年11月22日

(45)【発行日】2021年12月15日

(54)【発明の名称】システム同定装置、システム同定方法及びコンピュータプログラム

(51)【国際特許分類】

   G05B 23/02 20060101AFI20211202BHJP

   G01M 7/02 20060101ALI20211202BHJP

【ＦＩ】

   G05B23/02 P

   G01M7/02 Z

【請求項の数】10

【全頁数】15

(21)【出願番号】特願2020-501746(P2020-501746)

(86)(22)【出願日】2019年2月18日

(86)【国際出願番号】JP2019005805

(87)【国際公開番号】WO2019163701

(87)【国際公開日】20190829

【審査請求日】2020年8月13日

(31)【優先権主張番号】特願2018-29218(P2018-29218)

(32)【優先日】2018年2月21日

(33)【優先権主張国】JP

(73)【特許権者】

【識別番号】000004237

【氏名又は名称】日本電気株式会社

(74)【代理人】

【識別番号】100109313

【弁理士】

【氏名又は名称】机 昌彦

(74)【代理人】

【識別番号】100124154

【弁理士】

【氏名又は名称】下坂 直樹

(72)【発明者】

【氏名】高田 宗一朗

(72)【発明者】

【氏名】井上 裕文

(72)【発明者】

【氏名】篠田 茂樹

(72)【発明者】

【氏名】菊池 克

【審査官】 黒田 暁子

(56)【参考文献】

【文献】 特開平４−７７７９８（ＪＰ，Ａ）

【文献】 特開平３−２１７９０１（ＪＰ，Ａ）

【文献】 国際公開第２０１５／０５９９５６（ＷＯ，Ａ１）

【文献】 国際公開第２０１６／１１４１３６（ＷＯ，Ａ１）

【文献】 国際公開第２０１５／１１８７３７（ＷＯ，Ａ１）

【文献】 小松 正貴 ほか，実現理論による近接固有値を有する構造物の振動特性推定，構造工学論文集，土木学会，2013年，Vol. 59A，p. 340-352，ISSN 1881-820X

【文献】 中溝 高好，信号解析とシステム同定，東京: コロナ社，1988年，p. 22-24, 49-53, 121-127，ISBN 4-339-03081-3

(58)【調査した分野】（Int.Cl.，ＤＢ名）

Ｇ０５Ｂ ２３／０２

Ｇ０１Ｍ ７／０２

(57)【特許請求の範囲】

【請求項1】

分析対象を励振させた位置において計測された入力信号及び出力信号に基づいて自己周波数応答関数を算出し、算出された前記自己周波数応答関数から得られるインパルス応答関数と、前記分析対象がモデル化された仮想二自由度モデルのインパルス応答関数とを用いて前記分析対象のシステム同定を行う分析部、
を備えるシステム同定装置。

【請求項2】

前記分析部は、多変数ニュートン法により仮想二自由度モデルの前記インパルス応答関数を推定し、推定により得られた前記インパルス応答関数に基づいて前記分析対象のシステム同定を行う、
請求項１に記載のシステム同定装置。

【請求項3】

前記分析対象の物理データに基づいて多変数ニュートン法に用いられる初期値を算出する初期値設定部をさらに備える、
請求項２に記載のシステム同定装置。

【請求項4】

前記分析対象は、管路である、
請求項１に記載のシステム同定装置。

【請求項5】

前記分析対象を励振させる加振部と、
前記加振部により前記分析対象を励振させた位置における入力信号及び出力信号を計測する計測部と、
前記加振部が前記分析対象を励振させる位置と前記計測部が前記分析対象を計測する位置とを揃えるように前記加振部及び前記計測部を設置する設置位置決め部とをさらに備える、
請求項１に記載のシステム同定装置。

【請求項6】

前記加振部は、力センサを内蔵したインパルスハンマー又は電磁式加振器である、
請求項５に記載のシステム同定装置。

【請求項7】

前記計測部は、加速度ピックアップ、レーザー変位計、レーザードップラー型速度計又は接触式変位計である、
請求項５に記載のシステム同定装置。

【請求項8】

前記分析部は、前記システム同定により前記分析対象の系の質量、剛性定数及び減衰係数を得る、請求項１に記載のシステム同定装置。

【請求項9】

分析対象を励振させる加振ステップと、
前記加振ステップにおいて前記分析対象を励振させた位置における入力信号及び出力信号を計測する計測ステップと、
前記計測ステップにおいて計測された前記入力信号及び前記出力信号に基づいて自己周波数応答関数を算出し、算出された前記自己周波数応答関数から得られるインパルス応答関数と、前記分析対象がモデル化された仮想二自由度モデルのインパルス応答関数とを用いて前記分析対象のシステム同定を行う分析ステップと、
を有するシステム同定方法。

【請求項10】

コンピュータに、
分析対象を励振させた位置において計測された入力信号及び出力信号に基づいて自己周波数応答関数を算出し、算出された前記自己周波数応答関数から得られるインパルス応答関数と、前記分析対象がモデル化された仮想二自由度モデルのインパルス応答関数とを用いて前記分析対象のシステム同定を行う分析ステップ、
を実行させるコンピュータプログラム。

【発明の詳細な説明】

【技術分野】

【0001】

本発明は、システム同定装置、システム同定方法及び記録媒体に関する。

【背景技術】

【0002】

石油、ガス、水道などのプラントシステムや、産業用ロボット等の物理システムを対象として、ＩｏＴ（Internet of Things）による監視や制御を行う際には、対象とするシステムのモデリング技術が重要である。このモデリング技術として、観測データに基づいて物理システムの数理的モデリングを行うものがある（例えば、特許文献１〜４、非特許文献１参照）。

【先行技術文献】

【特許文献】

【0003】

【特許文献1】国際公開２０１５／１１８７３７号

【特許文献2】国際公開２０１５／０５９９５６号

【特許文献3】特開平４−７７７９８号公報

【特許文献4】特開平３−２１７９０１号公報

【非特許文献】

【0004】

【非特許文献1】中溝 高好、「信号解析とシステム同定」、コロナ社、ｐ．２２−２４、４９−５３，１２１−１２７、１９８８年

【発明の概要】

【発明が解決しようとする課題】

【0005】

しかし、モデリングの対象システムが、うなり現象や共鳴が発生するような、近接固有値を有する系である場合、特許文献１〜４および非特許文献１に記載の技術では、システム同定が困難な場合があった。なお、近接固有値を有する系は、固有値が縮退した重根を有する系も含むものとする。
本発明は、近接固有値を有する系のシステム同定を行うことができるシステム同定装置等を提供することを目的としている。

【課題を解決するための手段】

【0006】

本発明の第１の態様によれば、システム同定装置は、分析対象を励振させた位置において計測された入力信号及び出力信号に基づいて自己周波数応答関数を算出し、算出された前記自己周波数応答関数から得られるインパルス応答関数と、前記分析対象がモデル化された仮想二自由度モデルのインパルス応答関数とを用いて前記分析対象のシステム同定を行う分析部と、を備える。

【0007】

本発明の第２の態様によれば、システム同定方法は、分析対象を励振させる加振ステップと、前記加振ステップにおいて前記分析対象を励振させた位置における入力信号及び出力信号を計測する計測ステップと、前記計測ステップにおいて計測された前記入力信号及び前記出力信号に基づいて自己周波数応答関数を算出し、算出された前記自己周波数応答関数から得られるインパルス応答関数と、前記分析対象がモデル化された仮想二自由度モデルのインパルス応答関数とを用いて前記分析対象のシステム同定を行う分析ステップと、を有する。

【0008】

本発明の第３の態様によれば、記録媒体は、コンピュータに、分析対象を励振させた位置において計測された入力信号及び出力信号に基づいて自己周波数応答関数を算出し、算出された前記自己周波数応答関数から得られるインパルス応答関数と、前記分析対象がモデル化された仮想二自由度モデルのインパルス応答関数とを用いて前記分析対象のシステム同定を行う分析ステップ、を実行させるプログラムを記録する。

【発明の効果】

【0009】

本発明によれば、近接固有値を有する系のシステム同定を行うことができる。

【図面の簡単な説明】

【0010】

【図1】本発明の第１の実施形態によるシステム同定装置の構成を示すブロック図である。

【図2】同実施形態によるシステム同定装置の処理を示すフローチャートである。

【図3】第２の実施形態によるシステム同定装置の構成を示すブロック図である。

【図4】同実施形態によるシステム同定装置の処理を示すフローチャートである。

【図5】第３の実施形態によるシステム同定装置の構成を示すブロック図である。

【図6】同実施形態による消火栓カプラの状況を示す図である。

【図7】同実施形態による仮想の仮想二自由度モデルを示す図である。

【図8】同実施形態によるシステム同定結果を示す図である。

【図9】システム同定実験の結果を示す図である。

【図10】本発明の実施形態によるシステム同定装置の最小構成を示すブロック図である。

【発明を実施するための形態】

【0011】

以下、図面を参照して本発明の一実施の形態を説明する。
観測データに基づく物理システムの数理的なモデリング手法は、「システム同定問題」と呼ばれている。本問題では大別して、（１）システムの入力および出力の信号が既知の場合と、（２）入力が未知の場合とに分類される。さらには、時間領域信号又は周波数領域信号を用いた技術も知られている。

【0012】

時間領域情報を用いた技術では、ＡＲモデル(Autoregressive model：自己回帰モデル)やＭＡモデル(Moving Average model：移動平均モデル)、ＡＲＭＡモデル(Autoregressive moving average model：自己回帰移動平均モデル)、ＡＲＸモデル（Auto-Regressive eXogeneous model）などの多項式モデルが用いられる。多項式モデルでは、周波数領域が平坦化され、異なる固有値が近接している近接固有値を有する系への適用が困難であった。一方で、周波数領域情報を用いた場合には、固有値双方のピーク位置が不明瞭となる場合が多く、曲線適合法を適用することができない。特に、フーリエ変換による周波数領域のサンプル数によっては、ピーク位置を視認することすらできない。

【0013】

上記のことから、近接固有値を有する系では固有値が近接するため、周波数領域同定法および時間領域同定法では、システム同定が困難となる。このように、近接固有値を有する系のシステム同定問題は未だ十分に解決されていない問題である。本実施形態のシステム同定装置等は、このような問題を解決し、近接固有値を有する系（固有値が重なる重根を有する系も含む）のシステム同定を行う。

【0014】

［第１の実施形態］
図１は、第１の実施形態によるシステム同定装置１の構成を示すブロック図である。システム同定装置１は、設置位置決め部１０１、加振部１０２、計測部１０３、信号収集部１０４及び分析部１０５を有する。対象物理システム１０６は、システム同定装置１による同定対象である。設置位置決め部１０１は、対象物理システム１０６に、加振部１０２及び計測部１０３を設置する。加振部１０２は、設置位置決め部１０１を介して対象物理システム１０６を励振させる。計測部１０３は、加振部１０２が設置位置決め部１０１を介して対象物理システム１０６を励振させたときの、対象物理システム１０６への入力信号と出力信号を検知する。信号収集部１０４は、計測部１０３が検知した入力信号及び出力信号をデータ化する。分析部１０５は、信号収集部１０４により得られたデータを分析し、対象物理システム１０６のシステム同定を行う。

【0015】

図２は、システム同定装置１の処理を示すフローチャートである。測定者は、対象物理システム１０６に対し、設置位置決め部１０１を介して加振部１０２及び計測部１０３を設置する（ステップＳ１１０）。このとき、入力位置と出力位置とを一致させるようにする。入力位置とは、加振部１０２が設置される、すなわち、振動の原因が入力される位置である。出力位置とは、計測部１０３が設定される、すなわち、対象物理システム１０６の振動を計測する位置である。

【0016】

自己周波数応答関数を計測するため、加振部１０２により、設置位置決め部１０１を介して対象物理システム１０６を励振させる。計測部１０３は、対象物理システム１０６への振動の入力信号と対象物理システム１０６からの振動の出力信号を検知する。信号収集部１０４は、計測部１０３が検知した入力信号及び出力信号をデータ化し、分析部１０５に出力する。分析部１０５は、得られたデータを分析する。具体的には、分析部１０５は、最初に、入力信号と出力信号をそれぞれ高速フーリエ変換（ＦＦＴ：fast fourier transform）する。分析部１０５は、周波数領域で出力信号を入力信号により除算することにより、自己周波数応答関数を求める（ステップＳ１２０）。

【0017】

次いで、分析部１０５は、自己周波数応答関数において、対象とする近接固有値が存在する周波数帯域のみズーミングを行う（ステップＳ１３０）。分析部１０５は、ズーミングした自己周波数応答関数に逆フーリエ変換（Inverse fourier transform）を施すことにより、自己周波数応答関数のインパルス応答関数を求める（ステップＳ１４０）。

【0018】

分析部１０５は、次のステップＳ１６０において用いられる初期値及びステップ幅の入力を受ける（ステップＳ１５０）。分析部１０５は、ステップＳ１４０で得られたインパルス応答関数に対し、後述する仮想二自由度系のインパルス応答関数を用いた多変数ニュートン法を適用する（ステップＳ１６０）。分析部１０５は、多変数ニュートン法の実行時に、ステップＳ１５０において入力された初期値及びステップ幅を用いる。

【0019】

分析部１０５は、解が収束しないと判定した場合（ステップＳ１７０：ＮＯ）、ステップＳ１５０に戻って、新たに使用する初期値及びステップ幅の入力を受け、ステップＳ１６０の処理を行う。分析部１０５は、解が収束したと判定した場合（ステップＳ１７０：ＹＥＳ）、その収束したときの仮想二自由度系のインパルス応答関数に基づいて対象物理システム１０６の系の質量、剛性定数及び減衰係数を求め、システム同定を行う（ステップＳ１８０）。
本実施形態によれば、近接固有値を有する系のシステム同定が可能である。

【0020】

［第２の実施形態］
本実施形態は、第１の実施形態を管路のシステム同定に適用する。
図３は、第２の実施形態によるシステム同定装置３の構成を示すブロック図である。システム同定装置３は、設置位置決め部３０１、加振部３０２、計測部３０３、信号収集部３０４、分析部３０５、記憶部３０６及び初期値設定部３０７を有する。対象物理システム３０８は、管路であり、システム同定装置３による同定対象である。

【0021】

設置位置決め部３０１、加振部３０２、計測部３０３及び信号収集部３０４はそれぞれ、第１の実施形態の設置位置決め部１０１、加振部１０２、計測部１０３及び信号収集部１０４と同様の機能を有する。記憶部３０６は、管路管理台帳データを記憶する。管路管理台帳データは、対象物理システム３０８の物理データである口径、材種、管厚の情報を含む。初期値設定部３０７は、記憶部３０６から読み出した口径、材種、管厚のデータに基づいて、分析部３０５において用いられる多変数ニュートン法のパラメータの初期値を算出する。分析部３０５は、初期値設定部３０７が算出した初期値を多変数ニュートン法において用いる点以外は、第１の実施形態の分析部１０５と同様の機能を有する。

【0022】

図４は、システム同定装置３の処理を示すフローチャートである。測定者は、対象物理システム３０８に対し、設置位置決め部３０１を介して加振部３０２及び計測部３０３を設置する（ステップＳ３１０）。加振部３０２により、設置位置決め部３０１を介して対象物理システム３０８を励振させる。計測部３０３は、加振位置における対象物理システム３０８への入力信号と対象物理システム３０８からの出力信号を検知する。信号収集部３０４は、計測部３０３が検知した入力信号及び出力信号をデータ化する。分析部３０５は、入力信号と出力信号を用いて自己周波数応答関数を求める（ステップＳ３２０）。分析部３０５は、自己周波数応答関数において、対象とする近接固有値が存在する周波数帯域のみズーミングを行う（ステップＳ３３０）。分析部３０５は、ズーミング部分を逆フーリエ変換し、自己周波数応答関数のインパルス応答関数を求める（ステップＳ３４０）。

【0023】

初期値設定部３０７は、記憶部３０６から対象物理システム３０８の口径、材種、管厚のデータを読み出す（ステップＳ３５０）。初期値設定部３０７は、読み出したデータに基づいて多変数ニュートン法のパラメータの初期値を算出する（ステップＳ３６０）。分析部３０５は、ステップ幅の入力を受ける（ステップＳ３７０）。分析部３０５は、ステップＳ３４０において得られたインパルス応答関数に対し、後述する仮想二自由度系のインパルス応答関数を用いた多変数ニュートン法を適用する（ステップＳ３８０）。多変数ニュートン法においては、ステップＳ３６０において算出された初期値と、ステップＳ３７０において入力されたステップ幅とを用いる。

【0024】

分析部３０５は、解が収束しないと判定した場合（ステップＳ３９０：ＮＯ）、ステップＳ３７０に戻って、新たに使用するステップ幅の入力を受け、ステップＳ３８０の処理を行う。分析部３０５は、解が収束したと判定した場合（ステップＳ３９０：ＹＥＳ）、その収束したときの仮想二自由度系のインパルス応答関数に基づいて対象物理システム３０８の系の質量、剛性定数及び減衰係数を求め、システム同定を行う（ステップＳ４００）。
本実施形態によれば、管路のシステム同定が可能である。

【0025】

［第３の実施形態］
ここでは、第１の実施形態を水道管路のシステム同定に適用する。
図５は、本実施形態のシステム同定装置５の構成を示すブロック図である。システム同定装置５は、消火栓カプラ５０１、ハンマー５０２、センサ５０３、データロガー５０４及び同定処理部５０５を有する。管路５０６は、システム同定装置５によるシステム同定対象の水道管路である。ハンマー５０２として、力センサを内蔵したインパルスハンマー、市販ハンマーに加速度ピックアップを設置したもの、電磁式加振器などが例示できる。センサ５０３として、加速度ピックアップ、レーザードップラー型速度計、レーザー変位計、接触式変位計などが例示できる。同定処理部５０５は、例えば、プロセッサ、メモリ及びＨＤＤ（Hard disk drive）により実現される。プロセッサは、処理をコンピュータに実行させるための同定処理プログラムをＨＤＤから読み出して実行することにより同定処理部５０５として動作する。

【0026】

図６は、管路５０６に設置された消火栓カプラ５０１の状況を示す図である。計測者は、管路５０６に消火栓カプラ５０１及びセンサ５０３を設置する（図２のステップＳ１１０）。計測者は、ハンマー５０２により図６に示す消火栓カプラ５０１を打撃し、管路５０６を励振する。センサ５０３は、打撃位置（Tapping Point）と同じ測定位置（Measurements point for vibration response）における励振後の出力信号を検出する。データロガー５０４は、ハンマー５０２の入力信号及びセンサ５０３の出力信号を収集する。データロガー５０４は、収集した入力信号及び出力信号をデータ化して同定処理部５０５に出力する。

【0027】

同定処理部５０５は、入力信号と出力信号それぞれに高速フーリエ変換（ＦＦＴ）処理を行う。ＦＦＴにより得られた入力信号のスペクトルをＸ（ω）、出力信号のスペクトルをＹ（ω）とそれぞれ表す。ωは周波数である。同定処理部５０５は、各周波数領域のスペクトルを除し、自己周波数応答関数Ｌ（ω）＝Ｙ（ω）／Ｘ（ω）を求める（図２のステップＳ１２０）。ここで、自己周波数応答関数の算出にあたり、Ｌ（ω）＝（Ｙ（ω）・Ｘ^＊（ω））／（Ｘ（ω）・Ｘ^＊（ω））をＨ_１推定、Ｌ（ω）＝（Ｙ（ω）・Ｙ^＊（ω））／（Ｘ（ω）・Ｙ^＊（ω））をＨ_２推定と呼ぶが、どちらを用いても良い。なお、Ｘ^＊（ω）はＸ（ω）の複素共役、Ｙ^＊（ω）はＹ（ω）の複素共役である。

【0028】

同定処理部５０５は、得られた自己周波数応答関数Ｌ（ω）において、着目する近接固有値が存在する周波数帯域のみをズーミングする（図２のステップＳ１３０）。着目する近接固有値が存在する周波数帯域にはピークが現れる。そこで、同定処理部５０５は、自己周波数応答関数Ｌ（ω）におけるピークを検出することにより、着目する近接固有値が存在する周波数帯域を特定する。あるいは、同定処理部５０５は、自己周波数応答関数Ｌ（ω）をシステム同定装置５が備える表示装置に表示し、表示を確認したユーザが、ピークが現れる周波数帯域を入力してもよい。同定処理部５０５は、特定された周波数帯域の自己周波数応答関数Ｌ（ω）を抽出し、閾値よりも小さい値については０に置き換えるズーミングを行う。同定処理部５０５は、ズーミングした結果を逆ＦＦＴすることにより、インパルス応答関数ｇ_ｅ（ｔ）を得る。なお、ｔは時間を表す（図２のステップＳ１４０）。

【0029】

ここで、システム同定用モデルとして仮想二自由度モデルを仮想する。なお、仮想二自由度モデルは、対称型二自由度ばね質量系とも呼ばれる。
図７は、仮想二自由度モデルを示す図である。仮想二自由度モデルは、等しい質量、ばね定数、及び、減衰係数からなる一自由度ばね質量系が、ばね及びダッシュポットによって接続された系である。Ｍは質量、Ｋはばね定数、Ｃは減衰係数、Ｆは外力ベクトル、ｘ_１、ｘ_２は変位ベクトルである。また、Δ_Ｋはばね定数の変化、Δ_Ｃは減衰係数の変化を示す。仮想二自由度モデルは、運動方程式を表す質量行列、剛性行列、減衰行列が対称行列となるシステムであり、固有ベクトルが対称となる性質を有しており、近接固有値を有する系のシステム同定に好適となる。

【0030】

図７に示す仮想二自由度モデルの自己周波数応答関数Ｌ_１１を求めると、以下の式（１）となる。ｓは複素数を表す。ｓはｓ＝ｉ・ωで表される（ただし、ωは周波数、ｉは虚数単位）。

【0031】

【数1】

【0032】

式（１）をラブラス逆変換し、仮想二自由度モデルのインパルス応答関数ｇ_１１を求めると、式（２）となる。なお、ｓ＝ｉ・ωを式（１）に代入した場合、これを逆フーリエ変換すると、式（２）となる。

【0033】

【数2】

【0034】

δ（ｔ）はディラックのデルタ関数である。第１項の（１／Ｍ）・δ（ｔ）は、無視できる程度に他の項よりも小さな値である。
なお、ここで、各パラメータω_ｄ１、ω_ｄ２、α、β、γを式（３）とした。

【0035】

【数3】

【0036】

これより、未知パラメータは計５つとなる。実験値のインパルス応答関数ｇ_ｅ（ｔ）と式（２）の差の二乗和Ｊを目的関数とした多変数ニュートン法によって、パラメータ推定を行う更新式を求める。目的関数Ｊを式（４）に示し、更新式を式（５）に示す。ｉはサンプリングの番号、ｔ_ｉはｉ番目のサンプリングを行った時間を表す。

【0037】

【数4】

【0038】

【数5】

【0039】

ここで、λはステップ調整用パラメータであり、パラメータ推定アルゴリズムが発散する場合に０．００１〜０．１の間で調整すると、収束性が改善される。特に、管路系のシステム同定では、λを０．０１程度に設定することが好ましい。ｇ＾_１１は、現在のパラメータ値を用いて式（２）を算出したものである。図２のステップＳ１５０においては、各パラメータ（ω_ｄ１、ω_ｄ２、α、β、γ）の初期値、及び、ステップ幅となるλの値が入力される。なお、同定処理部５０５は、初期値の算出に用いられる情報の入力を受け、入力された情報に基づいて初期値を算出してもよい。

【0040】

同定処理部５０５は、ステップＳ１４０により得られたインパルス応答関数ｇ_ｅ（ｔ_ｉ）と、現在の各パラメータの値を用いて式（２）により算出したｇ＾_１１（ｔ_ｉ）とに基づいて、式（４）により二乗和Ｊを算出する。同定処理部５０５は、二乗和Ｊが閾値以下となるように、式（５）により各パラメータの値を更新する（図２のステップＳ１６０）。

【0041】

同定処理部５０５は、パラメータの更新を繰り返し、Ｊの値やその値の変化によって収束したか否かを判断する。同定処理部５０５は、収束しないと判断した場合（図２のステップＳ１７０：ＮＯ）、新たな初期値及びλの入力を受ける（ステップＳ１５０）。同定処理部５０５は、収束したと判断した場合、そのときのパラメータ値を用い、式（５）の関係に基づいて、質量Ｍ、ばね定数Ｋ及び減衰係数Ｃを算出する（ステップＳ１８０）。

【0042】

数値実験により、本実施形態によるシステム同定法の動作検証をおこなった。口径１００ｍｍのダクタイル鋳鉄管で構成された管路を模擬するため、真値がＭ＝１３．３０７０ｋｇ、Ｋ＝２．８９９４×１０^９Ｎ／ｍ、Ｃ＝１０００Ｎｓ／ｍ、Δ_Ｋ＝５．７９８９×１０^７Ｎ／ｍ、Δ_Ｃ＝６６６．６６６７Ｎｓ／ｍとした。本条件よりサンプリング周波数５０ｋＨｚで２０ｍｓのインパルス応答関数を生成し、さらに平均０、分散１の正規性白色雑音をインパルス応答関数に加算し、試験用のテストデータとした。初期値は各真値に０．９５を乗じた値を用い、ステップ調整用パラメータは０．０１を用いた。

【0043】

図８は、上記の条件によるシステム同定結果を示す図である。同図において横軸は周波数を表し、縦軸はアクセレランス（Accelerance）を表す。同図には、同定したパラメータを用いて計算した自己周波数応答関数（Identified）と、真値の自己周波数応答関数（Experiment）をそれぞれ示している。同図によれば、真値と、同定結果の両者の良い一致が確認できる。なお、推定値は、Ｍ＝１３．３０７３ｋｇ、Ｋ＝２．８９８０×１０^９Ｎ／ｍ、Ｃ＝９９９．９３１２Ｎｓ／ｍ、Δ_Ｋ＝６．０６５２×１０^７Ｎ／ｍ、Δ_Ｃ＝６６６．６７３０Ｎｓ／ｍであった。以上により、本実施形態による同定アルゴリズムの動作が確認できる。

【0044】

なお、第２の実施形態のように管路のシステム同定を行う場合、パラメータω_ｄ１、ω_ｄ２、α、β、γの初期値は以下のように算出される。これらの初期値を決める主たるパラメータである質量Ｍ及びばね定数Ｋは、以下の式（６）を用いて算出される。

【0045】

【数6】

【0046】

ここで、Ｒは半径、Ａは断面積である。管長をＬ、管厚をｈとした場合、Ａ＝ｈＬの関係がある。半径Ｒ、管厚ｈは管路管理台帳データから読み出した口径、管厚から求められる。管長Ｌは、管路管理台帳データから読み出してもよく、ユーザが入力してもよい。Ｅは管の弾性係数、Ｉは断面二次モーメントであり、Ｉ＝Ｌ・ｈ^３／１２である。ρは、管の密度である。弾性係数Ｅ及び管の密度ρは、管路管理台帳データから読み出した材種に応じた値としてもよく、予め決められた値としてもよい。Δ_Ｋはばね定数Ｋの１／１００、Δ_Ｃは減衰係数Ｃの１／３程度が望ましい。初期値の算出に用いる減衰係数Ｃ、減衰係数の変化Δ_Ｃは、管路管理台帳データから読み出した口径、管厚、材種のうち一以上に応じた値としてもよく、予め決められた値としてもよい。初期値設定部３０７は、これらの値を用いて式（３）により各パラメータω_ｄ１、ω_ｄ２、α、β、γの初期値を算出する。
以上が初期値設定部の具体的な表式である。

【0047】

［実験］
供用後の水道管を試験用管路に設置し、通水環境下でシステム同定実験を実施した。供試管路は口径１００ｍｍ、肉厚１０ｍｍの普通鋳鉄管で供用済みの管を用いた。管路の上部には消火栓を設置し、カプラ上に加速度センサを設置した。

【0048】

図９は、システム同定実験の結果を示す図である。同図において横軸は周波数を表し、縦軸はアクセレランス（Accelerance）を表す。また図中の丸のプロットは自己周波数応答関数の実験値（Experiment）を表し、実線は同定結果（Identified）を表している。同図に示すように、実験値と推定値のピーク位置やスペクトル形状について、良好な一致が確認された。なお、得られた推定値はＭ＝１４．８４４６ｋｇ、Ｋ＝３．４３４６×１０^９Ｎ／ｍ、Ｃ＝９０３．５４９１Ｎｓ／ｍ、Δ_Ｋ＝１．２８７７×１０^８Ｎ／ｍ、Δ_Ｃ＝９０３．５４９１Ｎｓ／ｍであった。

【0049】

関連技術と比較した優位性を示すため、特許文献１で用いられているＡＲモデルを用いた同定を行った。そして、ＡＲモデルを用いた同定方法と本実施形態による同定方法とを比較した。

【0050】

その結果、ＡＲ法による同定結果は、実験値との大きな乖離が確認され、同定困難であることがわかった。これは、インパルス応答関数がうなり波形をともなっており、時間領域同定法で用いる多項式モデルでは、本波形の特徴を記述するのに限界が生じていることを示している。以上の検討により、本実施形態が優れた効果を奏することがわかった。

【0051】

図１０は、本発明の実施形態によるシステム同定装置の最小構成を示すブロック図である。同図に示す最小構成のシステム同定装置１ａは、上述した第１の実施形態の分析部１０５を少なくとも備えればよい。分析部１０５は、分析対象を励振させた位置において計測された入力信号及び出力信号に基づいて自己周波数応答関数を算出する。そして、分析部１０５は、算出された自己周波数応答関数から得られるインパルス応答関数と、分析対象がモデル化された仮想二自由度モデルのインパルス応答関数とを用いて、分析対象のシステム同定を行う。

【0052】

本実施形態によれば、仮想二自由度モデルの自己周波数応答関数に対応するインパルス応答関数を用いて、近接固有値を有する系のシステム同定を行うことができる。

【0053】

上述した実施形態におけるシステム同定装置１、１ａ、３、５は、バスで接続されたＣＰＵ（Central Processing Unit）やメモリや補助記憶装置などを備え、システム同定プログラムを実行することによって上述した実施形態におけるシステム同定装置１、１ａ、３、５の一部の機能を実現する。なお、システム同定装置１、１ａ、３、５の一部の機能は、ＡＳＩＣ（Application Specific Integrated Circuit）やＰＬＤ（Programmable Logic Device）やＦＰＧＡ（Field Programmable Gate Array）等のハードウェアを用いて実現されても良い。システム同定プログラムは、コンピュータ読み取り可能な記録媒体に記録されても良い。コンピュータ読み取り可能な記録媒体とは、例えばフレキシブルディスク、光磁気ディスク、ＲＯＭ（Read Only Memory）、ＣＤ−ＲＯＭ（Compact Disc Read only memory）等の可搬媒体、コンピュータシステムに内蔵されるハードディスク等の記憶装置である。システム同定プログラムは、電気通信回線を介して送信されても良い。

【0054】

上記の実施形態の一部又は全部は、以下の付記のようにも記載されうるが、以下には限られない。

【0055】

（付記１）分析対象を励振させた位置において計測された入力信号及び出力信号に基づいて自己周波数応答関数を算出し、算出された前記自己周波数応答関数から得られるインパルス応答関数と、前記分析対象がモデル化された仮想二自由度モデルのインパルス応答関数とを用いて前記分析対象のシステム同定を行う分析部、を備えるシステム同定装置。

【0056】

（付記２）前記分析部は、多変数ニュートン法により仮想二自由度モデルの前記インパルス応答関数を推定し、推定により得られた前記インパルス応答関数に基づいて前記分析対象のシステム同定を行う、付記１に記載のシステム同定装置。

【0057】

（付記３）前記分析対象の物理データに基づいて多変数ニュートン法に用いられる初期値を算出する初期値設定部をさらに備える、付記２に記載のシステム同定装置。

【0058】

（付記４）前記分析対象は、管路である、付記１に記載のシステム同定装置。

【0059】

（付記５）前記分析対象を励振させる加振部と、前記加振部により前記分析対象を励振させた位置における入力信号及び出力信号を計測する計測部と、前記加振部が前記分析対象を励振させる位置と前記計測部が前記分析対象を計測する位置とを揃えるように前記加振部及び前記計測部を設置する設置位置決め部とをさらに備える、付記１に記載のシステム同定装置。

【0060】

（付記６）前記加振部は、力センサを内蔵したインパルスハンマー又は電磁式加振器である、付記５に記載のシステム同定装置。

【0061】

（付記７）前記計測部は、加速度ピックアップ、レーザー変位計、レーザードップラー型速度計又は接触式変位計である、付記５に記載のシステム同定装置。

【0062】

（付記８）前記設置位置決め部は、消火栓カプラである、付記５に記載のシステム同定装置。

【0063】

（付記９）前記分析部は、前記システム同定により前記分析対象の系の質量、剛性定数及び減衰係数を得る、付記１に記載のシステム同定装置。

【0064】

（付記１０）分析対象を励振させる加振ステップと、前記加振ステップにおいて前記分析対象を励振させた位置における入力信号及び出力信号を計測する計測ステップと、前記計測ステップにおいて計測された前記入力信号及び前記出力信号に基づいて自己周波数応答関数を算出し、算出された前記自己周波数応答関数から得られるインパルス応答関数と、前記分析対象がモデル化された仮想二自由度モデルのインパルス応答関数とを用いて前記分析対象のシステム同定を行う分析ステップと、を有するシステム同定方法。

【0065】

（付記１１）コンピュータに、分析対象を励振させた位置において計測された入力信号及び出力信号に基づいて自己周波数応答関数を算出し、算出された前記自己周波数応答関数から得られるインパルス応答関数と、前記分析対象がモデル化された仮想二自由度モデルのインパルス応答関数とを用いて前記分析対象のシステム同定を行う分析ステップ、を実行させるプログラム。
以上、実施形態（及び実施例）を参照して本願発明を説明したが、本願発明は上記実施形態（及び実施例）に限定されるものではない。本願発明の構成や詳細には、本願発明のスコープ内で当業者が理解し得る様々な変更をすることができる。
この出願は、２０１８年２月２１日に出願された日本出願特願２０１８−０２９２１８を基礎とする優先権を主張し、その開示の全てをここに取り込む。

【産業上の利用可能性】

【0066】

近接固有値を有する系のシステム同定に適用可能である。よってその工業的価値は高い。

【符号の説明】

【0067】

１、１ａ、３、５…システム同定装置
１０１、３０１…設置位置決め部
１０２、３０２…加振部
１０３、３０３…計測部
１０４、３０４…信号収集部
１０５、３０５…分析部
１０６、３０８…対象物理システム
３０６…記憶部
３０７…初期値設定部
５０１…消火栓カプラ
５０２…ハンマー
５０３…センサ
５０４…データロガー
５０５…同定処理部
５０６…管路

【図1】

【図2】

【図3】

【図4】

【図5】

【図6】

【図7】

【図8】

【図9】

【図10】