Integrated energy system applications can significantly improve energy efficiency. In this paper

we establish an integrated energy system containing heat

electricity and gas. The existing power flow (PF) calculation method applied to integrated energy systems (IESs) does not consider non-smooth constraints

such as the piecewise pipeline friction coefficient and generator buses reactive power limits

etc. Mixed integer nonlinear programming (MINLP) is conventionally used to deal with piecewise pipeline friction coefficients in gas network parts

but it is both complex and inefficient. Hence

we develop a piecewise linear function-based fitting method that can reduce the number of integer variables and enhanced the computational efficiency. In the electric network part

if the reactive power of the PV bus violates limits

it will be converted into a PQ bus

which is a non-differentiable and non-smooth constraint. Mixed complementarity problems are conventionally introduced to represent the PV-PQ buses type switching relationship and are addressed by the Newton-Raphson (NR) method. However

the above method is sensitive to the initial point. Here

we introduce a robust projected Levenberg-Marquardt (PLM) algorithm to cope with this issue. We demonstrate the advantages of our method and validate it both in a small-scale system and large-scale network test cases.