# Possible to solve Algebraic Riccati Equation through ODE45?

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.

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

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.

• I have only read the article "schur method for solving algebraic riccati equations" from 1978. That method is not good enough according to me. Sure, it finds the solution. But develop a LQR control law with that method is not what I recommending. – Daniel Mårtensson Nov 21 '17 at 11:45
• @DanielMårtensson That IS what matlab is using inside LQR. – percusse Nov 21 '17 at 11:52
• I have tried both the schur method and solve the solution in Matlab via Control system Toolbox. Did not give the same results. – Daniel Mårtensson Nov 21 '17 at 17:38
• @DanielMårtensson read the source code. It uses gcare.m – percusse Nov 21 '17 at 19:31
• So schur's method is gcare.m ? – Daniel Mårtensson Nov 21 '17 at 21:36