Solving Differential Riccati Equations A Nonlinear
Solving Differential Riccati Equations A Nonlinear
Solving Differential Riccati Equations A Nonlinear
Abstract. Differential algebraic Riccati equations are at the heart of many applications in
control theory. They are time-depent, matrix-valued, and in particular nonlinear equations that
require special methods for their solution. Low-rank methods have been used heavily computing
a low-rank solution at every step of a time-discretization. We propose the use of an all-at-once
arXiv:1912.06944v1 [math.NA] 14 Dec 2019
space-time solution leading to a large nonlinear space-time problem for which we propose the use
of a Newton–Kleinman iteration. Approximating the space-time problem in low-rank form requires
fewer applicatons of the discretized differential operator and gives a low-rank approximation to the
overall solution.
Key words. Optimal Control, Low-rank methods, Riccati equations, Non-linear problems
(1.2) Ṗ + AT P + P A − P BB T P + C T C = 0, P (tf ) = M,
where M ∈ Rn×n is associated with a terminal penalty term at the final time horizon
tf , see the details below. From a computational point of view, equation (1.2) poses
many numerical challenges. Indeed, in this case we have to solve for the unknown
P (·) : [0, tf ] → Rn×n , i.e., a time-varying matrix with n2 entries in each step. If the
system (1.1) results from a spatial semi-discretization of a partial differential equation
(PDE), a computationally efficient method is crucial. A key ingredient in this context
is that of a low numerical rank. For the algebraic counterpart of (1.2), it is known,
see, e.g., [9, 30, 32, 40] that P ≈ LLT can be well approximated by low-rank factors
L ∈ Rn×k with k n if at least the control or the observation matrices B and C
correspond to finite-dimensional operators. In particular, this holds true for many
relevant PDEs of parabolic type (cf. [32]) and practically realisable controllers.
In the time-varying case, only a few theoretical results have been obtained. For
a recent discussion on low-rank solutions of differential Riccati equations, we refer to
[41]. Nevertheless, many numerical approaches exist and most of them rely on a time
∗
† Institute of Mathematics and Scientific Computing, University of Graz, 8010 Graz, Austria
(tobias.breiten@uni-graz.at )
‡ University of Bath, Department of Mathematical Sciences, BA2 7AY Bath, United Kingdom
(s.dolgov@bath.ac.uk)
§ Technische Universität Chemnitz, Department of Mathematics, Scientific Computing Group,
where the operators G and H depend on the chosen time stepping scheme with s
stages. In the particular case of the DRE, an implicit Euler scheme for example
results in a series of algebraic Riccati equations
>
1 1
τA − I Pi + Pi τ A − I − τ Pi BB > Pi + τ C > C + Pi+1 = 0,
2 2
with τ denoting the time step size. Note that this equation runs backward in time with
Pi+1 already known at step i. For such equations, efficient numerical methods that
are based on a low-rank Newton-Kleinman iteration (more details are given below)
[31, 16] are known to perform well.
A slightly different approach is given by splitting methods [21] that divide (1.2)
into two parts
The benefit of this decomposition is that the first equation is affine and the second
equation can be solved exactly, see [21, Lemma 3.2]. One can thus apply splitting
schemes such as, e.g., Lie splitting or Strang splitting. As for the time stepping
schemes, the key idea is to approximate the unknown P (·) for each time step ti via
low-rank factors.
Let us emphasize that the previous methods all rely on storing the information
for each time step separately. Hence, if each matrix P (ti ) ≈ Li LTi , with Li ∈ Rn×ki
Ptf
one still has to store n · i=1 k entries. While this is reasonable for problems with
few time steps, it might cause difficulties if long time horizons tf or fine time grids
have to be considered. The strategy we pursue in this manuscript is different in the
sense that we proceed by discretizing the DRE in space and time simultaneously.
Of course, we then have a nonlinear problem of vast dimensionality defined on the
space-time cylinder. We here propose the use of multivariate tensor decompositions
SOLVING DRES: A NONLINEAR SPACE-TIME METHOD 3
that allow us to easily approximate the solution of a linear problem defined over a
high-dimensional space. This means that we propose the use of an outer nonlinear
solver of Newton-type that at its core needs to solve a linearized problem. For this
we propose low-rank tensor methods that depend on the rank r of the solution over
the whole time-interval and assumes this rank to be (approximately) constant. As
our numerical examples confirm, this seems to be a reasonable assumption and could
lead to a significant reduction in assembly time and storage requirements.
Structure of the paper. In the next section we introduce the main idea of
low-rank methods for all-at-once approaches. These methods tackle the solution of
the matrix equation in a holistic way, where both space and time are discretized
simultaneously and then solved in a coupled fashion. We also motivate that the rank
of a space-time-solution is often found to be small so that our approach seems a feasible
alternative. The heart of this paper is given in Section 3, where we introduce the low-
rank method for solving the differential Riccati equation. The method of choice relies
on the tensor train format, which we introduce in Section 4. In the last section we
present results for various setups and illustrate that our suggested all-at-once scheme
provides a viable alternative to other methods.
Notation. We start by recalling several common notation from numerical-
(multi)linear algebra. Given two matrices V ∈ Rp×q and W ∈ Rn×m , the Kronecker
product is defined by
w11 V . . . w1m V
W ⊗ V := .. .. ..
. . . .
wm
and remind the reader of the following relation between these two operators
W T ⊗ V vec(Y ) = vec(V Y W ).
Since our approach for (1.2) crucially depends on the above ideas, we provide a
more detailed introduction into the topic on the example of open-loop control. For
this purpose, let us consider a control system of the form
E ẋ(t) = Ax(t) + Bu(t), x(0) = x0 ,
y(t) = Cx(t),
where E, A ∈ Rn×n , B ∈ Rn×m denotes a control operator and C ∈ Rp×n is an
observation operator. Given a desired trajectory xd : [0, tf ] → Rn , we want to solve
the optimal control problem
(2.1)
Ztf Ztf
1 2 β 1
min n J(x, u) : = ky(t) − Cxd (t)k dt + ku(t)k2 dt + ky(tf ) − Cxd (tf )k2
2
x∈L (0,tf ;R ) 2 2 2
u∈L2 (0,tf ;Rm ) 0 0
with k1,2,3 being small in comparison to nt . One can then implement the iterative
solver to maintain the low-rank style of the solution in combination with an additional
truncation scheme based on a truncated SVD or a QR reduction [42, 27, 8].
Example 2.1. Let us consider the one-dimensional heat equation
∂ ∂2
x(ξ, t) = 2 x(ξ, t) + χω u(ξ, t) in (0, 1) × (0, 2),
∂t ∂ξ
∂ ∂
x(0, t) = 0 = x(1, t) in (0, 2),
∂ξ ∂ξ
(ξ − 0.25)2
1
x(ξ, 0) = √ exp − in (0, 1),
0.05π 0.05
with control domain ω = (0.1, 0.4)∪(0.6, 0.9) and desired trajectory (see Fig. 1, middle)
(ξ − 0.75)2
1
xd (ξ, t) = (t − 2)x(ξ, 0) + t √ exp − .
0.02π 0.02
In Figure 1 (left), we show the actual state trajectory obtained by solving the discrete
optimality system (2.3) corresponding to a spatio-temporal discretization with n =
1000 and nt = 2000 grid points. Figure 1 (right) shows a rapid (exponential) decay
Singular values
0
10
X̄ = [x̄1 , . . . , x̄nt ]
P̄ = [p̄1 , . . . , p̄nt ]
10−8
σi
10−16
0 20 40
i
Fig. 1: 1D heat equation with n = 1000 and nt = 2000. Left. Controlled and desired
states. Right. Singular values of solutions.
of the singular values of the state and adjoint matrices, indicating that a low-rank
approximation of the form (2.4) can be computed for any desired accuracy.
We now may utilize the approximations (2.4) within our favourite and most suit-
able Krylov-subspace solver. In the case of linear PDE-constraints this method per-
forms very well. As the Riccati equation is nonlinear one needs to investigate how the
6 TOBIAS BREITEN AND SERGEY DOLGOV AND MARTIN STOLL
procedure changes for nonlinear equations. Our approach here is motivated by recent
work in [14], where one performs a nonlinear iteration as an outer solver, e.g., a New-
ton or Picard iteration, and then solves a linearized space-time discretized equation in
low-rank form. Of course the structure of the linearized systems is more complicated
than the one given above. Typically, the nonlinearity adds more terms to the sum
of the discretized matrices. The number of terms is depending on the rank of the
solution from the previous step of the nonlinear iteration. We derive the structure
of the linearized equations for our problem in Section 3.2. The goal of our solver
is to decouple the spatial and temporal degrees of freedom since every increase in n
and nt will result in a dramatic increase of the computational cost. To break the
curse-of-dimensionality we will use an outer non-linear solver and a space-time inner
solver that is based on the tensor train format introduced in Section 4.
3. A low-rank all-at-once method for the differential Riccati equation.
In this section, we extend the discussion from Section 2 on open-loop control to the
case of the closed-loop variant of problem (2.1). While the open-loop control problem
requires the solution of a linear system in saddle point form we here face the more
challenging differential Riccati equations. We begin with a summary of the theoretical
foundation of optimal closed-loop control via differential Riccati equations. For a
numerical example from PDE boundary control, we visualize the approximability of
the time-varying solution P (·) of (1.2) in terms of the singular values of corresponding
matrix unfoldings. Based on this observation, we derive a novel tensor-based Newton-
Kleinman method.
3.1. Optimal feedback control: low-rank or not low-rank. Let us again
consider an optimal control problem of the form
(3.1)
Ztf Ztf
1 2 β 1
min n J(x, u) : = ky(t) − Cxd (t)k dt + ku(t)k2 dt + ky(tf ) − Cxd (tf )k2
2
x∈L (0,tf ;R ) 2 2 2
u∈L2 (0,tf ;Rm ) 0 0
∂
x(ξ, t) = ∆ξ x(ξ, t) in Ω × (0, tf ),
∂t
m
∂ X
x(ξ, t) = 0 = αi (ξ)u(t) on Γ × (0, tf ),
∂ν i=1
x(ξ, 0) = x0 (ξ) in Ω,
Ztf Ztf
1 β
min 12 J(x, u) : = ky(t) − ek2R9 dt + ku(t)k2R12 dt.
2
u∈L (0,tf ;R ) 2 2
0 0
with e = (1, . . . , 1)T and yi (t) = |ω1i | ωi x(·, t) dt, i = 1, . . . , 9. We have implemented
R
a spatial discretization of the problem with n = 1089 piecewise linear finite elements
and a temporal implicit Euler scheme with nt = 1000 points. In Figure 2, we show
the singular values of the matrix unfoldings
P (t1 )
P1 = [vec(P (t1 )), . . . , vec(P (tnt ))] ∈ Rn ×nt , P2 = ... ∈ Rn·nt ×n .
2
P (tnt )
In both cases, the numerical ranks r(P1 ) = 252 and r(P2 ) = 54 (for machine preci-
sion) are significantly smaller than the maximum rank of r = 1000, allowing us to
approximate the tensor P ∈ Rn×n×nt by a low-rank representation P̃ . Some results
for such an approximation are given in Figure 3 where a full representation of the
solution P is compared to a tensor train approximation P̃ (cf. Section 4).
σi (P (·))
Γ4 Γ5 Γ6 σ1 (P (·)) P1 ∈ Rn ×nt
2
10−2
Γ3 Γ7 P2 ∈ Rn·nt ×n
10−5
Γ2 Γ8 10−8
ω7 ω8 ω9
ω4 ω5 ω6
10−11
ω1 ω2 ω3
Γ1 Γ9
Γ12 Γ11 Γ10 10−14
50 100 150 200 250
Fig. 2: Left. Control and observation domains. Right. Singular values of matrix
unfoldings.
8 TOBIAS BREITEN AND SERGEY DOLGOV AND MARTIN STOLL
kP −P̃ k #kB(P )
kP k r(P1 ) r(P2 ) #kB(P̃ ) #kB(P̃ )
ε 252 54 1.21 · 105 78
10 −14
250 52 1.16 · 105 82
10−12 216 37 7.18 · 104 132
10−10 176 29 4.62 · 104 205
10−8 133 22 2.70 · 104 354
10−6 85 15 1.20 · 104 793
10−4 45 8 3.59 · 103 2642
10−2 19 3 6.87 · 102 13810
for Pi+1 , where A(Pi ) = A−BB T Pi E denotes the closed-loop system operator associ-
ated with the current iteration. Note that each step still requires solving a differential
matrix equation for the time-varying unknown Pi+1 : [0, t] → Rn×n . However, in con-
trast to (DRE) the resulting differential Lyapunov equation is linear. Applying the
vectorization operator to both sides of the above equation, we obtain the equivalent
equation
E = ET ⊗ ET , L = E T ⊗ AT + AT ⊗ E T ,
M(Pi ) = E T ⊗ E T Pi BB T + E T Pi BB T ⊗ E T
C = vec(C T C), G(Pi ) = (E T Pi ⊗ E T Pi )B, B = vec(BB T ).
Let us emphasize that the subscript i indicates the iteration of the outer non-linear
solver and the superscript j enumerates the time-step of the current vector. Note that
the above scheme indeed is implicit since the equation is running backwards in time.
SOLVING DRES: A NONLINEAR SPACE-TIME METHOD 9
Similar as before, we represent all time steps in a single vector that we denote by
1
Pi
.. 2
pi := . ∈ Rn nt .
Pint
where
G(Pi1 )
−1 1 C
.. .. .. ..
1 . . .
.
D= , c = − , g(pi ) = −
. .
τ −1 1 C ..
−1 C + τ1 C G(Pint )
For the main idea, we may think of pi as being a tensor which we approximate. For
this we realize that we can write
Pij = eTj ⊗ In ⊗ In pi .
Pij ≈ eTj ⊗ In ⊗ In (u ⊗ v ⊗ w) = uj (v ⊗ w) .
For obtaining a matrix valued expression, let us interpret Pij = vec(Pij ) as the vec-
torization of a matrix. It then follows that
We use the latter expression for a derivation of the term M(Pij ) which is of the form
≈ uj E T ⊗ E T wv T BB T + E T wv T BB T ⊗ E T .
With this in mind, let us now move beyond the case of a simple rank-one approxima-
tion and utilize the tensor train format illustrated in Section 4. This format allows
the approximation of the state tensor via
rX
1 ,r2
with r1 and r2 the ranks of the tensor train approximation. As a result, we can
approximate the operator M via
rX
1 ,r2 T T
(1) T (1)
M(pi ) ≈ blkdiag (E T ⊗ E T us1 ,k u(3)
s2 u(2)
s1 ,s2 BB T
+ E u u(3)
s1 ,k s2 u(2)
s1 ,s2 BB T ⊗ E T )
s1 ,s2 =1 k=1,...,nt
rX1 ,r2 T T
= diag(u(1)
s1 ) ⊗(E
T
⊗ E T u(3)
s2 u(2)
s1 ,s2 BB T + E T u(3)
s2 u(2)
s1 ,s2 BB T ⊗ E T )
s1 ,s2 =1
| {z }
=:Uj
Note that we have included only one sum in the derivation of g(pi ) and have collected
T
(3) (2)
the sum over s2 into the matrix products with u1:r2 uj,1:r2 , where 1 : r2 indicates
r2 vectors collected into a n × r2 matrix, which is standard Matlab notation. Alto-
gether, we obtain a linear system of tensor product structure of the form Api+1 = f ,
where
(3.2) A = Int ⊗ E T ⊗ AT + Int ⊗ AT ⊗ E T + D ⊗ E T ⊗ E T + M(pi ),
rf
X
(3.3) f =− eτ ⊗ cTk ⊗ cTk + g(pi )
k=1
where eτ = [1, . . . , 1, 1 + τ1 ]T ∈ Rnt . The resulting system is now solved with the
Alternating Minimal Energy (AMEn) solver, presented next. Note that the number of
terms in M(pi ) depends on the previous solution and we will illustrate in the numerical
experiments how this behaves with respect to changes in the system parameters and
system size. We will use a tolerance εnewton to assess the convergence of the Newton
iterates via kpkp
i+1 −pi k
i+1 k
< εnewton .
3.3. Nested approach. As the nonlinear iteration can become quite expensive
with a possibly high rank of the state variables we aim to drastically reduce this cost
by using a nested approach [22, 37]. There the authors solve a nonlinear problem
on a coarse mesh and then transfer the solution to the next finer mesh as the initial
guess for the nonlinear iteration. This process is continued until the final mesh-size
is solved. Typically, this reduces the number of nonlinear iterations drastically. Our
aim is to apply this approach here as well. We will solve the low-rank DRE on a
coarse mesh to obtain a solution that we then transfer to the next finer mesh as the
intial guess for the Newton–Kleinman iteration.
SOLVING DRES: A NONLINEAR SPACE-TIME METHOD 11
we then obtain
rX
1 ,r2
This will now be a low-rank solution on the fine grid and it will be the starting guess
for the Newton–Kleinman iteration on the finer mesh. It is obvious that this process
can be continued further.
4. The tensor train decomposition and algorithms. The efficient solution
of tensor-valued equations has recently seen much progress [43, 3, 33, 19, 26, 20].
While there is a variety of available tensor formats we decide to use the so-called tensor
train (TT) representation [34, 33, 36]. As our problem essentially consists of tensors
of order three, the Tucker format [19, 26] would also be appropriate. The availability
of tailored solvers along with the possibility of designing and using preconditioners
makes the TT-format and its toolbox1 [35] an ideal candidate for our purposes.
In the following we will briefly introduce the TT format, while we point to the
2
literature for details. To this end, suppose p ∈ Rn nt is the approximate solution
vector of a space-time matrix equation as discussed in Section 3. However, its ele-
ments can be also naturally enumerated by three indices i1 , i2 , i3 , corresponding to
the discretization in time and the two spatial dimensions, respectively. Introducing a
multi-index
i1 i2 i3 = (i1 − 1)n2 + (i2 − 1)n + i3 ,
nt ,n,n
we can denote p = p(i1 i2 i3 ) i1 ,i2 ,i3 =1 , and consider p as a three-dimensional tensor
The ranges r1 , r2 are called TT ranks and the u(m) , m = 1, 2, 3 are the so-called
TT blocks, with u(1) ∈ Rnt ×r1 , u(2) ∈ Rr1 ×n×r2 and u(3) ∈ Rr2 ×n . When fixing
(2)
the indices we get for u(2) (i2 ) ∈ Rr1 ×r2 a matrix slice, for us1 ,s2 ∈ Rn a vector,
(2)
and for us1 ,s2 (i2 ) a scalar. The values of r1 , r2 depend on the accuracy enforced in
1 https://github.com/oseledets/TT-Toolbox
12 TOBIAS BREITEN AND SERGEY DOLGOV AND MARTIN STOLL
Note that the Kronecker product matrix assembly (3.2) is a particular case of (4.2).
For the numerical solution of Ap = f we use the Alternating Minimal Energy
(AMEn) algorithm [13], which is an enhanced version of the Alternating Linear
Scheme (ALS) [24]. The initial idea is to rewrite the TT decomposition (4.1) as
a linear map from the elements of a single TT block, and consider Ap = f as an
overdetermined system on these elements. Specifically for the three-dimensional TT
format, we introduce frame matrices
r2 T
X (2) 2
nt ×nt r1
(4.3) U6=1 = Int ⊗ u1:r1 ,s2 ⊗ u(3)
s2 ∈ Rn ,
s2 =1
T 2
nt ×r1 nr2
(4.4) U6=2 = u(1) ⊗ In ⊗ u(3) ∈ Rn ,
r1
X (2) 2
nt ×r2 n
(4.5) U6=3 = u(1)
s1 ⊗ us1 ,1:r2 ⊗ In ∈ Rn ,
s1 =1
(u(2) (n))T
One can notice that, given (4.1), the identity p = U6=k vecT (u(k) ) holds for any k =
1, 2, 3. Now we can plug this reduced basis ansatz into the original system Ap = f ,
and resolve it by projecting onto the same frame matrix,
(4.6) T
(U6= k AU6=k )vecT (u
(k) T
) = U6= kf .
Now, the ALS algorithm iterates over k = 1, 2, 3, solving (4.6) in each step, and
updating the TT block u(k) . This simple algorithm requires an initial guess in the
TT format (4.1) with fixed TT ranks, which might be difficult to guess a priori. To
circumvent this issue, AMEn algorithm computes additionally a TT approximation
of the current residual f − Ap in a form similar to (4.1) with (smallish) TT ranks
ρ1 , ρ2 , and expands the TT blocks of the solution with the TT blocks of this residual
approximation [13]. This allows us to increase the TT ranks of the solution from
r1 , r2 to r1 + ρ1 , r2 + ρ2 in each iteration and improve convergence. The iteration
continues until the accuracy and the TT ranks of the solution reach their desired
SOLVING DRES: A NONLINEAR SPACE-TIME METHOD 13
values. At the same time, we can descrease the ranks using SVD if they become too
large. Specifically, let u(1) , u(2) , u(3) be the TT blocks at the given AMEn iteration,
and let û(1) , û(2) , û(3) be the TT blocks after one full sweep over k = 1, 2, 3. We
choose a stopping tolerance εamen > 0 and stop the algorithm when
The same threshold εamen is used for SVD within the AMEn iteration when the TT
ranks need to be reduced.
For the practical efficiency we can notice that the reduced matrix Bk = U6=
T
k AU6=k
in (4.6) can be written in a Kronecker product form that inherits the original matrix
TT decomposition (4.2). In particular, we can write
RX
1 ,R2
(1) (2) (3)
Bk = B`1 ⊗ B`1 ,`2 ⊗ B`2 ,
`1 ,`2 =1
where B (k) = A(k) (that corresponds to the identity factor in the frame matrix U6=k ),
and hence is large but sparse, while the other matrices are dense but small, of sizes
r1 × r1 , r2 × r2 or even 1 × 1 if k = 1 or k = 3 is considered. This allows us to design
an efficient block Jacobi preconditioner for a GMRES solver of (4.6). Without loss
(2)
of generality, we can assume that k = 2, in which case B`1 ,`2 ∈ Rn×n is sparse, and
(1) (3)
B`1 ∈ Rr1 ×r1 and B`2 ∈ Rr2 ×r2 are dense. Then we construct the block diagonal
preconditioner as follows,
RX
1 ,R2
(1) (2) (3)
(4.7) P2 = diag(B`1 ) ⊗ B`1 ,`2 ⊗ diag(B`2 ).
`1 ,`2 =1
2 https://github.com/oseledets/TT-Toolbox
14 TOBIAS BREITEN AND SERGEY DOLGOV AND MARTIN STOLL
As illustrated in Figure 4 the control consists of five spatially constant control patches
ω1 , . . . , ω5 , given by
ω1 = {ξ ∈ Ω | kξ − ξM1 k2 ≤ 0.01},
ωi = {ξ ∈ Ω | kξ − ξMi k2 ≤ 0.0025}, i = 2, 3, 4, 5,
with centres defined
! as
1
2 0.1 0.78 0.2 0.9
ξM1 = 1 , ξM2 = , ξM3 = , ξM4 = , ξM5 = .
0.2 0.23 0.7 0.87
2
ω5
ω̃3 ω̃6 ω̃9
ω4
ω3
ω2
ω̃1 ω̃4 ω̃7
80
10−1
r1 (pi ), r2 (pi )
n = 289
kpi+1 −pi k
60 n = 1089
kpi+1 k
n = 289 n = 4225
−4
10 n = 1089 n = 16641
n = 4225 40
n = 16441
20
10−7
2 4 6 8 10 5 10
Iteration number i Iteration number i
in the previous case. However, the penalization parameter of the optimal control
problem is set to β = 1. As for the case of an interior control, we consider tracking
of a temporally constant state xd that is visualized in Figure 8. The results shown in
Figure 7 indicate a similar behaviour to the the previously computed case. Namely,
that eventually the quadratic convergence of the Newton method kicks in and that
again relying on the nester approach saves many nonlinear iterations and thus more
solves of the linearized space-time problem. Moreover, we do not observe an increase
of the TT-ranks during the Newton-Kleinman iteration.
80
10−1
r1 (pi ), r2 (pi )
n = 289
kpi+1 −pi k
60 n = 1089
kpi+1 k
n = 4225
−7
20
10
2 4 6 8 10 2 4 6 8 10
Iteration number i Iteration number i
80
10−1 n = 289
r1 (pi ), r2 (pi )
n = 1089
n = 289
kpi+1 −pi k
n = 4225
60 n = 1089
kpi+1 k
n = 16441
n = 4225
10−4 n = 16641
40
20
10−7
2 4 6 8 2 4 6 8
Iteration number i Iteration number i
3 https://www.mpi-magdeburg.mpg.de/projects/mess
18 TOBIAS BREITEN AND SERGEY DOLGOV AND MARTIN STOLL
1.2 1.2
1 1
1 1
0.8 0.8
0.4 0.4
0 0
1 1
0.2 0.2
1 1
0.5 0.5
0.5 0.5
0 0
0 0 0 0
Fig. 10: Left. Desired state xd (ξ1 , ξ2 ). Right. Controlled state x(ξ1 , ξ2 , tf ).
REFERENCES
[1] H. Abou-Kandil, G. Freiling, V. Ionescu, and G. Jank, Matrix Riccati equations, Sys-
tems & Control: Foundations & Applications, Birkhäuser Verlag, Basel, 2003. In control
and systems theory.
[2] N. Ahmed, G. Matthies, L. Tobiska, and H. Xie, Discontinuous Galerkin time step-
ping with local projection stabilization for transient convection–diffusion-reaction problems,
Computer Methods in Applied Mechanics and Engineering, 200 (2011), pp. 1747–1756.
[3] J. Ballani and L. Grasedyck, A projection method to solve linear systems in tensor format,
Numerical Linear Algebra with Applications, 20 (2013), pp. 27–43.
[4] R. Bellman, Dynamic programming, Courier Corporation, 2013.
[5] P. Benner, S. Dolgov, A. Onwunta, and M. Stoll, Low-rank solvers for unsteady Stokes-
Brinkman optimal control problem with random data, Computer Methods in Applied Me-
chanics and Engineering, 304 (2016), pp. 26–54.
[6] P. Benner and H. Mena, BDF methods for large-scale differential riccati equations, in in
Proceedings of Mathematical Theory of Network and Systems, MTNS, 2004.
[7] , Rosenbrock methods for solving Riccati differential equations, IEEE Transactions on
Automatic Control, 58 (2013), pp. 2950–2956.
[8] P. Benner, A. Onwunta, and M. Stoll, Low rank solution of unsteady diffusion equations
with stochastic coefficients, SIAM Journal on Uncertainty Quantification, 3 (2015), pp. 622–
649.
[9] P. Benner and J. Saak, Numerical solution of large and sparse continuous time algebraic
matrix Riccati and Lyapunov equations: a state of the art survey, GAMM-Mitteilungen,
36 (2013), pp. 32–52.
[10] A. Brooks and T. Hughes, Streamline upwind/Petrov-Galerkin formulations for convection
dominated flows with particular emphasis on the incompressible Navier-Stokes equations,
Computer Methods in Applied Mechanics and Engineering, 32 (1982), pp. 199–259.
SOLVING DRES: A NONLINEAR SPACE-TIME METHOD 19
[11] S. Dolgov, TT-GMRES: solution to a linear system in the structured tensor format, Russian
Journal of Numerical Analysis and Mathematical Modelling, 28 (2013), pp. 149–172.
[12] S. Dolgov and B. Khoromskij, Two-level QTT-Tucker format for optimized tensor calculus,
SIAM Journal on Matrix Analysis and Applications, 34 (2013), pp. 593–623.
[13] S. Dolgov and D. Savostyanov, Alternating minimal energy methods for linear systems in
higher dimensions, SIAM Journal on Scientific Computing, 36 (2014), pp. A2248–A2271.
[14] S. Dolgov and M. Stoll, Low-rank solutions to an optimization problem constrained by
the Navier-Stokes equations, SIAM Journal on Scientific Computing, 39 (2017), pp. A255–
A280.
[15] D. Donoho et al., High-dimensional data analysis: The curses and blessings of dimension-
ality, AMS math challenges lecture, 1 (2000), p. 32.
[16] F. Feitzinger, T. Hylla, and E. Sachs, Inexact Kleinman–Newton method for Riccati
equations, SIAM Journal on Matrix Analysis and Applications, 31 (2009), pp. 272–288.
[17] C. Fu and J. Pfaendtner, Lifting the curse of dimensionality on enhanced sampling of
reaction networks with parallel bias metadynamics, Journal of Chemical Theory and Com-
putation, 14 (2018), pp. 2516–2525.
[18] L. Grasedyck and W. Hackbusch, A multigrid method to solve large scale sylvester equa-
tions, SIAM Journal on Matrix Analysis and Applications, 29 (2007), pp. 870–894.
[19] L. Grasedyck, D. Kressner, and C. Tobler, A literature survey of low-rank tensor ap-
proximation techniques, GAMM-Mitteilungen, 36 (2013), pp. 53–78.
[20] W. Hackbusch, Tensor Spaces And Numerical Tensor Calculus, Springer–Verlag, Berlin,
2012.
[21] E. Hansen and T. Stillfjord, Convergence analysis for splitting of the abstract differential
Riccati equation, SIAM Journal on Numerical Analysis, 52 (2014), pp. 3128–3139.
[22] R. Herzog and E. Sachs, Preconditioned conjugate gradient method for optimal control
problems with control and state constraints, SIAM Journal on Matrix Analysis and Appli-
cations, 31 (2010), pp. 2291–2317.
[23] M. Hinze, Optimal and instantaneous control of the instationary Navier-Stokes equations,
habilitation, Technisches Universität Dresden, 2002.
[24] S. Holtz, T. Rohwedder, and R. Schneider, The alternating linear scheme for tensor
optimization in the tensor train format, SIAM Journal on Scientific Computing, 34 (2012),
pp. A683–A713.
[25] B. N. Khoromskij and V. Khoromskaia, Multigrid accelerated tensor approximation of
function related multidimensional arrays, SIAM Journal on Scientific Computing, 31
(2009), pp. 3002–3026.
[26] T. Kolda and B. Bader, Tensor decompositions and applications, SIAM Review, 51 (2009),
pp. 455–500.
[27] D. Kressner and C. Tobler, Krylov subspace methods for linear systems with tensor product
structure, SIAM Journal on Matrix Analysis and Applications, 31 (2010), pp. 1688–1714.
[28] N. Lang, Numerical methods for large-scale linear time-varying control systems and related
differential matrix equations, PhD thesis, Technische Universität Chemnitz, 2018.
[29] N. Lang, H. Mena, and J. Saak, On the benefits of the LDLT factorization for large-
scale differential matrix equation solvers, Linear Algebra and its Applications, 480 (2015),
pp. 44–71.
[30] Y. Lin and V. Simoncini, A new subspace iteration method for the algebraic Riccati equation,
Numerical Linear Algebra with Applications, 22 (2015), pp. 26–47.
[31] H. Mena, A. Ostermann, L.-M. Pfurtscheller, and C. Piazzola, Numerical low-rank
approximation of matrix differential equations, Journal of Computational and Applied
Mathematics, 340 (2018), pp. 602–614.
[32] M. R. Opmeer, Decay of singular values of the Gramians of infinite-dimensional systems, in
2015 European Control Conference (ECC), 2015, pp. 1183–1188.
[33] I. Oseledets, DMRG approach to fast linear algebra in the TT-format, Computational Meth-
ods in Applied Mathematics, 11 (2011), pp. 382–393.
[34] , Tensor-train decomposition, SIAM Journal on Scientific Computing, 33 (2011),
pp. 2295–2317.
[35] I. Oseledets, S. Dolgov, V. Kazeev, D. Savostyanov, O. Lebedeva, P. Zhlobich,
T. Mach, and L. Song, TT-Toolbox. https://github.com/oseledets/TT-Toolbox.
[36] I. Oseledets and E. Tyrtyshnikov, Breaking the curse of dimensionality, or how to use
SVD in many dimensions, SIAM Journal on Scientific Computing, 31 (2009), pp. 3744–
3759.
[37] J. W. Pearson, M. Stoll, and A. J. Wathen, Preconditioners for state-constrained optimal
control problems with Moreau-Yosida penalty function, Numerical Linear Algebra with
20 TOBIAS BREITEN AND SERGEY DOLGOV AND MARTIN STOLL