Implementation and Performance Evaluation of a Bivariate Cut-HDMR Metamodel for Semiconductor Packaging Design Problems with a Large Number of Input Variables

A metamodeling technique based on Bivariate Cut High Dimensional Model Representation (Bivariate Cut HDMR) is implemented for a semiconductor packaging design problem with 10 design variables. Bivariate Cut-HDMR constructs a metamodel by considering only up to second-order interactions. The implementation uses three uniformly distributed sample points (s = 3) with quadratic spline interpolation to construct the component functions of Bivariate Cut-HDMR, which can be used to make a direct comparison with a metamodel based on Central Composite Design (CCD). The performance of Bivariate Cut-HDMR is evaluated by two well-known error metrics: R-squared and Relative Average Absolute Error (RAAE). The results are compared with the performance of CCD. Bivariate Cut HDMR does not compromise the accuracy compared to CCD, although the former uses only one-fifth of sample points (201 sample points) required by the latter (1045 sample points). The sampling schemes and the predictions of cut-planes and boundary-planes are discussed to explain possible reasons for the outstanding performance of Bivariate Cut HDMR.


Introduction
Numerous metamodeling techniques (also known as response surface methods, surrogate models, or reduced-order models) have been developed and implemented for engineering design optimization [1]. Metamodeling includes two parts: generation of discrete sample points and connection of the discrete sample points. Each metamodeling technique possesses its own characteristics that can be suited for certain applications.
For a typical engineering system, a metamodel considering up to second-order interactions is often sufficient to describe system responses [2,3]. For example, a metamodeling technique called central composite design (CCD) has been implemented widely in the field of semiconductor packaging design community, which uses quadratic polynomial functions for fitting sample points [4][5][6]. It was implemented for commercial software such as optiSLang [7], Design-Expert [8], etc. The CCD metamodeling technique requires P number of sample points to produce the metamodel for N number of input variables, defined as [3]: As the number of input variables increases, the computational cost may become prohibitively high due to an extremely large number of sample points required (this is the well-known "curse of dimensionality"). This situation will be exacerbated when the modeling requires a computationally expensive analysis such as time-dependent and nonlinear analysis that is routinely encountered in complex semiconductor packaging architectures [9][10][11].
In order to build accurate and efficient metamodels for high dimensional input-output systems, numerous advanced metamodeling techniques such as the high-dimensional model representation (HDMR) technique [12][13][14], reduced dimensional polynomial chaos expansion [15], and active and rank-adaptive tensor regression [16] have been developed to enhance the efficiency of metamodeling in various engineering fields.
Among these techniques, one family of HDMR, called Cut-HDMR, possesses two unique practical features: (1) it involves function evaluations only at sample points, and, more importantly, (2) it determines the number of sample points from a pre-defined function of the number of input variables regardless of the nature of engineering applications, i.e., selection of sample points is simple and straightforward [14,17,18]. Based on these features, numerous metamodeling techniques based on Cut-HDMR have been developed such as RBF-HDMR [19], Adaptive MLS-HDMR [20], and Kriging-HDMR [21].
HDMR decomposes a multivariate function into multiple lower-order component functions, based on the hierarchical structure of interaction effects of the input variables. The high performance of some of the metamodeling techniques based on Cut-HDMR considering up to second-order component functions (this will be referred to as Bivariate Cut-HDMR) has been confirmed for nonlinear numerical test functions [19,20] and statistical analysis of multiconductor transmission line networks [22].
The objectives of this paper are (1) to introduce the cut-HDMR to the semiconductor packaging design community and to help implement the Bivariate Cut-HDMR for those who are not familiar with the HDMR, and (2) to investigate the performance of Bivariate Cut-HDMR for a complex semiconductor packaging problem (10 design input valuables). The result is compared with the performance of CCD, which has been utilized widely in the semiconductor packaging industry.

Background: Bivariate Cut-HDMR
The fundamentals of HDMR are described first. A specific HDMR that uses the Dirac measure located at a cut center, called Cut-HDMR, is presented together with its approximated version, Bivariate Cut-HDMR, which considers up to second-order component functions.

High-Dimensional Model Representation (HDMR)
The concept of high-dimensional expansion was implemented originally to estimate the sensitivity of a function with respect to arbitrary groups of variables [23]. Later, the term, HDMR, was first introduced by Rabitz and Alis [12]. They detailed and completed the general foundations of HDMR.
The HDMR expansion is performed based on the interaction effects of input variables. The term "interaction" employed here means that more than one variable act together to affect the performance function. This is distinctly different from the term "correlation" employed in statistics, which represents whether and how strongly a pair of random variables are related.
A general form of HDMR is defined as [12]: where y(x) and the bold letter, x, represent the performance function and the vector of input variables, (x 1 , x 2 , . . . , x N ), respectively; y 0 is a constant representing the mean of the performance function, which is called zeroth-order effect or mean effect; y i (x i ) represents the effect when the variable x i acts independently on y(x), which is called first-order effect or main effect; y ij x i , x j is the effect on y(x) when the variables x i and x j act together, which is called second-order effect or bivariate interaction effect. It should be noted that y ij x i , x j excludes the main effects of x i and x j as well as the mean effect. The subsequent terms indicate the higher order interaction effects of more variables acting together on y(x). The last term y 12...N (x 1 , x 2 , . . . , x N ) represents the residual influence. Each effect in the general form of HDMR is called a component function. The general form of component functions can be expressed as [12]: The functions in the above equation are defined as: . . , n} is an n-dimensional unit cube and γ is a measure [24]. A measure is a function that quantifies the size of sets. A measure assigns a non-negative real number or +∞ to subsets of a certain set. Each distinct measure embodies a different way to assess how big a set is.
There is no unique decomposition of the model output y(x 1 , x 2 , . . . , x N ); all HDMR expansions follow the general form in Equation (2). The choice of a particular HDMR expansion depends on the application and the nature of any constraints in sampling input variables. For example, for the uncertainty analysis of a model output (e.g., an analysis of the variance of an output), the component functions in the HDMR should be chosen to represent the independent contributions of input variables to the overall uncertainty of the output. It is known as ANOVA-HDMR [12,25].
The ANOVA-HDMR is typically carried out by multi-dimensional Monte Carlo integration due to its complexity. The Monte Carlo integration needs a large number of sample points to attain good accuracy. It is impractical for the advanced semiconductor packaging applications where computational cost of each sample point is high. However, another approach of HDMR, Cut-HDMR, can tackle the challenge.

Cut-HDMR and Bivariate Cut-HDMR
Cut-HDMR uses the Dirac measure [26] located at a point m = (m 1 , m 2 , . . . , m n ) (also known as cut center): By combining it with Equations (3) and (4), the component functions of Cut-HDMR can be expressed as: where y(m 1 , . . . , m i−1 , x i , m i+1 , . . . , m n ) is a 1D performance function along the x i direction that passes through m; y m 1 , . . . , x i , . . . , x j , . . . , m n is a 2D performance function of the x i , x j plane that passes through m, and so on. Equation (6) shows that Cut-HDMR is an expression as a superposition of its values on lines, planes, and hyperplanes of higher orders passing through the cut center, m. The expansions of Cut-HDMR do not contain any integral. Cut-HDMR uses only arithmetic computation to determine the component functions, and thus it requires the least amount of computational cost compared to other HDMRs [12,27].
For most well-defined physical systems, the high-order interactions are negligible [14,28], and thus the multivariate performance function of such a physical system can be approximated well by the sum of low-order component functions. Experience shows that an HDMR expansion up to the second-order often provides a satisfactory description of the function for many high-dimensional systems when the input variables are properly chosen [14].
It has been proven that the mean values of input variables, µ, are the optimal cut center m when only the terms up to the second-order are considered [25]. Accordingly, the metamodel based on Bivariate Cut-HDMR can be obtained by substituting Equation (6) into Equation (2) with the cut center being the mean values of input variables. This Bivariate Cut-HDMR metamodel is written as [29]: where y 0 = µ = [µ 1 , µ 2 , . . . , µ N ] T is the vector of the mean values of N input variables (cut center); µ ∼i is µ without the element µ i ; µ ∼ij is µ without the elements µ i and µ j ; y x i , µ ∼i is the 1D performance function along the x i direction that passes through µ (cut-line); and y x i , x j , µ ∼ij is a 2D performance function on the (x i , x j ) plane that passes through µ (cut-plane). Figure 1 illustrates the concept of Bivariate Cut-HDMR using an arbitrary 2D function, which is decomposed into four component functions. Figure 1a shows the 2D function, x 2 + y 2 + xy − 14x − 16y + 122 = 0, as the black meshed surface, and the dot in the figure represents the zeroth-order effect (i.e., a constant). In Figure 1b, the blue curve is the 1D performance function along the x 1 direction, in which x 2 is kept as µ 2 . The green line is the zeroth-order effect along the x 1 direction. The main effect of x 1 is the red curve, obtained by subtracting the green line from the blue curve.
The same procedure can be applied to obtain the main effect of x 2 , as shown in Figure 1c. In Figure 1d, the blue surface is obtained by the superposition of the red curves in Figure 1b,c, which represents the performance function without any interaction effects. The green plane is the zeroth-order effect. By subtracting the blue surface and the green plane from the black surface, the interaction effect of the (x 1 , x 2 ) pair is obtained, which is shown as the red surface. ,  Figure 1d, the blue surface is obtained by the superposition of the red curves in Figure 1b,c, which represents the performance function without any interaction effects. The green plane is the zeroth-order effect. By subtracting the blue surface and the green plane from the black surface, the interaction effect of the 12 ( , ) xx pair is obtained, which is shown as the red surface.

Implementation for Semiconductor Packaging Application
The Bivariate Cut-HDMR technique is implemented to construct a metamodel for a semiconductor packaging application. The application involves warpage prediction of a thin flat ball grid array (TFBGA) package with 10 design input variables. Figure 2 shows the schematic diagram of a TFBGA package. The first chip is attached to a substrate by the first die attach film (DAF). The second chip is attached to the first chip by the second DAF. Then, they were encapsulated by epoxy molding compound (EMC). A stacked die TFBGA package is often used as the top package of a Package-on-Package (PoP). Warpage at solder pad areas is one of the most critical factors to high PoP stacking yield [30].

Implementation for Semiconductor Packaging Application
The Bivariate Cut-HDMR technique is implemented to construct a metamodel for a semiconductor packaging application. The application involves warpage prediction of a thin flat ball grid array (TFBGA) package with 10 design input variables. Figure 2 shows the schematic diagram of a TFBGA package. The first chip is attached to a substrate by the first die attach film (DAF). The second chip is attached to the first chip by the second DAF. Then, they were encapsulated by epoxy molding compound (EMC). A stacked die TFBGA package is often used as the top package of a Package-on-Package (PoP). Warpage at solder pad areas is one of the most critical factors to high PoP stacking yield [30].   A finite element (FE) model was constructed for warpage prediction. Figure 3 shows details of the FE model built by a commercial FE analysis package (ANSYS ® ). The quarter symmetry model of boundary conditions and the die stack configuration are shown in (a) and (b); and the enlarged view of cross-section is shown in (c). The material properties and the nominal dimensions used in the model are summarized in Tables 1 and 2. The nominal dimensions of the TFBGA package are adopted from the design in Refs. [10,31]. The TFBGA package was subjected to the EMC molding process at 175 °C, which was used as a stress-free temperature. The conventional lead-free solder reflow profile with the peak temperature of 260 °C was considered [32].

Description of TFBGA Package
In this implementation, 10 design variables were considered for the warpage prediction of solder pad areas. The details of design variables are summarized in Table 3. The design spaces of the package dimensions and the material properties were defined by the values found in the literature: package dimensions in [31,[33][34][35][36][37][38][39] and material properties in [40][41][42][43][44][45].  The TFBGA package was subjected to the EMC molding process at 175 • C, which was used as a stress-free temperature. The conventional lead-free solder reflow profile with the peak temperature of 260 • C was considered [32]. In this implementation, 10 design variables were considered for the warpage prediction of solder pad areas. The details of design variables are summarized in Table 3. The design spaces of the package dimensions and the material properties were defined by the values found in the literature: package dimensions in [31,[33][34][35][36][37][38][39] and material properties in [40][41][42][43][44][45].

Sample Points
The number of sample points to construct a Bivariate Cut-HDMR metamodel can be generally expressed as [12]: where N is the number of input variables and s is the number of sample points taken along the direction of each input variable. N(s − 1) points are used to construct 1D performance functions, and N(N − 1)(s − 1) 2 /2 points are used to construct the 2D performance functions. For the univariate terms (i.e., s number of sample points distributed along each input variable), the center becomes the reference point, and the remaining (s − 1) sample points are evenly distributed on two sides with respect to the reference point. For the bivariate terms, the sample points form a uniform gird on a plane with the center as a reference point.
Cut-HDMR, in its original form [14], states that a set of sample points can be selected to calculate the values of corresponding component functions and to form a look-up table that can be used to interpolate component functions at an arbitrary point in the design domain. There has been no universally accepted sampling strategy and interpolation algorithm. The implementation of this study uses three uniformly distributed sample points (s = 3) with quadratic spline interpolation to construct the component functions of Bivariate Cut-HDMR. In this way, the Bivariate Cut-HDMR metamodel can be compared directly with the CCD metamodel. Figure 4 and Table 4 show the number of sample points required for the CCD and Bivariate Cut-HDMR metamodels. After N = 7, the number of sample points for CCD becomes more than double the number of sample points for Bivariate Cut-HDMR. Considering only the number of sample points, Bivariate Cut-HDMR has a significant advantage over CCD when a metamodeling problem has a large number of input variables.
riate Cut-HDMR. In this way, the Bivariate Cut-HDMR metamodel can be compared directly with the CCD metamodel. Figure 4 and Table 4 show the number of sample points required for the CCD and Bivariate Cut-HDMR metamodels. After 7 N  , the number of sample points for CCD becomes more than double the number of sample points for Bivariate Cut-HDMR. Considering only the number of sample points, Bivariate Cut-HDMR has a significant advantage over CCD when a metamodeling problem has a large number of input variables.

Construct Performance Functions
After the warpage values at the 201 sample points are obtained, the metamodel can be constructed by applying Equation (7). The quadratic spline interpolation scheme was adopted with all sample points to form the 1D performance functions (cut-lines), y x i , µ ∼i , and the 2D performance functions (cut-planes), y x i , x j , µ ∼ij by following the procedures below: • 1D performance functions: 1. Select a design variable.

2.
Find the three sample points along the design variable that was obtained earlier, i.e., high, mid, and low values of the design variable and other design variables keep the mean values.

3.
Construct the 1D function of the design variable with the three sample points using quadratic spline interpolation. This can be done by using the built-in function that is available in commercial software (e.g., MATLAB).

4.
Select another design variable and repeat steps 2-3 until all 1D performance functions along each design variable are built.
• 2D performance functions: 1. Select a pair of design variables.

2.
Find the nine sample points along two design variables that were obtained earlier (other design variables keep the mean values) as shown in the figure.

3.
Construct the 2D function of the design variable with the nine sample points using quadratic spline interpolation. This can be done by using the built-in function that is available in commercial software (e.g., MATLAB).

4.
Select another pair of design variables and repeat steps 2-3 until all 2D performance functions of each pair of design variables are built.
The performance functions are illustrated in Figures 5 and 6. Figure 5 shows the 1D performance functions of EMC thickness and substrate thickness, and Figure 6 shows the 2D performance functions of two pairs of design variables. The pair of substrate thickness and EMC CTE has the strongest second-order interaction effect among other pairs. In contrast, the pair of package width and length and 1st chip thickness has the weakest second-order interaction effect. Red dots indicate the sample points used to construct the cut-lines and the cut-planes. The performance functions are illustrated in Figures 5 and 6. Figure 5 shows the 1D performance functions of EMC thickness and substrate thickness, and Figure 6 shows the 2D performance functions of two pairs of design variables. The pair of substrate thickness and EMC CTE has the strongest second-order interaction effect among other pairs. In contrast, the pair of package width and length and 1st chip thickness has the weakest secondorder interaction effect. Red dots indicate the sample points used to construct the cut-lines and the cut-planes.   Following is the example of determining the response of a random input by using the constructed metamodel. Assuming that a random input x is (0.81, 0. 25 where the warpage at the cut center, y 0 = µ = [µ 1 , µ 2 , . . . , µ 11 ] T , is 0.9 µm; y x i , µ ∼i and y x i , x j , µ ∼ij are the values of x on all known 1D performance functions and 2D performance functions that were constructed earlier. Thus, the warpage value at the random input x can be calculated; it was −49.2 µm. The above Bivariate Cut-HDMR procedure was integrated in MATLAB (R2020b) codes, and they are available at https://www.mathworks.com/matlabcentral/fileexchange/9289 0-bivariate-cut-hdmr (accessed on 25 May 2021). Those who are interested in implementing Bivariate Cut-HDMR metamodeling can run the script readily by following the instructions.

Performance Evaluation
The performance of Bivariate Cut-HDMR is evaluated using two well-known error metrics. The performance of CCD is also evaluated for comparison.

Error Metrics
Two error metrics employed to evaluate the performance are: [46] • Metric 1: R-squared where m is the number of total test sample points; y(x i ) is a performance function at the ith new sample point used for validity check;ŷ(x i ) is an approximated performance function at the ith new sample point; and y(x i ) is the mean of all y(x i ). R-squared indicates the overall accuracy of a metamodel, and its maximum value is 1.
• Metric 2: Relative average absolute error (RAAE) where STD is the standard deviation of all y(x i ). Similar to R-squared, RAAE quantifies the overall accuracy of a metamodel. The closer a value of RAAE is to zero, the more accurate a metamodel is.
Monte Carlo simulation (MCS) was performed to produce 1000 additional sample points. They were used to evaluate the performance of Bivariate Cut-HDMR using the above metrics. The results of the performance metrics are summarized in Table 5. The values of R-squared and RAAE are 0.9855 and 0.0880, respectively. A metamodel based on CCD was also constructed for comparison. A total of 1045 sample points were required for the CCD metamodel, which built a 10D quadratic function to define the warpage behavior. The additional sample points obtained from MCS were utilized again to evaluate the performance of CCD metamodel. The results are also shown in Table 5. The values of R-squared and RAAE are 0.9662 and 0.1472, respectively.
More direct and quantitative comparisons are shown in Figure 7, where the absolute errors of the MCS sample points are compared. The absolute errors of half the MCS sample points of Bivariate Cut-HDMR are less than 5 µm. The outcome is remarkable. Bivariate Cut HDMR used only one-fifth of sample points (201 sample points) required by CCD (1045 sample points). However, Bivariate Cut-HDMR does not compromise the accuracy when compared to CCD. The following section is intended to provide some insight into this performance of Bivariate Cut HDMR. A metamodel based on CCD was also constructed for comparison. A total of 1045 sample points were required for the CCD metamodel, which built a 10D quadratic function to define the warpage behavior. The additional sample points obtained from MCS were utilized again to evaluate the performance of CCD metamodel. The results are also shown in Table 5. The values of R-squared and RAAE are 0.9662 and 0.1472, respectively.
More direct and quantitative comparisons are shown in Figure 7, where the absolute errors of the MCS sample points are compared. The absolute errors of half the MCS sample points of Bivariate Cut-HDMR are less than 5 µ m. The outcome is remarkable. Bivariate Cut HDMR used only one-fifth of sample points (201 sample points) required by CCD (1045 sample points). However, Bivariate Cut-HDMR does not compromise the accuracy when compared to CCD. The following section is intended to provide some insight into this performance of Bivariate Cut HDMR.  Figure 8 shows the sampling schemes of Bivariate Cut-HDMR and CCD for a threevariable ( 3 N  , 3 s  ) example. In the figure, the red point is the mean point for both Bivariate Cut-HDMR and CCD; the blue points are used to construct the functions of three lines in the X-, Y-and Z-directions for Bivariate Cut-HDMR and the axial points for CCD; and the yellow points together with the blue points are used to construct the functions of three planes (X-Z plane, Y-Z plane, and X-Y plane) for Bivariate Cut-HDMR and the factorial points for CCD. It also illustrates one of the cut-planes (green planes) and one of the boundary-planes (magenta planes) of both metamodels.  Figure 8 shows the sampling schemes of Bivariate Cut-HDMR and CCD for a threevariable (N = 3, s = 3) example. In the figure, the red point is the mean point for both Bivariate Cut-HDMR and CCD; the blue points are used to construct the functions of three lines in the X-, Y-and Z-directions for Bivariate Cut-HDMR and the axial points for CCD; and the yellow points together with the blue points are used to construct the functions of three planes (X-Z plane, Y-Z plane, and X-Y plane) for Bivariate Cut-HDMR and the factorial points for CCD. It also illustrates one of the cut-planes (green planes) and one of the boundary-planes (magenta planes) of both metamodels.

variable (
3 N  , 3 s  ) example. In the figure, the red point is the mean point for both Bivariate Cut-HDMR and CCD; the blue points are used to construct the functions of three lines in the X-, Y-and Z-directions for Bivariate Cut-HDMR and the axial points for CCD; and the yellow points together with the blue points are used to construct the functions of three planes (X-Z plane, Y-Z plane, and X-Y plane) for Bivariate Cut-HDMR and the factorial points for CCD. It also illustrates one of the cut-planes (green planes) and one of the boundary-planes (magenta planes) of both metamodels. The sampling points of Bivariate Cut-HDMR are utilized to construct the first-order and second-order component functions, i.e., every sample point is used to construct the 1D and 2D performance functions (as illustrated in Figures 5 and 6). On the other hand, the sample points of CCD are aimed to cover the boundaries of a design domain.

Prediction of Cut-Planes
As mentioned earlier, the sampling scheme of Bivariate Cut-HDMR is designed to construct the cut-lines and cut-planes. The prediction on the cut-planes performed by both metamodels are compared. As shown in Figure 8, Bivariate Cut-HDMR has more sample points (9) than CCD (5) on the cut-planes (green planes). Figures 9 and 10 show the two predicted surfaces (cut-planes), which were studied in the TFBGA application. Each figure has the identical nine dots (warpage values obtained from the FE model) in (a) and (b). Red dots are the sample points used to construct for each metamodel. Blank dots are the sample points that were used to construct the Bivariate Cut-HDMR metamodel but not used to construct the CCD metamodel. The sampling points of Bivariate Cut-HDMR are utilized to construct the first-order and second-order component functions, i.e., every sample point is used to construct the 1D and 2D performance functions (as illustrated in Figures 5 and 6). On the other hand, the sample points of CCD are aimed to cover the boundaries of a design domain.

Prediction of Cut-Planes
As mentioned earlier, the sampling scheme of Bivariate Cut-HDMR is designed to construct the cut-lines and cut-planes. The prediction on the cut-planes performed by both metamodels are compared. As shown in Figure 8, Bivariate Cut-HDMR has more sample points (9) than CCD (5) on the cut-planes (green planes). Figures 9 and 10 show the two predicted surfaces (cut-planes), which were studied in the TFBGA application. Each figure has the identical nine dots (warpage values obtained from the FE model) in (a) and (b). Red dots are the sample points used to construct for each metamodel. Blank dots are the sample points that were used to construct the Bivariate Cut-HDMR metamodel but not used to construct the CCD metamodel.  The surfaces (cut-planes) in Figures 9a and 10a were constructed by Bivariate Cut-HDMR (quadratic spline interpolation) with two sets of nine sample points shown in the The surfaces (cut-planes) in Figures 9a and 10a were constructed by Bivariate Cut-HDMR (quadratic spline interpolation) with two sets of nine sample points shown in the figures. There is no error between warpage values obtained from FE (dots) modeling and the predicted surfaces.
The surfaces in Figures 9b and 10b were plotted by the CCD metamodel obtained from 1045 sample points. The five sample points shown in Figures 9b and 10b were just a small portion of the total 1045 sample points used to construct the CCD metamodel (a second-order polynomial function). This attempt for CCD to fit all 1045 sample points inevitably produces the discrepancy between true warpage values (dots) and predicted surfaces in the entire design domain, especially in the corners, as shown in Figures 9b and 10b.

Prediction of Boundary-Planes
The example in Figure 8 (N = 3) shows five sample points on the boundary-planes of both metamodels. It is important, however, to note that there are lesser or no sample points on the boundary-planes of both metamodels when the number of input variables (N) increases. On the boundary-planes of the TFBGA application (N = 10), there were no sample point for Bivariate Cut-HDMR and only four sample points for CCD. Figure 11 shows two predicted surfaces (boundary-planes) of the TFBGA application. Variables other than the two variables shown in the plots were kept at their maximum values, i.e., it represents one of the boundary-planes in the design domain. Red dots in (b) are the sample points used to construct the CCD metamodel. They also appear in (a), although they are not used for Bivariate Cut-HDMR.  Figures 9b and 10b were just a small portion of the total 1045 sample points used to construct the CCD metamodel (a second-order polynomial function). This attempt for CCD to fit all 1045 sample points inevitably produces the discrepancy between true warpage values (dots) and predicted surfaces in the entire design domain, especially in the corners, as shown in Figures 9b and 10b.

Prediction of Boundary-Planes
The example in Figure 8 ( 3 N  ) shows five sample points on the boundary-planes of both metamodels. It is important, however, to note that there are lesser or no sample points on the boundary-planes of both metamodels when the number of input variables (N) increases. On the boundary-planes of the TFBGA application ( 10  N ), there were no sample point for Bivariate Cut-HDMR and only four sample points for CCD. Figure 11 shows two predicted surfaces (boundary-planes) of the TFBGA application. Variables other than the two variables shown in the plots were kept at their maximum values, i.e., it represents one of the boundary-planes in the design domain. Red dots in (b) are the sample points used to construct the CCD metamodel. They also appear in (a), although they are not used for Bivariate Cut-HDMR. The Bivariate Cut-HDMR surface (boundary-planes) of Figure 11a are plotted by 201 sample points, while the CCD surface of Figure 11b are plotted by 1045 sample points including the four sample points on the boundary-plane. It is noteworthy that the predicted 2D performance function of Bivariate Cut-HDMR is similar to the CCD surface, The Bivariate Cut-HDMR surface (boundary-planes) of Figure 11a are plotted by 201 sample points, while the CCD surface of Figure 11b are plotted by 1045 sample points including the four sample points on the boundary-plane. It is noteworthy that the predicted 2D performance function of Bivariate Cut-HDMR is similar to the CCD surface, despite the fact that CCD utilizes four sample points on the boundary-plane, but Bivariate Cut-HDMR does not.

Conclusions
Bivariate Cut-High Dimensional Model Representation (Bivariate Cut-HDMR) was implemented successfully for the warpage problem of a thin flat ball grid array package with 10 design variables. The implementation with three uniformly distributed sample points (s = 3) in conjunction with quadratic spline interpolation allowed for comparing its performance with a metamodel based on Central Composite Design (CCD).
The performance of both metamodels were evaluated by two well-known error metrics: R-squared and Relative Average Absolute Error (RAAE). The results were compared with the performance of CCD: the R-squared values of CCD and Cut-HDMR were 0.9662 and 0.9855, respectively; the RAAE values of CCD and Cut-HDMR were 0.1472 and 0.0880, respectively.
The outcome was remarkable. Bivariate Cut HDMR used only one-fifth of sample points (201 sample points) required by CCD (1045 sample points); however, Bivariate Cut-HDMR did not compromise the accuracy in both error metrics compared to CCD, which was confirmed by more direct and quantitative comparisons using the absolute errors of the Monte Carlo simulation (MCS) sample points.
Two technical reasons for the outstanding performance of Bivariate Cut-HDMR were discussed: (1) Sampling scheme: the sample points of Bivariate Cut-HDMR were utilized to construct the first-order and second-order component functions, while the sample points of CCD were aimed to cover the boundaries of a design domain. (2) Predictions of cut-planes and boundary-planes: Bivariate Cut-HDMR predicted cutplanes more accurately despite the smaller number of sample points, while both techniques produced similar accuracy for boundary-plane predictions.