Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I am trying to find the values of two parameters that allow for a specific result of my differential equations, under different initial conditions.

Given a system of equations:

TestSystem = {H'[t] == p - H[t] (k1 sh[t]), H[0] == H0,
   sh'[t] == -k1 H[t] sh[t] + k2 ss[t], sh[0] == sh0,

   st == sh[t] + ss[t]};

I use ParametricNDSolve to compute the values of the dependent variables, while leaving k1 and k2 unspecified:

TestSolver[H0e_, st0_] :=
 ParametricNDSolve[
  TestSystem /. p -> 0 /. {H0 -> H0e}
      /. sh0 -> st0 /. st -> st0
    /. k2 -> k2estim /. k1 -> k1estim, {h, sh, st}, {t, 0, 
   300}, {k1estim, k2estim}];

The conditions that I want to evaluate are then specified by m,d,hh, at t=300:

m = 100 sh[k1estim, k2estim][300]/st;
d = 100 ss[k1estim, k2estim][300]/st;
hh = h[k1estim, k2estim][300];

cond1a[k1estim_, k2estim_] := 
 Evaluate[m /. st -> 5.*^-6 /. TestSolver[1.*^-5, 5.*^-6]]
cond1b[k1estim_, k2estim_] := 
 Evaluate[d /. st -> 5.*^-6 /. TestSolver[1.*^-5, 5.*^-6]]
cond1c[k1estim_, k2estim_] := 
 Evaluate[hh /. TestSolver[1.*^-5, 5.*^-6]]

cond2a[k1estim_, k2estim_] := 
 Evaluate[m /. st -> 5.*^-6 /. TestSolver[5.*^-5, 5.*^-6]]
cond2b[k1estim_, k2estim_] := 
 Evaluate[d /. st -> 5.*^-6 /. TestSolver[5.*^-5, 5.*^-6]]
cond2c[k1estim_, k2estim_] := 
 Evaluate[hh /. TestSolver[5.*^-5, 5.*^-6]]

Now, say that I wanted to find the following

1 < cond1a < 100,
0 < cond1b < 1,
1.*^-8 < cond1c < 1*^-5,
1 < cond2a < 50,
50 < cond2b < 100,
1.*^-8 < cond2c < 1*^-5,

How could I find the values of k1 and k2 that would satisfy such a thing?

I was considering using NMinimize like this:

NMinimize[{1/cond1a+cond1b+cond1c+cond2a+1/cond2b+cond2c,cons},{k1,k2}]

(where cons would be the required constrains).

Yet, this does not seem like an efficient and plausible way to solve this problem. Does anyone have any ideas on how I could approach this problem? Given that this is a multi-objective optimization problem, how would you suggest that I could specify the objective function?

Thank you so much!

share|improve this question
1  
It appears to be a constraint satisfaction problem rather than optimization per se. You might thus use a simple objective, say (the constant function) 1, and pass the constraints in a way that they get evaluated for numeric values of the parameters k1,k2. –  Daniel Lichtblau Aug 20 '13 at 21:07
    
@DanielLichtblau I see! You are suggesting considering any objective function of the above (e.g. minimize cond1b, since I want it to be close to 0), and apply all the constraints, correct? –  Sosi Aug 21 '13 at 12:53
    
Yes, that's correct. –  Daniel Lichtblau Aug 21 '13 at 14:46
    
Well, indeed that seems like a nice suggestion! Now my problem is having a very wide span for both parameters. Evaluating everything yields some NMinimize::incst: error. I'll try to work around this. Also, would you like to post your suggestion as answer? –  Sosi Aug 21 '13 at 15:16

1 Answer 1

up vote 2 down vote accepted

I had to make a number of corrections and adjustments to the code, but here is the idea.

TestSystem = {h'[t] == p - h[t] (k1 sh[t]), h[0] == H0, 
   sh'[t] == -k1 h[t] sh[t] + k2 ss[t], sh[0] == sh0, 
   st == sh[t] + ss[t]};

TestSolver[H0e_, st0_] := 
  ParametricNDSolve[
   TestSystem /. p -> 0 /. {H0 -> H0e} /. sh0 -> st0 /. st -> st0 /. 
     k2 -> k2estim /. k1 -> k1estim, {h, sh, ss}, {t, 0, 
    300}, {k1estim, k2estim}];

m[k1estim_?NumberQ, k2estim_?NumberQ] := 
  100 sh[k1estim, k2estim][300]/st;
d[k1estim_?NumberQ, k2estim_?NumberQ] := 
  100 ss[k1estim, k2estim][300]/st;
hh[k1estim_?NumberQ, k2estim_?NumberQ] := h[k1estim, k2estim][300];

tsres1 = TestSolver[1.*^-5, 5.*^-6];
tsres5 = TestSolver[5.*^-5, 5.*^-6];

cond1a[k1estim_?NumberQ, k2estim_?NumberQ] := 
 m[k1estim, k2estim] /. st -> 5.*^-6 /. tsres1
cond1b[k1estim_?NumberQ, k2estim_?NumberQ] := 
 d[k1estim, k2estim] /. st -> 5.*^-6 /. tsres1
cond1c[k1estim_?NumberQ, k2estim_?NumberQ] := 
 hh[k1estim, k2estim] /. tsres1


cond2a[k1estim_?NumberQ, k2estim_?NumberQ] := 
 m[k1estim, k2estim] /. st -> 5.*^-6 /. tsres5
cond2b[k1estim_?NumberQ, k2estim_?NumberQ] := 
 d[k1estim, k2estim] /. st -> 5.*^-6 /. tsres5
cond2c[k1estim_?NumberQ, k2estim_?NumberQ] := 
 hh[k1estim, k2estim] /. tsres5

constraints[k1_, k2_] = {1 <= cond1a[k1, k2] <= 100, 
   0 <= cond1b[k1, k2] <= 1, 1.*^-8 <= cond1c[k1, k2] <= 1*^-5, 
   1 <= cond2a[k1, k2] <= 50, 50 <= cond2b[k1, k2] <= 100, 
   1.*^-8 <= cond2c[k1, k2] <= 1*^-5};

Now can do e.g.

Timing[{min, vals} = 
  NMinimize[{1/cond1a[k1, k2]^2 + cond1b[k1, k2]^2 + 
     cond1c[k1, k2]^2 + 1/cond2a[k1, k2]^2 + 1/cond2b[k1, k2]^2 + 
     cond2c[k1, k2]^2, 
    constraints[k1, k2]}, {{k1, -100, 100}, {k2, -100, 100}}]]

I'm not sure how good a job it will do though. I played a bit more and got a possibly viable result by changing the initial ranges.

Timing[{min, vals} = 
  NMinimize[{1/cond1a[k1, k2]^2 + cond1b[k1, k2]^2 + 
     cond1c[k1, k2]^2 + 1/cond2a[k1, k2]^2 + 1/cond2b[k1, k2]^2 + 
     cond2c[k1, k2]^2, 
    constraints[k1, k2]}, {{k1, -1000, 1000}, {k2, -10, 10}}]]

[some warning messages omitted]

(* {17.440000, {0.235039300883, {k1 -> 5534.72060398, 
   k2 -> 0.0377418151532}}} *)

Now check for constraint satisfaction.

In[314]:= {cond1a[k1, k2], cond1b[k1, k2], cond1c[k1, k2], 
  cond2a[k1, k2], cond2b[k1, k2], cond2c[k1, k2]} /. vals

(* Out[314]= {99.5161215794, 0.483878420575, 
 1.0140995435*10^-8, 49.9885948839, 50.0114051161, 
 5.75765494239*10^-6} *)

In[315]:= constraints[k1, k2] /. vals

(* Out[315]= {True, True, True, True, True, True} *)
share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.