# Method for constrained NonlinearModelFit with NDSolve

I'm trying to fit some data with a system of differential equations :

data={{6., 60.1536}, {9., 57.9807}, {12., 60.9089}, {15., 59.4291}, {18., 61.3227}, {21., 61.8788}, {24., 67.2192}, {27., 66.2767}, {30., 68.705}, {33., 68.2107}, {36., 70.8731}, {39., 68.7269}, {42.,
73.2306}, {45., 72.3068}, {48., 74.8006}, {51., 72.1975}, {54., 76.577}, {57., 75.5894}, {60., 76.342}, {63., 75.5134}, {66.,
77.47}, {69., 76.6854}, {72., 78.7422}, {75., 78.6074}, {78., 81.0158}, {81., 82.8521}, {84., 85.1395}, {87., 85.211}, {90.,
84.5157}, {93., 83.622}, {96., 88.1703}, {99., 85.6195}, {102., 86.8345}, {105., 86.5568}, {108., 87.5942}, {111., 88.3053}, {114.,
88.3475}, {117., 89.3993}, {120., 91.7091}, {123., 89.7268}, {126., 90.6704}, {129., 89.7999}, {132., 90.369}, {135., 88.7787}, {138.,
90.3022}, {141., 89.8267}, {144., 91.2241}, {147., 91.2859}, {150., 92.992}, {153., 91.0079}, {156., 93.0784}, {159., 90.8868}, {162.,
92.7426}, {165., 92.757}, {168., 94.4202}, {171., 92.2914}, {174.,
90.3876}, {177., 89.3376}, {180., 89.814}, {183., 88.9134}, {186.,
89.7058}, {189., 91.642}, {192., 90.3205}, {195., 87.8566}, {198., 87.6065},            {201., 87.0403}, {204., 87.3344}, {207., 87.2313}, {210.,
87.3705}, {213., 86.9135}, {216., 87.2684}, {219., 87.2989}, {222.,
85.4766}, {225., 85.3534}, {228., 86.535}, {231., 86.1929}, {234.,
86.089}, {237., 85.9466}, {240., 85.1389}, {243., 85.0242}, {246.,
84.4313}, {249., 83.7604}, {252., 81.9419}, {255., 83.773}, {258.,
82.7046}, {261., 84.7331}, {264., 86.0393}, {267., 84.7472}, {270.,
79.1677}, {273., 80.9426}, {276., 79.9624}, {279., 75.5272}, {282.,
79.3103}, {285., 80.8015}, {288., 81.3927}, {291., 80.1678}, {294.,
80.268}, {297., 79.9067}, {300., 76.9766}, {303., 81.8132}, {306.,
73.6449}, {309., 76.4059}, {312., 76.4056}, {315., 81.7311}, {318.,
80.8468}, {321., 80.958}, {324., 86.9248}, {327., 78.3434}, {330.,
74.8752}, {333., 78.0912}, {336., 81.5165}, {339., 72.7919}, {342.,
74.2966}, {345., 79.2233}, {348., 81.9791}, {351., 74.3276}, {354.,
85.1221}, {357., 78.6944}, {360., 75.8183}, {363., 75.6696}, {366.,
75.9147}, {369., 76.3326}, {372., 80.0048}, {375., 79.8311}, {378.,
79.0427}, {381., 81.8084}, {384., 73.5742}, {387., 84.2291}, {390.,
84.9122}, {393., 82.6657}, {396., 78.2888}, {399., 90.0235}, {402.,
83.3667}, {405., 81.7737}, {408., 81.19}, {411., 82.3131}, {414.,
79.8072}, {417., 74.4822}, {420., 75.6291}, {423., 82.2655}, {426.,
73.704}, {429., 81.4184}, {432., 72.1127}, {435., 74.7053}, {438.,
79.4664}, {441., 86.4491}, {444., 79.5096}, {447., 77.1761}, {450.,
83.082}, {453., 80.3418}, {456., 85.3873}, {459., 85.7409}, {462.,
73.3735}, {465., 72.2276}, {468., 82.7752}, {471., 71.6917}, {474.,
78.5233}, {477., 82.4042}, {480., 83.8073}, {483., 91.5845}, {486.,
82.8906}, {489., 87.3935}, {492., 89.9856}, {495., 74.1819}, {498.,
77.5752}, {501., 82.6796}, {504., 79.2659}, {507., 81.5865}, {510.,
82.709}, {513., 88.4083}, {516., 81.7317}, {519., 76.2638}, {522.,
86.2863}, {525., 93.2163}, {528., 82.6943}}


And now the equations :

Rc = 50

model[R0_?NumberQ, c_?NumberQ, d_?NumberQ] :=  Module[{y, m, x},
First[y /. NDSolve[{y'[x] == (0.20 - 0.005*(y[x] - Rc) - m[x]*y[x]),
m'[x] == c - d*m[x], y[0] == R0, m[0] == 0}, {y, m}, {x, 0,
530}, Method -> {"EquationSimplification" -> "Residual"}]]]

nlm = NonlinearModelFit[data, {model[59, c, d][x], 0 < c, 0 < d }, {{c, 0.5}, {d, 1}}, x]


But I've been waiting for a very long time and Mathematica didn't give me the solutions.

So I wanted to ask you if there was some "Method-> " I could add in my constrained NonlinearModelFit in order to get faster the result. I looked for such a method but what I read left me very unclear on how to do it exactly... Apparently most of the methods are not ok for constrained problems, and I tried Method->"NMinimize" which did not change anything. But apparently it has to be written : Method-> {"NMinimize",Method->...} and I don't know what to write exactly....

• I think you need a better model or better starting values. If you plot the predictions from your current starting values (ListPlot[Transpose[{data[[All, 1]], model[59, 0.5, 1][#] & /@ data[[All, 1]]}]]), you won't get a curve that looks anything like the data. Also, for down the road, the variance about the curve is not constant (one of the assumptions with the defaults of NonlinearModelFit) but that can be examined after a reasonable fit (assuming constant variance) is found.
– JimB
Commented Apr 15, 2019 at 15:38
• @JimB Thank you for those explanations. Could you explain me please what : ListPlot[Transpose[{data[[All, 1]], model[59, 0.5, 1][#] & /@ data[[All, 1]]}]] exactly means ?
– J.A
Commented Apr 15, 2019 at 15:46
• The predictions for your starting values are found with model[59, 0.5, 1][#] & /@ data[[All, 1]]. The # symbol is replaced with each value of the predictor (data[[All,1]]). Then Transpose[{x,y}] constructs the predicted dataset.
– JimB
Commented Apr 15, 2019 at 16:01
• Note that @CarlWoll 's answer suggests the current model is not a very good fit. In addition, looking at fit["CorrelationMatrix"] gives a result of $\left( \begin{array}{cc} 1. & 1. \\ 1. & 1. \\ \end{array} \right)$ which can suggest an overparameterized model or that lots of combinations of $c$ and $d$ give the same (poor) fit.
– JimB
Commented Apr 15, 2019 at 16:11

I would use ParametricNDSolveValue to create the model, because then parameter derivatives of the model can be used when creating the fit:

Rc = 50;
model = ParametricNDSolveValue[
{
y'[x] == 0.20 - 0.005 (y[x] - Rc) - m[x] y[x],
m'[x] == c - d m[x],
y[0] == R0,
m[0] == 0
},
y,
{x, 0, 530},
{R0, c, d}
];


Then, use NonlinearModelFit:

fit = NonlinearModelFit[
data,
{model[59, c, d][x], 0 < c, 0 < d},
{{c, .5}, {d, 1}},
x
]; //AbsoluteTiming


{1.17995, Null}

Finally, plot the data and the fit:

Show[
ListPlot[data],
Plot[fit[x], {x, 0, 550}]
]