Response surface methodology (RSM) relies on the design of experiments and empirical modelling techniques to find the optimum of a process when the underlying fundamental mechanism of the process is largely unknown. This paper proposes an iterative RSM framework, where Gaussian process (GP) regression models are applied for the approximation of the response surface. GP regression is flexible and capable of modelling complex functions, as opposed to the restrictive form of the polynomial models that are used in traditional RSM. As a result, GP models generally attain high accuracy of approximating the response surface, and thus provide great chance of identifying the optimum. In addition, GP is capable of providing both prediction mean and variance, the latter being a measure of the modelling uncertainty. Therefore, this uncertainty can be accounted for within the optimization problem, and thus the process optimal conditions are robust against the modelling uncertainty. The developed method is successfully applied to the optimization of trans-stilbene conversion in the epoxidation of trans-stilbene over cobalt ion-exchanged faujasite zeolites (Co2+-NaX) catalysts using molecular oxygen. (C) 2009 Elsevier B.V. All rights reserved.