[Math] Solving continuous algebraic Riccati equation using generalized eigenproblem algorithm

control theorylinear-controlmatricesoptimal control

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.

Best Answer

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.

Related Question