Shape with area $a$ and shortest average distance between any two points
Solution 1:
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.
Solution 2:
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".
Solution 3:
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$.