Gaussian process (GP) model based optimization is widely applied in simulation and machine learning. In general, it first estimates a GP model based on a few observations from the true response and then uses this model to guide the search, aiming to quickly locate the global optimum. Despite its successful applications, it has several limitations that may hinder its broader use. First, building an accurate GP model can be difficult and computationally expensive, especially when the response function is multimodal or varies significantly over the design space. Second, even with an appropriate model, the search process can be trapped in suboptimal regions before moving to the global optimum because of the excessive effort spent around the current best solution. In this work, we adopt the additive global and local GP (AGLGP) model in the optimization framework. The model is rooted in the inducing points based GP sparse approximations and is combined with independent local models in different regions. With these properties, the AGLGP model is suitable for multimodal responses with relatively large data sizes. Based on this AGLGP model, we propose a combined global and local search for optimization (CGLO) algorithm. It first divides the whole design space into disjoint local regions and identifies a promising region with the global model. Next, a local model in the selected region is fit to guide detailed search within this region. The algorithm then switches back to the global step when a good local solution is found. The global and local natures of CGLO enable it to enjoy the benefits of both global and local search to efficiently locate the global optimum. Summary of Contribution: This work proposes a new Gaussian process based algorithm for stochastic simulation optimization, which is an important area in operations research. This type of algorithm is also regarded as one of the state-of-the-art optimization algorithms for black-box functions in computer science. The aim of this work is to provide a computationally efficient optimization algorithm when the baseline functions are highly nonstationary (the function values change dramatically across the design space). Such nonstationary surfaces are very common in reality, such as the case in the maritime traffic safety problem considered here. In this problem, agent-based simulation is used to simulate the probability of collision of one vessel with the others on a given trajectory, and the decision maker needs to choose the trajectory with the minimum probability of collision quickly. Typically, in a high-congestion region, a small turn of the vessel can result in a very different conflict environment, and thus the response is highly nonstationary. Through our study, we find that the proposed algorithm can provide safer choices within a limited time compared with other methods. We believe the proposed algorithm is very computationally efficient and has large potential in such operational problems.