(57)【要約】【課題】本発明の目的は、油圧システムの動特性に影響を与える物理量をパラメータとして推定するパラメータ推定方法の適用範囲を広げることである。【解決手段】油圧システムのパラメータ推定方法は、油圧システムの動特性に影響を与える物理量を推定対象パラメータとして推定するものである。油圧システムのパラメータ推定方法は、油圧システムの状態量と推定対象パラメータとを含む状態ベクトルを設定する工程と、状態ベクトルを用いて、油圧システムの状態方程式を設定する工程と、ベイズ推定に基づく状態推定手法を状態方程式に適用して、状態ベクトルを推定することにより、推定対象パラメータを推定する工程とを備える。【選択図】図１
2019057198.tif 000002
【特許請求の範囲】
【請求項１】
油圧システムの動特性に影響を与える物理量を推定対象パラメータとして推定する油圧システムのパラメータ推定方法であって、
前記油圧システムの状態量と前記推定対象パラメータとを含む状態ベクトルを設定する工程と、
前記状態ベクトルを用いて、前記油圧システムの状態方程式を設定する工程と、
ベイズ推定に基づく状態推定手法を前記状態方程式に適用して、前記状態ベクトルを推定することにより、前記推定対象パラメータを推定する工程とを備える、油圧システムのパラメータ推定方法。
【請求項２】
請求項１に記載の油圧システムのパラメータ推定方法であって、
前記推定対象パラメータは、前記油圧システムにおける作動油の体積弾性率である、油圧システムのパラメータ推定方法。
【請求項３】
請求項１又は２に記載の油圧システムのパラメータ推定方法であって、
前記状態方程式が非線形であり、
前記推定対象パラメータを推定する工程では、前記状態方程式に非線形の状態推定手法を適用して、前記状態ベクトルを推定することにより、前記推定対象パラメータを推定する、油圧システムのパラメータ推定方法。
【請求項４】
請求項１〜３の何れか１項に記載の油圧システムのパラメータ推定方法であって、
前記状態方程式が雑音項を含み、
前記推定対象パラメータを推定する工程では、前記状態方程式にカルマンフィルタを適用して、前記状態ベクトルを推定することにより、前記推定対象パラメータを推定する、油圧システムのパラメータ推定方法。
【請求項５】
油圧システムの動特性に影響を与える物理量を推定対象パラメータとして推定する油圧システムのパラメータ推定方法であって、
前記油圧システムの動特性を示す第１方程式を導出する工程と、
前記推定対象パラメータを含む状態ベクトルを用いた第２方程式を前記第１方程式から導出する工程と、
逐次推定手法を前記第２方程式に適用して、前記状態ベクトルを推定することにより、前記推定対象パラメータを推定する工程とを備える、油圧システムのパラメータ推定方法。
【請求項６】
請求項５に記載の油圧システムのパラメータ推定方法であって、
前記推定対象パラメータを推定する工程では、前記第２方程式に逐次最小二乗法を適用して、前記状態ベクトルを推定することにより、前記推定対象パラメータを推定する、油圧システムのパラメータ推定方法。
【請求項７】
請求項６に記載の油圧システムのパラメータ推定方法であって、さらに、
前記油圧システムの状態量を測定する工程を備え、
前記推定対象パラメータを推定する工程では、前記状態量の測定値に含まれるノイズに起因する前記推定対象パラメータの推定誤差を小さくするために、前記逐次最小二乗法を前記状態方程式に適用するに際して、前記状態量の測定値を補正する、油圧システムのパラメータ推定方法。
【請求項８】
請求項７に記載の油圧システムのパラメータ推定方法であって、
前記推定対象パラメータを推定する工程では、前記状態量の測定値の絶対値が大きいほど、その重みを大きくするように、前記状態量の測定値を補正する、油圧システムのパラメータ推定方法。
【請求項９】
請求項７に記載の油圧システムのパラメータ推定方法であって、
前記推定対象パラメータを推定する工程では、前記状態量の測定値の変化量が大きいほど、その重みを小さくするように、前記状態量の測定値を補正する、油圧システムのパラメータ推定方法。
【請求項１０】
請求項７に記載の油圧システムのパラメータ推定方法であって、
前記推定対象パラメータを推定する工程では、前記ノイズをフィルタによって除去することにより、前記状態量の測定値を補正する、油圧システムのパラメータ推定方法。
【請求項１１】
請求項５に記載の油圧システムのパラメータ推定方法であって、
前記第２方程式が雑音項を含み、
前記推定対象パラメータを推定する工程では、前記第２方程式にカルマンフィルタを適用して、前記状態ベクトルを推定することにより、前記推定対象パラメータを推定する、油圧システムのパラメータ推定方法。
【請求項１２】
請求項５〜１１の何れか１項に記載の油圧システムのパラメータ推定方法であって、
前記推定対象パラメータは、前記油圧システムにおける作動油の体積弾性率である、油圧システムのパラメータ推定方法。
【請求項１３】
請求項５〜１１の何れか１項に記載の油圧システムのパラメータ推定方法であって、
前記推定対象パラメータは、前記油圧システムにおける作動油の体積弾性率と縮流係数である、油圧システムのパラメータ推定方法。
【請求項１４】
請求項５〜１１の何れか１項に記載の油圧システムのパラメータ推定方法であって、
前記推定対象パラメータは、前記油圧システムにおける作動油の体積弾性率の圧力依存性を示すものである、油圧システムのパラメータ推定方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、油圧システムの動特性に影響を与える物理量をパラメータとして推定する油圧システムのパラメータ推定方法に関する。
【背景技術】
【０００２】
油圧システムの制御において、油圧システムの動特性を正確にモデル化することは重要である。油圧システムの動特性は、作動油の体積弾性率や粘度、レイノルズ数等の様々な物理量に大きく依存する。そして、これらの物理量には、油圧システムの稼働中に動的に変化するパラメータに依存するものがある。例えば、作動油の体積弾性率は、作動油への気泡混入率、作動油の温度や圧力等に依存する。そのため、油圧システムの動特性に影響を与える物理量については、その値をリアルタイムで取得することが望ましい。
【０００３】
例えば、作動油の体積弾性率をリアルタイムで取得する方法として、油圧システムの作動油の体積変化とそれに対する発生圧力を測定し、その比率（つまり、体積弾性率の定義そのもの）を計算することが考えられる。しかしながら、一般的な油圧システムにおいて、作動油の体積変化を正確に測定するためには、通常の作業を中断したり、作動油の体積変化とそれに対する発生圧力を測定するための回路を別途設ける必要がある。そのため、油圧システムの運用効率の低下や運用コストの増大を招いてしまう。
【０００４】
作動油の体積弾性率をリアルタイムで取得する他の方法として、体積弾性率を数学関数でモデル化し、その関数の変数である物理パラメータをセンサにて測定し、代入することが考えられる。しかしながら、作動油の体積弾性率が依存する作動油への気泡混入率は、センサによって直接測定することが難しく、そのため、全ての状態に対して十分な精度を保証する数学モデルを構築することは困難である。
【０００５】
そこで、作動油の体積弾性率をリアルタイムで取得する方法として、システム同定理論の適用が考えられる。システム同定理論を用いたパラメータ同定法として、例えば、最小二乗法がある。
【０００６】
下記特許文献１には、最小二乗法に基づくパラメータ同定法が開示されている。特許文献１に記載のパラメータ同定法では、鉛直多関節油圧マニピュレータの非線形モデルを表す状態空間方程式を定め、当該状態空間方程式に、従来の逐次同定法を修正した新たな逐次同定法を適用する。この逐次同定法では、各ステップにおいて、測定したパラメータ及び既知のパラメータを用いて、状態空間方程式に基づき、体積弾性率及び流量係数を含む未知のパラメータについての線形方程式を立てるが、当該線形方程式を用いた同定計算は行わない。各ステップの終了後に、未知のパラメータについての線形方程式を、全ステップを横断して２つの線形方程式に厳密に統合し、これらを各時刻に亘って連立することで２つの正規方程式を導く。そして、これら２つの正規方程式を解くことで、体積弾性率及び流量係数を含む未知のパラメータを同定する。
【先行技術文献】
【特許文献】
【０００７】
【特許文献１】特開２０１５−７７６４３号公報
【発明の概要】
【発明が解決しようとする課題】
【０００８】
しかしながら、特許文献１に記載のパラメータ同定法は、状態の全てが観測可能であるとともに状態方程式が同定パラメータに対して線形に変形できる場合にしか、適用できない。加えて、同定パラメータの数が増えると、計算コストが指数関数的に増加してしまい、リアルタイムでの推定が困難になる。したがって、パラメータ同定法の適用範囲が著しく限定される。
【０００９】
本発明の目的は、油圧システムの動特性に影響を与える物理量をパラメータとして推定するパラメータ推定方法の適用範囲を広げることである。
【課題を解決するための手段】
【００１０】
本願の発明者は、上記の目的を達成するために、推定対象の物理量を推定対象パラメータとして含む状態ベクトルを設定することに着目して検討を進めた。そして、当該状態ベクトルを用いた方程式に対して、ベイズ推定に基づく状態推定手法又は逐次推定手法を適用することにより、パラメータ推定手法の適用範囲を広げることができるという新たな知見を得るに至った。本発明は、このような知見に基づいて完成されたものである。
【００１１】
本発明の第１の態様による油圧システムのパラメータ推定方法は、油圧システムの動特性に影響を与える物理量を推定対象パラメータとして推定する油圧システムのパラメータ推定方法であって、前記油圧システムの状態量と前記推定対象パラメータとを含む状態ベクトルを設定する工程と、前記状態ベクトルを用いて、前記油圧システムの状態方程式を設定する工程と、ベイズ推定に基づく状態推定手法を前記状態方程式に適用して、前記状態ベクトルを推定することにより、前記推定対象パラメータを推定する工程とを備える。
【００１２】
上記第１の態様による油圧システムのパラメータ推定方法においては、ベイズ推定に基づく状態推定手法を適用するので、状態ベクトルに含まれる油圧システムの状態量の全てが観測可能でなくてもよい。その結果、パラメータ推定方法の適用範囲を広げることができる。
【００１３】
上記第１の態様による油圧システムのパラメータ推定方法において、好ましくは、前記推定対象パラメータは、前記油圧システムにおける作動油の体積弾性率である。この場合、油圧システムの動作を停止することなく、作動油の体積弾性率を得ることができる。
【００１４】
上記第１の態様による油圧システムのパラメータ推定方法において、前記状態方程式が非線形であってもよい。この場合、前記推定対象パラメータを推定する工程では、好ましくは、前記状態方程式に非線形の状態推定手法を適用して、前記状態ベクトルを推定することにより、前記推定対象パラメータを推定する。このような態様においては、状態方程式が非線形であっても、推定対象パラメータを推定することができるので、パラメータ推定方法の適用範囲を広げることができる。
【００１５】
上記第１の態様による油圧システムのパラメータ推定方法において、好ましくは、前記状態方程式が雑音項を含み、前記推定対象パラメータを推定する工程では、前記状態方程式にカルマンフィルタを適用して、前記状態ベクトルを推定することにより、前記推定対象パラメータを推定する。この場合、油圧システムが有する雑音等の影響を考慮して、推定対象パラメータを推定することができる。
【００１６】
本発明の第２の態様による油圧システムのパラメータ推定方法は、油圧システムの動特性に影響を与える物理量を推定対象パラメータとして推定する油圧システムのパラメータ推定方法であって、前記油圧システムの動特性を示す第１方程式を導出する工程と、前記推定対象パラメータを含む状態ベクトルを用いた第２方程式を前記第１方程式から導出する工程と、逐次推定手法を前記第２方程式に適用して、前記状態ベクトルを推定することにより、前記推定対象パラメータを推定する工程とを備える。
【００１７】
上記第２の態様による油圧システムのパラメータ推定方法においては、逐次推定手法を適用するので、推定対象パラメータを推定する際の計算コストの増加を抑制することができる。そのため、パラメータ推定手法の適用範囲を広げることができる。
【００１８】
上記第２の態様による油圧システムのパラメータ推定方法において、好ましくは、前記推定対象パラメータを推定する工程では、前記第２方程式に逐次最小二乗法を適用して、前記状態ベクトルを推定することにより、前記推定対象パラメータを推定する。この場合、あるタイムステップでの推定対象パラメータの推定値を、その前のタイムステップでの推定対象パラメータの推定値の関数として逐次的に（少ない計算コストで）導出することができる。
【００１９】
上記第２の態様による油圧システムのパラメータ推定方法は、好ましくは、さらに、前記油圧システムの状態量を測定する工程を備え、前記推定対象パラメータを推定する工程では、前記状態量の測定値に含まれるノイズに起因する前記推定対象パラメータの推定誤差を小さくするために、前記逐次最小二乗法を前記第２方程式に適用するに際して、前記状態量の測定値を補正する。この場合、補正した状態量の測定値を用いて、推定対象パラメータを推定することができるので、推定対象パラメータの推定精度を向上させることができる。
【００２０】
上記第２の態様による油圧システムのパラメータ推定方法において、前記推定対象パラメータを推定する工程では、前記状態量の測定値の絶対値が大きいほど、その重みを大きくするように、前記状態量の測定値を補正するようにしてもよい。この場合、状態量の測定値の絶対値の大きさに応じて、状態量の測定値が推定対象パラメータの推定誤差に与える影響を調整することができるため、推定対象パラメータの推定精度を向上させることができる。
【００２１】
上記第２の態様による油圧システムのパラメータ推定方法において、前記推定対象パラメータを推定する工程では、前記状態量の測定値の変化量が大きいほど、その重みを小さくするように、前記状態量の測定値を補正するようにしてもよい。この場合、状態量の測定値の変化量の大きさに応じて、状態量の測定値が推定対象パラメータの推定誤差に与える影響を調整することができるので、推定対象パラメータの推定精度を向上させることができる。
【００２２】
上記第２の態様による油圧システムのパラメータ推定方法において、前記推定対象パラメータを推定する工程では、前記ノイズをフィルタによって除去することにより、前記状態量の測定値を補正するようにしてもよい。この場合、状態量の測定値に含まれるノイズがフィルタによって除去されているので、推定対象パラメータの推定精度を向上させることができる。
【００２３】
上記第２の態様による油圧システムのパラメータ推定方法において、好ましくは、前記第２方程式が雑音項を含み、前記推定対象パラメータを推定する工程では、前記第２方程式にカルマンフィルタを適用して、前記状態ベクトルを推定することにより、前記推定対象パラメータを推定する。この場合、油圧システムが有する雑音等の影響を考慮して、推定対象パラメータを推定することができる。
【００２４】
上記第２の態様による油圧システムのパラメータ推定方法において、好ましくは、前記推定対象パラメータは、前記油圧システムにおける作動油の体積弾性率である。この場合、油圧システムの動作を停止しなくても、作動油の体積弾性率を得ることができる。
【００２５】
上記第２の態様による油圧システムのパラメータ推定方法において、好ましくは、前記推定対象パラメータは、前記油圧システムにおける作動油の体積弾性率と縮流係数である。この場合、作動油の体積弾性率だけでなく、縮流係数も推定することができる。
【００２６】
上記第２の態様による油圧システムのパラメータ推定方法において、好ましくは、前記推定対象パラメータは、前記油圧システムにおける作動油の体積弾性率の圧力依存性を示すものである。この場合、作動油の体積弾性率が作動油の圧力に依存する度合いを推定することができる。
【発明の効果】
【００２７】
本発明の油圧システムのパラメータ推定方法によれば、パラメータ推定手法の適用範囲を広げることができる。
【図面の簡単な説明】
【００２８】
【図１】本発明の第１の実施の形態による油圧システムのパラメータ推定方法を示すフローチャートである。
【図２】本発明の第１の実施の形態による油圧システムのパラメータ推定方法が適用される油圧システムを示す回路図である。
【図３】本発明の第１の実施の形態で用いるｆｈｉｎ（ｘｖ）、ｆｈｏｕｔ（ｘｖ）、ｆｒｉｎ（ｘｖ）及びｆｒｏｕｔ（ｘｖ）を示すグラフである。
【図４】油圧システムへの入力としての方向流量制御弁のスプール位置を示すグラフである。
【図５】本発明の第１の実施の形態での推定結果を示すグラフであって、作動油の体積弾性率Ｋの推定を開始してからの時間変化を示すグラフである。
【図６】本発明の第１の実施の形態での推定結果を示すグラフであって、圧力Ｐｐの推定を開始してからの時間変化を示すグラフである。
【図７】本発明の第１の実施の形態での推定結果を示すグラフであって、圧力Ｐｈの推定を開始してからの時間変化を示すグラフである。
【図８】本発明の第１の実施の形態での推定結果を示すグラフであって、圧力Ｐｒの推定を開始してからの時間変化を示すグラフである。
【図９】本発明の第１の実施の形態での推定結果を示すグラフであって、スプール位置ｘｃの推定を開始してからの時間変化を示すグラフである。
【図１０】本発明の第１の実施の形態での推定結果を示すグラフであって、スプール速度ｘｃの推定を開始してからの時間変化を示すグラフである。
【図１１】本発明の第２の実施の形態による油圧システムのパラメータ推定方法を示すフローチャートである。
【図１２】逐次最小二乗法によって作動油の体積弾性率Ｋを推定した結果を示すグラフである。
【図１３】逐次最小二乗法によって作動油の体積弾性率Ｋを推定した結果を示すグラフである。
【図１４】逐次最小二乗法によって縮流係数Ｃを推定した結果を示すグラフである。
【図１５】重み付き逐次最小二乗法によって作動油の体積弾性率Ｋを推定した結果を示すグラフである。
【図１６】重み付き逐次最小二乗法によって縮流係数Ｃを推定した結果を示すグラフである。
【図１７】重み付き逐次最小二乗法によって作動油の体積弾性率Ｋを推定した結果を示すグラフである。
【図１８】重み付き逐次最小二乗法によって縮流係数Ｃを推定した結果を示すグラフである。
【図１９】周波数フィルタと逐次最小二乗法によって作動油の体積弾性率Ｋを推定した結果を示すグラフである。
【図２０】周波数フィルタと逐次最小二乗法によって縮流係数Ｃを推定した結果を示すグラフである。
【図２１】大気圧における体積弾性率Ｋ０の推定結果を示すグラフである。
【図２２】比例係数Ｋｐの推定結果を示すグラフである。
【図２３】大気圧における体積弾性率Ｋ０の推定結果を示すグラフである。
【図２４】比例係数Ｋｐの推定結果を示すグラフである。
【発明を実施するための形態】
【００２９】
以下、添付図面を参照しながら、本発明の実施の形態について詳述する。
【００３０】
［第１の実施の形態］
図１を参照して、本発明の第１の実施の形態による油圧システムのパラメータ推定方法について説明する。図１は、第１の実施の形態による油圧システムのパラメータ推定方法を示すフローチャートである。
【００３１】
第１の実施の形態による油圧システムのパラメータ推定方法は、油圧システムの動特性に影響を与える物理量を推定対象パラメータとして推定するものである。第１の実施の形態による油圧システムのパラメータ推定方法は、油圧システムの動特性を示す方程式を設定する工程（ステップＳ１０）と、油圧システムの状態量と推定対象パラメータとを含む状態ベクトルを設定する工程（ステップＳ１１）と、状態ベクトルを用いて、油圧システムの状態方程式と観測方程式を設定する工程（ステップＳ１２）と、ベイズ推定に基づく状態推定手法を状態方程式に適用して、状態ベクトルを推定することにより、推定対象パラメータを推定する工程（ステップＳ１３）とを備える。
【００３２】
このようなパラメータ推定手法は、例えば、制御ＣＰＵやシミュレータによって実行される。制御ＣＰＵは、例えば、実機システムを制御する。
【００３３】
第１の実施の形態による油圧システムのパラメータ推定方法は、例えば、図２に示す油圧システム１０に適用して、油圧システム１０の作動油の体積弾性率Ｋを推定対象パラメータとして推定するために用いられる。なお、第１の実施の形態による油圧システムのパラメータ推定方法が適用される油圧システムは、図２に示す油圧システム１０に限定されるものではない。また、第１の実施の形態による油圧システムのパラメータ推定方法によって推定される推定対象パラメータは、油圧システム１０の作動油の体積弾性率Ｋに限定されない。
【００３４】
図２を参照して、第１の実施の形態による油圧システムのパラメータ推定方法が適用される油圧システム１０について説明する。図２は、油圧システム１０を示す回路図である。
【００３５】
油圧システム１０は、スプール式の方向流量制御弁１４のスプール変位ｘｖを制御入力とし、油圧シリンダ１６の動作を制御するものである。油圧システム１０は、油圧ポンプ１２と、方向流量制御弁１４と、油圧シリンダ１６と、作動油タンク１８と、リリーフ弁２０とを備える。
【００３６】
油圧ポンプ１２は、作動油を吐出する。方向流量制御弁１４は、油圧ポンプ１２から油圧シリンダ１６に供給される作動油の方向及び流量を変化させる。油圧シリンダ１６は、油圧ポンプ１２が吐出する作動油（つまり、油圧ポンプ１２から供給される作動油）によって駆動される。作動油タンク１８は、油圧ポンプ１２から吐出される作動油を貯留する。リリーフ弁２０は、油圧ポンプ１２から吐出される作動油の圧力が所定の圧力を超えないように、開閉動作を行う。油圧ポンプ１２と方向流量制御弁１４とは、配管２２によって接続されている。方向流量制御弁１４と油圧シリンダ１６のシリンダヘッド側の液室とは、配管２４によって接続されている。油圧シリンダ１６のロッド側の液室と方向流量制御弁１４とは、配管２６によって接続されている。
【００３７】
再び、図１を参照しながら、第１の実施の形態による油圧システムのパラメータ推定方法の各工程について説明する。
【００３８】
先ず、ステップＳ１０において、油圧システム１０の動特性を示す方程式を設定する。油圧システム１０の動特性は、以下に示す方程式によって表される。
【００３９】
【数１】
2019057198.tif 000003
【００４０】
ここで、ｍは油圧シリンダ１６及び負荷対象の質量を示し、ｘｃはシリンダ位置（油圧シリンダ１６内でのシリンダヘッドの位置）を示し、Ａｈは油圧シリンダ１６のシリンダヘッド側の液室１６１の断面積を示し、Ａｒは油圧シリンダ１６のロッド側の液室１６２の断面積を示し、Ｐｐは配管２２内の圧力を示し、Ｐｈは配管２４内の圧力を示し、Ｐｒは配管２６内の圧力を示し、Ｐｔは作動油タンク１８内の圧力を示し、Ｋは作動油の体積弾性率を示し、Ｖｐは配管２２の容積を示し、Ｖｈは配管２４の容積を示し、Ｖｒは配管２６の容積を示し、ｌは油圧シリンダ１６の全長を示す。Ｑｐは油圧ポンプ１２から吐出される油量を示し、Ｑｈｉｎは方向流量制御弁１４を通過して油圧シリンダ１６の液室１６１へ流入する油量を示し、Ｑｈｏｕｔは油圧シリンダ１６の液室１６１から排出されて方向流量制御弁１４を通過する油量を示し、Ｑｒｉｎは方向流量制御弁１４を通過して油圧シリンダ１６の液室１６２へ流入する油量を示し、Ｑｒｏｕｔは油圧シリンダ１６の液室１６２から排出されて方向流量制御弁１４を通過する油量を示し、ＱｐＲはリリーフ弁２０を通過する油量を示す。
【００４１】
Ｑｐ、Ｑｈｉｎ、Ｑｈｏｕｔ、Ｑｒｉｎ、Ｑｒｏｕｔ、ＱｐＲは、以下に示す方程式によって表される。
【００４２】
【数２】
2019057198.tif 000004
【００４３】
ここで、ｆｈｉｎ（ｘｖ）はＡｈｉｎとｘｖとの関係を示す任意の関数であり、ｆｈｏｕｔ（ｘｖ）はＡｈｏｕｔとｘｖとの関係を示す任意の関数であり、ｆｒｉｎ（ｘｖ）はＡｒｉｎとｘｖとの関係を示す任意の関数であり、ｆｒｏｕｔ（ｘｖ）はＡｒｏｕｔとｘｖとの関係を示す任意の関数であり、ｆｐＲ（Ｐｐ−Ｐｔ）は、ＡｐＲとＰｐ−Ｐｔとの関係を示す任意の関数である。また、Ｃは縮流係数を示し、ρは作動油の密度を示し、Ａｈｉｎは方向流量制御弁１４において作動油が液室１６１に向かって流れる流路での開口面積を示し、Ａｈｏｕｔは方向流量制御弁１４において液室１６１からの作動油が流れる流路での開口面積を示し、Ａｒｉｎは方向流量制御弁１４において作動油が液室１６２に向かって流れる流路での開口面積を示し、Ａｒｏｕｔは方向流量制御弁１４において液室１６２からの作動油が流れる流路での開口面積を示し、ＡｐＲはリリーフ弁２０の開口面積を示し、ｘｖは方向流量制御弁１４のスプール位置を示す。
【００４４】
油圧システム１０の動特性を示す方程式を設定した後、ステップＳ１１において、油圧システム１０の状態量と推定対象パラメータとを含む状態ベクトルを設定する。
【００４５】
油圧システム１０の状態量は、圧力Ｐｐ、圧力Ｐｈ、圧力Ｐｒ、スプール位置ｘｃ及びスプール速度ｘｃである。また、前述のとおり、推定対象パラメータは、体積弾性率Ｋである。油圧システム１０の状態ベクトルｘは、以下のように表される。
【００４６】
【数３】
2019057198.tif 000005
【００４７】
油圧システム１０の状態ベクトルｘを設定した後、ステップＳ１２において、状態ベクトルｘを用いて、油圧システム１０の状態方程式と観測方程式を設定する。
【００４８】
油圧システム１０の状態方程式は、非線形である。油圧システム１０の状態方程式は、以下のように表される。
【００４９】
【数４】
2019057198.tif 000006
【００５０】
つまり、油圧システム１０の状態方程式は、１階の微分方程式の形で表される。
【００５１】
なお、上記状態方程式におけるｖ∈Ｒ６×１は、平均ベクトルが０∈Ｒ６×１であり、分散共分散行列がＱ∈Ｒ６×６である白色雑音ベクトルを表す。つまり、上記状態方程式は、雑音項を含む。
【００５２】
ここで、圧力センサや位置センサにより、油圧システム１０の状態量のうち、圧力Ｐｐ、圧力Ｐｈ、圧力Ｐｒ及びスプール位置ｘｃが観測可能であるとする。このとき、油圧システム１０の観測方程式は、以下のように表される。
【００５３】
【数５】
2019057198.tif 000007
【００５４】
なお、上記観測方程式におけるｗ∈Ｒ４×４は、平均ベクトルが０∈Ｒ４×４であり、分散共分散行列がＲ∈Ｒ４×４である白色雑音ベクトルを表す。
【００５５】
油圧システム１０の状態方程式と観測方程式を設定した後、ステップＳ１３において、ベイズ推定に基づく状態推定手法としてのUnscentedカルマンフィルタ（つまり、非線形の状態推定手法）を状態方程式に適用して、状態ベクトルｘを推定することにより、推定対象パラメータとしての作動油の体積弾性率Ｋを推定する。
【００５６】
［第１の実施の形態での推定に用いる条件］
以下、本実施の形態での推定に用いる条件について説明する。また、当該条件を用いて推定した結果についても説明する。
【００５７】
本実施の形態で用いる物理定数を、以下の表１に示す。
【００５８】
【表１】
2019057198.tif 000008
【００５９】
また、本実施の形態で用いるｆｈｉｎ（ｘｖ）、ｆｈｏｕｔ（ｘｖ）、ｆｒｉｎ（ｘｖ）及びｆｒｏｕｔ（ｘｖ）を、図３に示す。ｆｈｉｎ（ｘｖ）とｆｒｏｕｔ（ｘｖ）は、それぞれ、スプール位置ｘｖが負である場合は、開口面積がゼロであり、スプール位置ｘｖが正である場合は、スプール位置ｘｖが大きくなるに従って、開口面積が大きくなる。ｆｒｉｎ（ｘｖ）とｆｈｏｕｔ（ｘｖ）は、それぞれ、スプール位置ｘｖが負である場合は、スプール位置ｘｖの絶対値が小さくなるに従って、開口面積が小さくなり、スプール位置ｘｖが正である場合は、開口面積がゼロである。油圧システム１０への入力として、方向流量制御弁１４のスプール位置ｘｖを、図４に示すように与える。
【００６０】
作動油の体積弾性率Ｋの時間変化ｆＫは、ゼロであってもよいし、正の値（例えば、５０）であってもよい。作動油の体積弾性率Ｋの時間変化ｆＫがゼロである場合、作動油の体積弾性率Ｋが定数であると仮定した推定が行われる。一方、作動油の体積弾性率Ｋの時間変化ｆＫが正の値である場合、作動油の体積弾性率Ｋは変動パラメータとなり、モデル化誤差を模擬することができる。なお、本実施の形態では、作動油の体積弾性率Ｋの時間変化ｆＫを５０にしている。
【００６１】
このような条件の下で、作動油の体積弾性率Ｋを推定し、図５〜図１０に示すような結果を得た。
【００６２】
図５に示すように、推定を開始してから０．１秒以内に、作動油の体積弾性率Ｋの推定値が、初期値（５００ＭＰａ）から真値（１５００ＭＰａ）付近に移動していることが確認できる。作動油の体積弾性率Ｋ及びシリンダ速度が直接計測されていないにも関わらず、これらの数値を推定できていることが判る。
【００６３】
また、図５に示すように、推定モデルでは考慮されていない真値の時間変化に対して、推定値が追従できていることも確認できる。これにより、状態方程式にモデル化誤差が含まれていても、雑音ベクトルvの分散共分散行列Qを適切に設定することで実用的な推定が実現されていることが判る。
【００６４】
また、本実施の形態では、推定値の標準偏差を指標として、推定の信頼度を評価することができる。図５では、油流の状態に応じて推定値の標準偏差が時々刻々変化しており、推定誤差を考慮に入れた制御器を設計する際の定量的な指標として利用できる。
【００６５】
また、図６〜図１０に示すように、作動油の体積弾性率Ｋの推定のみならず、雑音ベクトルwの分散共分散行列Rを適切に設定することにより、雑音の入った測定値から状態ベクトルxの真値も推定できていることが確認できる。
【００６６】
［第１の実施の形態の応用例］
第１の実施の形態では、ベイズ推定に基づいた状態推定手法として、Unscentedカルマンフィルタを用いていたが、例えば、線形カルマンフィルタや拡張カルマンフィルタ、粒子フィルタ、アンサンブルカルマンフィルタ等を用いてもよい。なお、線形カルマンフィルタを用いる場合には、状態方程式は線形でなければならないが、例えば、動作範囲が限られた油圧システムであれば、その範囲内で非線形な状態方程式を線形な状態方程式に近似することにより、線形カルマンフィルタを適用することができる。また、限定的な仮定の下では、油圧システムの状態方程式を線形に表すこともできる。例えば、作動油の流れを層流であると仮定したり、油圧シリンダの代わりに、油圧モータを採用するのであれば、油圧システムの状態方程式を線形に表すことができる。
【００６７】
推定対象パラメータが動的に変化する場合、当該推定対象パラメータに関する既知の特性を状態方程式に反映させてもよい。例えば、作動油の体積弾性率Ｋは、等温条件において、圧力に比例するので、これを状態方程式に反映させてもよい。
【００６８】
状態方程式を第１原理モデリングではなく、ブラックボックスモデリングで導出し、その係数を推定するようにしてもよい。
【００６９】
第１の実施の形態による油圧システムのパラメータ推定方法によって得られる推定値は、油圧システムの制御に利用してもよいし、油圧システムのオペレータへ提示することにより、オペレータによる油圧システムの操作の支援及び／又はガイダンスを行うようにしてもよい。
【００７０】
［第２の実施の形態］
図１１を参照しながら、本発明の第２の実施の形態による油圧システムのパラメータ推定方法について説明する。図１１は、第２の実施の形態による油圧システムのパラメータ推定方法を示すフローチャートである。
【００７１】
第２の実施の形態による油圧システムのパラメータ推定方法は、第１の実施の形態による油圧システムのパラメータ推定方法と比べて、ベイズ推定に基づく状態推定手法の代わりに、逐次推定手法を用いている点で異なる。また、状態ベクトルが推定対象パラメータのみを含む場合もあり得る点で異なる。
【００７２】
第２の実施の形態による油圧システムのパラメータ推定方法は、油圧システムの動特性を示す第１方程式を設定する工程（ステップＳ２０）と、推定対象パラメータを含む状態ベクトルを用いた第２方程式を第１方程式から導出する工程（ステップＳ２１）と、逐次推定手法を第２方程式に適用して、状態ベクトルを推定することにより、推定対象パラメータを推定する工程（ステップＳ２２）とを備える。
【００７３】
このようなパラメータ推定方法は、例えば、制御ＣＰＵやシミュレータによって実行される。制御ＣＰＵは、例えば、実機システムを制御する。
【００７４】
ここで、ステップＳ２０の工程は、第１の実施の形態におけるステップＳ１０の工程と同じであるから、その詳細な説明は省略する。ステップＳ２１及びステップＳ２２の工程の詳細については、後述する。
【００７５】
第２の実施の形態による油圧システムのパラメータ推定方法は、第１の実施の形態による油圧システムのパラメータ推定方法と同様に、例えば、図２に示す油圧システム１０に適用して、油圧システム１０の作動油の体積弾性率Ｋを推定するために用いられる。なお、第２の実施の形態による油圧システムのパラメータ推定方法が適用される油圧システムは、図２に示す油圧システム１０に限定されるものではない。また、第２の実施の形態による油圧システムのパラメータ推定方法によって推定される推定対象パラメータは、油圧システム１０の作動油の体積弾性率Ｋに限定されない。
【００７６】
逐次推定手法は、各タイムステップでの推定値がその前のタイムステップでの推定値の関数として導出できるものであれば、特に限定されない。逐次推定手法は、逐次最小二乗法やカルマンフィルタを含む。逐次推定手法を適用することにより、計算コストを少なくすることができる。
【００７７】
【表２】
2019057198.tif 000009
【００７８】
表２は、本実施の形態にて採用される逐次推定手法と推定対象パラメータとの組み合わせを示す。例えば、（１）の組み合わせに係る態様は、逐次最小二乗法によって作動油の体積弾性率Ｋを推定する場合を示す。
【００７９】
以下、表２に示す（１）〜（８）の組み合わせに係る態様について説明する。なお、表２における空欄は、他の組み合わせに係る態様にて説明する逐次推定手法を適用することで実施できるため、その詳細な説明を省略していることを示しているだけであり、実施が不可能であることを意味しているのではない。例えば、（５）の組み合わせに係る態様にて説明する周波数フィルタ付の逐次最小二乗法を適用しても、作動油の体積弾性率Ｋを推定することができる。
【００８０】
また、以下では、表２に示す（１）〜（８）の組み合わせに係る態様の推定結果についても説明する。なお、推定に用いた条件は、第１の実施の形態での条件、つまり、表１、図３及び図４に示すものと同じであるから、その詳細な説明は省略する。
【００８１】
（１）の組み合わせに係る態様について
先ず、逐次最小二乗法によって作動油の体積弾性率Ｋを推定する場合について説明する。
【００８２】
この場合、ステップＳ２１において、油圧システム１０の動特性を表す方程式のうち、Ｐｐ、Ｐｈ及びＰｒを微分したものについての方程式を推定対象パラメータとしての作動油の体積弾性率Ｋに対して線形な形で離散化する。
【００８３】
先ず、Ｐｐ、Ｐｈ及びＰｒを微分したものについての方程式の右辺を、Ｋとそれ以外の項ａｐ、ａｈ、ａｒの積として陽に表すと、以下のようになる。
【００８４】
【数６】
2019057198.tif 000010
【００８５】
これらの式の左辺、つまり、Ｐｐ、Ｐｈ及びＰｒを微分したものを時間区切りΔｔで差分近似すると、以下のようになる。
【００８６】
【数７】
2019057198.tif 000011
【００８７】
ここで、ｋは時間区切りΔｔで離散化されたタイムステップを表す。Ｋ以外のパラメータは、観測可能であるとする。
【００８８】
各タイムステップｎの時点での上記式の両辺の二乗誤差の総和は、以下のようになる。
【００８９】
【数８】
2019057198.tif 000012
【００９０】
このような二乗誤差の総和を最小化するＫ（ｎ）を逐次最小二乗法によって導出する。これにより、ステップＳ２２が実行される。なお、上記式において、計算開始タイムステップａは、油圧システム１０の特性に応じて調整するパラメータである。
【００９１】
図１２は、作動油の体積弾性率Ｋの推定結果を示すグラフである。図１２に示すように、時間が経過するに従って、推定値が真値に近づいており、推定の精度が向上していることが判る。なお、推定値が真値に一致していないのは、Ｐｐ、Ｐｈ及びＰｒを微分したものを１次のオイラー法で差分近似したことで生じる近似誤差が原因である。
【００９２】
（２）の組み合わせに係る態様について
続いて、逐次最小二乗法によって作動油の体積弾性率Ｋと縮流係数Ｃを推定する場合について説明する。
【００９３】
この場合、ステップＳ２１において、油圧システム１０の動特性を表す方程式のうち、Ｐｐ、Ｐｈ及びＰｒを微分したものについての方程式を推定対象パラメータとしてのＫ及びＫＣに対して線形な形で離散化する。
【００９４】
先ず、Ｐｐ、Ｐｈ及びＰｒを微分したものについての方程式の右辺のうち、Ｋの項とＫＣの項とを陽に表すと、以下のようになる。
【００９５】
【数９】
2019057198.tif 000013
【００９６】
これらの式の左辺、つまり、Ｐｐ、Ｐｈ及びＰｒを微分したものを時間区切りΔｔで差分近似すると、以下のようになる。
【００９７】
【数１０】
2019057198.tif 000014
【００９８】
ここで、ｋは時間区切りΔｔで離散化されたタイムステップを表す。
【００９９】
左辺のベクトルをｙＬＳ（ｋ）、右辺のベクトルをｘＬＳ、行列をＡ（ｋ）とすると、上記式は、以下のようになる。
【０１００】
【数１１】
2019057198.tif 000015
【０１０１】
ここで、ｘＬＳ以外のパラメータは観測可能であるとする。
【０１０２】
各タイムステップｎの時点での上記式の両辺の二乗誤差の総和は、以下のようになる。
【０１０３】
【数１２】
2019057198.tif 000016
【０１０４】
このような二乗誤差の総和を最小化するｘＬＳ（ｎ）を逐次最小二乗法によって導出する。その後、ＫＣをＫで除することにより、縮流係数Ｃを推定する。これにより、ステップＳ２２が実行される。なお、上記式において、計算開始タイムステップａは、油圧システム１０の特性に応じて調整するパラメータである。
【０１０５】
図１３は、作動油の体積弾性率Ｋの推定結果を示すグラフである。図１４は、縮流係数Ｃの推定結果を示すグラフである。なお、図１３及び図１４では、計測値Ｐｐ、Ｐｈ、Ｐｒにノイズがない場合の推定結果だけでなく、計測値Ｐｐ、Ｐｈ、Ｐｒにノイズがある場合の推定結果も示している。ノイズは、計測値Ｐｐ、Ｐｈ、Ｐｒに対して平均がゼロで標準偏差が０．０４［ＭＰａ］であるホワイトノイズが加わり、且つ、計測値ｘｃに対して平均がゼロで標準偏差が１０［ｍｍ］のホワイトノイズが加わったものである。
【０１０６】
図１３に示すように、作動油の体積弾性率Ｋは、（１）の場合と同様に、時間が経過するに従って、推定値が真値に近づいており、推定の精度が向上していることが判る。縮流係数Ｃについても、図１４に示すように、正しく推定されていることが確認できる。
【０１０７】
これに対して、計測値にノイズがある場合は、Ｋ及びＣの推定結果の精度が低下している。これは、計測値Ｐｐ、Ｐｈ、Ｐｒに含まれるノイズがＱｐ、Ｑｈｉｎ、Ｑｈｏｕｔ、Ｑｒｉｎ、Ｑｒｏｕｔ、ＱｐＲについての方程式中の非線形関数（平方根）を通ることで、流量の平均値にずれを生み出すことが主な原因である。推定精度の低下は、油圧システム１０の動特性を示す方程式の非線形性が大きくなるほど、また、計測データの真値に対するノイズの大きさが大きくなるほど、顕著になる。計測値にノイズがある油圧システムに対しては、単純に逐次最小二乗法を適用することが好ましくないことが判る。
【０１０８】
（３）の組み合わせに係る態様について
続いて、計測値に含まれるノイズによる推定精度の低下を軽減するために、重み付き逐次最小二乗法によって作動油の体積弾性率Ｋと縮流係数Ｃを推定する場合について説明する。なお、この態様では、ステップＳ２１の工程は、（２）の組み合わせに係る態様と同じである。
【０１０９】
この態様では、以下に示す二乗誤差の重み付き総和を最小化するｘＬＳ（ｎ）を導出する。その後、ＫＣをＫで除することにより、縮流係数Ｃを推定する。これにより、ステップＳ２２が実行される。
【０１１０】
【数１３】
2019057198.tif 000017
【０１１１】
ここで、Ｗ（ｋ）はタイムステップｋにおける重み行列であり、計測値の信頼度等、既知の性質に応じて設定される。計測値の真値に対するノイズの大きさが大きくなるほど、推定の精度が顕著に低下することから、計測値の絶対値が大きいほど、重みを大きくして、その計測値に対する信頼度を上げるような推定を行う。Ｗ（ｋ）は、例えば、以下のように設定される。
【０１１２】
【数１４】
2019057198.tif 000018
【０１１３】
ここで、αｐ、αｈ、αｒは、重みの大きさを調整するチューニングパラメータである。なお、Ｗ（ｋ）の設定はこれに限定されない。
【０１１４】
図１５は、作動油の体積弾性率Ｋの推定結果を示すグラフである。図１６は、縮流係数Ｃの推定結果を示すグラフである。図１３及び図１４と比べて、計測値にノイズが含まれているにも関わらず、推定精度が向上していることが判る。
【０１１５】
（４）の組み合わせに係る態様について
続いて、計測値に含まれるノイズによる推定精度の低下を軽減するために、重み付き逐次最小二乗法によって作動油の体積弾性率Ｋと縮流係数Ｃを推定する場合について説明する。（４）の組み合わせに係る態様では、（３）の組み合わせに係る態様と比べて、重み行列Ｗ（ｋ）が異なっている。この重み行列Ｗ（ｋ）は、１タイムステップにおける圧力の変化が大きいほど、重みが小さくなるように、設定されている。つまり、圧力の変化が油圧システム１０の動作として想定される値よりも大きい場合、その計測値をノイズとみなし、推定への寄与を低減させる。重み行列Ｗ（ｋ）は、例えば、以下のように設定される。
【０１１６】
【数１５】
2019057198.tif 000019
【０１１７】
ここで、αは、想定変化量を調整するチューニングパラメータである。なお、Ｗ（ｋ）の設定はこれに限定されない。
【０１１８】
図１７は、作動油の体積弾性率Ｋの推定結果を示すグラフである。図１８は、縮流係数Ｃの推定結果を示すグラフである。図１３及び図１４と比べて、計測値にノイズが含まれているにも関わらず、推定精度が向上していることが判る。
【０１１９】
（５）の組み合わせに係る態様について
続いて、計測値に含まれるノイズによる推定精度の低下を軽減するために、周波数フィルタと逐次最小二乗法によって作動油の体積弾性率Ｋと縮流係数Ｃを推定する場合について説明する。
【０１２０】
この場合、圧力計測値の周波数成分のうち、油圧システム１０の動作として想定される周波数成分以外をノイズとみなしてフィルタにより除去し、フィルタ処理後の計測値に対して逐次最小二乗法を適用する。使用する周波数フィルタは、以下のような移動平均フィルタである。
【０１２１】
【数１６】
2019057198.tif 000020
【０１２２】
ここで、Ｐｆｉｌｔｅｒｅｄはフィルタ適用後の圧力値であり、Ｐは圧力計測値（Ｐｐ、Ｐｈ、Ｐｒ）である。
【０１２３】
なお、使用する周波数フィルタは、移動平均フィルタに限定されず、例えば、線形フィルタ、バターワースフィルタ、チェビシェフフィルタ等であってもよい。
【０１２４】
図１９は、作動油の体積弾性率Ｋの推定結果を示すグラフである。図２０は、縮流係数Ｃの推定結果を示すグラフである。図１９及び図２０では、上記式においてＮ＝２１とした移動平均フィルタを用いた推定結果を示している。図１３及び図１４と比べて、計測値にノイズが含まれているにも関わらず、推定精度が向上していることが判る。
【０１２５】
（６）の組み合わせに係る態様について
続いて、（６）の組み合わせに係る態様について説明する。この態様は、計測値に含まれるノイズによる推定精度の低下を軽減するために、Unscentedカルマンフィルタによって作動油の体積弾性率Ｋを推定する場合である。なお、Unscentedカルマンフィルタによって作動油の体積弾性率Ｋを推定する手法は、第１の実施の形態で説明したものと同様であるから、その詳細な説明は省略する。
【０１２６】
この手法によれば、油圧システム１０の状態量と推定対象パラメータ（作動油の体積弾性率Ｋ）を同時に推定することができるので、油圧システム１０の状態量の全てを計測する必要がなくなる。また、ノイズの大きさ等、既知の特性を考慮した推定を行うので、推定精度に優れる。
【０１２７】
（７）の組み合わせに係る態様について
続いて、（７）の組み合わせに係る態様について説明する。この態様は、逐次最小二乗法に基づく作動油の体積弾性率Ｋの圧力依存性を推定するものである。
【０１２８】
作動油の体積弾性率Ｋは、作動油の圧力に応じて変化することが知られている。例えば、等温且つ低圧領域では、作動油の体積弾性率Ｋは、以下のように、圧力に比例する関数として近似することができる。
【０１２９】
【数１７】
2019057198.tif 000021
【０１３０】
ここで、Ｋ０は大気圧における体積弾性率であり、Ｋｐは比例係数であり、Ｐはゲージ圧（圧力計測値）である。Ｋ０及びＫｐは作動油の温度や作動油の気泡含有率等、動的に変化するパラメータに依存する。したがって、逐次推定手法によって推定する。
【０１３１】
ステップＳ２１において、油圧システム１０の動特性を表す方程式のうち、Ｐｐ、Ｐｈ及びＰｒを微分したものについての方程式を推定対象パラメータとしてのＫ０及びＫｐに対して線形な形で離散化する。
【０１３２】
先ず、油圧システム１０の動特性を表す方程式のうち、Ｐｐ、Ｐｈ及びＰｒを微分したものについての方程式に対して、上記式を代入し、Ｋ０の項とＫｐの項を陽に示すと、以下のようになる。
【０１３３】
【数１８】
2019057198.tif 000022
【０１３４】
左辺、つまり、Ｐｐ、Ｐｈ及びＰｒを微分したものを時間区切りΔｔで差分近似すると、以下のようになる。
【０１３５】
【数１９】
2019057198.tif 000023
【０１３６】
ここで、ｋは時間区切りΔｔで離散化されたタイムステップを表す。
【０１３７】
左辺のベクトルをｙＬＳ（ｋ）、右辺のベクトルをｘＬＳ、行列をＡ（ｋ）とすると、上記式は、以下のようになる。
【０１３８】
【数２０】
2019057198.tif 000024
【０１３９】
ここで、ｘＬＳ以外のパラメータは観測可能であるとする。
【０１４０】
各タイムステップｎの時点での上記式の両辺の二乗誤差の総和は、以下のようになる。
【０１４１】
【数２１】
2019057198.tif 000025
【０１４２】
このような二乗誤差の総和を最小化するｘＬＳ（ｎ）を逐次最小二乗法によって導出する。これにより、ステップＳ２２が実行される。なお、上記式において、計算開始タイムステップａは、油圧システム１０の特性に応じて調整するパラメータである。
【０１４３】
図２１は、大気圧における体積弾性率Ｋ０の推定結果を示すグラフである。図２２は、比例係数Ｋｐの推定結果を示すグラフである。なお、図２１及び図２２では、計測値Ｐｐ、Ｐｈ、Ｐｒにノイズがない場合の推定結果だけでなく、計測値Ｐｐ、Ｐｈ、Ｐｒにノイズがある場合の推定結果も示している。ノイズは、計測値Ｐｐ、Ｐｈ、Ｐｒに対して平均がゼロで標準偏差が０．０２［ＭＰａ］であるホワイトノイズが加わり、且つ、計測値ｘｃに対して平均がゼロで標準偏差が１０［ｍｍ］のホワイトノイズが加わったものである。
【０１４４】
図２１及び図２２に示すように、ノイズがない場合は、十分な推定精度を確保できていることが判る。
【０１４５】
これに対して、計測値にノイズがある場合は、Ｋ０及びＫｐの推定結果の精度が低下しており、単純に逐次最小二乗法を適用することが好ましくないことが判る。
【０１４６】
（８）の組み合わせに係る態様について
続いて、（８）の組み合わせに係る態様について説明する。この態様は、（７）と同様に、作動油の体積弾性率Ｋの圧力依存性を推定するものであるが、逐次最小二乗法ではなく、Unscentedカルマンフィルタを適用する。
【０１４７】
先ず、ステップＳ２１において、油圧システム１０の状態方程式及び観測方程式を設定する。具体的には、以下のとおりである。
【０１４８】
油圧システム１０の状態量（圧力Ｐｐ、圧力Ｐｈ、圧力Ｐｒ、シリンダ位置ｘｃ、シリンダ速度ｘｃに推定対象パラメータであるＫ０及びＫｐを追加したものを、油圧システム１０の状態ベクトルとして、以下のように設定する。
【０１４９】
【数２２】
2019057198.tif 000026
【０１５０】
このとき、油圧システム１０の状態方程式は、以下のように表される。
【０１５１】
【数２３】
2019057198.tif 000027
【０１５２】
つまり、油圧システム１０の状態方程式は、１階の微分方程式の形で表される。
【０１５３】
なお、上記状態方程式におけるｖ∈Ｒ７×１は、平均ベクトルが０∈Ｒ７×１であり、分散共分散行列がＱ∈Ｒ７×７である白色雑音ベクトルを表す。
【０１５４】
ここで、Ｋ０の時間変化ｆＫ０及びＫｐの時間変化ｆＫｐは、ゼロとする。
【０１５５】
また、圧力センサや位置センサにより、油圧システム１０の状態量のうち、圧力Ｐｐ、圧力Ｐｈ、圧力Ｐｒ及びスプール位置ｘｃが観測可能であるとする。このとき、油圧システム１０の観測方程式は、以下のように表される。
【０１５６】
【数２４】
2019057198.tif 000028
【０１５７】
なお、上記観測方程式におけるｗ∈Ｒ４×４は、平均ベクトルが０∈Ｒ４×４であり、分散共分散行列がＲ∈Ｒ４×４である白色雑音ベクトルを表す。
【０１５８】
油圧システム１０の状態方程式と観測方程式を設定した後、ステップＳ２２において、Unscentedカルマンフィルタを状態方程式に適用して、状態ベクトルｘを推定することにより、油圧システム１０の状態量及び推定対象パラメータを推定する。
【０１５９】
図２３は、大気圧における体積弾性率Ｋ０の推定結果を示すグラフである。図２４は、比例係数Ｋｐの推定結果を示すグラフである。図２３及び図２４では、計測値Ｐｐ、Ｐｈ、Ｐｒにノイズがない場合の推定結果だけでなく、計測値Ｐｐ、Ｐｈ、Ｐｒにノイズがある場合の推定結果も示している。ノイズは、計測値Ｐｐ、Ｐｈ、Ｐｒに対して平均がゼロで標準偏差が０．０２［ＭＰａ］であるホワイトノイズが加わり、且つ、計測値ｘｃに対して平均がゼロで標準偏差が１０［ｍｍ］のホワイトノイズが加わったものである。
【０１６０】
図２３及び図２４では、図２１及び図２２と比べて、推定対象パラメータＫ０、Ｋｐ及びシリンダ速度ｘｃが直接計測されていないにも関わらず、推定精度が向上していることが判る。
【０１６１】
［第２の実施の形態の応用例］
第２の実施の形態では、過去の計測値への依存性を調整するパラメータとして、計算開始ステップａを用いたが、忘却係数を用いて指数関数的に過去の計測値の重みを減少させてもよい。
【０１６２】
第２の実施の形態では、測定値の絶対値に応じた重みを設定したが、測定値の絶対値が所定の閾値よりも小さい場合に、その測定値を推定に用いないようにしてもよい。
【０１６３】
第２の実施の形態では、作動油の体積弾性率を圧力に関して線形な関数で近似したが、圧力に関して非線形な関数で近似してもよい。
【０１６４】
第２の実施の形態では、作動油の体積弾性率を圧力の関数としてモデル化したが、任意のパラメータの関数としてモデル化し、その係数を求めるようにしてもよい。例えば、作動油の体積弾性率を、温度の関数として、或いは、温度及び圧力の関数として、モデル化してもよい。
【０１６５】
第２の実施の形態では、Unscentedカルマンフィルタを用いたが、例えば、線形カルマンフィルタや拡張カルマンフィルタ、粒子フィルタ等を用いてもよい。なお、線形カルマンフィルタを用いる場合には、状態方程式は線形でなければならないが、例えば、動作範囲が限られた油圧システムであれば、その範囲内で非線形な状態方程式を線形な状態方程式に近似することにより、線形カルマンフィルタを適用することができる。また、限定的な仮定の下では、油圧システムの状態方程式を線形に表すこともできる。例えば、作動油の流れを層流であると仮定したり、油圧シリンダの代わりに、油圧モータを採用するのであれば、油圧システムの状態方程式を線形に表すことができる。
【０１６６】
状態方程式を第１原理モデリングではなく、ブラックボックスモデリングで導出し、その係数を推定するようにしてもよい。
【０１６７】
第２の実施の形態による油圧システムのパラメータ推定方法によって得られる推定値は、油圧システムの制御に利用してもよいし、油圧システムのオペレータへ提示することにより、オペレータによる油圧システムの操作の支援及び／又はガイダンスを行うようにしてもよい。
【０１６８】
以上、本発明の実施の形態について詳述してきたが、これらはあくまでも例示であって、本発明は、上述の実施の形態の記載によって、何等、限定的に解釈されるものではない。
【符号の説明】
【０１６９】
１０  油圧システム
１２  油圧ポンプ
１４  方向流量制御弁
１６  油圧シリンダ
１８  作動油タンク
２０  リリーフ弁
【図１】
2019057198.tif 000029
【図２】
2019057198.tif 000030
【図３】
2019057198.tif 000031
【図４】
2019057198.tif 000032
【図５】
2019057198.tif 000033
【図６】
2019057198.tif 000034
【図７】
2019057198.tif 000035
【図８】
2019057198.tif 000036
【図９】
2019057198.tif 000037
【図１０】
2019057198.tif 000038
【図１１】
2019057198.tif 000039
【図１２】
2019057198.tif 000040
【図１３】
2019057198.tif 000041
【図１４】
2019057198.tif 000042
【図１５】
2019057198.tif 000043
【図１６】
2019057198.tif 000044
【図１７】
2019057198.tif 000045
【図１８】
2019057198.tif 000046
【図１９】
2019057198.tif 000047
【図２０】
2019057198.tif 000048
【図２１】
2019057198.tif 000049
【図２２】
2019057198.tif 000050
【図２３】
2019057198.tif 000051
【図２４】
2019057198.tif 000052