Depending on your integrals (dimension, smoothness, etc), consider approximating with a Gaussian quadrature, eg as in
and from then on you effectively have a sum of function values with (potentially) parameter-dependent weights, in any case it should be easy to handle with any nonlinear solver.
If this is about structural estimation, I recently asked a similar question
and got some good suggestions, which I am working through. Currently I find NLopt (Augmented Lagrangian for the constraint) very robust, but that’s on a toy problem.