# Matrix riccati differential equation using matlab

Tags:
1. Aug 5, 2016

### Jeffrey Eiyike

1. The problem statement, all variables and given/known data

2. Relevant equations

3. The attempt at a solution

2. Aug 5, 2016

### Jeffrey Eiyike

I can upload the code if you want to see my attempt.

3. Aug 5, 2016

### Staff: Mentor

Yes, please do. It's best if you post it inline using code tags. It should look like this:
[code]
[/code]

4. Aug 5, 2016

### Jeffrey Eiyike

Code (Text):
function dPdt = mRiccati(t,P, B, Q, R)
%computes Phi
P = reshape(P, size(Q)); % Convert from "n^2"-by-1 to "n"-by-"n"
dPdt = -Q + P' * B * R^-1 * B.' * P; % Determine derivative
dPdt = dPdt(:); % Convert from "n"-by-"n" to "n^2"-by-1 or converting to column vector
end

function dhdt = mRiccati2(t, P0, B, h, R)
%compute h
h = reshape(h, size(B)); % Convert from "n^2"-by-1 to "n"-by-"n"
dhdt = 0.5* P0' * B * R^-1 * B.' * h; % Determine derivative
dhdt = dhdt(:); % Convert from "n"-by-"n" to "n^2"-by-1
end

function dXdt = mRiccati3(t,X, h0, P0, B, R, sigma)
% Compute Chi
X = reshape(X, size(R)); % Convert from "n^2"-by-1 to "n"-by-"n"
dXdt = (1/2)* h0' * B * R^-1 * B.' * h0- ((sigma^2)/2)*trace(-P0) ; % Determine derivative
dXdt = dXdt(:); % Convert from "n"-by-"n" to "n^2"-by-1
end

clear all; close all
% sigma=[1 1 1 0.5 0.5 0.5 0.33];
sigma=[1 1 1 1.4142 1.4142 1.4142 1.7321];
q=[1 1 1 1 1 1 1];
C=[5; 5; 5; 10; 10; 10; 15];
e=[4; 3; 2; 2; 2; 2; 3];

dt=0.1;        %step size
B = [1; 1];    % Input vector Matrix
R = 1;          %penalty on the optimal  U that updates the contract
N = [0 3]; % Final condition Linear term
t_end = 20;
tspan = t_end:(-0.1):0;
%
%Solves the Matrix Differential equation Backwards

for ii=1:length(sigma);%:length(sigma)

Q = [q(ii) 0; 0 0]; %Final condition on the quadratic term

%Solution for the first Matrix Differential equation
[T1 P] = ode45(@(t,P)mRiccati(t,P, B, Q, R), tspan, Q);

for i= 1:length(P)
P0 = reshape(P(i,:), size(Q));

%Solution for the second Matrix Differential equation (Phi)
[T2, h] = ode45(@(t,h) mRiccati2(t, P0, B, h, R), tspan, N);

% obtain matrix P from last entry (initial value)(h)

h0 = reshape(h(i,:), size(B));

%Solution for the Third Matrix Differential equation(X)
[T3, X] = ode45(@(t,X) mRiccati3(t,X, h0, P0, B, R, sigma(ii)), tspan, 0);
end

% Flips the parameter
h=flipud(h);
P=flipud(P);
X=flipud(X);

% Compute the value function, error and contract integrand
for j= 1: length(P)
%         display(num2str(length(P)))
ZZ= [e(ii, j)  C(ii, j)]; % State
V(ii, j)= -0.5* ZZ*reshape(P(j,:), size(Q))*ZZ' + reshape(h(j,:), size(B))'*ZZ'- X(j);  % Value function
U(ii, j)= inv(R)*B'*(-reshape(P(j,:), size(Q))*ZZ'+ reshape(h(j,:), size(B)));
C(ii, j+1) =  C(ii, j) + U(ii, j)*dt;    %Contract
e(ii, j+1)=e(ii, j)+ U(ii, j)*dt ;%-sigma(ii)*sqrt(dt)*normrnd(0,1);   %error
end

Last edited by a moderator: Aug 10, 2016
5. Aug 7, 2016

### Staff: Mentor

So does your code work or are you having problems with it? If there are problems or errors, tell us about them. Don't make us guess at what might be going wrong.

6. Aug 7, 2016

### Jeffrey Eiyike

My error e is surpose to converge to zero an remain there and the contract C is suppose to converge to 1 and remain there. But mine converges and later starts increasing

7. Aug 7, 2016

### Staff: Mentor

Your three functions don't seem to me to be working correctly.
Code (Matlab M):
function dPdt = mRiccati(t,P, B, Q, R)
%computes Phi
P = reshape(P, size(Q)); % Convert from "n^2"-by-1 to "n"-by-"n"
...

If Q is an n2 by 1 vector, then size(Q) evaluates to [n2, 1], not [n, n]. See size in the matlab documentation. The call to the reshape function doesn't actually change P, as far as I can see. The same problem exists for your other three functions. There might be other problems, but this is the first one I noticed.

Where is the code for the mRiccati function? It doesn't appear to be something provided by MathWorks, but I could be wrong.