I tried to work on the problem and I think I am able to resolve some points. Here is my work:

Consider the Lagrangians: $$L_0=\int_{\mathbb{R}} g^u{f}^{1-u}\mathrm{d}\mu+\int_{\mathbb{R}}\left(\lambda_0(f-f_L)+\lambda_{00}(f_U-f)\right)\mathrm{d}\mu+\mu_0\left(\int_{\mathbb{R}} f\mathrm{d}\mu-1\right)$$

$$L_1=\int_{\mathbb{R}} g^u{f}^{1-u}\mathrm{d}\mu+\int_{\mathbb{R}}\left(\lambda_1(g-g_L)+\lambda_{11}(g_U-g)\right)\mathrm{d}\mu+\mu_1\left(\int_{\mathbb{R}} g\mathrm{d}\mu-1\right)$$

Taking the Gateux derivatives of the Lagrangians, at the direction of funtions $\psi_0$ and $\psi_1$, respectively, leads to

$$\frac{\partial L_0}{\partial f}=\int\left((1-u)\left(\frac{g}{f}\right)^u+\lambda_0-\lambda_{00}+\mu_0\right)\psi_0\mathrm{d}\mu=0\quad\quad(1)$$

$$\frac{\partial L_1}{\partial g}=\int\left(u\left(\frac{g}{f}\right)^{u-1}+\lambda_1-\lambda_{11}+\mu_1\right)\psi_1\mathrm{d}\mu=0\quad\quad\quad\,\,\,(2)$$

Here according to Gateux derivative, $\psi_0$ and $\psi_1$ are arbitrary functions. I take them as integrable functions with $\int \psi_0 \mathrm{d}\mu=1$ and $\int \psi_1 \mathrm{d}\mu=1$

There are actually $3$ cases for each Lagrangian $L_0$ and $L_1$. For $L_0$ we have $$f=f_L, \quad f=f_u, \quad f_L<f<f_U$$ and for $L_1$ we have $$g=g_L, \quad g=g_u, \quad g_L<g<g_U$$

The conditions above $\partial L_0/\partial f=0$ and $\partial L_1/\partial g=0$ make sense only for the conditions $f_L<f<f_U$ and $g_L<g<g_U$.

Hence, I can write the maximizing functions as

$$f(y)=\begin{cases}f_L\quad y\in E_0\\ f_U\quad y\in E_1\\h_0\quad y\in E_2\\\end{cases}\quad g(y)=\begin{cases}g_L\quad y\in E_3\\ g_U\quad y\in E_4\\h_1\quad y\in E_5\\\end{cases}$$

Based on this result I have following conclusions:

$1.$ For the general case, namely if $g_L\neq \infty$ and $g_U\neq \infty$, the sets $E_k$ and functions $h_0$ and $h_1$ are necessarily dependent on $u$, although there may be special cases. This is because there is no simplification which can show that $E_k$ and functions $h_0$ and $h_1$ may be written independent of $u$.

 

$2.$ if $g_L= \infty$ and $g_U= \infty$, the results are independent of $u$ and I can confirm this with experiments. In order to show this analytically, first, I should be able to show that $g/f$ is constant in (1) and (2).

 

I am trying to explain this: $\lambda_{00}$ and $\lambda_{11}$ are zero, $\lambda_{0}$ and $\lambda_{1}$ are some positive functions. $\mu_0$ and $\mu_1$ must be negative constants, if not the equations (1) and (2) above wont be $0$. My problem is that $(1−u)(g/f)^u+\lambda_0$ can also be a constant function which is equal to $-\mu_0$, if $\lambda_0$ is carefully chosen. Then, the integral equation will hold while $g/f$ is not a constant. However, this is never the case. Am I wrong?

Please feel free to comment and post a new anwer based on mine. Because still there are missing points.

I tried to work on the problem and I think I am able to resolve some points. Here is my work:

Consider the Lagrangians: $$L_0=\int_{\mathbb{R}} g^u{f}^{1-u}\mathrm{d}\mu+\int_{\mathbb{R}}\left(\lambda_0(f-f_L)+\lambda_{00}(f_U-f)\right)\mathrm{d}\mu+\mu_0\left(\int_{\mathbb{R}} f\mathrm{d}\mu-1\right)$$

$$L_1=\int_{\mathbb{R}} g^u{f}^{1-u}\mathrm{d}\mu+\int_{\mathbb{R}}\left(\lambda_1(g-g_L)+\lambda_{11}(g_U-g)\right)\mathrm{d}\mu+\mu_1\left(\int_{\mathbb{R}} g\mathrm{d}\mu-1\right)$$

Taking the Gateux derivatives of the Lagrangians, at the direction of funtions $\psi_0$ and $\psi_1$, respectively, leads to

$$\frac{\partial L_0}{\partial f}=\int\left((1-u)\left(\frac{g}{f}\right)^u+\lambda_0-\lambda_{00}+\mu_0\right)\psi_0\mathrm{d}\mu=0\quad\quad(1)$$

$$\frac{\partial L_1}{\partial g}=\int\left(u\left(\frac{g}{f}\right)^{u-1}+\lambda_1-\lambda_{11}+\mu_1\right)\psi_1\mathrm{d}\mu=0\quad\quad\quad\,\,\,(2)$$

Here according to Gateux derivative, $\psi_0$ and $\psi_1$ are arbitrary functions. I take them as integrable functions with $\int \psi_0 \mathrm{d}\mu=1$ and $\int \psi_1 \mathrm{d}\mu=1$

There are actually $3$ cases for each Lagrangian $L_0$ and $L_1$. For $L_0$ we have $$f=f_L, \quad f=f_u, \quad f_L<f<f_U$$ and for $L_1$ we have $$g=g_L, \quad g=g_u, \quad g_L<g<g_U$$

The conditions above $\partial L_0/\partial f=0$ and $\partial L_1/\partial g=0$ make sense only for the conditions $f_L<f<f_U$ and $g_L<g<g_U$.

Hence, I can write the maximizing functions as

$$f(y)=\begin{cases}f_L\quad y\in E_0\\ f_U\quad y\in E_1\\h_0\quad y\in E_2\\\end{cases}\quad g(y)=\begin{cases}g_L\quad y\in E_3\\ g_U\quad y\in E_4\\h_1\quad y\in E_5\\\end{cases}$$

Based on this result I have following conclusions:

$1.$ For the general case, namely if $g_L\neq \infty$ and $g_U\neq \infty$, the sets $E_k$ and functions $h_0$ and $h_1$ are necessarily dependent on $u$, although there may be special cases. This is because there is no simplification which can show that $E_k$ and functions $h_0$ and $h_1$ may be written independent of $u$.

$2.$ if $g_L= \infty$ and $g_U= \infty$, the results are independent of $u$ and I can confirm this with experiments. In order to show this analytically, first, I should be able to show that $g/f$ is constant in (1) and (2).

I am trying to explain this: $\lambda_{00}$ and $\lambda_{11}$ are zero, $\lambda_{0}$ and $\lambda_{1}$ are some positive functions. $\mu_0$ and $\mu_1$ must be negative constants, if not the equations (1) and (2) above wont be $0$. My problem is that $(1−u)(g/f)^u+\lambda_0$ can also be a constant function which is equal to $-\mu_0$, if $\lambda_0$ is carefully chosen. Then, the integral equation will hold while $g/f$ is not a constant. However, this is never the case. Am I wrong?

Please feel free to comment and post a new anwer based on mine. Because still there are missing points.