Given a bounded $\Omega\subset R^d$, $N$ functions $F_1,\dots,F_N\in L^\infty(\Omega)$, and a positive integrable $\omega\in L^1_+(\Omega)$, define the following function from $R^N$ to $R$: $$ f(\alpha)=\int_\Omega \max\limits_i \{F_i(x)+\alpha_i\}\omega(x) dx. $$ Note that $f$ is a convex function.

Question:I would like to show that the subdifferential \begin{equation} \partial f(\alpha)=\Bigg\{P=\int \theta(x)\omega(x)d x:\qquad \theta_i(x)\geq 0,\,\sum_{i=1}^N\theta_i(x)=1 \,\text{a.e.},\qquad \text{and } \theta_i(x)=0\text{ a.e. in }[F_i(x)+\alpha_i<\max_j(F_j(x)+\alpha_j)] \Bigg\}\qquad\qquad(1) \end{equation}

This may be elementary, but my convex analysis is quite rusty...

Here is how I guessed this result: fix any $x\in\Omega$, and look at the pointwise function of $\alpha$ above $x$ $$ f_x(\alpha)=\max\limits_i \{F_i(x)+\alpha_i\}. $$ By standard properties of the $\max$ function the subdifferential is classically obtained as the convex hull of all the active gradients, i-e $$ \partial f_x(\alpha)=Co\Big(\nabla_{\alpha}(F_i(x)+\alpha_i):\qquad i\in\{1...N\}\text{ s.t. } F_i(x)+\alpha_i=\max(F_j(x)+\alpha_j)\Big). $$ Since the functions in the max are here linear in $\alpha$, we immediately get that those can be identified as convex coefficients with zeros on the inactive components, i-e $$ \partial f_x(\alpha)=\Big\{\theta(x)=(\theta_1(x)...\theta_N(x)): \qquad\sum\limits_{i=1}^N\theta_i(x)=1,\theta_i(x)\geq 0,\text{ and }\theta_i(x)=0 \text{ if } F_i(x)+\alpha_i<\max(F_j(x)+\alpha_j)\Big\}. $$ Then by linearity of the subdifferential and because $f$ is the infinitesimal sum (i-e the integral) $f(\alpha)=\int f_x(\alpha)\omega(x)dx$ we expect that $$ P\in\partial f(\alpha) \quad\Leftrightarrow\quad P=\int p(x)\omega(x)dx\text{ for some vector-valued function }p(x)\in \partial f_x(\alpha), $$ which simply means that $P$ is of the form (1). However, I am facing the usual issue in subdifferential calculus: it is easy to check that any $P=\int p(x)dx $ as above belongs to $\partial f(\alpha)$, but are there any other subgradients of a different form? I claim that the answer is negative (by analogy of integrals with finite sums), but right now I'm stuck. In order words, how can I carefully justify the above "linearity" argument? Is it always true that if $g(\alpha)=\int g_x(\alpha)dx$ then its subdifferential can be recoverd from all the integrated subgradients $P=\int p_xdx\in \partial g(\alpha)$, where $p_x\in \partial g_x(\alpha)$ is any subgradient selection in the "frozen" fiber above $x$?

side-note on the context:for my research I need some kind of generalized vector-valued bathtub principle. More precisely, I am trying to characterize the minimizers of $s=(s_1,...s_N)\in (L^1(\Omega))^N\mapsto \mathcal F(s)=\sum\limits_{i=1}^N \int F_i(x)s_i(x)dx$ with the constraints that $s_i(x)\geq 0$, fixed mass $\int s_i(x)dx=m_i$, and a total saturation $\sum s_i(x)=\omega(x)$ a.e. for a given function $\omega\in L^1_+$. The above problem arose when I tried brute-force, i-e writing the Lagrange multipliers associated to the constraints and then apply a min-max argument. I couldn't find any reference on any kind of multicomponent bathtub principle so if anyone has any suggestion I would greatly appreciate.