WO2020208762A1

WO2020208762A1 - 蓄電池の特性推定装置および特性推定方法

Info

Publication number: WO2020208762A1
Application number: PCT/JP2019/015727
Authority: WO
Inventors: 智己竹上; 敏裕和田
Original assignee: 三菱電機株式会社
Priority date: 2019-04-11
Filing date: 2019-04-11
Publication date: 2020-10-15
Also published as: US20220155375A1; DE112019007176T5; US12072383B2; CN113661399B; CN113661399A; JPWO2020208762A1; JP7069412B2

Abstract

蓄電池の特性推定装置は、蓄電池の端子電圧を検出する電圧検出部と、蓄電池の休止状態における端子電圧の時系列データを取得する時系列データ取得部と、蓄電池モデルを提供するモデル提供部と、時系列データ取得部により取得された端子電圧の時系列データと、モデル提供部により提供された蓄電池モデルとに基づいて、蓄電池モデルのモデルパラメータを推定する、モデルパラメータ推定部とを備えている。蓄電池モデルは、蓄電池のＯＣＶを表現するＯＣＶ項、および蓄電池の電極を構成する粒子間のイオン拡散を表現する１次元の拡散方程式の基本解に基づく粒子間拡散項を含んでいる。

Description

蓄電池の特性推定装置および特性推定方法

　本発明は、蓄電池の特性推定装置および特性推定方法に係り、特に蓄電池の電圧緩和過程の推定装置および推定方法に関する。

　蓄電池の特性を表現する様々なモデルが提案されている。例えば、蓄電池の電圧緩和過程を表現するモデルとしては、ＯＣＶを含む一般的な蓄電池の等価回路をベースとした、ＣＲ並列素子の多直列モデルが知られている（例えば、特許文献１参照）。

　特許文献１では、蓄電池の電気化学的な応答を一次遅れ要素の重ね合わせによって表現し、ＣＲ並列素子の休止時の応答に対応する指数関数項によって、蓄電池の電圧緩和過程をモデル化している。指数関数項を表現するモデルパラメータは、蓄電池の休止状態において検出される電圧緩和過程の時系列データに基づいて推定される。

特許第５１７０８７６号

　しかしながら、特許文献１に記載のモデルによって、蓄電池の電圧緩和過程を高精度に推定することには限界がある。例えば、蓄電池の休止状態が短時間である場合には、短い時定数の電圧緩和過程が支配的になるが、蓄電池の休止状態が長時間である場合には、より長い時定数の電圧緩和過程が支配的になる。そのため、特許文献１のモデルでは、蓄電池の短時間の休止状態において検出される時系列データに基づいて、蓄電池の長時間の休止状態における電圧緩和過程を表現するモデルパラメータを推定することは困難である。

　本発明は、上記のような課題を解決するためのものであり、蓄電池の電圧緩和過程を高精度に推定することができる、蓄電池の特性推定装置および特性推定方法を提供することを目的とする。

　上記の課題を解決するために、本発明に係る蓄電池の特性推定装置は、蓄電池の端子電圧を検出する電圧検出部と、蓄電池の休止状態における端子電圧の時系列データを取得する時系列データ取得部と、蓄電池モデルを提供するモデル提供部であって、当該蓄電池モデルは、蓄電池のＯＣＶを表現するＯＣＶ項、および蓄電池の電極を構成する粒子間のイオン拡散を表現する１次元の拡散方程式の基本解に基づく粒子間拡散項を含む、モデル提供部と、時系列データ取得部により取得された端子電圧の時系列データと、モデル提供部により提供された蓄電池モデルとに基づいて、蓄電池モデルのモデルパラメータを推定する、モデルパラメータ推定部とを備える。

　また、本発明に係る蓄電池の特性推定方法は、蓄電池の端子電圧を検出するステップと、蓄電池の休止状態における端子電圧の時系列データを取得するステップと、蓄電池モデルを提供するステップであって、当該蓄電池モデルは、蓄電池のＯＣＶを表現するＯＣＶ項、および蓄電池の電極を構成する粒子間のイオン拡散を表現する１次元の拡散方程式の基本解に基づく粒子間拡散項を含む、蓄電池モデルを提供するステップと、端子電圧の時系列データと、蓄電池モデルとに基づいて、蓄電池モデルのモデルパラメータを推定するステップとを含む。

　本発明に係る蓄電池の特性推定装置および特性推定方法によれば、蓄電池の電圧緩和過程を高精度に推定することができる。

本発明の実施の形態１に係る蓄電池の特性推定装置の構成を示すブロック図である。図１の特性推定装置における時系列データ取得部の内部構成を示すブロック図である。蓄電池内部の電極のイメージ図である。蓄電池内部の電極の単粒子モデルによるイメージとその等価回路を示す図である。蓄電池内部の電極の１次元近似モデルによるイメージとその等価回路を示す図である。１次元近似等価回路モデルの等価回路を示す図である。図６の等価回路を近似変換した等価回路を示す図である。図７の等価回路をさらに近似変換した等価回路を示す図である。図８の等価回路をさらに近似変換した等価回路を示す図である。図８から図９の近似変換において、分布定数回路から拡散方程式が導出されることを説明する図である。電圧緩和過程における電位分布のイメージ図である。正極と負極を有する蓄電池の内部の１次元等価回路を示す図である。図１の特性推定装置におけるモデルパラメータ推定部の内部構成を示すブロック図である。半電池の電位緩和曲線を示す図である。半電池で取得した３種類のＳＯＣ－ＯＣＰデータを示す図である。図１５の全データに対して、本発明の実施の形態１に係る手法を適用して、ＯＣＰを推定した結果が示す図である。本発明の実施の形態１に係る特性推定装置の各機能を専用のハードウェアである処理回路で実現する場合を示した構成図である。本発明の実施の形態１に係る特性推定装置の各機能をプロセッサおよびメモリを備えた処理回路より実現する場合を示した構成図である。

　以下、添付図面を参照して、本願が開示する発明の実施の形態について、詳細に説明する。ただし、以下に示す実施の形態は一例であり、これらの実施の形態によって、本発明が限定されるものではない。

　実施の形態１．
　図１は、本発明の実施の形態１に係る蓄電池の特性推定装置１００の構成を示す図である。なお、図１には、特性推定装置１００に接続される蓄電池１０１も併せて示されている。

　本実施の形態１では、蓄電池１０１として、リチウムイオン蓄電池を想定する。ただし、蓄電池１０１は、鉛蓄電池、ニッケル水素蓄電池、全個体電池、ナトリウム硫黄電池、ナトリウム硫黄電池、マグネシウムイオン電池、電気二重層等であってもよい。

　特性推定装置１００は、電圧検出部１０２と、時系列データ取得部１０３と、モデル提供部１０４と、モデルパラメータ推定部１０５とを備えている。

　電圧検出部１０２は、蓄電池１０１の端子電圧Ｖを検出する。

　時系列データ取得部１０３は、電圧検出部１０２によって検出される蓄電池１０１の端子電圧Ｖに基づいて、蓄電池１０１の休止状態における端子電圧Ｖの時系列データ｛（ｐ，Ｖｐ）｜ｐ＝０，１，・・・ｍ－１｝を取得する。

　より詳細には、図２に示されるように、時系列データ取得部１０３は、電流検出部１０３ａと、休止状態検出部１０３ｂと、時系列データ記憶部１０３ｃと、時系列データ抽出部１０３ｄとを含んでいる。

　電流検出部１０３ａは、蓄電池１０１の電流Ｉを検出する。

　休止状態検出部１０３ｂは、電流検出部１０３ａによって検出される蓄電池１０１の電流Ｉに基づいて、蓄電池１０１の休止状態を検出する。ただし、休止状態検出部１０３ｂは、蓄電池１０１が必ずしも完全に休止している状態のみを検出する必要はなく、休止に近い状態を含めて検出してもよい。これ以降、蓄電池１０１の休止に近い状態を「準休止状態」と称することにする。

　蓄電池１０１の休止状態を検出する方法としては、例えば、予め決定された微小値εおよびδに対して、｜Ｉ｜＜εかつ｜ΔＩ｜＜δを満たしている状態を検出する方法を採用することができる。ここで、ΔＩは、ある時刻ｔにおける電流Ｉとその直前の時刻ｔ－Δｔにおける電流Ｉとの差分値である。

　この検出方法では、第１の不等式｜Ｉ｜＜εによって、電流検出部１０３ａの検出誤差に起因して電流Ｉが厳密にゼロとならない場合でも、休止状態を検出することができる。また、εを大きく設定することによって、上述の「準休止状態」を検出することもできる。

　また、第２の不等式｜ΔＩ｜＜δによって、電流Ｉの大きさだけでなく、電流Ｉの変化も考慮に入れるため、休止状態をより高精度に検出することができる。

　なお、休止状態検出部１０３ｂは、蓄電池１０１を搭載するシステムが蓄電池１０１の休止状態を別途検出している場合には、その検出信号を受信することによって、蓄電池１０１の状態時間を検出してもよい。

　時系列データ記憶部１０３ｃは、休止状態検出部１０３ｂによって蓄電池１０１の休止状態が検出されている間、電圧検出部１０２によって検出される蓄電池１０１の端子電圧Ｖの時系列データ｛（ｋ，Ｖｋ）｜ｋ＝０，１，・・・ｎ｝を記憶する。ここで、ｋはサンプリング時刻であり、サンプリング周期をｔｓとすると、時刻ｔ＝ｋ＊ｔｓの関係が成立する。なお、時系列データ記憶部１０３ｃは、必ずしも時系列データ取得部１０３の内部に設けられる必要はない。時系列データ記憶部１０３ｃは、例えばクラウド上等に設けられてもよい。

　なお、端子電圧Ｖの時系列データを記憶する際のサンプリング周期ｔｓは、時系列データ記憶部１０３ｃが決定してもよいし、電圧検出部１０２が端子電圧Ｖを検出する際に決定してもよい。

　時系列データ抽出部１０３ｄは、時系列データ記憶部１０３ｃに記憶されている時系列データ｛（ｋ，Ｖｋ）｜ｋ＝０，１・・・，ｎ－１｝から、不必要なデータを除外して、必要な時系列データ｛（ｐ，Ｖｐ）｜ｐ＝０，１，・・・，ｍ－１｝を抽出する。

　ここで、ｐはｋの時系列データから抽出された値であることから、ｐ＝０，１，・・・，ｍ－１は、ｔ＝ｒ₀＊ｔｓ，ｒ₁＊ｔｓ，・・・，ｒ_m-1＊ｔｓに対応する。すなわち、あるｐに対して、ｔ＝ｒｐ＊ｔｓとなる。ただし、ｒｐは抽出されたｐ個目の時系列データである。

　時系列データ抽出部１０３ｄにおける具体的な抽出方法としては、例えば、元の時系列データから、最初の一定時間のデータを除外することが考えられる。これは蓄電池１０１の休止状態の初期には、過去の充放電履歴に依存して蓄電池１０１内部のイオン分布に乱れが存在すると考えられるため、電圧緩和過程の推定に適さないデータとなる可能性があるからである。

　また、明らかなデータの飛びが存在する場合、これを除外してもよい。また、電圧変化に基づくサンプリングを行って、電圧変化の小さいデータを除外してもよい。

　このようにして、元の時系列データ｛（ｋ，Ｖｋ）｜ｋ＝０，１，・・・，ｎ－１｝から、不必要なデータを除外した時系列データ｛（ｐ，Ｖｐ）｜ｐ＝０，１，・・・，ｍ－１｝を抽出する。ただし、このようなデータ抽出の処理は、特性の推定方法、保存可能なデータ数等の制約に依存するものであり、必須のものではない。

　モデル提供部１０４は、拡散方程式の基本解に基づく蓄電池モデルを少なくとも１つ提供する。これ以降、このようなモデルを「拡散基本解モデル」と称することにする。

　モデルパラメータ推定部１０５は、時系列データ取得部１０４により取得された、蓄電池１０１の端子電圧Ｖの時系列データ｛（ｐ，Ｖｐ）｜ｐ＝０，１，・・・，ｍ－１｝と、モデル提供部１０４により提供された蓄電池モデルとに基づいて、蓄電池モデルのモデルパラメータを推定する。

　［蓄電池モデルの技術的な詳細］
　以下、本実施の形態１に係る蓄電池モデルの技術的な詳細について説明する。

　図３は、蓄電池内部の電極のイメージ図である。集電体３０１に多数の粒子３０２が固着してバルク３０３を形成しており、それらの全体が電解液３０４に浸されている。なお、正確には、電子伝導性を高めるための導電助剤、多数の粒子３０２を結着させるためのバインダー等も付着しているが、ここでは簡単のため省略している。

　蓄電池の充電時には、個々の粒子３０２に対して、電解液３０４からリチウムイオン３０５が供給されると共に、集電体３０１から電子３０６が供給される。蓄電池の放電時には、個々の粒子３０２から、リチウムイオン３０５および電子３０６が放出される。

　なお、先述したように、本実施の形態１では、蓄電池としてリチウムイオン蓄電池を想定している。そのため、図３は、リチウムイオン蓄電池の場合のイメージ図である。しかしながら、蓄電池は、鉛蓄電池、ニッケル水素蓄電池、全個体電池、ナトリウム硫黄電池、ナトリウム硫黄電池、マグネシウムイオン電池、電気二重層等であってもよい。

　例えば、蓄電池が全個体電池である場合、上述の「電解液」を「固体電解質」に読み替える。また、マグネシウムイオン電池の場合、「リチウムイオン」を「マグネシウムイオン」に読み替える。

　図４は、蓄電池内部の電極の単粒子モデルによるイメージと、その等価回路を示す図である。図４の等価回路において、Ｒｅは電子伝導抵抗、Ｒｃｔは反応抵抗、Ｃｄｌは電気二重層容量、Ｚｗは粒子３０７内部の拡散インピーダンス、ＯＣＶは粒子３０７の開回路電圧、Ｒｉはイオン伝導抵抗である。

　図４の単粒子モデルは、本来は図３のように多数存在する粒子３０２が全て同じ挙動を示すと仮定して、全粒子の挙動を単一粒子３０７に集約した理想モデルと捉えることができる。このような単粒子モデルは、蓄電池のモデルとして広く用いられているが、粒子間のリチウムイオン分布の偏りを表現できないという欠点を有する。

　図５は、蓄電池内部の電極の１次元近似モデルによるイメージ図と、その等価回路を示す図である。図５では簡単のため、粒子数を４で近似表現した場合を示している。この１次元近似モデルでは、図４の単粒子モデルとは異なり、粒子間のリチウムイオン分布の偏りを表現することが可能である。

　図５の等価回路は、図４の単粒子モデルの等価回路が並列に多段接続されたものとなっている。ただし、各粒子３０７の間に電子伝導抵抗とイオン伝導抵抗が存在するため、集電体３０１に近い粒子３０７ほど、電子伝導抵抗が小さくなり、電解液中のイオン伝導抵抗は大きくなる。反対に、集電体３０１から遠い粒子３０７ほど、電子伝導抵抗が大きくなり、電解液中のイオン伝導抵抗は小さくなる。

　このような電子伝導抵抗とイオン伝導抵抗の影響により、電流レートが大きくなるほど、各粒子３０７に供給されるイオン量に差がつき、粒子間でＳＯＣ分布の偏りが生じる。従来の蓄電池の休止状態における電圧緩和過程を推定する技術では、このような粒子間のＳＯＣ分布の偏りの影響を反映しながら実用に耐えるモデルが構築されてこなかった。そのため、粒子間のＳＯＣ分布の偏りの解消過程を反映した長時間の電圧緩和過程を表現することが困難であった。

　これに対して、本実施の形態１では、図５の１次元近似モデルから導出された、粒子間のＳＯＣ分布の偏りの解消過程を反映した長時間の電圧緩和過程を高精度に表現するモデルを用いる。

　なお、以降の説明において「イオン」とは、蓄電池の充放電に寄与するイオンのことである。例えば、リチウムイオン電池の場合はリチウムイオン、マグネシウムイオン電池の場合はマグネシウムイオン、ナトリウムイオン電池の場合はナトリウムイオンである。

　［１次元近似モデルの等価回路の近似変換］
　まず、図６の１次元近似モデルの等価回路を近似変換すると、図７のような等価回路が得られる。ただし、蓄電池は休止状態にあり、かつ休止状態となってから数秒～数百秒が経過しているものとする。また、電極電位ｖは、集電体側の粒子の電位に基づいている。

　図６から図７の近似変換では、電子伝導抵抗Ｒｅはイオン伝導抵抗Ｒｉと比較して非常に小さいため、これらを無視している。また、粒子内拡散を表現するワールブルグインピーダンスＺｗを分離して直列接続している。ワールブルグインピーダンスＺｗによって表現される粒子内拡散の過程は、数秒～数百秒程度の長い時定数を有する。これに対して、反応抵抗Ｒｃｔと電気二重層容量Ｃｄｌによる反応過程の時定数τ：＝Ｒｃｔ＊Ｃｄｌは、数十～数百ミリ秒程度の短い時定数である。そのため、ワールブルグインピーダンスＺｗを近似的に分離して、直列接続可能であることが一般的に知られている。

　次に、図７の等価回路をさらに近似変換すると、図８のような等価回路が得られる。図７から図８の近似変換では、各粒子の反応過程における過電圧を無視し、かつ粒子内拡散を表現するワールブルグインピーダンスＺｗを１つに統合して近似表現している。

　蓄電池が休止状態となった直後は、過去の充放電履歴に依存して各粒子の過電圧にばらつきが存在すると考えられる。しかしながら、反応過程の過電圧は、反応過程の時定数に従って、数十～数百ミリ秒程度でほぼゼロに収束すると考えられる。

　また、粒子内拡散の過電圧についても、粒子内拡散の時定数に対応した数秒～数百秒程度の時間が経過して粒子内拡散が収束に近づくと、粒子間の過電圧のばらつきは小さくなっていく。そのため、そこから先の過電圧の変化は、１つのワールブルグインピーダンスで代表して近似表現可能であると考えられる。

　また、各粒子のＯＣＶについては、粒子ＳＯＣのばらつきに応じたＯＣＶのばらつきを表現可能であればよいので、線形近似してキャパシタで表現したものを用いれば十分であると考えられる。

　次に、図８の等価回路をさらに近似変換すると、図９のような等価回路が得られる。図９の等価回路は、図８の等価回路におけるＯＣＶを表現するキャパシタとイオン抵抗の梯子の連なりを分布定数回路とみなして、拡散インピーダンスで表現したものである。

　図１０は、図８から図９の近似変換において、分布定数回路から拡散方程式が導出されることを説明する図である。図１０は、分布定数回路の微小距離Δｘの一部分を表現する回路である。ただし、ＲとＣはそれぞれ単位長さ当たりの抵抗と容量であり、ＲΔｘとＣΔｘはそれぞれ微小距離Δｘ当たりの抵抗と容量を表す。

　まず、電流に関する等式を考えると、

が成り立つから、両辺をΔｘで割ってΔｘ→０とすると、

となる。

　同様に、電圧に関する等式を考えると、

が成り立つから、両辺をΔｘで割ってΔｘ→０とすると、

となる。ここで、両辺をさらにｘで偏微分すると、

となる。

　最後に、式（２）と式（５）から∂ｉ（ｘ，ｔ）／∂ｘを消去し、Ｄ：＝１／ＲＣとすると、

となる。これは拡散方程式にほかならない。

　以上のようにして、蓄電池の休止状態における電極電位モデルとして、図９に示すような拡散インピーダンスを直列接続したものが得られる。ただし、粒子内拡散は粒子が理想的には球状であることから３次元拡散に従い、粒子間拡散は深さｘの方向、すなわち集電体からセパレータへの方向の拡散であることから１次元拡散に従う。

　［電圧緩和過程の数理モデルの導出］
　続いて、上記の等価回路モデルから、電圧緩和過程の数理モデルを導出する。

　まず、－∞＜ｘ＜＋∞の領域において、初期分布をｖ（ｘ，０）＝ｖ₀（ｘ）、境界条件をｖ（±∞，ｔ）＝０とする。なお、実際のリチウムイオンの拡散を考えるとき、リチウムイオンの位置ｘの領域は０≦ｘ＜＋∞となるが、この場合も、０≦ｘ＜＋∞の領域と左右対称となるように－∞＜ｘ≦０の領域を仮想的にとれは、原点で常に濃度勾配０となる。そのため、原点を境にリチウムイオンの流出入はマクロな視点では生じないため、条件を－∞＜ｘ＜＋∞の領域とした場合の議論を同様に適用することが可能である。

　また、リチウムイオンの位置ｘの領域は、厳密にはある値ｘ_maxによって、有限領域０≦ｘ≦ｘ_maxに制限されるはずである。しかしながら、一般的に、拡散層が非常に厚い場合、領域を半無限の０≦ｘ＜＋∞とみなすことも多い。そのため、本実施の形態１でもこれに従う。このとき、式（６）の拡散方程式の解は、

と表現される。しかしながら、式（７）の解析解を得ることは困難であり、得られるのは特殊な初期分布の場合のみに限られる。

　そこで、初期分布がｖ₀（ｘ）＝δ（ｘ）となる場合を考える。ただし、δ（ｘ）はデルタ関数であり、

と定義されるようなものである。初期分布がデルタ関数に従うというのは、ｘ＝０において濃度１の物質を滴下した状況と捉えることが可能である。すると、拡散方程式の解は、

と表わされる。Ｇ（ｘ，ｔ）はグリーン関数とも呼ばれ、拡散方程式の基本解として知られている。

　これを一般化すると、１次元拡散方程式の基本解に基づく項は、パラメータｋと関数ｆを用いて

と記述することができる。ここで、ｋは大きさを調整するパラメータである。

　本実施の形態１では、電極電位、あるいはリチウムイオン濃度の初期分布を、一定時間ｚ秒経過後の基本解の分布で仮定する。すなわち、ｖの推定値をハット付きのｖで表現したとき、

とする。ここで、ｋは大きさを調整するパラメータ、ｚは時間移動パラメータである。

　これを一般化すると、時間移動パラメータをさらに含む１次元拡散方程式の基本解に基づく項は、

と記述できる。以下に、時間移動パラメータを導入した式（１１）、式（１２）を用いることの妥当性とその原理を説明する。

　まず、蓄電池の充電時、電極内部の電位分布ｖ（ｘ，ｔ）が任意の時刻ｔ＞０で不等式

を満たすような関数ｆ_minとｆ_maxを見つけることを考える。関数ｆ_maxの初期分布をｆ_0,maxとして、ｆ_maxとｖとの差を具体的に計算してみると、

となるから、ｆ_max（ｘ，ｔ）≧ｖ（ｘ，ｔ）が成り立つためには、ｆ_0,max（ｘ）≧ｖ₀（ｘ）が成り立つようなｆ_0,max（ｘ）を見つければよい。

　ｖ₀が有限であること、すなわちｖ₀≦ｖ_maxであることと、ｖ₀（∞）＝０であることを踏まえると、

として、ｆ_0,max（ｘ）≧ｖ₀（ｘ）を満たすように、Ｍとｔ_Mを選ぶことができる。

　また、同様にして、ｆ_minについても、

として、ｆ_0,min（ｘ）≦ｖ₀（ｘ）を満たすように、ｍとｔ_mを選ぶことができる。このとき、式（１３）は、

となる。

　式（１７）は、電極の電位分布が過去の充電電流に依存して波打っていたとしても、電極界面から深部に向けて漸近的に低下していく形状となっていることを意味している。

　さらに、電極界面、すなわちｘ＝０において、ｔ＞＞ｔ_M，ｔ_mであるような状況を考えると、

となる。誤差εは時間の経過とともに小さくなっていき、ｔ→∞でε→０となる。つまり、電位の緩和挙動を見ているため当然であるが、ＭＧ（０，ｔ＋ｔ_M）とｍＧ（０，ｔ＋ｔ_m）、およびこれら２つに挟まれたｖ（ｘ，ｔ）の値は、時間の経過とともに同じ値に収束してゆき、これらの値同士のずれは小さくなっていく。

　したがって、式（１７）の２つの基本解モデルの間に存在するような基本解モデルの集合を考え、それらの中で実際の電位挙動に最も一致するようなモデルを選び出すことができれば、過去に検出した電位挙動から将来の電位挙動を推定することが可能となる。

　図１１は、電圧緩和過程における電位分布のイメージ図である。実線は実際の電位分布、上の点線は上限基本解モデル、下の点線は下限基本解モデル、それらの間の破線は推定基本解モデルを表わしている。

　実際の電位分布および推定基本解モデルは、ともに深さｚ方向のすべての位置において上限基本解モデルと下限基本解モデルとの間に収まっており、長時間経過後もそれらの位置関係は変わっていない。そして、上限基本解モデルと下限基本解モデルは実際の分布との誤差が比較的大きいのに対して、推定基本解モデルは実際の分布の低周波挙動をよく捉えられていることが分かる。

　リチウムイオン濃度は過去の充放電履歴などに依存し、局所的に濃度に高低が存在する。しかしながら、マクロな傾向としては、電極界面は濃度が高く、深くなるにつれて濃度が低くなる、あるいはその逆となるはずである。

　本実施の形態１に係る基本解モデルは、このような大域的な挙動を捉えたものとなっている。なお、局所的な電位の高低、すなわちイオン濃度の濃淡は、拡散現象の進行とともに消失していくものであるため、実測データが長時間であるほど初期分布の影響を受けにくくなり、モデルの精度は高くなると考えられる。

　式（１１）の推定基本解モデルの別の視点での解釈として、電極内部のイオン濃度を、z秒過去の時点で電極界面にもたらされた高濃度のイオンが拡散するときの分布で近似表現していると捉えることもできる。つまり、過去の充放電履歴に応じた局所的なイオン分布の濃淡が時間の経過とともに解消されていった後に最後に残るような、最も支配的で低周波であるような拡散過程を、パラメータｋとｚによって、大きさと位置を適切にずらした基本解によって、高精度に捉えることができる。

　ここまでの議論と同様にして、電極粒子内部の拡散についても数理モデルを導出する。電極粒子は理想的には球状であるため、粒子内部の拡散は３次元拡散に従うと考える。そこで、３次元球対称の場合を考えると、拡散方程式は、

となる。ただし、ｒは球面座標上での位置を表わす。これに対する基本解は、

となる。

　これを一般化すると、３次元拡散方程式の基本解に基づく項は、パラメータｋと関数ｆを用いて、

と記述することができる。

　１次元拡散方程式のときと同様に、初期分布がｚ_p秒後の基本解に従うとすると、粒子表面上ｒ＝ｒ_dでの電位は、

となる。

　これを一般化すると、

と記述することができる。

　［蓄電池への適用］
　図１２は、正極と負極を有する蓄電池の内部の１次元等価回路である。ここまでの説明は、正極または負極の一方のみの挙動を説明したものであったが、図１２では、セパレータを介して、正極と負極の両方の１次元等価回路が接続されている。

　図１２の等価回路から電圧緩和モデルを構築すると、一例として、

となる。ここで、右肩の添字ｐとｎは、それぞれ正極と負極のパラメータであることを意味している。

　式（２４）は必ずしもそのまま用いる必要はなく、一部の項や一部のパラメータを省略してもよい。

　例えば、指数関数の部分は、時間ｔの経過とともに急速に減衰することから、Ｗ^e＝０、つまり指数関数部分を１としてもよい。ただし、右肩の添字ｅは、ｐまたはｎを意味する。

　さらに、正極と負極の項を１つにまとめてもよい。粒子内拡散については、一般的に時定数が数秒～数百秒のオーダであることから、条件に応じて以下の近似式が成り立つ。

　ここで、右辺の第１条件は、正負極の拡散時定数が近い値である、あるいは、電極の一方の時定数が他方の時定数に比べて非常に大きいということを意味している。また、右辺の第２条件は、推定に用いる時刻ｔがｚ₁とｚ₂よりも十分に大きいときを表わしている。

　また、式（２４）では明示していないが、第３条件として、時間ｔがさらに大きいとき、粒子内拡散そのものをゼロと近似することも可能である。

　一方、粒子間拡散については、粒子内拡散と比較してｋｉが小さいことから、以下のような近似式が成り立つ。

　ただし、ｉ≠ｊである。ここで、右辺の第１条件は、正負極の粒子間拡散の時定数が近い値である、あるいは、一方の電極の過電圧が他方よりも非常に大きいことを意味している。また、右辺の第２条件は、推定に用いる時刻ｔが一方の電極の粒子間拡散の時定数よりも十分に大きいことを意味している。

　このように、正負極の特性、フィッティングするデータ範囲等によって、モデルはさまざまに形に簡単化される。

　本実施の形態１のモデルの本質は、これまで考慮されてこなかった、あるいは実用に耐えうるモデルの構築が困難であった、粒子間拡散を含む電圧緩和過程を、拡散方程式の基本解に基づくモデルで表現し、かつ、その基本解モデルに対して、時間移動パラメータｚを導入した点にある。これをより一般的なモデル式で記述すると、

となる。

　ここで、関数ｆ（ｔ）は、比較的短時間の直流抵抗、反応抵抗、粒子内拡散抵抗等に起因する過電圧を表現する関数である。この関数には、一般的に用いられるような、直流抵抗素子と１つ以上のＣＲ並列素子とによる等価回路を想定した関数を用いてもよい。

　一方、関数ｇ（ｔ）は、粒子間拡散を表現する拡散基本解項を補正する関数であり、過去の充放電履歴に基づく電極内部の電位分布の局所的な高低などに由来するような、時間ｔが比較的小さいときの電位変動をよく説明するためのものである。

　式（２７）の関数ｆ（ｔ）の一例として、一般的な電池等価回路に基づく項、例えば、直流抵抗と多段直列接続されたＣＲ並列素子とに基づく項を用いてもよい。すなわち、

である。

　この場合、等価回路のパラメータとして、Ｒ、ｋ_i、τ_i（ｉ＝１，２，・・・ｎ_CR）の一部または全部を推定パラメータに含むことになる。

　関数ｇ（ｔ）は、その目的上、比較的短時間、具体的には数秒～数百秒程度で１に収束し、拡散基本解項への寄与がなくなるようなものが望ましい。例えば、

とする、あるいは、

とする等、さまざまな形が考えられる。

　なお、関数ｇ（ｔ）は、通常は上述の通り粒子間拡散項の補正のためのものであるが、関数ｆ（ｔ）の代わりの役目を果たすようなものであってもよい。

　あるいは、本実施の形態１のもう一つの重要な本質である、粒子内拡散を３次元拡散方程式の基本解で表現したモデルを明示すると、

となる。

　モデル提供部１０４は、以上のようにして構築される、拡散方程式の基本解に基づくモデルの少なくとも１つを提供する。拡散基本解モデルは、モデル提供部１０４の内部に記憶されていてもよいし、例えばクラウド上等のモデル提供部１０４の外部に記憶されていてもよい。

　モデル提供部１０４は、次に述べるモデルパラメータ推定部１０５からのモデル変更を指示する制御信号を受信する毎に、自身が提供可能な拡散基本解モデルの集合の中から、これまでに出力していないモデルを１つ選択して、モデルパラメータ推定部１０５に出力する。

　また、モデル提供部１０４は、全てのモデルを出力してしまい、これまでに出力していないモデルが存在しなくなった場合には、モデルパラメータ推定部１０５に対して、モデルパラメータ推定の終了を指示する制御信号を送信する。

　［モデルパラメータの推定］
　先述したように、モデルパラメータ推定部１０５は、時系列データ取得部１０３により取得された、蓄電池１０１の端子電圧Ｖの時系列データ｛（ｐ，Ｖｐ）｜ｐ＝０，１，・・・，ｍ－１｝と、モデル提供部１０４により提供された拡散基本解モデルとに基づいて、蓄電池１０１の拡散基本解モデルのモデルパラメータを推定する。

　図１３は、モデルパラメータ推定部１０５の内部構成を示す図である。モデルパラメータ推定部１０５は、最適化問題構成部１０５ａと、推定演算部１０５ｂと、判定部１０５ｃとを備えて構成されている。

　最適化問題構成部１０５ａは、解くべき最適化問題を構成する。詳細には、最適化問題構成部１０５ａは、時系列データ取得部１０３により取得された時系列データ｛（ｐ，Ｖｐ）｜０，１，・・・，ｍ－１｝と、モデル提供部１０４により提供された拡散基本解モデルと、推定演算部１０５ｂからの最適化問題の構成の変更を指示する制御信号とに基づいて、最適化問題を構成して、評価関数Ｌおよび制約不等式を出力する。

　評価関数Ｌは、最も典型的には、以下のように二乗誤差の総和によって構成される。

　ただし、θは推定すべきモデルパラメータを縦に並べたベクトルであり、例えば、

という拡散基本解モデルを用いる場合には、

となる。

　また、制約不等式は、

と表現される。例えば、式（３４）のｋ₁とｋ₂が非負であり、また、ＯＣＶがＶｍｉｎ以上かつＶｍａｘ以下であることが既知であるならば、

として、制約不等式を構成すればよい。

　以上のようにして、最適化問題

を構成する。ただし、制約不等式は必須のものではなく、制約なしの最適化問題を構成してもよい。

　最適化問題構成部１０５ａは、推定演算部１０５ｂからの最適化問題の構成の変更を指示する制御信号を受信した場合には、最適化問題の構成をこれまでと異なるものに変更する。具体的には、評価関数を変更する、制約不等式を追加する等の操作を行う。

　また、最適化問題構成部１０５ａは、推定演算部１０５ｂからの制御信号を受信した際に、これまでと異なる構成の候補が存在しない場合には、モデル提供部１０４に対して、モデルの変更を指示する制御信号を送信する。

　推定演算部１０５ｂは、最適化問題構成部１０５ａから出力される評価関数Ｌと制約不等式ｈ_i≦（ｉ＝１，・・・ｎ_h）とに基づいて、ある最適化手法を用いて推定演算を行なうことによって最適化問題の最適解を求め、推定したモデルパラメータθと評価関数Ｌの値を出力する。

　例えば、最適化問題が制約不等式を含まない場合には、最急降下法、Newton法、Gauss-Newton法、Quasi-Newton法、Levenberg-Marquardt法等を用いることができる。

　また、最適化問題が制約不等式を含む場合には、ペナルティ関数法、乗数法、逐次二次計画法、内点法、GRG（Generalized Reduced Gradient）法等を用いることができる。

　最適化手法は、問題の規模、変数の数、性質、計算資源等に応じて、適切なものを使い分ければよい。また、最適化手法は上述の方法に限定されるものではなく、公知の他の最適化手法を用いてもよい。

　あるいは、特定のモデルパラメータをある値に固定すれば他のモデルパラメータに関して線形となるような場合には、当該特定のモデルパラメータ値を固定したうえで、当該他のモデルパラメータに対して線形最小二乗法を適用して、当該他のモデルパラメータの値を求めた後、当該特定のモデルパラメータ値を変え、当該他のモデルパラメータに対して線形最小二乗法を適用して、当該他のモデルパラメータの値を再び求める、というような操作を繰り返すことで、評価関数を最小化するような特定のモデルパラメータと他のモデルパラメータの値の組合せを見つけてもよい。

　推定演算部１０５ｂは、判定部１０５ｃから推定演算のやり直しを指示する制御信号を受信した場合には、設定を変更して推定演算のやり直しを行なう。これは、初期値の変更、推定手法の変更等を意味する。ただし、やり直しの可能性がなくなった場合、やり直しの回数が設定値を超えた場合等には、最適化問題構成部１０５ａに対して、最適化問題の変更を指示する制御信号を送信する。

　判定部１０５ｃは、推定演算部１０５ｂから出力される評価関数Ｌの値と推定モデルパラメータθと、モデル提供部１０６からの制御信号とに基づいて、制御信号または推定モデルパラメータθを出力する。

　判定部１０７ｃは、推定演算部１０５ｂから出力される評価関数Ｌの値が設定値以下の場合には、モデルパラメータθの推定が成功したと判定し、推定されたモデルパラメータθを出力する。

　反対に、判定部１０５ｃは、評価関数Ｌの値が設定値より大きい場合には、推定演算部１０５ｂに対して、モデルパラメータ推定のやり直しを指示する制御信号を送信する。この際、推定演算部１０５ｂは、過去の評価関数Ｌの値と比較して、Ｌの値が小さい方の推定モデルパラメータθを保存する。

　判定部１０５ｃは、モデル提供部１０４からモデルパラメータ推定の終了を指示する制御信号を受信した場合には、モデルパラメータ推定部１０５の処理を終了させ、保存している推定モデルパラメータθを出力する。

　以上説明したように、本実施の形態１に係る蓄電池の特性推定装置は、蓄電池の端子電圧を検出する電圧検出部と、蓄電池の休止状態における端子電圧の時系列データを取得する時系列データ取得部と、蓄電池モデルを提供するモデル提供部と、時系列データ取得部により取得された端子電圧の時系列データと、モデル提供部により提供された蓄電池モデルとに基づいて、蓄電池モデルのモデルパラメータを推定する、モデルパラメータ推定部とを備えている。蓄電池モデルは、蓄電池のＯＣＶを表現するＯＣＶ項、および蓄電池の電極を構成する粒子間のイオン拡散を表現する１次元の拡散方程式の基本解に基づく粒子間拡散項を含んでいる。これにより、蓄電池の休止状態における電圧緩和過程を高精度に推定することができる、

　［モデルパラメータの推定結果］
　以下に、拡散基本解モデルとして式（２４）を用いて、実際の蓄電池の時系列データから、モデルパラメータを推定した結果を示す。なお、推定したモデルパラメータはｋ₁、ｋ₂、ｚ₂、ＯＣＶの４つであり、推定手法としてＧＲＧ法を用いた。

　図１４には、リチウムイオン蓄電池を解体して、正極材料を用いて作製した半電池の電位緩和曲線が示されている。ただし、電位変化の小さい方が新品の電位変化であり、電位変化の大きい方が劣化品の電位変化である。また、比較のため、７２００秒時点における両者の電位をゼロに合わせてプロットしている。

　明らかに、劣化品の電位変化が大きく、７２００秒時点でも電位が収束していないことが分かる。このように、劣化した蓄電池では、一般に電極電位の緩和に要する時間が非常に長くなる。そのため、高精度なＳＯＣ－ＯＣＶデータを取得するためには、電圧緩和過程を予測することが重要となる。

　図１５は、作製した半電池で取得した３種類のＳＯＣ－ＯＣＰデータが示されている。示す。３種類のデータは、それぞれ概ねＳＯＨ＝１００％、８０％、６０％に対応している。一般的に、ＳＯＣ－ＯＣＰ曲線は劣化しても形状が変化しないと言われている。しかしながら、図１５では、ＳＯＨ＝８０％と６０％のデータは、ＳＯＣ＝１００％の新品のデータと比較して、高ＳＯＣ領域でやや高電位となっている。この原因としては、電位緩和が十分に収束していないという可能性が考えられる。

　図１６には、図１５の全データに対して、本実施の形態１に係る手法を適用して、ＯＣＰを推定した結果が示されている。高ＳＯＣ領域における劣化品の高電位傾向がほぼ消失し、全区間に渡って新品時の形状とほぼ一致していることがわかる。つまり、真のＯＣＰを精度良く推定できている蓋然性が高い。

　［まとめ］
　本実施の形態１に係る電圧緩和過程のモデルの利点は、次のようにまとめられる。

　第１に、粒子間拡散という、従来考慮されてこなかった物理現象を反映したモデルを用いる点である。これにより、数時間～十時間以上かかる長時間の電圧緩和過程について、物理的根拠のある高精度なモデル化を可能とした。

　第２に、拡散方程式の基本解に基づく項を用いるというアイデアにより、少ないモデルパラメータ数と高い精度を両立するモデルの構築を可能とし、その理論的根拠も示した。

　第３に、上記の拡散方程式の基本解に基づく項を用いるうえで、通常の拡散方程式の基本解をそのまま用いるのではなく、時間移動パラメータを導入することで、高精度なモデル化を可能とし、その理論的根拠も示した。

　第４に、上述の粒子間拡散の項に加えて、一般的な電池モデルと同様に粒子内拡散の項を加えたモデルを構築する際、３次元拡散方程式の基本解に基づく項で表現することで少ないモデルパラメータ数で高精度なモデル化を可能とし、その理論的根拠を示した。電圧緩和過程には粒子内拡散と粒子間拡散が混在しているため、上述の第１から第３までの考えに基づく粒子間拡散項のみで電圧緩和を高精度にモデル化し、予測することができるとは限らない。つまり、粒子間拡散を高精度にフィッティングするためには、粒子内拡散を高精度にモデル化およびフィッティングすることもまた重要となる。それゆえに、この第４の考えについても、第１から第３の考えに基づく長時間の電圧緩和過程を表現する項のモデルパラメータを高精度にフィッティングする上で重要であり、第１から第３のアイデアの効果を強化するアイデアである。

　第５に、本実施の形態１に係る粒子間拡散モデルは、従来の直流抵抗とＣＲ並列素子などによる等価回路モデルのモデルパラメータをより高精度に推定することを可能とする。一般的な等価回路モデルは長くて数百秒程度までの時定数の過電圧をモデル化しているが、数百秒以上の時定数の電圧緩和過程である粒子間拡散を考慮していない分、原理的に精度が低いものとなっていた。これに対して、本実施の形態１では、長い時定数の粒子間拡散を高精度に推定・予測できるので、必然的な結果として、より短い過程、すなわち電解液抵抗、電荷移動抵抗、粒子内拡散抵抗等を反映する等価回路パラメータも、精度よく推定することが可能となる。

　以上より、従来技術と比較して、物理モデルに基づきながら、条件に応じてモデルを簡単化してモデルパラメータを削減することも可能であり、かつ、モデルパラメータを推定するので、事前にモデルパラメータ値を取得しておく必要もない。

　最後に、本実施の形態１の効果として、電圧緩和過程を高精度にモデル化可能となり、ＯＣＶを高精度に予測できるようになった結果、蓄電池のＳＯＣリセット、ＳＯＣ－ＯＣＶデータの取得等に必要な時間が短縮されることになる。あるいは、蓄電池のＳＯＣリセット、ＳＯＣ－ＯＣＶデータの取得を従来よりも高精度に行なうことが可能となる。これは、実用上、蓄電池搭載システムの高効率利用、稼働率向上、劣化蓄電池の迅速診断等の意味で産業応用上の利点となる。また、本実施の形態１により、粒子内拡散と粒子間拡散を高精度に分離推定可能となることから、GITT法（Galvanostatic Intermittent Titration Technique）等の休止を含む公知の測定技術において、高精度に蓄電池モデルを構築可能となる。

　なお、上述した本発明の実施の形態１に係る特性推定装置における各機能は、処理回路によって実現される。各機能を実現する処理回路は、専用のハードウェアであってもよく、メモリに格納されるプログラムを実行するプロセッサであってもよい。図１７は、本発明の実施の形態１に係る特性推定装置の各機能を専用のハードウェアである処理回路４０００で実現する場合を示した構成図である。また、図１８は、本発明の実施の形態１に特性推定装置の各機能をプロセッサ５００１およびメモリ５００２を備えた処理回路５０００により実現する場合を示した構成図である。

　処理回路が専用のハードウェアである場合、処理回路４０００は、例えば、単一回路、複合回路、プログラム化したプロセッサ、並列プログラム化したプロセッサ、ＡＳＩＣ（Ａｐｐｌｉｃａｔｉｏｎ　Ｓｐｅｃｉｆｉｃ　Ｉｎｔｅｇｒａｔｅｄ　Ｃｉｒｃｕｉｔ）、ＦＰＧＡ（Ｆｉｅｌｄ　Ｐｒｏｇｒａｍｍａｂｌｅ　Ｇａｔｅ　Ａｒｒａｙ）、またはこれらを組み合わせたものが該当する。特性推定装置の各部の機能それぞれを個別の処理回路４０００で実現してもよいし、各部の機能をまとめて処理回路４０００で実現してもよい。また、時系列データ、拡散基本解モデル等のデータは、処理回路４０００の内部に記憶されてもよいし、例えばクラウド上等の処理回路４０００の外部に記憶されてもよい。

　一方、処理回路がプロセッサ５００１の場合、特性推定装置の各部の機能は、ソフトウェア、ファームウェア、またはソフトウェアとファームウェアとの組み合わせにより実現される。ソフトウェアおよびファームウェアは、プログラムとして記述され、メモリ５００２に格納される。プロセッサ５００１は、メモリ５００２に記憶されたプログラムを読み出して実行することにより、各部の機能を実現する。すなわち、特性推定装置は、処理回路５０００により実行されるときに、上述した各制御が結果的に実行されることになるプログラムを格納するためのメモリ５００２を備える。

　これらのプログラムは、上述した各部の手順あるいは方法をコンピュータに実行させるものであるともいえる。ここで、メモリ５００２とは、例えば、ＲＡＭ（Ｒａｎｄｏｍ　Ａｃｃｅｓｓ　Ｍｅｍｏｒｙ）、ＲＯＭ（Ｒｅａｄ　Ｏｎｌｙ　Ｍｅｍｏｒｙ）、フラッシュメモリ、ＥＰＲＯＭ（Ｅｒａｓａｂｌｅ　Ｐｒｏｇｒａｍｍａｂｌｅ　Ｒｅａｄ　Ｏｎｌｙ　Ｍｅｍｏｒｙ）、ＥＥＰＲＯＭ（Ｅｌｅｃｔｒｉｃａｌｌｙ　Ｅｒａｓａｂｌｅ　ａｎｄ　Ｐｒｏｇｒａｍｍａｂｌｅ　Ｒｅａｄ　Ｏｎｌｙ　Ｍｅｍｏｒｙ）等の、不揮発性または揮発性の半導体メモリが該当する。また、磁気ディスク、フレキシブルディスク、光ディスク、コンパクトディスク、ミニディスク、ＤＶＤ等も、メモリ２００２に該当する。また、時系列データ、拡散基本解モデル等のデータは、メモリ５００２に記憶されてもよいし、例えばクラウド上等に記憶されてもよい。

　なお、上述した各部の機能について、一部を専用のハードウェアで実現し、一部をソフトウェアまたはファームウェアで実現するようにしてもよい。

　このように、処理回路は、ハードウェア、ソフトウェア、ファームウェア、またはこれらの組み合わせによって、上述した各部の機能を実現することができる。

　１００　特性推定装置、１０１　蓄電池、１０２　電圧検出部、１０３　時系列データ取得部、１０４　モデル提供部、１０５　モデルパラメータ推定部。

Claims

　蓄電池の端子電圧を検出する電圧検出部と、
　前記蓄電池の休止状態における前記端子電圧の時系列データを取得する時系列データ取得部と、
　蓄電池モデルを提供するモデル提供部であって、該蓄電池モデルは、前記蓄電池のＯＣＶを表現するＯＣＶ項、および前記蓄電池の電極を構成する粒子間のイオン拡散を表現する１次元の拡散方程式の基本解に基づく粒子間拡散項を含む、モデル提供部と、
　前記時系列データ取得部により取得された前記端子電圧の時系列データと、前記モデル提供部により提供された前記蓄電池モデルとに基づいて、前記蓄電池モデルのモデルパラメータを推定する、モデルパラメータ推定部と
を備える、蓄電池の特性推定装置。
　前記粒子間拡散項は、時間移動パラメータをさらに含む、請求項１に記載の特性推定装置。
　前記蓄電池モデルは、前記蓄電池の電極を構成する粒子内のイオン拡散を表現する、３次元の拡散方程式の基本解に基づく粒子内拡散項をさらに含む、請求項１または２に記載の特性推定装置。
　前記蓄電池モデルは、電池等価回路に基づく項をさらに含む、請求項１～３のいずれか一項に記載の特性推定装置。
　蓄電池の端子電圧を検出するステップと、
　前記蓄電池の休止状態における前記端子電圧の時系列データを取得するステップと、
　蓄電池モデルを提供するステップであって、該蓄電池モデルは、前記蓄電池のＯＣＶを表現するＯＣＶ項、および前記蓄電池の電極を構成する粒子間のイオン拡散を表現する１次元の拡散方程式の基本解に基づく粒子間拡散項を含む、蓄電池モデルを提供するステップと、
　前記端子電圧の時系列データと、前記蓄電池モデルとに基づいて、前記蓄電池モデルのモデルパラメータを推定するステップと
を含む、蓄電池の特性推定方法。
　前記粒子間拡散項は、時間移動パラメータをさらに含む、請求項５に記載の特性推定方法。
　前記蓄電池モデルは、前記蓄電池の電極を構成する粒子内のイオン拡散を表現する、３次元の拡散方程式の基本解に基づく粒子内拡散項をさらに含む、請求項５または６に記載の特性推定方法。
　前記蓄電池モデルは、電池等価回路に基づく項をさらに含む、請求項５～７のいずれか一項に記載の特性推定方法。