(57)【要約】【課題】多次元のパラメータセットの時系列データからなる移動軌跡データを適切に成形できる移動軌跡成形装置を提供する。【解決手段】移動軌跡成形装置は、多次元のパラメータセットの時系列データである移動軌跡データを取得する移動軌跡データ取得部と、予め規定された平滑化フィルタを適用して前記移動軌跡データを平滑化する平滑化処理部と、を備える。前記平滑化フィルタは、前記パラメータセットどうしの時間の差が大きくなるにつれて平滑化の影響度が小さくなる重み付けを規定する第１の重み付け関数と、前記パラメータセットどうしの多次元空間内の位置が離れるにつれて平滑化の影響度が小さくなる重み付けを規定する第２の重み付け関数と、の積からなる。【選択図】図２０
2019021036.tif 000002
【特許請求の範囲】
【請求項１】
多次元のパラメータセットの時系列データである移動軌跡データを取得する移動軌跡データ取得部と、
予め規定された平滑化フィルタを適用して前記移動軌跡データを平滑化する平滑化処理部と、
を備え、
前記平滑化フィルタは、
前記パラメータセットどうしの時間の差が大きくなるにつれて平滑化の影響度が小さくなる重み付けを規定する第１の重み付け関数と、前記パラメータセットどうしの多次元空間内の位置が離れるにつれて平滑化の影響度が小さくなる重み付けを規定する第２の重み付け関数と、の積からなる
移動軌跡成形装置。
【請求項２】
前記第１の重み付け関数及び前記第２の重み付け関数は、ガウス分布関数である
請求項１に記載の移動軌跡成形装置。
【請求項３】
前記多次元のパラメータセットを、次元別に正規化する正規化処理部を更に備える
請求項１又は請求項２に記載の移動軌跡成形装置。
【請求項４】
平滑化された前記移動軌跡データによって特定される移動軌跡曲線を、多次元空間内において等間隔で刻むとともに、当該等間隔で刻まれた前記移動軌跡曲線上の各点の座標のデータ系列を取得する等間隔化処理部
を更に備える請求項１から請求項３の何れか一項に記載の移動軌跡成形装置。
【請求項５】
前記多次元のパラメータセットは、
モータと、負荷と、前記モータ及び前記負荷を連結する連結部材と、を具備する機械系における、前記モータのトルク、前記モータの角速度、前記負荷の角度、前記負荷の角速度、前記連結部材のねじれ角度、及び、前記連結部材のねじれ角速度のうちの少なくとも何れか一つを含む
請求項１から請求項４の何れか一項に記載の移動軌跡成形装置。
【請求項６】
請求項１から請求項５の何れか一項に記載の移動軌跡成形装置と、
当該移動軌跡成形装置によって成形された移動軌跡データである評価対象移動軌跡データと予め用意された基準移動軌跡データとの対比に基づいて、前記移動軌跡データに異常があるか否かを判定する異常判定部と、
を備える異常診断システム。
【請求項７】
前記異常判定部は、
前記評価対象移動軌跡データと前記基準移動軌跡データとの間のＤＴＷ距離を演算し、算出された前記ＤＴＷ距離が所定の判定閾値を上回ったか否かに基づいて、前記移動軌跡データに異常があるか否かを判定する
請求項６に記載の異常診断システム。
【請求項８】
多次元のパラメータセットの時系列データである移動軌跡データを取得する移動軌跡データ取得ステップと、
予め規定された平滑化フィルタを適用して前記移動軌跡データを平滑化する平滑処理ステップと、
を有し、
前記平滑化フィルタは、
前記パラメータセットどうしの時間の差が大きくなるにつれて平滑化の影響度が小さくなる重み付けを規定する第１の重み付け関数と、前記パラメータセットどうしの多次元空間内の位置が離れるにつれて平滑化の影響度が小さくなる重み付けを規定する第２の重み付け関数と、の積からなる
移動軌跡成形方法。
【請求項９】
コンピュータに、
多次元のパラメータセットの時系列データである移動軌跡データを取得する移動軌跡データ取得ステップと、
予め規定された平滑化フィルタを適用して前記移動軌跡データを平滑化する平滑処理ステップと、
を実行させるプログラムであって、
前記平滑化フィルタは、
前記パラメータセットどうしの時間の差が大きくなるにつれて平滑化の影響度が小さくなる重み付けを規定する第１の重み付け関数と、前記パラメータセットどうしの多次元空間内の位置が離れるにつれて平滑化の影響度が小さくなる重み付けを規定する第２の重み付け関数と、の積からなる
プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、移動軌跡成形装置、異常診断システム、移動軌跡成形方法及びプログラムに関する。
【背景技術】
【０００２】
負荷を駆動させるモータの制御系（モータ制御システム）の設計において制御性能を向上させるためにはフィードバック制御では限界があり、指令信号に対する応答特性が高いフィードフォワード制御を行う必要がある。また、フィードフォワード制御を実現するためには、制御システム側が、制御対象（即ち、モータ及び負荷を有する機械系）を模したモデルの逆モデルを持つ必要があるが、このフィードフォワード制御用のモデルと実際の制御対象の特性との間に誤差があると、制御性能の向上は図れない。
【０００３】
したがって、制御性能の高いフィードフォワード制御を実現するためには、制御対象を正確に模したモデルを得る必要がある。ここで、当該モデルをなすモデルパラメータである制御対象全体の慣性モーメントＪ、粘性係数Ｄ及びクーロン摩擦係数τｆは、制御対象全体を一慣性系モデルとして模した場合に、所定のエラーシステムによって算出される誤差信号とその内部信号とに基づいて同定することが可能である。
【０００４】
また、特許文献１では、モデルパラメータの推定手段として、変位−変形量の位相面図（ロストモーション補正モデル）から、剛性と飽和時の摩擦の大きさを推定し、ロストモーション補正量を算出する方法が開示されている。
また、特許文献２では、取り付け時の測定、力測定器測定による力と変位の関係の測定結果から、各モデルパラメータを求め、ロストモーション補正のタイミングに使用する方法が開示されている。
【先行技術文献】
【特許文献】
【０００５】
【特許文献１】特開平８−１５２９１０号公報
【特許文献２】特開２００２−０２３８５２号公報
【発明の概要】
【発明が解決しようとする課題】
【０００６】
また、上述の制御対象（モータ及び負荷を有する機械系）の異常診断やヘルスモニタリングを行う目的で、当該制御対象の機械特性を表す複数種類のパラメータの時系列を組み合わせて分析することが想定される。このような多次元のパラメータセットの時系列データである移動軌跡データを精度良く分析できるようにするために、上記移動軌跡データを適切に成形する手法の確立が望まれている。
【０００７】
本発明の目的は、上記課題に鑑みてなされたものであって、多次元のパラメータセットの時系列データからなる移動軌跡データを適切に成形できる移動軌跡成形装置、異常診断システム、移動軌跡成形方法及びプログラムを提供することにある。
【課題を解決するための手段】
【０００８】
本発明の第１の態様によれば、移動軌跡成形装置は、多次元のパラメータセットの時系列データである移動軌跡データを取得する移動軌跡データ取得部と、予め規定された平滑化フィルタを適用して前記移動軌跡データを平滑化する平滑化処理部と、を備える。前記平滑化フィルタは、前記パラメータセットどうしの時間の差が大きくなるにつれて平滑化の影響度が小さくなる重み付けを規定する第１の重み付け関数と、前記パラメータセットどうしの多次元空間内の位置が離れるにつれて平滑化の影響度が小さくなる重み付けを規定する第２の重み付け関数と、の積からなる。
【０００９】
また、本発明の第２の態様によれば、前記第１の重み付け関数及び前記第２の重み付け関数は、ガウス分布関数である。
【００１０】
また、本発明の第３の態様によれば、上述の移動軌跡成形装置は、前記多次元のパラメータセットを、次元別に正規化する正規化処理部を更に備える。
【００１１】
また、本発明の第４の態様によれば、上述の移動軌跡成形装置は、平滑化された前記移動軌跡データによって特定される移動軌跡曲線を、多次元空間内において等間隔で刻むとともに、当該等間隔で刻まれた前記移動軌跡曲線上の各点の座標のデータ系列を取得する等間隔化処理部を更に備える。
【００１２】
また、本発明の第５の態様によれば、前記多次元のパラメータセットは、モータと、負荷と、前記モータ及び前記負荷を連結する連結部材と、を具備する機械系における、前記モータのトルク、前記モータの角速度、前記負荷の角度、前記負荷の角速度、前記連結部材のねじれ角度、及び、前記連結部材のねじれ角速度のうちの少なくとも何れか一つを含む。
【００１３】
また、本発明の第６の態様によれば、異常診断システムは、上述の移動軌跡成形装置と、当該移動軌跡成形装置によって成形された移動軌跡データである評価対象移動軌跡データと予め用意された基準移動軌跡データとの対比に基づいて、前記移動軌跡データに異常があるか否かを判定する異常判定部と、を備える。
【００１４】
また、本発明の第７の態様によれば、前記異常判定部は、前記評価対象移動軌跡データと前記基準移動軌跡データとの間のＤＴＷ距離を演算し、算出された前記ＤＴＷ距離が所定の判定閾値を上回ったか否かに基づいて、前記移動軌跡データに異常があるか否かを判定する。
【００１５】
また、本発明の第８の態様によれば、移動軌跡成形方法は、多次元のパラメータセットの時系列データである移動軌跡データを取得する移動軌跡データ取得ステップと、予め規定された平滑化フィルタを適用して前記移動軌跡データを平滑化する平滑処理ステップと、を有する。前記平滑化フィルタは、前記パラメータセットどうしの時間の差が大きくなるにつれて平滑化の影響度が小さくなる重み付けを規定する第１の重み付け関数と、前記パラメータセットどうしの多次元空間内の位置が離れるにつれて平滑化の影響度が小さくなる重み付けを規定する第２の重み付け関数と、の積からなる。
【００１６】
また、本発明の第９の態様によれば、プログラムは、コンピュータに、多次元のパラメータセットの時系列データである移動軌跡データを取得する移動軌跡データ取得ステップと、予め規定された平滑化フィルタを適用して前記移動軌跡データを平滑化する平滑処理ステップと、を実行させる。前記平滑化フィルタは、前記パラメータセットどうしの時間の差が大きくなるにつれて平滑化の影響度が小さくなる重み付けを規定する第１の重み付け関数と、前記パラメータセットどうしの多次元空間内の位置が離れるにつれて平滑化の影響度が小さくなる重み付けを規定する第２の重み付け関数と、の積からなる。
【発明の効果】
【００１７】
上述の移動軌跡成形装置、異常診断システム、移動軌跡成形方法及びプログラムによれば、多次元のパラメータセットの時系列データからなる移動軌跡データを適切に成形することができる。
【図面の簡単な説明】
【００１８】
【図１】第１の実施形態に係るモータ制御システムが制御対象とする機械系の概要を示す図である。
【図２】第１の実施形態に係る機械系のブロック線図を示す図である。
【図３】第１の実施形態に係る機械系のモータ側及び負荷側の摩擦特性を説明する図である。
【図４】第１の実施形態に係る機械系の連結部材における不感帯の特性を説明する図である。
【図５】第１の実施形態に係るモータ制御システムの機能構成を示す図である。
【図６】第１の実施形態に係るパラメータ同定装置の機能構成を示す図である。
【図７】第１の実施形態に係るパラメータ同定装置の処理フローを示す図である。
【図８】第１の実施形態に係るパラメータ同定装置の処理を説明する図である。
【図９】第１の実施形態に係る複数の位相面図を説明する第１の図である。
【図１０】第１の実施形態に係る複数の位相面図を説明する第２の図である。
【図１１】第１の実施形態に係る複数の位相面図を説明する第３の図である。
【図１２】第１の実施形態に係る複数の位相面図を説明する第４の図である。
【図１３】第１の実施形態に係る複数の位相面図を説明する第５の図である。
【図１４】第１の実施形態に係る複数の位相面図を説明する第６の図である。
【図１５】第２の実施形態に係るモータ制御システム、機械系、及び、異常診断システムの全体構成を示す図である。
【図１６】第２の実施形態に係る異常診断システムの機能構成を示す図である。
【図１７】第２の実施形態に係る異常診断システムの処理フローを示す図である。
【図１８】第２の実施形態に係る異常診断システムの各処理を詳細に説明するための第１の図である。
【図１９】第２の実施形態に係る異常診断システムの各処理を詳細に説明するための第２の図である。
【図２０】第２の実施形態に係る異常診断システムの各処理を詳細に説明するための第３の図である。
【図２１】第２の実施形態に係る異常診断システムの各処理を詳細に説明するための第４の図である。
【図２２】第２の実施形態に係る異常診断システムの各処理を詳細に説明するための第５の図である。
【発明を実施するための形態】
【００１９】
＜第１の実施形態＞
以下、第１の実施形態に係るモータ制御システムについて、図１〜図１４を参照しながら詳細に説明する。
【００２０】
（制御対象とする機械系の概要）
図１は、第１の実施形態に係るモータ制御システムが制御対象とする機械系の概要を示す図である。
図１に示すように、第１の実施形態に係るモータ制御システム１は、モータ２０と、負荷２１と、モータ２０及び負荷２１を機械的に連結する軸２２（連結部材）と、を有してなる機械系２の動作を制御するシステムである。ここで、本実施形態において、モータ２０は、内部に回転角度、回転速度を検出可能な回転検出器（エンコーダ）を有し、精密な位置決めを実現可能なサーボモータである。また、本実施形態において、軸２２は、例えば、ボールねじ、ギア等の剛性要素、ガタ要素を有する部材である。
【００２１】
図１に示すように、モータ制御システム１は、図示しない上位機器（ホストコントローラ）から負荷２１の目標回転角度（以下、「目標角度θｔ」と記載する。）を示す角度指令値を受け付ける。そして、モータ制御システム１は、受け付けた目標角度θｔに応じてモータ２０に発生すべきトルクτを算出し、当該トルクτを示すトルク指令値を出力する。モータ２０は、受け付けたトルク指令値に基づいてトルクτで回転駆動し、そのトルクτがモータ２０側から負荷２１側へと軸２２を通じて伝達する。その結果、伝達されたトルクτ’に応じて負荷２１が回転駆動する。このようにして、モータ制御システム１により、負荷２１の回転角度θＬが所望の目標角度θｔに一致するような制御が実現される。
【００２２】
本実施形態に係るモータ制御システム１は、いわゆるクローズドシステムであって、モータ２０の回転速度（角速度）と、負荷２１の回転角度と、を示す実測データを入力し、これらの実測データに基づいてフィードバック制御を行う。
なお、以下の説明において、モータ２０の回転速度（角速度）を「モータ角速度ωＭ」とも記載し、モータ２０の回転角度を「モータ角度θＭ」とも記載する。また、負荷２１の回転速度（角速度）を「負荷角速度ωＬ」とも記載し、負荷２１の回転角度を「負荷角度θＬ」とも記載する。
また、モータ制御システム１は、更に、機械系２を２慣性系と見なして模した２慣性系モデル（後述）を予め有しており、当該２慣性系モデルに基づくフィードフォワード制御を行う。
なお、機械系２における負荷２１が、回転系ではなく直動系の負荷（テーブル負荷等）の場合、負荷２１に対する制御対象パラメータは、厳密には“角度”、“角速度”ではなく、“位置”、“速度”となる。しかし、この場合、モータ制御システム１は、２慣性系として、モータ２０の“角度”、“角速度”と同じ次元で取り扱うために、負荷の“位置”、“速度”を、モータ軸換算値としての“角度”、“角速度”に逐次変換して各種制御を行うものとする。
【００２３】
機械系２を２慣性系と見なした場合において、機械系２固有の特性を表すパラメータ群は、モータ２０側の特性を示すパラメータ、負荷２１側の特性を示すパラメータ、及び、モータ２０と負荷２１とを機械的に連結して動力を伝達する軸２２（連結部材）の特性を示すパラメータに大別される。具体的には、モータ２０側の特性を示すパラメータとして、モータ側慣性モーメントＪＭ、モータ側粘性係数ＤＭ、モータ側摩擦係数τｆＭがある。また、負荷２１側の特性を示すパラメータとして、負荷側慣性モーメントＪＬ、負荷側粘性係数ＤＬ、負荷側摩擦係数τｆＬがある。更に、軸２２の特性を示すパラメータとして、ねじれ剛性係数ＫＲ、ねじれ粘性係数ＤＲ、不感帯幅ＢＬがある。
【００２４】
（２慣性系モデルの概要）
図２は、第１の実施形態に係る機械系のブロック線図を示す図である。
本実施形態に係る機械系２（図１参照）をモータ２０、負荷２１及び軸２２からなる２慣性系と見なすことで、当該機械系２における入力と出力の関係を、図２に示すような伝達関数を用いたブロック線図で表すことができる。
【００２５】
図２に示すように、モータ２０は、モータ制御システム１（図１）から受け付けたトルク指令値に基づいて発生させたトルクτを入力とし、モータ角速度ωＭ及びモータ角度θＭを出力とする。モータ２０に入力されたトルクτは、モータ２０に生じる摩擦であるモータ側摩擦係数τｆＭ、及び、軸２２を通じて負荷２１に伝達するトルクτ’が差し引かれた後、伝達要素１／（ＪＭｓ＋ＤＭ）を通じてモータ角速度ωＭに変換され、更に伝達要素１／ｓを通じてモータ角度θＭに変換される。
ここで、モータ側摩擦係数τｆＭは、モータ２０の速度反転時（モータ角速度ωＭの符号反転時）に符号のみが反転するクーロン摩擦成分だけでなく、モータ２０の速度反転後の変位（モータ角速度ωＭのゼロからの積分値）に依存して非線形に変化する非線形摩擦成分を含んでいる。モータ側摩擦係数τｆＭは、この非線形摩擦成分をモデル化して規定し、モータ角速度ωＭを入力変数とするモータ側摩擦関数ＧＭを通じて得られる。モータ側摩擦関数ＧＭの詳細については後述する。
【００２６】
同様に、負荷２１は、モータ２０から軸２２を通じて伝達されたトルクτ’（＝（ＫＲ＋ＤＲｓ）（θＭ−θＬ−ＢＫＬＳ））を入力とし、負荷角速度ωＬ及び負荷角度θＬを出力とする。負荷２１に入力されたトルクτ’は、負荷２１に生じる摩擦である負荷側摩擦係数τｆＬが差し引かれた後、伝達要素１／（ＪＬｓ＋ＤＬ）を通じて負荷角速度ωＬに変換され、更に伝達要素１／ｓを通じて負荷角度θＬに変換される。
ここで、負荷側摩擦係数τｆＬも同様に、負荷２１の速度反転後の変位（負荷角速度ωＬのゼロからの積分値）に依存して非線形に変化する非線形摩擦成分を含んでいる。負荷側摩擦係数τｆＬは、この非線形摩擦成分をモデル化して規定し、負荷角速度ωＬを入力変数とする負荷側摩擦関数ＧＬを通じて得られる。負荷側摩擦関数ＧＬの詳細については後述する。
【００２７】
なお、機械系２をなすモータ２０及び負荷２１には、それぞれ、図示しない回転検出器（エンコーダ）が設置されている。モータ制御システム１は、当該回転検出器を通じて、モータ角度θＭ、負荷角度θＬの経時的変化を示す実測データ、及び、これらの時間微分より導出可能なモータ角速度ωＭ、負荷角速度ωＬの経時的変化を示す実測データを取得する。
【００２８】
また、図２に示すように、軸２２は、モータ角度θＭと負荷角度θＬとの偏差（以下、「ねじれ角度（θＭ−θＬ）」とも記載する。）を入力とし、ねじれ角度（θＭ−θＬ）に応じたトルクτ’を出力とする伝達系である。入力されたねじれ角度（θＭ−θＬ）は、不感帯（ガタ）における非線形特性を表す不感帯特性関数Ｆ１と、伝達要素ＫＲ及び伝達要素ＤＲｓと、を通じて、モータ２０及び負荷２１に印加されるトルクτ’に変換される。
ここで、不感帯特性関数Ｆ１とは、軸２２のねじれ角度（θＭ−θＬ）を入力変数とし、当該ねじれ角度（θＭ−θＬ）からガタ変位ＢＫＬＳを差し引いたものを出力する非線形関数である。不感帯特性関数Ｆ１の詳細については後述する。
また、ねじれ剛性係数ＫＲとは、軸２２のねじれ方向についての剛性の度合いを示すパラメータであって、軸２２のばね定数に相当する。即ち、ねじれ剛性係数ＫＲは、負荷２１に印加されるトルクτ’のうち、軸２２のねじれ角度（θＭ−θＬ）に比例する成分を与える。
また、ねじれ粘性係数ＤＲとは、軸２２のねばりの度合いを示すパラメータであって、負荷２１に印加されるトルクτ’のうち、軸２２におけるモータ角速度ωＭと負荷角速度ωＬとの偏差（以下、「ねじれ速度（ωＭ−ωＬ）」とも記載する。）に比例する成分を与える。
【００２９】
図２に示す機械系２において、モータ側慣性モーメントＪＭ、モータ側粘性係数ＤＭ、モータ側摩擦係数τｆＭ、負荷側慣性モーメントＪＬ、負荷側粘性係数ＤＬ、負荷側摩擦係数τｆＬ、ねじれ剛性係数ＫＲ、ねじれ粘性係数ＤＲ及び不感帯幅ＢＬは、機械系２の現実の特性を表すパラメータ群であり、いずれも個別には観測が困難な未知のパラメータである。
一方、機械系２についての入力であるトルクτ、及び、出力であるモータ角速度ωＭ、モータ角度θＭ、負荷角速度ωＬ及び負荷角度θＬは、上述の回転検出器を通じて観測可能なパラメータである。
【００３０】
図３は、第１の実施形態に係る機械系のモータ側及び負荷側の摩擦特性を説明する図である。
図３（ａ）は、モータ側摩擦係数τｆＭ及び負荷側摩擦係数τｆＬの各々に生じる非線形摩擦特性を示すグラフである。
良く知られているクーロン摩擦は、物体の速度（モータ角速度ωＭ、負荷角速度ωＬ）の方向（正、負の符号）に依存してその方向（正、負の符号）のみが変化し、その量は、速度（ω）、変位（θ）に対しては変動しないものとして知られている。しかしながら、軸２２において、例えば、ボールねじ、ボールベアリング等の転がり要素が含まれる場合には、通常のクーロン摩擦とは特性が異なる“転がり摩擦”を考慮する必要がある。
ここで、転がり摩擦は、速度反転直後の転がり要素が転動しない“微動領域”においては、速度反転後の変位（モータ角速度ωＭ、負荷角速度ωＬのゼロからの積分値）によって見かけのばね定数が動的に変化する非線形ばね特性を有しており、図３（ａ）に示すようなヒステリシスカーブを描く。これは、転がり要素と軌道面の接触部における弾性変形やすべりによるものと考えられている。
また、転がり要素が有効に転動する“粗動領域”においては、クーロン摩擦により摩擦の速度に対する静的特性を示す。即ち、図３（ａ）に示すように、モータ２０の摩擦（モータ側摩擦係数τｆＭ）は、速度反転後の変位が所定以上となった時点で、モータ側クーロン摩擦係数τｆＭｃで飽和する。また、負荷２１の摩擦（負荷側摩擦係数τｆＬ）は、速度反転後の変位が所定以上となった時点で、負荷側クーロン摩擦係数τｆＬｃで飽和する。
【００３１】
本実施形態においては、モータ側摩擦特性関数ＧＭ及び負荷側摩擦特性関数ＧＬは、図３（ａ）に示す非線形摩擦特性を折れ線近似してなるＧＭＳモデルに基づいて規定される。ここで、ＧＭＳモデルとは、Ｎ個の特性の異なるブロック（転がり要素）とばねとが並列に接続されている状態を模したものである。図３（ｂ）に示すグラフ（折れ線）を構成する複数の直線の各々が、Ｎ個のブロックのうちのいくつまでが転がりきっている（有効に転動している）状態か、を表現している。
図３（ｂ）に示すように、モータ側摩擦特性関数ＧＭは、変位０〜θ１、θ１〜θ２、θ２〜θ３、θ３〜θ４の各々において採用すべき折れ線を規定する傾きＫ１、Ｋ２、Ｋ３、Ｋ４、及び、オフセット（変位０における切片）τｆＭ１（＝−τｆＭｃ）、τｆＭ２、τｆＭ３、τｆＭ４、τｆＭ５（＝τｆＭｃ）によって規定される。
なお、図示を省略するが、負荷側摩擦特性関数ＧＬも、図３（ｂ）に示すモータ側摩擦特性関数ＧＭと同様に規定される。
【００３２】
図４は、第１の実施形態に係る機械系の連結部材における不感帯の特性を説明する図である。
【００３３】
図４（ａ）に示す不感帯特性関数Ｆ１は、モータ２０と負荷２１との間に設けられた軸２２の“遊び”（不感帯）の幅を示す不感帯幅ＢＬに基づく伝達特性を表している。
不感帯特性関数Ｆ１は、不感帯幅ＢＬに依存する非線形特性であって、軸２２のねじれ角度（θＭ−θＬ）を入力変数とする。図４（ａ）に示すように、軸２２のねじれ角度（θＭ−θＬ）の絶対値が不感帯幅ＢＬ以下の場合（−ＢＬ≦θＭ−θＬ≦＋ＢＬ）、ガタ出力はゼロとなる。即ち、この場合、モータ２０から負荷２１へトルクが伝達されない。
他方、軸２２のねじれ角度（θＭ−θＬ）の絶対値が不感帯幅ＢＬよりも大きい場合（θＭ−θＬ＜−ＢＬ，θＭ−θＬ＞＋ＢＬ）、ガタ出力は、ねじれ角度（θＭ−θＬ）から不感帯幅ＢＬだけ小さい値（θＭ−θＬ＋ＢＬ（θＭ−θＬ＜−ＢＬ）、又は、θＭ−θＬ−ＢＬ（θＭ−θＬ＞＋ＢＬ））となる。
【００３４】
また、図４（ｂ）に示すガタ変位関数Ｆ１’は、ガタ変位ＢＫＬＳとねじれ角度（θＭ−θＬ）との関係をしている。
ガタ変位ＢＫＬＳは、不感帯幅ＢＬ（−ＢＬ〜＋ＢＬ）の幅を有する不感帯における変位量を示すパラメータであって、ねじれ角度（θＭ−θＬ）に対し、図４（ｂ）に示すような特性を有している。即ち、図４（ａ）に示すガタ出力は、ガタ変位ＢＫＬＳを用いて“θＭ−θＬ―ＢＫＬＳ”と表すことができる。
図４（ｂ）に示すように、ガタ変位ＢＫＬＳは、ねじれ角度（θＭ−θＬ）が不感帯の最小値（−ＢＬ）よりも小さい範囲では、ガタ変位ＢＫＬＳは最小値（−ＢＬ）をとり、ねじれ角度（θＭ−θＬ）が不感帯の最大値（＋ＢＬ）よりも大きい範囲では、ガタ変位ＢＫＬＳは最大値（＋ＢＬ）をとる。また、ねじれ角度（θＭ−θＬ）が不感帯の最小値（−ＢＬ）以上かつ最大値（＋ＢＬ）以下の範囲においては、ガタ変位ＢＫＬＳは、ねじれ角度（θＭ−θＬ）と同一の値を有する特性となる。
【００３５】
（機能構成）
図５は、第１の実施形態に係るモータ制御システムの機能構成を示す図である。
図５に示すように、第１の実施形態に係るモータ制御システム１は、フィードバック制御部１０と、フィードフォワード制御部１１と、パラメータ同定装置１２と、を備えている。
【００３６】
フィードバック制御部１０は、回転検出器を通じて観測されるモータ角速度ωＭと、角度指令値により指定される負荷２１の目標角度θｔ及び観測される現在の負荷角度θＬの偏差（θｔ−θＬ）と、に基づいて、モータ２０の制御を行う。
具体的には、フィードバック制御部１０は、目標角度θｔと負荷角度θＬとの偏差（θｔ−θＬ）をゼロとするためのトルクを算出し、その算出結果を示すトルク指令値を出力する。その際、フィードバック制御部１０は、回転検出器を通じて観測されたモータ角速度ωＭを参照して、適切かつ迅速なフィードバック制御がなされるようなトルクを算出する。
【００３７】
フィードフォワード制御部１１は、機械系２を模した２慣性系モデルＭＯＤを内部に有している。この２慣性系モデルＭＯＤは、複数のモデルパラメータ群からなり、図２に示す機械系２の逆モデルとなるように構成されている。
具体的には、複数のモデルパラメータ群とは、モデルモータ側慣性モーメントＪＭ０、モデルモータ側粘性係数ＤＭ０、モデルモータ側摩擦係数τｆＭ０、モデル負荷側慣性モーメントＪＬ０、モデル負荷側粘性係数ＤＬ０、モデル負荷側摩擦係数τｆＬ０、モデル不感帯幅ＢＬ０、及び、モデルねじれ剛性係数ＫＲ０である。
【００３８】
２慣性系モデルＭＯＤは、機械系２を模して構築された逆モデルであり、当該２慣性系モデルＭＯＤをなす上記モデルパラメータ群（ＪＭ０、ＤＭ０、τｆＭ０、ＪＬ０、ＤＬ０、τｆＬ０、ＢＬ０、ＫＲ０）は、機械系２の実際の特性を表す未知のパラメータ群（ＪＭ、ＤＭ、τｆＭ、ＪＬ、ＤＬ、τｆＬ、ＢＬ、ＫＲ）の各々に対応している。
【００３９】
ここで、モデルパラメータ群（ＪＭ０、ＤＭ０、τｆＭ０、ＪＬ０、ＤＬ０、τｆＬ０、ＢＬ０、ＫＲ０）と、機械系２の実際の特性を示すパラメータ群（ＪＭ、ＤＭ、τｆＭ、ＪＬ、ＤＬ、τｆＬ、ＢＬ、ＫＲ）との間に誤差が全くないと仮定すると、２慣性系モデルＭＯＤは、現実の機械系２の逆モデルに完全に一致するものとなる。そうすると、フィードフォワード制御部１１が目標角度θｔ及び２慣性系モデルＭＯＤに基づいて算出したトルクを機械系２に出力して駆動した結果、機械系２を通じて実際に得られる負荷角度θＬは、完全に目標角度θｔに一致するはずである（θＬ＝θｔ）。このように、フィードフォワード制御部１１は、予め機械系２を模して規定された２慣性系モデルＭＯＤに基づいて、機械系２固有の特性を織り込みながらトルクを算出することで、フィードバック制御よりも応答性が高いフィードフォワード制御において、より高精度な制御を実現することができる。
【００４０】
なお、図５に示すように、モータ制御システム１は、フィードバック制御部１０により算出されたトルクと、フィードフォワード制御部１１により算出されたトルクと、を加算して得られるトルクτ（トルク指令値）を機械系２（モータ２０（図１、図２））に向けて出力する。
【００４１】
以上のように、本実施形態に係るモータ制御システム１は、モータ２０、負荷２１及び軸２２を具備する機械系２を模した２慣性系モデルＭＯＤをなす複数のモデルパラメータ群（ＪＭ０、ＤＭ０、τｆＭ０、ＪＬ０、ＤＬ０、τｆＬ０、ＢＬ０、ＫＲ０）に基づいて、モータ２０に発生させるべきトルクτを算出する。
【００４２】
フィードフォワード制御部１１において、モデルパラメータ群（ＪＭ０、ＤＭ０、τｆＭ０、ＪＬ０、ＤＬ０、τｆＬ０、ＢＬ０、ＫＲ０）は、予め、機械系２の特性を表すパラメータとして想定される値で規定されている。しかしながら、実際には、モデルパラメータ群（ＪＭ０、ＤＭ０、τｆＭ０、ＪＬ０、ＤＬ０、τｆＬ０、ＢＬ０、ＫＲ０）と、機械系２固有のパラメータ群（ＪＭ、ＤＭ、τｆＭ、ＪＬ、ＤＬ、τｆＬ、ＢＬ、ＫＲ）との間には、それぞれ誤差が存在している。このような誤差があると、精度の高いフィードフォワード制御を実現することができない。
そこで、本実施形態に係るモータ制御システム１は、上述の複数のモデルパラメータ群のうち、特に、モデルモータ側摩擦係数τｆＭ０、モデル負荷側摩擦係数τｆＬ０、モデル不感帯幅ＢＬ０及びモデルねじれ剛性係数ＫＲ０を別個独立に、精度良く同定可能なパラメータ同定装置１２を備えている。
【００４３】
図５に示すように、パラメータ同定装置１２は、機械系２におけるモータ角度θＭ、モータ角速度ωＭ、負荷角度θＬ及び負荷角速度ωＬの経時的変化を示す実測データを取得する。パラメータ同定装置１２は、上記４個の実測データに基づいて、モデルモータ側摩擦係数τｆＭ０、モデル負荷側摩擦係数τｆＬ０、モデル不感帯幅ＢＬ０及びモデルねじれ剛性係数ＫＲ０を同定する。
以下、パラメータ同定装置１２の機能について詳細に説明する。
【００４４】
図６は、第１の実施形態に係るパラメータ同定装置の機能構成を示す図である。
図６に示すように、パラメータ同定装置１２は、データ取得部１２０と、位相面図生成部１２１と、パラメータ同定部１２２と、を備えている。
【００４５】
データ取得部１２０は、モータ２０に対するトルク指令値（トルクτ）と、モータ２０の角度（モータ角度θＭ）及び角速度（モータ角速度ωＭ）の経時的変化を示す実測データと、負荷２１の角度（負荷角度θＬ）及び角速度（負荷角速度ωＬ）の経時的変化を示す実測データと、を取得する。
【００４６】
位相面図生成部１２１は、データ取得部１２０によって取得されたトルク指令値（トルクτ）と、モータ２０の角度及び角速度の実測データと、負荷２１の角度及び角速度の実測データと、に基づいて複数の位相面図を生成する。
ここで、図６に示すように、位相面図生成部１２１は、連結部材（軸２２）のねじれ角度（θＭ−θＬ）とトルク指令値（トルクτ）との関係を示す「ねじれ角度−トルク位相面図Ｐ１」を生成する。
また、位相面図生成部１２１は、トルク指令値（トルクτ）とモータ２０の角速度（モータ角速度ωＭ）及び負荷２１の角速度（負荷角速度ωＬ）との関係を示す「トルク−機械角速度位相面図Ｐ２」を生成する。
また、位相面図生成部１２１は、軸２２のねじれ角度（θＭ−θＬ）とモータ２０の角速度（モータ角速度ωＭ）及び負荷２１の角速度（負荷角速度ωＬ）との関係を示す「ねじれ角度−機械角速度位相面図Ｐ３」を生成する。
また、位相面図生成部１２１は、トルク指令値（トルクτ）と軸２２のねじれ速度（モータ角速度ωＭと負荷角速度ωＬとの差。以下、“ωＭ−ωＬ”と記載する。）との関係を示す「トルク−ねじれ速度位相面図Ｐ４」を生成する。
また、位相面図生成部１２１は、軸２２のねじれ角度（θＭ−θＬ）と軸２２のねじれ速度（ωＭ−ωＬ）との関係を示す「ねじれ角度−ねじれ速度位相面図Ｐ５」を生成する。
また、位相面図生成部１２１は、負荷２１の角度（負荷角度θＬ）とトルク指令値（トルクτ）との関係を示す「負荷角度−トルク位相面図Ｐ６」を生成する。
【００４７】
パラメータ同定部１２２は、位相面図生成部１２１によって生成された複数の位相面図に基づいて、モータ２０の摩擦（モータ側摩擦係数τｆＭ）、負荷２１の摩擦（負荷側摩擦係数τｆＬ）、軸２２の剛性（ねじれ剛性係数ＫＲ）、及び、軸２２の不感帯幅ＢＬを示すモデルパラメータ（τｆＭ０、τｆＬ０、ＫＲ０、ＢＬ０）を同定する。
【００４８】
（パラメータ同定装置の処理フロー）
図７は、第１の実施形態に係るパラメータ同定装置の処理フローを示す図である。
また、図８は、第１の実施形態に係るパラメータ同定装置の処理を説明する図である。
また、図９〜図１４は、第１の実施形態に係る複数の位相面図を説明する第１から第６の図である。
図７に示す処理フローは、例えば、機械系２の実運転の開始前等において、パラメータ同定装置１２（図６）が機械系２についてのパラメータ同定を行う際に実行される。
【００４９】
まず、モータ制御システム１のオペレータは、上位機器を操作して所定の角度指令値を出力させることで、負荷２１の角度（負荷角度θＬ）を小振幅かつ低周波数で反復動作させる（ステップＳ０１）。
【００５０】
次に、パラメータ同定装置１２のデータ取得部１２０は、小振幅かつ低周波数で反復動作中の機械系２から、モータ角速度ωＭ、モータ角度θＭ、負荷角速度ωＬ及び負荷角度θＬの経時的変化を示す実測データを取得する（ステップＳ０２）。
【００５１】
次に、パラメータ同定装置１２の位相面図生成部１２１は、ステップＳ０２で取得された実測データから複数の位相面図を生成する（ステップＳ０３）。
【００５２】
ここで、位相面図生成部１２１は、図８に示すように、負荷角度θＬの方向反転前後に取得されたねじれ角度（θＭ−θＬ）の実測データ、及び、トルクτの実測データから同一時刻にサンプリングされた実測値どうしを紐付けて抽出する。ここで、「ねじれ角度（θＭ−θＬ）の実測データ」は、モータ角度θＭの実測データから負荷角度θＬの実測データを差し引くことで得られる。そして、位相面図生成部１２１は、抽出した実測値に基づいて、図９に示すような、ねじれ角度（θＭ−θＬ）を横軸にとり、トルクτを縦軸にとる位相面図である「ねじれ角度−トルク位相面図Ｐ１」を生成する。
なお、位相面図生成部１２１は、負荷角度θＬの正方向から負方向への反転前後に取得される実測データと、負方向から正方向への反転前後に取得される実測データと、の各々について上記実測値の紐付けを行う。このようにして生成されたねじれ角度−トルク位相面図Ｐ１には、主に、機械系２の非線形摩擦特性（図３（ａ）参照）と不感帯特性関数Ｆ１に起因するヒステリシス曲線が表れる。
【００５３】
同様に、位相面図生成部１２１は、負荷角度θＬの方向反転前後におけるトルクτの実測データ、及び、モータ角速度ωＭ、負荷角速度ωＬの実測データから同一時刻にサンプリングされた実測値どうしを紐付けて抽出する。そして、位相面図生成部１２１は、抽出した実測値に基づいて、図１０に示すような、トルクτを横軸にとり、モータ角速度ωＭ及び負荷角速度ωＬを縦軸にとる位相面図である「トルク−機械角速度位相面図Ｐ２」を生成する。
同様に、位相面図生成部１２１は、負荷角度θＬの方向反転前後におけるねじれ角度（θＭ−θＬ）の実測データ、及び、モータ角速度ωＭ、負荷角速度ωＬの実測データから同一時刻にサンプリングされた実測値どうしを紐付けて抽出する。そして、位相面図生成部１２１は、抽出した実測値に基づいて、図１１に示すような、ねじれ角度（θＭ−θＬ）を横軸にとり、モータ角速度ωＭ及び負荷角速度ωＬを縦軸にとる位相面図である「ねじれ角度−機械角速度位相面図Ｐ３」を生成する。
更に、位相面図生成部１２１は、負荷角度θＬの方向反転前後におけるトルクτの実測データ、及び、ねじれ速度（ωＭ−ωＬ）の実測データから同一時刻にサンプリングされた実測値どうしを紐付けて抽出する。ここで、「ねじれ速度（ωＭ−ωＬ）の実測データ」は、モータ角速度ωＭの実測データから負荷角速度ωＬの実測データを差し引いて得られる。そして、位相面図生成部１２１は、抽出した実測値に基づいて、図１２に示すような、トルクτを横軸にとり、ねじれ速度（ωＭ−ωＬ）を縦軸にとる位相面図である「トルク−ねじれ速度位相面図Ｐ４」を生成する。
更に、位相面図生成部１２１は、負荷角度θＬの方向反転前後におけるねじれ角度（θＭ−θＬ）の実測データ、及び、ねじれ速度（ωＭ−ωＬ）の実測データから同一時刻にサンプリングされた実測値どうしを紐付けて抽出する。そして、位相面図生成部１２１は、抽出した実測値に基づいて、図１３に示すような、ねじれ角度（θＭ−θＬ）を横軸にとり、ねじれ速度（ωＭ−ωＬ）を縦軸にとる位相面図である「ねじれ角度−ねじれ速度位相面図Ｐ５」を生成する。
更に、位相面図生成部１２１は、負荷角度θＬの方向反転前後における負荷角度θＬの実測データ、及び、トルクτの実測データから同一時刻にサンプリングされた実測値どうしを紐付けて抽出する。そして、位相面図生成部１２１は、抽出した実測値に基づいて、図１４に示すような、負荷角度θＬを横軸にとり、トルクτを縦軸にとる位相面図である「負荷角度−トルク位相面図Ｐ６」を生成する。
【００５４】
次に、パラメータ同定部１２２は、ステップＳ０３で生成された各種位相面図に基づいて、２慣性系モデルＭＯＤをなす各種モデルパラメータのうち、モデルモータ側摩擦係数τｆＭ０、モデル負荷側摩擦係数τｆＬ０、モデル不感帯幅ＢＬ０及びモデルねじれ剛性係数ＫＲ０を同定する（ステップＳ０４）。具体的には、パラメータ同定部１２２は、各種位相面図における折れ点の位置等に基づいて上記モデルパラメータを同定する。
【００５５】
ここで、図２に示すブロック線図に基づいて、モータ２０についての運動方程式である式（１）、及び、負荷２１についての運動方程式である式（２）を得ることができる。
【００５６】
【数１】
2019021036.tif 000003
【００５７】
【数２】
2019021036.tif 000004
【００５８】
式（１）、（２）において、“τｆＭ［ωＭ］”は、モータ側摩擦特性関数ＧＭ（図３（ｂ））に基づいてモータ角速度ωＭの変数として規定されるモータ側摩擦係数τｆＭを示している。また、“τｆＬ［ωＬ］”は、負荷側摩擦特性関数ＧＬに基づいて負荷角速度ωＬの変数として規定される負荷側摩擦係数τｆＬを示している。
また、“ＢＫＬＳ［θＭ−θＬ］”は、ガタ変位関数Ｆ１’（図４（ｂ））に基づいてねじれ角度（θＭ−θＬ）の変数として規定されるガタ変位ＢＫＬＳを示している。
【００５９】
なお、現実の機械系２において、軸２２のねじれ粘性係数ＤＲの影響は小さいため、式（１）、式（２）においては、ねじれ粘性係数ＤＲを無視している。
また、上記運動方程式（式（１）、式（２））において、機械系２の反復動作を小振幅かつ低周波数にすると、モータ２０の慣性トルクＪＭｓ２θＭ、負荷２１の慣性トルクＪＬｓ２θＬ、モータ２０の粘性トルクＤＭｓθＭ、及び、負荷２１の粘性トルクＤＬｓθＬの影響も小さくなる。そこで、上記運動方程式におけるこれらの項を削除することで、以下に示す式（３）、式（４）が導出される。
【００６０】
【数３】
2019021036.tif 000005
【００６１】
【数４】
2019021036.tif 000006
【００６２】
式（３）より、機械系２の反復動作が小振幅かつ低周波数の条件下において、トルク指令値に従ってモータ２０が出力するトルクτは、モータ側摩擦係数τｆＭと負荷側摩擦係数τｆＬとの合計値となることが示される。
また、式（４）より、機械系２の反復動作が小振幅かつ低周波数の条件下において、負荷側摩擦係数τｆＬをねじれ剛性係数ＫＲで除した値とガタ変位ＢＫＬＳとの和がねじれ角度（θＭ−θＬ）となることが示される。即ち、軸２２におけるねじれ角度（θＭ−θＬ）からガタ変位ＢＫＬＳを差し引いた角度（変位）に、軸２２のばね定数（ねじれ剛性係数ＫＲ）を乗じた値が、負荷側摩擦係数τｆＬに一致する。
【００６３】
以下、式（３）、（４）の関係を考慮しながら、図９に示すねじれ角度−トルク位相面図Ｐ１における各折れ点Ｔ１１〜Ｔ１４について説明する。
【００６４】
ねじれ角度−トルク位相面図Ｐ１における折れ点Ｔ１１〜折れ点Ｔ１２は、負方向に飽和していたトルクτが正方向に向けて反転し、これに応じてモータ２０が正方向への回転駆動を開始するまでの軌跡を示している。即ち、折れ点Ｔ１１において最小値（負方向に飽和状態）となっていたトルクτが正方向に向けて徐々に増加していく過程において、当該トルクτがモータ２０に生じる摩擦の飽和値（モータ側クーロン摩擦係数τｆＭｃ）を上回った時点でモータ２０が正方向に回転駆動し始める。他方、折れ点Ｔ１２の時点では、不感帯（−ＢＬ〜＋ＢＬ）の存在により負荷２１の摩擦の影響は生じない。したがって、折れ点Ｔ１２に至った時点におけるトルクτは、式（３）より、モータ２０に生じる摩擦の飽和値（モータ側クーロン摩擦係数τｆＭｃ）に一致する。
【００６５】
また、ねじれ角度−トルク位相面図Ｐ１における折れ点Ｔ１２〜折れ点Ｔ１３は、軸２２のねじれ角度（θＭ−θＬ）が不感帯の最小値（−ＢＬ）から最大値（＋ＢＬ）に至るまでの軌跡を示している。即ち、折れ点Ｔ１２から折れ点Ｔ１３にかけては、ガタ変位ＢＫＬＳのみが変化し、負荷２１には動力が伝達されない。したがって、モータ２０のみが駆動するため、トルクτがモータ側クーロン摩擦係数τｆＭｃ一定を維持しながら、ねじれ角度（θＭ−θＬ）が変化する。したがって、折れ点Ｔ１３に至った時点におけるねじれ角度（θＭ−θＬ）は、式（４）より、不感帯の最大値（＋ＢＬ）に一致する。
【００６６】
また、ねじれ角度−トルク位相面図Ｐ１における折れ点Ｔ１３〜折れ点Ｔ１４は、ガタ変位ＢＫＬＳが正方向に飽和して、不感帯幅ＢＬを上回るねじれ角度（θＭ−θＬ）を通じて、負荷２１に印加されるトルクτ’が徐々に増加し、これに応じて負荷２１が正方向への回転駆動を開始するまでの軌跡を示している。即ち、折れ点Ｔ１３においてモータ２０のみを回転駆動させていたトルクτが正方向に向けて更に増加していく過程において、当該トルクτが、モータ２０に生じる摩擦の飽和値と負荷２１に生じる摩擦の飽和値（負荷側クーロン摩擦係数τｆＬｃ）との合計値と釣り合った時点で負荷２１が正方向に回転駆動し始める。したがって、折れ点Ｔ１４に至った時点におけるトルクτは、式（３）より、モータ側クーロン摩擦係数τｆＭｃと負荷側クーロン摩擦係数τｆＬｃとの合計値に一致する。
また、折れ点Ｔ１４に至った時点におけるねじれ角度（θＭ−θＬ）は、式（４）より、“τｆＬｃ／ＫＲ＋ＢＬ”に一致する。
【００６７】
以上の特性に基づき、パラメータ同定部１２２は、各折れ点Ｔ１１〜Ｔ１４の位置に基づいて各種モデルパラメータを同定する。
具体的には、パラメータ同定部１２２は、ねじれ角度−トルク位相面図Ｐ１における折れ点Ｔ１２（又は折れ点Ｔ１３）の縦軸の値を参照してモデルモータ側クーロン摩擦係数τｆＭｃ０を同定する。また、パラメータ同定部１２２は、折れ点Ｔ１４における縦軸の値から、先に同定したモデルモータ側クーロン摩擦係数τｆＭｃ０の値を差し引いてモデル負荷側クーロン摩擦係数τｆＬｃ０を同定する。更に、パラメータ同定部１２２は、折れ点Ｔ１３（又は折れ点Ｔ１２）の横軸の値を参照してモデル不感帯幅ＢＬ０を同定する。更に、パラメータ同定部１２２は、折れ点Ｔ１４における横軸の値と、先に同定したモデル不感帯幅ＢＬ０及びモデル負荷側クーロン摩擦係数τｆＬｃ０とに基づいて式（４）を解くことで、モデルねじれ剛性係数ＫＲ０を同定する。
【００６８】
次に、式（３）、（４）の関係を考慮しながら、図１０に示すトルク−機械角速度位相面図Ｐ２における各折れ点Ｔ２１〜Ｔ２４について説明する。
トルク−機械角速度位相面図Ｐ２における折れ点Ｔ２１〜折れ点Ｔ２２は、負方向に飽和していたトルクτが正方向に向けて反転し、これに応じて、トルクτが最小値（負の飽和状態）からほとんど変化しないまま、モータ２０及び負荷２１の角速度がゼロとなるまでの軌跡を示している。
また、折れ点Ｔ２２〜折れ点Ｔ２３は、モータ２０及び負荷２１の角速度がゼロの状態から徐々に正方向にトルクτが上昇して行く中で、不感帯（−ＢＬ〜＋ＢＬ）の存在により、モータ角速度ωＭのみが上昇する軌跡を示している。
また、折れ点Ｔ２３〜折れ点Ｔ２４は、トルクτが、モータ２０のみの駆動に要するトルク（モータ側クーロン摩擦係数τｆＭｃ）を上回ることで負荷２１にトルクτ’が印加され、モータ角速度ωＭと負荷角速度ωＬとが等しく上昇する軌跡を示している。
【００６９】
ここで、折れ点Ｔ２２の横軸の値は、モータ角速度ωＭ、負荷角速度ωＬがゼロであることに基づき、式（３）より“−τｆＭｃ−τｆＬｃ”と求めることができる。
また、折れ点Ｔ２３の横軸の値は、負荷２１の角速度が上昇を開始する時点（モータ２０のみが回転駆動している時点）におけるトルクである。したがって、式（３）より、折れ点Ｔ２３の横軸の値は、モータ側クーロン摩擦係数τｆＭｃに一致する。
【００７０】
パラメータ同定部１２２は、トルク−機械角速度位相面図Ｐ２における折れ点Ｔ２２の横軸の値から、先に同定したモデルモータ側クーロン摩擦係数τｆＭｃ０の値を差し引いてモデル負荷側クーロン摩擦係数τｆＬｃ０を同定する。
ここで、ねじれ角度−トルク位相面図Ｐ１の折れ点Ｔ１４では、トルクτ、ねじれ角度（θＭ−θＬ）が完全には飽和せず、モータ側粘性係数ＤＭ及び負荷側粘性係数ＤＬに基づく誤差成分が生じることが分かっている。そうすると、ねじれ角度−トルク位相面図Ｐ１の折れ点Ｔ１４の位置に基づいてモデル負荷側クーロン摩擦係数τｆＬｃ０を同定しようとすると、当該誤差成分に起因して精度が低下することが考えられる。
他方、トルク−機械角速度位相面図Ｐ２の折れ点Ｔ２２では、モータ角速度ωＭ、負荷角速度ωＬがゼロとなっている。したがって、折れ点Ｔ２２では、モータ角速度ωＭ、負荷角速度ωＬに比例する成分を与えるモータ側粘性係数ＤＭ及び負荷側粘性係数ＤＬの影響を受けていない。そのため、折れ点Ｔ２２の横軸の値から、より精度の高いモデル負荷側クーロン摩擦係数τｆＬｃ０を同定することができる。
【００７１】
また、パラメータ同定部１２２は、トルク−機械角速度位相面図Ｐ２における折れ点Ｔ２３の横軸の値を参照してモデルモータ側クーロン摩擦係数τｆＭｃ０を同定する。
【００７２】
次に、式（３）、（４）の関係を考慮しながら、図１１に示すねじれ角度−機械角速度位相面図Ｐ３における各折れ点Ｔ３１〜Ｔ３４について説明する。
ねじれ角度−機械角速度位相面図Ｐ３における折れ点Ｔ３１〜折れ点Ｔ３２は、負方向に飽和していたトルクτが正方向に向けて反転し、これに応じて、ねじれ角度（θＭ−θＬ）が最小値（負の飽和状態）からほとんど変化しないまま、モータ２０及び負荷２１の角速度がゼロとなるまでの軌跡を示している。
また、折れ点Ｔ３２〜折れ点Ｔ３３は、不感帯（−ＢＬ〜＋ＢＬ）の存在により負荷角速度ωＬがゼロのままモータ角速度ωＭのみが上昇することで、ねじれ角度（θＭ−θＬ）が徐々に増大していく軌跡を示している。
また、折れ点Ｔ３３〜折れ点Ｔ３４は、ねじれ角度（θＭ−θＬ）が不感帯の最大値（＋ＢＬ）を上回ることで負荷２１にトルクτ’が印加され、モータ角速度ωＭと負荷角速度ωＬとが等しく上昇する軌跡を示している。
【００７３】
ここで、折れ点Ｔ３２の横軸の値は、負荷角速度ωＬがゼロであること、及び、ガタ変位ＢＫＬＳが“−ＢＬ”（θＭ−θＬ＜−ＢＬ）であることに基づき、式（４）より“−τｆＬｃ／ＫＲ−ＢＬ”と求めることができる。
また、折れ点Ｔ３３の横軸の値は、負荷２１の角速度が上昇を開始する時点（負荷角速度ωＬがゼロの時点）におけるねじれ角度（θＭ−θＬ）である。したがって、折れ点Ｔ３３の横軸の値は、式（４）より、不感帯の最大値（＋ＢＬ）に一致する。
【００７４】
パラメータ同定部１２２は、ねじれ角度−機械角速度位相面図Ｐ３における折れ点Ｔ３２の横軸の値と、先に同定したモデル不感帯幅ＢＬ０及びモデル負荷側クーロン摩擦係数τｆＬｃ０とに基づいて式（４）を解くことで、モデルねじれ剛性係数ＫＲ０を同定する。
また、パラメータ同定部１２２は、ねじれ角度−機械角速度位相面図Ｐ３における折れ点Ｔ３３の横軸の値を参照してモデル不感帯幅ＢＬ０を同定する。
【００７５】
次に、図１２に示すトルク−ねじれ速度位相面図Ｐ４から導出される各折れ点Ｔ４１〜Ｔ４３について説明する。
図１２に示すトルク−ねじれ速度位相面図Ｐ４は、トルク−機械角速度位相面図Ｐ２（図１０）に示されるモータ角速度ωＭと負荷角速度ωＬとの差を示す位相面図である。したがって、トルク−ねじれ速度位相面図Ｐ４における折れ点Ｔ４１の位置は、トルク−機械角速度位相面図Ｐ２における折れ点Ｔ２１の位置に対応する。また、トルク−ねじれ速度位相面図Ｐ４における折れ点Ｔ４２の位置は、トルク−機械角速度位相面図Ｐ２における折れ点Ｔ２３の位置に対応する。更に、トルク−ねじれ速度位相面図Ｐ４における折れ点Ｔ４３の位置は、トルク−機械角速度位相面図Ｐ２における折れ点Ｔ２４の位置に対応する。
【００７６】
パラメータ同定部１２２は、トルク−ねじれ速度位相面図Ｐ４における折れ点Ｔ４２の横軸の値を参照してモデルモータ側クーロン摩擦係数τｆＭｃ０を同定する。
ここで、図１２に示す通り、トルク−ねじれ速度位相面図Ｐ４における折れ点Ｔ４２は、トルク−機械角速度位相面図Ｐ２（図１０）における折れ点Ｔ２３よりも明確に示されている。したがって、より精度よくモデルモータ側クーロン摩擦係数τｆＭｃ０を同定することができる。
【００７７】
次に、図１３に示すねじれ角度−ねじれ速度位相面図Ｐ５から導出される各折れ点Ｔ５１〜Ｔ５３について説明する。
図１３に示すねじれ角度−ねじれ速度位相面図Ｐ５は、ねじれ角度−機械角速度位相面図Ｐ３（図１１）に示されるモータ角速度ωＭと負荷角速度ωＬとの差を示す位相面図である。したがって、ねじれ角度−ねじれ速度位相面図Ｐ５における折れ点Ｔ５１の位置は、ねじれ角度−機械角速度位相面図Ｐ３における折れ点Ｔ３１の位置に対応する。また、ねじれ角度−ねじれ速度位相面図Ｐ５における折れ点Ｔ５２の位置は、ねじれ角度−機械角速度位相面図Ｐ３における折れ点Ｔ３２の位置に対応する。また、ねじれ角度−ねじれ速度位相面図Ｐ５における折れ点Ｔ５３の位置は、ねじれ角度−機械角速度位相面図Ｐ３における折れ点Ｔ３３の位置に対応する。更に、ねじれ角度−ねじれ速度位相面図Ｐ５における折れ点Ｔ５４の位置は、ねじれ角度−機械角速度位相面図Ｐ３における折れ点Ｔ３４の位置に対応する。
【００７８】
パラメータ同定部１２２は、ねじれ角度−ねじれ速度位相面図Ｐ５における折れ点Ｔ５３の横軸の値を参照して、モデル不感帯幅ＢＬ０を同定する。
また、パラメータ同定部１２２は、ねじれ角度−ねじれ速度位相面図Ｐ５における折れ点Ｔ５２の横軸の値と、先に同定したモデル不感帯幅ＢＬ０及びモデル負荷側クーロン摩擦係数τｆＬｃ０とに基づいて式（４）を解くことで、モデルねじれ剛性係数ＫＲ０を同定する。
ここで、図１３に示す通り、ねじれ角度−ねじれ速度位相面図Ｐ５における折れ点Ｔ５２及び折れ点Ｔ５３は、ねじれ角度−機械角速度位相面図Ｐ３（図１１）における折れ点Ｔ３２、Ｔ３３よりも明確に示されている。したがって、より精度よくモデル不感帯幅ＢＬ０及びモデルねじれ剛性係数ＫＲ０を同定することができる。
【００７９】
次に、図１４に示す負荷角度−トルク位相面図Ｐ６から導出される折れ点Ｔ６１及び折れ点Ｔ６３について説明する。
折れ点Ｔ６１及び折れ点Ｔ６３の時点においては、モータ側摩擦係数τｆＭ及び負荷側摩擦係数τｆＬの両方が飽和状態である。したがって、式（３）より、折れ点Ｔ６１の縦軸の値は、“−τｆＭｃ−τｆＬｃ”、折れ点Ｔ６３の縦軸の値は、“τｆＭｃ＋τｆＬｃ”と求めることができる。
【００８０】
パラメータ同定部１２２は、負荷角度−トルク位相面図Ｐ６における折れ点Ｔ６３の縦軸の値から、先に同定したモデルモータ側クーロン摩擦係数τｆＭｃ０の値を差し引いてモデル負荷側クーロン摩擦係数τｆＬｃ０を同定する。
ここで、図１４に示す通り、負荷角度−トルク位相面図Ｐ６における折れ点Ｔ６３は、ねじれ角度−トルク位相面図Ｐ１（図９）における折れ点Ｔ１４よりも明確に示されている。したがって、より精度よくモデル負荷側クーロン摩擦係数τｆＬｃ０を同定することができる。
【００８１】
本実施形態に係るパラメータ同定部１２２は、例えば、トルク−ねじれ速度位相面図Ｐ４の折れ点Ｔ４２から同定されたモデルモータ側摩擦係数τｆＭ０、トルク−機械角速度位相面図Ｐ２の折れ点Ｔ２２から同定されたモデル負荷側摩擦係数τｆＬ０を採用する。
また、パラメータ同定部１２２は、例えば、ねじれ角度−ねじれ速度位相面図Ｐ５の折れ点Ｔ５２及び折れ点Ｔ５３の各々から同定されたモデル不感帯幅ＢＬ０及びモデルねじれ剛性係数ＫＲ０を採用する。
これにより、いずれのモデルパラメータについても高い精度で同定することができる。
【００８２】
なお、パラメータ同定部１２２は、同定されたモデルモータ側クーロン摩擦係数τｆＭｃ０をモータ側摩擦特性関数ＧＭに代入する（図３（ｂ）参照）。そして、パラメータ同定部１２２は、モデルモータ側クーロン摩擦係数τｆＭｃ０に対し予め定められた係数を乗じることで、モータ側摩擦特性関数ＧＭを構成する各種定数（傾きＫ１、Ｋ２、Ｋ３、Ｋ４、オフセットτｆＭ１、τｆＭ２、τｆＭ３、τｆＭ４、τｆＭ５）を同定する。これにより、パラメータ同定部１２２は、モデルモータ側摩擦係数τｆＭ０を同定する。
同様に、パラメータ同定部１２２は、同定されたモデル負荷側クーロン摩擦係数τｆＬｃ０を付加側摩擦特性関数ＧＬに代入する。そして、パラメータ同定部１２２は、モデル負荷側クーロン摩擦係数τｆＬｃ０に対し予め定められた係数を乗じることで、負荷側摩擦特性関数ＧＬを構成する各種定数（傾き、オフセット）を同定する。これにより、パラメータ同定部１２２は、モデル負荷側摩擦係数τｆＬ０を同定する。
【００８３】
なお、本実施形態に係るパラメータ同定部１２２は、生成された各種位相面図に対し、ＰＳＯ（Particle Swarm Optimization）等を施し、各位相面図に示されたヒステリシス曲線との誤差面積が最小となるような折れ線近似演算を行うことで、各折れ点を検出する。
【００８４】
（作用、効果）
以上、第１の実施形態に係るパラメータ同定装置１２は、上述のデータ取得部１２０と、位相面図生成部１２１と、パラメータ同定部１２２と、を備える態様とする。
このようにすることで、複数の位相面図に表された複数の折れ点等から、機械系２のモータ側摩擦係数τｆＭ、負荷側摩擦係数τｆＬ、不感帯幅ＢＬ及びねじれ剛性係数ＫＲを別個独立に評価して、上記各種パラメータに対応するモデルパラメータを精度良く同定することができる。
したがって、モータ２０及び負荷２１を有する機械系２を精度よく模した２慣性系モデルＭＯＤを得ることができる。
【００８５】
また、第１の実施形態に係る位相面図生成部１２１は、位相面図として、トルクτと軸２２のねじれ速度（ωＭ−ωＬ）との関係を示すトルク−ねじれ速度位相面図Ｐ４を生成する。そして、パラメータ同定部１２２は、トルク−ねじれ速度位相面図Ｐ４における折れ点の位置に基づいてモータ２０の摩擦を示すモデルパラメータ（モデルモータ側摩擦係数τｆＭ０）を同定する。
このようにすることで、モデルモータ側摩擦係数τｆＭ０を一層精度よく同定することができる。
【００８６】
また、第１の実施形態に係る位相面図生成部１２１は、位相面図として、軸２２のねじれ角度（θＭ−θＬ）と軸２２のねじれ速度（ωＭ−ωＬ）との関係を示すねじれ角度−ねじれ速度位相面図Ｐ５を生成する。そして、パラメータ同定部１２２は、ねじれ角度−ねじれ速度位相面図Ｐ５における折れ点の位置に基づいて軸２２の剛性、及び、不感帯幅を示すモデルパラメータ（モデルねじれ剛性係数ＫＲ０、モデル不感帯幅ＢＬ０）を同定する。
このようにすることで、モデルねじれ剛性係数ＫＲ０、モデル不感帯幅ＢＬ０を一層精度よく同定することができる。
【００８７】
また、第１の実施形態に係る位相面図生成部１２１は、位相面図として、トルクτとモータ２０の角速度又は負荷２１の角速度との関係を示すトルク−機械角速度位相面図Ｐ２を生成する。そして、パラメータ同定部１２２は、トルク−機械角速度位相面図Ｐ２における折れ点の位置に基づいて負荷２１の摩擦を示すモデルパラメータ（モデル負荷側摩擦係数τｆＬ０）を同定する。
このようにすることで、モデル負荷側摩擦係数τｆＬ０を一層精度よく同定することができる。
【００８８】
また、第１の実施形態に係る位相面図生成部１２１は、位相面図として、負荷２１の角度（負荷角度θＬ）とトルクτとの関係を示す負荷角度−トルク位相面図Ｐ６を生成する。そして、パラメータ同定部１２２は、負荷角度−トルク位相面図Ｐ６における折れ点の位置に基づいて負荷２１の摩擦を示すモデルパラメータ（モデル負荷側摩擦係数τｆＬ０）を同定する。
このようにすることで、モデル負荷側摩擦係数τｆＬ０を一層精度よく同定することができる。
【００８９】
また、第１の実施形態に係る位相面図生成部１２１は、位相面図として、軸２２のねじれ角度（θＭ−θＬ）とトルクτとの関係を示すねじれ角度−トルク位相面図Ｐ１を生成する。そして、パラメータ同定部１２２は、ねじれ角度−トルク位相面図Ｐ１における折れ点の位置に基づいて、モータ２０の摩擦、負荷２１の摩擦、軸２２の剛性、及び、軸２２の不感帯幅のうちの少なくとも何れか一つを示すモデルパラメータを同定する。
このようにすることで、一つの位相面図から、モデルモータ側摩擦係数τｆＭ０、モデル負荷側摩擦係数τｆＬ０、モデルねじれ剛性係数ＫＲ０及びモデル不感帯幅ＢＬ０を同定することができる。
【００９０】
以上、第１の実施形態に係るモータ制御システム１及びパラメータ同定装置１２について詳細に説明したが、モータ制御システム１、パラメータ同定装置１２の具体的な態様は、上述のものに限定されることはなく、要旨を逸脱しない範囲内において種々の設計変更等を加えることは可能である。
【００９１】
例えば、他の実施形態に係るパラメータ同定部１２２は、複数の位相面図（Ｐ１〜Ｐ６）から同定可能なモデルパラメータのうち、同種のものについて複数の値が同定された場合、その平均値をモデルパラメータとしてもよい。例えば、パラメータ同定部１２２は、ねじれ角度−トルク位相面図Ｐ１から同定されたモデル不感帯幅ＢＬ０と、ねじれ角度−ねじれ速度位相面図Ｐ５から同定されたモデル不感帯幅ＢＬ０と、の平均値を２慣性系モデルＭＯＤに採用するモデル不感帯幅ＢＬ０としてもよい。
【００９２】
また、他の実施形態に係る位相面図生成部１２１は、一つの位相面図（ねじれ角度−トルク位相面図Ｐ１）のみを生成する態様であってもよい。この場合、パラメータ同定部１２２は、生成されたねじれ角度−トルク位相面図Ｐ１からモデルモータ側摩擦係数τｆＭ０、モデル負荷側摩擦係数τｆＬ０、モデルねじれ剛性係数ＫＲ０及びモデル不感帯幅ＢＬ０の全てを同定する。
このようにすることで、各モデルパラメータの同定処理を簡素化することができる。
【００９３】
また、他の実施形態に係る位相面図生成部１２１は、機械系２の反復動作ごとに実測データの取得を複数回繰り返し、これを平均化したものに基づいて各位相面図を生成してもよい。
このようにすることで、実測データのばらつき誤差が低減されるため、モデルパラメータの同定精度を更に高めることができる。
【００９４】
また、第１の実施形態に係るパラメータ同定装置１２は、モータ制御システム１におけるフィードフォワード機能に対してだけではなく、モータ制御システム１の異常診断や機械調整にも利用することができる。
【００９５】
また、第１の実施形態において、パラメータ同定装置１２は、上位機器から入力されるトルク指令値を「トルクτ」として、上述の各種処理を行うものとして説明したが、他の実施形態においてはこの態様に限定されない。
例えば、他の実施形態に係るパラメータ同定装置１２は、“トルク指令値”ではなく、モータ２０に流れる検出電流値から求まる“実測トルク”を「トルクτ」として用いて、上述の各種処理を行う態様であってもよい。
【００９６】
また、重力等の定常外乱がある場合には、トルクτやねじれ角度（θＭ−θＬ）に一定のオフセットがかかる。そこで、他の実施形態に係るパラメータ同定装置１２は、更に、トルクτ、ねじれ角度（θＭ−θＬ）の平均を求めてオフセット分を除去する機能、折れ点を正負で平均をとる機能を有していてもよい。
【００９７】
＜第２の実施形態＞
機械制御において、第１の実施形態で説明した位相面図（Ｐ１〜Ｐ６）によれば、機械系２の特性を示す複数の信号間の関係を可視化することができ、時間応答だけでは見つけにくい機械特性を詳細に把握することができる。このような位相面図を、基準となる位相面図と比較することで、機械系２の異常検知やヘルスモニタリングのための機械診断を行うことができる。
【００９８】
しかしながら、機械系２から取得した位相面図に基づいて異常検知、ヘルスモニタリング等を行う場合、人の目視による比較（判断）に頼る必要があり、取得した位相面図と基準位相面図との比較に基づく異常診断を自動ですることは困難である。即ち、「位相面図における２つの移動軌跡の差が大きければ異常である。」との判断は、人の目視によって行うことは容易であるが、これでは効率的に行うことができないし、類否判断が曖昧となり得る。すなわち、装置（コンピュータ等）を用いて、診断対象とする機械系２から取得した位相面図と、基準位相面図との比較に基づく異常の有無を自動的に判断する手法を確立することが望まれる。
【００９９】
ところで、地図上の移動軌跡データの類似度の検索手法として、動的時間伸縮法（ＤＴＷ：Dynamic Time Warping）がある。機械系２の機械特性を示す二次元のパラメータセット（例えば、トルクτとねじれ角度（θＭ-θＬ）との組み合わせ）の時系列データを空間上（ねじれ角度−トルク位相面図Ｐ１上）にプロットし、これを当該位相面図上における移動軌跡データとみなしてＤＴＷを適用することで、２つの移動軌跡データの類似度に基づく異常診断を行うことができる。
【０１００】
しかしながら、ＤＴＷはノイズに弱いため、位相面図に用いた場合、信号のノイズにより精度が悪くなりやすい。
また、ガタや非線形摩擦の影響を把握可能な微小応答を位相面図に表すと、位相面図上においてサンプルポイントが著しく偏在し、そのままではＤＴＷを適用できない。
更に、移動軌跡データのデータ数が大きい場合、ＤＴＷは計算量が過大となる。
【０１０１】
第２の実施形態に係る移動軌跡成形装置、及び、当該移動軌跡成形装置を備える異常診断システムは、上記課題に鑑みてなされたものである。以下、第２の実施形態に係る移動軌跡成形装置、及び、当該移動軌跡成形装置を備える異常診断システムについて、図１５〜図２２を参照しながら詳細に説明する。
【０１０２】
（全体構成）
図１５は、第２の実施形態に係るモータ制御システム、機械系、及び、異常診断システムの全体構成を示す図である。
【０１０３】
モータ制御システム１及び機械系２については、第１の実施形態と同様であるため、詳細な説明を省略する。
異常診断システム３は、機械系２の異常診断、ヘルスモニタリングを行うシステムである。異常診断システム３は、一般に用いられる汎用のコンピュータ等であってよい。
異常診断システム３は、モータ制御システム１の位相面図生成部１２１（図６参照）を通じて、移動軌跡データを取得する。ここで、「移動軌跡データ」とは、多次元のパラメータセット（多次元ベクトル）の時系列データであって、本実施形態においては、ねじれ角度−トルク位相面図Ｐ１上のプロットで示される、トルクτとねじれ角度（θＭ−θＬ）との組み合わせの時系列データ（実測データ）である（図８、図９等を参照）。
なお、「移動軌跡データ」は、「トルクτとねじれ角度（θＭ−θＬ）との組み合わせの時系列データ」に限定されない。他の実施形態においては、移動軌跡データは、他の位相面図（Ｐ２〜Ｐ６）上に示されるパラメータセット（例えば、トルクτとねじれ速度（ωＭ−ωＬ）との組み合わせ、負荷角度θＬとトルクτとの組み合わせ等）の時系列データであってもよい。
【０１０４】
第２の実施形態に係る異常診断システム３は、内部に移動軌跡成形装置３０を備えている。移動軌跡成形装置３０は、モータ制御システム１を通じて得られた移動軌跡データ（移動軌跡）を適切に成形する機能を有している。
【０１０５】
（機能構成）
図１６は、第２の実施形態に係る異常診断システムの機能構成を示す図である。
図１６に示すように、異常診断システム３は、ＣＰＵ３０Ａと、メモリ３１と、表示部３２と、操作部３３と、接続インターフェイス３４と、記録媒体３５と、を備えている。
ＣＰＵ３０Ａは、異常診断システム３の動作全体を司るプロセッサである。ＣＰＵ３０Ａは、記録媒体等に格納されたプログラムやデータをメモリ３１上に読み出し、当該プログラムに規定される処理を実行することで、後述の各機能を実現する。
メモリ３１は、ＣＰＵ３０Ａのワークエリア等として用いられる揮発性のメモリ（ＲＡＭ）である。
表示部３２は、例えば、液晶ディスプレイ、有機ＥＬディスプレイ等により実現され、ＣＰＵ３０Ａによる処理結果を表示する。
操作部３３は、例えば、マウス、タッチパネル及びキーボード等で構成され、操作者（ユーザ）の指示を受けてＣＰＵ３０Ａに各種操作等を入力する。
接続インターフェイス３４は、外部装置とのインターフェイスである。特に、本実施形態においては、接続インターフェイス３４は、モータ制御システム１と接続される。
記録媒体３５は、例えば、ＨＤＤ（Hard Disk Drive）、ＳＳＤ（Solid State Drive）等の大容量記録デバイスにより実現され、ＯＳ（Operation System）、アプリケーションプログラム、及び、各種データ等を記憶する。本実施形態に係る記録媒体３５は、例えば、後述する「基準移動軌跡データ」などを予め記録保存している。
【０１０６】
第２の実施形態に係るＣＰＵ３０Ａは、上述のプログラムに基づいて動作することで、図２に示す移動軌跡データ取得部３００、平滑化処理部３０１、正規化処理部３０２、等間隔化処理部３０３、及び、異常判定部３０４としての機能を発揮する。
【０１０７】
移動軌跡データ取得部３００は、接続インターフェイス３４を通じて、モータ制御システム１から移動軌跡データを取得する。上述した通り、本実施形態の例としては、移動軌跡データ取得部３００は、トルクτとねじれ角度（θＭ−θＬ）との組み合わせである、二次元のパラメータセット（二次元ベクトル）の時系列データを取得する。
【０１０８】
平滑化処理部３０１は、予め規定された平滑化フィルタ（後述する「バイラテラル軌跡フィルタｆ」）を適用して、移動軌跡データ取得部３００が取得した移動軌跡データを平滑化する。
【０１０９】
正規化処理部３０２は、上記二次元のパラメータセット（トルクτとねじれ角度（θＭ−θＬ）との組み合わせ）を、次元別（トルクτ、ねじれ角度（θＭ−θＬ）別）に正規化する。
【０１１０】
等間隔化処理部３０３は、平滑化処理部３０１によって平滑化された移動軌跡データによって移動軌跡曲線を特定する。そして、等間隔化処理部３０３は、当該移動軌跡曲線を二次元空間内（位相面図内）において等間隔で刻むとともに、当該等間隔で刻まれた移動軌跡曲線上の各点の座標のデータ系列を取得する。
【０１１１】
上述の移動軌跡データ取得部３００、平滑化処理部３０１、正規化処理部３０２及び異常判定部３０４は、移動軌跡成形装置３０を構成する。
【０１１２】
異常判定部３０４は、等間隔化処理部３０３が取得したデータ系列からなる評価対象移動軌跡データ（移動軌跡成形装置３０によって成形された移動軌跡データ）と、予め用意された基準移動軌跡データとの対比に基づいて、移動軌跡データに異常があるか否かを判定する。
【０１１３】
（異常診断システムの処理フロー）
図１７は、第２の実施形態に係る異常診断システムの処理フローを示す図である。
図１８〜図２２は、それぞれ、第２の実施形態に係る異常診断システムの各処理を詳細に説明するための第１の図〜第５の図である。
図１７に示す一連の処理フローは、例えば、機械系２に対して定期的に行われるヘルスチェック（異常診断）を行う際に実施される。
【０１１４】
まず、ＣＰＵ３０Ａの移動軌跡データ取得部３００は、接続インターフェイス３４を通じて、モータ制御システム１から移動軌跡データを取得する（ステップＳ１１）。移動軌跡データ取得部３００が取得する移動軌跡データの特徴について、図１８、図１９を参照しながら説明する。
【０１１５】
図１８は、トルクτ（モータトルク）の時系列データと、ねじれ角度（θＭ−θＬ）の時系列データの例を示している。また、図１９は、トルクτの時系列データと、ねじれ角度（θＭ−θＬ）の時系列データと、に基づいて作成される移動軌跡データＥ１（ねじれ角度−トルク位相面図Ｐ１上にプロットされるサンプルポイントの群）を示している。
【０１１６】
ここで、第１の実施形態では、ガタや非線形摩擦の特性を位相面図上に顕在化させる目的で、機械系２の反復動作を小振幅かつ低周波数（微小応答）とし、モータ２０の慣性トルクＪＭｓ２θＭ、負荷２１の慣性トルクＪＬｓ２θＬ、モータ２０の粘性トルクＤＭｓθＭ、及び、負荷２１の粘性トルクＤＬｓθＬの影響を小さくすることを説明した。
この場合、図１８に示すように、各時系列データは、負荷角度θＬの方向反転時付近の時間領域（領域Ｒ１、Ｒ２）に属するサンプルポイントの数よりも、負荷角度θＬが同一方向への移動を続ける時間領域（トルクτ、ねじれ角度（θＭ−θＬ）が安定する時間領域）（領域Ｑ１、Ｑ２）に属するサンプルポイントの数の方が大きくなる。
また、図１９に示すように、ねじれ角度−トルク位相面図Ｐ１上に描かれる移動軌跡の大部分は、領域Ｒ１、Ｒ２に属するサンプルポイントによって形成され、領域Ｑ１、Ｑ２に属する多数のサンプルポイントは、移動軌跡の形成にほとんど寄与しない。即ち、微小応答時に取得された移動軌跡データＥ１の各サンプルポイントは、ねじれ角度−トルク位相面図Ｐ１上に描かれる移動軌跡の一部の領域に著しく偏在する。
【０１１７】
図１７に示す処理フローにおいて、次に、ＣＰＵ３０Ａの正規化処理部３０２は、移動軌跡データＥ１の正規化処理を行う（ステップＳ１２）。ここで、移動軌跡データＥ１の各次元（トルクτ、ねじれ角度（θＭ−θＬ））は、スケールや単位が異なるため、それぞれが属する数値範囲が大きく異なる。そこで、正規化処理部３０２は、トルクτ及びねじれ角度（θＭ−θＬ）を同一のスケール、単位として扱えるように、“−１〜＋１”の範囲に収まるように正規化を行う。
【０１１８】
次に、ＣＰＵ３０Ａの平滑化処理部３０１は、ステップＳ１２で正規化された移動軌跡データＥ１に対し、「エッジ保存の平滑化処理」を行う（ステップＳ１３）。「エッジ保存の平滑化処理」とは、具体的には、移動軌跡データＥ１のうち大きく変化するエッジ領域（即ち、図１８、図１９に示す領域Ｒ１、Ｒ２）のサンプルポイントを動かすことなく保存して、それ以外の領域（図１８、図１９に示す領域Ｑ１、Ｑ２）のサンプルポイントのみを選択的に平滑化する処理を指す。
この「エッジ保存の平滑化処理」の詳細については、図２０を参照しながら説明する。
【０１１９】
図２０は、エッジ保存の平滑化処理を実現する平滑化フィルタであるバイラテラル軌跡フィルタｆを示している。
バイラテラル軌跡フィルタｆは、２つの重み付け関数（第１の重み付け関数Ｗ１、第２の重み付け関数Ｗ２）の積によって規定される。
第１の重み付け関数Ｗ１は、サンプルポイントどうしの時間の差が大きくなるにつれて平滑化の影響度が小さくなる重み付けを規定する。具体的には、第１の重み付け関数Ｗ１は、横軸を“サンプリングされた時間の差”とする、平均μ＝０、標準偏差σｇのガウス分布関数である。
第２の重み付け関数Ｗ２は、サンプルポイントどうしの位相面（ねじれ角度−トルク位相面図Ｐ１）上の位置が離れるにつれて平滑化の影響度が小さくなる重み付けを規定する。具体的には、第２の重み付け関数Ｗ２は、横軸を“位相面上の位置の差”とする、平均μ＝０、標準偏差σｃのガウス分布関数である。
【０１２０】
バイラテラル軌跡フィルタｆは、例えば、式（５）のように規定される。
【０１２１】
【数５】
2019021036.tif 000007
【０１２２】
ここで、式（５）の“Ｙ（ｉ）”は、移動軌跡データのｉ番目のサンプルポイントの二次元ベクトル（トルクτとねじれ角度（θＭ−θＬ）との組み合わせ）である。
【０１２３】
バイラテラル軌跡フィルタｆによれば、以下のようなフィルタリング効果が得られる。
まず、第１の重み付け関数Ｗ１により、主に、サンプリングされた時間の差が小さいサンプルポイントどうしで平滑化（平均化）がなされる。これにより、短い時間間隔で細かく変動するノイズ成分を除去（平滑化）することができる。
他方、第２の重み付け関数Ｗ２によれば、隣り合うサンプルポイント間の距離（位相面上の位置の差）が大きいほど、それらのサンプルポイントの平滑化への寄与度（重み）は小さくなる。したがって、時間間隔が小さいサンプルポイントどうしであっても、当該時間間隔で大きな移動軌跡の変化があった場合には、それらのサンプルポイントどうしでの平滑化はなされなくなる。これにより、エッジ保存の効果が得られる。
以上より、移動軌跡データＥ１にバイラテラル軌跡フィルタｆを適用することで、領域Ｑ１、Ｑ２に属するサンプルポイントのみが選択的に平滑化され、小さい時間間隔で大きな変化が生じたエッジ部分（領域Ｒ１、Ｒ２）に属するサンプルポイントは平滑化されず、そのまま保存されることとなる。
以下、バイラテラル軌跡フィルタｆにより「エッジ保存の平滑化」がなされた移動軌跡データを、「平滑化移動軌跡データＥ２」と表記して説明する。
【０１２４】
なお、第１の重み付け関数Ｗ１の標準偏差σｇと第２の重み付け関数Ｗ２の標準偏差σｃとは、平滑化処理の対象とする移動軌跡データの、平滑化したいノイズの大きさ、及び、保存したい（平滑化したくない）エッジの大きさ等に応じて適宜調整される。
【０１２５】
図１７に示す処理フローにおいて、次に、ＣＰＵ３０Ａの等間隔化処理部３０３は、平滑化移動軌跡データＥ２に対し、「等間隔化処理」を行う（ステップＳ１４）。「等間隔化処理」とは、平滑化移動軌跡データＥ２に示される移動軌跡の形を維持したまま、当該移動軌跡の一部分に偏在するサンプルポイントを均等化させる処理である。
この「等間隔化処理」の詳細については、図２１を参照しながら説明する。
【０１２６】
等間隔化処理部３０３は、図２１に示す３つの処理ステップＳ１４１、Ｓ１４２、Ｓ１４３を実行する。
【０１２７】
まず、等間隔化処理部３０３は、平滑化移動軌跡データＥ２において時系列上で隣り合うサンプルポイント間の距離を計測するとともに、各距離の総和を演算することで移動軌跡の全長を特定する（ステップＳ１４１）。
具体的には、図２１のステップＳ１４１に示すように、等間隔化処理部３０３は、平滑化移動軌跡データＥ２を構成する各サンプルポイントｅ１、ｅ２、・・・、ｅＭ（“Ｍ”は、平滑化移動軌跡データＥ２の総サンプルポイント数）のうち、隣り合うものどうしの各距離ｄ１、ｄ２、・・ｄ（Ｍ−１）を演算する。ここで、ｉ番目のサンプルポイントｅｉは、ねじれ角度−トルク位相面図Ｐ１上の座標（ｘｉ、ｙｉ）にプロットされるものとする。この場合、ｉ−１番目のサンプルポイントｅ（ｉ−１）からｉ番目のサンプルポイントｅｉまでの距離ｄｉ、及び、１番目のサンプルポイントｅ１からｉ番目のサンプルポイントｅｉまでの距離の総和ｓｉは、式（６）によって求められる。
【０１２８】
【数６】
2019021036.tif 000008
【０１２９】
等間隔化処理部３０３は、式（６）に基づいて平滑化移動軌跡データＥ２の最後のサンプルポイントｅＭまでの距離の総和ｓＭ（移動軌跡の全長）を算出する。
【０１３０】
次に、等間隔化処理部３０３は、平滑化移動軌跡データＥ２を構成する各サンプルポイントｅ１、ｅ２、・・、ｅＭに基づく曲線（以下、移動軌跡曲線Ｈと表記する。）を特定する。ある座標空間上にプロットされる複数のポイントに基づいて曲線を特定する手法は、既存技術を適用できる。例えば、等間隔化処理部３０３は、各サンプルポイントｅ１、ｅ２、・・、ｅＭを滑らかに接続するスプライン曲線によって移動軌跡曲線Ｈを特定してもよいし、曲線に振動を生じさせないようにする区分的三次エルミート内挿多項式で移動軌跡曲線Ｈを特定してもよい。
【０１３１】
次に、等間隔化処理部３０３は、ステップＳ１４２で求めた移動軌跡曲線Ｈを、ねじれ角度−トルク位相面図Ｐ１上において等間隔（距離ｄｃ）で刻み、刻まれた各点を新たなサンプルポイントｅ１’、ｅ２’、・・として特定する（ステップＳ１４３）。ここで、移動軌跡曲線Ｈを均等に１００分割する場合、等間隔化処理部３０３は、距離の総和ｓＭを分割数（１００）で割ることで距離ｄｃを算出する（ｄｃ＝ｓＭ／１００）。
以下、移動軌跡曲線Ｈが等間隔で刻まれて規定される新たなサンプルポイントｅ１’、ｅ２’、・・のデータ系列を、「評価対象移動軌跡データＥ２’」と表記して説明する。
【０１３２】
図１７において、次に、ＣＰＵ３０Ａの異常判定部３０４は、ステップＳ１４で取得された評価対象移動軌跡データＥ２’と、予め用意された基準移動軌跡データとのＤＴＷ距離を演算する。ここで、「基準移動軌跡データ」とは、機械系２が正常な状態に取得される、代表的な移動軌跡データである。即ち、評価対象移動軌跡データＥ２’が基準移動軌跡データと類似していれば、診断の対象である機械系２は「異常なし」と判断することができる。また、「ＤＴＷ距離」とは、動的時間伸縮法に基づいて２つの移動軌跡データの類似度を示す指標となるパラメータである。
以下、図２２を参照しながら、ＤＴＷ距離の求め方について説明する。
【０１３３】
図２２には、ねじれ角度−トルク位相面図Ｐ１上にプロットされた２つのデータ系列Ｘ、Ｙを示している。図２２に示すデータ系列Ｘが評価対象移動軌跡データＥ２’であり、データ系列Ｙが基準移動軌跡データであるものとする。
異常判定部３０４は、データ系列Ｘのサンプルポイントごとに、データ系列Ｙの各サンプルポイントのうち最も距離が短いサンプルポイントを特定する。例えば、データ系列Ｘのうちのｉ番目のサンプルポイントｘ（ｉ）と最も距離が短いデータ系列Ｙのサンプルポイントは、ｊ番目のサンプルポイントｙ（ｊ）となる。また、データ系列Ｘのうちの（ｉ＋１）番目のサンプルポイントｘ（ｉ＋１）と最も距離が短いデータ系列Ｙのサンプルポイントも、ｊ番目のサンプルポイントｙ（ｊ）となる。
【０１３４】
異常判定部３０４は、特定されたデータ系列Ｘのサンプルポイントとデータ系列Ｙのサンプルポイントとの組（例えば｛ｘ（ｉ），ｙ（ｊ）｝の組）の距離（ｃｏｓｔ）を算出し、更に、全ての組のｃｏｓｔの総和（ＤＴＷ距離）を算出する。
なお、ＤＴＷ距離の具体的な算出式は、式（７）に示す通りである。
【０１３５】
【数７】
2019021036.tif 000009
【０１３６】
式（７）において、ＤＴＷ（ｉ，ｊ）は、｛ｘ（ｉ），ｙ（ｊ）｝の組までのＤＴＷ距離である。
【０１３７】
図１７において、最後に、異常判定部３０４は、ステップＳ１５で算出したＤＴＷ距離が予め規定された所定の判定閾値を上回ったか否かの判定を行う（ステップＳ１６）。ステップＳ１５で算出したＤＴＷ距離が所定の判定閾値を上回っている場合、評価対象移動軌跡データＥ２’が基準移動軌跡データと大きく異なっていることから、異常判定部３０４は、「異常あり」との判定をする。他方、ステップＳ１５で算出したＤＴＷ距離が所定の判定閾値以下であった場合、評価対象移動軌跡データＥ２’が基準移動軌跡データと類似していることから、異常判定部３０４は、「異常なし」との判定をする。
【０１３８】
（作用・効果）
以上の通り、第２の実施形態に係る異常診断システム３（移動軌跡成形装置３０）は、多次元のパラメータセット（トルクτとねじれ角度（θＭ−θＬ）との組み合わせ）の時系列データである移動軌跡データを取得する。また、異常診断システム３は、予め規定された平滑化フィルタ（バイラテラル軌跡フィルタｆ）を適用して、取得した移動軌跡データを平滑化する。そして、上述のバイラテラル軌跡フィルタｆは、パラメータセットどうしの時間の差が大きくなるにつれて平滑化の影響度が小さくなる重み付けを規定する第１の重み付け関数Ｗ１と、パラメータセットどうしの多次元空間（位相面図）内の位置が離れるにつれて平滑化の影響度が小さくなる重み付けを規定する第２の重み付け関数Ｗ２と、の積からなる。
このようにすることで、移動軌跡の特徴が強く現れる“エッジ”（短時間内に位置が大きく変化する領域）を平滑化させることなく保存し、ノイズ成分のみを平滑化（除去）することができる。即ち、異常診断システム３によれば、移動軌跡データのノイズを適切に除去することができる。
なお、ＤＴＷはノイズの影響を受けて類否判断を誤りやすい特徴を有している。しかし、上述のとおり、エッジ保存の平滑化により移動軌跡データのノイズを適切に除去することで、ＤＴＷに基づく移動軌跡の類否判断の信頼性を高めることができる。
【０１３９】
また、第２の実施形態に係る異常診断システム３（移動軌跡成形装置３０）は、多次元のパラメータセット（トルクτとねじれ角度（θＭ−θＬ）との組み合わせ）を、次元別に正規化する。
このようにすることで、異なる単位、異なるスケールのパラメータの組み合わせを同系列のパラメータとして扱うことができるので、ＤＴＷに基づく移動軌跡の類否判断（異常診断）をより適切に行うことができる。
【０１４０】
更に、第２の実施形態に係る異常診断システム３（移動軌跡成形装置３０）は、平滑化移動軌跡データＥ２によって特定される移動軌跡曲線Ｈを、多次元空間（位相面図）内において等間隔で刻むとともに、当該等間隔で刻まれた移動軌跡曲線Ｈ上の各点の座標のデータ系列（評価対象移動軌跡データＥ２’）を取得する。
このようにすることで、移動軌跡上におけるサンプルポイントの偏在を解消することができる。これにより、ＤＴＷに基づく移動軌跡の類否判断（異常診断）をより適切に行うことができる。
また、移動軌跡データのサンプルポイント数を適切に（例えば１００ポイント程度）に削減することができるので、ＤＴＷの計算量を大幅に削減することができる。
【０１４１】
更に、第２の実施形態に係る異常診断システム３は、評価対象移動軌跡データＥ２’と、予め用意された基準移動軌跡データとの対比に基づいて、移動軌跡データに異常があるか否かを判定する。
また、異常診断システム３は、評価対象移動軌跡データＥ２’と基準移動軌跡データとの間のＤＴＷ距離を演算する。そして、異常診断システム３は、算出されたＤＴＷ距離が所定の判定閾値を上回ったか否かに基づいて、移動軌跡データに異常があるか否かを判定する。
これにより、位相面図に描かれる移動軌跡の類否判断を、ＤＴＷ距離に基づいて精度良く行うことができる。したがって、位相面図を用いた機械系２の異常診断を実現することができる。
【０１４２】
なお、同一空間（位相面図）内にプロットされる２つのデータ系列（移動軌跡）の類否判断の手法は、第２の実施形態で説明したＤＴＷに限定されない。他の実施形態に係る異常診断システム３は、公知の手法であってＤＴＷ以外の類否判断の手法を用いてもよい。
【０１４３】
また、第２の実施形態においては、移動軌跡データは、二次元のパラメータセット（トルクτとねじれ角度（θＭ−θＬ）との組み合わせ）の時系列であるものとして説明したが、他の実施形態においてはこの態様に限定されない。
他の実施形態においては、移動軌跡データは、三次元のパラメータセット（例えば、トルクτ、ねじれ角度（θＭ−θＬ）及びねじれ速度（ωＭ−ωＬ）の組み合わせ等）であってもよい。この場合、異常診断システム３による各処理（正規化処理、エッジ保存の平滑化処理、等間隔化処理、及び、ＤＴＷ距離の演算処理）は、三次元空間内にプロットされる移動軌跡データについて実行される。
更に、移動軌跡データは、四次元以上のパラメータセットの時系列データであってもよい。
【０１４４】
また、第２の実施形態に係る移動軌跡成形装置３０は、異常診断システム３に組み込まれ、機械系２の異常診断、ヘルスモニタリングの用途に用いるものとして説明したが、他の実施形態に係る移動軌跡成形装置３０は、他の用途に用いられてもよい。
【０１４５】
例えば、他の実施形態に係る移動軌跡成形装置３０は、小型飛行体（ドローン等）が実空間（三次元空間）内を移動する際に取得される三次元移動軌跡データに対し、ノイズ除去、等間隔化処理等を実行し、移動軌跡を成形する用途に用いられてもよい。
また、他の実施形態に係る移動軌跡成形装置３０は、アウトレットモール、百貨店、コンビニエンスストア等を徒歩移動する利用客の移動軌跡データに対し、ノイズ除去、等間隔化処理等を実行し、移動軌跡を成形する用途に用いられてもよい。
以上のような用途に用いることで、移動軌跡の類否判断、統計的分析を高確度、かつ、低負荷で実施することができる。
【０１４６】
また、上述の各実施形態においては、異常診断システム３及び移動軌跡成形装置３０の機能を実現するためのプログラムをコンピュータ読み取り可能な記録媒体に記録して、この記録媒体に記録されたプログラムをコンピュータシステムに読み込ませ、実行することにより各手順を行うものとしている。ここで、上述した異常診断システム３及び移動軌跡成形装置３０の各処理の過程は、プログラムの形式でコンピュータ読み取り可能な記録媒体に記憶されており、このプログラムをコンピュータが読み出して実行することによって上記各種処理が行われる。ここで、コンピュータ読み取り可能な記録媒体とは、磁気ディスク、光磁気ディスク、ＣＤ−ＲＯＭ、ＤＶＤ−ＲＯＭ、半導体メモリ等をいう。また、このコンピュータプログラムを通信回線によってコンピュータに配信し、この配信を受けたコンピュータが当該プログラムを実行するようにしても良い。
また、異常診断システム３及び移動軌跡成形装置３０の各機能が、ネットワークで接続される複数の装置に渡って具備される態様であってもよい。
【０１４７】
以上、本発明のいくつかの実施形態を説明したが、これらの実施形態は、例として示したものであり、発明の範囲を限定することは意図していない。これら実施形態は、その他の様々な形態で実施されることが可能であり、発明の要旨を逸脱しない範囲で種々の省略、置き換え、変更を行うことができる。これら実施形態やその変形は、発明の範囲や要旨に含まれると同様に、特許請求の範囲に記載された発明とその均等の範囲に含まれるものとする。
【符号の説明】
【０１４８】
１  モータ制御システム
１０  フィードバック制御部
１１  フィードフォワード制御部
１２  パラメータ同定装置
１２０  データ取得部
１２１  位相面図生成部
１２２  パラメータ同定部
２  機械系
２０  モータ
２１  負荷
２２  軸
３  異常診断システム
３０  移動軌跡成形装置
３０Ａ  ＣＰＵ
３００  移動軌跡データ取得部
３０１  平滑化処理部
３０２  正規化処理部
３０３  等間隔化処理部
３０４  異常判定部
３１  メモリ
３２  表示部
３３  操作部
３４  接続インターフェイス
３５  記録媒体
θｔ  目標角度（目標回転角度）
τ、τ’  トルク
θＭ  モータ角度（モータの角度）
ωＭ  モータ角速度（モータの角速度）
ＪＭ  モータ側慣性モーメント
ＤＭ  モータ側粘性係数
τｆＭ  モータ側摩擦係数
τｆＭｃ  モータ側クーロン摩擦係数
θＬ  負荷角度（負荷の角度）
ωＬ  負荷角速度（負荷の角速度）
ＪＬ  負荷側慣性モーメント
ＤＬ  負荷側粘性係数
τｆＬ  負荷側摩擦係数
τｆＬｃ  負荷側クーロン摩擦係数
ＫＲ  ねじれ剛性係数
ＤＲ  ねじれ粘性係数
ＢＬ  不感帯幅
ＭＯＤ  ２慣性系モデル
ＪＭ０  モデルモータ側慣性モーメント
ＤＭ０  モデルモータ側粘性係数
τｆＭ０  モデルモータ側摩擦係数
τｆＭｃ０  モデルモータ側クーロン摩擦係数
ＪＬ０  モデル負荷側慣性モーメント
ＤＬ０  モデル負荷側粘性係数
τｆＬ０  モデル負荷側摩擦係数
τｆＬｃ０  モデル負荷側クーロン摩擦係数
ＢＬ０  モデル不感帯幅
ＫＲ０  モデルねじれ剛性係数
Ｆ１  不感帯特性関数
Ｆ１’  ガタ変位関数
ＧＭ  モータ側摩擦特性関数
ＧＬ  負荷側摩擦特性関数
Ｐ１  ねじれ角度−トルク位相面図
Ｐ２  トルク−機械角速度位相面図
Ｐ３  ねじれ角度−機械角速度位相面図
Ｐ４  トルク−ねじれ速度位相面図
Ｐ５  ねじれ角度−ねじれ速度位相面図
Ｐ６  負荷角度−トルク位相面図
ｆ  バイラテラル軌跡フィルタ
Ｗ１  第１の重み付け関数
Ｗ２  第２の重み付け関数
Ｈ  移動軌跡曲線
【図１】
2019021036.tif 000010
【図２】
2019021036.tif 000011
【図３】
2019021036.tif 000012
【図４】
2019021036.tif 000013
【図５】
2019021036.tif 000014
【図６】
2019021036.tif 000015
【図７】
2019021036.tif 000016
【図８】
2019021036.tif 000017
【図９】
2019021036.tif 000018
【図１０】
2019021036.tif 000019
【図１１】
2019021036.tif 000020
【図１２】
2019021036.tif 000021
【図１３】
2019021036.tif 000022
【図１４】
2019021036.tif 000023
【図１５】
2019021036.tif 000024
【図１６】
2019021036.tif 000025
【図１７】
2019021036.tif 000026
【図１８】
2019021036.tif 000027
【図１９】
2019021036.tif 000028
【図２０】
2019021036.tif 000029
【図２１】
2019021036.tif 000030
【図２２】
2019021036.tif 000031