UFDC Home  myUFDC Home  Help 



Full Text  
MULTIPLE SURROGATES AND ERROR MODELING INT OPTIMIZATION OF LIQUID ROCKET PROPULSION COMPONENTS By TUSHAR GOEL 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 2007 O 2007 Tushar Goel To my parents Sushil and Ramesh, sister Manj ari, and brother Arun. ACKNOWLEDGMENTS This dissertation could not be completed without enormous help from my teachers, family, and friends. While I feel that words would never be sufficient to adequately reflect their contributions, I give a try. I am incredibly grateful to my advisors Prof. Raphael Haftka and Prof. Wei Shyy for their continuous encouragement, very generous support, patience, and peerless guidance. Both Prof. Shyy and Prof. Haftka provided me numerous opportunities to develop and to hone my research and personal skills. I am amazed by their neverending enthusiasm and depth of knowledge, and feel extremely fortunate to have been taught by them. I would like to especially thank my advisory committee members, Prof. NamHo Kim, Prof. Jacob N. Chung, and Prof. Andre I. Khuri, for their willingness to serve on my committee, for evaluating my dissertation, and for offering constructive criticism that has helped improved this work. I particularly thank Dr Kim, for many discussions during our weekly group meetings and afterwards. I feel deeply indebted to Prof. Nestor V. Queipo for a very fruitful collaboration. Not only did he significantly contribute to my work but also he was very supportive and helpful during the entire course of my graduate studies. I express my sincere gratitude to Dr Layne T Watson and Dr Daniel J Dorney for their collaboration and help in my research. I thank Dr Siddharth Thakur, who provided a huge assistance with the STREAM code and suggestions related to numerical aspects in my research work. I also thank Prof. Peretz P. Friedmann and Prof. KwangYong Kim for the opportunities to test some of our ideas. I thank my collaborators Dr Raj Vaidyanathan, Ms Yolanda Mack, Dr Melih Papila, Dr Yogen Utturkar, Dr Jiongyang Wu, Mr Abdus Samad, Mr Bryan Glaz, and Dr Li Liu. I learnt a lot from you and I feel sincerely indebted for your help, both personally and academically. I thank the staff of Mechanical Engineering department, particularly Jan, Pam, and David, for their help with administrative and technical support. I also am thankful to the staff at International Center, library, CIRCA, and ETD for their help with this thesis and other administrative details. I sincerely acknowledge the financial support provided by NASA Constellation University Institute Program (CUIP). I duly thank my colleagues in the Structural and Multidisciplinary Optimization Group and Computational Thermofluids Laboratory for their assistance and many fruitful discussions about all the worldly issues related to academics and beyond. I am highly obliged to Prof. Kalyanmoy Deb and Prof. Prashant Kumar at IIT Kanpur, who gave me very sage advice at different times in my life. They indeed have played a big role in shaping my career. I also thank my colleagues Emre, Eric, Pat, Nick, and Amor, who made my stay at the University of Michigan a memorable one. I am grateful to have true friends in Ashish, Jaco, Erdem, Murali, Siva, Ashwin, Girish, Saurabh, Ved, Priyank, Satish, Tandon, Sudhir, Kale, Dragos, Christian, Victor, Ben, and Palani for lending me a shoulder, when I had a bad day and for sharing with me the happy moments. These memories will remain etched for life. Finally, but not at all the least, I must say that I would never have completed this work, had it not been the unconditional love, appreciation, and understanding of my family. Despite the fact that we missed each other very much, they always motivated me to take one more step forward throughout my life and rej oiced all my achievements. To you, I dedicate this dissertation! TABLE OF CONTENTS page ACKNOWLEDGMENTS .............. ...............4..... LIST OF TABLES ............ ...... ._ ._ ...............12... LIST OF FIGURES .............. ...............16.... LIST OF ABBREVIATIONS ........._.._ ..... .___ ...............20.... AB S TRAC T ......_ ................. ............_........2 CHAPTER 1 INTRODUCTION AND SCOPE ................ ...............31................ Space Propulsion Sy stems ................. ............. ...............3.. 1.... Design Requirements of Propulsions Systems ................. ...............32........... ... System Identification and Optimization: Case Studies .............. .... .......... ..............3 Sensitivity Evaluation and Model Validation for a Cryogenic Cavitation Model ..........34 Shape Optimization of Diffuser Vanes ................. ...............34........... ... Surrogate M odeling .............. ...............35.... Issues with Surrogate Modeling .............. ...............37.... Sampling Strategies ................. ...............37................. Type of Surrogate M odel ................... ............. ...............38...... Estimation of Errors in Surrogate Predictions ................. ...............38............... Scope of Current Research .............. ...............39.... 2 ELEMENTS OF SURROGATE MODELING ................. ...............43........... ... Steps in Surrogate M odeling .............. ...............44.... Design of Experiments ................... .. ......... ...............45...... Numerical Simulations at Selected Locations ................. ...............45........... ... Construction of Surrogate Model .............. ...............45.... M odel Validation............... ......... ....................4 Mathematical Formulation of Surrogate Modeling Problem............... ...............46 Design of Experiments .............. ...............48.... Factorial Designs ................. ...............49................. Central Composite Designs ............... ........ ....... .... ........5 Variance Optimal DOEs for Polynomial Response Surface Approximations ................51 Latin Hypercube Sampling .............. ...............51.... Orthogonal Arrays .................. ............. .....................5 Optimal LHS, OAbased LHS, Optimal OAbased LHS ................. .......................53 Construction of Surrogate Model .............. .... ...............54. Polynomial Response Surface Approximation .............. ...............54.... Kriging M odeling ............. ...............56..... Radial Basis Functions ............. ...............58..... Kernelbased Regression ................. ...............61.__._...... Model Selection and Validation .............. ...............62.... Split Sam ple .............. ...............62.... Cross Validation ............. ...............63..... Bootstrapping .............. ...............64.... 3 PITFALLS OF USING A SINGLE CRITERION FOR SELECTING EXPERIMENTAL DE SIGNS ........._._ .......__. ...............72... Introducti on ........._..... ......_.. .. .... ...............72 Error Measures for Experimental Designs .............. ...............76.... Test Problems and Results ............ ... .... ..._.. ...............80... Comparison of Different Experimental Designs ......._. .......... .. ........._._.....81 Space filling characteristics of Doptimal and LHS designs................. ...............8 Tradeoffs among various experimental designs ........._...... ..... ...._._..................82 Extreme example of risks in single criterion based design: Minmax RMS bias C C D .............. .............. ....... ... ..... .......8 Strategies to Address Multiple Criteria for Experimental Designs .............. .... ........._...89 Combination of modelbased Doptimality criterion with geometry based LHS criterion ................... ............ ...... .. ............. ...............9 Multiple experimental designs combined with pointwise errorbased fi1tering.......92 Concluding Remarks .............. ...............94.... 4 ENSEMBLE OF SURROGATES ................. ...............105............... Conceptual Framework................ .. ... ...............10 Identification of Region of Large Uncertainty ................. ........._._ ...._._ .....10 Weighted Average Surrogate Model Concept............... ...............107 Nonparametric surrogate fi1ter .............. ...............109.... Best PRESS for exclusive assignments ........___..........___ ............... 109 .... Parametric surrogate fi1ter .............. ..... .. ......... ..........10 Test Problems, Numerical Procedure, and Prediction Metrics ......_._ ........... ........ .......112 Test Problems ................ ...............112................ BraninHoo function ................. ...............112......... ...... Camelb ack functi on ........._.._ ........... ...............112... GoldsteinPrice function ................. ...............112......... ...... Hartman functions ................. ......... ...............112...... Radial turbine design for space launch ................. ...............113.............. Numerical Procedure ................. ...............114................ Prediction M etrics ................. ...............115......... ...... Correlation coefficient ................. ...............115................ RMS error............... ...............116. M aximum error ................. ...............116......... ...... Results and Discussion ................... ... .. ...... ...............117..... Identification of Zones of High Uncertainty ................. ...............117.............. Robust Approximation via Ensemble of Surrogates .......... ................ ...............119 Correlations ............ ............ ...............120...... RM S errors ................. ................. 12......... 1.... Maximum absolute errors............... ..... ...... .........2 Studying the role of generalized crossvalidation errors ................. ................ ..123 Effect of sampling density................ ...............12 Sensitivity analysis of PWS parameters ................. ...............126........... ... Conclusions............... ..............12 5 ACCURACY OF ERROR ESTIMATES FOR SURROGATE APPROXIMATION OF NOISEFREE FUNCTIONS ................. ...............144......... ...... Introducti on ................. ...............144................ Error Estimation Measures ................. ............. ... ..........4 Error Measures for Polynomial Response Surface Approximation .............. ................146 Estimated standard error ................. ...............148................ Root mean square bias error............... ...............148. Error Measures for Kriging .............. .... ............. ...............149..... Model Independent Error Estimation Models ................ ...............................15 Generalized crossvalidation error ............. ...............151.... Standard deviation of responses ................ ................. .................. ..152 Ensemble of Error Estimation Measures ................. ...............154........... ... Averaging of multiple error measures ................. ...............154........... ... Identification of best error measure .............. ...............154.... Simultaneous application of multiple error measures ................ .....................155 Global Prediction Metrics............... ...............155 Root mean square error .............. ... ............. ...............155..... Correlation between predicted and actual errors ................. .........................156 Maximum absolute error ................. ...............157............... Test Problems and Testing Procedure .............. ...............157.... Test Problem s ................. ...............157......... ...... BraninHoo function .............. ...............158.... Camelb ack functi on ................. ...............158............... GoldsteinPrice function .............. ...............158.... Hartman functions ................ ...............158................ Radial turbine design problem .............. ...............159.... Cantilever beam design problem .............. ...............159.... Testing Procedure ................. ...............160................ Design of experiments .............. .....................160 Test points .............. ...............161.... Surrogate construction................ ............16 Error estim ation................ ..............16 Results: Accuracy of Error Estimates ................. ...............162........... ... Global Error Measures .............. ...............163.... Pointwise Error Measures ................. ...............165............... Root mean square errors............... ........... ..........16 Correlations between actual and predicted errors .............. ....................16 Maximum absolute errors............... ...............169 Ensemble of Multiple Error Estimators ............ .....___ ...............170. Averaging of Errors ............... ...... ..._ .... ... ..............17 Identification of Suitable Error Estimator for Kriging ............_ ..... .. ..............171 Detection of High Error Regions using Multiple Errors Estimators ................... ..........173 Conclusions............... .... ...........17 Global Error Estimators ............_ ......___ ...............175... Pointwise Error Estimation Models ............__......_ __ ...............176. Simultaneous Application of Multiple Error Measures. ....._____ ...... ...___...........177 6 CRYOGENIC CAVITATION MODEL VALIDATION AND SENSITIVITY EVALUATION ................. ...............199......... ...... Introducti on ................. ... ......... .. .. ....... .. ........ .............9 Cavitating Flows: Significance and Previous Computational Efforts ................... ........199 Influence of Thermal Environment on Cavitation Modeling ................. ................ ..201 Experimental and Numerical Modeling of Cryogenic Cavitation ............... ............._...203 Surrogate Modeling Framework ................. ...............204................ Scope and Organization................. ...... ........0 Governing Equations and Numerical Approach ................. ...............206........... ... Transportbased Cavitation Model .............. ...............207.... Thermodynamic Effects .............. ...............208.... Speed of Sound Model .............. ...............210.... Turbulence M odel ................. ...............211.............. Numerical Approach .............. ...............212.... Results and Discussion ............... .......... .. ........ .............1 Test Geometry, Boundary Conditions, and Performance Indicators.............................21 Surrogatesbased Global Sensitivity Assessment and Calibration ............... .... ........._..214 Global Sensitivity Assessment ................. ...............215.......... .... Surrogate construction..................... .. ... ........1 Main and interaction effects of different variables .............. ....................21 Validation of global sensitivity analysis ................ ............ ........_.........218 Calibration of Cryogenic Cavitation Model .............. ...............219.... Surrogate modeling of responses ................. ....._.. ...._.. ......_.._......220 Multiobjective optimization............... .............22 Optimization outcome for hydrogen .............. ...............222.... Validation of the calibrated cavitation model ........._.._.. ...._... ........_.......222 Investigation of Thermal Effects and Boundary Conditions .............. ....................23 Influence of Thermosensitive Material Properties ........._.._.. ....._.._ ........._......223 Impact of Boundary Conditions .............. ...............225.... Conclusions............... .. .........................2 6 Influence of Turbulence Modeling on Predictions ....._._._ ......._.. .. ...._._.........243 7 IMPROVING HYDRODYNAMIC PERFORMANCE OF DIFFUSER VIA SHAPE OPTIMIZ ATION ................. ...............245................ Introducti on ................. ...............245................ Problem Description ................ ...............247................ Vane Shape Definition .............. ....... .... .... ..................24 Mesh Generation, Boundary Conditions, and Numerical Simulation .........................250 SurrogateBased Design and Optimization ............ ......___....._ ...........25 Surrogate M odeling ............ ..... ._ ...............251... Global Sensitivity Assessment .............. ...............255.... Optimization of Diffuser Vane Performance ............_...... .__ ........._........257 Design Space Refinement Dimensionality Reduction .............. .....................5 Final Optimization............... ..... ..........26 Analysis of Optimal Diffuser Vane Shape .............. ...............260.... Flow Structure .............. ...............260.... Vane Loadings ................. ...............261................ Empirical Considerations .............. ...............261.... Summary and Conclusions .............. ...............262.... 8 SUMMARY AND FUTURE WORK .............. ...............277.... Pitfalls of Using a Single Criterion for Experimental Designs ................. .....................277 Summary and Learnings ................. ...............277................ Future W ork ................. ...............278................ Ensemble of Surrogates .............. ...............278.... Summary and Learnings ................. ...............278................ Future W ork. .............. .... .... ........ ... ......... ............27 Accuracy of Error Estimates for Noisefree Functions ................. ................ ........ .279 Summary and Learnings ................. ...............279................ Future W ork ................. .... ............ .. ........... ..... ..........28 System Identification of Cryogenic Cavitation Model ................ ................. ..........280 Summary and Learnings ....__. ................. ........_.._.........28 Future W ork. ...._.._.................. ... ......_.._..........2 1 Shape Optimization of Diffuser Vanes ....__. ................. ........._.._ ....... 28 Summary and Learnings ....__. ................. ........_.._.........28 Future W ork. ...._.._................. ......._ _. ..........28 APPENDIX A THEORETICAL MODELS FOR ESTIMATING POINTWISE BIAS ERRORS............. .283 DataIndependent Error Measure s ....__. ................. ........__. ........ 28 DataIndependent Bias Error Bounds ....__. ................. ............... 286 .... DataIndependent RMS Bias Error .............. ...............286.... DataDependent Error Measures .............. ...............287.... Bias Error Bound Formulation .............. ...............288.... Root Mean Square Bias Error Formulation .................. ...............289.............. Determining the Distribution of Coefficient Vector (3 .............. .....................8 Analytical Expression for Pointwise Bias Error Bound ...._ .............. ..... ..........291 Analytical Estimate of Root Mean Square Bias Error. ......____ ..... .. .............292 B APPLICATIONS OF DATAINDEPENDENT RMS BIAS ERROR MEASURES ..........294 Construction of Experimental Designs ............_..._ .. ........_._ ... .. ........ .............9 Why Minmax RMS Bias Designs Place Points near Center for Fourdimensional Space? ............ ........ ... ............29 Verification of Experimental Designs ................. ......... ...............296 .... Comparison of Experimental Designs ................. .......... ...............297 .... RMS Bias Error Estimates for Trigonometric Example............... ...............298 C GLOBAL SENSITIVITY ANALYSIS ................. ...............304............... D LACKOFFIT TEST WITH NONREPLICATE DATA FOR POLYNOMIAL RESPONSE SURFACE APPROXIMATION .............. ...............307.... LI ST OF REFERENCE S ................. ...............3.. 10......... ... BIOGRAPHICAL SKETCH .............. ...............329.... LIST OF TABLES Table page 21. Summary of main characteristics of different DOEs. ................ ............. ........ .......70 22. Examples of kernel functions and related estimation schemes. ................ ......................70 23. Summary of main characteristics of different surrogate models............._ ........._._......71 31. Doptimal design (25 points, 4dimensional space) obtained using JMP. ..........................102 32. LHS designs (25 points, 4dimensional space) obtained using MATLAB. .......................102 33. Comparison of RMS bias CCD, FCCD, Doptimal, and LHS designs for 4dimensional space (all designs have 25 points). ............. ...............102.... 34. Prediction performance of different 25point experimental designs in approximation of example functions F;, Fz, and F3 in fourdimensional spaces. ............. ....................10 35. Minmax RMS bias central composite designs for 25 dimensional spaces and corresponding design metrics. ............ ...............103..... 36. Mean and coefficient of variation (based on 100 instances) of different error metrics for various experimental designs in fourdimensional space (30 points). ................... ..........104 37. Reduction in errors by considering multiple experimental designs and picking one experimental design using appropriate criterion (filtering). ............. ......................0 41. Parameters used in Hartman function with three variables. ............. ....................13 42. Parameters used in Hartman function with six variables. ............. ......................3 43. Mean, coefficient of variation (COV), and median of different analytical functions. .........137 44. Range of variables for radial turbine design problem. ...........__.....___ ...............137 45. Numerical setup for the test problems. ........._.._... ...............138._...._... 46. Median, 1st, and 3rd quartile of the maximum standard deviation and actual errors in predictions of different surrogates at the location corresponding to maximum standard deviation over 1000 DOEs for different test problems. ........._.... ........._.....13 8 47. Median, 1st, and 3rd quartile of the minimum standard deviation and actual errors in predictions of different surrogates at the location corresponding to minimum standard deviation over 1000 DOEs for different test problems. ........._... ................1 39 48. Median, 1st, and 3rd quartile of the maximum standard deviation and maximum actual errors in predictions of different surrogates over 1000 DOEs for different test problems............... ...............13 49. Effect of design of experiment: Number of cases when an individual surrogate model yielded the least PRESS error. ............. ...............140.... 410. Opportunities of improvement via PWS: Number of points when individual surrogates yield errors of opposite signs. ................ ......... ........ ......... ...............140 411i. Mean and coefficient of variation of correlation coefficient between actual and predicted response for different surrogate models ................. .....__. ........._.._.. .140 412. Mean and coefficient of variation of RMS errors in design space for different surrogate m odels. .............. ...............141.... 413. Mean and coefficient of variation of maximum absolute error in design space. ................141 414. Mean and coefficient of variation of the ratio of RMS error and PRESS over 1000 DOEs ................. ...............142................ 415. The impact of sampling density in approximation of BraninHoo function. ...................142 416. The impact of sampling density in approximation of Camelback function. ...................... 143 417. Effect of parameters in parametric surrogate filter used for PWS. ............ ...................143 51. Summary of different error measures used in this study. ............ ...... ............... 187 52. Parameters used in Hartman function with three variables. ............ ...... ...............18 53. Parameters used in Hartman function with six variables. ............ ......................8 54. Range of variables for radial turbine design problem. ......................__ ...............188 55. Ranges of variables for cantilever beam design problem. ............. ......................8 56. Numerical setup for different test problems. ............. ...............188.... 57. Mean and coefficient of variation of normalized actual RMS error in the entire design space. ............. ...............189.... 58. Mean and coefficient of variation of ratio of global error measures and corresponding actual RMS error in design space. ............ ...............190..... 59. Mean and COV of ratio of root mean squared predicted and actual errors for different test problem s. ............ ...............191..... 510. Mean and COV of correlations between actual and predicted errors for different test problem s. ............. ...............192.... 511i. Mean and COV of ratio of maximum predicted and actual errors for different test problems............... ...............19 512. Mean and COV of ratio of root mean square average error and actual RMS errors for different test problems. ............ ...............194..... 513. Comparison of performance of individual error measures and GCV chosen error measure for kriging. ............. ...............194.... 514. Number of cases out of 1000 for which error estimators failed to detect high error regions. ............. ...............195.... 515. Number of cases out of 1000 for which error estimators failed to detect maximum error regions. ............. ...............196.... 516. Number of cases out of 1000 for which different error estimators wrongly marked low error regions as high error. ............. ...............197.... 517. Number of cases out of 1000 for which different error estimators wrongly marked low error regions as the maximum error region. ............. ...............198.... 518. High level summary of performance of different pointwise error estimators. ................... 198 61. Summary of a few relevant numerical studies on cryogenic cavitation ............... ...............239 62. Ranges of variables for global sensitivity analyses. ............ ...............240..... 63. Performance indicators and corresponding weights in surrogate approximations of prediction metrics P,f and T,f ............ ...............240 64. Performance indicators and corresponding weights in surrogate approximations of prediction metrics P,f and T,f in modelparameter space. ............. .....................24 65. Predicted and actual P,f and T,f at bestcompromise model parameter for liquid N2 (Case 290C). ............ ...............241..... 66. Description of flow cases chosen for the validation of the calibrated cryogenic cavitation model ................. ...............242....._._. ..... 67. Model parameters in LaunderSpalding and nonequilibrium k E turbulence models.....243 71. Design variables and corresponding ranges. ............ ...............273..... 72. Summary of pressure ratio on data points and performance metrics for different surrogate models fitted to Set A. ............. ...............273.... 73. Range of data, quality indicators for different surrogate models, and weights associated with the components of PWS for Set B and Set C data, respectively. .................. ...........274 74. Optimal design variables and pressure ratio obtained using different surrogates constructed using Set C data. ............. ...............274.... 75. Comparison of actual and predicted pressure ratio of optimal designs obtained from multiple surrogate models (Set C). ............ ...............275..... 76. Modified ranges of design variables and fixed parameters in refined design space. ...........275 77. Range of data, summary of performance indicators, and weighted associated with different surrogate models in the refined design space. ............. ......................7 78. Design variables and pressure ratio at the optimal designs predicted by different surrogates. ............ ...............276..... 79. Actual and empirical ratios of gaps between adj acent diffuser vanes. ..........._..._ ..............276 B1. Design variables and maximum RMS bias errors for minmax RMS bias central composite designs in N,=25 dimensional spaces. ............. ...............303.... B2. Comparison of different experimental designs for two dimensions. ............ ..................303 B3. Comparison of actual and predicted RMS bias errors for minmax RMS bias central composite experimental designs in fourdimensional space. ............. .....................30 LIST OF FIGURES Figure page 11. Schematic of liquid fuel rocket propulsion system. ............. ...............41..... 12. Classifieation of propulsion systems according to power cycles. ................ ................ ..42 21. Key stages of the surrogatebased modeling approach. ............. ...............66..... 22. Anatomy of surrogate modeling: model estimation + model appraisal. ............. .................66 23. A surrogate modeling scheme provides the expected value of the prediction and the uncertainty associated with that prediction ................. ...............67........... ... 24. Alternative loss functions for the construction of surrogate models. ............. ...................67 25. A twolevel full factorial design of experiment for three variables. ............ ....................68 26. A central composite design for threedimensional design space. ............. .....................6 27. A representative Latin hypercube sampling design with Ns = 6, N~ = 2 for uniformly distributed variables in the unit square. ............. ...............69..... 28. LHS designs with significant differences in terms of uniformity. ............ ....................69 31. Boxplots of radius of the largest unoccupied sphere inside the design space [1, 1]N". ........97 32. Illustration of the largest spherical empty space inside the 3D design space [1, 1]3 (20 points). ............ ...............97..... 33. Tradeoffs between different error metrics. ............ ...............98..... 34. Comparison of 100 Doptimal, LHS, and combination (Doptimality + LHS) experimental designs in fourdimensional space (30 points) using different metrics. .....99 35. Simultaneous use of multiple experimental designs concept, where one out of three experimental designs is selected using appropriate criterion (filtering). ........................100 41. Boxplots of weights for 1000 DOE instances (Camelback function). ............ .................128 42. Contour plots of two variable test functions. ............. ...............129.... 43. Boxplots of function values of different analytical functions. ................ .....................130 44. Contour plots of errors and standard deviation of predictions considering PRS, KRG, and RBNN surrogate models for BraninHoo function. ............ .....................13 45. Standard deviation of responses and actual errors in prediction of different surrogates at corresponding locations (b oxplots of 1000 DOEs using B raninHoo functi on)..............13 1 46. Correlations between actual and predicted response for different test problems. ...............132 47. Normal distribution approximation of the sample mean correlation coefficient data obtained using 1000 bootstrap samples (kriging, BraninHoo function). .....................133 48. RMS errors in design space for different surrogate models. ............ .....................13 49. Maximum absolute error in design space for different surrogate models. ..........................135 410. Boxplots of ratio of RMS error and PRESS over 1000 DOEs for different problems. .....136 51. Contour plots of two variable analytical functions. ............. ...............178.... 52. Cantilever beam subj ected to horizontal and vertical random loads. .........._... ........._.....178 53. Ratio of global error measures and relevant actual RMS error. ............. .....................18 54. Ratio of root mean square values of pointwise predicted and actual errors for different problems, as denoted by predicted error measure. ............. ..... ............... 18 55. Correlation between actual and predicted error measures for different problems. .............184 56. Ratio of maximum predicted and actual absolute errors in design space for different problem s. ............. ...............186.... 61. Variation of physical properties for liquid nitrogen and liquid hydrogen with tem perature. ............ ...............228..... 62. Experimental setup and computational geometries. ............ ...............229..... 63. Sensitivity indices of main effects using multiple surrogates of prediction metric (liquid N2, Case 290C). ............ ...............230..... 64. Influence of different variables on performance metrics quantified using sensitivity indices of main and total effects (liquid N2, Case 290C). ................ ......................23 1 65. Validation of global sensitivity analysis results for main effects of different variables (liquid N2, Case 290C). ............ ...............232..... 66. Surface pressure and temperature predictions using the model parameters for liquid N2 that minimized P,,, and Td,, respectively (Case 290C). ............. ....................23 67. Location of points ( C,, ) and corresponding responses used for calibration of the cryogenic cavitation model (liquid N2, Case 290C). ............. ..... ............... 23 68. Pareto optimal front and corresponding optimal points for liquid N2 (Case 290C) using different surrogates. ............ ...............233..... 69. Surface pressure and temperature predictions on benchmark test cases using the model parameters corresponding to original and bestcompromise values for different fluids. ........... ...............234..... 610. Surface pressure and temperature predictions using the original parameters and best compromise parameters for a variety of geometries and operating conditions. .............236 611. Surface pressure and temperature profie on 2D hydrofoil for Case 290C where the cavitation is controlled by, (1) temperaturedependent vapor pressure, and (2) zero latent heat, and hence isothermal flow field. ............ ...............237..... 612. Impact of different boundary conditions on surface pressure and temperature profie on 2D hydrofoil (Case 290C, liquid N2) and predictions on first computational point next to boundary. ............ ...............238..... 613. Influence of turbulence modeling on surface pressure and temperature predictions in cryogenic cavitating conditions. ............ ...............244..... 71. A representative expander cycle used in the upper stage engine. ............ .....................26 72. Schematic of a pump. ............. ...............264.... 73. Meanline pump flow path. ................ ............ ........ ......... ........ .........264 74. Baseline diffuser vane shape and timeaveraged flow. ............ ...............265..... 75. Definition of the geometry of the diffuser vane. ............ ...............265..... 76. Parametric Bezier curve............... ...............265. 77. A combination of H and Ogrids to analyze diffuser vane. ............ .....................26 78. Surrogate based design and optimization procedure. ............ ...............266..... 79. Surrogate modeling............... ...............26 710. Sensitivity indices of main effect using various surrogate models. ............ ..................267 711. Sensitivity indices of main and total effects of different variables using PWS. ...............268 712. Actual partial variance of different design variables. ............. ...............268.... 713. Baseline and optimal diffuser vane shape obtained using different surrogate models. .....269 714. Comparison of instantaneous and timeaveraged flow Hields of intermediate optimal (PRS) and baseline designs. ............. ...............270.... 715. Instantaneous and timeaveraged pressure for the final optimal diffuser vane shape. ......271 716. Pressure loadings on different vanes. ............ ...............271..... 717. Gaps between adjacent vanes. ............ ...............272..... B1. Twodimensional illustration of central composite experimental design constructed using two parameters al and at. ............. ...............301.... B2. Contours of scaled predicted RMS bias error and actual RMS error when assumed true model to compute bias error was quintic while the true model was trigonometric. .......301 Figure B3. Contours of scaled predicted RMS bias error and actual RMS error when different distributions of p(2) were specified. ............ ...............302..... LIST OF ABBREVIATIONS Alias matrix Constants in Hartman functions Estimated coefficients associated with ith basis function Vector of estimated coefficients of basis functions Covariance matrix in kriging Bounds on coefficient vectors Cavitation model parameters Constants in Hartman functions Specific heat at constant pressure Turbulence model parameters Defficiency Expected value of random variable x Average of surrogate models Error associated with i~th SUTTOgate model Approximation error at design point x Bias error at design point x Bias error bound at design point x Dataindependent bias error bound at design point x Root mean square bias error at design point x A a b b c` C Cdest, >prod CP Der E(x) Eavg E, e(x) eb (x) ebeb X) e (x) e "(x) ees(x) f(x) h h(x) h I K k L viz vi N DOE Ne Nlhs NRBF N, Ns N Ntest Standard error at design point x Function of design variables and decomposed functions Vector of basis functions in polynomial response surface model Vapor mass fraction Sensible enthalpy Radial basis function Vector of radial basis functions Identity matrix Thermal conductivity Turbulent kinetic energy Loss function, Latent heat Moment matrix Cavitation source terms Number of design of experiments Number of eigenvectors Number of Latin hypercube samples Number of radial basis functions Number of symbols for orthogonal arrays Number of sampled data points Number of surrogate models Number of test points N,, Np Np P p . p R R ad r r(x) sresp T T ~/ t u,v,w,ui,u ,uk,U Number of variables Number of basis functions in approximation model Number of basis functions missing from the approximation model Lz norm of the difference between predicted and benchmark experimental pressure data Turbulence production term Location of midpoint on the lower side of diffuser vane Pressure Constants in Hartman functions Correlation matrix Adjusted coefficient of multiple determination Strength of orthogonal array Radius of largest unoccupied sphere Vector of correlation between prediction and data points (kriging) Sensitivity indices, main, interaction, and total effects Standard deviation of responses Temperature Lz norm of the difference between predicted and benchmark experimental temperature data Time, magnitude of tangents of Bezier curves Velocity components V Null eigenvector V, V, V; Variance and its components (partial variances), volume w; Weight associated with ith SUTTOgate model X Gramian design matrix x Vector of design variables xi, x ,xk Space variables (coordinates) y Vector of responses y, y(x) Function or response value at a design point x 9, f(x) Predicted response at a design point x Z(x) Systematic departure term in kriging Z Subset of design variables for global sensitivity analysis a Volume fraction, Thermal diffusivity, Parameter in weighted average model a ,a2 Parameters used to define vertex and axial point locations in a central composite design P Vector of coefficients of true basis functions P* Estimated coefficient vector in kriging /7 Coefficient associated with a basis function in polynomial response surface approximation, coefficient of thermal expansion, parameter in weighted average model 7 Constant used to estimate root mean square bias error 3 Kroneckerdelta E, E(x) Turbulent dissipation, error in surrogate model r(x) True function or response B Probability density function, variable defining diffuser vane shape 0 Vector of parameters in the Gaussian correlation function /Z Regularization parameter pu Mean of responses at sampled points, dynamic viscosity Weights used in numerical integration p Density 2 Variance of noise, estimated process variance in kriging ~a Adjusted root mean square error Degree of correlation, flow variable 1 Vector of ones Sub scripts c Cavity, candidate true polynomial dopt Doptimal design i, j, kIndices krg Kriging I Liquid phase lhs Latin hypercube sampling design m Mixture of liquid and vapor max Maximum of the quantity min Minimum of the quantity prs Polynomial response surface pws PRESSbased weighted average surrogate rbnn Radial basis neural network RM~S Root mean square value t Turbulent wta Weighted average surrogate E Solution with least deviation from the data v Vapor phase oo Reference conditions Superscripts/Overhead Symbols I Data independent error measure R Reynolds stress (i) ith design point (i) All points except ith design point (1) Lower bound T Transpose (u) Upper bound (1) Terms in the response surface model (2) Terms missing from the response surface model ^ Predicted value Average of responses in surrogate models Nondimensional Numbers CFL Pr Re Courant, Freidricks and Levy number Prandtl number Reynolds number Cavitation number Acronyms ADS AIAA ASME ASP CCD CFD COV CV DOE ED EGO ESE FCCD GCV GMSE GSA KRG IGV LHS LH2 LOX All possible data sets American Institute of Aeronautics and Astronautics American Society of Mechanical Engineers Alkaline surfactantp olymer Central composite design Computational fluid dynamics Coefficient of variation Cross validation Design of experiment Experimental design Efficient global optimization Estimated standard error Facecentered central composite design Generalized cross validation error Generalized mean squared error Global sensitivity analysis Kriging Inlet guide vane Latin hypercube sampling Liquid Hydrogen Liquid Oxygen MSE Mean squared error NASA National Aeronautics and Space Administration NS NavierStokes NIST National Institute of Standards and Technology NPSF Nonparametric surrogate filter OA Orthogonal array POF Pareto optimal front PRESS Predicted residual sum of squares PRS Polynomial response surface PSF Parametric surrogate filter PWS PRESSbased weighted average surrogate RBF Radial basis function RBNN Radial basis neural network RMS Root mean square RMSBE Root mean square bias error RMSE Root mean squared error RP1 Refined petroleum SoS Speed of sound SQP Sequential quadratic programming SS Split sample SVR Support vector regression Operators E(a) Expected value of the quantity a max(a, b) Maximum of a and b min(a, b) Minimum of a and b r(a, b) Correlation between vectors a and b V(a) Variance of the quantity a cr(a) Standard deviation of the quantity a (a~iavg Spaceaveraged value of the quantity a (a>ma Maximum of the quantity a a Absolute value of the quantity a a (L2) Norm of the vector a 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 MULTIPLE SURROGATES AND ERROR MODELING IN OPTIMIZATION OF LIQUID ROCKET PROPULSION COMPONENTS By Tushar Goel May 2007 Chair: Raphael T. Haftka Cochair: Wei Shyy Major: Mechanical Engineering Design of space propulsion components is extremely complex, expensive, and involves harsh environments. Coupling of computational fluid dynamics (CFD) and surrogate modeling to optimize performance of space propulsion components is becoming popular due to reduction in computational expense. However, there are uncertainties in predictions using this approach, like empiricism in computational models and surrogate model errors. We develop methods to estimate and to reduce such uncertainties. We demonstrate the need to obtain experimental designs using multiple criteria by showing that using a singlecriterion may lead to high errors. We propose using an ensemble of surrogates to reduce uncertainties in selecting the best surrogate and sampling strategy. We also develop an averaging technique for multiple surrogates that protects against poor surrogates and performed at par with best surrogate for many problems. We assess the accuracy of different error estimation models, including an error estimation model based on multiple surrogates, used to quantify prediction errors. While no single error model performs well for all problems, we show possible advantage of combining multiple error models. We apply these techniques to two problems relevant to space propulsion systems. First, we employ surrogatebased strategy to understand the role of empirical model parameters and uncertainties in material properties in a cryogenic cavitation model, and to calibrate the model. We also study the influence of thermal effects on predictions in cryogenic environment in detail. Second, we use surrogate models to improve the hydrodynamic performance of a diffuser by optimizing the shape of diffuser vanes. For both problems, we observed improvements using multiple surrogate models. While we have demonstrated the approach using space propulsion components, the proposed techniques can be applied to any largescale problem. CHAPTER 1 INTRODUCTION AND SCOPE Liquid rocket propulsion systems are the most popular form of space propulsion systems for high thrust and specific impulse applications as required for space applications (Humble et al., 1995, Chapter 5). Unlike propulsion systems used in aircraft, space propulsion systems carry both fuel and oxidizer with the vehicle. This poses additional requirements on the selection of suitable propellants and design of propulsion systems. Apart from high energy density, the choice of propellants is also affected by the ease of storage and handling, mass or volume of propellant, and nature of products of combustion. Typical propellants used for liquid space propulsion are refined petroleum (RPl) with liquid oxygen (LOX), hypergolic propellants (monomethyl hydrazine with nitrogen tetroxide) and cryogens (liquid hydrogen LH2 and LOX). Cryogens (LH2 and LOX) are most popular due to higher power/gallon ratio and specific thrust, lower weight of LH2, and cleaner combustion products (water). Despite difficulties in storage (due to the tendency of cryogens to return to gaseous state unless supercooled, the boiling point of LH2 and LOX at standard conditions is  423o F and 298o F, respectively) and safety considerations, the rewards of using cryogens as space propellants are significant (NASA online facts, 1991). Space Propulsion Systems A conceptual schematic of a typical bipropellant space propulsion system is shown in Figure 11. There are five major components of the propulsion system: fuel and oxidizer storage tanks, fuel and oxidizer pumps, gas turbine, combustion chamber, and nozzle. Based on the type of power cycle, different propulsion systems are classified as follows: *GasGenerator Cycle (Figure 12(A)) A small amount of fuel and oxidizer is fed to gas generator where the fuel is burnt at less than optimal ratio to keep the temperature in turbine low. The hot gases produced in gas generator drive the turbine to produce power required to run the fuel and oxidizer pumps. Thrust of the engine is regulated by controlling the amount of propellants through gas generator. This is an open cycle since the hot gas from turbines is either dumped overboard or sent into the main nozzle downstream. This configuration is useful for moderate power requirements but not good for the applications which require high power. * Staged Combustion Cycle (Figure 12(B)) This is a closed cycle system in which an enriched mixture of fuel and oxidizer is generated in preburner, and after passing through the turbine this vaporized mixture is fed to the main combustion chamber. No fuel or oxidizer is wasted in this cycle as complete combustion takes place in the main combustion chamber. This cycle is used for high power applications but the engine development cost is high and components of propulsion system are subj ected to harsh conditions. * Expander Cycle (Figure 12(C)) In this closed cycle, the main combustion chamber is cooled by passing liquid fuel that in turn gets vaporized by exchange of heat. This fuel vapor runs the turbine and is fed back to the main combustion chamber, where complete combustion takes place. The limitation on heat transfer to the fuel limits the power available in the turbine. This configuration is more suitable to small/mid size engines. * PressureFed Cycle (Figure 12(D)) This is the simplest configuration which does not require any pump or turbine. Instead, the fuel and the oxidizer are fed to the combustion chamber by the tank pressure. This cycle can be applied only to relatively lowchamberpressure applications because higher pressure requires bulky storage tanks. Design Requirements of Propulsions Systems The design of an optimal propulsion system requires efficient performance of all components. As can be seen from Figure 12, pumps, turbine, combustion chamber, and nozzle are integral part of almost all space propulsion systems. The system level goal of designing a space propulsion system is to obtain the highest possible thrust with the lowest weight (Griffin and French, 1991), but the requirements of individual components are also governed by corresponding operating conditions. Storage tanks are required to withstand high pressure in cryogenic environments while keeping the weight of tank low. Nozzles impart high velocities to high temperature combustion products to produce maximum possible thrust. The requirements for the design of the turbine and the combustion chamber are high efficiency, compact design, and ability to withstand high pressures and temperatures. Similarly, pumps are required to supply the propellants to the combustion chamber at a desired flow rate and inj section pressure. Maj or issues in the design of pumps include harsh environments, compact design requirements, and cavitation under cryogenic conditions. Each subsystem has numerous design options for example, number of stages, number and geometry of blades in turbines and pumps, different configurations and geometry of inj ectors and combustion chambers etc., that relate to the subsystem level requirements. System Identification and Optimization: Case Studies The design of a propulsion system is extremely complex due to the often conflicting requirements posed on individual components and interaction of different subsystems. Experimental design of propulsion systems is very expensive, time consuming, and involves harsh environments. With improvements in numerical algorithms and increase in computational power, the role of computational fluid dynamics (CFD) for design of complex systems with various levels of complexities has grown manyfold. Computer simulation of design problems has not only reduced the cost of developing the designs but also has reduced risks and design cycle time. The flexibility of trying alternative design options also has increased extensively compared to the experiments. With improvements in computer hardware and CFD algorithms, the complexity of the simulation model is increasing in an effort to capture physical phenomenon more accurately. Our current efforts in the design of liquid rocket propulsion systems are focused on two distinct classes of problems. First, we try to understand the roles of thermal effects and model parameters used for numerical modeling of cryogenic cavitation; and second, we aim to carry out shape optimization of diffuser vanes to maximize diffuser efficiency. While our interests in the design of diffuser are motivated by the ongoing interests in the Moon and Mars exploration endeavors, the study of cryogenic cavitation model validation and sensitivity analysis is relatively more generic, but very relevant to the design of liquid rocket propulsion systems. We discuss each problem in more detail as follows. Sensitivity Evaluation and Model Validation for a Cryogenic Cavitation Model Though a lot of improvements have been made over the years in the numerical modeling of fluid physics, areas like cavitation, both in normal as well as cryogenic environment, are still continuously developing. Complex physics, phase changes and resulting variations in material properties, interaction of convection, viscous effects and pressure, time dependence, multiple time scales, interaction between different phases of fluids and many fluids, temperature dependent material properties, and turbulence make these problems difficult (Utturkar, 2005). Numerous algorithms (Singhal et al. 1997, Merkle et al. 1998, Kunz et al. 2000, Senocak and Shyy 2004ab, Utturkar et al. 2005ab) have been developed to capture the complex behavior of cavitating flows, which have serious implications on the performance of the propulsion components . We study one cryogenic cavitation model in detail to assess the role of thermal boundary conditions, thermal environment, uncertainties in material properties, and empirical model parameters on the prediction of pressure and temperature field. Finally, we calibrate the cryogenic cavitation model parameters and validate the outcome using several benchmark data. Shape Optimization of Diffuser Vanes The diffuser is a critical component in liquid rocket turbomachinery. The high velocity fluid from turbopump is passed through the diffuser to partially convert kinetic energy of the fluid into pressure. The efficiency of a diffuser is determined by its ability to induce pressure recovery that is characterized by the ratio of outlet to inlet pressure. While exploring different concepts for diffuser design, Dorney et al. (2006a) observed that diffusers with vanes are more effective than vaneless diffusers. Consequently, we seek further improvements in diffuser efficiency through shape optimization of vanes. Since computational cost of simulations is high, we adopt a surrogate modelbased framework for optimization and sensitivity analysis that is briefly introduced in next section and discussed in detail in a following chapter. Together, the present effort demonstrates that the same technical framework can be equally adapted to treat hardware and computational modeling development, with information to allow one to inspect the overall design space characteristics, to facilitate quantitative tradeoff considerations between of multiple and competing obj ectives, and to rank order the importance of the various design variables. Surrogate Modeling High computational cost of simulations involved in performance evaluation of propulsion components makes direct coupling of optimization tools and simulations infeasible for most practical problems. To alleviate the problems associated with the optimization of such complex and computationally expensive components, surrogate models based on a limited amount of data have been frequently used. Surrogate models offer a lowcost alternative to evaluate a large number of designs and are amenable to the optimization process. Surrogate models also can be used to assess the trends in the design space as well as help identify the problems associated with numerical simulations. A sample of recent application of surrogate models for the design of spacepropulsion system is given as follows. Lepsch et al. (1995) used polynomial response surface approximations to minimize the emptyweight of the dual fuel vehicles by considering propulsion systems and vehicle design parameters. Rowell et al. (1999) and the references within discuss the application of regression techniques for the design of singlestagetoorbit vehicles. Madsen et al. (2000) used polynomial response surface models for the design of diffusers. Gupta et al. (2000) used response surface methodology to improve the lifttodrag ratio of artificially blunted leading edge spherical cone that is a representative geometry for reentry vehicle while constraining the heat transfer rates. Chung and Alonso (2000) estimated boom and drag for low boom supersonic business jet design problem via kriging. Shyy et al. (2001ab) employed global design optimization techniques for the design of rocket propulsion components; Papila et al. (2000, 2001, 2002), Papila 2001, approximated the design obj ectives for supersonic turbines using polynomial response surface approximations and neural networks; Vaidyanathan et al. (2000, 2004ab), Vaidyanathan (2004) modeled performance indicators of liquid rocket inj ectors using polynomial response surface approximations; Simpson et al. (2001b) used kriging and polynomial response surface approximations for the design of aerospike nozzles; Steffen (2002a) developed an optimal design of a scramj et inj ector to simultaneously improve efficiency and pressure loss characteristics using polynomial response surfaces. In a followup work, Steffen et al. (2002b) used polynomial response surface approximations to analyze the parameter space for the design of a combined cycle propulsion components, namely, mixedcompression inlet, hydrogen fueled scramjet combustor, a ductedrocket nozzle. Charania et al. (2002) used response surface methodology to reduce computational cost in systemlevel uncertainty assessment for reusable launch vehicle design. Huque and Jahingir (2002) used neural networks in optimization of integrated inlet/ej ector system of an axisymmetric rocketbased combined cycle engine. Qu et al. (2003) applied different polynomial response surface approximations for the structural design of hydrogen tanks; Keane (2003) employed response surfaces for the optimal design of wing shapes; Umakant et al. (2004) used kriging surrogate models to compute probability density functions for robust configuration design of airbreathing hypersonic vehicle. Baker et al. (2004ab) used response surface methodology for system level optimization of booster and ramj et combustor while considering multiple constraints. Levy et al. (2005) used support vector machines and neural networks to approximate obj ectives temperature field and pressure loss while searching Pareto optimal solutions for the design of a combustor. Issues with Surrogate Modeling The accuracy of surrogate models is an important factor for its effective application in design and optimization process. Some of the issues that influence the accuracy of surrogate models are: (1) number and location of sampled data points, (2) numerical simulations, and (3) choice of the surrogate model. The characterization of uncertainty in predictions via different error estimation measures is also useful for optimization algorithms such as, EGO (Jones et al., 1998) etc. These issues are widely discussed in literature (Li and Padula 2005, Queipo et al. 2005). While a detailed review of different issues and surrogate models is provided in next chapter, we briefly discuss the intent of current research in context of addressing the issues related to surrogate prediction accuracy. Sampling Strategies Typically the amount of data used for surrogate model construction is limited by the availability of computational resources. Many design of experiments (DOE) techniques are developed to sample the locations for conducting simulations such that the errors in approximation are reduced. However, these strategies mostly optimize a single criterion that caters to an assumed source of error in approximation. The dominant sources of errors in approximation are rarely known a priori for practical engineering problems. This renders the selection of appropriate DOE technique a very difficult task because an unsuitable DOE may lead to poor approximation. We demonstrate this issue of high errors in approximation due to a singlecriterion based DOE with the help of simple examples, and highlight the need to simultaneously consider multiple criteria. Type of Surrogate Model The influence of the choice of surrogate model on prediction accuracy is widely explored in literature. Many researchers have compared different surrogates for various problems and documented their recommendations. The general conclusion is that no single surrogate model may prove adequate for all problems, and the selection of a surrogate for a given problem is often influenced by the past experience. However, as we demonstrate in this work, the suitability of any surrogate model also depends on the sampling density and sampling scheme. Then, the selection of an appropriate surrogate model becomes more complicated exercise. Here, we present methods to exploit available information by simultaneous use of multiple surrogate models. Specifically, we use multiple surrogate models to assess the regions of high uncertainty in response predictions. Further we develop a weighted average surrogate model, which is demonstrated to be more robust than the individual surrogate models for a wide variety of problems and sampling schemes. Estimation of Errors in Surrogate Predictions Since surrogate models only approximate the response, there are errors in predictions. An estimate of prediction errors is beneficial in determining the sampling locations in many surrogatebased optimization methods such as EGO (Jones et al., 1998) and for adaptive sampling. The effectiveness of these methods depends on the accuracy of error estimation measures. Mostly, these error estimation measures are based on statistical assumptions. For example, prediction variance for polynomial response surface approximation is developed assuming that errors are exclusively due to noise that follows a normal distribution of zero mean and cr2 Variance independent of locations. When these assumptions are not satisfied, the accuracy of error estimation measures is questionable. We compare different error estimation measures using a variety of test problems and give recommendations. We also explore the idea of simultaneously using multiple error estimation measures. While we use the presented surrogatebased framework for design and optimization of diffuser design and cryogenic cavitation model, we employ analytical examples, which are primarily used to test optimization algorithms, to exhibit the key concepts relating to different issues with surrogate modeling. Scope of Current Research In short, the goal of present work is to develop methodologies for designing optimal propulsion systems while addressing issues related to numerical uncertainties in surrogate modeling. The scope of the current work can be summarized as follows: 1) To illustrate risks in using a single criterion based experimental design for approximation and the need to consider multiple criteria. 2) To explore the use of an ensemble of surrogates to help identify regions of high uncertainties in predictions and to possibly provide a robust prediction method. 3) To appraise different error estimation measures and to present methods to enhance the error detection capability by combining multiple error measures. 4) To demonstrate the proposed surrogatemodel based approach to liquid rocket propulsion problems dealing with (a) cryogeniccavitation model validation and sensitivity study to appraise the influence of different parameters on performance, and (b) shape optimization of diffuser vanes to maximize diffuser efficiency. The organization of this work is as follows. We review different surrogate models and relevant issues associated with surrogate modeling in Chapter 2. In Chapter 3, we evidence risks in using a single criterion for constructing design of experiments. Methods to harness the potential of multiple surrogate models are presented in Chapter 4. The performance of different error estimation measures is compared in Chapter 5 and we propose methods to simultaneously use multiple error measures. This surrogatebased analysis and optimization framework is applied to two problems related to liquid rocket propulsion, model validation and sensitivity analysis of a cryogenic cavitation model in Chapter 6, and shape optimization of diffuser vanes in Chapter 7. We recapitulate maj or conclusions of the current work and delineate the scope of future work in Chapter 8. Oxid izer Figure 11. Schematic of liquid fuel rocket propulsion system. A man. uno B GARpap Om~Lhm RIUp~ C 5=* "" D Figure 12. Classification of propulsion systems according to power cycles. A) Gasgenerator cycle. B) Staged combustion cycle. C) Expander cycle. D) Pressurefed cycle. CHAPTER 2 ELEMENTS OF SURROGATE MODELING Surrogate models are widely accepted for the design and optimization of components with high computational or experimental cost (typically encountered in CFD simulations based designs), as they offer a computationally less expensive way of evaluating designs. Surrogate models are constructed using the limited data generated from the analysis of carefully selected designs. Numerous successful applications of surrogate models for design and optimization of aerospace systems, automotive components, electromagnetic applications and chemical processes etc. are available in literature. A few examples are given as follows. Kaufman et al. (1996), Balabanov et al. (1996, 1999), Papila and Haftka (1999, 2000), and Hosder et al. (2001) constructed polynomial response surface approximations (PRS) for structural weight based on structural optimizations of high speed civil transport. Hill and Olson (2004) applied PRS to approximate noise models in their effort to reduce the noise in the conceptual design of transport aircraft. Madsen et al. (2000), Papila et al. (2000, 2001), Shyy et al. (2001ab), Vaidyanathan et al. (2000, 2004ab), Goel et al. (2004), and Mack et al. (2005b, 2006) used polynomial and neural networksbased surrogate models as design evaluators for the optimization of propulsion components including turbulent flow diffuser, supersonic turbine, swirl coaxial inj ector element, liquid rocket inj ector, and radial turbine designs. Burman et al. (2002), Goel et al. (2005), and Mack et al (2005a) used different surrogates to maximize the mixing efficiency facilitated by a trapezoidalshaped bluff body in the time dependent Navier Stokes flow while minimizing the resulting drag coefficient. Knill et al. (1999), Rai and Madavan (2000, 2001) and, Madavan et al. (2001) used surrogate models for airfoil shape optimization. Dornberger et al. (2000) used neural networks and polynomial response surfaces for design of turbine blades. Kim et al. (2002), Keane (2003), Jun et al. (2004) applied PRS to optimize wing designs. Ong et al. (2003) used radial basis functions to approximate the obj ective function and constraints of an aircraft wing design. Bramantia et al. (2001) used neural networkbased models to approximate the design obj ectives in electromagnetic problems. Farina et al. (2001) used multiquadrics interpolation based response surface approximations to optimize the shape of electromagnetic components like, Ccore and magnetizer. Wilson et al. (2001) used response surface approximations and kriging for approximating the obj ectives while designing piezomorph actuators. Redhe et al. (2002ab), Craig et al. (2002), and Stander et al. (2004) used PRS, kriging and neural networks in design of vehicles for crashworthiness. RaisRohani and Singh (2002) used PRS to approximate limit state functions for estimating the reliability of composite structures. Rikards and Auzins (2002) developed PRS to model buckling and axial stiffness constraints while minimizing the weight of composite stiffened panels. Vittal and Haj ela (2002) proposed using PRS to estimate the statistical confidence intervals on the reliability estimates. Zerpa et al. (2005) used kriging, radial basis functions, and PRS to optimize the cumulative oil recovery from a heterogene ou s, multi pha se re servoi r subj ect to an A SP (al kal ine surfactantp olymer) flooding . Steps in Surrogate Modeling Li and Padula (2005) and Queipo et al. (2005) have given a comprehensive review of the relevant issues in surrogate modeling. In this chapter, we discuss the key steps in the surrogate modeling as explained with the help of Figure 21. Discussion of maj or issues and the most prominent approaches followed in each step of the surrogate modeling process, briefly described below, lays the outline of this chapter. Design of Experiments (DOEs) The design of experiment is the sampling plan in design variable space. Other common names for DOEs are experimental designs or sampling strategies. The key question in this step is how we assess the goodness of such designs, considering the number of samples is severely limited by the computational expense of each sample. We discuss the most prominent approaches related to DOE in a subsequent section. Later in Chapter 3, we demonstrate some practical issues with the construction of DOEs. Numerical Simulations at Selected Locations Here, the computationally expensive model is executed for all the designs selected using the DOE specified in the previous step. In context of the present work, the details of relevant numerical simulation tools used to evaluate designs were briefly discussed in appropriate chapters . Construction of Surrogate Model Two questions are of interest in this step: 1) what surrogate models) should we use (model selection) and, 2) how do we Eind the corresponding parameters (model identification)? A formal description of the problem of interest is discussed in next section. A framework for the discussion and mathematical formulation of alternative surrogatebased modeling approaches is outlined in next section and the section on construction of surrogate models. Model Validation The purpose of this step is to establish the predictive capabilities of the surrogate model away from the available data (generalization error). Different schemes to estimate generalization error for model validation are discussed in this chapter. We compare different error estimation measures specific to a few popular surrogates in a following chapter. Mathematical Formulation of Surrogate Modeling Problem With reference to Figure 22, surrogate modeling can be seen as a nonlinear inverse problem for which one aims to determine a continuous function ( y(x)) of a set of design variables from a limited amount of available data (y). The available data y while deterministic in nature can represent exact evaluations of the function y(x) or noisy observations; and in general cannot carry sufficient information to uniquely identify y(x). Thus, surrogate modeling deals with the twin problems of: 1) constructing a model f(x) from the available data y (model estimation), and 2) assessing the errors E attached to it (model appraisal). A general description of the anatomy of inverse problems can be found in Snieder (1998). Using the surrogate modeling approach the prediction of the simulationbased model output is formulated as y(x) = f(x)+ e(x) The prediction expected value and its variance V(y) are illustrated in Figure 23, with B being a probability density function. Different model estimation and model appraisal components of the prediction have been shown to be effective in the context of surrogate based analysis and optimization (see for example, McDonald et al., 2000; Chung and Alonso, 2000; Simpson et al., 2001a; Jin et al., 2001), namely polynomial response surface approximation (PRS), Gaussian radial basis functions (GRF) (also referred as radial basis neural networks RBNN), and (ordinary) kriging (KRG) as described by Sacks et al. (1989). Model estimation and appraisal components of these methods are presented in a following section. A good paradigm to illustrate how particular solutions (9) to the model estimation problem can be obtained is provided by regularization theory (see for example, Tikhonov and Arsenin (1977), and Morozov (1984)), which imposes additional constraints on the estimation. More precisely, jl can be selected as the solution to the following Tikhonov regularization problem: minI(9),=h ~Lv y, (x"))+AnD"'@sdx (2.1) where S is the family of surrogate models under consideration, L(x) is a loss or cost function used to quantify the so called empirical error (e.g., L(x) = x ), Ai is a regularization parameter, and D"'f represents the value of the mderivative of the proposed model at location x. Note that D"'f represents a penalty term; for example, if m is equal to two, it penalizes high local curvature. Hence, the first term enforces closeness to the data (goodness of fit), while the second term addresses the smoothness of the solution with Ai (a real positive number) establishing the tradeoff between the two. Increasing values of Ai provide smoother solutions. The purpose of the regularization parameter Ai is, hence, to help implement Occam's razor principle (Ariew, 1976), which favors parsimony or simplicity in model construction. A good discussion on statistical regularization of inverse problems can be found in Tenorio (2001). The quadratic loss function (i.e., Lz norm) is most commonly used in part because it typically allows easy estimation of the parameters associated with the surrogate model; however, it is very sensitive to outliers. The linear (also called Laplace) loss function takes the absolute value of its argument (i.e., L1 norm); on the other hand, the Huber loss function is defined as quadratic for small values of its argument and linear otherwise. The so called r loss function has received considerable attention in the context of the support vector regression surrogate (Vapnik, 1998; Girosi, 1998), and assigns an error equal to zero if the true and estimated values are within an r distance. Figure 24 illustrates the cited loss functions. Design of Experiments As stated earlier, the design of experiment is the sampling plan in design variable space and the key question in this step is how we assess the goodness of such designs. In this context, of particular interest are sampling plans that provide a unique value (in contrast to random values) for the input variables at each point in the input space, and are modelindependent; that is, they can be efficiently used for fitting a variety of models. Typically, the primary interest in surrogate modeling is minimizing the error and a DOE is selected accordingly. Two maj or components of the empirical error and corresponding expressions for average error formulation considering quadratic loss functions are given as follows. * Variance: measures the extent to which the surrogate model f(x) is sensitive to particular data set D. Each data set D corresponds to a random sample of the function of interest. This characterizes the noise in the data, Evar (x) = E4DS [f E4DS [~X)1]2 (2.2) * Bias: quantifies the extent to which the surrogate model outputs (i.e., f(x)) differ from the true values (i.e., y(x)) calculated as an average over all possible data sets D (ADS), ELm,, (x) = CE4D iSX)] y(x) (2.3) In both expressions, E4DS denotes the expected value considering all possible data sets. There is a natural tradeoff between bias and variance. A surrogate model that closely fits a particular data set (lower bias) will tend to provide a larger variance and viceversa. We can decrease the variance by smoothing the surrogate model but if the idea is taken too far then the bias error becomes significantly higher. In principle, we can reduce both bias (can choose more complex models) and variance (each model more heavily constrained by the data) by increasing the number of points, provided the latter increases more rapidly than the model complexity. In practice, the number of points in the data set is severely limited (e.g., due to computational expense) and often during the construction of surrogate model, a balance between bias and variance errors is sought. This balance can be achieved, for example, by reducing the bias error while imposing penalties on the model complexity (e.g., Tikhonov regularization). With reference to most applications, where the actual model is unknown (see previous and next sections), and data is collected from deterministic computer simulations, bias error is the dominant source of error because the numerical noise is small, and a DOE is selected accordingly. When response surface approximation is used, there is good theory for obtaining minimum bias designs (Myers and Montgomery, 1995, Section 9.2) as well as some implementations in low dimensional spaces (Qu et al., 2004). For the more general case, the bias error can be reduced through a DOE that distributes the sample points uniformly in the design space (Box and Draper, 1959; Sacks and Ylvisaker, 1984; as referenced in Tang, 1993). However, fractional factorial designs replace dense full factorial designs for computationally expensive problems. The uniformity property in designs is sought by, for example, maximizing the minimum distances among design points (Johnson et al., 1990), or by minimizing correlation measures among the sample data (Iman and Conover, 1982; Owen, 1994). Practical implementations of a few most commonly used DOEs are discussed next. Factorial Designs Factorial designs are one of the simplest DOEs to investigate the effect of main effects and interactions of variables on the response for boxshaped design domain. 2N factorial design where each design variable is allowed to take two extreme levels is often used as a screening DOE to eliminate the unimportant variables. Qualitative and binary variables can also be used for this DOE. A typical twolevel full factorial DOE for three variables is shown in Figure 25. These designs can be used to create a linear polynomial response surface approximation. For higher order approximations, the number of levels of each variable is increased for example a quadratic polynomial response surface approximation can be fitted by a threelevel full factorial design (3N' designs). Some times the number of experiments is reduced by using 2hi (or 1/2p) fractional factorial designs, which require the selection of p independent design generators (least influential interactions). Typically, these designs are classified according to resolution number: * Resolution III: Main effects are not aliased with other main effects, but are confounded with one or more twoway interactions. * Resolution IV: Main effects are not aliased with other main effects or twoway interactions. Two factor interactions are confounded with other twoway interactions. * Resolution V: Main effects and twoway interactions are not confounded with one another. More details can be obtained from Myers and Montgomery (1995, Chapter 4, pp. 156 179). Factorial designs produce orthogonal designs for polynomial response surface approximation. However, for higher dimensions (N~ > 6 ), factorial designs requires a large number of experiments making them particularly unattractive for computationally expensive problems. Central Composite Designs These include designs on 2 '~ vertices, 2N~ axial points, and N~ repetitions of central point. The distance of axial point a is varied to generate, facecentered, spherical design, or orthogonal design. A typical central composite design for threevariable problem is shown in Figure 26. These designs reduce variance component of the error. The repetitions at the center reduce the variance, improve stability (defined as the ratio of maximum variance to minimum variance in the entire design space), and give an idea of magnitude of the noise, but are not useful for the computer simulations. These designs are also not practical for higher dimension spaces (N, > 8) as the number of simulations becomes very high. For N, > 3 when the designs on the vertices of the design spaces are not feasible, BoxBehnken designs can be used for quadratic polynomial response surface approximation. These are spherical designs (sampling designs at these locations enables us to exactly determine the function value on the points equi distant from the center) but these designs introduce higher uncertainty near the vertices. Variance Optimal DOEs for Polynomial Response Surface Approximations Moment matrix M = XIX/N 1 for PRS (see next section) is very important quantity as this affects the prediction variance, and the confidence on the coefficients, hence, used to develop different variance optimal DOEs. Doptimal design maximizes the determinant of the moment matrix M~to maximize the confidence in the coefficients of polynomial response surface approximation. Aoptimal design minimizes the trace of inverse of2~to minimize the sum of variances of the coefficients. Goptimal design minimizes the maximum of prediction variance. Ioptimal design minimizes the integral of the prediction variance over the design domain. All these DOEs require the solution of difficult optimization problem, which is solved heuristically in higher dimensional spaces. Latin Hypercube Sampling (LHS) Stratified sampling ensures that all portions of a given partition are sampled. LHS (McKay et al., 1979) is a stratified sampling approach with the restriction that each of the input variables (xk) has all portions of its distribution represented by input values. A sample of size Ns can be constructed by dividing the range of each input variable into N, strata of equal marginal SAn example of matrix X is given in Equation (2.6). probability 1/Ns and sampling once from each stratum. Let us denote this sample by x j = 1, 2,..., Ns; k = 1, 2,.., N .The sample is made of components of each of the x 's matched at random. Figure 27 illustrates a LHS design for two variables, when six designs are selected. While LHS represents an improvement over unrestricted stratified sampling (Stein, 1987), it can provide sampling plans with very different performance in terms of uniformity measured by, for example, maximum minimumdistance among design points, or by correlation among the sample data. Figure 28 illustrates this shortcoming; the LHS plan in Figure 28(B) is significantly better than that in Figure 28(A). Orthogonal Arrays (OA) These arrays were introduced by C. R. Rao in the late 40's (Rao, 1947), and can be defined as follows. An OA of strength r is a matrix of N, rows and N, columns with elements taken from a set of N, symbols, such that in any Ns x r submatrices each of the (N, ), possible rows occur the same number 7 (index) of times. The number of rows ( Ns) and columns ( N, ) in the OA definition represents the number of samples and variables or factors under consideration, respectively. The N, symbols are related to the levels defined for the variables of interest, and the strength r is an indication of how many effects can be accounted for (to be discussed later in this section) with values typically between two and four for reallife applications. Such an array is denoted by OA( Ns, N,, N r ). Note that, by definition, a LHS is an OA of strength one, OA (Ns, N,, N ,1). There are two limitations on the use of OA for DOE: *Lack of flexibility: Given desired values for the number of rows, columns, levels, and strength, the OA may not exist. For a list of available orthogonal arrays, theory and applications, see, for example, Owen (1992), Hedayat et al. (1999), and references therein. *Point replicates: OA designs proj ected onto the subspace spanned by the effective factors (most influential design variables) can result in replication of points. This is undesirable for deterministic computer experiments where the bias of the proposed model is the main concern. Optimal LHS, OAbased LHS, Optimal OAbased LHS Different approaches have been proposed to overcome the potential lack of uniformity of LHS. Among those, most of them adjust the original LHS by optimizing a spread measure (e.g., minimum distance or correlation) of the sample points. The resulting designs have been shown to be relatively insensitive to the optimal design criterion (Ye et al., 2000). Examples of this strategy can be found in the works of Iman and Conover (1982), Johnson et al. (1990), and Owen (1994). Tang (1993) and Ye (1998) presented the construction of strength r OAbased LHS which stratify each rdimensional margin, and showed that they offer a substantial improvement over standard LHS. Another strategy optimizes a spread measure of the sample points, but restricts the search of LHS designs, which are orthogonal arrays, resulting in so called optimal OAbased LHS (Leary et al., 2003). Adaptive DOE, in which model appraisal information is used to place additional samples, have also been proposed (Jones et al., 1998, Sasena et al., 2000, Williams et al., 2000). A summary of main characteristics and limitations of different DOE techniques is listed in Table 21. If feasible, two sets of DOE are generated, one (so called training data set) for the construction of the surrogate (next section), and second for assessing its quality (validation as discussed in a later section). Given the choice of surrogate, the DOE can be optimized to suit a particular surrogate. This has been done extensively for minimizing variance in polynomial RSA (e.g., D and A optimal designs, Myers and Montgomery, 1995, Chapter 8) and to some extent for minimizing bias (e.g., Qu et al., 2004). However, for nonpolynomial models, the cost of the optimization of a surrogatespecific DOE is usually prohibitive, and so is rarely attempted. Construction of Surrogate Model There are both parametric (e.g., polynomial response surface approximation, kriging) and nonparametric (e.g., proj ectionpursuit regression, radial basis functions) alternatives for constructing surrogate models. The parametric approaches presume the global functional form of the relationship between the response variable and the design variables is known, while the non parametric ones use different types of simple, local models in different regions of the data to build up an overall model. This section discusses the estimation and appraisal components of the prediction of a sample of both parametric and nonparametric approaches. Specifically, the model estimation and appraisal components corresponding to polynomial response surface approximation (PRS), kriging (KRG), and radial basis functions (RBF) surrogate models are discussed next, followed by a discussion of a more general nonparametric approach called kernelba~sed regression. Throughout this section a square loss function is assumed unless otherwise specified, and given the stochastic nature of the surrogates, the available data is considered a sample of a population. Polynomial Response Surface Approximation (PRS) The regression analysis is a methodology to study the quantitative relation between a function of interest y and Npl basis functions f,, using Ns sample values of the response y, a for a set of basis functions ~f z) (Draper and Smith, 1998, Section 5.1). Monomials are the most preferred basis functions by practitioners. For each observation i, a linear equation is formulated: y, (f (x)) = ff '/~) + EI; E(ez) = 0, V'(ez) = cr2, (2.4) where the errors E, are considered independents with expected value equal to zero and variance 02 ,and 13 represents the quantitative relation between basis functions. The set of equations specified in Equation (2.4) can be expressed in matrix form as: y = XI) + E; E(E) = 0, V(E) = a 2I, (2.5) where Xis a Ns x Np matrix of basis functions, also known as Gramian design matrix, with the design variable values for sampled points. A Gramian design matrix for a quadratic polynomial in two variables (N, = 2; Npl = 6) is shown in Equation (2.6) 1xx ') x,(1)2 X1(1) (1) (2) 1 2) 1 (2) (2) 1(2)2 X1(2) (2) (2) X = 1 ,22 12 2 22(2.6) 1 x" "1 x(N,)2 1(N) (Ns) (Ns)2 J(x) = [ b f, (x), (2.7) where bj is the estimated value of the coefficient associated with the jth basis function fi (x). Then, the error in approximation at a design point x is given as e(x) = y(x) (x) The coefficient vector b can be obtained by minimizing a loss function L, defined as L = [ez ,(2.8) I=1 where el is the error at ith data point, p is the order of loss function, and N, is the number of sampled design points. A quadratic loss function (p = 2), that minimizes the variance of the error in approximation, is mostly used because we can estimate coefficient vector b using an analytical expression, as follows, b =(XTX) 1X y. (2.9) The estimated parameters b (by least squares) are unbiased estimates of 3 that minimize variance. At a new set of basis function vector f for design point P, the predicted response and the variance of the estimation are given as f(x)= bl f,(x;), an~d V (f(x))= cr2f T(XI'X)'f). (2.10) Kriging Modeling (KRG) Kriging is named after the pioneering work of D.G. Krige (a South African mining engineer), and was formally developed by Matheron (1963). More recently, Sacks et al. (1989, 1993), and Jones et al. (1998) made it wellknown in the context of the modeling, and optimization of deterministic functions, respectively. The kriging method in its basic formulation estimates the value of a function (response) at some unsampled location as the sum of two components: the linear model (e.g., polynomial trend), and a systematic departure representing low (large scale) and high frequency (small scale) variation components, respectively. The systematic departure component represents the fluctuations around the trend, with the basic assumption being that these are correlated, and the correlation is a function of distance between the locations under consideration. More precisely, it is represented by a zero mean, secondorder, stationary process (mean and variance constant with a correlation depending on a distance) as described by a correlation model. Hence, these models (ordinary kriging) suggest estimating deterministic functions as: y(x) = pu + E(x), E(E) = 0, cov(E(x(I), E(x(")) f 0 Vi, j, (2.11) where pu is the mean of the response at sampled design points, and E is the error with zero expected value, and with a correlation structure that is a function of a generalized distance between the sample data points. A possible correlation structure (Sacks et al., 1989) is given by: cvex) x ) = a2 O 1 f k x ) (2.12) where N_ denotes the number of dimensions in the set of design variables x; a identifies the standard deviation of the response at sampled design points, and, A is a parameter, which is a measure of the degree of correlation among the data along the kth direction. Specifically, the parameters pu, E, and # are estimated using a set of N, samples (x, y), such that a likelihood function is maximized (Sacks et al., 1989). Given a probability distribution and the corresponding parameters, the likelihood function is a measure of the probability of the sample data being drawn from it. The model estimates at unsampled points is: E~fx)) p r"R(f p),(2.13) where ^` above the letters denotes estimates, r identifies the correlation vector between the set of prediction points and the points used to construct the model, R is the correlation matrix among the Ns sample points, and 1 denotes an Ns vector of ones. On the other hand, the estimation variance at unsampled design points is given by: 1R (2r)4 V(y(x))= a2 rTR'r+ (.4 1T R 1 Gaussian processes (Williams and Rasmussen, 1996, Gibbs, 1997), another wellknown approach to surrogate modeling, can be shown to provide identical expressions for the prediction and prediction variance to those provided by kriging, under the stronger assumption that the available data (model responses) is a sample of a multivariate normal distribution (Rao, 2002, Section 4a). Radial Basis Functions (RBF) Radial basis functions have been developed for the interpolation of scattered multivariate data. The method uses linear combinations of NRB radially symmetric functions hr (x), based on Euclidean distance or other such metric, to approximate response functions as, NRBF y(x)= w~lh,(x)+4, (2.15) where w represents the coefficients of the linear combinations, h~ (x) the radial basis functions, and E, independent errors with variance cr2 The flexibility of the model, that is its ability to fit many different functions, derives from the freedom to choose different values for the weights. Radial basis functions are a special class of functions with their main feature being that their response decreases (or increases) monotonically with distance from a central point. The center, the distance scale, and the precise shape of the radial function are parameters of the model. A typical radial function is the Gaussian, which in the case of a scalar input, is expressed as, h, (x)= exp .) (2.16) The parameters are its center c and its radius 3. Note that the response of the Gaussian RBF decreases monotonically with the distance from the center, giving a significant response only in the center neighborhood. Given a set of N, input/output pairs (sample data), a radial basis function model can be expressed in matrix form as, f = Hw, (2.17) where H is a matrix given by, hdx (X) hdx (g  h xB (l)) H = z (2.18) Similar to polynomial response surface approximation method, by solving Equation (2. 17), the optimal weights (in the least squares sense) can be found to be, w = A H'y, (2.19) where A' is a matrix given by, A = (H H) (2.20) The error variance estimate can be shown to be given by, 2 yTr2y cr (2.21) tr ( .) where P, is a proj section matrix, , = I HA H (2.22) The RBF model estimate for a new set input values is given by, J(x) = hl v~, (2.23) where h is a column vector with the radial basis functions evaluations, h, (x) h = h (x)(2.24) On the other hand, the prediction variance is the variance of the estimated model f(x) plus the error variance and is given by: PVy(x))= V(h~ri)+ V(g)= (h7(HrH) h+1) y Pr(225 Ns NRB Radial basis function is also known as radial basis neural networks (RBNN) as described by Orr (1996, 1999ab). The MATLAB implementation of radial basis functions or RBNN (function 'newrb'), used in this study, is described as follows. Radial basis neural networks are twolayer networks consisting of a radialbasis function and a linear output layer. The output of each neuron is given by f =radbujas wxx0.832 prad (2.26) radba~s(x) = exp(x2), (2.27) where w is the weight vector associated with each neuron, x is the input design vector, 'spread' is a user defined value that controls the radius of influence for each neuron where the radius is half the value of parameter 'spread'. Specifically, the radius of influence is the distance at which the influence reaches a certain small value. If 'spread' is too small, the prediction will be poor in regions that are not near the position of a neuron; and if 'spread' is too large, the sensitivity of the neurons will be small. Neurons are added to the network one by one until a specified mean square error goal 'goal' is reached. If the error goal is set to zero, neurons will be added until the network exactly predicts the input data. However, this can lead to overfitting of the data, which may result in poor prediction between data points. On the other hand, if error goal is large, the network will be undertrained and predictions even on data points will be poor. For this reason, an error goal is judiciously selected to prevent overfitting while keeping the overall prediction accuracy high. Kernelbased Regression The basic idea of RBF can be generalized to consider alternative loss functions and basis functions in a scheme known as kernelbased regression. With reference to Equation (2.1), it can be shown that independent of the form of the loss function L, the solution of the variational problem has the form (the Representer Theorem: see Girosi, 1998; Poggio and Smale, 2003): j)(x) = az G(x, x '') >+b, (2.28) where G(x, x ')) is a (symmetric) kernel function that determines the smoothness properties of the estimation scheme. Table 22 shows the kernel functions of selected estimation schemes with the kernel parameters being estimated by model selection approaches (see next section for details). If the loss function L is quadratic, the unknown coefficients in Equation (2.28) can be obtained by solving the linear system, (NsAI~+G)a = y, (2.29) where l is the identity matrix, and G a square positive definite matrix with elements Ggl = G;((xm, xij) Note that, the linear system i s well posed, since (NsAIll+ G) i s strictly positive and well conditioned for large values of NsA If loss function L is nonquadratic, the solution of the variational problem still has the form of Equation (2.28) but the coefficients a~ are found by solving a quadratic programming problem in what is known as 'support vector regression' (Vapnik, 1998). Maj or characteristics of different surrogate models are summarized in Table 23. Comparative studies have shown that depending on the problem under consideration, a particular modeling scheme (e.g., polynomial response surface approximation, kriging, radial basis functions) may outperform the others and in general, it is not known a priori which one should be selected. See for example, the works of Friedman and Stuetzle (1981), Yakowitz & Szidarovsky (1985), Laslett (1994), Giunta and Watson (1998), Simpson et al. (2001ab), Jin et al. (2001). Considering, plausible alternative surrogate models can reasonably fit the available data, and the cost of constructing surrogates is small compared to the cost of the simulations, using multiple surrogates may offer advantages compared to the use of a single surrogate. Recently, multiple surrogatebased analysis and optimization approaches have been suggested by Zerpa et al. (2005) and Goel et al. (2006b) based on the model averaging ideas of Perrone and Cooper (1993), and Bishop (1995). The multiple surrogatebased analysis approach is based on the use of weighted average models, which can be shown to reduce the prediction variance with respect to that of the individual surrogates. The idea of multiple surrogatebased approximations is discussed in Chapter 4. Model Selection and Validation Generalization error estimates assess the quality of the surrogates for prediction and can be used for model selection among alternative models; and establish the adequacy of surrogate models for use in analysis and optimization studies (validation). This section discusses the most prominent approaches in the context of surrogate modeling. Split Sample (SS) In this scheme, the sample data is divided into training and test sets. The former is used for constructing the surrogate while the latter, if properly selected, allows computing an unbiased estimate of the generalization error. Its main disadvantages are, that the generalization error estimate can exhibit a high variance (it may depend heavily on which points end up in the training and test sets), and that it limits the amount of data available for constructing the surrogates. Cross Validation (CV) It is an improvement on the split sample scheme that allows the use of the most, if not all, of the available data for constructing the surrogates. In general, the data is divided into k subsets (kfold crossvalidation) of approximately equal size. A surrogate model is constructed k times, each time leaving out one of the subsets from training, and using the omitted subset to compute the error measure of interest. The generalization error estimate is computed using the k error measures obtained (e.g., average). If k equals the sample size, this approach is called leaveone out crossvalidation (known also as PRESS in the polynomial response surface approximation terminology). Equation (2.30) represents a leaveoneout calculation when the generalization error is described by the mean square error (GMSE). GM~SE = k (7)2,(230 where g(') represents the prediction at x ') using the surrogate constructed from all sample points except (x ') y, ). Analytical expressions are available for leaveoneout case for the GMSE without actually performing the repeated construction of the surrogates for both polynomial response surface approximation (Myers and Montgomery, 1995, Section 2.7) and kriging (Martin and Simpson, 2005). The advantage of crossvalidation is that it provides nearly unbiased estimate of the generalization error, and the corresponding variance is reduced (when compared to splitsample) considering that every point gets to be in a test set once, and in a training set k1 times (regardless of how the data is divided); the variance of the estimation though may still be unacceptably high in particular for small data sets. The disadvantage is that it requires the construction of k surrogate models; this is alleviated by the increasing availability of surrogate modeling tools. A modified version of the CV approach called GCVgeneralized cross validation, which is invariant under orthogonal transformations of the data (unlike CV) is also available (Golub et al., 1979). If the Tikhonov regularization approach for regression is adopted, the regularization parameter ii can be identified using one or more of the following alternative approaches: CV cross validation (leaveoneout estimates), GCV (smoothed version of CV), or the Lcurve (explained below). While CV and GCV can be computed very efficiently (Wahba, 1983; Hutchinson and de Hoog, 1985), they may lead to very small values of ii even for large samples (e.g., very flat GCV function). The Lcurve (Hansen, 1992) is claimed to be more robust and have the same good properties of GCV. The Lcurve is a plot of the residual norm (first term) versus the norm nII)m jH of the solution for different values of the regularization parameter and displays the compromise in the minimization of these two quantities. The best regularization parameter is associated with a characteristic Lshaped "corner" of the graph. Bootstrapping This approach has been shown to work better than crossvalidation in many cases (Efron, 1983). In its simplest form, instead of splitting the data into subsets, subsamples of the data are considered. Each subsample is a random sample with replacement from the full sample, that is, it treats the data set as a population from which samples can be drawn. There are different variants of this approach (Hall, 1986; Efron and Tibshirani, 1993; Hesterberg et al., 2005) that can be used for model identification as well as for identifying confidence intervals for surrogate model outputs. However, this may require considering several dozens or even hundreds of subsamples. For example, in the case of polynomial response surface approximation (given a model), regression parameters can be estimated for each of the subsamples and a probability distribution (and then confidence intervals) for the parameters can be identified. Once the parameter distributions are estimated, confidence intervals on model outputs of interest (e.g., mean) can also be obtained. Bootstrapping has been shown to be effective in the context of neural network modeling; recently, its performance in the context of model identification in regression analysis is also being explored (Ohtani, 2000, Kleijnen and Deflandre 2004). I Design of Experiments  Numerical Simulations at Selected Locations Model Validation Figure 21. Key stages of the surrogatebased modeling approach. Simulation Analysis problem based model y Appraisal Data yi problem If necessary Construction of Surrogate Models (Model Selection and Identification) Estimation problem Figure 22. Anatomy of surrogate modeling: model estimation + model appraisal. The former provides an estimate of function while the latter forecasts the associated error. I I B V B(E(y), V(y)) E(y) Figure 23. A surrogate modeling scheme provides the expected value of the prediction E(y) (solid line) and the uncertainty associated with that prediction, illustrated here using a probability density function 8 .  C  D Figure 24. Alternative loss functions for the construction of surrogate models. A) Quadratic. B) Laplace. C) Huber. D) r loss function. Figure 25. A twolevel full factorial design of experiment for three variables. Figure 26. A central composite design for threedimensional design space. e Y. b b b b b Figure 27. A representative Latin hypercube sampling design with Ns = 6, N~ = 2 for uniformly distributed variables in the unit square. Figure 28. LHS designs with significant differences in terms of uniformity (Leary et al., 2003). A) Random LHS. B) Correlation optimized LHS. C) OAbased LHS. (Figure reprinted with kind permission of Taylor and Francis group, Leary et al., 2003, Figure Table 21. Summary of main characteristics of different DOEs. Limitations Table 22. Examples of kernel functions and related estimation schemes. Kernel function Estimation scheme Gr(x, x,)= (1 + xxl)d Polynomial of degree d (PRD) G~xx)=xX x, Linear splines (LSP) G~x, xl ex i x 2 Gaussian radial basis function (GRF) Main features Used to investigate main effect and interaction of variables for boxshaped domains, gives orthogonal designs, caters to noise Cater to noise, applicable for boxshaped domains, repetition of points on center to improve stability Cater to noise, Applicable to irregular domains too  Maximize confidence in coefficients  Minimize the sum of variances of coefficients  Minimize the maximum of prediction variance  Minimize the integral of prediction variance over design domain Caters to bias error, stratified sampling method, good for high number of variables Name Factorial designs CCD Variance optimal designs Doptimal A optimal Goptimal Ioptimal Latin hypercube sampling (LHS) Orthogonal arrays (OA) OAbased LHS Irregular domains not good for N, > 6 Irregular domains, not good for N, > 8, repetition of points not useful for simulations High computational cost, not good when noise is low Not good when noise is significant, occasional poor DOE due to random components Limited number of orthogonal arrays, difficult to create OAs Limited OAs, may leave large holes in design space Boxshaped domains, the moment matrix is diagonal for monomial basis functions so coefficients of approximation are uncorrelated Combine OA and LHS designs to improve distribution of points Table 23. Summary of main characteristics of different surrogate models. Name Main features Polynomial response Global parametric approach, good for slow varying functions, easy to surface approximation construct, good to handle noise, not very good for simulations based data Kriging Global parametric approach, handles smooth and fast varying functions, computationally expensive for large amount of data Radial basis function Local, nonparametric approach, computationally expensive, good for fast varying functions Kernelbased Global, nonparametric approach, uses different loss functions, functions relatively new approach CHAPTER 3 PITFALLS OF USING A SINGLE CRITERION FOR SELECTING EXPERIMENTAL DESIGNS Introduction Polynomial response surface (PRS) approximations are widely adopted for solving optimization problems with high computational or experimental cost as they offer a computationally less expensive way of evaluating designs. It is important to ensure the accuracy of PRSs before using them for design and optimization. The accuracy of a PRS, constructed using a limited number of simulations, is primarily affected by two factors: (1) noise in the data; and (2) inadequacy of the fitting model (called modeling error or bias error). In experiments, noise may appear due to measurement errors and other experimental errors. Numerical noise in computer simulations is usually small, but it can be high for illconditioned problems, or if there are some unconverged solutions such as those encountered in computational fluid dynamics or structural optimization. The true model representing the data is rarely known, and due to limited data available, usually a simple model is fitted to the data. For simulationbased PRS, modeling/bias error due to an inadequate model is mainly responsible for the error in prediction. In design of experiments techniques, sampling of the points in design space seeks to reduce the effect of noise and reduce bias errors simultaneously. However, these obj ectives (noise and bias error) often conflict. For example, noise rejection criteria, such as Doptimality, usually produce designs with more points near the boundary, whereas the bias error criteria tend to distribute points more evenly in design space. Thus, the problem of selecting an experimental design (also commonly known as design of experiment or DOE) is a multiobj ective problem with conflicting obj ectives (noise and bias error). The solution to this problem would be a Pareto optimal front of experimental designs that yields different tradeoffs between noise and bias errors. Seeking the optimal experimental designs considering only one criterion, though popular, may yield minor improvements in the selected criterion with significant deterioration in other criteria. In the past, the maj ority of the work related to the construction of experimental designs is done by considering only one design obj ective. When noise is the dominant source of error, there are a number of experimental designs that minimize the effect of variance (noise) on the resulting approximation, for example, the Doptimal design, that minimizes the variance associated with the estimates of coefficients of the response surface model. Traditional variance based designs minimize the effect of noise and attempt to obtain uniformity (ratio of maximum to minimum error in design space) over the design space, but they do not address bias errors. Classical minimum bias designs consider only spaceaveraged or integrated error measures (Myers and Montegomery, 1995, pp. 208279) in experimental designs. The bias component of the averaged or integrated mean squared error is minimized to obtain socalled minimum bias designs. The fundamentals of minimizing integrated mean squared error and its components can be found in Myers and Montgomery (1995, Chapter 9), and Khuri and Cornell (1996, Chapter 6). Venter and Haftka (1997) developed an algorithm implementing a minimumbias criterion for irregularly shaped design spaces where no closed form solution exists for experimental design. They compared minimumbias and Doptimal experimental designs for two problems with two and three variables. The minimumbias experimental design was found to be more accurate than Doptimal for average error but not for maximum error. Qu et al. (2004) implemented Gaussian quadraturebased minimum bias design and presented minimum bias central composite designs for up to six variables. There is some work done on developing experimental designs by minimizing the integrated mean squared error accounting for both variance and bias errors. Box and Draper (1963) minimized integrated mean squared errors averaged over the design space by combining average weighted variance and average bias errors. Draper and Lawrence (1965) minimized the integrated mean square error to account for model inadequacies. Kupper and Meydrech (1973) specified bounds on the coefficients associated with the assumed true function to minimize integrated mean squared error. Welch (1983) used a linear combination of variance and bias errors to minimize mean squared error. Montepiedra and Fedorov (1997) investigated experimental designs minimizing the bias component of the integrated mean square error subject to a constraint on the variance component or viceversa. Fedorov et al. (1999) later studied design of experiments via weighted regression prioritizing regions where the approximation is needed to predict the response. Their approach considered both variance and bias components of the estimation error. Bias error averaged over the design space has been studied extensively, but there is a relatively small amount of work to account for pointwise variation of bias errors because of inherent difficulties. An approach for estimating bounds on bias errors in PRS by a pointwise decomposition of the mean squared error into variance and the square of bias was developed by Papila and Haftka (2001). They used the bounds to obtain experimental designs (EDs) that minimize the maximum absolute bias error. Papila et al. (2005) extended the approach to account for the data and proposed datadependent bounds. They assumed that the true model is a higher degree polynomial than the approximating polynomial, and that it satisfies the given data exactly. Goel et al. (2006a) generalized this bias error bounds estimation method to account for inconsistencies between the assumed true model and actual data. They demonstrated that the bounds can be used to develop adaptive experimental designs to reduce the effect of bias errors in the region of interest. Recently, Goel et al. (2006c) presented a method to estimate pointwise root mean square (RMS) bias errors in approximation prior to the generation of data. They applied this method to construct experimental designs that minimize maximum RMS bias error (minmax RMS bias designs). Since minimumbias designs do not achieve uniformity, designs that distribute points uniformly in design space (space filling designs like Latin hypercube sampling) are popular even though these designs have no claim to optimality. Since Latin hypercube sampling (LHS) designs can create poor designs, as illustrated by Leary et al. (2003), different criteria like, maximization of minimum distance between points, or minimization of correlation between points are used to improve its performance. We will demonstrate in this chapter that even optimized LHS designs can occasionally leave large holes in design space, which may lead to poor predictions. Thus, there is a need to consider multiple criteria. Some previous efforts of considering multiple criteria are as follows. In an effort to account for variance, Tang (1993) and Ye (1998) presented orthogonal array based LHS designs that were shown to be better than the conventional LHS designs. Leary et al. (2003) presented strategies to find optimal orthogonal array based LHS designs in a more efficient manner. Palmer and Tsui (2001) generated minimumbias Latin hypercube experimental designs for sampling from deterministic simulations by minimizing integrated squared bias error. Combination of facecentered cubic design and LHS designs is quite widely used (Goel et al., 2006d). The primary obj ective of this work is to demonstrate the risks associated with using a single criterion to construct experimental designs. Firstly, we compare LHS and Doptimal designs, and demonstrate that both these designs can leave large unsampled regions in design space that may potentially yield high errors. In addition, we illustrate the need to consider multiple criteria to construct experimental designs, as singlecriterion based designs may represent extreme tradeoffs among different criteria. Minmax RMS bias design, which yields a small reduction in the maximum bias error at the cost of a huge increase in the maximum variance, is used as an example. While the above issue of tradeoff among multiple criteria requires significant future research effort, we explore several strategies for the simultaneous use of multiple criteria to guard against selecting experimental designs that are optimal according to one criterion but may yield very poor performance on other criteria. In this context, we firstly discuss which criteria can be simultaneously used meaningfully; secondly, we explore how to combine different criteria. We show that complimentary criteria may cater to competing needs of experimental designs. Next, we demonstrate improvements by combining a geometrybased criterion LHS and a modelbased Doptimality criterion, to obtain experimental designs. We also show that poor experimental designs can be filtered out by creating multiple experimental designs and selecting one of them using an appropriate errorbased (pointwise) criterion. Finally, we combine the above mentioned strategies to construct experimental designs. The chapter is organized as follows: Different error measures used in this study are summarized in the next section. Following that, we show maj or results of this study. We illustrate the issues associated with single criterionbased experimental designs, and show a few strategies to accommodate multiple criteria. We close the chapter by recapitulating maj or findings. Error Measures for Experimental Designs Let the true response r(x) at a design point x be represented by a polynomial f "(x)P, where f(x) is the vector of basis functions and p is the vector of coefficients. The vector f(x) has two components: f a)(x) is the vector of basis functions used in the PRS or fitting model, and f (2)(X) iS the vector of additional basis functions that are missing in the linear regression model (assuming that the true model is a polynomial). Similarly, the coefficient vector P can be written as a combination of vectors P'' and p'2 that represent the true coefficients associated with the basis function vectors f '(x) and f (2x), respectively. Precisely, qx) =c f *(x: = "' f ~ ~ (f cz](x)~()" p. + f il' *()" p (3.1) Assuming normally distributed noise E with zero mean and variance G2 (N(0, a2)), the observed response y(x) at a design point x is given as y(x) = r(x)+ E. (3.2) The predicted response f(x) at a design point x is given as a linear combination of approximating basis functions vector f"')(x) with corresponding estimated coefficient vector b : J(x) = (f '(x))T b. (3.3) The estimated coefficient vector b is evaluated using the data (y) for Ns design points as (Myers and Montgomery, 1995, Chapter 2): b = (X() X(1 )1 X() y, (3.4) where X'' is the Gramian matrix constructed using f '(x) (refer to Appendix A). The error at a general design point x is the difference between the true response and the predicted response, e(x) = r(x) f(x) When noise is dominant, estimated standard error ees (x) , used to appraise error, is given as (Myers and Montogomery, 1995) ees(x)Ji~i= j\ar[f~)] of(x)(X X ) f (x) (3.5) where 0,2 is the estimated variance of the noise, and as is the standard error in approximation. When bias error is dominant, the root mean square of bias error e("'"(x) at design point x can be obtained as (Goel et al. 2006c, and Appendix A) e"(;;sa)= E ()) f x)f '"x) E(pp f x)A' "(x (3.6) where Ep(g(x)) is the expected value of g(x) with respect to a, and A is the alias matrix A = (X ')X ')1 X!'!TX '. However, the true model may not be known, and Equation (3.6) is only satisfied approximately for the assumed true model. Prior to generation of data, all components of the coefficient vector p (2) are assumed to have a uniform distribution between 7 and 7 (7 is a constant) such that E,; (b" ) = Y I, where l is an NP2 x NP2 identity matrix. Substituting this in Equation (3.6), the pointwise RMS bias error (Goel et al. 2006c, and Appendix A) is e(""(x) () ff11 wr'" (x)A 7~I f '(x) ' A'f'(x)] (3.7) = df (x)f (x)A f (x)A'f "f(x) Since 7 has a uniform scaling effect, prior to generation of data it can be taken as unity for the sake of simplicity. It is clear from Equation (3.7) that the RMS bias error at a design point x can be obtained from the location of data points (defines alias matrix A) used to fit the response surface, the form (f"') and f 2) of the assumed true function (which is a higher order polynomial than the approximating polynomial), and the constant 7. Goel et al. (2006c) demonstrated with examples that this error prediction model gives good estimates of actual error fields both when the true function is polynomial and when it is nonpolynomial. Two representative examples of a polynomial true function and a nonpolynomial true function trigonometricc function with multiple frequencies) are presented in Appendix B, respectively. Many criteria used to construct/compare different EDs are presented in the literature. A few commonly used error metrics as well as new bias error metrics are listed below. *Maximum estimated standard error in design space (Equation (3.5) with "a = 1) (eema ) max e, (x)). (3.8) *Spaceaveraged estimated standard error (Equation (3.5) with "a = 1) (e le s(x)dx. (3.9) es avg vol (V) , *Maximum absolute bias error bound (Papila et al., 2005) in design space (eI )max = mx e (x)), (3.10) where c(2) = 1 and (e)a ,2)XA/ x 2 (3.11) *Maximum RMS bias error in design space (Equation (3.7) with y = 1) (ebr~ms) = max (erm cx) (3.12) *Spaceaveraged RMS bias error (Equation (3.7) with y = 1) (e re) er:" (x)dx. (3.13) This criterion is the same as spaceaveraged bias error. Among all the above criteria, the standard error based criteria are the most commonly used. For all test problems, a design space coded as an N,,dimensional cube V = [1,1]"' is used and bias errors are computed following the common assumption that the true function and the response surface model are cubic and quadratic polynomials, respectively. Besides the above errormetric based criteria, the following criteria are also frequently used. * Defficiency (Myers and Montgomery, 1995, pp. 393) D, = maxm M) /a; M = XmrXN sy (3.14) Here, max MI in Equation (3.14) is taken as the maximum of all experimental designs. This criterion is primarily used to construct Doptimal designs. A high value of D efficiency is desired to minimize the variance of the estimated coefficients b. * Radius of the largest unoccupied sphere (rmax) We approximate the radius of the largest sphere that can be placed in the design space such that there are no experimental design points inside this sphere. A large value of rmax indicates large holes in the design space and hence a potentially poor experimental design. This criterion is not used to construct experimental designs, but this allows us to measure the spacefilling capability of any experimental design. Test Problems and Results This section is divided into two parts. In the first subsection, we compare widely used experimental designs, like LHS designs, Doptimal designs, central composite designs and their minimum bias error counterparts. We show that different designs offer tradeoffs among multiple criteria and experimental designs based on a single error criterion may be susceptible to high errors on other criteria. In the second subsection, we discuss a few possible strategies to simultaneously accommodate multiple criteria. Specifically, we present two strategies, (1) combination of a geometrybased criterion (LHS) with a modelbased criterion (Doptimality), and (2) simultaneous use of multiple experimental designs combined with pointwise error estimates as a filtering criterion to seek protection against poor designs. Comparison of Different Experimental Designs Space filling characteristics of Doptimal and LHS designs Popular experimental designs, like LHS designs that cater to bias errors by evenly distributing points in design space or numerically obtained Doptimal designs that reduce the effect of noise by placing the design points as far apart as possible, can occasionally leave large holes in the design space due to the random nature of the design (Doptimal) or due to convergence to local optimized LHS designs. This may lead to poor approximation. Firstly, we demonstrate that for Doptimal and optimized LHS designs, a large portion of design space may be left unsampled even for moderate dimensional spaces. For demonstration, we consider two to fourdimensional spaces V' = [1, 1]1.. The number of points in each experimental design is twice the number of coefficients in the corresponding quadratic polynomial, that is, 12 points for two dimensions, 20 points for three dimensions, and 30 points for fourdimensional design spaces. We also create experimental designs with 40 points in fourdimensional design space. We generate 100 designs in each group to alleviate the effect of randomness. Doptimal designs were obtained using the MATLAB routine 'candexch' such that duplicate points are not allowed (duplicate points are not useful to approximate data from deterministic functions or numerical simulations). We supplied a grid of points (grid includes corner/face points and points sampled at a grid spacing randomly selected between 0.15 and 0.30 units) and allocated a maximum of 50000 iterations to find a Doptimal design. LHS designs were constructed using the MATLAB routine 'lhsdesign' with the 'maximin' criterion that maximizes the minimum distance between points. We allocated 40 iterations for optimization. For each experimental design, we estimated the radius (rmax of the largest unsampled sphere that fits inside the design space and summarized the results with the help ofboxplots in Figure 31. The box encompasses the 25th to 75th percentiles and the horizontal line within the box shows the median value. The notches near the median represent the 95% confidence interval of the median value. It is obvious from Figure 31 that rmax increases with the dimensionality of the problem, i.e., the distribution of points in high dimensional spaces tends to be sparse. As expected, an increase in the density of points reduced rmax (compare fourdimensional space with 30 and 40 points). The reduction in rmax was more pronounced for Doptimal designs than LHS designs. LHS designs had a less sparse distribution compared to Doptimal designs, however, the median rmax of approximately 0.75 units in fourdimensional space for LHS designs indicated that a very large region in the design space remained unsampled and data points are quite far from each other. The sparse distribution of points in the design space is illustrated with the help of a three dimensional example with 20 points in Figure 32, where the largest unsampled sphere is shown. For both Doptimal and LHS designs, the large size of the sphere clearly demonstrates the presence of large gaps in the design space that makes the surrogate predictions susceptible to errors. This problem is expected to become more severe in high dimensional spaces. The results indicate that a single criterion (Doptimality for Doptimal designs, and maxmin distance for LHS designs) based experimental design may lead to poor performance on other criteria. Tradeoffs among various experimental designs Next, we illustrate tradeoffs among different experimental designs by comparing minmax RMS bias design (refer to Appendix B), facecentered cubic design (FCCD), Doptimal design (obtained using JMP, Table 31), and LHS design (generated using MATLAB routine 'lhsdesign' with 'maximin' criterion, and allocating 1000 iterations to get a design, Table 32) in fourdimensional space. Note that all experimental designs, except FCCD, optimize a single criterion, i.e., Doptimal designs optimize Defficiency, LHS designs maximize the minimum distance between points, and minmax RMS bias designs minimize the influence of bias errors. On the other hand, FCCD is an intuitive design obtained by placing the points on the faces and vertices. The designs were tested using a uniform 114 grid in the design space V = [1,1I]4 and different metrics are documented in Table 33. We observed that no single design (used in the generic sense, meaning a class of designs) outperformed other designs on 'all' criteria. The D optimal and the facecentered cubic design had high Defficiency; the minmax RMS bias design and the LHS design had low Defficiency. The minmax RMS bias design performed well on bias error based criteria but caused a significant deterioration in standard error based criteria, due to the peculiar placement of axial points near the center. While the Doptimal design was a good experimental design according to standard error based criteria, large holes in the design space (rmax = 1) led to poor performance on bias error based criteria. Since LHS designs neglect boundaries, they resulted in very high maximum standard error and bias errors. However, LHS designs yielded the least spaceaveraged bias error estimate. The FCCD design, which does not optimize any criterion, performed reasonably on all the metrics. However, we note that FCCD designs in high dimensional spaces are not practical due to the high ratio of the number of simulations to the number of polynomial coefficients. We used polynomial examples to illustrate the risks in using experimental designs constructed with a single criterion. Firstly, we considered a simple quadratic function F,(x) (Equation (3.15) in fourdimensional space) with normally distributed noise e (zero mean and unit variance), F, (x) = r(x) + E, 77(x) = 10(1+ xl + x3 (, X2~(.5 We construct a quadratic polynomial response surface approximation using the data at 25 points sampled with different experimental designs (minmax RMS bias design; FCCD; D optimal, Table 31; and LHS, Table 32) and compute actual absolute error in approximation at a uniform grid of 1 14 pOints in the design space V = [1,1I]4 The accuracy measures based on the data(, thatL are most3 commonI~lly used, are theC adjuse 3coefcin ofCIC~CI deeriato RCCIIIIII 'andth standard error normalized by the range of the function (RMSE) in Table 34. Very high values of normalized maximum, root mean square, and average absolute errors (normalized by the range of the data for the respective experimental design) in Table 34 indicate that the minmax RMS bias design (also referred to as RMS bias CCD) is indeed a poor choice of experimental design whenI theC error is due to nisel, thoVughl all approximlation ac~curac~y melasures (Rd sadr error) suggested otherwise. That is, the high errors come with no warning from the fitting process! High values of the ratio of spaceaveraged, root mean square, or maximum actual errors to sandrd rrorindcat th risks,+ associated,,,, wi,:th relying, on measures such as R determine the accuracy of approximations (We pursue this issue in detail in a Chapter 5). Among other experimental designs, LHS design had a high normalized maximum error near the corners, where no data is sampled. FCCD and Doptimal designs performed reasonably, with FCCD being the best design. Secondly, we illustrate errors due to large holes in design space observed in the previous section. A simple function that is likely to produce large errors would be a cosine with maximum at the center of the sphere. However, to ensure a reasonable approximation of the true function with polynomials, we used a truncated Maclaurin series expansion of a translated radial cosine function cos(k Ix xfhs l),namely F(x)=20j 1r pr l; r=k xx d,l (3.16) where Xed is a fixed point in design space, and k is a constant. We considered two instances of Equation (3.16) by specifying the center of the largest unoccupied sphere associated with LHS design (Table 32) and Doptimal design (Table 31) as the fixed point. F2 X)= 20 1~r r %);r = k2XX hsl ,X hs = [0. 168, 0. 168, 0. 141, 0. 167], (.7 F3~X) = 20 1 %lir r ) k3XXot, X optx = [0.0, 0.0, 0.0, 0.0], (.8 The constants k2 and k3 WeTO Obtained by maximizing the normalized maximum actual absolute error in approximation of F2 X) USing LHS experimental design and approximation of F3 X) USing Doptimal experimental design, respectively, subj ect to a reasonable approximation (determined by the condition, R~ > 0.90) of the two functions by all considered experimental designs (FCCD, Doptimal, LHS, and RMS bias CCD). As earlier, the errors were normalized by dividing the actual absolute errors by the range of data values used to construct the experimental design. Subsequently, the optimal value of the constants were k2 = 1.13 and k3 = 1.18 . We used quadratic polynomials to approximate F(x) and errors are evaluated at a uniform grid of 114 pOints in the design space V = [1,1]4 The quality of fit, maximum, root mean square, and average actual absolute errors in approximation for each experimental design are summarized in Table 34. We observed that despite a good quality of fit (high R~ and low normalized standard error), the normalized maximum actual absolute errors were high for all experimental designs. In particular, the approximations constructed using data sampled at the D optimal and the LHS designs performed very poorly. This means that the accuracy metrics, though widely trusted, can mislead the actual performance of the experimental design. The high maximum error in approximations using the LHS designs occurred at one of the corners that was not sampled (thus extrapolation error), however, we note that LHS designs yielded the smallest normalized spaceaveraged and root mean square error in approximation. On the other hand, the maximum error in approximations using Doptimal designs appeared at a test point closest to the center x, in the case of F,(x), and near x ig in the case of F,(x). Besides, high normalized spaceaveraged errors indicated poor approximation of the true function F(x). The other two experimental designs, FCCD and RMS bias CCD, performed reasonably on maximal errors. The relatively poor performance of RMS bias CCD for average and RMS errors is explained by recalling that the experimental design was constructed by assuming the true function to be a cubic polynomial, whereas F(x) was a quartic polynomial. An important characteristic for all experimental designs is the ratio of spaceaveraged, root mean square, or maximum actual absolute error to estimated standard error. When this ratio is large, the errors are unexpected and therefore, potentially damaging. The FCCD design provided a reasonable value of the ratio of actual to estimated standard errors, however, RMS bias design performed very poorly as the actual errors were much higher than the standard estimated error. This means that the estimated standard error is misleading about the actual magnitude of error that cannot be detected in an engineering example where we do not have the luxury of using a large number of data points to test the accuracy of approximation. Similarly, for all functions, the ratio of maximum actual absolute error to standard error for LHS designs (2952) was much higher than for Doptimal designs (about 9). The surprise element is also evident by the excellent values of R~ of 0.99 and 1.00 compared to 0.90 for the Doptimnal design. The results presented here clearly suggest that different experimental designs were non dominated with respect to each other and offered multiple (sometimes extreme) tradeoffs, and that it might be dangerous to use a single criterion based experimental design without thorough knowledge of the problem (which is rarely the case in practice). Extreme example of risks in single criterion based design: Minmax RMS bias CCD A more pathological case, demonstrating the risks in developing experimental designs using a single criterion, was encountered for moderate dimensional cases while developing the central composite design counterpart of the minimum bias design, i.e., minimizing the maximum RMS bias error. The performance of the minmax RMS bias designs constructed using two parameters (refer to Appendix B) for two to fivedimensional spaces on different metrics is given in Table 35. For two and threedimensional spaces, the axial points (given by at) were located at the face and the vertex points (given by al) were placed slightly inwards to minimize the maximum RMS bias errors. The RMS bias designs performed very reasonably on all the error metrics. A surprising result was obtained for optimal designs for four and fivedimensional spaces: while the parameter corresponding to vertex points (al) was at its upper limit (1.0), the parameter corresponding to the location of axial points (az) hit the corresponding lower bound (0.1). This meant that to minimize maximum RMS bias error, the points should be placed near the center. The estimated standard error was expectedly very high for this design. Contrasting facecentered cubic design for fourdimensional cases with three and fourdimensional minmax RMS bias designs (Table 35) isolated the effect of dimensionality and the change in experimental design (location of axial points) on different error metrics. The increase in bias errors (bounds and RMS error) was attributed to increase in dimensionality (small variation in bias errors with different experimental designs in fourdimensional design space), and the increase in standard error for minmax RMS bias design was the outcome of the change in experimental design (the location of axial points given by a2). This unexpected result for four and higher dimensional cases is supported by theoretical reasoning (Appendix B), and very strong agreement between the predicted and the actual RMS bias errors for the minmax RMS bias design and the facecentered central composite design (Appendix B). To further illustrate the severity of the risks in using a single criterion, we show the tradeoffs among the maximum errors (RMS bias error, estimated standard error, and bias error bound) for a fourdimnensional design space [1,I]4, Obtained by varying the location of the axial points (a2) fTOm near the center (a2=0.1i, minmax RMS bias design) to the face of the design space (a2=1.0, central composite design), while keeping the vertex locations (al=1.0) fixed. The tradeoff between maximum RMS bias error and maximum estimated standard error is shown in Figure 33(A), and the tradeoff between maximum RMS bias error and maximum bias error bound is shown in Figure 33(B). Moving the axial point away from the center reduced the maximum bias error bound and the maximum estimated standard error but increased the maximum RMS bias error. The relatively small variation in maximum RMS bias error compared to the variation in maximum estimated standard error and maximum bias error bound demonstrated the low sensitivity of maximum RMS bias error with respect to the location of axial points (a2), and explains the success of the popular central composite designs (a2=1.0) in handling problems with bias errors. While we noted that each design on the curves in Figure 33 corresponds to a nondominated (tradeoff) point, a small increase in maximum RMS bias error permits a large reduction in maximum estimated standard error, or in other words, the minimization with respect to a single criterion (here maximum RMS bias error) may lead to small gains at a cost of significant loss with respect to other criteria. Tradeoff between maximum bias error bound and maximum RMS bias error also reflected similar results, though the gradients were relatively small. The most important implication of the results presented in this section is that it may not be wise to select experimental designs based on a single criterion. Instead, tradeoff between different metrics should be explored to find a reasonable experimental design. While detailed exploration of this issue requires significantly more research, our initial attempts to simultaneously accommodate multiple criteria are illustrated next. Strategies to Address Multiple Criteria for Experimental Designs As discussed in the previous section, the experimental designs optimized using a single criterion may perform poorly on other criteria. While a bad experimental design can be identified by visual inspection in low dimensional spaces, we need additional measures to filter out bad designs in high dimensional spaces (Goel et al., 2006e). We explored several strategies to simultaneously accommodate multiple criteria in an attempt to avoid poor experimental designs. In this context, we discuss two issues: * Which criteria are meaningful for different experimental designs? and * How can we combine different criteria? Since the experimental designs are constructed to minimize the influence of bias error and noise, a sensible choice of suitable criteria for any experimental design should seek balance among the two sources of errors, i.e., bias and noise. Consequently, if we select an experimental design that primarily caters to one source of error, for example, noise, the secondary criterion should be introduced to address the other source of error, bias error in this case, and viceversa. We elaborate on this idea in a following subsection. Once we have identified criteria to construct experimental designs, we seek ways to combine different criteria. Taking inspiration from multiobj ective optimization problems, we can accommodate multiple criteria according to several methods, for example, * Optimize the experimental design to minimize a composite function that represents the weighted sum of criteria, * Optimize the experimental design to minimize the primary criterion while satisfying constraints on the secondary criteria, and * Solve a multiobj ective optimization problem to identify different tradeoffs and then select a design that suits the requirements the most. Here, we show two ways to avoid poor experimental designs using a fourdimensional example. Firstly, we present a method to combine the modelbased Doptimality criterion that caters to noise with the geometrybased LHS criterion that distributes points evenly in design space and reduces spaceaveraged bias errors. Secondly, we demonstrate that selecting one out of several experimental designs according to an appropriate pointwise errorbased criterion reduces the risk of obtaining poor experimental designs. Further, we show that the coupling of multiple criteria and multiple experimental designs may be effective to avoid poor designs. Combination of modelbased Doptimality criterion with geometry based LHS criterion We used an example of constructing an experimental design for a fourdimensional problem with 30 points (response surface model and assumed true model were quadratic and cubic polynomials, respectively). Three sets of experimental designs were generated as follows. The first set comprised 100 LHS experimental designs generated using the MATLAB routine 'lhsdesign' with the 'maximin' criterion (a maximum of 40 iterations were assigned to find an experimental design). The second set comprised 100 Doptimal experimental designs generated using the MATLAB routine 'candexch' with a maximum of 40 iterations for optimization. A grid of points, (with grid spacing randomly selected between 0. 15 and 0.30) including face and corner points, was used to search for the Doptimal experimental designs. The third set of (combination) experimental designs was obtained by combining Doptimal (modelbased criterion) and LHS designs (geometrybased criterion). We selected 30 design points from a 650 point LHS design ('lhsdesign' with 'maximin' criterion and a maximum of 100 'iterations' for optimization) using the Doptimality criterion ('candexch' with a maximum of 50000 iterations for optimization). For each design, we computed the radius rmax of the largest unsampled sphere, Defficiency, maximum and spaceaveraged RMS bias and estimated standard error using a uniform 114" grid in the design space [1,I]4 We show the tradeoff among different criteria for Doptimal, LHS, and combination designs in Figure 34. As can be seen from Figure 34(A), the Doptimal designs were the best and LHS designs were the worst with respect to the maximum estimated standard and RMS bias error. Compared to the LHS designs, the combination designs significantly reduced the maximum estimated standard error with marginal improvement on the maximum RMS bias error criterion (Figure 34(A)), and improved Defficiency without sacrificing rmax (Figure 34(D)). The advantages of using combination designs were more obvious in Figure 34(B), where we compared spaceaveraged bias and estimated standard errors. We see that Doptimal designs performed well on spaceaveraged estimated standard errors but yielded high spaceaveraged RMS bias errors. On the other hand, the LHS designs had low spaceaveraged RMS bias errors but high spaceaveraged estimated standard errors. The combination designs simultaneously yielded low spaceaveraged RMS bias and estimated standard errors. This result was expected because the Latin hypercube sampling criterion allows relatively uniform distribution of the 2 The average number of points in the uniform grid used to generate Doptimal designs was 1300. So to provide a fair comparison while keeping the computational cost low, we obtain 650 points using LHS and use this set of points to develop combination experimental designs. points by constraining the location of points that are used to generate combination designs using the Doptimality criterion. Similarly, we observed that unlike Doptimal designs, combination experimental designs performed very well on the spaceaveraged RMS bias error and the r;;2a criterion (refer to Figure 34(C)), and the performance was comparable to that of the LHS designs. Mean and coefficient of variation (COV) of different metrics for the three sets of experimental designs are tabulated in Table 36. Doptimal designs outperformed LHS designs in terms of the ratio of maximum to average error (stability), Defficiency, maximum RMS bias error, and maximum estimated standard error. Also, for most metrics, the variation in results due to sampling (COV) was the least among the three. As seen before, LHS designs performed the best on spaceaveraged RMS bias errors. The designs obtained by combining two criteria (D optimality and LHS), were substantially closer to the best of the two except for (ens")m Thus, they reduced the risk of large errors. Furthermore, the variation with samples (COV) is also reduced. The results suggested that though different experimental designs were nondominated (tradeoffs) with respect to each other, simultaneously considering multiple criteria by combining the modelbased Doptimality criterion and the geometrybased LHS criterion may be effective in producing more robust experimental designs with a reasonable tradeoff between bias errors and noise. Multiple experimental designs combined with pointwise errorbased filtering Next, we demonstrate the potential of using multiple experimental designs to reduce the risk of finding poor experimental designs. The main motivation is that the cost of generating experimental designs is not high so we can construct two or three experimental designs using LHS, or Doptimality, or a combination of the two criteria, and pick the best according to an appropriate criterion. To illustrate the improvements by using three EDs over a single ED, each of the two criteriamaximum RMS bias error and maximum estimated standard errorwere used to select the best (least error) of the three EDs. For illustration, 100 such experiments were conducted with LHS designs, Doptimal designs, and the combination of LHS and Doptimal designs (as described above). Actual magnitudes of maximum RMS bias error and maximum estimated standard error for all 300 designs and the 100 designs obtained after filtering using minmax RMS bias or maximum estimated standard error criteria are plotted in Figure 35 for three sets of (100) experimental designs. As is evident by the shortening of the upper tail and the size of the boxplots in Figure 35, both the minmax RMS bias and maximum estimated standard error criteria helped eliminate poor experimental designs for all three sets. Filtering had more impact on the maximum error estimates than the spaceaveraged error estimates. The numerical quantifieation of improvements in actual magnitude of maximum and spaceaveraged error based on 100 experiments is summarized in Table 37. We observed that the pointwise errorbased (minmax RMS bias or estimated standard error) Eiltering significantly reduced the mean and COV of maximum errors. We also noted improvements in the individual experimental designs using multiple criteria. LHS designs were most significantly improved by picking the best of three based on estimated maximum standard error. Doptimal designs combined with the min max RMS bias error based filtering criterion helped eliminate poor designs according to the RMS bias error criterion. It can be concluded from this exercise that potentially poor designs can be filtered out by considering a small number of experimental designs with an appropriate (min max RMS bias or maximum estimated standard) error criterion. The filtering criterion should be complimentary to the criterion used for construction of the experimental design, i.e., if a group of EDs are constructed using a variance based criterion, then the selection of an ED from the group should be based on bias error criterion, and viceversa. Results presented in this section indicate that use of multiple criteria (LHS and D optimality) and multiple EDs help reduce maximum and spaceaveraged bias and estimated standard errors. Implementing the above findings, we can obtain experimental designs with reasonable tradeoff between bias error and noise in three steps as follows: * Generate a large number of LHS experimental design points, * Select a Doptimal subset within the LHS design (combine modelbased and geometry based criteria), * Repeat first two steps three times and select the design that is the best according to one of the minmax RMS bias or maximum estimated standard error criteria (filtering using pointwise errorbased criterion). Concluding Remarks In this chapter, we demonstrated the risks of using a single criterion to construct experimental designs. We showed that constructing experimental designs by combining multiple (model, geometry, and error based) criteria and/or using multiple experimental designs reduces the risk of using a poor experimental design. For fourdimensional space, comparison of computed LHS and Doptimal designs, that involve random components and may yield poor approximation due to random components or convergence to local optima, revealed that the Doptimal designs were better for maximum errors, and LHS designs were better for spaceaveraged RMS bias errors. Both designs were susceptible to leaving large spheres in design space unsampled. A comparison of popular experimental designs (facecentered cubic design, minmax RMS bias design, Doptimal design, and LHS design) revealed the nondominated (tradeoff among different criteria) nature of different designs. The minmax RMS bias design, obtained by placing the axial points close to the center, performed the best in reducing maximum RMS bias error, but was the worst design for estimated standard error metrics and Defficiency. LHS design gave the best performance in terms of spaceaveraged bias errors. However, facecentered cubic design that is an intuitive design yielded a reasonable tradeoff between bias error and noise reduction on all metrics. The same conclusions were supported by approximation of three example polynomials that highlighted the susceptibility of different experimental designs to the nature of the problem, despite the fact that the accuracy metrics suggested a very good fit for each example. We concluded that different experimental designs, constructed using one error criterion, do not perform the best on all criteria. Instead, they offer tradeoffs. In moderate dimensional spaces these single criterionbased designs can often lead to extreme tradeoffs, particularly by using the maximum RMS bias error measure as a design criterion, such that small gains in the desired criterion are achieved at the cost of significant deterioration of performance in other criteria. A tradeoff study, conducted to study the variation of different error metrics with the location of axial points in centralcomposite designs, illustrated the perils of using a single criterion to construct experimental designs and emphasized the need to consider multiple criteria to tradeoff bias error and noise reduction. To address the risk of using a poor experimental design by considering a single criterion, we explored a few strategies to accommodate multiple criteria. We demonstrated that the experimental design obtained by combining two criteria, the Doptimality criterion with LHS design, offered a reasonable tradeoff between spaceaveraged RMS bias and estimated standard error, and spacefi11ing criteria. Specifically, combination designs significantly improved the poor experimental designs. We showed that the risk of getting a poor experimental design could be further reduced by choosing one out of three experimental designs using a pointwise error based criterion, e.g., minmax RMS bias or maximum estimated standard error criterion. The combination of Doptimal designs and minmax RMS bias error was particularly helpful in reducing bias errors. Finally, we adopted selection of experimental designs by combining the D optimality criterion with LHS design and selecting one out of three such combination designs to cater to both bias error and noise reduction. However, since these results are based on a limited number of examples, we note the need of future research to address the issues related to accommodating multiple criteria while constructing experimental designs. A 2D1 02 D30 4022 32 D3 D4 Fiue . oplt (aedo 10deins frais rna)ofte agetuncupe shr inside~~~~~~ th einsae[1 ]'(hr ,i tenme fvrals.xai hw th dmesinliy f h dsin pcean cresonin umerofp int inth expeimetaldesgn.Smalerr~n i desredto voi lage nsaple re ion. D... opia dein are selected usn MALA rotn 'cnech w peiyardo pint wihgi pcn ewe .5ad030 n aiu f500ieain Figure for Boxptimizatidon) LHS designs are gnratd usin MATLAB routine 'lgs uocpshsesign ponswith gi pcn ewe .1 n .0 n a maximum of 4000 iterations frotmzto.A pia ein.B H designs. A 11 1 B Figure 32. Illustration of the largest spherical empty space inside the 3D design space [1, 1]3 (20 points). A) Doptimal designs. B) LHS designs. 1.186, , , a Axial point~ 'nearceiffer~ 1.175 1.175 1.16 1.155 0 10 20 30 40 50 60 70 80 2 4( 86 68 7 7.2 A max(e ~ ~maxceb. B Figure 33. Tradeoffs between different error metrics. A) Maximum estimated standard error (es)max and maximum RMS bias error (e nn)max B) Maximum bias error bound (eb)ma and maximum RMS bias error (e nn)max for fourdimensional space. 25 points were used to construct central composite experimental designs with vertex location fixed at al= 1.0. Axial point on thei~e au2 21 0.95 0.9 0.75 0.65 0.5 0.5 *  o * Doptimal * LHS SCormbination ** = On 0 1 2 34 max(ee,(x)) 5 6 0.6 0.7 0.8 0.9 1 1.1 avg(e (x)) 1.2 1.3 B 4 * * ** * q* ** * ** I. *Doptimnal *LHS u combination 0O 0.7 Doptimal LHS 165 Combination 0965 0.75 0 8 0.85 0 9 0.9 1 0 a 0 5 0 6 0; ne 0 80 Defliciency Figure 34. Comparison of 100 Doptimal, LHS, and combination (Doptimality + LHS) experimental designs in fourdimensional space (30 points) using different metrics. rmax: radius of the largest unsampled sphere, e(""(x): RMS bias error, e,s(x): estimated standard error, (.)max: maximum of the quantity inside parenthesis, (.)avg: Spatial average of the quantity inside parenthesis. All metrics, except Defficiency (Dyf), are desired to be low. A) (ees)max vs. (e ""s)max. B) (ees)civg VS. (e tas)civg C) rmax VS. (e tas)civg D) D~f vs. rmax r Doptimal LHS o Combins~ion 0.95 0 9 0.851 " 0.8 62 0 3 7 6.5 5.5 4.5 3. 2.5 mBE mSE A AIILHS designs 2.2t ~ 1.1 1 aBE aSE AII LHS designs 1.1C aBE aSE F.srering using SE 1 1t 1 5I mBE mSE Filtering using SE 5 mBE mSE Filtering using BE I I aBE aSE F.Isrerig ursing BE 1.11 1 2.2 ~ 2.2C 0g~ 8 0 7C 1.4) . 1.41 0.7) 01.8 0.8) 0.8) 0 o5~ aBE aSE aBE aSE Filtering using BE Filtering using SE 0.6t , mBE mSE 0.6~ 0.6C 05~ _ mBE mSE mBE mSE aBE aSE B AII Deoptimal designs Filtering using BE Filtering using SE AII Doptimnal destgns mBE mSE Al P1Doptimal designs mBE mSE Filtering using BE mBE mSE Filtering using SE aBE aSE aBE aSE All Doptimaldesigns Filtering using BE aBE aSE Filtering using SE Figure 35. Simultaneous use of multiple experimental designs concept, where one out of three experimental designs is selected using appropriate criterion (filtering). Boxplots show maximum and spaceaveraged RMS bias, and standard estimated errors in four dimensional design space; considering all 300 designs, 100 designs filtered using minmax RMS bias error as a criterion and 100 designs filtered using min max estimated standard error as a criterion (mBE denotes maximum RMS bias error, mSE 