(57)【要約】【課題】  パーティクルフィルタを用いて制御対象となるシステムのパラメータを同定するに際して、パラメータの同定精度を向上できる技術を開示する。【解決手段】  システムパラメータ同定装置１０では、パーティクルフィルタ実行部２４が、複数のパーティクルのそれぞれに、システムの状態ベクトルを割り当て、観測値と、各パーティクルの推定値とに基づいて、各パーティクルの重みを算出し、その重みに基づいてパーティクルを更新することで状態ベクトルを推定してパラメータを同定する。パーティクルフィルタ実行部２４では、パーティクルは、任意の期間内の最初のステップからｎステップ目毎に（ｎ：２以上の整数）更新され、かつ、各パーティクルの重みは、パーティクルが更新される１個の更新周期内に含まれるｎ個のステップのうちの少なくとも２個のステップにおいて取得された観測値と推定値とに基づいて算出される。【選択図】  図１
2017083922.tif 000002
【特許請求の範囲】
【請求項１】
パーティクルフィルタを用いて制御対象のシステムのパラメータを同定するシステムパラメータ同定装置であって、
入力値取得部と、観測値取得部と、第１関数記憶部と、第２関数記憶部と、状態ベクトル算出部と、推定値算出部と、パーティクルフィルタ実行部と、を備えており、
前記入力値取得部は、前記システムへの入力値を、任意の期間においてステップ毎に取得可能であり、
前記観測値取得部は、前記システムからの出力の観測値を、任意の期間においてステップ毎に取得可能であり、
前記第１関数記憶部は、前記任意の期間内の任意のステップにおける前記システムの状態ベクトルを、前記任意のステップの１つ前のステップにおける前記システムの状態ベクトルに基づいて算出する第１関数を記憶しており、前記システムの状態ベクトルは、前記システムの状態量と、前記システムのパラメータと、前記入力値に関連するシステムノイズとを要素に含んでおり、
前記第２関数記憶部は、前記任意のステップにおける前記システムの状態ベクトルに基づいて、前記任意のステップにおける前記システムの出力の推定値を算出する第２関数を記憶しており、
前記状態ベクトル算出部は、前記任意のステップの１つ前のステップにおける前記システムの状態ベクトルを前記第１関数記憶部に記憶されている前記第１関数に入力して、前記任意のステップにおける前記システムの状態ベクトルを算出し、
前記推定値算出部は、前記任意のステップにおける前記システムの状態ベクトルを前記第２関数記憶部に記憶されている前記第２関数に入力して、前記任意のステップにおける前記システムの出力の推定値を算出し、
前記パーティクルフィルタ実行部は、複数のパーティクルのそれぞれに、前記システムの状態ベクトルを割り当て、前記観測値取得部から取得した観測値と、パーティクル毎に前記推定値算出部から取得した推定値とに基づいて、各パーティクルの重みを算出し、その重みに基づいてパーティクルを更新することで前記システムの状態ベクトルを推定して前記システムのパラメータを同定し、
前記パーティクルフィルタ実行部では、
パーティクルは、前記任意の期間内の最初のステップからｎステップ目毎に（ｎ：２以上の整数）更新され、
各パーティクルの重みは、パーティクルが更新される１個の更新周期内に含まれるｎ個のステップのうちの少なくとも２個のステップにおいて取得された観測値と推定値とに基づいて算出される、システムパラメータ同定装置。
【請求項２】
前記システムのパラメータにパラメータノイズを付与するノイズ付与部をさらに備えており、
前記ノイズ付与部は、パーティクルが更新されるステップにおいて、前記システムのパラメータに前記パラメータノイズを付与する、請求項１に記載のシステムパラメータ同定装置。
【請求項３】
前記パラメータノイズの分散が可変である、請求項２に記載のシステムパラメータ同定装置。
【請求項４】
前記任意の期間内においてパーティクルがｍ回（ｍ：２以上の整数）更新されるとすると、ｊ回目（ｊ：１〜ｍ−１までの任意の整数）に更新されるときのステップで付与される前記パラメータノイズの前記分散は、ｎ・ｊ＋１ステップ目からｎ・ｊ＋ｎステップ目までのｎ個のステップのうちの少なくとも１個のステップにおいて取得された入力値又は観測値に基づいて決定される、請求項３に記載のシステムパラメータ同定装置。
【請求項５】
各パーティクルの重みは、前記少なくとも２個のステップにおける同時確率から算出される、請求項１〜４のいずれか一項に記載のシステムパラメータ同定装置。
【請求項６】
請求項１〜５のいずれか一項に記載のシステムパラメータ同定装置を備える、前記システムを制御する制御器であって、
前記制御器の制御パラメータは、前記パーティクルフィルタ実行部により更新されたパーティクルの分散及び平均に基づいて決定される、制御器。
【請求項７】
パーティクルフィルタを用いて制御対象のシステムのパラメータを同定するシステムパラメータ同定方法であって、
入力値取得工程と、観測値取得工程と、状態ベクトル算出工程と、推定値算出工程と、パーティクルフィルタ実行工程と、を備えており、
前記入力値取得工程では、前記システムへの入力値を、任意の期間においてステップ毎に取得し、
前記観測値取得工程では、前記システムの出力の観測値を、任意の期間においてステップ毎に取得し、
前記状態ベクトル算出工程では、前記システムの状態量と、前記システムのパラメータと、入力値に関するシステムノイズとを要素に含む、前記システムの状態ベクトルであって、前記任意の期間内の任意のステップにおける前記システムの状態ベクトルを、前記任意のステップの１つ前のステップにおける前記システムの状態ベクトルに基づいて算出する第１関数に、前記１つ前のステップにおける前記システムの状態ベクトルを入力して前記任意のステップにおける前記システムの状態ベクトルを算出し、
前記推定値算出工程では、前記任意のステップにおける前記システムの状態ベクトルに基づいて前記任意のステップにおける前記システムの出力の推定値を算出する第２関数に、前記任意のステップにおける前記システムの状態ベクトルを入力して前記任意のステップにおける前記システムの出力の推定値を算出し、
前記パーティクルフィルタ実行工程では、複数のパーティクルのそれぞれに、前記システムの状態ベクトルを割り当て、前記観測値取得工程で取得した観測値と、パーティクル毎に前記推定値算出工程で取得した推定値とに基づいて、各パーティクルの重みを算出し、その重みに基づいてパーティクルを更新することで前記システムの状態ベクトルを推定して前記システムのパラメータを同定し、
前記パーティクルフィルタ実行工程では、
パーティクルは、前記任意の期間内の最初のステップからｎステップ目毎に（ｎ：２以上の整数）更新され、
各パーティクルの重みは、パーティクルが更新される１個の更新周期内に含まれるｎ個のステップのうちの少なくとも２個のステップにおいて取得された観測値と推定値とに基づいて算出される、システムパラメータ同定方法。
【請求項８】
パーティクルフィルタを用いて制御対象のシステムのパラメータを同定するためのコンピュータプログラムであって、
コンピュータに、入力値取得処理と、観測値取得処理と、状態ベクトル算出処理と、推定値算出処理と、パーティクルフィルタ実行処理と、を実行させ、
前記入力値取得処理では、前記システムへの入力値を、任意の期間においてステップ毎に取得し、
前記観測値取得処理では、前記システムの出力の観測値を、任意の期間においてステップ毎に取得し、
前記状態ベクトル算出処理では、前記システムの状態量と、前記システムのパラメータと、入力値に関連するシステムノイズとを要素に含む、前記システムの状態ベクトルであって、前記任意の期間内の任意のステップにおける前記システムの状態ベクトルを、前記任意のステップの１つ前のステップにおける前記システムの状態ベクトルに基づいて算出する第１関数に、前記１つ前のステップにおける前記システムの状態ベクトルを入力して前記任意のステップにおける前記システムの状態ベクトルを算出し、
前記推定値算出処理では、前記任意のステップにおける前記システムの状態ベクトルに基づいて前記任意のステップにおける前記システムの出力の推定値を算出する第２関数に、前記任意のステップにおける前記システムの状態ベクトルを入力して前記任意のステップにおける前記システムの出力の推定値を算出し、
前記パーティクルフィルタ実行処理では、複数のパーティクルのそれぞれに、前記システムの状態ベクトルを割り当て、前記観測値取得処理で取得した観測値と、パーティクル毎に前記推定値算出処理で取得した推定値とに基づいて、各パーティクルの重みを算出し、その重みに基づいてパーティクルを更新することで前記システムの状態ベクトルを推定して前記システムのパラメータを同定し、
前記パーティクルフィルタ実行処理では、
パーティクルは、前記任意の期間内の最初のステップからｎステップ目毎に（ｎ：２以上の整数）更新され、
各パーティクルの重みは、パーティクルが更新される１個の更新周期内に含まれるｎ個のステップのうちの少なくとも２個のステップにおいて取得された観測値と推定値とに基づいて算出される、コンピュータプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本明細書に開示する技術は、システムパラメータ同定装置、システムパラメータ同定方法、及びそのためのコンピュータプログラムに関する。詳しくは、制御対象となるシステムが有するパラメータを同定するための技術に関する。
【背景技術】
【０００２】
非特許文献１には、パーティクルフィルタを用いて制御対象となるシステムのパラメータを同定する技術が開示されている。パーティクルフィルタとは、所与の観測値と各パーティクルの推定値に基づいて各パーティクルの重みを求め、その重みに基づいてパーティクルを更新する手順を繰り返すことで、推定したい値を求める公知の手法である。
【先行技術文献】
【非特許文献】
【０００３】
【非特許文献１】益田哲也、杉江俊治、“粒子フィルタによる線形近似システムのパラメータ同定”、システム制御情報学会論文誌、ｖｏｌ．２４、Ｎｏ．１０、ｐ２５０−２５４、２０１１
【発明の概要】
【発明が解決しようとする課題】
【０００４】
非特許文献１では、システムの出力がステップ毎に観測され、パーティクルは、ステップ毎に更新される。即ち、パーティクルの更新周期は、システムの出力のステップ周期と等しい。この場合、各パーティクルの重みは、１個のステップにおける観測値と推定値とに基づいて算出される。このような手法では、実際のパラメータ（以下、真値とも称する）とは異なる値を有するパラメータに対してパーティクルの重みが大きくなる場合がある。このため、真値ではないパラメータを同定することが起こり得ることとなり、真値の同定精度が低い。
【０００５】
本明細書では、パーティクルフィルタを用いて制御対象となるシステムのパラメータを同定するに際して、パラメータの同定精度を向上できる技術を開示する。
【課題を解決するための手段】
【０００６】
本明細書は、パーティクルフィルタを用いて制御対象のシステムのパラメータを同定するシステムパラメータ同定装置を開示する。このシステムパラメータ同定装置は、入力値取得部と、観測値取得部と、第１関数記憶部と、第２関数記憶部と、状態ベクトル算出部と、推定値算出部と、パーティクルフィルタ実行部と、を備える。入力値取得部は、システムへの入力値を、任意の期間においてステップ毎に取得可能である。観測値取得部は、システムからの出力の観測値を、任意の期間においてステップ毎に取得可能である。第１関数記憶部は、任意の期間内の任意のステップにおけるシステムの状態ベクトルを、任意のステップの１つ前のステップにおけるシステムの状態ベクトルに基づいて算出する第１関数を記憶している。システムの状態ベクトルは、システムの状態量と、システムのパラメータと、入力値に関連するシステムノイズとを要素に含んでいる。第２関数記憶部は、任意のステップにおけるシステムの状態ベクトルに基づいて、任意のステップにおけるシステムの出力の推定値を算出する第２関数を記憶している。状態ベクトル算出部は、任意のステップの１つ前のステップにおけるシステムの状態ベクトルを第１関数記憶部に記憶されている第１関数に入力して、任意のステップにおけるシステムの状態ベクトルを算出する。推定値算出部は、任意のステップにおけるシステムの状態ベクトルを第２関数記憶部に記憶されている第２関数に入力して、任意のステップにおけるシステムの出力の推定値を算出する。パーティクルフィルタ実行部は、複数のパーティクルのそれぞれに、システムの状態ベクトルを割り当て、観測値取得部から取得した観測値と、パーティクル毎に推定値算出部から取得した推定値とに基づいて、各パーティクルの重みを算出し、その重みに基づいてパーティクルを更新することでシステムの状態ベクトルを推定してシステムのパラメータを同定する。パーティクルフィルタ実行部では、パーティクルは、任意の期間内の最初のステップからｎステップ目毎に（ｎ：２以上の整数）更新され、かつ、各パーティクルの重みは、パーティクルが更新される１個の更新周期内に含まれるｎ個のステップのうちの少なくとも２個のステップにおいて取得された観測値と推定値とに基づいて算出される。
【０００７】
上記のシステムパラメータ同定装置では、パーティクルの更新周期が、ステップ周期のｎ倍（ｎ：２以上の整数）である。そして、各パーティクルの重みは、パーティクルが更新される１個の更新周期内に含まれるｎ個のステップのうちの少なくとも２個のステップにおいて取得された観測値と推定値とに基づいて算出される。このため、各パーティクルの重みを１個のステップにおける観測値と推定値とに基づいて算出する場合と比較して、真値とは異なるパラメータに対してパーティクルの重みが大きくなってしまうことが起こり難くなる。パーティクルの重みは、より真値に近いパラメータに対して大きくなり易くなる。この結果、制御対象となるシステムのパラメータの同定精度を向上できる。
【０００８】
また、本明細書は、パーティクルフィルタを用いて制御対象のシステムのパラメータを同定する新規なシステムパラメータ同定方法を開示する。このシステムパラメータ同定方法は、入力値取得工程と、観測値取得工程と、状態ベクトル算出工程と、推定値算出工程と、パーティクルフィルタ実行工程とを備える。入力値取得工程では、システムへの入力値を、任意の期間においてステップ毎に取得する。観測値取得工程では、システムからの出力の観測値を、任意の期間においてステップ毎に取得する。状態ベクトル算出工程では、システムの状態量と、システムのパラメータと、入力値に関連するシステムノイズとを要素に含む、システムの状態ベクトルであって、任意の期間内の任意のステップにおけるシステムの状態ベクトルを、任意のステップの１つ前のステップにおけるシステムの状態ベクトルに基づいて算出する第１関数に、１つ前のステップにおけるシステムの状態ベクトルを入力して任意のステップにおけるシステムの状態ベクトルを算出する。推定値算出工程では、任意のステップにおけるシステムの状態ベクトルに基づいて任意のステップにおけるシステムの出力の推定値を算出する第２関数に、任意のステップにおけるシステムの状態ベクトルを入力して任意のステップにおけるシステムの出力の推定値を算出する。パーティクルフィルタ実行工程では、複数のパーティクルのそれぞれに、システムの状態ベクトルを割り当て、観測値取得工程で取得した観測値と、パーティクル毎に推定値算出工程で取得した推定値とに基づいて、各パーティクルの重みを算出し、その重みに基づいてパーティクルを更新することでシステムの状態ベクトルを推定してシステムのパラメータを同定する。パーティクルフィルタ実行工程では、パーティクルは、任意の期間内の最初のステップからｎステップ目毎に（ｎ：２以上の整数）更新され、かつ、各パーティクルの重みは、パーティクルが更新される１個の更新周期内に含まれるｎ個のステップのうちの少なくとも２個のステップにおいて取得された観測値と推定値とに基づいて算出される。このシステムパラメータ同定方法によると、制御対象となるシステムのパラメータの同定精度を適切に向上できる。
【０００９】
また、本明細書は、パーティクルフィルタを用いて制御対象のシステムのパラメータを同定するための新規なコンピュータプログラムを開示する。このコンピュータプログラムは、コンピュータに、入力値取得処理と、観測値取得処理と、状態ベクトル算出処理と、推定値算出処理と、パーティクルフィルタ実行処理と、を実行させる。入力値取得処理では、システムへの入力値を、任意の期間においてステップ毎に取得する。観測値取得処理では、システムからの出力の観測値を、任意の期間においてステップ毎に取得する。状態ベクトル算出処理では、システムの状態量と、システムのパラメータと、入力値に関連するシステムノイズとを要素に含む、システムの状態ベクトルであって、任意の期間内の任意のステップにおけるシステムの状態ベクトルを、任意のステップの１つ前のステップにおけるシステムの状態ベクトルに基づいて算出する第１関数に、１つ前のステップにおけるシステムの状態ベクトルを入力して任意のステップにおけるシステムの状態ベクトルを算出する。推定値算出処理では、任意のステップにおけるシステムの状態ベクトルに基づいて任意のステップにおけるシステムの出力の推定値を算出する第２関数に、任意のステップにおけるシステムの状態ベクトルを入力して任意のステップにおけるシステムの出力の推定値を算出する。パーティクルフィルタ実行処理では、複数のパーティクルのそれぞれに、システムの状態ベクトルを割り当て、観測値取得処理で取得した観測値と、パーティクル毎に推定値算出処理で取得した推定値とに基づいて、各パーティクルの重みを算出し、その重みに基づいてパーティクルを更新することでシステムの状態ベクトルを推定してシステムのパラメータを同定する。パーティクルフィルタ実行処理では、パーティクルは、任意の期間内の最初のステップからｎステップ目毎に（ｎ：２以上の整数）更新され、かつ、各パーティクルの重みは、パーティクルが更新される１個の更新周期内に含まれるｎ個のステップのうちの少なくとも２個のステップにおいて取得された観測値と推定値とに基づいて算出される。このコンピュータプログラムによると、コンピュータに、パラメータを同定するための処理を適切に実行させることができる。
【００１０】
本明細書が開示する技術の詳細、及び、さらなる改良は、発明を実施するための形態及び実施例にて詳しく説明する。
【図面の簡単な説明】
【００１１】
【図１】実施例１のシステムパラメータ同定装置と外部記憶媒体のブロック図。
【図２】システムパラメータの同定工程を示すフローチャート。
【図３】ロボットアームへの入力トルクの時系列データ。
【図４】ロボットアームからの出力角度の時系列データ。
【図５】クーロン摩擦力γのパラメータノイズζ３の分散β３２の分散計算式のグラフ。
【図６】粘性摩擦係数ｄのパラメータノイズζ４の分散β４２の分散計算式のグラフ。
【図７】慣性モーメントＪのパラメータノイズζ５の分散β５２の分散計算式のグラフ。
【図８】最初のステップ及びパーティクル更新時のステップにおける各パラメータノイズζ１〜ζ５の分散β１２〜β５２を表したグラフ。
【図９】実験１のパーティクル更新時のステップにおける各パラメータのパーティクルの分布を表したグラフ。
【図１０】実験２のパーティクル更新時のステップにおける各パラメータのパーティクルの分布を表したグラフ。
【図１１】実施例２の制御器とシステムのブロック図。
【発明を実施するための形態】
【００１２】
以下に説明する実施例の主要な特徴を列記しておく。なお、以下に記載する技術要素は、それぞれ独立した技術要素であって、単独であるいは各種の組合せによって技術的有用性を発揮するものであり、出願時請求項記載の組合せに限定されるものではない。
【００１３】
（特徴１）  本明細書が開示するシステムパラメータ同定装置は、システムのパラメータにパラメータノイズを付与するノイズ付与部をさらに備えていてもよい。このノイズ付与部は、パーティクルが更新されるステップにおいて、システムのパラメータにパラメータノイズを付与してもよい。この構成によると、パーティクルが局所解（ローカルミニマム）に収束することを抑制でき、同定精度をさらに向上できる。
【００１４】
（特徴２）  本明細書が開示するシステムパラメータ同定装置では、パラメータノイズの分散が可変であってもよい。この構成によると、パーティクルの拡散の度合いを変更できる。
【００１５】
（特徴３）  本明細書が開示するシステムパラメータ同定装置では、任意の期間内においてパーティクルがｍ回（ｍ：２以上の整数）更新されるとすると、ｊ回目（ｊ：１〜ｍ−１までの任意の整数）に更新されるときのステップで付与されるパラメータノイズの分散が、ｎ・ｊ＋１ステップ目からｎ・ｊ＋ｎステップ目までのｎ個のステップのうちの少なくとも１個のステップにおいて取得された入力値又は観測値に基づいて決定されてもよい。この構成によると、パーティクルが不要に拡散してしまうことを抑制でき、効率的にパラメータを同定できる。
【００１６】
（特徴４）  本明細書が開示するシステムパラメータ同定装置では、各パーティクルの重みは、上記少なくとも２個のステップにおける同時確率から算出されてもよい。この構成によると、パーティクルの重みは、より真値に近いパラメータに対して大きくなるため、同定精度をさらに向上できる。
【００１７】
（特徴５）  システムを制御する制御器が、上記に記載のシステムパラメータ同定装置を備えていてもよい。この制御器の制御パラメータは、パーティクルフィルタ実行部により更新されたパーティクルの分散及び平均に基づいて決定されてもよい。この構成によると、制御パラメータを効率的に決定することができる。
【実施例１】
【００１８】
図１〜３を参照してシステムパラメータ同定装置１０について説明する。同定装置１０は、パーティクルフィルタを用いて制御対象となるシステム（例えば、ロボットアーム）のパラメータＰを同定する。同定装置１０には、システムのパラメータＰを同定するための各種処理を実行するコンピュータが搭載されている。コンピュータは、演算処理を行うＣＰＵ、演算処理のデータが一時的に記憶されるＲＡＭ、及びＣＰＵによって実行される演算プログラムが記憶されたＲＯＭを備えている。ＣＰＵがＲＯＭに記憶された演算プログラムを実行することで、ＣＰＵは、後述する状態ベクトル算出部、推定値算出部等として機能する。
【００１９】
同定装置１０のコンピュータは、入力値取得部１２と、観測値取得部１４と、第１関数記憶部１６と、第２関数記憶部１８と、状態ベクトル算出部２０と、推定値算出部２２と、パーティクルフィルタ実行部２４と、分散計算式記憶部２６と、ノイズ決定部２８と、を備える。
【００２０】
同定装置１０は、オフラインで使用される。即ち、外部記憶媒体３０には、所定期間（ｔ＝０〜Ｔ）における入力値データ３２と観測値データ３４とが予め記憶されている。入力値データ３２は、システムへの入力値を示す時系列データである。観測値データ３４は、入力に対するシステムからの出力の観測値を示す時系列データである。
入力値取得部１２は、外部記憶媒体３０に記憶されている入力値データ３２から、ステップ毎に（即ち、一定の周期で）入力値ｕを取得する。取得された入力値は、ノイズ決定部２８及び状態ベクトル算出部２０に出力される。この処理は、ＣＰＵで実施される。以下では、ｋステップ目（ｋ＝０〜Ｔ／Ｔｓ、Ｔｓ：ステップ周期）の入力値ｕをｕｋと表し、他の値の表記についても同様とする。
観測値取得部１４は、外部記憶媒体３０に記憶されている観測値データ３４から、ステップ毎に観測値Ｙｋｒを取得する。取得された観測値Ｙｋｒは、ノイズ決定部２８及びパーティクルフィルタ実行部２４に出力される。この処理は、ＣＰＵで実施される。なお、入力値取得部１２及び観測値取得部１４は、所定期間の一部の期間における入力値データ３２及び観測値データ３４を取得してもよい。
【００２１】
第１関数記憶部１６は、ＲＯＭに設けられており、式３で表される第１関数を記憶している。第１関数は、制御対象となるシステムの状態方程式と、隣接ステップ間のパラメータＰの関係式を統合した関数であり、詳細には、以下のように定義される。
【００２２】
制御対象となるシステムの状態方程式は、式１によって定義される。
【数１】
2017083922.tif 000003
ｆは既知の非線形関数、Ｘｋはｋステップ目のシステムの状態量、ηｋはｋステップ目のシステムノイズである。式１から明らかなように、ｋ＋１ステップ目の状態量Ｘｋ＋１は、ｋステップ目における状態量Ｘｋ及びシステムノイズηｋに基づいて算出される。なお、本実施例では、システムノイズηｋとして、平均が入力値ｕｋ、分散がα２（定数）であるシステムノイズが用いられる。
【００２３】
隣接ステップ間のパラメータＰの関係式は、式２によって定義される。
【数２】
2017083922.tif 000004
ζｋは、ｋステップ目のパラメータＰｋに対するパラメータノイズである。ζｋの平均はゼロ、分散はβ２である。本実施例では分散β２は可変であるが、定数であってもよい。分散β２については後で詳述する。
【００２４】
式１と式２を統合することにより、式３で表される新たな状態方程式（第１関数）が定義される。
【数３】
2017083922.tif 000005
Ｚｋは、ｋステップ目のシステムの状態ベクトルであり、Ｚｋ＝（ＸｋＴ，ＰｋＴ）Ｔとして定義される。ｆａは式１、２から導出される既知の非線形関数である。状態ベクトルＺｋは、状態量Ｘｋ、パラメータＰｋ、システムノイズηｋ、パラメータノイズζｋを要素に含んでいる。また、式３から明らかなように、ｋ＋１ステップ目の状態ベクトルＺｋ＋１は、ｋステップ目における状態ベクトルＺｋ、システムノイズηｋ、及びパラメータノイズζｋに基づいて算出される。
【００２５】
第２関数記憶部１８は、ＲＯＭに設けられており、式４で表される出力方程式（以下、第２関数とも称する）を記憶している。
【数４】
2017083922.tif 000006
Ｙｋはｋステップ目のシステムの出力の推定値、Ｈｋはｋステップ目の既知の行列、σｋはｋステップ目の観測ノイズである。σｋの平均は予め設定されており、分散はρ２（定数）である。式４から明らかなように、ｋステップ目の推定値Ｙｋは、ｋステップ目における状態ベクトルＺｋ及び観測ノイズσｋに基づいて算出される。
【００２６】
状態ベクトル算出部２０は、後述するパーティクルフィルタ実行部２４で発生させたＳ個のパーティクルのそれぞれについて、ｋステップ目における状態ベクトルＺｋｉ（ｉ＝１〜Ｓ）と、入力値取得部１２から取得した入力値ｕｋに所定の分散α２を持たせたシステムノイズηｋと、ノイズ決定部２８から取得したパラメータノイズζｋを、第１関数記憶部１６に記憶されている第１関数（式３参照）に入力して、ｋ＋１ステップ目の状態ベクトルＺ(〜)ｋ＋１ｉを算出する（なお、実際の表記ではチルダ「〜」は文字の上に付される）。状態ベクトル算出部２０で算出された状態ベクトルＺ(〜)ｋ＋１ｉは、推定値算出部２２に出力される。この処理はＣＰＵで実施される。なお、本実施例では、第１関数に入力して算出された状態ベクトルにはチルダ「〜」を付すことで、パーティクルフィルタ実行部２４で推定された状態ベクトルと区別する。
【００２７】
推定値算出部２２は、各パーティクルについて、状態ベクトル算出部２０から取得したｋ＋１ステップ目の状態ベクトルＺ(〜)ｋ＋１ｉを、第２関数記憶部１８に記憶されている第２関数（式４参照）に入力して、ｋ＋１ステップ目のシステムの出力の推定値Ｙ(〜)ｋ＋１ｉを算出する。このとき、観測ノイズσｋはゼロとして推定値Ｙ(〜)ｋ＋１ｉが算出される。推定値算出部２２で算出された推定値Ｙ(〜)ｋ＋１ｉは、パーティクルフィルタ実行部２４に出力される。この処理はＣＰＵで実施される。なお、本実施例では、第２関数に入力して算出された推定値にはチルダ「〜」を付す。
【００２８】
パーティクルフィルタ実行部２４は、最初のステップ（即ち、ｋ＝０）において、Ｓ個のパーティクルを発生させる。パーティクルは所定の範囲で一様に分布している。ここで、「所定の範囲」とは、同定すべきパラメータＰの真値を含む範囲である。各パーティクルには、所定の状態ベクトルＺ０ｉが割り当てられる。以下では、最初のステップにおける状態ベクトルＺ０ｉを、初期状態ベクトルとも称する。
また、パーティクルフィルタ実行部２４は、最初のステップからｎステップ目毎に、パーティクルを更新する。具体的には、パーティクルフィルタ実行部２４は、以下の手順でパーティクルを更新する。（１）まず、パーティクルの更新がｊ回目であるｋ＝ｎ・ｊステップ目（ｊ＝１〜ｍ、ｍ：所定期間Ｔ内でパーティクルが更新される回数、ｎ・ｍ＝Ｔ／Ｔｓ）において、観測値取得部１４から、ｊ−１回目のパーティクルの更新後に取得されたｎ個の観測値Ｙｒ（即ち、Ｙｎ（ｊ−１）＋１ｒからＹｎ（ｊ−１）＋ｎｒまでのｎ個の観測値Ｙｒ）を取得する。また、推定値算出部２２から、各パーティクルについて、ｊ−１回目のパーティクルの更新後に取得されたｎ個の推定値Ｙ(〜)ｉ（即ち、Ｙ(〜)ｎ（ｊ−１）＋１ｉからＹ(〜)ｎ（ｊ−１）＋ｎｉまでのｎ個の推定値Ｙ(〜)ｉ）を取得する。（２）次に、これらｎ個の観測値Ｙｒ及び推定値Ｙ(〜)ｉを次の式５で表される重み計算式に入力して、各パーティクルの重みｗｉを算出する。
【数５】
2017083922.tif 000007
なお、式５のＷは次の式６で定義される。
【数６】
2017083922.tif 000008
（３）続いて、パーティクル数が、重みの高い領域で多く、重みの低い領域で少なくなるようにパーティクルを更新して、ｋ＝ｎ・ｊステップ目の状態ベクトルＺｎ・ｊｉを推定して、ｋ＝ｎ・ｊステップ目のシステムパラメータＰｎ・ｊを同定する。なお、この状態ベクトルＺｎ・ｊｉは、ｎ個の観測値Ｙｒ及び推定値Ｙ(〜)ｉに基づいて推定されたものであり、第１関数から算出されたｋ＝ｎ・ｊステップ目の状態ベクトルＺ(〜)ｎ・ｊｉとは異なる。
パーティクルフィルタ実行部２４で推定された状態ベクトルＺｎ・ｊｉは、状態ベクトル算出部に出力される。これらの処理は、ＣＰＵで実施される。
【００２９】
分散計算式記憶部２６は、ＲＯＭに設けられており、パラメータノイズζｋの分散β２を記述する計算式を記憶している。この計算式は、入力値ｕ又は観測値Ｙｒの関数であり、パラメータ毎に予め設定されている（後述）。なお、この計算式は、入力値ｕ又は観測値Ｙｒの関数に限られず、入力値ｕ又は観測値Ｙｒから導出される値（例えばＹｒ）の関数であってもよい。
【００３０】
ノイズ決定部２８は、ｋ＝ｎ・ｊステップ目（但し、ｊ≠ｍ）において、入力値取得部１２又は観測値取得部１４から、ｊ回目のパーティクル更新後の最初のステップからｊ＋１回目のパーティクル更新時のステップまでのｎ個の入力値ｕ又はｎ個の観測値Ｙｒを取得する（即ち、ｕｎ・ｊ＋１からｕｎ・ｊ＋ｎまでのｎ個の入力値ｕ又はＹｎ・ｊ＋１ｒからＹｎ・ｊ＋ｎｒまでのｎ個の観測値Ｙｒを取得する）。そして、分散式記憶部から取得した分散式にｎ個の入力値ｕ又はｎ個の観測値Ｙｒから選択された１個の入力値ｕ又は観測値Ｙｒを入力して、分散β２を算出し、ｋ＝ｎ・ｊステップ目のパラメータノイズζｎ・ｊを決定する。一方、ｋ＝ｎ・ｊ以外のステップでは、上記処理は実行されず、パラメータノイズζは一律にゼロと設定される。このようにして決定されたパラメータノイズζｋは、状態ベクトル算出部２０に出力される。この処理はＣＰＵで実施される。
上述したことから明らかなように、パーティクルの更新が行われないステップ（ｋ＝ｎ・ｊ以外のステップ）では、ζｋ＝０である。即ち、パーティクルの更新が行われないステップでは、Ｐｋ＋１＝Ｐｋが成立する。このため、ノイズ決定部２８は、パーティクルの更新を行うステップにおいてのみパラメータノイズを決定し、状態ベクトル算出部２０に出力するということもできる。
なお、上記の構成では、ノイズ決定部２８はパーティクルの更新を行うステップにおいてのみパラメータノイズを決定したが、この構成に限られない。ノイズ決定部２８は、最初のステップ（ｋ＝０）においても、パラメータノイズζ０を決定してもよい。この場合、パラメータノイズζ０の分散β２は、ｋ＝１〜ｎステップ目までのｎ個の入力値ｕ又はｎ個の観測値Ｙｒから選択された１個の入力値ｕ又は観測値Ｙｒを分散式に入力することにより算出される。なお、ノイズ決定部２８が「ノイズ付与部」の一例に相当する。
【００３１】
次に、図２を参照して同定装置１０のコンピュータがシステムのパラメータＰを同定する処理について説明する。まず、ステップＳ２の入力値取得工程では、入力値取得部１２が、外部記憶媒体３０に記憶されている入力値データ３２から、ステップ毎に全期間の入力値ｕｋを取得する。次に、ステップＳ４の観測値取得工程では、観測値取得部１４が、外部記憶媒体３０に記憶されている観測値データ３４から、ステップ毎に全期間の観測値Ｙｋｒを取得する。続いて、ステップＳ６の初期状態ベクトル決定工程では、パーティクルフィルタ実行部２４が、最初のステップ（ｋ＝０）においてＳ個のパーティクルを発生させ、各パーティクルに初期状態ベクトルＺ０ｉ（ｉ＝１〜Ｓ）を割り当てる。なお、ステップＳ２〜Ｓ６の工程は、順不同である。
【００３２】
続いて、ステップＳ８の状態ベクトル算出工程では、状態ベクトル算出部２０が、ステップＳ６で割り当てられたパーティクル毎の初期状態ベクトルＺ０ｉと、システムノイズη０（即ち、平均が入力値ｕ０、分散がα２となるノイズ）と、ζ０＝０を第１関数に入力して、状態ベクトルＺ(〜)１ｉを算出する。次いで、状態ベクトル算出部２０は、Ｚ(〜)１ｉと、η１と、ζ１＝０を第１関数に入力して、Ｚ(〜)２ｉを算出する。状態ベクトル算出部２０は、この処理をｋ＝ｎステップ目までｎ回繰り返して、パーティクル毎にｎ個の状態ベクトルＺ(〜)１ｉ〜Ｚ(〜)ｎｉを算出する。なお、ｋ＝０〜ｎ−１ステップ目までのζｋは０である。続いて、ステップＳ１０の推定値算出工程では、推定値算出部２２が、ステップＳ８で算出されたｎ個の状態ベクトルＺ(〜)１ｉ〜Ｚ(〜)ｎｉを第２関数に入力して、パーティクル毎にｎ個の推定値Ｙ(〜)１ｉ〜Ｙ(〜)ｎｉを算出する。
【００３３】
続いて、ステップＳ１２のパーティクルフィルタ実行工程では、１回目のパーティクルの更新を行う。具体的には、パーティクルフィルタ実行部２４が、ステップＳ４で取得した全期間の観測値Ｙｋｒのうちの、ｋ＝１〜ｎステップ目までのｎ個の観測値Ｙ１ｒ〜Ｙｎｒと、ステップＳ１０で取得したｎ個の推定値Ｙ(〜)１ｉ〜Ｙ(〜)ｎｉを重み計算式（式５参照）に入力して、各パーティクルの重みｗｉを算出する。そして、パーティクルフィルタ実行部２４は、重みｗｉに基づいてパーティクルを更新して、パーティクル毎にｋ＝ｎステップ目の状態ベクトルＺｎｉを推定する。このときのパーティクルの分布により、ｋ＝ｎステップ目のパラメータＰが同定される。
【００３４】
続いて、ステップＳ１４では、パーティクルの更新が最後であるか（即ち、ｊ＝ｍであるか）否かを判定する。パーティクルの更新が最後ではない場合（ステップＳ１４でＮＯ）は、ステップＳ１６に進む。ステップＳ１６の分散決定工程では、ノイズ決定部２８が、ｋ＝ｎ＋１〜２ｎステップ目までのｎ個の入力値ｕｎ＋１〜ｕ２ｎ、又は、ｎ個の観測値Ｙｎ＋１ｒ〜Ｙ２ｎｒ（或いは、入力値ｕｋ又は観測値Ｙｋｒから導出される値）から選択された１個の入力値ｕｋ又は観測値Ｙｋｒを、分散計算式記憶部２６に記憶されている計算式に入力して分散β２を算出し、ｋ＝ｎステップ目のパラメータノイズζｎを決定する。
【００３５】
続いて、ステップＳ１８の状態ベクトル算出工程では、状態ベクトル算出部２０が、ステップＳ１２で推定されたｋ＝ｎステップ目の状態ベクトルＺｎｉと、ηｎと、ステップＳ１６で決定されたパラメータノイズζｎを、第１関数に入力してｋ＝ｎ＋１ステップ目の状態ベクトルＺ(〜)ｎ＋１ｉを算出する。ステップＳ１８の処理が終了すると、同定装置１０のコンピュータは、ステップＳ８の処理に戻る。ステップＳ８では、状態ベクトル算出部２０が、ｋ＝ｎ＋２〜２ｎステップ目までのｎ−１個の状態ベクトルＺ(〜)ｎ＋２ｉ〜Ｚ(〜)２ｎｉを算出する。続いて、ステップＳ１０では、ステップＳ１８及びステップＳ８で算出されたｎ個の状態ベクトルＺ(〜)ｎ＋１ｉ〜Ｚ(〜)２ｎｉから、ｎ個の推定値Ｙ(〜)ｎ＋１ｉ〜Ｙ(〜)２ｎｉが算出され、ステップＳ１２では、２回目のパーティクルの更新が行われる。以下、上述した処理を繰り返し、ステップＳ１２でｍ回目（即ち、最後）のパーティクルの更新が行われてｋ＝ｎ・ｍステップ目の状態ベクトルＺｎ・ｍｉが推定されると、同定装置１０のコンピュータは、ステップＳ１４でＹＥＳと判定し、パラメータＰの同定処理を終了する。このときのパーティクルの分布により、ｋ＝ｎ・ｍステップ目のパラメータＰが同定される。
【００３６】
上記の説明から明らかなように、パーティクルが更新されるステップ（ｋ＝ｎ・ｊ、ｊ＝１〜ｍ）では、状態ベクトルは、状態ベクトル算出部２０によって算出される状態ベクトルＺ(〜)ｎｉと、パーティクルフィルタ実行部２４によって推定される状態ベクトルＺｎｉの２種類が存在する。そして、ステップＳ１８の工程では、第１関数に入力される状態ベクトルとして、後者の状態ベクトルＺｎｉが用いられる。なお、上記のようにパラメータの同定がオフラインで行われる場合は、ステップＳ１６の工程のように都度パラメータノイズζｎ・ｊを決定する代わりに、予めｍ−１個の分散β２を算出してパラメータノイズζｎ・ｊを決定しておいてもよい。
【００３７】
上記の同定装置１０では、パーティクルの更新周期が、ステップ間隔Ｔｓのｎ倍である。そして、各パーティクルの重みｗｉは、パーティクルが更新される１個の更新周期内に含まれるｎ個のステップにおいて取得された観測値Ｙｒと推定値Ｙ(〜)ｉとに基づいて算出される。このため、各パーティクルの重みを１個のステップにおける観測値Ｙｒと推定値Ｙ(〜)ｉとに基づいて算出する場合と比較して、パーティクルの重みｗｉは、より真値に近いパラメータに対して大きくなり易くなる。この結果、制御対象となるシステムのパラメータの同定精度を向上できる。
【００３８】
また、上記の同定装置１０は、ノイズ決定部２８を備える。ノイズ決定部２８は、パーティクルが更新されるステップ（ｋ＝ｎ・ｊ、但し、ｊ≠ｍ）において、パラメータノイズζを決定して状態ベクトル算出部２０に出力する。この構成によると、パーティクルが局所解（ローカルミニマム）に収束することを抑制でき、同定精度をさらに向上できる。
【００３９】
また、上記の同定装置１０では、ノイズ決定部２８はパラメータノイズζの分散β２を、分散計算式記憶部２６に記憶されている分散計算式を用いて決定するため、分散β２は可変となる。この構成によると、分散β２に応じてパーティクルの拡散の度合いを変更できる。
【００４０】
特に、上記の同定装置１０では、パーティクルのｊ回目の更新時のステップで第１関数に入力されるパラメータノイズζｎ・ｊの分散β２が、ｋ＝ｎ・ｊ＋１ステップ目からｎ・ｊ＋ｎステップ目までのｎ個のステップのうちの１個のステップにおいて取得された入力値ｕ又は観測値Ｙｒに基づいて決定される。この構成によると、パラメータを同定するための情報の多寡に応じて分散β２を決定できる。このため、情報が少ないときに大きな分散β２を設定することに起因してパーティクルが不要に拡散してしまうという事態の発生を抑制でき、効率的にパラメータを同定できる。
【００４１】
また、上記の同定装置１０では、各パーティクルの重みｗｉが、ｎ個のステップにおける同時確率から算出される。このため、パーティクルの重みｗｉは、より真値に近いパラメータに対して大きくなり、同定精度をさらに向上できる。
【００４２】
次に、ロボットアームのシミュレーションモデルを用いて同定装置のパラメータ同定精度を検証した具体例について説明する。この具体例では、同定装置は、実験１と、実験１の比較例としての実験２を実施する。実験１では、上述したノイズ決定処理が行われる（即ち、パラメータノイズζの分散β２が可変である）。実験２では、ノイズ決定処理が行われない（即ち、パラメータノイズζの分散β２が定数である）。このロボットアームは、単軸運動モデルであり、ロボットアームの軸を回転軸として回転運動する。同定装置１０が同定するロボットアームのパラメータＰは、ｍ、φ、γ、ｄ、Ｊの５個である。ｍは転がり摩擦のヒステリシスの膨らみを表す係数、φは転がり出し変位領域幅、γはクーロン摩擦力、ｄは粘性摩擦係数、Ｊは慣性モーメントを表す。
【００４３】
（実験１）
まず、ロボットアームの状態方程式と、隣接ステップ間のパラメータＰの関係式を統合して、ロボットアームの第１関数を求める。ロボットアームの状態方程式は、次の式７で示すロボットアームの運動方程式を以下の手順で離散化することにより求められる。
【数７】
2017083922.tif 000009
θはロボットアームの回転角度、τはロボットアームに入力されるトルク、ｖはクーロン摩擦を包含する転がり摩擦力である。ここで、ｖ（ｔ）は次の式８によって定義される。
【数８】
2017083922.tif 000010
ωはロボットアームの角速度（＝θ）、ν（ニュー）は速度反転後の転がり摩擦力、δは速度反転後の移動距離、ξ＝δ／φである。ここで、ｇ（ξ）は次の式９によって定義される。
【数９】
2017083922.tif 000011
【００４４】
ここで、上述したように、式７〜９をパーティクルフィルタに実装するためには、これらを離散系に変換する必要がある。離散系に変換するためには、まず、速度反転を判定して、次の式１０、１１に従って速度反転後のころがり摩擦力νｋと速度反転後の移動距離δｋを定める。
【数１０】
2017083922.tif 000012
【数１１】
2017083922.tif 000013
式１０、１１を式８に代入すると、離散系での転がり摩擦力ｖｋが次の式１２の通りに導出できる。但し、ｇ（ξ）には式９、ξにはδｋ／φｋが用いられる。
【数１２】
2017083922.tif 000014
【００４５】
一方、式７を状態方程式に書き換えると、次の式１３となる。
【数１３】
2017083922.tif 000015
式１３を簡略化すると、次の式１４となる。
【数１４】
2017083922.tif 000016
式１４を、入力τ（ｔ）にゼロ次ホールドを仮定して離散化すると、次の式１５が得られる。但し、ｖ（ｔ）の離散化は、式１２に従っている。
【数１５】
2017083922.tif 000017
式１５を一般式に置き換えると次の式１６で示すロボットアームの状態方程式が求められる。
【数１６】
2017083922.tif 000018
ここで、ｆは非線形な転がり摩擦を包含した式１５を表す既知の非線形関数、Ｘｋ’はＸｋと転がり摩擦の状態量ωｋ−１、θｋ−１、δｋ−１、νｋ−１、ｖｋ−１を含有した状態行列、ηｋは平均τｋ、分散α２（＝２．５×１０−３）の入力信号であり、いわゆるシステムノイズに相当する。
【００４６】
一方、隣接ステップ間のパラメータＰの関係式は、式２で示した通りである。但し、Ｐｋは、このロボットアームのシミュレーションモデルで同定すべき５個のパラメータからなる行列であり、次の式１７で表される。
【数１７】
2017083922.tif 000019
【００４７】
式１６と式２を統合することにより、式３で表される新たな状態方程式（第１関数）が定義される。但し、Ｚｋ＝（Ｘｋ’Ｔ，ＰｋＴ）Ｔとして定義される。ｆａは式１６、式２から導出される既知の非線形関数である。
【００４８】
一方、観測量は角度θｋであるため、出力方程式（第２関数）は線形となり、次の式１８で表される。
【数１８】
2017083922.tif 000020
本モデルでは、観測ノイズσｋは、平均ゼロ、分散ρ２＝４．０×１０−６に設定されている。
【００４９】
図３は、ロボットアームに入力するトルクτ（ｔ）の時系列データを示し、図４は、ロボットアームから出力される角度θ（ｔ）の時系列データを示す。前者が入力値データ３２に相当し、後者が観測値データ３４に相当する。時系列データの期間ＴはＴ＝４５秒であり、ステップ間隔ＴｓはＴｓ＝０．０１秒である。同定装置１０の入力値取得部１２及び観測値取得部１４は、時系列データの全期間の入力値及び観測値を取得する。パーティクルは、ｎ＝２５０ステップ目に更新される。即ち、パーティクルの更新周期は２．５秒であり、パーティクルは計１８回更新される（ｊ＝１〜１８）。また、パーティクルの総数ＳはＳ＝５００００とする。また、パラメータの真値は、ｍ＝０．８、φ＝１．５、γ＝５．０、ｄ＝２．０、Ｊ＝４．０に設定されている。
【００５０】
次に、分散計算式記憶部２６が記憶している分散計算式について説明する。本モデルでは、同定すべきパラメータが５個あるため、パラメータノイズζｋは、各パラメータに対応する５個のパラメータノイズ及び分散を有する。以下では、パラメータｍｋ、φｋ、γｋ、ｄｋ、Ｊｋに対応するパラメータノイズ及び分散を、それぞれζ１ｋ及びβ１２、ζ２ｋ及びβ２２、ζ３ｋ及びβ３２、ζ４ｋ及びβ４２、ζ５ｋ及びβ５２とする。
【００５１】
（分散β１２、β２２の分散計算式）
転がり摩擦力ｖｋは、速度反転後に転がり出し変位領域幅φｋの距離を移動すると、クーロン摩擦力γと等価となる。このため、転がり摩擦のヒステリシスの膨らみを表す係数ｍｋとφｋは、速度反転後にのみ転がり摩擦力ｖｋに影響する同定パラメータである。従って、パラメータｍｋ、φｋのパラメータノイズζ１ｋ、ζ２ｋの分散β１２、β２２の分散計算式は、それぞれ次の式１９、式２０で表される。
【数１９】
2017083922.tif 000021
【数２０】
2017083922.tif 000022
なお、速度反転の有無は、観測値データθ（ｔ）（図３参照）から判定される。
【００５２】
（分散β３２の分散計算式）
パラメータγｋ（クーロン摩擦力）のパラメータノイズζ３ｋの分散β３２の分散計算式は、式２１で表される（図５参照）。なお、β３２ｍａｘ＝２．５×１０−３である。ｋ＝ｎ・ｊステップ目の分散β３２は、ｋ＝ｎ・ｊ＋１〜ｎ・ｊ＋ｎステップ目までの観測値データθ（ｔ）の最大角速度の絶対値ｍａｘ（｜ω｜）を式２１に入力することにより算出される。ロボットアームの高速回転時はクーロン摩擦力に対して粘性摩擦力が支配的になるため、クーロン摩擦力γｋのパラメータノイズζ３ｋの分散β３２が小さくなるように設定している。また、ロボットアームの極低速回転時は、速度反転直後である可能性が高いと判断し、この場合も分散β３２が小さくなるように設定している。
【数２１】
2017083922.tif 000023
【００５３】
（分散β４２の分散計算式）
粘性摩擦力は角速度に比例した摩擦力である。このため、角速度が小さいときは、ロボットアームの挙動に対する粘性摩擦力の影響が小さくなり、システム同定が困難となると考えられる。従って、パラメータｄｋ（粘性摩擦係数）のパラメータノイズζ４ｋの分散β４２が、最大角速度の絶対値ｍａｘ（｜ω｜）に比例するように、分散β４２の分散計算式を次の式２２のように定める（図６参照）。なお、β４２ｍａｘ＝２．５×１０−３である。
【数２２】
2017083922.tif 000024
【００５４】
（分散β５２の分散計算式）
パラメータＪｋ（慣性モーメント）は、ロボットアームが十分な加減速動作をするときにシステム同定が可能である。別言すれば、ロボットアームの定速回転時や停止時には同定ができない。このため、慣性モーメントＪｋのシステムノイズζ５ｋの分散β５２がトルクτ（入力値）に比例するように、分散β５２の分散計算式を次の式２３のように定める（図７参照）。ｋ＝ｎ・ｊステップ目の分散β５２は、ｋ＝ｎ・ｊ＋１〜ｎ・ｊ＋ｎステップ目までの入力値データτ（ｔ）の最大トルクの絶対値ｍａｘ（｜τ｜）を式２３に入力することにより算出される。なお、β５２ｍａｘ＝１．０×１０−２である。
【数２３】
2017083922.tif 000025
【００５５】
図８（ａ）〜（ｅ）は、最初のステップ（ｋ＝０）と、パーティクル更新時のステップ（ｋ＝ｎ・ｊ（ｊ＝１〜１７））におけるパラメータノイズζの分散を表す。横軸Ｎは、最初のステップから数えたときの分散β１２〜β５２の個数を表し、縦軸β１２〜β５２は分散の値を表す。本実験では、最初のステップにおいてもパラメータノイズζの分散が算出されるため、計１８個の分散が算出される。図８（ａ）〜（ｅ）から明らかなように、分散β１２〜β５２は可変である。例えば、図８（ａ）において、Ｎ＝１の分散β１２は、ｋ＝０ステップ目のパラメータノイズζ１・０の分散であり、Ｎ＝１８の分散β１２は、ｋ＝２５０・１７ステップ目のパラメータノイズζ１・２５０・１７の分散である。図９（ａ）〜（ｅ）は、上述した条件で同定装置１０がロボットアームのパラメータ同定処理を実行したときの各パラメータのパーティクルの分布を示す。図１０において、濃色はパーティクルの数が少ない領域を示し、淡色はパーティクルの数が多い領域を示している。また、破線はパラメータの真値を示し、丸印は、パーティクル更新時のステップにおけるパーティクルの平均値を示す。
【００５６】
（実験２）
図１０（ａ）〜（ｅ）は、各パラメータのパラメータノイズζｋの分散β２を定数とした場合の各パラメータのパーティクルの分布を示す。なお、この場合、各パラメータｍｋ、φｋ、γｋ、ｄｋ、Ｊｋのパラメータノイズζｋの分散を、それぞれβ６２、β７２、β８２、β９２、β１０２とすると、各分散は、β６２＝β１ｍａｘ２＝０．０２２、β７２＝β２ｍａｘ２＝０．０５２、β８２＝β３ｍａｘ２＝２．５×１０−３、β９２＝β４ｍａｘ２＝２．５×１０−３、β１０２＝β５ｍａｘ２＝１．０×１０−２と設定されている。
【００５７】
表１は、実験１、２を実施したときの最後のパーティクル更新時のステップにおけるパーティクルの平均と分散を各パラメータについて示している。表２は、実験１、２を各１０回実施したときの最後のパーティクル更新時のステップにおけるパーティクルの平均の誤差率（％）と分散の平均を各パラメータについて示している。
【表１】
2017083922.tif 000026
【表２】
2017083922.tif 000027
図９、１０、表１、２によると、実験１、２共にパーティクルの平均はパーティクルを更新する度に徐々に真値に近い値に収束している。しかしながら、実験１のほうが、パラメータＪを除いてパーティクルの平均が真値により近く、かつ、平均の誤差率も小さい。パラメータＪに関しても、実験１と２の平均の差及び平均の誤差率の差はごく僅かである。また、パラメータφを除いて、実験１のほうが、パーティクルの分散の平均は大幅に小さくなっている。パラメータφに関しても、実験１と２の分散の平均の差はごく僅かである。このことから、パラメータノイズの分散を可変としたほうが、分散を定数とする場合よりも、各回のシミュレーションの結果がばらつくことを抑制しながら、パラメータの同定精度をより向上できることが分かる。
【００５８】
また、実験１に関して、表１、表２の結果によると、パラメータφ、γ、ｄ、Ｊの平均の誤差率はいずれも０．５％未満であるのに対し、パラメータｍの平均の誤差率は１．８２９％と他のパラメータと比べてやや大きい。これは、パラメータｍの変化に対してθｋが不感であることが原因と考えられる。このため、同定結果をセンサレス力制御等の制御に応用する場合には、それほど大きな影響を与えないものと思われる。
【００５９】
また、図９、１０を比較すると、実験２では、実験１よりも、パーティクルの更新を繰り返す過程でパーティクルの平均が真値から大きく逸脱していることが分かる。これは、パラメータを同定するための十分な情報がない場合であっても、一定の分散を有するパラメータノイズζが同定すべきパラメータに与えられるためである。特に、実験２では、パラメータＪのパーティクルの平均が、ｊ＝１３、１４で真値から大幅に逸脱している。この原因は、次のように推測される。即ち、図３、４によると、ｊ＝１３、１４のステップ（即ち、ｔ＝３０〜３２．５秒）では、その約７．５秒前から入力トルクτがほぼゼロであり、出力角度θの角加速度がほぼゼロである。このように、ロボットアームの運動が角加速度を殆ど発生しない場合に、実験２のように一定の分散を有するパラメータノイズζをパラメータに与えると、パーティクルが不要に拡散して、真値ではない値にパーティクルが集合してしまうことが原因であると推測される。実験１のように、パラメータを同定するための十分な情報がある場合とない場合でパラメータノイズζの分散β２を可変とすることにより、パーティクルの更新の終了時だけではなく、パーティクルの更新を繰り返す過程においても、高い精度でパラメータを同定できることが分かる。
【実施例２】
【００６０】
次に、図１２を参照して実施例２について説明する。本実施例では、制御器４０が、制御対象となるシステム４２と通信可能に接続されている。制御器４０は、システムパラメータ同定装置４４と制御部４６を備える。制御部４６は、システム４２からの制御出力値を目標値に近づけるように、システム４２からの制御出力値を用いたフィードバック制御を行うことによってシステム４２への制御入力値を決定する。制御部４６は、システム４２への制御入力値を決定するための制御パラメータを有しており、制御入力値は、制御パラメータに基づいて決定される（後述）。
【００６１】
同定装置４４は、実施例１の同定装置１０と異なり、システムパラメータをオンラインで同定する。即ち、同定装置４４の入力値取得部及び観測値取得部は、システム４２の制御入力値及び制御出力値を、システム４２の動作中にリアルタイムで取得する。同定装置４４は、ノイズ決定処理を実施しない点を除いて、実施例１の同定装置１０と同様の処理を実施する。即ち、本実施例では、ｋ＝ｎ・ｊステップ目（ｊ＝１〜ｍ−１）におけるシステムパラメータのパラメータノイズζの分散β２は、可変ではなく一定とされる。同定装置４４のパーティクルフィルタ実行部２４により、ｋ＝ｎ・ｊステップ目（ｊ＝１〜ｍ）毎にパーティクルの更新が行われると、更新後のパーティクルの平均及び分散が、制御部４６に出力される。
【００６２】
本実施例では、制御パラメータを算出するための制御パラメータ計算式が、パーティクルの平均、パーティクルの分散、及び制御出力値と制御入力値との偏差を変数とする関数となるように設計されている。このため、制御部４６は、ｋ＝ｎ・ｊステップ目において同定装置４４から出力されたパーティクルの平均及び分散と、制御出力値と制御入力値との偏差とを制御パラメータ計算式に入力して、ｋ＝ｎ・ｊステップ目の制御パラメータを算出する。制御部４６は、この制御パラメータに基づいて、ｋ＝ｎ・ｊ＋１ステップ目の制御入力値を決定する。制御器４０は、上記の処理をｋ＝ｎ・ｊステップ目毎に繰り返すことにより、制御パラメータを同定する。
【００６３】
この構成によると、制御器４０は、システム４２のフィードバック制御とシステムパラメータの同定を同時に行う。制御器４０は、同定装置４４によるオンラインでのシステムパラメータの同定結果（具体的には、パーティクルの平均及び分散）を利用して、制御パラメータを決定する。この構成によると、制御部４６が試行錯誤的に制御パラメータを決定する構成と比較して、効率的に正確な制御パラメータを決定できる。
【００６４】
以上、本発明の具体例を詳細に説明したが、これらは例示にすぎず、特許請求の範囲を限定するものではない。特許請求の範囲に記載の技術には、以上に例示した具体例を様々に変形、変更したものが含まれる。
【００６５】
例えば、同定装置は、各パーティクルの重みｗｉを、１個の更新周期内に含まれるｎ個のステップのうちの少なくとも２個のステップにおいて取得された観測値Ｙｒと推定値Ｙ(〜)ｉとに基づいて算出してもよい。また、各パーティクルの重みｗｉは、同時確率以外の方法（即ち、式５以外の重み計算式）で求めてもよい。また、分散α２、ρ２は、ステップ毎に変化するように予め設定してもよい。また、パラメータノイズζ及び観測ノイズσの平均はゼロでなくてもよい。また、ノイズ決定部２８は、パーティクルが更新される全てのステップにおいてパラメータノイズζを決定する構成でなくてもよい。例えば、ノイズ決定部２８は、パーティクルが更新される全てのステップのうち、１つおきのステップにおいて（即ち、周期的に）パラメータノイズζを決定する構成であってもよい。
【００６６】
また、本明細書または図面に説明した技術要素は、単独であるいは各種の組合せによって技術的有用性を発揮するものであり、出願時請求項記載の組合せに限定されるものではない。また、本明細書または図面に例示した技術は複数目的を同時に達成するものであり、そのうちの一つの目的を達成すること自体で技術的有用性を持つものである。
【符号の説明】
【００６７】
１０：システムパラメータ同定装置、１２：入力値取得部、１４：観測値取得部、１６：第１関数記憶部、１８：第２関数記憶部、２０：状態ベクトル算出部、２２：推定値算出部、２４：パーティクルフィルタ実行部、２６：分散計算式記憶部、２８：ノイズ決定部、３０：外部記憶媒体、３２：入力値データ、３４：観測値データ、４０：制御器、４４：同定装置
【図１】
2017083922.tif 000028
【図２】
2017083922.tif 000029
【図３】
2017083922.tif 000030
【図４】
2017083922.tif 000031
【図５】
2017083922.tif 000032
【図６】
2017083922.tif 000033
【図７】
2017083922.tif 000034
【図８】
2017083922.tif 000035
【図９】
2017083922.tif 000036
【図１０】
2017083922.tif 000037
【図１１】
2017083922.tif 000038