Results 1 to 7 of 7

Math Help - Matlab - Running fmincon with ODE45 in seperate M-file

  1. #1
    Junior Member
    Joined
    Apr 2009
    Posts
    45

    Matlab - Running fmincon with ODE45 in seperate M-file

    Hi,

    I've been having some problems writing this program and i'm not sure where I am going wrong. It seems as though the variabe I want to change (a) does not.

    The problem is 4 sets of differentials, and the aim is to minimise the value a (which is a feedrate) to minimise a cost function (using fmincon)

    The first m-file is a simulation (which contain the differentials):

    The equations i am trying to represent are:

    dE/dt = -b1*E + F
    dG/dt = r1 - (1/b1)*r2
    dP/dt = r2
    dZ/dt = F

    F = a

    My F (flow rate) for now is simply F = a, i plan to make this more complicated alter

    ------------------------------------------------------------------------
    function dxdt=Feed(t,x);

    global Feed dxdt a b fd cost_function feedrate E G P Z F r2 x1 x2 x3 x4;

    dxdt = zeros(size(x));



    b = [0.019 0.95 50 43 6 35 0.05 0.04 30 1 178 100 5];


    E = x(1);
    G = x(2);
    P = x(3);
    Z = x(4);
    F = a


    r1 = b(10)*b(3)*x(1)/(b(4)*(1 + x(2)/b(5) + x(3)/b(6)) + b(3));
    r2 = b(8)*b(9)*x(2)/(b(7) + x(2));


    dxdt(1) = -b(1)*E + F; %dE/dt
    dxdt(2) = r1 - (1/b(2))*r2; %dG/dt
    dxdt(3) = b(8)*b(9)*G/(b(7) + G); %dP/dt
    dxdt(4) = F; %dz/dt


    if x(2)<=0
    dxdt(3)=0;
    dxdt(2)=r1;

    x(1) = dxdt(1)
    x(2) = dxdt(2)
    x(3) = dxdt(3)
    x(4) = dxdt(4)

    end
    ------------------------------------------------------------------------

    The cost fucntion file is then, which also solves the ODE:

    -------------------------------------------------------------------------------
    function J = myfun(a)



    global Feed dxdt a b F r2 x1 x2 x3 x4;


    [t,x] = ode45('Feed',[0 200],[16.35 0 0 0]);



    J = (dxdt(3) - b(11)*exp(-dxdt(3)/b(12)) - b(13)*dxdt(4));




    plot(t,[x(:,1) x(:,2) x(:,3)/10 x(:,4)]);

    end
    --------------------------------------------------------------------------

    Ths file to run all this is then:
    ---------------------------------------------------------------------------

    clear all
    clc

    global Feed dxdt E F x1 x2 x3 x4;


    a0 = 1;
    lb = 0;
    ub = 100;
    [a,fval] = fmincon(@myfun,a0,lb,ub)
    --------------------------------------------------------------------------


    What my aim is, that the run file runs the ODE, then calculates J, then changes the value of a, re runs the ODE and keeps doing this to find a value of a to minimise J (in myfun).

    But what ever value starting a value I chose this does not change. (I believe it is because the dxdt(4) does not change).

    Can anyone see anything obviously wrong with it?

    Thanks,

    Rem
    Follow Math Help Forum on Facebook and Google+

  2. #2
    Grand Panjandrum
    Joined
    Nov 2005
    From
    someplace
    Posts
    14,972
    Thanks
    4

    Re: Matlab - Running fmincon with ODE45 in seperate M-file

    Quote Originally Posted by rem88 View Post
    Hi,

    I've been having some problems writing this program and i'm not sure where I am going wrong. It seems as though the variabe I want to change (a) does not.

    The problem is 4 sets of differentials, and the aim is to minimise the value a (which is a feedrate) to minimise a cost function (using fmincon)

    The first m-file is a simulation (which contain the differentials):

    The equations i am trying to represent are:

    dE/dt = -b1*E + F
    dG/dt = r1 - (1/b1)*r2
    dP/dt = r2
    dZ/dt = F

    F = a

    My F (flow rate) for now is simply F = a, i plan to make this more complicated alter

    ------------------------------------------------------------------------
    function dxdt=Feed(t,x);

    global Feed dxdt a b fd cost_function feedrate E G P Z F r2 x1 x2 x3 x4;

    dxdt = zeros(size(x));



    b = [0.019 0.95 50 43 6 35 0.05 0.04 30 1 178 100 5];


    E = x(1);
    G = x(2);
    P = x(3);
    Z = x(4);
    F = a


    r1 = b(10)*b(3)*x(1)/(b(4)*(1 + x(2)/b(5) + x(3)/b(6)) + b(3));
    r2 = b(8)*b(9)*x(2)/(b(7) + x(2));


    dxdt(1) = -b(1)*E + F; %dE/dt
    dxdt(2) = r1 - (1/b(2))*r2; %dG/dt
    dxdt(3) = b(8)*b(9)*G/(b(7) + G); %dP/dt
    dxdt(4) = F; %dz/dt


    if x(2)<=0
    dxdt(3)=0;
    dxdt(2)=r1;

    x(1) = dxdt(1)
    x(2) = dxdt(2)
    x(3) = dxdt(3)
    x(4) = dxdt(4)

    end
    ------------------------------------------------------------------------

    The cost fucntion file is then, which also solves the ODE:

    -------------------------------------------------------------------------------
    function J = myfun(a)



    global Feed dxdt a b F r2 x1 x2 x3 x4;


    [t,x] = ode45('Feed',[0 200],[16.35 0 0 0]);



    J = (dxdt(3) - b(11)*exp(-dxdt(3)/b(12)) - b(13)*dxdt(4));




    plot(t,[x(:,1) x(:,2) x(:,3)/10 x(:,4)]);

    end
    --------------------------------------------------------------------------

    Ths file to run all this is then:
    ---------------------------------------------------------------------------

    clear all
    clc

    global Feed dxdt E F x1 x2 x3 x4;


    a0 = 1;
    lb = 0;
    ub = 100;
    [a,fval] = fmincon(@myfun,a0,lb,ub)
    --------------------------------------------------------------------------


    What my aim is, that the run file runs the ODE, then calculates J, then changes the value of a, re runs the ODE and keeps doing this to find a value of a to minimise J (in myfun).

    But what ever value starting a value I chose this does not change. (I believe it is because the dxdt(4) does not change).

    Can anyone see anything obviously wrong with it?

    Thanks,

    Rem
    What do you think :

    Code:
    x(1) = dxdt(1)
    x(2) = dxdt(2)
    x(3) = dxdt(3)
    x(4) = dxdt(4)
    does in the derivative function?

    CB
    Follow Math Help Forum on Facebook and Google+

  3. #3
    Junior Member
    Joined
    Apr 2009
    Posts
    45

    Re: Matlab - Running fmincon with ODE45 in seperate M-file

    Quote Originally Posted by CaptainBlack View Post
    What do you think :

    Code:
    x(1) = dxdt(1)
    x(2) = dxdt(2)
    x(3) = dxdt(3)
    x(4) = dxdt(4)
    does in the derivative function?

    CB
    I thought it replaced the values of x(1) in dxdt(1) etc for the next calculation?

    I put this in because dxdt(4) was not changing and I wasnt sure why.

    dxdt(4) should increase with time becuase it is the total amount of F (=a) fed, but it stays as dxdt(4) = a
    Follow Math Help Forum on Facebook and Google+

  4. #4
    Grand Panjandrum
    Joined
    Nov 2005
    From
    someplace
    Posts
    14,972
    Thanks
    4

    Re: Matlab - Running fmincon with ODE45 in seperate M-file

    Quote Originally Posted by rem88 View Post
    I thought it replaced the values of x(1) in dxdt(1) etc for the next calculation?
    You should not do this, ode45 is going to update x using dxdt. The function has one purpose; to calculate the derivatives at (t,x), it should do nothing more. Also it is bad practice to modify the arguments in a function.

    CB
    Follow Math Help Forum on Facebook and Google+

  5. #5
    Junior Member
    Joined
    Apr 2009
    Posts
    45

    Re: Matlab - Running fmincon with ODE45 in seperate M-file

    OK, thanks. I wasn't sure if they were updated. It must be the way I have written it that dxdt(4) doesn't increase
    Follow Math Help Forum on Facebook and Google+

  6. #6
    Junior Member
    Joined
    Apr 2009
    Posts
    45

    Re: Matlab - Running fmincon with ODE45 in seperate M-file

    I was wrong dxdt(4) does increase, I just was not seeing it. When I plot a graph I can see that it increases.

    But I dont see why the value for a does not change.
    Follow Math Help Forum on Facebook and Google+

  7. #7
    Junior Member
    Joined
    Apr 2009
    Posts
    45

    Re: Matlab - Running fmincon with ODE45 in seperate M-file

    I have figure out that all the derivatives update, but the a value does not change.

    I think I know the problem, because I am using global a, a is being replaced each time and a0 in fminunc is not used.

    But if I remove a, then the program doesnt work.

    I thought setting the inital condition in fminunc would update the a value, as a is a function of myfun
    Follow Math Help Forum on Facebook and Google+

Similar Math Help Forum Discussions

  1. Matlab ode45
    Posted in the Math Software Forum
    Replies: 8
    Last Post: October 24th 2009, 12:22 PM
  2. Matlab ode45
    Posted in the Differential Equations Forum
    Replies: 0
    Last Post: October 24th 2009, 06:23 AM
  3. Solving IVP in Matlab (ode45)
    Posted in the Math Software Forum
    Replies: 2
    Last Post: December 15th 2008, 10:01 PM
  4. ode45 Matlab
    Posted in the Math Software Forum
    Replies: 6
    Last Post: November 29th 2008, 09:21 PM
  5. ode45 MATLAB
    Posted in the Math Software Forum
    Replies: 2
    Last Post: September 24th 2008, 05:02 AM

Search Tags


/mathhelpforum @mathhelpforum