Philippe Callac

Rana Kabalan

Cyril Férandon

Gérard Barroso

Mitochondrial Inheritance During Outcrossing In a Substrate Simultaneously Inoculated With Spores And Mycelium Of Agaricus bisporus

Edinbrugh

Mitochondria inheritance in hybrid strains is uniparental and appears to depend on the used outcrossing method. When two homokaryons are mated in vitro, one of the two parental haplotypes is preferentially inherited. When crosses are performed in vivo, Buller phenomenon occurs and data show that one nucleus from the heterokaryotic spores (or the mycelium issued from them) migrates alone in the homokaryotic mycelium.

In a simple model where mycelia from spores (6% homokaryotic) would cross at random, 97% of the crosses would occur between homokaryons and heterokaryons. This would lead to hybrids without preferential haplotype inheritance (when homoA x heteroB and heteroA x homoB equally occurred). This agrees with the numerous different haplotypes found in wild populations.

DISCUSSION

sporeprint

PARENT A strain Bs 518 (brown cultivar C9) haplotype mt1 (genotype Prm1-1)

Agaricus bisporus is amphithallic: among the offspring, about 94% of the single spore isolates were heterokaryotic, while 6% were homokaryotic.

Spore suspension (about one billion

  INRA UR1264 Mycologie et sécurité des aliments BP 81,33883 Villenave d'Ornon, France email: callac@bordeaux.inra .fr 2-Université Victor Ségalen Bordeaux 2 Two primers: Mu and Mr Allele Prm1-1 : one intron I1 (≈ 1kb) Allele Pm1-2 : two introns I1'+I1'' (


Claire Rempenault

Julie Mielle

Kristina Schreiber

Pierre Corbeau

Laurence Macia

Bernard Combe

Jacques Morel

Claire Immediato Daien

Rachel Audo

CXCR5/CXCL13 pathway, a key driver for migration of regulatory B10 cells, is defective in patients with rheumatoid arthritis

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.


domain: spi phys math

text:
J H Bussemaker
email: jasper.bussemaker@dlr.de

N Bartoli

T Lefebvre
email: thierry.lefebvre@onera.fr

P D Ciampa
email: pier.ciampa@dlr.de

B Nagel

Effectiveness of Surrogate-Based Optimization Algorithms for System Architecture Optimization

The design of complex system architectures brings with it a number of challenging issues, among others large combinatorial design spaces. Optimization can be applied to explore the design space, however gradient-based optimization algorithms cannot be applied due to the mixed-discrete nature of the design variables. It is investigated how effective surrogate-based optimization algorithms are for solving the black-box, hierarchical, mixed-discrete, multiobjective system architecture optimization problems. Performance is compared to the NSGA-II multi-objective evolutionary algorithm. An analytical benchmark problem that exhibits most important characteristics of architecture optimization is defined. First, an investigation into algorithm effectiveness is performed by measuring how accurately a known Pareto-front can be approximated for a fixed number of function evaluations. Then, algorithm efficiency is investigated by applying various multi-objective convergence criteria to the algorithms and establishing the possible trade-off between result quality and function evaluations needed. Finally, the impact of hidden constraints on algorithm performance is investigated. The code used for this paper has been published.

Nomenclature

I. Introduction

T increasing complexity of aerospace products brings several challenges to the design of new systems [1]. Many requirements and constraints on the transportation system, the manufacturing system, the aircraft itself, individual systems and components have to be taken into account. In addition, ever stricter regulations are applied to emission and noise levels. An important solution for overcoming such challenges lies in the ability to take important decisions earlier on in the design process. Such decisions have a high impact on the performance of the final design, and suffer from uncertainty owed to the lack of knowledge at early design stages. Additionally, due to the complexity of aerospace systems, many interacting decisions have to be taken into account at the same time, leading to extremely large combinatorial design spaces. To be better able to take such architecting decisions with an increased level of certainty, it is needed to apply physics-based systematic design space exploration techniques [START_REF] Bussemaker | System Architecture Design Space Exploration: An Approach to Modeling and Optimization[END_REF]. This also reduces the reliance on expert judgment, which, while quick, can also suffer from expert bias, subjectivity, conservatism, and overconfidence.

First, it should be made possible to quantitatively evaluate the performance of candidate architectures, so they can be objectively compared to each other in terms of design goals. This can be a complicated multi-disciplinary, multi-organizational process involving many communication and socio-technical challenges [START_REF] Ciampa | A MBSE Approach to MDAO Systems for the Development of Complex Products[END_REF]. It requires the ability of quickly formulating and deploying design systems to support the quantitative multi-disciplinary analysis of candidate architectures and make sure such analysis is coherent, consistent, and meaningful [START_REF] Ciampa | AGILE Paradigm: The next generation of collaborative MDO for the development of aeronautical systems[END_REF].

Once it is possible to evaluate architecture candidates, it is then needed to model the system architecture design space in such a way as to explicitly identify all architecture design choices that need to be evaluated, and yield function-to-form mappings for the justification of the system architectures with respect to system stakeholders and their needs. In addition to enabling integration within the wider MBSE process, such a model offers the formalization needed to construct an optimization problem based on the identified design choices. An optimization algorithm can then explore the design space to find the most optimal system architecture, or the Pareto-optimal set of system architectures if multiple objectives are optimized for.

System architecture optimization problems are challenging to solve due to their black-box, hierarchical, mixed-integer and multi-objective nature (see Sec. II.B). Research is needed into what optimization algorithms can deal with this kind of problems, and to identify areas where further research might be necessary. This paper addresses this need, by comparing the performance of different optimization algorithms for solving a system architecture optimization problem. An overview is presented on system architecting, its place in the systems engineering process, and behavior of system architecture optimization problems. Then, a literature review of currently existing optimization algorithm components is given. Investigations are performed into the effectiveness (how well the "real" Pareto front is approximated) and the efficiency (how many function evaluations are needed for a certain approximation of the Pareto front). Finally, an investigation into architecture optimization subject to hidden constraints is presented.

II. System Architecting

A. Complex System Architecting

The development of modern complex systems needs to account for an ever increasing number of capabilities to be delivered, as well as for organizational boundaries, integration and communication challenges, and constraints stemming from all the stages of the product's life-cycle. Ideally every decision taken at each stage of the development should be evaluated along the entire life-cycle. The management of such development complexity requires a shift to a novel system development paradigm.

In this context, the DLR "Institute of System Architectures in Aeronautics" is developing a novel "model based conceptual framework" for architecting, designing and optimizing complex aeronautical systems. The expected impact is a drastic reduction in time and costs associated with the development, via an increased transparency, efficiency, and traceability of the design and decision making processes. The conceptual framework, introduced in [START_REF] Ciampa | Accelerating the Development of Complex Systems in Aeronautics via MBSE and MDAO: a Roadmap to Agility[END_REF] extends the scope of design and optimization methods to all the phases of the development life cycle of complex systems. The implementation of the concept is supported by the development of novel design methods and approaches, leveraging digital design engineering and modeling technologies. The work presented in this paper focuses on solving system architecture optimization problems, and is part of the European Commission funded project AGILE 4.0 (2019-2022) [START_REF]AGILE4.0[END_REF]. The project will act as test environment for architecture optimization, and several realistic system architecting problems will be defined and explored within the scope of the project.

System architecting can be seen as a step in the systems engineering process [START_REF] Crawley | System architecture: strategy and product development for complex systems[END_REF]. The system of interest consists of a set of interconnected components that together fulfill the functions that justify the existence of the system. A system architecture describes how these functions are fulfilled, by iteratively mapping elements of form (i.e. architecture components) to function, and defining connections among elements of form. This paper builds on previous work on modeling the design space of system architectures as presented in [START_REF] Bussemaker | System Architecture Design Space Exploration: An Approach to Modeling and Optimization[END_REF]. Modeling the system architecture design space enables the semantic representation of architecture choices involved in the construction of architectures. This enables integration in the MBSE process, and offers the formalization needed for formulating an optimization problem from the design space model.

B. System Architecture Optimization

Optimization is the formalization of a design task: an automated way of finding the best design to solve a given problem [START_REF] Eligius | Introduction to Nonlinear and Global Optimization[END_REF]. In optimization, an optimization algorithm systematically varies a set of design variables (a design vector), and observes the change in one or more merit functions (or objective functions) this change of design variables elicits. Its goal is then to find a design vector that minimizes (or maximizes) the values of these objective functions, while at the same time making sure that the values of any constraint functions remain within limits. The challenge of solving such a problem is increased by non-linear and a-priori unknown objective and constraint function behavior. Many different algorithms exist for solving optimization problem, with various degrees of effectiveness for different types of optimization problems. Here, effective means either to have a high probability of finding the optimal point at all, or finding (an approximation of) the optimal point in the least amount of steps [START_REF] Eligius | Introduction to Nonlinear and Global Optimization[END_REF]. This last point is important, as the evaluation of one design vector to find out its objective and constraint function values, usually takes several orders of magnitude longer than the internal calculations of the optimization algorithms. This is easy to confirm by thinking of aerodynamic wing shape optimization: one evaluation of a possible shape might involve running high-fidelity CFD (Computational Fluid Dynamics) analysis on a supercomputer.

In this paper we are dealing with the optimization of system architectures, which poses several difficulties for optimization algorithms [START_REF] Bussemaker | System Architecture Design Space Exploration: An Approach to Modeling and Optimization[END_REF]. In optimization, a distinction is made between design variables of continuous type and of discrete type. Discrete design variables can be further broken down into integer (ordering and distance have meaning) and categorical (no notion of order or distance) design variables. System architecture optimization problems might have all of these types of design variables: an example of an integer design variable might be the number of engines per wing, a categorical design variable might be whether to add a winglet or not, a continuous design variable might be the bypass ratio of the engine. This mixed-discrete nature of the design variables makes it more difficult to solve the optimization problem, as it excludes the use of gradient-based optimization algorithms: algorithms that use local derivative information to find the "direction" towards the optimum. Another characteristic of system architecture optimization is that because of conflicting stakeholder needs, there are in general multiple objectives that a system architecture might be optimized for. This means that rather than finding one optimal architecture, the goal is to find a set of Pareto-optimal architectures: architectures that are better than others for some objectives, but worse for other objectives [START_REF] Miettinen | Nonlinear Multiobjective Optimization[END_REF]. Within the Pareto set, no architecture can be said to be better than the other, and therefore an additional decision-making step that chooses the final architecture based on additional considerations needs to be performed after the optimization has been completed. Solving these so-called multi-objective optimization problems is hard, because instead of simply moving towards the optimal point, there is also a need to find design points sufficiently spread-out along the Pareto-front (i.e. the set of Pareto points in the solution-space). This is called the balance between exploitation (i.e. finding the best design) and exploration (i.e. finding new regions in the design space of potential best designs), the principle is visualized in Fig. 1.

Next, system architecture optimization problems feature decision hierarchy: design variables can be conditionally active based on other design variables. This means that there are regions in the design space where one or more of the design variables have no influence of the performance of the architecture. There are several ways of dealing with this problem [START_REF] Zaefferer | A First Analysis of Kernels for Kriging-Based Optimization in Hierarchical Search Spaces[END_REF] (see Fig. 2 for a visual explanation): ignore the effects, imputation, and explicit consideration. Ignoring the effects potentially confuses the optimization algorithm, because multiple different design vectors might actually map to the same design point (because of inactive design variables). This is solved by imputation: inactive design variables are replaced by some predefined value, so that no duplicate design vectors are evaluated. Explicitly considering hierarchy effects can yield the most effective optimization algorithms, however this requires special-purpose optimization algorithms, and information about which design variables are active at locations in the design space (the 𝛿-function of [START_REF] Hutter | A Kernel for Hierarchical Parameter Spaces[END_REF]). It should be noted that regardless of the strategy chosen, there will always be a difference between the apparent and feasible design space [START_REF] Bussemaker | System Architecture Design Space Exploration: An Approach to Modeling and Optimization[END_REF]: the apparent design space spans all combinations of all design variables (i.e. what the optimization algorithm sees), the feasible design space takes out combinations leading to infeasible or duplicate design points due to inactive design variables (i.e. what the architecture evaluator sees).

Finally, because quantitative performance evaluation of systems depends on numerical simulation techniques, architecture optimization problems may also be subject to hidden constraints. Hidden constraints are constraints that are not known to the optimizer, and only the status (i.e. satisfied or violated) is known qualitatively after completing an evaluation [START_REF] Mãijller | Surrogate Optimization of Computationally Expensive Black-Box Problems with Hidden Constraints[END_REF]. They are especially common in simulation-based optimization, because it can happen that simulations do not converge to a meaningful result [START_REF] Digabel | A Taxonomy of Constraints in Simulation-Based Optimization[END_REF]. This can also mean that if such constraint is violated, none of the other objective and constraint values has a valid result. Up to 60% of evaluations in a design of experiments may result in violated hidden constraints (i.e. infeasible results) [START_REF] Mãijller | Surrogate Optimization of Computationally Expensive Black-Box Problems with Hidden Constraints[END_REF]. In general, system architecture optimization problems can therefore be considered to have the following features:

• The design space is non-linear and its shape is unknown a-priori (the evaluation function is a "black-box"),

• Design variables are of mixed-discrete nature,

• A hierarchy between design variables may be present,

• The optimization problem may have multiple objectives to optimize for, and • Hidden constraints may exist in the design space.

III. Optimization Algorithm Components

This section will provide an overview of existing types of optimization algorithms that in principle can tackle mixed-integer, multi-objective, hierarchical, black-box optimization problems: the type of problem involved in system architecture optimization. Additionally, multi-objective convergence criteria are discussed, as detecting the convergence of the Pareto front towards the real Pareto front is a non-trivial task compared to detecting single-objective convergence.

A. Evolutionary Algorithms

Generally a good choice for solving mixed-discrete optimization problems are Evolutionary Algorithms (EA's), originally known as the Genetic Algorithm (GA). EA's are based on the natural selection process as observed by Darwin [START_REF] Glover | [END_REF]: an initial population of individuals (design points) is evolved (improved) by generating new generations (new design points) to replace parts of the current generation: parents are selected from the current population (current best design points) and offspring (new design points) is created by applying different kinds of genetic operators, like crossover and mutation. The offspring is then compared to the current population and used to replace all or parts of the current generation, based on comparing their fitness (objective) values and optionally retaining the current best individuals from the current population (elitism). This process is visualized in Fig. 3. Evolutionary algorithms are well suited for finding global optima, as long as the selection and configuration of the different genetic operators allow for sufficient exploitation (improving already-good points) and exploration (trying out new points in the design space).

Evolutionary Algorithms use crossover and mutation to generate offspring from selected parents. Different crossover and mutation operators are available. Additionally, different crossover and mutation operators can be used for continuous than for discrete design variables. Since in architecture optimization, both types of optimization variables will be present in the optimization problem, different types of evolutionary operators will be used on different design variable types. For the experiments in this study, the Simulated Binary Crossover and Uniform Crossover [START_REF] Blank | Pymoo: Multi-Objective Optimization in Python[END_REF] are used for continuous and discrete design variables, respectively. For the mutation operation, Polynomial Mutation and Bitflip mutation [START_REF] Blank | Pymoo: Multi-Objective Optimization in Python[END_REF] are used for continuous and discrete design variables, respectively. Multi-Objective Evolutionary Algorithms (MOEA's) are adaptations of the EA process where mainly the selection and replacement steps are modified to result in better exploration of the multi-objective search space [START_REF] Zitzler | Comparison of multiobjective evolutionary algorithms: empirical results[END_REF]. Instead of optimizing for one objective directly, rather the population is improved by binary comparison according to some Pareto dominance relation (see Fig. 4), which results in an improvement of design points towards to Pareto-front, whereas at the same time a proper spread along the Pareto-front is ensured. Selecting individuals for creating offspring is done based on Pareto-rank as well. The most used MOEA is the Nondominated Sorting Genetic Algorithm II (NSGA-II) [START_REF] Deb | A fast and elitist multiobjective genetic algorithm: NSGA-II[END_REF]. This algorithm is based on ranking the design points based on their distance from the Pareto-front (exploitation), closer is better, and the distance to other neighboring points (exploration), further is better.

Evolutionary Algorithms are effective at dealing with hierarchical optimization problems using the hidden genes approach [START_REF] Nyew | Structured-Chromosome Evolutionary Algorithms for Variable-Size Autonomous Interplanetary Trajectory Planning Optimization[END_REF]. The definition of the design space (the chromosomes) includes the maximum number of design variables (genes). Inactive design variables (hidden genes) are excluded from the cost function evaluation, however are still visible and actively take part in the design space exploration (e.g. evolutionary operators). There is therefore a difference between what the evolutionary algorithm operators on and what is used for cost function evaluation: this is the same phenomena of apparent vs feasible design spaces as discussed before.

Since Evolutionary Algorithms depend on comparison between results, for example through non-dominated sorting in the case of NSGA-II, the extreme barrier approach can be used to deal with hidden constraints [START_REF] Gratton | A Merit Function Approach for Direct Search[END_REF]: objectives and constraints are assigned a value of +∞ in case of violated hidden constraints, so that selection and replacement are guided towards feasible design points.

B. Surrogate-Based Optimization

Evolutionary algorithms are known for their robust effectiveness for finding global optima, but do require a relatively high number of function evaluations (however still many less than design space enumeration). This is a potentially problematic point, as it is expected that the main computation cost lies exactly there: at the function evaluations. Additionally, because simulation becomes more high-fidelity, this problem will only grow in the future. One class of optimization algorithms specifically designed to mitigate this problem are Surrogate-Based Optimization (SBO) algorithms [START_REF] Chugh | A survey on handling computationally expensive multiobjective optimization problems with evolutionary algorithms[END_REF]. Efficient Global Optimization (EGO) was one of the first SBO [START_REF] Jones | Efficient Global Optimization of Expensive Black -Box Functions[END_REF], based on Kriging model and the Expected Improvement (EI) infill criterion. The process, visualized in Fig. 5, is the same for all SBO algorithms: create an initial sample to build the surrogate model, sample new points based on an infill criterion and update the surrogate model, check if our convergence criterion has been reached, and sample new points if convergence has not been reached yet. The infill criterion is used to select new points to sample, and usually includes some trade-off between exploration and exploitation. For example, if Kriging surrogate models are used, model uncertainty information is used to aid in exploration (to explore areas to improve the current estimate of the surrogate model, in the hope that a new area of global optimum is found [START_REF] Tian | Multiobjective Infill Criterion Driven Gaussian Process-Assisted Particle Swarm Optimization of High-Dimensional Expensive Problems[END_REF]). Exploitation is aided by simply selecting points to evaluate where the surrogate model predicts a global optimum. The great advantage of this approach is that because evaluating a point in the surrogate model is extremely cheap compared to evaluating the real function, an exhaustive search of the design space to find the best infill points can be performed. Training a surrogate model, however, can take a non-trivial amount of time and therefore should be taken into account.

The two most popular interpolating surrogate models used for SBO are Radial Basis Function (RBF) and Kriging surrogate models [START_REF] Chugh | A survey on handling computationally expensive multiobjective optimization problems with evolutionary algorithms[END_REF]. The main difference is that next to the function estimate itself, Kriging models (also known as Gaussian process models [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF]) also calculate the uncertainty of the estimate. This allows the development of infill criteria that not only attempt to find the global minima (exploitation), but also attempt to improve the estimate of the surrogate model itself as to better be able to predict areas with global minima (exploration). Kriging models, however, have difficulty training for high-dimensional design spaces, which is less of a problem for RBF models [START_REF] Datta | A surrogate-assisted evolution strategy for constrained multi-objective optimization[END_REF]. Kriging modified with a Partial Least Squares estimation (KPLS) increases the training speed of Kriging models and might be a powerful mechanism for dealing with high-dimensional spaces [START_REF] Bouhlel | Efficient global optimization for high-dimensional constrained problems by using the Kriging models combined with the partial least squares method[END_REF]. SBO algorithms are very promising for efficiently finding the global optimum, however due to the nature of the used surrogate models are usually better suited for modeling continuous design spaces. A Kriging surrogate model uses kernels to calculate correlations between sample points, such that points that lie closer to each other are assumed to be more correlated to each other than points that lie further away. A common kernel is the exponential kernel, where correlation exponentially reduces with distance. The problem with mixed-integer optimization problems is that for the categorical variables, there is no notion of distance. Several different ways to solve this problem have been studied in the past, including using hybrid branch-and-bound algorithms [START_REF] Roy | A Mixed Integer Efficient Global Optimization Algorithm for the Simultaneous Aircraft Allocation-Mission-Design Problem[END_REF], treating discrete variables as continuous variables [START_REF] Müller | MISO: mixed-integer surrogate optimization framework[END_REF], continuous relaxation or dummy coding [START_REF] Garrido-Merchán | Dealing with categorical and integer-valued variables in Bayesian Optimization with Gaussian processes[END_REF], and developed special-purpose mixed-integer kernels [START_REF] Pelamatti | Efficient global optimization of constrained mixed variable problems[END_REF][START_REF] Zuniga | Global optimization for mixed categorical-continuous variables based on Gaussian process models with a randomized categorical space exploration step[END_REF].

Infill Criteria Infill criteria determine the place in the design space where new points are sampled to increase the prediction quality of the model, and to confirm regions with global minima. The simplest infill criteria are based on directly finding either the design point where the surrogate model predicts the best objective function value ŷ (bias towards exploitation, e.g. [START_REF] Gutmann | A Radial Basis Function Method for Global Optimization[END_REF]), or where the model uncertainty (variance ŝ) is highest (bias towards exploration, e.g. [START_REF] Gramacy | Parameter space exploration with Gaussian process trees[END_REF]) [START_REF] Horn | Model-Based Multi-objective Optimization: Taxonomy, Multi-Point Proposal, Toolbox and Benchmark[END_REF]. Much effort has been invested in finding infill criteria with more natural balances between exploration and exploitation. The Expected Improvement (EI) [START_REF] Jones | Efficient Global Optimization of Expensive Black -Box Functions[END_REF] is a more balanced metric and represents the amount the function can be improved with respect to the current best point given that an improvement is achieved, and is calculated by:

𝐸 𝐼 (𝑥) = ( 𝑓 𝑚𝑖𝑛 -ŷ(𝑥)) • Φ 𝑓 𝑚𝑖𝑛 -ŷ(𝑥) ŝ(𝑥) + ŝ(𝑥) • 𝜙 𝑓 𝑚𝑖𝑛 -ŷ(𝑥) ŝ(𝑥)
, where 𝑓 𝑚𝑖𝑛 is the function value of the current best sample, ŷ and ŝ are the function and variance estimates, respectively, Φ and 𝜙 the standard normal density and distribution functions, respectively. The infill criterion then selects the point with the highest 𝐸 𝐼 as new point to sample. Another criterion uses the Lower Confidence Bound (LCB) [START_REF] Cox | A statistical method for global optimization[END_REF], which is the lowest expected value to be observed at a point given its normal distribution, and is calculated by: where 𝛼 is a scaling parameter (a typical value is 2 [START_REF] Cox | A statistical method for global optimization[END_REF]). The LCB then replaces the objective function value and can be directly optimized for (minimized) to find a point with a potential new best function value. Another metric, the Probability of Improvement (PoI) [START_REF] Hawe | An Enhanced Probability of Improvement Utility Function for Locating Pareto Optimal Solutions[END_REF], gives the probability that some value with an offset to the current best value 𝑇 < 𝑓 𝑚𝑖𝑛 is improved upon given the model estimate and its normal distribution (visualized in Fig. 6):

𝐿𝐶 𝐵(𝑥) = ŷ(𝑥) -𝛼 • √ ŝ(𝑥),
𝑃𝑜𝐼 (𝑥) = Φ 𝑇 -ŷ(𝑥) ŝ(𝑥)
The point with the maximum PoI can then be used for sampling the function.

If an RBF surrogate model is used, only ŷ is available and therefore such balanced infill criteria as discussed before cannot be used. One strategy to still enable exploration of the design space is to use Coordinated Perturbation (CP) [START_REF] Regis | Combining radial basis function surrogates and dynamic coordinate search in high-dimensional expensive black-box optimization[END_REF], a method similar to Simulated Annealing (SA), but where the probability of evaluating a new point is based on a trade-off between the distance to the current point (exploration) and the predicted function value (exploitation). CP modified for use in multi-objective optimization will be called the ŷ-Dist criterion, and implements the following differences with respect to the single-objective CP criterion:

1) Objective values are normalized and the euclidean distance is used to calculate the minimum distance (see [START_REF] Regis | Combining radial basis function surrogates and dynamic coordinate search in high-dimensional expensive black-box optimization[END_REF]);

2) The RBF criterion and minimum distance are defined as separate objectives in the infill problem, to naturally balance between exploration and exploitation; this also removes the need for normalization between the RBF criterion and the minimum distance.

Constraint Handling Constraint handling in Surrogate-Based Optimization can be approached in several different ways [START_REF] Priem | Upper trust bound feasibility criterion for mixed constrained Bayesian optimization with application to aircraft design[END_REF]. One common way is to transform the constrained problem into an unconstrained problem by modifying the objective to be penalized due to constraint violation, thereby driving the optimizer towards the feasible region. For this approach it is difficult to find a balance between the objective value and the penalty multiplier. Another way is to drive the infill criterion to ignore points where the surrogate model for the constraint functions expect the constraint to be violated. This, however, can suffer from estimation errors and potentially ignore points lying in a feasible but unexplored area of the design space. One way to solve this is the Expected Violation (EV) metric, that includes information about the model variance to only ignore new infill points if this value (analogue to EI) is above some predefined threshold. The Probability of Feasibility (PoF) metric works according to the same principle [START_REF] Schonlau | Global versus local search in constrained optimization of computer models[END_REF].

Dealing with hidden constraints is more problematic, since violated hidden constraints result in all outputs (i.e. objectives and constraints) being invalid. The extreme barrier approach (replacing invalid values with +∞) is not applicable, as infinite values cannot be modeled in surrogate models. One way of dealing with this is by replacing invalid values with the maximum of valid values, so that the surrogate model will predict unattractive objective and constraint values in the region of violated hidden constraints [START_REF] Mãijller | Surrogate Optimization of Computationally Expensive Black-Box Problems with Hidden Constraints[END_REF]. Another approach is by removing invalid design points from the surrogate model training set, and creating a second model that predicts whether a design point will be feasible or not [START_REF] Lee | Optimization Subject to Hidden Constraints via Statistical Emulation[END_REF]. The second predictor might then be based on some classification model, for example random trees.

Multi-Objective, Mixed-Integer and Hierarchical Problems

Dealing with multi-objective optimization problems has also been an important topic of research for SBO algorithms [START_REF] Rojas-Gonzalez | A Survey on Kriging-Based Infill Algorithms for Multiobjective Simulation Optimization[END_REF]. Many strategies have been studied, including multi-objective variations of the Expected Improvement [START_REF] Keane | Statistical Improvement Criteria for Use in Multiobjective Design Optimization[END_REF] and Probability of Improvement (PoI) [START_REF] Hawe | An Enhanced Probability of Improvement Utility Function for Locating Pareto Optimal Solutions[END_REF] criteria, the Expected Hypervolume Improvement (EHVI) [START_REF] Dos Passos | Multi-objective optimization with Kriging surrogates using "moko", an open source package[END_REF], or simply optimizing for the estimated objective values themselves [START_REF] Tian | Multiobjective Infill Criterion Driven Gaussian Process-Assisted Particle Swarm Optimization of High-Dimensional Expensive Problems[END_REF]. Infill criteria for multi-objective optimization are discussed in more details in a later section.

Extending Kriging surrogate models to support mixed-integer and hierarchical (i.e. variable-size) design spaces is an ongoing area of research, Tab. 1 provides an overview of existing approaches, and good reviews are provided by [START_REF] Zaefferer | A First Analysis of Kernels for Kriging-Based Optimization in Hierarchical Search Spaces[END_REF][START_REF] Pelamatti | Overview and Comparison of Gaussian Process-Based Surrogate Models for Mixed Continuous and Discrete Variables: Application on Aerospace Design Problems[END_REF]. One method is to simply train a surrogate model of the continuous design variables for each of the combinations of the discrete variables, known as Category-Wise (CW) Kriging [START_REF] Pelamatti | Surrogate model based optimization of constrained mixed variable problems: application to the design of a launch vehicle thrust frame[END_REF]. The problem with this approach is that for the optimization of expensive black-box functions there are often not enough samples available to construct an accurate surrogate model for each of the discrete variable combinations [START_REF] Pelamatti | Bayesian optimization of variable-size design space problems[END_REF]. An alternative approach is to model the design space including the continuous and discrete variables in one surrogate model. An obvious possibility is to treat the discrete variables as continuous variables, also called continuous relaxation; this method might be more appropriate for integer design variables, and less for categorical design variables. Kriging kernels capable of modeling mixed-integer spaces include the Hamming distance [START_REF] Hutter | Automated Configurations of Algorithms for Solving Hard Computational Problems[END_REF], Gower distance [START_REF] Halstrup | Black-box optimization of mixed discrete-continuous optimization problems[END_REF], and Symbolic Covariance (SC) [START_REF] Mccane | Distance functions for categorical and mixed variables[END_REF], Compound Symmetry (CS) [START_REF] Roustant | Group kernels for Gaussian process metamodels with categorical inputs[END_REF], and Latent Variables (LV) [START_REF] Zhang | A Latent Variable Approach to Gaussian Process Modeling with Qualitative and Quantitative Factors[END_REF] kernels. Kernels that can also handle hierarchical design spaces include the Arc [START_REF] Hutter | A Kernel for Hierarchical Parameter Spaces[END_REF], Indefinite Conditional [START_REF] Zaefferer | A First Analysis of Kernels for Kriging-Based Optimization in Hierarchical Search Spaces[END_REF], Imputation [START_REF] Zaefferer | A First Analysis of Kernels for Kriging-Based Optimization in Hierarchical Search Spaces[END_REF], Wedge [START_REF] Horn | Surrogates for hierarchical search spaces[END_REF], Sub-Problem-Wise (SPW) decomposition [START_REF] Pelamatti | Bayesian optimization of variable-size design space problems[END_REF], and Dimensional Variable-Wise (DVW) decomposition [START_REF] Pelamatti | Bayesian optimization of variable-size design space problems[END_REF] kernels. It is noted that Baert et al. [START_REF] Baert | Surrogate-Based Optimisation for a Mixed-Variable Design Space: Proof of Concept and Opportunities for Turbomachinery Applications[END_REF] mention that they have implemented mixed-integer Kernels inspired by Value Different Metrics (VDM) [START_REF] Wilson | Improved Heterogeneous Distance Functions[END_REF]. However, the VDM and derived mixed-integer metrics depend on the relative number of occurrences of output labels, and are therefore only applicable to classification problems (discrete output), and not to regression problems (continuous output) as with SBO.

An overview of SBO approaches including research into applying SBO for mixed-integer or multi-objective optimization is shown in Tab. 2. As can be seen, there has been no research on applying SBO for mixed-integer multi-objective problems simultaneously. Due to the modular nature of SBO algorithms, it should be possible to combine mixed-integer surrogate models with multi-objective infill criteria to enable this.

Infill Criteria for Multi-Objective Optimization Problems

Two main approaches can be identified for creating infill criteria for use with multi-objective optimization problems: modification of single-objective (SO) infill criteria, and the formulation of infill criteria specifically made for multi-objective (MO) optimization.

Single-objective Infill Criteria in Multi-Objective Optimization

Existing single-objective infill criteria, like the Lower Confidence Bounds (LCB) criterion, can be used in MO optimization by running a multi-objective optimization on infill criteria evaluated for each objective separately [START_REF] Tian | Multiobjective Infill Criterion Driven Gaussian Process-Assisted Particle Swarm Optimization of High-Dimensional Expensive Problems[END_REF]. This would then result in a Pareto-set of points balancing the infill criteria values for the objectives, naturally resulting in a set of infill points that balance the exploration and exploitation over the difference objectives. This principle can be used for any SO infill criterion that does not depend on the notion of current best objective value.

If an infill criterion depends on the current best objective values, for example the Expected Improvement (EI) or Probability of Improvement (PoI), a modification is needed to make them usable in MO optimization. This modification is needed, because if these criteria would simply be evaluated for each objective separately, only points near the edges of the current Pareto front would be selected for infill, as these points represent the respective best values of the objectives. One way to solve this problem is by scalarizing the optimization problem for various weighting vectors and apply the SO infill criteria on the scalarized problem. This approach is for example used by ParEGO [START_REF] Knowles | ParEGO: a hybrid algorithm with on-line landscape approximation for expensive multiobjective optimization problems[END_REF]: objective weightings are randomly selected from evenly distributed weightings, and scalarization is then done using the augmented Tchebycheff function.

Another way to solve this problem is to instead of using the current best objective value, to use the objective values of the closest point on the Pareto front. This principle is used in the experiments in this paper, and it is shown in Fig. 7: it can be seen that around each point currently on the Pareto front, a slope towards the improvement of the Pareto front is constructed. For the shown hypothetical bi-objective problem, these two infill criterion objectives will therefore result in a well-diversified selection of infill points along the existing Pareto front. Compared to ParEGO, this approach prevents the scalarization of the multi-objective problem.

Special-Purpose Multi-Objective Infill Criteria

Another strategy is by using infill criteria specifically designed for multi-objective (MO) problems. Most of these MO infill criteria attempt to transform the MO objective space into a single-objective space related to the distance to the current Pareto front. The goal is to find infill points that maximize some probability of finding an improvement over the current Pareto front. Table 3 compares the different MO infill criteria. Figure 8 visualizes the infill criteria objective spaces.

One of the earliest available criteria was the Expected Hypervolume Improvement (EHVI), as the hypervolume is a natural measure of the extent of the Pareto front. The problem with the EHVI is that there is no closed form for computing the expected improvement of the hypervolume, and that for dealing with the stochastic nature of the objective space (in terms of the mean predicted value and the variance), a Monte Carlo sampling approach is used. Computing the hypervolume itself is already expensive, and coupled with a Monte Carlo approach this leads to a factor (b) Euclidean Expected Improvement (EEI) [START_REF] Keane | Statistical Improvement Criteria for Use in Multiobjective Design Optimization[END_REF] (c) Enhanced Probability of Improvement (EPoI) [START_REF] Hawe | An Enhanced Probability of Improvement Utility Function for Locating Pareto Optimal Solutions[END_REF] (d) Expected Maximin Fitness Improvement (EMFI) [START_REF] Svenson | Multiobjective optimization of expensive-to-evaluate deterministic computer simulator models[END_REF] (e) Minimum Probability of Improvement (MPoI) [START_REF] Rahat | Alternative infill strategies for expensive multi-objective optimisation[END_REF] (f) Minimum Euclidean Probability of Improvement (MEPoI) of approximately 500 between the computation of the EHVI and the MPoI. Another difficulty of the EHVI is that since the calculation of the HV does not extend beyond the edges of the currently known Pareto front, the EHVI criterion does not indicate improvement potential near these edges. This effect can be seen in Fig. 8a. An additional consideration is the dependence of the computation time on the number of points in the Pareto front, as this might vary a lot between optimization runs, and the size of the Pareto front can grow very large if it is sufficiently diversified. Table 3 shows that over time, better MO infill criteria have been found that are quicker to compute, and have less dependence on the number of points in the existing Pareto front. The Minimum Variance of Pareto Front (MVPF) infill criterion derives its advantage from the selection of infill points after running a normal infill optimization based on the mean predicted objective values ( ŷ). Therefore, the computation of the MVPF infill criterion actually costs no time other than querying the surrogate model.

MO extensions of the PoI and EI criteria are the Enhanced Probability of Improvement (EPoI) and Euclidean Expected Improvement (EEI) criteria, respectively. Similarly to their SO criteria, the EPoI and EEI are both based on the probability of improvement over the current Pareto front. The EEI then takes this probability and multiplies it with the centroid of the probability integral. The original derivation by Keane [START_REF] Keane | Statistical Improvement Criteria for Use in Multiobjective Design Optimization[END_REF] shows a closed-form integral for the calculation of the EEI. However, due to computational cost concerns, the implementation used here replaces this integral with a Monte Carlo sampling approach, based on [START_REF] Hawe | An Enhanced Probability of Improvement Utility Function for Locating Pareto Optimal Solutions[END_REF]. This also enables the use of these criteria for MO optimization problems with more than two objectives. Additionally, for the EEI, the closed-form centroid calculation has been replaced by the euclidean-based improvement measure from [START_REF] Parr | Improvement Criteria for Constraint Handling and Multiobjective Optimization[END_REF] to enable use in higher-dimensional MO problems. Both the EPoI and EEI criteria suffer from the same limitation as the EHVI criterion, in that near the edges of the Pareto front no improvement is indicated, as shown in Fig. 8b and8c.

The implementation of the Expected Maximin Fitness Improvement EMFI) criterion is based on [START_REF] Svenson | Multiobjective optimization of expensive-to-evaluate deterministic computer simulator models[END_REF]. It deals with the predicted variance ŝ by performing Monte Carlo sampling, leading to the non-smooth infill objective space as shown in Fig. 8d. In this case it is important to select an algorithm that can handle non-smooth and/or discontinuous objective spaces for solving the infill optimization problem.

The Minimum Probability of Improvement (MPoI) criteria has been [START_REF] Rahat | Alternative infill strategies for expensive multi-objective optimisation[END_REF] extended using the euclidean-based improvement measure from [START_REF] Parr | Improvement Criteria for Constraint Handling and Multiobjective Optimization[END_REF] as the Minimum Euclidean Probability of Improvement (MEPoI) metric. The objective spaces of the MPoI and MEPoI criteria are shown in Fig. 8e and8f.

One issue with MO infill criteria, is that the MO infill search space is single-objective: there is only one point that maximizes a given MO infill criterion. It is possible to only search for one infill point per algorithm step, however then there is no guarantee that these new infill points are well-diversified along the existing Pareto front, especially since it is recognized that the infill criterion search landscape often is highly multi-model [START_REF] Rahat | Alternative infill strategies for expensive multi-objective optimisation[END_REF]. To enable the selection of multiple MO infill points per algorithm step, the single MO infill objective can be modulated by the predicted objective values to turn the single-objective MO infill criterion problem into a multi-objective infill problem. This principle is visualized in Fig. 9: an infill criterion (shown is the Enhanced Probability of Improvement [START_REF] Hawe | An Enhanced Probability of Improvement Utility Function for Locating Pareto Optimal Solutions[END_REF]) is modulated over the two underlying problem objectives to turn the single-objective EPoI infill problem into a bi-objective problem, corresponding to the number of objectives of the underlying problem.

C. Multi-Objective Convergence Criteria

Normally, EA's and SBO's are executing with an evaluation budget [START_REF] Martí | A stopping criterion for multi-objective optimization evolutionary algorithms[END_REF]. This means they run for a given amount of function evaluations and stop once the budget has been exceeded, and then choose the best design point(s) from its search history to represent the best design points. Computational cost, however, is potentially wasted if the optimization algorithm converges before the budget has been exhausted: the algorithm then continues to search the design space, but no new better design points are discovered. Detecting convergence in single-objective optimization is usually done based on the value of the objective function itself: convergence is then assumed once objective function improvement per iteration falls below some threshold [START_REF] Rudenko | A Steady Performance Stopping Criterion for Pareto-based Evolutionary Algorithms[END_REF]. Applying the same principle to multi-objective optimization is not possible, as there the shape of the Pareto-front is the relevant aspect, not the current best values of the objectives. Note that convergence criteria can be independently applied to both Evolutionary Algorithms and Surrogate Based Optimization algorithms.

Common multi-objective convergence criteria include the Mutual Domination Rate (MDR) [START_REF] Martí | A stopping criterion for multi-objective optimization evolutionary algorithms[END_REF], Fitness Homogeneity Indicator (FHI) [START_REF] Marti | A progress indicator for detecting success and failure in evolutionary multiobjective optimization[END_REF], Generational Distance (GD) [START_REF] Van Veldhuizen | Evolutionary Computation and Convergence to a Pareto Front[END_REF], Spread [START_REF] Deb | A fast and elitist multiobjective genetic algorithm: NSGA-II[END_REF], Hypervolume [START_REF] Zitzler | Multiobjective Evolutionary Algorithms: A Comparative Case Study and the Strength Pareto Approach[END_REF], Consolidation Ratio (CR) [START_REF] Goel | A non-dominance-based online stopping criterion for multi-objective evolutionary algorithms[END_REF], Maximal Crowding Distance (MCD) and Steady Performance Indicator (SPI) [START_REF] Rudenko | A Steady Performance Stopping Criterion for Pareto-based Evolutionary Algorithms[END_REF].

IV. Investigation of Optimization Algorithm Effectiveness

In the previous section evolutionary algorithms, surrogate-based algorithms, and multi-objective convergence criteria that may be effective for solving system architecture optimization problems: black-box, mixed-discrete, hierarchical, multi-objective optimization problems. The first aspect to be investigated in this work is the effectiveness of different algorithms. Effectiveness is seen as distinct from efficiency, and relates to how well an optimization algorithm is able to find the global optimum (for single-objective optimization) or the Pareto-front (for multi-objective optimization), regardless of how quick it converges to the final result (i.e. efficiency) [START_REF] Eligius | Introduction to Nonlinear and Global Optimization[END_REF].

The optimization algorithms and convergence criteria discussed in the previous section have been implemented in Python using the multi-objective optimization framework pymoo * [START_REF] Blank | Pymoo: Multi-Objective Optimization in Python[END_REF]. Surrogate-based optimization capabilities have been implemented using scikit-learn † [START_REF] Pedregosa | Scikit-learn: Machine Learning in Python[END_REF] for Kriging models with mixed-integer and/or hierarchical kernels, and using SMT ‡ [START_REF] Bouhlel | A Python surrogate modeling framework with derivatives[END_REF] for RBF and continuous Kriging surrogate models. The optimization algorithm code and the code for performing the optimization experiments described in this section is published at [78] § . Readers should be able to reproduce the results presented in this paper.

All optimization concepts discussed in this section are implemented in Python using the pymoo multi-objective optimization framework [START_REF] Blank | Pymoo: Multi-Objective Optimization in Python[END_REF]. The code for performing the optimization experiments that forms the basis for this paper is published online, to improve reproducibility.

A. Analytical Benchmark Problem

The analytical benchmark problem that will be used in this paper is based on the hierarchical Rosenbrock problem described by Pelamatti et al. [53]. This problem is based on the well-known high-dimensional Rosenbrock optimization problem, modified to have mixed-discrete and hierarchical design variables: two categorical variables determine the sub-problem the other variables act on (these are called dimensional variables in [START_REF] Pelamatti | Bayesian optimization of variable-size design space problems[END_REF]). Based on which of the 4 sub-problems are selected, different continuous variables are active. In addition, several integer variables are added that map to continuous values for evaluation in by the original Rosenbrock function. The hierarchical Rosenbrock problem features 8 continuous, 3 integer, and 2 categorical design variables, 1 objective, and 2 inequality constraints.

Two modifications are applied to the hierarchical Rosenbrock problem of [START_REF] Pelamatti | Bayesian optimization of variable-size design space problems[END_REF]: the problem is made multi-objective by adding a second objective, and another layer of categorical variables are added to increase the average amount of inactive variables. The second objective is calculated using

𝑓 2 (𝒙) = 10 - 𝑓 1 (𝒙) 40 2 + 4 ∑ 𝑖=1 200 (𝑥 𝑖 + 1) 2 , ( 1
)
where 𝑓 1 is the first objective, 𝒙 is the design vector, 𝑥 𝑖 are the first four continuous design variable as selected in the sub-problem. The Pareto front of the multi-objective hierarchical Rosenbrock problem is shown in Fig. 10a. As can be seen, it is composed of design points in two of the four original sub-problems.

The second modification is the repetition of the design variables, the creation of several random mappings from repeated variables to original variables, and the addition of a categorical variables selecting which mapping to use. The number of repetitions 𝑛 𝑟 𝑒 𝑝 and number of mappings 𝑛 𝑚𝑎 𝑝 are variable. Random mappings are created by randomly selecting which repetition an original design variable should map from for each of the design variables. The random state is fixed based on 𝑛 𝑟 𝑒 𝑝 and 𝑛 𝑚𝑎 𝑝 : this gives the impression of a random mapping, however it fixes the behavior of the problem between optimization runs, as needed for applying statistics on optimization algorithm performance. Additionally, each repetition modifies the objective values: they are moved along the existing Pareto front dimensions, so that a combined Pareto front is constructed that is constructed across repetitions. The result of this modification to the multi-objective hierarchical Rosenbrock problem is shown in Fig. 10b: it can be seen that the Pareto front consists of three parts, across two repetitions and two sub-problems.

The analytical architecture optimization benchmark problem used in this paper is therefore the repeated, multiobjective, hierarchical Rosenbrock problem, constructed with 𝑛 𝑟 𝑒 𝑝 = 2 and 𝑛 𝑚𝑎 𝑝 = 2. In total, this benchmark problem has 2 objectives, 2 inequality constraints, and 27 design variables, of which 16 continuous, 6 integer, and 5 categorical.

B. Experimental Setup

Different optimization algorithms will be compared with each other for their effectiveness at finding the Pareto front in the engine architecting problem. The different Kriging kernels and infill criteria for the Surrogate-Based Optimization (SBO) algorithm make for a large number of possible optimization algorithms. Therefore, before the actual algorithms will be compared with each other, the following two questions will be answered in order to down-select the SBO algorithms:

1) Which Kriging kernel performs best for system architecture optimization problems? 2) Should multi-objective infill criteria be used in their original form or their MO-modulated form?

The down-selected SBO algorithms will then be compared to each other and with the other algorithms for their effectiveness at finding the Pareto-front in the engine architecting design space. Each algorithm will be executed 8 times with the same parameter to obtain a feeling for the impact of stochastic elements in the algorithms (i.e. initial DOE, evolutionary operators, infill search). The initial Design of Experiments (DOE) for constructing the first surrogate models will be based on a Latin Hypercube Sampling (LHS) strategy. Algorithm effectiveness will be measured using the following metrics: Δ𝐻𝑉 [START_REF] Palar | On multi-objective efficient global optimization via universal Kriging surrogate model[END_REF] to measure how closely the Pareto-front is approximated, and spread [START_REF] Deb | A fast and elitist multiobjective genetic algorithm: NSGA-II[END_REF] to measure how evenly points are spread out along the Pareto-front. Constraints are handled by either rejecting newly found points with a constraint violation (for MOEA's), rejecting search points with predicted constraint violation (for RBF SBO), or by using the Probability of Feasibility (PoF) indicator (for Kriging SBO).

Table 4 presents the different algorithms that will be compared: an MOEA and multiple SBO algorithms. The algorithms are chosen such that the impact of different choices can be investigated: MOEA vs SBO, SBO surrogate model type, and infill criterion type. Note that the EHVI (Expected Hypervolume Improvement) infill criterion is not taken into account, due to its excessive computational cost (see also Table 3).

Kriging Surrogate Selection

The Kriging surrogate used in Surrogate Based Optimization (SBO) algorithms supports different kernels for representing similarity between two sampling points. These kernels are then used to construct the covariance matrix used in forming the prediction. As shown in Tab. 1, there are several ways of using Kriging for mixed-integer, hierarchical design spaces. This section aims to answer the questions of which Kriging kernel is the most appropriate for use for system architecting optimization problems. This will be done by comparing the different approaches shown in Tab. 1. Category-Wise Kriging is shown to be ineffective in SBO due to the requirement of enough samples for the remaining continuous design variables for each of the discrete combinations [START_REF] Pelamatti | Bayesian optimization of variable-size design space problems[END_REF], and will therefore not be taken into account. Each of the Kriging surrogates with the different MI and hierarchical design space modeling approaches will be executed on the analytical benchmark problem with a budget of 300 function evaluations, with the Expected Improvement (EI) and Minimum Probability of Improvement (MPoI) infill criteria, and otherwise the same parameters as for the main comparison experiment. Kernels will be compared using the Δ𝐻𝑉 and spread metrics. The DOE size is 135 points (5 times 27 design variables).

Figure 11 shows the result of this experiment, only for the continuous relaxation and dummy coding Kriging kernels. It shows that the continuous relaxation approach, i.e. treating discrete variables as continuous variables during training and infill optimization, results in the fastest convergence to the known Pareto front. The spread metrics show large differences between the selection of the infill criterion, which is investigated in a later experiment. Other than the continuous relaxation and dummy coding approaches, all special-purpose mixed-integer and hierarchical Kriging kernels did not manage to meaningfully improve the Δ𝐻𝑉 and spread as obtained by the initial DOE. Implementation of these kernels, as also available in the published code, has been validated using relevant analytical test problems used in the development of these kernels (e.g. from [START_REF] Zaefferer | A First Analysis of Kernels for Kriging-Based Optimization in Hierarchical Search Spaces[END_REF][START_REF] Zuniga | Global optimization for mixed categorical-continuous variables based on Gaussian process models with a randomized categorical space exploration step[END_REF][START_REF] Pelamatti | Bayesian optimization of variable-size design space problems[END_REF]. Unfortunately, all these test problems are single-objective problems, and of a relatively low dimensionality. An investigation into infill behavior has been performed, of which some results are shown in Fig. 18 (in the appendix). It shows that mixed-integer and hierarchical Kriging kernels result in less correlation between distance to the known Pareto front and the objective value of the infill criterion. Similar results have been demonstrated in [START_REF] Karlsson | Continuous surrogate-based optimization algorithms are well-suited for expensive discrete problems[END_REF]: they conclude that continuous SBO algorithms are well capable of, and often better at, solving discrete optimization problems.

Finally, it should be noted that compared to continuous relaxation, mixed-integer and hierarchical kernels take up to 100 times longer to train the model, and up to 20 times longer to evaluate an infill point. It is therefore concluded that the continuous relaxation Kriging approach is to be selected for high-dimensional architecture optimization problems, as this results in superior performance from a prediction, correlation between infill objective and distance to Pareto front, and training and infill time point of view.

Infill Batch Size Investigation

The modulation of multi-objective (MO) infill criteria (see Fig. 9) offers potential improvements over directly using MO infill criteria due to the selection of multiple infill points per iteration. Selecting multiple infill points enables opportunities for a better-spread selection of infill points, requires less surrogate model training iterations, and allows time-savings in the expensive black-box evaluation if parallelization is possible. An experiment is performed to investigate the difference between one infill point per iteration and multiple infill points, and to give a suggestion on the amount of infill points per iteration. A comparison is made using the Expected Improvement (EI), Minimum Probability of Improvement (MPoI), and Minimum Variance of Pareto Front (MVPF) infill criteria with a Kriging surrogate model, and the RBF ŷ-Dist algorithm. EI is an infill criterion originally developed for use in single-objective optimization problems, and it naturally extended for use in multi-objective optimization by considering the EI criterion for each objective dimension separately (see Fig. 7). The MPoI infill criterion is an MO infill criterion, and therefore needs modulation if multiple infill points per iteration are needed (see Fig. 9). MVPF is an approach for selecting infill points from a population of candidate infill points after infill optimization, and naturally can be applied for selecting multiple infill points that are well-spread along the Pareto front. RBF ŷ-Dist attempts to trade-off exploration and exploitation by solving a multi-objective problem with exploration (finding the best expected value) and exploitation (finding the largest distance to existing design points) as objectives, and therefore also naturally lends itself to the selection of multiple infill points per iteration.

The main reason for using SBO is the reduction of the number of function evaluations to reduce optimization time. In this experiment it is assumed that perfect parallelization is possible, and that therefore the optimization time is only driven by the number of infill iterations, rather than the number of function evaluations. The continuous relaxation Kriging surrogate will be used, and otherwise all evaluation metrics are similar as for the experiments for the selection of the Kriging surrogate model.

The first result shows that if no parallel evaluation is possible, choosing 1 infill point per iteration converges to the Pareto front in the least amount of function evaluations needed, see Fig. 12. However, if the performance is compared by infill iterations, these differences disappear, and a larger number of infill points can lead to similar or faster convergence. An infill batch size of 5 is selected as the best size if multiple infill points are selected, due to the similar convergence speed (in iterations) as for 1 infill point, and a good final convergence to the Pareto front. For RBF ŷ-Dist, an infill batch size of 20 is selected, as there a larger infill size leads to a faster convergence (in terms of iterations), with less influence on the spread. The infill batch size has no influence on convergence if seen in terms of number of evaluations.

Figure 13 shows a comparison of the three infill criteria (EI, MPoI, and MVPF) for infill sizes of 1 and 5, and the RBF ŷ-Dist algorithm with infill sizes 1 and 20. It shows that in terms of Δ𝐻𝑉 and spread, the larger infill sizes converge faster to the known Pareto front, however result in a worse spread of points along the Pareto front. Therefore there exists a trade-off between exploration and exploitation with regard to the infill batch size. An additional issue noted during the experiment was that Kriging surrogate models can have difficulty training when there are too many training points, and therefore it might be better to keep the amount of training points as low as possible. This effect does not apply to the RBF surrogate model. To conclude, selecting multiple infill points per iteration may result in a closer approximation of the Pareto front (Δ𝐻𝑉) in less number of iterations, however selecting 1 infill point per iteration:

• Converges faster in terms of function evaluations;

• Results in a better diversification (spread);

• Reduces the amount of training points, which improves training time and reduces convergence problems for Kriging surrogate models.

C. Results and Conclusions

Results on the effectiveness of the considered optimization algorithms (see Table 4) for system architecture optimization will be presented. NSGA2 is compared with various SBO algorithms. The Kriging SBO algorithms use the continuous relaxation kernel. One new point is selected per infill iteration for SBO (it is assumed no parallel evaluation of design points is possible), and an evaluation budget of 400 evaluations in total is used. NSGA2 is run with an offspring size of 10 per iteration. The DOE (or initial population) size is 𝑛 𝑑𝑣 • 5 (where 𝑛 𝑑𝑣 is the number of design variables) individuals, which amounts to 135 points for the analytical benchmark problem.

Results of the algorithm comparison using the analytical benchmark problem are shown in Fig. 14. Best optimization algorithms in the four compared algorithm classes are shown: NSGA2, RBF ŷ-Dist, Kriging ŷ, and Kriging MVPF. More detailed views on all compared algorithms are shown in the appendix in Fig. 19. Results show that SBO converges faster than NSGA2. This is expected, as NSGA2 essentially is a randomized guided search algorithm [START_REF] Glover | [END_REF]; these kinds of algorithms are powerful at finding a global optimum or Pareto front, however need many function evaluations for this. Within the SBO algorithms, the Kriging SBO algorithms ( ŷ and MVPF) initially converge the fastest, however are later overtaken by RBF SBO. All four best algorithms achieve low spreads, indicating a good distribution of points along the Pareto front. Finally, it should be noted that training and querying RBF surrogate models is much faster than Kriging: the infill time varies between 1 and 2 ms for Kriging, whereas it stays around 0.2 ms for RBF; the training time varies between 10 and 20 seconds for Kriging and is about 50 ms for RBF.

V. Investigation of Optimization Algorithm Efficiency

Next, optimization algorithm efficiency is investigated by comparing the amount of architecture evaluations needed for convergence towards the benchmark problem Pareto-front. Efficiency can also be defined in terms of other metrics, like memory usage or wall time, however since it is expected that one function evaluation takes significantly more time and resources than the optimization algorithm steps, this is seen as the most relevant parameter.

A. Experimental Setup

Optimization algorithm efficiency will be tested by applying several multi-objective convergence criteria, shown in Tab. 5. The results of the previous effectiveness investigation will be used to investigate the impact of the convergence criteria on the different optimization runs. Convergence criteria are evaluated every 10 function evaluations, to allow for iterations where no better points are found. The Δ𝐻𝑉 and spread metrics are used to compare results. The efficiency experiment is applied to the best-performing algorithms from the effectiveness investigation: RBF ŷ-Dist, Kriging ŷ, and Kriging MVPF.

B. Results and Conclusions

Results on the efficiency experiment applied on the analytical benchmark problem are presented in Fig. 15. It shows that there is a clear trade-off between the number of function evaluations and the achieved quality of the multi-objective results in terms of distance to the Pareto front (Δ𝐻𝑉). The MDR (Mutual Domination Rate) and HV (Hypervolume) convergence metrics offer the largest reduction in function evaluations for the least reduction in Δ𝐻𝑉, most consistently. An example of the convergence zone in an optimization run is shown in Fig. 15d, which indeed seems to be shortly before the Δ𝐻𝑉 remains steady. The spread is not affected significantly by a reduction in function evaluations. This can be explained by the lack of a clear trend in the spread metric throughout an optimization run, as can be seen in Fig. 14b.

VI. Surrogate-Based Optimization Subject to Hidden Constraints

Hidden constraints, also known as unknown, unspecified, and forgotten constraints, exhibit themselves when a design evaluations could not complete correctly and yields a meaningless results [START_REF] Mãijller | Surrogate Optimization of Computationally Expensive Black-Box Problems with Hidden Constraints[END_REF]. This is a common occurrence in simulation-based optimization due to the non-convergence of numerical solvers, for example because the solution is not physically possible [START_REF] Digabel | A Taxonomy of Constraints in Simulation-Based Optimization[END_REF]. Up to 60% of evaluations in an optimization run may result in violated hidden constraints. A realistic engine architecting problem published in [START_REF] Bussemaker | System Architecture Optimization: An Open Source Multidisciplinary Turbofan Architecting Problem[END_REF] shows similar results: in the initial design of experiments of 205 design points, approximately 67% did not yield a meaningful result due to an unconverged design loop. The number of infeasible design points may be reduced by choosing a design space parameterization that reduces the chance for representing infeasible designs: for example, much effort has been invested in developing parametric geometry representations in the aircraft design discipline [START_REF] Kulfan | Fundamental" Parameteric Geometry Representations for Aircraft Component Shapes[END_REF]. Despite these efforts, hidden constraints will still play a role in some simulation-based optimization problems.

Evolutionary algorithms like NSGA-II are able to solve problems where large parts of the design space are subject to hidden constraints, however need many function evaluations to do so. For example, in [START_REF] Bussemaker | System Architecture Optimization: An Open Source Multidisciplinary Turbofan Architecting Problem[END_REF] a Pareto front of optimal engine architectures is found in 4000 function evaluations. Applying the SBO algorithms investigated in previous sections to the engine architecting problem shows that these algorithms are unable to produce a large improvement over the initial DOE. To support the claim that this is because of hidden constraints, the analytical benchmark problem used in previous sections is modified to also include hidden constraints. A hidden constraint is added that is violated if:

• 𝑔 2 is violated, or This results in a design space where in an initial design of experiments approximately 60% of design points yield an infeasible result, and a Pareto front with several gaps due to these hidden constraints, see Fig. 16. Running the same experiment as for determining algorithm effectiveness (i.e. 135 DOE points, 400 evaluations in total, 1 infill per iteration for SBO) shows that SBO algorithms are not significantly better than NSGA2 when exposed to hidden constraints, see Fig. 17. Kriging SBO algorithms are not able to improve the Δ𝐻𝑉 metric at all. The RBF SBO algorithms are able to find an improvement, although one that stabilizes at Δ𝐻𝑉 = 0.35, whereas without hidden constraints the Pareto front is approached much more closely. It is also shown that the number of infeasible infill evaluations is approximately 30% for NSGA2 and Kriging SBO algorithms, whereas it is up to 45% for RBF SBO algorithms.

At this point it can be concluded that the previously investigate SBO algorithms do not perform significantly better than NSGA2 for problems with hidden constraints. Hidden constraints need to be considered when developing and selecting algorithms for system architecture optimization.

VII. Conclusion and Outlook

System architecture optimization problems are hard to solve, because they in general suffer from a hierarchical design space, mixed-integer design variables, and multiple design objectives. A literature review discussing different types of optimization algorithms that are able to solve these kind of optimization problems is presented. Most suitable for this are Multi-Objective Evolutionary Algorithms (MOEA's). However, MOEA's suffer from needing many function evaluations for finding a Pareto-front. Surrogate Based Optimization (SBO) algorithms are specifically designed to reduce the number of function evaluations needed by constructing mathematical models of the design space. SBO algorithms use infill criteria to determine places in the design space to sample, with a trade-off between exploration and exploitation. In the past, research has been performed to enable SBO's to solve mixed-integer problems and multi-objective problems.

Optimization algorithms considered in the effectiveness (ability to find Pareto front) and efficiency (number of function evaluations needed to find Pareto front) experiments include NSGA2 (a MOEA) and various RBF-and Kriging-based SBO algorithms. RBF SBO does not expose information about estimation uncertainty, and therefore only directly sampling the surrogate ( ŷ) or applying a distance-based exploration penalty ( ŷ-Dist) can be used as infill. For Kriging SBO, many different single-objective and multi-objective infill criteria are investigated. Infill criteria are modified to be used in multi-objective optimization. The benchmark problem used in the experiments is based on a modification of a hierarchical Rosenbrock problem. All code implementing the various optimization algorithms, test problems, and experiments has been published at [START_REF] Bussemaker | jbussemaker/ArchitectureOptimizationExperiments: Release for AIAA paper[END_REF].

An investigation is performed into the best Kriging kernel for Kriging SBO. Next to continuous relaxation and dummy coding, several special-purpose mixed-integer and hierarchical Kriging kernels are investigated. Only the continuous relaxation and dummy coding kernels were able to approximate the known Pareto front, due to bad performance in modeling the high-dimensional hierarchical search space of the analytical benchmark problem. Continuous relaxation provided the best performance, both in optimization effectiveness and training and infill time. It should be investigated in more details why the special-purpose mixed-integer and hierarchical Kriging kernels were not effective for mixed-integer and hierarchical architecting problems. One area of focus should be the high-dimensionality of architecting problems in general. Additionally, it should be investigated whether machine learning models might be used as surrogate models for SBO [START_REF] Chugh | A survey on handling computationally expensive multiobjective optimization problems with evolutionary algorithms[END_REF], even though normally they might require a large number of samples for precise estimates [START_REF] Goodfellow | Deep Learning[END_REF].

Next, it is investigated how many infill points should be generated at every infill iteration. It is concluded that selecting multiple points per infill may result in a closer approximation of the Pareto front, however selecting 1 infill point per iteration: converges faster in terms of function evaluations, results in a better diversification along the Pareto front, and reduces the amount of training points in surrogate model.

Algorithm effectiveness is investigated by comparing NSGA2 with RBF and Kriging SBO. In case of SBO, 1 infill point per iteration is used, and for Kriging SBO the continuous relaxation kernel is used. It is shown that for SBO performs significantly better than NSGA2 in terms of function evaluations needed to find the Pareto front. For RBF SBO, the ŷ-Dist infill provides the best performance; for Kriging ŷ (direct evaluation) and MVPF (Minimum Variance of the Pareto Front) provide best performance. Algorithm efficiency is investigated by applying various convergence criteria on the effectiveness results. It is shown that the MDR (Mutual Domination Rate) and HV (Hypervolume) convergence metrics are best able to reduce the number of function evaluations without reducing distance to the Pareto front too much.

Even though it is shown that SBO algorithms perform well on the mixed-discrete, hierarchical, multi-objective analytical test problem, they do not offer any performance improvement over NSGA2 when hidden constraints are added. Future research will focus on hidden constraints in architecture optimization, to enable the application on realistic simulation-based optimization problems that may be subject to analysis convergence issues. The engine benchmark problem from [START_REF] Bussemaker | System Architecture Optimization: An Open Source Multidisciplinary Turbofan Architecting Problem[END_REF] will be used to test specific approaches for dealing with hidden constraints in surrogate-based optimization. Surrogate-based algorithms that can deal with mixed-discrete, hierarchical and multi-objective design space will be applied in various industry-provided use-cases in the AGILE 4.0 project. These use-cases feature evaluation times of multiple minutes to hours per design point, and will therefore serve as a good example of why surrogate-based optimization that can deal with architecture design spaces is needed. This work considered no multi-level optimization approach, where for example continuous variables are optimized at an inner-loop compared to an outer-loop optimizing the discrete design variables [START_REF] Barjhoux | A bi-level methodology for solving large-scale mixed categorical structural optimization[END_REF]. Scalarization of the multi-objective problem was also not considered (e.g. [START_REF] Désidéri | Multiple-Gradient Descent Algorithm (MGDA)[END_REF]), and although often weights involved in scalarization lead to an element of arbitrariness, it might speed up convergence for certain classes of optimization problems.

Figure 1

 1 Figure 1 The goals of multi-objective optimization: progress towards the Pareto-front while maintaining sufficient diversity, adapted from [8].

Figure 2

 2 Figure 2 Three strategies for dealing with hierarchical optimization problems [10]: naive (ignore the effects), imputation, and explicit consideration. The corrector converts a design vector to its canonical form (𝑥 𝑖𝑚 𝑝 ) by imputing inactive design variables, and optionally returns which design variables are active (𝛿). This makes sure design points are only evaluated once. The contour plots on the right give a notional view of how accurate the optimizer can model the design space, based on the test problem from [10]: on the left side of the domain, only one design variable is active.

Figure 3

 3 Figure 3 General procedure of an Evolutionary Algorithm, adapted from [8].

Figure 4

 4 Figure 4 Pareto dominance relation, from [16]. In this example, point D and E are both in the Pareto set, since they dominate all other points, however not each other.

Figure 5

 5 Figure 5 Process of the Efficient Global Optimization (EGO) algorithm, adapted from [21]. Surrogate-Based Optimization (SBO) algorithms are based on the same process, with variations in the initial sampling strategy, the underlying surrogate model, the infill criterion, and the stopping criterion.

Figure 6

 6 Figure 6Illustration of the Probability of Improvement (PoI), from[START_REF] Hawe | An Enhanced Probability of Improvement Utility Function for Locating Pareto Optimal Solutions[END_REF]: the probability that the function value will be better than some value 𝑇 < 𝑓 𝑚𝑖𝑛 given the surrogate model estimate ŷ and variance ŝ.

  (a) First infill criterion objective. (b) Second infill criterion objective.

Figure 7

 7 Figure 7 Modification of the Probability of Improvement (PoI) for use in multi-objective optimization, visualized for a normalized bi-objective solution space with 5 hypothetical points on the Pareto front, and a constant variance. The shown criterion objectives are to be minimized for finding infill points.

Figure 8

 8 Figure 8 Visualization of multi-objective infill criteria for a hypothetical normalized bi-objective solution space ( 𝑓 0 and 𝑓 1 ) with 5 points on the Pareto front, and a constant variance.

( a )

 a Infill criterion objective. (b) Modulated over 𝑓 0 . (c) Modulated over 𝑓 1 .

Figure 9

 9 Figure 9 Modulation of a multi-objective infill criterion over the two underlying function objectives. This guides the selection of infill points along the currently existing Pareto front. Shown is the Enhanced Probability of Improvement infill criterion [33] for a hypothetical normalized bi-objective solution space ( 𝑓 0 and 𝑓 1 ) with 5 points on the Pareto front, and a constant variance.

  (a) Addition of a second objective 𝑓 2 . (b) Repetition of the problem.

Figure 10

 10 Figure 10 Modifications of the hierarchical Rosenbrock problem [53]. "SP" refers to the sub-problem index (see [53]), "rep." to the replication index.

( a )

 a Δ𝐻 𝑉 metric (lower is better). (b) Spread metric (lower is better).

Figure 11

 11 Figure 11 Comparison of Kriging with continuous relaxation and dummy coding kernels on the analytical benchmark problem. Both Kriging models are compared using the EI (Expected Improvement) and MPoI (Minimum Probability of Improvement) infill criteria. Continuous relaxation provides faster convergence towards the Pareto front (Δ𝐻𝑉). However, the MPoI infill criterion with continuous relaxation results in a worse spread.

  (a) Comparison by function evaluations. (b) Comparison by infill iterations.

Figure 12

 12 Figure 12 Comparison of Kriging with continuous relaxation and MPoI (Minimum Probability of Improvement) infill criterion for different infill batch sizes. For an infill batch larger than 1, the modulated infill criterion is used (see Fig. 9). In terms of function evaluations, choosing only 1 infill point per iteration results in fastest convergence. If parallel evaluation is possible, infill sizes larger than 1 can result in faster convergence.

  (a) Δ𝐻 𝑉 metric.(b) metric.

Figure 13

 13 Figure 13 Comparison between 1 or 5 design points per infill iteration for different Kriging infill criteria, and the RBF ŷ-Dist algorithm. Larger infill batch sizes result in faster convergence (Δ𝐻𝑉) in terms of number of iterations, however with a worse spread along the Pareto front.

  (a) Δ𝐻 𝑉 metric. (b) Spread metric.

Figure 14

 14 Figure 14 Comparison between algorithm effectiveness on the analytical benchmark problem. Shown are the best algorithm in the four compared classes: MOEA, RBF SBO, Kriging SBO with SO-infill criteria, and Kriging SBO with MO-infill criteria. SBO converges faster than MOEA (NSGA2). Within the SBO algorithms, Kriging SBO with MO-infill criteria converges the fastest. Spread values of all algorithms are comparable, however RBF SBO and Kriging SBO with MO-infill result in the highest spreads.

  (a) Reduction in function evaluations. (b) Resulting Δ𝐻 𝑉 . (c) Resulting spread. (d) HV convergence applied to the RBF ŷ-Dist algorithm.

Figure 15 Figure 16

 1516 Figure 15Comparison between algorithm efficiency on the analytical benchmark problem. Shown are average resulting number of evaluations, Δ𝐻𝑉, and spread for different algorithms and convergence criteria. There is a clear trade-off between remaining Δ𝐻𝑉 and number of evaluations; the effect on the resulting spread is not significant.

Figure 17

 17 Figure 17 Comparison between optimization algorithms on the analytical benchmark problem with hidden constraints. It shows that SBO algorithms are not able to find a significant improvement towards the Pareto front, if they are at all able to find an improvement.

  (a) RBF SBO (Δ𝐻 𝑉 ). (b) RBF SBO (spread).(c) Kriging SBO with SO-infill (Δ𝐻 𝑉 ).(d) Kriging SBO with SO-infill (spread).(e) Kriging SBO with MO-infill (Δ𝐻 𝑉 ).(f) Kriging SBO with MO-infill (spread).

Figure 19

 19 Figure 19 Comparison between algorithm effectiveness on the analytical benchmark problem. Focus on different SBO algorithm classes. The Kriging SBO algorithm are all based on continuous relaxation.

Table 1 Approaches for applying surrogate models on mixed-integer (MI) and hierarchical (H) design spaces. 𝑛 𝑐 represents the number of combinations of the discrete variables. Note that the MI kernels imply that continuous kernels are used for continuous design variables.

 1

	Nr of Models Input space Input pre-processing	Kriging Kernel	Reference
	𝑛 𝑐	Continuous	N/A	Continuous	e.g. Pelamatti et al. [45]
			Continuous	Continuous	e.g. Garrido-Merchán et al. [30]
			relaxation
			Dummy coding	Continuous	Herrera et al. [46]
		MI		Hamming distance (Ham)	e.g. Zaefferer et al. [47]
				Gower distance (Gow)	Halstrup [48]
				Symbolic Covariance (SC)	McCane et al. [49]
	1			Compound Symmetry (CS)	Roustant et al. [50]
				Latent Variables (LV)	Zhang et al. [51]
			N/A	Arc	Hutter et al. [11]
				Indefinite Conditional (Ico)	Zaefferer et al. [10]
		MI + H		Imputation (Imp)	Zaefferer et al. [10]
				Wedge (Wed)	Horn et al. [52]
				SPW decomposition	Pelamatti et al. [53]
				DVW decomposition	Pelamatti et al. [53]

Table 2 Overview of Surrogate-Based Optimization (SBO) approaches, with a focus on Mixed-Integer (MI), Multi-Objective (MO), Hierarchical (H) and constraint handling (𝑔) approaches. Abbreviations: Expected Improvement (EI), Expected Violation (EV), Coordinate Perturbation (CP), Kriging with Partial Least Squares (KPLS), Probability of Feasibiliy (PoF), Probability of Improvement (PoI), Expected Hyper-Volume Improve- ment (EHVI), Euclidean EI (EEI), Expected Maximin Fitness Improvement (EMFI), Minimum Variance of Pareto Front (MVPF), Lower Confidence Bound (LCB), Upper Trust Bound (UTB). ŷ and ŝ refer to the function and variance estimates, respectively.

 2

	Publication	Year	Name	Surrogate	Infill criteria	MI H MO	𝑔
	Jones et al. [23]	1998	EGO	Kriging	EI
	Gutmann [34]	2001		RBF	ŷ
	Sasena [21]	2002	SuperEGO	Kriging	EI			EV
	Gramacy et al. [35] Regis et al. [38]	2004 2013	DYCORS	Kriging RBF	ŝ CP			!
	Yi et al. [57]	2014		Kriging	ŷ vs ŝ			ĝ
	Bartoli et al. [58] Hemker [60] Chen et al. [61]	2016 2008 SurOpt MINLP SEGOMOE 2015	KPLS, MOE Kriging Kriging	ŷ vs EI (WB2S [59]) ŝ ŷ vs ŝ	! !		! ! !
	Baert et al. [55] Müller [29] Roy et al. [28]	2015 2016 2017	MV-SBO MISO AMIEGO	Kriging w/ MI kernels RBF, Kriging KPLS	CP, EI, ŷ EI	! ! !		! EV
	Zaefferer et al. [10]	2018		Kriging w/ MI kernels	EI	! !
	Horn et al. [52]	2019		Kriging w/ MI kernels	EI	! !
	Pelamatti et al. [45] Priem et al. [39] Jeong et al. [62] Knowles [63] Emmerich et al. [64] 2006 2019 2020 2005 2006 Keane [43] 2006 Hawe et al. [33] 2007 Svenson et al. [65] 2016 Datta et al. [26] 2016 Palar et al. [66] 2017 Rahat et al. [67] 2017 Mueller [68] 2017 dos Passos et al. [44] 2018 Tian et al. [24] 2019	Multi-EGO ParEGO SMES-RBF SOCEMO	Kriging w/ MI kernels Kriging Kriging Kriging Kriging Kriging Kriging Kriging RBF Kriging Kriging RBF Kriging Kriging	EI WB2S MO EI Scalarized EI MO EI, PoI EEI Enhanced PoI EMFI ŷ EHVI, EEI Minimum PoI ŷ EHVI, MVPF ŷ vs ŝ	! ! !	PoF UTB ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! !
	Tran et al. [69]	2020 srMO-BO-3GP	Multi-level Kriging	Scalarized EI		! !

Table 3 Comparison of implemented multi-objective infill criteria. 𝑛 𝑝 𝑓 represents the number of points in the existing Pareto front.

 3

	Name	Publication	Year	Computation Cost
				Relative 𝑛 𝑝 𝑓 dependence
	Expected Hypervolume Improvement (EHVI)	Emmerich et al. [64] 2006	~500	N/A
	Euclidean Expected Improvement (EEI)	Keane [43]	2006	~10	𝑛 2 𝑝 𝑓
	Enhanced Probability of Improvement (EPoI) Expected Maximin Fitness Improvement (EMFI)	Hawe et al. [33] Svenson et al. [65]	2007 2016	~10 ~10	𝑛 2 𝑝 𝑓 √ 𝑛 𝑝 𝑓
	Minimum Probability of Improvement (MPoI)	Rahat et al. [67]	2017	1	N/A
	Minimum Variance of Pareto Front (MVPF)	dos Passos et al. [44] 2018	~0	N/A

(a) Expected Hypervolume Improvement (EHVI)

[START_REF] Emmerich | Single-and multiobjective evolutionary optimization assisted by Gaussian random field metamodels[END_REF]

Table 4 Compared optimization algorithms.

 4

	Class	Name	Surrogate Model Infill Criterion Constraints
	MOEA	NSGA2	N/A	N/A	Reject
	RBF SBO	RBF ŷ	RBF	ŷ	Reject ĝ
		RBF ŷ-Dist		ŷ-Dist
		Krg ŷ		ŷ
	Kriging	Krg ŝ		ŝ
	SBO with	Krg LCB		LCB
	SO infill	Krg EI		EI
		Krg PoI		PoI
		Krg EEI	Kriging	EEI	PoF
	Kriging SBO with MO infill	Krg EPoI Krg MPoI Krg MEPoI Krg EMFI		EPoI MPoI MEPoI EMFI
		Krg MVPF		MVPF

Table 5 Investigated convergence criteria. Δ in the convergence limit indicates that the rate of change of the (filtered) convergence criterion is checked, rather than the absolute value. EMA represents the Exponential Moving Average filter.

 5

	Name	Convergence criterion	Filter	Settings	Convergence limit
	Budget	Evaluation budget			> 400
	MDR	Mutual Domination Rate	EMA(𝑛 = 2)		< 10%
	MGBM Martí-García-Berlanga-Molina Kalman(𝑟 = 0.1, 𝑞 = 0.1)		< 10%
	FHI	Fitness Homogeneity Indicator	EMA(𝑛 = 2)		Δ < 1e-4
	GD	Generational Distance	EMA(𝑛 = 2)		Δ < 1e-3
	IGD	Inverse Generational Distance	EMA(𝑛 = 2)		Δ < 1e-3
	Spread	Spread	EMA(𝑛 = 5)		Δ < 1e-4
	HV	Hypervolume	EMA(𝑛 = 2)		Δ < 1e-3
	CR	Consolidation Ratio	EMA(𝑛 = 2)	𝑛 𝑑𝑒𝑙𝑡 𝑎 = 1	> 80%
	MCD	Maximal Crowding Distance	EMA(𝑛 = 2)		Δ < 5e-4
	SPI	Steady Performance Indicator	EMA(𝑛 = 2)	𝑛 𝑙𝑎𝑠𝑡 = 4	< 2%

  • 𝑓 1 is within 3 of a multiple of 20 (i.e. |0.5 -( 𝑓 1 /20) mod 1| > 0.35), or • 𝑓 2 is larger than 1000 and within 15 of a multiple of 100 (i.e. |0.5 -( 𝑓 2 /100) mod 1| > 0.35).

* pymoo: https://pymoo.org/ † scikit-learn: https://scikit-learn.org/ ‡ SMT (Surrogate Modeling Toolbox): https://smt.readthedocs.io/ § jbussemaker/ArchitectureOptimizationExperiments

Acknowledgments

The research presented in this paper has been performed in the framework of the AGILE 4.0 project (Towards Cyber-physical Collaborative Aircraft Development) and has received funding from the European Union Horizon 2020 Programme under grant agreement n• 815122.

Appendix: Extended Results