Understanding wood characteristics is essential for enhancing the economic sustainability of the kraft process. In a recent work, we included wood characteristics in the phenomenological modeling of digesters. The developed model allows the optimization of the process. Here, we propose a constrained multi-objective optimization approach to minimize the residual lignin content while maximizing the carbohydrate content in the cellulose pulp, using the Non-dominated Sorting Genetic Algorithm II. In order to reduce computational costs during the optimization, we developed a surrogate model, based on simulated data, to predict key performance indicators. We evaluated two machine learning techniques: Multilayer Perceptron (MLP) and decision tree-based methods, considering single and multiple outputs. The combination between the multiple-output approach and MLP performed well and resulted in a more simplified surrogate model. The results for the simulation and the surrogate model were similar, and the computational cost was approximately 50 times lower for the surrogate model. An uncertainty was associated with the multiplicative factor of the delignification rate as a way to observe how its estimate interferes with the Pareto curve. The present work presents contributions to data-driven modeling and multi-objective optimization of digesters, allowing us to determine the importance of wood characteristics in operation.