As a continuation of the work done by the QSE (Electrical Supply Quality) research group at the UPC (Polytechnic University of Catalonia) on harmonic load flow in electric power networks, this thesis aims to study existing harmonic load flow formulations, as well as the numerical resolution of the nonlinear equation systems derived from these formulations, in order to propose improvements for the former and compare performances of numerical methods for the latter. The improvements in the harmonic load flow formulations are related to a reduction in the number of iterations, for which an improved formulation is proposed. The comparison of numerical resolution methods is focused on analysing harmonic load flow formulation convergences and accuracies. The specific goals of the thesis are: (1) To propose an improved formulation for the harmonic load flow problem. This formulation should be applicable to electrical networks with highly distorted voltages. (2) To analyse the numerical resolution of all the considered harmonic load flow formulations (existing and improved) in terms of convergence and accuracy, by using a well-known numerical method (Newton-Raphson) and an alternative numerical method (Levenberg-Marquardt).