Method for constrained NonlinearModelFit with NDSolve

Question

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....

Thank you in advance

Answer 1

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}]
]