∎

11institutetext: G. Kirsten 22institutetext: Dipartimento di
Matematica,Università di Bologna,Piazza di Porta S. Donato, 5, Bologna,
I-40127, Italy.
22email<EMAIL_ADDRESS>33institutetext: L. Saluzzi 44institutetext:
Department of Mathematics, Imperial College London, South Kensington Campus,
London, SW7 2AZ, UK.
44email<EMAIL_ADDRESS>
# A multilinear HJB-POD method for
the optimal control of PDEs In loving memory of Maurizio Falcone

Gerhard Kirsten  Luca Saluzzi

(Received: date / Accepted: date)

###### Abstract

Optimal control problems driven by evolutionary partial differential equations
arise in many industrial applications and their numerical solution is known to
be a challenging problem. One approach to obtain an optimal feedback control
is via the Dynamic Programming principle. Nevertheless, despite many
theoretical results, this method has been applied only to very special cases
since it suffers from the curse of dimensionality. Our goal is to mitigate
this crucial obstruction developing a new version of dynamic programming
algorithms based on a tree structure and exploiting the compact representation
of the dynamical systems based on tensors notations via a model reduction
approach. Here, we want to show how this algorithm can be constructed for
general nonlinear control problems and to illustrate its performances on a
number of challenging numerical tests. Our numerical results indicate a large
decrease in memory requirements, as well as computational time, for the
proposed problems. Moreover, we prove the convergence of the algorithm and
give some hints on its implementation.

###### Keywords:

dynamic programming, optimal control, tree structure, model order reduction,
error estimates

###### MSC:

49L20, 49J15, 49J20, 93B52

##  1 Introduction

Feedback control is a fundamental concept in engineering and applied
mathematics, where the goal is to design a system that can regulate a process
to achieve a desired behavior. One of the most powerful tools in feedback
control is the Hamilton Jacobi Bellman (HJB) equation, which provides a
framework for optimal control of dynamical systems. The HJB equation is a
partial differential equation that arises from the calculus of variations and
has a wide range of applications, including $e.g.$ robotics, aerospace and
finance. The main disadvantage of this approach comes in the form of the so-
called curse of dimensionality; the phenomenon for which the complexity of a
problem increases exponentially as the number of variables or dimensions
involved in the problem grows. In real applications, the dynamical system may
be described by a large number of state variables either because the
continuous problem is in high-dimension or because the dynamical system is
obtained via a discretization in space of a Partial Differential Equation
(PDE). In the context of linear dynamics and quadratic cost functional, the
HJB is equivalent to the Differential Riccati Equation for the finite horizon
control problem and to the Algebraic Riccati Equation in the infinite horizon
case. This setting has been widely researched, leading to several promising
high-dimensional solvers kirsten2019 ; BBKS_2020 .

For general nonlinear problems, such a reformulation does not exist and the
HJB equation must be tackled directly. In the recent years several efforts
have been employed in the mitigation of the curse of dimensionality arising in
optimal control, among those we mention sparse grids GK16 , max-plus algebra
McEneaney_2007 ; Akian_Gaubert_Lakhoua_2008 ; maxplusdarbon ;
akian2023adaptive , artificial neural networks Han_Jentzen_E_2018 ;
Darbon_Langlois_Meng_2020 ; Kunisch_Walter_2021 ; sympocnet ; Zhou_2021 ;
Onken2021 ; ruthotto2020machine , the application of tensor formats
dolgov2022data ; oster2022approximating ; richter2021solving  and radial basis
functions alla2021hjb .

In this paper we aim to mitigate the curse of dimensionality via a graph-based
optimization algorithm, the Tree Structure Algorithm (TSA) for the resolution
of the finite horizon HJB problem AFS19 . The TSA leads to the construction of
a tree in the direction of all the possible controlled trajectories. Due to
its flexible structure, this technique has been already applied in different
settings, $e.g.$ state constraints problems afs20  and high-dimensional
semidiscrete PDEs Alla_Saluzzi_2020  and its convergence is ensured by
rigorous error estimates saluzzi2022error . Furthermore, a geometrical pruning
based on the distance of the nodes has been introduced to avoid the
exponential growth of the tree and obtain a quadratic growth rate in the
context of LQR problems saluzzi2022error . Unfortunately, for general
nonlinear problems this criterion may be not effective since the tree nodes
may spread out faster, leading to a further curse of dimensionality. This is
the first shortcoming of the TSA that we will aim to address in this paper. We
will investigate:

  * •

a bilinear setting where the application of the geometrical pruning also
yields a good reduction in the cardinality of the tree,

  * •

an optimal control problem based on monotone controls with an efficient tree-
based data structure,

  * •

a statistical pruning rule based on the iterative knowledge of the value
function on the tree nodes.

A second shortcoming of the TSA that we address in this paper is related to
the computational cost of evaluating and constructing full-dimensional tree
nodes. Given the exponential growth in the cardinality of the tree, that can
merely be mitigated by pruning techniques, a massive computational effort may
be required to construct and evaluate the discrete problem on the tree nodes,
as the dimension of the discrete problem is increased.

A first attempt to address this issue was proposed in Alla_Saluzzi_2020  where
the authors applied a combination of the Proper Orthogonal Decomposition (POD)
volkwein2011model  for the linear part of the problem and Discrete Empirical
Interpolation method (DEIM) chaturantabut2010nonlinear  for the nonlinear
terms, to reduce the complexity of the problem. The coupling of the POD
technique and the HJB equation dates back to the pioneering paper by Kunisch
and co-authors KVX2004 , which was then further developed in a series of works
KX2005 ; kunisch2010optimal ; HV2005 . Nevertheless, in this setting, the POD-
DEIM algorithm itself has some computational drawbacks. More precisely, the
dimension of the vectors that need to be stored for constructing the tree and
the POD and DEIM spaces increase exponentially to the order of the dimension
of the underlying dynamical system. This may lead to a large bottleneck in the
computation of the reduced spaces for larger problems.

Instead, in this paper, we take advantage of the composite structure of a
subset of semilinear PDEs, where the underlying PDEs can be written and
discretized in Array form, leading to discrete semilinear Matrix and Tensor
equations Simoncini2017 ; Autilia2019matri ; palitta2016 ; kirsten.22 . Given
this particular structure, we aim to show how the tree can be constructed in
low-dimension by applying a Higher-Order POD-DEIM (HO-POD-DEIM) model order
reduction kirsten.22  to the discrete problem and then solving the HJB
equation on the low-dimensional tree. Our computational results on several
benchmark problems indicate that the new algorithm leads to memory
requirements that increase linearly in the dimension of the underlying
dynamical system, instead of exponentially. Furthermore, the convergence of
the proposed technique is established by the derivation of rigorous error
estimates for the reduced discrete dynamical system and for the discrete value
function computed on the reduced tree. These theoretical results extend the
error bounds obtained in sorensen2016  to semi-implicit schemes as well as the
proposed HO-POD-DEIM technique.

Our construction focuses on general high-dimensional semidiscretized PDEs. A
simplified matrix-oriented version of our framework is experimentally explored
for the Navier-Stokes (NS) equation in the companion manuscript KSF2023 ,
where the application to systems of differential equations is also discussed.
This is a classical example known to be computationally expensive, where the
application of MOR techniques helps in the computation of the solution (see
$e.g.$ QR2007 ; SR2018 ; pichi2022driving ). Here we consider complex problems
in high dimension, showing the promising numerical results on benchmark
problems. Furthermore, we deepen the analysis of all the ingredients of this
new methodology, including pruning techniques, error estimates and important
implementational nuances. We believe that the numerical simulations presented
in the last section illustrate that DP is now also feasible for more complex,
higher-dimensional problems from a computational point of view, and we hope
that this brings it closer to the application of challenging industrial
problems.

The paper is organized as follows. In the second section we introduce the
optimal control framework and the Tree Structure Algorithm. Section 3 is
devoted to the Model Order Reduction setting and its coupling with the TSA,
whereas in Section 4 we present some hints for an efficient implementation of
the proposed algorithm. In Section 5 we examine different pruning criteria for
the TSA showing some results in the reduction of the cardinality of the tree,
and Section 6 presents an error bound for the approximation of the value
function via the reduced order model algorithm. Finally, in the last section
we present some numerical experiments to show the effectiveness of the
proposed method.

##  2 The optimal control problem

Let us consider the classical finite horizon optimal control problem that we
use as a model problem. The system is driven by

$\left\\{\begin{array}[]{l}\dot{y}(s)=f(y(s),u(s),s),\;\;s\in(t,T],\\\
y(t)=x\in\mathbb{R}^{N}.\end{array}\right.$ (1)

Here, $y:[t,T]\rightarrow\mathbb{R}^{N}$ is the solution,
$u:[t,T]\rightarrow\mathbb{R}^{m}$ is the control,
$f:\mathbb{R}^{N}\times\mathbb{R}^{m}\times[t,T]\rightarrow\mathbb{R}^{N}$ is
the dynamics and

$\mathcal{U}=\\{u:[t,T]\rightarrow U,\mbox{measurable}\\}$

is the set of admissible controls where $U\subset\mathbb{R}^{m}$ is a compact
set. We define the cost functional for the finite horizon optimal control
problem as

$J_{x,t}(u):=\int_{t}^{T}L(y(s,u),u(s),s)\,ds+g(y(T)),$ (2)

where $L:\mathbb{R}^{N}\times\mathbb{R}^{m}\times[t,T]\rightarrow\mathbb{R}$
is the running cost and $g:\mathbb{R}^{N}\rightarrow\mathbb{R}$ is the final
cost. In the present work we will assume that the functions $f,L$ and $g$ are
bounded:

$\displaystyle\begin{aligned} |f(x,u,s)|\leq M_{f},&\quad|L(x,u,s)|\leq
M_{L},\quad|g(x)|\leq M_{g},\cr&\forall\,x\in\mathbb{R}^{N},u\in
U\subset\mathbb{R}^{m},s\in[t,T],\end{aligned}$ (3)

the functions $f$ and $L$ are Lipschitz-continuous with respect to the first
variable

$\displaystyle\begin{aligned} &|f(x,u,s)-f(y,u,s)|\leq
L_{f}|x-y|,\quad|L(x,u,s)-L(y,u,s)|\leq
L_{L}|x-y|,\cr&\qquad\qquad\qquad\qquad\forall\,x,y\in\mathbb{R}^{N},u\in
U\subset\mathbb{R}^{m},s\in[t,T],\end{aligned}$ (4)

and finally the cost $g$ is also Lipschitz-continuous:

$|g(x)-g(y)|\leq L_{g}|x-y|,\quad\forall x,y\in\mathbb{R}^{N}.$ (5)

Note that these assumptions guarantee uniqueness for the trajectory $y(t)$ by
the Carathéodory theorem (we refer to e.g. BC08  for a precise statement).

The aim is to construct a state-feedback control law $u(t)=\Phi(y(t),t),$ in
terms of the state equation $y(t),$ where $\Phi$ is the feedback map. The
optimality conditions are derived via the well-known Dynamic Programming
Principle (DPP) introduced by R. Bellman. We first introduce the value
function for an initial datum $(x,t)\in\mathbb{R}^{N}\times[t,T]$:

$v(x,t):=\inf\limits_{u\in\mathcal{U}}J_{x,t}(u)$ (6)

which can be represented via the DPP, i.e. for every $\tau\in[t,T]$:

$v(x,t)=\inf_{u\in\mathcal{U}}\left\\{\int_{t}^{\tau}L(y(s),u(s),s)ds+v(y(\tau),\tau)\right\\}.$
(7)

Due to (7) the HJB can be derived for every $x\in\mathbb{R}^{N}$, $s\in[t,T)$:

$\left\\{\begin{array}[]{ll}&-\dfrac{\partial v}{\partial
s}(x,s)+\max\limits_{u\in U}\left\\{-L(x,u,s)-\nabla v(x,s)\cdot
f(x,u,s)\right\\}=0,\\\ &v(x,T)=g(x).\end{array}\right.$ (8)

Once the value function is known, by e.g. solving (8), then the optimal
feedback control can be obtained as:

$u^{*}(t):=\operatorname*{arg\,max}_{u\in U}\left\\{-L(x,u,t)-\nabla
v(x,t)\cdot f(x,u,t)\right\\}.$ (9)

###  2.1 Dynamic programming on a tree structure

We briefly sketch the essential features of the dynamic programming approach
on a tree based on the discrete approximation of the dynamical system. More
details on the tree structure algorithm can be found in AFS19  where the
algorithm and several tests have been presented.
It is hard to find analytical solutions of the HJB equation (8) due to the
nonlinearity and classical approximation methods, e.g. finite difference or
semi-Lagrangian schemes, need a space discretization that is impossible to
manage in high-dimension (see the book FF13  for a comprehensive analysis of
approximation schemes for Hamilton-Jacobi equations). This has motivated
different approaches to mitigate the ”curse of dimensionality”.

We consider the discretized problem with a time step $\Delta t:=[(T-t)/N_{t}]$
where $N_{t}$ is the number of temporal time steps

$\left\\{\begin{array}[]{ll}V^{n}(x)=\min\limits_{u\in U}[\Delta
t\,L(x,u,t_{n})+V^{n+1}(x+\Delta tf(x,u,t_{n}))],\quad n=N_{t}-1,\dots,0,\\\
V^{N_{t}}(x)=g(x),\hskip 213.39566ptx\in\mathbb{R}^{N},\end{array}\right.$
(10)

where $t_{n}=t+n\Delta t,\,t_{N_{t}}=T$, and $V^{n}(x):=V(x,t_{n}).$ The
classical approach computes the solution through the application of an
interpolation operator to obtain the term $V^{n+1}(x+\Delta tf(x,u,t_{n}))$
based on the values sitting on the grid nodes. This direction will be
abandoned to build a tree structure and computing (10) only on a tree
structure. Starting from the initial condition $x$, we consider all the nodes
obtained following the discrete dynamics, e.g. for the explicit Euler scheme
with different discrete controls $u_{j}$. This gives in one step the points

$\zeta_{j}^{1}=x+\Delta t\,f(x,u_{j},t_{0}),\qquad j=1,\ldots,M.$ (11)

We assume that the control set $U$ is a bounded subset in $\mathbb{R}^{m}$ and
we discretize the control domain $U$ with constant step-size $\Delta u$
obtaining a discrete control set with a finite number of points $U^{\Delta
u}=\\{u_{1},...,u_{M}\\}$ that in the sequel we continue to denote by $U$
(with a slight abuse of notation).
Therefore, from every point $x$ we can reach $M$ points by (11). Identifying
the root of the tree with $\mathcal{T}^{0}=\\{x\\}$ we obtain the first level
of the tree $\mathcal{T}^{1}=\\{\zeta_{1}^{1},\ldots,\zeta^{1}_{M}\\}$. We can
proceed in this way so that all the nodes at the $n-$th time level, will be
given by

$\mathcal{T}^{n}=\\{\zeta^{n-1}_{i}+\Delta
tf(\zeta^{n-1}_{i},u_{j},t_{n-1})\\}_{j=1}^{M}\quad i=1,\ldots,M^{n-1}$

and all the nodes belonging to the tree can be shortly defined as

$\mathcal{T}:=\\{\zeta_{j}^{n}\\}_{j=1}^{M^{n}},\quad n=0,\ldots N_{t},$

where the nodes $\zeta^{n}_{i}$ are the result of the dynamics at time $t_{n}$
with the controls $\\{u_{j_{k}}\\}_{k=0}^{n-1}$:

$\zeta_{i_{n}}^{n}=\zeta_{i_{n-1}}^{n-1}+\Delta
tf(\zeta_{i_{n-1}}^{n-1},u_{j_{n-1}},t_{n-1})=x+\Delta
t\sum_{k=0}^{n-1}f(\zeta^{k}_{i_{k}},u_{j_{k}},t_{k}),$

with $\zeta^{0}=x$, $i_{k}=\left\lfloor\dfrac{i_{k+1}}{M}\right\rfloor$ and
$j_{k}\equiv i_{k+1}\mbox{mod }M$, where
$\zeta_{i}^{k}\in\mathbb{R}^{N},i=1,\ldots,M^{k}$ and
$\left\lfloor\cdot\right\rfloor$ is the ceiling function.

Despite the fact that the tree structure allows the resolution of high
dimensional problems, the construction may be expensive since
$|\mathcal{T}|=O(M^{N_{t}})$, where $N_{t}$ the number of time steps and $M$
is the number of controls. Whenever $M$ or $N_{t}$ are too large, the
construction turns out to be infeasible due to the memory allocation. In
Section 5 we will introduce two pruning criteria and theoretical results on
the reduction of the cardinality, showing their efficiency in avoiding the
allocation memory problem.

Once the tree structure $\mathcal{T}$ has been constructed, we compute the
numerical value function $V(x,t)$ on the tree nodes as

$V(x,t_{n})=V^{n}(x),\quad\forall x\in\mathcal{T}^{n},$ (12)

where $t_{n}=t+n\Delta t$. It is now straightforward to evaluate the value
function. Since the TSA defines a grid
$\mathcal{T}^{n}=\\{\zeta^{n}_{j}\\}_{j=1}^{M^{n}}$ for $n=0,\ldots,N_{t}$, we
can approximate (8) as follows:

$\begin{cases}V^{n}(\zeta^{n}_{i})=\min\limits_{u\in
U}\\{V^{n+1}(\zeta^{n}_{i}+\Delta tf(\zeta^{n}_{i},u,t_{n}))+\Delta
t\,L(\zeta^{n}_{i},u,t_{n})\\},\\\
\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\zeta^{n}_{i}\in\mathcal{T}^{n}\,,n=N_{t}-1,\ldots,0,\\\
V^{N_{t}}(\zeta^{N_{t}}_{i})=g(\zeta_{i}^{N_{t}}),\qquad\qquad\qquad\qquad\qquad\qquad\zeta_{i}^{N_{t}}\in\mathcal{T}^{N_{t}},\end{cases}$
(13)

where the minimization is computed by comparison on the discretized set of
controls $U$.

##  3 Reduced order models on a tree structure

Despite the fact that the tree structure algorithm avoids the construction of
a grid in high dimensions, the resulting memory requirements can still be
overwhelming. A first step towards relieving this computational demand via
model order reduction was presented in Alla_Saluzzi_2020 . More precisely, the
POD-DEIM algorithm from chaturantabut2010nonlinear  is used to reduce the
dimension of the discrete dynamical system, so that the the tree construction
is performed in low dimension. Nevertheless, the POD-DEIM algorithm itself has
some computational drawbacks. Firstly, if the discrete dynamical system from
the finite difference semi-discretization of a PDE in dimension $d$, the
memory requirements in both the offline and online phases of POD-DEIM are of
${\cal O}(N)$, where $N=\prod_{i=1}^{d}n_{i}$, where $n_{i}$ is the number of
discretization nodes in the $i$th spatial direction. A similar increase in
memory requirements is experienced for other discretization techniques.

Instead, for high-dimensional semi-discrete PDEs, we couple the tree structure
algorithm with the multilinear POD-DEIM algorithm presented in
Kirsten.Simoncini.arxiv2020  for the 2D case and in kirsten.22  for higher
dimensions. This will decrease the memory requirements to ${\cal
O}(\widetilde{N})$, with $\widetilde{N}=\sum_{i=1}^{d}n_{i}$.

To this end, we will first discuss some basic tensor notation required for the
new algorithm, after which we will review the standard HJB-POD algorithm,
before introducing the new multilinear one.

###  3.1 Notation and tensor basics

The third mode of a third-order tensor ${\boldsymbol{\cal
T}}\in\mathbb{R}^{n_{1}\times n_{2}\times n_{3}}$, is given by (see e.g.,
kolda2009 )

${\boldsymbol{\cal
T}}_{(3)}=\begin{pmatrix}{\bm{T}}_{1},{\bm{T}}_{2},\cdots,{\bm{T}}_{n_{2}}\end{pmatrix},$

where ${\bm{T}}_{i}\in\mathbb{R}^{n_{3}\times n_{1}},i=1,2,\ldots,n_{2}$ is
referred to as a lateral slice, and ${\boldsymbol{\cal T}}_{(3)}$ is a matrix
in $\mathbb{R}^{n_{3}\times n_{1}n_{2}}$. Multiplication between a tensor and
a matrix, is done via the $m-$mode product, which, for a tensor
${\boldsymbol{\cal T}}\in\mathbb{R}^{n_{1}\times n_{2}\times n_{3}}$ and a
matrix ${\bm{M}}\in\mathbb{R}^{n\times n_{m}}$, we express as

${\boldsymbol{\cal Q}}={\boldsymbol{\cal
T}}\times_{m}{\bm{M}}\quad\iff\quad{\boldsymbol{\cal
Q}}_{(m)}={\bm{M}}{\boldsymbol{\cal T}}_{(m)},$

in the $m$-th mode. The Kronecker product of two matrices
${\bm{M}}\in\mathbb{R}^{m_{1}\times m_{2}}$ and
${\bm{N}}\in\mathbb{R}^{n_{1}\times n_{2}}$ is defined as

${\bm{M}}\otimes{\bm{N}}=\begin{pmatrix}M_{1,1}{\bm{N}}&\cdots&M_{1,m_{2}}{\bm{N}}\\\
\vdots&\ddots&\vdots\\\
M_{m_{1},1}{\bm{N}}&\cdots&M_{m_{2},m_{2}}{\bm{N}}\end{pmatrix}\in\mathbb{R}^{m_{1}n_{1}\times
m_{2}n_{2}},$

and the vec$(\cdot)$ operator stacks the columns of a matrix one after the
other to form a long vector. For a third order tensor, the vectorization is
applied via the first mode unfolding. Furthermore,

$({\bm{M}}\otimes{\bm{N}})\mbox{vec}({\bm{X}})=\mbox{vec}({\bm{N}}{\bm{X}}{\bm{M}}^{\top}).$
(14)

As a result, if ${\boldsymbol{\cal X}}\in\mathbb{R}^{n_{1}\times n_{2}\times
n_{3}}$, and ${\bm{X}}={\boldsymbol{\cal X}}_{(3)}^{\top}$, then

$({\bm{L}}\otimes{\bm{M}}\otimes{\bm{N}})\mbox{vec}\left({\boldsymbol{\cal
X}}\right)=\mbox{vec}\left(({\bm{M}}\otimes{\bm{N}}){\bm{X}}{\bm{L}}^{\top}\right).$
(15)

More important properties include (see, e.g., golub13 ): (i)
$({\bm{M}}\otimes{\bm{N}})^{\top}={\bm{M}}^{\top}\otimes{\bm{N}}^{\top}$; (ii)
$({\bm{M}}_{1}\otimes{\bm{N}}_{1})({\bm{M}}_{2}\otimes{\bm{N}}_{2})=({\bm{M}}_{1}{\bm{M}}_{2}\otimes{\bm{N}}_{1}{\bm{N}}_{2})$;
and (iii) $\|{\bm{M}}\otimes{\bm{N}}\|_{2}=\|{\bm{M}}\|_{2}\|{\bm{N}}\|_{2}$.

###  3.2 POD-DEIM reduced dynamics

Consider the nonlinear dynamical system (1). In what follows we will assume
without loss of generality that the linear and nonlinear terms on the right-
hand side can be explicitly separated to yield a semilinear system of the form

$\left\\{\begin{array}[]{l}\dot{\bm{y}}(s)=f({\bm{y}}(s),u(s),s):={\bm{L}}{\bm{y}}(s)+{\bm{f}}({\bm{y}}(s),u(s),s),\;\;s\in(t,T],\\\
y(t)=x\in\mathbb{R}^{N},\end{array}\right.$ (16)

where ${\bm{L}}\in\mathbb{R}^{N\times N}$ and
${\bm{f}}:\mathbb{R}^{N}\times\mathbb{R}^{c}\times[t,T]\rightarrow\mathbb{R}^{N}$
is a continuous function in all arguments and locally Lipschitz-type with
respect to the first variable.

In Alla_Saluzzi_2020  the authors reduce the number of variables involved in
the system (16) by means of POD-DEIM. More precisely, consider a set of
solution snapshots $y_{i}=y(s_{i})$ collected at $n_{s}$ time instances
$s_{i}$ in the timespan $(t,T]$, and consider the snapshot matrix

${\bm{S}}=[y_{1},y_{2},\ldots,y_{n_{s}}]\in\mathbb{R}^{N\times
n_{s}},\qquad{\cal S}=\mbox{Range}({\bm{S}}).$

A POD basis of dimension $k\leq n_{s}$ is obtained by orthogonal reduction of
the matrix ${\bm{S}}$. That is, given the Singular Value Decomposition (SVD)

${\bm{S}}={\bm{V}}{\bm{\Sigma}}{\bm{W}}^{\top},\qquad{\bm{V}},{\bm{W}}\in\mathbb{R}^{N\times
n_{s}},{\bm{\Sigma}}\in\mathbb{R}^{n_{s}\times n_{s}},$

the POD basis is given by $\\{{\bm{v}}_{1},\ldots,{\bm{v}}_{k}\\}$, where
${\bm{V}}_{k}=[{\bm{v}}_{1},\ldots,{\bm{v}}_{k}]\in\mathbb{R}^{N\times k}$ is
the matrix of truncated left singular vectors related to the $k$ largest
singular values contained on the diagonal of ${\bm{\Sigma}}$.

Given the matrix ${\bm{V}}_{k}$, the state vector ${\bm{y}}(s)$ can be
approximated as ${\bm{y}}(s)\approx{\bm{V}}_{k}\widehat{\bm{y}}(s)$, for all
$s\in(t,T]$, where $\widehat{\bm{y}}(s)\in\mathbb{R}^{k}$ solves the reduced
dynamical system

$\left\\{\begin{array}[]{l}\dot{\widehat{y}}(s)={\bm{V}}^{\top}_{k}f({\bm{V}}_{k}\widehat{\bm{y}}(s),u(s),s):={\bm{V}}_{k}^{\top}{\bm{L}}{\bm{V}}_{k}\widehat{\bm{y}}(s)+{\bm{V}}_{k}^{\top}{\bm{f}}({\bm{V}}_{k}\widehat{\bm{y}}(s),u(s),s),\\\
\widehat{y}(t)={\bm{V}}_{k}^{\top}x.\end{array}\right.$ (17)

To ensure that the reduced model can be simulated with a computational cost
independent of $N$, we need to avoid lifting the nonlinear term before
projection onto the low-dimensional space. Consequently, the Discrete
Empirical Interpolation Method (DEIM) from chaturantabut2010nonlinear  is used
to interpolate the nonlinear function.

To this end we consider an approximation of the form

${\bm{f}}({\bm{V}}_{k}\widehat{\bm{y}}(s),u(s),s)\approx{\bm{\varPhi}}_{p}\widehat{\bm{f}}({\bm{V}}_{k}\widehat{\bm{y}}(s),u(s),s),\qquad\widehat{\bm{f}}({\bm{V}}_{k}\widehat{\bm{y}}(s),u(s),s)\in\mathbb{R}^{p},$

where
${\bm{\varPhi}}_{p}=[{\bm{\varphi}}_{1},\ldots,{\bm{\varphi}}_{p}]\in\mathbb{R}^{N\times
p}$, with $p\ll N$, and $\\{{\bm{\varphi}}_{1},\ldots,{\bm{\varphi}}_{p}\\}$
is a POD basis of dimension $p$ obtained from the set of snapshots
$\\{{\bm{f}}({y}_{i},u(s_{i}),s_{i})\\}_{i=1}^{n_{s}}$. The overdetermined
system is solved by interpolation, ensuring that the left and right side of
the equation is equal at $p$ selected points. That is, given the matrix
${\bm{P}}=[{\bm{e}}_{\rho_{1}},\ldots,{\bm{e}}_{\rho_{p}}]\in\mathbb{R}^{p}$
containing a subset of columns of the identity matrix, we ensure that
${\bm{P}}^{\top}{\bm{f}}({\bm{V}}_{k}\widehat{\bm{y}}(s),u(s),s)={\bm{P}}^{\top}{\bm{\varPhi}}_{p}\widehat{\bm{f}}({\bm{V}}_{k}\widehat{\bm{y}}(s),u(s),s)$,
so that

${\bm{f}}({\bm{V}}_{k}\widehat{\bm{y}}(s),u(s),s)\approx\widetilde{\bm{f}}(\widehat{\bm{y}}(s),u(s),s):={\bm{\varPhi}}_{p}({\bm{P}}^{\top}{\bm{\varPhi}}_{p})^{-1}{\bm{P}}^{\top}{\bm{f}}({\bm{V}}_{k}\widehat{\bm{y}}(s),u(s),s).$

Throughout this paper we deal with nonlinear functions that are evaluated
element-wise, so that

$\widetilde{\bm{f}}(\widehat{\bm{y}}(s),u(s),s)={\bm{\varPhi}}_{p}({\bm{P}}^{\top}{\bm{\varPhi}}_{p})^{-1}{\bm{f}}({\bm{P}}^{\top}{\bm{V}}_{k}\widehat{\bm{y}}(s),u(s),s),$

and the nonlinear term is only evaluated at $p$ entries.

###  3.3 A Multilinear HJB-POD-DEIM algorithm on a tree structure

In this section we illustrate how, under certain hypotheses, the discrete
system (16) can be expressed, integrated and reduced in terms of multilinear
arrays; see e.g., Simoncini2017 ; Autilia2019matri ; palitta2016 ; kirsten.22
. We focus specifically on the case where the discrete system (16) stems from
the space discretization of a semilinear PDE of the form

$\left\\{\begin{array}[]{l}\small\partial_{s}y(s,x)={\cal
L}\left(y(s,x)\right)+{\bm{f}}(\nabla{{y}(s,x)},{y}(s,x),u(s),s),\\\
\qquad\,\,\,s\in(t,T],x\in\Omega,\\\
y(t,x)=\tilde{y}(x),\;\;x\in\Omega,\end{array}\right.$ (18)

with ${\cal L}$ a linear differential operator, ${\bm{f}}$ a generic nonlinear
operator and $\Omega\subset\mathbb{R}^{d}$, for $d=2,3$.

####  3.3.1 Discretization in terms of multilinear arrays

Consider the operator $\mathcal{L}$ to be a second order differential operator
with separable coefficients. Then, the physical domain can be mapped to a
hypercubic domain
${\bm{\Omega}}=[a_{1},b_{1}]\times\cdots\times[a_{d},b_{d}]$, if the operator
is discretized via a tensor basis. Examples of such discretizations include,
but are not limited to, spectral methods, and finite differences on
parallelepipedal domains. See, e.g., Simoncini2017 ; palitta2016 ; kirsten.22
for more information regarding the assumptions on the operators, domains and
discretization techniques. Here we consider $\mathcal{L}$ as a $d-$dimensional
Laplace operator for illustration purposes, but more general operators can
also be treated. Under these conditions, it holds (from (16)) that

${\bm{L}}=\sum_{m=1}^{d}{\bm{I}}_{n_{d}}\otimes\cdots\otimes\overset{m}{\bm{A}}_{m}\otimes\cdots\otimes{\bm{I}}_{n_{1}}\in\mathbb{R}^{N\times
N},$

where ${\bm{A}}_{m}\in\mathbb{R}^{n_{m}\times n_{m}}$ contains the
approximation of the second derivative in the $x_{m}$ direction. We will also
consider problems where the nonlinear term $\bm{f}$ depends on the first
derivative of the state vector, so that we also define the matrix

${\bm{D}}=\sum_{m=1}^{d}{\bm{I}}_{n_{d}}\otimes\cdots\otimes\overset{m}{\bm{B}}_{m}\otimes\cdots\otimes{\bm{I}}_{n_{1}}\in\mathbb{R}^{N\times
N},$

where ${\bm{B}}_{m}\in\mathbb{R}^{n_{m}\times n_{m}}$ contains the
approximation of the first derivative in the $x_{m}$ direction. The vectors
${\bm{y}}(t)\in\mathbb{R}^{N}$ from (16) then represent the vectorization of
the elements of a tensor
$\boldsymbol{\mathcal{Y}}(t)\in\mathbb{R}^{n_{1}\times\cdots\times n_{d}}$,
such that ${\bm{y}}(t)={\tt vec}(\boldsymbol{\mathcal{Y}}(t))$,
${\bm{L}}{\bm{y}}=\mbox{vec}\left(\mathcal{A}(\boldsymbol{\mathcal{Y}})\right)$
and
${\bm{D}}{\bm{y}}=\mbox{vec}\left(\mathcal{D}(\boldsymbol{\mathcal{Y}})\right)$,
where111For the case $d=2$, (19) are Sylvester operators of the form
${\bm{A}}_{1}{\bm{Y}}+{\bm{Y}}{\bm{A}}_{2}^{\top}$ and
${\bm{B}}_{1}{\bm{Y}}+{\bm{Y}}{\bm{B}}_{2}^{\top}$ respectively Simoncini2017
.

$\mathcal{A}(\boldsymbol{\mathcal{Y}}):=\sum_{m=1}^{d}\,\boldsymbol{\mathcal{Y}}\times_{m}{\bm{A}}_{m}\quad\mbox{and}\quad\mathcal{D}(\boldsymbol{\mathcal{Y}}):=\sum_{m=1}^{d}\,\boldsymbol{\mathcal{Y}}\times_{m}{\bm{B}}_{m}.$
(19)

Moreover, if the function $\mathcal{F}:{\cal
S}\times[0,t_{f}]\rightarrow\mathbb{R}^{n_{1}\times\cdots\times n_{d}}$
represents the function ${\bm{f}}$ evaluated at the entries of the array
$\boldsymbol{\mathcal{Y}}$ and $\mathcal{D}(\boldsymbol{\mathcal{Y}})$, then
it holds that ${\bm{f}}({\bm{D}}{\bm{y}},{\bm{y}}(s),u(s),s)={\tt
vec}\left(\mathcal{F}\left(\mathcal{D}(\boldsymbol{\mathcal{Y}}),\boldsymbol{\mathcal{Y}},u(s),s\right)\right)$,
and (16) can be written in the form

$\begin{cases}\dot{\boldsymbol{\mathcal{Y}}}(s)&=\mathcal{A}(\boldsymbol{\mathcal{Y}}(s))+\mathcal{F}\left(\mathcal{D}(\boldsymbol{\mathcal{Y}}(s)),\boldsymbol{\mathcal{Y}}(s),u(s),s\right),\\\
{\boldsymbol{\mathcal{Y}}}(t)&=\boldsymbol{\mathcal{X}}\in\mathbb{R}^{n_{1}\times\cdots\times
n_{d}}.\\\ \end{cases}$ (20)

The boundary conditions are contained in the matrices ${\bm{A}}_{m}$ and
${\bm{B}}_{m}$, $m=1,\ldots,d$; see e.g., Autilia2019matri ; palitta2016 .
From here on forward we consider the case where $n_{1}=\cdots=n_{d}=n$, so
that $N=n^{d}$.

####  3.3.2 Higher-Order POD (HO-POD) model reduction

As it has been shown in kirsten.22 , great savings in terms of memory
requirements and computational time can be obtained by applying model order
reduction directly to the system (20) instead of first vectorizing and
applying model reduction to the vectorized system (16). To this end, we
consider an approximation of the form

$\boldsymbol{\mathcal{Y}}(s)\approx\widetilde{\boldsymbol{\mathcal{Y}}}(s):=\widehat{\boldsymbol{\mathcal{Y}}}(s)\bigtimes_{m=1}^{d}{\bm{V}}_{m},$

where ${\bm{V}}_{m}\in\mathbb{R}^{n_{m}\times k_{m}}$ are tall matrices with
orthonormal columns and
$\widehat{\boldsymbol{\mathcal{Y}}}(s)\in\mathbb{R}^{k_{1}\times\cdots\times
k_{d}}$ $(k_{m}\ll n)$ satisfies the low-dimensional equation

$\begin{cases}\dot{\widehat{\boldsymbol{\mathcal{Y}}}}(s)&=\widehat{\mathcal{A}}(\widehat{\boldsymbol{\mathcal{Y}}}(s))+\widehat{\mathcal{F}}\left(\widehat{\mathcal{D}}(\widehat{\boldsymbol{\mathcal{Y}}}(s)),\widehat{\boldsymbol{\mathcal{Y}}}(s),u(s),s\right),\\\
\widehat{\boldsymbol{\mathcal{Y}}}(t)&=\boldsymbol{\mathcal{X}}\bigtimes_{m=1}^{d}{\bm{V}}_{m}^{\top}\in\mathbb{R}^{k_{1}\times\cdots\times
k_{d}},\\\ \end{cases}$ (21)

where

$\widehat{\mathcal{F}}\left(\widehat{\mathcal{D}}(\widehat{\boldsymbol{\mathcal{Y}}}(s)),\widehat{\boldsymbol{\mathcal{Y}}}(s),u(s),s\right)={\mathcal{F}}\left({\mathcal{D}}(\widetilde{\boldsymbol{\mathcal{Y}}}(s)),\widetilde{\boldsymbol{\mathcal{Y}}}(s),u(s),s\right)\bigtimes_{m=1}^{d}{\bm{V}}_{m}^{\top}$
(22)

and

$\widehat{\mathcal{A}}(\widehat{\boldsymbol{\mathcal{Y}}}):=\sum_{m=1}^{d}\,\widehat{\boldsymbol{\mathcal{Y}}}\times_{m}\widehat{\bm{A}}_{m},\quad\widehat{\bm{A}}_{m}={\bm{V}}_{m}^{\top}{\bm{A}}_{m}{\bm{V}}_{m}.$
(23)

The matrices ${\bm{V}}_{m}\in\mathbb{R}^{n_{m}\times k_{m}}$ can be obtained
via the HO-POD algorithm described in kirsten.22 . That is, given a set of
snapshots $\\{\boldsymbol{\mathcal{Y}}(s_{i})\\}_{i=1}^{n_{s}}$, each matrix
${\bf V}_{m}$ is constructed in order to approximate the left range space of
the matrix

$\boldsymbol{\mathcal{S}}_{(m)}=\begin{pmatrix}\boldsymbol{\mathcal{Y}}_{(m)}(s_{1}),\ldots,\boldsymbol{\mathcal{Y}}_{(m)}(s_{n_{s}})\end{pmatrix}\in\mathbb{R}^{n\times
n^{d-1}n_{s}},\quad\mbox{for}\quad m=1,\ldots,d,$

where $m$ represents the mode along which the tensor is unfolded, and
$\boldsymbol{\mathcal{S}}\in\mathbb{R}^{n\times\cdots\times n}$ is a tensor of
order $d$ containing the snapshots. Note that neither the matrices
$\boldsymbol{\mathcal{S}}_{(m)}$ or the tensor $\boldsymbol{\mathcal{S}}$ is
ever explicitly constructed or stored. Instead we follow the dynamic algorithm
initially introduced in Kirsten.Simoncini.arxiv2020  for approximating the
left range space of $\boldsymbol{\mathcal{S}}_{(m)}$. In kirsten.22  a simpler
algorithm was used to construct the approximation space in the tensor setting.
Here we implement the more refined dynamic algorithm for the tensor setting;
the inclusion of snapshot information is discussed here, whereas snapshot
selection will be presented in section 4.

Suppose $\kappa$222We refer the reader to Kirsten.Simoncini.arxiv2020  for a
detailed experimental analysis on the role of the parameter $\kappa$. is the
maximum admissible dimension for the reduced space in all modes, selected
a-priori, and consider the initial condition $\boldsymbol{\mathcal{Y}}(t)$.
Let

$\boldsymbol{\mathcal{Y}}(t)\approx\boldsymbol{\mathcal{C}}(t)\bigtimes_{m=1}^{d}{\bm{U}}_{m}^{(0)},$

represent the sequentially truncated higher order SVD333For the case $d=2$,
however, we just use the standard MATLAB SVD function. (STHOSVD)
vannieuwenhoven2012  of $\boldsymbol{\mathcal{Y}}(t)$, where
${\bm{U}}_{m}^{(0)}$ contains the first $\kappa$ dominant left singular
vectors of $\boldsymbol{\mathcal{Y}}_{(m)}(t)$. For each mode these left
singular vectors are collected into the matrix
$\widetilde{\bm{V}}_{m}={\bm{U}}_{m}^{(0)}$, $m=1,\ldots,d$.

Subsequently, suppose the snapshot at time instance $s_{j}$ has been selected
for inclusion into the approximation space and let
$\boldsymbol{\mathcal{Y}}(s_{j})\approx\boldsymbol{\mathcal{C}}(s_{j})\bigtimes_{m=1}^{d}{\bm{U}}_{m}^{(j)},$
represent the STHOSVD of the selected snapshot and let
${\bm{\Sigma}}_{m}^{(j)}$ contain the first $\kappa$ singular values of
$\boldsymbol{\mathcal{Y}}_{(m)}(s_{j})$ on the main diagonal. The
approximation spaces are updated by appending the new singular values and
vectors, so that

$\widetilde{\bm{V}}_{m}\leftarrow[\widetilde{\bm{V}}_{m},{\bm{U}}_{m}^{(j)}],\quad{\rm
and}\quad\widetilde{\bm{\Sigma}}_{m}\leftarrow{\tt
blkdiag}(\widetilde{\bm{\Sigma}}_{m},{\bm{\Sigma}}_{m}^{(j)}).$

Eventually the diagonal entries of $\widetilde{\bm{\Sigma}}_{m}$ are reordered
decreasingly and truncated so that the largest $\kappa$ values are retained,
with the vectors in $\widetilde{\bm{V}}_{m}$ reordered and truncated
accordingly.

At the end of the procedure, when all snapshots have been processed, the final
basis vectors are obtained by orthogonal reduction of the matrices
$\widetilde{\bm{V}}_{m}$. More precisely, let
$\widetilde{\bm{V}}_{m}=\overline{\bm{V}}_{m}\,\overline{\bm{\Sigma}}_{m}\,\overline{\bm{W}}_{m}^{\top}$
be the SVD of $\widetilde{\bm{V}}_{m}$. The final basis matrices
${\bm{V}}_{m}$ are obtained by truncating the first $k_{m}$ dominant singular
vectors of $\overline{\bm{V}}_{m}$ according to the criterion

${\sqrt{\sum_{i=k_{m}+1}^{\kappa}({\sigma}_{m}^{(i)})^{2}}}<\tau\,{\sqrt{\sum_{i=1}^{\kappa}({\sigma}_{m}^{(i)})^{2}}},$
(24)

for some $\tau\in(0,1)$, where ${\sigma}_{m}^{(i)}$ is the $i$-th diagonal
element of $\overline{\bm{\Sigma}}_{m}$.

####  3.3.3 Higher-Order DEIM (HO-DEIM) approximation of the nonlinear term

It is clear from (22) that a bottleneck forms around the reduced nonlinear
term, similar to the vector setting. To this end, we consider the HO-DEIM
algorithm from kirsten.22  to circumvent the issue. Consequently, suppose the
tall matrices ${\bm{\Phi}}_{m}\in\mathbb{R}^{n_{m}\times p_{m}}$, for
$m=1,2,\ldots,d$, have been constructed as the output of the HO-POD method
described above for the snapshots
$\\{\mathcal{F}\left(\mathcal{D}(\boldsymbol{\mathcal{Y}}(s_{i})),\boldsymbol{\mathcal{Y}}(s_{i}),u(s_{i}),s_{i}\right)\\}_{i=1}^{n_{s}}$.
Furthermore, consider $d$ matrices ${\bm{P}}_{m}\in\mathbb{R}^{n_{m}\times
p_{m}}$ each containing a subset of columns of the $n_{m}\times n_{m}$
identity matrix. The matrices ${\bm{P}}_{m}$ are each respectively obtained as
the output of the q-deim algorithm gugercin2018  with input ${\bm{\Phi}}_{m}$.
The ho-deim approximation of 22 is then given by

$\begin{split}\widehat{\mathcal{F}}\left(\widehat{\mathcal{D}}(\widehat{\boldsymbol{\mathcal{Y}}}(s)),\widehat{\boldsymbol{\mathcal{Y}}}(s),u(s),s\right)&\approx{\mathcal{F}}\left({\mathcal{D}}(\widetilde{\boldsymbol{\mathcal{Y}}}(s)),\widetilde{\boldsymbol{\mathcal{Y}}}(s),u(s),s\right)\bigtimes_{m=1}^{d}{\bf
F}_{m}\\\ &=\widehat{\mathcal{F}}^{\mbox{\tiny\sc
deim}}\left(\widehat{\mathcal{D}}(\widehat{\boldsymbol{\mathcal{Y}}}(s)),\widehat{\boldsymbol{\mathcal{Y}}}(s),u(s),s\right),\end{split}$
(25)

where

${\bf F}_{m}={\bf V}_{m}^{\top}{\bf\Phi}_{m}({\bf
P}_{m}^{\top}{\bf\Phi}_{m})^{-1}{\bf P}_{m}^{\top}.$

If $\mathcal{F}$ is evaluated element-wise at the components of
$\widetilde{\boldsymbol{\mathcal{Y}}}(s)$ and
${\mathcal{D}}(\widetilde{\boldsymbol{\mathcal{Y}}}(s))$, then it holds that

$\begin{split}\savestack{\tmpbox}{\stretchto{\scaleto{\scalerel*[width("{\mathcal{F}}\left({\mathcal{D}}(\widetilde{\boldsymbol{\mathcal{Y}}}(s)),\widetilde{\boldsymbol{\mathcal{Y}}}(s),u(s),s\right)")]{\kern-0.6pt\bigwedge\kern-0.6pt}{\rule[-505.89pt]{4.30554pt}{505.89pt}}}{}}{0.5ex}}\ensurestackMath{\stackon[1pt]{{\mathcal{F}}\left({\mathcal{D}}(\widetilde{\boldsymbol{\mathcal{Y}}}(s)),\widetilde{\boldsymbol{\mathcal{Y}}}(s),u(s),s\right)}{\tmpbox}}:=&{\mathcal{F}}\left({\mathcal{D}}(\widetilde{\boldsymbol{\mathcal{Y}}}(s)),\widetilde{\boldsymbol{\mathcal{Y}}}(s),u(s),s\right)\bigtimes_{m=1}^{d}{\bf
P}_{m}^{\top}\\\
=&{\mathcal{F}}\left({\mathcal{D}}(\widetilde{\boldsymbol{\mathcal{Y}}}(s))\bigtimes_{m=1}^{d}{\bf
P}_{m}^{\top},\,\,\widetilde{\boldsymbol{\mathcal{Y}}}(s)\bigtimes_{m=1}^{d}{\bf
P}_{m}^{\top},\,u(s),s\right).\end{split}$ (26)

In this case the nonlinear term is evaluated at only $p_{1}p_{2}\cdots p_{m}$
entries.

####  3.3.4 The reduced optimal control problem on a tree structure

In this section we explore how the full procedure combining the tree structure
algorithm and the HO-POD-DEIM Model reduction technique is split into an
offline and online stage to solve the HJB equation (8) and determine the
optimal control (9).

##### Offline Stage

The offline stage consists of two important steps, namely snapshot collection
and basis construction. To this end, we select a coarse time step
$\widehat{\Delta t}$ and control set $\widehat{U}$. The basis is constructed
on the fly following Sections 3.3.2-3.3.3, on the nodes of the tree, which is
constructed following Section 2.1. This is a computationally expensive step,
as the full-dimensional space is explored in this phase. To this end, we
discuss a collection of nuances related to the implementation in Section 4.

##### Online Stage

At this stage, the computed basis vectors are exploited to construct a reduced
dimensional tree, approximate the reduced value function and compute the
optimal trajectory at a fraction of the inital cost.

  * •

Construction of the reduced tree. Here we fix the desired wider discrete
control set $\widetilde{U}\subset\widehat{U}$ and/or a smaller time step
$\Delta t\leq\widehat{\Delta t}$ for the resolution of the HO-POD-DEIM reduced
dynamical system

$\begin{cases}\dot{\widehat{\boldsymbol{\mathcal{Y}}}}(s)&=\widehat{\mathcal{A}}(\widehat{\boldsymbol{\mathcal{Y}}}(s))+\widehat{\mathcal{F}}^{\mbox{\tiny\sc
deim}}\left(\widehat{\mathcal{D}}(\widehat{\boldsymbol{\mathcal{Y}}}(s)),\widehat{\boldsymbol{\mathcal{Y}}}(s),u(s),s\right),\\\
\widehat{\boldsymbol{\mathcal{Y}}}(t)&=\boldsymbol{\mathcal{X}}\bigtimes_{m=1}^{d}{\bm{V}}_{m}^{\top}\in\mathbb{R}^{k_{1}\times\cdots\times
k_{d}}.\\\ \end{cases}$ (27)

Following Section 2.1, we build the reduced tree $\widehat{\mathcal{T}}$ as
done for the offline stage. The cardinality of tree, however, still grows
exponentially, despite the reduced dimension of the dynamical system. As a
result we analyze a collection of important pruning criteria in Section 5.1 in
an attempt to reduce the cardinality of the tree dynamically during the
construction.

  * •

Approximation of the reduced value function. The value function computed in
the reduced space will be denoted by $\widehat{V}(\widehat{x},t)$ and its
approximation at time $t_{n}$ as $\widehat{V}^{n}(\widehat{x})$. Its
resolution will follow the classical scheme introduced in Section 2.1 :

$\begin{cases}\widehat{V}^{n}(\widehat{\zeta}^{n}_{i})=\min\limits_{u\in\widetilde{U}}\\{\widehat{V}^{n+1}(T^{n+1}(\widehat{\zeta}^{n}_{i},u)+\Delta
t\,\widehat{L}(\widehat{\zeta}^{n}_{i},u,t_{n})\\},\\\
\qquad\qquad\qquad\qquad\qquad\widehat{\zeta}^{n}_{i}\in\widehat{\mathcal{T}}^{n}\,,n=N_{t}-1,\ldots,0,\\\
\widehat{V}^{N_{t}}(\widehat{\zeta}^{N_{t}}_{i})=\widehat{g}(\widehat{\zeta}_{i}^{N_{t}}),\qquad\qquad\qquad\qquad\qquad\qquad\widehat{\zeta}_{i}^{N_{t}}\in\widehat{\mathcal{T}}^{N_{t}},\end{cases}$
(28)

where

$\widehat{L}(\widehat{\zeta},u,t)=L(\widehat{\zeta}\bigtimes_{m=1}^{d}{\bm{V}}_{m}^{\top},u,t),$
$\widehat{g}(\widehat{\zeta})=g(\widehat{\zeta}_{i}^{N_{t}}\bigtimes_{m=1}^{d}{\bm{V}}_{m}^{\top})$

and $T^{n+1}(\widehat{\zeta},u)$ stands for the time evolution of the node
$\widehat{\zeta}$ with control $u$ at time $t_{n+1}$.

  * •

Computation of the optimal trajectory.

The optimal trajectory can be seen as a specific path in the tree structure.
For this reason during the computation of the value function we store the
minimizing indices in (28). Once completed the computation of the value
function, the optimal path will be given just following the tree branches
which returns the minimum index.

##  4 Hints for the implementation

In both the offline and online phases of the procedure, great savings in terms
of CPU time and memory requirements can be obtained if implemented in an
efficient way. Consequently, in this section we discuss how the snapshots are
selected and how the simulated tree nodes can be efficiently stored to save on
memory requirements in the offline phase. Moreover, we discuss how the reduced
model can be efficiently simulated at many time steps and control inputs, to
avoid high computational costs in the online phase.

###  4.1 Snapshot selection

The full order model is simulated on a coarse timegrid with two control inputs
which are the two extremes of the control domain, as discussed in
Alla_Saluzzi_2020 . To avoid excessive computational work we only include
information from snapshots that are not yet well approximated in the current
basis. The condition for snapshot inclusion is given by the projection error,
that is:

$\mbox{if:}\quad\frac{\|{\boldsymbol{\cal Y}}(s_{i})-{\boldsymbol{\cal
Y}}(s_{i})\bigtimes_{m=1}^{d}{\bm{V}}_{m}{\bm{V}}_{m}^{\top}\|_{F}}{\|{\boldsymbol{\cal
Y}}(s_{i})\|_{F}}>\tau\quad\mbox{then:}\quad\mbox{\bf Include},$ (29)

for some $\tau\in(0,1)$. Here the matrices ${\bm{V}}_{m}$ contain the current
basis vectors in all modes, updated dynamically with snapshot information from
the previous selected snapshots as described in section 3.3.2.

###  4.2 Efficient memory allocation by low-rank storage of tree nodes

One challenge of the method presented in Alla_Saluzzi_2020  is in terms of
memory in the offline phase, since the high fidelity solutions need to be
calculated and stored for several time steps and control inputs in order to
form a reduced order model. In particular, the nodes of the tree $\mathcal{T}$
are vectorized and stored in a matrix
${\bm{T}}\in\mathbb{R}^{N\times|\mathcal{T}|}$, where
$|\mathcal{T}|=O(M^{N_{t}})$, with $M$ fixed as the number of control inputs
and $N_{t}$ as the number of time steps. The exponential growth of the second
dimension greatly limits the number of snapshots that can be stored.

In this paper, we suggest the following improvement. Since the full order
model is simulated in array form, the snapshots at the resulting tree nodes
are either matrices or tensors. Consequently, we can take advantage of the
(possible) low-rank structure of each node. That is, we compute the STHOSVD of
each computed nodal value, truncated to the first $\kappa$ singular vectors in
each mode, so that
$\boldsymbol{\mathcal{Y}}(s_{j})\approx\boldsymbol{\mathcal{C}}(s_{j})\bigtimes_{m=1}^{d}{\bm{U}}_{m}^{(j)}$,
with $\kappa$ selected a-priori as discussed in section 3.3.2. As a result,
the node can be stored in low-rank form to be recalled for later computations.
More precisely, we collect and store only the dominant singular vectors in
each mode and the low-dimensional core tensors such that

$\overline{\bm{U}}_{m}\leftarrow[\overline{\bm{U}}_{m},{\bm{U}}_{m}^{(j)}]\quad{\rm
for}\,\,\,m=1,\ldots,d\,\,\,{\rm
and}\quad\overline{\bm{c}}\leftarrow[\overline{\bm{c}},{\rm
vec}(\boldsymbol{\mathcal{C}}(s_{j}))].$

When required at the next time level, the snapshots can easily be computed
from its Tucker decomposition (kolda2009, , Section 4). This process allows us
to store vectors of length $n_{1},\ldots,n_{d}$ instead of $N=n_{1}n_{2}\cdots
n_{d}$. The number of vectors stored depends on the rank of the considered
snapshots. A further advantage is that when a snapshot is selected for
inclusion into the approximation space, a HOSVD is required as discussed in
section 3.3.2, which will be readily available thanks to this procedure.

Furthermore, it has been observed that only the nodes from the previous level
of the tree need to be stored, since the snapshots from the earlier levels are
automatically processed and discarded during the HO-POD basis construction.
Finally, we observe that the computation of the value function does not
require the knowledge of the nodes, but only of the corresponding cost
evaluation. In this way we are going to store only the corresponding scalar
cost and the nodes will be erased after the computation of its tree sons.

###  4.3 Efficient simulation of the reduced model (27)

An important ingredient in the success of the HO-POD-DEIM reduction procedure
is the ability to integrate the reduced model (27) in array form without
vectorization. In this paper we consider the semi-implicit Euler scheme, given
that the considered model is typically associated with a stiff linear term and
a nonstiff nonlinear term, but several alternatives can be considered
Autilia2019matri ; kirsten.22 . More precisely, suppose
${\widehat{\boldsymbol{\cal Y}}}^{(j)}$ is an approximation of
${\widehat{\boldsymbol{\cal Y}}}(s_{j})$, then the linear system

$(\widehat{{\cal I}}-\Delta t\widehat{\cal{A}}){\widehat{\boldsymbol{\cal
Y}}}^{(j)}={\widehat{\boldsymbol{\cal Y}}}^{(j-1)}+\Delta
t\widehat{\mathcal{F}}^{\mbox{\tiny\sc
deim}}\left(\widehat{\mathcal{D}}({\widehat{\boldsymbol{\cal
Y}}}^{(j-1)}),{\widehat{\boldsymbol{\cal
Y}}}^{(j-1)},u(s_{j-1}),s_{j-1}\right)$ (30)

needs to evaluated at each time level $s_{j}$. Once again, vectorizing the
linear system (30) at each time step will reduce the computational gains
related to the reduction in Array form. Instead, (30) can be solved in array
form using the direct method presented in simoncini.2020.boll ; kirsten.22 .

##  5 Pruning techniques

Although theoretically the tree structure enables to compute the solution for
arbitrary high dimensional problems, since we are not restricted to the direct
discretization of a domain, its construction turns to be computationally
expensive, due to the exponential growth of its cardinality, For this reason
in this section we are going to introduce and analyse different $pruning$
$criteria$ able to reduce the growth of the tree, but keeping the same
accuracy.

###  5.1 Geometric pruning

A pruning criterion based on a comparison of the nodes in euclidean norm has
been introduced in saluzzi2022error . More precisely, two given nodes
$\zeta^{n}_{i}$ and $\zeta^{n}_{j}$ will be merged if

$\|\zeta^{n}_{i}-\zeta^{n}_{j}\|\leq\epsilon,\quad\mbox{ with }i\neq j\mbox{
and }n=0,\ldots,N_{t},$ (31)

for a given threshold $\epsilon>0$. To ensure first order convergence, the
threshold $\epsilon$ must scale as $\Delta t^{2}$ (we refer to
saluzzi2022error  for more details about the error estimates). This pruning
criterion has been successfully applied to low and high dimensional problems,
but the main drawback is the expensive computation of distances in high
dimension. One possible solution relies on the projection of the data onto a
lower dimension minimizing the variance of the data. This procedure is already
encoded in the above described algorithm, since we are reducing the dimension
of the problem keeping the main features.

###  5.2 Statistical pruning

In this section we introduce a new iterative pruning criterion based on
statistical information about the value function. We suppose we are starting
with a certain control set $U_{1}$. First, we construct the tree
$\mathcal{T}_{1}$ based on the control set $U_{1}$ and the value function
computed on the tree will be denoted by $V_{1}(x,t)$. Afterwards, we refine
the constructed tree based on the information on the value function: fixing a
ratio $\rho\in(0,1]$ of the nodes, we retain just those with the lowest value
function, obtaining a new tree $\widetilde{\mathcal{T}}_{1}$. More precisely,
we have that $|\widetilde{\mathcal{T}}_{1}|=\rho|{\mathcal{T}}_{1}|$ and for
every time level $t_{n}\in[t,T]$ and every node $\zeta\in\mathcal{T}_{1}$
there exists a node $\widetilde{\zeta}\in\widetilde{\mathcal{T}}_{1}$ such
that $V_{1}(\widetilde{\zeta},t_{n})\leq V_{1}(\zeta,t_{n})$. Hence, we can
start with the construction of a new tree $\mathcal{T}_{2}$ with a wider
control set $U_{2}\supset U_{1}$ such that the nodes are constrained in the
zones where the previous value function had the lowest values, $i.e.$

$\min_{\widetilde{\mathcal{T}}^{n}_{1}}\widetilde{\zeta}\leq\zeta^{n}_{j}\leq\max_{\widetilde{\mathcal{T}}^{n}_{1}}\widetilde{\zeta},\quad\forall\zeta^{n}_{j}\in\mathcal{T}^{n}_{2},\,n\in\\{N_{start},\ldots,N_{t}\\}$

where the minimum and the maximum are computed element-wise, as well the
inequalities. The statistical pruning is applied starting from an arbitrary
time $t_{N_{start}}$ since the first levels contain few nodes. We usually will
fix $N_{start}=3$. The entire procedure can be iterated doubling the number of
controls in each step. Computed the tree $\widetilde{\mathcal{T}}_{k}$ at the
$k$-th iteration, the subsequent tree $\mathcal{T}_{k+1}$ will satisfy the
constraint

$\min_{\widetilde{\mathcal{T}}^{n}_{k}}\widetilde{\zeta}\leq\zeta^{n}_{j}\leq\max_{\widetilde{\mathcal{T}}^{n}_{k}}\widetilde{\zeta},\quad\forall\zeta^{n}_{j}\in\mathcal{T}^{n}_{k+1},\,n\in\\{N_{start},\ldots,N_{t}\\}.$
(32)

Since we neglect the nodes which do not satisfy (32), the problem can be
regarded as a state-constrained problem where the constraint is given by the
relation (32). We refer to afs20  for more details about the coupling of the
tree with state-constrained problems. The ratio $\rho$ is fixed such that it
still retains the optimal trajectory from the previous iteration. In this way
we can ensure that the value function is not increasing during the iterative
procedure. Therefore, we denote by $V^{n}_{k}(x)$ the value function obtained
in the $k$-th iteration at the point $x$ at time $t_{n}$. By construction, we
can notice that the iterative value function at the initial time is non
increasing, $i.e.$

$V^{0}_{k+1}(x)\leq V^{0}_{k}(x),\quad\forall k\geq 0$

and bounded from below since

$V^{n}_{k}(x)\geq-TM_{f}-M_{g},\quad\forall k\geq
0,\,n\in\\{0,\ldots,N_{t}\\},\,x\in\mathcal{T}^{n}_{k},$

using the hypothesis (3). Hence, the iterative scheme is convergent and it can
repeated until we reach a maximum number of iterations or it satisfies a
stopping criterion. In Algorithm 1 the method is sketched.

Algorithm 1 Statistical pruning

1:Choose an initial condition $x_{0}$, a ratio $\rho$, a starting time
$t_{N_{start}}$, a tolerance $tol$, a maximum number of iteration $k_{max}$
and a initial number of discrete controls $M$

2:Build a tree $\mathcal{T}_{1}$ with $M$ controls and compute the value
function $V^{n}_{1}(x)$

3:while $res>tol$ and $k\leq k_{max}$ do

4: $M:=2M-1$

5: Construct $\widetilde{\mathcal{T}}_{k}$ retaining a ratio $\rho$ of
$\mathcal{T}_{k}$ with the lowest value function

6: Construct $\mathcal{T}_{k+1}$ under the constraint (32)

7: Compute the value function $V^{n}_{k+1}(x)$

8: $res=|V^{0}_{k}(x_{0})-V^{0}_{k+1}(x_{0})|$

9: $k=k+1$

10:end while

In Figure 1 we show an application of the statistical pruning under the Van
der Pol dynamics:

$\left\\{\begin{array}[]{l}\dot{y_{1}}(t)=y_{2}(t),\;\;t\in(0,T],\\\
\dot{y_{2}}(t)=\omega(1-y^{2}_{1}(t))y_{2}(t)-y_{1}(t)+u(t),\;\;t\in(0,T],\\\
(y_{1}(0),y_{2}(0))=(0.4,-0.3),\end{array}\right.$

where $\omega=0.15$, $T=1.4$ and $u:[0,T]\rightarrow[0,1].$ Fixing a time step
$\Delta t=0.2$, we display the initial full tree $\mathcal{T}_{1}$ with
discrete controls $\\{0,1\\}$, its refinement $\widetilde{\mathcal{T}}_{1}$
with $\rho=0.3$ and the new tree $\mathcal{T}_{2}$ with discrete controls
$\\{0,0.5,1\\}$.

Figure 1: Application of the statistical pruning to Van der Pol oscillator
with $\rho=0.3$.

###  5.3 Monotone control

In this section we restrict the admissible set of controls to monotone
controls, $e.g.$

$\overline{\mathcal{U}}=\\{u:[0,T]\rightarrow U,\,u(\cdot)\mbox{ monotone in
}[0,T]\\}.$

In economy different problems can be formulated as optimal control problems
with monotone controls ($e.g.$ adjustment theory of investment problems).
Under this constraint, in barron1985viscosity  Barron proved that the value
function is a generalized solution of the quasi-variational inequality and its
numerical treatment has been investigated in a series of papers
philipp2015discrete ; aragone2018fully . We consider the non decreasing case
without loss of generality. Let us introduce the notation which will be useful
in this section. We define as $\mathcal{T}_{M,N}$ the tree obtained using $M$
discrete controls and $N$ time steps, while we denote as
$\overline{\mathcal{T}}_{M,N}$ the tree constructed via monotone controls. In
Figure 2 we show the structure of the tree $\overline{\mathcal{T}}_{2,N}$. In
this case we are using $2$ discrete controls $u_{1}<u_{2}$. When we apply
$u_{2}$, the corresponding subtree will have just one node for each level,
since the control cannot decrease by hypothesis.

$x$$\zeta^{1}_{1}$$\zeta^{2}_{3}$$\zeta^{2}_{4}$$\ldots$$\zeta^{1}_{2}$$\ldots$
Figure 2: Example of the tree $\overline{\mathcal{T}}_{2,n}$

In this framework we have a great improvement in terms of the cardinality of
the tree, as stated in the following proposition.

###### Proposition 1

Given $M$ discrete controls and $N$ time steps, the cardinality of the tree
based on monotone controls is given by

$\left|\overline{\mathcal{T}}_{M,N}\right|=\frac{(M+N)!}{M!N!}$ (33)

Proof. We will proceed by induction on the pair $(M,N)$. It is easy to check
that $\left|\overline{\mathcal{T}}_{1,N}\right|=N+1$ and
$\left|\overline{\mathcal{T}}_{M,1}\right|=M+1$. Now let us suppose (33) holds
for a pair $(M,N)$. First, we are going to prove that the result holds for the
pair $(M,N+1)$. Given the particular structure of the tree, we can write

$\left|\overline{\mathcal{T}}_{M,N+1}\right|=\sum_{k=1}^{M}\left|\overline{\mathcal{T}}_{k,N}\right|+1=\frac{(M+N+1)!-M!(N+1)!}{M!(N+1)!}+1=\frac{(M+N+1)!}{M!(N+1)!},$

obtaining the result. Afterwards, let us demonstrate it for the pair
$(M+1,N)$. In this case we can split the tree in the following way

$\left|\overline{\mathcal{T}}_{M+1,N}\right|=\left|\overline{\mathcal{T}}_{M,N}\right|+\left|\overline{\mathcal{T}}_{M+1,N-1}\right|=\sum_{k=2}^{N}\left|\overline{\mathcal{T}}_{M,k}\right|+\left|\overline{\mathcal{T}}_{M+1,1}\right|$
$=\frac{(M+N+1)!-(M+2)!N!}{(M+1)!N!}+M+2=\frac{(M+N+1)!}{(M+1)!N!}$

and this completes the proof.

In general the cardinality of the tree grows as $O(M^{N})$ which is infeasible
due to the huge amount of memory allocations. Fixing the number of discrete
control $M$, the cardinality of the tree based on a monotone control grows as
$O(N^{M}/M!)$, yielding an affordable algorithm for the computation of the
optimal control with a high number of time steps.

###  5.4 Bilinear control

Let us consider the following bilinear dynamical system:

$\dot{y}=Ly+uy,\quad u\in U\subset\mathbb{R}.$ (34)

Discretizing (34) via a semi-implicit scheme, we obtain

$y^{n}=(I-\Delta tL)^{-1}y^{n-1}(1+\Delta tu^{n-1})=(I-\Delta
tL)^{-n}y_{0}\prod_{i=0}^{n-1}(1+\Delta tu^{i}).$ (35)

Let us consider now a new evolution $\tilde{y}^{n}$ of the discrete scheme at
time $t_{n}$ with controls $\\{\tilde{u}^{i}\\}_{i=0}^{n-1}$. Then the
distance between the two dynamics is given by

$\|y^{n}-\tilde{y}^{n}\|\leq\|(I-\Delta
tL)^{-n}y_{0}\|\left|\prod_{i=0}^{n-1}(1+\Delta
tu^{i})-\prod_{i=0}^{n-1}(1+\Delta t\tilde{u}^{i})\right|.$ (36)

Let us introduce a definition which will be useful in this section.

###### Definition 1

A discrete system with discrete controls satisfies the sum-based pruning
property if for every pair of vectors $(u^{0},\ldots,u^{n-1})$ and
$(\tilde{u}^{0},\ldots,\tilde{u}^{n-1})$ such that

$\sum_{i=0}^{n-1}u^{i}=\sum_{i=0}^{n-1}\tilde{u}^{i}$ (37)

the corresponding discrete solution $y^{n}$ and $\tilde{y}^{n}$ satisfy the
geometrical pruning rule, $e.g.$ $\|y^{n}-\tilde{y}^{n}\|\leq C\Delta t^{2}$.

This class of discrete systems benefits from an important improvement in terms
of the cardinality of the corresponding tree, as stated in the following
proposition. For the proof we refer to Proposition 3.12 in saluzzi2022error .

###### Proposition 2

The cardinality of the tree based on a system with $M$ discrete controls and
$N$ time steps satisfying the sum-based pruning property is at most
$\frac{N(N-1)}{2}(M-1)+N+1$.

The discrete dynamics (35) with two discrete controls belongs to the class of
the system verifying the sum-based pruning property as stated in the next
proposition.

###### Proposition 3

The discrete system (35) satisfies the sum-based pruning with $2$ discrete
controls. Hence, the cardinality of the corresponding tree is at most
$\frac{N(N+1)}{2}+1$.

###### Proof

Let us consider two pair of vectors $(u^{0},\ldots,u^{n-1})$ and
$(\tilde{u}^{0},\ldots,\tilde{u}^{n-1})$ verifying the sum-based pruning
property (37). Since we are considering two discrete control $(u_{1},u_{2})\in
U\times U$, the upper bound for distance (36) between the two corresponding
dynamics becomes

$\|(I-\Delta tL)^{-n}y_{0}\|\left|(1+\Delta tu_{1})^{k_{1}}(1+\Delta
tu_{2})^{n-k_{1}}-(1+\Delta tu_{1})^{k_{2}}(1+\Delta tu_{2})^{n-k_{2}}\right|$

with $k_{1},k_{2}\in\\{0,\ldots,n\\}$. By property (37) we immediately see
that either $k_{1}=k_{2}$ or $u_{1}=u_{2}$, which implies that the two
corresponding solutions coincide.

In this case we can directly construct the tree based on this structure,
without implementing any pruning criterion. The construction of the tree based
on two discrete controls may be used as a fast and cheap procedure to get
information about the full dimensional system. Once constructed the basis and
projected the system onto the lower dimensional space, it is possible to
consider an higher number of discrete controls.

##  6 An error bound for the multilinear HJB-POD-DEIM algorithm

The aim of this section is to derive an error estimate for the approximation
of value function with the HO-POD-DEIM algorithm applied to the tree
structure. The main reference of this section is chat2012 , where the authors
obtain a state space error bounds for the solutions of the reduced systems via
a POD-DEIM approach and the application of an implicit scheme for the time
integration. Following their proof, we are going to extend the result to semi-
implicit schemes in our multilinear setting.

First of all, we consider the vectorized form of dynamical system (16), whose
semi-implicit discretization with stepsize $\Delta t$ and discrete controls
$\\{u^{j}\\}_{j=0}^{N_{t}-1}$ reads

$\frac{y^{j}-y^{j-1}}{\Delta
t}=Ly^{j}+{\bm{f}}(y^{j-1},u^{j-1},t^{j-1}),\qquad j=1,\ldots,N_{t}.$ (38)

Taking into account the basis in vector form

$V_{Y}=V_{d}\otimes\cdots\otimes V_{1},\quad
V_{F}={\bf\Phi_{d}}\otimes\cdots\otimes{\bf\Phi_{1}},$
$\mathbb{P}=V_{F}\left(({\bf P_{d}}\otimes\cdots\otimes{\bf
P_{1}})^{\top}V_{F}\right)^{-1}\left({\bf P_{d}}\otimes\cdots\otimes{\bf
P_{1}}\right)^{\top},$

the vectorized form of the semi-implicit scheme for the reduced dynamics (27)
reads

$\frac{\hat{y}^{j}-\hat{y}^{j-1}}{\Delta
t}=V_{Y}^{\top}L\,V_{Y}\hat{y}^{j}+V_{Y}^{\top}\mathbb{P}{\bm{f}}(V_{Y}\hat{y}^{j-1},u^{j-1},t^{j-1}),\qquad
j=1,\ldots,N_{t}.$ (39)

Our aim is to prove an error estimate between full order scheme (38) and
reduced one (39).

For this purpose we introduce the logarithmic norm of matrix
$A\in\mathbb{C}^{n\times n}$ defined as

$\mu(A)=\sup_{x\in\mathbb{C}^{n}\setminus\\{0\\}}\frac{Re<Ax,x>}{\|x\|^{2}}.$
(40)

The logarithmic norm plays an important role for the stability analysis for
continuous and discrete linear dynamical systems. Indeed, it is possible to
prove that $\|e^{tA}\|\leq e^{t\mu(A)}$ $\forall t\geq 0$ (see $e.g.$
soderlind2006logarithmic ) and by this inequality we can state that the
dynamical system is stable if $\mu(A)\leq 0$. The definition of this norm will
be fundamental in the treatment of the implicit part of the scheme, while the
Lipschitz-continuity of f will be employed for the estimation of the explicit
part. In the following proposition we prove that the error between the full
order model (38) and the lifted reduced order model (39) depends on the
accuracy of the HO-POD and HO-DEIM basis. The proof can be found in Appendix
A.

###### Proposition 1

Given $\\{y^{k}\\}_{k=0}^{N_{t}}$ the solution of the (38) and
$\\{\hat{y}^{k}\\}_{k=0}^{N_{t}}$ solution of (39) with controls
$\\{u^{j}\\}_{j=0}^{N_{t}-1}$ and time step $\Delta t$ satisfying $\Delta
t\,\mu(V_{Y}^{\top}LV_{Y})<1$, then

$\sum_{k=0}^{N_{t}}|y^{k}-V_{Y}\hat{y}^{k}|^{2}\leq
C(T)\left(\mathcal{E}_{y}+\mathcal{E}_{f}\right)$ (41)

with

${\cal E}_{y}=\sum_{j=0}^{N_{t}}|y^{j}-V_{Y}V_{Y}^{\top}y^{j}|^{2},\quad{\cal
E}_{f}=\sum_{j=0}^{N_{t}-1}|{\bm{f}}(y^{j},t^{j},u^{j})-V_{F}V_{F}^{\top}{\bm{f}}(y^{j},t^{j},u^{j})|^{2}.$

Finally, we are ready to prove a convergence result for the continuous value
function $v(x,t)$, solution of the HJB equation (8), and the discrete value
function solution $\\{\widehat{V}^{n}(\widehat{x})\\}_{n}$, solution of the
scheme (28). For this purpose, let us define the continuous version of the DDP
for the full model

$\begin{cases}V(x,s)=\min\limits_{u\in
U}\\{V(x+(t_{n+1}-s)f(x,u,s),t_{n+1})+(t_{n+1}-s)\,L(x,u,s)\\},\\\
V(x,T)=g(x),\hskip
170.71652ptx\in\mathbb{R}^{d},s\in[t_{n},t_{n+1}),\end{cases}$ (42)

and its reduced version which reads:

$\begin{cases}\widehat{V}(\widehat{x},s)=\min\limits_{u\in
U}\\{\widehat{V}(\widehat{x}+(t_{n+1}-s)f^{\ell}(x^{\ell},u,s),t_{n+1})++(t_{n+1}-s)\,\widehat{L}(\widehat{x},u,s)\\},\\\
\widehat{V}(\widehat{x},T)=\widehat{g}(\widehat{x}),\hskip
170.71652pt\widehat{x}\in\mathbb{R}^{\ell},s\in[t_{n},t_{n+1}).\end{cases}$
(43)

Given the exact value function $v(x,s)$ and its continuous reduced
approximation $\widehat{V}(V_{Y}x,s)$, the following theorem provides an error
estimate for the approximation of the HJB equation by the HO-POD-DEIM
approach. The assumptions and the main procedure for the following result can
be found in Theorem 5.1 in Alla_Saluzzi_2020 .

###### Theorem 6.1

Given $v(x,s)$ the solution of the HJB equation (8) and its reduced
approximation $\widehat{V}(V_{Y}x,s)$ solution of the scheme (43), the
following estimate holds

$|v(x,s)-\widehat{V}(V_{Y}x,s)|\leq C(T)\left(\Delta
t+\mathcal{E}_{y}+\mathcal{E}_{f}\right).$ (44)

###### Proof

The proof follows closely the procedure adopted for Theorem 5.1 in
Alla_Saluzzi_2020 . The only difference arises in the estimation of the
projection error between the FOM and the lifted ROM solutions and in this case
we apply Proposition 1 to obtain the result.

##  7 Numerical tests

In this section we test the proposed technique in different frameworks. In the
first numerical test we consider a bilinear advection-diffusion equation,
comparing the vector and matrix cases for the construction of the reduced
basis. The second test is devoted to a nonlinear reaction-diffusion PDE where
we show the efficiency of the statistical pruning coupled with the MOR
technique. Finally, in the third test we consider a more challenging problem:
the control of the 3D viscous Burgers’ equations. We use this final example to
indicate the power of the proposed algorithm in terms of CPU time and memory
requirements with respect to the vector construction of the problem. The
numerical tests are performed on a Dell XPS 13 with Intel Core i7, 2.8GHz and
16GB RAM. The codes are written in Matlab R2022a.

###  7.1 Test 1: Advection-diffusion equation

In the first numerical test we conside the following bilinear advection-
diffusion equation:

$\begin{cases}\partial_{s}y+(c_{1},c_{2})\cdot\nabla y=\sigma\Delta
y+yu(s)&(x,s)\in\Omega\times[0,T],\\\
y(x,s)=0&(x,s)\in\partial\Omega\times[0,T],\\\
y(x,0)=y_{0}(x)&x\in\Omega,\end{cases}$ (45)

with

$u\in\mathcal{U}=\\{u:[0,T]\rightarrow U,u\,measurable\\}.$

The aim of the optimal control is to drive the solution to the equilibrium
$\overline{y}(x)\equiv 0$ and to this end we introduce the following cost
functional:

$J_{y_{0},t}(u)=\int_{t}^{T}\left(\int_{\Omega}|y(x,s)|^{2}dx+|u(s)|^{2}\right)ds+\int_{\Omega}|y(x,T)|^{2}dx.$

Since we are considering a bilinear optimal control problem, we can benefit of
the results presented in Section 5.4, $e.g.$ discretizing the control set with
two discrete controls, the total cardinality of the tree is order
$O(N_{t}^{2})$. We fix $T=1$, $U=[-3,-1]$, $y_{0}(x)=\max(2-\|x\|^{2},0)$,
$\Omega=[-5,5]^{2}$, $\widehat{\Delta t}=\Delta t=0.05$,
$\widehat{U}=\\{-3,-1\\}$, $c_{1}=0.5$, $c_{2}=0$ and $\sigma=0$. Later on we
will discuss the behaviour of the algorithm considering different choices for
the coefficients $c_{1},c_{2}$ and $\sigma$. Furthermore, we impose
$\tau=10^{-4}$ for both methods to obtain the same projection error. In this
setting the cardinality of the tree is $3321$. In Table 1 we show the
dimension of the basis varying the number of the grid points in each
direction. Since the system is driven along an axes, the HO-POD procedure
requires more basis in one direction with respect to the other. We note that
the maximum of the dimensions of the HO-POD basis is equal to the number of
POD basis for any choice of $n$.

$n$ | $101$ | $121$ | $141$ | $161$ | $181$ | $201$
---|---|---|---|---|---|---
POD | 7 | 7 | 7 | 7 | 8 | 8
HO-POD | $(3,7)$ | $(3,7)$ | $(3,7)$ | $(3,7)$ | $(3,8)$ | $(3,8)$
Table 1: Dim. of basis for POD and HO-POD varying the number of the grid
points per dimension with $\tau=10^{-4}$, $\Delta t=0.05$ and two discrete
controls.

In the top left panel of Figure 3 we compare the CPU time for the offline
phase for the POD and HO-POD algorithms. As stated previously, HO-POD requires
less storage and enables to treat with very high dimensional problems. In
particular, we note a difference of almost two orders of magnitude between the
POD and HO-POD offline stages for $n=201$. In the top right and bottom panels
of Figure 3 a comparison of the computational times and cardinality of the
pruned trees varying the number of discrete controls for the online phase and
fixing $n=161$ is presented. The HO-POD is again performing better than the
POD algorithm since the geometrical pruning turns to be more efficient in the
HO-POD setting. Indeed, fixing $M=7$ discrete controls, the cardinality of the
HO-POD tree reaches almost order $10^{5}$, against the order $\approx 10^{7}$
for the POD tree and $10^{16}$ for the unpruned tree.

Figure 3: Test 1: CPU time for the offline stage (top left), CPU time for the
online stage (top right) and cardinality of the pruned tree (bottom) for POD
and HO-POD techniques.

The optimal trajectory computed via HO-POD with $7$ discrete controls for
different time instances is displayed in Figure 4, noting that the solution is
getting closer to the stationary solution.

Figure 4: Test 1: Optimal trajectory for $t=0$ (left), $t=1$ (central) and
$t=2$ (right) for HO-POD and $n=161$.

Finally, Table 2 shows the different CPU times for the offline phase for the
two methods. First of all, we notice that HO-POD is faster in all cases, but
we obtain a particular speed-up in presence of a one-direction convection,
since the construction of the basis operates separately in each direction.

$(c_{1},c_{2},\sigma)$ | POD | HO-POD
---|---|---
$(1,0,0)$ | $17.78s$ | 0.67s
$(1,1,0)$ | 18s | 1.77s
$(0,0,1)$ | 16.9s | 0.84s
$(1,1,1)$ | 34s | 1.73s
Table 2: Offline CPU times for POD and HO-POD varying the coefficient
$(c_{1},c_{2},\sigma)$ with $n=601$ and $\widehat{\Delta t}=0.05$.

###  7.2 Test 2: Allen-Cahn equation

We consider the following nonlinear PDE with homogeneous Neumann boundary
conditions:

$\begin{cases}\partial_{s}y=\sigma\Delta
y+y\left(1-y^{2}\right)+y_{0}(x)u(s)&(x,s)\in\Omega\times[0,T],\\\
\partial_{n}y(x,s)=0&(x,s)\in\partial\Omega\times[0,T],\\\
y(x,0)=y_{0}(x)&x\in[-1,1]^{2}.\end{cases}$ (46)

Our aim is to steer the solution to the unstable equilibrium
$\overline{y}\equiv 0$ minimizing the following cost functional

$J_{y_{0},t}(u)=\int_{t}^{T}\left(\int_{\Omega}|y(x,s)|^{2}dx+\gamma|u(s)|^{2}\right)ds+\int_{\Omega}|y(x,T)|^{2}dx.$

We fix $T=1$, $\gamma=0.01$, $\sigma=0.1$, $\Omega=[-1,1]^{2}$ $U=[-2,0]$ and
$y_{0}(x)=2+\cos(2\pi x_{1})\cos(2\pi x_{2})$. Furthermore, we set $\Delta
t=0.1$, $\tau=10^{-3}$ and we discretize the domain $[-1,1]^{2}$ with $601$
equidistant points, obtaining a grid of $361201$ points. In the offline phase
we consider $2$ discrete controls and we construct a rough tree with $2047$
nodes. Since the problem is not linear, we apply the HO-POD-DEIM strategy and
the dimensions of the basis turns to be $k_{1}=k_{2}=p_{1}=p_{2}=5$, hence the
low dimension solution lives in $\mathbb{R}^{5\times 5}$. The offline stage
took $174$ seconds. In this example the dynamical system is nonlinear and we
do not have any a priori estimate for the introduced pruning criteria. Hence,
we are going to apply the statistical pruning discussed in Section 5.2. We fix
the ratio $\rho=0.2$ and we iterative the statistical pruning strategy
explained in Algorithm 1 with a stopping tolerance $tol=10^{-4}$. In the left
panel of Figure 5 we show the cardinality of the low dimensional tree in
logarithm scale. We recall that the cardinality of the full tree would be
$O(M^{11})$, where $M$ is the number of discrete controls, reaching an order
of $\approx$ $10^{17}$ in the case of $33$ controls. The application of the
statistical pruning achieves a great improvement in terms of memory storage,
gaining almost 12 orders of magnitudes with respect to the full tree with $33$
discrete controls. The total cost varying the number of controls is displayed
in the left panel of Figure 5. The cost functional shows a decreasing
behaviour as expected and the algorithms stops with $33$ controls since the
stopping rule has been satisfied. In Figures 6-7 the optimal trajectories at
time instances $t\in\\{0,0.5,1\\}$ and the control signal are displayed. We
note that the control signal is driving the solution to the unstable
equilibrium $\tilde{y}\equiv 0$.

Figure 5: Test 2: Cardinality of the tree in logarithmic scale (left) and
total cost (right) for $\Delta t=0.1$ and varying the number of controls.

Figure 6: Test 2: Optimal trajectory at $t=0$ (left) and $t=0.5$ (right) for
$\Delta t=0.1$ and $M=33$.

Figure 7: Test 2: Optimal trajectory at $t=1$ (left) and control signal
(right) for $\Delta t=0.1$ and $M=33$.

###  7.3 Test 3: 3D viscous Burgers’ equation

Here we consider the nonlinear 3D viscous Burgers’ equation (see, e.g.,
GAO2017 ) given by

$\begin{cases}\partial_{t}y_{1}&=\frac{1}{r}\Delta
y_{1}-\underline{y}\cdot\nabla y_{1}+y_{1}u\\\
\partial_{t}y_{2}&=\frac{1}{r}\Delta y_{2}-\underline{y}\cdot\nabla
y_{2}+y_{2}u,\\\ \partial_{t}y_{3}&=\frac{1}{r}\Delta
y_{3}-\underline{y}\cdot\nabla y_{3}+y_{3}u,\\\ \end{cases}$ (47)

where $y_{1}(x_{1},x_{2},x_{3},t)$, $y_{2}(x_{1},x_{2},x_{3},t)$ and
$y_{3}(x_{1},x_{2},x_{3},t)$ are the three velocities to be determined, with
$x=(x_{1},x_{2},x_{3})\in[0,1]^{3}$, $t\in[0,1]$ and the Reynold’s number
$r=100$. Furthermore, the system is subject to homogeneous Dirichlet boundary
conditions and initial states

$\begin{split}y_{1}(x,y,z,0)&=\frac{1}{10}\sin(2\pi x_{1})\sin(2\pi
x_{2})\cos(2\pi x_{3})\\\ y_{2}(x,y,z,0)&=\frac{1}{10}\sin(2\pi
x_{1})\cos(2\pi x_{2})\sin(2\pi x_{3})\\\
y_{3}(x,y,z,0)&=\frac{1}{10}\cos(2\pi x_{1})\sin(2\pi x_{2})\sin(2\pi
x_{3}).\\\ \end{split}$

A finite difference space discretization in the cube yields a system of ODEs
of the form (20), with nonlinear functions given by

$\begin{split}{\cal F}_{i}({\cal D}_{i}({\boldsymbol{\cal
Y}}_{i}),{\boldsymbol{\cal Y}}_{1},{\boldsymbol{\cal Y}}_{2},{\boldsymbol{\cal
Y}}_{3},t)&=\sum_{k=1}^{3}({\boldsymbol{\cal
Y}}_{i}\times_{k}{B}_{ki})\circ{\boldsymbol{\cal Y}}_{k},\end{split}$

for $i=1,2,3$, where ${B}_{1i}\in\mathbb{R}^{n\times n}$,
${B}_{2i}\in\mathbb{R}^{n\times n}$ and ${B}_{3i}\in\mathbb{R}^{n\times n}$
contain the coefficients for a first order centered difference space
discretization in the $x_{1}-$, $x_{2}-$ and $x_{3}-$ directions respectively,
and $n$ is the dimension of the discretized tensor in each spatial direction.
For a more detailed discussion on the space discretization and HO-POD-DEIM
model reduction of systems of ODEs in array form, we point the reader to
kirsten.22 , as well as the companion manuscript KSF2023 .

We consider the following cost functional

$J_{y_{0},t}(u)=\int_{t}^{T}\left(\int_{\Omega}\sum_{i=1}^{3}|y_{i}(x,s)|^{2}dx+\dfrac{1}{10}|u(s)|^{2}\right)ds+\int_{\Omega}\sum_{i=1}^{3}|y_{i}(x,T)|^{2}dx.$
(48)

The control $u(t)$ will be taken in the following admissible set of controls

$\mathcal{U}=\\{u:[0,T]\rightarrow[-2,0]\\}.$

We therefore construct one tree for the control $u$ containing the approximate
solution of each of the three equations at its nodes. Constructing and storing
this tree of course leads to extremely demanding memory requirements and
computational effort to construct the approximations spaces for the reduced
models.

We therefore use this experiment to illustrate the massive computational gain
of the HO-POD-DEIM method, in combination with the snapshot selection
algorithm and the low-rank storage algorithm. We first investigate the
computational load required in the offline phase by the HO-POD-DEIM method as
well as standard POD-DEIM applied to the system (47) discretized in vector
form. For the vectorized system we also apply a semi-implicit Euler time
discretization to each of the three equations, and each linear system is
solved using the Matlab function pcg preconditioned with an incomplete
Cholesky factorization with a drop tolerance of $10^{-4}$.

Below we illustrate the computational load both in terms of CPU time and
memory requirements. On the left of Figure 8 we plot the computational time
required by both methods to construct the full-dimensional tree, with
$N_{t}=10$ and two controls, whose nodal values are used to construct either
the HO-POD-DEIM basis or the standard POD-DEIM basis. To construct the tree,
all nodal values from the previous time level need to be stored. To this end,
one of the computational bottlenecks in the construction of the reduced model
is memory requirements. Consequently, we further illustrate the power of the
proposed algorithm on the right of Figure 8, where we plot the maximum memory
required (in mb) at any point in the offline phase of the respective
algorithms. The plot indicates a massive difference in memory requirements,
mainly related to the low-rank basis construction used in the HO-POD-DEIM
algorithm, as well as the low-rank memory allocation method discussed in
Section 4.2. Furthermore we notice, that no data points are plotted in the
vector case for $n>60$, as this is where the computer ran out of its available
computational memory. Both plots are with respect to increasing $n$ and we
select $\tau=10^{-2}$ and $\kappa=20$ a priori.

Figure 8: Test: CPU time (left) and memory requirements (right) for both
methods applied to (47).

In what follows we consider the reduced model constructed by the HO-POD-DEIM
method and investigate the efficiency of the reduction. In Table 3 we indicate
the dimensions of the reduced approximation spaces determined to comply with
$\tau=10^{-2}$. To construct a tree with ten time steps and two controls,
pruned by the standard geometric pruning technique, the reduced model with
dimensions as in Table 3 required merely 19 seconds in comparison to the 361
seconds required by the full order model.

$y$ | $k_{1}$ | $k_{2}$ | $k_{3}$ | $p_{1}$ | $p_{2}$ | $p_{3}$ | error
---|---|---|---|---|---|---|---
$y_{1}$ | 6 | 11 | 12 | 10 | 18 | 19 | $2\cdot 10^{-2}$
$y_{2}$ | 6 | 15 | 13 | 6 | 20 | 17 | $3\cdot 10^{-2}$
$y_{3}$ | 6 | 12 | 13 | 5 | 18 | 18 | $2\cdot 10^{-2}$
Table 3: Dim. of basis and the average relative error compared to the full
order model with dimension $n=60$ and $\tau=10^{-2}$.

Finally, we also plot the cost functional in Figure 9 for both the full and
reduced order models as well as the optimal trajectories (unfolded along the
first mode) for all three equations at $t=0$, $t=0.5$ and $t=1$ in Figure 10.
Both these plots, indicate the convergence to the equilibrium of the reduced
model. We note that there is a visual superposition of the two curves of the
cost functional, demonstrating the effectiveness of the proposed methodology
for determining the optimal trajectory.

Figure 9: Cost functional for the optimal control.

Figure 10: All three controlled trajectories from top to bottom respectively,
unfolded along the first mode. We plot $t=0$ (left), $t=0.5$ (middle) and
$t=1$ (right).

##  8 Conclusions

In this paper we have introduced a new algorithm for approximating optimal
feedback controls related to optimal control problems driven by evolutionary
partial differential equations. The new algorithm is based on a tree structure
to avoid the construction of a grid in the solution of the HJB equations, and
exploits the compact representation of the dynamical systems based on tensor
notations via a higher-order model reduction approach. We have shown how the
algorithm can be constructed for general nonlinear control problems and given
some crucial hints on its implementation. Furthermore, we have studied the
existing pruning techniques for reducing the cardinality of the constructed
tree, and introduced a new statistical pruning technique for a further
reduction in the cardinality of the tree. To guarantee the convergence of the
method, we derived an error estimate depending on the time step and on the
accuracy of the HO-POD-DEIM basis. Finally, numerical tests on a number of
challenging benchmark problem have been discussed, indicating the power of the
method, with large savings both in terms of computational time and,
especially, memory. We believe that these promising numerical results brings
us one step closer to the application of DP in challenging, industrial
settings.

To this end, we plan to, in the near future, explore more challenging
industrial problems where the combination of the compact tensor representation
of the problem and the tree structure algorithm can give a competitive
advantage to DP for feedback control problems over possible competitors.

## Acknowledgements

A large part of this project was started together with Prof. Maurizio Falcone,
shortly before he passed away. We greatly acknowledge his contributions,
inspirations and leadership in this project. May his memory live a long time
in the mathematical community.

Funding. The second author is a member of INDAM GNCS (Gruppo Nazionale di
Calcolo Scientifico). This research has been partially supported by the INdAM-
GNCS project ”Finanziamento Giovani Ricercatori 2020-2021”.

Data access. Matlab codes implementing the numerical examples are available at
https://github.com/saluzzi/Multilinear_HJB_POD.

## Declarations

Conflict of interest. The authors declare that they have no conflict of
interest.

##  Appendix A Proof of Proposition 1

###### Proof

For the sake of simplicity we are going to define
$\hat{{L}}=V_{Y}^{\top}LV_{Y}$, $f_{k}={\bm{f}}(y^{k},t^{k},u^{k})$,
$\hat{f}_{k}=V_{Y}^{\top}\mathbb{P}{\bm{f}}(V_{Y}\hat{y}^{k},u^{k},t^{k})$ and
$f_{k,V}=V_{Y}^{\top}\mathbb{P}{\bm{f}}(V_{Y}V^{\top}_{Y}y^{k},u^{k},t^{k})$,
where we are considering the same control sequence
$\\{u_{k}\\}_{k=0}^{N_{t}-1}$.

We consider the error at time $t_{j}$ between the full model and the lifted
reduced model as

$E_{j}=y^{j}-V_{Y}\hat{y}^{j}$

and we rewrite it as a sum of two quantities

$E_{j}=\rho_{j}+\theta_{j}$

where

$\rho_{j}=y^{j}-V_{Y}V_{Y}^{\top}y^{j},\quad\theta_{j}=V_{Y}V_{Y}^{\top}y^{j}-V_{Y}\hat{y}^{j}.$

Multiplying (38) by $V_{Y}^{\top}$ and adding and subtracting
$\hat{L}V_{Y}^{\top}y^{j}+f_{j-1,V}$ we get

$V_{Y}^{\top}\frac{y^{j}-y^{j-1}}{\Delta
t}=\hat{L}V_{Y}^{\top}y^{j}+f_{j-1,V}+\hat{R}_{j},$

with

$\hat{R}_{j}=V_{Y}^{\top}Ly^{j}+V_{Y}^{\top}f_{j-1}-\hat{L}V_{Y}^{\top}y^{j}-f_{j-1,V}.$

Defining $\hat{\theta}_{j}=V_{Y}^{\top}\theta_{j}$, we obtain

$\frac{\hat{\theta}_{j}-\hat{\theta}_{j-1}}{\Delta
t}=V_{Y}^{\top}\frac{y^{j}-y^{j-1}}{\Delta
t}-\frac{\hat{y}^{j}-\hat{y}^{j-1}}{\Delta
t}=\hat{L}V_{Y}^{\top}y^{j}+f_{j-1,V}+\hat{R}_{j}-\hat{L}\hat{y}^{j}-\hat{f}_{j-1}.$

Since
$\left<\hat{\theta}_{j},\hat{\theta}_{j-1}\right>\leq\|\hat{\theta}_{j}\|\|\hat{\theta}_{j-1}\|$,
we get

$\frac{\|\hat{\theta}_{j}\|-\|\hat{\theta}_{j-1}\|}{\Delta
t}\leq\frac{1}{\|\hat{\theta}_{j}\|}\left<\hat{\theta}_{j},\frac{\hat{\theta}_{j}-\hat{\theta}_{j-1}}{\Delta
t}\right>$
$=\frac{1}{\|\hat{\theta}_{j}\|}\left(\left<\hat{\theta}_{j},\hat{L}\left(V_{Y}^{\top}y^{j}-\hat{y}^{j}\right)\right>+\left<\hat{\theta}_{j},f_{j-1,V}-\hat{f}_{j-1}+\hat{R}_{j}\right>\right)$
$\leq\mu(\hat{L})\|\hat{\theta}_{j}\|+\gamma\|\hat{\theta}_{j-1}\|+\|\hat{R}_{j}\|,$

where $\gamma=L_{\bf f}\|V_{Y}^{\top}\mathbb{P}\|$ and we used the definition
of logarithmic norm (40) and the Lipschitz-continuity of the function ${\bf
f}$. Defining $\zeta=\frac{1}{1-\Delta t\mu(\hat{L})}$ and $\eta=1+\Delta
t\gamma$ and by the fact that $\|\theta_{j}\|=\|\hat{\theta}_{j}\|$, it
follows

$\|\theta_{j}\|\leq\zeta\eta\|\theta_{j-1}\|+\Delta
t\,\zeta\|\hat{R}_{j}\|\leq(\zeta\eta)^{j}\|\theta_{0}\|+\Delta
t\sum_{k=1}^{j}\zeta^{k}\eta^{k-1}\|\hat{R}_{j-k+1}\|$ $\leq\Delta
t\zeta\left(\sum_{k=0}^{j-1}(\zeta\eta)^{2k}\sum_{k=1}^{j}\|\hat{R}_{k}\|^{2}\right)^{1/2},$

where we note that $\theta_{0}=0$ and $\zeta$ is positive due to the
assumption on the time step $\Delta t$. Let us define
$q=\sum_{k=0}^{N_{t}-1}(\zeta\eta)^{2k}$. Recalling the definition of
$\hat{R}_{k}$

$\hat{R}_{k}=V_{Y}^{\top}L\left(y^{k}-V_{Y}V_{Y}^{\top}y^{k}\right)+V_{Y}^{\top}\left(f_{k-1}-\mathbb{P}{\bf
f}(V_{Y}V^{\top}_{Y}y^{k-1},u^{k-1},t^{k-1})\right),$

we note that

$\|V_{Y}^{\top}\left(f_{k-1}-\mathbb{P}{\bf
f}(V_{Y}V^{\top}_{Y}y^{k-1},u^{k-1},t^{k-1})\right)\|$
$=\|V_{Y}^{\top}\left(f_{k-1}-\mathbb{P}f_{k-1}+\mathbb{P}f_{k-1}-\mathbb{P}{\bf
f}(V_{Y}V^{\top}_{Y}y^{k-1},u^{k-1},t^{k-1})\right)\|$ $\leq{\bf
c}\|V_{Y}^{\top}\|\|f_{k-1}-V_{F}V_{F}^{\top}f_{k-1}\|+\gamma\|\rho_{k-1}\|$

where we applied Proposition 1 from kirsten.22  with ${\bf
c}=\prod_{m=1}^{d}\|({\bf P_{m}}^{\top}{\bf\Phi_{m}})^{-1}\|$ and the
Lipschitz-continuity of the function ${\bf f}$ . Therefore, we obtain the
following upper bound for the term $\hat{R}_{k}$

$\|\hat{R}_{k}\|\leq\alpha\|\rho_{k}\|+\beta\|w_{k-1}\|$

where $\alpha=\|V_{Y}^{\top}L\|+\gamma$, $\beta={\bf c}\|V_{Y}^{\top}\|$ and
$w_{k-1}=f_{k-1}-V_{F}V_{F}^{\top}f_{k-1}$.

From these results we can get the following estimate for the generic term
$\theta_{j}$

$\|\theta_{j}\|^{2}\leq(\Delta
t\zeta)^{2}q\sum_{k=1}^{j}\|\hat{R}_{k}\|^{2}\leq 2(\Delta
t\zeta)^{2}q\sum_{k=1}^{j}(\alpha^{2}\|\rho_{k}\|^{2}+\beta^{2}\|w_{k-1}\|^{2})$

and finally

$\sum_{j=0}^{N_{t}}\|E_{j}\|^{2}=\sum_{j=0}^{N_{t}}\|\rho_{j}\|^{2}+\sum_{j=1}^{N_{t}}\|\theta_{j}\|^{2}\leq
C(T)\left(\mathcal{E}_{y}+\mathcal{E}_{f}\right)$

where

$C(T)=\max\\{1+2q\zeta^{2}T\Delta t\,\alpha^{2},2q\zeta^{2}T\Delta
t\beta^{2}\\}.$

###### Remark A.1

Supposing that $\gamma\leq-\mu(\hat{L})$, then $\zeta\eta<1$ and we obtain the
following upper bound for the quantity $q$

$q=\sum_{k=0}^{N_{t}-1}(\zeta\eta)^{2k}\leq\frac{1}{1-(\zeta\eta)^{2N_{t}}}.$

###### Remark A.2

The constant $C(T)$ depends on the coefficient ${\bf c}=\prod_{m=1}^{d}\|({\bf
P_{m}}^{\top}{\bf\Phi_{m}})^{-1}\|,$ which is minimized applying the q-deim
procedure, we refer to gugercin2018  for more details.

## References

  * (1) Akian, M., Gaubert, S., Lakhoua, A.: 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), 817–848 (2008).  DOI 10.1137/060655286.  URL http://dx.doi.org/10.1137/060655286
  * (2) Akian, M., Gaubert, S., Liu, S.: An adaptive multi-level max-plus method for deterministic optimal control problems.  arXiv preprint arXiv:2304.10342 (2023)
  * (3) Alla, A., Falcone, M., Saluzzi, L.: An Efficient DP Algorithm on a Tree-Structure for Finite Horizon Optimal Control Problems.  SIAM Journal on Scientific Computing 41(4), A2384–A2406 (2019)
  * (4) Alla, A., Falcone, M., Saluzzi, L.: A tree structure algorithm for optimal control problems with state constraints.  Rendiconti di Matematica e delle Sue Applicazioni 41, 193–221 (2020)
  * (5) Alla, A., Oliveira, H., Santin, G.: HJB-RBF based approach for the control of PDEs.  arXiv preprint arXiv:2108.02987 (2021)
  * (6) Alla, A., Saluzzi, L.: A HJB-POD approach for the control of nonlinear PDEs on a tree structure.  Applied Numerical Mathematics 155, 192–207 (2020).  DOI 10.1016/j.apnum.2019.11.023.  URL http://dx.doi.org/10.1016/j.apnum.2019.11.023
  * (7) Aragone, L.S., Parente, L.A., Philipp, E.A.: Fully discrete schemes for monotone optimal control problems.  Computational and Applied Mathematics 37, 1047–1065 (2018)
  * (8) Bardi, M., Capuzzo-Dolcetta, I.: Optimal Control and Viscosity Solutions of Hamilton-Jacobi-Bellman Equations.  Modern Birkhäuser Classics. Birkhäuser Boston (2008)
  * (9) Barron, E.N.: Viscosity solutions for the monotone control problem.  SIAM journal on control and optimization 23(2), 161–171 (1985)
  * (10) Benner, P., Bujanović, Z., Kürschner, P., Saak, J.: A numerical comparison of different solvers for Large-Scale, Continuous-Time Algebraic Riccati Equations and LQR problems.  SIAM Journal on Scientific Computing 42(2), A957–A996 (2020).  DOI 10.1137/18m1220960.  URL http://dx.doi.org/10.1137/18M1220960
  * (11) Chaturantabut, S., Sorensen, D.C.: Nonlinear model reduction via discrete empirical interpolation.  SIAM J. Sci. Comput. 32(5), 2737–2764 (2010)
  * (12) Chaturantabut, S., Sorensen, D.C.: A state space error estimate for POD-DEIM nonlinear model reduction.  SIAM J Numer Anal 50(1), 46–63 (2012)
  * (13) Darbon, J., Dower, P.M., Meng, T.: Neural network architectures using min-plus algebra for solving certain high-dimensional optimal control problems and Hamilton-Jacobi PDEs.  Mathematics of Control, Signals, and Systems 35(1), 1–44 (2023)
  * (14) Darbon, J., Langlois, G.P., Meng, T.: 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
  * (15) D’Autilia, M.C., Sgura, I., Simoncini, V.: Matrix-oriented discretization methods for reaction–diffusion PDEs: Comparisons and applications.  Computers & Mathematics with Applications pp. 2067–2085 (2020)
  * (16) Dolgov, S., Kalise, D., Saluzzi, L.: Data-driven tensor train gradient cross approximation for Hamilton-Jacobi-Bellman equations.  arXiv preprint arXiv:2205.05109 (2022)
  * (17) Drmač, Z., Gugercin, S.: A new selection operator for the discrete empirical interpolation method—improved a priori error bound and extensions.  SIAM J. Sci. Comput. 38(2), A631–A648 (2016)
  * (18) Falcone, M., Ferretti, R.: Semi-Lagrangian approximation schemes for linear and Hamilton-Jacobi equations.  SIAM (2013)
  * (19) Falcone, M., Kirsten, G., Saluzzi, L.: Approximation of optimal control problems for the Navier-Stokes equation via multilinear HJB-POD.  Applied Mathematics and Computation 442, 127722 (2023).  DOI https://doi.org/10.1016/j.amc.2022.127722.  URL https://www.sciencedirect.com/science/article/pii/S0096300322007901
  * (20) Gao, Q., Zou, M.: An analytical solution for two and three dimensional nonlinear Burgers’ equation.  Appl. Math. Modell. 45, 255 – 270 (2017).  DOI https://doi.org/10.1016/j.apm.2016.12.018.  URL http://www.sciencedirect.com/science/article/pii/S0307904X16306710
  * (21) Garcke, J., Kröner, A.: Suboptimal feedback control of PDEs by solving HJB equations on adaptive sparse grids.  Journal of Scientific Computing 70(1), 1–28 (2016).  DOI 10.1007/s10915-016-0240-7.  URL http://dx.doi.org/10.1007/s10915-016-0240-7
  * (22) Golub, G.H., van Loan, C.F.: Matrix Computations, fourth edn.  Johns Hopkins University Press, Baltimore (2013).  URL http://www.cs.cornell.edu/cv/GVL4/golubandvanloan.htm
  * (23) Han, J., Jentzen, A., E, W.: Solving high-dimensional partial differential equations using deep learning.  Proceedings of the National Academy of Sciences 115(34), 8505–8510 (2018).  DOI 10.1073/pnas.1718942115.  URL http://dx.doi.org/10.1073/pnas.1718942115
  * (24) Hinze, M., Volkwein, S.: Proper orthogonal decomposition surrogate models for nonlinear dynamical systems: Error estimates and suboptimal control.  In: Dimension reduction of large-scale systems, pp. 261–306. Springer (2005)
  * (25) Kirsten, G.: Multilinear POD-DEIM model reduction for 2D and 3D nonlinear systems of differential equations.  Journal of Computational Dynamics 9(2), 159–183 (2022)
  * (26) Kirsten, G., Simoncini, V.: A matrix-oriented POD-DEIM algorithm applied to nonlinear differential matrix equations (2020).  ArXiv 2006.13289
  * (27) Kirsten, G., Simoncini, V.: Order reduction methods for solving large-scale differential matrix Riccati equations.  SIAM J. Sci. Comput. 42(4), A2182–A2205 (2020)
  * (28) Kolda, T.G., Bader, B.W.: Tensor decompositions and applications.  SIAM Rev 51(3), 455–500 (2009)
  * (29) Kunisch, K., Volkwein, S.: Optimal snapshot location for computing POD basis functions.  ESAIM: Mathematical Modelling and Numerical Analysis 44(3), 509–529 (2010)
  * (30) Kunisch, K., Volkwein, S., Xie, L.: HJB-POD based feedback design for the optimal control of evolution problems.  SIAM J. on Applied Dynamical Systems 4, 701–722 (2004)
  * (31) Kunisch, K., Walter, D.: Semiglobal optimal feedback stabilization of autonomous systems via deep neural network approximation.  ESAIM: Control, Optimisation and Calculus of Variations 27, 16 (2021).  DOI 10.1051/cocv/2021009.  URL http://dx.doi.org/10.1051/cocv/2021009
  * (32) Kunisch, K., Xie, L.: POD-based feedback control of burgers equation by solving the evolutionary HJB equation.  Computers and Mathematics with Applications 49, 1113–1126 (2005)
  * (33) McEneaney, W.M.: A Curse-of-Dimensionality-Free Numerical Method for Solution of Certain HJB PDEs.  SIAM Journal on Control and Optimization 46(4), 1239–1276 (2007).  DOI 10.1137/040610830.  URL http://dx.doi.org/10.1137/040610830
  * (34) Meng, T., Zhang, Z., Darbon, J., Karniadakis, G.E.: Sympocnet: Solving optimal control problems with applications to high-dimensional multi-agent path planning problems (2022).  URL https://arxiv.org/abs/2201.05475.  Doi: 10.48550/ARXIV.2201.05475
  * (35) Onken, D., Nurbekyan, L., Li, X., Fung, S.W., Osher, S., Ruthotto, L.: A neural network approach applied to multi-agent optimal control.  In: 2021 European Control Conference (ECC). IEEE (2021).  DOI 10.23919/ecc54610.2021.9655103
  * (36) Oster, M., Sallandt, L., Schneider, R.: Approximating optimal feedback controllers of finite horizon control problems using hierarchical tensor formats.  SIAM Journal on Scientific Computing 44(3), B746–B770 (2022)
  * (37) Palitta, D., Simoncini, V.: Matrix-equation-based strategies for convection–diffusion equations.  BIT Numerical Mathematics 56(2), 751–776 (2016)
  * (38) Philipp, E.A., Aragone, L.S., Parente, L.A.: Discrete time schemes for optimal control problems with monotone controls.  Computational and Applied Mathematics 34(3), 847–863 (2015)
  * (39) Pichi, F., Strazzullo, M., Ballarin, F., Rozza, G.: Driving bifurcating parametrized nonlinear PDEs by optimal control strategies: application to Navier–Stokes equations with model order reduction.  ESAIM: Mathematical Modelling and Numerical Analysis 56(4), 1361–1400 (2022)
  * (40) Quarteroni, A., Rozza, G.: Numerical solution of parametrized Navier-Stokes equations by reduced basis methods.  Numerical Methods for Partial Differential Equations 23(4), 923–948 (2007)
  * (41) Richter, L., Sallandt, L., Nüsken, N.: Solving high-dimensional parabolic PDEs using the tensor train format.  In: International Conference on Machine Learning, pp. 8998–9009 (2021)
  * (42) Ruthotto, L., Osher, S.J., Li, W., Nurbekyan, L., Fung, S.W.: 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), 9183–9193 (2020)
  * (43) Saluzzi, L., Alla, A., Falcone, M.: Error estimates for a tree structure algorithm solving finite horizon control problems.  ESAIM: Control, Optimisation & Calculus of Variations 28 (2022)
  * (44) Simoncini, V.: Computational methods for linear matrix equations.  SIAM Rev 58(3), 377–441 (2016)
  * (45) Simoncini, V.: Numerical solution of a class of third order tensor linear equations.  Bollettino dell’Unione Matematica Italiana 13(3), 429–439 (2020)
  * (46) Söderlind, G.: The logarithmic norm. history and modern theory.  BIT Numerical Mathematics 46, 631–652 (2006)
  * (47) Sorensen, D.C., Embree, M.: A DEIM induced CUR factorization.  SIAM J. Sci. Comput. 38(3), A1454–A1482 (2016)
  * (48) Stabile, G., Rozza, G.: Finite volume POD-Galerkin stabilized reduced order methods for the parametrized incompressible Navier-Stokes equations.  Computers & Fluids 173, 923–948 (2018)
  * (49) Vannieuwenhoven, N., Vandebril, R., Meerbergen, K.: A new truncation strategy for the higher-order singular value decomposition.  SIAM J. Sci. Comput. 34(2), A1027–A1052 (2012)
  * (50) Volkwein, S.: Model reduction using proper orthogonal decomposition.  Lecture Notes, Institute of Mathematics and Scientific Computing, University of Graz 1025 (2011)
  * (51) Zhou, M., Han, J., Lu, J.: Actor-Critic Method for High Dimensional Static Hamilton–Jacobi–Bellman Partial Differential Equations based on Neural Networks.  SIAM Journal on Scientific Computing 43(6), A4043–A4066 (2021).  DOI 10.1137/21m1402303.  URL https://doi.org/10.1137\%2F21m1402303