Penalization method to convert Bayesian optimization methods into batch multi-objective Bayesian optimization methods

Figures
Abstract
Bayesian optimization is a surrogate-based global optimization method that is increasingly being used for engineering design. However, most methods are designed to be sequential, and in optimization problems where the functional evaluation can be easily parallelized, batch methods are more effective at reducing real-time computation. While existing batch multi-objective Bayesian optimization (B-MOBO) methods achieve strong performance, they are typically purpose-built from scratch, limiting ability to leverage the extensive library of proven sequential acquisition functions for parallel settings. In this paper, a generalizable transformation methodology is developed that converts compatible single- or multi-objective sequential Bayesian optimization methods into B-MOBO methods. The key innovation is a Euclidean distance-based composite acquisition function with multi-objective penalization averaging, which combines multiple sequential acquisition functions while incorporating objective-wise penalization information from all objectives and preventing redundant sampling in batch selection. To demonstrate the schema, this methodology is applied to two representative sequential acquisition strategies, the expected improvement method (single-objective) and the quality metrics method (multi-objective), to develop two new B-MOBO variants. These new methods are then compared against their sequential counterparts and established B-MOBO methods using both computation time and real-time analysis. Results show that the new B-MOBO methods are as effective as existing methods at solving optimization problems in terms of solution quality and show improvements when real-time is considered.
Citation: Holder A, DeBruin H, Sestito JM (2026) Penalization method to convert Bayesian optimization methods into batch multi-objective Bayesian optimization methods. PLoS One 21(7): e0354346. https://doi.org/10.1371/journal.pone.0354346
Editor: Yuanchao Liu, Northeastern University, CHINA
Received: November 7, 2025; Accepted: July 6, 2026; Published: July 21, 2026
Copyright: © 2026 Holder et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Data Availability: All programming files are available from GitHub (https://github.com/sestito/BatchMOBO_Public).
Funding: The author(s) received no specific funding for this work.
Competing interests: The authors have declared that no competing interests exist.
Introduction
To solve engineering design optimization problems efficiently, surrogate-based global optimization methods have been developed. Different from population based global optimization, such as particle swarms [1] and genetic algorithms [2], in surrogate-based optimization, a surrogate of the objective is constructed and updated during the searching process. The surrogate is used to identify the most promising region and guide the search. Since the surrogate is much cheaper to evaluate than the original objective, the overall sampling cost is reduced.
One such surrogate-based method is Bayesian optimization (BO). The typical surrogate used in BO is the Gaussian process (GP) model, which also quantifies the uncertainty. In BO, an acquisition function is constructed from the surrogate to identify the candidate sample most likely to improve the objective value, guiding the sequential sampling process. This surrogate-based approach is ideal for optimizing black-box objectives where evaluations are time-consuming and costly [3–8].
Solving single-objective problems with sequential Bayesian optimization (S-BO) has been well explored in design optimization [9–12]. However, in real-world applications, design problems usually involve optimizing multiple competing objectives, such that improving one objective often degrades another. This requires trade-offs between different objectives meaning there is no one optimal value. Instead, there is a set of optimal solutions, or values for the input parameters, that are optimal and form what is called the Pareto front, separating the feasible region from the infeasible one. Multi-objective Bayesian optimization (MOBO) methods have been developed to obtain the Pareto front. In sequential sampling MOBO (S-MOBO), a set of initial samples and their corresponding solutions are passed into the S-MOBO algorithm which recommends a sample that is most likely to exist on the Pareto front. S-MOBO algorithms can be used to discover the Pareto front of the objective space quickly and efficiently [13,14].
One major limitation of S-MOBO methods is their inherently sequential nature. Only one iteration of the design can be built at a time, which increases the real-time required for optimization. In many modern engineering applications, computational resources support parallel evaluation of objective functions through high-performance computing clusters, multi-core processors, or distributed experimental facilities [15], and recent surrogate-assisted methods have been developed to exploit these resources efficiently across a range of optimization problem classes [16,17]. This capability motivated the development of batch multi-objective Bayesian optimization (B-MOBO), which can recommend multiple sets (a batch) of design parameters simultaneously, whose solutions are expected to exist along the Pareto front [3,18,19]. By evaluating multiple candidates in parallel, B-MOBO methods can dramatically reduce real-time computation costs even though the total number of function evaluations may remain similar or slightly increase compared to sequential methods. This makes parallelization essential for time-critical engineering design problems where waiting for sequential evaluations is impractical.
The development of batch methods has progressed through several stages. Early batch approaches for single-objective optimization focused on strategies such as constant liar, Kriging believer [20], and q-Expected Improvement [21]. These methods addressed the fundamental challenge of batch selection on how to choose multiple points when the outcomes of other points in the batch are unknown. Gonzalez et al. introduced local penalization as an efficient heuristic for batch generation, using Lipschitz continuity estimates to create local repulsion around selected batch members, thereby preventing redundant sampling [15]. This approach achieved large computational savings compared to more elaborate alternatives. Subsequently, Alvi et al. extended local penalization to asynchronous settings, demonstrating that asynchronous methods can outperform synchronous batch approaches in both real-time and sample efficiency [22]. These strategies have since been adapted for multi-objective settings.
While existing B-MOBO methods have been developed, most approaches are either purpose-built B-MOBO methods or individual adaptions of S-MOBO [23–26]. Recent B-MOBO methods employ various strategies including ensemble acquisition functions, hypervolume-based criteria, and diversity-guided selection, as well as advanced surrogate modeling techniques for high-dimensional problems [27,28]. For instance, Daulton et al. developed noisy Expected Hypervolume Improvement method which handles noisy objectives and scales batch selection from exponential to polynomial complexity through cached box decomposition [29]. Lin et al. proposed parallel multi-fidelity MOBO approaches incorporating strategies for variable-fidelity optimization with correlated objectives [30]. Luković et al. introduced diversity-guided B-MOBO [3]. Recently, Ngo et al. presented MOBO-OSD, which generates diverse Pareto solutions via orthogonal search directions combined with Pareto front estimation for batch optimization [31].
A wide array of S-BO and S-MOBO methods already exist and are highly effective at discovering the Pareto front [7,32–34]. S-BO methods are effective for varying optimization problems, where selecting the best acquisition function for the problem at hand can reduce the time to discover the true Pareto front. However, while existing B-MOBO methods employ strategies such as local penalization and ensemble acquisition functions, a fundamental gap remains. These methods are typically designed from scratch for batch multi-objective optimization rather than providing a systematic framework to convert existing sequential methods. If compatible S-BO and S-MOBO methods can be turned into B-MOBO methods, these extensive libraries of proven sequential acquisition functions can be leveraged and parallelized for engineering design problems, combining their effectiveness with the advantage of batch selection.
This work addresses this gap by developing a Euclidean distance-based composite acquisition function with multi-objective penalization averaging that transforms compatible single-objective and multi-objective sequential acquisition functions into batch multi-objective methods. Unlike existing local penalization methods that primarily focus on preventing redundant sampling through Lipschitz-based exclusion zones, our penalization averaging approach incorporates objective-wise penalization information from all objectives while maintaining computational efficiency. The resulting schema enables systematic conversion of compatible S-BO and S-MOBO methods into B-MOBO methods without problem-specific modifications. To demonstrate the generality and effectiveness of this schema, we apply it to two well-established methods, the expected improvement method [34] and quality metrics method [35], converting them both to B-MOBO methods. These new methods are then compared to their S-MOBO counterparts as well as against an established B-MOBO method to evaluate the performance gains from parallelization and the competitiveness of the proposed transformation approach.
Methods
Proposed batch multi-objective Bayesian optimization transformation algorithm
The goal of the proposed B-MOBO transformation algorithm is to convert compatible S-BO and S-MOBO methods into B-MOBO methods. This process first involves the conversion of a S-BO method into a S-MOBO method, followed by the conversion of a S-MOBO method into a B-MOBO method. The transformation algorithm is outlined in Fig 1. A compatible method is defined as one whose acquisition function is non-negative, is computable pointwise from the trained surrogate, and is formulated such that larger acquisition values indicate more promising samples. Acquisition functions that do not satisfy these conditions require modification, such as shifting, before the transformations can be applied.
The acquisition function is assumed to be the combined acquisition function for S-MOBO. (a) illustrates the overall algorithm while (b) is the Batch Recommendation Loop.
To understand how the MOBO formation works, it is important to understand the BO formation. The goal of BO is to select a sample ( for n input parameters) that minimizes or maximizes, depending on the problem, an objective (f(x)). To do this, a surrogate that is cheap to evaluate is fit to the existing data to allow for predictions for the objective based on a set of input parameters. This surrogate is then used to create an acquisition function (a1(x)), that is maximized to yield the sample most likely to improve the objective. In this work, all objective functions are minimized while all acquisition functions are maximized, such that larger acquisition values indicate more promising samples.
To convert single-objective optimization methods into multi-objective optimization methods, a new acquisition function must be defined that can be used to select the next sample. First, an acquisition function is developed for each of m objectives. Next, the new combined acquisition function (a(x)) is formulated as the Euclidean distance defined as
(1)where is the acquisition function built for the ith objective. The new acquisition function is sufficient to convert compatible S-BO methods into corresponding S-MOBO methods. This formulation can be used for any number of objectives, even one. Hereafter, the following calculations are derived assuming the combined acquisition function.
For the conversion of S-MOBO methods into B-MOBO methods, a penalization method is developed based on a penalization method by González et al. which was originally written for single-objective maximization [15]. The González penalty method is adopted and modified to work for the minimization of objective functions as well as multiple objectives. The penalization method for single-objective optimization is defined such that the first new sample x1 is selected as normal. Each subsequent sample is selected from a penalized acquisition function defined as
(2)where is the ith new recommended sample. The net penalty function () is then the average of each objective’s penalty function such that
(3)where is the penalty function for each objective defined as
(4)(5)where and are the mean and standard deviation for predicted by the surrogate trained on the known samples and solutions. is the minimum observed value of the kth objective and serves as the reference optimum in the penalty. is a Lipschitz constant calculated as the maximum gradient magnitude of the kth objective’s trained surrogate within the input parameter bounds. Together, and set the radius of the exclusion zone applied around each previously selected batch member.
To adopt this penalization method for multi-objective optimization, the penalization function can be computed for each objective, such that there are m penalization functions. The penalization functions are then combined using a simple average, and then multiplied with the acquisition function. The average is selected over alternatives such as product, minimum, or weighted rules for three reasons. First, the average is invariant to the number of objectives, whereas the product compounds the penalties such that the net penalization strengthens as objectives are added, and the minimum retains only the single most-penalizing objective. As such, a region strongly penalized by one objective but promising for the others is over-excluded under the product and minimum rules, while the average preserves the penalization contribution of every objective and limits the influence of any single one. Second, a weighted average would require objective-specific weights to be specified or tuned for each problem, reintroducing the problem-specific adjustment that the transformation is designed to avoid. Third, for a single objective the average reduces exactly to the penalization of González et al. [15], preserving the original single-objective behavior. The simple average is also the most straightforward rule to apply to any acquisition function. L becomes and M becomes . As such, z becomes , becomes and the penalized acquisition function becomes
(6)In the above, i indexes the previously selected batch members that penalize the candidate, j indexes the batch candidate currently being selected, and k indexes the m objectives.
A visual representation of S-BO to B-MOBO transformation algorithm is demonstrated in Fig 2. This new penalization acquisition function can now be applied to existing S-BO and S-MOBO methods to create B-MOBO methods.
a), c), and e) outline how a S-BO method is converted into a B-MOBO method through the formation of a new acquisition function. b), d), and f) are the corresponding penalty functions for the a), c), and e) acquisition functions. g) is the product of e) and f) to demonstrate the new acquisition function that the next batch sample will be selected from.
Existing Bayesian optimization methods to transform
This methodology is applied to the expected improvement (EI) method [34], which is a S-BO method, and the quality metrics method [35], which is a S-MOBO method, to convert them both into B-MOBO methods.
For the quality metrics method, the relative improvement (I) for the quality metrics (QM) was changed from the relative change to the relative difference method defined as
(7)where are the current n sample points and x is the added sample set. It was changed to the relative difference method to account for quality metrics that can be negative. The quality metrics used for the optimization are the dimensional overall spread (DOS) and integral hyperarea difference (IHD). The DOS is defined by the minimum bounding hyperrectangle that encompasses the current estimated Pareto front. This is also known as a hypervolume with a reference point of the objective-wise minimum of the Pareto points. The IHD is defined as the area under a set of simplexes that estimates the Pareto front [35].
Results and discussions
To demonstrate the effectiveness of the proposed algorithm, the algorithm is applied to the multi-objective QM optimization method and single-objective EI method to create two new B-MOBO methods. To express the EI method as a multi-objective method in the Euclidean domain, it is henceforth denoted as EIMe. To quantify these new B-MOBO methods, the numerical examples ZDT1, ZDT2, ZDT3, FON [36,37], and DTLZ2 [38] and a welded beam engineering design problem [39] are adopted for optimization. To evaluate the effectiveness of the new methods, several comparisons will be made. First, each new B-MOBO method will be compared to their S-MOBO counterparts, henceforth known as a batch size of 1. Second, the new B-MOBO methods will be compared to the B-MOBO Ensemble method by Lyu et al. [14] and the q-Expected Hypervolume Improvement (qEHVI) method by Daulton et al. [29]. All Bayesian optimization methods use a Gaussian surrogate with a Matérn 5/2 kernel to balance smoothness and modeling flexibility. Kernel hyperparameters are optimized via maximum likelihood estimation with ten random restarts. These settings follow standard practices in the Bayesian optimization literature [40,41] and are held fixed across all experiments to ensure fair comparisons. The qEHVI method is implemented using the BoTorch library. Its acquisition function is optimized using multi-start L-BFGS-B with 20 restarts and 1024 raw samples, and the reference point is set to the objective-wise minimum of the observed data offset by 10% of each objective’s range. All methods are executed on comparable computational resources using identical initial sample sets and evaluation budgets.
A performance test is developed to systematically quantify each method for each optimization example and is executed for each optimization method, at each batch size (1, 2, and 3), and for each numerical example. To account for statistical randomness, the performance test is the optimization routine executed 30 times, and the optimization routine is defined as follows: An initial set of 30 samples are created using the Latin Hypercube Sampling [42,43] method and the samples’ corresponding solutions are computed from the functional evaluation. Next, the MOBO method is used to select the next sets of samples corresponding to its batch size: a batch size of one results in one set, a batch size of two results in two sets, and a batch size of three results in three sets. The corresponding solution is calculated for each sample set resulting in one, two, or three new samples and solutions. This data is combined with the original data, the optimization method is used again to recommend the next sample or samples, and the process repeats until a total of 330 samples and solutions are created.
To quantify the performance test, distinct performance metrics are calculated for each of the 30 runs. These metrics are the number of points on the Pareto front (nP), the Pareto quality metric dimensional overall spread (DOS), and the Pareto quality metric integral hyperarea difference (IHD). The mean, standard deviation, and range are calculated for each performance metric. In general, a higher nP is desired and DOS and IHD values closest to the ideal values are desired.
To compare batch and sequential methods, the first type of comparison that will be used is a computation time comparison. To keep computation time constant, function evaluations must be consistent requiring the exact same number of samples to be compared. To illustrate this comparison, let us look at a sequential method versus a batch method of size two. For the sequential method, 30 initial samples are generated, then the optimization is executed 300 times to create a total of 330 samples. For the batch method of size two, 30 initial samples are generated, then the optimization is executed 150 times to create a total of 330 samples. In this instance, both algorithms created 330 samples with 330 functional evaluations and took an equivalent amount of computation time. However, the sequential method required the 300 functional evaluations to be executed in series while the batch method allowed the 300 functional evaluations to be executed in parallel for sets of two. This type of analysis is good for examining the computational time benefits between methods, but falls short when looking at the real-time benefits which batch methods thrive in. The computation time analysis also isolates the two components of the transformation. At a batch size of 1, the penalization product in Eq. 2 is empty, so the results reflect the composite acquisition function alone. Comparing a batch size of 1 against batch sizes of 2 and 3 under the computation time analysis, where the total evaluation budget is matched, therefore isolates the effect of the multi-objective averaged penalization.
To examine the real-time benefits between batch and sequential methods, the number of optimizations is held constant to keep the real-time constant. To illustrate this comparison, a sequential method versus a batch method of size two is used. For the sequential method, 30 initial samples are generated, then the optimization is executed 100 times to create a total of 130 samples. For the batch method of size two, 30 initial samples are generated, then the optimization is executed 100 times to create a total of 230 samples. In this instance, assuming parallelization of functional evaluations, both methods will take an equivalent amount of real-time. It is important to note that the real-time analysis is a parallel-resource budget comparison rather than a direct algorithmic fairness comparison. Under this framing, batch methods with a batch size greater than one access more total functional evaluations within the same wall-clock window, which is the intended advantage of batch methods in applications where evaluations parallelize. The purpose of the real-time analysis is therefore to compare methods at matched batch sizes, that is, the new B-MOBO methods against existing B-MOBO methods, and to show how each method exploits available parallel compute. Differences observed across batch sizes within a single method should not be interpreted as intrinsic methodological superiority. To look at the benefits of the new batch methods, both the computation time and real-time comparisons will be utilized.
Two-dimensional numerical examples
The two-dimensional numerical examples are defined in Table 1. All numerical examples are executed with the dimension of input variables as n = 3.
The ZDT1 problem is first examined to demonstrate the different Bayesian optimization methods. The ZDT1 test problem’s Pareto front can also be analytically solved as [44]
(8)From the analytical solution of the Pareto front, the ideal quality metrics are determined to be DOS = 1 and IHD = 1/3.
To examine the ZDT1 performance test data, the performance metrics for each of the 30 samples are plotted as box-whisker’s plots with outliers removed. An outlier is determined to be outside a ± 2.698 times standard deviation. Fig 3 a – c has the computation time analysis while Fig 3 d – f has the real-time analysis.
a., b., and c. are generated using a computation time analysis while d., e., and f. are generated using a real-time analysis. A larger number of Pareto points indicates better performance, and dimensional overall spread (DOS) and integral hyperarea difference (IHD) values closer to the ideal values for the problem indicate better performance.
For the computation time analysis of the ZDT1 problem (Fig 3 a – c), we would expect the batch size of 1 to outperform the batch size of 2 and 3. This is because each new sample in the batch size of 1 uses all previous data to inform the next selected sample, which is not the case for the other batch sizes. For ZDT1 case, this is not observed as all batch methods appear to perform similarly with their respective optimization method. This could be due to the simplicity of the ZDT1 example. It is also apparent that each optimization method captures the ZDT1 Pareto front with different success. When examining the number of Pareto points, the QM method is the best, followed by Ensemble, followed by EIMe and qEHVI which show similar results. If looking at DOS, the QM method is the best where the qEHVI overestimates, the EIMe overestimates, and the Ensemble method underestimates. Overall, under the tested ZDT1 setting, the proposed batch variants did not show a clear performance disadvantage relative to their sequential counterparts when computation time is considered.
For the real-time analysis of the ZDT1 problem (Fig 3 d – f), the higher batch size produces a better approximated Pareto front than the lower batch size, as expected when parallel evaluation is available. For the number of Pareto points for the QM method, a batch size of 3 has more Pareto points than a batch size of 1. For all three methods, the DOS and IHD values all converge to the ideal values as batch size is increased. When considering real-time under parallel evaluation, the higher batch size makes more effective use of the available parallel evaluation resource within the same wall-clock window, and an optimization method more suited to solving the ZDT1 problem, such as the QM method, can be selected using the new schema.
The ZDT2 problem also has an exact solution for the Pareto front as
(9)This yields an optimal DOS that equals 1 and the optimal IHD of 2/3. The ZDT3 solution [36] is a bit more complex, where and
(10)and
(11)This yields an optimal DOS of approximately 1.51 and an optimal IHD of 0.73. The ideal parameters for FON are DOS = 1 and IHD = 2/3 [37].
For FON, a more complicated two-dimensional Pareto front, similar results as the ZDT1 test problem are observed. When examining the computation time analysis (Fig 4 a – c), it is again observed that the batch methods perform similarly to the sequential method for their respective optimization method. Specifically, for the EIMe, Ensemble, and qEHVI methods, the batch methods are slightly worse. Again, this is expected since the sequential method uses all previous information to guide the sample selection process. However, when going from a computation time analysis to a real-time analysis (Fig 4 d – f), all methods produce better approximated Pareto fronts as larger batches complete more functional evaluations under the parallel-resource budget, with the QM method seeing the largest improvement.
a., b., and c. are generated using a computation time analysis while d., e., and f. are generated using a real-time analysis. A larger number of Pareto points indicates better performance, and dimensional overall spread (DOS) and integral hyperarea difference (IHD) values closer to the ideal values for the problem indicate better performance.
Additionally, when looking at the computation time analysis (Fig 4 a – c), there are large discrepancies between the approximated Pareto fronts calculated using different methods. Even though the EIMe and QM methods have similar DOS and IHD values, the number of Pareto points the QM method has is over double the amount the EIMe method calculates. Additionally, the Ensemble and qEHVI methods perform poorly compared to both the EIMe and QM methods. This indicates that it is important to choose the best BO method for the problem.
The observations on the ZDT2, ZDT3, and FON problems further suggest that when a functional evaluation can be parallelized, the final approximated Pareto front can see improvements over traditional sequential optimization within the same wall-clock window. Additionally, the type of optimization method plays a large role in the established approximated Pareto front, as shown with FON. The ability to select the optimization method most optimal for a specific problem plays a key role in finding strong Pareto fronts. Coupling the correct optimization method with parallelization can result in an improved approximated Pareto front. This is why the ability to turn any compatible sequential optimization method into a B-MOBO method instead of relying on only developed B-MOBO methods can be a powerful tool for optimization problems.
Higher-dimensional numerical examples
Multi-objective optimization is not just for two-dimensional optimization. These methods must be able to scale to higher dimensions. In this work, these methods are tested with up to five objectives. The QM has difficulties with four or more objectives and is partially excluded from this analysis. The qEHVI method has difficulties with three or more objectives and is excluded from this analysis. To accomplish three to five objective problems, the numerical example DTLZ2 is used. Here, DTLZ2 is formulated to minimize m objectives using n input parameters [38]
(12)where and . The ideal Pareto solutions exists when through equal 0.5. The Pareto front is the surface of a unit radius m-sphere centered at the origin existing in the first orthant. This results in an optimal DOS of 1 and IHD of of the hypervolume of the n-sphere of radius 1. The worst possible guesses results in an n-sphere of radius , and the worst possible DOS and IHD are and of the hypervolume of the n-sphere of radius . If the Pareto front has a low diversity, the DOS and IHD calculations can be lower than the optimal values.
The DTLZ2 test problem is executed with m = 3 and n = 4 following the same execution settings as the two objective numerical examples. In three objective space, the optimal DOS and IHD are 1 and respectively. The worst possible DOS and IHD are 3.375 and respectively. The EIMe method performs the best. The QM method and Ensemble method both perform well. QM method calculates the quality metrics better than the Ensemble method and the Ensemble method gets more Pareto points than the QM method as seen in Fig 5. Like the two-dimensional problems, different optimization methods are effective for different problems. When looking at the real-time analysis (Fig 5 d – f), the EIMe method, which performed the best, sees improvements from the additional parallel evaluations. The QM and Ensemble methods see slight improvements. When the number of objectives is increased to 5 (Fig 6), similar trends are observed. The EIMe method performs better than the Ensemble method and both methods see improvements under the parallel-resource budget while showing minimal losses in the computation time analysis.
a., b., and c. are generated using a computation time analysis while d., e., and f. are generated using a real-time analysis. A larger number of Pareto points indicates better performance, and dimensional overall spread (DOS) and integral hyperarea difference (IHD) values closer to the ideal values for the problem indicate better performance.
a., b., and c. are generated using a computation time analysis while d., e., and f. are generated using a real-time analysis. A larger number of Pareto points indicates better performance, and dimensional overall spread (DOS) and integral hyperarea difference (IHD) values closer to the ideal values for the problem indicate better performance.
Engineering example
To compare the methods to an engineering problem, the engineering design of a welded beam that has been employed previously by Ray et al. [39] is used. This is a four input variables (n = 4) and two objective problem. With a load of 6,000 pounds, length of 14 inches, modulus of elasticity of 30 Msi, the objectives to minimize fabrication cost (f1) and end deflection (f2) and constraints shear stress in weld (g1), bending stress in beam (g2), geometric constraint of weld thickness less than beam width (g3), and buckling load constraint (g4) are
(13)where which correspond to weld thickness, weld length, beam height, and beam width all in inches and , , and are the shear stress in psi, bending stress in psi, and buckling load in pounds respectively. The constraints are implemented as a penalty () function as adapted from Coello et al. [45] as
(14)The engineering application yields results consistent with the numerical examples as seen in Fig 7. Real-time analysis reveals comparable performance across batch sizes, with smaller batches showing marginal advantages. The computation time analysis, however, demonstrates an advantage for larger batch sizes. Notably, when applying the qEHVI method to the welded beam problem, the acquisition optimization failed for all batch sizes greater than one. Because the penalty is added equally to both objectives, infeasible designs collapse onto a near-degenerate front that ill-conditions the qEHVI hypervolume box decomposition for batch sizes greater than one. These observations underscore the key advantages of the proposed methodology. By enabling the conversion of a compatible S-BO or S-MOBO method into a B-MOBO method, users gain the flexibility to parallelize sampling when functional evaluations support parallelization, while simultaneously selecting the optimization approach best suited to their specific application.
a., b., and c. are generated using a computation time analysis while d., e., and f. are generated using a real-time analysis. A larger number of Pareto points indicates better performance, and dimensional overall spread (DOS) and integral hyperarea difference (IHD) values closer to the ideal values for the problem indicate better performance.
Batch generation time analysis
To examine the computational efficiency of the transformation algorithm, an experiment is designed to evaluate each method on the ZDT1 benchmark problem with batch sizes ranging from 1 to 3. Additionally, the EIMe method is evaluated with batch sizes from 1 to 10 to assess scalability. For each method and batch size combination, 30 runs were conducted. Each run began with 30 initial samples with an additional 300 samples calculated using the method. Average time per new sample is then calculated. The results in Table 2 reveal distinct computational scaling behaviors across methods. For the EIMe method, the computation time per sample appears to decrease with larger batches, while qEHVI shows an increase as the batch size increases. The other two methods show relatively stable per-sample times. One-way ANOVA confirms that there are significant differences between batch sizes on time per sample for all methods (p < 0.001 for all methods). However, upon further examination, linear regression reveals more distinct patterns. Time per sample decreases significantly with batch size for the EIMe method (p = 0.004), shows no significant relationship for the QM and Ensemble methods (p = 0.494 and 0.178 respectively), and increases significantly for the qEHVI method (p = 0.010).
The decreasing or stable per-sample time observed for the transformation-based methods (EIMe and QM) is consistent with an observed empirical scaling trend on the order of O(b) or better, where b is the number of samples in the batch. The trends reported here are empirical observations from the timing experiment rather than proven complexity bounds. The per-batch-member computational cost does not appear to grow with batch size for the transformation-based methods, in contrast to the increasing per-sample cost observed empirically for qEHVI. When examining batch sizes up to 10 for the EIMe method, the empirical trend persists, indicating that the transformation approach remains efficient at larger batches in the timing experiment. The decreasing per-sample time observed for EIMe is consistent with amortizing certain calculations across batch members, which makes the schema attractive for applications requiring large batch sizes.
Conclusion
In this work, a methodology was developed to convert compatible sequential single-objective and multi-objective Bayesian optimization methods into batch multi-objective Bayesian optimization methods. The methodology was demonstrated on two representative sequential acquisition strategies, the expected improvement method and the quality metrics method, and the resultant batch methods were compared to their sequential counterparts and established B-MOBO baselines. When computation time is considered, the batch methods perform slightly worse than their sequential counterparts under the tested settings but remain competitive. When real-time is considered and functional evaluations can be parallelized, the batch methods access more total functional evaluations within the same wall-clock window and produce better approximated Pareto fronts than their sequential counterparts under this parallel-resource budget. This trend was observed for 2, 3, 4, and 5 objectives under the tested settings.
Depending on the functional evaluations and the problem being solved, different Bayesian optimization methods are recommended. As such, it is important to select the best Bayesian optimization method for the problem. Since this new methodology can convert compatible sequential Bayesian optimization methods to batch multi-objective methods, the most appropriate method for a given problem can be converted into a batch method using the new schema. Batch methods should only be used when the functional evaluations can be easily parallelized. When working with functional evaluations which must be sequential, a non-batch method is better.
Acknowledgments
The authors would like to thank Dr. Nicholas Rosasco for allowing them to use the Valparaiso University Compute Resources to run simulations.
References
- 1. Omkar SN, Khandelwal R, Ananth TVS, Narayana Naik G, Gopalakrishnan S. Quantum behaved Particle Swarm Optimization (QPSO) for multi-objective design optimization of composite structures. Exp Syst Appl. 2009;36(8):11312–22.
- 2. Deb K, Pratap A, Agarwal S, Meyarivan T. A fast and elitist multiobjective genetic algorithm: NSGA-II. IEEE Trans Evol Computat. 2002;6(2):182–97.
- 3.
Lukovic MK, Tian Y, Matusik W. Diversity-guided multi-objective Bayesian optimization with batch evaluations. In: Larochelle H, Ranzato M, Hadsell R, Balcan MF, Lin H, editors. Advances in Neural Information Processing Systems. Curran Associates, Inc.; 2020. p. 17708–20.
- 4. Liu K, Wu T, Detwiler D, Panchal J, Tovar A. Design for Crashworthiness of Categorical Multimaterial Structures Using Cluster Analysis and Bayesian Optimization. J Mech Design. 2019;141(12).
- 5. Pandita P, Bilionis I, Panchal J. Bayesian Optimal Design of Experiments for Inferring the Statistical Expectation of Expensive Black-Box Functions. J Mech Design. 2019;141(10).
- 6. Sarkar S, Mondal S, Joly M, Lynch ME, Bopardikar SD, Acharya R, et al. Multifidelity and Multiscale Bayesian Framework for High-Dimensional Engineering Design and Calibration. J Mech Design. 2019;141(12).
- 7. Choi J, An D, Won J. Bayesian Approach for Structural Reliability Analysis and Optimization Using the Kriging Dimension Reduction Method. J Mech Design. 2010;132(5).
- 8. Chen W, Xiong Y, Tsui K-L, Wang S. A Design-Driven Validation Approach Using Bayesian Prediction Models. J Mech Design. 2008;130(2).
- 9.
Močkus J. On bayesian methods for seeking the extremum. In: Lecture Notes in Computer Science (including subseries Lecture Notes in Artificial Intelligence and Lecture Notes in Bioinformatics). vol. 27 LNCS. Springer Verlag; 1975. p. 400–4.
- 10. Kushner HJ. A New Method of Locating the Maximum Point of an Arbitrary Multipeak Curve in the Presence of Noise. J Basic Eng. 1964;86(1):97–106.
- 11. Yang BS, Yeun Y-S, Ruy W-S. Managing approximation models in multiobjective optimization. Struct Multidiscip Optim. 2002;24(2):141–56.
- 12. Sacks J, Welch WJ, Mitchell TJ, Wynn HP. Design and analysis of computer experiments. Stat Sci. 1989;4(11):409–23.
- 13. Blank J, Deb K. Pymoo: Multi-Objective Optimization in Python. IEEE Access. 2020;8:89497–509.
- 14.
Lyu W, Yang F, Yan C, Zhou D, Zeng X. Multi-objective bayesian optimization for analog/RF circuit synthesis. In: Proceedings of the 55th Annual Design Automation Conference. 2018. p. 1–6.
- 15.
González J, Dai Z, Hennig P, Lawrence N. Batch Bayesian optimization via local penalization. In: Proceedings of the 19th International Conference on Artificial Intelligence and Statistics, AISTATS 2016. 2016. p. 648–57.
- 16. Liu Y, Liu J, Tang L. Expensive Constrained Optimization via Adaptive Switching Between Surrogate-Assisted Constrained and Unconstrained Search. IEEE Trans Evol Computat. 2025;:1–1.
- 17. Li F, Liu Y, Shen H, Pan A, Du W, Jin Y. Robust Multiobjective Evolutionary Algorithm Based on Surrogate-Assisted Robust Distance Metric. IEEE Trans Cybern. 2026;56(5):2902–15.
- 18. Zhang S, Yang F, Yan C, Zhou D, Zeng X. An Efficient Batch-Constrained Bayesian Optimization Approach for Analog Circuit Synthesis via Multiobjective Acquisition Ensemble. IEEE Trans Comput-Aided Des Integr Circuits Syst. 2022;41(1):1–14.
- 19.
Holder A, DeBruin H. Development of an Efficient Batch Multi-objective Bayesian Optimization Method for Engineering Design. 2022. Available from: https://scholar.valpo.edu/cus/1083
- 20.
Ginsbourger D, Riche RL, Carraro L. Kriging Is Well-Suited to Parallelize Optimization. Berlin, Heidelberg: Springer; 2010. p. 131–62.
- 21.
Chevalier C, Ginsbourger D. Fast Computation of the Multi-Points Expected Improvement with Applications in Batch Selection. vol. 7997 LNCS. Berlin, Heidelberg: Springer; 2013. p. 59–69.
- 22.
Alvi AS, Ru B, Calliess J, Roberts SJ, Osborne MA. Asynchronous Batch Bayesian Optimisation with Improved Local Penalisation. 36th International Conference on Machine Learning, ICML 2019. 2019. p. 373–87.
- 23.
Huang J, Yang F, Yan C, Zhou D, Zeng X. A Robust Batch Bayesian Optimization for Analog Circuit Synthesis via Local Penalization. In: Proceedings of the Asia and South Pacific Design Automation Conference, ASP-DAC. Institute of Electrical and Electronics Engineers Inc.; 2021. p. 146–51.
- 24.
Oh C, Bondesan R, Gavves E, Welling M. Batch Bayesian Optimization on Permutations using Acquisition Weighted Kernels. 2021;2. Available from: http://arxiv.org/abs/2102.13382
- 25. Liu J, Zhang D, Wu C, Tang H, Tian C. A multi-finger robot system for adaptive landing gear and aerial manipulation. Robot Auton Syst. 2021;146:103878.
- 26. Wang H, Xu H, Yuan Y, Zhang Z. An adaptive batch Bayesian optimization approach for expensive multi-objective problems. Inform Sci. 2022;611:446–63.
- 27. Li F, Shang Z, Liu Y, Shen H, Jin Y. Inverse distance weighting and radial basis function based surrogate model for high-dimensional expensive multi-objective optimization. Appl Soft Comput. 2024;152:111194.
- 28. Li F, Tian H, Shen H, Zhang X, Liu J, Gong Z. A preference modified inverted generational distance indicator guided algorithm for evolutionary multi-objective optimization. Swarm Evol Comput. 2025;99:102169.
- 29.
Daulton S, Balandat M, Bakshy E. Differentiable Expected Hypervolume Improvement for Parallel Multi-Objective Bayesian Optimization. 2020;10. Available from: http://arxiv.org/abs/2006.05078
- 30. Lin Q, Hu J, Zhou Q. Parallel multi-objective Bayesian optimization approaches based on multi-fidelity surrogate modeling. Aerosp Sci Technol. 2023;143:108725.
- 31.
Ngo L, Ha H, Chan J, Zhang H. MOBO-OSD: Batch Multi-Objective Bayesian Optimization via Orthogonal Search Directions. 2025;10. Available from: http://arxiv.org/abs/2510.20872
- 32. Shu L, Jiang P, Shao X, Wang Y. A New Multi-Objective Bayesian Optimization Formulation With the Acquisition Function for Convergence and Diversity. J Mech Design. 2020;142(9):1–35.
- 33. Couckuyt I, Deschrijver D, Dhaene T. Fast calculation of multiobjective probability of improvement and expected improvement criteria for Pareto optimization. J Glob Optim. 2014;60(3):575–94.
- 34. Zhan D, Cheng Y, Liu J. Expected Improvement Matrix-Based Infill Criteria for Expensive Multiobjective Optimization. IEEE Trans Evol Computat. 2017;21(6):956–75.
- 35. Sestito JM, Kempner M, Harris TAL, Zarkadoula E, Wang Y. Development of Aluminum Scandium Nitride Molecular Dynamics Force Fields with Scalable Multi-Objective Bayesian Optimization. JOM. 2022;74(9):3487–97.
- 36. Zitzler E, Deb K, Thiele L. Comparison of multiobjective evolutionary algorithms: empirical results. Evol Comput. 2000;8(2):173–95. pmid:10843520
- 37. Fonseca CM, Fleming PJ. Multiobjective optimization and multiple constraint handling with evolutionary algorithms - Part I: A unified formulation. IEEE Trans Syst Man Cybern A. 1998;28(1):26–37.
- 38.
Deb K, Thiele L, Laumanns M, Zitzler E. Scalable test problems for evolutionary multiobjective optimization. Evolutionary Multiobjective Optimization. 2005. p. 105–45.
- 39. Ray T, Liew KM. A Swarm Metaphor for Multiobjective Design Optimization. Eng Optim. 2002;34(2):141–53.
- 40.
Frazier PI. A Tutorial on Bayesian Optimization. 2018;7. Available from: http://arxiv.org/abs/1807.02811
- 41. Snoek J, Larochelle H, Adams RP. Practical Bayesian Optimization of Machine Learning Algorithms. Adv Neural Inform Process Syst. 2012;4:2951–9.
- 42. Park J-S. Optimal Latin-hypercube designs for computer experiments. J Stat Plan Infer. 1994;39(1):95–111.
- 43. Wang GG. Adaptive Response Surface Method Using Inherited Latin Hypercube Design Points. J Mech Design. 2003;125(2):210–20.
- 44.
Wang L, Ng AHC, Deb K. Multi-objective optimisation using evolutionary algorithms: An introduction. Multi-objective Optimisation Using Evolutionary Algorithms. Springer London; 2011. p. 3–34.
- 45. Coello CA. Use of a self-adaptive penalty approach for engineering optimization problems. Comput Ind. 2000;41(3):113–27.