0
$\begingroup$

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};
NMinimize[{fcost, Thread[paralist > 0]}, paralist]

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.

$\endgroup$
  • 1
    $\begingroup$ Consider using ParametricNDSolveValue[] along with the usual fcost[parameters_ /; VectorQ[parameters, NumericQ]] := (* stuff *). $\endgroup$ – J. M. will be back soon Jun 15 '15 at 19:29
  • $\begingroup$ @Guess who it is Thank you very much for your earnest edit and valuable response to my question. I will try it right now. $\endgroup$ – Zihu Guo Jun 15 '15 at 23:56
  • $\begingroup$ 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 $\endgroup$ – Zihu Guo Jun 16 '15 at 2:35
  • $\begingroup$ Can you edit your question to show the changes you did? $\endgroup$ – J. M. will be back soon Jun 16 '15 at 2:48
  • $\begingroup$ @Guesswhoitis. ok, I have edit my question. Thank you very much. $\endgroup$ – Zihu Guo Jun 16 '15 at 2:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.