Crystallization process operated in MSMPR mode is of great relevance in manufacturing because of its low operational requirements and easy maintenance. However, to enhance product purity and downstream operations, one needs to handle trade-off among residence time, mean and spread of crystal size distribution. In this paper, we, therefore, study this multi-objective optimization problem (MOOP) to obtain the Pareto solutions using high-fidelity (HF) population balance equations. We parallely aim to develop an optimal Pareto front using Multi-objective Bayesian Optimization (MOBO) technique which exploits the qEHVI as the acquisition function to identify the optimization path. The results show that MOBO requires only 2.4% of HF functional calls compared to evolutionary approaches to generate similar quality results. This indicates the potential of MOBO in optimizing the time expensive MOOPs with minimum number of function calls, which can be easily extended for response surface-based optimization studies involving design of experiments. © 2023 Taylor & Francis.