567提交

This commit is contained in:
ChuXun
2026-01-30 20:40:15 +08:00
parent 13151b06de
commit dbb62485a7
3 changed files with 961 additions and 0 deletions

View File

@@ -0,0 +1,176 @@
\section{Model Formulation}\label{sec:model}
We develop a mechanism-driven continuous-time model for smartphone battery drain that couples
(i) component-level power mapping from user/device inputs,
(ii) an equivalent-circuit battery model (ECM) with polarization memory,
(iii) a constant-power-load (CPL) algebraic closure for the discharge current,
(iv) lumped thermal dynamics, and
(v) slow health degradation (SOH).
All symbols are used consistently throughout.
\subsection{Total Power Decomposition $P_{\rm tot}$ (Screen/CPU/Network)}\label{sec:ptot}
Let the state vector be
\begin{equation}
\mathbf{x}(t)=\big[z(t),\,v_p(t),\,T_b(t),\,S(t),\,w(t)\big]^\top,
\end{equation}
where $z$ is the state-of-charge (SOC), $v_p$ is the polarization voltage, $T_b$ is the battery temperature,
$S$ is the state-of-health (SOH, capacity fraction), and $w$ is a continuous radio-tail state.
The exogenous input vector is
\begin{equation}
\mathbf{u}(t)=\big[L(t),\,C(t),\,N(t),\,\Psi(t),\,T_a(t)\big]^\top,
\end{equation}
where $L$ is screen brightness, $C$ is CPU load, $N$ is network activity intensity,
$\Psi$ is signal quality (larger is better), and $T_a$ is ambient temperature.
We model the instantaneous total power demand as an additive decomposition
\begin{equation}\label{eq:ptot_def}
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)),
\end{equation}
where $P_{\mathrm{bg}}$ is background/baseline power. The component mappings are chosen to be explicit and
mechanism-consistent:
\begin{align}
P_{\mathrm{scr}}(L)&=P_{\mathrm{scr},0}+k_L L^\gamma,\qquad \gamma>1,\label{eq:pscr}\\
P_{\mathrm{cpu}}(C)&=P_{\mathrm{cpu},0}+k_C C^\eta,\qquad \eta>1,\label{eq:pcpu}\\
P_{\mathrm{net}}(N,\Psi,w)&=P_{\mathrm{net},0}+k_N\frac{N}{(\Psi+\varepsilon)^\kappa}+k_{\mathrm{tail}}w,
\qquad \kappa>0,\ \varepsilon>0.\label{eq:pnet}
\end{align}
Here $(\Psi+\varepsilon)^{-\kappa}$ captures the increased radio power required under poor signal quality,
and $k_{\mathrm{tail}}w$ represents residual ``tail'' consumption after network bursts.
\subsection{Continuous Radio-Tail Dynamics $w(t)$}\label{sec:tail}
Instead of a discrete finite-state-machine tail model, we introduce a continuous tail state $w(t)\in[0,1]$:
\begin{equation}\label{eq:w_dyn}
\dot w(t)=\frac{\sigma(N(t))-w(t)}{\tau(N(t))},
\end{equation}
where
\begin{equation}\label{eq:sigma_tau}
\sigma(N)=\min(1,N),\qquad
\tau(N)=
\begin{cases}
\tau_\uparrow, & \sigma(N)\ge w,\\
\tau_\downarrow,& \sigma(N)< w,
\end{cases}
\qquad \tau_\uparrow\ll\tau_\downarrow.
\end{equation}
This formulation yields fast engagement of the tail state during activity increases and slow decay after activity
drops, while maintaining continuity and numerical robustness.
\subsection{ECM Terminal Voltage Equation}\label{sec:ecm}
We adopt a first-order Thevenin ECM with an ohmic resistance and one polarization branch:
\begin{equation}\label{eq:vterm}
V_{\mathrm{term}}(t)=V_{\mathrm{oc}}(z(t)) - v_p(t) - I(t)\,R_0(T_b(t),S(t)),
\end{equation}
where $V_{\mathrm{oc}}(z)$ is the open-circuit voltage (OCV) as a function of SOC, and
$R_0(T_b,S)$ is the temperature- and SOH-dependent ohmic resistance.
\subsection{CPL Closure: Quadratic Current and Discriminant $\Delta$}\label{sec:cpl}
Smartphone loads are well-approximated as constant-power over short time scales.
We therefore impose a CPL constraint:
\begin{equation}\label{eq:cpl}
P_{\mathrm{tot}}(t)=V_{\mathrm{term}}(t)\,I(t)
=\big(V_{\mathrm{oc}}(z)-v_p-I R_0(T_b,S)\big)I.
\end{equation}
This yields a quadratic equation in $I$ with discriminant
\begin{equation}\label{eq:delta}
\Delta(t)=\big(V_{\mathrm{oc}}(z)-v_p\big)^2-4R_0(T_b,S)P_{\mathrm{tot}}(t).
\end{equation}
Feasibility requires $\Delta(t)\ge 0$. When feasible, the physically consistent branch is
\begin{equation}\label{eq:I_cpl}
I_{\mathrm{CPL}}(t)=\frac{V_{\mathrm{oc}}(z)-v_p-\sqrt{\Delta(t)}}{2R_0(T_b,S)}.
\end{equation}
If $\Delta(t)<0$, the demanded power is not deliverable under the CPL assumption, indicating voltage-collapse risk.
\subsection{Coupled ODEs: SOC--Polarization--Thermal--SOH--Tail}\label{sec:odes}
Given $I(t)$, the coupled state dynamics are
\begin{align}
\dot z(t)&=-\frac{I(t)}{3600\,Q_{\mathrm{eff}}(T_b(t),S(t))},\label{eq:dz}\\
\dot v_p(t)&=\frac{I(t)}{C_1}-\frac{v_p(t)}{R_1C_1},\label{eq:dvp}\\
\dot T_b(t)&=\frac{1}{C_{\mathrm{th}}}\Big(I(t)^2R_0(T_b,S)+\frac{v_p(t)^2}{R_1}-hA\big(T_b(t)-T_a(t)\big)\Big),\label{eq:dTb}\\
\dot S(t)&=-\lambda_{\mathrm{sei}}|I(t)|^{m}\exp\!\left(-\frac{E_{\mathrm{sei}}}{R_gT_b(t)}\right),\qquad 0\le m\le 1,\label{eq:dS}\\
\dot w(t)&=\frac{\sigma(N(t))-w(t)}{\tau(N(t))}.\label{eq:dw}
\end{align}
Equation \eqref{eq:dTb} uses a nonnegative polarization dissipation term $v_p^2/R_1$ for energetic consistency.
\subsection{Constitutive Relations: OCV, $R_0(T_b,S)$, and $Q_{\rm eff}(T_b,S)$}\label{sec:constitutive}
\paragraph{OCV (modified Shepherd).}
We use a modified Shepherd form:
\begin{equation}\label{eq:voc_raw}
V_{\mathrm{oc}}(z)=E_0-K\Big(\frac{1}{z}-1\Big)+A e^{-B(1-z)}.
\end{equation}
\paragraph{Ohmic resistance with Arrhenius temperature dependence and SOH correction.}
\begin{equation}\label{eq:R0}
R_0(T_b,S)=R_{\mathrm{ref}}\exp\!\Big[\frac{E_a}{R_g}\Big(\frac{1}{T_b}-\frac{1}{T_{\mathrm{ref}}}\Big)\Big]\big(1+\eta_R(1-S)\big).
\end{equation}
\paragraph{Effective capacity.}
\begin{equation}\label{eq:Qeff}
Q_{\mathrm{eff}}(T_b,S)=Q_{\mathrm{nom}}\,S\Big[1-\alpha_Q(T_{\mathrm{ref}}-T_b)\Big]_+,
\qquad [x]_+=\max(x,0).
\end{equation}
\subsection{Incorporating Three Lightweight Refinements}\label{sec:refinements}
To improve robustness while preserving the mechanistic structure, we incorporate three ``micro-refinements.''
\paragraph{(i) Low-SOC singularity protection in $V_{\mathrm{oc}}$.}
The term $1/z$ in \eqref{eq:voc_raw} is numerically singular as $z\to 0$.
We introduce an effective SOC
\begin{equation}\label{eq:zeff}
z_{\mathrm{eff}}(t)=\max\{z(t),z_{\min}\},
\end{equation}
with a small reserve threshold $z_{\min}\in(0,1)$ (e.g., $z_{\min}=0.02$) reflecting a practical BMS ``unavailable''
low-SOC region. We then evaluate OCV using $z_{\mathrm{eff}}$:
\begin{equation}\label{eq:voc}
V_{\mathrm{oc}}(z)=E_0-K\Big(\frac{1}{z_{\mathrm{eff}}}-1\Big)+A e^{-B(1-z_{\mathrm{eff}})}.
\end{equation}
\paragraph{(ii) Nonnegative polarization heating.}
Thermal generation is written as $I^2R_0+v_p^2/R_1$, which is always nonnegative and aligns with resistive dissipation
in the polarization branch. This choice avoids sign ambiguities that can arise with alternative $Iv_p$ forms.
\paragraph{(iii) Lightweight current saturation (throttling/PMIC limiting).}
Real devices may throttle performance or limit current under low voltage or high temperature.
We model this with a temperature-dependent current cap:
\begin{equation}\label{eq:I_sat}
I(t)=\min\big(I_{\mathrm{CPL}}(t),\,I_{\max}(T_b(t))\big),
\end{equation}
where a simple continuous form is
\begin{equation}\label{eq:Imax}
I_{\max}(T_b)=I_{\max,0}\Big[1-\rho_T\,(T_b-T_{\mathrm{ref}})\Big]_+,\qquad \rho_T\ge 0.
\end{equation}
When $I_{\mathrm{CPL}}>I_{\max}$, the device operates in a degraded regime with delivered power
$P_{\mathrm{del}}(t)=V_{\mathrm{term}}(t)I(t)\le P_{\mathrm{tot}}(t)$, corresponding to throttling.
\subsection{Initial Conditions and Termination Definitions (TTE and optional $t_\Delta$)}\label{sec:ic_tte}
We use
\begin{equation}\label{eq:ic}
z(0)=z_0,\qquad v_p(0)=0,\qquad T_b(0)=T_a(0),\qquad S(0)=S_0,\qquad w(0)=0.
\end{equation}
We define the time-to-end (time-to-empty / time-to-shutdown) as
\begin{equation}\label{eq:TTE}
\mathrm{TTE}=\inf\Big\{t>0:\ V_{\mathrm{term}}(t)\le V_{\mathrm{cut}}\ \ \text{or}\ \ z(t)\le 0\Big\}.
\end{equation}
Optionally, to quantify CPL infeasibility as a voltage-collapse risk indicator, we define
\begin{equation}\label{eq:tDelta}
t_{\Delta}=\inf\Big\{t>0:\ \Delta(t)\le 0\Big\}.
\end{equation}
With throttling \eqref{eq:I_sat}, $t_\Delta$ is interpreted as the onset time at which pure CPL operation becomes
infeasible, even if the system may continue operating in a degraded mode.
\subsection{Closed-Loop Structure Summary}\label{sec:summary_loop}
The model forms a closed-loop chain:
\begin{equation}\label{eq:loop}
\mathbf{u}(t)\ \Rightarrow\ P_{\mathrm{tot}}(t)\ \Rightarrow\
\big(V_{\mathrm{oc}}(z_{\mathrm{eff}}),R_0(T_b,S),\Delta(t)\big)\ \Rightarrow\
I(t)\ \Rightarrow\ \dot{\mathbf{x}}(t)\ \Rightarrow\ \big(V_{\mathrm{term}}(t),z(t),\mathrm{TTE}\big).
\end{equation}
Nonlinear feedback arises because $P_{\mathrm{tot}}$ is enforced via CPL, while $R_0$ and $Q_{\mathrm{eff}}$
depend on $(T_b,S)$, which in turn evolve under the dissipated power.
\subsection{(Optional) Scaling and Time-Scale Discussion}\label{sec:scaling}
Although not required for computation, a brief scale analysis clarifies stiffness and numerical choices.
Let $\tau_p=R_1C_1$ denote the polarization time constant, and $\tau_{\mathrm{th}}=C_{\mathrm{th}}/(hA)$ the thermal
time constant. Typically $\tau_p\ll \tau_{\mathrm{th}}$, implying fast electrical transients and slower thermal drift.
Moreover, the tail dynamics introduce $\tau_\uparrow\ll \tau_\downarrow$.
These separated time scales motivate a time step that resolves $\tau_p$ and $\tau_\uparrow$ in explicit integration,
as enforced later in the numerical method.

View File

@@ -0,0 +1,392 @@
% =========================
% 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}

View File

@@ -0,0 +1,393 @@
%===========================================================
\section{Uncertainty Quantification and Statistical Inference}
\label{sec:uq}
%===========================================================
This section extends the deterministic continuous-time framework in
Sections~\ref{sec:model_formulation}--\ref{sec:numerics} by modeling future
usage inputs as continuous-time stochastic processes and propagating the
resulting uncertainty through the mechanistic battery model. The objective is
to obtain a \emph{distribution} of time-to-end (TTE) rather than a single-point
estimate, and to quantify the global sensitivity of TTE to key parameters via
variance-based indices. Importantly, the underlying electro-thermal-aging
dynamics and the constant-power-load (CPL) closure are unchanged; randomness
enters only through exogenous inputs and (optionally) uncertain parameters.
%-----------------------------------------------------------
\subsection{Motivation and Model Choices for Random Inputs (OU / Regime Switching)}
\label{subsec:uq_motivation}
%-----------------------------------------------------------
Smartphone usage is intrinsically uncertain beyond a short forecasting horizon:
screen brightness $L(t)$, CPU load $C(t)$, network activity $N(t)$, and signal
quality $\Psi(t)$ exhibit mean-reverting fluctuations, cross-correlations, and
occasional abrupt changes (e.g., screen-off $\to$ gaming; good $\to$ poor
coverage). A purely deterministic extrapolation of $\mathbf{u}(t)$ therefore
tends to understate variability and cannot support probabilistic statements
(e.g., ``runtime exceeds $t$ with $90\%$ probability'').
We model the future input vector
\begin{equation}
\mathbf{u}(t)=[L(t),C(t),N(t),\Psi(t),T_a(t)]^\top
\end{equation}
as a continuous-time stochastic process, while preserving the mechanistic
mapping
$\mathbf{u}(t)\mapsto P_{\mathrm{tot}}(t)\mapsto I(t)\mapsto \dot{\mathbf{x}}(t)$.
Two choices are considered:
\begin{enumerate}
\item \textbf{Bounded multivariate OU (Option U1).}
A multivariate Ornstein--Uhlenbeck process provides mean reversion and
cross-correlation in a continuous-time setting. Smooth bounding transforms
ensure physical admissibility ($L,C,N\in[0,1]$ and $\Psi$ in a prescribed
range).
\item \textbf{Regime-switching OU (Option U2).}
A continuous-time Markov chain $r(t)$ captures discrete ``modes'' (idle,
browsing, video, gaming; good/poor coverage). Within each regime, an OU process
drives the latent inputs. This yields bursty but still continuous trajectories.
\end{enumerate}
Both options are mechanism-compatible and avoid black-box regression: the
battery physics remain deterministic conditional on the sampled input path.
%-----------------------------------------------------------
\subsection{Mathematical Definitions and Bounding Maps}
\label{subsec:uq_definitions}
%-----------------------------------------------------------
\paragraph{Option U1: Bounded multivariate OU.}
Let $\mathbf{y}(t)\in\mathbb{R}^4$ denote latent (unbounded) Gaussian processes
associated with $[L,C,N,\Psi]$. We define
\begin{equation}
d\mathbf{y}(t)=\mathbf{K}\big(\boldsymbol{\mu}-\mathbf{y}(t)\big)\,dt
+\mathbf{\Sigma}\,d\mathbf{W}(t),
\label{eq:mvou}
\end{equation}
where $\mathbf{K}\succ 0$ controls correlation times, $\boldsymbol{\mu}$ is the
long-run mean, $\mathbf{\Sigma}$ sets diffusion intensity, and $\mathbf{W}(t)$
is a standard $4$-dimensional Brownian motion. Cross-channel correlations are
encoded in $\mathbf{\Sigma}\mathbf{\Sigma}^\top$.
Ambient temperature is modeled separately as a scalar OU process:
\begin{equation}
dT_a(t)=k_a\big(\mu_a-T_a(t)\big)\,dt+\sigma_a\,dW_a(t).
\label{eq:ou_ta}
\end{equation}
To enforce physical bounds, we map latent variables to admissible inputs using
a smooth logistic transform $\sigma(s)=(1+e^{-s})^{-1}$:
\begin{align}
L(t)&=\sigma\!\big(y_L(t)\big),\qquad
C(t)=\sigma\!\big(y_C(t)\big),\qquad
N(t)=\sigma\!\big(y_N(t)\big), \label{eq:bound_lcn}\\
\Psi(t)&=\Psi_{\min}+(\Psi_{\max}-\Psi_{\min})\,\sigma\!\big(y_\Psi(t)\big).
\label{eq:bound_psi}
\end{align}
This choice yields continuous trajectories and avoids nonphysical discontinuous
jumps that could artificially trigger the CPL infeasibility condition.
\paragraph{Option U2: Regime-switching OU.}
Let $r(t)\in\{1,\dots,R\}$ be a continuous-time Markov chain with generator
matrix $\mathbf{Q}=[q_{ij}]$, where $q_{ij}\ge 0$ for $j\neq i$ and
$q_{ii}=-\sum_{j\neq i}q_{ij}$. Conditional on $r(t)$, we define
\begin{equation}
d\mathbf{y}(t)=\mathbf{K}_{r(t)}\big(\boldsymbol{\mu}_{r(t)}-\mathbf{y}(t)\big)\,dt
+\mathbf{\Sigma}_{r(t)}\,d\mathbf{W}(t),
\label{eq:rsou}
\end{equation}
and map $\mathbf{y}(t)$ to $\{L,C,N,\Psi\}$ using
Eqs.~\eqref{eq:bound_lcn}--\eqref{eq:bound_psi}. Ambient temperature can also be
regime dependent:
\begin{equation}
dT_a(t)=k_{a,r(t)}\big(\mu_{a,r(t)}-T_a(t)\big)\,dt+\sigma_{a,r(t)}\,dW_a(t).
\label{eq:rsou_ta}
\end{equation}
This formulation captures abrupt mode changes while keeping inputs continuous
between switching times.
%-----------------------------------------------------------
\subsection{Discrete-Time Input Generation (Update Equations)}
\label{subsec:uq_generation}
%-----------------------------------------------------------
For Monte Carlo simulation, we require discrete-time updates over time step
$\Delta t$. For a scalar OU process
\begin{equation}
dy=k(\mu-y)\,dt+\sigma\,dW,
\label{eq:ou_scalar}
\end{equation}
the exact (in distribution) update is
\begin{equation}
y_{n+1}=\mu+(y_n-\mu)e^{-k\Delta t}
+\sigma\sqrt{\frac{1-e^{-2k\Delta t}}{2k}}\,\xi_n,
\qquad \xi_n\sim\mathcal{N}(0,1).
\label{eq:ou_exact}
\end{equation}
For the multivariate OU \eqref{eq:mvou}, one may use the matrix-exponential form
\begin{equation}
\mathbf{y}_{n+1}=\boldsymbol{\mu}+\mathbf{A}\big(\mathbf{y}_n-\boldsymbol{\mu}\big)
+\mathbf{B}\,\boldsymbol{\xi}_n,
\qquad \boldsymbol{\xi}_n\sim\mathcal{N}(\mathbf{0},\mathbf{I}),
\label{eq:mvou_exact}
\end{equation}
where $\mathbf{A}=e^{-\mathbf{K}\Delta t}$ and $\mathbf{B}$ satisfies
$\mathbf{B}\mathbf{B}^\top=\int_0^{\Delta t}e^{-\mathbf{K}s}\mathbf{\Sigma}\mathbf{\Sigma}^\top e^{-\mathbf{K}^\top s}\,ds$.
In practice, choosing $\mathbf{K}$ diagonal yields a simple componentwise
update using \eqref{eq:ou_exact}, while correlations can be retained through
$\mathbf{\Sigma}$.
For regime switching, over sufficiently small $\Delta t$ we approximate
\begin{equation}
\mathbb{P}\big(r_{n+1}=j\,\big|\,r_n=i\big)\approx
\begin{cases}
q_{ij}\Delta t, & j\neq i,\\
1+q_{ii}\Delta t, & j=i,
\end{cases}
\label{eq:ctmc_step}
\end{equation}
then update $\mathbf{y}$ using \eqref{eq:mvou_exact} with parameters associated
with the realized regime $r_n$.
%-----------------------------------------------------------
\subsection{Monte Carlo Propagation and TTE Distribution}
\label{subsec:uq_mc}
%-----------------------------------------------------------
Let $\omega$ denote the randomness driving $\mathbf{u}(t,\omega)$ (and, if
included, uncertain parameters). For each sampled input path $\omega_m$, the
battery dynamics are integrated using the deterministic solver from
Section~\ref{sec:numerics}: RK4 with nested CPL current evaluation at each
substep, including low-SOC OCV protection $z_{\mathrm{eff}}=\max\{z,z_{\min}\}$,
nonnegative polarization heating $v_p^2/R_1$, and the lightweight current cap
$I=\min(I_{\mathrm{CPL}},I_{\max}(T_b))$.
The runtime endpoint is defined by
\begin{equation}
\mathrm{TTE}(\omega)=\inf\left\{t>0:\; V_{\mathrm{term}}(t,\omega)\le V_{\mathrm{cut}}
\ \text{or}\ z(t,\omega)\le 0\right\}.
\label{eq:tte_uq}
\end{equation}
(Optionally, the CPL infeasibility risk time
$t_{\Delta}=\inf\{t>0:\Delta(t,\omega)\le 0\}$ may be recorded as a separate
diagnostic.)
Given $M$ independent sample paths $\{\omega_m\}_{m=1}^M$, we obtain
$\mathrm{TTE}_m=\mathrm{TTE}(\omega_m)$ and form the empirical CDF
\begin{equation}
\widehat{F}_{\mathrm{TTE}}(t)=\frac{1}{M}\sum_{m=1}^M \mathbf{1}\{\mathrm{TTE}_m\le t\}.
\label{eq:emp_cdf}
\end{equation}
The empirical mean and variance are
\begin{equation}
\widehat{\mu}_{\mathrm{TTE}}=\frac{1}{M}\sum_{m=1}^M \mathrm{TTE}_m,\qquad
\widehat{\sigma}^2_{\mathrm{TTE}}=\frac{1}{M-1}\sum_{m=1}^M(\mathrm{TTE}_m-\widehat{\mu}_{\mathrm{TTE}})^2.
\label{eq:tte_moments}
\end{equation}
\paragraph{Monte Carlo error.}
For standard Monte Carlo estimators of smooth functionals of TTE, the
statistical error decays as $O(M^{-1/2})$. We therefore increase $M$ until key
summaries (mean and selected quantiles) stabilize under doubling $M$.
%-----------------------------------------------------------
\subsection{Confidence Intervals, Quantiles, and Survival Curves}
\label{subsec:uq_inference}
%-----------------------------------------------------------
\paragraph{Confidence interval for the mean.}
By the central limit theorem, an approximate $95\%$ confidence interval for the
mean TTE is
\begin{equation}
\widehat{\mu}_{\mathrm{TTE}}\ \pm\ 1.96\,\frac{\widehat{\sigma}_{\mathrm{TTE}}}{\sqrt{M}}.
\label{eq:ci_mean}
\end{equation}
When $M$ is moderate and the distribution is skewed, a nonparametric bootstrap
over $\{\mathrm{TTE}_m\}$ can be used to obtain robust confidence bounds.
\paragraph{Quantiles.}
Let $\mathrm{TTE}_{(1)}\le \cdots \le \mathrm{TTE}_{(M)}$ denote the ordered
samples. The empirical $p$-quantile is
\begin{equation}
\widehat{q}_p=\mathrm{TTE}_{(\lceil pM\rceil)},\qquad p\in(0,1).
\label{eq:quantile}
\end{equation}
In particular, the median is $\widehat{q}_{0.5}$, and the lower-tail quantile
$\widehat{q}_{0.1}$ supports conservative ``guaranteed runtime'' statements.
\paragraph{Survival function.}
A reliability-style summary is the survival curve
\begin{equation}
\widehat{S}(t)=\mathbb{P}(\mathrm{TTE}>t)\approx 1-\widehat{F}_{\mathrm{TTE}}(t).
\label{eq:survival}
\end{equation}
This directly answers: ``what is the probability the phone remains operational
beyond time $t$?''
%-----------------------------------------------------------
\subsection{Variance-Based Global Sensitivity (Sobol Indices)}
\label{subsec:uq_sobol}
%-----------------------------------------------------------
We quantify global parameter importance via variance-based sensitivity indices
for the scalar quantity of interest (QoI)
\begin{equation}
Y=\mathrm{TTE}.
\end{equation}
Let $\boldsymbol{\xi}=(\xi_1,\dots,\xi_d)$ denote uncertain factors (e.g.,
$k_L,\gamma,k_N,\kappa,\mu_a$ and other parameters as needed), assumed
independent with prescribed prior distributions. Because usage randomness
$\omega$ also contributes variance, we recommend defining the QoI as the
\emph{conditional expectation} over usage paths:
\begin{equation}
Y(\boldsymbol{\xi})=\mathbb{E}_{\omega}\big[\mathrm{TTE}(\boldsymbol{\xi},\omega)\big],
\label{eq:qoi_condexp}
\end{equation}
which yields stable and actionable sensitivities to design/physics parameters.
In computations, \eqref{eq:qoi_condexp} is approximated by an inner Monte Carlo
average over $M_{\omega}$ usage realizations.
The first-order Sobol index of factor $\xi_i$ is defined as
\begin{equation}
S_i=\frac{\mathrm{Var}\big(\mathbb{E}[Y\mid \xi_i]\big)}{\mathrm{Var}(Y)},
\label{eq:sobol_first}
\end{equation}
and the total-effect index is
\begin{equation}
S_{T_i}=1-\frac{\mathrm{Var}\big(\mathbb{E}[Y\mid \boldsymbol{\xi}_{\sim i}]\big)}{\mathrm{Var}(Y)},
\label{eq:sobol_total}
\end{equation}
where $\boldsymbol{\xi}_{\sim i}$ denotes all factors except $\xi_i$. Large
$S_i$ indicates a strong main effect, while a large gap $S_{T_i}-S_i$ indicates
substantial interaction and/or nonlinearity (expected here due to CPL feedback
and electro-thermal coupling).
%-----------------------------------------------------------
\subsection{Saltelli Sampling and Estimation}
\label{subsec:uq_saltelli}
%-----------------------------------------------------------
We employ the Saltelli sampling scheme for efficient estimation of Sobol
indices. Let $\mathbf{A},\mathbf{B}\in\mathbb{R}^{N\times d}$ be two independent
sample matrices of $\boldsymbol{\xi}$. For each $i\in\{1,\dots,d\}$, construct
$\mathbf{A}^{(i)}_B$ by replacing the $i$-th column of $\mathbf{A}$ with the
$i$-th column of $\mathbf{B}$. Denote the corresponding model evaluations by
\begin{equation}
Y_A^{(n)}=Y(\mathbf{A}_n),\quad
Y_B^{(n)}=Y(\mathbf{B}_n),\quad
Y_{A^{(i)}_B}^{(n)}=Y(\mathbf{A}^{(i)}_{B,n}),
\qquad n=1,\dots,N.
\end{equation}
We estimate $\mathrm{Var}(Y)$ from the pooled samples and compute Sobol
estimators in the following commonly used form:
\begin{equation}
\widehat{S}_i=
\frac{\frac{1}{N}\sum_{n=1}^N Y_B^{(n)}\left(Y_{A^{(i)}_B}^{(n)}-Y_A^{(n)}\right)}
{\widehat{\mathrm{Var}}(Y)},
\label{eq:saltelli_first}
\end{equation}
\begin{equation}
\widehat{S}_{T_i}=
\frac{\frac{1}{2N}\sum_{n=1}^N \left(Y_A^{(n)}-Y_{A^{(i)}_B}^{(n)}\right)^2}
{\widehat{\mathrm{Var}}(Y)}.
\label{eq:saltelli_total}
\end{equation}
\paragraph{Nested averaging over usage paths.}
Each $Y(\cdot)$ above is computed as
\begin{equation}
Y(\boldsymbol{\xi})\approx \frac{1}{M_{\omega}}\sum_{m=1}^{M_{\omega}}
\mathrm{TTE}(\boldsymbol{\xi},\omega_m),
\label{eq:nested_mc}
\end{equation}
where $\{\omega_m\}$ are i.i.d.\ usage realizations generated by
Option~U1/U2. This inner average reduces the Monte Carlo noise in $Y$ so that
the outer Saltelli estimators converge reliably in $N$.
%-----------------------------------------------------------
\subsection{Optional: Variance Reduction (LHS / Quasi-Monte Carlo)}
\label{subsec:uq_varred}
%-----------------------------------------------------------
While plain Monte Carlo converges at rate $O(M^{-1/2})$, variance reduction can
improve efficiency when computational budgets are tight.
\paragraph{Latin hypercube sampling (LHS).}
For estimating the TTE distribution under uncertain inputs/parameters, LHS can
replace i.i.d.\ sampling of low-dimensional uncertain parameters
$\boldsymbol{\xi}$ to reduce estimator variance without changing the model. LHS
is especially effective when the dominant uncertainty is parameter-driven.
\paragraph{Quasi-Monte Carlo (QMC).}
For Sobol estimation (outer sampling), low-discrepancy sequences (e.g., Sobol
sequences) can improve convergence of integral estimates in moderate
dimensions. In this work, QMC can be applied to generate $\mathbf{A},\mathbf{B}$
before constructing $\mathbf{A}^{(i)}_B$. Because our QoI involves a nested
average \eqref{eq:nested_mc}, QMC primarily benefits the outer parameter
integration, while the inner usage randomness still scales as
$O(M_\omega^{-1/2})$.
\paragraph{Control variates (conceptual).}
If a simplified surrogate (e.g., the same model with fixed $T_b=T_a$ or without
aging) is available, it may serve as a control variate to reduce variance of
$\mathrm{TTE}$. We do not rely on this technique in the baseline pipeline.
%-----------------------------------------------------------
\subsection{Optional: Unified Two-Level Uncertainty (Inputs and Parameters)}
\label{subsec:uq_twolevel}
%-----------------------------------------------------------
When both usage inputs and physical/power parameters are uncertain, the full
QoI can be viewed hierarchically as
\begin{equation}
\mathrm{TTE}=\mathrm{TTE}(\boldsymbol{\xi},\omega),
\end{equation}
with $\boldsymbol{\xi}$ representing uncertain parameters (e.g., $k_L,\gamma,
k_N,\kappa,\mu_a,hA$) and $\omega$ representing stochastic input realizations
from Option~U1/U2. Two complementary summaries are useful:
\paragraph{Unconditional runtime distribution.}
The overall distribution integrates over both sources of uncertainty:
\begin{equation}
F_{\mathrm{TTE}}(t)=\mathbb{P}(\mathrm{TTE}\le t)=
\int \mathbb{P}\!\left(\mathrm{TTE}(\boldsymbol{\xi},\omega)\le t\ \big|\ \boldsymbol{\xi}\right)
\,p(\boldsymbol{\xi})\,d\boldsymbol{\xi}.
\label{eq:unconditional}
\end{equation}
This is estimated by outer sampling of $\boldsymbol{\xi}$ and inner sampling of
$\omega$.
\paragraph{Sensitivity of conditional mean runtime.}
For design guidance, sensitivities are computed for
$Y(\boldsymbol{\xi})=\mathbb{E}_{\omega}[\mathrm{TTE}(\boldsymbol{\xi},\omega)]$
as in \eqref{eq:qoi_condexp}, yielding Sobol indices that reflect how parameter
variation shifts \emph{expected} runtime under random usage.
\paragraph{Practical computation.}
A computationally efficient compromise is to (i) propagate usage uncertainty
with a large $M$ at nominal parameters to obtain $F_{\mathrm{TTE}}$, and (ii)
compute Sobol indices with moderate inner averaging $M_\omega$ and outer sample
size $N$ to rank parameter importance.
%-----------------------------------------------------------
\subsection*{Algorithmic Summary}
%-----------------------------------------------------------
For completeness, the full UQ pipeline used in subsequent sections can be
summarized as follows:
\begin{itemize}
\item Generate stochastic input paths $\mathbf{u}(t,\omega)$ using
Eqs.~\eqref{eq:mvou}--\eqref{eq:bound_psi} (Option~U1) or
Eqs.~\eqref{eq:rsou}--\eqref{eq:rsou_ta} (Option~U2), with discrete updates
given by \eqref{eq:ou_exact}--\eqref{eq:ctmc_step}.
\item For each path, solve the mechanistic battery model using RK4 with nested
CPL current evaluation (Section~\ref{sec:numerics}) and record
$\mathrm{TTE}$ from \eqref{eq:tte_uq}.
\item Construct the empirical distribution \eqref{eq:emp_cdf}, compute moments
\eqref{eq:tte_moments}, confidence intervals \eqref{eq:ci_mean}, quantiles
\eqref{eq:quantile}, and survival curve \eqref{eq:survival}.
\item For global sensitivity, evaluate $Y(\boldsymbol{\xi})$ via nested averaging
\eqref{eq:nested_mc} and estimate Sobol indices with Saltelli sampling
\eqref{eq:saltelli_first}--\eqref{eq:saltelli_total}.
\end{itemize}