There are a large number of limiting non-smooth constraints in multi-region interconnected power grids. These constraints increase the non-convexity of the power flow model. Given that the traditional distributed power flow algorithm, based on a heuristic rule processing limit, is prone to convergence problems in the calculation process, a fully distributed power flow calculation method for an interconnected power grid that can robustly handle non-smooth constraints is proposed. First, the grid is partitioned according to the hierarchical and partitioned dispatch mode, and the non-smooth constraints in the model are smoothly processed. Then, based on the bi-level alternate direction inexact Newton method with a second-order convergence rate, the power flow problem is transformed into a problem of solving the optimal step increment. Based on the null space method, the coefficient matrix of the state variable is reduced in dimension, and the conjugate gradient algorithm is used to update the dual multiplier of each partition, and the second-order information is used in the multiplier update process to improve the convergence of the algorithm. The coordination layer does not need to be involved in distributed computing between multiple regions, and only a small amount of boundary information needs to be communicated, so the communication burden is light. Finally, the 30-bus system and 182-bus system are used as test examples to verify that the proposed method has higher accuracy and better convergence when setting bad initial values and dealing with non-smooth constraints. This work is supported by the National Natural Science Foundation of China (No. 52177125 and No. 51707196).