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

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

Quote:

Originally Posted by

**rem88** 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

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

Quote:

Originally Posted by

**CaptainBlack** 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

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

Quote:

Originally Posted by

**rem88** 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

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

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.

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