# optimization problem with NDSolve and Import cannot get any results, even without error message

I want to minimize the function ff. When ff is calculated for a specified point the answer is correct:

ff[4*10^(-4), 6*10^(-4)]
ff[3*10^(-4), 5*10^(-4)]
Out=-0.916638
Out=-0.91548


but when I want to optimize ff with NMinimize, it just keeps running, but don't show any results, even without error message, which is the most strange thing. First I wonder if something is wrong of the function ff. I also wonder if it's because my codes are too complicated due to the use of Import of data. I have tried to use different Methods in Nminimize, but it doesn't work. I'm new using Mathematica. Could anyone give me some advice, please?

resNMinimize =
Reap[NMinimize[{ff[q1, q2],
1*10^(-5) <= q1 <= 0.001 && 1*10^(-5) <= q2 <= 0.001}, {q1, q2},
(* Method -> "SimulatedAnnealing", AccuracyGoal -> 9,
PrecisionGoal -> 8,EvaluationMonitor\[RuleDelayed]Sow[{q1,q2}]*)(*,
StepMonitor\[RuleDelayed]Sow[{q1,q2}],MaxIterations\[Rule]40*)]]


and this is the ff function:

Clear["Global*"];
n = 6; ke = 0; W = 1.22; wh = 0.12; Wd = 0.005; l = 4.125; h = 0.072;
a0 = -1.112 ; a1 = -0.8115 ; b1 = 2.52 ; a2 = 1.609 ; b2 = 1.156 ; a3 \
= 0.9558 ; b3 = -0.7058 ; a4 = -0.1763 ; b4 = -0.5241 ; a5 = -0.1802 \
; b5 = 0.00232 ; a6 = -0.01007 ; b6 = 0.0327 ; w = 0.03138 ;
u = a0 + a1*Cos[y*w*W] + b1*Sin[y*w*W] + a2*Cos[2*y*w*W] +
b2*Sin[2*y*w*W] + a3*Cos[3*y*w*W] + b3*Sin[3*y*w*W] +
a4*Cos[4*y*w*W] + b4*Sin[4*y*w*W] + a5*Cos[5*y*w*W] +
b5*Sin[5*y*w*W] + a6*Cos[6*y*w*W] + b6*Sin[6*y*w*W];
uav = Integrate[u, {y, 0, W}]/W;
U = u/uav;
v = -Integrate[D[u, x], {y, 0, y}];
V = v/uav;
t = ke*W/uav; h1 = wh/W; h2 = (Wd + wh)/W; theta = 0; Q = 1;
fai = Table[Cos[m*Pi*y], {m, 0, n}];
NN = N[Integrate[fai^2, {y, 0, 1}]];
xm = l/W;
obsdata = Import["https://pastebin.com/raw/Fz38738r", "Table"];
ff[e1_?NumericQ, e2_?NumericQ] := Module[{},
(*dy1=e1/(uav*W);*)
eqns = Table[{FullSimplify[
Sum[cc[j]'[x]/Sqrt[NN[[j + 1]]*NN[[m + 1]]]*
Chop[Integrate[fai[[j + 1]]*fai[[m + 1]]*U, {y, 0, 1}]], {j,
0, n}]] ==
FullSimplify[
Sum[cc[j][x]/
Sqrt[NN[[j + 1]]*NN[[m + 1]]]*(Chop[
Integrate[
fai[[j + 1]]*D[e1/(uav*W)*Dt[fai[[m + 1]], y], y], {y,
0, 0.5}]] +
Chop[Integrate[
fai[[j + 1]]*D[e2/(uav*W)*Dt[fai[[m + 1]], y], y], {y,
0.5, 1}]] -
Integrate[
V*fai[[m + 1]]*Dt[fai[[j + 1]], y], {y, 0, 1}]), {j, 0,
n}] - t*cc[m][x]],
cc[m] ==
N[Integrate[fai[[m + 1]]/Sqrt[NN[[m + 1]]]*theta, {y, 0, h1}]] +
N[Integrate[fai[[m + 1]]/Sqrt[NN[[m + 1]]]*Q, {y, h1, h2}]] +
N[Integrate[
fai[[m + 1]]/Sqrt[NN[[m + 1]]]*theta, {y, h2, 1}]]}, {m, 0,
n}];
vars = Table[cc[j], {j, 0, n}];
s = NDSolve[eqns, vars, {x, 0, xm},
Method -> {"IndexReduction" -> Automatic,
"EquationSimplification" -> {"Residual", "TimeConstraint" -> 10,
"SimplifySystem" -> True}}];
s1 = Flatten[s];
prediction =
Table[Sum[
fai[[m + 1]]/Sqrt[NN[[m + 1]]]*(cc[m][x] /. s1), {m, 0, n}], {y,
1/2*1/122, (1 - 1/2*1/122), 1/122}, {x, 1/2*25/122,
400/122 + 25/122/2, 25/122}];
-Correlation[Flatten[obsdata], Flatten[prediction]]]


And below is the content of conc_D72_Inj12.mat:

• I can't test your function because I do not have conc_D72_Inj12.mat. But you should certainly import the data once for all, i.e. outside the definition of ff. Here you import it everytime you evaluate ff! – anderstood Sep 11 '17 at 23:32
• Thanks for your advice. I will change to import the data outside the function of ff. – M.L. Sep 12 '17 at 17:13
• If you need further help, you should include this missing file (or replace it with dummy data). – anderstood Sep 12 '17 at 17:28
• Hi,Anderstood. How should I include the data? I don't find the place to attach files. – M.L. Sep 12 '17 at 17:52
• I would recommend using pastebin.com and providing the command to include the data (something like Import[www.pastebin.com/raw/...,"Table"]). Don't forget to update your code accordingly. As a general rule, anyone should be able to observe your problem and evaluate ff by simply copying and pasting your code, without having to spend time fiddling to reproduce your issue. – anderstood Sep 12 '17 at 18:03

You really need to understand that everything which can be evaluated once, should be evaluated once. So I started to rewrite your function ff by defining the set of equations, as a function of e1,e2, "globally" (outside of the definition of ff), same for prediction which took almost one second to evaluate every time:

rhs1[e1_, e2_] :=
Sum[cc[j][x]/Sqrt[NN[[j + 1]]*NN[[m + 1]]]*(Integrate[
fai[[j + 1]]*D[e1/(uav*W)*Dt[fai[[m + 1]], y], y], {y, 0,
0.5}] + Integrate[
fai[[j + 1]]*D[e2/(uav*W)*Dt[fai[[m + 1]], y], y], {y, 0.5,
1}] - Integrate[
V*fai[[m + 1]]*Dt[fai[[j + 1]], y], {y, 0, 1}]), {j, 0, n}] - t*cc[m][x];
lhs1 := Sum[cc[j]'[x]/Sqrt[NN[[j + 1]]*NN[[m + 1]]]*
NIntegrate[fai[[j + 1]]*fai[[m + 1]]*U, {y, 0, 1}], {j, 0, n}];
rhs2 := Integrate[fai[[m + 1]]/Sqrt[NN[[m + 1]]]*theta, {y, 0, h1}] +
N[Integrate[fai[[m + 1]]/Sqrt[NN[[m + 1]]]*Q, {y, h1, h2}]] +
N[Integrate[fai[[m + 1]]/Sqrt[NN[[m + 1]]]*theta, {y, h2, 1}]];
eqns[e1_, e2_] =
Table[{lhs1 == rhs1[e1, e2], cc[m] == rhs2}, {m, 0, n}]

vars = Table[cc[j], {j, 0, n}]
prediction = Flatten@Table[Sum[fai[[m + 1]]/Sqrt[NN[[m + 1]]]*(cc[m][x]),
{m, 0, n}],{y, 1/2*1/122, (1 - 1/2*1/122), 1/122}, {x, 1/2*25/122,
400/122 + 25/122/2, 25/122}]


Then, ff becomes much more readable:

ff[e1_?NumericQ, e2_?NumericQ] := Module[{},
s = Flatten@NDSolve[eqns[e1, e2], vars, {x, 0, xm}];
-Correlation[Flatten[obsdata], prediction /. s]]


Evaluating one ff takes about a tenth of a second. You can plot it:

Plot3D[ff[q1, q2], {q1, 0.00001, 0.001}, {q2, 0.00001, 0.001},
PlotPoints -> 2]


From the plot you can guess that the constrained minimum is close to $(0.001,0.00001)$ which you can confirm using FindMinimum:

FindMinimum[{ff[q1, q2],
1*10^(-5) <= q1 <= 0.001 && 1*10^(-5) <= q2 <= 0.001}, {{q1,
0.0005}, {q2, 0.0005}}]
(* {-0.920901, {q1 -> 0.000839124, q2 -> 0.0000100353}} *)


You can plot the minimum together with the surface:

Show[Plot3D[ff[q1, q2], {q1, 0.00001, 0.001}, {q2, 0.00001, 0.001},
PlotPoints -> 4],
Graphics3D[{Red, PointSize[0.05],
Point[({q1, q2} /. min[])~Join~{min[]}]}]] Note that FindMinimum returns a local minimum, but your function seems to be convex so it should not matter.

Take-home message Optimizing procedures simply need to evaluate a cost function (here, ff). The cost function just needs to output a real number, from values of input parameters. What is inside it really does not matter, you can import data, do whatever you need, as long as in the end it returns a real number. Of course, if the cost function is very long to evaluate, the minimization procedure might take an unreasonable amount of time, because it has to call the cost function many times (in particular if you have many parameters and need to compute the Jacobian of the cost function...).

• Wow! You are brilliant. I will study it carefully. Thank you very much. – M.L. Sep 13 '17 at 18:53
• Sorry to disturb you again. I have another question. Actually, I have large numbers of parameters, like e1,e2,...e100. How should I include them in my function ff? I have used Table, but the effect is not good. – M.L. Sep 14 '17 at 17:44
• @J.L. I would start by extending the code exactly as it is now, just by adding more variables. But optimizing over 100 variables is not an easy task... If you have more questions, you should start a new thread. – anderstood Sep 14 '17 at 17:53
• Yes, it should be not easy to optimize over 100 variables. Ok, I will start with a new question. – M.L. Sep 14 '17 at 17:57
• @Mr.Wizard Apparently the OP reverted again. The line eq = ` is not valid. I personally don't mind much, but that's not really the purpose of the website ("help me and then I'll remove details in my question"). – anderstood Dec 22 '17 at 18:05