17 KiB
%=========================================================== \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} (OptionU1) or
Eqs.\eqref{eq:rsou}--\eqref{eq:rsou_ta} (OptionU2), 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}