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.
Warning The 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)
Note The type of the initial value,
X0
orX0_lr
, dictates the type used for the whole trajectory,sol.X
andsol_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 Penzl1999Projection
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.
Warning The output matrix
C
of the included configuration differs from all the other configurations hosted at MOR Wiki by a factor of 10.