domain: spi phys math info

text:
Mohamed Amine Bouhlel

John T Hwang

Nathalie Bartoli

Rémi Lafage

Joseph Morlier

Joaquim R R A Martins

A Python surrogate modeling framework with derivatives

The surrogate modeling toolbox (SMT) is an open-source Python package that contains a collection of surrogate modeling methods, sampling techniques, and benchmarking functions. This package provides a library of surrogate models that is simple to use and facilitates the implementation of additional methods. SMT is different from existing surrogate modeling libraries because of its emphasis on derivatives, including training derivatives used for gradient-enhanced modeling, prediction derivatives, and derivatives with respect to training data. It also includes unique surrogate models: kriging by partial least-squares reduction, which scales well with the number of inputs; and energyminimizing spline interpolation, which scales well with the number of training points. The efficiency and effectiveness of SMT are demonstrated through a series of examples. SMT is documented using custom tools for embedding automatically tested code and dynamically generated plots to produce high-quality user guides with minimal effort from contributors. SMT is maintained in a public version control repository 1 .

Motivation and significance

In the last few decades, numerical models for engineering have become more complex and accurate, but the computational time for running these models has not necessarily decreased. This makes it difficult to complete engineering tasks that rely on these models, such as design space exploration and optimization. Surrogate modeling is often used to reduce the computational time of these tasks by replacing expensive numerical simulations with approximate functions that are much faster to evaluate. Surrogate models are constructed by evaluating the original model at a set of points, called training points, and using the corresponding evaluations to construct an approximate model based on mathematical functions.

Surrogate modeling is often used in the context of design optimization because of the repeated model evaluations that are required. Derivatives play an important role in optimization, because problems with a large number of optimization variables require gradient-based algorithms for efficient scalability. Therefore, situations frequently arise where there are requirements for surrogate models associated with the computation or use of derivatives.

There are three types of derivatives in surrogate modeling: prediction derivatives, training derivatives, and output derivatives. Prediction derivatives are the derivatives of the surrogate model outputs with respect to the inputs, and they are the derivatives required when using a surrogate model in gradient-based optimization. Training derivatives are derivatives of the training outputs with respect to the training inputs that provide additional training data that increase the accuracy of the surrogate model. Output derivatives are derivatives of the prediction outputs with respect to the training outputs, which are required if a surrogate model is reconstructed within a gradientbased optimization process. We describe these derivatives in more detail in Section 3. [START_REF] Bartoli | Adaptive modeling strategy for constrained global optimization with application to aerodynamic wing design[END_REF].

Various packages that build surrogate models have been developed using different programming languages, such as Scikit-learn in Python [START_REF] Pedregosa | Scikit-learn: Machine learning in Python[END_REF], SUMO in MATLAB [START_REF] Gorissen | A surrogate modeling and adaptive sampling toolbox for computer based design[END_REF], and GPML in MATLAB and Octave [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF]. However, these packages do not handle the derivatives.

In this paper, we introduce a new Python package called the surrogate modeling toolbox (SMT). SMT is different from existing surrogate modeling libraries because of its emphasis on derivatives, including all three types of derivatives described above. SMT also includes newly developed surrogate models that handle derivatives and do not exist elsewhere: partial least-squares (PLS)-based surrogate models [START_REF] Bouhlel | Gradient-enhanced kriging for highdimensional problems[END_REF][START_REF] Bouhlel | An improved approach for estimating the hyperparameters of the kriging model for high-dimensional problems through the partial least squares method[END_REF][START_REF] Bouhlel | Improving kriging surrogates of high-dimensional design models by partial least squares dimension reduction[END_REF], which are suitable for high-dimensional problems, and the regularized minimal-energy tensorproduct spline (RMTS), which is suitable for low-dimensional problems with up to hundreds of thousands of sampling points [START_REF] Hwang | A fast-prediction surrogate model for large datasets[END_REF]. To use SMT, the user should first provide a set of training points. This could be done either by using the sampling techniques and benchmarking functions implemented in SMT or by directly importing the data. Then, the user can build a surrogate model based on the training points and make predictions for the function values and derivatives.

The main goal of this work is to provide a simple Python toolbox that contains a set of surrogate modeling functions and supports different kinds of derivatives that are useful for many applications in engineering. SMT includes sampling techniques, which are necessary for the construction of a surrogate model. Various sample functions are also included to facilitate the benchmarking of different techniques and for reproducing results. SMT is suitable for both novice and advanced users and is supported by detailed documentation available online with examples of each implemented surrogate modeling approach 2 . It is hosted publicly in a version-controlled repository 3 , and is easily imported and used within Python scripts. It is released under the New BSD License and runs on Linux, macOS, and Windows operating systems. Regression tests are run automatically on each operating system whenever a change is committed to the repository.

The remainder of this paper is organized as follows. First, we present the architecture and the main implementation features of SMT in Section 2 and then we describe the methods implemented in SMT in Section 3. Section 4 gives an example of SMT usage and presents a set of benchmarking functions implemented within SMT. We apply SMT to two engineering problems in Section 5 and present the conclusions in Section 6.

2 Software architecture, documentation, and automatic testing SMT is composed of three main modules (sampling methods, problems, and surrogate models) that implement a set of sampling techniques, benchmarking functions, and surrogate modeling techniques, respectively. Each module contains a common interface inherited by the corresponding methods, and each method implements the functions required by the interface, as shown in Figure 1. SMT's documentation is written using reStructuredText and is generated using the Sphinx package for documentation in Python, along with custom extensions 4 . The documentation pages include embedded code snippets that demonstrate the usage. These code snippets are extracted dynamically from actual tests in the source code, ensuring that the code snippets are always up to date. The print output and plots from the code snippets are also generated dynamically by custom Sphinx extensions and embedded in the documentation page. This leads to high-quality documentation with low effort, requiring only the addition of a custom directive and the path to locating the test code. Similarly, another custom directive embeds a table of options, default values, valid types, valid values, and descriptions in the documentation for the surrogate model, sampling method, or benchmarking problem. The documentation also uses existing Sphinx directives to embed descriptions of the user-callable methods in the classes.

In addition to the user documentation, we also provide developer documentation that explains how to contribute code to SMT. The developer documentation includes a different list of application programming interface (API) methods for the Surrogate-Model, SamplingMethod, and Problem classes, which are classes that must be imple- When a developer issues a pull request, the request is merged once the automatically triggered tests run successfully and at least one reviewer approves it. The repository on GitHub5 is linked to two continuous integration testing services, Travis CI (for Linux and macOS) and AppVeyor (for Windows), which trigger test suite execution whenever code is committed and prevent changes from being merged if any of the tests fail.

Surrogate modeling methods

To build a surrogate model, two main steps are necessary. First, we generate a set of training points from the input space where the quantity of interest is computed. This step can be done by using one of the sampling techniques implemented in SMT (Section 3.1), or by uploading an existing training dataset. Second, we train the desired surrogate model on those points and make a prediction of the output values and derivatives (Section 3.2).

Sampling methods

SMT contains a library of sampling methods used to generate sets of points in the input space, either for training or for prediction.

Random sampling: this class creates random samples from a uniform distribution over the design space.

Latin hypercube sampling: this is a statistical method for generating a quasi-random sampling distribution. It is among the most popular sampling techniques in computer experiments thanks to its simplicity and projection properties in highdimensional problems. Five construction criteria are available in SMT: four criteria defined in the pyDOE package 6 and the enhanced stochastic evolutionary criterion [START_REF] Jin | An efficient algorithm for constructing optimal design of computer experiments[END_REF].

Full-factorial sampling: this is a common sampling method where all input variables are set at two levels each. These levels are usually denoted by +1 and -1 for high and low levels, respectively. Full-factorial sampling computes all possible combinations of these levels across all such input variables.

Figure 2 shows the implementation of the Random class, which inherits from the SamplingMethod class.

Surrogate models

Table 1 lists the surrogate modeling methods currently available in SMT and summarizes the advantages and disadvantages of each method. The methods include both well-established methods and methods recently developed by the authors that use derivative information. Among the well-established methods, we implement kriging [START_REF] Sacks | Designs for computer experiments[END_REF], radial basis functions (RBF) [START_REF] Powell | The Theory of Radial Basis Function Approximation in[END_REF], inverse distance weighting (IDW) [START_REF] Shepard | A two-dimensional interpolation function for irregularly-spaced data[END_REF], least squares (LS) [16, ch. 3], and quadratic polynomials (QP) [16, ch. 3]. In Figure 3, we show the implementation of the RBF class, which inherits from the SurrogateModel class.

The methods that we recently developed are kriging combined with partial least squares (KPLS) [START_REF] Bouhlel | Improving kriging surrogates of high-dimensional design models by partial least squares dimension reduction[END_REF], KPLSK for the construction of a standard kriging model in highdimensional problems [START_REF] Bouhlel | An improved approach for estimating the hyperparameters of the kriging model for high-dimensional problems through the partial least squares method[END_REF], gradient-enhanced KPLS (GE-KPLS) [START_REF] Bouhlel | Gradient-enhanced kriging for highdimensional problems[END_REF], and RMTS [START_REF] Hwang | A fast-prediction surrogate model for large datasets[END_REF]. Surrogate models provide a vector of prediction outputs y for a given vector of prediction inputs x ∈ R nx , and can be expressed as

y = f (x, x t , y t ) , (1)
where x t ∈ R nx and y t are the vectors of training inputs and outputs, respectively, which are used to build the surrogate model a priori, and x is the unknown point to predict with the surrogate model. We now describe the surrogate modeling methods available in SMT.

LS and QP

The LS method [START_REF] Hastie | The Elements of Statistical Learning[END_REF] fits a linear model with coefficients β = (β 0 , β 1 , . . . , β nx ), where n x is the number of dimensions, to minimize the residual sum of the squares between the observed responses in the dataset, and the responses predicted by the linear approximation, i.e., min

β ||Xβ -y|| 2 2 , (2)
where X = 1, x (1) T , . . . , x (nt) T T . The vectors in X have (n t × n x + 1) dimensions and n t is the number of training points. The square polynomial model is given by [START_REF] Hastie | The Elements of Statistical Learning[END_REF]

y = Xβ + , (3)
where is a vector of random errors and

X =     1 x (1) 1 . . . x (1) nx x (1) 1 x
(1) 2

. . . x

(1) . . . x

nx-1 x (1) nx x (1) 1 2 . . . x (
(nt) nx x (nt) 1 x (nt) 2 . . . x (nt) nx-1 x (nt) nx x (nt) 1 2 . . . x (nt) nx 2     . ( 4
)
The vector of estimated polynomial regression coefficients using ordinary LS estimation is β = X T X -1 X T y.

(5)

IDW

The IDW model is an interpolating method where the unknown points are calculated with a weighted average of the sampling points using [START_REF] Shepard | A two-dimensional interpolation function for irregularly-spaced data[END_REF]

y =          nt i β x, x (i) t y
(i) t nt i β x, x (i) t , if x = x (i) t ∀i y (i) t , if x = x (i) t for some i , (6)
where x ∈ R nx is the prediction input vector, y ∈ R is the prediction output, x

t ∈ R nx is the input vector for the ith training point, and y (i) t ∈ R is the output value for the ith training point. The weighting function β is defined as

β x (i) , x (j) = ||x (i) -x (j) || -p 2 , (7)
where p is a positive real number called the power parameter, which must be strictly greater than one for the derivatives to be continuous.

RBF

The RBF surrogate model [START_REF] Powell | The Theory of Radial Basis Function Approximation in[END_REF] represents the interpolating function as a linear combination of basis functions, one for each training point. RBF are named as such because the basis functions depend only on the distance from the prediction point to the training point for the basis function. The coefficients of the basis functions are computed during the training stage. RBF are frequently augmented to global polynomials to capture the general trends. The prediction equation for RBF is given by

y = p (x) w p + nt i φ x, x (i) t w r , (8)
where x ∈ R nx is the prediction input vector, y ∈ R is the prediction output, x

t ∈ R nx is the input vector for the ith training point, p (x) ∈ R np is the vector mapping the polynomial coefficients to the prediction output, φ x, x (i) t ∈ R nt is the vector mapping the RBF coefficients to the prediction output, w p ∈ R np is the vector of polynomial coefficients, and w r ∈ R nt is the vector of RBF coefficients. The coefficients, w p and w r , are computed by solving the following linear augmented Gram system:

        φ x (1) t , x (1) t . . . φ x (1) t , x (nt) t p x (1) t T . . . . . . . . . . . . φ x (nt) t , x (1) t
. . . φ x (nt) t , x (nt) t p x (nt) t T p x (1) t . . . p x (nt) t 0              w r1 . . . w rn t w p      =      y (1) t . . . y (nt) t 0      . ( 9
)
Only Gaussian basis functions are implemented currently. These can be written as follows

φ x (i) , x (j) = exp - ||x (i) -x (j) || 2 2 d 2 0 , (10)
where d 0 is a scaling parameter.

Kriging-based models

Kriging [START_REF] Sacks | Designs for computer experiments[END_REF] is an interpolating model that is a linear combination of a known function f i (x) added to a realization of a stochastic process, Z (x), to obtain

ŷ = k i=1 β i f i (x) + Z (x) , (11)
where β i is the ith linear regression coefficients and k is the number of linear regression coefficients to be determined. Here, Z (x) has a mean of zero and a spatial covariance function given by cov Z

x (i) , Z x (j) = σ 2 R x (i) , x (j) , (12)
where σ 2 is the process variance and R is the correlation. Two types of correlation functions are available in SMT: the exponential (Ornstein-Uhlenbeck process)

R x (i) , x (j) = nx l=1 exp -θ l x (i) l -x (j) l
, and the Gaussian correlation function

R x (i) , x (j) = nx l=1 exp -θ l x (i) l -x (j) l 2 ,
where θ l ∈ R + . The number of hyperparameters θ is equal to the number of variables n x . They are estimated by maximizing the likelihood function using the gradientfree optimization algorithm COBYLA [START_REF] Powell | A Direct Search Optimization Method That Models the Objective and Constraint Functions by Linear Interpolation[END_REF]. These two correlation functions are called abs exp (exponential) and squar exp (Gaussian) in SMT. The deterministic term, the sum in the kriging prediction [START_REF] Deb | An efficient constraint handling method for genetic algorithms[END_REF], is replaced by a constant, a linear model, or a quadratic model, all of which are available in SMT. KPLS [START_REF] Bouhlel | Improving kriging surrogates of high-dimensional design models by partial least squares dimension reduction[END_REF] is a kriging model that uses the PLS method. KPLS is faster than kriging because fewer hyperparameters need to be estimated to achieve high accuracy. This model is suitable for high-dimensional problems owing to the kernel constructed through the PLS method. The PLS method is a well-known tool for high-dimensional problems that searches the direction that maximizes the variance between the input and output variables. This is done by a projection in a smaller space spanned by the principal components. The PLS information is integrated into the kriging correlation matrix to scale the number of inputs by reducing the number of hyperparameters. The number of principal components, h, which corresponds to the number of hyperparameters in KPLS is much lower than n x . For example, the PLS-Gaussian correlation function is

R x (i) , x (j) = nx l=1 h k=1 exp -θ k w (k) 2 l x (i) l -x (j) l 2 . ( 13
)
We provide more details on the KPLS method in previous work [START_REF] Bouhlel | Improving kriging surrogates of high-dimensional design models by partial least squares dimension reduction[END_REF].

The KPLSK model is built in two steps [START_REF] Bouhlel | An improved approach for estimating the hyperparameters of the kriging model for high-dimensional problems through the partial least squares method[END_REF]. The first step is to run KPLS and estimate the hyperparameters expressed in the reduced space with h dimensions. The second step is to express the estimated hyperparameters in the original space with n x dimensions, and then use this as a starting point to optimize the likelihood function of a standard kriging model. The idea here is to guess a "good" initial hyperparameter and apply a gradient-based optimization using a classic kriging kernel. This guess is provided by the KPLS construction: the solutions (θ * 1 , . . . , θ * h ) and the PLS-coefficients w

η l = h k=1 θ * k w (k) l 2 (14)
for l = 1, . . . , n x , we express the initial hyperparameters point in the original space. We demonstrate this using a KPLS-Gaussian correlation function R KPLS in the following example [START_REF] Bouhlel | An improved approach for estimating the hyperparameters of the kriging model for high-dimensional problems through the partial least squares method[END_REF]:

R KPLS x (i) , x (j) = h k=1 nx l=1 exp -θ k w (k) l 2 x (i) l -x (j) l 2 = exp nx l=1 h k=1 -θ k w (k) l 2 x (i) l -x (j) l 2 = exp nx l=1 -η l x (i) l -x (j) l 2 = nx l=1 exp -η l x (i) l -x (j) l 2 , (15)
where the change of variables is performed on the third line, and the last line is a standard Gaussian correlation function. The hyperparameters point ( 14) provides a starting point for a gradient-based optimization applied on a standard kriging method.

3.2.5 Gradient-enhanced KPLS models GE-KPLS [START_REF] Bouhlel | Gradient-enhanced kriging for highdimensional problems[END_REF] is a gradient-enhanced kriging (GEK) model with the PLS approach.

GEK is an extension of kriging that exploits gradient information. GEK is usually more accurate than kriging; however, it is not computationally efficient when the number of inputs, the number of sampling points, or both, are high. This is primarily owing to the size of the corresponding correlation matrix that increases in proportion to both the number of inputs and the number of sampling points.

To address these issues, GE-KPLS exploits the gradient information with a slight increase in the size of the correlation matrix and reduces the number of hyperparameters. The key idea of GE-KPLS is to generate a set of approximating points around each sampling point using a first-order Taylor approximation. Then, the PLS method is applied several times, each time on a different number of sampling points with the associated sampling points. Each PLS provides a set of coefficients that contribute to each variable nearby the associated sampling point to the output. Finally, an average of all PLS coefficients is computed to estimate the global influence to the output. Denoting these coefficients by w

(k) 1 , . . . , w (k)
nx , the GE-KPLS Gaussian correlation function is given by

R GE-KPLS x (i) , x (j) = σ nx l=1 h k=1 exp -θ k w (k) l x (i) l -w (k) l x (j) l 2 . ( 16
)
This approach reduces the number of hyperparameters (reduced dimension) from n x to h, where h << n x . As mentioned previously, PLS is applied several times with respect to each sampling point, which provides the influence of each input variable around that point. The idea here is to add only m approximating points (m ∈ [1, n x ]) around each sampling point. Only the m highest coefficients given by the first principal component are considered, which usually capture the most useful information. Bouhlel and Martins [START_REF] Bouhlel | Gradient-enhanced kriging for highdimensional problems[END_REF] provide more details on this approach.

RMTS method

The RMTS model is a type of surrogate model for low-dimensional problems with large datasets that has fast prediction capability [START_REF] Hwang | A fast-prediction surrogate model for large datasets[END_REF]. The underlying mathematical functions are tensor-product splines, which limits RMTS to up to four-dimensional problems, or five-dimensional problems in certain cases. On the other hand, tensorproduct splines enable a fast prediction time that does not increase with the number of training points. Unlike other methods, such as kriging and RBF, RMTS is not susceptible to numerical issues when there is a large number of training points or when points are too close together. The prediction equation for RMTS is given by

y = F (x) w, (17)
where x ∈ R nx is the prediction input vector, y ∈ R is the prediction output, w ∈ R nw is the vector of spline coefficients, and F (x) ∈ R nw is the vector mapping the spline coefficients to the prediction output. RMTS computes the coefficients of the splines, w, by solving an energy minimization problem subject to the conditions that the splines pass through the training points. This is formulated as an unconstrained optimization problem where the objective function consists of a term with the second derivatives of the splines, a term that represents the approximation error for the training points, and another term for regularization. Thus, this optimization problem can be written as min

w 1 2 w T Hw + 1 2 βw T w + 1 2 1 α nt i F x (i) t w -y (i) t 2 , (18)
where

x (i)
t ∈ R nx is the input vector for the ith training point, y

t ∈ R is the output value for the ith training point, H ∈ R nw×nw is the matrix of second derivatives, F x (i) t ∈ R nw is the vector mapping the spline coefficients to the ith training output, and α and β are the regularization coefficients.

In problems with a large number of training points relative to the number of spline coefficients, the energy minimization term is not necessary and can be set to zero by setting the reg cons option to zero. In problems with a small dataset, the energy minimization is necessary. When the true function has high curvature, the energy minimization can be counterproductive. This can be addressed by increasing the quadratic approximation term to one of higher order and using Newton's method to solve the resulting nonlinear system. The nonlinear formulation is given by min

w 1 2 w T Hw + 1 2 βw T w + 1 2 1 α nt i F x (i) t w -y (i) t p , ( 19
)
where p is the order set by the approx order option. The number of Newton iterations is specified in the nonlinear maxiter option. RMTS is implemented in SMT with two choices of splines: B-splines and cubic Hermite splines. RMTB uses B-splines with a uniform knot vector in each dimension. The number of B-spline control points and the B-spline order in each dimension are options that trade off efficiency and precision of the interpolant. For the cubic Hermite splines, RMTC divides the domain into tensor-product cubic elements. For adjacent elements, the values and derivatives are continuous. The number of elements in each dimension is an option that trades off efficiency and precision. B-splines are usually the better choice when training time is the most important factor, whereas cubic Hermite splines are the better choice when the accuracy of the interpolant is most important [START_REF] Hwang | A fast-prediction surrogate model for large datasets[END_REF].

Derivatives

There are three types of derivatives in SMT.

Prediction derivatives ( dy/ dx) are derivatives of predicted outputs with respect to the inputs at which the model is evaluated. These are computed together with the prediction outputs when the surrogate model is evaluated [START_REF] Li | Data-based approach for fast airfoil analysis and optimization[END_REF]. These are required for gradient-based optimization algorithms based on surrogate models [START_REF] Bartoli | An adaptive optimization strategy based on mixture of experts for wing aerodynamic design optimization[END_REF][START_REF] Bartoli | Adaptive modeling strategy for constrained global optimization with application to aerodynamic wing design[END_REF].

Training derivatives ( dy t / dx t ) are derivatives of the training outputs with respect to the corresponding inputs. These are provided by the user and are used to improve the model accuracy in GE-KPLS.

When the adjoint method is used to compute training derivatives, a high-quality surrogate model can be constructed with a low relative cost, because the adjoint method computes these derivatives at a cost independent of the number of inputs.

Output derivatives ( dy/ dy t ) are derivatives of predicted outputs with respect to training outputs, which is a measure of how the prediction changes with a change in training outputs, accounting for the retraining of the surrogate model. These post-training derivatives are used when the surrogate model is trained within an optimization iteration. This feature is not commonly available in other frameworks; however, it is required when the training of the surrogate model is embedded in a gradient-based optimization. In this case, derivatives of the prediction outputs with respect to the training outputs must be computed and combined with derivatives from other parts of the model using, for example, the chain rule.

Given its focus on derivatives, SMT is synergistic with the OpenMDAO framework [START_REF] Gray | OpenMDAO: An open-source framework for multidisciplinary design, analysis, and optimization[END_REF], which is a software framework for gradient-based multidisciplinary analysis and optimization [START_REF] Hwang | A computational architecture for coupling heterogeneous numerical models and computing coupled derivatives[END_REF][START_REF] Martins | Review and unification of methods for computing derivatives of multidisciplinary computational models[END_REF][START_REF] Martins | Multidisciplinary design optimization: A survey of architectures[END_REF]. An SMT surrogate model can be a component that is part of a larger model developed in OpenMDAO and can provide the derivatives that OpenMDAO requires from its components to compute the coupled derivatives of the multidisciplinary model.

Additional surrogate modeling methods

To extend surrogate modeling to higher-level methods, we implemented a component within SMT named extensions. These methods require additional and sometimes different steps than the usual surrogate modeling. For example, multi-fidelity methods combine data generated from different sources, such as coarse and fine mesh solutions. Another example is the mixture of experts (MoE) class of methods, which linearly combines several surrogates. Three such methods are available within SMT.

MoE modeling performs a weighted sum of local surrogate models (experts) instead

of one single model. It is based on splitting the input space into several subspaces via clustering algorithms and training a surrogate model within each subspace. Hastie et al. [START_REF] Hastie | The Elements of Statistical Learning[END_REF] provide a general introduction to the MoE method. The implementation of MoE within SMT is based on Gaussian mixture models and expectation maximization [START_REF] Bettebghor | Surrogate modeling approximation using a mixture of experts based on em joint estimation[END_REF].

Variable-fidelity modeling (VFM) samples points using both low-fidelity function evaluations that are cheap to evaluate and more costly high-fidelity evaluations.

It then uses the differences between the high-and low-fidelity evaluations to construct a bridge function that corrects the low-fidelity model [START_REF] Han | Improving variable-fidelity surrogate modeling via gradient-enhanced kriging and a generalized hybrid bridge function[END_REF]. SMT implements additive and multiplicative bridge functions.

Multi-fidelity kriging (MFK) uses a correlation factor ρ(x) (constant, linear, or quadratic) and a discrepancy function δ(x). The high-fidelity model is given by

y high (x) = ρ(x)y low (x) + δ(x).
SMT follows the formulation proposed by Le Gratiet [START_REF] Gratiet | Multi-fidelity Gaussian process regression for computer experiments[END_REF], which is recursive and can be easily extended to multiple levels of fidelity.

Example and benchmark functions

SMT contains a library of analytical and engineering problems that can be used for instructional or benchmarking purposes. The analytical functions currently implemented in SMT are the sphere, Branin, L p norm, Rosenbrock, and tensor-product functions, all of which are detailed in A. The engineering functions currently implemented in SMT are the cantilever beam, robot arm, torsion vibration, water flow, welded beam, and wing weight problems (described in B).

Figure 4 shows the implementation of the Sphere class, which inherits from the Problem class. This function is a simple example that serves as a good first test for newly developed surrogate modeling and surrogate-based optimization methods. This function is a continuous and convex function where the global minimum is at the origin. Figure 5 shows an example of the use of RMTC within SMT for the robot arm function [START_REF] An | Quasi-regression[END_REF]. This function is scalable; however, the number of dimensions must be an even number (we use two dimensions in this example). This function gives the position of a robot arm, which is made by multiple segments, in a two-dimensional space where the shoulder of the arm is fixed at the origin. This function is highly nonlinear, and the use of training derivatives samples is particularly beneficial in this case [START_REF] Bouhlel | Gradient-enhanced kriging for highdimensional problems[END_REF].

Applications

In this section, we describe two applications that highlight the unique features of SMT. We also provide an overview of the main previous applications realized using SMT. The first application is the computation and validation of the prediction derivatives for an airfoil analysis tool, where we apply GE-KPLS. The second application describes a practical use for output derivatives, where we use RMTS to compute the outputs derivatives for a surrogate model that is dynamically trained with an optimization iteration.

Airfoil analysis and shape optimization tool

Li et al. [START_REF] Li | Data-based approach for fast airfoil analysis and optimization[END_REF] used the GE-KPLS implementation in SMT with an MoE technique to develop a data-driven approach to airfoil and shape optimization based on computational fluid dynamics (CFD) simulations 7 . They used a database of 1100 existing Functions values if kx = None or derivative values if kx is an int . """ ne , nx = x . shape y = np . zeros (( ne , 1) , complex ) if kx is None : y [: , 0] = np . sum ( x **2 , 1). T else : y [: , 0] = 2 * x [: , kx ] return y airfoils8 and enriched the database with 100,000 more generated airfoils. The data was used to create a surrogate model of force coefficients (lift, drag, and moment) with respect to flight speed (Mach number), the angle of attack, and airfoil shape variables. The shape design variables consisted of shape modes obtained by singular-value decomposition of the existing airfoil shapes. The approach proved successful, enabling shape optimization in 2 sec using a personal computer with an error of less than 0.25% for subsonic flight conditions. This is thousands of times faster than optimization using direct calls to CFD.

To validate the prediction derivatives provided by SMT, we built a surrogate model of the airfoil drag coefficient similar to the one above using the following four inputs: first thickness mode t 1 , first camber mode c 1 , Mach number M , and angle of attack α. We compute the partial derivatives of the drag coefficient (C d ) with respect to these To train the surrogate model, we use 20 training points with their respective training derivatives. We use the GE-KPLS model, where we provide training derivatives using an adjoint approach [START_REF] Lyu | Automatic differentiation adjoint of the Reynolds-averaged Navier-Stokes equations with a turbulence model[END_REF][START_REF] Mader | ADjoint: An approach for the rapid development of discrete adjoint solvers[END_REF]. To validate the prediction derivatives of the constructed surrogate model, we generate 100 new points different from the training points and compute the relative error on those points using

ε = || ∂C d ∂x -∂ C d ∂x || 2 || ∂C d ∂x || 2 , (20)
where ∂ Ĉd /∂x is the vector of 100 surrogate model prediction derivatives with respect to one of the four variables (x = t 1 , c

Aircraft design optimization considering dynamics

The design optimization of commercial aircraft using CFD is typically done by modeling the aircraft performance at a small number of representative flight conditions [START_REF] Kenway | Multipoint high-fidelity aerostructural optimization of a transport aircraft configuration[END_REF][START_REF] Kenway | Multipoint aerodynamic shape optimization investigations of the Common Research Model wing[END_REF]. This simplification is made because it would be prohibitively expensive to simulate the full, discretized mission profile using CFD at all points. However, full-mission simulation is sometimes necessary, such as when the flight is short or when considering morphing aircraft designs [START_REF] Burdette | Impact of morphing trailing edge on mission performance for the common research model[END_REF][START_REF] Martins | Encyclopedia of Aerospace Engineering, volume Green Aviation, chapter Fuel burn reduction through wing morphing[END_REF]. This requires a surrogate model of the aerodynamic performance as a function of a small number of parameters such as flight speed and altitude. In a design optimization context, the surrogate model must be retrained each optimization iteration because the aircraft design changes from iteration to iteration as the shape is being optimized. This is the situation in an aircraft allocation-mission-design optimization problem recently solved using RMTS [START_REF] Hwang | Solution of ordinary differential equations in gradient-based multidisciplinary design optimization[END_REF]. In this work, the optimization problem maximized airline profit by optimizing the twist, span, sweep, and shape of the wing of a nextgeneration commercial aircraft. The profit was computed by simulating the fuel efficiency and flight time for the aircraft on a set of routes operated by the hypothetical airline. To do this, a surrogate model was generated for the wing lift and drag coefficients (C L and C D ) as a function of the angle of attack (α), Mach number (M ), and altitude (h). In each optimization iteration, the training C L and C D values were computed at a series of points in M -α-h space using CFD. The training outputs were used to retrain the RMTS surrogate model, and the mission simulations for all the airline routes were performed using inexpensive evaluations of the trained surrogate model, as shown in Figure 7. However, because the overall optimization problem was solved using a gradient-based algorithm, we required derivatives of the prediction outputs (C L and C D values at the M -α-h points at which we evaluated the surrogate) with respect to the training outputs (C L and C D values at the fixed M -α-h points where the training points were located).

Figure 8 shows the surrogate model with the altitude axis eliminated by projecting onto the other two axes, where C L replaces α as one of the inputs. The horizontal and vertical axes represent the inputs. The red points are the training points and the black points are the prediction points. Therefore, the output derivatives are the derivatives of the outputs at the black points with respect to the training outputs at the red points.

Other applications

In addition to the two applications we just described, SMT has been used to solve other engineering problems. We summarize these applications in Table 3. The number of input variables ranges from 2 to 99. This demonstrates SMT's ability to solve different engineering problems of various complexities. SMT arose from our own research needs, namely the development of the GE-KPLS and RMTS surrogate models, which require the handling of derivatives. GE-KPLS is a method we recently developed that scales much better with the number of inputs than the other methods, thanks in part to the use of the derivative information. RMTS is another method that we recently developed that is best suited to low-dimensional problems and handles both structured and unstructured training points.

SMT is distributed through an open-source license and is freely available online9 . We provide documentation that caters to both users and potential developers 10 . SMT allows users collaborating on the same project to have a common surrogate modeling tool that facilitates the exchange of methods and reproducibility of results. SMT is a recently developed tool and, so far, it has primarily been used in aerospace engineering applications. However, this tool is useful to anyone needing to use or develop surrogate modeling techniques, regardless of the application. In particular, given the advantages of the newer surrogate modeling techniques (GE-KPLS and RMTS), using SMT is particularly attractive.

A Analytical functions implemented in SMT

In the following descriptions, n x is the number of dimensions.

Sphere The sphere function is quadratic, continuous, convex, and unimodal. It is given by nx i=1

x 2 i , -10 ≤ x i ≤ 10, for i = 1, . . . , n x .

Branin [START_REF] Forrester | Engineering Design via Surrogate Modelling-A Practical Guide[END_REF] The Branin function is commonly used in optimization and has three global minima. It is given by Rosenbrock [START_REF] Shang | A note on the extended Rosenbrock function[END_REF] The Rosenbrock function is a continuous, nonlinear, and nonconvex function and usually used in optimization. The minimum of this function is situated in a parabolic-shaped flat valley and is given by nx-1 i=1

f (x) = x 2 - 5
x i+1 -x 2 i 2 + (x i -1) 2 , -2 ≤ x i ≤ 2, for i = 1, . . . , n x .

Tensor-product The tensor-product function approximates a step function, which causes oscillations with some surrogate models [START_REF] Hwang | A fast-prediction surrogate model for large datasets[END_REF], and is given by

B Engineering problems implemented in SMT

Cantilever beam [START_REF] Cheng | Trust region based mode pursuing sampling method for global optimization of high dimensional design problems[END_REF] This function models a simple uniform cantilever beam with vertical and horizontal loads 50 600 Robot arm [START_REF] An | Quasi-regression[END_REF] This function gives the position of a robot arm Torsion vibration [START_REF] Liping | A comparison of metamodeling methods using practical industry requirements[END_REF] This function gives the low natural frequency of a torsion problem

1 2π -b - √ b 2 -4ac 2a ,
where K i = πG i d i 32L i , M j = ρ j πt j D j 4g

, J j = 0.5M j Water flow [START_REF] Morris | Bayesian design and analysis of computer experiments: Use of derivatives in surface prediction[END_REF] This function characterizes the flow of water through a borehole that is drilled from the ground surface through two aquifers Welded beam [START_REF] Deb | An efficient constraint handling method for genetic algorithms[END_REF] The shear stress of a welded beam problem is given by Wing weight [START_REF] Forrester | Engineering Design via Surrogate Modelling-A Practical Guide[END_REF] The estimate of the weight of a light aircraft wing is given by 0.036S

D j 2 , a = 1, b = -K 1 +K2 J 1 + K 2 +K3 J 2 , c = K 1 K 2 +K 2 K 3 +K 3 K 1 J 1 J 2 ,
τ 2 + τ 2 +

Figure 1 :

 1 Figure 1: Architecture of SMT.

Figure 2 :

 2 Figure 2: Implementation of the Random class, which inherits from the Sampling-Method class. The compute function evaluates the requested number of sampling points uniformly over the design space.

Figure 3 :

 3 Figure 3: Implementation of the RBF class, which inherits from the SurrogateModel class.

  nx for k = 1, . . . , h. Using the change of variables

Figure 4 :

 4 Figure 4: Implementation of the Sphere class, which inherits from the Problem class.The main function is evaluate, which computes either the output values or the derivatives depending on the variable kx.

Figure 5 :

 5 Figure 5: Example of use of RMTC within SMT on the robot arm function [1].

Figure 6 :

 6 Figure6: Top left: drag and lift coefficients for the baseline, the solution of Optimization 1, and the solution of Optimization 2. Bottom left: airfoil shapes for all three cases; the shapes for both optimizations is indistinguishable. Right: pressure distributions. Optimization 1 achieved a 0.03% lower drag compared with Optimization 2 using a fraction of the evaluations.

Fig. 8

 8 Fig. 8 The contours show L/D, the red markers indicate training point locations, and the black markers and lines indicate the path of the aircraft through C L -M space for the longest-range mission. The L/D = 26 contour is highlighted in grey.

Figure 8 :

 8 Figure 8: RMTS surrogate model for the performance of an aircraft, represented by the lift-to-drag ratio (L/D), as a function of Mach number (M ), and lift coefficient (C L ) [19]. In this context, an output derivative is the derivative of the L/D values at the prediction points (shown in black) with respect to the L/D values at the training points (shown in red). The C L -M values of the red points are the training inputs, which are fixed, and the C L -M values of the black points are the prediction inputs. The path of the black points represents the flight of a commercial aircraft, from climb with C L ∼ 0.5 to cruise at M > 0.8 to descent with C L ∼ 0.3. The red points appear to be duplicated because there is a third input, altitude, which is eliminated by projecting onto the C L -M plane.

  x = (x 1 , x 2 ) with -5 ≤ x 1 ≤ 10 and 0 ≤ x 2 ≤ 15.

  i ) , -1 ≤ x i ≤ 1, for i = 1, . . . , n x , or nx i=1 exp (x i ) , -1 ≤ x i ≤ 1, for i = 1, . . . , n x , or nx i=1 tanh (x i ) , -1 ≤ x i ≤ 1, for i = 1, . . . , n x , or nx i=1 exp -2x 2 i , -1 ≤ x i ≤ 1, for i = 1, . . . , n x .

  where b i ∈ [0.01, 0.05], h i ∈ [0.3, 0.65], l i ∈ [0.5, 1].

4 i=1 L i cos i j=1 θ j 2 + 4 i=1 L i sin i j=1 θ j 2 ,

 4242 where L i ∈ [0, 1] for i = 1, . . . , 4 and θ j ∈ [0, 2π] for j = 1, . . . , 4.

  2πT u (H u -H l ) 05 ≤ r w ≤ 0.15, 100 ≤ r ≤ 50000, 63070 ≤ T u ≤ 115600, 990 ≤ H u ≤ 1110, 63.1 ≤ T l ≤ 116, 700 ≤ H l ≤ 820, 1120 ≤ L ≤ 1680, and 9855 ≤ K w ≤ 12045.

2 ,

 2 for h ∈ [0.125, 1], and l, t ∈ [5, 10].

Table 1 :

 1 Surrogate modeling methods provided by SMT.

	Method	Advantages (+) and disadvantages (-)		Derivatives		References
			Train. Pred. Out.
	Kriging	+ Prediction variance, flexible	No	Yes	No	Sacks et al. [39]
		-Costly if number of inputs or training points is large
		-Numerical issues when points are too close to each other
	KPLS	+ Prediction variance, fast construction	No	Yes	No	Bouhlel et al. [7]
		+ Suitable for high-dimensional problems
		-Numerical issues when points are too close to each other
	KPLSK	+ Prediction variance, fast construction	No	Yes	No	Bouhlel et al. [6]
		+ Suitable for high-dimensional problems
		-Numerical issues when points are too close to each other
	GE-KPLS + Prediction variance, fast construction	Yes	Yes	No	Bouhlel and Martins [5]
		+ Suitable for high-dimensional problems
		+ Control of the correlation matrix size
		-Numerical issues when points are too close to each other
		-Choice of step parameter is not intuitive
	RMTS	+ Fast prediction	Yes	Yes	Yes	Hwang and Martins [18]
		+ Training scales well up to 10 5 training points
		+ No issues with points that are too close to each other
		-Poor scaling with number of inputs above 4
		-Slow training overall
	RBF	+ Simple, only a single tuning parameter	No	Yes	Yes	Powell [36]
		+ Fast training for small number of training points
		-Susceptible to oscillations
		-Numerical issues when points are too close to each other
	IDW	+ Simple, no training required	No	Yes	Yes	Shepard [41]
		-Derivatives are zero at training points
		-Poor overall accuracy
	LS	+ Simple, fast construction	No	Yes	No	Hastie et al. [16]
		-Accurate only for linear problems
	QP	+ Simple, fast construction	No	Yes	No	Hastie et al. [16]
		-Large number of points required for large number of inputs

Table 2 :

 2 1 , M, α), ∂C d /∂x is the vector of corresponding reference derivative values with respect to x, and ||.|| 2 is the L 2 norm. The resulting relative errors corresponding to each variable (listed in Table 2) show that the surrogate model provides an accurate estimate of the derivatives with a relative error less or equal to 10 -2 . Relative error of the prediction derivative of C d with respect to t 1 , c 1 , M , and α. The surrogate model yields an accurate prediction derivative with an error less than or equal to 10 -2 . ∂C d /∂t 1 ∂C d /∂c 1 ∂C d /∂M ∂C d

	/∂α

  Extended design structure matrix[START_REF] Lambe | Extensions to the design structure matrix for the description of multidisciplinary design, analysis, and optimization processes[END_REF] for the aircraft allocation-missiondesign optimization problem that used RMTS[START_REF] Hwang | Solution of ordinary differential equations in gradient-based multidisciplinary design optimization[END_REF]. The aerodynamic surrogate model predicts wing lift and drag coefficient as a function of the angle of attack, Mach number, and altitude. This surrogate model is retrained each iteration because the lift and drag training data change each optimization iteration as the aircraft wing design changes. Since a gradient-based algorithm is used for this large-scale optimization problem, we require output derivatives: derivatives of the lift and drag prediction outputs with respect to the lift and drag training data. These are provided by RMTS.

	Optimizer	Area, sweep, twist, shape	Area, sweep, twist, shape	Altitude, cruise Mach	Altitude	Altitude		Allocation variables
	Geometric	Geometric
	constraints	model
			CFD solver		CFD training data
				Atmospheric models	Mach	Mach	Velocity, density
					Aerodynamic surrogate		Drag
						Propulsion surrogate	SFC
	Thrust constraints				Lift	Thrust	Mission analysis eq.	Fuel burn, block time
	Profit, alloc.							Allocation
	constraints							model
	Figure 7: information includes training derivatives used for gradient-enhanced modeling, pre-
	diction derivatives used for surrogate-model-based optimization, and derivatives with
	respect to the training data used when the optimization loop includes reconstructing
	the surrogate model. However, SMT does not need to involve derivatives and pro-
	vides a simple general interface to various surrogate modeling techniques. Together
	with the provided problems, SMT is a convenient framework to develop new methods
	for sampling and surrogate modeling, and then benchmarking them against existing
	methods.

Table 3 :

 3 Summary of SMT applications to engineering design problems.

	Problem	Surrogate	Objective of the study	Number of	Reference
		used		design variables
	3D blade	KPLS	Compare the accuracy and efficiency of the	Up to 99	Bouhlel et al. [7]
			KPLS model and the Optimus [34]
			implementation of the kriging model
	3D blade	KPLSK	Compare the accuracy and efficiency of the	99	Bouhlel et al. [6]
			KPLSK model and the KPLS model
	Automotive	KPLS,	Minimize the mass of a vehicle subject to	50	Bouhlel et al. [8]
		KPLSK	68 constraints
	Eight different GE-KPLS Compare the accuracy and efficiency of the	Up to 15	Bouhlel and Martins [5]
	engineering		GE-KPLS model and the indirect GEK
	problems
	Aircraft wing	MoE	Minimize the drag of aircraft wings subject	Up to 17	Bartoli et al. [2]
			to lift constraints
	2D airfoils	GE-KPLS Build a fast interactive airfoil analysis and	Up to 16	Li et al. [26]
			design optimization tool
	Aircraft	RMTS	Optimize an aircraft design while analyzing 3	Hwang and Munster [19]
	performance		the full mission using a CFD surrogate
	Rotor	RMTS	Optimize a rotor using a surrogate model	2	Hwang and Ning [20]
	analysis		for the lift and drag of blade airfoils

  for d 1 ∈ [1.8, 2.2], L 1 ∈ [9, 11], G 1 ∈ [105300000, 128700000], d 2 ∈ [1.638, 2.002], L 2 ∈ [10.8, 13.2], G 2 ∈ [5580000, 6820000], d 3 ∈ [2.025, 2.475], L 3 ∈ [7.2, 8.8], G 3 ∈ [3510000, 4290000], D 1 ∈ [10.8, 13.2], t 1 ∈ [2.7, 3.3], ρ 1 ∈ [0.252, 0.308], D 2 ∈ [12.6, 15.4], t 2 ∈ [3.6, 4.4], and ρ 1 ∈ [0.09, 0.11].

  (N z W dg ) 0.49 + S w W p , where 150 ≤ S w ≤ 200, 220 ≤ W f w ≤ 300, 6 ≤ A ≤ 10, -10 ≤ Λ ≤ 10, 16 ≤ q ≤ 45, 0.5 ≤ λ ≤ 1, 0.08 ≤ t c ≤ 0.18, 2.5 ≤ N z ≤ 6, 1700 ≤ W dg ≤ 25000, and 0.025 ≤ W p ≤ 0.08.

	0.758 w	W 0.0035 f w	A cos 2 Λ	cos Λ q 0.006 λ 0.04 100tc	-0.3

http://smt.readthedocs.io/en/latest

https://github.com/SMTorg/smt

https://smt.readthedocs.org

https://github.com/SMTorg/smt

https://pythonhosted.org/pyDOE/randomized.html

https://github.com/mdolab/adflow

http://webfoil.engin.umich.edu

ConclusionsSMT is unique compared with existing surrogate modeling libraries in that it is designed from the ground up to handle derivative information effectively. The derivative

https://github.com/SMTorg/SMT

http://smt.readthedocs.io/en/latest/

Acknowledgments

This work was partially supported by the Air Force Office of Scientific Research (AFOSR) MURI on "Managing multiple information sources of multi-physics systems," Program Officer Jean-Luc Cambier, Award Number FA9550-15-1-0038. This work is part of the activities of ONERA-ISAE-ENAC joint research group and was partially supported by an ONERA internal project on multi-fidelity methods (MUFIN). The authors would like to thank Rémi Vauclin and Mostafa Meliani for the implementation of MFK, and Rémy Priem for the implementation of MoE.

//github.com/SMTorg/SMT