Exponential integrators for large-scale stiff Riccati differential equations
Introduction
In this paper we are concerned with numerical methods for large-scale systems of stiff Riccati differential equations (RDEs) of the following form where are given matrices and is the unknown matrix-valued function. RDEs of this form occur in many important applications such as optimal control, -control, filtering, boundary value problems for systems of ODEs and many others (see e.g., [1], [2], [3], [4]). In most control problems, the coefficient matrices and are obtained from the discretization of operators defined on infinite dimensional spaces, and the fast and slow modes exist. This means that the associated RDEs will be fairly large and stiff.
An important special case of (1) is the symmetric RDEs here, and . It is obvious that the solution of symmetric RDEs is symmetric as is also a solution. Symmetric RDEs are possibly the most widely studied equations due to their importance in linear–quadratic optimal control problems. Another special mention should be made of RDEs (1) with , yielding the so-called matrix Sylvester differential equations (SDEs). For a thorough description of these equations and some qualitative issues, we refer the reader to [1], [4], [5], [6], [7] and the references appearing therein.
Many numerical methods have been developed in the past for solving RDEs. Perhaps the most natural numerical technique is to rewrite (1) as an -vector ODEs based on Kronecker product, and then to use a standard numerical integrator such as Runge–Kutta or linear multi-step solver [8]. However, these approaches are not suitable for the solution of large stiff RDEs. They are generally computationally expensive and it is hard to exploit the structure inherited in some large practical problems. For stiff RDEs, some matrix-valued versions of implicit time integration schemes, such as BDF, Rosenbrock methods have been explored through a direct time discretization of (1), see e.g., [9], [10], [11]. Recently, some other unconventional numerical methods have been also developed for RDEs and related problems, including splitting methods [12], [13], [14] and projection methods [15], [16], [17], etc.
The aim of this paper is to introduce exponential integrators for large-scale stiff RDEs. In the past two decades exponential integrators have become a popular tool for solving large-scale stiff semi-linear systems of ODEs A general derivation of exponential integrators is based on the variation-of-constants formula By approximating the nonlinear terms by an appropriate algebra polynomial, various types of exponential integrators can be exploited. Different approximations result in different types of exponential integrators of either multi-step type or Runge–Kutta type, see e.g., [18], [19], [20], [21]. A main advantage of exponential integrators with stiff order is that they do not suffer from an order reduction even if the matrix is a discretization of a unbounded linear operator. For a full overview of exponential integrators and associated software, we refer the readers to the reviews [22], [23] and references therein. Although the matrix differential equations (1) can be reformulated as the form (3) and solved by an exponential integrator, this approach will generate very large and not be appropriate.
In the present paper we propose matrix-valued versions of exponential integrators for stiff RDEs (1). The methods provide an efficient alternative to implicit integrators for computing the solution of RDEs. For large-scale system, in many applications it is often observed, both practical and theoretical, the solution has low numerical rank and can be approximated by products of low-rank matrices [24], [25]. To utilize such structure, we introduce how low-rank approximations can be applied to exponential integrators. Thus we are able to save computational and memory storage requirements compared to a straightforward application of exponential integrators.
The remainder of the paper is organized as follows: In Section 2, we give some basic results and properties of RDEs. In Section 3, the exponential Rosenbrock-type methods are introduced for the application to the RDEs. In Section 4 we show some issues of implementation and the low-rank approximations for both typical exponential integration schemes are exploited. Section 5 is devoted to some numerical examples and comparisons with splitting methods of similar orders. Finally, we draw some conclusions in Section 6.
Section snippets
Preliminaries
We start with recalling a general result on the solution of the RDEs. The following result shows that the RDEs (1) can be equivalently written in an integral form (see e.g., [26]).
Theorem 1 The exact solution of the RDEs (1) can be written in integral form
Proof The proof can be done directly by differentiating both sides. □
In particular, over the time interval , by using a change of variables in (5) gives
Exponential Rosenbrock-type integrators for RDEs
In this section, we consider the time discretization of RDEs (1). Rewrite Eqs. (1) as where denotes the Fréchet derivative of and the nonlinear remainder at , respectively: with and .
It is obvious that is a Sylvester operator. Formally, by the variation of constants formula, the exact solution of (18) can be written as follows: The above expression has a
General implementation
For exponential integrators, the main computational cost is to approximate the exponential and exponential-type operators at each time step. To our knowledge, there is no explicit method to evaluate the Sylvester operator exponential in the literature. For the computation of the first -function, the following formula gives an indirect way. Define the augmented matrix by Using the formula (10.40) arising in [29], we have
Numerical experiments
In this section, we present some numerical experiments to assess the performance of exponential integration methods. We report the numerical experience of the second-order scheme ExpEuler and the third-order scheme Erow3. All experiments are performed under Windows 10 and MATLAB R2018b running on a laptop with an Intel Core i7 processor with 1.8 GHz and RAM 8 GB. Unless otherwise stated, we use the relative errors at the final time, measured in the Frobenius norm.
For the sake of simplicity we
Conclusion
In this paper we have introduced matrix versions of exponential integrators for approximating the solution of the RDEs. We have discussed the construction of the schemes and studied their performance on several test problems. We have presented an efficient implementation based on and factorizations to deal with large-scale symmetric and nonsymmetric problems. The schemes possess several favorable properties. They preserve the equilibria of the RDEs and are explicit, involving only
Acknowledgments
We thank the referees for suggestions that lead to an improved paper and the author of [14] for sharing his MATLAB codes. This work was supported in part by the Jilin Scientific and Technological Development Program, PR China (Grant Nos. 20180101224JC and 20200201276JC) and the Natural Science Foundation of Jilin Province, PR China (Grant Nos. 20190499kJ and 20200822KJ), and the Scientific Startup Foundation for Doctors of Changchun Normal University, PR China (Grant No. 002006059).
References (42)
- et al.
Remarks on the time-varying Riccati equations
Systems Control Lett.
(1999) - et al.
Convergence of the solution of a nonsymmetric matrix Riccati differential equation to its stable equilibrium solution
J. Math. Anal. Appl.
(2006) Global existence and stability of solutions of matrix Riccati equations
J. Math. Anal. Appl.
(2001)- et al.
Numerical low-rank approximation of matrix differential equations
J. Comput. Appl. Math.
(2018) - et al.
On the benefits of the factorization for large-scale differential matrix equation solvers
Linear Algebra Appl.
(2015) - et al.
Low-rank Newton-ADI methods for large nonsymmetric algebraic Riccati equations
J. Franklin Inst. B
(2016) - et al.
Solution of a nonsymmetric algebraic Riccati equation from a two- dimensional transportmodel
Linear Algebra Appl.
(2011) - et al.
Matrix Riccati Equations in Control and Systems Theory
(2003) - et al.
Numerical Solution of Boundary Value Problems for Ordinary Differential Equations
(1988) Introduction to Control Theory
(1993)
Riccati Differential Equations
Numerical Methods for Ordinary Differential Equations
Rosenbrock methods for solving Riccati differential equations
IEEE Trans. Automat. Control
Efficient matrix-valued algorithms for solving stiff Riccati differential equations
IEEE Trans. Automat. Control
Numerical integration of the differential Riccati equation and some related issues
SIAM J. Numer. Anal.
Low-rank second-order splitting of large-scale differential Riccati equations
IEEE Trans. Automat. Control
Adaptive high-order splitting schemes for large-scale differential Riccati equations
Numer. Algorithms
Low rank approximate solutions to large-scale differential matrix Riccati equations
Approximate solution to large nonsymmetric differential Riccati problems
A structure preserving Krylov subspace method for large scale differential Riccati equations
Explicit exponential Runge–Kutta methods for semilinear parabolic problems
SIAM J. Numer. Anal.
Cited by (11)
Computing the Lyapunov operator φ-functions, with an application to matrix-valued exponential integrators
2022, Applied Numerical MathematicsCitation Excerpt :Although MDEs (7) can be reformulated as a standard (vector-valued) ordinary differential equation and solved by a standard exponential integrator, this approach will usually be memory consuming as well as computationally expensive. Recently, in [30], the matrix-valued exponential Rosenbrock-type methods are proposed for solving DREs. The important ingredient to implementation of vector-valued exponential integrators is the computation of the matrix φ-functions.
Low-Rank Exponential Integrators for Differential Riccati Equation
2023, Frontiers in Artificial Intelligence and Applications