I want to solve the Algebraic Riccati Equation via MATLAB or Octave.

The first thing I do is to create the Hamiltonian matrix of A, B, Q and R.

> H = [A -(B.*inv(R).*B'); -Q -A']

Then I do the Schur decomposition of H.

> [U, S] = schur(H)

Then I will find U21 and U11 because U is a square matrix.

> [m,n] = size(U)
> U11 = U(1:(m/2), 1:(n/2))
> U21 = U((m/2+1):m, 1:(n/2))

Then I find the solution to Algebraic Riccati Equation.

> X = U21*inv(U11)

The problem is this method does not give the same solution values for matrix X, when using MATLAB's function X= care(A, B, Q, R).

Why? Have I done this method wrong ?

Source: http://dspace.mit.edu/bitstream/handle/1721.1/1301/R-0859-05666488.pdf

  • 2
    $\begingroup$ From a first glance the following should be corrected: in the definition of the Hamiltonian $H$ the (1,2) block is $-B*inv(R)*B'$ and does not involve an element by element multiplication as described by $-B.*inv(R).*B'$. However, this might only be a typo in your question since $-B.*inv(R).*B'$ is not defined in most cases (dimensions do not agree). Maybe if you provide your example we can get more insight. $\endgroup$ – RTJ Oct 24 '17 at 22:41
  • 2
    $\begingroup$ @CTNT Is it also faster and more accurate to use B/R instead of B*inv(R). $\endgroup$ – Kwin van der Veen Oct 25 '17 at 0:19
  • $\begingroup$ @KwinvanderVeen I have changed the question now. $\endgroup$ – Daniel Mårtensson Oct 25 '17 at 10:08
  • $\begingroup$ Type edit care if you want to compare to how Matlab solves this. $\endgroup$ – horchler Oct 25 '17 at 13:38
  • 1
    $\begingroup$ You can always check if your X is correct by plugging it into the CARE you are trying to solve. Also normally Q and R are symmetric, in which case X should be symmetric as well. $\endgroup$ – Kwin van der Veen Oct 25 '17 at 14:20

I wanted to add a little bit of information from my own experience using this method in Matlab for others who find this. The problem is that Matlab's schur command doesn't order the eigenvalues properly for the method you are describing. Take a look at the diagonal of the S matrix you are outputting and make sure the eigenvalues that correspond to the stable (negative) eigenvalues are in the upper-most diagonal elements of S. If they are not, you can use Matlab's ordschur command to put them there.

Looking at your problem, the stable eigenvalues appear in the middle two elements of the diagonal in the S matrix. This means the eigenvectors in U that you are using are incorrect.

See the ordshur documentation for more info. You can use the 'lhp' (left-hand plane) keyword to place the generalized eigenvectors corresponding to the stable eigenvalues in the first columns of U and then your method should work. Note that this method won't work if your matrices are singular or if R violates some conditions. You'll need to use the extended Hamiltonian pencil to solve if that is the case.

Alternatively, you can also play with the 'udi' keyword for ordering the eigenvalues which lie inside the unit circle for the discrete Hamiltonian method and the extended Hamiltonian pencil corresponding to the DARE.

H = [A, -B*inv(R)*B.'; -Q, -A.'];

[U, S] = schur(H);

[U, S] = ordschur(U, S, 'lhp');

U = mat2cell(U, size(H)/2, size(H)/2);
X = U{2, 1}/U{1, 1};

Here is the answer. It will give the same results as MATLAB's command care(A, B, Q, R). But MATLAB's command dare(A, B, Q, R) will not give the same results if we follow the PDF source in the question.

A = [0 1; -3 -4];
B = [1; 1];
C = eye(2);
D = [0; 0];
R = 4;
Q = [1 1; 3 0];

H = [A -(B*inv(R)*B'); -Q -A'];

[U, S] = schur(H);

[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)

>> X = 
      -6.4282  -1.6249
      -7.1621  -2.2296

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.