Optimal Bounds for POD Approximations of Infinite Horizon Control Problems Based on Time Derivatives

Journal of Scientific Computing (2025) 103:19 https://doi.org/10.1007/s10915-025-02833-0 Optimal Bounds for POD Approximations of Infinite Horizon Control Problems Based on Time Derivatives Javier de Frutos1 ·Bosco García-Archilla2 ·Julia Novo3 Received: 30 October 2023 / Revised: 30 January 2025 / Accepted: 17 February 2025 / Published online: 1 March 2025 © The Author(s) 2025 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 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 allows 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. Keywords 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. BJulia Novo julia.nov[email protected] Javier de Frutos [email protected]a.es Bosco García-Archilla [email protected] 1Instituto de Investigación en Matemáticas (IMUVA), Universidad de Valladolid, Valladolid, Spain 2Departamento de Matemática Aplicada II, Universidad de Sevilla, Seville, Spain 3Departamento de Matemáticas, Universidad Autónoma de Madrid, Madrid, Spain 123 19 Page 2 of 28 Journal of Scientific Computing (2025) 103 :19 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 [10], 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 [14], a domain decomposition technique is considered. In [26] semi-Lagrangian methods are studied. The authors in [22] 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 [19]. In [9] a tensor decomposition approach is presented. In [10] an approach based on low-rank tensor train decompositions is applied. Methods using sparse grids for HJB equations are presented in [6]. The solution of HJB equations on a tree structure was presented in [2]. The author of [23,24] discusses an approach to certain nonlinear HJB PDEs which is not subject to the curse of dimensionality. The approach utilizes the max-plus algebra. In [8] a data-driven approach based on the knowledge of the value function and its gradient on sample points is developed. The authors of [3]presentanew approach where the value function is computed using radial basis functions. Expanded literature on the control of partial differential equations using dynamic programming approach can be found in the last two references. 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 ⊂Rnand project then onto a reduced space r⊂Rrwith r<nto 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 r⊂Rr. In the second approach, the authors use a uniform mesh over the reduced space 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 [7], 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, his the time step while kis the spatial mesh diameter. This error bound improves existing results in the literature, where only O(k/h) error bounds are proved, see [12,13]. To bound the error in the first method in [1] the authors follow the technique in [12, Corollary 2.4], [13, Theorem 1.3] obtaining a bound for the error of size O(k/h).Forthe second method in [1], the factor 1/halso 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 [15] 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.The idea of using snapshots approaching the time derivatives is not new, although most of the references in the literature employ first order difference quotients (DQs) (i.e. first order divided finite 123 Journal of Scientific Computing (2025) 103 :19 Page 3 of 28 19 differences) instead of Galerkin time derivatives, as in [15] and the present paper. In [21] the set of snapshots (at different times) is increased with DQs to carry out the error analysis for the case in which projections respect to the H1 0norm are considered. In a more recent paper, [20], the authors show that the use of DQs has the added property of allowing to prove pointwise estimates in time. In a later paper, [11], the authors prove that one does not need to double the set of snapshots with values at different times plus DQs since only DQs plus a single initial value are enough to get pointwise estimates. This is a very interesting result because one can work with the same number of snapshots as in the standard case (the one with only values of the states at different times). In [15], the authors prove that this is also the case with time derivatives. A set of snapshots based on time derivatives plus the snapshot at the initial time (or the mean value of the states) is able to provide pointwise in time error estimates. This is the idea we apply in the present paper. Moreover, we carry out a different error analysis based on the recent results obtained in [7] that allow us to get sharper error bounds free of 1/hfactors. This is in agreement with the numerical investigations in the literature where the 1/hbehaviour in the error bounds of fully discrete methods has never been observed. Also, the use of snapshots based on time derivatives allows 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 [7] together with the use of snapshots based on time derivatives, are the key ingredients to get error bounds for the new method that are optimal in terms of the time step hand the mesh diameter of the reduced space kr.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 Sect.2we state the model problem and some preliminary results. In Sect.3we introduce the POD approximation and carry out the error analysis of the method. Finally, in Sect.4we show some numerical experiments in which we implement the method we propose in the paper. In the experiments of Sect.4we choose the same numerical tests as in [1] to compare our results with 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, ·denotes any norm associated to an inner product and ·∞denotes the maximum norm for vectors in Rn,n≥1. We will also denote by ·2the standard euclidean norm. In particular, in the numerical experiments, we use a weighted norm ·slightly different from the standard euclidean norm ·2. For a nonlinear mapping f:Rn×Rm→Rn, and a given initial condition y0∈Rnlet us consider the controlled nonlinear dynamical system ˙y(t)=f(y(t), u(t)) ∈Rn,t>0,y(0)=y0∈Rn,(1) together with the infinite horizon cost functional J(y,u)=∞ 0 g(y(t), u(t))e−λtdt.(2) 123 19 Page 4 of 28 Journal of Scientific Computing (2025) 103 :19 In (2)λ>0 is a given weighting parameter and g:Rn×Rm→R. The set of admissible controls is Uad ={u∈U|u(t)∈Uad for almost all t≥0}, where U=L2(0,∞;Rm)and Uad ⊂Rmis a compact convex subset. As in [1, Assumption 2.1] we assume the following hypotheses: •The right-hand side fin (1) is continuous and globally Lipschitz-continuous in both the first and second arguments; i.e., there exists a constant Lf>0 satisfying f(y,u)−f(˜y,u)≤Lfy−˜y,∀y,˜y∈Rn,u∈Uad,(3) f(y,u)−f(y,˜u)≤Lfu−˜u,∀u,˜u∈Uad,y∈Rn.(4) •The right-hand side fin (1) satisfies that there exists a constant Mf>0 such that the following bound holds f(y,u)∞≤Mf,∀y∈⊂Rn,u∈Uad,(5) where is a bounded polyhedron such that for sufficiently small h>0 the following inward pointing condition on the dynamics holds y+hf(y,u)∈, ∀y∈, u∈Uad.(6) •The running cost gis continuous and globally Lipschitz-continuous in both the first and second arguments; i.e., there exists a constant Lg>0 satisfying |g(y,u)−g(˜y,u)|≤Lgy−˜y,∀y,˜y∈Rn,u∈Uad,(7) |g(y,u)−g(y,˜u)|≤Lgu−˜u,∀u,˜u∈Uad,y∈Rn.(8) •Moreover, there exists a constant Mg>0 such that |g(y,u)|≤Mg,∀(y,u)∈×Uad.(9) From the assumptions made on fthere exists a unique solution of (1)y=y(y0,u)defined on [0,∞)for every admissible control u∈Uad and for every initial condition y0∈Rn,see [4, Chapter 3]. We define the reduced cost functional as follows: ˆ J(y0,u)=J(y(y0,u), u), ∀u∈Uad,y0∈Rn,(10) where y(y0,u)solves (1). Then, the optimal control can be formulated as follows: for given y0∈Rnwe consider min u∈Uad ˆ J(y0,u). The value function of the problem is defined as v:Rn→Ras follows: v(y)=inf ˆ J(y,u)|u∈Uad,y∈Rn.(11) This function gives the best value for every initial condition, given the set of admissible controls Uad. It is characterized as the viscosity solution of the HJB equation corresponding to the infinite horizon optimal control problem: λv(y)+sup u∈Uad {−f(y,u)·∇v(y)−g(y,u)}=0,y∈Rn.(12) 123 Journal of Scientific Computing (2025) 103 :19 Page 5 of 28 19 The solution of (12) is unique for sufficiently large λ,λ>max(Lg,Lf),[4]. Let us consider first a time discretization where his a strictly positive step size. We consider the following semidiscrete scheme for (12): vh(y)=min u∈Uad {(1−λh)vh(y+hf(y,u)) +hg(y,u)},y∈Rn.(13) As it is well-known equation (13) represents a numerical approximation related to the HJB equation (12) (see Remark 7). The following convergence result for the semidiscrete approximation [12, Theorem 2.3] requires that for (y,˜y,u)∈Rn×Rn×Uad f(y+˜y,u)−2f(y,u)+f(y−˜y,u)≤Cf˜y2,(14) g(y+˜y,u)−2g(y,u)+g(y−˜y,u)≤Cg˜y2.(15) Theorem 1 Let assumptions (3),(5),(6),(7),(9),(14)and (15)hold and let λ> max(2Lg,Lf).Letvand vhbe the solutions of (12)and (13), respectively. Then, there exists a constant C ≥0, that can be bounded explicitly, such that the following bound holds sup y∈Rn|v(y)−vh(y)|≤Ch,h∈[0,1/λ). (16) As in [1] let us suppose that there exists a bounded polyhedron ⊂Rnsuch that for h>0 small enough (6) holds. We consider a fully-discrete approximation to (12). Let Sjms j=1be a family of simplices which defines a regular triangulation of  = ms  j=1 Sj,k=max 1≤j≤ms (diam Sj). We assume we have nsvertices/nodes ˆy1,..., ˆynsin the triangulation. Let Vkbe the space of piecewise affine functions from to Rwhich are continuous in having constant gradients in the interior of any simplex Sjof the triangulation. Then, a fully discrete scheme for the HJB equations is given by vh,k(ˆyi)=min u∈Uad (1−λh)vh,k(ˆyi+hf(ˆyi,u)) +hg(ˆyi,u),(17) for any vertex ˆyi∈. There exists a unique solution of (17) in the space Vk,see[4, 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 Lu>0 such that u(t)−u(s)2≤Lu|t−s|,(18) then first order of convergence both in time and space is proved in [7, Theorem 6]. Theorem 2 Assume conditions (3)–(5),(7)–(9)and (18)hold. Assume λ>L with L= CnL f. Then, for 0≤h≤1/(2λ) there exist positive constants C1=C1(λ, Mf,Mg,Lf,Lg) and C2=C2(λ, Lf,Lg,Lu)such that |v(y)−vh,k(y)|≤C1(h+k)+C2h,y∈. Condition (18) can be weakened and one can still get convergence as proved in [7, Theorem 7]. Assume the following convexity assumption introduced in [5, (A4)] and denoted by (CA) as in [5,7], 123 19 Page 6 of 28 Journal of Scientific Computing (2025) 103 :19 •(CA) For every y∈Rn, {f(y,u), g(y,u), u∈Uad} is a convex subset of Rn+1. Theorem 3 Assume conditions (3),(4),(5),(7),(8),(9)and (CA) hold. Assume λ>L with L defined as in Theorem 2. Then, for 0≤h≤1/(2λ) there exist positive constants C1=C1(λ, Mf,Mg,Lf,Lg)and C2=C2(λ, Mf,Mg,Lf,Lg)such that for y ∈ |v(y)−vh,k(y)|≤C1(h+k)+C2 1 (1+β)2λ2(log(h))2h1 1+β,β=√nL f λ.(19) Let us observe that since β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 [7]and[15]. 3.1 POD Approximation Based on Time Derivatives For p∈Nlet us choose different pairs (uν,yν 0)p ν=1in U×.SinceU=L2(0,∞;Rm) the controls do not need to be constants, as those taken in the numerical experiments. By yν=y(uν;yν 0),ν=1,...,p, we denote the solutions of (1) corresponding to those chosen initial conditions and controls. Let us fix T>0andM>0andtaket=T/Mand tj=jt,j=0,...,M.For N=M+1 we define the following space V=span zν 1,zν 2,...,zν Np ν=1, with zν 1=√N yν,yν=1 N M  j=0 yν(tj) zν j=τyν t(tj−1), j=2,...,N, so that V=span √N yν,τyν t(t1),...,τyν t(tN)p ν=1, where the factor τin front of the temporal derivatives is a time scale and it makes the snapshots dimensionally correct. In the numerical experiments we take τ=1. The correlation matrix corresponding to the snapshots is given by K=((ki,j)) ∈RpN×pN, with the entries ki,j=1 pN zi k,zj l,k,l=1,...,N,i,j=1,...,p, 123 Journal of Scientific Computing (2025) 103 :19 Page 7 of 28 19 and where here, and in the sequel, (·,·)denotes the inner product in Rnto which the norm ||·||is associated. Let us denote for simplicity V=span w1,w 2,...,wpN:= z1 1,...z1 N,...,zp 1,...,zp N. Following [21], we denote by λ1≥λ2,... ≥λd>0 the positive eigenvalues of Kand by v1,...,vd∈RpN its associated eigenvectors of euclidean norm 1. Then, the (orthonormal) POD basis functions of Vare given by ϕk=1 √pN 1 √λk pN  j=1 vj kwj,k=1,...,d,(20) where vj kis the j-th component of the eigenvector vk. The following error estimate is known from [21, Proposition 1] 1 pN pN  j=1      wj− r  k=1 (wj,ϕ k)ϕk     2 = d  k=r+1 λk,(21) from which one can deduce for ν=1,...,p      yν− r  k=1 (yν,ϕ k)ϕk     2 +τ2 M+1 M  j=1      yν t(tj)− r  k=1 (yν t(tj), ϕk)ϕk     2 ≤p d  k=r+1 λk.(22) In the sequel, we will denote by Vr=span{ϕ1,ϕ 2,...,ϕ r},1≤r≤d,(23) and by PrRn→Vr, the orthogonal projection onto Vr.Then(21) can be written as 1 pN pN  j=1 wj−Prwj  2= d  k=r+1 λk. The following lemma is proved in [15, Lemma 3.2]. Lemma 1 Let T >0,t=T/M, tn=nt, n =0,1,...M, let X be a Banach space, z∈H2(0,T;X). Then, the following estimate holds max 0≤k≤Nzk2 X≤3z2 X+12T2 M M  n=1zn t2 X+16T 3(t)2T 0ztt(s)2 Xds,(24) where z=1 M+1M j=0zj. Using Lemma 1we can prove pointwise estimates for the projections onto Vr. Lemma 2 The following bounds hold for ν=1,...,p max 0≤j≤Myν(tj)−Pryν(tj)2≤3+24 T2 τ2p d  k=r+1 λk+16T 3(t)2T 0yν tt(s)2ds. (25) 123 19 Page 8 of 28 Journal of Scientific Computing (2025) 103 :19 Proof We argue as in [15, Lemma 3.4]. Taking z=yν(tj)−Pryν(tj)in (24) and applying (22) (taking into account that (M+1)/M≤2) yields max 0≤n≤Myν(tj)−Pryν(tj)2≤3+24 T2 τ2p d  k=r+1 λk +16T 3(t)2T 0yν tt(s)−Pryν tt(s)2ds. Now, since Pris an orthogonal projection, we have yν tt(s)−Pryν tt(s)2≤yν tt(s)2and the proof is finished  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 rwith r<n. To start we need to introduce some notation.We use a slightly different notation from the one used in [1]. In particular, as stated below, Pr c.isalwaysused to denote coefficients and ϕis always used for the linear combination based on the POD basis functions of the reduced order space. More precisely: for any y∈⊂Rnlet us denote by Pr cy∈Rrthe coefficients of the projection of yonto Vr Pr cy={(y,ϕ k)}r k=1.(26) For any yr∈Rrlet us denote by ϕyr∈Rnthe vector whose coefficients in the POD basis are the components of yr, i.e., ϕyr= r  j=1 yr jϕj,(27) where yr jis the jcomponent of the vector yr. For fand gin (1), (2)and(yr,u)∈Rr×Uad we define fr(yr,u)=Pr cf(ϕyr,u)∈Rr, gr(yr,u)=g(ϕyr,u)∈R.(28) 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 r⊂Rr satisfying Pr cy∈r,∀y∈. (29) The following lemma proves that the inward pointing condition for rfollows from (29). Lemma 3 Condition (29)implies that yr+hfr(yr,u)∈r,yr=Pr cy,y∈, provided the step size h or Pry−yis sufficiently small. Proof We follow [1, Remark 4.5] for the proof. We first observe that yr+hfr(yr,u)=Pr cy+hPr cf(ϕyr,u). Adding and subtracting hPr cf(y,u)we get yr+hfr(yr,u)=Pr c(y+hf(y,u)) +hPr c(f(ϕyr,u)−f(y,u)). (30) 123 Journal of Scientific Computing (2025) 103 :19 Page 9 of 28 19 Applying condition (6)y+hf(y,u)∈and applying (29) the first term on the right-hand side of (30) verifies Pr c(y+hf(y,u)) ∈r. Then, we only need to show that the second term on the right-hand side of (30) is small enough for hor Pry−ysufficiently small. Let us denote by z=f(ϕyr,u)−f(y,u)∈Rn.SincePr cz={(z,ϕ k)}r k=1∈Rr,taking into account that the functions ϕkdefine an orthonormal basis and that Pris a projection, we have Pr cz2=Prz≤z. Applying the above inequality together with (3), we get Pr cz2 2=Pr c(f(ϕyr,u)−f(y,u))2 2≤f(ϕyr,u)−f(y,u)2 ≤L2 fϕyr−y2=L2 fPry−y2,(31) so that the proof is concluded.  We can now define the reduced order problem we solve in practice. For frand grdefined in (28) and a given initial condition yr 0∈Rrlet us consider the controlled nonlinear dynamical system ˙yr(t)=fr(yr(t), u(t)) ∈Rr,t>0,yr(0)=yr 0∈Rr,(32) together with the infinite horizon cost functional Jr(yr,u)=∞ 0 gr(yr(t), u(t))e−λtdt.(33) As in (10), we define the reduced cost functional ˆ Jr(yr 0,u)=Jr(yr(yr 0,u), u), ∀u∈Uad,yr 0∈Rr,(34) where yr(yr 0,u)solves (32). Then, the POD optimal control can be formulated as follows: for given yr 0∈Rrwe consider min u∈Uad ˆ Jr(yr 0,u). The value function of the problem vr:Rr→Ris defined as follows: vr(yr)=inf ˆ J(yr,u)|u∈Uad,yr∈Rr.(35) Remark 1 It is easy to check that the regularity assumptions for frand granalogous to those for fand g,(3), (4), (5), (7), (8)and(9), hold from the definition of frand grand the properties being true for fand g. To get in practice a fully discrete approximation in the reduced space let us define Sr jmr s j=1 a family of simplices which defines a regular triangulation of r.Weassumewehavens vertices/nodes in the triangulation ˆy1 r,..., ˆyns r∈rand r= mr s  j=1 Sr j,kr=max 1≤j≤mr s (diam Sr j). Let Vkrbe the space of piecewise affine functions from rto Rwhich are continuous in rhaving constant gradients in the interior of any simplex Sr jof the triangulation. As in [1, 123 19 Page 16 of 28 Journal of Scientific Computing (2025) 103 :19 Section 1.2], the controls could not be unique but one can select the control with minimum norm. Now, let us observe that for the discrete value function it holds vh,k(ˆyi+hf(ˆyi,ui h,k)) =vh,k(ˆyi)+hf(yi,ui h,k)·∇vh,k(ˆyi)+O(h2). (52) Taking into account that vh,k(ˆyi)=(1−λh)vh,k(ˆyi+hf(ˆyi,ui h,k)) +hg(ˆyi,ui h,k), (53) and inserting (52)into(53)weget vh,k(ˆyi)=(1−λh)vh,k(ˆyi)+hf(yi,ui h,k)·∇vh,k(ˆyi)+O(h2) +hg(ˆyi,ui h,k). And then λhvh,k(ˆyi)=hf(ˆyi,ui h,k)·∇vh,k(ˆyi)+hg(ˆyi,ui h,k)+O(h2). From which λvh,k(ˆyi)=f(ˆyi,ui h,k)·∇vh,k(ˆyi)+g(ˆyi,ui h,k)+O(h). Now, since λv( ˆyi)=f(ˆyi,ui)·∇v(ˆyi)+g(ˆyi,ui), (54) and vh,k(ˆyi)→v(ˆyi), for h,k→0 we obtain f(ˆyi,ui h,k)·∇vh,k(ˆyi)+g(ˆyi,ui h,k)→f(ˆyi,ui)·∇v(ˆyi)+g(ˆyi,ui). (55) Arguing as in [13, Section 1.2], let us define L(y,u)=1 λ(f(y,u)·∇v(y)+g(y,u)), and let us associate with ya (unique) control u(y)such that L(y,u(y)) =min u∈Uad L(y,u)=v(y). Assume ∇vh,k(ˆyi)→∇v(ˆyi)(56) (which we have not proved) and ui h,k→uifor h,k→0. Then, on the one hand, from (55), L(ˆyi,ui h,k)→L(ˆyi,ui), and, on the other L(ˆyi,ui h,k)→L(ˆyi,ui) which implies ui=uiand ui h,k→uifor h,k→0. Finally, let us observe that the argument in [13, Section 1.2] had already proved that for any fixed hand k→0 the fully discrete controls ui h,kconverge to the corresponding semi-discrete time control defined in (13), for that value of h. 123 Journal of Scientific Computing (2025) 103 :19 Page 17 of 28 19 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]. The authors in [1] apply state snapshots in the reduced order method instead of snapshots based on time derivatives. We observe that, as explained in detail in the introduction, in the last case it is not necessary to consider both, state snapshots and time derivatives, since it has already been proved that only with time derivatives optimal bounds can be obtained. We observe that we have chosen the closed-loop control type approach instead of the open-loop control type approach in the present paper. Also, the theory of the present paper develops the first approach. We do not compare the present method with methods based on the first approach since our aim is just to propose, analyze and check in practice a new method that could be better or not (probably depending on the examples) than other methods in the literature. The numerical experiments of this section show that our method works fine in practice and is able to provide accurate approximations. We first notice that due to numerical reasons we have to choose a finite time horizon, so we select a sufficiently large te>0, which, in the experiments that follow, it was fixed to te=3. As in [1], we consider the following convection-reaction-diffusion equation zt−εzxx +γzx+μ(z3−z)=ub in I×(0,te), z(·,0)=z0in I, z(·,t)=0in∂I×(0,te), (57) with ε=1/10, and where I=(0,a)is an open interval, z:I×[0,te]→Rdenotes the state, and γand μare positive constants. The controls ubelong to the closed, convex, bounded set Uad =L2(0,te,[ua,ub]), for real values ua<ub. The cost functional to minimize is given by te 0 e−λtz(·,t,u)2 L2(I)+1 100 |u(t)|2dt,(58) wherewesetλ=1. Notice then that in (58) the aim is to drive the state to zero. We use a finite-difference method on a uniform grid of size x=l/Nwith N=100 on the interval I=(0,l)to discretize (57) 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) [25], 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) vr h,k(yi r)=min u∈Uad (1−λh)vr h,k(yi r+hfr(yi r,u)) +hgr(yi r,u),i=1,...,ns, was solved by fixed-point iteration, the stopping criterium being that two consecutive iterates differ in the maximum norm in less than a given tolerance TOLv, initially set to TOLv= 5×10−4. For the first iterate we choose a family of constant controls ul,l=1,...,pand at any point of the mesh yi r(initial condition) we compute the approximate solution of (32) corresponding to this initial condition and control ul. Then, we compute the value of the functional cost (33). Finally, the value of the initial iterate at yi ris the minimum between the values of the functional cost for l=1,...,p. 123 19 Page 18 of 28 Journal of Scientific Computing (2025) 103 :19 Once (36) is solved, the optimal control ur h,k(yi r)=argminu∈Uad (1−λh)vr h,k(yi r+hfr(yi r,u)) +hgr(yi r,u),(59) is obtained at any mesh point yi r,i=1,...,ns. Then, the suboptimal feedback operator r(y)is computed by interpolation. This means that for y∈we project onto the POD space to get Pryand then Pry= ns  i=1 μiyi r, r(y)= ns  i=1 μiur h,k(yi r), where the coefficients μisatisfy 0 ≤μi≤1, ns i=1μi=1. With this, the closed-loop system y(t)=f(y(t), r(y(t))), y(0)=y0,(60) 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. We will see below that very different approximations to the solution of (60) can be obtained with different values of the tolerance TOLvfor the fixed-point iteration solving (36)(seeFig.3), so that we solved this equation for decreasing values of TOLv, each one 5 times smaller than the previous one until the relative error between the solutions of (60) corresponding to two consecutive values of TOLvwas below 10% (it usually turned out to drop dramatically from above 10% to less than 0.01%). Here and in the sequel, by the relative error of a quantity ˆywith respect to y we mean y−ˆy/max(|y|,10−3). For the optimal HJB states, for every value of time tfor which the solution of (60) was computed, we computed the maximum or the relative errors of the components of y(t). For Test 2 in Sect. 4.2, due to the discontinous initial datum, it proved impossible to drive the relative error of two optimal HJB states computed with two different tolerances TOLvbelow 10%, so that we checked that the value of the relative errors measured in the norm (63) below was smaller than 10%. With respect to the computational cost of solving (59) by fixed-point iteration, it is obviously proportional to the number of iterations. In the experiments below, these values were 1382 for all values of rin Test 1 below, 362, 263 and 242 for r=2,3,4, respectively, in Test 2 below, and 365, 1103 and 1539 for r=3,4,5 in Test 3 in Sect. 4.3 below. On each iteration, the bulk of the cost is finding the nonnegative scalars μi j,j=1,...,ns, such that yi r+hfr(yi r,u)=μi 1y1 r+···+μi nsyns r, which was 70% of the cost of the iteration for r=5 in Test 3 in Sect.4.3, 79% for r=4 in Test 2, and 95% for r=4 in Test 1 in Sect. 4.1, followed by the cost of obtaining fr(yi r,u),i=1,...,ns, which was 4% for r=4inTests 1 and 2 below to 28% for r=5 in Test 3 in Sect.4.3. We note that cost of obtaining fr(yi r,u) can be substantially diminished using appropriate tensors or by means of techniques like discrete empirical inerpolation, which, for simplicity, we did not use in our codes. 4.1 Test 1: Semilinear Equation As in [1], we consider (57) with γ=0andμ=1, a=1andb(x)=z0(x)=2x(1−x). It is easy to check that the uncontrolled solution converges, as t→∞to a non-null steady state (see also [1, Fig. 6.1]), and that the null solution is unstable. 123 Journal of Scientific Computing (2025) 103 :19 Page 19 of 28 19 For the finite-difference approximation, we consider y:[0,te]→RN−1with components yj(t)≈z(xj,t),xj=jx,j=1,...,N−1, x=1/N, solution of Cyt=1 10 Ay +C(F(y)+uB)(61) where the components of Fand Bare, respectively Fj=yj(1−y2 j),Bj=2xj(1−xj), j=1,...,N−1, and Aand Care (N−1)×(N−1)tridiagonal matrices matrices given by A=1 (x)2 ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ −21 1−21 ......... 1−21 1−2 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ,C=1 12 ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 10 1 1101 ......... 1101 110 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ,(62) so that the finite-difference discretization (61) is fourth-order convergent. The norm we consider in RN−1is given by y2=x N−1  j=1 y2 j.(63) Let as observe that this norm is an approximation to the integral 1 0y(x)2dx of a function with values yjat the spatial mesh nodes. To compute the snapshots, as in [1], for constant controls u∈Usnap ={−1,0,1},we obtained the solutions y(n)=y(tn)of (61)every1/20 time units, that is, for tn=n/20, n=1,...20te, and then the time derivatives y(n) twere computed from identity (61). For the reduced spaces, we consider the cases of POD basis with only r=2,3 and 4 elements, also as in [1]. The POD approximation yrwas then the mean of the snapshots plus a linear combination of the POD basis. The control set Uad is given by 41 controls equally distributed in [−1,1]. As in [1], to define the domain r, we compute the projections of all the snapshots. With this procedure we obtain a set of points in Rr. Then, we define an hypercube containing this set of points. The aim of this procedure, in view of Lemma 3, is that the set rdefined in this way satisfies the invariance condition yr+hfr(yr,u)∈r,yr∈r,u∈Uad.(64) The set rfor r=4wasgivenby r=[−0.87,0.41]×[−0.01,0.02]×[−0.01,0.01]×[−0.01,0.01]. For this set we checked that condition (64) holds. We notice that our set ris considerable smaller than the corresponding set in [1](see [1, 6.1. Test 1]) where the authors use the standard euclidean norm in Rnrather than the norm (63) we use here. Since the domain is smaller we also consider partitions of rsmaller than those in [1]. We take maximum diameter kr=0.01, and, as in [1], we choose h=0.1kr. In Fig. 1we 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. 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 123 19 Page 20 of 28 Journal of Scientific Computing (2025) 103 :19 Fig. 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). The red crosses correspond to the values of the controls that we have joined by a blue line Fig. 2 Test 1: Value of the cost functional (58) on the optimal HJB states for r=2,3,4. The red crosses correspond to the values of the cost values that are joined by a pdf line Fig. 3 Test 1: Relative error between the optimal HJB states with r=4 corresponding to solving (36)by fixed point iteration with tolerances TOLv=5×10−4and TOLv=1×10−4(left), TOLv=1×10−4 and TOLv=2×10−5(centre), and optimal HJB controls (right) integrator that we use are more accurate than those in [1]. We also notice that there is little discrepancy between the values of the optimal HJB states for the different values of rthat we tried. We also computed the values of the cost functional (58) on the optimal HJB states for the three values of r. The values are shown in Fig. 2. It can be seen that the values decrease with rand that they differ in the ninth significant digit. As mentioned above there can be a significant difference between the optimal HJB states computed with solutions obtained by solving (36) with different tolerances TOLv.InFig.3we show the relative error between the optimal HJB states corresponding to tolerances TOLv= 5×10−4and TOLv=1×10−4(left), and between this one and that corresponding to TOLv= 2×10−5(centre). The right plot shows the corresponding optimal HJB controls. Figure 3 123 Journal of Scientific Computing (2025) 103 :19 Page 21 of 28 19 Table 1 Relative errors of the optimal HJB states and controls for r=3 computed with x=1/N,N=25 and N=50, with respect to those computed with N=100 Ny Rate r(y)rate 25 7.24 ×10−51.99 ×10−5 50 4.25 ×10−64.09 1.17 ×10−64.09 Fig. 4 Test 1: Results for different values of kr; relative errors between HJB states (top left) and controls (bottom left) with respect to to kr=0.005; HBJ controls (centre) and values of the cost functional (58) (right) shows the importance of solving (36) accurately in order to obtain good optimal HJB states, thus, justifying that we computed the (approximations to the) solution of (36) with decreasing values of TOLvuntil the relative error of the corresponding optimal HJB estates was below 10%. Following the suggestion of one of the reviewers, we checked if the optimal control computed with x=1/100 stabilizes the same PDE computed on finer meshes. We tried this for r=4andx=1/400 and 1/800, and the controlled solutions converged to zero as fast as in the case x=1/100. Given the accuracy with which we had computed the controls for x=1/100 and the accuracy of the discretization itself, we did not expect otherwise. The results above suggest that, for this problem, it is enough with r=3. For this value of rwe now check the effect of the finite-difference mesh in the optimal HJB states. In Table 1we show the relative errors of the optimal HJB states computed with N=25 and N=50 with respect to that computed with N=100, as well as the relative errors of the corresponding controls. For the controls, we show in Table 1the maximum for all values of t∈{0,0.05,0.1,...,3}of the relative errors, and for the states we show the maximum onthesamevaluesoftof the maximum of the relative errors of the state on all the values of the corresponding spatial grid. They confirm that the finite-difference discretization is of order 4. Due to the excellent accuracy obtained with x=1/50, the results that follow are done with that value of x. We now check the effect of different values krof the diameter of the partition of r.Todo this, we compare the results obtained with r=3, x=1/50 and kr=0.02,0.01,0.005. In order not to spoil the better accuracy obtained with the smaller value of krwe took Uad with 161 controls equally distributed in [−1,1]for kr=0.005, and, to simplify computations with only 11 controls for kr=0.02 (we also try with 21 and 41 controls, but, although we do not have at present an explanation for it, using only 11 controls with kr=0.02 gave somewhat better results). The results can be seen in Fig. 4. The plots on the left show the (maximum of the 51 points of the spatial grid of the) relative errors of the optimal HJB states (top) and their controls (bottom) for kr=0.02 and kr=0.01 with respect to those of kr=0.005, while the plot in the centre shows the HJB controls. The errors, as expected, are smaller for kr=0.01 than for kr=0.02, except for the controls for t∈[1.8,2.55]where they 123 19 Page 22 of 28 Journal of Scientific Computing (2025) 103 :19 Fig. 5 Test 1: Results for POD basis extracted from snapshots (x=1/50, r=3); Relative errors of HJB state (left) and control (centre) with respect to the case where POD basis is taken from time derivatives; HBJ control (right) are slightly larger. We also notice that, for kr=0.01, the relative errors remain below 10% except for t∈[1.5.2.2](where they remain below 18%) in the case of the errors in the optimal HJB states and t∈[1.8,3]for the controls. Notice, however, that the largest errors take place where both the states and the controls are close to zero (recall the plots in Fig. 1), were it is difficult to obtain small relative errors. Maybe this is the reason for the similar values of the cost functional (58) for the three values of kr, which are shown on the right plot in Fig. 4; the relative errors (with respect to kr=0.005) for kr=0.02 and kr=0.01 are 0.081% and 0.0023%, respectively. One may wonder what is the result if the snapshots are used to obtain the POD basis as in [1] instead of the time derivatives as in the present paper. Thus, we repeated our computations but replacing the time derivatives by the snapshots minus their mean. We did not find any significant difference. In Fig. 5we show the results corresponding to x=1/50, r=3, the set rbeing r=(−0.42,0.9)×(−0.01,0.02)×(−0.01,0.01). The optimal HJB control is shown on the right-plot, while the other two plots show the relative errors of the HJB state (left) and its control (centre) with respect to the results when the POD basis is taken from the time derivatives, r=3andx=1/100. We see that the relative errors are below 0.1%, and thus, no difference can be seen between the right-plot in Fig. 5and the centre plot in Fig. 1. We also notice that our results when the POD basis is taken from the snapshots are better than those in [1]. We believe that this is due to the higher accuracy of our computations (fourth-order convergent finite-difference method instead of a second-order convergent one, NDF with small tolerances to compute the snapshots instead of implicit Euler method, denser sets Uad for the control variable, smaller tolerance TOLvin the fixed point method to solve (36), etc). The fact that very similar results are obtained when the POD basis is taken from the snapshots or the time derivatives should not be surprising. As shown in [18], wether better results are obtained if the POD basis is extracted from the snapshots or from their difference quotients is case-dependent and, as shown in [16], very similar results are usually obtained when the POD basis is taken from the time derivatives or the snapshots difference quotients. The advantage of using time derivatives or difference quotients for the POD basis is more from the theoretical side, since it allows to prove optimal convergence of the POD methods with less assumptions than when the POD basis is extracted from the snapshots. More recently, in [17], it has been proved that using only snapshots for the POD basis, it is possible to prove error estimates for the corresponding POD methods with convergence rates as close to optimal as the smoothness of the solution from where the snapshots are taken allows. In 123 Journal of Scientific Computing (2025) 103 :19 Page 23 of 28 19 view of the recent results in [17], the analysis in the present paper can be easily adapted to cover also the case where POD basis is taken from the snapshots. 4.2 Test 2: Advection–Diffusion Equation As in [1], we now consider (57) with γ=1andμ=0, I=(0,2)and z0(x)= max(0,0.5sin(π x)).Wetakebas 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 dyj dt =yj+1−2yj+yj−1 10(x)2−yj+1−yj−1 2x+b(xj), j=1,...,N−1,y0=yN=0, which is second order convergent in problems with sufficiently smooth solutions. Since the initial state z0does not possess second-order derivatives in L2, we notice then that the time derivative ztblows up when t→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))/tof states at t=0andt=t. Perhaps also for lack bounded time derivatives at t=0 and the more dissipative nature of the implicit Euler method, we found that, in the computation of the optimal HBJ states and controls, better results were obtained if the implicit Euler method with t=1/20 was used instead of the NDF with small tolerances. Also, for reasons that we do not understand at present, we found that better results were obtained when the POD approximation was a linear combination of the POD basis plus the initial condition y0, instead of a linear combination of the POD basis plus the mean as in the previous section. In the previous test we had an invariance set so that we did not need to impose any boundary condition for solving (36). In this test we found it impossible to find a set r satisfying condition (64) both when the POD basis is taken from the states and from their time derivatives. In this last case the set rfor r=4 we used in the experiments was r=(−0.5,0.7)×(−0.3,1.5)×(−0.3,0.2)×(−0.05,0.15). To overcome the lack of invariance of this set, whenever for a vertex yi rwe had yi r+ hfr(yi r,u)/∈r, we simply replaced yi r+hfr(yi r,u)by its closest point on ∂r.This resulted in changing the value of yi r+hfr(yi 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 yi r+hfr(yi 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 Riccati equation. In Fig. 6we show the uncontrolled solution (left), the optimal LQR state (middle) and the optimal LQR control (right). In Fig.7we 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 Fig.7 compare very well with the optimal LQR control of Fig.6even for the case with only r=2 basis functions in our POD method. 123 19 Page 24 of 28 Journal of Scientific Computing (2025) 103 :19 Fig. 6 Test 2: Uncontrolled solution (left), the optimal LQR state (middle) and the optimal LQR control (right) Fig. 7 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 Fig. 8(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 uHJB −uLQR/max(10−3,uLQR)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=3andr=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 Fig.8and that restricting ourselves to the time interval in which the optimal control is sufficiently away from zero the errors for r=2are also below 35%. Again, the results here (which correspond to kr=0.1) compare favourably with those in the literature. As in the previouis test, we checked for r=4 if the optimal control computed with x=1/100 stabilizes the same PDE discretized with =1/400 and 1/800 with similar results as in the previous section. 4.3 Test 3: A Two-Dimensional Reaction–Diffusion Equation We extend (57) to two dimensions. In particular, we consider, 123 Journal of Scientific Computing (2025) 103 :19 Page 25 of 28 19 Fig. 8 Test 2: relative errors uHJB −uLQR/max(10−3,uLQR)of the optimal HJB controls with respect to the optimal LQR control (left) and difference between the optimal LQR state and the optimal HJB state computed with r=4 POD basis functions (right) Fig. 9 Test 3: Uncontrolled solution at t=0,1.5,3 zt−εz+(z3−z)=ub in ×(0,te), z(·,0)=z0in , z(·,t)=0in∂ ×(0,te), (65) with ε=1/10, =[0,1]×[0,1]and z:×[0,te]denotes the state. The control u belongs to Uad =L2(0,te,[ua,ub]), with ua=−1andub=1. The cost function is as (58) but with the state measured in L2() instead of L2(I),thatis te 0 e−λtz(·,t,u)2 L2() +1 100 |u(t)|2dt, with λ=1 as before. Similarly to Sect.4.1,wetakeb(x,y)=z0(x,y)=4x(1−x)y(1−y) and te=3. In Fig. 9we show the uncontrolled solution at the initial time, at t=te/2 and t=te. For the finite-difference approximation, we consider y:[0,te]→R(N−1)2with components yk(t)≈z(xk,t),where,fork=(j−1)(N−1)+i,i,j=1,...,N−1, xk=(xi,yj), and xi=ix,yj=jy,x=y=1/N, solution of ˆ Cy t=1 10 ˆ Ay +ˆ C(ˆ F(y)+uˆ B)(66) where the components of ˆ Fand ˆ Bare, respectively ˆ Fk=yk(1−y2 k),ˆ Bk=4xi(1−xi)yj(1− yj),k=(j−1)(N−1)+i,i,j=1,...,N−1, and ˆ Aand ˆ Care (N−1)2×(N−1)2 matrices given by ˆ A=I⊗A+A⊗I,ˆ C=I⊗C+C⊗I,whereIis the identity of order N−1, ⊗represents the Kronecker product of matrices, and Aand Care the matrices in (62) so that, as in Sect.4.1, the finite-difference discretization (66) is fourth-order convergent. The 123