# PGDoubling

## a MATLAB package to solve algebraic Riccati equations and optimal control problems using permuted graph bases

© Federico Poloni, 2012 -- 2014

## Purpose

This package uses the algorithms described in [MehP12, MehP13] to
compute Lagrangian invariant subspaces of Hamiltonian and "extended
Hamiltonian" pencils. The resulting algorithm is more accurate than
Matlab's care on several tricky test problems [carex].

It is the right choice for you if:

• you want to compute a
linear-quadratic regulator, solve a Riccati equation, or a
generalized Riccati equation (those that correspond to a matrix pencil
in the form:
    [ 0  I 0]    [0  A  B]
[ -I 0 0]s - [A' Q  S]
[ 0  0 0]    [B' S' R]


);

• you are not satisfied with the accuracy obtained by Matlab's care, or wish to test a new algorithm;

• your problem is small-scale (up to a few hundred degrees of freedom);

• your system is not a descriptor system (i.e., no E coefficient,
just \dot{x} = Ax + Bu).

• Note: your problem could also be a singular control problem (Lur'e
equations, singular R); there is also a function to compute the
optimal parameter gamma in H_infinity control;

The current implementation is in pure Matlab (no mex-files), so it
is slower that care. Rewriting the more intensive parts in Fortran
or C will very likely result in a huge speedup.

## Tests

The files named test*.m are tests. In order to run them, you need to install and have in your Matlab path (change it using pathtool) the following two libraries:

The command runtests should take about one minute and return success. You will see some warnings about ill-conditioned matrices being inverted, but it is ok.

## Terminology

A CanBasis is a basis matrix for an n-dimensional subspace in the form
P*[I(n);X], where P is a permutation matrix and X is arbitrary.
Usually the routines return CanBases in which each element of X is
smaller (in modulus) than a certain configurable threshold.
A CanBasis is stored as a struct can with fields can.X and can.p (permutation vector).

A SymBasis is the structured version of a CanBasis for a Lagrangian
subspace. P is a "symplectic swap matrix" (a symplectic analogous of a
permutation), and X is Hermitian.
A SymBasis is stored as a struct sym with fields sym.X and sym.v (boolean vector).

You can represent subspaces, matrices and matrix pencils through
canBases and symBases. See the basic theory in [MehP12]. (all will be
clearer if you have read this paper.)

## Documentation

All functions contain help; type help functionName in Matlab to access
it. You can look in the "test" functions for some more usage examples.
In the following there is a quick introduction.

## Basic usage: matrix equations

    [X,Y,U,V]=solveCARE(A,G,Q,...options...)


Returns the symmetric (semi)stabilizing solution of a CARE
A'X+XA+Q=XGX. Also returns the (semi)stabilizing solution of the "dual
CARE" YA'+AY+YQY=G. More important, U and V are basis matrices for the
(semi)stabilizing and (semi)antistabilizing Lagrangian subspaces of the
Hamiltonian, i.e., subspaces U=[U1;U2] and V=[V1;V2] such that
X=U2/U1, Y=V1/V2. You should always try to restate your successive
computations in order to use U and V rather than X and Y
, as X and Y
may be ill-conditioned while U and V are tame.

Recommended options for robustness:
solveCARE(A,G,Q,'type','sign','safer',true). If your problem has large
size, you may consider increasing the threshold with
(...'threshold',sqrt(n),'diagonalThreshold',sqrt(n)). The impact of
the thresholds hasn't been investigated thoroughly yet.

[X,Y,U,V]=solveECARE(A,B,Q,R,S,...options...)


solve a (possibly singular) generic control problem in the form:

    [ 0  I 0]    [0  A  B]
[ -I 0 0]s - [A' Q  S],
[ 0  0 0]    [B' S' R]


returning the (semi)stabilizing X and Y. In addtion, returns bases U and V
for the stabilizing subspace [X;I] and the antistabilizing one [I;Y]. Note that the position of X and I in the subspace is different from the one in the 2x2 problem. I am just following the standard notation for each of the two problems. In any case, you should work with U and V rather than with X and Y, as argued above.

    gamma=gammaIteration(A1,B1,B2,C1,C2,D11,D12,D21,D22,tol,...options...);


run the gamma iteration [BenBMX07,MehP13] on a H_inf control problem,
with tolerance tol.

Some more matrix equations could be added --- These techniques work also
for DARE, X=Q+A*inv(X)A', X=Q-Ainv(X)*A', and 0=P+QY+RY\^2 (with
some modifications). If you would like to solve one of these with this package, contact the author.

## Basic usage: building blocks

canBasisFromSubspace takes a subspace and returns X and a permutation P
defining its CanBasis. symBasisFromSymplecticSubspace takes a Lagrangian
subspace and returns its SymBasis. The associated permutation matrix
P is obtained by rowSwap(eye(2*n),sym.v,'N'), if needed.

You may also be interested in symBasisFromSymplecticPencil,
symBasisFromHamiltonianPencil, symBasisFromEvenPencil (actually only works
for even pencils in the form above; a generic even pencil does not have
a symBasis).

Some random stuff that an end user might find interesting:
checkCAREInvariantSubspaceResidual (computes and returns several error
measures), riccatiCoefficientsFromHamiltonian, hamiltonian (should be
named hamiltonianFromRiccatiCoefficients), pirq ((symplectic swap
P)RQ factorization), ChuLM07Carex (returns the results of the
experiments shown in [ChuLM07]), randomLagrangianSubspace.

## Feedback

This library does not have a large user base; I am happy to hear about end users' experiences ("I tried it and it works great; I tried it and it doesn't
work") and to receive feedback and bug reports.

If you use it in a paper, please cite us when it is relevant. In the academia we basically live off citation counts.

Check COPYING.txt to see if you qualify to use, modify and redistribute the library in other projects for free.

## FAQ

Q: why are the function names so long?

A: I find it more expressive than the default 'shorten everything' Matlab syntax. You can configure Matlab's editor for tab completion and the fact that the names are long will not slow you down.

Q: why is it so slow?

A: mainly, because the PiRQ factorization is implemented with a loop in pure Matlab. Rewriting it in a MEX-file would result in a massive speed improvement on large problems. In general, the functions are not aggressively optimized; some more performance could be squeezed out.

Q: why do you roll your own option parsing functions? OOP is slow in Matlab!

A: I haven't found yet an option parser for Matlab that satisfies some basic
requirements, namely: (1) can pass options unchanged to subroutines (2)
complains at the end if there are unused parameters (3) has an
acceptable syntax which does not require me to specify option names in
18 different places. If you know of one, let me know.

Q: does it work when there is a mass matrix E?

A: no, and it would require significant changes to the algorithm to
support it. We'll be working on that hopefully.

## References

[MehP12]
Mehrmann, Volker; Poloni, Federico Doubling algorithms with
permuted Lagrangian graph bases.
SIAM J. Matrix Anal. Appl. 33
(2012), no. 3, 780–805.
[MehP13]
Mehrmann, Volker; Poloni, Federico Using permuted graph bases in H∞
control.
Automatica J. IFAC 49 (2013), no. 6, 1790–1797.
[ChuLM07]
Chu, Delin; Liu, Xinmin; Mehrmann, Volker A numerical method for
computing the Hamiltonian Schur form.
Numer. Math. 105 (2007), no.
3, 375–412.
[BenBMX07]
Benner, Peter; Byers, Ralph; Mehrmann, Volker; Xu, Hongguo A robust
numerical method for the γ-iteration in H∞ control.
Linear Algebra
Appl. 425 (2007), no. 2-3, 548–570.
[carex]
Jörn Abels , Peter Benner CAREX - A Collection of Benchmark Examples for Continuous-Time Algebraic Riccati Equations (Version 2.0) (1999),
http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.40.4899