Convex optimization problem outputs a unique solution if $f_l\leq q_o$ but not if $f_l\leq q_o\leq f_u$, why? I have the following convex optimization problem ($q_0$ and $q_1$ are unkown):
$$\begin{array}{ll} \text{maximize} & \displaystyle\int_{\Omega} q_1^u{q_0}^{1-u}\mathrm{d}\mu\\ \text{subject to} & \displaystyle\int_{\Omega} q_0 = 1,\quad \displaystyle\int_{\Omega} q_1 =1 \\  & f_l \leq {q_0} \leq f_u\\ & g_l \leq q_1 \leq g_u\end{array}$$
where $u\in(0,1) $ and
$$\int_{\Omega}f_l \mathrm{d}\mu\leq 1,\quad\int_{\Omega}g_l \mathrm{d}\mu\leq 1$$
$$\int_{\Omega}f_u \mathrm{d}\mu\geq 1,\quad\int_{\Omega}g_u \mathrm{d}\mu\geq 1$$
Here, $q_0,q_1$ are density functions on $\Omega$, and $f_l,f_u,g_l,g_u$ are some known positive functions on $\Omega$.

I try to solve this optimization problem and I get exactly the same solution for every $u$, if I set $f_u=\infty$ and $g_u=\infty$, i.e., if there are only lower bounds for $q_0$ and $q_1$. If I also put $g_u$ and $f_u$ as some functions which are bounded above by $B<\infty$, the solution is not independent of $u$, I get different solutions for every $u$. Can one show this using math?

P.S: In the examples I had, $g_u$ and $f_u$ were also integrable over $\Omega$.
Addendum: I tried with $f_l=0$ and $g_l=0$ and for this case, namely with only upper bounding functions, I also get results independent of $u$, namely the same result for every $u$.
 A: 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.
