Let $Z$ be a random variable with finite mean $\mathbb E[Z]$ and let $\phi:(-\infty,\infty) \rightarrow (-\infty,\infty]$ be a convex l.s.c function which is differentiable at $z=1$ with gradient $\phi'(1)=0$.
Question
How to find / characterize all the numbers $c\in \mathbb R$ verying the "stochastic" variational inequality $$ \begin{split} &1 \in \mathbb E[\partial \phi^*(Z-c)],\\ &\phi'(1) = 0. \end{split} $$
N.B.: Here $\phi^*(z) := \sup_t zt - \phi(t)$ defines the usual convex conjugate and $$\partial \phi^*(z) := \{g \in \mathbb R | \phi^*(t) \ge \phi^*(z) + g(z-t)\;\forall t \in \mathbb R\} $$ is the subdifferential of $\phi^*$ at $z$.
Some concrete examples
Z takes on finitely many values. Suppose $Z$ can only take one of $n$ distinct values $z_1,\ldots,z_n$ with probabilities $p_1,\ldots,p_n$.
Then the above problem reads:
$$ \begin{split} &1 \in \sum_{i=1}^np_i\partial \phi^*(z_i-c),\\ &\phi'(1) = 0, \end{split} $$ which is just a classical variational inclusion problem. The particular case $n=1$ gives
$$ \begin{split} &1 \in \partial \phi^*(z_1-c) \iff z_1-c \in \partial \phi(1),\\ &\phi'(1) = 0, \end{split} $$ giving the unique solution $c=z_1$.
Quadratic potential. Suppose $\phi(z) = \alpha(z-1)^2/2$, for some $\alpha > 0$. A simple computation gives, $\phi^*(z) = \dfrac{z^2}{2\alpha} + z$, from which it follows that the unique solution for our problem is $c=\mathbb E[Z]$.