Nowadays, solving of many engineering problems, such as the development of transport vehicles, involves the use of computationally expensive simulations, such as computational fluid dynamics (CFD), the finite element method (FEM), or others, as well as their combinations, e.g. the fluid-structure interaction (FSI). These will be further addressed as simulations. In practice, typically, such simulations suffer from two major drawbacks. Firstly, even thought the computational power available to an engineer for running such simulations has been gradually increasing over years, in practice the requirements on the model accuracy had increased as well [48], and as such, there is still a large quantity of simulations that take hours, days and occasionally even more, just to evaluate a single design [43]. Secondly, they are often of a black-box nature, meaning that given an input X, an output Y is returned to the user without an accompanying explicit relation between the two [41].
Within a practically infinite space of design choices for the system under study, the task of an engineer is to identify the most influential design parameters and understand their degree of influence. Through the use of simulations, often the next goal is to meet a set of threshold requirements on selected targets, while minimizing or maximizing a certain property of the overall system, such as cost, weight, etc. Mathematically, this can be formulated as a multi-objective optimization problem with the goal to: minimize f m ( X) m = 1, ..., P
subject to g j ( X) ≤ 0 j = 1, ..., Q
where f are the objective functions, g and h are the inequality and equality constraints, respectively, and X is the vector of design variables in the design search space bounded by X L i and X U i from the lower and upper side, respectively [47]. In case of multi-objective optimization when m ≥ 2, the notion of singular optimality needs to be augmented as there is typically not a single solution X * minimizing all of the objectives of Equation 1 at the same time. In such cases, the concept of so-called Pareto optimality is considered [36], which can be loosely described as: "for each solution that is contained in the Pareto set, one can only improve one objective by accepting a tradeoff in at least one other objective" [39].
Due to said high computational costs and the black-box nature of simulations, multi-objective parametric studies and optimizations often become too lengthy and impractical for actual application [42,41]. The lack of derivatives required for derivative-based optimization methods, such as in Liu and Reynolds [35] or Peitz and Dellnitz [39], can be addressed by the use of derivative-free optimization (DFO) methods, such as genetic, evolutionary or swarm algorithms [27]. However, these methods compensate for the lack of gradients by evaluating on larger sets of candidate solutions, which leaves the issue of high computational cost unresolved. One way to resolve that is by using a model of lower fidelity requiring less computational resources [33]. However, the accuracy of such lower fidelity models can become unacceptably low. As an alternative, a surrogate-based approach can be adopted [20,5,12].
A surrogate model, also addressed as a metamodel, is a mathematical model that upon training approximates the response of the original model based on a sample set from the design space [31]. Such computationally cheap surrogate can then be called during the optimization defined by Equation 1-4 instead of the original, expensive simulation. In addition, the use of a surrogate model opens up several new opportunities, which include:
1. using the surrogate as a substructure of a larger model [34] 2. cheap numerical evaluation of gradients using the surrogate 3. repeating the optimization for various formulations of the problem (within the bounds of validity of the surrogate) -compare various formulations of the same car side impact problem: Youn et al. [54], Guo, Wang, and Wang [26], and Tanabe and Ishibuchi [46] 4. design exploration through sensitivity studies and visualization [49,38] On the other hand, a disadvantage is the added complexity of the overall task due to introduction of extra tuning parameters related to the surrogate model and its training. Due to the often black-box nature of simulations, it is not possible to a-priori quantify the sample set size required for an accurate training of the surrogate [4]. Therefore, for a purely optimization-focused approach, in general it is not guaranteed that such surrogate-based approach is indeed computationally cheaper than optimizing directly using the original model. Nevertheless, it has been established that for budget-limited tasks with expensive simulations, it is indeed the case [44,24,38]. This is the case especially for multi-objective optimization, when the exploration of the full Pareto front is required. Confirming results of an empirical comparison can be found e.g. in Voutchkov and Keane [50] on benchmark problems and in Bessa and Pellegrino [3] on a practical problem of the design of ultra-thin carbon fiber deployable shells.
Example uses of surrogate-based optimization (SBO) from the area of structural optimization include the design of 3D weaving composite stiffened panels by Fu, Ricci, and Bisagni [21], hypersonic vehicle's metallic thermal protection by Guo et al. [25] or a wellhead connector for a subsea Christmas tree by Zeng et al. [55]. The SBO methodology has been already integrated into various commercial software packages, such as OptiSLang 1 1
OptiSLang, https://www.dynardo.de/en/ software/optislang/ansys-optislang.html [Accessed on 03/08/2020] or HEEDS 2 . This can be often practical, but not always. Firstly, the user must rely only on the available capabilities of such software packages, which can sometimes be insufficient for the application at hand. Secondly, especially for individuals or smaller businesses, the license price for such packages can be too expensive, thus inaccessible. Last but not least, for research purposes, commercial solutions often offer only a limited insight into the used methods and source code, leading to limited customabilizity, a feature often required for research.
In this paper, the core of a developed Python opensource SBO framework is presented. In Section 2, the different elements of the proposed framework are presented. Next, in Section 3, the use of the framework is demonstrated on a benchmark problem. Finally, the capabilities of the framework are summarized in Section 4 together with an outline of suggested further developments.
This framework integrates elements of surrogate-based modeling and automated model selection together with derivative-free optimization. One of the driving ideas behind it is modularity, such that it could be readily customized for the user's particular needs. Upon completion of the master thesis project, the framework will be made publicly available at github.com/apanzo/ optimization.
The top-level illustration of the framework is shown in Figure 1. As a starting point, the simulation is prepared, and the framework settings, such as the selection of a surrogate model and optimization algorithm, are set. Both the use of a surrogate and performing optimization are optional. The surrogate loop contains sample selection, model evaluation, results storage and retrieval, optimization and training of the surrogate, and finally accessing its convergence. The optimization loop consists of the actual optimization, and, in case that a surrogate model is used, a verification of the obtained results against the original model. In the following subsections, key parts of the framework are discussed in the order from optimization start to end as indicated in Figure 1.
Considering the expensiveness of the simulations addressed by this framework, one of the keys of success is the smart selection of the sample points within the explored design space, such that the surrogate model obtains a sufficiently representative sample to be trained on. For this purpose, simpler sampling strategies such as the grid search or random search are not suitable. The former is expensive as it does not benefit from sampling in multiple dimensions, as the sample projection on each of the input's axes is independent of the number of inputs dimensions, and in the case of re-sampling with a finer sample, either the new sample discards all of the previous samples, or the re-sampling is to be done on an integer multiple of samples per input dimension of the current sample. The random search does not suffer from those, but its weakness lies in its poor space-filling property, thus a large amount of samples is required to generate a sufficiently representative sample.
As a baseline within this framework, quasi-random sampling strategies are implemented. In particular, the two available sampling methods are the cell centered Latin hypercube sampling (LHS) [14,18] and the original Halton sequence sampling [13]. The principal logic behind the LHS method is that it is a sparse grid method where the projection of the sample on each of the input's axes contains only one sample per uniform interval. As such, the sample size is significantly reduced compared to the full grid. On the contrary, the disadvantage of the poor resampling pertains. This is not the case for the Halton sequence, which is a multi-dimensional version of the van Der Cor-put sequence. It is defined as P = {φ q 1 (Z), ..., φ q n (Z)} with bases q 1 , ..., q n that are mutually coprime and in practice taken as the first s primes [13,22]. φ q (Z), the inverse radix number, is obtained from an integer
where q is the basis and the largest exponent is r = int ln Z ln q , and 0 ≤ a i < q are the coefficients [13,11], as
The advantage is that each point of the sequence is generated independently of the previous point. An example is shown in Figure 2. It is noted however that the original method is suitable only up to 8 input dimensions, as in higher dimensions, spurious correlations between the inputs occur, considerably deteriorating the sample's quality [13]. This can be leveraged by using its modified versions that address this issue [19,37]. To take it one step further, considering that these, so called static sampling methods, only take into account the space-filling criteria, that is the input's quality, and not the nature of the response, that is the output's quality, an adaptive sampling method that does so is implemented as well. As an example, if S is the full design space and the response of a model is flat everywhere apart from the region D ∈ S, where there is a valley in the response, only a few samples are required outside D, while in D more samples are required to train the surrogate. The selected method is from Eason and Cremaschi [17], which determines the new samples based on a balance of the variance in predictions of the K models from cross-validation, explained further in subsection 2.4, for exploitation and the distance to the nearest sample for exploration, each normalized by its maximum. The new sample is chosen where the sum of the two criteria is the largest.
For the practical use of the framework, it has to be able to interact with external solvers that evaluate the simulations. This independence is one of the main strong points of the framework, as it allows to integrate results from different solvers. As an example, for a FSI study, different FEM and CFD solvers can be used and their results integrated within this framework.
The implementation of different evaluators is application dependent, therefore the interaction with each new solver has to be customized according to its specific interface. In the course of development of this framework, a custom ANSYS evaluator has been integrated. For its use, the procedure is to firstly build a parametric model inside the ANSYS Workbench environment 3 , with the input and output parameters defined. These parameters are passed over by name to the framework together with the path to the prepared model. The rest of the interaction is automated. After determining the sample points, firstly the framework checks that there are available free licenses for carrying out the computations. If there are available licenses, a request is sent to ANSYS Workbench to evaluate the selected sample points. Once all the samples have been evaluated, the framework retrieves the defined output parameters and integrates them within its internal database.
Among others, the most commonly used surrogatemodeling techniques are kriging, radial basis functions (RBF), artificial neural networks (ANN) or support vector machines (SVM) [17,12]. The former two are integrated within the framework using the Surrogate Modeling Toolbox (SMT) library [10] natively, while the ANN is integrated as a custom class into SMT using the TensorFlow library [1]. In the course of the reported research, the ANN has been used due to previous experience in the research group, therefore it is elaborated upon closer in this section.
In the general sense, ANNs are one of the subclasses of the broader field of machine learning (ML), which in turn is a subclass of artificial intelligence (AI). Most commonly, the ANNs are trained using a form of supervised learning, that is by providing both the input and output data during the training process. They can be used both for classification, that is categorization of the input content into different sub-classes with common features, as well as regression. For the purpose of optimization, regression is typically used and the task of learning is to obtain the unknown relation between input and output data.
The core concept of an ANN is a single neuron. Inspired by the function of a real neuron inside the living brain, the mathematical model from input x to output y consists of performing a weighted summation of all the inputs
and applying a non-linear transformation
with g the so-called activation function. Common activation functions are the logistic sigmoid
To provide a flexible range of outputs based on the neuron's activations, the neurons are organized into layers. The typical neural network architecture for regression is the multi-layer perceptron (MLP), consisting of the first and last layers containing the amount of neurons equal to the number of input and output dimensions, respectively, and a selected number of intermediate, hidden layers, which can each contain an arbitrary amount of neurons. Building upon the Kolmogorov's theorem [32], Hecht-Nielsen, Ne, and Corporation [28] proved that using specific activation functions, even an ANN with a single hidden layer of neurons can approximate any continuous function. However, it is not guaranteed that such ANN can actually learn such representation [6], and in practice, using multiple hidden layers of neurons can simplify the learning process [29].
The training of the network is equivalent to the determination of appropriate weights w from Equation 5for each neuron, such that the desired mapping x → y from the input data to the output data is obtained. The common training method is the backpropagation method, where the prediction error E after the forward pass is propagated back to the contribution to that error from each neuron, and their weights are updated in each training iteration. This is schematically illustrated in Figure 3.
Independent of the selected surrogate model, to obtain an accurate model, a step of crucial importance is the selection of the model's hyperparameters, that is the parameters that are not directly learned through the training process. This is specific for each of the surrogate models. As an example, in the case of the ANN, these are Figure 3: Illustration of backpropragation the amount of hidden layers and neurons in each of them, the learning rate parameter, regularization parameters or the activation function choice. For this purpose, modelspecific methods, such as the neuron pruning [56], where the neurons of low significance are removed, as well as model-independent, such as a random search, Bayesian optimization or meta-heuristic optimization approaches, can be used.
At this stage, a random search and Bayesian optimization (BO) [45] are included for the ANN using Keras Tuner 4 . A hypermodel is defined with default hyperparameters stored in a configuration file. Within, these hyperparameters can be either changed manually, or specified to be subject of automatic tuning. For ANNs, the tunable hyperparameters are: number of hidden layers, number of neurons in each layer, initial learning rate, activation function and the regularization parameter. The current implementation of the BO allows only setting a fixed optimization budget, therefore the selected option is to start with an initial random sample of 3x times the number of optimized hyperparameters, and stop at 10x times the number. The acquisition function is the upper-confidence bound [45]. The frequency of hyperparameter optimization is set as at the first iteration, and then every time the sample set size increases by 50%.
A metric of the surrogate's accuracy of choice is to be tracked to determine the convergence of the sample set size. In our approach, the tracked metric is the mean absolute error (MAE)
where y pred is the prediction of the surrogate, y train is the training value and N is the total number of samples. Since the amount of data is relatively low, to perform validation of the trained surrogate model, the K-fold cross-validation approach is adopted [9]. In each iteration, the data is split into a set of K-folds, and K surrogates are trained each time leaving out a different fold of the data. This holdout set is used to calculate the generalization error of the model on previously unseen data. With increasing size of the data set, leaving out the holdout set will affect the surrogate's prediction less and less, so the average E MAE over the K surrogates will decrease. The resampling loop is terminated when a satisfactory level of E MAE is attained.
Once the surrogate model has been trained to a desired level of accuracy, the optimization run can be started. Within this framework, various population-based optimization algorithms are provided by the Pymoo library [7]. Pymoo provides implementations of a series of population-based algorithms, including:
• NSGA-II and NSGA-III
• RNSGA-II and RNSGA-III
Among those, the most commonly used in the industry is the NSGA-II method [15] and its updated NSGA-III version [16] adapted for the many-objective cases of m > 2, where an even exploration of the Pareto-front is accomplished by the use of reference directions. A generic scheme of a genetic algorithm is shown in Figure 4. The idea behind is similar to Darwinian evolution, that is that once the initial population of solutions is generated (Step 1), the fitness of the individuals is evaluated (Step 2) and only the fittest survive in each generation (Step 3). From those, a selection parent selection takes places (Step 4), from which new offsprings are generated through crossover (Step 5). Additionally, a random mutation is applied to their genome (Step 6) to maintain diversity within the population and explore the design space better. At the end of each iteration, the convergence criterion is evaluated, and either the optimization stops, or continues with a new generation [52].
The available termination criteria are the maximum:
• number of evaluations
• time
• design space tolerance
• objective space tolerance.
Due to the computational cheapness of the surrogate, the optimization process is not limited by the computational cost of evaluating the simulation, such as the maximal number of its evaluations. Therefore, the selection of the termination criterion can focus purely on the convergence of the optimization. In this case the last two criteria are of main importance. A threshold tolerance on the metrics tracking the evolution of solutions in the design or objective space is defined, and if the solutions in the specified amount of last iterations do not improve over this threshold, the optimization is terminated. For robustness purposes, also the maximal amount of evaluations or generations can be specified, e.g. for cases when the optimal solutions oscillate around some value, and thus in the defined window would otherwise never have converged.
It is also possible to enter the optimization directly using the original model without constructing the surrogate. In this case however, for an expensive simulation, the choices on the optimization algorithm and its components, especially the population size and number of generated offsprings in each generation, have to be more carefully selected with respect to the number of required simulations.
Once an optimum (set) has been found, given that this optimum has been found based on the surrogate, it is important to verify it's accuracy. In particular, it could be that the algorithm has found an optimum that lies in a spurious extremum caused by the surrogate's inaccuracy which does not pertain to the original model. Therefore, by submitting a sample from the Pareto optimal set back to an evaluation using the original model, it can be verified whether the obtained solutions indeed pertain to the original model. If the mean or maximal error within this verification exceeds a pre-defined tolerance, the surrogate is retrained on an enlarged sample, and the optimization occurs again, as seen in Figure 1.
As an important implementation note, the data is normalized within the internal data flow. This is done such that common settings can be used across different problems, as well as for a single problem to treat inputs and outputs of different magnitudes comparably. Not doing so would mean that for example when the prediction error is summed, the weight in kilograms would completely dominate deformation in millimeters, which is undesirable. Therefore, all the inputs and outputs are normalized by dividing each quantity with its absolute maximal value. In such a way, all the data is guaranteed to lie within the [-1,1] range, while maintaining its sign, conforming to the constraint violation formulation in Equation 2.
As an illustrative example of a practical problem, the results on a parametrized nonlinear FEM simulation of a car side impact problem [30] ran using the proposed framework are presented in this section.
The design space is defined by 7 input parameters
that represent the thicknesses of structural members such as B-pillars, beams or the floor. A single solution was, at the time of the problem's publication, reported to take about 20h [54]. Even thought this time has likely reduced with today's hardware, it is still a good illustrative example of surrogate-based modeling.
For the purpose of benchmarking, the response has been parametrized by analytical expressions, which are presented hereafter [30]. The 3 objectives of the design study are to minimize the structural weight, impact force on the passenger and average velocity of the side members that absorbs the impact load, represented as f 1 ( x) ≡ W = 1.98 + 4.9x 1 + 6.67x 2 + 6.98x 3 + 4.01x 4 +1.78x 5 + 0.00001x 6 + 2.73x 7
respectively. In addition, there are 10 constraints that take into account criteria such as the abdomen load
The format of the user input using the JSON formatting 5 is shown below:
" evaluator ": " benchmark " ,
" sampling ": " halton " ,
" d e fa u lt _ s am p le _ co e f ": 10 ,
" resampling ": " geometric " , To translate it, the carside problem is evaluated from the benchmark problems set. The initial sample is 10x the number of input dimensions, thus 70, and in every additional iteration, the sample set is increased by 20%. The convergence metric for the model is the MAE metric, with a maximal value of 0.03. The ANN is the used surrogate model and its performance metric is calculated using 5-fold cross validation. The hyperparameters are re-optimized when the amount of samples increases by 50%. The optimization uses the NSGA-III algorithm and its termination is based upon the objective space tolerance, which terminates, if there is no change by more than 0.001 in the tracked metrics: changes of the objective functions, inverted generational distance (IGD) [2] and the nadir point, over the last 5 generations. This convergence check is performed every 5 generations. Additionally, a maximum of 150 generations is allowed during a single optimization run. The verification mean error limit is set as 5%. Additionally, not present in the input file, the ANN's hyperparameters are defined as in Table 1. The output layer function is linear. The number of hidden layers and neurons per hidden layer are optimized using BO between [2,10] and [6,200], respectively. For optimization, the algorithm is initiated with default settings 7 . The reference directions are obtained using an energy method from Blank et al. [8] and their number, which also defines the population size, is set at 30x the number of objectives, that is 90 in this case. The amount of generated offsprings is equal to the population size.
In Figure 5, the comparison of training data with the prediction of the final trained ANN is shown. The sample set converged after 3 surrogate training iterations, that is with 101 samples evaluated. This is just over the defined population size in one generation of the genetic algorithm. The optimization took 70 generations to converge, which amount to 9000 evaluations of the model. That as a factor of 90x more than the required surrogate training sample size, confirming the advantage of using the surrogate. It is noted however, that thus far the selected optimization parameters are determined empirically on a trial-and-error basis.
As a validation of the result, the numerical values of the Pareto optimal solutions are unfortunately not presented in Jain and Deb [30] nor other paper, however, the qualitative comparison with their solutions in the objective space as shown in Figure 6 indicates that the ranges of all 3 objectives fall within nearly identical ranges. The final mean verification error between the Pareto optimal set and the original model, as discussed in subsection 2.6, was 2.10%.
In this paper, a working core example of a modular surrogate-based optimization framework was presented. A more detailed elaboration on the ANN surrogate was drawn, together with a top-level description of each of the key submodules. An illustrative working example was shown, that demonstrated the interaction and use of the framework together with the efficiency of the surrogate based approach. Compared to commercial solutions, the modularity makes it suitable for research purposes, while the open-source nature is beneficiary for small companies or users with advanced customization needs.
On top of the correct implementation of the methods, the main challenges related to its use include the selection of a particular surrogate and optimization algorithm or the termination, model selection and convergence criteria. On a fully general scale, the no free lunch (NFL) theorems [51] prove that an universal optimal choice is impossible. However, if the set of problems is limited to an expected class of practical engineering problems, at least a good baseline starting point is possible [53]. Thus, the main aim of the follow-up work is to perform testing on a wide range of problems with a range of framework setups to thoroughly establish the influence of the tunable parameters of the surrogate model as well as of the optimization algorithm.
Further suggestions include a quantitative comparative study of the performance of different proposed adaptive sampling methods, which show the largest potential on the overall SBO performance, as the expensive simulations remains main bottleneck overall. Beyond that, softwarewise, additional suggested features include a graphical user interface (GUI), inclusion of more surrogate models from SMT, such as the RBF, extending the amount of supported evaluators e.g. for Abaqus, and others. Finally, for practice oriented research, a composite layup optimization study in terms of the number of plies and their ply angles is an interesting area of possible application of this framework within the field of advanced structures.