is this sort of maximization problem understood?
Take a function $f:\mathbb{R}_+^n\to \mathbb{R}_+$ increasing in each dimension, differentiable, bounded above by some constant $C$. (possibly not concave). Fix $v\in\mathbb{R}^n_+$. \begin{align} \text{maximize }\ & \int f(x)\mu(dx) \\ \text{over }\ & \mu \text{, a distribution on } \mathbb{R}_+^n \\ \text{subject to }\ & E_\mu[x]=v. \end{align}
Here is a finite point mass result that uses Caratheodory's theorem.
Define $$A = \{(x, f(x)) \in \mathbb{R}^{n+1} : (x_1, ..., x_n) \geq 0\}$$ Let $X$ be any nonnegative random vector and assume $E[X]$ and $E[f(X)]$ are both finite. Then $(X,f(X)) \in A$ always, and so, $$ E[(X, f(X))] \in Conv(A) $$ where $Conv(A)$ is the convex hull of the set $A$.
Fix $v=(v_1, ..., v_n)$ as a vector with positive entries. We want to maximize $E[f(X)]$ subject to $E[X]=v$. This is equivalent to \begin{align} \mbox{Maximize:} \quad & y \quad (Eq. 1) \\ \mbox{Subject to:} \quad & x=v \quad (Eq. 2)\\ & (x,y) \in Conv(A) \quad (Eq. 3) \end{align} For simplicity assume a solution $(x^*,y^*)$ exists. Since $(x^*, y^*) \in Conv(A)$, and $A$ has $n+1$ dimensions, Caratheodory's theorem ensures $(x^*,y^*)$ can be written as the convex combination of at most $n+2$ points in $A$. However, since $(x^*, y^*)$ is on the boundary of $Conv(A)$, a simple extension to Caratheodory ensures $(x^*,y^*)$ can be written as a convex combination of at most $n+1$ points in $A$, so there are probabilities $p_1, ..., p_{n+1}$ that sum to 1, and nonnegative vectors $x_1, ..., x_{n+1} \in \mathbb{R}^{n}$, such that: $$(x^*, y^*) = \sum_{i=1}^{n+1} p_i(x_i, f(x_i)) $$ That is, we define an optimal random variable $X^*$ to take the value $x_i$ with probability $p_i$, for $i \in \{1, ..., n+1\}$, and this solves the problem. $\Box$