Efficient global optimization of expensive black-box functions is a challenging problem in engineering optimization where the number of function evaluations is limited by time or cost. This paper introduces a response surface methodology that is particularly effective at modeling nonlinear, multimodal functions common in engineering. The approach involves fitting response surfaces to data collected from evaluating objective and constraint functions at a few points. These surfaces are then used for visualization, tradeoff analysis, and optimization. The key to using response surfaces for global optimization is balancing the need to exploit the approximating surface (by sampling where it is minimized) with the need to improve the approximation (by sampling where prediction error may be high). This balance requires solving auxiliary problems that have previously been considered intractable, but we show how these computational obstacles can be overcome. The response surface methodology used here is based on modeling the objective and constraint functions with stochastic processes, which may seem complex but is intuitively based on calibrating a model that summarizes how the function typically behaves. This approach has a long history in mathematical geology, global optimization, and statistics. In mathematical geology, it is called 'kriging' and dates back to the early 1960s. In global optimization, it is called 'Bayesian global optimization' or the 'random function approach', dating back to a seminal article by Harold Kushner in 1964. In statistics, it began in the early 1970s with a focus on approximating integrals and other hard-to-compute functionals of functions. Most recently, the focus has been on developing accurate approximations to expensive computer codes for visualization and optimization. In our work, we take the stochastic process model commonly used in the statistics literature and apply it to global optimization. Two things set us apart from previous work in Bayesian global optimization.Efficient global optimization of expensive black-box functions is a challenging problem in engineering optimization where the number of function evaluations is limited by time or cost. This paper introduces a response surface methodology that is particularly effective at modeling nonlinear, multimodal functions common in engineering. The approach involves fitting response surfaces to data collected from evaluating objective and constraint functions at a few points. These surfaces are then used for visualization, tradeoff analysis, and optimization. The key to using response surfaces for global optimization is balancing the need to exploit the approximating surface (by sampling where it is minimized) with the need to improve the approximation (by sampling where prediction error may be high). This balance requires solving auxiliary problems that have previously been considered intractable, but we show how these computational obstacles can be overcome. The response surface methodology used here is based on modeling the objective and constraint functions with stochastic processes, which may seem complex but is intuitively based on calibrating a model that summarizes how the function typically behaves. This approach has a long history in mathematical geology, global optimization, and statistics. In mathematical geology, it is called 'kriging' and dates back to the early 1960s. In global optimization, it is called 'Bayesian global optimization' or the 'random function approach', dating back to a seminal article by Harold Kushner in 1964. In statistics, it began in the early 1970s with a focus on approximating integrals and other hard-to-compute functionals of functions. Most recently, the focus has been on developing accurate approximations to expensive computer codes for visualization and optimization. In our work, we take the stochastic process model commonly used in the statistics literature and apply it to global optimization. Two things set us apart from previous work in Bayesian global optimization.