# ParametricNDSolve and Constraint satisfaction problem

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!

-
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

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

-
Sorry to bump this old thread, but some time ago when I was no longer working with this I thought about a different approach that possibly could give what I wanted. I was thinking that maybe using CountourPlot[] possibly would be more efficient, with the countours specifying which areas would satisfy the constraints? Anyways, this is just an idea. If I get some spare time I'll try it out just for fun – Sosi Oct 2 '14 at 13:14