Multi-objective optimization of complex reaction network models is often demanding due to the large computational expense for such calculations. We consider one such model for long-chain branched polyvinyl acetate (PVAc-LCB) production consisting of a large set of stiff differential equations. In such a case, where the evaluation of the solution becomes expensive, optimization of the problem using a relatively inexpensive surrogate is a feasibility to explore. We propose a Bayesian optimization framework that balances exploitation (maximizing information from the current best solution from a set of candidates) and exploration (minimizing the uncertainty in unexplored landscape) to solve this extremely complex problem. The Gaussian Process-based surrogate model provides the landscape that is evaluated using a suitable acquisition function for sampling the next best location towards the global optima. We also present a comparison study between the results of the proposed Bayesian approach and that obtained using the Non-Dominated Sorting Genetic Algorithm-II (NSGA-II) without the use of a surrogate. The proposed strategy builds a Pareto front that is comparable to the high-fidelity Pareto front with a computational gain of almost 100-fold. © 2022 IEEE.