[Math] Solve Algebraic Riccati Equation via MATLAB

control theorylinear-controlMATLABoptimal control

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

Best Answer

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};
Related Question