Harmonic Resonance Identification and Mitigation in Power System Using Modal Analysis

Abstract

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 identification 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 sufficient 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 first part of the paper, dominant network nodes in terms of resonance amplification 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 specific resonance are determined next. This is performed by calculating the critical mode (i.e., the mode that experiences resonance) sensitivity coefficients 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.

1. 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 both existing networks when topology is changing or new equipment is being connected to the system, as well as in the network planning procedure.

The most frequently used harmonic resonance analysis method is an impedance frequency scan [5,6]. It is a relatively straightforward analysis technique, where the impedance seen from a certain bus is calculated as a function of frequency. This technique shows if parallel and/or series resonance exist in a system. However, it does not give any answers as to which network components are involved in the resonance, which network nodes are most affected, where is the best location, and what is the best approach to mitigate the problem [7].

These questions have been addressed partially in the literature by several authors. The idea of the harmonic resonance mode analysis (HRMA) has been introduced by Xu et al. [7] as a tool for parallel harmonic resonance analysis. Compared to the conventional impedance frequency scan analysis method, the harmonic resonance mode analysis enables individual analysis of resonance phenomena. The method is based on admittance matrix eigen-decomposition, which yields eigenvalues and eigenvectors used to determine bus participation in a particular resonance phenomenon. Obtained results can be further used to determine the modal sensitivity. The concept of the modal impedance magnitude sensitivity was presented by Huang et al. [8]. Furthermore, modal impedance resonance frequency can also be calculated using the method presented by Cui and Wang [9] and altered for use with the complex admittance matrix by Hu et al. [10]. All of the aforementioned calculations can be performed for series harmonic resonance analysis with modifications presented by Zhou et al. [11]. With all mentioned calculations, parameters of the specific network elements with the highest influence on the resonance frequencies can be selected. Using the Newton–Rhapson method presented by Hu et al. for use with the HRMA method in [10], the values of parameters can be altered in such a way that resonance frequencies are shifted to the values, at which the harmonic voltage amplification is minimized.

While the principles of the harmonic resonance mode analysis were thoroughly described in the aforementioned literature, it was noted that there is a lack of systematic approach to performing the analysis comprehensively. Additionally, two more issues, that were not addressed in the literature, were encountered during the implementation of the method. The first issue was the dependence of the modal sensitivity on the selected frequency calculation step. Specifically, by decreasing the calculation step the sensitivity indices converge to the actual values. The second issue was noticed during the modal resonance frequency shift. It was observed, that the numbering of the modal impedances may change after the new value of the adjusted parameter is calculated during a Newton–Rhapson method iteration. This makes it difficult to keep track of the adjusted modal impedance resonance frequency properly without any further measures.

The main objective of this paper is therefore to propose a systematic approach to harmonic resonance (parallel and series) analysis and mitigation utilizing the HRMA method. Additionally, the algorithm for calculation of more accurate resonance frequency values of the resonance modes is presented to address the issues with calculation step dependence of modal sensitivity indices. The final proposed improvement is the algorithm for tracking the adjusted modal impedance resonance frequency. In described approach, the harmonic resonance mode analysis through admittance matrix eigen-decomposition is first presented. Modal sensitivity analyses of impedance magnitude and resonance frequency are performed next to determine the influence of network parameters on both quantities. To improve the accuracy of the modal resonance frequencies calculation, an accurate resonance frequency calculation procedure is proposed. Based on the results of all the previous calculations, the harmonic resonance mitigation method is presented next, to adjust the resonance frequencies of the modal impedances. Parameters with the highest influence on resonance frequencies were 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.

2. 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.

2.1. Resonance Mode Analysis

The power network can be represented using the admittance matrix $\mathbf{Y}$. In this case, the nodal voltage vector $\mathbf{U}$ is equal to

$$\mathbf{U}={\mathbf{Y}}^{-1}\mathbf{I},$$

(1)

where $\mathbf{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 $\mathbf{Y}$ being close to singularity has been discovered by Xu et al. [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 $\mathbf{Y}$, which can be decomposed into [7]:

$$\mathbf{Y}=\mathbf{L}\Lambda \mathbf{T},$$

(2)

where $\mathbf{L}$ and $\mathbf{T}$ are the left and the right eigenvectors matrix, respectively, and $\Lambda$ is the eigenvalues matrix, which is diagonal. The elements of the eigenvalue matrix are eigenvalues ${\lambda }_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 ${\Lambda }^{-1}$. Since the eigenvalue matrix is diagonal, the elements of the ${\Lambda }^{-1}$ matrix are equal to ${\lambda }_i^{-1}$. 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

$$\mathbf{V}={\Lambda }^{-1}\mathbf{J},$$

(3)

where $\mathbf{V}$ and $\mathbf{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 $\mathbf{T}$ matrix can be used to determine modal current injection vector $\mathbf{J}$ as

$$\mathbf{J}=\mathbf{TI},$$

(4)

where $\mathbf{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 $\mathbf{L}$ matrix determines the modal voltage effect on the nodal voltages $\mathbf{U}$ as

$$\mathbf{U}=\mathbf{LV}.$$

(5)

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 $\mathbf{Y}$ is symmetrical, the relationship between the right and left eigenvalue matrix is $\mathbf{T}={\mathbf{L}}^{-1}$ and ${\mathbf{L}}^{-1}={\mathbf{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

$$P{F}_{i,j}=L_{j,i}\cdot T_{i,j},$$

(6)

where the subscripts $i$ and $j$ denote mode number and bus number, respectively.

2.2. 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 $\alpha$ and represents either the inductance, capacitance or resistance of a network element. Figure 2 presents the effect of changing the parameter $\alpha$ 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.

2.2.1. 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

$$\frac{\partial \left|{\lambda }_i\right|}{\partial \alpha }=\frac{\partial \left|{\lambda }_i\right|}{\partial G}\cdot \frac{\partial G}{\partial \alpha }+\frac{\partial \left|{\lambda }_i\right|}{\partial B}\cdot \frac{\partial B}{\partial \alpha }$$

(7)

where $G$ and $B$ represent the real and imaginary parts of the network branch admittance, respectively. In Equation (7), $\frac{\partial G}{\partial \alpha }$ and $\frac{\partial B}{\partial \alpha }$ are derivatives of the series RLC branch admittance with respect to one of the parameters $\alpha$ ($R$, $X_L$ or $B_C$). The derivatives $\frac{\partial \left|{\lambda }_i\right|}{\partial G}$ and $\frac{\partial \left|{\lambda }_i\right|}{\partial B}$ are equal to

$$\frac{\partial \left|{\lambda }_i\right|}{\partial G}=u_i=\frac{S_{i,\mathrm{Re}}\cdot {\lambda }_{i,\mathrm{Re}}+S_{i,\mathrm{Im}}\cdot {\lambda }_{i,\mathrm{Im}}}{\sqrt{{\lambda }_{i,\mathrm{Re}}^2+{\lambda }_{i,\mathrm{Im}}^2}}$$

(8)

And

$$\frac{\partial \left|{\lambda }_i\right|}{\partial B}=v_i=\frac{S_{i,\mathrm{Re}}\cdot {\lambda }_{i,\mathrm{Im}}-S_{i,\mathrm{Im}}\cdot {\lambda }_{i,\mathrm{Re}}}{\sqrt{{\lambda }_{i,\mathrm{Re}}^2+{\lambda }_{i,\mathrm{Im}}^2}},$$

(9)

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,

$${\mathbf{S}}_i={\mathbf{L}}_{:,i}\cdot {\mathbf{T}}_{i,:},$$

(10)

where the operator ${\mathbf{L}}_{:,i}$ and ${\mathbf{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={\mathbf{S}}_i\left(j,j\right)$, 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 $S_i={\mathbf{S}}_i\left(j,j\right)+{\mathbf{S}}_i\left(k,k\right)-{\mathbf{S}}_i\left(j,k\right)-{\mathbf{S}}_i\left(k,j\right)$.

The modal impedance magnitude sensitivity is equal to

$$\frac{\partial \left|z_{\mathrm{m},i}\right|}{\partial \alpha }=-{\left|{\lambda }_i\right|}^{-2}\cdot \frac{\partial \left|{\lambda }_i\right|}{\partial \alpha }.$$

(11)

In order to achieve the comparability between parameter changes done with different kinds of elements, the result from Equation (11) is normalized using equation

$${\frac{\partial \left|z_{\mathrm{m},i}\right|}{\partial \alpha }|}_{\mathrm{norm}}=\frac{\partial \left|z_{\mathrm{m},i}\right|}{\partial \alpha }\cdot \frac{\alpha }{\left|z_{\mathrm{m},i}\right|}.$$

(12)

During the calculation process, it was observed that the value of the modal impedance magnitude sensitivity depends on the calculation step $\mathsf{\Delta }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.

• Import $i$-th mode data, identify the resonance frequency $f_{\mathrm{r},i}$, and set the counter to $t=1$.
• Set $a^{\left(t-1\right)}=f_{\mathrm{r},i}-\mathsf{\Delta }f$, $b^{\left(t-1\right)}=f_{\mathrm{r},i}$ and $c^{\left(t-1\right)}=f_{\mathrm{r},i}+\mathsf{\Delta }f$, then calculate $d^{\left(t\right)}=\frac{a^{\left(t-1\right)}+b^{\left(t-1\right)}}{2}$, $e^{\left(t\right)}=\frac{b^{\left(t-1\right)}+c^{\left(t-1\right)}}{2}$, and modal impedance amplitudes $\left|z_{\mathrm{m},i}\left(d^{\left(t\right)}\right)\right|$ and $\left|z_{\mathrm{m},i}\left(e^{\left(t\right)}\right)\right|$.
• Find at which frequency the modal impedance magnitude is higher (point $d^{\left(t\right)}$ or point $e^{\left(t\right)}$) and set it as $b^{\left(t\right)}$. If $\left|z_{\mathrm{m},i}\left(d^{\left(t\right)}\right)\right|>\left|z_{\mathrm{m},i}\left(e^{\left(t\right)}\right)\right|$, set $a^{\left(t\right)}=a^{\left(t-1\right)}$ and $c^{\left(t\right)}=b^{\left(t-1\right)}$, else $a^{\left(t\right)}=b^{\left(t-1\right)}$ and $c^{\left(t\right)}=c^{\left(t-1\right)}$.
• Set the new resonance frequency approximation for $t$-th iteration $f_{\mathrm{r},i}^{\left(t\right)}=b^{\left(t\right)}$.
• If $\left|\left|z_{\mathrm{m},i}\left(b^{\left(t\right)}\right)\right|-\left|z_{\mathrm{m},i}\left(b^{\left(t-1\right)}\right)\right|\right|>{\epsilon }_{\mathrm{stop}}$, where ${\epsilon }_{\mathrm{stop}}$ is the stopping criterion, continue to step 2, else construct an array with elements $\left[f_{\mathrm{r},i}-2\delta f,f_{\mathrm{r},i}-\delta f,f_{\mathrm{r},i},f_{\mathrm{r},i}+\delta f,f_{\mathrm{r},i}+2\delta f\right]$, where $\delta f$ is the diminished calculation interval, and finish the calculation.

2.2.2. 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

$$\frac{\partial f_{\mathrm{r},i}}{\partial \alpha }=\frac{\frac{{\partial }^2\left|{\lambda }_i\right|}{\partial f\partial \alpha }}{\frac{{\partial }^2\left|{\lambda }_i\right|}{\partial f^2}}.$$

(13)

The first part of the Equation (13), $\frac{{\partial }^2\left|{\lambda }_i\right|}{\partial f\partial \alpha }$, is equal to

$$\frac{{\partial }^2\left|{\lambda }_i\right|}{\partial f\partial \alpha }=\frac{\partial u_i}{\partial f}\frac{\partial G}{\partial \alpha }+u_i\frac{{\partial }^{2}G}{\partial \alpha \partial f}+\frac{\partial v_i}{\partial f}\frac{\partial B}{\partial \alpha }+v_i\frac{{\partial }^{2}B}{\partial \alpha \partial f},$$

(14)

where $u_i$ is evaluated in Equation (8) and $v_i$ in Equation (9). The derivatives of $u_i$ and $v_i$ are

$$\frac{\partial u_i}{\partial f}=\frac{u_i\left(f_{\mathrm{r},i}+\mathsf{\Delta }f\right)-u_i\left(f_{\mathrm{r},i}-\mathsf{\Delta }f\right)}{2\mathsf{\Delta }f}$$

(15)

and

$$\frac{\partial v}{\partial f}=\frac{v_i\left(f_{\mathrm{r},i}+\mathsf{\Delta }f\right)-v_i\left(f_{\mathrm{r},i}-\mathsf{\Delta }f\right)}{2\mathsf{\Delta }f}.$$

(16)

The expressions $\frac{\partial G}{\partial \alpha }$, $\frac{\partial B}{\partial \alpha }$, $\frac{{\partial }^{2}G}{\partial \alpha \partial f}$ and $\frac{{\partial }^{2}B}{\partial \alpha \partial f}$ from the Equation (14) are derivatives of the series RLC branch impedance with respect to parameter $\alpha$ and to frequency $f$.

The second part of the Equation (13), $\frac{{\partial }^2\left|{\lambda }_i\right|}{\partial f^2}$, is equal to

$$\frac{{\partial }^2\left|{\lambda }_i\right|}{\partial f^2}=\frac{\left|{\lambda }_i\left(f_{\mathrm{r},i}+2\mathsf{\Delta }f\right)\right|+\left|{\lambda }_i\left(f_{\mathrm{r},i}-2\mathsf{\Delta }f\right)\right|-2\left|{\lambda }_i\left(f_{\mathrm{r},i}\right)\right|}{4{\left(\mathsf{\Delta }f\right)}^2}.$$

(17)

The resonance frequency sensitivity, evaluated with the Equation (13) is further normalized [9,10] using the equation

$${\frac{\partial f_{\mathrm{r},i}}{\partial \alpha }|}_{\mathrm{norm}}=\frac{\partial f_{\mathrm{r},i}}{\partial \alpha }\cdot \alpha .$$

(18)

2.3. Series Resonance Modal Analysis

It has been shown that applying the resonance mode analysis method to the inverse nodal admittance matrix ${\mathbf{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 ${\mathbf{Z}}_{\mathrm{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.

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 $\left|z_{\mathrm{m},i}\left(f,\alpha \right)\right|$ 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 $\mathbf{h}\mathbf{\left(}\mathbf{a}\mathbf{\right)}$ is calculated using the equation

$$\mathbf{h}\mathbf{\left(}\mathbf{a}\mathbf{\right)}={\mathbf{f}}_{\mathrm{r},\mathrm{current}}\mathbf{\left(}\mathbf{a}\mathbf{\right)}-{\mathbf{f}}_{\mathrm{r},\mathrm{desired}},$$

(19)

where $\mathbf{a}$ is the vector of changed parameters $\alpha$, ${\mathbf{f}}_{\mathrm{r},\mathrm{current}}\left(\mathbf{a}\right)$ is the current resonance frequencies vector at parameter vector $\mathbf{a}$ and ${\mathbf{f}}_{\mathrm{r},\mathrm{desired}}$ is the vector of desired resonance frequencies. The calculation is performed using the equation

$${\mathbf{a}}^{\left(t+1\right)}={\mathbf{a}}^{\left(t\right)}-{\mathbf{J}}^{-1}\left({\mathbf{a}}^{\left(t\right)}\right)\cdot \mathbf{h}\left({\mathbf{a}}^{\left(t\right)}\right)$$

(20)

until the stopping criterion $\left|{\mathbf{a}}^{\left(t+1\right)}-{\mathbf{a}}^{\left(t\right)}\right|<{\epsilon }_{\mathrm{stop}}$ is not satisfied. In Equation (20), $\mathbf{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 ${\widehat{f}}_{\mathrm{r},i}^{\left(t+1\right)}$ after the parameter changes in $t$-th iteration is calculated with equation

$${\widehat{f}}_{\mathrm{r},i}^{\left(t+1\right)}=f_{\mathrm{r},i}^{\left(t\right)}+\mathbf{J}\left(a^{\left(t\right)}\right)\cdot \Delta \mathbf{a}$$

(21)

while the modal impedance magnitude is calculated using the equation

$$\left|z_{\mathrm{m},i}^{\left(t+1\right)}\right|=\left|z_{\mathrm{m},i}^{\left(t\right)}\right|+{\Delta }_{\mathrm{zm}}\left({\mathbf{a}}^{\left(t\right)}\right)\cdot \Delta \mathbf{a},$$

(22)

where ${\Delta }_{\mathrm{zm}}$ is modal impedance sensitivity coefficients matrix for which ${\mathsf{\Delta }}_{\mathrm{zm}}\left(i,p\right)=\frac{\partial \left|z_{\mathrm{m},i}^{\left(t+1\right)}\right|}{\partial {\alpha }_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.

4. 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.

4.1. 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. Industrial network parameters.

Component	Parameter	Value
TR 1	S	20 MVA
	uk	10.36%
	R/X	1/10
TR 2	S	20 MVA
	uk	10.36%
	R/X	1/10
TR 3	S	40/30/10 MVA
	uk	10.05/11.89/7.84%
	R/X	1/10
TR 4	S	36 MVA
	uk	4.11%
	R/X	1/10
TR 5	S	8 MVA
	uk	8.07%
	R/X	1/10
TR 6	S	1 MVA
	uk	5%
	R/X	1/5
K 1	S	6.3 MVAr
	ft	395 Hz
	R/XL	1/10
	R/XC	1/200
K 2	S	2.1 MVAr
	ft	264 Hz
	R/XL	1/10
	R/XC	1/200
K 3	S	0.6 MVAr
	ft	2542 Hz
	R/XL	1/10
	R/XC	1/200
K 4	S	2 MVAr
	ft	980 Hz
	R/XL	1/10
	R/XC	1/200
Load @ 3	P	54 MW
	Q	46 MVAr
Load @ 4	P	6 MW
	Q	3.5 MVAr
Load @ 7	P	10 MW
	Q	10 MVAr
Load @ 9	P	0.8 MW
	Q	0.5 MVAr
Load @ 6	P	3 MW
	Q	2 MVAr
Network	Sk”	1100 MVA
	R/X	1/10
Line 7–8	R	0.162 Ω
	X	0.185 Ω

.

4.2. 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).

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. Eigenvectors and participation factors of the industrial system.

Quantity	Node (j)	${\mathit{f}}_{\mathbf{r}\mathbf{,}\mathbf{5}}$ 340 Hz	Node (j)	${\mathit{f}}_{\mathbf{r}\mathbf{,}\mathbf{7}}$ 290 Hz
${\lambda }_i$/pu	-	0.00747	-	0.00343
${\mathit{PF}}_{i,j}$	9	0.8350	9	0.9825
	6	0.1616	6	0.0074
	5	0.0433	8	0.0046
	7	0.0401	5	0.0041
	8	0.0370	7	0.0036
	1	0.0015	1	0.0001
	2	0.0009	2	<0.0001
	3	0.0007	3	<0.0001
	4	0.0007	4	<0.0001

. 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.

4.2.1. 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. Resonance sensitivity coefficients of mode 5.

Parameter	$\frac{\partial {\mathit{f}}_{\mathbf{r},5}}{\partial \mathit{\alpha }}/\mathbf{pu}/100\%$	$\frac{\partial \left|{\mathit{z}}_{\mathbf{m},5}\right|}{\partial \mathit{\alpha }}/\%/100\%$
K 1–BC	−1.3231	1.0199
K 2–XL	−0.9089	1.3848
K 4–BC	−0.9050	0.6603
K 1–XL	−0.6476	0.5047
TR 6–XL	−0.6031	−2.3231
K 2–BC	−0.5606	0.8360
K 3–BC	−0.5420	−3.4423

.

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.

4.2.2. Parallel Resonance Frequency Shift Mitigation

In Figure 8 and Figure 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 Figure 8 and Figure 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.

4.2.3. 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 sensitivity coefficients for relevant modes and busses.

Mode 5 (bus 6)
Parameter	$\frac{\partial f_{\mathrm{r},5}}{\partial \alpha }$/pu/100%	$\frac{\partial \left|z_{\mathrm{m},5}\right|}{\partial \alpha }$/%/100%
K 1–BC	−3.5135	−1.1661
K 1–XL	−2.0975	−0.6944
TR 6–XL	−0.4423	−0.9174
TR 3 (primary winding)–XL	−0.4190	0.7630
K 2–XL	−0.3549	1.2650
Mode 7 (bus 7)
Parameter	$\frac{\partial f_{\mathrm{r},5}}{\partial \alpha }$/pu/100%	$\frac{\partial \left|z_{\mathrm{m},5}\right|}{\partial \alpha }$/%/100%
K 3–BC	−3.0065	−0.9836
TR 6–XL	−2.7491	0.5279
Load, bus 9–XL	−0.2114	0.4512
Mode 7 (bus 8)
Parameter	$\frac{\partial f_{\mathrm{r},5}}{\partial \alpha }$/pu/100%	$\frac{\partial \left|z_{\mathrm{m},5}\right|}{\partial \alpha }$/%/100%
K 3–BC	−3.0180	−0.9837
TR 6–XL	−2.7845	0.5303
Load, bus 9–XL	−0.2106	0.4534
Mode 1 (bus 8)
Parameter	$\frac{\partial f_{\mathrm{r},5}}{\partial \alpha }$/pu/100%	$\frac{\partial \left|z_{\mathrm{m},5}\right|}{\partial \alpha }$/%/100%
K 1–BC	−4.5736	−0.9496
K 1–XL	−3.9764	−0.8222
Line 7−8–XL	−0.3360	1.5542
K 1–R	−0.1426	−0.7339
Mode 5 (bus 9)
Parameter	$\frac{\partial f_{\mathrm{r},5}}{\partial \alpha }$/pu/100%	$\frac{\partial \left|z_{\mathrm{m},5}\right|}{\partial \alpha }$/%/100%
K 1–BC	−1.5664	−1.6133
K 2–XL	−1.1790	1.1064
K 4–BC	−0.9425	−0.0948
K 2–BC	−0.7719	−0.7339
K 1–XL	−0.7282	0.7179
TR 3 (primary winding)–XL	−0.4908	−0.0825

.

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.

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.

4.2.4. 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 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.

The relations between the impedances of busses 6, 7, 8 and 9, shown in Figure 13, Figure 14, Figure 15 and Figure 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.

5. 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.