A novel optimisation-based model of the power flow (PF) problem is proposed using complementarity conditions to properly represent generator bus voltage controls, including reactive power limits and voltage recovery processes. This model is then used to prove that the Newton–Raphson (NR) solution method for solving the PF problem is basically a step of the generalised reduced gradient algorithm applied to the proposed optimisation problem. To test the accuracy, flexibility and the numerical robustness of the proposed model, the IEEE 14-bus, 30-bus, 57-bus, 118-bus and 300-bus test systems and large real 1211-bus and 2975-bus systems are used, benchmarking the results of the proposed PF model against the standard NR method. It is shown that the proposed model yields adequate solutions, even in the case when the NR method fails to converge.