# How can I minimize a target function that has to solve differential equations numerically?

My problem is to minimize a target function with differential equations in it.

The target function is

$$\min \mathrm{J}\left( p \right) =\sum_{i=1}^5{|x_i\left( p \right) -y_i|},$$

where $$x_i(p)$$ are differential equations like this:

stateFunctionV[v1_, v2_, v3_, v4_, v5_, v6_, v7_, v8_] :=
{μ =
(67/100)(x2[t]/(x2[t] + 28/100))(1 - x3[t]/940)(1 - x4[t]/1050)^3 (1 - x5[t]/361);
x1'[t] == μ*x1[t],
x2'[t] == -(v1 + μ/v5)* x1[t],
x3'[t] == (v2 + μ*v6)*x1[t],
x4'[t] == (v3 + μ*v7)* x1[t],
x5'[t] == (v4 + μ*v8)* x1[t],
x1 == 0.2025, x2 == 441.337, x3 == 0, x4 == 0, x5 == 0}


and $$y_i$$ is the experiment data.

The code:

PIP1[value_] :=
{Sum[
Sum[
Abs[
NDSolveValue[{
μ = (67/100)(x2[t]/(x2[t] + 28/100))(1 - x3[t]/940)(1 -       x4[t]/1050)^3(1 - x5[t]/361);
x1'[t] == μ*x1[t],
x2'[t] == -(value[] + μ/value[])x1[t],
x3'[t] == (value[] + μ value[])x1[t],
x4'[t] == (value[] + μ value[])x1[t],
x5'[t] == (value[] + μ value[])x1[t],
x1 == 0.2025, x2 == 441.337, x3 == 0, x4 == 0,
x5 == 0}, {x1, x2, x3, x4, x5}, {t, 0, 6.5}][[i]][data2[[i, j, 1]]] - data2[[i, j, 2]]],
{j, 1, 7}],
{i, 1, 5}]}

NMinimize[PIP1[{v1, v2, v3, v4, v5, v6, v7, v8}[]], {v1, v2, v3, v4, v5, v6, v7, v8}]


doesn't work. I also tried to take the differential functions out as a new function and it also failed.

So I wonder if there's a way to minimize this problem with NMinimize or not? I would be grateful for your advice on how to couple NMinimize and NDSolveValue.

The experiment data are below:

data2 =
{{{0, 0.2025}, {2, 0.76}, {3, 1.28}, {4, 2.2}, {5, 2.75}, {6, 2.95}, {6.5, 3.055}},
{{0, 441.337}, {2, 400.435}, {3, 300.326}, {4, 219.348}, {5, 120.109}, {6, 50}, {6.5, 25}},
{{0, 0}, {2, 42.2895}, {3, 78.6579}, {4, 117.026}, {5, 152.158}, {6, 229.553}, {6.5, 224.1205}},
{{0, 0}, {2, 1.66667}, {3, 13.55}, {4, 16.6667}, {5, 22.7333}, {6, 34.4333}, {6.5, 71.2333}},
{{0, 0}, {2, 6.73913}, {3, 7.86957}, {4, 9.3913}, {5, 8.34783}, {6, 10.9565}, {6.5, 11.8261}}}; The target function is to minimize the sum of differences between the experimental data at a given time (0,1,2,3,4,5,6,6.5h) and the solutions of the differential equations at that moment.

• You might want to look at using ParametricNDSolveValue[] instead. Also, FindFit[] should be able to handle this if you change the "NormFunction" option. – J. M.'s torpor Jan 27 at 16:24
• As defined, IP1 returns a list, but it should return a number. Therefore, replace the curly brackets by round brackets – Daniel Huber Jan 27 at 17:18
• @J.M. Thank you so much! Could I use FindFit with the dataset that has 5 outputs at a given time? Should I manipulate the dataset and how? :) – Charmbracelet Jan 28 at 7:46
• @DanielHuber Thank you, I edited the question. – Charmbracelet Jan 28 at 10:49

## 1 Answer

Clear[obj]
data2 = {{{0, 0.2025}, {2, 0.76}, {3, 1.28}, {4, 2.2}, {5, 2.75}, {6, 2.95}, {6.5, 3.055}}, {{0, 441.337}, {2, 400.435}, {3, 300.326}, {4, 219.348}, {5, 120.109}, {6, 50}, {6.5, 25}}, {{0, 0}, {2, 42.2895}, {3, 78.6579}, {4, 117.026}, {5, 152.158}, {6, 229.553}, {6.5, 224.1205}}, {{0, 0}, {2, 1.66667}, {3, 13.55}, {4, 16.6667}, {5, 22.7333}, {6, 34.4333}, {6.5, 71.2333}}, {{0, 0}, {2, 6.73913}, {3, 7.86957}, {4, 9.3913}, {5, 8.34783}, {6, 10.9565}, {6.5, 11.8261}}};

mu = (67/100) (x2[t]/(x2[t] + 28/100)) (1 - x3[t]/940) (1 - x4[t]/1050)^3 (1 - x5[t]/361);
ODE = {x1'[t] == mu*x1[t],
x2'[t] == -(v1 + mu/v5)*x1[t],
x3'[t] == (v2 + mu*v6)*x1[t],
x4'[t] == (v3 + mu*v7)*x1[t],
x5'[t] == (v4 + mu*v8)*x1[t],
x1 == 0.2025, x2 == 441.337, x3 == 0, x4 == 0, x5 == 0};
pfun = ParametricNDSolveValue[ODE, {x1, x2, x3, x4, x5}, {t, 0, 6.5}, {v1, v2, v3, v4, v5, v6, v7, v8}];
weight = {1, 1/400, 1/40, 1, 1/6};
obj[v1_?NumberQ, v2_?NumberQ, v3_?NumberQ, v4_?NumberQ, v5_?NumberQ, v6_?NumberQ, v7_?NumberQ, v8_?NumberQ] := Sum[weight[[i]]* Sum[Abs[pfun[v1, v2, v3, v4, v5, v6, v7, v8][[i]][data2[[i]][[j, 1]]] - data2[[i]][[j, 2]]], {j, 1, 7}], {i, 1, 5}]
solvs = NMinimize[obj[v1, v2, v3, v4, v5, v6, v7, v8], {v1, v2, v3, v4, v5, v6, v7, v8}, Method -> "DifferentialEvolution"]
pfun0 = pfun[v1, v2, v3, v4, v5, v6, v7, v8] /. solvs[];

grs = Table[Plot[pfun0[[i]][t], {t, 0, 6.5}], {i, 1, 5}];
pts = Table[ListPlot[data2[[i]], PlotStyle -> Red], {i, 1, 5}];
Show[grs, pts, PlotRange -> All] 