## 1 Introduction  
Conventionally, electricity, gas and heating networks are designed and operated independently as individual systems. The increasing application of combined heat and power units (CHP), gas-fired units and other energy conversion facilities (e.g., power to gas (P2G) and power to heat (P2H) equipment) enhance the degree of interdependence of different energy systems [1; 2]. There is a growing interest in considering these multiple energy systems as an integrated whole, and the coordinated planning and operation of integrated energy systems (IES) to achieve the technical, economic, and environmental benefits from multiple energy utilization [3; 4; 5]. Steady-state modeling and analysis of an IES are usually considered to be the basis for its planning and operation. Many researchers have studied the steady-state analysis of the electricity-gas IES (i.e., IES composed of electricity and gas networks), the electricity-heat IES and the electricity-gas-heat IES. [PERSON] et al. [6] conducted a unified power- and gas-flow analysis of electricity-gas IES considering bi-directional energy conversion. [PERSON] et al. [7] proposed a decomposed method for power and gas flow analysis to explore the impacts of gas composition on both electricity and gas networks. Considering the influence of environmental temperature on the pipelines and the gas network, [PERSON] et al. [8] proposed viewing the gas temperature as a state variable in modeling power and gas flow problems. [PERSON] et al. [9; 10] demonstrated the necessity for combined analysis of electricity and heating networks. [PERSON] et al. [11] studied the physical interaction between electricity and heating networks. References [12; 13; 14] performed a steady-state energy-flow analysis of IEScontaining electricity, gas and heat based on the Newton-Raphson method. [PERSON] et al. [15] presented an integrated electricity-gas-heat transmission network model and applied it to the system in Great Britain. [PERSON] et al. [16] aimed for an urban EGH-IES and proposed an operational multi-objective optimization model. An integrated method for optimizing electricity, gas and heat in an EGH-IES was investigated in [17]. Although the references mentioned above have studied the steady-state analysis of electricity-gas IES [6; 7; 8], electricity-heat IES [9; 10; 11] and electricity-gas-heat IES [12; 13; 14; 15; 16; 17], they are all deterministic analyses, in essence, and cannot address the effects of the uncertainties existing in an IES (e.g., loads and renewable energies) on the operation of IES. Compared with conventional separated electricity, gas or heating networks, IES is exposed to more uncertain sources, such as intermittent renewable energy sources (wind, solar, etc.) and electricity/gas/heat loads. Moreover, uncertainties in IES are more complicated than in individual energy networks due to the mutual interaction among multiple energy networks. Any uncertain variation in one network may be transmitted to the other via coupling facilities (e.g., CHP unit, etc. ), which may finally deteriorate the overall performance of the IES and even pose risks to the security of each energy network. It is therefore becoming a kind of necessity to perform the analysis of IES with the uncertainties considered carefully. Probabilistic power flow analysis, which has been studied extensively, is a powerful tool for recognizing and evaluating the effects of uncertainties on electricity networks [18; 19]. It is therefore logical and desirable to extend the probabilistic power flow technique to the multi-energy flow analysis of IES. [PERSON] et al. [20; 21] presented a probabilistic power and gas flow model considering uncertainties in wind speeds and electricity/gas/heat loads, and then solved the problem by Monte Carlo simulation (MCS). [PERSON] et al. [22] applied the cumulant method and Gram-Charlier expansion to conduct probabilistic power and gas flow analysis. [PERSON] et al.