Power systems consist of generators, transformers, loads, and distributed power sources interconnected through lines. The reliable operation of the system requires that voltage and current magnitudes, angles, and power flow are within allowable ranges. Several methodologies have been proposed to solve the power flow consisting of nonlinear algebraic equations. Among them, the Newton-Raphson method is widely used due to its simplicity in building the bus admittance matrix, high accuracy, and fast convergence. However, the method has several limitations, such as singularity issues that occur when some or all Jacobian matrix elements become zero. This prevents the power flow calculation from converging since the inversion of the Jacobian matrix cannot be obtained. To address this limitation, this study proposes a novel and robust power flow calculation methodology that can solve singularities for various ill conditions using the Moore-Penrose Pseudo-inverse. This method improves the classical Newton-Raphson method by addressing its shortcomings. The performance of the proposed algorithm was evaluated using various testbeds, including balanced and unbalanced radial systems, a large-scale power grid, and systems with ungrounded transformers. The accuracy of the proposed algorithm was verified by comparing the power flow calculation with DIgSILENT Power Factory. The testbed consisted of modified IEEE 4-, 13-, 37-, and 69-bus systems.