【Pythonでベイズ推定】PyStanでリチウムイオン電池を劣化推定

スポンサーリンク

こんにちは。電池研究のかたわら、ベイズ統計やPythonの勉強をしています。
今回は、ベイズ推定の計算を簡単に実装できる確率的プログラミング言語の一つであるStanを使った回帰モデルの実装例を紹介します。Stanは、PythonやR、Matlabなどの環境から動かすことができるようになっており、私はPython環境で使っています。

ただの直線を予測しても面白くないので、今回はリチウムイオン電池の劣化推移をイメージした非線形なデータを例題として扱います。
ベイズ的な回帰モデルとすることでデータの不確かさを加味しやすいので、「データが少ないけど予測したい」「予測がどの程度ブレる可能性があるのか知りたい」といった場合に有用です。

①リチウムイオン電池の劣化モデル【ルート則】

リチウムイオン電池の劣化推移を予測する一般的なモデルとして、ルート則というモデルが一般的です。
これは、リチウムイオン電池の容量低下が、時間の√に比例するというモデルです。実際、特定の劣化メカニズムは、確かにルート則に従うことは知られています。(例えばこちらの文献:負極表面のSEI被膜の成長は固相の一次元拡散が律速になっているため、被膜成長速度が√(時間)に比例することを実験的に確認した文献です。)
https://www.gs-yuasa.com/jp/technology/technical_report/pdf/vol10_2/010_02_008.pdf

ルート則に従うこの劣化は、リチウムイオン電池の最も一般的な(かつ理想的な)劣化メカニズムのため、ルート則による寿命予測とは「電池が理想的な劣化推移を辿ると仮定した時の予測」と言い換えることができます。
なお、ルート則の拡張版として、べき乗則もよく用いられます。ルート則は時間の0.5乗とも表すことができますが、このべき定数0.5を変数としたモデルです。劣化前の電池の容量yを100として、使用時間に伴って容量が低下していくとすると、
【ルート則】
$$ y = 100 \ – \ α\sqrt{x} \ = \ 100 \ – \ α \cdot x^{0.5}$$
【べき乗則】
$$y = \ 100 \ – \ α \cdot x^{β}$$
例えばべき定数β=0.3で、劣化が非常に緩やか場合や、べき定数β=1.0で直線的に容量が低下してしまう場合など、「電池が理想的な劣化推移」とならない様子を表現することができます。

○検証に用いるデータ

今回はルート則にノイズを加えたシミュレーションデータをNumpyで生成します。なお、定数は適当です。
予測の際には、べき乗数(真値 0.5)も未知パラメータとしたべき乗則モデルでフィッティングしてみたいと思います。

[Python code]

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import pystan
%matplotlib inline

##適当なパラメータを生成
x = np.arange(0.5, 5., 0.5)
alpha = 5.
beta = 0.5
sigma = 0.2

##初期容量を100とし、べき乗則に従って容量が減衰していくモデルに正規分布に従うノイズを加える
y_obs = 100 - alpha * x ** beta + np.random.normal(0, sigma, len(x))

##可視化
fig, ax1 = plt.subplots(figsize=(7, 5))
ax1.scatter(x, y_obs)
plt.xlim(0, 10)
plt.ylim(80, 100)

②Stanによる基本の実装

基本的には、線形モデルの実装と大差はありません。data{ }、parameters{ } でそれぞれインプットデータと推定したいパラメータを定義し、model{ }に仮定したモデルを定義します。tranceformed parameters{ }など、その他のブロックも必要に応じて使います。今回は、tranceformed parameters{ }ブロックでxのべき乗値(pow_X)を生成しました。

なお、データやパラメータが配列の場合(例: [x1, x2, x3 …xn])、n個の数値の配列(real x[n])ではなくvector型(vector[n] X)として定義すると、確率分布からサンプリングする際の計算がシンプルになるため高速化できるとのことです。またコードが見やすくなるものメリットです。(Stanのマニュアル:ループ処理が遅いということではないそうです。)

べき乗の計算($x^β$の部分)をしているpow()関数には、ベクトルには対応していないようだったのでループ処理をしています。(ここはただの代入なので、ベクトル化できたとしても速度は変わらないんですかね。)

stan_model = """
data {
   int<lower=1> N;
   vector[N] X_obs;
   vector[N] Y_obs;
}


parameters {
   real<lower=0> alpha;
   real<lower=0> beta;
   real<lower=0> sigma;
}

transformed parameters {
   vector[N] pow_X;
   for (i in 1:N) pow_X[i] = pow(X_obs[i], beta);
}


model {
   Y_obs ~ normal(100 - alpha * pow_X , sigma);
}
"""
sm = pystan.StanModel(model_code=stan_model)

モデルを定義したら、コンパイルします。ここは少し時間がかかるので待ちましょう。

モデルに渡すデータをPythonの辞書型で定義します。
コンパイル済みのモデルにデータを渡して解析を実行します。繰り返し回数は3000回、そのうち最初の500回をウォームアップとして棄却します。これは、マルコフ連鎖モンテカルロ法を用いてパラメータ推定しているためで、乱数を発生させて、パラメータを少しずつ変化させながら、データとモデルに適合するパラメータを取得していきます。このような逐次計算を、初期値を変えて3回行い(chains=3の部分)、それらを集計します。それぞれのパラメータについて、繰り返し回数分((3000-500)×3 = 7500)の数値が取得され、この数値の分布がベイズ推定したパラメータの確率分布(事後分布)になります。
モデルに問題がなければ、数秒で計算が終わるはずです。

stan_data = {
    "N":len(x),
    "X_obs":x,
    "Y_obs":y_obs
}

fit = sm.sampling(data = stan_data, iter=3000, warmup=500, chains=3)

print()関数で計算結果を表示してみます。
Rhatが1.0となっていることから、計算が収束していることがわかります。(pow_X[2]は、「x=1.0のべき乗は、βの値に依らず常に1.0だから、ばらつきは計算できないよ!」となりnanと表示されていますが問題ありません。)

print(fit)
##以下は出力結果
"""
Inference for Stan model: anon_model_528e88c91936a01c8ba5483b18b51d5a.
3 chains, each with iter=3000; warmup=500; thin=1; 
post-warmup draws per chain=2500, total post-warmup draws=7500.

           mean se_mean     sd   2.5%    25%    50%    75%  97.5%  n_eff   Rhat
alpha      5.01  3.5e-3   0.15   4.71   4.93   5.01    5.1   5.31   1830    1.0
beta       0.49  5.8e-4   0.03   0.44   0.47   0.49    0.5   0.54   1899    1.0
sigma      0.26  2.2e-3    0.1   0.15    0.2   0.24    0.3    0.5   1985    1.0
pow_X[1]   0.71  2.9e-4   0.01   0.69   0.71   0.71   0.72   0.74   1907    1.0
pow_X[2]    1.0     nan    0.0    1.0    1.0    1.0    1.0    1.0    nan    nan
pow_X[3]   1.22  2.9e-4   0.01   1.19   1.21   1.22   1.23   1.24   1894    1.0
pow_X[4]    1.4  5.7e-4   0.02   1.35   1.39    1.4   1.42   1.45   1891    1.0
pow_X[5]   1.56  8.4e-4   0.04   1.49   1.54   1.56   1.58   1.64   1888    1.0
pow_X[6]   1.71  1.1e-3   0.05   1.62   1.68   1.71   1.73   1.81   1886    1.0
pow_X[7]   1.84  1.4e-3   0.06   1.73   1.81   1.84   1.87   1.97   1884    1.0
pow_X[8]   1.96  1.6e-3   0.07   1.83   1.92   1.96    2.0   2.11   1882    1.0
pow_X[9]   2.08  1.8e-3   0.08   1.93   2.03   2.08   2.12   2.25   1881    1.0
lp__       7.89    0.04   1.55    3.9   7.21    8.3   9.02    9.6   1616    1.0

Samples were drawn using NUTS at Tue Oct 13 17:09:56 2020.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).
"""

③予測分布のグラフ化

推定したα、βを用いて電池劣化推移をグラフ化してみます。α、βが確率分布を表す、7500個のサンプルであることから、計算される容量Yも7500個の数値になります。
Numpyで計算したあと、7500個のサンプルの50%ベイズ信頼区間、95%ベイズ信頼区間を計算し、グラフに描写します。なお、信頼区間の計算にはScipyを使いました。


x_sim = np.arange(0, 20, 0.1)
a = fit['alpha']
b = fit['beta']

y_calc = pd.DataFrame()
for i in np.arange(len(x_sim)):
    y_calc[i] = 100 - a*x_sim[i]**b

from scipy.stats import mstats
low_y50, high_y50 = mstats.mquantiles(y_calc, [0.25, 0.75], axis=0)   ##配列y_calcから50%ベイズ信頼区間を計算
low_y95, high_y95 = mstats.mquantiles(y_calc, [0.025, 0.975], axis=0)  ##配列y_calcから95%ベイズ信頼区間を計算

fig, ax1 = plt.subplots(figsize=(7,5))
ax1.scatter(x, y_obs, color = "r", label="observed") ##測定データ
ax1.plot(x_sim, y_calc.mean(), linestyle="dashed", color="black", label = "mean") ##推定した分布の中央値
ax1.fill_between(x_sim, low_y50, high_y50, alpha=0.6, color="darkgray", label = "50%CI")
ax1.fill_between(x_sim, low_y95, high_y95, alpha=0.3, color="gray", label = "95%CI")

ax1.legend()
plt.xlim(0, 10)
plt.ylim(80, 100)

④シミュレーションの実装【generated quantities{}ブロック】

上記では、容量の予測値YをPythonで計算しましたが、Stanのモデル内で計算することもできます。また正規分布に従う測定誤差 Normal(0, sigma)を考慮したサンプリングも行えます。
改変したモデルは下記です。


stan_model2 = """
data {
   int<lower=1> N;
   vector[N] X_obs;
   vector[N] Y_obs;
   int<lower=1> N_sim;
   vector[N_sim] X_sim;
}


parameters {
   real<lower=0> alpha;
   real<lower=0> beta;
   real<lower=0> sigma;
}


model {
   vector[N] pow_X;
   for (i in 1:N) pow_X[i] = pow(X_obs[i], beta);
   Y_obs ~ normal(100 - alpha * pow_X , sigma);
}

generated quantities {
   vector[N_sim] Mu;
   vector[N_sim] Y_calc;
   {
      vector[N_sim] pow_X_sim;
      for (i in 1:N_sim) {
         pow_X_sim[i] = pow(X_sim[i], beta);
         Mu[i] = 100 - alpha * pow_X_sim[i];
         Y_calc[i] = normal_rng(Mu[i], sigma);
      }
   }
}

"""

なお、最初のモデルではxのべき乗値pow(x)をtransformed parameters {}ブロックで生成しましたが、今回はmodel{}ブロックで生成させました。
こうすることで、計算結果の出力時に値がpow(X)が出力されなくなります。

解析を実行するPythonコードは基本的に同様ですが、今回はシミュレーションデータもStan内で生成するため、シミュレーションするx値もStanに渡します。


sm2 = pystan.StanModel(model_code=stan_model2)

x_sim = np.arange(0.1, 10, 0.1)

stan_data = {
    "N":len(x),
    "X_obs":x,
    "Y_obs":y_obs,
    "N_sim":len(x_sim),
    "X_sim":x_sim
}

fit2 = sm2.sampling(data = stan_data, iter=3000, warmup=500, chains=3)

モデルに含めていない将来のデータ(とはいえ仮想データですが)とともに、予測結果を表示してみます。


y_calc = fit2['Y_calc']
low_y50, high_y50 = mstats.mquantiles(y_calc, [0.25, 0.75], axis=0)
low_y95, high_y95 = mstats.mquantiles(y_calc, [0.025, 0.975], axis=0)

x_2 = np.arange(5, 20, 0.5)

alpha = 5.
beta = 0.5
sigma = 0.2

y_2 = 100 - alpha*x_2**beta + np.random.normal(0, sigma, len(x_2))

fig, ax1 = plt.subplots(figsize=(7,5))
ax1.scatter(x, y_obs, label="Observed") ##実測データ(Stanに渡した分)
ax1.scatter(x_2, y_2, label="Observed") ##実測データ(将来のデータ)

ax1.plot(x_sim, np.median(y_calc, axis=0), linestyle="dashed", color="black", label = "median")
ax1.fill_between(x_sim, low_y50, high_y50, alpha=0.5, color="darkgray", label = "50%CI")
ax1.fill_between(x_sim, low_y95, high_y95, alpha=0.3, color="gray", label = "95%CI")

ax1.legend()
plt.xlim(0, 10)
plt.ylim(80, 100)

ベイズ推定を使えば、データの不確かさを考慮した上で、予測結果を確率分布として得られるため、例えば「データが少ないけど予測したい」「予測がどの程度ブレる可能性があるのか知りたい」といった際に有用です。
ただし、あくまで解析者が仮定したモデルを前提としているため、誤ったモデルを仮定すれば誤った予測分布になりうるという点には注意が必要です。

また、今回は用いませんでしたが、パラメータに関する事前知識がある場合、その事前知識を確率分布(事前分布)としてモデルに反映することができます。(例えば、「過去のデータからすると、αは3.0〜7.0くらいの間には収まるだろう」、など。緩い制約条件のようなものと考えて良いでしょうか。)
事前分布を応用した、階層ベイズモデルというのもあります。そのうち書きたいと思います。

今回は以上となります!

タイトルとURLをコピーしました