Shape with area $a$ and shortest average distance between any two points Consider the collection $\mathcal{X}$ of 


*

*compact and connected subsets of $\mathbb{R}^2$,  

*that are the closure of some open subset of $\mathbb{R}^2$, and

*have area $a$.


For any $\Omega \in \mathcal{X}$ and any $x,y \in \Omega$, let $d(x,y,\Omega)$ be the length of the shortest path lying entirely inside $\Omega$ and connecting $x$ to $y$.
Also, let $\mathbb{1}(x,y,\Omega)$ be the indicator function that equals $1$ if $x$ and $y$ are both in $\Omega$, and zero otherwise.
I am interested in the following problem:

$$\min_{\Omega \in \mathcal{X}} \int_{x\in \mathbb{R}^2}\int_{y\in \mathbb{R}^2} d(x,y,\Omega) \mathbb{1}(x,y,\Omega) dx dy$$



*

*Is there a well-known solution to this problem?

*Is it possible to show that a solution to this problem must be convex?

*Is the problem even well-defined (e.g., does $d(x,y,\Omega)$ always exist? I tried to help guaranteeing it by including only closed and connected sets in $\mathcal{X}$, but I am not 100% sure that's enough)?

*Is the problem guaranteed to have a solution without more regularity conditions?


Some related problems are described here: https://link-springer-com.dist.lib.usu.edu/content/pdf/10.1007%2Fs10958-012-0717-3.pdf. Unfortunately, these related problems either (1) focus on "networks" rather than "thick" sets, (2) focus right away on convex sets, or (3) do not fix the area and instead impose an additional "cost" to the minimization problem for increasing the area. 
 A: I don't know if this is well known (I guess it is), but here is an argument showing that a disc is the unique minimizer. Since the minimizer will end up being convex, we can consider straight-line distances instead of distances within the set.
Given such a set $\Omega$ with centre of mass $(0,0),$ consider the set $\Omega^{\mathrm{new}}$ defined by
$$(x,y)\in\Omega^{\mathrm{new}} \iff |y|\leq\frac12 \int 1(x,y',\Omega)dy'.$$
So in each vertical line segment $\Omega_x=\{y\mid (x,y)\in\Omega\},$ we replace the set by a centred version.
By Fubini's theorem, the areas of $\Omega^{\mathrm{new}}$ and $\Omega$ are equal. For a fixed pair of distinct $x$-coordinates $x,x'$ consider the effect on the distance between points on those lines,
$$w(\Omega_x,\Omega_{x'})=\int_{y\in \Omega_x} \int_{y'\in \Omega_{x'}} |(x,y)-(x',y')| dy dy'$$
I claim this decreases when replacing $\Omega$ by $\Omega^{\mathrm{new}},$ except in the case where $\Omega$ and $\Omega^{\mathrm{new}}$ are a.e. equal or one is a null set.
To see this, first consider replacing $\Omega_{x'}$ by the interval (not necessarily centred at zero) of the same length but of the form $\{y\mid \Phi(y)<c\},$ where $\Phi(y)=\int_{y'\in \Omega_{x'}} |(x,y)-(x',y')| dy'.$ Note $\Phi$ is strictly convex so this does indeed produce an interval. This operation does not increase $w(\Omega_x,\Omega_{x'}).$ We then do the same the other way round. Both slices are now intervals, and a vertical translation takes them to $\Omega^{\mathrm{new}}.$ Doing this a bit more carefully would show that the weight strictly decreases for sets that actually change.
This shows that a minimizer must have reflection symmetry through every line through its centre of mass, and must therefore be rotationally invariant, and a bit of thought shows it has to be a disk.
The existence of a minimizer can be shown by compactness. The set of measures $\mu$ of total measure $a,$ contained in some fixed large disk, subject to the additional restriction $\mu(S)\leq\lambda(S)$ (bounded by Lebesgue measure), is compact in the vague (weak-*) topology. (I'm using measures here, but you can use something like $L^2$ if that's more intuitive.) Any indicator function of a compact set with area $a$ can be put into such a space, then we can invoke compactness to pick a minimizer $\mu.$
I will now argue that the minimizer has to be a set. (Actually this isn't really needed for the above argument - you could just replace sets by functions everywhere.)
 The potential $\Phi_\mu(x)=\int_{x'\in\mathbb R^2} |x-x'|d\mu(x')$ is convex, and in particular continuous. If there were sets $S_1,S_2$ with $\max \Phi_\mu(S_1)<\min \Phi_\mu(S_2)$ and $\mu(S_1)<\lambda(S_1)$ then we could move some mass from $S_2$ to $S_1$ and get a smaller total potential, a contradiction. This means $\mu$ is actually uniform distribution on a set of the form $\{\Phi_\mu<c\}.$
It is possible to have a compact connected closed set that is regular (equal to the closure of its interior) and not path-connected. Just thicken the sine part of the topologist's sine curve, and stick two of them back to back (e.g. the closure of the set of $(x,y)$ with $0<|x|<1$ and $|y-\sin 1/x|<\tfrac 1 4$), then scale to get the area correct. So the integral may not be finite. You can modify this example to get path-connected but a diverging integral.
A: Not a positive answer, but from here we find that the average distance of two random points in a $w \times h$ rectangle is:
$$
\frac1{15}
\left(
\frac{w^3}{h^2}+\frac{h^3}{w^2}+d
\left(
3-\frac{w^2}{h^2}-\frac{h^2}{w^2}
\right)
+\frac52
\left(
\frac{h^2}{w}\log\frac{w+d}{h}+\frac{w^2}{h}\log\frac{h+d}{w}
\right)
\right)\;,
$$
where $d=\sqrt{w^2+h^2}$.
Now WLOG let $a = 1$, $w = x$ and $h = 1/x$. We get the formula:
$$
\frac1{15}
\left(
x^5+x^{-5}+d
\left(
3-x^4-x^{-4}
\right)
+\frac52
\left(
x^{-3}\log(x^2+xd)+x^3\log(x^{-2}+d/x)
\right)
\right)\;,
$$
where $d=\sqrt{x^2+x^{-2}}$.
For $x > 0$ this has a minimum at $x = 1$, ruling out "long thin rectangles".  
A: This is by no means a full answer, but points out some structure.
I assume that $\Omega$ is convex and regular enough.
Let us denote
$$
F(\Omega)
=
\int_{x\in\Omega}\int_{y\in\Omega}|x-y|\,dy\,dx.
$$
In the convex setting, you are minimizing $F$.
Let $\ell(x,v)$ denote the distance from $x\in\Omega$ to $\partial\Omega$ in the direction $v\in S^1$.
We use the convention that if $x\in\partial\Omega$ and $v$ points inward, then $\ell(x,v)>0$.
More analytically, we can define $\ell(x,v)=\sup\{t\geq0;x+tv\in\Omega\}$.
The inner integral defining $F(\Omega)$ can be written in polar coordinates.
The radial integral is easy to compute, and we find
$$
F(\Omega)
=
\frac13\int_{x\in\Omega}\int_{v\in S^1}\ell(x,v)^3\,dv\,dx.
$$
This is an integral over the ring bundle $S\Omega=\Omega\times S^1$ (also called the sphere bundle — $S^1$ is just the 1D sphere).
One can change integration from the sphere bundle an integral over all lines and the space of lines using the so-called Santaló formula (see e.g. proposition 8.2 in these lecture notes for a proof in $\mathbb R^2$).
This leads to
$$
F(\Omega)
=
\frac1{3}
\int_{x\in\partial\Omega}
\int_{v\in S^1}
|\langle v,\nu_x\rangle|
\int_0^{\ell(x,v)}\ell(x+tv,v)^3
\,dt\,dv\,dx.
$$
Here $\nu_x$ is the unit normal at $x$.
The innermost integral is again an explicit 1D integral, and we get
$$
F(\Omega)
=
\frac1{12}
\int_{x\in\partial\Omega}
\int_{v\in S^1_{x,in}}
|\langle v,\nu_x\rangle|
\ell(x,v)^4
\,dv\,dx.
$$
On the other hand, the area of $\Omega$ is
$$
|\Omega|
=
\frac1{2\pi}
\int_{x\in\partial\Omega}
\int_{v\in S^1_{x,in}}
|\langle v,\nu_x\rangle|
\ell(x,v)
\,dv\,dx.
$$
See exercise 96 in the linked notes.
Here $S^1_{x,in}$ is the set of $v\in S^1$ which point towards the interior of $\Omega$ from $x\in\partial\Omega$.
If $\partial_{in}S\Omega$ denotes the inward pointing boundary of the sphere bundle (the whole boundary is $\partial S\Omega=\partial\Omega\times S^1$) and we denote by $\sigma$ the measure corresponding to $|\langle v,\nu_x\rangle|\,dv\,dx$, we have
$$
F(\Omega)
=
\frac1{12}
\int_{\partial_{in}S\Omega}\ell(x,v)^4\,d\sigma(x,v)
$$
and
$$
|\Omega|
=
\frac1{2\pi}
\int_{\partial_{in}S\Omega}\ell(x,v)\,d\sigma(x,v).
$$
This puts $F(\Omega)$ and $|\Omega|$ in a very similar neat form.
The suspected extremal case is a disc of some radius $R>0$.
In this case $\ell(x,v)=2R|\langle v,\nu_x\rangle|$ when $v$ points inward.
Therefore it might be convenient to use the measure $\tilde\sigma$ corresponding to $dv\,dx$ instead.
Let us denote $u(x,v)=|\langle v,\nu_x\rangle|$ for brevity.
Let me denote the length of the perimeter by $P=|\partial\Omega|$ and the diameter of $\Omega$ by $D$.
Denoting $L^p=L^p(\partial_{in}S\Omega,\tilde\sigma)$, we have (for any $\alpha\geq0$ and $p\geq1$)
$$
\begin{split}
F(\Omega)&=\frac1{12}\|\ell^4u\|_{L^1},\\
|\Omega|&=\frac1{2\pi}\|\ell u\|_{L^1},\\
P&=\frac1\pi\|1\|_{L^p},\\
D&=\|\ell u^\alpha\|_{L^\infty}.
\end{split}
$$
These with the isoperimetric inequality $4\pi A\leq P^2$ and the isodiametric inequality $4A\leq\pi D^2$ and Hölder's inequality give something to play with.
To give anything sharp, one has to apply Hölder's inequality in a way that gives equality when $\ell$ is a constant multiple of $u$.
