This doctoral dissertation focuses on multiobjective stochastic simulation optimization, and has been developed on the on the interface of two key research fields: operations research and artificial intelligence. Operations research models are commonly applied to support decisions in complex business and industrial systems (e.g., operations management, transportation science, or healthcare management). Often, the decision maker has different objectives to optimize simultaneously (e.g., minimize the waiting time of patients in the emergency department, while also minimizing personnel costs; or, alternatively, minimize the cost of a production system, while also minimizing the profit risk and maximizing the service level). This problem is known as multiobjective optimization (MO), and the goal is to find a solution set that reflects the optimal trade-offs between the different conflicting objectives (referred to as the Pareto-optimal set). Moreover, the objectives also tend to be stochastic (i.e., they contain uncertainty), and they cannot be directly observed (i.e., the performance can only be estimated). Often, stochastic simulation is used as a tool to study such systems: the simulation model numerically evaluates the objectives for a given decision vector. Despite its relevance in practice, the current literature on multiobjective stochastic simulation optimization is scarce. In particular, we focus on settings where the simulation is computationally expensive (i.e., it takes a long time to run a single replication) and the decision maker has limited simulation budget. Ideally, an algorithm should search the solution space in a way that is both effective (i.e., having high probability of finding the Pareto-optimal solutions), and efficient (using limited computational budget). To circumvent the expensive simulation runs, different techniques can be used to approximate the simulation outcomes through metamodels. In particular, as it is common in Bayesian optimization and machine learning, we have exploited the use of Gaussian Processes, which subsume kriging metamodels as predictive models, to account for the intrinsic uncertainty perturbing the objectives during optimization. In Chapter 1 we explicitly establish the link between mathematical optimization (subject of study in the operations research field), and kriging metamodels (subject of study in the artificial intelligence domain); we then discuss the positioning of multiobjective optimization from a Bayesian rationale, and the role of machine learning in this dissertation. In Chapter 2 we outline the essential theoretical background. In Chapter 3 we survey the most relevant kriging-based algorithms for multiobjective simulation optimization. These algorithms perform a sequential search of so-called infill points, used to update the kriging metamodel at each iteration. An infill criterion helps to balance local exploitation and global exploration during this search by using the information provided by the kriging metamodels. Most research has been done on algorithms for deterministic problem settings; only very recently, algorithms for stochastic simulation outputs have been proposed. The current literature for stochastic problems shows a significant shortcoming: all algorithms either ignore the intrinsic noise in the outcomes that stems from the stochastic simulation (by fitting deterministic kriging metamodel to the mean observed outcomes, treating this mean as if it were sampled with infinite precision), or they assume that the noise is homogeneous (i.e., that it has the same strength over the search space). In practical settings, though, the noise is usually heterogeneous. In Chapter 4, we propose a multiobjective simulation optimization algorithm that contains two crucial elements: the search phase implements stochastic kriging to account for the inherent noise in the outputs when constructing the metamodel, and the accuracy phase uses a well-known multiobjective ranking and selection procedure in view of maximizing the probability of selecting the true Pareto-optimal points by allocating extra replications on competitive solutions. We evaluate the impact of these elements on the search and identification effectiveness, for a set of test functions with different Pareto front geometries, and varying levels of heterogeneous noise. Our results show that the use of stochastic kriging is essential in improving the search efficiency; yet, the allocation procedure appears to lose effectiveness in settings with high noise. This emphasizes the need for further research on multiobjective ranking and selection methods. In Chapter 5 we consider multiobjective ranking and selection problems with heterogeneous noise and correlation between the mean values of alternatives. From a Bayesian perspective, we propose a sequential sampling technique that uses a combination of screening, stochastic kriging metamodels, and hypervolume estimates to decide how to allocate samples. Extensive experimental results show the proposed sampling criteria to perform relatively well on their own, but the best performance was obtained when combined. Extensive computational experiments show that the proposed method only requires a small fraction of samples compared to the standard equal allocation method, and it is competitive to the state-of-the-art, with the exploitation of the correlation structure being the dominant contributor to the improvement. In Chapter 6 we summarize the most relevant findings of the dissertation, and put forward a considerable number of challenging research questions that present valuable opportunities for further work.