- Email: [email protected]
Construction of response surfaces based on progressive-lattice-sampling experimental designs with application to uncertainty propagation
Structural Safety 26 (2004) 201–219 www.elsevier.com/locate/strusafe
Construction of response surfaces based on progressivelattice-sampling experimental designs with application to uncertainty propagation V.J. Romero*, L.P. Swiler, A.A. Giunta Sandia National Laboratories, Albuquerque, NM 87185-0828, USA
Abstract Response surface functions are often used as simple and inexpensive replacements for computationally expensive computer models that simulate the behavior of a complex system over some parameter space. Here we examine several data fitting and interpolation techniques (finite-element interpolation, kriging, and polynomial regression) that can be used to construct a sequence of progressively upgraded response surface approximations based on Progressive Lattice Sampling (PLS) incremental experimental designs. PLS is a paradigm for structured sampling of a hypercube parameter space by placing and incrementally adding samples at each level of the design in a manner intended to efficiently leverage the samples at all previous levels. When combined with compatible interpolation methods, PLS can be used to construct efficiently upgradable response surface approximations. Upon upgrading, convergence heuristics are gained that can be used to estimate the magnitude of the approximation error entailed in the response surface. The three interpolation methods tried here are examined for fitting and convergence behavior in several test problems. Published by Elsevier Ltd. Keywords: Response surfaces; Experimental design; Fitting and interpolation; Uncertainty analysis
1. Introduction Response surface functions are often used as simple and inexpensive replacements of complex and computationally expensive computer models that simulate the behavior of a system over some parameter space. As such, these ‘‘surrogate’’ response surface models fit and interpolate data that is generated by running the simulation model at selected points in the parameter space, e.g. based on some experimental design. The response surface is then used to predict system response over the parameter space at reduced cost, say for uncertainty analysis where the * Corresponding author. Tel.: +1-505-844-5890; fax: +1-505-844-8251. E-mail address: [email protected] (V.J. Romero). 0167-4730/$ - see front matter Published by Elsevier Ltd. doi:10.1016/j.strusafe.2003.03.001
202
V.J. Romero et al. / Structural Safety 26 (2004) 201–219
response surface may be evaluated hundreds or thousands of times. Besides expense, problems associated with numerical noise in the underlying simulation models, and model coupling in multidisciplinary analysis also motivate the use of response surfaces. Examples of response surface usage to facilitate large-scale optimization and uncertainty analysis are cited in references [1–4]. Two issues that arise when using response surface functions are accuracy and the cost of procuring the data samples needed to create the functions. With a sufficiently compliant global fitting/interpolating function over the parameter space, response surface accuracy generally increases as the number of data points increases (if the points are appropriately placed throughout the parameter space), until the essential character of the function is effectively mapped out. Thereafter, it is not cost effective to continue adding samples. Since a single high-fidelity physics simulation (i.e., one data sample) can take hours to compute, it is desirable to minimize the number of simulations that are needed to construct an accurate response surface. A method is presented in [5] to address the accuracy and computational cost associated with using response surface functions under the following circumstances: (1) the physics simulation model is relatively expensive to evaluate; (2) the parameter space is a hypercube or can be accurately and inexpensively mapped into one; (3) the sampled or ‘‘target’’ function is a continuous, deterministic function over the parameter space; (4) reasonably general, arbitrary target functions are to be fitted; and (5) approximate response values are desired over the entire parameter space—i.e., for uncertainty propagation problems or global optimization. In the Finite-Element/Progressive-Lattice-Sampling method, the response surface is a continuous, piecewise-smooth global interpolation function comprised of a patchwork of linear to quadratic polynomials over adjoining subspaces, in the manner of ‘‘finite element’’ representation of spatial functions or ‘‘fields’’ in computational continuum mechanics. The data samples fit and interpolated by the finite-element (FE) response surfaces are arranged on a structured lattice. The so-called ‘‘Progressive Lattice Sampling’’ (PLS) method of generating the lattice points may be viewed as an ‘‘incremental’’ experimental design method. The PLS design attempts to place and add sets of samples such that all samples are efficiently leveraged as the design progresses from one level to the next. Hence, the FE/PLS combination can be used to construct efficiently upgradable response surfaces. A mathematical analysis of finite element interpolation [6] shows that for a continuous and infinitely differentiable function over the parameter space, the domain integral of the pointwise absolute error goes to zero as the spacing between samples goes to zero. This analysis can be applied to the FE/PLS method. Hence, in the limit, FE/PLS response surfaces converge everywhere in the domain to target functions in this class (though the convergence rate is generally not uniform over the domain). Upon upgrading a response surface to the next level of the PLS design, the resulting change in predicted response at any point in the domain is a heuristic indicator of the magnitude of the local error in the response surface approximation. When the incremental change goes to zero, this indicates that the local approximation error has become negligible. The purpose of the current study is to compare the performance of the FE interpolation method against two other response surface construction methods (polynomial regression and kriging) that are compatible with upgradable PLS experimental designs. The convergence behaviors of the three interpolation schemes are examined on several 2-D test problems as the number of data samples is increased in passing through the various levels of the PLS design.
V.J. Romero et al. / Structural Safety 26 (2004) 201–219
203
2. Progressive lattice sampling (PLS) experimental designs As presented here, the implementation of Progressive Lattice Sampling has changed somewhat since it was first introduced in [5], but the essence is still the same. The basic concept (which is also utilized in e.g., upgradable quadrature [7,8] and adaptive nodalization for error reduction in the solution of partial differential equations, [9]) is to progressively sample the parameter space in a manner that efficiently builds on previous information (i.e., fully leverages previous samples). Under the present circumstances this means placing and adding samples in a manner that preserves uniformity of sampling coverage over the hypercube parameter space in the various stages or levels of the incremental experimental design (IED). Uniform coverage over the parameter space is desirable for general response surface construction as will be discussed later. A 2-D Progressive Lattice Sampling IED sequence is shown in Fig. 1 for variables denoted P1 and P2. In this example, the perimeter boundary of the square represents the bounds of the P1–P2
Fig. 1. 2-D Progressive Lattice Sampling (PLS) incremental experimental design.
204
V.J. Romero et al. / Structural Safety 26 (2004) 201–219
parameter space. Level 1 of the PLS design consists of three samples. One sample is located at the center of the parameter space. Assuming no a priori knowledge of the target function (or joint probability density function if an uncertainty propagation problem) exists, the two other samples are located at arbitrarily picked mid-sides of the perimeter boundary. For n parameters, Level 1 requires n+1 samples. Level 2 adds n samples to complete a 2n+1 face-centered layout. Level 3 adds a 2n factorial design to arrive at a ‘‘Central Composite’’ design. Level 4 adds a Box–Behnken design to complete a 3-level full factorial 3n design. (In 2-D, Levels 3 and 4 have the same layout.) Level 5 adds a sub-scaled 2n full factorial design as shown in the figure. Level 6 adds the appropriate samples to complete a 5n full factorial design. Level 7 adds a sub-scaled 4n full factorial design in the interior of the parameter space as shown. Note that Level 4 requires 9 samples for n=2 parameters, 27 for n=3, and 81 for n=4. Level 7 requires 41 samples for n=2, 189 for n=3, and 881 for n=4. Thus, PLS becomes increasingly expensive as the dimensionality of the parameter space grows. The strength of PLS is that it provides an efficient way to add sample sites so that uniform distribution of the samples over the hypercube parameter space is maintained at the various design levels. Such progressive sampling is one element of an evolving strategy for converging, and assessing convergence of, response surface approximations to the target function being sampled. An assessment of convergence or approximation error is essential in estimating and limiting errors arising from using response surface approximations (RSA) in ‘‘meta procedures’’ such as optimization and uncertainty analysis. Progressing through the first four PLS levels for such convergence assessment can be done for n=5 dimensions with 43 samples, for n=7 with 143 samples, and for n=9 with 531 samples.
3. Compatible interpolation schemes Some compatible fitting/interpolation schemes that can be used with PLS to construct progressively upgradable response surfaces are described here. These are piecewise finite-element interpolation, polynomial regression, and kriging. 3.1. Piecewise finite-element (FE) interpolation In conjunction with Progressive Lattice Sampling, piecewise finite element (FE) interpolation was originally used [5]. The arrangement of samples in each PLS level allows the parameter space to be subdivided into a regular pattern of adjacent polygons, which for two parameter space dimensions results in triangular and quadrilateral 2-D finite elements (FEs). For the 2-D parameter space, a combination of 3-, 4-, and 6-node triangular FEs and 9-node square FEs exist at the various PLS levels. Low-order polynomial interpolation (see below) exists over each element. At the boundary between any two adjacent elements the polynomial functions are automatically matched along the inter-element boundary because the polynomial orders of each element are automatically constrained to be that of their mutual boundary. The order is dictated by the number of nodes or sample points along the boundary: linear for two end nodes, and quadratic when a third node exists at the boundary midpoint. The collection of all the elements together creates a C0-continuous global function over the parameter space that possesses local low-order
V.J. Romero et al. / Structural Safety 26 (2004) 201–219
205
polynomial behavior over each element. As such, the global RSA has considerable freedom to locally conform to the data. Details on the implementation of the FE interpolation technique used here are documented in [5]. In the PLS progression shown in Fig. 1, Levels 1, 2 and 3(4) each have a single finite element and therefore a single polynomial interpolation function over the entire parameter space. These single-polynomial global RSA have the same number of ascending monomials or terms in the polynomial (beginning with the constant term) as the number of sample points over the parameter space, so are sometimes referred to as saturated polynomials. These saturated polynomials replicate the data values at the sample points, whereas the regression polynomials discussed in the next section have fewer terms than the number of samples in the parameter space and therefore generally cannot replicate all of the data exactly, but rather strive to best fit the data in a leastsquares sense. For Level 1, a linear three-term saturated-polynomial RSA is created from three data samples. For Level 2 a ‘‘simple quadratic’’ five-term saturated polynomial (having no twoparameter interaction terms) is created from five data samples. Level 3(4) uses a bi-quadratic saturated polynomial that includes all two-parameter interaction terms plus several higher-order monomials (see next section). For Levels 5 and higher, the parameter space is partitioned into multiple finite elements that together constitute a global interpolation function. This partitioning is denoted with dashed lines in Fig. 1. Level 5 contains four 6-node quadratic triangular elements (first six ascending monomials), and Level 6 contains four 9-node bi-quadratic square elements. Level 7 is a combination of 4-node linear-to-quadratic (see [5]) transition triangular elements at the corners of the parameter space, quadratic triangles at the mid-sides, and five bi-quadratic squares filling the interior. 3.2. Polynomial regression Polynomial regression methods are commonly used to create response surface functions from a set of data samples. Regression is popular since the calculations are simple and the resulting function is a closed-form algebraic expression. For example, a quadratic polynomial has the form Y^ ðpÞ ¼ C0 þ
n X i¼1
Ci pi þ
n X n X Ci ;j pi pj
ð1Þ
i¼1 j¼1
where is Y^ ðpÞ the estimate of the target function at a point in the parameter space having n independent coordinates =(p=(p1, p2, . . ., pn), and the C0, Ci, Ci,j are constant coefficients. To calculate the values of the coefficients, a system of linear equations is formed by applying the above polynomial model at each of K sampling points where the kth point has coordinates p(k) and target function value Y(k). The number of sampling points must be greater than or equal to the number of ascending monomials (terms) in the polynomial. If equal, then the solution of the system of equations determines the coefficients of a saturated polynomial. The saturated polynomial RSA exactly matches the target values at all of the sampling points. When more sampling points than terms in the polynomial are available, then the system of equations is overdetermined
206
V.J. Romero et al. / Structural Safety 26 (2004) 201–219
and a regression procedure is invoked to solve for the coefficients. Because the regression polynomial does not have as many terms as there are data samples to fit, it is over-constrained and cannot in general match the true Y(p) values at the sample points. Instead, it can only be made to fit the data in some ‘‘best’’ sense. Most commonly, the method of least squares [10] is used. For PLS Levels 1, 2, and 3(4), the polynomials used in the regression method are the same saturated polynomials that the FE method uses. For Levels 5 and above, regression is used because, although a global polynomial of sufficiently high order can always be found to pass through all data samples in the parameter space, higher-order polynomials that would enable this for the more densely populated Levels 5–7 would tend to exhibit undesirable oscillations between the data points. The regression approach circumvents this by generally limiting the order of the fitting polynomial to quadratic or cubic order and using regression to minimize global fitting error. In this work the regression polynomial used for Level 5 was a 9-term Lagrange bi-quadratic polynomial having the standard six ascending monomials of a quadratic polynomial in two dimensions [n=2 in Eq. (1)], plus (C7(p1)2p2+C8p1(p2)2+C9(p1)2(p2)2. This was fit to the 13 data samples of Level 5. The functions for Levels 6 and 7 were 16-term bi-cubic Lagrange polynomials (see [6]) which were fit to 25 and 41 data samples, respectively. Note that all of the terms in these polynomial functions were retained in the regression procedure used here. Under some conditions this is not desirable, and methods such as stepwise regression (see [10]) exist to find a quasioptimal subset of monomials. However, these issues were not explored here. 3.3. Kriging Kriging is an interpolation method that originated in the mining and geological spatial statistics community [11]. The premise that underlies kriging is that there is spatial correlation among the response values, Y(p), obtained from the data samples. Because of the spatial correlation assumption, a kriging interpolation response value, Y^ ðpÞ, is biased by the Y(p) response values at ‘‘nearby’’ data points, rather than giving all response values equal weight as is done in polynomial regression. Kriging interpolation is somewhat similar to a weighted least squares approach in traditional statistics, where higher weights are given to response values at nearby data points. The kriging method used in this study produces a C2-continuous interpolating function over the entire parameter space. The kriging method does not assume an underlying global trend (and associated global functional form) as is done in polynomial regression. Thus, it is often useful in approximating arbitrary, potentially multimodal, functions. The form of the kriging function used in this study is Y^ ðpÞ ¼ þ rT R1 ðY eÞ
ð2Þ
where r is a vector of local correlation data; R is a matrix of spatial correlation values among the K data samples; and (Ye) is a vector that expresses the deviation of the vector of responses, Y, from the estimate of the mean, , where e is a K1 vector with all entries equal to one. Both r and R are computed using an exponential correlation function [12], and are dependent on a vector of correlation parameters, ={1, . . ., yn} where n is the dimension of the n-tuple p. The i,jth entry
V.J. Romero et al. / Structural Safety 26 (2004) 201–219
in R is computed as " # n 2 X ðiÞ ðjÞ Ri;j ¼ exp t pt pt
207
ð3Þ
t¼1
for i=1, . . ., K, and j=1, . . ., K. Similarly, the ith entry of the local correlation vector, r, is computed as " # n 2 X ðiÞ ri ¼ exp t pt pt ð4Þ t¼1
for i=1, . . ., K. The unknown parameters in the kriging model are the constant term , and the elements of the correlation vector . The term is computed using a generalized least squares approach and is computed using a maximum likelihood estimation (MLE) procedure. References [12–15] describe the mathematical and statistical methods used to estimate values for and . The matrix algebra and MLE procedure employed in kriging make this approach the most complex and computationally expensive of the three fitting/interpolation schemes examined in this study. However, the computational expense of the fitting/interpolation mathematics in kriging will usually be significantly less than the cost of obtaining even a single response value in engineering applications that involve expensive computer simulations. For this reason, computational expense of kriging versus piecewise finite element interpolation and polynomial regression is not considered to be an issue in this study.
4. Application examples Fig. 2 plots the 2-D target function used in this study. The response value of this multimodal function, Y(p), is defined as:
Fig. 2. Target function to be approximated (‘‘exact function’’).
208
V.J. Romero et al. / Structural Safety 26 (2004) 201–219
r Yðp1; p2Þ ¼ 0:8r þ 0:35sin 2:4 pffiffiffi ½1:5sinð1:3Þ 2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi p2 on the domain 04p141, 04p241. where r ¼ ðp1Þ2 þðp2Þ2 and ¼ arctan p1
ð5Þ
Exact values (samples) of this function are obtained at the indicated lattice points in Fig. 1. The data of the various levels is then fitted and interpolated with piecewise finite element methods, kriging, and regression polynomials. Fig. 3 plots response surfaces for finite element interpolation of PLS Levels 1, 3(4), and 6. This illustrates the convergence of the FE RSA to the target function as samples are added in progressing through various levels of the PLS design. To examine the fitting performance (within the Progressive Lattice Sampling framework) of the various response surface construction methods utilized here, global and local measures of fitting error over the parameter space are investigated next. 4.1. Global fitting performance A simple measure of quality of global fit over the parameter space can be calculated at each PLS level as: 441 P
exacti predictedi
approximate spatial average absolute error ¼ i¼1
441
ð6Þ
where exact and predicted values in the summation come from respective evaluation of the exact function and the particular response surface approximation at 441 equally spaced points on a 2121 square grid overlaid on the domain. This measure is an expedient approximation to the global average integrated absolute error over the domain, which would require a much more involved calculation. The value of the current metric can depend strongly on the grid used if it
Fig. 3. Successive response surface approximations to the exact function by finite-element interpolation of PLS designs: (A)—Level 1 having 3 samples; (B)—Level 3(4) having 9 samples; (C)—Level 6 having 25 samples.
V.J. Romero et al. / Structural Safety 26 (2004) 201–219
209
underresolves the functions. As Figs. 2 and 3 suggest, a 2121 plotting grid appears to be sufficiently dense to achieve adequate representation of the target and approximation functions, though a formal grid sensitivity study was not run. Fig. 4 presents the results from the various interpolation schemes at PLS Levels 1–7. All interpolation methods produce a reduction of global fitting error as the number of Lattice samples increases in progressing to each new level of the PLS IED. Through the first three levels the loworder saturated polynomial approximations perform noticeably better than kriging. With increasing sample density in the higher levels (Levels 5–7) the relative performance of kriging improved quickly, eventually giving the best global fit of any method at Level 7 (where the target function was essentially exactly replicated with kriging). At these higher levels, where polynomial regression and piecewise finite-element (FE) interpolation are used instead of saturated polynomials, kriging performed as well as or better than regression. FE interpolation performed consistently better than regression, and somewhat better than kriging except at Level 7 where kriging performed somewhat better than FE. Global fitting quality as it impacts uncertainty propagation will be investigated next. In a probabilistic uncertainty propagation problem, the system response function maps uncertain (random variable) inputs of the system into an output uncertainty distribution. Concretely, let the values of two uncertain inputs p1 and p2 come from independent normal distributions having means 0.5 and standard deviations =0.5/3. The corresponding joint probability density function (JPDF) of these input random variables is represented in Fig. 5 after truncation of the JPDF beyond the unit p1–p2 parameter space depicted in Fig. 2. The JPDF is then renormalized to integrate to unity over this space. The JPDF defines the relative propensity for attaining specific combinations of input values marked by any given point (location) in the joint parameter space. This likelihood function for attaining given input combinations maps through the applicable input/output functional relationship of the system (such as that shown in Fig. 2) into a corresponding likelihood function for attaining various system output values.
Fig. 4. Level-wise convergence of spatially averaged absolute error.
210
V.J. Romero et al. / Structural Safety 26 (2004) 201–219
Operationally, the JPDF of the input random variables can be mapped through the system input/output response function into an output response distribution via standard Monte Carlo (MC) sampling techniques. Toward this purpose, one million Latin-Hyper-cube- stratified MC samples or (p1,p2) parameter combinations (points in the parameter space) were generated by the LHS [16] code according to the JPDF in Fig. 5. The 2-D target function and the various approximations to it introduced earlier in this section were then evaluated at all 10 6 (p1,p2) parameter sets, yielding a population of 106 response values for each function. The impact of the various response surface approximations on the accuracy of global metrics of the mapped response distribution is shown in Figs. 6 and 7 below. Fig. 6 plots, at each PLS level, the calculated means of the output distributions obtained from using the various RSA construction techniques. Also plotted is the mean of the response distribution obtained from the exact 2-D function. As the number of PLS samples over the parameter space increases, each RSA construction scheme generally yields improving distribution means which approach the exact value. However, no RSA scheme yields unbroken monotonic convergence through the succession of PLS Levels. This is in contrast to the previous global indicator of fitting quality (spatially averaged absolute error), where all RSA methods exhibited monotonic convergence. Fig. 7 plots another global indicator of the impact of RSA fitting error over the entire parameter space. The variance of each mapped response distribution is calculated and then the square root is taken to retain the same physical units as the two preceding global indicators. The kriging approach to RSA construction exhibits large oscillation (the largest of all methods) under the sparse parameter space sampling at initial PLS levels. Ultimately, however, at Level 7 (41 samples) kriging produces the most accurate results of any construction method. The level-wise convergence of the other interpolation schemes is more orderly. The performance of polynomial regression levels off after Level 6 (25 samples), after which the constraints of the bi-cubic polynomial global functional form apparently cannot be significantly overcome by adding more PLS data samples to the regression. Such leveling off is also seen in Figs. 4 and 6.
Fig. 5. Joint Probability Density Function describing the random variables in the problem: normally distributed parameters p1 and p2 with means 0.5, std. deviations =0.167, and truncation of the unit square parameter space at 3 above and below the mean values.
V.J. Romero et al. / Structural Safety 26 (2004) 201–219
211
Fig. 6. Level-wise convergence of mean of mapped response distribution.
Fig. 7. Level-wise convergence of the square root of the variance of the mapped response distribution.
4.2. Local fitting performance Frequently in probabilistic analysis, values of the cumulative distribution function (CDF) of the mapped response distribution are desired. The CDF is simply an integral function of the response distribution (or more properly the response probability density function, PDF) such that the value of the CDF function for any response value r is the integrated probability of attaining a response of magnitude r or less, according to the propensities indicated by the response PDF: ðr PDFðr0 Þdr: ð7Þ CDFðrÞ ¼ 1
212
V.J. Romero et al. / Structural Safety 26 (2004) 201–219
A complementary question is, ‘‘What is the probability that response exceeds some threshold value r?’’ Here the PDF is integrated between the limits of r and 1. For a given threshold value r, this ‘‘exceedence probability’’ EP(r) and the CDF probability CDF(r) are complements: each is 1 minus the other. Either quantity can serve as a mechanism for measuring localized response surface approximation error as explained next. We use the aforementioned MC populations of output response values arrived at using the exact and approximate response functions. For each population, the EP is determined as the ratio of the number of calculated response values at or above some specified threshold value r, to the total number (106) of samples. We then note that the following congruency relationship exists: for a response threshold value of say r=1.5, the MC quotient for exceedence probability, is equal (to within applicable confidence intervals1) to the integral of the JPDF over the exceedence (shaded) region of the parameter space shown in Fig. 8. The figure shows a cutting plane passed through the exact function at a response threshold r=1.5. The shaded region of the threshold plane marks the region of the parameter space where the function exceeds the stipulated threshold value. Accordingly, a multidimensional integration of the JPDF (Fig. 5) over this shaded exceedence region of the parameter space also gives the applicable EP. Through this congruence relationship, we can attribute differences between exact- and RSA-mapped MC results to local RSA errors responsible for altering the shape of the exceedence region over which the JPDF is integrated. We develop this connection further in the following. For a given threshold value r and distribution of (p1, p2) MC sampling points in the parameter space, exceedence probability depends only on the form of the function being sampled. Observing that the numerator in the MC ratio corresponds to the number of sampling points that lie in the shaded region of the parameter space where response exceeds the stipulated threshold value, it is apparent that only the local functional form that governs the threshold-specific contour which separates exceedence and complement regions of the parameter space is important. That is, it is
Fig. 8. Cutting plane through exact function at a response value of 1.5. 1
Actually, MC sampling errors are immaterial here because all functions (exact and RSAs) are evaluated on the same set of 106 MC parameter sets. Any small bias error in the results due to finite MC sampling shows up as a common bias among all the results. Hence, the plots in this paper accurately reveal the relative accuracy of the various RSAs.
V.J. Romero et al. / Structural Safety 26 (2004) 201–219
213
only necessary to accurately approximate the exact function in the vicinity of the threshold contour that demarcates the exceedence and non-exceedence regions of the parameter space. Such local errors between exact and approximate exceedence regions will now be examined for several stipulated threshold levels. Fig. 9 shows several other threshold values to be considered here, and the associated exact exceedence and non-exceedence (‘‘complement’’) regions of the parameter space. Fig. 10 illustrates the impact of response surface approximations on the projected exceedence and complement regions of the parameter space. As the figure shows, for a threshold value of 1.5, a FE interpolation of the Level 1 PLS design (based on 3 samples of the exact function) yields no exceedence region at all. FE response surfaces based on the Level 3(4) 9-sample and Level 6 25sample designs do progressively better in approximating the exact exceedence and complement regions. Figs. 11–13 depict analogous results for threshold values of 1.0, 0.5, and 0.2, respectively. As the figures suggest, it is generally the case that as the RSA are upgraded with more Lattice points, the approximate exceedence and non-exceedence regions approach the exact regions of the target function. These convergence trends will be quantitatively characterized next. Fig. 14 plots results for the 1.5 threshold level. Fairly large initial oscillations in the kriging results are evident as Lattice samples are added to the progressive experimental design. However, eventual convergence to virtually the exact probability occurs by Level 7 (41 PLS samples). The level-wise convergence of the other two interpolation schemes is somewhat more orderly. Although the exact probability here is 0.000252, all methods yield zero probability at Levels 1 and 2 because the underlying designs do not sample the corners of the hypercube. Any significant nonlinear interaction effects in the actual function as the vertices of the hyperspace are approached will not accurately show up in response surfaces based on these designs. The 9-sample Level 3(4) design does sample the corners, but even though the FE interpolation of this design yields a reasonably representative exceedence region (see Fig. 10C), the magnitude of the exceedence probability is so small in this problem that the tiny 0.00015 error in calculated probability represents a 60% deviation from the exact result. Thus, when small exceedence probabilities are involved, calculated results are very sensitive to even small errors in local representation of the
Fig. 9. Cutting planes through exact function showing associated exceedence and complement regions for response threshold values of 1.0, 0.5, and 0.2, respectively.
214
V.J. Romero et al. / Structural Safety 26 (2004) 201–219
Fig. 10. Shaded region of the p1–p2 parameter space where a threshold value of 1.5 is exceeded by: (A), the exact function; and (B)–(D), response surfaces from FE interpolation of Level 1, 3(4), and 6 PLS designs, respectively.
Fig. 11. Shaded region of the p1–p2 parameter space where a threshold value of 1.0 is exceeded by: (A), the exact function; and (B)–(D), response surfaces from FE interpolation of Level 1, 3(4), and 6 PLS designs, respectively.
Fig. 12. Shaded regions of the p1–p2 parameter space where a threshold value of 0.5 is exceeded by: (A), the exact function; and (B)– (D), response surfaces from FE interpolation of Level 1, 3(4), and 6 PLS designs, respectively.
true function. Ultimately, at Level 7 the FE result is about 11% in error, with the cubic regression polynomial yielding a 25% error, and kriging yielding effectively no error. Fig. 15 presents results for the threshold level of 1.0. The trends here are similar to those of Fig. 14, with kriging exhibiting largest initial oscillation but eventually yielding the best (almost exact) result at Level 7. (The remarkable value yielded by kriging at Level 1 is specious because the underlying experimental design is clearly too sparse to justify that degree of accuracy.) The 9sample Level 3(4) design visually appears to begin resolving reasonably well the exceedence region for this threshold (see Fig. 11C), but the FE-derived approximation shown in the figure still results in a large 84% error in cal-culated probability. Again, the magnitude of the exceed-
V.J. Romero et al. / Structural Safety 26 (2004) 201–219
215
Fig. 13. Shaded region of the p1–p2 parameter space where a threshold value of 0.2 is exceeded by: (A), the exact function; and (B)–(D), response surfaces from FE interpolation of Level 1, 3(4), and 6 PLS designs, respectively.
Fig. 14. Level-wise convergence of calculated exceedence probability for a threshold of 1.5.
ence probability is fairly small here (exact result=0.00756), so that the small error magnitude of 0.0064 in the calculated result represents a large relative error. Ultimately, at Level 7 (41 samples) the FE result is only about 4.6% in error, with PR about 11% in error, and negligible error for kriging. Fig. 16 presents results for a threshold value of 0.5. In contrast to the behavior at other threshold levels, here the early convergence behavior of kriging is less erratic than that of loworder saturated polynomials. Still, in absolute terms, kriging exhibits greater error over the first three levels. However, the apparent good performance of saturated polynomials over Levels 1–3 is misleading. Fig. 12B shows the saturated polynomial exceedence region at Level 1 having 3 PLS samples. Although this region is completely different from the exact one in Fig. 12A, only a 3% error exists in the calculated probability. The 9-sample Level 3 result corresponding to Fig. 12C yields only a 6.2% error, yet the approximated exceedence region is still qualitatively very different from the exact one. Therefore, the apparent high relative accuracy of these results is purely coincidental. One tip-off to this fact is that results actually get significantly worse in progressing through the first three levels, suggesting that something is amiss. However, in a real application where the exact result is not known, these same results will not necessarily send a
216
V.J. Romero et al. / Structural Safety 26 (2004) 201–219
Fig. 15. Level-wise convergence of calculated exceedence probability for a threshold of 1.0.
Fig. 16. Level-wise convergence of calculated exceedence probability for a threshold of 0.5.
warning signal: it might appear from looking at the kriging and polynomial interpolations through the first three levels that results are converging to an apparent asymptotic probability of about 0.48, and that the incremental cost of going to the next level of the PLS design can be foregone because the RSA has already apparently closely converged to the target function. Thus, it is difficult in practice to determine when an experimental design and accompanying response surface truly resolve the relevant characteristics of a target function. Therefore, it is important to go as far up the progressive design as affordable, at least (apparently from the problems in this paper) until three consecutive PLS levels yield a stable result. Even then there are no guarantees of accuracy—but one can only do the best that one can afford to do.
V.J. Romero et al. / Structural Safety 26 (2004) 201–219
217
The fact that kriging and FE interpolation yield essentially exact results at Level 7 suggests that the underlying 41-sample PLS design adequately captures the complex variation of the exact function, at least in the vicinity of the exact 0.5 contours (see Fig. 12). That the bi-cubic polynomial did comparatively poorly at Levels 6 and 7 when regressed to the 25-sample and 41sample data sets reflects the conformance constraints of the imposed bi-cubic form of the fitting function. These constraints do not allow adequate conformance to the particularly complex exceedence regions shown in Fig. 12(A) for the 0.5 threshold value. Fig. 17 presents the results for the 0.2 threshold, which yields an exact probability of 0.98439. Again, kriging results oscillate the most over the first few PLS levels, but ultimately become the most accurate (essentially exact) at Level 7. By Level 5 (13 samples) each of the interpolation methods does reasonably well, but before this level there is insufficient sampling of the parameter space to resolve the rather complex target region shown in Fig. 13.
5. Summary and conclusions Over the seven metrics studied here (three global and four local indicators of fitting quality), general reduction of error in the calculated results is exhibited by all interpolation methods as more PLS samples are added to the upgradable PLS experimental design. However, while all schemes exhibited increasing level-wise accuracy for the global metric of spatially-averaged fitting error over the parameter space, some unexpected oscillations where exhibited in the level-wise convergence of the mean and variance global metrics. Similar oscillations were observed in calculated local metrics (exceedence probabilities). These oscillations tended to be worse at the first three PLS levels, with kriging results generally oscillating the most under the sparse sampling at these levels. The erratic behavior of kriging at lower PLS levels may be a result of highly inaccurate (due to insufficient sampling of the target function initially) estimates of the correct correlation length needed to accurately represent the target function.
Fig. 17. Level-wise convergence of calculated exceedence probability for a threshold of 0.2.
218
V.J. Romero et al. / Structural Safety 26 (2004) 201–219
Polynomial regression was sometimes superior to FE and kriging, but FE and kriging usually outperformed polynomial regression. FE generally performed better than kriging except at Level 7 where kriging always performed best (converging essentially to exact results in all cases). The performance of regression generally leveled off after Level 5, after which apparently the constraints of the bi-cubic polynomial global functional form could not be significantly overcome by adding more data samples to the regression. Certainly, many more diverse and higher-dimensional functions must be tested before more firm conclusions can be drawn about the relative applicability of the PLS-compatible response surface construction schemes tried here, but this paper demonstrates a framework and process for such studies. By default, users may be initially constrained to using polynomial regression because of the general availability of software packages for this, but general kriging and finiteelement interpolation schemes may also become widely available in the future. Over the seven cases, no method always performed best or worst, and all methods were best at least once and worst at least once. This is perhaps an indication that several interpolation methods should be tried when constructing a response surface, with intimate knowledge of their tendencies, constraints, biases, and idiosyncrasies being used to ultimately select a method. Even though the greatest share of the expense and difficulty in constructing a response surface normally lies in procuring the data points in the experimental design, all of this effort can be for naught if the interpolation of the data is faulty. Therefore, it is absolutely critical to robustly fit and interpolate the data. In addition, the number and location of data samples must be sufficient to adequately resolve the function for the purposes at hand. As seen here, such adequacy can depend strongly on the particular purpose of the interpolation (what global and/or local fitting needs/metrics apply), and on the interpolation scheme being used. It can be very difficult to determine when a particular sampling design and interpolation scheme sufficiently resolve a function, but this must be done if the response surface is to be used as an effective inexpensive replacement for the actual function. Monitoring convergence heuristics of progressive convergent response surface approximations can help in this regard, and PLS appears to be an efficient approach for generating suitable underlying incremental experimental designs.
Acknowledgements This work was supported by the United States Department of Energy under Contract DE-AC04-94AL8500. Sandia is a multiprogram laboratory operated by Sandia Corporation, a Lockheed Martin Company, for the United States Department of Energy. The authors thank the Sandia reviewers and the Structural Safety referees for their contributions to the improvement of this paper. References [1] Unal R, Lepsch RA, Engelund W, Stanley DO. Approximation model building and multidisciplinary design optimization using response surface methods. Paper AIAA-1996-4044-CP presented at the 6th Multi-Disciplinary Analysis and Optimization Symposium, Belleview, WA, 4–6 September 1996.
V.J. Romero et al. / Structural Safety 26 (2004) 201–219
219
[2] Roux WJ, Stander N, Haftka RT. Response surface approximations for structural optimization. International J Numerical Mthds Engr 1998;42:517–34. [3] Venter G, Haftka RT, Starnes Jr JH. Construction of response surface approximation for design optimization. AIAA Journal 1998;36(12):2242–9. [4] Romero VJ. Optimization and non-deterministic analysis with large simulation models: issues and directions. ASME HTD-Vol. 357-2, Proceedings of the 7th AIAA/ASME Joint Thermophysics and Heat Transfer Conference, Albuquerque, NM, 15–18 June, 1998. [5] Romero VJ, Bankston SD. Finite-element/progressive-lattice-sampling response surface methodology and application to benchmark probability quantification problems. Sandia National Laboratories report SAND98-0567 printed March 1998. [6] Strang G, Fix GJ. An analysis of the finite element method. Englewood Cliffs (NJ): Prentice-Hall; 1973. [7] Patterson TNL. The optimum addition of points to quadrature formulae. Math Comp 1968;22:847–56. [8] Genz AC, Malik AA. An imbedded family of fully symmetric numerical integration rules. SIAM J Numer Anal 1983;20(3). [9] Baden SB, Chrisochoides NP, Gannon DB, Norman ML, editors. IMA Volume on Structured Adaptive Mesh Refinement (SAMR) Grid Methods 2000. Springer-Verlag; January. [10] Myers RH, Montgomery DC. Response surface methodology: process and product optimization using designed experiments. NY: Wiley; 1995. [11] Cressie N. Statistics for spatial data. New York: John Wiley and Sons, Inc; 1991. [12] Koehler JR, Owen AB. Computer experiments. In: Ghosh S, Rao CR, editors. Handbook of statistics, vol. 13. New York: Elsevier Science; 1996. p. 261–308. [13] Giunta AA. Aircraft multidisciplinary design optimization using design of experiments methods and response surface modeling methods. PhD dissertation, Virginia Polytechnic Institute and State University, Blacksburg, VA, 1997. [14] Giunta AA, Watson LT. A comparison of approximation modeling techniques: polynomial versus interpolating models. Paper AIAA-1998-4758 presented at the 7th AIAA/USAF/NASA/ISSMO Symposium on Multidisciplinary Analysis and Optimization, St. Louis, MO, 2–4 September 1998. [15] Sacks J, Welch WJ, Mitchell TJ, Wynn HP. Design and analysis of computer experiments. Statistical Science 1989;4(4):409–35. [16] Iman RL, Shortencarier MJ. A FORTRAN77 program and user’s guide for the generation of latin hypercube and random samples to use with computer models. Sandia National Laboratories report SAND83-2365 (RG), printed March 1984.
Construction of response surfaces based on progressivelattice-sampling experimental designs with application to uncertainty propagation V.J. Romero*, L.P. Swiler, A.A. Giunta Sandia National Laboratories, Albuquerque, NM 87185-0828, USA
Abstract Response surface functions are often used as simple and inexpensive replacements for computationally expensive computer models that simulate the behavior of a complex system over some parameter space. Here we examine several data fitting and interpolation techniques (finite-element interpolation, kriging, and polynomial regression) that can be used to construct a sequence of progressively upgraded response surface approximations based on Progressive Lattice Sampling (PLS) incremental experimental designs. PLS is a paradigm for structured sampling of a hypercube parameter space by placing and incrementally adding samples at each level of the design in a manner intended to efficiently leverage the samples at all previous levels. When combined with compatible interpolation methods, PLS can be used to construct efficiently upgradable response surface approximations. Upon upgrading, convergence heuristics are gained that can be used to estimate the magnitude of the approximation error entailed in the response surface. The three interpolation methods tried here are examined for fitting and convergence behavior in several test problems. Published by Elsevier Ltd. Keywords: Response surfaces; Experimental design; Fitting and interpolation; Uncertainty analysis
1. Introduction Response surface functions are often used as simple and inexpensive replacements of complex and computationally expensive computer models that simulate the behavior of a system over some parameter space. As such, these ‘‘surrogate’’ response surface models fit and interpolate data that is generated by running the simulation model at selected points in the parameter space, e.g. based on some experimental design. The response surface is then used to predict system response over the parameter space at reduced cost, say for uncertainty analysis where the * Corresponding author. Tel.: +1-505-844-5890; fax: +1-505-844-8251. E-mail address: [email protected] (V.J. Romero). 0167-4730/$ - see front matter Published by Elsevier Ltd. doi:10.1016/j.strusafe.2003.03.001
202
V.J. Romero et al. / Structural Safety 26 (2004) 201–219
response surface may be evaluated hundreds or thousands of times. Besides expense, problems associated with numerical noise in the underlying simulation models, and model coupling in multidisciplinary analysis also motivate the use of response surfaces. Examples of response surface usage to facilitate large-scale optimization and uncertainty analysis are cited in references [1–4]. Two issues that arise when using response surface functions are accuracy and the cost of procuring the data samples needed to create the functions. With a sufficiently compliant global fitting/interpolating function over the parameter space, response surface accuracy generally increases as the number of data points increases (if the points are appropriately placed throughout the parameter space), until the essential character of the function is effectively mapped out. Thereafter, it is not cost effective to continue adding samples. Since a single high-fidelity physics simulation (i.e., one data sample) can take hours to compute, it is desirable to minimize the number of simulations that are needed to construct an accurate response surface. A method is presented in [5] to address the accuracy and computational cost associated with using response surface functions under the following circumstances: (1) the physics simulation model is relatively expensive to evaluate; (2) the parameter space is a hypercube or can be accurately and inexpensively mapped into one; (3) the sampled or ‘‘target’’ function is a continuous, deterministic function over the parameter space; (4) reasonably general, arbitrary target functions are to be fitted; and (5) approximate response values are desired over the entire parameter space—i.e., for uncertainty propagation problems or global optimization. In the Finite-Element/Progressive-Lattice-Sampling method, the response surface is a continuous, piecewise-smooth global interpolation function comprised of a patchwork of linear to quadratic polynomials over adjoining subspaces, in the manner of ‘‘finite element’’ representation of spatial functions or ‘‘fields’’ in computational continuum mechanics. The data samples fit and interpolated by the finite-element (FE) response surfaces are arranged on a structured lattice. The so-called ‘‘Progressive Lattice Sampling’’ (PLS) method of generating the lattice points may be viewed as an ‘‘incremental’’ experimental design method. The PLS design attempts to place and add sets of samples such that all samples are efficiently leveraged as the design progresses from one level to the next. Hence, the FE/PLS combination can be used to construct efficiently upgradable response surfaces. A mathematical analysis of finite element interpolation [6] shows that for a continuous and infinitely differentiable function over the parameter space, the domain integral of the pointwise absolute error goes to zero as the spacing between samples goes to zero. This analysis can be applied to the FE/PLS method. Hence, in the limit, FE/PLS response surfaces converge everywhere in the domain to target functions in this class (though the convergence rate is generally not uniform over the domain). Upon upgrading a response surface to the next level of the PLS design, the resulting change in predicted response at any point in the domain is a heuristic indicator of the magnitude of the local error in the response surface approximation. When the incremental change goes to zero, this indicates that the local approximation error has become negligible. The purpose of the current study is to compare the performance of the FE interpolation method against two other response surface construction methods (polynomial regression and kriging) that are compatible with upgradable PLS experimental designs. The convergence behaviors of the three interpolation schemes are examined on several 2-D test problems as the number of data samples is increased in passing through the various levels of the PLS design.
V.J. Romero et al. / Structural Safety 26 (2004) 201–219
203
2. Progressive lattice sampling (PLS) experimental designs As presented here, the implementation of Progressive Lattice Sampling has changed somewhat since it was first introduced in [5], but the essence is still the same. The basic concept (which is also utilized in e.g., upgradable quadrature [7,8] and adaptive nodalization for error reduction in the solution of partial differential equations, [9]) is to progressively sample the parameter space in a manner that efficiently builds on previous information (i.e., fully leverages previous samples). Under the present circumstances this means placing and adding samples in a manner that preserves uniformity of sampling coverage over the hypercube parameter space in the various stages or levels of the incremental experimental design (IED). Uniform coverage over the parameter space is desirable for general response surface construction as will be discussed later. A 2-D Progressive Lattice Sampling IED sequence is shown in Fig. 1 for variables denoted P1 and P2. In this example, the perimeter boundary of the square represents the bounds of the P1–P2
Fig. 1. 2-D Progressive Lattice Sampling (PLS) incremental experimental design.
204
V.J. Romero et al. / Structural Safety 26 (2004) 201–219
parameter space. Level 1 of the PLS design consists of three samples. One sample is located at the center of the parameter space. Assuming no a priori knowledge of the target function (or joint probability density function if an uncertainty propagation problem) exists, the two other samples are located at arbitrarily picked mid-sides of the perimeter boundary. For n parameters, Level 1 requires n+1 samples. Level 2 adds n samples to complete a 2n+1 face-centered layout. Level 3 adds a 2n factorial design to arrive at a ‘‘Central Composite’’ design. Level 4 adds a Box–Behnken design to complete a 3-level full factorial 3n design. (In 2-D, Levels 3 and 4 have the same layout.) Level 5 adds a sub-scaled 2n full factorial design as shown in the figure. Level 6 adds the appropriate samples to complete a 5n full factorial design. Level 7 adds a sub-scaled 4n full factorial design in the interior of the parameter space as shown. Note that Level 4 requires 9 samples for n=2 parameters, 27 for n=3, and 81 for n=4. Level 7 requires 41 samples for n=2, 189 for n=3, and 881 for n=4. Thus, PLS becomes increasingly expensive as the dimensionality of the parameter space grows. The strength of PLS is that it provides an efficient way to add sample sites so that uniform distribution of the samples over the hypercube parameter space is maintained at the various design levels. Such progressive sampling is one element of an evolving strategy for converging, and assessing convergence of, response surface approximations to the target function being sampled. An assessment of convergence or approximation error is essential in estimating and limiting errors arising from using response surface approximations (RSA) in ‘‘meta procedures’’ such as optimization and uncertainty analysis. Progressing through the first four PLS levels for such convergence assessment can be done for n=5 dimensions with 43 samples, for n=7 with 143 samples, and for n=9 with 531 samples.
3. Compatible interpolation schemes Some compatible fitting/interpolation schemes that can be used with PLS to construct progressively upgradable response surfaces are described here. These are piecewise finite-element interpolation, polynomial regression, and kriging. 3.1. Piecewise finite-element (FE) interpolation In conjunction with Progressive Lattice Sampling, piecewise finite element (FE) interpolation was originally used [5]. The arrangement of samples in each PLS level allows the parameter space to be subdivided into a regular pattern of adjacent polygons, which for two parameter space dimensions results in triangular and quadrilateral 2-D finite elements (FEs). For the 2-D parameter space, a combination of 3-, 4-, and 6-node triangular FEs and 9-node square FEs exist at the various PLS levels. Low-order polynomial interpolation (see below) exists over each element. At the boundary between any two adjacent elements the polynomial functions are automatically matched along the inter-element boundary because the polynomial orders of each element are automatically constrained to be that of their mutual boundary. The order is dictated by the number of nodes or sample points along the boundary: linear for two end nodes, and quadratic when a third node exists at the boundary midpoint. The collection of all the elements together creates a C0-continuous global function over the parameter space that possesses local low-order
V.J. Romero et al. / Structural Safety 26 (2004) 201–219
205
polynomial behavior over each element. As such, the global RSA has considerable freedom to locally conform to the data. Details on the implementation of the FE interpolation technique used here are documented in [5]. In the PLS progression shown in Fig. 1, Levels 1, 2 and 3(4) each have a single finite element and therefore a single polynomial interpolation function over the entire parameter space. These single-polynomial global RSA have the same number of ascending monomials or terms in the polynomial (beginning with the constant term) as the number of sample points over the parameter space, so are sometimes referred to as saturated polynomials. These saturated polynomials replicate the data values at the sample points, whereas the regression polynomials discussed in the next section have fewer terms than the number of samples in the parameter space and therefore generally cannot replicate all of the data exactly, but rather strive to best fit the data in a leastsquares sense. For Level 1, a linear three-term saturated-polynomial RSA is created from three data samples. For Level 2 a ‘‘simple quadratic’’ five-term saturated polynomial (having no twoparameter interaction terms) is created from five data samples. Level 3(4) uses a bi-quadratic saturated polynomial that includes all two-parameter interaction terms plus several higher-order monomials (see next section). For Levels 5 and higher, the parameter space is partitioned into multiple finite elements that together constitute a global interpolation function. This partitioning is denoted with dashed lines in Fig. 1. Level 5 contains four 6-node quadratic triangular elements (first six ascending monomials), and Level 6 contains four 9-node bi-quadratic square elements. Level 7 is a combination of 4-node linear-to-quadratic (see [5]) transition triangular elements at the corners of the parameter space, quadratic triangles at the mid-sides, and five bi-quadratic squares filling the interior. 3.2. Polynomial regression Polynomial regression methods are commonly used to create response surface functions from a set of data samples. Regression is popular since the calculations are simple and the resulting function is a closed-form algebraic expression. For example, a quadratic polynomial has the form Y^ ðpÞ ¼ C0 þ
n X i¼1
Ci pi þ
n X n X Ci ;j pi pj
ð1Þ
i¼1 j¼1
where is Y^ ðpÞ the estimate of the target function at a point in the parameter space having n independent coordinates =(p=(p1, p2, . . ., pn), and the C0, Ci, Ci,j are constant coefficients. To calculate the values of the coefficients, a system of linear equations is formed by applying the above polynomial model at each of K sampling points where the kth point has coordinates p(k) and target function value Y(k). The number of sampling points must be greater than or equal to the number of ascending monomials (terms) in the polynomial. If equal, then the solution of the system of equations determines the coefficients of a saturated polynomial. The saturated polynomial RSA exactly matches the target values at all of the sampling points. When more sampling points than terms in the polynomial are available, then the system of equations is overdetermined
206
V.J. Romero et al. / Structural Safety 26 (2004) 201–219
and a regression procedure is invoked to solve for the coefficients. Because the regression polynomial does not have as many terms as there are data samples to fit, it is over-constrained and cannot in general match the true Y(p) values at the sample points. Instead, it can only be made to fit the data in some ‘‘best’’ sense. Most commonly, the method of least squares [10] is used. For PLS Levels 1, 2, and 3(4), the polynomials used in the regression method are the same saturated polynomials that the FE method uses. For Levels 5 and above, regression is used because, although a global polynomial of sufficiently high order can always be found to pass through all data samples in the parameter space, higher-order polynomials that would enable this for the more densely populated Levels 5–7 would tend to exhibit undesirable oscillations between the data points. The regression approach circumvents this by generally limiting the order of the fitting polynomial to quadratic or cubic order and using regression to minimize global fitting error. In this work the regression polynomial used for Level 5 was a 9-term Lagrange bi-quadratic polynomial having the standard six ascending monomials of a quadratic polynomial in two dimensions [n=2 in Eq. (1)], plus (C7(p1)2p2+C8p1(p2)2+C9(p1)2(p2)2. This was fit to the 13 data samples of Level 5. The functions for Levels 6 and 7 were 16-term bi-cubic Lagrange polynomials (see [6]) which were fit to 25 and 41 data samples, respectively. Note that all of the terms in these polynomial functions were retained in the regression procedure used here. Under some conditions this is not desirable, and methods such as stepwise regression (see [10]) exist to find a quasioptimal subset of monomials. However, these issues were not explored here. 3.3. Kriging Kriging is an interpolation method that originated in the mining and geological spatial statistics community [11]. The premise that underlies kriging is that there is spatial correlation among the response values, Y(p), obtained from the data samples. Because of the spatial correlation assumption, a kriging interpolation response value, Y^ ðpÞ, is biased by the Y(p) response values at ‘‘nearby’’ data points, rather than giving all response values equal weight as is done in polynomial regression. Kriging interpolation is somewhat similar to a weighted least squares approach in traditional statistics, where higher weights are given to response values at nearby data points. The kriging method used in this study produces a C2-continuous interpolating function over the entire parameter space. The kriging method does not assume an underlying global trend (and associated global functional form) as is done in polynomial regression. Thus, it is often useful in approximating arbitrary, potentially multimodal, functions. The form of the kriging function used in this study is Y^ ðpÞ ¼ þ rT R1 ðY eÞ
ð2Þ
where r is a vector of local correlation data; R is a matrix of spatial correlation values among the K data samples; and (Ye) is a vector that expresses the deviation of the vector of responses, Y, from the estimate of the mean, , where e is a K1 vector with all entries equal to one. Both r and R are computed using an exponential correlation function [12], and are dependent on a vector of correlation parameters, ={1, . . ., yn} where n is the dimension of the n-tuple p. The i,jth entry
V.J. Romero et al. / Structural Safety 26 (2004) 201–219
in R is computed as " # n 2 X ðiÞ ðjÞ Ri;j ¼ exp t pt pt
207
ð3Þ
t¼1
for i=1, . . ., K, and j=1, . . ., K. Similarly, the ith entry of the local correlation vector, r, is computed as " # n 2 X ðiÞ ri ¼ exp t pt pt ð4Þ t¼1
for i=1, . . ., K. The unknown parameters in the kriging model are the constant term , and the elements of the correlation vector . The term is computed using a generalized least squares approach and is computed using a maximum likelihood estimation (MLE) procedure. References [12–15] describe the mathematical and statistical methods used to estimate values for and . The matrix algebra and MLE procedure employed in kriging make this approach the most complex and computationally expensive of the three fitting/interpolation schemes examined in this study. However, the computational expense of the fitting/interpolation mathematics in kriging will usually be significantly less than the cost of obtaining even a single response value in engineering applications that involve expensive computer simulations. For this reason, computational expense of kriging versus piecewise finite element interpolation and polynomial regression is not considered to be an issue in this study.
4. Application examples Fig. 2 plots the 2-D target function used in this study. The response value of this multimodal function, Y(p), is defined as:
Fig. 2. Target function to be approximated (‘‘exact function’’).
208
V.J. Romero et al. / Structural Safety 26 (2004) 201–219
r Yðp1; p2Þ ¼ 0:8r þ 0:35sin 2:4 pffiffiffi ½1:5sinð1:3Þ 2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi p2 on the domain 04p141, 04p241. where r ¼ ðp1Þ2 þðp2Þ2 and ¼ arctan p1
ð5Þ
Exact values (samples) of this function are obtained at the indicated lattice points in Fig. 1. The data of the various levels is then fitted and interpolated with piecewise finite element methods, kriging, and regression polynomials. Fig. 3 plots response surfaces for finite element interpolation of PLS Levels 1, 3(4), and 6. This illustrates the convergence of the FE RSA to the target function as samples are added in progressing through various levels of the PLS design. To examine the fitting performance (within the Progressive Lattice Sampling framework) of the various response surface construction methods utilized here, global and local measures of fitting error over the parameter space are investigated next. 4.1. Global fitting performance A simple measure of quality of global fit over the parameter space can be calculated at each PLS level as: 441 P
exacti predictedi
approximate spatial average absolute error ¼ i¼1
441
ð6Þ
where exact and predicted values in the summation come from respective evaluation of the exact function and the particular response surface approximation at 441 equally spaced points on a 2121 square grid overlaid on the domain. This measure is an expedient approximation to the global average integrated absolute error over the domain, which would require a much more involved calculation. The value of the current metric can depend strongly on the grid used if it
Fig. 3. Successive response surface approximations to the exact function by finite-element interpolation of PLS designs: (A)—Level 1 having 3 samples; (B)—Level 3(4) having 9 samples; (C)—Level 6 having 25 samples.
V.J. Romero et al. / Structural Safety 26 (2004) 201–219
209
underresolves the functions. As Figs. 2 and 3 suggest, a 2121 plotting grid appears to be sufficiently dense to achieve adequate representation of the target and approximation functions, though a formal grid sensitivity study was not run. Fig. 4 presents the results from the various interpolation schemes at PLS Levels 1–7. All interpolation methods produce a reduction of global fitting error as the number of Lattice samples increases in progressing to each new level of the PLS IED. Through the first three levels the loworder saturated polynomial approximations perform noticeably better than kriging. With increasing sample density in the higher levels (Levels 5–7) the relative performance of kriging improved quickly, eventually giving the best global fit of any method at Level 7 (where the target function was essentially exactly replicated with kriging). At these higher levels, where polynomial regression and piecewise finite-element (FE) interpolation are used instead of saturated polynomials, kriging performed as well as or better than regression. FE interpolation performed consistently better than regression, and somewhat better than kriging except at Level 7 where kriging performed somewhat better than FE. Global fitting quality as it impacts uncertainty propagation will be investigated next. In a probabilistic uncertainty propagation problem, the system response function maps uncertain (random variable) inputs of the system into an output uncertainty distribution. Concretely, let the values of two uncertain inputs p1 and p2 come from independent normal distributions having means 0.5 and standard deviations =0.5/3. The corresponding joint probability density function (JPDF) of these input random variables is represented in Fig. 5 after truncation of the JPDF beyond the unit p1–p2 parameter space depicted in Fig. 2. The JPDF is then renormalized to integrate to unity over this space. The JPDF defines the relative propensity for attaining specific combinations of input values marked by any given point (location) in the joint parameter space. This likelihood function for attaining given input combinations maps through the applicable input/output functional relationship of the system (such as that shown in Fig. 2) into a corresponding likelihood function for attaining various system output values.
Fig. 4. Level-wise convergence of spatially averaged absolute error.
210
V.J. Romero et al. / Structural Safety 26 (2004) 201–219
Operationally, the JPDF of the input random variables can be mapped through the system input/output response function into an output response distribution via standard Monte Carlo (MC) sampling techniques. Toward this purpose, one million Latin-Hyper-cube- stratified MC samples or (p1,p2) parameter combinations (points in the parameter space) were generated by the LHS [16] code according to the JPDF in Fig. 5. The 2-D target function and the various approximations to it introduced earlier in this section were then evaluated at all 10 6 (p1,p2) parameter sets, yielding a population of 106 response values for each function. The impact of the various response surface approximations on the accuracy of global metrics of the mapped response distribution is shown in Figs. 6 and 7 below. Fig. 6 plots, at each PLS level, the calculated means of the output distributions obtained from using the various RSA construction techniques. Also plotted is the mean of the response distribution obtained from the exact 2-D function. As the number of PLS samples over the parameter space increases, each RSA construction scheme generally yields improving distribution means which approach the exact value. However, no RSA scheme yields unbroken monotonic convergence through the succession of PLS Levels. This is in contrast to the previous global indicator of fitting quality (spatially averaged absolute error), where all RSA methods exhibited monotonic convergence. Fig. 7 plots another global indicator of the impact of RSA fitting error over the entire parameter space. The variance of each mapped response distribution is calculated and then the square root is taken to retain the same physical units as the two preceding global indicators. The kriging approach to RSA construction exhibits large oscillation (the largest of all methods) under the sparse parameter space sampling at initial PLS levels. Ultimately, however, at Level 7 (41 samples) kriging produces the most accurate results of any construction method. The level-wise convergence of the other interpolation schemes is more orderly. The performance of polynomial regression levels off after Level 6 (25 samples), after which the constraints of the bi-cubic polynomial global functional form apparently cannot be significantly overcome by adding more PLS data samples to the regression. Such leveling off is also seen in Figs. 4 and 6.
Fig. 5. Joint Probability Density Function describing the random variables in the problem: normally distributed parameters p1 and p2 with means 0.5, std. deviations =0.167, and truncation of the unit square parameter space at 3 above and below the mean values.
V.J. Romero et al. / Structural Safety 26 (2004) 201–219
211
Fig. 6. Level-wise convergence of mean of mapped response distribution.
Fig. 7. Level-wise convergence of the square root of the variance of the mapped response distribution.
4.2. Local fitting performance Frequently in probabilistic analysis, values of the cumulative distribution function (CDF) of the mapped response distribution are desired. The CDF is simply an integral function of the response distribution (or more properly the response probability density function, PDF) such that the value of the CDF function for any response value r is the integrated probability of attaining a response of magnitude r or less, according to the propensities indicated by the response PDF: ðr PDFðr0 Þdr: ð7Þ CDFðrÞ ¼ 1
212
V.J. Romero et al. / Structural Safety 26 (2004) 201–219
A complementary question is, ‘‘What is the probability that response exceeds some threshold value r?’’ Here the PDF is integrated between the limits of r and 1. For a given threshold value r, this ‘‘exceedence probability’’ EP(r) and the CDF probability CDF(r) are complements: each is 1 minus the other. Either quantity can serve as a mechanism for measuring localized response surface approximation error as explained next. We use the aforementioned MC populations of output response values arrived at using the exact and approximate response functions. For each population, the EP is determined as the ratio of the number of calculated response values at or above some specified threshold value r, to the total number (106) of samples. We then note that the following congruency relationship exists: for a response threshold value of say r=1.5, the MC quotient for exceedence probability, is equal (to within applicable confidence intervals1) to the integral of the JPDF over the exceedence (shaded) region of the parameter space shown in Fig. 8. The figure shows a cutting plane passed through the exact function at a response threshold r=1.5. The shaded region of the threshold plane marks the region of the parameter space where the function exceeds the stipulated threshold value. Accordingly, a multidimensional integration of the JPDF (Fig. 5) over this shaded exceedence region of the parameter space also gives the applicable EP. Through this congruence relationship, we can attribute differences between exact- and RSA-mapped MC results to local RSA errors responsible for altering the shape of the exceedence region over which the JPDF is integrated. We develop this connection further in the following. For a given threshold value r and distribution of (p1, p2) MC sampling points in the parameter space, exceedence probability depends only on the form of the function being sampled. Observing that the numerator in the MC ratio corresponds to the number of sampling points that lie in the shaded region of the parameter space where response exceeds the stipulated threshold value, it is apparent that only the local functional form that governs the threshold-specific contour which separates exceedence and complement regions of the parameter space is important. That is, it is
Fig. 8. Cutting plane through exact function at a response value of 1.5. 1
Actually, MC sampling errors are immaterial here because all functions (exact and RSAs) are evaluated on the same set of 106 MC parameter sets. Any small bias error in the results due to finite MC sampling shows up as a common bias among all the results. Hence, the plots in this paper accurately reveal the relative accuracy of the various RSAs.
V.J. Romero et al. / Structural Safety 26 (2004) 201–219
213
only necessary to accurately approximate the exact function in the vicinity of the threshold contour that demarcates the exceedence and non-exceedence regions of the parameter space. Such local errors between exact and approximate exceedence regions will now be examined for several stipulated threshold levels. Fig. 9 shows several other threshold values to be considered here, and the associated exact exceedence and non-exceedence (‘‘complement’’) regions of the parameter space. Fig. 10 illustrates the impact of response surface approximations on the projected exceedence and complement regions of the parameter space. As the figure shows, for a threshold value of 1.5, a FE interpolation of the Level 1 PLS design (based on 3 samples of the exact function) yields no exceedence region at all. FE response surfaces based on the Level 3(4) 9-sample and Level 6 25sample designs do progressively better in approximating the exact exceedence and complement regions. Figs. 11–13 depict analogous results for threshold values of 1.0, 0.5, and 0.2, respectively. As the figures suggest, it is generally the case that as the RSA are upgraded with more Lattice points, the approximate exceedence and non-exceedence regions approach the exact regions of the target function. These convergence trends will be quantitatively characterized next. Fig. 14 plots results for the 1.5 threshold level. Fairly large initial oscillations in the kriging results are evident as Lattice samples are added to the progressive experimental design. However, eventual convergence to virtually the exact probability occurs by Level 7 (41 PLS samples). The level-wise convergence of the other two interpolation schemes is somewhat more orderly. Although the exact probability here is 0.000252, all methods yield zero probability at Levels 1 and 2 because the underlying designs do not sample the corners of the hypercube. Any significant nonlinear interaction effects in the actual function as the vertices of the hyperspace are approached will not accurately show up in response surfaces based on these designs. The 9-sample Level 3(4) design does sample the corners, but even though the FE interpolation of this design yields a reasonably representative exceedence region (see Fig. 10C), the magnitude of the exceedence probability is so small in this problem that the tiny 0.00015 error in calculated probability represents a 60% deviation from the exact result. Thus, when small exceedence probabilities are involved, calculated results are very sensitive to even small errors in local representation of the
Fig. 9. Cutting planes through exact function showing associated exceedence and complement regions for response threshold values of 1.0, 0.5, and 0.2, respectively.
214
V.J. Romero et al. / Structural Safety 26 (2004) 201–219
Fig. 10. Shaded region of the p1–p2 parameter space where a threshold value of 1.5 is exceeded by: (A), the exact function; and (B)–(D), response surfaces from FE interpolation of Level 1, 3(4), and 6 PLS designs, respectively.
Fig. 11. Shaded region of the p1–p2 parameter space where a threshold value of 1.0 is exceeded by: (A), the exact function; and (B)–(D), response surfaces from FE interpolation of Level 1, 3(4), and 6 PLS designs, respectively.
Fig. 12. Shaded regions of the p1–p2 parameter space where a threshold value of 0.5 is exceeded by: (A), the exact function; and (B)– (D), response surfaces from FE interpolation of Level 1, 3(4), and 6 PLS designs, respectively.
true function. Ultimately, at Level 7 the FE result is about 11% in error, with the cubic regression polynomial yielding a 25% error, and kriging yielding effectively no error. Fig. 15 presents results for the threshold level of 1.0. The trends here are similar to those of Fig. 14, with kriging exhibiting largest initial oscillation but eventually yielding the best (almost exact) result at Level 7. (The remarkable value yielded by kriging at Level 1 is specious because the underlying experimental design is clearly too sparse to justify that degree of accuracy.) The 9sample Level 3(4) design visually appears to begin resolving reasonably well the exceedence region for this threshold (see Fig. 11C), but the FE-derived approximation shown in the figure still results in a large 84% error in cal-culated probability. Again, the magnitude of the exceed-
V.J. Romero et al. / Structural Safety 26 (2004) 201–219
215
Fig. 13. Shaded region of the p1–p2 parameter space where a threshold value of 0.2 is exceeded by: (A), the exact function; and (B)–(D), response surfaces from FE interpolation of Level 1, 3(4), and 6 PLS designs, respectively.
Fig. 14. Level-wise convergence of calculated exceedence probability for a threshold of 1.5.
ence probability is fairly small here (exact result=0.00756), so that the small error magnitude of 0.0064 in the calculated result represents a large relative error. Ultimately, at Level 7 (41 samples) the FE result is only about 4.6% in error, with PR about 11% in error, and negligible error for kriging. Fig. 16 presents results for a threshold value of 0.5. In contrast to the behavior at other threshold levels, here the early convergence behavior of kriging is less erratic than that of loworder saturated polynomials. Still, in absolute terms, kriging exhibits greater error over the first three levels. However, the apparent good performance of saturated polynomials over Levels 1–3 is misleading. Fig. 12B shows the saturated polynomial exceedence region at Level 1 having 3 PLS samples. Although this region is completely different from the exact one in Fig. 12A, only a 3% error exists in the calculated probability. The 9-sample Level 3 result corresponding to Fig. 12C yields only a 6.2% error, yet the approximated exceedence region is still qualitatively very different from the exact one. Therefore, the apparent high relative accuracy of these results is purely coincidental. One tip-off to this fact is that results actually get significantly worse in progressing through the first three levels, suggesting that something is amiss. However, in a real application where the exact result is not known, these same results will not necessarily send a
216
V.J. Romero et al. / Structural Safety 26 (2004) 201–219
Fig. 15. Level-wise convergence of calculated exceedence probability for a threshold of 1.0.
Fig. 16. Level-wise convergence of calculated exceedence probability for a threshold of 0.5.
warning signal: it might appear from looking at the kriging and polynomial interpolations through the first three levels that results are converging to an apparent asymptotic probability of about 0.48, and that the incremental cost of going to the next level of the PLS design can be foregone because the RSA has already apparently closely converged to the target function. Thus, it is difficult in practice to determine when an experimental design and accompanying response surface truly resolve the relevant characteristics of a target function. Therefore, it is important to go as far up the progressive design as affordable, at least (apparently from the problems in this paper) until three consecutive PLS levels yield a stable result. Even then there are no guarantees of accuracy—but one can only do the best that one can afford to do.
V.J. Romero et al. / Structural Safety 26 (2004) 201–219
217
The fact that kriging and FE interpolation yield essentially exact results at Level 7 suggests that the underlying 41-sample PLS design adequately captures the complex variation of the exact function, at least in the vicinity of the exact 0.5 contours (see Fig. 12). That the bi-cubic polynomial did comparatively poorly at Levels 6 and 7 when regressed to the 25-sample and 41sample data sets reflects the conformance constraints of the imposed bi-cubic form of the fitting function. These constraints do not allow adequate conformance to the particularly complex exceedence regions shown in Fig. 12(A) for the 0.5 threshold value. Fig. 17 presents the results for the 0.2 threshold, which yields an exact probability of 0.98439. Again, kriging results oscillate the most over the first few PLS levels, but ultimately become the most accurate (essentially exact) at Level 7. By Level 5 (13 samples) each of the interpolation methods does reasonably well, but before this level there is insufficient sampling of the parameter space to resolve the rather complex target region shown in Fig. 13.
5. Summary and conclusions Over the seven metrics studied here (three global and four local indicators of fitting quality), general reduction of error in the calculated results is exhibited by all interpolation methods as more PLS samples are added to the upgradable PLS experimental design. However, while all schemes exhibited increasing level-wise accuracy for the global metric of spatially-averaged fitting error over the parameter space, some unexpected oscillations where exhibited in the level-wise convergence of the mean and variance global metrics. Similar oscillations were observed in calculated local metrics (exceedence probabilities). These oscillations tended to be worse at the first three PLS levels, with kriging results generally oscillating the most under the sparse sampling at these levels. The erratic behavior of kriging at lower PLS levels may be a result of highly inaccurate (due to insufficient sampling of the target function initially) estimates of the correct correlation length needed to accurately represent the target function.
Fig. 17. Level-wise convergence of calculated exceedence probability for a threshold of 0.2.
218
V.J. Romero et al. / Structural Safety 26 (2004) 201–219
Polynomial regression was sometimes superior to FE and kriging, but FE and kriging usually outperformed polynomial regression. FE generally performed better than kriging except at Level 7 where kriging always performed best (converging essentially to exact results in all cases). The performance of regression generally leveled off after Level 5, after which apparently the constraints of the bi-cubic polynomial global functional form could not be significantly overcome by adding more data samples to the regression. Certainly, many more diverse and higher-dimensional functions must be tested before more firm conclusions can be drawn about the relative applicability of the PLS-compatible response surface construction schemes tried here, but this paper demonstrates a framework and process for such studies. By default, users may be initially constrained to using polynomial regression because of the general availability of software packages for this, but general kriging and finiteelement interpolation schemes may also become widely available in the future. Over the seven cases, no method always performed best or worst, and all methods were best at least once and worst at least once. This is perhaps an indication that several interpolation methods should be tried when constructing a response surface, with intimate knowledge of their tendencies, constraints, biases, and idiosyncrasies being used to ultimately select a method. Even though the greatest share of the expense and difficulty in constructing a response surface normally lies in procuring the data points in the experimental design, all of this effort can be for naught if the interpolation of the data is faulty. Therefore, it is absolutely critical to robustly fit and interpolate the data. In addition, the number and location of data samples must be sufficient to adequately resolve the function for the purposes at hand. As seen here, such adequacy can depend strongly on the particular purpose of the interpolation (what global and/or local fitting needs/metrics apply), and on the interpolation scheme being used. It can be very difficult to determine when a particular sampling design and interpolation scheme sufficiently resolve a function, but this must be done if the response surface is to be used as an effective inexpensive replacement for the actual function. Monitoring convergence heuristics of progressive convergent response surface approximations can help in this regard, and PLS appears to be an efficient approach for generating suitable underlying incremental experimental designs.
Acknowledgements This work was supported by the United States Department of Energy under Contract DE-AC04-94AL8500. Sandia is a multiprogram laboratory operated by Sandia Corporation, a Lockheed Martin Company, for the United States Department of Energy. The authors thank the Sandia reviewers and the Structural Safety referees for their contributions to the improvement of this paper. References [1] Unal R, Lepsch RA, Engelund W, Stanley DO. Approximation model building and multidisciplinary design optimization using response surface methods. Paper AIAA-1996-4044-CP presented at the 6th Multi-Disciplinary Analysis and Optimization Symposium, Belleview, WA, 4–6 September 1996.
V.J. Romero et al. / Structural Safety 26 (2004) 201–219
219
[2] Roux WJ, Stander N, Haftka RT. Response surface approximations for structural optimization. International J Numerical Mthds Engr 1998;42:517–34. [3] Venter G, Haftka RT, Starnes Jr JH. Construction of response surface approximation for design optimization. AIAA Journal 1998;36(12):2242–9. [4] Romero VJ. Optimization and non-deterministic analysis with large simulation models: issues and directions. ASME HTD-Vol. 357-2, Proceedings of the 7th AIAA/ASME Joint Thermophysics and Heat Transfer Conference, Albuquerque, NM, 15–18 June, 1998. [5] Romero VJ, Bankston SD. Finite-element/progressive-lattice-sampling response surface methodology and application to benchmark probability quantification problems. Sandia National Laboratories report SAND98-0567 printed March 1998. [6] Strang G, Fix GJ. An analysis of the finite element method. Englewood Cliffs (NJ): Prentice-Hall; 1973. [7] Patterson TNL. The optimum addition of points to quadrature formulae. Math Comp 1968;22:847–56. [8] Genz AC, Malik AA. An imbedded family of fully symmetric numerical integration rules. SIAM J Numer Anal 1983;20(3). [9] Baden SB, Chrisochoides NP, Gannon DB, Norman ML, editors. IMA Volume on Structured Adaptive Mesh Refinement (SAMR) Grid Methods 2000. Springer-Verlag; January. [10] Myers RH, Montgomery DC. Response surface methodology: process and product optimization using designed experiments. NY: Wiley; 1995. [11] Cressie N. Statistics for spatial data. New York: John Wiley and Sons, Inc; 1991. [12] Koehler JR, Owen AB. Computer experiments. In: Ghosh S, Rao CR, editors. Handbook of statistics, vol. 13. New York: Elsevier Science; 1996. p. 261–308. [13] Giunta AA. Aircraft multidisciplinary design optimization using design of experiments methods and response surface modeling methods. PhD dissertation, Virginia Polytechnic Institute and State University, Blacksburg, VA, 1997. [14] Giunta AA, Watson LT. A comparison of approximation modeling techniques: polynomial versus interpolating models. Paper AIAA-1998-4758 presented at the 7th AIAA/USAF/NASA/ISSMO Symposium on Multidisciplinary Analysis and Optimization, St. Louis, MO, 2–4 September 1998. [15] Sacks J, Welch WJ, Mitchell TJ, Wynn HP. Design and analysis of computer experiments. Statistical Science 1989;4(4):409–35. [16] Iman RL, Shortencarier MJ. A FORTRAN77 program and user’s guide for the generation of latin hypercube and random samples to use with computer models. Sandia National Laboratories report SAND83-2365 (RG), printed March 1984.