Here, we propose a strategy for the global optimization of process flowsheets, a fundamental problem in process systems engineering, based on algebraic surrogates that are built from rigorous simulations via Bayesian symbolic regression. The applied method provides a closed-form expression that can be optimized to global optimality using state-of-the-art solvers, where BARON or ANTIGONE were the solvers of choice. When predicting unseen test data, the algebraic models show a similar accuracy level compared to traditional surrogates based on Gaussian processes. However, they can be more easily optimized to global optimality due to their analytical closed-form structure, which allows the user to apply well-established global deterministic solvers. We show the capabilities of our approach in several case studies, ranging from process units to full flowsheets. The performance of our approach is assessed by comparing the CPU time for model building, the prediction accuracy of the identified model, and the CPU time for the subsequent optimization with a proven benchmark.