3S LiPo Battery Pack Thevenin ECM (12-state)

def _mission_current(t: float) -> float:
    """Pack current draw (A) for a typical survey mission.

    0–60 s:    hover (15 A)
    60–300 s:  cruise (8 A)
    300–360 s: hover (15 A)
    360–420 s: descent (5 A)
    420+ s:    idle (0.5 A)
    """
    if t < 60.0:
        return 15.0
    elif t < 300.0:
        return 8.0
    elif t < 360.0:
        return 15.0
    elif t < 420.0:
        return 5.0
    else:
        return 0.5

def _r0_effective(soc: float, temp: float, variability: float) -> float:
    """R0 as a function of SOC, temperature, and cell variability."""
    soc_factor = 1.0 + _R0_SOC_FACTOR * (1.0 - np.clip(soc, 0.0, 1.0))
    temp_factor = 1.0 + _ALPHA_R * (temp - 25.0)
    return _R0_BASE * soc_factor * temp_factor * variability

def rhs_battery_3s(t, y):
    """3S LiPo pack: 12-state 2nd-order Thevenin ECM with thermal coupling."""
    dy = np.zeros(12)
    i_pack = _mission_current(t)

    for cell in range(3):
        base = cell * 4
        soc = y[base]
        v_rc1 = y[base + 1]
        v_rc2 = y[base + 2]
        t_cell = y[base + 3]

        r0 = _r0_effective(soc, t_cell, _R0_VARIABILITY[cell])

        # SOC dynamics: dSOC/dt = -I / Q
        dy[base] = -i_pack / _Q_COULOMBS

        # RC pair dynamics
        dy[base + 1] = (i_pack / _C1) - (v_rc1 / _TAU1)
        dy[base + 2] = (i_pack / _C2) - (v_rc2 / _TAU2)

        # Thermal: lumped node with ohmic + polarisation heat generation
        p_ohmic = i_pack**2 * r0
        p_rc1 = v_rc1**2 / _R1
        p_rc2 = v_rc2**2 / _R2
        p_total = p_ohmic + p_rc1 + p_rc2
        q_conv = (t_cell - _T_AMB) / _R_TH_CONV
        dy[base + 3] = (p_total - q_conv) / _C_TH

    return dy