Decisions on risk mitigation measures, such as reducing the impacts of natural hazards like floods, are often made in the presence of significant uncertainty. Ideally, various risk mitigation strategies are compared to identify the optimal combination of mitigation measures that minimizes an objective function involving risk. Most commonly, a probabilistic, sampling-based model is used to estimate the risk entering this objective function. Due to the significant randomness associated with the outcome of hazard processes and model uncertainties, sampling-based risk estimates exhibit large variance. Consequently, the estimates of the objective function show large variance too, which makes this a challenging stochastic optimization problem. This challenge is aggravated by computational constraints, as risk assessment models are often expensive to run, which limits the number of model evaluations that can be obtained. Modeling and accounting for this sampling uncertainty in the optimization process is a central challenge that we address in our work. We present a Bayesian optimization framework tailored to these settings. The surrogate is a heteroscedastic Gaussian process that explicitly represents input-dependent observation variance through a separate Gaussian process. We apply a Monte Carlo–based Knowledge Gradient as acquisition function, which accounts for uncertainty in the estimated means. Practical challenges and identified solutions arising in the optimization are presented.