知財求人 - 知財ポータルサイト「IP Force」
▶ カルソニックカンセイ株式会社の特許一覧 ▶ 学校法人慶應義塾の特許一覧 ▶ 国立大学法人京都大学の特許一覧 ▶ 国立大学法人群馬大学の特許一覧
(19)【発行国】日本国特許庁(JP)
(12)【公報種別】特許公報(B1)
(11)【特許番号】
(24)【登録日】2023-03-31
(45)【発行日】2023-04-10
(54)【発明の名称】推定装置および推定方法
(51)【国際特許分類】
G01R 31/367 20190101AFI20230403BHJP
G01R 31/389 20190101ALI20230403BHJP
G01R 31/388 20190101ALI20230403BHJP
G01R 31/3828 20190101ALI20230403BHJP
H01M 10/48 20060101ALI20230403BHJP
H02J 7/00 20060101ALI20230403BHJP
【ＦＩ】
G01R31/367
G01R31/389
G01R31/388
G01R31/3828
H01M10/48 P
H02J7/00 Y
【請求項の数】
7
(21)【出願番号】P 2022172918
(22)【出願日】2022-10-28
【審査請求日】2022-12-16
【早期審査対象出願】
(73)【特許権者】
【識別番号】000004765
【氏名又は名称】マレリ株式会社
(73)【特許権者】
【識別番号】899000079
【氏名又は名称】慶應義塾
(73)【特許権者】
【識別番号】504132272
【氏名又は名称】国立大学法人京都大学
(73)【特許権者】
【識別番号】504145364
【氏名又は名称】国立大学法人群馬大学
(74)【代理人】
【識別番号】100148301
【弁理士】
【氏名又は名称】竹原 尚彦
(74)【代理人】
【識別番号】100176991
【弁理士】
【氏名又は名称】中島 由布子
(74)【代理人】
【識別番号】100217696
【弁理士】
【氏名又は名称】川口 英行
(72)【発明者】
【氏名】長村 謙介
(72)【発明者】
【氏名】川口 貴弘
(72)【発明者】
【氏名】足立 修一
(72)【発明者】
【氏名】丸田 一郎
【審査官】小川 浩史
(56)【参考文献】
【文献】特開２０２０－１１２５５１（ＪＰ，Ａ）
【文献】川口貴弘；大山隆景；丸田一郎；長村謙介；片芝惇平；足立修一，「μ-マルコフモデルを用いた二次電池の実効抵抗の逐次推定」，自動車技術会論文集，2021年11月25日，vol. 52, No. 6，pp. 1323-1328，DOI： https://doi.org/10.11351/jsaeronbun.52.1323
【文献】大山隆景；足立修一；丸田一郎；片芝惇平；長村謙介，「ベクトル型忘却要素を用いた時変インピーダンスの逐次推定」，第６３回自動制御連合講演会講演論文集，2020年12月14日，pp. 351-354，DOI： https://doi.org/10.11511/jacc.63.0_351
【文献】佐々木理沙子；川口貴弘；丸田一郎；足立修一；片芝惇平；長村謙介，「μ-マルコフモデルを用いた自動車用二次電池の内部インピーダンスの推定」，計測自動制御学会論文集，2020年02月15日，Vol. 56, No. 2，pp. 67-69，DOI： https://doi.org/10.9746/sicetr.56.67
(58)【調査した分野】(Int.Cl.，ＤＢ名)
Ｇ０１Ｒ ３１／３６－３１／３９６
Ｈ０１Ｍ １０／４２－１０／４８
Ｈ０２Ｊ ７／００
(57)【特許請求の範囲】
【請求項1】
バッテリの所定時間における過電圧変化を示す実効抵抗を推定する推定装置であって前記バッテリを流れる電流と、前記バッテリの端子電圧の測定値に基づいて、前記バッテリの過電圧を算出する過電圧算出部と、
前記電流と前記過電圧の時系列データに基づいて、ＦＩＲモデルを含むモデルを構成する多項式の各項の係数を推定することで、前記バッテリのシステムを同定するシステム同定部と、
前記バッテリのシステムに基づいて、前記実効抵抗を算出する実効抵抗算出部と、を備え、
前記システム同定部は、前記ＦＩＲモデルを構成する多項式を複数の時間区分で分割し、各時間区分に含まれる項に対して共通の係数を推定することで、前記ＦＩＲモデルの多項式の項数よりも少ない数の係数を推定する、推定装置。
【請求項2】
前記電流と前記過電圧の時系列データに対して、白色化フィルタによりフィルタリング処理を行い、前記電流が白色信号である前記時系列データを前記システム同定部に出力するフィルタリング部を備える、請求項１記載の推定装置。
【請求項3】
前記実効抵抗算出部は、前記時系列データの入力開始時刻から前記所定時間経過後の最後の時刻までの、各時間区分で推定された前記係数に、各時間区分の項数を掛けた数値を積算することで前記実効抵抗を算出し、前記システム同定部において、前記時間区分のいずれかの境界と前記最後の時刻が一致するように、前記時間区分の長さが設定される、請求項１記載の推定装置。
【請求項4】
前記システム同定部において、前記時系列データにおける入力開始時刻に近い時間区分ほど長さが短くなるように、各時間区分の長さが設定される、請求項１記載の推定装置。
【請求項5】
前記システム同定部は、前記各時間区分に含まれる項に対して、前記共通の係数として定数を推定する、請求項１記載の推定装置。
【請求項6】
前記システム同定部は、前記ＦＩＲモデルを含むモデルとして、前記ＦＩＲモデルにＡＲＸモデルを組み合わせたμマルコフモデルを適用し、前記時系列データの入力開始時刻から前記所定時間経過後の最後の時刻までの前記時系列データを前記ＦＩＲモデルに適用する、請求項１記載の推定装置。
【請求項7】
バッテリの所定時間における過電圧変化を示す実効抵抗を推定する推定方法であって、前記バッテリを流れる電流と、前記バッテリの端子電圧の測定値に基づいて、前記バッテリの過電圧を算出する過電圧算出ステップと、
前記電流と前記過電圧の時系列データに基づいて、ＦＩＲモデルを含むモデルを構成する多項式の各項の係数を推定することで、前記バッテリのシステムを同定するシステム同定ステップと、
前記バッテリのシステムに基づいて、前記実効抵抗を算出する実効抵抗算出ステップと、を備え、
前記システム同定ステップでは、前記ＦＩＲモデルを構成する多項式を複数の時間区分で分割し、各時間区分に含まれる項に対して共通の係数を推定することで、前記ＦＩＲモデルの多項式の項数よりも少ない数の係数を推定する、推定方法。
【発明の詳細な説明】
【技術分野】
【0001】
本発明は、推定装置および推定方法に関する。
【背景技術】
【0002】
バッテリは、例えば、モータ等の車両の駆動源に電力を供給するため用いられる。
バッテリの充電制御や車両の運転制御に際して、バッテリのＳＯＰ(State of Power)を推定することが求められる。ＳＯＰは、バッテリが所定時間に充電または放電可能電力を意味するものである。ＳＯＰは、例えば、バッテリの実効抵抗から推定することができる。
バッテリの実効抵抗は、バッテリの内部インピーダンスの特性を表現する特徴量である。実効抵抗を推定する方法として、ＦＩＲ(Finite Impulse Response)モデルや、ＦＩＲモデルとＡＲＸ（Auto-Regressive eXogeneous）モデルを組み合わせたμマルコフモデルにより、バッテリのシステムを同定する方法が提案されている（例えば、特許文献１参照）。
バッテリの充電制御や車両の運転制御に際して、バッテリのＳＯＰ(State of Power)を推定することが求められる。ＳＯＰは、バッテリが所定時間に充電または放電可能電力を意味するものである。ＳＯＰは、例えば、バッテリの実効抵抗から推定することができる。
バッテリの実効抵抗は、バッテリの内部インピーダンスの特性を表現する特徴量である。実効抵抗を推定する方法として、ＦＩＲ(Finite Impulse Response)モデルや、ＦＩＲモデルとＡＲＸ（Auto-Regressive eXogeneous）モデルを組み合わせたμマルコフモデルにより、バッテリのシステムを同定する方法が提案されている（例えば、特許文献１参照）。
【先行技術文献】
【特許文献】
【0003】
【文献】特許第６７３７４９０号公報
【発明の概要】
【発明が解決しようとする課題】
【0004】
実効抵抗は、ＳＯＰの推定の他にも様々な処理に用いることができ、より長い時間の実効抵抗を推定するニーズがある。実効抵抗を求める時間が長くなると、ＦＩＲモデルにおいて推定するパラメータの数が増加し、演算負荷の増大に繋がる。実効抵抗を推定する際に、推定するパラメータの数を削減することが求められている。
【課題を解決するための手段】
【0005】
本発明の一態様は、バッテリの所定時間における過電圧変化を示す実効抵抗を推定する推定装置であって、
前記バッテリを流れる電流と、前記バッテリの端子電圧の測定値に基づいて、前記バッテリの過電圧を算出する過電圧算出部と、
前記電流と前記過電圧の時系列データに基づいて、ＦＩＲモデルを含むモデルを構成する多項式の各項の係数を推定することで、前記バッテリのシステムを同定するシステム同定部と、
前記バッテリのシステムに基づいて、前記実効抵抗を算出する実効抵抗算出部と、を備え、
前記システム同定部は、前記ＦＩＲモデルを構成する多項式を複数の時間区分で分割し、各時間区分に含まれる項に対して共通の係数を推定することで、前記ＦＩＲモデルの多項式の項数よりも少ない数の係数を推定する。
前記バッテリを流れる電流と、前記バッテリの端子電圧の測定値に基づいて、前記バッテリの過電圧を算出する過電圧算出部と、
前記電流と前記過電圧の時系列データに基づいて、ＦＩＲモデルを含むモデルを構成する多項式の各項の係数を推定することで、前記バッテリのシステムを同定するシステム同定部と、
前記バッテリのシステムに基づいて、前記実効抵抗を算出する実効抵抗算出部と、を備え、
前記システム同定部は、前記ＦＩＲモデルを構成する多項式を複数の時間区分で分割し、各時間区分に含まれる項に対して共通の係数を推定することで、前記ＦＩＲモデルの多項式の項数よりも少ない数の係数を推定する。
【0006】
また、本発明の一態様は、バッテリの所定時間における過電圧変化を示す実効抵抗を推定する推定方法であって
前記バッテリを流れる電流と、前記バッテリの端子電圧の測定値に基づいて、前記バッテリの過電圧を算出する過電圧算出ステップと、
前記電流と前記過電圧の時系列データに基づいて、ＦＩＲモデルを含むモデルを構成する多項式の各項の係数を推定することで、前記バッテリのシステムを同定するシステム同定ステップと、
前記バッテリのシステムに基づいて、前記実効抵抗を算出する実効抵抗算出ステップと、を備え、
前記システム同定ステップでは、前記ＦＩＲモデルを構成する多項式を複数の時間区分で分割し、各時間区分に含まれる項に対して共通の係数を推定することで、前記ＦＩＲモデルの多項式の項数よりも少ない数の係数を推定する。
前記バッテリを流れる電流と、前記バッテリの端子電圧の測定値に基づいて、前記バッテリの過電圧を算出する過電圧算出ステップと、
前記電流と前記過電圧の時系列データに基づいて、ＦＩＲモデルを含むモデルを構成する多項式の各項の係数を推定することで、前記バッテリのシステムを同定するシステム同定ステップと、
前記バッテリのシステムに基づいて、前記実効抵抗を算出する実効抵抗算出ステップと、を備え、
前記システム同定ステップでは、前記ＦＩＲモデルを構成する多項式を複数の時間区分で分割し、各時間区分に含まれる項に対して共通の係数を推定することで、前記ＦＩＲモデルの多項式の項数よりも少ない数の係数を推定する。
【発明の効果】
【0007】
本発明によれば、推定するパラメータの数を削減することができる。
【図面の簡単な説明】
【0008】
【図1】本発明の実施形態に係る推定装置の構成を示すブロック図である。
【図2】実効抵抗を説明する図である。
【図3】過電圧算出部の構成を示すブロック図である。
【図4】推定部の構成を示すブロック図である。
【図5】バッテリシステムを表現する等価回路を示す図である。
【図6】本実施形態のＰＣＦＩＲモデルと、従来のＦＩＲモデルの概念図である。
【図7】ＰＣＦＩＲモデルにおける時間区分を説明する図である。
【図8】本実施形態における推定装置の処理を示すフローチャートである。
【図9】数値シミュレーションに用いた回路構成を示す図である。
【図10】（ａ）は数値シミュレーションにおける入力データを示す図であり、（ｂ）は出力データを示す図である。
【図11】本実施形態のＰＣＦＩＲモデルと従来のＦＩＲモデルによる実効抵抗の推定結果を示すグラフである。（ａ）は、各モデルにおけるインパルス応答である電圧を示す図であり、各モデルで推定されるパラメータに対応するものである。（ｂ）は、各モデルにおけるステップ応答である実効抵抗を示す図である。
【図12】各時間区分の長さを一定とした場合の推定結果であり、（ａ）は、インパルス応答を示し、（ｂ）はステップ応答である実効抵抗を示している。
【図13】各時間区分の長さを段階的に長くした場合の推定結果である。（ａ）は、インパルス応答を示し、（ｂ）はステップ応答である実効抵抗を示している。
【図14】変形例１に係る区分線形ＦＩＲモデルと、従来のＦＩＲモデルの概念図である。
【図15】変形例２に係る推定部の構成を示すブロック図である。
【発明を実施するための形態】
【0009】
以下、本発明の実施形態に係る推定装置を、図面を参照して説明する。
図１は、実施形態に係る推定装置の構成を示すブロック図である。
図１に示すように、推定装置１は、電気自動車またはハイブリッド電気自動車等の車両に設けられたバッテリ５０の実効抵抗を推定するものである。バッテリ５０には、電圧センサ６０および電流センサ７０が接続されている。バッテリ５０は、充電可能な二次電池であり、例えばリチウムイオンバッテリを用いることができるが、他の種類のバッテリであっても良い。バッテリ５０は、車両に設けられている。バッテリ５０は放電することにより、車両を駆動する電気モータへ電力を供給する。また、車両の制動時には、電気モータからの回生エネルギーにより充電される。バッテリ５０はまた、急速充電器や家庭用コンセント等の外部の充電設備によっても充電される。
図１は、実施形態に係る推定装置の構成を示すブロック図である。
図１に示すように、推定装置１は、電気自動車またはハイブリッド電気自動車等の車両に設けられたバッテリ５０の実効抵抗を推定するものである。バッテリ５０には、電圧センサ６０および電流センサ７０が接続されている。バッテリ５０は、充電可能な二次電池であり、例えばリチウムイオンバッテリを用いることができるが、他の種類のバッテリであっても良い。バッテリ５０は、車両に設けられている。バッテリ５０は放電することにより、車両を駆動する電気モータへ電力を供給する。また、車両の制動時には、電気モータからの回生エネルギーにより充電される。バッテリ５０はまた、急速充電器や家庭用コンセント等の外部の充電設備によっても充電される。
【0010】
電圧センサ６０はバッテリ５０の端子電圧（terminal voltage）を測定する。電流センサ７０はバッテリ５０を流れる電流を測定する。電圧センサ６０および電流センサ７０は、例えば、車両のイグニッションＯＮからイグニッションＯＦＦまでの間、所定のサンプリング周期で端子電圧および電流の測定を行う。推定装置１は、電圧センサ６０および電流センサ７０に有線または無線により接続され、それぞれが測定した端子電圧および電流の測定値が、順次入力される。
【0011】
推定装置１は、バッテリ５０の端子電圧と電流の測定値を用いて、バッテリ５０の実効抵抗(effective resistance)を推定する。
実効抵抗は、二次電池であるバッテリ５０の内部インピーダンスの特性を表現する特徴量であり、直流抵抗や通電抵抗とも呼ばれる。
図２は、実効抵抗を説明する図である。
図２に示すように、実効抵抗は、具体的には、大きさ１Ａのステップ信号の入力に対する過電圧変化を意味するものである。以降の説明では、所定時間（ｘ秒）経過後の実効抵抗を、ｘ秒実効抵抗という。推定装置１は、ユーザによって指定されるｘ秒実効抵抗を推定する。
実効抵抗は、例えば、バッテリ５０のＳＯＰ（State of Power）の推定に用いることができる。ＳＯＰは充電可能電力と放電可能電力の総称であるが、以下の説明において、ＳＯＰは、端子電圧の下限が与えられたときの放電可能電力を意味するものとして用いる。すなわち、ｘ秒実効抵抗から、バッテリ５０がｘ秒間継続して放電可能なパワーを推定することができる。
実効抵抗は、二次電池であるバッテリ５０の内部インピーダンスの特性を表現する特徴量であり、直流抵抗や通電抵抗とも呼ばれる。
図２は、実効抵抗を説明する図である。
図２に示すように、実効抵抗は、具体的には、大きさ１Ａのステップ信号の入力に対する過電圧変化を意味するものである。以降の説明では、所定時間（ｘ秒）経過後の実効抵抗を、ｘ秒実効抵抗という。推定装置１は、ユーザによって指定されるｘ秒実効抵抗を推定する。
実効抵抗は、例えば、バッテリ５０のＳＯＰ（State of Power）の推定に用いることができる。ＳＯＰは充電可能電力と放電可能電力の総称であるが、以下の説明において、ＳＯＰは、端子電圧の下限が与えられたときの放電可能電力を意味するものとして用いる。すなわち、ｘ秒実効抵抗から、バッテリ５０がｘ秒間継続して放電可能なパワーを推定することができる。
【0012】
バッテリ５０のＳＯＰは、バッテリ５０の充放電制御や車両の運転制御に利用することができる。
例えば、ハイブリッド車両において、エンジン走行中に必要に応じてモータを駆動してエンジンをアシストすることがある。この場合、バッテリ５０は、モータがエンジンをアシストする時間（例えば１０秒）分、モータを駆動するための電流を継続して放電する必要がある。
また、例えば、バッテリ５０は、エンジンがアイドリングストップから再始動する際に、クランキングを行うための時間（例えば、０．１秒）分、モータを駆動するための電流を継続して放電する必要がある。
すなわち、１０秒実効抵抗や０．１秒実効抵抗等のｘ秒実効抵抗から推定されたＳＯＰを用いることで、車両がモータアシスト走行やクランキングを行うことができるかを判定することができる。
例えば、ハイブリッド車両において、エンジン走行中に必要に応じてモータを駆動してエンジンをアシストすることがある。この場合、バッテリ５０は、モータがエンジンをアシストする時間（例えば１０秒）分、モータを駆動するための電流を継続して放電する必要がある。
また、例えば、バッテリ５０は、エンジンがアイドリングストップから再始動する際に、クランキングを行うための時間（例えば、０．１秒）分、モータを駆動するための電流を継続して放電する必要がある。
すなわち、１０秒実効抵抗や０．１秒実効抵抗等のｘ秒実効抵抗から推定されたＳＯＰを用いることで、車両がモータアシスト走行やクランキングを行うことができるかを判定することができる。
【0013】
推定装置１は、例えば、車両に設けられたＥＣＵ（Electronic Control Unit）から構成することができる。
ＥＣＵは、図示は省略するが、ＣＰＵ等のプロセッサと、ＲＯＭおよびＲＡＭ等のメモリから構成される。メモリには、推定装置１で実行される各種プログラムが格納されており、プロセッサがプログラムを実行することで、図１に示す機能構成が実現される。メモリには、また、推定装置１が行う処理に必要なデータが格納され、さらに推定装置１の処理結果が一時的に記憶される。
詳細な説明は省略するが、ＥＣＵは、実効抵抗の推定処理に加えて、例えば、実効抵抗を用いたＳＯＰの推定処理や、ＳＯＰに基づいたバッテリ５０の充放電制御および車両の運転制御を行っても良い。あるいは、ＥＣＵは推定した実効抵抗を外部に出力し、外部の制御装置において実効抵抗を用いた処理を行っても良い。
あるいは、推定装置１は、車両の外部に設けられたコンピューターとしても良い。この場合、車両外部の推定装置１は、例えば、車両に設けられたＥＣＵと通信を行って電圧センサ６０および電流センサ７０の測定値を取得しても良い。
ＥＣＵは、図示は省略するが、ＣＰＵ等のプロセッサと、ＲＯＭおよびＲＡＭ等のメモリから構成される。メモリには、推定装置１で実行される各種プログラムが格納されており、プロセッサがプログラムを実行することで、図１に示す機能構成が実現される。メモリには、また、推定装置１が行う処理に必要なデータが格納され、さらに推定装置１の処理結果が一時的に記憶される。
詳細な説明は省略するが、ＥＣＵは、実効抵抗の推定処理に加えて、例えば、実効抵抗を用いたＳＯＰの推定処理や、ＳＯＰに基づいたバッテリ５０の充放電制御および車両の運転制御を行っても良い。あるいは、ＥＣＵは推定した実効抵抗を外部に出力し、外部の制御装置において実効抵抗を用いた処理を行っても良い。
あるいは、推定装置１は、車両の外部に設けられたコンピューターとしても良い。この場合、車両外部の推定装置１は、例えば、車両に設けられたＥＣＵと通信を行って電圧センサ６０および電流センサ７０の測定値を取得しても良い。
【0014】
図１に示すように、推定装置１は、過電圧算出部２０および推定部３０を備える。
過電圧算出部２０は、電圧センサ６０および電流センサ７０によって所定周期でサンプリングされた、バッテリ５０の端子電圧ｖ（ｋ）および電流ｉ（ｋ）を取得し、過電圧η（ｋ）を算出する。
推定部３０は、バッテリ５０の電流ｉ（ｋ）を、後記するＰＣＦＩＲ（Piecewise Constant FIR）モデルにおける入力ｕ（ｋ）とし、過電圧η（ｋ）を出力ｙ（ｋ）として、バッテリ５０のシステムを同定することで、実効抵抗Ｒ（Ｋ）を推定する。
過電圧算出部２０は、電圧センサ６０および電流センサ７０によって所定周期でサンプリングされた、バッテリ５０の端子電圧ｖ（ｋ）および電流ｉ（ｋ）を取得し、過電圧η（ｋ）を算出する。
推定部３０は、バッテリ５０の電流ｉ（ｋ）を、後記するＰＣＦＩＲ（Piecewise Constant FIR）モデルにおける入力ｕ（ｋ）とし、過電圧η（ｋ）を出力ｙ（ｋ）として、バッテリ５０のシステムを同定することで、実効抵抗Ｒ（Ｋ）を推定する。
【0015】
図３は、過電圧算出部２０の構成を示すブロック図である。
図３に示すように、過電圧算出部２０は、ＳＯＣ算出部２１、ＳＯＣ－ＯＣＶ変換部２２および演算器２３を備える。
ＳＯＣ算出部２１は、端子電圧ｖ（ｋ）および電流ｉ（ｋ）を用いて、バッテリ５０の充電率であるＳＯＣ（State of Charge）ｘc（ｋ）を算出する。ＳＯＣ算出部２１は、電流積算法またはカルマンフィルタ等の既知の手法により、ＳＯＣｘc（ｋ）を算出することができる。
ＳＯＣ－ＯＣＶ変換部２２は、既知のＳＯＣ－ＯＣＶ特性ｆocv（・）に基づいて、ＳＯＣ算出部２１が算出したＳＯＣｘc（ｋ）を、開回路電圧（Open Circuit Voltage: ＯＣＶ）ｖoc（ｋ）に変換する。開回路電圧ｖoc（ｋ）は、バッテリ５０に電流が流れていない際の電圧を意味する。
演算器２３は、バッテリ５０の端子電圧ｖ（ｋ）と、開回路電圧ｖoc（ｋ）の差を演算し、バッテリ５０の過電圧η（ｋ）として推定部３０に出力する。
図３に示すように、過電圧算出部２０は、ＳＯＣ算出部２１、ＳＯＣ－ＯＣＶ変換部２２および演算器２３を備える。
ＳＯＣ算出部２１は、端子電圧ｖ（ｋ）および電流ｉ（ｋ）を用いて、バッテリ５０の充電率であるＳＯＣ（State of Charge）ｘc（ｋ）を算出する。ＳＯＣ算出部２１は、電流積算法またはカルマンフィルタ等の既知の手法により、ＳＯＣｘc（ｋ）を算出することができる。
ＳＯＣ－ＯＣＶ変換部２２は、既知のＳＯＣ－ＯＣＶ特性ｆocv（・）に基づいて、ＳＯＣ算出部２１が算出したＳＯＣｘc（ｋ）を、開回路電圧（Open Circuit Voltage: ＯＣＶ）ｖoc（ｋ）に変換する。開回路電圧ｖoc（ｋ）は、バッテリ５０に電流が流れていない際の電圧を意味する。
演算器２３は、バッテリ５０の端子電圧ｖ（ｋ）と、開回路電圧ｖoc（ｋ）の差を演算し、バッテリ５０の過電圧η（ｋ）として推定部３０に出力する。
【0016】
図４は、推定部３０の構成を示すブロック図である。
図４に示すように、推定部３０は、フィルタリング部３１、特徴ベクトル算出部３２、パラメータ推定部３３および実効抵抗算出部３４を備える。
前記したように、推定部３０は、バッテリ５０の電流ｉ（ｋ）を、ＰＣＦＩＲモデルにおける入力ｕ（ｋ）として、過電圧η（ｋ）を出力ｙ（ｋ）として扱う。
フィルタリング部３１は、入力である電流ｕ（ｋ）と、出力である過電圧ｙ（ｋ）に対して、白色化フィルタを用いたフィルタリング処理を行い、フィルタリングされた電流ｕf（ｋ）と過電圧ｙf（ｋ）を特徴ベクトル算出部３２に出力する。白色化フィルタは、例えば、予め設定され推定装置１のメモリに格納されたものを用いることができる。
図４に示すように、推定部３０は、フィルタリング部３１、特徴ベクトル算出部３２、パラメータ推定部３３および実効抵抗算出部３４を備える。
前記したように、推定部３０は、バッテリ５０の電流ｉ（ｋ）を、ＰＣＦＩＲモデルにおける入力ｕ（ｋ）として、過電圧η（ｋ）を出力ｙ（ｋ）として扱う。
フィルタリング部３１は、入力である電流ｕ（ｋ）と、出力である過電圧ｙ（ｋ）に対して、白色化フィルタを用いたフィルタリング処理を行い、フィルタリングされた電流ｕf（ｋ）と過電圧ｙf（ｋ）を特徴ベクトル算出部３２に出力する。白色化フィルタは、例えば、予め設定され推定装置１のメモリに格納されたものを用いることができる。
【0017】
特徴ベクトル算出部３２は、フィルタリングされた電流ｕf（ｋ）を入力ｕ（ｋ）とし、過電圧ｙf（ｋ）を出力ｙ（ｋ）として、ＰＣＦＩＲモデルを構成する回帰ベクトルφを算出する。
パラメータ推定部３３は、回帰ベクトルφを用いて、ＰＣＦＩＲモデルを構成するパラメータベクトルθを推定する。
すなわち、特徴ベクトル算出部３２は、明らかになっている入力ｕ（ｋ）および出力ｙ（ｋ）から、ＰＣＦＩＲモデルの特徴ベクトルである回帰ベクトルφを算出する。パラメータ推定部３３は、算出された回帰ベクトルφから、ＰＣＦＩＲモデルのパラメータ（係数）を推定する。これによって、バッテリ５０のシステムが同定される。すなわち、特徴ベクトル算出部３２およびパラメータ推定部３３は、システム同定部として機能する。
実効抵抗算出部３４は、同定されたバッテリシステムを用いて、バッテリ５０の実効抵抗Ｒ（Ｋ）を算出する。
推定部３０は、実効抵抗算出部３４が算出した実効抵抗を外部に出力しても良く、あるいは、ＥＣＵのメモリに格納して、バッテリの充放電制御等の他の制御に用いても良い。
パラメータ推定部３３は、回帰ベクトルφを用いて、ＰＣＦＩＲモデルを構成するパラメータベクトルθを推定する。
すなわち、特徴ベクトル算出部３２は、明らかになっている入力ｕ（ｋ）および出力ｙ（ｋ）から、ＰＣＦＩＲモデルの特徴ベクトルである回帰ベクトルφを算出する。パラメータ推定部３３は、算出された回帰ベクトルφから、ＰＣＦＩＲモデルのパラメータ（係数）を推定する。これによって、バッテリ５０のシステムが同定される。すなわち、特徴ベクトル算出部３２およびパラメータ推定部３３は、システム同定部として機能する。
実効抵抗算出部３４は、同定されたバッテリシステムを用いて、バッテリ５０の実効抵抗Ｒ（Ｋ）を算出する。
推定部３０は、実効抵抗算出部３４が算出した実効抵抗を外部に出力しても良く、あるいは、ＥＣＵのメモリに格納して、バッテリの充放電制御等の他の制御に用いても良い。
【0018】
［バッテリシステム］
図５は、バッテリシステムを表現する等価回路を示す図である。
図５において，ｉ（ｋ）はバッテリ５０に流れる電流であり、放電方向を正とする。また、ｖ（ｋ）は、端子電圧であり、ｖoc（ｋ）は、開回路電圧である。ここで、ｋは離散時刻を表し、たとえばｉ（ｋ）は、時刻ｋＴsにおける電流ｉの値を意味する。ただし、Ｔｓはサンプリング周期である。
図５に示すように、バッテリシステムは、起電力と内部インピーダンスからなる等価回路として表現することができる。起電力における電圧が開回路電圧ｖoc（ｋ）であり、内部インピーダンスにおける電圧が過電圧η（ｋ）である。図５のグラフに示すように、開回路電圧は、バッテリ５０の充電率を示すＳＯＣ（State of Charge）に対して、いわゆるＳＯＣ－ＯＣＶ特性と呼ばれる対応関係を有する。
図５は、バッテリシステムを表現する等価回路を示す図である。
図５において，ｉ（ｋ）はバッテリ５０に流れる電流であり、放電方向を正とする。また、ｖ（ｋ）は、端子電圧であり、ｖoc（ｋ）は、開回路電圧である。ここで、ｋは離散時刻を表し、たとえばｉ（ｋ）は、時刻ｋＴsにおける電流ｉの値を意味する。ただし、Ｔｓはサンプリング周期である。
図５に示すように、バッテリシステムは、起電力と内部インピーダンスからなる等価回路として表現することができる。起電力における電圧が開回路電圧ｖoc（ｋ）であり、内部インピーダンスにおける電圧が過電圧η（ｋ）である。図５のグラフに示すように、開回路電圧は、バッテリ５０の充電率を示すＳＯＣ（State of Charge）に対して、いわゆるＳＯＣ－ＯＣＶ特性と呼ばれる対応関係を有する。
【0019】
一般に、バッテリ５０の放電時には端子電圧は開回路電圧よりも低くなり、充電時には端子電圧は開回路電圧よりも高くなる。これは電池の内部インピーダンスの存在によるものであり、このような端子電圧と開回路電圧の差が、過電圧（overpotential）と呼ばれる。
すなわち、過電圧η(k)は、以下の式（１）で表される。
ただし、ｖ（ｋ）、ｖoc（ｋ）、η（ｋ）の方向は、図５に示したとおりである。
式（１）において、電流ｉ（ｋ）と開回路電圧ｖoc（ｋ）は、電池のＳＯＣを介して対応づけられることが知られている。ＳＯＣは、例えば、電流積算法により求めることができる。電池のＳＯＣをｘc（ｋ）とすると、ＳＯＣｘc（ｋ）は、電流ｉ（ｋ）を用いて以下の式（２）により求められる。
ただし、ＦＣＣは満充電容量である。また、電流ｉ（ｋ）はサンプル点間で０次ホールドされているものとした。ＳＯＣは、他にも、カルマンフィルタ等を用いて推定可能である。
開回路電圧ｖoc（ｋ）は、ＳＯＣｘc（ｋ）に依存して、以下の式（３）により表される。
ここで、ｆocv（・）は、図５のグラフに示す、ＳＯＣ－ＯＣＶ特性と呼ばれる非線形関数である。
すなわち、過電圧η(k)は、以下の式（１）で表される。
【数1】
式（１）において、電流ｉ（ｋ）と開回路電圧ｖoc（ｋ）は、電池のＳＯＣを介して対応づけられることが知られている。ＳＯＣは、例えば、電流積算法により求めることができる。電池のＳＯＣをｘc（ｋ）とすると、ＳＯＣｘc（ｋ）は、電流ｉ（ｋ）を用いて以下の式（２）により求められる。
【数2】
開回路電圧ｖoc（ｋ）は、ＳＯＣｘc（ｋ）に依存して、以下の式（３）により表される。
【数3】
【0020】
電流ｉ（ｋ）と過電圧η（ｋ）の関係、すなわち内部インピーダンスの特性を記述することを考える。内部インピーダンスのモデルは、大きく静的モデルと動的モデルに分けられる。静的モデルでは、電流ｉ（ｋ）と過電圧η（ｋ）の関係は以下の式（４）により表される。
ただし、Ｒは内部抵抗の値である。
一方、動的モデルでは電流ｉ（ｋ）と過電圧η（ｋ）の関係は、伝達関数Ｇ（ｑ）を用いて以下の式（５）により記述される。
ただし、ｑは、ｑｘ（ｋ）＝ｘ（ｋ＋１）を満たすシフトオペレータである。
【数4】
一方、動的モデルでは電流ｉ（ｋ）と過電圧η（ｋ）の関係は、伝達関数Ｇ（ｑ）を用いて以下の式（５）により記述される。
【数5】
【0021】
［放電可能電力と実効抵抗］
前記したように、バッテリ５０の充放電制御や車両の運転制御において、ＳＯＰ（放電可能電力）を推定することが求められている。
図５に示したバッテリシステムの内部インピーダンスが、式（４）の静的モデルで表される場合、端子電圧を下限値ｖminに保つと、バッテリ５０を流れる電流は、ｉ（ｋ）＝（ｖoc（ｋ）－ｖmin）／Ｒとなる。したがって、バッテリ５０から外部回路に供給される電力は、以下の式（６）で表される。
放電可能電力は、区間ｋ∈［０，Ｋ］の間で継続して取り出せる電力、すなわち，ｋ∈［０，Ｋ］におけるｐ（ｋ）の最小値である。区間ｋ∈［０，Ｋ］は、推定が要求される時間（ｘ秒）に応じてあらかじめ決定される。放電可能電力ｐdisは、以下の式（７）で表される。
前記したように、バッテリ５０の充放電制御や車両の運転制御において、ＳＯＰ（放電可能電力）を推定することが求められている。
図５に示したバッテリシステムの内部インピーダンスが、式（４）の静的モデルで表される場合、端子電圧を下限値ｖminに保つと、バッテリ５０を流れる電流は、ｉ（ｋ）＝（ｖoc（ｋ）－ｖmin）／Ｒとなる。したがって、バッテリ５０から外部回路に供給される電力は、以下の式（６）で表される。
【数6】
【数7】
【0022】
前記した式（５）の動的モデルを用いた場合にも、式（６）と同様の式を用いて放電可能電力を評価できると便利である。そのために、バッテリ５０の実効抵抗Ｒ（ｋ）を、以下の式（８）の通り定義する。
ただし、ηdis（ｋ）は大きさｉdis＞０の一定電流を時刻０からｋまで流したときの過電圧である。この実効抵抗を用いると、以下の式（９）に示すように、近似することができる。
式（９）において、実効抵抗Ｒ（ｋ）は本来、電流ｉ（ｋ）に依存して決まる量であるが、電流ｉ（ｋ）が一定値であるときの式（８）で代用していることに注意が必要である。したがって、この近似は、Ｋが十分に短い等、電流ｉ（ｋ）が一定とみなせる場合に妥当なものとなる。このとき、放電可能電力は、以下の式（１０）に示す通りとなる。
式（５）の動的モデルに対応する実効抵抗Ｒ（ｋ）を求めることを考える。システムの線形性に注意すると，大きさｉdisの一定値電流を流すとき、以下の式（１１）となる。
ただし、ｕstep（ｋ）はｋ＜０で０、ｋ≧０で１の値を取る単位ステップ信号である。式（１１）を式（８）と比較することにより、以下の式（１２）が得られる。
【数8】
【数9】
【数10】
【数11】
【数12】
【0023】
多くのバッテリ５０において、式（１２）の実効抵抗は、時刻ｋに対して単調増加である。区間ｋ∈［０，Ｋ］においてｖoc（ｋ）が一定値であるとみなせるとき、式（１０）の放電可能電力は、以下の式（１３）となる。
このことから、区間ｋ∈［０，Ｋ］の最後の時刻Ｋにおける実効抵抗Ｒ（Ｋ）、すなわちｘ秒実効抵抗を求めることが、放電可能電力の推定に重要であることがわかる。なお、以下の説明では、実効抵抗をＲK：＝Ｒ（Ｋ）とも表記する。
【数13】
【0024】
前記したように、バッテリ５０のＳＯＣｘc（ｋ）は電流ｉ（ｋ）と端子電圧ｖ（ｋ）から推定可能である。また、ＳＯＣ－ＯＣＶ特性ｆocv（・）が既知であるとすると、開回路電圧ｖoc（ｋ）を計算することができ、その結果、過電圧η（ｋ）の値を算出することができる。このとき、実効抵抗Ｒ（Ｋ）の推定は、以下の問題として考えられる。
・問題：各時刻ｋにおいて入出力データｕ（ｋ）＝電流ｉ（ｋ）と、ｙ（ｋ）＝過電圧η（ｋ）が与えられている。この場合に、本実施形態では、与えられたＫに対する実効抵抗ＲKを推定せよ。
・問題：各時刻ｋにおいて入出力データｕ（ｋ）＝電流ｉ（ｋ）と、ｙ（ｋ）＝過電圧η（ｋ）が与えられている。この場合に、本実施形態では、与えられたＫに対する実効抵抗ＲKを推定せよ。
【0025】
［推定方法の従来例］
本実施形態における実効抵抗の推定方法を説明する前に、理解を深めるために、従来の実効抵抗ＲKの推定方法を説明する。
以下では、システム同定分野での慣例に従い、入力をｕ（ｋ）、出力をｙ（ｋ）と記述する。すなわち、ｕ（ｋ）＝電流ｉ（ｋ）であり、ｙ（ｋ）＝過電圧η（ｋ）である。
式（５）に示す内部インピーダンスＧ（ｑ）のインパルス応答をｇ（ｋ）とすると、実効抵抗ＲKは以下の式（１４）で表される。
すなわち、実効抵抗ＲKは、内部インピーダンスＧ（ｑ）のインパルス応答ｇ（ｋ）を計算して、そのＫ番目までの和を算出することで求めることができる。すなわち、ｘ秒実効抵抗は、入力開始時刻からｘ秒が経過した最後の時刻Ｋまでの各時刻における入力に対するインパルス応答の和として表される。
このことから、内部インピーダンスＧ（ｑ）の表現として、インパルス応答を正確かつ効率よく表現できるモデルを用いることが望ましい。この観点から、内部インピーダンスＧ（ｑ）を表す数学モデルの有力な候補のひとつとして、インパルス応答を直接パラメータ化する有限インパルス応答(Finite Impulse Response)モデルがある。有限インパルス応答モデルは、以降「ＦＩＲモデル」と表記する。ＦＩＲモデルは、以下の式（１５）に表される。
ただし、ｂ0，・・・，ｂnはインパルス応答であり、ｗ（ｋ）は式誤差である。また、ｎはＦＩＲモデルの次数であり、インパルス応答が０に収束する時間に応じて選ばれる。式（１５）の右辺の、誤差ｗ（ｋ）を除く各項が、電流ｕ（ｋ）のサンプリング周期（例えば０．１秒）に対応している。
ＦＩＲモデルにおいて、パラメータであるｂ0，・・・，ｂnが推定できると、実効抵抗ＲKはＫ≦ｎに対して、以下の式（１６）により計算することができる。
次数ｎを十分に大きく取ることによって、ＦＩＲモデルは任意の安定な線形システムを記述することができる。一方で、実効抵抗ＲKの計算には、ｂ0，ｂ1，・・・，ｂKのみが用いられるため、式（１５）において推定されたｂK+1，・・・，ｂnは、実効抵抗ＲKの計算には利用されない。前記したように、ｘ秒実効抵抗として、例えば０．１秒実効抵抗や１０秒実効抵抗の推定が求められることが多い。すなわち、実効抵抗ＲKを推定するとき、最後の時刻Ｋがインパルス応答が０に収束する時間に応じて選ばれる次数ｎに比べて、非常に小さいことが多い。すなわち、ＦＩＲモデルを用いて実効抵抗を推定した場合、多くの冗長なパラメータを推定することになる。このことから、推定すべきパラメータを削減したモデルが求められる。
本実施形態における実効抵抗の推定方法を説明する前に、理解を深めるために、従来の実効抵抗ＲKの推定方法を説明する。
以下では、システム同定分野での慣例に従い、入力をｕ（ｋ）、出力をｙ（ｋ）と記述する。すなわち、ｕ（ｋ）＝電流ｉ（ｋ）であり、ｙ（ｋ）＝過電圧η（ｋ）である。
式（５）に示す内部インピーダンスＧ（ｑ）のインパルス応答をｇ（ｋ）とすると、実効抵抗ＲKは以下の式（１４）で表される。
【数14】
このことから、内部インピーダンスＧ（ｑ）の表現として、インパルス応答を正確かつ効率よく表現できるモデルを用いることが望ましい。この観点から、内部インピーダンスＧ（ｑ）を表す数学モデルの有力な候補のひとつとして、インパルス応答を直接パラメータ化する有限インパルス応答(Finite Impulse Response)モデルがある。有限インパルス応答モデルは、以降「ＦＩＲモデル」と表記する。ＦＩＲモデルは、以下の式（１５）に表される。
【数15】
ＦＩＲモデルにおいて、パラメータであるｂ0，・・・，ｂnが推定できると、実効抵抗ＲKはＫ≦ｎに対して、以下の式（１６）により計算することができる。
【数16】
【0026】
そこで、実効抵抗ＲKの推定にμマルコフモデルを用いることが提案されている。μマルコフモデルは、ＦＩＲモデルにＡＲＸモデルを組み合わせたものである。μマルコフモデルは、以下の式（１７）で表される。
ただし、ｂl，（ｌ＝０，１，・・・，Ｋ）は、離散時刻Ｋまでのインパルス応答である。αl，βl，（ｌ＝１，２，・・・，ｐ）は、システムの特性に依存して決まるパラメータである。ｐはパラメトリックモデル部分の次数である。ｗ（ｋ）は式誤差である。
式（１７）の右辺第１項がＦＩＲモデルを表し、右辺第２項および右辺第３項が、ＡＲＸモデルを表している。
このμマルコフモデルは、ＦＩＲモデルのうち、実効抵抗ＲKの計算に不要なｋ＞Ｋの項を有理伝達関数で置き換えたものと解釈することができる。すなわち、μマルコフモデルは、ＦＩＲモデルにおける実効抵抗の推定に不要なパラメータの数を削減することができる。
【数17】
式（１７）の右辺第１項がＦＩＲモデルを表し、右辺第２項および右辺第３項が、ＡＲＸモデルを表している。
このμマルコフモデルは、ＦＩＲモデルのうち、実効抵抗ＲKの計算に不要なｋ＞Ｋの項を有理伝達関数で置き換えたものと解釈することができる。すなわち、μマルコフモデルは、ＦＩＲモデルにおける実効抵抗の推定に不要なパラメータの数を削減することができる。
【0027】
ここで、次数ｐを真のシステムの次数と等しく選ぶと、式（１７）が真のシステムと等価になるパラメータの組が存在する。一方、次数ｐが真の次数と異なる場合は、このようなパラメータの組は存在しない。しかしながら、インパルス応答の推定値ｂlに関しては、ｐの選び方に依らず有効推定値を与えることが示されており、これがμマルコフモデルの優れた性質である。
【0028】
式（１７）において、パラメータベクトルθを、以下の式（１８）とする。
また、回帰ベクトルφ(k)を、以下の式（１９）とする。
この場合、式（１７）は、以下の式（２０）の線形回帰モデルで記述することができる。
このとき、最小二乗法によるパラメータ推定値θｈａｔは、以下の式（２１）と計算される。なお、数式中では、パラメータ推定値θｈａｔは、符号θの上に「＾」を付けて示している。
ただし、Ｎは推定に用いる入出力データ数である。
【数18】
【数19】
【数20】
【数21】
【0029】
前記した式（２０）のμマルコフモデルに基づく線形回帰モデルを逐次推定する方法を説明する。
線形回帰モデルに対する最小二乗法は、以下の式（２２）、式（２３）、式（２４）に示すように逐次化できることが知られている。これを逐次最小二乗法という。
ただし、Ｐ（ｋ）は、以下の式（２５）に表される共分散行列である。
このアルゴリズムは、θ（０）＝０、Ｐ-1（０）＝０のときに、式（２１）の最小二乗法に一致する。実用上は、θｈａｔ（０）を適当な初期値θｈａｔ0とし、共分散行列の初期値Ｐ（０）は適当な正定数γを用いてγＩと置くことが多い。
線形回帰モデルに対する最小二乗法は、以下の式（２２）、式（２３）、式（２４）に示すように逐次化できることが知られている。これを逐次最小二乗法という。
【数22】
【数23】
【0030】
式（２２）～式（２４）の逐次最小二乗法では、時刻ｋ以前のデータに対してすべて同じ重みをつけてパラメータを推定する。しかし、実際の内部インピーダンスは電池のＳＯＣや温度によって時々刻々と動特性が変化する時変システムである。このため、過去になるほど小さな重みをつけて推定を行うことが適切である。このような考え方に基づいた推定法として忘却要素付き逐次最小二乗法が知られている。忘却要素付き逐次最小二乗法における評価関数は、以下の式（２６）で与えられる。
ただし、0＜λ≦1は忘却要素であり、ユーザによって調整されるパラメータである。これに対応する推定法は、以下の式（２７）～式（２９）となる。
ただし、Ｐ（ｋ）は以下の式（３０）に表される。
これは、逐次最小二乗法における共分散行列に対応する量である。明らかに，λ＝１のとき、このアルゴリズムは逐次最小二乗法に一致する。忘却要素λの値を小さくすると，パラメータ変動に対する追従性が向上する。その一方で、実効的に利用できるデータが少なくなり、雑音に対する感度が上がってしまうため、推定値の変動が大きくなる欠点もある。したがって、適切な大きさのλを設定することが、逐次推定において重要である。
【数24】
【数25】
【数26】
【0031】
このように、μマルコフモデルは、ＦＩＲモデルにおけるｋ＞Ｋの項に関するパラメータ数を削減するものである。しかしながら、ｘ秒実効抵抗として、推定が求められる時間が長くなる場合、Ｋが大きな値を取ることになる。このような場合、ＦＩＲモデルのｋ＜Ｋの項で推定するパラメータｂl，（ｌ＝０，１，・・・，Ｋ）の数が多くなる。このような場合は、μマルコフモデルを用いても、十分にパラメータ数を削減できるとは言えない。
長い時間の実効抵抗の推定が要求される場合でもパラメータ数を削減できるように、ｋ＜Ｋの項に対するパラメータ数を削減することが求められている。
長い時間の実効抵抗の推定が要求される場合でもパラメータ数を削減できるように、ｋ＜Ｋの項に対するパラメータ数を削減することが求められている。
【0032】
［本実施形態のモデル構造］
本実施形態では、ＦＩＲモデルを修正したモデルを提案する。本実施形態で用いるモデルは、従来のＦＩＲモデルと区別して、区分定数ＦＩＲ（Piecewise Constant FIR）モデルという。区分定数ＦＩＲモデルを、以降、「ＰＣＦＩＲモデル」という。
ＰＣＦＩＲモデルは、ＦＩＲモデルを構成する項を複数の時間区分で分割し、各時間区分において共通のパラメータを用いるものである。本実施形態では、共通のパラメータの一例として定数を用いているため、区分定数ＦＩＲモデルと称しているが、各時間区分におけるパラメータは共通の関数であれば良く、定数に限定されるものではない。共通の関数として、他にも、１次関数、２次関数、指数関数等を用いることができる。
本実施形態では、ＦＩＲモデルを修正したモデルを提案する。本実施形態で用いるモデルは、従来のＦＩＲモデルと区別して、区分定数ＦＩＲ（Piecewise Constant FIR）モデルという。区分定数ＦＩＲモデルを、以降、「ＰＣＦＩＲモデル」という。
ＰＣＦＩＲモデルは、ＦＩＲモデルを構成する項を複数の時間区分で分割し、各時間区分において共通のパラメータを用いるものである。本実施形態では、共通のパラメータの一例として定数を用いているため、区分定数ＦＩＲモデルと称しているが、各時間区分におけるパラメータは共通の関数であれば良く、定数に限定されるものではない。共通の関数として、他にも、１次関数、２次関数、指数関数等を用いることができる。
【0033】
図６は、本実施形態のＰＣＦＩＲモデルと従来のＦＩＲモデルを示す概念図である。図６では、それぞれのモデルで推定されるインパルス応答、すなわちパラメータを示したものである。実線が本実施形態のＰＣＦＩＲモデルを示し、破線が従来のＦＩＲモデルを示す。
図６に示すように、従来のＦＩＲモデルでは、項ごとに個別のパラメータが推定されるため、インパルス応答は滑らかに下降する曲線を形成する。一方、ＰＣＦＩＲモデルでは、各時間区分に含まれる項のパラメータを定数としているため、インパルス応答は階段状に下降する直線を形成する。図６からわかるように、ＰＣＦＩＲモデルは、従来のＦＩＲモデルと比較して、推定するパラメータ数が削減されている。
図６に示すように、従来のＦＩＲモデルでは、項ごとに個別のパラメータが推定されるため、インパルス応答は滑らかに下降する曲線を形成する。一方、ＰＣＦＩＲモデルでは、各時間区分に含まれる項のパラメータを定数としているため、インパルス応答は階段状に下降する直線を形成する。図６からわかるように、ＰＣＦＩＲモデルは、従来のＦＩＲモデルと比較して、推定するパラメータ数が削減されている。
【0034】
ＰＣＦＩＲモデルは、ＦＩＲモデルの式（１５）を修正した、以下の式（３１）により表される。
ただし、ｂ1，ｂ2，・・・，ｂmはインパルス応答であり、本実施形態で推定されるパラメータである。ｍはパラメータ数であり、すなわち、項を分割する時間区分の数を意味する。ｎ1，ｎ2，・・・，ｎmは、共通のパラメータとする項の数であり、すなわち、各時間区分に含まれる項の数を意味し、言い換えると各時間区分の長さを示すものである。パラメータ数ｍと、項数ｎ1，ｎ2，・・・，ｎmは、それぞれユーザによって指定され、ＰＣＦＩＲモデルの区分情報として、推定装置１のメモリに記憶される。
【数27】
【0035】
図７は、ＰＣＦＩＲモデルにおける時間区分を説明する図である。
図７に示すように、ＰＣＦＩＲモデルにおいては、多項式に含まれる項が、複数の時間区分Ｄ1，Ｄ2・・・Ｄmに分割されている。
最初の時間区分Ｄ1においては、ｎ1の数の項が含まれ、定数ｂ1が各項のパラメータとして用いられている。次の時間区分Ｄ2においては、ｎ2の数の項が含まれ、定数ｂ2が各項のパラメータとして用いられている。最後の時間区分Ｄmにおいては、ｎmの数の項が含まれ、定数ｂmが各項のパラメータとして用いられている。
図７に示すように、ＰＣＦＩＲモデルにおいては、多項式に含まれる項が、複数の時間区分Ｄ1，Ｄ2・・・Ｄmに分割されている。
最初の時間区分Ｄ1においては、ｎ1の数の項が含まれ、定数ｂ1が各項のパラメータとして用いられている。次の時間区分Ｄ2においては、ｎ2の数の項が含まれ、定数ｂ2が各項のパラメータとして用いられている。最後の時間区分Ｄmにおいては、ｎmの数の項が含まれ、定数ｂmが各項のパラメータとして用いられている。
【0036】
式（３１）において，ｍ＝ｎとし、ｎ1＝ｎ2＝・・・＝ｎm＝１とすると、１つの項に１つのパラメータが推定されることになり、ＦＩＲモデルと一致するものとなる。したがって、推定パラメータ数の削減という観点から、ＰＣＦＩＲモデルにおいては、ｍ＜ｎとし、ｎ1，ｎ2，・・・，ｎmは、少なくとも一つを２以上の数とすると良い。すなわち、推定されるパラメータの数ｍが、ＰＣＦＩＲモデルを構成する多項式の項数ｎよりも少なくなることが望ましい。これによって、ＰＣＦＩＲモデルを用いた場合に、従来のＦＩＲモデルと比較して、推定するパラメータの数を減らすことができる。
ｎ1，ｎ2，・・・，ｎmは、例えば、全て同じ数に設定しても良い。この場合、各時間区分に含まれる項の数が一定となる。すなわち、図６に示すように、各時間区分が同じ長さとなる。
また、例えば、ｎ1，ｎ2，・・・，ｎmの少なくとも一つを他と異なる数に設定しても良い。を含むようにしても良い。一例として、ｎ1＜ｎ2＜・・・＜ｎmと、段階的に各時間区分に含まれる項数が大きくなるように設定しても良い。
また、いずれかの時間区分の境界に、推定が求められるｘ秒実効抵抗の最後の時刻Ｋが一致するように、ｍおよびｎ1，ｎ2，・・・，ｎmを設定すると良い。
ｎ1，ｎ2，・・・，ｎmは、例えば、全て同じ数に設定しても良い。この場合、各時間区分に含まれる項の数が一定となる。すなわち、図６に示すように、各時間区分が同じ長さとなる。
また、例えば、ｎ1，ｎ2，・・・，ｎmの少なくとも一つを他と異なる数に設定しても良い。を含むようにしても良い。一例として、ｎ1＜ｎ2＜・・・＜ｎmと、段階的に各時間区分に含まれる項数が大きくなるように設定しても良い。
また、いずれかの時間区分の境界に、推定が求められるｘ秒実効抵抗の最後の時刻Ｋが一致するように、ｍおよびｎ1，ｎ2，・・・，ｎmを設定すると良い。
【0037】
さらに、本実施形態のＰＣＦＩＲモデルは、前記したμマルコフモデルに適用しても良い。すなわち、ＰＣＦＩＲモデルを適用することで、μマルコフモデルの式（１７）は、以下の式（３２）に修正することができる。
式（１７）の従来のμマルコフモデルでは、ＦＩＲモデルで推定されるパラメータｂl，（ｌ＝０，１，・・・，Ｋ）の数は、実効抵抗を求める通電時間Ｋに依存していたが、ＰＣＦＩＲモデルを適用した式（３２）においては、パラメータ数ｍはユーザが任意に設定可能である。
ＰＣＦＩＲモデルをμマルコフモデルに適用することで、ｋ＜Ｋの項についてはＰＣＦＩＲモデルを適用することでパラメータ数が削減され、ｋ＞Ｋの項についてはＡＲＸモデルを適用することでさらにパラメータ数が削減される。
【数28】
ＰＣＦＩＲモデルをμマルコフモデルに適用することで、ｋ＜Ｋの項についてはＰＣＦＩＲモデルを適用することでパラメータ数が削減され、ｋ＞Ｋの項についてはＡＲＸモデルを適用することでさらにパラメータ数が削減される。
【0038】
［ＰＣＦＩＲモデルによる実効抵抗の演算方法］
前記した式（３１）において実効抵抗ＲKを演算する場合、項の数ｎ1，ｎ2，・・・，ｎmを、以下の式（３３）を満たす整数ｍKが存在するように選んだとする。
このとき、ＰＣＦＩＲモデルから算出される実効抵抗ＲKは、以下の式（３４）となる。
すなわち、ｘ秒実効抵抗は、入力開始時刻からｘ秒が経過した最後の時刻Ｋまでの各時間区分で推定されたパラメータｂiに、各時間区分の項数ｎiを掛けた数値（ｎiｂi）を積算することで算出される。
ＰＣＦＩＲモデルの式（３１）において、回帰ベクトルφ（ｋ）およびパラメータベクトルθは、以下の式（３５）および式（３６）と置くことができる。
この場合、式（３１）は、以下の式（３７）の線形回帰モデルで記述することができる。
したがって、本実施形態で用いるＰＣＦＩＲモデルに対しても、μマルコフモデルに適用した一括最小二乗法（式（２１））や逐次最小二乗法（式（２７）～（２９））を適用して、実効抵抗ＲKを算出することができる。
なお、ＰＣＦＩＲモデルをμマルコフモデルに適用したものを用いる場合、式（３２）のＡＲＸモデルの特徴量を式（３５）の回帰ベクトルφ（ｋ）に加え、推定するパラメータベクトルθに、ＡＲＸモデルのパラメータαlとβlを加えれば良い。
前記した式（３１）において実効抵抗ＲKを演算する場合、項の数ｎ1，ｎ2，・・・，ｎmを、以下の式（３３）を満たす整数ｍKが存在するように選んだとする。
【数29】
【数30】
ＰＣＦＩＲモデルの式（３１）において、回帰ベクトルφ（ｋ）およびパラメータベクトルθは、以下の式（３５）および式（３６）と置くことができる。
【数31】
【数32】
なお、ＰＣＦＩＲモデルをμマルコフモデルに適用したものを用いる場合、式（３２）のＡＲＸモデルの特徴量を式（３５）の回帰ベクトルφ（ｋ）に加え、推定するパラメータベクトルθに、ＡＲＸモデルのパラメータαlとβlを加えれば良い。
【0039】
［ＰＣＦＩＲモデルの最小二乗推定値の性質］
ＰＣＦＩＲモデルの最小二乗推定値については、次の性質が成り立つことを証明することができる。
定理１．入出力データは定常であるとし、入力ｕ（ｋ）は白色信号であり、また、入力ｕ（ｋ）と雑音ｗ（ｋ）は無相関であると仮定する。ＰＣＦＩＲモデルから前記した式（３４）によって求めた実効抵抗ＲKは、データ数Ｎが∞のとき、ＦＩＲモデルやμマルコフモデルから式（１６）によって求めた実効抵抗ＲKと一致する。
この定理により、入力信号が白色であり、データ数が十分に多ければ、ＰＣＦＩＲモデルは、従来のＦＩＲモデルまたはμマルコフモデルと等しい実効抵抗ＲKの推定結果を与えることがわかる。
一方、データ数が有限である場合には、パラメータ数ｍとモデルの精度のあいだにはトレードオフの関係がある。パラメータ数ｍを大きくするほど推定精度が上がるが、推定装置１の演算負荷も増加する。そのため、実際の運用に際して、ユーザは、推定装置１の演算負荷と、要求される推定精度のバランスを考慮して、パラメータ数ｍを設定する必要がある。
ＰＣＦＩＲモデルの最小二乗推定値については、次の性質が成り立つことを証明することができる。
定理１．入出力データは定常であるとし、入力ｕ（ｋ）は白色信号であり、また、入力ｕ（ｋ）と雑音ｗ（ｋ）は無相関であると仮定する。ＰＣＦＩＲモデルから前記した式（３４）によって求めた実効抵抗ＲKは、データ数Ｎが∞のとき、ＦＩＲモデルやμマルコフモデルから式（１６）によって求めた実効抵抗ＲKと一致する。
この定理により、入力信号が白色であり、データ数が十分に多ければ、ＰＣＦＩＲモデルは、従来のＦＩＲモデルまたはμマルコフモデルと等しい実効抵抗ＲKの推定結果を与えることがわかる。
一方、データ数が有限である場合には、パラメータ数ｍとモデルの精度のあいだにはトレードオフの関係がある。パラメータ数ｍを大きくするほど推定精度が上がるが、推定装置１の演算負荷も増加する。そのため、実際の運用に際して、ユーザは、推定装置１の演算負荷と、要求される推定精度のバランスを考慮して、パラメータ数ｍを設定する必要がある。
【0040】
ここで、入力ｕ（ｋ）が白色信号ではなく、以下の式（３８）に示すように、整形フィルタＨ（ｑ）によって整形された有色信号であるとする。
ただし、ｖ（ｋ）は白色雑音である。
この場合は、前記の定理１の仮定が成り立たないため、有色信号の入力ｕ（ｋ）をそのまま用いて、ＰＣＦＩＲモデルにより実効抵抗を推定することができない。
【数33】
この場合は、前記の定理１の仮定が成り立たないため、有色信号の入力ｕ（ｋ）をそのまま用いて、ＰＣＦＩＲモデルにより実効抵抗を推定することができない。
【0041】
そこで、本実施形態では、入力ｕ（ｋ）と出力ｙ（ｋ）に対して、白色化フィルタＨ-1（ｑ）を用いたフィルタリング処理を行うことで、入力ｕ（ｋ）を白色化させる。白色化フィルタＨ-1（ｑ）は、前記した整形フィルタＨ-1（ｑ）の逆フィルタを意味するものである。
フィルタリング処理は、以下の式（３９）および式（４０）で表される。
このように、Ｈ（ｑ）の逆フィルタである白色化フィルタＨ-1（ｑ）を通した入出力ｕf（ｋ）とｙf（ｋ）を用いることにより、入力ｕ（ｋ）が白色信号であるときの同定問題に帰着することができるため、正しく実効抵抗ＲKを求めることができる。白色化フィルタＨ-1（ｑ）が既知である場合は、予め設定して推定装置１のメモリに格納し、フィルタリング処理に用いることができる。
フィルタリング処理は、以下の式（３９）および式（４０）で表される。
【数34】
【0042】
［実効抵抗の推定処理］
図８は、本実施形態における推定装置１の処理を示すフローチャートである。
図８に示すように、推定装置１の過電圧算出部２０は、電圧センサ６０および電流センサ７０によって測定された、端子電圧ｖ（ｋ）および電流ｉ（ｋ）を取得する（ステップＳ０１）。
過電圧算出部２０のＳＯＣ算出部２１は、例えば、式（２）の電流積算法や、カルマンフィルタ等を用いて、バッテリ５０のＳＯＣｘｃ（ｋ）を算出する（ステップＳ０２）。
ＳＯＣ－ＯＣＶ変換部２２は、メモリに格納されている既知のＳＯＣ－ＯＣＶ特性（式（３）参照）に基づいて、ＳＯＣｘｃ（ｋ）を、開回路電圧ｖoc（ｋ）に変換する（ステップＳ０３）。
演算器２３は、開回路電圧ｖoc（ｋ）と、端子電圧ｖ（ｋ）との差を演算して、過電圧η（ｋ）として推定部３０に出力する（ステップＳ０４）。
図８は、本実施形態における推定装置１の処理を示すフローチャートである。
図８に示すように、推定装置１の過電圧算出部２０は、電圧センサ６０および電流センサ７０によって測定された、端子電圧ｖ（ｋ）および電流ｉ（ｋ）を取得する（ステップＳ０１）。
過電圧算出部２０のＳＯＣ算出部２１は、例えば、式（２）の電流積算法や、カルマンフィルタ等を用いて、バッテリ５０のＳＯＣｘｃ（ｋ）を算出する（ステップＳ０２）。
ＳＯＣ－ＯＣＶ変換部２２は、メモリに格納されている既知のＳＯＣ－ＯＣＶ特性（式（３）参照）に基づいて、ＳＯＣｘｃ（ｋ）を、開回路電圧ｖoc（ｋ）に変換する（ステップＳ０３）。
演算器２３は、開回路電圧ｖoc（ｋ）と、端子電圧ｖ（ｋ）との差を演算して、過電圧η（ｋ）として推定部３０に出力する（ステップＳ０４）。
【0043】
推定部３０では、電流センサ７０によって測定された電流ｉ（ｋ）と、過電圧算出部２０で算出された過電圧η（ｋ）を時系列データとして取得し、電流ｉ（ｋ）を入力ｕ（ｋ）、過電圧η（ｋ）を出力ｙ（ｋ）とする（ステップＳ０５）。
推定部３０のフィルタリング部３１は、メモリから白色化フィルタを取得し、入力ｕ（ｋ）と出力ｙ（ｋ）に対して、式（３９）および式（４０）に示すフィルタリング処理を行う（ステップＳ０６）。
入力ｕ（ｋ）が有色信号である場合に、このフィルタリング処理によって白色信号とすることができる。フィルタリング部３１は、ｕf（ｋ）とｙf（ｋ）の時系列データを特徴ベクトル算出部３２に出力する。
推定部３０のフィルタリング部３１は、メモリから白色化フィルタを取得し、入力ｕ（ｋ）と出力ｙ（ｋ）に対して、式（３９）および式（４０）に示すフィルタリング処理を行う（ステップＳ０６）。
入力ｕ（ｋ）が有色信号である場合に、このフィルタリング処理によって白色信号とすることができる。フィルタリング部３１は、ｕf（ｋ）とｙf（ｋ）の時系列データを特徴ベクトル算出部３２に出力する。
【0044】
特徴ベクトル算出部３２は、フィルタリング処理されたｕf（ｋ）とｙf（ｋ）を、それぞれ入力ｕ（ｋ）および出力ｙ（ｋ）とする（ステップＳ０７）。特徴ベクトル算出部３２は、推定装置１のメモリから、ユーザによって指定された、ＰＣＦＩＲモデルの区分情報であるパラメータ数ｍおよび項数ｎ1，ｎ2，・・・，ｎmを取得する。特徴ベクトル算出部３２は、入力ｕ（ｋ）および出力ｙ（ｋ）を用いてＰＣＦＩＲモデルの式（３５）を構成し、特徴ベクトルである回帰ベクトルφを算出する（ステップＳ０８）。
【0045】
パラメータ推定部３３は、回帰ベクトルφを用いて、一括最小二乗法（式（２１））や逐次最小二乗法（式（２７）～（２９））により、式（３６）に示すパラメータベクトルθを推定する（ステップＳ０９）。
ステップＳ０８、Ｓ０９の特徴ベクトル算出部３２およびパラメータ推定部３３の処理によって、ＰＣＦＩＲモデルを適用した式（３１）または式（３２）によるバッテリシステムが同定される。
実効抵抗算出部３４は、推定されたパラメータベクトルθを用いて、式（３４）により、実効抵抗ＲKを算出する（ステップＳ１０）。
ステップＳ０８、Ｓ０９の特徴ベクトル算出部３２およびパラメータ推定部３３の処理によって、ＰＣＦＩＲモデルを適用した式（３１）または式（３２）によるバッテリシステムが同定される。
実効抵抗算出部３４は、推定されたパラメータベクトルθを用いて、式（３４）により、実効抵抗ＲKを算出する（ステップＳ１０）。
【0046】
［数値シミュレーション］
本実施形態のＰＣＦＩＲモデルの有用性を、数値シミュレーションを通して説明する。
図９は、数値シミュレーションに用いた回路構成を示す図である。
図１０の（ａ）は、数値シミュレーションにおける入力データを示す図であり、図１０の（ｂ）は出力データを示す図である。
図１１は、本実施形態のＰＣＦＩＲモデルと従来のＦＩＲモデルによる実効抵抗の推定結果を示すグラフである。図１１の（ａ）は、各モデルにおけるインパルス応答である電圧（ｖ）を示す図であり、各モデルで推定されるパラメータに対応するものである。図１１の（ｂ）は、各モデルにおけるステップ応答である実効抵抗を示す図である。また、実線がＰＣＦＩＲモデルを示し、破線が従来のＦＩＲモデルを示す。
本実施形態のＰＣＦＩＲモデルの有用性を、数値シミュレーションを通して説明する。
図９は、数値シミュレーションに用いた回路構成を示す図である。
図１０の（ａ）は、数値シミュレーションにおける入力データを示す図であり、図１０の（ｂ）は出力データを示す図である。
図１１は、本実施形態のＰＣＦＩＲモデルと従来のＦＩＲモデルによる実効抵抗の推定結果を示すグラフである。図１１の（ａ）は、各モデルにおけるインパルス応答である電圧（ｖ）を示す図であり、各モデルで推定されるパラメータに対応するものである。図１１の（ｂ）は、各モデルにおけるステップ応答である実効抵抗を示す図である。また、実線がＰＣＦＩＲモデルを示し、破線が従来のＦＩＲモデルを示す。
【0047】
図９は、図４のバッテリシステムの等価回路の内部インピーダンスの動特性を、ＲＣ並列回路で表現したものである。図９において、Ｒ0，…，Ｒ4は抵抗値、Ｃ1，…，Ｃ4はコンデンサの静電容量である。
図９の回路に対して、以下の式（４１）に示す整形フィルタで整形された有色信号である入力ｕ（ｋ）を印加し、雑音ｗ（ｋ）を加えて、出力ｙ（ｋ）を得た。
数値シミュレーションでは０．１秒周期でサンプリングを行った。図１０は、数値シミュレーションで得られた入出力データの一部を示すものである。図１０の（ａ）が入力である電流ｕ（ｋ）を示すデータであり、図１０の（ｂ）が出力である過電圧ｙ（ｋ）を示すデータである。
図１０の入出力データを用いて、本実施形態のＰＣＦＩＲモデルと、従来のＦＩＲモデルとのそれぞれによって実効抵抗の推定を行い、図１１に示す推定結果を得た。
図９の回路に対して、以下の式（４１）に示す整形フィルタで整形された有色信号である入力ｕ（ｋ）を印加し、雑音ｗ（ｋ）を加えて、出力ｙ（ｋ）を得た。
【数35】
図１０の入出力データを用いて、本実施形態のＰＣＦＩＲモデルと、従来のＦＩＲモデルとのそれぞれによって実効抵抗の推定を行い、図１１に示す推定結果を得た。
【0048】
従来のＦＩＲモデルでは、パラメータ数が８００であるのに対し、本実施形態におけるＰＣＦＩＲモデルでは、パラメータ数は４であり、大幅にパラメータ数が削減されていることがわかる。
さらに図１１の（ｂ）に示すように、ＰＣＦＩＲモデルは、各時間区分内の実効抵抗の数値はＦＩＲモデルに対してずれが生じているが、各時間区分の境界では、ＰＣＦＩＲモデルの実効抵抗はＦＩＲモデルの実効抵抗と一致している。
さらに図１１の（ｂ）に示すように、ＰＣＦＩＲモデルは、各時間区分内の実効抵抗の数値はＦＩＲモデルに対してずれが生じているが、各時間区分の境界では、ＰＣＦＩＲモデルの実効抵抗はＦＩＲモデルの実効抵抗と一致している。
【0049】
実効抵抗ＲKの推定においては、推定したいｘ秒実効抵抗の精度が確保できれば良い。そのため、ｘ秒実効抵抗の最後の時刻Ｋを、時間区分の境界に一致させることで、ＰＣＦＩＲモデルを用いた場合の実効抵抗の推定精度を確保することができる。具体的には、式（３１）または式（３２）における区分数ｍおよび各区分の項数ｎ1，ｎ2，・・・，ｎmを調整することで、時間区分の境界と最後の時刻Ｋを一致させることができる。
【0050】
図１１は、データ数Ｎが十分に多い場合の推定結果である。すなわち、前記した定理１．で説明したように、データ数Ｎが十分に多い場合には、ＰＣＦＩＲモデルは、ＦＩＲモデルと等しい実効抵抗の推定結果を得ることができる。
【0051】
ここで、データ数Ｎが少ない場合についても検討する。
図１２および図１３は、データ数Ｎが少ない場合の推定結果である。
図１２は各時間区分の長さを一定（ｎ1＝ｎ2・・・＝ｎm）とした場合の推定結果を示している。図１２の（ａ）は、インパルス応答であり、図１２の（ｂ）は、ステップ応答である実効抵抗を示している。
図１３は、各時間区分の長さを段階的に長くした場合（ｎ1＜ｎ2＜・・・＜ｎm）の推定結果である。図１３の（ａ）は、インパルス応答であり、図１３の（ｂ）は、ステップ応答である実効抵抗を示している。
図１２および図１３では、ＰＣＦＩＲモデルにより、雑音ｗ（ｋ）を変化させながら複数回推定を行っている。この場合、推定結果にはバラつきが出る。図１２および図１３で示した２本の点線の間に、複数回の推定結果の９５％が入っている。すなわち、２本の点線は、推定結果のバラつきの大きさを示している。
また、図１２および図１３の実線は、複数回の推定結果の平均値である。一点鎖線は、ＰＣＦＩＲモデルによる推定の正解値であり、図１１の実線で示した値に対応する。破線は、従来のＦＩＲモデルによる推定結果である。
図１２および図１３は、データ数Ｎが少ない場合の推定結果である。
図１２は各時間区分の長さを一定（ｎ1＝ｎ2・・・＝ｎm）とした場合の推定結果を示している。図１２の（ａ）は、インパルス応答であり、図１２の（ｂ）は、ステップ応答である実効抵抗を示している。
図１３は、各時間区分の長さを段階的に長くした場合（ｎ1＜ｎ2＜・・・＜ｎm）の推定結果である。図１３の（ａ）は、インパルス応答であり、図１３の（ｂ）は、ステップ応答である実効抵抗を示している。
図１２および図１３では、ＰＣＦＩＲモデルにより、雑音ｗ（ｋ）を変化させながら複数回推定を行っている。この場合、推定結果にはバラつきが出る。図１２および図１３で示した２本の点線の間に、複数回の推定結果の９５％が入っている。すなわち、２本の点線は、推定結果のバラつきの大きさを示している。
また、図１２および図１３の実線は、複数回の推定結果の平均値である。一点鎖線は、ＰＣＦＩＲモデルによる推定の正解値であり、図１１の実線で示した値に対応する。破線は、従来のＦＩＲモデルによる推定結果である。
【0052】
図１２および図１３を比較すると、各時間区分の長さを段階的に長くした図１３の方が、２本の破線の間の範囲が狭く、推定結果のバラつきが小さいことがわかる。また、図１２、図１３のいずれの場合も、推定結果の平均（実線）を取ることで正解値（一点鎖線）に近づくが、図１２の（ｂ）では、時間が経過するごとに、平均値（実線）と正解値（一点鎖線）の乖離が大きくなっている。一方、図１３の（ｂ）では、平均値（実線）と正解値（一点鎖線）の乖離が小さく抑えられている。
図１２の（ａ）および図１３の（ａ）のＦＩＲモデル（破線）からわかるように、インパルス応答は、入力開始時刻（時間０）に近い方が変化が大きく、入力開始時刻から離れるにつれて、変化が収束していく。そのため、図１３のように各時間区分の長さを段階的に長くした場合、インパルス応答の変化が大きい箇所において推定されるパラメータが多くなるため、推定精度を向上させることができる。
このように、本実施形態におけるＰＣＦＩＲモデルは、各時間区分の項数ｎiを任意に設定することができるため、例えばデータ数Ｎが少ないような場合でも、各時間区分の長さを調整して、推定精度の向上に寄与することができる。
図１２の（ａ）および図１３の（ａ）のＦＩＲモデル（破線）からわかるように、インパルス応答は、入力開始時刻（時間０）に近い方が変化が大きく、入力開始時刻から離れるにつれて、変化が収束していく。そのため、図１３のように各時間区分の長さを段階的に長くした場合、インパルス応答の変化が大きい箇所において推定されるパラメータが多くなるため、推定精度を向上させることができる。
このように、本実施形態におけるＰＣＦＩＲモデルは、各時間区分の項数ｎiを任意に設定することができるため、例えばデータ数Ｎが少ないような場合でも、各時間区分の長さを調整して、推定精度の向上に寄与することができる。
【0053】
以上の通り、本実施形態の推定装置１は、以下の構成を有する。
（１）推定装置１は、バッテリ５０の所定時間（ｘ秒）における過電圧変化を示す実効抵抗を推定する。
推定装置１は、過電圧算出部２０と、特徴ベクトル算出部３２およびパラメータ推定部３３（システム同定部）と、実効抵抗算出部３４と、を備える。
過電圧算出部２０は、バッテリ５０を流れる電流ｉ（ｋ）と、バッテリ５０の端子電圧ｖ（ｋ）の測定値に基づいて、バッテリ５０の過電圧η（ｋ）を算出する。
特徴ベクトル算出部３２およびパラメータ推定部３３は、電流と過電圧の時系列データ（入力ｕ（ｋ）、出力ｙ（ｋ））に基づいて、ＦＩＲモデルを含むモデルを構成する多項式の各項のパラメータｂ1，ｂ2，・・・，ｂm（係数）を推定することで、バッテリ５０のシステムを同定する。
実効抵抗算出部３４は、バッテリ５０のシステムに基づいて、実効抵抗ＲKを算出する。
特徴ベクトル算出部３２およびパラメータ推定部３３は、ＦＩＲモデルを構成する多項式を複数の時間区分で分割し、各時間区分に含まれる項に対して共通のパラメータを推定することで、ＦＩＲモデルの多項式の項数よりも少ない数のパラメータを推定する。
（１）推定装置１は、バッテリ５０の所定時間（ｘ秒）における過電圧変化を示す実効抵抗を推定する。
推定装置１は、過電圧算出部２０と、特徴ベクトル算出部３２およびパラメータ推定部３３（システム同定部）と、実効抵抗算出部３４と、を備える。
過電圧算出部２０は、バッテリ５０を流れる電流ｉ（ｋ）と、バッテリ５０の端子電圧ｖ（ｋ）の測定値に基づいて、バッテリ５０の過電圧η（ｋ）を算出する。
特徴ベクトル算出部３２およびパラメータ推定部３３は、電流と過電圧の時系列データ（入力ｕ（ｋ）、出力ｙ（ｋ））に基づいて、ＦＩＲモデルを含むモデルを構成する多項式の各項のパラメータｂ1，ｂ2，・・・，ｂm（係数）を推定することで、バッテリ５０のシステムを同定する。
実効抵抗算出部３４は、バッテリ５０のシステムに基づいて、実効抵抗ＲKを算出する。
特徴ベクトル算出部３２およびパラメータ推定部３３は、ＦＩＲモデルを構成する多項式を複数の時間区分で分割し、各時間区分に含まれる項に対して共通のパラメータを推定することで、ＦＩＲモデルの多項式の項数よりも少ない数のパラメータを推定する。
【0054】
バッテリ５０の実効抵抗ＲKは、バッテリ５０の内部インピーダンスの特性を表現する特徴量であり、入力開始時刻からｘ秒が経過した最後の時刻Ｋまでの各時刻における入力に対するインパルス応答の和として表される。
推定装置１は、インパルス応答をパラメータとしてバッテリ５０のシステムを表現するＰＣＦＩＲモデルを用いて、実効抵抗ＲKを推定する。特徴ベクトル算出部３２およびパラメータ推定部３３は、バッテリ５０の電流と過電圧の時系列データｕ（ｋ）、ｙ（ｋ）に基づいてパラメータｂ1，ｂ2，・・・，ｂmを推定してバッテリ５０のシステムを同定することで、実効抵抗ＲKを算出することができる。
ここで、従来のＦＩＲモデルによりバッテリ５０のシステムを同定する場合、多項式を構成する項の数ｎと同数のパラメータｂ0，・・・，ｂnを推定する必要がある。しかしながら、実効抵抗の推定に必要なのは、入力開始時刻から最後の時刻Ｋまでのインパルス応答ｂ0，ｂ1，・・・，ｂKのみである。すなわち、従来のＦＩＲモデルでは、実効抵抗ＲKの推定に不要な多くのパラメータｂK+1，・・・，ｂnを推定することになる。
従来のＦＩＲモデルとＡＲＸモデルを組み合わせたμマルコフモデルは、実効抵抗ＲKの推定に不要な項（ｋ＞Ｋの項）について推定するパラメータの数を削減することができる。しかしながら、推定する実効抵抗ＲKの時間（ｘ秒）が大きくなると、実効抵抗ＲKの推定に必要なｋ＜Ｋの項の数が多くなるため、十分に推定するパラメータの数を削減することができない。
推定装置１は、インパルス応答をパラメータとしてバッテリ５０のシステムを表現するＰＣＦＩＲモデルを用いて、実効抵抗ＲKを推定する。特徴ベクトル算出部３２およびパラメータ推定部３３は、バッテリ５０の電流と過電圧の時系列データｕ（ｋ）、ｙ（ｋ）に基づいてパラメータｂ1，ｂ2，・・・，ｂmを推定してバッテリ５０のシステムを同定することで、実効抵抗ＲKを算出することができる。
ここで、従来のＦＩＲモデルによりバッテリ５０のシステムを同定する場合、多項式を構成する項の数ｎと同数のパラメータｂ0，・・・，ｂnを推定する必要がある。しかしながら、実効抵抗の推定に必要なのは、入力開始時刻から最後の時刻Ｋまでのインパルス応答ｂ0，ｂ1，・・・，ｂKのみである。すなわち、従来のＦＩＲモデルでは、実効抵抗ＲKの推定に不要な多くのパラメータｂK+1，・・・，ｂnを推定することになる。
従来のＦＩＲモデルとＡＲＸモデルを組み合わせたμマルコフモデルは、実効抵抗ＲKの推定に不要な項（ｋ＞Ｋの項）について推定するパラメータの数を削減することができる。しかしながら、推定する実効抵抗ＲKの時間（ｘ秒）が大きくなると、実効抵抗ＲKの推定に必要なｋ＜Ｋの項の数が多くなるため、十分に推定するパラメータの数を削減することができない。
【0055】
本実施形態では、システム同定部である特徴ベクトル算出部３２およびパラメータ推定部３３ＦＩＲモデルは、従来のＦＩＲモデルを修正したＰＣＦＩＲモデルを用いている。ＰＣＦＩＲモデルでは、ＦＩＲモデルを構成する多項式を複数の時間区分で分割し、各時間区分に含まれる項に対して共通のパラメータを推定する。これによって、推定するパラメータの数が、ＰＣＦＩＲモデルを構成する多項式の項数よりも少なくなる（ｍ＜ｎ）ため、推定する実効抵抗ＲKの秒数（ｘ秒）が大きい場合でも、十分に推定するパラメータの数を削減することができる。これによって、推定装置１の演算負荷を低減することができる。
【0056】
（２）推定装置１は、フィルタリング部３１を備える。
フィルタリング部３１は、電流と過電圧の時系列データ（ｕ（ｋ）、ｙ（ｋ））に対して、白色化フィルタＨ-1（ｑ）によりフィルタリング処理を行い、電流が白色信号である時系列データ（ｕf（ｋ）、ｙf（ｋ））を、システム同定部である特徴ベクトル算出部３２に出力する。
フィルタリング部３１は、電流と過電圧の時系列データ（ｕ（ｋ）、ｙ（ｋ））に対して、白色化フィルタＨ-1（ｑ）によりフィルタリング処理を行い、電流が白色信号である時系列データ（ｕf（ｋ）、ｙf（ｋ））を、システム同定部である特徴ベクトル算出部３２に出力する。
【0057】
ＰＣＦＩＲモデルの入力である電流ｕ（ｋ）が有色信号である場合には、前記したＰＣＦＩＲモデルについての定理１．が成立しないため、実効抵抗ＲKの推定精度に影響を与える可能性がある。本実施形態では、フィルタリング部３１が白色化フィルタＨ-1（ｑ）によりフィルタリング処理を行う。電流を白色信号とすることで、定理１．が成立する状態となり、ＰＣＦＩＲモデルを用いた精度の高い推定を行うことができる。
なお、入力の電流ｕ（ｋ）として白色信号を用いることができる場合には、フィルタリング部３１は省略しても良い。
なお、入力の電流ｕ（ｋ）として白色信号を用いることができる場合には、フィルタリング部３１は省略しても良い。
【0058】
（３）実効抵抗算出部３４は、時系列データの入力開始時刻から所定時間（ｘ秒）経過後の最後の時刻Ｋまでの、各時間区分で推定されたパラメータｂiに、各時間区分の項数ｎiを掛けた数値（ｎiｂi）を積算することで、実効抵抗ＲKを算出する。
特徴ベクトル算出部３２およびパラメータ推定部３３（システム同定部）において、時間区分のいずれかの境界と最後の時刻Ｋが一致するように、時間区分の長さが設定される。
特徴ベクトル算出部３２およびパラメータ推定部３３（システム同定部）において、時間区分のいずれかの境界と最後の時刻Ｋが一致するように、時間区分の長さが設定される。
【0059】
ＰＣＦＩＲモデルによる推定結果を、従来のＦＩＲモデルによる推定結果と比較した場合（図１１参照）、各時間区分内の実効抵抗の数値はずれが生じるが、各時間区分の境界の数値は一致する傾向がある。そのため、推定装置１では、時間区分の境界と最後の時刻Ｋが一致するように時間区分の長さを設定することで、ｘ秒実効抵抗の推定精度を高めることができる。
【0060】
（４）特徴ベクトル算出部３２およびパラメータ推定部３３（システム同定部）において、時系列データにおける入力開始時刻に近い時間区分ほど長さが短くなるように、各時間区分の長さが設定される。
【0061】
ＰＣＦＩＲモデルのパラメータによって示されるインパルス応答は、入力開始時刻に近いほど変化が大きい。そこで、推定装置１では、入力開始時刻に近いほど時間区分が長くなるように各時間区分の長さを設定する。これによって、インパルス応答の変化が大きい箇所において推定されるパラメータが多くなるため、実効抵抗ＲKの推定精度を向上させることができる。
【0062】
（５）特徴ベクトル算出部３２およびパラメータ推定部３３（システム同定部）は、各時間区分に含まれる項に対して、共通のパラメータとして定数を推定する。
【0063】
各時間区分で推定するパラメータを定数とすることで、各時間区分で１つのパラメータのみが推定対象となる。そのため、ＰＣＦＩＲモデルにおいて推定されるパラメータ数を大きく削減することができる。
【0064】
（６）特徴ベクトル算出部３２およびパラメータ推定部３３（システム同定部）は、ＦＩＲモデルを含むモデルとして、ＦＩＲモデルにＡＲＸモデルを組み合わせたμマルコフモデルを適用し、時系列データの入力開始時刻から所定時間経過後の最後の時刻Ｋまでの時系列データをＦＩＲモデルに適用する。
【0065】
μマルコフモデルにＰＣＦＩＲモデルを適用することで、多項式のｋ＞Ｋの項にはＰＣＦＩＲモデルが適用され、ｋ＜Ｋの項にＡＲＸモデルが適用されるため、推定するパラメータ数をさらに削減することができる。
【0066】
［変形例１］
前記した実施形態では、ＰＣＦＩＲモデルにおいて、各時間区分における共通するパラメータとして定数を用いる例を説明したが、この態様に限定されない。前記したように、共通するパラメータは共通する関数であれば良く、例えば、１次関数、２次関数、指数関数等を用いることができる。変形例１では、共通するパラメータとして１次関数を用いた区分線形ＦＩＲモデルを説明する。
図１４は、変形例１に係る区分線形ＦＩＲモデルと、従来のＦＩＲモデルの概念図である。
図１４に示すように、区分線形モデルでは、各時間区分に含まれる項のパラメータを１次関数としているため、各時間区分において推定されるパラメータは、下向きに傾斜する直線を形成し、さらに時間区分ごとに異なる傾きを示す。図１４からわかるように、区分線形ＦＩＲモデルも、従来のＦＩＲモデルと比較して、推定するパラメータ数が大幅に削減されている。
区分線形モデルは、ＦＩＲモデルの式（１５）を修正した、以下の式（４１）により表される。
ただし、ｎ1，ｎ2・・・，ｎmは、ユーザが指定する項数であり、β11，β12，・・・，βm1，βm2が推定されるパラメータである。すなわち、変形例１では、一次関数の傾きと切片を構成するパラメータを推定することになる。例えば、最初の時間区分において、最初の時刻ｕ（ｋ）のパラメータがβ11であり、最後の時刻ｕ（ｋ－ｎ1＋１）の時刻のパラメータがβ12とであり、最初の時刻と最後の時刻の間のパラメータは、β11とβ12の内分点として設定される。すなわち、各時間区分において、２つのパラメータが推定される。
前記した実施形態では、ＰＣＦＩＲモデルにおいて、各時間区分における共通するパラメータとして定数を用いる例を説明したが、この態様に限定されない。前記したように、共通するパラメータは共通する関数であれば良く、例えば、１次関数、２次関数、指数関数等を用いることができる。変形例１では、共通するパラメータとして１次関数を用いた区分線形ＦＩＲモデルを説明する。
図１４は、変形例１に係る区分線形ＦＩＲモデルと、従来のＦＩＲモデルの概念図である。
図１４に示すように、区分線形モデルでは、各時間区分に含まれる項のパラメータを１次関数としているため、各時間区分において推定されるパラメータは、下向きに傾斜する直線を形成し、さらに時間区分ごとに異なる傾きを示す。図１４からわかるように、区分線形ＦＩＲモデルも、従来のＦＩＲモデルと比較して、推定するパラメータ数が大幅に削減されている。
区分線形モデルは、ＦＩＲモデルの式（１５）を修正した、以下の式（４１）により表される。
【数36】
【0067】
区分線形ＦＩＲモデルの式（４２）において、回帰ベクトルφ（ｋ）およびパラメータベクトルθは、それぞれ以下の式（４３）および式（４４）により表される。
区分線形ＦＩＲモデルの式（４２）は、ＰＣＦＩＲモデルと同様に、式（３７）の線形回帰モデルで記述することができる。
したがって、変形例１の区分線形ＦＩＲデルに対しても、一括最小二乗法（式（２１））や逐次最小二乗法（式（２７）～（２９））を適用して、実効抵抗ＲKを算出することができる。
なお、変形例１の区分線形ＦＩＲモデルは、前記したμマルコフモデルに適用しても良い。
【数37】
したがって、変形例１の区分線形ＦＩＲデルに対しても、一括最小二乗法（式（２１））や逐次最小二乗法（式（２７）～（２９））を適用して、実効抵抗ＲKを算出することができる。
なお、変形例１の区分線形ＦＩＲモデルは、前記したμマルコフモデルに適用しても良い。
【0068】
以上のように、変形例１に係る推定装置１は、以下の構成を備える。
（ｉ）特徴ベクトル算出部３２とパラメータ推定部３３（システム同定部）は、ＦＩＲモデルを構成する多項式を複数の時間区分で分割し、各時間区分に含まれる項に対して、共通の係数として、一次関数を推定する。
（ｉ）特徴ベクトル算出部３２とパラメータ推定部３３（システム同定部）は、ＦＩＲモデルを構成する多項式を複数の時間区分で分割し、各時間区分に含まれる項に対して、共通の係数として、一次関数を推定する。
【0069】
ＰＣＦＩＲモデルと同様に、各時間区分に含まれる項において共通の関数を用いることで、全体で推定されるパラメータ数を削減することができる。これによって、推定装置１の演算負荷を低減することができる。さらに、係数を一次関数とすることで、各時間区分で推定されるパラメータが２倍となる分、実効抵抗の推定精度を向上させることができる。
【0070】
［変形例２］
図１５は、変形例２に係る推定部３０Ａの構成を示すブロック図である。
前記した実施形態では、白色化フィルタが既知である場合の例を説明した。この場合、推定部３０のフィルタリング部３１（図４参照）は、予め設定された白色化フィルタをメモリから取得して、フィルタリング処理を行う。
変形例２では、白色化フィルタが未知である場合を説明する。変形例２では、実効抵抗の推定を行う際に、電流センサ７０で測定される電流ｉ（ｋ）の時系列データに対して時系列解析を行うことで、フィルタリング処理で用いる白色化フィルタＨ-1（ｑ）を作成する。
図１５は、変形例２に係る推定部３０Ａの構成を示すブロック図である。
前記した実施形態では、白色化フィルタが既知である場合の例を説明した。この場合、推定部３０のフィルタリング部３１（図４参照）は、予め設定された白色化フィルタをメモリから取得して、フィルタリング処理を行う。
変形例２では、白色化フィルタが未知である場合を説明する。変形例２では、実効抵抗の推定を行う際に、電流センサ７０で測定される電流ｉ（ｋ）の時系列データに対して時系列解析を行うことで、フィルタリング処理で用いる白色化フィルタＨ-1（ｑ）を作成する。
【0071】
図１５に示すように、変形例２において、推定部３０Ａは、フィルタリング部３１の前段に、時系列解析部３５と、フィルタ作成部３６と、を更に備える。
時系列解析部３５は、電流センサ７０で測定される電流ｉ（ｋ）の時系列データを用いて、既知の時系列解析手法により、式（３８）に定義される整形フィルタＨ（ｑ）を推定する。時系列解析に用いるモデルは限定されるものではないが、例えば、移動平均モデル（moving average model：ＭＡモデル）、自己回帰モデル（autoregressive model：ＡＲモデル）等を用いることができる。
時系列解析部３５は、電流センサ７０で測定される電流ｉ（ｋ）の時系列データを用いて、既知の時系列解析手法により、式（３８）に定義される整形フィルタＨ（ｑ）を推定する。時系列解析に用いるモデルは限定されるものではないが、例えば、移動平均モデル（moving average model：ＭＡモデル）、自己回帰モデル（autoregressive model：ＡＲモデル）等を用いることができる。
【0072】
フィルタ作成部３６は、時系列解析部３５で推定された整形フィルタＨ（ｑ）の逆フィルタである白色化フィルタＨ-1（ｑ）を作成する。フィルタ作成部３６は、作成した白色化フィルタＨ-1（ｑ）をフィルタリング部３１に出力する。フィルタリング部３１は、前記した実施形態と同様に、白色化フィルタＨ-1（ｑ）を用いて、入力ｕ（ｋ）と出力ｙ（ｋ）に対して、式（３９）および式（４０）によるフィルタリング処理を行う。
【0073】
以上のように、変形例２に係る推定装置１は、以下の構成を備える。
（ｉｉ）推定装置１は、時系列解析部３５と、フィルタ作成部３６と、を備える。
時系列解析部３５は、バッテリ５０を流れる電流ｉ（ｋ）の時系列データを時系列解析することで、電流の整形フィルタＨ（ｑ）を推定する。
フィルタ作成部３６は、整形フィルタＨ（ｑ）の逆フィルタである、白色化フィルタＨ-1（ｑ）を作成し、フィルタリング部３１に出力する。
（ｉｉ）推定装置１は、時系列解析部３５と、フィルタ作成部３６と、を備える。
時系列解析部３５は、バッテリ５０を流れる電流ｉ（ｋ）の時系列データを時系列解析することで、電流の整形フィルタＨ（ｑ）を推定する。
フィルタ作成部３６は、整形フィルタＨ（ｑ）の逆フィルタである、白色化フィルタＨ-1（ｑ）を作成し、フィルタリング部３１に出力する。
【0074】
入力ｕ（ｋ）が有色信号であり、白色化フィルタが未知である場合にも、時系列解析部３５およびフィルタ作成部３６の処理で作成した白色化フィルタＨ-1（ｑ）により、フィルタリング部３１がフィルタリング処理を行うことができる。これによって、前記した定理１．が成り立つ状態とすることができ、推定装置１の実効抵抗ＲKの推定精度を向上させることができる。
【符号の説明】
【0075】
１ 推定装置
２０ 過電圧算出部
３０ 推定部
３１ フィルタリング部
３２ 特徴ベクトル算出部（システム同定部）
３３ パラメータ推定部（システム同定部）
３４ 実効抵抗算出部
５０ バッテリ
６０ 電圧センサ
７０ 電流センサ
２０ 過電圧算出部
３０ 推定部
３１ フィルタリング部
３２ 特徴ベクトル算出部（システム同定部）
３３ パラメータ推定部（システム同定部）
３４ 実効抵抗算出部
５０ バッテリ
６０ 電圧センサ
７０ 電流センサ
【要約】
【課題】推定するパラメータの数を削減する。
【解決手段】推定装置１の過電圧算出部２０は、バッテリ５０を流れる電流と、バッテリ５０の端子電圧の測定値に基づいて、バッテリ５０の過電圧を算出する。特徴ベクトル算出部３２およびパラメータ推定部３３は、電流と過電圧の時系列データに基づいて、ＦＩＲモデルを含むモデルを構成する多項式の各項のパラメータを推定することで、バッテリ５０のシステムを同定する。実効抵抗算出部３４は、バッテリ５０のシステムに基づいて、実効抵抗を算出する。特徴ベクトル算出部３２およびパラメータ推定部３３は、ＦＩＲモデルを構成する多項式を複数の時間区分で分割し、各時間区分に含まれる項に対して共通のパラメータを推定することで、ＦＩＲモデルの多項式の項数よりも少ない数のパラメータを推定する。
【選択図】図１
【解決手段】推定装置１の過電圧算出部２０は、バッテリ５０を流れる電流と、バッテリ５０の端子電圧の測定値に基づいて、バッテリ５０の過電圧を算出する。特徴ベクトル算出部３２およびパラメータ推定部３３は、電流と過電圧の時系列データに基づいて、ＦＩＲモデルを含むモデルを構成する多項式の各項のパラメータを推定することで、バッテリ５０のシステムを同定する。実効抵抗算出部３４は、バッテリ５０のシステムに基づいて、実効抵抗を算出する。特徴ベクトル算出部３２およびパラメータ推定部３３は、ＦＩＲモデルを構成する多項式を複数の時間区分で分割し、各時間区分に含まれる項に対して共通のパラメータを推定することで、ＦＩＲモデルの多項式の項数よりも少ない数のパラメータを推定する。
【選択図】図１
【図1】
【図2】
【図3】
【図4】
【図5】
【図6】
【図7】
【図8】
【図9】
【図10】
【図11】
【図12】
【図13】
【図14】
【図15】