UFDC Home  myUFDC Home  Help 



Full Text  
ACCURACY OF RESPONSE SURFACE APPROXIMATIONS FOR WEIGHT EQUATIONS BASED ON STRUCTURAL OPTIMIZATION By MELIH PAPILA A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2001 Aileme... ACKNOWLEDGEMENTS I would like to express my special thanks and appreciation to Dr. Raphael T. Haftka, chairman of my advisory committee. He has been a continuous and excellent source of guidance for my Ph.D. research. I have always been amazed and motivated by his neverending enthusiasm to explore and willingness to share his knowledge and experience. Working and interacting with him have broadened my vision not only professionally, but also personally. I also would like to thank the members of my advisory committee, Dr. Andre I. Khuri, Dr. Andrew J. Kurdila, Dr. Bhavani V. Sankar, Dr. Peter G. Ifju and Dr. Wei Shyy. I am grateful for their willingness to review my Ph.D. research and to provide me their constructive comments that helped me to complete this dissertation. I must express my gratitude to my former advisor, Dr. Mehmet A. Akgun, who made me realize how much I like research. I have always enjoyed interacting with him in my professional and personal life. I also thank the colleagues with whom I interacted during this study. In particular, I thank Dr. Roberto Vitali, Dr. Satchi Venkataraman, Dr. Raluca Rosca and Dr. Ibrahim Ebcioglu at the University of Florida, and Dr. Hongman Kim at Virginia Tech. The financial support provided by NASA grants NAG12000 and NAG12177 is gratefully acknowledged. I thank my dear friends from home; Murat (MK), Canan, Selami (SYA), Senol, Tarkan, Sartuk, Bulent, Yusuf and Merve who kept constant contact with me. I am grateful for the friendship of Kamal and Alpay who always inspired me with their positive attitudes. I am also grateful to have my best friend Guray and Elif in my life. They deserve special thanks for their constant support and true friendship. The best thing happened to me during my study was getting to know wonderful people, Farshad, Sibel, and their family, and having the chance to make them an important part of my life. I thank my parents, Gonul and Mithat Papila; Meric, Ozer, Cagil and Irmak Torgal; Hurriyet, Ibrahim and Eray Uzgoren; Emine, Bulent and Yigit Topoglu: without the love and support of my family I could not have finished this dissertation. Finally, my deepest appreciation goes to my beautiful and loving wife, Nilay. Her love and presence have been and will be a tremendous inspiration for every challenge in my life. Nilay has given me her support in every respect: emotionally as my better half, professionally as my colleague. I only hope that I have been able to give her the same support in return. TABLE OF CONTENTS page A C K N O W LE D G E M EN T S ......... ................................................................................. iii LIST OF TABLES ................... ............................ .. vii LIST OF FIGU RES ......... .... .. ...... ........ .... .... ................... ..... ix ABSTRACT ............................................. xiii CHAPTERS 1. IN T R O D U C T IO N ........... ................................................. .... ........... ........ F ocus ........................................................ .............. 1 A p p ro ach ........................ ........................................................................ . 7 Scope and O objectives .................. .................... ............... 9 2. LITERA TURE SU RVEY ........................................................................ 12 W ing (System level) Structural W eight Equations ................................. ................. 13 Panel (Component Level) Weight Equations in MultiLevel Procedures ................. 18 A accuracy of R response Surfaces ......................................... ................ .............. 23 Use of RS in Uncertainty Quantification .............. ................ ....... .............. 30 3. RESPONSE SURFACE METHODOLOGY........................ ..................... 32 D design of E xperim ents ..........3......... .. ............ .... .................... ....... ............. 33 Fitting Procedure: Standard Least Squares .......................................................... 36 Adequacy and Accuracy of RS ....................................................... .......... .... 39 Hypothesis Testing in RS M ethodology .............. ............................. ....... ....... 43 Standard L ackofFit T est .................................................................. ... 48 LackofFit with NonReplicated Design Points....................................... 52 Iteratively Reweighted Least Squares.............. ................ ...... .............. 58 4. BIAS AND NOISE ERRORS IN RESPONSE SURFACE APPROXIMATIONS ... 63 M ean Squared E rror C riterion ............................... .. ... .. .............................. 63 Estimated Standard Prediction Error and Residual Based Error Bound ................... 69 Polynomial Examples.......................................... ............ ........ 71 Coefficient of Correlation ............................................................................. 72 Cubic Polynomial in Two Variables without Noise ................ ............... 74 Cubic Polynomial in Two Variables with Noise.............................................. 81 Cubic Polynomial in Five Variables without Noise........................................... 87 Cubic Polynomial in Five Variables with Noise............................................. 89 Effect of Experimental Design: Eigenvalues with a Minimum Bias Design........ 91 Discussion ............................................ ............... 94 5. HIGH SPEED CIVIL TRANSPORT (HSCT) WING WEIGHT EQUATIONS ....... 95 H SC T D design D definition .................................................................... .............. 95 Structural Optimization and Finite Element M odel........................................ ...... 98 Construction of Response Surface Approximations .............................................. 103 Five Configuration Variables ............................... .............. 103 Numerical noise and outlier analysis .................................................... 104 Lackoffit ................................................... 108 Error analysis based on M SEP eigenvalues .......... ................................ 111 Cubic fitting model in the approximation ............. ..... ................. 113 Data vs. RS ....... ........ ......... ................ 117 Repair vs. IRLS .................................. ..................... 118 Ten Configuration V ariables................................................ 121 Design for Uncertainty ....................... .. .. .................... 123 Parameter Uncertainty: Material Properties and Geometry ............................ 125 R S U uncertainty .................. ............................. ................ .. 130 Discussion ............................................ 132 6. CRACKED COMPOSITE PANEL WEIGHT EQUATION................................. 134 Stress Intensity Factor C alculation .................................................. .................... 137 Panel Problem D definition ............... .................. ....................................... 142 Optimization and Response Surface Comparisons .............................................. 145 Low Fidelity R results ..................................... ...... .. ....... .. .......... .. 146 H ighFidelity R results ................... .................. ...... ................ .............. .. 147 Comparison of LowFidelity and HighFidelity Results................................... 150 Discussion ............................................ 153 7. CON CLUD ING REM ARK S .............. ........................................................... 155 APPENDICES A NOISE IN OPTIMIZATION RESULTS DUE TO ROUNDOFF ERRORS.......... 157 B FU Z Z Y SE T TH E O R Y ................................................................... .................... 165 R E FER EN CE S .......... .............. ..... .............. ................ 168 BIOGRAPHICAL SKETCH .................................... ............................ 177 LIST OF TABLES Table page 1. Comparison of actual and predicted weights from Eq. (1.1) ................................ 3 2. Experimental designs used in the present work ................................ .............. 36 3. Verification of the lackoffit test with nonreplicated design points .................. 57 4. Standard weighting functions for IRLS ................................... .............. 60 5. Coefficient of correlation between the absolute residual eTj and the squareroot of maximum eigenvalues m from Eq. (4.27) for two variable cubic polynomial example without noise................................... 77 6. Coefficient of correlations between the absolute residual eTj and the squareroot of maximum eigenvalues x from Eq. (4.27) for two variable cubic polynomial example with noise.................................................. 82 7. Coefficients of correlation between the estimated standard error and absolute true residuals e j for cubic polynomial example with noise............. 85 8. Facecentered central composite design (FCCD) and minimumbias central composite design (MBCCD) for twovariable polynomial example ........ 92 9. Load cases in H SCT design study................................ ........................ .. ....... 96 10. Configuration design variables for HSCT with corresponding value ranges........ 97 11. Statistics for fivevariable full quadratic response surface approximations for HSCT wing bending material weight......... ........... ........................... 105 12. The near replicate design points................ ...... .... ................................. 108 13. Results of the lackoffit test for quadratic RS fitted to 126 design points of fivevariable HSCT wing problem .................................................................. 109 14. Results of the lackoffit test for quadratic RS fitted to 115 design points of fivevariable HSCT wing problem ............... ......... ........................ 110 15. Coefficients of correlation r between Gma and s for fivevariable HSCT wing problem at 3146 design points ............................................. 111 16. Statistics for fivevariable cubic response surface approximations for wing bending material weight ........... ...... ......... .... .............. 114 17. Results of the lackoffit test for cubic RS fitted to data of fivevariable HSCT wing problem ......... ............................................ .............. 115 18. Average Residual Percentage for IRLS fit and RS after Outlier Repair, fivevariable HSCT problem................................................... 120 19. Statistics for 10variable response surface approximations for HSCT wing bending material weight ........... ............. ............... .............. 122 20. RScbc term estimates and statistics, fivevariable HSCT problem ................ 124 21. Uncertainty levels and ranges for fuzzy variables in fivevariable HSCT wing problem...................................... .............. 126 22. Uncertainty for blue and red design: Combined effect variable fuzziness and R S m odeling uncertainty ......... ................................................ .............. 132 23. M material Properties for A S4/35016........... ................................. .............. 144 24. Configurations for threelevel full factorial experimental design.................... 145 25. Summary of RS approximation (weight equation) based on results of low fidelity optim ization ........................................... ......................... ... 147 26. Summary for quadratic RS approximation (weight equation) based on highfidelity results ...... ..... ........................................ ................ .. 148 27. Summary for first order RS approximation (weight equation) based on highfidelity results ...... ..... ........................................ ................ .. 149 28. Plythicknesses of lowfidelity optima for configurations deviating from the general trend ................................................ ........ .. .......... 152 29. Summary of RS approximation (weight equation) based on results of low fidelity optimization after repair of the three outliers by highfidelity ............. 152 30. Convergence history for Configuration8: Convergence setting Case 1............. 160 31. Objective function with different initial structural designs for Config. 8.......... 161 LIST OF FIGURES Figure page 1. Aircraft gross weight variation with wing span .................................. ................ 2 2. Threelevel full factorial (left) and face centered central composite (right) experimental designs for three codedconfiguration variables ............................ 35 3. Sample Fdistributions ................. .............................................. 45 4. The tdistributions of various degrees of freedom df ................ ............... 48 5. 3D representation of clustering design points to generate nearreplicate data for lackoffit test............... ...................... ................... 56 6. 1D example of the effect of an outlier in the data .................. .......... .. 58 7. 1D example of IRLS procedure ....... ........... .. ............. ................. 59 8. W fighting functions for IRLS.................. ......... ........................ .............. 61 9. Linear correlation examples .............. .................................. 73 10. Twovariable threelevel factorial design points............................ ......... ... .. 75 11. Squareroot of maximum eigenvalue max contours from Eq. (4.27) for a twovariable cubic polynomial, fitted with a quadratic model................. 75 12. Absolute residual ej contours twovariable cubic polynomial example (polynomials defined in Table 5), without noise. ............................................ 77 13. Comparison between squareroot of eigenvalues max and absolute residual ed twovariable cubic polynomial example ...................... .............. 78 14. Comparison of squareroot of eigenvalues max and maximum absolute residual e j max contours (maximum of residuals among the four polynomials defined in Table 5) twovariable cubic polynomial example, without noise. ............................................ 79 15. Comparison between squareroot of eigenvalues max and maximum of true residual errors ej max due to the four polynomial in Table 5 two variable cubic polynomial example, without noise................... ................. 79 16. Conservative error bound s contours twovariable cubic polynomial example (polynomials defined in Table 5), without noise.................. ............... 80 17. Comparison between conservative error bound s, and absolute residual ej twovariable cubic polynomial example (polynomials defined in Table 5), without noise..................... ........................ .............. ...... .... ..... 81 18. Absolute residual ej contours twovariable cubic polynomial example (polynomial defined in Table 5 and Table 6), with noise ................................... 82 19. Comparison between squareroot of eigenvalues Gmax and absolute residual ej twovariable cubic polynomial example (polynomials defined in T able 5), w ith noise.................................................. .. ..... .............. 83 20. Comparison of squareroot of eigenvalues max and maximum absolute residual erj, max contours (maximum of residuals among the four polynomials defined in Table 6) twovariable cubic polynomial example, with noise ..... .. ...................................... ............... 84 21. Comparison between squareroot of eigenvalues max and maximum of residual errors erj, max due to the four polynomial in Table 5 two variable polynom ial exam ple, with noise.................................. ..................... 84 22. Normalized estimated standard error ee, / a for a twovariable........................ 86 23. Comparison between conservative error bound s, and absolute residual ej I twovariable cubic polynomial example (polynomials defined in Table 7), w ith noise..................... ...... ...... .... ........ ......................... ...... .. 86 24. Comparison between squareroot of eigenvalues Gmax and absolute true residuals eT and the trend line manually drawn for largest errors characterized by the eigenvalues fivevariable cubic polynomial example, w without noise ........... ..... ...... ......... ..... ....... ..........88 25. Comparison between conservative error bound s, and absolute true residual eg fivevariable cubic polynomial example, without noise............... 89 26. Comparison between squareroot of eigenvalues x = J max and absolute residuals y = e fivevariable cubic polynomial example, with noise.............. 90 27. Comparison between conservative error bound s, and absolute residual e1 fivevariable cubic polynomial example, with noise ..................... 91 28. Eigenvaluecontours ( kmax ). a) Facecentered central composite design (FCCD); b) minimumbias central composite design (MBCCD) for two variable polynom ial exam ple ................................................................. ....... 93 29. Absolute residualcontours in twovariable polynomial example (for Polynomial 1 in Table 5). a) Facecentered central composite design (FCCD); b) Minimumbias central composite design (MBCCD)....................... 93 30. HSCT wing planform and five & 10configuration variables. .......................... 97 31. Flow chart for the RS based wing weight equation ....................................... 99 32. HSCT finite element model and structural design variables............................. 99 33. Closeup view for the typical wing box cell ........................................... .... 100 34. Flow chart of the optimization in GENESIS.......................... 101 35. Effect of overly tight optimization parameters: Difference in optimization results between Case 1 and Case 3 ............................................. 102 36. Outliers and effect of repair on quadratic response surface, fivevariable H SCT problem .............. ...... .............................. ............ .. 107 37. Comparison between squareroot of eigenvalues max and conservative error bound s, quadratic RS of HSCT 126 design points.......... 112 38. Comparison between squareroot of eigenvalues max and conservative error bound s, quadratic RS of HSCT 115 design points.......... 112 39. Comparison of predictions of squareroot of MSEP, s, and s, at 115 design points quadratic RS................................... .............. 113 40. Comparison between squareroot of eigenvalues max and absolute error between reduced cubic and quadratic RS fivevariable HSCT wing problem........................................... ............ 117 41. Noisy wing bending material weight Wb compared to smoother objective function (total wing structural weight)........................ ...... .............. 118 42. Extreme value design line. ................................... .............. 119 43. Comparison of RS with outlier repair and with outlier elimination over the design line, fivevariable HSCT problem. ............. ........................... ........ 120 44. Description of vertex method for threevariable case............................ ... 127 45. Fivevariable HSCT case: Contours in plane with almost constant wing bending. a) Nominal wing bending material weight contours; b) parameter uncertainty band width (upper bound minus lower bound) contours; c) Modeling (RS) uncertainty band width contours........ .................................. 129 46. Comparison of predictions of squareroot ofMSEP, s, (=x) and s, (=y) at 152 design points of fivevariable HSCT wing problem .............................. 131 47. Blue and Red design planform s .............. ...... ......................................... 131 48. Fit over stress values to find K ......... .. .. ... .. ........................ .... .......... 140 49. Stress distribution in the vicinity of the crack tip............................................. 141 50. Composite BladeStiffened Panel with a Crack ....... .... ................................ 143 51. Comparison of prediction by RS approximation and optimization data based on lowfidelity optimization for stiffened panel with crack ................... 146 52. Comparison of prediction by RS approximation and optimization data based on highfidelity optimization for stiffened panel with crack ...... .......... 149 53. Ratio of highfidelity and lowfidelity optimal weights............................ 150 54. Comparison of relative errors wrt highfidelity results obtained by low fidelity RS before and after repair.............................. ................. 153 55. Difference in results on PC and UNIX: Convergence setting Case 1............... 159 56. Difference in results on PC and UNIX: Convergence setting Case 2................. 159 57. M em bership functions ............................ ............ ...................... .............. 166 Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy ACCURACY OF RESPONSE SURFACE APPROXIMATIONS FOR WEIGHT EQUATIONS BASED ON STRUCTURAL OPTIMIZATION By Melih Papila December 2001 Chair: Raphael T. Haftka, Distinguished Professor Major Department: Aerospace Engineering, Mechanics and Engineering Science Accurate weight prediction methods are vitally important for aircraft design optimization. Therefore, designers seek weight prediction techniques with low computational cost and high accuracy, and usually require a compromise between the two. The compromise can be achieved by combining stress analysis and response surface (RS) methodology. While stress analysis provides accurate weight information, RS techniques help to transmit effectively this information to the optimization procedure. The focus of this dissertation is structural weight equations in the form of RS approximations and their accuracy when fitted to results of structural optimizations that are based on finite element analyses. Use of RS methodology filters out the numerical noise in structural optimization results and provides a smooth weight function that can easily be used in gradientbased configuration optimization. In engineering applications RS approximations of low order polynomials are widely used, but the weight may not be modeled well by loworder polynomials, leading to bias errors. In addition, some structural optimization results may have highamplitude errors (outliers) that may severely affect the accuracy of the weight equation. Statistical techniques associated with RS methodology are sought in order to deal with these two difficulties: (1) highamplitude numerical noise (outliers) and (2) approximation model inadequacy. The investigation starts with reducing approximation error by identifying and repairing outliers. A potential reason for outliers in optimization results is premature convergence, and outliers of such nature may be corrected by employing different convergence settings. It is demonstrated that outlier repair can lead to accuracy improvements over the more standard approach of removing outliers. The adequacy of approximation is then studied by a modified lackoffit approach, and RS errors due to the approximation model are reduced by using higher order polynomials. In addition, remaining error in the RS weight equation is characterized, and two error measures are introduced. One of the measures is qualitative, and it is used to identify regions of design space where RS accuracy may be poor. The second measure is quantitative, but conservative. Its use is demonstrated for incorporating the uncertainty due to RS weight equation with parameter uncertainty in order to compare competing designs for robustness. CHAPTER 1 INTRODUCTION An aircraft designer would probably place weight at the top of the list of important concerns in the design of an aircraft. This is because all disciplines in the design process such as aerodynamics, structure, performance and propulsion interact with each other through the weight. The tradeoff between aerodynamic and structural efficiency is a common example of such a discipline interaction. Structural weight, a big proportion of the total weight, affects the lift that must be generated, and consequently the drag. In turn aerodynamic design aspects affect the loads, and hence structural design, and the structural weight. Focus The focus of this dissertation is accurate and affordable weight estimation. The impact of the weight starts from the very beginning in aircraft design process. Definition of conceptual characteristics initiates the process, and it is followed by the preliminary design phase where prediction of the weight is one of the first steps. The final size and cost of the aircraft are strongly influenced by the estimated weight at the preliminary design phase. The more reliable and precise starting data are, the faster is the convergence to the promising configurations. Empirical weight equations. It is a common practice in conceptual design to estimate the weight by empirical equations derived from historical aircraft data. They are one of the earliest approaches used for weight prediction. Early wing weight prediction methods, for instance, were of a purely statistical nature and included only simple relationships between wing structure weight and the most significant design parameters such as wing span and wing gross area. If only a nominal wing span is available to a designer, a purely statistical approach is simply to collect and plot weight versus wing span data for comparable existing aircraft, and then try to make a fit to obtain an expression for weight as a function of wing span. For example Figure 1 presents statistical/historical data (source: Aviation Week & Space Technology, March 16, 1992, pp. 71119) for the gross weight versus wing span on a logarithmic scale for existing aircraft. A first order polynomial curve fit on this plot gives an empirical equation for aircraft gross weight. 7.00 Data 6.00 Linear fit 5.00 y = 2.0186x+ 1.1524 S4.00 0 3.00 o 2.00 1.00 0.00 0.00 0.50 1.00 1.50 2.00 2.50 x=log(Span, ft) Figure 1: Aircraft gross weight variation with wing span determined by statistical evaluation of existing aircraft Empirical equations are still widely used for first order estimates or as reasonability checks for more detailed analysis. Most airplane manufacturers have developed their own methods based on such empirical relations (Niu, 1988; Loretti 1997). An example of how to develop an empirical equation can be found in Staton (1996). Another example of empirical equations presented in Torenbeek (1992) is a wing weight equation, Ww for subsonic transports and executive jets: Ww =17bS WZ (1.1) where b is wing span, S is wing gross area, WMZF is maximum zero fuel weight, and WMTO is maximum takeoff weight. Table 1 compares the actual and estimated weight for subsonic transport and executive jet airplanes, from Eq. (1.1). Table 1: Comparison of actual and predicted weights from Eq. (1.1) from Torenbeek (1992) Aircraft WMTo, (kN) Ww, estimated (kN) Ww, actual (kN) Error, % Airbus A300/B2 1343.6 188.43 196.31 4.0 Airbus A340 2487.0 295.32 340.85 13.4 BAC 111/300 387.0 38.60 42.90 10.0 Cessna Citation II 59.16 6.79 5.73 +18 Boeing 727100 71.8 75.89 79.02 4.0 Boeing 737200 513.8 39.79 47.21 15.7 Boeing 747100 3158.4 446.15 383.35 +16.1 Fokker F28 Mk 4000 315.8 30.22 33.28 9.2 Lockheed 10111 1912.8 224.59 210.86 +6.5 McD. Douglas DC9/30 480.4 40.60 50.71 19.9 McD. Douglas MD80 622.8 60.95 69.22 11.9 McD. Douglas DC10/10 2024.0 250.46 217.97 +14.9 McD. Douglas DC10/30 2468.9 252.74 261.83 5.3 The main advantage of empirical weight equations is that they are easy to use. They do not require expertise in weight engineering, so aerodynamics or performance experts can easily use them. They do not require sophisticated computational resources, but they also have some limitations that should be considered. Empirical equations are generally a function of exterior geometric parameters, such as wingspan and thickness to chord ratio, and they may not reflect interior structural details such as frame, stiffener shapes and spacing. Since they are based on existing aircraft data, empirical equations may be misleading for unconventional designs or for new concepts that require extrapolation instead of interpolation from the data. This may happen, for instance, when new materials or technology appear (advanced composite materials, bonded construction), or when operating temperatures become much higher than those for previous structures. When conceptual characteristics are not far beyond the available knowledge and experience (in other words when interpolation is possible) accuracy by the empirical equations may be of the order of 25% (Niu, 1988). However, as Schmidt et al. (1997) suggested an error of 2% in fuselage weight prediction may be equivalent to about 1.5 tons of freight or 20 passengers. Alternative to empirical weight equations. Historically, structural weight estimation of aircraft that lie outside the range of applicability of empirical methods has been determined by stressanalysis based approaches or by construction of prototypes. The latter approach is very expensive and not practical unless the design is converged, and verification is needed rather than the prediction. Stressanalysis based structural weight prediction procedures are widely used throughout the design phase. Howe (1996) presented a comparison between an empirical and a stressanalysis based wing weight equation. Data for 118 aircraft ranging from the Cessna 150 to the Boeing 747/200 were used for comparison. Howe reported that both equations gave reasonably good results, but the latter were more accurate. However, it was also noted that the test data were limited to aircraft of light alloy construction. Predictions by empirical equations can be misleading for aircraft of different and advanced materials such as composite materials. That is because the use of such materials has been dealt with somewhat arbitrary factoring in empirical equations as historical data includes primarily the aircraft of traditional materials. Weight determination by stress analysis can be done at different levels of complexity ranging from closed form beam and plate solutions early in the design process to advanced finite element (FE) analysis at later phases of the process. Designers seek the lightest practical arrangement of material and structural details capable of withstanding the required loads. Depending on the complexity of the analysis and weight prediction, the computational time required for a lightweight design may be substantial. With advances in computational power, the level of complexity that can be afforded in the preliminary design phase is increasing. For instance, FE analysis usage in aerospace weight engineering has progressed to the point that preliminary structural designs may be analyzed with low levels of uncertainty when combined with a corresponding knowledge of the asbuilt structural design. However, the multi disciplinary nature of the modern design process cannot afford long cycle times, since it entails frequent changes and extensive tradein studies. In addition, conceptual designers are not usually experts in FE modeling. That is why these designers still demand computationally affordable, easy to use, but also reliable approaches. As a consequence, the idea of combining the reliability of stressanalysis based prediction and the ease of manipulation of analytical expressions in multidisciplinary studies has attracted interest. In general, structural weight equations have been developed for various types of wing and fuselage structures since they are primary loadcarrying elements. Wing and fuselage structures are composed of a large number of curved or flat, usually stiffened panels. Design optimization of all panels simultaneously, in other words the complete structure, requires detailed structural modeling and appears to be beyond the present computational capabilities since the panels have complex geometry and failure modes. The current practice is to design wing or fuselage structures at two levels, global and panel, and this can increase the required time for weight estimation. In addition, other problems that affect stressanalysis based weight prediction and design optimization are the following Numerical simulations based on stress analysis such as structural optimizations are prone to numerical noise due to reasons including discretization errors, incomplete convergence of iterative procedures, and roundoff errors. When these simulations are used in a multidisciplinary optimization procedure employing gradientbased optimization algorithms the procedure suffers from the noisy calculations, and may fail to find the optimum. As a consequence, smooth responses are desired in optimization problems. To reduce the errors and the noise generated by the computational models the complexity of the model may be increased. A finer finite element mesh, for example, can be used in structural analyses. This, however, increases the computational cost and design cycle time. *Oftentimes interactive use of different structural codes is required for capturing various failure modes, such as buckling, delamination, or crack propagation. Integration of the different codes is not an easy task. Approach The approach taken into consideration is weight equations by response surfaces. Response surface (RS) techniques (Myers and Montgomery, 1995: Khuri and Cornell, 1996), which fit a large number of analyses with simple functions, offer an attractive set of features that address the problems mentioned in previous section: They filter out the noise and help to identify outliers They may require large computational effort beforehand to construct, but this pays off during the design process since they are easy to compute They can act as interfaces between the different levels of analysis or between different disciplines. Response surface methodology can be used to construct weight equations. Approximations by RS are similar to empirical relations in the sense that they both employ statistical fitting procedures. The source data used, however, are different. RS based equations utilize the data created by stressanalysis based structural optimizations at preselected design points, providing an overall perspective within the design space (design of experiments) whereas empirical relations employ historical aircraft data. Applications for RS based weight equations. Multidisciplinary optimization (MDO) has become increasingly important with the increasing complexity of engineering systems. SobieszczanskiSobieski and Haftka (1997) presented a survey for the recent developments in this field. They indicated the two main challenges of MDO are organizational complexity and computational expense. Response surface techniques are particularly suitable for MDO applications. They can be used for approximating the objective functions and constraints. For the organizational challenge these techniques provide a convenient representation of the data from one discipline to other disciplines and to the system. For the challenge of computational expense RS techniques may first seem inappropriate since they require a large number of analyses for data generation beforehand. However, the simplicity and insignificant cost of calculations by RS once constructed make them attractive in multidisciplinary design optimization studies. The organizational challenge exists also in a singlediscipline application. Oftentimes design optimization uses procedures employing different levels of analysis and optimization as in global/local structural design. In these applications, integration of different codes is generally required. Approximations by RS techniques may provide an easy means to connect such codes. They may also be used to correct the lowcost analysis results with the help of a small number of highfidelity results whose computational cost is often prohibitive to use at each design throughout the design process. Another important application of RS techniques is design against uncertainty. The random nature of this problem makes it computationally intensive, and it typically entails at least an order of magnitude higher computational cost than deterministic design (Venter and Haftka, 1998). Response surface approximations can be used to reduce this high computational cost since they are easy to calculate and implement. The designer, however, should keep in mind that the RS approximation itself may impose another type of uncertainty due to modeling error that needs to be addressed in uncertainty studies. Challenges associated with RS based weight equations. Structural optimizations that are often used to generate data for RS based weight equations are prone to numerical noise. Optimization based on finite element analyses, for example, may introduce noise due to discretization errors affecting the analysis results or even due to roundoff error characteristics of the computation platform affecting progress by the optimizer (Papila and Haftka, 1999). Although one of the most appealing features of RS approximations is the ability to filter the data, highamplitudenoise outliers that may be experienced in structural optimizations may cause the loss of accuracy of the approximation or weight equation. Therefore, it is important to make sure the data for constructing the weight equation are free from severe flaws before relying on the weight equation. In case of data generated by structural optimizations, for instance, optimum weight data free from flaws such as incomplete convergence would increase the reliability of the weight equation. Even when the data generated are free from severe flaws or outliers, RS based weight equations are still approximations. This is because the number of the data is usually limited and it does not cover every possible design point. In addition, traditional polynomial fitting models may not represent equally the true nature or complexity of the weight function for every region in the design space. Therefore, it is also important to study the accuracy of the weight equation and to pinpoint the design regions where the approximation may be poor so that designer may steer away while using the weight equation or at least be cautious. Scope and Objectives The main thrust of this dissertation is to address the challenges associated with the RS based weight equations as mentioned in the previous section, that is, how to deal with the errors that are introduced by the process of fitting structural optimizations with RS for structural weight equations. This includes reducing error by identifying and repairing outliers, studying the accuracy of approximation, reducing RS errors due to the approximation model by using higher order polynomials, and characterizing the remaining error for use in design for uncertainty. Therefore, the two main objectives of the work are the following: Reduce errors associated with RS constructed on the basis of structural optimizations: To reduce noise error the aim is to identify the highly erroneous results and label them as outliers within the data. For the data generated by structural optimizations outliers may be due to reasons such as failure of convergence and local optima. Once detected they may be studied further for possible correction. The corrected or repaired data set is expected to increase the accuracy and reliability of the RS. To deal with modeling (bias) errors, higherorder polynomials are used. Characterize the remaining error and investigate statistical tools to determine the regions of design space where RS accuracy may be poor. A literature survey of wing and panel weight equations, use of RS techniques for accurate weight estimations and for analysis/optimization code integration, and the use of RS based weight equations in design sensitivity to uncertainty are presented in Chapter 2. Chapter 3 provides the background of RS methodology. Uncertainty due to use of RS is given in Chapter 4 by deriving a measure for high modeling error locations in the design space. Chapter 5 presents an application of the methodology for weight equations of a high speed civil transport wing weight problem, and its use in design for uncertainty. 11 Chapter 6 summarizes similar weight equation construction for a stiffened composite panel with a crack problem and Chapter 7 presents the concluding remarks for this dissertation. CHAPTER 2 LITERATURE SURVEY Weight equations allow designers to design systems without getting bogged down in the redesign of components as they change system characteristics. In aircraft design the upper level system is the aircraft itself and the main structural components are wing and fuselage. These in turn can be considered as systems of components such as panels, frames, etc. The flight optimization system (FLOPS) wing weight equation (McCullers, 1984), for example, predicts how the aircraft structural weight changes with use of composite materials without making panel and wing level detailed and computationally burdensome structural analyses or optimization. Therefore, accurate and easytouse weight equations are a useful and important tool for the designer. Shanley (1960) published weight prediction methods based on elementary strength/stiffness considerations, augmented by experimental results and statistical data. This approach enabled weight engineers to obtain more accurate weight estimates than traditional empirical weight equations and to consider sensitivity of weight to the design. Moreover, it also paved the way to the 'station analysis' methods and use of more involved analysis models. In general, structural weight equations have been developed for various types of wing and fuselage structures since they are primary loadcarrying elements. Shanley (1960) reported that existing empirical formulas were sufficiently efficient for the other portions of the airplane structure. This appears to hold even for recent aircraft concepts. Balabanov et al. (1996, 1998 and 1999) investigated weight due to wing bending resistant structural elements, wing bending material weight, of a futuristic highspeed civil transport (HSCT). They used empirical formulas from the FLOPS program (McCullers, 1984; 1997) for the entire structure except for the wing structural weight as in Grasmeyer et al. (1998) and Gern et al. (1999) for another unconventional concept, transport aircraft with strutbraced wing. The focus of this work is on accurate weight equations for wings and for panels. The former are used by the aircraft designer to design the entire airplane; the latter are used by the wing or fuselage designer to design the entire wing or the entire fuselage. Wing (System level) Structural Weight Equations Wing weight equations started with fitting of historical data using specially tailored expressions with variables raised to various powers. They progressed to equations based on stress analysis of simple beam models of the wing and fuselage that were adjusted to fit historical data. Shanley (1960) reviewed the work by Lipp and Kelley. Lipp (1938) developed formulas for the structural weight of wings by actual integration of assumed lift distribution in the basic equations for the conditions of shear and bending. The final equations were expressed in terms of load factor, span, root thickness, allowable stress and taper ratios for chord and thickness. Kelley (1944) used more accurate spanwise lift distributions that increased the difficulties in the integration and presented results in chart form for engineering use. Shanley (1960) modeled the wing as a box beam whose flanges carried the bending load whereas webs carried shear loads. He indicated that the most important terms are span, mean aerodynamic chord, and the ratio of wing span to depth. Torenbeek (1992) presented a method providing estimations for the primary parts of a wing structure. The method makes use of elementary stress analysis combined with historical data. The wing group weight is expressed as the sum of the primary weight (top and bottom covers, spars, ribs, and attachments) and the secondary weight comprised of the weight of components in front of the front spar, components behind the rear spar, plus any miscellaneous weight. The basic weight is the minimum weight of the primary structure required to resist bending and shear loads where stresses are taken into account. Formulas based on historical data were used in calculations of the secondary weight and for the correction to the primary weight due to factors such as joints and mountings for engines. Prediction errors of less than 4% were found for four representative transport aircraft with takeoff weights between 50 and 3500 kN. Weight equations also progressed to using a multitude of structural design optimizations (a review for structural optimization was presented by Vanderplaats, 1999) based on simple plate models. McCullers (1984) developed the weight equations in FLOPS by combinations of historical data and structural optimizations based on equivalent plate models for wings using the wing aeroelastic synthesis procedure (WASP) and the aeroelastic tailoring and structural optimization (ATSO) program (McCullers and Lynch, 1972; 1974). Optimum wing designs data were used to account for aeroelastic tailoring via composite materials, forward sweep angle effects and flutter and divergence weight penalties for very high aspect ratio wings on the wing weight. The equations are based on analytical expressions of primary wing configuration parameters such as aspect ratio, thicknesstochord ratio, sweep angle, etc. correlated by constant factors based on historical data and results of structural optimizations. For different categories of aircraft the constants and the expressions vary, so that equations are applicable for a range of aircraft types, from military fighters to civil transonic and even supersonic transport aircraft. These and similar weight equations attempt to fit any airplane, or at least a broad class of airplanes. These types of equations are referred to as quasianalytical estimation in Staton (1996). They offer very good estimation accuracy for existing technology aircraft, and through their sensitivity to engineering parameters allow accurate trade studies. Their advantages over pure statistical equations lie in their use in extrapolating outside parameter ranges in the database of aircraft. One modem trend has been to replace wing weight equations with structural optimization based on beam models inside the aircraft optimization problem. For example, Kroo et al. (1991) used this approach in aerodynamicstructural design studies of joinedwing aircraft. They combined a vortexlattice aerodynamic model with fully stressed beam sizing routine. Two structural box models: a symmetrical box of uniform cap and web thickness, and an asymmetrical arrangement of cap thickness so that diagonal corners of the box have the same thickness were used. They investigated selected joinedwing designs and studied the effects of design parameters on drag and structural weight. They employed a simple structural optimization following the aerodynamic analysis for each design. In later studies (Kroo et al., 1994; Braun et al., 1996), they used twolevel collaborative optimization that allows the designer to include other disciplines besides structures and aerodynamics. Collaborative optimization allows some autonomy for each discipline to optimize on its own, but the twolevel formulation can be very ill conditioned (Alexandrov and Lewis, 2000). Gern et al. (2000) also embedded a beammodel structural optimization for weight calculations of a strutbraced transport wing inside the aerodynamic optimization. They employed a fully stressed design procedure accounting for aeroelastic load distribution to estimate the wing bending material weight. The remaining portion of the wing structural weight due to flaps, slats, spoilers, ribs etc. was estimated by FLOPS weight equations. Finite element (FE) analysis based procedures have also been used in the preliminary design phase. Advances in the computational resources allow the use of at least coarse FE models in this phase. Huang et al. (1996) performed simultaneous structuralaerodynamic optimization by combining a FE based structural optimization for wing bending material weight (structural weight needed to resist bending) with a simple weight equation. Weight prediction from structural optimization was first compared with two quasiempirical wing weight equations for ten high speed civil transport (HSCT) wing planforms. The equations were from FLOPS (McCullers, 1984; 1997) and by Grumman Aerospace Corporation (York and Labell, 1980). The comparison indicated that the weight equations predicted well the average structural weight and effects of airfoil thickness on structural weight, but did not model accurately all of the effects of planform changes on the structural weight. The FLOPS weight equation was chosen since it agreed better with the structural optimization results than Grumman equation. The structural optimization and weight equations were combined through the following interlacing procedure: structural optimization and FLOPS weight estimations were found for the initial design and their ratio was calculated. Five cycles of aerodynamic design optimization were performed by using the FLOPS weight estimation corrected by the ratio. After five cycles, the structural optimization was repeated for the new aerodynamic design and the ratio was updated. These steps were repeated until convergence. The interlacing procedure of the FLOPS weight equation and the FE based structural optimization results was successful in terms of obtaining designs reflecting realistic weight estimation balanced with computational cost. It did not protect, however, the optimization from getting stuck at a design where the simple equation was too optimistic. It was noted that derivative of the ratio or scale factor may be helpful for this problem. Performing repeated structural optimizations during the aerodynamic design process has the advantage that use of a general structural weight equation that must fit all transports is not required. However, this approach has three main disadvantages: noise, illconditioning, and code/discipline integration. Giunta et al. (1994) showed that numerical noise caused convergence problems in design optimization for HSCT employing gradientbased optimization techniques. Balabanov et al. (1996) had to deal with numerical noise of several sources in their HSCT design problem. Venter and Haftka (1998) experienced noise associated with the dependence of discretization error on their shape design variables. Multilevel treatments, such as collaborative optimization can lead to illconditioning and computational difficulties (Alexandrov and Lewis, 2000). In addition, the integration of structural optimization software with the other disciplines also requires substantial effort (e.g. Rohl et al., 1994). An alternate approach that avoids the above difficulties is to generate a weight equation for a particular aircraft by fitting the data generated via finite element based optimizations within the context of RS methodology. This approach has the advantages of the traditional weight equation, but without its disadvantages. The main difficulty is obtaining adequate accuracy, and this is one of the thrusts of the present work. Response surface techniques have become an important tool in multidisciplinary optimization (MDO) studies where numerical noise is often an issue besides the organizational difficulty of coupling simulations from different disciplines. They filter out numerical noise, and provide a convenient representation of data from one discipline to another and are easy to interface with an optimizer due to their ease of implementation. Kaufman et al. (1996) fitted quadratic polynomial RS approximation to the structural weight obtained with multiple structural optimizations with GENESIS (VMA, 1997). The HSCT configuration optimizations performed with the RS produced superior results than with the FLOPS weight equation. Balabanov et al. (1996, 1998 and 1999) investigated RS construction for wing bending material weight of HSCT based on structural optimization results of number of different configurations and its use in the configuration optimization problem. They used a larger number of lowfidelity (LF) models and tuned the results based on small number of highfidelity (HF) analyses (Balabanov et al., 1998) in order to handle the computational cost problem associated with FE based structural optimization. Papila and Haftka (1999, 2000a) also constructed RS weight equations based on structural optimizations (see chapter 5). Panel (Component Level) Weight Equations in MultiLevel Procedures Wing and fuselage structures are composed of large number of curved or flat, usually stiffened panels. Designing all the wing panels simultaneously constitutes a complex optimization problem requiring detailed structural modeling of the entire wing and appears to be beyond present computational capabilities. For integrated structural wing (as system) and panel (as component) design, there are similar approaches to those described in the previous section. In the past, there were simple, stressanalysis based panel weight equations like Gerard's (1960). He performed minimum weight analysis of various forms of stiffenedplate construction considered as the compression covers of multicell wing box beams. This, however, did not take into account the effects of panel stiffness on load redistribution (designer can help one panel by beefing up another) and overall stiffness constraints. The current practice is to design the system, wing or fuselage structures at two levels; global (system) and panel (component) levels. At the global level, the structure is modeled with moderate level of detail for individual panels, and the internal load distribution is obtained. With these loads, panellevel design optimization for detailed geometry is then performed by programs such as Panel Analysis and Sizing Code (PASCO) by Stroud and Anderson (1981) or PANDA2 by Bushnell (1987). Depending of the complexity of the problem, integration of global and panel levels may become difficult. Schmit and Mehrinfar (1981) showed that the integration of the levels can be handled well for simple structural models by a twolevel decomposition similar to collaborative optimization (CO) approach. One difficulty associated with the decomposition; elimination of the component level (local) variables in terms of system level (global) variables was mentioned in Haftka and Giirdal (1992). For panel problems, as well as for complex truss and frame crosssectional forms, it is impossible to find analytical expressions for eliminating local variables and replacing them n il/ global variables. It is possible to keep both local and global; variables, and supplement the problem ii i/h equality constraints that guarantee the consistency of the global variables i/ i/h the local variables (Haftka and Giirdal, 1992). However, this approach often tends to make the optimization problem more illconditioned (e.g. Thareja and Haftka, 1986). As noted in the previous section twolevel optimization suffers from ill conditioning and the fact that panel (local) level optima are not usually smooth function of the global level design variables. A successful integration of wing (system or global) and panel (component or local) level structural optimizations was performed by Rohl et al. (1994) in threelevel decomposition approach for the preliminary aerodynamicstructural design of an HSCT wing. The levels of the approach from top to bottom, made use of FLOPS as a general aircraft sizing and performance code for aerodynamic optimization, ASTROS as a structural optimization tool based on FE analyses for wing box (wing level), and PASCO for panellevel optimization. The main purpose of the paper was to demonstrate the integration of PASCO into ASTROS to supply buckling constraint evaluation for wing level optimization, and its application in a multilevel wing design procedure. The overall integration and execution of the procedure were performed by UNIX shell scripts calling specifically designed finite element preprocessor for ASTROS using the FLOPS output, PASCO preprocessor, and postprocessor. The results of this integration prototype without RS application indicated that it worked, but required the designer in the loop making decisions in buckling optimization. Besides the modifications on the codes and customized interfaces, it was noted that fully automatic run requires a logic implementation. Response surface approximations offer an attractive means also for integrating the design of a single panel with the overall structure and for integrating different codes. Ragon et al. (1997) introduced a methodology for global/local design optimization of large wing structures. They used ADOP, Aeroelastic Design Optimization Program of McDonnel Douglas in global/local optimization, and PASCO, Panel Analysis and Sizing Code as the local design code. The methodology starts with generating many optimal designs for each panel type for a range of loading and inplane stiffnesses (All and A66). The data generated were then utilized to construct an RS weight equation as a function of the inplane loads and stiffness parameters for the optimal panel weight satisfying all the local failure constraints. However, for some values of the stiffness parameters no feasible designs exist. When the panel RS weight equation was interfaced with ADOP, during the global iterations ADOP tended to move into infeasible portions of the design space. This was because no data corresponding to the infeasible design points was included in the RS model, and the design regions bordering the infeasible regions tended to be those with low structural weights. This problem was corrected by constructing a new RS using the data that included penalty terms for those original design points for which PASCO could not generate a feasible design. The penalty terms required higher order (cubic) polynomial for accurate fitting. The use of RS in terms of stiffness constraints is an advance over traditional weight equations. It allows the global design to control the load carried by the panel by specifying the panel inplane thickness. In other global/local procedures responses other than weight were approximated by RS easing the code integration. Venkataraman and Haftka (1997), for example, utilized separate local and global analysis codes for a liquidhydrogen tank structure design made of laminatedcomposite material. They used NASTRAN for the global model of irregular geometry and PANDA2 for the local regular stiffened panel analysis. An iterative design procedure between NASTRAN and PANDA2 failed to converge since the latter did not have any information of the global constraints. They overcame this problem by a RS approximation for the global response used in PANDA2. Liu et al. (2000) aimed to solve the infeasibility problem faced by Ragon et al. (1997) who needed to satisfy both stiffness and strength constraints at the panel level. Instead of weight equation Liu et al. developed failure load equation. Their procedure was applied to laminated composite wing panels when the design involves discrete and combinatorial optimization. Global and local design processes were coordinated through an equation that predicts the buckling load multiplier. The equation was a cubic RS as a function of number of plies of each orientation and the loads on the panel. Response surface was fitted to a large number of panel stacking optimization based on analytical solution for simply supported unstiffened panel in the configuration space of number of plies and the loads. The reason for the cubic model is that buckling load is proportional to the cube of the thickness. They used Doptimal experimental design points where panel level optimizations by permutation genetic algorithm (GA) for the stacking sequence. The objective function was maximum buckling load. Then winglevel optimization was performed by Genesis using the cubic RS in the buckling constraint. They applied the methodology for 6, 18, and 54variable cases. It was shown that a cubic RS can fit to the buckling load of the optimal panel stacking sequence as a function of loading on the panel and the given number of plies, and it can be used in order to avoid the problem of infeasible designs faced by Ragon et al. (1997). Ragon et al. (1997) and Liu et al. (2000) generated RS based equations with panel (component) design based on simple models. One goal of the present work is to push the procedure further in the context of general FE modeling of stiffened panels. Accuracy of Response Surfaces Increasing use of RS techniques in design optimization studies has brought attention to ways of increasing the accuracy of the RS approximations. Two sources of error have been mainly investigated: numerical noise and modeling error due to approximation functions. Efforts for identifying the cause of the numerical noise, for improving the data by reducing the noise and by proper selection of design points Design of Experiments (DOE) and for reducing modeling errors pay back in the accuracy of the RS and in its reliable application. Numerical experiments may be noisy due to several reasons including discretization errors, incomplete convergence of iterative procedures, and round off errors. In order to reduce the effect of noise in RS accuracy DOE offers variance optimal designs. Standard designs such as central composite design (CCD), for example, are quite efficient against noise effect (Montgomery and Myers, 1995). However, there are highdimensional problems where more economical designs are generally required. D optimal design criterion (Montgomery and Myers, 1995; Khuri and Cornell, 1996), for example, is commonly used in such problems. Balabanov et al. (1996) used Doptimal design for 29variable HSCT design problem in order to reduce the number of structural optimizations. Response surface of an effective experimental design providing low variance level can be used as a datafiltering tool. Venter et al. (1998a and 1998b) made use ofRS approximations for filtering noise in finite element analyses. In their case the noise was associated with the dependence of discretization error on shape design variables. Giunta et al. (1994) found numerical noise in their calculation of aerodynamic drag components for the aircraft that caused convergence problems in design optimization for HSCT. They demonstrated the noise in a two dimensional example that showed variations were small so that at all points the accuracy of the drag was acceptable, but oscillatory behavior created difficulties for the gradientbased optimization technique. Quadratic RS approximation was used to filter out the noise in aerodynamic simulations. Giunta et al. (1997) applied the RS approach for 10variable HSCT configuration optimization problem. They reported improvement in designs obtained based on the smooth RS approximation compared to the original noisy simulations. Use of RS resulted in virtually identical optimal design starting from two different initial designs whereas using the original noiseproducing aerodynamic model method resulted in different optimal designs. Balabanov et al. (1996) investigated the noisy behavior of wing bending material weight of HSCT designs (that was indicated by Kaufman et al., 1996, but without identifying its source). They found that a large portion of the noise was due to incomplete optimization of the wing camber that affects the wing bending material weight via the aerodynamic loads. Much of the remaining numerical noise was due to the structural optimization process itself through incomplete convergence. Accuracy of bending material weight RS after correcting the camber optimization problem increased substantially compared to one in Kaufman et al. (1996). To improve the accuracy further, minimum bias and minimum variance design strategies were tried. They worked comparably well in terms of rootmeansquareerror (RMSE) and average error. The maximum error, however, is about 30% less in minimum bias design based RS. While lowamplitude, random noise is filtered well by RS approximations, data with large errorsoutliersmay cause significant loss of accuracy. Therefore, robust regression employs techniques for statistical methods, such as iteratively reweighted least square (IRLS) fitting (Holland and Welsch, 1977), to detect and remove or weight down outliers. Once detected, outliers can be investigated further for possible mistakes. These techniques help also to identify flaws in analysis/optimization procedures so that improvement on the results can be achieved. Optimization procedures often produce poor results due to algorithmic difficulties, software problems, or local optima. When a single optimization is flawed, it may be difficult to detect the problem, because the ill conditioning responsible for the problem may make it difficult to apply optimality criteria unambiguously. However, in many applications a large number of optimizations are performed for a range of problem parameters (e.g. Balabanov et al., 1999). When such multiple optimizations are available, statistical analysis can be done to detect incorrect optimal results. For example, Papila and Haftka (1999 and 2000a) and Kim et al. (2000 and 2001) used IRLS procedures to detect such outlier points. More details on these references are presented in Chapter 5. Researchers have also found ways for dealing with modeling error. For instance, it can be decreased by proper selection/definition of variables/factors. Concept of intervening design variables and intervening functions was introduced by Schmit and his coworkers for accurate approximations (e.g. Schmit and Farshi, 1974). They are transformations of the original variables and functions or arguments of the functions. Finding/defining intervening variables as functions of the original variables would help to model response with lower order polynomials (or lower order derivatives in case of local Taylor series approximations), and to decrease computational cost to obtain the approximation (Haftka and Girdal, 1992). Another way of decreasing modeling error and increase the RS accuracy is use of reasonable design space approach (Balabanov et al., 1997). The approach seeks inexpensive constraints such as simple geometric constraints that prevent combinations of design variables resulting in unreasonable geometry configurations. Eliminating large portions of the design space occupying those unreasonable configurations renders the design space more similar to a simplex. The approach improves the modeling accuracy of the RS by shrinking the region where the RS is fitted. Kaufman et al. (1996) investigated the ways of improving the RS accuracy for weight equations. They looked for intervening variables that may reduce the number of data points or increase the accuracy. They also used inexpensive approximate analysis methods and geometric constraints to find the regions of the design space where reasonable HSCT configurations are likely to appear. Unreasonable designs were not eliminated, but moved towards the reasonable region so that they reside on the edges of the reasonable design space. This approach did not reduce the number of structural optimizations, but caused increase in accuracy. The reasonable design space approach was also applied by Roux et al. (1998) to truss test problems and by Balabanov et al. (1999) in fitting RS weight equation for HSCT wing. Knill et al. (1999) demonstrated a method enabling the efficient implementation of accurate computational fluid dynamics (CFD) predictions into high dimensional, highly constrained MDO problems like HSCT. They utilized simple conceptual level models to select an appropriate set of intervening functions for which to construct RS models since accuracy with loworder polynomials can be improved by transformations of the real function or its arguments. Supersonic linear theory aerodynamics was used as lowfidelity model to generate data for full term RS approximation. Stepwise regression was then used with this data, and RS models in order to determine the significant terms/variables in calculation of the desired aerodynamic quantities. Computationally expensive Euler solutions were performed next only for set of configurations reflecting the effects of the significant terms. This reduced substantially the number of highfidelity Euler analyses required. They used the reducedterm RS models with Euler analyses as an additive correction RS. Another way for reducing modeling error is to use higher order polynomials or more complex functions in RS. Venter et al. (1996) made use of dimensional analysis for reducing the number of variables and the design space for minimum weight design of plate with a step. The reduction in the number of variables made use of higher order polynomials in RS model affordable. They used cubic and quartic polynomials and improved the accuracy for the stress and buckling constraints. Papila and Haftka (1999 and 2000a) also achieved substantial improvement in accuracy by using cubic polynomials in RS approximation for HSCT wing bending material weight (see in Chapter 5). Neural network (NN) based RS approximations have also been investigated as alternative to traditional polynomialbased RS approximations within the context of design optimization. They allow fitting highly nonlinear response. Their main drawbacks are the cost of fitting the approximation and lack of statistical tools for effective design of experiments. There are comparisons in the literature of the performance of NNbased and polynomialbased RS approximations. For example, Carpenter and Barthelemy (1993) used NNbased and polynomialbased approximations to develop RS for several test problems. They found that two methods perform comparable based on the number of undetermined parameters. Papila (N.) et al. (1999) investigated the effect of data size and relative merits between polynomial and NNbased RS in handling varying data characteristics. They demonstrated that as the nature of the experimental data becomes complex NNbased approximations perform better than polynomialbased. The neural network technique and the polynomial RS were integrated to offer enhanced optimization capabilities by Shyy et al. (1999 and 2001) and Papila (N.) et al. (2001). Vaidyanathan et al. (2000) applied NN and polynomialbased RS approximations in the preliminary design of two rocket engine components, gasgas injector and supersonic turbine. They demonstrated that NNbased and polynomialbased approximations can perform comparably for modest data sizes. Papila (N.) et al. (2000) investigated the performance of both polynomial and NNbased RS in their multi objective turbine blade design optimization problem. Oftentimes, designers sacrifice the fidelity of the analysis model used in data generation for RS approximation due to computational cost. This may cause decrease in accuracy of the data and approximation. There is a simple way to improve accuracy in such cases: correction response surfaces (CRS) that make use of lowfidelity (LF) results with small number of highfidelity (HF) results. Mason et al. (1994), for example, studied design optimization of curved composite frames typical of a support structure for semimonocoque aircraft structures and applied twolevels of fidelity and correction of HF solutions onto the LF ones. Vitali et al. (1997) also used two levels of fidelity and RS technique in design of hatstiffened panel for blended wing body aircraft. Venkataraman et al. (1998) demonstrated effective combination of codes/models at different fidelity levels by CRS for structural optimization of a ringstiffened cylindrical shell of revolution. These works were examples of CRS in component design. Balabanov et al. (1998) used the CRS approach to improve the accuracy of RS based wing weight equations. They used coarse FE models for thousands of HSCT configurations, and a quadratic RS was constructed for wing bending material weight. Then about a hundred of the configurations were optimized with a refined FE model to construct a linear CRS. Both additive and multiplicative corrections were applied. The linear additive correction was found slightly better for their design problem. They tried to correct either the prediction of the RS or to correct the data themselves by the high fidelity optimizations. The correction reduced errors by more than half. Another important relevant work in the sense that LF/HF coupling of numerical simulations applied was by Knill et al. (1999). Their main contribution is that besides CRS they used the RS fitted on LF results to identify the important polynomial terms, so that the RS model of HF had fewer terms and required fewer analyses. Design of experiments offers also minimumbias criterion for reducing modeling (bias) error. Venter and Haftka (1997) developed an algorithm implementing minimum bias based criterion, necessary for an irregularly shaped design space where no closed form solution exist for minimumbias design. They compared minimumbias and D optimal designs for two problems with two and three variables. Minimumbias design based RS was found more accurate than Doptimal for the two problems except for the percent maximum error. Percent maximum error was much bigger by minimumvariance design in Balabanov et al. (1996). Errors in RS approximations represent a part of the uncertainty in the design process. The studies discussed above attempt to increase the accuracy and to decrease the level of uncertainty associated with the use of RS. Another important objective of the present work is try to characterize further the errors in RS predictions so that the uncertainty introduced by them can be quantified. Use of RS in Uncertainty Quantification During the design of a new aircraft concept like HSCT, modifications and changes in geometry and material properties are likely. Therefore, it is important to determine how sensitive to change the design is. Uncertainty can be studied by probabilistic or by fuzzy set models. In both cases, design subject to uncertainty is computationally expensive. Response surface approximations alleviate the computational burden because they allow for an inexpensive evaluation of the effect of changes although they cause uncertainty of some degree themselves (Papila and Haftka, 2001a). Reliability based design studies made use of RS approximations (Schueller et al., 1989; Rajashekhar and Ellingwood, 1993; Romero and Bankston, 1998). Xiao et al. (1999) used a reliability based MDO for shape optimization of a transport aircraft wing, treating material properties, geometric variables and loading conditions as random variables with assumed probability distributions. They constructed RS approximations for the objective function and all the constraints and applied Monte Carlo simulation for calculating reliabilities. Qu et al. (2000) investigated the use of RS combined with efficient Monte Carlo simulation techniques in reliabilitybased design optimization of angleply laminates subjected to mechanical loads at cryogenic temperatures. Uncertainties for the problem were in material properties due to short of characterization through the entire operating temperature range, its failure modes at cryogenic temperatures, stiffness and strength properties manufacturing dependent. Response surface approximations were used to reduce computational cost and filter out random sampling noise in the Monte Carlo simulation. Fuzzy set model of uncertainty and RS approximations, and within the vertex method (Dong and Shah, 1987) were employed by Venter and Haftka (1998b), for calculating the possibility of failure in droppedply composite laminate design. Same approach was also used by Papila and Haftka (1999) for sensitivity to uncertainty of wing bending material weight in HSCT designs (See Chapter 5). CHAPTER 3 RESPONSE SURFACE METHODOLOGY Response surfaces (RS) approximate numerical or physical experimental data by an analytical expression that is usually a loworder polynomial. Khuri and Cornell (1996, p. 3) described the three key steps within the context of the methodology: Preselecting the pointsdesign pointswhere experiments are performed. This step is called design of experiments to represent the design space so that the points selected will yield adequate and reliable measurements/calculations of the response of interest. Determining a mathematical approximation model that best fits the data collected/generated from the set of design points of the previous step. This is performed by conducting appropriate statistical test of hypotheses concerning the model's parameters. Using the RS approximation for predicting the response for a given set of the experimental factors/variables. A common application is to determine the optimal settings of the factors that produce optimum (maximum or minimum) value of the response. This chapter summarizes first the design of experiments used in this dissertation. Section starting on page 36 presents the standard leastsquares fitting technique. In the section starting on page 39, statistical tools to check the accuracy of the RS are given including a customized lackoffit test (page 52). Then, section on page 58 presents Iteratively Reweighted Least Squares (IRLS) accommodated with the standard weighting function. Design of Experiments Experimental designs possess the following properties required to fit a pth order polynomial RS for k variables/factors: at leastp+l levels of each factor (k +l)(k + 2).....(k + p) at least (k + 1)(k + 2) (k + distinct design points p! Montgomery and Myers (1995, p. 280) summarized important characteristics that make choice of experimental design effective. An efficient design will Result in a good fit of the model to the data Give sufficient information to allow a test for lack of fit Allow models of increasing order to be constructed sequentially Provide an estimate of "pure" experimental error Be insensitive (robust) to the presence of outliers in the data Be robust to errors in control of design levels Be cost effective Provide a check on the homogeneous variance assumption Provide a good distribution of prediction variance Although not all of the above properties are required in every RS experience, most of them must be given serious consideration on each occasion where one designs experiments. In many practical situations the scientist or engineer specifies strict ranges on the design/configuration variables (experimental factors). That is region of interest and region of operability are the same and a cuboidal. In general, factors of the experiments or configuration variables v, in this work have different order of magnitudes and different units of measurement. To avoid illconditioning and simplify the numerical calculations used in the parameter estimates for the response surface approximation, variables are coded by Eq. (3.1) so that the smallest and the largest limits for configuration or control variables are assigned as 1 and +1, respectively. v, [max(v,) + min(v, )] / 2 x = (3.1) [max(v ) min(v )] / 2 For the application problems of cuboidal design regions in this study, different experimental designs and their combinations were considered based on the number of factors or configuration variables and the degree of polynomial used. As observed from the literature for the engineering design problems studied in this work quadratic polynomial RS approximations are commonly used (e.g. Balabanov et al. (1996, 1998 and 1999). Cubic or even higher order polynomials, on the other hand, were also applied (e.g. Venter et al., 1996); Papila (N.) et al. 1999). In this work, quadratic polynomials are used as an initial model and based on the accuracy achieved first order and cubic polynomials are also studied. Threelevel full factorial designs (FFD) and facecentered central composite designs (FCCD) are typically applied for quadratic models depending on the number of factors involved. These experimental designs are schematically shown in Figure 2 for the case of three coded factors or configuration variables. In the cuboidal region, it is important that the points be "pushed to the extreme" of the region and it is important for the region to be covered in a symmetric fashion. These properties result in attractive distribution of prediction variance (Montgomery and Myers, 1995, p. 312313). Threelevel FFD and FCCD accomplish this. ,, < ; ". ' S A ' O ,   r  SI I I x X2 x X2 x1 X2 X1 Figure 2: Threelevel full factorial (left) and face centered central composite (right) experimental designs for three codedconfiguration variables With the increase in the number of factors, initial experimental designs considered here move from the full factorial design to FCCD and with more factors added from FCCD to a Doptimal (Khuri and Cornell, 1996; SAS, 1998) subsets of FCCD. Moreover, the feasibility of the data points and their usefulness in terms of design performance is checked by quick, but representative evaluations where possible. Unreasonable designs are removed in order to reduce the number of designs from the original FCCD to provide feasible or reasonable designs for the Doptimal selection. When higher order (e.g. cubic) polynomials are needed, orthogonal arrays (Owen, 1994) with four or more levels or the addition of a subbox in the original design space is used to provide the required number of levels. Table 2 presents a summary of the experimental designs used in this study for number of factors and degree of polynomials employed. Table 2: Experimental designs used in the present work # of factors Quadratic polynomial Cubic polynomial 2 Threelevel full factorial design Threelevel full factorial design 3 Threelevel full factorial design + subFCCD 5 FCCD + orthogonal array Doptimal subset of 10Doptmal subset of + orthogonal array (FCCD infeasible designs) Fitting Procedure: Standard Least Squares The response at jh design point x,j is denoted as y, [= y(x,)] and given by Eq. (3.2) y, = r7 + C (3.2) where ir [ = (x ) ] is the true mean response and e, represents other random sources of variation such as measurement error and noise not accounted for in r7. In other words if the experimenter has the luxury of performing the experiment many times at x,, the average of the observations will tend to r, despite the random errors j In RS technique the true mean response is assumed to be given in terms of coefficients A s and shape functions f (x ) where i=1, nb as nb rf(x ) = f (x)= fT (3.3) 1=1 where fT ={ (x, ) (X ..) f J (x)} and p = t f, * Bold face is used for vector and matrix representations, e.g. x = x, X2 T for twovariable case where superscript T stands for transpose operation The shape functions are usually assumed as monomials. For instance, for quadratic RS in twovariables Eq. (3.4) shows the shape functions and their vector representation f1(x ) 1 f2(Xi) j f = f, f3(xi) x2j (3.4) A(x,) x f,=f(,) f4 (X ) (34) f6 (X) X2 The methodology (Myers and Montgomery, 1995; Khuri and Cornell, 1996) makes assumptions, socalled ideal error assumptions (or independent and identically distributed, i.i.d. assumptions) The RS expression is exact and represents the true mean response The error e, in the experiments is uncorrelated, normally distributed random noise. This random variable has expected value of zero, E(, )= 0, and constant variance, Var(e,) =a2 throughout the design space. Or, N(0, '21) where E is error vector and I is the identity matrix. The method of least squares is used to estimate the coefficients s. The least squares function L is the sum of squares of the difference between collected response data yj and the assumed true mean response r77 [Eq. (3.3)] L= y,ff (x,) (3.5) 7=1 The function L is minimized with respect to f, s, and leastsquare estimators b,s of f, s satisfy L =2 f(x) y bk( ) =0 i ,2,....,b (3.6) aA k=l This yields a system of linear equations that can be solved for b,s. Then the RS approximation is given as 17b j(x,)= [ bf(xj) 1=1 (3.7) = fTb The difference (residual) between the data y, and the estimate defined in Eq. (3.7) is given for the th point x, as e, = y A(x) (3.8) In general case j(xj) may not be the true mean response and the difference from the data may not be only the random error, therefore the residual error is denoted by e, instead of e . Matrix forms of Eq. (3.5), (3.6), (3.7) and (3.8) for N data points are L = yTy 2pTXTy + PTXTX (3.9) S= 2XTy +2XTXb = 0 (3.10) y =Xb (3.11) e =y Xb (3.12) where X is the matrix whose terms in the row associated with the point x, are formed by shape functions f (x,). For instance, for a quadratic model in twovariables with N data points [see Eqs. (3.3) and (3.4)] X is given as fT 1 x1 X21 fT 1 12 X22 X= fT 1 X JN 1 XlN X2N The system of equations, Eq. (3.10) may XTXb= XTy 2 Xl 1] ]X2 X12 X12 22 2 X2N XN X2N be written as and the coefficient vector b can now be expressed as b = (XTX)XTy (3.15) If Eq. (3.14) is substituted in Eq. (3.9) the remaining error vector e, is found to satisfy ee, = yTy TXTy that is called as sum of squares of the residuals SSE. (3.16) Adequacy and Accuracy of RS Adequacy and accuracy of the RS approximation are mainly affected by the following three factors whose effects also interact with one another Use of finite number of data points Noise in the data Inadequacy of the fitting model (3.13) (3.14) If the true mean response is exactly represented by the RS form as assumed, the components b,s of vector b in Eq. (3.15) are the unbiased estimates of the A s, in other words the expected value of b is equal to P as can be shown as E(b) = E[(XTX) XTy] =E[(XTX) XT(X + E)] (3.17) = E[(XTX) XTX + (XTX)' XTE] Since E(E) = 0 and (XTX)XTX = I (I is the identity matrix) E(b)= P (3.18) However, when a finite number of design points are used the b,s will not be exactly equal 8,s due to effect of random error in the data. The variance of a random variable Y, Var(Y), is the expected value of the squared difference between the variable and its expected value: Var(Y) = E[(Y E(Y))2] (3.19) Similarly the variancecovariance matrix of a random vector Y, Var(Y) is defined as Var(Y) = E[(Y E(Y))(Y E(Y))T ] (3.20) Following Eq. (3.20) the variancecovariance of the coefficient estimator vector b is derived in Eq. (3.21) Var(b) = E{[b E(b)][b E(b)]T = E[(XTX)'XTy _ ][(XTX) XTy ]T} = E{[(XTX) XT(Xp + ) p][(XTX)IXT(Xp + E) p]T = E{(XTX) XT X(XTX)1 (3.21) = (XTX) XTE(ET)X(XTX)I = (XTX)1XTVar(E)X(XTX)1 = (XTX) XTo2INX(XTX) = 02(XTX)1 The variancecovariance matrix is symmetric and its ith diagonal (as i =1,2...n ) characterizes the variance of estimator b, that causes variance of j(x,) socalled prediction variance. The prediction variance can be expressed as [see also Eq. (3.7)] Var L(x,)] = Var[fb] = ff Var(b)f, (3.22) where f, = f(x) Substitution of Eq. (3.21) into Eq. (3.22) gives prediction variance in Eq. (3.23). Var[I(xj)]= 2fT (XTX)lf, (3.23) The random noise error variance o2 is not usually available and needs also to be estimated. Sum of squares of the residuals, SSE from Eq. (3.16) is used to develop an estimator s2 (error mean square) for '2. s2 =MSE SSE ee (3.24) Nn, Nn, Sum of squares of the residuals SSE and evidently s2 represents variation in the data not accounted for by the fitting model. The error mean square s2 is an unbiased estimate of the noise error variance U2 provided that fitting model exactly models the true mean response 17. Replacing 72 with s2 estimates for the variancecovariance matrix [Eq. (3.21)] and prediction variance [Eq. (3.23)] can be given as Var(b) = s2(XTX)1 Var[5(xj)], s2ff (XTX) 1 Positive square root of the prediction variance is usually used as an estimate of the prediction error at a design point x,, also called estimated standard error as given in Eq. (3.26). ee (x)= Var[(x )] sff(X TX) f (3.26) where estimator s is positive square root of MSE [Eq. (3.24)] often called in engineering literature as rootmeansquareerror (RMSE). The estimated standard error shows dependence on the location of the design point. Furthermore, as Eq. (3.3) is, in general, only an approximation to the true function, s will contain not only noise error, but also modeling (bias) error. The estimator s is also used to calculate coefficient of variation as given in Eq. (3.27) s cV = (3.27) V y N zYJ where y =  N Besides s and c,, the quality of the approximation is often measured by comparing the sum of squares of the residuals SSE (sum of squares unaccounted for by the fitted model) and the sum of squares due to regression SSREG (sum of squares explained by the fitted model). As SSE is given again by Eq. (3.28) the latter sum of squares is defined by Eq. (3.29) N SSE = (Y, )2 (3.28) J=1 N SSRE = 2 (3.29) J=1 where yj = )(x ) and y is the average of the data as used in Eq. (3.27). Their sum shows the total variation in the data set. It is called as total sum of squares SSr and can also be expressed as N SST = (y, y)2 (3.30) J=1 The matrix forms of the sum of squares are given in Eq. (3.31) SSE =yTy bTXTy SSREG =bTXTy ( )2 (3.31) T (1Ty)2 SS, =y y N where 1 is a vector of N ones. Hypothesis Testing in RS Methodology Statistical hypothesis testing (Stark, 2001) is formalized as making a decision between rejecting or not rejecting a null hypothesis on the basis of a set of observations/calculations. Two types of errors can result from any decision rule (test): rejecting the null hypothesis when it is true (a "Type I error"), and failing to reject the null hypothesis when it is false (a "Type II error"). For any hypothesis, it is possible to develop various decision rules. Typically, the chance of a "Type I error" that researcher is willing to allow is specified ahead of time. The chance that the test erroneously rejects the null hypothesis when true, is called the significance level of the decision rule (a). The observations are then used to calculate a test statistic. Researcher compares the test statistic (such as Fstatistic) and its statistical table value (Fdistribution) associated with the selected significance level a If the statistic is larger than the table value, the test is considered as "significant" at level a, and null hypothesis is rejected. The most commonly used level of significance is 0.05 by which researcher accepts 5% probability of making a mistake by rejecting the null hypothesis. There are certain statistical hypotheses and relevant tests to measure the efficiency and reliability of the RS approximation model. These tests also presume the ideal error conditions. The tests are performed on a null hypothesis Ho against an alternative hypothesis Ha. A statistic customized for the testing is calculated and checked against a relevant statistical table of the statistic as a function of degrees) of freedom that describes how many redundant results/measurements exist in an overdetermined system. One important test is the test of the significance of the fitted regression equation as formalized by the following hypotheses: Ho :All coefficients (excluding the intercept) = 0 (3.32) H : At least one coefficient (other than intercept) is not zero The test determines if there is a relationship between the response variable and a subset of the regression variables. Rejection of Ho implies that at least one of the variables play role on the response values. The test statistic for this test is the Fstatistic as given in Eq. (3.33). F = SSREG(nb1) (3.33) SSE (N n ) Decision based on this test at a significance level a is made by comparing the F statistic calculated by Eq. (3.33) with the Fdistribution value ., where dfn and dfd are degrees of freedom for numerator and degrees of freedom for denominator, respectively. For the test in Eq. (3.32) dfn=nb1 and dfd=Nnb. The F distribution is the distribution of the ratio of two estimates of variance. The F distribution has two parameters dfn and dfd and may be denoted as F .,. The dfn is the number of degrees of freedom that the estimate of variance used in the numerator is based on. The dfd is the number of degrees of freedom that the estimate used in the denominator is based on. The shape of the Fdistribution F ., depends on values of dfn and dfd. Figure 3 shows two examples for Fdistribution. ., is the Fvalue where the rightmost tail area under distribution F. is equal to a. For example, F005,4,12 and F.05,10,100 are equal to 3.26 and 1.93, respectively (associated tail areas equal to 0.05 are as shown on Figure 3). F4,12 5% of the distribution . is greater than 3.26 O 1 2 3 4 5 F SF'10,100 S\ 5 of the distribution iCL is greater than 1.93 0 1 2 F Figure 3: Sample Fdistributions The value F, can also be defined as the upper 100a percentile of the distribution with degrees of freedom at numerator dfn and at denominator dfd. If the F statistic by Eq. (3.33) is beyond the upper 100a percent point of the distribution with dfn = nb 1 and dfd = N n,, that is if F > Fn _, _n then the null hypothesis is in Eq. (3.32) rejected at the a level of significance. This is interpreted that the variation accounted for by the model (through the values of b, other than intercept) is significantly greater than the unexplained variation (Khuri and Cornell, 1996, p. 31). In general, the test statistic calculated such as F in Eq. (3.33) is different for a null hypothesis than the table value of the selected significance level a. The statistic corresponds to a level in the table calledpvalue. Thepvalue is the smallest significance level for which any of the tests would have rejected the null hypothesis (Stark, 2001). Therefore, it provides a sense of significance of the evidence against the null hypothesis. The lower the pvalue, the more significant the evidence and the test. If the significance level a is set as 0.05, the test resulting in apvalue under 0.05 would be significant, and the null hypothesis is rejected in favor of the alternative hypothesis. In the test of the significance of the fitted regression equation formalized in Eq. (3.32), the pvalue is determined as the rightmost tail area under the distribution Fn _,N beyond the F statistic value by Eq. (3.33). For instance, if test statistic F is 2.50 on FI0,00 (Figure 3), then pvalue would be 0.01 (FOO10,100o =2.50). This means the researcher can reject the null hypothesis at any significance level a larger than 0.01 and conclude in favor of alternative hypothesis. Coefficients of determination R2 and R2 (R2adjusted), accompanying statistics to Fstatistic, are respectively a measure and an adjusted measure of the proportion of total variation of the values ofy, about the mean value of the response y explained by the fitted model (Khuri and Cornell, 1996, p. 3132) and given as SSR R = REG ST (3.34) R21 _SSE (Nnb) a SS /(N 1) An R value larger than 0.9 is typically required for an adequate approximation. Another hypothesis testing concerns the individual coefficients in the model: Ho :p, =0 (3.35) H, :A 0 While the test for Eq. (3.32) is performed for the comparison of the two types of variance, the test for Eq. (3.35) is to check if the mean value for individual coefficients of the model is different than zero. Such a test would be useful in determining the significance of the each variable and regression coefficient in the fitting model. The test is performed by comparing the coefficient estimates to their respective estimated standard errors. The test statistic, tstatistic is defined by t, = b (3.36) S where s, is the estimated standard deviation of ith coefficient b, (called standard error) and is obtained as the positive square root of the estimated ith term of the diagonal (as i= 1,2...nb) of variancecovariance matrix given in Eq. (3.21). For the alternative hypothesis in Eq. (3.35) the test is twosided. Absolute value of the calculated tstatistic is compared with the value obtained from the Student's tdistribution table (denoted as tf ) relevant to the degree of freedom, df=N nb. Student's tdistributions of different degree of freedom df are shown in Figure 4. The tdistribution is a symmetric, bell shaped probability distribution, similar to the standard normal curve. It differs from the standard normal curve, however, in that it has an additional parameter, called degrees of freedom, which changes its shape. The table value t, for a level of significance a is the upper 100(a/2) S,Nnb percentile of the distribution td For instance, the table values of different degree of freedom indicated on Figure 4 are values concerning level of significance a = 0.05. The area beyond 3.18 for df3, for instance, is equal to a 2 = 0.025. df= C, df= =6 df= 3 a/2=0.025 3.18 196 196 3.18 4.00 2A5 trit 2A5 pvalue/2=0.014 Figure 4: The tdistributions of various degrees of freedom df The statistical significance or pvalue for this test is determined as two times the area under the tdistribution tdf and beyond the absolute value of tstatistic calculated by Eq. (3.36). For example, if one tests the null hypothesis H: i = 0 and calculates statistic t=4.0 with df3, from the relevant tdistribution the pvalue is found 0.028 (= 2x0.014) as shown on Figure 4. This means that, there is a 2.8% chance of making an error to reject the null and say blz0. If researcher chooses to work with a level of significance a more than or equal to 0.028 then the test concludes that there is no significant evidence to accept the coefficient being zero. Standard LackofFit Test Lackoffit test is used as a tool to check if the fitting model is not adequate and it does not contain a sufficient number of terms. In order to distinguish the fitting model and the missing terms subscripts, 1 and 2, respectively, will be used in matrix and 49 function forms in the formulation. The fitting model and the true model are represented as in Eqs. (3.37) and (3.38), respectively. y = Xl1P + E (3.37) E(y) = X101 + X2,2 (3.38) where X1 and P1 respectively denote the model previously denoted by X and P in the second section of this chapter (starting at page 36), and X2, P2 are similar to X1 and P1, but reflect terms present in the true response that are not included in the fitting model. For example, for two variables with a quadratic fitting model, and a cubic true model, response and its expectation at design point x1 = [xs x2s J are given as Y, = A1 + A2X, + A3X2 +A 14 + A5X ~X2 + 2 + (3.39) E(y,) = A1 + A")12 + A31 + A4 X32 + A14 5 X X2 + 2 16X2 + + (3.40) 2 I 1 2 13 2 2 ]2 2 0 23 jly 2 24 2] Then, the coefficient vectors (P1 and P2), matrices X1 and X2 are given in Eq. (3.41) (3.42), and (3.43), respectively. 1 = 1, A12 A13 A14 ]1A5 A16]T (3.41) 2= iA2 g (3.41) P2 = L[21 022 23 24 T 1 X11 X21 X X11X21 X2 11 11 21 21 12 2 1 X12 X22 X12 X12 X22 X22 X = 2 2 (3.42) 1 xlJ x2 X xJ XlX X2 X2N X 2N X1 x2N 2 1 X2A 2N 3 2 2 3 x11 11 21 X11X21 21 3 2 2 3 X12 1222 x12 22 X22 X2 3 2 2 3 (3.43) xJ3 x x2J lJ X2J X2J 3 2 2 3 X1N lXN 2N 2N X2N The lack of fit test supplies statistical evidence that the fitting model is unable to account adequately for the variation in the observed response. The test tries to isolate the portion of the residual variation unaccounted for by the fitted model. For that purpose, evidence is sought based on estimator of pure error variance that does not depend on the form of the fitted model required. An ideal scenario for obtaining such an estimate of pure experimental error variance is to have n, replications at several design points x,. Therefore, response forth design point among M distinct design points at its kth replication can now be represented as yk= y(x )k = 7(xJ )+ k (3.44) (wherej=l, M and k=, n,). The total number of design points is still denoted by N and equal to M N = n (3.45) J=1 For the fh design point the average yj is an estimate of the true mean response r/(x ) nj ZYjk y, k=l (3.46) n J that is equal to the true mean response when n > oo. By these estimates, the sum of squares due to pure error SSpE can be calculated as for the design points M nj SSPE = (Jk )2 (3.47) j=1 k=1 It is compared to the sum of squares due to lackoffit, SSLOF using the difference between the estimate y given by the model and y . SSLOF = i ,() 2 (3.48) Equations (3.47) and (3.48) can also be rewritten in matrix forms M 1 1T SSPE nI ni nJ (3.49) y=l n] M 11T SSLOF =yT[IN XI(XTX1) 1X]y y ni y (3.50) where yn is the response vector of n, replications atjth design point, and 1, is a (n, x 1) vector of ones. Note that Eq. (3.50) is the difference between total sum of squares SSE [Eq. (3.31)] and pure error sum of squares SSPE [Eq. (3.49)]. Therefore, an Fratio is formed by using these partitions of residual sum of square F = LOF dLOF (3.51) SSP / d PE PE where dLOF = M nb and dpE = N M are degrees of freedom associated with SSLOF and SSpE, respectively. Lack of fit can be detected at the a level of significance if the F ratio exceeds the table value, Fa,Mn_,NM where the latter quantity is the upper 100a percentile of the central Fdistribution (Khuri and Cornell, 1996, p. 3841). Similar to hypothesis testing for Eq. (3.32) the pvalue is determined as the rightmost tail area under the distribution FMnb,_, beyond the F statistic value by Eq. (3.51). LackofFit with NonReplicated Design Points Standard lackoffit test described in previous section requires exact replications at several of the design points. Unlike physical experiments, exact replications by numerical experiments/simulations do not supply additional information since they result in exactly the same response when repeated. It is still possible, however, to obtain an estimator that approximates the notion of pure experimental error (Hart, 1997, p. 123). The idea is to treat the observations at neighboring design points as "near" replicates. Then test procedures are based on a pseudo error estimator of the error variance. Neill and Johnson (1985) generalized the standard pureerror lack of fit test to accommodate the case of nonreplication for first order fitting model. Here, the method is extended to a higherorder fitting model. The method considers a design point as near replicate design point Xjk within a distance 6jk from the associated reference point x, subscriptt 'jk' stands for kth nearreplicate around xj) xjk = X + 6jk (3.52) where distance or disturbance vector (for a general, nvariable case) is 6jk =[ clk 2j k pjk n. jk]T (3.53) The monomial terms at design point xjk for the approximation model are calculated using Eq. (3.52) and (3.53). In a quadratic fitting model, for instance, they are given as in Eq. (3.54) XPJ = X2 p + {&k} x = x, {22p x k + = 2x k +k } (3.54) XpkX,k = XjX, + {Xfpj ,k + x,Ypjk + cpjkC~jk } where, r =1, 2,..., n (number of design variables). The design matrix X1 for the data set augmented by the nearreplicate designs is shown by Eq. (3.55) X, = X + A (3.55) where X1 is the design matrix as if the nearreplicate designs are replaced by their associated reference design points xi so that it represents an experimental design with replications, and A, is the disturbance design matrix that turns X1 into X1. Entries in a row of the disturbance matrix Al correspond to the fitting model coefficient vector P1 in Eq. (3.56). P1 = l1 A2 A3 "' A, "' A,, ]T (3.56) For example the (j+k)th row of the disturbance matrix corresponding to kth near replication of the /h design point is given as (assuming th design point is followed by its k near replicates) A1 j+k (A (A A)j2 (jk3 (Ajk ), (Ajk) b] (3.57) For design points that are not near replicates the entries in the associated row are all equal to zero. The nonzero entries in Eq. (3.57) are formed by the corresponding monomial terms within the brackets [quadratic fitting model case, Eq. (3.54)]. 54 For instance, in a twovariable problem where the data has two near replicates (k =1,2) for design point x, and the fitting model is quadratic, design matrices, coefficient vector and disturbance design matrix are given as in the following equations. 1 ,,1 21 11 X1 1 1 21 1 x12 22 <12 12 22 22 1 2 2 X, = xlJ X2J lJ xlj]2] 2 (3.58) 1 Xll X2]j (XlJ1)2 (xljlx2j1) (X2 1)2 1 x12 2j2 (X1)2 (x2 12 2j2 (XC22)2 1 2 2 XIN X2N X X X X 2N N lN 2N 2N 1 x,,1 1 1 1 1 21 1 1 2 2 1 2 2 X12 22 X2 x1222 x22 XX = 1 l x2j 1 X X2 2 (3.59) 1 2 2 1 X, X2, X X1,X], X2] 1 2 2 1 xN x, 2N X2N X1N X2N X 2N P1 = 11 2 A13 A14 A15 16 ] (3.60) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 Am = (3.61) 0 (A 1)12 (A,1)13 (A,1)14 (Aj1)15 (A1)16 0 (A2 )12 (Aj2)13 (A2 )14 (Aj2)15 (Aj2)16 0 0 0 0 0 0 where for the kth near replication nonzero terms in A, are given by Eq. (3.62). 0 0 0 00 0 O (A jk)12 1 { jk} (AJk )13 {2jk} (AJ k)l4 = {2xl Jlk + k } (3.62) (A k )15 = {Xlj,2k + X2 1,&k + d~lk 2 2jk} (Ak )16 = 2x2 &2Jk + 2jk For the data with near replicates, the response vector can now be represented in the following matrix form after substituting Eq. (3.55) into Eq. (3.37). y = X p + A p + E (3.63) Now a corrected response vector, y can be defined by using the fitting model to project the nearreplicate data to the reference points y = y A i (3.64) y = X11 + (3.65) Equation (3.65) allows a standard lackoffit test since the design matrix X1 now is with replicated design points that allow calculating the pure experimental error variance. The corrected response vector, y, however, is not observable. Neill and Johnson (1985) proposed to replace it with an observable response vector y as given in Eq. (3.66) that, at least asymptotically, is comparable to y. y = yAlb, (3.66) The test formalized in this approach (Neill and Johnson, 1985) is given in Eq. (3.67) H : E(y) = (X + A)p (3.67) H, :E(y) = (X1 + A))Pl + (X2 + A,2)p where X2 and A2 are respectively similar to X1 and A1, but corresponding to the missing term coefficient vector p2. Similar to standard lackoffit test, the test statistic in Eq. (3.51), now denoted as F, is based on sums of square error SSPE and SSLOF [Eqs. (3.49) and (3.50), respectively], but with y, y, and X1 replaced with y, y5 and X1, respectively. Neill and Johnson (1985) claimed that F is asymptotically comparable to the test statistic obtained when replication actually exists. A fivevariable problem was used in order to test the lackoffit test with non replicated data. First a quadratic function in fivevariables is calculated in 126 distinct design points (coded variables between 1 and +1) that will be used for the fivevariable HSCT wing problem (see Chapter 5 for details of the experimental design). The data is contaminated with a random, normally distributed noise (o = 3000). Then, a face centered composite design (FCCD) around the central design with a disturbance size of 0.1 in coded space (42 near replicate designs) was generated, and contaminated results were calculated similarly. Figure 5 shows the clustering nearreplicate designs around the center design in 3D case in order to visualize the procedure. (0.1 0.1,0.1) 0.1,0. ,0.1) (1 1 ,1) 1 ,1) (0 1,0 .1 ) (1, 1 i ' I I I cV3 (0.1,0.1 0.1) 0 CV0 .CV1_. (1,1,1)0C( 0 (1,1,1) Figure 5: 3D representation of clustering design points to generate nearreplicate data for lackoffit test The left cubicle represents the whole design space with distinct points (equivalent to 126 design points in fivevariable example) and the right cubicle represents the 42 nearreplicates around the center design. The lackoffit test for this set of data obtained from a quadratic function should give no significant indication for the inadequacy of fitting model. The same approach was repeated by generating data from a cubic function while keeping the fitting model quadratic. This latter case should indicate significant statistical evidence for the inadequacy of the fitting model. Table 3 shows the results of the two tests. The second column presents the results when the true function is also quadratic polynomial of five variables. For a significance level of a = 0.05, the calculated statistic F is smaller than the relevant value from the F distribution table. That is, there is no significant statistical evidence for rejecting the null hypothesis; the quadratic fitting model is the true model. The pvalue of almost one indicates that the chance of making an error by rejecting the null hypothesisa quadratic model is adequateis about 100%. In contrast, results with a cubic true function indicate that there is a statistical evidence to reject the null hypothesis, and the chance of rejecting it erroneously (pvalue) is less than 1%. Table 3: Verification of the lackoffit test with nonreplicated design points Statistical parameters Quadratic true function Cubic true function Quadratic fitting model Quadratic fitting model M 126 126 N 168 168 nb 21 21 F 0.6176 2.1530 FO.5, nbM 1.5704 1.5704 pvalue 0.9747 0.0030 Iteratively Reweighted Least Squares As mentioned in previous section, one of the factors affecting the accuracy of the RS is noise in the data. While lowamplitude, random noise is filtered well by RS approximations, data with large errorsoutliersmay cause significant loss of accuracy. An example in onevariable case is shown graphically in Figure 6 where two noisy data set and true mean response data are compared. Data labeled as "data with outlier" is the same as the data labeled as "noisy data", but with a larger error at the rightmost design point. While the fit of the original "noisy data" is very close to true mean response function, the slope of the function has changed substantially due to the outlier. 9 S true function 8 O x noisy data 7 0 data with outlier fit on noisy data 6  fit on data with outlier. * 4 3 5 ". / 5 I t^J I 2 0 0 0 2 4 6 8 10 12 x Figure 6: 1D example of the effect of an outlier in the data Outliers are atypical (by definition), infrequent observations. In the example, because of the way in which the regression line is determined (especially the fact that it is based on minimizing not the sum of simple distances, but the sum of squares of distances of data points from the line), the outlier has a profound influence on the slope of the regression line. A single outlier is capable of considerably changing the slope of the regression line as shown in Figure 6. Robust regression employs techniques for statistical methods, such as iteratively reweighted least square (IRLS) fitting (Holland and Welsch, 1977), to detect and remove or weight down outliers. Once detected, outliers can be investigated further for possible mistakes. In order to detect outliers, IRLS was used in the present work. Iteratively re weighted least square procedures start with an initial response surface fitted to data and then assign low weights to points with large errors and refit the data using a weighted least square procedure. The process is repeated until convergence, and weights of the response data that do not fit the underlying model (outliers) tend to converge to small values. This effectively eliminates these points from the fitting process. A graphical representation of the procedure is shown in Figure 7 for a 1D case. 1.4 5 [ Data Data 4.5 1.2 LS IRLS .4 1. I weighting 3.5 O i 3.5 o) Outliers 3 0.8 CM .2.5 ) 0.6 0.4 '1.5  0 0.2 0.4 0.6 0.8 1 X Figure 7: 1D example of IRLS procedure In IRLS procedure, the weight w assigned to a data point can be determined by several weighting functions (Myers and Montgomery, 1995, pp. 667673). Table 4 presents the standard weighting functions considered in this dissertation and collaborative research, Kim et al. (2001) focused on efficient IRLS application for results from optimizations. Table 4: Standard weighting functions for IRLS Name W(ej) Range Tuning Constant Huber's minmax 1 e/s H H=1.0 H s/e I e,/s > H Beaton and Tukey's biweight [1( e /B)22 eI < B B=1.0 1.9 0 e/s > B Biweight function by Beaton and Tukey (1974) is preferred to minmax weighting function by Huber (1964) because it gives zero weight to the outliers and thus the outliers are clearly detected. The weight w assigned to a data point can be formalized as 1 l if e,/s w = B (3.68) 10 otherwise where e, is the residual, see Eq. (3.8), s is the estimate of noise, see Eq. (3.24), and B is tuning constant, usually 1 used. The shape of the biweight function in Figure 8 clearly shows that it penalizes outliers with zero or low weight. Points of r (= ej / s) within the supports of the bell shaped biweight function are called inliers, and outliers are the points outside the range of the supports. 1.6 Huber 1.4 Biweight S Biased 1.2 0 1 0.8 S 0.6 0.4  0.2  I U 0 2 0 4 Residual, r Figure 8: Weighting functions for IRLS (Kim et al., 2001) The weighting function called "Biased" in Figure 8, introduced in Kim et al. (2001), is an alternative weighting function to standard biweight function that takes into account onesided errors due to incomplete convergence in minimization problems. After selecting weighting function the coefficients b,s for IRLS approximation can be found by solving Eq. (3.69). XTWX1bi = XTWy (3.69) where W is a diagonal weighting matrix using the weights from Eq. (3.68) w(el) 0 .. 0 0 w(e2) "'. : W = (3.70) 0 ... 0 w(e )_ When the weights are all one, the solution of Eq. (3.69) reduces to ordinary least square solution. However, in general Eqs. (3.69) and (3.70) are nonlinear equations, and the iterative IRLS procedure employs Eq.(3.71). b b) bz+ X 1 X X1\ '(yXb') (3.71) Normally, the IRLS procedure is used only to eliminate or reduce the effect of errors by assigning them low weights (Holland and Welsch, 1977) that basically leaves them out and is referred as robust regression. However, the procedure is used here for the purpose of identifying outliers and factors causing them so that the repair of erroneous data can be sought rather than removing the erroneous points. The JMP (SAS, 1998) statistical software was used in this work for construction of RS approximations and IRLS procedure. CHAPTER 4 BIAS AND NOISE ERRORS IN RESPONSE SURFACE APPROXIMATIONS Applications making use of RS approximations may suffer in regions of design space where the approximation may be poor. The traditional loworder fitting models such as quadratic polynomials may not always represent the behavior of complex engineering systems. On the other hand higher order models are generally impractical to use in high dimensional problems and designers may have to live with lowerorder models. In this situation, information concerning the design regions where the loworder models may suffer the most would be valuable. In this chapter a measure that can provide such information is derived. Mean Squared Error Criterion As a measure of the error in the approximation, the mean squared error of prediction MSEP defined as in Eq. (4.1) is used in this dissertation MSEP(x) = Ei(x) 7(x)]2 (4.1) where 77(x) and j(x) are the true mean response and the prediction by the fitted model, respectively at x. MSEP is by definition an expected value that would be reached if the number of data used in approximation were unlimited. Equation (4.1) may be rewritten as given in Eq. (4.2). MSEP(x) = EUf(x) Ej(x)]+ [E)(x) (x)]}2 = E[(x) ()]2 (4.2) 2Effj(x) Ej)(x)][E)(x) 7(x)]}+ E[Ej(x) 7(x)]2 The crossproduct term in the above squaring operation is EU (x) Ej(x)][Ej(x) 7(x)]} E (x)j)(x) )(x)77(x)1 = (E)(x))2 + Ej(x)17(x) (4.3) = (Ep(x))2 EF(x)77(x) (EJ(x))2 + Ep(x)77(x) =0 since E7(x) = 7(x) and E(E.(x)) = Ej(x). Therefore MSEP = E[i(x) Ey(x)]2 + [E(x) 7/(x)]2 (4.4) The first term in Eq. (4.4) represents the variance error due to random noise [see Eq. (3.19)] and the second term represents bias error due to inadequate modeling. This error expression is usually integrated over the design space and the integral is minimized by choosing the experimental designs that control the effect of one or both types of error (Khuri and Cornell, 1996, pp. 207247; Venter and Haftka, 1997). Here instead it is attempted to characterize the error in the predictions of an RS approximation already constructed and to determine the design regions where RS prediction may suffer due to either or both types of error. Therefore Eq. (4.4) is used to investigate the variation of MSEP from pointtopoint. The expectation of the predicted response at a given design point x can be expressed as EF(x) = fTE(bl) (4.5) where f, is the vector of shape functions f [see Eq. (3.3)] calculated at point x. The mean and the variation for the coefficient estimates are given as [see Eq.s (3.15) and (3.21)] E(b,)= (XX, 1 X1E(y) (4.6) T ~ Var(b)= (x )'X Var(y)X (X X, )' 2 tXX ( 1 (4.7) where a is the standard deviation of the noise error. The first part of the mean squared error in Eq. (4.4) is equal to the prediction variance at x [due to Eq. (3.19)]. E[(x) Ej(2 Varli(x) = fjVar(bl)fi (4.8) = 2fT (XTXI) f The second part of Eq. (4.4) is the squared error due to the inadequacy of the model used [Ei(x) r/(x)]2 = (Bias[(x)])2 (4.9) Therefore mean squared error of prediction in Eq. (4.4) can be rewritten as MSEP = Var[l(x)]+ (Bias[(x)])2 (4.10) The form of the prediction variance, Var[4(x)], in Eq (4.8) does not change even in the presence of inadequate modeling, since variance is only caused by random noise. However, when bias error is present, the residual error mean square s2 is not an unbiased estimate of a2. Using Eqs. (3.16) and (3.24) we get * Subscripts 1 and 2 assign matrices and vectors for the fitting model and the missing terms, respectively. 2 yT(N P) Np, (4.11) where P = X(XlTX 1 (4.12) The accuracy of this estimate depends on factors such as the available data or number of points and also the adequacy of the fitting model that determines whether the estimate is unbiased or biased. The true response at x can be written as 17(x)= = l 2 (4.13) 1(X) =f 2 0+f2P, (4.13) where f2 are terms missing from the assumed model. Since one usually does not know the true response, it is often assumed to be a higher order polynomial when monomials are used as shape functions. The mean of the true response is given as E(y)= X1 P +X2p2 (4.14) where X2 is similar to X1, but due to the terms (shape functions) present in the true response that are not included in the fitting model. After substitution of Eq. (4.14), Eq. (4.6) becomes E(bl)= 0P +Ap2 (4.15) whereA = (XTX1)1 XX2 is called the alias matrix. Substitution of Eq. (4.15) into Eq. (4.5) yields Ej(x) = flr(p + Ap2) (4.16) (Bias[(x)])2 ={fT (pl + Ap2) fT fj 2 (4.17) = pT[ATfc f2 [fTA fT]p2 The MSEP at a given point can now be estimated as by using s2 instead of 02 MSEP(x)= s2fT(XX1)f1 +p Mp2 (4.18) where M = [ATf1 f2 fTA fT] (4.19) Note that MSEP is an expectation by definition [Eq. (4.1)], and Eq. (4.18) is its estimate by the data available and associated s2. Using E(y) from Eq. (4.14), the expected value of biased error mean square s2 given as (Seber, 1977) 2rP X (IN P)X2 2 E(s2) = 2 (N P) Npl (4.20) since (IN P) is an idempotent matrix, that is a matrix whose square is equal to itself. It is also positive semidefinite matrix, so E(s2) >2 provided that (IN P)X2P2 0. Equation (4.20) means that when bias errors are present, s2 is a biased estimate of 072. The standard lackoffit test detects, when significant, the presence of the second terms in Eqs. (4.18) and (4.20). If this expected value, E(s2), is substituted in prediction variance contribution in Eq. (4.18) expected value for the estimate MSEP can be expressed as E[MSEP(x)] T (X )'lf 02 + 2TK2 + pM2 I ()N p (4.21) 2f1T (XTX,) fi + p Gp2 where K= Xj (X2 XA) (4.22) flT (XTX )lfi G= fT(xTx) K+M (4.23) N p, The aim of this chapter is to identify points where the bias error is large. It is assumed that the terms missing from the fitting model are known, but there is not enough data to calculate the corresponding coefficients P2. If one can estimate the size of p2, it is possible to formulate a constrained maximization problem for the largest magnitude of the mean squared error predictor that may be experienced at any given design point for the worst possible P2 of that magnitude. max E[MSEP(x)] P2 (4.24) such that C 2 = C The Lagrangian for this optimization problem can be written as L(2, A) = E[MSEP(x)]+ A(p 2 c) (4.25) Differentiating the Lagrangian with respect to P2 V[o2frT(XX1)f i+ V(TGi2)+ V T2 c) 0 (4.26) where V = a a] Equation (4.26) yields the following I 3/21 3f22 32 p, . eigenvalue problem at a design point x ; G02 +/V2 =0 or Gp2 GP2 =0 (4.27) for which the maximum eigenvalue /Gmax characterizes the maximum possible mean squared and bias error associated with the assumed true model that includes shape functions missing in the fitting model. The corresponding eigenvector defines the coefficients of the missing shape functions that results in the largest bias error when fitted only with the assumed model. The eigenvectors and the function experiencing the worst possible bias error may be different pointtopoint although the magnitude of the missing coefficient vector is constrained. So the eigenvalue calculated does not reflect the true function corresponding to the data (as the data is insufficient to calculate P2 ). It reflects instead, the assumed form of true function with the shape function coefficients P2 (among all the possible combinations such that PIT2 = c2) causing the largest error. Estimated Standard Prediction Error and Residual Based Error Bound The eigenvalue measure characterizes the largest possible mean squared error at a design point for which noise also contributes as well as the insufficient fitting model. The measure, on the other hand, does not pay any attention to the actual data, and so does not include actual noise in the data. The matrix G defined in Eq. (4.23) suggests that associated eigenvalues primarily characterize the error due to missing terms (shape functions) of the true function and the effect of the experimental design. As a result, the characterization is qualitative only, and may be used to locate design points of possible large approximation error. The success of the eigenvalue measure to pinpoint the locations of high error regions may be determined by checking its correlation with the residuals. The actual error at a design point x, associated with the RS approximation y(x,) needs the true mean response r7(x,) that is not always available even for data points due to noise. For engineering problems, such as wing structural weight as studied in Chapter 5, the true function is also not known. In order to test the use of eigenvalue measure for such problems, another measure based on the available data, approximation and associated estimated standard prediction error was used. This alternate measureconservative error boundwas also used as an approximate tool to quantify uncertainty due to the use of the RS approximation )(x ). Recalling MSEP expression in Eq. (4.4) and recognizing the first term as being the prediction variance, the problem of calculating an estimate of the MSEP at any design point stands on the second term described in Eq. (4.9) as also repeated in the following (Bias [y(x )])2 = [E5(x,) r(x, )]2 (4.28) Equation (4.28) cannot be calculated exactly since neither terms on right hand side are available. To be able to determine approximately, Eq. (4.28) was replaced by (Bias (x)])2 = [(x,) y(x )]2 (4.29) where right hand side is square of residual e, also defined in Eq. (3.8). While y(x ) is available at any design points, y(x ) is only available at data points, and additional experiment/simulation is needed for any other design point. Since the alternate error measure was also intended to be available for any design point, the right hand side was replaced by an RS approximation fitted on available data of residual es. Two measures of error are collectively used for the data of residual related RS, Absolute value of residual relative to predicted response at a design point x, eR(xJ) = e (4.30) Equation (4.30) measures the combined effects of modeling and noise error. Average of the absolute relative residuals in Eq. (4.30) N SeR (x ) eR = (4.31) R N ' where N is number of data points. Equation (4.31) is calculated to compare with error data obtained by Eq. (4.30) at the design points. Since the error is not completely characterized by these two measures first the maximum of the two was used as the residual error measure in order to be more conservative in modeling uncertainty analysis eRmax (X,) = max{eR (x ), eR (4.32) Then an RS approximation eR (x) was constructed for eRmax(x,) to predict the result of Eq. (4.29) by Eq. (4.33) (Biask(xj) = [ R(Xj)(xj)]2 (4.33) An estimate of error at a design point x, was then defined by squareroot of an alternate estimate of MSEP that is sum of prediction variance as first term in Eq. (4.18) and prediction by Eq. (4.33) at x, s = AMSEP(x,) = s2f(XXTi) fl + R(xJ)xJ)]2 (4.34) Equation (4.34) was treated as a conservative error or residual bound that may be suffered in the predictions by f(x ), and can be calculated pointtopoint. Polynomial Examples The eigenvalue measure was derived as a tool to determine the regions of potential high bias error when the fitting model is missing some of the shape functions of the true model. For instance, the traditional loworder polynomial fitting models are widely used since the higherorder models may not be applicable due to insufficient data. Therefore, the ideal case for demonstrating the use of eigenvalue estimate of bias error is to check the measure for problems where the true function is a polynomial of a higher degree than the fitting model used. The true function for the examples was chosen as a cubic polynomial while the fitting model was kept as quadratic. In other words for each example, the data were generated by a cubic polynomial and then fitted by quadratic model. In order to see the effect of noise, data contaminated with random and normally distributed noise were also studied. Two and fivevariable polynomials were used to demonstrate the use of the eigenvalue measure for different dimensional problems. Coefficient of Correlation As a statistical indication for the success of the measure to pinpoint the high bias error regions in the design space, coefficient of correlation between the squareroot of eigenvalues calculated and the bias error (or residual calculated as the absolute difference between RS prediction and true response value) and its graphical representation can be used. The most widelyused type of correlation coefficient is Pearson r, also called simple linear correlation. The Pearson coefficient of correlation, r between two variables zi (e.g. eigenvalues) and z2 (e.g. absolute bias error or residual) is calculated as S(ziz2), Z1Z2 r = n (4.35) 1 2 where z1 and z2 are the averages, and oa1 and o'2 are the standard deviations in the set of values for variables. It is a measure of the linear relation between two variables and can range from 1 to +1. A value of 0 represents a lack of correlation between the variables. The value of 1 represents a perfect negative correlation while a value of +1 represents a perfect positive correlation. The correlation coefficient r represents the linear relationship between two variables, and if squared, the resulting value r2 is equal to coefficient of determination [R2 in Eq. (3.34)] as if one treats one of the variables e.g. z2 being fitted as a linear function of the other, zl. The sample plots for r are shown in Figure 9. r=0 Z2 S Z2 Z2a r < 0 Zl r=1 Zl Z1 r>0 r=+1 z2i z2 Figure 9: Linear correlation examples (Lowry, 2001) The coefficient of determination r2 represents the proportion of common variation in the two variables (i.e., the "strength" or "magnitude" of the relationship). On the other hand the correlation is not intended to explain or explore cause and effect type relation, and so the fit is not used to predict one from the other. In order to evaluate the correlation between variables, it is important to know "magnitude" or "strength" as well as the significance of the correlation that is the significance level at which one can reject the null hypothesis that there is zero correlation, H : r = 0 and accept that calculated nonzero correlation is not by pure chance. The tstatistic for this test is given as t = (4.36) (1 r2)/(N 2) where N is the number of data (or sample size) to check on the correlation between the variables. It is compared with the ttable value corresponding the significance level, a and degree of freedom (N2) as explained in section concerning hypothesis testing (page 43). Thepvalue can also be determined from the ttable. The significance orpvalue in this dissertation was calculated by a JAVA script posted on the web by Richard Lowry (2001). High positive correlation between eigenvalues and residuals indicates that if designer moves from a design point to another one where eigenvalue is higher, the error due to approximation may also increase. In other words, eigenvalue measure can warn and make the designer cautious against inaccurate regions of approximation. Cubic Polynomial in Two Variables without Noise First a cubic polynomial in two variables was used as the true model, and a quadratic polynomial as the fitted model with no noise in the data. y(x) = (x) = fTl + fT2 (4.37) where x = [ x, ]T f T 3 2 2 ] 1 f X = [x x~2 x x2 x2] (4.38) P, = A i A2 t A A4 A5 A6 i]T 02 = 21 322 )323 24 i]T Data for a threelevel design as shown in Figure 10 is assumed to be available to form the design matrices X1 and X2 [Eq.s (3.42) and (3.43), respectively] that is, there is not enough data points to characterize all ten coefficients. 1,1 1,1 X1, 2 1,1 1,1 Figure 10: Twovariable threelevel factorial design points Equation (4.27) depends only on the experimental design and the fitting and assumed true models, not the response data. The eigenvalue problem in Eq. (4.27) with P2 2 =1 was solved, at 21x21 mesh points over the design region. The maximum eigenvalues and corresponding eigenvectors were determined. Figure 11 shows a contour plot of the squareroot of the eigenvalues. It is seen that high errors are expected at center of the boundaries of the design region: (1,0), (1,0), (0,1) and (0,1), and low values at (0,0), (0.8,0.8), (0.8,0.8), (0.8,0.8) and (0.8,0.8). 1 2VGmax 0.8289 0.5 0.8034 0.7778 0.7522 0.7266 S0.7011 0.6755 ," C : 0.6499 0.6243 0.5987 0.5732 0.5476 0.5 0.5220 S0.4964 0.4709 0.4453 tT .5 0'.5 x1 Figure 11: Squareroot of maximum eigenvalue VGmax contours from Eq. (4.27) for a twovariable cubic polynomial, fitted with a quadratic model (p, =l11 1,12 13 1814 15 J16]T) Next, a specific example with P = [50 3 5 1 2 1]T was used, and sets of 02 that are the eigenvectors associated with the maximum eigenvalues obtained at nine data points in Figure 10 were studied. Note that such an eigenvector at a design point (with P2T2 = 1) defines the cubic polynomial that introduces the largest bias error at that point, when fitting model is quadratic. Here, without noise, absolute true residual at design point x,, le1 is the bias error as given in Eq. (4.39) leT, = 7/(x,) )(x,) (4.39) There are four eigenvectors to investigate since the other five are negative of either one of the four eigenvectors, and since the absolute residuals are identical when only the sign changes. Since the eigenvalue contours consist of different third order polynomials instead of a single one, exact agreement cannot be expected for any single polynomial. Instead, a given polynomial is expected to have a positive correlation between the absolute true residuals e and the squareroot of the maximum eigenvalues Umax Table 5 shows the polynomial coefficients and the associated coefficient of correlation between eD where i denoting the polynomial/eigenvector and ma of Eq. (4.27). Figure 12 shows the absolute true residual plot for the polynomials studied. Comparing Figure 11 and Figure 12, it is seen that each polynomial achieves maximal errors in one of the regions predicted by Figure 11 and at associated design point. For instance, Polynomial 1 is the eigenvector determined at design points (1,1) and (+1,+1), where errors are high as shown in Figure 12a. 77 Table 5: Coefficient of correlation between the absolute residual e T and the squareroot of maximum eigenvalues mx from Eq. (4.27) for twovariable cubic polynomial example without noise (j = [50 3 5 1 2 1]T) i Corresponding design point (X1, x2) f21 022 23 024 1 (1,1) and (+1,+1) 0 0.707 0.707 0 0.797 2 (1,+1) and (+1,1) 0 0.707 0.707 0 0.797 3 (1,0) and (+1,0) 0 0 1 0 0.406 4 (0,1) and (0,+1) 0 1 0 0 0.406 Maximum of four 0.973 0 E a) (a) e T I 04870 04548 o' 04225 03903 03580 0 3258 0 2935 02613 ' 0 2290 01967 0 1645 S01322 01000 00355 0 00032 i 3 06623 06184 05746 0 5307 04868 04430 0 3991 0 3553 03114 0 2675 0 2237 01798 01360 00921 S00482 00044 Vq e^ T, 2 04870 04225 03903 0 3580 0 3258 0 2935 02613 0 2290 0 1967 0 1645 0 1322 00677 0 00355 0 0032 _0'5 015 (b) eTJ 4 06623 06184 05746 05307 04868 04430 03991 03553 03114 02675 02237 01798 01360 00921 00482 00044 1 5 ... I s 0 1 5' (c) (d) Figure 12: Absolute residual ej contours twovariable cubic polynomial example (polynomials defined in Table 5), without noise. a) Polynomial 1; b)Polynomial 2; c) Polynomial 3; d) Polynomial 4 Comparison of absolute residual for the polynomials and the squareroot of eigenvalues can also be shown by Figure 13. The correlations in Table 5 and Figure 13 show that even for a single polynomial the eigenvalue distribution can be used to identify locations of possible high bias error. 0.8 Polynomial 1 0.7 Polynomial 2 A Polynomial 3 Polynomial 4 , 0.6 , 0.5 Ml 0.3 0.2  0.1 .. 0 A . : ;_ 0.4 0.5 0.6 0.7 0.8 0.9 VTGmax Figure 13: Comparison between squareroot of eigenvalues 2max and absolute residual ej twovariable cubic polynomial example (polynomials defined in Table 5), without noise. Figure indicates that eigenvalues correlate well with maximal errors. See correlation coefficients in Table 5 When the maximum of the absolute residual among the four polynomials, e max = max(e1 , e e e2 3 ) was used, the correlation is 0.973 as reported in last row of Table 5. With the maximum of all four, contours given in Figure 14.b are very similar to those of Figure 11 (also repeated as Figure 14.a). Highcorrelation r=0.973 between mx and e max can also be visualized by Figure 15. max wMPP" 05 05 0 8289 08034 07522 07266 07011 06755 0 6499 06243 05987 05732 05476 05220 04964 S 04709 S04453 NMw 05 0 5 ^eT, max 06184 05746 04868 04430 03991 0 3553 03114 02675 02237 0 1798 01360 00921 00482 00044 V IV__. W, 1t1 05 0 05 i1 5 0 '5 xl (a) X(b) Figure 14: Comparison of squareroot of eigenvalues Jmax and maximum absolute residual ej m contours (maximum of residuals among the four polynomials defined in Table 5) twovariable cubic polynomial example, without noise. a) Eigenvalue field; b)Maximum true residual field 0.8 0.7 0.6 0.5 0.1 . 0 0.4 0.5 0.6 0.7 0.8 0.9 vAG max Figure 15: Comparison between squareroot of eigenvalues max and maximum of true residual errors er, max due to the four polynomial in Table 5 twovariable cubic polynomial example, without noise. Coefficient of correlation is 0.973 Similar contour plots (Figure 16) were also created for the conservative error bound s, of each polynomial in Table 5 as defined in Eq. (4.34). It seems that conservative error bound also shows high error at the design points where eigenvectors 80 were determined to define the four polynomials. For instance, Polynomial 3 is the eigenvector determined at design points (1,0) and (+1,0), where residual and conservative error bound are both high compared to other design points as shown in Figure 12c and Figure 16c. 1 r S1 2 S0 7800 0 7764 05 07646 05 07610 07492 0 7455 07338 07301 07184 07146 07029 06991 06875 0 6837 S0 06721 0 06682 06567 06528 06413 06373 06258 06219 0 6104 06064 05 05950 05909 05796 05755 05641 05600 0 5487 05446 1 5 0 5 1 1 1 ' 5 0 5 . . X, X, (a) (b) 1 S3 1 S S S J 4 S07830 07866 05 07657 07691 07484 07516 07311 07341 07138 07166 0 6965 0 6990 06793 06815 06620 ' 06640 06447 0 6465 06274 0 6290 06101 06115 05928 05940 0 5 05755 0 5764 S05582 05589 S05409 05414 0 5236 0 5239 11 015 X, X, (c) (d) Figure 16: Conservative error bound s contours twovariable cubic polynomial example (polynomials defined in Table 5), without noise. a) Polynomial 1; b) Polynomial 2; c) Polynomial 3; d) Polynomial 4 Variation of s within each distribution is smaller compared to corresponding eT but s, is higher than e for all design points. This can be seen better in Figure 17 that indicates s, determined a bound for the error conservatively for nonoise example although it did not require the true function values. 0.8 Polynomial 1 0.7 Polynomial 2 A Polynomial 3 0.6 Polynomial 4 0.5 eTj, 0.4  0.3 0.2 '. . 0.1  0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 S, Figure 17: Comparison between conservative error bound and absolute residual ej twovariable cubic polynomial example (polynomials defined in Table 5), without noise Cubic Polynomial in Two Variables with Noise The twovariable polynomial example presented in the previous section was repeated, with random, normally distributed noise E of zero mean and standard deviation o = 0.3 in order to simulate more general problems with both noise and modeling error. Table 6 summarizes the coefficient of correlations for the four polynomials used also in Table 5 when noise introduced as given in Eq. (4.40). y(x) = flT + fTp + (4.40) Figure 18 shows the absolute true residual contours when noise is present in the data. Comparison of Figure 12 and Figure 18 shows that the effect of the noise is to shift slightly the location of high residual regions. 82 Table 6: Coefficient of correlations between the absolute residual e 1, and the squareroot of maximum eigenvalues mx from Eq. (4.27) for twovariable cubic polynomial example with noise (o =0.3), (P1 = [50 3 5 1 2 ]T) i Corresponding design point (X1, x2) f21 022 23 P24 1 (1,1) and (+1,+1) 0 0.707 0.707 0 0.590 2 (1,+1) and (+1,1) 0 0.707 0.707 0 0.370 3 (1,0) and (+1,0) 0 0 1 0 0.482 4 (0,1) and (0,+1) 0 1 0 0 0.299 Maximum of four 0.859 1 ' eTj 1 I 08118 07581 07044 06507 05970 05433 04896 04359 03822 03285 02747 S 02210 01673 0 1136 00599 0 0062 " 015 x (a) e T 2 08988 08394 07800 07206 06611 06017 05423 04829 04235 03641 03046 02452 01858 01264 00670 00076 x (b) * eIj 07806 07290 06774 06258 05742 05226 04710 04195 03679 03163 02647 02131 01615 01099 00583 00067 er, 4 0 7617 07112 06608 06104 05600 05095 04591 04087 03583 03078 02574 02070 01566 0 1062 00557 00053 1 J5 05 x (c) x (d) Figure 18: Absolute residual e T contours twovariable cubic polynomial example (polynomial defined in Table 5 and Table 6), with noise (o = 0.3). a) Polynomial 1; b) Polynomial 2; c) Polynomial 3; d) Polynomial 4 1 5 I 05 Comparison curves for absolute residual and the squareroot of eigenvalues for the polynomials are shown in Figure 19. Figure 13 and Figure 19 together show that noise caused the residuals spread more at an eigenvalue level calculated. 1 Polynomial 1 0.9 Polynomial 2 SPolynomial 3 0.8 Polynomial 4 0.7  0.6 . e 1. 0 . 0.4 . 0f 4 0.3 . O.2 I ; ... .. ,. :? *... V ... 0.2 : . 0.1 0  0.4 0.5 0.6 0.7 0.8 0.9 2G max Figure 19: Comparison between squareroot of eigenvalues max and absolute residual e j twovariable cubic polynomial example (polynomials defined in Table 5), with noise (o = 0.3) With the maximum of all four absolute residuals, contours given in Figure 20.b are very similar to eigenvalue contours of Figure 11 (also repeated as Figure 20.a). The slight shift on the locations of maximum residual errors can also been seen in Figure 20.b. This shift is attributed to the noise in the data. Coefficient of correlation r=0.859 between and er max is still high. Figure 21 shows the linear correlation between those two quantities and the increased spread due to the noise effect compared to Figure 15. In spite of the cloudy data it reflects the benefit of eigenvalue measure since high residuals are mostly mapped qualitatively by the magnitude of the eigenvalues. w 05 0 5 _1 i5 0L o 5 S0 8289 08034 05 0 7778 0 7522 07011 0 6755 0 6499  0 6243 05987 05732 0 5476 05220 05 04964 04709 0 4453 1 I 5 015 (a) (b) Figure 20: Comparison of squareroot of eigenvalues max and maximum absolute residual e max contours (maximum of residuals among the four polynomials defined in Table 6) twovariable cubic polynomial example, with noise (o = 0.3). a) Eigenvalue field; b) Maximum residual field 1 0.9 0.8 0.7 0.6 eg 0.5 Smax 0.5 0.4 0.3 0.2 0.1 0 S' Ii: ., , I * *3a * ,?.. ' * s1 0.4 0.5 0.6 0.7 0.8 0.9 AG max Figure 21: Comparison between squareroot of eigenvalues max and maximum of residual errors e, max due to the four polynomial in Table 5 twovariable polynomial example, with noise (o = 0.3). Coefficient of correlation is 0.859 eJ max 08990 08417 07843 07270 06697 06123 05550 04976 04403 03829 03256 02683 S2109 0 1536 0 0962 0 0389 The absolute residuals in Figure 20 are due also to noise, so they should also be compared to the estimated standard error e,, (square root of prediction variance). Table 7 summarizes the coefficient of correlation between the estimated standard error and the absolute residuals for four polynomials. Table 7: Coefficients of correlation between the estimated standard error and absolute true residuals e,1 for cubic polynomial example with noise (o = 0.3), (P1 =[50 3 5 1 2 1]T) i Corresponding design point (x1, x2) i21 Pi22 23 24 S r 1 (1,1) and (+1,+1) 0 0.707 0.707 0 1.127 0.043 2 (1,+1) and (+1,1) 0 0.707 0.707 0 0.543 0.363 3 (1,0) and (+1,0) 0 0 1 0 0.793 0.020 4 (0,1) and (0,+1) 0 1 0 0 0.171 0.433 Maximum of four 0.126 Figure 22 shows the field of estimated standard error [square root of Eq. (4.8)] normalized by a. The rootmeansquareerror estimator s for a is different for each polynomial as reported in Table 7. The distribution of the estimated standard error of the Polynomial i ee, can be pictured as Figure 22 scaled by relevant s from Table 7. The coefficient of correlation between the estimated standard error and the maximum eigenvalues is only 0.097 for all polynomials in this example as both distributions depend only on the experimental design Similar to previous example without noise, Figure 23 shows s, determining error bound conservatively for residuals er The effect of noise on the conservative representation of erT, by sj can be seen by the comparison of Figure 17 and Figure 23. 86 The different magnitudes of biased estimates s for a from Table 7 caused s, being conservative at different levels. Error bound s 1 (Polynomial 1), for instance is the most conservative if used to determine residuals for Polynomial 1 e1 . p. 05 S0 O 8956 08759 08562 08364 08167 07970 07773 07576 07379 07182 06985 06788 06591 06394 06196 05999 05 1 d ,A5 Xi x Figure 22: Normalized estimated standard error e, / a for a twovariable quadratic polynomial model with experimental design of Figure 10 1.2 SPolynomial 1 D Polynomial 2 1 Polynomial 3 Polynomial 4 0.8 " as. .. * 0.4 ; . 0.2  0 . 0 0.2 0.4 0.6 0.8 1 1.2 Figure 23: Comparison between conservative error bound sJ and absolute residual eT twovariable cubic polynomial example (polynomials defined in Table 7), with noise (o = 0.3) 