3
$\begingroup$

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.

$\endgroup$
2
$\begingroup$

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.

$\endgroup$
5
  • $\begingroup$ Is that a better solution ? I tried ordschur in Octave but it's not exactly same command as in MATLAB. MATLAB has very smart built in flags which Octave don't have. $\endgroup$ – Daniel Mårtensson Dec 4 '17 at 22:22
  • $\begingroup$ @DanielMårtensson It's not better, it is necessary otherwise the order is random. I don't know which flags you mean. $\endgroup$ – percusse Dec 5 '17 at 0:09
  • $\begingroup$ If you know an answer with MATLAB code. You can post it here so I can accept is a answer :) Flags = keywords se.mathworks.com/help/matlab/ref/ordschur.html. Octave has no keywords. $\endgroup$ – Daniel Mårtensson Dec 5 '17 at 14:03
  • $\begingroup$ @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. $\endgroup$ – percusse Dec 5 '17 at 14:09
  • $\begingroup$ I don't have matlab. Sorry. $\endgroup$ – Daniel Mårtensson Dec 5 '17 at 14:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.