(57)【要約】【課題】実車環境下において二次電池の等価回路モデルのＳＯＣ推定精度を向上できるようにした二次電池の等価回路モデルのパラメータ推定装置およびパラメータ推定方法を提供する。【解決手段】演算部は、等価回路モデルのパラメータのうちまず例えば抵抗Ｒ０の抵抗値を例えばナイキスト線図上などで周波数領域により算出する。そして、演算部は、この算出された抵抗Ｒ０の抵抗値を除く残りのパラメータを時間領域の例えば差分方程式などを用いて算出する。演算部が抵抗Ｒ０の抵抗値を周波数領域により算出するときには、予め定められた周波数範囲内における２つの周波数のインピーダンスＺ１、Ｚ２を算出し、これらのインピーダンスＺ１、Ｚ２を直線近似することで虚部Ｉｍがゼロとなる条件を満たす実部Ｒｅを抵抗Ｒ０の抵抗値として算出する。【選択図】図６
2015135277.tif 000002
【特許請求の範囲】
【請求項１】
ハイブリッド自動車又は電気自動車などの車両に搭載される二次電池（３）が等価回路モデル（９）に置換された状態で、前記二次電池（３）の検出電流及び検出電圧に応じて前記二次電池（３）の等価回路モデル（９）のパラメータを推定する推定部（８）を備え、
前記推定部（８）は、前記二次電池（３）が等価的に評価された等価回路モデル（９）のパラメータのうち一部の第１パラメータを周波数領域により算出した後、この算出された第１パラメータを除く残りの第２パラメータを時間領域により算出することを特徴とする二次電池の等価回路のパラメータ推定装置。
【請求項２】
請求項１記載の二次電池の等価回路のパラメータ推定装置において、
前記第１パラメータは第１抵抗（Ｒ０）を含み、
前記推定部（８）は、前記周波数領域において得られる複数のインピーダンス（Ｚ１、Ｚ２）を直線近似することで虚部がゼロとなる条件を満たす実部を前記第１抵抗（Ｒ０）の抵抗値として算出することを特徴とする二次電池の等価回路のパラメータ推定装置。
【請求項３】
ハイブリッド自動車又は電気自動車などの車両に搭載される二次電池（３）が等価回路モデル（９）に置換された状態で、前記二次電池（３）の検出電流及び検出電圧に応じて前記二次電池（３）の等価回路モデル（９）のパラメータを推定する推定部（８）を備え、
前記推定部（８）は、前記二次電池（３）が等価的に評価された等価回路モデル（９）のパラメータのうち一部の第１抵抗（Ｒ０）を含む第１パラメータについて、周波数領域において得られる複数のインピーダンス（Ｚ１、Ｚ２）を直線近似することで虚部がゼロとなる条件を満たす実部を前記第１抵抗（Ｒ０）の抵抗値として算出した後、この算出された第１抵抗（Ｒ０）を除く残りの第２パラメータを周波数領域において最小２乗法により算出することを特徴とする二次電池の等価回路のパラメータ推定装置。
【請求項４】
請求項１〜３の何れか一項に記載の二次電池の等価回路のパラメータ推定装置において、
前記等価回路モデル（９）は、第１抵抗（Ｒ０）と、前記第１抵抗（Ｒ０）に１又は複数個直列接続され、コンデンサ（Ｃ１、Ｃ２…Ｃｎ）及び第２抵抗（Ｒ１、Ｒ２…Ｒｎ）が並列接続されたＲＣ直並列回路（１２）と、を備えて構成され、
前記推定部（８）は、前記第１パラメータとして前記第１抵抗（Ｒ０）の抵抗値を算出することを特徴とする二次電池の等価回路のパラメータ推定装置。
【請求項５】
請求項２〜４の何れか一項に記載の二次電池の等価回路のパラメータ推定装置において、
前記推定部（８）は、前記第１抵抗（Ｒ０）の抵抗値を周波数領域により算出するときに、予め定められた周波数範囲内における第１及び第２周波数の第１及び第２インピーダンス（Ｚ１、Ｚ２）を算出し、前記第１及び第２インピーダンス（Ｚ１、Ｚ２）を直線近似することで虚部がゼロとなる条件を満たす実部を前記第１抵抗（Ｒ０）の抵抗値として算出することを特徴とする二次電池の等価回路のパラメータ推定装置。
【請求項６】
請求項２〜５の何れか一項に記載の二次電池の等価回路のパラメータ推定装置において、
前記推定部（８）は、前記第１抵抗（Ｒ０）の抵抗値を周波数領域により算出するための周波数として、前記二次電池（３）に接続された主機以外のコンポーネント（６）から発生され重畳された周波数成分の信号を用いて推定することを特徴とする二次電池の等価回路のパラメータ推定装置。
【請求項７】
ハイブリッド自動車又は電気自動車などに搭載される二次電池（３）が等価回路モデル（９）に置換された状態で、前記二次電池（３）の検出電流及び検出電圧に応じて前記二次電池の等価回路モデル（９）のパラメータを推定部（８）により推定する推定方法であって、
前記推定部（８）は、前記二次電池（３）が等価的に評価された等価回路モデル（９）のパラメータのうち一部の第１パラメータを周波数領域により算出し、
前記推定部（８）は、前記算出された第１パラメータを除く残りの第２パラメータを時間領域により算出することを特徴とする二次電池の等価回路のパラメータ推定方法。
【請求項８】
請求項７記載の二次電池の等価回路のパラメータ推定方法において、
前記第１パラメータは第１抵抗（Ｒ０）を含み、
前記推定部（８）は、前記周波数領域において得られる複数のインピーダンス（Ｚ１、Ｚ２）を直線近似することで虚部がゼロとなる条件を満たす実部を前記第１抵抗（Ｒ０）の抵抗値として算出することを特徴とする二次電池の等価回路のパラメータ推定方法。
【請求項９】
ハイブリッド自動車又は電気自動車などの車両に搭載される二次電池（３）が等価回路モデル（９）に置換された状態で、前記二次電池（３）の検出電流及び検出電圧に応じて前記二次電池（３）の等価回路モデル（９）のパラメータを推定部（８）により推定する推定方法であって、
前記推定部（８）は、前記二次電池（３）が等価的に評価された等価回路モデル（９）のパラメータのうち一部の第１抵抗（Ｒ０）を含む第１パラメータについて、周波数領域において得られる複数のインピーダンス（Ｚ１、Ｚ２）を直線近似することで虚部がゼロとなる条件を満たす実部を前記第１抵抗（Ｒ０）の抵抗値として算出し、
前記推定部（８）は、前記算出された第１抵抗（Ｒ０）を除く残りの第２パラメータを周波数領域において最小２乗法により算出することを特徴とする二次電池の等価回路のパラメータ推定方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、二次電池の等価回路のパラメータ推定装置及びパラメータ推定方法に関する。
【背景技術】
【０００２】
例えば、ハイブリッド自動車（ＨＶ）又は電気自動車（ＥＶ）などの車両には二次電池が搭載されている。このような二次電池の充電状態（ＳＯＣ：State of Charge）を逐次推定することは、バッテリパックの容量削減に大きく寄与するため重要である。
【０００３】
二次電池の等価回路モデルを使用する方法が二次電池の充電状態を推定するために用いられる（例えば特許文献１参照）。特許文献１記載の技術では、二次電池を等価回路モデルで模擬し、等価回路モデルの各パラメータを逐次最小２乗法により推定している。そして、その算出された各パラメータを用いて開回路電圧（ＯＣＶ：Open Circuit Voltage）を算出している。さらに、予め記録された開回路電圧（ＯＣＶ）と充電状態（ＳＯＣ）との関係に基づいて、その開回路電圧（ＯＣＶ）に対応する充電状態の値（ＳＯＣ）を算出している。一般に、ＳＯＣを推定するため等価回路モデルを用いた方法が知られている。
【０００４】
他方、周波数応答分析（ＦＲＡ：Frequency Response Analyzer）の計測分野においては、インピーダンスのナイキスト線図を描くことで等価回路モデル中の各パラメータを推定する方法が知られている（例えば、特許文献２参照）。この特許文献２記載の技術では、電流または電圧を階段状または矩形波状に変化させたときの各時間の電圧または電流の変化を用いてインピーダンスの周波数特性を測定し、ナイキスト線図を描いている。
【先行技術文献】
【特許文献】
【０００５】
【特許文献１】特開２０１１−１２２９５１号公報
【特許文献２】特開２００６−１６２２８３号公報
【発明の概要】
【発明が解決しようとする課題】
【０００６】
特許文献１記載の技術を用いてＲＣ並列回路をｎ個直列に接続したＲＣ直並列回路を等価回路モデルとすると、この等価回路モデルにおいてＳＯＣ推定精度を向上するためには直列接続個数を増加する、すなわちｎを大きくすることが考えられるが、高次の計算を必要とするため演算量が大幅に増加してしまう。例えば、ハイブリッド自動車（ＨＶ）又は電気自動車（ＥＶ）などの車両バッテリ制御技術にこの演算方法を用いようとしても、例えば演算装置の処理能力及び搭載メモリ量に制約があることが多いため、実写環境下において充分な精度でＳＯＣを推定することが困難である。
【０００７】
特許文献２記載の技術を用いてインピーダンスのナイキスト線図に基づき等価回路モデルのパラメータを推定したとしても、ナイキスト線図の全体を精度良く描こうとした場合に必要な周波数成分数が多くなり演算量が膨大になる。また、ＦＲＡの計測器は、ユーザが指示できれば必要な周波数成分を含んだ信号を発生させることができるものの、例えば前述の車両内環境下では二次電池の電流、電圧は主に主機モータの電流に応じて変化するため、周波数成分がユーザによる運転パターンに大きく依存する。このためＦＲＡ計測器に周波数成分を指示したとしても実車環境に応じた周波数成分を網羅していない場合も考えられる。また電流、電圧を矩形波状に変化させることも困難である。したがって、ナイキスト線図全体を精度良く描くことが困難である。
【０００８】
さらに高周波数領域では、二次電池の等価回路モデルの容量成分が０に近づき、ほぼ抵抗成分Ｒ０のみとなりやすいが、この高周波領域ではインピーダンスの絶対値が小さくなってしまい、電流変化又は／及び電圧変化が生じにくい。このため、実車環境下ではノイズの影響が大きく測定が困難である。
【０００９】
本発明の目的は、推定演算量を大幅に削減し、実車環境下においても二次電池の等価回路モデルのＳＯＣを精度良く推定できるようにした二次電池の等価回路モデルのパラメータ推定装置およびパラメータ推定方法を提供することにある。
【課題を解決するための手段】
【００１０】
請求項１記載の発明によれば、推定部は、二次電池が等価的に評価された等価回路モデルのパラメータのうち一部の第１パラメータを周波数領域により算出した後、この算出された第１パラメータを除く残りの第２パラメータを時間領域により算出するため、大幅に推定演算量を削減できる。これにより、実車環境下の限られたリソースの中でも、精度良くＳＯＣを推定できる。
【００１１】
請求項２記載の発明によれば、推定部は、等価回路モデルのパラメータのうち第１抵抗の抵抗値を周波数領域により算出した後、この算出された第１抵抗以外の第２パラメータを時間領域により算出するため、大幅に推定演算量を削減できる。これにより、実車環境下の限られたリソースの中でも、精度良くＳＯＣを推定できる。
【００１２】
請求項３記載の発明によれば、推定部は、等価回路モデルのパラメータのうち第１抵抗の抵抗値を周波数領域により算出した後、この算出された第１抵抗以外の第２パラメータを周波数領域において最小２乗法を用いて算出するため、大幅に推定演算量を削減できる。これにより、実車環境下の限られたリソースの中でも、精度良くＳＯＣを推定できる。
【図面の簡単な説明】
【００１３】
【図１】本発明の一実施形態に係るシステムブロック構成例を概略的に示す電気的構成図
【図２】二次電池の等価回路モデルの一例
【図３】二次電池の充電状態の推定処理の流れの一例を概略的に示すフローチャート
【図４】二次電池の等価回路モデルにおいて二次成分まで考慮した場合のインピーダンスのナイキスト線図の一例
【図５】第１抵抗のパラメータを推定するための処理内容の一例を概略的に示すフローチャート
【図６】（ａ）〜（ｄ）は処理内容の一例をナイキスト線図上で示す説明図
【発明を実施するための形態】
【００１４】
以下、本発明の一実施形態を説明する。図１は、本実施形態に係る二次電池の等価回路パラメータの推定装置１の電気的構成を概略的なブロック図により示す。
車両内には主機モータによる主機負荷２が搭載されている。この負荷２には二次電池３のセル群から電力供給されている。電圧値測定部４が二次電池３のセルの電圧値を測定するために設けられており、温度測定部５が二次電池３の温度を測定するために設けられる。主機モータによる負荷２のほかには他の負荷（他負荷）６が接続されており、この他負荷６にも二次電池３のセル群から電力供給されている。
【００１５】
他負荷６としては、例えばＤＣＤＣコンバータ、エアコンディショナのコンプレッサ、インバータなどが挙げられる。電流値測定部７は、二次電池３の電流値を測定するために設けられる。
【００１６】
演算部８は、例えばマイクロコンピュータを用いて構成され、この演算部８には、電圧値測定部４の電圧値測定情報、温度測定部５の温度測定情報、及び、電流値測定部７の電流値測定情報が与えられる。演算部８は、これらの測定情報の少なくとも一部以上を用いて二次電池３の等価回路のパラメータの値を演算する。
【００１７】
二次電池３は、図２に示す等価回路モデル９に模擬できる。この図２に示す二次電池３の等価回路モデル９は、二次電池３の純電圧要素１０と内部インピーダンス要素１１とを端子間に備える。この内部インピーダンス要素１１は、抵抗（第１抵抗相当）Ｒ０と、この抵抗Ｒ０に複数個直列接続されコンデンサＣ１…Ｃｎ及び抵抗Ｒ１…Ｒｎがそれぞれ並列接続されたＲＣ直並列回路１２と、を備える。抵抗Ｒ０は０次成分となる。
【００１８】
なお、ＲＣ直並列回路１２は、そのＲＣ並列回路の直列接続個数を１個としても良い。ここで、電流値測定部７は、二次電池３から負荷２及び６に流れる電流値Ｉを測定することで二次電池３に流れる電流を測定でき、電圧値測定部４は、二次電池３の両端子間の電圧値を測定できる。
【００１９】
この図２において、ＶＯＣＶ（ＯＣＶ：Open Circuit Voltage）が開回路電圧を示しており、Ｖｚを過電圧と定義すると、開回路電圧ＶＯＣＶは、電気化学的並平衡状態における電極間の電位差を示しており、過電圧Ｖｚは二次電池３に電流が流れたときの二次電池３の内部インピーダンス要素１１による電圧降下を示している。なお、厳密に言えば、インダクタ成分なども存在するが、例えば無視可能な程度となるため図２には図示していない。二次電池３の充電状態（ＳＯＣ）を知るため、このような等価回路モデル９を使用する。
【００２０】
図３は二次電池の充電状態を推定するための演算部による処理を概略的に示すフローチャートである。演算部８は、二次電池３の電流値ｉ（ｔ）、電圧値ｖ（ｔ）、温度Ｔ（ｔ）の情報をそれぞれ電流値測定部７、電圧値測定部４、温度測定部５から取得する（Ｓ１）。演算部８は、まず一次成分となる抵抗Ｒ０のパラメータ（抵抗値）を推定する（Ｓ２）。このステップＳ２の詳細は後述する。演算部８は、この算出された抵抗Ｒ０のパラメータを用いて、その他の抵抗Ｒ１…Ｒｎ及びコンデンサＣ１…Ｃｎのパラメータ（各インピーダンス値）を推定する（Ｓ３）。また、このステップＳ３の方法は、時間領域と周波数領域における演算処理を組み合わせた演算処理方法、または、周波数領域の演算処理方法となるが、その詳細については、演算例１、演算例２を挙げてそれぞれ後述する。
【００２１】
そして、演算部８は、これらのパラメータ（各インピーダンス値）に応じて開回路電圧ＶＯＣＶを演算する（Ｓ４）。さらに、演算部８は、ステップＳ４において算出された開回路電圧ＶＯＣＶ、事前に算出されていた開回路電圧ＶＯＣＶ−充電状態（ＳＯＣ）のマップ、温度測定部５の温度情報、に応じて二次電池３の充電状態（ＳＯＣ）を推定する（Ｓ５）。演算部８は、これらのステップＳ１〜Ｓ５の処理を一定周期Ｔ０毎に実施する。
【００２２】
図４は等価回路モデル９において二次成分まで考慮した場合のインピーダンスのナイキスト線図を示している。また図５は、演算部８が図３のステップＳ２において抵抗Ｒ０のパラメータ（抵抗値）を推定するための処理内容をフローチャートにより概略的に示している。この図５に示す処理は実車環境にて行なわれる。また図６（ａ）〜図６（ｄ）は処理内容をナイキスト線図で表している。
【００２３】
図４に示すインピーダンスのナイキスト線図において、抵抗Ｒ０付近では実部Ｒｅがほぼ一定となる。すなわち、高周波数領域において虚部Ｉｍのみが変化する特性となる。なお、インダクタンス成分が含まれるため、厳密には高周波数領域において虚部Ｉｍが正になることもある。そして、低周波数になるに従い、コンデンサＣ１、Ｃ２の一次成分、二次成分の影響が大きくなり、図示破線に沿うような特性となる。ただし、この図４に示す破線は理論的に考慮した場合の理想的な波形を示すが、実測定上では図４の実線に示すように破線の一部に沿った特性となる。
【００２４】
二次電池３およびその周辺回路の設計者は、特に虚部Ｉｍの成分が０に近づく周波数範囲を実験的手法などで予め算出できる。すなわち、設計者は抵抗Ｒ０を跨いで虚部Ｉｍが変化する周波数ｆ０を含む周波数範囲を経験的に求めておくことができる。この周波数ｆ０は例えば数百Ｈｚ〜数ｋＨｚの範囲内に予め設定することができる。
【００２５】
実車環境下では、図５に示すように、演算部８は、この予め求められた周波数ｆ０の第１のインピーダンスＺ１を、例えばウェーブレット変換を用いて周波数領域で算出する（Ｓ１１）。なお、ここでインピーダンスＺ１を周波数領域で算出する手法はウェーブレット変換処理に限られるものではない。
【００２６】
次に、演算部８は、周波数ｆ０におけるインピーダンスＺ１の虚部Ｉｍに注目し、インピーダンスＺ１の虚部Ｉｍが０であれば（Ｓ１２：ＹＥＳ）、抵抗Ｒ０の抵抗値を、算出された第１のインピーダンスＺ１の実部Ｒｅに設定する（Ｓ１３）。逆に、演算部８は、第１のインピーダンスＺ１の虚部Ｉｍが０でなく負であれば（Ｓ１２：ＮＯ、Ｓ１４：ＹＥＳ）、周波数ｆ０より高周波数ｆ０＋Δｆ（但しΔｆ＞０）において第２のインピーダンスＺ２を算出する（Ｓ１５）。この第２のインピーダンスＺ２はウェーブレット変換処理の手法を用いて算出できる。なお、ここでもインピーダンスＺ２を周波数領域で算出する手法はウェーブレット変換処理に限られるものではない。ここで、Δｆは数百Ｈｚ程度の値を採用することができる。
【００２７】
また、演算部８は、インピーダンスＺ１が０でなく正であれば（Ｓ１２：ＮＯ、Ｓ１４：ＮＯ）、周波数ｆ０より低周波数ｆ０−Δｆ（但しΔｆ＞０）において第２のインピーダンスＺ２を算出する（Ｓ１６）。
【００２８】
例えば、図６（ａ）に示すように、演算部８は、第１のインピーダンスＺ１の虚部Ｉｍが負であれば、高周波数ｆ０＋Δｆにおいて第２のインピーダンスＺ２を算出する。高周波数ｆ０＋Δｆで求めれば、第２のインピーダンスＺ２は第１のインピーダンスＺ１に比較して虚部Ｉｍが正方向に遷移することになり、第２のインピーダンスＺ２はナイキスト線図上の実軸との交点方向に動く。
【００２９】
また、図６（ｂ）に示すように、演算部８は、第１のインピーダンスＺ１の虚部Ｉｍが正であれば、その低周波数ｆ０−Δｆにおいて第２のインピーダンスＺ２を算出する。低周波数ｆ０−Δｆで求めれば、第２のインピーダンスＺ２は第１のインピーダンスＺ１に比較してその虚部Ｉｍが負方向に遷移することになり、第２のインピーダンスＺ２はナイキスト線図上で実軸との交点方向に動く。
【００３０】
そして、演算部８は、これらの第１及び第２のインピーダンスＺ１及びＺ２の虚部Ｉｍに注目し、これらの虚部Ｉｍの符号が同じであるか判定する（Ｓ１７）。演算部８は、虚部Ｉｍの符号が同一であれば（Ｓ１７：ＹＥＳ）、これらの第１及び第２のインピーダンスＺ１及びＺ２を例えば直線近似処理して外方に挿引（外挿）し（Ｓ１８）、ナイキスト線図の実軸との交点を求め虚部Ｉｍが０となる実部Ｒｅを算出する。このとき、この算出された実部Ｒｅを抵抗Ｒ０の抵抗値とする（Ｓ２０）。
【００３１】
例えば、図６（ｃ）に示すように、第１及び第２のインピーダンスＺ１及びＺ２の虚部Ｉｍが共に負であるときには、これらのインピーダンスＺ１及びＺ２を外挿することでナイキスト線図の実軸との交点を求めることで抵抗Ｒ０の抵抗値を求めることができる。これは、第１及び第２のインピーダンスＺ１及びＺ２の虚部Ｉｍが共に正であるときも同様となる。
【００３２】
逆に、演算部８は、ステップＳ１７において虚部Ｉｍの符号が互いに異なれば（Ｓ１７：ＮＯ）、第１及び第２のインピーダンスＺ１及びＺ２を例えば直線近似処理して内方に挿引（内挿）し（Ｓ１９）、ナイキスト線図の実軸との交点を求め虚部Ｉｍが０となる実部Ｒｅを算出する（Ｓ２０）。
【００３３】
例えば、図６（ｄ）に示すように、第１のインピーダンスＺ１の虚部Ｉｍが負であり、第２のインピーダンスＺ２の虚部Ｉｍが正であるときには、これらのインピーダンスＺ１及びＺ２を内挿することでナイキスト線図の実軸との交点を求め、抵抗Ｒ０の抵抗値を求めることができる。これは、第１のインピーダンスＺ１の虚部Ｉｍが正であり、第２のインピーダンスＺ２の虚部Ｉｍが負であるときも同様となる。
【００３４】
以下、図３に示したステップＳ３のパラメータ推定方法の詳細例（演算例１、演算例２）を挙げ、図３のステップＳ２に示したように抵抗Ｒ０の抵抗値（インピーダンスＺ０）を予め算出した場合と、抵抗Ｒ０の抵抗値を予め算出しない場合とで、その演算量を比較する。
【００３５】
＜演算例１：ステップＳ３に時間領域の差分方程式を用いた場合の演算例＞
図２の等価回路モデルにおいてｎ＝２とすると、インピーダンスＺ（ｓ）は下記の（１）式のような伝達関数により求めることができる。
【００３６】
【数１】
2015135277.tif 000003
このインピーダンスＺ（ｓ）は、その分子が２次式であり分母も２次式となる。この（１）式では、分子の２次、１次、０次の各係数をａ、ｂ、ｃとし、分母の２次、１次、０次の各係数を「１」、ｄ、ｅとしている。またこのインピーダンスＺ（ｓ）に流れる電流をＩ（ｓ）とすると、このインピーダンスＺ（ｓ）にかかる電圧Ｖ（ｓ）は、下記の（２）式のようになる。
【００３７】
【数２】
2015135277.tif 000004
また、Ｔをサンプリング周期として下記の（３）式のように定義する。
【００３８】
【数３】
2015135277.tif 000005
この（３）式を用いて電圧Ｖ（ｓ）を離散化すると、下記の（４）式のように求めることができる。
【００３９】
【数４】
2015135277.tif 000006
この（４）式を時間領域の差分方程式によって表わすと、下記の（５）式のようになる。  
【００４０】
【数５】
2015135277.tif 000007
ここで、最小二乗法を用いて計測値にフィットする関数の係数を求めることを考える。最小二乗法は、測定値にフィットするように関数の係数を求める手法であり、測定値と関数の残差の二乗和を最小とするように係数を決定する。例えば、関数の残差をｅとすれば、下記の（６）式のように残差ｅを定義できる。
【００４１】
【数６】
2015135277.tif 000008
また、このとき残差ｅの二乗ノルムは、下記の（７）式のように計算できる。
【００４２】
【数７】
2015135277.tif 000009
この（７）式が最小となる条件は、下記の（８）式を満たす場合であり、この条件を満たすｘの値は（９）式のように求められる。すなわち、いわゆる擬似逆行列（（９）式のｙの係数部）を求めることで、最小二乗法による残差ｅの二乗ノルムが最小となる条件を満たす係数を求めることができる。
【００４３】
【数８】
2015135277.tif 000010
そこで、前述の（５）式を行列式に展開すると、下記の（１０）式のように求められる。
【００４４】
【数９】
2015135277.tif 000011
最小二乗法による残差ｅの二乗ノルムが最小となる条件を満たすように（１０）式を変換すると、（１１）式のように求められる。
【００４５】
【数１０】
2015135277.tif 000012
抵抗Ｒ０の抵抗値が予め求められていない場合、ｘが５×１の行列式となるため、この行列式を解くために最低５つの方程式が必要となる。この場合、７サンプリング分の電圧測定値Ｖ、電流測定値Ｉを必要とする。
【００４６】
しかし、抵抗Ｒ０の抵抗値が予め求められていると、（１）式に代わるインピーダンスＺ（ｓ）の伝達関数は、下記の（１２）式のようにＲ０の項がなくなる。
【００４７】
【数１１】
2015135277.tif 000013
このインピーダンスＺ（ｓ）は、その分子が１次式であり分母も２次式となる。この（１２）式では、分子の１次、０次の各係数をｆ、ｇとし、分母の２次、１次、０次の各係数を「１」、ｈ、ｉとしている。（５）式と同様に（１２）式を時間領域の差分方程式によって表現すると、下記の（１３）式のように表すことができる。
【００４８】
【数１２】
2015135277.tif 000014
この（１３）式を行列式に展開すると、下記の（１４）式のように求められる。
【００４９】
【数１３】
2015135277.tif 000015
最小二乗法による残差ｅの二乗ノルムが最小となる条件を満たすように（１４）式を変換すると、下記の（１５）式のように求められる。
【００５０】
【数１４】
2015135277.tif 000016
抵抗Ｒ０の抵抗値が予め求められていると、ｘが４×１の行列式となるため、この行列式を解くため４つの方程式で解くことができる。したがって６回分の電圧測定値Ｖ、電流測定値Ｉのサンプリング値のみで時間領域の差分方程式を用いて解くことができる。未知数が一つ減少すると、必要な連立方程式の数が減少することになり、行列式の次数が下がる。すると、必要な電流測定値Ｉ、電圧測定値Ｖのサンプリング数を減少させることができる。
【００５１】
以下、抵抗Ｒ０の抵抗値を予め算出しない場合と、抵抗Ｒ０の抵抗値を予め算出した場合とについて、擬似逆行列（ｙの係数）の計算量を詳細に考察する。
「Ａ」がＮ行Ｍ列とすると、
（Ｘ１）Ｎ行Ｍ列の行列「Ａ」について、（ｉ，ｊ）要素と（ｊ，ｉ）要素を入れ替えたＭ行Ｎ列の転置行列「ＡＴ」を作成する。
（Ｘ２）転置行列「ＡＴ」と、行列「Ａ」とを乗算する。
（Ｘ３）乗算結果「ＡＴＡ」の逆行列を計算する。
（Ｘ４）逆行列（ＡＴＡ）−１と転置行列「ＡＴ」を乗算する。
【００５２】
これらの（Ｘ１）〜（Ｘ４）の処理が必要となる。（Ｘ１）では行列の入替処理となり加減乗除処理は必要としない。（Ｘ２）ではＮ×Ｍ×Ｍ回の乗算と（Ｎ−１）×Ｍ×Ｍ回の加減算処理が必要となる。（Ｘ３）ではＭ行×Ｍ列の行列の逆行列を計算するため、ガウスの消去法（掃き出し法）を用いるとＯ（Ｍ＾３）の処理が必要となる。（Ｘ４）ではＮ×Ｍ×Ｍ回の乗算処理とＮ×Ｍ×（Ｍ−１）回の加減算処理が必要となる。
【００５３】
例えばｎ＝２とすると、
抵抗Ｒ０の抵抗値を予め求めない場合にはＭ＝６となる。このため、（Ｘ２）では３６Ｎ回の乗算処理と、３６（Ｎ−１）回の加減算処理を必要とする。（Ｘ３）では数百回の処理（例えば２５２回の乗算処理と２１０回の加減算処理）を必要とする。（Ｘ４）では数十回の処理（例えば３６Ｎ回の乗算処理と３０Ｎ回の加減算処理）を必要とする。（Ｘ２）〜（Ｘ４）を合計すると、乗算処理が例えば７２Ｎ＋２５２回、加減算処理が例えば６６Ｎ＋１７４回となる。
【００５４】
逆に、抵抗Ｒ０の抵抗値を予め求めている場合にはＭ＝４となる。このため、（Ｘ２）では１６Ｎ回の乗算と１６（Ｎ−１）回の加減算処理を必要とする。（Ｘ３）ではガウスの消去法を用いると数十回（例えば８０回の乗算処理と６０回の加減算処理）で良い。（Ｘ４）では十数回（例えば１６Ｎ回の乗算処理と１２Ｎ回の加減算処理）で良い。これらの（Ｘ２）〜（Ｘ４）を合計しても、乗算処理が例えば３２Ｎ＋８０回、加減算処理が例えば２８Ｎ＋４４回となる。これにより、抵抗Ｒ０の抵抗値が予め算出されていれば、大幅に処理回数を削減できることがわかる。
【００５５】
なお、抵抗Ｒ０の抵抗値を予め算出する際には、前述したように周波数領域において直線近似する必要があるが、この近似処理に必要な処理数は、周波数成分数×（加減算処理１回+乗算処理１回）×２（実部と虚部）となり、抵抗Ｒ０の抵抗値を予め算出しない場合に比較してこの分だけ演算量が増加する。しかしながら、この演算量は、前述の処理削減量に比較すれば微々たるものであり、この増加分を考慮したとしても減少分の寄与度が大きく、大幅に処理回数を削減できる。
【００５６】
また、抵抗Ｒ０の抵抗値が予め算出されている場合には、これらの処理時間は必要な電流測定値Ｉ、電圧測定値Ｖのサンプリング数の減少分だけ減少することになる。これにより、さらに処理時間を低減できることになる。
【００５７】
（演算例２：他の演算例）
ｓ→ｊωとした他の演算例について演算量を考慮する。図２の等価回路モデルにおいてｎ＝２とすると、インピーダンスＺ（ｊω）は下記の（１６）式のように表すことができる。
【００５８】
【数１５】
2015135277.tif 000017
ここで、ａ＝１／（Ｒ１×Ｃ１）、ｂ＝１／（Ｒ２×Ｃ２）、ｃ＝１／Ｃ１、ｄ＝１／Ｃ２、ｅ＝Ｒ０、とすると、インピーダンスＺ（ｊω）は下記の（１７）式のように表すことができる。
【００５９】
【数１６】
2015135277.tif 000018
このインピーダンスＺ（ｊω）はその分子が２次式であり分母も２次式となる。この（１７）式を行列式に展開すると、下記の（１８）式のように求められる。
【００６０】
【数１７】
2015135277.tif 000019
前述したように、最小二乗法による残差ｅの二乗ノルムが最小となる条件を満たすように（１８）式を変換すると、下記の（１９）式のように求められる。
【００６１】
【数１８】
2015135277.tif 000020
抵抗Ｒ０の抵抗値が予め求められていない場合、ｘが６×１の行列式となるため、この行列式を解くために最低でも３つの周波数においてインピーダンスＺ（ｊω）を推定しなければならない。したがって、行列Ａは最低６×６行列となる。
【００６２】
他方、抵抗Ｒ０の抵抗値が予め求められていると、（１６）式に代わるインピーダンスＺ（ｊω）の伝達関数は、下記の（２０）式のようにＲ０の項がなくなる。
【００６３】
【数１９】
2015135277.tif 000021
ここで、ａ＝１／（Ｒ１×Ｃ１）、ｂ＝１／（Ｒ２×Ｃ２）、ｃ＝１／Ｃ１、ｄ＝１／Ｃ２、とすると、インピーダンスＺ（ｊω）は下記の（２１）式のように表すことができる。
【００６４】
【数２０】
2015135277.tif 000022
このインピーダンスＺ（ｊω）はその分子が１次式であり、分母が２次式となる。この（２１）式を行列式に展開すると、下記の（２２）式のように求められる。
【００６５】
【数２１】
2015135277.tif 000023
さらに前述したように、最小二乗法による残差ｅの二乗ノルムが最小となる条件を満たすように（２２）式を変換すると、下記の（２３）式のように求められる。
【００６６】
【数２２】
2015135277.tif 000024
抵抗Ｒ０の抵抗値が予め求められている場合、ｘが４×１の行列式となるため、この行列式を解くためには、２つの周波数においてインピーダンスＺ（ｊω）を推定すれば良くなり、行列Ａの行数、列数は４×４行列で良い。
【００６７】
すなわち、抵抗Ｒ０の抵抗値が予め求められていると、伝達関数の分子の次数を低減でき、未知数が減少することで行列式の次数が下がり、必要な周波数成分も下がる。抵抗Ｒ０の抵抗値を予め算出せず例えば周波数領域で全てのパラメータを算出する場合と比較して演算量を削減できる。このような計算手法を用いても演算量を大幅に削減できる。
【００６８】
＜まとめ＞
本実施形態によれば、演算部８は、等価回路モデル９のパラメータのうちまず例えば抵抗Ｒ０の抵抗値を例えばナイキスト線図上などで周波数領域により算出した後、この算出された抵抗Ｒ０の抵抗値を除く残りの抵抗Ｒ１…Ｒｎの抵抗値、コンデンサＣ１…Ｃｎの容量値を時間領域の例えば差分方程式などを用いて算出している。すると、差分方程式を用いて算出するときに未知数が１つ減少するため演算量を大幅に削減できる。また、等価回路モデル９の全てのパラメータを時間領域で算出するよりも演算量を削減できる。これにより、実車環境下の限られたリソースの中でも、精度良くＳＯＣを推定できる。
【００６９】
また、演算部８が抵抗Ｒ０の抵抗値を周波数領域により算出するときには、予め定められた周波数範囲内における２つの周波数のインピーダンスＺ１、Ｚ２を算出し、これらのインピーダンスＺ１、Ｚ２を直線近似することで虚部Ｉｍがゼロとなる条件を満たす実部Ｒｅを抵抗Ｒ０の抵抗値として算出している。このため、演算量を極力少なくしながら抵抗Ｒ０の抵抗値を極力正確に求めることができる。全てのパラメータを推定するために、全ての周波数領域のナイキスト線図を描く必要がなくなり、抵抗Ｒ０付近のインピーダンスのナイキスト線図を描けばよいため、実車環境下の限られたリソースの中でも精度良く抵抗Ｒ０を容易に推定することができる。
【００７０】
また、別の例では、演算部８は、等価回路モデル９のパラメータのうちまず例えば抵抗Ｒ０の抵抗値を例えばナイキスト線図上などで周波数領域により算出した後、この算出された抵抗Ｒ０の抵抗値を除く残りの抵抗Ｒ１…Ｒｎの抵抗値、コンデンサＣ１…Ｃｎの容量値を周波数領域において最小２乗法を用いて算出している。すると、未知数が１つ減少するため演算量を大幅に削減できる。
【００７１】
（適用例）
図１において、主機負荷２は車両内の主機モータを示すが、この主機モータの駆動制御パターンによる変化だけでは、電流測定値Ｉ、電圧測定値Ｖに数ｋＨｚ以上の高周波成分は生じにくい。そこで、演算部８は、二次電池３以外に接続されている主機負荷２以外の他負荷６から発生され重畳された周波数成分の信号を用いて抵抗Ｒ０の抵抗値を推定しても良い。例えば、演算部８は、主機負荷２以外の負荷６から発生され重畳された周波数成分の一部又は全部におけるインピーダンスを前述の第１インピーダンスＺ１、第２インピーダンスＺ２とし、内挿又は外挿して虚部Ｉｍが０となる実部Ｒｅを求めることで抵抗Ｒ０の抵抗値を推定しても良い。前述したように、主機モータ以外の他負荷６としては、例えばＤＣＤＣコンバータ（スイッチングレギュレータ）、エアコンディショナのコンプレッサ、インバータなどが挙げられる。
【００７２】
ＤＣＤＣコンバータは例えば１００ｋＨｚ、エアコンディショナのコンプレッサは例えば２０ｋＨｚ、インバータは例えば５〜１０ｋＨｚのスイッチング周波数成分をそれぞれ発生するため、これらの周波数成分が電流値測定部７の電流測定値Ｉ、電圧値測定部４の電圧測定値Ｖにも含まれることになる。演算部８は、例えばウェーブレット変換処理などを用い、各他負荷６の発生成分がどのような周波数成分を含んでいるかを判定でき、さらに前述のインピーダンスＺ１、Ｚ２に対応するインピーダンスを算出できる。そして、演算部８は、例えば直線近似することで抵抗Ｒ０の抵抗値を精度良く求めることができる。
【００７３】
本適用例によれば、演算部８は他負荷６から発生され重畳された周波数成分の信号を用いて抵抗Ｒ０を推定するため、抵抗Ｒ０の抵抗値を精度良く求めることができ、Ｓ／Ｎの高いインピーダンス計測処理ができる。
【００７４】
（他の実施形態）
本発明は、前述実施形態に限定されるものではなく、例えば、以下に示す変形または拡張が可能である。
等価回路モデル９は、抵抗Ｒ０と、抵抗Ｒ０に直列接続されたＲＣ直並列回路１２とを備えた形態を示したが、等価回路モデル９はこの形態に限られるものではない。例えば、インダクタンス成分を含んで構成することもできる。
【００７５】
第１パラメータとして抵抗Ｒ０の抵抗値、第２パラメータとしてそれ以外の抵抗Ｒ１…Ｒｎ、コンデンサＣ１…Ｃｎの値に適用した実施形態を示したが、これに限定されるものではない。例えば、抵抗Ｒ０と共に抵抗Ｒ１の抵抗値を第１パラメータとして加入しても良いし、さらに抵抗Ｒ２の抵抗値などを第１パラメータとして加入しても良い。
【００７６】
なお、特許請求の範囲に付した括弧付き符号は本願明細書の構成要素に対応する符号を付したものであり構成要素の一例を挙げたものである。したがって、本願に係る発明は当該特許請求の範囲の構成要素に付した符号の要素に限られるわけではなく、特許請求の範囲内の用語又はその均等の範囲で様々な拡張が可能である。
【符号の説明】
【００７７】
図面中、２は主機による負荷、３は二次電池、６は他負荷（主機以外のコンポーネント）、８は演算部（推定部）、９は等価回路モデル、を示す。
【図１】
2015135277.tif 000025
【図２】
2015135277.tif 000026
【図３】
2015135277.tif 000027
【図４】
2015135277.tif 000028
【図５】
2015135277.tif 000029
【図６】
2015135277.tif 000030