Bivariate optimal density Consider any feasible $p:[0,1]^2\to [0,1]$ that allows discontinuities and the problem
$$\min_{p(.)} \int_0^1\int_0^1 p(x,y)^2 dF(x) dG(y)$$
s.t.  $$\int_0^1 p(x,y)dG(y)=k\phantom{0} for \phantom{0} a.e.\phantom{0} x\in [0,1]$$
where $F,G$ are distribution functions and $1/2<k<1$.
How could I solve the problem analytically?

By your condition on $p(x,y)$, we see that $p(x,y)=p(y)$ a.e.. Hence you just need to find
$$ \min_p \int_0^1 p(y)^2 dG(y) $$
s.t.
$$ \int_0^1 p(y) dG(y) =k $$
To minimize the functional, we perform a variation according to Lagrange multipliers. Take $p \to p + h \xi$ then we see first order terms in $h$ correspond to
$$ \delta^{(1)} =\int_0^1 [2p(y)  - \lambda ]\xi(y) dG(y)$$
To get a critical point, we need $p(y)=\lambda/2$ a.e. since we need $\delta^{(1)} =0 $ for all test functions $\xi$. But the constraint gives us that the solution is
$$ p(x,y) = k \quad a.e. $$
For reference : Calculus of variations: Lagrange multipliers
