[Math] analytic solution to structured algebraic Riccati equation

control theorylinear algebraoptimization

In solving a model I have written down for a research paper, I am left with two Algebraic Riccati Equations, that is I need to solve for $X$ in the equation

\begin{align*}
X = A^\top (X + XB(R + B^\top X B)^{-1}B^\top X))A + Q.
\end{align*}

These arise from solving a LQG problem. From what I have picked up over the last few days, there is no general solution to this equation, i.e. an $f$ such that $X = f(A,B,R,Q)$ for all $(A,B,R,Q)$. However, there seems to be many special instances in which substantial analytical progress can be made.

My problem has quite a bit of structure, and I wonder if I could be pointed to some results that would help me. My problem is characterized by:

(1) $B = (1,1,0,0)$ for one problem and $B = (0,0,1,0)$ for the other, (2) $R$ a positive scalar, (3) $Q$ symmetric, and (4) (for one problem) $A$ also symmetric and $Q$ diagonal.

If anyone is aware of general results or methods for approach that might be useful in my particular case, I would be very grateful. It would be nice if these could be stated with a minimum of engineering jargon, my background is more in math and some of the terminology is not helpful.

Best Answer

There is no analytic method and for real-world problems in control theory, you do not need an analytic solution.

Here is the solution suggested from Vasile Sima which I wrote the C++ codes based upon here.

/*
This file contains a continuous-time algebraic Riccati equation solver based on the explanations from Dr. Vasile Sima.
Further reference: Sima, Vasile. Algorithms for linear-quadratic optimization. Vol. 200. CRC Press, 1996.
*/

// A C++ implementation of a Continuous-time Algebraic Riccati Equation (CARE) solver based on Schur vectors.

// Solves a CARE, A'X + XA - XBR^(-1)B'X + Q = 0, using Schur vectors method, with A and Q n-by-n, B n-by-m, and R m-by-m real matrices, and Q and R symmetric, R positive definite. It is assumed, e.g., that the matrix pair (A,B) is stabilizable, and the matrix pair (C,A) is detectable, where Q = C'C, and rank(C) = rank(Q).

// Main steps:

// 1. Construct the Hamiltonian matrix H, H = [ A -B*(R\B'); Q -A' ]; (in MATLAB notation).

// 2. Reduce H to a real Schur form S and accumulate the transformations in U, S = U'HU, U orthogonal.

// 3. Reorder the real Schur form S so that all n stable eigenvalues are moved to the leading part of S, accumulating the transformations in U.

// 4. Compute the unique positive-semidefinite stabilizing solution X, X = U(n+1:2n,1:n)*U(1:n,1:n)^(-1).

// 5. Symmetrize X, X = ( X + X' )/2.

If you are using MATLAB, you can simply use care, dare and even lqg commands.

An ad-hoc solution is using Genetic Algorithm to solve these problems via C++ (solver, NLP demo) or MATLAB.

Related Solutions

[Math] A System of Matrix Equations (2 Riccati, 1 Lyapunov)

I posting it as another answer because I don't want to delete anything in previous yet and the other option is to make it too long.

4) The equation for $X$ is Sylvester equation, but you right this is closely related to Lyapunov. Lyapunov equation is a special case of Sylvester equation.

5) Remove linear algebra tag, the problem is non linear.

6) Add optimization and oprimal control tags.

7) Since we figured out that the problem with your code is well known numerical instability of matrix inverse (I should have noticed it instantly). I tried to write a code which bypasses the problem, i.e. do not invert matrices explicitly. I've ended with the following code which also give more hints on where to look for additional conditions about the matrices. Hopefully these conditions will help to find to prove the conjecture.

Briefly about the algorithm,

use matlab's algorithm dare to get the ricatti matrices $P$ and $W$
To obtain the values of $K$ and $L$ use zero finding algorithm (fsolve) using their formula in ($1$) and ($2$) respectively.
Use matlab algorithm dlyap to obtain $X$
To verify the result solve perturb the required equation so that there no matrix inverse, i.e. $${\bf K}=(\gamma{\bf I}_{n}+(1-\gamma) {\bf X} {\bf W}^{-1}){\bf L}$$ $${\bf K}=(\gamma{\bf I}_{n}+(1-\gamma) {\bf X} {\bf W}^{-1}) {\bf W} {\bf H}({\bf H}^\top {\bf W} {\bf H}+ {\bf R})^{-1} $$ $${\bf K}({\bf H}^\top {\bf W} {\bf H}+ {\bf R})=(\gamma{\bf I}_{n}+(1-\gamma) {\bf X} {\bf W}^{-1}) {\bf W} {\bf H} $$ $${\bf K}({\bf H}^\top {\bf W} {\bf H}+ {\bf R})=(\gamma{\bf W}+(1-\gamma) {\bf X} ) {\bf H} $$

possible hints for the conjecture The algorithm is more stable with respect to matrix inverse. The finding of $K,L$ is somewhat a problem thought because the initial guess is not guided. In general, it may be essential difficulty in arriving to correct values. However when I run the algorithm $1000$ times I get $>600$ fails of dare and only few fails of fsolve. I think dare may shed some light on the original problem.

Thus, I've found that dare is very capricious algorithm. More precisely it doesn't has a solution for any input, but for narrowed family of matrices. The exact conditions as described in matlab aren't clear to me, so I even asked another question with hope that somebody will give an answer. I've found several unanswered questions about algebraic Ricatti equations, so I would suggest to go to literature to discover under what conditions to $F,Q,H,R$ the condition that the solution to Ricatti exists. The same about Silverster\Lyapunov. I would also investigate the link between algebraic and differential equations, perhaps the reduction to DE make the prove easier, so worth a try.

The code

n=9;
r=5;

F=round(10*rand(n),0);
H=round(10*rand(n,r),0);
R= round(10*rand(r),0);
R=R+R';
Q=round(10*rand(n),0);
Q=Q+Q';

P = dare(F',H,Q,R/gamma);
W = dare(F',H,Q,R);

options=optimset('Algorithm',{'levenbergmarquardt',0.0005});%,'Display','iter');
[K,fvalK,exitflagK] = fsolve(@(x) F*(P-x*H'*P)*F'+Q-P,zeros(n,r),options);
[L,fvalL,exitflagL] = fsolve(@(x) F*(W-x*H'*W)*F'+Q-W,zeros(n,r),options);

if exitflagK<=0 || exitflagL<=0 
        error('try again') 
end

X = dlyap((I-K*H')*F,(I-H*L')*F',K*H'*W);

tmp = K*(H'*W*H+R)-(gamma*W+(1-gamma)*X)*H;   
err=norm(tmp(:),inf);
disp(['K(H''W*H+R) - (gamma*W+(1-gamma)*X)*H) = ', num2str(err)])

[Math] Algebraic Riccati Inequality Solution via LMI

Your question is incorrect to begin with, as the continuous-time ARE is

$\quad A^TP + P A - PBR^{-1}B^TP + Q = 0$

Hence, just as bad form.

The LMI formulations of LQ in both discrete-time and continuous-time would typically be done in both the Riccati matrix and the feedback matrix.

Find minimum-trace $P$ such that

$ (A-BK)^T P (A-BK) - P + K^TRK + Q \preceq 0$

Multiply with $P^{-1}$ from left and right, define $Y = P^{-1}$ and $F = KP^{-1}$, and apply suitable Schur complements to arrive at linear form.

Best Answer

Related Solutions

[Math] A System of Matrix Equations (2 Riccati, 1 Lyapunov)

[Math] Algebraic Riccati Inequality Solution via LMI

Related Question