# Parallel Interior-Point Solver for Block-Structured Nonlinear Programs on
SIMD/GPU Architectures

François Pacauda, Michel Schanenb, Sungho Shinb, Daniel Adrian Maldonadob,
Mihai Anitescub a Centre Automatique et Systèmes, Mines Paris - PSL, Paris,
France; b Mathematics and Computer Science Department, Argonne National
Laboratory, Lemont, USA

###### Abstract

We investigate how to port the standard interior-point method to new exascale
architectures for block-structured nonlinear programs with state equations.
Computationally, we decompose the interior-point algorithm into two successive
operations: the evaluation of the derivatives and the solution of the
associated Karush-Kuhn-Tucker (KKT) linear system. Our method accelerates both
operations using two levels of parallelism. First, we distribute the
computations on multiple processes using coarse parallelism. Second, each
process uses a SIMD/GPU accelerator locally to accelerate the operations using
fine-grained parallelism. The KKT system is reduced by eliminating the
inequalities and the state variables from the corresponding equations, to a
dense matrix encoding the sensitivities of the problem’s degrees of freedom,
drastically minimizing the memory exchange. We demonstrate the method’s
capability on the supercomputer Polaris, a testbed for the future exascale
Aurora system. Each node is equipped with four GPUs, a setup amenable to our
two-level approach. Our experiments on the stochastic optimal power flow
problem show that the method can achieve a 50x speed-up compared to the state-
of-the-art method.

##  1 Introduction

Solving complex engineering problems often resorts to the solution of large-
scale block-structured nonlinear programs. As such, there has been a long
interest in designing efficient nonlinear optimization algorithms,
particularly by using parallel computing. Parallelism can happen at two
levels. At first, _coarse parallelism_ splits the program into large
computational chunks, usually dispatched to multiple processors using a
message-passing interface in distributed memory. In this paradigm, the
parallel algorithm is designed to minimize the communication between the
different processes. In a complementary direction, _fine-grained parallelism_
breaks down the program into small tasks, fast to compute in shared memory.
This method requires a large number of processors to be efficient, and it is
usually better on SIMD architectures with low communication overhead, as
provided by Graphical Processing Units (GPUs). In the mathematical
optimization community, coarse parallelism has traditionally been used to
solve large-scale block-structured optimization problems, as encountered in
dynamic or stochastic nonlinear programs. On the contrary, fine-grained
parallelism has gained attraction only recently, with the renewed interests
for machine learning applications and stochastic gradient algorithms. In this
work, we combine coarse and fine-grained parallelism to solve block-structured
nonlinear problems on new exascale architectures, where the solution algorithm
is streamlined on different GPUs using CUDA-aware MPI.

###  1.1 Literature review

In his pioneering work [39, 40], Robert Schnabel identified three practical
approaches to run optimization algorithms in parallel: (i) parallelize the
function evaluations; (ii) parallelize the linear algebra; and (iii)
parallelize the optimization algorithm itself.

The first attempt to parallelize the evaluations has been to streamline the
computation of the derivatives using finite-differences [29]. Soon, it has
been noted that parallelizing the forward pass in automatic differentiation
(AD) is also straightforward, provided that we can propagate the tangents
(encoding the first-order sensitivity) in parallel [20]. Unfortunately, doing
the same in the reverse pass is not trivial, as adjoining a mutable code leads
to race conditions (e.g., every read becomes a write operation). This has led
to extensive research on adapting automatic differentiation to parallel
environments [4, 19, 27]. Now, most state-of-the-art differentiable tools
employ a Domain Specific Language (DSL) constraining the user to specific
differentiable operations. In particular, this approach has been adopted
mainly in machine learning, leading to the development of fast AD libraries
efficiently generating the derivatives efficiently on hardware accelerators
such as GPUs or TPUs [32, 3].

The parallelization of linear algebra is usually more involved, as most large-
scale optimization methods fall back on the solution of sparse indefinite
Karush-Kuhn-Tucker (KKT) systems [30]. In the 1980s, preliminary results were
obtained by running iterative methods in parallel, using block-Krylov [36] or
block-truncated Newton methods [28]. However, block iterative algorithms are
quickly limited by the lack of generic preconditioners for KKT systems. The
1990s witnessed the emergence of the interior-point methods (IPM), together
with the development of large-scale sparse direct linear solvers [12, 38]. In
IPM, a significant portion of the time is spent solving a sequence of
(indefinite) KKT systems, hence the method directly benefits from efficient
sparse linear solvers able to run in parallel [1, 13]. In the 2000s, it was
shown that, for block-structured optimization problems as we consider here,
the layout of the optimization problem can be exploited further in a Schur
complement approach to solve the Newton step in parallel [2, 9, 22, 17, 45,
44, 33]. These developments led to the development of mature decomposition-
based parallel nonlinear solvers for scenario-based problems in the 2010s [16,
8, 46, 34, 41].

Eventually, running an optimization algorithm fully in parallel generally
requires a subtle combination of (i) and (ii), often devolving to a software
engineering problem. The challenge is to evaluate the derivatives _and_ solve
the resulting KKT system each in parallel; all this while minimizing the
communication between the different processes. This has led to the development
of different prototypes for MPI-parallel modelers [10, 43, 21, 34], most of
them extending a specific AD backend [15, 5, 14]. Such approaches have been
successfully applied to solve large-scale block-structured nonlinear problems,
as encountered in stochastic programming and dynamic optimization.

###  1.2 Contributions

In this article, we introduce a new parallel algorithm to solve block-
structured nonlinear programs involving state equations on exascale
supercomputers. Our algorithm uses the parallel interior-point solver MadNLP
[41], using two layers of parallelism to streamline both the evaluation of the
derivatives and the solution of the KKT system. This framework targets new
exascale supercomputers, where each node is assigned to multiple GPUs
connected with a unified memory (designed to have fast memory exchange between
the different GPUs).

We demonstrate the capability of the algorithm on scenario-based power flow
problems (block-OPF), here formulated as two-stage stochastic nonlinear
programs. The scenarios can be stochastic or represent contingencies (which
can be interpreted as stochastic outcomes with uniform distribution), as is
the case of the very widely used security-constrained AC optimal power flow
(SC-ACOPF) problem [7]. SC-ACOPF is one of the core analyses undertaken in the
planning, operational planning, and real-time operation of transmission
systems [7]. SC-ACOPF is run several times a day by many operators in the US
and the world. For brevity, we will refer to such problems as stochastic.

The block structure of such problems is given by the different scenarios
associated with the stochastic problem, leading to potential parallelism in
both the evaluation of the derivatives and the solution of the resulting
block-angular KKT system. The parallel solution of the block-OPF problem with
a Schur complement approach has been studied extensively both with PIPS-NLP
[8, 37] (multiprocessing) and with Beltistos [23, 25] (multiprocessing \+
factorization of the dense Schur complement on the GPU). Compared to the
state-of-the-art solver Beltistos, our approach carries out almost all
computation on the GPUs including a global CUDA-aware MPI reduction, from the
evaluation of the derivatives to the assembling of the Schur complement. We
test our implementation on the pre-exascale supercomputer Polaris, where each
node is equipped with 4 A100 GPUs, and we solve block-OPF problems with up to
9,251 nodes.

##  2 Problem statement

In systems engineering, it is common to encounter optimization problems with
relatively few degrees of freedom – ”controls”. Then, the goal is to
appropriately fix the values for the degrees of freedom, e.g., by minimizing a
given operational cost while satisfying the physical equations of the problem.
In that context, the internal state of the system is described by a _state_
variable $x\in\mathbb{R}^{n_{x}}$, whose values depend on the current
_controls_ $u\in\mathbb{R}^{n_{u}}$ associated with the problem’s degrees of
freedom. If the problem is well-posed, this translates to the _state equation_
$g(x,u)=0$, where the function $g$ exhibits the physical structure of the
problem (e.g., a differential equation encoding a dynamics, or a nonlinear
network flow associated with static balance equations). When the system faces
uncertainties, it is often appropriate to choose a control $u$ feasible under
a finite set of conditions (or scenarios). That is, the control $u$ must
satisfy $N$ different state equations

$g_{i}(x_{i},u)=0\quad\text{for all}\quad i=1,\cdots,N,$ (1)

where the state $x_{i}$ now depend on the current scenario $i$. The variables
$x_{i}$ can be assimilated into a recourse variable. The $N$ functions
$g_{1},\cdots,g_{N}$ define the block structure of the problem.

###  2.1 Block-structured nonlinear programs

In addition to satisfying the $N$ state equations (1), we aim at minimizing
the average operating costs on the $N$ different scenarios. The corresponding
problem formulates as a two-stage nonlinear program, which, in our case, is a
nonlinear program with _partially separable structure_ [11]:

$\min_{\begin{subarray}{c}x_{1},\cdots,x_{N},\\\
u\end{subarray}}\;\sum_{i=1}^{N}f_{i}(x_{i},u)\quad\text{s.t.}\quad\left\\{\begin{aligned}
&x_{i}\geq 0\;,\quad u\geq 0\\\ &g_{i}(x_{i},u)=0\;,\\\ &h_{i}(x_{i},u)\leq
0\;,\end{aligned}\right.\quad\forall i=1,\cdots,N\;,$ (2)

with $f_{i}:\mathbb{R}^{n_{x}}\times\mathbb{R}^{n_{u}}\to\mathbb{R}$,
$g_{i}:\mathbb{R}^{n_{x}}\times\mathbb{R}^{n_{u}}\to\mathbb{R}^{n_{x}}$,
$h_{i}:\mathbb{R}^{n_{x}}\times\mathbb{R}^{n_{u}}\to\mathbb{R}^{m}$ smooth
functions encoding the objective, the state equations, and the operational
constraints, respectively. We note that the number of variables ($N\times
n_{x}+n_{u}$) and constraints ($N\times(m+n_{x})$) are linearly proportional
to the number of blocks $N$.

In addition, if we introduce local control variables $u_{1},\cdots,u_{N}$ with
the additional coupling constraint $u_{1}=\cdots=u_{N}=u$, we get a problem
with a separable structure, solvable using the primal decomposition method; at
the expense of increasing the search space [35, 11].

By introducing slack variables $s_{1},\cdots,s_{N}$, we rewrite (2) in
standard form:

$\min_{\begin{subarray}{c}x_{1},\cdots,x_{N},\\\ s_{1},\cdots,s_{N},\\\
u\end{subarray}}\;\sum_{i=1}^{N}f_{i}(x_{i},u)\quad\text{s.t.}\quad\left\\{\begin{aligned}
&u\geq 0\;,\quad x_{i}\geq 0\;,\quad s_{i}\geq 0\\\ &g_{i}(x_{i},u)=0\;,\\\
&h_{i}(x_{i},u)+s_{i}=0\;,\end{aligned}\right.\quad\forall i=1,\cdots,N\;.$
(3)

We define $y_{i}\in\mathbb{R}^{n_{x}}$ the multipliers (or _adjoints_)
associated to the equality constraints $g_{i}(x_{i},u)=0$,
$z_{i}\in\mathbb{R}^{m}$ the multipliers associated to the operational
constraints $h_{i}(x_{i},u)+s_{i}=0$, as well as $\lambda,\kappa_{i},\nu_{i}$
the three multipliers associated to the respective bound constraints $u\geq
0,x_{i}\geq 0,s_{i}\geq 0$. The Lagrangian associated to (3) is:

$L(x,u,s;y,z,\lambda,\mu,\nu):=\\\
\sum_{i=1}^{N}\Big{[}f_{i}(x_{i},u)+y_{i}^{\top}g_{i}(x_{i},u)+z_{i}^{\top}\big{(}h_{i}(x_{i},u)+s_{i}\big{)}-\kappa_{i}x_{i}-\nu_{i}s_{i}\Big{]}-\lambda
u\;,$ (4)

with $x:=(x_{1},\cdots,x_{N})$, $s:=(s_{1},\cdots,s_{N})$,
$y:=(y_{1},\cdots,y_{N})$, $z:=(z_{1},\cdots,z_{N})$. To simplify the
notations, we define the _extended_ objective function and the _extended_
constraints:

$f(x,u):=\sum_{i=1}f_{i}(x_{i},u)\;,\quad
g(x,u):=\begin{bmatrix}g_{1}(x_{1},u)\\\ \vdots\\\
g_{N}(x_{N},u)\end{bmatrix}\;,\quad h(x,u):=\begin{bmatrix}h_{1}(x_{1},u)\\\
\vdots\\\ h_{N}(x_{N},u)\end{bmatrix}.$

We assume the functions $f,g,h$ are twice differentiable. We denote

$\displaystyle H=\partial_{(x,u)}h(x,u)\in\mathbb{R}^{Nm\times(Nn_{x}+n_{u})}$
Jacobian of the inequality cons. $\displaystyle
G=\partial_{(x,u)}g(x,u)\in\mathbb{R}^{Nn_{x}\times(Nn_{x}+n_{u})}$ Jacobian
of the equality cons. $\displaystyle
W=\nabla_{(x,u)}^{2}L(x,u,s;\cdot)\in\mathbb{R}^{(Nn_{x}+n_{u})\times(Nn_{x}+n_{u})}$
Hessian of Lagrangian.

###  2.2 Interior-point method

The interior-point method (IPM) [30, Chapter 19] is a classical approach to
solve (3).

####  2.2.1 KKT system

The Karush-Kuhn-Tucker (KKT) equations associated to (3) can be expressed as

$\displaystyle\nabla_{x}f_{i}+(G_{x}^{i})^{\top}y_{i}+(H_{x}^{i})^{\top}z_{i}-\kappa_{i}=0,$
$\displaystyle\forall i=1,\cdots,N$ (5a)
$\displaystyle\sum_{i=1}^{N}\Big{(}\nabla_{u}f_{i}+(G_{u}^{i})^{\top}y_{i}+(H_{u}^{i})^{\top}z_{i}\Big{)}-\lambda=0,$
(coupling) (5b) $\displaystyle z_{i}-\nu_{i}=0,$ $\displaystyle\forall
i=1,\cdots,N$ (5c) $\displaystyle g_{i}(x_{i},u)=0,$ $\displaystyle\forall
i=1,\cdots,N$ (5d) $\displaystyle h_{i}(x_{i},u)+s_{i}=0,$
$\displaystyle\forall i=1,\cdots,N$ (5e) $\displaystyle
X_{i}\kappa_{i}=0,\;(x_{i},\kappa_{i})\geq 0,$ $\displaystyle\forall
i=1,\cdots,N$ (5f) $\displaystyle S_{i}\nu_{i}=0,\;(s_{i},\nu_{i})\geq 0,$
$\displaystyle\forall i=1,\cdots,N$ (5g) $\displaystyle
U\lambda=0,\;(u,\lambda)\geq 0,$ (5h)

where $U=\text{diag}(u)$, $X_{i}=\text{diag}(x_{i})$,
$S_{i}=\text{diag}(s_{i})$.

The interior-point method uses a homotopy parameter $\mu>0$ to replace the
complementarity constraints (5f)-(5g)-(5h) by the smooth approximations:
$X_{i}\kappa_{i}=\mu e_{n_{x}},S_{i}\nu_{i}=\mu e_{m},U\lambda=\mu e_{n_{u}}$
($e_{n}$ being the vector of all ones of dimension $n$). The resulting
(smooth) system of nonlinear equations can be solved iteratively using Newton
method, where at each iteration, the descent direction is updated by solving
the following _augmented_ linear system:

$\begin{bmatrix}W+\Sigma_{p}&0&G^{\top}&H^{\top}\\\ 0&\Sigma_{s}&0&I\\\
G&0&0&0\\\ H&I&0&0\end{bmatrix}\begin{bmatrix}p_{d}\\\ p_{s}\\\ p_{y}\\\
p_{z}\end{bmatrix}=-\begin{bmatrix}r_{1}\\\ r_{2}\\\ r_{3}\\\
r_{4}\end{bmatrix}$ (6)

with $r_{1}=\begin{bmatrix}\nabla_{x}f+G_{x}^{\top}y+H_{x}^{\top}z-\mu
X^{-1}e_{n_{x}}\\\ \nabla_{u}f+G_{u}^{\top}y+H_{u}^{\top}z-\mu
U^{-1}e_{n_{u}}\end{bmatrix}$, $r_{2}=z-\mu S^{-1}e_{m}$, $r_{3}=g(x,u)$,
$r_{4}=h(x,u)+s$. The primal descent direction $p_{d}$ decomposes as
$p_{d}=(p_{x_{1}},\cdots,p_{x_{N}},p_{u})$.

####  2.2.2 Block angular structure

The linear system (6) is sparse and symmetric indefinite, and can be
factorized using the Bunch-Kaufman algorithm. However, it is often beneficial
to exploit its block-angular structure. Indeed, both the Hessian of the
Lagrangian and the Jacobians have a block-angular structure, given as

$W=\begin{bmatrix}W_{x_{1}x_{1}}&&&W_{x_{1}u}\\\ &\ddots&&\vdots\\\
&&W_{x_{N}x_{N}}&W_{x_{N}u}\\\
W_{ux_{1}}&\ldots&W_{ux_{N}}&W_{uu}\end{bmatrix},\quad
G=\begin{bmatrix}G^{1}_{x_{1}}&&&G^{1}_{u}\\\ &\ddots&&\vdots\\\
&&G^{N}_{x_{N}}&G^{N}_{u}\end{bmatrix}\;.$

By reordering the linear system (6), we can expose the block-angular structure
of the KKT system as:

$\begin{bmatrix}A_{1}&&&B_{1}^{\top}\\\ &\ddots&&\vdots\\\
&&A_{N}&B_{N}^{\top}\\\ B_{1}&\ldots&B_{N}&A_{0}\end{bmatrix}$ (7)

with

$A_{0}=W_{uu},\quad
A_{i}=\begin{bmatrix}W_{x_{i}x_{i}}+\Sigma_{x_{i}}&0&G_{x_{i}}^{\top}&H_{x_{i}}^{\top}\\\
0&\Sigma_{s_{i}}&0&I\\\ G_{x_{i}}&0&0&0\\\ H_{x_{i}}&I&0&0\end{bmatrix},\quad
B_{i}=\begin{bmatrix}W_{x_{i}u}\\\ (G_{u}^{i})^{\top}\\\
(H_{u}^{i})^{\top}\end{bmatrix}^{\top}\;.$

The block-angular structure (7) can be exploited to solve the KKT linear
system in parallel using a Schur complement approach. In that case, the
submatrices $A_{i}$ can be factorized independently to assemble the Schur
complement in parallel [8].

###  2.3 Condensation and reduction

Instead of reordering the augmented KKT system (6) as a block angular matrix
(7), we propose an alternative approach based on successive condensation and
reduction of the KKT system, following the method introduced in [31]. If the
structure is well-defined, we show that we can condense the KKT system (6) to
a dense matrix with size $n_{u}\times n_{u}$ in two steps: first, by removing
the inequality constraints in (6), then by exploiting the structure of the
equality constraints to reduce the condensed system to a dense matrix. The
condensation and reduction steps are illustrated in Figure 1.

Figure 1: Successive reductions for a block-structured nonlinear problem with
$N=3$: Augmented system (6), Condensed system (8), Reduced system (11).

####  2.3.1 Condensation step

The condensation step allows reducing the size of the KKT system drastically
if the number of inequality constraints is large111It is equivalent to the
normal equations in linear programming [30, Chapter 16, p.412].

######  Proposition 2.1 (Condensed KKT system).

The linear system (6) is equivalent to

$\begin{bmatrix}K+\Sigma_{p}&G^{\top}\\\
G&0\end{bmatrix}\begin{bmatrix}p_{d}\\\
p_{y}\end{bmatrix}=-\begin{bmatrix}r_{1}+H^{\top}(\Sigma_{s}r_{4}-r_{2})\\\
r_{3}\end{bmatrix}\;,$ (8)

where $K\in\mathbb{R}^{(Nn_{x}+n_{u})\times(Nn_{x}+n_{u})}$ is the _condensed
matrix_ $K:=W+H^{\top}\Sigma_{s}H$. The descent directions $p_{s}$ and $p_{z}$
are recovered as

$\left\\{\begin{aligned}
&p_{z}=\Sigma_{s}\big{[}Hp_{d}+r_{4}\big{]}-r_{2}\;,\\\
&p_{s}=-\Sigma_{s}^{-1}\big{[}r_{2}+p_{z}\big{]}\;.\end{aligned}\right.$ (9)

###### Proof.

See [31, Theorem 2.2]. ∎

The condensed matrix $K$ inherits the block-angular structure of the Hessian
of the Lagrangian $W$.

######  Proposition 2.2.

The condensed matrix $K=W+H^{\top}\Sigma_{s}H$ has a block-angular structure,
given as

$K=\begin{bmatrix}K_{x_{1}x_{1}}&&&K_{x_{1}u}\\\ &\ddots&&\vdots\\\
&&K_{x_{N}x_{N}}&K_{x_{N}u}\\\
K_{ux_{1}}&\ldots&K_{ux_{N}}&K_{uu}\end{bmatrix}$ (10)

where we have defined the condensed blocks
$K_{x_{i}x_{i}}:=W_{x_{i}x_{i}}+(H_{x_{i}}^{i})^{\top}\Sigma_{s_{i}}H^{i}_{x_{i}}$,
$K_{ux_{i}}:=W_{ux_{i}}+(H^{i}_{u})^{\top}\Sigma_{s_{i}}H^{i}_{x_{i}}$ and
$K_{uu}:=W_{uu}+\sum_{i=1}^{N}(H^{i}_{u})^{\top}\Sigma_{s_{i}}H^{i}_{u}$.

###### Proof.

This is proved by induction. ∎

####  2.3.2 Reduction step

In addition, we can exploit the structure of the equality constraints
$g_{1},\cdots,g_{N}$ to further reduce the size of the linear system (8) down
to a dense matrix with size $n_{u}\times n_{u}$. Equation (10) exhibits the
structure w.r.t. the state $x$ and the control $u$, we rewrite as such the
condensed KKT system (8) as

$\begin{bmatrix}K_{x_{1}x_{1}}&&&K_{x_{1}u}&(G^{1}_{x_{1}})^{\top}&&\\\
&\ddots&&\vdots&&\ddots&\\\
&&K_{x_{N}x_{N}}&K_{x_{N}u}&&&(G^{N}_{x_{N}})^{\top}\\\
K_{ux_{1}}&\ldots&K_{ux_{N}}&K_{uu}&(G^{1}_{u})^{\top}&\ldots&(G^{1}_{u})^{\top}\\\
G^{1}_{x_{1}}&&&G^{1}_{u}&&&\\\ &\ddots&&\vdots&&&\\\
&&G^{N}_{x_{N}}&G^{N}_{u}&&&\\\ \end{bmatrix}\begin{bmatrix}p_{x_{1}}\\\
\vdots\\\ p_{x_{N}}\\\ p_{u}\\\ p_{y}^{1}\\\ \vdots\\\
p_{y}^{N}\end{bmatrix}=-\begin{bmatrix}\hat{r}_{1}^{1}\\\ \vdots\\\
\hat{r}_{1}^{N}\\\ \hat{r}_{2}\\\ \hat{r}_{3}^{1}\\\ \vdots\\\
\hat{r}_{3}^{N}\end{bmatrix}\;,$

where we have renamed the right-hand-side in (8) as $\hat{r}$.

######  Proposition 2.3 (Reduction).

Assume that for all $i=1,\cdots,N$ the Jacobian matrices
$G_{x}^{i}\in\mathbb{R}^{n_{x}\times n_{x}}$ are invertible. Then the linear
system (8) is equivalent to

$\hat{K}_{uu}\,p_{u}=-\hat{r}_{2}+\sum_{i=1}^{N}\Big{[}(G_{u}^{i})^{\top}(G_{x}^{i})^{-\top}\hat{r}_{1}^{i}+\big{[}K_{ux_{i}}-(G_{u}^{i})^{\top}(G_{x}^{i})^{-\top}K_{x_{i}x_{i}}\big{]}(G_{x}^{i})^{-1}\hat{r}_{3}^{i}\Big{]}$
(11)

with $\hat{K}_{uu}:=Z^{\top}KZ$ and $Z\in\mathbb{R}^{(n_{u}+Nn_{x})\times
n_{u}}$ is the reduction operator defined as

$Z=\begin{bmatrix}-(G_{x}^{1})^{-1}G_{u}^{1}\\\ \vdots\\\
-(G_{x}^{N})^{-1}G_{u}^{N}\\\ I\end{bmatrix}\;.$ (12)

The descent directions $p_{x}$ and $p_{y}$ are recovered as

$\left\\{\begin{aligned}
&p_{x}^{i}=-(G_{x}^{i})^{-1}\big{[}\hat{r}_{3}^{i}+G_{u}^{i}p_{u}\big{]}\\\
&p_{y}^{i}=-(G_{x}^{i})^{-\top}\big{[}\hat{r}_{1}^{i}+K_{x_{i}x_{i}}p_{x}^{i}+K_{x_{i}u}p_{u}\big{]}\;.\end{aligned}\right.$
(13)

###### Proof.

See [31, Theorem 2.1]. ∎

The reduction (11) is equivalent to a Schur complement approach applied to the
condensed KKT system (8). In Proposition (2.1), we have shown that the
condensed matrix $K$ has a block-angular structure. The associated condensed
KKT system (8) is also inheriting a block-angular structure in the form of
(7), where the blocks are given by

$A_{0}=K_{uu}\;,\quad
A_{i}=\begin{bmatrix}K_{x_{i}x_{i}}&(G_{x}^{i})^{\top}\\\
G_{x}^{i}&0\end{bmatrix}\;,\quad B_{i}=\begin{bmatrix}K_{x_{i}u}\\\
G_{u}^{i}\end{bmatrix}^{\top}\;.$ (14)

######  Proposition 2.4.

Assume that for each $i=1,\cdots,N$ the Jacobian $G_{x}^{i}$ is invertible.
Let $S_{uu}=A_{0}-\sum_{i=1}^{N}B_{i}A_{i}^{-1}B_{i}^{\top}$ be the Schur
complement associated to the block-angular system (7) with the matrices
$(A_{i},B_{i})$ defined in (14). Then, the Schur complement $S_{uu}$ is equal
to the reduced matrix $\widehat{K}_{uu}$ defined in (11): $S_{uu}=Z^{\top}KZ$.

###### Proof.

First, note that if the Jacobian $G_{x}^{i}$ is invertible, then the block
matrix $A_{i}$ defined in (14) is also invertible, with

$A_{i}^{-1}=\begin{bmatrix}0&(G_{x}^{i})^{-1}\\\
(G_{x}^{i})^{-\top}&-(G_{x}^{i})^{-\top}K_{x_{i}x_{i}}(G_{x}^{i})^{-1}\end{bmatrix}\;.$
(15)

Using (14)-(15), we expand the expression of the terms in the sum constituting
the Schur complement $S_{uu}$:

$\displaystyle B_{i}A_{i}^{-1}B_{i}^{\top}$
$\displaystyle=\begin{bmatrix}K_{ux_{i}}&(G_{u}^{i})^{\top}\end{bmatrix}\begin{bmatrix}0&(G_{x}^{i})^{-1}\\\
(G_{x}^{i})^{-\top}&-(G_{x}^{i})^{-\top}K_{x_{i}x_{i}}(G_{x}^{i})^{-1}\end{bmatrix}\begin{bmatrix}K_{x_{i}u}\\\
G_{u}^{i}\end{bmatrix}\;,$
$\displaystyle=(G_{u}^{i})^{\top}(G_{x}^{i})^{-\top}K_{x_{i}u}+K_{ux_{i}}(G_{x}^{i})^{-1}(G_{u}^{i})-(G_{u}^{i})^{\top}(G_{x}^{i})^{-\top}K_{x_{i}x_{i}}(G_{x}^{i})^{-1}G_{u}^{i}\;.$

Hence, the Schur complement
$S_{uu}=A_{0}-\sum_{i=1}^{N}B_{i}A_{i}^{-1}B_{i}^{\top}$ expands as

$\displaystyle S_{uu}$
$\displaystyle=K_{uu}-\sum_{i=1}^{N}\big{[}(G_{u}^{i})^{\top}(G_{x}^{i})^{-\top}K_{x_{i}u}+K_{ux_{i}}(G_{x}^{i})^{-1}(G_{u}^{i})-(G_{u}^{i})^{\top}(G_{x}^{i})^{-\top}K_{x_{i}x_{i}}(G_{x}^{i})^{-1}G_{u}^{i}\big{]}$
$\displaystyle=Z^{\top}KZ\;.$

We recover the expression of the reduced matrix $\widehat{K}_{uu}$ in
Proposition 2.3. ∎

###  2.4 Discussion

Hence, we can interpret the reduction step as a Schur complement approach.
Forming the Schur complement has always been the bottleneck when solving
distributed block angular problems in parallel [26, 8]. Its reduction
operation involves large memory transfers between the processes, with the
number of transfers being on the order of $\mathcal{O}(log(p))$, where $p$ is
the number of processes. Due to the quasi-shared memory architecture on GPUs,
the reduction can be implemented efficiently [31]. In the next section, we
propose to extend [31] to assemble the reduced matrix $\widehat{K}_{uu}$ using
two levels of parallelism, using both MPI and CUDA, thus reducing the reliance
on distributed memory.

##  3 Parallel implementation

In the previous section, we have detailed the structure of block-angular
nonlinear programs and presented the condensation and reduction steps for the
KKT system. The loose coupling between the blocks is favorable for
parallelizing the evaluation of the derivatives and the solution of the block-
angular KKT system. _Globally_ , we can distribute the computation on
different processes using MPI (coarse parallelism). _Locally_ , we can further
streamline the computation using GPU accelerators (fine-grained parallelism).
This paradigm, with its two levels of parallelism, is directly in line with
what is currently offered by the new exascale architectures, where each node
has 4 to 8 GPUs, all sharing a unified memory for fast communication. We
present in §3.1 how we streamline the evaluation of the model using automatic
differentiation, and in §3.2 how we parallelize the solution of the KKT
system.

###  3.1 Parallel automatic-differentiation

First, we present how to evaluate the model in parallel using automatic
differentiation [18]. We illustrate the procedure in Figure 2. The goal of the
algorithm is to streamline the evaluation of the $N$ scenarios on $N/M$ GPUs,
$M$ being the number of scenarios evaluated locally on each GPU (we suppose
here that $N$ is a multiple of $M$).

root$g_{1},\cdots,g_{4}$$g_{13},\cdots,g_{16}$$g_{1}$$g_{4}$$\cdots$$g_{13}$$g_{16}$$\cdots$$\cdots$$\cdots$
Figure 2: Parallel evaluation of the derivatives for $g_{1},\cdots g_{N}$ on 4
GPUs: we have a total of $N=16$ scenarios, each GPU evaluating $M=16/4=4$
scenarios locally.

####  3.1.1 Local parallelism

The first level of parallelism streamlines the evaluation of the model on
SIMD/GPU devices. We have designed our implementation to run entirely on the
GPU device, to avoid any data transfer between the host and the device.

##### Block evaluation

We suppose that the nonlinear functions $(f_{i},g_{i},h_{i})$ share the same
structure, its expressions yielding the same Abstract Syntax Tree (AST) for
all $i=1,\cdots,M$. We illustrate the block evaluation on a simple abstract
tree, but the reasoning extends to more complicated structures. We suppose
that for all $i$, the functions $f_{i},g_{i},h_{i}$ depend linearly on a
nonlinear basis matrix
$\psi:\mathbb{R}^{n_{x}}\times\mathbb{R}^{n_{u}}\to\mathbb{R}^{n_{b}}$: that
is, there exists three _sparse_ matrices $L_{f},L_{g},L_{h}$ such that

$f_{i}(x_{i},u)=L_{f}\psi(x_{i},u)\;,\quad
g_{i}(x_{i},u)=L_{g}\psi(x_{i},u)\;,\quad
h_{i}(x_{i},u)=L_{h}\psi(x_{i},u)\;.$ (16)

Suppose we aim to evaluate the $M$ functions $g_{1},\cdots,g_{M}$ in batch for
the states $x_{1},\cdots,x_{M}$. The structure (16) is directly amenable for
SIMD evaluation. We denote by
$X_{M}=(x_{1},\cdots,x_{M})\in\mathbb{R}^{n_{x}\times M}$ the dense matrix
obtained by concatenating the $M$ states together. By using a proper GPU
kernel or a parallel modeler, we can evaluate the basis in a SIMD fashion and
build the matrix
$\Psi(X_{M},u):=\big{(}\psi(x_{1},u),\cdots,\psi(x_{M},u)\big{)}\in\mathbb{R}^{n_{b}\times
M}$. Then, evaluating the functions $g_{1},\cdots,g_{M}$ simultaneously
translates to the evaluation of one SpMM product:

$\big{(}g_{1}(x_{1},u),\cdots,g_{M}(x_{M},u)\big{)}=L_{g}\Psi(X_{M},u)\in\mathbb{R}^{n_{x}\times
M}\;.$ (17)

The total memory required in the two successive operations is
$\mathcal{O}((n_{x}+n_{b})\times M)$, and depends linearly on the number of
blocks $M$. We note the SpMM operations are generally implemented efficiently
in the vendor library (cusparse for CUDA, rocSPARSE for AMDGPU).

##### First-order derivatives

Suppose that for a given $i$ we have a differentiable implementation
$\tt{gb}_{i}:\mathbb{R}^{n_{d}}\to\mathbb{R}^{n_{x}}$ associated to the
function $g_{i}$. We aim to evaluate the Jacobian-matrix products $(\nabla
g_{i})D$ for $p$ tangents encoded in a matrix $D\in\mathbb{R}^{n_{d}\times p}$
using forward-mode AD and operator overloading. This operation translates to
propagating _forward_ a vector of dual numbers. Denoting by
$\underline{d}\in\mathbb{D}_{p}^{n_{d}}$ the dual number encoding the $p$
tangents stored in $D$, evaluating $(\nabla g_{i})D$ simply amounts to call
$\tt{gb}_{i}(\underline{d})$ and extract the results in the dual numbers
returned as a result. As $G^{i}$ is sparse, we can apply the technique of
Jacobian coloring [18] to compress the independent columns of the sparse
matrix $G^{i}$ and reduces the number of required seeding tangents $p$ needed
to evaluate the full Jacobian.

Suppose now we want to evaluate the sparse Jacobians $G^{1},\cdots,G^{M}$ in
batch. As the functions $g_{i}$ are based on the same AST, their respective
Jacobians $G^{1},\cdots,G^{M}$ are sharing the same sparsity pattern. By
seeding a matrix of dual numbers
$\underline{D}_{M}=(\underline{d}_{1},\cdots,\underline{d}_{M})\in\mathbb{D}_{p}^{n_{d}\times
M}$, we can use the same operation as (17) to streamline the evaluation of the
$M$ Jacobian-vector products using the SIMD kernel $\Psi(\cdot)$ and SpMM
operations:

$\big{(}\text{\tt gb}_{1}(\underline{d}_{1}),\cdots,\text{\tt
gb}_{M}(\underline{d}_{M})\big{)}:=L_{g}\Psi(\underline{D}_{M})\in\mathbb{D}_{p}^{n_{x}\times
M}\;.$ (18)

Once the results are evaluated, it remains to uncompress the dual outputs to
build the $M$ sparse Jacobians $G^{1},\cdots,G^{M}$. Hence, we can streamline
the evaluation of the Jacobian along with the number of tangents $p$ and the
number of blocks $M$. This comes at the expense of increasing memory usage to
$\mathcal{O}((n_{x}+n_{b}+n_{d})\times M\times p)$ (to store the dual matrices
associated to the input, the intermediate basis $\Psi$ and the output).

##### Second-order derivatives

The evaluation of the second-order derivatives follows the same procedure,
using _forward-over-reverse_ AD. For each $i$, we suppose available an adjoint
function $\text{\tt
adj\\_gb}_{i}:\mathbb{R}^{n_{d}}\times\mathbb{R}^{n_{x}}\to\mathbb{R}^{n_{d}}$
which for any primal $x\in\mathbb{R}^{n_{d}}$ and adjoint
$y\in\mathbb{R}^{n_{x}}$ evaluates the Jacobian-transpose vector product
$(G^{i}(x))^{\top}y$ (reverse-mode). Using forward-mode AD on top of
$\text{\tt adj\\_gb}_{i}$, we can compute the second-order derivatives
$y^{\top}\nabla^{2}g_{i}(x)V$ for $p$ directions $V$ by calling $\text{\tt
adj\\_gb}_{i}(\underline{x},\underline{y})$. Using Hessian coloring, we can
compress the independent columns of the sparse matrix $y^{\top}\nabla^{2}g(x)$
and reduce the number of seeding tangents $p$ required to evaluate the full
Hessian. We note that in general obtaining an adjoint $\text{\tt
adj\\_gb}_{i}$ running in parallel is nontrivial due to potential race
conditions incurred by the control flow reversal of the original code.

Computing the Hessian $y^{\top}\nabla^{2}g_{i}(x)$ in parallel for $i=1,\cdots
M$ amounts to defining two matrices of dual numbers
$\underline{X}_{M}=(\underline{X}_{1},\cdots,\underline{X}_{M})\in\mathbb{D}_{p}^{n_{d}\times
M}$,
$\underline{Y}_{M}=(\underline{y}_{1},\cdots,\underline{y}_{M})\in\mathbb{D}_{p}^{n_{x}\times
M}$ and evaluate
$\nabla\Psi(\underline{X}_{M})^{\top}L_{g}^{\top}\underline{Y}_{M}$. The dual
outputs are uncompressed to build the $M$ sparse Hessians (as the sparsity
pattern of the Hessians is different than those of the Jacobians, the matrix
$\underline{X}_{M}$ employed here is different than the one used in (18)). The
total memory required to store the duals is
$\mathcal{O}((2n_{x}+n_{d}+n_{b})\times M\times p)$. For more details, we
refer to the vector forward mode as described in [18].

####  3.1.2 Global parallelism

Now, if we have several GPUs at our disposal, we can push the parallelism
further by distributing the evaluations using multiprocessing and a Message
Passing Interface (MPI) library. Coming back at our original problem (2), we
illustrate in Figure 2 how to dispatch the evaluation of the $N$ nonlinear
constraints $g_{1},\cdots,g_{N}$ (the same reasoning applies to the objectives
$f_{1},\cdots,f_{N}$ and the inequality constraints $h_{1},\cdots,h_{N}$). We
use the streamlined implementation described in the previous subsection to
evaluate the constraints in a batch of size $M$: the first GPU evaluates the
constraints $g_{1},\cdots,g_{M}$, the second GPU evaluates
$g_{M+1},\cdots,g_{2M}$, and so on. In total, the evaluation of the $N$
constraints requires $N/M$ GPUs (if $M=1$, each GPU evaluate one constraint;
if $M=N$, we use only one GPU evaluating all the constraints).

The implementation has been designed to minimize the communication between the
different processes: each batch $g_{1},\cdots,g_{M}$ stores the data it needs
_locally_ , the only data exchange with the other processes being the vector
of input and the vector of output. In addition, we will see in the next
section we do not have to transfer the first- and second-order information if
a parallel linear solver is being used.

###  3.2 Parallel KKT solver

By exploiting the block-angular structure of the KKT system, we can solve the
Newton step in parallel using a Schur complement approach. The challenge lies
in the computation of the Schur complement matrix
$S=A_{0}-\sum_{i=1}^{N}B_{i}A_{i}^{-1}B_{i}^{\top}$. Each product
$B_{i}A_{i}^{-1}B_{i}^{\top}$ requires the factorization of the matrix $A_{i}$
and the solution of a linear system with multiple (sparse) right-hand-side
$A_{i}^{-1}B_{i}$. State-of-the-art methods are evaluating the Schur
complement using an _incomplete augmented factorization_ applied on the
auxiliary matrix $\begin{bmatrix}A_{i}&B_{i}^{\top}\\\ B_{i}&0\end{bmatrix}$,
as currently implemented in the Pardiso linear solver [33]. Here, we use an
alternative approach building on the reduced KKT system §2.3.2 (equivalent to
the Schur complement approach). As the reduction can be streamlined on GPU
accelerators [31], this approach can assemble the Schur complement in parallel
using CUDA-aware MPI. We illustrate the parallel computation of the Schur
complement in Figure 3.

root(Assembling)$K,G$$K,G$$K,G$$K,G$ADADADAD$\hat{K}_{uu}^{1:4}$$\hat{K}_{uu}^{5:8}$$\hat{K}_{uu}^{9:12}$$\hat{K}_{uu}^{13:16}$(Reduction)$Z^{\top}KZ$$Z^{\top}KZ$$Z^{\top}KZ$$Z^{\top}KZ$MPI_AllReduce+(Schur
compl.)$\hat{K}_{uu}$ Figure 3: Parallel computation of the Schur complement.

##### Assembling the sparse matrices.

Using the procedure introduced in §3.1.1, we evaluate _locally_ the Jacobians
$G^{1},\cdots,G^{M}$ the Jacobians $H^{1},\cdots,H^{M}$ and the Hessians
$W^{1},\cdots,W^{M}$. Using dedicated kernels, we uncompress the results in
the block-angular sparse Jacobians

$G_{x}^{1:M}=\begin{bmatrix}G_{x_{1}}^{1}&&\\\ &\ddots&\\\
&&G_{x_{M}}^{M}\end{bmatrix},~{}G_{u}^{1:M}=\begin{bmatrix}G_{u}^{1}\\\
\vdots\\\
G_{u}^{M}\end{bmatrix},~{}H^{1:M}=\begin{bmatrix}H_{x_{1}}^{1}&&&H_{u}^{1}\\\
&\ddots&&\vdots\\\ &&H_{x_{M}}^{M}&H_{u}^{M}\end{bmatrix}\;,$

and sparse Hessian

$W^{1:M}=\begin{bmatrix}W_{x_{1}x_{1}}&&&W_{x_{1}u}\\\ &\ddots&&\vdots\\\
&&W_{x_{M}x_{M}}&W_{x_{M}u}\\\
W_{ux_{1}}&\ldots&W_{ux_{M}}&W_{uu}\end{bmatrix}\;.$

Once the sparse matrices are obtained, we recover the condensed matrix
$K^{1:M}=W^{1:M}+(H^{1:M})^{\top}\Sigma(H^{1:M})$ (Proposition 2.1) using one
SpGEMM operation and we factorize the matrix $G_{x}^{1:M}$ using a sparse LU
factorization (potentially running in batch as the matrices
$G_{x_{1}}^{1},\cdots,G_{x_{M}}^{M}$ are sharing the same sparsity pattern).
Once the matrix $G_{x}^{1:M}$ is factorized as $PG_{x}^{1:M}Q=LU$ ($P$, $Q$
being two permutation matrices), computing $(G_{x}^{1:M})^{-1}b$ translates to
two backsolves (SpSV) and two matrix-vector multiplications (SpMV), as
$(G_{x}^{1:M})^{-1}b=QU^{-1}L^{-1}Pb$.

##### Local reduction.

Once the sparse matrices are built, we evaluate locally the reduced matrix
$\widehat{K}_{uu}^{1:M}$ on the GPU, using $\text{div}(n_{u},n_{batch})+1$
matrix-matrix product $\widehat{K}_{uu}^{1:M}V$ (with
$V\in\mathbb{R}^{n_{u}\times n_{batch}}$ a dense matrix encoding $n_{batch}$
vectors of the Cartesian basis of $\mathbb{R}^{n_{u}}$). The evaluation of one
batched matrix-matrix product $\widehat{K}_{uu}^{1:M}V=(Z^{\top}K^{1:M}Z)V$
proceeds in three steps

  1. 1.

Solve $T_{x}=-(G_{x}^{1:M})^{-1}(G_{u}^{1:M}V)$.

  2. 2.

Evaluate $\begin{bmatrix}L_{x}\\\
L_{u}\end{bmatrix}:=\begin{bmatrix}K_{xx}^{1:M}&K_{xu}^{1:M}\\\
K_{ux}^{1:M}&K_{uu}^{1:M}\end{bmatrix}\begin{bmatrix}T_{x}\\\ V\end{bmatrix}$.

  3. 3.

Set $\widehat{K}_{uu}^{1:M}V=L_{u}-G_{u}^{1:M}(G_{x}^{1:M})^{-\top}L_{x}$.

In total, we need 2 SpSM and 3 SpMM operations in the first step, 1 SpMM in
the second step, and 2 SpSM and 3 SpMM operations in the third step, giving a
total of 4 SpSM and 7 SpMM operations. More than the computation, the
reduction is limited by the memory, as we have to store the three buffers
$L_{x},T_{x},T_{u}$ with a total size of $(2M\times n_{x}+n_{u})\times
n_{batch}$. If $n_{x}$ is too large, it is in our interest to reduce $M$ (by
using more GPUs) or to reduce $n_{batch}$ (at the expense of computing more
matrix-matrix product $\widehat{K}_{uu}^{1:M}V$).

##### Global reduction.

Once we obtain the locally reduced matrices $\widehat{K}_{uu}^{nM+1:(n+1)M}$
for $n=0,\cdots,N/M-1$, we can assemble the global reduced matrix
$\widehat{K}_{uu}=\sum_{n=0}^{N/M-1}\widehat{K}_{uu}^{nM+1:(n+1)M}$ using one
_all reduce_ (MPI_Allreduce) operation. The size of the reduced matrix
$\widehat{K}_{uu}$ is $n_{u}\times n_{u}$, hence limiting the memory transfer
required in the algorithm.

###  3.3 Discussion

We have presented a practical way to assemble the Schur complement on multi-
GPU architectures. The parallelism occurs both at the _local_ level (SIMD
evaluations on the GPUs) and at the _global_ level (distributed computation
with MPI). The algorithm has the advantage of assembling the sparse Jacobians
and Hessians only locally, as the reduction occurs before proceeding to the
memory transfer with MPI_Allreduce. The reduced matrix has a dimension
$n_{u}\times n_{u}$, which compresses the memory transfer significantly if the
number of degrees of freedom $n_{u}$ is small. However, this comes at the
expense of storing a vector of dual numbers (whose memory is linearly
proportional to the number of blocks $M$ evaluated locally and the number of
tangents $p$ being employed to evaluate the sparse derivatives) and additional
buffers in the reduction algorithm. In the next section, we will test an
implementation of the algorithm on CUDA GPUs, and show that the algorithm is
practical.

##  4 Numerical results

We demonstrate the capabilities of the algorithm we introduced in Section §3
on the supercomputer Polaris, using CUDA-aware MPI to dispatch the solution on
multiple GPUs. We present in §4.1 the stochastic optimal power flow problem,
and give in §4.2 detailed assessments of the algorithms we have introduced
earlier in §3. Eventually, we present in §4.3 a benchmark comparing our
parallel solution algorithm with a state-of-the-art solution method running on
the CPU.

###  4.1 Settings

####  4.1.1 Case study: the block-structured optimal power flow

The stochastic optimal power flow problem aims at finding an optimal dispatch
for the generators $u$. The solution $u$ should minimize the operational costs
while satisfying the physical constraints (power flow equations $g(x,u)=0$,
here playing the role of the state equations) and operational constraints
(line flow constraints $h(x,u)\leq 0$) on a given set of scenarios. Each
scenario is assigned given load parameters (energy demands) and potential
contingencies (line tripping). The values of the state $x$ depend on the local
scenario we are in, the state $x$ being the _recourse_ variable in our case.
As such, the problem has a partially separable structure as introduced in
Problem (2), the control $u$ being shared across all scenarios. We refer to
[7] for the original presentation of the stochastic optimal power flow problem
and to [26, 8, 23, 24] for practical algorithms solving the stochastic optimal
power flow problem (some also focus on the multistage setting, which is not
covered in this article). For our benchmark, we look at reference instances
provided by MATPOWER [47], whose characteristics are detailed in Table 1. We
recall that in our case, the size of the Schur complement matrix
$\hat{K}_{uu}$ is given by the number of controls $n_{u}$.

Name | #bus | #lines | #gen | $n_{x}$ | $n_{u}$
---|---|---|---|---|---
case118 | 118 | 186 | 54 | 181 | 107
case1354pegase | 1,354 | 1,991 | 260 | 2,447 | 519
case2869pegase | 2,869 | 4,582 | 510 | 5,227 | 1,019
case9241pegase | 9,241 | 16,049 | 1,445 | 17,036 | 2,889
Table 1: MATPOWER instances used in the benchmark.

####  4.1.2 Implementation

The algorithm has been implemented entirely in Julia 1.8. The Schur complement
approach has been developed as an extension of the nonlinear optimization
solver MadNLP [41], using CUDA-aware MPI as provided in [6]. We have used the
package ExaPF as a nonlinear modeler for the optimal power flow problem. All
the results presented here have been generated on the supercomputer Polaris
equipped with a total of 560 nodes, each node having with 1 CPU and 4 A100
GPUs.

###  4.2 Assessment of the parallel implementation

####  4.2.1 Assessing the performance of the parallel automatic differentation

We first assess the performance of the parallel automatic differentiation we
introduced in §3.1 in a multi-GPU setting. We compare the performance we
obtain with a CPU implementation. We use case1354pegase as a representative
instance, and display the time spent in the automatic differentiation as we
increase the total number of scenarios $N$. The results are displayed in
Figure 4.

We observe that the computation time depends linearly on the number of
scenarios, as expected. For $N=8$, it is not worthwhile dispatching the
evaluation on multiple GPUs as the problem is small enough to be evaluated on
a single GPU. For $N=512$, the evaluation time is 12.3s on the CPU, compared
to 0.50, 0.41, 0.31, and 0.28s using 1, 2, 4 and 8 GPUs, respectively. Hence,
we get a 40x speed-up when evaluating the derivatives in a multi-GPU setting,
and it is not worthwhile to use more than 4 GPUs (one node).

Figure 4: Time spent to evaluate the model and its derivatives with automatic
differentiation.

####  4.2.2 Assessing the performance of the parallel KKT solver

We proceed to the same performance analysis to assess the performance of the
parallel KKT solver detailed in §3.2. We compare the time required to evaluate
the full solution of the KKT system afresh (including reduction time,
factorization time and backsolve time) on case1354pegase as we increase the
number of scenarios $N$. As a reference, we give the time taken by the sparse
linear solvers HSL MA27 (single-threaded) and HSL MA57 (multi-threaded). The
results are displayed in Figure 5.

On the left, we display the evolution of the time spent in the linear solver
as we increase the number of scenarios. For $N=512$, we observe that we get a
linear speed-up as we increase the number of GPUs: using 8 GPUs, the parallel
KKT solver is 40x faster than using HSL MA27 on the CPU. Interestingly, we
observe that HSL MA57 is not faster than HSL MA27, despite being
multithreaded. This is consistent with the observation made in [42], and
illustrates the difficulty of parallelizing effectively the sparse LDL
factorization (Bunch-Kaufman). On the right, we display a performance profile
detailing the time spent in MA27 and the parallel KKT solver on case1354pegase
with $N=512$ scenarios. We observe that most of the time in HSL MA27 is spent
on factorizing the sparse augmented KKT system (6). On the other side, the
factorization of the dense reduced matrix $\hat{K}_{uu}$ is trivial using
LAPACK on the GPU; the bottleneck in the parallel KKT solver is the reduction
algorithm itself. Fortunately, the reduction algorithm can run in parallel: we
get a linear speed-up as we increase the number of GPUs used in the reduction
algorithm.

Figure 5: Time spent to solve the KKT system for case1354pegase.

####  4.2.3 Assessing the memory consumption

We have observed in §3.1 that the total memory required to store the duals is
$\mathcal{O}((2n_{x}+n_{d}+n_{b})\times M\times p)$, with $M$ being the number
of scenarios stored locally ($M=N$ on 1 GPU, $M=N/2$ on 2 GPUs) and $p$ the
number of tangents. We display in Table 2 the memory taken by the automatic
differentiation backend and by the parallel KKT solver for case1354pegase as
we increase the number of scenarios $N$. We note that storing the duals is
expensive in terms of memory, with up to 10.9GB for $N=512$ on one GPU (as a
reference, each NVIDIA A100 GPU on Polaris has 40GB of memory available). By
evaluating the model on different processes with MPI, we can split the memory
consumption on the different GPUs we are using, leading to better use of the
resource at our disposal.

| 1 GPU | 2 GPUs
---|---|---
$N$ | AD | KKT solver | AD | KKT solver
8 | 171.1 | 92.3 | 85.5 | 48.1
16 | 342.2 | 181.5 | 171.1 | 93.1
32 | 684.3 | 360.0 | 342.2 | 183.2
64 | 1,368.7 | 716.8 | 684.3 | 363.2
128 | 2,737.3 | 1,430.5 | 1,368.7 | 723.4
256 | 5,474.7 | 2,858.0 | 2,737.3 | 1,443.6
512 | 10,949.3 | 5,712.8 | 5,474.7 | 2,884.1
Table 2: Memory consumption in MB

###  4.3 Parallel solution of the block-structured OPF problem

We analyze the parallel performance of our implementation on block-structured
OPF problems.

####  4.3.1 Assessing the parallel performance w.r.t. the number of scenarios

First, we are interested in the scaling of the parallel algorithm in relation
to the total number of scenarios $N$. We consider the case118 instance, and
increase the number of scenarios $N$ from 8 up to 2,048. For each $N$, we
solve the block-structured OPF problem with MadNLP using our parallel KKT
solver, and we compare with the performance we obtained with HSL MA27. The
results are displayed in Figure 6. We observe that the solver HSL MA27 is
initially faster than our parallel KKT solver, as the problem is too small to
benefit from parallelism. However, as soon as $N\geq 16$ the parallel KKT
solver becomes competitive with HSL MA27. The relative performance is
improving as we increase the number of scenarios $N$: for $N=512$, we get a
68x speed-up when using 8 GPUs, compared to the reference HSL MA27 (10.4s
versus 712s). Interestingly, using 2 nodes (=8 GPUs) does not lead to any
speed-up compared to a single node (=4 GPUs) if $N\leq 256$; this setting is
attractive only when the size of the problem becomes sufficiently large
($N\geq 1024$) to compensate for the additional memory exchange.

Figure 6: Time to solve the block-structured OPF problem case118 as a function
of the number of scenarios $N$.

####  4.3.2 Assessing the parallel performance w.r.t. the size of the problem

Second, we increase the size of the problems. We set a fixed number of
scenarios $N=8$, and look at the time to solution for case1354pegase,
case2869pegase and case9241pegase. We detail the respective dimension of each
problem in Table 3. We display the results in Figure 7, and give the detailed
benchmark in Table 4. On the left (a), we display the total time required to
find the solution of the three instances as a function of the number of GPUs;
on the right (b), we show the performance profile associated to
case9241pegase. In (a), we observe that overall the parallel algorithm is
faster than the CPU implementation. The parallel algorithm scales well as we
increase the number of GPUs we are using, the parallel algorithm being 35x
faster than the reference when using 8 GPUs to solve case9241pegase. In (b),
we detail the time spent in the different operations for case9241pegase: the
time spent to factorize the Schur complement with Lapack (using cusolve) is
constant as the size of the Schur complement remains the same as we increase
the number of GPUs. We observe that the time spent in the AD decreases
linearly with the number of GPUs exploited, but the relative time spent in AD
is negligible (less than 5% of the total time). Most of the time is spent in
the parallel reduction, as discussed earlier in §4.2.2.

| $N$ | nvar | ncon |  $\hat{K}_{uu}$ (mb)
---|---|---|---|---
1354pegase | 8 | 20,095 | 53,520 | 2.1
2869pegase | 8 | 42,835 | 119,216 | 7.9
9241pegase | 8 | 139,177 | 404,640 | 63.7
1354pegase | 512 | 1,253,383 | 4,425,280 | 2.1
Table 3: Dimension of the instances we have used in our benchmark. Figure 7:
For a fixed number of scenarios $N=8$, (a) total time spent solving the block-
OPF case1354pegase, case2869pegase and case9241pegase with MadNLP (b)
performance profile for case9241pegase with varying number of GPUs.

| 1354pegase | 2869pegase | 9241pegase
---|---|---|---
| #it | AD | KKT | Tot. | #it | AD | KKT | Tot. | #it | AD | KKT | Tot.
CPU | 44 | 2.6 | 4.2 | 7.0 | 77 | 11.9 | 27.4 | 40.3 | 136 | 205.6 | 771.8 | 984.1
1 GPU | 44 | 0.3 | 1.8 | 2.1 | 93 | 1.1 | 11.7 | 12.8 | 98 | 5.5 | 112.3 | 117.8
2 GPUs | 44 | 0.3 | 1.1 | 1.4 | 93 | 0.8 | 7.4 | 8.2 | 98 | 3.4 | 56.8 | 60.2
4 GPUs | 44 | 0.3 | 1.0 | 1.3 | 93 | 0.8 | 5.7 | 6.5 | 98 | 2.3 | 35.8 | 38.1
8 GPUs | 44 | 0.2 | 1.0 | 1.2 | 93 | 0.6 | 5.1 | 5.7 | 98 | 1.4 | 26.4 | 27.7

Table 4: Detailed results

####  4.3.3 Assessing the parallel performance on a very large-scale instance

We finish our numerical experiments by solving a very large-scale instance:
case1354pegase with $N=512$ scenarios. The dimension of the resulting
optimization problem is displayed in Table 3: the problem has more than 1
million variables, and 4 millions constraints. We solve this instance on resp.
1 node, 2, 4 and 8 nodes (resp. 4, 8, 16 and 32 GPUs). The results are
displayed in Figure 8. We observe that the scaling is almost perfect when we
use 2 nodes (8 GPUs) instead of a single node (4 GPUs) but we do not observe
the same behavior when we increase the number of nodes to 4 and 8. On that
instance, the gain we get when using 8 nodes (32 GPUs) is marginal compared to
when using 4 nodes (16 GPUs): the solving time only decreases from 67s to 58s.
This corroborate our observations: it is better to pack all the computation on
a single node to use four A100 GPUs connected together via unified memory
(NVLINK has a transfer rate of 600GB/s). When we have to use more than 2
nodes, the memory transfers are more involved as they have to pass through the
network of the supercomputer.

Figure 8: Solving case1354pegase with $N=512$

##  5 Conclusion

We show promising results for leveraging massively parallel SIMD architectures
like GPUs for block-structured nonlinear programs. The parallelism is applied
to both the derivative evaluation and the solution of the KKT linear system.
The main operation in the KKT algorithm is the assembling of the Schur
complement, the factorization of the dense Schur complement being fast to
carry on the GPU.

At all levels, the method benefits significantly from the massive parallelism,
achieving a speedup of around 40 for the derivatives compared to a sequential
CPU implementation. The speedup is very application dependent, not least on
the Hessian coloring and the problem’s structure. The assembling of the Schur
complement is bottlenecked by a distributed reduction operation bound by the
interconnect’s latency and throughput between GPUs. Current, so-called super
nodes with multiple GPUs connected via fast networks like NVLINK greatly
accelerate this operation. Lastly, our method is limited by the memory
capacity of the GPU accelerators as it grows linearly with the number of
problem blocks. In the context of ACOPF we are confident that upcoming GPUs
will provide enough memory to solve a large number of scenarios in parallel,
even for the largest grid instances (e.g., Eastern Interconnection with 70,000
nodes).

With the upcoming release of the Aurora supercomputer, these SIMD
architectures will allow new science in regimes that were impossible with
previous CPU architectures.

## Acknowledgment

This material was based upon work supported by the U.S. Department of Energy,
Office of Science, Office of Advanced Scientific Computing Research (ASCR)
under Contract DE-AC02-06CH11347 and by NSF through award CNS-1545046. The
authors gratefully acknowledge the funding support from the Applied
Mathematics Program within the U.S. Department of Energy’s (DOE) Office of
Advanced Scientific Computing Research (ASCR) as part of the project ExaSGD.
This research used resources of the Argonne Leadership Computing Facility,
which is a DOE Office of Science User Facility supported under Contract DE-
AC02-06CH11357.

## References

  * Amestoy et al., [2000] Amestoy, P. R., Duff, I. S., and L’excellent, J.-Y. (2000).  Multifrontal parallel distributed symmetric and unsymmetric solvers.  Computer methods in applied mechanics and engineering, 184(2-4):501–520.
  * Birge and Qi, [1988] Birge, J. R. and Qi, L. (1988).  Computing block-angular Karmarkar projections with applications to stochastic programming.  Management science, 34(12):1472–1479.
  * Bradbury et al., [2018] Bradbury, J., Frostig, R., Hawkins, P., Johnson, M. J., Leary, C., Maclaurin, D., Necula, G., Paszke, A., VanderPlas, J., Wanderman-Milne, S., and Zhang, Q. (2018).  JAX: composable transformations of Python+NumPy programs.
  * Bücker et al., [2001] Bücker, H. M., Lang, B., an Mey, D., and Bischof, C. H. (2001).  Bringing together automatic differentiation and OpenMP.  In Proceedings of the 15th international conference on Supercomputing, pages 246–251.
  * Bussieck and Meeraus, [2004] Bussieck, M. R. and Meeraus, A. (2004).  General algebraic modeling system (GAMS).  In Modeling languages in mathematical optimization, pages 137–157. Springer.
  * Byrne et al., [2021] Byrne, S., Wilcox, L. C., and Churavy, V. (2021).  MPI. jl: Julia bindings for the Message Passing Interface.  In Proceedings of the JuliaCon Conferences, volume 1, page 68.
  * Capitanescu et al., [2011] Capitanescu, F., Ramos, J. M., Panciatici, P., Kirschen, D., Marcolini, A. M., Platbrood, L., and Wehenkel, L. (2011).  State-of-the-art, challenges, and future trends in security constrained optimal power flow.  Electric power systems research, 81(8):1731–1741.
  * Chiang et al., [2014] Chiang, N., Petra, C. G., and Zavala, V. M. (2014).  Structured nonconvex optimization of large-scale energy systems using PIPS-NLP.  In 2014 Power Systems Computation Conference, pages 1–7. IEEE.
  * Choi and Goldfarb, [1993] Choi, I. C. and Goldfarb, D. (1993).  Exploiting special structure in a primal—dual path-following algorithm.  Mathematical Programming, 58(1):33–52.
  * Colombo et al., [2009] Colombo, M., Grothey, A., Hogg, J., Woodsend, K., and Gondzio, J. (2009).  A structure-conveying modelling language for mathematical and stochastic programming.  Mathematical Programming Computation, 1(4):223–247.
  * DeMiguel and Nogales, [2008] DeMiguel, V. and Nogales, F. J. (2008).  On decomposition methods for a class of partially separable nonlinear programs.  Mathematics of Operations Research, 33(1):119–139.
  * Duff, [2004] Duff, I. S. (2004).  MA57—a code for the solution of sparse symmetric definite and indefinite systems.  ACM Transactions on Mathematical Software (TOMS), 30(2):118–144.
  * Duff and Van Der Vorst, [1999] Duff, I. S. and Van Der Vorst, H. A. (1999).  Developments and trends in the parallel solution of linear systems.  Parallel Computing, 25(13-14):1931–1970.
  * Dunning et al., [2017] Dunning, I., Huchette, J., and Lubin, M. (2017).  JuMP: A modeling language for mathematical optimization.  SIAM review, 59(2):295–320.
  * Fourer et al., [1990] Fourer, R., Gay, D. M., and Kernighan, B. W. (1990).  A modeling language for mathematical programming.  Management Science, 36(5):519–554.
  * Gondzio and Grothey, [2009] Gondzio, J. and Grothey, A. (2009).  Exploiting structure in parallel implementation of interior point methods for optimization.  Computational Management Science, 6(2):135–160.
  * Gondzio and Sarkissian, [2003] Gondzio, J. and Sarkissian, R. (2003).  Parallel interior-point solver for structured linear programs.  Mathematical Programming, 96(3):561–584.
  * Griewank and Walther, [2008] Griewank, A. and Walther, A. (2008).  Evaluating derivatives: principles and techniques of algorithmic differentiation.  SIAM.
  * Hovland and Bischof, [1998] Hovland, P. and Bischof, C. (1998).  Automatic differentiation for message-passing parallel programs.  In Proceedings of the First Merged International Parallel Processing Symposium and Symposium on Parallel and Distributed Processing, pages 98–104. IEEE.
  * Hovland, [1997] Hovland, P. D. (1997).  Automatic differentiation of parallel programs.  University of Illinois at Urbana-Champaign.
  * Huchette et al., [2014] Huchette, J., Lubin, M., and Petra, C. (2014).  Parallel algebraic modeling for stochastic optimization.  In 2014 First Workshop for High Performance Technical Computing in Dynamic Languages, pages 29–35. IEEE.
  * Jessup et al., [1994] Jessup, E. R., Yang, D., and Zenios, S. A. (1994).  Parallel factorization of structured matrices arising in stochastic programming.  SIAM journal on Optimization, 4(4):833–846.
  * Kardoš et al., [2019] Kardoš, J., Kourounis, D., and Schenk, O. (2019).  Two-level parallel augmented Schur complement interior-point algorithms for the solution of security constrained optimal power flow problems.  IEEE Transactions on power systems, 35(2):1340–1350.
  * Kardoš et al., [2020] Kardoš, J., Kourounis, D., and Schenk, O. (2020).  Structure-exploiting interior point methods.  In Parallel Algorithms in Computational Science and Engineering, pages 63–93. Springer.
  * Kardoš et al., [2022] Kardoš, J., Kourounis, D., Schenk, O., and Zimmerman, R. (2022).  Beltistos: A robust interior point method for large-scale optimal power flow problems.  Electric Power Systems Research, 212:108613.
  * Lubin et al., [2011] Lubin, M., Petra, C. G., Anitescu, M., and Zavala, V. (2011).  Scalable stochastic optimization of complex energy systems.  In SC’11: Proceedings of 2011 International Conference for High Performance Computing, Networking, Storage and Analysis, pages 1–10. IEEE.
  * Moses et al., [2021] Moses, W. S., Churavy, V., Paehler, L., Hückelheim, J., Narayanan, S. H. K., Schanen, M., and Doerfert, J. (2021).  Reverse-mode automatic differentiation and optimization of GPU kernels via Enzyme.  In Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis, pages 1–16.
  * Nash and Sofer, [1989] Nash, S. G. and Sofer, A. (1989).  Block truncated-Newton methods for parallel optimization.  Mathematical Programming, 45(1):529–546.
  * Nash and Sofer, [1991] Nash, S. G. and Sofer, A. (1991).  A general-purpose parallel algorithm for unconstrained optimization.  SIAM Journal on Optimization, 1(4):530–547.
  * Nocedal and Wright, [2006] Nocedal, J. and Wright, S. J. (2006).  Numerical optimization.  Springer series in operations research. Springer, New York, 2nd edition.
  * Pacaud et al., [2022] Pacaud, F., Shin, S., Schanen, M., Maldonado, D. A., and Anitescu, M. (2022).  Accelerating condensed interior-point methods on SIMD/GPU architectures.  arXiv preprint arXiv:2203.11875.
  * Paszke et al., [2019] Paszke, A., Gross, S., Massa, F., Lerer, A., Bradbury, J., Chanan, G., Killeen, T., Lin, Z., Gimelshein, N., Antiga, L., et al. (2019).  Pytorch: An imperative style, high-performance deep learning library.  Advances in neural information processing systems, 32.
  * Petra et al., [2014] Petra, C. G., Schenk, O., Lubin, M., and Gärtner, K. (2014).  An augmented incomplete factorization approach for computing the Schur complement in stochastic optimization.  SIAM Journal on Scientific Computing, 36(2):C139–C162.
  * Rodriguez et al., [2021] Rodriguez, J. S., Parker, R., Laird, C. D., Nicholson, B., Siirola, J. D., and Bynum, M. (2021).  Scalable parallel nonlinear optimization with PyNumero and Parapint.  Optimization Online.
  * Ruszczynski, [1993] Ruszczynski, A. (1993).  Interior point methods in stochastic programming.  Working paper.
  * Saad, [1980] Saad, Y. (1980).  On the rates of convergence of the Lanczos and the block-Lanczos methods.  SIAM Journal on Numerical Analysis, 17(5):687–706.
  * Schanen et al., [2018] Schanen, M., Gilbert, F., Petra, C. G., and Anitescu, M. (2018).  Toward multiperiod ac-based contingency constrained optimal power flow at large scale.  In 2018 Power Systems Computation Conference (PSCC), pages 1–7. IEEE.
  * Schenk and Gärtner, [2004] Schenk, O. and Gärtner, K. (2004).  Solving unsymmetric sparse systems of linear equations with PARDISO.  Future Generation Computer Systems, 20(3):475–487.
  * Schnabel, [1985] Schnabel, R. B. (1985).  Parallel computing in optimization.  In Computational Mathematical Programming, pages 357–381. Springer.
  * Schnabel, [1995] Schnabel, R. B. (1995).  A view of the limitations, opportunities, and challenges in parallel nonlinear optimization.  Parallel computing, 21(6):875–905.
  * Shin et al., [2021] Shin, S., Coffrin, C., Sundar, K., and Zavala, V. M. (2021).  Graph-based modeling and decomposition of energy infrastructures.  IFAC-PapersOnLine, 54(3):693–698.
  * Tasseff et al., [2019] Tasseff, B., Coffrin, C., Wächter, A., and Laird, C. (2019).  Exploring benefits of linear solver parallelism on modern nonlinear optimization applications.  arXiv preprint arXiv:1909.08104.
  * Watson et al., [2012] Watson, J.-P., Woodruff, D. L., and Hart, W. E. (2012).  PySP: modeling and solving stochastic programs in python.  Mathematical Programming Computation, 4(2):109–149.
  * Word et al., [2014] Word, D. P., Kang, J., Akesson, J., and Laird, C. D. (2014).  Efficient parallel solution of large-scale nonlinear dynamic optimization problems.  Computational Optimization and Applications, 59(3):667–688.
  * Zavala et al., [2008] Zavala, V. M., Laird, C. D., and Biegler, L. T. (2008).  Interior-point decomposition approaches for parallel solution of large-scale nonlinear parameter estimation problems.  Chemical Engineering Science, 63(19):4834–4845.
  * Zhu et al., [2009] Zhu, Y., Word, D., Siirola, J., and Laird, C. D. (2009).  Exploiting modern computing architectures for efficient large-scale nonlinear programming.  In Computer Aided Chemical Engineering, volume 27, pages 783–788. Elsevier.
  * Zimmerman et al., [2010] Zimmerman, R. D., Murillo-Sánchez, C. E., and Thomas, R. J. (2010).  MATPOWER: Steady-state operations, planning, and analysis tools for power systems research and education.  IEEE Transactions on Power Systems, 26(1):12–19.

Government License: The submitted manuscript has been created by UChicago
Argonne, LLC, Operator of Argonne National Laboratory (“Argonne”). Argonne, a
U.S. Department of Energy Office of Science laboratory, is operated under
Contract No. DE-AC02-06CH11357. The U.S. Government retains for itself, and
others acting on its behalf, a paid-up nonexclusive, irrevocable worldwide
license in said article to reproduce, prepare derivative works, distribute
copies to the public, and perform publicly and display publicly, by or on
behalf of the Government. The Department of Energy will provide public access
to these results of federally sponsored research in accordance with the DOE
Public Access Plan. http://energy.gov/downloads/doe-public-access-plan.


[a]Gregory Foote

# Exploring the Potential of the Pulsed Laser onboard the CALIPSO Satellite to
Improve Calibration with VERITAS

###### Abstract

Imaging Atmospheric Cherenkov Telescopes (IACTs) are used to detect bright
nanosecond-duration flashes of optical light originating from interactions of
cosmic/gamma-rays in the atmosphere. A natural calibration source with similar
characteristics does not exist; however, satellite-based laser systems provide
a potential alternative. The CALIPSO satellite is one such facility which uses
a suite of instruments to gather information about the atmosphere. Of
particular interest is the CALIOP instrument, which emits 20-nanosecond laser
pulses at 1064 nm and 532 nm at a rate of 20 Hz towards the Earth. The TAIGA-
HiSCORE collaboration announced a detection of CALIOP laser pulses at the 37th
ICRC in 2021, demonstrating that the laser footprint extends to at least tens
of kilometers from the subsatellite point. We have used the VERITAS IACT to
observe CALIPSO, and show here the results of using these observations to help
to calibrate the array. We also discuss the potential of this technique for
cross-calibration between different IACT facilities and for relative
calibration between the telescopes of future large arrays.

##  1 Introduction

The Very Energetic Radiation Imaging Telescope Array System (VERITAS) is a set
of four imaging atmospheric Cherenkov telescopes located at the Fred Lawrence
Whipple Observatory (FLWO) in southern Arizona (31$\deg$ 40′N, 110$\deg$ 57′W,
1.3km a.s.l.), and is designed to observe the phenomenon of “Cherenkov
showers” where large bursts of Cherenkov light are produced high in the
atmosphere originating from the interactions between high energy gamma
rays/cosmic rays and the Earth’s atmosphere. Each telescope is a 12 m Davies-
Cotton telescope, with all four telescopes roughly 100 m apart in a
rectangular layout. Each 12 m reflector is made of 350 hexagonal mirrors which
focus light onto a camera composed of an array of 499 photo-multiplier tubes
(PMTs). Each PMT is fitted with a hexagonal light cone in order to remove the
gaps between PMTs to tile the camera plane [1]. In order to prevent being
overwhelmed with data and to only collect images of Cherenkov showers, VERITAS
employs a three-tiered trigger system. This trigger system begins with a
pixel-level trigger associated with the signal measured by a single PMT
exceeding a set threshold. The next level up is a telescope-level trigger
associated with a group of three or more neighboring pixel triggers within 5
ns of each other in the same camera. Finally, an array-level trigger is formed
if two or more telescope-level triggers enter a triggered state within 50 ns
of each other after correcting for delays of the Cherenkov light between
telescopes. Triggered events are read out using flash-ADCs which record 32-ns
signal traces for each PMT with 2-ns samples [2].

The Cloud-Aerosol LIDAR and Infrared Pathfinder Satellite Observations
(CALIPSO) satellite [3] is a joint atmospheric science mission built,
maintained, and launched by NASA and CNES, which operates at an altitude of
700 km. Its primary instrument, the Cloud-Aerosol LIdar with Orthogonal
Polarization (CALIOP) instrument, is of particular interest to this study
because it hosts two lasers: a 532 nm laser and 1064 nm laser, both of which
produce 20-ns pulses at a rate of 20.16 Hz directed towards the Earth . While
532 nm is within the detectable wavelength range of Cherenkov telescopes, the
footprint on the ground was expected to be too small (on the order of 100 m)
to be observable from anywhere but the direct satellite path. However, the
TAIGA-HiScore telescope announced in 2021 that they were able to detect
CALIPSO [4], which refuted the assumption about the footprint and spurred
VERITAS to begin observations. Using VERITAS, CALIPSO has been shown to be
consistently detectable. We have conducted an observation campaign of CALIPSO
passes, as well as developed analysis methods to exploit its potential as a
calibration source for Cherenkov telescopes. Potentially visible CALIPSO
passes occur approximately every three days and, as a stable, nanosecond-
pulsed light source located at large distance, its potential as a calibration
source is unique.

##  2 VERITAS Observations of CALIPSO

The observation campaign of CALIPSO with VERITAS covers a wide range of
azimuth angles, and elevation angles between zenith and $55\deg$ elevation. As
the dataset of CALIPSO observations grows, we are developing an understanding
of the correlation between the elevation and azimuth of the satellite and
VERITAS’s response. This allows us to predict the response at a given sky
location. This also provides an opportunity to study the laser footprint of
CALIPSO which could potentially help other observatories utilize CALIPSO for
internal or cross-calibrations, including with VERITAS.

We identify CALIPSO pulses in our data by exploiting two characteristics: the
lack of parallax and the pulse arrival time. First, we require the centers of
the images from each telescope (within a single event) to be within two pixel
widths ($0.3\deg$) of each other; as seen in figure 1. We also calculate the
phase of the pulse arrival time with respect to the 20.16 Hz laser pulse
frequency. If a large collection of events with phases within 0.007 of each
other appears in the observation run, they are classed as CALIPSO pulses. Once
the pulses are identified in our data, we perform analysis on them to get
calibration products.

Figure 1: Two events as seen by VERITAS are showcased, the top row produced by
CALIPSO, and the bottom row produced by a regular cosmic-ray, a) shows a
CALIPSO pulse as seen by each telescope, b) shows digital counts versus time
from the bright pixel in telescope 4 for this pulse, c) shows a high-altitude
cosmic-ray as seen by each telescope, and d) shows digital counts versus time
from the bright pixel in telescope 3 for this cosmic-ray. The typical bounds
of the pulse sum are shown with the blue boxes in c) and d).

##  3 Camera Rotations

During the construction of VERITAS, the cameras were placed on their support
structures with a slight rotation accidentally. It is therefore necessary to
measure and correct for this effect in our analysis. Initially this rotation
was measured using a plumb bob, but the measurement became simpler once CCD
cameras were installed which record images of both the PMT camera of each
telescope and the sky in the direction to which they are pointed [5]. The
current method positions the telescopes at $180\deg$ azimuth and angled
towards the Celestial Equator, such that a bright star travels through the
middle of all the cameras. When the star travels through the field-of-view,
the track across the middle of the camera is compared between the CCD camera’s
view and each individual telescope to determine the amount of rotation which
is needed to correct for the camera rotations. CALIPSO pulses can be used to
perform a similar measurement, because we know that the position and track of
CALIPSO should be the same in every camera. We can also potentially measure
changes in camera rotations at different elevations and across time, if such
changes exist.

Using data from an observation of CALIPSO taken on January 29${}^{\text{th}}$,
2023 at $58\deg$ elevation, we calculated the camera rotations by first
solving for the angle between the positive X-axis and CALIPSO’s path in each
camera. This angle is found by first performing a least-squares fit of a line
on the path CALIPSO makes in the camera, as seen in figure 2. We then use
equation 1 to find the angle and equation 2 to find its error.

$\theta_{i}=\arctan\left(m_{i\text{,fit}}\right)$ (1)

$\delta\theta_{i}=\left|\arctan\left(m_{i\text{,fit}}\right)\right|\cdot\left|\frac{\delta
m_{i\text{,fit}}}{m^{3}_{i\text{,fit}}+m_{i\text{,fit}}}\right|$ (2)

The relative rotation of each telescope is solved for by setting telescope 1
to zero. Doing this provides a comparison to the currently used values within
VERITAS, which used the traditional method. The comparison of this analysis
and the currently used values is found in Table 1.

Figure 2: The fitted path through the VERITAS camera (dotted line) along with each image center for every captured event (blue dots). Telescope ID | Current Rotations | This Analysis
---|---|---
1 | $0.$ | $0.\pm 0.077$
2 | $0.612$ | $0.692\pm 0.092$
3 | $-0.806$ | $-0.426\pm 0.088$
4 | $2.027$ | $2.086\pm 0.079$
Table 1: Table of the measured rotations from this analysis and the rotation
measurements currently used in VERITAS. The rotations have been modified from
their original values to be relative to Telescope 1 for ease of comparison;
all measurements are in degrees.

##  4 Relative Throughputs

The optical systems of VERITAS experience significant change over time. The
telescope mirrors suffer from degradation due to exposure to the southern
Arizona desert and climate throughout the year. These changes are corrected
for with maintenance tasks including mirror cleanings, swaps, replacements and
re-coatings. Alongside this, the PMTs also degrade over time from two effects:
operating at high voltages over long periods of time and charged positive ions
accumulating within them. Measuring these effects is crucial in maintaining an
accurate model of VERITAS within our simulations. A full discussion of current
VERITAS calibration techniques to measure the impacts on the throughputs of a
given telescope (defined as the product of the reflectivity, accounting for
the mirrors, and the gain, accounting for the PMTs) and the results thereof
can be found in [6].

We are exploring the possibility of utilizing CALIPSO for performing some
aspects of this calibration. Performing this calibration with CALIPSO is
possible because the spacing of the VERITAS telescopes relative to the
altitude of CALIPSO, 100 meters versus 700 kilometers, ensures that the amount
of energy received by each telescope should be approximately the same.

The total amount of light captured by the VERITAS camera is estimated to be
proportional to the sum of all traces from PMTs which pass standard image
cleaning procedures. The portion of the trace which is included in this sum
normally begins roughly halfway through the rising edge of the captured pulse
and covers 12 ns [7]. This method, called the double-pass method, works well
for the pulses produced from Cherenkov showers, which are usually entirely
captured within the 32-ns FADC window stored. However, CALIPSO pulses are
relatively slow, with a full-width at half-maximum of 20 ns, making this
method unreliable. The differences in the pulse characteristics are shown in
figure 1. For this analysis, the analysis method was adjusted to begin the
pulse sum at the start of the 32 ns captured window and end at the maximum of
the pulse. While normally this would significantly increase the amount of
noise, the uniformity and intensity of the captured CALIPSO pulses mitigates
this.

The above method yields the amount of light captured per pulse per telescope;
to provide a more useful comparison we find the average of the light intensity
across all pulses for a given telescope in a single transit; along with the
error of that average. We move to a per-transit basis as opposed to working
with individual pulses because across the entire transit the laser does not
significantly change. To get an individual telescope’s relative throughput for
a transit, we divide the response from each telescope by the error-weighted
average of all telescopes for that transit. These relative throughputs can
then be compared across different transits at different elevations. The
results from 7 CALIPSO transits are shown in figure 3. An early verification
of our method is to use it to show an increase of the relative throughput of
Telescope 1 after March 10${}^{\text{th}}$, 2023. Between March
1${}^{\text{st}}$ and 10${}^{\text{th}}$, 91 mirror facets were replaced and
aligned which increased the reflectivity of Telescope 1. While more data is
required to understand how well this method performs relative to traditional
methods, the method shows initial promise.

Figure 3: The relative throughputs for each telescope across time, calculated
using CALIPSO laser pulses. The completion date of mirror facet replacements
and alignment on Telescope 1, March 10${}^{\text{th}}$, 2023, is indicated by
the solid black line.

##  5 Conclusion

CALIPSO represents an unique calibration tool as it guarantees a point-like,
bright, pulsed, distant optical source whose location and energy are known.
This enables several calibrations of VERITAS, two of which are explored here,
with data taken under normal operating conditions. Additional calibration uses
of CALIPSO which we have not yet explored include measurements of the
telescope optical point spread functions and atmospheric model validation, as
well as studies of the stability of the telescope properties over time. As the
dataset of CALIPSO observations with VERITAS grows, our understanding of the
VERITAS response to the laser will improve, which will help with current and
future calibration studies.

##  6 Acknowledgements

This research is supported by grants from the U.S. Department of Energy Office
of Science, the U.S. National Science Foundation and the Smithsonian
Institution, by NSERC in Canada, and by the Helmholtz Association in Germany.
This research used resources provided by the Open Science Grid, which is
supported by the National Science Foundation and the U.S. Department of
Energy’s Office of Science, and resources of the National Energy Research
Scientific Computing Center (NERSC), a U.S. Department of Energy Office of
Science User Facility operated under Contract No. DE-AC02-05CH11231. We
acknowledge the excellent work of the technical support staff at the Fred
Lawrence Whipple Observatory and at the collaborating institutions in the
construction and operation of the instrument.

## References

  * [1] J. Holder, _VERITAS: Status and Highlights_ , in _International Cosmic Ray Conference_ , vol. 11 of _International Cosmic Ray Conference_ , p. 137, Jan., 2011, DOI [1111.1225].
  * [2] A. Weinstein, _The VERITAS Trigger System_ , in _International Cosmic Ray Conference_ , vol. 3 of _International Cosmic Ray Conference_ , pp. 1539–1542, Jan., 2008, DOI [0709.4438].
  * [3] D.M. Winker, M.A. Vaughan, A. Omar, Y. Hu, K.A. Powell, Z. Liu et al., _Overview of the CALIPSO Mission and CALIOP Data Processing Algorithms_ , _JAtOT_ 26 (2009) 2310.
  * [4] A. Porelli and Taiga Collaboration, _Satellite-based Calibration of the TAIGA-HiSCORE Cerenkov Array by the LIDAR on-board CALIPSO_ , in _37th International Cosmic Ray Conference_ , p. 876, Mar., 2022, DOI.
  * [5] J. Holder, V.A. Acciari, E. Aliu, T. Arlen, M. Beilicke, W. Benbow et al., _Status of the VERITAS Observatory_ , in _American Institute of Physics Conference Series_ , F.A. Aharonian, W. Hofmann and F. Rieger, eds., vol. 1085 of _American Institute of Physics Conference Series_ , pp. 657–660, Dec., 2008, DOI [0810.0474].
  * [6] C.B. Adams, W. Benbow, A. Brill, J.H. Buckley, J.L. Christiansen, A. Falcone et al., _The throughput calibration of the VERITAS telescopes_ , _A &A_ 658 (2022) A83 [2111.04676].
  * [7] M.K. Daniel, _The VERITAS standard data analysis_ , in _International Cosmic Ray Conference_ , vol. 3 of _International Cosmic Ray Conference_ , pp. 1325–1328, Jan., 2008, DOI [0709.4006].

## Full Author List: VERITAS Collaboration

A. Acharyya1, C. B. Adams2, A. Archer3, P. Bangale4, J. T. Bartkoske5, P.
Batista6, W. Benbow7, J. L. Christiansen8, A. J. Chromey7, A. Duerr5, M.
Errando9, Q. Feng7, G. M. Foote4, L. Fortson10, A. Furniss11,12, W. Hanlon7,
O. Hervet12, C. E. Hinrichs7,13, J. Hoang12, J. Holder4, Z. Hughes9, T. B.
Humensky14,15, W. Jin1, M. N. Johnson12, M. Kertzman3, M. Kherlakian6, D.
Kieda5, T. K. Kleiner6, N. Korzoun4, S. Kumar14, M. J. Lang16, M. Lundy17, G.
Maier6, C. E McGrath18, M. J. Millard19, C. L. Mooney4, P. Moriarty16, R.
Mukherjee20, S. O’Brien17,21, R. A. Ong22, N. Park23, C. Poggemann8, M.
Pohl24,6, E. Pueschel6, J. Quinn18, P. L. Rabinowitz9, K. Ragan17, P. T.
Reynolds25, D. Ribeiro10, E. Roache7, J. L. Ryan22, I. Sadeh6, L. Saha7, M.
Santander1, G. H. Sembroski26, R. Shang20, M. Splettstoesser12, A. K.
Talluri10, J. V. Tucci27, V. V. Vassiliev22, A. Weinstein28, D. A. Williams12,
S. L. Wong17, and J. Woo29

1Department of Physics and Astronomy, University of Alabama, Tuscaloosa, AL
35487, USA

2Physics Department, Columbia University, New York, NY 10027, USA

3Department of Physics and Astronomy, DePauw University, Greencastle, IN
46135-0037, USA

4Department of Physics and Astronomy and the Bartol Research Institute,
University of Delaware, Newark, DE 19716, USA

5Department of Physics and Astronomy, University of Utah, Salt Lake City, UT
84112, USA

6DESY, Platanenallee 6, 15738 Zeuthen, Germany

7Center for Astrophysics $|$ Harvard & Smithsonian, Cambridge, MA 02138, USA

8Physics Department, California Polytechnic State University, San Luis Obispo,
CA 94307, USA

9Department of Physics, Washington University, St. Louis, MO 63130, USA

10School of Physics and Astronomy, University of Minnesota, Minneapolis, MN
55455, USA

11Department of Physics, California State University - East Bay, Hayward, CA
94542, USA

12Santa Cruz Institute for Particle Physics and Department of Physics,
University of California, Santa Cruz, CA 95064, USA

13Department of Physics and Astronomy, Dartmouth College, 6127 Wilder
Laboratory, Hanover, NH 03755 USA

14Department of Physics, University of Maryland, College Park, MD, USA

15NASA GSFC, Greenbelt, MD 20771, USA

16School of Natural Sciences, University of Galway, University Road, Galway,
H91 TK33, Ireland

17Physics Department, McGill University, Montreal, QC H3A 2T8, Canada

18School of Physics, University College Dublin, Belfield, Dublin 4, Ireland

19Department of Physics and Astronomy, University of Iowa, Van Allen Hall,
Iowa City, IA 52242, USA

20Department of Physics and Astronomy, Barnard College, Columbia University,
NY 10027, USA

21 Arthur B. McDonald Canadian Astroparticle Physics Research Institute, 64
Bader Lane, Queen’s University, Kingston, ON Canada, K7L 3N6

22Department of Physics and Astronomy, University of California, Los Angeles,
CA 90095, USA

23Department of Physics, Engineering Physics and Astronomy, Queen’s
University, Kingston, ON K7L 3N6, Canada

24Institute of Physics and Astronomy, University of Potsdam, 14476 Potsdam-
Golm, Germany

25Department of Physical Sciences, Munster Technological University,
Bishopstown, Cork, T12 P928, Ireland

26Department of Physics and Astronomy, Purdue University, West Lafayette, IN
47907, USA

27Department of Physics, Indiana University-Purdue University Indianapolis,
Indianapolis, IN 46202, USA

28Department of Physics and Astronomy, Iowa State University, Ames, IA 50011,
USA

29Columbia Astrophysics Laboratory, Columbia University, New York, NY 10027,
USA