Multiparametric response surface construction by means of proper generalized decomposition: An extension of the PARAFAC procedure

Comput. Methods Appl. Mech. Engrg. 253 (2013) 543–557

Contents lists available at SciVerse ScienceDirect

Comput. Methods Appl. Mech. Engrg. journal homepage: www.elsevier.com/locate/cma

Multiparametric response surface construction by means of proper generalized decomposition: An extension of the PARAFAC procedure F. El Halabi a,b, D. González a, Agustí Chico-Roca c, M. Doblaré a,b,d,⇑ a

Group of Structural Mechanics and Materials Modelling (GEMM), Aragón Institute of Engineering Research (13A), Universidad de Zaragoza, Spain Centro de Investigación Biomédica en Red en Bioingeniería, Biomateriales y Nanomedicina (CIBER-BBN), Spain c Ascamm Technology Center, Ascamm Private Foundation, Barcelona, Spain d Abengoa Research, Seville, Spain b

a r t i c l e

i n f o

Article history: Received 4 October 2011 Received in revised form 7 August 2012 Accepted 8 August 2012 Available online 18 August 2012 Keywords: Proper generalized decomposition (PGD) Multidimensional least-square Response surface methodology (RSM)

a b s t r a c t Multidimensional problems are found in almost every scientific field. In particular, this is standard in parametric design, inverse analysis, in optimization and in metamodeling analysis, whereby statistical or deterministic approximations of multiparametric solutions are built from the results of experimental campaigns or computer simulations. Multidimensional fitting or approximation of response functions exponentially increase their complexity and computational cost with the number of dimensions responding to the well-known ‘‘curse of dimensionality’’. To reduce the order of complexity and make the solution of many-parameter problems affordable, we propose to combine the model reduction technique known as Proper Generalized Decomposition (PGD) and the response surface (RSM) methodology. As a proof of concept we have used a simple fitting procedure as it is least squares, although other more complex fitting procedures may be easily included. The combined algorithm is presented and its capabilities discussed in a set of multidimensional examples. The data samples to be fit in each of these examples are obtained by means of appropriate discretizing the interval of interest for each design factor and then generating the output values (exact or stochastically modified) by means of virtual experiments. To have an idea of the number of discretization points needed along each direction, the Taguchi’s design of experiments is used. The obtained results show evident improvements in computer time and accuracy when compared to other traditional multiparametric approximation techniques based on polynomial functions and the standard Levenberg–Marquardt algorithm, especially in problems with non-linear behavior and with high number of design parameters. Further comparison was done with the PARAFAC-ALS algorithm. The combination of PGD and RSM seems to be an appealing tool for accurately modelling multiparametric problems in almost real-time if a sufficient set of previous off-line results is available despite the intrinsic complexity of the problem and of the number of parameters involved. Ó 2012 Elsevier B.V. All rights reserved.

1. Introduction Many systems are modelled in such a way that an output (dependent) variable depends upon a big number of input (independent) variables. The usual approach to analyze this type of problems has been determining the effect of the independent variable on each output one, once at a time, while holding the rest of independent variables constant or, alternatively, to perform a set of physical or virtual experiments under specific combinations of the independent variables for getting an overall picture of the combined behavior of the output within the range of interest of

⇑ Corresponding author at: Group of Structural Mechanics and Materials Modelling (GEMM), Aragón Institute of Engineering Research (13A), Universidad de Zaragoza, Spain. Tel.: +34 976 761000; fax: +34 976 762578. E-mail address: [email protected] (M. Doblaré). 0045-7825/$ - see front matter Ó 2012 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.cma.2012.08.005

the input parameters. These two techniques are known as the ’’univariate’’ and ’’factorial’’ methods, respectively [1]. This situation is typical in parametric design, optimization, inverse analysis and statistical problems where the high number of initial parameters are usually reduced, fixing some of them or discarding their contribution, in order to make the problem manageable. In addition, the cost of getting a simple value of the output variable (mostly a time or spatial distribution) may be very expensive, either if they are derived from demanding experimental tests or from highly complex mathematical models [2]. A large variety of techniques has been implemented to deal with this kind of problems with the aim of accelerating the obtention of the dependent variable for specific combinations of independent factors by simplifying the model (model reduction approaches [3–6]), by approximating the overall solution from a limited set of input–output relations (functional fitting [7], Kriging [8], Response Surface approaches [9–11,3,12], design of

544

F. El Halabi et al. / Comput. Methods Appl. Mech. Engrg. 253 (2013) 543–557

experiments [13–16]) or, finally, by using heuristic metamodels (neural networks [4,5], inductive learning [17]). All of them have their respective advantages and drawbacks but in any of them we have to face the well-known ’’curse of dimensionality’’ [18] since the complexity and cost exponentially increase with the number of dimensions, making this problem impossible when the number of dimensions is over a certain number that depends on the computer capacity. Several strategies have been proposed to solve this problem. In particular, we use here a recent method called Proper Generalized Decomposition (PGD) [19–21] that approximates a multidimensional variable in terms of a finite sum of products of separated functions of each independent variable following the well-known Fourier scheme. As a proof of concept of the possibilities of this technique in metamodel approaches and, in particular, in approximating variables depending on many independent parameters, we discuss here the combination of the classical RSM [3] formulation with PGD, resulting into a methodology that will be here referred as RSPGD, although coupling PGD with any other approximation technique would be straightforward. Former works have been presented based in a separated representation of the approximating function to solve fitting problems, as shown in [22], where Proper Orthogonal Decomposition (POD) was used [23]. Both, PGD and POD are model reduction techniques based on separate representation of the variable, however the main difference lies in the way the reduced basis functions are defined and constructed. The POD is classically used as an ’’a posteriori’’ model reduction technique for long-time simulations. In this work we focus in a generalization of POD, as it is PGD that is considered an ’’a priori’’ construction of such separated representation [24]. Another separated representation technique is the so called Parallel Factor (PARAFAC) analysis. This procedure is an extension of low-rank matrix decomposition to higher way arrays [25]. It decomposes a given array in a sum of multilinear terms, analogous to the familiar bilinear vector outer products that appear in matrix decomposition. PARAFAC is applied in the analysis of data arrays indexed by three or more independent variables, just like singular value decomposition (SVD) is applied in ordinary matrix (two-way array) analysis. Unlike SVD, PARAFAC does not impose orthogonality constraints [26]. This method, in its standard form can not deal with sparse cloud of points, thus, a full factorial data distribution over mesh discretization is needed. A multitude of algorithms have been developed to fit a PARAFAC model [27,25], i.e., GRAM-DTLD, PARAFAC-ALS, ASD, SWATLD, PMF3 and dGN. In the case of a fitting problem like the one treated in this work, and once the problem has been discretized, the resulting problem is formally similar to PARAFAC, since it uses a sum of tensorial decompositions of the function in terms of parameter-associated vectors. However, some differences in concept and physical meaning can be mentioned: (i) it approximates a whole multidimensional function in terms of products of continuous functions which is not the objective of PARAFAC; (ii) in the proposed procedure, the functions of each approximation term are piecewise linear defined in all the parameter range and not only at nodal coordinates like in standard PARAFAC; (iii) the ’’factors’’ meaning in PARAFAC is completely different [28]. There is no meaning here except one summand of the decomposition; (iv) The number of summands (’’factors’’) is here determined without any physical interpretation, but only as the number of terms needed to converge and is not predefined. Therefore, the proposed procedure can be interpreted as a particular adaptation of PARAFAC that uses piecewise linear functions to weigh the nodal data points that do not lie on the multi-dimensional grid. For further comparison with the here proposed fitting procedure, the PARAFAC model was fit with an Alternative Least Square (ALS) procedure [25].

Between the most effective multidimensional approximation technique for fitting problems we can mention Sparse Grid [29], where the classification problem is discretized and solved in a sequence of conventional grids with uniform mesh sizes in each coordinate direction. The sparse grid solution is then obtained from the solutions on these different grids by linear combination. For a D-dimensional problem, the sparse grid approach employs only a complexity O(2n ðlog 2n ÞD1 ), where 2n gives the mesh size, then curse of the dimensionality of full grid methods affects sparse grids much less. However, they still depend exponentially on the dimensionality D. This limits the sparse grid method to at most 18 dimensions [30]. Meanwhile, the PGD approximation computational complexity is O(N  d), where N is the number of degrees of freedom on each dimensional axis and D is the dimension. Response surface methodology (RSM) was initially developed to model experimental tests [31] relating the output of a sequence of pre-designed experiments in terms of a set of independent (input) variables. Soon, it was extended to modelling numerical experiments [9–11,3,12]. Although RSM appeared in chemical problems, nowadays it may be found in a vast number of applications, such as, tool-life testing, multi component reactions, food studies or reverse engineering, among a crowd [3,32]. As most metamodeling techniques, the basic approach for RSM involves the definition of (a) an experimental (or virtual) design for generating data; (b) a model to represent the data; and (c) a fitting approach. There are a variety of options for each of these components as described in [6]. RSM usually employs central composite designs, second order approximation functions and least squares regression with their associated pros and cons [33]. For example, polynomials functions drive to high computational cost in manyparameter problems and non-physical oscillations in many occasions [33]. Local approximations (piecewise continuous cubic polynomial splines, B-splines, NURBS, etc.) have been also used in RSM, among a crowd of other possibilities [34–37]. As in any other approximation method, computational cost increases dramatically with the problem dimension as well as the number of samples needed to get a good accuracy. The former problem is solved here by using PGD. Following this method, the response is then approximated by a sum of products of one-dimensional functions, each depending on one independent variable and described, as in the standard Finite Element Method, by 1-D splines with small support. The latter is mitigated by using a good design of experiments, usually abbreviated as DoE. The objective of DoE is selecting the samples (combinations of independent factors) where the response should be evaluated to maintain the accuracy while reducing at a maximum the number of the samples needed and then the cost of constructing the response surface. Several DoEs have been proposed starting from the evident and most expensive (but most accurate) full factorial [38], fractional factorial, central composite [16], D-optimal, Taguchi’s [14], Latin hypercube [39], Audze-Eglais’ approach [40], etc. Finally, getting the control (nodal) values of these splines is here solved by least square minimization of the error between the approximation function and the initially known cloud of samples, although using any other fitting procedure is straightforward. The paper is established as follows. After this Introduction, a brief description of RSPGD for a general multidimensional case is presented in Section 2. Then several numerical examples are introduced in Section 3, discussing the possibilities of the method and its cost and accuracy when compared to other more traditional methods as functional fitting using polynomial approximation functions and the Levenberg–Marquardt fitting procedure. Section 4 shows a more complex and real problem of designing a cranial plate using seven design parameters. We finish the paper with some concluding remarks and additional discussion.

545

F. El Halabi et al. / Comput. Methods Appl. Mech. Engrg. 253 (2013) 543–557

2. The RSPGD approach

E¼

The RSM problem can be resumed as follows: Get the best (in a certain sense) approximation response function f ðx1 ; x2 ; . . . ; xD Þ partially defined by a cloud of values of that function coming from experiments, numerical calculation or any other possible source, and dependent on the values of D quantitative factors or input variables x ¼ fx1 ; x2 ; . . . ; xD g that are considered capable of exact control and that define a hyper surface in the bounded region X  RD (called the experimental region) in which we consider constrained the values of x by practical limitations. Polynomials of first and second order, exponential, sinusoidal, logarithmic, and many other types of functions have been used as approximating functions while their associated constants are determined by means of some kind of regression analysis, being least squares the most used of these strategies. Of course, other fitting procedures are found in the literature such as weighted least squares regression, best linear unbiased predictor used for Kriging, back-propagation mostly implemented in neural networks, and many others [7,8,12,16,5,17]. In particular, least-squares polynomial curve fitting is available in many standard computer programs solving the approximation problem in one and two dimensions. On the contrary, the extension to many dimensions is not frequently found and has not been implemented due to its intrinsic complexities [41]. One reason is the high number of parameters to handle for a high-dimensional fitting that, as commented above, increases exponentially with the number of dimensions. However, and as commented, the interest in solving high-dimensional problems has encouraged the development of model reductions techniques that try to avoid the so-called ‘‘curse of dimensionality’’ that appears in traditional strategies such as multifactorial mesh-based discretization or sparse grid methods [42]. The most important model reduction methods are based in the idea that the response of complex models can be approximated by the one of a surrogate model which represents the projection of the initial model on a low dimensional reduced basis [24]. The way of defining and constructing this basis distinguishes the specific model reduction technique. In this work, a model reduction method based on separation of variables named as Proper Generalized Decomposition (PGD) [19,43,20,21] has been used. PGD have been introduced in different contexts like parameterized PDEs [44–47], multiscale models [48,49], parametric modeling and structural optimization [50], multi-dimensional PDEs [19,43,21,24], etc. This technique constructs the approximation of the solution by means of a sum of products (sometimes called finite sum decomposition) of separated one-dimensional functions depending on each of the problem dimensions or parameters. Each of these functions is determined by an iterative method, with no initial assumption on their form, although most of the times, they are expressed as piecewise linear splines with small support as in standard linear one-dimensional finite elements. Let wðxÞ be a scalar function of x ¼ ðx1 ; x2 ; . . . ; xD Þ 2 X ¼ QD k¼1 ½lk ; Lk . In the PGD approach, this function is approximated as:

FðxÞ 

T D X Y

ai

i¼1

F ki ðxk Þ

k¼1

F ki ðxk Þ ¼

T X

ai F i

ð1Þ

i¼1

with the i  th one-dimensional functions of the k  th variable xk that have to be determined along the process, D the number of independent variables (dimensions), F i the product of D functions F ki and T the number of sums or terms for the approximation. The least square approach is then established as follows: given a set of pairs (xm ; wm ) with m ¼ 1; . . . ; N known values of the function w for N combinations of independent parameters, find the function FðxÞ expressed as in (1) which minimizes E defined as:

N X ½wm  Fðxm Þ2

ð2Þ

m¼1

F ki ðxk Þ ¼

Mk X Nkn ðxk Þzkin

ð3Þ

n¼1

with N kin ðxk Þ the standard linear one-dimensional spline shape function with value 1 at the associated node n and zero at the rest of the Mk nodes of the chosen discretization of the interval ½lk ; Lk ; zkin the nodal value vector of the function F ki at node n. In order to compute the minimum of E the partial derivatives with respect to the parameters of FðxÞ must be zero. The numerical scheme to solve the problem here presented consists in an iterative procedure as follows: let us assume T terms of the finite sum already computed, then, the following three stages proposed in [19,20] are followed for solving the problem: 1. Get the coefficients ai from the linear system of equations derived from the minimization of the functional E with respect to the coefficients ai N X @E @F ¼ 2 ½wm  Fðxm Þ ¼ 0; i ¼ 1; . . . ; T : @ ai @ ai m¼1

ð4Þ

The equivalent linear system may be expressed as:

Ka ¼ f

ð5Þ

For a better understanding of the procedure a system of only two terms in the sum (T ¼ 2) and two dimensions (D ¼ 2) will be considered in the following. In this simple case, we have for the matrix and vectors in (5) the following fully-developed expression:

K¼

 N  X F 1 ðxm ÞF 1 ðxm Þ F 1 ðxm ÞF 2 ðxm Þ m¼1

a¼



F 2 ðxm ÞF 1 ðxm Þ F 2 ðxm ÞF 2 ðxm Þ 



N X a1 wm F 1 ðxm Þ ; f¼ wm F 2 ðxm Þ a2 m¼1

ð6Þ

 ð7Þ

Note that the size of K is T  T, and F i ; I ¼ 1; . . . ; T represent the T Q Q functions F i ðxÞ ¼ Dk¼1 F ki ðxk Þ ¼ Dk¼1 Nk ðxk Þzki , being Nk the onedimensional shape functions for the each dimension xk and zki the vector of nodal values for the summand i and the dimension xk , both with Mk components. 2. A convergence check for the overall solution, that is, if for the already computed values of the basis functions F ki ðxk Þ, i ¼ 1; . . . ; T; k ¼ 1; . . . ; D and coefficients ai ; i ¼ 1; . . . ; T the value of the relative error  is below a certain predefined accuracy limit.

¼

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi E < TOL PN 2 m¼1 ½wm 

ð8Þ

with E given by (2). Then, the solving process finishes and if not, we move to the stage below. In some cases additional stopping criteria are needed, as shown later in the numerical examples section. These problems appears in cases with data containing random perturbations, data input distribution by LHS DoE and oversampled cases, where the residual (8) can not be improved by adding more terms into the approximation. In addition to the relative error stopping criterion (8), a maximum number of terms in the approximation is established. Also, the weighting value (aT Þ for the last computed term is compared to the one corresponding to the first one, and if aT < a1 with  ¼ 103 , a non-important improvement of the approximation generated by adding additional terms might be

546

F. El Halabi et al. / Comput. Methods Appl. Mech. Engrg. 253 (2013) 543–557

expected and the algorithm stops. Another, stopping criterion can be introduced by checking the relative error for each term:

er ¼

~ ðx0 ÞkL2 kuðx0 Þ  u kuðx0 ÞkL2

T D D X Y Y ai F ki ðxk Þ þ Rk ðxk Þ i¼1

k¼1

ð10Þ

k¼1

Here F kTþ1 ðxk Þ has been substituted by Rk ðxk Þ since the normalized basis functions F kTþ1 ðxk Þ will be obtained by normalizing the functions Rk ðxk Þ once the iterative process described below converges. To compute the enrichment functions, a non-linear problem with as many equations as the number of dimensions (k ¼ 1; 2; . . . ; D), obtained by replacing (10) into (2), must be solved. The non-linear solver implemented is an alternating directions (fixed-point) scheme. In previous works, Newton schemes were successfully applied [19]. However the fixed-point alternative based on the resolution of smaller linear systems within an iterative scheme seemed to be faster, simpler and more robust [43]. The non-linear system is then transformed into a set of D independent linear systems that are solved cyclically until convergence. Each basis function Rk ðxk Þ is computed from the linear system established by making null the partial derivative of the error function E with respect to the parameters of such function, assuming the rest of the functions known. To illustrate this technique let us consider again the previous simple case (T ¼ 2 and D ¼ 2), being r and s the nodal value vectors for the first and second dimension functions, respectively. The alternating direction procedure to get the enrichment functions is then expressed as: compute rðx1 Þ with sðx2 Þ known, by means of: " !# N T X X T @E 1T 1T 2T 1 1 2 2 2 1 2 ¼2 wm  ai N ðxm Þzi N ðxm Þzi þ N ðxm ÞrN ðxm Þs @r m¼1 i¼1   T  N1 ðx1m ÞN2 ðx2m Þs ¼ 0 ð11Þ After operating and rearranging the equation above, the following system of equations is obtained:

K1  r ¼ f 1  l1 K1 ¼

N X N1 ðx1m ÞN2T ðx2m ÞsN1T ðx1m ÞN2T ðx2m Þs

ð12Þ ð13Þ

m¼1

f1 ¼

N X wm N1 ðx1m ÞN2T ðx2m Þs

ð14Þ

m¼1

l1 ¼

N X T X

e enr 1 ¼ kriter  riter1 k2 < TOL1

ð9Þ

where uðx0 Þ represent the vector of the response values for param~ ðx0 Þ its respective approximation. This error shows eter set x0 and u to converge to a fix value in problems where the input data values do not match with nodal distribution of the discretization, i.e., random or LHS distributions. 3. A new term T þ 1 is added to the finite sum, so the new basis functions F kTþ1 ðxk Þ, k ¼ 1; . . . ; D have to be obtained. In this ’’enrichment stage’’ [19,21] the response function is then written as:

Fðx1 ; . . . ; xD Þ ¼

the stopping criterion is activated here when the following conditions are both fulfilled:

ai N1 ðx1m ÞN1T ðx1m Þz1i ðsT N2 N2T ðx2m Þz2i Þ

ð15Þ

m¼1 i¼1

Once r is computed, s is obtained from the associated linear system assuming r known. After solving the two linear systems (for the bi-dimensional case), the convergence is checked verifying that both functions reach a fixed point. Considering riter1 ; siter1 and riter ; siter the computed functions r; s at the previous and present iterations respectively,

e enr 2 ¼ ksiter  siterþ1 k2 < TOL2

ð16Þ

At this point, it is important to remark that, even when the non-linear solver converges sufficiently fast, this step consumes the main part of the global computing time. The vectors for starting this stage (enrichment) when seeking the first sum of the approximation are random, subsequent additional starting vectors are defined as the previous converged sum. In previous tests, this procedure provided good convergence rate. This algorithm has several advantages: it is easy to implement, it is robust and simple to extend to high order dimensions. However, the residual decreases almost linearly with the iterations while other methods can provide, at least in principle, superlinear or even quadratic convergence rate. It is also worth mentioning, even if it is not relevant for the purposes of this work, that several types of constraints can be imposed to the loading vectors in a relatively straightforward way. Furthermore, the convergence to a fix point of a basis function within this enrichment stage is plotted in Fig. 1, for different dimensions. The error checked was the associated to (9) where almost linear convergence can be observed. For a better comprehension of the process above, a summary is presented in algorithm 1. Algorithm 1. RSPGD Methodology 1: 2: 3: 4: 5: 6: 7: 8: 9: 10: 11: 12: 13: 14: 15:

Load xm and wm SET i = 1 Number of sums of the approximation While  > TOLls do Initialize rk for k ¼ 2; . . . ; D While e enr k > TOLk do for k = 1, . . ., D do Solve (1) and Compute rk using (9) Compute e enr k using (15) end for end while pffiffiffiffiffiffiffiffiffiffiffiffiffiffi T Normalize rk by rk = rk  rk Compute ai by (5) Compute  by (8) Set i ¼ i þ 1 end while

3. Validation and numerical examples A generic (D-dimensional) algorithm following the procedure described in the previous section was programmed within MATLAB R2010a and a set of numerical examples were solved to ensure the proper working of the method. All the examples were run in a standard CPU with the following specification: processor Intel (R) core I5 of 3.20 GHz, 4.00 GB in RAM and a 64 bits operative system. Different aspects were inspected while solving these examples such as, relative error with different number of approximation terms, time of execution and relation between accuracy and discretization. Furthermore, comparisons were done with the multidimensional response surface Levenberg–Marquardt algorithm with polynomials functions [51]. Also, the three and nine factors problems solution with RSPGD were compared with the CANDECOMP/PARAFAC algorithm [27,25]. The inputs needed for the surface fitting problem solved by RSPGD is the cloud of sample points (values of the response

547

F. El Halabi et al. / Comput. Methods Appl. Mech. Engrg. 253 (2013) 543–557

Table 1 Results of the 2 factors fitting problem for the sixth-order polynomial function data (17).

1

Fix Point Error

10

Discretization n Sample distribution n Number of sample

−2

10

−5

4 Elements n Full factorial n 25

10

3 Dimensions −7

7 Dimensions

10

1

2

4

RSPGD 3th D-PF

9 Dimensions

4 Elements n LHS n 50

−9

10

Method

6

8

10

12

Iteration

RSPGD 3th D-PF

Fig. 1. Convergence curve for the fix point alternative directions algorithm of the first basis function in different dimensional problems.

function for a predefined set of combinations of the input parameters), discretization (number and location of the nodal points along each dimension), experimental domain (discretization intervals for each input parameter) and the maximum error desired for the approximation. As for other fitting techniques, the number of input data and its distribution represents an important issue in the final accuracy and computational cost of the proposed methodology. Even more, the approximation error values will tend asymptotically to a fix value governed by the discretization used, where no improvement will be observed by adding a new term to the approximation. Fully factorial and Latin Hypercube Sample (LHS) [39] designs have been used to establish the input data distributions for the next numerical examples. The proposed methodology allows each parameter to have a different discretization. Due to this, a previous study of the effect of each factor onto the response is recommended, with the purpose of optimizing the discretization for each input parameter. A previous coarse study of the parameter effect on the response is then proposed for high dimensional problems based on the Taguchi DoE. Even more, it is interesting to explore the possibility of building a discretization directly derived from the given points; this could improve CPU time due to the option to incorporate a different discretization size in each dimension. Furthermore, non-uniform discretization can be straightforward introduced, enhancing the approximation where more sample points are found. However, in the next numerical examples we show that the discretization is not critical for CPU time, but it is for the number of sample points needed when designing an experiment in predefined ranges of parameters. CPU time increased linearly with the number of sample points. 3.1. Two factors problems For a simple visualization of the possibilities of the approach, preliminary examples with two and three factors are solved. Different response functions with two dimensions were studied, namely a six degree polynomial and an exponential function. The template functions implemented in the two dimensions fitting problems were:

f ðx; yÞ ¼ ðx6 þ y6 þ x  yÞ=100 f ðx; yÞ ¼ eðx

2 þy2 Þ

þxy

ð17Þ ð18Þ

For the polynomial structure (17) two cases were analyzed with different discretizations, four and six elements, respectively. The intervals of variation were equal for both variables, x; y 2 ½0; 3. A least square multidimensional curve fitting based on the Levenberg– Marquardt algorithm was used for comparative purposes. This curve fitting strategy was implemented for a third order polynomial

6 Elements n Full factorial n 49

RSPGD 3th D-PF

6 Elements n LHS n 70

RSPGD 3th D-PF

Approximation terms

Execution time, s

a

1 2 3 –

0.013 0.041 0.051 0.285

25.40 7.85 0.02 –

1 2 5 –

0.017 0.050 0.144 0.347

20.8 7.5 0.76 –

1 2 3 –

0.017 0.048 0.067 0.283

34.0 10.9 0.03 –

1 2 5 –

0.02 0.08 0.28 0.36

32.31 9.78 0.74 –

Value

approximation function, which from now on will be denoted by 3th D-PF. A full factorial and LHS distribution were used for input data. Table 1 shows the Design of Experiments (DOE) used, number of samples, computer time and a’s values (weighting coefficients) for each RSPGD approximation term. Approximated surfaces for the six element per direction discretization using a full factorial and a 70 points LHS distribution are plotted in Fig. 2(a) and (b), respectively. Two relative errors obtained by the approximation obtained through the RSPGD and the standard regression approximation were computed as in (9), the first of these errors refers to the input data values and the second to response values associated to a number of test parameter sets chosen randomly. These errors will be referred from now on as local and global errors, respectively. Local and global errors against number of approximation terms for the RSPGD procedure for different discretizations and data distribution can be observed in Fig. 2(c) and (d), respectively. A second study was accomplished for the input data function (17) with a random perturbation or ‘‘noise’’ between 0% and 10% of its response value. The sample distribution analyzed for this case was a full factorial design resulting in 49 sample points (6 elements per parameter). Results were compared to the surface obtained by 3th D-PF. Relative errors computed by (9), execution time and a values are shown in Table 2. Comparing this result with the previous case without the random perturbation, it was shown that for the latter case with two terms in the approximation, with a full factorial sample distribution, a local relative error of 104 was achieved, where the a value of the third term is three orders of magnitude less than of the second term. The basis vectors for this two terms are plotted in Fig. 3(a). Instead, the case with a random noise for the same problem configuration needs five terms to reach a relative error of 1004 and an a value three orders of magnitude lower than the one corresponding to the second term (see Table 2). These extra terms represent the fitting of the approximate surface to the random noise introduced. Note that the first two basis vectors have similar behavior to the previous case (without perturbation) as can be observed in Fig. 3(b). However, neither the basis functions values and the corresponding a of these first two terms are identical to the case without perturbation, concluding that the method tries to fit the noise with all its approximation terms as shown in Fig. 4. The third theoretical example with two factors corresponds to the response function (18). Computer time, relative error values

548

F. El Halabi et al. / Comput. Methods Appl. Mech. Engrg. 253 (2013) 543–557

Approximate Surface Input Points Test Points

Approximate Surface Input Points Test Points

15

Response

Response

15 10 5 0

10 5 0

3 2 Y

1

1

0 0

3

3

2

2 Y

X

1

0 0

0.2

0.4 4 Elements (Factorial) 4 Elements (LHS) 6 Elements (Factorial) 6 Elements (LHS) 3th D−PF (LHS 70 points)

X

4 Elements (Factorial) 4 Elements (LHS) 6 Elements (Factorial) 6 Elements (LHS) 10 Elements (Factorial) 3th D−PF (LHS 70 points)

0.05

Relative Error

Relative Error

3

2

(b)

(a)

0.3

1

0.1

0.005

0.05 0

0.001 1

2

3

4

5

Approximation Terms

1

2

3

4

5

6

7

Approximation Terms

(d)

(c)

Fig. 2. Results corresponding to the response function (17). (a) Approximated 3-terms RSPGD response surface with six elements per direction and factorial data distribution. (b) Approximated 5-terms RSPGD response surface with six elements per direction and 70 data points with LHS distribution. (c) Local relative errors curves comparison. (d) Global relative errors curves comparison for 50 test points.

Table 2 Results of the 2 factors fitting problem solved by RSPGD and 3th order polynomial approximation function for the response function (17) with perturbation in the input data. Method

RSPGD

3th D-PF

Execution time, s

Local relative error (approximation terms)

Alpha value

0.021 0.061

0.295 (1) 1:2  102 (2)

32.48 10.56

0.089

2:4  103 (3)

0.43

0.11

1:2  103 (4)

0.07

0.16

3:1  104 (5)

0.04

0.31

5:9  102

–

and relation between the number of samples and computer time were checked, again compared with a 3th D-PF. In this example, two discretization cases with equal number of elements along x and y were evaluated. The corresponding interval was set to x; y 2 ½0; 1 for both parameters. Table 3 summarizes the results obtained with the RSPGD approximation for the different cases and number of sample points. Again, local and global errors were computed against number of approximation terms for different discretization and input data distributions. These errors are summarized in Fig. 5(a) and (b) for local and global errors, respectively. To study the behaviour with different number of data we take fifty, hundred, thousand, and ten-thousand data points. In Fig. 5(c), we show the relation of the global error obtained through a 3-terms RSPGD approximation for 100 random test points and the different number of data point cases. Computer time for computing a 3-terms RSPGD approximation with different number of data points is presented in Fig. 5(d). For the 6th order polynomial function, a 2-terms RSPGD approximation was sufficient to get local residuals in the order of 1009 , for a full factorial DoE of both discretizations cases,

lower than the one of the 3th D-PF, due to the coincidence of the samples points with discretization nodes, a direct consequence of the full factorial DoE. For the cases with a LHS DoE, where sample data did not match discretization nodes, four terms were needed to get similar local errors than the 3th D-PF. However, for LHS distribution the approximation with 5-terms presents lower global errors than the full factorial case, as shown in Fig. 2(d). In the same figure the convergence of the approximation to a fixed global error value can be clearly observed, depending on the discretization. Another important checking parameters are the a’s values that can be taken as an indicator of the contribution of each new term to the total approximation. Observing these a’s values for this case it is immediate to notice its decreasing behavior. With respect to CPU time, for the second discretization and full factorial DoE with 49 sample points to get a 3-terms approximation, 0.067 s were needed, which represents approximately one quarter the cost to get the 3th D-PF approximation. For the 70 LHS DoE case, a 5-terms approximation leads to a time cost of 0.28 s, also lower than for the 3th D-PF (77%). Fig. 2(a) and (b) show the improvement of the approximation, when increasing the number of terms and the discretizations. The second analysis, where a random perturbation was added to the input values, shows lower local errors values than for the 3th D-PF, specially near the samples points, due to the local nature of RSPGD. Similar results for the example with the sample function (18) can be observed in Fig. 5(a) and (b). Computing the approximation for this last case for different number of samples corresponding to a LHS distribution, we can observe how the global errors of the approximation tend to a fix value while increasing the number of data points. Also a linear increment in CPU time can be observed when increasing the number of data points.

549

F. El Halabi et al. / Comput. Methods Appl. Mech. Engrg. 253 (2013) 543–557 0.8

0.8

F1 F2 F3

0.6

G1 G2 G3

0.6

0.4

0.2

0.2

F

G

0.4

0

0

−0.2

−0.2

−0.4

−0.4

−0.6

−0.6 0

1

2

3

0

0.8 0.6

2

3

2

3

(a)

G1 G2 G3 G4 G5

1

1

y

x

F1 F2 F3 F4 F5

1 0.8

0.4

0.4

F

G

0.2 0 0 −0.2 −0.4 −0.4

−0.6 −0.8

−0.8

−1 0

1

2

3

0

y

1

x

(b) Fig. 3. Resulting basis functions for x parameter (F) and y parameter (G) of the two factors fitting problem corresponding to the response function (17). (a) without random perturbation, (b) with a random perturbation between 0% and 10% of the response value.

3.2. Three factors problems As for the two parameters problems, RSPGD is here applied to problems with the response function dependent of three factors. Two different input data functions have been analyzed:

f ðx; y; zÞ ¼ 0:1y  0:3x4 þ z4 ;

x; y; z 2 ½0; 3

f ðx; y; zÞ ¼ sinðxyÞ þ cosð2yÞ þ

10 ; zþxþ1

ð19Þ

h pi x; y 2 0; ; z 2 ½0; 5 4 ð20Þ

The results obtained through these example are also compared to those computed by the CANDECOMP/PARAFAC algorithm. For the first 3-D example (19), all cases studied were compared to a 3th D-PF (30 terms!) and to the results obtained by means of the PARAFAC-ALS algorithm implemented using the optimized MATLAB function parafac als [52]. Note that for the PARAFAC-ALS procedure, the response data needs to be known at every xijk i ¼ 1; . . . ; I; j ¼ 1; . . . ; J; k ¼ 1; . . . ; K, where i, j and k are nodes of the unidimensional discretization of parameter x; y and z, respectively. In other words, in its standard form this procedure is restricted to fully factorial distribution. Full factorial and LHS DoE were used, the corresponding results for different discretizations and with

different numbers of terms are summarized in Table 4. Fig. 6(a) and (b) show the effect of number of approximation terms and discretization onto the accuracy with 4-terms in the RSPGD approximation. The second uses (20) as response function; two discretizations were run and the approximation evolution compared with a 3th D-PF and PARAFAC-ALS. Discretization, DoE, sample points, computational time and a’s values results are shown in Table 5. Local and global errors for (5  2  5) elements for the different approximation methods are plotted in Fig. 7(a) and (b), respectively. Observing the results obtained for the 4th order polynomial function (19), local error curves representing results of a 4terms RSPGD for two different discretizations with a LHS DoE, again show lower values than for the full factorial case. Solving the problem with the PARAFAC-ALS the computed local errors for each approximation term were similar to those obtained through the RSPGD. Global errors (computed for 1000 test points) for a discretization of 10, 5 and 10 elements for parameter x; y and z, respectively, show slightly lower values solving the fitting problem with RSPGD than those obtained by PARAFAC-ALS. To reach global errors obtained by the 3th D-PF, only a 2-term RSPGD approximation with the discretization mentioned before, and 4-term RSPGD for a discretization of 5  2  5 were needed.

550

F. El Halabi et al. / Comput. Methods Appl. Mech. Engrg. 253 (2013) 543–557

Fig. 4. Contribution to the approximation of the first two terms and the rest.

Table 3 Results of the 2 factor problem with a sinusoidal response function defined by (18) for different discretizations obtained by RSPGD and 3th order polynomial approximation function. Discretization n

Method Approximation terms Execution a Value time, s

Sample distribution n Number of sample 4 Elements n Full factorial n 25

1 2 3

0.014 0.04 0.052

3th D-PF –

0.285

1:2  1016 –

1 2 5 3th D-PF –

0.014 0.049 0.140 0.3

4.05 0.81 0.09 –

1 2 3

0.016 0.042 0.061

6.36 1.24

RSPGD

4 Elements n LHS n 50

RSPGD

6 Elements n Full factorial n 49

RSPGD

6 Elements n LHS n 70

3.92 0.89

3th D-PF –

0.3

3:0  1016 –

1 2 5 3th D-PF –

0.022 0.067 0.21 0.32

5.55 1.14 0.11 –

RSPGD

Computational time increment is more remarked in the RSPGD by adding a new term, due to the progressive construction of the basis function. For the PARAFAC-ALS technique the approximation terms must be established ‘a priori’, since all basis functions are computed at once, justifying the low variation in CPU time between a 1-term and 4-terms approximation, as shown in Tables 4 and 5. For large number of samples, lower CPU times were observed when using PARAFAC-ALS procedure. Even more, no significant increment of the PARAFAC-ALS time was registered between  200 and  725 samples. However, for the RSPGD a three times increment in CPU time was observed between these two number of sample points. For the sinusoidal response function (20), local errors of about 103 for a 5-terms approximation function with two different

discretizations were obtained for the RSPGD and the PARAFACALS approximations. Also, the improvement in local approximation when a full factorial DoE is used has to be remarked. For an 8terms RSPGD and PARAFAC-ALS approximations, this error reaches a value of 107 . However, in the case of LHS distribution with approximately the same number of data points (200), this local error converges to  8  103 . Comparing with the approximation with 3th D-PF, local errors of  3  102 , were obtained. Again, lower global errors were obtained with an LHS distribution. For full factorial DoE, global errors by RSPGD and PARAFAC-ALS, were very similar. We also provide the CPU time related to the online computations, i.e., the time required to compute an approximate solution for a new design point. For the 3-D example, it can be seen that generating 1000 approximate solutions using a 5-terms RSPGD or PARAFAC-ALS for the 5  2  5 case discretization, would take around 0.1 s, which is faster than generating the same 1000 points by means of the 3th D-PF approximation, that consumes a total time of around 1.5 s. 3.3. Nine factors problem In the following subsection, RSPGD will be applied to a theoretical moderately high dimensional problem. The function considered was: 2

f ¼ A2 þ EðB

þC 2 Þ

þ DEFG þ K 2 L;

A; B; C; D; E; F; G; K; L 2 ½0; 1 ð21Þ

An important task to complete the setup of the RSPGD problem is the discretization size of each problem parameter. For this example, the Taguchi’s Method was used to achieve such task. Taguchi’s Method for designs of experiments has become increasingly popular for its simplicity and effectiveness. For additional details of the method the reader is referred to [14,13,15,53]. A 9 three-level factor experimental design was implemented, from which the effect of each factor on the main output variable was obtained. Fig. 8 shows the effect of each factor on the response. With this, a guide to

551

F. El Halabi et al. / Comput. Methods Appl. Mech. Engrg. 253 (2013) 543–557 0.2

0.01 0.001

Relative Error

Relative Error

0.25 0.1

4 Elements (Factorial) 4 Elements (LHS) 6 Elements (Factorial) 6 Elements (LHS) 3th D−PF (LHS 70 points)

0 1

2

3

4

5

6

1

2

3

4

5

Approximation Terms

(a)

(b) 10

Computational Time, s

Relative Error

0.005

Approximation Terms

0.01

0.005

0.002 3 Elements per direction 6 Elements per direction 3rd Degree Polinomial Fit

0.001 50

0.05

0.001

7

4 Elements (Factorial) 4 Elements (LHS) 6 Elements (Factorial) 6 Elements (LHS) 10 Elements (Factorial) 3th D−PF (LHS 70 points)

100

1000

10000

6

7

1

3 Elements per direction 6 Elements per direction

10

10

10

0

−1

−2

50

100

1000

Nº. Sample Points

Nº. Sample Points

(c)

(d)

10000

Fig. 5. Comparison between RSPGD and a 3rd order polynomial fitting corresponding to the response function (18): (a) local relative errors curves, (b) global relative errors curves, (c) global relative errors curve vs. number of samples with an LHS distribution, (d) CPU time vs. number of samples with an LHS distribution.

Table 4 Results of the 3 factor problem with a polynomial response function defined by (19) for different discretizations obtained by RSPGD, PARAFAC-ALS and the 3th D-PF. Discretization n Sample distribution n Number of sample

Method

RSPGD 525 n Full factorial n 196

PARAFAC-ALS

3th D-PF 525 n LHS n 200

RSPGD

3th D-PF RSPGD 10  5  10 n Full factorial n 726

PARAFAC-ALS

3th D-PF

a

Approximation terms

Execution time, s

1 2 3 4 1 2 3 4 –

0.06 0.157 0.220 0.330 0.17 0.18 0.18 0.20 0.441

443.8 103.7 1.35 0.09

1 2 3 4 –

0.09 0.25 0.43 0.59 0.45

430.9 100.8 10.9 3.8 –

1 2 3 4 1 2 3 4 –

0.22 0.51 0.78 1.07 0.19 0.21 0.21 0.23 0.42

637.1 148.0 2.1 0.15

Value

DoE was tested for both discretization cases, resulting in 7776 and 26244 points for discretizations 1 and 2, respectively. Then, an LHS DoE was created and analyzed with 3000 and 5000 points for discretizations 1 and 2, respectively. Results obtained were compared to the PARAFAC-ALS algorithms for the full factorial case. A quadratic polynomial regression fitting for comparison was also built. The polynomial fitting expression writes,

~¼ u

9 9 9 X 9 X X X C i a2i K i ai þ Pij ai aj þ G; i¼1

–

–

choose the discretization size for each parameter is obtained; for example, more elements must be used for factor A than for factor B, due to the stronger non-linear behaviour shown in Fig. 8a. The interval for each problem parameter x was set to xi 2 ½0; 1 with i ¼ 1; . . . ; 9. Following the parameter effect curves, two discretizations were analyzed, as detailed in Table 6. A full factorial

i¼1

ð22Þ

i¼1 j¼iþ1

where 55 coefficients are determined by the Levenberg–Marquardt algorithm. As for the previous examples, local and global errors were computed, being the latter obtained from 5000 test points. In terms of local errors for both discretization cases with full factorial DoE, lower values and faster convergence rate were registered using the PARAFAC-ALS procedure as shown in Fig. 9(a). RSPGD shows relatively slow local convergence, i.e., a 40-terms approximation was needed to get similar local error than the 15-terms PARAFAC-ALS approximation. This considerable difference was not observed when global errors were computed. Fig. 9(b) shows the convergence to a global error value associated to a specific discretization. Here again a faster convergence to the fix global error is observed using the PARAFAC-ALS technique, observing that for the second discretization with a full factorial distribution a 5terms PARAFAC-ALS approximation leads to a relative error of 2:9  102 , while an 11-terms RSPGD approximation was needed to reach such error value. The 9-D quadratic polynomial fitting registered relatively high global errors ( 0:3) in comparison with the other fitting procedures. Running the same discretization cases with an LHS DoE, a similar rate of convergence than the full factorial cases was observed, as mentioned previously. As for the 3-D examples, again global

552

F. El Halabi et al. / Comput. Methods Appl. Mech. Engrg. 253 (2013) 543–557

Relative Error

10 10 10

10

RSPGD−3x1x3(Factorial) RSPGD−3x1x3(LHS) RSPGD−5x2x5(Factorial) RSPGD−5x2x5(LHS) RSPGD−10x5x10(Factorial) 3th D−PF (LHS 200 points) PARAFAC−10x5x10(Factorial)

1

0

−1

−2

Relative Error

10

−3

RSPGD−3x1x3(Factorial) RSPGD−3x1x3(LHS) RSPGD−5x2x5(Factorial) RSPGD−5x2x5(LHS) 3th D−PF (LHS 200 points) PARAFAC−5x2x5(Factorial)

−6

1

2

0.1

0.01 0.005 3

4

1

2

3

Approximation Terms

Approximation Terms

(a)

(b)

4

Fig. 6. Comparison between RSPGD, PARAFAC-ALS and a 3rd order polynomial fitting corresponding to the response function (19): (a) local relative errors curves, (b) global relative errors curves.

Table 5 Results of the 3 factor problem with a response function defined by (20) for different discretizations obtained by RSPGD, PARAFAC-ALS and the 3th D-PF. Discretization n Sample distribution n Number of sample

Method

RSPGD 636 n Full factorial n 196

PARAFAC-ALS

3th D-PF 636 n LHS n 200

RSPGD 3th D-PF

RSPGD 848 n Full factorial n 405

PARAFAC-ALS

3th D-PF 848 n LHS n 400

RSPGD 3th D-PF

Approximation terms

Execution time, s

a

1 2 3 4 5 1 2 3 4 5 –

0.05 0.15 0.25 0.35 0.47 0.15 0.19 0.23 0.24 0.25 0.44

63.6 4.18 1.72 0.59 0.15

1 2 3 5 –

0.09 0.22 0.43 0.75 0.45

62.5 4.09 1.87 0.92 –

1 2 3 4 5 1 2 3 4 5 –

0.09 0.27 0.42 0.63 0.79 0.17 0.19 0.23 0.19 0.22 0.44

89.3 5.51 2.26 0.81 0.02

1 2 3 5 –

0.09 0.27 0.39 0.72 0.45

62.5 4.09 1.87 0.92 –

Value

–

–

errors achieved with an LHS distribution are slightly lower than for full factorial cases, even using lower numbers of samples, as shown in Fig. 9(c). With respect to the CPU time, a notable difference can be observed in Fig. 9(d) between the RSPGD and the PARAFAC-ALS, where for the first discretization with full factorial DoE (7776) points, 15-terms are computed in around 4 min and with 3000 points with an LHS distribution the same RSPGD approximation consume around 1.7 min getting better global error values. The CPU time to generate 1000 new design points with a 15terms RSPGD approximation for the 9-D problem was around 0.7 s, that represents an attractive feature of the proposed method.

4. Application example In this section, one of the many possible applications for this proposed surface fitting procedure is presented. Basically, we propose to get the specific response of a complex finite element model for any combination of its design factors, by means of RSPGD. A cranial implant manufactured by Rapid Prototyping (RP) has been chosen, Fig. 10. At the moment RP technology is available for a number of biocompatible thermoplastics and metal alloys suitable for prosthesis fabrication. Among the different RP technologies, Selective Laser Sintering (SLS) is the most versatile in terms of the variety of materials that can be used for fabricating the prototypes. First, a CAD representation of the desired part to be fabricated (i.e., the STL model of the prototype) is transformed into thin layers. Following, a high power laser selectively fuses small particles of material on a fixtureless table to produce the desired product layer by layer by scanning each of the virtual layers from the transformed CAD model on the surface of a powder bed. After a given layer is completed, the fixtureless base is move down one layer thickness and a new layer of material is applied on top. The process is repeated until the component is completed [54,55]. In most cases this manufacturing technique induces a transversely isotropic final product material, where significant differences in mechanical properties appear between the growth direction and its perpendicular plane. Anisotropic mechanical properties of the final product can be controlled somehow by manufacturing parameters such as: camera temperature, velocity of building directions, angle of the fixtureless table, etc. Clinically, two basic variables are the main ones to be checked in a cranial implant computational model: maximum displacements and stresses against an applied load, ensuring brain tissue and implant integrity. The finite element model of the implant used in this work consists of an implant shell model with a discretization of 1875 quadrilateral elements. Contact conditions are introduced to simulate the interaction with cranium. A seven factor response surface will be generated implementing the RSPGD approximation. The parameters chosen are: two coordinates, Lx and Ly , representing the applied load location on the implant, (see Fig. 10), two angles that define the applied load direction in 3D, axyz and ayz , two elastic moduli, Eyz and Ex , for plane yz and x directions, respectively and the angle of the fixtureless table (hm ) was considered as another parameter design. For this example, the discretization size of each parameter is chosen according to the effect of each variable on the response, which was obtained through a previous Taguchi’s approach. A 7 three-level factor experimental design was implemented, from which the effect of each factor on the main output variable was obtained.

553

F. El Halabi et al. / Comput. Methods Appl. Mech. Engrg. 253 (2013) 543–557 −1

0.1

Relative Error

Relative Error

10

−2

10

RSPGD (Factorial) −3

PARAFAC (Factorial)

10

3th D−PF (Factorial)

RSPGD (Factorial) PARAFAC (Factorial) 3th D−PF (Factorial) RSPGD (LHS 200 points) 3th D−PF (LHS 200 points)

0.05

0.02

RSPGD (LHS 200 points) 3th D−PF (LHS 200 points) −4

10

0.01 1

2

3

4

5

1

2

3

Approximation Terms

Approximation Terms

(a)

(b)

4

5

Fig. 7. Data function (20) for a 5-terms RSPGD and PARAFAC-ALS and a 3rd order polynomial fitting with a discretization of 5, 5 and 2 elements in directions x, y and z, respectively: (a) local relative errors curves, (b) global relative errors curves.

Fig. 8. Effect of each factor from a 9 three-level factor Taguchi experiment design on the response following function (21).

Table 6 Elements in each parameter discretization of the 9-D example. Discretization case

A

B

C

D

E

F

G

K

L

1 2

3 3

2 2

1 2

2 2

2 2

1 2

2 2

2 2

1 2

Fig. 11 shows the effect of each factor on the maximum displacement. According to these curves, the response has an important

correlation with Lx , Ly and (Ex ). Also note that for Eyz the variation in the response is almost linear. A more exhaustive Taguchi experimental design can be made, i.e., 7 four-level factor, to have a better understanding of the behavior of the response. According to this analysis, the parameters intervals and discretization were chosen as shown in Table 7. With the parameters and their intervals defined as in Table 7, the discretization for each parameter must be established. A full factorial design was implemented to establish the sample points and its distribution. 3888 sample resulted for the 7

F. El Halabi et al. / Comput. Methods Appl. Mech. Engrg. 253 (2013) 543–557

Relative Error

10

10

10

RSPGD Discretization 1 PARAFAC Discretization 1 RSPGD Discretization 2 PARAFAC Discretization 2 2th D−PF Discretization 2

−1

−2

−3

1

3

5

7

9

11

13

15

0.3

0.1

0.05

0.02

1

3

5

7

9

11

Approximation Terms

Approximation Terms

(a)

(b)

0.3

13

15

400 200

RSPGD Local Discr. 1 RSPGD Local Discr. 2 RSPGD Global Discr. 1 RSPGD Global Discr. 2 2th D−PF Global Discr. 2

0.1

100

CPU Time, s

0.2

Relative Error

RSPGD Discretization 1 PARAFAC Discretization 1 RSPGD Discretization 2 PARAFAC Discretization 2 2th D−PF Discretization 2

0.2

Relative Error

554

0.05

RSPGD Discr. 1 (Fac) RSPGD Discr. 2(Fac) PARAFAC−ALS Discr. 1(Fac) PARAFAC−ALS Discr. 2(Fac) RSPGD Discr. 1 (LHS) RSPGD Discr. 2(LHS)

10

1 0.02

1

3

5

7

9

11

13

15

1

2

4

6

8

Approximation Terms

Approximation Terms

(c)

(d)

10

12

Fig. 9. Comparison between RSPGD and PARAFAC-ALS corresponding to the 9-d response function (21). (a) Local errors curves for different discretizations with a full factorial sample distribution. (b) Global errors curves for different discretizations with a full factorial sample distribution. Local and global errors curves for different discretizations with an LHS sample distribution. (d) CPU time for different discretization and samples distribution.

Fig. 10. Cranial implant geometry.

dimension mesh according to the discretization shown in Table 7. As mentioned before two response variables were analyzed: the maximum displacement and the maximum principal stress. The output of the RSPGD algorithm for this problem are the vectors that represent the basis functions of the approximation and its respective a values, with which it is possible to evaluate the response variable at any point of the 7-dimensional domain. Note that the RSPGD algorithm shown in Algorithm 1 will be inside another loop over the nodes of the finite element mesh of the cranial implant model. A script within MATLAB R2010a that compute the response approximation (on-line) for each node corresponding to a specific combination of the parameters was programmed. The input for this script are the desired factor combinations and the basis functions of the approximation terms with their respective a values obtained by the RSPGD procedure (off-line). Our cranial implant finite element model contains 1958 nodes and it takes less than 2 s to compute the response approximation

(on-line) for all nodes and write them in a visualization file for a new design set of parameters. Stopping criteria used to run the RSPGD approximation for all model nodes were: (a) Maximum number of terms for the RSPGD (15 for this example); (b) the last a value computed is less than 1% of the a corresponding to the first term of basis function; (c) local relative error defined by (9) reaches a tolerance of 104 . Each FEM simulation of the cranial implant model takes approximately 10 min. The nodes response was finally approximated by different numbers of terms. For example, CPU time for a node (central node 263) with 15-terms RSPGD is summarized in Table 8. Note that each RSPGD analysis is well suited for parallel computing since nodal computations are completely independent of each other. Fig. 12 shows a distribution of the displacement in the x direction computed by the finite element method and the RSM by the proposed technique (RSPGD) for the following factors combination: Lx ¼ 22:5 mm, Ly ¼ 27:0 mm, axyz ¼ 0:7854 rad, ayz ¼ 3:1416 rad, Eyz ¼ 2500:0 MPa, Ex ¼ 1000:0 MPa and hm ¼ 0:5236 rad. Errors between 20 samples points used as sample data and its approximation by RSPGD were computed. The mean of the error and the standard deviations registered were 1:17%  0:92 and 1:9%  1:1 for the displacement and maximum stress response, respectively. Also the response of 10 random factor combinations by FEM simulation and by RSM were computed and compared. The mean error and standard deviation resulted 1:59%  1:02 and 2:79%  1:25 for the displacement and stress response, respectively. An observation on the evolution of the approximation by RSPGD through the number of terms for a central node (263) was made. a’s values effectively decrease with increasing the number of terms as shown in Table 8.

555

F. El Halabi et al. / Comput. Methods Appl. Mech. Engrg. 253 (2013) 543–557

Fig. 11. Effect of each factor from a 7 three-level factor Taguchi experiment design on the maximum displacement response.

Table 7 Intervals and discretization of each parameter for the cranial implant problem. Parameter

Initial value

Final value

Number of elements

Lx (mm) Ly (mm) axyz (rad) ayz (rad) Eyz (MPa) Ex (MPa) hm (rad)

0 0

45 55 p=2 2p 4000 2000 p=3

3 3 2 2 2 3 2

p=4 0 1000 1000 0

Table 8 Results of the RSPGD approximation of the displacement in the x-direction for the central node 263. Method RSPGD Fitting

1-term 7-terms 10-terms 12-terms 15-terms

Local error

Execution time, s

a Value

0.31 0.10 0.052 0.023 0.011

0.72 30.5 54.8 90.2 143.3

0.1190 0.0357 0.0263 0.0122 0.0098

556

F. El Halabi et al. / Comput. Methods Appl. Mech. Engrg. 253 (2013) 543–557

Fig. 12. Direction x displacements distribution for factors combination: Lx ¼ 22:5 mm, Ly ¼ 27:0 mm, axyz ¼ 0:7854 rad, ayz ¼ 3:1416 rad , Eyz ¼ 2500:0 MPa, Ex ¼ 1000:0 MPa and hm ¼ 0:5236 rad: (a) FEM simulation, (b) RSPGD.

5. Discussion and conclusions

References

A methodology has been presented that integrates the proper generalized decomposition (PGD) method with the response surface methodology (RMS), here referred to as RSPGD. It was shown through some analytical examples with two and three dimensions that the response surface approximation obtained by RSPGD offers a clear advantage in cases where the sample points correspond to complex functions, i.e., sinusoidal and 6th or higher-order polynomials. Reductions in computer time and residual magnitude are clearly observed when comparing with the traditional fitting procedure based in polynomial approximation functions solved by the Levenberg–Marquardt algorithm. Unlike traditional fitting procedures where prescribed approximation function types are used, in RSPGD no previous constraint in the approximation is established which means higher flexibility to fit functions with high variations. However, due to its local character, discretization has to be made in such a manner that sufficient information is provided in each element to ensure convergency. Full factorial and LHS DoE have shown to be suitable for this fitting procedure ensuring good global and local accuracy. This technique can be also used in problems with data containing random errors and oversampling. Even more, linear increase in computational cost was observed with the number of data points, although obtaining such sample points could be expensive. Further optimizations can be considered for CPU time reductions, such as, advanced programming and implementation in a more optimal programming language (FORTRAN, C, C++, etc.). The principal advantage of this method appears in high dimensional problems (>5); in this work, an application with seven design factors was analyzed. Once the response surface was built (on-line) for each node, which takes an average of 90 s per node, any new response evaluation for any specific combination of factors is obtained in about 1.8 s. These results can be improved by several ways, one of them can be by adding more enrichments basis functions to the approximation. However, we have to keep in mind the compromise between computer time and the accuracy of the approximation. Another way is decreasing the size of the discretization for each factor. A study to optimize these factor’s discretizations is proposed, where each factor effect over the response is obtained through a previous Taguchi’s design of experiments. Then, the factors with higher effect may be discretized finer.

[1] M.B. Wilk, O. Kempthorne, Some aspects of the analysis of factorialexperiments in a completely randomized design, Annals of Mathematical Statistics 27 (4) (1956) 950–985. [2] G.E.P. Box, K.B. Wilson, On the experimental attainment of optimum conditions, Journal of the Royal Statistical Society Series B – Statistical Methodology 13 (1) (1951) 1–45. [3] W.J. Hill, W.G. Hunter, A review of response surface methodology – a literature survey, Technometrics 8 (4) (1966) 571–590. [4] S.N. Patnaik, J.D. Guptill, D.A. Hopkins, Subproblem optimization with regression and neural network approximators, Computer Methods in Applied Mechanics and Engineering 194 (30–33) (2005) 3359–3373. [5] P. Koutsourelakis, Design of complex systems in the presence of large uncertainties: A statistical approach, Computer Methods in Applied Mechanics and Engineering 197 (49–50) (2008) 4092–4103. [6] R. Jin, W. Chen, T.W. Simpson, Comparative studies of metamodeling techniques under multiple modeling criteria, Structural and Multidisciplinary Optimization 23 (2000) 1–13. [7] M. Nakashima, Variable coefficient a-stable explicit Runge-Kutta, methods (1995). [8] S. Sakata, F. Ashida, M. Zako, On applying kriging-based approximate optimization to inaccurate data, Computer Methods in Applied Mechanics and Engineering 196 (13–16) (2007) 2055–2069. [9] O.L. Davies, The Design and Analysis of Industrial Experiments, Hafner, 1954, p. 495 (Chapter 11). [10] D.R. Read, The design of chemical experiments, Biometrics 10 (1) (1954) 1–15. [11] H. Grohskopf, Statistics in the chemical process industries – present and future, Industrial and Engineering Chemistry 52 (6) (1960) 497–499. [12] H. Arellano-garcia, J. Schoneberger, S. Korkel, Optimal design of experiments in the chemical engineering, Chemie Ingenieur Technik 79 (10) (2007) 1625– 1638. [13] V.N. Nair, B. Abraham, J. Mackay, A. Box, R.N. Kacker, T.J. Lorenzen, J.M. Lucas, R.H. Myers, G.G. Vining, J.A. Nelder, M.S. Phadke, J. Sacks, W.J. Welch, A.C. Shoemaker, K.L. Tsui, S. Taguchi, C.F.J. Wu, Taguchi parameter design – a panel discussion, Technometrics 34 (2) (1992) 127–161. [14] K.L. Tsui, An overview of taguchi method and newly developed statisticalmethods for robust design, Iie Transactions 24 (5) (1992) 44–57. [15] K.N. Otto, E.K. Antonsson, Extensions to the taguchi method of product design, Journal of Mechanical Design 115 (1) (1993) 5–13. [16] H.J. Park, S.H. Park, Extension of central composite design for second-order response surface model building, Communications in Statistics – Theory and Methods 39 (7) (2010) 1202–1211. [17] P. Langley, H.A. Simon, Applications of machine learning and rule induction, Communications of the ACM 38 (11) (1995) 55–64. [18] P. Ladevèze, L. Chamoin, On the verification of model reduction methods based on the proper generalized decomposition, Computer Methods in Applied Mechanics and Engineering 200 (23–24) (2011) 2032–2047. [19] A. Ammar, B. Mokdad, F. Chinesta, R. Keunings, A new family of solvers for some, classes of multidimensional partial differential equations encountered in kinetic theory modeling of complex fluids, Journal of Non-Newtonian Fluid Mechanics 139 (3) (2006) 153–176. [20] F. Chinesta, A. Ammar, E. Cueto, Recent advances and new challenges in the use of the proper generalized decomposition for solving multidimensional models, Archives of Computational Methods in Engineering 17 (4) (2010) 327– 350. [21] D. Gonzalez, A. Ammar, F. Chinesta, E. Cueto, Recent advances on the use of separated representations, International Journal for Numerical Methods in Engineering 81 (5) (2010) 637–659. [22] C. Audouze, F. De Vuyst, P.B. Nair, Reduced-order modeling of parameterized pdes using time-space-parameter principal component analysis, International Journal for Numerical Methods in Engineering 80 (8) (2009) 1025–1057.

Acknowledgements This work has been partially financed by Ascamm Technology Center in Barcelona. The authors also thank the Instituto de Salud Carlos III (ISCIII) through the Ciber-BBN initiative.

F. El Halabi et al. / Comput. Methods Appl. Mech. Engrg. 253 (2013) 543–557 [23] Y. Liang, H. Lee, S. Lim, W. Lin, K. Lee, C. Wu, Proper orthogonal decomposition and its applications part i: Theory, Journal of Sound and Vibration 252 (3) (2002) 527–544. [24] A. Nouy, A priori model reduction through proper generalized decomposition for solving time-dependent partial differential equations, Computer Methods in Applied Mechanics and Engineering 199 (23–24) (2010) 1603–1626. [25] G. Tomasi, R. Bro, A comparison of algorithms for fitting the parafac model, Computational Statistics & Data Analysis 50 (7) (2006) 1700–1734. [26] N.D. Sidiropoulos, G.B. Giannakis, R. Bro, Blind PARAFAC receivers for DSCDMA systems, IEEE Transactions on Signal Processing 48 (3) (2000) 810–823. [27] N. Faber, R. Bro, P.K. Hopke, Recent developments in candecomp/parafac algorithms: a critical review, Chemometrics and Intelligent Laboratory Systems 65 (1) (2003) 119–137. [28] R. Harshman, Foundations of the parafac procedure: Models and conditions for an explanatory multi-modal factor analysis, UCLA Working Papers in Phonetics 16 (1) (1970) 1–84. [29] J. Garcke, M. Hegland, Fitting multidimensional data using gradient penalties and the sparse grid combination technique (Apr. 2009). [30] H.-J. Bungartz, M. Griebel, Sparse grids, Acta Numerica 13 (2004) 147–269. [31] G.E.P. Box, N.R. Draper, Empirical model-building and response surface, John Wiley & Sons, Inc., New York, NY, USA, 1986. [32] U. Alibrandi, N. Impollonia, G. Ricciardi, Probabilistic eigenvalue buckling analysis solved through the ratio of polynomial response surface, Computer Methods in Applied Mechanics and Engineering 199 (9–12) (2010) 450–464. [33] H. Akima, A new method of interpolation and smooth curve fitting based on local procedures, Journal of the ACM 17 (4) (1970) 589–602. [34] C.D. Boor, Bicubic spline interpolation, Journal of Mathematics and Physics 41 (3) (1962) 212–&. [35] C.A. Hall, Natural cubic and bicubic spline interpolation, SIAM Journal on Numerical Analysis 10 (6) (1973) 1055–1060. [36] K.L. Rasmussen, P.V. Sharma, Bicubic spline interpolation – quantitative test of accuracy and efficiency, Geophysical Prospecting 27 (2) (1979) 394–408. [37] H. Inoue, A least-squares smooth fitting for irregularly spaced data – finiteelement approach using the cubic b-spline basis, Geophysics 51 (11) (1986) 2051–2066. [38] D.J. Street, L. Burgess, Factorial Designs, John Wiley & Sons, Inc., 2007. [39] R.J. McKay, Adequacy of variable subsets in multivariate regression, Technometrics 21 (4) (1979) 475–479. [40] S.J. Bates, J. Sienz, D.S. Langley, Formulation of the audze-eglais uniform latin hypercube design of experiments, Advances in Engineering Software 34 (8) (2003) 493–506. [41] F.H. Lesh, Multi-dimensional least-squares polynomial curve fitting, Communications of the ACM 2 (9) (1959) 29–30.

557

[42] B. Khoromskij, O(dlog n)-quantics approximation of n-d tensors in highdimensional numerical modeling, Constructive Approximation 34 (2) (2011) 257–280. [43] A. Ammar, B. Mokdad, F. Chinesta, R. Keunings, A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modelling of complex fluids: Part ii: Transient simulation using space-time separated representations, Journal of Non-Newtonian Fluid Mechanics 144 (2–3) (2007) 98–121. [44] A. Nouy, A. Clement, F. Schoefs, N. Moes, An extended stochastic finite element method for solving stochastic partial differential equations on random domains, Computer Methods in Applied Mechanics and Engineering 197 (51–52) (2008) 4663–4682. [45] A. Nouy, Recent developments in spectral stochastic methods for the numerical solution of stochastic partial differential equations, Archives of Computational Methods in Engineering 16 (3) (2009) 251–285. [46] A. Nouy, O.P. Le Maitre, Generalized spectral decomposition for stochastic nonlinear problems, Journal of Computational Physics 228 (1) (2009) 202–235. [47] A. Doostan, G. Iaccarino, A least-squares approximation of partial differential equations with high-dimensional random inputs, Journal of Computational Physics 228 (12) (2009) 4332–4345. [48] F. Chinesta, A. Ammar, E. Cueto, Proper generalized decomposition of multiscale models, International Journal for Numerical Methods in Engineering 83 (8–9) (2010) 1114–1132. [49] D. Neron, P. Ladeveze, Proper generalized decomposition for multiscale and multiphysics problems, Archives of Computational Methods in Engineering 17 (4) (2010) 351–372. [50] A. Leygue, E. Verron, A first step towards the use of proper general decomposition method for structural optimization, Archives of Computational Methods in Engineering 17 (4) (2010) 465–472. [51] M. Lampton, Damping-undamping strategies for the levenberg-marquardt nonlinear least-squares method, Computers in Physics 11 (1) (1997) 110–115. [52] E. Acar, D.M. Dunlavy, T.G. Kolda, Cpopt: Optimization for fitting candecomp/ parafac models (extended abstact), in: CASTA 2008: Workshop on Computational Algebraic Statistics, Theories and Applications, 2008. [53] S. Kumar, P. Kumar, H.S. Shan, Optimization of tensile properties of evaporative pattern casting process through taguchi’s method, Journal of Materials Processing Technology 204 (1–3) (2008) 59–69. [54] C. Chua, K. Leong, C. Lim, Rapid prototyping: Principles and applications, in: World Scientific, 2003. [55] R. Bibb, D. Eggbeer, P. Evans, A. Bocca, A. Sugar, Rapid manufacture of customfitting surgical guides, Rapid Prototyping Journal 15 (5) (2009) 346–354.