A multivariate hybrid approach applied to AISI 52100 hardened steel turning optimization

A multivariate hybrid approach applied to AISI 52100 hardened steel turning optimization

Journal of Materials Processing Technology 189 (2007) 26–35 A multivariate hybrid approach applied to AISI 52100 hardened steel turning optimization ...

210KB Sizes 0 Downloads 35 Views

Journal of Materials Processing Technology 189 (2007) 26–35

A multivariate hybrid approach applied to AISI 52100 hardened steel turning optimization Anderson P. Paiva, Jo˜ao Roberto Ferreira ∗ , Pedro P. Balestrassi Federal University of Itajuba, Industrial Engineering, Av. BPS, 1303 CP 50 Pinheirinho, 37500-903 Itajuba, Minas Gerais, Brazil Received 13 October 2006; received in revised form 5 December 2006; accepted 27 December 2006

Abstract This paper presents an alternative hybrid approach, combining response surface methodology (RSM) and principal component analysis (PCA) to optimize multiple correlated responses in a turning process. Since a great number of manufacturing processes present sets of correlated responses, this approach could be extended to many applications. As a case study, the turning process of the AISI 52100 hardened steel is examined considering three input factors: cutting speed (Vc ), feed rate (f) and depth of cut (d). The outputs considered were: the mixed ceramic tool life (T), processing cost per piece (Kp ), cutting time (Ct ), the total turning cycle time (Tt ), surface roughness (Ra ) and the material removing rate (MRR). The aggregation of these targets into a single objective function is conducted using the score of the first principal component (PC1 ) of the responses’ correlation matrix and the experimental region (Ω) is used as the main constraint of the problem. Considering that the first principal component cannot be enough to represent the original data set, a complementary constraint defined in terms of the second principal component score (PC2 ) is added. The original responses have the same weights and the multivariate optimization lead to the maximization of MRR while minimize the other outputs. The kind of optimization assumed by the multivariate objective function can be established examining the eigenvectors of the correlation matrix formed with the original outputs. The results indicate that the multiresponse optimization is achieved at a cutting speed of 238 m/min, with a feed rate of 0.08 mm/rev and at a depth of cut of 0.32 mm. © 2007 Elsevier B.V. All rights reserved. Keywords: Response surface methodology (RSM); Principal component analysis (PCA); Generalized reduced gradient (GRG); Hard turning; Mixed ceramic

1. Introduction A great number of manufacturing processes present sets of correlated responses. It is mainly the continuous processes which are naturally multivariate, and for a variety of reasons there is a lack of efficiency of their representative models. This paper considers the hard turning process, but its findings can be applied to many other processes. As with most parts of machining processes, hard turning is specific, which means it requires a unique model for every application, raw material, and cutting condition. In looking for the best process operating condition, one must evaluate the behavior of some desired features that are, at first, considered significant as a function of the factor’s increment. This is the traditional experimental strategy. In many situations the multiple quality characteristics measured are highly correlated. There is plenty of literature that

Corresponding author. Tel.: +55 35 36291150; fax: +55 35 36291148. E-mail address: [email protected] (J.R. Ferreira).

0924-0136/$ – see front matter © 2007 Elsevier B.V. All rights reserved. doi:10.1016/j.jmatprotec.2006.12.047

evaluates machining processes with multiple responses [1–9]. Even though the responses present a highly correlated structure, there are no major concerns by the authors regarding the influence of such structures on the models coefficient estimation or to the residual independence, as referenced by Box et al. [10]. The best-known optimization method for multiple responses is the so-called Harrington’s desirability function [11]. According to Khuri and Conlon [12], however, this method is unable to incorporate the correlation structure that exists in the original data set. This means that the interrelationship between the many responses can lead the experimenter to inconclusive results when the analysis is univariate. This kind of solution can differ widely from the simultaneous optimum solution. Many authors try to surmount this problem by using principal component analysis (PCA). Bratchell [13], however, warns of two important drawbacks in doing so: (a) the conflict that exists between maximum and minimum values in a group of variables that needs to be simultaneously optimized and (b) the lower significance of principal components for data with lower intercorrelation. In both cases PCA is not adequate for optimizing

A.P. Paiva et al. / Journal of Materials Processing Technology 189 (2007) 26–35

the responses. This work proposes an alternative procedure of optimization that is able to diminish the drawbacks of the multivariate approach without many constraints for the responses and presents a solution for the case where just one principal component is insufficient to represent the major part of the variation structure. To establish a machining process of high quality means finding out the important parameters and levels that are related to certain quality characteristics. These characteristics are generally nonlinear and sometimes expressed by stochastic models. In the turning process of AISI 52100 steel, the quality characteristics used in this work are the mixed ceramic tool life (T), processing cost per piece (Kp ), cutting time (Ct ), the total turning cycle time (Tt ), surface roughness (Ra ) and material removal rate (MRR). The experimental region (Ω) is considered to be the main constraint. Since these responses exhibit a high degree of dependency, a multivariate strategy can be employed to overcome the influence of these structures on the machining characteristics model building. According to Box et al. [10], ignoring the dependence among the errors, linear dependencies among the expected values of the responses or linear dependencies in the data, can lead to difficulties and wrong conclusions. In the present study, the multivariate response surface approach is adopted to optimize the output variables affecting the machining quality characteristics. First, the experimentation is conducted in the region of interest according to a Box–Wilson central composite design (CCD) when the six responses are recorded and calculated [14]. Second, the regression analysis is conducted using the first principal component score of the original output variables. In spite of the adequate representation proportioned by the latent variables, an alternative procedure is proposed employing a multivariate constraint formed by the second largest principal component score (PC2 ) and an objective function written in terms of the first principal component (PC1 ). Finally, the generalized reduced gradient (GRG) optimization algorithm is used to determine the optimal parameters of the turning process. The large number of experiments necessary to establish an adequate functional relationship between the observed responses and the cutting parameters (cutting speed, feed rate and depth of cut), usually make the experimentation costprohibitive. According to Choudhury and El-Baradie [15], many researchers have been investigating the effects of the cutting parameters – by varying a single parameter per experiment – on the machining outputs, such as tool life and surface roughness. The current study, on the other hand, takes into consideration the simultaneous variation of factors to build forecasting models for all relevant outputs. This approach is known as design of experiment methodology (DOE) and consists of planning experiments capable of generating appropriate data for efficient statistical analysis, which in turn produces valid and objective conclusions [14]. From the many experimental strategies available, this study will make use of the response surface methodology (RSM). The RSM arrangements, such as central composite design (CCD), for instance, use a combination of a factorial arrangement (complete or fractional), average points of the factor levels


(center points) and the axial points (extremes) to adjust, when it is necessary, a second-order polynomial model. Many researchers have employed these methodologies to study a material machinability. Beauchamp et al. [16] employed a complete factorial of six factors to investigate the influence of the cutting parameters on the surface roughness in the turning process. Noordin et al. [17] applied the RSM to describe the performance of carbide tools in the AISI 1045 steel turning process. Choudhury and El-Baradie [15] also used the RSM to model tool life in the turning of high resistance steel. Alauddin et al. [18] developed a similar work. Choudhury and Bartarya [19] employed factorials at three levels to study the influence of temperature on tool wear. Other researchers [20,21] have implemented the Taguchi methodology, along with such other techniques as dual optimization, artificial neural networks (ANN), Monte Carlo simulation and fuzzy logic. The objective of all these works is to optimise – through a small, but efficient number of experiments – the machining process. 2. Response surface methodology overview According to Montgomery [14], RSM is a collection of mathematical and statistical tools used to model and analyze problems whose desired responses are influenced by many variables. In general, the relationship between dependent and independent variables is unknown. One must thus find a reasonable approximation for the real relationship between the response (y) and the set of independent variables (x). Usually a lower-order polynomial in some region of interest is employed. However, if there is curvature in the system, then the approximating function must be a polynomial of higher order, like the second-order model seen in Eq. (1). y = β0 +


βi xi +


k  i=1

βii xi2 +

βij xi xj + ε


where β is the polynomial coefficient, k = 2 and ε is the error term. Montgomery [14] considers it unlikely that a specific polynomial model approximates a real model for the whole experimental space covered for the independent variables. For a specific region, however, the approximation is usually efficient. The ordinary least squares (OLS) method is used to estimate the parameters (β) that in matrix form can be written as: −1 βˆ = (XT X) XT Y


where X is the matrix of factors level and Y is the response. The evaluation of the presence of curvature in the model is based on the analysis of center points for the factors level. Derringer and Suich [11], dealing with multiple response problems, improved the algorithm of desirability function. In this method, the statistical model is first obtained using OLS. Then using a set of transformations based on the limits imposed on the responses, a conversion is conducted for each one of the responses resulting in an individual desirability function di , with 0 ≤ di ≤ 1. These individual values are then combined using a


A.P. Paiva et al. / Journal of Materials Processing Technology 189 (2007) 26–35

geometrical average, such as: D = (d1 (Y1 )d2 (Y2 )· · ·dk (Yk ))1/k


This value of D gives a solution of commitment and is restricted to the interval [0,1]. D is close to 1 when the responses are close to its specification. The type of transformation depends on the desired optimization direction. The desirability function approach to a problem of optimization is simple, easy to apply, and allows the user to judge the importance of each response. However, this methodology does not take into account the variance and correlations of the responses [12]. Ignoring these correlations, one can modify the structure of the overall desirability function, which in turn may jeopardize the establishing of an optimum operating condition. 3. Multivariate response surface approach The principal component analysis is one of the most widely applied tools used to summarize common patterns of variation among variables. Moreover, PCA is able to retain meaningful information in the early axes while summarizing in later axes variations associated with experimental error, measurement inaccuracy, and rounding. According to Johnson and Wichern [22], the PCA method is algebraically a linear combination of p random variables X1 , X2 , . . ., Xp . Geometrically, these combinations represent a selection of a new system of coordinates obtained from an original system rotation. The coordinative axes now have the variables X1 , X2 , . . ., Xp . The new axes represent the direction of maxima. The principal components are uncorrelated and depend only on the matrix of covariance  (or on the matrix of correlation ρ) of the variables X1 , X2 , . . ., Xp and its development does not require the multivariate normality assumption. Assuming that  is the covariance matrix associated to the random vector XT = [X1 , X2 , . . ., Xp ] and that this matrix has pairs of eigenvalues–eigenvectors (λ1 , e1 ), (λ2 , e2 ), . . ., ≥(λp , ep ), where λ1 ≥ λ2 ≥ · · · ≥ λp ≥ 0, then the ith principal component is given by: Yi =

eTi X

= e1i X1 + e2i X2 + · · · + epi Xp ,

i = 1, 2, . . . , p (4)

If the eigenvectors are perpendicular, the ith component will be the result found in Eqs. (5) and (6): Maximize :

Var(Ti X)


Subject to : Ti i = 1 Cov(Ti X, Tk X) = 0,


Sometimes it is useful to write the linear combination in a form of principal component score. In this way, xpn being a ran√ dom observation, x¯ p the pth response average, spp the response standard deviation, p the response and [E] being the eigenvectors matrix of the multivariate set, results in: PCscore = [Z][E]


There are a variety of stopping rules to estimate the adequate number of nontrivial axes or, in other words, the number of significant principal components. The most popular methods are those based on the Kaiser’s criteria [22]. According to this rule, only the principal components whose eigenvalues are greater than 1 should be kept to represent the original set. Moreover, the accumulated variance explained should be greater than 80%. This criterion is adequate when it is used with the correlation matrix. Otherwise, the covariance matrix should only be used for a set of original responses written in the same scale. There are several difficulties when dealing with multivariate responses. Independently modeling each response variable does not take into account the relationships or correlations among the variables. Hence, special care must be taken to avoid misleading interpretations when analyzing multiresponse data. The basic problem relates to the fitting of multiresponse models while ignoring the three kinds of dependencies that can occur: (i) dependence among the errors, (ii) linear dependencies among the expected value of the responses and (iii) linear dependencies in the original data. To overcome these difficulties, a hybrid strategy based on multivariate statistics for summarizing and reducing the dimensionality of the data can be employed [10]. Once the PCA factorizes the multivariate data into a number of independent factors, a natural formulation for the multiresponse problem is to change the original response variables by a principal component score equation modeled through OLS algorithm. This factorization take into account the variances and correlations among the original variables. To force the solution to fall inside the experimental region, a constrained nonlinear programming problem written in terms of principal components could be expressed as shown in Eqs. (8) and (9). Minimize :

k k    βij xi xj PC1 =β0 + βi xi + βii xi2 + i=1

Subject to :

xT x ≤ ρ 2


(8) (9)

Optimum values can be obtained by locating the stationary point of the multivariate fitted surface. The objective is to find the settings of x’s that can optimize the multivariate objective function (PC1 ) subject only to the constraint that defines the region of interest Ω. In other words, the appropriate optimum value of the fitted multivariate model is located by using some optimization algorithm through a constrained procedure to force the optimum to lie within the experimental region. There are two different experimental regions of interest in optimization: spherical and cuboidal. For cuboidal designs, the constraint is written as −1 ≤ xi ≤ 1, i = 1, 2, . . ., k (k is the number of control variables), and for spherical designs the constraint is defined by (xT x) ≤ ρ2 , where ρ is the design radius. The value of ρ should be chosen in order to avoid solutions that are too far outside the experimental region that led to response surfaces established in Eq. (1). For a central composite design, a logical choice is ρ = α, where α is the axial distance. In the case of cuboidal designs (such as Box–Behnken and factorial or fractional factorial designs), a

A.P. Paiva et al. / Journal of Materials Processing Technology 189 (2007) 26–35

natural choice for the lower and upper bounds on the x’s are the experimental low and high coded levels, respectively. To solve the proposed constrained nonlinear optimization system there are several available methods such as gradientbased, sequential quadratic programming (SQP), genetic algorithms (GA), simulated annealing and ant colony [24,25]. According to K¨oksoy and Doganaksoy [23], the generalized reduced gradient (GRG) is one of the most robust and efficient methods of constrained nonlinear optimization. The expression “reduced gradient” comes from the substitution of the constraints on the objective function, decreasing then the number of variables and consequently reducing the number of present gradients. When making the partition of the original variables into basics (Z) (or dependents) and nonbasics (Y) (or independent), one can write F(X) = F(Z,Y) and h(X) = h(Z,Y). In order to meet the condition of optimality it is necessary that dhj (X) = 0. Making A = z hj (X)T and B = Y hj (X)T , then dY = −B−1 A dZ. Consequently the GRG can be defined as: GR =

d T F (X) = ∇z F (X) − [B−1 A] ∇Y F (X)T dZ


The searching direction is SX = [−GR dY]T . For the iterations one can use Xk+1 = Xk + αSk+1 , verifying at each step if Xk+1 is adequate and h(Xk+1 ) = 0. The final step consists of solving F(X) as a function of α, using a one-dimensional algorithm of search like the Newton method. The main goal of optimization is to find the best possible combination of factors – termed design variables – to maximize a given objective function. Calculus-based solution methods are most commonly used. These methods require that the optimization start from a number of initial points to avoid converging with a local minimum. Thus, calculus-based methods are best when the solution space is convex. For nonconvex solution spaces, the results might not be globally optimal [25]. Duffuaa et al. [24] have compared the results of a number of gradient-based optimization algorithms with different machining models. Their approach has limitations because of its use of gradient-based methods, which are not ideal for nonconvex problems. They have concluded that the generalized reduced gradient method is the most suitable for solving machining optimization models. Treating the principal components rather than the original machining multiple responses has several advantages. If the first principal component represents a high proportion of the total variance in the data, it provides a univariate summary of the multivariate responses. Inspection of the loadings (eigenvectors) reveals the kind of relationship between the ith principal component score equation and the original responses. According to Bratchell [13], a response surface model for principal component provides a model of the overall response, which takes into account the correlations among the response variables and their relative importance. Linear relationships among the response variables can be immediately identified by zero eigenvalues and omitted from further consideration, which avoids unnecessary work when the number of responses exceeds the number of observations. Even though a set of variables is well represented by the principal components, a single principal component is not always


Table 1 Chemical composition of the AISI 52100 steel (%) C
















enough for this representation. Moreover, this does not occur in the majority of the complex manufacturing processes. To represent adequately the responses associated to different principal components it is necessary to add multivariate constraints, written in terms of the smallest principal components. An eigenvector or correlation analysis among responses and principal component scores can reveal which components must be chosen to create the multivariate optimization system. Another necessary aspect in the association of PCA and RSM is related to the type of optimization that the multivariate objective functions must follow. It is worth noting that is the main reason why PCA is more widely used with Taguchi than with response surface designs. The analysis in Taguchi designs is done employing the concept of loss function. Specifically that means each kind of optimization (maximization, minimization or normalization) can be represented by a proper signal-to-noise relation applied to the original responses. Due to the mathematical nature of this relation, the signal-to-noise ratio must always be maximized [2–8]. In RSM, however, the approach is totally different. In surmounting this barrier, this article proposes to find out the kind of optimization analyzing the correlation observed between the significant principal components and each original response. If there is a positive correlation between a principal component score and a specific original response, then they will have the same direction of optimization. If the correlation between the two output variables is negative, the maximization of one variable implies the minimization of the other and vice versa. 4. Experimental procedure To comply with the objectives of this research the work pieces used in the turning process were made with dimensions of Ø 49 mm × 50 mm. All of them were quenched and tempered. After this heat treatment, their hardness was between 53 and 55 HRC, up to a depth of 3 mm below the surface. The work piece material was AISI 52100 steel, with the chemical composition shown in Table 1. The machine tool used was a CNC lathe with a 5.5 kW spindle motor with conventional roller bearings. The mixed ceramic (Al2 O3 + TiC) inserts used were coated with a very thin layer of titanium nitride (TiN) presenting a chamfer on the edges. Their ISO code was CNGA 120408 S01525 and they were made by Sandvik Coromant (Sandvik class CC6050). The tool holder presented negative geometry with ISO code DCLNL 1616H12 and entering angle χr = 95◦ . Tool flank wear measurements (VBmax ) were taken through an optical microscope. The end of life criteria was the breaking of the tool point. Adopting this experimental condition, the work pieces were machined using the range of parameters as defined by Table 2. A sequential set of experimental runs was established using a blocked central composite design built according to the design shown in Table 3. The blocked CCD was adopted in this work because the experiments were done in 2 different weeks, a considerable time elapsed between the running of the two portions of the entire design. To assess the interference of this nuisance variable over the responses, the ANOVA of the blocks was used. In the set of recorded responses, tool life (T), surface roughness (Ra ) and cutting time (Ct ), were observed, while total cost (Kp ), the total turn-


A.P. Paiva et al. / Journal of Materials Processing Technology 189 (2007) 26–35

Table 2 CCD factor levels Parameter


Cutting speed Feed Depth of cut


V f d

Levels (coded)

m/min mm/rev mm






187.34 0.0342 0.1025

200 0.050 0.150

220 0.075 0.225

240 0.100 0.300

252.66 0.1158 0.3475

Table 3 Parameters symbols and values adopted in the study

[26] is determined as follows:




Batch size (units) Secondary time (min) Tool approximation and retreat time (min) Set-up time (min) Insert changing time (min) Machine and labor costs (US$) Tool holder price (US$) Average tool holder life (number of edges) Insert price (US$) Number of cutting edges on the insert Section length Initial diameter (mm) Final diameter (mm) Average diameter (mm)

Z ts ta tp tft Sm + Sh Vsi Nfp Kpi Ns lf D d Dm

1.000 0.5 0.1 60 1 80 200 1.000 50 4 50 49 46 47.5

Tt = (1 + tft )

 l πd   f 1000fVc

+ ts + ta +

tp 1 − tft Z Z


As defined in ref. [26], the total cost of the turning process (Kp ) can be described as: Kp =




1 Z

(Sh + Sm ) +

Ct Ct (Sh + Sm ) + [Kft + tft (Sh + Sm )] (13) 60 T

For interchangeable inserts, the tool cost per life (Kft ) can be represented by: Kft =

Kpi Vsi + Nfp Ns


The symbols used in Eqs. (11)–(14) and respective values adopted in this study are shown in Table 3. Initially, the experimental region (Ω) was considered as the basic constraint for the optimization problem.

5. Results and discussion ing cycle time (Tt ) and material removal rate were calculated using Eqs. (11)–(14). According to Cauchick-Miguel and Coppini [26], the cutting time (Ct ), also called machining time, is the time that the tool actually spends in the feed mode or cutting and removing chips. Mathematically, this variable can be described in cylindrical turning as: Ct =

lf πDm 1000fVc

The RSM model, as described in the above, was applied using an adequate central composite design to collect the data of the six responses, representing the information necessary to build a second-order model, as can be seen in Table 4. Experimentally, two least squares-based models can be obtained: coded and uncoded based parameters. According to some researchers [14,23], the coded units approach is better than uncoded because it provides resources to eliminate any spurious statistical results due to different measurement scales for


where lf is the length of part, Dm the work piece average diameter, f the feed rate and Vc is the adopted cutting speed. The total turning cycle time (Tt ) in minutes

Table 4 Cutting parameters and responses for the CCD n













1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18

1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2

200 240 200 240 200 240 200 240 220 220 187.34 252.66 220 220 220 220 220 220

0.05 0.05 0.1 0.1 0.05 0.05 0.1 0.1 0.075 0.075 0.075 0.075 0.0342 0.1158 0.075 0.075 0.075 0.075

0.15 0.15 0.15 0.15 0.3 0.3 0.3 0.3 0.225 0.225 0.225 0.225 0.225 0.225 0.1025 0.3475 0.225 0.225

16.75 11.50 9.85 8.50 11.50 7.45 8.20 6.25 8.60 6.80 10.10 7.60 17.50 7.20 12.00 6.70 7.20 9.10

7.70 6.41 3.85 3.21 3.85 3.21 1.92 1.60 3.11 3.10 3.65 2.71 6.82 2.01 6.82 2.01 3.09 3.11

8.82 7.63 4.90 4.24 4.84 4.30 2.82 2.52 4.13 4.23 4.67 3.72 7.87 2.95 8.05 2.97 4.20 4.11

17.59 17.26 11.49 10.45 10.71 11.20 6.74 6.62 10.10 11.44 10.82 9.49 15.45 7.49 17.96 7.78 11.09 9.82

1.50 1.80 3.00 3.60 3.00 3.60 6.00 7.20 3.71 3.71 3.16 4.26 1.69 5.73 1.69 5.73 3.71 3.71

0.330 0.280 0.695 0.565 0.245 0.420 0.565 0.610 0.360 0.420 0.340 0.450 0.320 0.720 0.360 0.310 0.370 0.290

4.27 3.01 −0.22 −0.74 0.70 −0.50 −2.50 −3.31 −0.48 −0.63 0.23 −1.18 3.64 −2.70 3.24 −1.97 −0.54 −0.33

−0.59 0.24 −1.79 −0.73 1.10 0.25 −0.41 −0.55 0.63 0.29 0.58 0.18 −0.36 −1.41 −0.40 1.30 0.61 1.07

A.P. Paiva et al. / Journal of Materials Processing Technology 189 (2007) 26–35

the factors. In addition, uncoded units often lead to co-linearity among the terms in the model. This inflates the variability in the coefficients estimates and makes them difficult to interpret. For these reasons, the model based on the coded units was employed. Table 6 presents the coefficients of the model obtained by the employment of the ordinary least square algorithm. Another important aspect in the statistical model-building process is associated with the amount of explicability of the dependent variables y by the predictors x. Table 6 shows the adjusted R-squared. This expression is of the larger-is-better type. A larger adjusted R2 indicates a high degree of explanation of the interest response. The adjusted R2 takes into account the fact that R2 tends to overestimate the actual amount of variation accounted for in the sample analysis. That is, if one applies the regression equation derived from a sample to another independent sample, one almost always gets a smaller R2 in the new sample than in the original. It can be noted that the models found in the present work are adequate once all of them exhibit a large adjusted R2 . The analysis of variance (ANOVA) was usually used to formally test for significance of the main effects and interactions. To refine the model, a common approach consists of removing any nonsignificant term from the full model. As a decision rule, if P-value is lower than 0.05 the corresponding term will be considered significant to the model. Otherwise, if P-value is greater than 0.05 the term will be excluded. According to Montgomery [14], this procedure is convenient in obtaining a simplified model but could decrease the coefficient of determination R2 and increase the error term S. Moreover, the exclusion of any term should follow the hierarchy principle. This is a modelbuilding principle that suggests when a particular polynomial term is included in a model, all lower-order polynomial terms should also included, even those that do not individually exhibit significance. The hierarchy principle promotes an internal consistency in the model. Table 6 shows the complete model for each response. In spite of some nonsignificant terms being found, their exclusion from the complete model increased the error S and reduced R2 (adj.). A complete second-order model was then considered to overcome this problem. The first two principal components shown in Table 4 were obtained by applying Eq. (7) to the response data. Hence, after determining the principal components scores PC1 and PC2 , the regression method was applied to create the related equa-


Fig. 1. First principal component (PC1 ) surface plot in function of f and d.

Fig. 2. First principal component (PC1 ) surface plot in function of d and V.

tions. Figs. 1 and 2 represent the fitted second-order response surfaces for the first principal component score in terms of the cutting parameters. Table 5 shows that the first principal component represents 81.6% of variation in the responses, and sufficient variance–covariance explanation. Moreover, the eigenvectors show a high positive correlation between PC1 and the original responses. Tool life (T), cutting time (Ct ), total turning cycle time (Tt ) and process cost (Kp ), while showing a high negative correlation between PC1 and material removing rate. Although there is a good explanation observed in first principal component, a poor correlation between PC1 and surface roughness (Ra ) is detected. Otherwise, it could be noted that PC2 and Ra have a strong and negative correlation. For the process

Table 5 Principal component analysis of the original turning responses Eigenvalue Proportion Cumulative

4.897 0.816 0.816

0.720 0.120 0.936

0.267 0.044 0.981

0.116 0.019 1.000

0.001 0.000 1.000

0.000 0.000 1.000








T Ct Tt Kp MRR Ra

0.403 0.445 0.446 0.436 −0.424 −0.265

−0.180 −0.160 −0.154 −0.107 −0.018 −0.952

−0.806 −0.028 0.032 0.419 −0.4 0.112

0.233 −0.284 −0.296 −0.344 −0.805 0.105

0.317 −0.548 −0.319 0.697 0.105 0.003

−0.015 −0.627 0.767 −0.135 −0.007 0.000


A.P. Paiva et al. / Journal of Materials Processing Technology 189 (2007) 26–35

Table 6 Full quadratic models for each response Term









Constant V f d V2 f2 d2 Vf Vd fd R2 (adj.) (%)

−0.4758 −0.4569 −1.8452 −1.5328 −0.0430 0.3113 0.3732 0.1413 −0.0288 0.2788 99.20

0.672 0.019 −0.465 0.454 −0.154 −0.626 −0.127 0.116 −0.361 −0.016 85.00

7.9680 −1.2510 −2.3410 −1.6390 0.2340 1.5470 0.4220 0.7500 0.0750 0.6750 85.00

3.1160 −0.3320 −1.3830 −1.3830 −0.0060 0.4570 0.4570 0.1210 0.1210 0.4390 99.10

4.1800 −0.3180 −1.4360 −1.4550 −0.0230 0.4330 0.4700 0.0960 0.1260 0.4390 99.30

10.6220 −0.2380 −2.5840 −2.8610 −0.1960 0.2970 0.8220 −0.1650 0.2180 0.5450 97.20

3.7130 0.3380 1.2380 1.2380 0.0000 0.0000 0.0000 0.1130 0.1130 0.4130 99.90

0.3560 0.0160 0.1360 −0.0080 0.0230 0.0700 0.0000 −0.0260 0.0500 −0.0180 89.10

Table 7 Significance of the coefficients of full quadratic models for first and second principal component scores Term

Constant V f d V2 f2 d2 Vf Vd fd











−0.476 −0.457 −1.845 −1.533 −0.043 0.311 0.373 0.141 −0.029 0.279

0.0985 0.0543 0.0543 0.0543 0.0580 0.0580 0.0580 0.0701 0.0701 0.0701

−4.8300 −8.4200 −34.0030 −28.2460 −0.7420 5.3680 6.4350 2.0160 −0.4100 3.9790

0.0010 0.0000 0.0000 0.0000 0.4790 0.0010 0.0000 0.0790 0.6920 0.0040

0.669 0.019 −0.465 0.453 −0.151 −0.626 −0.125 0.118 −0.360 −0.018

0.1634 0.0900 0.0900 0.0900 0.0962 0.0962 0.0962 0.1162 0.1162 0.1162

4.0960 0.2060 −5.1630 5.0300 −1.5750 −6.5070 −1.3020 1.0110 −3.0990 −0.1510

0.0030 0.8420 0.0001 0.0010 0.1540 0.0000 0.2290 0.3410 0.0150 0.8840

improvement, the responses Ra , Tc , Tt and Kp must be minimized while MRR must be maximized. Examining the eigenvectors in Table 5, it is possible to note that is impossible to maximize the tool life (T) and minimize MRR simultaneously. Probably, the multivariate optimization approach can leads to tool life (T) minimization. By analyzing the correlation between the principal components and the responses (expressed in terms of eigenvectors) it is observed that when minimizing the model built with the first principal component score PC1 , the responses MRR will be maximized while Tc , Tt , Kp and T will be minimized. Table 6 presents the full quadratic models of each response of interest. Tables 7–9 show the significance of the coefficients and the ANOVA table of the full quadratic models of PC1 and PC2 . The well established polynomial of second order can be justified observing the relation between P-value and the significance

level, which in this research was adopted as 0.05. Full quadratic models were used for all responses where there was no lack of fit in this kind of model. By analyzing the regression established between PC2 and Ra , it can be observed that Ra is less than 0.3 ␮m when PC2 assumes values greater than 0.6. To address the entire set of optimization directions, a nonlinear optimization algorithm system can be written in terms of the multivariate response surface models of PC1 and PC2 , as described by the coefficients in Table 7. In this way, using the GRG method, a spherical constraint will be imposed on the factor levels, i.e., the values that optimize the responses of interest must belong to an experimental interval −r ≤ xi ≤ +r. On this case, where was used a blocked central composite design for three factors, the natural choice recommended by Montgomery

Table 8 ANOVA for first principal component PC1

Table 9 ANOVA for second principal component PC2



Regression Linear Square Interaction

9 3 3 3

Residual error Lack-of-fit Pure error

8 5 3








82.912 79.511 2.614 0.788

9.212 26.503 0.871 0.263

234.62 674.97 22.19 6.69

0.000 0.000 0.000 0.014

Regression Linear Square Interaction



Residual error Lack-of-fit Pure error

0.3141 0.266 0.0477 83.226

0.039 0.0533 0.0159







9 3 3 3

11.356 5.615 4.592 1.149

1.262 1.872 1.530 0.383

11.69 17.33 14.17 3.55

0.001 0.001 0.001 0.067

8 5 3

0.864 0.556 0.308

0.108 0.111 0.103





A.P. Paiva et al. / Journal of Materials Processing Technology 189 (2007) 26–35


Table 10 Multivariate optimization results Cutting parameters Multivariate optimization V Coded units Uncoded units


Multiple optimization f

0.9127 238.25



0.3202 0.083


1.3158 0.3237

V (m/min)a

f (mm/rev)a

d (mm)a

0.4001 228.0010

0.3679 0.0842

1.5400 0.3405

Optimal values of responses Kp (US$ per piece) T (min) Ct (min) Tt (min) MRR (cm3 /min) Ra (mm) a

7.284 5.637 1.602 2.567 6.251 0.415

7.410 6.000 1.750 2.715 6.458 0.400

Units (SI).

[14] for the radius is 1.633. Minimize :

(8 )


Subject to :

xT x ≤

k(1 + (ns0 /ns )) (1 + (nc0 /nc ))

1/2  ; PC2 ≥ 0.6 (9 )

To compare the results obtained with the multivariate response surface optimization, an alternative constrained and nonlinear optimization procedure can be proposed. For this approach the following system of equations is suggested. If the multivariate optimization is adequate, similar results can be expected. Using the fitted equations of the machining outputs shown in Table 6 as constraints and choosing the cost per piece (Kp ) as an objective function, the optimization system can be written as: Min :



xT x ≤ ρ 2 ;

S.T. :

T ≥ 6;

MRR ≥ 6;

Ra ≤ 0.4

Ct ≤ 2;

Tt ≤ 3; (16)

The constraint in Eq. (16) represents the desired values for the decision variables. In this equation, xT x ≤ ρ2 represents an experimental hyperesphere. Table 10 shows the results obtained using the multivariate multiresponse optimization and traditional constrained nonlinear optimization approaches. It is possible to note that the multivariate approach presented results slightly different of those obtained with the traditional optimization, although the process cost have been considering more appropriated. However, in order to optimize the original set of responses, the minimization of PC1 provided a larger value for cutting speed (V) than the traditional approach. A reasonable explanation to this event could be stated in terms of the significance of the cutting parameters to each response. While V, f and d are equally important to explain the models of PC1 , T, Ct , Tt and MRR, for Kp and Ra the cutting speed is not significative. Then, when the objective function is Kp , the GRG algorithm tries to achieve the ceramic tool life constraint (T ≥ 6 min), reducing in

this case the value of the cutting speed to obtain the desired target. Otherwise, when Kp is replaced by PC1 in the optimization system, the cutting speed must be increased because the minimization of PC1 implies the optimization of the set of responses. Note that the coefficients for the cutting speed parameter are negatives in all models (Table 6). The same fact can be observed with MRR. For the maximization of the material removal rate, a larger value of V produces the expected optimization. In this model, the coefficient of V is positive. Once the results are compatible with the expected values and the hard turning theory, the multivariate method may be considered suitable for improving the machining process, mainly when a large set of correlated responses are employed. Although in this case the first principal component alone was not enough to represent an adequate optimization set, an auxiliary multivariate constraint based on the second principal component was used. Forcing the solution to fall within the experimental region through the means of a spherical constraint also leads the optimization problem to a feasible solution (Table 11). 6. Sensitivity analysis of the turning parameters To assess the sensitivity of the turning parameters in an optimal condition, there are two possible strategies: (a) change the constraint values (ρ2 , PC2 ) or (b) change the optimum value x* = (V* , f* , d* ). The adequate approach to develop the strategy (a) is to use the concept of Lagrange multipliers. According to Nash and Sofer [27], the Lagrange multipliers express the gradient at the optimum as a linear combination of the rows of the constraint matrix and can indicate the sensitivity of the optimal objective value to changes in the data. Assuming that the objective function is twice continuously differentiable, and considering small perturbations (δ) in the right side of the constraints it is possible to use a Taylor series to obtain the approximation: f (¯x) = f (x∗ ) +


δi λ∗i



where x* represents a local minimum, such as x∗ = arg minf (x). x In particular, Eq. (17) is valid if x¯ is the minimizer of the


A.P. Paiva et al. / Journal of Materials Processing Technology 189 (2007) 26–35

Table 11 Sensitivity analysis of the turning parameters

Optimal −2% +2% +δ −δ

Optimal −2% +2% +δ −δ








−2.733 −2.489 −2.971 −2.808 −2.533

5.637 5.807 5.543 5.461 5.928

1.602 1.712 1.509 1.589 1.694

2.567 2.691 2.457 2.552 2.667

7.285 7.597 6.961 7.286 7.484

6.252 5.930 6.586 6.381 5.922

0.415 0.397 0.436 0.410 0.424




g1 (x)

g2 (x)

λ1(xt x)

λ2(PC2 )

238.250 233.577 243.014 243.139 231.666

0.083 0.081 0.085 0.080 0.086

0.324 0.317 0.330 0.341 0.301

2.667 2.041 3.432 3.767 1.567

0.600 0.720 0.450 0.700 0.500



perturbed problem. If the right-hand side of the ith constraint changes by δi then the optimal objective value changes by approximately δi λ*i . Hence, λ*i represents the change in the optimal objective per unit change in the ith right-hand side. For this reason, the Lagrange multpliers are also called shadow prices [27]. To write a nonlinear constrained optimization problem in terms of an unconstrained form and considering the inequality constraints, it is plausible to consider the Lagragian function as follow [28]:  p  n   L(xi , λj , sj ) = f (xi ) − λj aji xi − bj j=1

− p






aji xi − bj − sj2



where i=1 aji xi − bj are the constraints established in the optimization problem and sj are the slack variables used to transform the inequalities in equality constraints. To solve Eq. (18) it is necessary to use the Karush– Kuhn–Tucker conditions (KKT) [26,27], what can be done using the gradient of the Lagrangian function. Then, upon elimination of sj , it is possible to write: ∂L(X, λ) ∂L(X, λ) = 0, =0 (19) ∂X ∂λ The Lagrange multipliers can be determined using the optimization routine available in the Microsoft Excel Solver. The two Lagrange multipliers obtained with the present data were λ1 = −0.1868 (referring to the first constraint) and λ2 = −0.8188. Substituting these values in Eq. (17), it follows that: PC1 (¯x) = PC1 (x∗ ) +


δi λ∗i


= −2.733 − 0.1868δ(xT x) + 0.8188δ(PC2 )


Adopting δ(xT x) = 1.1 and δ(PC2 ) = 0.1 as the values of respective perturbations associated to the constraints, the new optimal value for PC1 is −2.808. Using the opposite values of the considered perturbations, PC1 becomes −2.533. The values assumed

by the original set of responses are presented in Table 8. The perturbations values are different because the constraints are not represented in the same scale. In order to indicate the range of the parameters which would give the same performance, another possible analysis can be done changing the uncoded parameters levels. Table 8 also describes the new values assumed by the original responses when the optimum is changed in ±2%. In general, this modification produces 9% of variation in the PC1 value, what corresponds approximately to 5% of variation in the most of the original optimal values. 7. Conclusions This research shows that by using a multivariate approach a multiresponse manufacturing process with a moderate to high correlation structure can be adequately represented. In the present study the results indicate that the approach is adequate in situations where the multiple responses exhibit some correlation. The first principal component is responsible for most of the variance–covariance present in the original data. Because of this, a regression equation built using ordinary least square method over the score of these latent variables seemed a good way to reduce the dimensionality of the response set and to form a singular and representative objective function. It was observed that the multivariate constraint based on the second principal component and a spherical constraint was efficient in promoting the adequate optimum for the problem. This paper also shows that the results obtained from the experimental investigation over the influence of the cutting parameters are consistent with the theory. In this case, to maximize the material removal rate while minimizing the cutting times, costs and surface quality simultaneously, a cutting speed of approximately 238 m/min, with feed rate of 0.08 mm/rev and a depth of cut of 0.32 mm are the values that attained the desired quality conditions. A variation of ±2% in these values does not modify significantly the multiresponse optimum. The negative correlation observed between PC1 and the ceramic tool life (T), however, lead to the minimization of this response. It was possible to state that the response surface methodology combined with principal component analysis are very useful techniques when modeling and creating equations

A.P. Paiva et al. / Journal of Materials Processing Technology 189 (2007) 26–35

to forecast and to optimize any dependent variable from the machining process. Regarding other approaches, a small number of experiments were necessary to generate a useful data pool for those interested in the turning process of the AISI 52100 hardened steel using mixed ceramic tools. For the six responses of interest, a second-order model proved more appropriate. These conclusions cannot be extrapolated to different materials, tools, or machine tools and they are valid only in the adopted range level. It can, nonetheless, be recommended to fit the methodology of any production processes. References [1] F.P. Chang, K.T. Chiang, Optimization of the WEDM process of particle-reinforced material with multiple performance characteristics using grey relational analysis, J. Mater. Process. Technol. 180 (2006) 96– 101. [2] T. Yih-Fong, A hybrid approach to optimize multiple performance characteristics of high-speed computerized numerical control milling tool steels, Mater. Des. 28 (2007) 36–46. [3] C.P. Fung, P.C. Kang, Multi-response optimization in friction properties of PBT composites using Taguchi method and principal component analysis, J. Mater. Process. Technol. 170 (3) (2005) 602–610. [4] Y.A. Song, S. Park, S.W. Chae, 3D welding and milling. Part II. Optimization of the 3D welding process using experimental design approach, Int. J. Mach. Tools Manuf. 45 (2005) 1063–1069. [5] J.L. Lin, K.S. Wang, Y.S. Tarng, Optimization of the electrical discharge machining process based on the Taguchi method with fuzzy logics, J. Mater. Process. Technol. 102 (3) (2000) 48–55. [6] H.C. Liao, Multi-response optimization using weighted principal components, Int. J. Adv. Manuf. Technol. 27 (7) (2005) 720–725. [7] H.C. Liao, Chen F.Y.K., Optimizing multi-response problem in the Taguchi method by DEA based ranking method, Int. J. Qual. Reliab. Manage. 19 (7) (2002) 825–837. [8] C.Y. Nian, W.W. Yang, Y.S. Tarng, Optimization of turning operations with multiple performance characteristics, J. Mater. Process. Technol. 95 (1999) 90–96. [9] M.S. Kulkarni, V. Mariappan, Multiple response optimization for improved machined surface quality, J. Mater. Process. Technol. 141 (2003) 174–180. [10] G.E.P. Box, W.G. Hunter, J.F. MacGregor, J. Erjavec, Some problems associated with the analysis of multiresponse data, Technometrics 15 (1) (1973) 33–51. [11] G. Derringer, R. Suich, Simultaneous optimization of several response variables, J. Qual. Technol. 12 (4) (1980) 214–219.


[12] A.I. Khuri, M. Conlon, Simultaneous optimization of multiple responses represented by polynomial regression functions, Technometrics 23 (4) (1981) 363–375. [13] N. Bratchell, Multivariate response surface modeling by principal components analysis, J. Chemom. 3 (1989) (1989) 579–588. [14] D.C. Montgomery, Design and Analysis of Experiments, fourth ed., Wiley, New York, 2001. [15] I.A. Choudhury, M.A. El-Baradie, Tool-life prediction model by design of experiments for turning high strength steel (290 BHN), J. Mater. Process. Technol. 77 (1998) 319–326. [16] Y. Beauchamp, M. Thomas, A.Y. Yossef, J. Masounave, Investigation of cutting parameter effects on surface roughness in lathe boring operation by use of a full factorial design, Comp. Ind. Eng. 31 (3/4) (1996) 645–651. [17] M.Y. Noordin, V.C. Venkatesh, S. Sharif, S. Elting, A. Abdullah, Application of response surface methodology in describing the performance of coated carbide tools when turning AISI 1045 steel, J. Mater. Process. Technol. 145 (1.) (2004). [18] M. Alauddin, M.A. El-Baradie, M.S.J. Hashmi, Prediction of tool life in end milling by response surface methodology, J. Mater. Process. Technol. 71 (1997) 457–465. [19] S.K. Choudhury, G. Bartarya, Role of temperature and surface finish in prediction tool wear using neural network and design of experiments, Int. J. Mach. Tools Manuf. 43 (2003) 747–753. [20] M.N. Dhavlikar, M.S. Kulkarni, V. Mariappan, Combined Taguchi and dual response method for optimization of a centerless grinding operation, J. Mater. Process. Technol. 132 (2003) 90–94. [21] P.G. Benardos, G.C. Vosniakos, Prediction of surface roughness in CNC face milling using neural networks and Taguchi’s design of experiments, Robot. Comp. Integ. Manuf. 18 (2002) 343–354. [22] R.A. Johnson, D. Wichern, Applied Multivariate Statistical Analysis, Prentice-Hall Inc., Englewood Cliffs, New Jersey, 2002. [23] O. K¨oksoy, N. Doganaksoy, Joint optimization of mean and standard deviation using response surface methods, J. Qual. Technol. 35 (3) (2003) 237–334. [24] S.O. Duffuaa, A.N. Shuaib, A. Alam, Evaluation of optimization methods for machining economics models, Comp. Oper. Res. 20 (1992) 227–237. [25] Z. Khan, B. Prasad, T. Singh, Machining conditions optimization by genetic algorithms and simulated annealing, Comp. Oper. Res. 24 (7) (1997) 647–657. [26] P.A. Cauchick-Miguel, N.L. Coppini, Cost per piece determination in machining process: an alternative approach, Int. J. Mach. Tools Manuf. 36 (8) (1996) 939–946. [27] S.G. Nash, A. Sofer, Linear and Nonlinear Programming, first ed., McGraw-Hill Companies Inc., 1996, p. 692. [28] W.S. Dorn, On Lagrange multipliers and inequalities, Oper. Res. 9 (1) (1961) 95–104.