# Optimal bounds for POD approximations of infinite horizon control problems
based on time derivatives

Javier de Frutos Instituto de Investigación en Matemáticas (IMUVA),
Universidad de Valladolid, Spain. Research supported by Spanish MCIN/AEI under
grants PID2019-104141GB-I00 and PID2022-136550NB-I00 co-finanzed by FEDER (EU)
funds and by Junta de Castilla y León under grant VA169P20 co-finanzed by
FEDER (EU) funds<EMAIL_ADDRESS>Bosco García-Archilla Departamento de
Matemática Aplicada II, Universidad de Sevilla, Sevilla, Spain. Research is
supported by Spanish MCIN/AEI under grants PID2019-104141GB-I00 and
PID2022-136550NB-I00 co-finanzed by FEDER (EU) funds<EMAIL_ADDRESS>Julia
Novo Departamento de Matemáticas, Universidad Autónoma de Madrid, Spain.
Research supported Spanish MCIN/AEI under grants PID2019-104141GB-I00 and
PID2022-136550NB-I00 co-finanzed by FEDER (EU) funds and by Junta de Castilla
y León under grant VA169P20 co-finanzed by FEDER (EU) funds
<EMAIL_ADDRESS>
###### Abstract

In this paper we consider the numerical approximation of infinite horizon
problems via the dynamic programming approach. The value function of the
problem solves a Hamilton-Jacobi-Bellman (HJB) equation that is approximated
by a fully discrete method. It is known that the numerical problem is
difficult to handle by the so called curse of dimensionality. To mitigate this
issue we apply a reduction of the order by means of a new proper orthogonal
decomposition (POD) method based on time derivatives. We carry out the error
analysis of the method using recently proved optimal bounds for the fully
discrete approximations. Moreover, the use of snapshots based on time
derivatives allow us to bound some terms of the error that could not be
bounded in a standard POD approach. Some numerical experiments show the good
performance of the method in practice.

Key words. Dynamic programming, Hamilton-Jacobi-Bellman equation, optimal
control, proper orthogonal decomposition, snapshots based on time derivatives,
error analysis.

##  1 Introduction

In this paper we consider the numerical approximation of optimal control
problems. The subject is of importance for many applications such as aerospace
engineering, chemical processing and resource economics, among others.

The value function of an optimal control problem is obtained in terms of a
first-order nonlinear Hamilton-Jacobi-Bellman (HJB) partial differential
equation. A bottleneck in the computation of the value function comes from the
need to approach a nonlinear partial differential equation in dimension $n$,
which is a challenging problem in high dimensions.

Several methods have been studied in the literature trying to mitigate the so
called curse of dimensionality although it is still a difficult task. As
stated in [5], the relevance of efficient numerical methods can be seen by the
fact that methods solving the HJB equation are rarely used in practice due to
the necessary computational effort. We mention some related references that
are not intended to be a complete list. In [8], a domain decomposition
technique is considered. In [13] semi-Lagrangian methods are studied. The
authors in [12] apply data-based approximate policy iteration methods. A
procedure for the numerical approximation of high-dimensional HJB equations
associated to optimal feedback control problems for semilinear parabolic
equations is proposed in [3]. In [4] a tensor decomposition approach is
presented. In [5] an approach based on low-rank tensor train decompositions is
applied, see also the references therein.

In the present paper we concentrate on reduced order models based on proper
orthogonal decomposition (POD) methods. Our work is related to [1]. In this
reference the authors propose two different ways to apply POD methods in the
numerical approximation of the fully-discrete value function. In the first
approach, the authors choose a set of nodes in the original domain
$\Omega\subset{\mathbb{R}}^{n}$ and project then onto a reduced space
$\Omega^{r}\subset{\mathbb{R}}^{r}$ with $r<n$ to get a new set of nodes. The
problem in this procedure is that it produces a nonuniform grid in which the
mesh diameter cannot be predicted a priori. Consequently, the method is not
suitable to implement in practice. Furthermore, although this is not reflected
in the error bounds in [1], the error also depends on the interpolation
properties of the a priori unknown reduced mesh in
$\Omega^{r}\subset{\mathbb{R}}^{r}$. In the second approach, the authors use a
uniform mesh over the reduced space $\Omega^{r}$. This second method can be
implemented in practice (the numerical experiments in [1] are carried out with
this method). However, as the authors state, the computation of an upper bound
for the error in this case is much more involved and the error bound proved in
[1] has some drawbacks see [1, Remark 4.7, Remark 4.9].

Recently in [9], a new error analysis is introduced in which a bound of size
$O(h+k)$ is obtained for the fully discrete approximations to infinite horizon
problems via the dynamic programming approach. In this error bound, $h$ is the
time step while $k$ is the spatial mesh diameter. This error bound improves
existing results in the literature, where only $O(k/h)$ error bounds are
proved.

To bound the error in the first method in [1] the authors follow the technique
in [6, Corollary 2.4], [7, Theorem 1.3] obtaining a bound for the error of
size $O(k/h)$. For the second method in [1], the factor $1/h$ also multiplies
all the terms on the right-hand side of the a priori error bound.

In this paper we present a new approach, similar to the second method in [1],
but with snapshots based on the value at different times of the time
derivative of the state of the controlled nonlinear dynamical system, instead
of values of the state at different times. This new approach is inspired in
the recent results in [10] where the authors prove that the use of snapshots
based on time derivatives has the advantage of providing pointwise estimates
for the error between a function and its projection onto the POD space.
Moreover, we carry out a different error analysis based on the recent results
obtained in [9] that allow us to get sharper error bounds free of $1/h$
factors. This is in agreement with the numerical investigations in the
literature where the $1/h$ behaviour in the error bounds of fully discrete
methods has never been observed. Also, the use of snapshots based on time
derivatives allow us to give a bound for some terms that could not be bounded
with the standard approach. Both facts, the new technique used to bound the
error that follows ideas in [9], together with the use of snapshots based on
time derivatives, allow us to get error bounds for the new method that are
optimal in terms of the time step $h$ and the mesh diameter of the reduced
space $k_{r}$. As usual, our error bounds for the POD method depend also on
the size of the tail of eigenvalues in the singular value decomposition.

The outline of the paper is as follows. In Section 2 we state the model
problem and some preliminary results. In Section 3 we introduce the POD
approximation and carry out the error analysis of the method. Finally, in
Section 4 we show some numerical experiments in which we implement the method
we propose in the paper. In the experiments of Section 4 we choose the same
numerical tests as in [1] to compare our results with the those in this
related reference. The method introduced in the present paper seems to produce
better results than those shown in [1]. We finish the paper with some
conclusions.

##  2 Model problem and standard numerical approximation

In the sequel, $\|\cdot\|$ denotes any norm associated to an inner product and
$\|\cdot\|_{\infty}$ denotes the maximum norm for vectors in
${\mathbb{R}}^{n}$, $n\geq 1$. We will also denote by $\|\cdot\|_{2}$ the
standard euclidean norm.

For a nonlinear mapping

$f:{\mathbb{R}}^{n}\times{\mathbb{R}}^{m}\rightarrow{\mathbb{R}}^{n},$

and a given initial condition $y_{0}\in{\mathbb{R}}^{n}$ let us consider the
controlled nonlinear dynamical system

$\dot{y}(t)=f(y(t),u(t))\in{\mathbb{R}}^{n},\quad t>0,\quad
y(0)=y_{0}\in{\mathbb{R}}^{n},$ (1)

together with the infinite horizon cost functional

$J(y,u)=\int_{0}^{\infty}g(y(t),u(t))e^{-\lambda t}~{}dt.$ (2)

In (2) $\lambda>0$ is a given weighting parameter and

$g:{\mathbb{R}}^{n}\times{\mathbb{R}}^{m}\rightarrow{\mathbb{R}}.$

The set of admissible controls is

${\mathbb{U}}_{\rm ad}=\left\\{u\in{\mathbb{U}}\mid u(t)\in U_{\rm ad}\ {\rm
for}\ {\rm almost}\ {\rm all}\ t\geq 0\right\\},$

where ${\mathbb{U}}=L^{2}(0,\infty;{\mathbb{R}}^{m})$ and $U_{\rm
ad}\subset{\mathbb{R}}^{m}$ is a compact convex subset.

As in [1, Assumption 2.1] we assume the following hypotheses:

  * •

The right-hand side $f$ in (1) is continuous and globally Lipschitz-continuous
in both the first and second arguments; i.e., there exists a constant
$L_{f}>0$ satisfying

$\displaystyle\|f(y,u)-f(\tilde{y},u)\|$ $\displaystyle\leq$ $\displaystyle
L_{f}\|y-\tilde{y}\|,\quad\forall y,\tilde{y}\in{\mathbb{R}}^{n},u\in U_{\rm
ad},$ (3) $\displaystyle\|f(y,u)-f(y,\tilde{u})\|$ $\displaystyle\leq$
$\displaystyle L_{f}\|u-\tilde{u}\|,\quad\forall u,\tilde{u}\in U_{\rm
ad},y\in{\mathbb{R}}^{n}.$ (4)

  * •

The right-hand side $f$ in (1) satisfies that there exists a constant
$M_{f}>0$ such that the following bound holds

$\|f(y,u)\|_{\infty}\leq M_{f},\quad\forall
y\in\overline{\Omega}\subset{\mathbb{R}}^{n},u\in U_{\rm ad},$ (5)

where $\overline{\Omega}$ is a bounded polyhedron such that for sufficiently
small $h>0$ the following inward pointing condition on the dynamics holds

$y+hf(y,u)\in\overline{\Omega},\quad\forall y\in\overline{\Omega},u\in U_{\rm
ad}.$ (6)

  * •

The running cost $g$ is continuous and globally Lipschitz-continuous in both
the first and second arguments; i.e., there exists a constant $L_{g}>0$
satisfying

$\displaystyle|g(y,u)-g(\tilde{y},u)|$ $\displaystyle\leq$ $\displaystyle
L_{g}\|y-\tilde{y}\|,\quad\forall y,\tilde{y}\in{\mathbb{R}}^{n},u\in U_{\rm
ad},$ (7) $\displaystyle|g(y,u)-g(y,\tilde{u})|$ $\displaystyle\leq$
$\displaystyle L_{g}\|u-\tilde{u}\|,\quad\forall u,\tilde{u}\in U_{\rm
ad},y\in{\mathbb{R}}^{n}.$ (8)

  * •

Moreover, there exists a constant $M_{g}>0$ such that

$|g(y,u)|\leq M_{g},\quad\forall(y,u)\in\overline{\Omega}\times U_{\rm ad}.$
(9)

From the assumptions made on $f$ there exists a unique solution of (1)
$y=y(y_{0},u)$ defined on $[0,\infty)$ for every admissible control
$u\in{\mathbb{U}}_{\rm ad}$ and for every initial condition
$y_{0}\in{\mathbb{R}}^{n}$, see [2, Chapter 3]. We define the reduced cost
functional as follows:

$\displaystyle\hat{J}(y_{0},u)=J(y(y_{0},u),u),\quad\forall
u\in{\mathbb{U}}_{\rm ad},\quad y_{0}\in{\mathbb{R}}^{n},$ (10)

where $y(y_{0},u)$ solves (1). Then, the optimal control can be formulated as
follows: for given $y_{0}\in{\mathbb{R}}^{n}$ we consider

$\min_{u\in{\mathbb{U}}_{\rm ad}}\hat{J}(y_{0},u).$

The value function of the problem is defined as
$v:{\mathbb{R}}^{n}\rightarrow{\mathbb{R}}$ as follows:

$v(y)=\inf\left\\{\hat{J}(y,u)\mid u\in{\mathbb{U}}_{\rm ad}\right\\},\quad
y\in{\mathbb{R}}^{n}.$ (11)

This function gives the best value for every initial condition, given the set
of admissible controls $U_{\rm ad}$. It is characterized as the viscosity
solution of the HJB equation corresponding to the infinite horizon optimal
control problem:

$\lambda v(y)+\sup_{u\in U_{\rm ad}}\left\\{-f(y,u)\cdot\nabla
v(y)-g(y,u)\right\\}=0,\quad y\in{\mathbb{R}}^{n}.$ (12)

The solution of (12) is unique for sufficiently large $\lambda$,
$\lambda>\max(L_{g},L_{f})$, [2].

Let us consider first a time discretization where $h$ is a strictly positive
step size. We consider the following semidiscrete scheme for (12):

$v_{h}(y)=\min_{u\in U_{\rm ad}}\left\\{(1-\lambda
h)v_{h}(y+hf(y,u))+hg(y,u)\right\\},\quad y\in{\mathbb{R}}^{n}.$ (13)

The following convergence result for the semidiscrete approximation [6,
Theorem 2.3] requires that for
$(y,\tilde{y},u)\in{\mathbb{R}}^{n}\times{\mathbb{R}}^{n}\times U_{\rm ad}$

$\displaystyle\|f(y+\tilde{y},u)-2f(y,u)+f(y-\tilde{y},u)\|$
$\displaystyle\leq$ $\displaystyle C_{f}\|\tilde{y}\|^{2},$ (14)
$\displaystyle\|g(y+\tilde{y},u)-2g(y,u)+g(y-\tilde{y},u)\|$
$\displaystyle\leq$ $\displaystyle C_{g}\|\tilde{y}\|^{2}.$ (15)

###### Theorem 1

Let assumptions (3), (5), (6), (7), (9), (14) and (15) hold and let
$\lambda>\max(2L_{g},L_{f})$. Let $v$ and $v_{h}$ be the solutions of (12) and
(13), respectively. Then, there exists a constant $C\geq 0$, that can be
bounded explicitly, such that the following bound holds

$\sup_{y\in{\mathbb{R}}^{n}}|v(y)-v_{h}(y)|\leq Ch,\quad h\in[0,1/\lambda).$
(16)

As in [1] let us suppose that there exists a bounded polyhedron
$\Omega\subset{\mathbb{R}}^{n}$ such that for $h>0$ small enough (6) holds. We
consider a fully-discrete approximation to (12). Let
$\left\\{S_{j}\right\\}_{j=1}^{m_{s}}$ be a family of simplices which defines
a regular triangulation of $\Omega$

$\overline{\Omega}=\bigcup_{j=1}^{m_{s}}S_{j},\quad k=\max_{1\leq j\leq
m_{s}}({\rm diam}\ S_{j}).$

We assume we have $n_{s}$ vertices/nodes $y^{1},\ldots,y^{n_{s}}$ in the
triangulation. Let $V^{k}$ be the space of piecewise affine functions from
$\overline{\Omega}$ to ${\mathbb{R}}$ which are continuous in
$\overline{\Omega}$ having constant gradients in the interior of any simplex
$S_{j}$ of the triangulation. Then, a fully discrete scheme for the HJB
equations is given by

$v_{h,k}(y^{i})=\min_{u\in{U}_{\rm ad}}\left\\{(1-\lambda
h)v_{h,k}(y^{i}+hf(y^{i},u))+hg(y^{i},u)\right\\},$ (17)

for any vertex $y^{i}\in\overline{\Omega}$. There exists a unique solution of
(17) in the space $V^{k}$, see [2, Theorem 1.1, Appendix A].

For the fully discrete method if we assume that the controls are Lipschitz-
continuous; i.e., there exists a positive constant $L_{u}>0$ such that

$\displaystyle\|u(t)-u(s)\|_{2}\leq L_{u}|t-s|,$ (18)

then first order of convergence both in time and space is proved in [9,
Theorem 6].

###### Theorem 2

Assume conditions (3), (4), (5), (7), (8), (9) and (18) hold. Assume
$\lambda>\overline{L}$ with $\overline{L}=CnL_{f}$. Then, for $0\leq h\leq
1/(2\lambda)$ there exist positive constants
${C_{1}=C_{1}(\lambda,M_{f},M_{g},L_{f},L_{g})}$ and
${C_{2}=C_{2}(\lambda,L_{f},L_{g},L_{u})}$ such that

$|v(y)-v_{h,k}(y)|\leq C_{1}(h+k)+{C_{2}h},\quad y\in\overline{\Omega}.$

Condition (18) can be weakened and one can still get convergence as proved in
[9, Theorem 7] Assume

  * •

(CA) For every $y\in{\mathbb{R}}^{n}$,

$\displaystyle\left\\{f(y,u),g(y,u),\quad u\in U_{\rm ad}\right\\}$

is a convex subset of ${\mathbb{R}}^{n+1}$.

###### Theorem 3

Assume conditions (3), (4), (5), (7), (8), (9) and (CA) hold. Assume
$\lambda>\overline{L}$ with $\overline{L}$ defined as in Theorem 2. Then, for
$0\leq h\leq 1/(2\lambda)$ there exist positive constants
$C_{1}=C_{1}(\lambda,M_{f},M_{g},L_{f},L_{g})$ and
$C_{2}=C_{2}(\lambda,M_{f},M_{g},L_{f},L_{g})$ such that for
$y\in\overline{\Omega}$

$|v(y)-v_{h,k}(y)|\leq
C_{1}(h+k)+C_{2}\frac{1}{(1+\beta)^{2}\lambda^{2}}(\log(h))^{2}h^{\frac{1}{1+\beta}},\quad\beta=\frac{\sqrt{n}L_{f}}{\lambda}.$
(19)

Let us observe that since $\beta$ is smaller than 1, by weakening the
regularity requirements we loose at most half an order in the rate of
convergence in time of the method up to a logarithmic term.

##  3 POD approximation of the optimal control problem based on time
derivatives

In this section we present a new approach, similar to the second method in
[1], but with snapshots based on time derivatives at different times. We also
perform a completely different error analysis to the one appearing in [1],
inspired in the results in [9] and [10].

###  3.1 POD approximation based on time derivatives

For $p\in\mathbb{N}$ let us choose different pairs
$\left\\{(u^{\nu},y_{0}^{\nu})\right\\}_{\nu=1}^{p}$ in
${\mathbb{U}}\times\overline{\Omega}$. By $y^{\nu}=y(u^{\nu};y_{0}^{\nu})$,
$\nu=1,\ldots,p$, we denote the solutions of (1) corresponding to those chosen
initial conditions and controls.

Let us fix $T>0$ and $M>0$ and take $\Delta t=T/M$ and $t_{j}=j\Delta t$,
$j=0,\ldots,M$. For $N=M+1$ we define the following space

${\cal\boldsymbol{U}}=\mbox{span}\left\\{z_{1}^{\nu},z_{2}^{\nu},\ldots,z_{N}^{\nu}\right\\}_{\nu=1}^{p},$

with

$\displaystyle z_{1}^{\nu}$ $\displaystyle=$
$\displaystyle\sqrt{N}\overline{y}^{\nu},\quad\overline{y}^{\nu}=\frac{1}{N}\sum_{j=0}^{M}y^{\nu}(t_{j})$
$\displaystyle z_{j}^{\nu}$ $\displaystyle=$ $\displaystyle\tau
y_{t}^{\nu}(t_{j-1}),\ j=2,\ldots,N,$

so that

${\cal\boldsymbol{U}}=\mbox{span}\left\\{\sqrt{N}\overline{y}^{\nu},\tau
y_{t}^{\nu}(t_{1}),\ldots,\tau y_{t}^{\nu}(t_{N})\right\\}_{\nu=1}^{p},$

where the factor $\tau$ in front of the temporal derivatives is a time scale
and it makes the snapshots dimensionally correct. The correlation matrix
corresponding to the snapshots is given by
${K}=((k_{i,j}))\in{\mathbb{R}}^{pN\times pN}$, with the entries

$k_{i,j}=\frac{1}{pN}\left(z_{k}^{i},z_{l}^{j}\right),\quad
k,l=1,\ldots,N,\quad i,j=1,\ldots,p.$

and where here, and in the sequel, $(\cdot,\cdot)$ denotes the inner product
in ${\mathbb{R}}^{n}$ to which the norm $||\cdot||$ is associated. Let us
denote for simplicity

${\cal\boldsymbol{U}}=\mbox{span}\left\\{w_{1},w_{2},\ldots,w_{pN}\right\\}:=\left\\{z_{1}^{1},\ldots
z_{N}^{1},\ldots,z_{1}^{p},\ldots,z_{N}^{p}\right\\}.$

Following [11], we denote by
$\lambda_{1}\geq\lambda_{2},\ldots\geq\lambda_{d}>0$ the positive eigenvalues
of $K$ and by
$\boldsymbol{v}_{1},\ldots,\boldsymbol{v}_{d}\in{\mathbb{R}}^{pN}$ its
associated eigenvectors of Euclidean norm $1$. Then, the (orthonormal) POD
basis functions of $\cal\boldsymbol{U}$ are given by

$\varphi_{k}=\frac{1}{\sqrt{pN}}\frac{1}{\sqrt{\lambda_{k}}}\sum_{j=1}^{pN}v_{k}^{j}w_{j},\
k=1,\ldots,d,$ (20)

where $v_{k}^{j}$ is the $j$-th component of the eigenvector
$\boldsymbol{v}_{k}$. The following error estimate is known from [11,
Proposition 1]

$\displaystyle\frac{1}{pN}\sum_{j=1}^{pN}\left\|w_{j}-\sum_{k=1}^{r}(w_{j},\varphi_{k})\varphi_{k}\right\|^{2}=\sum_{k=r+1}^{d}\lambda_{k},$
(21)

from which one can deduce for $\nu=1,\ldots,p$

$\left\|\overline{y}^{\nu}-\sum_{k=1}^{r}(\overline{y}^{\nu},\varphi_{k})\varphi_{k}\right\|^{2}+\frac{\tau^{2}}{M+1}\sum_{j=1}^{M}\left\|y_{t}^{\nu}(t_{j})-\sum_{k=1}^{r}(y_{t}^{\nu}(t_{j}),\varphi_{k})\varphi_{k}\right\|^{2}\leq
p\sum_{k=r+1}^{d}\lambda_{k}.$ (22)

In the sequel, we will denote by

$\displaystyle{\cal\boldsymbol{U}}^{r}=\mbox{span}\\{\varphi_{1},\varphi_{2},\ldots,\varphi_{r}\\},\quad
1\leq r\leq d,$ (23)

and by $P^{r}\ :\ {\mathbb{R}}^{n}\to{\cal\boldsymbol{U}}^{r}$, the orthogonal
projection onto ${\cal\boldsymbol{U}}^{r}$. Then (21) can be written as

$\frac{1}{pN}\sum_{j=1}^{pN}\left\|w_{j}-P^{r}w_{j}\right\|^{2}=\sum_{k=r+1}^{d}\lambda_{k}.$

The following lemma in proved in [10, Lemma 3.2].

###### Lemma 1

Let $T>0$, $\Delta t=T/M$, $t^{n}=n\Delta t$, $n=0,1,\ldots M$, let $X$ be a
Banach space, $\boldsymbol{z}\in H^{2}(0,T;X)$. Then, the following estimate
holds

$\max_{0\leq k\leq
N}\|\boldsymbol{z}^{k}\|_{X}^{2}\leq{3}\|\overline{\boldsymbol{z}}\|_{X}^{2}+\frac{12T^{2}}{M}\sum_{n=1}^{M}\|\boldsymbol{z}_{t}^{n}\|_{X}^{2}+\frac{16T}{3}(\Delta
t)^{2}\int_{0}^{T}\|\boldsymbol{z}_{tt}(s)\|_{X}^{2}\ ds,$ (24)

where
$\overline{\boldsymbol{z}}=\frac{1}{M+1}\sum_{j=0}^{M}\boldsymbol{z}^{j}$.

Using Lemma 1 we can prove pointwise estimates for the projections onto
${\cal\boldsymbol{U}}^{r}$.

###### Lemma 2

The following bounds hold for $\nu=1,\ldots,p$

$\max_{0\leq j\leq
M}\|y^{\nu}(t_{j})-P^{r}y^{\nu}(t_{j})\|^{2}\leq\left(3+24\frac{T^{2}}{\tau^{2}}\right)p\sum_{k={r+1}}^{d}\lambda_{k}+\frac{16T}{3}(\Delta
t)^{2}\int_{0}^{T}\|y^{\nu}_{tt}(s)\|^{2}\ ds.$ (25)

###### Proof

We argue as in [10, Lemma 3.4]. Taking
$\boldsymbol{z}=y^{\nu}(t_{j})-P^{r}y^{\nu}(t_{j})$ in (24) and applying (22)
(taking into account that $(M+1)/M\leq 2$) yields

$\displaystyle\max_{0\leq n\leq M}\|y^{\nu}(t_{j})-P^{r}y^{\nu}(t_{j})\|^{2}$
$\displaystyle\leq$
$\displaystyle\left(3+24\frac{T^{2}}{\tau^{2}}\right)p\sum_{k={r+1}}^{d}\lambda_{k}$
$\displaystyle{}+\frac{16T}{3}(\Delta
t)^{2}\int_{0}^{T}\|y^{\nu}_{tt}(s)-P^{r}y^{\nu}_{tt}(s)\|^{2}\ ds.$

Now, since $P^{r}$ is an orthogonal projection, we have
$\|y^{\nu}_{tt}(s)-P^{r}y^{\nu}_{tt}(s)\|^{2}\leq\|y^{\nu}_{tt}(s)\|^{2}$ and
the proof is finished $\Box$

###  3.2 The POD control problem

To mitigate the curse of dimensionality, the idea of the POD method is to work
on a space of dimension $r$ with $r<n$. To start we need to introduce some
notation.

For any $y\in\overline{\Omega}\subset{\mathbb{R}}^{n}$ let us denote by
$P_{c}^{r}y\in{\mathbb{R}}^{r}$ the coefficients of the projection of $y$ onto
${\cal\boldsymbol{U}}^{r}$

$\displaystyle P_{c}^{r}y=\left\\{(y,\varphi_{k})\right\\}_{k=1}^{r}.$ (26)

For any $y^{r}\in{\mathbb{R}}^{r}$ let us denote by
$\boldsymbol{\varphi}y^{r}\in{\mathbb{R}}^{n}$ the vector whose coefficients
in the POD basis are the components of $y^{r}$, i.e.,

$\displaystyle\boldsymbol{\varphi}y^{r}=\sum_{j=1}^{r}y^{r}_{j}\varphi_{j},$
(27)

where $y^{r}_{j}$ is the $j$ component of the vector $y^{r}$.

For $f$ and $g$ in (1), (2) and $(y^{r},u)\in{\mathbb{R}}^{r}\times U_{\rm
ad}$ we define

$\displaystyle f^{r}(y^{r},u)$ $\displaystyle=$ $\displaystyle
P^{r}_{c}f(\boldsymbol{\varphi}y^{r},u)\in{\mathbb{R}}^{r},$ (28)
$\displaystyle g^{r}(y^{r},u)$ $\displaystyle=$ $\displaystyle
g(\boldsymbol{\varphi}y^{r},u)\in{\mathbb{R}}.$

To have an inward pointing condition on the dynamics in the reduced space,
analogous to (6), following [1, Section 4.2], we assume that there exists a
bounded polyhedron $\overline{\Omega}^{r}\subset{\mathbb{R}}^{r}$ satisfying

$P_{c}^{r}y\in\Omega^{r},\quad\forall y\in\overline{\Omega}.$ (29)

The following lemma proves that the inward pointing condition for $\Omega^{r}$
follows from (29).

###### Lemma 3

Condition (29) implies that

$y^{r}+hf^{r}(y^{r},u)\in\overline{\Omega}^{r},\quad y^{r}=P_{c}^{r}y,\
y\in\overline{\Omega},$

provided the step size $h$ or $\|P^{r}y-y\|$ is sufficiently small.

###### Proof

We follow [1, Remark 4.5] for the proof. We first observe that

$y^{r}+hf^{r}(y^{r},u)=P_{c}^{r}y+hP^{r}_{c}f(\boldsymbol{\varphi}y^{r},u).$

Adding and subtracting $hP^{r}_{c}f(y,u)$ we get

$\displaystyle
y^{r}+hf^{r}(y^{r},u)=P_{c}^{r}(y+hf(y,u))+hP_{c}^{r}(f(\boldsymbol{\varphi}y^{r},u)-f(y,u)).$
(30)

Applying condition (6) $y+hf(y,u)\in\overline{\Omega}$ and applying (29) the
first term on the right-hand side of (30) verifies
$P_{c}^{r}(y+hf(y,u))\in\Omega^{r}$. Then, we only need to show that the
second term on the right-hand side of (30) is small enough for $h$ or
$\|P^{r}y-y\|$ sufficiently small.

Let us denote by $z=f(\boldsymbol{\varphi}y^{r},u)-f(y,u)\in{\mathbb{R}}^{n}$.
Since
$P_{c}^{r}z=\left\\{(z,\varphi_{k})_{X}\right\\}_{k=1}^{r}\in{\mathbb{R}}^{r}$,
taking into account that the functions $\varphi_{k}$ define an orthonormal
basis and that $P^{r}$ is a projection, we have

$\|P_{c}^{r}z\|_{2}=\|P^{r}z\|\leq\|z\|.$

Applying the above inequality together with (3), we get

$\displaystyle\|P_{c}^{r}z\|_{2}^{2}$ $\displaystyle=$
$\displaystyle\|P_{c}^{r}(f(\boldsymbol{\varphi}y^{r},u)-f(y,u))\|_{2}^{2}\leq\|f(\boldsymbol{\varphi}y^{r},u)-f(y,u)\|^{2}$
(31) $\displaystyle\leq$ $\displaystyle
L_{f}^{2}\|\boldsymbol{\varphi}y^{r}-y\|^{2}=L_{f}^{2}\|P^{r}y-y\|^{2},$

so that the proof is concluded. $\Box$

We can now define the reduced order problem we solve in practice. For $f^{r}$
and $g^{r}$ defined in (28) and a given initial condition
$y_{0}^{r}\in{\mathbb{R}}^{r}$ let us consider the controlled nonlinear
dynamical system

$\dot{y}^{r}(t)=f^{r}(y^{r}(t),u(t))\in{\mathbb{R}}^{r},\quad t>0,\quad
y^{r}(0)=y_{0}^{r}\in{\mathbb{R}}^{r},$ (32)

together with the infinite horizon cost functional

$J^{r}(y^{r},u)=\int_{0}^{\infty}g^{r}(y^{r}(t),u(t))e^{-\lambda t}~{}dt.$
(33)

As in (10), we define the reduced cost functional

$\displaystyle\hat{J}^{r}(y_{0}^{r},u)=J^{r}(y^{r}(y_{0}^{r},u),u),\quad\forall
u\in{\mathbb{U}}_{\rm ad},\quad y_{0}^{r}\in{\mathbb{R}}^{r},$ (34)

where $y^{r}(y_{0}^{r},u)$ solves (32). Then, the POD optimal control can be
formulated as follows: for given $y_{0}^{r}\in{\mathbb{R}}^{r}$ we consider

$\min_{u\in{\mathbb{U}}_{\rm ad}}\hat{J}^{r}(y_{0}^{r},u).$

The value function of the problem
$v^{r}:{\mathbb{R}}^{r}\rightarrow{\mathbb{R}}$ is defined as follows:

$v^{r}(y^{r})=\inf\left\\{\hat{J}(y^{r},u)\mid u\in{\mathbb{U}}_{\rm
ad}\right\\},\quad y^{r}\in{\mathbb{R}}^{r}.$ (35)

Remark 1 It is easy to check that the regularity assumptions for $f^{r}$ and
$g^{r}$ analogous to those for $f$ and $g$, (3), (4), (5), (7), (8) and (9),
hold from the definition of $f^{r}$ and $g^{r}$ and the properties being true
for $f$ and $g$.

To get in practice a fully discrete approximation in the reduced space let us
define $\left\\{S_{j}^{r}\right\\}_{j=1}^{m_{s}^{r}}$ a family of simplices
which defines a regular triangulation of $\Omega^{r}$. We assume we have
$n_{s}$ vertices/nodes in the triangulation
$y^{1}_{r},\ldots,y^{n_{s}}_{r}\in\overline{\Omega}^{r}$ and

$\overline{\Omega}^{r}=\bigcup_{j=1}^{m_{s}^{r}}S_{j}^{r},\quad
k_{r}=\max_{1\leq j\leq m_{s}^{r}}({\rm diam}\ S_{j}^{r}).$

Let $V^{k_{r}}$ be the space of piecewise affine functions from
$\overline{\Omega}^{r}$ to ${\mathbb{R}}$ which are continuous in
$\overline{\Omega}^{r}$ having constant gradients in the interior of any
simplex $S_{j}^{r}$ of the triangulation. As in [1, (4.15)] we introduce the
following POD fully discrete scheme for the HJB equations

$v_{h,k}^{r}(y^{i}_{r})=\min_{u\in{U}_{\rm ad}}\left\\{(1-\lambda
h)v_{h,k}^{r}(y^{i}_{r}+hf^{r}(y^{i}_{r},u))+hg^{r}(y^{i}_{r},u)\right\\},\
i=1,\ldots,n_{s},$ (36)

for any vertex $y^{i}_{r}\in\overline{\Omega}^{r}$. As in (17), there exists a
unique solution of (36) in the space $V^{k_{r}}$ defined by its nodal values
(36), see [2, Theorem 1.1, Appendix A].

The key point to carry out the error analysis is that (36) is the fully
discrete approximation to the continuous problem with value function defined
in (35). Moreover, we can apply Theorems 2 and 3 with $v$ and $v_{h,k}$
replaced by $v^{r}$ and $v^{r}_{h,k}$.

As in [1], for any node $y^{i}_{r}\in\overline{\Omega}_{r}$ we set

$y^{i}=\boldsymbol{\varphi}y^{i}_{r},\quad i=1,\ldots n_{s},$

and define

$\tilde{v}_{h,k}^{r}(y)=v_{h,k}^{r}(P^{r}_{c}y),\quad\forall
y\in\overline{\Omega}.$ (37)

For $y^{i}$, $i=1,\ldots n_{s}$, by definition, we have

$\displaystyle\tilde{v}_{h,k}^{r}(y^{i})=v_{h,k}^{r}(P^{r}_{c}y^{i})=v_{h,k}^{r}(y^{i}_{r}),$

since $y^{i}_{r}\in{\mathbb{R}}^{r}$ are the coordinates of
$y^{i}\in{\mathbb{R}}^{n}$ respect to the basis functions of
${\cal\boldsymbol{U}}^{r}$ (23), see (26), (27).

Taking into account that
$y^{i}_{r}+hf^{r}(y^{i}_{r},u)$=$P_{c}^{r}(y^{i}+hf(y^{i},u))$ and
$g^{r}(y_{r}^{i},u)=g(y^{i},u)$ then (36) can also be written as (see [1,
(4.17)])

$\tilde{v}_{h,k}^{r}(y^{i})=\min_{u\in{U}_{\rm ad}}\left\\{(1-\lambda
h)\tilde{v}_{h,k}^{r}(y^{i}+hf(y^{i},u))+hg(y^{i},u)\right\\},\
i=1,\ldots,n_{s}.$

Nevertheless, we do not use the above characterization of the fully discrete
value function to bound the error.

###  3.3 Error analysis of the method

To prove the main results of the paper we need a previous lemma that we now
state and prove. Next lemma bounds the difference between the value function
solving the original problem (11) and the value function solving the reduced
order problem (35).

###### Lemma 4

Let $v$ and $v^{r}$ be the solutions of (11) and (35), respectively. For
$y\in\overline{\Omega}$, let $P^{r}y\in{\mathbb{R}}^{n}$ be the projection of
$y$ onto ${\cal\boldsymbol{U}}^{r}$ and let $P^{r}_{c}y\in{\mathbb{R}}^{r}$ be
as defined in (26). Then, the following bound holds

$\displaystyle|v(P^{r}y)-v^{r}(P^{r}_{c}y)|$ $\displaystyle\leq$
$\displaystyle\ L_{g}\int_{0}^{\infty}te^{(L_{f}-\lambda)t}\max_{0\leq s\leq
t}\|(I-P^{r})f(y(s),u(s))\|dt$
$\displaystyle\quad+L_{g}\int_{0}^{\infty}te^{(L_{f}-\lambda)t}\max_{0\leq
s\leq t}\|(I-P^{r})f(y_{r}(s),u^{r}(s))\|dt,$

where $u,u^{r}:[0,\infty)\rightarrow{\mathbb{R}}^{m}$ are the controls giving
the minimum in (11) and (35) (with initial conditions $P^{r}y$ and
$P^{r}_{c}y$, respectively), $y(t)$ is the solution of (1) with $y(0)=P^{r}y$
and control $u(t)$ and $y_{r}(t)$ is the solution of (1) with $y(0)=P^{r}y$
and control $u^{r}(t)$.

###### Proof

We argue as in [9, Lemma 2]. Let $w:[0,\infty)\rightarrow{\mathbb{R}}^{m}$ be
a given control that we now fix and let $y(t)$ be the solution of (1) with
$y(0)=P^{r}y$ and control $w(t)$. Then

$y(t)=P^{r}y+\int_{0}^{t}f(y(s),w(s))\,ds.$ (39)

Let $y^{r}(t)$ be the solution of (32) with control $w(t)$ and
$y^{r}(0)=P^{r}_{c}y$. Then, recalling definitions (26) and (27) from which we
obtain $\boldsymbol{\varphi}P^{r}_{c}y=P^{r}y$, we can write

$\boldsymbol{\varphi}y^{r}(t)=P^{r}y+\int_{0}^{t}\boldsymbol{\varphi}P^{r}_{c}f(\boldsymbol{\varphi}y^{r}(s),w(s))\,ds.$
(40)

Subtracting (40) from (39) we get

$\displaystyle y(t)-\boldsymbol{\varphi}y^{r}(t)$ $\displaystyle=$
$\displaystyle\int_{0}^{t}(f(y(s),w(s))-\boldsymbol{\varphi}P^{r}_{c}f(\boldsymbol{\varphi}y^{r}(s),w(s)))\,ds$
$\displaystyle=$
$\displaystyle\int_{0}^{t}(f(y(s),w(s))-\boldsymbol{\varphi}P^{r}_{c}f(y(s),w(s)))\,ds$
$\displaystyle\quad+\int_{0}^{t}(\boldsymbol{\varphi}P^{r}_{c}f(y(s),w(s))-\boldsymbol{\varphi}P^{r}_{c}f(\boldsymbol{\varphi}y^{r}(s),w(s)))\,ds.$

Taking norms

$\displaystyle\|y(t)-\boldsymbol{\varphi}y^{r}(t)\|$ $\displaystyle\leq$
$\displaystyle\int_{0}^{t}\|(I-P^{r})f(y(s),w(s))\|\,ds$
$\displaystyle\quad+\int_{0}^{t}\|\boldsymbol{\varphi}P^{r}_{c}f(y(s),w(s))-\boldsymbol{\varphi}P^{r}_{c}f(\boldsymbol{\varphi}y^{r}(s),w(s))\|\,ds.$

Using again that $\boldsymbol{\varphi}P^{r}_{c}=P^{r}$ and applying that
$P^{r}$ is a projection together with (3) we get

$\displaystyle\|\boldsymbol{\varphi}P^{r}_{c}f(y(s),w(s))-\boldsymbol{\varphi}P^{r}_{c}f(\boldsymbol{\varphi}y^{r}(s),w(s))\|^{2}$
$\displaystyle\quad=\|P^{r}\left(f(y(s),w(s))-f(\boldsymbol{\varphi}y^{r}(s),w(s))\right)\|^{2}$
$\displaystyle\quad\leq\|f(y(s),w(s))-f(\boldsymbol{\varphi}y^{r}(s),w(s)\|^{2}\leq
L_{f}^{2}\|y(s)-\boldsymbol{\varphi}y^{r}(s)\|^{2}.$

Going back to (Proof)

$\|y(t)-\boldsymbol{\varphi}y^{r}(t)\|\leq
L_{f}\int_{0}^{t}\|y(s)-\boldsymbol{\varphi}y^{r}(s)\|ds+t\max_{0\leq s\leq
t}\|(I-P^{r})f(y(s),w(s))\|.$

Applying Gronwall lemma we get

$\displaystyle\|y(t)-\boldsymbol{\varphi}y^{r}(t)\|\leq{e^{L_{f}t}}\left(t\max_{0\leq
s\leq t}\|(I-P^{r})f(y(s),w(s))\|\right).$ (42)

We now observe that from definitions (10) and (34) we get

$\displaystyle|\hat{J}(P^{r}y,w)-\hat{J}^{r}(P^{r}_{c}y,w)|\leq\int_{0}^{\infty}|g(y(t),w(t))-g(\boldsymbol{\varphi}y^{r}(t),w(t))|e^{-\lambda
t}dt.$

Applying then the Lipschitz-continuity property of $g$, (7), together with
(42) we get

$\displaystyle|\hat{J}(P^{r}y,w)-\hat{J}^{r}(P^{r}_{c}y,w)|\leq$
$\displaystyle\quad L_{g}\int_{0}^{\infty}t{e^{(L_{f}-\lambda)t}}\max_{0\leq
s\leq t}\|(I-P^{r})f(y(s),w(s))\|\ dt.$ (43)

To conclude we will argue similarly as in in [9, Theorem 4].

For any $y\in\overline{\Omega}$, let
$u^{r}:[0,\infty)\rightarrow{\mathbb{R}}^{m}$ be the control giving the
minimum in (35) with initial condition $P^{r}_{c}y$. Since by definition of
$v$, $v(P^{r}y)\leq\hat{J}(P^{r}y,u^{r})$ and
$v^{r}(P^{r}_{c}y)=\hat{J}^{r}(P^{r}_{c}y,u^{r})$, applying (Proof) with
$w=u^{r}$ we get

$\displaystyle v(P^{r}y)-v^{r}(P^{r}_{c}y)$ $\displaystyle\leq$
$\displaystyle\hat{J}(P^{r}y,u^{r})-\hat{J}^{r}(P^{r}_{c}y,u^{r})$ (44)
$\displaystyle\leq$ $\displaystyle
L_{g}\int_{0}^{\infty}te^{(L_{f}-\lambda)t}\max_{0\leq s\leq
t}\|(I-P^{r})f(y_{r}(s),u^{r}(s))\|\ dt,$

where $y_{r}(t)$ is the solution of (1) with $y(0)=P^{r}y$ and control
$u^{r}(t)$.

On the other hand, let $u:[0,\infty)\rightarrow{\mathbb{R}}^{m}$ be the
control giving the minimum in (11) with initial condition $P^{r}y$. Arguing as
in (44) and applying (Proof) again, with $w=u$, we get

$\displaystyle v^{r}(P^{r}_{c}y)-v(P^{r}y)$ $\displaystyle\leq$
$\displaystyle\hat{J}^{r}(P^{r}_{c}y,u)-\hat{J}(P^{r}y,u)$ (45)
$\displaystyle\leq$ $\displaystyle
L_{g}\int_{0}^{\infty}te^{(L_{f}-\lambda)t}\max_{0\leq s\leq
t}\|(I-P^{r})f(y(s),u(s))\|\ dt,$

where $y(t)$ is the solution of (1) with $y(0)=P^{r}y$ and control $u$. From
(44) and (45) we conclude (4). $\Box$

In next theorem we bound the difference between the value function of the
original problem and our fully discrete approximation based on POD. Let
$u^{r}$ be the control giving the minimum in (35). For the proof of next
theorem we need to assume that there exists a positive constant $L_{u^{r}}>0$
such that

$\displaystyle\|u^{r}(t)-u^{r}(s)\|_{2}\leq L_{u^{r}}|t-s|.$ (46)

###### Theorem 4

Let $v$ be the solution of (11) and let $\tilde{v}_{h,k}^{r}$ be its fully
discrete POD approximation defined in (36)-(37). Assume conditions (3), (4),
(5), (7), (8), (9), (14), (15) and (46) hold. Assume
$\lambda>\max(2L_{g},L_{f},\overline{L}_{r})$ with $\overline{L}_{r}=CrL_{f}$.
Then, for any $y\in\overline{\Omega}$ and $0\leq h\leq 1/(2\lambda)$ there
exist positive constants $C_{1}$ and $C_{2}$ such that the following bound
holds for $y\in\overline{\Omega}$

$\displaystyle|v(y)-\tilde{v}_{h,k}^{r}(y)|$ $\displaystyle\leq$
$\displaystyle C_{1}(h+k_{r})+{C_{2}h}+\frac{L_{g}}{\lambda-
L_{f}}\|y-P^{r}y\|$ $\displaystyle\
+L_{g}\int_{0}^{\infty}te^{(L_{f}-\lambda)t}\max_{0\leq s\leq
t}\|(I-P^{r})f(y(s),u(s))\|dt$ $\displaystyle\
+L_{g}\int_{0}^{\infty}te^{(L_{f}-\lambda)t}\max_{0\leq s\leq
t}\|(I-P^{r})f(y_{r}(s),u^{r}(s))\|dt,$

where $u,u^{r}:[0,\infty)\rightarrow{\mathbb{R}}^{m}$ are the controls giving
the minimum in (11) and (35) (with initial conditions $P^{r}y$ and
$P^{r}_{c}y$, respectively), $y(t)$ is the solution of (1) with $y(0)=P^{r}y$
and control $u(t)$ and $y_{r}(t)$ is the solution of (1) with $y(0)=P^{r}y$
and control $u^{r}(t)$.

###### Proof

We first observe that adding and subtracting terms we get

$\displaystyle|v(y)-\tilde{v}_{h,k}^{r}(y)|$ $\displaystyle\leq$
$\displaystyle|v(y)-v(P^{r}y)|+|v(P^{r}y)-v^{r}(P^{r}_{c}y)|+|v^{r}(P^{r}_{c}y)-\tilde{v}_{h,k}^{r}(y)|.$

To bound the first term we write

$|v(y)-v(P^{r}y)|\leq|v(y)-v_{h}(y)|+|v_{h}(y)-v_{h}(P^{r}y)|+|v_{h}(P^{r}y)-v(P^{r}(y)|,$

and then apply (16) to the first and third terms and the Lipschitz-continuity
of $v_{h}$ that holds for $\lambda>L_{f}$ (see [7, p. 473]) to the second
term. Then

$|v(y)-v(P^{r}y)|\leq Ch+\frac{L_{g}}{\lambda-L_{f}}\|y-P^{r}y\|.$

To bound the second term we apply Lemma 4. To conclude we need to bound the
third term. To this end we observe that

$v^{r}(P^{r}_{c}y)-\tilde{v}_{h,k}^{r}(y)=v^{r}(P^{r}_{c}y)-v_{h,k}^{r}(P^{r}_{c}y)$

so that we can apply Theorem 2 to $v^{r}$ and $v^{r}_{h,k}$ to reach (4).
$\Box$

Remark 2 Let us observe that the first two terms on the right-hand side of (4)
give the rate of convergence of the method in terms of the time step $h$ and
mesh diameter $k_{r}$. The other three terms come from the POD approximation
and will be bounded at the end on this section. These terms depend on the set
of snapshots and the tail of the eigenvalues in the singular value
decomposition.

To apply Theorem 2 in the proof of Theorem 4 we use the properties of $f^{r}$
and $g^{r}$ that, as commented in Remark Proof are inherited from the assumed
hypothesis made on $f$ and $g$, (3), (4), (5), (7), (8), (9). We also need to
assume condition (46) holds for the control of the reduced order problem.
Condition (46) can be weakened and one can still get convergence assuming
instead

  * •

(CAr) For every $y^{r}\in{\mathbb{R}}^{r}$,

$\displaystyle\left\\{f^{r}(y^{r},u),g^{r}(y^{r},u),\quad u\in U_{\rm
ad}\right\\}$

is a convex subset of ${\mathbb{R}}^{r+1}$.

This result is stated in Theorem 5 below. In the proof of the theorem we apply
Theorem 3 instead of Theorem 2.

###### Theorem 5

Let $v$ be the solution of (11) and let $\tilde{v}_{h,k}^{r}$ be its fully
discrete POD approximation defined in (36)-(37). Assume conditions (3), (4),
(5), (7), (8), (9), (14), (15) and (CAr) hold. Assume
$\lambda>\max(2L_{g},L_{f},\overline{L}_{r})$ with $\overline{L}_{r}=CrL_{f}$.
Then, for any $y\in\overline{\Omega}$ and $0\leq h\leq 1/(2\lambda)$ there
exist positive constants $C_{1}$ and $C_{2}$ such that the following bound
holds for $y\in\overline{\Omega}$

$\displaystyle|v(y)-v_{h,k}(y)|$ $\displaystyle\leq$ $\displaystyle
C_{1}(h+k_{r})+C_{2}\frac{1}{(1+\beta)^{2}\lambda^{2}}(\log(h))^{2}h^{\frac{1}{1+\beta}}$
$\displaystyle\ +\frac{L_{g}}{\lambda-L_{f}}\|y-P^{r}y\|$ $\displaystyle\
+L_{g}\int_{0}^{\infty}te^{(L_{f}-\lambda)t}\max_{0\leq s\leq
t}\|(I-P^{r})f(y(s),u(s))\|dt$ $\displaystyle\
+L_{g}\int_{0}^{\infty}te^{(L_{f}-\lambda)t}\max_{0\leq s\leq
t}\|(I-P^{r})f(y_{r}(s),u^{r}(s))\|dt,$

where $\beta=\frac{\sqrt{n}L_{f}}{\lambda}$ and
$u,u^{r}:[0,\infty)\rightarrow{\mathbb{R}}^{m}$ are the controls giving the
minimum in (11) and (35) (with initial conditions $P^{r}y$ and $P^{r}_{c}y$,
respectively), $y(t)$ is the solution of (1) with $y(0)=P^{r}y$ and control
$u(t)$ and $y_{r}(t)$ is the solution of (1) with $y(0)=P^{r}y$ and control
$u^{r}(t)$.

###### Proof

The proof is the same than the proof of Theorem 4 applying Theorem 3 instead
of Theorem 2. $\Box$

Remark 3 The same comments as in Remark Proof apply with the difference that
the rate of convergence in terms of the time step $h$ is reduced due to the
weaker regularity requirements. Since $\beta$ is smaller than $1$ we loose at
most half an order in the rate of convergence in time up to a logarithmic
term.

To conclude we will give an estimation of the last three terms in (4) and (5).
The first term is bounded in the following lemma, where, recall, $p$ is then
number of trajectories in the set of snapshots (see Section 3.1).

###### Lemma 5

For $y\in\overline{\Omega}$ and $P^{r}$ the $X$-orthogonal projection onto
${\cal\boldsymbol{U}}^{r}$ the following bound holds

$\displaystyle\left\|y-P^{r}y\right\|\leq$
$\displaystyle\min_{\genfrac{}{}{0.0pt}{2}{1\leq l\leq p}{0\leq n\leq
M}}\left\|y-y^{l}(t_{n})\right\|+\Biggl{[}\left(3+24\frac{T^{2}}{\tau^{2}}\right)p\sum_{k={r+1}}^{d}\lambda_{k}$
$\displaystyle{}+\frac{16T}{3}(\Delta t)^{2}\max_{1\leq l\leq
p}\int_{0}^{T}\|y^{l}_{tt}(s)-P^{r}y^{l}_{tt}(s)\|^{2}\ ds\Biggr{]}^{1/2}.$
(49)

###### Proof

Let $k$ and $j$ be such that

$\left\|y-y^{j}(t_{k})\right\|=\min_{\genfrac{}{}{0.0pt}{2}{1\leq l\leq
p}{0\leq n\leq M}}\left\|y-y^{l}(t_{n})\right\|.$

Then, we can write

$y-P^{r}y=(I-P^{r})\left(y-y^{j}(t_{k})\right)+y^{j}(t_{k})-P^{r}y^{j}(t_{k}).$
(50)

Now noticing that $\left\|I-P^{r}\right\|\leq 1$ and recalling Lemma 2, from
(50) it follows (5) $\Box$

Remark 4 The first term on the right hand-side in (5) reflects the closeness
of the data $y$ to the set of snapshots while the other term is the projection
error onto the POD basis. Let us observe that using time derivatives in the
set of snapshots allow us to get a bound for this projection error in the
discrete maximum norm in time, see Lemma 2.

We bound the second term on the right-hand side of (4) in the following lemma.

###### Lemma 6

For $y\in\overline{\Omega}$ and $P^{r}$ the orthogonal projection onto
${\cal\boldsymbol{U}}^{r}$ let $u:[0,\infty)\rightarrow{\mathbb{R}}^{m}$ be
the control giving the minimum in (11) with initial condition $P^{r}y$. Let
$y(t)$ be the solution of (1) with $y(0)=P^{r}y$ and control $u(t)$. Then, for
any fixed $s\in[0,\infty)$ the following bound holds

$\displaystyle\|(I-P^{r})f(y(s),u(s))\|$ $\displaystyle\leq$
$\displaystyle\min_{\genfrac{}{}{0.0pt}{2}{1\leq l\leq p}{0\leq n\leq
M}}\|f(y(s),u(s))-f(y^{l}(t_{n}),u^{l}(t_{n}))\|$ (51)
$\displaystyle\quad+\frac{(M+1)}{\tau^{2}}p\sum_{k=r+1}^{d}\lambda_{k}.$

###### Proof

We argue as in the proof of Lemma 5. Let $k$ and $j$ be such that

$\|f(y(s),u(s))-f(y^{j}(t_{k}),u^{j}(t_{k}))\|=\min_{\genfrac{}{}{0.0pt}{2}{1\leq
l\leq p}{0\leq n\leq M}}\|f(y(s),u(s))-f(y^{l}(t_{n}),u^{l}(t_{n}))\|.$

Then,

$\displaystyle(I-P^{r})f(y(s),u(s))=(I-P^{r})(f(y(s),u(s))-f(y^{j}(t_{k}),u^{j}(t_{k})))$
$\displaystyle\quad+(I-P^{r})f(y^{j}(t_{k}),u^{j}(t_{k})),$

so that applying (22) and $\|I-P^{r}\|\leq 1$ we get (51). $\Box$

Remark 5 An error bound for $\|(I-P^{r})f(y_{r}(s),u^{r}(s))\|$ can be obtain
arguing exactly in the same way.

Remark 6 As in (5), the first term on the right-hand side of (51) reflects the
closeness of $f(y(s),u(s))$ to the set of snapshots and the second one is the
projection error onto the POD basis. Let us observe that the use of temporal
derivates in the set of snapshots is essential to get the bound (51).

In case one has a uniform distribution along the discrete times of the errors
in (22), as it is often the case (at least in our experience with numerical
computations concerning POD methods, see [10, Figure 1]), one would expect for
the second term on the right-hand side of (51) a behaviour as
$\frac{p}{\tau^{2}}\sum_{k=r+1}^{d}\lambda_{k}$ instead of the rude bound
$\frac{(M+1)}{\tau^{2}}p\sum_{k=r+1}^{d}\lambda_{k}$.

##  4 Numerical Experiments

We now present some numerical experiments. We closely follow those in [1] so
that the new method we propose can be compared with the method in [1]. We
first notice that due to numerical reasons we have to choose a finite time
horizon, so we select a sufficiently large $t_{e}>0$, which, in the
experiments that follow, it was fixed to $t_{e}=3$. As in [1], we consider the
following convection-reaction-diffusion equation

$\begin{array}[]{rcll}z_{t}-\varepsilon z_{xx}+\gamma
z_{x}+\mu(z^{3}-z)&=&ub&\hbox{\rm in $I\times(0,t_{e})$},\\\
z(\cdot,0)&=&z_{0}&\hbox{\rm in $I$},\\\ z(\cdot,t)&=&0&\hbox{\rm
in~{}$\partial I\times(0,t_{e})$},\end{array}$ (52)

with $\varepsilon=1/10$, and where $I=(0,a)$ is an open interval,
$z:I\times[0,t]\rightarrow{\mathbb{R}}$ denotes the state, and $\gamma$ and
$\mu$ are positive constants. The controls $u$ belong to the closed, convex,
bounded set ${\mathbb{U}}_{\rm ad}=L^{2}(0,t_{e},[u_{a},u_{b}])$, for real
values $u_{a}<u_{b}$. The cost functional to minimize is given by

$\int_{0}^{t_{e}}e^{-\lambda
t}\left(\left\|z(\cdot,t,u)\right\|_{L^{2}(I)}^{2}+\frac{1}{100}\left|u(t)\right|^{2}\right)\,dt,$
(53)

where we set $\lambda=1$. Notice then that in (53) the aim is to drive the
state to zero.

We use a finite-difference method on a uniform grid of size $\Delta x=l/N$
with $N=100$ on the interval $I=(0,l)$ to discretize (52) in space to obtain a
system of ordinary differential equations (ODEs). To obtain the snapshots, the
ODE system is integrated in time using Matlab’s command ode15s, which uses the
numeric differentiation formulae (NDF) [14], with sufficiently small
tolerances for the local errors (below $10^{-12}$). The snapshots were
obtained on a uniform (time) grid of diameter $1/20$. The time derivatives
were obtained by evaluating the right-hand-side of the system of ODEs.

As in [1], equation (36)

$v_{h,k}^{r}(y^{i}_{r})=\min_{u\in{U}_{\rm ad}}\left\\{(1-\lambda
h)v_{h,k}^{r}(y^{i}_{r}+hf^{r}(y^{i}_{r},u))+hg^{r}(y^{i}_{r},u)\right\\},\
i=1,\ldots,n_{s},$

was solved by fixed-point iteration, the stopping criterium being that two
consecutive iterates differ in the maximum norm in less that $5\times
10^{-4}$. For the first iterate we choose a family of constant controls
$u^{l}$, $l=1,\ldots,p$ and at any point of the mesh $y_{r}^{i}$ (initial
condition) we compute the approximate solution of (32) corresponding to this
initial condition and control $u^{l}$. Then, we compute the value of the
functional cost (33). Finally, the value of the initial iterate at $y_{r}^{i}$
is the minimum between the values of the functional cost for $l=1,\ldots,p$.

Once (36) is solved, the optimal control

$u_{h,k}^{r}(y^{i}_{r})=\operatorname*{argmin}_{u\in{U}_{\rm
ad}}\left\\{(1-\lambda
h)v_{h,k}^{r}(y^{i}_{r}+hf^{r}(y^{i}_{r},u))+hg^{r}(y^{i}_{r},u)\right\\},$
(54)

is obtained at any mesh point $y^{i}_{r}$, $i=1,\ldots,n_{s}$. Then, the
suboptimal feedbak operator $\Phi^{r}(y)$ is computed by interpolation. This
means that for $y\in\overline{\Omega}$ we project onto the POD space to get
$P^{r}y$ and then

$P^{r}y=\sum_{i=1}^{n_{s}}\mu_{i}y^{i}_{r},\quad\Phi^{r}(y)=\sum_{i=1}^{n_{s}}\mu_{i}u_{h,k}^{r}(y^{i}_{r}),$

where the coefficients $\mu_{i}$ satisfy $0\leq\mu_{i}\leq 1$,
$\sum_{i=1}^{n_{s}}\mu_{i}=1$.

With this, the closed-loop system

$y^{\prime}(t)=f(y(t),\Phi^{r}(y(t))),\qquad y(0)=y_{0},$ (55)

is integrated, again, using the NDF formulae as implemented in Matlab’s
command ode15s with the same tolerances as in the computation of the
snapshots.

###  4.1 Test 1: Semilinear equation

As in [1], we consider (52) with $\gamma=0$ and $\mu=1$, $l=1$ and
$b(x)=z_{0}(x)=2x(1-x)$. It is easy to check that the uncontrolled solution
converges, as $t\rightarrow\infty$ to a non-null steady state (see also [1,
Fig. 6.1]), and that the null solution is unstable.

For the finite-difference approximation, we consider
$y:[0,t_{e}]\rightarrow{\mathbb{R}}^{N-1}$ with components $y_{j}(t)\approx
z(x_{j},t)$, $x_{j}=j\Delta x$, $j=1,\ldots,N-1$, solution of

$Cy_{t}=\frac{1}{10}Ay+C(F(y)+uB)$ (56)

where the components of $F$ and $B$ are, respectively
$F_{j}=y_{j}(1-y_{j}^{2})$, $B_{j}=2x_{j}(1-x_{j})$, $j=1,\ldots,N-1$, and $A$
and $C$ are $(N-1)\times(N-1)$ tridiagonal matrices matrices given by

$A=\frac{1}{(\Delta x)^{2}}\left[\begin{array}[]{cccccc}-2&\ 1\\\ \ 1&-2&\
1\\\ &\ddots&\ddots&\ddots\\\ &&\ 1&-2&\ 1\\\ &&&\
1&-2\end{array}\right],\qquad
C=\frac{1}{12}\left[\begin{array}[]{cccccc}10&1\\\ 1&10&1\\\
&\ddots&\ddots&\ddots\\\ &&1&10&1\\\ &&&1&10\end{array}\right],$

so that the finite-difference discretization (56) is fourth-order convergent.
The norm we consider in ${\mathbb{R}}^{N-1}$ is given by

$\left\|y\right\|^{2}=\Delta x\sum_{j=1}^{N-1}y_{j}^{2}.$ (57)

Let as observe that this norm is an approximation to the integral
$\int_{0}^{1}y(x)^{2}dx$ of a function with values $y_{j}$ at the spatial mesh
nodes.

To compute the snapshots, as in [1], for constant controls $u\in U_{\rm
snap}=\\{-1,0,1\\}$, we obtained the solutions $y^{(n)}=y(t_{n})$ of (56)
every $1/20$ time units, that is, for $t_{n}=n/20$, $n=1,\ldots 20t_{e}$, and
then the time derivatives $y_{t}^{(n)}$ were computed from identity (56). For
the reduced spaces, we consider the cases of POD basis with only $r=2,3$ and
$4$ elements, also as in [1]. The control set ${\mathbb{U}}_{\rm ad}$ is given
by 21 controls equally distributed in $[-1,1]$.

As in [1], to define the domain $\overline{\Omega}^{r}$, we compute the
projections of all the snapshots. With this procedure we obtain a set of
points in $\mathbb{R}^{r}$. Then, we define an hypercube containing this set
of points. The aim of this procedure, in view of Lemma 3, is that the set
$\overline{\Omega}^{r}$ defined in this way satisfies the invariance condition

$y^{r}+hf^{r}(y^{r},u)\in\overline{\Omega}^{r},\quad
y^{r}\in\overline{\Omega}^{r},\quad u\in{\mathbb{U}}_{\rm ad}.$ (58)

The set $\overline{\Omega}^{r}$ for $r=4$ was given by

$\overline{\Omega}^{r}=[-0.9,0.425]\times[-0.075,0.07]\times[-0.04,0.04]\times[-0.025,0.25].$

For this set we checked that condition (58) holds.

We notice that our set $\overline{\Omega}^{r}$ is considerable smaller than
the corresponding set in [1] (see [1, 6.1. Test 1]) where the authors use the
standard euclidean norm in ${\mathbb{R}}^{n}$ rather than the norm (57) we use
here. Since the domain is smaller we also consider partitions of $\Omega^{r}$
smaller than those in [1]. We take maximum diameter $k_{r}=0.02$ and
$k_{r}=0.01$, and, as in [1], we choose $h=0.1k_{r}$.

In Fig. 1 we have represented on top the optimal solution (left) for $r=4$,
the difference between optimal solution with 4 and 2 POD basis functions (top
middle) and the difference between optimal solution with 4 and 3 POD basis
functions (top right). On bottom we have represented the optimal controls for
$r=2$, $3$ and $4$.

Figure 1:  Test 1: Optimal HJB states computed with r=4 POD basis functions
(top-left), difference between optimal solution with 4 and 2 POD basis
functions (top-middle), difference between optimal solution with 4 and 3 POD
bases (top-right). Optimal HJB controls with r=4,3,2 (bottom).

We observe that we get much better results than those in [1], although (apart
from using a different set of snapshots) in our method, both the finite-
difference method and the time integrator that we use are more accurate than
those in [1].

###  4.2 Test 2: Advection-diffusion equation

As in [1], we now consider (52) with $\gamma=1$ and $\mu=0$, $I=(0,2)$ and
$z_{0}(x)=\max(0,0.5\sin(\pi x))$. We take $b$ as the characteristic function
of the inverval $(1/2,1)$. To compute the POD basis we compute the time
derivatives of the states for constant controls $u=-2.2,-1.1,0$. The
semidiscretization was done with a standard finite difference method

$\frac{dy_{j}}{dt}=\frac{y_{j+1}-2y_{j}+y_{j-1}}{10(\Delta
x)^{2}}-\frac{y_{j+1}-y_{j-1}}{2\Delta x}+b(x_{j}),\quad j=1,\ldots,N-1,\qquad
y_{0}=y_{N}=0,$

which is second order convergent in problems with sufficiently smooth
solutions.

Since the initial state $z_{0}$ does not possess second-order derivatives in
$L^{2}$, we notice then that the time derivative $z_{t}$ blows up when
$t\rightarrow 0$. For this reason, after the spatial discretization by finite
differences, we replaced the time derivative at $t=0$ by the difference
quotient $(y^{(1)}-y^{(0)})/\Delta t$ of states at $t=0$ and $t=\Delta t$,

In this test we found it impossible to find a set $\Omega^{r}$ satisfying
condition (58) both when the POD basis is taken from the states and from their
time derivatives. In this last case the set $\Omega^{r}$ for $r=4$ we used in
the experiments was

$\Omega^{r}=(-0.5,0.7)\times(-0.3,1.5)\times(-0.3,0.2)\times(-0.05,0.15).$

To overcome the lack of invariance of this set, whenever for a vertex
$y^{i}_{r}$ we had
$y^{i}_{r}+hf^{r}(y^{i}_{r},u)\not\in\overline{\Omega}^{r}$, we simply
replaced $y^{i}_{r}+hf^{r}(y^{i}_{r},u)$ by its closest point on
$\partial\Omega^{r}$. This resulted in changing the value of
$y^{i}_{r}+hf^{r}(y^{i}_{r},u)$ in less than 2% in the first two coordinates
in the POD basis and $15\%$ in the remaining ones, except for some negative
values of the fourth coordinate where errors up to $60\%$ were encountered.
For example, an error of $15\%$ in the third coordinate means that for some
values of $y^{i}_{r}+hf^{r}(y^{i}_{r},u)$ the third coordinate could be in the
set $[-0.345,0.23]$ instead of $[-0.3,0.2]$. Nevertheless, as we will see
below, the results obtained with the POD approximation in this test were
excellent.

Since this problem is linear-quadratic, the solution of HJB equation can be
computed by solving Ricatti equation. In Fig. 2 we show the uncontrolled
solution (left), the optimal LQR state (middle) and the optimal LQR control
(right).

Figure 2:  Test 2: Uncontrolled solution (left), the optimal LQR state
(middle) and the optimal LQR control (right).

In Figure 3 we have represented on top the optimal solution (left) for $r=4$,
the difference between optimal solution with 4 and 2 POD basis functions (top
middle) and the difference between optimal solution with 4 and 3 POD basis
functions (top right). On the bottom part we have represented the optimal
controls for $r=2$, $3$ and $4$. Also in this case, the improvement with
respect to the results in [1] is remarkable. In particular, the optimal
controls in Figure 3 compare very well with the optimal LQR control of Figure
2 even for the case with only $r=2$ basis functions in our POD method.

Figure 3:  Test 2: Optimal HJB states computed with $r=4$ POD basis functions
(top-left), difference between optimal solution with 4 and 2 POD basis
functions (top middle), difference between optimal solution with 4 and 3 POD
basis functions (top-right). Optimal HJB controls with $r=4,3,2$ (bottom).

To conclude, in Figure 4 (left) we show the difference between the optimal LQR
state and the optimal HJB state computed with $r=4$ POD basis functions. On
the right, we show the relative errors $\left|u_{\rm
HJB}-u_{LQR}\right|/\max(10^{-3},\left|u_{LQR}\right|)$ of the optimal HJB
controls with respect to the optimal LQR control for $r=2,3$ and $4$. It can
be seen a very good agreement between HJB and LQR optimal states. With respect
to the optimal controls, we notice that whereas with $r=3$ and $r=4$ POD basis
functions the errors do not exceed 30% and, indeed, they stay below 10% most
of the time, this is not the case of $r=2$ POD basis functions, where errors
are above 100% for more than half the time interval. However, let us observe
that we are considering relative errors on the right of Figure 4 and that
restricting ourselves to the time interval in which the optimal control is
sufficiently away from zero the errors for $r=2$ are also below 35%.

Figure 4:  Test 2: difference between the optimal LQR state and the optimal
HJB state computed with $r=4$ POD basis functions (left) and relative errors
$\left|u_{\rm HJB}-u_{LQR}\right|/\max(10^{-3},\left|u_{LQR}\right|)$ of the
optimal HJB controls with respect to the optimal LQR control (right).

Again, the results here (which correspond to $k_{r}=0.1$) compare favourably
with those in the literature.

##  5 Conclusions

In this paper we introduce a reduced order method based on POD to mitigate the
curse of dimensionality in the numerical approximation of HJB equations. The
novelty of the method is the use of snapshots based on temporal derivatives of
the controlled nonlinear dynamical system.

We carry out the error analysis of the method based on the recent results
obtained in [9] that allow us to get sharper error bounds than those appearing
in the literature. In particular, the factor $1/h$ where $h$ is the time step
of the fully discrete method does not appear in our error bounds. Our error
bounds are optimal in terms of the time step $h$ and the mesh diameter of the
reduced space $k_{r}$ and as usual depend also on the size of the tail of
eigenvalues in the singular value decomposition. The use of snapshots based on
time derivatives allow us to give a bound for some of the terms in the error
that could not be bounded with the standard approach.

Some numerical experiments are performed in which we show that the new method
we propose obtains better results than a similar POD method presented in [1].

## References

  * [1] A. Alla, M. Falcone & S. Volkwein, Error analysis for POD approximations of infinite horizon problems via the dynamic programming approach, SIAM J. Control Optim., 55 (2017), 3091–3115.
  * [2] M. Bardi & I. Capuzzo-Dolcetta, Optimal Control and Viscosity Solutions of Hamilton-Jacobi-Belmann Equations, Springer Science+Business Media, LLC, New York, 1997.
  * [3] D. Kalise & K. Kunisch, Polynomial approximation of high-dimensional Hamilton-Jacobi-Bellman equations and applications to feedback control of semilinear parabolic PDEs. SIAM J Sci Comput. 40 (2018), A629-A652.
  * [4] S. Dolgov, D. Kalise & K. Kunisch. Tensor decomposition methods for high-dimensional Hamilton-Jacobi-Bellman Equations. SIAM J Sci Comput. 43 (2021), A1625-A650.
  * [5] M. Eigel, R. Schneider & D. Sommer, Dynamical low-rank approximations of solutions to the Hamilton-Jacobi-Bellman equation, Numer. Linear Algebra Appl., 30 (2023) Paper No. e2463, 20 pp.
  * [6] M. Falcone, A numerical approach to the infinite horizon problem of deterministic control theory, Appl. Math. Optim., 15 (1987), 1–13.
  * [7] M. Falcone, Numerical solution of dynamic programming equations, in Optimal Control and Viscosity Solutions of Hamilton-Jacobi-Bellman Equations, Birkhäuser Boston, Boston, MA, 1997, pp. 471-504
  * [8] M. Falcone, P. Lacunara & Seghini, A splitting algorithm for Hamilton-Jacobi-Bellman equations, Appl. Numer. Math., 15 (1994), 207–218.
  * [9] J. de Frutos & J. Novo, Optimal bounds for numerical approximations of infinite horizon problems based on dynamic programming approach, SIAM. J. Control Optim. 61 (2023) 415–433.
  * [10] B. García-Archilla, V. John & J. Novo, POD-ROMs for incompressible flows including snapshots of the temporal derivative of the full order solution, SIAM J. Numer. Anal. 61 (2023) 1340–1368.
  * [11] K. Kunisch & S. Volkwein, Galerkin proper orthogonal decomposition methods for parabolic problems, Numer. Math., 90 (2001), 117-148.
  * [12] B. Luo, W. Huai-Ning, T. Huang & D. Liu Data-based approximate policy iteration for affine nonlinear continuous-time optimal control design. Automatica, 50 (2014), 3281–3290.
  * [13] D. Tonon, M.S. Aronna & D. Kalise, Optimal control: novel directions and applications. Vol 1. Cham: Springer; 2017.
  * [14] L. F. Shampine & M. W. Reichelt, The Matlab ODE Suite, SIAM J. Sci. Comput., 18 (1997), pp. 1–22.