%===========================================================
\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}
