Abstract State engineering of quantum objects is a central requirement for precision sensing and quantum computing implementations. When the quantum dynamics can be described by analytical solutions or simple approximation models, optimal state preparation protocols have been theoretically proposed and experimentally realized. For more complex systems such as interacting quantum gases, simplifying assumptions do not apply anymore and the optimization techniques become computationally impractical. Here, we propose Bayesian optimization based on multi-output Gaussian processes to learn the physical properties of a Bose-Einstein condensate within few simulations only. We evaluate its performance on an optimization study case of diabatically transporting the quantum gas while keeping it in its ground state. Within a few hundred executions, we reach a competitive performance to other protocols. While restricting this benchmark to the well known Thomas-Fermi approximation for straightforward comparisons, we expect a similar performance when employing more complex theoretical models, which would be computationally more challenging, rendering standard optimal control theory protocols impractical. This paves the way for efficient state engineering of complex quantum systems including mixtures of interacting gases or cold molecules.