At the distribution network level, integrated electricity-gas systems (IEGSs) exhibit significant non-convexity due to the non-negligible impedance of power grids and the complexity of dynamic gas transmission. Moreover, considering the stochastic nature of disaster-inducing mechanisms of extreme events, resilience assessment requires load-shedding calculations under multiple fault scenarios, significantly increasing evaluation time. To accelerate the resilience assessment process, this chapter and the subsequent one propose a fast resilience assessment method for distribution-level IEGS that accounts for the non-convexity of dynamic gas flow. This chapter focuses on the optimal energy flow problem in distribution-level IEGS and introduces an efficient non-convex optimization algorithm. First, based on the physical characteristics of distribution-level IEGS, a non-convex optimal energy flow model is formulated, incorporating alternating current (AC) power flow equations for the electricity network and quasi-dynamic energy flow equations for the gas network. Second, an equivalent transformation is applied to standardize nonlinear and non-convex constraints of different mathematical forms into a standard bilinear equality constraint format. Next, within the framework of the spatial branch-and-bound (SB&B) algorithm, the concept, generation method, and proof process of the bilinear branching feasible region are proposed to guide branching decisions and accelerate algorithm convergence. Then, the SB&B algorithm is embedded using the callback mechanism of the CPLEX solver. Finally, by comparing the computational performance of 12 solution algorithms for mixed-integer bilinear optimization problems on two IEGS benchmark test cases of different scales, the accuracy and efficiency of the proposed non-convex optimization algorithm are validated.