Simulation and optimization of chemical flowsheets rely on the solution of a large number of non-linear equations. Finding such solutions can be supported by constructing machine-learning based surrogates, relating features and outputs by simple, explicit functions. In order to generate training data for those surrogates computationally efficiently, schemes to adaptively sample the feature space are mandatory. In this article, we present a novel family of utility functions to favor an adaptive, Bayesian exploration of the feature space in order to identify regions that are convergent, fulfill customized inequality constraints and are Pareto-optimal with respect to conflicting objectives. The benefit is illustrated by small toy-examples as well as by industrially relevant chemical flowsheets.