(57)【特許請求の範囲】
【請求項１】
蓄電池の端子間電圧を計測する電圧計測部と、
予め測定した、容量誤差を含む前記蓄電池の定格満充電容量および前記端子間電圧に依存するように前記蓄電池のＳｏＣ推定値を推定するＳｏＣ推定部と、
前記蓄電池を流れる、オフセット誤差を含む電流を計測する電流計測部と、
前記ＳｏＣ推定値の、電流測定値の前記オフセット誤差に関する偏微分係数である第１係数を計算する第１係数計算部と、
前記ＳｏＣ推定値の、予め測定された前記蓄電池の前記定格満充電容量と前記蓄電池の使用時における実際の満充電容量との差である容量誤差に関する偏微分係数である第２係数を計算する第２係数計算部と、
前記第１係数および前記第２係数を含む係数情報と、前記蓄電池を流れる電流と、前記ＳｏＣ推定値とから、電流測定値の前記オフセット誤差および前記満充電容量の容量誤差を推定する誤差推定部と、
を備える
蓄電池パラメータ推定装置。
【請求項２】
前記電流計測部によって計測された前記電流を積算して電流積算値を計算する電流積算部を更に備え、
前記誤差推定部は、前記係数情報と、前記電流積算値と、前記ＳｏＣ推定値とから、前記オフセット誤差および容量誤差を推定する
請求項１に記載の蓄電池パラメータ推定装置。
【請求項３】
前記誤差推定部は、前記係数情報の今回の計算値と前回の計算値との差分と、前記ＳｏＣ推定値の今回の計算値と前回の計算値との差分と、記蓄電池を流れる電流とから、前記オフセット誤差および容量誤差を推定する
請求項１に記載の蓄電池パラメータ推定装置。
【請求項４】
前記ＳｏＣ推定部は、前記定格満充電容量と、前記端子間電圧と、予め測定した前記蓄電池の等価回路の特性パラメータとから、前記蓄電池の前記ＳｏＣ推定値を推定し、
前記ＳｏＣ推定値の、前記蓄電池の前記等価回路の特性パラメータの定格値と実際値との差に関する偏微分係数である第３係数を計算する第３係数計算部を更に備え、
前記誤差推定部が参照する前記係数情報は、前記第１係数、前記第２係数、および前記第３係数を含む
請求項１から３のいずれか１項に記載の蓄電池パラメータ推定装置。
【請求項５】
前記ＳｏＣ推定部は、Ｋａｌｍａｎフィルタを用いて、前記蓄電池の前記ＳｏＣ推定値を推定する、
請求項１から４のいずれか１項に記載の蓄電池パラメータ推定装置。
【請求項６】
前記誤差推定部は、
前記係数情報またはその今回の計算値と前回の計算値との差分と、前記ＳｏＣ推定値またはその今回の計算値と前回の計算値との差分とから、説明変数を計算する説明変数計算部と、
前記係数情報またはその今回の計算値と前回の計算値との差分と、前記ＳｏＣ推定値またはその今回の計算値と前回の計算値との差分とから、観測変数を計算する観測変数計算部と、
前記説明変数および前記観測変数に基づいて、前記オフセット誤差および前記容量誤差に関するパラメータの推定値であるパラメータ推定ベクトルを計算するパラメータ推定部と、
前記パラメータ推定ベクトルに基づいて、前記オフセット誤差を補正するオフセット計算部と、
前記パラメータ推定ベクトルに基づいて、容量誤差を推定する容量計算部と、
を備える
請求項１から５のいずれか１項に記載の蓄電池パラメータ推定装置。
【請求項７】
前記誤差推定部は、
前記説明変数と、前周期の逆分散とから、ＲＬＳゲインを計算するＲＬＳゲイン計算部と、
前周期の前記逆分散と、前記ＲＬＳゲインと、前記説明変数とから、前記逆分散を計算する逆分散計算部と、
を更に備え、
前記パラメータ推定部は、前記説明変数と、前記観測変数と、前記ＲＬＳゲインと、前周期のパラメータ推定ベクトルとから、前記パラメータ推定ベクトルを計算する
請求項６に記載の蓄電池パラメータ推定装置。
【請求項８】
前記誤差推定部は、
前記説明変数と、前記観測変数と、前周期の共分散行列とから、前記共分散行列を計算する共分散計算部と、
前記共分散行列に関する固有値問題を解くことにより、前記共分散行列の最小固有値に対応する固有ベクトルを計算する固有ベクトル計算部と、
を更に備え、
前記パラメータ推定部は、前記固有ベクトルから、前記パラメータ推定ベクトルを計算する
請求項６に記載の蓄電池パラメータ推定装置。
【請求項９】
前記誤差推定部は、
前記説明変数と、前記観測変数と、前周期の共分散行列とから、前記共分散行列を計算する共分散計算部と、
前記説明変数と、前記共分散行列と、前周期のパラメータ推定ベクトルとから、前記パラメータ推定ベクトルを求めるための２次式の係数を計算する係数計算部と、
前記２次式の係数を持つ２次式の実数根を求める２次式求解部と、
を更に備え、
前記パラメータ推定部は、前記２次式の係数と、前記２次式の実数根とから、前記パラメータ推定ベクトルを計算する
請求項６に記載の蓄電池パラメータ推定装置。
【請求項１０】
蓄電池の端子間電圧を計測する電圧計測ステップと、
予め測定した、容量誤差を含む前記蓄電池の定格満充電容量および前記端子間電圧に依存するように前記蓄電池のＳｏＣ推定値を推定するＳｏＣ推定ステップと、
前記蓄電池を流れる、オフセット誤差を含む電流を計測する電流計測ステップと、
前記ＳｏＣ推定値の、電流測定値の前記オフセット誤差に関する偏微分係数である第１係数を計算する第１係数計算ステップと、
前記ＳｏＣ推定値の、予め測定された前記蓄電池の前記定格満充電容量と前記蓄電池の使用時における実際の満充電容量との差である容量誤差に関する偏微分係数である第２係数を計算する第２係数計算ステップと、
前記第１係数および前記第２係数を含む係数情報と、前記蓄電池を流れる電流と、前記ＳｏＣ推定値とから、電流測定値の前記オフセット誤差および前記満充電容量の容量誤差を推定する誤差推定ステップと、
を有する
蓄電池パラメータ推定方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、蓄電池を流れる電流測定値のオフセット誤差、および予め測定された蓄電池の定格満充電容量（ＦＣＣ：Ｆｕｌｌ  Ｃｈａｒｇｅ  Ｃａｐａｃｉｔｙ）と蓄電池の使用時における実際のＦＣＣとの差である容量誤差を推定する蓄電池パラメータ推定装置および蓄電池パラメータ推定方法に関するものである。
【背景技術】
【０００２】
リチウムイオン蓄電池やニッケル水素電池、鉛蓄電池等の蓄電池を効率的に使用するためには、蓄電池の充電状態（ＳｏＣ：  Ｓｔａｔｅ  ｏｆ  Ｃｈａｒｇｅ）を高精度に推定することが不可欠である。しかし、蓄電池のＳｏＣは直接計測できる物理量ではないため、蓄電池を流れる電流や蓄電池の端子間電圧、蓄電池温度の測定値、更にこれらの物理量の履歴からＳｏＣを推定する必要がある。蓄電池のＳｏＣを推定する最も基本的な手法としては、開回路電圧法および電流積算法が知られている。
【０００３】
開回路電圧法とは、蓄電池を静定状態に置き、蓄電池内部の熱力学的状態が十分に平衡に達した蓄電池の端子間電圧、すなわち開回路電圧（ＯＣＶ：  Ｏｐｅｎ  Ｃｉｒｃｕｉｔ  Ｖｏｌｔａｇｅ）を計測し、事前に決定されるＳｏＣに対するＯＣＶの依存性に基づいて、蓄電池のＳｏＣの値を求めるものである。このため、開回路電圧法によるＳｏＣの推定中は、長時間にわたって蓄電池を使用することができない。
【０００４】
一方、電流積算法とは、蓄電池を流れる電流を積算し、蓄電池のＦＣＣで除算することで、積算開始時刻からのＳｏＣの変化を計算するものである。電流積算法では、ＳｏＣが明らかな状態、すなわち蓄電池が満充電状態か完全放電状態に有る場合を基準にとることで、蓄電池を使用しながらＳｏＣを推定することができる。しかし、電流積算法では、広く用いられているホール効果素子では避けることのできない電流値のわずかなオフセット誤差（以下、「オフセット誤差」と略す）によって、蓄電池のＳｏＣの推定精度が大きく悪化してしまうという問題がある。
【０００５】
電流積算法を用いた従来の装置として、蓄電池の等価回路を基に、蓄電池を流れる電流を入力とし、蓄電池の端子間電圧を出力とする線形システムを考え、この線形システムの状態に、電流測定値のオフセット誤差を追加した拡大系をもとに、Ｋａｌｍａｎフィルタを構成することで、蓄電池のＳｏＣと電流測定値のオフセット誤差を精度よく推定するバッテリの充電率推定装置が開示されている（例えば、特許文献１参照）。
【０００６】
また、開回路電圧法と電流積算法を組合せた従来の装置として、蓄電池を流れる電流が一定以上の絶対値を持つような、オフセット誤差の影響が少ない期間に限定して蓄電池を流れる電流を積算し、同期間におけるＳｏＣの変化を比較することで、蓄電池のＦＣＣを推定する電池容量算出装置が開示されている（例えば、特許文献２参照）。
【０００７】
特許文献２では、ＳｏＣはＯＣＶから推定され、ＯＣＶは蓄電池の等価回路に基づく再帰的演算によって求められる。等価回路の特性パラメータは、蓄電池を流れる電流と蓄電池の端子間電圧の推移から適応機構によって調整される。この結果、電流測定値のオフセット誤差は、電流積算期間を電流が一定以上の絶対値を持つときに限ることで、その影響を減じられるとしている。
【先行技術文献】
【特許文献】
【０００８】
【特許文献１】特開２０１２−５７９６４号公報
【特許文献２】特開２０１２−５８０２８号公報
【非特許文献】
【０００９】
【非特許文献１】Ｄ．−Ｚ．Ｆｅｎｇ，Ｘ．−Ｄ．Ｚｈａｎｇ，Ｄ．−Ｘ．Ｃｈａｎｇ  ａｎｄ  Ｗ．Ｘ．Ｚｈｅｎｇ，“Ａ  ｆａｓｔ  ｒｅｃｕｒｓｉｖｅ  ｔｏｔａｌ  ｌｅａｓｔ  ｓｑｕａｒｅｓ  ａｌｇｏｒｉｔｈｍ  ｆｏｒ  ａｄａｐｔｉｖｅ  ＦＩＲ  ｆｉｌｔｅｒｉｎｇ，”ＩＥＥＥ  Ｔｒａｎｓａｃｔｉｏｎ  ｏｎ  Ｓｉｇｎａｌ  Ｐｒｏｃｅｓｓｉｎｇ，ｖｏｌ．５２，Ｎｏ．１０，２００４．
【発明の概要】
【発明が解決しようとする課題】
【００１０】
しかしながら、従来技術には、以下のような課題がある。
特許文献１では、蓄電池の定格ＦＣＣおよび等価回路の特性パラメータを既知としているために、蓄電池の定格ＦＣＣおよび等価回路の特性パラメータに誤差が含まれる場合には、これらの誤差がＳｏＣや電流測定値のオフセット誤差の推定値に伝播してしまうという課題がある。
【００１１】
一方、特許文献２では、等価回路の特性パラメータを適合させる適応機構が参照する電流値には、オフセット誤差が含まれている。従って、等価回路の特性パラメータは、オフセット誤差に起因する誤差を持つ。この誤差は、等価回路の特性パラメータに依存して算出されるＯＣＶの推定値、さらにはＳｏＣの推定値に伝播し、最終的には蓄電池のＦＣＣの推定値まで伝播する。
【００１２】
更に、例えば、特許文献２に開示された技術を用いて蓄電池のＦＣＣを推定し、この推定したＦＣＣを用いて、特許文献１に開示された手法により、電流測定値のオフセット誤差を推定したとしても、このオフセット誤差には、実際の電流測定値のオフセット誤差の推定値に起因する誤差が重畳していることとなり、正しい電流オフセット値を求めることができない。
【００１３】
このように、従来の装置では、オフセット誤差に起因する誤差が蓄電池のＦＣＣの推定値まで伝播するために、蓄電池を流れる電流の測定値がオフセット誤差を有している限り、蓄電池のＦＣＣを正しく推定することができないという課題があった。
【００１４】
本発明は、上記のような課題を解決するためになされたものであり、電流測定値がオフセット誤差を有する場合でも、蓄電池を流れる電流測定値のオフセット誤差、および予め測定された蓄電池の定格満充電容量（ＦＣＣ）と蓄電池の使用時における実際のＦＣＣとの差である容量誤差を高精度に推定することができる蓄電池パラメータ推定装置および蓄電池パラメータ推定方法を得ることを目的とする。
【課題を解決するための手段】
【００１５】
本発明に係る蓄電池パラメータ推定装置は、蓄電池の端子間電圧を計測する電圧計測部と、予め測定した、容量誤差を含む蓄電池の定格ＦＣＣおよび端子間電圧に依存するように蓄電池のＳｏＣ推定値を推定するＳｏＣ推定部と、蓄電池を流れる、オフセット誤差を含む電流を計測する電流計測部と、ＳｏＣ推定値の、電流測定値のオフセット誤差に関する偏微分係数である第１係数を計算する第１係数計算部と、ＳｏＣ推定値の、予め測定された蓄電池の定格ＦＣＣと蓄電池の使用時における実際のＦＣＣとの差である容量誤差に関する偏微分係数である第２係数を計算する第２係数計算部と、第１係数および第２係数を含む係数情報と、蓄電池を流れる電流と、ＳｏＣ推定値とから、電流測定値のオフセット誤差およびＦＣＣの容量誤差を推定する誤差推定部とを備えるものである。
【００１６】
本発明に係る蓄電池パラメータ推定方法は、蓄電池の端子間電圧を計測する電圧計測ステップと、予め測定した、容量誤差を含む蓄電池の定格ＦＣＣおよび端子間電圧に依存するように蓄電池のＳｏＣ推定値を推定するＳｏＣ推定ステップと、蓄電池を流れる、オフセット誤差を含む電流を計測する電流計測ステップと、ＳｏＣ推定値の、電流測定値のオフセット誤差に関する偏微分係数である第１係数を計算する第１係数計算ステップと、ＳｏＣ推定値の、予め測定された蓄電池の定格ＦＣＣと蓄電池の使用時における実際のＦＣＣとの差である容量誤差に関する偏微分係数である第２係数を計算する第２係数計算ステップと、第１係数および第２係数を含む係数情報と、蓄電池を流れる電流と、ＳｏＣ推定値とから、電流測定値のオフセット誤差およびＦＣＣの容量誤差を推定する誤差推定ステップとを有するものである。
【発明の効果】
【００１７】
本発明では、電流測定値および定格ＦＣＣの何れかまたは両方に誤差が含まれ、この誤差は、蓄電池の端子間電圧の測定値から推定されるＳｏＣ推定値に基づいて推定される。この結果、電流測定値に誤差が含まれる場合でも、電流測定値のオフセット誤差およびＦＣＣの容量誤差を高精度に推定することができる蓄電池パラメータ推定装置および蓄電池パラメータ推定方法を得ることができる。
【図面の簡単な説明】
【００１８】
【図１】本発明の実施の形態１に係る蓄電池パラメータ推定装置の全体構成図である。
【図２】本発明の実施の形態１に係る蓄電池の等価回路である。
【図３】本発明の実施の形態１に係るＳｏＣ推定部の全体構成図である。
【図４】本発明の実施の形態１に係る第１係数計算部の全体構成図である。
【図５】本発明の実施の形態１に係る第２係数計算部の全体構成図である。
【図６】本発明の実施の形態１に係る誤差推定部の全体構成図である。
【図７】本発明の実施の形態１に係る説明変数計算部のブロック線図である。
【図８】本発明の実施の形態１に係る観測変数計算部ブロック線図である。
【図９】本発明の実施の形態１に係る蓄電池パラメータ推定方法のフローチャートである。
【図１０】本発明の実施の形態１に係るＳｏＣ推定部が行う処理の１周期のフローチャートである。
【図１１】本発明の実施の形態１に係る第１係数計算部が行う処理の１周期のフローチャートである。
【図１２】本発明の実施の形態１に係る第２係数計算部が行う処理の１周期のフローチャートである。
【図１３】本発明の実施の形態１に係る誤差推定部が行う処理の１周期のフローチャートである。
【図１４】本発明の実施の形態２に係る蓄電池パラメータ推定装置の全体構成図である。
【図１５】本発明の実施の形態２に係る蓄電池パラメータ推定方法のフローチャートである。
【図１６】本発明の実施の形態３に係る蓄電池パラメータ推定装置の全体構成図である。
【図１７】本発明の実施の形態３に係る誤差推定部の全体構成図である。
【図１８】本発明の実施の形態３に係る説明変数計算部のブロック線図である。
【図１９】本発明の実施の形態３に係る観測変数計算部のブロック線図である。
【図２０】本発明の実施の形態３に係る蓄電池パラメータ推定方法のフローチャートである。
【図２１】本発明の実施の形態３に係る誤差推定部が行う処理の１周期のフローチャートである。
【図２２】本発明の実施の形態４に係る蓄電池パラメータ推定装置の全体構成図である。
【図２３】本発明の実施の形態４に係る誤差推定部の全体構成図である。
【図２４】本発明の実施の形態４に係る蓄電池パラメータ推定方法のフローチャートである。
【図２５】本発明の実施の形態４に係る誤差推定部が行う処理の１周期のフローチャートである。
【図２６】本発明の実施の形態５に係る蓄電池パラメータ推定装置の全体構成図である。
【図２７】本発明の実施の形態５に係る誤差推定部の全体構成図である。
【図２８】本発明の実施の形態５に係る蓄電池パラメータ推定方法のフローチャートである。
【図２９】本発明の実施の形態５に係る誤差推定部が行う処理の１周期のフローチャートである。
【発明を実施するための形態】
【００１９】
以下、本発明における蓄電池パラメータ推定装置および蓄電池パラメータ推定方法の好適な実施の形態について図面を用いて説明する。なお、各図において同一、または相当する部分については、同一符号を付して説明する。
【００２０】
実施の形態１．
図１は、本発明の実施の形態１に係る蓄電池パラメータ推定装置１０１の全体構成図である。電流計測部１０２は、蓄電池１０３を流れる電流１０４を計測する。また、電圧計測部１０５は、蓄電池１０３の正極および負極の端子間の電圧１０６を計測する。ここで、電流計測部１０２が測定する電流１０４の測定値には、オフセット誤差１１６が含まれている。
【００２１】
ＳｏＣ推定部１０７は、予め測定した蓄電池１０３の定格満充電容量（ＦＣＣ）および電圧計測部１０５によって計測された電圧１０６に依存した方法で、蓄電池１０３の等価回路に基づいて、蓄電池１０３のＳｏＣ推定値１０８を推定する。なお、図１に示すように、ＳｏＣ推定部１０７は、電流計測部１０２によって計測された電流１０４を更に参照するようにしてもよい。ここで、蓄電池１０３の定格ＦＣＣは、予め測定された蓄電池１０３の定格ＦＣＣと蓄電池１０３の使用時における実際のＦＣＣとの差である容量誤差１１７を有している。
【００２２】
第１係数計算部１０９は、ＳｏＣ推定部１０７の内部演算を参照しつつ、ＳｏＣ推定値１０８の、電流計測部１０２によって測定された電流１０４のオフセット誤差１１６に関する偏微分係数である第１係数１１０を、蓄電池１０３の等価回路に基づいて計算する。
【００２３】
第２係数計算部１１１は、ＳｏＣ推定部１０７の内部演算を参照しつつ、ＳｏＣ推定値１０８の、予め測定された蓄電池１０３の定格ＦＣＣと実際のＦＣＣとの容量誤差１１７に関する偏微分係数である第２係数１１２を、蓄電池１０３の等価回路に基づいて計算する。
【００２４】
電流積算部１１３は、電流計測部１０２によって計測された電流１０４を積算することで、電流積算値１１４を計算する。
【００２５】
誤差推定部１１５は、第１係数計算部１０９によって計算された第１係数１１０と、第２係数計算部１１１によって計算された第２係数１１２と、電流積算部１１３によって計算された電流積算値１１４と、ＳｏＣ推定部１０７によって推定されたＳｏＣ推定値１０８とから、オフセット誤差１１６および容量誤差１１７を推定して出力する。
【００２６】
図２は、本発明の実施の形態１に係る蓄電池１０３の等価回路である。図２に示す電圧Ｅ［Ｖ］を有する電圧源２０１は、蓄電池１０３のＯＣＶ（開回路電圧）を表す。ＯＣＶは、ＳｏＣの関数であり、蓄電池１０３毎に定められる。抵抗Ｒｄ［Ω］の電気抵抗２０２と、容量Ｃｄ［Ｆ］のキャパシタンス２０３とで構成される並列回路２０４は、蓄電池１０３の緩和過程を表す。抵抗Ｒ０［Ω］の電気抵抗２０５は、緩和過程以外の電池内部抵抗を表す。
【００２７】
ここで、電気量ｑｂ［Ｃ］を蓄電池１０３に蓄えられている電気量とし、蓄電池１０３の定格ＦＣＣをＦｃｃ［Ｃ］とし、ｓ［％］を蓄電池１０３のＳｏＣとすると、電気量ｑｂ［Ｃ］は、
【数１】
0006388563.tif 000002
で表される。
【００２８】
等価回路に関するＫｉｒｃｈｈｏｆｆの電流法則、および電気抵抗に関するＯｈｍの法則によれば、Ｖ［Ｖ］で表される蓄電池１０３の電圧１０６は、以下の微分方程式系、
【数２】
0006388563.tif 000003
に従う。ただし、
【数３】
0006388563.tif 000004
である。また、Ｉは蓄電池１０３を流れる電流１０４、Ｃｄ、Ｒｄ、Ｒ０は、図２に示す蓄電池１０３の等価回路の特性パラメータである。
【００２９】
今、電流計測部１０２および電圧計測部１０５のサンプリング周期を、ｔｓ［秒］とすると、微分方程式系（２）は、以下の差分方程式系に離散化される。
【数４】
0006388563.tif 000005
ただし、
【数５】
0006388563.tif 000006
である。
【００３０】
図３は、本発明の実施の形態１に係るＳｏＣ推定部１０７の全体構成図である。本実施の形態１におけるＳｏＣ推定部１０７は、図２の等価回路、すなわち差分方程式系（５）に基づき、拡張Ｋａｌｍａｎフィルタを用いて、蓄電池１０３のＳｏＣ推定値１０８を計算する。
【００３１】
ＳｏＣ推定部１０７は、電流計測部１０２、電圧計測部１０５による計測ごとに周期的に動作する。ある周期における動作において、状態予測部３０１は、電流計測部１０２によって計測されて電流記憶部３０２に記憶された１周期前の電流１０４である電流Ｉｋ−１と、状態推定部３１０で計算されて状態記憶部３０３に記憶された１周期前の推定状態ｘｋ−１とから、
【数６】
0006388563.tif 000007
に従って予測状態ｘ＾ｋを予測する。
【００３２】
測定残差計算部３０４は、状態予測部３０１によって計算された予測状態ｘ＾ｋと、電流計測部１０２によって周期的に計測された蓄電池１０３の電流１０４である電流Ｉｋと、電圧計測部１０５によって計測された蓄電池１０３の電圧１０６である電圧Ｖｋとから、測定残差ｅｋを、
【数７】
0006388563.tif 000008
に従って計算する。
【００３３】
共分散予測部３０５は、共分散推定部３１１によって計算され、共分散記憶部３０６に記憶された１周期前の共分散推定値Ｐｋ−１から、予測共分散Ｐ＾ｋを、
【数８】
0006388563.tif 000009
に従って計算する。ただし、Ｑとは、
【数９】
0006388563.tif 000010
によって予め計算され、Φ１２およびΦ２２とは、
【数１０】
0006388563.tif 000011
によって定義される行列である。また、σＩ２は、電流１０４の観測誤差の分散を表す所定の正の実数定数である。
【００３４】
観測行列計算部３０７は、状態予測部３０１によって計算された予測状態ｘ＾ｋから、
【数１１】
0006388563.tif 000012
に従って観測行列Ｈｋを計算する。
【００３５】
残差共分散計算部３０８は、観測行列計算部３０７によって計算された観測行列Ｈｋと、共分散予測部３０５によって計算された予測共分散Ｐ＾ｋとから、残差共分散Ｓｋを、
【数１２】
0006388563.tif 000013
に従って計算する。ただし、σＶ２は、電圧１０６の観測誤差の分散を表す所定の正の実数定数である。
【００３６】
Ｋａｌｍａｎゲイン計算部３０９は、共分散予測部３０５によって計算された予測共分散Ｐ＾ｋと、観測行列計算部３０７によって計算された観測行列Ｈｋと、残差共分散計算部３０８によって計算された残差共分散Ｓｋとから、ＫａｌｍａｎゲインＫｋを、
【数１３】
0006388563.tif 000014
によって計算する。
【００３７】
状態推定部３１０は、状態予測部３０１によって計算された予測状態ｘ＾ｋと、Ｋａｌｍａｎゲイン計算部３０９によって計算されたＫａｌｍａｎゲインＫｋと、測定残差計算部３０４によって計算された測定残差ｅｋとから、推定状態ｘｋを、
【数１４】
0006388563.tif 000015
によって計算し、結果を状態記憶部３０３に記憶させる。
【００３８】
共分散推定部３１１は、Ｋａｌｍａｎゲイン計算部３０９によって計算されたＫａｌｍａｎゲインＫｋと、観測行列計算部３０７によって計算された観測行列Ｈｋと、共分散予測部３０５によって計算された予測共分散Ｐ＾ｋとから、共分散推定値Ｐｋを、
【数１５】
0006388563.tif 000016
によって計算し、共分散記憶部３０６に記憶させる。
【００３９】
ＳｏＣ計算部３１２は、状態推定部３１０によって計算された推定状態ｘｋの第２要素である充電量ｑｂ，ｋを、予め測定した蓄電池１０３の定格ＦＣＣＦｃｃで除することによって、ＳｏＣ推定値１０８を計算する。
【００４０】
ＳｏＣ推定部１０７は、ＳｏＣ推定値１０８を算出するために電流１０４の測定値を参照する。もし、電流１０４の測定値が実際に蓄電池１０３を流れている電流と厳密に等しく、また、蓄電池１０３の特性が、図２に示した等価回路と厳密に一致するなら、ＳｏＣ推定部１０７は、平均二乗誤差の観点から最良のＳｏＣ推定値１０８を与える。
【００４１】
しかしながら、電流１０４の測定値が、オフセット誤差ΔＩを含むならば、ＳｏＣ推定値１０８は、電流測定値のオフセット誤差に起因する偏り誤差を持つ。本実施の形態１における第１係数計算部１０９は、このオフセット誤差ΔＩに起因するＳｏＣ推定値１０８の偏り誤差を近似的に定量する。
【００４２】
図４は、本発明の実施の形態１に係る第１係数計算部１０９の全体構成図である。予測状態微分係数計算部４０１は、１周期前に推定状態微分係数計算部４０４で計算され、推定状態微分係数記憶部４０２に記憶された、パラメータｐ１に関する推定状態の偏微分係数ｄｘｋ−１／ｄｐ１から、パラメータｐ１に関する予測状態の偏微分係数ｄｘ＾ｋ／ｄｐ１を、
【数１６】
0006388563.tif 000017
に従い計算する。ここで、オフセット誤差ΔＩは、ΔＩ＝Ｉｔｙｐ×ｐ１と表され、Ｉｔｙｐは、パラメータｐ１の絶対値を１よりも十分に小さくするための定数である。
【００４３】
測定残差微分係数計算部４０３は、予測状態微分係数計算部４０１によって計算された、パラメータｐ１に関する予測状態の偏微分係数から、パラメータｐ１に関する測定残差の偏微分係数ｄｅｋ／ｄｐ１を、
【数１７】
0006388563.tif 000018
に従い計算する。ただし、観測行列Ｈｋは、状態予測部３０１によって計算される予測状態ｘ＾ｋから、式（１３）によって計算される。
【００４４】
推定状態微分係数計算部４０４は、予測状態微分係数計算部４０１によって計算された、パラメータｐ１に関する予測状態の偏微分係数ｄｘ＾ｋ／ｄｐ１と、Ｋａｌｍａｎゲイン計算部３０９によって計算されたＫａｌｍａｎゲインＫｋと、測定残差微分係数計算部４０３によって計算された、パラメータｐ１に関する測定残差の偏微分係数ｄｅｋ／ｄｐ１とから、パラメータｐ１に関する推定状態の偏微分係数ｄｘｋ／ｄｐ１を、
【数１８】
0006388563.tif 000019
に従い計算する。
【００４５】
ＳｏＣ微分係数計算部４０５は、推定状態微分係数計算部４０４によって計算された推定状態の偏微分係数ｄｘｋ／ｄｐ１の第２要素ｄｑｂ，ｋ／ｄｐ１を、予め測定した蓄電池１０３の定格ＦＣＣＦｃｃで除することによって第１係数１１０を計算する。
【００４６】
ＳｏＣ推定部１０７では、ＳｏＣ推定値１０８を算出するために蓄電池１０３の定格ＦＣＣＦｃｃを用いている。もし、ＳｏＣ推定部１０７で用いた定格ＦＣＣＦｃｃが、蓄電池１０３の実際のＦＣＣと一致しており、電流測定値がオフセット誤差を持たず、また、蓄電池１０３の他の特性が図２に示した等価回路と一致するなら、図３に示したＳｏＣ推定部１０７は、二乗誤差の意味で最良のＳｏＣ推定値１０８を与える。
【００４７】
しかしながら、ＳｏＣ推定部１０７で用いた定格ＦＣＣＦｃｃが、容量誤差を持っている場合、ＳｏＣ推定値１０８は偏り誤差を持つ。第２係数計算部１１１は、この定格ＦＣＣＦｃｃの容量誤差に起因するＳｏＣ推定値１０８の偏り誤差を、近似的に定量する。
【００４８】
図５は、本発明の実施の形態１に係る第２係数計算部１１１の全体構成図である。以下、蓄電池１０３の実際のＦＣＣをＦｃｃとし、ＳｏＣ推定部１０７で用いた定格ＦＣＣをＦ〜ｃｃとする。ここで、Ｆ〜ｃｃがＦｃｃに対し、
【数１９】
0006388563.tif 000020
のようにずれていたとする。ここで、ずれはＦｃｃに対して十分小さいとし、パラメータｐ２の絶対値は１より十分小さいとする。
【００４９】
予測状態微分係数計算部５０１は、１周期前に推定状態微分係数計算部５０９で計算され、推定状態微分係数記憶部５０２に記憶された、パラメータｐ２に関する推定状態の偏微分係数ｄｘｋ−１／ｄｐ２から、パラメータｐ２に関する予測状態の偏微分係数ｄｘ＾ｋ／ｄｐ２を、
【数２０】
0006388563.tif 000021
によって計算する。
【００５０】
観測行列微分係数計算部５０３は、状態予測部３０１によって計算された予測状態ｘ＾ｋから、パラメータｐ２に関する観測行列の偏微分係数ｄＨｋ／ｄｐ２を計算する。
【００５１】
測定残差微分係数計算部５０４は、観測行列微分係数計算部５０３によって計算された、パラメータｐ２に関する測定行列の偏微分係数ｄＨｋ／ｄｐ２と、状態予測部３０１によって計算された予測状態ｘ＾ｋと、観測行列計算部３０７によって計算された観測行列Ｈｋと、予測状態微分係数計算部５０１によって計算された、パラメータｐ２に関する予測状態の偏微分係数ｄｘ＾ｋ／ｄｐ２とから、パラメータｐ２に関する測定残差の偏微分係数ｄｅｋ／ｄｐ２を、
【数２１】
0006388563.tif 000022
によって計算する。
【００５２】
予測共分散微分係数計算部５０５は、１周期前に推定共分散微分係数計算部５１０で計算され、推定共分散微分係数記憶部５０６に記憶された、パラメータｐ２に関する推定共分散の偏微分係数ｄＰｋ−１／ｄｐ２から、パラメータｐ２に関する予測共分散の偏微分係数ｄＰ＾ｋ／ｄｐ２を、
【数２２】
0006388563.tif 000023
に従って計算する。
【００５３】
残差共分散微分係数計算部５０７は、観測行列微分係数計算部５０３によって計算された観測行列の偏微分係数ｄＨｋ／ｄｐ２と、共分散予測部３０５によって計算された予測共分散Ｐ＾ｋと、観測行列計算部３０７によって計算された観測行列Ｈｋと、予測共分散微分係数計算部５０５によって計算された、パラメータｐ２に関する予測共分散の偏微分係数ｄＰ＾ｋ／ｄｐ２とから、パラメータｐ２に関する残差共分散の偏微分係数ｄＳｋ／ｄｐ２を、
【数２３】
0006388563.tif 000024
によって計算する。
【００５４】
Ｋａｌｍａｎゲイン微分係数計算部５０８は、予測共分散微分係数計算部５０５によって計算された、パラメータｐ２に関する予測共分散の偏微分係数ｄＰ＾ｋ／ｄｐ２と、観測行列微分係数計算部５０３によって計算された、パラメータｐ２に関する観測行列の偏微分係数ｄＨｋ／ｄｐ２と、残差共分散計算部３０８によって計算された残差共分散Ｓｋと、共分散予測部３０５によって計算された予測共分散Ｐ＾ｋと、観測行列微分係数計算部５０３によって計算された、パラメータｐ２に関する観測行列の偏微分係数ｄＨｋ／ｄｐ２と、残差共分散微分係数計算部５０７によって計算された、パラメータｐ２に関する測定残差共分散の偏微分係数ｄＳｋ／ｄｐ２とから、パラメータｐ２に関するＫａｌｍａｎゲインの偏微分係数ｄＫｋ／ｄｐ２を、
【数２４】
0006388563.tif 000025
に従い計算する。
【００５５】
推定状態微分係数計算部５０９は、予測状態微分係数計算部５０１によって計算された予測状態の偏微分係数ｄｘ＾ｋ／ｄｐ２と、Ｋａｌｍａｎゲイン微分係数計算部５０８によって計算された、パラメータｐ２に関するＫａｌｍａｎゲインの偏微分係数ｄＫｋ／ｄｐ２と、測定残差計算部３０４によって計算された測定残差ｅｋと、Ｋａｌｍａｎゲイン計算部３０９によって計算されたＫａｌｍａｎゲインＫｋと、測定残差微分係数計算部５０４によって計算された測定残差の偏微分係数ｄｅｋ／ｄｐ２とから、パラメータｐ２に関する推定状態の偏微分係数ｄｘｋ／ｄｐ２を、
【数２５】
0006388563.tif 000026
に従って計算し、結果を推定状態微分係数記憶部５０２に記憶させる。
【００５６】
推定共分散微分係数計算部５１０は、予測共分散微分係数計算部５０５によって計算された、パラメータｐ２に関する予測共分散の偏微分係数ｄＰ＾ｋ／ｄｐ２と、Ｋａｌｍａｎゲイン微分係数計算部５０８によって計算されたＫａｌｍａｎゲインの偏微分係数ｄＫｋ／ｄｐ２と、観測行列計算部３０７によって計算された観測行列Ｈｋと、共分散予測部３０５によって計算された予測共分散Ｐ＾ｋと、Ｋａｌｍａｎゲイン計算部３０９によって計算されたＫａｌｍａｎゲインＫｋと、観測行列微分係数計算部５０３によって計算された、パラメータｐ２に関する観測行列の偏微分係数ｄＨｋ／ｄｐ２とから、パラメータｐ２に関する推定共分散の偏微分係数ｄＰｋ／ｄｐ２を、
【数２６】
0006388563.tif 000027
に従い計算し、結果を推定共分散微分係数記憶部５０６に記憶させる。
【００５７】
ＳｏＣ微分係数計算部５１１は、予め測定した蓄電池１０３の定格ＦＣＣと、推定状態微分係数計算部５０９によって計算された、パラメータｐ２に関する推定状態の偏微分係数ｄｘｋ／ｄｐ２の第２要素ｄｑｂ，ｋ／ｄｐ２と、状態推定部３１０によって推定された推定状態ｘｋの第２要素である充電量ｑｂ，ｋとから、第２係数ｄｓｋ／ｄｐ２を、
【数２７】
0006388563.tif 000028
に従い計算する。
【００５８】
電流積算部１１３は、初期時刻からの電流１０４の積算値を、
【数２８】
0006388563.tif 000029
に従って積算する。
【００５９】
さて、もし、蓄電池１０３のＳｏＣが正しく推定され、電流測定値がオフセット誤差を持たず、定格ＦＣＣが蓄電池１０３の実際のＦＣＣに等しいならば、
【数２９】
0006388563.tif 000030
が成り立つ。ただし、ｑ０は初期充電量である。
【００６０】
しかし、実際には、電流測定値はオフセット誤差を持ち、定格ＦＣＣは実際のＦＣＣからずれており、さらに、ＳｏＣ推定値１０８は、電流測定値のオフセット誤差およびＦＣＣの容量誤差に起因する偏り誤差を持っている。ゆえに、式（３１）は、
【数３０】
0006388563.tif 000031
と補正される。パラメータｐ１，ｐ２に関する高次の項を無視して整理すると、式（３２）は、
【数３１】
0006388563.tif 000032
と変形される。
【００６１】
ここで（３３）の左辺は各時刻において定数であり、右辺は各時刻においてパラメータｐ１，ｐ２および充電量ｑ０の１次式である。ゆえに、目的関数
【数３２】
0006388563.tif 000033
を最小にするようなＲｅｃｕｒｓｉｖｅ  ｌｅａｓｔ  ｓｑｕａｒｅｓ  ｆｉｌｔｅｒ  を用いて、平均二乗誤差の観点から、最適なパラメータ推定値ｐ＾ｋを推定することができる。ただし、λは、０より大きく１以下の実定数である忘却係数である。
【００６２】
本実施の形態１において、ベクトルｕｋは説明変数であり、ＳｏＣ推定値ｓｋ、第１係数ｄｓｋ／ｄｐ１、第２係数ｄｓｋ／ｄｐ２から、
【数３３】
0006388563.tif 000034
によって定められる。また、スカラーｙｋは観測変数であり、本実施の形態１においては、
【数３４】
0006388563.tif 000035
で定義される。
【００６３】
また、ベクトルｐ＾ｋは、本実施の形態１において、パラメータｐ１の推定値を第１要素とし、パラメータｐ２の推定値を第２要素とし、初期充電量ｑ０の推定値を第３要素とするパラメータ推定ベクトルである。この事実を利用し、蓄電池のオフセット誤差および容量誤差を補正することができるパラメータｐ１，ｐ２を推定することが、本発明の要点である。
【００６４】
図６は、本発明の実施の形態１に係る誤差推定部１１５の全体構成図である。説明変数計算部６０１は、ＳｏＣ推定部１０７によって推定されたＳｏＣ推定値１０８と、第１係数計算部１０９によって計算された第１係数１１０と、第２係数計算部１１１によって計算された第２係数１１２とから、式（３５）に基づいて説明変数ｕｋを計算する。
【００６５】
観測変数計算部６０２は、ＳｏＣ推定値１０８と、電流積算部１１３によって計算された電流積算値１１４とから、式（３６）に従って観測変数ｙｋを計算する。
【００６６】
ＲＬＳゲイン計算部６０３は、説明変数計算部６０１によって計算された説明変数ｕｋと、逆分散計算部６０５によって計算され、逆分散記憶部６０４に記憶されている１周期前の逆分散Ｙｋ−１とから、ＲＬＳゲインｇｋを、
【数３５】
0006388563.tif 000036
によって計算する。
【００６７】
逆分散計算部６０５は、逆分散計算部６０５によって計算され、逆分散記憶部６０４に記憶されている１周期前の逆分散Ｙｋ−１と、ＲＬＳゲイン計算部６０３によって計算されたＲＬＳゲインｇｋと、説明変数計算部６０１によって計算された説明変数ｕｋとから、逆分散Ｙｋを、
【数３６】
0006388563.tif 000037
に従い計算する。ただし、Ｉｄは、単位行列を表す。
【００６８】
パラメータ推定部６０６は、説明変数計算部６０１によって計算された説明変数ｕｋと、観測変数計算部６０２によって計算された観測変数ｙｋと、ＲＬＳゲイン計算部６０３によって計算されたＲＬＳゲインｇｋと、パラメータ記憶部６０７に記憶されている１周期前のパラメータ推定ベクトルｐ＾ｋ−１とから、パラメータ推定ベクトルｐ＾ｋを、
【数３７】
0006388563.tif 000038
に従い計算し、パラメータ記憶部６０７に記憶する。
【００６９】
電流測定値のオフセット誤差ΔＩは、ΔＩ＝Ｉｔｙｐ×ｐ１で求められ、ここで、パラメータｐ１はパラメータ推定ベクトルｐ＾ｋの第１要素である。また、蓄電池のＦＣＣは、Ｆｃｃ＝Ｆ〜ｃｃ（１−ｐ２）で求められ、ここで、パラメータｐ２はパラメータ推定ベクトルｐ＾ｋの第２要素である。
【００７０】
オフセット計算部６０８および容量計算部６０９は、この計算を実行し、それぞれ、電流測定値のオフセット誤差１１６および蓄電池の容量誤差１１７の推定値を計算する。
【００７１】
図７は、本発明の実施の形態１に係る説明変数計算部６０１のブロック線図である。時刻ｔｋは、１周期前の時刻ｔｋ−１にサンプリング周期ｔｓを加算することによって計算される。
【００７２】
説明変数ｕｋの第１要素は、時刻ｔｋに定数Ｉｔｙｐをかけたものと、第１係数計算部１０９によって計算された第１係数１１０に、定格ＦＣＣＦｃｃをかけたものとの差をとることにより計算される。説明変数ｕｋの第２要素は、ＳｏＣ推定部１０７によって推定されたＳｏＣ推定値１０８に、第２係数計算部１１１によって計算された第２係数１１２を足しあわせたものに定格ＦＣＣＦｃｃをかけ、符号を逆転することによって計算される。説明変数ｕｋの第３要素は定数−１に設定されている。
【００７３】
図８は、本発明の実施の形態１に係る観測変数計算部６０２のブロック線図である。観測変数ｙｋは、電流積算部１１３によって計算された電流積算値１１４と、ＳｏＣ推定部１０７によって推定されたＳｏＣ推定値ｓｋに、定格ＦＣＣＦｃｃをかけたものとの差を取ることによって計算される。
【００７４】
図９は、本発明の実施の形態１に係る蓄電池パラメータ推定方法のフローチャートである。フローチャート中の各ステップは、図１に示す蓄電池パラメータ推定装置１０１の各部が行う処理に対応する。ただし、本実施の形態１における蓄電池パラメータ推定方法の実行順序は、図９に示したものに限定されない。各ステップの依存関係を壊さない限りステップの実行順序の入れ替えが許容される。
【００７５】
図１０は、本発明の実施の形態１に係るＳｏＣ推定部１０７が行う処理の１周期のフローチャートである。フローチャート中の各ステップは、図３に示すＳｏＣ推定部１０７の各部が行う処理に対応する。ただし、本実施の形態１におけるＳｏＣ推定部１０７の実行順序は図１０に示したものに限定されない。各ステップの依存関係を壊さない限りステップの実行順序の入れ替えが許容される。
【００７６】
図１１は、本発明の実施の形態１に係る第１係数計算部１０９が行う処理の１周期のフローチャートである。フローチャート中の各ステップは、図４に示す第１係数計算部１０９の各部が行う処理に対応する。ただし、本実施の形態１における第１係数計算部１０９の実行順序は図１１に示したものに限定されない。各ステップの依存関係を壊さない限りステップの実行順序の入れ替えが許容される。
【００７７】
図１２は、本発明の実施の形態１に係る第２係数計算部１１１が行う処理の１周期のフローチャートである。フローチャート中の各ステップは、図５に示す第２係数計算部１１１の各部が行う処理に対応する。ただし、本実施の形態１における第２係数計算部１１１の実行順序は図１２に示したものに限定されない。各ステップの依存関係を壊さない限りステップの実行順序の入れ替えが許容される。
【００７８】
図１３は、本発明の実施の形態１に係る誤差推定部１１５が行う処理の１周期のフローチャートである。フローチャート中の各ステップは、図６に示す誤差推定部１１５の各部が行う処理に対応する。ただし、本実施の形態１における誤差推定部１１５の実行順序は図１３に示したものに限定されない。各ステップの依存関係を壊さない限りステップの実行順序の入れ替えが許容される。
【００７９】
また、本実施の形態１の詳細な計算方法は、本実施の形態１に示されるものに限定されない。例えば、式（３８）は、説明変数ｕｋの重み付き共分散行列Ｒｋの逆行列である逆分散Ｙｋを再帰的に計算する式であるが、計算量は増大するものの通常の行列反転を用いて、
【数３８】
0006388563.tif 000039
と計算してもよい。同様に、パラメータ推定行列ｐ＾ｋは、
【数３９】
0006388563.tif 000040
で定義されるベクトルｒｋを用いて、
【数４０】
0006388563.tif 000041
と計算してもよい。あるいは、本実施の形態１において忘却係数λは定数としているが、これをｋに応じて可変としてもよい。
【００８０】
更に、本実施の形態１において、ＳｏＣ推定部１０７は、Ｋａｌｍａｎフィルタを用いて、蓄電池１０３を流れる電流１０４と、蓄電池１０３の端子間電圧１０６とから、ＳｏＣ推定値１０８を求めているが、本実施の形態１によるＳｏＣ推定部１０７は、Ｋａｌｍｎａｎフィルタを用いた実装に限定されるものではない。電流測定値のオフセット誤差および定格ＦＣＣの容量誤差に関するＳｏＣ推定値１０８の偏微分係数が計算できる限りにおいて、本発明は適用可能である。
【００８１】
本実施の形態１における蓄電池パラメータ推定装置１０１は、例えば、メモリと入出力機能を備えたマイクロコントローラとして実現して良い。その場合、電流計測部１０２、電圧計測部１０５を除く各部は、マイクロコントローラに組み込まれたソフトウェアとして実現される。
【００８２】
以上のように、実施の形態１によれば、蓄電池の電流測定値および予め測定された蓄電池の定格ＦＣＣの何れかまたは両方に誤差が含まれ、この誤差は、蓄電池の等価回路を用いて、蓄電池の端子間電圧の測定値から推定されるＳｏＣ推定値に基づいて推定される。この結果、電流測定値に誤差が含まれる場合でも、電流測定値のオフセット誤差およびＦＣＣの容量誤差を高精度に推定することができる蓄電池パラメータ推定装置を得ることができる。
【００８３】
これにより、より効率のよいエネルギー管理が可能となり、省エネルギーに資する。また、高精度な蓄電池パラメータの推定により、蓄電池の劣化状態を継続的に把握することができるので、蓄電池の寿命の管理や、蓄電システムのメンテナンスコスト削減等にも資する。
【００８４】
実施の形態２．
図１４は、本発明の実施の形態２に係る蓄電池パラメータ推定装置１０１の全体構成図である。図１４に示す本実施の形態２の蓄電池パラメータ推定装置１０１は、先の実施の形態１の図１と比較して、図２に示した蓄電池１０３の等価回路の特性パラメータＲｄ，Ｃｄ，Ｒ０の誤差を追加考慮して電流測定値のオフセット誤差１１６およびＦＣＣの容量誤差１１７を推定する点が異なっている。
【００８５】
具体的には、ＳｏＣ推定部１０７において用いる、予め測定した標準パラメータＲｄ，Ｃｄ，Ｒ０が、蓄電池の使用時における実際の電池と比較して誤差を持つと考え、これらを、
【数４１】
0006388563.tif 000042
のようにパラメトライズする。その上で、パラメータｐ３、ｐ４、ｐ５にそれぞれ対応する第３係数計算部１４０１、第４係数計算部１４０２、第５係数計算部１４０３により、先の実施の形態１と同様に、ＳｏＣ推定値ｓｋの、パラメータｐ３，ｐ４，ｐ５に対するに関する変化率である第３係数、第４係数、第５係数をそれぞれ計算する。
【００８６】
誤差推定部１４０４は、式（３１）に加え、
【数４２】
0006388563.tif 000043
の関係を利用する。先の実施の形態１の第１係数および第２係数を含む係数情報に加えて、第３係数、第４係数、第５係数を追加することにより、式（３２）は、
【数４３】
0006388563.tif 000044
のように拡張される。また、式（４４）、式（４５）に同様の変形を行うことにより、
【数４４】
0006388563.tif 000045
を得る。ただし、
【数４５】
0006388563.tif 000046
である。パラメータに関する高次の項を無視すると、式（４６）、式（４７）および式（４８）は、パラメータｐｊに関し線形な関係に書き換えられる。
【００８７】
本実施の形態２における誤差推定部１４０４は、式（３４）に代えて目的関数Ｊ、
【数４６】
0006388563.tif 000047
を最小化する。ただし、Σは、所定の正定対称な重み行列である。また、観測変数ｙｌは、
【数４７】
0006388563.tif 000048
と定義される。また、説明変数ｕｋは、
【数４８】
0006388563.tif 000049
とおいて、行列［ｕ１，ｋｕ２，ｋｕ３，ｋ］となる。式（３７）に代えて、
【数４９】
0006388563.tif 000050
と定義されるＲＬＳゲインｇｋを、式（３８）、式（３９）と合わせて、先の実施の形態１と全く同様にパラメータ推定ベクトルｐ＾ｋを推定できる。
【００８８】
図１５は、本発明の実施の形態２に係る蓄電池パラメータ推定方法のフローチャートである。フローチャート中の各ステップは、図１４に示す蓄電池パラメータ推定装置１０１の各部が行う処理に対応する。ただし、本実施の形態２における蓄電池パラメータ推定装置１０１の実行順序は、図１５に示したものに限定されない。各ステップの依存関係を壊さない限りステップの実行順序の入れ替えが許容される。
【００８９】
以上のように、実施の形態２によれば、図２に示した蓄電池１０３の等価回路の特性パラメータＲｄ，Ｃｄ，Ｒ０の誤差を考慮して電流測定値のオフセット誤差１１６およびＦＣＣの容量誤差１１７を推定することができる。
【００９０】
なお、本実施の形態２による蓄電池パラメータ推定装置１０１は、図２に示した等価回路の特性パラメータを補正する構成に限られるものではない。例えば、複数のＣＲ回路や、拡散インピーダンス等を含む回路を等価回路として考え、それに含まれるパラメータを補正する構成をとることも可能である。
【００９１】
また、抵抗や拡散係数等のパラメータの温度依存性をＡｒｒｈｅｎｉｕｓの式などで表現し、そこに含まれる頻度因子や活性化エネルギーなどをパラメータとして補正することも可能である。等価回路の特性パラメータＲｄ，Ｃｄ，Ｒ０の一部のみを推定するような構成も可能である。
【００９２】
実施の形態３．
先の実施の形態１では、式（３５）および式（３６）によって定められる説明変数ｕｋ、および観測変数ｙｋを基に、パラメータ推定ベクトルｐ＾ｋを推定している。ここで、説明変数ｕｋの第１要素は、ｔｋに比例する項を含んでおり、これは時間の経過に従って際限なく増大する。このような際限なく増大する量を計算過程に含むことは、組み込みシステムにおけるソフトウェアの実装に適したものではない。
【００９３】
本実施の形態３は、このような問題を解決するためのものであり、観測変数ｙｋと説明変数ｕｋをアフィン空間上で考えることにより、際限なく増大する量を直接取り扱うこと無く、先の実施の形態１および２と同様の効果を得るものである。以下、本実施の形態３による蓄電池パラメータ推定装置１０１を、本実施の形態１の蓄電池パラメータ推定装置１０１との差異を中心に説明する。
【００９４】
まず、観測変数ｙｋ、説明変数ｕｋの重み付き時間平均を、
【数５０】
0006388563.tif 000051
と定義する。ただし、Ｓｋは、初項１、公比λの等比数列の、最初のｋ項分の和である。ここで、組（ｙｋ，ｕｋ）が属するベクトル空間をアフィン空間とみなし、時刻ｔｋにおける局所座標（ｚｌ，ｖｌ）を、ｚｌ：＝ｙｌ−ｙｋ，ｖｌ：＝ｕｌ−ｕｋととる。すると、局所座標表示した観測変数ｙｋと説明変数ｕｋの重み付き時間平均は、
【数５１】
0006388563.tif 000052
で計算できる。差分ｙｋ−ｙｋ−１、ｕｋ−ｕｋ−１はともに有限の値であり、λＳｋ−１／Ｓｋは１未満であるため、ｚｋ，ｖｋは常に有限の値を取る。
【００９５】
ここで、観測変数を、ｙｋの代わりに＜ｚｋ＞とし、説明変数を、ｕｋの代わりに＜ｖｋ＞とすることで、先の実施の形態１の誤差推定部１１５と同様に、パラメータ推定ベクトルｐ＾ｋを計算できる。ただし、本実施の形態３においては、パラメータをアフィン空間で考えたことから、ｑ０を推定することはできない。従って、説明変数ｕｋは、最後の要素−１を除いた２次元ベクトルとし、パラメータ推定ベクトルｐ＾ｋも２次元ベクトルとする。
【００９６】
図１６は、本発明の実施の形態３に係る蓄電池パラメータ推定装置１０１の全体構成図である。本実施の形態１の蓄電池パラメータ推定装置１０１は、電流積算部１１３が不要となり、電流計測部１０２によって計測された電流１０４が、誤差推定部１６０１に直接入力されることを除き、先の実施の形態１に記載の蓄電池パラメータ推定装置１０１と同様の構成を取る。
【００９７】
図１７は、本発明の実施の形態３に係る誤差推定部１６０１の全体構成図である。本実施の形態３の誤差推定部１６０１は、観測変数計算部１７０２に電流積算値１１４が入力される代わりに電流１０４が入力されることを除いて、先の実施の形態１の誤差推定部１１５と同様の構成を取る。
【００９８】
図１８は、本発明の実施の形態３に係る説明変数計算部１７０１のブロック線図である。説明変数ｕｋ−ｕｋ−１の第１要素は、第１係数計算部１０９によって計算された第１係数１１０から、１周期前の第１係数１１０を引き、定格ＦＣＣＦｃｃを掛けた値を、サンプリング周期ｔｓに定数Ｉｔｙｐをかけた値から引くことで計算される。
【００９９】
説明変数ｕｋ−ｕｋ−１の第２要素は、ＳｏＣ推定部１０７によって計算されたＳｏＣ推定値１０８と、第２係数計算部１１１によって計算された第２係数１１２との和を計算し、計算した和の前回サンプル時の値との差分を取り、さらに定格ＦＣＣＦｃｃをかけ、符号を逆転することで計算される。
【０１００】
一方、λｋは、前回サンプル時に計算されたλｋ−１に所定の忘却係数λを掛けることで計算され、初項１、公比λの等比数列の、最初のｋ項分の和である級数Ｓｋは、前回サンプル時に計算された級数Ｓｋ−１にλｋを加えることで計算され、さらにＳｋとＳｋ−１の比をとり、この比に係数λを掛けることでλＳｋ−１／Ｓｋが計算される。説明変数＜ｖｋ＞は、前回サンプル時の＜ｖｋ−１＞と、説明変数ｕｋ−ｕｋ−１との差分を取り、さらに、λＳｋ−１／Ｓｋを掛けあわせることで計算される。
【０１０１】
図１９は、本発明の実施の形態３に係る観測変数計算部１７０２のブロック線図である。ＳｏＣ推定部１０７によって計算されたＳｏＣ推定値１０８に、定格ＦＣＣＦｃｃをかけ、前回サンプル時との差分を取り、この差分を、電流１０４にサンプリング周期ｔｓをかけたものから引くことで、観測変数ｙｋ−ｙｋ−１が計算される。観測変数＜ｚｋ＞は、前回サンプル時の＜ｚｋ−１＞と、観測変数ｙｋ−ｙｋ−１との差分を取り、さらに、λＳｋ−１／Ｓｋを掛けあわせることで計算される。
【０１０２】
本実施の形態３を先の実施の形態２と組み合わせる場合、値ｙｋをベクトルと読み替え、ベクトルｕｋを行列と読み替えればよい。また、説明変数、観測変数の計算において、時刻ｔｋおよび電流積算値ｑｃｃ，ｋに係る項以外は、定義通りにｕｋ−ｕｋ−１やｙｋ−ｙｋ−１等が計算される必要がある。
【０１０３】
図２０は、本発明の実施の形態３に係る蓄電池パラメータ推定方法のフローチャートである。フローチャート中の各ステップは、図１６に示す蓄電池パラメータ推定装置１０１の各部が行う処理に対応する。ただし、本実施の形態３における蓄電池パラメータ推定装置１０１の実行順序は、図２０に示したものに限定されない。各ステップの依存関係を壊さない限りステップの実行順序の入れ替えが許容される。
【０１０４】
図２１は、本発明の実施の形態３に係る誤差推定部１６０１が行う処理の１周期のフローチャートである。フローチャート中の各ステップは、図１７に示す誤差推定部１６０１の各部が行う処理に対応する。ただし、本実施の形態３における誤差推定部１６０１の実行順序は、図２１に示したものに限定されない。各ステップの依存関係を壊さない限りステップの実行順序の入れ替えが許容される。
【０１０５】
以上のように、実施の形態３によれば、時間の経過に従って際限なく増大する量を直接取り扱うこと無く、先の実施の形態１と同様の効果を得ることができる。
【０１０６】
実施の形態４．
先の実施の形態１では、電流測定値のオフセット誤差１１６、およびＦＣＣの容量誤差１１７の推定にＲｅｃｕｒｓｉｖｅ  Ｌｅａｓｔ  Ｓｑｕａｒｅｓ  Ｆｉｌｔｅｒ  （ＲＬＳフィルタ）を用いている。この手法は、観測変数ｙｋにのみランダムな誤差が重畳し、説明変数ｕｋのランダム誤差が無視できるほど小さい場合に、平均二乗誤差の観点から最適な推定を与える。
【０１０７】
しかし、説明変数ｕｋが大きなランダム誤差を持っていた場合、推定値がバイアスを持つことが知られている。本実施の形態４は、電流測定値のオフセット誤差１１６およびＦＣＣの容量誤差１１７の推定にＲｅｃｕｒｓｉｖｅ  Ｔｏｔａｌ  Ｌｅａｓｔ  Ｓｑｕａｒｅｓ  Ｆｉｌｔｅｒ（ＲＴＬＳフィルタ）を用いることで、この問題を解決するものである。
【０１０８】
本実施の形態４において、蓄電池パラメータ推定装置１０１は、式（５１）に代えて目的関数Ｊ、
【数５２】
0006388563.tif 000053
を最小にするパラメータ推定ベクトルｐ＾ｋを求める。ただし、δｋ、εｋは、任意のｋに対して、説明変数をｕｋ、観測変数をｙｋとして、
【数５３】
0006388563.tif 000054
を満たす。今、ｙ＾ｌ：＝ｙｌ−δｌ、ｕ＾ｌ：＝ｕｌ−εｌとおくと、式（６１）から、
【数５４】
0006388563.tif 000055
が成り立つ。これを用いると、式（６０）は、
【数５５】
0006388563.tif 000056
と変形することができる。式（６３）を最小化するベクトルｕ＾ｌは、
【数５６】
0006388563.tif 000057
の解として得られる。
【０１０９】
今、式（６４）の解である
【数５７】
0006388563.tif 000058
は、
【数５８】
0006388563.tif 000059
と表現される。これを式（６３）に代入して整理すると、
【数５９】
0006388563.tif 000060
と書き換えられる。さらに、Λ１／２を、行列ΛのＣｈｏｌｅｓｋｙ分解とし、Ｂｋ：＝Λ１／２Ａｋ、ζｌ：＝Λ１／２ξｌとおくと、
【数６０】
0006388563.tif 000061
を得る。
【０１１０】
ここで、ＡｋＴの零空間が１次元であることから、ある零でないベクトルｂｋが存在して、ＢｋＴｂｋ＝０が成り立つ。このようなｂｋは、スカラー倍の自由度を除いて一意である。射影行列の性質から、このようなベクトルｂｋに関して、
【数６１】
0006388563.tif 000062
が成り立つ。従って、式（６８）は、
【数６２】
0006388563.tif 000063
と書き換えられる。さて、ＢｋＴｂｋ＝０、すなわちＡｋＴΛ１／２ｂｋ＝０であるから、ａｋ：＝Λ１／２ｂｋとおくと、式（７０）はさらに、
【数６３】
0006388563.tif 000064
と変形できる。
【０１１１】
ゆえに、ａｋは、組（Ｒｋ、Λ）に関する一般化固有値問題の解である最小固有値に対応する固有ベクトルとして求められる。また、ＡｋＴａｋ＝０であるから、ベクトルａｋの第一要素をａ０，ｋとおき、残りの要素をａ＾ｋとおくと、ｐ＾ｋａ０，ｋ＝ａ＾ｋである。ゆえに、パラメータ推定ベクトルｐ＾ｋは、ａ＾ｋ／ａ０，ｋと表現される。
【０１１２】
固有ベクトルの求め方としては、公知の手法を用いてよい。例えば、Ｈｏｕｓｅｈｏｌｄｅｒ変換法、Ｇｉｖｅｎｓ回転法、Ａｒｎｏｌｄｉ法などが知られている。今、共分散行列Ｒｋは、正定対称行列であるから、Ｌａｎｃｚｏｓ法も適用可能である。さらに最小固有値に対応する固有ベクトルさえ求まればよいことから、逆べき乗法やＲａｙｌｅｉｇｈ  ｑｕｏｔｉｅｎｔ  ｉｔｅｒａｔｉｏｎ法も利用できる。
【０１１３】
本実施の形態４における蓄電池パラメータ推定装置１０１を、主に、先の実施の形態１の蓄電池パラメータ推定装置１０１との差異を中心に説明する。
【０１１４】
図２２は、本発明の実施の形態４に係る蓄電池パラメータ推定装置１０１の全体構成図である。本実施の形態４における蓄電池パラメータ推定装置１０１は、誤差推定部２２０１を除き、先の実施の形態１の蓄電池パラメータ推定装置１０１と同様の構成を取る。
【０１１５】
図２３は、本発明の実施の形態４に係る誤差推定部２２０１の全体構成図である。説明変数計算部６０１と観測変数計算部６０２は、先の実施の形態１と同様である。
【０１１６】
共分散計算部２３０１は、共分散記憶部２３０４に記憶された前回サンプル時の共分散行列Ｒｋ−１と、説明変数計算部６０１によって計算された説明変数ｕｋと、観測変数計算部６０２によって計算された観測変数ｙｋとから、共分散行列Ｒｋを次の再帰式、
【数６４】
0006388563.tif 000065
に従って計算し、結果を共分散記憶部２３０４に記憶させる。
【０１１７】
固有ベクトル計算部２３０２は、共分散計算部２３０１によって計算された共分散行列Ｒｋと、所定の重み行列Λとに基づいて、式（７１）の一般化固有値問題を解き、最小固有値に対応する固有ベクトルａｋを求める。
【０１１８】
パラメータ計算部２３０３は、ａｋの第１要素で第２要素以降の残りの要素を割ることで、パラメータ推定値ｐ＾ｋを算出する。オフセット計算部６０８と容量計算部６０９は、先の実施の形態１と同様である。
【０１１９】
図２４は、蓄電池パラメータ推定方法のフローチャートである。フローチャート中の各ステップは、図２２に示す蓄電池パラメータ推定装置１０１の各部が行う処理に対応する。ただし、本実施の形態４における蓄電池パラメータ推定装置１０１の実行順序は、図２４に示したものに限定されない。各ステップの依存関係を壊さない限りステップの実行順序の入れ替えが許容される。
【０１２０】
図２５は、本発明の実施の形態４に係る誤差推定部２２０１が行う処理の１周期のフローチャートである。フローチャート中の各ステップは、図２３に示す誤差推定部２２０１の各部が行う処理に対応する。ただし、本実施の形態４における誤差推定部２２０１の実行順序は、図２５に示したものに限定されない。各ステップの依存関係を壊さない限りステップの実行順序の入れ替えが許容される。
【０１２１】
本実施の形態４を、先の実施の形態２と組み合わせることも可能である。このためには、Σを正定対称な重み行列として、共分散行列Ｒｋを、
【数６５】
0006388563.tif 000066
のように計算すればよい。
【０１２２】
本実施の形態４を、先の実施の形態３と組み合わせることも可能である。このことは、先の実施の形態１と先の実施の形態３の差異が、説明変数計算部６０１、観測変数計算部６０２に限定されることを考慮すれば明らかである。当然ながら、先の実施の形態３に記載したように、先の実施の形態２と組み合わせた先の実施の形態３を、本実施の形態４とさらに組み合わせることも可能である。
【０１２３】
以上のように、実施の形態４によれば、電流１０４の測定値や電圧１０６が大きなランダム誤差を持っていた場合においても、先の実施の形態１と同様の効果を得ることができる。
【０１２４】
実施の形態５．
先の実施の形態４の蓄電池パラメータ推定装置１０１は、固有ベクトル計算部２３０２において、一般化固有値問題を解いており、一般的には、これは計算量が大きい処理である。そこで、本実施の形態５では、非特許文献１に開示された手法を応用し、共分散行列Ｒｋがサンプル周期ごとに大きく変化しないことを利用して、計算量を削減する構成を開示する。
【０１２５】
まず、パラメータ推定ベクトルｐ＾ｋは、ａ＾ｋ／ａ０，ｋと表現され、パラメータ推定ベクトルｐｋはスカラー倍の自由度が有ることから、ｐｋ＝［−１ａ＾ｋＴ］Ｔとパラメトライズできる。
【０１２６】
ここで、以下の更新計画
【数６６】
0006388563.tif 000067
を導入する。ただし、θｋはスカラーであり、ｗｋは、列｛ｗｋ，ｗｋ−１，…｝が効率よくＡｋの像を張るような時系列である。
【０１２７】
そこで、非特許文献１では、ｗｋとして説明変数ｕｋを用いることとしているが、実施の形態４では、ランダムな零でないベクトルが用いられている。他の選択肢としては、例えば、過去ｎ−１個のｗｋと直交する方向ベクトルを生成して利用する、等が考え得る。
【０１２８】
式（７４）に含まれるスカラーθｋは、式（７１）に示した目的関数Ｊを最小化するようなスカラーθｋを計算する。式（７１）に式（７４）を代入すると、
【数６７】
0006388563.tif 000068
のように書き換えられる。ただし、
【数６８】
0006388563.tif 000069
である。
【０１２９】
式（７５）を最小化するスカラーθｋは、方程式、
【数６９】
0006388563.tif 000070
の解として与えられる。ここで、分母Ｄｋ２は常に正であるから、分子が０となるようなθｋを求めればよい。分子Ｎｋおよび分母Ｄｋを展開すると、
【数７０】
0006388563.tif 000071
である。
【０１３０】
簡単のため、分子Ｎｋ＝Ｎ２，ｋθｋ２＋２Ｎ１，ｋθｋ＋Ｎ０，ｋとおき、分母Ｄｋ＝Ｄ２，ｋθｋ２＋２Ｄ１，ｋθｋ＋Ｄ０，ｋとおくと、それぞれのスカラーθｋに関する微分は、ｄＮｋ／ｄθｋ＝２Ｎ２，ｋθｋ＋２Ｎ１，ｋ、ｄＤｋ／ｄθｋ＝２Ｄ２，ｋθｋ＋２Ｄ１，ｋと表現される。ゆえに式（７７）の分子は、
【数７１】
0006388563.tif 000072
と書ける。
【０１３１】
簡単のため、式（８０）をさらに２αｋθｋ２＋２βｋθｋ＋２γｋθｋとおく。これは２次式であるから、おそらく相異なる２つの実数根
【数７２】
0006388563.tif 000073
を有する。
【０１３２】
目的関数Ｊの増減表を書くと、表１、表２のようになるから、明らかに、θ＋ｋがＪを最小にする解である。なお、増減表においてＣは、定数Ｎ２，ｋ／Ｄ２，ｋである。
【０１３３】
【表１】
0006388563.tif 000074
【０１３４】
【表２】
0006388563.tif 000075
【０１３５】
増減表から、２次式αｋθｋ２＋βｋθｋ＋γｋθｋは、必ず相異なる２つの実数根を持つ。しかし、αｋが、０または０に非常に近い時には、根を求めることは数値的に不安定になる。この場合、θｋ＝０とすればよい。
【０１３６】
本実施の形態５における蓄電池パラメータ推定装置１０１を、主に先の実施の形態４の蓄電池パラメータ推定装置１０１との差異を中心に説明する。
【０１３７】
図２６は、本発明の実施の形態５に係る蓄電池パラメータ推定装置１０１の全体構成図である。誤差推定部２６０１を除き、先の実施の形態４の蓄電池パラメータ推定装置１０１と同様の構成を取る。
【０１３８】
図２７は、本発明の実施の形態５に係る誤差推定部２６０１の全体構成図である。説明変数計算部６０１と観測変数計算部６０２は、先の実施の形態１と同様である。共分散計算部２３０１と共分散記憶部２３０４は、先の実施の形態４と同様である。
【０１３９】
係数計算部２７０１は、説明変数計算部６０１によって計算された説明変数ｕｋと、共分散計算部２３０１によって計算された共分散行列Ｒｋと、パラメータ記憶部２７０４に記憶された前回サンプル時のパラメータ推定値ｐ＾ｋ−１とから、式（７８）、（７９）、（８０）により、２次式の係数αｋ、βｋ、γｋを計算する。
【０１４０】
２次式求解部２７０２は、係数計算部２７０１によって計算された係数を持つ２次式の実数根θ＋ｋを、式（８１）の第２式に基づいて求める。
【０１４１】
パラメータ計算部２７０３は、係数計算部２７０１によって計算された２次式の係数と、２次式求解部２７０２によって求められた２次式の実数根θ＋ｋとから、式（７４）により、パラメータ推定ベクトルｐ＾ｋを計算し、パラメータ記憶部２７０４に記憶させる。ただし、２次式の２次の係数がほとんど０であるときは、ｐ＾ｋ＝ｐ＾ｋ−１とする。オフセット計算部６０８と容量計算部６０９は、先の実施の形態１と同様である。
【０１４２】
図２８は、本発明の実施の形態５に係る蓄電池パラメータ推定方法のフローチャートである。フローチャート中の各ステップは、図２６に示す蓄電池パラメータ推定装置１０１の各部が行う処理に対応する。ただし、本実施の形態５における蓄電池パラメータ推定装置１０１の実行順序は、図２８に示したものに限定されない。各ステップの依存関係を壊さない限りステップの実行順序の入れ替えが許容される。
【０１４３】
図２９は、本発明の実施の形態５に係る誤差推定部２６０１が行う処理の１周期のフローチャートである。フローチャート中の各ステップは、図２７に示す誤差推定部２６０１の各部が行う処理に対応する。ただし、本実施の形態５における誤差推定部２６０１の実行順序は、図２９に示したものに限定されない。各ステップの依存関係を壊さない限りステップの実行順序の入れ替えが許容される。
【０１４４】
本実施の形態５を、先の実施の形態２と組み合わせることも可能である。
【０１４５】
以上のように、実施の形態５によれば、パラメータ推定ベクトルｐ＾ｋを、一般化固有値問題を解くこと無く、先の実施の形態４と同様の効果を得ることができる。これにより、計算量が削減でき、組み込みシステムにおけるソフトウェアの実装に適する。
【符号の説明】
【０１４６】
１０１  蓄電池パラメータ推定装置、１０２  電流計測部、１０３  蓄電池、１０４  電流、１０５  電圧計測部、１０６  蓄電池電圧、１０７  ＳｏＣ推定部、１０８  ＳｏＣ推定値、１０９  第１係数計算部、１１０  第１係数、１１１  第２係数計算部、１１２  第２係数、１１３  電流積算部、１１４  電流積算値、１１５  誤差推定部、１１６  オフセット誤差、１１７  容量誤差、２０１  電圧源、２０２  電気抵抗、２０３  キャパシタンス、２０４  並列回路、２０５  電気抵抗、３０１  状態予測部、３０２  電流記憶部、３０３  状態記憶部、３０４  観測残差計算部、３０５  共分散予測部、３０６  共分散記憶部、３０７  観測行列計算部、３０８  観測残差共分散計算部、３０９  Ｋａｌｍａｎゲイン計算部、３１０  状態推定部、３１１  共分散推定部、３１２  ＳｏＣ計算部、４０１  予測状態微分係数計算部、４０２  推定状態微分係数記憶部、４０３  観測残差微分係数計算部、４０４  推定状態微分係数計算部、４０５  ＳｏＣ微分係数計算部、５０１  予測状態微分係数計算部、５０２  推定状態微分係数記憶部、５０３  観測行列微分係数計算部、５０４  測定残差微分係数計算部、５０５  予測共分散微分係数計算部、５０６  推定共分散微分係数記憶部、５０７  観測残差共分散微分係数計算部、５０８  Ｋａｌｍａｎゲイン微分係数計算部、５０９  推定状態微分係数計算部、５１０  推定共分散微分係数計算部、５１１  ＳｏＣ微分係数計算部、６０１  説明変数計算部、６０２  観測変数計算部、６０３  ＲＬＳゲイン計算部、６０４  逆分散記憶部、６０５  逆分散計算部、６０６  パラメータ推定部、６０７  パラメータ記憶部、６０８  オフセット計算部、６０９  容量計算部、１４０１  第３係数計算部、１４０２  第４係数計算部、１４０３  第５係数計算部、１４０４  誤差推定部、１６０１  誤差推定部、１７０１  説明変数計算部、１７０２  観測変数計算部、２２０１  誤差推定部、２３０１  共分散計算部、２３０２  固有ベクトル計算部、２３０３  パラメータ計算部、２３０４  共分散記憶部、２６０１  誤差推定部、２７０１  係数計算部、２７０２  ２次式求解部、２７０３  パラメータ計算部、２７０４  パラメータ記憶部。
【図１】
0006388563.tif 000076
【図２】
0006388563.tif 000077
【図３】
0006388563.tif 000078
【図４】
0006388563.tif 000079
【図５】
0006388563.tif 000080
【図６】
0006388563.tif 000081
【図７】
0006388563.tif 000082
【図８】
0006388563.tif 000083
【図９】
0006388563.tif 000084
【図１０】
0006388563.tif 000085
【図１１】
0006388563.tif 000086
【図１２】
0006388563.tif 000087
【図１３】
0006388563.tif 000088
【図１４】
0006388563.tif 000089
【図１５】
0006388563.tif 000090
【図１６】
0006388563.tif 000091
【図１７】
0006388563.tif 000092
【図１８】
0006388563.tif 000093
【図１９】
0006388563.tif 000094
【図２０】
0006388563.tif 000095
【図２１】
0006388563.tif 000096
【図２２】
0006388563.tif 000097
【図２３】
0006388563.tif 000098
【図２４】
0006388563.tif 000099
【図２５】
0006388563.tif 000100
【図２６】
0006388563.tif 000101
【図２７】
0006388563.tif 000102
【図２８】
0006388563.tif 000103
【図２９】
0006388563.tif 000104