Power distribution systems are experiencing a largescale integration of Converter-Interfaced Distributed Energy Resources (CIDERs). This complicates the analysis and mitigation of harmonics, whose creation and propagation are facilitated by the interactions of converters and their controllers through the grid. In this paper, a method for the calculation of the so-called Harmonic Power-Flow (HPF) in three-phase grids with CIDERs is proposed. The distinguishing feature of this HPF method is the generic and modular representation of the system components. Notably, as opposed to most of the existing approaches, the coupling between harmonics is explicitly considered. The HPF problem is formulated by combining the hybrid nodal equations of the grid with the closed-loop transfer functions of the CIDERs, and solved using the Newton-Raphson method. The grid components are characterized by compound electrical parameters, which allow to represent both transposed or non-transposed lines. The CIDERs are represented by modular linear time-periodic systems, which allows to treat both grid-forming and gridfollowing control laws. The method’s accuracy and computational efficiency are confirmed via time-domain simulations of the CIGRÉ low-voltage benchmark microgrid. This paper is divided in two parts, which focus on the development (Part I) and the validation (Part II) of the proposed method.