Firstly, you need also to assume that $\int u \mathrm{d}V_g$ and $\int k e^{4u}\mathrm{d} V_g$ are $C^1$ functionals on your function space.
Then: let $\mathcal{H}$ be a real Hilbert space. Let $S, K_1, K_2$ be $C^1$ functionals such that $K_1$ and $K_2$ intersect transversely, then the set $\mathcal{N} = \{K_1 = k_1\}\cap\{K_2 = k_2\}$ is a $C^1$ Hilbert manifold, and we can apply the theory of Lagrange Multiplier as follows: if $u\in \mathcal{N}$ is a constrained critical point of $S$, we must have $DS(u) = \lambda_1 DK_1(u) + \lambda_2 DK_2(u)$ where $D$ is the Frechet derivative on $\mathcal{H}$.
Applying to the calculus of variations, you have that Frechet derivatives $DS(u):\mathcal{H}\to\mathbb{R}$ satisfy the Gateaux differentiation condition $ DS(u) \cdot h = \lim_{t\to 0} \frac{1}{t}[ S(u+th) - S(u)] $ which using calculus of variations and applying to your case where $\mathcal{H} = L^2(M,\mathrm{d}V_g)$ with $\langle f,g\rangle = \int fg \mathrm{d}V_g$, and $S$ is the energy functional, $K_1$ and $K_2$ are the given constraints gives $ DS[u] \cdot h = \int h P(u) + u P(h) \mathrm{d}V_g $ where using that $P$ is self-adjoint gives you $ = \int 2h P(u) \mathrm{d}V_g ~.$ Similarly you get $ DK_1[u] \cdot h = \int h \mathrm{d}V_g $ and $ DK_2[u] \cdot h = \int 4 k e^{4u} h \mathrm{d}V_g~. $
Putting this together tells you that for a constrained critical point of $S$ under constraints $K_1,K_2$, the solution $u$ must satisfies that for all $h\in\mathcal{H}$
$ \int h\left( 2P(u) - \lambda_1 - 4\lambda_2k e^{4u}\right) \mathrm{d}V_g = 0$
and hence the term inside the parentheses must be 0.
For the case in the paper you are reading, some further justification probably is involved. (For example, you hadn't told us on what function space is $P$ a self-adjoint operator).) If it is $L^2$ as I inferred, $\int ke^{4u}\mathrm{d}V_g$ is not actually $C^1$ (it is not even bounded). Nevertheless, one can usually apply some renormalisation to replace $K_2$ by another functional that is identical to $K_2$ when $|K_2[u]|\leq \epsilon$ and suitably cut-off outside. Slightly more worrisome is that $P$ is likely to be only densely defined on $L^2$ for most geometric applications, which throws some more caveats into the above.
But on a formal level, the above derivation (and more simply speaking, Rahul's answer) explains what you seek.