파라미터 추정과 하이브리드 모델링

파라미터 추정

왜 필요한가. 같은 모델의 셀이라도 내부 저항은 제조 편차로 몇 % 다르고, 온도가 내려가면 몇 배로 커지며, 노화하면서 서서히 증가한다. 충전 상태(state of charge, SOC) 추정이나 출력 한계 계산은 이 값을 그대로 쓰므로, 파라미터가 틀리면 그 위에 쌓은 모든 계산이 틀린다. 파라미터의 변화 자체가 노화의 지표로 쓰이기도 한다.

모델. 30.2장의 1-RC 등가회로 모델(equivalent circuit model, ECM)은 개방 전압(open-circuit voltage, OCV) 원천, 직렬 저항 $R_0$, 그리고 저항 $R_1$과 커패시터 $C_1$의 병렬 쌍으로 이루어진다. 방전 전류를 양수로 두면 단자 전압은

이다. $R_0$는 전해질·집전체 저항처럼 전류 변화에 즉시 반응하는 부분이고, RC 쌍은 전하 이동과 확산처럼 시상수 $\tau = R_1C_1$로 천천히 쌓이고 풀리는 과전압을 흉내 낸다. 전류가 한 샘플 간격 $\Delta t$ 동안 일정하면(영차 유지, zero-order hold) $V_1$은 정확히

로 갱신된다. OCV 곡선은 아주 낮은 전류로 따로 측정했다고 가정한다.

펄스 응답에서 읽는 것. 휴지 상태의 셀에 전류 펄스를 걸면 전압 응답이 세 부분으로 나뉜다. 펄스 시작 순간 전압이 $I R_0$만큼 즉시 떨어지고, 이어서 RC 쌍이 시상수 $\tau$로 $I R_1$까지 더 떨어지며, 펄스가 끝나면 같은 순서로 회복된다. 즉시 변화에서 $R_0$, 천천히 변하는 부분의 크기에서 $R_1$, 그 속도에서 $\tau$를 읽는다. 아래 위젯에서 세 슬라이더를 움직이며 응답의 어느 부분이 어느 파라미터에 반응하는지 확인한다.

비선형 최소제곱. 손으로 읽는 대신 측정 전체를 한꺼번에 맞춘다. 측정 전압 $y_k$와 모델 전압 $V_k(\boldsymbol\theta)$의 잔차 제곱합

을 최소화한다. 잡음이 독립 가우시안이면 이것이 최우추정이다(12.2장). 세 파라미터는 단위와 크기가 크게 달라서($R_0 \approx 0.02\ \Omega$, $C_1 \approx 2000$ F) 원래 값으로 최적화하면 척도가 나쁘다. 로그 파라미터로 최적화하면 척도가 상대 변화로 통일되고 양수 조건도 저절로 지켜진다.

가상의 셀(용량 3 Ah, $R_0 = 20$ mΩ, $R_1 = 15$ mΩ, $C_1 = 2$ kF, $\tau = 30$ s)에 3 A 방전 펄스를 300 s 걸고 앞뒤로 휴지를 둔 10분짜리 실험을 1 Hz로 기록한다. 측정 잡음은 1 mV다. 파라미터 값은 설명용이고 특정 셀의 측정값이 아니다.

import numpy as np
from scipy.optimize import least_squares

Q_AH = 3.0                                           # 용량 (설명용 값)
def ocv_v(soc):                                      # 개방 전압 곡선 (설명용, 저율 시험으로 따로 쟀다고 가정)
    return 3.5 + 0.55 * soc + 0.04 * np.log((soc + 0.02) / (1.02 - soc))

def simulate_v(theta, t_s, i_a, soc0=0.8):           # 1-RC ECM, I > 0이 방전. theta = (R0 [Ω], R1 [Ω], C1 [F])
    r0_ohm, r1_ohm, c1_f = theta                     # 스칼라 또는 (G,) 배열(격자 계산용)
    dt_s = t_s[1] - t_s[0]
    a = np.exp(-dt_s / (r1_ohm * c1_f))              # step 동안 전류가 일정하면(영차 유지) 정확한 갱신
    soc = soc0 - np.concatenate([[0], np.cumsum(i_a[:-1])]) * dt_s / (3600 * Q_AH)   # (T,)
    v1, v1_hist = 0.0 * a, []
    for i_k in i_a:
        v1_hist.append(v1)
        v1 = a * v1 + r1_ohm * (1 - a) * i_k
    return ocv_v(soc) - i_a * r0_ohm - np.stack(v1_hist, axis=-1)   # (T,) 또는 (G, T)

t_s = np.arange(0, 600.0)                            # (600,) 1 Hz 샘플
i_a = np.where((t_s >= 60) & (t_s < 360), 3.0, 0.0)  # 3 A 방전 펄스 300 s, 앞뒤 휴지
theta_true = np.array([0.020, 0.015, 2000.0])        # R0 20 mΩ, R1 15 mΩ, C1 2 kF → τ = 30 s
rng = np.random.default_rng(0)
v_meas = simulate_v(theta_true, t_s, i_a) + rng.normal(0, 1e-3, t_s.shape)   # 측정 잡음 1 mV

def residual(log_theta):                             # 로그 파라미터: 양수 보장, 척도가 다른 파라미터를 상대 변화로
    return simulate_v(np.exp(log_theta), t_s, i_a) - v_meas
fit = least_squares(residual, np.log([0.01, 0.01, 1000.0]))
r0, r1, c1 = np.exp(fit.x)
print(f"추정: R0 = {r0*1e3:.2f} mΩ, R1 = {r1*1e3:.2f} mΩ, C1 = {c1:.0f} F, τ = {r1*c1:.1f} s")
print(f"잔차 RMS = {np.sqrt(np.mean(fit.fun**2))*1e3:.2f} mV (측정 잡음 1 mV)")

추정: R0 = 19.95 mΩ, R1 = 15.07 mΩ, C1 = 1987 F, τ = 30.0 s
잔차 RMS = 1.00 mV (측정 잡음 1 mV)

세 파라미터가 1% 안쪽으로 맞고, 잔차의 RMS가 측정 잡음과 같은 1 mV다. 잔차가 잡음 수준까지 내려갔다는 것은 모델이 데이터를 설명하는 데 필요한 구조를 갖췄다는 필요조건이다. 충분조건은 아니다. 다음 코드가 그 이유를 보인다.

식별 가능성은 실험이 정한다. 24.6장에서 본 것처럼, 추정값을 해석하려면 그 값이 데이터로 정해져야 한다. 참값 근처에서 로그 파라미터에 대한 민감도 행렬 $\mathbf{S}$(측정 시각 × 파라미터)를 계산하면 근사 공분산 $\sigma^2(\mathbf{S}^\top\mathbf{S})^{-1}$에서 각 파라미터의 상대 표준오차와 파라미터 사이의 상관을 얻는다. 같은 모델에 세 실험 설계를 비교한다.

def sensitivity(t_s, i_a, idx):                      # 참값에서 로그 파라미터 민감도 행렬 (len(idx), 3)
    f = lambda lt: simulate_v(np.exp(lt), t_s, i_a)[idx]
    lt0, h = np.log(theta_true), 1e-6
    return np.column_stack([(f(lt0 + h * e) - f(lt0 - h * e)) / (2 * h) for e in np.eye(3)])

t_short = np.arange(0, 24.0)                         # 펄스 3 s만, 휴지 없음
i_short = np.where(t_short >= 20, 3.0, 0.0)
designs = {"1 Hz, 펄스 300 s + 휴지": (t_s, i_a, np.arange(600)),
           "같은 실험, 60 s 간격 샘플": (t_s, i_a, np.arange(30, 600, 60)),
           "1 Hz, 펄스 3 s, 휴지 없음": (t_short, i_short, np.arange(24))}
for name, (tt, ii, idx) in designs.items():
    S = sensitivity(tt, ii, idx)
    cov = 1e-3**2 * np.linalg.inv(S.T @ S)           # 잡음 1 mV에서의 근사 공분산 (로그 파라미터)
    sd = np.sqrt(np.diag(cov)); corr = cov / np.outer(sd, sd)
    print(f"{name:18s} 상대오차 R0 {sd[0]:6.3f} R1 {sd[1]:6.3f} C1 {sd[2]:6.3f}"
          f" | 상관 R0–R1 {corr[0, 1]:+.3f}, R1–C1 {corr[1, 2]:+.3f}")

1 Hz, 펄스 300 s + 휴지 상대오차 R0  0.004 R1  0.006 C1  0.013 | 상관 R0–R1 -0.968, R1–C1 -0.865
같은 실험, 60 s 간격 샘플  상대오차 R0  0.132 R1  0.172 C1  0.335 | 상관 R0–R1 -0.998, R1–C1 -0.991
1 Hz, 펄스 3 s, 휴지 없음 상대오차 R0  0.016 R1 19.999 C1  1.077 | 상관 R0–R1 +0.494, R1–C1 +0.952

• 1 Hz, 펄스 300 s + 휴지. 세 파라미터 모두 상대오차가 1.3% 이하다. 그래도 $R_0$와 $R_1$의 상관이 −0.968로 강하다. 펄스 동안 전압이 머무는 높이는 합 $R_0 + R_1$로 정해지고, 둘을 나누는 정보는 펄스 시작과 끝의 몇 초에만 있기 때문이다.
• 같은 실험, 60 s 간격 기록. 실험은 같은데 데이터 로거가 1분에 한 번만 기록했다. $\tau = 30$ s보다 기록 간격이 길어 즉시 변화와 느린 변화를 구분할 샘플이 거의 없다. $R_0$–$R_1$ 상관이 −0.998이 되고, 상대오차가 13–34%로 커진다. 데이터가 정하는 것은 사실상 합 $R_0 + R_1$뿐이다.
• 1 Hz, 펄스 3 s, 휴지 없음. 펄스가 $\tau$보다 훨씬 짧으면 $V_1 \approx I t / C_1$로 거의 직선이라 $R_1$이 거의 나타나지 않는다. $R_1$의 상대오차가 약 20, 즉 2000%다. 이 실험으로는 $R_1$을 사실상 알 수 없다.

잔차가 작아도 식별되지 않을 수 있다. 셋째 설계에서 최적화를 돌리면 잔차는 잡음 수준까지 내려가고 $R_1$에는 어떤 숫자가 나온다. 그 숫자는 초기값과 최적화기가 멈춘 위치가 정한 것이다. 그래서 파라미터 추정 결과는 항상 불확실성(또는 상관)과 함께 보고하고, 상관이 ±0.99에 가까운 쌍은 개별 값이 아니라 정해진 조합(여기서는 $R_0 + R_1$)만 해석한다. 민감도 행렬은 참값 근처의 국소 선형화이므로, 비선형성이 강하면 24.6장의 프로파일 우도로 확인한다.

베이지안 추정 맛보기. 같은 문제를 사후분포로 풀면 "데이터가 무엇을 정했는가"가 더 직접 보인다. 같은 모델의 다른 셀들에서 얻은 값으로 사전분포를 둔다. $R_1$은 중앙값 10 mΩ, $C_1$은 중앙값 1500 F이고 둘 다 로그 표준편차 0.5인 로그정규분포다(참값과 일부러 다르게 두었다). $R_0$는 즉시 변화로 잘 정해지므로 참값에 고정하고, $(\log R_1, \log C_1)$ 격자 위에서 사후분포를 계산해 $R_1$의 주변분포를 본다.

lr1 = np.linspace(np.log(2e-3), np.log(0.1), 400)     # log R1 격자
lc1 = np.linspace(np.log(200.0), np.log(2e4), 150)    # log C1 격자
LR, LC = [g.ravel() for g in np.meshgrid(lr1, lc1, indexing="ij")]   # (G,) G = 400 × 150
log_prior = -0.5 * ((LR - np.log(0.010)) / 0.5) ** 2 - 0.5 * ((LC - np.log(1500.0)) / 0.5) ** 2  # 사전: 같은 모델 셀들의 값

def r1_summary(log_w):                               # 격자 가중치 → R1 주변분포의 중앙값과 90% 구간 [mΩ]
    w = np.exp(log_w - log_w.max()).reshape(400, 150).sum(axis=1)
    cdf = np.cumsum(w) / w.sum()
    return [np.exp(np.interp(q, cdf, lr1)) * 1e3 for q in (0.05, 0.5, 0.95)]

print("사전분포   R1 = {1:5.1f} mΩ  [{0:5.1f}, {2:5.1f}]".format(*r1_summary(log_prior)))
data = {"펄스 300 s, 60 s 간격": (t_s[::30], i_a[::30], slice(1, None, 2)),   # 30 s step으로 계산, t = 30, 90, ... s만 기록
        "펄스 3 s, 휴지 없음  ": (t_short, i_short, slice(None))}
for name, (tt, ii, keep) in data.items():
    y = (simulate_v(theta_true, tt, ii) + rng.normal(0, 1e-3, tt.shape))[keep]
    pred = simulate_v((theta_true[0], np.exp(LR), np.exp(LC)), tt, ii)[:, keep]   # (G, T), R0는 알려졌다고 고정
    log_post = log_prior - 0.5 * np.sum((pred - y) ** 2, axis=1) / 1e-3**2
    print(f"사후({name}) R1 = {{1:5.1f}} mΩ  [{{0:5.1f}}, {{2:5.1f}}]".format(*r1_summary(log_post)))

사전분포   R1 =  10.0 mΩ  [  4.4,  22.7]
사후(펄스 300 s, 60 s 간격) R1 =  14.8 mΩ  [ 14.6,  15.1]
사후(펄스 3 s, 휴지 없음  ) R1 =   9.8 mΩ  [  4.3,  22.3]

60 s 간격으로만 기록했어도 $R_0$를 알면 $R_1$은 잘 정해진다. 사후 90% 구간이 14.6–15.1 mΩ로 좁고, 사전 중앙값 10 mΩ에서 데이터 쪽(참값 15 mΩ)으로 옮겨 갔다. 3 s 펄스에서는 사후분포가 사전분포와 거의 같다(중앙값 9.8 vs 10.0 mΩ, 구간 4.3–22.3 vs 4.4–22.7 mΩ). 사후분포가 사전분포와 같다는 것은 데이터가 그 파라미터에 대해 아무것도 말하지 않았다는 뜻이다. 이때 사후 평균을 "추정값"으로 보고하면 사전분포를 데이터로 위장한 셈이 된다. 반대로 이 비교는 식별 가능성을 점검하는 좋은 도구다. 사전분포와 사후분포를 나란히 그리는 것만으로 어떤 파라미터가 데이터로 정해졌는지 보인다(24.7장). 파라미터가 많은 전기화학 모델에서는 격자 대신 MCMC(13.5장)를 쓴다.

대리모델

왜 필요한가. 등가회로 모델은 빠르지만 내부 상태(전극의 리튬 농도 분포, 전해질 농도)를 표현하지 못한다. 그래서 리튬 석출 위험이나 고율 방전 용량처럼 내부 상태가 정하는 양을 예측하려면 30.2장의 전기화학 모델(Doyle–Fuller–Newman 계열의 P2D 모델이나 단일 입자 모델)이 필요하다. 이런 모델은 편미분방정식과 대수방정식이 섞인 DAE라서 한 번 푸는 데 수 ms에서 수 초가 걸린다. 파라미터 추정, 설계 최적화, 불확실성 전파처럼 모델을 수천에서 수백만 번 부르는 작업에는 너무 느리다.

정의. 대리모델(surrogate model)은 고충실도 시뮬레이터의 입력–출력 관계를 시뮬레이션 실행 결과로 학습해 빠르게 근사하는 모델이다. 배터리에서 입력은 물성(확산계수, 반응 속도 상수, 입자 크기)과 운전 조건(C-rate, 온도, 초기 SOC)이고, 출력은 스칼라(방전 용량, 최고 온도)이거나 시계열(전압 곡선)이다. GP, 신경망, 축소 모델(POD)의 일반론과 검증 원칙은 24.5장에서 다뤘다. 여기서는 배터리에서 특히 주의할 두 가지를 정리한다.

• 운전 조건 일반화. 대리모델은 학습한 조건의 범위 안에서만 믿을 수 있다. 배터리는 저온·고율처럼 위험하고 흥미로운 조건이 대개 학습 범위의 가장자리에 있다. 그래서 대리모델을 학습 범위 밖에서 따로 평가해야 한다. 0.5–2 C로 학습한 모델이 3 C에서 맞는지는 학습 범위 안의 시험 오차로 알 수 없다.
• 물성 일반화. 물성을 입력으로 넣은 대리모델은 "다른 셀"에 대한 예측이 된다. 이 경우 시험 시뮬레이션은 학습에 쓴 물성 조합과 겹치지 않게 따로 뽑고, 물성 공간의 어느 부분에서 오차가 큰지를 함께 보고한다.

연구에서의 의미. 대리모델이 고충실도 모델을 대신하면, 그 위에서 얻은 결론(최적 설계, 파라미터 추정값, 신뢰 구간)에는 대리모델의 오차가 섞인다. 대리모델로 찾은 최적점이나 추정값은 마지막에 고충실도 모델로 다시 계산해 확인하는 것이 기본이다.

하이브리드 모델

왜 필요한가. 물리 모델은 구조가 있어 학습 범위 밖으로 외삽할 수 있지만, 모든 현상을 담지 못해 편향이 있다. 데이터 기반 모델은 편향을 줄이지만 학습 범위 밖에서 근거 없이 행동한다. 하이브리드 모델(hybrid model)은 물리 모델이 아는 부분은 물리 모델에 맡기고, 모르는 부분만 데이터로 배운다. 결합 방식(잔차 보정, 파라미터 학습, 직렬 결합)의 일반론은 24.5장에서 다뤘다.

물리 구조가 외삽을 정한다. 가장 단순한 하이브리드는 물리 법칙의 함수 꼴을 정하고 계수만 데이터로 맞추는 것이다. 셀 저항의 온도 의존성은 흔히 Arrhenius 꼴 $R_0(T) = R_{\text{ref}}\exp\!\big[\tfrac{E_a}{R_g}\big(\tfrac{1}{T} - \tfrac{1}{T_{\text{ref}}}\big)\big]$로 모델링한다. $E_a$는 활성화 에너지, $R_g$는 기체 상수, $T$는 절대온도다. 15–35 °C에서 펄스 시험으로 추정한 $R_0$ 다섯 개(추정 오차 2%)로 같은 자유도의 두 모델을 맞추고, 학습 범위 밖의 −10 °C, 0 °C, 45 °C에서 비교한다. 참 $R_0(T)$를 Arrhenius 꼴로 만든 가상 셀이므로, 이 실험은 "구조가 맞을 때"의 이점을 보인다.

T_meas_c = np.array([15.0, 20.0, 25.0, 30.0, 35.0])  # 측정 온도 (°C)
def r0_true(T_c):                                    # 가상 셀의 참 R0(T): Arrhenius, Ea = 30 kJ/mol (설명용)
    return 0.020 * np.exp(30e3 / 8.314 * (1 / (T_c + 273.15) - 1 / 298.15))
r0_meas = r0_true(T_meas_c) * (1 + 0.02 * np.random.default_rng(1).normal(size=5))   # 추정 오차 2%

poly = np.polyfit(T_meas_c, r0_meas, 2)              # 데이터 기반: T의 2차 다항식
arr = np.polyfit(1 / (T_meas_c + 273.15), np.log(r0_meas), 1)   # 물리 구조: ln R0 = a + (Ea/R)(1/T)
for T_c in [-10.0, 0.0, 25.0, 45.0]:
    p, a = np.polyval(poly, T_c), np.exp(np.polyval(arr, 1 / (T_c + 273.15)))
    print(f"T = {T_c:5.1f} °C  참값 {r0_true(T_c)*1e3:5.1f} mΩ | 2차 다항식 {p*1e3:5.1f} | Arrhenius {a*1e3:5.1f}")
print(f"Arrhenius 적합의 활성화 에너지 = {arr[0] * 8.314 / 1e3:.1f} kJ/mol (참값 30)")

T = -10.0 °C  참값 100.0 mΩ | 2차 다항식  76.6 | Arrhenius 102.2
T =   0.0 °C  참값  60.5 mΩ | 2차 다항식  55.0 | Arrhenius  61.5
T =  25.0 °C  참값  20.0 mΩ | 2차 다항식  20.0 | Arrhenius  20.1
T =  45.0 °C  참값   9.3 mΩ | 2차 다항식  11.5 | Arrhenius   9.3
Arrhenius 적합의 활성화 에너지 = 30.3 kJ/mol (참값 30)

학습 범위 안인 25 °C에서는 두 모델이 같다. 범위 밖에서는 2차 다항식이 −10 °C에서 참값 100 mΩ를 76.6 mΩ로 23% 낮게 예측한다. 저온에서 저항을 낮게 보면 출력 한계를 과대평가해 전압이 하한 아래로 떨어질 수 있으므로, 이것은 안전 방향의 오류가 아니다. Arrhenius 적합은 −10 °C에서도 2% 오차이고, 적합한 활성화 에너지 30.3 kJ/mol 자체도 물리적으로 해석할 수 있는 양이다.

잔차 보정. 물리 모델이 빠뜨린 현상(LFP 셀의 OCV 히스테리시스, 온도 의존 확산, 노화에 따른 저항 증가)은 물리 모델 출력과 측정의 차이를 데이터 기반 모델로 배워 보정한다. 실습에서 확산을 무시한 저충실도 모델에 GP 잔차 보정을 붙인 하이브리드가 데이터만 쓴 GP보다 저온 외삽 오차를 크게 줄이는 것을 확인한다.

열화 모델링. 하이브리드 구조는 열화에서도 자주 쓰인다. 용량 감소를 SEI(solid electrolyte interphase) 성장처럼 물리적 근거가 있는 함수 꼴(예: 시간의 제곱근에 비례하는 성장)로 두고, 그 계수가 온도·SOC·전류에 어떻게 의존하는지를 데이터로 배우는 방식이다. 이렇게 하면 열화 모델의 파라미터가 사이클마다 다시 추정되는 등가회로 파라미터($R_0$의 증가)와 연결되어, 30.4장의 수명 예측에 물리적 해석을 줄 수 있다. 다만 짧은 기간의 열화 데이터로 함수 꼴을 고르면 여러 꼴이 같은 데이터를 똑같이 잘 설명하는 경우가 많다. 서로 다른 함수 꼴의 예측이 장기에서 얼마나 갈라지는지를 함께 보고한다.

평가

대리모델과 하이브리드 모델은 세 축으로 평가한다. 한 축만 보고하면 모델의 쓸모를 판단할 수 없다.

• 정확도. 학습 범위 안의 시험 오차와 학습 범위 밖의 오차를 나누어 보고한다. 배터리에서 범위 밖은 대개 저온, 고율, 노화가 진행한 셀이다. 평균 오차(RMSE)만이 아니라 최대 오차도 쓴다. BMS는 최악의 경우를 견뎌야 하기 때문이다.
• 물리적 일관성. 예측이 물리 법칙과 상식을 어기지 않는지 따로 검사한다. 방전 용량은 0과 정격 용량 사이여야 하고, 같은 온도에서 C-rate가 커지면 줄어야 하며, 저항은 양수이고 온도가 내려가면 커져야 한다. 이런 검사는 학습 데이터가 없는 조건까지 촘촘한 격자로 한다. 정확도가 좋은 모델도 일관성을 어길 수 있고, 일관성 위반은 그 모델을 최적화나 제어에 넣었을 때 최적화기가 파고드는 구멍이 된다.
• 속도. 평가 1회의 시간을 고충실도 모델과 같은 조건(같은 기계, batch 크기 명시)에서 잰다. 학습에 쓴 고충실도 실행의 비용도 함께 쓴다. 대리모델을 만드는 데 고충실도 실행 1만 번이 들었는데 실제 사용에서 1000번만 부른다면 대리모델은 이득이 아니다.

세 축 사이에는 대개 trade-off가 있다. 실습의 결과처럼 하이브리드는 정확하지만 저충실도 물리 모델을 함께 계산하므로 순수 데이터 모델보다 느리다. 어느 쪽을 고를지는 쓰임새(실시간 BMS, 오프라인 설계 최적화)가 정한다.

참고 문헌

• Plett (2015). Battery Management Systems, Volume II: Equivalent-Circuit Methods. Artech House. — 등가회로 모델과 그 파라미터 추정, BMS에서의 쓰임
• Raue, Kreutz, Maiwald, Bachmann, Schilling, Klingmüller, Timmer (2009). Structural and practical identifiability analysis of partially observed dynamical models by exploiting the profile likelihood. Bioinformatics, 25(15), 1923–1929. — 구조적·실용적 식별 가능성과 프로파일 우도