Matrix riccati differential equation using matlab

[Jeffrey Eiyike]

Homework Statement

Homework Equations

The Attempt at a Solution

 

[Jeffrey Eiyike]

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

 

[Mark44]

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

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

 

[Jeffrey Eiyike]

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
endfunction 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

 

[Mark44]

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.

 

[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

 

[Mark44]

Your three functions don't seem to me to be working correctly.

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 n² by 1 vector, then size(Q) evaluates to [n², 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.