Optimized data exploration applied to the simulation of a chemical process

TL;DR

This paper introduces a novel iterative algorithm for exploring complex parameter spaces in chemical process simulations, improving feasibility classification and guiding exploration towards regions of interest using machine learning techniques.

Contribution

The paper presents a new exploration algorithm combining machine learning methods to efficiently identify feasible regions in unknown parameter spaces, with improvements over existing Kriging-based approaches.

Findings

Outperforms Kriging-based exploration in binary feasibility classification

Successfully applied to industrial chemical process simulation

Achieves good data space approximation with few data points

Abstract

In complex simulation environments, certain parameter space regions may result in non-convergent or unphysical outcomes. All parameters can therefore be labeled with a binary class describing whether or not they lead to valid results. In general, it can be very difficult to determine feasible parameter regions, especially without previous knowledge. We propose a novel algorithm to explore such an unknown parameter space and improve its feasibility classification in an iterative way. Moreover, we include an additional optimization target in the algorithm to guide the exploration towards regions of interest and to improve the classification therein. In our method we make use of well-established concepts from the field of machine learning like kernel support vector machines and kernel ridge regression. From a comparison with a Kriging-based exploration approach based on recently published results we can show the advantages of our algorithm in a binary feasibility classification scenario with a discrete feasibility constraint violation. In this context, we also propose an improvement of the Kriging-based exploration approach. We apply our novel method to a fully realistic, industrially relevant chemical process simulation to demonstrate its practical usability and find a comparably good approximation of the data space topology from relatively few data points.

Tables

Table 1. Table 1: Toy data exploration setups (1) to (5) with their respective exploration hyperparameters. We list the total number of samples N 𝑁 N and the parameter G 𝐺 G controlling the number of initial samples G 2 superscript 𝐺 2 G^{2} , which are placed on a regular grid in 𝝌 G = 𝝌 toy subscript 𝝌 G subscript 𝝌 toy \bm{\chi}_{\mathrm{G}}=\bm{\chi}_{\mathrm{toy}} , Eq.   20 . The weights s 𝑠 s , o 𝑜 o and r 𝑟 r constituting the weight vector 𝐰 𝐰 \mathbf{w} , Eq.   11 , determine the influence of the three components of the utility function, Eq.   8 . Specifically, s 𝑠 s controls the outcome prediction uncertainty impact, o 𝑜 o the optimization target prediction impact and r 𝑟 r the distance to the nearest neighbor impact, respectively, Eq.   12 . Setups (1) to (4) use our algorithm, whereas setup (5) represents a pure grid approach in the framework of our algorithm, which is finished after the initial sampling step since G 2 = N superscript 𝐺 2 𝑁 G^{2}=N .

Setup	$N$	$G$	$s$	$o$	$r$
(1)	$16$	$3$	1	0	1
(2)	$25$	$3$	1	0	1
(3)	$64$	$3$	1	0	1
(4)	$25$	$3$	1	3	1
(5)	$25$	$5$	0	0	0

Table 2. Table 2: Toy data exploration characteristics for each of the setups (1) to (5) from table   1 . We list the success rate R 𝑅 R , Eq.   31 , the score σ 𝜎 \sigma , Eq.   23 , the ratio of false positives r fp subscript 𝑟 fp r_{\mathrm{fp}} , Eq.   25a , the ratio of false negatives r fn subscript 𝑟 fn r_{\mathrm{fn}} , Eq.   25b , the validity ratio α 𝛼 \alpha , Eq.   27 , the best optimization score t best subscript 𝑡 best t_{\mathrm{best}} , Eq.   30a , and the worst optimization score t worst subscript 𝑡 worst t_{\mathrm{worst}} , Eq.   30b . The arrows indicate whether high values ( ↑ ↑ \uparrow ) or low values ( ↓ ↓ \downarrow ) are considered more favorable for each characteristic. Depending on the context, either a low value of t worst subscript 𝑡 worst t_{\mathrm{worst}} (e. g., for o = 0 𝑜 0 o=0 ) or a high value (e. g., for o > 0 𝑜 0 o>0 ) can be desired.

Setup	$R$ $\uparrow$	$\sigma$ $\uparrow$	$r_{\mathrm{fp}}$ $\downarrow$	$r_{\mathrm{fn}}$ $\downarrow$	$\alpha$ $\uparrow$	$t_{\mathrm{best}}$ $\uparrow$	$t_{\mathrm{worst}}$
(1)	$0.807\left(8\right)$	$1.000$	$0.167\left(7\right)$	$0.026\left(3\right)$	$0.455$	$1.000$	$-1.000$
(2)	$0.857\left(7\right)$	$1.000$	$0.102\left(4\right)$	$0.041\left(4\right)$	$0.389$	$1.000$	$-1.568$
(3)	$0.952\left(4\right)$	$1.000$	$0.039\left(4\right)$	$0.009\left(2\right)$	$0.524$	$1.561$	$-1.568$
(4)	$0.761\left(8\right)$	$1.000$	$0.220\left(7\right)$	$0.019\left(3\right)$	$0.786$	$1.716$	$-1.000$
(5)	$0.797\left(8\right)$	$1.000$	$0.178\left(7\right)$	$0.025\left(3\right)$	$0.389$	$1.000$	$-1.000$

Table 3. Table 3: Benchmark exploration characteristics for the three setup collections ( 6 ) n subscript 6 𝑛 (6)_{n} , [ 7 ] n subscript delimited-[] 7 𝑛 [7]_{n} , and { 8 } n subscript 8 𝑛 \{8\}_{n} , where n 𝑛 n represents the total number of sampled points N 𝑁 N . The setups ( 6 ) 25 subscript 6 25 (6)_{25} , ( 6 ) 50 subscript 6 50 (6)_{50} , ( 6 ) 100 subscript 6 100 (6)_{100} and ( 6 ) 150 subscript 6 150 (6)_{150} use our algorithm, whereas [ 7 ] 25 subscript delimited-[] 7 25 [7]_{25} , [ 7 ] 50 subscript delimited-[] 7 50 [7]_{50} , [ 7 ] 100 subscript delimited-[] 7 100 [7]_{100} and [ 7 ] 150 subscript delimited-[] 7 150 [7]_{150} use the Kriging-based algorithm. The setups { 8 } 25 subscript 8 25 \{8\}_{25} , { 8 } 50 subscript 8 50 \{8\}_{50} , { 8 } 100 subscript 8 100 \{8\}_{100} and { 8 } 150 subscript 8 150 \{8\}_{150} use the modified Kriging-based algorithm. For each setup, we show the success rate R 𝑅 R , Eq.   31 , the score σ 𝜎 \sigma , Eq.   23 , the ratio of false positives r fp subscript 𝑟 fp r_{\mathrm{fp}} , Eq.   25a , the ratio of false negatives r fn subscript 𝑟 fn r_{\mathrm{fn}} , Eq.   25b , the validity ratio α 𝛼 \alpha , Eq.   27 , the best optimization score t best subscript 𝑡 best t_{\mathrm{best}} , Eq.   30a , and the worst optimization score t worst subscript 𝑡 worst t_{\mathrm{worst}} , Eq.   30b . The arrows indicate whether high values ( ↑ ↑ \uparrow ) or low values ( ↓ ↓ \downarrow ) are considered more favorable for each characteristic. Depending on the context, either a low or a high value of t worst subscript 𝑡 worst t_{\mathrm{worst}} might be preferable. According to these characteristics, our algorithm is comparable to the Kriging-based approach for few samples, but superior for many. The modified Kriging-based approach generally performs better than the original version. A visualization of exploration characteristics for different sample sizes can be found in Figure   3 .

Setup	$R$ $\uparrow$	$\sigma$ $\uparrow$	$r_{\mathrm{fp}}$ $\downarrow$	$r_{\mathrm{fn}}$ $\downarrow$	$\alpha$ $\uparrow$	$t_{\mathrm{best}}$ $\uparrow$	$t_{\mathrm{worst}}$
${\left(6\right)}_{25}$	$0.816\left(8\right)$	$1.000$	$0.029\left(3\right)$	$0.155\left(7\right)$	$0.136$	$0.667$	$-1.333$
${\left(6\right)}_{50}$	$0.924\left(5\right)$	$1.000$	$0.026\left(3\right)$	$0.050\left(4\right)$	$0.667$	$1.435$	$-1.606$
${\left(6\right)}_{100}$	$0.975\left(3\right)$	$1.000$	$0.011\left(2\right)$	$0.014\left(2\right)$	$0.887$	$1.580$	$-1.723$
${\left(6\right)}_{150}$	$0.988\left(2\right)$	$0.973$	$0.005\left(1\right)$	$0.006\left(2\right)$	$0.875$	$1.696$	$-1.724$
${\left[7\right]}_{25}$	$0.841\left(7\right)$	$1.000$	$0.125\left(6\right)$	$0.034\left(4\right)$	$0.471$	$1.229$	$-1.225$
${\left[7\right]}_{50}$	$0.894\left(6\right)$	$1.000$	$0.074\left(5\right)$	$0.032\left(3\right)$	$0.471$	$1.558$	$-1.568$
${\left[7\right]}_{100}$	$0.950\left(4\right)$	$1.000$	$0.018\left(3\right)$	$0.032\left(3\right)$	$0.351$	$1.558$	$-1.568$
${\left[7\right]}_{150}$	$0.960\left(4\right)$	$1.000$	$0.014\left(2\right)$	$0.026\left(3\right)$	$0.339$	$1.558$	$-1.568$
${\left\{8\right\}}_{25}$	$0.880\left(6\right)$	$1.000$	$0.076\left(5\right)$	$0.044\left(4\right)$	$0.471$	$1.220$	$-1.226$
${\left\{8\right\}}_{50}$	$0.922\left(5\right)$	$1.000$	$0.053\left(4\right)$	$0.025\left(3\right)$	$0.613$	$1.481$	$-1.524$
${\left\{8\right\}}_{100}$	$0.957\left(4\right)$	$1.000$	$0.021\left(3\right)$	$0.022\left(3\right)$	$0.587$	$1.571$	$-1.530$
${\left\{8\right\}}_{150}$	$0.975\left(3\right)$	$0.993$	$0.013\left(2\right)$	$0.012\left(2\right)$	$0.596$	$1.631$	$-1.645$

Table 4. Table 4: Data exploration setups for the chemical simulation. We list the setups (1) to (4), which use our algorithm, the setups (5) and (6), which represent a pure grid approach in the framework of our algorithm, and the setups ⟨ ⟨ \langle 7 ⟩ ⟩ \rangle to ⟨ ⟨ \langle 9 ⟩ ⟩ \rangle , which are based on a LHS approach, with their respective exploration hyperparameters. For each setup, we show the total number of samples N 𝑁 N . For setups (1) to (6) we also show and the parameter G 𝐺 G controlling the number of initial samples G 4 superscript 𝐺 4 G^{4} , which are placed on a regular grid in 𝝌 G subscript 𝝌 G \bm{\chi}_{\mathrm{G}} . The latter either corresponds to 𝝌 sim subscript 𝝌 sim \bm{\chi}_{\mathrm{sim}} , Eq.   39 , or 𝝌 win subscript 𝝌 win \bm{\chi}_{\mathrm{win}} , Eq.   41 . The weights s 𝑠 s , o 𝑜 o and r 𝑟 r determine the explorative behavior by changing the influence of the three components of the utility function, Eq.   8 . Specifically, s 𝑠 s controls the outcome prediction uncertainty impact, o 𝑜 o the optimization target prediction impact and r 𝑟 r the distance to the nearest neighbor impact, respectively, Eq.   12 . All of these exploration hyperparameters have no meaning for the LHS setups.

Setup	$N$	$G$	${𝝌}_G$	$s$	$o$	$r$
(1)	$40$	$2$	${𝝌}_{\mathrm{win}}$	1	0	1
(2)	$40$	$2$	${𝝌}_{\mathrm{win}}$	1	1	1
(3)	$81$	$2$	${𝝌}_{\mathrm{win}}$	1	0	1
(4)	$81$	$2$	${𝝌}_{\mathrm{win}}$	1	1	1
(5)	$81$	$3$	${𝝌}_{\mathrm{sim}}$	0	0	0
(6)	$256$	$4$	${𝝌}_{\mathrm{sim}}$	0	0	0
$⟨$7$⟩$	$40$	LHS in ${𝝌}_{\mathrm{sim}}$
$⟨$8$⟩$	$81$	LHS in ${𝝌}_{\mathrm{sim}}$
$⟨$9$⟩$	$256$	LHS in ${𝝌}_{\mathrm{sim}}$

Table 5. Table 5: Data exploration characteristics for the chemical simulation for each of the setups (1) to (4), which use our algorithm, the setups (5) and (6), which represent a pure grid approach in the framework of our algorithm, and the setups ⟨ ⟨ \langle 7 ⟩ ⟩ \rangle to ⟨ ⟨ \langle 9 ⟩ ⟩ \rangle , which are based on a LHS approach. The respective hyperparameters are shown in  table   4 . We list the success rate R 𝑅 R , Eq.   31 , the score σ 𝜎 \sigma , Eq.   23 , the ratio of false positives r fp subscript 𝑟 fp r_{\mathrm{fp}} , Eq.   43a , the ratio of false negatives r fn subscript 𝑟 fn r_{\mathrm{fn}} , Eq.   43b , the validity ratio α 𝛼 \alpha , Eq.   27 , the best optimization score t best subscript 𝑡 best t_{\mathrm{best}} , Eq.   30a , and the worst optimization score t worst subscript 𝑡 worst t_{\mathrm{worst}} , Eq.   30b . The arrows indicate whether high values ( ↑ ↑ \uparrow ) or low values ( ↓ ↓ \downarrow ) are considered more favorable for each characteristic. Depending on the context, either a low or a high value of t worst subscript 𝑡 worst t_{\mathrm{worst}} might be preferable. These characteristics show that our algorithm provides us with a very good estimation of the feasibility region from relatively few data points in comparison with grid-based or LHS approaches.

Setup	$R$ $\uparrow$	$\sigma$ $\uparrow$	$r_{\mathrm{fp}}$ $\downarrow$	$r_{\mathrm{fn}}$ $\downarrow$	$\alpha$ $\uparrow$	$t_{\mathrm{best}}$ $\uparrow$	$t_{\mathrm{worst}}$
(1)	$0.913\left(6\right)$	$1.000$	$0.069\left(5\right)$	$0.018\left(3\right)$	$0.739$	$0.841$	$0.588$
(2)	$0.863\left(7\right)$	$1.000$	$0.134\left(6\right)$	$0.003\left(1\right)$	$1.105$	$0.894$	$0.584$
(3)	$0.921\left(5\right)$	$1.000$	$0.047\left(4\right)$	$0.033\left(3\right)$	$0.841$	$0.833$	$0.573$
(4)	$0.724\left(9\right)$	$1.000$	$0.276\left(7\right)$	$0.000\left(1\right)$	$1.314$	$0.930$	$0.582$
(5)	$0.479\left(10\right)$	$1.000$	$0.000\left(1\right)$	$0.521\left(7\right)$	$0.066$	$0.725$	$0.725$
(6)	$0.485\left(10\right)$	$1.000$	$0.000\left(1\right)$	$0.515\left(7\right)$	$0.138$	$0.817$	$0.633$
$⟨$7$⟩$	$0.480\left(10\right)$	$0.875$	$0.000\left(1\right)$	$0.521\left(7\right)$	$0.481$	$0.847$	$0.567$
$⟨$8$⟩$	$0.521\left(10\right)$	$0.988$	$0.479\left(7\right)$	$0.000\left(1\right)$	$0.350$	$0.856$	$0.585$
$⟨$9$⟩$	$0.667\left(9\right)$	$0.984$	$0.333\left(8\right)$	$0.000\left(1\right)$	$0.326$	$0.902$	$0.556$

Equations

$$\mathcal{S}:\mathbf{\chi }\to S$$

$$S\equiv \eta \otimes \tau .$$

$$\eta \equiv \{\mathrm{valid},\mathrm{invalid}\}$$

$$d(\mathbf{x})\equiv (\mathbf{x},y,t)$$

$$D\equiv D({\mathbf{x}}_1,\dots ,{\mathbf{x}}_n)\equiv \{d({\mathbf{x}}_1),\dots ,d({\mathbf{x}}_n)\}.$$

$$D_x\equiv \{{\mathbf{x}}_1,\dots ,{\mathbf{x}}_n\}$$

$$\mathcal{E}:\mathcal{D}\otimes \mathbf{\chi }\to S\otimes {\eta }_p$$

$$U(\mathcal{E},D,\mathbf{x},\mathbf{w})$$

$${\text{x}}_{\mathrm{new}}=\underset{\text{x}\in \mathbf{\chi }}{\mathrm{argmax}}U(\mathcal{E},D,\mathbf{x},\mathbf{w})$$

$$\mathbf{u}\equiv \mathbf{u}(\mathcal{E},D,\mathbf{x})\equiv \left(\begin{array}{c}{U}_s(\mathcal{E},D,\mathbf{x})\\ U_o(\mathcal{E},D,\mathbf{x})\\ U_r(D,\mathbf{x})\end{array}\right)$$

$$\mathbf{w}\equiv \mathbf{w}(s,o,r)\equiv \left(\begin{array}{c}s\\ o\\ r\end{array}\right),$$

$$\mathbf{u}\sim \left(\begin{array}{c}\text{Outcome prediction uncertainty}\\ \text{Optimization target prediction}\\ \text{Distance to nearest neighbor}\end{array}\right).$$

$$U_s(\mathcal{E},D,\mathbf{x})$$

$$S({\hat{p}}_{\hat{y}})$$

$$U_o(\mathcal{E},D,\mathbf{x})$$

$${\hat{t}}_{\mathrm{r}}(\mathcal{E},D,\mathbf{x})$$

$$t_{\max }(D)$$

$$t_{\min }(D)$$

$$U_r(D,\mathbf{x})$$

$$\mathbf{x}\equiv \left(\begin{array}{c}{x}_1\\ x_2\end{array}\right)$$

$${\mathbf{\chi }}_{\mathrm{toy}}\equiv [-2,2]\otimes [-2,2].$$

$$y({\mathcal{S}}_{\mathrm{toy}},\mathbf{x})\equiv \{\begin{array}{ll}\mathrm{valid}& \text{if}||\mathbf{x}|{|}_2\le \sqrt{3}\wedge x_1{x}_2\ge 0\\ \mathrm{invalid}& \text{otherwise}\end{array}$$

$$t({\mathcal{S}}_{\mathrm{toy}},\mathbf{x})\equiv x_1,$$

$$R\equiv R(D_{\mathrm{expl}})$$

$$\sigma \equiv \sigma (D_{\mathrm{expl}})$$

$$g\equiv \frac{\sigma }{R}$$

$$r_{\mathrm{fp}}\equiv r_{\mathrm{fp}}(D_{\mathrm{expl}})$$

$$r_{\mathrm{fn}}\equiv r_{\mathrm{fn}}(D_{\mathrm{expl}})$$

$$r_{\mathrm{fp}}+r_{\mathrm{fn}}=1-R$$

$$\alpha \equiv \alpha (D_{\mathrm{expl}})$$

$$\displaystyle\alpha_{\infty}=\frac{3\pi}{32-3\pi}\approx$0.417$$

$

$$\displaystyle r_{\mathrm{fp}\infty}=\frac{3\pi}{64}\approx$0.147$$

$

Full text

Optimized data exploration applied to the simulation of a chemical process

Raoul Heesea,b,111Corresponding author: raoul.heese@itwm.fraunhofer.de, Michał Walczaka,b, Tobias Seidelb,

Norbert Asprionc, Michael Bortza,b

*a**Fraunhofer Center for Machine Learning

bFraunhofer ITWM Optimization Department, Fraunhofer-Platz 1, 67663 Kaiserslautern, Germany

cChemical and Process Engineering BASF SE, Carl-Bosch-Str. 38, 67056 Ludwigshafen, Germany*

Abstract

In complex simulation environments, certain parameter space regions may result in non-convergent or unphysical outcomes. All parameters can therefore be labeled with a binary class describing whether or not they lead to valid results. In general, it can be very difficult to determine feasible parameter regions, especially without previous knowledge. We propose a novel algorithm to explore such an unknown parameter space and improve its feasibility classification in an iterative way. Moreover, we include an additional optimization target in the algorithm to guide the exploration towards regions of interest and to improve the classification therein. In our method we make use of well-established concepts from the field of machine learning like kernel support vector machines and kernel ridge regression. From a comparison with a Kriging-based exploration approach based on recently published results we can show the advantages of our algorithm in a binary feasibility classification scenario with a discrete feasibility constraint violation. In this context, we also propose an improvement of the Kriging-based exploration approach. We apply our novel method to a fully realistic, industrially relevant chemical process simulation to demonstrate its practical usability and find a comparably good approximation of the data space topology from relatively few data points.

1 Introduction

Chemical process design is in its nature a multicriteria optimization (MCO) task [16, 4]. To find a Pareto-optimal design [10], an engineer needs to explore multiple trade-offs between competing (design) parameters and objectives. Recently, powerful decision support tools have been developed for usage in industrial contexts to aid an engineer in this design process. These tools enable him to navigate through Pareto-efficient solutions from simulation runs performed for different parameter variations [4, 7].

In a typical industrial application, a flowsheet simulation of a chemical process involves solving several hundreds to several thousands of nonlinear equations for each specified parameter combination. Only certain combinations of the parameters lead to numerically converging, physically reasonable results. Therefore, in order to estimate the Pareto frontier it is necessary to determine valid parameter combinations, i. e., to find the operation window or feasibility region in the parameter space.

A binary feasibility classification task in this context can be challenging due to an often vast parameter space and the computational effort involved in solving the system of nonlinear equations. There exist different approaches from various branches of engineering how to solve this problem. Many studies propose to reduce the computational cost by approximating the feasible region. For example, in Ref. [12] parameters are classified according to their feasibility using a support vector machine (SVM) [9] combined with k-means clustering to reduce the model training cost.

An alternative to binary classification of feasible points is to train a surrogate model of a continuous feasibility function, which reflects by how much the constraints which define the feasibility range are violated. The surrogate model is usually trained adaptively by sampling regions with high model uncertainty and the vicinity of the predicted boundary. Suggested ways to construct a surrogate model include high dimensional model representation techniques [1], Kriging-based modelling [5, 6], and most recently, radial basis function (RBF)-based modelling [30]. The RBF-based modeling outperforms Kriging-based modeling in terms of achieving better accuracy with fewer sampling points, however, both approaches exhibit limitations in five- and higher-dimensional problems [30]. Moreover, in certain applications it is impossible to access the information about maximum constraints violation.

In this manuscript, we present a novel method to sequentially explore feasibility classifications. Inspired by Ref. [6] our proposed algorithm does not only provide us with an adaptive estimator for feasibility, but can also take an optimization target into account. This allows us to focus the exploration on parameter regions of interest. We therefore consider our approach as a method for optimized data exploration, which exceeds bare feasibility classification. Briefly put, we propose an adaptive screening of an operation window combining feasibility classification and optimization.

For demonstration purposes we apply our method to the simulation of a realistic chemical process. In particular, we consider the simulation of two serialized distillation columns to separate an azeotropic mixture of chloroform and acetone [3]. We show that we can improve the estimation of feasible regions in comparison with a uniform grid approach and a latin hypercube sampling (LHS) [15] approach. We use Chemasim, a BASF in-house flowsheet simulator, to run the simulations. By design, Chemasim, like other flowsheet simulators, does not provide a complete quantification of equation violations for unfeasible parameters, hence a surrogate model of a continuous feasibility function can not be applied. Our algorithm, however, is based entirely on a binary feasibility classification, which we can directly extract from Chemasim.

In the following, we will first present a formal description of our algorithm. Subsequently, we will explain its functionality with the help of a toy example. Using a benchmark, we will show the strengths of our method in comparison with a Kriging-based exploration approach for a binary feasibility classification scenario. We also propose an improvement of this Kriging-based approach for a discrete feasibility constraint violation. Finally, we will outline the chemical process simulation to which we have applied our algorithm and present the results. We will conclude with a brief summary and outlook.

2 Optimized data exploration

Our method of optimized data exploration can be considered as a sequential design of experiments, where each new experiment (i. e., each new simulation evaluation) is chosen based on a utility function. Summarized, the method consists of four major steps:

1. (i)

Start-up: Ensure that an initial data set is available. This preliminary step is required to ensure at least basic knowledge about the data topology. Therefore, one may use previously obtained data or evaluate a set of parameters, e. g., on a regular grid. 2. (ii)

Choice: Choose the parameter for which the utility function is maximized. The utility function makes use of the previously obtained data to interpolate or extrapolate missing information. 3. (iii)

Evaluation: Evaluate the simulation for the chosen parameter. 4. (iv)

Repetition: Repeat from (ii) or finish the exploration based on a termination criterion.

The choice of a suitable utility function and a sufficient initial knowledge about the data are of course crucial for the success of our approach.

From a practical perspective, data exploration is often driven by an optimization, hence we consider in the following that we will not only evaluate a simulation outcome for its validity but will also extract some kind of quantitative optimization target from the result. In order to emphasize on the main aspects of our data exploration method, we will only consider a single optimization target. A MCO problem can, however, always be reduced to this one-dimensional case by an appropriate scalarization of its targets.

In this section, we will first describe the general framework of our method in which we introduce the main formal ingredients. We will subsequently explain the exploration algorithm itself in more detail. A toy example will help us to demonstrate our method.

2.1 Framework

We use simulation evaluations and estimators to generate data. Therefore, we will first formally define those concepts. Furthermore, we will discuss our choice of a utility function. These definitions will serve as a framework for the following studies.

2.1.1 Data generation

We consider a simulation $\mathcal{S}$ which is described by a mapping

[TABLE]

of parameters $\mathbf{x}$ from the compact $p$-dimensional parameter space $\bm{\chi}\subset\mathbb{R}^{p}$ onto the solution space

[TABLE]

This solution space consists of two parts. First, a classification space

[TABLE]

containing two classes which describe whether the simulation outcome $y\equiv y(\mathcal{S},\mathbf{x})\in\eta$ was valid (i. e., numerically convergent and physically reasonable in the sense that the numerically obtained simulation result fulfills a set of predefined conditions) or invalid. And second, the optimization target space $\tau\subset\mathbb{R}$ which contains all possible results for the optimization target $t\equiv t(\mathcal{S},\mathbf{x})\in\tau$. The optimization target is meaningful only for valid simulation outcomes. Without loss of generality we assume in the following that a maximization of $t$ is considered optimal. The symbol $\otimes$ in Eq. 2 denotes a tensor product.

Each evaluation of the simulation therefore leads to a data point

[TABLE]

given by the evaluated parameter $\mathbf{x}$, the corresponding outcome $y$ and the optimization target $t$ from the mapping $\mathcal{S}$, Eq. 1. Although the optimization target is meaningless for invalid outcomes, we still include it in $d(\mathbf{x})$ to achieve a unified notation for valid and invalid data points. Note that we assume that the result of a simulation is purely deterministic and completely defined by the choice of parameters $\mathbf{x}$. The collection of results from $n$ evaluations consequently allows us to define a data set

[TABLE]

Data exploration is achieved by adding new elements to such a set. For reasons of convenience we use

[TABLE]

to denote the corresponding collection of evaluated parameters.

2.1.2 Prediction

A sufficiently large data set allows us to perform model predictions for solutions of yet unevaluated parameters with the help of an estimator

[TABLE]

which maps from the space of possible data sets $\mathcal{D}$ and the parameter space $\bm{\chi}$ onto the solution space $S$, Eq. 2, and the predicted outcome probability space $\eta_{p}\equiv[0,1]$. The solution space contains the predicted outcome $\hat{y}\equiv\hat{y}(\mathcal{E},D,\mathbf{x})\in\eta$ and the predicted optimization target $\hat{t}\equiv\hat{t}(\mathcal{E},D,\mathbf{x})\in\tau$, whereas the predicted outcome probability space contains the probability of the predicted outcomes $\hat{p}_{\hat{y}}\equiv\hat{p}_{\hat{y}}(\mathcal{E},D,\mathbf{x})\in[0,1]$. Briefly put, $\mathcal{E}$ allows us to predict simulation outcomes with their associated probability and optimization targets for arbitrary parameters based on previous simulation results.

So far, $\mathcal{E}$ is considered completely general and our algorithm is not limited to a specific choice. However, it has turned out in practice that $\mathcal{E}$ is best represented by two independent estimators: First, a kernel SVM used for a classification of the outcome $\hat{y}$ and second, a kernel ridge regression (RR) [17] to predict the optimization target $\hat{t}$. The probability of the predicted outcomes $\hat{p}_{\hat{y}}$ is then obtained from Platt scaling [20]. The kernel RR can be understood as a surrogate model for the optimization target. For both of these estimations we use the well-known kernel method [11], which is discussed in A in more detail.

2.1.3 Utility function

At the heart of our data exploration method lies the utility function

[TABLE]

which describes the estimated exploration benefit of evaluating the parameter $\mathbf{x}$ based on a previously obtained data set $D$. Therefore, in each sequential exploration step a new evaluation of the simulation is chosen for the parameter

[TABLE]

with the best utility score. The two ingredients of the utility function are the utility vector

[TABLE]

and the weight vector

[TABLE]

which both consist of three components. Each component of the utility vector can be assigned a straightforward interpretation:

[TABLE]

The components (or weights) $s\geq 0$, $o\geq 0$ and $r\geq 0$ of the weight vector determine the influence of $U_{s}(\mathcal{E},D,\mathbf{x})$, $U_{o}(\mathcal{E},D,\mathbf{x})$ and $U_{r}(D,\mathbf{x})$, respectively, on the utility function. In other words, $\mathbf{w}$ determines the explorative behavior. The expression $||\mathbf{w}||_{1}=s+o+r$ in Eq. 8 represents the 1-norm of $\mathbf{w}$ and since by definition $\mathbf{u}\in[0,1]^{3}$, as we will see below, one has $U(\mathcal{E},D,\mathbf{x},\mathbf{w})\in[0,1]$. In the following, we will formally define the components of $\mathbf{u}$, Eq. 10.

The first component

[TABLE]

represents the estimated outcome prediction uncertainty based on the Shannon information entropy [27]

[TABLE]

in bits. Here we have recalled the estimated probability of predicting a valid or invalid outcome $\hat{p}_{\hat{y}}(\mathcal{E},D,\mathbf{x})$ from the estimator mapping, Eq. 7.

The second component

[TABLE]

represents the estimated optimization score. It is based on

[TABLE]

which makes use of the extremal optimization targets

[TABLE]

respectively, to rescale the predicted optimization target $\hat{t}(\mathcal{E},D,\mathbf{x})$, from Eq. 7, in such a way that $\hat{t}_{\mathrm{r}}(\mathcal{E},D,\mathbf{x})\in[0,1]$.

As a last component, the utility vector contains the classification feature space distance [25]

[TABLE]

Here, the expression $||\cdot||_{2}$ stands for the 2-norm distance and $\gamma_{\mathrm{C}}$ represents a hyperparameter of the kernel SVM estimator. A derivation of Eq. 18 can be found in B.

By definition, $U_{r}(D,\mathbf{x})$ becomes smaller the more similar the parameter $\mathbf{x}$ is to its nearest neighbor in $D_{x}$, Eq. 6. Therefore, this expression can also be understood as an artificial repulsion of data points which ensures that newly suggested parameters $\textbf{x}_{\mathrm{new}}$ explore unknown regions of the parameter space $\bm{\chi}$. The amount of utility reduction with decreasing distance is controlled by $\gamma_{\mathrm{C}}$. In A this hyperparameter is discussed in more detail. Note that it is determined during the training of $\mathcal{E}$, Eq. 7, as explained further below.

Summarized, we have introduced the utility function $U(\mathcal{E},D,\mathbf{x},\mathbf{w})$, Eq. 8, based on the utility vector $\mathbf{u}$, Eq. 10, and the weight vector $\mathbf{w}$, Eq. 11. The utility vector consists of three components with a distinct meaning, Eq. 12. These components are weighted by the weights $s$, $o$ and $r$ contained in the weight vector. Consequently, we can control the explorative behavior by tuning the weights.

2.2 Algorithm

Our data exploration method is outlined in Algorithm 1. As described in the introduction of this section, it can be partitioned in four major steps. These steps correspond to the following lines of the algorithm:

1. (i)

Start-up: Line 2 to 11 2. (ii)

Choice: Line 13 3. (iii)

Evaluation: Line 14 4. (iv)

Repetition: Line 12 and 16

As a termination criterion we use a desired number of newly evaluated points $N$. We assume here that the simulation $\mathcal{S}$, Eq. 1, is defined by the application and is not modified during the exploration process. By contrast, the estimator $\mathcal{E}$, Eq. 7, is retrained in each iteration step. For the sake of simplicity we omit $\mathcal{S}$ and $\mathcal{E}$ in the notation.

2.2.1 Functions

The entry point to the algorithm is the Exploration function. Its arguments represent an initial data set from previous evaluations $D_{\mathrm{init}}$ (which might also be an empty set), the parameter space to explore $\bm{\chi}$, the parameter $G\geq 0$ controlling the number of initial evaluations $G^{p}$ in the hyperrectangular compact set $\bm{\chi}_{\mathrm{G}}\subseteq\bm{\chi}$, the desired total number of points to evaluate $N\geq G^{p}$ and the weight vector $\mathbf{w}$, Eq. 11, which controls the behavior of the utility function, Eq. 8. Furthermore, the Exploration function contains three implicit functions:

• •

Simulation: Performs the mapping $\mathcal{S}$ and returns the corresponding data point $d(\mathbf{x})$.

• •

Suggestion: Performs two steps in order to obtain the next parameter to evaluate. First, the estimator $\mathcal{E}$ for the current data set $D_{\mathrm{expl}}$ is trained. Second, Eq. 9 is evaluated and the parameter $\textbf{x}_{\mathrm{new}}$ with the best utility is returned.

• •

Grid: Returns a set of $G^{p}\leq N$ parameters which constitute a regular grid in the parameter subspace $\bm{\chi}_{\mathrm{G}}$. If previous knowledge about the data topology is available, it is generally reasonable to choose a starting grid in such a way that both valid and invalid solutions are sampled in a region of an expectably good optimization target.

After $N$ points have been evaluated, the explored data set $D_{\mathrm{expl}}\supset D$ is returned.

The simulation $\mathcal{S}$, Eq. 1, is considered completely general up to this point. Therefore, our algorithm is universal and can be applied in many different scenarios. We will specify simulation mappings further below in explicit examples for data exploration.

2.2.2 Estimator training

During the training of $\mathcal{E}$ we tune the hyperparameters of the SVM and the RR independently in each iteration step by cross validation [11] as soon as the size of the data set $D_{\mathrm{expl}}$ allows it. For training and prediction we use standardized parameters by removing the mean and scaling to unit variance. All available data at each iteration step is considered as training data.

The optimization target $t$ corresponding to an invalid outcome $y$ is meaningless, but it might nevertheless be of practical use to be able to assign some numerical value to it, e. g., to train an estimator. Therefore, we assume in the following that invalid optimization targets taken from a data set $D$ correspond to the worst valid optimization target of this data set $t_{\mathrm{min}}(D)$, Eq. 17b, for all practical purposes.

Moreover, we assume that the initial data set $D_{\mathrm{init}}$ or the data set obtained from the initial grid sampling is sufficiently large so that a suitable estimator of our choice is well-defined. Otherwise, either the start-up step has to be changed accordingly or a different estimator has to be chosen. We will not further discuss such pathological cases.

3 Demonstration

To illustrate our data exploration method from the previous section, we will first present a two-dimensional toy example (i. e., $p=2$). Specifically, we consider parameters

[TABLE]

in the parameter space

[TABLE]

Thus, the parameter space is a simple square of edge length 4.

3.1 Toy simulation

The chosen toy simulation $\mathcal{S}_{\mathrm{toy}}$, Eq. 1, is given by

[TABLE]

respectively, where $||\mathbf{x}||_{2}$ denotes the 2-norm of $\mathbf{x}$. Thus, simulation outcomes $y$ are considered valid within and on a circle of radius $\sqrt{3}$ in the first and third quadrant and invalid otherwise as shown in Figure 1. This particular example is interesting to study data exploration behavior because it has both sharp and round edges and two distinct feasibility areas only connected by a single point. The optimization target $t$ directly corresponds to $x_{1}$.

We assume no previous knowledge about the data so that $D_{\mathrm{init}}=\{\}$. The initial grid is chosen to expand uniformly across the whole parameter space $\bm{\chi}_{\mathrm{G}}=\bm{\chi}_{\mathrm{toy}}$, Eq. 20. To solve Eq. 9 numerically in each iteration step we use a differential evolution approach [28]. Although our proposed algorithm, Algorithm 1, is completely deterministic in its general form, the usage of cross validation for the estimator training, the statistics involved in the Platt scaling and a differential evolution approach to solve the optimization problem introduce a certain degree of randomness. The presented results are therefore chosen as representative examples.

3.2 Results

The results are shown in Figure 2. Each row, (1) to (5), corresponds to a different set of exploration hyperparameters $N$, $G$, $\bm{\chi}_{\mathrm{G}}$, $s$, $o$ and $r$ as summarized in table 1. We will hereafter refer to a set of such hyperparameters as a setup. Each setup yields a different explored data set $D_{\mathrm{expl}}$. Column (a) shows the evaluated outcomes $y$, Eq. 21a, and labels them as valid ( ) or invalid ( ). The contours ( ) separate the valid from the invalid predicted outcomes $\hat{y}$, Eq. 7, as given by the final estimator. The corresponding regions are shaded accordingly ( and ). Column (b) shows the optimization target $t$, Eq. 21b, for valid outcomes as color-interpolated markers from the best possible result $t=\sqrt{3}$ ( ) to the worst possible result $t=-\sqrt{3}$ ( ). Also shown are the invalid outcomes ( ). For the first and second row we also mark the iteration number $n+1$ from Algorithm 1, i. e., the chronological order in which the points have been evaluated after the initial grid had been set up.

Apparently, setup (1) already provides us with a rough estimate of the parameter topology with 16 evaluations, which is then refined by setups (2) and (3) with 9 and 39 additional evaluations, respectively. For all of these first three setups, the optimization target has no effect on the utility. This situation changes for setup (4), where we make it the main contribution to the utility function. Consequently, the evaluations are concentrated in an area with high values of $t$ and the lower left part of the topology is hardly explored. Finally, setup (5) shows a pure grid approach in the framework of our algorithm for comparison, which is finished after the initial sampling step since $G^{2}=N$.

3.2.1 Relative success rate and score

In C we define the relative success rate

[TABLE]

and the score

[TABLE]

as a quality measure for our estimator. The score represents the fraction of correct outcome predictions performed for all parameters in the explored data set $D_{\mathrm{expl}}$, whereas $R$ is defined as the fraction of correct outcome predictions for almost all parameters in the parameter space $\bm{\chi}$. In other words, $\sigma$ is a local and $R$ a global quality measure of our estimator.

The fraction

[TABLE]

consequently tells us how well the outcome predictions from the estimator obtained from the explored data set generalize to the whole parameter space. If $D_{\mathrm{expl}}$ is a representative subset of $\bm{\chi}$, $g$ can consequently be used to measure whether the estimator is overfitted ($g>1$) or underfitted ($g<1$). However, since the aim of our exploration algorithm is to find a representative subset of the parameter space in the first place, we can instead consider $g$ as an estimated quality measure for our training set itself. For $g>1$ the local prediction is better than the global prediction, hence the explored data can be seen as an oversimplified subset. For $g<1$, conversely, the local prediction is worse than the global prediction and the explored data can be seen as overcomplicated subset. If $g=1$, the local and global predictions are equally good. In this case the explored data set can be considered as perfectly representative.

The kernel SVM we use for classification allows us to almost always achieve a perfect score $\sigma$ and therefore we remain in the realm $g\geq 1$. Practice has shown that it seems to be a good approach to use such oversimplified data subsets during exploration. This observation could be due to the fact that our utility function is in such cases mostly based on the main features of the classification border and tends to neglect minor details, which are usually not important for all but the very last exploration steps.

3.2.2 Ratios of false positives and false negatives

Two additional global quality measures for our estimator are given by the ratio of false positives

[TABLE]

defined in D. We use the convention that positive results correspond to valid outcomes and negative results to invalid outcomes. Consequently, $r_{\mathrm{fp}}$ represents the fraction of wrongly predicted valid outcomes and $r_{\mathrm{fn}}$ the fraction of wrongly predicted invalid outcomes, respectively, performed for almost all parameters in the parameter space $\bm{\chi}$.

An estimator of high quality is indicated by a high relative success rate $R$, Eq. 22, a low ratio of false positives $r_{\mathrm{fp}}$ and a low ratio of false negatives $r_{\mathrm{fn}}$. Depending on the application, either false positives or false negatives might be considered far more adverse. However, if the two ratios are to be valuated equally, considering $R$ might be a sufficient quality measure since by definition

[TABLE]

holds true.

3.2.3 Validity ratio

The calculation of valid data points yields meaningful optimization targets and therefore provides us with more information than the calculation of invalid data points. Hence, the fraction of valid to invalid data points can serve as a measure for the usefulness of a sampling. For this purpose we have defined the validity ratio

[TABLE]

as the fraction of valid to invalid data points in the explored data set $D_{\mathrm{expl}}$ in E. Moreover, we have defined its reference limit

[TABLE]

as the fraction of the total volumes for valid and invalid outcomes, respectively, in the whole parameter space. The validity ratio $\alpha$ of a uniform random sampling in the whole parameter space will eventually converge to $\alpha_{\infty}$ for a sufficient number of samples. Therefore, this value represents a reasonable reference scale for the validity rate $\alpha$ which other sampling approaches have to be measured up to.

The calculation of $\alpha_{\infty}$ also allows us to specify a worst-case reference limit for the ratios of false positives and false negatives, Eq. 25. Specifically, we consider a completely randomized estimator that predicts valid and invalid outcomes with equal chance. Performing such random predictions for almost all parameters in $\bm{\chi}$ lead us to the reference limit for the ratio of false positives

[TABLE]

respectively as explained in E. An estimator exceeding these limits is consequently worse than a “random coin-tossing” estimator.

3.2.4 Exploration characteristics

We list the characteristic values of exploration $R$, Eq. 22, $\sigma$, Eq. 23, $r_{\mathrm{fp}}$, Eq. 25a, $r_{\mathrm{fn}}$, Eq. 25b, and $\alpha$, Eq. 27 together with the best and worst optimization targets

[TABLE]

respectively, for each of the explored data sets, (1) to (5), in table 2. In Eq. 30 we have recalled the extremal optimization targets from Eq. 17.

A Monte Carlo approach is used to calculate

[TABLE]

with $\mathcal{N}=$10,000$$ data points; see C. It is important to emphasize that global quality measures like $R$ are only possible because our considerations are not limited to predefined data sets, but rather make use of the fact that we can calculate new simulation data on demand. This enables us to create the randomly chosen samples necessary for the Monte Carlo approach.

3.2.5 Discussion

A comparison of the setups (2) and (5) in table 2 shows that our data exploration approach increases the highest relative success rate $R$ by almost $8\%$ in comparison with a regular grid sampling with the same number of sampling points. As expected, the highest value for $R$ is given by setup (3), followed by the setups (2) and (1). The worst performance is given by setup (4) and the grid approach, setup (5). However, while setup (4) has a poor value of $R$, it leads to the best optimization target $t_{\mathrm{best}}\approx$1.716, which almost reaches the theoretical limit of $\sqrt{3}\approx$1.732. This is no surprise given our choice of the weight vector $\mathbf{w}$, Eq. 11, which enforces parameter sampling near the best optimization target. Summarized, we see from our toy example that by tuning the weight vector we can intuitively control the exploration behavior.

The score $\sigma$ is perfect for all setups, which can be expected from such small training sets. Thus, we remain in the realm $g>1$, Eq. 24.

Only the setups (2) and (3) have ratios $r_{\mathrm{fp}}$ below the worst-case reference limit $r_{\mathrm{fp}\infty}$, Eq. 29a, while the ratios of all other setups exceed it. The best ratio is achieved for setup (3) and the worst for setup (4), which incorporates the optimization target. The ratios $r_{\mathrm{fn}}$, on the other hand, are almost the same for all setups and are all much smaller than the worst-case reference limit $r_{\mathrm{fn}\infty}$, Eq. 29b. We assume that this result is due to the fact that invalid outcomes are mostly found in the outer realm of the parameter space, a topological behavior that can be uncovered with almost any sampling method even with only a few samples because of the low dimensionality of the parameter space.

A comparison of the validity ratios $\alpha$ reveals that the best result is achieved for setup (4), followed by (3) and (1). Setups (2) and (5) even fail to beat the reference limit $\alpha_{\infty}$, Eq. 28. As we will see further below, these relatively bad validity ratios are a consequence of the very small number of samples considered here and will improve with ongoing exploration. The fact that setup (1) has a better validity ratio than setup (2) is also a result of the sparse sampling.

For an exploration approach with a weight $o=0$, for which the optimization target is ignored, a lower value of $t_{\mathrm{worst}}$ can be considered more favorable since it indicates a more complete parameter space exploration. On the other hand, if the optimization target is of importance by choosing a weight $o>0$, exploration of such uninteresting regions should rather be avoided and a higher value of $t_{\mathrm{worst}}$ can be considered more favorable. The latter is the case for setup (4) and we find that it shares the same value for $t_{\mathrm{worst}}$ with setup (1) and the grid approach, setup (5).

3.2.6 Summary

The toy example clearly shows how the weights in the utility function can be used to control the explorative behavior in an intuitive way. For an accomplished data exploration, we consider a high value of $R$, $\alpha$ and $t_{\mathrm{best}}$ as desirable results, whereas $r_{\mathrm{fp}}$ and $r_{\mathrm{fn}}$ should be small. In other words, we seek (i) a good outcome prediction while (ii) evaluations in regions with a good optimization score should be preferred and (iii) the explored data set has a good ratio of valid to invalid outcomes. In this sense, these five exploration characteristics can be seen as objectives of a MCO problem with possibly conflicting goals. Relative importance of these objectives varies by the application.

4 Benchmark

To demonstrate the strengths of our algorithm, we will in the following briefly present a benchmark with a Kriging-based exploration approach, which has been described in Refs. [5, 30]. We have already briefly mentioned this alternative strategy in section 1. Summarized, the Kriging-based exploration works in an iterative way similar to our novel method: Starting from an initial data set, a feasibility estimator is trained each step with the currently explored data points. However, in contrast to our algorithm, the estimator is continuous and therefore relies on a continuous feasibility function which reflects the degree of feasibility constraint violation. Based on the continuous estimator, a new parameter is suggested and the simulation is evaluated for the new parameter. The resulting data point is included in the set of explored data points and the next iteration begins.

A more detailed explanation of the Kriging-based exploration can be found in Refs. [5, 30]. Our algorithm is described in section 2. For the benchmark we will make use of our toy simulation, Eq. 21, from section 3.

4.1 Competing algorithms

The main purpose of the benchmark is to study the behavior of the two algorithms of interest – our proposed method and the Kriging-based approach – with respect to the number of sampled points $N$. Therefore, we run our algorithm on a collection of setups $(6)_{n}$ defined by the hyperparameters $G=4$, $\bm{\chi}_{\mathrm{G}}=\bm{\chi}_{\mathrm{toy}}$, $s=2$, $o=0$ and $r=1$. By setting $s>r$ we suppressed the spreading of data points in favor of a more precise feasibility border sampling. Each setup in the collection only differs by its number of samples $N=n$.

The notation of the setups follows section 3. Furthermore, we use an additional convention: The type of brackets we use to denote an exploration setup stands for the type of algorithm this setup is using to determine the sampling. Setups with round brackets refer to our algorithm, whereas non-round brackets represent different algorithms, as we will see in the following.

4.1.1 Kriging-based approach

For comparison, we run the Kriging-based exploration algorithm on a collection of setups $[7]_{n}$, where $n$ stands for the total number of sampled points in the same sense as for the setups $(6)_{n}$. For the estimator we use a Gaussian Process Regression (GPR) [22] and achieved the best results using a Matérn kernel [14] with smoothness parameter $\nu=$1.5$$. In each iteration the hyperparameters of the GPR kernel are optimized by maximizing the log-marginal-likelihood of the GPR model with the help of the L-BFGS-B algorithm [8]. The training data is standardized by removing the mean and scaling to unit variance. Because the toy simulation, Eq. 21, is defined as a binary classification problem and by design no explicit information about feasibility constraint violation is available, we use

[TABLE]

as target values for the GPR model. In each iteration step, a new parameter $\textbf{x}^{\mathrm{GPR}}_{\mathrm{new}}\in\bm{\chi}$ is sampled where the expected improvement [30]

[TABLE]

becomes maximal so that

[TABLE]

holds true. Here we have made use of the GPR model prediction $\hat{y}^{\mathrm{GPR}}(\textbf{x})$ and its corresponding standard error $s^{\mathrm{GPR}}(\textbf{x})$ evaluated for the parameter x. The initial data set for the training of the GPR model is chosen as a regular grid of $4\times 4$ data points in complete analogy to the initial data set for our algorithm. The numerical solution of Eq. 34 is obtained with a differential evolution approach.

4.1.2 Modified Kriging-based approach

The expected improvement, Eq. 33, is supposed to focus the sampling on the feasibility border where the predicted violation $\hat{y}^{\mathrm{GPR}}(\textbf{x})$ vanishes. However, in our binary classification scenario with the constraint violation described by Eq. 32, all points are penalized by either one of the two outcomes regardless of their distance to the border. This lack of knowledge about a continuous distance measure might significantly worsen the performance of the Kriging-based approach. Therefore, we suggest a modification which can compensate this deficiency. Specifically, we consider a reformulation

[TABLE]

of Eq. 32, where

[TABLE]

represents an upper bound of the closest Euclidean distance of the parameter $\mathbf{x}$ to the feasibility border. According to Eq. 36, this upper bound corresponds to the Euclidean distance of $\mathbf{x}$ to its nearest neighbor of opposing feasibility $\mathbf{x}^{\prime}$ in the currently explored data set $D_{\mathrm{expl}}$. In particular, additional samples in $D_{\mathrm{expl}}$ can only improve $\mathrm{d_{fb}}(\mathbf{x})$, which converges to a tight bound for the theoretical limit of an infinite number of samples.

Apart from the major difference of using $y^{\mathrm{GPR}}_{\mathrm{continuous}}(\mathbf{x})$ instead of $y^{\mathrm{GPR}}(\mathbf{x})$ we leave the Kriging-based algorithm unchanged. For comparison, we run this modified algorithm on a collection of setups $\{8\}_{n}$, where $n$ stands for the total number of sampled points in the same sense as for the setups $(6)_{n}$ and $[7]_{n}$.

4.2 Results

For the benchmark we compare the exploration characteristics introduced in section 3.2. Specifically, we consider the success rate $R$, Eq. 31, the score $\sigma$, Eq. 23, the ratio of false positives $r_{\mathrm{fp}}$, Eq. 25a, the ratio of false negatives $r_{\mathrm{fn}}$, Eq. 25b, and the validity ratio $\alpha$, Eq. 27, together with the best and worst optimization scores $t_{\mathrm{best}}$ and $t_{\mathrm{worst}}$, respectively, Eq. 30. The exploration characteristics are summarized in table 3 for four different numbers of samples $n\in\{25,50,100,150\}$ for each of the three candidate algorithms. Additionally, we show the progression of certain characteristics for an increasing number of samples $n\in[20,150]$ in Figure 3. The vertical lines ( ) correspond to the four chosen numbers of samples from table 3. Finally, snapshots of the explored outcomes for these samples can be found in Figure 4 analogously to the left column in Figure 2.

4.2.1 Discussion

The top plot in Figure 3 shows that after about 35 samples our algorithm ( ) steadily exceeds the success rate $R$ of the Kriging-based method ( ) by roughly $0.02$, and similarly, by roughly $0.01$ the success rate $R$ of the modified Kriging algorithm ( ), albeit after about 80 samples. Both Kriging models exhibit higher $R$ values initially, however at the cost of increased false positive rates $r_{\mathrm{fp}}$. The false positive rate is consistently better for our algorithm over all number of samples. Both Kriging-based algorithms exhibit a peak in the false positive rate in the initial sampling stage up, even up to the reference limit $r_{\mathrm{fp}\infty}$, Eq. 29a, indicating overestimation of the feasible range as illustrated in the top row of Figure 4. On the other hand, the third plot in Figure 3 shows that our algorithm tends to underestimate the feasible region, especially for the first 35 samples. However, the major difference between benchmarked algorithms is highlighted in the bottom plot in Figure 3. Our algorithm exceeds the reference validity ratio $\alpha_{\infty}$, Eq. 28 already after about 35 samples. After about 70 samples the validity ratio of our algorithm is twice as high as that of the Kriging-based algorithm, and about $4/3$ times higher than that of the modified Kriging algorithm. Such high validity ratio values indicate that majority of new points are placed within feasible region and thus relevant parameter space is sampled more efficiently. Without modification, the Kriging-based algorithm reaches the validity ratio of a uniform sampling, thus confirming that the binary constraint violation formulation in Eq. 32 does not penalize sufficiently sampling far from the feasibility border. The modified constraint violation function, Eq. 35, appears to partially address this issue as it results in an increased validity ratio, albeit it does not reach the level of our algorithm.

Our algorithm constantly pushes $t_{\mathrm{best}}$ and $t_{\mathrm{worst}}$ to their limits with an increasing number of samples, whereas the Kriging-based algorithm already reaches fixed values after about 25 sampled points. Those fixed values are far from the theoretical limits given by $\pm\sqrt{3}\approx\pm$1.732$$. This indicates that an exploration of the feasibility region is rather coarse with the Kriging-based approach. The modified Kriging algorithm reaches similar $t_{\mathrm{best}}$ and $t_{\mathrm{worst}}$ values indicating better sampling at the feasible boundary than the unmodified Kriging-based algorithm. Note that we have not included the optimization target in our utility function, which would lead to a better value of $t_{\mathrm{best}}$ for much less samples as shown in section 3.2. Moreover, the optimization target is not explicitly considered in the Kriging-based algorithms in the first place. The best optimization target can therefore be expected to be much worse if the target maximum is not directly located on the feasibility border.

From the progression of $R$ in Figure 3 we see that the prediction quality of all compared algorithms appears to be almost saturated for 150 samples. Although additional data points can still lead to an improvement of the estimators’ precision, this improvement is expected to be rather small and might not be worth the additional calculation effort. According to table 3 we have $\sigma<R$ so that $g<1$, Eq. 24, for $(6)_{150}$, which indicates that the local prediction is worse than the global prediction. Thus, the explored data set can be seen as an overcomplicated subset of the complete parameter space. Such a behavior could in fact be used as a possible quantitative stopping criterion for our algorithm in agreement with the qualitatively observed saturation.

4.2.2 Summary

The benchmark has shown that our method is comparable to the Kriging-based approach for a very small number of samples, but outperforms it for an increasing number of iterations. It is particularly remarkable that we achieve a much higher ratio of valid to invalid data points with our algorithm, which makes our sampling more efficient. Moreover, we can clearly consider our modification of the Kriging-based approach as an improvement for the case of a discrete feasibility constraint violation.

5 Application

In the current section we will apply the proposed algorithm from section 2 to the simulation of a realistic chemical process. This process describes the inner workings of a production plant consisting of two connected distillation columns. The task of the plant is to separate a two-component mixture consisting of chloroform and acetone by means of a so-called pressure swing distillation. Such kind of operations are very common in chemical engineering and their physical description is well-known [3].

As shown in Figure 5 each column has one input stream (feed) and two output streams (distillate and bottom stream), consisting of chloroform and acetone mixtures. Different parameters such as the operating pressure $P$ influence the concentrations of the components in the output streams. The two columns are connected in such a way that the bottom stream of the first column constitutes the feed for the second column and the bottom stream of the second column is recycled and mixed with the educt of the process to form the feed of the first column. The two distillate streams constitute the products of the process, which are desired to have high purities.

The distillation process we consider is not straightforward since the mixture of chloroform and acetone exhibits an azeotropic behavior. This means that at the azeotropic point, the composition in the vapor phase equals the composition in the liquid phase, which is sensitive to $P$ in this particular case. One possibility to surpass this distillation limit and to separate the azeotropic mixture is to use a special setup in which the first column operates at low pressure ($P=$1\text{,}\mathrm{bar}) and the second column operates at high pressure ($P=$10\text{\,}\mathrm{bar}). As a result, a high concentration of acetone in the distillate stream of the first column can be achieved and the concentrations in the corresponding bottom stream are close to the azeotrope. The higher pressure in the second column changes the azeotrope composition in such a way that a high concentration of chloroform can be achieved in the distillate stream of the second column.

5.1 Chemical simulation

The most common approach for modeling distillation columns is the equilibrium stage model (ESM) [3], which is based on a cascade of interconnected equilibrium stages. Every stage has a vapor (boiling) and a liquid (condensing) output stream, where the vapor stream of each stage raises to the stage above and the liquid output stream of each stage flows to the stage below. The reflux ratio of a column represents the ratio between the internal liquid streams of the last two stages and the distillation stream. A higher reflux ratio means that more energy is needed for cooling and heating. For every stage certain physically motivated equations (conservation of mass and conservation of enthalpy together with thermodynamic equilibrium and closure conditions) must hold. The ESM is used in all commercially available flow sheet simulators and is used wold-wide to simulate chemical distillation processes.

We model each column with an ESM consisting of $28$ stages (not shown in Figure 5). The non-random two-liquid (NRTL) model [23] is used to describe the interactions between the substances in each stage. The complete chemical process is then represented by a system of about $400$ coupled linearly independent equations, some of which are highly nonlinear. We also have a comparable number of independent internal variables describing, e. g., the concentrations, temperatures and flow rates of the internal streams between the stages.

As already mentioned in section 1, we use the flowsheet simulator Chemasim to perform the calculations. Given the parameters $\mathbf{x}$, we formally define the function

[TABLE]

to describe the outcome of a Chemasim evaluation in the classification space $\eta$, Eq. 3. Convergence and physical feasibility are exclusively decided by Chemasim-internal criteria and we have no knowledge about the degree of a violation, i. e., we have no access to a continuous feasibility function.

We consider a four-dimensional parameter space (i. e., $p=4$) for the chemical process simulation. It is constituted by the mass fractions of acetone at the distillate stream of column one $m_{\mathrm{ac}}\in[0.1,1.0]$ and chloroform at the distillate stream of column two $m_{\mathrm{cl}}\in[0.8,1.0]$, as well as the reflux ratios of the two columns $r_{1}\in[5,35]$ and $r_{2}\in[5,35]$, respectively. Using the parameter vector

[TABLE]

our parameter space can consequently be written as

[TABLE]

and therefore corresponds to a four dimensional hyperrectangle.

The simulation $\mathcal{S}_{\mathrm{sim}}$, Eq. 1, can formally be expressed by

[TABLE]

respectively, where we have recalled Eq. 37. Our optimization target $t(\mathcal{S}_{\mathrm{sim}},\mathbf{x})$ is chosen as the average of the two distillate mass fractions, which we aim to maximize. Recall that by definition, the optimization target is only meaningful for valid simulation outcomes $y(\mathcal{S}_{\mathrm{sim}},\mathbf{x})$. We specify sufficient internal variables of the simulation with a suitably chosen but fixed value so that given the parameters $\mathbf{x}$ the system of equations is well-defined and can be solved by Chemasim.

We assume that we have no previous knowledge about the data (i. e., $D_{\mathrm{init}}=\{\}$) except for the fact that there is a parameter space window

[TABLE]

in which both valid and invalid parameters can be found. For the start-up step of our algorithm we choose a regular grid which expands inside of a chosen hyperrectangle $\bm{\chi}_{\mathrm{G}}$. Eq. 9 is again solved numerically using a differential evolution approach.

5.2 Results

The results are shown in Figure 6 analogously to Figure 2. To achieve a two dimensional illustration of the four dimensional parameter space $\bm{\chi}_{\mathrm{sim}}$, Eq. 39, we show three different planar cuts through $\bm{\chi}_{\mathrm{sim}}$ in the first three columns and project all data points onto the respective planes. Specifically, we choose the planes for which

[TABLE]

holds true with $c=10$ in column ($\mathrm{a_{1}}$), $c=20$ in column ($\mathrm{a_{2}}$) and $c=30$ in column ($\mathrm{a_{3}}$), respectively. As a consequence, the outcomes of the projected data points do not necessarily have to correspond to the estimated feasibility borders ( ) on each plane. It is also important to emphasize that the exact theoretical feasibility border is unknown so we cannot plot it. For the optimization target visualization in column (b), the specific choice of the projection plane does not affect the plots.

Each of the five rows represents a different setup as summarized in table 4. Specifically, the first two setups are based on $N=40$ evaluations and while setup (1) ignores the optimization target, it is taken into account by setup (2) with an equal weight. The third and fourth setup also contain these two cases but for $N=81$ evaluations. All of these four setups have an initial grid parameter $G=2$ and make use of a grid in the parameter space window $\bm{\chi}_{\mathrm{win}}$, Eq. 41. Finally, the setups (5) and (6) represent pure grid approaches in the framework of our algorithm, which are finished after the initial sampling step since $G^{4}=N$. In table 4 we also list three additional setups, $\langle$7$\rangle$ to $\langle$9$\rangle$, which are not shown in Figure 6. Each of these setups represents a typical LHS of $N$ parameter points in $\bm{\chi}_{\mathrm{sim}}$, Eq. 39. The other exploration hyperparameters have no meaning for these setups.

As one would expect, it becomes apparent from columns ($\mathrm{a_{1}}$) to ($\mathrm{a_{3}}$) that the estimated feasible regions expand and include higher distillate mass fractions $m_{\mathrm{ac}}$ and $m_{\mathrm{cl}}$, respectively, with increased reflux ratios $r_{1}=r_{2}$. Interestingly, this observation can be made for all depicted setups. The only exception can be found in ($1\mathrm{a_{2}}$), which covers a smaller validity area than ($1\mathrm{a_{1}}$), but still includes higher distillate mass fractions. Since $r_{1},r_{2}\in[5,35]$, column ($\mathrm{a_{2}}$) in fact shows a projection onto the central reflux plane.

5.2.1 Exploration characteristics

Figure 6 serves as an illustration of the exploration behavior, however, to quantify the success of individual setups, we compare the exploration characteristics introduced in section 3.2. The characteristic values of exploration are given by the success rate $R$, Eq. 22, the score $\sigma$, Eq. 23, the ratio of false positives $r_{\mathrm{fp}}$, Eq. 25a, the ratio of false negatives $r_{\mathrm{fn}}$, Eq. 25b, and the validity ratio $\alpha$, Eq. 27. We show these values together with the best and worst optimization scores $t_{\mathrm{best}}$ and $t_{\mathrm{worst}}$, respectively, Eq. 30, in table 5 for all nine setups.

For the calculation of $R$, $r_{\mathrm{fp}}$, $r_{\mathrm{fn}}$ and the respective reference limits $\alpha_{\infty}$, $r_{\mathrm{fp}\infty}$ and $r_{\mathrm{fn}\infty}$ we use a Monte Carlo approach with $\mathcal{N}=$10,000$$ data points. Specifically, we make use of the success rate approximation, Eq. 31, the approximation of the ratio of false positives

[TABLE]

as explained in C and D, respectively. Moreover, the above-mentioned reference limits are approximated by

[TABLE]

[TABLE]

as explained in E and D, respectively. A numerical evaluation of Eq. 44 yields

[TABLE]

which indicates a close balance between valid and invalid points in the whole parameter space. Using Eq. 46 we can directly determine

[TABLE]

which are also very similar due to this balancing.

5.2.2 Discussion

We find that according to table 5 setup (3) reaches the best value for the relative success rate $R$, closely followed by setup (1). Both results have overlapping confidence intervals. A slightly worse value for $R$ is achieved by the setups (2) and (4). By contrast, remarkably bad results can be found for the grid and LHS approaches. A comparison from setup (6) with setup (3) reveals that our data exploration approach results in an improvement of $R$ by almost $90\%$ with only one third of the sampling points ($81$ instead of $256$). This result shows that our algorithm is way more efficient especially for higher dimensional spaces than a uniform grid sampling.

The scores $\sigma$ show that the estimators trained on the LHS samples also perform bad on the training set itself. Since $\sigma$ increases with a higher number of samples, we assume that this behavior might be a consequence of the sparsity of the sample in the four dimensional parameter space, which can lead to an overly complex prediction of the feasibility border. It might therefore be possible that using a different estimator may result in a better overall performance for the LHS approaches. For all setups, we have $g>1$, Eq. 24.

When we compare the ratios of false positives and false negatives for the first four setups, we find that setup (3) has the lowest value of $r_{\mathrm{fp}}$, followed by (1), (2) and (4). The lowest value of $r_{\mathrm{np}}$ is, on the other hand, achieved by (4), followed by (2), (1) and (3). It is remarkable that the setups (4), $\langle$8$\rangle$ and $\langle$9$\rangle$ exceed the reference limit $r_{\mathrm{fp}\infty}$, Eq. 47a, but have almost no false negatives, whereas the setups (5), (6) and $\langle$7$\rangle$ exceed the reference limit $r_{\mathrm{np}\infty}$, Eq. 47b, but have almost no false positives. Clearly, all of the associated estimators generalize very badly due to the unrepresentative sampling used for the training.

The validity ratio $\alpha$ is best for setup (4), followed by setup (2). Both of these ratios clearly exceed the reference limit $\alpha_{\infty}$, Eq. 46. The other setups, including the three LHS approaches, only have validity ratios smaller than $\alpha_{\infty}$ and can therefore be considered less useful than a typical random sampling. This result is no surprise since respecting the optimization score will guide the exploration towards parameter space realms with valid outcomes.

As expected, we see from a comparison between (1) and (2) or (3) and (4) in Figure 6 that omission of the optimization target leads to a more regular spreading of the evaluations in the parameter space. Again, it is no surprise that incorporating the optimization score can lead to a worse outcome classification.

The highest value for $t_{\mathrm{best}}$ is reached by setup (4), closely followed by setup $\langle$9$\rangle$. All other setups achieve slightly worse values. Therefore, we find that even a small number of evaluations provides us with a reasonably good optimization target. From a direct optimization of the simulation by means of a MISQP procedure [24] we find an optimization target of $t\approx 0.950$, which is only slightly better than our best exploration result of $t\approx 0.930$. Its calculation can, however, require a few hundred parameter evaluations. Moreover, the Chemasim-internal optimizer has difficulties solving this specific problem due to the fact that the optimization target, Eq. 40b, has vanishing derivatives with respect to $r_{1}$ and $r_{2}$, which makes them insignificant for a gradient-based optimization approach. The feasibility outcome, on the other hand, explicitly depends on those parameters. Since the optimizer ignores this relation, it can be very difficult to find an optimal target in Chemasim without injecting expert knowledge.

The value of $t_{\mathrm{worst}}$ shows us in how far regions of a rather uninteresting optimization score have been explored. Exploration of such regions might be necessary to improve $R$, but should be avoided in favor of more relevant regions if the optimization target is of interest. From a comparison between (3) and (4) we find, without surprise, that incorporating the optimization target into the utility function leads to a more favorable value of $t_{\mathrm{worst}}$. However, a comparison between (1) and (2) shows the opposite effect. We assume that this behavior is a result of an unfavorably trained estimator for the optimization target during one iteration of the exploration process.

5.2.3 Summary

Our data exploration method has proven to be clearly superior to grid-based or LHS approaches. From the calculation of comparably few data points we already get a very good estimation of the feasibility region in the whole parameter space. Moreover, by incorporating the optimization target in the utility function we can achieve a sampling which contains data points very close to the optimum. Such data points can serve as suitable starting points for a rigorous optimization algorithm.

6 Conclusions and Outlook

From our benchmark, we have found that our method yields better exploration characteristics than a previously suggested Kriging-based approach [5, 30] for a binary feasibility classification scenario as soon as a critical number of sampling points has been exceeded. We have also seen that the ratio of valid to invalid data points in the sampling is much higher with our strategy, which can be important if the data set should be further used, e. g., to train a shortcut model or for optimization purposes. The performance of the Kriging-based approach might have been worsened by the lack of knowledge about a continuous feasibility constraint violation in our example. Therefore, we have suggested an improvement of the original approach for such cases.

The chemical process simulation has shown that our data exploration method provides us with an excellent approximation of the data space topology from a relatively small number of data points in comparison with grid-based or LHS approaches. By tuning the explorative hyperparameters we can intuitively control the exploration behavior. In this way we can concentrate the exploration on parameter regions with a relatively good optimization target. Hereby, we have discovered an almost perfect optimization target which could be used as a very suitable starting point for rigorous optimization strategies to speed up the optimization process. Since the simulated chemical process is fully realistic and industrially relevant we have demonstrated that our method is applicable to a real-world problem.

It is important to emphasize that some industrially relevant applications certainly require more than four design parameters. In such cases both our data exploration method and the previously suggested Kriging-based approach might be challenged by the exploration of a high dimensional parameter space [30]. To circumvent this curse of dimensionality, we propose to couple the data exploration strategy with a dimensionality reduction method [29] in order to restrict the search to a lower dimensional manifold. However, an extended discussion of such an approach is beyond the scope of this manuscript.

Naturally, our method introduces an additional computational overhead in comparison with conventional data exploration strategies like a regular grid or a randomized sampling. Therefore, it is best suited for simulation environments where suggesting the next data point to evaluate takes significantly less time than the actual evaluation. Simulations of chemical processes constitute perfect candidates for this requirement due to their computational complexity and the difficulty of predicting their operation window.

As a stopping criterion for exploration we have used a fixed number of evaluations. This has allowed us to compare our method with grid-based and random sampling approaches. In practice one might instead want to make use of a suitable precision measure for the estimator to stop the evaluation at a sufficient accuracy. In the scope of our benchmark we have suggested to use the fraction of the relative success rate to the score as a possible stopping criterion comparable to a saturated exploration.

We have shown that the utility function strongly influences the outcome of our algorithm. It seems natural to ask in how far this function could be modified or generalized. In the following we will therefore briefly discuss open questions for future research related to this topic.

First of all, from a more general point of view, maximizing the utility function can also be regarded as a MCO problem, where each element of the utility vector corresponds to an objective function. In this sense the complete parameter space represents the feasible set of decision vectors. We use the weight vector to reduce this problem to a single objective problem. However, it could be insightful to treat the data exploration multicriterially. In a similar way it would also be possible to take a non-scalarized MCO target into account (e. g., in our application from section 5 both mass fractions could act as optimization targets). As indicated in the introduction, such a scenario is common in the context of chemical process engineering and therefore of great practical interest.

Furthermore, we always assume a constant weight vector for the entire duration of the exploration. One possible modification of our algorithm would be a dynamic approach in which the coefficients of the weight vector change in each step of the iteration. For example, an exploration which starts with a strong weight of the outcome prediction uncertainty and ends with a strong weight of the optimization target prediction would first focus on the classification boundary and later on the region of interest. Such a dynamic approach could also be combined with a MCO strategy.

In this manuscript we have focused on a binary classification approach to separate valid from invalid simulation outcomes. However, for certain applications it could be beneficial to distinguish different causes of invalidity. One could for example use a ternary classification approach to separate valid, physically unfeasible and numerically divergent outcomes. Depending on the simulation, even more different invalidity classes might be of interest. Our exploration method could be straightforwardly generalized to incorporate such an extension. SVMs could still be used to calculate prediction probabilities for such a multi-class classification problem [31]. Since the task of classification itself is still an active area of research with connections to many other scientific fields (such as quantum mechanics [26] with promising results [2, 13]), it can be expected that the quality and performance of classification methods is subject to future improvements.

In real-world applications the outcomes of a simulation may crucially depend on a large number of different simulation parameters (e. g., algorithmic parameters or initialization values) and their mutual interactions. This means that divergent outcomes could turn into feasible outcomes for different simulation parameters and vice versa. However, in a complex simulation environment the reasons for a numerical divergence can become practically untraceable. To take this behavior into account, it would be possible to make use of a stochastic perspective, where numerically divergent outcomes are only considered invalid with a certain probability. This invalidity probability would consequently reflect the ignorance about the influence of the simulation parameters on the outcomes.

Prediction probabilities for the (multi-class) classification of outcomes already provide a statistical framework that can be exploited in this context. It would for example allow to identify numerically divergent data points of high uncertainty, i. e., data points that result from a divergent simulation run but have a comparably low probability of belonging to the class of divergent outcomes. Such data points could then be re-evaluated (and possibly re-labeled if a different outcome occurs) using differently tuned simulation parameters. If the simulation environments allows it, a higher computational effort (e. g., by increasing the numerical precision or the iteration steps of the underlying equation solver) could be used for each revision in the hope that convergence can eventually be achieved. The allocation of suitable sample weights [18] would be an intrinsic way of SVMs to assign a higher certainty to re-calculated data points with the same divergent outcome.

Although the proposed statistical perspective fits rather naturally into our new method of optimized data exploration, we only consider it a conceptional idea that goes beyond the scope of this manuscript. Furthermore, the benefit of re-evaluating divergent data points harshly depends on the specific simulation environment and requires full control over the simulation parameters.

Summarized, we consider the straightforward variability of our algorithm through a modification of the utility function as a conceptional strength, which allows us to study different approaches in a single framework. Therefore, our method can also be seen as a very versatile starting point for further research in the fields of data exploration, feasibility classification and optimization. Moreover, its demonstrated practicability gives way for different kind of applications in the realm of chemical process engineering and beyond.

7 Acknowledgements

This work was developed in the Fraunhofer Cluster of Excellence “Cognitive Internet Technologies”. The authors would like to thank Christian Bauckhage, Jürgen Franke and Marius Kloft for their helpful and constructive comments. Our numerical examples were realized with the help of scikit-learn [19].

Appendix A Kernel methods

As explained in section 2.1.2 we use a kernel SVM and a kernel RR to perform predictions, Eq. 7. Both estimators rely on the kernel method [11]. Specifically, we assume that there exists a mapping

[TABLE]

from the parameter space $\bm{\chi}$ to a feature space $\mathcal{F}$ in such a way that the inner product

[TABLE]

represents a Gaussian kernel

[TABLE]

for all parameters $\mathbf{x}$ and $\mathbf{x}^{\prime}$ in $\bm{\chi}$, where $\gamma$ is a hyperparameter of the feature space metric and $||\cdot||_{2}$ stands for the 2-norm distance.

It is important to emphasize that the feature space for the SVM is not necessarily the same as for the RR. Therefore, we assume in fact two feature space mappings, namely $\phi_{\mathrm{C}}$ for the SVM, which maps to the classification feature space $\mathcal{F}_{\mathrm{C}}$, and $\phi_{\mathrm{R}}$ for the RR, which maps to the regression feature space $\mathcal{F}_{\mathrm{R}}$. Both mappings are defined in analogy to Eq. 48 and are associated with the Gaussian kernels $k_{\mathrm{C}}(\mathbf{x},\mathbf{x}^{\prime})$ and $k_{\mathrm{R}}(\mathbf{x},\mathbf{x}^{\prime})$, respectively, of the form given by Eq. 50 with hyperparameters $\gamma_{\mathrm{C}}$ and $\gamma_{\mathrm{R}}$, respectively.

Appendix B Classification feature space distance

The feature space distance between two parameters $\mathbf{x_{1}}$ and $\mathbf{x_{2}}$ with the common feature mapping $\phi$, Eq. 48, can be written as [25]

[TABLE]

based on an arbitrary kernel function $k$. We make use of this expression in section 2.1.3 to define the third component of the utility vector, Eq. 10, as the classification feature space distance

[TABLE]

between the parameter $\mathbf{x}$ and its nearest neighbor $\mathbf{x_{\mathrm{NN}}}$ from the data set $D_{x}$, Eq. 6. Here, $\phi_{\mathrm{C}}$ represents the classification feature space mapping of the kernel SVM estimator described in A. The nearest neighbor is obtained from solving

[TABLE]

The choice of a Gaussian kernel $k_{\mathrm{C}}(\mathbf{x},\mathbf{x}^{\prime})$, Eq. 50, allows us to rewrite Eqs. 52 and 53 as Eq. 18.

By definition, $U_{r}(D,\mathbf{x})$ reduces the utility the closer the parameter $\mathbf{x}$ is located to its nearest neighbor $\mathbf{x_{\mathrm{NN}}}$ in the classification feature space $\mathcal{F}_{\mathrm{C}}$, i. e., the more similar the parameter is to its nearest neighbor in $D_{x}$. For neighbors with identical features, $U_{r}(D,\mathbf{x})$ becomes [math] and for neighbors with fundamentally different features it can asymptotically approach $1$.

Appendix C Relative success rate and score

To quantify the quality of the exploration we make use of the relative success rate

[TABLE]

It contains the volume of correct outcome predictions

[TABLE]

with the indicator function

[TABLE]

and the total volume of the parameter space

[TABLE]

We can also write

[TABLE]

where the score

[TABLE]

of a test data set $D_{\mathrm{t}}$ represents the fraction of correct outcome predictions performed for all parameters $\mathbf{x}\in D_{\mathrm{t}}$. The absolute values in Eq. 59 denote a set cardinality, a notation which we will also use in the following sections. The specific test data set

[TABLE]

from Eq. 58 is required to contain data points $d(\mathbf{x})$, Eq. 4, for all or almost all parameters $\mathbf{x}$ from $\bm{\chi}$.

Summarized, the relative success rate $R(D)\in[0,1]$ is a quality measure for an outcome estimator taken across the whole parameter space, whereas the score $\sigma(D_{\mathrm{t}})\in[0,1]$ may also be used for a local quality measure depending on the chosen test data set $D_{\mathrm{t}}$. A larger value of $R(D)$ and $\sigma(D_{\mathrm{t}})$ indicates an estimator of higher quality.

In practice, the relative success rate can be determined by evaluting the integral in Eq. 55 with the help of a standard Monte Carlo approach [21], which also provides us with an estimated error. Given the collection of evaluated parameters

[TABLE]

in the sense of Eq. 6, which have been chosen randomly from $\bm{\chi}$, we can make use of Eq. 31 with the approximation

[TABLE]

and its estimated error

[TABLE]

respectively. Here we have introduced the Monte Carlo expectation value

[TABLE]

and the Monte Carlo error estimate

[TABLE]

which can both be considered a functional with respect to a function $f:\bm{\chi}\rightarrow\mathbb{R}$. The two expressions additionally depend on a collection of evaluated parameters $D_{x}$, Eq. 6. The symbol $z_{\alpha_{\mathcal{P}}/2}$ refers to the standard $\alpha_{\mathcal{P}}/2$ normal quantile associated with the probability $\mathcal{P}=1-\alpha_{\mathcal{P}}$ of the precise integration result with the approximation $E^{\mathrm{MC}}[f;D_{x})$ being within the confidence interval $[E^{\mathrm{MC}}[f;D_{x})-S^{\mathrm{MC}}[f;D_{x},\alpha_{\mathcal{P}}),E^{\mathrm{MC}}[f;D_{x})+S^{\mathrm{MC}}[f;D_{x},\alpha_{\mathcal{P}})]$. For all of our numerical examples we set $\mathcal{P}=95\%$ so that $\alpha_{\mathcal{P}}=5\%$ and |z_{0.025}|\approx1.960$$. We introduce the appropriate abbreviation

[TABLE]

to simplify our notation.

Appendix D Ratios of false positives and false negatives

A common quality measure for a binary classification method is the number of false positives (for which the classification method improperly predicts a positive result) and the number of false negatives (for which the classification method improperly predicts a negative result), respectively. Following this concept, we can define the ratio of false positives

[TABLE]

for the estimator $\mathcal{E}$ trained with respect to a test data set $D_{\mathrm{t}}$. By definition, a positive result corresponds to a valid outcome, whereas a negative result corresponds to an invalid outcome.

If we choose the specific test set $D_{\chi}(\mathcal{S})$, Eq. 60, we find

[TABLE]

respectively, in analogy to Eq. 58. Here we have recalled total parameter space volume $V_{\mathrm{total}}(\bm{\chi})$, Eq. 57. Moreover, we have introduced the volume of false positives

[TABLE]

which are based on the indicator functions

[TABLE]

respectively.

Summarized, the ratios $r_{\mathrm{fp}}(D)\in[0,1]$ and $r_{\mathrm{fn}}(D)\in[0,1]$ can be considered as a quality measure for an outcome estimator taken across the whole parameter space. A worst-case reference value for the ratios of false positives and false negatives will be introduced in E based on evaluating randomized predictions for almost all parameters in $\bm{\chi}$. Note that an estimator of high quality is indicated by a large relative success rate $R(D)$, Eq. 54, but small ratios of false positives and false negatives; also see Eq. 26.

In complete analogy to C, we can use a Monte Carlo approach to numerically calculate $V_{\mathrm{fp}}(\mathcal{S},\mathcal{E},D,\bm{\chi})$ and $V_{\mathrm{fn}}(\mathcal{S},\mathcal{E},D,\bm{\chi})$ in Eq. 68 using a collection $D^{\mathrm{MC}}_{x}$ of $\mathcal{N}$ evaluated parameters, Eq. 61, which have been chosen randomly from $\bm{\chi}$. This leads us to Eq. 43 with the approximations

[TABLE]

as well as their estimated errors

[TABLE]

respectively. Here we have recalled the Monte Carlo expectation value, Eq. 64, and the Monte Carlo error estimate, Eq. 66, to simplify our notation.

Appendix E Validity ratio

In general, valid data points are more insightful than invalid data points because they contain information about the optimization target. To quantify the usefulness of a sampling, we therefore introduce the validity ratio

[TABLE]

The higher the fraction $\alpha(D)$ of valid to invalid data points in the data set $D$, the more useful the corresponding sampling.

However, to provide a clear scale for $\alpha$, some kind of reference value is necessary. Such a reference value can be obtained from the idea of an infinite sample: For a given simulation $\mathcal{S}$, a fully dense sampling in the parameter space $\bm{\chi}$ with infinitely many data points results in the infinitely large data set $D_{\infty}(\bm{\chi})$. For this reference data set the validity rate converges to the reference limit

[TABLE]

given by the fractions

[TABLE]

of the total volumes for valid and invalid outcomes in $\bm{\chi}$, respectively. Here we have introduced the valid volume

[TABLE]

with the indicator function

[TABLE]

and the invalid volume

[TABLE]

based on the total parameter space volume $V_{\mathrm{total}}(\bm{\chi})$, Eq. 57.

For the two-dimensional toy example from section 3 the reference limit $\alpha_{\infty}(\mathcal{S}_{\mathrm{toy}},\bm{\chi}_{\mathrm{toy}})$, Eq. 75, can be calculated explicitly using the definitions of the parameter space $\bm{\chi}_{\mathrm{toy}}$ and the toy simulation $\mathcal{S}_{\mathrm{toy}}$, Eqs. 20 and 21, respectively. It is straightforward to see from Figure 1 that the valid volume

[TABLE]

corresponds to the area of two quarters of a circle of radius $\sqrt{3}$. On the other hand, the total parameter space volume

[TABLE]

is equivalent to the area of a square of side length $4$. According to Eq. 78 the invalid volume

[TABLE]

is the difference of the two aforementioned volumes. Thus, the reference limit

[TABLE]

is given by Eq. 28.

For the chemical process simulation from section 5, an explicit expression of the reference limit of the validity rate can not be obtained since we use a flowsheet simulator to perform calculations. Therefore, we have to fall back to a numerical approximation using a Monte Carlo approach, Eq. 44. We will discuss this method further below in more detail.

Our previous considerations also allow us to define a reference value for the ratio of false positives and false negatives, Eq. 68, from D in a straightforward way. Suppose a completely randomized estimator, who predicts valid and invalid outcomes with the same chance. For a fully dense sampling $D_{\infty}(\bm{\chi})$, the ratios with respect to this randomized estimator converge to the worst-case reference limits

[TABLE]

respectively. One has

[TABLE]

so that $r_{\mathrm{fp}\infty}+r_{\mathrm{fn}\infty}=\frac{1}{2}$ in general and $r_{\mathrm{fp}\infty}=r_{\mathrm{fn}\infty}=\frac{1}{4}$ for $\alpha_{\infty}=1$, as expected for random outcome predictions.

For the toy example from section 3, these terms can be evaluated using the explicit expression of $\alpha_{\infty}$, Eq. 82, which results in Eq. 29. In contrast, we only have an approximation of $\alpha_{\infty}$, Eq. 44, for the chemical process simulation from section 5 and can therefore only use the approximations given by Eq. 45. This approach is explained in the following.

Analogously to the considerations from C and D, a Monte Carlo approach allows a numerical evaluation of $V_{\mathrm{valid}}(\mathcal{S},\bm{\chi})$ in Eq. 75 when we presume a collection $D^{\mathrm{MC}}_{x}$ of $\mathcal{N}$ evaluated parameters, Eq. 61, which have been chosen randomly from $\bm{\chi}$. As a result, we find Eq. 44 with the approximation

[TABLE]

and its estimated error

[TABLE]

respectively. Here we have made use of the abbreviations

[TABLE]

and

[TABLE]

which are based on the Monte Carlo expectation value and the Monte Carlo error estimate, Eqs. 64 and 66, respectively.

Using these results, we can also numerically evaluate Eq. 84. Specifically, one has Eq. 45 with an approximation only depending on Eq. 85. The respective estimated error is given by

[TABLE]

and depends on Eqs. 85 and 86, respectively.

Summarized, for the chemical process simulation from section 5 the evaluation of Eq. 44 based on Eqs. 85 and 86 leads us to Eq. 46. Moreover, from the evaluation of Eqs. 45 and 89 we find Eq. 47.