Possible to solve Algebraic Riccati Equation through ODE45?

Question

The Algebraic Riccati Equation is:

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

And the Algebraic Riccati Difference Equation is:

$$A'X + XA - XBR^{-1}B'X + Q = \frac{dX}{dt}$$

Some times $\frac{dX}{dt}$ is just called $S$ in books.

But is it not possible to solve the Algebraic Riccati Difference Equation through ODE45 and break the simulation when $-tol <= \frac{dX}{dt} <= tol$ ?

You may wonder "Why not use Schur's method to solve the Riccati equations?". Well, Schur's method is a very bad method and there is a reason why MATLAB and Octave are not using it.

Sure, the method find the solution $X$ to the riccati equation, but build a control law $L$ for the LQR controller by using that solution is not what I recommending.

Right now I using fsolve, but that's a built in library in Octave and a function in Optimization package for MATLAB.

Answer 1

Yes it is. Here is the solution:

Step 1: - Create the function. In this case it's CARE

function [value, isterminal, direction] = func(t, X, A, B, Q, R)

 X = reshape(X, size(A)); % Vector to matrix
 value = A'*X + X*A - X*B*inv(R)*B'*X + Q; % Values is the derivative of X
 value = value(:); % Turn value to a vector

 isterminal = 1;   % Stop the integration - An event!  
 direction  = 0;

Step 2: - Create your matrices

A =

     1     1
     2     1

>> B

B =

     1
     1

>> Q

Q =

     1     0
     0    1

>> R

R =

     5

Step 3: - Simulate

 >> [T X] = ode45(@(t,X)func(t, X, A, B, Q, R), [0 100], X0) ;

Step 4: - Find the last values for X and generate a matrix x

 >> x = reshape(X(size(X,1), :), size(A))

 x =

    8.7541    5.4230
    5.4230    5.0536

Step 5: - Check if the solution is almost zero

 >> A'*x + x*A - x*B*inv(R)*B'*x + Q

 ans =

     0.0022    0.0017
     0.0017    0.0013

Yes it is.

Here is the plot how we can se that after a few seconds, the derivative of $X$ is almost zero.

I use this command

 plot(T, X)

....This is a knif....solution :)

Answer 2

You may wonder "Why not use Schur's method to solve the Riccati equations?". Well, Schur's method is a very bad method and there is a reason why MATLAB and Octave are not using it.

Nope, that's precisely what matlab is using. It first checks whether $R$ is well-conditioned and if so forms the Hamiltonian and uses the QZ decomposition. If $R$ is not well-conditioned then it forms the Extended Hamiltonian and applies the Van Dooren's deflation.

You can follow the logic from SciPy's implementation of solve_continuous_are which always uses the extended Hamiltonian instead of an $R$ dependent solution.