I am not really answering your question, so don't worry about the bounty. Just some random thoughts. Don't have time to do more. I am using the Lebesgue measure.
If $\int f_{u}=1$ and $\int g_{u}=1$ then since the functions $s\in
(0,\infty)\mapsto s^{u}$ and $s\in(0,\infty)\mapsto s^{1-u}$ are increasing,
you have that $q_{0}^{u}q_{1}^{u}\leq f_{u}^{u}g_{u}^{1-u}$ and so $f_{u}$ and
$g_{u}$ are the unique solution.
If $\int f_{u}>1$ and $\int g_{u}=1$, then the solution $q_{u}$ cannot always
coincide with $f_{u}$ since $\int q_{0}=1$ and so $q_{0}<f_{u}$ in some set
$E\subseteq\Omega$ with positive measure. If $\int f_{l}=1$, then since
$f_{l}\leq q_{0}$ and $1=\int f_{l}=\int q_{0}$, you have $\int(q_{0}%
-f_{l})=0$ and so $q_{0}-f_{l}=0$ for $\mu$ a.e. $x\in\Omega$. Thus $f_{l}$
and $g_{u}$ are the unique solution.
If $\int f_{l}<1$, assume that there exists a solution $q_{0}$ and assume that
for some $n$ the set $G_{n}=\{x\in G:f_{l}(x)+\frac{1}{n}<q_{0}(x)<f_{u}%
(x)-\frac{1}{n}\}$ has positive measure. Take $\varphi=0$ outside $G_{n}$,
$|\varphi|\leq\frac{1}{2n}$ in $G_{n}$, and $\int_{G_{n}}\varphi=0$. Setting $F(q)=\int q^{u}g_{u}^{1-u}$, you have%
$$
g(t):=F(q_{0}+t\varphi)\geq g(0)=F(q_{u})
$$
for every $t$ small and so
$$
0=\frac{dg}{dt}(0)=\int uq_{0}^{u-1}g_{u}^{1-u}\varphi=u\int_{G_{n}}%
q_{0}^{u-1}g_{u}^{1-u}\varphi.
$$
In turn, given $\psi$ with $|\psi|\leq\frac{1}{4n(1+\mu(G_{0}))}$ and $\psi=0$
outside $G_{n}$, we have that $\varphi:=\psi-\int_{G_{n}}\psi$ satisfies
$|\varphi|\leq\frac{1}{2n}$ in $G_{n}$, and $\int_{G_{n}}\varphi=0$ and so by
Fubini's theorem$$
0=\int_{G_{n}}q_{0}^{u-1}g_{u}^{1-u}\left( \psi-\int_{G_{n}}\psi
\,d\mu\right) \,d\mu=\int_{G_{n}}\psi\left( q_{0}^{u-1}g_{u}^{1-u}
-\int_{G_{n}}q_{0}^{u-1}g_{u}^{1-u}\,d\mu\right) \,d\mu.
$$
Since this is true for every $\psi$ with $|\psi|\leq\frac{1}{4n(1+\mu(G_{n}%
))}$, it should imply that $q_{0}^{u-1}g_{u}^{1-u}-\int_{G_{n}}q_{0}
^{u-1}g_{u}^{1-u}\,d\mu=0$ for $\mu$ a.e. $x\in G_{n}$. So $q_{0}^{u-1}
g_{u}^{1-u}$ is constant in $G_{n}$, which implies that $q_{0}=c_{n}g_{u}$ in
$G_{n}$ for some $c_{n}>0$. But $G_{n}\subseteq G_{n+1}$ and so $c_{n}
=c_{n+1}$. Thus $q_{0}=cg_{u}$ in $\bigcup G_{n}=\{x\in G:f_{l}(x)<q_{0}
(x)<f_{u}(x)\}$.
This shows that if $\int g_{u}=1$ then either there is no solution $q_{0}$ or
there is and it has the form$$
q_{0}(x)=\left\{
\begin{array}
[c]{ll}
f_{u}(x) & x\in E_{1},\\
cg_{u}(x) & x\in E_{2},\\
f_{l}(x) & x\in E_{3}.
\end{array}
\right.
$$
I guess that taking simple choices of functions, such as, $f_{l}=g_{l}=0$,
$g_{u}=1$, $f_{u}=2$ and $\Omega=[0,1]$, one could try to see if there is
uniqueness or not.
I don't have time to do the interesting case $\int g_{u}>1$, $\int f_{u}>1$, $\int g_{l}<1$,$\int f_{l}<1$.
Edit
The solution always exist. provided
$$
\left( s,t\right) \mapsto s^{u}t^{1-u}=f(s,t)
$$
is concave.
So if
$$
\ell=\sup\int q_{1}^{u}q_{0}^{1-u}$$
and you construct a sequence $(q_{1,n},q_{0,n})$ such that
$$
\ell=\lim_{n\rightarrow\infty}\int q_{1,n}^{u}q_{0,n}^{1-u}$$
then if $g_{u}$ and $f_{u}$ are bounded then you can find a subsequence
$(q_{1,n_{k}},q_{0,n_{k}})\overset{\ast}{\rightharpoonup}(q_{1},q_{0})$ in
$L^{\infty}$. Since $f$ is concave
$$
\ell=\limsup_{k\rightarrow\infty}\int q_{1,n_{k}}^{u}q_{0,n_{k}}^{1-u}\leq\int
q_{1}^{u}q_{0}^{1-u}.
$$
Again by weak star convergence, if $\Omega$ has finite measure $1=\int
q_{1,n_{k}}\rightarrow\int q_{1}$ and $1=\int q_{0,n_{k}}\rightarrow\int
q_{0}$. The bounds should also be satisfied. $\int_{E}f_{l}\leq\int_{E}%
q_{0,k}\rightarrow\int_{E}q_{0}$ so $\int_{E}f_{l}\leq\int_{E}q_{0}$ for every
measurable set $E$ so $f_{l}\leq q_{0}$ and the same is true for the other
constraints. So you do have a maximum.
