Harmonic Resonance Identiﬁcation and Mitigation in Power System Using Modal Analysis

: Due to a rising share of power electronic devices in power networks and the consequent rise in harmonic distortion, impedance resonances are an important issue. Nowadays, the frequency scan method is used for resonance phenomena identiﬁcation and analysis. The main disadvantage of the method is its inability to decouple different resonance phenomena. This means that the method is also unable to provide sufﬁcient information about the effects that the parameters of network elements have on different resonance phenomena. Furthermore, it was also noted that despite the fact that the harmonic resonance mode analysis is well described in the literature, there is a lack of systematic approach to the analysis procedure. Thus the main objective of this paper is to address this disadvantage and to propose a systematic approach to harmonic resonance analysis and mitigation, utilizing modal analysis. In the ﬁrst part of the paper, dominant network nodes in terms of resonance ampliﬁcation of harmonics are determined. This is done by analysis of the eigenvalues of the network admittance matrix. Using the eigenvalue analysis results, key parameters of network elements involved in a speciﬁc resonance are determined next. This is performed by calculating the critical mode (i.e., the mode that experiences resonance) sensitivity coefﬁcients with respect to network parameters. In the second part of the paper, the procedure for modal resonance frequency shift is presented. The shift is performed by changing the value of a selected parameter so that the modal resonance frequency matches the desired resonance frequency value. The parameter value is calculated with the Newton–Rhapson method. Presented analysis considers both parallel and series resonances. The effectiveness of the proposed method is demonstrated on an actual industrial-network model.


Introduction
Harmonic distortion of voltages and currents in power systems is an inevitable phenomenon, due to the nonlinear nature of loads and elements constantly present across different voltage levels. In particular, this applies to industrial networks that have witnessed a large increase in the number of connected power converters, supplying various loads [1]. Excessive levels of harmonic distortion may be detrimental to the network and connected devices, as it can disrupt their operation and, in extreme cases, even cause their failure. In order to limit these negative effects, recommended guidelines determine the levels of permitted harmonic distortion in the supply voltage [2].
Under normal operating conditions, harmonic distortion does not generally exceed limits set by the standards [3]. However, with harmonic resonance conditions present in the network, voltages and/or currents can get amplified significantly. Resonance has been reported as one of the main harmonic-related causes for problems with electrical equipment [4]. To prevent equipment failure, appropriate methods for harmonic resonance identification and mitigation are essential. Having said that, these methods need to be used In the last part of the paper, the proposed approach to harmonic resonance analysis is tested by means of simulations on an actual industrial network model.

Harmonic Resonance Identification
The presented resonance mode analysis method consists of several steps summarized in Figure 1. The most important are network admittance matrix calculation, admittance matrix eigen-decomposition, modal sensitivity indices calculation and modal resonance frequency shift. In order to perform the series resonance analysis, the admittance matrix needs to be modified as it is presented in Figure 1. The modification is comprised of identification of the network bus, for which the series resonance analysis is performed, and elimination of the row and the column related to the analyzed network bus. The steps presented in Figure 1 are explained in the next subsections as follows. The basic theory of modal analysis and the meaning of transformation of the network is explained in Section 2.1. Next, the modal and the resonance frequency indices are presented in Section 2.2. Together with the sensitivity indices, the algorithm for calculation of the accurate value of the resonance frequency is presented in the subsection. Finally, modifications needed for series resonance analysis are presented in Section 2.3. The resonance frequency shift is addressed later, in Section 3.
Energies 2021, 14, x FOR PEER REVIEW 3 of 21 selected and, using the Newton method, the values of selected parameters at which resonance frequencies reach the desired values were determined. In the last part of the paper, the proposed approach to harmonic resonance analysis is tested by means of simulations on an actual industrial network model.

Harmonic Resonance Identification
The presented resonance mode analysis method consists of several steps summarized in Figure 1. The most important are network admittance matrix calculation, admittance matrix eigen-decomposition, modal sensitivity indices calculation and modal resonance frequency shift. In order to perform the series resonance analysis, the admittance matrix needs to be modified as it is presented in Figure 1. The modification is comprised of identification of the network bus, for which the series resonance analysis is performed, and elimination of the row and the column related to the analyzed network bus. The steps presented in Figure 1 are explained in the next subsections as follows. The basic theory of modal analysis and the meaning of transformation of the network is explained in SubSection 2.1. Next, the modal and the resonance frequency indices are presented in SubSection 2.2. Together with the sensitivity indices, the algorithm for calculation of the accurate value of the resonance frequency is presented in the subsection. Finally, modifications needed for series resonance analysis are presented in SubSection 2.3. The resonance frequency shift is addressed later, in Section 3.  Figure 1. Flowchart of the presented analysis method.

Resonance Mode Analysis
The power network can be represented using the admittance matrix Y. In this case, the nodal voltage vector U is equal to

Resonance Mode Analysis
The power network can be represented using the admittance matrix Y. In this case, the nodal voltage vector U is equal to where I denotes the nodal current injections vector. The admittance matrix representation of the network is chosen because of the simplicity of its construction. Moreover, the parallel resonance's relation with admittance matrix Y being close to singularity has been discov-  [7]. It has also been noted by the authors that one of the eigenvalues of the singular matrix is equal to zero. That is why the admittance matrix eigen-decomposition needs to be performed in order to do a more detailed study of the resonances.
Admittance Matrix Eigen-Decomposition The harmonic resonance mode analysis is based on presentation of electrical network with the admittance matrix Y, which can be decomposed into [7]: where L and T are the left and the right eigenvectors matrix, respectively, and Λ is the eigenvalues matrix, which is diagonal. The elements of the eigenvalue matrix are eigenvalues λ i , also named modes, where subscript i denotes mode number, and its dimension is equal to the dimension of the admittance matrix. Calculation of the inverse value of the eigenvalue matrix yields the modal impedances matrix Λ −1 . Since the eigenvalue matrix is diagonal, the elements of the Λ −1 matrix are equal to λ −1 i . They are named modal impedances. Each modal impedance resonance determines one nodal impedance parallel frequency. The voltage equation of the network presented in the modal domain can be written as where V and J are modal voltages vector and modal current injections vector, respectively. The corresponding eigenvectors matrices determine the parallel resonance effect on the network nodes. The T matrix can be used to determine modal current injection vector J as where I denotes the nodal current injections vector. The influence of the right eigenvector matrix on the degree of effect of a current injection into a certain network bus to the particular modal current can be seen from Equation (4). The right eigenvectors are the rows of the right eigenvector matrix. That is why the right eigenvectors give information on modal excitability [7]. The L matrix determines the modal voltage effect on the nodal voltages U as U = LV.
Equation (5) describes the influence of the left eigenvector matrix to the degree of effect of a certain modal voltage expressed in the bus voltages. Namely, it gives information on modal observability [7]. The left eigenvectors are the columns of the left eigenvector matrix.
As the admittance matrix Y is symmetrical, the relationship between the right and left eigenvalue matrix is T = L −1 and L −1 = L T , where the operator T indicates matrix transposition. This means that the bus with the highest observability level also has the highest excitability level. Thus only one of the two eigenvectors gives enough information for the harmonic resonance analysis [7].
Since the right and left eigenvectors essentially give the same information, they can be combined into the participation factors [7]. They are calculated using the equation where the subscripts i and j denote mode number and bus number, respectively.

Modal Sensitivity Analysis
To shift the modal resonance frequency as effectively as possible, the effect of each network parameter on the modal impedance needs to be evaluated. The general parameter is marked with the letter α and represents either the inductance, capacitance or resistance of Energies 2021, 14, 4017 5 of 20 a network element. Figure 2 presents the effect of changing the parameter α on the modal impedance characteristic. The full red line represents a modal impedance characteristic before the parameter change and the dotted red line represents a modal impedance characteristic after the change. It can be observed that the modal impedance magnitude and the modal resonance frequency may change.

Modal Sensitivity Analysis
To shift the modal resonance frequency as effectively as possible, the effect of each network parameter on the modal impedance needs to be evaluated. The general parameter is marked with the letter and represents either the inductance, capacitance or resistance of a network element. Figure 2 presents the effect of changing the parameter on the modal impedance characteristic. The full red line represents a modal impedance characteristic before the parameter change and the dotted red line represents a modal impedance characteristic after the change. It can be observed that the modal impedance magnitude and the modal resonance frequency may change. To fully understand the particular parameter effect, both, the modal impedance magnitude and the resonance frequency sensitivity need to be calculated.

Modal Impedance Magnitude Sensitivity
The modal impedance magnitude together with the participation factor determines the effect of a nodal harmonic current injection on the network nodal voltages magnitude. That means, in order to mitigate a harmonic voltage magnitude, the modal impedance magnitude needs to be minimized. The modal impedance magnitude sensitivity is calculated through the eigenvalue magnitude sensitivity using the process thoroughly described in [8].
In summary, the eigenvalue sensitivity is calculated using the equation where and represent the real and imaginary parts of the network branch admittance, respectively. In Equation (7)  To fully understand the particular parameter effect, both, the modal impedance magnitude and the resonance frequency sensitivity need to be calculated.

Modal Impedance Magnitude Sensitivity
The modal impedance magnitude together with the participation factor determines the effect of a nodal harmonic current injection on the network nodal voltages magnitude. That means, in order to mitigate a harmonic voltage magnitude, the modal impedance magnitude needs to be minimized. The modal impedance magnitude sensitivity is calculated through the eigenvalue magnitude sensitivity using the process thoroughly described in [8].
In summary, the eigenvalue sensitivity is calculated using the equation where G and B represent the real and imaginary parts of the network branch admittance, respectively. In Equation (7), ∂G ∂α and ∂B ∂α are derivatives of the series RLC branch admittance with respect to one of the parameters α (R, X L or B C ). The derivatives And where S i denotes sensitivity coefficient for analyzed element. The sensitivity coefficients for the i-th mode are calculated from entries of the sensitivity matrix, where the operator L :,i and T i,: are left and right eigenvector of i-th mode, respectively. The calculation of S i for the element to which the analyzed parameter relates is performed depending on how the element is connected. In case it is connected as a shunt element, S i is equal to S i = S i (j, j), where j is the node that the element is connected to. If the element is connected in series between nodes j and k, then the sensitivity coefficient is equal to The modal impedance magnitude sensitivity is equal to In order to achieve the comparability between parameter changes done with different kinds of elements, the result from Equation (11) is normalized using equation During the calculation process, it was observed that the value of the modal impedance magnitude sensitivity depends on the calculation step ∆ f . The smaller the calculation step, the better the sensitivity coefficient approximation is. This is because the density of the calculation points on the observed frequency interval increases, and the identified resonance frequency converges to the real resonance frequency. To avoid the dependence on the calculation step size, an accurate resonance frequency calculation procedure is proposed. The calculation is performed using the modified idea of bisection, which is an iterative numerical method for calculation of zeros of a function. The algorithm is based on the assumption that the modal impedance function is continuous and concave in the surroundings of the resonance frequency. The difference between the bisection and the proposed method is the latter being used to determine the frequency at which the maximum of the modal impedance characteristic is located.
The calculation can be summarized in the following steps.
1. Import i-th mode data, identify the resonance frequency f r,i , and set the counter to t = 1.

2.
Set , and modal impedance amplitudes z m,i d (t) and z m,i e (t) .

3.
Find at which frequency the modal impedance magnitude is higher (point d (t) or point Set the new resonance frequency approximation for t-th iteration f (t) where ε stop is the stopping criterion, continue to step 2, else construct an array with elements where δ f is the diminished calculation interval, and finish the calculation.

Resonance Frequency Sensitivity
To obtain full information on the system parameter effect on modal impedance, the modal resonance frequency sensitivity also needs to be calculated. Namely, the example with a parallel RLC circuit has shown, that the inductor and the capacitor do not have any effect on the modal impedance magnitude. However, they have an effect on its resonance frequency [9]. The resonance frequency sensitivity has been calculated using the method described in [9,10]. It is calculated using the equation The first part of the Equation (13), where u i is evaluated in Equation (8) and v i in Equation (9). The derivatives of u i and v i are and The expressions ∂G ∂α , ∂B ∂α , ∂ 2 G ∂α∂ f and ∂ 2 B ∂α∂ f from the Equation (14) are derivatives of the series RLC branch impedance with respect to parameter α and to frequency f . The second part of the Equation (13), The resonance frequency sensitivity, evaluated with the Equation (13) is further normalized [9,10] using the equation

Series Resonance Modal Analysis
It has been shown that applying the resonance mode analysis method to the inverse nodal admittance matrix Y −1 is not appropriate for the series resonance analysis [7,11]. Since series resonance relates to high loop current values, a loop impedance matrix Z loop should be constructed and analyzed. As harmonic voltage is applied between a network node and the reference node to excite a harmonic current, the observed node should be connected to the reference node through the dummy branch (a short circuit connection) [11]. Such connection is necessary because the impedance of an ideal voltage source is equal to zero.
The main drawback of the loop impedance matrix is the complexity of its construction, which is more complex in comparison with the construction of the nodal admittance matrix. Because of that, Neufeld et al. [12] proposed the representation of the network, whose topology is changed with the dummy branch method, using the nodal admittance matrix.
The implemented network nodal admittance matrix modification for series resonance modal analysis was performed by deleting the row and the column referring to the analyzed node. Because of that, the original node numbers were saved for the later analysis of the results. The network modification algorithm is summarized in Figure 3. the analyzed node. Because of that, the original node numbers were saved for the later analysis of the results. The network modification algorithm is summarized in Figure 3.

Harmonic Resonance Mitigation
Based on the results of resonance mode and modal sensitivity analysis from the previous section, the resonance frequency shift method is introduced below to shift critical modes to safe frequencies.

Resonance Frequency Shift
Since the high values of parallel resonance modal impedances at their resonance frequencies coincide with high values of certain bus impedances, it is necessary to shift these resonance frequencies as far as possible from the frequencies at which harmonic current injections are expected. Alternatively, the bus impedance magnitude may be lowered by the shift of the resonance frequencies of the series resonance modal impedances. In this case, the resonance frequencies of the modal impedances coincide with low values of certain bus impedance which means that the resonance frequencies of the series resonances modal impedance should be shifted as close as possible to the frequencies of the expected harmonic current injections.
The function m, ( , ) describing the dependence between the modal impedance magnitude, frequency and parameter size, is calculated numerically. Consequentially, the chosen parameter value at which the resonance frequency is equal to the desired value is calculated using the Newton-Raphson method [10]. The zero of function h(a) is calculated using the equation

Harmonic Resonance Mitigation
Based on the results of resonance mode and modal sensitivity analysis from the previous section, the resonance frequency shift method is introduced below to shift critical modes to safe frequencies.

Resonance Frequency Shift
Since the high values of parallel resonance modal impedances at their resonance frequencies coincide with high values of certain bus impedances, it is necessary to shift these resonance frequencies as far as possible from the frequencies at which harmonic current injections are expected. Alternatively, the bus impedance magnitude may be lowered by the shift of the resonance frequencies of the series resonance modal impedances. In this case, the resonance frequencies of the modal impedances coincide with low values of certain bus impedance which means that the resonance frequencies of the series resonances modal impedance should be shifted as close as possible to the frequencies of the expected harmonic current injections.
The function |z m,i ( f , α)| describing the dependence between the modal impedance magnitude, frequency and parameter size, is calculated numerically. Consequentially, the chosen parameter value at which the resonance frequency is equal to the desired value is calculated using the Newton-Raphson method [10]. The zero of function h(a) is calculated using the equation where a is the vector of changed parameters α, f r,current (a) is the current resonance frequencies vector at parameter vector a and f r,desired is the vector of desired resonance frequencies.
The calculation is performed using the equation until the stopping criterion a (t+1) − a (t) < ε stop is not satisfied. In Equation (20), J represents the Jacobi matrix, which consists of frequency sensitivity coefficients. During the resonance frequency shift calculations, it was observed that resonance mode numbering may change. Because of this, a procedure of the expected modal impedance resonance frequency and magnitude assessment has been implemented. The modal resonance frequency estimationf (t+1) r,i after the parameter changes in t-th iteration is calculated with equationf while the modal impedance magnitude is calculated using the equation where ∆ zm is modal impedance sensitivity coefficients matrix for which ∂α p , where p denotes a parameter number. After the evaluation of the Equations (21) and (22) the modal impedance closest to the estimated values is chosen for resonance frequency shift in the next iteration.

Application Example
To illustrate the practical implementation of the proposed method and to evaluate its efficiency, the analysis is demonstrated on a real industrial network model. A simplified scheme of the analyzed system is shown in Figure 4. where a is the vector of changed parameters , f r,current (a) is the current resonance frequencies vector at parameter vector a and f r,desired is the vector of desired resonance frequencies. The calculation is performed using the equation until the stopping criterion a ( ) − a ( ) < stop is not satisfied. In Equation (20), J represents the Jacobi matrix, which consists of frequency sensitivity coefficients.
During the resonance frequency shift calculations, it was observed that resonance mode numbering may change. Because of this, a procedure of the expected modal impedance resonance frequency and magnitude assessment has been implemented. The modal resonance frequency estimation r, where zm is modal impedance sensitivity coefficients matrix for which Δ zm ( , ) = m,

( )
, where denotes a parameter number. After the evaluation of the Equations (21) and (22) the modal impedance closest to the estimated values is chosen for resonance frequency shift in the next iteration.

Application Example
To illustrate the practical implementation of the proposed method and to evaluate its efficiency, the analysis is demonstrated on a real industrial network model. A simplified scheme of the analyzed system is shown in Figure 4.

System Data
As can be observed from Figure 4, there are three 110/20 kV transformers (TR) in the network: TR 1 and TR 2 supplying system II and TR 3 supplying system I. Through tertiary winding of transformer TR 3, a 5 kV system is supplied. Altogether, there are more

System Data
As can be observed from Figure 4, there are three 110/20 kV transformers (TR) in the network: TR 1 and TR 2 supplying system II and TR 3 supplying system I. Through tertiary winding of transformer TR 3, a 5 kV system is supplied. Altogether, there are more than 20 passive reactive power compensators connected to the low voltage level 0.4 kV of the system, which are, due to the simplicity, excluded from the analysis. Three central compensators K 1, K 2 and K 4, connected to the medium-voltage level (5 kV and 20 kV) and one (K 3), connected to the low voltage have been considered. Due to the poorly designed compensation scheme, which caused rather complex circumstances in terms of resonant amplification of harmonics present in the network, this network represents an interesting example for application of the proposed resonance analysis method from Figure 1. The network parameters are given in Table 1.

Simulation Results
The analysis of the system resonances was performed according to the algorithm, presented in Figure 1. The parallel and series resonance frequencies are identified from the system bus frequency scan results. Further analyses are performed for both parallel and series resonances. They consist of modal eigen-decomposition, modal resonance frequency and modal impedance magnitude sensitivity calculation and finally the bus impedance mitigation with modal resonance frequency shift. Figure 5 presents the frequency scan results of the industrial network. There are two different frequency scan representations, one with absolute per unit impedance values and one with values normalized to the impedance at the nominal network frequency. Two bus impedance frequency scan representations are given because, firstly, modal decomposition is performed on the network admittance matrix, which yields the absolute impedance values when bus current injection is performed, and consequentially modal impedances, eigenvalues and eigenvectors relate to them. On the other hand, it is important to know the relative impedance magnitude increase in comparison to the value at the nominal frequency and, as can be seen from Figure 5, the resulting information about harmonic resonance severity at different busses may change dramatically.

Simulation Results
The analysis of the system resonances was performed according to the algorithm, presented in Figure 1. The parallel and series resonance frequencies are identified from the system bus frequency scan results. Further analyses are performed for both parallel and series resonances. They consist of modal eigen-decomposition, modal resonance frequency and modal impedance magnitude sensitivity calculation and finally the bus impedance mitigation with modal resonance frequency shift. Figure 5 presents the frequency scan results of the industrial network. There are two different frequency scan representations, one with absolute per unit impedance values and one with values normalized to the impedance at the nominal network frequency. Two bus impedance frequency scan representations are given because, firstly, modal decomposition is performed on the network admittance matrix, which yields the absolute impedance values when bus current injection is performed, and consequentially modal impedances, eigenvalues and eigenvectors relate to them. On the other hand, it is important to know the relative impedance magnitude increase in comparison to the value at the nominal frequency and, as can be seen from Figure 5, the resulting information about harmonic resonance severity at different busses may change dramatically.   From Figure 5, several frequencies can be seen at which the potential for parallel resonances exists. The most critical are resonances at 5th and 7th harmonic components. Due to the limited space, the results are shown for the resonance at the 7th harmonic frequency only.
Modal impedance characteristics of the industrial system are shown in Figure 6. Modes with the highest value and, consequentially, the biggest influence on the bus parallel resonances are mode 8 (resonance frequency 220 Hz), mode 7 (resonance frequency 291 Hz), mode 5 (resonance frequency 340 Hz) and mode 6 (resonance frequency 435 Hz). frequency only.
Modal impedance characteristics of the industrial system are shown in Figure 6. Modes with the highest value and, consequentially, the biggest influence on the bus parallel resonances are mode 8 (resonance frequency 220 Hz), mode 7 (resonance frequency 291 Hz), mode 5 (resonance frequency 340 Hz) and mode 6 (resonance frequency 435 Hz). As shown in Figure 5, there is a significant possibility for parallel resonance at the 7th harmonic frequency at busses 1, 5, 6, 7, 8 and 9. As stated before, the modal impedances (shown in Figure 6) are related to the absolute per unit bus impedance values, which is also confirmed by the eigenvectors and the participation factors shown in Table  2. Among the mentioned 6 busses, bus 5 was not analyzed as it is a mathematical accessory and thus does not physically exist, and bus 1 whose impedance magnitude increase at the 7th harmonic frequency is much lower than for the other busses in question. The modal impedance of mode 5 has the resonance frequency closest to the frequency of the 7th harmonic (350 Hz) as well as the highest value at that frequency. This means As shown in Figure 5, there is a significant possibility for parallel resonance at the 7th harmonic frequency at busses 1, 5, 6, 7, 8 and 9. As stated before, the modal impedances (shown in Figure 6) are related to the absolute per unit bus impedance values, which is also confirmed by the eigenvectors and the participation factors shown in Table 2. Among the mentioned 6 busses, bus 5 was not analyzed as it is a mathematical accessory and thus does not physically exist, and bus 1 whose impedance magnitude increase at the 7th harmonic frequency is much lower than for the other busses in question.

Parallel Resonance Analysis
The modal impedance of mode 5 has the resonance frequency closest to the frequency of the 7th harmonic (350 Hz) as well as the highest value at that frequency. This means that mode 5 needs to be shifted in order to mitigate the potential parallel resonance at the 7th harmonic component.
From the eigenvectors and participation factors presented in Table 2, it can be seen that the highest impedance magnitude is in bus 9, followed by the impedance magnitude of bus 6, which is consistent with per unit absolute impedance values from Figure 5. The normalized nodal impedance characteristics, however, show that the highest relative impedance magnitude increase is in bus 6, followed by impedances of the buses 9, 7, 5 and 8. Since there is the highest impedance magnitude increase, it can be concluded that the most severe 7th harmonic parallel resonance consequences can be expected in bus 6. The modal impedance magnitude and resonance frequency sensitivity coefficients of the parameters with the biggest influence on resonance frequency are presented in Figure 7 and in Table 3.
that the highest impedance magnitude is in bus 9, followed by the impedance magnitude of bus 6, which is consistent with per unit absolute impedance values from Figure 5. The normalized nodal impedance characteristics, however, show that the highest relative impedance magnitude increase is in bus 6, followed by impedances of the buses 9, 7, 5 and 8. Since there is the highest impedance magnitude increase, it can be concluded that the most severe 7th harmonic parallel resonance consequences can be expected in bus 6.
The modal impedance magnitude and resonance frequency sensitivity coefficients of the parameters with the biggest influence on resonance frequency are presented in Figure  7 and in Table 3.  From Figure 8 and Table 3 can be observed that the elements with the highest frequency sensitivity coefficients are located in the nodes from 6 to 9, where the reactive power compensators are located. The susceptance of compensator K 1 has the highest resonance frequency sensitivity. That is why it was chosen for performing the 7th harmonic parallel resonance mitigation.

Parallel Resonance Frequency Shift Mitigation
In Figures 8 and 9, impedances seen from busses 6 and 7 before shift, after shift and without compensator in operation are shown, respectively. As it can be seen, the parallel  From Figure 8 and Table 3 can be observed that the elements with the highest frequency sensitivity coefficients are located in the nodes from 6 to 9, where the reactive power compensators are located. The susceptance of compensator K 1 has the highest resonance frequency sensitivity. That is why it was chosen for performing the 7th harmonic parallel resonance mitigation.

Parallel Resonance Frequency Shift Mitigation
In Figures 8 and 9, impedances seen from busses 6 and 7 before shift, after shift and without compensator in operation are shown, respectively. As it can be seen, the parallel resonance frequency in the vicinity of the 7th harmonic has been shifted to a lower frequency. The frequency of the mode 5 resonance has been shifted from 340 Hz to 320 Hz. The shift was achieved by changing the compensator K 1 capacitive susceptance from 0.0415 pu to 0.0671 pu. The new parameter value was calculated in 5 iterations. Due to the limited space, only results for busses 6 and 7 are shown, however, similar results have been achieved also at other busses of interest (busses 8 and 9).  The relations between bus 6 and 7 impedances shown in Figures 8 and 9 are presented in Figure 10 at the 7th harmonic frequency (350 Hz). The values presented are normalized to the impedance magnitudes at the nominal network frequency. The impedance after the parallel resonance frequency shift with parameter change is lower, which means that the parallel resonance mitigation at the 7th harmonic frequency was successful. The relations between bus 6 and 7 impedances shown in Figures 8 and 9 are presented in Figure 10 at the 7th harmonic frequency (350 Hz). The values presented are normalized to the impedance magnitudes at the nominal network frequency. The impedance after the parallel resonance frequency shift with parameter change is lower, which means that the parallel resonance mitigation at the 7th harmonic frequency was successful.

Series Resonance Analysis
The series resonance analysis studies for the 7th harmonic were performed for each bus individually because of the network topology changes, needed to be done for every bus. The results of the series resonance analyses (for all buses of interest) are shown in Table 4.

Series Resonance Analysis
The series resonance analysis studies for the 7th harmonic were performed for each bus individually because of the network topology changes, needed to be done for every  Table 4.   Figure 11 shows the modal impedances frequency characteristics for the analyzed network when bus 6 is short-circuited. The number of modes is decreased by one, comparing to the unmodified system (see Figure 6). That is a consequence of a decreased number of rows and columns of the modified system admittance matrix in accordance to the original system. The excluded row and column belong to the short-circuited bus. Figure 11 shows that mode 5 has the resonance frequency closest to the series resonance of bus 6 impedance at 374 Hz, which was decided to be shifted. The series resonance (sensitivity) analysis results for bus 6 are presented in Figure 12. The results show that the series resonance frequency of bus 6 impedance at 325 Hz is determined by mode 5. The compensator K 1 has the biggest effect on its resonance frequency, 373 Hz, namely, its capacitance and inductance. These elements would thus be the most appropriate choice for parallel resonance mitigation.
of the parameters of compensator K 1.
The series resonance of node 9 at frequency 325 Hz is determined by resonance mode 5. The resonance frequency of mode 5 is 328 Hz, and, as the results show, the capacitance of compensator K 1 and inductance K 2 have the biggest effect on the series resonance. Consequentially, one of the two elements is the most appropriate for the 7th harmonic parallel resonance mitigation in bus 9.   The series resonance of node 9 at frequency 325 Hz is determined by resonance mode 5. The resonance frequency of mode 5 is 328 Hz, and, as the results show, the capacitance of compensator K 1 and inductance K 2 have the biggest effect on the series resonance. Consequentially, one of the two elements is the most appropriate for the 7th harmonic parallel resonance mitigation in bus 9.   Series resonance analysis results for bus 7 show that the closest series resonance frequency to the 7th harmonic is at 300 Hz, and that the capacitance of compensator K 3 and inductance of transformer TR 6 have the biggest influence on it. Choosing between these two elements, it would be more appropriate to select the capacitor K 3 capacitance. Alternatively, series filters K 1 and K 2 are connected to bus 7, so they could be tuned to frequency enabling as efficient parallel resonance mitigation as possible.
The series resonance of bus 8 that is closest to 7th harmonic frequency is at 302 Hz, which is determined by mode 7 with resonance frequency at 301 Hz. The results show the capacitance of compensator K 3 and the inductance of transformer TR 6 having the biggest influence.
There is another series resonance in bus 8 impedance that could be useful for 7th harmonic mitigation. It is located at a frequency of 461 Hz and is determined by mode 1. According to the presented results it can be said that it is determined by compensator K 1, connected directly to bus 7. Since both series resonances seem to be connected, it appears that the bus 8 parallel resonance at 7th harmonic should be mitigated using one of the parameters of compensator K 1.
The series resonance of node 9 at frequency 325 Hz is determined by resonance mode 5. The resonance frequency of mode 5 is 328 Hz, and, as the results show, the capacitance of compensator K 1 and inductance K 2 have the biggest effect on the series resonance. Consequentially, one of the two elements is the most appropriate for the 7th harmonic parallel resonance mitigation in bus 9.

Series Resonance Frequency Shift Mitigation
The 7th harmonic parallel resonance mitigation by shifting series resonances has been performed for each bus individually.
The resonance frequency of mode 5 has been shifted from 373 Hz to 345 Hz using the capacitive susceptance of compensator K 1. The parameter has changed from 0.0415 pu to 0.0556 pu and series resonance frequency of bus 6 impedance characteristic has changed from 374 Hz to 336 Hz. The calculation was completed in 4 iterations. Figure 13 shows the comparison between the bus 6 impedance characteristics before and after the change and the characteristic without any compensation.  For bus 8, the 7th parallel harmonic resonance mitigation has been achieved by mode 1 resonance frequency shift from 482 Hz to 330 Hz. The series resonance observed at bus 8 was shifted from 461 Hz to 345 Hz. The shift was performed by increasing the capacitive susceptance of the compensator K 1 from 0.00415 pu to 0.00721 pu. The calculation was For bus 7 impedance it has been decided that the most appropriate way to mitigate the parallel resonance at 7th harmonic frequency is to tune either compensator K 1 or K 2 to frequency close to the 7th harmonic. Between the two compensators, K 1 was chosen. The shift presented in Figure 14 has been achieved by changing the tuning frequency of the compensator K 1 from 483 Hz to 345 Hz.
For bus 8, the 7th parallel harmonic resonance mitigation has been achieved by mode 1 resonance frequency shift from 482 Hz to 330 Hz. The series resonance observed at bus 8 was shifted from 461 Hz to 345 Hz. The shift was performed by increasing the capacitive susceptance of the compensator K 1 from 0.00415 pu to 0.00721 pu. The calculation was finished in 3 iterations. The comparison of the bus 8 frequency impedance characteristic is shown in Figure 15.
In Figure 16, the results of the 7th harmonic parallel resonance mitigation of bus 9 with series resonance frequency shift are presented. This was achieved by shifting the resonance frequency of mode 5 from 328 Hz to 340 Hz with change of the capacitor K 2 inductance. The bus 9 impedance series resonance has changed from 325 Hz to 357 Hz and the inductive reactance has decreased from 1.772 pu to 1.423 pu. The calculation was finished in 5 iterations.  For bus 8, the 7th parallel harmonic resonance mitigation has been achieved by mode 1 resonance frequency shift from 482 Hz to 330 Hz. The series resonance observed at bus 8 was shifted from 461 Hz to 345 Hz. The shift was performed by increasing the capacitive susceptance of the compensator K 1 from 0.00415 pu to 0.00721 pu. The calculation was finished in 3 iterations. The comparison of the bus 8 frequency impedance characteristic is shown in Figure 15.
In Figure 16, the results of the 7th harmonic parallel resonance mitigation of bus 9 with series resonance frequency shift are presented. This was achieved by shifting the resonance frequency of mode 5 from 328 Hz to 340 Hz with change of the capacitor K 2 inductance. The bus 9 impedance series resonance has changed from 325 Hz to 357 Hz and the inductive reactance has decreased from 1.772 pu to 1.423 pu. The calculation was finished in 5 iterations.     The relations between the impedances of busses 6, 7, 8 and 9, shown in Figures 13-16, are presented in Figure 17 at the 7th harmonic frequency. The values presented are normalized to the impedance magnitudes at the nominal network frequency. The parallel resonances of the busses in question have been successfully mitigated with the series resonance frequency shift too. The biggest mitigation effect was achieved for busses 7 and 8, while for the other two busses, the impedance is lowered but is still larger comparing with the case when the compensation in the network is absent.

Conclusions
The main objective of this paper is to propose a systematic approach to harmonic resonance identification and mitigation in power system using modal analysis. The main steps in the proposed calculation procedure are as follows. First, the network is transformed into the modal domain. In this way, individual resonance phenomena are separated and represented with modal impedances, enabling their individual analysis. The influence of elements and their parameters at the nominal frequency was included in the analysis. This was achieved by calculating the sensitivity of the amplitude of modal impedances and resonant frequencies of modal impedances to changes in parameters. The analysis of series resonances was included in the analysis with appropriate adjustments of the admittance matrix, thus enabling the use of the same computational procedures as for parallel resonances. Finally, with the help of sensitivity coefficients, the parameter values of network elements at which the adjusted resonance frequency of the modal impedance reaches the desired value were determined using the Newton method. The effectiveness of the proposed method was evaluated on a real industrial network model.
As highlighted in the paper, the following solutions are presented.
• When calculating the sensitivity coefficients, it was observed that sensitivity coefficients of modal impedances are dependent on the size of a calculation step Δf. The problem was solved by proposing an additional calculation of a more accurate value of the resonance frequency. • Another problem, that was encountered, was the question of how to follow the shifted mode when changing the resonance frequency. Namely, it was observed that the frequency characteristic of one modal impedance can have several resonant frequencies, and the fact that the labels of some resonant modes changed in a way that was not exactly consistent when the parameters were changed made it even more difficult. This obstacle was solved by estimation of expected resonance frequency and modal impedance magnitude in every iteration of the resonance frequency shift calculation, so that after changing the parameters in the new iteration of the calculation, we chose the resonance mode that was closest to the estimated value.

Conclusions
The main objective of this paper is to propose a systematic approach to harmonic resonance identification and mitigation in power system using modal analysis. The main steps in the proposed calculation procedure are as follows. First, the network is transformed into the modal domain. In this way, individual resonance phenomena are separated and represented with modal impedances, enabling their individual analysis. The influence of elements and their parameters at the nominal frequency was included in the analysis. This was achieved by calculating the sensitivity of the amplitude of modal impedances and resonant frequencies of modal impedances to changes in parameters. The analysis of series resonances was included in the analysis with appropriate adjustments of the admittance matrix, thus enabling the use of the same computational procedures as for parallel resonances. Finally, with the help of sensitivity coefficients, the parameter values of network elements at which the adjusted resonance frequency of the modal impedance reaches the desired value were determined using the Newton method. The effectiveness of the proposed method was evaluated on a real industrial network model.
As highlighted in the paper, the following solutions are presented.
• When calculating the sensitivity coefficients, it was observed that sensitivity coefficients of modal impedances are dependent on the size of a calculation step ∆f. The problem was solved by proposing an additional calculation of a more accurate value of the resonance frequency. • Another problem, that was encountered, was the question of how to follow the shifted mode when changing the resonance frequency. Namely, it was observed that the frequency characteristic of one modal impedance can have several resonant frequencies, and the fact that the labels of some resonant modes changed in a way that was not exactly consistent when the parameters were changed made it even more difficult. This obstacle was solved by estimation of expected resonance frequency and modal impedance magnitude in every iteration of the resonance frequency shift calculation, so that after changing the parameters in the new iteration of the calculation, we chose the resonance mode that was closest to the estimated value.
Overall, the results show that the transformation of a network, presented with the admittance matrix, into modal domain is a useful tool for individual analysis of each individual resonance phenomenon. In relation to the modal sensitivity indices, it has been noted that together with the enhanced resonance frequency calculation, the proposed approach reliably determines the degree of effect that every individual parameter has on modal impedance characteristics. Finally, it can be seen from the results, that using the proposed resonance frequency shift method, the bus impedances were successfully lowered at the frequency in question (the 7th harmonic). Consequently, the degree of bus voltage harmonic distortion, caused by nodal harmonic current injections, is successfully mitigated.