(57)【特許請求の範囲】
【請求項１】
過電圧モデルとヒステリシスモデルとを備えるバッテリの等価回路モデルのパラメータを推定するバッテリのパラメータ推定装置において、
前記バッテリに第１の振幅を有する第１の電流と、前記第１の振幅よりも小さい第２の振幅を有する第２の電流とをそれぞれ入力し、
前記第１の電流に応じた前記バッテリの出力に基づいて、前記過電圧モデルに係るパラメータを推定し、
前記第２の電流に応じた前記バッテリの出力に基づいて、前記ヒステリシスモデルに係るパラメータを推定することを特徴とするバッテリのパラメータ推定装置。
【請求項２】
請求項１に記載のバッテリのパラメータ推定装置において、
前記第１の電流に応じた前記バッテリの出力に基づいて、前記ヒステリシスモデルに係るパラメータを固定しつつ、前記過電圧モデルに係るパラメータを推定することを特徴とするバッテリのパラメータ推定装置。
【請求項３】
請求項１又は２に記載のバッテリのパラメータ推定装置において、
前記第２の電流に応じた前記バッテリの出力に基づいて、前記過電圧モデルに係るパラメータを固定しつつ、前記ヒステリシスモデルに係るパラメータを推定することを特徴とするバッテリのパラメータ推定装置。
【請求項４】
請求項１乃至３いずれか一項に記載のバッテリのパラメータ推定装置において、
前記ヒステリシスモデルに係るパラメータの推定と、前記過電圧モデルに係るパラメータの推定とを交互に繰り返すことを特徴とするバッテリのパラメータ推定装置。
【請求項５】
請求項１乃至４いずれか一項に記載のバッテリのパラメータ推定装置において、
前記バッテリに入力される電流は、Ｍ系列信号であることを特徴とするバッテリのパラメータ推定装置。
【請求項６】
請求項１乃至５いずれか一項に記載のバッテリのパラメータ推定装置において、
前記ヒステリシスモデルに係るパラメータは、前記ヒステリシスモデルにおける、電圧降下の最大範囲及び入力電流の特性により関連付けられるヒステリシス抵抗を含むことを特徴とするバッテリのパラメータ推定装置。
【請求項７】
請求項１乃至６いずれか一項に記載のバッテリのパラメータ推定装置において、
前記ヒステリシスモデルに係るパラメータは、前記ヒステリシスモデルにおける、電圧降下の最大範囲及び電圧降下の速さの特性により関連付けられるヒステリシス容量を含むことを特徴とするバッテリのパラメータ推定装置。
【請求項８】
請求項６に記載のバッテリのパラメータ推定装置において、
前記電圧降下の最大範囲をＭ（ｔ）、前記入力電流をｕ（ｔ）とした場合に、前記ヒステリシス抵抗Ｒh（ｔ）を式
【数１】
0006711981.tif 000036
で表すことを特徴とするバッテリのパラメータ推定装置。
【請求項９】
請求項７に記載のバッテリのパラメータ推定装置において、
前記電圧降下の最大範囲をＭ（ｔ）、前記電圧降下の速さをΓ（ｔ）とした場合に、前記ヒステリシス容量Ｃh（ｔ）を式
【数２】
0006711981.tif 000037
で表すことを特徴とするバッテリのパラメータ推定装置。
【請求項１０】
請求項８に記載のバッテリのパラメータ推定装置において、
前記ヒステリシス抵抗Ｒh（ｔ）を式
【数３】
0006711981.tif 000038
で表すことを特徴とするバッテリのパラメータ推定装置。
【請求項１１】
請求項８に記載のバッテリのパラメータ推定装置において、
前記ヒステリシス抵抗Ｒh（ｔ）を式
【数４】
0006711981.tif 000039
で表すことを特徴とするバッテリのパラメータ推定装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、バッテリのパラメータ推定装置に関する。
【背景技術】
【０００２】
バッテリの内部状態又は等価回路モデルのパラメータを推定する装置として、例えば特許文献１に記載のものが知られている。特許文献１に記載のバッテリのパラメータ推定装置は、バッテリの充放電電流および端子電圧を検出し、これらを入力として、抵抗と容量を含むバッテリの等価回路モデルを用いてカルマンフィルタでそのパラメータやバッテリの内部状態、及び開回路電圧を推定（算出）する。
【０００３】
上述のバッテリの等価回路モデルにおいて、過電圧に関するモデル（過電圧モデル）を用いてバッテリのＳＯＣ−ＯＣＶ特性を表すことができる。しかしながら、実際のバッテリでは充電後と放電後とでＳＯＣ−ＯＣＶ特性が異なるヒステリシス現象が発生することがある。この場合にはバッテリのＳＯＣ−ＯＣＶ特性を正確に表すことができない。ヒステリシス現象は電極の材料により発生し、特にリン酸リチウムを用いた場合にはヒステリシス現象の影響が大きく出る。
【０００４】
ここで、バッテリのヒステリシス現象を取り扱うためにバッテリの等価回路にヒステリシスによる電圧降下を表すヒステリシス素子を付け加えたモデル（ヒステリシスモデル）が提案されている。バッテリのヒステリシス現象とは、バッテリの充放電に伴う状態の変動において、その変動の履歴によってバッテリの平衡状態が変動することを言う。ヒステリシスのないバッテリでは、バッテリの充放電の履歴にかかわらず、ある程度の時間放置することによって元の平衡状態に戻る。一方で、ヒステリシスのあるバッテリでは、バッテリの充放電の履歴によっては、いくら放置しても元の平衡状態に戻らないことがある。
【０００５】
ヒステリシスモデルを考慮して、バッテリの内部状態又は等価回路モデルのパラメータを推定する装置として、例えば特許文献２、非特許文献１に記載のものが知られている。特許文献２、非特許文献１では、バッテリのヒステリシスモデルに係るパラメータを、ＳＯＣ等とあわせてカルマンフィルタで同時推定している。
【先行技術文献】
【特許文献】
【０００６】
【特許文献１】特開２０１４−７４６８２号公報
【特許文献２】特許第４５１１６００号公報
【非特許文献】
【０００７】
【非特許文献１】G. L. Plett: "Extended Kalman filtering for battery management systems of LiPB-based HEV battery packs Part 2. Modeling and identification", Journal of Power Sources 134 (2004) 262-276
【発明の概要】
【発明が解決しようとする課題】
【０００８】
しかし、特許文献２、非特許文献１では、ヒステリシスの性質を考慮せずに任意の入力をバッテリの等価回路に加えて、ヒステリシスモデルに係るパラメータを推定している。この場合、ヒステリシスによって生じる電圧降下が相対的に小さい場合にもヒステリシスモデルに係るパラメータを推定することで、推定精度が落ちることがある。一方で、過電圧モデルによって生じる電圧降下が相対的に小さい場合に過電圧モデルに係るパラメータを推定することで、推定精度が落ちることがある。また、任意の入力から推定するために、カルマンフィルタとして、ＵＫＦ（Ｕｎｓｃｅｎｔｅｄ  Ｋａｌｍａｎ  Ｆｉｌｔｅｒ）等の高度なものを用いる必要がある。
【０００９】
かかる事情に鑑みてなされた本発明の目的は、ヒステリシスモデルに係るパラメータの推定精度と過電圧モデルに係るパラメータの推定精度とをそれぞれ向上させることができるバッテリのパラメータ推定装置を提供することにある。
【課題を解決するための手段】
【００１０】
上記課題を解決するために、第１の観点に係るバッテリのパラメータ推定装置は、
過電圧モデルとヒステリシスモデルとを備えるバッテリの等価回路モデルのパラメータを推定するバッテリのパラメータ推定装置において、
前記バッテリに第１の振幅を有する第１の電流と、前記第１の振幅よりも小さい第２の振幅を有する第２の電流とをそれぞれ入力し、
前記第１の電流に応じた前記バッテリの出力に基づいて、前記過電圧モデルに係るパラメータを推定し、
前記第２の電流に応じた前記バッテリの出力に基づいて、前記ヒステリシスモデルに係るパラメータを推定することを特徴とする。
【００１１】
上記課題を解決するために、第２の観点に係るバッテリのパラメータ推定装置は、
前記第１の電流に応じた前記バッテリの出力に基づいて、前記ヒステリシスモデルに係るパラメータを固定しつつ、前記過電圧モデルに係るパラメータを推定することを特徴とする。
【００１２】
上記課題を解決するために、第３の観点に係るバッテリのパラメータ推定装置は、
前記第２の電流に応じた前記バッテリの出力に基づいて、前記過電圧モデルに係るパラメータを固定しつつ、前記ヒステリシスモデルに係るパラメータを推定することを特徴とする。
【００１３】
上記課題を解決するために、第４の観点に係るバッテリのパラメータ推定装置は、
前記ヒステリシスモデルに係るパラメータの推定と、前記過電圧モデルに係るパラメータの推定とを交互に繰り返すことを特徴とする。
【００１４】
上記課題を解決するために、第５の観点に係るバッテリのパラメータ推定装置は、
前記バッテリに入力される電流は、Ｍ系列信号であることを特徴とする。
【００１５】
上記課題を解決するために、第６の観点に係るバッテリのパラメータ推定装置は、
前記ヒステリシスモデルに係るパラメータは、前記ヒステリシスモデルにおける、電圧降下の最大範囲及び入力電流の特性により関連付けられるヒステリシス抵抗を含むことを特徴とする。
【００１６】
上記課題を解決するために、第７の観点に係るバッテリのパラメータ推定装置は、
前記ヒステリシスモデルに係るパラメータは、前記ヒステリシスモデルにおける、電圧降下の最大範囲及び電圧降下の速さの特性により関連付けられるヒステリシス容量を含むことを特徴とする。
【００１７】
上記課題を解決するために、第８の観点に係るバッテリのパラメータ推定装置は、
前記電圧降下の最大範囲をＭ（ｔ）、前記入力電流をｕ（ｔ）とした場合に、前記ヒステリシス抵抗Ｒh（ｔ）を式
【数１】
0006711981.tif 000002
で表すことを特徴とする。
【００１８】
上記課題を解決するために、第９の観点に係るバッテリのパラメータ推定装置は、
前記電圧降下の最大範囲をＭ（ｔ）、前記電圧降下の速さをΓ（ｔ）とした場合に、前記ヒステリシス容量Ｃh（ｔ）を式
【数２】
0006711981.tif 000003
で表すことを特徴とする。
【００１９】
上記課題を解決するために、第１０の観点に係るバッテリのパラメータ推定装置は、
前記ヒステリシス抵抗Ｒh（ｔ）を式
【数３】
0006711981.tif 000004
で表すことを特徴とする。
【００２０】
上記課題を解決するために、第１１の観点に係るバッテリのパラメータ推定装置は、
前記ヒステリシス抵抗Ｒh（ｔ）を式
【数４】
0006711981.tif 000005
で表すことを特徴とする。
【発明の効果】
【００２１】
第１の観点に係るバッテリのパラメータ推定装置によれば、ヒステリシスモデルに係るパラメータの推定精度と過電圧モデルに係るパラメータの推定精度とをそれぞれ向上させることができる。
【００２２】
第２の観点に係るバッテリのパラメータ推定装置によれば、パラメータの推定値をより効率よく真値に近づけることができる。
【００２３】
第３の観点に係るバッテリのパラメータ推定装置によれば、パラメータの推定値をより効率よく真値に近づけることができる。
【００２４】
第４の観点に係るバッテリのパラメータ推定装置によれば、パラメータの推定値をより速く真値に近づけることができる。
【００２５】
第５の観点に係るバッテリのパラメータ推定装置によれば、ヒステリシスモデルを線形モデルとして取り扱うことができ、且つ、広い周波数帯域を取り扱うことができる。
【００２６】
第６の観点に係るバッテリのパラメータ推定装置によれば、バッテリの等価回路にヒステリシス素子を増やすことなくヒステリシス現象を扱える。
【００２７】
第７の観点に係るバッテリのパラメータ推定装置によれば、バッテリの等価回路にヒステリシス素子を増やすことなくヒステリシス現象を扱える。
【００２８】
第８の観点に係るバッテリのパラメータ推定装置によれば、バッテリの等価回路にヒステリシス素子を増やすことなくヒステリシス現象を扱える。
【００２９】
第９の観点に係るバッテリのパラメータ推定装置によれば、バッテリの等価回路にヒステリシス素子を増やすことなくヒステリシス現象を扱える。
【００３０】
第１０の観点に係るバッテリのパラメータ推定装置によれば、より正確な推定値をより速く算出できる。
【００３１】
第１１の観点に係るバッテリのパラメータ推定装置によれば、より簡易なモデル構成とすることができ、ヒステリシス現象の取り扱いが容易になる。
【図面の簡単な説明】
【００３２】
【図１】実施形態１に係るパラメータ推定装置の構成例の機能ブロック図である。
【図２】ｎ次のＲＣ梯子回路（Ａ）フォスタ型（Ｂ）カウエル型である。
【図３】リン酸鉄リチウムイオン電池のＳＯＣ−ＯＣＶ特性の測定結果の例である。
【図４】Plettによるヒステリシスモデルの等価回路である。
【図５】可変抵抗Ｒhと可変容量Ｃhによって構成されるＲＣ並列回路である。
【図６】パラメータ感度プロットの一例である。
【図７】入力電流の振幅の大きさとパラメータの推定精度との関係の例である。
【図８】入出力データに基づくパラメータ推定方法の一例である。
【図９】入力信号に応じたパラメータに限定した推定手順の一例である。
【図１０】図８の方法を繰り返したときの推定値の変動の一例である。
【発明を実施するための形態】
【００３３】
（実施形態１）
以下、実施形態１に係るバッテリのパラメータ推定装置について、図面を参照しながら詳細に説明する。本実施形態に係るバッテリのパラメータ推定装置は、電気自動車やハイブリッド電気自動車などの車両に用いられる。このような車両には、車両を駆動する電気モータ、バッテリ、これらのコントローラなどが搭載され、電気モータへの電力の供給（放電）や制動時における電気モータからの制動エネルギの回生、地上充電設備からのバッテリへの電力回収（充電）が行われる。このような充放電電流のバッテリへの出入りがあると、バッテリの内部状態が変化していき、この内部状態をバッテリのパラメータ推定装置で推定しながらモニタしていくことで、バッテリの残量など必要な情報を収集している。
【００３４】
［機能ブロック］
図１に示されるように、パラメータ推定装置１は、パラメータ推定部１０と、記憶部２０と、充電率算出部３０と、電荷量算出部４０と、健全度算出部５０とを備える。パラメータ推定部１０は、電流センサ（入力検出部）２及び電圧センサ（出力検出部）３を介してバッテリ４に接続される。パラメータ推定部１０は、電源装置５に接続される。
【００３５】
電流センサ２は、バッテリ４に入力される電流（充放電電流）を検出して、検出した充放電電流をパラメータ推定部１０に対して出力する。電圧センサ３は、バッテリ４から出力される電圧（端子電圧）を検出して、検出した端子電圧をパラメータ推定部１０に対して出力する。
【００３６】
バッテリ４は、例えばリチャージャブル・バッテリ（二次電池）である。バッテリ４は、本実施の形態においてリチウム・イオン・バッテリであるものとして説明するが、他の種類のバッテリであってもよい。
【００３７】
電源装置５は、例えば電流源であり、パラメータ推定部１０からの指示に応じて、バッテリ４に所定の波形を有する電流を入力する。
【００３８】
パラメータ推定部１０は、出力予測部１１と、誤差算出部１２と、パラメータ更新部１３とを備える。パラメータ推定部１０は、例えばプロセッサ又はマイクロコンピュータ等で構成される。出力予測部１１、誤差算出部１２及びパラメータ更新部１３は、それぞれ別個のチップとして構成されてもよいし、１チップにまとめて構成されてもよい。
【００３９】
記憶部２０は、例えば半導体メモリ又は磁気記憶装置等で構成される。記憶部２０は、バッテリ等価回路モデル２１を格納する。記憶部２０は、パラメータ推定装置１が取り扱うデータ又は情報等を格納してもよい。
【００４０】
バッテリ等価回路モデル２１は、バッテリ４の内部状態を等価回路で表したモデル（等価回路モデル）である。バッテリ等価回路モデル２１は、パラメータとして、バッテリ４の開回路電圧（ＯＣＶ：Ｏｐｅｎ  Ｃｉｒｃｕｉｔ  Ｖｏｌｔａｇｅ）、並びに、バッテリ４の内部抵抗及び容量を有する。バッテリ等価回路モデル２１のパラメータは、記憶部２０に格納される。
【００４１】
出力予測部１１は、電流センサ２からバッテリ４の充放電電流（ｉ）を取得し、記憶部２０からバッテリ等価回路モデル２１のパラメータを取得する。出力予測部１１は、バッテリ４の充放電電流とバッテリ等価回路モデル２１のパラメータとに基づいて、バッテリ４の端子電圧の予測値を算出する。
【００４２】
誤差算出部１２は、電圧センサ３からバッテリ４の端子電圧（ｖ）を取得し、出力予測部１１からバッテリ４の端子電圧の予測値を取得する。誤差算出部１２は、電圧センサ３から取得した端子電圧と、出力予測部１１から取得した端子電圧の予測値との誤差を算出する。
【００４３】
パラメータ更新部１３は、誤差算出部１２から端子電圧の誤差を取得する。パラメータ更新部１３は、端子電圧の誤差に基づいて、記憶部２０に格納されているバッテリ等価回路モデル２１のパラメータを更新する。
【００４４】
バッテリ等価回路モデル２１のパラメータは、バッテリ４の内部状態を表す。つまりパラメータ推定部１０は、出力予測部１１、誤差算出部１２及びパラメータ更新部１３を動作させることにより、バッテリ４の内部状態を推定する。
【００４５】
充電率算出部３０は、パラメータ推定部１０で推定されたバッテリ４の内部状態に基づいて、バッテリ４の充電率（ＳＯＣ：Ｓｔａｔｅ  Ｏｆ  Ｃｈａｒｇｅ）を算出する。ＳＯＣは、バッテリ４の充電容量に対する充電量の比である。バッテリ４のＯＣＶとＳＯＣとの関係は、バッテリ４の温度又は劣化度合い等に影響されにくい。よって、バッテリ４のＳＯＣは、予め実験等によって求められたバッテリ４のＯＣＶとＳＯＣとの関係、及び、パラメータ推定部１０で推定されたＯＣＶに基づいて算出されうる。本実施形態において、充電率算出部３０は、上述の方法によりバッテリ４のＳＯＣを算出する。充電率算出部３０で算出されるＳＯＣは、バッテリ４のバッテリ・マネジメントに利用される。
【００４６】
電荷量算出部４０は、電流センサ２からバッテリ４の充放電電流値（ｉ）を取得する。電荷量算出部４０は、充放電電流値（ｉ）を逐次積算していくことでバッテリ４に出入りする電荷量を算出する。電荷量算出部４０は、充放電電流値（ｉ）から算出した電荷量を、逐次積算演算を開始する前にバッテリ４に充電されている電荷量（残存電荷量ともいう）から減算することで、バッテリ４に現在充電されている電荷量（現存電荷量ともいう）を算出する。
【００４７】
健全度算出部５０は、バッテリ４の健全度（ＳＯＨ：Ｓｔａｔｅ  Ｏｆ  Ｈｅａｌｔｈ）を算出する。ＳＯＨは、バッテリ４の現在の満充電容量（ＦＣＣ：Ｆｕｌｌ  Ｃｈａｒｇｅ  Ｃａｐａｃｉｔｙ）が、バッテリ４の新品時のＦＣＣに対して何％であるかを示す。バッテリ４の新品時のＦＣＣは、既知の値である。バッテリ４の新品時のＦＣＣは、健全度算出部５０に格納されてもよいし、記憶部２０に格納されてもよい。健全度算出部５０は、充電率算出部３０からバッテリ４のＳＯＣを取得し、電荷量算出部４０からバッテリ４の現存電荷量を取得する。健全度算出部５０は、バッテリ４のＳＯＣと現存電荷量とに基づいて、バッテリ４の現在のＦＣＣを算出する。健全度算出部５０は、バッテリ４の新品時のＦＣＣとバッテリ４の現在のＦＣＣとに基づいて、バッテリ４のＳＯＨを算出する。
【００４８】
健全度算出部５０は、バッテリ４のＳＯＨと、バッテリ４の現存電荷量と、バッテリ４のＯＣＶとの関係に基づいて、バッテリ４のＳＯＨを算出してもよい。健全度算出部５０は、例えば、本出願人の出願による特開２０１２−５７９５６号公報に開示されているような、所定幅で区分けしたＳＯＨごとに現存電荷量とＯＣＶとの関係を表わす特性表に基づいて、ＳＯＨを算出してもよい。この場合、健全度算出部５０は、記憶部２０に格納されたＯＣＶの推定値と、電荷量算出部４０で算出された現存電荷量とが特性表のいずれのＳＯＨの範囲に入るのか判定し、バッテリ４のＳＯＨを算出する。
【００４９】
［バッテリ等価回路モデル］
本実施形態においては、バッテリ４の内部状態を推定するためのバッテリ等価回路モデル２１として、ＯＣＶ、過電圧パラメータ及びヒステリシスパラメータの３つの要素を含むモデルを考える。バッテリ等価回路モデル２１は、例えば、ＯＣＶに対応する電圧源と、過電圧パラメータ及びヒステリシスパラメータにそれぞれ対応するインピーダンスとが直列に接続する開回路である。開回路に流れる電流は、バッテリ４の充放電電流（ｉ）に対応する。開回路の両端の電圧は、バッテリ４の端子電圧（ｖ）に対応する。
【００５０】
＜過電圧モデル＞
過電圧パラメータは、過電圧に関するモデル（過電圧モデル）を表すパラメータである。過電圧は、バッテリ４のＯＣＶと、バッテリ４に充放電電流が入力される際の端子電圧との差であり、バッテリ４の内部インピーダンスに起因する。バッテリ４の内部インピーダンスは、バッテリ４の電解液の電気抵抗、及び、バッテリ４の電極反応に起因するインピーダンスを含む。
【００５１】
バッテリ４の電極反応は、電解液と活物質との界面における電荷移動過程と、電解液又は活物質におけるイオンの拡散過程とを含む。例えば、リチウム・イオン・バッテリ等の物理過程（ｎｏｎ−Ｆａｒａｄａｉｃ  ｐｒｏｃｅｓｓ）バッテリ、つまり、拡散過程が支配的なバッテリにおいては、拡散過程に起因するインピーダンスであるワールブルグインピーダンスの影響が支配的となる。
【００５２】
本実施形態においては、過電圧パラメータに対応するインピーダンスとして、電解液等の抵抗を含む内部抵抗（Ｒ0）と、ワールブルグインピーダンス（Ｚw）とが直列に接続されるインピーダンスを想定する。
【００５３】
ＳＯＣは、バッテリ４の充放電電流（ｉ）と、バッテリ４の新品時の満充電容量（ＦＣＣ0）と、バッテリ４のＳＯＨとを用いて、式（１）で表される。
【数５】
0006711981.tif 000006
【００５４】
また、ワールブルグインピーダンス（Ｚw）の伝達関数は、式（２）により表される。
【数６】
0006711981.tif 000007
ただし、ｓはラプラス演算子であり、ＲdはＺw（ｓ）の低周波極限（ω→０）で算出される拡散抵抗であり、τdは、拡散反応の速度に対応する拡散時定数である。また、拡散抵抗及び拡散時定数を用いた式（３）により拡散容量（Ｃd）が定義される。
【数７】
0006711981.tif 000008
【００５５】
式（２）において、ワールブルグインピーダンス（Ｚw）は、ラプラス演算子（ｓ）の平方根を含む。ワールブルグインピーダンス（Ｚw）を時間領域へ変換しやすくするために、ワールブルグインピーダンス（Ｚw）を近似する。
【００５６】
＜＜無限級数展開＞＞
ワールブルグインピーダンス（Ｚw）は、例えば式（４）に示すように、無限級数の和として表されうる。
【数８】
0006711981.tif 000009
ただし、
【数９】
0006711981.tif 000010
である。上述の式（４）は、図２（Ａ）に示されるような、抵抗（Ｒ1〜Ｒn）とコンデンサ（Ｃ1〜Ｃn）との並列回路が直列にｎ個接続されたｎ次フォスタ型回路に対応する。式（５）（６）から明らかなように、ワールブルグインピーダンス（Ｚw）を近似したｎ次のフォスタ型等価回路モデルによれば、拡散容量（Ｃd）及び拡散抵抗（Ｒd）を用いて、バッテリ等価回路モデル２１の他のパラメータ（Ｒn又はＣn等）を算出可能である。
【００５７】
＜＜連分数展開＞＞
ワールブルグインピーダンス（Ｚw）は、式（７）に示すように、連分数展開によっても表されうる。
【数１０】
0006711981.tif 000011
ただし、
【数１１】
0006711981.tif 000012
である。上述の式（７）は、図２（Ｂ）に示されるような、並列に接続されたｎ個の抵抗（Ｒ1〜Ｒn）のそれぞれが、直列に接続されたｎ個のコンデンサ（Ｃ1〜Ｃn）の間に接続されたｎ次カウエル型回路に対応する。式（８）（９）から明らかなように、ワールブルグインピーダンス（Ｚw）を近似したｎ次のカウエル型等価回路モデルによれば、拡散容量（Ｃd）及び拡散抵抗（Ｒd）を用いて、バッテリ等価回路モデル２１の他のパラメータ（Ｒn又はＣn等）を算出可能である。
【００５８】
＜ヒステリシスモデル＞
ヒステリシスパラメータは、バッテリ４のヒステリシス現象に関するモデル（ヒステリシスモデル）を表すパラメータである。バッテリ４において、充電後と放電後とでＳＯＣ−ＯＣＶ特性が異なるヒステリシス現象が発生することがある。ヒステリシス現象は電極の材料により発生し、特にリン酸リチウムを用いた場合にはヒステリシス現象の影響が大きく出る。
【００５９】
図３はリン酸鉄リチウムイオン電池のＳＯＣ−ＯＣＶ特性の測定結果の一例である。図３（Ａ）によれば充電時の特性と放電時の特性との間でＯＣＶの差が生じていることが分かる。また図３（Ａ）の破線囲み部を拡大した図３（Ｂ）において、ＳＯＣが約３０％の時点で放電するようにしてもヒステリシス特性を示すことが分かる。ヒステリシス現象が発生する場合、上述の過電圧モデルだけでは、バッテリ４のＳＯＣ−ＯＣＶ特性を正確に表すことができない。
【００６０】
このようなヒステリシス現象を表すモデルの一つであるPlettによるヒステリシスモデルは、図４に示される等価回路で表される。ここでＶHがヒステリシス電圧を表す素子である。このヒステリシスモデルは、以下の式（１０）で表される。
【数１２】
0006711981.tif 000013
ここで、ｖh（ｔ）はヒステリシス電圧、Γ（ｔ）はヒステリシスモデルの電圧降下の速さ（ＳＯＣ−ＯＣＶ曲線の傾きに相当）、Ｍ（ｔ）はヒステリシスモデルの電圧降下の最大範囲、ｕ（ｔ）は入力電流を表すパラメータである。
【００６１】
本来、バッテリ４のヒステリシスはバッテリ内部の電気化学反応の結果としてあらわれてくるもので、バッテリ４内部の電荷移動過程やイオンの拡散過程と密接な関係が有る。しかしながらPlettによるヒステリシスモデルは、電荷移動過程やイオンの拡散過程とは独立した反応を追加してｖh（ｔ）を表している。そのため、式（１０）に基づいてバッテリ４の内部状態を推定する場合には、ＲＣ並列回路の抵抗及び容量の推定に加え、ヒステリシス電圧を表すΓ（ｔ）及びＭ（ｔ）を推定する必要がある。つまりヒステリシスを考えないモデルによる推定と比較して、推定すべきパラメータが２つ（ΓとＭ）増加する。
【００６２】
ここで、本実施形態においては、かかる独立した反応を電荷移動過程やイオンの拡散過程に対応したモデルに統合する。ヒステリシスモデルを表す式（１０）は、以下の式（１１）（１２）を用いて、以下の式（１３）のように変形される。
【数１３】
0006711981.tif 000014
【００６３】
式（１３）は、図５に示されるような可変抵抗（Ｒh）と可変容量（Ｃh）によって構成されるＲＣ並列回路を表す式と同等であると解釈できる。特に式（１１）で表されるように、モデルの抵抗は電流の大きさによって可変となる可変抵抗であることが特徴である。
【００６４】
以上説明してきたように、ヒステリシスモデルを式（１１）〜（１３）で表すことができる。このようにヒステリシスモデルを表すことで、ヒステリシス素子を増やすことなく、ヒステリシス現象を取り扱うことができる。
【００６５】
ヒステリシスモデルに係るパラメータは、可変抵抗（Ｒh）及び可変容量（Ｃh）である。可変抵抗（Ｒh）は、ヒステリシス抵抗ともいう。ヒステリシス抵抗は、ヒステリシスモデルの電圧降下の最大範囲（Ｍ（ｔ））及び入力電流（ｕ（ｔ））の特性により関連付けられる。可変容量（Ｃh）は、ヒステリシス容量ともいう。ヒステリシス容量は、ヒステリシスモデルの電圧降下の最大範囲（Ｍ（ｔ））及びヒステリシスモデルの電圧降下の速さ（Γ（ｔ））の特性により関連付けられる。
【００６６】
［連続時間システム同定法］
本実施形態においては、バッテリ等価回路モデル２１で表されるシステムに連続時間システム同定法を適用して、バッテリ等価回路モデル２１のパラメータが推定される。連続時間システム同定法による推定の対象となるパラメータは、過電圧モデルのＲ0、Ｒd及びτd（Ｃd）、並びに、ヒステリシスモデルのＲh及びＣhである。
【００６７】
連続時間システム同定法によれば、バッテリ等価回路モデル２１で表されるシステムの入出力データの関係を表す式に基づいて、各パラメータが推定される。以下、連続時間システム同定法によるパラメータ推定方法について説明する。
【００６８】
＜伝達関数＞
ｕ（ｔ）が入力された場合にｙ（ｔ）が出力される線形時不変システムを考える。該システムの伝達関数がＧ（ｐ，θ）と表される場合、ｕ（ｔ）とｙ（ｔ）との関係は、以下の式（１４）により表される。
【数１４】
0006711981.tif 000015
ｐは、微分演算子（ｄ／ｄｔ）である。θは、ｍ個のパラメータの集合（θ＝［θ1・・・θm］）である。θは、ｍ次元の列ベクトルとしても表される。
【００６９】
Ｇ（ｐ，θ）は、パラメータ（θ）に関して線形な関数であるＮ（ｐ，θ）及びＤ（ｐ，θ）を用いて、以下の式（１５）のようにも表される。
【数１５】
0006711981.tif 000016
Ｎ（ｐ，θ）は、ｍ次元行ベクトルとして表されるｎ（ｐ）＝［ｎ1（ｐ）・・・ｎm（ｐ）］に対してｍ次元列ベクトルとして表されるθを掛けて算出される内積である。Ｄ（ｐ，θ）は、ｍ次元行ベクトルとして表されるｄ（ｐ）＝［ｄ1（ｐ）・・・ｄm（ｐ）］に対してｍ次元列ベクトルとして表されるθを掛けて算出される内積である。
【００７０】
＜パラメータ誤差＞
伝達関数（Ｇ）に含まれるパラメータ（θ）は、連続時間システム同定法により推定される。パラメータの推定値をθ＾とした場合、システムの出力の推定値（ｙ＾（ｔ））は、以下の式（１６）のように表される。
【数１６】
0006711981.tif 000017
なお、数式内で用いられる、文字の上に「＾」が付される記号は、明細書における表記の制約上、「ｙ＾」のように表すことがある。
【００７１】
連続時間システム同定法においては、システムの出力（ｙ（ｔ））と、出力の推定値（ｙ＾（ｔ））との誤差が最小となるようにパラメータ（θ）が推定される。ｙ＾（ｔ）とｙ（ｔ）との誤差は、推定誤差信号（ｅ（ｔ））として、以下の式（１７）により定義される。
【数１７】
0006711981.tif 000018
【００７２】
ここでパラメータの真値（θ0）が存在してｙ（ｔ）＝Ｇ（ｐ，θ0）ｕ（ｔ）が成り立つ場合、式（１７）は、式（１４）〜（１６）を用いて、以下の式（１８）のように変形される。
【数１８】
0006711981.tif 000019
【００７３】
さらに式（１８）は、Ｎ（ｐ，θ）＝ｎ（ｐ）・θ、及び、Ｄ（ｐ，θ）＝ｄ（ｐ）・θの関係を用いて、以下の式（１９）のように変形される。
【数１９】
0006711981.tif 000020
【００７４】
式（１９）に含まれるθ0−θ＾は、パラメータ誤差を表す。パラメータ誤差をεとすると、ε＝θ0−θ＾となる。この場合、式（１９）は、以下の式（２０）のように変形される。
【数２０】
0006711981.tif 000021
【００７５】
ここで、部分空間の基底ベクトル（Ｖ）を式（２１）のように定義する。
【数２１】
0006711981.tif 000022
【００７６】
式（２１）のように部分空間が定義される場合、式（２０）によれば、パラメータ誤差（ε）は、推定誤差信号（ｅ（ｔ））が基底ベクトル（Ｖ）を有する部分空間に射影されたときの射影ベクトルとなる。
【００７７】
式（２１）に示される基底ベクトル（Ｖ）に含まれるパラメータの真値（θ0）を、パラメータの推定値（θ＾）に置き換える場合、パラメータ誤差（ε）は、以下の式（２２）のように表される。
【数２２】
0006711981.tif 000023
Ｉは単位行列である。ΔはＧ（ｐ，θ0）とＧ（ｐ，θ＾）との差に起因する項である。
【００７８】
＜パラメータ推定＞
式（２０）〜（２２）により示される推定誤差信号（ｅ（ｔ））が最小化されるようにパラメータの推定値（θ＾）を算出することにより、パラメータの推定値を真値に近づけることができる。推定誤差信号（ｅ（ｔ））は、例えばガウス・ニュートン法を用いた反復計算を行うことによって最小化されうる。
【００７９】
ガウス・ニュートン法を用いて推定誤差信号を最小化するようにパラメータの推定値を算出するための反復計算は、例えば以下の式（２３）のように表される。
【数２３】
0006711981.tif 000024
上添え字のｋは、ｋ回目の反復計算で得られる値であることを示している。つまりθ＾kは、ｋ回目の反復計算で得られるパラメータの推定値である。αは、反復計算によるパラメータ推定のアルゴリズムの安定性を決めるステップ幅パラメータであり、０＜α＜１である。
【００８０】
式（２３）は、式（２２）を用いて、以下の式（２４）に示されるように変形される。
【数２４】
0006711981.tif 000025
【００８１】
式（２４）で表されるシステムの平衡点は、θ＾k＝θ0である。該システムの伝達関数をインパルス応答のｚ変換の形式で表した場合、Ｈ（ｚ）＝［−α／｛ｚ−（１−α）｝］×Ｉと表される。０＜α＜１である場合、Ｈ（ｚ）の極（１−α）が単位円の内側となるので、該システムは安定であるといえる。つまり、反復計算においてθ＾kが発散してしまうことなく、真値であるθ0に収束するといえる。
【００８２】
＜バッテリ等価回路モデルへの適用＞
バッテリ４から取得される入出力データがサンプリング時刻をｔi（ｉ＝０、１、・・・、Ｎ）として離散的に取得される場合を考える。この場合における基底ベクトル（ＶD）及び推定誤差信号（ｅD）はそれぞれ、以下の式（２５）及び（２６）のように表される。
【数２５】
0006711981.tif 000026
【００８３】
パラメータ誤差（ε）は、式（２０）で表されるｅ（ｔ）及び式（２１）で表されるＶをそれぞれｅD及びＶDに置換して変形することで、例えば式（２７）のように表される。
【数２６】
0006711981.tif 000027
ここで、Ｔは転置行列を意味する。
【００８４】
バッテリ等価回路モデル２１のパラメータ（θ）は、ＯＣＶ、過電圧パラメータ及びヒステリシスパラメータであり、以下の式（２８）で表される。
【数２７】
0006711981.tif 000028
【００８５】
バッテリ等価回路モデル２１のパラメータの推定に、式（２４）で表されるシステムを適用するために、パラメータの各要素と基底ベクトルとを対応づける。ここで、式（１６）の両辺をパラメータ（θ）の要素であるθjで偏微分して、以下の式（２９）が導かれる。
【数２８】
0006711981.tif 000029
【００８６】
さらに、式（２９）と式（２１）とに基づいて、以下の式（３０）が導かれる。
【数２９】
0006711981.tif 000030
【００８７】
式（３０）は、基底ベクトル（Ｖ）の要素（Ｖj）が、出力の推定値（ｙ＾）をパラメータの各要素（θj）で偏微分したものとして表されることを示している。式（３０）によって、バッテリ等価回路モデル２１のパラメータと式（２４）で表されるシステムの基底ベクトルとが対応づけられる。式（３０）をパラメータの感度関数ともいう。
【００８８】
＜パラメータ感度プロット＞
式（３０）で示されるパラメータの感度関数から、パラメータの感度が定義される。ラプラス演算子（ｓ）で表される伝達関数（Ｇ（ｓ））を用いて、パラメータの各要素（θj）に対応する正規化パラメータ感度（以下、単に感度又はパラメータ感度ともいう）が式（３１）のように定義される。
【数３０】
0006711981.tif 000031
【００８９】
パラメータが式（２８）に示されるものである場合の正規化パラメータ感度は式（３２）のように表される。
【数３１】
0006711981.tif 000032
【００９０】
正規化パラメータ感度は、周波数特性を有する。図６に例示されるように、横軸に周波数をとり、縦軸に正規化パラメータ感度をとったグラフを、パラメータ感度プロットという。
【００９１】
バッテリ等価回路モデル２１のパラメータを推定する場合、パラメータの各要素を同時に推定してもよい。この際、感度が高い要素ほど、精度よく推定されうる。例えば図６（Ａ）において、Ｒ0は、広い周波数帯域において高い感度を有している。一方でＣhは、他の要素と比較して感度が低くなっている。この場合、Ｒ0はＣhよりも精度よく推定されうる。
【００９２】
図６（Ａ）において、τdは、表示されている周波数帯域の大部分で感度が低くなっているが、Ωτdで示される周波数帯域においては、比較的高い感度を有している。この場合、τdは、入出力データのうちΩτdで示される周波数帯域のデータに基づいて、比較的精度よく推定されうる。
【００９３】
＜入力信号の使い分け＞
以上説明してきたように、連続時間システム同定法を用いて、バッテリ等価回路モデル２１のパラメータを推定できる。パラメータの推定の精度を向上させるためには、正規化パラメータ感度が高い周波数帯域のデータを用いる方がよい。よって、パラメータの各要素を全て同時に推定するのではなく、感度が高いパラメータに限定して推定してもよい。
【００９４】
ここで、バッテリ４に入力される電流の振幅が比較的大きい場合と、比較的小さい場合とで、それぞれ異なるパラメータ感度プロットが得られる。図６（Ａ）には、入力電流の振幅が比較的大きい場合のパラメータ感度プロットが示される。一方で図６（Ｂ）には、入力電流の振幅が比較的小さい場合のパラメータ感度プロットが示される。
【００９５】
入力電流の振幅が比較的大きい場合、図６（Ａ）に示されるように、表示されている周波数の全帯域にわたって、Ｒ0及びＲdのパラメータ感度は、比較的高い。よって、Ｒ0及びＲdについては比較的精度よく推定されうる。τdのパラメータ感度は、Ωτdで示される周波数帯域において、比較的高くなっている。よって、τdについては、Ωτdで示される周波数帯域のデータに基づいて、比較的精度よく推定されうる。一方で、Ｒh及びＣhのパラメータ感度は、表示されている周波数の全帯域にわたって比較的低い。よって、Ｒh及びＣhについては、推定精度を向上させることが難しい。
【００９６】
入力電流の振幅が比較的小さい場合、図６（Ｂ）に示されるように、表示されている周波数の帯域の一部において、Ｒh及びＣhのパラメータ感度は、比較的高い。例えば、Ｒhのパラメータ感度は、ΩRhで示される周波数帯域において、比較的高くなっている。また、Ｃhのパラメータ感度は、ΩChで示される周波数帯域において、比較的高くなっている。よって、Ｒh及びＣhについては、それぞれΩRh及びΩChで示される周波数帯域のデータに基づいて、比較的精度よく推定されうる。一方で、Ｒd及びτdのパラメータ感度は、表示されている周波数の全帯域にわたって比較的低い。よって、Ｒd及びτdについては、推定精度を向上させることが難しい。Ｒ0のパラメータ感度は、ΩR0で示される周波数帯域において、比較的高くなっている。よって、Ｒ0については、ΩR0で示される周波数帯域のデータに基づいて、比較的精度よく推定されうる。
【００９７】
以上のように、入力電流の振幅が比較的大きい場合と比較的小さい場合とで、それぞれ異なるパラメータが比較的精度よく推定されうる（図７参照）。入力電流の振幅が比較的小さい場合（図７の２行目）、ヒステリシスパラメータ（Ｒh及びＣh）の推定精度は比較的高い。一方で、過電圧パラメータ（Ｒd及びτd）の推定精度は比較的低い。入力電流の振幅が比較的大きい場合（図７の３行目）、ヒステリシスパラメータ（Ｒh及びＣh）の推定精度は比較的低い。一方で、過電圧パラメータ（Ｒd及びτd）の推定精度は比較的高い。
【００９８】
［パラメータ推定方法］
以上説明してきたように、バッテリ４の内部状態は、連続時間システム同定法を用いて推定されうる。本実施形態に係るパラメータ推定装置１は、パラメータの各要素を同時に推定するのではなく、入力信号に応じた感度が高い要素を選択的に推定する。このようにすることで、パラメータの推定精度を向上させることができる。以下、本実施形態に係るパラメータ推定方法について説明する。
【００９９】
＜入力電流＞
本実施形態において、バッテリ等価回路モデル２１は、ヒステリシスモデルを含む。ヒステリシスモデルは、電流の大きさによって可変となる可変抵抗を含む。バッテリ等価回路モデル２１に入力される電流の大きさ（絶対値）が一定である場合、式（１１）で表されるヒステリシスモデルの抵抗（Ｒh）は時不変となる。この場合、ヒステリシスモデルは、線形時不変のモデルとして取り扱うことができるものとなる。
【０１００】
入力される電流は、例えばＭ系列信号（Ｍａｘｉｍｕｍ  Ｌｅｎｇｔｈ  Ｓｅｑｕｅｎｃｅ）である。Ｍ系列信号は、振幅が一定の周期信号であり、且つ、１周期内では擬似的な乱数としての性質を示すものである。Ｍ系列信号は、その性質上、広帯域の周波数成分を含む。
【０１０１】
一般的なＭ系列信号は、０と１との２値をとるものである。本実施形態において入力されるＭ系列信号は、絶対値が同一である正の値と負の値との２値をとるものとする。つまりＭ系列信号が入力される場合、Ｒhは時不変となる。
【０１０２】
入力される電流は、Ｍ系列信号に限られず、振幅が一定の矩形波の信号であってもよい。入力される電流は、広帯域の周波数成分を含む信号であってもよいし、パラメータの推定に用いられる一部の周波数成分を含む信号であってもよい。
【０１０３】
＜入出力データに基づくパラメータの推定＞
パラメータ推定部１０は、図８のフローチャートを実行して、１組の入出力データ（ｕ（ｔ）及びｙ（ｔ））に基づく反復計算によってパラメータの推定値を算出する。
【０１０４】
まず、パラメータ推定部１０は、電源装置５を用いて、バッテリ４に対してＭ系列信号の電流を入力する（ステップＳ１）。入力される電流は、Ｍ系列信号に限られず、他の信号であってもよい。
【０１０５】
続いて、パラメータ推定部１０は、電流センサ２から入力電流（ｕ（ｔ））を取得し、電圧センサ３から端子電圧（ｙ（ｔ））を取得する（ステップＳ２）。パラメータ推定部１０は、ｕ（ｔ）及びｙ（ｔ）を所定のサンプリング時刻で取得してもよい。本実施形態においては、所定のサンプリング時刻をｔi（ｉ＝０、１、・・・、Ｎ）とする。パラメータ推定部１０は、取得したｕ（ｔ）及びｙ（ｔ）のデータを記憶部２０に格納してもよい。
【０１０６】
続いて、パラメータ推定部１０は、パラメータの推定値の初期値（θ＾0）を設定する（ステップＳ３）。ガウス・ニュートン法は、勾配法の一種であり、初期値の選択によっては局所最適解に陥るおそれもある。よって初期値は、例えばバッテリ４に関する既知の事前情報に基づいて設定される。
【０１０７】
続いて、パラメータ推定部１０は、反復計算の回数のカウンタ（ｋ）を０に設定する（ステップＳ４）。つまり、ステップＳ３で設定された初期値（θ＾0）は、０回目の推定値である。
【０１０８】
続いて、パラメータ推定部１０は、ｋ回目の推定値（θ＾k）に基づいて、基底（ＶD）と推定誤差信号（ｅD）とを設定する（ステップＳ５）。基底及び推定誤差信号は、上述の式（２５）、（２６）、（２８）及び（３０）に基づいて設定される。
【０１０９】
続いて、パラメータ推定部１０は、（ｋ＋１）回目の反復計算（ｋ＝０の場合、１回目の反復計算）を行い、パラメータの推定値（θ＾k+1）を算出する（ステップＳ６）。この算出は、式（２２）〜（２４）に基づいて行われる。パラメータ推定部１０は、算出した推定値を記憶部２０に格納してもよい。
【０１１０】
続いて、パラメータ推定部１０は、パラメータの推定値が収束したか判定する（ステップＳ７）。収束の判定は、例えば以下の式（３３）を用いて行われる。
【数３２】
0006711981.tif 000033
ζは、収束判定値であり、適宜設定されうる。パラメータ推定部１０は、式（３３）が満たされる場合に推定値が収束したと判定する。
【０１１１】
パラメータの推定値が収束していない場合（ステップＳ７：ＮＯ）、パラメータ推定部１０は、反復計算の回数のカウンタ（ｋ）に１を加える（ステップＳ８）。その後、パラメータ推定部１０は、ステップＳ５に戻って、パラメータ推定の反復計算を続ける。
【０１１２】
パラメータの推定値が収束した場合（ステップＳ７：ＹＥＳ）、パラメータ推定部１０は、図８のフローチャートの処理を終了する。図８に示される反復計算による推定結果は、ステップＳ６で算出された（ｋ＋１）回目の推定値（θ＾k+1）である。
【０１１３】
＜入力信号に応じたパラメータの推定＞
図９は、入力信号に応じたパラメータに限定して推定するパラメータ推定方法の一例である。パラメータ推定部１０は、図９のフローチャートに沿って、振幅が比較的大きい入力信号と振幅が比較的小さい入力信号とをそれぞれ使い分けて、各入力信号に対応するパラメータを推定する。
【０１１４】
まず、パラメータ推定部１０は、電源装置５を用いて、振幅が比較的大きいＭ系列信号をバッテリ４に対して入力する（ステップＳ２１）。入力される電流は、Ｍ系列信号に限られず、他の信号であってもよい。ステップＳ２１で入力される振幅が比較的大きい電流のことを、第１の電流ともいう。第１の電流は、第１の振幅を有するものとする。
【０１１５】
続いて、パラメータ推定部１０は、電流センサ２から入力電流（ｕ１（ｔ））を取得し、電圧センサ３から端子電圧（ｙ１（ｔ））を取得する（ステップＳ２２）。パラメータ推定部１０は、ｕ１（ｔ）及びｙ１（ｔ）を所定のサンプリング時刻で取得してもよい。ステップＳ２２では、図８のステップＳ２と同様に、所定のサンプリング時刻をｔi（ｉ＝０、１、・・・、Ｎ）とする。パラメータ推定部１０は、取得したｕ１（ｔ）及びｙ１（ｔ）のデータを記憶部２０に格納してもよい。
【０１１６】
続いて、パラメータ推定部１０は、ステップＳ２２で取得したｕ１（ｔ）及びｙ１（ｔ）に基づいて、図８のステップＳ３〜Ｓ８の処理を実行して、バッテリ等価回路モデル２１のパラメータを推定する（ステップＳ２３）。ｕ１（ｔ）及びｙ１（ｔ）は、バッテリ４に入力される電流の振幅が比較的大きい場合に対応する入出力データである。この場合、パラメータ感度プロットは、図６（Ａ）に例示されるものと同様の傾向を有するものとなる。図６（Ａ）によれば、過電圧パラメータ（Ｒ0、Ｒd、τd）の感度が比較的高い。よって、パラメータの推定精度が高められるように、パラメータ推定部１０は、図８のステップＳ６で算出するパラメータを過電圧パラメータ（Ｒ0、Ｒd、τd）に限定する。一方で、パラメータ推定部１０は、ヒステリシスパラメータ（Ｒh、Ｃh）を図８のステップＳ３で設定された初期値のまま固定する。
【０１１７】
続いて、パラメータ推定部１０は、電源装置５を用いて、振幅が比較的小さいＭ系列信号をバッテリ４に対して入力する（ステップＳ２４）。入力される電流は、ステップＳ２１と同様にＭ系列信号に限られず、他の信号であってもよい。ステップＳ２４で入力される振幅が比較的小さい電流のことを、第２の電流ともいう。第２の電流は、第２の振幅を有するものとする。第２の振幅は、上述の第１の振幅よりも小さい。
【０１１８】
続いて、パラメータ推定部１０は、電流センサ２から入力電流（ｕ２（ｔ））を取得し、電圧センサ３から端子電圧（ｙ２（ｔ））を取得する（ステップＳ２５）。パラメータ推定部１０は、ｕ２（ｔ）及びｙ２（ｔ）を所定のサンプリング時刻で取得してもよい。ステップＳ２５では、ステップＳ２２と同様に、所定のサンプリング時刻をｔi（ｉ＝０、１、・・・、Ｎ）とする。パラメータ推定部１０は、取得したｕ２（ｔ）及びｙ２（ｔ）のデータを記憶部２０に格納してもよい。
【０１１９】
続いて、パラメータ推定部１０は、ステップＳ２５で取得したｕ２（ｔ）及びｙ２（ｔ）に基づいて、図８のステップＳ３〜Ｓ８の処理を実行して、バッテリ等価回路モデル２１のパラメータを推定する（ステップＳ２６）。ｕ２（ｔ）及びｙ２（ｔ）は、バッテリ４に入力される電流の振幅が比較的小さい場合に対応する入出力データである。この場合、パラメータ感度プロットは、図６（Ｂ）に例示されるものと同様の傾向を有するものとなる。図６（Ｂ）によれば、ヒステリシスパラメータ（Ｒh、Ｃh）の感度が比較的高い。よって、パラメータの推定精度が高められるように、パラメータ推定部１０は、図８のステップＳ６で算出するパラメータをヒステリシスパラメータ（Ｒh、Ｃh）に限定する。一方で、パラメータ推定部１０は、過電圧パラメータ（Ｒ0、Ｒd、τd）をステップＳ２３で算出された推定値のまま固定する。
【０１２０】
続いて、パラメータ推定部１０は、ステップＳ２２で取得したｕ１（ｔ）及びｙ１（ｔ）に基づいて、図８のステップＳ３〜Ｓ８の処理を実行して、バッテリ等価回路モデル２１の過電圧パラメータを推定する（ステップＳ２７）。この際、図８のステップＳ３で設定される初期値として、ステップＳ２３で算出された過電圧パラメータの推定値が用いられる。ヒステリシスパラメータ（Ｒh、Ｃh）は、ステップＳ２６で算出された推定値のまま固定される。
【０１２１】
続いて、パラメータ推定部１０は、ステップＳ２６及びＳ２７においてそれぞれ算出された過電圧パラメータ及びヒステリシスパラメータの推定値が収束したか判定する（ステップＳ２８）。パラメータ推定部１０は、例えば、ステップＳ２６及びＳ２７の処理の前後の推定値の差が、所定の値未満となった場合に、推定値が収束したと判定する。収束判定の方法は、これに限られない。例えば、推定値の変化から、最終的な収束値を予測して、収束したものとみなしてもよい。
【０１２２】
パラメータの推定値が収束していない場合（ステップＳ２８：ＮＯ）、パラメータ推定部１０は、ステップＳ２６に戻って、パラメータ推定を続ける。この場合、過電圧パラメータの推定とヒステリシスパラメータの推定とが交互に行われることとなる。パラメータの一部が収束した場合、収束したパラメータを除いて推定を続けてもよいし、収束したパラメータも含めて推定を続けてもよい。
【０１２３】
パラメータの推定値が収束した場合（ステップＳ２８：ＹＥＳ）、パラメータ推定部１０は、図９のフローチャートの処理を終了する。図９に示される処理によるパラメータの推定結果は、ステップＳ２６及びＳ２７で算出された推定値である。
【０１２４】
図９のフローチャートの処理において、パラメータ推定部１０は、ステップＳ２６における過電圧パラメータの推定を繰り返す際に、既に取得した入出力データ（ｕ１（ｔ）及びｙ１（ｔ））を再取得せずにそのまま利用している。パラメータ推定部１０は、再度バッテリ４に電流信号を入力して、入出力データ（ｕ１（ｔ）及びｙ１（ｔ））を再取得してもよい。
【０１２５】
図９のフローチャートの処理において、パラメータ推定部１０は、過電圧パラメータとヒステリシスパラメータとを交互に推定する。パラメータ推定部１０は、過電圧パラメータとヒステリシスパラメータとのいずれか一方を連続して推定してもよい。
【０１２６】
図９のフローチャートの処理において、パラメータ推定部１０は、ステップＳ２７におけるヒステリシスパラメータの推定を繰り返す際に、既に取得した入出力データ（ｕ２（ｔ）及びｙ２（ｔ））を再取得せずにそのまま利用している。パラメータ推定部１０は、再度バッテリ４に電流信号を入力して、入出力データ（ｕ２（ｔ）及びｙ２（ｔ））を再度取得してもよい。パラメータ推定部１０は、振幅が比較的大きい電流信号と振幅が比較的小さい電流信号とを交互にバッテリ４に入力してもよい。パラメータ推定部１０は、振幅が比較的大きい電流信号と振幅が比較的小さい電流信号とのいずれか一方を連続してバッテリ４に入力してもよい。
【０１２７】
図９のフローチャートの処理において、パラメータ推定部１０は、第２の電流よりも先に第１の電流をバッテリ４に入力している。パラメータ推定部１０は、第１の電流よりも先に第２の電流をバッテリ４に入力してもよい。
【０１２８】
図９のフローチャートの処理において、パラメータ推定部１０は、ヒステリシスパラメータよりも先に過電圧パラメータを推定している。パラメータ推定部１０は、過電圧パラメータよりも先にヒステリシスパラメータを推定してもよい。
【０１２９】
図９のフローチャートの処理において、ＯＣＶの推定は、ステップＳ２６で行われてもよいし、ステップＳ２７で行われてもよいし、両方のステップでそれぞれ行われてもよい。Ｒ0の推定についても、ステップＳ２６で行われてもよいし、ステップＳ２７で行われてもよいし、両方のステップでそれぞれ行われてもよい。
【０１３０】
＜推定精度及び収束の速さ＞
本実施形態に係るパラメータ推定装置１は、図８及び図９に示されるフローチャートの処理を実行して、バッテリ等価回路モデル２１のパラメータを推定できる。本実施形態に係るパラメータ推定方法によれば、パラメータの推定精度を高くでき、且つ、推定値の収束を速くすることができる。
【０１３１】
図１０は、図８のパラメータ推定方法を繰り返したときの推定値の変動の一例である。図１０（Ａ）は、ヒステリシスパラメータの一つであるＣhの推定値の変動を示す。図１０（Ｂ）は、過電圧パラメータの一つであるＲdの推定値の変動を示す。図１０において、パラメータの真値は、一点鎖線で示される。本実施形態に係るパラメータ推定方法による推定値の変動は、実線で示される。比較例に係る推定値の変動は、破線で示される。
【０１３２】
＜＜ヒステリシスパラメータの推定値＞＞
図１０（Ａ）を参照して、Ｃhの推定値の変動について説明する。比較例の一つは、振幅が比較的小さい電流信号を入力した場合に取得される入出力データのみに基づいてＣhを推定する場合である。この場合、図１０（Ａ）によれば、Ｃhの推定値（破線で示される）は、パラメータ推定を繰り返すことで真値に徐々に近づく。
【０１３３】
他の比較例は、振幅が比較的大きい電流信号を入力した場合に取得される入出力データのみに基づいてパラメータを推定する場合である。この場合も、図１０（Ａ）によれば、Ｃhの推定値（破線で示される）は、パラメータ推定を繰り返すことで真値に徐々に近づく。
【０１３４】
一方で、本実施形態に係るパラメータ推定方法により算出されたＣhの推定値（実線で示される）は、図１０（Ａ）によれば、比較例に係るＣhの推定値（破線で示される）よりも真値に近く、比較例よりも真値に速く近づく。
【０１３５】
＜＜過電圧パラメータの推定値＞＞
図１０（Ｂ）を参照して、Ｒdの推定値の変動について説明する。比較例の一つは、振幅が比較的小さい電流信号を入力した場合に取得される入出力データのみに基づいてパラメータを推定する場合である。この場合、図１０（Ｂ）によれば、Ｒdの推定値（破線で示される）は、パラメータ推定を繰り返すことで真値に徐々に近づくが、真値から大きく異なる値に収束する。
【０１３６】
他の比較例は、振幅が比較的大きい電流信号を入力した場合に取得される入出力データのみに基づいてパラメータを推定する場合である。この場合、図１０（Ｂ）によれば、Ｒdの推定値（破線で示される）は、パラメータ推定を繰り返すことで真値に徐々に近づく。
【０１３７】
一方で、本実施形態に係るパラメータ推定方法により算出されたＲdの推定値（実線で示される）は、図１０（Ｂ）によれば、比較例に係るＲdの推定値よりも真値に速く近づく。
【０１３８】
以上説明してきたように、本実施形態に係るパラメータ推定装置１によれば、推定しようとするパラメータに応じて入力電流を選択して、パラメータの推定精度を向上させることができる。また、本実施形態に係るパラメータ推定装置１によれば、パラメータの推定値を効率よく真値に近づけることができる。また、本実施形態に係るパラメータ推定装置１によれば、パラメータの推定値をより速く真値に近づけることができる。
【０１３９】
（実施形態２）
実施形態１では、ヒステリシスモデルのパラメータの一つであるＲhを式（１１）の形式で表した。実施形態２では、Ｒhを表す形式として式（１１）以外の形式を用いる場合について説明する。実施形態１と重複する説明については省略する。
【０１４０】
式（１１）ではヒステリシスモデルの可変抵抗（Ｒh）が入力電流（ｕ（ｔ））の絶対値｜ｕ（ｔ）｜を分母とする関数で表される。一方、ヒステリシス現象をさらに正確に取り扱うために、可変抵抗（Ｒh）を表す式の分母を入力電流（ｕ（ｔ））の関数として拡張して表すことが考えられる。ここで、入力電流（ｕ（ｔ））の関数として、ｆ（｜ｕ（ｔ）｜）を定義する。関数ｆ（ｘ）は入力ｘに対する出力を表し、入力ｘと出力ｆ（ｘ）の関係は任意に定められる。関数ｆ（｜ｕ（ｔ）｜）を用いて可変抵抗（Ｒh）を表すと、以下の式（３４）のようになる。
【数３３】
0006711981.tif 000034
式（３４）においても、モデルの抵抗は電流の大きさによって可変となる可変抵抗であることが特徴である。
【０１４１】
関数ｆ（ｘ）の形式の一つに一次関数がある。これは、ｆ（ｘ）＝αｘ＋βと表す形式であり、α及びβは定数である。この形式で可変抵抗（Ｒh）を表すと、式（３５）のようになる。
【数３４】
0006711981.tif 000035
式（３５）においてβ＞０であれば、ｕ（ｔ）＝０の場合に可変抵抗（Ｒh）は有限の値をとる。すなわち可変抵抗（Ｒh）を推定するときに無限大に発散することがなくモデルが安定し、実際のモデルに近づく。また式（３５）においてβ＝１であれば、ｕ（ｔ）＝０の場合にＲh（ｔ）＝Ｍ（ｔ）となるため、モデルのパラメータの理解が容易になる。このように本実施形態に係るパラメータ推定装置１によれば、ヒステリシスを正確に扱えるため、より早く正確な推定値を得ることができる。特に、ｆ（｜ｕ（ｔ）｜）＝α｜ｕ（ｔ）｜＋βと表す形式とした式（３５）のモデルによれば、より早く正確な推定値を得ることができる。
【０１４２】
さらに式（３５）においてα＝１、β＝０とすれば、ｆ（｜ｕ（ｔ）｜）＝｜ｕ（ｔ）｜であり、式（１１）の形式でモデルを表すこととなる。これによりモデル構成がより簡易にできるため、より容易にヒステリシス現象を取り扱うことができる。
【０１４３】
以上、関数ｆ（ｘ）の形式として一次関数について説明したが、これに限られない。二次関数などの多項式であってもよいし、有理関数、無理関数、対数関数や指数関数などであってもよい。いずれの形式であっても、バッテリ４の等価回路にヒステリシス素子を増やすことなくヒステリシスを扱うことができ、ヒステリシス特性に則したモデルとして係数を選ぶことにより正確な推定値を得ることができる。
【０１４４】
本発明を諸図面および実施例に基づき説明してきたが、当業者であれば本開示に基づき種々の変形または修正をおこなうことが容易であることに注意されたい。従って、これらの変形または修正は本発明の範囲に含まれることに留意されたい。例えば、各構成部、各ステップなどに含まれる機能などは論理的に矛盾しないように再配置可能であり、複数の構成部およびステップなどを１つに組み合わせたり、或いは分割したりすることが可能である。
【０１４５】
例えば、上述の実施形態において、ワールブルグインピーダンス（Ｚw）を無限級数展開又は連分数展開により近似したが、任意の方法で近似してもよい。例えば、無限乗積展開を用いて近似することが考えられる。
【符号の説明】
【０１４６】
１  パラメータ推定装置
２  電流センサ（入力検出部）
３  電圧センサ（出力検出部）
４  バッテリ
５  電源装置
１０  パラメータ推定部
１１  出力予測部
１２  誤差算出部
１３  パラメータ更新部
２０  記憶部
２１  バッテリ等価回路モデル
３０  充電率算出部
４０  電荷量算出部
５０  健全度算出部
【図１】
0006711981.tif 000040
【図２】
0006711981.tif 000041
【図３】
0006711981.tif 000042
【図４】
0006711981.tif 000043
【図５】
0006711981.tif 000044
【図６】
0006711981.tif 000045
【図７】
0006711981.tif 000046
【図８】
0006711981.tif 000047
【図９】
0006711981.tif 000048
【図１０】
0006711981.tif 000049