End-to-End Reinforcement Learning of Koopman Models for Economic Nonlinear
Model Predictive Control

Daniel Mayfranka,d, Alexander Mitsosc,a,b, Manuel Dahmena,∗

  * a

Forschungszentrum Jülich GmbH, Institute of Climate and Energy Systems, Energy
Systems Engineering (ICE-1), Jülich 52425, Germany

  * b

RWTH Aachen University, Process Systems Engineering (AVT.SVT), Aachen 52074,
Germany

  * c

JARA-ENERGY, Jülich 52425, Germany

  * d

RWTH Aachen University, Aachen 52062, Germany

(Economic) nonlinear model predictive control ((e)NMPC) requires dynamic
models that are sufficiently accurate and computationally tractable. Data-
driven surrogate models for mechanistic models can reduce the computational
burden of (e)NMPC; however, such models are typically trained by system
identification for maximum prediction accuracy on simulation samples and
perform suboptimally in (e)NMPC. We present a method for end-to-end
reinforcement learning of Koopman surrogate models for optimal performance as
part of (e)NMPC. We apply our method to two applications derived from an
established nonlinear continuous stirred-tank reactor model. The controller
performance is compared to that of (e)NMPCs utilizing models trained using
system identification, and model-free neural network controllers trained using
reinforcement learning. We show that the end-to-end trained models outperform
those trained using system identification in (e)NMPC, and that, in contrast to
the neural network controllers, the (e)NMPC controllers can react to changes
in the control setting without retraining.

Keywords: Economic model predictive control; Koopman; Reinforcement learning;
End-to-end learning .

End-to-End Reinforcement Learning of Koopman Models for Economic Nonlinear
Model Predictive Control

Daniel Mayfranka,d, Alexander Mitsosc,a,b, Manuel Dahmena,∗

  * a

Forschungszentrum Jülich GmbH, Institute of Climate and Energy Systems, Energy
Systems Engineering (ICE-1), Jülich 52425, Germany

  * b

RWTH Aachen University, Process Systems Engineering (AVT.SVT), Aachen 52074,
Germany

  * c

JARA-ENERGY, Jülich 52425, Germany

  * d

RWTH Aachen University, Aachen 52062, Germany

(Economic) nonlinear model predictive control ((e)NMPC) requires dynamic
models that are sufficiently accurate and computationally tractable. Data-
driven surrogate models for mechanistic models can reduce the computational
burden of (e)NMPC; however, such models are typically trained by system
identification for maximum prediction accuracy on simulation samples and
perform suboptimally in (e)NMPC. We present a method for end-to-end
reinforcement learning of Koopman surrogate models for optimal performance as
part of (e)NMPC. We apply our method to two applications derived from an
established nonlinear continuous stirred-tank reactor model. The controller
performance is compared to that of (e)NMPCs utilizing models trained using
system identification, and model-free neural network controllers trained using
reinforcement learning. We show that the end-to-end trained models outperform
those trained using system identification in (e)NMPC, and that, in contrast to
the neural network controllers, the (e)NMPC controllers can react to changes
in the control setting without retraining.

Keywords: Economic model predictive control; Koopman; Reinforcement learning;
End-to-end learning

##  1 Introduction

Data-driven surrogate models present a promising avenue for rendering economic
nonlinear model predictive control (eNMPC) tractable by reducting the
computational demand associated with solving the underlying optimal control
problem (McBride and Sundmacher, (2019)). System identification (SI)
represents the predominant methodology employed for training data-driven
dynamic surrogate models; however, SI narrowly focuses on optimizing average
prediction accuracy on a collection of simulation samples. In contrast, recent
studies have demonstrated the superiority of dynamic surrogate models, trained
explicitly for optimal performance in control applications via reinforcement
learning (RL) techniques, over SI-trained models (Amos et al., (2018); Yin et
al., (2022); Iwata and Kawahara, (2022); Chen et al., (2019); Gros and Zanon,
(2019)). Nonetheless, these findings were limited to the learning of linear
models (Chen et al., (2019)), applications devoid of state constraints (Amos
et al., (2018); Yin et al., (2022); Iwata and Kawahara, (2022); Chen et al.,
(2019)), or primarily concentrated on the adaptation of bounds and cost
function rather than the dynamic MPC model itself (Gros and Zanon, (2019)).

RL research primarily focuses on learning model-free control policies that
eschew the utilization of system state predictions to dictate control actions.
In contrast, learning a dynamic model and employing its predictions to derive
a control law, as exemplified in eNMPC, has crucial advantages: (i) retraining
is rendered unnecessary should constraints or objective functions be altered,
provided that the system dynamics persist unaltered; (ii) learning the system
dynamics may be more sample-efficient than learning a model-free control
policy (Amos et al., (2018); Chen et al., (2019)); (iii) MPC has a well-
established theory regarding performance and stability guarantees, especially
for linear models and recent publications aim to extend this theory to include
(learned) eNMPC (Gros and Zanon, (2019); Angeli et al., (2011)).

Using RL to train dynamic surrogate models promises to combine the
aforementioned advantages of model-based policies with the typical advantage
of end-to-end learning over SI, i.e., the task-optimal performance of the
resulting model. To train dynamic surrogate models using RL algorithms,
(eN)MPC is viewed as a policy whose control outputs can be differentiated with
respect to the parameters of the surrogate model. This requires
differentiating a solution point of a parametric optimization problem with
respect to the parameters (Fiacco and Kyparisis, (1985); Fiacco and Ishizuka,
(1990); Ralph and Dempe, (1995)). Lewis and Vrabie, (2009) and Lewis et al.,
(2012) show that RL can be used to learn goal-directed dynamic models by
connecting RL to the linear quadratic regulator (Brunton and Kutz, (2022)), a
model predictive feedback controller for control settings with linear dynamics
and quadratic cost functions. Amos et al., (2018) and Gros and Zanon, (2019)
study how to learn MPC controllers end-to-end by viewing these as
differentiable policies. Amos et al., (2018) use imitation learning (Hussein
et al., (2017)), a supervised learning approach aiming to train a control
policy through expert demonstrations, instead of RL to learn an NMPC but do
not consider control settings with inequality constraints on state variables.
Moreover, they learn a parameterized version of the true underlying dynamic
equations, which requires an a priori known structure of the system dynamics.
Gros and Zanon, (2019) use RL in the form of deep Q-learning, a popular RL
algorithm for environments with discrete action spaces and the first RL
algorithm to claim superhuman performance on an Atari game (Mnih et al.,
(2013)), to learn linear MPC and eNMPC controllers. However, in the nonlinear
case, they focus on learning advantageous parameterizations of the constraints
and cost function and only parameterize the dynamic model with a learnable
constant added to the transition function. Chen et al., (2019) use both
imitation learning and RL together with the differentiable MPC solver by Amos
et al., (2018) to learn task-optimal linear dynamic models for HVAC control
without constraints on state variables.

Given assumptions we describe in Subsection 2.3, a solution point $\bm{a}^{*}$
of a parametric optimization problem is differentiable with respect to the
parameters $\bm{\theta}$ locally at a given $\bm{\theta}^{*}$. One of the
assumptions is that the second-order sufficient conditions (SOSC) for a local
minimum hold at $\bm{a}^{*}$, which, for convex parametric programs with
strictly convex objective functions, is easily verifiable through the Karush-
Kuhn-Tucker (KKT) conditions (Karush, (1939); Boyd et al., (2004)). Therefore,
for end-to-end learning of dynamic models for MPC, it is desirable to learn
dynamic models that lead to convex optimal control problems (OCPs). A simple
condition for convex OCPs are formulations with linear dynamics and convex
constraints. Koopman theory (Koopman, (1931)) aims to find globally-valid
linear approximations of nonlinear process dynamics through a nonlinear
coordinate transformation into a higher-dimensional linear state space.
Koopman theory and its extension to actuated systems for control applications
(Korda and Mezić, (2018), therefore, provides a powerful framework to obtain
surrogate models that can accurately represent continuous nonlinear dynamical
systems while still resulting in convex OCPs, especially when combined with
modern deep learning algorithms (Lusch et al., (2018)). Iwata and Kawahara,
(2022) and Yin et al., (2022) both use RL for end-to-end learning of task-
optimal Koopman models to obtain gain matrices that act on the Koopman
invariant subspaces, minimizing the squared distance of the system state from
a desired target state. However, similarly to linear-quadratic regulator
control (Kalman et al., (1960)), this approach cannot easily be extended to
(e)NMPC applications or applications with hard bounds on state variables.

Koopman-based MPCs have seen application in various domains due to their
favorable balance between high representational capacity and comparatively low
computational demand in solving the resulting OCPs. These applications range
from ground vehicle (Cibulka et al., (2020)) to air vehicle control (Folkestad
and Burdick, (2021)), flow control (Arbabi et al., (2018)), and soft robotics
(Bruder et al., (2019)). Some contributions have also studied the application
of Koopman-based MPCs to control of chemical processes. Narasingam and Kwon,
(2019) integrate Koopman models with Lyapunov-based MPC, thus guaranteeing
controller feasibility and closed-loop stability. Narasingam and Kwon, (2020)
integrate a priori system knowledge in the choice of the Koopman basis
functions and achieve successful data-driven control of a hydraulic fracturing
process. Using a learned disturbance estimator, Son et al., 2022b  develop
offset-free Koopman Lyapunov-based MPC to compensate for the unavoidable
plant-model-mismatch resulting from a finite-dimensional approximation of the
Koopman operator. They showcase their approach on the control of a batch pulp
digester (Son et al., (2021)). For hybrid systems with differing local
dynamics, Son et al., 2022a  develop a hybrid Koopman MPC approach by
clustering the training data and by training one local Koopman model for each
cluster. Furthermore, Albalawi and Hameed, (2023) demonstrate the efficacy of
Koopman-based economic MPC on a continuous stirred-tank reactor (CSTR) case
study.

We present a method for end-to-end RL of Koopman models for optimal
performance in (e)NMPC applications with hard constraints on states. Using the
cvxpylayers package (Agrawal et al., 2019a ) to construct differentiable
Koopman-MPC policies, we use Proximal Policy Optimization (PPO) (Schulman et
al., (2017)), a prominent RL algorithm for tasks with continuous action
spaces, to learn task-optimal Koopman models. We test our approach on two
distinct case studies derived from a continuous stirred-tank reactor model
from the literature (Flores-Tlacuahuac and Grossmann, (2006); Du et al.,
(2015)): (i) an NMPC case study, wherein the controller objective is to
stabilize the product concentration in the presence of an externally varying
product flow rate, and (ii) a demand response case study, in which eNMPC seeks
to minimize electricity costs without violating constraints on state
variables. We assess the resulting control performance by comparing it to that
of dynamic models trained exclusively with SI as well as model-free policies
trained by RL.

Our findings demonstrate that end-to-end learning consistently outperforms SI
in learning dynamic surrogate models of the CSTR model in our control
applications. Additionally, we investigate to which extent the learned models
and policies still function in the case of modified constraints. Such a change
in the control setting is an example of so-called distribution shift (Zhang et
al., (2021)), which is a major challenge in machine learning applications
(Quinonero-Candela et al., (2008); Zhang et al., (2021); Yao et al., (2022)).
We show that, unlike model-free policies, MPCs employing end-to-end learned
dynamic surrogate models may successfully adapt to such a change in the
control setting without a need for re-training. To our knowledge, this work is
the first to connect end-to-end learning of Koopman models to control
applications with constraints on system variables and the first to study the
ability of end-to-end learned dynamic models to adapt to a distribution shift
in control.

The remainder of this paper is organized as follows: Section 2 provides the
theoretical background to our work and presents our method. Section 3 provides
numerical examples on a case study. Section 4 discusses the conclusions and
directions for future work.

##  2 Method

Subsections 2.1 \- 2.3 provide the necessary theory upon which we build our
approach for end-to-end RL of dynamic Koopman surrogate models for MPC
applications (Subsection 2.4).

###  2.1 Koopman theory for control

The objective of applied Koopman theory is to obtain globally valid linear
representations of nonlinear dynamics by lifting the system states into a
higher dimensional Koopman space through a nonlinear coordinate transformation
(Koopman, (1931)). Multiple ways to extend Koopman theory to controlled
systems have been proposed (e.g., Proctor et al., (2018); Williams et al.,
(2016); Korda and Mezić, (2018)). We use the approach by Korda and Mezić,
(2018) as it results in convex OCPs, under the easily verifiable condition
that the objective function and inequality constraint functions are convex.
The resulting models are of the form

$\displaystyle\bm{z}_{0}$ $\displaystyle=\bm{\psi}_{\bm{\theta}}(\bm{x}_{0}),$
(1a) $\displaystyle\bm{z}_{t+1}$
$\displaystyle=\bm{A}_{\bm{\theta}}\bm{z}_{t}+\bm{B}_{\bm{\theta}}\bm{u}_{t},$
(1b) $\displaystyle\hat{\bm{x}}_{t}$
$\displaystyle=\bm{C}_{\bm{\theta}}\bm{z}_{t},$ (1c)

where $\bm{z}_{t}\in\mathbb{R}^{N}$ is the lifted latent space vector,
$\hat{\bm{x}}_{t}$ is the prediction of the states
$\bm{x}_{t}\in\mathbb{R}^{n}$, and $\bm{u}_{t}\in\mathbb{R}^{m}$ is the
control input vector. $\bm{A}_{\bm{\theta}}\in\mathbb{R}^{N\times N}$ and
$\bm{B}_{\bm{\theta}}\in\mathbb{R}^{N\times m}$ together form an approximation
to the Koopman operator for an actuated system, and
$\bm{C}_{\bm{\theta}}\in\mathbb{R}^{n\times N}$ is a linear decoder. The
initial condition of the predictor is obtained by passing the initial
condition $\bm{x}_{0}$ through a nonlinear encoder
$\bm{\psi}_{\bm{\theta}}\colon\mathbb{R}^{n}\mapsto\mathbb{R}^{N}$, where
typically $N\gg n$. All model components can be learned from data by adjusting
the model parameters $\bm{\theta}$. When using models of the form given by
Equations (1a) - (1c) as part of MPC, the nonlinear encoding of the initial
state into the higher dimensional Koopman space can be performed before
solving the OCP and, therefore, outside of the scope of the OCP optimizer.
Given an initial Koopman state
$\bm{z}_{0}=\bm{\psi}_{\bm{\theta}}(\bm{x}_{0})$, the OCP optimizer sees only
the linear components of the Koopman model, i.e., $\bm{A}_{\bm{\theta}}$ and
$\bm{B}_{\bm{\theta}}$ for the dynamics, and $\bm{C}_{\bm{\theta}}$ to
evaluate the objective function or constraints in case they depend on state
variables. Thus, all equality constraints imposed by the Koopman model are
linear. Therefore, the resulting OCPs are convex if the objective function and
constraints on control and system state variables are convex.

###  2.2 Deep reinforcement learning with continuous action spaces

Deep RL utilizes nonlinear function approximators, such as deep artificial
neural networks, to learn optimal control policies through sequential feedback
from trial and error actuation (Sutton and Barto, (2018)). The optimal control
problem is typically represented by a Markov Decision Process (MDP), wherein
the next state $\bm{x}_{t+1}$ and reward $r_{t+1}$ depend only on the current
state $\bm{x}_{t}$ and action $\bm{u}_{t}$. The MDP of a specific RL problem
is referred to as the environment of that problem. An episode is defined as a
single run of an agent in an environment, from the initial state to a terminal
state. In engineering applications with continuous action spaces, RL methods
that involve agents learning parameterized policies
$\bm{\pi}_{\bm{\theta}}(\bm{u}_{t}|\bm{x}_{t})\colon\mathbb{R}^{n}\mapsto\mathbb{R}^{m}$,
directly mapping states to actions, are the most successful (Sutton and Barto,
(2018)). These methods can be classified as either policy-gradient or actor-
critic methods. In the last decade, actor-critic methods have been the focus
of intense research efforts (Silver et al., (2014); Fujimoto et al., (2018);
Schulman et al., 2015a ; Schulman et al., (2017); Haarnoja et al., (2018)) and
today’s state-of-the-art algorithms typically vastly outperform pure policy-
gradient methods.

###  2.3 Post-optimal sensitivity analysis of convex problems

Obtaining post-optimal sensitivities of solution mappings of optimization
problems has been a topic of study for several decades (e.g., Fiacco and
Kyparisis, (1985); Ralph and Dempe, (1995); Agrawal et al., 2019b ). Given a
parametric nonlinear program with parameters $\bm{\theta}$, and a solution
point $\bm{a}^{*}$ for a given $\bm{\theta}^{*}$, the solution point is
differentiable with respect to the parameters in a neighborhood of
$\bm{\theta}^{*}$ under the following assumptions (Pistikopoulos et al.,
(2020); Fiacco and Ishizuka, (1990)): (i) The objective function and
constraints of the problem are twice continuously differentiable in $\bm{a}$.
(ii) The gradients of objective function and constraints with respect to
$\bm{a}$ and the constraints are once continuously differentiable in
$\bm{\theta}$. (iii) The second-order sufficient conditions (SOSC) for a local
minimum hold at $\bm{a}^{*}$. (iv) The linear independence constraint
qualification (LICQ) is satisfied at $\bm{a}^{*}$. (v) Strict complementary
slackness (SCS) holds at $\bm{a}^{*}$. Note that the above conditions also
imply that the KKT conditions hold (Karush, (1939); Boyd et al., (2004)), and
that $\bm{a}^{*}$ is an isolated local minimum (Fiacco and Ishizuka, (1990)).
For a convex parametric problem this means that the solution map of the
optimization problem becomes single-valued, i.e., the solution map is locally
an implicit function of the problem parameters.

A straightforward but costly way to obtain the derivatives of a solution point
$\bm{a}^{*}$ that satisfies conditions (i-v) with respect to the parameters
$\bm{\theta}$ is to unroll the computational graph of an iterative
optimization algorithm and backpropagate through every step of the entire
algorithm using algorithmic differentiation (e.g., Domke, (2012)).
Alternatively, it is possible to apply the implicit function theorem (Nocedal
and Wright, (1999)) to the KKT conditions of the optimization problem to
obtain backpropagation gradients for a small fraction of the computational
cost required to solve the optimization problem (e.g., Fiacco and Ishizuka,
(1990)). Using the implicit function theorem, Amos and Kolter, (2017) and
Agrawal et al., 2019a  connected solvers for quadratic and convex problems,
respectively, to the popular deep learning frameworks PyTorch (Paszke et al.,
(2019)) and Tensorflow (Abadi et al., (2016)). In doing so, they substantially
lowered the required effort to use post-optimal sensitivity analysis in a deep
learning project.

###  2.4 End-to-end learning of Koopman models for MPC

Figure 1: Workflow from mechanistic model to task-optimal dynamic Koopman
surrogate model.

Figure 1 summarizes our proposed workflow for training task-optimal Koopman
surrogate models. We start with the standard SI procedure for dynamic
surrogate models, given a mechanistic model: First, we generate a data set
$\mathcal{D}$ from which the system dynamics can be learned using supervised
learning, by (randomly) actuating the mechanistic model. Second, we fit a
dynamic surrogate model, in our case a Koopman model with learnable parameters
$\bm{\theta}$, to the data. As described in Section 2.1, our Koopman models
consist of a nonlinear encoder $\bm{\psi_{\theta}}$, the matrices
$\bm{A}_{\bm{\theta}}$ and $\bm{B}_{\bm{\theta}}$ that represent the linear
dynamics inside the Koopman space, and a decoder matrix
$\bm{C}_{\bm{\theta}}$. Lusch et al., (2018) identify three high-level
requirements for autoregressive dynamic Koopman models that correspond to
three types of loss functions that need to be combined when performing system
identification. The extension to Koopman models for controlled systems and the
model structure by Korda and Mezić, (2018) is straightforward and results in
the following requirements and loss functions:

  1. 1.

Identification of nonlinear lifting functions $\bm{\psi_{\theta}}$ that allow
for a linear reconstruction through $\bm{C}_{\bm{\theta}}$. The associated
loss function is the autoencoder loss given by

$||\bm{C_{\theta}}\bm{\psi_{\theta}}(\bm{x}_{t})-\bm{x}_{t}||.$ (2)

  2. 2.

Identification of linear latent space dynamics, using the prediction loss:

$||\bm{A}_{\bm{\theta}}\bm{\psi_{\theta}}(\bm{x}_{t})+\bm{B}_{\bm{\theta}}\bm{u}_{t}-\bm{\psi_{\theta}}(\bm{x}_{t+1})||$
(3)

  3. 3.

State prediction in the original system state, combining all elements of the
model. The associated loss function is given by

$||\bm{C_{\theta}}(\bm{A}_{\bm{\theta}}\bm{\psi_{\theta}}(\bm{x}_{t})+\bm{B}_{\bm{\theta}}\bm{u}_{t})-\bm{x}_{t+1}||.$
(4)

The norm $||\cdot||$ is the mean-squared error, averaging over the number of
samples. Depending on the system, it can be beneficial to use multi-step
versions of the prediction loss (Equation 3) and the combined loss (Equation
4) to improve model accuracy when predicting over multiple time steps in
closed-loop fashion. We follow standard SI practices, such as using a separate
validation data set for early stopping and normalizing the inputs and outputs
of the model.

The main contribution of this work is our method for the third and final step
depicted in Figure 1, the end-to-end refinement of Koopman models for optimal
performance as part of MPC in a specific control task.

Figure 2: Method for end-to-end refinement of dynamic Koopman surrogate model.
The RL agent consists of a stochastic actor and a critic. The actor is an MPC
policy utilizing a dynamic Koopman surrogate model. The critic is a
feedforward neural network. The environment consists of the mechanistic model
of the system that is to be controlled, and a reward function that depends
upon the controllers task.

Figure 2 summarizes our approach towards end-to-end refinement of Koopman
models for optimal performance in (e)NMPC applications through RL. We use the
actor-critic PPO algorithm (Schulman et al., (2017)) to learn the optimal
model parameters $\bm{\theta}$ for a given control task. Note that the
refinement applies to all elements of the Koopman model, i.e., the predictor
consisting of $\bm{A}_{\bm{\theta}}$ and $\bm{B}_{\bm{\theta}}$, and the
autoencoder given by $\bm{\psi_{\theta}}$ and $\bm{C}_{\bm{\theta}}$. Using
the cvxpylayers package (Agrawal et al., 2019a ), we construct automatically
differentiable stochastic Koopman-MPC policies
$\bm{\pi}_{\bm{\theta}}(\bm{u}_{t}|\bm{x}_{t})\colon\mathbb{R}^{n}\mapsto\mathbb{R}^{m}$
that serve as actors. These policies take a state $\bm{x}_{t}$ as input and
compute the Koopman state $\bm{z}_{t}$ through a multilayer perceptron encoder
$\bm{\psi_{\theta}}$. Then, by evaluating objective function and constraints
through the decoder, any (e)MPC problem with a convex objective function and
convex constraints (in state space $\bm{x}$) can be formulated and solved
using cvxpylayers to obtain the optimal control solution $\bm{u}^{*}_{t}$. As
is common in applying MPC and suitable for an RL approach, our MPC policies
are implemented using a receding horizon approach and only return the control
solution for the first control step in the MPC horizon. Given the current time
step $t$ and an MPC horizon of $t_{f}+1$ steps with the corresponding sets
$\mathrm{T}_{+1}=\\{t,\dots,t+t_{f}\\}$ and
$\mathrm{T}=\\{t,\dots,t+t_{f}-1\\}$, the OCPs that are solved inside the
Koopman-MPC policies take the form

$\displaystyle\underset{(\bm{u}_{t})_{t\in\mathrm{T}}}{\operatorname{min}}$
$\displaystyle\sum_{t\in\mathrm{T}_{+1}}\Phi(\bm{C_{\theta}}\bm{z}_{t},\bm{u}_{t})+M\bm{s}^{\intercal}_{t}\bm{s}_{t},$
(5a) $\displaystyle\text{s.t. }\bm{z}_{t+1}$
$\displaystyle=\bm{A_{\theta}}\bm{z}_{t}+\bm{B_{\theta}}\bm{u}_{t}\quad\forall
t\in\mathrm{T},$ (5b) $\displaystyle\underline{\bm{x}}-\bm{s}_{t}$
$\displaystyle\leq\bm{C_{\theta}}\bm{z}_{t}\leq\bar{\bm{x}}+\bm{s}_{t}\quad\forall
t\in\mathrm{T}_{+1},$ (5c) $\displaystyle\underline{\bm{u}}$
$\displaystyle\leq\bm{u}_{k}\leq\bar{\bm{u}}\quad\forall t\in\mathrm{T},$ (5d)
$\displaystyle\bm{g}(\bm{C_{\theta}}\bm{z}_{t},\bm{u}_{t},\bm{s}_{t})$
$\displaystyle\leq 0\quad\forall t\in\mathrm{T}_{+1}.$ (5e)

Both state variables $\bm{x}_{t}$ and control variables $\bm{u}_{t}$ can have
upper and lower bounds (see Equations (5c) and (5d)). To guarantee the
feasibility of the OCPs, we add slack variables $\bm{s}_{t}$ to the state
bounds and penalize their use quadratically in the objective function
(Equation (5a)) with a penalty factor $M$. The stage cost $\Phi$ of the
objective function can be any convex function and depends on the aim of the
policy. The linear transition function in Koopman space is given by Equation
(5b). Depending on the application, the OCPs can also include further convex
inequality constraints (see Equation 5e).

As PPO is an on-policy RL method, the actor must be stochastic during training
to ensure exploration. Therefore, during training, we sample the action that
the actor will take from a normal distribution
$\bm{u}_{t}\sim\bm{\mathcal{N}}(\bm{u}^{*}_{t},\bm{\sigma}^{2})$ with the mean
$\bm{u}^{*}_{t}$ being the deterministic solution of the Koopman-MPC (see
Figure 2). Depending on the application, it can be beneficial to use a
constant standard deviation $\bm{\sigma}$, or start with a large standard
deviation and gradually reduce the amount of exploration, e.g., by reducing
$\bm{\sigma}$ linearly up to a specific value, as the policy becomes more
competent (e.g., Chen et al., (2019)).

We use a critic in the form of a multilayer perceptron (MLP) with learnable
parameters $\bm{\phi}$ and generalized advantage estimation (Schulman et al.,
2015b ) to calculate the advantage estimates during training. These are used
to compute the actor loss via the clipped PPO loss function and to update the
actor parameters $\bm{\theta}$ accordingly (Schulman et al., (2017)). We also
use global gradient clipping to stabilize the training, as described by
Engstrom et al., (2020). We use the Adam optimizer (Kingma and Ba, (2014)) to
train both the actor and the critic.

We use the running average cumulative reward of a fixed number of preceding
episodes as a performance metric during training. We choose the number of
episodes to ensure a relatively small variance while limiting the influence of
outdated policy parameters due to training updates (cf. Sections 3.3 and 3.4).
At the end of a training run, we choose the agent that obtained the highest
average cumulative reward as the final agent. We then evaluate the performance
of that agent by testing it without any exploration, i.e.,
$\bm{u}_{t}=\bm{u}^{*}_{t}$.

Due to the convex nature of the OCPs, our Koopman-MPC policies are not only
automatically differentiable, but also computationally cheap compared to many
alternatives that utilize general nonlinear dynamic models. We expect that
Koopman-MPC policies are real-time feasible for a potentially vast range of
systems.

##  3 Numerical experiments

###  3.1 Case study description

We base our numerical examples on a dimensionless benchmark continuous
stirred-tank reactor (CSTR) model (Flores-Tlacuahuac and Grossmann, (2006); Du
et al., (2015)). The dynamics of the CSTR are given by two nonlinear ordinary
differential equations:

$\displaystyle\dot{c}(t)=$
$\displaystyle\>(1-c(t))\dfrac{\rho(t)}{V}-c(t)ke^{-\frac{N}{T(t)}},$ (6a)
$\displaystyle\begin{split}\dot{T}(t)=&\>(T_{f}-T(t))\dfrac{\rho(t)}{V}+c(t)ke^{-\frac{N}{T(t)}}\\\
&\>-F(t)\alpha_{c}(T(t)-T_{c})\end{split}$ (6b)

The system states $\bm{x}$ are the product concentration $c$ and the
temperature $T$. The control inputs $\bm{u}$ are the production rate $\rho$
and the coolant flow rate $F$. All parameters are given in Table 1.

Table 1: CSTR model parameters (Flores-Tlacuahuac and Grossmann, (2006); Du et al., (2015)) | symbol | value
---|---|---
volume | $V$ | $20$
reaction constant | $k$ | $300\frac{1}{h}$
activation energy | $N$ | $5$
feed temperature | $T_{f}$ | $0.3947$
heat transfer coefficient | $\alpha_{c}$ | $1.95\cdot 10^{-4}$
coolant temperature | $T_{c}$ | $0.3816$

Based on the CSTR model, we construct both NMPC and eNMPC applications: In
NMPC (Section 3.3), the controller shall stabilize the product concentration
$c$ and temperature $T$ at given target values, given random perturbations to
the production rate $\rho$. The goal in our eNMPC case (Section 3.4) is to
optimize process economics while satisfying process constraints, given
electricity price predictions. Table 2 presents the bounds of the system
states and control inputs as well as the target values of the system states in
NMPC.

Table 2: Lower (lb) and upper (ub) bounds of system states and control inputs, target values of system states in NMPC, and steady-state (ss) values used to evaluate the economic benefit of flexible production in eNMPC. | NMPC (Section 3.3) | eNMPC (Section 3.4)
---|---|---
variable | lb | ub | target | lb | ub | ss
$c$ | - | - | 0.1367 | 0.1231 | 0.1504 | 0.1367
$T$ | - | - | 0.7293 | 0.6 | 0.8 | 0.7293
$\rho$ | $0.8\frac{1}{h}$ | $1.2\frac{1}{h}$ | - | $0.8\frac{1}{h}$ | $1.2\frac{1}{h}$ | $1.0\frac{1}{h}$
$F$ | $0.0\frac{1}{h}$ | $700.0\frac{1}{h}$ | - | $0.0\frac{1}{h}$ | $700.0\frac{1}{h}$ | $390.0\frac{1}{h}$

We use discrete time steps of length $\Delta
t_{\text{discr}}=$15\text{\,}\text{min}$$ and control steps of length $\Delta
t_{\text{ctrl}}=$60\text{\,}\text{min}$$. Given a current state $\bm{x}_{t}$
and control action $\bm{u}_{t}$, we use Equations (6a) - (6b) and the RK45
solver (Dormand and Prince, (1980)) in SciPy (Virtanen et al., (2020)) to
calculate the next state $\bm{x}_{t+1}$ in the environment.

We compare the performances of three different controller types:

  1. 1.

Koopman-SI: MPC controller using a Koopman model to approximate the CSTR
dynamics. The model is trained using SI on a data set generated using the
mechanistic model (Equations (6a) - (6b)) and random control inputs. We
provide a more detailed explanation of the data sampling and SI procedure in
Subsection 3.2.

  2. 2.

Koopman-RL: MPC controller using a Koopman model to approximate the CSTR
dynamics. We use the model of the Koopman-SI controller as an initial guess
and tune it for task-optimal MPC performance using the approach described in
Section 2.4.

  3. 3.

MLP: A model-free neural network policy
$\bm{\pi}_{\bm{\theta}}(\bm{u}_{t}|\bm{x}_{t})$ in the form of a multilayer
perceptron (MLP), trained by RL. Similarly to the Koopman-RL, each action is
sampled from a normal distribution
$\bm{u}_{t}\sim\bm{\mathcal{N}}(\bm{u}^{*}_{t},\bm{\sigma}^{2})$, where the
neural network produces $\bm{u}^{*}_{t}$ given the current state of the
environment.

Matching our proposed workflow (see Figure 1), the remainder of this section
is structured as follows: Subsection 3.2 elaborates on the data sampling
procedure using the mechanistic model and discusses the subsequent SI process
using the Koopman model architecture. Subsections 3.3 and 3.4 describe the
end-to-end refinement of the Koopman models and report the performances of all
three controller types for the NMPC and eNMPC applications, respectively.
Finally, Subsection 3.5 presents our findings on a modified eNMPC, which we
use to study the adaptability of end-to-end refined Koopman MPCs to
modifications of constraints after training has been completed.

###  3.2 Data sampling and system identification

Given randomly sampled control inputs $\bm{u}$ within their respective bounds,
the product concentration $c$ typically quickly moves away from its feasible
range in eNMPC and also away from its target value in NMPC (see Table 2).
Therefore, such a sampling strategy would not produce a data set particularly
suitable for learning a good dynamic surrogate model, as most data would be in
irrelevant parts of the state space. Instead of actuating the mechanistic
model by randomly sampling both control inputs $\rho$ and $F$, we, therefore,
take a different approach: We generate multiple trajectories using the
mechanistic model, each trajectory starting from the target state of the NMPC,
which is also at the center of the feasible range of the eNMPC. For each
trajectory, we randomly generate a series of control inputs for either $\rho$
or $F$. Given this series of inputs, we then solve an OCP using the
mechanistic model to determine a corresponding series of values for the
respective other control variable, which leads to $c$ not diverging much from
its feasible range in eNMPC. As in NMPC and eNMPC, we use control steps of
length $\Delta t_{\text{ctrl}}=$60\text{\,}\text{min}$$, i.e., $\rho$ and $F$
change every four time steps.

Given an OCP horizon of $t_{f}$ time steps with the corresponding set
$\mathrm{T}=\\{0,\dots,t_{f}-1\\}$, and a randomly generated trajectory for
the production rate $\rho_{\text{target},t}\forall t\in\mathrm{T}$, the OCP
minimizes

$\sum_{t\in\mathrm{T}}(c_{\text{target}}-c_{t})^{2}+w_{\rho}(\rho_{\text{target},t}-\rho_{t})^{2},$
(7)

subject to the mechanistic model (Equations (6a) - (6b)), and subject to the
given bounds on control variables (see Table 2), with $c_{\text{target}}$
being the target value for $c$ in NMPC (see Table 2) and $w_{\rho}=10$.
Equivalently, given a sampled trajectory for the coolant flow rate
$F_{\text{target},t}\forall t\in\mathrm{T}$, the OCP minimizes

$\sum_{t\in\mathrm{T}}(c_{\text{target}}-c_{t})^{2}+w_{F}\cdot(F_{\text{target},t}-F_{t})^{2},$
(8)

with $w_{F}=0.1$. We choose $w_{\rho}$ and $w_{F}$ so that the resulting data
set has good coverage of the feasible part of the state space (see Table 2,
eNMPC). We use the GEKKO Optimization Suite (Beal et al., (2018)) with default
settings to formulate and solve the OCPs. Given the resulting trajectories for
the control inputs $\bm{u}$, we use the RK45 solver (Dormand and Prince,
(1980)) in SciPy (Virtanen et al., (2020)) to compute the system response and
generate the data set for SI. We generate 84 trajectories for training and
validation purposes. Each of those trajectories has a length of 5 days, i.e.,
480 time steps. We use 63 of those trajectories for training and the remaining
21 for validation.

We perform SI by minimizing the sum of three loss functions (see Equations 2
\- 4), corresponding to the three requirements to a dynamic Koopman model (see
Section 2.4). To obtain Koopman models that perform well over multiple time
steps of closed-loop prediction, we take a stochastic curriculum learning
approach (Bengio et al., (2015); Zaremba and Sutskever, (2014)), i.e., we
start the training procedure by minimizing the mean-squared error (MSE) of
predictions over a single time step, and, progress onto minimizing the average
MSE over multiple time steps of continuous closed-loop prediction in later
epochs. Such a curriculum learning procedure is designed to result in models
that make accurate long-term predictions but to avoid exploding/vanishing
gradients resulting from concatenating predictions of an untrained model in
the first epochs. Specifically, in the first epoch, we use the one-step MSE as
the loss function. After 250 epochs, we use the average MSE over 240 time
steps, i.e., 2.5 days, of closed-loop prediction in every epoch. Between the
first and the 250th epoch, we decrease the probability of taking the one-step
loss linearly from $100\text{\,}\%$ to $0\text{\,}\%$ in favor of the 240-step
loss. We use the Adam optimizer (Kingma and Ba, (2014)) with a learning rate
of $0.5\cdot 10^{-4}$, a mini-batch size of 64 samples, and a maximum number
of 5,000 epochs. We stop training early after at least 350 epochs if the
validation loss has not reached a new minimum for 100 consecutive epochs.

A pragmatic and common practice when training ML models is to oversize the
models and then use regularization techniques and early stopping to prevent
overfitting (Goodfellow et al., (2016)). However, we aim to obtain models that
strike a good balance between accuracy and computational tractability.
Therefore, we avoid unnecessary oversizing of our Koopman surrogate models by
starting with a small model and repeatedly training larger models until the
performance gains become negligible. All results shown in the following are
obtained using a Koopman model with a latent space dimensionality of eight,
i.e., $\bm{A}_{\bm{\theta}}\in\mathbb{R}^{8\times
8},\bm{B}_{\bm{\theta}}\in\mathbb{R}^{8\times
2},\bm{C}_{\bm{\theta}}\in\mathbb{R}^{2\times 8}$, and an MLP encoder
$\bm{\psi}\colon\mathbb{R}^{2}\mapsto\mathbb{R}^{8}$ with two hidden layers
(four and six neurons, respectively), and hyperbolic tangent activation
functions.

###  3.3 NMPC

Here, the controller task is to stabilize the product concentration $c$ and
temperature $T$ at constant target values (see Table 2). Every eight hours,
the production rate $\rho$ is set to a random value within its bounds for the
next eight hours. The controller can respond by varying the coolant flow rate
$F$ within its bounds to stabilize $c$ and $T$. Thus, the production rate
becomes an external input $\rho_{\text{ext}}$ and the coolant flow rate $F$ is
the control input.

During RL, we define the reward $r_{t}$ for each control action as
proportional to the negative squared relative deviation between $c$, $T$, and
their respective steady-state values:

$\begin{split}r_{t}=&\>-((c_{t}-c_{\text{target}})/(c_{\text{ub}}-c_{\text{lb}}))^{2}\\\
&\>-((T_{t}-T_{\text{target}})/(T_{\text{ub}}-T_{\text{lb}}))^{2}\end{split}$
(9)

Relative deviations are achieved by normalizing each variable using the size
of its feasible range (see Table 2). Aptly, using the feasible range of each
variable, we also normalize the inputs and outputs of the Koopman models to
values between zero and one. At each control step, the MPC controller solves
an OCP in Koopman space, minimizing the squared distance from the target state
$\bm{x}_{\text{target}}$:

$\displaystyle\underset{(F_{t})_{t\in\mathrm{T}}}{\operatorname{min}}$
$\displaystyle\sum_{t\in\mathrm{T}_{+1}}||C_{\bm{\theta}}\bm{z}_{t}-\bm{x}_{\text{ss}}||,$
(10a) $\displaystyle\text{s.t. }\bm{z}_{0}$
$\displaystyle=\bm{\psi}_{\bm{\theta}}(\bm{x}_{0}),$ (10b)
$\displaystyle\bm{z}_{t+1}$
$\displaystyle=A_{\bm{\theta}}\bm{z}_{t}+B_{\bm{\theta}}\begin{pmatrix}\rho_{\text{ext}}\\\
F_{t}\end{pmatrix}\;$ $\displaystyle\forall t\in\mathrm{T},$ (10c)
$\displaystyle\underline{F}_{t}$ $\displaystyle\leq F_{t}\leq\bar{F}_{t}\;$
$\displaystyle\forall t\in\mathrm{T}.$ (10d)

We use an MPC horizon of $t_{f}=12$ quarter hours, i.e., three hours, and
episodes with a duration of 288 time steps, i.e., three days. As the
production rate $\rho_{\text{ext}}$ is adjusted every eight hours, each
episode includes eight step changes.

We repeat the RL training of the Koopman-RL and MLP controllers over ten
random seeds. The values of training hyperparameters are given in the
supplementary material. The MLP controllers have two fully connected hidden
layers, each with 256 neurons, followed by an output layer with one neuron for
$F$. All layers have hyperbolic tangent activation functions. As a performance
indicator during the training process, we use the running average sum of
rewards during the previous 30 episodes, i.e., 2,160 control steps. As the
policies are updated every 2,048 control steps, taking the running average of
30 episodes has a relatively low variance while not being overly strongly
influenced by actions resulting from outdated policy parameters. In the
following, we call the undiscounted cumulative reward of an episode the score
of that episode.

Figure 3: Summary of training runs in NMPC. Each training configuration is run
ten times with different random seeds. We average the running score over the
last 30 episodes. The lines represent individual training runs. Additionally,
we highlight the highest average running score achieved by a policy type (star
for Koopman-RL, dot for MLP).

Figure 3 summarizes the RL training results for both RL-trained controller
types. We show the running average of training scores of each training run.
Due to using a pretrained Koopman model, the Koopman-RL controllers start the
RL training at comparatively high scores. Still, after around 5,000 episodes,
most Koopman-RL controllers have significantly increased their average scores
compared to the SI baseline, and the performance reaches a plateau. Some
training runs show rather steadily improving scores from the start, reaching
the plateau after about 1,000 episodes. Other trainings show an initial drop
in performance, followed by a recovery. The MLP controllers are randomly
initialized and obtain low scores at the beginning of the training. The
performance improves gradually and plateaus after about 3,000 episodes at
average scores below those of the SI-pretrained Koopman controllers.

After completing all training runs, we test the policies that achieved the
highest average running scores without exploration, i.e., we set
$\bm{u}_{t}=\bm{u}^{*}_{t}$, on 100 test episodes. Table 3 shows the results
of this test. A direct comparison of the three controllers, given the same
randomly generated production rate trajectory, is given in Figure 4. Note that
as exploration is turned off in these tests, the performance may be higher
than during training (cf. Figure 3).

Table 3: NMPC: test results over 100 episodes. Scores represent the negative summed squared relative deviation from the target values, thus higher is better and a score of 0 is the upper performance bound. Score | avg. | std. | min. | max.
---|---|---|---|---
MLP | -2.21 | 0.73 | -4.40 | -0.43
Koopman-SI | -1.35 | 0.62 | -3.06 | -0.29
Koopman-RL | -0.50 | 0.19 | -1.27 | -0.12
Figure 4: NMPC: comparison of the controller behavior, given the same randomly
generated production rate trajectory. Best viewed in color.

All three controllers exhibit sensible behavior. However, Table 3 shows that
the Koopman-RL controller manages to keep the state of the CSTR on average
closer to the target than the other two controllers. This finding is confirmed
by Figure 4, which shows that the Koopman-RL controller manages to stabilize
the product concentration $c$ quicker and much closer to the target than the
other two controllers. The MLP controller does not stabilize $c$ at or near
the target; instead, it makes $c$ oscillate around its target. Regarding the
temperature variable $T$, all three controllers produce an average deviation
from $T_{\text{target}}$ that is substantially larger than the average
deviation of $c$ from $c_{\text{target}}$. This has two reasons: Firstly, the
feasible range of $T$ is larger than that of $c$ (see Table 2), which affects
the reward calculation (Eq. 9) as we normalize deviations using the feasible
range of each variable. Secondly, the objective function in the Koopman OCPs
(Eq. 10a) is also affected due to the normalization of the inputs of the
Koopman models. The overall performance regarding $T$ is very similar for all
three controllers. Again, the $MLP$ controller produces oscillating behavior,
whereas both Koopman controllers establish steady states at some distance to
the target after each step change.

###  3.4 eNMPC

Here, the controller task is to optimize process economics through demand
response, i.e., shifting electricity consumption in time to benefit from
intervals with relatively low prices, given electricity price profiles. The
state variables $c$ and $T$ are subject to box constraints, with upper and
lower bound values given in Table 2. The control inputs
$\bm{u}_{t}=(\rho_{t},F_{t})^{\intercal}$ can be varied freely within their
bounds. For simplicity, we assume that electric power consumption is
proportional to the coolant flow rate $F$. Thus, the controller can improve
the process economics by shifting process cooling to low-price intervals. To
enable flexible operation, we assume a product storage with filling level
$l_{t}$ and a maximum capacity of six hours of steady-state production.
Because the number of variables in an OCP scales with the prediction horizon,
RL training with a differentiable MPC controller and a long time horizon can
quickly become computationally infeasible. We identify a price prediction
horizon that results in a performant eNMPC by repeatedly solving eNMPC control
problems using the mechanistic model (Equations (6a) - (6b)). Based on the
results, we choose a price prediction horizon of nine hours, as longer
horizons yield no substantial improvement in control performance.

We use historic day-ahead electricity prices from the Austrian market (Open
Power System Data, (2020)) from March 29, 2015 to March 25, 2018 to train the
controllers. The final testing is performed on Austrian electricity price data
from March 26, 2018 to September 30, 2018. Our goal is to train controllers
that perform well when used continuously over an arbitrarily long time
horizon. In addition, when testing controllers on short episodes in the eNMPC
setting, the economic performance of a controller can be heavily influenced by
how the controller utilizes the product storage. For example, a controller
that fills up the storage initially to obtain flexibility might incur higher
costs in a short episode than a steady-state production but might achieve
economic savings in the long run. Therefore, we test the final performance of
the controllers in a single continuous episode that covers half a year.
However, during training, it makes sense to use shorter episodes to allow for
more frequent restarts. This is especially important at the beginning of the
training procedure when controllers are likely to steer the system to
undesirable states from which little can be learned. Therefore, as in the NMPC
case study, we choose an episode length of 288 control steps, i.e., three
days, for training. At the start of each training episode, the state variables
get reset to their steady-state values, the storage level is initialized
randomly between one and two hours of steady-state production volume, and a
connected series of electricity prices of length $\mathrm{T}$ is sampled from
the data. In the roughly half-year-long test episode, we start with an empty
storage to have a consistent starting point for the three different
controllers. The product concentration $c$ and the temperature $T$ are
initialized at their steady-state values.

For RL, we choose a reward that depends on the electricity cost savings
compared to steady-state production (see Table 2) and, on whether any bounds
on state variables were violated during the previous control interval. First,
we check for constraint violations, and compute the steady-state costs

$c_{\text{ss}}=F_{\text{ss}}\cdot p_{t}\cdot\Delta t_{\text{ctrl}},$ (11)

and the costs induced by the current control action

$c_{t}=F_{t}\cdot p_{t}\cdot\Delta t_{\text{ctrl}},$ (12)

using the electricity price $p_{t}$. Then, we calculate the reward

$r_{t}=\begin{cases}-1,&\text{if constraint violation}\\\
(c_{\text{ss}}-c_{t})\beta,&\text{otherwise},\end{cases}$ (13)

wherein $\beta$ is a hyperparameter used to balance the influence of cost
savings and constraint violations.

At each control step, the eMPC controller solves an OCP in Koopman space,
aiming to minimize the operational costs while satisfying all constraints:

$\displaystyle\underset{(\rho_{t},F_{t})_{t\in\mathrm{T}}}{\operatorname{min}}\sum_{t\in\mathrm{T}_{+1}}(F_{t}p_{t}\Delta
t_{\text{ctrl}}+M\bm{s}^{\intercal}_{t}\bm{s}_{t}),$ (14a)
$\displaystyle\text{s.t. }\bm{z}_{0}$
$\displaystyle=\bm{\psi}_{\bm{\theta}}\begin{pmatrix}c_{0}\\\
T_{0}\end{pmatrix},$ (14b) $\displaystyle\bm{z}_{t+1}$
$\displaystyle=\bm{A}_{\bm{\theta}}\bm{z}_{t}+\bm{B}_{\bm{\theta}}\bm{u}_{t}\quad\forall
t\in\mathrm{T},$ (14c) $\displaystyle l_{t+1}$
$\displaystyle=l_{t}+(\rho_{t}-\rho_{\text{ss}})\Delta
t_{\text{discr}}\quad\forall t\in\mathrm{T},$ (14d) $\displaystyle\bm{x}_{t}$
$\displaystyle=\bm{C}_{\bm{\theta}}\bm{z}_{t}\quad\forall
t\in\mathrm{T}_{+1},$ (14e)
$\displaystyle\underline{\bm{x}}_{t}-\bm{s}_{\bm{x},t}$
$\displaystyle\leq\bm{x}_{t}\leq\bar{\bm{x}}_{t}+\bm{s}_{\bm{x},t}\quad\forall
t\in\mathrm{T}_{+1},$ (14f) $\displaystyle 0-s_{l,t}$ $\displaystyle\leq
l_{t}\leq 6.0+s_{l,t}\quad\forall t\in\mathrm{T}_{+1},$ (14g)
$\displaystyle\bm{s}_{t}$ $\displaystyle=\begin{pmatrix}\bm{s}_{\bm{x},t}\\\
s_{l,t}\end{pmatrix}\quad\forall t\in\mathrm{T}_{+1},$ (14h)
$\displaystyle\bm{0}$ $\displaystyle\leq\bm{s}_{t}\quad\forall
t\in\mathrm{T}_{+1},$ (14i) $\displaystyle\underline{\bm{u}}_{t}$
$\displaystyle\leq\bm{u}_{t}\leq\bar{\bm{u}}_{t}\quad\forall t\in\mathrm{T}.$
(14j)

As explained in Section 2, we introduce slack variables $\bm{s}_{\bm{x},t}$
for the state variables and $s_{l,t}$ for the product storage to ensure
feasibility of the policy but penalize their use in the objective function.

The values of the training hyperparameters are given in the supplementary
material. The MLP controllers have separate input layers for $\bm{x}_{t}$,
$l_{t}$, and the trajectory of future electricity prices. Each of those input
layers is followed by a hidden layer, with 100, 56, and 100 neurons,
respectively. The outputs of those layers are concatenated and passed through
two fully connected layers, each of size 256 neurons. The output layer has two
neurons, one each for $\rho$ and $F$. All layers have hyperbolic tangent
activation functions.

As in NMPC, we repeat the RL training of the Koopman-RL and MLP controllers
over ten random seeds. Again, we use the average running score during the
previous 30 episodes, i.e., 2,160 control steps, as a performance indicator
during the training process. Figure 5 summarizes the RL training results for
both controller types. We show the running average of training scores of each
training run. Unlike in NMPC, the SI pretraining of the Koopman models does
not lead to high initial scores, as the resulting Koopman MPCs frequently
produce minor constraint violations resulting in low rewards (see Equation
13). The performance of the Koopman-RL controllers improves and saturates
after about 2,500 episodes. After 5,000 episodes, the maximum running scores
of the Koopman-RL controllers have reached values between $1.4$ and $3.8$. The
running scores of the MLP controllers rise to around zero after only about 750
episodes. Afterward, the running scores continue to rise slowly, saturating
only after about 20,000 episodes at maximum running scores between $4.6$ and
$5.7$. We stopped the RL refinement of the Koopman-RL controllers at 5,000
episodes because of the high computational cost of the training procedure
compared to the MLP controllers, which is due to the necessity to solve and
backpropagate through OCPs repeatedly. Because our case study is relatively
small, analyzing training runtimes would give little insight into the expected
runtime for a practical control problem. Thus, we defer this analysis for
future work on more realistic systems.

Figure 5: Summary of RL training runs in eNMPC. Each training configuration is
run ten times with different random seeds. We average the running score over
the previous 30 episodes. The lines represent individual training runs. The
constant dashed green line represents the median maximum performance when
training the MLP policies for 20,000 episodes. Additionally, we highlight the
highest average running score achieved by a policy type (star for Koopman-RL,
dot for MLP). Best viewed in color.

We observe that convergence to high average scores is unstable, especially for
the Koopman-RL policies. This is a common finding in policy optimization
algorithms (Schulman et al., 2015a ; Schulman et al., (2017)), stemming from
the fact that small perturbations to the parameters $\bm{\theta}$ of a policy
can have outsized effects on the resulting actions, and thus also the rewards.
We hypothesize that this effect is aggravated when learning the parameters of
an MPC policy, as small changes in trajectory predictions can lead to very
different active sets of constraints, resulting in vastly different control
outputs.

Similar to our analysis in Section 3.3, after completing the controller
training runs, we test the policies which achieved the highest ever running
scores during training. We perform this testing without exploration, i.e., we
set $\bm{u}_{t}=\bm{u}^{*}_{t}$, and we report the performance metrics in
Table 4. Additionally, we show the control behavior of the three controllers
in a three-day test episode using a randomly sampled continuous electricity
price profile from the test set in Figure 6 for illustration purposes.

Figure 6: eNMPC: comparison of the controller behavior, given a randomly sampled continuous electricity price trajectory from the test set. The dotted lines represent the bounds (Table 2). Best viewed in color. Table 4: eNMPC: test results over a single roughly half-year long test episode with electricity price data from March 26, 2018 to September 30, 2018. Costs are stated relative to steady-state production. The constraint violation column reports the percentage of control steps that result in constraint violations. | Cost | Constr. viols. [%] | Avg. storage level [h]
---|---|---|---
MLP | 0.91 | 0.24 | 0.80
Koopman-SI | 0.90 | 8.84 | 0.28
Koopman-RL | 0.94 | 0.22 | 5.32

Table 4 shows that even though the Koopman-SI controller achieves the highest
cost savings, it does so at the expense of causing many constraint violations.
The Koopman-RL controller and the MLP controller cause few constraint
violations. In case of the Koopman-RL controller, all violations are related
to the bounds of the concentration variable $c$. The MLP controller achieves
higher cost savings than the Koopman-RL controller. Interestingly, all
constraint violations by the MLP controller are violations of the non-
negativity constraint of the product storage. Apparently, the controllers
utilize the product storage quite differently. While the Koopman-SI controller
mostly keeps the storage near-empty, the Koopman-RL controller exhibits the
opposite behavior, generally operating close to the maximum storage capacity
of six hours, which could partially explain the difference in economic
performance compared to the MLP controller. Most importantly, end-to-end
training clearly improved the overall performance of the Koopman MPC in this
application.

####  3.4.1 Analysis of Koopman embedding before and after end-to-end training

A central challenge in applied Koopman theory is identifying a suitable set of
functions for mapping between the original state space and the lifted space
where the dynamics can be approximated linearly with high accuracy. Here, we
analyze how the proposed end-to-end refinement of a Koopman model for task-
optimal performance in eNMPC affects the autoencoder part of the model, which
defines the mapping between the original and lifted domains and which is pre-
trained by SI and then fine-tuned by RL. To this end, we generate a data set
of 100,000 steps in the demand response environment using the full mechanistic
model by applying the Koopman-SI controller. This data set contains 8,302
steps with constraint violations in either $c$ or $T$, and 91,698 steps
without constraint violations. We use this data set to calculate the
autoencoder MSE for both the Koopman-SI and the Koopman-RL model as a measure
for the reconstruction loss. Furthermore, we pass the 100,000 points through
the autoencoders of both models. Then, we calculate the rates at which points
with violations of state bounds get correctly mapped to the exterior of the
feasible region defined by the state bounds (sensitivity). To determine the
specificity, we proceed analogously for the points without constraint
violations.

Table 5: Analysis of autoencoder performance and constraint violation detection rates for the Koopman-SI and Koopman-RL models. The table displays the autoencoder mean-squared error (MSE) as a measure for the reconstruction loss, sensitivity (true positive rate), and specificity (true negative rate) for different combinations of encoder and decoder from Koopman-SI (SI) and Koopman-RL (RL) models. Herein, positive refers to a point with a constraint violation in $c$ or $T$, and negative to a point without a constraint violation. Encoder model |  Decoder model |  Autoencoder MSE |  True pos. (sensitivity) |  True neg. (specificity)
---|---|---|---|---
SI | SI | $3.2\cdot 10^{-5}$ | 0.844 | 0.995
RL | RL | $3.5\cdot 10^{-4}$ | 0.980 | 0.947
SI | RL | $3.5\cdot 10^{-5}$ | 0.844 | 0.995
RL | SI | $4.1\cdot 10^{-4}$ | 0.971 | 0.947

As can be seen from the upper part of Table 5, the end-to-end refinement leads
to a tenfold increase in the reconstruction loss; still, the loss remains at a
relatively low level. Importantly, the autoencoder of the Koopman-RL model
shows a much higher sensitivity than the autoencoder of the Koopman-SI model,
i.e., the rate at which a point that violates a constraint on $c$ or $T$ still
does so after being passed through the autoencoder rises from $84.4\;\%$ to
$98.0\;\%$. When using the Koopman models as part of an eNMPC, this increase
in sensitivity and the related decrease in specificity leads to more
conservative controller behavior. The higher sensitivity and specificity are
expected as the Koopman model was refined using a reward function that heavily
punishes constraint violations (see Eq. 13).

To identify which part of the autoencoder model was modified to what extent
during the refinement, we combine the encoder of one model with the decoder of
the other model, i.e., the encoder of the Koopman-SI model with the decoder of
the Koopman-RL model and vice versa. As can be seen from the lower part of
Table 5, the combinations of encoders and decoders from different models still
yield low reconstruction losses, suggesting that the Koopman models before and
after the refinement operate with very similar lifted spaces. Interestingly,
using the encoder of the Koopman-SI model together with the decoder of the
Koopman-RL model produces almost identical results to using the full Koopman-
SI autoencoder. Conversely, using the Koopman-RL encoder in conjunction with
the Koopman-SI decoder produces very similar results to using the full
Koopman-RL autoencoder. These observations indicate that the end-to-end
refinement affected the encoder to a larger extent than the decoder.

###  3.5 eNMPC with adapted bounds

One key advantage of MPC controllers over model-free policies is that they may
be able to adapt to changing control settings, such as shifted bounds or a
change in the objective function. Model-free policies cannot adapt to such
changes, as constraints and objective function are learned implicitly by the
policy through the reward signal of the environment. Therefore, any change to
the environment would require a retraining of the policy.

We shift the bounds of the product concentration variable $c$ and rerun the
test (cf. Subsection 3.4) without retraining. We examine three different
constraint adaptions: relaxed bounds, tightened bounds, and a shifted feasible
region of $c$, i.e., a different target product. Table 6 shows the adapted
bounds of $c$ in all three cases and reports the frequency of constraint
violations that result from applying the controllers. Additionally, Figure 7
illustrates the resulting trajectories of $c$ when using the controllers on a
three-day test episode for the different bound adaptations. We waive reporting
the production costs in Table 6 as they become meaningless in the presence of
vast differences in the rates of constraint violations.

Table 6: Adapted lower (lb) and upper (ub) bounds of $c$ and constraint violations resulting from applying the controllers without retraining. The constraint violation column reports the percentage of control steps that result in constraint violations. | Adapted bounds of $c$ | Constr. viols. [%]
---|---|---
| lb | ub | MLP | Koopman-SI | Koopman-RL
Tightened | 0.1299 | 0.1435 | 49.20 | 19.53 | 1.73
Relaxed | 0.1162 | 0.1572 | 0.24 | 4.41 | 0.13
Shifted | 0.1504 | 0.1777 | 100.00 | 14.94 | 5.83
Figure 7: eNMPC with adapted $c$ bounds: comparison of the controller
behavior, given a randomly sampled continuous electricity price trajectory.
The original bounds used during RL training are shown in dotted black lines,
whereas the shifted bounds are shown in dashed black lines. Best viewed in
color.

Both Table 6 and Figure 7 show that the Koopman-RL controller is best at
adapting to a change in the control setting, even though RL trained it for
optimal performance in a somewhat different setting. Even though the tightened
bounds and shifted bounds cases lead to more frequent constraint violations
than in the original setting (cf. Table 4), the Koopman-RL controller
apparently still profits from its end-to-end refinement, as constraint
violations are markedly rarer than for the Koopman-SI controller. In the
tightened bounds and relaxed bounds cases, the relative difference in the
frequency of constraint violations compared to the Koopman-SI controller stays
similar to the original case (cf. Table 4). In the shifted bounds case, the
Koopman-RL controller produces more constraint violations than in the other
two cases, which is expected, as the RL training discouraged ever entering the
shifted bounds range.

Unsurprisingly, the MLP controller is agnostic to a change in the control
setting and fails in the tightened bounds and shifted bounds cases. In the
relaxed bounds case, it causes few constraint violations due to its unaltered
behavior. Again, all constraint violations are violations of the non-
negativity constraint of the product storage, which we did not relax. As the
MLP controller respects the original training bounds of the state variables
$c$ and $T$, it does not take advantage of the enhanced flexibility potential
of the process in the relaxed bounds case.

##  4 Conclusion

This article presents a method for end-to-end RL of Koopman models for optimal
performance in (e)NMPC applications. We exploit the modeling approach by Korda
and Mezić, (2018) to construct Koopman-MPC controllers with convex underlying
OCPs, and cvxpylayers (Agrawal et al., 2019a ) for post-optimal sensitivity
analysis, thus turning the Koopman-MPC controllers into automatically
differentiable policies that can be trained using state-of-the-art RL methods
like PPO (Schulman et al., (2017)). Applying a trained Koopman-MPC controller
has a comparatively low computational cost due to the convexity of the OCPs.
Therefore, our method should be able to produce predictive controllers that
are real-time capable for potentially large control systems.

We demonstrate the effectiveness of our method by comparing the control
performance to those of MPC utilizing models trained by system identification
and model-free neural network controllers trained by RL. Specifically, we
investigate two control applications derived from a CSTR model: (i) a setpoint
tracking problem (NMPC, Section 3.3), and (ii) a demand response problem with
hard constraints on state variables (eNMPC, Section 3.4). Our RL-trained
Koopman models consistently outperform their counterparts trained by system
identification.

Additionally, we alter the feasible region of a state variable in eNMPC and
test the controllers without retraining (Section 3.5). Such changes in the
control setting present a significant challenge to data-driven control. Due to
the explicit constraints in the Koopman-eNMPCs, end-to-end training using RL
preserves the controllers inherent adaptability to changes in the environment.
Thus, even though the Koopman-eNMPCs lose some of the gains in performance
from the RL refinement procedure, they can adapt better to the changes to the
environment than the Koopman-eNMPCs trained by system identification alone. As
expected, the neural network controllers cannot adapt as they implicitly learn
the constraints of their environment through the rewards received during
training.

We validated our method on nonlinear but relatively small case studies, far
from the complexity of many real-world systems. Future work should therefore
test the approach on more complex systems and, eventually, real-world systems.
However, applying our method to more complex control problems would likely
necessitate training for more iterations. Additionally, it would increase the
computational burden associated with each iteration due to larger Koopman
surrogate models and more expensive environment steps. Both factors would
increase the overall computational burden of the training process and could
render our approach computationally intractable for very large problems.
Future work should therefore explore ways to increase the data efficiency and
stability of the learning process. To this end, differentiable Koopman (e)NMPC
could be combined with model-based RL algorithms (Janner et al., (2019)) or
policy optimization algorithms that exploit the differentiability of
environments that are based on mechanistic simulation models (Xu et al.,
(2022); Mayfrank et al., (2024)). Another avenue of possible future research
is combining RL and differentiable Koopman (e)NMPC to learn task-optimal
disturbance estimators (Son et al., 2022b ) instead of refining the dynamic
model using RL.

## Declaration of Competing Interest

We have no conflict of interest.

## Acknowledgements

This work was performed as part of the Helmholtz School for Data Science in
Life, Earth and Energy (HDS-LEE) and received funding from the Helmholtz
Association of German Research Centres.

We thank the anonymous reviewer for raising the idea of using differentiable
Koopman (e)NMPC to learn task-optimal disturbance estimators for offset-free
MPC (see Conclusion). We also want to thank our former ICE-1 colleagues
Danimir T. Doncevic and Florian J. Baader for fruitful discussions and
literature recommendations.

## Declaration of generative AI and AI-assisted technologies in the writing
process

During the preparation of this work Daniel Mayfrank used ChatGPT and Grammarly
in order to correct grammar and spelling and to improve style of writing.
After using these tools, all authors reviewed and edited the content as needed
and take full responsibility for the content of the publication.

## Nomenclature

### Abbreviations

CSTR | continuous stirred-tank reactor
---|---
(e)(N)MPC | (economic) (nonlinear) model
| predictive control
HVAC | heating, ventilation and air
| conditioning
KKT | Karush-Kuhn-Tucker
LICQ | linear independence constraint
| qualification
MDP | Markov decision process
MLP | multilayer perceptron
MSE | mean squared error
OCP | optimal control problem
PPO | proximal policy optimization
RL | reinforcement learning
SCS | strict complementary slackness
SI | system identification
SOSC | second-order sufficient condition

### Greek Symbols

$\alpha$ | learning rate
---|---
$\beta$ | reward calculation hyperparameter
$\gamma$ | reward discount factor
$\epsilon$ | clipping hyperparameter
$\bm{\theta}$ | learnable parameters of actor/predictor
$\lambda$ | generalized advantage estimation
| hyperparameter
$\bm{\mu}$ | expected value for action selection
$\bm{\pi}$ | policy
$\rho$ | CSTR production rate
$\bm{\sigma}$ | standard deviation for action selection
$\bm{\phi}$ | learnable parameters of critic
$\Phi$ | MPC stage cost
$\bm{\psi}$ | encoder MLP

### Latin Symbols

$\bm{a}$ | solution point of an optimization
---|---
| problem
$\bm{A}$ | autoregressive part of Koopman
| dynamics matrix
$\bm{B}$ | external input part of Koopman
| dynamics matrix
$c$ | CSTR product concentration
$\bm{C}$ | decoder matrix of Koopman model
$F$ | CSTR coolant flow rate
$\bm{g}$ | inequality constraints
$K_{\text{PPO}}$ | number of epochs per update
$M$ | penalty factor for slack variables
$M_{\text{PPO}}$ | minibatch size
$\bm{\mathcal{N}}$ | normal distribution
$N_{\text{PPO}}$ | number of parallel actors
$p$ | electricity price
$r$ | reward
$\bm{s}$ | slack variables
$t$ | indicates time-dependency of variable
$T$ | CSTR temperature
$\mathrm{T}$ | set of discrete time steps
$T_{\text{PPO}}$ | control steps between updates to
| actor and critic
$\bm{u}$ | control variables
$\bm{x}$ | system state variables
$\bm{z}$ | Koopman state variables

### Subscripts

ext | indicates external input
---|---
lb | lower bound
ss | steady-state
$t$ | discrete time step
ub | upper bound

### Superscripts

$*$ | indicates optimality
---|---

## Author contributions

Daniel Mayfrank: Conceptualization, Methodology, Software, Investigation,
Writing - original draft, Visualization

Alexander Mitsos: Conceptualization, Writing - review & editing, Supervision,
Funding acquisition

Manuel Dahmen: Conceptualization, Methodology, Writing - review & editing,
Supervision, Funding acquisition

## Bibliography

  * Abadi et al., (2016) Abadi, M., Barham, P., Chen, J., Chen, Z., Davis, A., Dean, J., Devin, M., Ghemawat, S., Irving, G., Isard, M., et al. (2016).  TensorFlow: a system for large-scale machine learning.  In 12th USENIX Symposium on Operating Systems Design and Implementation (OSDI 16), pages 265–283.
  * (2) Agrawal, A., Amos, B., Barratt, S., Boyd, S., Diamond, S., and Kolter, J. Z. (2019a).  Differentiable convex optimization layers.  Advances in Neural Information Processing Systems, 32:9558–9570.
  * (3) Agrawal, A., Barratt, S., Boyd, S., Busseti, E., and Moursi, W. M. (2019b).  Differentiating through a cone program.  arXiv preprint arXiv:1904.09043.
  * Albalawi and Hameed, (2023) Albalawi, F. and Hameed, S. W. (2023).  Koopman-based economic model predictive control for nonlinear systems.  In 2023 9th International Conference on Control, Decision and Information Technologies (CoDIT), pages 822–828.
  * Amos et al., (2018) Amos, B., Jimenez, I., Sacks, J., Boots, B., and Kolter, J. Z. (2018).  Differentiable MPC for end-to-end planning and control.  Advances in Neural Information Processing Systems, 31:8299–8310.
  * Amos and Kolter, (2017) Amos, B. and Kolter, J. Z. (2017).  Optnet: Differentiable optimization as a layer in neural networks.  In International Conference on Machine Learning, pages 136–145.
  * Angeli et al., (2011) Angeli, D., Amrit, R., and Rawlings, J. B. (2011).  On average performance and stability of economic model predictive control.  IEEE Transactions on Automatic Control, 57(7):1615–1626.
  * Arbabi et al., (2018) Arbabi, H., Korda, M., and Mezic, I. (2018).  A data-driven Koopman model predictive control framework for nonlinear flows.  arXiv preprint arXiv:1804.05291.
  * Beal et al., (2018) Beal, L., Hill, D., Martin, R., and Hedengren, J. (2018).  Gekko optimization suite.  Processes, 6(8):106.
  * Bengio et al., (2015) Bengio, S., Vinyals, O., Jaitly, N., and Shazeer, N. (2015).  Scheduled sampling for sequence prediction with recurrent neural networks.  Advances in Neural Information Processing Systems, 28:1171–1179.
  * Boyd et al., (2004) Boyd, S., Boyd, S. P., and Vandenberghe, L. (2004).  Convex Optimization.  Cambridge University Press.
  * Bruder et al., (2019) Bruder, D., Gillespie, B., Remy, C. D., and Vasudevan, R. (2019).  Modeling and control of soft robots using the Koopman operator and model predictive control.  arXiv preprint arXiv:1902.02827.
  * Brunton and Kutz, (2022) Brunton, S. L. and Kutz, J. N. (2022).  Data-driven Science and Engineering: Machine Learning, Dynamical Systems, and Control.  Cambridge University Press.
  * Chen et al., (2019) Chen, B., Cai, Z., and Bergés, M. (2019).  Gnu-RL: A precocial reinforcement learning solution for building HVAC control using a differentiable MPC policy.  In Proceedings of the 6th ACM International Conference on Systems for Energy-Efficient Buildings, Cities, and Transportation, pages 316–325.
  * Cibulka et al., (2020) Cibulka, V., Haniš, T., Korda, M., and Hromčík, M. (2020).  Model predictive control of a vehicle using Koopman operator.  IFAC-PapersOnLine, 53(2):4228–4233.  21st IFAC World Congress.
  * Domke, (2012) Domke, J. (2012).  Generic methods for optimization-based modeling.  In Artificial Intelligence and Statistics, pages 318–326.
  * Dormand and Prince, (1980) Dormand, J. R. and Prince, P. J. (1980).  A family of embedded Runge-Kutta formulae.  Journal of Computational and Applied Mathematics, 6(1):19–26.
  * Du et al., (2015) Du, J., Park, J., Harjunkoski, I., and Baldea, M. (2015).  A time scale-bridging approach for integrating production scheduling and process control.  Computers & Chemical Engineering, 79:59–69.
  * Engstrom et al., (2020) Engstrom, L., Ilyas, A., Santurkar, S., Tsipras, D., Janoos, F., Rudolph, L., and Madry, A. (2020).  Implementation matters in deep policy gradients: A case study on PPO and TRPO.  arXiv preprint arXiv:2005.12729.
  * Fiacco and Ishizuka, (1990) Fiacco, A. V. and Ishizuka, Y. (1990).  Sensitivity and stability analysis for nonlinear programming.  Annals of Operations Research, 27(1):215–235.
  * Fiacco and Kyparisis, (1985) Fiacco, A. V. and Kyparisis, J. (1985).  Sensitivity analysis in nonlinear programming under second order assumptions.  In Systems and Optimization, pages 74–97. Springer.
  * Flores-Tlacuahuac and Grossmann, (2006) Flores-Tlacuahuac, A. and Grossmann, I. E. (2006).  Simultaneous cyclic scheduling and control of a multiproduct CSTR.  Industrial & Engineering Chemistry Research, 45(20):6698–6712.
  * Folkestad and Burdick, (2021) Folkestad, C. and Burdick, J. W. (2021).  Koopman NMPC: Koopman-based learning and nonlinear model predictive control of control-affine systems.  In 2021 IEEE International Conference on Robotics and Automation (ICRA), pages 7350–7356.
  * Fujimoto et al., (2018) Fujimoto, S., Hoof, H., and Meger, D. (2018).  Addressing function approximation error in actor-critic methods.  In International Conference on Machine Learning, pages 1587–1596.
  * Goodfellow et al., (2016) Goodfellow, I., Bengio, Y., and Courville, A. (2016).  Deep Learning.  MIT Press.
  * Gros and Zanon, (2019) Gros, S. and Zanon, M. (2019).  Data-driven economic NMPC using reinforcement learning.  IEEE Transactions on Automatic Control, 65(2):636–648.
  * Haarnoja et al., (2018) Haarnoja, T., Zhou, A., Abbeel, P., and Levine, S. (2018).  Soft actor-critic: Off-policy maximum entropy deep reinforcement learning with a stochastic actor.  In International Conference on Machine Learning, pages 1861–1870.
  * Hussein et al., (2017) Hussein, A., Gaber, M. M., Elyan, E., and Jayne, C. (2017).  Imitation learning: A survey of learning methods.  ACM Computing Surveys (CSUR), 50(2):1–35.
  * Iwata and Kawahara, (2022) Iwata, T. and Kawahara, Y. (2022).  Data-driven end-to-end learning of pole placement control for nonlinear dynamics via Koopman invariant subspaces.  arXiv preprint arXiv:2208.08883.
  * Janner et al., (2019) Janner, M., Fu, J., Zhang, M., and Levine, S. (2019).  When to trust your model: Model-based policy optimization.  Advances in Neural Information Processing Systems, 32:12498–12509.
  * Kalman et al., (1960) Kalman, R. E. et al. (1960).  Contributions to the theory of optimal control.  Boletín de la Sociedad Matemática Mexicana, 5(2):102–119.
  * Karush, (1939) Karush, W. (1939).  Minima of functions of several variables with inequalities as side constraints.  M. Sc. Dissertation. Department of Mathematics, University of Chicago.
  * Kingma and Ba, (2014) Kingma, D. P. and Ba, J. (2014).  Adam: A method for stochastic optimization.  arXiv preprint arXiv:1412.6980.
  * Koopman, (1931) Koopman, B. O. (1931).  Hamiltonian systems and transformation in Hilbert space.  Proceedings of the National Academy of Sciences, 17(5):315–318.
  * Korda and Mezić, (2018) Korda, M. and Mezić, I. (2018).  Linear predictors for nonlinear dynamical systems: Koopman operator meets model predictive control.  Automatica, 93:149–160.
  * Lewis and Vrabie, (2009) Lewis, F. L. and Vrabie, D. (2009).  Reinforcement learning and adaptive dynamic programming for feedback control.  IEEE Circuits and Systems Magazine, 9(3):32–50.
  * Lewis et al., (2012) Lewis, F. L., Vrabie, D., and Vamvoudakis, K. G. (2012).  Reinforcement learning and feedback control: Using natural decision methods to design optimal adaptive controllers.  IEEE Control Systems Magazine, 32(6):76–105.
  * Lusch et al., (2018) Lusch, B., Kutz, J. N., and Brunton, S. L. (2018).  Deep learning for universal linear embeddings of nonlinear dynamics.  Nature Communications, 9(1):1–10.
  * Mayfrank et al., (2024) Mayfrank, D., Ahn, N. Y., Mitsos, A., and Dahmen, M. (2024).  Task-optimal data-driven surrogate models for eNMPC via differentiable simulation and optimization.  arXiv preprint arXiv:2403.14425.
  * McBride and Sundmacher, (2019) McBride, K. and Sundmacher, K. (2019).  Overview of surrogate modeling in chemical process engineering.  Chemie Ingenieur Technik, 91(3):228–239.
  * Mnih et al., (2013) Mnih, V., Kavukcuoglu, K., Silver, D., Graves, A., Antonoglou, I., Wierstra, D., and Riedmiller, M. (2013).  Playing Atari with deep reinforcement learning.  arXiv preprint arXiv:1312.5602.
  * Narasingam and Kwon, (2019) Narasingam, A. and Kwon, J. S.-I. (2019).  Koopman Lyapunov-based model predictive control of nonlinear chemical process systems.  AIChE Journal, 65(11):e16743.
  * Narasingam and Kwon, (2020) Narasingam, A. and Kwon, J. S.-I. (2020).  Application of Koopman operator for model-based control of fracture propagation and proppant transport in hydraulic fracturing operation.  Journal of Process Control, 91:25–36.
  * Nocedal and Wright, (1999) Nocedal, J. and Wright, S. J. (1999).  Numerical Optimization.  Springer.
  * Open Power System Data, (2020) Open Power System Data (2020).  Open power system data.  https://data.open-power-system-data.org/time_series/ (accessed on 2022-08-29).
  * Paszke et al., (2019) Paszke, A., Gross, S., Massa, F., Lerer, A., Bradbury, J., Chanan, G., Killeen, T., Lin, Z., Gimelshein, N., Antiga, L., et al. (2019).  PyTorch: An imperative style, high-performance deep learning library.  Advances in Neural Information Processing Systems, 32:8024–8035.
  * Pistikopoulos et al., (2020) Pistikopoulos, E. N., Diangelakis, N. A., and Oberdieck, R. (2020).  Multi-Parametric Optimization and Control.  John Wiley & Sons.
  * Proctor et al., (2018) Proctor, J. L., Brunton, S. L., and Kutz, J. N. (2018).  Generalizing Koopman theory to allow for inputs and control.  SIAM Journal on Applied Dynamical Systems, 17(1):909–930.
  * Quinonero-Candela et al., (2008) Quinonero-Candela, J., Sugiyama, M., Schwaighofer, A., and Lawrence, N. D. (2008).  Dataset Shift in Machine Learning.  MIT Press.
  * Ralph and Dempe, (1995) Ralph, D. and Dempe, S. (1995).  Directional derivatives of the solution of a parametric nonlinear program.  Mathematical Programming, 70(1):159–172.
  * (51) Schulman, J., Levine, S., Abbeel, P., Jordan, M., and Moritz, P. (2015a).  Trust region policy optimization.  In International Conference on Machine Learning, pages 1889–1897.
  * (52) Schulman, J., Moritz, P., Levine, S., Jordan, M., and Abbeel, P. (2015b).  High-dimensional continuous control using generalized advantage estimation.  arXiv preprint arXiv:1506.02438.
  * Schulman et al., (2017) Schulman, J., Wolski, F., Dhariwal, P., Radford, A., and Klimov, O. (2017).  Proximal policy optimization algorithms.  arXiv preprint arXiv:1707.06347.
  * Silver et al., (2014) Silver, D., Lever, G., Heess, N., Degris, T., Wierstra, D., and Riedmiller, M. (2014).  Deterministic policy gradient algorithms.  In International Conference on Machine Learning, pages 387–395.
  * Son et al., (2021) Son, S. H., Choi, H.-K., and Kwon, J. S.-I. (2021).  Application of offset-free Koopman-based model predictive control to a batch pulp digester.  AIChE Journal, 67(9):e17301.
  * (56) Son, S. H., Choi, H.-K., Moon, J., and Kwon, J. S.-I. (2022a).  Hybrid Koopman model predictive control of nonlinear systems using multiple EDMD models: An application to a batch pulp digester with feed fluctuation.  Control Engineering Practice, 118:104956.
  * (57) Son, S. H., Narasingam, A., and Kwon, J. S.-I. (2022b).  Development of offset-free Koopman Lyapunov-based model predictive control and mathematical analysis for zero steady-state offset condition considering influence of Lyapunov constraints on equilibrium point.  Journal of Process Control, 118:26–36.
  * Sutton and Barto, (2018) Sutton, R. S. and Barto, A. G. (2018).  Reinforcement Learning: An Introduction.  MIT press.
  * Virtanen et al., (2020) Virtanen, P., Gommers, R., Oliphant, T. E., Haberland, M., Reddy, T., Cournapeau, D., Burovski, E., Peterson, P., Weckesser, W., Bright, J., et al. (2020).  SciPy 1.0: fundamental algorithms for scientific computing in Python.  Nature Methods, 17(3):261–272.
  * Williams et al., (2016) Williams, M. O., Hemati, M. S., Dawson, S. T., Kevrekidis, I. G., and Rowley, C. W. (2016).  Extending data-driven Koopman analysis to actuated systems.  IFAC-PapersOnLine, 49(18):704–709.
  * Xu et al., (2022) Xu, J., Makoviychuk, V., Narang, Y., Ramos, F., Matusik, W., Garg, A., and Macklin, M. (2022).  Accelerated policy learning with parallel differentiable simulation.  arXiv preprint arXiv:2204.07137.
  * Yao et al., (2022) Yao, H., Choi, C., Cao, B., Lee, Y., Koh, P. W. W., and Finn, C. (2022).  Wild-time: A benchmark of in-the-wild distribution shift over time.  Advances in Neural Information Processing Systems, 35:10309–10324.
  * Yin et al., (2022) Yin, H., Welle, M. C., and Kragic, D. (2022).  Embedding Koopman optimal control in robot policy learning.  In 2022 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), pages 13392–13399. IEEE.
  * Zaremba and Sutskever, (2014) Zaremba, W. and Sutskever, I. (2014).  Learning to execute.  arXiv preprint arXiv:1410.4615.
  * Zhang et al., (2021) Zhang, M., Marklund, H., Dhawan, N., Gupta, A., Levine, S., and Finn, C. (2021).  Adaptive risk minimization: Learning to adapt to domain shift.  Advances in Neural Information Processing Systems, 34:23664–23678.