特許 6815299 油圧システムのパラメータ推定方法

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

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

(11)【特許番号】6815299

(24)【登録日】2020年12月24日

(45)【発行日】2021年1月20日

(54)【発明の名称】油圧システムのパラメータ推定方法

(51)【国際特許分類】

   G05B 23/02 20060101AFI20210107BHJP

【ＦＩ】

   G05B23/02 P

【請求項の数】12

【全頁数】29

(21)【出願番号】特願2017-182179(P2017-182179)

(22)【出願日】2017年9月22日

(65)【公開番号】特開2019-57198(P2019-57198A)

(43)【公開日】2019年4月11日

【審査請求日】2019年9月30日

(73)【特許権者】

【識別番号】000001199

【氏名又は名称】株式会社神戸製鋼所

(74)【代理人】

【識別番号】100067828

【弁理士】

【氏名又は名称】小谷 悦司

(74)【代理人】

【識別番号】100115381

【弁理士】

【氏名又は名称】小谷 昌崇

(74)【代理人】

【識別番号】100109058

【弁理士】

【氏名又は名称】村松 敏郎

(72)【発明者】

【氏名】山田 崇

(72)【発明者】

【氏名】西田 吉晴

【審査官】 堀内 亮吾

(56)【参考文献】

【文献】 特開２０１５−０７７６４３（ＪＰ，Ａ）

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

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

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

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

Ｇ０５Ｂ ２３／００ −２３／０２

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

【請求項1】

作動油を吐出する油圧ポンプと、シリンダへッド側の液室及びロッド側の液室を有して前記油圧ポンプから吐出される作動油により駆動される油圧シリンダと、前記油圧ポンプから前記油圧シリンダに供給される作動油の方向及び流量を変化させる方向流量制御弁と、前記油圧ポンプと前記方向流量制御弁とを接続するポンプ配管と、前記方向流量制御弁と前記シリンダヘッド側の液室とを接続するへッド側配管と、前記ロッド側の液室と前記方向流量制御弁とを接続するロッド側配管と、を備えた油圧システムの動特性に影響を与える物理量として前記油圧システムにおける作動油の体積弾性率を推定対象パラメータとして推定する油圧システムのパラメータ推定方法であって、
前記油圧システムの状態量と前記体積弾性率とを含み、前記油圧システムの状態量として前記ポンプ配管内の圧力、前記へッド側配管内の圧力、前記ロッド側配管内の圧力、前記油圧シリンダにおけるスプール位置及びスプール速度を含む状態ベクトルを設定する工程と、
前記状態ベクトルを用いて、前記油圧システムの非線形の状態方程式を設定する工程と、
前記ポンプ配管内の圧力、前記へッド側配管内の圧力、前記ロッド側配管内の圧力、及び前記スプール位置を観測可能な物理量とし、ベイズ推定に基づく非線形の状態推定手法を前記状態方程式に適用して、前記状態ベクトルを推定することにより、前記体積弾性率を推定する工程とを備える、油圧システムのパラメータ推定方法。

【請求項2】

請求項１に記載の油圧システムのパラメータ推定方法であって、
前記状態方程式が雑音項を含み、
前記推定対象パラメータを推定する工程では、前記状態方程式にカルマンフィルタを適用して、前記状態ベクトルを推定することにより、前記推定対象パラメータを推定する、油圧システムのパラメータ推定方法。

【請求項3】

作動油を吐出する油圧ポンプと、シリンダへッド側の液室及びロッド側の液室を有して前記油圧ポンプから吐出される作動油により駆動される油圧シリンダと、前記油圧ポンプから前記油圧シリンダに供給される作動油の方向及び流量を変化させる方向流量制御弁と、前記油圧ポンプと前記方向流量制御弁とを接続するポンプ配管と、前記方向流量制御弁と前記シリンダヘッド側の液室とを接続するへッド側配管と、前記ロッド側の液室と前記方向流量制御弁とを接続するロッド側配管と、を備えた油圧システムの動特性に影響を与える物理量として前記油圧システムにおける作動油の体積弾性率及び当該体積弾性率の圧力依存性を示すパラメータの少なくとも一方を推定対象パラメータとして推定する油圧システムのパラメータ推定方法であって、
前記油圧システムの動特性を示す第１方程式を導出する工程と、
前記推定パラメータを含む状態ベクトルを用いた第２方程式であって前記第１方程式のうち前記ポンプ配管内の圧力、前記へッド側配管内の圧力及び前記ロッド側配管内の圧力をそれぞれ微分したものについての方程式を前記推定対象パラメータに対して線形な形で離散化した第２方程式を導出する工程と、
逐次推定手法を前記第２方程式に適用することにより、前記推定対象パラメータを推定する工程とを備える、油圧システムのパラメータ推定方法。

【請求項4】

請求項３に記載の油圧システムのパラメータ推定方法であって、
前記推定対象パラメータを推定する工程では、前記第２方程式に逐次最小二乗法を適用することにより、前記推定対象パラメータを推定する、油圧システムのパラメータ推定方法。

【請求項5】

請求項４に記載の油圧システムのパラメータ推定方法であって、さらに、
前記油圧システムの状態量を測定する工程を備え、
前記推定対象パラメータを推定する工程では、前記状態量の測定値に含まれるノイズに起因する前記推定対象パラメータの推定誤差を小さくするために、前記逐次最小二乗法を前記状態方程式に適用するに際して、前記状態量の測定値を補正する、油圧システムのパラメータ推定方法。

【請求項6】

請求項５に記載の油圧システムのパラメータ推定方法であって、
前記推定対象パラメータを推定する工程では、前記状態量の測定値の絶対値が大きいほど、その重みを大きくするように、前記状態量の測定値を補正する、油圧システムのパラメータ推定方法。

【請求項7】

請求項５に記載の油圧システムのパラメータ推定方法であって、
前記推定対象パラメータを推定する工程では、前記状態量の測定値の変化量が大きいほど、その重みを小さくするように、前記状態量の測定値を補正する、油圧システムのパラメータ推定方法。

【請求項8】

請求項５に記載の油圧システムのパラメータ推定方法であって、
前記推定対象パラメータを推定する工程では、前記ノイズをフィルタによって除去することにより、前記状態量の測定値を補正する、油圧システムのパラメータ推定方法。

【請求項9】

請求項３に記載の油圧システムのパラメータ推定方法であって、
前記第２方程式が雑音項を含み、
前記推定対象パラメータを推定する工程では、前記第２方程式にカルマンフィルタを適用することにより、前記推定対象パラメータを推定する、油圧システムのパラメータ推定方法。

【請求項10】

請求項３〜９の何れか１項に記載の油圧システムのパラメータ推定方法であって、
前記推定対象パラメータは、前記油圧システムにおける作動油の体積弾性率である、油圧システムのパラメータ推定方法。

【請求項11】

請求項３〜９の何れか１項に記載の油圧システムのパラメータ推定方法であって、
前記推定対象パラメータは、前記油圧システムにおける作動油の体積弾性率と縮流係数である、油圧システムのパラメータ推定方法。

【請求項12】

請求項３〜９の何れか１項に記載の油圧システムのパラメータ推定方法であって、
前記推定対象パラメータは、前記油圧システムにおける作動油の体積弾性率の圧力依存性を示すものである、油圧システムのパラメータ推定方法。

【発明の詳細な説明】

【技術分野】

【0001】

本発明は、油圧システムの動特性に影響を与える物理量をパラメータとして推定する油圧システムのパラメータ推定方法に関する。

【背景技術】

【0002】

油圧システムの制御において、油圧システムの動特性を正確にモデル化することは重要である。油圧システムの動特性は、作動油の体積弾性率や粘度、レイノルズ数等の様々な物理量に大きく依存する。そして、これらの物理量には、油圧システムの稼働中に動的に変化するパラメータに依存するものがある。例えば、作動油の体積弾性率は、作動油への気泡混入率、作動油の温度や圧力等に依存する。そのため、油圧システムの動特性に影響を与える物理量については、その値をリアルタイムで取得することが望ましい。

【0003】

例えば、作動油の体積弾性率をリアルタイムで取得する方法として、油圧システムの作動油の体積変化とそれに対する発生圧力を測定し、その比率（つまり、体積弾性率の定義そのもの）を計算することが考えられる。しかしながら、一般的な油圧システムにおいて、作動油の体積変化を正確に測定するためには、通常の作業を中断したり、作動油の体積変化とそれに対する発生圧力を測定するための回路を別途設ける必要がある。そのため、油圧システムの運用効率の低下や運用コストの増大を招いてしまう。

【0004】

作動油の体積弾性率をリアルタイムで取得する他の方法として、体積弾性率を数学関数でモデル化し、その関数の変数である物理パラメータをセンサにて測定し、代入することが考えられる。しかしながら、作動油の体積弾性率が依存する作動油への気泡混入率は、センサによって直接測定することが難しく、そのため、全ての状態に対して十分な精度を保証する数学モデルを構築することは困難である。

【0005】

そこで、作動油の体積弾性率をリアルタイムで取得する方法として、システム同定理論の適用が考えられる。システム同定理論を用いたパラメータ同定法として、例えば、最小二乗法がある。

【0006】

下記特許文献１には、最小二乗法に基づくパラメータ同定法が開示されている。特許文献１に記載のパラメータ同定法では、鉛直多関節油圧マニピュレータの非線形モデルを表す状態空間方程式を定め、当該状態空間方程式に、従来の逐次同定法を修正した新たな逐次同定法を適用する。この逐次同定法では、各ステップにおいて、測定したパラメータ及び既知のパラメータを用いて、状態空間方程式に基づき、体積弾性率及び流量係数を含む未知のパラメータについての線形方程式を立てるが、当該線形方程式を用いた同定計算は行わない。各ステップの終了後に、未知のパラメータについての線形方程式を、全ステップを横断して２つの線形方程式に厳密に統合し、これらを各時刻に亘って連立することで２つの正規方程式を導く。そして、これら２つの正規方程式を解くことで、体積弾性率及び流量係数を含む未知のパラメータを同定する。

【先行技術文献】

【特許文献】

【0007】

【特許文献1】特開２０１５−７７６４３号公報

【発明の概要】

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

【0008】

しかしながら、特許文献１に記載のパラメータ同定法は、状態の全てが観測可能であるとともに状態方程式が同定パラメータに対して線形に変形できる場合にしか、適用できない。加えて、同定パラメータの数が増えると、計算コストが指数関数的に増加してしまい、リアルタイムでの推定が困難になる。したがって、パラメータ同定法の適用範囲が著しく限定される。

【0009】

本発明の目的は、油圧システムの動特性に影響を与える物理量をパラメータとして推定するパラメータ推定方法の適用範囲を広げることである。

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

【0010】

本願の発明者は、上記の目的を達成するために、推定対象の物理量を推定対象パラメータとして含む状態ベクトルを設定することに着目して検討を進めた。そして、当該状態ベクトルを用いた方程式に対して、ベイズ推定に基づく状態推定手法又は逐次推定手法を適用することにより、パラメータ推定手法の適用範囲を広げることができるという新たな知見を得るに至った。本発明は、このような知見に基づいて完成されたものである。

【0011】

本発明の第１の態様による油圧システムのパラメータ推定方法は、油圧システムの動特性に影響を与える物理量を推定対象パラメータとして推定する油圧システムのパラメータ推定方法であって、前記油圧システムの状態量と前記推定対象パラメータとを含む状態ベクトルを設定する工程と、前記状態ベクトルを用いて、前記油圧システムの状態方程式を設定する工程と、ベイズ推定に基づく状態推定手法を前記状態方程式に適用して、前記状態ベクトルを推定することにより、前記推定対象パラメータを推定する工程とを備える。

【0012】

上記第１の態様による油圧システムのパラメータ推定方法においては、ベイズ推定に基づく状態推定手法を適用するので、状態ベクトルに含まれる油圧システムの状態量の全てが観測可能でなくてもよい。その結果、パラメータ推定方法の適用範囲を広げることができる。

【0013】

上記第１の態様による油圧システムのパラメータ推定方法において、好ましくは、前記推定対象パラメータは、前記油圧システムにおける作動油の体積弾性率である。この場合、油圧システムの動作を停止することなく、作動油の体積弾性率を得ることができる。

【0014】

上記第１の態様による油圧システムのパラメータ推定方法において、前記状態方程式が非線形であってもよい。この場合、前記推定対象パラメータを推定する工程では、好ましくは、前記状態方程式に非線形の状態推定手法を適用して、前記状態ベクトルを推定することにより、前記推定対象パラメータを推定する。このような態様においては、状態方程式が非線形であっても、推定対象パラメータを推定することができるので、パラメータ推定方法の適用範囲を広げることができる。

【0015】

上記第１の態様による油圧システムのパラメータ推定方法において、好ましくは、前記状態方程式が雑音項を含み、前記推定対象パラメータを推定する工程では、前記状態方程式にカルマンフィルタを適用して、前記状態ベクトルを推定することにより、前記推定対象パラメータを推定する。この場合、油圧システムが有する雑音等の影響を考慮して、推定対象パラメータを推定することができる。

【0016】

本発明の第２の態様による油圧システムのパラメータ推定方法は、油圧システムの動特性に影響を与える物理量を推定対象パラメータとして推定する油圧システムのパラメータ推定方法であって、前記油圧システムの動特性を示す第１方程式を導出する工程と、前記推定対象パラメータを含む状態ベクトルを用いた第２方程式を前記第１方程式から導出する工程と、逐次推定手法を前記第２方程式に適用して、前記状態ベクトルを推定することにより、前記推定対象パラメータを推定する工程とを備える。

【0017】

上記第２の態様による油圧システムのパラメータ推定方法においては、逐次推定手法を適用するので、推定対象パラメータを推定する際の計算コストの増加を抑制することができる。そのため、パラメータ推定手法の適用範囲を広げることができる。

【0018】

上記第２の態様による油圧システムのパラメータ推定方法において、好ましくは、前記推定対象パラメータを推定する工程では、前記第２方程式に逐次最小二乗法を適用して、前記状態ベクトルを推定することにより、前記推定対象パラメータを推定する。この場合、あるタイムステップでの推定対象パラメータの推定値を、その前のタイムステップでの推定対象パラメータの推定値の関数として逐次的に（少ない計算コストで）導出することができる。

【0019】

上記第２の態様による油圧システムのパラメータ推定方法は、好ましくは、さらに、前記油圧システムの状態量を測定する工程を備え、前記推定対象パラメータを推定する工程では、前記状態量の測定値に含まれるノイズに起因する前記推定対象パラメータの推定誤差を小さくするために、前記逐次最小二乗法を前記第２方程式に適用するに際して、前記状態量の測定値を補正する。この場合、補正した状態量の測定値を用いて、推定対象パラメータを推定することができるので、推定対象パラメータの推定精度を向上させることができる。

【0020】

上記第２の態様による油圧システムのパラメータ推定方法において、前記推定対象パラメータを推定する工程では、前記状態量の測定値の絶対値が大きいほど、その重みを大きくするように、前記状態量の測定値を補正するようにしてもよい。この場合、状態量の測定値の絶対値の大きさに応じて、状態量の測定値が推定対象パラメータの推定誤差に与える影響を調整することができるため、推定対象パラメータの推定精度を向上させることができる。

【0021】

上記第２の態様による油圧システムのパラメータ推定方法において、前記推定対象パラメータを推定する工程では、前記状態量の測定値の変化量が大きいほど、その重みを小さくするように、前記状態量の測定値を補正するようにしてもよい。この場合、状態量の測定値の変化量の大きさに応じて、状態量の測定値が推定対象パラメータの推定誤差に与える影響を調整することができるので、推定対象パラメータの推定精度を向上させることができる。

【0022】

上記第２の態様による油圧システムのパラメータ推定方法において、前記推定対象パラメータを推定する工程では、前記ノイズをフィルタによって除去することにより、前記状態量の測定値を補正するようにしてもよい。この場合、状態量の測定値に含まれるノイズがフィルタによって除去されているので、推定対象パラメータの推定精度を向上させることができる。

【0023】

上記第２の態様による油圧システムのパラメータ推定方法において、好ましくは、前記第２方程式が雑音項を含み、前記推定対象パラメータを推定する工程では、前記第２方程式にカルマンフィルタを適用して、前記状態ベクトルを推定することにより、前記推定対象パラメータを推定する。この場合、油圧システムが有する雑音等の影響を考慮して、推定対象パラメータを推定することができる。

【0024】

上記第２の態様による油圧システムのパラメータ推定方法において、好ましくは、前記推定対象パラメータは、前記油圧システムにおける作動油の体積弾性率である。この場合、油圧システムの動作を停止しなくても、作動油の体積弾性率を得ることができる。

【0025】

上記第２の態様による油圧システムのパラメータ推定方法において、好ましくは、前記推定対象パラメータは、前記油圧システムにおける作動油の体積弾性率と縮流係数である。この場合、作動油の体積弾性率だけでなく、縮流係数も推定することができる。

【0026】

上記第２の態様による油圧システムのパラメータ推定方法において、好ましくは、前記推定対象パラメータは、前記油圧システムにおける作動油の体積弾性率の圧力依存性を示すものである。この場合、作動油の体積弾性率が作動油の圧力に依存する度合いを推定することができる。

【発明の効果】

【0027】

本発明の油圧システムのパラメータ推定方法によれば、パラメータ推定手法の適用範囲を広げることができる。

【図面の簡単な説明】

【0028】

【図1】本発明の第１の実施の形態による油圧システムのパラメータ推定方法を示すフローチャートである。

【図2】本発明の第１の実施の形態による油圧システムのパラメータ推定方法が適用される油圧システムを示す回路図である。

【図3】本発明の第１の実施の形態で用いるｆ_ｈｉｎ（ｘ_ｖ）、ｆ_ｈｏｕｔ（ｘ_ｖ）、ｆ_ｒｉｎ（ｘ_ｖ）及びｆ_ｒｏｕｔ（ｘ_ｖ）を示すグラフである。

【図4】油圧システムへの入力としての方向流量制御弁のスプール位置を示すグラフである。

【図5】本発明の第１の実施の形態での推定結果を示すグラフであって、作動油の体積弾性率Ｋの推定を開始してからの時間変化を示すグラフである。

【図6】本発明の第１の実施の形態での推定結果を示すグラフであって、圧力Ｐ_ｐの推定を開始してからの時間変化を示すグラフである。

【図7】本発明の第１の実施の形態での推定結果を示すグラフであって、圧力Ｐ_ｈの推定を開始してからの時間変化を示すグラフである。

【図8】本発明の第１の実施の形態での推定結果を示すグラフであって、圧力Ｐ_ｒの推定を開始してからの時間変化を示すグラフである。

【図9】本発明の第１の実施の形態での推定結果を示すグラフであって、スプール位置ｘ_ｃの推定を開始してからの時間変化を示すグラフである。

【図10】本発明の第１の実施の形態での推定結果を示すグラフであって、スプール速度ｘ_ｃの推定を開始してからの時間変化を示すグラフである。

【図11】本発明の第２の実施の形態による油圧システムのパラメータ推定方法を示すフローチャートである。

【図12】逐次最小二乗法によって作動油の体積弾性率Ｋを推定した結果を示すグラフである。

【図13】逐次最小二乗法によって作動油の体積弾性率Ｋを推定した結果を示すグラフである。

【図14】逐次最小二乗法によって縮流係数Ｃを推定した結果を示すグラフである。

【図15】重み付き逐次最小二乗法によって作動油の体積弾性率Ｋを推定した結果を示すグラフである。

【図16】重み付き逐次最小二乗法によって縮流係数Ｃを推定した結果を示すグラフである。

【図17】重み付き逐次最小二乗法によって作動油の体積弾性率Ｋを推定した結果を示すグラフである。

【図18】重み付き逐次最小二乗法によって縮流係数Ｃを推定した結果を示すグラフである。

【図19】周波数フィルタと逐次最小二乗法によって作動油の体積弾性率Ｋを推定した結果を示すグラフである。

【図20】周波数フィルタと逐次最小二乗法によって縮流係数Ｃを推定した結果を示すグラフである。

【図21】大気圧における体積弾性率Ｋ０の推定結果を示すグラフである。

【図22】比例係数Ｋｐの推定結果を示すグラフである。

【図23】大気圧における体積弾性率Ｋ０の推定結果を示すグラフである。

【図24】比例係数Ｋｐの推定結果を示すグラフである。

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

【0029】

以下、添付図面を参照しながら、本発明の実施の形態について詳述する。

【0030】

［第１の実施の形態］
図１を参照して、本発明の第１の実施の形態による油圧システムのパラメータ推定方法について説明する。図１は、第１の実施の形態による油圧システムのパラメータ推定方法を示すフローチャートである。

【0031】

第１の実施の形態による油圧システムのパラメータ推定方法は、油圧システムの動特性に影響を与える物理量を推定対象パラメータとして推定するものである。第１の実施の形態による油圧システムのパラメータ推定方法は、油圧システムの動特性を示す方程式を設定する工程（ステップＳ１０）と、油圧システムの状態量と推定対象パラメータとを含む状態ベクトルを設定する工程（ステップＳ１１）と、状態ベクトルを用いて、油圧システムの状態方程式と観測方程式を設定する工程（ステップＳ１２）と、ベイズ推定に基づく状態推定手法を状態方程式に適用して、状態ベクトルを推定することにより、推定対象パラメータを推定する工程（ステップＳ１３）とを備える。

【0032】

このようなパラメータ推定手法は、例えば、制御ＣＰＵやシミュレータによって実行される。制御ＣＰＵは、例えば、実機システムを制御する。

【0033】

第１の実施の形態による油圧システムのパラメータ推定方法は、例えば、図２に示す油圧システム１０に適用して、油圧システム１０の作動油の体積弾性率Ｋを推定対象パラメータとして推定するために用いられる。なお、第１の実施の形態による油圧システムのパラメータ推定方法が適用される油圧システムは、図２に示す油圧システム１０に限定されるものではない。また、第１の実施の形態による油圧システムのパラメータ推定方法によって推定される推定対象パラメータは、油圧システム１０の作動油の体積弾性率Ｋに限定されない。

【0034】

図２を参照して、第１の実施の形態による油圧システムのパラメータ推定方法が適用される油圧システム１０について説明する。図２は、油圧システム１０を示す回路図である。

【0035】

油圧システム１０は、スプール式の方向流量制御弁１４のスプール変位ｘｖを制御入力とし、油圧シリンダ１６の動作を制御するものである。油圧システム１０は、油圧ポンプ１２と、方向流量制御弁１４と、油圧シリンダ１６と、作動油タンク１８と、リリーフ弁２０とを備える。

【0036】

油圧ポンプ１２は、作動油を吐出する。方向流量制御弁１４は、油圧ポンプ１２から油圧シリンダ１６に供給される作動油の方向及び流量を変化させる。油圧シリンダ１６は、油圧ポンプ１２が吐出する作動油（つまり、油圧ポンプ１２から供給される作動油）によって駆動される。作動油タンク１８は、油圧ポンプ１２から吐出される作動油を貯留する。リリーフ弁２０は、油圧ポンプ１２から吐出される作動油の圧力が所定の圧力を超えないように、開閉動作を行う。油圧ポンプ１２と方向流量制御弁１４とは、配管２２によって接続されている。方向流量制御弁１４と油圧シリンダ１６のシリンダヘッド側の液室とは、配管２４によって接続されている。油圧シリンダ１６のロッド側の液室と方向流量制御弁１４とは、配管２６によって接続されている。

【0037】

再び、図１を参照しながら、第１の実施の形態による油圧システムのパラメータ推定方法の各工程について説明する。

【0038】

先ず、ステップＳ１０において、油圧システム１０の動特性を示す方程式を設定する。油圧システム１０の動特性は、以下に示す方程式によって表される。

【0039】

【数1】

【0040】

ここで、ｍは油圧シリンダ１６及び負荷対象の質量を示し、ｘ_ｃはシリンダ位置（油圧シリンダ１６内でのシリンダヘッドの位置）を示し、Ａ_ｈは油圧シリンダ１６のシリンダヘッド側の液室１６１の断面積を示し、Ａ_ｒは油圧シリンダ１６のロッド側の液室１６２の断面積を示し、Ｐ_ｐは配管２２内の圧力を示し、Ｐ_ｈは配管２４内の圧力を示し、Ｐ_ｒは配管２６内の圧力を示し、Ｐ_ｔは作動油タンク１８内の圧力を示し、Ｋは作動油の体積弾性率を示し、Ｖ_ｐは配管２２の容積を示し、Ｖ_ｈは配管２４の容積を示し、Ｖ_ｒは配管２６の容積を示し、ｌは油圧シリンダ１６の全長を示す。Ｑ_ｐは油圧ポンプ１２から吐出される油量を示し、Ｑ_ｈｉｎは方向流量制御弁１４を通過して油圧シリンダ１６の液室１６１へ流入する油量を示し、Ｑ_ｈｏｕｔは油圧シリンダ１６の液室１６１から排出されて方向流量制御弁１４を通過する油量を示し、Ｑ_ｒｉｎは方向流量制御弁１４を通過して油圧シリンダ１６の液室１６２へ流入する油量を示し、Ｑ_ｒｏｕｔは油圧シリンダ１６の液室１６２から排出されて方向流量制御弁１４を通過する油量を示し、Ｑ_ｐＲはリリーフ弁２０を通過する油量を示す。

【0041】

Ｑ_ｐ、Ｑ_ｈｉｎ、Ｑ_ｈｏｕｔ、Ｑ_ｒｉｎ、Ｑ_ｒｏｕｔ、Ｑ_ｐＲは、以下に示す方程式によって表される。

【0042】

【数2】

【0043】

ここで、ｆ_ｈｉｎ（ｘ_ｖ）はＡ_ｈｉｎとｘ_ｖとの関係を示す任意の関数であり、ｆ_ｈｏｕｔ（ｘ_ｖ）はＡ_ｈｏｕｔとｘ_ｖとの関係を示す任意の関数であり、ｆ_ｒｉｎ（ｘ_ｖ）はＡ_ｒｉｎとｘ_ｖとの関係を示す任意の関数であり、ｆ_ｒｏｕｔ（ｘ_ｖ）はＡ_ｒｏｕｔとｘ_ｖとの関係を示す任意の関数であり、ｆ_ｐＲ（Ｐ_ｐ−Ｐ_ｔ）は、Ａ_ｐＲとＰ_ｐ−Ｐ_ｔとの関係を示す任意の関数である。また、Ｃは縮流係数を示し、ρは作動油の密度を示し、Ａ_ｈｉｎは方向流量制御弁１４において作動油が液室１６１に向かって流れる流路での開口面積を示し、Ａ_ｈｏｕｔは方向流量制御弁１４において液室１６１からの作動油が流れる流路での開口面積を示し、Ａ_ｒｉｎは方向流量制御弁１４において作動油が液室１６２に向かって流れる流路での開口面積を示し、Ａ_ｒｏｕｔは方向流量制御弁１４において液室１６２からの作動油が流れる流路での開口面積を示し、Ａ_ｐＲはリリーフ弁２０の開口面積を示し、ｘ_ｖは方向流量制御弁１４のスプール位置を示す。

【0044】

油圧システム１０の動特性を示す方程式を設定した後、ステップＳ１１において、油圧システム１０の状態量と推定対象パラメータとを含む状態ベクトルを設定する。

【0045】

油圧システム１０の状態量は、圧力Ｐ_ｐ、圧力Ｐ_ｈ、圧力Ｐ_ｒ、スプール位置ｘ_ｃ及びスプール速度ｘ_ｃである。また、前述のとおり、推定対象パラメータは、体積弾性率Ｋである。油圧システム１０の状態ベクトルｘは、以下のように表される。

【0046】

【数3】

【0047】

油圧システム１０の状態ベクトルｘを設定した後、ステップＳ１２において、状態ベクトルｘを用いて、油圧システム１０の状態方程式と観測方程式を設定する。

【0048】

油圧システム１０の状態方程式は、非線形である。油圧システム１０の状態方程式は、以下のように表される。

【0049】

【数4】

【0050】

つまり、油圧システム１０の状態方程式は、１階の微分方程式の形で表される。

【0051】

なお、上記状態方程式におけるｖ∈Ｒ^６×１は、平均ベクトルが０∈Ｒ^６×１であり、分散共分散行列がＱ∈Ｒ^６×６である白色雑音ベクトルを表す。つまり、上記状態方程式は、雑音項を含む。

【0052】

ここで、圧力センサや位置センサにより、油圧システム１０の状態量のうち、圧力Ｐ_ｐ、圧力Ｐ_ｈ、圧力Ｐ_ｒ及びスプール位置ｘ_ｃが観測可能であるとする。このとき、油圧システム１０の観測方程式は、以下のように表される。

【0053】

【数5】

【0054】

なお、上記観測方程式におけるｗ∈Ｒ^４×４は、平均ベクトルが０∈Ｒ^４×４であり、分散共分散行列がＲ∈Ｒ^４×４である白色雑音ベクトルを表す。

【0055】

油圧システム１０の状態方程式と観測方程式を設定した後、ステップＳ１３において、ベイズ推定に基づく状態推定手法としてのUnscentedカルマンフィルタ（つまり、非線形の状態推定手法）を状態方程式に適用して、状態ベクトルｘを推定することにより、推定対象パラメータとしての作動油の体積弾性率Ｋを推定する。

【0056】

［第１の実施の形態での推定に用いる条件］
以下、本実施の形態での推定に用いる条件について説明する。また、当該条件を用いて推定した結果についても説明する。

【0057】

本実施の形態で用いる物理定数を、以下の表１に示す。

【0058】

【表1】

【0059】

また、本実施の形態で用いるｆ_ｈｉｎ（ｘ_ｖ）、ｆ_ｈｏｕｔ（ｘ_ｖ）、ｆ_ｒｉｎ（ｘ_ｖ）及びｆ_ｒｏｕｔ（ｘ_ｖ）を、図３に示す。ｆ_ｈｉｎ（ｘ_ｖ）とｆ_ｒｏｕｔ（ｘ_ｖ）は、それぞれ、スプール位置ｘ_ｖが負である場合は、開口面積がゼロであり、スプール位置ｘ_ｖが正である場合は、スプール位置ｘ_ｖが大きくなるに従って、開口面積が大きくなる。ｆ_ｒｉｎ（ｘ_ｖ）とｆ_ｈｏｕｔ（ｘ_ｖ）は、それぞれ、スプール位置ｘ_ｖが負である場合は、スプール位置ｘ_ｖの絶対値が小さくなるに従って、開口面積が小さくなり、スプール位置ｘ_ｖが正である場合は、開口面積がゼロである。油圧システム１０への入力として、方向流量制御弁１４のスプール位置ｘ_ｖを、図４に示すように与える。

【0060】

作動油の体積弾性率Ｋの時間変化ｆ_Ｋは、ゼロであってもよいし、正の値（例えば、５０）であってもよい。作動油の体積弾性率Ｋの時間変化ｆ_Ｋがゼロである場合、作動油の体積弾性率Ｋが定数であると仮定した推定が行われる。一方、作動油の体積弾性率Ｋの時間変化ｆ_Ｋが正の値である場合、作動油の体積弾性率Ｋは変動パラメータとなり、モデル化誤差を模擬することができる。なお、本実施の形態では、作動油の体積弾性率Ｋの時間変化ｆ_Ｋを５０にしている。

【0061】

このような条件の下で、作動油の体積弾性率Ｋを推定し、図５〜図１０に示すような結果を得た。

【0062】

図５に示すように、推定を開始してから０．１秒以内に、作動油の体積弾性率Ｋの推定値が、初期値（５００ＭＰａ）から真値（１５００ＭＰａ）付近に移動していることが確認できる。作動油の体積弾性率Ｋ及びシリンダ速度が直接計測されていないにも関わらず、これらの数値を推定できていることが判る。

【0063】

また、図５に示すように、推定モデルでは考慮されていない真値の時間変化に対して、推定値が追従できていることも確認できる。これにより、状態方程式にモデル化誤差が含まれていても、雑音ベクトルvの分散共分散行列Qを適切に設定することで実用的な推定が実現されていることが判る。

【0064】

また、本実施の形態では、推定値の標準偏差を指標として、推定の信頼度を評価することができる。図５では、油流の状態に応じて推定値の標準偏差が時々刻々変化しており、推定誤差を考慮に入れた制御器を設計する際の定量的な指標として利用できる。

【0065】

また、図６〜図１０に示すように、作動油の体積弾性率Ｋの推定のみならず、雑音ベクトルwの分散共分散行列Rを適切に設定することにより、雑音の入った測定値から状態ベクトルxの真値も推定できていることが確認できる。

【0066】

［第１の実施の形態の応用例］
第１の実施の形態では、ベイズ推定に基づいた状態推定手法として、Unscentedカルマンフィルタを用いていたが、例えば、線形カルマンフィルタや拡張カルマンフィルタ、粒子フィルタ、アンサンブルカルマンフィルタ等を用いてもよい。なお、線形カルマンフィルタを用いる場合には、状態方程式は線形でなければならないが、例えば、動作範囲が限られた油圧システムであれば、その範囲内で非線形な状態方程式を線形な状態方程式に近似することにより、線形カルマンフィルタを適用することができる。また、限定的な仮定の下では、油圧システムの状態方程式を線形に表すこともできる。例えば、作動油の流れを層流であると仮定したり、油圧シリンダの代わりに、油圧モータを採用するのであれば、油圧システムの状態方程式を線形に表すことができる。

【0067】

推定対象パラメータが動的に変化する場合、当該推定対象パラメータに関する既知の特性を状態方程式に反映させてもよい。例えば、作動油の体積弾性率Ｋは、等温条件において、圧力に比例するので、これを状態方程式に反映させてもよい。

【0068】

状態方程式を第１原理モデリングではなく、ブラックボックスモデリングで導出し、その係数を推定するようにしてもよい。

【0069】

第１の実施の形態による油圧システムのパラメータ推定方法によって得られる推定値は、油圧システムの制御に利用してもよいし、油圧システムのオペレータへ提示することにより、オペレータによる油圧システムの操作の支援及び／又はガイダンスを行うようにしてもよい。

【0070】

［第２の実施の形態］
図１１を参照しながら、本発明の第２の実施の形態による油圧システムのパラメータ推定方法について説明する。図１１は、第２の実施の形態による油圧システムのパラメータ推定方法を示すフローチャートである。

【0071】

第２の実施の形態による油圧システムのパラメータ推定方法は、第１の実施の形態による油圧システムのパラメータ推定方法と比べて、ベイズ推定に基づく状態推定手法の代わりに、逐次推定手法を用いている点で異なる。また、状態ベクトルが推定対象パラメータのみを含む場合もあり得る点で異なる。

【0072】

第２の実施の形態による油圧システムのパラメータ推定方法は、油圧システムの動特性を示す第１方程式を設定する工程（ステップＳ２０）と、推定対象パラメータを含む状態ベクトルを用いた第２方程式を第１方程式から導出する工程（ステップＳ２１）と、逐次推定手法を第２方程式に適用して、状態ベクトルを推定することにより、推定対象パラメータを推定する工程（ステップＳ２２）とを備える。

【0073】

このようなパラメータ推定方法は、例えば、制御ＣＰＵやシミュレータによって実行される。制御ＣＰＵは、例えば、実機システムを制御する。

【0074】

ここで、ステップＳ２０の工程は、第１の実施の形態におけるステップＳ１０の工程と同じであるから、その詳細な説明は省略する。ステップＳ２１及びステップＳ２２の工程の詳細については、後述する。

【0075】

第２の実施の形態による油圧システムのパラメータ推定方法は、第１の実施の形態による油圧システムのパラメータ推定方法と同様に、例えば、図２に示す油圧システム１０に適用して、油圧システム１０の作動油の体積弾性率Ｋを推定するために用いられる。なお、第２の実施の形態による油圧システムのパラメータ推定方法が適用される油圧システムは、図２に示す油圧システム１０に限定されるものではない。また、第２の実施の形態による油圧システムのパラメータ推定方法によって推定される推定対象パラメータは、油圧システム１０の作動油の体積弾性率Ｋに限定されない。

【0076】

逐次推定手法は、各タイムステップでの推定値がその前のタイムステップでの推定値の関数として導出できるものであれば、特に限定されない。逐次推定手法は、逐次最小二乗法やカルマンフィルタを含む。逐次推定手法を適用することにより、計算コストを少なくすることができる。

【0077】

【表2】

【0078】

表２は、本実施の形態にて採用される逐次推定手法と推定対象パラメータとの組み合わせを示す。例えば、（１）の組み合わせに係る態様は、逐次最小二乗法によって作動油の体積弾性率Ｋを推定する場合を示す。

【0079】

以下、表２に示す（１）〜（８）の組み合わせに係る態様について説明する。なお、表２における空欄は、他の組み合わせに係る態様にて説明する逐次推定手法を適用することで実施できるため、その詳細な説明を省略していることを示しているだけであり、実施が不可能であることを意味しているのではない。例えば、（５）の組み合わせに係る態様にて説明する周波数フィルタ付の逐次最小二乗法を適用しても、作動油の体積弾性率Ｋを推定することができる。

【0080】

また、以下では、表２に示す（１）〜（８）の組み合わせに係る態様の推定結果についても説明する。なお、推定に用いた条件は、第１の実施の形態での条件、つまり、表１、図３及び図４に示すものと同じであるから、その詳細な説明は省略する。

【0081】

（１）の組み合わせに係る態様について
先ず、逐次最小二乗法によって作動油の体積弾性率Ｋを推定する場合について説明する。

【0082】

この場合、ステップＳ２１において、油圧システム１０の動特性を表す方程式のうち、Ｐ_ｐ、Ｐ_ｈ及びＰ_ｒを微分したものについての方程式を推定対象パラメータとしての作動油の体積弾性率Ｋに対して線形な形で離散化する。

【0083】

先ず、Ｐ_ｐ、Ｐ_ｈ及びＰ_ｒを微分したものについての方程式の右辺を、Ｋとそれ以外の項ａ_ｐ、ａ_ｈ、ａ_ｒの積として陽に表すと、以下のようになる。

【0084】

【数6】

【0085】

これらの式の左辺、つまり、Ｐ_ｐ、Ｐ_ｈ及びＰ_ｒを微分したものを時間区切りΔｔで差分近似すると、以下のようになる。

【0086】

【数7】

【0087】

ここで、ｋは時間区切りΔｔで離散化されたタイムステップを表す。Ｋ以外のパラメータは、観測可能であるとする。

【0088】

各タイムステップｎの時点での上記式の両辺の二乗誤差の総和は、以下のようになる。

【0089】

【数8】

【0090】

このような二乗誤差の総和を最小化するＫ（ｎ）を逐次最小二乗法によって導出する。これにより、ステップＳ２２が実行される。なお、上記式において、計算開始タイムステップａは、油圧システム１０の特性に応じて調整するパラメータである。

【0091】

図１２は、作動油の体積弾性率Ｋの推定結果を示すグラフである。図１２に示すように、時間が経過するに従って、推定値が真値に近づいており、推定の精度が向上していることが判る。なお、推定値が真値に一致していないのは、Ｐ_ｐ、Ｐ_ｈ及びＰ_ｒを微分したものを１次のオイラー法で差分近似したことで生じる近似誤差が原因である。

【0092】

（２）の組み合わせに係る態様について
続いて、逐次最小二乗法によって作動油の体積弾性率Ｋと縮流係数Ｃを推定する場合について説明する。

【0093】

この場合、ステップＳ２１において、油圧システム１０の動特性を表す方程式のうち、Ｐ_ｐ、Ｐ_ｈ及びＰ_ｒを微分したものについての方程式を推定対象パラメータとしてのＫ及びＫＣに対して線形な形で離散化する。

【0094】

先ず、Ｐ_ｐ、Ｐ_ｈ及びＰ_ｒを微分したものについての方程式の右辺のうち、Ｋの項とＫＣの項とを陽に表すと、以下のようになる。

【0095】

【数9】

【0096】

これらの式の左辺、つまり、Ｐ_ｐ、Ｐ_ｈ及びＰ_ｒを微分したものを時間区切りΔｔで差分近似すると、以下のようになる。

【0097】

【数10】

【0098】

ここで、ｋは時間区切りΔｔで離散化されたタイムステップを表す。

【0099】

左辺のベクトルをｙ_ＬＳ（ｋ）、右辺のベクトルをｘ_ＬＳ、行列をＡ（ｋ）とすると、上記式は、以下のようになる。

【0100】

【数11】

【0101】

ここで、ｘ_ＬＳ以外のパラメータは観測可能であるとする。

【0102】

各タイムステップｎの時点での上記式の両辺の二乗誤差の総和は、以下のようになる。

【0103】

【数12】

【0104】

このような二乗誤差の総和を最小化するｘ_ＬＳ（ｎ）を逐次最小二乗法によって導出する。その後、ＫＣをＫで除することにより、縮流係数Ｃを推定する。これにより、ステップＳ２２が実行される。なお、上記式において、計算開始タイムステップａは、油圧システム１０の特性に応じて調整するパラメータである。

【0105】

図１３は、作動油の体積弾性率Ｋの推定結果を示すグラフである。図１４は、縮流係数Ｃの推定結果を示すグラフである。なお、図１３及び図１４では、計測値Ｐ_ｐ、Ｐ_ｈ、Ｐ_ｒにノイズがない場合の推定結果だけでなく、計測値Ｐ_ｐ、Ｐ_ｈ、Ｐ_ｒにノイズがある場合の推定結果も示している。ノイズは、計測値Ｐ_ｐ、Ｐ_ｈ、Ｐ_ｒに対して平均がゼロで標準偏差が０．０４［ＭＰａ］であるホワイトノイズが加わり、且つ、計測値ｘ_ｃに対して平均がゼロで標準偏差が１０［ｍｍ］のホワイトノイズが加わったものである。

【0106】

図１３に示すように、作動油の体積弾性率Ｋは、（１）の場合と同様に、時間が経過するに従って、推定値が真値に近づいており、推定の精度が向上していることが判る。縮流係数Ｃについても、図１４に示すように、正しく推定されていることが確認できる。

【0107】

これに対して、計測値にノイズがある場合は、Ｋ及びＣの推定結果の精度が低下している。これは、計測値Ｐ_ｐ、Ｐ_ｈ、Ｐ_ｒに含まれるノイズがＱ_ｐ、Ｑ_ｈｉｎ、Ｑ_ｈｏｕｔ、Ｑ_ｒｉｎ、Ｑ_ｒｏｕｔ、Ｑ_ｐＲについての方程式中の非線形関数（平方根）を通ることで、流量の平均値にずれを生み出すことが主な原因である。推定精度の低下は、油圧システム１０の動特性を示す方程式の非線形性が大きくなるほど、また、計測データの真値に対するノイズの大きさが大きくなるほど、顕著になる。計測値にノイズがある油圧システムに対しては、単純に逐次最小二乗法を適用することが好ましくないことが判る。

【0108】

（３）の組み合わせに係る態様について
続いて、計測値に含まれるノイズによる推定精度の低下を軽減するために、重み付き逐次最小二乗法によって作動油の体積弾性率Ｋと縮流係数Ｃを推定する場合について説明する。なお、この態様では、ステップＳ２１の工程は、（２）の組み合わせに係る態様と同じである。

【0109】

この態様では、以下に示す二乗誤差の重み付き総和を最小化するｘ_ＬＳ（ｎ）を導出する。その後、ＫＣをＫで除することにより、縮流係数Ｃを推定する。これにより、ステップＳ２２が実行される。

【0110】

【数13】

【0111】

ここで、Ｗ（ｋ）はタイムステップｋにおける重み行列であり、計測値の信頼度等、既知の性質に応じて設定される。計測値の真値に対するノイズの大きさが大きくなるほど、推定の精度が顕著に低下することから、計測値の絶対値が大きいほど、重みを大きくして、その計測値に対する信頼度を上げるような推定を行う。Ｗ（ｋ）は、例えば、以下のように設定される。

【0112】

【数14】

【0113】

ここで、α_ｐ、α_ｈ、α_ｒは、重みの大きさを調整するチューニングパラメータである。なお、Ｗ（ｋ）の設定はこれに限定されない。

【0114】

図１５は、作動油の体積弾性率Ｋの推定結果を示すグラフである。図１６は、縮流係数Ｃの推定結果を示すグラフである。図１３及び図１４と比べて、計測値にノイズが含まれているにも関わらず、推定精度が向上していることが判る。

【0115】

（４）の組み合わせに係る態様について
続いて、計測値に含まれるノイズによる推定精度の低下を軽減するために、重み付き逐次最小二乗法によって作動油の体積弾性率Ｋと縮流係数Ｃを推定する場合について説明する。（４）の組み合わせに係る態様では、（３）の組み合わせに係る態様と比べて、重み行列Ｗ（ｋ）が異なっている。この重み行列Ｗ（ｋ）は、１タイムステップにおける圧力の変化が大きいほど、重みが小さくなるように、設定されている。つまり、圧力の変化が油圧システム１０の動作として想定される値よりも大きい場合、その計測値をノイズとみなし、推定への寄与を低減させる。重み行列Ｗ（ｋ）は、例えば、以下のように設定される。

【0116】

【数15】

【0117】

ここで、αは、想定変化量を調整するチューニングパラメータである。なお、Ｗ（ｋ）の設定はこれに限定されない。

【0118】

図１７は、作動油の体積弾性率Ｋの推定結果を示すグラフである。図１８は、縮流係数Ｃの推定結果を示すグラフである。図１３及び図１４と比べて、計測値にノイズが含まれているにも関わらず、推定精度が向上していることが判る。

【0119】

（５）の組み合わせに係る態様について
続いて、計測値に含まれるノイズによる推定精度の低下を軽減するために、周波数フィルタと逐次最小二乗法によって作動油の体積弾性率Ｋと縮流係数Ｃを推定する場合について説明する。

【0120】

この場合、圧力計測値の周波数成分のうち、油圧システム１０の動作として想定される周波数成分以外をノイズとみなしてフィルタにより除去し、フィルタ処理後の計測値に対して逐次最小二乗法を適用する。使用する周波数フィルタは、以下のような移動平均フィルタである。

【0121】

【数16】

【0122】

ここで、Ｐ_{ｆｉｌｔｅｒｅｄ}はフィルタ適用後の圧力値であり、Ｐは圧力計測値（Ｐ_ｐ、Ｐ_ｈ、Ｐ_ｒ）である。

【0123】

なお、使用する周波数フィルタは、移動平均フィルタに限定されず、例えば、線形フィルタ、バターワースフィルタ、チェビシェフフィルタ等であってもよい。

【0124】

図１９は、作動油の体積弾性率Ｋの推定結果を示すグラフである。図２０は、縮流係数Ｃの推定結果を示すグラフである。図１９及び図２０では、上記式においてＮ＝２１とした移動平均フィルタを用いた推定結果を示している。図１３及び図１４と比べて、計測値にノイズが含まれているにも関わらず、推定精度が向上していることが判る。

【0125】

（６）の組み合わせに係る態様について
続いて、（６）の組み合わせに係る態様について説明する。この態様は、計測値に含まれるノイズによる推定精度の低下を軽減するために、Unscentedカルマンフィルタによって作動油の体積弾性率Ｋを推定する場合である。なお、Unscentedカルマンフィルタによって作動油の体積弾性率Ｋを推定する手法は、第１の実施の形態で説明したものと同様であるから、その詳細な説明は省略する。

【0126】

この手法によれば、油圧システム１０の状態量と推定対象パラメータ（作動油の体積弾性率Ｋ）を同時に推定することができるので、油圧システム１０の状態量の全てを計測する必要がなくなる。また、ノイズの大きさ等、既知の特性を考慮した推定を行うので、推定精度に優れる。

【0127】

（７）の組み合わせに係る態様について
続いて、（７）の組み合わせに係る態様について説明する。この態様は、逐次最小二乗法に基づく作動油の体積弾性率Ｋの圧力依存性を推定するものである。

【0128】

作動油の体積弾性率Ｋは、作動油の圧力に応じて変化することが知られている。例えば、等温且つ低圧領域では、作動油の体積弾性率Ｋは、以下のように、圧力に比例する関数として近似することができる。

【0129】

【数17】

【0130】

ここで、Ｋ_０は大気圧における体積弾性率であり、Ｋ_ｐは比例係数であり、Ｐはゲージ圧（圧力計測値）である。Ｋ_０及びＫ_ｐは作動油の温度や作動油の気泡含有率等、動的に変化するパラメータに依存する。したがって、逐次推定手法によって推定する。

【0131】

ステップＳ２１において、油圧システム１０の動特性を表す方程式のうち、Ｐ_ｐ、Ｐ_ｈ及びＰ_ｒを微分したものについての方程式を推定対象パラメータとしてのＫ_０及びＫ_ｐに対して線形な形で離散化する。

【0132】

先ず、油圧システム１０の動特性を表す方程式のうち、Ｐ_ｐ、Ｐ_ｈ及びＰ_ｒを微分したものについての方程式に対して、上記式を代入し、Ｋ_０の項とＫ_ｐの項を陽に示すと、以下のようになる。

【0133】

【数18】

【0134】

左辺、つまり、Ｐ_ｐ、Ｐ_ｈ及びＰ_ｒを微分したものを時間区切りΔｔで差分近似すると、以下のようになる。

【0135】

【数19】

【0136】

ここで、ｋは時間区切りΔｔで離散化されたタイムステップを表す。

【0137】

左辺のベクトルをｙ_ＬＳ（ｋ）、右辺のベクトルをｘ_ＬＳ、行列をＡ（ｋ）とすると、上記式は、以下のようになる。

【0138】

【数20】

【0139】

ここで、ｘ_ＬＳ以外のパラメータは観測可能であるとする。

【0140】

各タイムステップｎの時点での上記式の両辺の二乗誤差の総和は、以下のようになる。

【0141】

【数21】

【0142】

このような二乗誤差の総和を最小化するｘ_ＬＳ（ｎ）を逐次最小二乗法によって導出する。これにより、ステップＳ２２が実行される。なお、上記式において、計算開始タイムステップａは、油圧システム１０の特性に応じて調整するパラメータである。

【0143】

図２１は、大気圧における体積弾性率Ｋ０の推定結果を示すグラフである。図２２は、比例係数Ｋｐの推定結果を示すグラフである。なお、図２１及び図２２では、計測値Ｐ_ｐ、Ｐ_ｈ、Ｐ_ｒにノイズがない場合の推定結果だけでなく、計測値Ｐ_ｐ、Ｐ_ｈ、Ｐ_ｒにノイズがある場合の推定結果も示している。ノイズは、計測値Ｐ_ｐ、Ｐ_ｈ、Ｐ_ｒに対して平均がゼロで標準偏差が０．０２［ＭＰａ］であるホワイトノイズが加わり、且つ、計測値ｘ_ｃに対して平均がゼロで標準偏差が１０［ｍｍ］のホワイトノイズが加わったものである。

【0144】

図２１及び図２２に示すように、ノイズがない場合は、十分な推定精度を確保できていることが判る。

【0145】

これに対して、計測値にノイズがある場合は、Ｋ_０及びＫ_ｐの推定結果の精度が低下しており、単純に逐次最小二乗法を適用することが好ましくないことが判る。

【0146】

（８）の組み合わせに係る態様について
続いて、（８）の組み合わせに係る態様について説明する。この態様は、（７）と同様に、作動油の体積弾性率Ｋの圧力依存性を推定するものであるが、逐次最小二乗法ではなく、Unscentedカルマンフィルタを適用する。

【0147】

先ず、ステップＳ２１において、油圧システム１０の状態方程式及び観測方程式を設定する。具体的には、以下のとおりである。

【0148】

油圧システム１０の状態量（圧力Ｐ_ｐ、圧力Ｐ_ｈ、圧力Ｐ_ｒ、シリンダ位置ｘ_ｃ、シリンダ速度ｘ_ｃに推定対象パラメータであるＫ_０及びＫ_ｐを追加したものを、油圧システム１０の状態ベクトルとして、以下のように設定する。

【0149】

【数22】

【0150】

このとき、油圧システム１０の状態方程式は、以下のように表される。

【0151】

【数23】

【0152】

つまり、油圧システム１０の状態方程式は、１階の微分方程式の形で表される。

【0153】

なお、上記状態方程式におけるｖ∈Ｒ^７×１は、平均ベクトルが０∈Ｒ^７×１であり、分散共分散行列がＱ∈Ｒ^７×７である白色雑音ベクトルを表す。

【0154】

ここで、Ｋ_０の時間変化ｆ_Ｋ０及びＫ_ｐの時間変化ｆ_Ｋｐは、ゼロとする。

【0155】

また、圧力センサや位置センサにより、油圧システム１０の状態量のうち、圧力Ｐ_ｐ、圧力Ｐ_ｈ、圧力Ｐ_ｒ及びスプール位置ｘ_ｃが観測可能であるとする。このとき、油圧システム１０の観測方程式は、以下のように表される。

【0156】

【数24】

【0157】

なお、上記観測方程式におけるｗ∈Ｒ^４×４は、平均ベクトルが０∈Ｒ^４×４であり、分散共分散行列がＲ∈Ｒ^４×４である白色雑音ベクトルを表す。

【0158】

油圧システム１０の状態方程式と観測方程式を設定した後、ステップＳ２２において、Unscentedカルマンフィルタを状態方程式に適用して、状態ベクトルｘを推定することにより、油圧システム１０の状態量及び推定対象パラメータを推定する。

【0159】

図２３は、大気圧における体積弾性率Ｋ０の推定結果を示すグラフである。図２４は、比例係数Ｋｐの推定結果を示すグラフである。図２３及び図２４では、計測値Ｐ_ｐ、Ｐ_ｈ、Ｐ_ｒにノイズがない場合の推定結果だけでなく、計測値Ｐ_ｐ、Ｐ_ｈ、Ｐ_ｒにノイズがある場合の推定結果も示している。ノイズは、計測値Ｐ_ｐ、Ｐ_ｈ、Ｐ_ｒに対して平均がゼロで標準偏差が０．０２［ＭＰａ］であるホワイトノイズが加わり、且つ、計測値ｘ_ｃに対して平均がゼロで標準偏差が１０［ｍｍ］のホワイトノイズが加わったものである。

【0160】

図２３及び図２４では、図２１及び図２２と比べて、推定対象パラメータＫ_０、Ｋ_ｐ及びシリンダ速度ｘ_ｃが直接計測されていないにも関わらず、推定精度が向上していることが判る。

【0161】

［第２の実施の形態の応用例］
第２の実施の形態では、過去の計測値への依存性を調整するパラメータとして、計算開始ステップａを用いたが、忘却係数を用いて指数関数的に過去の計測値の重みを減少させてもよい。

【0162】

第２の実施の形態では、測定値の絶対値に応じた重みを設定したが、測定値の絶対値が所定の閾値よりも小さい場合に、その測定値を推定に用いないようにしてもよい。

【0163】

第２の実施の形態では、作動油の体積弾性率を圧力に関して線形な関数で近似したが、圧力に関して非線形な関数で近似してもよい。

【0164】

第２の実施の形態では、作動油の体積弾性率を圧力の関数としてモデル化したが、任意のパラメータの関数としてモデル化し、その係数を求めるようにしてもよい。例えば、作動油の体積弾性率を、温度の関数として、或いは、温度及び圧力の関数として、モデル化してもよい。

【0165】

第２の実施の形態では、Unscentedカルマンフィルタを用いたが、例えば、線形カルマンフィルタや拡張カルマンフィルタ、粒子フィルタ等を用いてもよい。なお、線形カルマンフィルタを用いる場合には、状態方程式は線形でなければならないが、例えば、動作範囲が限られた油圧システムであれば、その範囲内で非線形な状態方程式を線形な状態方程式に近似することにより、線形カルマンフィルタを適用することができる。また、限定的な仮定の下では、油圧システムの状態方程式を線形に表すこともできる。例えば、作動油の流れを層流であると仮定したり、油圧シリンダの代わりに、油圧モータを採用するのであれば、油圧システムの状態方程式を線形に表すことができる。

【0166】

状態方程式を第１原理モデリングではなく、ブラックボックスモデリングで導出し、その係数を推定するようにしてもよい。

【0167】

第２の実施の形態による油圧システムのパラメータ推定方法によって得られる推定値は、油圧システムの制御に利用してもよいし、油圧システムのオペレータへ提示することにより、オペレータによる油圧システムの操作の支援及び／又はガイダンスを行うようにしてもよい。

【0168】

以上、本発明の実施の形態について詳述してきたが、これらはあくまでも例示であって、本発明は、上述の実施の形態の記載によって、何等、限定的に解釈されるものではない。

【符号の説明】

【0169】

１０ 油圧システム
１２ 油圧ポンプ
１４ 方向流量制御弁
１６ 油圧シリンダ
１８ 作動油タンク
２０ リリーフ弁

【図1】

【図2】

【図3】

【図4】

【図5】

【図6】

【図7】

【図8】

【図9】

【図10】

【図11】

【図12】

【図13】

【図14】

【図15】

【図16】

【図17】

【図18】

【図19】

【図20】

【図21】

【図22】

【図23】

【図24】