State-dependent Riccati equation feedback stabilization for nonlinear PDEs

The synthesis of suboptimal feedback laws for controlling nonlinear dynamics arising from semi-discretized PDEs is studied. An approach based on the State-dependent Riccati Equation (SDRE) is presented for 2 and ∞ control problems. Depending on the nonlinearity and the dimension of the resulting problem, offline, online, and hybrid offline-online alternatives to the SDRE synthesis are proposed. The hybrid offline-online SDRE method reduces to the sequential solution of Lyapunov equations, effectively enabling the computation of suboptimal feedback controls for two-dimensional PDEs. Numerical tests for the Sine-Gordon, degenerate Zeldovich, and viscous Burgers’ PDEs are presented, providing a thorough experimental assessment of the proposed methodology.


Introduction
Feedback control laws are ubiquitous in modern science and engineering and can be found in autonomous vehicles, fluid flow control, and network dynamics, among many others.A distinctive feature of feedback laws is their ability to compensate external perturbations in real time.While an offline training phase is often affordable, an operational feedback law must be able to provide a control signal at a rate that is determined by the underlying time scales of the physical system to be controlled.This requirement poses a formidable computational constraint for the synthesis of feedback controls which require an online optimization procedure.
A natural approach to generate an optimal feedback law for real-time control is to follow a dynamic programming approach.Here, we solve a nonlinear Hamilton-Jacobi-Bellman (HJB) partial differential equation (PDE) for the value function of the optimal control problem under study.This is done in an offline phase, and once the value function has been computed, a feedback law is obtained as by-product.Once online, the cost of implementing an HJBbased control in real time, assuming the current state of the system is available, is reduced to the evaluation of a nonlinear feedback map.Unfortunately, the dynamic programming approach is not suitable for systems described by large-scale dynamics, as the computational complexity of approximating the associated high-dimensional HJB PDEs goes beyond the reach of traditional computational methods.Only very recently, the use of effective computational approaches such as sparse grids [19,31], tree structure algorithms [2], polynomial approximation [28,29,4] tensor decomposition methods [47,21,18,39], and representation formulas [13,12] have addressed the solution of highdimensional HJB PDEs.Recent works making use of deep learning [24,17,38,26,34,37,30] anticipate that the synthesis of optimal feedback laws for largescale dynamics can be a viable path in the near future.
An alternative to the dynamic programming approach is the use of a Nonlinear Model Predictive Control (NMPC) framework [22].Here, an optimal open-loop control law is computed at every sampling instant of the dynamics.However, the control action is optimized over a prediction horizon, which is sufficiently large to guarantee asymptotic stability of the closed-loop, but short enough to ensure a computing time that is compatible with the rate at which the control law is called.The NMPC framework embodies the trade-off between optimality and real-time computability.It has been shown [23] that the NMPC corresponds to a relaxation of the dynamic programming approach, in the sense that the NMPC law is suboptimal with respect to the optimal stabilizing feedback law provided by the dynamic programming approach, although its suboptimality can be controlled by increasing the prediction horizon.
There exists a third synthesis alternative, which incorporates elements from both dynamic programming and NMPC, known as the State-Dependent Riccati Equation (SDRE) approach [15,16].The SDRE method originates from the dynamic programming and the HJB PDE associated to infinite horizon optimal stabilization, however, it circumvents its solution by reformulating the feedback synthesis as the sequential solution of state-dependent Algebraic Riccati Equations (ARE), which are updated online along a trajectory.The SDRE feedback is implemented similarly as in NMPC, but the online solution of an optimization problem is replaced by a nonlinear matrix equation.Alternatives to the online formulation include the use of neural networks [48,1] in supervised learning, and reformulating the SDRE synthesis as an optimization problem [27].However, in this work we will focus on addressing the SDRE synthesis from a numerical linear algebra perspective.The efficient solution of matrix equations arising in feedback control has been subject of extensive research over the last decades, leading to the design of solvers which can effectively be applied to large-scale dynamics such as those arising in the control of systems governed by PDEs (see, e.g., [3,8,9,10,44,46]), and agent-based models [25].Moreover, under certain stabilizability conditions, this feedback law generates a locally asymptotically stable closed-loop and approximates the optimal feedback law from the HJB PDE.
The purpose of the present work is to study the design of SDRE feedback laws for the control of nonlinear PDEs, including nonlinear reaction and transport terms.This is a class of problems that are inherently large-scale, and where the presence of nonlinearities renders linear controllers underperformant.In this framework, the use of dynamic programming or NMPC methods is often prohibitively expensive, as for instance, in feedback control for multidimensional PDEs.Here instead, we propose and assess different alternatives for control design based on the SDRE approach which can be effectively implemented for two and three dimensional PDEs.The methodology studied in this paper is based on the parametrization of the SDRE synthesis proposed in [6,5].In these works, different SDRE feedback laws are developed based on the representation of the nonlinear terms in the system dynamics.This is particularly relevant to PDE dynamics, as unlike the nonlinear ODE world, there exists a clear taxonomy of physically meaningful nonlinearities.In some particular cases, the SDRE approach is reduced to a series of offline computations, and a real-time controller only requires a nonlinear feedback evaluation.We explore the limitations of such a parametrization.In the more general case, the SDRE synthesis requires the sequential solution of AREs at a high rate, a task that can is demanding for large-scale dynamics.Here, we propose a variant requiring the solution of a Lyapunov equation instead, thus mitigating the computational cost associated to the online synthesis.
The paper is structured as follows.In section 2 and its subsections we review the H 2 /H ∞ optimal feedback synthesis problem.In section 3 and its subsections we present the suboptimal approximation to these feedback laws by means of the SDRE approach, including offline, online, and hybrid offlineonline implementations, illustrated with an application to the control of the Sine-Gordon equation.Section 4 discusses numerical linear algebra aspects of the solution of the Algebraic Riccati and Lyapunov equations which constitute the core building blocks of the SDRE feedback synthesis.Finally, in section 5 we perform a computational assessment of the proposed methodology on the two-dimensional degenerate Zeldovich and viscous Burgers PDEs.

Optimal feedback synthesis for nonlinear dynamics
In this section we revisit the use of dynamic programming and Hamilton-Jacobi-Bellman/Isaacs PDEs for the computation of optimal feedback controls for nonlinear dynamics.We begin by stating the problem of optimal feedback stabilization via H 2 synthesis, to then focus on robustness under perturbation in the framework of H ∞ control.
2.1 The H 2 synthesis and the Hamilton-Jacobi-Bellman PDE We consider the following infinite horizon optimal control problem: where y(t) = (y 1 (t), . . ., y d (t)) ⊤ ∈ R d denotes the state of the system, and the control signal u(•) belongs to U := L ∞ (R + ; R m ).The running cost is given by y R ≻ 0. We assume the system dynamics f (y) : R d → R d to be such that f (0) = 0, and the control operator g(y) : R d → R d×m to be C 1 (R d ).This formulation of the H 2 synthesis corresponds to the asymptotic stabilization of nonlinear dynamics towards the origin.By the application of the Dynamic Programming Principle, it is well-known that the optimal value function characterizing the solution of this infinite horizon control problem is the unique viscosity solution of the Hamilton-Jacobi-Bellman equation The explicit minimizer u * of eq. ( 2) is given by By inserting (3) into (2), we obtain the HJB equation to be understood in the classical sense.We recall that under the additional linearity assumptions f (x) = Ax with A ∈ R d×d and g(x) = B ∈ R d×m , the value function is known to be a quadratic form, V (x) = x ⊤ Πx, with Π ∈ R d×d positive definite, and eq.(HJB) becomes an Algebraic Riccati Equation for Solving for V (x) in eq.(HJB) globally in R d allows an online synthesis of the optimal feedback law (3) by evaluating the gradient ∇V (x) and g(x) at the current state x.This leads to an inherently robust control law in the sense that if the state of the system is perturbed to x ′ = x+δx, there still exists a stabilizing feedback action departing from the perturbed state, namely, u * (x ′ ).However, this control design neglects the modelling of the disturbance/uncertainties affecting the dynamics, with no general stabilization guarantees.We overcome this limitation by formulating an H ∞ synthesis, which we describe in the following.

2.2
The H ∞ problem and the Hamilton-Jacobi-Isaacs PDE We extend the previous formulation by considering a nonlinear dynamical system of the form where an additional disturbance signal w(•) ∈ W, with W = L 2 (R + ; R p ) enters the system through h(y) : R d → R d×p .We assume that y = 0 is an equilibrium of the system for u = w = 0.The H ∞ control objective is to achieve both internal stability of the closed-loop dynamics and disturbance attenuation through the design of a feedback law u = u(y) such that for a given γ > 0, and for all T ≥ 0 and w ∈ L 2 ([0, T [, R p ), the inequality holds.Here, P ∈ R p×p , P ≻ 0 , and y is the solution to (4).The parameter γ plays a crucial role in H ∞ control design.We say that the control system (4) has L 2 -gain not greater than γ, if (5) holds.Finding the smallest γ * for which this inequality is verified, also known as the H ∞ -norm of the system, is a challenging problem in its own right.The simplest yet computationally expensive method to find the H ∞ -norm of a system is through a bisection algorithm, as described in [11].Applying dynamic programming techniques leads to a characterization of the value function for this problem in terms of the solution of a Hamilton-Jacobi-Isaacs PDE valid for all γ ≥ γ * .The unconstrained minimizer and maximizer of (6), u * γ and w * γ respectively, are explicitly computed as where V γ (x) solves the Hamilton-Jacobi-Isaacs equation with and V γ (0) = 0.For an initial condition which is not a steady state we have see e.g.[43,Theorem 16].When there is no confusion, we denote V γ (x) by V .Note that if the disturbance attenuation is neglected by taking γ → ∞, we recover the solution of (HJB) instead.Moreover, under the linearity assumptions , where Π ∈ R d×d is positive definite, and eq.(HJI) becomes the following Algebraic Riccati Equation for solving the H ∞ problem for full state feedback.In the following section we discuss the construction of computational methods to synthesize nonlinear feedback laws inspired by the solution of HJB/HJBI PDEs.
3 Sub-optimal feedback control laws Despite the extensive parametrization of the control problem, eqns.(HJB) and (HJI) remain first-order, stationary nonlinear PDEs whose numerical approximation is a challenging task.These nonlinear PDEs are cast over R d , where d relates to the dimension of the state space of the dynamics, which can be arbitrarily large.In particular, if the system dynamics correspond to nonlinear PDEs, the direct solution of the resulting infinite-dimensional eqns.
(HJB) and (HJI) remains an open computational problem.We explore different alternatives to circumvent the solution of high-dimensional HJB PDEs by resorting to the sequential solution of the Algebraic Riccati Equations (ARE) and (ARE ∞ ) respectively, providing a tractable alternative for feedback synthesis in large-scale nonlinear dynamics.The different techniques we propose trade the optimality associated to the HJB-based control for computability, and hence will be referred to as suboptimal feedback laws.
For the sake of simplicity, we focus on the H ∞ synthesis.The H 2 feedback follows directly choosing h(x(t)) := 0 in (4).We assume a quadratic cost, i.e. ℓ(x) = x ⊤ Qx with Q ∈ R d×d symmetric positive semidefinite.

Linear-Quadratic Regulator (LQR) control law
The simplest suboptimal control law for nonlinear systems uses the optimal feedback law for the linear-quadratic control problem arising from linearization around an equilibrium, which we assume to be x = 0.For γ > γ * , we solve eq.
, and H = h(0).From the solution Π, we obtain the linear feedback law Such a feedback law, applied to the original nonlinear system, can only guarantee stabilization in a neighbourhood around the origin.

State-Dependent Riccati Equation (SDRE)
If we express the system dynamics through a space-dependent representation of the form ẋ we can synthesize a suboptimal feedback law inherited from (HJI) by following an approach known as the State-dependent Riccati Equation (SDRE).This method has been thoroughly analyzed in the literature, see, e.g.[14,5], and despite being purely formal, it is extremely effective for stabilizing nonlinear dynamics.The SDRE approach is based on the idea that infinite horizon optimal feedback control for systems of the form ( 12) are linked to a state-dependent ARE: where Solving this equation leads to a state-dependent Riccati operator Π(x), from where we obtain a nonlinear feedback law given by It is important to observe that the operator equation ( 13) admits an analytical solution only in a limited number of cases.More importantly, even if this solution is computed for every state x, the closed-loop differs from the optimal feedback obtained from solving (HJI), as the SDRE approach assumes that the value function is always locally approximated as V (x) = x ⊤ Π(x)x.From an optimal control perspective the SDRE can be interpreted as a model predictive control loop where at a given instant, the dynamics (A(x), B(x), H(x)) are frozen at the current state and an LQR feedback is numerically approximated.The procedure is illustrated in Algorithm 3.1.The resulting feedback is naturally different from the optimal nonlinear feedback law, and will remain different regardless of how fast the SDRE feedback is updated along a trajectory.The latter is also observed in [27], where it is shown that a direct derivation of the SDRE approach departing from the HJB PDE would lead to additional terms in (12).Notwithstanding, it is possible to show local asymptotic stability for the SDRE feedback.We recall the following result [5] concerning asymptotic stability of the SDRE closed-loop in the H 2 case.
Proposition 1 Assume a nonlinear system where x, and the pair (A(x), B(x)) is stabilizable for every x in a non-empty neighbourhood of the origin.Then, the closed-loop dynamics generated by the feedback law ( 14) are locally asymptotically stable.

Algorithm 3.1 SDRE-MPC loop
Require: Set u(t) := −K(x(t k ))x(t) 5: Integrate system dynamics to x(t k+1 ) 6: end for Assuming the stabilizability hypothesis above, the main bottleneck in the implementation of Algorithm 3.1 is the high rate of calls to an ARE solver for (13).Moreover, these ARE calls are expected to be sufficiently fast for real-time feedback control.This is a demanding computational task for the type of largescale dynamics arising in optimal control of PDEs.In the following, we discuss two variations the SDRE-MPC algorithm which mitigate this limitation by resorting to offline and more efficient online computations.

Offline approximation of the SDRE
A first alternative for more efficient SDRE computations was proposed in [6], inspired by a power series argument first discussed in [49].Assuming that the state operator A(x) can be decomposed into A(x) = A 0 + f 1 (x)A 1 , where A 0 and A 1 are constant matrices and f (x) is a scalar function, B(x) = B and H(x) = H, then the Riccati operator Π(x) solving ( 13) is approximated by where the matrices L n ∈ R d×d solve After solving a single ARE and N Lyapunov equations, an N −order approximation of Π(x) yields the feedback Unfortunately, the reduction above is only possible for a scalar nonlinearity.
If the nonlinearity that can be expressed as where A j ∈ R d×d and the state dependence is restricted to r scalar functions f j (x) : R d → R, then a first order approximation of Π(x) is given by where Π 0 solves (ARE ∞ ) for A 0 and the remaining Π j satisfy the Lyapunov equations with C 0 = A 0 − SΠ 0 , and Overall, this approach requires the computation of the ARE associated to A 0 in addition to r Lyapunov equations whose solution is fully parallelizable.Its implementation is summarized below.

Algorithm 3.2 Offline SDRE
1: Offline computations: 2: Compute Π 0 from (ARE∞) 3: for j = 1, . . ., r do 4: Compute Π j from the Lyapunov equation ( 22) 5: end for 6: Online process: 7: Set u(x(t)) := −R −1 B ⊤ Π(x(t))x(t) 8: Integrate system dynamics for x(t) An offline-online SDRE approach Although the previous approach is a valid variant to circumvent the online solution of Riccati equations at a high rate, it becomes unfeasible in cases where both d and r are large, as it requires the solution of r Lyapunov equations (22) and storage of the solution matrix with d 2 entries each.Such a large-scale setting arises naturally in feedback control of dynamics from semidiscretization of nonlinear PDEs and agent-based models.We present a variant of the offline SDRE approach which circumvents this limitation by resorting to an online phase requiring the solution of a single Lyapunov equation per step.
Let us define the quantity W (x) = r j=1 Π j f j (x).Multiplying each equation in (22) by its corresponding f j (x), it follows that W (x) satisfies the Lyapunov equation Therefore, the feedback law can be expressed as The feedback law (25) can be computed by solving an offline Riccati equation for Π 0 and an online Lyapunov equation for W (x); see section 4 for a discussion on the computational costs of the two approaches.The offline-online SDRE approach is summarized in Algorithm 3.3 below.

Algorithm 3.3 Offline-online SDRE
1: Offline phase: 2: Compute Π 0 from (ARE∞) 3: Online phase: 4: Integrate system dynamics to x(t k+1 ) 9: end for Remark 1 The approximation of the state space solution x(t k+1 ) at time t k+1 can be performed with both explicit or implicit time-stepping schemes.In both cases, the feedback gain remains frozen at K(x(t n )).

A preliminary test: the damped Sine-Gordon equation
We present a preliminary assessment of the effectiveness of Algorithm 3.1 and Algorithm 3.3 for the H 2 case.In section 5 we will focus on large-scale, twodimensional PDEs and H ∞ control.Given a domain Ω ⊂ R, we consider the control of the damped Sine-Gordon equation (see e.g.[42]) with homogeneous Dirichlet boundary conditions over Ω × R + 0 : where the control variable u(t) acts through an indicator function χ ωc (ξ) supported over ω c ⊂ Ω.The cost functional to be minimized is given by: where ω o := ∪ z i=1 ω oi ⊂ Ω represents a collection of local patches where we average the state.Defining y(t) = (X(•, t), Ẋ(•, t)) ⊤ , we write the dynamics as a first-order abstract evolution system ẏ(t) = Ay(t) + f (y(t)) + Bu(t) , where We approximate the operators above using a finite difference discretization in space.Given the domain Ω = [ξ L , ξ R ], we construct the uniform grid . We define the discrete state X i (t) := X(ξ i , t), and the discrete augmented state Y (t) = (X 1 (t), . . ., X d (t), Ẋ1 (t), . . ., Ẋd (t)) ⊤ .The discrete operators read where I d is the d × d identity matrix, and ∆ d is the discrete Laplace operator1 , To express the nonlinearity in a semilinear form consistent with (12) we define In our test we consider the following values for the given parameters, and In the cost functional (27) we set R = 1, ω o (x) = [−2.5,−1.5] ∪ [1.5, 2.5] and z = 2.In this test we take d = 402 nodes in the finite difference discretization.Time-stepping is performed with an implicit Euler method with a time step of 0.1.The small size Riccati and Lyapunov equations are solved using the Matlab functions icare and lyap, respectively.Controlled dynamics with different feedback controls are shown in Figure 1.The presence of the damping term α∂ t X(ξ, t) generates a stable trajectory for both the uncontrolled and LQR-controlled dynamics using the linearized feedback (11).However, we still observe differences in the state and control variables with respect to the SDRE controllers, SDRE-MPC and SDRE offline-online, described in Algorithms 3.1 and 3.3, respectively.The accumulated running costs in Figure 1 (top-left) indicate that the SDRE-MPC implementation achieves the best closed-loop performance, followed by SDRE offline-online, both outperforming linearized LQR and uncontrolled trajectories.However, the main difference between the SDRE-MPC and SDRE offline-online closed-loops is related to computational time.The SDRE-MPC solver requires the solution of multiple AREs in sequence, taking a total of 23 minutes of CPU time for this test, whereas the online solution of Lyapunov equations of the SDRE offline-online implementation reduces this computation time to 45 seconds.A deeper study on the methods performance in the large scale setting will be reported on in section 5, while implementation aspects associated with the solution of these algebraic equations are discussed in the next section.

Solving large-scale Algebraic Riccati/Lyapunov equations
The time steps discussed in Section 3.2 all use matrices that stem from the solution of algebraic matrix equations, and more precisely the (quadratic) Riccati and the (linear) Lyapunov equations.The past few years have seen a dramatic improvement in the effectiveness of numerical solution strategies for solving these equations in the large scale setting.For a survey in the linear case we refer the reader to the recent article [44], while for the algebraic Riccati equation we point the reader to, e.g., [7,45,10].
In our derivation we found projection methods to be able to adapt particularly well to the considered setting, with a similar reduction framework for both linear and quadratic problems; other approaches are reviewed for instance in [9,44].We emphasize that all considered methods require that the zero order term in the matrix equation, e.g., matrix Q in the Riccati equation (13), be low rank.
The general idea consists of first determining an approximation space K k that can be naturally expanded if needed, and then seeking an approximate solution in this space, by imposing a Galerkin condition on the matrix residual for computing the projected approximate solution.
Let V k be such that K k = range(V k ), with V k having orthonormal columns.Recalling that in both the Riccati and Lyapunov equations the solution matrix is symmetric, the approximate solution can be written as V k Y V ⊤ k .Consider the Riccati equation ( 13) for a fixed x = x(t * ), so that we set Imposing the Galerkin condition on R means that the residual matrix R be orthogonal to the approximation space, in the matrix sense, that is where the orthogonality of the columns of V k was used.If V k has small dimensions, the reduced Riccati matrix equation on the right also has small dimensions and can be solved by a "dense" method to determine Y X ; see, e.g., [10].The cost of solving the reduced quadratic equation with coefficient matrices of size k is at least 63 k3 floating point operations with an invariant subspace approach [35].Note that the large matrix V k Y X V ⊤ k is never constructed explicitly, since it would be dense even for sparse data.
Analogously, for the Lyapunov equation in (24), we can write This reduced Lyapunov equation can be solved by means of a "dense" method at a cost of about 15 k3 floating point operations for coefficient matrices of size k, if the real Schur decomposition is employed; see, e.g., [44].Note that the computational cost is significantly lower than that of solving the reduced Riccati equation with matrices of the same size.

On the selection of the approximation space
Choices as approximation space explored in the literature include polynomial and rational Krylov subspaces [44].They both enjoy the property of being nested as they enlarge, that is where k is associated with the space dimension.Rational Krylov subspaces have emerged as the key choice because they are able to deliver accurate approximate solutions with a relatively small space dimension, compared with polynomial spaces.Given a starting tall matrix V 0 and an invertible stable coefficient matrix A * , we have used two distinct rational spaces: the Extended Krylov subspace, ), which only involves matrix-vector products and solves with A * , and the (fully) Rational Krylov subspace, where σ j can be computed a-priori or adaptively.In both cases, the space is expanded iteratively, one block of vectors at the time, and systems with A * or with (A * − σ j I) are solved by fast sparse methods.For A * real valued and stable, the σ j 's are selected to be in C + , so that A * − σ j I is nonsingular.The actual choice of the shifts is a key step and a rich literature is available, yielding theoretically grounded effective strategies; see, e.g., the discussion in [44].
In our implementation we used the Extended Krylov subspace for solving the Lyapunov equation in Algorithm 3.3, which has several advantages, such as the computation of the sparse Cholesky factorization of A 0 once for all.On the other hand, we used the Rational Krylov subspace for the Riccati equation, which has been shown to be largely superior over the Extended Krylov on this quadratic equation, in spite of requiring the solution of a different (shifted) sparse coefficient matrix at each iteration [9,45].Nonetheless, in Section 4.2 we report an alternative approach that makes the Extended Krylov subspace competitive again for the Riccati problem with A 0 symmetric.Except for the operations associated with the reduced problems, the computational costs per iteration of the Riccati and Lyapunov equation solvers are very similar, if the same approximation space is used.
Although we refer to the specialized literature for the algorithmic details2 , we would like to include some important implementation details that are specific to our setting.In particular, the matrix C 0 employed in Algorithm 3.3 is given by C 0 = A 0 − BR −1 B ⊤ − 1 2γ 2 HP −1 H ⊤ Π 0 , which is not easily invertible if explicitly written, since it is dense in general.Note that A 0 is in general sparse, as it stems from the discretization of a partial differential operator.We can write Using the classical Sherman-Morrison-Woodbury formula, the product C −1 0 V for some tall matrix V can be obtained as which is assumed to be nonsingular.Therefore, C −1 0 V can be obtained by first solving sparse linear systems with A 0 , and then using matrix-matrix products.More precisely, the following steps are performed: We also recall that Π 0 , the solution to the initial Riccati equation, is stored in factored form, and this should be taken into account when computing matrix-matrix products with Π 0 .
While trying to employ the Rational Krylov space, we found that the structure of C 0 made the selection of optimal shifts {σ j } particularly challenging, resulting in a less effective performance of the method.Hence, our preference went for the Extended Krylov space above, with the above enhancement associated with solves with C 0 .

Feedback matrix oriented implementation
In Algorithm 3.1 the Riccati equation needs to be solved at each time step t n .However, its solution Π(x n ) is only used to compute the feedback matrix Hence, it would be desirable to be able to immediately compute K(x n ) without first computing Π(x n ).This approach has been explored in the Riccati equation literature but also for other problems based on Krylov subspaces, see, e.g., [41,33].
In this section, in the case where the matrix A(x n ) is symmetric, we describe the implementation of one of the projection methods described above, that is able to directly compute K(x n ) without computing Π(x n ) (in factored form), and more importantly, without storing and computing the whole approximation basis.The latter feature is particularly important for large scale problems, for which dealing with the orthogonal approximation basis represents one of the major computational and memory costs.To the best of our knowledge, this variant of the Riccati solver is new, while it is currently explored in [40] for related control problems and the rational Krylov space.
Here we consider the Extended Krylov subspace.For A(x n ) symmetric, the orthonormal basis of EK k can be constructed by explicitly orthogonalizing only with respect to the previous two basis blocks.Hence, only two previous blocks of vectors need to be stored in memory, and require explicit orthogonalization when the new block of vectors is added to the basis [41].This is also typical of polynomial Krylov subspaces constructed for symmetric matrices, giving rise to the classical Lanczos three-term recurrence [33].
With this procedure, in the reduced equation ( 31) the matrices computed as k grows by updating the new terms at each iteration, and the solution Y X can be obtained without the whole matrix V k being available.Note that the stopping criterion does not require the computation of the whole residual matrix, so that also in the standard solver the full matrix V k Y X V ⊤ k is never explicitly accessed.However, to be able to compute the basis V k appearing on the right still seems to be required.As already done in the literature, this problem can be overcome by a so-called "twostep" procedure: at completion, once the final Y X is available, the basis V k is computed again one block at the time, and the corresponding terms in the product Y X V ⊤ k are updated.Since A(x n ) is already factorized and the orthogonalization coefficients are already available (they correspond to the non-zero entries of A k ), then the overall computational cost is feasible; we refer the reader to [41] and its references for additional details for the two-step procedure employed for different purposes.

Large-scale nonlinear dynamical systems
In this section we present a numerical assessment of the proposed methodology applied to the synthesis of feedback control for two-dimensional nonlinear PDEs.The first test is a nonlinear diffusion-reaction equation, known as the degenerate Zeldovich equation, where the origin is an unstable equilibrium and traditional linearization-based controllers fail.The second test studies the viscous Burgers' equation with a forcing term.We discretize the control problem in space using finite differences, similarly as in Section 3.4.Controlled trajectories are integrated in time using an implicit Euler method, which is accelerated using a Jacobian-Free Newton Krylov method (see e.g.[32]).The goal of all our tests is the optimal and robust stabilization of the dynamics to the origin, encoded in the optimization of the following cost: This expression is similar to (27), but includes the The reported numerical simulations were performed on a MacBook Pro with CPU Intel Core i7-6, 2,6GHz and 16GB RAM, using Matlab [36].

Case study 1: the degenerate Zeldovich equation
We consider the control of a Zeldovich-type equation arising in combustion theory [20] over Ω × R + 0 , with Ω ⊂ R 2 and Neumann boundary conditions: The scalar control and disturbance act, respectively, through functions χ ωc (ξ) and χ ω d (ξ) with support ω c , ω d ⊂ Ω.The uncontrolled dynamics have three equilibrium points: Our goal is to stabilize the system to X ≡ 0, which is an unstable equilibrium point.A first step towards the application of the proposed framework is the space discretization of the system dynamics, leading to a finite-dimensional state-space representation.Following the setting presented in section 3.4, using a finite difference discretization leads to (34) where the discrete state X(t) = (X 1 (t), . . ., X d (t)) ⊤ ∈ R d corresponds to the approximation of X(ξ, t) at the grid points and the symbol • denotes the Hadamard or component-wse product.The matrix ∆ d ∈ R d×d is the finite difference approximation of the Neumann Laplacian and the matrices B, H ∈ R d are the discretization of the indicator functions supported over ω c and ω d , respectively.The discretization of (32) follows similarly as in section 3.4.Once the finite-dimensional state-space representation is obtained, we proceed to express the system in semilinear form (20) and implement the proposed algorithms.To set Algorithm 3.1 and Algorithm 3.3, from (34) we define where diag(v) denotes a diagonal matrix with the components of the vector v on the main diagonal, and decompose A(X) as where ∆ d is the two-dimensional discrete Laplacian and δ i,j denotes the Kronecker delta.In this test we set and the initial condition x(ξ, 0) = sin(ξ 1 ) sin(ξ 2 ), on a discretized space grid of n ξ1 × n ξ2 nodes with n ξ1 = n ξ2 = 101 (d = 10201).For the matrices B, H and C we considered a collection of patches depicted in Figure 2, and given by and In the following, we analyse results for the H 2 and H ∞ controls, i.e. γ = 0 and γ = 0, respectively, in (32).
Test 1. Experiments for H 2 -control.We start by presenting results for H 2control, i.e.P ≡ 0 in (32) and H ≡ 0 in (34).In Figure 3 we show a snapshot of the controlled trajectories at t = 3, a horizon sufficiently large for the dynamics to approach a stationary regime.In the top-left panel the uncontrolled problem reaches the stable equilibrium X ≈ 1.02.In the top-right panel we show the results of the LQR control computed by linearizing equation (34) around the origin, which also fails to stabilize around the unstable equilibrium X ≡ 0. The controlled solutions with Algorithm 3.1 and Algorithm 3.3 are shown at the bottom of the same figure: both algorithms reach the desired configuration.
The corresponding control input is shown in the top-right panel of Figure 4. We observe that the LQR control has a completely different behavior with respect to the control computed by Algorithm 3.1 or Algorithm 3.3.
The performance results of the different controlled trajectories is presented in Figure 4 where we show the evaluation of the cumulative cost functional in the top-left panel.As expected, Algorithm 3.1 provides the best closed-loop performance among the proposed algorithms.However, in terms of efficiency Algorithm 3.3 is faster than Algorithm 3.1 when increasing the dimension of the problem as shown in the bottom-left panel of Figure 4.When the dimension d increases (x-axis in the plot), the cost functional converges to 0.6 for Algorithm 3.1 and to 1 for Algorithm 3.3.Both methods are able to stabilize the problem.We omit reporting the behavior of the LQR-based control as it fails to stabilize the dynamics.To compare the proposed approaches we compute the H 2 − and H ∞ − controls with the same disturbance using both Algorithm 3.1 and Algorithm 3.3.The results are presented in Figure 5 and Figure 6.In every test case, the SDRE-based methodologies effectively stabilize the perturbed dynamics to a small neighbourhood around X ≡ 0.
A quantitative study is proposed in Figure 7 where we show the evaluation of the cumulative H 2 cost functional for both disturbances in the left panels.As expected, Algorithm 3.1 with H ∞ −control exhibits the best performance, closely followed by Algorithm 3.3.The right panels of Figure 7 show the different control inputs, which reflect the observed differences in closed-loop performance.

Test 3. On the use of a feedback matrix oriented implementation
To conclude the first case study we provide a numerical example where we synthesize the feedback operator K(x) = −R −1 B ⊤ Π(x) directly, circumventing the computation of Π(x), as explained in Section 4.2.For this test we The performance of the different feedback synthesis methods is presented in Figure 9.Here we show the evaluation of the cumulative cost functional in the top-left panel.For completeness we provide the control inputs on the top-right panel.We compare the performance of Algorithm 3.1 using Rational Krylov (RK) and Extended Krylov (EKSM) subspaces as the problem dimension in-   ary conditions: In this case, the scalar control and disturbance act, respectively, through the indicator function χ ωc (ξ), χ ω d (ξ) with ω c , ω d ⊂ Ω.A finite difference discretization of the space of the system dynamics leads to a state space representation of the form where the matrices ∆ d , D ∈ R d×d and B, H ∈ R d are finite-dimensional approximations of the Laplacian, gradient, control and disturbance operators, respectively, and the exponential term is understood component-wise.In particular, D is obtained using a backward finite difference discretization where B n := tridiag([ −1 1 0 ], n), and ⊗ denotes the Kronecker product.We proceed to express the semi-discretized dynamics in semilinear form.For this, we define and then For our numerical experiments we set Ω = [0, 1] × [0, 1], ǫ = 0.1, R = 0.05, and initial condition x 0 (ξ) = sin(ξ 1 ) sin(ξ 2 ), on a discretized space grid of n ξ1 × n ξ2 nodes with n ξ1 = n ξ2 = 101 (d = 10201).The matrices B and C are defined as in the previous study case (see Figure 2).In the following we discuss the results for H 2 and H ∞ control, i.e., P = 0 and P = 0, respectively, in (32).Finally, we discuss the results for P = 1, γ = 0.1 and disturbance w(t) = {0.1 sin(2t)} in (35).The results presented in Figure 12 are in line with our first case study.In this example, Algorithm 3.1 stabilizes the solution faster.Since it is difficult to visualize differences in the controlled state variables, we provide a qualitative analysis through Figure 12.We show the evaluation of the cost functional (32) on the left panel and the control inputs on the right panel.Again, we find that Algorithm 3.1 with γ = 0 has the lowest values for the cost functional.

Conclusions and future work
In this work we have discussed different alternatives for the synthesis of feedback laws for stabilizing nonlinear PDEs.In particular, we have studied the use of state-dependent Riccati equation methods, both for H 2 and H ∞ synthesis.Implementing an SDRE feedback law requires expressing the dynamics in semilinear form and the solution of algebraic Riccati equations at an arbitrarily high rate.This is a stringent limitation in PDE control, where highdimensional dynamics naturally emerge from space discretization.Hence, we study offline and offline-online synthesis alternatives which circumvent or mitigate the computational effort required in the SDRE synthesis.Most notably, we have proposed an offline-online method which replaces the sequential solution of algebraic Riccati equations by Lyapunov equations.Through extensive computational experiments, including two-dimensional nonlinear PDEs, we have assessed that the SDRE offline-online method provides a reasonable approximation of purely online SDRE synthesis, yielding similar performance results at a reduced computational cost.Moreover, the nonlinearities arising in nonlinear reaction and nonlinear advection PDE models can be easily represented within the semilinearization framework required by SDRE methods.In conclusion, SDRE-based feedback laws constitute a reasonable alternative for suboptimal feedback synthesis for large-scale, but well structured, nonlinear dynamical systems.Future research directions include the study of the SDRE methodology for high-dimensional systems arising from interacting particle systems, and the interplay with deep learning techniques to lower the computational burden associated to a real-time implementation.

3 Fig. 4
Fig. 4 Test 1. Top: Cumulative cost functional functional (left) with H 2 −control and corresponding control inputs (right).Bottom: CPU time for Algorithm 3.1 and Algorithm 3.3 (left), and convergence of the cost functional with respect to the dimension of the problem d (x-axis) (right).
creases.Specifically, we recall that the Extended Krylov subspace computes directly the gain matrix.As expected, Algorithm 3.1 provides the best closedloop performance among the proposed algorithms.However, in terms of CPU time Algorithm 3.3 is faster than Algorithm 3.1 when increasing the problem dimension d as shown in the bottom-left panel of Figure9.When the problem dimension d increases, the cost functional converges to 0.9 for Algorithm 3.1 and to 1.3 for Algorithm 3.3.The two different Krylov subspace methods in Algorithm 3.1 lead exactly to the same solution.In the plot legend we refer to Alg 3.1 RK for Rational Krylov subspaces and to Alg 3.1 EKSM for Extended Krylov subspaces.5.2 Case study 2: the viscous Burgers equation with exponential forcing termThe second experiment deals with the control of a viscous Burgers equation with exponential forcing term over Ω × R + 0 , with Ω ⊂ R 2 and Dirichlet bound-

Test 4 .
Experiments for H 2 −control.The trajectories of the controlled system problem with P = 0 in(32) are shown in Figure10.The uncontrolled solution tends to move towards the topright corner of Ω.All algorithms tend to control the solution to zero for large time, but at a different rate.The control inputs are then shown in the top-right of Figure11.Algorithm 3.1 has the largest control input, leading to the smallest cumulative

Fig. 9
Fig. 9 Test 3: cumulative functional (top-left) with H 2 −control and corresponding control input (top-right).Evaluation of the cost functional with respect to the dimension of the problem d (x-axis) (bottom-right) and CPU time (bottom-left) for both Algorithm 3.1 using Rational Krylov subspaces and Extended Krylov subspaces and Algorithm 3.3.

3 Fig. 11
Fig. 11 Test 4. Top: evaluation of the cumulative cost functional with H 2 −control (left) and control inputs (right).Bottom: cost functional for dynamics of increasing dimension for different algorithms (left) and CPU time for both Algorithm 3.1 and Algorithm 3.3 (right).