# Solving continuous algebraic Riccati equation using generalized eigenproblem algorithm

I want to solve the continuous algebraic riccati equation:

$$A^TX + XA - XBR^{-1}B^TX + Q = 0$$

To solve this, I have been using Schur's method to solve algebratic riccati equations.

Schur's method to solve algebratic riccati equations

The problem is the solution $X$ is not the same solution when I use MATLAB's command X = CARE(A, B, Q, R). MATLAB is using Generalized Eigenproblem Algorithms and Software for Algebraic Riccati Equations.

Generalized Eigenproblem Algorithms and Software for Algebraic Riccati Equations,'' Proc. IEEE, 72(1984), 1746--1754.

So my question is how I can solve the algebraic riccati equation using the generalized eigenproblem algorithm?

Assume that we have those matrices.

A =

0     1
-2    -3

>> B

B =

0
1

>> Q

Q =

1     0
0     2

>> R

R =

1


Using Schur's method. We will found the solution $X$ by using this:

>> H = [A -B*inv(R)*B'; -Q -A']; % Hamiltonian matrix
>> [U, S] = schur(H); % Schur decomposition
>> [m,n] = size(U);
>> U11 = U(1:(m/2), 1:(n/2));
>> U21 = U((m/2+1):m, 1:(n/2));
>> X = U21*inv(U11) % The solution

X =

-6.0000   -2.8074
-8.1926   -3.0000

>> A'*X + X*A - X*B*inv(R)*B'*X + Q % Check if this equation is near 0

ans =

1.0e-13 *

0.0711    0.0888
0.1776    0.1155


Yes it is! Now we use MATLAB's command CARE which use the generalized eigenproblem algorithm.

>> X = care(A, B, Q, R)

X =

1.5737    0.2361
0.2361    0.3871

>> A'*X + X*A - X*B*inv(R)*B'*X + Q  % Check if this equation is near 0

ans =

1.0e-14 *

-0.1332   -0.0208
-0.0208    0.0888


Let's try Octave/MATLAB's command fsolve.

>> fun = @(P, A, B, Q, R) P*A+A'*P-P*B*inv(R)*B'*P + Q;
>> p0 = 1*ones(2);
>> X = fsolve(@(P) fun(P, A, B, Q, R), p0)
X =

1.57368   0.23607
0.23607   0.38705

>> A'*X + X*A - X*B*inv(R)*B'*X + Q
ans =

-1.7419e-08  -1.1153e-08
-1.1153e-08  -6.9926e-09


So both give good solutions. The only difference was that Schur's method give the solution $X$ as matrix only, but with the generalized eigenproblem algorithm, we got a positive semi-definite solution $X$. The fsolve command works very good and give a equal solution as MATLAB's command CARE do.

Let's change matrix A.

>> A = [0 1; -3 -5]; % The other matrices are the same
>> X = care(A, B, Q,R)

X =

0.9118    0.0215
0.0215    0.1994


Let's use fsolve with the new matrix A.

X = fsolve(@(P) fun(P, A, B, Q, R), p0)
X =

1.53014   0.16228
0.16228   0.22729

>> A'*X + X*A - X*B*inv(R)*B'*X + Q
ans =

-9.9554e-08  -5.1336e-08
-5.1336e-08  -2.1732e-08


So here we got a change of the solution X, but the solution X has still positive values.

The problem with your solution is that you don't sort the stable eigenbasis and pick that up for the solution of $X$. This step is crucial for detecting whether there is a viable solution for ARE or not.
Hence either you add an ordschur after the Schur decomposition or you perform the same with a $QZ$ decomposition via ordqz. Otherwise you are not guaranteed to pick up the stabilizing solution by no means.
• @DanielMårtensson It is in gcare.m lines 116-117. Type edit gcare.m in the editor. I've mentioned this to you before. I don't use matlab myself anymore. – percusse Dec 5 '17 at 14:09