Simulation Optimization is the name of the field you are looking for. It is closely related to, and overlaps with. Machine Learning, such as training Deep Learning models.
In Simulation Optimization, there is no such thing as training on all the data, because there are an infinite number of data points which could be generated, not only at different values of the optimization variables, but the number of replications increasing to infinity.
There are a number of general approaches, including stochastic gradient descent (originally called Kiefer-Wolfowitz when it was invented 70 years ago, although that was the finite difference version, and only in one dimension), Gradient based methods, to include Quasi-Newton, or even Newton (simulating for Hessian).
Any algorithm which uses averages of sample for objective function and gradients (or Hessian) can be called Sample Average Approximation (SAA). The term is generic enough that it can refer to anything. But SAA is often used as a synonym for any algorithm using common random numbers, which is also called Sample Path Optimization, although that too can mean almost anything.
It is noteworthy, especially in conjunction with derivative based methods, but also for zeroth order (Objective function only) methods, use of common random numbers can be extremely important for performance. If calculating a finite difference derivative at two nearby points, if the errors are independent, the variances add, and get divided by a small number, which creates a huge variance in the finite difference derivative; that's not a good thing. However, if common random numbers are used, and particularity if there is good synchronization in the implementation so that very high positive correlation is achieved at nearby points, the individual data point errors might cancel out to a large extent, and finite difference derivatives could potentially be more accurate than evaluation of the objective at an individual point Similarly for differencing gradients from adjacent iterations in order to feed a Quasi-Newton Hessian update - this concept seems foreign to many Machine Learning people.
Surrogate models (sometime known as Response Surface Methodology) fit a model to objective values at sampled points, and then optimize the surrogate function; then additional points ere sampled to improve the model, then re-optimized, etc. One manifestation is fitting a Gaussian Process, then optimizing that. (a form of Bayesian Optimization) Stochastic Kriging is another term for essentially the same thing. The uncertainty in the fitted model expands between sampled data points. There are also a variety of smoothing methods, involving some form of splines, or otherwise.
Evolutionary algorithms and various "animal" algorithms are used by some people, but I don't have a very positive view of them. Perhaps they have their use. You seem to have a handle on those, so I won't further elaborate.
Going back to derivative based methods, there are several possibilities. Finite difference gradient are always applicable if the objective function is theoretically differentiable. A method called Likelihood Ratio Method (it was my question over 40 years ago to the then student who developed LRM, of how you could simulate for a derivative, which I wanted for use in optimization, that stimulated him to invent the method) is fairly broadly applicable to directly simulate for derivatives, but is usually not very fast (has high variance), and not used that widely - it basically involves differentiating under the expectation, but instead of differentiating the (unintegrated) objective function, i.e., the integrand, it differentiates the probability measure driving the simulation) A closely related method, known as week derivatives, is essentially a more modern version of Likelihood Ratio Method, but as of yet, no one has created an Automatic Differentiation (AD) like machinery to make it "easy" to use on complicated models - I think that it might be doable by marrying the expertise of an AD expert to modify existing AD tools to incorporate different end node rules needed for weak derivatives (there might be some theoretical applicability issues though).
If the objective function is continuous (smooth) per sample path, not just in expected value, then it is generally valid to interchange differentiation and expectation (simulation evaluates expectation, even when estimating probabilities). This is known as Infinitesimal Perturbation Analysis (IPA). The objective function is differentiated, hen simulated. Automatic Differentiation (a.k.a. Algorithmic Differentiation) can be used to directly simulate for the gradient, especially using the so-called Reverse Mode of Automatic Differentiation, which computes the entire gradient, no matter how many components, at a fixed overhead ratio to just the objective function (presuming three is enough memento). The Reverse Mode of AD applied to Neural Networks is better known as Backpropagation. Sometimes Backpropagation is used to refer just to the Reveres Mode AD, while sometime it is used sued to refer to the whole stochastic gradient descent algorithm using Reverse Mode AD.
IPA using the Reverse Mode of Automatic Differentiation can compute (simulate) an entire gradient "at once". By contrast, LRM calculates the gradient one component at a time, by separately simulating with respect to the probability measure (partially) differentiated with respect to one component at a time. This can be a big disadvantage relative to IPA in models having a large number of optimization variables.
Complex Step Differentiation is a method of doing finite differences using complex numbers. Despite looking very different from Automatic Differentiation, it essentially amounts to the same thing (once the simulation is suitably modified to handle complex inputs), and achieves the same accuracy as Automatic Differentiation,. Unfortunately, its scope of applicability is the same as IPA, not that of finite differences, i.e., it requires each sample path be differentiable, not just in expectation. Further, Complex Step Differentiation really corresponds to the Forward Mode of Automatic Differentiation; and therefore, computes only one component of the gradient at a time. This places it at a disadvantage to RM Automatic Differentiation when the number of optimization variables sis large, but might be advantageous for a small number of variables, due to avoiding the overhead of RM Automatic Differentiation.
If some or all of the optimization variables are discrete, that's a whole different ball game which I won't get into here.
