Overview
Atlassian Sourcetree is a free Git and Mercurial client for Windows.
Atlassian Sourcetree is a free Git and Mercurial client for Mac.
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 linearquadratic 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 smallscale (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 mexfiles), 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:
 MatlabXUnit from https://github.com/psexton/matlabxunit.
 CAREX from https://wwwuser.tuchemnitz.de/~benner/pub/Software/CAREX/carex_m.tar.gz.
The command runtests
should take about one minute and return success. You will see some warnings about illconditioned matrices being inverted, but it is ok.
Terminology
A CanBasis is a basis matrix for an ndimensional 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 illconditioned 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=QAinv(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.
License
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 MEXfile 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. 23, 548–570.
 [carex]
 Jörn Abels , Peter Benner CAREX  A Collection of Benchmark Examples for ContinuousTime Algebraic Riccati Equations (Version 2.0) (1999), http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.40.4899