## 1 Introduction  
The increasing demand for energy has resulted in increased concerns regarding energy sources, and the exploration of new energy supply modes has attracted considerable interest. The concept of integrated energy systems (IESs) challenges that of traditional energy systems. IESs, which can fully utilize energy in various forms and optimize resource allocation, have rapidly developed in recent years [1]. Compared with traditional energy systems, these systems have significant advantages [2; 3]. For example, all types of energy are closely coupled, which makes them complementary and mutually beneficial, and they can realize cascade utilization and collaborative optimization of energy. Presently, research on IESs is mostly focused on three aspects: energy flow analysis and calculation [4], multi-energy system optimization scheduling [2; 3], and multi-energy system planning and multi-market games [5; 6]. The ability to integrate various types of energy structures is a compelling advantage of IESs [7; 8]. However, owing to the significant amount of renewable energy sources (RESs) involved, crucial challenges in terms of energy uncertainties arise. The uncertainties come from the intermittency of RESs; additionally, the volatility of the loads in different energy structures aggravate the energy uncertainties in IESs [9], which creates significant challenges for the operation and maintenance of IESs. We can determine the probabilistic distribution of each node state in the IESs by calculating the probabilistic energy flow, which provides guidance to address possible risks. Therefore, the probabilistic energy flow (PEF) in IESs must be studied. The basis of PEF is the calculation of steady-state energy flow (SSEF), which has been studied extensively. Because of the similarity between the SSEF model of the power system and the gas network, [PERSON]'s method can be used to solve the SSEF of the gas network. When the initial point of [PERSON]'s method is near the equilibrium point, the method has good and fast convergence. Moreover, the IES mostly works near the equilibrium point. Therefore, [PERSON]'s method is a commonly used method for solving the energy flow problem. Moreover, the solutions obtained using this method for multiple energy systems can be divided into two forms according to the Jacobian matrices used in the solution, i.e., separate Newton iterations for each system [10] or a unified Newton iteration for all systems [11]. In [10], each energy system (e.g., power system or natural gas system) is represented by its own Jacobian matrix and independent iterations are performed using [PERSON]'s method. The unbalanced quantity is iterated using the energy coupling unit until the whole energy system achieves convergence. This method is relatively simple, but the computation speed is somewhat slow. In [11], the power and natural gas networks are represented by a single Jacobian matrix for the unified Newton iterations. Then, the integrated energy flow of the entire network can be obtained. This method is more complicated in cases where there is a compressor in the gas grid, but the iterations are faster. Presently, research on the impact of uncertainties in IESs is limited. However, the calculations performed for these systems are similar to the probabilistic power flow (PPF) calculations in traditional power systems. They include the Monte Carlo method [12], analytical method [13], stochastic spectral analysis method [14], and approximation method [15]. Generally, a Monte Carlo simulation (MCS) requires several samples to obtain good precision according to the law of large numbers; this simulation requires considerable time to compute. However, it is possible to achieve the most accurate results with an MCS; therefore, it is used as reference to measure the precision of other methods. The studies on MCSs mainly focus on improving sampling efficiency to reduce the number of sampling representatives, for instance stratified sampling [16; 17] and importance sampling [18].