Modern power systems face a grand challenge in grid management due to increased electricity demand, imminent disturbances, and uncertainties associated with renewable generation, which can undermine grid security. The security assessment is directly connected to the robustness of the operating condition and is evaluated by analyzing proximity to the Power Flow (PF) solution space’s boundary. Calculating the location of such a boundary is a computationally challenging task, linked to the PF equations’ non-linear nature, presence of technical constraints, and complex network topology. This paper introduces a general framework to characterize points on the PF solution space boundary in terms of auxiliary variables subject to algebraic constraints. Then we develop an adaptive continuation algorithm to trace 1-dimensional sections of boundary curves which exhibits robust performance and computational tractability. Implementation of the algorithm is described in detail, and its performance is validated on different test networks.