1. Limited time only! Sign up for a free 30min personal tutor trial with Chegg Tutors
    Dismiss Notice
Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Homework Help: Matrix riccati differential equation using matlab

  1. Aug 5, 2016 #1
    aassignment.png 1. The problem statement, all variables and given/known data

    2. Relevant equations

    3. The attempt at a solution
  2. jcsd
  3. Aug 5, 2016 #2
    I can upload the code if you want to see my attempt.
  4. Aug 5, 2016 #3


    Staff: Mentor

    Yes, please do. It's best if you post it inline using code tags. It should look like this:
    << your matlab code>>
  5. Aug 5, 2016 #4
    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

    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

    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

    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);
    % Flips the parameter

    % 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    
    Last edited by a moderator: Aug 10, 2016
  6. Aug 7, 2016 #5


    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.
  7. Aug 7, 2016 #6
    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
  8. Aug 7, 2016 #7


    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.
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook

Have something to add?
Draft saved Draft deleted