Analysis of Grid-Connected Photovoltaic Generation Systems in the Harmonic Domain

: In this contribution, a harmonic domain (HD) model of grid-connected photovoltaic (PV) generation systems is proposed. The topology of the PV generation system contains a PV array connected to a DC / DC converter and a DC / AC converter, coupled to the grid through a ﬁlter. The individual components of the PV system are modeled in the HD. The PV array is represented as a Thevenin equivalent, whereas the DC / DC converter and the DC / AC converter are represented in the HD using switching function matrices. The periodic steady-state solution of the PV system is obtained through simple matrix / vector operations. The results obtained by the proposed model are validated against the response obtained with the PSCAD / EMTDC ® simulator, widely accepted by the power industry.


Introduction
Worldwide societies are increasingly concerned with the already powerful effects of climate warming. By the end of 2018 and over the course of 2019, atypical weather phenomena have taken place around the world, e.g., powerful hurricanes, droughts, floods, forest fires, and massive kelp arrival in coasts, all of them related to climate change due to the irrational consumption of fossil fuels and their polluting subproducts. Electrical energy generated in power plants using such fuels is one of the main causes of pollution and climate change around the world. Therefore, a fast transition to renewable energy sources is imperative, while the demand for electrical energy is continuously increasing to satisfy the multiple needs of growing societies and their associated industrial infrastructure.
Among the alternatives of renewable energy sources, photovoltaic (PV) technology is fast developing. It includes several technical aspects related to the dynamic operation, control, and power quality issues, among others. The above can be adequately assessed through the development and use of appropriate simulation models [1].
Successful examples of development of large PV systems have been reported [2,3] with their integration with the power network considered. An area of opportunity deals with the development of small PV energy systems and secondary energy sources, such as capacitor/battery banks and diesel generators [4][5][6]. PV systems are usually conveniently interfaced to the grid/load through converter arrangements; i.e., for their isolated operation [7,8] or grid-connected performance [9]. Efficiency and safety are important issues to be considered during the operation of PV systems. In particular, it has been observed that the efficiency may decline due to losses associated with the transformer operation and other auxiliary components. Alternative arrangements based on the elimination of transformers and isolated sensors may contribute to solve this problem and to reduce costs [10].
The reasons mentioned above have encouraged the modeling and analysis of electrical grids of generation with renewable energy sources under transient and/or periodic steady-state operation. This research is in particular centered on the periodic steady-state modeling and analysis of grids with PV sources in the harmonic domain (HD).
Traditionally, time-domain (TD) methods have been used to simulate grid-connected PV systems where a postprocessing stage is needed for the extraction of harmonics. This, in fact, is performed after the transient has died out and the periodic steady-state solution has been obtained. However, for lighted damped systems, achieving the periodic steady-state solution may require long simulation times [11]. Newton methods have been proposed to speed up the solution in time domain in electronic circuits [12] and in electrical networks [13]. The method described in [12] is referred to as the shooting method, and those proposed in [13] are referred to as Newton methods of extrapolation to the limit cycle based on Poincaré maps [14]. Their use and application are at present consolidated, based on the several applications reported for the efficient periodic steady-state solution of nonlinear and time-varying components and electrical networks [15]. In particular, in a recent contribution [16], the periodic steady state of a microgrid with PV generation is obtained in the time domain using a Newton method [13] and cubic splines.
As an alternative, the HD framework [17][18][19] is a consolidated powerful alternative for the determination of the periodic steady-state solution of nonlinear power networks. The harmonic domain concept was formally introduced in [17] to describe a general frame of reference for harmonic analysis with phases, harmonics, and harmonic cross-coupling, where linear and nonlinear components and effects of the power network were explicitly represented in a unified framework representation for their periodic steady-state solution of the entire electric network. Later references [18,19] reported the progress achieved in the harmonic domain in power network modeling and analysis. A more recent contribution [15] concisely details the procedure of linearization of nonlinear power components and the unified solution of linear and nonlinear networks in the harmonic domain. This is the approach followed in the research reported in this contribution. In recent contributions [20][21][22] the HD concept is used to report the modeling of aspects related to the connection of PV systems to the grid under harmonic distortion conditions. In another recent contribution in the frequency domain [23], the harmonics produced by the voltage ripple in a PV inverter with feedback control are analyzed. In [24], a model is proposed to assess the harmonic distortion produced from grid-connected PV inverters.
The paper is detailed in the sections to follow as: Section 2 describes the modeling of grid-connected PV generation system, Section 3 details the state-space TD representation of the total PV system, Section 4 presents the HD model, Section 5 shows a flowchart of the overall solution process, Section 6 shows a test-case solution with the proposed HD methodology and its validation against the simulator PSCAD/EMTDC ® , widely accepted by the power industry, and Section 7 gives the overall conclusions of this investigation.

System Configuration
In Figure 1, a single-phase grid-connected PV system is connected to a utility grid through a PV array, a capacitor link, a DC/DC and DC/AC converter, and a filter.
The decoupling AC-system dynamics from the PV array is achieved through the capacitor connected at the PV array output. The voltage level and maximum power point tracking (MPPT) of the PV array is achieved with the DC/DC converter; the DC/AC converter allows to obtain the AC power. The grid-connected PV inverters should satisfy the grid code, i.e., the network connection criteria, construction of the power generation system/network, system protection, and operation of the system [25]. A filter is used to mitigate the total harmonic distortion (THD) at the point of common coupling (PCC) with the PV system. The HD model of a single-phase PV system connected to a grid is described next.

PV System Representation
A PV device is an element that converts sunlight into electricity. A PV cell is the basic unit of a PV device; an arrangement of series-parallel connected cells form a module. Figure 2 shows the equivalent circuit of a PV module. The equation that mathematically describes the characteristic of a PV module is [26,27]: where and are the saturation and generated PV currents of the PV module, respectively, VT = / is the thermal voltage of the array with cells in series, is the Shockley diode current, being the second term in Equation (1), q is the electron charge (1.60217646 × 10 −23 J/K), T is the p-n junction temperature (Kelvin), a is the ideal diode constant, is the equivalent parallel resistance, and is the equivalent series resistance of the PV array. In this research work, a PV array is represented as a PV Thevenin equivalent, obtained with the calculation of the and parameters of Figure 2. The above is achieved by solving the nonlinear relation, Equation (1), with the PV array working at maximum power point (MPP). The PV Thevenin equivalent is calculated as follows: Step 1. In a PV array, and are given by the number of series and parallel cells. The total PV array characteristics are calculated as: = .

PV System Representation
A PV device is an element that converts sunlight into electricity. A PV cell is the basic unit of a PV device; an arrangement of series-parallel connected cells form a module. Figure 2 shows the equivalent circuit of a PV module. The equation that mathematically describes the characteristic of a PV module is [26,27]: where I 0 and I pv are the saturation and generated PV currents of the PV module, respectively, V T = N S kT/q is the thermal voltage of the array with N S cells in series, I d is the Shockley diode current, being the second term in Equation (1), q is the electron charge (1.60217646 × 10 −23 J/K), T is the p-n junction temperature (Kelvin), a is the ideal diode constant, R P is the equivalent parallel resistance, and R S is the equivalent series resistance of the PV array.

PV System Representation
A PV device is an element that converts sunlight into electricity. A PV cell is the basic unit of a PV device; an arrangement of series-parallel connected cells form a module. Figure 2 shows the equivalent circuit of a PV module. The equation that mathematically describes the characteristic of a PV module is [26,27]: where and are the saturation and generated PV currents of the PV module, respectively, VT = / is the thermal voltage of the array with cells in series, is the Shockley diode current, being the second term in Equation (1), q is the electron charge (1.60217646 × 10 −23 J/K), T is the p-n junction temperature (Kelvin), a is the ideal diode constant, is the equivalent parallel resistance, and is the equivalent series resistance of the PV array. In this research work, a PV array is represented as a PV Thevenin equivalent, obtained with the calculation of the and parameters of Figure 2. The above is achieved by solving the nonlinear relation, Equation (1), with the PV array working at maximum power point (MPP). The PV Thevenin equivalent is calculated as follows: Step 1. In a PV array, and are given by the number of series and parallel cells. The total PV array characteristics are calculated as: = .
(5) In this research work, a PV array is represented as a PV Thevenin equivalent, obtained with the calculation of the R P and R S parameters of Figure 2. The above is achieved by solving the nonlinear relation, Equation (1), with the PV array working at maximum power point (MPP).
The PV Thevenin equivalent is calculated as follows: Step 1. In a PV array, V OC and I SC are given by the number of series and parallel cells. The total PV array characteristics are calculated as: Step 2. It is assumed that I pv is equal to the maximum generated current I sc , and I 0 is calculated as: Step 3. With the PV array working at MPP, I MPP is given by: from the P-V curve of a PV module, the derivative of the power with respect to voltage is zero at MPP, i.e., by applying the Newton-Raphson (NR) method to the simultaneous solution of Equations (7) and (8), R S and R P can be obtained.
Step 4. With R S and R P known, and using the Thevenin theorem, the equivalent voltage source is calculated as: The equivalent resistance is: Equations (9) and (10) represent the resultant Thevenin equivalent to be interfaced.

Boost Converter
The boost converter (BC) is shown in Figure 3. It maintains the output voltage v out at the desired level, based on the input voltage v in and the load, which may fluctuate. The output voltage is higher than the input voltage.
Step 2. It is assumed that is equal to the maximum generated current , and is calculated as: Step 3. With the PV array working at MPP, IMPP is given by: from the P-V curve of a PV module, the derivative of the power with respect to voltage is zero at MPP, i.e., by applying the Newton-Raphson (NR) method to the simultaneous solution of Equations (7) and (8), and can be obtained. Step 4. With and known, and using the Thevenin theorem, the equivalent voltage source is calculated as: The equivalent resistance is: Equations (9) and (10) represent the resultant Thevenin equivalent to be interfaced.

Boost Converter
The boost converter (BC) is shown in Figure 3. It maintains the output voltage at the desired level, based on the input voltage and the load, which may fluctuate. The output voltage is higher than the input voltage. The diode is reverse biased when the switch is "on" and the output stage is isolated. The input supplies energy to the inductor. If the switch is "off", the output stage receives energy from the inductor and from the input. The switch is controlled by a pulse-width modulation (PWM) technique with a duty ratio d, determined as, which is adjusted as needed [26]. In Equation (8), TS is the BC switching period, being of the order of kHz [28]; is the time the switch remains in "on" state.

DC/AC Converter
The single-phase full-bridge inverter DC/AC shown in Figure 4 generates a sinusoidal AC output voltage. Switching of the inverter voltage is based on the PWM technique, where two sinusoidal waveforms and a high-frequency triangular carrier waveform are used to generate the The diode is reverse biased when the switch is "on" and the output stage is isolated. The input supplies energy to the inductor. If the switch is "off", the output stage receives energy from the inductor and from the input. The switch is controlled by a pulse-width modulation (PWM) technique with a duty ratio d, determined as, which is adjusted as needed [26]. In Equation (8), T S is the BC switching period, being of the order of kHz [28]; T on is the time the switch remains in "on" state.

DC/AC Converter
The single-phase full-bridge inverter DC/AC shown in Figure 4 generates a sinusoidal AC output voltage. Switching of the inverter voltage is based on the PWM technique, where two sinusoidal waveforms and a high-frequency triangular carrier waveform are used to generate the PWM signal [26]. The PWM signal is obtained through comparison of the triangular waveform V tri against V control and −V control as shown in Figure 5a. The reference voltage and carrier intersection through the peak value of V control (and −V control ) is shown in Figure 5b. The modulation ratio m a is calculated as: Energies 2019, 12, 4785 5 of 14 m a is also the ratio between the magnitude of the fundamental component of the output voltage V ac and the input voltage V dc . The frequency of V tri is on the order of several kHz and defines the order of the harmonics of V ac [27]. Harmonics in V ac appear as sidebands, each centered at twice the frequency modulation ratio m f , i.e., where fs is the switching frequency and f 0 is the AC-system frequency. For the unipolar PWM voltage switching in the voltage source inverter (VSI), the control of leg A is independent of the control for leg B. For instance, in leg A, when V control < V tri , S 1 is turned on, otherwise it is turned off (S 2 is complementary to S 1 ). The same applies for the control of leg B, but using −V control instead of V control .
Energies 2019, 12, x FOR PEER REVIEW 5 of 14 PWM signal [26]. The PWM signal is obtained through comparison of the triangular waveform against and − as shown in Figure 5a. The reference voltage and carrier intersection through the peak value of (and − ) is shown in Figure 5b. The modulation ratio is calculated as: is also the ratio between the magnitude of the fundamental component of the output voltage and the input voltage . The frequency of is on the order of several kHz and defines the order of the harmonics of [27]. Harmonics in appear as sidebands, each centered at twice the frequency modulation ratio , i.e., where fs is the switching frequency and f0 is the AC-system frequency. For the unipolar PWM voltage switching in the voltage source inverter (VSI), the control of leg A is independent of the control for leg B. For instance, in leg A, when < , is turned on, otherwise it is turned off ( is PWM signal [26]. The PWM signal is obtained through comparison of the triangular waveform against and − as shown in Figure 5a. The reference voltage and carrier intersection through the peak value of (and − ) is shown in Figure 5b. The modulation ratio is calculated as: is also the ratio between the magnitude of the fundamental component of the output voltage and the input voltage . The frequency of is on the order of several kHz and defines the order of the harmonics of [27]. Harmonics in appear as sidebands, each centered at twice the frequency modulation ratio , i.e., where fs is the switching frequency and f0 is the AC-system frequency. For the unipolar PWM voltage switching in the voltage source inverter (VSI), the control of leg A is independent of the control for leg B. For instance, in leg A, when < , is turned on, otherwise it is turned off ( is  Figure 6 shows the single-phase diagram of a grid-connected PV system. It consists of a voltage source, five transmission lines, a capacitor bank, two linear loads, and the PV array. The associated equivalent circuit is shown in Figure 7. The PV source has been replaced by the corresponding Thevenin equivalent, as explained in Section 2.  Figure 6 shows the single-phase diagram of a grid-connected PV system. It consists of a voltage source, five transmission lines, a capacitor bank, two linear loads, and the PV array. The associated equivalent circuit is shown in Figure 7. The PV source has been replaced by the corresponding Thevenin equivalent, as explained in Section 2.  Twelve differential algebraic equations (DAE), i.e., six differential equations and six algebraic equations represent the dynamic operation of the grid-connected PV system. The voltage at the capacitors and the currents in the inductors are selected as the state variables. The DAE set is detailed next.

State-Space Representation of the Grid-Connected PV Generation System
Differential equations: = .
(19) Figure 6. Grid-connected PV system. Figure 6 shows the single-phase diagram of a grid-connected PV system. It consists of a voltage source, five transmission lines, a capacitor bank, two linear loads, and the PV array. The associated equivalent circuit is shown in Figure 7. The PV source has been replaced by the corresponding Thevenin equivalent, as explained in Section 2.  Twelve differential algebraic equations (DAE), i.e., six differential equations and six algebraic equations represent the dynamic operation of the grid-connected PV system. The voltage at the capacitors and the currents in the inductors are selected as the state variables. The DAE set is detailed next.

State-Space Representation of the Grid-Connected PV Generation System
Differential equations: = . Twelve differential algebraic equations (DAE), i.e., six differential equations and six algebraic equations represent the dynamic operation of the grid-connected PV system. The voltage at the capacitors and the currents in the inductors are selected as the state variables. The DAE set is detailed next.

Harmonic Domain Model
The grid-connected PV system of Figure 7 can be solved in the TD with a state-space matrix equation representation based on nonautonomous ordinary differential equations of the form: where [A] is the state matrix of size n × n, [B] is the control or input matrix of size n × r, and u is the input vector of dimension r. Based on [6], the set of ordinary differential equations given by Equation (32) is conveniently converted into the HD following the procedure detailed in [18] resulting in: where: and: In Equation (35) the differentiation matrix D has been arranged in block-diagonal form to account for the twelve HD variables. Also, the switching functions S i = S 1 − S 2 and S 2 are represented in the HD as Toeplitz-type matrices S i and S 2 , respectively, and introduced into matrix [A] (Equation (37)). Note that in Equation (37) I d represents an identity matrix of appropriate dimensions.
Alternatively, for the HD solution of Equation (33), the steady-state solution of the PV system for a given operating point can be obtained as: In practical applications, the PV system must operate around the MPP. This can be achieved by an appropriate duty ratio of the BC. The switching function of the BC has been included in Equation (39) as a Toeplitz-type matrix for a given duty ratio. The solution of Equation (39) with any duty ratio leads to an operating point possibly different than the MPP. Then, an iterative solution scheme is needed. Figure 8 shows the proposed scheme. Generally speaking, Figure 8 shows that the HD system Equation (39) is solved for a given duty ratio and the MPP operating condition is assessed. In case that the response is negative, the duty ratio is increased or decreased based on a perturb and observe algorithm to achieve the MPP. In the flow diagram of Figure 8, V pv0 and I pv0 correspond to the DC components of the HD vectors V pv and I pv , respectively. Also, ∆d represents an increment (or decrement) of the duty ratio of S 1 . In our research we assumed d = 0.25 and a tolerance for convergence of 0.001 p.u.
Energies 2019, 12, x FOR PEER REVIEW 9 of 14 (39) as a Toeplitz-type matrix for a given duty ratio. The solution of Equation (39) with any duty ratio leads to an operating point possibly different than the MPP. Then, an iterative solution scheme is needed. Figure 8 shows the proposed scheme. Generally speaking, Figure 8 shows that the HD system Equation (39) is solved for a given duty ratio and the MPP operating condition is assessed. In case that the response is negative, the duty ratio is increased or decreased based on a perturb and observe algorithm to achieve the MPP. In the flow diagram of Figure 8, and correspond to the DC components of the HD vectors and , respectively. Also, Δd represents an increment (or decrement) of the duty ratio of . In our research we assumed d = 0.25 and a tolerance for convergence of 0.001 p.u.

Flowchart of HD Solution
In practical applications, the PV system must operate around the MPP. This can be achieved by an appropriate duty ratio of BC. In this research work, the perturb and observe algorithm is used to reach the MPP. Furthermore, the modeled PV system is a non-feedback system, i.e., the output of the system has no influence or effect on the control action of the input signal. In other words, it is an open-loop system. The switching function of the BC has been included in Equation (28) as a Toeplitztype matrix for a given duty ratio.
The flowchart of the solution process to obtain the periodic steady-state solution of gridconnected PV sources in the HD is shown in Figure 9. It has five blocks, where in the first block the parameters of the network, including the PV system, are read; in the second block, the Thevenin equivalent circuit of the PV array is calculated; in the third block, the set of ODEs is obtained from the circuit of Figure 7 and converted to the HD. In the fourth block the switching functions are represented in the HD as Toeplitz-type matrices and introduced into the matrix [A]. Finally, in the fifth block, the solution of the PV generation system for a given operating point is calculated. Through this procedure, the periodic steady-state solution of the entire grid is obtained via simple matrix/vector operations.

Flowchart of HD Solution
In practical applications, the PV system must operate around the MPP. This can be achieved by an appropriate duty ratio of BC. In this research work, the perturb and observe algorithm is used to reach the MPP. Furthermore, the modeled PV system is a non-feedback system, i.e., the output of the system has no influence or effect on the control action of the input signal. In other words, it is an open-loop system. The switching function of the BC has been included in Equation (28) as a Toeplitz-type matrix for a given duty ratio.
The flowchart of the solution process to obtain the periodic steady-state solution of grid-connected PV sources in the HD is shown in Figure 9. It has five blocks, where in the first block the parameters of the network, including the PV system, are read; in the second block, the Thevenin equivalent circuit of the PV array is calculated; in the third block, the set of ODEs is obtained from the circuit of Figure 7 and converted to the HD. In the fourth block the switching functions are represented in the HD as Toeplitz-type matrices and introduced into the matrix [A]. Finally, in the fifth block, the solution of the PV generation system for a given operating point is calculated. Through this procedure, the periodic steady-state solution of the entire grid is obtained via simple matrix/vector operations.

Simulation Results
The grid-connected PV system shown in Figure 7 is used as a test system. The corresponding parameters are given in Table A1 of Appendix A. Standard test conditions (STC) in the PV array are assumed, i.e., 1000 W/m 2 and 25 °C, with the PV system operating at MPP. Figure 10 shows the resultant TD waveforms and the corresponding harmonic content of the voltage at PCC. It is remarked that to obtain the TD waveforms presented in this section, Equation (28) is evaluated with a sampling time-step of 0.1 μs. The results of Figure 10 have been obtained by considering the explicit representation of 50 and 220 harmonics, respectively, in the HD model of the PV system. It can be noticed that higher order harmonics produced by power electronics components are not expected to appear beyond the twentieth harmonic. ®. Figure 9. Flowchart of the grid-connected PV system solution in the HD.

Simulation Results
The grid-connected PV system shown in Figure 7 is used as a test system. The corresponding parameters are given in Table A1 of Appendix A. Standard test conditions (STC) in the PV array are assumed, i.e., 1000 W/m 2 and 25 • C, with the PV system operating at MPP. Figure 10 shows the resultant TD waveforms and the corresponding harmonic content of the voltage at PCC. It is remarked that to obtain the TD waveforms presented in this section, Equation (28) is evaluated with a sampling time-step of 0.1 µs. The results of Figure 10 have been obtained by considering the explicit representation of 50 and 220 harmonics, respectively, in the HD model of the PV system. It can be noticed that higher order harmonics produced by power electronics components are not expected to appear beyond the twentieth harmonic.

Simulation Results
The grid-connected PV system shown in Figure 7 is used as a test system. The corresponding parameters are given in Table A1 of Appendix A. Standard test conditions (STC) in the PV array are assumed, i.e., 1000 W/m 2 and 25 °C, with the PV system operating at MPP. Figure 10 shows the resultant TD waveforms and the corresponding harmonic content of the voltage at PCC. It is remarked that to obtain the TD waveforms presented in this section, Equation (28) is evaluated with a sampling time-step of 0.1 μs. The results of Figure 10 have been obtained by considering the explicit representation of 50 and 220 harmonics, respectively, in the HD model of the PV system. It can be noticed that higher order harmonics produced by power electronics components are not expected to appear beyond the twentieth harmonic.  By comparing the individual harmonics obtained in the HD with the explicit representation of 50 harmonics, it is clear that all meet the IEEE 1547-2018 standard [25]. In the absence of a harmonic resonance, it is sufficient to explicitly represent in the HD solution the presence of the first 25 or at most 50 harmonics, thus dramatically reducing the required computational effort needed without losing accuracy in the HD solution.
The calculated THD and powers at PCC are listed in Table 1. The maximum error between responses obtained with the HD method and the PSCAD/EMTDC ® simulator is negligible, i.e., it is only 0.02% for the THD in the variable . The current waveform at PCC is illustrated in Figure 11a and the corresponding comparison in harmonic content is shown in Figure 11b. Again, a close agreement between responses can be observed. It is shown that the harmonic distortion of the waveform is mainly due to the presence of third, fifth, and seventh harmonics, these being of 5.12%, 0.15%, and 0.26%, respectively. The third harmonic (5.12%) exceeds the permissible limit set by the IEEE 1547-2018 standard [25], which is 4%.  By comparing the individual harmonics obtained in the HD with the explicit representation of 50 harmonics, it is clear that all meet the IEEE 1547-2018 standard [25]. In the absence of a harmonic resonance, it is sufficient to explicitly represent in the HD solution the presence of the first 25 or at most 50 harmonics, thus dramatically reducing the required computational effort needed without losing accuracy in the HD solution.
The calculated THD and powers at PCC are listed in Table 1. The maximum error between responses obtained with the HD method and the PSCAD/EMTDC ® simulator is negligible, i.e., it is only 0.02% for the THD in the variable i PCC . The current waveform at PCC is illustrated in Figure 11a and the corresponding comparison in harmonic content is shown in Figure 11b. Again, a close agreement between responses can be observed. It is shown that the harmonic distortion of the waveform is mainly due to the presence of third, fifth, and seventh harmonics, these being of 5.12%, 0.15%, and 0.26%, respectively. The third harmonic (5.12%) exceeds the permissible limit set by the IEEE 1547-2018 standard [25], which is 4%.

Conclusions
A model of a grid-connected PV generation system in the harmonic domain has been proposed. It has been evaluated for standard test conditions (STC). The proposed model can include any number of harmonics depending on the desired resolution. However, according to the obtained results, it should be noticed that in the absence of higher frequency resonances, using a smaller harmonic spectrum may represent a considerable reduction in the required computational resources to obtain the HD solution.
A major characteristic of the proposed model is that the periodic steady state can be obtained in a direct way through a solution process in the HD. Here, switching functions (of the BC and inverter) are represented in the HD as Toeplitz-type matrices.
This contribution is relevant since the developed method allows the determination of the periodic steady-state solution in networks with photovoltaic generation systems in the harmonic domain, which is a frame of reference that allows an accurate prediction of individual harmonics and total harmonic distortion (THD) of voltage and current through a unified solution where harmonics and harmonic cross-coupling effects are explicitly represented.
The reported results allow a precise assessment of individual harmonic magnitudes and THD, which indicates that the continuous operation of PV units in the electrical network is related with harmonics that may exceed permissible limits set by established standards.
The obtained results have been successfully validated against the response obtained with the PSCAD/EMTDC ® simulator, widely accepted by the power industry. A close agreement between the obtained responses has been obtained in all cases.

Acknowledgments:
The authors want to acknowledge CONACYT for finance assistance to conduct this investigation and to the Facultad de Ingeniería Eléctrica, División de Estudios de Posgrado, UMSNH, Morelia, México for the access and use of its facilities.

Conflicts of Interest:
The authors declare no conflict of interest.

Conclusions
A model of a grid-connected PV generation system in the harmonic domain has been proposed. It has been evaluated for standard test conditions (STC). The proposed model can include any number of harmonics depending on the desired resolution. However, according to the obtained results, it should be noticed that in the absence of higher frequency resonances, using a smaller harmonic spectrum may represent a considerable reduction in the required computational resources to obtain the HD solution.
A major characteristic of the proposed model is that the periodic steady state can be obtained in a direct way through a solution process in the HD. Here, switching functions (of the BC and inverter) are represented in the HD as Toeplitz-type matrices.
This contribution is relevant since the developed method allows the determination of the periodic steady-state solution in networks with photovoltaic generation systems in the harmonic domain, which is a frame of reference that allows an accurate prediction of individual harmonics and total harmonic distortion (THD) of voltage and current through a unified solution where harmonics and harmonic cross-coupling effects are explicitly represented.
The reported results allow a precise assessment of individual harmonic magnitudes and THD, which indicates that the continuous operation of PV units in the electrical network is related with harmonics that may exceed permissible limits set by established standards.
The obtained results have been successfully validated against the response obtained with the PSCAD/EMTDC ® simulator, widely accepted by the power industry. A close agreement between the obtained responses has been obtained in all cases.