Here is the lower bound of $0.49$ for all $\alpha\ge 1$. Note that $\min_\pi F(\pi)$ is a non-decreasing function of $\alpha$, so it is enough to consider $\alpha=1$. Also, the truth is about $0.55$ for $\alpha=1$ and $0.62$ for $\alpha=2$, so, if you care about entropy, the uniform distribution that gives $\frac 23$ for all $\alpha$ is not too far from the truth.
If you do it numerically and try just to consider functions constant on short intervals, then do not forget the self-action on each interval in the sum to be optimized. It looks small compared to the whole for spread distributions but without it everything will quickly converge to a singleton. I made this error when trying to optimize numerically and it took me some time to realize and fix it.
Now the bound. I'll just consider $G(\pi)=\frac{F(\pi)}{2}=\iint \pi(x)^2\pi(y)|x-y|\,dx\,dy$. Note that this expression is invariant with respect to the scaling $\pi\mapsto c\pi(c\cdot)$, so the support restriction is not important. What matters is only that $\int\pi=1$.
Write the Lebesgue decomposition $\pi(x)=\int_{t>0}\chi_{E_t}(x)\,dt$ where $E_t$ is a decreasing nested family of measurable sets. Then
$$
G(\pi)=\iint_{s,t>0}2sQ(s,t)\,ds\,dt
$$
where
$$
Q(s,t)=\iint_{E_s\times E_t}|x-y|\,dx\,dy\,.
$$
It is not hard to see that if $0<\ell<L$ are lengths of $E_s,E_t$ (not necessarily in this order), then the least value of $Q$ is given by 2 concentric intervals, in which case it is
$$
\frac{\ell L^2}{4}+\frac{\ell^3}{12}\ge \frac{\ell L(\ell+L)}{6}\,.
$$
The RHS is just made up to get a symmetric expression in $s,t$, which leads to a nice optimization problem but costs about $12\%$ in the precision of the bound. Also it is clear from here that the only real competitors are symmetric unimodal distributions.
Now, putting $f(t)=|E_t|$, we see that we need to estimate
$$
\frac 16\iint 2s f(s)f(t)[f(s)+f(t)]\,ds\,dt=\frac 13\iint (s+t) f(s)^2f(t)\,ds\,dt\\
=\frac 13\int (s+a)f(s)^2\,ds
$$
from below under the conditions $\int f(s)\,ds=1,\int sf(s)\,ds=a>0$.
This is an almost pure Hilbert space question except for the positivity constraint. We'll just define $g_b(s)=\frac {(1-bs)_+}{a+s}$ for $b>0$ and use Cauchy-Schwarz:
$$
\left[\int (s+a)f(s)^2\,ds\right]\left[\int (s+a)g_b(s)^2\,ds\right]\ge
\left[\int (s+a)f(s)g_b(s)\,ds\right]^2
$$
The integral on the right is
$$
\int f(s)(1-bs)_+\ge \int f(s)(1-bs)=1-ab\,,
$$
so we get the bound
$$
\int (s+a)f(s)^2\,ds\ge (1-ab)^2\left[\int (s+a)g_b(s)^2\,ds\right]^{-1}
$$
as long as $ab<1$. Choosing $ab=\frac 14$, we get the right hand side independent of $a$ and equal approximately $0.735$ (the best choice is the solution of some transcendental equation, but we have lost $12\%$ already, so who cares now) giving $G(\pi)\ge 0.245$, as promised.
Simple functions that give values of $F(\pi)$ close to the true minimum (on $[-0.5,0.5]$ instead of $[0,1]$) are (properly normalized) $(1-2x)^2$ and $1-x$ for $\alpha=1$ and $\alpha=2$ respectively. Those are not exact optimizers but you may find it hard to beat them when doing numeric optimization.