# PEMWE: 0D + Step1b.1 統合ベイズ状態空間モデル レポート
## 1. 目的
本レポートの目的は、KIT の incremental PEMWE 実測 EIS データ（Figure 5 の 4 条件）を用いて、
1. Step1b.1 で得られる **EIS 等価回路パラメータ**を、
2. 0D の **電圧分解モデル**へ接続し、
3. さらにそれらを **ベイズ状態空間モデル**として統合する
ことです。
狙いは、EIS の各スペクトルを個別にフィッティングして終わるのではなく、
**「EIS で見えている抵抗・緩和を、セル電圧と両立する latent health state に変換する」**ことにあります。
## 2. 今回の入力データ
今回使用したデータは、ユーザー提供の KIT raw txt 4 本です。
- `Fig5_750mAcm-2.txt`
- `Fig5_1000mAcm-2.txt`
- `Fig5_1500mAcm-2.txt`
- `Fig5_2000mAcm-2.txt`
各 txt には `Frequency/Hz`, `Impedance/Ohm`, `Phase/deg`, およびヘッダ内の `Potential`, `Current` が含まれます。ここから `Zre`, `Zim` を復元し、Step1b.1 の回路でフィッティングしました。
## 3. モデルの全体構造
今回の統合モデルは 2 層です。
### 3.1 Step1b.1: EIS 観測モデル
各スペクトルに対して
$$
Z(\omega) = j\omega L_{par} + R_{ohm} + ZARC_1 + ZARC_2
$$
を当てています。
ここで `L_par` は高周波の寄生インダクタンス、`R_ohm` は高周波オーム抵抗、`ZARC_1`, `ZARC_2` は非理想 semicircle です。
ただし `ZARC_1`, `ZARC_2` をそのまま物理機構に一対一対応させず、**時定数 `tau` の大小で HF / LF に並べ替えた上で、操作的 proxy として使う**方針を採っています。
## 4. なぜ `R_HF` 単独ではなく `R_non = R_HF + R_LF` を使うのか
今回の 4 条件では、HF / LF の厳密な分離は条件によって安定性が異なりました。特に 1000 mA cm⁻² 付近では、2-ZARC 分解をそのまま kinetics / transport に一対一対応させると HF arc が極端に小さくなり、解釈が不安定になります。
そこで、統合モデルでは
- `ASR = R_ohm * A`
- `R_non = (R_HF + R_LF) * A`
- `R_LF = R_LF * A`
を使いました。
この設計の意味は、`ASR` をオーム損失、`R_non` を全非オーム損失（主に kinetics 側）、`R_LF` を低周波輸送寄与の proxy として使うことです。つまり、**kinetics は total non-ohmic で表し、transport は LF branch を追加的に参照する**構造です。
## 5. 0D 電圧モデル
統合モデルの 0D 側は、次の compact な電圧式です。
$$
V_t = v_{base} + \eta_{act,t} + \eta_{ohm,t} + \eta_{diff,t}
$$
各項は以下です。
### 5.1 オーム損失
$$
\eta_{ohm,t} = j_t \cdot ASR_t
$$
### 5.2 活性化損失
$$
\eta_{act,t} = B_{act}\,\mathrm{asinh}\left(\frac{j_t}{2i_{0,eff,t}}\right)
$$
ただし
$$
i_{0,eff,t} = \frac{\kappa_{act}}{R_{non,t}}
$$
と置いています。`R_non` が大きいほど実効交換電流密度 `i0_eff` は小さくなり、活性化損失は大きくなります。
### 5.3 輸送損失
$$
\eta_{diff,t} = B_{diff}\ln\left(\frac{i_{L,eff,t}}{i_{L,eff,t}-j_t}\right)
$$
ただし
$$
i_{L,eff,t} = j_t + \frac{\kappa_{diff}}{R_{LF,t}}
$$
と置いています。`R_LF` が大きいほど `iL_eff` が operating current に近づき、輸送損失が増えます。
### 5.4 `v_base`
`v_base` は、理想可逆電圧と、ここでは明示分離していない基準項をまとめた lumped baseline です。今回の Figure 5 データは 4 operating points のみで、温度・圧力の時系列情報がないため、まずは `v_base` として吸収し、将来の熱収支付き 0D へ拡張できる形にしています。
## 6. 状態空間モデル
latent state は各 operating point ごとに
$$
x_t = [\log ASR_t,\; \log R_{non,t},\; \log R_{LF,t}]^\top
$$
です。
### 6.1 状態遷移
今回は真の長期劣化時系列ではなく、同日の 4 operating points しかないため、最小構成として **local level random walk** を使いました。
$$
x_t \sim \mathcal{N}(x_{t-1}, Q)
$$
ここで `Q = diag(q_ASR, q_Rnon, q_RLF)` です。
重要なのは、この `q` は **今のデータでは「劣化速度」ではなく、「状態の滑らかさ」**を表していることです。将来、真の PoC 時系列（起動停止・累積通電量・水質・圧力差など）を入れる段階で、`x_t = x_(t-1) + B s_t + η_t` へ自然に拡張できます。
### 6.2 観測方程式
EIS 観測は log-space で
$$
\log y^{(EIS)}_t \sim \mathcal{N}(x_t, \Sigma_t)
$$
とし、`Σ_t` は Step1b.1 の posterior SD から近似した log-space の観測分散です。電圧観測は
$$
V^{obs}_t \sim \mathcal{N}(V_{0D}(j_t, x_t), \sigma_V^2)
$$
としています。
## 7. 推定手順
今回の integrated layer は **完全な HMC ではなく MAP + Laplace 近似**です。
理由は、(i) Step1b.1 自体はすでにベイズ単時点フィットで不確実性を持っている、(ii) 今回の統合対象は 4 点のみ、(iii) まずは構造を明示し、latent state と電圧式の整合を見る段階だからです。
実装としては、
1. MAP 推定（`scipy.optimize.minimize`）
2. MAP 近傍の有限差分 Hessian
3. Hessian 逆行列を posterior covariance とみなす
4. そこから多変量正規近似で posterior draws を生成
としています。
## 8. 結果
### 8.1 Stage1b.1 から抽出した feature
| file | j [A cm^-2] | V [V] | ASR median [Ω cm²] | R_non median [Ω cm²] | R_LF median [Ω cm²] |
|---|---:|---:|---:|---:|---:|
| Fig5_750mAcm-2.txt | 0.75 | 1.6615 | 0.1243 | 0.0434 | 0.0137 |
| Fig5_1000mAcm-2.txt | 1.00 | 1.7015 | 0.1243 | 0.0337 | 0.0335 |
| Fig5_1500mAcm-2.txt | 1.50 | 1.7790 | 0.1246 | 0.0264 | 0.0070 |
| Fig5_2000mAcm-2.txt | 2.00 | 1.8540 | 0.1264 | 0.0231 | 0.0040 |
ASR はおおむね **0.124–0.126 Ω cm²** でかなり安定です。一方で `R_non` は **0.044 → 0.023 Ω cm²** 程度へ低下し、高電流密度ほど total non-ohmic 抵抗が小さく見えています。`R_LF` は 1000 mA cm⁻² で一時的に大きく、その後 1500–2000 mA cm⁻² で小さくなっています。これは LF branch の解釈が operating point に依存して動くことを示しており、**HF/LF の生 split をそのまま物理量と断定しない**という今回の設計判断を支持します。
### 8.2 latent state の平滑化
この図では、× が Step1b.1 観測量、線が統合状態空間モデルの latent state median、帯が 90% credible band です。ASR はもともと安定で、state-space を通しても大きくは動きません。いっぽう `R_non`, `R_LF` は、EIS 単独ではややばらつきがあるものの、電圧観測を同時に使うことで、**セル電圧と整合する範囲へ再配置**されています。
### 8.3 電圧 fit と分解
この図では、measured voltage、posterior predictive band、posterior median、および `v_base`, `v_base + η_act`, `v_base + η_act + η_ohm` を重ねています。今回の 4 点データでは posterior median はほぼ観測点を通っています。ただしこれは「4 点に 0D を当てたから当然ぴったり」というより、**EIS 観測で拘束された latent state のもとで 0D 電圧式が矛盾なく成立する**ことの確認として読むべきです。
### 8.4 implied effective states
この図は、統合モデルが暗黙に持つ `i0_eff` と `iL_eff` の posterior を示します。ここで重要なのは、**これらの絶対値は calibration-dependent** だということです。今回の `κ_act`, `κ_diff` は compact proxy の scale parameter なので、この `i0_eff`, `iL_eff` は“真の材料定数”というより **EIS とセル電圧を両立させる実効状態量**です。
## 9. グローバル・パラメータの結果
| parameter | median | q05 | q95 |
|---|---:|---:|---:|
| v_base | 1.43165 | 1.34391 | 1.52345 |
| kappa_act | 0.00324919 | 0.000745111 | 0.0147433 |
| kappa_diff | 0.108293 | 0.0252385 | 0.435654 |
| q_asr | 0.033301 | 0.00869816 | 0.119183 |
| q_rnon | 0.200952 | 0.0934225 | 0.433643 |
| q_rlf | 0.655766 | 0.435852 | 0.96582 |
| sigma_v | 0.00013856 | 1.79599e-05 | 0.00103577 |
`v_base` は約 1.432 V、`q_asr < q_rnon < q_rlf` で、ASR は安定、LF 側は最も変動しやすい、という傾向です。`sigma_v` は約 1.39e-04 V で、今回の compact 0D + latent-state coupling で電圧残差はかなり小さいです。
## 10. 物理的解釈
### 10.1 何が比較的強く言えるか
- `ASR` は、膜＋接触＋電子経路を含んだ **高周波オーム損失**の proxy としてかなり安定
- `R_non` は、主として **非オーム損失の総量**として扱うのが妥当
- `R_LF` は、**輸送・気泡・水管理寄りの遅い過程**を含む proxy と読むのが自然
- `L_par` はセル内部ではなく **測定系寄生**として扱うべき
### 10.2 何をまだ言い切れないか
- HF / LF の 2 arc を、それぞれ単一の素過程へ一対一対応させること
- `i0_eff`, `iL_eff` を材料固有の真値とみなすこと
- 今回の `q` をそのまま degradation rate とみなすこと
## 11. なぜこの integrated model が有効か
このモデルの本質は、**「EIS をそのまま電圧に入れた」のではなく、EIS を latent health state に変換してから 0D と結合した」**ところにあります。Step1b.1 は観測モデル、0D は状態モデルです。
- Step1b.1: スペクトル → 抵抗・緩和の proxy
- 0D: proxy → 電圧・実効速度論・実効輸送制約
- state-space: operating point 間で一貫した latent state として平滑化
この役割分担があるので、将来、真の時系列 EIS、起動停止回数、累積通電量、水質、温度、圧力差、crossover / gas purity を入れても自然に拡張できます。
## 12. 限界と次の一手
### 限界
1. 今回のデータは 4 operating points で、真の degradation 時系列ではない
2. integrated layer は Laplace 近似であり、HMC ではない
3. 温度・圧力は明示的な状態としては入れていない
4. `η_par` は baseline へ吸収しており、別観測で分離していない
### 次の一手
1. **真の longitudinal EIS + polarization 時系列**へ適用
2. 状態遷移を `x_t = x_(t-1) + B s_t + η_t` へ拡張
3. `T_t` を latent state に加えて 0D 熱収支を入れる
4. `R_non` を 0D / 1D 物理モデルと結び、`ECSA_eff`, `κ_eff`, `iL_eff` の同定へ進む
5. 将来的には HMC / particle MCMC へ置換