Computing expectation of nonlinear function of multinomial logit

Suppose I have constants \{\delta_j\}_{j \in J} and have random variables \{\varepsilon_j\}_{j\in J} distributed iid according to type I extreme value distribution as in a standard multinomial logit model. Let F: \mathbb{R}^J \rightarrow \mathbb{R} be the joint distribution \varepsilon_j's. Suppose further that I have a concave function u: \mathbb{R} \rightarrow \mathbb{R}.

I want to compute the expectation \int u \left( \max{\{\delta_j + \varepsilon_j \}_{j \in J}} \right) dF(\varepsilon).

Further, I want to make the \delta_j's be functions of some other vector of parameters x, and then maximize the function H(x) = \int u \left( \max{\{\delta_j(x) + \varepsilon_j \}_{j \in J}} \right) dF(\varepsilon) over x.

Is there a recommended method for doing this? To compute the integral I could use monte carlo simulation I guess. Toward that end, what is the best way to sample from type I extreme value?

I want to make sure that whatever way I use to compute the integral will make it easy to optimize the integral over the parameters as described in my “Further” sentence above. Any recommendations on this?