# DifferentialRiccatiEquations.jl

This package provides algorithms to solve autonomous Generalized Differential Riccati Equations (GDRE)

```
\left\{
\begin{aligned}
E^T \dot X E &= C^T C + A^T X E + E^T X A - E^T X BB^T X E,\\
X(t_0) &= X_0.
\end{aligned}
\right.
```

More specifically:

- Dense Rosenbrock methods of orders 1 to 4
- Low-rank symmetric indefinite (LRSIF) Rosenbrock methods of order 1 and 2, $X = LDL^T$

In the latter case, the (generalized) Lyapunov equations arizing in the Rosenbrock stages are solved using a LRSIF formulation of the Alternating-Direction Implicit (ADI) method, as described by LangEtAl2015. The ADI uses the self-generating parameters described by Kuerschner2016.

WarningThe low-rank 2nd order Rosenbrock method suffers from the same problems as described by LangEtAl2015.

The user interface hooks into CommonSolve.jl by providing the `GDREProblem`

problem type
as well as the `Ros1`

, `Ros2`

, `Ros3`

, and `Ros4`

solver types.

# Getting started

The package can be installed from Julia's REPL:

```
pkg> add git@gitlab.mpi-magdeburg.mpg.de:jschulze/DifferentialRiccatiEquations.jl.git
```

To run the following demos, you further need the following packages and standard libraries:

```
pkg> add LinearAlgebra MAT SparseArrays UnPack
```

What follows is a slightly more hands-on version of `test/rail.jl`

.
Please refer to the latter for missing details.

## Dense formulation

The easiest setting is perhaps the dense one,
i.e. the system matrices `E`

, `A`

, `B`

, and `C`

as well as the solution trajectory `X`

are dense.
First, load the system matrices from e.g. `test/Rail371.mat`

(see License section below)
and define the problem parameters.

using DifferentialRiccatiEquations
using LinearAlgebra
using MAT, UnPack
P = matread("Rail371.mat")
@unpack E, A, B, C, X0 = P
tspan = (4500., 0.) # backwards in time

Then, instantiate the GDRE and call `solve`

on it.

prob = GDREProblem(E, A, B, C, X0, tspan)
sol = solve(prob, Ros1(); dt=-100)

The trajectories $X(t)$, $K(t) := B^T X(t) E$, and $t$ may be accessed as follows.

sol.X # X(t)
sol.K # K(t) := B^T X(t) E
sol.t # discretization points

By default, the state $X$ is only stored at the boundaries of the time span `tspan`

,
as one is mostly interested only in the feedback matrices $K$.
To store the full state trajectory, pass `save_state=true`

to `solve`

.

sol_full = solve(prob, Ros1(); dt=-100, save_state=true)

## Low-rank formulation

Continuing from the dense setup, assemble a low-rank variant of the initial value, $X_0 = LDL^T$ where $E^T X_0 E = C^T C / 100$ in this case. Both dense and sparse factors are allowed for $D$.

using SparseArrays
q = size(C, 1)
L = E \ C'
D = sparse(0.01I(q))
X0_lr = LDLᵀ(L, D)
Matrix(X0_lr) ≈ X0

Passing this low-rank initial value to the GDRE instance
selects the low-rank algorithms and computes the whole trajectories in $X$ that way.
Recall that these trajectories are only stored iff one passes the keyword argument `save_state=true`

to `solve`

.

prob_lr = GDREProblem(E, A, B, C, X0_lr, tspan)
sol_lr = solve(prob_lr, Ros1(); dt=-100)

NoteThe type of the initial value,`X0`

or`X0_lr`

, dictates the type used for the whole trajectory,`sol.X`

and`sol_lr.X`

.

## Solver introspection / Callbacks

To record information during the solution process,
e.g. the residual norms of every ADI step at every GDRE time step,
define a custom observer object and associated callback methods.
Refer to the documentation of the `Callbacks`

module for further information.

julia> import DifferentialRiccatiEquations.Callbacks
help?> Callbackss

Note that there are currently no pre-built observers.

## ADI shift parameter selection

The ADI shifts may be configured using keyword arguments of `solve`

.

shifts = Shifts.Projection(2)
solve(::GALEProblem, ::ADI; shifts)
adi_kwargs = (; shifts)
solve(::GDREProblem, ::Ros1; adi_kwargs)
solve(::GAREProblem, ::NewtonADI; adi_kwargs)

Pre-built shift strategies include:

`Heuristic`

shifts described by Penzl1999`Projection`

shifts described by BennerKuerschnerSaak2014- User-supplied shifts via the
`Cyclic`

wrapper

Refer to the documentation of the `Shifts`

module for further information.

julia> import DifferentialRiccatiEquations.Shifts
help?> Shiftss

# Acknowledgments

I would like to thank the code reviewers:

- Jens Saak (https://github.com/drittelhacker)
- Martin Köhler (https://github.com/grisuthedragon)
- Fan Wang (https://github.com/FanWang00)

# License

The DifferentialRiccatiEquations package is licensed under MIT, see `LICENSE`

.

The `test/Rail371.mat`

data file stems from BennerSaak2005 and is licensed under CC-BY-4.0.
See MOR Wiki for further information.

WarningThe output matrix`C`

of the included configuration differs from all the other configurations hosted at MOR Wiki by a factor of 10.