AbstractA new probabilistic method is presented for the harmonic load‐flow analysis in electrical power systems. In the method, the network is modelled using the bus impedance matrix at each harmonic order. The load flow is solved by a current injection method. For harmonic sources, the amplitude variation and the phase angle variation are taken into account. For the summation of the harmonics originating from different sources, a probability density function, derived from the jointly normal distribution is used. As a final result, statistical estimates are obtained for the node voltages and branch currents of each harmonic frequency considered.