The inherent intermittency and variability of renewable energy sources present significant challenges to the optimal design and implementation of integrated energy systems (IES). This paper introduces a novel stochastic optimization model that integrates advanced scenario generation and clustering algorithm for renewable energy sources within a multi-objective, bi-level optimization framework. Specifically, the clearness index is employed to represent the stochastic distribution of solar radiation intensity by beta distribution, while wind speed uncertainty is modeled seasonally using the Weibull distribution. Monte Carlo sampling with synchronous back substitution is applied for scenario generation and reduction of solar radiation and wind speed. To address the multi-objective evaluation, the analytic hierarchy process is utilized, and the joint optimization is achieved by combining a region contraction algorithm with stochastic programming. The proposed methodology is validated on an IES featuring various heating devices, incorporating uncertainties in both wind and solar energy. The results indicate that the absorption heat pump-based scheme achieves superior energy-saving performance, achieving an energy rate of 0.4724. Additionally, the compression heat pump-based scheme exhibits excellent economic efficiency and environmental sustainability, with a cost of energy of 0.3639 and a renewable fraction of 0.5536.