Many global optimization tasks have partially or completely shared search spaces. This is particularly common in computational heterogeneous catalysis, where catalysts are optimized to their lowest-energy configurations under various idealized chemical conditions to construct phase diagrams. If each task represents the global optimization under a unique reaction condition, all tasks essentially shares the same configurational search space as well as the same black-box function that takes as input the configuration and returns the relaxed energy obtained by an expensive electronic structure calculation. Bayesian optimization represents an efficient paradigm for solving a single global optimization task of an expensive black-box function. Evolutionary multitasking is a newly emerging paradigm for solving multiple global optimization tasks simultaneously. Here we present a Bayesian evolutionary multitasking (BEM) framework combining the best of the two paradigms. As a case study, we aim to derive the phase diagram of Pt-doped Ni nanoparticles covering a wide variety of steam methane reforming conditions. By integrating knowledge such as graph theory and lattice symmetry, the BEM framework as a whole can significantly accelerate the mapping of phase diagrams for very complex heterogeneous systems. The efficient calculation of surface free energies of low index facets in the reactive catalyst environment provides insight into catalyst deactivation through coking and sulfur poisoning.