特許 7612124 モデルパラメータ推定装置およびモデルパラメータ推定方法

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

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

(11)【特許番号】

(24)【登録日】2024-12-26

(45)【発行日】2025-01-10

(54)【発明の名称】モデルパラメータ推定装置およびモデルパラメータ推定方法

(51)【国際特許分類】

   H01M 10/48 20060101AFI20241227BHJP

   G01R 31/382 20190101ALI20241227BHJP

   G01R 31/385 20190101ALI20241227BHJP

   G01R 31/389 20190101ALI20241227BHJP

   G01R 31/392 20190101ALI20241227BHJP

   G01R 31/367 20190101ALI20241227BHJP

【ＦＩ】

H01M10/48 P

G01R31/382

G01R31/385

G01R31/389

G01R31/392

G01R31/367

【請求項の数】 20

(21)【出願番号】P 2024551013

(86)(22)【出願日】2022-10-14

(86)【国際出願番号】 JP2022038310

(87)【国際公開番号】W WO2024079867

(87)【国際公開日】2024-04-18

【審査請求日】2024-08-29

【早期審査対象出願】

(73)【特許権者】

【識別番号】000006013

【氏名又は名称】三菱電機株式会社

(74)【代理人】

【識別番号】110002941

【氏名又は名称】弁理士法人ぱるも特許事務所

(72)【発明者】

【氏名】竹上 智己

【審査官】辻丸 詔

(56)【参考文献】

【文献】国際公開第２０１８／２２９８９８（ＷＯ，Ａ１）

【文献】特開２０１８－０７７０７６（ＪＰ，Ａ）

【文献】特開２０１５－０８１８００（ＪＰ，Ａ）

【文献】特開２００６－１０５８２３（ＪＰ，Ａ）

【文献】特開２０１４－１８２０７２（ＪＰ，Ａ）

【文献】国際公開第２０１６／０３１１９１（ＷＯ，Ａ１）

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

【文献】特開２０１５－２２４９２７（ＪＰ，Ａ）

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

Ｈ０１Ｍ １０／４８

Ｇ０１Ｒ ３１／３６７

Ｇ０１Ｒ ３１／３８２

Ｇ０１Ｒ ３１／３８５

Ｇ０１Ｒ ３１／３８９

Ｇ０１Ｒ ３１／３９２

Ｈ０２Ｊ ７／００

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

【請求項1】

非線形パラメータと線形パラメータを用いて対象システムを表現した状態空間モデルに対して前記非線形パラメータに値を代入して得られた状態方程式と、前記対象システムの入力と出力それぞれの測定値の時系列データに基づき、前記入力の測定値に対する前記対象システムの状態を示す状態量を算出する状態量算出部、
前記状態空間モデルと前記状態量と前記入力の測定値に基づき算出される出力の推定値と前記出力の測定値との誤差を最小化する前記線形パラメータを推定する線形パラメータ推定部、および
最小化された前記誤差が小さくなるように、所定の収束条件を満たすまで前記非線形パラメータの値の更新を繰り返す非線形パラメータ更新部、
を備えたことを特徴とするモデルパラメータ推定装置。

【請求項2】

前記状態量算出部は、前記対象システムが蓄電池であり、前記入力が前記蓄電池の電流であり、前記出力が前記蓄電池の端子電圧であるとして、前記状態量を算出することを特徴とする請求項１に記載のモデルパラメータ推定装置。

【請求項3】

前記非線形パラメータ更新部は、前記出力の推定値の前記非線形パラメータに対する勾配または勾配近似値を利用して、あるいは勾配フリーの非線形最適化手法を用いて、前記非線形パラメータの値を更新することを特徴とする請求項１または２に記載のモデルパラメータ推定装置。

【請求項4】

前記非線形パラメータ更新部は、前記非線形パラメータの値の更新にあたり、当該更新前の値である前回値を摂動させた前記非線形パラメータそれぞれの摂動値を生成し、
前記状態量算出部は、前記それぞれの摂動値に対する状態量を算出し、
前記線形パラメータ推定部は、前記それぞれの摂動値に対する前記状態量に対する線形パラメータを推定し、
前記非線形パラメータ更新部は、少なくとも前記それぞれの摂動値に対する前記状態量に対する前記線形パラメータの推定値に基づき、前記出力の推定値の前記非線形パラメータに対する勾配近似値を算出し、算出した勾配近似値に基づき前記前回値を更新する
ことを特徴とする請求項１または２に記載のモデルパラメータ推定装置。

【請求項5】

前記非線形パラメータは、前記蓄電池の充電率と満充電容量と拡散時定数と前記電流の測定に用いるセンサのオフセット電流のうちの少なくとも１つを含み、
前記線形パラメータは、前記蓄電池の直流抵抗と拡散抵抗とコンデンサ容量のうちの少なくとも１つを含むことを特徴とする請求項２に記載のモデルパラメータ推定装置。

【請求項6】

前記線形パラメータ推定部は、線形最小二乗問題の線形最小二乗解を用いて前記線形パラメータを推定することを特徴とする請求項２に記載のモデルパラメータ推定装置。

【請求項7】

前記状態空間モデルには、電気量と開回路電圧との関係を表わし、前記線形パラメータを含むＯＣＶ関数を含むことを特徴とする請求項６に記載のモデルパラメータ推定装置。

【請求項8】

前記ＯＣＶ関数は前記電気量の多項式であり、
前記線形パラメータは、前記多項式の係数を含むことを特徴とする請求項７に記載のモデルパラメータ推定装置。

【請求項9】

前記状態空間モデルは、前記線形パラメータに対する線形等式制約を含み、
前記線形パラメータ推定部は、前記線形等式制約を含む線形最小二乗問題をラグランジュ未定乗数法により解くことで前記線形パラメータを推定することを特徴とする請求項２または５に記載のモデルパラメータ推定装置。

【請求項10】

前記状態空間モデルには、電気量と開回路電圧との関係を表わし、前記線形パラメータを含む区分多項式のＯＣＶ関数を含み、
前記線形パラメータは、前記区分多項式の線形パラメータを含み、
前記線形等式制約は、前記区分多項式の節点における線形等式制約を含むことを特徴とする請求項９に記載のモデルパラメータ推定装置。

【請求項11】

前記状態空間モデルは、前記線形パラメータに対する線形不等式制約を含み、
前記線形パラメータ推定部は、前記線形不等式制約を含む線形最小二乗問題を解くことで前記線形パラメータを推定することを特徴とする請求項１または２に記載のモデルパラメータ推定装置。

【請求項12】

前記状態空間モデルは、前記蓄電池の充電率と開回路電圧との関係を表わすＯＣＶ特性情報を含み、
前記非線形パラメータは、満充電容量を含むことを特徴とする請求項２、５、６のいずれか１項に記載のモデルパラメータ推定装置。

【請求項13】

前記状態空間モデルは、前記蓄電池の正極電気量と正極電位との関係を表わす正極電位特性情報と、前記蓄電池の負極電気量と負極電位との関係を表わす負極電位特性情報とを含み、
前記非線形パラメータは、前記満充電容量にかえて、正極容量と負極容量と正極電気量と負極電気量のうちの少なくとも１つを含むことを特徴とする請求項１２に記載のモデルパラメータ推定装置。

【請求項14】

非線形パラメータと線形パラメータを用いて対象システムを表現した状態空間モデルの前記非線形パラメータの初期値を設定するステップ、
前記状態空間モデルに対して前記非線形パラメータに値を代入して得られた状態方程式と、前記対象システムの入力と出力それぞれの測定値の時系列データに基づき、前記入力の測定値に対する前記対象システムの状態を示す状態量を算出する状態量算出ステップ、
前記状態空間モデルと前記状態量と前記入力の測定値に基づき算出される前記出力の推定値と前記出力の測定値との誤差を最小化する前記線形パラメータを推定する線形パラメータ推定ステップ、および
最小化された前記誤差が小さくなるように、所定の収束条件を満たすまで前記非線形パラメータの値の更新を繰り返す非線形パラメータ更新ステップ、
を含むことを特徴とするモデルパラメータ推定方法。

【請求項15】

前記初期値を設定するステップでは、前記対象システムが蓄電池であるとして前記初期値を設定し、前記状態量算出ステップでは、前記入力が前記蓄電池の電流であり、前記出力が前記蓄電池の端子電圧であるとして前記状態量を算出する、ことを特徴とする請求項１４に記載のモデルパラメータ推定方法。

【請求項16】

前記非線形パラメータ更新ステップでは、前記出力の推定値の前記非線形パラメータに対する勾配または勾配近似値を利用して、あるいは勾配フリーの非線形最適化手法を用いて、前記非線形パラメータの値を更新することを特徴とする請求項１４または１５に記載のモデルパラメータ推定方法。

【請求項17】

前記非線形パラメータは、前記蓄電池の充電率と満充電容量と拡散時定数と前記電流の測定に用いるセンサのオフセット電流のうちの少なくとも１つを含み、
前記線形パラメータは、前記蓄電池の直流抵抗と拡散抵抗とコンデンサ容量のうちの少なくとも１つを含むことを特徴とする請求項１５に記載のモデルパラメータ推定方法。

【請求項18】

前記状態空間モデルは、前記線形パラメータに対する線形不等式制約を含み、
前記線形パラメータ推定ステップでは、前記線形不等式制約を含む線形最小二乗問題を解くことで前記線形パラメータを推定することを特徴とする請求項１４または１５に記載のモデルパラメータ推定方法。

【請求項19】

前記状態空間モデルは、前記蓄電池の充電率と開回路電圧との関係を表わすＯＣＶ特性情報を含み、
前記非線形パラメータは、満充電容量を含むことを特徴とする請求項１５に記載のモデルパラメータ推定方法。

【請求項20】

前記状態空間モデルは、前記蓄電池の正極電気量と正極電位との関係を表わす正極電位特性情報と、前記蓄電池の負極電気量と負極電位との関係を表わす負極電位特性情報とを含み、
前記非線形パラメータは、前記満充電容量にかえて、正極容量と負極容量と正極電気量と負極電気量のうちの少なくとも１つを含むことを特徴とする請求項１９に記載のモデルパラメータ推定方法。

【発明の詳細な説明】

【技術分野】

【0001】

本願は、モデルパラメータ推定装置およびモデルパラメータ推定方法に関するものである。

【背景技術】

【0002】

環境負荷低減のため、再生可能エネルギーを活用するための定置用蓄電システムが普及している。また、ＥＶ（電気自動車：Electric Vehicle）、ＨＥＶ（ハイブリッド電気自動車：Hybrid Electric Vehicle）、ＰＨＶ（プラグインハイブリッド自動車：Plug-in Hybrid Vehicle）などの電動車両が実用化され、さらには電動航空機などの開発も進んでいる。

【0003】

一方、これらの機器には、リチウムイオン電池などの蓄電池が用いられているが、蓄電池は使用とともに劣化が進行し、性能が低下することが知られている。それゆえ、蓄電池の高精度なＳＯＣ（充電率：State of Charge）推定、劣化診断、寿命予測などのため、蓄電池を高精度にモデリングするためのモデルパラメータの推定技術が求められている。

【0004】

そこで、出力誤差法に基づく連続時間システム同定法によるモデルパラメータ推定技術が提案されている（例えば、特許文献１参照。）。具体的には、蓄電池充放電時の電流・電圧の時系列データに対し非線形最適化手法を適用することで、連続時間システムとして構築された蓄電池モデルの出力電圧が計測電圧と一致するようにモデルパラメータを推定している。この方法によれば、適切に推定の初期値を選ぶことで、ＯＣＶ（開回路電圧：Open Circuit Voltage）特性の非線形性を考慮しつつ高精度にモデルパラメータを推定することが可能となる。

【先行技術文献】

【特許文献】

【0005】

【文献】特開２０１４－８６３１３号公報（段落００１７～００１８、図１、段落００４４～００４５、図５）

【発明の概要】

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

【0006】

しかしながら、特許文献1で提案された方法では、推定に先立ち全パラメータの初期値を設定する必要があるが、初期値の与え方によっては推定値が発散、あるいは真値とは異なる局所解に収束することがある。それゆえ、対象システムの事前情報に基づく初期値の設定が必須であり、事前情報が少ないシステムへの適用が困難であった。

【0007】

本願は、上記のような課題を解決するための技術を開示するものであり、少ない情報でパラメータを推定することができるモデルパラメータ推定装置およびモデルパラメータ推定方法を得ることを目的とする。

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

【0008】

本願に開示されるモデルパラメータ推定装置は、非線形パラメータと線形パラメータを用いて対象システムを表現した状態空間モデルに対して前記非線形パラメータに値を代入して得られた状態方程式と、前記対象システムの入力と出力それぞれの測定値の時系列データに基づき、前記入力の測定値に対する前記対象システムの状態を示す状態量を算出する状態量算出部、前記状態空間モデルと前記状態量と前記入力の測定値に基づき算出される出力の推定値と前記出力の測定値との誤差を最小化する前記線形パラメータを推定する線形パラメータ推定部、および最小化された前記誤差が小さくなるように、所定の収束条件を満たすまで前記非線形パラメータの値の更新を繰り返す非線形パラメータ更新部、を備えたことを特徴とする。

【0009】

本願に開示されるモデルパラメータ推定方法は、非線形パラメータと線形パラメータを用いて対象システムを表現した状態空間モデルの前記非線形パラメータの初期値を設定するステップ、前記状態空間モデルに対して前記非線形パラメータに値を代入して得られた状態方程式と、前記対象システムの入力と出力それぞれの測定値の時系列データに基づき、前記入力の測定値に対する前記対象システムの状態を示す状態量を算出する状態量算出ステップ、前記状態空間モデルと前記状態量と前記入力の測定値に基づき算出される前記出力の推定値と前記出力の測定値との誤差を最小化する前記線形パラメータを推定する線形パラメータ推定ステップ、および最小化された前記誤差が小さくなるように、所定の収束条件を満たすまで前記非線形パラメータの値の更新を繰り返す非線形パラメータ更新ステップ、を含むことを特徴とする。

【発明の効果】

【0010】

本願に開示されるモデルパラメータ推定装置あるいはモデルパラメータ推定方法によれば、非線形パラメータのみの初期値設定でモデルパラメータを推定することができるので、少ない情報でパラメータを推定することが可能となる。

【図面の簡単な説明】

【0011】

【図1】実施の形態１にかかるモデルパラメータ推定装置の構成を説明するためのブロック図である。

【図2】実施の形態１にかかるモデルパラメータ推定装置の演算処理を実行する部分のハードウェア構成例を示すブロック図である。

【図3】モデルパラメータの推定対象である蓄電池の等価回路を示す図である。

【図4】モデルパラメータの推定対象である蓄電池の別モデルとしてフォスター型等価回路を示す図である。

【図5】実施の形態１にかかるモデルパラメータ推定装置の動作、あるいはモデルパラメータ推定方法を示すフローチャートである。

【図6】実施の形態１にかかるモデルパラメータ推定装置、あるいはモデルパラメータ推定方法をある蓄電池システムに適用した際の電圧の計測値、電流の計測値、およびモデルによる電流の推定値それぞれの経時変化を示すグラフ形式の図である。

【図7】実施の形態１にかかるモデルパラメータ推定装置、あるいはモデルパラメータ推定方法をある蓄電池システムに適用した際に推定されたＯＣＶ関数として電気量と開回路電圧の関係を示すグラフ形式の図である。

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

【0012】

実施の形態１．
図１～図７は、実施の形態１にかかるモデルパラメータ推定装置の構成と動作、およびモデルパラメータ推定方法について説明するためのものであり、図１はモデルパラメータ推定装置の構成を説明するための推定対象である蓄電池とモデルパラメータ推定装置を含めた蓄電池システムのブロック図、図２はモデルパラメータ推定装置の演算処理を実行する部分のハードウェア構成例を示すブロック図である。

【0013】

また、図３はモデルパラメータの推定対象である蓄電池を過電圧の直流成分と緩和成分で表現する等価回路を示す図であり、図４は蓄電池の別モデルとしてＣＲ並列素子を多直列接続して拡散インピーダンスを表現するフォスター型等価回路を示す図である。そして、図５はモデルパラメータ推定装置におけるモデルパラメータ推定動作、つまりモデルパラメータ推定方法を示すフローチャートである。

【0014】

さらに図６はモデルパラメータ推定装置、あるいはモデルパラメータ推定方法をある蓄電池システムに適用した際の、蓄電池の電圧の計測値の経時変化、蓄電池の電流の計測値の経時変化、およびモデルによる電流の推定値の経時変化をそれぞれ示す３つのグラフを縦に並べた図であって、横軸が共通（同期）した時間を示し、縦軸がそれぞれ電圧の計測結果、電流の計測結果、および電流の推計値を示している。また、図７はこのときに推定されたＯＣＶ関数として、ある規格容量によって規格化した電気量を横軸に、開回路電圧を縦軸にしたグラフ形式の図である。

【0015】

以下、本願の実施の形態にかかるモデルパラメータ推定装置、およびモデルパラメータ推定方法について、図面を適宜参照しつつ詳細に説明する。なお、各図において同一符号は同一もしくは相当部分を示している。

【0016】

実施の形態１にかかるモデルパラメータ推定装置３は、図１に示すように、蓄電池１を診断するため、蓄電池１を有する蓄電池システム１００に設けられたものである。そして、蓄電池１の出力特性である電流を検出する電流検出装置２１、電圧を検出する電圧検出装置２２を有する検出部２に接続され、蓄電池１のモデルパラメータを推定するように構成している。

【0017】

蓄電池１は、典型的にはリチウムイオン蓄電池を想定しているが、リチウムイオン蓄電池に限らず、鉛蓄電池、ニッケル水素蓄電池、ニッケルカドミウム蓄電池、全固体蓄電池など他の種類の蓄電池であってもよい。また、診断対象の蓄電池１は、単一セルの蓄電池の他に、複数のセルを任意の直列接続と並列接続の組み合わせで接続した蓄電池モジュールであってもよい。また、モジュールをさらに任意の直列接続と並列接続の組み合わせで接続したパックなどであってもよく、どの単位でひとかたまりの蓄電池１と捉えるかは任意である。

【0018】

モデルパラメータ推定装置３は、検出部２が検出した電流値Iと電圧値Ｖから計測値の時系列データを取得する時系列データ取得部３１と、蓄電池モデルを取得するとともに、後述する非線形パラメータφの初期値を設定するモデル取得部３２を備える。さらに、時系列データ取得部３１が取得した時系列データと、モデル取得部３２が設定した蓄電池モデルと非線形パラメータφの初期値からモデルパラメータを推定するパラメータ推定部３３を備えている。なお、ここでは、モデル取得部３２が初期値の設定を行う例を示したが、別途、あるいはパラメータ推定部３３内に初期値設定部を設けて初期値を設定（あるいは読み込み）するようにしてもよい。

【0019】

パラメータ推定部３３は、非線形パラメータφを更新する非線形パラメータ更新部３３１と、非線形パラメータφを固定した際の蓄電池モデルの状態方程式に基づき、蓄電池モデルの状態量を算出する状態量算出部３３２を備える。さらに、蓄電池モデルと状態量に基づいて線形パラメータψを推定する線形パラメータ推定部３３３と、収束判定条件に基づいてパラメータ推定が収束しているか否かを判定する判定部３３４を備えている。

【0020】

なお、モデルパラメータ推定装置３は、図２に示すように、プロセッサ３０ａと記憶装置３０ｂを備えた一つのハードウェア３０によって構成することが考えられる。記憶装置３０ｂは、図示していないが、ランダムアクセスメモリ等の揮発性記憶装置と、フラッシュメモリ等の不揮発性の補助記憶装置とを具備する。また、フラッシュメモリの代わりにハードディスクの補助記憶装置を具備してもよい。プロセッサ３０ａは、記憶装置３０ｂから入力されたプログラムを実行する。この場合、補助記憶装置から揮発性記憶装置を介してプロセッサ３０ａにプログラムが入力される。また、プロセッサ３０ａは、演算結果等のデータを記憶装置３０ｂの揮発性記憶装置に出力してもよいし、揮発性記憶装置を介して補助記憶装置にデータを保存してもよい。

【0021】

つまり、モデルパラメータ推定装置３を構成する各部の機能、例えば、モデル取得部３２と、状態量算出部３３２と、線形パラメータ推定部３３３と、非線形パラメータ更新部３３１の機能は、ソフトウェア、ファームウェア、またはそれらの組み合わせにより実現される。ソフトウェアおよびファームウェアは、プログラムとして記述されており、記憶装置３０ｂに格納されている。プロセッサ３０ａは、記憶装置３０ｂに記憶されたプログラムを読み出して、そのプグラムを実行することにより、モデルパラメータ推定装置３の各部の機能を実現する。

【0022】

＜本願の技術の一般的な説明＞
本願のモデルパラメータ推定装置３およびモデルパラメータ推定方法の技術内容について、はじめに、対象を蓄電池に限定しない一般的な問題設定で説明する。
まず、対象システムの状態空間モデルは式（１）のように表されるものとする。

【0023】

【数1】

ただし、ドットｘは状態量を表わすベクトル、ｕは入力、ｙは出力であり、ｆとｇはそれぞれある非線形なベクトル値関数である。また、φ_ｘとφ_ｙは、それぞれ状態方程式と出力方程式における非線形パラメータベクトルであり、ψは線形パラメータベクトルである。

【0024】

ここで解きたいのは、等間隔または不等間隔で任意にサンプリングされた時刻｛ｔ_ｋ｜ｋ＝０，・・・，Ｎ｝での時系列入出力データ｛ｕ（ｔ_ｋ）｜ｋ＝０，・・・，Ｎ｝,｛ｙ（ｔ_ｋ）｜ｋ＝０，・・・，Ｎ｝の計測値が与えられたとき、式（２）に示すモデルの出力ハットｙ（推定値）と計測値ｙとの誤差εの二乗和に基づく評価関数Ｊ（θ）を最小化するようなパラメータθを求めることである。

【0025】

【数2】

ただし、θは線形パラメータψと非線形パラメータφを並べたベクトルであり、ψ：＝［ψ_ｘ ^Τ，ψ_ｙ ^Τ］^Τに対し、θ：＝［φ^Τ，ψ^Τ］^Τと定義されるようなものである。

【0026】

これは非線形最適化問題であることから、公知の非線形最適化手法により局所最適解を得ることが可能である。例えば、評価関数のパラメータに関する勾配、または推定出力 ハットｙ（ｔ_ｋ，θ）のパラメータに関する勾配、に基づくガウス＝ニュートン法（Gauss-Newton method）、あるいはレーベンバーグ＝マーカート法（Levenberg-Marquardt method）などを用いることが可能である。また、勾配フリー、つまり勾配情報に基づかないネルダー＝ミード法（Nelder-Mead method）などの手法を用いることも可能であり、非線形最適化の具体的手法に関して限定されることはない。

【0027】

パラメータに関する勾配を求めるに当たっては、特許文献１に記載されているような感度方程式の計算に基づく方法を用いることができる。また、片側差分または両側差分などによる勾配近似手法を用いてもよく、随伴変数を導入することによる随伴変数法（アジョイント法：Adjoint method）を用いてもよい。

【0028】

なお、式（１）では連続時間の状態空間モデルを記載しているが、対象システムを離散時間で考えてもよい。このとき、状態空間モデルは式（３）のように記述される。

【0029】

【数3】

ただし、ｘ_ｋは時刻ｋでの状態量を表わすベクトル、ｕ_ｋとｙ_ｋはそれぞれ時刻ｋでの入力と出力であり、ｆ_ｄとｇ_ｄはそれぞれある非線形なベクトル値関数である。連続時間システムを離散時間システムに変換するには、公知の方法、例えばゼロ次ホールドあるいは一次ホールドなどを用いることが可能である。

【0030】

一方、連続時間システムをそのまま扱う場合には、もとの状態方程式に対し公知の微分方程式の数値解法を適用することで、入力された時系列データから任意の時刻の状態量を計算することが可能となる。微分方程式の数値解法として、典型的にはオイラー法（Euler method）またはルンゲクッタ法（Runge-Kutta method）を用いればよい。

【0031】

しかしながら、非線形最適化問題を扱う出力誤差法一般の課題として、パラメータ推定値が初期値依存性をもつことが挙げられる。それゆえ、背景技術で述べたように、推定パラメータの初期値を適切に設定しないと、大域的最適解が得られないという問題がある。そこで、式（１）の状態空間モデルの構造に着目し、初期値依存性を減らすことを考える。まず、モデルの出力を時系列に従い縦に並べると、式（４）のようになる。

【0032】

【数4】

つまり、線形パラメータψと非線形パラメータφを分離した形で表現できる。このとき、出力ベクトルＹ：＝［ｙ（ｔ_０），・・・，ｙ（ｔ_Ｎ）］^Τと誤差ベクトルＥ：＝［ε（ｔ_０），・・・，ε（ｔ_Ｎ）］^Τに対し、式（５）関係式が成り立つ。
Ｙ＝Ｇ（φ)ψ＋Ｅ ・・・（５）

【0033】

そのため、式（５）における非線形パラメータφを固定したとき、線形パラメータψの線形最小二乗解は、線形最小二乗法により式（６）のように求まる。
ψ(φ)＝(Ｇ^Τ(φ)Ｇ(φ))^－１Ｇ^Τ(φ)Ｙ ・・・（６）
式（６）の線形最小二乗解を利用すると、解くべき問題は、元の最適化問題が式（７）の問題Ｐ１であったのに対し、式（８）の問題Ｐ２に変形される。

【0034】

【数5】

【0035】

つまり、元々線形パラメータψと非線形パラメータφを全て推定する必要があった問題Ｐ１が、非線形パラメータφのみの推定問題（問題Ｐ２）に変形された。問題Ｐ２の非線形パラメータφの解が得られれば（固定できれば）、これを式（６）に代入することで、線形パラメータψの解も即座に得られる。

【0036】

なお、問題Ｐ２もまた非線形最適化問題であることから、問題Ｐ１を解く場合と同様に、上述のような公知の様々な非線形最適化手法を用いることが可能である。問題Ｐ１を問題Ｐ２に変換すると、推定パラメータが非線形パラメータφのみになることで、以下の利点が生じる。

【0037】

第一に、初期値依存性を低下させることができる。元々の問題Ｐ１を解く場合には線形パラメータψと非線形パラメータφ全ての初期値を設定する必要があったのに対し、問題Ｐ２においては、線形パラメータψは線形最小二乗解として初期値なしに求まる。そのため、非線形パラメータφの初期値のみを設定すれば問題を解くことが可能となっている。したがって、パラメータの真値を推定するうえで問題がより簡単なものとなっている。

【0038】

第二に、多くの非線形最適化手法の適用において、計算量が削減される。例えば、モデル出力のパラメータに関する勾配情報を利用する手法を利用する場合には、線形パラメータψの勾配情報が不要となる分、勾配ベクトル、あるいはヤコビ行列（Jacobian matrix）などのサイズが小さくなる。また、別の観点として、推定パラメータ数が減ることにより、パラメータが収束するまでの繰り返し計算の回数が少なくて済むことが期待される。

【0039】

第三に、数値的安定性が向上する。一般に推定パラメータ間で値の大きさのオーダーが大きく異なる場合には桁落ちなどの問題が生じやすい。ゆえに、推定パラメータ数が減ることでこうした問題が生じにくくなると期待される。

【0040】

つぎに、線形パラメータψが、線形等式制約を含む場合と、線形等式制約と線形不等式制約を含む場合に分けて、解を得るための具体的な手法について説明する。

【0041】

＜１：線形パラメータが線形等式制約を含む場合＞
線形パラメータψが線形等式制約を含む場合、線形パラメータψの等式制約付き最小二乗問題は式（９）の問題Ｐ３のように定式化される。

【0042】

【数6】

【0043】

ここで、Ｃ_ｅｑとｄ_ｅｑは、それぞれ線形パラメータψに対する線形等式制約を記述する行列とベクトルである。問題Ｐ３を先程と同様の思想で解くため、非線形パラメータφをある値で固定したＧ＝Ｇ（φ）に対し、問題Ｐ３の副問題ＳＰ３を考えると式（１０）のようになる。

【0044】

【数7】

副問題ＳＰ３は等式制約付き線形最小二乗問題であるから、ラグランジュ未定乗数法（method of Lagrange multiplier）で解くことができる。具体的には、副問題ＳＰ３はラグランジュ乗数ベクトルλを用いた式（１１）に示す評価関数Ｌを最小化するようなψを求めればよい。

【0045】

【数8】

そこで、評価関数Ｌのψとλそれぞれによる勾配が０になる点を考えると式（１２）となる。

【0046】

【数9】

まとめると式（１３）となり、解を即座に求めることが可能である。

【0047】

【数10】

それゆえ、問題Ｐ３を解くためには、非線形最適化手法による非線形パラメータの更新と、ラグランジュ未定乗数法による線形パラメータの算出を交互に繰り返せばよい。

【0048】

＜２：線形パラメータが線形等式制約と線形不等式制約を含む場合＞
線形パラメータが線形等式制約と線形不等式制約を含む場合、線形パラメータの不等式制約付き最小二乗問題は式（１４）の問題Ｐ４のように定式化される。

【0049】

【数11】

ここで、Ｃ_ｉｎｅｑとｄ_ｉｎｅｑはそれぞれ、線形パラメータψに対する線形不等式制約を記述する行列とベクトルである。問題Ｐ３のときと同様にして問題Ｐ４の副問題ＳＰ４を考えると式（１５）のようになる。

【0050】

【数12】

【0051】

副問題ＳＰ４を解く方法として、内点法、逐次２次計画法などを用いることが可能である。この場合、繰り返し計算が必要となるが、副問題ＳＰ４は凸最適化問題であることから、大域的最適解が得られる。もちろん、副問題ＳＰ３は副問題ＳＰ４に含まれるため、副問題ＳＰ３に対してもこのような繰り返し計算による解法を適用することが可能である。それゆえ、問題Ｐ４を解くためには、非線形最適化手法による非線形パラメータの更新と、凸最適化問題の解法による線形パラメータの算出を繰り返せばよい。

【0052】

ここから、問題Ｐ１、問題Ｐ３、問題Ｐ４に対する非線形最適化手法の具体的な適用方法について、問題Ｐ４を対象に説明する。なお、問題Ｐ１および問題Ｐ３は問題Ｐ４に含まれるため、問題Ｐ４を対象に説明すれば十分である。実際、問題Ｐ１は問題Ｐ４において不等式制約と等式制約が存在しないという特殊な場合である。また、問題Ｐ３も、問題Ｐ４において不等式制約が存在しないという特殊な場合である。

【0053】

副問題ＳＰ４におけるある非線形パラメータφに対する線形パラメータψの最適解をψ＝ψ^＊（φ）としたとき、問題Ｐ４で最小化したい評価関数は、式（１６）となる。

【0054】

【数13】

【0055】

ψ^＊（φ）は上述の方法で副問題ＳＰ４を解けば得られることから、問題Ｐ４を解くためにはＬ（φ）を最小化する非線形パラメータφを見つければよい。以下では、非線形最適化手法として、勾配法の一種である最急降下法、ガウス＝ニュートン法、および勾配フリーの方法の一つであるネルダー＝ミード法について順に説明する。

【0056】

＜最急降下法＞
最急降下法を用いる場合は、非線形パラメータφの更新式は式（１７）となる。

【0057】

【数14】

ただし、αは所定の正の値であり、更新のたびに直線探索などにより決定してもよい。

【0058】

勾配を直接求めるのが困難な場合は、勾配近似値として片側差分、あるいは両側差分などの値を用いることができる。たとえば、片側差分においては、φ_ｊを十分に小さなΔ_ｊだけ摂動させることで偏微分を式（１８）のように近似的に計算することができる。

【0059】

【数15】

ただし、ｅ_ｊは非線形パラメータφと同じ長さでｊ番目の要素のみ１であるような単位ベクトルである。つまり、φ＋Δ_ｊｅ_ｊは非線形パラメータφの第ｊ要素のみを摂動させた摂動値である。このとき、非線形パラメータφの更新式は式（１９）となる。

【0060】

【数16】

なお、評価関数Ｌの計算においては、非線形パラメータφの値ごとに状態量を算出し、線形パラメータψを推定し、出力推定値を算出することが必要となる。

【0061】

＜ガウス＝ニュートン法＞
ガウス＝ニュートン法の更新式は、式（２０）となる。ただし、Λ（ψ）は式（２１）で表され、∂ハットｙ（ｔ_ｋ，φ，ψ^＊（φ））／∂φ_ｊが感度関数、すなわち出力推定値のある非線形パラメータφ_ｊに関する偏微分である。

【0062】

【数17】

【0063】

感度関数を直接求めることが困難な場合には、感度関数に代えて感度関数の近似値を用いる。つまり、勾配そのものではなく勾配近似値を用いる。この場合にも、例えば片側差分を用いる場合の計算式は式（２２）となる。

【0064】

【数18】

【0065】

＜ネルダー＝ミード法＞
一例として、ネルダーミード法を用いる場合、１回のパラメータ更新の内部で、繰り返し回数ｌ（エル）、非線形パラメータ数ｎ_φに対してｎ_φ＋１個の変数φ^ｌ，１，φ^ｌ，２，・・・，φ^{ｌ，ｎφ＋１}を用い、評価関数値Ｌ（φ^ｌ，１），Ｌ（φ^ｌ，２），・・・，Ｌ（φ^{ｌ，ｎφ＋１}）を計算する。さらに、鏡像、拡大、収縮、縮小という４種類の操作を行なったパラメータ値に対しても評価関数値を計算する。そして、このように算出した評価関数値を利用しつつ、公知のアルゴリズムに従いパラメータをφ^ｌからφ^ｌ＋１に更新する。

【0066】

このように、非線形最適化の方法によっては、非線形パラメータ推定値を更新するために、異なる複数の非線形パラメータ推定値に対する評価関数または出力推定値を求めることが必要となる。

【0067】

なお、本願の主旨から外れるため説明を省略していたが、元々の問題Ｐ１が非線形パラメータφに対して等式制約と不等式制約の少なくとも一方を含む場合においても、非線形パラメータφの更新において公知の制約条件付き非線形最適化手法の適用が可能である。例えば上述の内点法、逐次２次計画法などを適用することで局所解を得ることができる。

【0068】

以上が一般問題を対象として説明した本願のパラメータ推定技術の内容である。この内容を前提とし、本技術を蓄電池１に適用した場合について説明する。

【0069】

＜蓄電池への適用＞
図１で説明したように、電流検出装置２１は、蓄電池１の電流を検出し、検出した電流値Iを出力する。電圧検出装置２２は、蓄電池１の端子電圧を検出し、検出した電圧値Ｖを出力する。ただし、蓄電池１が複数セルの直並列接続で構成されている場合、検出電流および検出電圧は、蓄電池１を構成する一部である単一または複数のセルのものであってもよい。そして、以下では、時系列データのサンプリング周期をｔ_ｓ秒であるとして説明する。ただし、サンプリング周期ｔ_ｓは固定されている必要はなく、可変であってもよい。

【0070】

時系列データ取得部３１は、ある時刻の時系列｛ｔ_ｋ｜ｋ＝０,１，・・・，Ｎ｝に対する検出電流の時系列データ｛Ｉ（ｔ_ｋ）｜ｋ＝０,１，・・・，Ｎ｝および検出電圧の時系列データ｛Ｖ（ｔ_ｋ）｜ｋ＝０，１，・・・，Ｎ｝を取得する。簡単のためサンプリング周期ｔ_ｓが固定の場合を考えると、初期時刻ｔ_０に対しｔ_ｋ＝ｔ_０＋ｋｔ_ｓの関係が成り立つ。ただし、時系列データは、時系列データ取得部３１が内部に保持していてもよいし、外部のＰＣ、サーバ、クラウドなどを経由して取得してもよい。

【0071】

モデル取得部３２は、推定対象となる蓄電池１に対する蓄電池モデルを保有しており、パラメータ推定部３３において利用される。そして、初期値を含む蓄電池モデルのパラメータは、パラメータ推定部３３で推定されたパラメータに基づき適宜更新される。なお、モデル取得部３２が保有する蓄電池モデルは、モデルパラメータ推定装置３が元々内部に保有していてもよいし、外部のＰＣ、サーバ、クラウドなどから取得したものであってもよい。

【0072】

蓄電池モデルとして、例えば図３に示すような等価回路モデルを用いる。図３で抵抗素子Ｒ_０は、蓄電池１の電解液抵抗、集電体金属抵抗、接触抵抗、電荷移動抵抗などに由来する過電圧の直流成分を表現する。また、抵抗素子Ｒ_ｄとコンデンサ素子Ｃ_ｄは、蓄電池１の拡散成分に由来する過電圧の緩和成分を表現する。

【0073】

このとき、蓄電池１の状態空間モデルは、例えば式（２３）～式（２５）のように記述される。

【0074】

【数19】

【0075】

ここで、ｑは電気量、Ｉは検出した電流値（検出電流）、Ｉ_ｏｆｆは検出電流に含まれる電流オフセット誤差、ｑ_ｄはコンデンサ素子Ｃ_ｄに蓄えられる電荷、ｆ_ＯＣＶは蓄電池１の電気量ｑに依存するＯＣＶ特性を表わす関数（ＯＣＶ関数ｆ_ＯＣＶ）である。ただし、電気量ｑは規格化された値であってもよく、その場合はＯＣＶ関数ｆ_ＯＣＶもまた規格化された電気量ｑに対する関数とする。

【0076】

離散時間の状態空間モデルを用いる場合には、例えば、ゼロ次ホールドにより、式（２３）～式（２５）を式（２６）～式（２８）のように変換したものを用いる。

【0077】

【数20】

【0078】

別の蓄電池モデルとして、図４に示すようなＣＲ並列素子を多直列接続したフォスター型の等価回路モデルを用いてもよい。図４の等価回路モデルは、図３の等価回路モデルと比較して、ＣＲ並列素子が３直列となっているため、より詳細な蓄電池モデルとなっている。もちろん、ＣＲ並列素子の直列接続数は任意であってよい。

【0079】

また、拡散インピーダンスをフォスター（Foster）型回路で近似した場合、各ＣＲ並列素子は拡散抵抗Ｒ_ｄと拡散容量Ｃ_ｄのみを用いて以下のように表わすことができる（例えば、Kuhn, Estelle, et al. "Modelling Ni-mH battery using Cauer and Foster structures." Journal of power sources 158.2 (2006): 1490-1497.参照。）。

【0080】

【数21】

ここで、ｎ_ｄは近似次数である。例えば近似次数ｎ_ｄ＝３のとき、図４の等価回路モデルのように３段のＣＲ並列素子で拡散インピーダンスを近似したことになる。

【0081】

このとき、式（２７）と式（２８）は、式（３０）と式（３１）のようになる。

【0082】

【数22】

ただし、ｋ_ｉは式（３２）に記載の通りである。
ｋ_ｉ＝４／（（２ｉ－１）^２π^２） ・・・（３２）

【0083】

なお、蓄電池モデルとして他のものを用いてもよく、例えば、カウエル（Cauer）型回路を用いることも可能である。また、等価回路を用いずに蓄電池１の内部現象を直接記述することを試みたいわゆる電気化学モデルを用いてもよい。蓄電池１のモデルとしては公知の様々なものを利用することが可能であり、式（１）または式（３）の枠組みで表現可能であれば限定されることはない。ただし、以下では簡単のため、図３で説明した等価回路モデルを仮定して説明する。

【0084】

ＯＣＶ特性を表わすＯＣＶ関数ｆ_ＯＣＶには、典型的には、電気量ｑに対し線形または区分線形な関数を用いる。例えば、多項式を用いる場合は式（３３）のように表わされ、パラメータ線形な表現が可能である。

【0085】

【数23】

ただし、ｍ_ｐは多項式の次数を表わし、ｃ_ｉはｑ^ｉの係数である。

【0086】

あるいは、式（３４）に示す区分線形関数（１次スプライン曲線）を用いる。

【0087】

【数24】

ただし、ｍ_ｓ（≧２）は節点の個数であり、式（３４）においては関数がｑ＝ｐ_１，ｐ_２，・・・，ｐ_ｍｓのｍ_ｓ点によりｍ_ｓ―１個の区間に分割されており、ｍ_ｓ≧３のとき等式制約として式（３５）を満たす。
a_ｉp_ｉ＋１＋b_ｉ＝a_ｉ＋１p_ｉ＋１＋b_ｉ＋１ for i＝1,…,ｍ_ｓ－2 ・・（３５）

【0088】

あるいは、式（３６）の３次スプライン曲線を用いる。

【0089】

【数25】

ただし、ｍ_ｓ≧３のとき等式制約として式（３７）を満たす。

【0090】

【数26】

【0091】

あるいは、他の任意の次数のスプライン曲線を用いる。

【0092】

スプライン曲線のような区分多項式においても、パラメータ線形な表現が可能である。一例として、区分線形関数の場合、式（３８）のように表現することが可能である。ただし、δ_ｉ，ｋは式（３９）に示す通りである。

【0093】

【数27】

【0094】

同様にして、３次スプライン曲線など他の場合においてもパラメータ線形な表現が可能である。また、他のＯＣＶ関数ｆ_ＯＣＶとして、線形パラメータψと非線形パラメータφを両方含む場合を考えることも可能である。その場合は、関数形が式（４０）のようになる。

【0095】

【数28】

ただし、ｇ_ＯＣＶはあるベクトル値関数であり、φ_ＯＣＶとψ_ＯＣＶはそれぞれ、ＯＣＶ関数ｆ_ＯＣＶに含まれる非線形パラメータφと線形パラメータψのベクトルである。

【0096】

他のＯＣＶ関数ｆ_ＯＣＶとして、非線形パラメータφのみを含む場合を考えることも可能である。つまり、ＯＣＶ関数ｆ_ＯＣＶの関数形は必ずしも線形パラメータψを含む形に制限されることはない。ただし、線形パラメータψを数多く含む場合において、とくに、本願における問題Ｐ１を問題Ｐ２に変換するアプローチの利点が強く生かされることになる。ただし、以下では式（３４）をＯＣＶ関数ｆ_ＯＣＶとして用いた場合を説明する。

【0097】

非線形パラメータ更新部３３１は、検出電流の時系列データ｛Ｉ（ｔ_ｋ）｜ｋ＝０，１，・・・，Ｎ｝と検出電圧の時系列データ｛Ｖ（ｔ_ｋ）｜ｋ＝０，１，・・・，Ｎ｝と蓄電池モデルに基づき、蓄電池モデルの非線形パラメータφを更新する。ただし、ここでの非線形パラメータφは式（４１）のようになる。

【0098】

【数29】

つまり、蓄電池１の状態空間モデルにおいて、状態方程式に含まれるパラメータと、出力方程式に含まれる非線形パラメータ（パラメータ線形表現できないパラメータ）を指す。ただし、特定の非線形パラメータ値が既知である場合など、必ずしも全ての非線形パラメータφを推定パラメータに含める必要はない。

【0099】

なお、蓄電池１の状態空間モデルの作り方に応じて、非線形パラメータφの要素は様々に異なりうる。例えば、ＯＣＶ関数ｆ_ＯＣＶに非線形パラメータφ_ＯＣＶが含まれない場合、必然的に、非線形パラメータφにもφ_ＯＣＶは含まれない。また、式（２４）、式（３０）の状態方程式に含まれるパラメータは時定数τ_ｄと電流オフセット誤差Ｉ_ｏｆｆのみであるが、モデル化の方法に応じて、他のパラメータも状態方程式に含まれうる。

【0100】

例えば、多段ＣＲ並列素子を式（２９）のように拡散インピーダンスの近似により表現した場合の時定数はτ_ｄのみの１パラメータとなるが、そうでない場合はＣＲ並列素子の数だけ時定数のパラメータが独立して存在することになる。また、例えば式（２４）のＣＲ過電圧の状態方程式に代えて式（４２）のＣＲ過電圧の状態方程式を用い、式（２５）に代えて出力方程式を式（４３）とした場合、線形パラメータψではなく非線形パラメータφにＣ_ｄを含める。

【0101】

【数30】

Ｖ＝Ｒ_０（Ｉ－Ｉ_ｏｆｆ）＋Ｖ_ｄ＋ｆ_ＯＣＶ（ｑ） ・・（４３）

【0102】

ただし、このようにした場合、τ_ｄとＣ_ｄが両方とも非線形パラメータφに含まれてしまうことを考えると、本願の技術思想の長所を活かすには式（２４）を用いたほうがよい。

【0103】

非線形パラメータφの具体的な更新方法は、上述したように公知の繰り返し計算に基づく非線形最適化手法を利用することが可能である。

【0104】

非線形パラメータφの初期値には、所定の値を用いる。所定の値は、可能であれば非線形パラメータφの真値に関する事前知識を用いる。例えば、拡散時定数τ_ｄはリチウムイオン蓄電池においては通常数十秒～数百秒程度であることが知られているため、初期値として例えば１００秒などと設定する。また、電流オフセット誤差Ｉ_ｏｆｆは微小な値であるため、初期値として例えば０Ａなどと設定する。

【0105】

これらの設定は一般的なリチウムイオン蓄電池および電流センサに対し妥当性の高い値であり、蓄電池１の直並列接続構成などにも依存しないことから、多くの場合において初期値チューニング作業が不要となるという意味で製品搭載上の利点を有する。また、蓄電池１の新品時の特性評価結果、またはスペックシート上のデータなどが利用できる場合は、これを利用してもよい。また、過去にパラメータ推定を実施している場合には、過去の推定結果を利用してもよい。

【0106】

なお初期値は、モデル取得部３２において設定してもよいが、非線形パラメータ更新部３３１において設定してもよく、別途設けた図示しない初期値設定部で設定するようにしてもよい。その際、上述した事前知識に基づき、自動設定するようにしてもよいが、図示しない入力Ｉ／Ｆで事前情報、初期値候補等を提示したうえで入力を促し、入力に基づいて初期値を設定するようにしてもよい。

【0107】

なお、非線形パラメータφの更新には、後述する状態量算出部３３２と線形パラメータ推定部３３３の動作を内部に含んでいてもよい。例えば、非線形パラメータφの更新に評価関数の非線形パラメータφに関する勾配情報を用いる場合には、状態量算出部３３２と線形パラメータ推定部３３３と同様の計算が必要となる。

【0108】

状態量算出部３３２は、検出電流の時系列データ｛Ｉ（ｔ_ｋ）｜ｋ＝０，１，・・・，Ｎ｝と蓄電池モデルとに基づき、蓄電池モデル（すなわち蓄電池１の状態空間モデル）における状態量の時系列値を算出する。例えば、式（２３）～式（２５）の状態空間モデルにおける状態量とは、電気量ｑと電荷ｑ_ｄを指す。

【0109】

状態空間モデルに基づく状態量の時系列値の計算には、上述したように連続時間と離散時間のいずれの状態空間モデルを用いてもよい。前者の場合には、例えば４次のルンゲクッタ法などの公知の数値解法を利用すればよく、後者の場合には、モデルに従って逐次計算していけばよい。あるいは、状態方程式の解に基づいて計算してもよい。

【0110】

線形パラメータ推定部３３３は、検出電流の時系列データ｛Ｉ（ｔ_ｋ）｜ｋ＝０，１，・・・，Ｎ｝と検出電圧の時系列データ｛Ｖ（ｔ_ｋ）｜ｋ＝０，１，・・・，Ｎ｝と蓄電池モデルと状態量とに基づき、線形パラメータψを推定し、電圧のモデル出力を算出する。具体的には、式（４４）のように出力方程式において線形パラメータψを分離して表現できるため、式（４５）のように表わすことができる。

【0111】

【数31】

【0112】

ゆえに、ある非線形パラメータφの推定値に対する線形パラメータψの線形最小二乗解は、式（４６）と推定できる。
ψ（φ）＝（Ｇ^Τ（φ）Ｇ（φ））^－１Ｇ^Τ（φ）ν ・・（４６）
よって、蓄電池１における出力誤差法に基づく電圧誤差最小化問題である問題Ｐ１は、式（４６）により問題Ｐ２に変換できる。また、この線形パラメータψの推定値を用いれば、電圧のモデル出力も算出できる。

【0113】

等式制約を含む場合は、等式制約を考慮して問題Ｐ３を考えればよい。例えば、ＯＣＶ関数ｆ_ＯＣＶが１次スプライン曲線（区分線形関数）の場合は、式（３５）を満たさなければならないから、ψ_{ＯＣＶ，ｉ}＝［ａ_ｉ ｂ_ｉ］^Τとしたとき、式（４７）に対し、等式制約式は式（４８）となる。
χ_ｉ＝［ｐ_ｉ＋１ １］^Τ ・・（４７）

【0114】

【数32】

【0115】

また、ＯＣＶ関数ｆ_ＯＣＶが３次スプライン曲線の場合は、式（３７）を満たさなければならないから、ψ_{ＯＣＶ，ｉ}＝［ａ_ｉ ｂ_ｉ ｃ_ｉ ｄ_ｉ］^Τとしたとき、式（４９）に対し、等式制約式は式（５０）となる。

【0116】

【数33】

【0117】

このように定式化することで、解きたい問題を問題Ｐ３および副問題ＳＰ３として扱うことが可能となる。このときの解法は既に説明したとおりである。なお、制約等式は上述の例で説明したものに限られず、線形制約等式で表現される単一または複数の任意の等式であってよい。

【0118】

また、線形パラメータψに対する不等式制約を含む場合を考えることも可能である。例えば、Ｒ_０とＣ_ｄに関しＲ_{０，ｍｉｎ}≦Ｒ_０≦Ｒ_{０，ｍａｘ}，Ｃ_{ｄ，ｍｉｎ}≦Ｃ_ｄ≦Ｃ_{ｄ，ｍａｘ}という上下限制約に関する事前知識を反映させたい場合は、不等式制約を式（５１）のように表現すればよい。

【0119】

【数34】

【0120】

これにより、解きたい問題を問題Ｐ４および副問題ＳＰ４として扱う事が可能となる。このときの解法は既に説明したとおりである。なお、制約不等式は上述の例で説明したものに限られず、線形制約不等式で表現される単一または複数の任意の不等式であってよい。

【0121】

なお、ここではＯＣＶ関数ｆ_ＯＣＶの多項式近似方法およびスプライン近似方法を説明したが、同様の方法でＲ_０，Ｒ_ｄ，Ｃ_ｄ，τ_ｄなどの回路パラメータの電気量依存性をモデル化してもよい。また、同様の方法で回路パラメータの温度依存性をモデル化してもよい。ただし、温度依存性をモデル化する場合には蓄電池１の温度を別途計測する必要があり、その場合にはモデルパラメータ推定装置３が温度検出器の出力を取得する必要がある。なお、温度モデルとしては、公知のアレニウス式（Arrhenius equation）などが利用可能である。

【0122】

判定部３３４は、検出電圧の時系列データ｛Ｖ（ｔ_ｋ）｜ｋ＝０，１，・・・，Ｎ｝と、電圧のモデル出力と、線形パラメータ推定値と、非線形パラメータ推定値とに基づいて、パラメータ推定値および電圧推定値の収束を判定する。所定の判定基準を満たせば推定を終了し、モデルパラメータを外部に出力する。判定基準を満たさなければ、非線形パラメータ更新部３３１に戻り、非線形パラメータφの更新を繰り返す。

【0123】

具体的には、パラメータθのｌ（エル）回目の推定値の第ｊ要素をθ_ｊ ^ｌとしたとき、特許文献１と同様に式（５２）を満たしているかどうかを判定する。

【0124】

【数35】

ただし、ｎ_θはθの要素数、ε_θは所定の十分に小さな値である。

【0125】

あるいは、非線形パラメータφのみに対して式（５３）を満たしているかどうかを判定する。

【0126】

【数36】

ただし、ｎ_φはφの要素数、ε_φは所定の十分に小さな値である。

【0127】

あるいは、電圧に関する二乗平均平方根誤差（ＲＭＳＥ:Root-Mean-Square Error）を基準にして、式（５４）を満たしているかどうかを判定する。

【0128】

【数37】

ただし、ε_Ｖは所定の十分に小さな値である。

【0129】

あるいは、電圧に関する平均絶対値誤差（ＭＡＥ:Mean Absolute Error）を基準にして判定しても良い。あるいは、非線形パラメータφの推定の繰り返し回数が所定の上限値に達しているかどうかを収束の条件として判定する。あるいは、上記の複数の方法による複合的な判定を行なう。あるいは、他の任意の収束判定方法を用いてもよい。非線形最適化における推定パラメータの収束判定方法として、公知の様々なものが知られているため、収束判定方法に関して限定されることはない。

【0130】

なお、パラメータ推定部３３の最終的なパラメータ推定値は、モデル取得部３２に保存して次回のモデルパラメータ推定の初期値に利用してもよいし、外部のＰＣ、サーバ、クラウド上などに保存してもよい。

【0131】

＜モデルパラメータ推定装置の処理手順＞
つぎに、実施の形態１にかかるモデルパラメータ推定装置３の処理手順、つまりモデルパラメータ推定方法について、図５のフローチャートを参照して説明する。ただしここでは、最も典型的な場合、すなわち勾配近似値を用いた非線形最適化手法を用いた場合を説明する。

【0132】

はじめに、時系列データ取得部３１が、電流検出装置２１が出力する蓄電池１の電流値Iと、電圧検出装置２２が出力する蓄電池１の電圧値Ｖと、をそれぞれ計測値の時系列データとして取得する（ステップＳ１３１０）。つぎに、モデル取得部３２が、外部から、または内部に保有する蓄電池モデルを取得し、非線形パラメータφの初期値の設定または読み込みを行う（ステップＳ１３２０）。蓄電池モデルとしては、式（３）に対応する例えば、式（２３）～式（２５）、あるいは式（２６）～式（２８）に示した状態空間モデルを取得する。

【0133】

取得した時系列データ、蓄電池モデル、および初期値はパラメータ推定部３３に出力され、パラメータ推定部３３内の各部で以下のような繰り返し演算（ステップＳ２３３１～ステップＳ２４００）が実行される。

【0134】

ステップＳ２３３１では、非線形パラメータ更新部３３１が、繰り返し計算に基づく勾配法の一種であるような非線形最適化手法に従い、線形パラメータψに対して最小化された誤差が小さくなるように、非線形パラメータφを更新する。勾配の計算において、例えば片側差分による勾配近似値を用いる場合、後述する前回の非線形パラメータ更新値（前回値）の要素毎の摂動値に対する出力推定値に基づき、勾配近似値を計算する。ただし、片側差分などによる勾配近似値の計算においては、前回値に対する出力推定値も利用する。なお、初回は、非線形パラメータの初期値として所定の値をそのまま設定する（例えばτ^d＝１００、Ｉ_ｏｆｆ＝０など）。

【0135】

ステップＳ２３３２では、状態量算出部３３２が、更新された非線形パラメータφによる蓄電池モデルの状態方程式に基づき、蓄電池モデルの状態量を算出する。ステップＳ２３３３では、線形パラメータ推定部３３３が、状態量算出部３３２から出力された状態量に基づき、線形パラメータψおよび蓄電池モデルの出力値を推定する。

【0136】

そして、ステップＳ２３３４では、判定部３３４が、所定の収束判定条件により、パラメータ推定部３３におけるパラメータ推定の収束を判定する。収束判定条件を満たすと判定した場合（ステップＳ２４００で「Ｙｅｓ」）は、パラメータ推定を終了し、推定結果のパラメータを外部に出力する。このとき、推定結果のパラメータはモデル取得部３２で保持してもよい。

【0137】

一方、収束判定条件を満たさないと判定した場合（ステップＳ２４００で「Ｎｏ」）は、ステップＳ３３３１に進み、ステップＳ３３３１以下の各ステップの実行を継続する。

【0138】

ステップＳ３３３１では、非線形パラメータ更新部３３１が、非線形パラメータ要素毎の摂動値を設定する。ステップＳ３３３２では、状態量算出部３３２が、非線形パラメータ更新部３３１で設定した非線形パラメータ要素毎の摂動値それぞれによる蓄電池モデルの状態方程式に基づき、非線形パラメータ要素毎の摂動値それぞれに対する蓄電池モデルの状態量を算出する。

【0139】

ステップＳ３３３３では、線形パラメータ推定部３３３が、状態量算出部３３２から出力された非線形パラメータ要素毎の摂動値それぞれに対する状態量に基づき、非線形パラメータ要素毎の摂動値それぞれに対する線形パラメータψおよび蓄電池モデルの出力値を推定する。そして、上述したステップＳ２３３１以下の各ステップの実行を継続する。

【0140】

以上が、実施の形態1にかかるモデルパラメータ推定装置３の処理手順、つまりモデルパラメータ推定方法における動作の一例である。なお、図５で説明した処理の内容および順序は一例であり、この内容および順序に限定されない。

【0141】

つぎに、実施の形態１にかかるモデルパラメータ推定装置３あるいはモデルパラメータ推定方法をある蓄電池システム１００の動作時の電流と電圧の時系列データに対して適用した結果について、図６を用いて説明する。なおこの適用では、直並列接続構成による蓄電池システム１００における電流と電圧を蓄電池１としてセル換算している。

【0142】

ここでは、蓄電池１のモデルとして図４で示したフォスター型回路で近似し、Ｒ_１、Ｒ_２、Ｒ_３およびＣ_１、Ｃ_２、Ｃ_３は拡散インピーダンスを有限近似する式（２９）に従うものとした。また、ＯＣＶ関数ｆ_ＯＣＶとしてはｍ_ｓ＝１０で等間隔に分割した３次スプライン曲線を用いた。簡単のため、非線形パラメータφはτ_ｄのみとし、線形パラメータψはＲ_０、Ｃ_ｄおよび３次スプライン曲線の全パラメータとした。

【0143】

図６上段の電流値Ｉの時系列データでは、前半が自由な充放電、後半が定電力充電となっている。いずれにおいても中段の実データ（計測値：実線）と下段のモデル出力（推定値：破線）が精度よく一致していることがわかる。

【0144】

また、このとき推定されたＯＣＶ関数ｆ_ＯＣＶを図７に示す。ただし、横軸にはある規格容量ｑ_ｔｙｐによって規格化した電気量ｑ／ｑ_ｔｙｐを用いている。推定結果は各点をなめらかにつなぎつつ、完全に線形ではない電気量ｑ／ｑ_ｔｙｐに対するＯＣＶの非線形な変動をモデル化できている。この効果により、図６の電圧（中段：実データ、下段：モデル出力）に関し、充電中の微妙な凹凸形状を含めて精度よく推定できている。

【0145】

実施の形態２．
実施の形態１においては、ＯＣＶ関数が未知の場合の処理について説明した。本実施の形態２においては、ＯＣＶ関数が既知である場合について説明する。なお、モデルパラメータ推定装置の構成、およびモデルパラメータ推定方法の基本動作は実施の形態１と同様であり、同様部分の説明を省略するとともに実施の形態１で用いた図を援用する。

【0146】

ＯＣＶ関数ｆ_ＯＣＶが既知の場合、式（２３）～式（２５）で表わされる蓄電池１の状態空間モデルにおいて、式（２５）は満充電容量ｑ_ｍａｘを用いて式（５５）のように表すことができる。

【0147】

【数38】

このとき、ＯＣＶ関数ｆ_ＯＣＶは既知であるため、非線形パラメータφは式（５６）のようになる。
φ＝［τ_ｄ Ｉ_ｏｆｆ ｑ_ｍａｘ］^Τ ・・（５６）

【0148】

また、出力方程式を変換すれば、式（５７）のようになり、実施の形態1と同様にパラメータ線形で表現できるから、線形パラメータψは式（５８）のようになる。

【0149】

【数39】

この場合、実施の形態1と異なり、満充電容量ｑ_ｍａｘも含めて推定することができるという利点がある。

【0150】

また、ＯＣＶ関数ｆ_ＯＣＶのより詳細な情報として、正極電位関数ｆ_ｐと負極電位関数ｆ_ｎが既知であり、ＯＣＶ関数ｆ_ＯＣＶが電極パラメータを用いて式（５９）のように表わせるとする。

【0151】

【数40】

ただし、電極パラメータのうち、ｑ_ｐ，０とｑ_ｎ，０はそれぞれ正極と負極の初期充電量であり、ｑ_{ｐ，ｍａｘ}とｑ_{ｎ，ｍａｘ}はそれぞれ正極と負極の満充電容量である。このように表現したとき、非線形パラメータφは式（６０）のようになる。
φ＝［τ_ｄ Ｉ_ｏｆｆ ｑ_ｐ，０ ｑ_{ｐ，ｍａｘ} ｑ_ｎ，ｏ ｑ_{ｎ，ｍａｘ}］^Τ ・・（６０）
線形パラメータψについては式（５８）と同様である。

【0152】

この場合、ＯＣＶ関数ｆ_ＯＣＶをより詳細にモデル化し、蓄電池１の劣化をより詳細に診断することができる。具体的には、蓄電池１の劣化に伴う正極容量の低下、負極容量の低下、正極と負極の電気量のずれを区別して把握することができる。なお、正極と負極の電気量のずれは、おもに負極表面に形成される被膜（ＳＥＩ：Solid Electrolyte Interface）の成長とリチウムの析出に由来して発生する現象と言われており、蓄電池１の満充電容量の低下を引き起こすものである。

【0153】

上記のように本願のモデルパラメータ推定装置３あるいはモデルパラメータ推定方法では、非線形パラメータφの更新と線形パラメータψの推定を分離する。あるいは問題Ｐ２のように線形パラメータの線形最小二乗解が埋め込まれた問題を考えて非線形パラメータφを更新していくようにした。これによる最適化計算時の利点は、上述したように大別して３つ挙げられる。

【0154】

第一に、初期値設定が非線形パラメータφのみでよいことによる、初期値依存性の低下が挙げられる。第二に、計算量が削減されることが挙げられる。第三に、数値的安定性の向上が挙げられる。

【0155】

線形パラメータψの推定においては、制約なし線形最小二乗問題と線形等式制約付き線形最小二乗問題と線形不等式制約問題のいずれにおいても最適解を得る事が可能である。とくに、制約なし線形最小二乗問題においては、線形最小二乗解が繰り返し計算なしで得られる。また、線形等式制約付き線形最小二乗問題においても、ラグランジュ未定乗数法により大域的な最適解が繰り返し計算なしで得られる。

【0156】

こうした利点は、蓄電池モデルを複雑化し、推定すべきパラメータ数が増えた場合においてより生かされることになる。それゆえ、事前に保有するＯＣＶ特性として、正極電位特性と負極電位特性を別々に保有したうえで、電極パラメータの一部または全部を推定することにより、より詳細なモデル化および劣化診断をするような場合にも有用となる。

【0157】

ここで特許文献１に記載された出力誤差法に基づく連続時間システム同定法について再度検討すると、ＯＣＶ特性データが必須とはされていないものの、ＯＣＶ特性データを保有していない場合のモデルパラメータ推定方法は明記されていない。実用においては、ＯＣＶ特性データを仮定しないモデルパラメータ推定技術も重要となる。例えば、ＯＣＶ特性は劣化などにより変動することが知られているため、推定の高精度化と劣化診断への応用のためにはＯＣＶ特性を含むモデルパラメータの推定が重要となる。

【0158】

また、蓄電池のリユースにおいては、ＯＣＶ特性を含むモデルパラメータの真値に関する事前情報がないまたは少ない状況下でのモデル化が必要となりうる。ＯＣＶ特性をある関数でパラメータ表現し、モデルパラメータに含めて一括推定する場合、特許文献１に記載の方法では上述の初期値設定と収束性の問題があるため、最適値に収束させることがより難しくなる。

【0159】

これに対し、本願では、電気量（ＳＯＣまたは規格化された電気量であってもよい）とＯＣＶとの関係であるＯＣＶ特性に関し、パラメータを含む関数で表現したうえで該パラメータを推定することで、ＯＣＶ特性が未知であってもＯＣＶ特性を推定することが可能となる。

【0160】

とくに、多項式、区分多項式などの多数の線形パラメータψを含む関数で表現すれば、非線形パラメータφの更新と線形パラメータψの推定のプロセスを分離した本願の特徴が活かされる。これにより、多数の線形パラメータψの初期値を繊細に設定する必要なく、全パラメータの推定を容易に実行可能となる。また上述のとおり、線形パラメータψの推定においては線形等式制約が含まれていても大域的な解が得られるため、線形等式制約を含む区分多項式としてＯＣＶ関数ｆ_ＯＣＶを扱った場合にも容易にパラメータ推定することが可能である。なお、多項式ではなく区分多項式を用いることで、曲線の局所的な変動も含めて精度良く表現しやすくなる。

【0161】

なお、本願は、様々な例示的な実施の形態および実施例が記載されているが、１つ、または複数の実施の形態に記載された様々な特徴、態様、および機能は特定の実施の形態の適用に限られるのではなく、単独で、または様々な組み合わせで実施の形態に適用可能である。従って、例示されていない無数の変形例が、本願明細書に開示される技術の範囲内において想定される。例えば、少なくとも１つの構成要素を変形する場合、追加する場合または省略する場合、さらには、少なくとも１つの構成要素を抽出し、他の実施の形態の構成要素と組み合わせる場合が含まれるものとする。

【0162】

以上のように、本願のモデルパラメータ推定装置３によれば、非線形パラメータφと線形パラメータψを用いて蓄電池１を表現した状態空間モデルに対して非線形パラメータφに値（例えば、初期値）を代入して得られた状態方程式と、蓄電池１の電流（電流値Ｉ）と端子電圧（電圧値Ｖ）それぞれの測定値の時系列データに基づき、電流値Ｉの測定値に対する蓄電池１の状態（例えば、電気量ｑ、電荷ｑ_ｄ）を示す状態量を算出する状態量算出部３３２、状態空間モデルと状態量と電流値Ｉの測定値に基づき算出される電圧値Ｖの推定値と電圧値Ｖの測定値との誤差を最小化する線形パラメータψを推定する線形パラメータ推定部３３３、推定した線形パラメータψが収束しているか否かを判定する判定部３３４、および判定部３３４が収束したと判定するまで、最小化された誤差が小さくなるように、所定の収束条件を満たすまで非線形パラメータφの値の更新を繰り返す非線形パラメータ更新部３３１、を備えるように構成した。そのため、初期値設定が非線形パラメータφのみに限られることで、少ない情報でパラメータを推定することができる。

【0163】

さらに、初期値依存性が低下でき、計算量も削減され、数値的安定性も向上する。その際、初期値としては蓄電池に対して一般的に知られた妥当性の高い値を用いることができるので、システムを更新する際にも追加情報なしに推定が可能となる。

【0164】

とくに、非線形パラメータ更新部３３１は、電圧値Ｖの推定値の非線形パラメータφに対する勾配または勾配近似値を利用して、あるいは勾配フリーの非線形最適化手法を用いて、非線形パラメータφの値を更新するようにすれば、容易に収束させることができる。

【0165】

あるいは、非線形パラメータ更新部３３１は、非線形パラメータφの値の更新にあたり、当該更新前の値である前回値を摂動させた非線形パラメータφそれぞれの摂動値を生成し、状態量算出部３３２は、それぞれの摂動値に対する状態量を算出し、線形パラメータ推定部３３３は、それぞれの摂動値に対する状態量に対する線形パラメータψを推定し、非線形パラメータ更新部３３１は、少なくともそれぞれの摂動値に対する状態量に対する線形パラメータψの推定値に基づき、端子電圧（電圧値Ｖ）の推定値の非線形パラメータφに対する勾配近似値を算出し、算出した勾配近似値に基づき前回値を更新するように構成すれば、より確実に収束する。

【0166】

それらの際、非線形パラメータφは、蓄電池１の充電率と満充電容量ｑ_ｍａｘと拡散時定数τ_ｄと電流（電流値Ｉ）の測定に用いるセンサのオフセット電流（電流オフセット誤差Ｉ_ｏｆｆ）のうちの少なくとも1つを含み、線形パラメータψは、蓄電池１の直流抵抗と拡散抵抗Ｒ_ｄとコンデンサ容量のうちの少なくとも１つを含むようにすれば、蓄電池１の状態を確実に評価できる。

【0167】

それらの際、線形パラメータ推定部３３３は、線形最小二乗問題の線形最小二乗解を用いて線形パラメータψを推定するようにすれば、計算量が少なくなる。

【0168】

その際、状態空間モデルには、電気量ｑと開回路電圧との関係を表わし、線形パラメータψを含むＯＣＶ関数ｆ_ＯＣＶを含むようにすれば、正確に蓄電池１をモデル化できる。

【0169】

さらに、ＯＣＶ関数ｆ_ＯＣＶは電気量ｑの多項式であり、線形パラメータψは、多項式の係数を含むようにすれば、電気量ｑに対するＯＣＶの非線形な変動をモデル化できる。

【0170】

状態空間モデルは、線形パラメータψに対する線形等式制約を含み、線形パラメータ推定部３３３は、線形等式制約を含む線形最小二乗問題をラグランジュ未定乗数法により解くことで線形パラメータψを推定するように構成すれば、非線形パラメータφの更新と、ラグランジュ未定乗数法による線形パラメータψの算出を交互に繰り返すことで、確実に収束可能である。

【0171】

状態空間モデルには、電気量ｑと開回路電圧との関係を表わし、線形パラメータψを含む区分多項式のＯＣＶ関数ｆ_ＯＣＶを含み、線形パラメータψは、区分多項式の線形パラメータを含み、線形等式制約は、区分多項式の節点における線形等式制約を含むようにすれば、非線形パラメータφの更新と線形パラメータψの推定のプロセスを分離した本願の特徴を十分に活かすことができる。

【0172】

あるいは、状態空間モデルは、線形パラメータψに対する線形不等式制約を含み、線形パラメータ推定部３３３は、線形不等式制約を含む線形最小二乗問題を解くことで線形パラメータψを推定するように構成しても、最適解を得る事が可能である。

【0173】

また、状態空間モデルは、蓄電池１の充電率と開回路電圧との関係を表わすＯＣＶ特性情報を含み、非線形パラメータφは、満充電容量ｑ_ｍａｘを含むようにすれば、蓄電池１の劣化を詳細に診断することができる。

【0174】

その際、状態空間モデルは、蓄電池１の正極電気量と正極電位との関係を表わす正極電位特性情報と、蓄電池１の負極電気量と負極電位との関係を表わす負極電位特性情報とを含み、非線形パラメータφは、満充電容量ｑ_ｍａｘにかえて、正極容量と負極容量と正極電気量と負極電気量のうちの少なくとも１つを含むように構成すれば、蓄電池１の劣化をより詳細に診断することができる。

【0175】

また、本願のモデルパラメータ推定方法によれば、非線形パラメータφと線形パラメータψを用いて蓄電池１を表現した状態空間モデルの非線形パラメータφの初期値を設定するステップ（Ｓ１３２０）、状態空間モデルに対して非線形パラメータφに値を代入して得られた状態方程式と、蓄電池１の電流（電流値Ｉ）と端子電圧（電圧値Ｖ）それぞれの測定値の時系列データに基づき、電流値Ｉの測定値に対する蓄電池の状態（例えば、電気量ｑ、電荷ｑ_ｄ）を示す状態量を算出する状態量算出ステップ（Ｓ２３３２）、状態空間モデルと状態量と電流値Ｉの測定値に基づき算出される電圧値Ｖの推定値と電圧値Ｖの測定値との誤差を最小化する線形パラメータψを推定する線形パラメータ推定ステップ（Ｓ２３３３）、推定した線形パラメータψが収束しているか否かを判定する判定ステップ（Ｓ２３３４）、および線形パラメータψが収束したと判定されるまで、最小化された誤差が小さくなるように、所定の収束条件を満たすまで非線形パラメータφの値の更新を繰り返す非線形パラメータ更新ステップ（Ｓ２３３１）、を含むように構成した。そのため、初期値設定が非線形パラメータφのみに限られることで、少ない情報でパラメータを推定することができる。

【0176】

さらに、初期値依存性が低下でき、計算量も削減され、数値的安定性も向上する。その際、初期値としては蓄電池に対して一般的に知られた妥当性の高い値を用いることができるので、システムを更新する際にも追加情報なしに推定が可能となる。

【0177】

非線形パラメータ更新ステップ（Ｓ２３３１）では、電圧値Ｖの推定値の非線形パラメータφに対する勾配または勾配近似値を利用して、あるいは勾配フリーの非線形最適化手法を用いて、非線形パラメータφの値を更新するようにすれば、容易に収束させることができる。

【0178】

あるいは、非線形パラメータφの値の更新にあたり、当該更新前の値である前回値を摂動させた非線形パラメータφ（の要素）それぞれの摂動値を生成する摂動値設定ステップ(Ｓ３３３１)、それぞれの摂動値に対する状態量を算出する摂動毎状態量算出ステップ（Ｓ３３３２）、およびそれぞれの摂動値に対する状態量に対する線形パラメータψを推定する摂動毎線形パラメータ推定ステップ（Ｓ３３３３）、をさらに含み、非線形パラメータ更新ステップ（Ｓ２３３１）では、少なくともそれぞれの摂動値に対する状態量に対する線形パラメータの推定値に基づき、端子電圧(電圧値Ｖ)の推定値の非線形パラメータφに対する勾配近似値を算出し、算出した勾配近似値に基づき前回値を更新するように構成すれば、より確実に収束する。

【0179】

それらの際、非線形パラメータφは、蓄電池１の充電率と満充電容量ｑ_ｍａｘと拡散時定数τ_ｄと電流（電流値Ｉ）の測定に用いるセンサのオフセット電流（電流オフセット誤差Ｉ_ｏｆｆ）のうちの少なくとも1つを含み、線形パラメータψは、蓄電池１の直流抵抗と拡散抵抗Ｒ_ｄとコンデンサ容量のうちの少なくとも１つを含むようにすれば、蓄電池１の状態を確実に評価できる。

【0180】

線形パラメータ推定ステップ（Ｓ２３３３）では、線形最小二乗問題の線形最小二乗解を用いて線形パラメータψを推定するようにすれば、計算量が少なくなる。

【0181】

その際、状態空間モデルには、電気量ｑと開回路電圧との関係を表わし、線形パラメータψを含むＯＣＶ関数ｆ_ＯＣＶを含むようにすれば、正確に蓄電池１をモデル化できる。

【0182】

さらに、ＯＣＶ関数ｆ_ＯＣＶは電気量ｑの多項式であり、線形パラメータψは、多項式の係数を含むようにすれば、電気量ｑに対するＯＣＶの非線形な変動をモデル化できる。

【0183】

状態空間モデルは、線形パラメータψに対する線形等式制約を含み、線形パラメータ推定ステップ（Ｓ２３３３）では、線形等式制約を含む線形最小二乗問題をラグランジュ未定乗数法により解くことで線形パラメータψを推定するように構成すれば、非線形パラメータφの更新と、ラグランジュ未定乗数法による線形パラメータψの算出を交互に繰り返すことで、確実に収束可能である。

【0184】

状態空間モデルには、電気量ｑと開回路電圧との関係を表わし、線形パラメータψを含む区分多項式のＯＣＶ関数ｆ_ＯＣＶを含み、線形パラメータψは、区分多項式の線形パラメータを含み、線形等式制約は、区分多項式の節点における線形等式制約を含むようにすれば、非線形パラメータφの更新と線形パラメータψの推定のプロセスを分離した本願の特徴を十分に活かすことができる。

【0185】

あるいは状態空間モデルは、線形パラメータψに対する線形不等式制約を含み、線形パラメータ推定ステップ（Ｓ２３３３）では、線形不等式制約を含む線形最小二乗問題を解くことで線形パラメータψを推定するように構成しても、最適解を得る事が可能である。

【0186】

また、状態空間モデルは、蓄電池１の充電率と開回路電圧との関係を表わすＯＣＶ特性情報を含み、非線形パラメータφは、満充電容量ｑ_ｍａｘを含むようにすれば、蓄電池１の劣化を詳細に診断することができる。

【0187】

その際、状態空間モデルは、蓄電池１の正極電気量と正極電位との関係を表わす正極電位特性情報と、蓄電池１の負極電気量と負極電位との関係を表わす負極電位特性情報とを含み、非線形パラメータφは、満充電容量ｑ_ｍａｘにかえて、正極容量と負極容量と正極電気量と負極電気量のうちの少なくとも１つを含むように構成すれば、蓄電池１の劣化をより詳細に診断することができる。

【符号の説明】

【0188】

１：蓄電池、 ２１：電流検出装置、 ２２：電圧検出装置、 ３：モデルパラメータ推定装置、 ３１：時系列データ取得部、 ３２：モデル取得部、 ３３：パラメータ推定部、 ３３１：非線形パラメータ更新部、 ３３２：状態量算出部、 ３３３：線形パラメータ推定部、 ３３４：判定部、 Ｉ：電流値、 Ｉ_ｏｆｆ：電流オフセット誤差、 ｑ：電気量、 ｑ_ｍａｘ：満充電容量、 Ｒ_ｄ：拡散抵抗、 Ｖ：電圧値、 τ_ｄ：拡散時定数。

【図1】

【図2】

【図3】

【図4】

【図5】

【図6】

【図7】