Numerical optimization yields strange result with imaginary number

Question

Edit: I changed the range for $r$ to $d/s \leq r \leq 1/2$ from the original one, $0 \leq r \leq 1/2$ and now it works! But there has emerged another issue, which is presented below.

Consider the following function: $$V=f(i(r),r;q,c,d,s)$$ That is, $V$ is a function of $i$ and $r$ with parameters $q$, $c$, $d$, $s$, and $i$ is a function of $r$. So ultimately it narrows down to $V$ as a function of $r$ along with the parameters. Given certain parameter values, I would like to numerically find $r$ that maximizes $V$ under the condition for the range of $r$ as $d/s \leq r \leq 1/2$.

Here is the code:

Clear["Global`*"];
c = 0.8; d = 0.8; q = 1.5; s = 4;
i = -0.1875` + 0.3125` r + (-10.5625` + (391.875` - 401.5625` r) r)/(-2197.` + r (-31335.` + (75825.` - 48125.` r) r) + 320.` \[Sqrt](r (6591.` + r (-193686.` + r (2.76172`*^6 + r (-7.51875`*^6 + (7.511625`*^6 - 2.5675`*^6 r) r))))))^(1/3) - 0.0625` (-2197.` + r (-31335.` + (75825.` - 48125.` r) r) + 320.` \[Sqrt](r (6591.` + r (-193686.` + r (2.76172`*^6 + r (-7.51875`*^6 + (7.511625`*^6 -2.5675`*^6 r) r))))))^(1/3);
v = (1 - r) y + r x;
V[r_?NumericQ] := NIntegrate[If[v - d (1 + i) > 0, v - d (1 + i), 0]/s^2, {x, 0, (d (1 + i))/r}, {y, (d (1 + i))/(1 - r) - r/(1 - r) x, s}] + NIntegrate[If[v - d (1 + i) > 0, v - d (1 + i), 0]/s^2, {x, (d (1 + i))/r, s}, {y, 0, s}]
NMaximize[{V[r], d/s <= r <= 1/2}, r]

which produces {1.13693, {r -> 0.40297}} as expected.

But to confirm the solution, when I evaluate V[] at the optimal r:

V[0.40297]

I get an error as follows:

Why this??

Answer 1

The accident in OP's code is interesting enough, so let me write an answer.

First let me explain why V[0.40297] "doesn't work". The issue can be boiled down to the following:

Why doesn't

aaa = 2 b; 
f[b_] := aaa; 
f[343]

output 686? How to fix this sample?

Well actually I've written a Chinese tutorial elaborating this issue and have been thinking about translating it to English for years but I'm just too lazy to do it. Now we have a detailed post discussing the topic, please have a look:

Why doesn't b = a; f[a_] := b; f[2] return 2?

If you don't bother to, then, in short, f[343] doesn't evaluate to 686 because the 2 b in aaa isn't explicitly there when the function f is defined, so it's not noticed. The most basic fix for this issue is to define aaa as a function of b to make the b in aaa explicit:

Clear[aaa, b, f]
aaa[b_] = 2 b;
f[b_] := aaa[b];
f[343]
(* 686 *)

And your V can be fixed in exactly the same manner i.e. make the r in i and v explicit:

Clear[i, v, r, x, y];
c = 0.8; d = 0.8; q = 1.5; s = 4;
i[r_] = -0.1875` + 
   0.3125` r + (-10.5625` + (391.875` - 401.5625` r)  r)/(-2197.` + 
      r  (-31335.` + (75825.` - 48125.` r)  r) + 
      320.`  √(r  (6591.` + 
            r (-193686.` + 
               r (2.76172`*^6 + 
                  r (-7.51875`*^6 + (7.511625`*^6 - 2.5675`*^6 r) r))))))^(1/3) - 
   0.0625`  (-2197.` + r  (-31335.` + (75825.` - 48125.` r)  r) + 
      320.`  √(r  (6591.` + 
            r (-193686.` + 
               r (2.76172`*^6 + 
                  r (-7.51875`*^6 + (7.511625`*^6 - 2.5675`*^6 r) r))))))^(1/3);
v[r_] = (1 - r) y + r x;
V[r_?NumericQ] := 
 With[{v = v[r], i = i[r]}, 
  NIntegrate[If[v - d (1 + i) > 0, v - d (1 + i), 0]/
    s^2, {x, 0, (d (1 + i))/r}, {y, (d (1 + i))/(1 - r) - (r x)/(1 - r), s}] + 
   NIntegrate[If[v - d (1 + i) > 0, v - d (1 + i), 0]/
    s^2, {x, (d (1 + i))/r, s}, {y, 0, s}]]

This isn't the only fix, of course. For example, the original i and v can be used if we control the evaluation order a bit:

(* The following is for the original v and i *)
(V[r_?NumericQ] := 
    With[{v = #, i = #2}, 
     NIntegrate[
       If[v - d  (1 + i) > 0, v - d  (1 + i), 0]/s^2, {x, 
        0, (d  (1 + i))/r}, {y, (d  (1 + i))/(1 - r) - (r  x)/(1 - r), s}] + 
      NIntegrate[
       If[v - d  (1 + i) > 0, v - d  (1 + i), 0]/s^2, {x, (d  (1 + i))/r, s}, {y, 0, 
        s}]]) &[v, i]

"OK, then why does the improper V in original code happen to work?" As mentioned above, this accident is the interesting part of this question, and it seems to be an issue that's never discussed in this site before! For better illustration, I'll again use a toy example:

Clear[f, r]

expr = 2 r^2; f[r_?NumericQ] := expr

NMinimize[f[r], r]
(* {2.79991*10^-17, {r -> -3.7416*10^-9}} *)

Clearly, f isn't proper, but NMinimize still manages to "work". What's happening here? The reason seems to be, an Experimental`NumericalFunction has been created internally:

Trace[NMinimize[f[r], r], _Experimental`CreateNumericalFunction, 
    TraceInternal -> True] // Flatten // Union // First
(* HoldForm[Experimental`CreateNumericalFunction[{r}, f[r], {}, 
     {WorkingPrecision -> MachinePrecision, Compiled -> Automatic}]] *)

And Experimental`NumericalFunction owns a really striking feature: it penetrates the black-box function defined with _?NumericQ!:

tst = Experimental`CreateNumericalFunction[{r}, f[r], {}];
tst[{3}]
(* 18. *)

I've added this to this community wiki.

Remark

It's reasonable to guess NMinimize has essentially called Block to penetrate the _?NumericQ. Indeed, Block is mentioned in document of numeric functions like NIntegrate, etc.:

NIntegrate has attribute HoldAll and effectively uses Block to localize variables. FindMinimum has attribute HoldAll, and effectively uses Block to localize variables.

But since Trace[NMinimize[f[r], r], _Block, TraceInternal -> True] doesn't find anything related, this is merely a guess.