Main Content

劣化するバッテリー システムの非線形の状態推定

この例では次を使用します。

スクリプトを開く

この例では、Simulink® のアンセンテッド カルマン フィルターを使用して非線形システムの状態を推定する方法を示します。また、イベントベースのカルマン フィルターを開発してシステム パラメーターを更新し、状態推定の精度を高める方法についても説明します。この例は、Control System Toolbox™ または System Identification Toolbox™ のいずれかで実行されます。この例では Predictive Maintenance Toolbox™ は必要ありません。

概要

次のような等価回路をもつバッテリー モデルについて考えます [1]。

モデルは電圧源 $E_m$、直列抵抗器 $R_0$、および 1 つの RC ブロック $R_1$ と $C_1$ で構成されています。バッテリーは充電と放電のサイクルを交互に繰り返します。この例では、測定された電流、電圧およびバッテリー温度の値を使用してバッテリー モデルの充電状態 (SOC) を推定します。バッテリーは非線形システムであると仮定し、アンセンテッド カルマン フィルターを使って SOC を推定します。バッテリーの容量は各放電サイクルごとに低下するため、SOC の推定値が不正確になります。イベントベースの線形カルマン フィルターを使用して、バッテリーが充電状態から放電状態に切り替わるときにバッテリー容量を推定します。推定した容量は、バッテリーの健康状態を示すために使用できます。

Simulink モデルは、バッテリー モデル、Unscented Kalman Filter ブロック、イベントベースの Kalman Filter ブロックの 3 つの主要コンポーネントを含みます。さらに詳しい説明については、以下の節を参照してください。

open_system('BatteryExampleUKF/')

バッテリー モデル

Simscape™ 言語を使用して熱効果のあるバッテリー モデルを実装します。

バッテリー モデルの状態遷移の方程式は次のようになります。

$$ \frac{d}{dt} \left(
\begin{array}{cc}
 SOC \\
 U_{1}
\end{array} \right) = \left(
\begin{array}{cc}
 0 \\
 -\frac{1}{R_1(SOC,T_b)*C_1(SOC,T_b)}U_1
\end{array} \right) + \left(
\begin{array}{cc}
 -\frac{1}{3600*C_q} \\
 \frac{1}{C_1(SOC,T_b)}\end{array} \right)I
+ W
$$

ここで、$R_1(SOC,T_b)$ と $C_1(SOC,T_b)$ は RC ブロック内の熱および SOC 依存の抵抗器とコンデンサ、$$ U_1 $$ はコンデンサ $C_1$ にかかる電圧、$I$ は入力電流、$T_b$ はバッテリーの温度、$C_q$ はバッテリーの容量 (単位: Ah)、$W$ はプロセス ノイズです。

次の図に示すように、入力電流は、バッテリーの放電時はランダムに生成されたパルス、バッテリーの充電時は一定です。

測定方程式は以下で与えられます。

$$ \begin{array} {ll}
E = E_m(SOC,T_b) - U_1 - IR_0(SOC,T_b) + V
\end{array} $$

ここで、$E$ は測定された電圧出力、$R_0(SOC,T_b)$ は直列抵抗器、$E_m = E_m(SOC,T_b)$ は電圧源からの起電力、$V$ は測定ノイズです。

モデル内では $R_0, R_1, C_1$ および $E_m$ が 2D ルックアップ テーブルで、SOC とバッテリーの温度に依存します。ルックアップ テーブルのパラメーターは、実験データを使用して特定されます [1]。

充電状態 (SOC) の推定

Unscented Kalman Filter ブロックを使用するには、状態方程式と測定方程式の MATLAB® 関数または Simulink 関数が定義されていなければなりません。この例では Simulink 関数を使用する方法を示します。アンセンテッド カルマン フィルターは離散時間フィルターであるため、最初に状態方程式を離散化します。この例ではオイラー法による離散化を使用します。サンプリング時間を $$ Ts $$ とします。一般的な非線形システム $\dot{x} = f(x,u)$ に対して、システムを $x_{T+1} = x_{T} + f(x_{T},u_{T})*Ts$ として離散化できます。

非線形バッテリー システムの状態ベクトルは次のとおりです。

$$ x_T = \left( \begin{array}{cc} SOC_T \\ U_{1_T} \end{array} \right) $$.

オイラー法による離散化を適用すると、次の方程式が得られます。

$$
\left(
\begin{array}{cc}
 SOC_{T+1} \\
 U_{1_{T+1}}
\end{array} \right) = \left(
\begin{array}{cc}
 SOC_{T} \\
 U_{1_{T}}
\end{array} \right) + \left(
\begin{array}{cc}
 -\frac{1}{3600*C_q}I \\
 -\frac{1}{R_1(SOC_T,T_b)*C_1(SOC_T,T_b)}U_1+\frac{1}{C_1(SOC_T,T_b)} I
\end{array} \right)Ts
+ W_T
$$

離散化された状態遷移方程式は、以下に示す "batteryStateFcn" という Simulink 関数として実装されます。関数入力 x は状態ベクトルで、関数出力 xNext は、離散化した状態遷移方程式を使って求められた次のステップでの状態ベクトルです。x と xNext の信号の次元とデータ型を指定する必要があります。この例では、x と xNext の信号の次元は 2、データ型は double です。追加の入力は、温度、推定容量、および電流です。追加の入力は状態遷移方程式への入力であり、UKF ブロックでは必須でない点に注意してください。

測定関数も、以下に示す "batteryMeasurementFcn" という名前の Simulink 関数として実装されます。

ブロック パラメーターを以下のように設定します。

[システム モデル] タブで、ブロック パラメーターを次のように指定します。

次のパラメーターを指定します。

  • 状態遷移の関数: batteryStateFcn。

離散化された状態遷移方程式を実装する、以前に定義した Simulink 関数の名前。

  • プロセス ノイズ: 時変共分散 $\left[ \begin{array}{cc} 2e-8 & 0 \\ 0 & 3e-7 \end{array} \right]$ の Additive。Additive は、ノイズ項が最終信号に直接加算されることを意味します。

SOC および $U_1$ のプロセス ノイズは、バッテリー システムの動的な特性に基づいて推定されます。バッテリーは定格容量が 30 Ah で、放電/充電サイクルにおける電流の振幅が平均して 15 A です。したがって、1 度の放電処理または充電処理に約 2 時間 (7200 秒) かかります。最大変化量は SOC で 100%、$U_1$ で約 4 ボルトです。

SOC と $U_1$ のステップあたりの最大変化量は、それぞれ $max(|dSOC|) \approx \frac{100\%}{3600*2}*Ts$ と $max(|dU_1|) \approx \frac{4}{3600*2}*Ts$ です。ここで、$T_s$ はフィルターのサンプリング時間です。この例では $T_s$ が 1 秒に設定されます。

プロセス ノイズ $W$ は次のようになります。

$$ W = \left[ \begin{array}{cc} (max(|dSOC|))^2 & 0 \\ 0 & (max(|dU_1|))^2 \end{array} \right] \approx \left[ \begin{array}{cc} 2e-8 & 0 \\ 0 & 3e-7 \end{array} \right] $$.

  • 初期状態: $$ \left( \begin{array}{cc} 1 \\ 0 \end{array} \right) $$。

SOC の初期値は 100% (満充電のバッテリー) と仮定され、$U_1$ の初期値は 0 に設定されます。これは、$U_1$ の事前情報がないためです。

  • 初期の共分散:

初期の共分散は初期推定値の正確さと信頼度を示します。初期推定値の最大誤差は、SOC で 10%、$U_1$ で 1 V であると仮定します。初期共分散行列は $$ \left[ \begin{array}{cc} 0.01 & 0 \\ 0 & 1 \end{array} \right] $$ に設定されます。

  • Unscented 変換のパラメーター:設定は [2] に基づきます。

     - Alpha: 1. Determine the spread of sigma points around x. Set Alpha to be 1 for larger spread.
     - Beta: 2. Used to incorporate prior knowledge of the distribution. The nominal value for Beta is 2.
     - Kappa: 0. Secondary scaling parameter. The nominal value for Kappa is 0.

  • 測定の関数: batteryMeasurementFcn。

測定関数を実装する、以前に定義した Simulink 関数の名前。

  • 測定ノイズ: 時不変共分散 1e-3 の Additive。

測定ノイズは、測定機器の精度に基づいて推定されます。バッテリー電圧を測定する電圧計の精度は約 1% です。バッテリーの電圧は約 4V です。同じことが、$max(dE_m) \approx 4*1\% = 0.04$ で表されます。したがって、$V = (max(dE_m))^2 \approx 1e-3$ と設定します。

  • サンプル時間: $Ts$。

バッテリーの劣化の推定

バッテリーの劣化は減少する容量 $C_q$ によってモデル化されます。この例では、劣化の影響を示すため、バッテリーの容量が放電/充電サイクルごとに 1 Ah 減少するように設定されています。容量の劣化率が事前にわかっていないので、$$ C_q $$ の状態方程式をランダム ウォークになるよう設定します。

$C_{q_{k+1}} = C_{q_k} + W_{C_q}$

ここで、$k$ は放電/充電サイクルの数、$W_{C_q}$ はプロセス ノイズです。

バッテリーは、バッテリーの SOC が 30% になると自動的に充電され、SOC が 90% になると放電に切り替わるように構成されています。この情報を使用して、電流 $I$ を 1 回の充電サイクルまたは放電サイクルで積分し、バッテリーの容量を測定します (クーロン カウンター)。

$C_q$ の測定方程式は次のとおりです。

$$ C_{q_k}^{Measured} = C_{q_k} + V_{C_q} =
\frac{\int_{t_{k-1}}^{t_k}{I}dt}{(\Delta SOC)_{nominal}}
=\frac{\int_{t_{k-1}}^{t_k}{I}dt}{|0.9-0.3|}=
\frac{\int_{t_{k-1}}^{t_k}{I}dt}{0.6}$$

ここで、$V_{C_q}$ は測定ノイズです。

バッテリーの劣化の状態と測定方程式は、次の状態空間形式にすることができます。

$$ C_{q_{k+1}} = A_{C_q} C_{q_k} + W_{C_q} $$

$$ C_{q_k}^{Measured} = C_{C_q} C_{q_k} + V_{C_q} $$

ここで、$A_{C_q}$ と $C_{C_q}$ は 1 に等しくなります。

上記の線形システムについて、カルマン フィルターを使ってバッテリーの容量を推定します。線形カルマン フィルターから推定された $C_q$ を使用して、SOC の推定を改善します。この例では、イベントベースの線形カルマン フィルターを使って $C_q$ を推定します。$C_q$ は充電サイクルまたは放電サイクルごとに 1 回測定されるので、線形カルマン フィルターが有効になるのは充電または放電の終了時のみです。

ブロック パラメーターとオプションを次のように設定します。

[モデル パラメーター] をクリックして、プラント モデルとノイズの特性を指定します。

  • モデルソース: Input Port。

イベントベースのカルマン フィルターを実装するため、イベントが発生した場合にのみ状態方程式が有効になります。つまり、状態方程式もイベントベースです。線形システム $x_{t+1} = Ax_t + Bu_t + w_t$ に対し、状態方程式を次のように設定します。

$$ x_{t+1} = \left \{ \begin{array}{cc} Ax_t + Bu_t + w_t, t=t_{enabled} \\ x_t, t \neq t_{enabled} \end{array} \right. $$.

  • A: $$ \left \{ \begin{array}{cc} A_{C_q}, t=t_{enabled} \\ 1, t \neq t_{enabled} \end{array} \right. $$。この例では、$A_{C_q}=1$ となります。その結果、$A$ は常に 1 に等しくなります。

  • C: 1。ここで、$C_{q_k}^{Measured} = C_{q_k} + V_{C_q} = \frac{\int_{t_{k-1}}^{t_k}{I}dt}{0.6}$。

  • 初期推定の [ソース]: Dialog。Initial state x[0] に初期状態を指定します。

  • 初期状態 x[0]: 30。これはバッテリーの定格容量 (30 Ah) です。

  • Q: $$ \left \{ \begin{array}{cc} 1, t=t_{enabled} \\ 0, t \neq t_{enabled} \end{array} \right. $$

これはプロセス ノイズ $W_{C_q}$ の共分散です。容量の低下率は放電/充電サイクルごとに約 1 Ah なので、プロセス ノイズを 1 に設定します。

  • R: 0.1。これは測定ノイズ $V_{C_q}$ の共分散です。容量の測定誤差は 1% 未満であると仮定します。バッテリーの容量が 30 Ah の場合、測定ノイズは $V_{C_q} \approx (0.3)^2 \approx 0.1$ です。

  • サンプル時間: Ts。

[オプション] をクリックし、測定値の更新を制御する入力端子 Enable を追加します。イネーブル端子は、バッテリー容量の推定値を連続して更新するのではなく、充電/放電イベントにおいて更新するために使用されます。

Enable を 0 に設定しても、状態方程式を使用した予測はディセーブルにならない点に注意してください。これが、状態方程式もイベントベースとなるよう設定される理由です。Kalman Filter ブロックにイベントベースの A と Q を設定することにより、Enable が 0 に設定されると、状態方程式を使用する予測がディセーブルになります。

結果

システムをシミュレートするには、バッテリーのパラメーターを読み込みます。このファイルには $E_m(SOC,T)$、$R_0(SOC,T)$、$R_1(SOC,T)$ などのバッテリー パラメーターが含まれています。

load BatteryParameters.mat

システムをシミュレートします。

sim('BatteryExampleUKF')

各タイム ステップで、アンセンテッド カルマン フィルターが電圧測定値 $E_m$ に基づいて SOC の推定値を提供します。実際の SOC、推定された SOC、およびその差異をプロットします。

% Synchronize two time series
[RealSOC, EstimatedSOC] = synchronize(RealSOC, EstimatedSOC, 'intersection');

figure;
subplot(2,1,1)
plot(100*RealSOC,'b','LineWidth',1.5);
hold on
plot(100*EstimatedSOC,'r--','LineWidth',1);
title('State of Charge');
xlabel('Time (s)');
ylabel('SOC (%)');
legend('Actual','UKF estimate','Location','Best','Orientation','horizontal');
axis tight

subplot(2,1,2)
DiffSOC = 100*(RealSOC - EstimatedSOC);
plot(DiffSOC.Time, DiffSOC.Data, 'LineWidth', 1.5);
xlabel('Time(s)');
ylabel('\Delta SOC (%)','Interpreter','Tex');
legend('Difference between Real SOC and Estimated SOC','Location','Best')
axis tight

初期推定誤差の後、SOC は実際の SOC へと迅速に収束します。最終の推定誤差は 0.5% 内です。アンセンテッド カルマン フィルターにより SOC を正確に推定できます。

放電/充電の各遷移において、SOC の推定を改善するためにバッテリーの容量が推定されます。バッテリー システムはインジケーター信号を出力してバッテリーの現在の処理を通知します。放電処理はインジケーター信号で -1 として表され、充電処理は 1 で表されます。この例では、容量の推定にカルマン フィルターをイネーブルにするかディセーブルにするかを判定するためにインジケーター信号の変化が使用されます。実際の容量と推定された容量、および充電/放電のインジケーター信号をプロットします。

figure;
subplot(2,1,1);
plot(RealCapacity,'b','LineWidth',1.5);
hold on
plot(EstimatedCapacity,'r--','LineWidth',1.5);
xlabel('Time (s)');
ylabel('Capacity (Ah)');
legend('Actual','KF estimate','Location','Best');

subplot(2,1,2);
plot(DischargeChargeIndicator.Time,DischargeChargeIndicator.Data,'b','LineWidth',1.5);
xlabel('Time(s)');
ylabel('Indicator Signal');

一般に、カルマン フィルターを使用して実際の容量を追跡することが可能です。容量の推定値と実際の値の間には 1/2 サイクルの遅延があります。この遅延は、測定のタイミングに起因します。バッテリーの容量低下の計算は、完全な放電/充電サイクルが 1 度終了したときに行われます。クーロン カウンターにより、前回の放電または充電サイクルに基づき容量測定が行われます。

まとめ

この例では、Simulink の Unscented Kalman Filter ブロックを使用してリチウム電池の非線形の状態推定を行う方法を説明しました。また、バッテリーの容量を推定するためのイベントベースのカルマン フィルターを開発する手順も示しています。新しく推定された容量は、アンセンテッド カルマン フィルターでの SOC の推定を改善するために使用されます。

参考文献

[1] Huria, Tarun, et al."High fidelity electrical model with thermal dependence for characterization and simulation of high power lithium battery cells."Electric Vehicle Conference (IEVC), 2012 IEEE International.IEEE, 2012.

[2] Wan, Eric A., and Rudolph Van Der Merwe."The unscented Kalman filter for nonlinear estimation."Adaptive Systems for Signal Processing, Communications, and Control Symposium 2000.AS-SPCC.Ieee, 2000.

関連するトピック