# Graph-accelerated non-intrusive polynomial chaos expansion using partially
tensor-structured quadrature rules

Bingran Wang111Ph.D Candidate, Department of Mechanical and Aerospace
Engineering University of California, San Diego, La Jolla, CA, USA  Nicholas
C. Orndorff222Ph.D Student, Department of Mechanical and Aerospace Engineering
and John T. Hwang333Assistant Professor, Department of Mechanical and
Aerospace Engineering University of California, San Diego, La Jolla, CA, USA

###### Abstract

Recently, the graph-accelerated non-intrusive polynomial chaos (NIPC) method
has been proposed for solving uncertainty quantification (UQ) problems. This
method leverages the full-grid integration-based NIPC method to address UQ
problems while employing the computational graph transformation approach,
AMTC, to accelerate the tensor-grid evaluations. This method exhibits
remarkable efficacy on a broad range of low-dimensional (three dimensions or
less) UQ problems featuring multidisciplinary models. However, it often does
not scale well with problem dimensions due to the exponential increase in the
number of quadrature points when using the full-grid quadrature rule. To
expand the applicability of this method to a broader range of UQ problems,
this paper introduces a new framework for generating a tailored, partially
tensor-structured quadrature rule to use with the graph-accelerated NIPC
method. This quadrature rule, generated through the designed quadrature
approach, possesses a tensor structure that is tailored for the computational
model. The selection of the tensor structure is guided by an analysis of the
computational graph, ensuring that the quadrature rule effectively capitalizes
on the sparsity within the computational graph when paired with the AMTC
method. This method has been tested on one 4D and one 6D UQ problem, both
originating from aircraft design scenarios and featuring multidisciplinary
models. Numerical results show that, when using with graph-accelerated NIPC
method, our approach generates a partially tensor-structured quadrature rule
that outperforms the full-grid Gauss quadrature and the designed quadrature
methods (more than 40% reduction in computational costs) in both of the test
problems.

##  1 Introduction

Engineers rely heavily on computational models to design, analyze, and
optimize engineering systems. However, such computational models are often
deterministic and inevitably affected by uncertainties inherent in real-world
scenarios. Uncertainty quantification (UQ) plays a pivotal role in ensuring
the reliability of predictions made by these computational models. Its
applications span various scientific and engineering domains, including
weather forecasting [1, 2], machine learning [3, 4]), structural analysis [5,
6], and aircraft design [7, 8, 9], where uncertainties can significantly
impact system behavior and performance.

Forward UQ, also known as uncertainty propagation, seeks to assess how
uncertainties associated with model inputs affect a computational model’s
outputs or quantities of interest (QoIs). These input uncertainties can stem
from variations in operating conditions and model parameters, referred to as
aleatoric uncertainties, or from model errors due to incomplete knowledge,
termed epistemic uncertainties. By quantifying the impact of input
uncertainties, forward UQ plays a crucial role in supporting decision-making
and risk assessment in engineering system design. In this paper, we address a
UQ problem within a probabilistic framework, representing uncertain inputs as
continuous probability density distributions. The objective is to estimate the
statistical moments or risk measures of the QoI.

We focus on employing the non-intrusive polynomial chaos (NIPC) method [10] to
solve the UQ problem, a popular approach especially suitable for UQ problems
involving fewer than 10 dimensions. The NIPC method approximates the model
response through a linear combination of polynomial chaos expansion (PCE)
basis functions, determined by the distributions of the random inputs. The PCE
coefficients are computed using either sampling methods [11, 12, 13] or
integration methods [14, 15]. In [16], Eldred and Burkardt demonstrated that
for low-dimensional UQ problems, the integration-based and regression-based
NIPC methods show comparable performance, with significantly higher
convergence rates than the Monte Carlo method.

The integration approach makes use of the orthogonality properties of the PCE
basis functions and transforms the UQ problems into multidimensional
integration problems. Common methods for solving the integration problems
include Gauss quadrature [17], Smolyak sparse-grid [18], and designed
quadrature [19] methods. The Gauss quadrature method formulates the quadrature
rule using a tensor product of univariate quadrature rules. This approach is
easy to use but the number of quadrature points increases exponentially with
the number of dimensions. The Smolyak sparse-grid approach breaks the tensor
structure of the Gauss quadrature rule by dropping the higher-order cross
terms while achieving minimal sacrifice on its accuracy. In comparison, the
design quadrature method generates the quadrature rule by solving an
optimization problem to find approximate solutions to the moment-matching
equations. The designed quadrature method results in a quadrature rule that
does not possess a tensorial structure but can typically achieve better
performance than the previous two methods.

Wang et al. recently introduced a new method called Accelerated Model
evaluations on Tensor grids using Computational graph transformations (AMTC)
[20]. This method reduces the computational cost of model evaluations on
tensor-grid inputs by modifying the computational graph of the model within
the middle end of a modeling language’s three-stage compiler. This
modification eliminates redundant evaluations at the operation level caused by
the tensor structure of the inputs. Their previous work applies AMTC with
integration-based NIPC to solve UQ problems, under the assumption of a
specific type of sparsity within the computational graphs of the computational
models. This integration is termed the graph-accelerated NIPC method and has
demonstrated superior efficiency in various UQ and optimization under
uncertainty problems involving multidisciplinary systems [8]. However, the
graph-accelerated NIPC method has so far only been implemented with a full-
grid Gauss quadrature rule, targeting low-dimensional UQ problems.

This paper introduces a novel framework aimed at producing a quadrature rule
possessing a specified tensor structure, tailored for use with the graph-
accelerated NIPC method in addressing the UQ problem. Within this framework,
we first identify a suitable tensor structure option for the quadrature rule
through an analysis of the computational model. Subsequently, we employ the
designed quadrature method to construct the quadrature rule that adheres to
this tensor structure. The resulting quadrature rule is strategically designed
to leverage the sparsity inherent in the computational graph of the model,
thus maximizing efficiency when employing the AMTC method for accelerating
tensor-grid evaluations.

This method has been applied to two UQ problems, one with four uncertain
inputs and one with six uncertain inputs, both derived from practical aircraft
design problems involving multidisciplinary systems. The proposed method
produces new partially tensor-structured quadrature rules distinct from
existing ones for both test problems. These partially structured quadrature
rules also outperform existing methods when employed with the graph-
accelerated NIPC method for solving UQ problems.

This paper is structured as follows. Section 2 provides background information
on NIPC, quadrature rules, and AMTC. In Section 3, we elaborate on the
proposed framework. Section 4 presents the numerical results obtained from the
test problems. Finally, Section 5 summarizes the work and provides concluding
remarks.

##  2 Background

###  2.1 Non-intrusive polynomial chaos

Polynomial chaos, also known as polynomial chaos expansion (PCE), traces its
origins back to Wiener, who utilized Hermite polynomials in developing a
theoretical framework for Gaussian random processes [21]. Ghanem further
expanded upon this concept, employing Hermite polynomials as orthogonal bases
to represent random processes [22]. However, this method suffers from
convergence issues for problems with non-Gaussian uncertain inputs. This issue
was addressed by Xiu and Karniadakis in their work [10], wherein they
introduced the generalized polynomial chaos (gPC) method. The gPC method
strategically employs different orthogonal polynomials tailored to uncertain
inputs with different distributions, thereby achieving optimal convergence in
addressing UQ problems.

We examine a general UQ problem involving a function represented as

$f=\mathcal{F}(u),$ (1)

where $\mathcal{F}:\mathbb{R}^{d}\to\mathbb{R}$ represents the model
evaluation function, $u\in\mathbb{R}^{d}$ denotes the input vector, and
$f\in\mathbb{R}$ represents a scalar output. The UQ problem aims to compute
the statistical moments or risk measures of the output random variable $f(U)$
under the effect of the input uncertainties, denoted as a stochastic vector,
$U:=[U_{1},\ldots,U_{d}]$. The input uncertainties are assumed to be mutually
independent and adhere to the probability distribution $\rho(u)$ with support
$\Gamma$.

Using the gPC method, the stochastic output can be represented as a linear
combination of orthogonal polynomials:

$f(U)=\sum_{i=0}^{\infty}\alpha_{i}\Phi_{i}(U),$ (2)

where $\Phi_{i}(U)$ denotes the PCE basis functions that are chosen based on
the uncertain inputs’ probability distribution $\rho(u)$, and $\alpha_{i}$
represents the PCE coefficients that need to be determined. In practical
applications, it is common to truncate the PCE basis functions up to a
specific polynomial order and approximate the stochastic output as

$f(U)\approx\sum_{i=0}^{q}\alpha_{i}\Phi_{i}(U).$ (3)

With $p$ as the upper bound for the total polynomial order, the resultant
number of PCE basis functions, $q+1$, satisfies

$q+1=\frac{(d+p)!}{d!p!}.$ (4)

The PCE basis functions are chosen to satisfy the orthogonality property as

$\left<\Phi_{i}(U),\Phi_{j}(U)\right>=\delta_{ij},$ (5)

where $\delta_{ij}$ is the Kronecker delta and the inner product is defined as

$\left<\Phi_{i}(U),\Phi_{j}(U)\right>=\int_{\Gamma}\Phi_{i}(u)\Phi_{j}(u)\rho(u)du.$
(6)

In the one-dimensional case, the PCE basis functions consist of univariate
orthogonal polynomials chosen according to the distribution of the uncertain
input. Table 1 presents common types of continuous random variables alongside
their corresponding orthogonal polynomials and support ranges. In multi-
dimensional cases, the PCE basis functions are formed as tensor products of
the univariate orthogonal polynomials. For instance, in the 2D case, the first
six PCE basis functions are expressed as

$\displaystyle\Phi_{0}(U)=\phi_{0}(U_{1})\phi_{0}(U_{2})$
$\displaystyle\Phi_{1}(U)=\phi_{1}(U_{1})\phi_{0}(U_{2})$
$\displaystyle\Phi_{2}(U)=\phi_{0}(U_{1})\phi_{1}(U_{2})$
$\displaystyle\Phi_{3}(U)=\phi_{2}(U_{1})\phi_{0}(U_{2})$
$\displaystyle\Phi_{4}(U)=\phi_{1}(U_{1})\phi_{1}(U_{2})$
$\displaystyle\Phi_{5}(U)=\phi_{0}(U_{1})\phi_{2}(U_{2}),$

where $\phi_{i}(U_{j})$ represents the $i$-th univariate orthogonal polynomial
corresponding to the uncertain input, $U_{j}$. If the uncertain inputs are
dependent, the PCE basis functions can be generated numerically following the
approach outlined in [23].

Using the gPC approximation in 3, the statistical moments of the QoI can be
analytically computed as

$\mu_{f}=\alpha_{0},$ (7)
$\sigma_{f}=\sum_{i=i}^{d}\alpha_{i}^{2}\left<\Phi_{i}^{2}\right>.$ (8)

For risk measures like probability of failure or Conditional Value at Risk
(CVaR), the gPC model can be used as an inexpensive-to-evaluate surrogate
model, and one can apply sampling methods to approximate them.

Table 1: Orthogonal polynomials for common types of continuous random variables Distribution | Orthogonal polynomials | Support range
---|---|---
Normal | Hermite | $(-\infty,\infty)$
Uniform | Legendre | $[-1,1]$
Exponential | Laguerre | $[0,\infty)$
Beta | Jacobi | $(-1,1)$
Gamma | Generalized Laguerre | $[0,\infty)$

To solve the PCE coefficients $\alpha_{i}$ in (3), the non-intrusive
polynomial chaos (NIPC) method applies either the integration or regression
approach. The integration approach makes use of the orthogonality property in
(5) and projects the QoI onto each basis function, which results in a multi-
dimensional integration problem:

$\displaystyle\alpha_{i}$
$\displaystyle=\frac{\left<f(U),\Phi_{i}\right>}{\left<\Phi_{i}^{2}\right>}$
(9)
$\displaystyle=\frac{1}{\left<\Phi_{i}^{2}\right>}\int_{\Gamma}f(u)\Phi_{i}(u)\rho(u)du$
$\displaystyle=\frac{1}{\left<\Phi_{i}^{2}\right>}\int_{\Gamma_{1}}\ldots\int_{\Gamma_{d}}f(u_{1},\ldots,u_{d})\Phi_{i}(u_{1},\ldots,u_{d})\rho(u_{1},\ldots,u_{d})du_{1}\ldots
du_{d}.$

In practical problems, this integration problem is often approximated using
numerical quadrature rules. Conversely, the regression approach computes the
PCE coefficients by solving a linear least-squares problem:

$\begin{bmatrix}\Phi_{1}(u^{(1)})&\ldots&\Phi_{q}(u^{(1)})\\\
\vdots&&\vdots\\\ \Phi_{1}(u^{(n)})&\ldots&\Phi_{q}(u^{(n)})\\\
\end{bmatrix}\begin{bmatrix}\alpha_{1}\\\ \vdots\\\
\alpha_{q}\end{bmatrix}=\begin{bmatrix}f(u^{(1)})\\\ \vdots\\\
f(u^{(n)})\end{bmatrix},$ (10)

where $(u^{(1)},\ldots,u^{(n)})$ represents a set of $n$ sample points
generated from the distributions of the uncertain inputs. A general rule of
thumb dictates that the number of samples $n$ should be 2-3 times the number
of coefficients. Consequently, the cost of the regression-based NIPC method
can be reduced by employing sparse PCE basis functions, resulting in fewer
coefficients to be estimated. For those interested in sparse PCE, a recent
review paper can be found in [24].

###  2.2 Numerical quadrature rules

This paper focuses on employing the integration-based NIPC method for solving
multi-dimensional UQ problems. In this context, the integration-based NIPC
method is essentially solving the integration problem as in (9). We denote the
integration problem as

$I(f)=\int_{\Gamma}f(u)\rho(u)du,$

where $u=(u_{1},\ldots,u_{d})\in\mathbb{R}^{d}$. Solving this integration
problem often involves utilizing the numerical quadrature approach, which
approximates the integral as a weighted sum of function evaluations at
specific quadrature points:

$I(f)\approx\sum_{i=1}^{n}w^{(i)}f(u^{(i)}),$ (11)

where $u^{(i)}\in\Gamma$ and $w^{(i)}>0$ represent the nodes and weights,
respectively, to be determined by the quadrature rule.

The objective of the quadrature rule is to effectively approximate a large set
of functions with a minimal number of function evaluations. This is typically
achieved by ensuring equality conditions hold for all functions $f$ within a
polynomial subspace $\Pi$, as expressed by:

$\int_{\Gamma}f(u)\rho(u)du=\sum_{i=1}^{n}w^{(i)}f(u^{(i)})\quad\text{for
all}\quad f\in\Pi_{r}.$ (12)

Here, we consider the cross-polynomial subspace with a total polynomial order
of $r$, defined as:

$\Pi_{r}=\text{span}\\{u^{\alpha}\ |\ \alpha\in\Lambda_{r}\\},$ (13)

where

$\alpha=(\alpha_{1},\ldots,\alpha_{d}),\quad
u^{\alpha}=\prod_{j=1}^{d}(u_{j})^{\alpha_{j}}\quad$ (14)

and

$\Lambda_{r}=\left\\{\alpha\mid\sum_{i=1}^{d}\alpha_{i}\leq r\right\\}.$ (15)

In many numerous integration problems, the evaluation functions exhibit smooth
behavior and can be accurately approximated by the polynomials within the
cross-polynomial subspace. In such instances, the quadrature rule enforcing
Equation (12) can prove to be highly effective.

####  2.2.1 Gauss quadrature rule

One approach to generating the quadrature points is through the use of the
Gauss quadrature rule [17]. In the one-dimensional case, the quadrature points
$u^{(1)},\ldots,u^{(k)}$ are the roots of the orthogonal polynomial $p_{k}(u)$
with degree $k$. The polynomials here are defined in the same way as the PCE
basis functions, and they also satisfy orthogonality:

$\left<p_{i},p_{j}\right>=\int_{\Gamma}p_{i}(u)p_{j}(u)\rho(u)du=\delta_{ij}.$
(16)

The weights of the Gauss quadrature rule, $w^{(1)},\ldots,w^{(k)}$ can be
computed by solving the moment-matching equations:

$\sum_{i=1}^{k}p_{j}\left(u^{(i)}\right)w^{(i)}=\left\\{\begin{array}[]{ll}1&\text{
if }j=0\\\ 0&\text{ for }j=1,\ldots,k-1.\end{array}\right.$ (17)

In this approach, the univariate Gauss quadrature rule with $k$ nodes and
weights can exactly integrate the univariate polynomials up to degree $2k-1$.

For multi-dimensional integrations, the Gauss quadrature rule extends the
univariate rule to form a tensor structure, Assuming $k$ quadrature points in
each dimension, the integral is approximated as

$I\approx\sum_{i_{1}=0}^{k}\ldots\sum_{i_{d}=0}^{k}w_{1}^{(i_{1})}\ldots
w_{d}^{(i_{d})}f(u_{1}^{(i_{1})},\ldots,u_{d}^{(i_{d})})\rho(u_{1}^{(i_{1})},\ldots,u_{n}^{(i_{d})}),$
(18)

where $(u_{i}^{(1)},\ldots,u_{i}^{(k)})$ and
$(w_{i}^{(1)},\ldots,w_{i}^{(k)})$ are the nodes and weights for the 1D Gauss
quadrature rule in the $u_{i}$ dimension, respectively. We represent the full-
grid quadrature points as

$\boldsymbol{u}=\boldsymbol{u}^{k}_{1}\times\cdots\times\boldsymbol{u}^{k}_{d},$
(19)

where $\boldsymbol{u}^{k}_{i}:=\\{u_{i}^{(j)}\\}_{j=1}^{k}$. Indeed, while
this method allows exact integration of cross-polynomials up to a total order
of $2k-1$, it suffers from the curse of dimensionality. This is because the
total number of quadrature points increases exponentially with the number of
dimensions, as $n=k^{d}$.

To mitigate the curse of dimensionality the Gauss quadrature rule suffers
from, the Smolyak sparse grid method [18] offers a solution. This method
strategically drops higher-order cross terms in the full-grid Gauss quadrature
points, reducing the number of required quadrature points while preserving the
accuracy of the quadrature rule to a significant extent. The sparse grid
quadrature points can be expressed as:

$\boldsymbol{u}=\bigcup_{\ell-d+1\leq|\mathbf{i}|\leq\ell}\left(\boldsymbol{u}^{i_{1}}_{1}\times\cdots\times\boldsymbol{u}^{i_{n}}_{d}\right),$
(20)

where $\ell$ is the level of construction.

####  2.2.2 Designed quadrature method

For high-dimensional integration problems, a more effective quadrature rule is
offered by the designed quadrature method proposed by Keshavarzzadeh et al.
[19]. The core idea of the designed quadrature method is to numerically
generate a quadrature rule that satisfies the multi-variate moment-matching
equations:

$\sum_{i=1}^{n}\Phi_{\boldsymbol{i^{\prime}}}\left(u^{(i)}\right)w^{(i)}=\left\\{\begin{array}[]{ll}1&\text{
if }|\boldsymbol{i^{\prime}}|=0\\\ 0&\text{ for }\alpha\in
0<|\boldsymbol{i^{\prime}}|<k,\end{array}\right.$ (21)

where

$\pi_{\alpha}=\prod_{i=1}^{d}p_{\alpha_{i}}^{(i)}(u_{i})$ (22)

and $p^{(i)}_{j}$ is the $j$-th orthogonal polynomial in the $u_{i}$
dimension. If the nodes and weights of the quadrature rule satisfy the
equations in (21), then the quadrature rules can exactly integrate cross-term
polynomials up to $(2k-1)$th order. We denote the quadrature points and the
weights as $\boldsymbol{u}:=(u^{(1)},\ldots,u^{(n)})$ and the weights
$\boldsymbol{w}:=(w^{(1)},\ldots,w^{(n)})$ respectively. Then, the moment-
matching equations in (21) can be written as residual equations:

$\mathcal{R}(\boldsymbol{u},\boldsymbol{w})=0,$ (23)

where $\mathcal{R}:\mathbb{R}^{d\times
n}\times\mathbb{R}^{n}\to\mathbb{R}^{n}$.

The designed quadrature method approximates the solution of the moment-
matching equations by solving an optimization problem formulated as

$\displaystyle\min_{\boldsymbol{u},\boldsymbol{w}}$
$\displaystyle\|\mathcal{R}(\boldsymbol{u},\boldsymbol{w})\|_{2}$ (24) subject
to $\displaystyle u^{(j)}\in\Gamma,\quad j=1,\ldots,n,$ $\displaystyle
w^{(j)}>0,\quad j=1,\ldots,n.$

The objective of this optimization, as shown in Equation (24), is to minimize
the $L^{2}$ norm of the residual equations while constraining the nodes to lie
within the support range and ensuring that the weights are positive. This
method has demonstrated superior effectiveness compared to the sparse-grid
Gauss quadrature rule for many high-dimensional integration problems. The
number of quadrature points required is typically chosen as
$n=\eta\frac{(2d)^{k-1}}{(k-1)!}$, where $\eta\in[0.5,0.9]$ serves as a tuning
parameter.

###  2.3 Computational graph transformations for efficient tensor-grid
evaluations

Figure 1: Computational graph of function $f=cos(u_{1})+exp(-u_{2})$[20]

(a)

(b)

Figure 2: Computational graphs with data size for full-grid input points
evaluation on $f=cos(u_{1})+exp(-u_{2})$ [20]

Accelerated Model evaluations on Tensor grids using Computational graph
transformations (AMTC) is a computational graph transformation method,
recently introduced by Wang et al. in [20], with its core idea initially
proposed in [25]. AMTC accelerates the evaluation time of a computational
model on tensor-grid inputs by eliminating repeated evaluations on the
operation level. Any computational model can be represented as a computational
graph with basic operations. In a conventional for-loop approach, when
evaluating a computational model on tensor-grid inputs, each operation on the
computational graph needs to be evaluated at every input point. However, with
tensor-grid inputs, where inputs possess a fully tensorial structure, only a
small set of unique values exists in each input space. Consequently, many
operations within the computational model might depend on only a few inputs,
leading to many redundant repeated evaluations at the operation level with a
for-loop approach. To address this inefficiency, AMTC eliminates the redundant
evaluations by generating a modified computational graph. The modified
computational graph is partitioned into smaller sub-graphs based upon on which
inputs each operation depends. Operations within each sub-graph share the same
input space and are solely evaluated on distinct nodes within their input
space. To ensure proper data flow between these sub-graphs, Einstein summation
(Einsum) operations are inserted at each connection joint between sub-graphs.

To consider an illustrative example, assume we want to evaluate a simple
function,

$f=\cos(u_{1})+\exp(-u_{2}),$ (25)

at full-grid input points with $k$ input points at each dimension. The given
function can be represented as a sequence of elementary operations:

$\displaystyle\xi_{1}=\varphi_{1}(u_{1})=\cos(u_{1});$ (26)
$\displaystyle\xi_{2}=\varphi_{2}(u_{2})=-u_{2};$
$\displaystyle\xi_{3}=\varphi_{3}(\xi_{2})=\exp(\xi_{2});$ $\displaystyle
f=\varphi_{4}(\xi_{1},\xi_{3})=\xi_{1}+\xi_{3},$

and its computational graph is depicted in Fig.1. The tensor-grid input points
are denoted as

$\boldsymbol{u}=\begin{Bmatrix}(u_{1}^{(1)},u_{2}^{(1)})&\ldots&(u_{1}^{(k)},u_{2}^{(1)})\\\
\vdots&\ddots&\vdots\\\
(u_{1}^{(1)},u_{2}^{(k)})&\ldots&(u_{1}^{(k)},u_{2}^{(k)})\\\ \end{Bmatrix}.$
(27)

The computational graphs to perform these full-grid evaluations with and
without the AMTC method are shown in Fig. 2. In comparison, without AMTC, each
operation in the computational graph is evaluated $k^{2}$ times, corresponding
to the total number of input points. However, by applying AMTC, operations
that depend only on $u_{1}$ and $u_{2}$ are only evaluated $k$ times since
their outputs have only $k$ distinct values.

AMTC has been integrated into the middle end of the compiler for the
Computational System Design Language (CSDL) [26], a recently developed domain-
specific language tailored for tackling multidisciplinary design, analysis,
and optimization problems. CSDL has gained significant popularity in
addressing large-scale multidisciplinary design optimizations, with
applications including designing novel air-taxi [27, 28] and spacecraft [29]
concepts. Within the CSDL compiler, the AMTC serves as a middle-end algorithm,
automatically generating the modified computational graph to minimize the
model evaluation cost on tensor-grid inputs.

This method has been used in combination with the integration-based NIPC
method, which forms full-grid quadrature rules for low-dimensional UQ
problems. We refer to the combination of integration-based NIPC and AMTC as
graph-accelerated NIPC, and this method has demonstrated exceptional
effectiveness in solving certain low-dimensional UQ problems involving
multidisciplinary or multi-point systems [20].

##  3 Methodology

This paper proposes a new framework that generates a tensor-structured
quadrature rule that is specifically tailored for a given computational model.
This quadrature rule is intended to be used with the graph-accelerated NIPC
method to efficiently solve the UQ problem.

###  3.1 Generation of tensor-structured quadrature rules

We use the degree of polynomials that a quadrature rule can exactly integrate
as the measure of its accuracy. Under this criterion, if we maintain a non-
tensorial structure for the quadrature rule, the most effective strategy is to
utilize the designed quadrature method. In this scenario, the quadrature
points can be expressed as:

$\boldsymbol{u}=\boldsymbol{u}^{n}_{\\{1,\ldots,d\\}}.$ (28)

The nodes and weights in the quadrature rule are determined by solving the
optimization problem outlined in (24), and the number of quadrature points is
chosen as

$n=\eta\frac{(2d)^{k-1}}{(k-1)!},\quad\eta\in[0.5,0.9].$ (29)

Adopting this method ensures that the resulting quadrature rule can
effectively integrate polynomials up to $(2k-1)$th order. If a fully tensorial
structure for the quadrature points is necessary, the optimal method to employ
is the full-grid Gauss quadrature method. In this approach, the quadrature
points are expressed as:

$\boldsymbol{u}=\boldsymbol{u}^{k}_{\\{1\\}}\times\cdots\times\boldsymbol{u}^{k}_{\\{d\\}},$
(30)

where $\boldsymbol{u}^{k}_{\\{i\\}}$ represents the 1D Gauss quadrature points
in $u_{i}$ dimension. The sets of 1D Gauss quadrature points in each dimension
are determined separately to ensure the nodes and weights can exactly
integrate the univariate polynomials up to $(2k-1)$th order. Maintaining a
fully tensorial structure of the 1D Gauss quadrature points and weights
ensures that the full-grid quadrature rule can exactly integrate cross-
polynomials up to the $(2k-1)$th order, albeit while requiring $n=k^{d}$
quadrature points. In fact, in each dimension, the 1D Gauss quadrature points
and weights correspond to the global minimum of the optimization problem
outlined in (24). Assuming that the optimization problem in (24) always finds
the global minimum, employing the full-grid Gauss quadrature rule is
equivalent to applying the designed quadrature method in each dimension and
subsequently forming a fully tensorial structure between the 1D solutions.

In addition to the fully tensorial and non-tensorial structure options, we can
also preserve a partially tensorial structure for the quadrature points while
ensuring accurate integration of all polynomials up to a specified order. For
instance, in a 4D integration problem, we could maintain a partially tensorial
structure as:

$\boldsymbol{u}=\boldsymbol{u}^{n_{1}}_{\\{1,2\\}}\times\boldsymbol{u}^{n_{2}}_{\\{3,4\\}},$
(31)

where a tensorial structure is retained between the quadrature points in the
space of $u_{1},u_{2}$ and the quadrature points in the space of
$u_{3},u_{4}$. If the quadrature rules in each space are selected such that
the quadrature points can precisely integrate polynomials up to the same
order, then this partially tensorial quadrature rule can also precisely
integrate the polynomials to the same order in the space of $u$. The total
number of tensorial structure options for the quadrature rule can be
represented by the Bell number from combinatorial mathematics, which is given
by:

$B_{d}=\sum_{i=0}^{d}\left\\{\begin{array}[]{l}i\\\ d\end{array}\right\\}.$
(32)

For example, for $d=3$, $B_{3}=5$, the five options of the tensorial structure
are:

$\displaystyle\boldsymbol{u}$
$\displaystyle=\boldsymbol{u}^{n}_{\\{1,2,3\\}},$ (33)
$\displaystyle\boldsymbol{u}$
$\displaystyle=\boldsymbol{u}^{n_{1}}_{\\{1,2\\}}\times\boldsymbol{u}^{n_{2}}_{\\{3\\}},$
$\displaystyle\boldsymbol{u}$
$\displaystyle=\boldsymbol{u}^{n_{1}}_{\\{1,3\\}}\times\boldsymbol{u}^{n_{2}}_{\\{2\\}},$
$\displaystyle\boldsymbol{u}$
$\displaystyle=\boldsymbol{u}^{n_{1}}_{\\{2,3\\}}\times\boldsymbol{u}^{n_{2}}_{\\{1\\}},$
$\displaystyle\boldsymbol{u}$
$\displaystyle=\boldsymbol{u}^{n_{1}}_{\\{1\\}}\times\boldsymbol{u}^{n_{2}}_{\\{2\\}}\times\boldsymbol{u}^{n_{3}}_{\\{3\\}}.$

(a)

(b)

(c)

Figure 3: Visualizations of quadrature points with different tensorial
structure options

Visualizations of quadrature points with different tensorial structures can be
found in Fig. 3.

Given a specific tensorial structure option for the quadrature rule, we need
to decide on the method for generating the quadrature points within each space
and determine the requisite number of quadrature points to attain the desired
level of accuracy. In our method, for a particular tensorial structure of the
quadrature rule, we employ the designed quadrature method in each space. This
method determines the nodes and weights to accurately integrate polynomials in
that space up to a specified order. Adapted from the Gauss quadrature rule, we
use the level of accuracy $k$ for a quadrature rule to indicate its capability
of precisely integrating polynomials up to the $(2k-1)$th order in the
integration space. The number of quadrature points utilized in the designed
quadrature method is given by (29), tailored for high-dimensional integration
problems. In our setting, we also utilize the designed quadrature method for
relatively low-dimensional spaces. Therefore, we choose the required number of
quadrature points using:

$n=\left\\{\begin{array}[]{ll}\frac{k^{d}}{d}&\text{ for }d\leq 2\\\
0.9\frac{(2d)^{k-1}}{(k-1)!}&\text{ for }d>2.\end{array}\right.$ (34)

When the number of dimensions is greater than two, we follow the designed
quadrature rule in (29) with $\eta=0.9$. However, for 1D and 2D cases, we
choose $n=\frac{k^{d}}{d}$. In this way, in the 1D space, it chooses $k$
quadrature points to achieve $k$ level of accuracy, which matches the Gauss
quadrature rule. Consequently, for the 1D space, we can avoid solving the
optimization problem in the designed quadrature method and directly use the 1D
Gauss quadrature rule, as it represents the exact solution of the designed
quadrature method. This approach results in generating the full-grid Gauss
quadrature rule if a fully tensorial structure option is chosen, and
generating the designed quadrature rule if a non-tensorial structure option is
chosen.

###  3.2 Tensor-grid evaluations with AMTC

When not using the AMTC method for model evaluations, tensor-structured
quadrature rules often underperform compared to non-tensorial quadrature
rules. This is because tensor-structured quadrature rules always require
evaluations of the model at more input points than non-tensorial quadrature
rules to achieve the same level of accuracy. However, when utilizing the
computational graph transformation method, AMTC, redundant operations incurred
by the tensorial structure of the input points are eliminated at the operation
level. Consequently, the overall cost of evaluating the model is no longer
directly tied to the number of quadrature points utilized. This gives us the
reason for devising a custom tensor-structured quadrature rule tailored for a
given computational model.

For instance, consider a 3D UQ problem involving two computationally expensive
solvers. The computational model can be written as

$\displaystyle f_{1}=\mathcal{F}_{1}(u_{1},u_{2}),$ (35) $\displaystyle
f_{2}=\mathcal{F}_{2}(f_{1},u_{3}),$

where $\mathcal{F}_{1}$ and $\mathcal{F}_{2}$ represent the first and second
solver, respectively, and the UQ problem aims to compute the risk measures of
$f_{2}$ under the effect of the three uncertain inputs. In this case, we have
two uncertain inputs associated with solver 1, and one uncertain input
associated with solver 2. Consequently, the operations in solver 1 are
dependent on $u_{1}$ and $u_{2}$ while the operations in solver 2 are
dependent on $u_{1}$, $u_{2}$, and $u_{3}$. In such scenario, an optimal
quadrature rule can be easily formulated following a partially tensorial
structure:

$\boldsymbol{u}=\boldsymbol{u}^{k_{1}}_{\\{1,2\\}}\times\boldsymbol{u}^{k_{2}}_{\\{3\\}}.$
(36)

The corresponding model evaluations can be significantly accelerated using the
AMTC method. To demonstrate, we show the computational graphs when performing
model evaluations on fully tensorial quadrature points
$\boldsymbol{u}^{k}_{\\{1\\}}\times\boldsymbol{u}^{k}_{\\{2\\}}\times\boldsymbol{u}^{k}_{\\{3\\}}$
and partially tensorial quadrature points
$\boldsymbol{u}^{k_{1}}_{\\{1,2\\}}\times\boldsymbol{u}^{k_{2}}_{\\{3\\}}$
after applying AMTC in Fig. 4. We observe that, in both cases, AMTC generates
modified computational graphs that reduce the number of evaluations on solver
1. In comparison, the partially tensor-structured quadrature rule is more
favorable than the fully tensor-structured quadrature rule, as both the
$u_{1}$ and $u_{2}$ affect the same operations on the computational graph, so
there is no need to maintain a tensor structure between $u_{1}$ and $u_{2}$
like in the fully tensor-structured quadrature rule. For the partially tensor-
structured quadrature rule, by grouping the space of $u_{1}$ and $u_{2}$, the
total number of quadrature points in this space, $k_{1}$, is smaller than the
full-grid quadrature rule which has $k^{2}$ quadrature points in this space.
Consequently, utilizing the partially tensor-structured quadrature rule
results in fewer evaluations for both solver 1 and solver 2.

(a)

(b)

Figure 4: Computational graphs after AMTC on fully and partially tensor
structured input points

###  3.3 Finding a desirable tensor structure option

For specific UQ problems, particularly those involving multidisciplinary
systems, certain partially tensor-structured quadrature rules demonstrate
superior performance compared to existing ones when utilized with the graph-
accelerated NIPC method. However, identifying an optimal tensor-structure
option may be challenging.

Theoretically, we can estimate the model evaluation costs after applying AMTC
for each tensor-structured rule. To do that, for a computational model
$f=\mathcal{F}(u)$, we denote the data and operations nodes as $\xi_{i}$ and
$\varphi_{j}$, respectively. Here $l$ represents the number of operations in
the computational graph. For example, consider the function

$f=cos(u_{1})+exp(-u_{2}),$ (37)

with its corresponding computational graph described as

$\displaystyle\xi_{1}=\varphi_{1}(u_{1})=cos(u_{1});$ (38)
$\displaystyle\xi_{2}=\varphi_{2}(u_{2})=-u_{2};$
$\displaystyle\xi_{3}=\varphi_{3}(\xi_{1})=exp(\xi_{2});$ $\displaystyle
f=\varphi_{4}(\xi_{2},\xi_{3})=\xi_{2}+\xi_{3}.$

From the computational graph, we can generate the dependency information which
stores whether the output of one operation is dependent on one input as a
binary number, written as $D(\varphi_{i},u_{j})$. The dependency information
can be viewed in a matrix. The dependency matrix for the function (37) is
presented in Tab. 2, which indicates the dependencies between operations and
input variables.

Table 2: Operations dependency information for function $f=cos(u_{1})+exp(-u_{2})$ $D(\varphi_{i},u_{j})$ | $u_{1}$ | $u_{2}$
---|---|---
$\varphi_{1}$ | 1 | 0
$\varphi_{2}$ | 0 | 1
$\varphi_{3}$ | 0 | 1
$\varphi_{4}$ | 1 | 1

To estimate the evaluation cost of each tensorial structure for the quadrature
rule to achieve the $k$-level of accuracy, we also need to store the
information regarding the evaluation cost of each operation on one quadrature
point. This information can be represented either by the total number of
floating-point operations (FLOPs) or the evaluation time. For the non-
tensorial structure, if we assume all operations are dependent on at least one
input, no repeated evaluations can be reduced by the AMTC method. Assume that
the operations’ cost on all of the quadrature points is roughly the same, the
total model evaluation cost can be approximated as:

$\text{Cost}\approx nO\left(f(u)\right)\approx n\sum_{i=1}^{l}O(\varphi_{i}),$
(39)

where $O(\cdot)$ represents the evaluation cost of the function/operation and
the number of quadrature points $n$ is chosen as in (34). If we maintain a
full-grid structure with $k$ quadrature points in each dimension, the
evaluation cost without the AMTC method is

$\text{Cost}\approx k^{d}O\left(\mathcal{F}(u^{(1)})\right)\approx
k^{d}\sum_{i=1}^{l}O(\varphi_{i}).$ (40)

However, with the AMTC method, each operation in the graph is evaluated only
for distinct values. In the full-grid case, if the output of an operation is
dependent on $\tilde{d}$ number of uncertain inputs, there are only
$k^{\tilde{d}}$ quadrature points in its dependent inputs space, and that
operation will be evaluated for $k^{\tilde{d}}$ number of times using the AMTC
method. In this case, we omit the cost of the Einstein summation nodes added
to the computational graph by the AMTC method, the evaluation cost can be
estimated as

$\text{Cost}\approx\sum_{i=1}^{l}k^{\sum_{j=1}^{d}D(\varphi_{i},u_{j})}O(\varphi_{i}).$
(41)

Using the same idea, we can also estimate the model evaluations cost for any
tensorial structure option of the quadrature points. For example, consider the
setting where we have a partially tensorial structure for the quadrature
points, such as

$\boldsymbol{u}=\boldsymbol{u}^{k_{1}}_{\\{1,2\\}}\times\boldsymbol{u}^{k_{2}}_{\\{3\\}}.$
(42)

In our method, $k_{1}$ and $k_{2}$ are chosen from (34) using $d=1$ and $d=2$
respectively. In this case, there are $k_{1}$ quadrature points in the space
of $u_{1}$ and $u_{2}$, and the $k_{2}$ quadrature points in the space of
$u_{3}$. With the AMTC method, for the operations that are dependent on only
$u_{1}$, only $u_{2}$ or only $u_{1}$ and $u_{2}$, they will be evaluated
$n_{1}$ times. Operations dependent on only $u_{3}$, will be evaluated $n_{2}$
times. Operations not dependent on any input are evaluated once. The rest of
the operations are evaluated $n_{1}n_{2}$ times. The model evaluations cost
for this tensor-structured quadrature rule can be estimated as

$\text{Cost}\approx\sum_{i=1}^{l}k_{1}^{D(\varphi_{i},u_{1})D(\varphi_{i},u_{2})}k_{2}^{D(\varphi_{i},u_{3})}O(\varphi_{i}).$
(43)

Estimating the total model evaluation cost for all tensor structure options
enables us to choose the optimal one with the lowest estimated cost. However,
the challenge lies in the fact that the total number of tensor structure
options follows the Bell number as shown in (32), which scales significantly
with the problem dimensions. To avoid this significant overhead cost
associated with numerically finding the optimal tensor structure option, we
propose an alternative approach. This approach involves selecting a desirable
tensor structure option based on our understanding of the computational model
and an analysis of the computational graph. By leveraging information from the
computational graph and characteristics of the computational model, we can
make informed decisions to choose a tensor structure that is likely to lead to
efficient model evaluations. This approach balances computational efficiency
with accuracy, making it a practical solution for many UQ problems with
complicated computational models.

When dealing with computational models involving multiple disciplines or
solvers, identifying a desirable tensor-structure option can sometimes be
achieved by examining the computational graph at the solver level. For
instance, in the example provided in (35), we can easily deduce the optimal
tensor structure based on insights gained from the computational graph. In
other cases, we can perform some analysis on the computational graph to help
us decide on the desired tensor-structure option. Our rationale here is based
on the observation that complex computational models typically contain one or
two computationally expensive operations, such as linear or non-linear
solvers, which significantly contribute to the overall computational cost. We
aim to identify uncertain inputs that these expensive operations do not depend
on and then establish a tensor structure to minimize the number of evaluations
required for these operations, with the utilization of AMTC. To do that, we
define the sparsity ratio (SR) of an uncertain input as the ratio of its
dependent operations’ cost over the total model evaluation cost which can be
approximated as

$\text{SR}(u_{i})\approx\frac{\sum_{i=1}^{l}O(\varphi_{i})D(\varphi_{i},u_{j})}{O(f(u))},$
(44)

where $\text{SR}(u_{i})$ represents the sparsity ratio of $u_{i}$ and it
measures the percentage of the model evaluation cost that is dependent on
$u_{i}$. We identify the sparse uncertain inputs based on the condition that
their sparsity ratio is $<5\%$. We then partition the uncertain inputs as

$u=(u_{s},u_{ns}),$ (45)

where $u_{s}=(u_{s_{1}},\ldots,u_{s_{d_{1}}})\in\mathbb{R}^{d_{1}}$ is the set
of sparse uncertain inputs and
$u_{ns}=(u_{ns_{1}},\ldots,u_{ns_{d_{2}}})\in\mathbb{R}^{d_{2}}$ is the set of
non-sparse uncertain inputs with $d=d_{1}+d_{2}$.

In this case, if all of the uncertain inputs are sparse uncertain inputs, we
may want to form a full-grid quadrature rule following the Gauss quadrature
method, such that the quadrature points can be written as

$\boldsymbol{u}=\boldsymbol{u}^{k}_{\\{1\\}}\times\ldots\times\boldsymbol{u}^{k}_{\\{d\\}},$
(46)

where $\boldsymbol{u}^{k}_{\\{i\\}}$ represents the $k$ quadrature points in
$u_{i}$ space. In this way, even though the number of quadrature points
increases exponentially as $n=k^{d}$, the model evaluation cost can be
significantly accelerated using the AMTC method which only evaluates each
operation on the distinct quadrature points in the space of the uncertain
inputs on which it depends. However, in most cases, for most UQ problems,
there often only exist a few sparse uncertain inputs. In this case, we suggest
choosing the quadrature rule that maintains a tensor structure between the
quadrature points in the sparse uncertain input space and quadrature points in
the space of non-sparse uncertain inputs, such that the quadrature points can
be written as

$\boldsymbol{u}=\boldsymbol{u}^{k_{1}}_{\\{ns\\}}\times\boldsymbol{u}^{k_{2}}_{\\{s\\}},$
(47)

where $\boldsymbol{u}^{k_{1}}_{\\{ns\\}}$ represents quadrature points in
$u_{ns}$ space with $k_{1}$ points and $\boldsymbol{u}^{k_{2}}_{\\{s\\}}$
represent the quadrature points in $u_{{s}}$ space with $k_{2}$ points. This
results in a total of $k_{1}k_{2}$ quadrature points in the quadrature rule.
However, when using the AMTC method to modify the computational graph, the
computational cost on the modified computational graph would scale roughly
linearly with $k_{1}$ as most of the expensive operations in the computational
graph would only depend on $u_{ns}$ and would only be evaluated $n_{1}$ times.
An alternative option is to formulate the quadrature rule as

$\boldsymbol{u}=\boldsymbol{u}^{k_{1}}_{\\{ns\\}}\times\boldsymbol{u}^{k}_{\\{s_{1}\\}}\times\ldots,\boldsymbol{u}^{k}_{\\{s_{d_{1}}\\}},$
(48)

which forms a tensor structure between the space of each sparse uncertain
input. This option results in more quadrature points but it is easier to
generate the quadrature rule and can be more effective to use if we know the
sparse uncertain inputs are associated with different solvers.

We show the specific steps of using this method with graph-accelerated NIPC in
Alg. 1.

Algorithm 1 Partially tensor-structured quadrature rule for graph-accelerated
NIPC

1:Specify a computational graph for a numerical model and the uncertain
inputs, $U$

2:Specify the level of accuracy of the quadrature rule, $k$

3:Choosing a desired tensor structure option based on the analysis of the
computational graph

4:Compute the required number of quadrature points in each space following
(29)

5:Compute the quadrature points and weights in each space by solving (24)

6:Formulate the quadrature rule using the chosen tensorial structure

7:Perform the model evaluations on the tensor-structured quadrature rule using
AMTC

8:Post-process the evaluations using NIPC to solve the UQ problem

##  4 Numerical Results

###  4.1 Laser-beam powered UAV multidisciplinary model

Figure 5: A circular cruise mission around a ground station Uncertain inputs | Distributions
---|---
Payload weight $W_{p}$ ($kg$) | $\mathcal{N}(90,10)$
Atmospheric extinction coefficient $\sigma$ | $\mathcal{N}(0.2,0.02)$
Flight altitude $h$ (m) | $\mathcal{N}(10,000,1000)$
Flight velocity $v$ ($m/s$) | $\mathcal{N}(100,10)$
Table 3: Uncertain inputs’ distributions for the UAV problem Figure 6:
Multidisciplinary structure of the UAV model

The first test problem is a 4D UQ problem adapted from a UAV design scenario
[30]. This computational model assesses the total energy stored during a
cruise mission undertaken by a laser-beam-powered UAV as it orbits a ground
station for one complete cycle while being recharged by the laser beam emitted
from the ground station. The mission is demonstrated in Fig. 5. We consider
four uncertain inputs associated with the operating conditions, which are
flight velocity, flight altitude, payload weight, and atmospheric extinction
coefficient. Their distributions are shown in Tab. 3. The objective of this UQ
problem is to compute the expectation of the stored energy considering the
influence of these uncertain inputs.

The computational model used in this problem involves four disciplines:
aerodynamics, laser beams, weights, and performance. The multidisciplinary
structure of the computational model is shown in Fig. 6. For the laser-beam
model, we apply the Beers–Lambert law for optical power transfer [31] to
compute the power received by the aircraft as

$P_{\text{rec}}=\eta P_{\text{tra}},\quad\eta=e^{-\sigma R},$ (49)

where $P_{\text{rec}}$ and $P_{\text{tra}}$ denote the power received and
transmitted, respectively; $\sigma$ is the atmospheric extinction coefficient
which differs across weather conditions; and $R$ is the distance between the
aircraft and the power source. In the aerodynamic model, inputs such as flight
velocity, flight altitude, and the total aircraft weight are considered. This
model is responsible for computing the total drag force exerted on the wing to
sustain steady-state flight. Additionally, it determines the lift coefficient
and drag coefficient by employing a linear model and a quadratic polar model,
respectively, derived from airfoil data. The weights model estimates the
weight of each aircraft component using the statistical equations from [32]
and outputs the total aircraft weight. The performance model takes inputs such
as flight velocity and total drag force to compute the total energy stored
during the mission. This computation is based on the following equations:

$\displaystyle n=\sqrt{\left(\frac{v}{rg}\right)^{2}+1},$ (50) $\displaystyle
P_{\text{req}}=\frac{D}{n},$ $\displaystyle
E=(P_{\text{rec}}-P_{\text{req}})t,$

where $n$ represents the turn load, $r$ is the radius of the mission, $g$ is
the gravitational acceleration, $D$ denotes the total drag force,
$P_{\text{req}}$ is the required power to maintain steady flight, and $t$
represents the mission time.

On this test problem, following our proposed method, sparsity ratios of the
uncertain inputs are measured, which are shown in Tab. 4.

Table 4: Sparsity ratio of the uncertain inputs in the UAV model Uncertain inputs | Sparsity ratio
---|---
$v$ | 92%
$h$ | 96%
$W_{p}$ | 4%
$\sigma$ | 2%

Using the sparsity ratios, two sparse uncertain inputs are identified: $W_{p}$
and $\sigma$, as both of them have a sparsity ratio less than 5%. Combining
this information with our understanding of the computational model, we
constructed the partially tensor-structured quadrature rule as:

$\boldsymbol{u}=\boldsymbol{u}^{n_{1}}_{\\{h,v\\}}\times\boldsymbol{u}^{n_{2}}_{\\{W_{p}\\}}\times\boldsymbol{u}^{n_{3}}_{\\{\sigma\\}}.$
(51)

Our proposed method is compared with the designed quadrature method and the
full-grid Gauss quadrature method, all of which are used with the NIPC method
to solve the UQ problem. The UQ convergence plots, with and without using
AMTC, are shown in Fig. 7. From the results, we observe that the partially
tensor-structured quadrature rule is outperformed by the designed quadrature
method when AMTC is not used to accelerate the model evaluations. This finding
aligns with our theoretical analysis, as the non-tensorial quadrature rule
requires fewer quadrature points than the tensor-structured quadrature rules
to achieve the same level of accuracy. However, when we utilize AMTC to
accelerate the model evaluations, the advantages of the structured quadrature
rules become evident. Both the model evaluations associated with the full-grid
quadrature rule and the partially tensor-structured quadrature rule are
significantly accelerated by using the AMTC method. In this scenario, the
partially tensor-structured quadrature rule achieves the best performance,
with more than 50% reduction in evaluation costs than the other two methods.
This is because, our method constructed a tailored tensorial structure that
makes the most effective use of the inherent sparsity within the computational
graph of this multidisciplinary system, resulting in the lowest model
evaluation costs under the same level of accuracy.

(a)

(b)

Figure 7: Convergence plots with and without AMTC for the UAV problem

###  4.2 Air-taxi trajectory design multidisciplinary model

The second test problem is a 6D UQ problem adapted from an air-taxi trajectory
optimization scenario. The computational models compute the average sound
pressure level (SPL) at the ground when the air taxi flies along a given
trajectory. For the UQ problem, we aim to compute the expectation of the
average sound pressure level under the effect of six uncertain inputs. Three
uncertain inputs are associated with the initial conditions of the trajectory,
and the other three uncertain inputs are associated with the parameters in the
acoustic model. The distributions of the uncertain inputs can be found in Tab.
5. Notably, the control inputs used in this problem are generated from solving
a deterministic trajectory optimization problem, which aims to minimize the
total energy expended throughout the trajectory. Details of the computational
model along with the trajectory optimization formulation can be found in [33].

Uncertain inputs | Distributions
---|---
Initial velocity $v_{0}$ ($m/s$) | $\mathcal{U}(30,35)$
Initial flight path angel $\gamma_{0}$ ($\deg$) | $\mathcal{U}(0,2)$
Initial altitude $h_{0}$ ($m$) | $\mathcal{U}(0,20)$
$\beta_{1}$ | $\mathcal{N}(0.0209,0.001)$
$\beta_{2}$ | $\mathcal{N}(18.2429,0.5)$
$\beta_{3}$ | $\mathcal{N}(6.729,0.2)$
Table 5: Uncertain inputs’ distributions for the air-taxi problem

In this problem, the measured sparsity ratios are shown in Tab. 6. Based on
this information, three uncertain inputs are identified as sparse uncertain
inputs which are the three uncertain inputs associated with the acoustic
model. This finding aligns with the multidisciplinary structure of the
computational model, depicted in Fig. 8. From this figure, we observe that the
computational model comprises four disciplines: flight dynamics, aerodynamics,
propulsion, and acoustics, with two-way coupling among the first three
disciplines. In this case, the coupled disciplines functions as a nonlinear
operation and is solved iteratively. Due to its dominant role in computational
costs, the three uncertain inputs associated with the nonlinear solver are the
dense uncertain inputs in this computational model. Conversely, the parameter
uncertain inputs are the sparse uncertain inputs as they do not affect the
upstream nonlinear solver in this computational model. Based on this
reasoning, we choose the partially tensor-structured option as

$\boldsymbol{u}=\boldsymbol{u}^{n_{1}}_{\\{v_{0},h_{0},\gamma_{0}\\}}\times\boldsymbol{u}^{n_{2}}_{\\{\beta_{1},\beta_{2},\beta_{3}\\}},$
(52)

in which we form a tensor structure between the space of dense uncertain
inputs and the space of sparse uncertain inputs.

Figure 8: Multidisciplinary structure of the air-taxi model Table 6: Sparsity ratio of the uncertain inputs in the air-taxi model Uncertain inputs | Sparsity ratio
---|---
$v_{0}$ | 99%
$h_{0}$ | 99%
$\gamma_{0}$ | 99%
$\beta_{1}$ | 2%
$\beta_{2}$ | 2%
$\beta_{3}$ | 2%

The UQ convergence plots with and without using AMTC are shown in Fig. 9.
Similar to the first test problem, the partially tensor-structured quadrature
rule is outperformed by the designed quadrature method when AMTC is not used.
However, it becomes the most efficient method when using AMTC to accelerate
the model evaluations. The reduction in computational costs is more than 40%
in most cases. This can be explained by the fact that the tailored tensorial
structure option takes full advantage of the inherent sparsity of the
computational graph, resulting in the most efficient performance with the AMTC
method. However, we also observe that when the desired relative error is
exceptionally low (1e-6 and lower), the performance of the partially tensor-
grid quadrature rule deteriorates on this test problem. This phenomenon can be
attributed to the fact that the designed quadrature rule method can achieve
similar performance to the Gauss quadrature rule while requiring a
significantly smaller number of model evaluations. However, as the dimension
of the optimization problem increases (as we have more optimization variables
and more constraints), it becomes more challenging to solve the optimization
problem to a tight tolerance using this method. Consequently, the designed
quadrature method may struggle to achieve the same level of accuracy as the
Gauss quadrature rule when an extremely high level of accuracy is required.

(a)

(b)

Figure 9: Convergence plots with and without AMTC for the airtaxi problem

##  5 Conclusion

In this paper, we introduced a novel approach for generating partially tensor-
structured quadrature rules tailored for utilization with the graph-
accelerated NIPC method to tackle UQ problems. Our method entails initially
choosing an appropriate tensorial structure option for the quadrature rule
through an analysis of the computational graph. Subsequently, the quadrature
rule possessing the chosen tensorial structure is generated using the designed
quadrature method. This method has been tested on two UQ problems: one with
four uncertain inputs and the other with six uncertain inputs, both involving
multidisciplinary systems and are derived from aircraft design scenarios. The
numerical results show that our method generates a new, partially tensor-
structured quadrature rule for both problems. Notably, this new rule exhibits
superior performance compared to both the existing designed quadrature method
and the full-grid Gauss quadrature method when used with the graph-accelerated
NIPC method.

One limitation of this proposed method is its reliance on the inherent
sparsity of the computational graph, which restricts the generation of new
quadrature rules for every problem. However, our method broadens the scope of
applicability of the graph-accelerated NIPC method, particularly to medium-
dimensional UQ problems featuring multidisciplinary systems. This is achieved
by generating partially tensor-structured rules to more efficiently leverage
the inherent sparsity within the computational graph. Another limitation lies
in the potential performance decline of our method for high-dimensional UQ
problems (ten dimensions or higher), where alternative approaches such as
Monte Carlo methods may outperform it. A promising direction for future
research involves integrating this method with dimension-reduction techniques
to address higher-dimensional UQ problems.

## Acknowledgments

The material presented in this paper is, in part, based upon work supported by
DARPA under grant No. D23AP00028-00.

## References

  * Joslyn and Savelli [2010] Joslyn, S., and Savelli, S., “Communicating forecast uncertainty: Public perception of weather forecast uncertainty,” _Meteorological Applications_ , Vol. 17, No. 2, 2010, pp. 180–195.  10.1002/met.190.
  * Pappenberger et al. [2005] Pappenberger, F., Beven, K. J., Hunter, N. M., Bates, P. D., Gouweleeuw, B. T., Thielen, J., and de Roo, A. P. J., “Cascading model uncertainty from medium range weather forecasts (10 days) through a rainfall-runoff model to flood inundation predictions within the European Flood Forecasting System (EFFS),” _Hydrology and Earth System Sciences_ , Vol. 9, No. 4, 2005, pp. 381–393.  10.5194/hess-9-381-2005.
  * Hüllermeier and Waegeman [2021] Hüllermeier, E., and Waegeman, W., “Aleatoric and epistemic uncertainty in machine learning: An introduction to concepts and methods,” _Machine Learning_ , Vol. 110, 2021, pp. 457–506.  10.1007/s10994-021-05946-3.
  * Abdar et al. [2021] Abdar, M., Pourpanah, F., Hussain, S., Rezazadegan, D., Liu, L., Ghavamzadeh, M., Fieguth, P., Cao, X., Khosravi, A., Acharya, U. R., et al., “A review of uncertainty quantification in deep learning: Techniques, applications and challenges,” _Information Fusion_ , Vol. 76, 2021, pp. 243–297.
  * Wan et al. [2014] Wan, H.-P., Mao, Z., Todd, M. D., and Ren, W.-X., “Analytical uncertainty quantification for modal frequencies with structural parameter uncertainty using a Gaussian process metamodel,” _Engineering Structures_ , Vol. 75, 2014, pp. 577–589.  10.1016/j.engstruct.2014.06.028.
  * Hu et al. [2018] Hu, Z., Mahadevan, S., and Ao, D., “Uncertainty aggregation and reduction in structure–material performance prediction,” _Computational Mechanics_ , Vol. 61, No. 1, 2018, pp. 237–257.  10.1007/s00466-017-1448-6.
  * Ng and Willcox [2016] Ng, L. W., and Willcox, K. E., “Monte Carlo information-reuse approach to aircraft conceptual design optimization under uncertainty,” _Journal of Aircraft_ , Vol. 53, No. 2, 2016, pp. 427–438.  10.2514/1.C033352.
  * Wang et al. [2024a] Wang, B., Orndorff, N. C., Joshy, A. J., and Hwang, J. T., “Graph-accelerated large-scale multidisciplinary design optimization under uncertainty of a laser-beam-powered aircraft,” _AIAA SCITECH 2024 Forum_ , 2024a, p. 0169.  10.2514/6.2024-0169.
  * Lim et al. [2022] Lim, D., Kim, H., and Yee, K., “Uncertainty propagation in flight performance of multirotor with parametric and model uncertainties,” _Aerospace Science and Technology_ , Vol. 122, 2022, p. 107398.  10.1016/j.ast.2022.107398.
  * Xiu and Karniadakis [2002] Xiu, D., and Karniadakis, G. E., “The Wiener–Askey polynomial chaos for stochastic differential equations,” _SIAM journal on scientific computing_ , Vol. 24, No. 2, 2002, pp. 619–644.  10.1137/S1064827501387826.
  * Jones et al. [2013] Jones, B. A., Doostan, A., and Born, G. H., “Nonlinear propagation of orbit uncertainty using non-intrusive polynomial chaos,” _Journal of Guidance, Control, and Dynamics_ , Vol. 36, No. 2, 2013, pp. 430–444.  10.2514/1.57599.
  * Hosder et al. [2007] Hosder, S., Walters, R., and Balch, M., “Efficient sampling for non-intrusive polynomial chaos applications with multiple uncertain input variables,” _48th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference_ , 2007, p. 1939.  10.2514/6.2007-1939.
  * Blatman and Sudret [2011] Blatman, G., and Sudret, B., “Adaptive sparse polynomial chaos expansion based on least angle regression,” _Journal of computational Physics_ , Vol. 230, No. 6, 2011, pp. 2345–2367.  10.1016/j.jcp.2010.12.021.
  * Keshavarzzadeh et al. [2017] Keshavarzzadeh, V., Fernandez, F., and Tortorelli, D. A., “Topology optimization under uncertainty via non-intrusive polynomial chaos expansion,” _Computer Methods in Applied Mechanics and Engineering_ , Vol. 318, 2017, pp. 120–147.  10.1016/j.cma.2017.01.019.
  * Xiaojing et al. [2018] Xiaojing, W., Zhang, W., Shufang, S., and Zhengyin, Y., “Sparse grid-based polynomial chaos expansion for aerodynamics of an airfoil with uncertainties,” _Chinese Journal of Aeronautics_ , Vol. 31, No. 5, 2018, pp. 997–1011.  10.1016/j.cja.2018.03.011.
  * Eldred and Burkardt [2009] Eldred, M., and Burkardt, J., “Comparison of non-intrusive polynomial chaos and stochastic collocation methods for uncertainty quantification,” _47th AIAA aerospace sciences meeting including the new horizons forum and aerospace exposition_ , 2009, p. 976.  10.2514/6.2009-976.
  * Golub and Welsch [1969] Golub, G. H., and Welsch, J. H., “Calculation of Gauss quadrature rules,” _Mathematics of computation_ , Vol. 23, No. 106, 1969, pp. 221–230.  10.1090/S0025-5718-69-99647-1.
  * Smolyak [1963] Smolyak, S. A., “Quadrature and interpolation formulas for tensor products of certain classes of functions,” _Doklady Akademii Nauk_ , Vol. 148, Russian Academy of Sciences, 1963, pp. 1042–1045.
  * Keshavarzzadeh et al. [2018] Keshavarzzadeh, V., Kirby, R. M., and Narayan, A., “Numerical integration in multiple dimensions with designed quadrature,” _SIAM Journal on Scientific Computing_ , Vol. 40, No. 4, 2018, pp. A2033–A2061.  10.1137/17M1137875.
  * Wang et al. [2024b] Wang, B., Sperry, M., Gandarillas, V. E., and Hwang, J. T., “Accelerating model evaluations in uncertainty propagation on tensor grids using computational graph transformations,” _Aerospace Science and Technology_ , Vol. 145, 2024b, p. 108843.  10.1016/j.ast.2023.108843.
  * Wiener [1938] Wiener, N., “The homogeneous chaos,” _American Journal of Mathematics_ , Vol. 60, No. 4, 1938, pp. 897–936.  10.2307/2371268.
  * Ghanem [1999] Ghanem, R., “Ingredients for a general purpose stochastic finite elements implementation,” _Computer methods in applied mechanics and engineering_ , Vol. 168, No. 1-4, 1999, pp. 19–34.  10.1016/S0045-7825(98)00106-6.
  * Lee and Rahman [2020] Lee, D., and Rahman, S., “Practical uncertainty quantification analysis involving statistically dependent random variables,” _Applied Mathematical Modelling_ , Vol. 84, 2020, pp. 324–356.  10.1016/j.apm.2020.03.041.
  * Lüthen et al. [2021] Lüthen, N., Marelli, S., and Sudret, B., “Sparse polynomial chaos expansions: Literature survey and benchmark,” _SIAM/ASA Journal on Uncertainty Quantification_ , Vol. 9, No. 2, 2021, pp. 593–649.  10.1137/20M1315774.
  * Wang et al. [2022] Wang, B., Sperry, M., Gandarillas, V. E., and Hwang, J. T., “Efficient uncertainty propagation through computational graph modification and automatic code generation,” _AIAA AVIATION 2022 Forum_ , 2022, p. 3997.  10.2514/6.2022-3997.
  * Gandarillas et al. [(accepted] Gandarillas, V., Joshy, A. J., Sperry, M. Z., Ivanov, A. K., and Hwang, J. T., “A graph-based methodology for constructing computational models that automates adjoint-based sensitivity analysis,” _Structural and Multidisciplinary Optimization_ , (accepted).
  * Sarojini et al. [2023] Sarojini, D., Ruh, M. L., Joshy, A. J., Yan, J., Ivanov, A. K., Scotzniovsky, L., Fletcher, A. H., Orndorff, N. C., Sperry, M., Gandarillas, V. E., et al., “Large-scale multidisciplinary design optimization of an evtol aircraft using comprehensive analysis,” _AIAA SciTech 2023 Forum_ , 2023, p. 0146.  10.2514/6.2023-0146.
  * Ruh et al. [2024] Ruh, M. L., Fletcher, A., Sarojini, D., Sperry, M., Yan, J., Scotzniovsky, L., van Schie, S. P., Warner, M., Orndorff, N. C., Xiang, R., et al., “Large-scale multidisciplinary design optimization of a NASA air taxi concept using a comprehensive physics-based system model,” _AIAA SCITECH 2024 Forum_ , 2024, p. 0771.  10.2514/6.2024-0771.
  * Gandarillas and Hwang [2023] Gandarillas, V., and Hwang, J. T., “TALOS: A toolbox for spacecraft conceptual design,” _arXiv preprint arXiv:2303.14936_ , 2023.  10.48550/arXiv.2303.14936.
  * Orndorff et al. [2023a] Orndorff, N. C., Wang, B., Ruh, M. L., Fletcher, A., and Hwang, J. T., “Gradient-based sizing optimization of power-beaming-enabled aircraft,” _AIAA AVIATION 2023 Forum_ , 2023a, p. 4019.  10.2514/6.2023-4019.
  * Kim et al. [2001] Kim, I. I., McArthur, B., and Korevaar, E. J., “Comparison of laser beam propagation at 785 nm and 1550 nm in fog and haze for optical wireless communications,” _Optical wireless communications III_ , Vol. 4214, Spie, 2001, pp. 26–37.  10.1117/12.417512.
  * Raymer [2012] Raymer, D., _Aircraft design: a conceptual approach_ , American Institute of Aeronautics and Astronautics, Inc., 2012.  10.2514/4.869112.
  * Orndorff et al. [2023b] Orndorff, N. C., Sarojini, D., Scotzniovsky, L., Gill, H., Lee, S., Cheng, Z., Zhao, S., Mi, C., and Hwang, J. T., “Air-taxi transition trajectory optimization with physics-based models,” _AIAA SCITECH 2023 Forum_ , 2023b, p. 0324.  10.2514/6.2023-0324.


# Real-Time Electric Vehicle Smart Charging at Workplaces: A Real-World Case
Study

Nathaniel Tucker2, Gustavo Cezar3, and Mahnoosh Alizadeh2  2Department of
Electrical and Computer Engineering, University of California, Santa Barbara,
California, 93106, USA  3SLAC National Accelerator Laboratory, GISMo Group,
California, 94025, USA

###### Abstract

We study a real-time smart charging algorithm for electric vehicles (EVs) at a
workplace parking lot in order to minimize electricity cost from time-of-use
electricity rates and demand charges while ensuring that the owners of the EVs
receive adequate levels of charge. Notably, due to real-world constraints, our
algorithm is agnostic to both the state-of-charge and the departure time of
the EVs and uses scenario generation to account for each EV’s unknown future
departure time as well as certainty equivalent control to account for the
unknown EV arrivals in the future. Real-world charging data from a Google
campus in California allows us to build realistic models of charging demand
for each day of the week. We then compare various results from our smart
charging algorithm to the status quo for a two week period at a Google parking
location.

##  I Introduction

Many large companies across the United States are installing Electric Vehicle
Supply Equipment (EVSEs) within their parking lots to allow for employees and
visitors to charge their Electric Vehicles (EVs) during their stay [1]. There
are numerous benefits that can come from workplace charging stemming from the
highly flexible nature of EV charging at workplaces [2, 3] (i.e., EVs are
plugged in for long durations and their charging can be shifted depending on
other conditions including behind-the-meter solar generation or electricity
rates). For example, Google (who has installed EVSEs at several of its Bay
Area locations for employee EV charging) allows employees to leave their EVs
connected to the EVSE for long durations and does not force them to unplug the
EV after reaching full charge, whereas other workplaces might require
employees to move EVs after charging is complete. The approach of not
enforcing EVs to be removed once fully charged may require more investment in
EVSEs to satisfy the needs of the employees, but the investment would allow
for greater flexibility that can be harnessed by a smart charging algorithm to
benefit the workplace (by providing financial savings) and the utility (by
enabling its customers to respond to price signals) [4, 5].

The Grid Integration, Systems, and Mobility (GISMo) group within the SLAC
National Accelerator Laboratory and the Smart Infrastrucuture Systems Lab
(SISL) at UCSB are working on a project to study the impacts of electric
vehicles on the distribution system and designing tools and algorithms to be
used by workplace charging station operators to efficiently manage asset
operations and minimize electricity charges. In collaboration with Google, the
project has access to historical charging data from previous years at multiple
Bay Area workplace parking locations for Google (within PG&E’s service
territory). The locations exhibit typical workplace charging behavior in that
most session start times occur around 8 to 9 A.M. and most session end times
occur between 4 and 6 P.M. The dataset contains 15-minute interval data for
over 10,000 charging sessions and between 700 and 1000 sessions per month. The
datasets include pertinent information such as charging session start times,
end times, 15 minute average power delivered, total energy delivered, and many
others.

Challenges: There are several key challenges to designing an algorithm to
schedule the charging profiles of numerous EVs. First and foremost, the
algorithm must run in real-time without knowledge of the future EV arrivals.
The algorithm must adapt its planned power schedules as more information is
revealed (i.e., as more EVs arrive to the parking lot). Second, contrary to
most smart charging algorithms presented in the literature, the algorithm must
be able to function with limited information from each EV [6], [7].
Specifically, when an EV plugs in, the algorithm does not get access to the
EV’s State of Charge (SoC) nor does it know the EV’s future departure time
(most level 2 chargers do not sense EV SoC nor do they request user input for
future departure times). As such, our smart charging algorithm must predict
how much maximum charge an EV may consume as well as the EV’s future departure
time. Such challenges have been acknowledged in past papers including [8, 9,
10, 11]. Third, all of the EV charging schedules within a parking lot are
coupled due to the local transformer capacity constraint [12]. As such, the
algorithm cannot over-allocate power at any given time; therefore, the
algorithm should make use of a model of the future EV arrivals to avoid over
allocating power due to unexpected arrivals.

Contribution: In this paper, we present an EV smart charging algorithm for
workplace parking lots equipped with EVSEs that operates in real-time to
minimize electricity cost from time-of-use electricity rates and demand
charges while ensuring that the owners of the EVs receive adequate levels of
charge and the entire system safely operates within the local transformer
capacity constraint. Notably, our algorithm is SoC and departure time agnostic
and uses both scenario generation to account for each EV’s unknown future
departure time as well as certainty equivalent control to account for the
unknown EV arrivals in the future. We build models from the Google dataset for
each day of the week and our algorithm uses these models as the expected
future when optimizing the EV charging schedules.

##  II Problem Description

###  II-A Base Optimization

Before we discuss the smart charging algorithm that operates in real-time with
a predicted model of the future, let us consider the simpler problem of
scheduling the charging power to a single EV if we know the arrival time,
departure time, and maximum energy. Specifically, let us break one day into 15
minute segments (96 segments total) and let us denote $t=1,\dots,T$ as the
time steps each day and $T=96$. Furthermore, let $x$ be a Tx1 indicator vector
where the $t$-th element is 1 if the EV is plugged in at time $t$ and 0 if
not. Additionally, let $p$ be a Tx1 vector where the $t$-th element
corresponds to the electricity rate at time $t$ from the local grid ($/kWh).
Last, let us denote our decision variable $e$ as a Tx1 vector where the $t$-th
element indicates how much energy will be delivered to the EV in time epoch
$t$. With this notation, we can formulate the charge scheduling optimization
problem for the singular EV as follows:

$\displaystyle\hskip 8.61108pt\max_{e}$ $\displaystyle
w_{1}\log(e^{T}x+1)-w_{2}p^{T}e$ (1a) subject to: $\displaystyle 0\leq e\leq
e_{max},$ (1b) $\displaystyle e_{del\\_min}\leq e^{T}x\leq e_{del\\_max}.$
(1c)

The objective function (1a) has two weight coefficients, $w_{1}$ and $w_{2}$,
for each of the two terms. The first term corresponds to the utility that the
owner of the EV receives for the energy that their EV receives. The
logarithmic utility term was chosen to model the diminishing returns in user
utility for EVs receiving excessive amounts of energy (e.g., the first 20kWh
charged is more valuable to the EV owner than the second 20kWh). The second
term corresponds to the cost of the energy that is purchased from the local
grid (we will add demand charges in the next section). The weights, $w_{1}$
and $w_{2}$, let us adjust the relative importance of user utility from
receiving energy and energy cost. Constraint (1b) ensures that the energy
delivered at each time step is non-negative and not greater than the energy
$e_{max}$ that the EVSE can deliver in a single timestep. Constraint (1c)
ensures that the amount of energy the EV receives is greater than a minimum
amount $e_{del\\_min}$ and less than a maximum amount $e_{del\\_max}$.

##  III Real-Time Smart Charging Algorithm

In this section, we modify the base optimization presented in Section II that
will be solved at each 15 minute time step $t=1,\dots,96$ each day.
Specifically, assume that at an arbitrary time $t$, there are currenly $I$ EVs
plugged in at the location.

Departure time scenario generation: When an EV plugs in to an EVSE, we
generate $N$ potential departure times for that EV and create a scenario in
our optimization problem for each potential departure time and solve the
optimization across all scenarios. As time progresses, if a potential
departure time is no longer feasible (i.e., the potential departure time is
the current time step and the EV has not yet departed), then that scenario is
removed from the optimization via dynamic scenario weights (let $C_{n}$ be a
weight coefficient for each scenario that is set to 0 if the scenario is no
longer feasible.). Furthermore let $x_{i,n}$ be the Tx1 binary vector that
indicates when EV $i$ is available to charge in scenario $n$.

Certainty equivalent control for future EV arrivals: We make use of our
dataset to generate a model for an average day that consists of estimated
arrival times, departure times, and energy requests for each day of the week.
We then use these daily models in the real-time optimization to account for
the unknown future EV arrivals. Specifically, at time $t$, let us assume that
there are $J$ EVs in the certainty equivalent daily model that are expected to
arrive in the future. Let $x_{j}$ be the Tx1 binary vector indicating when EV
$j$ is available to charge. The decision variables that determine how much
energy is delivered at a given time $t$ are Tx1 vectors labeled as
$e_{1},\dots,e_{I}$ for the actual EVs plugged in and $e_{I+1},\dots,e_{I+J}$
for the future EV arrivals from the model. Additionally, let $p_{d}$ be the
demand charge ($/kW) to be assessed on the monthly peak load. The optimization
at time $t$ can be written as:

$\displaystyle\max_{\begin{subarray}{c}e_{1},\dots,e_{I}\\\
e_{I+1},\dots,e_{I+J}\end{subarray}}$
$\displaystyle\sum_{i=1}^{I}\sum_{n=1}^{N}\frac{1}{C_{n}}\Big{[}w_{1}\log(e_{i}^{\;T}x_{i,n}+1)-w_{2}p^{T}e_{i}\Big{]}$
(2a) $\displaystyle+$
$\displaystyle\sum_{j=I+1}^{I+J}\Big{[}w_{1}\log(e_{j}^{\;T}x_{j}+1)-w_{2}p^{T}e_{j}\Big{]}-w_{2}p_{d}\hat{e}_{inc}$
subject to: $\displaystyle 0\leq e_{k}\leq e_{max},\hskip 19.5pt\forall
k=1,\dots,I+J$ (2b) $\displaystyle e_{i}^{\;T}x_{i,n}\geq
e_{del\\_min},\;\;\forall i=1,\dots,I,$ (2c) $\displaystyle\hskip
83.95807pt\forall n=1,\dots,N,$ $\displaystyle e_{j}^{\;T}x_{j}\geq
e_{del\\_min},\hskip 6.45831pt\forall j=I+1,\dots,I+J,$ (2d)
$\displaystyle\sum_{k=1}^{I+J}e_{k}(t)\leq e_{trans},\;\;\hskip 20.0pt\forall
t=1,\dots,T,$ (2e)
$\displaystyle\hat{e}_{inc}\geq\sum_{k=1}^{I+J}e_{k}(t)-\hat{e}_{old},\hskip
7.0pt\forall t=1,\dots,T,$ (2f) $\displaystyle\hat{e}_{inc}\geq 0.$ (2g)

The first term of the objective function (2a) accounts for all $I$ EVs
currently plugged in and their $N$ potential departure times each while the
second term of the objective function accounts for all $J$ EVs in the future
model. The third term ($w_{2}p_{d}\hat{e}_{inc}$) accounts for the demand
charge from increasing the current month’s peak demand ($\hat{e}_{old}$) by
$\hat{e}_{inc}$. Constraint (2b) ensures that the energy delivered is
nonnegative and less than the EVSE max $e_{max}$. Constraint (2c) ensures a
minimum amount of energy is delivered to each EV currently plugged in,
constraint (2d) ensures a minimum amount of energy is delivered to each EV in
the future model, and constraint (2e) ensures that the sum of all energy
delivered by the EVSEs at each time $t$ does not exceed the transformer
constraint $e_{trans}$. Constraints (2f)-(2g) keep track of any increase to
the current month’s peak load for the demand charge. We note that
$\hat{e}_{old}$ corresponds to the previous peak energy demand that has been
observed during the month. The pseudocode for the daily algorithm can be
viewed in Algorithm 1 below.

Algorithm 1 Real-Time Smart Charging

1: for each day do

2: Update current parking lot state

3: for each 15 minute interval $t$ do

4: if new departure from parking lot then

5: Update parking lot state

6: end if

7: if new arrival to parking lot then

8: Generate $N$ potential departure times for new arrival

9: Update Parking lot state

10: end if

11: Formulate optimization for time $t$:

12: for each EV $i$ plugged in at time $t$ do

13: Add EV $i$ to total objective function (2a)

14: Add EV $i$ to active constraints (2b)-(2g)

15: end for

16: for each future EV $j$ in daily model $t_{model}>t$ do

17: Add EV $j$ to total objective function (2a)

18: Add EV $j$ to active constraints (2b)-(2g)

19: end for

20: Solve optimization (2a)-(2g) for time $t$

21: Store planned energy schedule for each EV $i$

22: Set each EVSE’s output power for the current 15 minute interval

23: Update peak load $\hat{e}_{old}$ for demand charge calculation (if a new
peak load is observed)

24: end for

25: end for

##  IV Case Study Results

We examine a two week period from June 17 - June 29 in 2019 at a Bay Area
workplace from our Google EV dataset. The location has 57 level 2 EVSEs with
50-100 EVs arriving each weekday and is under PG&E’s E-19 rate structure.

First, the EV charging session data was filtered by weekday and then filtered
again by arrival time. Namely, each charging session was put into one of 12
possible groups corresponding to 2-hour windows for the arrival times (e.g.,
an EV charging session that started at 9:48am would be stored in the
8:00am-10:00am group). Once this was done, daily arrival time histograms were
generated and the average stay duration and average energy consumption were
calculated for each of the 12 groups. The average arrivals per weekday, the
average arrivals per 2 hour window, the groups’ average stay durations, and
the groups’ average energy consumption were then used to create the
algorithm’s future model each day and to generate potential departure times
for each EV arrival. The daily average arrival time histogram can been viewed
in Figure 1. The 2 hour groups’ average energy consumption and average stay
duration can be viewed in Table I. Last, we note that all of these simulations
were done in Python with CVX and Mosek on a Laptop with an i7 processor and
16gb of RAM. Solving the daily smart charging optimization (2a)-(2g) at each
15 minute time step took less than a second, so the algorithm is fit to run in
real-time. Moreover, the optimization problem’s complexity is not affected by
the number of arriving EVs each day; rather, the problem size grows only as
the number of chargers increases. Additionally, for implementation, (2a)-(2g)
has to be solved every 15 minutes but for the 57 chargers in our case study,
(2a)-(2g) was solved in less than a second. Thus, the algorithm is scalable
and there is significant extra time for computation for a larger dataset
(i.e., more chargers at the parking lot).

Figure 1: Average arrival time histogram for weekdays in 2019, split into 2
hour blocks, to be used for the daily model.

Arrival Time Block | Avg Charge Amount (kWh) | Avg Stay Duration
---|---|---
12:00am-2:00am | 5.56 | 4 hrs 38 mins
2:00am-4:00am | 4.00 | 2 hrs 02 mins
4:00am-6:00am | 12.91 | 3 hrs 52 mins
6:00am-8:00am | 14.63 | 5 hrs 29 mins
8:00am-10:00am | 15.79 | 6 hrs 02 mins
10:00am-12:00pm | 9.27 | 6 hrs 02 mins
12:00pm-2:00pm | 7.41 | 11 hrs 15 mins
2:00pm-4:00pm | 6.80 | 16 hrs 06 mins
4:00pm-6:00pm | 7.14 | 16 hrs 27 mins
6:00pm-8:00pm | 6.61 | 23 hrs 19 mins
8:00pm-10:00pm | 6.78 | 25 hrs 54 mins
10:00pm-12:00am | 7.74 | 10 hrs 01 mins

TABLE I: Average charge amount (kwh) and average stay duration for different
arrival time blocks.

###  IV-A Example Day’s Charging Schedule Evolution

Figure 2 shows an example as to how the real-time optimization changes the
planned power output as time progresses and more EVs arrive. These plots
correspond to a Monday with 78 EVs arriving to 57 EVSEs with
$(w_{1},w_{2})=(2,1)$ and a transformer capacity constraint of 150kW. Starting
with the Top Left: This plot shows the planned charging schedule that is
calculated at 12:00am on the given Monday for the entire location. There are
no actual EV arrivals at this early time in the morning so the entire planned
charging schedule is created by looking at the future model for the day
(shaded blue region). Top Right: This shows the planned charging power as of
10:45am in the morning. Everything that is not in the shaded blue region
corresponds to time-steps in the past, meaning that the red plot in the white
region corresponds to actual charging power. However, everything in the blue
region is still estimated via the future model. Bottom Left: This plot shows
the planned charging schedule as of 1:30pm. Bottom Right Plot: This plot shows
the charging schedule as of 11:30pm, when there are no more arrivals for that
day. Additionally, the smart charging algorithm is able to ensure that the
actual charging schedule stays below the transformer capacity constraint.

Figure 2: Evolution of planned charging power. Top Left: At 12:00am the
algorithm is fully using the expected future as a model for what could happen
that day. Top Right: At 10:45am, many EVs have arrived and the power output
has started to deviate from the expected future. Bottom Left: Planned power
output at 1:30pm. Bottom Right: At 11:30pm, there are no more arrivals, Red
curve shows the daily power output of the 57 EVSEs. Additionally, the green
circles indicate spikes in total charging power in the evening due to the
algorithm waiting for cheap electricity rates to charge multiple EVs,
discussed in Section IV-D.

Figure 3 presents a comparison of a weekday’s predicted daily charging
schedule and the actual daily charging power that ocurred that day. As seen in
Figure 3 the daily model does a good job predicting the future load for this
location. Furthermore, note that most of the power is scheduled to be
delivered during the partial-peak electricity rates during the mid-morning
(8:30am-12:00pm) and the charging power decreases rapidly during the peak
electricity rates mid-day (12:00pm-6:00pm).

Figure 3: A comparison of the model’s predicted daily charging schedule vs.
the actual daily charging power.

###  IV-B Test Case Comparisons

Furthermore, in Table II, we present simulation results for several test
cases. In these tests, we vary $(w_{1},w_{2})$ in the objective function
between (2,1) and (10,1) to show how the algorithm swaps priority from
minimizing electricity cost to maximizing user utility for receiving energy.
Additionally, we varied the transformer constraint that couples all the EVs
together between 250kW and 100kW. In Table II, we also show the total energy
delivered over the 2 week period and the total cost of purchasing that energy
from the grid from energy rates as well as the total cost due to demand
charges (displayed as percentages compared to the status quo values). The last
2 columns in Table II show whether or not the test case included a constraint
that forces the EVSEs to charge at a certain rate for the first hour that a
new EV is plugged in. The idea behind this is to ensure that EVs will receive
some minimum amount of energy, even if they are only plugged in for a short
duration, or if they arrive during peak electricity rates. We experimented
with forcing the EVs to charge for the first hour of a new session at their
maximum rate, half their maximum rate, a quarter of their maximum rate, and
without this constraint altogether. The second to last column indicates
whether or not this constraint was included, and the last column indicates the
charge rate that was used for this first hour ($p_{max}$ = 6.6kW). The first
row of this table shows the total energy delivered and the total energy cost
for the status quo (i.e., no algorithm in place, just the energy consumed and
cost for June 17-June 29, 2019). The $(w_{1},w_{2})=(10,1)$ results in Tests
07-13 indicate that the smart charging algorithm was able to reduce energy
costs and demand charges while delivering adequate energy (greater than $80\%$
of the status quo energy) to the EVs.

Test # | $w_{1}$ | $w_{2}$ |  20cmCoupling |  |  |  |  |
---|---|---|---|---|---|---|---|---
Constraint |  |  |  |  |  |  |  |
(kW) |  20cmEnergy |  |  |  |  |  |  |
Delivered |  20cmElectricity |  |  |  |  |  |  |
Purchase Cost |  20cmDemand |  |  |  |  |  |  |
Charge Cost |  20cmForced |  |  |  |  |  |  |
Initial |  |  |  |  |  |  |  |
Charge |  20cmInitial |  |  |  |  |  |  |
Charge |  |  |  |  |  |  |  |
Rate |  |  |  |  |  |  |  |
20cmStatus Quo | n/a | n/a | n/a | 100% | 100% | 100% | n/a | n/a
01 | 2 | 1 | 250 | 50.13% | 45.78% | 64.50% | Yes | $\frac{1}{2}p_{max}$
02 | 2 | 1 | 150 | 50.01% | 45.83% | 88.76% | Yes | $\frac{1}{2}p_{max}$
03 | 2 | 1 | 125 | 50.68% | 46.33% | 73.96% | Yes | $\frac{1}{2}p_{max}$
04 | 2 | 1 | 110 | 50.44% | 45.94% | 65.09% | Yes | $\frac{1}{2}p_{max}$
05 | 2 | 1 | 110 | 44.20% | 36.40% | 65.09% | No | n/a
06 | 2 | 1 | 100 | 45.05% | 37.45% | 59.17% | No | n/a
07 | 10 | 1 | 100 | 81.93% | 81.36% | 59.17% | No | n/a
08 | 10 | 1 | 150 | 83.85% | 81.03% | 88.76% | No | n/a
09 | 10 | 1 | 250 | 84.21% | 80.98% | 89.94% | No | n/a
10 | 10 | 1 | 100 | Infeasible | Infeasible | Infeasible | Yes | $\frac{1}{2}p_{max}$
11 | 10 | 1 | 110 | 81.76% | 81.09% | 65.09% | Yes | $\frac{1}{2}p_{max}$
12 | 10 | 1 | 150 | 83.49% | 81.80% | 88.76% | Yes | $\frac{1}{2}p_{max}$
13 | 10 | 1 | 150 | Infeasible | Infeasible | Infeasible | Yes | $p_{max}$
13 | 10 | 1 | 250 | 87.25% | 84.81% | 97.63% | Yes | $p_{max}$

TABLE II: Results for 13 different test cases.

###  IV-C Transformer Capacity Constraints

As shown in Table II, we vary the transformer capacity constraint that couples
the charging power of all the EVs. In the status quo, there is no coupling
constraint and the total load peaks at 169 kW. However, with a smart charging
algorithm in place, we can constrain the total load. This would allow for the
location to use a smaller transformer capacity or increase their other non-EVs
loads (e.g., the nearby offices can safely use more power without worrying
about exceeding the transformer capacity due to the EVSEs). Figure 4 shows a 5
day comparison of the Test 08 load (see Table II) compared to the status quo
load. As seen in this plot, the smart charging algorithm was able to enforce a
transformer capacity limit at 150kW without sacrificing much energy delivered.
Furthermore, the red curve (the smart charging profile) drops below the blue
curve (the status quo) during the middle of each day to avoid the peak
electricity rates from 12:00pm-6:00pm.

Figure 4: Comparison of the Monday-Friday weekday load for the status quo
(i.e., no algorithm in place) vs. the smart charging algorithm with parameters
as listed in Test 08. Additionally, the smart charging algorithm enforces the
transformer constraint (green line).

###  IV-D Infeasibilities Due to Transformer Capacities

Additionally, we would like to discuss the effects of removing the constraint
forcing the EVSEs to charge at a certain rate for the first hour an EV is
plugged in. Recall, this constraint was added to ensure that EVs are charged
before they depart. If this constraint is removed, the smart charging
algorithm becomes overly optimistic about each EV’s departure time.
Specifically, the algorithm optimistically predicts that the EVs will stay
until the off-peak electricity rates in the evening and plans the daily
charging schedule as seen in Plots 3,4 of Figure 2 (circled in green).
Including the constraint that forces the EVSEs to charge a new arrival for the
first hour removes the second peak that is seen in Plots 3,4 of Figure 2.

The constraint that forces EVSEs to charge new EV arrivals for the first hour
of being plugged in also affects whether or not the optimization is feasible
each day. Namely, many EVs arrive around the start of the workday at 8:00am
and plug in to an EVSE. If many EVSEs are forced to charge near their maximum
rating at the same time, the total load of the parking lot can exceed the
transformer capacity. As such, in Figure 5 we show the feasible and infeasible
regions of operation for the smart charging algorithm with $w_{1}$=10 and
$w_{2}$=1.

Figure 5: The feasible and infeasible regions of operation for various
transformer capacities and forced initial charge rates for new EV arrivals.

##  V Conclusion

In this paper, we presented a smart charging algorithm for a workplace parking
lot equipped with EVSEs that operates in real-time to minimize electricity
cost from time-of-use electricity rates and demand charges while ensuring that
the owners of the EVs receive adequate levels of charge. Our algorithm uses
both scenario generation to account for each EV’s unknown future departure
time as well as certainty equivalent model predictive control to account for
the unknown EV arrivals in the future. We build models from our Google dataset
for each day of the week and our algorithm uses these models as the expected
future when optimizing the EVs charging schedules. In the future, we have
access to the meter data for loads other than the EVSEs for all locations in
our dataset as well as the single-line diagram specifications and we hope to
include more operational constraints and account for the non-EV loads in the
smart charging algorithm. Furthermore, some locations have energy storage and
solar generation which we hope to include into future algorithms. Finally,
future work will deploy this framework in real-time field operation.

Acknowledgment: This work was funded by the California Energy Commission under
grant EPC-17-020. SLAC National Accelerator Laboratory is operated for the US
Department of Energy by Stanford University under Contract DE-AC02-76SF00515.
Thank you to Rolf Schreiber from Google for providing data.

## References

  * [1] B. Kellison, “Workplace EV Charging is on the Rise,” Wood Mackenzie, https://greentechmedia.com.
  * [2] N. Tucker and M. Alizadeh, “Online pricing mechanisms for electric vehicle management at workplace charging facilities,” in _2018 56th Annual Allerton Conference_ , 2018, pp. 351–358.
  * [3] Z. Ji, X. Huang, Z. Zhang, M. Jiang, and Q. Xu, “Evaluating the vehicle-to-grid potentials by electric vehicles: A quantitative study in china by 2030,” in _2020 IEEE PESGM_ , 2020, pp. 1–5.
  * [4] B. Ferguson, V. Nagaraj, E. C. Kara, and M. Alizadeh, “Optimal planning of workplace electric vehicle charging infrastructure with smart charging opportunities,” in _2018 21st ITSC_. IEEE, 2018.
  * [5] N. Tucker and M. Alizadeh, “An online admission control mechanism for electric vehicles at public parking infrastructures,” _IEEE Transactions on Smart Grid_ , vol. 11, no. 1, pp. 161–170, 2019.
  * [6] P.-Y. Kong and G. K. Karagiannidis, “Charging schemes for plug-in hybrid electric vehicles in smart grid: A survey,” _IEEE Access_ , vol. 4, pp. 6846–6875, 2016.
  * [7] J. C. Mukherjee and A. Gupta, “A review of charge scheduling of electric vehicles in smart grid,” _IEEE Systems Journal_ , vol. 9, no. 4, 2015.
  * [8] L. Richard and M. Petit, “Fast charging station with battery storage system for ev: Optimal integration into the grid,” in _2018 IEEE Power Energy Society General Meeting (PESGM)_ , 2018, pp. 1–5.
  * [9] S. Wang, S. Bi, and Y.-J. A. Zhang, “A reinforcement learning approach for ev charging station dynamic pricing and scheduling control,” in _2018 IEEE Power Energy Society General Meeting (PESGM)_ , 2018, pp. 1–5.
  * [10] B. Khaki, Y.-W. Chung, C. Chu, and R. Gadh, “Nonparametric user behavior prediction for distributed ev charging scheduling,” in _2018 IEEE Power Energy Society General Meeting (PESGM)_ , 2018, pp. 1–5.
  * [11] Y. Xiong, B. Wang, C.-C. Chu, and R. Gadh, “Electric vehicle driver clustering using statistical model and machine learning,” in _2018 IEEE Power Energy Society General Meeting (PESGM)_ , 2018, pp. 1–5.
  * [12] S. Powell, E. C. Kara, R. Sevlian, G. V. Cezar, S. Kiliccote, and R. Rajagopal, “Controlled workplace charging of electric vehicles: The impact of rate schedules on transformer aging,” _Applied Energy_ , 2020.