# ParametricNDSolveValue, NMaximize: Dealing with complex results in maximization problem

I analyze a system of differential equations in several variables including mus[t], mui[t], and w[t]. Time runs from t=0 to t=T. The two initial values mus[0] and mui[0] can freely be chosen; their choice fully determines system dynamics. I want to find values for mus[0] and mui[0] that maximize w[T] subject to system dynamics. My strategy is to use the methods ParametricNDSolveValue and NMaximize.

My code

u[x_] = Log[x] - x + 1;
aformula[x_] = (-1 + (1 + 4 x)^0.5)/(2 x);

par = {beta -> 0.3 + gamma, gamma -> 1/7, nu -> 0.67/365,
rho -> 0.05/365, kappa -> 197};
sys = Simplify[{
ns'[t] == -beta  (a[t])^2  ns[t]  ni[t],
ni'[t] == beta  (a[t])^2  ns[t]  ni[t] - gamma  ni[t],
(rho + nu) mus[t] - mus'[t] ==
u[a[t]] + (mui[t] - mus[t]) beta (a[t])^2  ni[t],
(rho + nu) mui[t] - mui'[t] ==
u[a[t]] + (mui[t] - mus[t]) beta (a[t])^2  ns[t] -
gamma (kappa + mui[t]),
a[t] ==
aformula[2 beta (mus[t] - mui[t]) ns[t] ni[t]/(ns[t] + ni[t])],
w'[t] ==
Exp[-(rho + nu) t] (-kappa  gamma  ni[t] + (ni[t] + ns[t]) u[
a[t]]),
h'[t] == Exp[-(rho + nu) t] (-kappa  gamma  ni[t]),
ns[0] == 0.9999223,
ni[0] == 0.0000527,
w[0] == 0,
h[0] == 0
}];

wT = ParametricNDSolveValue[
Flatten[{sys //. par, mus[0] == mus0, mui[0] == mui0}],
w[T], {t, 0, T}, {mus0, mui0, T}]

Plot3D[wT[Mus0, Mui0, 500], {Mus0, -200, -180}, {Mui0, -23000, -18000}]
NMaximize[{wT[Mus0, Mui0, 500], -200 <= Mus0 <= -180, -23000 <=
Mui0 <= -18000}, {Mus0, Mui0}]


The issue The code works for the values given above (Mus0 in the range -200 to -180, Mui0 in the range -23000 to -18000). Here is part of the output: {-46.7287, {Mus0 -> -180., Mui0 -> -20639.9}}. But when I allow for a larger set of values for Mus0 and Mui0, then it doesn't. Error messages include ParametricNDSolveValue::mconly: For the method IDA, only machine real code is available. Unable to continue with complex values or beyond floating-point exceptions. and InterpolatingFunction::dmval: Input value {500} lies outside the range of data in the interpolating function. Extrapolation will be used. and NMaximize::cvdiv: Failed to converge to a solution. The function may be unbounded. I suspect that the source of the problem is that for ''bad'' choices of Mus0 or Mui0, some variables attain complex values and so does w[T]. I know that system dynamics with complex values cannot constitute the solution of my problem. So all I need to do is to find a way to ''discard'' complex system dynamics.

Options To deal with complex system dynamics, I considered several options. (1) Replace the function aformula with a function that is globally real-valued. (Note that aformula is real-valued over the relevant range, which is x>=0, but not for sufficiently negative x.) I considered functions that are weighted averages of Re[aformula] and, e.g., 1, but I didn't find such averages that wouldn't create numerical issues for very large or small x. (2) Interrupt ParametricNDSolveValue when complex values appear. I understand that WhenEvent might allow to do something like this when one uses NDSolve, but I didn't manage. (The relevant event to trigger such an interruption could be when mus[t] < mui[t] or a[t]>1.) (3) Impose in NMaximize that only paths with mus[t] > mui[t] should be considered. I don't know how to impose such a condition. (4) Others.

Questions Which option is preferable? How can I implement it? Thank you.

Post scriptum A proper solution should have the following properties: (a) All function values are real. (b) a[t]>0 (implied by (a)). (c) ns[t]>0 and ni[t]>=0 (implied by ODEs).

• Could you show desirable range for Mus0, Mui0? Commented Apr 11 at 10:39
• Don't know desirable range for Mus0 and Mui0. Commented Apr 11 at 13:06
• Do you mean $-\infty < Mus0< \infty$ and $-\infty < Mui0< \infty$? Commented Apr 11 at 13:16
• Yes, I do not impose any a-priori restrictions on Mus0 or Mui0 except that they are real. Commented Apr 11 at 13:32
• Please, see my answer with solution without restrictions. Commented Apr 11 at 14:13

One possible solution is to define real functions in a form

real[x_] := Sqrt[x^2]; u[x_] = Log[real[x]] - x + 1;
aformula[x_] = (-1 + real[(1 + 4  x)]^0.5)/(2  x);
a[t_] :=
aformula[2  beta  (mus[t] - mui[t])  ns[t]  ni[t]/(ns[t] + ni[t])];


With these definitions we can compute global maximum as follows

par = {beta -> 0.3 + gamma, gamma -> 1/7, nu -> 0.67/365,
rho -> 0.05/365, kappa -> 197};
sys = {ns'[t] == -beta   (a[t])^2   ns[t]   ni[t],
ni'[t] ==
beta   (a[t])^2   ns[t]   ni[t] -
gamma   ni[t], (rho + nu)  mus[t] - mus'[t] ==
u[a[t]] + (mui[t] - mus[t])  beta  (a[t])^2   ni[
t], (rho + nu)  mui[t] - mui'[t] ==
u[a[t]] + (mui[t] - mus[t])  beta  (a[t])^2   ns[t] -
gamma  (kappa + mui[t]),
w'[t] ==
Exp[-(rho + nu)  t]  (-kappa   gamma   ni[t] + (ni[t] + ns[t])  u[
a[t]]), h'[t] ==
Exp[-(rho + nu)  t]  (-kappa   gamma   ni[t]), ns[0] == 0.9999223,
ni[0] == 0.0000527, w[0] == 0, h[0] == 0};

wT = ParametricNDSolveValue[
Flatten[{sys //. par, mus[0] == mus0, mui[0] == mui0}],
w[T], {t, 0, T}, {mus0, mui0, T}]
NMaximize[wT[Mus0, Mui0, 500], {Mus0, Mui0}]

(*Out[]= {161.849, {Mus0 -> -14697.1, Mui0 -> 39820.9}}*)


Update 1. Taken into account constraint a[t]>0,ni[t]>0,ns[t]>0 we can define two functions with and without WhenEvent to compute extremum and solution as follows

wT2[mus0_?NumberQ, mui0_?NumberQ, T_?NumberQ] :=
Module[{a0, ni0 = 0.0000527, ns0 = 0.9999223, gamma = 1/7,
nu = 0.67/365, rho = 0.05/365, kappa = 197, beta, a, ns, ni, mus,
mui, h, w, t, sol, resol, tm}, beta = 0.3 + gamma;
real[x_] := Sqrt[x^2]; u[x_] = Log[real[x]] - x + 1;
aformula[x_] = (-1 + real[(1 + 4  x)]^0.5)/(2  x);
sys = {ns'[t] == -beta  (a)^2  ns[t]  ni[t],
ni'[t] ==
beta  (a)^2  ns[t]  ni[t] - gamma  ni[t], (rho + nu) mus[t] -
mus'[t] ==
u[a] + (mui[t] - mus[t]) beta (a)^2  ni[t], (rho + nu) mui[t] -
mui'[t] ==
u[a] + (mui[t] - mus[t]) beta (a)^2  ns[t] -
gamma (kappa + mui[t]),
w'[t] ==
Exp[-(rho + nu) t] (-kappa  gamma  ni[t] + (ni[t] + ns[t]) u[
a]), h'[t] == Exp[-(rho + nu) t] (-kappa  gamma  ni[t]),
ns[0] == ns0, ni[0] == ni0, w[0] == 0, h[0] == 0,
mus[0] == mus0, mui[0] == mui0} /.
a -> aformula[
2 beta (mus[t] - mui[t]) ns[t] ni[t]/(ns[t] + ni[t])];
a0 = aformula[2 beta (mus0 - mui0) ns0 ni0/(ns0 + ni0)];
eps = 10^-12;
pred[x1_?NumberQ, x2_?NumberQ, x3_?NumberQ] :=
x1 < eps || x2 < eps || x3 < eps; If[a0 > 0, sol =
Quiet@
NDSolveValue[{sys,
WhenEvent[
pred[aformula[
2 beta (mus[t] - mui[t]) ns[t] ni[t]/(ns[t] + ni[t])],
ns[t], ni[t]], "StopIntegration"]}, {ns, ni, mus, mui, h,
w}, {t, 0, T}], Nothing];
tm = If[a0 <= 0, 0, sol[[6]]["Domain"][[1, 2]]];
resol = If[(tm - T)^2 > 10^-2  T, -1.  10^10, sol[[6]][tm]]; resol];

wT1[mus0_?NumberQ, mui0_?NumberQ, T_?NumberQ] :=
Module[{s, a0, ni0 = 0.0000527, ns0 = 0.9999223, gamma = 1/7,
nu = 0.67/365, rho = 0.05/365, kappa = 197, beta, a, ns, ni, mus,
mui, h, w}, beta = 0.3 + gamma; u[x_] = Log[x] - x + 1;
aformula[x_] = (-1 + (1 + 4 x)^0.5)/(2 x);
sys = {ns'[t] == -beta  (a)^2  ns[t]  ni[t],
ni'[t] ==
beta  (a)^2  ns[t]  ni[t] - gamma  ni[t], (rho + nu) mus[t] -
mus'[t] ==
u[a] + (mui[t] - mus[t]) beta (a)^2  ni[t], (rho + nu) mui[t] -
mui'[t] ==
u[a] + (mui[t] - mus[t]) beta (a)^2  ns[t] -
gamma (kappa + mui[t]),
w'[t] ==
Exp[-(rho + nu) t] (-kappa  gamma  ni[t] + (ni[t] + ns[t]) u[
a]), h'[t] == Exp[-(rho + nu) t] (-kappa  gamma  ni[t]),
ns[0] == ns0, ni[0] == ni0, w[0] == 0, h[0] == 0, mus[0] == mus0,
mui[0] == mui0} /.
a -> aformula[
2 beta (mus[t] - mui[t]) ns[t] ni[t]/(ns[t] + ni[t])]; s =
NDSolve[
sys, {ns, ni, mus, mui, h, w}, {t, 0, T}]; {aformula[
2  beta  (mus[t] - mui[t])  ns[t]  ni[t]/(ns[t] + ni[t])], ns[t],
ni[t], mus[t], mui[t], h[t], w[t]} /. s[[1]]]


Using function wT2 we compute local and global maximum as

local =
FindMaximum[wT2[mus0, mui0, 500], {{mus0, -180}, {mui0, -20600}},
MaxIterations -> 1000]

Out[]= {-40.851, {mus0 -> -107.27, mui0 -> -24528.}}

global =
NMaximize[wT2[mus0, mui0, 500],
Element[{mus0, mui0},
Rectangle[{-200, -2.5  10^4}, {-100, -2  10^4}]]]

Out[]= {-40.8532, {mus0 -> -100.014, mui0 -> -24721.9}}


Visualization

smax = wT1[mus0, mui0, 500] /. local[[2]];

var = {a, ns, ni, mus, mui, h, w}; Table[
Plot[smax[[i]], {t, 0, 500}, PlotRange -> All,
AxesLabel -> {t, var[[i]]}], {i, Length[smax]}]


smax1 = wT1[mus0, mui0, 500] /. global[[2]];

var = {a, ns, ni, mus, mui, h, w}; Table[
Plot[smax1[[i]], {t, 0, 500}, PlotRange -> All,
AxesLabel -> {t, var[[i]]}], {i, Length[smax1]}]


Note, that solution computed with FindMaximum is not converges due to lack of precision. Therefore solutions shown above look different.

• Changing Log[x] to Log[real[x]] changes the problem, opening the way for negative a[t] without triggering complex values of other functions. Commented Apr 12 at 6:09
• @beginners Thank you very much for your explanation. Maybe you don't understand what your are looking for. If you are looking for global maximum, then you can't use u and aformula as it is and you should defined these functions for negative argument as well, for example as in my answer. If you are looking for a local maximum, then you should add constrain on solution, as for example in your PS. But in this case we can't use ParametricNDSolveValue and need some new method to solve the problem. Commented Apr 12 at 11:15
• Thank you. These are very helpful directions. Commented Apr 14 at 14:39
• @beginners Please, pay attention for eps in the function wT2 definition. This is computational error of the model which shows how close solution comes to singularity. Commented Apr 14 at 20:26