Thread: 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

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

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

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

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

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.

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