ParametricNDSolve and Constraint satisfaction problem

Question

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!

Answer 1

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} *)