I wish to extremise
$$ Q = \int_0^h u \, \,dy $$
with the following constraints
$$B = \int_0^h u g \, \, dy \\ M = \int_0^h u^2 +\left(\int_0^y g \, \,dy^*\right) \, \, dy$$
where $M,B$ are constant. So my idea is to set this up as a variational problem as such
$$\varepsilon \left(Q + \lambda B + \mu M\right) = 0$$
We assume
$$u = u_0 + \varepsilon u_1 \\ g = g_0 + \varepsilon g_1 $$
I also wish to constrain that $0 < g \leq g^*$ for some given $g^*$, but I will ignore that for the time being. Taking $O(\varepsilon)$ terms gives
$$\varepsilon \left[\int_0^h u_1 + \lambda(u_0g_1+u_1g_0) + \mu 2 u_0u_1 + \mu\left[y \int_0^h g_1 \right]_0^h - \mu \int_0^h y g_1 \, \,dy \right] = 0$$
Where I have used integration by parts on the $\int g$ term in $M$. This gives
$$ \int_0^h \left[ u_1(1 + \lambda g_0 + 2 \mu u_0) + g_1(\lambda u_0 + \mu(h-y) \right] = 0 $$
Taking the two terms to be independent gives
$$u_0 = -\frac{\mu}{\lambda}(h-y) \\ g_0 = -\frac{1}{\lambda} + 2\left(\frac{\mu}{\lambda}\right)^2(h-y)$$
All good so far I think? I expected $u_0,g_0$ to be linear. For ease of notation set $\omega = -\mu / \lambda$ and $-1/\lambda = \beta$.
Then we can arrive substituiting $u_0, g_0$ into the definitions of $M,B$ with some algebra, at
$$M = \omega^2 h^3 + \frac{\beta h^2}{2} \\ B = \frac{2}{3}\omega^3 h^3 + \frac{\beta \omega h^2}{2}$$
And hence, these two simultaneous equations give
$$Q = \frac{(3(M \omega - B)^{2/3}}{2 \omega} $$
So here is where I run out of steam a bit. I guess I want to extremise $Q$ so take
$$\frac{\partial Q}{\partial \omega} = 0 \Rightarrow \omega = \frac{3B}{M}$$
However, this $\omega$ implies that $\beta = -2M / h^2 < 0$, which means that the constant term in $g_0$ is negative, which is unphysical in my set up, so we must have $\beta = 0$ at the maximal value.
Proceeding with $\beta = 0$ gives us
$$\omega = \frac{3B}{2M}$$
which gives me a maximal
$$Q = \frac{M}{2^{2/3} 3^{1/3} B^{1/3}}$$.
However, I am not sure if this answer/my method is correct, nor how to incorporate the inequality constraint that $g_0$ is smaller than some fixed value.
I can get a measurement for maximal $Q$ via a different method, which provides me with the correct powers/dimension, but the prefactors (in particular the $2^{2/3}$) are incorrect.
Any pointers on how to answer this type of question? How to incorporate the inequality? Or if indeed my method is correct.
EDIT
To show why I select $\beta = 0$ as the maximal case. Consider the plot of $Q$ against $\omega$. If we input some reasonable (physically sound) values for $M = 4.5$ and $B = 2$ we have
So we can see that $Q$ has a maximal value. However, plotting $\beta$ against $\omega$ for the same values of $M,B$ shows us that $\beta$ gives a negative value at this maximal value, which is unphsyical. So the maximal $Q$ goes with $\beta = 0$.

