Optimal experimental design: from design point to design region

Bubel, Martin et al. Article — Published Version Optimal experimental design: from design point to design region Statistical Papers Provided in Cooperation with: Springer Nature Suggested Citation: Bubel, Martin et al. (2025) : Optimal experimental design: from design point to design region, Statistical Papers, ISSN 1613-9798, Springer, Berlin, Heidelberg, Vol. 66, Iss. 5, https://doi.org/10.1007/s00362-025-01725-7 This Version is available at: https://hdl.handle.net/10419/323273 Standard-Nutzungsbedingungen: Die Dokumente auf EconStor dürfen zu eigenen wissenschaftlichen Zwecken und zum Privatgebrauch gespeichert und kopiert werden. Sie dürfen die Dokumente nicht für öffentliche oder kommerzielle Zwecke vervielfältigen, öffentlich ausstellen, öffentlich zugänglich machen, vertreiben oder anderweitig nutzen. Sofern die Verfasser die Dokumente unter Open-Content-Lizenzen (insbesondere CC-Lizenzen) zur Verfügung gestellt haben sollten, gelten abweichend von diesen Nutzungsbedingungen die in der dort genannten Lizenz gewährten Nutzungsrechte. Terms of use: Documents in EconStor may be saved and copied for your personal and scholarly purposes. You are not to copy documents for public or commercial purposes, to exhibit the documents publicly, to make them publicly available on the internet, or to distribute or otherwise use the documents in public. If the documents have been made available under an Open Content Licence (especially Creative Commons Licences), you may exercise further usage rights as specified in the indicated licence. http://creativecommons.org/licenses/by/4.0/ REGULAR ARTICLE Statistical Papers (2025) 66:109 https://doi.org/10.1007/s00362-025-01725-7 Abstract Optimal experimental designs are used in chemical engineering to obtain precise mathematical models. The optimal design consists of design points with a maximal amount of information and thus lead to more precise models than statistical designs. In general, the optimal design depends on an uncertain estimate of unknown model parameters θ . The optimal designs are therefore also uncertain and continuously shift in the design space, as the value of θ changes. We present two approaches to capture this behavior when computing optimal designs, a global clustering approach and a local approximation of the confidence regions. Both methods find an optimal design and assign the optimal design points confidence regions which can be used by an experimenter to decide which design points to use. The clustering approach requires a Monte Carlo sampling of the uncertain parameters and then identifies regions of high weight density in the design space. The local approximation of the confidence regions is obtained via an error propagation using the derivatives of the optimal design points and weights. We apply the introduced approaches to mathematical examples as well as to an application example modeling vapor-liquid equilibria. Keywords Confidence regions · Robust optimization · Optimal experimental design Received: 12 July 2024 / Revised: 12 May 2025 © The Author(s) 2025 Optimal experimental design: from design point to design region MartinBubel1· PhilippSeufert1· GlebKarpov2· JanSchwientek1· MichaelBortz1· IvanOseledets2,3 Martin Bubel [email protected].de 1 Fraunhofer Institut für Technound Wirtschaftsmathematik, Fraunhofer-Platz 1, D-67663 Kaiserslautern, Germany 2 Skoltech, Moscow, Russia 3 Artificial Intelligence Research Institute AIRI, Moscow, Russia 1 3 M. Bubel et al. 1 Introduction In modern industry, decision-making has increasingly become model-based, particularly in the process industry, where mathematical models are employed for the design and optimization of processes. These models, whether fully empirical or physically motivated, typically contain unknown parameters that require calibration to align with the actual system. Model calibration relies on experimental data gathered from wet-lab experiments or simulations such as computational fluid dynamics (CFD) simulations. Although models are often utilized for optimization - such as minimizing energy costs while maintaining throughput - the critical task of model calibration must precede this. The term experimental design refers to a series of experiments conducted specifically for model calibration. The effectiveness of an experimental design directly impacts the reliability and uncertainty of the calibrated model Bubel et al. (2024). Given that experiments can be costly, there is a significant interest in optimal experimental designs (OED) specifically tailored for model calibration. While the computation of optimal experimental designs is straightforward for models that are linear in their parameters, see e.g. Fedorov and Leonov (2013), for nonlinear models, approximations are required, as there is no general closed form expression for model uncertainty Bubel et al. (2025). A common approach for computing optimal experimental designs for nonlinear regression models involves linearization of the model by differentiation w.r.t. the model parameters, which creates a dependency on those parameters’ values. Consequently, it is often more efficient to conduct experiments and calibrate the model iteratively, utilizing sequential updates on model parameters and the resulting optimal experimental designs. This workflow, as depicted in Fig. 1, begins with existing data or a set of initial experiments, followed by iterative model calibration on the available data and the computation of new optimal experimental designs, and has been applied in multiple works, e.g. Atkinson (2008); Vanaret et al. (2020); Bubel et al. (2024). The transition from linear to nonlinear models is thoroughly reviewed in Fedorov and Leonov Fig. 1 Iterative workflow for model validation and adjustment 1 3 109 Page 2 of 28 Optimal experimental design: from design point to design region Fedorov and Leonov (2013). For linear models, the optimal experimental design remains independent of the reference parameter value. In contrast, for nonlinear models, the optimal design is influenced by these parameters, necessitating the discussion of model-based optimal experimental designs. Many methods for computing those, e.g. those of Yang and Stufken Yang and Stufken (2011) or Duarte et al. Duarte et al. (2018), are grounded in the Kiefer-Wolfowitz equivalence theorem Kiefer and Wolfowitz (1960), which outlines necessary and sufficient conditions for OED optimization. Traditional algorithms for computing optimal designs include the vertex direction method (VDM) introduced by Fedorov Fedorov and Leonov (2013) and Wynn Wynn (1970), the vertex exchange method (VEM) Böhning (1986), and the multiplicative algorithm (MUL) Pukelsheim and Torsney (1991); Silvey et al. (1978). Recent adaptations have aimed to enhance convergence and reduce runtimes, giving rise to methods such as the cocktail algorithm Yu (2011), the Yang-Biedermann-Tang algorithm Yang et al. (2013), the Weighted-Discretization-Approach Vanaret et al. (2020), MaxVol Mikhalev and Oseledets (2018); Vanaret et al. (2020), and a randomized exchange algorithm Harman et al. (2020). When dealing with nonlinear regression models, the dependence of optimal experimental designs on reference parameter values, which may be uncertain, introduces uncertainty into the designs themselves. To address this, the objective function in the optimization of experimental design is modified to account for uncertainty. Two predominant approaches are the average case approach Fedorov and Leonov (2013); Asprey and Macchietto (2002), which replaces the objective function with its expected value concerning the uncertain parameters, and the worst case approach Körkel et al. (2004); Asprey and Macchietto (2002); Duarte and Wong (2014), which minimizes the worst case scenario, also known as the minimax approach. Several algorithmic schemes for the worst case approach have been developed Asprey and Macchietto (2002); Körkel et al. (2004); Duarte and Wong (2014), while recent work has explored the conditional value at risk in OED contexts Valenzuela et al. (2015); Kusumo et al. (2022). There are also methods that combine expectation and variance, as presented by Mesbah and Streif (2015). Alternative strategies to manage uncertainty in OED include a method proposed by Gottu Mukkula and Paulen (2019), which uses Monte-Carlo simulations to compute the confidence region of the parameter estimate ˜ θ , providing a more precise estimate compared to linearization techniques using Fisher information matrices. Additionally, multi-stage designs introduced by Gottu Mukkula et al. (2021) involve performing a subset of experiments before observing the parameter value θ , with the remaining experiments planned based on the assumed knowledge of θ . Dror and Steinberg Dror and Steinberg (2006) present an approach that computes several optimal designs for varying parameter values, subsequently applying a clustering algorithm to obtain optimal experiments. This method was adapted by Biedermann and Woods Biedermann and Woods (2011), with a summary available in Chapter 14.6 of the Handbook of Design and Analysis of Experiments Biedermann and Yang (2015). In this paper, we introduce two additional approaches to manage uncertainty in OED. These methods are founded on the concept of replacing optimal design points xi with optimal regions Xi⊂X within the design space. The first approach utilizes clustering to average optimal designs across several sampled scenarios, identifying regions in the design space with high weight 1 3 Page 3 of 28 109 M. Bubel et al. density. The second approach examines sensitivities in the optimal design relative to the uncertain parameters θ . This method offers the advantage of requiring only local information about θ . By employing a linear approximation, we can assign each optimal design point a confidence ellipsoid, which serves as supplementary information for determining which experiments to conduct. 2 Optimal experimental design In this section we give a brief introduction into the theory of optimal experimental design. We focus on the main definitions and results and also briefly discuss optimal experimental design under uncertainty. The definitions and results mainly follow Chapters 1 and 2 of Fedorov and Leonov (2013), but can also be found similarly in most OED literature Atkinson (2008); Duarte and Wong (2015); Vanaret et al. (2020). In design of experiments we consider a parameterized model function f(x, θ) , which maps from a design space X⊂Rd X to an output space Y⊂Rd y, with θ∈Θ⊂Rd are the unknown model parameters from the set Θ , i.e. f: X ×Θ→Rd y, ( x, θ )→ f ( x, θ ). (1) To estimate the unknown parameters we perform experiments for design points x1,...,x n and take measurements y1,...,y n . We assume that the existence of a set of true parameters θt such that the observations are given by yi=f(xi,θ t)+εi∈Rd y, where the εi∼N(0,Σ) are random independent identically distributed (i.i.d.) measurement errors with covariance matrix Σ∈Rd y ×d y. 2.1 Locally optimal experimental designs A design (or experimental plan) ξN for the model f consists of design points xi∈X and repetitions ri of the experiments, with ∑ri=N the total number of experiments. Alternatively we can also write the number of repetitions via p i =ri N . We observe that the values pi sum to 1, and the design ξN thus is a (discrete) probability measure on the design space X. As is typical in OED we relax the notion of a design and define a design measure ξ as a probability measure on X. We introduce Ξ(X) as the space of all design measures. For any finite discrete design measures ξ=∑N i=1 ω i δ xi , given by a finite set of design points x1,...,x N∈X and weights ω1,...,ω N , the support of ξ is given by supp (ξ)={x1,...,x N}. (2) The efficiency of an experimental design ξ for model calibration is equivalent to the uncertainty of the least-squares estimator θLSE . For nonlinear regression models, the uncertainty of the least-squares estimator is not known in closed form, but can be approximated using the Fisher information matrix, which is defined as 1 3 109 Page 4 of 28 Optimal experimental design: from design point to design region m( x, θ )= D θ f ( x, θ )T·Σ−1· D θ f ( x, θ ), (3) where Σ denotes the covariance matrix of the measurement errors and Dθf(x, θ) denotes the Jacobian matrix of the model f with respect to the model parameters θ . Consequently, the covariance the of least-squares estimator θLSE depends on the experimental design ξ as Cov θLSE ≈ (M(ξ,θ))−1= (∫X m(x, θ)ξ(dx) )−1 . (4) Commonly, the uncertainty of the least-squares estimator is quantified in terms of a suitable scalar function Ψ , usually referred to as design criteria, see for instance Example 2.1 in Fedorov and Leonov (2013). In this work, we use the log-D-criterion ΨD , which is defined as Ψ logD = log (det ( M −1)). (5) While, from a theoretical perspective, the optima of the log-D and the D-criterion are the same, we prefer the log-D-criterion for numerical reasons. The interested reader is referred to the work of Ament et al Ament et al. (2024), which features a rigorous comparison of objective function and their log-based alternatives. Using the approximation of parameter uncertainty (4) and the scalar mapping (5), the resulting optimization problem for OED is given by ξ∗= arg min ξΨ logD ( M ( ξ, ˜ θ )). (6) In (6), ˜ θ denotes the reference value for which we compute the Fisher information matrix M(ξ, ˜ θ) . Given a parameter estimate ˜ θ , the solution to (6) is considered a locally optimal design, as it is only valid for θ=˜ θ . 2.2 Optimal experimental designs under uncertainty In the previous section we have considered locally optimal designs, which depend on a reference parameter value ˜ θ . For the task of model calibration, the reference parameter value usually is chosen as the least-squares estimate ˜ θ=θLSE , which usually is uncertain. If, θLSE is far from the true parameter value θt , local optimal experimental designs can result in non-optimal design for the true parameters θt , resulting in a less efficient model calibration process. An example of this has been shown for instance by Körkel et al. in Körkel et al. (2004). Hedging against uncertainty in the reference parameter value includes the consideration of uncertainty in local optimal experimental designs. Two classical approaches to handle uncertainty in OED are the average caseand the worst case approach, presented e.g. in Fedorov and Leonov (2013); Asprey and Macchietto (2002); Duarte and Wong (2014, 2015). The average case design is the solution to 1 3 Page 5 of 28 109 M. Bubel et al. arg min ξ Eθ [Ψ( M ( ξ,θ ))] , (7) while worst-case optimal design are obtained by solving arg min ξ sup θ∈Θ Ψ( M ( ξ,θ )). (8) For both approaches, under mild assumptions, an equivalence theorem can be formulated similar to those for locally optimal designs, see for instance Chapter 2.9 in Fedorov and Leonov (2013) for the average case problem and Asprey and Macchietto (2002); Duarte and Wong (2014) for the worst case problem. Comparing both approaches, we find that the average case approach is known to be too optimistic as it tends to ignore single bad scenarios. This is particularly problematic, if the true scenario θ=θt is one of those bad scenarios. In contrast, the worst case approach tends to be too pessimistic as it only focuses on the worst case scenario, generally leaving too much room for improvement in favorable scenarios. While these approaches have all been successfully applied to optimal experimental design, they suffer from two problems. First, the mentioned approaches give little information on the individual scenarios and it is somewhat unclear, how an optimal design will perform in the unknown true scenario. Second, the approaches all require optimal designs in (almost every) parameter value θ , or similar distribution information, to be evaluated. This information can be computationally expensive to obtain. That’s why we propose two alternative approaches to handle the uncertainty in the model parameters θ in this work, to arrive at a better understanding of uncertainty contained in experimental designs. Additionally, note that while we focus on the uncertainty in the model parameters θ throughout this work, the model function f may contain further uncertain parameters. The presented methods can easily be extended to such parameters. 3 Computing optimal regions in the design space In this section we introduce two alternative approaches to handle the uncertainties of the parameters θ . Instead of focusing on the design criterion Ψ and applying an uncertainty approach to the objective function of the optimization problem, we want to focus on the optimal design points x∗ i(θ) . A similar approach has been presented in Dror and Steinberg (2006) by Dror and Steinberg and later been modified by Biedermann and Woods in Biedermann and Woods (2011). They propose a Monte-Carlo sampling of optimal designs in the design space and subsequent clustering algorithms to identify a robust design. In the following we assume the uncertain parameters θ to be random variables with a known probability distribution. This can for example be a prior distribution, as is used for Bayesian design criteria. For the two approaches developed in this work, we state the following assumptions Assumption 1 (i) The design space X is compact. 1 3 109 Page 6 of 28 Optimal experimental design: from design point to design region (ii) The mapping x→ m(x, ˜ θ) is continuous for the reference value ˜ θ∈Θ . (iii) There exists a design ξ∈Ξfin( X, ˜ θ ) for the reference value ˜ θ∈Θ . If the above assumptions hold for all θ∈Θ , we can find a finite discrete optimal design ξ∗(θ) for each value of θ . The design points x∗ i(θ) and weights ω∗ i(θ) are also random variables, as they depend on θ . We expect the design points to smear in the design space, as the value of θ changes. The approaches presented in this chapter aim at identifying these regions in the design space. As each region is associated with a random variable x∗ i(θ) , we want to approximate the confidence region of each design point x∗ i(θ) . Computing the exact confidence regions is overly ambitious, as the regions may overlap, and the number of design points in an optimal design can change in the value of θ . An exact representation can only be found under strong mathematical assumptions. However, a simple approximation is possible and gives the necessary insight, into how the uncertainty of θ effects each design point x∗ i(θ) and the weights ω∗ i(θ) , and by extension the optimal design ξ∗(θ) . In both presented approaches we discuss approximations for the mean µ and the covariance matrix Σ assigned to the design points x∗ i(θ) . The mean µ can then be used as representative experiment, whereas the covariance matrix gives information on the uncertainty assigned to the experiment. If the design point x∗ i(θ) is normally distributed, its confidence region is an ellipsoid spanned by the matrix Σ with center µ . For a level α∈[0,1] the confidence ellipsoid is described by the equation Cα= y ( y − µ )TΣ−1( y − µ )≤ χ 2 1−α( k ), (9) where χ2 1−α( k ) denotes the 1−α quantile of the χ -squared distribution with k degrees of freedom. In general, the design points will not be normally distributed. However, the ellipsoid Cα can be used to approximate the actual confidence regions, or help interpret the size of the uncertainty given by the covariance matrix. 3.1 Clustering of design points As first approach to identify the confidence regions in the design space we propose a Monte-Carlo based sampling approach. Afterwards we apply a clustering algorithm to identify the regions in the design space. In Dror and Steinberg (2006); Biedermann and Woods (2011) Dror, Steinberg, Biedermann and Woods also propose clustering approaches to identify robust designs from a set of sampled optimal designs. They suggest the k-means clustering algorithm, which identifies k clusters among a set of design points. We use an alternative clustering algorithm. We sample the uncertain parameters θ and compute optimal designs ξj=ξ∗(θj) for each sampled scenario θ1,...,θM . By using the equivalence theorem Kiefer and Wolfowitz (1960) we assume the designs ξj to be discrete designs with a bounded number of support points. With the optimal designs, we can compute the average design 1 3 Page 7 of 28 109 M. Bubel et al. ξ =1 M· M ∑ j =1 ξj . (10) This design is again a design measure, which may however consist of a vast amount of support points, each with very small weight. Such a design in general cannot be used to identify experiments in applications. Using the clustering algorithm, we find regions in the design space with a high average weight, or equivalently, regions with a high weight density. In order to identify these regions we have found the algorithm DBSCAN Ester et al. (1996); Schubert et al. (2017) to be a sensible choice. This algorithm takes a list of points with corresponding weights and identifies clusters in regions with high weight density. We briefly discuss how the algorithm works. The algorithm depends on two parameters, a search radius r>0 and a minimum weight ωmin . For each design point xi the algorithm computes the total weight found in a sphere of radius r around this point xi . If this weight exceeds the threshold ωmin , the point is assigned as cluster point. If the total weight is less then ωmin , the point is discarded. The algorithm is well documented and previously published, thus we refer to the literature, for instance Ester et al. (1996); Schubert et al. (2017), as well as the implementation found in scikitlearn Pedregosa and Varoquaux (2011) for details. Compared to the k-means algorithm used by Dror and Steinberg Dror and Steinberg (2006), the DBSCAN algorithm has the advantage, that it does not require the number of clusters to be specified. Additionally, it also considers the optimal weights when identifying the clusters, whereas the k-means algorithm only works on the support points, thereby assigning each point the same significance. After performing the clustering we obtain a number of clusters cj , each given by a list of indices I c j={ i 1 ,...,i N j} , such that the cluster cj consists of points xi,i∈Icj, and corresponding weights ωi,i∈Icj . These clusters describe regions with high weight-density. In average the optimal designs ξ1,...,ξM allocate weight to these regions and thus it appears sensible, to perform an experiment from within these regions. For each cluster cj we compute some significant parameters. ●First, we compute the total weight of the cluster, given by ω( cj )= ∑ i ∈ Ic j ωi . (11) This value indicates the importance of the cluster with respect to the design ξ . The more weight a cluster has, the higher the priority for performing experiments from this cluster. ●Then, we consider the weighted mean of the cluster, 1 3 109 Page 8 of 28 Optimal experimental design: from design point to design region In this second region there are many points with small weight and there is no clear point to prioritize. We now apply the clustering algorithm with a search radius of r=0.035 and a minimal weight threshold of ωmin =0.1 . The algorithm finds two clusters. The corresponding cluster means, weights and confidence intervals are plotted in Fig. 3aand the clusters are also given given in Table 2. The clusters have a total weight of 0.955. Some design points with positive average weights corresponding in sum to 0.045 are not covered by the two clusters. Last we use the implicit function theorem to compute the derivatives Dθx∗ i and Dθω∗ i and then approximate the confidence intervals. For the derivatives we obtain Dθx∗ 2≈(0,0)T , Dθω∗ 1≈0 and Dθω∗ 2≈0 . For the point x∗ 1 we compute Dθx∗ 1≈(0,0.0625)T . These all correspond to the analytical derivatives, which are 0 for the constant weights and the point x∗ 2 and Dθx∗ 1= (0,θ−2 2) . Via the assumed distribution of θ we approximate the confidence intervals. The results for the optimal design points are plotted in Fig. 3b We observe, that both the clustering approach as well as the approximation of the confidence regions correctly depict the effect of the uncertain parameters. Both methods identify the optimal point x∗ 2=1.0 , which does not change in the parameters θ . This behavior is mirrored in both methods, as both confidence intervals have length 0. Also the second optimal design region, centered around x∗ 1=0.75 , is identified. The corresponding design point changes in θ as is indicated by the computed confidence intervals. As second mathematical example we consider an exponential function with additional power term, given by f2:( x, ( θ 1 ,θ 2)) → θ 1·exp(− θ 1· x )+ x θ 2 , (24) with X= [0.8,5] ⊂R and θ∈R2 . For illustration purposes, we consider the special case where only the second parameter θ2 is of considerable uncertainty. This means, the parameters are samples for a normal distribution with mean θ0= (1.05,−1.28)T and covariance matrix Σ=(0 . 00 . 0 0.00.01 ). (25) This allows for examining the influence of the uncertainty of the parameter θ2 on the optimal design. Note, that the variance Var [θ1]=0 , whereby the parameter θ1 is deterministic and not uncertain. The presented approaches can also be applied in this setting, where some of the parameters θ are fixed at their mean θ0 . Using the mean parameter value θ0 as reference value, we obtain an optimal design with 3 optimal support points. To compute these values we replace the design space X with a fine grid of 841 candidate points. The grid is given by Representative Weight Variance 1 0.7486 0.455 0.052 2 1.0 0.5 0.0 Table 2 Computed clusters for the exponential function f1 with 100 sampled parameter values 1 3 Page 15 of 28 109 M. Bubel et al. Xgrid ={0.8+0.005 ·j|j=0,...,840}. (26) The values of the optimal weights and the optimal support points are given in Table 3. We see, that the optimal design for the reference value ˜ θ= (1.05,−1.28)T consists of three support points. This can however change, even for small changes in the value of the reference parameter. For example, the optimal designs for the reference parameters (1.05,−1.26)T and (1.05,−1.4)T consist of only 2 optimal design points. Corresponding results can be seen in Fig. 4. We sample M= 100 points θj from the distribution of θ and compute the optimal designs for the corresponding reference values. Again we use the grid Xgrid to compute the corresponding designs ξj=ξ∗(θj) . We average the designs to obtain ξ =1 M· M ∑ j =1 ξj . (27) The averaged weights are plotted in Fig. 5a. We observe in this figure that we have three regions in the design space which are of interest, the borders x=0.8 and x=5.0 , as well as a small region at approximately x≈1.85 . The location of the three optimal design points only varies slightly when the value of θ changes. We apply the clustering approach to the averaged weights, with a radius of r=0.05 and a minimum weight of ωmin =0.1 . Three clusters are identified, matching the observations from the average design. The cluster representatives, weights and confidence intervals are plotted in Fig. 5a and are also given in Table 4. Note, that the plotted cluster representatives overlap the average weights on the borders of the design space. The clusters have a total weight of 1.0 and thus cover all design points with positive average weight. Next we compute the derivatives Dθx∗ i and Dθω∗ i and then approximate the corresponding confidence intervals locally. For the computation of the confidence interFig. 4 Optimal designs for the exponential function f2 and reference parameters ˜ θ= (1.05,−1.28)T ( circles), ˜ θ= (1.05,−1.26)T ( lower triangles) and ˜ θ= (1.05,−1.4)T ( upper triangles) Design point x∗ i Weight ω∗ i 1 0.8 0.1298 2 1.855 0.4878 3 5.0 0.3824 Table 3 Optimal design for the exponential function f2 1 3 109 Page 16 of 28 Optimal experimental design: from design point to design region vals we only consider the derivatives with respect to the second parameter θ2 and the corresponding entries of the covariance matrix Σ . This is done, as θ1 has variance 0 and is uncorrelated to the parameter θ2 . The confidence intervals for the design points are given in Fig. 5b, the confidence intervals for the weights are given in Table 5. We observe, that both approaches identify three regions in the design space which are of interest. This matches the observation from the average design. The two points on the border of the design space, x∗ 1=0.8 and x∗ 3=5.0 , are correctly identified by both methods. Both methods assign these points confidence intervals of length 0. However, the local approximation assigns the optimal point x∗ 2=1.855 a large confidence interval. The average design and the clustering approach however indicate, that this point only changes slightly. This deviation is probably due to the local behavior of the implicit function approach, whereas the clustering approach has access to global data. From the localized approximation we also obtain confidence intervals for the optimal weights. The weights ω∗ 1 and ω∗ 3 have large confidence intervals. This indicates, that these weights vary strongly in the uncertain parameter θ . The weight ω∗ 2 on the other hand has a smaller confidence interval and we thus conclude it is more stable. Design point x∗ i Weight ω∗ Confidence interval for ω∗ 1 0.8 0.2245 [−6 . 0242 , 6 . 2837] 2 1.855 0.4978 [0.1213, 0.8543] 3 5.0 0.2777 [−5.4051,6.1698] Table 5 Confidence intervals of the optimal weights ω∗ for the mathematical exponentialpower model Representative Weight Variance 1 0.8 0.2245 0.0 2 1.8594 0.4978 0.014 3 5.0 0.2777 0.0 Table 4 Computed clusters for the exponential function f2 and 100 sampled parameter values Fig. 5 Resulting optimal designs from the Monte Carlo sampling, the clustering approach and the approximated confidence regions for the exponential function f2 1 3 Page 17 of 28 109 M. Bubel et al. We have seen, that the number of support points changes in the parameter value θ . The disappearance and in particular the appearance of such points is not a local behavior, and thus we cannot expect the local method to correctly predict this behavior or approximate it. However, the large confidence regions assigned to the weights ω∗ 1 and ω∗ 3 , indicate that these weights may go to 0. Using the approximated gradients we can even estimate for which parameter values this behavior occurs. 4.2 Application example – the flash As application example we consider a flash, which for example can be found in Vanaret et al. (2020); Seufert et al. (2021, 2024). We give a brief description of the model and the OED setup, where we closely follow Vanaret et al. (2020). The flash is a pressurized container, into which a liquid input stream ˙ F with composition vector xF enters. A heat duty ˙ Q is applied in the container so that vapor and liquid phase are in equilibrium. We obtain two output streams, a liquid output stream ˙ L at concentration yL and a vapor output stream ˙ V at concentration yV . The vapor-liquid equilibrium is described by the so called MESH equations. For an input stream consisting of two components (methanol and water) these equations are given as Biegler et al. (1997) ●Mass balances ˙ F ·x F m =˙ V·y V m +˙ L·y L m ˙ F·x F w= ˙ V·y V w+ ˙ L·y L w (28) ●Equilibrium P ·y V m = P 0 m ( T ) ·y L m·γm ( y L m,y L w,T,θ ) P·y V w=P 0 w(T)·y L w·γw(y L m,y L w,T,θ) (29) ●Summation xF m+xF w=yV m+yV w=yL m+yL w=1 (30) ●Heat balance ˙ Q+˙ F·HL ( x F m,x F w,T F )= V·HV(y V m,y V w,T)+ ˙ L·HL(y L m,y L w,T) (31) The temperature TF denotes the temperature of the input feed, the value T denotes the temperature at equilibrium in the flash. The functions P0 m,P0 w denote the vapor pressure of the pure components and the functions HL and HV denote the enthalpies of the molar liquid and vapor streams. The model parameters θ=(a12,a 21,b 12,b 21) denote parameters of the activity coefficient model γ Renon and Prausnitz (1968) used to describe the thermodynamic equilibrium. We can rewrite the MESH as given above as an implicit model 1 3 109 Page 18 of 28 Optimal experimental design: from design point to design region 0= g ( x, s, θ ), (32) where s:== (yV,yL,T,x F w) is a set of state variables, the input variable x:== (xF m,P) consists of the concentration of methanol in the feed and the pressure at which the flash operates, and g corresponds to the four Mesh equations. We define the mapping h:R6→R2,s→ (T,yV m) and the model function f( x, θ ) :== h ( s ( x, θ )), (33) where s(x, θ) denotes an solution of the implicit model (32) given the inputs x=(xF m,P) and the model parameters θ . An exhaustive discussion on the existence of derivatives for this system can be found in Schmid 2024 Schmid et al. (2024). The flash model example demonstrates, that the proposed methods can be straightforwardly extended to implicitly defined models, which are commonly encountered in practical case studies. Optimal experimental design for implicit model functions has also been considered by Duarte Duarte et al. (2021). As resulting design space we select X= [0.0,1.0] ×[0.5,5.0] , i.e. the molar fraction xF m is chosen from within the range (0.0 to 1 . 0) mol mol and pressure from the range (0.5 to 5.0) bar. As outputs of the model we measure the values ( T,y V m)= f ( x, θ ) , the temperature T at equilibrium and the molar fraction yV m of methanol in the vapor output stream, y. We assume the model parameters θ=( a 12 ,a 21 ,b 12 b 21)T to be normally distributed with mean E[θ]=(−3.9402,6.3037,1337.558,−1891.945)T. (34) The values for the standard deviations are obtained by assigning each parameter a 5 percent deviation from its nominal value, i.e. Cov [θ]ii = (0.05 ·θ0 i)2 , resulting in the covariance matrix Cov [θ]≈    0 . 0388 0 0 0 00.0993 0 0 0 0 4472.65 0 0 0 0 8948.64   . (35) Using the expectation θ0 as reference value we can compute an optimal design. For this purpose we replace the design space X with a uniform grid, with 101 grid points in the xF dimension and 91 grid points in the P dimension. The resulting grid is given by X grid =  i 100,j 20  T  i=0,...,100, j= 10,...,100  (36) and the optimal design ξ∗ is given in Table 6. For the flash we sample M= 200 different values θj of the model parameters from the described distribution N(θ0,Σθ) and compute the corresponding optimal 1 3 Page 19 of 28 109 M. Bubel et al. designs ξj on the specified grid Xgrid . In Fig. 6a the corresponding support points are plotted. The size of the points and the coloring hereby indicate the weight assigned to each point. As all computed optimal design points have a molar fraction xF m≤0.6 , we only plot the relevant part of the design space and exclude all values xF m>0.6 . Next, we apply the clustering algorithm to the averaged weights. The result of this method can be seen in Fig. 6b. Here, each cluster cj is represented via its mean point cj,mean its weight ω(cj) and its confidence ellipsoid. For the clustering we have used a radius r=0.04 and a minimal weight ωmin =0.04 . The results are also given in Table 7. The clusters have a total weight of 0.9015. Thus design points with positive average weights corresponding in sum to 0.0985 are not covered by the clusters. Last we compute the derivatives Dθx∗(θ0),D θω∗(θ0) and approximate the confidence regions via the implicit function theorem. The confidence intervals for the Representative Weight Cov11 Cov12 Cov22 1 (0 . 1533 , 0 . 5)T 0.2030 0.0023 0.0 0.0 2 (0 . 3310 , 1 . 53)T 0.1775 0.0005 − 0.0022 0.0475 3 (0 . 2224 , 1 . 65)T 0.0495 0.0002 − 0.0005 0.0056 4 (0 . 0746 , 5 . 0)T 0.2330 0.0001 0.0 0.0 5 (0 . 3464 , 5 . 0)T 0.2385 0.0014 0.0 0.0 Table 7 Computed clusters for the flash Fig. 6 Resulting optimal designs from the Monte Carlo sampling, the averaging of weights and the clustering approach for the flash Design point x∗ i Weight ω∗ i 1 (0.07,5.0)T 0.2317 2 (0.1,2.0)T 0.0493 3 (0.16,0.5)T 0.2497 4 (0.33,1.7)T 0.2257 5 (0.34,5.0)T 0.2435 Table 6 Optimal design for the flash on the grid Xgrid and reference parameters θ0 1 3 109 Page 20 of 28 Optimal experimental design: from design point to design region weights are given in Table 8. Unfortunately the approximated confidence regions are very large and cannot be plotted appropriately, thus we forego such a plot. From the results we observe that the clustering approach gives reasonable approximations of the confidence region. We validate these results by computing the values Y j :== ∑ i ∈ Ij ωi ( θj ), where θj,j =1,...,200, denote the sampled parameter values from the clustering approach and Ij denotes the indices of the optimal design points x∗ i(θj) , which are within one of the approximated confidence regions. For the values Yj some statistical figures are given in Table 9. We observe, that for half the samples we achieve a coverage of 94% , and we only have a coverage of less than 84% in one quarter of evaluated cases. In contrast to the mathematical examples, the local approach does not yield reasonable approximations, resulting in confidence ellipsoids that are excessively large and exceed the physically plausible range for the concentration values x1 . We attribute this to two main reasons. First, the variance considered in the flash example is substantial, leading to approximation errors in local approach. Second, the nonlinearity of the optimal experimental design in θ contributes to further approximation errors when employing linearization-based methods. Both of these effects were discussed in Sect. 3.3, where we noted that the implicit function approach is only locally accurate. For this specific example, we conclude that the primary reason for the failure of the local approach is the high nonlinearity present in the NRTL model and its parameterization. As analyzed by Labarta et al. Labarta et al. (2022), the NRTL model can exhibit highly nonlinear behavior, which appears to be the case with the chosen parameters in this example. While parameter transformations, as proposed by Bates and Watts Bates and Watts (1981), may mitigate nonlinearity, they are not universally applicable to all nonlinear regression models, including the NRTL model; Statistical figures mean 0.914175 standard deviation 0.095014 minimum 0.499907 25% quantile 0.848292 50% quantile 0.942115 75% quantile 1.000000 maximum 1.000000 Table 9 Mean, Standard Deviation and quantiles of the values Yj Design point x∗ i Weight ω∗ Confidence interval for ω∗ 1 (0.07,5.0)T 0.2317 [0.1807,0.2828] 2 (0.1,2.0)T 0.0493 [−0.1309,0.2296] 3 (0.16,0.5)T 0.2497 [0.2476,0.2518] 4 (0.33,1.7)T 0.2257 [0.1272,0.3243] 5 (0.34,5.0)T 0.2435 [0.2103,0.2767] Table 8 Confidence ellipsoids of the optimal weights ω∗ for the flash 1 3 Page 21 of 28 109 M. Bubel et al. hence, we cannot circumvent the nonlinearity in this instance. To derive reasonable information from the local approach, we recommend decreasing the covariance of the model parameters. For this particular example, we find that multiplying the diagonal entries in (35) by a factor of 0.1 to 0.01 is reasonable to achieve confidence ellipsoids comparable in size to those obtained from the clustering approach. Subsequently, the computed local derivatives (see Appendix A) of the optimal design points from the flash example can be utilized to determine how much and in which direction the optimal design points shift under uncertainty. Although some interference is necessary (scaling of parameter covariance), the experimenter may still gain insights into which experiments are more and less stable, respectively. Ultimately, the chosen noise level both reflects a realistic scenario from the process engineering domain and nicely demonstrates the superior stability of the clustering approach over the local approach, which is why we decided to keep the relatively large parameter covariances in (35) for the flash example. 5 Conclusion In this paper we have presented two approaches to handle uncertainty and unknown model parameters in optimal experimental design. The first approach is a clustering approach which utilizes the DBSCAN algorithm. This approach first averages the optimal weights and then utilizes the algorithm to identify regions in the design space with a high weight density. We identify each region with a representative point, weight and covariance matrix. The clusters correspond to design regions in which OED assigns a large weight over all sampled scenarios. Performing experiments in these regions thus seems reasonable, regardless of the true value of the uncertain parameters. Via the covariance matrix we also obtain information on how uncertain each cluster is and can approximate confidence regions. The main disadvantage of this approach is, that we need to sample the parameter space Θ and compute optimal designs for each of the samples. Similar to alternative approaches like the worst-case design or the average case designs this can be computationally expensive, in particular in high dimensional settings. The second approach approximates the confidence ellipsoids of the optimal design points x∗ i and optimal weights ω∗ i . This approach utilizes the implicit function theorem to compute the derivatives Dθx∗ i and Dθω∗ i . As this is a local approach in x∗ i,ω∗ i and in the nominal parameter θ0 , we require significantly less evaluations of the model f and its Jacobian matrix Dθf(x, θ) . In the examples we observe, that this approach provides additional information to each computed optimal design ξ∗ . Via the confidence ellipsoids an experimenter can make a more informed decision on which experiments to perform. However, this approach only is a local approach and may lead to wrong and misleading approximations. This is particularly the case, when the optimal design points and optimal weights are very non-linear in θ0 . Large variances in the uncertain parameters can also be problematic, as the linearizations applied get worse for 1 3 109 Page 22 of 28 Optimal experimental design: from design point to design region large deviations from the nominal (expected) scenario. This effect is seen in the flash example of this contribution. Combining both approaches, the data points from the clusters can be used to compute higher moments at each cluster. These moments can be used to investigate the non-linearity of each optimal design point. This can then indicate, how much we trust the local approximations at each point. Last, we have seen that in some examples the weight ω∗ i assigned to an optimal design point can go to 0 as θ changes. The optimal design point x∗ i is then no longer part of the design and is removed. Via the gradient and confidence intervals on the optimal weights ω∗ i , we can predict under which parameter changes the weight converges to 0. The opposite behavior, when a design point is suddenly assigned some positive weight, is more difficult to predict and can be part of future research. Derivatives of designs via the implicit function theorem In this Appendix we give details on the application of the implicit function theorem to the Karush-Kuhn-Tucker (KKT) conditions of the optimal experimental design problem. This follows ideas from the so-called parametric optimization problem Jarre and Stoer (2004). We consider the optimization problem min ω,x 1,...,xn ˜ Ψ( ω,x1,...,x n,θ 0) s.t. gi(ω)≤0for i=0,...,n h j (x i ) ≤ 0 for i=1,...,n, j =1,...,2 · d x, (37) introduced for the implicit function approach. To improve notation we denote x=(x1,...,x n)∈Rn×d x. The Lagrangian function for this problem is given by Lθ:( ω, x ,λ )→ L θ( ω, x ,λ ), (38) with L θ ( ω, x ,λ )= ˜ Ψ(ω,x,θ)+ n ∑ i=0 λi·gi(ω)+ n ∑ i=1 2·dx ∑ j=1 λij ·hj(xi) . (39) As the points x∗ i(θ0) and weights ω∗ i(θ0) are optimal for the value θ0 , we find Lagrange multipliers λ∗ such that the KKT conditions for Lθ 0 are fulfilled. Then D ωiLθ0 ( ω ∗ , x ∗,λ ∗)=0, D x i Lθ0 ( ω∗,x ∗ ,λ ∗ )=0. (40) 1 3 Page 23 of 28 109 M. Bubel et al. In the following we denote by I ∗=  (i, j)  hj(xi)=0,i =1 ,...,n, j=1,...,2 · dx  (41) the indices of the active constraints. The corresponding Lagrange multipliers are denoted by λI∗=( λ ij , ( i, j )∈ I ∗), (42) and the constraints by hI∗( x 1 ,...,x n)=( h j( x i) , ( i, j )∈ I ∗). (43) For the implicit function theorem we define the function F :(ω,x, λ, θ) → ( DωLθ ( ω, x ,λ ) DxLθ0(ω,x,λ) ) (44) We consider the derivative of F with respect to (ω∗,x∗,λ ∗ 0,λ ∗ I∗) , given by the matrix DF :== ( DF∗ ( Dh∗ ) T Dh∗0 ), (45) where we denote DF∗=( DωF ( ω ∗ , x∗ ,λ ∗ ,θ 0) ,D x F ( ω ∗ , x∗ ,λ ∗ ,θ 0)), (46) and Dh ∗= ( Dωg0 ( ω ∗) 0 0DxhI∗(x∗) ). (47) We then consider the derivative of F with respect to θ , given by DG :== ( DθF ( ω ∗ , x∗ ,λ ∗ ,θ 0) 0 0 ) (48) If the matrix DF is non-singular, we can apply the implicit function theorem and compute the derivatives    Dθω ∗( θ 0) Dθx∗(θ0) Dθλ∗ 0(θ0) D θ λ ∗ I∗( θ0 )    =−DF−1· DG (49) 1 3 109 Page 24 of 28