To efficiently tackle the optimal operation problem of multi-stakeholder integrated energy systems (IESs), this paper develops a multi-objective optimization and multi-attribute decision-making support method. Mathematically, The optimal operation of IESs interconnected with distributed district heating and cooling units (DHCs) via the power grid and gas network, can be formulated as a multi-objective optimization problem considering both economic, reliability and environment-friendly objectives with numerous constraints of each energy stakeholder. Firstly, a multi-objective group search optimizer with probabilistic operator and chaotic local search (MPGSO) is proposed to balance global and local optimality during the random search iteration. The MPGSO utilizes a crowding probabilistic operator to select producers to explore areas with higher potential but less crowding and reduce the number of fitness function calculations. Moreover, a new parameter selection strategy based on chaotic sequences with limited computational complexity is adopted to escape the local optimal solutions. Consequently, a set of superior Pareto-optimal fronts could be obtained by the MPGSO. Subsequently, a multi-attribute decision-making support method based on the interval evidential reasoning (IER) approach is used to determine a final optimal solution from the Pareto-optimal solutions, taking multiple attributes of each stakeholder into consideration. To verify the effectiveness of the MPGSO, the DTLZ suite of benchmark problems are tested compared with the original GSOMP, NSGA-II and SPEA2. Additionally, simulation studies are conducted on a modified IEEE 30-bus system connected with distributed DHCs and a 15-node gas network to verify the proposed approch. The quality of the obtained Pareto-optimal solutions is assessed using a set of criteria, including hypervolume (HV), generational distance (GD), and Spacing index, among others. Simulation results show that the number of Pareto-optimal solutions (NPS) of MPGSO are higher by about 32.6 %-62.1%, computation time (CT) are lower by about 2.94 %-46.1 % compared with other algorithms. Besides, to further evaluate the performance of the proposed approach in addressing larger-scale issues, the study employs the modified IEEE 118-bus system of greater magnitude. The proposed MPGSO algorithm effectively handles multi-objective and non-convex optimization problems with Pareto sets in terms of better convergence and distributivity.