Introduction

The presence of harmonics is unavoidable in a power system due to the nonlinear characteristics of equipment and connected loads [1]. Integration of solar PV system further intensifies the harmonics in the distribution network. These harmonics emanate within the inverter during the normal power conversion process [2]. The nonlinear behavior of solar PV with changing weather condition increases inverter switching, which injects additional harmonics in the system [3, 4]. The magnitude and order of harmonics depend on the power inverter technology used, adopted modulation technique, method of commutation, existence of high or low frequency coupling transformer, interconnection configuration and filters used at the output side of the inverter [5–9].

The interaction between grid components and number of solar PV units further amplifies the harmonic distortion [10].

The objective of this paper is to analyze the distribution network with solar PV integration, and to study the voltage/current harmonic distortions at various buses/branches. Paper also deals with the deviation of these distortions from the standard IEEE limits [11], thereby emphasizes the need for developing suitable control strategies to eliminate them in order to maintain the clean power supply. To examine these issues, the Newton-Raphson algorithm is reformulated and implemented by including the solar PV and nonlinear loads on the smart radial distribution feeder. The results of the harmonic power flow analysis are comprehensively analyzed and the observations are presented.

The paper is organized as follows: The modeling of various power system components in harmonic domain is presented in “Modeling of Power System Components”. The harmonic analysis and simulation algorithm is discussed in “Harmonic Analysis Algorithm”. Design of smart radial distribution grid is given in “Design of Smart Distribution Grid System”. “Simulation Results Analysis” presents simulation results analysis and associated discussion on the harmonic distortions. Finally, “Conclusion” concludes the paper by summarizing the work carried and investigations made.

Modeling of Power System Components

For analyzing the issues connected with the harmonics, we need to develop appropriate models of various power system components, used in our study. This section presents the harmonic domain modeling of various power system components.

Modeling of Solar PV

A PV system exhibits the nonlinear I-V characteristics and introduces harmonics in the system by itself. An ideal PV cell is modeled as a current source which represents the generated current due to incident light and a parallel diode for p-n junction of PV cell but the model of a real PV cell contains additional parallel resistance Rp, represents the current leakage through the cell and a series resistance Rs, accounts for an extra voltage drop between the junction as shown in Fig. 1a. In an ideal situation, to determine the upper limit of open circuit voltage and the efficiency of solar cell, it is assumed that only radiative recombination takes place. But In the actual solar cell, the recombination via impurities predominates as presented in Fig. 1b [12–15]. Furthermore, the presence of solar array capacitance, as presented in Fig. 2, also increases the voltage ripple and has the dominant effect on the performance of the switching regulators at higher switching frequencies [16, 17].

Fig. 1
figure 1

a Equivalent circuit, Homo junction solar cell. b Equivalent circuit, Solar cell with impurity recombination current

Full size image
Fig. 2
figure 2

Equivalent circuit with solar array capacitance

Full size image

The equation for current and voltage relationship of PV array is expressed in Eq. 1.

$$ I=I_{ph} -I_{0} \left[ \exp \left( \frac{q\left( V+R_{s} I \right)}{N_{s} KT} \right)-1 \right]-\left( \frac{V+R_{s} I}{R_{p} } \right) $$
(1)

Where, I ph : Photo generated current (A) i.e. current generated by the incident light (it is directly proportional to the sun irradiation), I d : Diode equation (A), I 0 : Reverse saturation or leakage current of the diode (A), q: Charge of electron (1.160217646 ×10−19 C), K: Boltzmann constant (1.3806503 ×10−23 J/K), T: Temperature of the p-n junction (∘Kelvin), N s : Number of cells, R s : Series resistance (Ω), R p : Parallel resistance (Ω), I d1: Recombination current.

Modeling of Transformer

The shunt branch is the main source of harmonics in transformer. These harmonics are produced due to nonlinear core loss resistance, nonlinear inductance and frequency dependent resistances of primary and secondary windings. Fig. 3 represents the transformer equivalent circuit containing phase shifter and electrical equivalents of harmonic producing elements. Moreover, the generated harmonics also depends on the phase shifting and transformer winding connections [18, 19].

Fig. 3
figure 3

Transformer equivalent circuit for harmonics studies

Full size image

Where,

R1, R2:

Primary and secondary winding resistance representing copper losses,

X1, X2:

Primary and secondary winding reactance representing leakage flux,

Lm:

Magnetizing reactance, representing core flux

Rm:

Magnetizing resistance representing Iron loss in the core

Rf1, Rf2:

Primary and secondary winding resistances depending on frequency, represents skin effect

N1: N2:

Ideal transformer turns ratio

The equivalent π harmonic model of transformer can be developed as shown in Fig. 4 [20–22].

Fig. 4
figure 4

Transformer π equivalent harmonic model

Full size image

Where, \(Y_{A}^{(h)},Y_{B}^{(h)}, Y_{C}^{(h)}\) are equivalent admittance of transformers after considering all the elements and harmonic producing branches. The harmonics admittance matrix for complete transformer is expressed in Eq. 2.

$$ [Y^{(h)}_{xy(tr)} ]=\left[\begin{array}{cc} (Y_{A}^{(h)} +Y_{B}^{(h)} ) & -Y_{A}^{(h)}\\ -Y_{A}^{(h)} & (Y_{C}^{(h)} +Y_{A}^{(h)} ) \end{array}\right] $$
(2)

Modeling of Distribution Lines and Cables

Mathematical modeling of power systems in frequency and time domain techniques employ an admittance matrix model of the system. For balanced system the admittance matrix is defined as [24, 25].

$$ \left[ Y_{\textit{bus}} \right] = \left[ A \right]. \left[ {Y_{\textit{prim}} } \right] \left[ A \right]^{T} $$
(3)

In Eq. 3, Incidence matrix [A] is used to represent network connectivity.

An un-transposed, LV distribution feeder feeds the loads which are highly nonlinear and contains harmonics. In order to analyze these conditions, a harmonic model of the power system network is required [23]. For non-sinusoidal operating conditions the admittance matrix will be defined at harmonic frequencies. One such LV three phase, four wire system is shown in Fig. 5, in which the conductors are connected between the buses x and y. The size of the harmonics admittance matrix for this system will be of 4 × 4 (4), consisting of self and transfer admittances between phases and neutral [26, 27, 29].

$$ [Y^{(h)}_{xy}] = \left[\begin{array}{cccc} y_{aa}^{(h)} & y_{ab}^{(h)} & y_{ac}^{(h)} & y_{an}^{(h)} \\ y_{ba}^{(h)} & y_{bb}^{(h)} & y_{bc}^{(h)} & y_{bn}^{(h)} \\ y_{ca}^{(h)} & y_{cb}^{(h)} & y_{cc}^{(h)} & y_{cn}^{(h)} \\ y_{na}^{(h)} & y_{nb}^{(h)} & y_{nc}^{(h)} & y_{nn}^{(h)} \end{array} \right] $$
(4)
Fig. 5
figure 5

Representation of three phase four wire distribution line

Full size image

Modeling of Load

The harmonic power flow analysis requires equivalent circuits of linear loads because they affect the resonance conditions and provide damping at higher frequencies [30]. At fundamental frequency, linear loads are modeled as conventional PQ and PV buses. However, shunt admittances are used to model them at harmonic frequencies [28]. The admittance of a linear load connected to bus y at the h th harmonic frequency is given in Eq. 5.

$$ y_{y}^{(h)} =\left( \frac{P_{y}^{(1)} -j\frac{Q_{y}^{(1)}}{h}}{\left| {V_{y}^{(1)} } \right|^{2}} \right) $$
(5)

At the harmonic frequencies, the nonlinear loads (such as line commutated converters, arc furnace etc.) can be modeled as harmonic current sources and shunt admittance. The harmonic current at any bus y will be function of its fundamental and harmonic voltages as expressed in Eq. 6.

$$ I^{\left( h \right)}=f_{y}^{h} \left( V_{y}^{(1)}, V_{y}^{(5)} ,V_{y}^{(7)}, V_{y}^{(9)} ,.....V_{y}^{(L)} ,\alpha_{y}, \beta_{y}\right) $$
(6)

Where, L is the maximum harmonic order, α k and β k are nonlinear load control parameters.

Harmonic Analysis Algorithm

The harmonic power flow is a reformulation of the conventional Newton-Raphson power flow algorithm to include non-linear loads [31]. In this section, the conventional Newton-Raphson algorithm is reformulated by including the solar PV and nonlinear loads. The solar PV can be modeled as a constant PQ model and negative load with current injecting into the node. For steady-state analysis, constant negative PQ model of DG resources can be used [32]. This method checks the power mismatch at the end of feeders and laterals, in order to find the harmonics and system operating variables of PV integrated distribution system [33, 34]. Initially, the variables are arranged into a real vector Z, which consists of fundamental, harmonic voltages and nonlinear solar PV elements.

$$ [Z]=\left[{V^{(1)}, V^{(5)},V^{(7)}.......V^{(L)},[{\Omega} ]} \right]^{T} $$
(7)

In Eq. 7, [Ω] are the nonlinear parameters of solar PV and L is the maximum harmonic order. The algorithm steps for power flow solution are given below:

Step 1::

Read the distribution networks line data and bus data. Implement the fundamental power flow analysis and compute an initial value for the fundamental bus voltage magnitudes and phase angles.

Step 2::

Calculate the fundamental and harmonic current mismatch vector for solar PV and various load buses in Fig. 6.

$$ {\Delta} I^{(1)}=\left[ \begin{array}{l} I^{(1)}_{r,1} ,I^{(1)}_{i,1} ,.........I^{(1)}_{r,x-1},I^{(1)}_{i,x-1} , \\ I^{(1)}_{r,x} +g^{(1)}_{r,x} ,I^{(1)}_{i,x} +g^{(1)}_{i,x} ,... \\ I^{(1)}_{r,y} +g^{(1)}_{r,y} ,I^{(1)}_{i,y} +g^{(1)}_{i,y} \end{array} \right]^{T} $$
(8)
$$ {\Delta} I^{(h)}=\left[ \begin{array}{l} I^{(h)}_{r,1} ,I^{(h)}_{i,1} ,.........I^{(h)}_{r,x-1} ,I^{(h)}_{i,x-1} , \\ I^{(h)}_{r,x} +g^{(h)}_{r,x} ,I^{(h)}_{i,x} +g^{(h)}_{i,x} , \\ I^{(h)}_{r,y} +g^{(h)}_{r,y} ,I^{(h)}_{i,y} +g^{(h)}_{i,y} \end{array} \right]^{T} $$
(9)

In Eq. 8, \(I^{(1)}_{r,x} \) and \(I^{(1)}_{i,x} \) are the real and imaginary bus injection currents at bus x. Superscript (1) is used for fundamental frequency. \(g^{(1)}_{r,x} \) and \(g^{(1)}_{i,x} \) are real and imaginary line currents. Subscript (y) is used for corresponding bus injection current at y bus. In Eq. 9, \(I^{(h)}_{r,x}\) and \(I^{(h)}_{i,x}\) are the real and imaginary bus injection currents at bus x. Superscript (h) is used for harmonic frequency. \(g^{(h)}_{r,x}\) and \(g^{(h)}_{i,x}\) are real and imaginary line currents.

Step 3::

Evaluate the mismatch vector for mismatch in power and current on each bus. Check the convergence.

$$ {\Delta} M=[{\Delta} W,{\Delta} I^{(5)},.......{\Delta} I^{(L)},{\Delta} I^{(1)}]^{T} $$
(10)

In Eq. 10, ΔW is the mismatch power vector and (ΔI (5),.......ΔI (L),ΔI (1)) is mismatch current vector for harmonics including fundamental.

Step 4::

Evaluate jacobian matrix and calculate correction vector.

$$ J^{(h)}=\left[\begin{array}{cccc} \frac{\partial {\Delta} P_{2}^{(h)}}{\partial \delta_{2}^{(h)}} & \frac{\partial {\Delta} P_{2}^{(h)}}{\partial V_{2}^{(h)}}...... & \frac{\partial {\Delta} P_{2}^{(h)}}{\partial \delta_{y}^{(h)}} & \frac{\partial {\Delta} P_{2}^{(h)}}{\partial V_{y}^{(h)}} \\ \frac{\partial {\Delta} Q_{2}^{(h)}}{\partial \delta_{2}^{(h)}} & \frac{\partial {\Delta} Q_{2}^{(h)}}{\partial V_{2}^{(h)}}...... & \frac{\partial {\Delta} Q_{2}^{(h)}}{\partial \delta_{y}^{(h)}} & \frac{\partial {\Delta} Q_{2}^{(h)}}{\partial V_{y}^{(h)}} \\ {..} & {..} & {..} & {..} \\ {..} & {..} & {..} & {..} \end{array} \right] $$
(11)

In Eq. 11, J (h) is the jacobian matrix at harmonic frequency which has all entries zero for linear buses. ΔP and ΔQ are the change in real and reactive powers, w.r.t voltage V and phase angle δ.

Step 5::

Compute the bus voltage correction vector as given in Eq. 12. η is the number of iterations.

$$ {\Delta} Z=Z^{(\eta )}-Z^{(\eta +1)} $$
(12)
Step 6::

Update the vector Zas given in Eq. 13.

$$ Z^{(\eta +1)}=Z^{(\eta )}-{\Delta} Z $$
(13)
Step 7::

Update the real and reactive power mismatch due to solar PV and nonlinear loads at given buses.

$$ {\Delta} P_{k}^{\textit{nonlinear}} =P_{k}^{(t)} +\sum\limits_{h} {F_{r,k}^{(h)}} $$
(14)
$$ {\Delta} Q_{k}^{nonlinear} =Q_{k}^{(t)} +\sum\limits_{h} {F_{i,k}^{(h)} } $$
(15)

Where, k = 1,2,3,4,......x, y.

In Eqs. 14–15, P k and Q k are the total real and reactive power at bus k.F r, k and F i, k are the total real and imaginary line powers.

Step 8::

Go to step 2. Calculate the mismatch vectors again and check the convergence.

Fig. 6
figure 6

Distribution grid with solar PV integration

Full size image

Design of Smart Distribution Grid System

The solar PV integrated smart distribution grid system is modeled in Fig. 6. Loads are fed from an 11KV grid through a transformer of rating 1.5 MVA, 11/0.433 KV. Laterals are feeding the loads of varying nature in the area. Table 1 presents the parameters of solar PV which are connected at bus 5 and 16 respectively.

Table 1 Solar PV parameters
Full size table

Simulation Results Analysis

This section presents the extensive analysis of harmonics with and without solar PV integration at the various buses and branches of the modeled smart distribution grid system. The branch wise comparison of total harmonic distortions in current, with and without solar PV integration, has been presented in Table 2.

Table 2 Current THD in various branches
Full size table

It is found that the harmonic distortion in current rises sharply after solar PV integration. The current harmonic distortion (%THD) in the branch connected between buses 2–3 increases to a very high value of 50.18 %. The nature of such distorted current is shown in Fig. 7. Analysis of harmonic amplitude spectrum, as shown in Fig. 8 depicts that it contains the 5 th, 7 th, 11 th, 13 th, 17 th, and 19 th order harmonics, which violates the acceptable limits as given in Table 6. The share of 5 th and 7 th order harmonic is more, which affects the system performance severely and needs immediate filtering.

Fig. 7
figure 7

Harmonic current waveform in the branch connected between buses 2 and 3

Full size image
Fig. 8
figure 8

Spectrum of harmonic current in the branch connected between buses 2 and 3

Full size image

Harmonics makes the distribution transformer current highly distorted as shown in Fig. 9. Amplitude spectrum in Fig. 10 contains the 5 th, 7 th, 11 th, 13 th, 17 th and 19 th order harmonics which violates the standard limits given in Table 6. Such current escalate the iron losses and reduces the transformer efficiency, consequently derating becomes essential.

Fig. 9
figure 9

Transformer harmonic current waveform

Full size image
Fig. 10
figure 10

Spectrum of transformer harmonic current

Full size image

The bus wise comparison of total harmonic distortions in voltage, with and without solar PV integration, has been presented in Table 3. Results show that with solar PV integration, the THD at bus 16 rises to a very high value equal to 13.67 % as compared to the THD value of 1.28 % without solar PV, Moreover, all other buses are also affected by harmonics. A clear violation of 5 % voltage standard limits, as given in Table 7 can be seen after solar PV integration.

Table 3 Voltage THD at various buses
Full size table

The distorted voltage waveform at bus 16 is presented in Fig. 11. The individual voltage harmonic contribution is depicted in concerned harmonic spectrum Fig. 12. It is found that, the harmonic spectrum contains 5 th, 7 th, 13 th, 17 th, and 19 th order harmonics which violates the individual harmonic distortion limits of 3 %. Moreover, these harmonic distortions alter the value of fundamental and RMS voltages as presented in Table 4. Therefore, a proper mitigation technique is essential to filter out these harmonics and resolve the related issues.

Fig. 11
figure 11

Voltage waveform at bus 16

Full size image
Fig. 12
figure 12

Spectrum of harmonic voltage at bus 16

Full size image
Table 4 Variation in voltage due to harmonics after solar PV integration
Full size table

Table 5 presents the comparative results of harmonic power flow variables. The variation in current is very high at buses 5 and 16 after solar PV integration, which poses burden on the connecting cables and requires higher size cables. Moreover, negative variation in power flow and current is observed at some buses. Such reverse power can affect the protection system and false tripping of protective gears can occur (Tables 6 and 7).

Table 5 System power flow variations after solar PV integration
Full size table
Table 6 Limits of maximum harmonic current distortion inPercent, (Based on IEEE Std. 519-1992, 1993) [11]
Full size table
Table 7 Limits of harmonic voltage distortion at the PCC (Based on IEEE Std. 519-1992, 1993) [11]
Full size table

Conclusion

The power distribution system has a large number of nonlinear loads, which generates harmonics in the system. Integration of solar PV further enhances the distortion and thereby degrades the quality of power supply. Mitigation of harmonics is possible only if their exact nature, individual percentage and THD are known. Normal power flow method does not provide above information. In this paper, harmonic power flow approach is implemented for a solar PV integrated smart distribution system. Comparison results of before and after solar PV integration has been analyzed. It is observed that integration of solar PV alters the system power flow variables and increases the harmonics in the system beyond the prescribed limit, which is a major challenge and needs to be addressed in the integration process. The grid must be able to absorb these challenges in order to achieve the smart grid objectives. The analysis presented in this paper will be helpful for designing an electronic filter as a control solution and make the system harmonic free.