Minimization with integration over region

Question

I'm having some difficulties implementing a minimization problem were the objective function involves a numeric integration. This simplified problem highlights one of the issues I'm currently having: how to find the location of the minimum average value of a function. For this problem assume the function can be described by

exactfun[x_?NumericQ, y_?NumericQ] := (x - .5)^4 + (y - .5)^4

I want to minimize the average value of the function over a disk region centered at $\left(a,b\right)$ with radius $r$. The exact value of $f_{avg} \left(a,b,r\right)$ can be determined explicitly for this particular function, though in general (and in my application), it is not possible. Here, $f_{avg} \left(a,b,r\right)$ is

FullSimplify[
 Integrate[(x - 1/2)^4 + (y - 1/2)^4, {x, y} \[Element] 
    Disk[{a, b}, r], 
   Assumptions -> {{a, b} \[Element] Reals && r > 0}]/
  Area@Disk[{a, b}, r]]
(*1/8 (1 + 4 (-1 + b) b (1 + 2 (-1 + b) b) + 6 r^2 + 
   2 (2 (-1 + a) a (1 + 2 (-1 + a) a) + 
      6 ((-1 + a) a + (-1 + b) b) r^2 + r^4))*)

Minimizing with respect to the disk center yields

Solve[{D[1/
     8 (1 + 4 (-1 + b) b (1 + 2 (-1 + b) b) + 6 r^2 + 
       2 (2 (-1 + a) a (1 + 2 (-1 + a) a) + 
          6 ((-1 + a) a + (-1 + b) b) r^2 + r^4)), a] == 0, 
  D[1/8 (1 + 4 (-1 + b) b (1 + 2 (-1 + b) b) + 6 r^2 + 
       2 (2 (-1 + a) a (1 + 2 (-1 + a) a) + 
          6 ((-1 + a) a + (-1 + b) b) r^2 + r^4)), b] == 0}, {a, 
  b}, Reals]
(*{{a -> 1/2, b -> 1/2}}*)

I'm interested in doing this same problem numerically. However, in the evaluation of NMinimize there is an error regarding integration limits that I can't sort out. I tried constructing the optimization function two ways, and each failed.

numericint[a_?NumericQ, b_?NumericQ, r_?NumericQ] := 
 NIntegrate[exactfun[x, y], {x, y} \[Element] Disk[{a, b}, r]]/
  Area@Disk[{a, b}, r]
numericint2[a_?NumericQ, b_?NumericQ, r_?NumericQ] := 
 Module[{disk = Disk[{a, b}, r]}, 
  NIntegrate[exactfun[x, y], {x, y} \[Element] disk, 
    Method -> {Automatic, "SymbolicProcessing" -> False}]/Area@disk]

Running the optimization for a particular $r$ and in a certain region fails

NMinimize[
 numericint[x, y, .05], {x, y} \[Element] Rectangle[{0, 0}, {1, 1}]]
(*NIntegrate::ilim: Invalid integration variable or limit(s) in True.*)
NMinimize[
 numericint2[x, y, .05], {x, y} \[Element] Rectangle[{0, 0}, {1, 1}]]
(*NIntegrate::ilim: Invalid integration variable or limit(s) in True.*)

Each of the functions appear to evaluate fine outside of the NMinimize environment, but are crashing and burning during the minimization. For reference:

numericint[.25, .25, .05]
(*0.00828281*)
numericint2[.25, .25, .05]
(*0.00828281*)
1/8 (1 + 4 (-1 + b) b (1 + 2 (-1 + b) b) + 6 r^2 + 
    2 (2 (-1 + a) a (1 + 2 (-1 + a) a) + 
       6 ((-1 + a) a + (-1 + b) b) r^2 + r^4)) /. {a -> .25, b -> .25,
   r -> .05}
(*0.00828281*)

Answer 1

The issue seems to be with the specification of the region of integration in the form {x, y} ∈ Disk[{a, b}, r], which is evaluating to True when being evaluated inside the NMinimize command. The problem persists when using the FindMinimum command as well. It can be worked around, however, by specifying the region of integration in terms of explicit bounds on x and y. For a circle of radius $r$ centered at $(a,b)$, we have \begin{align} a - r \leq &\: x \leq a + r, \\ b - \sqrt{ r^2 - (x - a)^2} \leq &\: y \leq b + \sqrt{ r^2 - (x - a)^2} \end{align} and we can program these bounds in explicitly:

numericint3[a_?NumericQ, b_?NumericQ, r_?NumericQ] := 
 NIntegrate[
   exactfun[x, y], {x, a - r, a + r}, {y, b - Sqrt[r^2 - (x - a)^2], 
    b + Sqrt[r^2 - (x - a)^2]}]/(π r^2)

FindMinimum[{numericint3[x, y, .05], 0 <= x <= 1, 0 <= y <= 1}, {x, y}]

(* {1.5625*10^-6, {x -> 0.5, y -> 0.5}} *)

Using FindMinimum takes about 6 seconds to evaluate on my machine. On the other hand, my attempt to use NMinimize to do the same thing is still running after about 5–10 minutes. I'll update this answer when & if it finishes.

UPDATE: NMinimize finished after 13m05s on my machine, yielding the same result:

NMinimize[
 numericint3[x, y, .05], {x, y} ∈ Rectangle[{0, 0}, {1, 1}]]

(* {1.5625*10^-6, {x -> 0.5, y -> 0.5}} *)