# DiffGEPCI: 3D MRI Synthesis from mGRE Signals
using 2.5D Diffusion Model

Yuyang Hu, Satya V. V. N. Kothapalli, Weijie Gan,
Alexander L. Sukstanskii, Gregory F. Wu, Manu Goyal,
Dmitriy A. Yablonskiy, Ulugbek S. Kamilov
Washington University in St.Louis

###### Abstract

We introduce a new framework called DiffGEPCI for cross-modality generation in
magnetic resonance imaging (MRI) using a 2.5D conditional diffusion model.
DiffGEPCI can synthesize high-quality Fluid Attenuated Inversion Recovery
(FLAIR) and Magnetization Prepared-Rapid Gradient Echo (MPRAGE) images,
without acquiring corresponding measurements, by leveraging multi-Gradient-
Recalled Echo (mGRE) MRI signals as conditional inputs. DiffGEPCI operates in
a two-step fashion: it initially estimates a 3D volume slice-by-slice using
the axial plane and subsequently applies a refinement algorithm (referred to
as 2.5D) to enhance the quality of the coronal and sagittal planes.
Experimental validation on real mGRE data shows that DiffGEPCI achieves
excellent performance, surpassing generative adversarial networks (GANs) and
traditional diffusion models.

##  1 Introduction

Magnetic Resonance Imaging (MRI) is a versatile medical imaging technique that
can produce multiple image contrasts that provide complementary information
about a tissue. There is a growing interest in MRI to synthesize multiple
naturally co-registered images corresponding to different MRI contrasts from a
single scan [1]. Gradient Echo Plural Contrast Imaging (GEPCI) [2] has been
previously proposed for generating FLAIR [3], MPRAGE [4], and other commonly
used contrasts images from a single Multi-Gradient-Recalled Echo (mGRE)
sequence data [5]. However, due to its purely physics-driven nature, original
GEPCI algorithms did not provide perceptually high-quality images similar to
those used in clinic. We seek to address this issue by proposing a new data-
driven method for generating FLAIR and MPRAGE images with high-perceptual
quality from mGRE.

Generative models are currently popular for image-to-image translation [6].
Conditional generative adversarial networks (cGAN) have been shown to create
realistic images [7, 8, 9]. The training of cGANs relies on adversarial
learning, where a discriminator captures target distribution information and
guides a mapping from source to target. Conditional diffusion models (CDMs)
have recently surpassed cGANs in terms of ability to generate high-quality and
realistic images [10, 11, 12, 13]. CDMs start with noise and progressively
refine output through iterative denoising while also integrating conditional
information. This advancement has enabled CDMs to achieve state-of-the-art
results in multi-contrast MR synthesis [14]. In the context of 3D multi-
contrast MRI synthesis, current diffusion models are limited to generating
images slice-by-slice. This approach can potentially compromise the
preservation of fine details and may introduce artifacts along the other two
dimensions. Directly implementing a 3D diffusion model to address this
challenge may not be feasible for two primary reasons: (a) it necessitates a
significantly larger dataset (comprising thousands of volumes as opposed to
thousands of slices) and (b) requires higher computational cost.

We introduce the Denoising Diffusion Model for Gradient Echo Plural Contrast
Imaging (DiffGEPCI) method that incorporates a 2.5D processing into diffusion-
model-based MRI synthesis. DiffGEPCI uses a conditional diffusion model
trained with mGRE signals as conditional input, sourced from three distinct
planes: Axial, Coronal, and Sagittal. During inference, DiffGEPCI generates a
primary translated volume from the Axial plane and incorporates a specialized
2.5D refinement module to improve performance across all three planes.
DiffGEPCI consistently yields state-of-the-art results for both mGRE
$\rightarrow$ FLAIR and mGRE $\rightarrow$ MPRAGE synthesis.

Figure 1: DiffGEPCI relies on a two-step process for enhancing the quality of
MRI synthesis. In the first step (a), a initial target modality volume is
generated from the axial plane. In the second step (b), this volume is split
into two additional plane volumes, and a 10-step inverse diffusion process is
applied to improve their quality. Finally, the three volumes are fused by
averaging their pixel values, leading to enhanced performance in both coronal
and sagittal planes.  Figure 2: Performance Comparison of several baselines
and our proposed method in Aixal planes. The top row illustrates the visual
performance mGRE $\rightarrow$ FLAIR task, while the bottom row presents mGRE
$\rightarrow$ MPRAGE. Our method estimates intricate structures while avoiding
troublesome artifacts, leading to the best PSNR and SSIM performance.

##  2 Proposed Method

###  2.1 Conditional Denoising Diffusion Probabilistic Model

Conditional Denoising Diffusion Probabilistic Model (cDDPM) is designed to
learn a parameterized Markov chain to convert Gaussian distribution to
conditional data distribution $p({\bm{x}},{\bm{y}})$, see [11, 9] for review
of DDPM. In this paper, ${\bm{x}}$ refers to our target modalities and the
condition ${\bm{y}}$ refers to the 10 mGRE signals (see [2] for information on
mGRE). The forward diffusion process starts from clean sample ${\bm{x}}_{0}$
and then gradually adds Gaussian noise to ${\bm{x}}_{0}$ according to the
following transition probability:

$q({\bm{x}}_{t}|{\bm{x}}_{t-1}):={\mathcal{N}}{({\bm{x}}_{t};\sqrt{1-\beta_{t}}{\bm{x}}_{t-1},\beta_{t}{\mathbf{I}}}),$
(1)

where $t\in[0,T]$, ${\mathcal{N}}(\cdot)$ denotes the Gaussian probability
density function (pdf), $\beta_{t}$ refers to the noise schedule of the
process. By parameter changes $\alpha_{t}:=1-\beta_{t}$ and
$\bar{\alpha}_{t}:=\prod^{t}_{s=1}\alpha_{s}$, we can rewrite ${\bm{x}}_{t}$
as a linear combination of noise $\bm{\epsilon}$ and ${\bm{x}}_{0}$

${\bm{x}}_{t}=\sqrt{\bar{\alpha}_{t}}{\bm{x}}_{0}+\sqrt{1-\bar{\alpha}_{t}}\bm{\epsilon},$
(2)

where $\epsilon\sim{\mathcal{N}}(0,{\mathbf{I}})$. This allows a closed-form
expression for the marginal distribution for sampling ${\bm{x}}_{t}$ given
${\bm{x}}_{0}$

$q({\bm{x}}_{t}|{\bm{x}}_{0}):={\mathcal{N}}({\bm{x}}_{t};\sqrt{\bar{\alpha}_{t}}{\bm{x}}_{0},(1-\bar{\alpha}_{t}){\mathbf{I}}).$
(3)

The goal of the reverse process is to generate a clean image ${\bm{x}}_{0}$
given a noise vector ${\bm{x}}_{T}$ and the condition ${\bm{y}}$. This is
achieved by training a convolutional neural network (CNN) to reverse the
Markov Chain from ${\bm{x}}_{T}$ to ${\bm{x}}_{0}$. The learning is formulated
as the estimation of a parameterized Gaussian process

$p({\bm{x}}_{t-1}|{\bm{x}}_{t}):={\mathcal{N}}({\bm{x}}_{t-1};\bm{\mu}_{\bm{\theta}}({\bm{x}}_{t},t,{\bm{y}}),\sigma^{2}_{t}{\mathbf{I}}),$
(4)

where the mean $\bm{\mu}_{\bm{\theta}}({\bm{x}}_{t},t)$ is only uncertain
variable we need to estimate and
$\sigma^{2}_{t}=\frac{1-\bar{a}_{t}}{1-\bar{a}_{t-1}}\beta_{t}$. cDDPM uses a
CNN trained to estimate $\bm{\epsilon}_{\bm{\theta}}$ to estimate added noise
$\bm{\epsilon}$ to reparametrize sampling from the conditional distribution
into

${\bm{x}}_{t-1}=\frac{1}{\sqrt{a_{t}}}({\bm{x}}_{t}-\frac{\beta_{t}}{\sqrt{1-\bar{a}_{t}}}\bm{\epsilon}_{\bm{\theta}}({\bm{x}}_{t},t,{\bm{y}}))+\sigma_{t}{\bm{z}},$
(5)

where ${\bm{z}}\sim{\mathcal{N}}(0,{\mathbf{I}})$. The CNN
$\bm{\epsilon}_{\bm{\theta}}$ used in our implementation is learned by
minimizing the following expected loss

$\mathcal{L}({\bm{\theta}})=\mathbb{E}_{\bm{\epsilon},{\bm{x}}_{t},t}[\left\lVert\bm{\epsilon}_{\bm{\theta}}({\bm{x}}_{t},t,{\bm{y}})-\bm{\epsilon}\right\rVert_{2}^{2}].$
(6)

###  2.2 DiffGEPCI for High-Quality Volumetric Synthesis

We now introduce DiffGEPCI, a new framework specifically designed for
translating 3D mGRE signals into the target modalities by using 2.5D
conditional diffusion models.

Our data preprocessing standardizes each target and condition volume by
padding them to a uniform size across all planes, ensuring that the volume
dimension is $N\times N\times N$. We use notations ${\bm{x}}^{A_{i}}$,
${\bm{x}}^{C_{j}}$, and ${\bm{x}}^{S_{k}}$ to denote the 2D slices from the
the Axial, Coronal, and Sagittal planes, indexed using $i$, $j$, and $k$.
Furthermore, we represent the 3D volumes derived from the Axial, Coronal, and
Sagittal planes as ${\bm{x}}^{A}$, ${\bm{x}}^{C}$, and ${\bm{x}}^{S}$,
respectively.

DiffGEPCI is trained with three distinct target planes using the same model
with the loss function defined in (6), where ${\bm{x}}_{0}\sim
q({\bm{x}}_{0})$ and $q$ is the probability distribution of target modalities
from three planes. DiffGEPCI uses the same cDDPM model as implemented in [10].
This approach enables the use of a single model $\bm{\epsilon}_{\bm{\theta}}$
to estimate all planes in the target modality.

As depicted in Figure 1(a), during the inference phase, we initiate the
prediction process by estimating each slice of the target modality along the
Axial plane in a 1000-step iterative process, as defined in Algorithm 1. Once
we predict all the axial slices of the 3D volume, we obtain the initial target
volume, denoted as
${\bm{x}}_{\text{inital}}^{\text{3D}}={\bm{x}}_{0}^{A}=({\bm{x}}_{0}^{A_{i}})^{N}_{i=1}$.
Although ${\bm{x}}_{\text{inital}}^{\text{3D}}$ provides reasonably good
performance, it still exhibits artifacts in the Coronal and Sagittal planes,
as shown in Fig. 3.

To address inconsistencies along Coronal and Sagital planes, we introduce a
subsequent refinement module illustrated in Fig. 1(b). First, we split the
${\bm{x}}_{\text{inital}}^{\text{3D}}$ into coronal and sagittal planes,
creating ${\bm{x}}^{C}=({\bm{x}}^{C_{j}})^{N}_{j=1}$ and
${\bm{x}}^{S}=({\bm{x}}^{S_{k}})^{N}_{k=1}$. Next, we apply 10 forward
diffusion steps to introduce a minor amount of noise to each slice, resulting
in ${\bm{x}}_{10}^{C}=({\bm{x}}_{10}^{C_{j}})^{N}_{j=1}$ and
${\bm{x}}_{10}^{S}=({\bm{x}}_{10}^{S_{k}})^{N}_{i=1}$. Following this, as
shown in Algorithm 2, the same model $\bm{\epsilon}_{\bm{\theta}}$ is used to
execute 10 steps of an inverse diffusion process for each slice, yielding
${\bm{x}}_{0}^{C}$ and ${\bm{x}}_{0}^{S}$. Finally, the three volumes are
fused to obtain the final output ${\bm{x}}_{\text{final}}^{\text{3D}}$.

##  3 Experimental Validation

###  3.1 Dataset

All human studies were approved by Washington University IRB. We collected
high-resolution 3D brain MRI data from 29 volunteers using a Siemens 3T Trio
MRI scanner and a 32-channel phased-array head coil, obtaining mGRE, FLAIR,
and MPRAGE images at ($1\times 1\times 1$mm) resolution for each subject to
create paired datasets. mGRE data were acquired using a three dimensional (3D)
multi-gradient-echo sequence with FOV 256 mm$\times$192 mm, repetition time TR
= 50 ms, flip angle 35°, 10 gradient echoes with first gradient echo time TE1
= 4 ms, echo spacing = 4 ms. Additional phase stabilization echo (the
navigator) is collected for each line in k-space to correct for image
artefacts due to the physiological fluctuations [15]. The Generalized
Autocalibrating Partially Parallel Acquisitions (GRAPPA) [16] algorithm was
used to reduce imaging time to 10 minutes. The 29 subjects were split into 23,
2, and 4 for training, validation, and testing, respectively. We extracted 85
slices from each subject’s brain volume, containing the most relevant brain
regions. We co-registered the corresponding FLAIR and MPRAGE data to the mGRE
sequence using a 3D Affine transformation from the Dipy package [17].

Algorithm 1 : Initial synthesis

1: input: $\bm{\epsilon}_{\bm{\theta}}$: pre-trained model, ${\bm{y}}$: mGRE
signal

2: for $i=N,N-1,\cdots,0$ do

3: Sample ${\bm{x}}_{1000}^{A_{i}}\sim{\mathcal{N}}(0,{\mathbf{I}})$

4: for $t=1000,999,\cdots,0$ do

5: ${\bm{z}}\sim{\mathcal{N}}(0,{\mathbf{I}})$

6:
${\bm{x}}_{t-1}^{A_{i}}=\frac{1}{\sqrt{a_{t}}}({\bm{x}}_{t}^{A_{i}}-\frac{\beta_{t}}{\sqrt{1-\bar{a}_{t}}}\bm{\epsilon}_{\bm{\theta}}({\bm{x}}_{t}^{A_{i}},t,{\bm{y}}^{A_{i}}))+\sigma_{t}{\bm{z}}$.

Algorithm 2 : 2.5D Refinement Module

1: input: $\bm{\epsilon}_{\bm{\theta}}$: pre-trained model; ${\bm{y}}$: mGRE
signal; ${\bm{x}}_{0}^{A}$, ${\bm{x}}_{10}^{C}$ and ${\bm{x}}_{10}^{S}$:
initialization volume,

2: for $j=N,N-1,\cdots,0$ do

3: for $t=10,9,\cdots,0$ do

4: ${\bm{z}}\sim{\mathcal{N}}(0,{\mathbf{I}})$

5:
${\bm{x}}_{t-1}^{C_{j}}=\frac{1}{\sqrt{a_{t}}}({\bm{x}}_{t}^{C_{j}}-\frac{\beta_{t}}{\sqrt{1-\bar{a}_{t}}}\bm{\epsilon}_{\bm{\theta}}({\bm{x}}_{t}^{C_{j}},t,{\bm{y}}^{C_{j}}))+\sigma_{t}{\bm{z}}$.

6: for $k=N,N-1,\cdots,0$ do

7: for $t=10,9,\cdots,0$ do

8: ${\bm{z}}\sim{\mathcal{N}}(0,{\mathbf{I}})$

9:
${\bm{x}}_{t-1}^{S_{k}}=\frac{1}{\sqrt{a_{t}}}({\bm{x}}_{t}^{S_{k}}-\frac{\beta_{t}}{\sqrt{1-\bar{a}_{t}}}\bm{\epsilon}_{\bm{\theta}}({\bm{x}}_{t}^{S_{k}},t,{\bm{y}}^{S_{k}}))+\sigma_{t}{\bm{z}}$.

10: return
${\bm{x}}_{\text{final}}^{\text{3D}}=\frac{1}{3}({\bm{x}}_{0}^{A}+{\bm{x}}_{0}^{C}+{\bm{x}}_{0}^{S}).$

Tasks |  mGRE $\rightarrow$ FLAIR
---|---
planes | Axial | Coronal | Sagitta
Metrics | PSNR | SSIM | PSNR | SSIM | PSNR | SSIM
U-Net [18] | 21.58 | 0.510 | 21.80 | 0.607 | 20.45 | 0.677
Pix2Pix [7] | 24.22 | 0.898 | 23.09 | 0.857 | 22.89 | 0.843
cDDPM [11] | 26.75 | 0.914 | 24.26 | 0.872 | 24.00 | 0.864
[5pt/5pt] DiffGEPCI | 26.92 | 0.922 | 24.88 | 0.889 | 24.38 | 0.885

Tasks |  mGRE $\rightarrow$ MPRAGE
---|---
planes | Axial | Coronal | Sagitta
Metrics | PSNR | SSIM | PSNR | SSIM | PSNR | SSIM
U-Net [18] | 25.70 | 0.597 | 25.45 | 0.696 | 26.06 | 0.784
Pix2Pix [7] | 28.31 | 0.948 | 26.48 | 0.912 | 26.45 | 0.922
cDDPM [11] | 29.07 | 0.950 | 27.57 | 0.917 | 27.17 | 0.925
[5pt/5pt] DiffGEPCI | 29.72 | 0.956 | 28.16 | 0.930 | 27.95 | 0.938

Table 1: Quantitatively comparison of DiffGEPCI and several baselines for mGRE
$\rightarrow$ FLAIR and mGRE $\rightarrow$ MPRAGE translation task on testset.
The best and  second best results are highlighted. Note the excellent
quantitative performance of DiffGEPCI, which suggests the potential of 2.5D
refinement can significantly improve the performance in all three planes.
Figure 3: Performance Comparison of various baselines and our proposed
DiffGEPCI for mGRE $\rightarrow$ MPRAGE translation in Coronal and Sagittal
planes. The top row shows the Coronal view, while the bottom row presents the
Sagittal plane.

###  3.2 Comparisons

We compared DiffGEPCI with three baseline methods. The first method, called
U-Net [18], corresponds to directly use of a U-Net model to translate the
input images into corresponding output images. The second method is pix2pix,
which is a cGAN model trained with architecture and loss functions adopted
from [7]. The third method is cDDPM [11], which is a 2D conditional DDPM
trained to generate modality outputs based on the condition provided by the
other. We used the peak signal-to-noise ratio (PSNR) and structural similarity
index (SSIM) as image quality metrics. Before we calculate these quantitative
metrics, a brain extraction tool, implemented in [19], is used to generate
REMs to mask out both skull and background in all modalities.

###  3.3 Results

Figure 2 shows the visual performance of mGRE image translation to FLAIR and
MPRAGE sequences. It is clear that U-Net introduces blurring artifacts and
detail loss. Pix2Pix outperform U-Net due to its adversarial nature, but still
results in blurry images. The traditional 2D diffusion model excels in
recovering details compared with the previous two methods. However, upon
closer examination of the generated volume in the other two planes, cross-
slice artifacts become evident, as depicted in Figure 3. By using the 2.5D
refinement module, DiffGEPCI successfully mitigates these artifacts, attaining
state-of-the-art performance across all three planes, as demonstrated in both
Figure 2 and Figure 3. The quantitative results over the whole testing set is
presented in Table 1, highlighting the quality improvements achieved by
DiffGEPCI.

##  4 Conclusion

DiffGEPCI addresses an important problem of generating translated 3D MRI
modalities from a given mGRE signal. Instead of relying solely on a 2D
diffusion model that introduces artifacts in orthogonal planes, our novel 2.5D
strategy has demonstrated the ability to significantly enhances the image
quality of the estimated volumes across all three imaging planes. Through
numerical validation using experimentally collected data, DiffGEPCI
outperformed several baseline algorithms both quantitatively and
qualitatively.

##  5 Compliance with Ethical Standards

This study was performed in line with the principles of the Declaration of
Helsinki. Approval was granted by the Ethics Committee of Washington
University (Date: 09/05/23/ No: 201906092-1001).

## Acknowledgment

This work was supported was supported by the NSF CAREER award CCF-2043134 and
the the National Institutes of Health (NIH) awards RF1 AG082030 and RF1
AG077658.

## References

  * [1] A. Panda, B. B Mehta, S. Coppo, Y. Jiang, D. Ma, N. Seiberlich, M. A Griswold, and V. Gulani,  “Magnetic resonance fingerprinting–an overview,”  Current opinion in biomedical engineering, vol. 3, pp. 56–66, 2017\.
  * [2] J. Luo, B. D Jagadeesan, A. H Cross, and D. A Yablonskiy,  “Gradient echo plural contrast imaging—signal model and derived contrasts: T2*, T1, phase, SWI, T1f, FST2* and T2*-SWI,”  Neuroimage, vol. 60, no. 2, pp. 1073–1082, 2012.
  * [3] R. Bakshi, S. Ariyaratana, R. HB Benedict, and L. Jacobs,  “Fluid-attenuated inversion recovery magnetic resonance imaging detects cortical and juxtacortical multiple sclerosis lesions,”  Archives of neurology, vol. 58, no. 5, pp. 742–748, 2001.
  * [4] John P Mugler III and J. R Brookeman,  “Rapid three-dimensional T1-weighted MR imaging with the MP-RAGE sequence,”  Journal of Magn. Reson. Med., vol. 1, no. 5, pp. 561–567, 1991\.
  * [5] A. D Elster,  “Gradient-echo MR imaging: techniques and acronyms.,”  Radiology, vol. 186, no. 1, pp. 1–8, 1993.
  * [6] Y. Pang, J. Lin, T. Qin, and Z. Chen,  “Image-to-image translation: Methods and applications,”  IEEE Trans. on Multimedia, vol. 24, pp. 3859–3881, 2021.
  * [7] P. Isola, J. Zhu, T. Zhou, and A. A Efros,  “Image-to-image translation with conditional adversarial networks,”  in Proc. IEEE Conf. Comput. Vis. Pattern Recognit., 2017, pp. 1125–1134.
  * [8] S. U. Dar, M. Yurt, L. Karacan, A. Erdem, E. Erdem, and T. Cukur,  “Image synthesis in multi-contrast MRI with conditional generative adversarial networks,”  IEEE Trans. Med. Imag., vol. 38, no. 10, pp. 2375–2388, 2019.
  * [9] K. Armanious, C. Jiang, M. Fischer, T. Küstner, T. Hepp, K. Nikolaou, S. Gatidis, and B. Yang,  “MedGAN: Medical image translation using GANs,”  Comput. med. imag. and graph., vol. 79, pp. 101684, 2020.
  * [10] P. Dhariwal and A. Nichol,  “Diffusion models beat gans on image synthesis,”  Proc. Ann. Conf. Neural Information Processing Systems, vol. 34, pp. 8780–8794, 2021.
  * [11] J. Ho, A. Jain, and P. Abbeel,  “Denoising diffusion probabilistic models,”  Proc. Ann. Conf. Neural Information Processing Systems, vol. 33, pp. 6840–6851, 2020.
  * [12] C. Saharia, W. Chan, H. Chang, C. Lee, J. Ho, T. Salimans, D. Fleet, and M. Norouzi,  “Palette: Image-to-image diffusion models,”  in ACM SIGGRAPH 2022 Conf. Proceedings, 2022, pp. 1–10.
  * [13] A. Kazerouni, E. K. Aghdam, M. Heidari, R. Azad, M. Fayyaz, I. Hacihaliloglu, and D. Merhof,  “Diffusion models in medical imaging: A comprehensive survey,”  Med. Image Analysis, p. 102846, 2023.
  * [14] M. Özbey, O. Dalmaz, S. U. Dar, H. A Bedel, Ş. Özturk, A. Güngör, and T. Çukur,  “Unsupervised medical image translation with adversarial diffusion models,”  IEEE Trans. Med. Imag., 2023.
  * [15] J. Wen, A. H Cross, and D. A Yablonskiy,  “On the role of physiological fluctuations in quantitative gradient echo MRI: implications for GEPCI, QSM, and SWI,”  Magn. Reson. Med., vol. 73, no. 1, pp. 195–203, 2015.
  * [16] M. A Griswold, P. M Jakob, R. M Heidemann, M. Nittka, V. Jellus, J. Wang, B. Kiefer, and A. Haase,  “Generalized autocalibrating partially parallel acquisitions (GRAPPA),”  Magn. Reson. Med., vol. 47, no. 6, pp. 1202–1210, 2002.
  * [17] E. Garyfallidis, M. Brett, B. Amirbekian, A. Rokem, Stefan Van D. W., Maxime Descoteaux, Ian Nimmo-Smith, and Dipy Contributors,  “Dipy, a library for the analysis of diffusion MRI data,”  Frontiers in neuroinformatics, vol. 8, pp. 8, 2014.
  * [18] O. Oktay, J. Schlemper, L. L. Folgoc, M. Lee, M. Heinrich, K. Misawa, K. Mori, S. McDonagh, N. Y Hammerla, B. Kainz, et al.,  “Attention U-net: Learning where to look for the pancreas,”  arXiv preprint arXiv:1804.03999, 2018.
  * [19] M. Jenkinson, M. Pechaud, and Smith S.,  “BET2: MR-based estimation of brain, skull and scalp surfaces,”  in Eleventh Annual Meeting of the Organization for Human Brain Mapping, 2005, p. 167.


# Zero-Order Optimization for
Gaussian Process-based Model Predictive Control

A. Lahr1, A. Zanelli1, A. Carron1 and M. N. Zeilinger1 1Institute for Dynamic
Systems and Control, ETH Zürich, Zürich CH-8092, Switzerland. E-mail:
<EMAIL_ADDRESS>project has received funding by the European Union’s
Horizon 2020 research and innovation programme, Marie Skłodowska-Curie grant
agreement No. 953348, ELO-X.

###### Abstract

By enabling constraint-aware online model adaptation, model predictive control
using Gaussian process (GP) regression has exhibited impressive performance in
real-world applications and received considerable attention in the learning-
based control community. Yet, solving the resulting optimal control problem in
real-time generally remains a major challenge, due to i) the increased number
of augmented states in the optimization problem, as well as ii)
computationally expensive evaluations of the posterior mean and covariance and
their respective derivatives. To tackle these challenges, we employ i) a
tailored Jacobian approximation in a sequential quadratic programming (SQP)
approach and combine it with ii) a parallelizable GP inference and automatic
differentiation framework. Reducing the numerical complexity with respect to
the state dimension $n_{x}$ for each SQP iteration from ${O}(n_{x}^{6})$ to
${O}(n_{x}^{3})$, and accelerating GP evaluations on a graphical processing
unit, the proposed algorithm computes suboptimal, yet feasible, solutions at
drastically reduced computation times and exhibits favorable local convergence
properties. Numerical experiments verify the scaling properties and
investigate the runtime distribution across different parts of the algorithm.

##  I Introduction

Real-world applications of model predictive control using Gaussian processes
(GP-MPC), such as vision-based robot path-tracking [1], trajectory tracking
using a robotic arm [2], autonomous racing [3, 4], or high-speed quadrotor
flight [5], have showcased its potential to leverage closed-loop data for
constraint-aware online model adaptation.

Yet, the computational cost associated with GP inference remains a major
challenge. To fully exploit the rich, state-dependent uncertainty description
induced by the posterior covariance of the GP, it should be included in the
optimal control problem (OCP) formulation and propagated through the dynamics
model over the prediction horizon [6]. However, incorporating the covariance
propagation into the OCP constitutes a major limitation for real-time
implementation of the algorithm on embedded hardware, both in terms of
additional optimization variables capturing the state covariance, as well as
computationally expensive evaluation and differentiation of the GP posterior
covariance.

As a result, practical implementations resort to various heuristics to speed
up the algorithm. Popular approaches include GP approximations with a fixed
number of basis functions [7, 8], reduction of optimization variables by
fixing the state covariances in the optimization problem based on their
predicted value at the last MPC iteration [3, 6], or to completely ignore the
uncertainty description provided by the GP posterior covariance in the control
algorithm [1], [5].

###  I-A Contributions

This paper addresses the aforementioned challenges by

  * i)

applying a tailored inexact sequential quadratic programming (SQP) algorithm
[9, 10] to solve the GP-MPC optimal control problems, returning suboptimal,
yet feasible solutions at convergence,

  * ii)

speeding up GP computations by efficient inference and automatic
differentiation (AD) [11, 12], as well as parallelization on a graphical
processing unit (GPU),

  * iii)

showing that the resulting algorithm maintains favorable local convergence
properties if the GP posterior covariance and process noise covariance at the
solution are sufficiently small.

###  I-B Related Work

A tailored Jacobian approximation for stochastic and robust nonlinear model
predictive control (NMPC) has been originally developed in [9, 10] where it is
referred to as a zero-order method. Based on this approach, [13] presents a
robust NMPC algorithm that improves the ellipsoidal uncertainty propagation by
efficiently optimizing for optimal linear feedback policies used in the
predictions. In the stochastic setting, in [14] the tailored Jacobian
approximation is coupled with a more accurate uncertainty propagation based on
linear-regression Kalman filtering, showing real-time feasible timings for an
automotive example. With respect to the aforementioned works, this paper
presents a direct extension of [9, 10] to the setting of Gaussian process-
based MPC. As such, the improvements presented in [13, 14] can be
straightforwardly applied to the results in this paper as well.

To the best of the authors’ knowledge, [8] is the only publication that
applies the tailored Jacobian approximation to a variant of GP-MPC. Therein,
the tailored stochastic NMPC algorithm is used in conjunction with an
approximate GP based on a finite set of basis functions for an autonomous
driving application, where the distributions of the weights are estimated
jointly with the system state using a particle filter. Conceptually, the
presented work differs from [8] insofar as we apply exact GP inference with an
approximate uncertainty propagation, while in [8], a finite-dimensional GP
approximation is used and the uncertainty propagation is performed using a
particle filter. Moreover, we additionally focus on a theoretical analysis of
the convergence properties and provide a high-performance implementation of
the proposed algorithm.

##  II Problem formulation

We consider discrete-time, nonlinear dynamics of the form

$\displaystyle x(k+1)=\psi(x(k),u(k))+B\left(\eta(x(k),u(k))+w(k)\right),$ (1)

where
$\psi:\mathbb{R}^{n_{x}}\times\mathbb{R}^{n_{u}}\rightarrow\mathbb{R}^{n_{x}}$,
$\eta:\mathbb{R}^{n_{x}}\times\mathbb{R}^{n_{u}}\rightarrow\mathbb{R}^{n_{w}}$
describe the known and unknown parts of the system dynamics as functions of
the state $x(k)\in\mathbb{R}^{n_{x}}$ and input $u(k)\in\mathbb{R}^{n_{u}}$ at
time step $k$, respectively, and are assumed to be twice continuously
differentiable. The process noise $w(k)\sim{N}(0,\Sigma^{w})$ is assumed to be
i.i.d. in each component, i.e., with diagonal covariance matrix
$\Sigma^{w}=\mathrm{diag}\left(\sigma_{1}^{2},\ldots,\sigma_{n_{w}}^{2}\right)$;
it affects the states via the matrix $B\in\mathbb{R}^{n_{x}\times n_{w}}$,
allowing to model the process noise in a lower-dimensional space. The system
is subject to $n_{h}$ individual chance constraints, i.e.,

$\displaystyle\mathrm{Pr}(h_{j}(x(k),u(k))\leq 0)\geq p_{j},$ (2)

for all constraints $j=1,\ldots,n_{h}$ and time steps $k$, with satisfaction
probability $0<p_{j}\leq 1$. This formalism also captures hard constraints,
for example on the inputs, by setting $p_{j}=1$ for the respective constraint.

The key idea of GP-MPC is to model the unknown dynamics $\eta$ as a Gaussian
process

$\displaystyle d(x,u)\sim{GP}\left(\mu^{d}(x,u),\Sigma^{d}(x,u)\right),$ (3)

with posterior mean
$\mu^{d}:\left(\mathbb{R}^{n_{x}}\times\mathbb{R}^{n_{u}}\right)\rightarrow\mathbb{R}^{n_{w}}$
and posterior covariance
$\Sigma^{d}:(\mathbb{R}^{n_{x}}\times\mathbb{R}^{n_{u}})\rightarrow\mathbb{R}^{n_{w}\times
n_{w}}$, respectively. Note that by $d(x,u)$ we denote a GP conditioned on
data points already; see e.g. [15, Chap. 2] for the inference formulas. At
every sampling time, the problem to be solved can be formulated as a
stochastic OCP111In (4), terminal constraints have been omitted for
simplicity., where a sequence of control polices
$\bm{\pi}:=\\{\pi_{i}\\}_{i=0}^{N-1}$,
$\pi:\mathbb{R}^{n_{x}}\rightarrow\mathbb{R}^{n_{u}}$, is determined over a
given set of admissible policies, that minimizes the expected value of a user-
defined cost over a finite time horizon of length $N$,

$\displaystyle\underset{\bm{\pi}}{\min}\quad$
$\displaystyle\mathbb{E}\Bigg{[}c_{f}(x_{N})+\sum_{i=0}^{N-1}c_{i}(x_{i},u_{i})\Bigg{]}$
(4a) $\displaystyle\mathrm{s.t.}\quad$ $\displaystyle\forall
i\in\\{0,\ldots,N-1\\},$ (4b) $\displaystyle u_{i}=\pi_{i}(x_{i}),$ (4c)
$\displaystyle x_{i+1}=\psi(x_{i},u_{i})+B\left(d(x_{i},u_{i})+w_{i}\right),$
(4d) $\displaystyle\mathrm{Pr}(h_{j}(x_{i},u_{i})\leq 0)\geq
p_{j},\>j=1,\ldots,n_{h},$ (4e) $\displaystyle x_{0}=x(k).$ (4f)

As computing exact solutions to the stochastic OCP (4) is generally
intractable, in the following we present a common approach to approximate (4)
by a deterministic problem, as previously proposed by [6].

###  II-A Mean and covariance propagation

Propagating the uncertainty introduced by the GP model through the nonlinear
dynamics model in a computationally efficient way generally proves to be very
challenging: As Gaussianity of the state distribution is lost after a
nonlinear transformation, so is the possibility to capture the entire
distribution by its first two moments. To retain computational tractability,
we use a common linearization-based approximation of the propagation of the
expected state and covariance [16], which results in the following
deterministic update equations [17], starting from
$x_{0}\sim{N}(\mu^{x}_{0},\Sigma^{x}_{0})$,

$\displaystyle\mu_{i+1}^{x}$
$\displaystyle=\psi(\mu_{i}^{x},u_{i})+B\mu^{d}(\mu_{i}^{x},u_{i}),$ (5)
$\displaystyle\Sigma_{i+1}^{x}$
$\displaystyle=\tilde{A}_{i}\Sigma_{i}^{x}\tilde{A}_{i}^{\top}+B\left(\Sigma^{d}(\mu_{i}^{x},u_{i})+\Sigma^{w}\right)B^{\top},$
(6) where $\displaystyle\tilde{A}_{i}$
$\displaystyle:=\left.\frac{\partial}{\partial
x}\left(\psi(x,u_{i})+B\mu^{d}(x,u_{i})\right)\right|_{x=\mu^{x}_{i}}.$ (7)

###  II-B Chance constraints

Given the covariance propagation above, it is possible to efficiently
formulate the individual chance constraints (4e) as deterministic constraints
on the mean prediction by a suitable tightening, i.e.,
$\bar{h}_{j}(\mu^{x}_{i},u_{i},\Sigma^{x}_{i})\leq 0$, where

$\displaystyle\bar{h}_{j}(\mu^{x}_{i},u_{i},\Sigma^{x}_{i}):=h_{j}(\mu^{x}_{i},u_{i})+\alpha_{j}\sqrt{C_{j}(\mu^{x}_{i},u_{i})\Sigma^{x}_{i}C_{j}(\mu^{x}_{i},u_{i})^{\top}}$

and $C_{j}(\mu^{x}_{i},u_{i}):=\frac{\partial h_{j}}{\partial
x}(\mu^{x}_{i},u_{i})$. For general probability distributions, the tightening
factor $\alpha:=\sqrt{\frac{p_{j}}{1-p_{j}}}$ can be chosen based on the
Chebyshev inequality; for Gaussian distributions, setting
$\alpha_{j}:=\Phi^{-1}(p_{j})$, where $\Phi^{-1}(\cdot)$ is the inverse
cumulative density function of a standard normal Gaussian variable, is a less
conservative choice.

###  II-C Optimal control problem

With the above simplifications and an approximate formulation of the expected
cost, we arrive at the following deterministic approximation of the stochastic
OCP (4) in terms of the predicted state mean
$\bm{\mu}:=\\{\mu^{x}_{i}\\}_{i=0}^{N}$, state covariance
$\bm{\Sigma}:=\\{\Sigma^{x}_{i}\\}_{i=0}^{N}$ and control input sequence
$\bm{u}:=\\{u_{i}\\}_{i=0}^{N-1}$, cf. [6],

$\displaystyle\underset{\bm{\mu},\bm{\Sigma},\bm{u}}{\min}\quad$
$\displaystyle c_{f}(\mu_{N}^{x})+\sum_{i=0}^{N-1}c_{i}(\mu_{i}^{x},u_{i})$
(8a) $\displaystyle\mathrm{s.t.}\quad$ $\displaystyle\forall
i\in\\{0,\ldots,N-1\\},$ (8b)
$\displaystyle\mu_{i+1}^{x}=\psi(\mu_{i}^{x},u_{i})+B\mu^{d}(\mu_{i}^{x},u_{i}),$
(8c)
$\displaystyle\Sigma_{i+1}^{x}=\tilde{A}_{i}\Sigma_{i}^{x}\tilde{A}_{i}^{\top}+B\left(\Sigma^{d}(\mu_{i}^{x},u_{i})+\Sigma^{w}\right)B^{\top},$
(8d) $\displaystyle\bar{h}_{j}(\mu^{x}_{i},u_{i},\Sigma^{x}_{i})\leq
0,\>j=1,\ldots,n_{h},$ (8e) $\displaystyle\mu_{0}^{x}=x(k),$ (8f)
$\displaystyle\Sigma_{0}^{x}=0.$ (8g)

###### Remark 1

While optimizing over control input sequences in (8) has computational
advantages compared to optimization over a sequence of feedback policies in
(4), it can lead to underconfident predictions and an overly conservative
controller. As a remedy, linear state feedback can be incorporated into the
controller [6]. In that case, hard input constraints can generally not be
satisfied and should be replaced with individual chance constraints; the
adaptation of the OCPs (4) and (8) is straightforward and can be found in [6].

Solving (8) in a receding horizon fashion yields a highly performant and
adaptive, yet uncertainty-aware, control strategy, demonstrated by real-world
applications such as [3, 5]. Nevertheless, for a real-time GP-MPC
implementation, there remain major computational challenges to be addressed.

First, the OCP (8) has $n_{x}+n_{u}+(n_{x}+n_{x}^{2})/2$ optimization
variables at each prediction stage, structurally equivalent to OCPs arising
from a direct multiple shooting discretization [18]. For this set of problems,
among the most competitive solvers at the time of this writing are interior
point algorithms exploiting the block-sparse structure of the multiple
shooting approach [19]. Still, these algorithms’ cubic computational
complexity in the number of shooting nodes leads to a total computational
complexity of ${O}(n_{x}^{6})$, which becomes prohibitive even for a moderate
number of states.

Second, even without updating the GP data points online, evaluating the GP
posterior mean, covariance, and their respective Jacobians generally scales
quadratically with the number of data points [15, p. 19], posing strict
limitations on the number of data points amenable for online inference. In
this regard, solving (8) is particularly demanding, as computing the
constraint Jacobians (8d) requires not only to compute the Jacobian of the
nominal and GP dynamics but also, their respective second-order derivatives.

A popular heuristic to ensure computational tractability despite the
limitations stated above has been to propagate the state covariance (8d)
outside the optimizer, based on input and state predictions from the last MPC
iteration [3], [6]. Fixing the covariance matrices $\Sigma^{x}_{i}$ at each
stage in (8) eliminates the corresponding $(n_{x}+n_{x}^{2})/2$ augmented
states as well as the need to compute the equality constraint Jacobian derived
from (8d). However, this heuristic generally does not lead to feasible
solutions for problem (8), as Fig. 1 illustrates.

$x_{2}$
---
$x_{1}$
---
Figure 1: Infeasibility arising from fixing the covariances based on the
previous MPC instance. The predicted state trajectory and covariances from the
previous time step are drawn with solid lines in blue; predicted covariances
around current predicted state trajectory in red; reference for
previous/current time step with dashed lines; infeasible region in light red.
When the linearized dynamics at the shooting nodes vary strongly from one time
step to another, in this example based on a reference change, fixing the
covariances based on the previous MPC instance might lead to significant
prediction errors.

To address the aforementioned challenges, we propose a tailored sequential
quadratic programming (SQP) algorithm that computes suboptimal, yet feasible,
solutions to (8) with drastically reduced computational footprint.

##  III Zero-order method for GP-MPC

In the following, we present a tailored optimization method for Gaussian
process-based MPC based on inexact sequential quadratic programming. To avoid
excessive notation in the following exposition, let us rewrite the GP-MPC
optimal control problem (8) in the following compact form,

$\displaystyle\underset{y,P}{\min}\quad$ $\displaystyle c(y)$ (9a)
$\displaystyle\mathrm{s.t.}\quad$ $\displaystyle f(y)=0,$ (9b) $\displaystyle
g(y,P)=0,$ (9c) $\displaystyle\bar{h}(y,P)\leq 0.$ (9d)

Thereby we have concatenated all mean states and inputs into
$y:=(y_{0},\ldots,y_{N-1},\mu^{x}_{N})$ and vectorized covariances into
$P:=(\overrightarrow{\Sigma^{x}_{0}},\ldots,\overrightarrow{\Sigma^{x}_{N}})$,
where $y_{i}:=(\mu_{i}^{x},u_{i})$ for all $i=0,\ldots,N-1$ and
$\overrightarrow{(\cdot)}$ is the column-wise vectorization operator.

###  III-A Sequential quadratic programming

Sequential quadratic programming [20] finds a solution to (9) by iteratively
solving quadratic programs that locally approximate the original NLP at the
current solution estimate $(\hat{y},\hat{P})$. At each iteration, the solution
$(\Delta y,\Delta P)$ of the quadratic subproblem

$\displaystyle\underset{\Delta y,\Delta P}{\min}\quad$
$\displaystyle\frac{1}{2}\Delta y^{\top}M_{yy}\Delta y+\left.\frac{\partial
c}{\partial y}\right|_{y=\hat{y}}\Delta y$ (10a)
$\displaystyle\mathrm{s.t.}\>\quad$ $\displaystyle\begin{bmatrix}f\\\
g\end{bmatrix}(\hat{y},\hat{P})+\left.\begin{bmatrix}\frac{\partial
f}{\partial y}&0\\\ \frac{\partial g}{\partial y}&\frac{\partial g}{\partial
P}\end{bmatrix}\right|_{\overset{\scriptstyle
y=\hat{y}}{P=\hat{P}}}\begin{bmatrix}\Delta y\\\ \Delta P\end{bmatrix}=0,$
(10b)
$\displaystyle\bar{h}(\hat{y},\hat{P})+\left.\begin{bmatrix}\frac{\partial\bar{h}}{\partial
y}&\frac{\partial\bar{h}}{\partial
P}\end{bmatrix}\right|_{\overset{\scriptstyle
y=\hat{y}}{P=\hat{P}}}\begin{bmatrix}\Delta y\\\ \Delta P\end{bmatrix}\leq 0,$
(10c)

is then used to update the current solution estimate, where $M_{yy}$ denotes
the chosen approximation of the Hessian of the Lagrangian of (9). The core SQP
algorithm is summarized in Alg. 1.

1 Initialize $\hat{y},\hat{P}$;

2 while  _termination criterion not met_ do

3 Solve QP (10) for $\Delta y,\Delta P$;

4 Set $y_{+}:=\hat{y}+\Delta y$, $P_{+}:=\hat{P}+\Delta P$;

5 Update $\hat{y}:=y_{+}$, $\hat{P}:=P_{+}$;

6

7 end while

Algorithm 1 Standard SQP iteration

Under standard regularity assumptions, it can be shown that the SQP iteration
locally converges to a Karush-Kuhn-Tucker (KKT) point of the original NLP (9)
if the Hessian approximation $M_{yy}$ is sufficiently accurate, cf. [20, Thm.
3.2].

Computationally, the SQP algorithm is ideally suited for NMPC, as the
computational load can be reduced significantly by running only a fixed number
of iterations before applying the control input to the plant, usually with
negligible effect on the control performance, c.f. the Real-Time Iteration
[21]. Still, as the number of optimization variables at each shooting node
carries over from (8) to the quadratic subproblems in (10) when employing a
sparsity-exploiting interior point solver, their cubic scaling with respect to
the number of augmented states (or ${O}(n_{x}^{6})$ with respect to the state
dimension) poses a major limitation on the nominal system dimension amenable
for a GP-MPC application.

###  III-B Zero-order algorithm

To alleviate the limited scalability of optimizing over the state covariance
matrices in GP-MPC, we employ an inexact SQP method that has been initially
presented for stochastic [9] and robust NMPC [10]. The core idea is that, by
using a tailored Jacobian approximation in the SQP algorithm, the equality
constraints in (10) corresponding to the covariance propagation (9c) can be
decoupled from the optimization problem. This allows for subsequent
elimination of the associated optimization variables $\Delta P$ from the QP
(10), which we present in the following.

The tailored Jacobian approximation is obtained by neglecting the Jacobian of
the covariance propagation (9c) with respect to the nominal variables $y$,
i.e., by setting $\frac{\partial g}{\partial y}:=0$ in (10b). This leads to a
zero-order approximation of the corresponding equality constraints at each SQP
iteration. Under the approximate Jacobian, rearranging (10b) yields an
approximation of $\Delta P$ based solely on the current linearization point
$(\hat{y},\hat{P})$, i.e.,

$\displaystyle\Delta\tilde{P}=-\left(\left.\frac{\partial g}{\partial
P}\right|_{\overset{\scriptstyle
y=\hat{y}}{P=\hat{P}}}\right)^{-1}g(\hat{y},\hat{P}),$ (11)

which can be used to eliminate $\Delta P\approx\Delta\tilde{P}$ from the QP
subproblems (10).

Instead of solving the linear system in (11), however, let us show how
$\Delta\tilde{P}$ may be obtained more efficiently and intuitively by making
use of the fact that $g(y,P)$ corresponds to a vectorized version of the
covariance propagation (6). We start by noting that the corresponding,
vectorized equality constraints (9c) may be written as

$\displaystyle g(y,P)=A(y)P+b(y)=0,$ (12)

where $A(y)\in\mathbb{R}^{n_{x}^{2}(N+1)\times n_{x}^{2}(N+1)}$ is invertible
and contains the system dynamics and $b(y)\in\mathbb{R}^{n_{x}^{2}(N+1)}$, the
vectorized process noise and GP posterior covariances for all prediction
steps; see Appendix -A for an explicit construction. Due to linearity of
$g(y,P)$ in $P$, obtaining $\Delta\tilde{P}$ in (11) based on the current
linearization point $(\hat{y},\hat{P})$ is equivalent to solving for

$\displaystyle\tilde{P}_{+}=-A(\hat{y})^{-1}b(\hat{y})$ (13)

and setting $\Delta\tilde{P}:=\tilde{P}_{+}-\hat{P}$. Since (12) is obtained
by vectorization and stacking of (6) for all stages, it can easily be verified
that solving for $\tilde{P}_{+}$ corresponds to propagating the predicted
covariance as given by equation (6), based on the current iterate $\hat{y}$
and initial covariance $\Sigma^{x}_{0}:=0$, followed by an application of the
vectorization operator.

After eliminating $\Delta P\approx\Delta\tilde{P}$, problem (10) reduces to

$\displaystyle\underset{\Delta y}{\min}\quad$ $\displaystyle\frac{1}{2}\Delta
y^{\top}M_{yy}\Delta_{y}+\left.\frac{\partial c}{\partial
y}\right|_{y=\hat{y}}\Delta y$ (14a) $\displaystyle\mathrm{s.t.}\quad$
$\displaystyle f(\hat{y})+\left.\frac{\partial f}{\partial
y}\right|_{\overset{\scriptstyle y=\hat{y}}{P=\hat{P}}}\Delta y=0,$ (14b)
$\displaystyle\bar{h}(\hat{y},\hat{P})+\left.\begin{bmatrix}\frac{\partial\bar{h}}{\partial
y}&\frac{\partial\bar{h}}{\partial
P}\end{bmatrix}\right|_{\overset{\scriptstyle
y=\hat{y}}{P=\hat{P}}}\begin{bmatrix}\Delta y\\\
\Delta\tilde{P}\end{bmatrix}\leq 0.$ (14c)

The modified iteration of the SQP algorithm alternates the solution of QP (14)
with the covariance propagation based on (6) until convergence, as summarized
in Alg. 2. Due to the reduced size of the QPs, we recover again the
computational complexity of ${O}(n_{x}^{3})$ when employing sparsity-
exploiting interior point QP solvers. This provides a drastic improvement
compared with the complexity of ${O}(n_{x}^{6})$ for the QPs in the original
problem (10), as also demonstrated by the numerical examples in Section V.

1 Initialize $\hat{y},\hat{P}$;

2 while  _termination criterion not met_ do

3 Obtain $\tilde{P}_{+}=A(\hat{y})^{-1}b(\hat{y})$ from (6);

4 Set $\Delta\tilde{P}:=\tilde{P}_{+}-\hat{P}$;

5 Solve QP (14) for $\Delta y$;

6 Set $y_{+}:=\hat{y}+\Delta y$;

7 Update $\hat{y}:=y_{+}$, $\hat{P}:=\tilde{P}_{+}$;

8

9 end while

Algorithm 2 Modified SQP iteration

##  IV Local convergence properties

At convergence, the solution obtained from Alg. 2 is guaranteed to be feasible
and suboptimal for the original NLP (8) [10], [21]. This is in contrast to
existing heuristics discussed in Fig. 1, where feasibility of the obtained
solution with respect to (8) cannot be guaranteed.

Whether Alg. 2 converges to a feasible point of (8), however, depends on the
error in the Jacobian approximation. In [10], this has been studied for
equality constraints of the form $0=A(y)P-\tilde{\sigma}^{2}b(y)$, with a
scalar uncertainty parameter $\tilde{\sigma}>0$. In particular, it has been
shown that, for $\tilde{\sigma}\rightarrow 0$, the tailored Jacobian
approximation does not deteriorate the local convergence properties of the SQP
algorithm. In the following, we translate these results to the case where
$b(y)$ is encoding uncertainty based on the GP posterior covariance by
defining a parameter $\sigma:=\sup_{y\in{B}(\bar{y},r_{\gamma})}\|b(y)\|$,
which implicitly bounds the process noise and GP posterior covariance in a
small neighborhood around a fixed point $(\bar{y},\bar{P})$ of Alg. 2. This
way, for GPs based on a twice continuously differentiable kernel, we can show
that, for small enough process noise and GP posterior covariance, applying the
tailored Jacobian approximation to GP-MPC preserves the local convergence
properties of the SQP algorithm.

To this end, we first provide a brief review of the standard arguments for
local convergence analysis of inexact Newton-type methods using strongly
regular generalized equations.

###  IV-A Generalized equations

Generalized equations [22] allow us to reformulate the KKT conditions of (9)
compactly as the set inclusion

$\displaystyle 0\in F(z)+{N}_{K}(z),$ (15)

where

$\displaystyle F(z):=-\begin{bmatrix}\nabla_{(y,P)}{L}(z)\\\ f(y)\\\
A(y)P+b(y)\\\ \bar{h}(y,P)\end{bmatrix}$ (16)

and $z:=(y,P,\lambda^{\mu},\lambda^{\Sigma},\nu)$ contains the primal
variables $y,P$ and the Lagrange multipliers
$\lambda^{\mu},\lambda^{\Sigma},\nu$ corresponding to the constraints
(9b)-(9d), respectively. Thereby, ${N}_{K}(z)$ denotes the normal cone to the
set
$K:=\mathbb{R}^{n_{y}+n_{P}}\times\mathbb{R}^{n_{f}}\times\mathbb{R}^{n_{g}}\times\mathbb{R}^{n_{h}}_{+}$
at $z$.

While solutions to (15) are KKT points of the original problem (8), the
suboptimal solution obtained by Alg. 2 will instead satisfy the perturbed
generalized equation

$\displaystyle 0\in\tilde{F}(z)+{N}_{K}(z),$ (17)

where

$\displaystyle\tilde{F}(z):=-\begin{bmatrix}\nabla_{(y,P)}\tilde{{L}}(z)\\\
f(y)\\\ A(y)P+b(y)\\\ \bar{h}(y,P)\end{bmatrix}.$ (18)

Evidently, the suboptimality is thereby caused by the inexact Jacobian
approximation used in the QP subproblems, leading to a perturbed Lagrangian
gradient

$\displaystyle\nabla_{(y,P)}\tilde{{L}}(z)=\nabla_{(y,P)}{L}(z)-\begin{bmatrix}\frac{\partial}{\partial
y}\left(A(y)P+b(y)\right)^{\top}\\\ 0\end{bmatrix}\lambda^{\Sigma}$ (19)

in the stationarity conditions in (16) and (18). Every iterate of Alg. 2
solves the linearized generalized equation

$\displaystyle 0\in\tilde{F}(\hat{z})+J(\hat{z})(z-\hat{z})+{N}_{K}(z),$ (20)

where $J(\hat{z})\approx\left.\frac{\partial\tilde{F}}{\partial
z}\right|_{z=\hat{z}}$ is the Jacobian approximation around the linearization
point $\hat{z}$.

To obtain a sufficient local convergence criterion for Alg. 2 to a solution
$\bar{z}:=(\bar{y},\bar{P},\bar{\lambda}^{\mu},\bar{\lambda}^{\Sigma},\bar{\nu})$
of (17), we need the following key requirements to be satisfied. Thereby,
${B}(z,r)$ denotes a ball of radius $r\in\mathbb{R}$ centered at
$z\in\mathbb{R}^{n_{z}}$.

######  Assumption 1 (cf. [23], Ass. 1)

Let (17) be strongly regular222Roughly speaking, the strong regularity
assumption implies that there exists a single-valued and Lipschitz-continuous
localization of the solution map (as defined in [24, p.4]) of (17) with
respect to small perturbations around a point $\bar{z}$, which can also be
stated in terms of a nonsingularity condition on the Jacobian of
$\tilde{F}(z)$, c.f. the implicit-function theorem in the case of a fixed
active set [24, Thm. 1B.1]. As such, by assuming strong regularity, we
implicitly also assume differentiability of $\tilde{F}(z)$ at $\bar{z}$. For
the NLP (9) in particular, this requires differentiability of the square-root
terms present in the tightened constraints (9d) at the solution, achievable
by, e.g., a suitable design of $h(y_{i})$, $B$ and $\Sigma^{w}$. at $\bar{z}$,
with Lipschitz constant $\gamma$ in the neighborhood ${B}(\bar{z},r_{\gamma})$
defined by $r_{\gamma}>0$.

######  Assumption 2 (cf. [23], Ass. 3)

Let there exist a neighborhood ${B}(\bar{z},r_{\tilde{\kappa}})$, with
$0<r_{\tilde{\kappa}}<r_{\gamma}$, and a positive constant $\tilde{\kappa}$,
with $\gamma\tilde{\kappa}<\frac{1}{2}$, such that, for any
$\hat{z}\in{B}(\bar{z},r_{\tilde{\kappa}})$, it holds that

$\displaystyle\left\|J(\hat{z})-\left.\frac{\partial\tilde{F}}{\partial
z}\right|_{z=\bar{z}}\right\|\leq\tilde{\kappa}.$ (21)

Under Assumptions 1 and 2, a sufficient local convergence criterion for Alg. 2
reads as follows.

######  Lemma 1 (cf. [23], Lemma 2)

Let Assumptions 1 and 2 hold. Denote by $z_{+}$ a solution to (20) constructed
at the linearization point $\hat{z}$. Then, there exist strictly positive
constants $\kappa<1$ and $r_{\kappa}$, such that, for any
$\hat{z}\in{B}(\bar{z},r_{\kappa})$, it holds that

$\displaystyle\|z_{+}-\bar{z}\|\leq\kappa\|\hat{z}-\bar{z}\|.$ (22)

Note that for Newton-type optimization, Assumption 1 is standard and could be
replaced by the stronger, but more frequently encountered, assumptions of
linear independence constraint qualification and strong second-order
sufficient condition, see [22, Thm. 4.1]. Regarding Assumption 2, the goal for
the following section is to show that it is satisfied for sufficiently small
process noise and GP posterior covariance matrices, ensuring local convergence
of Alg. 2 by Lemma 1.

###  IV-B Local convergence for small uncertainties

The difference between the approximate and exact Jacobian of (18), at the
linearization point $\hat{z}$ and a solution $\bar{z}$ of (17), respectively,
is given by

$\displaystyle J(\hat{z})-\left.\frac{\partial\tilde{F}}{\partial
z}\right|_{z=\bar{z}}=\begin{bmatrix}M_{yy}-\left.\frac{\partial^{2}\tilde{{L}}}{\partial
y^{2}}\right|_{z=\bar{z}}&0\\\ 0&0\\\ \left.\frac{\partial}{\partial
y}\left(A(y)P+b(y)\right)\right|_{\overset{\scriptstyle
y=\bar{y}}{P=\bar{P}}}&0\\\ 0&0\end{bmatrix}.$ (23)

Hence, from Lemma 1 and equation (23), we can deduce local convergence in a
neighborhood around $\bar{z}$ if the errors in the Hessian- and tailored
Jacobian approximations are sufficiently small. As only the latter is
essential to our method, we focus on the Jacobian difference induced by the
tailored Jacobian approximation; see e.g. [20] for a discussion of the Hessian
approximation’s role in SQP methods.

The main result makes use of the following regularity assumption.

###### Assumption 3

Let $b(y)$ be twice continuously differentiable for all
$y\in{B}(\bar{y},r_{\gamma})$.

Note that this assumption can also be phrased in terms of the GP kernel
function’s regularity: As $b(y)$ denotes the vectorized process noise and GP
posterior covariances, the latter of which are a linear combination of kernel
function evaluations [15], we restrict ourselves to twice continuously
differentiable kernels. This captures many of the kernels commonly used in
practice, such as squared exponential, linear or Matèrn kernels with $\nu\geq
5/2$.

###### Lemma 2

Let Assumptions 1 and 3 hold and define
$\sigma:=\sup_{y\in{B}(\bar{y},r_{\gamma})}{\|b(y)\|}$. Then, for any
$\epsilon>0$, there exists a $\delta\in\mathbb{R}$ such that if
$\sigma\leq\delta$, it holds that

$\displaystyle\left\|\left.\frac{\partial}{\partial
y}\left(A(y)P+b(y)\right)\right|_{\overset{\scriptstyle
y=\bar{y}}{P=\bar{P}}}\right\|\leq\epsilon.$ (24)

###### Proof:

We will show the above implication by considering both summands inside (24)
separately and applying the triangle inequality.

For the first term, rearranging and taking the norm of the equality constraint
(12), evaluated at the solution $\bar{z}$, leads to

$\displaystyle\|\bar{P}\|=\|A(\bar{y})^{-1}b(\bar{y})\|\leq\|A(\bar{y})^{-1}\|\|b(\bar{y})\|\leq\|A(\bar{y})^{-1}\|\sigma.$

Hence, by linearity of $\frac{\partial}{\partial y}A(y)P$ in $P$ and the
triangle inequality, for any $\epsilon>0$, there exists a
$\delta\in\mathbb{R}$ such that $\sigma\leq\delta$ implies
$\left\|\left.\frac{\partial}{\partial
y}A(y)\bar{P}\right|_{y=\bar{y}}\right\|\leq\epsilon$.

For the second term, consider the Taylor expansion of the $i$-th component of
$b(y)$ around the expansion point $\bar{y}$, evaluated at
$y_{t}\in{B}(\bar{y},r_{\gamma})$,

$\displaystyle b_{i}(y_{t})=b_{i}(\bar{y})+t\left.\frac{\partial
b_{i}(y)}{\partial y}\right|_{y=\bar{y}}\delta y+\frac{t^{2}}{2}\delta
y^{\top}\left.\frac{\partial^{2}b_{i}(y)}{\partial
y^{2}}\right|_{y=\xi_{i}}\delta y,$ (25)

where $t:=\|y_{t}-\bar{y}\|$ and $\delta y:=(y_{t}-\bar{y})/t$ (for
$y_{t}\neq\bar{y}$) is a unit vector. By the mean-value theorem, the Hessian
is evaluated at some $\xi_{i}\in\Xi$ for each component $i$, with the set
$\Xi:=\left\\{\bar{y}+s(y_{t}-\bar{y})\>\middle|\>s\in[0,1]\text{ and
}\bar{y}\in\mathbb{R}^{n_{y}},y_{t}\in{B}(\bar{y},r_{\gamma})\right\\}$. Since
$b_{i}(y)$ is twice continuously differentiable, there exists some
$M_{2}\in\mathbb{R}$ such that, for all $\xi_{i}\in\Xi$,
$i\in\\{1,\ldots,n_{x}^{2}(N+1)\\}$, it holds that

$\displaystyle\left\|\left.\frac{\partial^{2}b_{i}(y)}{\partial
y^{2}}\right|_{y=\xi_{i}}\right\|\leq M_{2}.$ (26)

From (25), by taking the absolute value, dividing by $t>0$ and inserting the
definitions of $\sigma$ and $M_{2}$, we obtain

$\displaystyle\left|\left.\frac{\partial b_{i}(y)}{\partial
y}\right|_{y=\bar{y}}\delta y\right|$
$\displaystyle=\left|\frac{b_{i}(y_{t})-b_{i}(\bar{y})}{t}-\frac{t}{2}\delta
y^{\top}\left.\frac{\partial^{2}b_{i}(y)}{\partial
y^{2}}\right|_{y=\xi_{i}}\delta y\right|$
$\displaystyle\leq\frac{2\sigma}{t}+\frac{M_{2}}{2}t.$ (27)

Since (25) holds for all $0<t<r_{\gamma}$ and $\delta y$ with $\|\delta
y\|=1$, taking the supremum over $\delta y$ on the left and the infimum over
$t$ on the right-hand side preserves inequality (27), i.e.,

$\displaystyle\sup_{\|\delta y\|=1}\left|\left.\frac{\partial
b_{i}(y)}{\partial y}\right|_{y=\bar{y}}\delta y\right|$
$\displaystyle\leq\inf_{t>0}\frac{2\sigma}{t}+\frac{M_{2}}{2}t.$ (28)

By applying the definition of the induced matrix norm on the left, and solving
for the value of the infimum on the right-hand side ($t<r_{\gamma}$ can always
be achieved by increasing $M_{2}$), (28) simplifies to

$\displaystyle\left\|\left.\frac{\partial b_{i}(y)}{\partial
y}\right|_{y=\bar{y}}\right\|$ $\displaystyle\leq 2\sqrt{\sigma M_{2}}.$ (29)

Thus, we have shown that the norms of both summands in (24) scale with
${O}(\sigma)$, which proves the assertion by the triangle inequality and the
definition of ${O}(\sigma)$. ∎

To summarize, in tandem with Lemma 1 and a sufficiently accurate Hessian
approximation $M_{yy}$ across all SQP iterates, Lemma 2 establishes guaranteed
convergence of Alg. 2 to a solution $\bar{z}$ of (17) for sufficiently small
process noise and GP covariances in a local neighborhood and, practical
convergence depending on the regularity of (17) in terms of the Lipschitz
constant $\gamma$.

##  V Numerical results

Alg. 2 has been prototyped in Python using the corresponding acados interface.
For the nominal dynamics, evaluation and sensitivity computation is thereby
performed using an acados integrator and just-in-time compiled CasADi
functions [25]; for the GP mean and covariance, the corresponding computations
are carried out with PyTorch [11] using the GPyTorch library [12].333An open-
source implementation of the following example is available at
https://gitlab.ethz.ch/ics/zero-order-gp-mpc, doi:10.3929/ethz-b-000611298. In
the following, we compare different variants of the proposed algorithm against
a “naïve” GP-MPC implementation and nominal MPC using a scalable benchmarking
example.

###  V-A Hanging chain example

To test the scaling properties of Alg. 2, we apply it on a slightly modified
variant of the hanging chain example, a popular benchmark for numerical
methods for NMPC [10], [19]. The system is defined by a chain of masses
$m:=0.033$\mathrm{k}\mathrm{g}$$ connected by linear springs with stiffness
$k:=30.3$\mathrm{N}\mathrm{/}\mathrm{m}$$. The mass at one end of the chain is
tied to the origin; the $n_{u}=3$ velocity components of the other end’s mass
constitute the control inputs of the system. The system state is given by the
position components of the controlled mass as well as the position and
velocity components of the intermediate masses, resulting in a state space
dimension of $n_{x}=6(n_{\text{mass}}-2)+3$. The initial state of the chain is
computed based on its resting position, where the controlled end is placed at
$(x_{\text{init}},y_{\text{init}},z_{\text{init}}):=(6l(n_{\text{mass}}-1),0,0)$,
with length $l:=0.033$\mathrm{m}$$. After applying a control input of
$u_{\text{init}}:=(1,1,1)$ for $T_{\text{init}}:=1$\mathrm{s}$$, the control
task is to restore the resting position while ensuring that none of the masses
violates a wall constraint at $y_{\text{wall}}:=-0.05$\mathrm{m}$$. The
example is modified by adding a latent force

$\displaystyle
f_{\text{lat}}(x,v_{x}):=\alpha_{\text{lat}}\left(v_{x}-\sin\left(\beta_{1}\frac{2\pi
x}{l}\right)-\sin\left(\beta_{2}\frac{2\pi x}{l}\right)^{2}\right)^{2}$

to the $y$-acceleration of each intermediate mass, where
$\alpha_{\text{lat}}:=-0.1$, $\beta_{1}:=2$, $\beta_{1}:=3$, and $x$ and
$v_{x}$ denote the position and velocity along the $x$-axis of the frame,
respectively. The continuous-time dynamics are discretized using an implicit
Runge-Kutta integrator, with time step $T_{s}:=0.2$\mathrm{s}$$. A Gauss-
Newton Hessian approximation $M_{yy}$ is employed. The model mismatch on each
of the $n_{w}:=3(n_{\text{mass}}-2)$ velocity states of the intermediate
masses is captured using independent GPs with squared exponential kernel.
Training data is generated by recording $D:=15N_{x_{0}}$ samples of the model
mismatch from closed-loop simulations with a nominal model predictive
controller that only considers the nominal dynamics and no constraint
tightening, starting from $N_{x_{0}}$ perturbed initial conditions. The
following experiments were performed using an Intel i9-7940X processor running
at 3.10 GHz and an NVIDIA GeForce RTX 2080 Ti GPU.

###  V-B Timings for increasing state and GP output dimension

Fig. 2 shows the scaling of the mean computation times per SQP iteration as
the number of states $n_{x}$ is increased by adding masses to the chain, for
the following controller implementations:

  * •

“nominal”: Using only the nominal model and no uncertainty description,

  * •

“naïve”: Solving (8) exactly without any data points (no GP-related
computations444 Due to software limitations, for a fair comparison the “naïve”
GP-MPC implementation only makes use of the GP prior, i.e., $\mu^{d}(y)\equiv
0$ and $\Sigma^{d}(y)\equiv const.$; it corresponds to a stochastic NMPC
implementation with the process noise covariance inflated by the GP prior
covariance. Therefore the present timings can be viewed as a lower bound for
the actual timings when solving (8) with a GP conditioned on data, which would
additionally require not only GP posterior evaluations and derivative
computations, but also expensive computations associated with the Hessian of
the GP posterior mean. ), including the covariances as optimization variables
into an augmented state,

  * •

“alg2-cpu$C$-$D$”: Alg. 2 with $D$ data points and GP inference and AD on $C$
CPU cores, and

  * •

“alg2-gpu-$D$”: Alg. 2 with $D$ data points and GP inference and AD performed
on GPU.

It becomes evident that the scaling of order ${O}(n_{x}^{6})$ for the “naïve”
method quickly becomes prohibitive in terms of computation time. In
comparison, for $n_{x}=39$, a roughly $1000$-fold speed-up by means of the
zero-order optimization strategy can be observed. Comparing the computation
times of Alg. 2 for different number of data points clearly shows the effects
of GPU acceleration, which accounts for a $5-10$-fold speed-up and better
scaling properties per SQP iteration for $n_{x}\geq 33$ and the associated GP
output dimension.

Figure 2: SQP timings comparison for increasing number of states $n_{x}$ and
GP dimension $n_{w}=(n_{x}-3)/2$.

###  V-C Timing profile of Alg. 2

Fig. 3 displays the runtime profile of Alg. 2 for $n_{x}=33$ states, divided
into the following components:

  * •

“acd_itf”: acados interface for get/set operations,

  * •

“acd_itg”: acados implicit Runge-Kutta integrator,

  * •

“acd_qp”: acados QP solve of (14) using HPIPM [26],

  * •

“con_tight”: NumPy covariance propagation (6) and constraint tightening for
(14c),

  * •

“gpytorch”: GPyTorch GP inference and AD.

For 150 data points, Fig. 3 shows that the overhead introduced by
parallelizing gpytorch operations is higher than the computational speed-ups,
especially on the GPU. For 1500 data points, however, the parallelization
leads to significantly lower computation times, saving about 60-70% of the
computational costs associated with the GPs, which at this point dominate the
computational footprint of the method. For 0 data points, potential for
improvement can be seen in terms of the constraint tightening computations,
which could be improved by switching to a C/C++ implementation.

Figure 3: Timing profile for Alg. 2 variants for $n_{x}=33$
($n_{\text{mass}}=7$).

##  VI Conclusions

To make GP-MPC computationally tractable, both the complexity to solve the
OCP, as well as evaluating and differentiating the GP posterior mean and
covariance need to be addressed. To tackle these challenges, this paper has
presented an inexact SQP approach with a tailored Jacobian approximation,
while parallelizing GP inference and differentiation routines on a GPU. The
results demonstrate that drastically reduced computation times can be achieved
while ensuring feasibility of the converged iterates and maintaining favorable
convergence properties, pushing the computational boundaries to apply GP-MPC
in real-world scenarios.

### Acknowledgments

We would like to thank Katrin Baumgärtner and Jonathan Frey for pointing us
towards the acados Cython interface, and Johannes Köhler for many helpful
discussions.

## References

  * [1] C. J. Ostafew, A. P. Schoellig, and T. D. Barfoot, “Learning-based nonlinear model predictive control to improve vision-based mobile robot path-tracking in challenging outdoor environments,” in _2014 IEEE International Conference on Robotics and Automation (ICRA)_ , 2014, pp. 4029–4036.
  * [2] A. Carron, E. Arcari, M. Wermelinger, L. Hewing, M. Hutter, and M. N. Zeilinger, “Data-Driven Model Predictive Control for Trajectory Tracking With a Robotic Arm,” _IEEE Robotics and Automation Letters_ , vol. 4, no. 4, pp. 3758–3765, 2019.
  * [3] J. Kabzan, L. Hewing, A. Liniger, and M. N. Zeilinger, “Learning-Based Model Predictive Control for Autonomous Racing,” _IEEE Robot. Autom. Lett._ , vol. 4, no. 4, pp. 3363–3370, 2019.
  * [4] J. Kabzan, M. Valls, V. Reijgwart, H. Hendrikx, C. Ehmke, M. Prajapat, A. Bühler, N. Gosala, M. Gupta, R. Sivanesan, A. Dhall, E. Chisari, N. Karnchanachari, S. Brits, M. Dangel, I. Sa, R. Dube, A. Gawel, M. Pfeiffer, and R. Siegwart, “AMZ Driverless: The full autonomous racing system,” _Journal of Field Robotics_ , 2020.
  * [5] G. Torrente, E. Kaufmann, P. Föhn, and D. Scaramuzza, “Data-Driven MPC for Quadrotors,” _IEEE Robotics and Automation Letters_ , vol. 6, no. 2, pp. 3769–3776, 2021.
  * [6] L. Hewing, J. Kabzan, and M. N. Zeilinger, “Cautious Model Predictive Control Using Gaussian Process Regression,” _IEEE Trans Control Syst Technol_ , vol. 28, no. 6, pp. 2736–2743, 2020.
  * [7] L. P. Fröhlich, C. Küttel, E. Arcari, L. Hewing, M. N. Zeilinger, and A. Carron, “Model Learning and Contextual Controller Tuning for Autonomous Racing,” _2022 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), to appear_ , 2022.
  * [8] S. Vaskov, R. Quirynen, M. Menner, and K. Berntorp, “Friction-Adaptive Stochastic Predictive Control for Trajectory Tracking of Autonomous Vehicles,” in _2022 American Control Conference (ACC)_ , 2022, pp. 1970–1975.
  * [9] X. Feng, S. D. Cairano, and R. Quirynen, “Inexact Adjoint-based SQP Algorithm for Real-Time Stochastic Nonlinear MPC,” _IFAC-PapersOnLine_ , vol. 53, no. 2, pp. 6529–6535, 2020.
  * [10] A. Zanelli, J. Frey, F. Messerer, and M. Diehl, “Zero-Order Robust Nonlinear Model Predictive Control with Ellipsoidal Uncertainty Sets,” _IFAC-PapersOnLine_ , vol. 54, no. 6, pp. 50–57, 2021.
  * [11] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, A. Desmaison, A. Kopf, E. Yang, Z. DeVito, M. Raison, A. Tejani, S. Chilamkurthy, B. Steiner, L. Fang, J. Bai, and S. Chintala, “PyTorch: An Imperative Style, High-Performance Deep Learning Library,” in _Advances in Neural Information Processing Systems 32_. Curran Associates, Inc., 2019, pp. 8024–8035.
  * [12] J. Gardner, G. Pleiss, K. Q. Weinberger, D. Bindel, and A. G. Wilson, “GPyTorch: Blackbox Matrix-Matrix Gaussian Process Inference with GPU Acceleration,” in _Advances in Neural Information Processing Systems_ , vol. 31. Curran Associates, Inc., 2018.
  * [13] F. Messerer and M. Diehl, “An Efficient Algorithm for Tube-based Robust Nonlinear Optimal Control with Optimal Linear Feedback,” in _2021 60th IEEE Conference on Decision and Control (CDC)_ , 2021, pp. 6714–6721.
  * [14] R. Quirynen and K. Berntorp, “Uncertainty Propagation by Linear Regression Kalman Filters for Stochastic NMPC,” _IFAC-PapersOnLine_ , vol. 54, no. 6, pp. 76–82, 2021.
  * [15] C. E. Rasmussen and C. K. I. Williams, _Gaussian Processes for Machine Learning_ , ser. Adaptive Computation and Machine Learning. Cambridge, Massachusetts: MIT Press, 2006.
  * [16] A. Girard, C. E. Rasmussen, and R. Murray-Smith, “Gaussian Process priors with Uncertain Inputs: Multiple-Step-Ahead Prediction,” Department of Computing Science, University of Glasgow, Technical Report TR-2002-119, 2002.
  * [17] L. Hewing, A. Liniger, and M. N. Zeilinger, “Cautious NMPC with Gaussian Process Dynamics for Autonomous Miniature Race Cars,” in _2018 European Control Conference (ECC)_ , 2018, pp. 1341–1348.
  * [18] H. G. Bock and K. J. Plitt, “A Multiple Shooting Algorithm for Direct Solution of Optimal Control Problems*,” _IFAC Proceedings Volumes_ , vol. 17, no. 2, pp. 1603–1608, 1984.
  * [19] D. Kouzoupis, G. Frison, A. Zanelli, and M. Diehl, “Recent Advances in Quadratic Programming Algorithms for Nonlinear Model Predictive Control,” _Vietnam J. Math._ , vol. 46, no. 4, pp. 863–882, 2018.
  * [20] P. T. Boggs and J. W. Tolle, “Sequential quadratic programming,” _Acta Numerica_ , vol. 4, pp. 1–51, 1995.
  * [21] H. G. Bock, M. Diehl, E. Kostina, and J. P. Schlöder, “1. Constrained Optimal Feedback Control of Systems Governed by Large Differential Algebraic Equations,” in _Real-Time PDE-Constrained Optimization_. Society for Industrial and Applied Mathematics, 2007, pp. 3–24.
  * [22] S. M. Robinson, “Strongly Regular Generalized Equations,” _Math. Oper. Res._ , vol. 5, no. 1, pp. 43–62, 1980.
  * [23] A. Zanelli, Q. Tran-Dinh, and M. Diehl, “Contraction Estimates for Abstract Real-Time Algorithms for NMPC,” in _2019 IEEE 58th Conference on Decision and Control (CDC)_ , 2019, pp. 8085–8092.
  * [24] A. L. Dontchev and R. T. Rockafellar, _Implicit Functions and Solution Mappings: A View from Variational Analysis_ , ser. Springer Monographs in Mathematics. New York, NY: Springer, 2009.
  * [25] J. A. E. Andersson, J. Gillis, G. Horn, J. B. Rawlings, and M. Diehl, “CasADi: A software framework for nonlinear optimization and optimal control,” _Math. Prog. Comp._ , vol. 11, no. 1, pp. 1–36, 2019.
  * [26] G. Frison and M. Diehl, “HPIPM: A high-performance quadratic programming framework for model predictive control,” _IFAC-PapersOnLine_ , vol. 53, no. 2, pp. 6563–6569, 2020.

###  -A Explicit form of (12)

By linearity of $g(y,P)$ in $P$, we can show that $\frac{\partial g}{\partial
P}(y,P)$ is indeed invertible for our problem setting. It holds that

$\displaystyle A(y)$ $\displaystyle:=\begin{bmatrix}I_{n_{x}^{2}}&&\\\
L(y_{0})&I_{n_{x}^{2}}&\\\ &\ddots&\ddots&\\\ &&L(y_{N-1})&I_{n_{x}^{2}}\\\
\end{bmatrix}$ (30) and $\displaystyle b(y)$
$\displaystyle:=\begin{bmatrix}0\\\ (B\otimes
B)\overrightarrow{\left(\Sigma^{d}(y_{0})+\Sigma^{w}\right)}\\\ \vdots\\\
(B\otimes
B)\overrightarrow{\left(\Sigma^{d}(y_{N-1})+\Sigma^{w}\right)}\end{bmatrix},$
(31)

where $L(y_{i}):=-\tilde{A}_{i}(y_{i})\otimes\tilde{A}_{i}(y_{i})$ and
“$\otimes$” denotes the Kronecker product. Differentiating equation (12), we
obtain that $\frac{\partial g}{\partial P}(y)=A(y)$; invertibility of
$\frac{\partial g}{\partial P}$ follows since $A(y)$ is lower triangular.