Optimization methods are essential to improve the operation of energy conversion systems including energy storage equipment and fluctuating renewable energy. Modern systems consist of many components, operating in a wide range of conditions and governed by nonlinear balance equations. Consequently, identifying their optimal operation (e.g. minimizing operational costs) requires solving challenging optimization problems, with the global optimum often hidden behind many local ones. In this work, we propose a hybrid method that advantageously combines Bayesian optimization (BO) and Interior Point OPTimizer (IPOPT). The BO is a global approach exploiting Gaussian process regression to build a surrogate model of the cost function to be optimized, while IPOPT is a local approach using quasi-Newton updates. The proposed BO-IPOPT combination allows leveraging the parameter space exploration of the BO with the quasi-Newton convergence of IPOPT once solution candidates are in the neighborhood of an optimum. Using a challenging constrained test function, we test BO-IPOPT in accuracy and computational efficiency. Finally, we showcase the proposed method in the optimal operation of a renewable steam generation system. The results show that BO-IPOPT combines high accuracy and computational efficiency, achieving up to 50% better objective function values at the same CPU time than other state-of-the-art methods.