Three-phase Harmonic Analysis Method for Unbalanced Distribution Systems

Due to the unbalanced features of distribution systems, a three-phase harmonic analysis method is essential to accurately analyze the harmonic impact on distribution systems. Moreover, harmonic analysis is the basic tool for harmonic filter design and harmonic resonance mitigation; therefore, the computational performance should also be efficient. An accurate and efficient three-phase harmonic analysis method for unbalanced distribution systems is proposed in this paper. The variations of bus voltages, bus current injections and branch currents affected by harmonic current injections can be analyzed by two relationship matrices developed from the topological characteristics of distribution systems. Some useful formulas are then derived to solve the three-phase harmonic propagation problem. After the harmonic propagation for each harmonic order is calculated, the total harmonic distortion (THD) for bus voltages can be calculated accordingly. The proposed method has better computational performance, since the time-consuming full admittance matrix inverse employed by the commonly-used harmonic analysis methods is not necessary in the solution procedure. In addition, the proposed method can provide novel viewpoints in calculating the branch currents and bus voltages under harmonic pollution which are vital for harmonic filter design. Test results demonstrate the effectiveness and efficiency of the proposed method.

Traditional harmonic analysis methods are developed based on single-phase representation, and they cannot be used straightforwardly in distribution systems due to the topological characteristics of distribution systems.The well-known topological characteristics of distribution systems include radial network structure, unbalanced loads, multi-phase and unbalanced operation, large number of branches and nodes, wide-ranging resistance and reactance values, etc. Owing to the unbalanced line sections and loads in distribution systems, a three-phase harmonic analysis method is essential to accurately analyze the harmonic impact on distribution systems [7][8][9][10][11][12].The harmonic analysis is also the basic tool for harmonic filter design and harmonic resonance mitigation and has to be done repeatedly for each harmonic order [13][14][15][16]; therefore, the accuracy and efficiency are both important issues for a well-designed harmonic analysis method.Conventional harmonic analysis methods [7][8][9][10][11] don't consider the special topological characteristics of distribution systems such as the unbalanced loads and radial network structure; therefore, the solution accuracy cannot be further improved and computational time cannot be effectively reduced.A harmonic analysis method was proposed based on frequency-scan formulation, exact three-phase models, and commonly used forward/backward sweep technique that can fully use the topological characteristics of distribution systems [12].However, some specific data formats, such as the coefficient vector defined to record whether the system harmonic currents flow through the branch, needed to be designed and integrated into the commonly used forward/backward sweep technique.Therefore, the method proposed in [12] cannot be used straightforwardly to solve other harmonic problems, i.e., optimal harmonic filter design, etc.
An accurate and efficient three-phase harmonic analysis method for unbalanced distribution systems is proposed in this paper.The variations of bus voltages, bus current injections and branch currents affected by harmonic currents of nonlinear loads and shunt equipments can be analyzed by two relationship matrices developed from the topological characteristics of distribution systems.Therefore, the topological characteristics of distribution systems can be completely integrated into the proposed method for computational performance improvement.Some useful formulas are derived to solve the three-phase harmonic propagation problem.After the harmonic propagation for each harmonic order was calculated, the total harmonic distortion (THD) for bus voltages could be calculated accordingly.Exact three-phase models are used for the proposed harmonic analysis method; therefore, accurate solutions can be obtained.The proposed method is also efficient, since the time-consuming full admittance matrix inverse employed by commonly-used harmonic analysis methods is not necessary in the solution procedure.Compared with the method proposed in [12], no specific data format needs to be designed and integrated in the proposed solution procedure.Moreover, the proposed method can provide novel viewpoints in calculating the branch currents and bus voltages under harmonic pollution.The harmonic currents aborted by the shunt equipments such as harmonic filters can also be easily quantified in the proposed solution procedure.This data is vital for harmonic filter design and harmonic resonance mitigation.Test results demonstrate the validity of the proposed method.

Commonly-Used Harmonic Analysis Method
Due to the unbalanced line sections and loads in distribution systems, three-phase models should be used in harmonic analysis.For example, the impedance matrix of a three-phase line section ([Z abc ]) can be expressed as: The commonly-used harmonic analysis method [1,6] is developed based on the full admittance matrix inverse of a distribution system and can be obtained by: where [Y (h) ] is the full admittance matrix of a distribution system for the h th -order harmonic; [V (h) ] and [I (h) ] are the vectors of bus voltage and current injection for the h th -order harmonic, respectively.Equation ( 2) can be solved by triangular factorization and forward/backward substitution or by an inverse matrix.The procedures are time-consuming for large-scale unbalanced distribution systems.Especially, a harmonic propagation calculation should be done repeatedly for each harmonic order.If the topological characteristics of distribution systems can be taken into account to avoid the triangular factorization of a full admittance matrix, computational time can be significantly reduced.The voltage THD after the harmonic propagation of each order has been solved can be represented as: (%) 100% where and VTHD i (%)are the h th -order harmonic voltage and voltage THD percentage for bus i, respectively.

Proposed Relationship Matrices
The proposed method is developed based on two relationship matrices, the bus-current-injection to branch-current (BIBC) matrix and branch-current to bus-voltage (BCBV) matrix.Detailed derivations of these two matrices can be found in [17,18].In this paper, only the building algorithms of these two matrices are shown.The relationship between bus current injections and branch currents for a radial distribution system can be expressed as follows: where [BIBC] is the bus-current-injection to branch-current matrix for a radial distribution system; [I] and [B] are the vectors for bus current injections and branch currents, respectively.The constant BIBC matrix is an upper triangular matrix and has non-zero entries of +1 only.The building algorithm of BIBC matrix from substation to downstream sequentially is: Procedure 1 For a distribution system with m-branch sections and n-bus, the dimension of the BIBC matrix is m × (n − 1); Procedure 2 If a line section (B k ) is located between bus i and bus j, copy the column of the i-th bus of the BIBC matrix to the column of the j-th bus and fill a +1 to the position of the k-th row and j-th bus column; Procedure 3 Repeat Procedure 2 until all line sections are included in the BIBC matrix.
The relationship between branch currents and bus voltages for a radial distribution system can be written as: where [BCBV] is the branch-current to bus-voltage matrix; [V 0 ] and [V] are the vectors of no-load bus voltages and bus voltages for distribution systems, respectively.For a distribution system with N buses, the dimension for [V 0 ] and [V] is (N − 1) by 1.
The constant BCBV matrix has non-zero entries consisted of complex line impedance values.The building algorithm of BCBV matrix from substation to downstream sequentially is: Procedure 4 For a distribution system with m-branch sections and n-bus, the dimension of the BCBV matrix is (n − 1) × m; Procedure 5 If a line section (B k ) is located between bus i and bus j, copy the row of the i-th bus of the BCBV matrix to the row of the j-th bus and fill the line impedance (Z ij ) to the position of the j-th bus row and k-th column; Procedure 6 Repeat Procedure 5 until all line sections are included in the BCBV matrix.
Note that if the line section between bus i and bus j as represented in Procedures 1-6 is a three-phase line section, the corresponding branch current B k will be a 3 × 1 vector, the +1 s in the BIBC matrix become a 3 × 3 identity matrix, and the Z ij in the BCBV matrix constitute a 3 × 3 impedance matrix.Equations ( 4) and ( 5) are very useful for derivations of the proposed harmonic analysis method.Detailed formulas will be derived in the next section.

Three-Phase Harmonic Analysis Method
Equations ( 4) and ( 5) provide a new perspective in observing the relationship between bus voltages, bus current injections and branch currents.As can be seen in Equation ( 4), the mathematical equation between bus current injections and branch currents, can be used to calculate the variations of the branch currents affected by the harmonic current injections.In Equation ( 5), the mathematical equation between the branch currents and bus voltages, can be used to calculate the variations of the bus voltages resulting from the harmonic branch currents.Figure 1 illustrates the concepts of the proposed harmonic propagation.The harmonic currents flowing in a distribution system can be divided into two categories: the harmonic currents contributed by nonlinear loads, Is (h) as denoted in Figure 1, and the harmonic currents absorbed by shunt equipments, ( ) h If as denoted in Figure 1, such as harmonic filters.
In general, harmonic currents contributed by nonlinear loads can be obtained by field measurements and/or mathematical formulas.Some research aimed to find the mathematic formulas for the harmonic currents contributed by nonlinear loads [7,19].Using the derived mathematical formulas, the harmonic impacts of nonlinear loads can be conducted only by simulations; therefore, the cost and time for power quality analyzer installation and measurement can be saved.The harmonic currents obtained by mathematical formulas can also be used in the proposed harmonic analysis method.The harmonic currents absorbed by shunt equipments need to be calculated according to the network topology and harmonic currents of nonlinear loads.If the harmonic currents absorbed by shunt equipments were calculated, the branch currents and bus voltages affected by these two types of harmonic currents could then be analyzed by Equations ( 4) and ( 5), respectively.The harmonic currents of nonlinear loads can be expressed in a vector form as: where [Is (h) ] is the h th -order harmonic current vector of nonlinear loads.The branch currents affected by the h th -order harmonic currents of nonlinear loads can be calculated by substituting Equation (6a) into Equation ( 4) and be expressed as: where [Bs (h) ] is the branch current vector affected by the h th -order harmonic currents of nonlinear loads.
Similarly, the h th -order harmonic current vectors for shunt equipments, [If (h) ], and branch currents affected by the h th -order harmonic currents of shunt equipments, [Bf (h) ], can be expressed as: The aggregate branch current vector affected by the h th -order harmonic currents of nonlinear loads and shunt equipments, [B (h) ], can be written as: Equation ( 8) can be rewritten as: ( ); ( ) x Isbus g y Ifbus h where ( ) Isbus ⋅ and ( ) Ifbus ⋅ are the buses relating to the nonlinear loads and shunt equipments, respectively.Ns and Nf are the numbers of nonlinear loads and shunt equipments, respectively.
[BIBC z ] is the column vector of [BIBC] relating to the bus of (h) z Is or (h) z If .
The bus voltages affected by the harmonic branch currents can be obtained by substituting Equation ( 9) into Equation ( 5) and be represented as: Equation ( 10) can be further expressed as: ( ) The voltage of bus i in Equation ( 11) can therefore be rewritten as: ( ) where [BCBV i ] is the row vector of [BCBV] relating to bus i.
The voltages of the shunt equipments such as harmonic filters can also be calculated by multiplying their harmonic currents absorbed and harmonic impedance.For example, the voltage of shunt equipment at bus y can be expressed as: where ( )   h y V and ( ) h y Zf are the bus voltage and shunt impedance of the h th -order harmonic at bus y, respectively.
The bus voltages of shunt equipments can then be expressed in vector form as Equation ( 14).Equation ( 14) can be rearranged and written as Equation ( 15) and then Equation ( 16) can be derived.Therefore, the harmonic currents absorbed by the shunt equipments can be calculated by Equation (17).After the harmonic currents absorbed by shunt equipments are calculated, the branch currents and bus voltages affected by the h th -order harmonic currents can be calculated by Equations ( 9) and (11), respectively.The main equipment that needs to be considered in the distribution harmonic analysis includes transformers, capacitors, inductors and line sections: From Equation (17), it can be illustrated that for a distribution system with Nf shunt equipment, [ZHA (h) ] is a square matrix with dimension Nf.Instead of inverting the full admittance matrix of a distribution system, only the matrix with dimension Nf needs to be inverted; therefore, the proposed method has better computational performance.Besides, the proposed method can provide novel viewpoints in calculating the branch currents and bus voltages under harmonic pollution.
In general, a swing bus should be assigned for load flow calculation.For distribution systems, the swing bus is the substation bus and its voltage will maintain constant under all conditions.In harmonic analysis, the harmonic propagation for the substation bus can also be calculated if necessary.Namely, the harmonic currents can propagate through the distribution substation and beyond.In this situation, the short-circuit capacity in a substation bus can be used to integrate the substation bus into harmonic analysis.The strength of the substation bus is directly proportional to its short-circuit capacity.The Thevenin equivalent impedance looking back into the transmission system from the substation bus can be calculated by its short-circuit capacity [20].The Thevenin equivalent impedance in p.u. can be expressed as: Figure 2 illustrates the concepts of including the substation bus into harmonic analysis, and it can be seen that the substation bus can be treated as the Thevenin equivalent impedance of its short-circuit capacity.The transmission system can then be considered as an infinite bus with infinite short-circuit capacity and zero equivalent impedance and will maintain its voltage under all conditions.Therefore, the substation bus can be integrated into the proposed method without difficulty.The frequency-based component models as presented in [2][3][4]6] are used for the proposed harmonic analysis method.In a distribution system, the main devices considered in the harmonic analysis include distribution cable, transformer, capacitor, inductor, etc.Using the frequency-based component models, the line impedances required for the BCBV matrix under different harmonic orders are scalable.The BIBC matrix is a constant matrix and will not need to be modified in the proposed harmonic analysis method.A large-scale induction motor can be treated as a harmonic source with impedance.Therefore, the frequency-based component models can also be utilized.Other models can also be integrated into the proposed method if the relationship between the component impedance and harmonic frequency can be obtained.

Test Results and Discussions
Many distribution systems have been tested; however, only several cases are shown here to explain the solution procedures and demonstrate the validity of the proposed method.The fundamental frequency and base voltage for the following test cases are 60 Hz and 11.4 kV, respectively.

Simple Seven-Bus System
A simple seven-bus system as shown in Figure 3 is used as an example to build the matrices and vectors in Equation ( 16).The simple seven-bus system (equivalent to an 18-node system) consists of the three-phase, double-phase and single-phase line sections and buses.The simple seven-bus system and the nonlinear loads are not based on a real-world situation and are only used to explain the solution procedures of the proposed harmonic analysis method.Besides, for explanation purposes, many harmonic filters are added to the simple seven-bus system.Figure 3 illustrates that nonlinear loads and harmonic filters are interconnected to Buses 3 and 6, respectively.The matrices and vectors used for the proposed method are shown in the Appendix.To simplify the case explained in the Appendix, only the situation of nonlinear loads in Bus 3 as illustrated in Figure 3 is tested.Tables 1 and 2 show the line data and bus data for the simple seven-bus system of Figure 3, respectively.
0.3477 1.0141 0.1565 0.4777 ( / ) 0.1565 0.4777 0.3375 1.0478 bc j j mi j j Table 3 lists the bus voltages solved by the unbalanced distribution load flow method proposed in [17], and it can be seen that the feeder is operated under unbalanced conditions.The harmonic current magnitudes of nonlinear loads as illustrated in Figure 3 are listed in Table 4.Note that nonlinear loads are interconnected in Buses 3 and 5 for the following tests.The parameters of single tune harmonic filters including resonance frequency installed in Bus 6 are as listed in Table 5. Substituting the data from Tables 4 and 5 and Equation (19) into Equations (A1) and (A2) and Equations ( 16) and ( 17), the harmonic currents absorbed for each harmonic order can be calculated.For example, the 5 th -order harmonic currents absorbed by the harmonic filters installed in Bus 6 are 5.49 A, 2.91 A and 4.28 A for phases a, b and c, respectively.The branch currents and bus voltages for the 5 th -order harmonic order as illustrated in Figures 4 and 5 can be calculated by Equation ( 9) and (11), respectively.From Figures 4 and 5, the effects of harmonic current injections on branch currents and bus voltages can be clearly observed.After the bus voltages for each harmonic order are calculated, and the voltage THD (%), as listed in Table 6, can be obtained.The proposed harmonic analysis method obtains the same solution as those obtained by [6,12]; therefore, the validity of the proposed can be demonstrated.Since the solutions of the proposed method and the conventional methods are same, the solutions aren't shown in the paper.From the simple example, it can be seen that the proposed method can be used to calculate the bus voltages for each harmonic order and the voltage THD (%) effectively.If the reactance and capacitance parameters of a harmonic filter installed in bus 6a with a resonance frequency 282 Hz are changed to 36.8331Ω and 866.4 Ω (resonance frequency 291 Hz), respectively, the 5 th -order harmonic currents absorbed by the harmonic filters installed in Bus 6 will be changed to 8.103 A, 2.268 A and 3.641 A for phases a, b and c, respectively.The branch currents after harmonic filter parameter adjustment are as shown in Figure 6.The effects of harmonic filter parameter adjustment on branch currents and bus voltages can therefore be effectively determined.The proposed method can provide novel viewpoints in calculating the branch currents and bus voltages under harmonic pollution.These novel viewpoints have great potential to be used for optimal harmonic filter design and harmonic resonance mitigation.Due to limited space, the applications will be investigated in future work.
If the substation bus with short-circuit capacity 250 MVA is included in the harmonic analysis, the bus voltages for the 5 th -order harmonic order can then be calculated and illustrated in Figure 7.The voltage THD (%) of each bus is listed in Table 7.Note that the Thevenin equivalent impedance is 0.0004 p.u. with respect to the short-circuit capacity 250MVA and the base MVA 0.1 MVA.Comparisons of Figures 5 and 7 and Tables 6 and 7, the effects with and without substation bus for harmonic analysis can be clearly observed.From the simple test cases, it can be seen that the substation bus can be integrated into the proposed method without difficulty.

Feasibility and Performance Tests
An actual three-phase distribution system, including 6 laterals, 105 segments and 123 load points in the suburban area [21], is also used to test the proposed method.Different test scenarios have been used to demonstrate the feasibility of the proposed method; however, only one test scenario is shown here.Figure 8 illustrates that ten three-phase nonlinear loads and five three-phase harmonic filters are interconnected to the actual distribution system.The harmonic currents for the 5 th -order harmonic order are also illustrated in Figure 8. Obviously, the harmonic currents absorbed by the harmonic filters can be effectively calculated.The voltage THD (%) for each harmonic load and harmonic filter is illustrated in Figure 9.Note that the number of harmonic filters is 15, because five three-phase harmonic filters are installed in the actual distribution system.Only the [ZHA (h) ] matrix with a dimension of 15 needs to be inverted; therefore, the proposed method has better computational performance and has the feasibility to be used in actual distribution systems.
In order to demonstrate the computational performance of the proposed method, two other methods are used for comparisons.The three methods include: Method 1: the proposed method; Method 2: the commonly-used admittance-matrix-based method [10] and; Method 3: the forward/backward sweep based method [16].In order to demonstrate the computational performance of the proposed method, four three-phase systems, 15-, 32-, 63-, and 88-bus systems represented in [16], are used for our tests.Table 8 shows the normalized execution time (NET) of these three methods.From Table 8, it can be observed that the NETs of 88-bus system for Methods 1 and 2 are 17.53 and 3020.91 (about 172 times faster), respectively.Since the full admittance matrix inverse for each harmonic order can be omitted, Method 1 outperforms Method 2 for large-scale systems.The computational performances for Methods 1 and 3 are close.The NETs as shown in Table 8 include the computational time for load flow solution.The computational time of the forward/backward sweep based load flow method is less; therefore, the performance of Method 3 is better than Method 1.The NETs for Methods 1 and 3 excluding the computational time for a load flow solution are shown in Table 9.Table 9 shows that, excluding the computational time of load flow, Method 1 is slightly faster than Method 3. Although, the performance differences are not very obvious, the proposed method provides novel viewpoints in observing the branch currents and bus voltages under harmonic pollution in the solution procedures.Besides, Reference [12] needs the specific data format such as the coefficient vector defined to record whether the system harmonic currents flow through the branch to be designed and built first.The specific data format is then integrated into the commonly used forward/backward sweep techniques to calculate harmonic propagation.Therefore, the method proposed in [12] cannot be straightforwardly used to solve other harmonic problems, i.e., optimal harmonic filter design, etc.No specific data format needs to be designed and integrated in the solution procedure of Method 1.Moreover, the novel viewpoints in observing the branch currents and bus voltages under harmonic pollution which are vital for harmonic filter design can also be obtained by in the solution procedures of the proposed method.Therefore, the effectiveness and efficiency of the proposed method can be demonstrated.The phase information of the harmonic current sources can enhance the accuracy of harmonic analysis.The phase information for the different harmonic current sources in a distribution system needs a reference angle and the conventional power quality analyzer cannot determine the reference angle.Phasor-measurement unit can be used to acquire the phase information [22]; however, it is too expensive to be used in distribution systems.If the phase information of the harmonic current sources is acquired, the harmonic currents with phase information can be integrated into the proposed harmonic analysis method without modifying the proposed method.Without the phase information of the harmonic current sources in the above tests, the results can be considered as worst cases for harmonic analysis.

Conclusions
Two relationship matrices derived from the topological characteristics of the distribution systems were used to analyze the variations in bus voltages, bus current injections and branch currents affected by harmonic current injections.Accordingly, some useful formulas were derived to analyze the three-phase harmonic propagation problem.The proposed method can provide novel viewpoints in calculating the branch currents and bus voltages under harmonic pollution.The harmonic currents aborted by the shunt equipments can also be easily quantified in the proposed solution procedure.This data is vital for harmonic filter design.Several tests were conducted to demonstrate the performance and feasibility of the proposed method.Moreover, the proposed method can obtain the same solutions as those obtained by the commonly-used methods with less computational time; therefore, the proposed method has great potential to be used in on-line operation.The applications of the proposed method for optimal harmonic filter design and harmonic resonance mitigation will be discussed in future works. 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 1 0 0

Figure 1 .
Figure 1.Concepts of the proposed harmonic propagation.
equivalent impedance of the substation bus with respect to the short-circuit capacity (MVA SC ) and base MVA (MVA B ).

Figure 2 .
Figure 2. Including the substation bus in harmonic analysis.

Figure 3 .
Figure 3.A simple seven-bus test system.

Figure 4 .
Figure 4. Branch currents of the 5 th -order harmonic order.

Figure 5 .
Figure 5. Bus voltage of the 5 th -order harmonic order.

Figure 6 .
Figure 6.Branch currents of the 5 th -order harmonic order after harmonic filter parameter adjustment.

Figure 7 .
Figure 7. Bus voltage of the 5 th -order harmonic order with substation bus.

Figure 8 .
Figure 8. Harmonic currents of the 5 th -order harmonic order.

Table 1 .
Line data of the simple seven-bus test system.

Table 2 .
Load data of a simple seven-bus test system.

Table 3 .
Load flow solution of the simple seven-bus test system.

Table 4 .
Harmonic currents of nonlinear loads.

Table 5 .
Parameters of installed harmonic filters.

Table 6 .
Voltage total harmonic distortion (THD) (%) of the simple seven-bus test system.

Table 7 .
Voltage THD (%) of a simple seven-bus test system with substation bus.

Table 8 .
Performance comparisons.NET 1.0 is set for the 15-bus system in Method 1; + Ns and Nf are the numbers of nonlinear loads and shunt equipments, respectively. *:

Table 9 .
Performance comparisons for Methods 1 and 3 excluding the computational time of load flow.NET 1.0 is set for the 15-bus system in Method 1. *: