MIDLMedical Imaging with Deep Learning 2022 Full Paper – MIDL 2022 submission
– Under Review Anna Yeaton1<EMAIL_ADDRESS>
Rahul G. Krishnan4<EMAIL_ADDRESS>
Rebecca Mieloszyk2<EMAIL_ADDRESS>
David Alvarez-Melis3<EMAIL_ADDRESS>
Grace Huynh2<EMAIL_ADDRESS>
1 Department of Pathology, NYU Grossman School of Medicine
4 Department of Computer Science, University of Toronto
2 Microsoft Research, Redmond
3 Microsoft Research, New England

# Hierarchical Optimal Transport for Comparing Histopathology Datasets

###### Abstract

Scarcity of labeled histopathology data limits the applicability of deep
learning methods to under-profiled cancer types and labels. Transfer learning
allows researchers to overcome the limitations of small datasets by pre-
training machine learning models on larger datasets _similar_ to the small
target dataset. However, similarity between datasets is often determined
heuristically. In this paper, we propose a principled notion of distance
between histopathology datasets based on a hierarchical generalization of
optimal transport distances. Our method does not require any training, is
agnostic to model type, and preserves much of the hierarchical structure in
histopathology datasets imposed by tiling. We apply our method to H&E stained
slides from The Cancer Genome Atlas from six different cancer types. We show
that our method outperforms a baseline distance in a cancer-type prediction
task. Our results also show that our optimal transport distance predicts
difficulty of transferability in a tumor vs. normal prediction setting.

###### keywords:

Histopathology, domain adaptation, transfer learning, optimal transport

††editors: Under Review for MIDL 2022

##  1 Introduction

Histopathology images are routinely used in the diagnostic workup of many
cancers. Beyond the standard identification of tumor grade and subtype,
histopathology images also contain an abundance of visual information that may
have predictive and prognostic value. Recent advances in deep learning for
histopathology are facilitating the extraction of this information, enabling
prediction of genetic alterations, treatment response and survival (Coudray et
al., 2018; Muhammad et al., 2021; Wulczyn et al., 2021; Echle et al., 2020).

Despite the promise of supervised deep learning, the large, labeled datasets
required to train complex networks on histopathology images are scarce,
especially in less common cancers and label types. To overcome data scarcity
during model training, transfer learning techniques can be utilized by pre-
training models on a larger—ideally similar—dataset. The model can then be
fine-tuned on the small target dataset of interest. Until now, determining
which similar dataset to use for pre-training in histopathology has been
driven by intuition or trial and error (Srinidhi et al., 2021).

In this paper, we introduce a principled approach to measuring similarities
between histopathology datasets. Specifically, we propose a novel distance
between histopathology datasets which we call Hierarchical Histopathology
Optimal Transport (HHOT) based on Optimal Transport (OT), a method to compare
distributions. Our method uses OT to compute the distance between individual
slides and to compute the distances between entire datasets (fig:Figure1). We
show that HHOT is highly predictive of transfer learning accuracy in a tumor
vs. normal prediction setting, and that it is significantly faster than a
naive (non-hierarchical) Optimal Transport approach.

fig:Figure1 [For a pair of slides, the OT distance is calculated using the OT
distance between all pairs of tiles. ][centred]  [For a pair of datasets, the
OT distance is calculated using the OT distance between all pairs of
slides][centred]

Figure 1: Schematic of how HHOT distances are calculated at the slide (a) and
dataset (b) level. Tiling of individual slides, done to overcome memory
limits, introduces hierarchical structure to the calculation of OT between
different datasets.

##  2 Background

Optimal Transport (OT) is a mathematical framework centered around the goal of
comparing probability distributions, with deep theory (Villani, 2003, 2008;
Santambrogio, 2015) and applications to various fields, ranging from economics
(Galichon, 2016) to meteorology (Cullen and Maroofi, 2003). Although it can be
formulated in more general settings, in this work we are interested in its
discrete Euclidean formulation, which considers two finite collections of
points $\\{\textbf{x}^{(i)}\\}_{i=1}^{n}$, $\\{\textbf{y}^{(j)}\\}_{j=1}^{m}$,
$\textbf{x}^{(i)},\textbf{y}^{(j)}\in\mathbb{R}^{d}$, represented as empirical
distributions: $\mu=\sum_{i=1}^{n}\textbf{p}_{i}\delta_{\textbf{x}^{(i)}}$,
$\nu=\sum_{j=1}^{m}\textbf{q}_{j}\delta_{\textbf{y}^{(j)}}$, where p and q are
probability vectors (non-negative and adding up to one).

At a high level, the goal of OT is to find an optimal correspondence between
these distributions, and in doing so, define a notion of similarity between
them. Given a cost function (often called the ground metric) between pairs of
points $c:\mathbb{R}^{d}\times\mathbb{R}^{d}\rightarrow\mathbb{R}+$, the goal
of OT is to find a correspondence between $\mu$ and $\nu$ with minimal cost.
Formally, the Kantorovich formulation of discrete OT seeks a coupling matrix
$\Gamma\in\mathbb{R}^{n\times m}$ that solves:

$\mathrm{OT}_{c}(\mu,\nu)\stackrel{{\scriptstyle\mathclap{\tiny\mbox{def.}}}}{{=}}\min_{\Gamma\in\Pi(\mu,\nu)}\langle\Gamma,\textbf{C}\rangle=\sum_{ij}\Gamma_{ij}\textbf{C}_{ij},$
(1)

where
$\textbf{C}_{ij}\stackrel{{\scriptstyle\mathclap{\tiny\mbox{def.}}}}{{=}}c(\textbf{x}^{(i)},\textbf{y}^{(j)})$.
The constraint set $\Pi(\mu,\nu)$ enforces $\Gamma$ to be measure-preserving,
i.e., to have $\mu$ and $\nu$ as its marginals:

$\Pi(\mu,\nu)\stackrel{{\scriptstyle\mathclap{\tiny\mbox{def.}}}}{{=}}\\{\Gamma\in\mathbb{R}^{n\times
m}_{+}\medspace|\medspace\Gamma\textbf{1}=\textbf{p},\medspace\Gamma^{\top}\textbf{1}=\textbf{q}\\}.$
(2)

It can be shown that for
$c(\textbf{x},\textbf{y})=\|\textbf{x}-\textbf{y}\|^{p}$,
$\mathrm{OT}_{c}(\mu,\nu)^{1/p}$ is a proper distance metric between
distributions (i.e., satisfies all metric axioms) (Peyré and Cuturi, 2019). As
noted above, the coupling matrix $\Gamma$ can be interpreted as a soft
matching or probabilistic correspondence between the elements of $\mu$ and
$\nu$, in the sense that $\Gamma_{ij}$ is high if
$\textbf{x}^{(i)},\textbf{y}^{(j)}$ are in ‘correspondence’, and low
otherwise. In some cases (such as the common case $n=m$ and
$\textbf{p},\textbf{q}$ uniform), the optimal coupling turns out to be sparse,
in which case $\Gamma$ defines a deterministic one-to-one mapping between the
$\textbf{x}^{(i)}$’s and $\textbf{y}^{(j)}$’s.

Problem (1) is a linear programming problem and thus solvable exactly in
$O(n^{3})$ time (Peyré and Cuturi, 2019). This makes it impractical even for
moderately sized collections of points. However, seminal work by Cuturi (2013)
showed that a regularized version of this problem can be solved much more
efficiently. The regularization consists of adding a entropy term on the
objective:

$\mathrm{OT}_{c,\varepsilon}(\mu,\nu)\stackrel{{\scriptstyle\mathclap{\tiny\mbox{def.}}}}{{=}}\min_{\Gamma\in\Pi(\mu,\nu)}\langle\Gamma,\textbf{C}\rangle+\varepsilon\mathrm{H}(\Gamma).$
(3)

This entropy-regularized OT problem can be solved very efficiently using the
Sinkhorn-Knopp algorithm (Cuturi, 2013). One downside of this regularization
is that $\mathrm{OT}_{c,\varepsilon}(\mu,\mu)\neq 0$, which in particular
implies this quantity is no longer a valid distance. To alleviate this, prior
work (Genevay et al., 2018; Salimans et al., 2018) has considered a debiased
version this quantity, also known as the Sinkhorn divergence:

$\mathrm{SD}_{c,\varepsilon}(\mu,\nu)\stackrel{{\scriptstyle\mathclap{\tiny\mbox{def.}}}}{{=}}\mathrm{OT}_{c,\varepsilon}(\mu,\nu)-\tfrac{1}{2}\bigl{(}\mathrm{OT}_{c,\varepsilon}(\mu,\mu)+\mathrm{OT}_{c,\varepsilon}(\nu,\nu)\bigr{)}$
(4)

In addition to satisfying $\mathrm{SD}_{c,\varepsilon}(\mu,\nu)\geq 0$, with
equality if and only if $\mu=\nu$, this divergence comes with many other
desirable theoretical properties: it is positive, convex, metrizes weak
converge in distribution (Feydy et al., 2019), and leads to faster statistical
rates of estimation of the exact OT problem (Chizat et al., 2020).

In Section 4.1, we will introduce our method using
$\mathrm{OT}_{c,\varepsilon}(\mu,\nu)$ for notational simplicity, noting that
it can naturally use $\mathrm{SD}_{c,\varepsilon}(\mu,\nu)$ instead. In
addition, we will drop $c$ and $\epsilon$ from the notation for $\mathrm{OT}$
when these are clear from the context.

##  3 Related work

Our HHOT method for histopathology builds on previous work in hierarchical OT
for other domains and previous work in non-hierarchical OT within
histopathology. For example, Yurochkin et al. (2019) have described a
hierarchical OT method for measuring distances between documents, using words
and topics as the hierarchical levels. In addition, specifically for
histopathology images, non-hierarchical OT has been used to compare individual
cell morphology (Basu et al., 2014; Wang et al., 2011) and to quantify domain
shift at the tile level (Stacke et al., 2021). In addition, the relationship
between OT-calculated dataset distances and the difficulty of transferability
has been previously described by (Alvarez-Melis and Fusi, 2020; Gao and
Chaudhari, 2021; Achille et al., 2021), although the notion of distance they
define is generic and thus does not leverage the hierarchical nature of
individual datasets.

##  4 Optimal transport between histopathology datasets

###  4.1 Hierarchical Histopathology Optimal Transport

We consider a pair of histopathology datasets, $\mathsf{D}_{a}$ and
$\mathsf{D}_{b}$, collected from different tissues, centers, or populations.
We seek a notion of distance that lets us assess their similarity. Each
dataset consists of slides s, which in turn are subdivided into tiles t. Let
$n$ and $m$ denote, respectively, the number of slides in each dataset. In
addition, we denote by $n_{i}$ the number of tiles in the $i$-th slide of the
first dataset, and analogously for $m_{j}$.

We can view $\mathsf{D}_{a}$ and $\mathsf{D}_{b}$ as point clouds or, more
formally, empirical distributions as described in section 2. From this
viewpoint, we can think of these datasets as samples of slide images sampled
from two different underlying distributions. Unless additional information is
provided, we can simply take the weights associated to each slide (p and q) to
be uniform, as is typically done in practical application of OT to point
clouds or images (Peyré and Cuturi, 2019). After defining a suitable notion of
distance between pairs of slides, one could in principle use problem (3) (or
its debiased counterpart, eq. (4)) to obtain a notion of similarity between
$\mathrm{D}_{a}$ and $\mathrm{D}_{b}$. However, this naive approach would
require loading entire slides into memory, which is computationally
infeasible—precisely the reason why these images are tiled in the first place.

Our proposed solution to this computational hurdle is to interpret the slides
themselves as collections (formally, distributions) of tiles, and use OT once
more, now to define a notion of distance between these. To this end, we first
define the cost between individual tiles as the distance between them.
Although we could in principle use the Euclidean distance between their raw
pixel representations, a more meaningful comparison can be obtained by first
embedding these images in some lower dimensional space (e.g., using a neural
network pre-trained on a large image dataset), and then computing a distance
between them. Hence, we define
$c(\textbf{t}_{u},\textbf{t}_{v})=\|\phi(\textbf{t}_{u})-\phi(\textbf{t}_{v})\|^{2}$,
where $\phi$ is a pre-trained encoder.

We collect all such pairwise costs in a matrix $\textbf{C}_{\text{tile}}$ of
size $n_{i}\times m_{j}$, and solve the corresponding OT problem:

$\mathrm{OT}_{\phi,\varepsilon}(\textbf{s}_{i},\textbf{s}_{j})=\min_{\Gamma\in\Pi(\textbf{p}_{i},\textbf{q}_{j})}\langle\Gamma_{\text{tile}},\textbf{C}^{ij}_{\text{tile}}\rangle+\varepsilon\textup{H}(\Gamma).$
(5)

As in section 2, we can also define a debiased version of this slide-to-slide
distance:

$\mathrm{SD}_{\phi,\varepsilon}(\textbf{s}_{i},\textbf{s}_{j})=\mathrm{OT}_{\phi,\varepsilon}(\textbf{s}_{i},\textbf{s}_{j})-\tfrac{1}{2}\bigl{(}\mathrm{OT}_{\phi,\varepsilon}(\textbf{s}_{i},\textbf{s}_{i})+\mathrm{OT}_{\phi,\varepsilon}(\textbf{s}_{j},\textbf{s}_{j})\bigr{)}$
(6)

Compared to other possible ways to compare slides using their tiles (like
using a mean or centroid tile), this OT-based approach is appealing because
(i) it does not lose information by aggregating the tiles, and (ii) it
recovers some of the global structure that it lost by tiling, i.e., the
relation between the tiles in the context of the slide. Specifically, in
operationalizing similarity through matching, problem (5) seeks to find
corresponding tiles across slides, and does it coherently as a result of the
marginal constraint (e.g., a single tile cannot be matched to all tiles in the
other slide).

Once we have computed (5) for every pair of slides
$(\textbf{s}_{i},\textbf{s}_{j})$ from the two datasets, we collect them in a
matrix $\textbf{C}_{\text{slide}}$ with entries
$[\textbf{C}_{\text{slide}}]_{ij}=\mathrm{OT}_{\varepsilon}(\textbf{s}_{i},\textbf{s}_{j})$.
With this, finally we have a ground cost between slides, which we can use to
compute the sought-after distance between datasets using OT once more:

$\mathrm{OT}_{\varepsilon}(\mathsf{D}_{a},\mathsf{D}_{b})=\min_{\Gamma\in\Pi(\textbf{p},\textbf{q})}\langle\Gamma,\textbf{C}_{\text{slide}}\rangle+\varepsilon\textup{H}(\Gamma),$
(7)

and its debiased counterpart:

$\mathrm{SD}_{\phi,\varepsilon}(\mathsf{D}_{a},\mathsf{D}_{b})=\mathrm{OT}_{\phi,\varepsilon}(\mathsf{D}_{a},\mathsf{D}_{b})-\tfrac{1}{2}\bigl{(}\mathrm{OT}_{\phi,\varepsilon}(\mathsf{D}_{a},\mathsf{D}_{b})+\mathrm{OT}_{\phi,\varepsilon}(\mathsf{D}_{a},\mathsf{D}_{b})\bigr{)}.$
(8)

###  4.2 Computational Implementation

We use the python optimal transport (POT) (Flamary et al., 2021) and geomloss
(Feydy et al., 2019) libraries to solve the individual OT problems (5) and
(7). Using the vanilla Sinkhorn algorithm, solving the first of these to
$\delta_{1}$-accuracy has $O(n_{i}m_{j}/\delta_{1})$ computational complexity
(Altschuler et al., 2017), and analogously $O(nm/\delta_{2})$ for the latter.
Taking $\delta_{1}=\delta_{2}$, the total complexity scales as
$O((nm+\sum_{ij}n_{i}m_{j})/\delta)$. fig:Figure4c shows empirical runtimes
for our method. Our implementation of HHOT can be found here:
https://github.com/ayeaton/HHOT

##  5 Methodology and Results

fig:Figure2 [Source tile visualizing cartilage within a LUAD slide.]  [Target
slide of LUAD. The orange highlighted region is cartilaginous.]  [Heatmap of
coupling between the source tile and target slide. High coupling is observed
in the cartilaginous region.]

Figure 2: Representative example of a source tile tightly coupling to target
tiles of the same tissue type, suggesting that HHOT between tiles implicitly
incorporates biological information.

We used whole slide images retrieved from the TCGA
(https://portal.gdc.cancer.gov/) for six common cancer types, including both
primary tumor samples and matched normal samples. Specifically we focused on
Stomach Adenocarcinoma (STAD), Bladder Urothelial Carcinoma (BLCA), Lung
Adenocarcinoma (LUAD), Lung Squamous-cell Carcinoma (LUSC), Colon
Adenocarcinoma (COAD), and Pancreatic Adenocarcinoma (PAAD). Slide images were
tiled into 512x512 pixel non-overlapping images at 20x magnification, and
tiles with more than 50 percent background were discarded as described in
Coudray et al. (2018); Noorbakhsh et al. (2020). We extracted a median of 474
tiles for BLCA, 777 for COAD, 445 for LUAD, 448 for LUSC, 278 for PAAD and 622
for STAD. We then used Inception-V3 pre-trained on ImageNet to compress
512x512 pixel images to 2048 length vectors. For comparison of collections of
tiles, we set the regularization parameter $\varepsilon$ to 0.25.

We calculated OT between slides, and found that tiles with the highest
coupling were those that were visually similar. As an example, we show the OT
coupling matrix between tiles from two representative slides in fig:Figure2.
We chose an example tile from the source slide which displays cartilage
tissue. We then show the coupling of this source tile to all tiles in the
target slide. We visualize the strength of coupling in fig:Figure2c; the
highest coupling to the cartilage source tile is tightly localized and
corresponds to a cartilaginous region in the target. These results demonstrate
that an OT solution implicitly incorporates biological structure.

fig:Figure3  [Centroid UMAP.] [OT UMAP.] [KNN performance. ] [HHOT distance
retains visual diversity within cancer datasets. ]

Figure 3: HHOT enables better clustering of slide images by cancer type, as
compared to the centroid tile distance method. HHOT also reflects the
diversity of images within a single cancer type. This suggests that HHOT
distances better preserve morphological information.

In slide-to-slide comparisons, HHOT preserves more relevant histomorphological
information than other methods, such as using a mean or centroid tile. We
demonstrate this with a task for clustering the slides by cancer type. We
created a matrix $\textbf{C}_{\text{slide}}$ with entries
$[\textbf{C}_{\text{slide}}]_{ij}=\mathrm{OT}_{\varepsilon}(\textbf{s}_{i},\textbf{s}_{j})$
as described above including all pairs of slides. For reference, we also
calculated a distance matrix between slides using the centroid tile of a slide
as described in Howard et al. (2021). We visualized the relationship between
slides using UMAP for both OT and centroid-tile distance (fig:Figure3a,
fig:Figure3b). Visually, we observe that HHOT distance enables better
clustering of the slide images by cancer type. Quantitatively, using K-nearest
neighbors over K ranging from two to eight, we performed a cancer-type
classification task with our two distance matrices as input. We show that OT
retains similarity within cancer-types and dissimilarity across cancer-types
as expected, and does so better than the centroid-tile distance baseline.
fig:Figure3e shows representative examples of how HHOT distances reflect the
diversity of images with a single cancer type. HHOT tightly groups PAAD
slides, and these slide images are visually very similar. In contrast, HHOT
shows more variability in the LUSC slides, and these images are more variable.

We observed that target datasets that are similar to the source dataset based
on HHOT distance are more improved by pre-training, i.e. a negative
correlation between HHOT distance and transferability. In this paper, we focus
on a tumor vs. normal task. For our pre-training tasks, we standardized the
total dataset size to 209 slides, and used a cross validation scheme to create
four datasets per cancer type with 169 slides for training and 40 for
validation. For our fine-tuning task, we standardized the dataset to 65
slides, and used a cross validation scheme to create four datasets per cancer
type with 25 slides for training and 40 for validation. We conducted our task
over all the pre-training datasets and all the fine-tuning datasets for 16
experiments per cancer type comparison. For each dataset pair, we pre-trained
a single-layer perceptron to predict tumor or normal status, using feature
vectors $\phi(\textbf{t})$ and a learning rate of 1e-2. We then fine-tuned
this pre-trained model for each of the other cancer types, using only 25
target slides and a learning rate of 1e-10.

We quantify transferability across tasks using the relative improvement in AUC
obtained by transfer learning:
$\textsf{Transferability}(D_{T},D_{S})=(\textsf{AUC}(D_{S}\rightarrow
D_{T})-\textsf{AUC}(D_{T}))/\textsf{AUC}(D_{T})$. We observed a negative
correlation between the HHOT distance and transferability most strongly for
PAAD, LUAD, LUSC, and STAD (fig:Figure4a). Consistently, visually similar
datasets, such as the two lung cancer datasets, show the smallest HHOT
distance and the highest transferability (green and purple squares and plus in
fig:Figure4a). For BLCA and COAD, we observe that the AUC of a model without
pre-training is already very high (fig:Figure4b); thus, pre-training on other
datasets cannot improve on this already high AUC. This is explicitly
quantified for these two cancers by the nearly zero slope of the best fit line
(fig:Figure4b, boxplot of no-pretaining AUCs aggregating four datasets of 25
target slides).

fig:example2

Figure 4: HHOT is predictive of model transferability across cancer types.
Datasets with large HHOT distances, representing more visually different data,
have worse transferability. The regression coefficients are [$\beta$,
$\beta_{0}$] [0.0024, $-5.47\text{\times}{10}^{-5}$] (COAD), [0.0016,
$-6\text{\times}{10}^{-5}$] (BLCA), [0.022, $-6.8\text{\times}{10}^{-4}$]
(PAAD), [0.042, $-1.7\text{\times}{10}^{-3}$] (LUAD), [0.043,
$-2\text{\times}{10}^{-3}$] (LUSC), [0.049, $-2.7\text{\times}{10}^{-3}$]
(STAD). For COAD and BLCA, the baseline AUC is already quite high (b), so
transfer learning with other cancer types results in little improvement. HHOT
is faster than non-hierarchical OT between slides (c).

[HHOT vs. transferability. ]  [No-transfer AUC.] [HHOT reduces runtime.]

Finally, we show that not only does HHOT preserve the hierarchical structure
of the data and correlate with difficulty of transfer learning, it is also
much faster than a flat (non-hierarchical) OT approach. We compare the time it
took to compute OT distance between four to eighteen slides, with 100 tiles
each and observed that the time it took to calculate OT distance between
slides increased in a polynomial order in the flat case (Figure 4).

##  6 Discussion

In this paper we introduced a principled approach to compare histopathology
datasets called HHOT. Our work adds to the OT and histopathology literature by
proposing a method to compare datasets while also preserving the structure
lost by standard tiling approaches. Specifically, we propose to solve first an
inner tile-to-tile OT problem for all pairs of slides, and then solve the
outer slide-to-slide OT problem between datasets. We first show that
correspondences in the OT coupling matrix $\Gamma_{tile}$ map the same type of
tissue across slides. We then show that OT distance performs better than a
naive, centroid-tile distance based method in a cancer-type prediction task.
We also show that HHOT distance between histopathology datasets correlates
with transferability. Furthermore, we find that the degree of correlation
between HHOT distance and transferability is associated with baseline AUC.
Simply, if the task is already close to optimal, it is difficult to improve
using pre-training. Finally, we show that HHOT distance is much faster than a
naive, flat approach to comparing histopathology datasets. In addition to
applications presented in this paper, promising applications of HHOT include
outlier detection, clustering analysis, dataset visualization, and
facilitating multi-modal integration of whole slide images and molecular data.

In conclusion, we have demonstrated that HHOT has many benefits for
researchers working with tiled histopathology data. Future work may focus on
using OT to direct dataset creation to optimize transferability, tune tile
size and normalization hyper-parameters, and compare feature vectors learned
from different models. R.G.K., R.M., D.A.M, and G.H were employed by Microsoft
corporation while performing this work. A.Y. was employed by Microsoft
corporation for a duration of this work. Part of this work has used computing
resources at the NYU School of Medicine High Performance Computing Facility.

## References

  * Achille et al. (2021) Alessandro Achille, Giovanni Paolini, Glen Mbeng, and Stefano Soatto.  The information complexity of learning tasks, their structure and their distance.  _Information and Inference: A Journal of the IMA_ , 10(1):51–72, 01 2021.  ISSN 2049-8772.  10.1093/imaiai/iaaa033.  URL https://doi.org/10.1093/imaiai/iaaa033.
  * Altschuler et al. (2017) Jason Altschuler, Jonathan Niles-Weed, and Philippe Rigollet.  Near-linear time approximation algorithms for optimal transport via sinkhorn iteration.  In I Guyon, U V Luxburg, S Bengio, H Wallach, R Fergus, S Vishwanathan, and R Garnett, editors, _Advances in Neural Information Processing Systems 30_ , pages 1964–1974. Curran Associates, Inc., 2017.
  * Alvarez-Melis and Fusi (2020) David Alvarez-Melis and Nicolo Fusi.  Geometric dataset distances via optimal transport.  In H. Larochelle, M. Ranzato, R. Hadsell, M. F. Balcan, and H. Lin, editors, _Advances in Neural Information Processing Systems_ , volume 33, pages 21428–21439. Curran Associates, Inc., 2020.  URL https://proceedings.neurips.cc/paper/2020/file/f52a7b2610fb4d3f74b4106fb80b233d-Paper.pdf.
  * Basu et al. (2014) Saurav Basu, Soheil Kolouri, and Gustavo K. Rohde.  Detecting and visualizing cell phenotype differences from microscopy images using transport-based morphometry.  _Proceedings of the National Academy of Sciences_ , 111(9):3448–3453, 2014.  ISSN 0027-8424.  10.1073/pnas.1319779111.  URL https://www.pnas.org/content/111/9/3448.
  * Chizat et al. (2020) Lenaic Chizat, Pierre Roussillon, Flavien Léger, François-Xavier Vialard, and Gabriel Peyré.  Faster wasserstein distance estimation with the sinkhorn divergence.  In _NeurIPS_ , 2020.  URL https://proceedings.neurips.cc/paper/2020/hash/17f98ddf040204eda0af36a108cbdea4-Abstract.html.
  * Coudray et al. (2018) Nicolas Coudray, Paolo Santiago Ocampo, Theodore Sakellaropoulos, Navneet Narula, Matija Snuderl, David Fenyö, Andre L. Moreira, Narges Razavian, and Aristotelis Tsirigos.  Classification and mutation prediction from non–small cell lung cancer histopathology images using deep learning.  _Nature Medicine_ , 2018.
  * Cullen and Maroofi (2003) Mike Cullen and Hamed Maroofi.  The fully compressible Semi-Geostrophic system from meteorology.  _Arch. Ration. Mech. Anal._ , 167(4):309–336, May 2003.  ISSN 0003-9527, 1432-0673.  10.1007/s00205-003-0245-x.
  * Cuturi (2013) Marco Cuturi.  Sinkhorn distances: Lightspeed computation of optimal transport.  In C J C Burges, L Bottou, M Welling, Z Ghahramani, and K Q Weinberger, editors, _Advances in Neural Information Processing Systems 26_ , pages 2292–2300. Curran Associates, Inc., 2013.  URL https://proceedings.neurips.cc/paper/2013/file/af21d0c97db2e27e13572cbf59eb343d-Paper.pdf.
  * Echle et al. (2020) Amelie Echle, Niklas Rindtorff, Titus Brinker, Tom Luedde, Alexander Pearson, and Jakob Kather.  Deep learning in cancer pathology: a new generation of clinical biomarkers.  _British Journal of Cancer_ , 124:1–11, 11 2020.  10.1038/s41416-020-01122-x.
  * Feydy et al. (2019) Jean Feydy, Thibault Séjourné, François-Xavier Vialard, Shun-ichi Amari, Alain Trouvé, and Gabriel Peyré.  Interpolating between optimal transport and mmd using sinkhorn divergences.  In _AISTATS_ , pages 2681–2690. PMLR, 2019.  URL http://proceedings.mlr.press/v89/feydy19a.html.
  * Flamary et al. (2021) Rémi Flamary, Nicolas Courty, Alexandre Gramfort, Mokhtar Z. Alaya, Aurélie Boisbunon, Stanislas Chambon, Laetitia Chapel, Adrien Corenflos, Kilian Fatras, Nemo Fournier, Léo Gautheron, Nathalie T.H. Gayraud, Hicham Janati, Alain Rakotomamonjy, Ievgen Redko, Antoine Rolet, Antony Schutz, Vivien Seguy, Danica J. Sutherland, Romain Tavenard, Alexander Tong, and Titouan Vayer.  Pot: Python optimal transport.  _Journal of Machine Learning Research_ , 22(78):1–8, 2021.  URL http://jmlr.org/papers/v22/20-451.html.
  * Galichon (2016) Alfred Galichon.  _Optimal Transport Methods in Economics_.  Princeton University Press, August 2016.  ISBN 9781400883592, 9781400883592.  10.1515/9781400883592.
  * Gao and Chaudhari (2021) Yansong Gao and Pratik Chaudhari.  An information-geometric distance on the space of tasks.  In Marina Meila and Tong Zhang, editors, _Proceedings of the 38th International Conference on Machine Learning_ , volume 139 of _Proceedings of Machine Learning Research_ , pages 3553–3563. PMLR, 18–24 Jul 2021.  URL https://proceedings.mlr.press/v139/gao21a.html.
  * Genevay et al. (2018) Aude Genevay, Gabriel Peyré, and Marco Cuturi.  Learning generative models with sinkhorn divergences.  In _International Conference on Artificial Intelligence and Statistics_ , pages 1608–1617. PMLR, 2018.  URL https://proceedings.mlr.press/v84/genevay18a.html.
  * Howard et al. (2021) Frederick M. Howard, James Dolezal, Sara Kochanny, Jefree Schulte, Heather Chen, Lara Heij, Dezheng Huo, Rita Nanda, Olufunmilayo I. Olopade, Jakob N. Kather, Nicole Cipriani, Grossmank Robert L., and Alexander T. Pearson.  The impact of site-specific digital histology signatures on deep learning model accuracy and bias.  _Nature Communications_ , 2021.
  * Muhammad et al. (2021) Hassan Muhammad, Chensu Xie, Carlie S Sigel, Michael Doukas, Lindsay Alpert, Amber Lea Simpson, and Thomas J Fuchs.  EPIC-survival: End-to-end part inferred clustering for survival analysis, with prognostic stratification boosting.  In Mattias Heinrich, Qi Dou, Marleen de Bruijne, Jan Lellmann, Alexander Schläfer, and Floris Ernst, editors, _Proceedings of the Fourth Conference on Medical Imaging with Deep Learning_ , volume 143 of _Proceedings of Machine Learning Research_ , pages 520–531. PMLR, 07–09 Jul 2021.  URL https://proceedings.mlr.press/v143/muhammad21a.html.
  * Noorbakhsh et al. (2020) Javad Noorbakhsh, Saman Farahmand, Ali Foroughi pour, Sandeep Namburi, Dennis Caruana, Mohammad Rimm, David Soltanieh-ha, Kourosh Zarringhalam, and Jeffrey H. Chuang.  Deep learning-based cross-classifications reveal conserved spatial behaviors within tumor histological images.  _Nature communications_ , 11:6367, 2020.  https://doi.org/10.1038/s41467-020-20030-5.
  * Peyré and Cuturi (2019) Gabriel Peyré and Marco Cuturi.  Computational optimal transport.  _Foundations and Trends in Machine Learning_ , 11(5-6):355–607, 2019.  ISSN 1935-8237.  10.1561/2200000073.
  * Salimans et al. (2018) Tim Salimans, Han Zhang, Alec Radford, and Dimitris Metaxas.  Improving GANs using optimal transport.  In _International Conference on Learning Representations_ , 2018.  URL https://openreview.net/forum?id=rkQkBnJAb.
  * Santambrogio (2015) Filippo Santambrogio.  _Optimal Transport for Applied Mathematicians: Calculus of Variations, PDEs, and Modeling_.  Birkhäuser, Cham, 2015.  ISBN 9783319208275.  10.1007/978-3-319-20828-2.
  * Srinidhi et al. (2021) Chetan L. Srinidhi, Ozan Ciga, and Anne L. Martel.  Deep neural network models for computational histopathology: A survey.  _Medical Image Analysis_ , 67:101813, 2021.  ISSN 1361-8415.  https://doi.org/10.1016/j.media.2020.101813.  URL https://www.sciencedirect.com/science/article/pii/S1361841520301778.
  * Stacke et al. (2021) Karin Stacke, Gabriel Eilertsen, Jonas Unger, and Claes Lundström.  Measuring domain shift for deep learning in histopathology.  _IEEE Journal of Biomedical and Health Informatics_ , 25(2):325–336, 2021.  10.1109/JBHI.2020.3032060.
  * Villani (2003) Cédric Villani.  _Topics in Optimal Transportation_.  American Mathematical Soc., 2003.  ISBN 9780821833124.
  * Villani (2008) Cédric Villani.  _Optimal transport, Old and New_ , volume 338.  Springer Science & Business Media, 2008.  ISBN 9783540710493.
  * Wang et al. (2011) Wei Wang, John A. Ozolek, Dejan Slepčev, Ann B. Lee, Cheng Chen, and Gustavo K. Rohde.  An optimal transportation approach for nuclear structure-based pathology.  _IEEE Transactions on Medical Imaging_ , 30(3):621–631, 2011.  10.1109/TMI.2010.2089693.
  * Wulczyn et al. (2021) Ellery Wulczyn, David F. Steiner Steiner, Melissa Moran, Markus Plass, Robert Reihs, Fraser Tan, Isabelle Flament-Auvigne, Trissia Brown, Peter Regitnig, Po-Hsuan Cameron Chen, Narayan Hegde, Apaar Sadhwani, Robert MacDonald, Benny Ayalew, Greg S. Corrado, Lily H. Peng, Daniel Tse, Heimo Müller, Zhaoyang Xu, Yun Liu, Martin C. Stumpe, Kurt Zatloukal, and Craig H. Mermel.  Interpretable survival prediction for colorectal cancer using deep learning.  _npj Digital Medicine_ , 2021.
  * Yurochkin et al. (2019) Mikhail Yurochkin, Sebastian Claici, Edward Chien, Farzaneh Mirzazadeh, and Justin M Solomon.  Hierarchical optimal transport for document representation.  In H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett, editors, _Advances in Neural Information Processing Systems_ , volume 32. Curran Associates, Inc., 2019.  URL https://proceedings.neurips.cc/paper/2019/file/8b5040a8a5baf3e0e67386c2e3a9b903-Paper.pdf.


# Reinforcement Learning based on MPC/MHE for Unmodeled and Partially
Observable Dynamics

Hossein Nejatbakhsh Esfahani, Arash Bahari Kordabad, Sébastien Gros Authors
are with Department of Engineering Cybernetics, Norwegian University of
Science and Technology (NTNU), Trondheim, Norway. {hossein.n.esfahani,
arash.b.kordabad<EMAIL_ADDRESS>
###### Abstract

This paper proposes an observer-based framework for solving Partially
Observable Markov Decision Processes (POMDPs) when an accurate model is not
available. We first propose to use a Moving Horizon Estimation-Model
Predictive Control (MHE-MPC) scheme in order to provide a policy for the POMDP
problem, where the full state of the real process is not measured and
necessarily known. We propose to parameterize both MPC and MHE formulations,
where certain adjustable parameters are regarded for tuning the policy. In
this paper, for the sake of tackling the unmodeled and partially observable
dynamics, we leverage the Reinforcement Learning (RL) to tune the parameters
of MPC and MHE schemes jointly, with the closed-loop performance of the policy
as a goal rather than model fitting or the MHE performance. Illustrations show
that the proposed approach can effectively increase the performance of close-
loop control of systems formulated as POMDPs.

##  I Introduction

Reinforcement Learning (RL) is a powerful tool for solving Markov Decision
Processes (MDP) problems [1]. RL methods often use Deep Neural Network (DNN)
to approximate either the optimal policy underlying the MDP directly or the
action-value function from which the optimal policy can be indirectly
extracted.

Recent publications are discussing RL for POMDPs. A neural network-based
computation of belief states (posterior distributions over states) was
proposed to aggregate historical information needed to estimate a belief state
[2, 3]. An RL algorithm tailored to POMDPs was proposed in [4] that
incorporated spectral parameter estimation within an exploration-exploitation
strategy. A data-driven algorithm based on approximate Dynamic Programming
(ADP) was proposed in [5] to stabilize a plant with partially observable
dynamics. The authors used an action-dependent heuristic dynamic programming
(ADHDP) algorithm including two neural networks as an actor-critic (AC) method
to estimate both the unmeasured state and the performance index. The proposed
ADP-based approach in [6] is similar to classic RL algorithms but requires
only measurements of the input/output data and not of the full system state.

In [7], a neural network-based actor-critic structure was proposed to
approximate the control policies where a full system state is not accessible.
In [8], a fuzzy neural network was used to find the local optimal policy. A
recurrent neural network (RNN) was proposed in [9] to learn and infer the true
state observations from the noisy and correlated observations in a POMDP. Most
of the proposed RL-based control techniques in the above literature are based
on DNN-based approximators.

Model predictive control (MPC) is a popular and widely used practical approach
to optimal control. MPC is often selected for its capability to handle both
input and state constraints [10]. At each time instant, MPC calculates the
input and corresponding state sequence minimizing a cost function while
satisfying the constraints over a given prediction horizon. The first input is
applied, and the optimal solution is recalculated at the next time instant
based on the latest state of the system.

In many practical applications, some states of the plant are estimated using
an observer since they can not be directly measured, and the plant is possibly
not fully observable. The Moving Horizon Estimation (MHE) is a well-known
model-based observer in order to estimate the states of processes. In this
paper, we use this type of observer as a natural choice for the MPC scheme
[11].

Recently, the integration of machine learning in model predictive control has
been presented, with the aim of learning the model of the system, the cost
function or even the control law directly [12, 13]. These approaches are based
on DNN-based approximation. The direct combination of RL and MPC has been
investigated in [14, 15, 16]. It is shown that a single MPC scheme can capture
the optimal value function, action-value function, and policy of an MDP, even
if the MPC model is inaccurate, hence providing a valid and generic function
approximator for RL. The applications of this new MPC-based RL framework have
been recently presented in [17, 18].

However, these approaches assume that the state of the process is known and
can be fully measured. For many applications this assumption is not fulfilled.
To address this issue, this paper proposes to use a state observer such as an
MHE combined with the MPC scheme to build a policy based on the historic of
the available measurements rather than on the full state of the system. MHE
delivers state estimations by fitting the process model trajectory to past
measurements obtained on the real system. We adopt an MPC-based Q-learning
algorithm to tune the parameters included in the MHE-MPC scheme for the
closed-loop performance of the resulting policy.

This paper is organised as follows. In Section II, some background material is
given. Then the MPC and MHE schemes are detailed. The implementation of the
Q-learning algorithm for tuning both the MPC and MHE schemes together is
detailed in Section III. An illustrative example is proposed in Section IV.
Finally, conclusions and future work are given in Section V.

##  II Preliminaries and problem formulation

In the context of reinforcement learning a partially observable real plant is
described by a discrete POMDP having (possibly) stochastic state transitions
as follows:

$\displaystyle\boldsymbol{\mathrm{x}}_{k+1}=f^{\text{plant}}(\boldsymbol{\mathrm{x}}_{k},\boldsymbol{\mathrm{u}}_{k},\boldsymbol{\mathrm{\zeta}})$
(1a)
$\displaystyle\boldsymbol{\mathrm{y}}_{k}=h(\boldsymbol{\mathrm{x}}_{k},\boldsymbol{\mathrm{\eta}})$
(1b)

where the full state $\boldsymbol{\mathrm{x}}_{k}$ is not measurable or not
even known and $\boldsymbol{\mathrm{x}}_{k+1}$ is the next plant state vector
under stochastic transition with some random disturbances
$\boldsymbol{\mathrm{\zeta}}$. The model outputs $\boldsymbol{\mathrm{y}}$ are
measured on the real system and delivered by the output function $h$
associated with some random measurement noises $\eta$. We present next the MHE
and MPC schemes and how they can be used to create the action-value function
approximation required in Q-learning.

###  II-A Parameterized MHE Formulation

For a POMDP, the measurements available from the real process at a given time
instant do not constitute a Markov state. As a result, the full history of the
measurements becomes possibly relevant to the optimal policy. However,
building a policy based on the complete measurement history to solve the
POMDPs is not realistic. The RL community either considers a limited sequence
of past observations as a sufficient history or estimates a belief state using
a recurrent neural network. In this paper, we propose a more structured
solution to address this issue, by using MHE as a model-based approach to
build a state from the measurement history. The complete measurement history
is then transformed into a (possibly small) model state that is compatible
with the selected policy.

The MHE-based observer at the physical time $k$ can be stated as the following
Nonlinear Least-Squares problem:

$\displaystyle\left\\{\hat{\boldsymbol{\mathrm{x}}}_{k-N_{\text{MHE}},\ldots,k},\hat{\boldsymbol{\mathrm{u}}}_{k-N_{\text{MHE}},\ldots,k-1}\right\\}$
$\displaystyle=\operatorname*{arg\,min}_{\boldsymbol{\mathrm{x}},\boldsymbol{\mathrm{u}}}\left\|\boldsymbol{\mathrm{x}}_{k-N_{\text{MHE}}}-\tilde{\boldsymbol{\mathrm{x}}}_{k-N_{\text{MHE}}}\right\|_{A^{\theta}_{r}}^{2}$
$\displaystyle\quad+\sum_{i=k-N_{\text{MHE}}}^{k}\left\|\bar{\boldsymbol{\mathrm{y}}}_{i}-\boldsymbol{\mathrm{y}}(\boldsymbol{\mathrm{x}}_{i})\right\|_{Q^{\theta}_{E}}^{2}+\phi_{\theta}(\boldsymbol{\mathrm{x}}_{i})$
$\displaystyle\quad+\sum_{i=k-N_{\text{MHE}}}^{k-1}\left\|\boldsymbol{\mathrm{u}}_{i}-\bar{\boldsymbol{\mathrm{u}}}_{i}\right\|_{R^{\theta}_{E}}^{2}+\phi_{\theta}(\boldsymbol{\mathrm{u}}_{i})$
(2a)
$\displaystyle\quad\mathrm{s.t.}\quad\boldsymbol{\mathrm{x}}_{i+1}=f_{\theta}^{\text{MHE}}(\boldsymbol{\mathrm{x}}_{i},\boldsymbol{\mathrm{u}}_{i})$
(2b)

where $k$ is the current time instant, $i$ is the time instant along the
estimation horizon window.
$\bar{\boldsymbol{\mathrm{y}}}_{i},\bar{\boldsymbol{\mathrm{u}}}_{i}$ are the
measurements available at the physical time $k$ while their corresponding
values obtained from the MHE model are
$\boldsymbol{\mathrm{y}}(\boldsymbol{\mathrm{x}}_{i}),\boldsymbol{\mathrm{u}}_{i}$,
respectively. Let us to consider the mismatch between the model (observer) and
the real plant measurements is explainable by normal centered output noise,
then matrices $Q^{\theta}_{E}$ and $R^{\theta}_{E}$ are the inverse of the
covariance matrices associated to these noises on the plant output and control
input measurements, respectively. The first term in (2a) is an arrival cost
weighted with matrix $A^{\theta}_{r}$, which aims at approximating the
information prior to $k-N_{\text{MHE}}$, where
$\tilde{\boldsymbol{\mathrm{x}}}$ is the available estimation for the state at
time $k-N_{\text{MHE}}$. In practice, since the MHE fitting error is not only
coming from some normal centered output noise but also model error, more
intricate noise, and possibly unmodelled dynamics, it is very difficult to
decide what symmetric positive semi-definite weighting matrices
$Q^{\theta}_{E}$, $R^{\theta}_{E}$, $A^{\theta}_{r}$ ought to be used to
obtain the best closed-loop performance. To address this issue, we propose to
adjust them using the RL algorithm. Moreover, as the Least-Squares cost as a
choice of penalty in the MHE are not necessarily sufficient, we introduce a
cost modification $\phi_{\theta}$ tuned by RL. Note that we consider a
gradient form of the cost modification in this paper
$\phi_{\theta}(\boldsymbol{\mathrm{x}}_{i})=f_{1}^{\top}\boldsymbol{\mathrm{x}}_{i}$
and
$\phi_{\theta}(\boldsymbol{\mathrm{u}}_{i})=f_{2}^{\top}\boldsymbol{\mathrm{u}}_{i}$,
where $f_{1}$ and $f_{2}$ are labeled as RL parameters $\theta$.

###  II-B Parameterized MPC Formulation

In this work we will consider the MPC scheme as a value function approximator
that can be formulated as:

$\displaystyle V_{\theta}(\boldsymbol{\mathrm{x}}_{k})=$
$\displaystyle\min_{\boldsymbol{\mathrm{x}},\boldsymbol{\mathrm{u}},\boldsymbol{\mathrm{\sigma}}}\quad\gamma^{k+N_{\text{MPC}}}\left(V^{f}_{\theta}(\boldsymbol{\mathrm{x}}_{k+N_{\text{MPC}}})+\boldsymbol{\mathrm{w}}_{f}^{\top}\boldsymbol{\mathrm{\sigma}}_{k+N_{\text{MPC}}}\right)$
$\displaystyle\quad+\sum_{i=k}^{k+N_{\text{MPC}}-1}\gamma^{i}\left(l_{\theta}(\boldsymbol{\mathrm{x}}_{i},\boldsymbol{\mathrm{u}}_{i})+\boldsymbol{\mathrm{w}}^{\top}\boldsymbol{\mathrm{\sigma}}_{i}\right)$
(3a) $\displaystyle\mathrm{s.t.}$
$\displaystyle\quad\boldsymbol{\mathrm{x}}_{i+1}=f_{\theta}^{\text{MPC}}(\boldsymbol{\mathrm{x}}_{i},\boldsymbol{\mathrm{u}}_{i}),$
(3b)
$\displaystyle\quad\boldsymbol{\mathrm{x}}_{k}=\hat{\boldsymbol{\mathrm{x}}}_{k},$
(3c)
$\displaystyle\quad\boldsymbol{\mathrm{g}}(\boldsymbol{\mathrm{u}}_{i})\leq
0,$ (3d)
$\displaystyle\quad\boldsymbol{\mathrm{h}}_{\theta}(\boldsymbol{\mathrm{x}}_{i},\boldsymbol{\mathrm{u}}_{i})\leq\boldsymbol{\mathrm{\sigma}}_{i},\quad\boldsymbol{\mathrm{h}}_{\theta}^{f}(\boldsymbol{\mathrm{x}}_{k+N_{\text{MPC}}})\leq\boldsymbol{\mathrm{\sigma}}_{k+N_{\text{MPC}}}$
(3e)
$\displaystyle\quad\boldsymbol{\mathrm{\sigma}}_{k,\ldots,k+N_{\text{MPC}}-1}\geq
0$ (3f)

where $l_{\theta}$ is the stage cost, $V^{f}_{\theta}$ the terminal cost,
$f_{\theta}^{\text{MPC}}$ the MPC model (possibly but not necessarily
different from the MHE model), $\boldsymbol{\mathrm{h}}_{\theta}$ the mixed
constraints, $\boldsymbol{\mathrm{g}}$ the pure input constraints,
$\boldsymbol{\mathrm{h}}_{\theta}^{f}$ the terminal constraints. The MPC
initial conditions in (3c) are delivered by MHE scheme at the current time
instant $k$. In many real processes, there are uncertainties and disturbances
that may cause an MPC scheme to become infeasible. Therefore, an $\ell_{1}$
relaxation of the mixed constraints (3f) is introduced. An exact penalty is
imposed on the corresponding slack variables
$\boldsymbol{\mathrm{\sigma}}_{k}$ with large enough weights
$\boldsymbol{\mathrm{w}},\boldsymbol{\mathrm{w}}_{f}$ such that the
trajectories predicted by the MPC scheme will respect the constraints. All
elements in the above MPC scheme are parameterized by $\theta$, which will be
adjusted by RL, as detailed in [14].

Let us consider the policy at the physical current time $k$ as:

$\displaystyle\pi_{\theta}(\boldsymbol{\mathrm{x}}_{k})=\boldsymbol{\mathrm{u}}_{k}^{*}$
(4)

where, $\boldsymbol{\mathrm{u}}_{k}^{\star}$ is the first element of the input
sequence
$\boldsymbol{\mathrm{u}}_{k}^{\star},\cdots,\boldsymbol{\mathrm{u}}_{k+N_{\text{MPC}}-1}^{\star}$
solution of (3). We next consider this optimal policy delivered by the MPC
scheme as an action $\boldsymbol{\mathrm{a}}_{k}$ in the context of
reinforcement learning where, it is selected according to the above policy
with the possible addition of exploratory moves [1]. Then, an action-value
function approximation $Q_{\theta}$ can be formulated as:

$\displaystyle
Q_{\theta}(\boldsymbol{\mathrm{x}}_{k},\boldsymbol{\mathrm{u}}_{k})=\min_{\boldsymbol{\mathrm{x}},\boldsymbol{\mathrm{u}},\boldsymbol{\mathrm{\sigma}}}$
$\displaystyle\quad\eqref{eq:v1}$ (5a) $\displaystyle\mathrm{s.t.}$
$\displaystyle\quad\eqref{eq:v2}-\eqref{eq:violation}$ (5b)
$\displaystyle\quad\boldsymbol{\mathrm{u}}_{k}=\boldsymbol{\mathrm{a}}_{k}$
(5c)

Note that the proposed approximations (3)-(5) satisfies the fundamental
equalities underlying the Bellman equations [19]:

$\displaystyle\pi_{\theta}(\boldsymbol{\mathrm{x}}_{k})=\arg\min_{\boldsymbol{\mathrm{u}}}Q_{\theta}(\boldsymbol{\mathrm{x}}_{k},\boldsymbol{\mathrm{u}}_{k}),\quad
V_{\theta}(\boldsymbol{\mathrm{x}}_{k})=\min_{\boldsymbol{\mathrm{u}}}Q_{\theta}(\boldsymbol{\mathrm{x}}_{k},\boldsymbol{\mathrm{u}}_{k})$
(6)

##  III MPC/MHE-based RL

In this section, we present the algorithmic details needed to implement a
classic Q-learning algorithm on the combination of MPC/MHE schemes.

###  III-A Q-Learning for MPC/MHE

A classical off-policy Q-Learning algorithm is based on the temporal-
difference learning procedure [1] in which the updating rule for the RL
parameters can be expressed as follows:

$\displaystyle\delta_{k}=L(\boldsymbol{\mathrm{x}}_{k},\boldsymbol{\mathrm{u}}_{k})$
$\displaystyle+\gamma
V_{\theta}(\boldsymbol{\mathrm{x}}_{k+1})-Q_{\theta}(\boldsymbol{\mathrm{x}}_{k},\boldsymbol{\mathrm{u}}_{k}),$
(7a) $\displaystyle\theta\leftarrow\theta$
$\displaystyle+\alpha\delta_{k}\nabla_{\theta}Q_{\theta}(\boldsymbol{\mathrm{x}}_{k},\boldsymbol{\mathrm{u}}_{k})$
(7b)

where scalar $\alpha>0$ is a step size, $0<\gamma\leq 1$ is a discount factor
and $\delta_{k}$ is the Temporal Difference (TD) error at the physical time
$k$. In the above TD learning algorithm, a baseline stage cost
$L(\boldsymbol{\mathrm{x}}_{k},\boldsymbol{\mathrm{u}}_{k})$ (reward in the
context of RL) is defined as a function of state-action pair in order to
provide an evaluation signal. Indeed, the baseline cost affects the agent
behavior and control policy via RL parameter updating, where the TD error is
appeared.

The gradient of function $Q_{\theta}$ needed in (7) requires one to compute
the sensitivities of the optimal value of Nonlinear Programming (NLP) (5).
This sensitivity ought to be computed with care since the RL parameters
$\theta$ impact $Q_{\theta}$ both directly via the MPC scheme and indirectly
via the MHE scheme, by modifying the state estimation
$\hat{\boldsymbol{\mathrm{x}}}_{k}$ at the physical current time instant $k$
that enters as an initial condition
$\boldsymbol{\mathrm{x}}_{k}=\hat{\boldsymbol{\mathrm{x}}}_{k}$ in the MPC
scheme. The gradient $\nabla_{\theta}Q_{\theta}$ associated to the proposed
MPC/MHE scheme is given by the following total derivative:

$\displaystyle\frac{dQ_{\theta}}{d\theta}=\frac{\partial
Q_{\theta}}{\partial\theta}+\frac{\partial
Q_{\theta}}{\partial\boldsymbol{\mathrm{x}}_{k}}\frac{\partial\hat{\boldsymbol{\mathrm{x}}}_{k}}{\partial\theta},\quad\boldsymbol{\mathrm{x}}_{k}=\hat{\boldsymbol{\mathrm{x}}}_{k}$
(8)

We detail next how to compute the above sensitivities.

###  III-B Sensitivities of the MPC/MHE scheme

Let us define the Lagrange functions
$\mathcal{L}_{\theta},\hat{\mathcal{L}_{\theta}}$ associated to the MPC and
MHE problems (2), (5) as follows:

$\displaystyle\mathcal{L}_{\theta}=\Phi_{\theta}+\boldsymbol{\mathrm{\lambda}}^{\top}G_{\theta}+\boldsymbol{\mathrm{\mu}}^{\top}H_{\theta}$
(9)
$\displaystyle\hat{\mathcal{L}}_{\theta}=\hat{\Phi}_{\theta}+\hat{\boldsymbol{\mathrm{\lambda}}}^{\top}\hat{G}_{\theta}$
(10)

where $H_{\theta}$ gathers the inequality constraints of (5) and
$\Phi_{\theta},\hat{\Phi}_{\theta}$ are the costs of the MPC and MHE
optimization problems, respectively. Variables
$\boldsymbol{\mathrm{\lambda}},\hat{\boldsymbol{\mathrm{\lambda}}}$ are the
Lagrange multipliers associated to the equality constraints
$G_{\theta},\hat{G}_{\theta}$ of the MPC and MHE, respectively. Variables
$\boldsymbol{\mathrm{\mu}}$ are the Lagrange multipliers associated to the
inequality constraints of the MPC scheme. Let us label the primal variables as
$\boldsymbol{\mathrm{p}}=\left\\{\boldsymbol{\mathrm{X}},\boldsymbol{\mathrm{U}}\right\\}$
and
$\hat{\boldsymbol{\mathrm{p}}}=\left\\{\hat{\boldsymbol{\mathrm{X}}},\hat{\boldsymbol{\mathrm{U}}}\right\\}$
for the MPC and MHE, respectively. The primal-dual variables of the MPC and
MHE schemes will be labeled as
$\boldsymbol{\mathrm{z}}=\left\\{\boldsymbol{\mathrm{p}},\boldsymbol{\mathrm{\lambda}},\boldsymbol{\mathrm{\mu}}\right\\}$
and
$\hat{\boldsymbol{\mathrm{z}}}=\left\\{\hat{\boldsymbol{\mathrm{p}}},\hat{\boldsymbol{\mathrm{\lambda}}}\right\\}$,
respectively.

The sensitivities of the MPC scheme (5) required in (8) can be obtained by the
sensitivity analysis detailed in [20] as follows:

$\displaystyle\frac{\partial
Q_{\theta}}{\partial\theta}=\frac{\partial\mathcal{L}_{\theta}(\boldsymbol{\mathrm{x}}_{k},\boldsymbol{\mathrm{z}}^{\star})}{\partial\theta},\quad\frac{\partial
Q_{\theta}}{\partial\boldsymbol{\mathrm{x}}_{k}}=\frac{\partial\mathcal{L}_{\theta}(\boldsymbol{\mathrm{x}}_{k},\boldsymbol{\mathrm{z}}^{\star})}{\partial\boldsymbol{\mathrm{x}}_{k}}$
(11)

where $\boldsymbol{\mathrm{z}}^{\star}$ is the primal-dual solution vector of
(5).

The sensitivity
$\frac{\partial\hat{\boldsymbol{\mathrm{x}}}_{k}}{\partial\theta}$ associated
to the MHE scheme can be obtained via using the Implicit Function Theorem
(IFT) on the Karush Kuhn Tucker (KKT) conditions underlying the parametric
NLP. Assuming that Linear Independence Constraint Qualification (LICQ) and
Second Order Sufficient Condition (SOSC) hold [21] at
$\hat{\boldsymbol{\mathrm{z}}}^{\star}$, then, the following holds:

$\displaystyle\frac{\partial\hat{\boldsymbol{\mathrm{z}}}^{\star}}{\partial\theta}=-\frac{\partial
R_{\theta}}{\partial\hat{\boldsymbol{\mathrm{z}}}}^{-1}\frac{\partial
R_{\theta}}{\partial\theta}$ (12)

where

$\displaystyle
R_{\theta}=\begin{bmatrix}\nabla_{\hat{\boldsymbol{\mathrm{p}}}}\hat{\mathcal{L}}_{\theta}\\\
\hat{G}_{\theta}\end{bmatrix}$ (13)

are the KKT conditions associated to the MHE scheme (2). As
$\hat{\boldsymbol{\mathrm{x}}}_{k}$ is part of
$\hat{\boldsymbol{\mathrm{z}}}^{\star}$, the sensitivity of the MHE solution
$\frac{\partial\hat{\boldsymbol{\mathrm{x}}}_{k}}{\partial\theta}$ required in
(8) can be extracted from matrix
$\frac{\partial\hat{\boldsymbol{\mathrm{z}}}^{\star}}{\partial\theta}$.

###  III-C Constrained RL steps

The adjustable weighting matrices in the proposed parameterization of both MPC
and MHE in (2), (3) and (5) are tuned using Q-learning. As a requirement, the
weighting matrices $Q^{\theta}_{E}$, $R^{\theta}_{E}$, $A^{\theta}_{r}$ must
be positive semidefinite. However, the RL steps delivered by Q-learning do not
necessarily respect this requirement, and we need to enforce it via
constraints on the RL steps throughout the learning process. To address this
requirement, we formulate a Semi-Definite Program (SDP) as a least squares
optimization problem:

$\displaystyle\min_{\Delta\theta}$
$\displaystyle\quad\frac{1}{2}\left\|\Delta\theta\right\|^{2}-\alpha\delta_{k}\nabla_{\theta}Q_{\theta}(\boldsymbol{\mathrm{x}}_{k},\boldsymbol{\mathrm{u}}_{k})^{\top}\Delta\theta$
(14a) $\displaystyle\mathrm{s.t.}$ $\displaystyle\quad
Q^{\theta}_{E}({\theta+\Delta\theta})\geq 0$ (14b) $\displaystyle\quad
R^{\theta}_{E}({\theta+\Delta\theta})\geq 0$ (14c) $\displaystyle\quad
A^{\theta}_{r}({\theta+\Delta\theta})\geq 0$ (14d)

where we assume that the weighting matrices
$Q^{\theta}_{E},R^{\theta}_{E},A^{\theta}_{r}$ are linear functions of
$\theta$. Then, these matrices are updated in each time instant due to
updating $\Delta\theta$, which is a solution of the above SDP scheme. The
proposed learning process is described in the Alg. 1.

Algorithm 1 (MPC+MHE)-Based RL

0:
$\alpha,tol>0,\theta=\theta_{0},\boldsymbol{\mathrm{x}}_{0},\boldsymbol{\mathrm{u}}_{0}$

while Iter do

1\. Measure output $\boldsymbol{\mathrm{y}}_{k}$ from (1b) at current time $k$

2\. Obtain
$\hat{\boldsymbol{\mathrm{x}}}_{k}$,$\frac{\partial\hat{\boldsymbol{\mathrm{x}}}_{k}}{\partial\theta}$
from (2) and (12)

3\. Obtain $\pi_{\theta}(\boldsymbol{\mathrm{x}}_{k})$,
$V_{\theta}(\boldsymbol{\mathrm{x}}_{k})$ from (3)

4\. Exploration:
$\boldsymbol{\mathrm{u}}_{k}=\pi_{\theta}(\boldsymbol{\mathrm{x}}_{k})+d$,
$d\sim\mathcal{N}\left(\mu,\sigma^{2}\right)$

5\. Obtain $Q_{\theta}$, $\frac{\partial Q_{\theta}}{\partial\theta}$ from
(5), (11)

6\. Assemble $\frac{dQ_{\theta}}{d\theta}$ from (8)

7\. Apply $\boldsymbol{\mathrm{u}}_{k}$ to the real plant: (1a)

8\. Evaluate baseline
$L(\boldsymbol{\mathrm{x}}_{k},\boldsymbol{\mathrm{u}}_{k})$

9\. RL update:

\- Obtain $\Delta\theta$ from (14)

\- $\theta\leftarrow\theta+\Delta\theta$

$k\leftarrow k+1$

end while

##  IV Numerical Example

In this section, we illustrate the performance of the proposed MPC/MHE-based
RL scheme, which is tested on a constrained two-mass-spring-damper system
shown in Figure 1, for which the MPC/MHE model ignores some of the dynamics.

Figure 1: Two-Mass-Spring-Damper

The control input acts on mass 1, and the position of mass 2 is measured. Let
us consider $m_{1}=0.8$ kg, $m_{2}=0.5$ kg, $k=25$ $\frac{N}{m}$, $b=3$
$\frac{Ns}{m}$ and define the plant dynamics as:

$\displaystyle\begin{bmatrix}\dot{x}_{1}\\\ \dot{x}_{2}\\\ \dot{x}_{3}\\\
\dot{x}_{4}\end{bmatrix}=\begin{bmatrix}\phantom{-}0&\phantom{-}0&\phantom{-}1&\phantom{-}0\\\
\phantom{-}0&\phantom{-}0&\phantom{-}0&\phantom{-}1\\\
-\frac{k}{m_{1}}&\phantom{-}\frac{k}{m_{1}}&-\frac{b}{m_{1}}&\phantom{-}\frac{b}{m_{1}}\\\
\phantom{-}\frac{k}{m_{2}}&-\frac{k}{m_{2}}&\phantom{-}\frac{b}{m_{2}}&-\frac{b}{m_{2}}\end{bmatrix}\begin{bmatrix}x_{1}\\\
x_{2}\\\ x_{3}\\\ x_{4}\end{bmatrix}+\begin{bmatrix}0\\\ 0\\\
\frac{1}{m_{1}}\\\ 0\end{bmatrix}u$

where $x_{1},x_{2}$ are positions of masses 1, 2, respectively, and
$x_{3},x_{4}$ are corresponding velocities. Variable $u$ is the control input
applied to the first mass. In this simulation, we propose to formulate the
MPC/MHE scheme based on a partially observable model. More precisely, the
adopted MPC scheme is presented based on a 2-states model, capturing only the
position and velocity of mass 1. The MHE scheme is based on the same model,
but is fed as measurements the position of the mass 2. Let
$f_{\theta}^{\text{MPC}}$ in (3b) be a partially observable and inaccurate MPC
model as follows:

$\displaystyle\begin{bmatrix}\dot{x}_{1}\\\
\dot{x}_{3}\end{bmatrix}=\left(\begin{bmatrix}0&1\\\
0&0\end{bmatrix}+A^{\text{bias}}\right)\begin{bmatrix}x_{1}\\\
x_{3}\end{bmatrix}+\left(\begin{bmatrix}0\\\
\frac{1}{m_{1}+m_{2}}\end{bmatrix}+B^{\text{bias}}\right)u$ (15)

where $A^{\text{bias}},B^{\text{bias}}$ are adjusted as model bias RL
parameters $\theta$ via Q-learning in order to tackle the inaccurate above MPC
model. In this example, the following optimization problem as a parameterized
MPC scheme is solved at each time instant $k$:

$\displaystyle\min_{\boldsymbol{\mathrm{x}},\boldsymbol{\mathrm{u}},\boldsymbol{\mathrm{\sigma}}}\quad$
$\displaystyle\theta_{c}+\frac{\gamma^{N^{f}}}{2}\left(\boldsymbol{\mathrm{x}}^{\top}_{N^{f}}M^{f}_{\theta}\boldsymbol{\mathrm{x}}_{N^{f}}+\boldsymbol{\mathrm{w}}^{\top}_{f}\boldsymbol{\mathrm{\sigma}}_{N^{f}}\right)$
$\displaystyle+\sum_{i=k}^{k+N_{\text{MPC}}-1}\frac{\gamma^{i}}{2}\left((c-\theta_{r})^{\top}M_{\theta}(c-\theta_{r})+\boldsymbol{\mathrm{w}}^{\top}\boldsymbol{\mathrm{\sigma}}_{i}\right)$
(16a) $\displaystyle\mathrm{s.t.}$
$\displaystyle\quad\boldsymbol{\mathrm{x}}_{i+1}=f_{\theta}^{\text{MPC}}(\boldsymbol{\mathrm{x}}_{i},\boldsymbol{\mathrm{u}}_{i}),$
(16b)
$\displaystyle\quad\boldsymbol{\mathrm{x}}_{k}=\hat{\boldsymbol{\mathrm{x}}}_{k},$
(16c) $\displaystyle\quad-1\leq\boldsymbol{\mathrm{u}}_{i}\leq 1,$ (16d)
$\displaystyle\quad\begin{bmatrix}0\\\
-10\end{bmatrix}+\underline{\theta}-\boldsymbol{\mathrm{\sigma}}_{i}\leq\boldsymbol{\mathrm{x}}_{i}\leq\begin{bmatrix}10\\\
10\end{bmatrix}+\bar{\theta}+\boldsymbol{\mathrm{\sigma}}_{i}$ (16e)
$\displaystyle\quad\boldsymbol{\mathrm{\sigma}}_{k,\ldots,k+N_{\text{MPC}}-1}\geq
0$ (16f)

where $N^{f}=k+N_{\text{MPC}}$,
$c=\left[\boldsymbol{\mathrm{x}}_{k},\boldsymbol{\mathrm{u}}_{k}\right]^{\top}$
and the MPC parameters subject to the RL scheme are:

$\displaystyle\theta=\left(\theta_{c},M^{f}_{\theta},\theta_{r},M_{\theta},A^{\text{bias}},B^{\text{bias}},\underline{\theta},\bar{\theta}\right)$
(17)

The positive semidefinite weighting matrices ($M^{f}_{\theta}$,
$M_{\theta}$,$Q^{\theta}_{E}$, $R^{\theta}_{E}$, $A^{\theta}_{r}$) in both MPC
and MHE schemes are adjusted using the constrained RL steps in (14). One can
choose a baseline stage cost used in the RL scheme (7) as:

$\displaystyle
L(\boldsymbol{\mathrm{x}}_{k},\boldsymbol{\mathrm{u}}_{k})=l(\boldsymbol{\mathrm{x}}_{k},\boldsymbol{\mathrm{u}}_{k})+\boldsymbol{\mathrm{w}}^{\top}\max(0,\boldsymbol{\mathrm{h}}(\boldsymbol{\mathrm{x}}_{k}))$
(18)

where $l(\boldsymbol{\mathrm{x}}_{k},\boldsymbol{\mathrm{u}}_{k})$ is adopted
as a quadratic function of the state and action deviations from their desired
values. The second term in the above baseline is considered to penalize the
constraint violations, where $\boldsymbol{\mathrm{h}}\geq 0$ is pure
inequality vector of constraints on the states and
$\boldsymbol{\mathrm{w}}=[10,10]$. Note that different step sizes $\alpha$
were used for the different parameters based on the problem scaling. The
desired values for the MPC model states (position and velocity of the first
body) $x_{1},x_{3}$ are chosen at $[0,0]^{T}$ , respectively. We apply a
process noise $\zeta\sim\mathcal{N}(\mu=0,\,(\sigma=0.02)^{2})$ on the
velocity of the second body and a measurement noise
$\eta\sim\mathcal{N}(\mu=0,\,(\sigma=0.05)^{2})$ on the position of second
body. In this simulation we choose $\gamma=0.9$, and
$N_{\text{MPC}}=N_{\text{MHE}}=8$.

As this simulation is considered as a POMDP and uncertain scenario and there
are both process and measurement Gaussian noises, the violations are observed
on the states in Figure 5. We demonstrate that the proposed MPC/MHE-based RL
can attenuate these violations and increase the control performance even if
the controller/observer models are unmodeled and partially observable.

In this example, we consider three different scenarios.

1) Without learning (MPC+MHE): In the first scenario, there are large
violations of the position and velocity constraints affecting the closed-loop
performance (large cost $J$) shown in Figure 5 and Figure 6.

2) MPC-based RL learning (MPC-RL+MHE): In the second scenario, the learning is
only performed on the MPC scheme. This MPC-based RL reduces the violations and
increases the closed-loop performance shown in Figure 5 by reducing the
discounted sum of the RL stage cost $J$ over a receding horizon. There is also
a decrease of the baseline cost shown in Figure 7 after starting the MPC
learning while the MHE learning is not still activated.

3) MPC/MHE-based RL learning (MPC-RL+MHE-RL): Finally in the third scenario,
the performance is improved after allowing the MHE to be adjusted using the
Q-learning algorithm and there is a solid decrease in the TD-error and
baseline cost and an increasing closed-loop performance (decrease of $J$). The
evolution of MHE parameters are illustrated in Figure 8 and Figure (9). The
evolution of MPC parameters are depicted in Figures 2, 3, and 4.

Figure 2: MPC adjustment: reference signals Figure 3: MPC adjustment:
Constraint and stage cost Figure 4: MPC adjustment: model bias Figure 5:
Positions of masses and closed-loop performance. The brown lines are shown as
lower bound (0 $mm$) and upper bound (10 $mm$) constraints on the positions.
Position references are (0 $mm$). Figure 6: Control input and velocities of
masses. The brown lines are shown as lower bound (-10 $mm/s$) and upper bound
(10 $mm/s$) constraints on the velocities. Figure 7: RL performance: baseline
cost and TD error Figure 8: MHE adjustment: Arrival matrix and penalizing
weights Figure 9: Evolution of arrival cost and gradients in MHE

##  V Conclusion

This paper proposed the combination of MPC-based Reinforcement Learning with
an MHE scheme to tackle POMDPs. The introduction of an MHE scheme allows to
deploy MPC-based Reinforcement Learning without a full state measurement, and
without necessarily holding a correct representation of the system state in
the MPC and MHE models. Furthermore, we propose to tune the MHE and MPC
schemes jointly, focusing directly on the closed-loop performance, as opposed
to using indirect criteria such as decreasing the MHE output error. We detail
the application of Q-learning to this approach, and test it in a simulated
spring mass example operating under constraints, where only a part of the real
system dynamics are modelled in the MPC and MHE schemes. We show that the
method manages to tune the MHE and MPC scheme to reduce the constraints
violations and improve the closed-loop performance. Future work will propose
an stability and feasibility analysis on the proposed MPC/MHE-based RL scheme.

## References

  * [1] R. S. Sutton and A. G. Barto, _Reinforcement learning: An introduction_. MIT press, 2018.
  * [2] Z. D. Guo, M. G. Azar, B. Piot, B. A. Pires, T. Pohlen, and R. Munos, “Neural predictive belief representations,” _CoRR_ , vol. abs/1811.06407, 2018.
  * [3] T. Gangwani, J. Lehman, Q. Liu, and J. Peng, “Learning belief representations for imitation learning in pomdps,” in _UAI_ , 2019.
  * [4] K. Azizzadenesheli, “Reinforcement learning in structured and partially observable environments,” Ph.D. dissertation, University of California, Irvine, 2019.
  * [5] X. Zhong, Z. Ni, Y. Tang, and H. He, ““data-driven partially observable dynamic processes using adaptive dynamic programming,” in _in Proc. IEEE Symposium on Adaptive Dynamic Programming and Reinforcement Learning_ , 2014, pp. 1–8.
  * [6] F. Lewis and K. Vamvoudakis, “‘reinforcement learning for partially observable dynamic processes: adaptive dynamic programming using measured output data,” _IEEE Trans. Syst. Man, Cybern. B, Cybern.,_ , vol. 41, no. 1, pp. 14–25, 2011.
  * [7] H. Jiang and H. He, “Data-driven distributed output consensus control for partially observable multiagent systems,” _IEEE Trans. Cybern._ , vol. 49, no. 3, pp. 848–858, 2019.
  * [8] A. Bahari Kordabad and M. Boroushaki, “Emotional learning based intelligent controller for mimo peripheral milling process,” _Journal of Applied and Computational Mechanics_ , vol. 6, no. 3, pp. 480–492, 2020.
  * [9] Y. Wang, K. Velswamy, and B. Huang, “A novel approach to feedback control with deep reinforcement learning,” _IFAC-PapersOnLine_ , vol. 51, no. 18, pp. 31 – 36, 2018, 10th IFAC Symposium on Advanced Control of Chemical Processes ADCHEM 2018.
  * [10] J. B. Rawlings, D. Q. Mayne, and M. Diehl, _Model predictive control: theory, computation, and design_. Nob Hill Publishing Madison, WI, 2017, vol. 2.
  * [11] P. Kuhl, M. Diehl, T. Kraus, J. P. Schloder, and H. G. Bock, “A real-time algorithm for moving horizon state and parameter estimation,” _Computers and Chemical Engineering_ , vol. 35, no. 1, pp. 71 – 83, 2011\.
  * [12] B. Karg and S. Lucia, “Efficient representation and approximation of model predictive control laws via deep learning,” _arXiv preprint arXiv:1806.10644_ , 2018.
  * [13] B. Karg and S. Lucia, “Learning-based approximation of robust nonlinear predictive control with state estimation applied to a towing kite,” in _2019 18th European Control Conference (ECC)_ , 2019, pp. 16–22.
  * [14] S. Gros and M. Zanon, “Data-driven economic nmpc using reinforcement learning,” _IEEE Transactions on Automatic Control_ , vol. 65, no. 2, pp. 636–648, 2019.
  * [15] M. Zanon, S. Gros, and A. Bemporad, “Practical reinforcement learning of stabilizing economic mpc,” in _2019 18th European Control Conference (ECC)_. IEEE, 2019, pp. 2258–2263.
  * [16] S. Gros and M. Zanon, “Reinforcement learning for mixed-integer problems based on mpc,” _arXiv preprint arXiv:2004.01430_ , 2020.
  * [17] A. Bahari Kordabad, W. Cai, and S. Gros, “MPC-based reinforcement learning for economic problems with application to battery storage,” in _2021 20th European Control Conference (ECC) (Accepted)_. IEEE, 2020.
  * [18] A. Bahari Kordabad, H. Nejatbakhsh Esfahani, A. M. Lekkas, and S. Gros, “Reinforcement learning based on scenario-tree MPC for ASVs (accepted),” in _2021 American Control Conference (ACC)_. IEEE, 2021.
  * [19] D. Bertsekas, “Dynamic programming and optimal control,” _Athena Scientific, 3rd edition_ , 2005.
  * [20] C. Buskens and H. Maurer., _Online Optimization of Large Scale Systems, chapter Sensitivity Analysis and Real-Time Optimization of Parametric Nonlinear Programming Problems, pages 3–16_. Berlin Heidelberg, Berlin, Heidelberg, 2001.
  * [21] J. Nocedal and S. Wright, _Numerical Optimization_. Springer, 2 edition, 2006.