(57)【要約】  【課題】電気化学反応系のインピーダンス・スペクトルから、化学的に意味のある等価回路を自動的に生成する。  【解決手段】インピーダンス・スペクトルを対数周波数軸上で等間隔に区分し、区分の端点のインピーダンスから夫々Ｒ／／Ｃ、Ｒ／／Ｌ等の等価回路要素を先ず割り当てる。これらの要素を直列接続した初期近似等価回路を、各要素中の抵抗成分を正値に拘束されたパラメーターとして実際のインピーダンス・スペクトルに対してフィッティングする。次に、各要素中の全ての成分を拘束条件なしでフィッティングし、その結果から、パラメーターが負になった要素を削除することなどにより、余分な要素の削除を行う。その後、全ての成分を正値に拘束して再度フィッティングを行う。【選択図】なし【特許請求の範囲】
【請求項１】
以下のステップ（ア）〜（コ）を設け、所定の周波数範囲にわたって測定された電気化学反応系のインピーダンスから、コイルからなる等価回路要素、キャパシターからなる等価回路要素、抵抗からなる等価回路要素、抵抗とコイルとの並列接続からなる等価回路要素、及び抵抗とキャパシターとの並列接続からなる等回路要素からなる群から選択された少なくとも一の等価回路要素の直列接続として構成される前記電気化学反応系の等価回路を推定する方法。
（ア）前記周波数範囲上の複数の周波数で測定されたインピーダンスの測定値の系列を受け取る。
（イ）周波数の順に並べられた前記インピーダンスの測定値の系列を複数の周波数区分に分割するとともに、前記周波数区分の各々から一つのインピーダンスを選択する。
（ウ）前記周波数区分毎に、前記選択された一つのインピーダンス及びそのインピーダンスが測定された周波数に基づいて、前記群から選択された一つの等価回路要素を対応付け、前記対応付けられた等価回路要素からなる等価回路要素リストを得る。
（エ）ステップ（ウ）で得られた前記等価回路要素リスト中の前記等価回路要素を直列接続した回路のインピーダンスが前記測定されたインピーダンスの系列にフィッティングされるように、前記等価回路要素リスト中の等価回路要素に含まれる抵抗の抵抗値を正の値の範囲で最適化する。
（オ）ステップ（エ）によって最適化された等価回路要素リスト中で、隣接する前記周波数区分に対応する２つの等価回路要素の何れかが前記抵抗とコイルとの並列接続からなる等価回路要素であり、かつ周波数が高い方の前記周波数区分から選択された前記インピーダンスの実数項の値が周波数が低い方の前記周波数区分から選択された前記インピーダンスの実数項の値以下である場合には、前記周波数が高い方の周波数区分に対応する前記等価回路要素を前記コイルからなる等価回路要素に変更する。
（カ）前記等価回路要素リスト中に含まれる前記等価回路要素に対応する前記測定されたインピーダンスの実数項から前記等価回路要素のインピーダンスの計算値を減算し、前記減算の結果を前記等価回路要素リスト中の全ての前記等価回路要素について平均して平均実数項差を得る。
（キ）前記平均実数項差の値が正の場合には、前記平均実数項差の値の抵抗値を有する抵抗からなる等価回路要素を前記等価回路要素リストに追加する。
（ク）ステップ（キ）で得られた前記等価回路要素リスト中の前記等価回路要素を直列接続した回路のインピーダンスが前記測定されたインピーダンスの系列にフィッティングされるように、前記等価回路要素リスト中の等価回路要素に含まれる抵抗の抵抗値、コイルのインダクタンス及びキャパシタのキャパシタンスを、前記抵抗値、インダクタンス及びキャパシタンスの値の正負に拘束条件を付けずに最適化する。
（ケ）ステップ（ク）の最適化によって、抵抗の抵抗値、コイルのインダクタンス及びキャパシタンスのうちの少なくとも一つが負の値になった前記等価回路要素を前記等価回路要素リストから除去する。
（コ）ステップ（ケ）によって得られた前記等価回路要素リスト中の前記等価回路要素を直列接続した回路のインピーダンスが前記測定されたインピーダンスの系列にフィッティングされるように、前記等価回路要素リスト中の等価回路要素に含まれる抵抗の抵抗値、コイルのインダクタンス及びキャパシタのキャパシタンスを前記抵抗値、インダクタンス及びキャパシタンスが全て正であるとの拘束条件を付けて最適化する。
【請求項２】
前記複数の周波数区分への分割は対数周波数軸上で等間隔となるように行われる、請求項１に記載の等価回路を推定する方法。
【請求項３】
前記インピーダンスの系列は、対数周波数軸上で等間隔で配置された複数の周波数について測定されたインピーダンスの系列である、請求項１または２に記載の等価回路を推定する方法。
【請求項４】
ステップ（イ）における前記選択された一つのインピーダンスは前記周波数区分の端点の周波数に対応するインピーダンスである、請求項１から３の何れかに記載の等価回路を推定する方法。
【請求項５】
前記ステップ（ウ）における前記一つの等価回路の対応付けは、
前記選択された一つのインピーダンスの虚数項が正の場合には前記抵抗とコイルとの並列接続からなる等価回路要素を対応付け、
前記選択された一つのインピーダンスの虚数項が負の場合には前記抵抗とキャパシターとの並列接続からなる等価回路要素を対応付ける、
請求項１から４の何れかに記載の等価回路を推定する方法。
【請求項６】
前記ステップ（ウ）における前記一つの等価回路の対応付けは、前記選択された一つのインピーダンスの虚数項がゼロである場合には行わない、請求項１から５の何れかに記載の等価回路を推定する方法。
【請求項７】
前記等価回路要素リスト中に含まれる前記コイルからなる等価回路要素の個数が複数になった場合には、前記コイルからなる複数個の等価回路要素を一つにまとめる、請求項１から６の何れかに記載の等価回路推定方法。
【請求項８】
前記最適化は全て最小二乗法によって行う、請求項１から７の何れかに記載の等価回路推定方法。
【請求項９】
前記周波数範囲上の複数の周波数で測定された電気化学反応系のインピーダンスの測定値の系列を受け取り、請求項１から８の何れかに記載の方法に基づいて前記電気化学反応系の等価回路を推定する装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は実験により得られた電気化学インピーダンス・スペクトルから簡便に適切な等価回路を推定する方法に関する。
【背景技術】
【０００２】
例えば電池材料の電気化学特性を評価する一般的な手法として、電気化学インピーダンス測定が知られている。電気化学インピーダンス測定によって直接的に得られたデータから測定対象の等価回路の構成及び等価回路中の素子のパラメーターを解析する。これは、電池等の内部で生起している化学反応機構等の解析に有用な手法である。
【０００３】
しかし、精度が有限である有限個の測定データだけからでは完全に自動的に測定対象の等価回路を求めるのは一般には容易ではない。更には、そもそも測定データから等価回路を求める目的は、燃料電池の解析を例にとれば、電池内部で起こる様々な反応により複素インピーダンスプロットが複雑な形状になっているところ、それから個々の反応の寄与を、それぞれ等価回路中の回路要素として分離することで、個々の反応機構の解析を行うことである。このような目的を無視して良いのであれば、自動的な計算によって要素の組み合わせを複雑化して等価回路を組み上げていくことで、等価回路の複素インピーダンスプロットをいくらでも実測データの複素インピーダンスプロットに近づけることは、技術的には不可能ではない。しかし、これでは電気化学インピーダンス測定結果から等価回路を求める意味がなくなってしまう。従って、従来の電気化学インピーダンス測定装置で使用されていた電気化学インピーダンス解析プログラムでは、測定されたインピーダンスデータを利用者が観察して、測定対象で起こっている反応等を考えながら推定した等価回路及び当該等価回路中の各素子要素の初期パラメーターを利用者の経験に基づいて入力または大域最適化法によって概算した後に、最小自乗法を用いてこれらパラメーターの最適化を行っていた。この種の先行技術については、例えば非特許文献１、２を参照されたい。
【０００４】
利用者が等価回路及びその中に含まれる各要素の可変パラメーターの初期値である初期パラメーターを入力すると、これらを利用して、解析プログラムが最小自乗法を適用して可変パラメーターの最適値、つまりこの等価回路及び可変パラメーターの値から計算される各周波数での複素インピーダンスが測定された複素インピーダンスプロットに最小自乗法の意味で最適フィッティングするような可変パラメーターの値の組を計算する。利用者が推定した等価回路が最適なものであれば、上述の計算された可変パラメーターの組を有する等価回路から計算された複素インピーダンスプロットは、実測された複素インピーダンスプロットと非常によく一致する。しかし、推定した等価回路が不適切なものであった場合には、非線形最小自乗法は初期パラメーターが最適値から離れすぎていると発散することから上記計算が発散して結果が得られなかったり、あるいは収束してもその結果の複素インピーダンスプロットは実測されたものとかけ離れたものとなる。また、初期パラメーターが最適な値と大きく異なっていた場合でも、上述の理由で計算が発散したり、収束してもその収束値は一見もっともらしいが実際には不適切なものである等の問題が起こることがある。
【０００５】
従来は、測定されたインピーダンスの解析に当たって、抵抗、コンデンサー、ＣＰＥ、コイル、ワールブルグインピーダンス及びゲリッシャーインピーダンスを素子として扱い、これらを任意に直並列接続して等価回路を設定できるようになっていた。このように、従来技術においては設定できる等価回路の自由度が非常に大きいため、最適化を行った結果を適用した等価回路全体の複素インピーダンスを計算して利用者に提示することはできても、等価回路中の意味のある部分的な構造を切り出すための統一的な手法を設定することは困難であり、従って等価回路について最適化を行った結果が当該等価回路中の部分的な構造について妥当であるかどうかを検証することは困難であった。
【０００６】
等価回路及びその初期パラメーターを利用者が推定する際には、通常は利用者が経験や知識から妥当と考えられるものを与えたり、あるいは測定されたデータの一部から別途簡単な理論計算を行って、その結果から推定したり経験・知識からの推定結果の妥当性をある程度検証していた。しかし、利用者にこのような判断や作業を要求することは、利用者に過大な負担を強いるだけではなく、利用者の経験や知識により測定装置から最終的に得られる結果の妥当性に大きな差が出たり、あるいは同じ結果が最終的に得られるとしても、利用者の習熟度の低い場合には過大な試行錯誤を要することがあった。
【０００７】
また、Boukampは等価回路解析が可能なインピーダンス・スペクトルは測定された対数周波数範囲を７つ程度に分け、分けられた対数周波数の境界値に対応する特性周波数を持つ、抵抗（Ｒ）と容量（Ｃ）の並列要素（以後Ｒ／／Ｃ要素と記載）を足し合わせた回路式を用い、各定数が正負どちらでも取り得ると設定して最小自乗法でパラメーターを最適化するとスペクトルが再現できることを報告している（非特許文献３、４）。この手法は測定スペクトルが等価回路解析できるか否かを判定するクラマース−クローニヒ変換機能として複数のインピーダンス解析ソフトに実装されている。この種のソフトウエアの例としては、例えば非特許文献５を参照されたい。非特許文献５の４ページ目，"Kramers-Kronig test settings"章の最後に，緩和時間を１０セットのデータあたり７つが良い旨の記載がある。ここで緩和時間は初期設定で固定している。しかし、この方法では本来正の値を取るべきパラメーターが負の値を取り得るため、等価回路の推定に用いることができない。
【０００８】
その他、測定インピーダンス・スペクトル解析では経験則パラメーターを用いた要素を使用して解析することが多い。しかし、経験則パラメーターを含む要素から得られる最適化パラメーターには物理的解釈が困難な値が含まれているため、スペクトルの再現ができても解釈ができない欠点がある。また、経験則パラメーターを用いる解析法では緩和時間が対数周波数に対して分布を持つと仮定して解析を行う。しかし分布を持つため、変化を捉え難い欠点がある。
【発明の概要】
【発明が解決しようとする課題】
【０００９】
本発明の課題は、電気化学インピーダンス・スペクトルに基づいて、利用者の知識や判断をできるだけ必要とすることなく適切な等価回路を推定することにある。
【課題を解決するための手段】
【００１０】
本発明の一側面によれば、以下のステップ（ア）〜（コ）を設け、所定の周波数範囲にわたって測定された電気化学反応系のインピーダンスから、コイルからなる等価回路要素、キャパシターからなる等価回路要素、抵抗からなる等価回路要素、抵抗とコイルとの並列接続からなる等価回路要素、及び抵抗とキャパシターとの並列接続からなる等回路要素からなる群から選択された少なくとも一の等価回路要素の直列接続として構成される前記電気化学反応系の等価回路を推定する方法が与えられる。
（ア）前記周波数範囲上の複数の周波数で測定されたインピーダンスの測定値の系列を受け取る。
（イ）周波数の順に並べられた前記インピーダンスの測定値の系列を複数の周波数区分に分割するとともに、前記周波数区分の各々から一つのインピーダンスを選択する。
（ウ）前記周波数区分毎に、前記選択された一つのインピーダンス及びそのインピーダンスが測定された周波数に基づいて、前記群から選択された一つの等価回路要素を対応付け、前記対応付けられた等価回路要素からなる等価回路要素リストを得る。
（エ）ステップ（ウ）で得られた前記等価回路要素リスト中の前記等価回路要素を直列接続した回路のインピーダンスが前記測定されたインピーダンスの系列にフィッティングされるように、前記等価回路要素リスト中の等価回路要素に含まれる抵抗の抵抗値を正の値の範囲で最適化する。
（オ）ステップ（エ）によって最適化された等価回路要素リスト中で、隣接する前記周波数区分に対応する２つの等価回路要素の何れかが前記抵抗とコイルとの並列接続からなる等価回路要素であり、かつ周波数が高い方の前記周波数区分から選択された前記インピーダンスの実数項の値が周波数が低い方の前記周波数区分から選択された前記インピーダンスの実数項の値以下である場合には、前記周波数が高い方の周波数区分に対応する前記等価回路要素を前記コイルからなる等価回路要素に変更する。
（カ）前記等価回路要素リスト中に含まれる前記等価回路要素に対応する前記測定されたインピーダンスの実数項から前記等価回路要素のインピーダンスの計算値を減算し、前記減算の結果を前記等価回路要素リスト中の全ての前記等価回路要素について平均して平均実数項差を得る。
（キ）前記平均実数項差の値が正の場合には、前記平均実数項差の値の抵抗値を有する抵抗からなる等価回路要素を前記等価回路要素リストに追加する。
（ク）ステップ（キ）で得られた前記等価回路要素リスト中の前記等価回路要素を直列接続した回路のインピーダンスが前記測定されたインピーダンスの系列にフィッティングされるように、前記等価回路要素リスト中の等価回路要素に含まれる抵抗の抵抗値、コイルのインダクタンス及びキャパシタのキャパシタンスを、前記抵抗値、インダクタンス及びキャパシタンスの値の正負に拘束条件を付けずに最適化する。
（ケ）ステップ（ク）の最適化によって、抵抗の抵抗値、コイルのインダクタンス及びキャパシタンスのうちの少なくとも一つが負の値になった前記等価回路要素を前記等価回路要素リストから除去する。
（コ）ステップ（ケ）によって得られた前記等価回路要素リスト中の前記等価回路要素を直列接続した回路のインピーダンスが前記測定されたインピーダンスの系列にフィッティングされるように、前記等価回路要素リスト中の等価回路要素に含まれる抵抗の抵抗値、コイルのインダクタンス及びキャパシタのキャパシタンスを前記抵抗値、インダクタンス及びキャパシタンスが全て正であるとの拘束条件を付けて最適化する。
ここで、前記複数の周波数区分への分割は対数周波数軸上で等間隔となるように行われてよい。
また、前記インピーダンスの系列は、対数周波数軸上で等間隔で配置された複数の周波数について測定されたインピーダンスの系列であってよい。
また、ステップ（イ）における前記選択された一つのインピーダンスは前記周波数区分の端点の周波数に対応するインピーダンスであってよい。
また、前記ステップ（ウ）における前記一つの等価回路の対応付けは、
前記選択された一つのインピーダンスの虚数項が正の場合には前記抵抗とコイルとの並列接続からなる等価回路要素を対応付け、
前記選択された一つのインピーダンスの虚数項が負の場合には前記抵抗とキャパシターとの並列接続からなる等価回路要素を対応付けるものであってよい。
また、前記ステップ（ウ）における前記一つの等価回路の対応付けは、前記選択された一つのインピーダンスの虚数項がゼロである場合には行わなくてよい。
また、前記等価回路要素リスト中に含まれる前記コイルからなる等価回路要素の個数が複数になった場合には、前記コイルからなる複数個の等価回路要素を一つにまとめてよい。
また、前記最適化は全て最小二乗法によって行ってよい。
本発明の他の側面によれば、前記周波数範囲上の複数の周波数で測定された電気化学反応系のインピーダンスの測定値の系列を受け取り、前記何れかの方法に基づいて前記電気化学反応系の等価回路を推定する装置が与えられる。
【発明の効果】
【００１１】
本発明によれば、利用者に多くの知識や経験を求めることなく、電気化学インピーダンス・スペクトルの測定結果からその等価回路を推定することができる。更に、このようにして得られた等価回路はそのパラメーターは物理的な意味を保持するものであり、その等価回路に基づいて測定対象の系で起こっている現象を理解するのに有用である。また、解析されたＲ／／Ｃ、Ｒ／／Ｌ要素は固有の緩和時間を持つため、スペクトル解析から緩和時間変化を正確に評価できるようになり、これにより電池や腐食の条件変化や劣化挙動の解析に有効である。
【図面の簡単な説明】
【００１２】
【図１Ａ】要素分解が難しい測定インピーダンス・スペクトルの対数周波数と虚数項との関係を示す図。
【図１Ｂ】要素分解が難しい測定インピーダンス・スペクトルの対数周波数と実数項との関係を示す図。
【図１Ｃ】要素分解が難しい測定インピーダンス・スペクトルの実数項と虚数項との関係を示す図。
【図２Ａ】ステップ［１］〜［１０］で自動推定したインピーダンスモデル値と測定値との比較（対数周波数と虚数項の関係）を示す図。区分数Ｎは１０に設定した。細い二重線はモデルからのインピーダンス計算値を示し、細い二重破線は部分インピーダンス計算値を示す。
【図２Ｂ】ステップ［１］〜［１０］で自動推定したインピーダンスモデル値と測定値との比較（対数周波数と実数項の関係）を示す図。区分数Ｎは１０に設定した。細い二重線はモデルからのインピーダンス計算値を示し、細い二重破線は部分インピーダンス計算値を示す。
【図２Ｃ】ステップ［１］〜［１０］で自動推定したインピーダンスモデル値と測定値との比較（実数項と虚数項の関係）を示す図。区分数Ｎは１０に設定した。細い二重線はモデルからのインピーダンス計算値を示し、細い二重破線は部分インピーダンス計算値を示す。
【図３Ａ】ステップ［１１］で最適化された結果（対数周波数と虚数項との関係）を示す図。細い二重線はインピーダンス計算値を示す。
【図３Ｂ】ステップ［１１］で最適化された結果（対数周波数と実数項との関係）を示す図。細い二重線はインピーダンス計算値を示す。
【図３Ｃ】ステップ［１１］で最適化された結果（実数項と虚数項との関係）を示す図。細い二重線はインピーダンス計算値を示す。
【図４Ａ】ステップ［１２］でモデルを修正後、ステップ［１３］、［１４］で最適化した結果（対数周波数と虚数項との関係）を示す図。細い二重線は計算値を示す。
【図４Ｂ】ステップ［１２］でモデルを修正後、ステップ［１３］、［１４］で最適化した結果（対数周波数と実数項との関係）を示す図。細い二重線はインピーダンス計算値を示す。
【図４Ｃ】ステップ［１２］でモデルを修正後、ステップ［１３］、［１４］で最適化した結果（実数項と虚数項との関係）を示す図。細い二重線はインピーダンス計算値である。
【図５Ａ】ノイズによる大きな値飛びを含む測定インピーダンス・スペクトル（対数周波数と虚数項との関係）を示す図。灰色の長方形枠で囲んだ部分にノイズによる値飛びが存在する。
【図５Ｂ】ノイズによる大きな値飛びを含む測定インピーダンス・スペクトル（対数周波数と実数項との関係）を示す図。
【図５Ｃ】ノイズによる大きな値飛びを含む測定インピーダンス・スペクトル（実数項と虚数項との関係）を示す図。
【図６Ａ】ステップ［１］〜［１０］で自動推定したインピーダンスモデルからのインピーダンス計算値と測定値との比較（対数周波数と虚数項との関係）示す図。区分数Ｎは８に設定した。細い二重線はモデルからのインピーダンス計算値を示し、細い二重破線は部分インピーダンス計算値を示す。
【図６Ｂ】ステップ［１］〜［１０］で自動推定したインピーダンスモデルからのインピーダンス計算値と測定値との比較（対数周波数と実数項との関係）を示す図。区分数Ｎは８に設定した。細い二重線はモデルからのインピーダンス計算値を、細い二重破線は部分インピーダンス計算値を示す。
【図６Ｃ】ステップ［１］〜［１０］で自動推定したインピーダンスモデルからのインピーダンス計算値と測定値の比較（実数項と虚数項との関係）を示す図。区分数Ｎは８に設定した。細い二重線はモデル値を、細い二重破線は部分インピーダンス計算値を示す。
【図７Ａ】ステップ［１１］で最適化した結果（対数周波数と虚数項との関係）を示す図。
【図７Ｂ】ステップ［１１］で最適化した結果（対数周波数と実数項との関係）を示す図。
【図７Ｃ】ステップ［１１］で最適化した結果（実数項と虚数項との関係）を示す図。
【図８Ａ】ステップ［１２］でモデルを修正後、ステップ［１３］、［１４］で最適化した結果（対数周波数と虚数項との関係）を示す図。細い二重線は計算値を示す。
【図８Ｂ】ステップ［１２］でモデルを修正後、ステップ［１３］、［１４］で最適化した結果（対数周波数と実数項との関係）を示す図。細い二重線は計算値を示す。
【図８Ｃ】ステップ［１２］でモデルを修正後、ステップ［１３］、［１４］で最適化した結果（実数項と虚数項との関係）を示す図。細い二重線は計算値を示す。
【発明を実施するための形態】
【００１３】
本発明では、今回は電気化学インピーダンスの測定データ数値から等価回路を推定するため、非特許文献３、４等に記載されたような従来の方法を改良し、測定モデルに用いられる要素の特徴を考慮することによって効率的に等価回路を推定する。ここにおいて、測定対象の系の性質から、その等価回路はＲ、Ｃ、Ｌ、ＲＣ並列接続（Ｒ／／Ｃ）、ＲＬ並列接続（Ｒ／／Ｌ）の等価回路要素の直列接続となるとしてよい。ここで、先ず、行われたインピーダンス測定周波数範囲を複数領域に区分し、各領域の端点における測定結果のインピーダンスを主に使用することで、当該周波領域についての等価回路の要素を得る。ここで、リアクタンス分（インピーダンスの虚数項）が０なら等価回路要素とはしない。このようにして得られた等価回路の全要素の直列回路を測定周波数の全範囲内の測定結果全体にフィッティングする。この段階ではフィッティングのために値を調節するパラメーターは抵抗分のみとし、またその値を正に拘束する。次に等価回路要素Ｒ／／ＬのうちでＬとすべきものを判定して必要な修正を行ってから、今度は抵抗値、インダクタンス、キャパシタンスを調節対象のパラメーターとし、これらパラメーターの正負についての拘束条件なしでフィッティングを行う。この結果中から、何れかのパラメーターが負になった等価回路要素をすべて削除した上で、今度はこれらのパラメーターを正の値に拘束してフィッティングする。ここで、過剰な回路単位をまとめるなど処理を行うことで、最終的な等価回路を得る。
【００１４】
本発明における等価回路の推定は、具体的には以下のようにして行われる。
【００１５】
測定モデルは、電気化学インピーダンス・スペクトルを抵抗（Ｒ）、キャパシター（Ｃ）、コイル（Ｌ）、抵抗とキャパシターとが並列接続した要素（Ｒ／／Ｃ）、抵抗とコイルとが並列接続した要素（Ｒ／／Ｌ）を適数直列接続した回路で再現するものである。
【００１６】
測定モデルに用いる要素のインピーダンスは下記式で表される：
・抵抗のインピーダンス（ＺＲ）：ＺＲ＝Ｒ・・・式（１）
・キャパシターのインピーダンス（ＺＣ）：ＺＣ＝１／（ｊωＣ）・・・式（２）
・コイルのインピーダンス（ＺＬ）：ＺＬ＝ｊωＬ・・・式（３）
・抵抗とキャパシターの並列接続要素（ＺＲ／／Ｃ）：ＺＲ／／Ｃ＝Ｒ／（１＋ＲＣ×ｊω）・・・式（４）
・抵抗とコイルの並列接続要素（ＺＲ／／Ｌ）：ＺＲ／／Ｌ＝Ｒ×（Ｌ／Ｒ）×ｊω／（１＋（Ｌ／Ｒ）×ｊω）・・・式（５）
ここで、ｊは単位虚数（ｊ２＝−１）。ωは角周波数であり、周波数（ｆ）との間にω＝２πｆ（πは円周率）の関係がある。
【００１７】
これら要素は下記の通り分類される
・［周波数依存性を持たない要素］：抵抗（Ｒ）
・［周波数依存性を持つ要素］且つ［虚数部が正の値を持つ要素］：コイル（Ｌ）、抵抗とコイルとの並列要素（Ｒ／／Ｌ）
・［周波数依存性を持つ要素］且つ［虚数部が負の値を持つ要素］：キャパシター（Ｃ）、抵抗とキャパシターとの並列要素（Ｒ／／Ｃ）
・［１つの要素のみ回路に含まれる要素］：抵抗（Ｒ）、コイル（Ｌ）、キャパシター（Ｃ）
・［１つの要素のみ回路に含まれる要素］且つ［高周波側で現れる要素］：抵抗（Ｒ）、コイル（Ｌ）
・［１つの要素のみ回路に含まれる要素］且つ［低周波側で現れる要素］：キャパシター（Ｃ）
・［時定数を持つ要素］：抵抗とキャパシターとの並列接続要素（Ｒ／／Ｃ）、抵抗とコイルとの並列接続要素（Ｒ／／Ｌ）
【００１８】
時定数をτとしたとき、ＺＲ／／Ｃ、ＺＲ／／Ｌとの関係は下記の通りとなる：
ＺＲ／／Ｃ＝Ｒ／（１＋τ×ｊω），τ＝１／（ＲＣ）・・・式（６）
ＺＲ／／Ｌ＝Ｒ×τ×ｊω／（１＋τ×ｊω），τ＝Ｒ／Ｌ・・・式（７）
【００１９】
更に、周波数変化によるインピーダンス実数項の増減に着目することで、以下の分類を行う：
・［周波数の減少によりインピーダンス実数項が減少する要素］：Ｒ／／Ｌ
・［周波数の減少によりインピーダンス実数項が増加する要素］：Ｒ／／Ｃ
【００２０】
電気化学反応系の等価回路の要素となり得るこれら要素の特徴を考慮すると、以下のステップによって等価回路を推定することができる。
［１］測定されたインピーダンス・スペクトルデータを周波数、実数項、虚数項を１行とするデータ配列にする（配列順序は周波数の昇順または降順）。データは対数周波数で等間隔にＭ分割されて昇順または降順に並べられたものを用いる。一般的にＭは周波数１桁あたり１０〜２０である。
［２］分割する周波数区分数を設定する。分割する周波数区分の数（Ｎ）の最大値はデータ行数／２であるが、３〜１２程度が望ましい。
［３］データ配列の開始行を１行目、最終行をＭ行目とすると、１行目とＭ行目を含めてＮ個のデータを、データ行間隔が等間隔になるように、このデータ配列中から抽出する。抽出されたデータをＮ行配列とする。ＮとＭの関係により単純に割り算して求まるデータ行間隔（Ｍ／Ｎ）が小数になる場合は、四捨五入または小数切り捨て操作により整数化して抽出する。抽出されたデータ行をＮ行配列中のＫ行目データ配列とする。
［４］等価回路推定の出発点として、Ｋ行目データ配列の各データ点（データ行）毎に一つの等価回路要素（Ｒ／／Ｌ要素、Ｒ／／Ｃ要素等）を割り当て、これら等価回路要素を直列接続した等価回路を得る。具体的には、Ｋ行目データ配列の各データ点毎に以下の操作を行う：
■Ｋ行目データ配列のデータ点中の虚数項が正の場合は当該データ点に対応する等価回路要素をＲ／／Ｌ要素とし、当該データ点の周波数及び式（７）からＫ番目のＲ値（ＲＫ）を計算し、さらにτ値（τＫ）及びＲＫからＬＫを計算する。なお、ここで添え字ＫはＮ行配列中のＫ行目データ配列から求めた値であるという意味で付与しており、実際にはＫ行目データ配列のデータ点毎に一つずつ求められる。以下同様である。
■Ｋ行目データ配列のデータ点中の虚数項が負の場合はＲ／／Ｃ要素とし、当該データ点の周波数及び式（７）からＲ値（ＲＫ）とＣ値（ＣＫ）を計算する。
■Ｋ行目データ配列のデータ点の虚数項が０の場合は操作なしとする。つまり、このデータ点に対応する等価回路要素は存在しないとする。
［５］ステップ［４］で求めた虚数項の関係及び実測された虚数項を用い、ＲＫ１〜ＲＫＮをパラメーターとして最小二乗を行い、ＲＫ１〜ＲＫＮの最適値を求める。最小二乗の際、ＲＫ１〜ＲＫＮは正の値に拘束する。
［６］１番目の要素タイプがＲ／／Ｌの場合、Ｋ行目データ配列中の１番目と２番目の行に対応するインピーダンスを測定した周波数での測定データ実数項を確認し、１番目の実数項値が２番目の実数項値以下である場合は１番目の要素タイプとＲ１値をそれぞれＬと０に変更する。１番目の実数項が２番目の実数項よりも大きい場合は変更無しとする。もちろん、この場合には抵抗Ｒ１の値である０それ自体には意味がない。単に、一見すると意味のありそうなランダムな値がここに入ると利用者が誤解するおそれがあるので、それを避けるために０を入れるにすぎない。
［７］上記操作を最大Ｎ番目の要素タイプまで繰り返す。
［８］ステップ［７］で決められた最大Ｎ個の要素について、式（６）、（７）から夫々のＣＫ、ＬＫを計算する。
［９］計算されたＲＫ、ＣＫ、ＬＫ値と要素タイプから複素インピーダンスを計算し、実数項の計算値を求める。
［１０］測定されたインピーダンス実数項列とステップ［９］で計算された実数項列の差を計算する。差の平均値が正の場合は回路要素Ｒを追加する。追加する回路要素Ｒの値は当該差の平均値に設定する。
［１１］推定された等価回路と実測データに対し、パラメーターの正負を拘束しない条件で複素非線形最小二乗法（Complex Nonlinear Least Squares）を適用する。複素非線形最小二乗は一つの入力変数（本願の場合は周波数）に対して２つの出力関数（実数項と虚数項）がある複素関数の最小二乗である。今現在，複素関数に対してパラメーターを実数のみに固定した条件での汎用最小二乗法は存在していない。複素関数に対しては実数関数に使用する非線形最小二乗法を特殊な方法で使用するため、複素非線形最小二乗と呼ばれている。さらに詳細には、例えば非特許文献６を参照されたい。なお、最小二乗近似によるフィッティングはこのステップだけではなく、ステップ［５］でも行っている。このように最小自乗近似を２回に分けて行う理由は、ステップ［５］で最初に最小自乗近似を行う段階では、パラメーターの初期値と実際のパラメーターの値とのかい離が大きいのが通常であるため、１ステップで複素非線形最小二乗を試みると特異行列エラーになる可能性が高いからである。そこで、先ずステップ［５］でエラーにならない程度のパラメーターの値に収束させ、次にステップ［１１］において弱い拘束条件下で複素非線形最小二乗を行う。ステップ［１１］においてエラーが生じた場合はステップ［３］のＮ値を１減じてからステップ［１１］までの操作を繰り返す。［１２］最適化によりパラメーターＬＫ、ＣＫ、ＲＫの何れか一つでも負に収束した要素は全て削除によってモデルを修正する。なお、このモデル修正を行った後では、前のステップに戻ることなく処理を進める。その他、測定データのノイズによる要素誤認と判断できるものについては削除する。誤認された要素の削除に当たっては、最適化された結果を利用者に提示し、利用者が要素誤認か否かを判断してもよい。あるいは、誤認された（あるいは誤認の可能性の高い）要素の特徴を予めコンピュータに蓄積したり、利用者による要素誤認か否かの判断をコンピュータに学習させることにより、この判断の一部あるいは全部を自動化することもできる。
［１３］最低周波数でのＲ／／Ｃ要素（最適化パラメーター値表の最下行Ｒ／／Ｃ要素）の抵抗値が大きく、その標準誤差（標準偏差を試料数の平方根で除算したもの）も大きい場合はＲ／／Ｃ要素をキャパシター要素に変更しても良い。なお、最適化パラメーター値表の構成を説明すれば、最適化パラメーター表の各行は一つの等価回路要素に対応する。また、一つの等価回路要素はＫ番目データ配列の一つのデータ点（行）に対応する（もちろん、処理過程で削除される等価回路要素もある）。従って、最適化パラメーター表中の行はＫ番目データ配列中の対応する行の順に配列されている。その結果、Ｋ番目データ配列が周波数の降順に配列されている場合には、最適化パラメーター表中の下にある行ほど、それに対応するＫ番目データ配列中の行のデータの測定周波数が低いことになる。「最低数波数でのＲ／／Ｃ要素」と言う時の「最低周波数での」とはこの意味である。
［１４］パラメーター値が全て正になるように拘束し、複素非線形最小二乗法でパラメーターの最適化を行う。複素非線形最小二乗を操作時にエラーが出る場合は、ステップ［１２］を行い、モデルの再修正を行う。
［１５］ステップ［１３］で最適化されたパラメーターの標準誤差値が最適化パラメーター値よりも小さいことを確認する。パラメーターの標準誤差値が最適化パラメーター値よりも大きい値になった場合は、ステップ［２］の区分数を減らすか、誤差値が大きい複数要素を１つにするよう、要素の一部を削除してモデルの再構築を行う。このステップの意義を説明すれば、ステップ［３］のＮ行配列から求められた最適化パラメーターは対数周波数に対して降順または昇順に並ぶ。近接する同じ要素タイプの誤差値がパラメーター値よりも大きい場合は分割数が多すぎることを意味する。従って、近接する同じ要素タイプを一つに統合するようにモデルを再構築するのである。
［１６］最適化モデルは、最小二乗残差和が小さく、全ての最適化パラメーター値がパラメーターの誤差値よりも大きいものを選択する。
【００２１】
上記一連のステップはプログラムコード化して半自動または全自動で行うことができる。
【００２２】
電気化学系の等価回路を推定するための上で説明した一連のステップは、以下の基本的な考え方に基づいている：
（１）ある周波数におけるインピーダンス虚数項が負の値である場合はＲ／／ＣまたはＣが優勢なインピーダンスである。
（２）ある周波数におけるインピーダンス虚数項が正の値である場合はＲ／／ＬまたはＬが優勢なインピーダンスである。ここで更に、
（２−１）Ｎ配列が対数周波数に対して降順に並べられ，かつＫ行目配列のインピーダンス虚数項が正の値であり，さらにＫ行目データ配列のインピーダンス実数項がＫ＋１行目データ配列のインピーダンス実数項よりも大きい場合にはＲ／／Ｌが優勢なインピーダンスである（Ｒ／／Ｌ要素インピーダンスは周波数の減少によりインピーダンス実数項が減少する性質のため）。
（２−２）Ｎ配列が対数周波数に対して降順に並べられ，かつＫ行目配列のインピーダンス虚数項が正の値であり，さらにＫ行目データ配列のインピーダンス実数項がＫ＋１行目データ配列のインピーダンス実数項と同一あるいは小さい場合はＬが優勢なインピーダンスである（Ｌ要素のインピーダンスのインピーダンス実数項は周波数によらず一定値を取る。ただし、Ｒ／／Ｃ要素と競合している場合は周波数の減少により実数項が増加する）。
（２−３）複数の選択周波数でＬが優勢になった場合、これらの要素は一つのＬにまとめる。
（３）インピーダンスの虚数項と実数項の間にはクラマース−クローニヒの関係が保持されているため、虚数項を最適化すると、実数項は曲線形状の等しい値が得られる。したがって測定されたインピーダンス実数項と、虚数項計算で得られたパラメーターから求めた計算インピーダンス実数項の差の平均値が抵抗（Ｒ）に対応する。この段階ではＲ値が正負に関わらずこの平均値を値として有する等価回路要素Ｒを追加する（ステップ［１０］に対応）。
（４）初期の周波数区分数が最適解を与える区分数よりも多い場合、パラメーターを非拘束した非線形最小二乗を施すと、過剰な要素のパラメーターが負の値に収束する。
（５）負のパラメーターに収束した要素は過剰分であるため、これらの要素を削除によってモデルを修正する。このモデル修正を行った後、前のステップに戻ることなく処理を進める。
（６）パラメーターが全て正の値に収束するモデルが得られた場合、最適化パラメーター値とその標準誤差値を比較し、Ｋ行目とＫ＋１行目の標準誤差がパラメーター値よりも大きい場合は過剰区分として、ＲＫとＲＫ＋１で値の小さい方を削除し一つにまとめる。まとめる際にはＲＫ値をＲＫ＋ＲＫ＋１に変更してからＫ＋１行目パラメーター行を削除する。
（７）最低周波数でのＲ／／Ｃ要素の抵抗値が他要素の抵抗値と比べて１桁以上大きく、その標準誤差値も大きい場合は、Ｃ要素をＲ／／Ｃで割り付けている可能性が高いため、Ｒ／／ＣをＣに置き換える。
【００２３】
区分する周波数区分数は最適解を与える区分数よりも多い値から試すのが良い。経験的には５〜１２程度が良い。ステップ［２］の周波数区分数決定については、複数通りの区分数を用いてステップ［３］〜［５］を行い、測定インピーダンスデータと計算値の二乗残差和が最も小さくなる区分数を自動で判定しても良い。また、ステップ［１２］はプログラムによる自動操作にしてもよい。
【００２４】
電気化学インピーダンスで現れうるワールブルグインピーダンスとゲリッシャーインピーダンス（これらのインピーダンスについては非特許文献７を参照）は上記回路でも再現できることから、本発明による解析を行った後、部分インピーダンス（最適化パラメーター値表の各行で表されるインピーダンス）の部分和を検証し、ワールブルグインピーダンスおよびゲリッシャーインピーダンスの存在有無を予測できる。
【００２５】
なお、上で説明したパラメーターの最適化では、具体的にはモデル値と測定値の自乗残差和を目的関数にして、近似ヘッセ行列の逆行列を計算する過程を含む方法を用いることができる。近似へッセ行列の逆行列を計算する手法としてはニュートン・ガウス法、レベンパーグ・マーカート法、直行距離回帰法などがある。ここで、最初期の区分数が最適化可能な区分数よりも多い場合、近似ヘッセ行列は逆行列を持たなくなるため、逆行列計算を試みた時点で計算エラーが発生する。この計算エラーが発生した場合は、区分数を一つ減じ、初期最適化を行う（パラメーターの正負を拘束しない最小自乗を行う）。
【００２６】
また、本発明による最適化はクラマース・クローニヒ変換が成り立つデータに対して有効である。従って、インピーダンスモデル推定前にクラマース・クローニヒ変換計算を行い、実測データが変換関係を満たしていることを確認するのが望ましい。ここにおいて、実測データの一部のみがクラマース・クローニヒ変換関係を満たしている場合、その関係が満たされている周波数領域のデータのみを用いてインピーダンスモデルを推定することが望ましい。従って、モデル推定を行う実測値データを、何らかの手段によって制限する機能を持つことが有効である。例えば、グラフ上のカーソルを用いて周波数領域を表示装置の画面上で選択できるようにすればよい。
【００２７】
以上のようにして推定された等価回路は、回路を表す図形の形態で表示装置やプリンタ等の出力装置から出力して良い。あるいは、図形以外の形態で出力してもよいし、以降のコンピュータプログラムでの処理に好都合な何らかのデータの形態で次段のプログラムに引き渡す、コンピュータ可読媒体への出力を行う、通信回線を介して送出する等の形態をとることもできる。
【実施例】
【００２８】
＜基本的な実施例＞
以下ではある電気化学系の電気化学インピーダンス・スペクトルの実測データに基づいて、その等価回路を推定する過程を説明する。実測データとしては要素分解が困難であるものを選んだ。この実測データの対数周波数（ｌｏｇｆ）とインピーダンス虚数項（Ｚｉｍｇ）との関係、対数周波数（ｌｏｇｆ）とインピーダンス実数項（Ｚｒｅａｌ）との関係及び実数項と虚数項との関係を夫々図１Ａ、図１Ｂ及び図１Ｃに示す。
【００２９】
なお、図１Ａ、図１Ｂ及び以降の図（枝番Ａ及びＢが付された図）は形式が統一されていて、対数周波数（ｌｏｇｆ）とインピーダンス虚数項（Ｚｉｍｇ）との関係を示す枝番Ａが付された図及び対数周波数（ｌｏｇｆ）とインピーダンス実数項（Ｚｒｅａｌ）との関係を示す枝番Ｂが付された図には、その下部にインピーダンス虚数項または実数項について、実測値と計算値との差であるΔＺｉｍｇ、ΔＺｒｅａｌも示す。ただし、図１Ａ及び図１Ｂは単に実測値を示し、この段階では計算値は存在していないので、ΔＺｉｍｇ及びΔＺｒｅａｌは何れもゼロ（０）であるとして図中に表記する。このゼロ表記は、計算値が存在しない図５Ａ及び図５Ｂでも行っている。
【００３０】
図１Ａ〜図１Ｃに示されたデータに基づいて、上記ステップ［１］〜［１０］に従って自動推定することで測定対象の電気化学系のインピーダンスモデル、つまり等価回路を得た。この実施例ではＮ＝１０と設定した。得られたインピーダンスモデルについてそのインピーダンスを計算し、これにより得られた計算値を図１Ａ〜図１Ｃに示された実測値のグラフにオーバーレイしたものを図２Ａ〜図２Ｃに示す。図２Ａ及び図２Ｂには実測と計算値との差も示した。
【００３１】
ステップ［１］〜［１０］に従って推定された等価回路と実測データに対し、ステップ［１１］においてパラメーターの正負を拘束しない条件で複素非線形最小二乗法を適用することで、最適化を行った。その結果得られた等価回路のインピーダンスの計算値を図２Ａ〜図２Ｃと同様にして、図３Ａ〜図３Ｃに示す。
【００３２】
ステップ［１１］で最適化されたパラメーター一覧を以下の表に示す。表中、w_element_label1項は要素タイプを、w_R1及びw_R1_std_uc1はＲＫ値及びその標準誤差値を、並びにw_T1及び w_T1_std_uc1はＣＫまたはＬＫ及びそれらの標準誤差値を夫々示す。また、灰色の枠で囲まれた要素（Point 0、Point 7、Point 9）は負の値に収束したパラメーターを含む要素を示す。
【００３３】
【表１】
2018055402.tif 000002
【００３４】
次に、ステップ［１２］において、最適化されたパラメーターが負の値に収束した要素を削除して、モデルの再構築を行った。再構築されたモデルに対してステップ［１３］及び［１４］により更に最適化を行った。その結果得られたインピーダンスモデル、即ち等価回路のインピーダンスの計算値を図２Ａ〜図２Ｃ等と同様にして、図４Ａ〜図４Ｃに示す。これらの図から見て取ることができるように、ステップ［１３］及び［１４］による最適化後においては、電気化学系のインピーダンスの実測値（図中「＋」マークで示す）と等価回路全体のインピーダンスの計算値（図中、連続した二重線で示す）とは非常によく一致することが確認された。
【００３５】
ステップ［１３］、［１４］で最適化されたパラメーター一覧を下表に示す。表中、w_element_label1は要素タイプを、w_R1及びw_R1_std_uc1はＲＫ値及びその標準誤差値を、並びにw_T1及び w_T1_std_uc1はＣＫまたはＬＫ及びそれらの標準誤差値を夫々示す。
【００３６】
【表２】
2018055402.tif 000003
【００３７】
＜実測されたインピーダンスの値に値飛び（ノイズ等による不連続点）が存在する場合の実施例＞
このような実測値の例を図５Ａ〜図５Ｃに示す。図５Ａからわかるように、低周波数側（ｌｏｇ ｆ＝１付近）で実測インピーダンスの実部の値が大きく飛んでいて、ノイズ等の影響を示唆している。この実測値に対して最初の実施例と同様な処理を行うことで等価回路を求めた。
【００３８】
すなわち、先ず上記ステップ［１］〜［１０］に従って自動推定することで図５Ａ〜図５Ｃに示す実測値に基づく等価回路を得た。この実施例ではＮ＝８と設定した。得られた等価回路についてそのインピーダンスを計算し、これにより得られた計算値を、図２Ａ〜図２Ｃと同様の態様で図６Ａ〜図６Ｃに示す。
【００３９】
上で得られた最初の等価回路について、ステップ［１１］においてパラメーターの正負を拘束しない条件で複素非線形最小二乗法を適用することで、最適化を行った。その結果得られた等価回路のインピーダンスの計算値を図３Ａ〜図３Ｃと同様の態様で図７Ａ〜図７Ｃに示す。
【００４０】
ステップ［１１］で最適化されたパラメーター一覧を下表に示す。表中、w_element_label1項は要素タイプを、w_R1及びw_R1_std_uc1はＲＫ値及びその標準誤差値を、並びにw_T1及び w_T1_std_uc1はＣＫまたはＬＫ及びそれらの標準誤差値をそれぞれ示す。また、灰色の枠で囲まれた要素（Point 6〜Point 8）は測定データに含まれるノイズにより誤認された要素を示す。
【００４１】
【表３】
2018055402.tif 000004
【００４２】
次に、ステップ［１２］において、最適化されたパラメーターで負の値に収束した要素を削除して、モデルの再構築を行った。再構築されたモデルに対してステップ［１３］及び［１４］により更に最適化を行った。その結果得られた等価回路のインピーダンスの計算値を図６Ａ〜図６Ｃ等と同様な態様で、図８Ａ〜図８Ｃに示す。これらの図から見て取ることができるように、実測値に値飛びがある図５Ａ〜図５Ｃの場合でも、ステップ［１３］及び［１４］による最適化後においては、電気化学系のインピーダンスの実測値（図中「＋」マークで示す）と等価回路全体のインピーダンスの計算値（図中、連続した二重線で示す）とは非常によく一致することが確認された。
【００４３】
ステップ［１３］、［１４］で最適化されたパラメーター一覧を下表に示す。表中、w_element_label1は要素タイプを、w_R1及び w_R1_std_uc1はＲＫ値及びその標準誤差値を、並びにw_T1及び w_T1_std_uc1はＣＫまたはＬＫ及びそれらの標準誤差値を夫々示す。
【００４４】
【表４】
2018055402.tif 000005
【００４５】
この実施例では実測値の値飛び等の明らかに不自然なデータ点をそのままにして等価回路推定を行った。しかし、等価回路推定を開始する前に、例えばこのような明らかに不自然なデータについては手作業で自然な値に修正したり、あるいはコンピュータによる学習によって不自然さの判別及び自動修正を行うことも可能である。また、値飛び等の不自然なデータのある周波数領域（図１Ａではグラフの左端付近の灰色のやや太い線で囲んだ領域）を等価回路推定の際の処理対象となる周波数範囲から除外してもよい。また、実測値に大きなノイズが乗るなど、インピーダンス測定の実測値が現実にあり得るデータから大きくかい離した場合には上述した等価回路推定の過程でエラーが起こり、推定が不可能となることも起こり得る。このようなエラーの発生を抑止するため、最初に実測されたインピーダンスデータに対してクラマース・クローニヒ変換を施してみて、この変換が成り立つ周波数領域だけについて上記の等価回路推定を行うようにしてもよい。
【００４６】
なお、最初の実施例の実測値においても、図１Ａにはｌｏｇ ｆ＝１．７付近で一点だけデータ飛びが見られる。それに対して、本実施例における図５Ａのノイズは図１Ａの単発ノイズと異なり，ある周波数帯に広がって現れている。図１Ａのようなノイズの場合には、元データ配列中に存在していてもそれがＫ行目データに当たる確率が低いため、等価回路推定への影響は現れ難い。これに対して、図５Ａの場合は多数のデータ点にノイズが乗っているため、等価回路推定に影響が現れやすくなる。そのため、図５Ａのような場合には、上でも言及したように、通常の推定のプロセスを開始する前にそのようなノイズが存在する領域を削除するなどの処置を取るのが好ましい。スペクトル形状からノイズ等の影響を完全に判定するための簡便かつ一般的なアルゴリズムがないため、起こり得るあらゆるノイズに対応するためにはこの段階に人間が介入して判断を行うことが必要となる。図５Ａは人による判定が好ましい場合の例である。
【００４７】
しかし、現実の電気化学反応系では通常は起こる可能性の少ないインピーダンススペクトルのパターンの一部を認識して、利用者に対してノイズの影響の可能性を指摘したり、あるいは、上述したように、最初に実測されたインピーダンスデータに対してクラマース・クローニヒ変換を施してみることで、この変換が成り立たない場合に利用者に対して警告を与えるなど、利用者による上記判断を支援することは実現可能である。また、ノイズの影響等の何らかの例外事象がインピーダンス測定対象の周波数範囲中のある区間で起こっている可能性を判定できる場合には、そのような区間を利用者に対して提示することで、利用者がその区間を外して推定処理を行うようにすることもできる。また、そのような区間を自動的に外して推定処理を進める自動判定処理を行ったり、あるいはこの自動判定を行うモードを事前に選択できるようにしてもよい。更に、通常の推定を行っても残存する誤差（ΔＺｉｍａｇ、ΔＺｒｅａｌ等）が妥当とみなすことができる値や事前に設定しておいた値よりも大きくなった場合に、ノイズの影響等の例外事象が起こった旨の警告を利用者に与えることもできる。
【先行技術文献】
【非特許文献】
【００４８】
【非特許文献１】J. R. Macdonald, "LEVM/ LEVMW MANUAL, Complex nonlinear least squares imittance, inversion, and simulation fitting programs for Windowsand MS-DOS"Ver. 8.12 (June 2013) 、[平成２８年９月２７日検索］インターネット＜URL：http://jrossmacdonald.com/jrm/wp-content/uploads/LEVMMANUAL-5-12-14.pdf＞
【非特許文献２】B. Boukamp, "A nonlinear least squares fit procedures for analysis of immittance data of electrochemical systems", Solid State Ionics, 20 (1986) 31-44.
【非特許文献３】B. A. Boukamp, "A Linear Kronig-Kramers Transform Test for Immittance Data Validation" J. Electrochem. Soc., 142 (1995) 1885.
【非特許文献４】B. A. Boukamp, "Electrochemical Impedance Spectroscopy"  Nano-Electrocatalysis, U. Lieden, 24-28, November 2008 、[平成２８年９月２７日検索］インターネット＜URL：http://www.lorentzcenter.nl/lc/web/2008/317/presentations/Boukamp.pdf＞
【非特許文献５】B.A. Boukamp, "Kramers-Kronig test for Windows", KK-Windowsマニュアル 、[平成２８年９月２７日検索］インターネット＜URL：https://www.utwente.nl/tnw/ims/publication/downloads/KK-windows.zip＞
【非特許文献６】K. Kobayashi, K. Terabe, T. Sukigara, Y. Sakka, Solid State Ionics, 249-250 (2013) 78-85.
【非特許文献７】[平成２８年９月２７日検索］インターネット＜URL：http://www.toyo.co.jp/material/products/detail/ZPLOT.html＞
【図１Ａ】
2018055402.tif 000006
【図１Ｂ】
2018055402.tif 000007
【図１Ｃ】
2018055402.tif 000008
【図２Ａ】
2018055402.tif 000009
【図２Ｂ】
2018055402.tif 000010
【図２Ｃ】
2018055402.tif 000011
【図３Ａ】
2018055402.tif 000012
【図３Ｂ】
2018055402.tif 000013
【図３Ｃ】
2018055402.tif 000014
【図４Ａ】
2018055402.tif 000015
【図４Ｂ】
2018055402.tif 000016
【図４Ｃ】
2018055402.tif 000017
【図５Ａ】
2018055402.tif 000018
【図５Ｂ】
2018055402.tif 000019
【図５Ｃ】
2018055402.tif 000020
【図６Ａ】
2018055402.tif 000021
【図６Ｂ】
2018055402.tif 000022
【図６Ｃ】
2018055402.tif 000023
【図７Ａ】
2018055402.tif 000024
【図７Ｂ】
2018055402.tif 000025
【図７Ｃ】
2018055402.tif 000026
【図８Ａ】
2018055402.tif 000027
【図８Ｂ】
2018055402.tif 000028
【図８Ｃ】
2018055402.tif 000029