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

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

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

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

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