Summary In well control (production) optimization, the computational cost of conducting a full-physics flow simulation on a 3D, rich grid-based model poses a significant challenge. This challenge is exacerbated in a robust optimization (RO) setting, where flow simulation must be repeated for numerous geological realizations, rendering RO impractical for many field-scale cases. In this paper, we introduce and discuss a new optimization workflow that addresses this issue by providing computational efficiency, i.e., achieving a near-global optimum of the predefined objective function with minimal forward model (flow-simulation) evaluations. In this workflow, referred to as “Bayesian optimization (BO),” the objective function for samples of decision (control) variables is first computed using a proper design experiment. Then, given the samples, a Gaussian process regression (GPR) is trained to mimic the surface of the objective function as a surrogate model. While balancing the dilemma to select the next control variable between high mean, low uncertainty (exploitation) and low mean, high uncertainty (exploration), a new control variable is selected, and flow simulation is run for this new point. Later, the GPR is updated, given the output of the flow simulation. This process continues sequentially until the termination criteria are satisfied. To validate the workflow and obtain a better insight into the detailed steps, we first optimized a 1D problem. The workflow is then implemented for a 3D synthetic reservoir model to perform RO in a realistic field scenario (8-dimensional and 45-dimensional optimization problems). The workflow is compared with two other commonly used gradient-free algorithms in the literature: particle swarm optimization (PSO) and genetic algorithm (GA). The main contributions are (1) developing a new optimization workflow to address the computational cost of flow simulation in RO, (2) demonstrating the effectiveness of the workflow on a 3D grid-based model, (3) investigating the robustness of the workflow against randomness in initiation samples and discussing the results, and (4) comparing the workflow with other optimization algorithms, showing that it achieves same near-optimal results while requiring only a fraction of the computational time.