FindMinimum running forever

Question

I cannot get any output from the below-mentioned optimization problem. The program keeps running forever. What can I do to solve this optimization problem successfully?

ClearAll["Global`*"];
b=19;
a=15;   
es=1;                                                 (*Element size for FEA*)
noofnode=IntegerPart[((b-a)/es)+1];
init=0&;
eb=380*10^9;     (*MPa*)
ea=71*10^9;       (*MPa*)
new=0.3;                             (*Poisson's Ratio*)                  
alphab=8.0*10^(-6);
alphaa=23.1*10^(-6);
rowb=0.96;
rowa=2.70;
ta=20;                      (*Degree C*)
tb=150;  
rpm=150.0;
wom=(2.0*\[Pi]*rpm)/60;
kb=35.0;
ka=204.0;
lam=(1/(b-a))*Log[kb/ka];               (*Constant for Exponential Temp Profile*)
f[ri_]:=ExpIntegralEi[-lam*ri];
c1=(tb-ta)/(f[b]-f[a]);
c2=((ta*f[b])-(tb*f[a]))/(f[b]-f[a]);
T[r_]:=c2+c1*ExpIntegralEi[-lam*r];
ed[v_]:=(v*ea)+((1-v)*eb);
alphad[v_]:=(v*alphaa)+((1-v)*alphab);
rowd[v_]:=(v*rowa)+((1-v)*rowb);
e[v1_,v2_,v3_,v4_,v5_]:=Interpolation[{ed[v1],ed[v2],ed[v3],ed[v4],ed[v5]}];
alpha[v1_,v2_,v3_,v4_,v5_]:=Interpolation[{alphad[v1],alphad[v2],alphad[v3],alphad[v4],alphad[v5]}];
row[v1_,v2_,v3_,v4_,v5_]:=Interpolation[{rowd[v1],rowd[v2],rowd[v3],rowd[v4],rowd[v5]}];
z[v1_?NumericQ,v2_?NumericQ,v3_?NumericQ,v4_?NumericQ,v5_?NumericQ]:=NDSolveValue[{(3+new) r^2 wom^2 row[v1,v2,v3,v4,v5][r-14]+F'[r]+(e[v1,v2,v3,v4,v5]'[r-14] (new F[r]-r (r^2 wom^2 row[v1,v2,v3,v4,v5][r-14]+F'[r])))/e[v1,v2,v3,v4,v5][r-14]+r^3 wom^2 row[v1,v2,v3,v4,v5]'[r-14]+r e[v1,v2,v3,v4,v5][r-14] (T[r] alpha[v1,v2,v3,v4,v5]'[r-14]+alpha[v1,v2,v3,v4,v5][r-14] T'[r])+r F''[r]==F[r]/r,F[a]==0,F[b]==0},F,{r,a,b}];
o[v1_,v2_,v3_,v4_,v5_]:=Sqrt[Sum[(((z[v1,v2,v3,v4,v5][j])*10^(-9))^2),{j,15,19}]/5];
FindMinimum[{o[v1,v2,v3,v4,v5],{0<=v1<=1,0<=v2<=1,0<=v3<=1,0<=v4<=1,0<=v5<=1}},{v1,v2,v3,v4,v5}]

Answer 1

This problem can be solved with using the Euler wavelets collocation method explained in our paper and on the page. First we map region $a\le r\le b$ on the unit interval by substitution r=a+(b-a) x in GE. We also define F2=F'',F1=F',v1=v', then we have

ClearAll["Global`*"];
b = 19;
a = 15;
eb = 380*10^9;(*MPa*)ea = 
 71*10^9;(*MPa*)new = 0.3;(*Poisson's Ratio*)alphab = 8.0*10^(-6);
alphaa = 23.1*10^(-6);
rowb = 0.96;
rowa = 2.70;
ta = 20;(*Degree C*)tb = 150;
rpm = 150.0;
wom = (2.0*\[Pi]*rpm)/60;
kb = 35.0;
ka = 204.0;
L = b - a;
lam = (1/(b - a))*Log[kb/ka];(*Constant for Exponential Temp Profile*)
f[ri_] := ExpIntegralEi[-lam*ri];
c1 = (tb - ta)/(f[b] - f[a]);
c2 = ((ta*f[b]) - (tb*f[a]))/(f[b] - f[a]);
T = c2 + c1*ExpIntegralEi[-lam*r]; T1 = c1 E^(-lam r)/r;
ed = (v*ea) + ((1 - v)*eb);
alphad = (v*alphaa) + ((1 - v)*alphab);
rowd = (v*rowa) + ((1 - v)*rowb);
e1 = D[ed, v] v1[x]/L;
alpha1 = D[alphad, v] v1[x]/L;
row1 = D[rowd, v] v1[x]/L;

 eq = {(3 + new) r^2 wom^2 rowd + 
     F1[x]/L + (e1 (new F[x] - r (r^2 wom^2 rowd + F1[x]/L)))/ed + 
     r^3 wom^2 row1 + r ed (T alpha1 + alphad T1/L) + r F2[x]/L^2 - 
     F[x]/r, F[a] == 0, F[b] == 0} /. {v -> v[x], r -> a + L x, 
    F[a] -> F[0], F[b] -> F[1]} // Simplify;

Second, we define wavelets and all functions to be computed

UE[m_, t_] := EulerE[m, t]
psi[k_, n_, m_, t_] := 
 Piecewise[{{2^(k/2) Sqrt[2/Pi] UE[m, 2^k t - 2 n + 1], (n - 1)/
      2^(k - 1) <= t < n/2^(k - 1)}, {0, True}}]
PsiE[k_, M_, t_] := 
 Flatten[Table[psi[k, n, m, t], {n, 1, 2^(k - 1)}, {m, 0, M - 1}]]
k0 = 2; M0 = 3; With[{k = k0, M = M0}, 
 var0 = Flatten[Table[cvar[n, m], {n, 1, 2^(k - 1)}, {m, 0, M - 1}]]];
nn = Length[var0];
dx = 1/(nn);  xl = Table[ l*dx, {l, 0, nn}]; zcol = 
 xcol = Table[(xl[[l - 1]] + xl[[l]])/2, {l, 2, 
    nn + 1}]; ycol = zcol; Psijk = 
 With[{k = k0, M = M0}, PsiE[k, M, t1]]; Int1 = 
 With[{k = k0, M = M0}, Integrate[PsiE[k, M, t1], t1]];
Int2 = Integrate[Int1, t1]; Psi[y_] := Psijk /. t1 -> y; 
int1[y_] := Int1 /. t1 -> y; int2[y_] := Int2 /. t1 -> y;
M = nn; U1 = Array[u1, {M}]; U2 = Array[u2, {M}];
F2[x_] := U1 . Psi[x]; F1[x_] := U1 . int1[x] + u0; 
F[x_] := U1 . int2[x] + u0 x + u00; v1[x_] := U2 . Psi[x]; 
v[x_] := U2 . int1[x] + v0;

Now we can solve GE and define F as function of v in collocation points

var = Join[U1, {u0, u00}]; eq1 = 
 Join[Table[eq[[1]] == 0, {x, xcol}], {eq[[2]], eq[[3]]}];
  

sol=Solve[eq1, var];
       
opt = Table[F[x] /. sol[[1]], {x, xcol}];

Finally we solve optimization problem with constraints $0\le v\le 1$ as follows

var2 = Join[U2, {v0}]; res = 
 Table[{v[x] >= 0, v[x] <= 1}, {x, Join[{0}, xcol, {1}]}]//Flatten; sol1 = NMinimize[{opt . opt, res}, var2];

Visualization

    {Plot[Evaluate[v[(r - a)/L] /. sol1[[2]]], {r, a, b}, 
  PlotRange -> All, Frame -> True, PlotLabel -> "v"], 
 Plot[Evaluate[F[(r - a)/L] /. sol[[1]] /. sol1[[2]]], {r, a, b}, 
  PlotRange -> All, Frame -> True, PlotLabel -> "F"]}

Update 1. We can force computation with option Method ->"DifferentialEvolution", for instance, for k0 = 2; M0 = 4;we define

scale = 
 Evaluate[opt[[1]] /. 
   Table[var2[[i]] -> RandomReal[{-1, 1}], {i, Length[var2]}]]

(*Out[]= -2.28631*10^7*)

sol1 = 
 NMinimize[{Norm[opt/scale], res}, var2, 
  Method -> "DifferentialEvolution"]

(*Out[]= {0.0000203025, {u2[1] -> -0.399257, u2[2] -> 0.107406, 
  u2[3] -> -0.0285942, u2[4] -> 0.0115987, u2[5] -> -0.267057, 
  u2[6] -> 0.042192, u2[7] -> -0.00704157, u2[8] -> 0.00169105, 
  v0 -> 1.}}*)

Visualization