Response Surface Models are often used as a surrogate for expensive black-box functions during optimization to reduce computational cost. Often, the CAE analysis models are highly nonlinear and multi-modal. A response surface approximation of such analysis as a result is highly multi-modal; i.e. it contains multiple local optima. A gradient-based optimizer working with such a response surface will often converge to the nearest local optimum. There does not exist any method to guarantee a global optima for non-convex multi-modal functions. For such problems, we propose an efficient algorithm to find as many distinct local optima as possible. The proposed method is specifically designed to work in large dimensions (about 100 ∼ 1000 design variables and similar number of constraints) and can identify most of the locally optimal solutions in a reasonable amount of time. We also present a method to analyze natural partitions in a noisy data set and use it to characterize the locally optimal basins found by our method.