Exponential integrators for large-scale stiff Riccati differential equations

https://doi.org/10.1016/j.cam.2020.113360Get rights and content

Highlights

  • The matrix-valued versions of exponential integrators for stiff RDEs are proposed.

  • Forward and backward recursion implementations are presented.

  • Two typical schemes are considered and the low-rank implementations are exploited to deal with large-scale symmetric and nonsymmetric problems.

  • The performance is illustrated using comparisons with typical integrators.

Abstract

Riccati differential equations arise in many different areas and are particularly important within the field of control theory. In this paper we consider numerical integration for large-scale systems of stiff Riccati differential equations. We show how to apply exponential Rosenbrock-type integrators to get approximate solutions. Two typical exponential integration schemes are considered. The implementation issues are addressed and some low-rank approximations are exploited based on high quality numerical algebra codes. Numerical comparisons demonstrate that the exponential integrators can obtain high accuracy and efficiency for solving large-scale systems of 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 X(t)=AX(t)+X(t)D+QX(t)GX(t)F(X(t)),X(0)=X0,where ARM×M,DRN×N,QRM×N,GRN×M are given matrices and X(t)RM×N is the unknown matrix-valued function. RDEs of this form occur in many important applications such as optimal control, H-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 A and D 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 X(t)=AX(t)+X(t)AT+QX(t)GX(t),X(0)=X0,here, Q=QT,G=GT and X0=X0T. It is obvious that the solution of symmetric RDEs is symmetric as X(t)T 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 G=0, 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 MN-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 y=Ly+f(y),yRN,LRN×N,y(0)=y0.A general derivation of exponential integrators is based on the variation-of-constants formula y(tn+hn)=ehnLy(tn)+0hne(hns)Lf(y(tn+s))ds.By approximating the nonlinear terms f(y(tn+s)) 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 L 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 L 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 X(t)=etAX0etD+0te(tτ)AQe(tτ)Ddτ0te(tτ)AX(τ)GX(τ)e(tτ)Ddτ.

Proof

The proof can be done directly by differentiating both sides. 

In particular, over the time interval [tn,tn+hn], by using a change of variables t=tn+shn in (5) gives X(tn+hn)=ehnAX(tn)ehnD+h

Exponential Rosenbrock-type integrators for RDEs

In this section, we consider the time discretization of RDEs (1). Rewrite Eqs. (1) as X(t)=Sn(X)+Gn(X),where Sn(X) denotes the Fréchet derivative of F(X) and Gn(X) the nonlinear remainder at Xn, respectively: Sn(X)=AnX+XDn,Gn(X)=F(X)Sn(X)with An=AXnG and Dn=DGXn.

It is obvious that Sn is a Sylvester operator. Formally, by the variation of constants formula, the exact solution of (18) can be written as follows: X(tn+hn)=ehnSn(X(tn))+hn01e(1s)hnSn(Gn(X(tn+shn)))ds.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 An by An=AnGn0Dn(M+N)×(M+N).Using the formula (10.40) arising in [29], we have eAn=eAn01e(1s)AnGnesDnds0eDn(M+N)×(M+N).

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 LDTT and LDUT 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)

  • ReidW.T.

    Riccati Differential Equations

    (1972)
  • ButcherJ.C.

    Numerical Methods for Ordinary Differential Equations

    (2008)
  • BennerP. et al.

    Rosenbrock methods for solving Riccati differential equations

    IEEE Trans. Automat. Control

    (2013)
  • ChoiC.H. et al.

    Efficient matrix-valued algorithms for solving stiff Riccati differential equations

    IEEE Trans. Automat. Control

    (1990)
  • DieciL.

    Numerical integration of the differential Riccati equation and some related issues

    SIAM J. Numer. Anal.

    (1992)
  • StillfjordT.

    Low-rank second-order splitting of large-scale differential Riccati equations

    IEEE Trans. Automat. Control

    (2015)
  • StillfjordT.

    Adaptive high-order splitting schemes for large-scale differential Riccati equations

    Numer. Algorithms

    (2018)
  • GüldoǧanY. et al.

    Low rank approximate solutions to large-scale differential matrix Riccati equations

    (2017)
  • HachedM. et al.

    Approximate solution to large nonsymmetric differential Riccati problems

    (2018)
  • KoskelaA. et al.

    A structure preserving Krylov subspace method for large scale differential Riccati equations

    (2017)
  • HochbruckM. et al.

    Explicit exponential Runge–Kutta methods for semilinear parabolic problems

    SIAM J. Numer. Anal.

    (2006)
  • Cited by (11)

    View all citing articles on Scopus
    View full text