# Statistical Proper Orthogonal Decomposition for model reduction in feedback
control

Sergey Dolgov  Department of Mathematical Sciences, University of Bath, United
Kingdom<EMAIL_ADDRESS>Dante Kalise  Department of Mathematics,
Imperial College London, United Kingdom<EMAIL_ADDRESS>Luca Saluzzi
Department of Mathematics, Scuola Normale Superiore, Pisa, Italy.
<EMAIL_ADDRESS>
###### Abstract

Feedback control synthesis for nonlinear, parameter-dependent fluid flow
control problems is considered. The optimal feedback law requires the solution
of the Hamilton-Jacobi-Bellman (HJB) PDE suffering the curse of
dimensionality. This is mitigated by Model Order Reduction (MOR) techniques,
where the system is projected onto a lower-dimensional subspace, over which
the feedback synthesis becomes feasible. However, existing MOR methods assume
at least one relaxation of generality, that is, the system should be linear,
or stable, or deterministic.

We propose a MOR method called Statistical POD (SPOD), which is inspired by
the Proper Orthogonal Decomposition (POD), but extends to more general
systems. Random samples of the original dynamical system are drawn, treating
time and initial condition as random variables similarly to possible
parameters in the model, and employing a stabilizing closed-loop control. The
reduced subspace is chosen to minimize the empirical risk, which is shown to
estimate the expected risk of the MOR solution with respect to the
distribution of all possible outcomes of the controlled system. This reduced
model is then used to compute a surrogate of the feedback control function in
the Tensor Train (TT) format that is computationally fast to evaluate online.
Using unstable Burgers’ and Navier-Stokes equations, it is shown that the SPOD
control is more accurate than Linear Quadratic Regulator or optimal control
derived from a model reduced onto the standard POD basis, and faster than the
direct optimal control of the original system.

##  1 Introduction

Nonlinear dynamics arising in fluid flow play a vital role in numerous
industrial applications, ranging from aerospace engineering to chemical
processing. Efficient management and control of these systems under random
state fluctuations or different parameters is crucial to ensure optimal
performance, stability, and safety. One key approach to achieving these
objectives is feedback control. This involves continuous monitoring of the
system’s behavior and applying corrective actions to minimize deviations from
a reference state. In the context of fluid flow problems, feedback control
enables engineers to actively manipulate various system parameters such as
flow rates, pressures, temperatures, and concentrations to stabilize the
dynamics. In the case of linear dynamics and quadratic cost functions, a
feedback controller that is relatively cheap computationally is the Linear
Quadratic Regulator (LQR), which reduces to the solution of a matrix Algebraic
Riccati equation [1, 2]. However, the complexity and nonlinearity inherent in
fluid dynamics pose significant challenges in the design and implementation of
feedback control strategies. Unlike other fields where linear models and
simplified assumptions may be sufficient, fluid systems often exhibit
intricate behaviors, including turbulence, flow separation, and time-dependent
phenomena. Consequently, developing effective feedback control techniques for
fluid problems turns out to be a very challenging task.

Potentially, the optimal feedback control can be computed by solving the
Hamilton-Jacobi-Bellman (HJB) nonlinear partial differential equation (PDE).
However, the HJB equation is usually of hyperbolic nature, and posed on the
state space of the dynamical system, which can be very high-dimensional. The
latter is typical when the system arises from a discretization of a PDE, such
as the Navier-Stokes equation in fluid dynamics. Moderately high-dimensional
HJB equations have been tackled with various methods: max-plus algebra [3, 4],
sparse grids and polynomials [5, 6, 7], tree-structure algorithms [8, 9], deep
neural networks [10, 11, 12, 13, 14, 15, 16, 17], low-rank tensor
decompositions [18, 19, 20] and kernel interpolation techniques [21, 22].
However, direct solution of the HJB equation for thousands of dimensions
(which is typical in Finite Element discretizations) remains out of reach.

The PDE formulation of the HJB equation can be bypassed by a data-driven
approach, where one computes the value function and/or control at each of the
given (usually random) states in a closed-loop fashion, and finds an
approximate feedback control in some ansatz by solving a regression problem by
minimizing an empirical risk of given samples of the control. Such regression
in a sparse polynomial basis was considered e.g. in [23, 24, 25], a tensor
format regression was employed in [26, 18, 27], and a neural network was
trained in [28, 29, 30, 31]. However, the dimension of the sought control
function is still that of the state space, requiring a sheer number of
unknowns in the approximation ansatz and training data.

Fortunately, not all states of the system are important or reachable,
especially in the controlled regime. Therefore, a promising approach to tackle
the curse of dimensionality is the Model Order Reduction (MOR). Here, the
dynamical system is projected onto some precomputed basis of lower dimension,
and the HJB equation (or a data-driven variant thereof) is solved for the
reduced system. The key question now is the algorithm to compute such a basis.

One of the simplest options is called Proper Orthogonal Decomposition (POD).
This method collects system states (”snapshots”) at certain time points, and
finds a basis that minimizes the total projection error of these snapshots
onto the basis. In the context of feedback control, the POD method was used to
deliver a reduced system (and hence a reduced HJB equation) in [32, 33].
However, the snapshots produced from an uncontrolled deterministic system may
be inaccurate, and even misleading for the controlled system at a different
realisation of random parameters or initial conditions. For instance, it may
be simply impossible to collect snapshots from an uncontrolled system that
exhibits a finite-time blowup.

A state of the art method to identify a basis that approximates well only
controllable and observable states is balanced truncation (BT) [34]. Note that
the original formulation of BT still uses the uncontrolled system, and thus
requires its stability, as well as linearity. This problem was mitigated by
the closed-loop balanced truncation [35] that incorporates the LQR,
stabilizing the system matrix. The linearity assumption aside, the reduction
of the closed-loop instead of uncontrolled system will be the first ingredient
of our paper.

Since BT uses the linear structure of the dynamical system quite explicitly,
its generalization to nonlinear systems is difficult, and still limited to
relatively simple cases. For example, a BT for bilinear systems was introduced
in [36], lifting of some nonlinear systems to quadratic-bilinear systems of
larger dimension was proposed in [37], and an $\mathcal{H}_{\infty}$-balanced
truncation for differential-algebraic systems, linearizing the system around
one state (the stationary solution), was presented in [38]. A more general
approach to nonlinear BT has been recently developed in [39, 40] however, the
nonlinear reduction also requires the solution of high-dimensional HJB type
equations, which is precisely the computational limitation we try to avoid in
this paper. As we will see in numerical examples, linearization around one
particular state may be inaccurate for more general nonlinearities and random
parameters, such as the Navier-Stokes equation with random inflow and boundary
control.

Random parameters (such as random coefficients or initial conditions) are
particularly difficult but relevant scenarios motivating the feedback control,
which needs to be robust with respect to this uncertainty. If the random
variables enter the model (bi)linearly, they can also be lifted to apply e.g.
bilinear BT [41]. Existence of a quadratic-bilinear structure suitable for BT
is not clear in general. On the other hand, POD methods require only that
random realizations of the model can be sampled: in this case, the
minimization of the projection error in POD can be straightforwardly extended
to the minimization of an empirical risk calculated using training samples of
all random variables in the model. A deep-learning/POD approach for parameter-
dependent nonlinear PDEs has been developed in [42].

Empirical balanced truncation using snapshots to approximate the Gramians
emerged soon after the first BT papers [43], and many more followed. Similarly
to [43], [44] proposes a balanced POD to compute the Gramians from numerical
simulations of state responses to unit impulses. Using similar ideas, [45]
shows that the POD kernel is an approximation to the controllability Gramian.
This is an important observation for this paper: a feedback control applies
usually to the fully observed state, in other terms the observation matrix is
the identity, which drops the main benefit of BT compared to POD: the extra
truncation of poorly observable states.

Thus, we propose a MOR method that resembles the POD to a large extent, except
that the snapshots are collected at random samples of time, parameters and
initial conditions. In contrast to particularly designed inputs (such as unit
impulses), this statistical approach makes the empirical risk a convergent
estimate of the expected risk. This method is motivated by Likelihood Informed
Subspace [46] and Active Subspace [47, 48] methods in statistics, where the
reduced basis is derived from an empirical mean Hessian of the log-likelihood,
or Gramian of the forward model. Here we propose a similar approach in the
feedback control context. Stabilization is achieved by collecting the
snapshots from a closed-loop system similarly to [35], but controlled with any
feasible (possibly suboptimal) regulator, such as the Pontryagin Maximum
Principle [49], or State-Dependent Riccati Equation [50]. The proposed method
is close to [51] which linearizes the system by estimating system matrices
from snapshots, followed by balanced truncation, and [52] where the POD is
extended to a low-rank tensor approximation of the full solution as a function
of state, parameters and time, followed by recompression of the basis at the
given parameter. However, in our statistical POD method we neither assume nor
seek any linearization, and in contrast to [52] we do this in the optimal
control setting.

Moreover, as soon as the statistical POD basis is ready, we pre-compute a low-
rank tensor approximation (specifically, the Functional Tensor Train (TT)
format [53, 54, 55]) of the closed-loop control of the reduced model as a
function of the reduced state. This TT approximation provides the desired
control in the feedback form that is fast to evaluate numerically, since it
needs only a modest amount of linear algebra operations with the TT
decomposition, in contrast to solving an optimal control problem from scratch.
This enables fast synthesis of a nearly-optimal control in the online regime.
A schematic of the entire procedure is shown in Figure 1.

Sample initial condition and parameters Solve closed-loop full systems trunc.
SVD of snapshot matrix Solve closed-loop reduced systems TT-Cross offline
Online system Interpolate TT control onlinebasisstate samplescontrol samples
project state$\begin{subarray}{c}\mbox{red.}\\\
\mbox{state}\end{subarray}$controlTT format Figure 1: Workflow of the
statistical POD and TT controller. In the offline stage, we draw (a modest
number of) random samples of the initial condition and parameters, solve the
full closed-loop systems on those parameters, and collect all snapshots into
one matrix. Computing the truncated SVD of this matrix delivers a reduced
basis minimizing the empirical risk. The TT-Cross algorithm is run to
approximate the optimal control of the reduced system in the feedback form in
the TT format. The TT-Cross draws adaptive samples of the reduced state, and
assimilates the corresponding control values. In the online stage, the system
state is measured and projected onto the reduced basis. The TT format of the
control is interpolated on this reduced state and fed back into the system.

The rest of the paper is organized as follows. In Section 2 we introduce
relevant background material on optimal control and tensor train
approximation. In Section 3 we introduce our statistical POD approach, and
present related error estimates in Section 4. The proposed methodology is
first assessed in Section 5 using Burgers’ equation, to conclude with a full
application to fluid flow control in Navier-Stokes in Section 6.

##  2 Background

###  2.1 Parameter-dependent optimal control problems

We are interested in a class of parameter-dependent optimal control problems
where a parameter $\mu\in\mathbb{R}^{M}$ accounts for uncertainties in
initial/boundary conditions of the control system propagating along state and
control trajectories. To introduce $\mu$ we assume a finite dimensional noise.
Namely, given a complete probability space $(\Omega,\mathcal{F},\mathbb{P})$
for random realisations of the system $\omega\in\Omega$, the $\sigma$-algebra
of events $\mathcal{F}$ and the probability measure $\mathbb{P}$, we assume
that any random variable $X(\omega)$ can be expressed as a deterministic
function $X(\omega)=x(\mu(\omega))$ of the random vector $\mu(\omega)$ with a
product density function $\pi(\mu)=\pi(\mu_{1})\cdots\pi(\mu_{M}),$ and
expectations with respect to $\mathbb{P}$ can be computed as integrals
weighted with $\pi$, $\mathbb{E}[X]=\int_{\mathbb{R}^{M}}x(\mu)\pi(\mu)d\mu.$
Now the system dynamics can be described as a parameter-dependent controlled
evolution equation:

$\left\\{\begin{array}[]{l}\frac{d}{dt}y(t;\mu)=f(y(t;\mu),u(t;\mu),\mu),\;\;t\in(0,T],\,\mu\in\mathbb{R}^{M},\\\
y(0;\mu)=x(\mu)\in\mathbb{R}^{d},\end{array}\right.$ (1)

where we denote by $y:[0,T]\times\mathbb{R}^{M}\rightarrow\mathbb{R}^{d}$ the
state of the system, by $u:[0,T]\times\mathbb{R}^{M}\rightarrow\mathbb{R}^{m}$
the control signal, and by
$\mathcal{U}=L^{\infty}([0,T]\times\mathbb{R}^{M};U_{ad}),U_{ad}\subseteq\mathbb{R}^{m}$,
the set of admissible controls.

Given the control system (1), we are concerned with the synthesis of a
feedback control law $u(y(t;\mu);\mu)$, that is, a control law that primarily
depends on the current state of the system $y(t;\mu)$. We design such a
control law by minimizing a finite horizon cost functional of the form

$J_{T}(u;x,\mu):=\int_{0}^{T}L(y(t;\mu),u(t;\mu))\,dt\,,\;\mu\in\mathbb{R}^{M},u\in\mathcal{U},$
(2)

where $L:\mathbb{R}^{d}\times\mathbb{R}^{m}\rightarrow\mathbb{R}$ is a
suitable running cost. Note that the dependence on the cost is with respect to
$x$, the initial condition, which generates the trajectory $y(t;\mu)$ For a
given parameter realization $\mu$ and initial condition $x(\mu)$, the optimal
control signal is given by the solution of

$\underset{u(\cdot)\in\mathcal{U}}{\inf}J_{T}(u;x,\mu)\,,\qquad\text{subject
to \eqref{eq}}\,.$ (3)

For a fixed $\mu$, we obtain deterministic dynamics, for which optimality
conditions for the dynamic optimization problem (3) are given by Pontryiagin’s
Maximum Principle (PMP)

$\displaystyle\frac{d}{dt}y(t;\mu)$ $\displaystyle=f(y(;\mu),u(t;\mu),\mu),$
(4) $\displaystyle y(0;\mu)$ $\displaystyle=x(\mu),$ (5)
$\displaystyle-\frac{d}{dt}p(t;\mu)$
$\displaystyle=\nabla_{y}f(y(t;\mu),u(t;\mu),\mu)^{\top}p(t;\mu)+\nabla_{y}L(y(t;\mu),u(t;\mu)),$
(6) $\displaystyle p(T;\mu)$ $\displaystyle=0,$ (7) $\displaystyle u(t;\mu)$
$\displaystyle=\underset{w\in
U_{ad}}{\operatorname*{arg\,min}}\\{L(y(t;\mu),w)+f(y(t;\mu),w,\mu)^{\top}p(t;\mu)\\}\,,$
(8)

where $p:[0,T]\times\mathbb{R}^{M}\rightarrow\mathbb{R}^{d}$ denotes the
adjoint variable. The PMP system can be solved using a reduced gradient method
(see $e.g.$ [6] for more details), which requires an initial guess for the
optimal control signal $u(t;\mu)$. In Section 3.1, we will explain how the
construction of reduced basis can improve the choice for the initial guess and
accelerate the entire algorithm.

###  2.2 Sub-optimal control laws using the State-Dependent Riccati Equation

The State-Dependent Riccati Equation (SDRE) is an effective technique for
feedback stabilization of nonlinear dynamics [56, 57]. The method is based on
the sequential solution of linear-quadratic control problems arising from
sequential linearization of the dynamics along a trajectory.

Given a realization of $\mu$, we consider an unconstrained
($U_{ad}=\mathbb{R}^{m}$) infinite-horizon quadratic cost functional

$J_{\infty}(u;x,\mu)=\int\limits_{0}^{+\infty}y(t;\mu)^{\top}Qy(t;\mu)+u(t)^{\top}Ru(t)\,dt\,,$
(9)

with $Q\in\mathbb{R}^{d\times d},Q\succeq 0\,,\;R\in\mathbb{R}^{m\times
m},R\succ 0\,,$ and system dynamics expressed in semilinear form

$\displaystyle\frac{d}{dt}y(t;\mu)$
$\displaystyle=A(y(t;\mu);\mu)y(t;\mu)+B(y(t;\mu);\mu)u(t)$ (10)
$\displaystyle y(0;\mu)$ $\displaystyle=x(\mu)\,.$ (11)

In the the case of linear dynamics in the state,
$A(y(t;\mu);\mu)=A(\mu)\in\mathbb{R}^{d\times d}$, and
$B(y(t;\mu);\mu)=B(\mu)\in\mathbb{R}^{d\times m}$, under standard
stabilizability assumptions, the optimal control is expressed in feedback form

$u(y;\mu)=-R^{-1}B^{\top}(\mu)P(\mu)y\,,$ (12)

where $P(\mu)\in\mathbb{R}^{d\times d}$ is the unique positive definite
solution of the Algebraic Riccati Equation (ARE)

$\displaystyle
A^{\top}(\mu)P(\mu)+P(\mu)A(\mu)-P(\mu)B(\mu)R^{-1}B^{\top}(\mu)P(\mu)+Q=0\,.$
(13)

In this case, the feedback operator $K=R^{-1}B^{\top}(\mu)P(\mu)$ does not
depend on the state. Formally, the SDRE method extends this logic to
semilinear systems by indexing with respect to the current state of the
trajectory, that is,

$u(y;\mu)=-R^{-1}B^{\top}(y;\mu)P(y;\mu)y\,,$ (14)

where $P(y;\mu)$ is the solution of the ARE

$\displaystyle
A^{\top}(y;\mu)P(y;\mu)+P(y;\mu)A(y;\mu)-P(y;\mu)B(y;\mu)R^{-1}B^{\top}(y;\mu)P(y;\mu)$
$\displaystyle=-Q$ (15)

where $A(y;\mu)$ and $B(y;\mu)$ are frozen at $y$. The implementation of the
SDRE control requires the sequential solution of eq. (15) as the state
$y(t;\mu)$ evolves in time. The SDRE feedback loop satisfies necessary
optimality conditions at a quadratic rate as the state is driven to zero (we
refer to [56] for more details and a complete statement of the result).

###  2.3 Tensor Train surrogate model for the control

Neither of the two previously presented control laws can be effectively
implemented for real-time stabilization of large-scale dynamics. Both the
numerical realization of the PMP optimality system and the computational
burden of solving AREs along a trajectory in the SDRE method surpass the time
scale required for real-time control. Hence, we propose to construct a
surrogate feedback law from state-parameter samples in a supervised learning
manner. As a surrogate ansatz we use the Tensor Train (TT) decomposition,
which can then be evaluated in real time along a controlled trajectory.

A scalar function $\tilde{u}(x):\mathbb{R}^{d}\rightarrow\mathbb{R}$ with
$x=(x_{1},\ldots,x_{d})\in\mathbb{R}^{d}$ is said to be represented in the
Functional Tensor Train (FTT) decomposition [54, 55] if it can be written as

$\tilde{u}(x)=\sum_{\alpha_{0},\ldots,\alpha_{d}=1}^{r_{0},\ldots,r_{d}}u^{(1)}_{\alpha_{0},\alpha_{1}}(x_{1})u^{(2)}_{\alpha_{1},\alpha_{2}}(x_{2})\cdots
u^{(d)}_{\alpha_{d-1},\alpha_{d}}(x_{d}),$ (16)

with some _cores_
$u^{(k)}(x_{k}):\mathbb{R}\rightarrow\mathbb{R}^{r_{k-1}\times r_{k}}$,
$k=1,\ldots,d$, and _ranks_ $r_{0},\ldots,r_{d}$. Without loss of generality,
we can let $r_{0}=r_{d}=1$, but the intermediate ranks can be larger than $1$
to account for the structure of the function. In practice, the exact
decomposition (16) may be not possible, and we seek to approximate a given
control function $u(x)$ by $\tilde{u}(x)$ in the form (16), minimising the
error $\|u-\tilde{u}\|$ in some (usually Euclidean) norm. We aim at scenarios
where $r_{k}$ stay bounded, or scale mildly. Quadratic value functions [20] or
weakly correlated Gaussian functions [58], for example, admit FTT
approximations with ranks at most polynomial in $d$ and poly-logarithmic in
the approximation error $\|u-\tilde{u}\|$. This allows one to avoid the curse
of dimensionality by replacing $u(x)$ by its approximation (16).

The FTT (16) has originally emerged as the TT decomposition [53] of tensors
(multiindex arrays). Indeed, the two come together. For practical computations
with (16) one needs to discretise each core using some basis functions
$\phi^{(k)}_{1}(x_{k}),\ldots,\phi^{(k)}_{n_{k}}(x_{k})$ in the $k$-th
variable. One can now write the $k$-th core using a 3-dimensional tensor
$\mathcal{U}^{(k)}\in\mathbb{R}^{r_{k-1}\times n_{k}\times r_{k}}$ of
expansion coefficients,

$u^{(k)}_{\alpha_{k-1},\alpha_{k}}(x_{k})=\sum_{i=1}^{n_{k}}\mathcal{U}^{(k)}(\alpha_{k-1},i,\alpha_{k})\phi_{i}^{(k)}(x_{k}).$
(17)

The multivariate function $\tilde{u}(x)$ becomes discretised in the Cartesian
basis, and defined by a $d$-dimensional tensor of expansion coefficients,

$\displaystyle\tilde{u}(x)$
$\displaystyle=\sum_{i_{1},\ldots,i_{d}=1}^{n_{1},\ldots,n_{d}}\mathcal{U}(i_{1},\ldots,i_{d})\phi^{(1)}_{i_{1}}(x_{1})\cdots\phi^{(d)}_{i_{d}}(x_{d}),$
$\displaystyle\mathcal{U}(i_{1},\ldots,i_{d})$
$\displaystyle=\sum_{\alpha_{0},\ldots,\alpha_{d}=1}^{r_{0},\ldots,r_{d}}\mathcal{U}^{(1)}(\alpha_{0},i_{1},\alpha_{1})\cdots\mathcal{U}^{(d)}(\alpha_{d-1},i_{d},\alpha_{d}).$
(18)

The tensor decomposition (18) is the original TT format [53]. Note that to
interpolate $\tilde{u}(x)$ on any given $x$ requires first $d$ univariate
interpolations of the cores as shown in (17), followed by the multiplication
of matrices $u^{(1)}(x_{1})\cdots u^{(d)}(x_{d})$. This can be implemented in
$2\sum_{k=1}^{d}n_{k}r_{k-1}r_{k}=\mathcal{O}(dnr^{2})$ operations, where we
define upper bounds $n:=\max_{k}n_{k}$ and $r:=\max_{k}r_{k}$.

To compute a TT approximation, instead of minimising the error
$\|u-\tilde{u}\|$ directly, one can use faster algorithms based on the so-
called _cross interpolation_ [59, 60, 27]. These methods consider sampling
sets of the Cartesian form $X_{<k}\oplus X_{k}\oplus X_{>k}$, for each
$k=1,\ldots,d,$ where $X_{<k}=\\{(x_{1},\ldots,x_{k-1})\\}$ of cardinality
$r_{k}$, $X_{k}=\\{x_{k}\\}$ of cardinality $n_{k}$, and
$X_{>k}=\\{(x_{k+1},\ldots,x_{d})\\}$ of cardinality $r_{k}$; $X_{<k}$ and
$X_{>k}$ are taken empty when $k=1$ and $k=d$, respectively. Note that the
number of samples in $X_{<k}\oplus X_{k}\oplus X_{>k}$ is $r_{k-1}n_{k}r_{k}$,
equal to the number of unknowns in $\mathcal{U}^{(k)}$. Therefore, one can
resolve the linear interpolation equations

$u(x)=\tilde{u}(x)\,,\quad\forall x\in X_{<k}\oplus X_{k}\oplus X_{>k}$

on the elements of $\mathcal{U}^{(k)}$ exactly, whenever $u^{(1)}(x_{1})\cdots
u^{(k-1)}(x_{k-1})$ for $(x_{1},\ldots,x_{k-1})\in X_{<k}$,
$\phi^{(k)}(x_{k})$ for $x_{k}\in X_{k}$ and $u^{(k+1)}(x_{k+1})\cdots
u^{(d)}(x_{d})$ for $(x_{k+1},\ldots,x_{d})\in X_{>k}$ are linearly
independent. Moreover, pivoting can be applied to thus computed
$\mathcal{U}^{(k)}$ to identify next sampling sets $X_{<k+1}$ and $X_{>k-1}$
as subsets of $X_{<k}\oplus X_{k}$ and $X_{k}\oplus X_{>k}$, respectively, to
make $u^{(1)}(x_{1})\cdots u^{(k)}(x_{k})$ and $u^{(k)}(x_{k})\cdots
u^{(d)}(x_{d})$ in the next step better conditioned. Iterating over
$k=1,\ldots,d$, the TT-Cross updates all TT cores until convergence. Vector
functions can be approximated component by component, since the number of
components in the control problems we consider is small. For a comprehensive
textbook on tensor methods see e.g. [61].

##  3 A statistical POD approach for optimal control problems

In this section we describe our main contribution, outlined in Figure 1. Our
ultimate goal is the fast synthesis of the feedback control that is resilient
to stochastic perturbations of the system. We split the entire procedure into
two stages. In the first (_offline_) stage, we construct a reduced state basis
that minimizes the average squared projection error for an ensemble of
controlled trajectories of random realizations of the dynamical system and/or
initial state. Due to this averaging of random samples of the error, we call
the span of thus obtained basis the Statistical POD subspace. It accommodates
states that are well controllable on average, and hence its dimension can be
much smaller than the original problem dimension. In the second (_online_)
stage, the system is projected onto the Statistical POD subspace, and an
optimal control problem is solved for the reduced system in the feedback form.
Due to its smaller dimension, the reduced optimal control problem is faster to
solve numerically, and can enable the real-time control synthesis. The
computations can be accelerated even further if the control function of the
reduced problem is approximated in the TT format. The latter step can be moved
into the offline stage, making the online stage a mere TT interpolation.

Offline Stage

The offline phase for the statistical POD method constitutes of the collection
of snapshots and the construction of the reduced basis and, hence, the reduced
dynamics. In this phase we rely on the application of the PMP system (4)-(8)
to obtain the optimal trajectory for each given sample of the parameter. We
fix $N$ different realizations
$\underline{\mu}=(\mu_{1},\ldots,\mu_{N})\in\mathbb{R}^{M\times N}$. For each
realization $\mu_{i}\in\underline{\mu}$, we solve the system (4)-(8) and we
collect the corresponding optimal trajectory
$(y^{*}(t_{1};\mu_{i}),\ldots,y^{*}(t_{n_{t}};\mu_{i}))$ at some $n_{t}$ time
instances $\\{t_{i}\\}_{i=1}^{n_{t}}$. The final snapshot matrix
$Y_{\underline{\mu}}$ will contain all the optimal trajectories at the given
time instances $t_{1},\ldots,t_{n_{t}}$ and realizations
$\mu_{1},\ldots,\mu_{N}$:

$Y_{\underline{\mu}}=[y^{*}(t_{1};\mu_{1}),\ldots,y^{*}(t_{n_{t}};\mu_{1}),y^{*}(t_{1};\mu_{2}),\ldots,y^{*}(t_{n_{t}};\mu_{N})]\in\mathbb{R}^{d\times
n_{t}N},$ (19)

where each $y^{*}(t_{i};\mu_{j})\in\mathbb{R}^{d}$ is a column vector.

At this point we perform a Singular Value Decomposition (SVD) of the snapshot
matrix
$Y_{\underline{\mu}}=U_{\underline{\mu}}\Sigma_{\underline{\mu}}V_{\underline{\mu}}^{\top}$
and we denote as $U_{\underline{\mu}}^{\ell}$ the first $\ell$ columns of the
orthogonal matrix $U_{\underline{\mu}}$. The number of basis $\ell$ can be
selected by making the truncation error

$\mathcal{E_{\underline{\mu}}}(\ell)=\frac{\sum_{i=1}^{\ell}\sigma^{2}_{i,{\underline{\mu}}}}{\sum_{i=1}^{\min\\{Nn_{t},d\\}}\sigma^{2}_{i,{\underline{\mu}}}}$
(20)

below a desired threshold. Here, $\\{\sigma_{i,{\underline{\mu}}}\\}_{i}$ are
the singular values of the matrix $Y_{\underline{\mu}}$. The quantity (20) is
related to the projection error arising in the truncation of the singular
value decomposition to the first $\ell$ basis. To shorten the notation, in
what follows, we will denote the reduced basis matrix as
$U^{\ell}:=U_{\underline{\mu}}^{\ell}$ .

The full order dynamics (1) is then reduced via projection onto the subspace
spanned by $U^{\ell}$ to the following reduced dynamics

$\left\\{\begin{array}[]{l}\frac{d}{dt}y^{\ell}(t;\mu)=(U^{\ell})^{\top}f(U^{\ell}y^{\ell}(t;\mu),u(t;\mu),\mu),\;\;t\in(0,+\infty),\mu\in\underline{\mu},\\\
y^{\ell}(0;\mu)=(U^{\ell})^{\top}x\in\mathbb{R}^{\ell}.\end{array}\right.$
(21)

To ease the notation, we define

$f^{\ell}(y^{\ell}(t;\mu),u(t;\mu),\mu)=(U^{\ell})^{\top}f(U^{\ell}y^{\ell}(t;\mu),u(t;\mu),\mu).$

The procedure for the offline stage is presented in Algorithm 1.

Algorithm 1 Offline Stage for the Statistical POD method

1:Fix the final time $T$, a time discretization $\\{t_{i}\\}_{i=1}^{n_{t}}$,
$N$ realizations $\\{\mu_{i}\\}_{i=1}^{N}$, truncation threshold
$\varepsilon\geq 0$ and $Y_{\underline{\mu}}=[\,]$

2:for $i=1,\ldots,N$ do

3: Solve the PMP system (4)-(8)

4: Compute the optimal trajectory
$Y_{\mu_{i}}=[y^{*}(t_{1};\mu_{i}),\ldots,y^{*}(t_{n_{t}};\mu_{i})]$

5: $Y_{\underline{\mu}}=[Y_{\underline{\mu}}\;Y_{\mu_{i}}]$

6:end for

7:Perform the SVD
$Y_{\underline{\mu}}=U_{\underline{\mu}}\Sigma_{\underline{\mu}}V_{\underline{\mu}}^{\top}$

8:Select $\min\ell:~{}\mathcal{E_{\underline{\mu}}}(\ell)<\varepsilon^{2}$
(20) and take the first $\ell$ columns of $U_{\underline{\mu}}$ into
$U^{\ell}.$

9:Construct the reduced dynamics (21)

######  Remark 3.1.

In general, $f(U^{\ell}y^{\ell}(t;\mu),u(t;\mu),\mu)$ is still computationally
expensive since it requires the evaluation of the nonlinearity on the lifted
variable $U^{\ell}y^{\ell}(t;\mu)\in\mathbb{R}^{d}$. For these cases one may
apply Empirical Interpolation Method (EIM, [62]) and Discrete Empirical
Interpolation Method (DEIM, [63]) to solve this issue.

######  Remark 3.2.

In some cases, such as problems arising from semidiscretization of nonlinear
PDEs, the dynamical system may present a quadratic-bilinear form:

$f(y(t;\mu),u(t;\mu),\mu)=Ay+{\bf T}(y\otimes
y)+\sum_{i=1}^{m}\theta_{k}(\mu)N_{k}y_{k}u_{k}+\theta_{0}(\mu)Bu,$ (22)

where $A,N_{k}\in\mathbb{R}^{d\times d}$, ${\bf T}\in\mathbb{R}^{d\times
d^{2}}$, $B\in\mathbb{R}^{d\times m}$ and
$\theta_{i}:\mu\rightarrow\mathbb{R}$. In this case, we can assemble the
reduced matrices

$A^{\ell}=(U^{\ell})^{\top}AU^{\ell},\;N_{k}^{\ell}=(U^{\ell})^{\top}N_{k}U^{\ell},\;{\bf
T}^{\ell}=(U^{\ell})^{\top}{\bf T}(U^{\ell}\otimes
U^{\ell}),\,B^{\ell}=(U^{\ell})^{\top}B$

and define the following reduced dynamics

$\frac{d}{dt}y^{\ell}(t;\mu)=A^{\ell}y^{\ell}+{\bf T}^{\ell}(y^{\ell}\otimes
y^{\ell})+\sum_{i=1}^{m}\theta_{k}(\mu)N^{\ell}_{k}y^{\ell}_{k}u_{k}+\theta_{0}(\mu)B^{\ell}u,$
(23)

avoiding the use of empirical interpolation techniques.

Online Stage

Given the reduced dynamics (21), the reduced PMP for a given $\mu$ reads

$\displaystyle\frac{d}{dt}y^{\ell}(t,\mu)$
$\displaystyle=f^{\ell}(y^{\ell}(t;\mu),u(t;\mu),\mu),$ (24) $\displaystyle
y^{\ell}(0;\mu)$ $\displaystyle=y^{\ell}_{0}(\mu),$ (25)
$\displaystyle-\frac{d}{dt}p^{\ell}(t;\mu)$
$\displaystyle=\nabla_{y^{\ell}}f^{\ell}(y^{\ell}(t;\mu),u(t;\mu),\mu)^{\top}p^{\ell}(t;\mu)+\nabla_{y^{\ell}}L(U^{\ell}y^{\ell}(t;\mu),u(t;\mu)),$
(26) $\displaystyle p^{\ell}(T;\mu)$ $\displaystyle=0,$ (27) $\displaystyle
u(t)$ $\displaystyle=\operatorname*{arg\,min}_{w\in
U_{ad}}\\{L(U^{\ell}y^{\ell}(t;\mu),w)+f^{\ell}(y^{\ell}(t;\mu),w,\mu)^{\top}p^{\ell}(t;\mu)\\}\,.$
(28)

The reduced-order state variable is $\ell$-dimensional and it can be exploited
to construct a faster approximation of the optimal control. The reduced-order
dynamics enables the synthesis of an optimal feedback law, now set in a lower
dimensional domain, helping in the mitigation of the curse of dimensionality.

Let us consider the feedback map $u^{\ell}:\mathbb{R}^{\ell}\rightarrow
U_{ad}$ obtained solving the reduced PMP (24)-(28). This controller can be
used in the full order model dynamics considering first the projection of the
state and then the evaluation of the feedback map:

$\dot{y}(t;\mu)=f(y(t;\mu),u^{\ell}((U^{\ell})^{\top}y(t;\mu)),\mu),$ (29)

obtaining a control problem where the computation of the optimal feedback is
independent from the original dimension of the dynamical system.

######  Remark 3.3.

The presented methodology is based on the construction of a snapshot set of
optimal trajectories obtained via the resolution of the PMP system (4)-(8).
The sampling of the controlled dataset may be derived also upon the SDRE
framework introduced in Section 2.2, solving the sequential AREs (15) along
the optimal trajectory. Since the number of variables in the dynamical system
may be arbitrary large, the SDRE approach requires efficient solvers for the
resolution of high-dimensional Riccati equations. We refer to [2] for a
comprehensive discussion of such techniques.

######  Remark 3.4.

The singular values for the Statistical POD method may present a slow decay
due to the nature of the problem and the size of the parameter domain. In this
context a good approximation requires a large number of POD basis and the
online stage will still be affected by the dimensionality problem. To this
end, an efficient surrogate model for the feedback control in the reduced
space accelerates the procedure, obtaining a fast and reliable method for the
computation of the optimal control solution. We consider the representation of
the feedback control $u^{\ell}(x^{\ell})$ in the Functional Tensor Train
format (16), obtained via a TT Cross procedure (see Section 2.3). The optimal
reduced trajectories obtained by solving the system (24)-(28) are employed as
sampling points for the Cross approximation.

######  Remark 3.5.

In contrast to existing POD techniques for the approximation of optimal
control problems, the proposed Statistical POD method introduces stochastic
terms in the dynamical system to explore the manifold of controlled solutions
for the parameterized problem. This enforces the robustness of the
corresponding reduced basis to perturbations, an essential feature for
feedback control. In particular, as shown in Section 4, the empirical risk of
the Statistical POD procedure converges to the expected risk. Thus, increasing
the amount of training data for constructing the reduced-order system, we
obtain a model that is more representative of the true system dynamics,
including unseen scenarios.

###  3.1 Efficient computation of the snapshots

The main building block in the offline phase is represented by the computation
of different optimal trajectories via the PMP system (4)-(8). The resolution
of this system relies strongly on the choice of the initial guess for the
optimal control $u^{*}(\cdot)$. In this section we introduce a procedure which
accelerates the high-dimensional computation exploiting the information about
the optimum for the reduced problem. Let us consider a prefixed initial guess
$u^{0}(\cdot)$. In the offline stage presented in Algorithm 1, the different
optimal trajectory realizations are stored in the snapshot matrix
$Y_{\underline{\mu}}$ and the basis are constructed only at the end of the
loop. An SVD decomposition can be added between Lines 2-6 of Alg. 1 to
construct a reduced basis containing partial information on the statistical
behaviour of the optimal trajectories and build a reduced dynamics. The
constructed reduced system is employed for a fast computation of the
approximated optimal control $u_{red}(\cdot)$, which can be used as initial
guess for the high-dimensional computation. In particular, at the $i$-th step
we consider as initial guess the control which achieves a lower cost
functional, $i.e.$

$u^{i}(\cdot)\in\operatorname*{arg\,min}_{u(\cdot)\in\\{u^{0}(\cdot),u_{red}^{i}(\cdot)\\}}J(u).$
(30)

The optimized procedure is sketched in Algorithm 2.

Algorithm 2 Optimized Offline Stage for the Statistical POD method

1:Fix the final time $T$, a time discretization $\\{t_{i}\\}_{i=1}^{n_{t}}$,
$N$ realizations $\\{\mu_{i}\\}_{i=1}^{N}$, truncation threshold
$\varepsilon\geq 0$, an initial control guess $u^{0}$ and
$Y_{\underline{\mu}}=[\,]$

2:for $i=1,\ldots,N$ do

3: if $i==1$ then

4: Solve the PMP system (4)-(8) with initial guess $u^{0}$

5: else

6: Solve the reduced PMP system (24)-(28) and compute $u_{red}^{i}$

7: Solve the PMP system (4)-(8) with initial guess given by (30)

8: end if

9: Compute the optimal trajectory
$Y_{\mu_{i}}=[y^{*}(t_{1};\mu_{i}),\ldots,y^{*}(t_{n_{t}};\mu_{i})]$

10: $Y_{\underline{\mu}}=[Y_{\underline{\mu}}\;Y_{\mu_{i}}]$

11: Perform the SVD
$Y_{\underline{\mu}}=U_{\underline{\mu}}\Sigma_{\underline{\mu}}V_{\underline{\mu}}^{\top}$

12: Select $\min\ell:~{}\mathcal{E_{\underline{\mu}}}(\ell)<\varepsilon^{2}$
and take the first $\ell$ columns of $U_{\underline{\mu}}$ into $U^{\ell}.$

13: Construct the reduced dynamics (21)

14:end for

##  4 Expected risk of statistical POD

We analyse the proposed technique providing a rigorous error estimate on the
expected risk. For this reason we first recall an eigenvalue perturbation
theory, the Bauer-Fike theorem, which will be employed for the proof of
Proposition 1.

######  Theorem 1 (Bauer-Fike).

Let $A$ be an $d\times d$ diagonalizable such that $A=V\Lambda V^{-1}$ and let
$E$ be an arbitrary $d\times d$ matrix. Then for every $\mu\in\sigma(A+E)$
there exists $\lambda\in\sigma(A)$ such that

$|\lambda-\mu|\leq\kappa_{p}(V)\|E\|_{p}$

where $\kappa_{p}(V)$ is the condition number in $p$-norm of the matrix $V$.

We are ready to formulate the error estimate on the expected projection error
characterizing the statistical POD approach.

######  Proposition 1.

Consider $t$ as a uniformly distributed random variable on $[0,T]$, and
consider $N$ iid realizations of the optimal trajectories
$y^{(i)}:=y(t_{i};\mu_{i})$, $i=1,\ldots,N$. Assume that
$\sqrt{\mathrm{Var}[y_{j}y_{k}]}\leq C<\infty$ for any $j,k=1,\ldots,d$, where
$y_{j}$ is the $j$th component of the random vector $y(t,\mu)$. Let

$G=\frac{1}{N}\sum_{i=1}^{N}y(t_{i};\mu_{i})y(t_{i};\mu_{i})^{\top}\quad\mbox{and}\quad
G^{*}=\mathbb{E}[yy^{\top}],$

and let $U^{\ell}\in\mathbb{R}^{d\times\ell}$ be an orthonormal matrix of
$\ell$ leading eigenvectors of $G$. Then

$\mathbb{E}[\|y(t;\mu)-U^{\ell}(U^{\ell})^{\top}y(t;\mu)\|_{2}^{2}]\leq\sum_{j=\ell+1}^{d}\lambda^{*}_{j}+2d\ell\frac{C}{\sqrt{N}}\leq\sum_{j=\ell+1}^{d}\lambda_{j}+(d+\ell)d\frac{C}{\sqrt{N}},$
(31)

where $\lambda^{*}_{j}$ are eigenvalues of $G^{*}$ sorted descending, and
$\lambda_{j}$ are eigenvalues of $G$ sorted descending.

###### Proof.

Let $G^{*}=U_{*}\Lambda_{*}U_{*}^{\top}$ and $G=U\Lambda U^{\top}$ be the
eigenvalue decompositions of $G^{*}$ and $G$, respectively, where eigenvalues
in $\Lambda$ and $\Lambda_{*}$ are sorted descending. Note that the matrix
$U^{\ell}$ contains the first $\ell$ columns of $U$. Moreover, let us
introduce $P^{\ell}=U^{\ell}(U^{\ell})^{\top}$ to shorten the notation.

Since $y(t;\mu)$ in (31) is independent of precomputed trajectories $y^{(i)}$,
and hence of $G$ and $U^{\ell}$, we can factorise the total expectation into
those over $y$ and $G$,

$\mathbb{E}[\|y-U^{\ell}(U^{\ell})^{\top}y\|_{2}^{2}]=\mathbb{E}_{G}[\mathbb{E}_{y}[\|y-U^{\ell}(U^{\ell})^{\top}y\|_{2}^{2}]],$

and hence express first

$\mathbb{E}_{y}[\|y-P^{\ell}y\|_{2}^{2}]=\mathbb{E}_{y}[y^{\top}y-y^{\top}P^{\ell}y]=\mathrm{tr}(\mathbb{E}_{y}[yy^{\top}-P^{\ell}yy^{\top}])$

due to orthogonality and the cyclic permutation under the trace of a matrix.
Now,

$\mathrm{tr}(\mathbb{E}_{y}[yy^{\top}])=\mathrm{tr}(\Lambda_{*})=\sum_{j=1}^{d}\lambda_{j}^{*},$

while

$\displaystyle\mathrm{tr}(\mathbb{E}_{y}[P^{\ell}yy^{\top}])$
$\displaystyle=\mathrm{tr}(P^{\ell}G^{*})=\mathrm{tr}(P^{\ell}G)+\mathrm{tr}(P^{\ell}(G^{*}-G))$
$\displaystyle=\mathrm{tr}((U^{\ell})^{\top}U\Lambda
U^{\top}U^{\ell}))+\mathrm{tr}(P^{\ell}(G^{*}-G)).$

Due to orthogonality, the first trace is just the trace of the leading
$\ell\times\ell$ submatrix of $\Lambda$, which is
$\sum_{j=1}^{\ell}\lambda_{j}$. Now taking also $\mathbb{E}_{G}$, we obtain

$\displaystyle\mathbb{E}[\|y-P^{\ell}y\|_{2}^{2}]$
$\displaystyle=\sum_{j=1}^{d}\lambda_{j}^{*}-\mathbb{E}_{G}[\sum_{j=1}^{\ell}\lambda_{j}]-\mathbb{E}_{G}[\mathrm{tr}((U^{\ell})^{\top}(G^{*}-G)U^{\ell})]$
$\displaystyle=\sum_{j=\ell+1}^{d}\lambda_{j}^{*}+\mathbb{E}_{G}[\sum_{j=1}^{\ell}(\lambda_{j}^{*}-\lambda_{j})]-\mathbb{E}_{G}[\mathrm{tr}((U^{\ell})^{\top}(G^{*}-G)U^{\ell})]$
$\displaystyle\leq\sum_{j=\ell+1}^{d}\lambda_{j}^{*}+\ell\cdot\mathbb{E}_{G}[\|G^{*}-G\|_{2}]+\ell\cdot\mathbb{E}_{G}[\|G^{*}-G\|_{2}],$

due to the Bauer-Fike Theorem 1 (second term), and

$|\mathrm{tr}(Q^{\top}AQ)|=|\sum_{j=1}^{\ell}q_{j}^{\top}Aq_{j}|\leq\ell\max_{j\in\\{1,\ldots,\ell\\}}|q_{j}^{\top}Aq_{j}|\leq\ell\|A\|_{2}$

for any $A\in\mathbb{R}^{d\times d}$ and orthonormal
$Q\in\mathbb{R}^{d\times\ell}$ for the third term.

By the Jensen’s inequality and norm equivalence,

$(\mathbb{E}_{G}[\|G^{*}-G\|_{2}])^{2}\leq\mathbb{E}_{G}[\|G^{*}-G\|_{2}^{2}]\leq\mathbb{E}_{G}[\|G^{*}-G\|_{F}^{2}]=\sum_{j,k=1}^{d}\mathbb{E}_{G}[(G^{*}_{j,k}-G_{j,k})^{2}].$

Since $\mathbb{E}_{G}[G]=G^{*}$,
$\mathbb{E}_{G}[(G^{*}_{j,k}-G_{j,k})^{2}]=\mathrm{Var}[G_{j,k}].$ In turn,
$G_{j,k}$ is a sum of iid random variables $y^{(i)}_{j}y^{(i)}_{k}/N$. For
those,

$\mathrm{Var}[G_{j,k}]=N\cdot\mathrm{Var}[y^{(i)}_{j}y^{(i)}_{k}/N]=\frac{1}{N}\mathrm{Var}[y_{j}y_{k}]\leq\frac{C^{2}}{N},$

where in the penultimate step we used that $y^{(i)}$ are samples from the same
distribution of $y$. Summing over $j,k$ and taking the square root gives the
first claim of the proposition. The second claim comes from rewriting

$\sum_{j=\ell+1}^{d}\lambda^{*}_{j}=\sum_{j=\ell+1}^{d}\lambda_{j}+\sum_{j=\ell+1}^{d}(\lambda^{*}_{j}-\lambda_{j}),$

and using again Theorem 1 to get
$|\lambda^{*}_{j}-\lambda_{j}|\leq\|G^{*}-G\|_{2}$. ∎

Note that the expected risk consists of the truncated eigenvalues (bias), as
for the classical POD technique, and a variance term decaying as $1/\sqrt{N},$
as usual for Monte-Carlo methods.

##  5 Neumann boundary control for 2D Burgers’ equation

To better illustrate the different blocks of our approach, in this section we
study a stabilization problem for the solution of the viscous Burgers’
equation in a bidimensional domain, where the control appears as a Neumann
boundary condition.

###  5.1 Problem formulation

More precisely, we consider the following state equation

$\begin{cases}\partial_{t}y-\nu\Delta y+y\cdot\nabla y=\alpha
y&(\xi,t,\mu)\in\Xi\times[0,T]\times\mathbb{R}^{M},\\\
\nu\partial_{n}y=u(t)&(\xi,t,\mu)\in\partial\Xi_{1}\times[0,T]\times\mathbb{R}^{M},\\\
\nu\partial_{n}y=0&(\xi,t,\mu)\in\partial\Xi_{2}\times[0,T]\times\mathbb{R}^{M},\\\
y(\xi,0;\mu)=\tilde{y}_{0}(\xi;\mu)&(\xi,\mu)\in\Xi\times\mathbb{R}^{M},\\\
\end{cases}$ (32)

where $\nu>0$ is the viscosity constant, $\alpha\geq 0$ regulates the unstable
term, $\xi=(\xi_{1},\xi_{2})\in\Xi=[0,1]^{2}$ is the spatial variable,
$\partial\Xi_{1}=\\{0\\}\times[0,1]$ and
$\partial\Xi_{2}=\partial\Xi\setminus\partial\Xi_{1}$, see Figure 2. Our aim
is to drive the dynamical system to the equilibrium $\overline{y}\equiv 0$ and
the corresponding cost functional we want to minimize reads:

$J_{T}(u;\mu)=\int_{0}^{T}\int_{\Xi}|y(\xi,t;\mu)|^{2}\,d\xi\,dt+\int_{0}^{T}|u(t)|^{2}\;dt.$
(33)

$\partial\Xi_{1}$$\partial\Xi_{2}$$\partial\Xi_{2}$$\partial\Xi_{2}$

Figure 2: Domain for the Burgers’ equation (32).

For the application of the statistical POD technique, we assume that the
initial condition is a realisation of the following random field,

$\tilde{y}_{0}(\xi;\mu)=y_{0}(\xi)+\sum_{i=1}^{M_{1}}\sum_{j=1}^{M_{2}}(i+j)^{-\gamma}\mu_{i+(j-1)M_{1}}\cos(i\pi\xi_{1})\cos(j\pi\xi_{2}),$
(34)

where $\mu_{i+(j-1)M_{1}}\sim\mathcal{N}(0,\sigma^{2})$ are normally
distributed random variables with zero mean and variance $\sigma^{2}$, $M_{1}$
and $M_{2}$ are the maximal frequencies in the first and second spatial
variables, and $y_{0}(\xi)$ is the mean initial state. We take
$\sigma=\sigma_{0}=0.05$ and $y_{0}\equiv 0.05$ by default. The parameter
$\gamma$ refers to the decay of the Fourier coefficients and it is related to
the regularity we want to assume. In Figure 3 we can observe two initial
conditions fixing $\gamma=3$ on the left and $\gamma=4$ on the right. Higher
values for $\gamma$ correspond to smoother initial conditions since the
Fourier coefficients decay more rapidly.

Figure 3: Observed initial conditions for $\gamma=3$ (left) and $\gamma=4$
(right).

First, we discretize the Burgers’ equation (32) via $P^{1}$ finite elements
$\\{\phi_{i}\\}_{i=1}^{d}$ in space centered at uniform grid points
$\\{(\xi^{i}_{1},\xi^{i}_{2})\\}_{i=1}^{d}$, obtaining the following
semidiscrete form

$E\dot{y}(t)=(C+\alpha E)y(t)+F(y(t)\otimes y(t))+Bu(t),$ (35)

where $E\in\mathbb{R}^{d\times d}$ is the mass matrix,
$C\in\mathbb{R}^{d\times d}$ is the stiffness matrix, $F\in\mathbb{R}^{d\times
d^{2}}$ arises from the nonlinear convective term and

$(B)_{i}=\begin{cases}1&(\xi^{i}_{1},\xi^{i}_{2})\in\partial\Xi_{1},\\\
0&\mbox{otherwise,}\\\ \end{cases}\quad i=1,\ldots,d.$

We apply an implicit Euler scheme for the time discretization. In particular,
we consider $d=289$ finite elements and $n_{t}=100$ time steps, fixing the
final time $T=50$. The computation of the optimal trajectories is based on the
Pontryagin’s Maximum Principle system (4)-(8), solved via a gradient-descent
method for each sample of $\mu$. More precisely, we consider the state-adjoint
system and the gradient of the cost functional following [64].

Once the reduced basis are computed, following the construction of the reduced
dynamics (23), the reduced system reads

$\displaystyle\frac{d}{dt}{y^{\ell}}(t;\mu)$
$\displaystyle=A^{\ell}y^{\ell}(t;\mu)+F^{\ell}(y^{\ell}\otimes
y^{\ell})+B^{\ell}u(t),$ (36) $\displaystyle y^{\ell}(0;\mu)$
$\displaystyle=x^{\ell}=(U^{\ell})^{\top}x,$

which can be rewritten in the following semilinear form

$\dot{y^{\ell}}(t;\mu)=\mathcal{A}^{\ell}(y^{\ell})y^{\ell}(t;\mu)+B^{\ell}u(t),$
(37)

where

$(\mathcal{A}^{\ell}(y))(i,j)=A^{\ell}(i,j)+\sum_{k=1}^{\ell}F^{\ell}(i,(j-1)\ell+k)\,y_{k},\quad
i,j\in\\{1,\ldots,\ell\\}.$

Note that the reduced system (36) is free from $\mu$ (which appears only the
full initial state $\tilde{y}_{0}$), hence (36) can be written in the usual
feedback form where the reduced initial state is a variable $x^{\ell}$. For
the remaining part of the section we fix default $\gamma=4$, $M_{1}=M_{2}=8$,
$\nu=0.02$ and $T=50$.

###  5.2 Error indicators.

Throughout the numerical experiments, we will measure the error of the reduced
model using the following error indicators:

$\mathcal{E}_{J}=|J(u^{*})-J(u_{red}^{*})|,\quad\mathcal{E}_{y}=\sum_{k=1}^{n_{t}}(t_{k}-t_{k-1})\|y(u^{*},t_{k})-y(u_{red}^{*},t_{k})\|_{2},$
(38)

where $J(u)$ is the total cost of the full model computed on the control $u$,
$u^{*}$ is the optimal control of the full model, $u_{red}^{*}$ is the optimal
control of the reduced model, while $y(u,t)$ is the full system state with
control $u$ at time $t$. In other terms, $\mathcal{E}_{J}$ and
$\mathcal{E}_{y}$ test respectively the cost and the trajectory in the
original model when the control signal is computed using the reduced model.
The TT performance will be assessed using

$\mathcal{E}_{TT}=|J(u_{red,TT})-J(u_{red})|,$ (39)

where $u_{red}$ is any feasible feedback law for the reduced system (not
necessarily optimal), and $u_{red,TT}$ is its TT approximation. For
simplicity, in the numerical tests we use the SDRE controller.

###  5.3 SPOD accuracy in the stable case ($\alpha=0$)

We begin considering the stable case $\alpha=0$, avoiding the exponential
growth for the solution, while in the second subsection we analyse the
behaviour for the unstable case. We apply the statistical POD strategy for the
reduction of (32)-(33) with $\alpha=0$ and we compare it with the POD methods
based on the information of two specific dynamics: the uncontrolled solution
and the optimal trajectory starting from the reference initial condition
$y_{0}$. The training for the statistical POD is based on $N=40$ independent
optimal trajectories each started from a sample of the initial condition (34).

In the left panel of Figure 5 we show the decay of singular values of
$Y_{\underline{\mu}}$ varying the number of samples $N\in\\{1,20,40\\}$. As
expected, inclusion of more snapshots from different model regimes produces a
slower decay.

However, singular vectors of $Y_{\underline{\mu}}$ with more snapshots produce
a more accurate reduced model. In the right panel of Figure 4 we compare the
model reduction errors (38) for three different sets of snapshots used for
POD: $Y_{\underline{\mu}}$ consisting of $N=40$ controlled trajectories
started from random samples of the initial condition (”statistical POD”), one
controlled trajectory started from the mean initial condition (”1 POD
controlled”), and one uncontrolled trajectory started from the mean initial
condition (”1 POD uncontrolled”). In all three cases $\ell=20$ dominant
singular vectors of the snapshot matrix are selected as the basis for model
reduction. We see that the uncontrolled snapshots give a very inaccurate basis
for reduction of the controlled system, using the controlled snapshots gives a
more accurate model, especially when multiple random samples of snapshots are
used.

Figure 4: Left: singular values of $Y_{\underline{\mu}}$ for different numbers
$N$ of realisations of initial conditions. Right: mean error in the cost
$\mathcal{E}_{J}$ depending on the reduced basis size $\ell$ for statistical
POD built upon $N=40$ controlled trajectory realisations started from random
iid samples from (34), and the POD built upon 1 realisation of controlled and
uncontrolled trajectory started from the mean $y_{0}$.

In Figure 5 we compare the cumulative CPU time for the computation of the
offline phase with Algorithm 1 and with its optimized version, Algorithm 2.
The ”not-optimized” algorithm has an almost constant increase in the CPU time
and at the first iterations it performs slightly better since the statistical
basis constructed with the first snapshots do not help in the construction of
a good initial guess. On the other hand, the increase of the knowledge on the
ensemble of reduced dynamics leads to a decrease of the CPU time at each step
for the optimized version, demonstrating its beneficial support for the high-
dimensional resolution of the PMP system.

Figure 5: Comparison of the CPU time for the offline phase and the optimized
offline phase increasing the number of samples.

In Table 1 we compare the accuracy of the three techniques changing the
parameters of the problem. In the first case we double the standard deviation
for the random variables, such that
$\mu_{i}\sim\mathcal{N}(0,4\sigma_{0}^{2})$, in the second case we consider
$\gamma=3$, obtaining a set of initial conditions less smooth. The dimension
of the reduced basis in each case is $\ell=6$. The uncontrolled dynamics is
not sufficient to build a good reduced model.

The POD technique based on controlled solutions obtains better results for
both indicators, especially for the statistical POD. In particular, fixing
$\gamma=3$ the statistical POD reveals an improvement of one order of
magnitude.

Error | Uncontrolled | 1 POD | SPOD
---|---|---|---
| $\sigma=2\sigma_{0}$ | $\sigma=\sigma_{0}$ | $\sigma=2\sigma_{0}$ | $\sigma=\sigma_{0}$ | $\sigma=2\sigma_{0}$ | $\sigma=\sigma_{0}$
| $\gamma=4$ | $\gamma=3$ | $\gamma=4$ | $\gamma=3$ | $\gamma=4$ | $\gamma=3$
$\mathcal{E}_{y}$ | 0.052 | 4.3e-3 | 1.11e-3 | 1.1e-3 | 3.17e-4 | 2.0e-4
$\mathcal{E}_{J}$ | 0.047 | 5.2e-3 | 1.65e-3 | 2.3e-4 | 1.58e-4 | 1.1e-4
Table 1: Mean errors in the cost functional and optimal trajectory started
from random iid samples from (34) with modified standard deviation or
$\gamma$. The reduced basis size is $\ell=6$ for all methods. Statistical POD
(SPOD) is computed from $N=20$ trajectories started from samples from the
original (34) with $\sigma=\sigma_{0}$ and $\gamma=4$.

In Table 2 we compare the 1 POD technique and the statistical POD with
completely different initial data which cannot be represented in the form
(34). For these cases we fix the reduced basis size $\ell=30$ to obtain
accurate approximations. We note that for all the test cases the statistical
approach is more accurate than 1 POD, with a difference of almost one order of
magnitude. Furthermore, we see that the accuracy is decreasing for the
different choices of $y_{0}$. This is due to the fact that the distance
between the initial conditions considered in Table 2 and the mean of the ones
represented in the form (34) is increasing, leading to the necessity of
increasing the basis size to obtain the same accuracy.

| SPOD | 1 POD
---|---|---
$y_{0}$ | $\mathcal{E}_{y}$ | $\mathcal{E}_{J}$ | $\mathcal{E}_{y}$ | $\mathcal{E}_{J}$
0.1 | 7.66e-5 | 3.14e-6 | 1.34e-4 | 3.11e-5
0.2 | 6.07e-5 | 1.81e-5 | 1.88e-4 | 1.52e-4
$\cos(\pi\xi_{1})\cos(\pi\xi_{2})$ | 4.72e-3 | 9.48e-3 | 1.87e-2 | 0.20
Table 2: Errors in the cost and controlled state started from different
initial conditions for the 1 POD and SPOD (with $N=20$) methods, both building
reduced bases of the same size $\ell=30$.

This demonstrates the higher accuracy of the statistical approach even for
initial conditions far from the set of training samples.

###  5.4 Faster control synthesis using a TT approximation

We now assess the pre-computation of the TT approximation of the control in
the feedback form, $u_{red,TT}(x^{\ell})\approx u_{red}(x^{\ell})$, aiming at
faster synthesis of $u_{red,TT}(x^{\ell})$ in the online regime by simply
interpolating the TT decomposition, in contrast to computing
$u_{red}(x^{\ell})$ via optimisation. Recall that as long as $x^{\ell}$ is
fixed, the rest of the model is independent of $\mu$, therefore,
$u_{red,TT}(x^{\ell};\mu)$ is actually just $u_{red,TT}(x^{\ell})$. We fix the
reduced model order $\ell=20$ for all methods in this subsection. Recall that
$u_{red,TT}(x^{\ell})$ is computed in the offline stage via the TT-Cross
algorithm sampling $u_{red}(x^{\ell})$ at certain states
$x^{\ell}\in[-1,1]^{20}$. We fix the stopping $tol=10^{-3}$ for the TT-Cross,
and $n=6$ Legendre basis functions for discretizing each component of
$x^{\ell}$. For faster computations in this offline stage, we compute
$u_{red}(x^{\ell})$ via the State-Dependent Riccati Equation (SDRE) applied to
the semilinear system (37) instead of PMP. This introduces a negligible
difference to the total cost, $e.g.$ choosing
$y_{0}=0.5\cos(\pi\xi_{1})\cos(\pi\xi_{2})$ the error between SDRE and PMP in
the computation of the total cost is order $10^{-3}$.

In Table 3 we notice that the SPOD method also gives a more structured reduced
model which lends itself to a more accurate TT approximation of the optimal
control.

$y_{0}$ | SPOD | 1 POD
---|---|---
0.05 | 2.47e-6 | 1.90e-3
0.1 | 7.69e-5 | 5.42e-3
$0.5\cos(\pi\xi_{1})\cos(\pi\xi_{2})$ | 6.64e-6 | 9.21e-3
Table 3: Average errors (39) due to the TT approximation for the 1 POD and
SPOD (computed from $N=20$ trajectories) methods.

Now we compare the SDRE controller (with and without TT approximation) and the
LQR controller. The advantage of the latter is that only one Riccati equation
(13) (that for the linearization at the origin) needs to be solved and used in
the computation of the feedback law (14), but we show that the absence of
nonlinear terms in the control strategy leads to a slower stabilization of the
system. In the left panel of Figure 6 we show the shape of the initial
condition $y_{0}=0.5\cos(\pi x)\cos(\pi y)$, and in the right panel we show
the running cost of the solution with the three controllers. We see that the
LQR controller gives a higher cost, which makes it less attractive despite the
computational simplicity. The reduction of the computing time can be achieved
using the TT interpolation instead.

In the same figure we see that the costs using the direct SDRE computation and
the TT interpolation are indistinguishable, indicating a negligible error due
to the TT approximation.

Figure 6: Initial condition $y_{0}=0.5\cos(\pi\xi_{1})\cos(\pi\xi_{2})$ (left)
and decay of the running cost with different techniques starting from this
$y_{0}$ (right).

CPU times of computing the control for one system state are compared in Table
4. The dimension of the reduced space is still fixed equal to $20$. First of
all, we note that the application of SPOD enables a speed-up of 50 times for
Pontryagin’s, while the use of SDRE achieves a speed up of two orders with
respect to the reduced PMP. Finally, the computation of a TT surrogate
function gains a further extra speed-up order of magnitude, achieving in the
end a final acceleration of 5 orders between the reduced TT-SDRE and the full
PMP.

Full PMP | Reduced PMP | Reduced SDRE | Reduced TT
---|---|---|---
41.22 s | 0.78 s | 4.25e-3 s | 4.57e-4 s
Table 4: Comparison in terms of the CPU time for the computation of a single
control with different FOM and ROM solver, fixing $\ell=20$ and
$y_{0}(\xi)=0.5\cos(\pi\xi_{1})\cos(\pi\xi_{2})$.

###  5.5 SPOD accuracy in the unstable case ($\alpha>0$)

In this subsection we move to the unstable case, fixing the parameter
$\alpha=0.2$ in the dynamical system (32). The addition of this term changes
drastically the behaviour of the uncontrolled dynamics, characterized now by
an exponential growth. It is clear that in this case the construction of the
POD basis based on the uncontrolled solution is meaningless, since it quickly
diverges far away from the desired state. In the left panel of Figure 7 we
show the decay of the singular values of different POD strategies varying the
number of sampled initial conditions. The decay of the singular values

is slower than that for the stable dynamics, reflecting the more complex
nature of the problem. In the right panel of Figure 7 we report the mean error
in the computation of the cost functional for 1 POD and the statistical
approach varying the basis size. The initial conditions are sampled again in
the form (34). The POD method based on the uncontrolled dynamics is not
reported since it achieves an error of order $\approx 10^{5}$, demonstrating
its inefficiency in this context. SPOD presents a better accuracy for all the
test cases.

Figure 7: Left: singular values of $Y_{\underline{\mu}}$ for different numbers
of offline realisations $N$. Right: mean error in the cost $\mathcal{E}_{J}$
(38) depending on the reduced basis size $\ell$ for SPOD built upon $N=40$
controlled trajectory realisations started from random iid samples from (34),
and the POD built upon 1 realisation of the controlled trajectory started from
the mean $y_{0}$.

In Table 5 we report the comparison of the two techniques varying the
parameters appearing in the definition of the initial conditions (34) with
$\ell=20$ reduced basis. Again, in the first columns we consider as random
variable $\mu_{i}\sim\mathcal{N}(0,4\sigma_{0}^{2})$, doubling the reference
standard deviation considered in (34), while in the last columns we fix the
decay exponent $\gamma=3$. The SPOD gets a better results for all the
indicators, especially for the case $\sigma=2\sigma_{0}$, where
$\mathcal{E}_{y}$ performs two order of magnitudes better than the standard
POD.

| 1 POD | SPOD
---|---|---
Error | $\sigma=2\sigma_{0}$ | $\sigma=\sigma_{0}$ | $\sigma=2\sigma_{0}$ | $\sigma=\sigma_{0}$
| $\gamma=4$ | $\gamma=3$ | $\gamma=4$ | $\gamma=3$
$\mathcal{E}_{y}$ | 3.5e-4 | 3.91e-4 | 2.1e-6 | 1.33e-4
$\mathcal{E}_{J}$ | 7.3e-5 | 7.57e-5 | 1.8e-6 | 9.15e-6
Table 5: Unstable Burgers’ example: mean errors in the cost and optimal
trajectory of the reduced model of order $\ell=20$ started from random iid
samples from (34) varying the standard deviation $\sigma$ or the decay
parameter $\gamma$. The SPOD basis is built from $N=20$ trajectories started
from samples from the original (34) with $\sigma=\sigma_{0}$ and $\gamma=4$.

Finally, we pass to the application of the SDRE considering the reduced
dynamical system (37) and we apply the Tensor Train Cross for the construction
of a surrogate model. We fix a number of basis $\ell=20$. In Table 6 we
compare the total cost obtained using LQR, SDRE and its approximation via TT.
We immediately notice that SDRE achieves better results than LQR for all the
study tests and the TT approximation obtains almost the same cost as SDRE for
the digits displayed.

$y_{0}$ | LQR | SDRE | TT-SDRE
---|---|---|---
0.05 | 0.0809 | 0.0808 | 0.0808
Random i.c. (34) | 0.0798 | 0.0756 | 0.0756
$0.5\cos(\pi\xi_{1})\cos(\pi\xi_{2})$ | 0.1592 | 0.0879 | 0.0878
Table 6: Unstable Burgers’ example: total cost $J$ with different initial
conditions for LQR, SDRE and TT-SDRE with SPOD (with $N=20$ samples) methods.

The left panel of Figure 8 shows the configuration of the controlled solution
via TT-SDRE at final time with infinity norm of order $10^{-2}$. In the right
panel of Figure 8 the running costs for LQR, SDRE and TT SDRE. There is a
visual superposition of the curves for SDRE and TT SDRE, reflected also in
Table 6. The SDRE controller is able to reduce the running cost faster, which
is reflected by the lower total cost in Table 6.

Figure 8: Unstable Burgers’ example: final configuration for the controlled
solution via TT-SDRE (left) and decay of the running cost with different
techniques (right) starting from $y_{0}=0.5\cos(\pi\xi_{1})\cos(\pi\xi_{2})$.

##  6 Dirichlet boundary control for the Navier-Stokes equations

We consider a more challenging problem: the optimal control of the 2D
incompressible Navier-Stokes (NS) equation via Dirichlet boundary control.

###  6.1 Problem formulation

The NS equation reads

$\begin{cases}\partial_{t}y-\nu\Delta y+y\cdot\nabla y+\nabla
p=0&(\xi,t,\mu)\in\Xi\times[0,T]\times\mathbb{R}^{M},\\\ \nabla\cdot
y=0&(\xi,t,\mu)\in\Xi\times[0,T]\times\mathbb{R}^{M},\\\
y=g(\xi;\mu)&(\xi,t,\mu)\in\Gamma_{in}\times[0,T]\times\mathbb{R}^{M},\\\
y=0&(\xi,t,\mu)\in\Gamma_{w}\times[0,T]\times\mathbb{R}^{M},\\\
y=u(t)&(\xi,t,\mu)\in\Gamma_{u}\times[0,T]\times\mathbb{R}^{M},\\\
\nu\partial_{n}y-p\vec{n}=0&(\xi,t,\mu)\in\Gamma_{out}\times[0,T]\times\mathbb{R}^{M},\\\
y(\xi,0;\mu)=y_{0}(\xi;\mu)&(\xi,\mu)\in\Xi\times\mathbb{R}^{M},\\\
\end{cases}$ (40)

where $\nu>0$ is the viscosity parameter and

$g(\xi;\mu)=4\xi_{2}(1-\xi_{2})+\frac{1}{2}\sum_{k=1}^{M}k^{-\gamma}\sin(2\pi
k\xi_{2})\mu_{k},$ (41)

is the uncertain inflow with $\mu=(\mu_{1},\ldots,\mu_{M})$ being independent
random variables each distributed uniformly on $[-c,c]$, where the half-width
$c>0$ will be varied in the numerical tests. Note that the mean inflow is
given by

$\overline{g}(\xi)=4\xi_{2}(1-\xi_{2}).$ (42)

The initial condition is set as

$y_{0}(\xi;\mu)=\begin{cases}g(\xi;\mu),&\xi\in\Gamma_{in},\\\
0,&\mbox{otherwise.}\end{cases}$

The equation is posed on a backward facing step domain, illustrated in Figure
9.

$\Gamma_{in}$$\Gamma_{w}$$\Gamma_{w}$$\Gamma_{u}$$\Gamma_{w}$$\Gamma_{out}$$g$$(-1,1)$$(5,1)$$(0,-1)$$(0,0)$

Figure 9: The backward step domain for the Navier-Stokes flow with uncertain
inflow

The control $u(t)$ is taken in the piecewise-constant form

$u(t)=\sum_{i=0}^{n_{t}}\begin{bmatrix}u_{1}(t_{i})\\\
u_{2}(t_{i})\end{bmatrix}\chi_{[t_{i},t_{i+1})}(t),$

where $\\{t_{i}\\}_{i=0}^{n_{t}}$ is a uniform discretization of the time
interval $[0,T]$.

We introduce the following cost functional

$J_{T}(y,u;\mu)=\int_{0}^{T}\int_{\Xi}|\nabla\times
y(\xi,t,\mu)|^{2}\,d\xi\,dt+\int_{0}^{T}\delta|u(t)|^{2}\;dt$ (43)

with the aim of minimizing the vorticity of the flow over a time interval
$[0,T]$ with a penalty cost for the control weighted by the parameter
$\delta>0$.

Since both the control and the random field appear in the boundary conditions,
it is convenient to split the solution in the form

$y=\widetilde{y}+\underline{y},$ (44)

where $\widetilde{y}$ takes into account the boundary conditions, and
$\underline{y}$ has homogeneous boundary conditions, and formulate the
feedback control problem on $\underline{y}$. In turn, $\widetilde{y}$ needs to
satisfy only the boundary and divergence-free conditions, so we choose it as
the solution of the stationary Stokes equation

$\begin{cases}-\nu\Delta\widetilde{y}+\nabla\widetilde{p}=0&(\xi,\mu)\in\Xi\times\mathbb{R}^{M},\\\
\nabla\cdot\widetilde{y}=0&(\xi,\mu)\in\Xi\times\mathbb{R}^{M},\\\
\widetilde{y}=g(\xi;\mu)&(\xi,\mu)\in\Gamma_{in}\times\mathbb{R}^{M},\\\
\widetilde{y}=0&(\xi,\mu)\in\Gamma_{w}\times\mathbb{R}^{M},\\\
\widetilde{y}=u&(\xi,\mu)\in\Gamma_{u}\times\mathbb{R}^{M},\\\
\nu\partial_{n}\widetilde{y}-p\vec{n}=0&(\xi,\mu)\in\Gamma_{out}\times\mathbb{R}^{M}.\\\
\end{cases}$ (45)

Note that $\widetilde{y}$ is a linear function of both $u$ and $g(\xi;\mu)$,
which is in turn a linear function of $\mu$ due to (41). This allows us to use
superposition and write the Stokes solution in the form of a linear map,

$\widetilde{y}(\widetilde{\mu},u)=\widetilde{Y}\widetilde{\mu}+Uu,\qquad\widetilde{Y}=\begin{bmatrix}\widetilde{Y}_{0}&\ldots&\widetilde{Y}_{M}\end{bmatrix},\quad
U=\begin{bmatrix}U_{1}&U_{2}\end{bmatrix},$ (46)

where $\widetilde{Y}_{0}$ is the solution of (45) with $\mu=0$ and $u=0$ (that
is, the uncontrolled Stokes solution with the mean inflow),
$\widetilde{Y}_{j}$ for $j=1,\ldots,M$ is the solution of (45) with
$\mu_{j}=1$, $\mu_{i}=0$ for $i\neq j$ and $u=0$,

$\widetilde{\mu}=\begin{bmatrix}1&\mu_{1}&\cdots&\mu_{M}\end{bmatrix}^{\top}$

is the random vector augmented by $1$ for brevity of the map (46), and $U_{j}$
is the solution of (45) with $g(\xi;\mu)=0$, $u_{j}=1$, $u_{i}=0$ for $i\neq
j$. As a by-product, this allows us to precompute (45) without time dependence
of $u(t)$, only for those $M+3$ initial inputs.

Plugging (46) into (40) we can write the following Navier-Stokes-type equation
on $\underline{y}$:

$\begin{cases}\partial_{t}\underline{y}-\nu\Delta\underline{y}+\widetilde{y}\cdot\nabla\underline{y}+\underline{y}\cdot\nabla\widetilde{y}+\underline{y}\cdot\nabla\underline{y}+\nabla\underline{p}=\nu\Delta\widetilde{y}-\widetilde{y}\cdot\nabla\widetilde{y},\\\
\nabla\cdot\underline{y}=0,\end{cases}$ (47)

with all boundary and initial conditions homogeneous.

###  6.2 Semi-discretization

We discretize (45) and (47) in space using the stable $P_{2}-P_{1}$ Taylor-
Hood finite elements pair, involving bilinear elements
$\\{\varphi_{i}\\}_{i=1}^{d_{p}}$ for the pressure and biquadratic elements
$\\{\phi_{i}\\}_{i=1}^{d_{v}}$ for the velocity, duplicated such that
$i=1,\ldots,d_{v}/2$ indexes the first component of the velocity, and
$i=d_{v}/2+1,\ldots,d_{v}$ indexes the second component. By convention,
$\phi_{i}$ and $\phi_{j}$ for any $i\leq d_{v}/2$ and $j>d_{v}/2$ are assumed
non-overlapping. The resulting semidiscretization of (47) reads

$\begin{cases}E\dot{\underline{y}}(t)+A\underline{y}(t)+F(\widetilde{y})\underline{y}(t)+F^{*}(\widetilde{y})\underline{y}(t)+F(\underline{y})\underline{y}(t)+D^{\top}p(t)=f_{0}(t),\\\
D\underline{y}(t)=0,\end{cases}$ (48)

where $E,A,F(\cdot),F^{*}(\cdot)$ and $D$ are mass and stiffness matrices with
elements

$\displaystyle(E)_{i,j}$ $\displaystyle=\int_{D}\phi_{i}\phi_{j}\,d\xi,$
$\displaystyle i,j$ $\displaystyle=1,\ldots,d_{v},$ $\displaystyle(A)_{i,j}$
$\displaystyle=\int_{D}\nu\nabla\phi_{i}\cdot\nabla\phi_{j}\,d\xi,$
$\displaystyle i,j$ $\displaystyle=1,\ldots,d_{v},$
$\displaystyle(F(v))_{i,j}$
$\displaystyle=\int_{D}\phi_{i}(v\cdot\nabla)\phi_{j}\,d\xi,$ $\displaystyle
i,j$ $\displaystyle=1,\ldots,d_{v},$ $\displaystyle(F^{*}(v))_{i,j}$
$\displaystyle=\int_{D}\phi_{i}\phi_{j}\partial_{\xi_{k}}v\,d\xi,$
$\displaystyle i,j$ $\displaystyle=1,\ldots,d_{v},\;k=\begin{cases}1,&j\leq
d_{v}/2,\\\ 2,&j>d_{v}/2,\end{cases}$ $\displaystyle(D)_{i,j}$
$\displaystyle=(D_{1})_{i,j}+(D_{2})_{i,j},$ where
$\displaystyle(D_{k})_{i,j}$
$\displaystyle=\int_{D}\varphi_{i}\partial_{\xi_{k}}\phi_{j}\,dx,$
$\displaystyle k$ $\displaystyle=1,2,\;i=1,\ldots,d_{p},\;j=1,\ldots,d_{v},$
$\displaystyle f_{0}(t)$
$\displaystyle=-A\widetilde{y}(\widetilde{\mu},u(t))-F(\widetilde{y})\widetilde{y}.$

Note that $\dot{\widetilde{y}}=0$ almost everywhere since $u(t)$ is piecewise-
constant. To obtain a feedback form, we can further separate $\widetilde{\mu}$
and $u$:

$\begin{cases}E\dot{\underline{y}}(t)+A\underline{y}(t)+C(\widetilde{\mu})\underline{y}(t)\\\
+F_{yu}(\underline{y}\otimes
u)+F(\underline{y})\underline{y}(t)+D^{\top}p(t)=f(\widetilde{\mu})-B(\widetilde{\mu})u-F_{uu}(u\otimes
u),\\\ D\underline{y}(t)=0,\end{cases}$ (49)

where

$\displaystyle C(\widetilde{\mu})$
$\displaystyle=\sum_{j=0}^{M}\widetilde{\mu}_{j}(F(\widetilde{Y}_{j})+F^{*}(\widetilde{Y}_{j})),$
$\displaystyle(F_{yu})_{i,k+2(j-1)}$
$\displaystyle=(F(U_{k}))_{i,j}+(F^{*}(U_{k}))_{i,j},$ $\displaystyle i,j$
$\displaystyle=1,\ldots,d_{v},\;k=1,2,$
$\displaystyle(B(\widetilde{\mu}))_{i,k}$
$\displaystyle=(AU)_{i,k}+\sum_{j=0}^{M}\widetilde{\mu}_{j}(F(\widetilde{Y}_{j})U_{k}+F(U_{k})\widetilde{Y}_{j})_{i},$
$\displaystyle i$ $\displaystyle=1,\ldots,d_{v},\;k=1,2,$
$\displaystyle(F_{uu})_{i,j+2(k-1)}$ $\displaystyle=(F(U_{j})U_{k})_{i},$
$\displaystyle i$ $\displaystyle=1,\ldots,d_{v},\;j,k=1,2,$ $\displaystyle
f(\widetilde{\mu})$
$\displaystyle=-A\widetilde{Y}\widetilde{\mu}-\sum_{j=0}^{M}\widetilde{\mu}_{j}F(\widetilde{Y}_{j})\widetilde{Y}\widetilde{\mu}.$

Applying the same semidiscretization to the cost functional, we obtain

$\displaystyle\tilde{J}_{T}(y,u)$
$\displaystyle=\int_{0}^{T}y(t)^{\top}D_{vort}y(t)+\delta|u(t)|^{2}\;dt$ (50)
$\displaystyle=\int_{0}^{T}\underline{y}(t)^{\top}D_{vort}\underline{y}(t)+u(t)^{\top}(\delta
I_{2}+U^{T}D_{vort}U)u(t)dt$
$\displaystyle+\int_{0}^{T}2\widetilde{\mu}^{\top}\widetilde{Y}^{\top}D_{vort}\underline{y}(t)+2u(t)^{\top}U^{\top}D_{vort}\underline{y}(t)dt,$

where $I_{2}\in\mathbb{R}^{2\times 2}$ is the identity matrix and

$D_{vort}=(D_{2}-D_{1})^{\top}(D_{2}-D_{1}).$

The system of ODEs (49) is approximated via an implicit Euler scheme with
$n_{t}$ time steps, and the resulting nonlinear algebraic system is solved via
a Newton’s method with stopping tolerance $10^{-6}$.

By default, we fix $M=8$, $\gamma=3$, $\nu=2\cdot 10^{-3}$, $\delta=10^{-3}$,
$T=20$, $n_{t}=80$, $d_{v}=10382$, $d_{p}=1340$ obtaining a problem of
dimension $d=11722$.

###  6.3 Model Order Reduction

Next, we pass to the construction of the reduced basis and the corresponding
reduced dynamical system. Fixing a realization $\mu$, the corresponding
snapshot matrix $Y_{\mu}$ will be formed just from the divergence-free
velocity snapshots $\\{\underline{y}^{0},\ldots,\underline{y}^{n_{t}}\\}$ with
homogeneous boundary conditions. Since the solution of the ROM dynamics is a
linear combination of the snapshots, the reduced trajectory directly benefits
from the divergence-free property and homogeneous boundary conditions and it
solves following reduced ODEs system

$\begin{cases}E^{\ell}\dot{y}^{\ell}(t)+A^{\ell}y^{\ell}(t)+C^{\ell}(\widetilde{\mu})y^{\ell}(t)\\\
+F^{\ell}_{yu}(y^{\ell}(t)\otimes
u)+F^{\ell}(y^{\ell}(t))y^{\ell}(t)=f^{\ell}({\widetilde{\mu}})-B^{\ell}(\widetilde{\mu})u-F^{\ell}_{uu}(u\otimes
u),\\\ y^{\ell}(0)=x^{\ell},\end{cases}$ (51)

where

$\displaystyle E^{\ell}$ $\displaystyle=(U^{\ell})^{\top}EU^{\ell},$
$\displaystyle A^{\ell}$ $\displaystyle=(U^{\ell})^{\top}AU^{\ell},$
$\displaystyle C^{\ell}(\widetilde{\mu})$
$\displaystyle=\sum_{j=0}^{M}\widetilde{\mu}_{j}\left[(U^{\ell})^{\top}F(\widetilde{Y}_{j})U^{\ell}+(U^{\ell})^{\top}F^{*}(\widetilde{Y}_{j})U^{\ell}\right],$
$\displaystyle F^{\ell}_{yu}$
$\displaystyle=(U^{\ell})^{\top}F_{yu}(U^{\ell}\otimes I_{2}),$ $\displaystyle
F^{\ell}(y^{\ell})$
$\displaystyle=\sum_{k=1}^{\ell}y^{\ell}_{k}\left[(U^{\ell})^{\top}F(U^{\ell}_{k})U^{\ell}\right],$
$\displaystyle(B^{\ell}(\widetilde{\mu}))_{k}$
$\displaystyle=(U^{\ell})^{\top}AU+\sum_{j=0}^{M}\widetilde{\mu}_{j}\left[(U^{\ell})^{\top}F(\widetilde{Y}_{j})U_{k}+(U^{\ell})^{\top}F(U_{k})\widetilde{Y}_{j}\right],$
$\displaystyle k=1,2,$ $\displaystyle F_{uu}^{\ell}$
$\displaystyle=(U^{\ell})^{\top}F_{uu},$ $\displaystyle
f^{\ell}({\widetilde{\mu}})$
$\displaystyle=-(U^{\ell})^{\top}A\widetilde{Y}\widetilde{\mu}-\sum_{j=0}^{M}\widetilde{\mu}_{j}(U^{\ell})^{\top}F(\widetilde{Y}_{j})\widetilde{Y}\widetilde{\mu},$

while the reduced cost functional reads

$\tilde{J}^{\ell}_{T}(y^{\ell},u)=\int_{0}^{T}y^{\ell}(t)^{\top}D^{\ell}_{vort}y^{\ell}(t)+u(t)^{\top}Ru(t)+2\widetilde{\mu}^{\top}Y^{\ell}_{vort}y^{\ell}(t)+2u(t)^{\top}U^{\ell}_{vort}y^{\ell}(t)dt,$

with

$\displaystyle D^{\ell}_{vort}$
$\displaystyle=(U^{\ell})^{\top}D_{vort}U^{\ell},$ $\displaystyle R$
$\displaystyle=\delta I_{2}+U^{T}D_{vort}U,$ $\displaystyle Y^{\ell}_{vort}$
$\displaystyle=\widetilde{Y}^{\top}D_{vort}U^{\ell}$ $\displaystyle
U^{\ell}_{vort}$ $\displaystyle=U^{\top}D_{vort}U^{\ell}.$

Note that coefficients involving $U^{\ell}$ can be precomputed, after which
the reduced ODEs system (51) can be both assembled and solved for each $\mu$
and $u$ with the complexity independent of the full dimension $d$.

Figure 10: Left: singular values of $Y_{\underline{\mu}}$ for different
numbers of offline realisations $N$. Right: model reduction errors (38) tested
on the mean inflow (42). The SPOD is built upon $N=15$ controlled trajectory
realisations considering random iid inflows from (41), and the 1 POD is built
upon the controlled trajectory with mean inflow (42).

###  6.4 Accuracy of the controlled reduced model

In the left panel of Figure 10 we show the decay of the singular values for 1
POD and for the SPOD using different number of realizations for the collection
of the snapshots ($N\in\\{5,10,15\\}$). We note that the singular values above
$10^{-3}$ are almost indistinguishable in SPOD using 10 and 15 realisations.
This shows that $N=15$ is sufficient for the statistical procedure to converge
for the most relevant modes. We now pass to study the performances of the SPOD
approach and its difference with the 1 POD technique. The statistical approach
has been constructed upon $N=15$ trajectories controlled by PMP considering
the random inflow (41) and fixing the parameter $c=1$. The right panel of
Figure 10 displays the behaviour of the different error indicators for a
system where the inflow is set to its mean value, $g=\overline{g}(\xi)$. It is
interesting to note that although 1 POD is built exactly upon the controlled
snapshots of the mean inflow, leading to a projection error of order $10^{-5}$
for $\ell=70$, the SPOD performs better for both error indicators for higher
reduced dimensions, reflecting the richness of the SPOD basis for
accommodating different trajectories.

Figure 11: Average errors over 10 iid random inflows with $c=1$ (left) and
$c=2$ (right) versus the reduced basis size $\ell$. The SPOD is built upon
$N=15$ controlled trajectory realisations with random iid inflows from (41).
The 1 POD is built upon 1 realisation of the controlled trajectory with mean
inflow $\overline{g}(\xi)$.

In Figure 11 we show the analysis of the mean of the errors (38) on 10 random
samples of $\mu$ changing the half-width $c$ of the interval of $\mu_{k}$. In
the left panel we consider $c=1$, employed also for the construction of the
statistical basis. The results are immediately evident: all the error
indicators for the 1 POD resolutions are stuck between order $10^{-2}$ and
$10^{-3}$, while SPOD presents in general a decreasing behaviour, reaching
order $10^{-5}$ for the error indicator $\mathcal{E}_{J}$ with 70 reduced
basis vectors. In the right panel the numerical experiments are run with
parameter $c=2$, introducing optimal trajectories possibly far from those
considered during the basis construction. This is reflected in the order of
the different error indicators, which achieve at most order $10^{-2}$, but
still yielding a better approximation compared to the 1 POD strategy.

###  6.5 Suboptimal controllers

Lastly, we compare different faster controllers: the PMP applied to the full
order model but with the mean inflow, the SDRE applied to the reduced order
model using 1 POD and SPOD bases, as well as the LQR applied to reduced
models. Firstly, we show the uncontrolled flow at the final time in Figure 12
with a prefixed random inflow (41) with

$\mu_{*}=\begin{bmatrix}\;\;\;0.0984\\\ -0.3838\\\ -0.1259\\\ \;\;\;0.0398\\\
-0.4314\\\ \;\;\;0.1589\\\ \;\;\;0.2323\\\ \;\;\;0.0947\end{bmatrix}.$ (52)

On the left panel we show the absolute velocity
$|v|(\xi)=\sqrt{y_{1}(\xi)^{2}+y_{2}(\xi)^{2}}$, while on the right panel we
show the velocity vector field $y=(y_{1}(\xi),y_{2}(\xi))$. It is possible to
notice the presence of vortexes due to re-circulation issues. The total cost
of the uncontrolled dynamics is equal to $4.6279$.

Figure 12: (Uncontrolled case) Absolute velocity (left) and velocity vector
field (right) at final time $T=20$. $J_{T}(y,u)=4.6279$.

Using PMP on the full system with the mean inflow $\overline{g}(\xi)$ to
compute the control signal, but applying this signal to the system with a
random realisation of the inflow gives the flow as shown in Figure 13. We see
that this control is unable to reduce the vortexes completely. This indicates
the need for controllers that are more specific and robust to random inputs to
the model.

Figure 13: (Deterministic controller) Mean velocity (left) and velocity vector
field (right) at final time $T=20$. The optimal control is computed via
Pontryagin for the Full Order dynamics in absence of noise and plugged into a
perturbed dynamical system with random inflow (41). $J_{T}(y,u)=4.9015$.

Now consider applying SDRE and LQR to the reduced system (51). Since SDRE does
not take into account quadratic terms in the control, we must omit them, and
write down the reduced dynamics in a semilinear form

$\dot{y}^{\ell}(t)=\mathcal{A}^{\ell}(y^{\ell}(t))y^{\ell}(t)+\mathcal{B}^{\ell}(y^{\ell}(t))u(t),\quad
y^{\ell}(0)=x^{\ell},$ (53)

where

$\mathcal{A}(y^{\ell}(t))=-(E^{\ell})^{-1}(A^{\ell}+C^{\ell}(\widetilde{\mu})+F^{\ell}(y^{\ell}(t))),$
$\mathcal{B}(y^{\ell}(t))=-(E^{\ell})^{-1}(F^{\ell}_{yu}(y^{\ell}(t)\otimes
I_{2})+B^{\ell}(\widetilde{\mu})).$

First we compute the solution controlled via a LQR feedback, solving the
Riccati equation (15) for $x=\underline{0}$, obtaining the matrix $P_{0}$. At
this point we consider the linearized feedback map

$u(x^{\ell})=-R^{-1}(\mathcal{B}(\underline{0})^{\top}P_{0}+2U^{\ell}_{vort})x^{\ell}$

and the resulting total cost is $3.3228$. The final configuration and the
velocity vector field is shown in Figure 14. We note that the resulting flow
is less turbulent than the uncontrolled case, but the solution is still far
from the laminar regime.

Figure 14: (LQR controller) Mean velocity (left) and velocity vector field
(right) at final time $T=20$. $J_{T}(y,u)=3.3228$.

###  6.6 Tensor Train approximation of the reduced SDRE control

Now we consider the TT approximation of the feedback control function computed
by SDRE on 1 POD and SPOD reduced models. Since the TT Cross approximates
scalar functions, it is applied twice, one per each component of the control
(14). Moreover, we fix $\mu=\mu_{*}$ as defined in (52). This allows us to
approximate again a function $u(x^{\ell})$ depending on the reduced state
only. In these state variables, we consider a domain
$X^{\ell}=\bigtimes_{i=1}^{\ell}[a_{i},b_{i}]$, where the interval ranges are
chosen such that the domain contains all the reduced snapshots. Furthermore,
we fix $tol=10^{-3}$ and we consider $n=6$ Legendre basis functions per
dimension.

Table 7 displays the approximation error $\mathcal{E}_{TT}$ (39). The error in
the statistical framework is in the order of $10^{-2}$, while the
approximation error with the 1 POD basis is order $10^{-1}$, reflecting the
fact that the projection onto the 1 POD basis is not able to approximate
perturbed trajectories.

$\ell$ |  $err_{TT}$ 1 POD |  $err_{TT}$ SPOD
---|---|---
5 | $2.49e-1$ | $3.60e-2$
10 | $3.50e-1$ | $6.21e-3$
20 | $1.19e-1$ | $1.62e-2$
Table 7: Error in the TT approximation using POD and SPOD.

Finally, in Figures 15-16 we show the final configuration of the controlled
solution respectively for TT-SDRE-1POD and TT-SDRE-SPOD. For the 1 POD
approach, it is possible to note by the right panel of Figure 15 that a
turbulent regime is still active, while for the statistical approach (in
Figure 16) the fluid presents a more laminar behaviour.

Figure 15: (1 POD controller) Mean velocity (left) and velocity vector field
(right) at final time $T=20$. The optimal control is computed via Tensor Train
Cross and 1 POD with $\ell=20$ basis. $J_{T}(y,u)=3.0880$.

Figure 16: (SPOD controller) Mean velocity (left) and velocity vector field
(right) at final time $T=20$. The optimal control is computed via Tensor Train
Cross and SPOD with $\ell=20$ basis. $J_{T}(y,u)=2.9527$.

## Concluding remarks

We have developed a model order reduction method for synthesis of feedback
control laws for nonlinear, parameter-dependent dynamics, including fluid flow
problems. The reduction phase is inspired by POD techniques, requiring
sampling of the (sub)optimal control problem solutions for different
realizations of the random variables. Snapshots are compressed for the
construction of a statistical POD basis which minimizes the empirical risk.

The resulting reduced order model facilitates the construction of a data-
driven stabilizing feedback law in the tensor train format. The low-rank
tensor train structure enables the real-time implementation of a feedback
control for high-dimensional problems such as vorticity minimization in 2D
Navier-Stokes.

Future research directions include the design of robust $\mathcal{H}_{\infty}$
controllers, the development of greedy sampling strategies for both parameters
and initial conditions which can alleviate the computational cost of the
offline phase, and the training of higher dimensional surrogates using
physics-informed neural networks.

## Acknowledgements

This research was supported by the UK Engineering and Physical Sciences
Research Council New Horizons Grant EP/V04771X/1, the New Investigator Award
EP/T031255/1 and the Standard Grant EP/T024429/1.

## References

  * [1] G. Kirsten, V. Simoncini, Order reduction methods for solving Large-Scale Differential Matrix Riccati Equations, SIAM Journal on Scientific Computing 42 (4) (2020) A2182–A2205.  doi:10.1137/19m1264217.
URL http://dx.doi.org/10.1137/19M1264217

  * [2] P. Benner, Z. Bujanović, P. Kürschner, J. Saak, A numerical comparison of different solvers for Large-Scale, Continuous-Time Algebraic Riccati Equations and LQR problems, SIAM Journal on Scientific Computing 42 (2) (2020) A957–A996.  doi:10.1137/18m1220960.
URL http://dx.doi.org/10.1137/18M1220960

  * [3] M. Akian, S. Gaubert, A. Lakhoua, The max-plus finite element method for solving deterministic optimal control problems: Basic properties and convergence analysis, SIAM Journal on Control and Optimization 47 (2) (2008) 817–848.  doi:10.1137/060655286.
URL http://dx.doi.org/10.1137/060655286

  * [4] J. Darbon, P. M. Dower, T. Meng, Neural network architectures using min plus algebra for solving certain high dimensional optimal control problems and Hamilton-Jacobi PDEs, doi: 10.48550/ARXIV.2105.03336 (2021).
URL https://arxiv.org/abs/2105.03336

  * [5] J. Garcke, A. Kröner, Suboptimal feedback control of PDEs by solving HJB equations on adaptive sparse grids, Journal of Scientific Computing 70 (1) (2016) 1–28.  doi:10.1007/s10915-016-0240-7.
URL http://dx.doi.org/10.1007/s10915-016-0240-7

  * [6] B. Azmi, D. Kalise, K. Kunisch, Optimal feedback law recovery by gradient-augmented sparse polynomial regression, Journal of Machine Learning Research 22 (48) (2021) 1–32.
  * [7] K. Kunisch, D. Vásquez-Varas, D. Walter, Learning optimal feedback operators and their sparse polynomial approximations, Journal of Machine Learning Research 24 (301) (2023) 1–38.
URL http://jmlr.org/papers/v24/22-1010.html

  * [8] A. Alla, M. Falcone, L. Saluzzi, High-order approximation of the finite horizon control problem via a tree structure algorithm, IFAC-PapersOnLine 52 (2) (2019) 19–24.
  * [9] M. Falcone, G. Kirsten, L. Saluzzi, Approximation of Optimal Control Problems for the Navier-Stokes equation via multilinear HJB-POD, Applied Mathematics and Computation 442 (2023) 127722.
  * [10] J. Han, A. Jentzen, W. E, Solving high-dimensional partial differential equations using deep learning, Proceedings of the National Academy of Sciences 115 (34) (2018) 8505–8510.  doi:10.1073/pnas.1718942115.
URL http://dx.doi.org/10.1073/pnas.1718942115

  * [11] J. Darbon, G. P. Langlois, T. Meng, Overcoming the curse of dimensionality for some Hamilton–Jacobi partial differential equations via neural network architectures, Research in the Mathematical Sciences 7 (3) (2020).  doi:10.1007/s40687-020-00215-6.
URL http://dx.doi.org/10.1007/s40687-020-00215-6

  * [12] K. Kunisch, D. Walter, Semiglobal optimal feedback stabilization of autonomous systems via deep neural network approximation, ESAIM: Control, Optimisation and Calculus of Variations 27 (2021) 16.  doi:10.1051/cocv/2021009.
URL http://dx.doi.org/10.1051/cocv/2021009

  * [13] T. Meng, Z. Zhang, J. Darbon, G. E. Karniadakis, Sympocnet: Solving optimal control problems with applications to high-dimensional multi-agent path planning problems, doi: 10.48550/ARXIV.2201.05475 (2022).
URL https://arxiv.org/abs/2201.05475

  * [14] M. Zhou, J. Han, J. Lu, Actor-critic method for high dimensional static hamilton–jacobi–bellman partial differential equations based on neural networks, SIAM Journal on Scientific Computing 43 (6) (2021) A4043–A4066.  doi:10.1137/21m1402303.
URL https://doi.org/10.1137%2F21m1402303

  * [15] D. Onken, L. Nurbekyan, X. Li, S. W. Fung, S. Osher, L. Ruthotto, A neural network approach applied to multi-agent optimal control, in: 2021 European Control Conference (ECC), IEEE, 2021.  doi:10.23919/ecc54610.2021.9655103.
  * [16] L. Ruthotto, S. J. Osher, W. Li, L. Nurbekyan, S. W. Fung, A machine learning framework for solving high-dimensional mean field game and mean field control problems, Proceedings of the National Academy of Sciences 117 (17) (2020) 9183–9193.
  * [17] G. Albi, S. Bicego, D. Kalise, Gradient-augmented Supervised Learning of Optimal Feedback Laws Using State-Dependent Riccat Equations, IEEE Control Systems Letters 6 (2022) 836–841.
  * [18] L. Sallandt, Computing high-dimensional value functions of optimal feedback control problems using the tensor-train format, Ph.D. thesis, Technische Universität Berlin (2022).
  * [19] L. Richter, L. Sallandt, N. Nüsken, Solving high-dimensional parabolic pdes using the tensor train format, in: International Conference on Machine Learning, PMLR, 2021, pp. 8998–9009.
  * [20] S. Dolgov, D. Kalise, K. K. Kunisch, Tensor Decomposition Methods for High-dimensional Hamilton–Jacobi–Bellman Equations, SIAM Journal on Scientific Computing 43 (3) (2021) A1625–A1650.  doi:10.1137/19m1305136.
  * [21] A. Alla, H. Oliveira, G. Santin, HJB-RBF based approach for the control of PDEs, Journal of Scientific Computing 96 (1) (2023) 25.
  * [22] T. Ehring, B. Haasdonk, Hermite kernel surrogates for the value function of high-dimensional nonlinear optimal control problems (2023).  arXiv:2305.06122.
  * [23] W. Kang, L. C. Wilcox, Mitigating the curse of dimensionality: sparse grid characteristics method for optimal feedback control and HJB equations, Computational Optimization and Applications 68 (2) (2017) 289–315.  doi:10.1007/s10589-017-9910-0.
URL http://dx.doi.org/10.1007/s10589-017-9910-0

  * [24] W. Kang, L. Wilcox, A causality free computational method for HJB equations with application to rigid body satellites, AIAA Guidance, Navigation, and Control Conference (2015).  doi:10.2514/6.2015-2009.
URL http://dx.doi.org/10.2514/6.2015-2009

  * [25] B. Azmi, D. Kalise, K. Kunisch, Optimal feedback law recovery by gradient-augmented sparse polynomial regression, Journal of Machine Learning Research 22 (2021).
URL
http://gateway.webofknowledge.com/gateway/Gateway.cgi?GWVersion=2&SrcApp=PARTNER_APP&SrcAuth=LinksAMR&KeyUT=WOS:000628529500048&DestLinkType=FullRecord&DestApp=ALL_WOS&UsrCustomerID=1ba7043ffcc86c417c072aa74d649202

  * [26] M. Oster, L. Sallandt, R. Schneider, Approximating optimal feedback controllers of finite horizon control problems using hierarchical tensor formats, SIAM Journal on Scientific Computing 44 (3) (2022) B746–B770.
  * [27] S. Dolgov, D. Kalise, L. Saluzzi, Data-driven tensor train gradient cross approximation for hamilton–jacobi–bellman equations, SIAM Journal on Scientific Computing 45 (5) (2023) A2153–A2184.
  * [28] T. Nakamura-Zimmerer, Q. Gong, W. Kang, Adaptive deep learning for high-dimensional Hamilton–Jacobi–Bellman equations, SIAM Journal on Scientific Computing 43 (2) (2021) A1221–A1247.  doi:10.1137/19m1288802.
URL https://doi.org/10.1137%2F19m1288802

  * [29] T. Nakamura-Zimmerer, Q. Gong, W. Kang, QRnet: Optimal regulator design with LQR-augmented neural networks, IEEE Control Systems Letters 5 (4) (2021) 1303–1308.  doi:10.1109/lcsys.2020.3034415.
URL https://doi.org/10.1109%2Flcsys.2020.3034415

  * [30] A. Bensoussan, J. Han, S. C. P. Yam, X. Zhou, Value-gradient based formulation of optimal control problem and machine learning algorithm, doi: 10.48550/arXiv.2103.09280 (2021).
URL https://arxiv.org/abs/2103.09280

  * [31] D. Onken, L. Nurbekyan, X. Li, S. W. Fung, S. Osher, L. Ruthotto, A neural network approach for high-dimensional optimal control, arXiv preprint arXiv:2104.03270 (2021).
  * [32] K. Kunisch, S. Volkwein, L. Xie, HJB-POD-based feedback design for the optimal control of evolution problems, SIAM Journal on Applied Dynamical Systems 3 (4) (2004) 701–722.
  * [33] K. Kunisch, L. Xie, POD-based feedback control of the Burgers equation by solving the evolutionary HJB equation, Computers & Mathematics with Applications 49 (7-8) (2005) 1113–1126.
  * [34] B. Moore, Principal component analysis in linear systems: Controllability, observability, and model reduction, IEEE Transactions on Automatic Control 26 (1) (1981) 17–32.  doi:10.1109/TAC.1981.1102568.
  * [35] P. M. R. Wortelboer, Frequency weighted balanced reduction of closed-loop mechanical servo-systems, Ph.D. thesis, delft University of Technology (1994).
  * [36] P. Benner, T. Damm, Lyapunov equations, energy functionals, and model order reduction of bilinear and stochastic systems, SIAM J. Control Optim. 49 (2) (2011) 686–711.  doi:10.1137/09075041X.
URL https://doi.org/10.1137/09075041X

  * [37] B. Kramer, K. Willcox, Balanced Truncation Model Reduction for Lifted Nonlinear Systems, Springer International Publishing, 2022, pp. 157–174.  doi:10.1007/978-3-030-95157-3_9.
  * [38] P. Benner, J. Heiland, S. W. R. Werner, Robust output-feedback stabilization for incompressible flows using low-dimensional $\mathcal{H}_{\infty}$-controllers, Computational Optimization and Applications 82 (1) (2022) 225–249.  doi:10.1007/s10589-022-00359-x.
  * [39] B. Kramer, S. Gugercin, J. Borggaard, L. Balicki, Nonlinear balanced truncation: Part 1-computing energy functions (2022).  arXiv:2209.07645.
  * [40] B. Kramer, S. Gugercin, J. Borggaard, Nonlinear balanced truncation: Part 2 – model reduction on manifolds (2023).  arXiv:2302.02036.
  * [41] P. Benner, M. Redmann, Model reduction for stochastic systems, Stochastic Partial Differential Equations: Analysis and Computations 3 (3) (2015) 291–338.  doi:10.1007/s40072-015-0050-1.
  * [42] S. Fresca, L. Dede, A. Manzoni, A comprehensive deep learning-based approach to reduced order modeling of nonlinear time-dependent parametrized pdes, Journal of Scientific Computing 87 (61) (2021) 1–36.
URL https://doi.org/10.1007/s10915-021-01462-7

  * [43] S. Lall, J. E. Marsden, S. Glavaški, Empirical model reduction of controlled nonlinear systems, IFAC Proceedings Volumes 32 (2) (1999) 2598–2603, 14th IFAC World Congress 1999, Beijing, Chia, 5-9 July.  doi:10.1016/S1474-6670(17)56442-3.
  * [44] C. W. Rowley, Model reduction for fluids, using balanced proper orthogonal decomposition, International Journal of Bifurcation and Chaos 15 (03) (2005) 997–1013.  doi:10.1142/S0218127405012429.
  * [45] K. Willcox, J. Peraire, Balanced model reduction via the proper orthogonal decomposition, AIAA Journal 40 (11) (2002) 2323–2330.  doi:10.2514/2.1570.
  * [46] T. Cui, K. J. H. Law, Y. M. Marzouk, Dimension-independent likelihood-informed mcmc, Journal of Computational Physics 304 (2016) 109–137.
  * [47] P. G. Constantine, E. Dow, Q. Wang, Active subspace methods in theory and practice: Applications to kriging surfaces, SIAM Journal on Scientific Computing 36 (4) (2014) A1500–A1524.
  * [48] O. Zahm, P. G. Constantine, C. Prieur, Y. Marzouk, Gradient-based dimension reduction of multivariate vector-valued functions, SIAM Journal on Scientific Computing 42 (1) (2020) A534–A558.
  * [49] D. E. Kirk, Optimal Control Theory, Dover Publications, 2004.
  * [50] H. T. Banks, B. M. Lewis, H. T. Tran, Nonlinear feedback controllers and compensators: a state-dependent Riccati equation approach, Computational Optimization and Applications 37 (2) (2007) 177–218.  doi:10.1007/s10589-007-9015-2.
URL http://dx.doi.org/10.1007/s10589-007-9015-2

  * [51] A. Yousefi, B. Lohmann, Balancing & optimization for order reduction of nonlinear systems, in: Proceedings of the 2004 American Control Conference, Vol. 1, 2004, pp. 108–112.  doi:10.23919/ACC.2004.1383588.
  * [52] A. V. Mamonov, M. A. Olshanskii, Tensorial parametric model order reduction of nonlinear dynamical systems, arXiv preprint 2302.08490 (2023).
URL http://arxiv.org/abs/2302.08490

  * [53] I. V. Oseledets, Tensor-train decomposition, SIAM J. Sci. Comput. 33 (5) (2011) 2295–2317.  doi:10.1137/090752286.
  * [54] D. Bigoni, A. P. Engsig-Karup, Y. M. Marzouk, Spectral Tensor-Train decomposition, SIAM Journal on Scientific Computing 38 (4) (2016) A2405–A2439.  doi:10.1137/15M1036919.
  * [55] A. Gorodetsky, S. Karaman, Y. Marzouk, A continuous analogue of the tensor-train decomposition, Computer Methods in Applied Mechanics and Engineering 347 (2019) 59–84.  doi:10.1016/j.cma.2018.12.015.
  * [56] T. Çimen, State-dependent Riccati equation (SDRE) control: a survey, IFAC Proceedings Volumes 41 (2) (2008) 3761–3775.
  * [57] A. Alla, D. Kalise, V. Simoncini, State-dependent Riccati equation feedback stabilization for nonlinear PDEs, doi: 10.48550/ARXIV.2106.07163 (2021).
URL https://arxiv.org/abs/2106.07163

  * [58] P. B. Rohrbach, S. Dolgov, L. Grasedyck, R. Scheichl, Rank bounds for approximating Gaussian densities in the Tensor-Train format, SIAM/ASA Journal on Uncertainty Quantification 10 (3) (2022) 1191–1224.  doi:10.1137/20M1314653.
  * [59] I. V. Oseledets, E. E. Tyrtyshnikov, TT-cross approximation for multidimensional arrays, Linear Algebra Appl. 432 (1) (2010) 70–88.  doi:10.1016/j.laa.2009.07.024.
  * [60] S. Dolgov, D. Savostyanov, Parallel cross interpolation for high–precision calculation of high–dimensional integrals, Comp. Phys. Comm. 246 (2020) 106869\.  doi:10.1016/j.cpc.2019.106869.
  * [61] W. Hackbusch, Tensor Spaces And Numerical Tensor Calculus, Springer–Verlag, Berlin, 2012.
  * [62] M. Barrault, Y. Maday, N. C. Nguyen, A. T. Patera, An ‘empirical interpolation’method: application to efficient reduced-basis discretization of partial differential equations, Comptes Rendus Mathematique 339 (9) (2004) 667–672.
  * [63] S. Chaturantabut, D. C. Sorensen, Nonlinear model reduction via discrete empirical interpolation, SIAM Journal on Scientific Computing 32 (5) (2010) 2737–2764.
  * [64] M. Heinkenschloss, Numerical solution of implicitly constrained optimization problems, Tech. rep. (2008).