Abstract: To improve the decomposition level of multiple loads in integrated energy systems and the overall performance of prediction models, a method for predicting multiple loads in integrated energy systems is proposed, considering composite index optimization modal decomposition and Stacking integration. Firstly, permutation entropy combined with mutual information is used as the fitness function, and the Golden Jackal Optimization algorithm is used to obtain the optimal parameter combination for variational modal decomposition adaptively. Then, the multiple load sequences are decomposed into a collection of intrinsic mode functions. Secondly, the Mean Impact Value (MIV) algorithm based on BP neural network perturbation is used to screen features related to meteorological, date, and load factors associated with multiple loads, thereby constructing a highly coupled feature matrix for multiple loads. Considering each single model's differences and advantages, a stacking ensemble learning model is constructed to predict multiple loads by reducing overfitting using the k-fold cross-validation method. Finally, the proposed method is validated using a multiple-load dataset from the Tempe campus of Arizona State University in the United States. The results show that the mean absolute percentage errors of the proposed method in predicting electrical, cooling, and heating loads are 0.903%, 2.713%, and 1.616%, respectively. The prediction accuracy is considerably improved compared to other prediction models.