# How do I solve the Riccati Differential Equation with fsolve in MATLAB / Octave?

Let's say that I have this equation. This is the Riccati Differential Equation:

$$P A + A^T P - P B R^{-1} B^T P + Q = 0$$

I know $A, B, Q, R$. My goal is to find $P$. How do I use the MATLAB / Octave command fsolve to find $P$ ?

EDIT: 

Here is the answer - An example to get LQR matrix.

>> fun = @(P, A, B, Q, R) P*A+A'*P-P*B*inv(R)*B'*P + Q
fun =

@(P, A, B, Q, R) P * A + A' * P - P * B * inv (R) * B' * P + Q
>> p0 = 0.1*ones(2)
p0 =

0.10000   0.10000
0.10000   0.10000

>> P = fsolve(@(P) fun(P, Amat, Bmat, Q, R), p0)
P =

55.67637    0.57241
0.57241    2.29935

>>

>> inv(R)*Bmat'*P
ans =

2.2269e-02   2.2895e-04
-1.4350e-03  -5.7645e-03

>> lqr(Amat, Bmat, Q, R)
ans =

2.2269e-02   2.2899e-04
-1.4353e-03  -5.7644e-03

>>

• fminsearch could also be possible to use. I have more experience using it anyway. A third approach could be to do a linearization and minimize a regularized least squares norm of vectorization of P using Kronecker products for matrix multiplication. – mathreadler Jul 2 '17 at 10:16
• Thank you for your reply. I don't think that fminsearch is the right solution. I'm using GNU Octave and Octave have fsolve as a built in function. The same as MATLAB has. – Daniel Mårtensson Jul 2 '17 at 10:18
• What do you mean by "right" solution? Is it part of a course where you have to use a particular function to pass? If so then you are supposed to show any own attempts at the problem in the question. – mathreadler Jul 2 '17 at 10:21
• So if I use the fminsearch, I can solve the Riccati equation? – Daniel Mårtensson Jul 2 '17 at 10:34
• Yes it should work. For example if you aim to minimize $\|PA+A^TP -PBR^{-1}B^TP+Q\|_2$. But whichever function you use you will need to get comfortable using the programming language so you know what you are doing. – mathreadler Jul 2 '17 at 11:22

Define a function, demofun:

function err = demofun(p, A, B, Q, R)
P = reshape(p, [3 3]); %% Reshape p to your desired matrix size
F = P*A+A'*P-P*B*inv(R)*B'*P + Q;

% Calculate the error value
err = sum(F(:).*F(:));

return


Then write the script that optimizes the function

A = randn(3);
B = randn(3);
Q = randn(3);
R = randn(3);

p0 = ones(9, 1)*0.1; % Initial guess for P, vectorized to p

tmp = fsolve(@(p) demofun(P, A, B, Q, R), p0);

P = reshape(tmp, [3 3]); % Reshape vectorized answer to matrix format

• No need to use reshape. Let's say that I got those pastebin.com/b9q9rcv1 This will give me L = [2.2269e-02 2.2899e-04; -1.4353e-03 -5.7644e-03] if I use the LQR command in MATLAB. – Daniel Mårtensson Jul 2 '17 at 20:32