Bayesian Optimization (BO) has recently gained popularity as an efficient derivative-free method for the global optimization of expensive noisy black-box objective functions. These characteristics render BO a promising tool to tackle optimization problems involving numerical simulations of complex unsteady flows at moderate-to-high Reynolds numbers. In this work, we assess the efficiency of Bayesian Optimization by considering two canonical flow problems: the drag reduction in the two-dimensional and three-dimensional flow around circular cylinders at, respectively, Re=500 and Re=3900, through tangential-velocity actuation at the cylinder wall. The root-mean-square of the drag coefficient with and without penalty terms is considered as the objective function. Several variants of Bayesian Optimization are assessed and compared against competing optimization algorithms such as Particle Swarm Optimization, CMA-ES, Nelder–Mead and the Explorative Gradient Method. Results show that in this case, the serial and the parallel BO techniques outperform other algorithms.