19 KiB
% ========================= % 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)<w.\end{cases}
\label{eq:wdot_sec6}
\end{align}
We discretize in time with RK4:
\begin{equation}
\mathbf{x}_{n+1}=\mathbf{x}_n+\frac{\Delta t}{6}\left(\mathbf{k}_1+2\mathbf{k}_2+2\mathbf{k}_3+\mathbf{k}_4\right),
\label{eq:rk4_update_sec6}
\end{equation}
where each stage \mathbf{k}_j=f(\mathbf{x},\mathbf{u},I(\mathbf{x},\mathbf{u})) is computed at
the corresponding intermediate state and time, and \emph{each} stage uses the nested algebraic
evaluation \eqref{eq:Delta_sec6}--\eqref{eq:I_limit_sec6}. Inputs \mathbf{u}(t) at intermediate
times are obtained by interpolation (piecewise-constant or piecewise-linear) or by direct
evaluation if \mathbf{u}(t) is generated procedurally.
\subsection{Projection, physical constraints, and robustness} \label{subsec:projection_robustness}
To prevent numerical drift outside admissible ranges, we apply a mild projection after each
successful time step:
\begin{equation}
z\leftarrow \min(1,\max(0,z)),\quad
S\leftarrow \min(1,\max(0,S)),\quad
w\leftarrow \min(1,\max(0,w)).
\label{eq:projection_sec6}
\end{equation}
This projection is used only to suppress floating-point accumulation and does not change the
continuous model definition. Additionally, the low-SOC protection is enforced solely in the OCV
calculation via z_{\mathrm{eff}}=\max\{z,z_{\min}\} to avoid the (1/z) singularity while preserving
the runtime termination criterion based on z(t) and V_{\mathrm{term}}(t).
Robustness safeguards include: (i) enforcing Q_{\mathrm{eff}}(T_b,S)\ge 0 via the [\cdot]_+ operator;
(ii) recording infeasibility when \Delta<0; and (iii) saturating current by
\eqref{eq:I_limit_sec6}, which prevents unrealistically large currents under low-voltage conditions.
\subsection{Step-size selection, stability, and convergence control} \label{subsec:stepsize_convergence}
The fastest electrical time scale is the polarization time constant
\begin{equation}
\tau_p=R_1C_1.
\label{eq:taup_sec6}
\end{equation}
To resolve the polarization dynamics and avoid stage-level oscillations in an explicit method,
we impose the time-step bound
\begin{equation}
\Delta t \le 0.05,\tau_p.
\label{eq:dt_bound_sec6}
\end{equation}
When tail dynamics are active, one may further restrict \Delta t \le 0.05\,\tau_\uparrow.
We enforce convergence through step-halving. Over a candidate step from t_n to t_{n+1}=t_n+\Delta t,
we compute two solutions: one using a single RK4 step of size \Delta t and another using two RK4
steps of size \Delta t/2. The step is accepted if
\begin{equation}
\left|z_{\Delta t}-z_{\Delta t/2}\right|\infty < 10^{-4},
\label{eq:step_halving_soc_sec6}
\end{equation}
and, for runtime outputs, the inferred TTE changes by less than 1\% under step-halving in the
same scenario:
\begin{equation}
\frac{\big|\mathrm{TTE}{\Delta t}-\mathrm{TTE}{\Delta t/2}\big|}{\mathrm{TTE}{\Delta t/2}} < 1%.
\label{eq:step_halving_tte_sec6}
\end{equation}
If the criterion fails, we set \Delta t\leftarrow \Delta t/2 and recompute the step.
\subsection{Event detection and TTE interpolation} \label{subsec:event_detection}
The runtime termination time is defined as
\begin{equation}
\mathrm{TTE}=\inf\left{t>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}