# How to optimize a cost function which involve the solution of a ODE system using the NMinmize?

I have created a cost function, the value of which depends on the solution of a ODE (ordinary differential equation). The parameters of the ODE are not determined, and my goal is to determine them which can minimize the cost function.

This is the ODE:

odeeqns[parameter_] :=
Module[{sa, sb, sc, sd, db, dc, dd, kab, kbc, kcd, kde, kef, kdg,
Jab, Jbc, Jcd, equations},
{sa, sb, sc, sd, db, dc, dd, kab, kbc, kcd, kde, kef, Jab, Jbc,
Jcd} = parameter;
equations = {molA[t] == sa,
molB'[t] ==
sb - db molB[t] + (kab molA[t] molB[t])/(molB[t] + Jab),
molC'[t] ==
sc - dc molC[t] + (kbc molB[t] molC[t])/(molC[t] + Jbc),
molD'[t] ==
sd - dd molD[t] + (kcd molC[t] molD[t])/(molD[t] + Jcd),
molB[0] == 0, molC[0] == 0, molD[0] == 0};
equations
]


and this is the cost function I have made, it depend on the solution at $t=28$, and fixed parameters $sa$. And the argument of the cost function is a list of numbers.

fcost[parameters_List] :=
Module[{sa0, sa1, tt, parameter0, parammeter1, equa0, equa1, sol0,
sol1, molC0, molB1, molC1, cc},
tt = 28;
sa0 = 0;
sa1 = 1;
parameter0 = Join[{sa0}, parameters];
parammeter1 = Join[{sa1}, parameters];
equa0 = odeeqns[parameter0];
equa1 = odeeqns[parammeter1];
sol0 = NDSolve[equa0, {molC, molD}, {t, 0, 200}] // First;
sol1 = NDSolve[equa1, {molC, molD}, {t, 0, 200}] // First;
{molC0, molD0} = {molC, molD} /. sol0;
{molC1, molD1} = {molC, molD} /. sol1;
cc = (molC0[tt]/molC1[tt] - 12)^2 + (molD0[tt]/molD1[tt] - 3)^2;
cc
]


Given a set of parameters, I can get the value of the cost function; for example

paralist = {0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1,
0.1, 0.1, 0.3, 0.1};
fcost[paralist]


I get the value

131.977

So I tried to solve the optimize problem using the fucntion NMinimize, but I do not know how to use it correctly. The examples on the documentation center are really simple.

I tried like this

paralist = {sa, sb, sc, sd, db, dc, dd, kab, kbc, kcd, kde, kef, Jab, Jbc, Jcd};


But I got an error. So can I solve this problem with this function NMinimize?

I changed the fcost function like this as @Guesswhoitis's comment:

fcost2[parameterlist_ /; VectorQ[parameterlist, NumericQ]] :=
Module[{sa0, sa1, tt, parameter0, parammeter1, equa0, equa1, sol0,
sol1, molC0, molB1, molC1, cc, myfunD0, myfunC0, myfunC1, myfunD1,
paralist, sa, sb, sc, sd, db, dc, dd, kab, kbc, kcd, kde, kef, Jab,
Jbc, Jcd, myfunD, myfunC},
tt = 28;
sa0 = 0;
sa1 = 1;
paralist = {sa, sb, sc, sd, db, dc, dd, kab, kbc, kcd, kde, kef,
Jab, Jbc, Jcd};
myfunD = ParametricNDSolveValue[odeeqns[paralist],
molD, {t, 0, 100}, paralist];
myfunC =
ParametricNDSolveValue[odeeqns[paralist], molC, {t, 0, 100},
paralist];
parameter0 = Join[{sa0}, parameterlist];
parammeter1 = Join[{sa1}, parameterlist];
myfunD0 = myfunD @@ parameter0;
myfunD1 = myfunD @@ parammeter1;
myfunC0 = myfunC @@ parameter0;
myfunC1 = myfunC @@ parammeter1;
cc = (myfunC0[tt]/myfunC1[tt] - 12)^2 + (myfunD0[tt]/myfunD1[tt] -
3)^2;
cc
]


and

parameterlist = {0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1,
0.1, 0.1, 0.3, 0.1};
fcost2[parameterlist]


I got the same result

131.977

but I still have no idea about how to optimize it with NMinimize correctly.

I tried like this

variableist = {sb, sc, sd, db, dc, dd, kab, kbc, kcd, kde, kef, Jab,
Jbc, Jcd};
NMinimize[
fcost2[variableist], {{sb, 0.1}, {sc, 0.1}, {sd, 0.1}, {db,
0.1}, {dc, 0.1}, {dd, 0.1}, {kab, 0.1}, {kbc, 0.1}, {kcd,
0.1}, {kde, 0.1}, {kef, 0.1}, {Jab, 0.1}, {Jbc, 0.3}, {Jcd,
0.1}}]


This is the output

NMinimize::lvar: Variables {{sb,0.1},{sc,0.1},{sd,0.1},{db,0.1},{dc,0.1},{dd,0.1},{kab,0.1},{kbc,0.1},{kcd,0.1},{kde,0.1},{kef,0.1},{Jab,0.1},{Jbc,0.3},{Jcd,0.1}} should be a list of variables, with each element being a variable, or a list containing a variable and lower and upper bounds for the starting region for that variable. >>

NMinimize[ fcost2[{sb, sc, sd, db, dc, dd, kab, kbc, kcd, kde, kef, Jab, Jbc, Jcd}], {{sb, 0.1}, {sc, 0.1}, {sd, 0.1}, {db, 0.1}, {dc, 0.1}, {dd, 0.1}, {kab, 0.1}, {kbc, 0.1}, {kcd, 0.1}, {kde, 0.1}, {kef, 0.1}, {Jab, 0.1}, {Jbc, 0.3}, {Jcd, 0.1}}]

I really do not know how to use this function when the parameter of the object function is a List.

• Consider using ParametricNDSolveValue[] along with the usual fcost[parameters_ /; VectorQ[parameters, NumericQ]] := (* stuff *). Commented Jun 15, 2015 at 19:29
• @Guess who it is Thank you very much for your earnest edit and valuable response to my question. I will try it right now. Commented Jun 15, 2015 at 23:56
• Hi, @Guesswhoitis. I have changed my fcost function like this fcost[parameters_ /; VectorQ[parameters, NumericQ]] := (* stuff *), which involved the funciton ParametricNDSolveValue[]. But I still have no idea to optimize it with the function NMinimize Commented Jun 16, 2015 at 2:35
• Can you edit your question to show the changes you did? Commented Jun 16, 2015 at 2:48
• @Guesswhoitis. ok, I have edit my question. Thank you very much. Commented Jun 16, 2015 at 2:50