% ========================= % Section 6: Numerical Solution & Identification % ========================= \section{Numerical Solution and Parameter Identification} \label{sec:numerics_id} This section describes a reproducible computational implementation of the coupled ODE--algebraic closure induced by the constant-power-load (CPL) constraint, and a mechanism-driven parameter identification workflow. We employ an explicit fourth-order Runge--Kutta integrator (RK4) with a nested algebraic evaluation of the discharge current at each substage, adaptive step-halving for convergence control, and event detection for the time-to-end (TTE). Parameter estimation is performed by targeted sub-experiments and log-based regressions that preserve physical interpretability. \subsection{RK4 with substage nested algebraic evaluation of $I$} \label{subsec:rk4_nestedI} The state vector is $\mathbf{x}(t)=[z(t),v_p(t),T_b(t),S(t),w(t)]^\top$ and the input vector is $\mathbf{u}(t)=[L(t),C(t),N(t),\Psi(t),T_a(t)]^\top$. For any $(\mathbf{x},\mathbf{u})$, we compute the total power demand \begin{equation} P_{\mathrm{tot}}(t)=P_{\mathrm{bg}}+P_{\mathrm{scr}}(L(t))+P_{\mathrm{cpu}}(C(t))+P_{\mathrm{net}}(N(t),\Psi(t),w(t)). \label{eq:Ptot_def_sec6} \end{equation} The terminal voltage satisfies the ECM relation \begin{equation} V_{\mathrm{term}}(t)=V_{\mathrm{oc}}(z_{\mathrm{eff}}(t))-v_p(t)-I(t)\,R_0(T_b(t),S(t)), \label{eq:Vterm_sec6} \end{equation} with the low-SOC protection $z_{\mathrm{eff}}(t)=\max\{z(t),z_{\min}\}$. Under the CPL assumption, \begin{equation} P_{\mathrm{tot}}(t)=V_{\mathrm{term}}(t)\,I(t) =\big(V_{\mathrm{oc}}(z_{\mathrm{eff}})-v_p-I R_0\big)\,I, \label{eq:CPL_sec6} \end{equation} which yields the discriminant \begin{equation} \Delta(t)=\big(V_{\mathrm{oc}}(z_{\mathrm{eff}})-v_p\big)^2-4R_0(T_b,S)\,P_{\mathrm{tot}}(t). \label{eq:Delta_sec6} \end{equation} If $\Delta(t)\ge 0$, the physically consistent branch of the quadratic solution is \begin{equation} I_{\mathrm{CPL}}(t)=\frac{V_{\mathrm{oc}}(z_{\mathrm{eff}}(t))-v_p(t)-\sqrt{\Delta(t)}}{2R_0(T_b(t),S(t))}. \label{eq:Icpl_sec6} \end{equation} To reflect device-side protection (PMIC current limiting / OS throttling), we apply a temperature-dependent saturation \begin{equation} I(t)=\min\big(I_{\mathrm{CPL}}(t),\,I_{\max}(T_b(t))\big), \qquad I_{\max}(T_b)=I_{\max,0}\big[1-\rho_T\,(T_b-T_{\mathrm{ref}})\big]_+. \label{eq:I_limit_sec6} \end{equation} When $\Delta(t)<0$, CPL delivery is infeasible. In such cases we record a ``collapse-risk'' event (see Section~\ref{subsec:event_detection}) and place the system in a strong degradation regime by taking $I(t)=I_{\max}(T_b(t))$ (equivalently, one may cap $P_{\mathrm{tot}}$), while the runtime termination (TTE) is still defined by voltage/SOC cutoffs. Given the current mapping $I=I(\mathbf{x},\mathbf{u})$, the ODE right-hand side is evaluated using the established dynamics: \begin{align} \dot z &= -\frac{I}{3600\,Q_{\mathrm{eff}}(T_b,S)}, \label{eq:zdot_sec6}\\ \dot v_p &= \frac{I}{C_1}-\frac{v_p}{R_1C_1}, \label{eq:vpdot_sec6}\\ \dot T_b &= \frac{1}{C_{\mathrm{th}}}\Big(I^2R_0(T_b,S)+\frac{v_p^2}{R_1}-hA\,(T_b-T_a)\Big), \label{eq:Tbdot_sec6}\\ \dot S &= -\lambda_{\mathrm{sei}}|I|^{m}\exp\!\left(-\frac{E_{\mathrm{sei}}}{R_gT_b}\right), \label{eq:Sdot_sec6}\\ \dot w &= \frac{\sigma(N)-w}{\tau(N)}, \quad \sigma(N)=\min(1,N), \quad \tau(N)=\begin{cases}\tau_\uparrow,&\sigma(N)\ge w,\\ \tau_\downarrow,&\sigma(N)0:\ V_{\mathrm{term}}(t)\le V_{\mathrm{cut}}\ \text{or}\ z(t)\le 0\right\}. \label{eq:TTE_def_sec6} \end{equation} During integration, we monitor $V_{\mathrm{term}}(t)-V_{\mathrm{cut}}$ and $z(t)$ for sign changes. If $V_{\mathrm{term}}(t_n)>V_{\mathrm{cut}}$ but $V_{\mathrm{term}}(t_{n+1})\le V_{\mathrm{cut}}$, we approximate the crossing time by linear interpolation: \begin{equation} t_\star \approx t_n + \Delta t\,\frac{V_{\mathrm{term}}(t_n)-V_{\mathrm{cut}}}{V_{\mathrm{term}}(t_n)-V_{\mathrm{term}}(t_{n+1})}. \label{eq:interp_voltage_event_sec6} \end{equation} Similarly, if $z(t_n)>0$ and $z(t_{n+1})\le 0$, \begin{equation} t_\star \approx t_n + \Delta t\,\frac{z(t_n)}{z(t_n)-z(t_{n+1})}. \label{eq:interp_soc_event_sec6} \end{equation} If both events occur within the same step, we take the earlier of the two interpolated times as TTE. In addition, we optionally record a CPL infeasibility (voltage-collapse risk) time \begin{equation} t_\Delta=\inf\{t>0:\ \Delta(t)\le 0\}, \label{eq:tDelta_def_sec6} \end{equation} which is useful for diagnosing ``sudden shutdown'' risk even when current limiting postpones the actual cutoff event. \subsection{Algorithm 3: Simulation procedure} \label{subsec:algorithm3} \begin{algorithm}[t] \caption{RK4 simulation with nested CPL current evaluation and event handling} \label{alg:rk4_cpl} \begin{algorithmic}[1] \REQUIRE Initial state $\mathbf{x}(0)=[z_0,0,T_a(0),S_0,0]^\top$, input trajectory $\mathbf{u}(t)$, parameters $\Theta$, cutoff $V_{\mathrm{cut}}$, step bound $\Delta t_{\max}$. \ENSURE Trajectories $\mathbf{x}(t)$, $V_{\mathrm{term}}(t)$, and $\mathrm{TTE}$ (and optionally $t_\Delta$). \STATE Set $t\leftarrow 0$, $\mathbf{x}\leftarrow \mathbf{x}(0)$, choose $\Delta t\le \Delta t_{\max}$. \STATE Initialize flags: $\texttt{risk\_recorded}\leftarrow \texttt{false}$. \WHILE{$V_{\mathrm{term}}(t)>V_{\mathrm{cut}}$ \AND $z(t)>0$} \STATE Evaluate $\mathbf{u}(t)$ and (if needed) $\mathbf{u}(t+\Delta t/2)$, $\mathbf{u}(t+\Delta t)$. \STATE Compute $z_{\mathrm{eff}}=\max\{z,z_{\min}\}$, then $V_{\mathrm{oc}}(z_{\mathrm{eff}})$, $R_0(T_b,S)$, $Q_{\mathrm{eff}}(T_b,S)$, and $P_{\mathrm{tot}}$ via \eqref{eq:Ptot_def_sec6}. \STATE Compute $\Delta$ via \eqref{eq:Delta_sec6}. \IF{$\Delta<0$} \IF{\NOT $\texttt{risk\_recorded}$} \STATE Record $t_\Delta\leftarrow t$; $\texttt{risk\_recorded}\leftarrow \texttt{true}$. \ENDIF \STATE Set $I\leftarrow I_{\max}(T_b)$ \COMMENT{strong degradation / protection} \ELSE \STATE Compute $I_{\mathrm{CPL}}$ via \eqref{eq:Icpl_sec6} and apply saturation \eqref{eq:I_limit_sec6}. \ENDIF \STATE Perform one RK4 step \eqref{eq:rk4_update_sec6} with nested current evaluation at each substage. \STATE Apply projection \eqref{eq:projection_sec6}. \STATE Step-halving check: compare $\Delta t$ vs.\ two half-steps; if \eqref{eq:step_halving_soc_sec6} fails, set $\Delta t\leftarrow \Delta t/2$ and recompute this step. \STATE Update $V_{\mathrm{term}}$ via \eqref{eq:Vterm_sec6}; test event conditions. \IF{event detected within this step} \STATE Interpolate event time using \eqref{eq:interp_voltage_event_sec6} or \eqref{eq:interp_soc_event_sec6}. \STATE Set $\mathrm{TTE}\leftarrow t_\star$ and \textbf{break}. \ENDIF \STATE Update $t\leftarrow t+\Delta t$. \ENDWHILE \RETURN $\mathrm{TTE}$, trajectories, and optionally $t_\Delta$. \end{algorithmic} \end{algorithm} \subsection{Overall strategy for parameter identification} \label{subsec:id_strategy} Parameters are grouped to enable targeted identification with minimal confounding: \begin{itemize} \item \textbf{Open-circuit voltage (OCV)} parameters $(E_0,K,A,B)$ are identified from an OCV--SOC curve. \item \textbf{ECM electrical parameters} $(R_0,R_1,C_1)$ are identified from current pulse tests by separating the instantaneous ohmic drop from the relaxation dynamics. \item \textbf{Thermal parameters} $(C_{\mathrm{th}},hA)$ are identified from heating and cooling transients. \item \textbf{Aging parameters} $(\lambda_{\mathrm{sei}},m,E_{\mathrm{sei}})$ are identified from capacity fade data under controlled current/temperature conditions. \item \textbf{Device power-mapping parameters} (screen/CPU/network) are identified from controlled workload logs by isolating each subsystem and fitting the prescribed mechanistic forms. \item \textbf{Tail parameters} $(k_{\mathrm{tail}},\tau_\uparrow,\tau_\downarrow)$ are identified from network-burst experiments by fitting the post-burst decay shape. \end{itemize} This staged approach preserves physical interpretability and avoids black-box regressions. \subsection{OCV fitting: $(E_0,K,A,B)$} \label{subsec:ocv_fit} Given OCV--SOC samples $\{(z_i,V_i)\}_{i=1}^M$ collected under quasi-equilibrium conditions, we estimate $(E_0,K,A,B)$ by least squares using the protected SOC $z_{i,\mathrm{eff}}=\max\{z_i,z_{\min}\}$: \begin{equation} \min_{E_0,K,A,B}\ \sum_{i=1}^M\left[ V_i-\left(E_0-K\left(\frac{1}{z_{i,\mathrm{eff}}}-1\right)+A e^{-B(1-z_{i,\mathrm{eff}})}\right) \right]^2. \label{eq:ocv_ls_sec6} \end{equation} The resulting OCV model is then used in the time-domain simulations through \eqref{eq:Vterm_sec6}. \subsection{Pulse-based identification: $(R_0,R_1,C_1)$} \label{subsec:pulse_id} \paragraph{Ohmic resistance $R_0$.} At fixed SOC and temperature, apply a current step of magnitude $\Delta I$ and measure the instantaneous voltage drop $\Delta V(0^+)$, yielding \begin{equation} R_0 \approx \frac{\Delta V(0^+)}{\Delta I}. \label{eq:R0_pulse_sec6} \end{equation} \paragraph{Polarization branch $(R_1,C_1)$.} After removing the ohmic drop, the remaining relaxation is approximately first-order with time constant $\tau_p=R_1C_1$. Denote the relaxation component by $V_{\mathrm{rel}}(t)=V_{\mathrm{term}}(t)-\big(V_{\mathrm{oc}}-\Delta I\,R_0\big)$. Then \begin{equation} V_{\mathrm{rel}}(t)\approx -\Delta I\,R_1\,e^{-t/\tau_p}, \label{eq:relax_exp_sec6} \end{equation} so that a linear fit of $\ln|V_{\mathrm{rel}}(t)|$ versus $t$ yields $\tau_p$ and $R_1$, and hence \begin{equation} C_1=\frac{\tau_p}{R_1}. \label{eq:C1_from_tau_sec6} \end{equation} \subsection{Temperature/aging coupling: $(R_{\mathrm{ref}},E_a,\eta_R,Q_{\mathrm{nom}},\alpha_Q)$} \label{subsec:temp_aging_coupling} \paragraph{Arrhenius temperature dependence for $R_0$.} Measure $R_0$ at multiple temperatures $T_b^{(j)}$ (e.g., by \eqref{eq:R0_pulse_sec6}) and fit \begin{equation} \ln R_0^{(j)}=\ln R_{\mathrm{ref}}+\frac{E_a}{R_g}\left(\frac{1}{T_b^{(j)}}-\frac{1}{T_{\mathrm{ref}}}\right), \label{eq:arrhenius_fit_sec6} \end{equation} to obtain $R_{\mathrm{ref}}$ and $E_a$. \paragraph{SOH correction for resistance.} Using measurements across different SOH levels $S$, fit \begin{equation} \frac{R_0(T_b,S)}{R_0(T_b,1)}\approx 1+\eta_R(1-S) \label{eq:etaR_fit_sec6} \end{equation} to obtain $\eta_R$. \paragraph{Effective capacity parameters.} From capacity tests across temperatures, estimate $Q_{\mathrm{nom}}$ and $\alpha_Q$ using \begin{equation} Q_{\mathrm{eff}}(T_b,S)=Q_{\mathrm{nom}}\,S\left[1-\alpha_Q(T_{\mathrm{ref}}-T_b)\right]_+. \label{eq:Qeff_fit_sec6} \end{equation} \subsection{Power mapping identification: $(k_L,\gamma,k_C,\eta,k_N,\kappa,\ldots)$} \label{subsec:power_mapping_id} \paragraph{Screen mapping.} Under controlled conditions with minimal CPU/network activity, vary brightness $L$ and measure total power. After subtracting background and CPU baseline, fit \begin{equation} P_{\mathrm{scr}}(L)=P_{\mathrm{scr},0}+k_L L^\gamma,\qquad \gamma>1. \label{eq:screen_fit_sec6} \end{equation} \paragraph{CPU mapping.} With fixed brightness and network conditions, apply controlled workloads to vary CPU load $C$ and fit \begin{equation} P_{\mathrm{cpu}}(C)=P_{\mathrm{cpu},0}+k_C C^\eta,\qquad \eta>1. \label{eq:cpu_fit_sec6} \end{equation} \paragraph{Network mapping and signal-quality penalty.} At fixed throughput proxy $N=N_0$, vary signal quality $\Psi$ and fit \begin{equation} P_{\mathrm{net}}(N_0,\Psi,w)\approx P_{\mathrm{net},0}+k_N\frac{N_0}{(\Psi+\varepsilon)^\kappa}+k_{\mathrm{tail}}w. \label{eq:net_fit_sec6} \end{equation} For steady experiments where $w$ is constant or negligible, define $\Delta P_{\mathrm{net}}(\Psi)=P_{\mathrm{net}}-P_{\mathrm{net},0}$ and fit in log space: \begin{equation} \ln \Delta P_{\mathrm{net}}(\Psi)\approx \ln(k_N N_0)-\kappa\ln(\Psi+\varepsilon), \label{eq:kappa_fit_sec6} \end{equation} yielding $\kappa$ (slope) and $k_N$ (intercept). \subsection{Tail parameter identification: $(k_{\mathrm{tail}},\tau_\uparrow,\tau_\downarrow)$} \label{subsec:tail_id} Conduct a network-burst experiment: drive $N(t)$ high for a short period and then reduce it rapidly. After the burst, $N\approx 0$ and the tail state decays approximately exponentially with time constant $\tau_\downarrow$: \begin{equation} w(t)\approx w(t_0)\,e^{-(t-t_0)/\tau_\downarrow},\qquad P_{\mathrm{tail}}(t)=k_{\mathrm{tail}}w(t). \label{eq:tail_decay_sec6} \end{equation} A linear fit of $\ln P_{\mathrm{tail}}(t)$ versus $t$ yields $\tau_\downarrow$, and the amplitude identifies $k_{\mathrm{tail}}$. The rise time $\tau_\uparrow$ is obtained by fitting the initial ramp-up segment during burst onset, consistent with $\tau_\uparrow\ll\tau_\downarrow$. \subsection{Thermal parameter identification: $(C_{\mathrm{th}},hA)$} \label{subsec:thermal_id} Using a heating--cooling experiment, identify the lumped thermal time constant. During the cooling phase where $I\approx 0$ and $v_p\approx 0$, \eqref{eq:Tbdot_sec6} reduces to \begin{equation} \dot T_b \approx -\frac{hA}{C_{\mathrm{th}}}(T_b-T_a), \label{eq:cooling_sec6} \end{equation} so that \begin{equation} T_b(t)-T_a \approx (T_b(t_0)-T_a)\,e^{-(hA/C_{\mathrm{th}})(t-t_0)}. \label{eq:cooling_exp_sec6} \end{equation} Fitting the exponential decay yields $hA/C_{\mathrm{th}}$. Then, using the heating phase with known heat generation $\dot Q \approx I^2R_0+v_p^2/R_1$, one can separate $C_{\mathrm{th}}$ and $hA$. \subsection{Aging parameter identification: $(\lambda_{\mathrm{sei}},m,E_{\mathrm{sei}})$} \label{subsec:aging_id} From controlled aging data providing $S(t)$ under known $(I,T_b)$ conditions, the SEI-driven degradation model \eqref{eq:Sdot_sec6} can be identified by log-linear regression. Approximating $\dot S$ via finite differences, \begin{equation} \ln(-\dot S)\approx \ln \lambda_{\mathrm{sei}} + m\ln|I| - \frac{E_{\mathrm{sei}}}{R_g}\frac{1}{T_b}. \label{eq:aging_loglin_sec6} \end{equation} A multi-condition fit across varying currents and temperatures yields $(\lambda_{\mathrm{sei}},m,E_{\mathrm{sei}})$. This procedure preserves the mechanistic form and avoids black-box regression. \subsection{Parameter table: nominal values, ranges, and sources} \label{subsec:param_table} Table~\ref{tab:param_summary} summarizes the parameters used in simulations, including nominal values and uncertainty ranges for sensitivity analysis. Nominal values are obtained via the identification procedures above or from manufacturer specifications / literature when direct measurements are unavailable. Ranges should be selected to reflect measurement uncertainty and device-to-device variability (e.g., $\pm 10\%$--$\pm 20\%$ for power-map gains, and temperature-dependent parameters constrained by Arrhenius fits). \begin{table}[t] \centering \caption{Parameter summary (to be finalized): nominal values, uncertainty ranges, and sources.} \label{tab:param_summary} \begin{tabular}{llll} \hline Category & Parameter & Nominal / Range & Source / Method \\ \hline OCV & $E_0,K,A,B$ & (fill) / (fill) & OCV--SOC LS fit \eqref{eq:ocv_ls_sec6} \\ ECM & $R_{\mathrm{ref}},E_a$ & (fill) / (fill) & Arrhenius fit \eqref{eq:arrhenius_fit_sec6} \\ ECM & $R_1,C_1$ & (fill) / (fill) & Pulse relaxation \eqref{eq:relax_exp_sec6} \\ SOH coupling & $\eta_R$ & (fill) / (fill) & Resistance vs.\ SOH \eqref{eq:etaR_fit_sec6} \\ Capacity & $Q_{\mathrm{nom}},\alpha_Q$ & (fill) / (fill) & Capacity tests \eqref{eq:Qeff_fit_sec6} \\ Thermal & $C_{\mathrm{th}},hA$ & (fill) / (fill) & Cooling/heating fits \eqref{eq:cooling_exp_sec6} \\ Aging & $\lambda_{\mathrm{sei}},m,E_{\mathrm{sei}}$ & (fill) / (fill) & Log-linear fit \eqref{eq:aging_loglin_sec6} \\ Screen & $P_{\mathrm{scr},0},k_L,\gamma$ & (fill) / (fill) & Screen power fit \eqref{eq:screen_fit_sec6} \\ CPU & $P_{\mathrm{cpu},0},k_C,\eta$ & (fill) / (fill) & CPU power fit \eqref{eq:cpu_fit_sec6} \\ Network & $P_{\mathrm{net},0},k_N,\kappa$ & (fill) / (fill) & Signal penalty \eqref{eq:kappa_fit_sec6} \\ Tail & $k_{\mathrm{tail}},\tau_\uparrow,\tau_\downarrow$ & (fill) / (fill) & Burst/decay \eqref{eq:tail_decay_sec6} \\ Protection & $I_{\max,0},\rho_T,z_{\min}$ & (fill) / (fill) & Device policy / assumption \\ \hline \end{tabular} \end{table}