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$.
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.