Metodoloski zvezki, Vol. 12, No. 1, 2015, 11-24 Asymptotic Normality of the Optimal Solution in Multiresponse Surface Mathematical Programming Jose A. Díaz-García1 Francisco J. Caro-Lopera2 Abstract An explicit form for the perturbation effect on the matrix of regression coefficients on the optimal solution in multiresponse surface methodology is obtained in this paper. Then, the sensitivity analysis of the optimal solution is studied and the critical point characterisation of the convex program, associated with the optimum of a multiresponse surface, is also analysed. Finally, the asymptotic normality of the optimal solution is derived by the standard methods. 1 Introduction The multiesponse surface methodology explores the relationships among several explanatory variables and more than one response variables. The addressed methodology considers a sequence of designed experiments in order to obtain a simultaneous optimal response. To reach this aim the method uses a second-degree polynomial model for each response variable. With that constraint the technique is just an approximation, but they are succeed in literature because such models can be easily interpreted, estimated and applied; moreover they perform well under the usual uncertainty about the process or phenomenon under consideration. In fact, a number of laws in sciences are usually explained with second-degree polynomial models, given that the first and second corresponding flows are well understood and they explain some crucial intrinsic property of the phenomenon. In recent decades, the multiresponse surface methodology has attracted the attention of many quality engineers in different industries. Quality improvement or optimization for such process or phenomenons, needs precise identification of the operation stages and their effects on the response variables. Therefore, multistage systems require special methods and solutions, since applying uni-response surface techniques may lead to suboptimal or even inaccurate results. Some examples of the mentioned industries are: agricultural, pharmaceutical, chemical, assembly, semiconductor, textile, and petroleum industries as well as queuing, healthcare, traffic control, and transportation systems. Start by assuming that a researcher knows a process or phenomenon and a corresponding set of observable responses variables Y1, ■ ■ ■ ,Yr which depends on some input Universidad Autónoma Agraria Antonio Narro, Calzada Antonio Narro 1923, Col. Buenavista, 25315 Saltillo, Coahuila, Mexico, jadiaz@uaaan.mx 2 Department of Basic Sciences, Universidad de Medellín, Medellín, Colombia, fjcaro@udem.edu.co 12 J. A. Diaz-Garcia and F. J. Caro-Lopera variables, x\,... xn. Suppose also that the input variables xi's can be controlled by the researcher with a minimum error. Typically, we have that Yk (x) = rqk (xi ,...xn), k = 1,..., r, and x = (xi,... Xn)', (1.1) where the form of the functions nk(-)'s, usually termed as the true response surface, are unknown and perhaps, very complex. The success of the response surface methodology depends on the approximation of nk(■) by a polynomial of low degree in some particular region. As it was mentioned, in the context of this paper we will assume that nk(■) can be soundly approximated by a polynomial of second order, that is n n n n Yk ( x) — Pok I ^ ^ fiikxi I ^ ^ ftiikxi + y] y /ijkxixj (1.2) i=1 i=1 i=1 j>i where the unknown parameters /js can be estimated via regression's techniques, as it shall be described in the subsequent section. We also study the levels of the input variables xi's such that the response variables Y1, ■ ■ ■ ,Yr are simultaneously minimal (optimal). This can be achieved if the following multiobjetive mathematical program is solved ( Y1(x) ^ Y2(x) mm x (1.3) V Yr (x) y subject to x E X, where X is certain operating region for the input variables xi's. Now, two questions, closely related, can be observed: 1. When the estimations of (1.2) for k = 1,..., r are considered into (1.3), the critical point x* obtained as solution shall be a function of the estimators /j s of the j. Thus, given that/j are random variables, then x* = x*(/3j's) is a random vector too. So, under the assumption that the distribution of 3 is known, then, we ask for the distribution of x* (/j's). 2. And, given that a point estimate of x*(3j's) should not be sufficient, then we would ask also for an estimated region or an estimated interval. In particular, the distribution of the critical point in a univariate response surface model was studied by (Diaz Garcia and Ramos-Quiroga, 2001, 2002), when y(x) is defined as an hyperplane. Now, in the context of the mathematical programming problems, the sensitivity analysis studies the effect of small perturbations in: (1) the parameters on the optimal objective Multiresponse Surface Mathematical Programming 13 function value and (2) the critical point. In general, these parameters shape the objective function and constraint the approach to the mathematical programming problem. In particular, (Jagannathan, 1977; Dupacova, 1984; Fiacco and Ghaemi, 1982) have studied the sensitivity analysis of the mathematical programming, among many other authors. As an immediate consequence of the sensitivity analysis, the corresponding asymptotic normality study of the critical point emerges naturally, which can be performed by standard methods of mathematical statistics (see similar results for the case of maximum likelihood estimates in (Aitchison and Silvey, 1958)). This last consequence makes the sensitivity analysis as an interesting source of statistical research. However, this approach must be fitted into the classical philosophy of the sensitivity analysis; i.e., we need to translate the general sensitivity analysis methodology into the statistical language. This involves, for example, to study the effect on the model estimators of adding and/or excluding variables and/or observations, see (Chatterjee and Hadi, 1988). This paper proposes a solution in order to establish the effect of perturbations of the matrix of regression parameters on the optimal solution of the multiresponse surface model and the asymptotic normality of the critical point. First, in Section 2 some notation is proposed. Then, the multiresponse surface mathematical program is set in Section 3 as a multiobjective mathematical programming problem and a general solution is considered in terms of a functional. Then, the characterisation of the critical point is given in Section 4 by stating the first-order and second-order Kuhn-Tucker conditions. Finally, the asymptotic normality of a critical point is established in Section 5 and for a particular form of the functional, the asymptotic normality of a critical point is also derived. 2 Notation For the sake of completeness, the main properties and usual notations are given here. But for a detailed discussion of the multiresponse surface methodology we recommend references (Khuri and Cornell, 1987; Khuri and Conlon, 1981, Chap. 7). Let N be the number of experimental runs and r be the number of response variables, which can be measured for each setting of a group of n coded variables x\,x2,... ,xn. We assume that the response variables can be modeled by a second order polynomial regression model in terms of Xj, i = 1,... ,n. Hence, the kth response model can be written as Yfc = Xfcpk + ek, k = 1,..., r, (2.1) where Yk is an N x 1 vector of observations on the kth response, Xk is an N x p matrix of rank p termed the design or regression matrix, p = 1 + n + n(n + 1)/2, fik is a p x 1 vector of unknown constant parameters, and ek is a random error vector associated with the kth response. For purposes of this study, it is assumed that X1 = • • • = Xr = X. Therefore, (2.1) can be written as Y = Xb + e (2.2) where Y Yl.Y2.---.Yr ,b A.^2. ••• moreover = (Pok, Pi k, . . . , Pnk, fill k, . . . , Pnnk, Pl2k, . . . , P(n-1)nk)' 14 J. A. Diaz-Garcia and F. J. Caro-Lopera and e = ei:e2. ■ ■ ■ .er , such that e — xr(0, IN ® S) i.e. e has an N x r matrix multivariate normal distribution with E(e) == 0 and Cov(vec e') = IN ® S, where S is a r x r positive definite matrix. Now, if A AilAo. • • • :Ar with Aj, j = 1, ■ ■ ■ , r the columns of A; then vec A = (Ai, A'2,..., A'r)/ and ® denotes the direct (or Kronecker) product of matrices, see (Muirhead, 1982, Theorem 3.2.2, p. 79). In addition denote • x = (xi , X2, ... , %n ) : The vector of controllable variables or factors. ß i.ß 2. ••• .ß r : The least squares estimator of B given by from where (X'X)-1X/Y, Y k = (X'X) 1X'Yfc = Wok ,filk, . . . , Pn, Pllk, . . . , Punk, @12k, . . . , Y(n_1)nk)' k = 1, 2,..., r. Moreover, under the assumption that e — xr(0, IN ® S), then b - Npxr(b, (X'X)-1 ® S), with Cov(vecb') = (X'X)-1 ® S. Z— (1, X1, X2, . . . , Xn, X1, X 2, . . . , Xn, X1X 2,X1X3 . . . , Xn_1 Xn ) . ßik = (ßik,..., ßnk)/ and B k ( 2ßi ik ßi 2k ß2 ik 2 ß22k \ ßnik ßn2k ßink \ ß2nk 2 ßnnk / Yk (x) = z'(x)3 k n n n n = p0k + y ] PikXi + y ] piikXi + ^ ] ^^ pijkXiXj i=1 i=1 i=1 j>i = p0k + /31kx + x' Bkx : The response surface or predictor equation at the point x for the kth response variable. For the aim of this paper it is assumed that Bk, k = 1,..., r, are positive definite matrices. Note that this last assumption is not always true and should be verified, for example via the canonical analysis, see (Khuri and Cornell, 1987, Subsection 5.5.1, pp. 180-186). 2 Y(x) = ^Y1(x), p2(x),..., Yr(x)J = b'z(x): The multiresponse surface or predicted response vector at the point x. 2Observe that, alternatively can be assumed that Bk, k = 1,..., r, are negative definite matrices and then in equation (3.1) the minimization should be replaced by maximization. 1 2 Multiresponse Surface Mathematical Programming 15 - Y'(In - X(X'X)-1X/)Y S = -—-: The estimator of the variance-covariance matrix N — p S such that (N — p)S has a Wishart distribution with (N — p) degrees of freedom and the parameter S; this fact is denoted as (N — p)S ~ Wr(N — p, S). Here, Im denotes an identity matrix of order m. Finally, notice that E (Y (x)) = E (b/z(x)) = b/z(x) (2.3) and Cov(Y (x)) = z/(x)(X/X)-1z(x)S. (2.4) An unbiased estimator of Cov(Y(x)) is given by CoV(Y (x)) = z/(x)(X/X)-1z(x)S. (2.5) 3 Multiresponse surface mathematical programming In the following sections, we make use of the multiresponse mathematical programming and multiobjective mathematical programming. For convenience, the concepts and notations required are listed below in terms of the estimated model of multiresponse surface mathematical programming, but further detailed properties can be found in (Khuri and Conlon, 1981; Khuri and Cornell, 1987; Rios etal., 1989; Steuer, 1986; Miettinen, 1999). The multiresponse mathematical programming or multiresponse optimisation (MRO) problem is proposed, in general, as follows min Y (x) = min I £(x) A £(x) V Y (x) ) (3.1) subject to x E X. It is a nonlinear multiobjective mathematical programming problem, see (Steuer, 1986; Ríos et al., 1989; Miettinen, 1999); and X denotes the experimental region, usually taken as a hypersphere X = {x|x'x < c2,c E where, c is set according to the experimental design model under consideration, see (Khuri and Cornell, 1987). Alternatively, the experimental region can be taken as a hypercube X = {x|/j < Xi 0, k = 1,... ,r, with y^k=1 Wk = 1, see (Rios et al., 1989). And observe that, as 13 k, k = 1,... ,r, are positive definite matrices, then Y^k=1 13k is a positive definite matrix too. Then the equivalent nonlinear scalar mathematical program defined in this manner is a quadratic program, and hence a convex program, see (Rao, 1979, p. 662). □ Let x* (iB) e Kn be the unique optimal solution of program (3.3) with the corresponding Lagrange multiplier A* (b) e K. The Lagrangian is defined by L(x, A; f Y(x) + A(||x||2 - c2). (4.1) Similarly, x*(b) G Kn denotes the unique optimal solution of program (1.3) with the corresponding Lagrange multiplier A* (b) G K. Now we establish the local Kuhn-Tucker conditions that guarantee that the Kuhn- Tucker point r* x* , A* G Kn+1 is a unique global minimum of convex df program (3.3). First recall that for f : Kn ^ K, — = Vx denotes the gradient of dx function f. Theorem 1. The necessary and sufficient conditions that a point x* (b ) G Kn for arbitrary fixed b G Kp, be a unique global minimum of the convex program (3.3) is that, x* (b) and the corresponding Lagrange multiplier A*(b) G K, fulfill the Kuhn-Tucker first order conditions (4.2) (4.3) (4.4) (4.5) with respect VxL(x, A; b ) = Vxf (Y (x) J + 2A(b )x = 0 VAL(x, A; b) = ||x||2 - c2 < 0 ~ |x||2 - c2) = 0 0 In addition, assume that strict complementarity slackness holds at x* to A*(b), that is A*(b) > 0 ^ ||x||2 - c2 = 0. Analogously, the Khun-Tucker condition (4.2) to (4.5) for b = b are stated next. (4.6) Multiresponse Surface Mathematical Programming 19 Corollary 1. The necessary and sufficient conditions that a point x*(b) e for arbitrary fixed b e be a unique global minimum of the convex program (3.2) is that, x* (b) and the corresponding Lagrange multiplier A* (b) e K, fulfill the Kuhn-Tuckerfirst order conditions VxL(x, A; b) = Vx/ (Y(x)) + 2A(b)x 0 VAL(x,A; \x\I2- c2 < 0 x||2 - c2) 0 (4.7) (4.8) (4.9) (4.10) and A(b) = 0 when ||x||2 - c2 < 0 at [x*(b),A* Observe that, due to the strict convexity of the constraint and objective function, the second-order sufficient condition is evidently fulfilled for the convex program (3.3). The next result states the existence of a once continuously differentiable solution to program (3.3), see (Fiacco and Ghaemi, 1982). Theorem 2. Assume that (4.6) holds and the second-order sufficient condition is satisfied by the convex program (3.3). Then is a unique global minimum of program (3.2) and A*(b) is also unique. 1. x* 2. For b e k(b) (is an £—neighborhood or open ball), there exist a unique once continuously differentiable vector function x* A* e 3fT+1 satisfying the second order sufficient conditions of problem (3.2), such that r* [x*(b), A*(b)]' andhence, x*(b) is a unique global minimum of problem (3.3) with associated unique Lagrange multiplier A* (iB). 3. For b e k(b), the status of the constraint is unchanged and A*(b) > 0 ^ ||x||2 — c2 = 0 holds. 0 * r 5 Asymptotic normality of the critical point This section considers the statistical and mathematical programming aspects in the sensitivity analysis of the optimum of a estimated multiresponse surface model. Theorem 3. Assume that: 1. For any 11$ e k(b), the second-order sufficient condition is fulfilled for the convex program (3.3) such that the second order derivatives d2L(x, A; b) d2L(x, A; b) d2L(x, A; b) d2L(x, A; b) d2L(x, A; b) dxdx' , dxd vec' b ' dxdA , dAdx' , dAd vec b 20 J. A. Diaz-Garcia and F. J. Caro-Lopera exist and are continuous in x*(b), À*(b) G Ve([x*(b), À*(b)]') and d2L(x, À; b) dxdx' is positive definite. 2. ibv, the estimator of the true parameter vector bv, is based on a sample of size Nv such that tnv(bv - bv) - npxr(b, ©), n© = (x'x)-1 0 s. N v 3. (4.6) holds for b = b. Then asymptotically vN x*(b ) - x* A Nn(0n, S), where the n x n variance-covariance matrix d x* d vec] © d x* 1 d vec B/ N —© = (X'X)-1 ® S such that all elements of ( dx furthermore vec b ) are continuous on any iB G V d x* d vec] [I - P-1Q(Q'P-1Q)-1 Q'] P-1G where P Q d2L(x, À; ] d xd x' d2L(x, À;~ dÀdx d2L(x, À; G dxd vec' 11$ Proof. According to Theorem 1 and Corollary 1, the Kuhn-Tucker conditions (4.2)-(4.5) at x* (b), À* (b) and the conditions (4.7)-(4.10) at [x* (b), À* (b)]' are fulfilled for mathematical programs (3.2) and (3.3), respectively. From conditions (4.7)-(4.10) of Corollary 1, the following system equation VxL(x,À; b) = Vx/ Y(x) + x VAL(x,À;B) = ||x||2 - c2 0 (5.1) (5.2) 0 Multiresponse Surface Mathematical Programming 21 has a solution x* (b), A*(b) > 0, b. The nonsingular Jacobian matrix of the continuously differentiable functions (5.1) and (5.2) with respect to x and A at x* A* is ( d2L(x, A; B) d2L(x, A; B) \ V dxdx' d2L(x, A;. d x'dA dAdx 0 / According to the implicit functions theorem, there is a neighborhood V^(b) such that for arbitrary b e V(b), the system (5.1) and (5.2) has a unique solution x*(b), A*(b), b and by Theorem 2, the components of x*(b), A*(b) are continuously differentiable function of b, see (Bigelow and Shapiro, 1974). Their derivatives are given by ( d. x* \ d vec dA* ' V d ( d2L(x, A; b) d2L(x, A; b) \ vec J d xd x' d2L(x, A;: V d x'dA dAd x 0 d2L(x, A;] dxd vec' ] 0 (.5.3) The explicit form of (dx* (b)/d vec b) follows from (5.3) and by the formula P Q Q' 0 i [I - P-1Q(Q'P-1 Q)-1Q']P-1 P-1Q(Q'P-1Q)-1 (Q'P-iQ)-iQ'P-i -(Q'P-iQ)-i where P is symmetric and nonsingular. Then from assumption 2, (Rao, 1973, (iii), p. 388) and (Bishop etal., 1991, Theorem 14.6-2, p. 493) (see also (Cramer, 1946, p. 353)) we have x x Nn[ 0ra, dx* d vec] © d x* d vec] (5.4) Finally note that all elements of (Sx*/db) are continuous on V^(b), so that the asymptotical distribution (5.4) can be substituted by vN x x d . , i n i dx 0n, d vec] © 5x* d vec] see (Rao, 1973, (iv), pp.388-389). As a particular case, assume that the functional in (3.3) is defined as □ f (Y(x))=£wfcffc(x), = !, with wk known constants. Then, k=1 k=1 * 22 J. A. Diaz-Garcia and F. J. Caro-Lopera Corollary 2. Suppose the hypothesis 1 to 3 of Theorem 3 are fulfilled. Then asymptotically vN x (b) — x A Nn(0n, S) where the n x n variance-covariance matrix d x* d vec] © d x* such that all elements of ( dx* more dx* where S d vec b / \vx*(b)'S-1x* d2L(x, A; f , -1 © = (X'X)-1 0 s d vec by N vecb) are continuous on any b e V£(B); further- - In I M (x x*(b)x* (b )'S-1 S dxdx' and M(x) = Vxz'(x) 2 J] wk B k - 2A*(b)In k= 1 d z'(x) d x (0:I„:2 diag(x):C1: ••• iC„_1) e»nxp, with ( 01 \ Ci ,i = 1,... ,n - 1, 0j e ^n-i,j = 1,...,i - 1; 0U x'Ai y xiIn_i observing that when i = 1 (i.e. j = 0), this row does not appear in C1; and 1 01 \ Ai 0i I n_ i 0k G Wl-i,k = 1,...,i. Proof. The required result follows from Theorem 3 and observing that in this particular case M(x) ^ wk0k + 2A(b)x VxL(x, A; b) = < k=1 or ^ Wk [^1k + 2Bkx] + 2A(b)x k=1 V\L(x, A; 0) = ||x||2 - c2 = 0 □ 0 Multiresponse Surface Mathematical Programming 23 Conclusions As the reader can check, the results of the paper can be computed easily from the estimates of the parameters obtained through multivariate multiple regression and a known explicit form of the functional f (•). A few basic routines in software R or MATLAB shall be sufficient for achievement this objective. In addition, as a consequence of Theorem 2 now is feasible to establish confidence regions and intervals and hypothesis tests on the critical point, see (Bishop et al., 1991, Section 14.6.4, pp. 498-500); it is also possible to identify operating conditions as regions or intervals, instead of isolated points. The results of this paper can be taken as a good first approximation to the exact problem. However, in some applications the number of observations can be relatively small and perhaps the results obtained in this work should be applied with caution. Acknowledgements The authors wish to thank the Editor and the anonymous reviewers for their constructive comments on the preliminary version of this paper. This paper was written during J. A. Díaz-García's stay as a professor at the Department of Statistics and O. R. of the University of Granada, España. F. J. Caro-Lopera was supported by project No. 158 of University of Medellin. References [1] Aitchison, J. and S. D. Silvey, S. D. (1958): Maximum likelihood estimation of parameters subject to restraints. Annals of Mathematical Statistics, 29, 813-828. [2] Biles, W. E. (1975): A response surface method for experimental optimization of multi-response process. Industrial & Engeneering Chemistry Process Design Development, 14, 152-158. [3] Gigelow, J. H. and Shapiro, N. Z. (1974): Implicit function theorem for mathematical programming and for systems of iniqualities. Mathematical Programming, 6(2), 141156. [4] Bishop, Y. M. M., Finberg, S. E. and Holland, P. W. (1991): Discrete Multivariate Analysis: Theory and Practice. The MIT press, Cambridge. [5] Chatterjee, S. and Hadi, A. S. (1988): Sensitivity Analysis in Linear Regression. John Wiley: New York. [6] Cramer, H. (1946): Mathematical Methods of Statistics. Princeton University Press, Princeton. [7] Díaz García, J. A. and Ramos-Quiroga, R. (2001): An approach to optimization in response surfaces. Communication in Statatistics, Part A- Theory and Methods, 30, 827-835. 24 J. A. Diaz-Garcia and F. J. Caro-Lopera [8] Díaz Garcia, J. A. and Ramos-Quiroga, R. (2002): Erratum. An approach to optimization in response surfaces. Communication in Statatistics, Part A- Theory and Methods, 31, 161. [9] Dupacova, J. (1984): Stability in stochastic programming with recourse-estimated parameters. Mathematical Programming, 28, 72-83. [10] Fiacco, A. V. and Ghaemi, A. (1982): Sensitivity analysis of a nonlinear structural design problem. Computers & Operations Research, 9(1), 29-55. [11] Jagannathan, R. (1977): Minimax procedure for a class of linear programs under uncertainty. Operations Research, 25, 173-177. [12] Kazemzadeh, R. B., Bashiri, M., Atkinson, A. C. and Noorossana, R. (2008): A General Framework for Multiresponse Optimization Problems Based on Goal Programming. European Journal of Operational Research, 189, 421-429. [13] Khuri, A. I. and Conlon, M. (1981): Simultaneous optimization of multiple responses represented by polynomial regression functions. Technometrics, 23, 363-375. [14] Khuri, A. I. and Cornell, J. A. (1987): Response Surfaces: Designs and Analysis. Marcel Dekker, Inc., NewYork. [15] Miettinen, K. M. (1999): Non linear multiobjective optimization. Kluwer Academic Publishers, Boston. [16] Muirhead, R. J. (1982): Aspects of multivariate statistical theory. Wiley Series in Probability and Mathematical Statistics, John Wiley & Sons, Inc., 1982. [17] Myers, R. H., Montgomery, D. C. and Anderson-Cook, C. M. (2009): Response surface methodology: process and product optimization using designed experiments. Third edition, Wiley, New York,. [18] Rao, C. R. (1973): Linear Statistical Inference and its Applications. (2nd ed.) John Wiley & Sons, New York. [19] Rao, S. S. (1979): Optimization Theory and Applications. Wiley Eastern Limited, New Delhi. [20] Ríos, S., Ríos Insua, S. and Ríos Insua, M. J. (1989): Procesos de decision Multicri-terio. EUDEMA, Madrid, (in Spanish). [21] Steuer, R. E. (1986): Multiple criteria optimization: Theory, computation and applications. John Wiley, New York.