How to solve a Riccati Control (differential) Equation?

Hi everyone! First time writing in this forum :)

My problem is that I have to solve the Riccati Control Equation, of the form:

dS/dt = -A' S -S A - Q + S G S

where A' is the transpose. A, S, Q and G are real matrices in R^(N x N). I need the steady state solution of this equation (within a real time simulink scheme, but this for now is not necessary).

Should I use a ODE solver? But how i define the function with matrices? and what about the timespan?

Thank you for your help!

Answer 1

[Giovanni Mottola]

You don't need to solve a differential equation. "Steady-state" means the solution (here, the matrix S) remains constant, which means that the derivative is zero. Substituting dS/dt=0 in your equation we get:

0= -A' S -S A - Q + S G S

which is an algebraic problem (no derivatives involved) in N x N equations and N x N unknowns.

This can be solved directly. You may define a function func.m as follows:

function diff = func( A, Q, G, x )
diff = A.'*x+x*a+Q-x*G*x;
end

Then from the command window define your initial guess x0 and call:

fsolve(@(x) func(A, Q, G, x), x0)

The alternative would be to use the "care" command ( http://it.mathworks.com/help/control/ref/care.html ), but you first have to find the vector B such that B*B.'=G.

[Davide Fabbroni]

Thank you very much for your answer!

However, I know that the steady state solution is the solution of the algebraic equation, but I needed the differential equation in order to keep track of that solution. My problem was that the care problem for my system and weights resulted in a Hamiltonian with eigenvalues on the imaginary axis.

While the differential equation can be always solved, even if there is no steady state solution, being that solution divergent.

I found this code that reshape the matrices in vectors so that I can use the ODE solver. I only had to modify it a little.

function dXdt = mRiccati(t, X, A, B, Q, R) 
X = reshape(X, size(A)); 

%Convert from "n^2"-by-1 to "n"-by-"n" 
dXdt = A.'*X + X*A - X*B*inv(R)*B.'*X + Q; 

%Determine derivative 
dXdt = dXdt(:); %Convert from "n"-by-"n" to "n^2"-by-1

Then, in the main code:

% Differential
X0 = ones((nx+nu)*(nx+nu),1); 
% You can use the following command to solve the system of 
% differential equations:
[T X] = ode45(@mRiccati, [0 100], X0, [], Ar, Br, Qr, P); 
% ODE45 returns "X" as a vector at each time step. 

% You may use the following code to reshape each row of "X" to 
% get the matrix and store it in a cell array:
[m n] = size(X); 
XX = mat2cell(X, ones(m,1), n); 
fh_reshape = @(x)reshape(x,size(Ar)); 
XX = cellfun(fh_reshape,XX,'UniformOutput',false);
XX_last = cell2mat(XX(end))

This code is not mine, but I can't find the answer on this forum anymore.

[jalal khodaparast]

I want to solve a complex differential Riccati equation numerically. I used ODE45 and ODE15s but The solution has unstable oscillation.

How can I solve a differential Riccati equation with complex values?

Answer 2

[LAKHLIFA SADEK]

<https://www.mathworks.com/matlabcentral/fileexchange/72853-symmetric-differential-matrix-riccati-equations?s_tid=prof_contriblnk>