Method of Lagrange multipliers with random variables I'll illustrate the issue I'm having with a simple problem.

Let $c_1, c_2 \in \mathbb{R}$, and $Z$ a real-valued random variable. Let $u:\mathbb{R} \rightarrow \mathbb{R} $ be a differentiable function, and $f(c_1,Z)$ be a real-valued function that is differentiable with respect to $c_1$.

The problem is:
Maximise: $u(c_1) + \mathbb{E}[u(c_2)]$ 
such that $c_2=f(c_1, Z)$

This problem is easily solved by direct substitution, and the answer is

$u'(c_1) + \mathbb{E}\left[u'(f(c_1,Z)) \frac{\partial f(c_1,Z)}{\partial c_1}\right] = 0 \label{answer}$

The issue is how to write down a Lagrangian whose extrema correspond to the solution of this problem in the normal way.

My first instinct was to write
$\mathcal{L} = u(c_1) + \mathbb{E}[u(c_2)] + \lambda(c_2 - f(c_1,Z))$. 

However, this doesn't seem to make sense.
If you treat $c_2$ as if it is a random variable, then the derivative of  $\mathbb{E}[u(c_2)]$ with respect to $c_2$ gives zero, and this can't possibly give the right answer.
On the other hand, it doesn't make sense to treat it as non-stochastic either, since it is 'forced' to be stochastic by the constraint.

Question: How do I write a Lagrangian whose extrema correpond to the solution of the above optimisation problem?
If you treat $c_2$ as if it is a random variable, then the derivative of
  $E[u(c_2)]$ with respect to $c_2$ gives zero.


Why? Such an assertion doesn't follow from anywhere. The subtleties lie elsewhere. The problem with the Lagrangian considered by the OP 

$$\mathcal{L} = u(c_1) + \mathbb{E}[u(c_2)] + \lambda(c_2 - f(c_1,Z))$$

is that it also is a random variable since now $Z$ appears outside the expected value (The direct substitution approach followed by maximization with respect to $c_1$ only, does not create any such issues).

Now, do we /  can we maximize a random variable? Well, no, because the essential characteristic of a random variable is that it is a function whose value cannot be set by command and control.

But one could say "ok, let's pretend that this Lagrangian is not a random variable, and just write down the conditions for maximization, even though we know that we can't force the solution".

But this won't work: if one attempts to do it one will eventually obtain

$$u'(c_1) + \mathbb{E}\left[u'(f(c_1,Z))\right] \cdot  \frac{\partial f(c_1,Z)}{\partial c_1} = 0 $$

which is not the same as the condition obtained through direct substitution, because here the partial derivative is outside the expected value.

(for those who may think "hey, then how do we apply maximization procedures in the maximum likelihood approach" the answer is that there, nothing is random anymore when we get to apply the maximization steps).