Frequency-dependent Parameter Identification for Improved Dynamic State Estimation Based Protection Based on Characteristic Signal Injection of Half-bridge MMC in Flexible DC Grids

Ming Nie; Jinghan He; Meng Li; Huiyuan Zhang; Keao Chen

网刊加载中。。。

使用Chrome浏览器效果最佳，继续浏览，你可能不会看到最佳的展示效果，

确定继续浏览么?

复制成功，请在其他浏览器进行阅读

Frequency-dependent Parameter Identification for Improved Dynamic State Estimation Based Protection Based on Characteristic Signal Injection of Half-bridge MMC in Flexible DC Grids PDF

- ORCID：
Ming Nie ¹
✉
- ORCID：
Jinghan He ¹ (Fellow, IEEE)
✉
- ORCID：
Meng Li ¹
✉
- ORCID：
Huiyuan Zhang ²
✉
- ORCID：
Keao Chen ¹
✉

1. School of Electrical Engineering, Beijing Jiaotong University, Beijing100044, China； 2. China Electric Power Research Institute, State Grid Corporation of China, Beijing100192, China

Updated：2024-12-18

DOI：10.35833/MPCE.2023.000800

OUTLINE

Abstract

In flexible DC grids, the rapid rise of fault current requires that the line protection must complete the fault identification within a few milliseconds. Dynamic state estimation based protection (DSEBP) provides a new idea for flexible DC line protection with good performance. However, the operating frequency in the DC grid is 0 Hz. When the DC grid is operating normally, it is difficult to identify the line parameters online to improve the performance of the protection method. This paper proposes a method to identify the frequency-dependent parameters of flexible DC grids based on the characteristic signal injection of half-bridge modular multilevel converter (HB-MMC). The characteristic signal is extracted by the Prony algorithm to calculate the line parameter under different frequencies. Afterwards, the number and position of residues and poles of frequency-dependent parameters are determined using the vector fitting method. Finally, an improved DSEBP is proposed. The simulation shows that the frequency-dependent parameters obtained by the proposed parameter identification method can be used in the improved DSEBP normally, and the identified parameters have better precision.

Keywords

Parameter identification; modular multilevel converter (MMC); state estimation; protection; flexible DC grid

I. Introduction

FLEXIBLE DC grid offers some advantages, including the absence of commutation failure, the ability to provide island power supply, and convenience for networking. These advantages have led to the widespread application of the flexible DC grid in power delivery and asynchronous network interconnection [

1]-[3]. However, the weakness of power electronics to withstand overcurrent is contradictory with the rapid rise of fault current. Improper handling of a single fault may even affect the normal operation of the whole DC grid. Hence, the effective implementation of transmission line protection in flexible DC grids plays a vital role in ensuring the secure operation of the power system.

In the current flexible DC grid projects, the main protection methods usually consist of the voltage derivative-based protection method and traveling wave protection method [

4]. The backup protection method usually involves the under-voltage protection method and current differential protection method [5]. Nevertheless, it should be noted that the efficacy of the main protection methods is impacted by fault resistance. The operation speed of the current differential protection is affected by the distribution capacitance, and the selectivity of the under-voltage protection needs to be improved. Therefore, it is a big challenge to coordinate operation speed and sensitivity of the protection [6], [7].

To improve the protection performance, the dynamic state estimation method is used for flexible DC protection [

8]. Dynamic state estimation based protection (DSEBP) considers the frequency-dependent parameter characteristics to weaken the effect of transient processes and improve the protection reliability. The DSEBP provides a new idea to solve the inherent contradiction in the coordination of protection methods. But the precise frequency-dependent parameters are difficult to obtain. Therefore, it is very important to obtain accurate parameters of flexible DC transmission lines to enhance the practicality of DSEBP.

Parameter identification methods can be divided into offline measurement methods and online measurement calculation methods [

9]. The principle of the offline measurement method changes the wiring of the line according to different identified parameters. Subsequently, the relevant identified parameters are calculated by the measured voltage, current, and other basic quantities [10]. The offline measurement method is practical and easy to operate, but it needs to build an independent measurement circuit or use an additional power supply at the test site, which leads to many operational steps, obvious errors, and other shortcomings. In [11], the open-circuit impedance and short-circuit impedance of the transmission line are measured offline and calculated by Taylor series expansion. Thereafter, the relationship between distributed conductance and distributed impedance is established and the relevant parameters are solved. The primary disadvantage of offline measurement is that the transmission line must be converted to the state of outage maintenance, which disrupts the normal operation of the system.

Currently, the majority of researchers are engaged in the study of online parameter identification methods. The classical theory method, based on Carson formula, is used to calculate the line parameters through geometrical and structural configurations of the towers, conductors, and ground wires [

12]. However, the accuracy of this method cannot meet the engineering requirements, as the calculation process depends on various simplifications. The ongoing developments in phasor measurement unit (PMU) technology and improvements in supervisory control and data acquisition (SCADA) systems facilitate further advances in line parameter identification [13]-[15]. In [16], a nonlinear network parameter identification method based on a hybrid measurement system combining wide-area measurement system (WAMS) and SCADA is proposed. Its principle is relatively simple, yet the line model employs an equivalent model, which results in a slight lack of accuracy. In [17], an online parameter estimation method based on measured PMU data is proposed for transmission lines. The accurate formula for calculating equivalent parameters of transmission lines is derived based on a uniform transmission line equation. In addition to PMU measurement, some researchers propose the use of fault records from the protection relays for parameter identification [18]. However, the current method of transmission line parameter identification focuses mainly on AC lines, whereas the research on DC line parameter identification is notably scarce. The main reason is that the operating frequency in the DC grid is 0 Hz, which makes it difficult to identify parameters in a steady state.

Fault identification based on the characteristic signal injection in DC grids is one of the research highlights [

19], [20]. It uses the controllable characteristics of modular multilevel converters (MMCs) to adjust the input quantity of submodules by changing the control strategy. Then, the generated characteristic signal of the specific frequency is injected into the transmission line. According to the principles of traveling wave reflection, the faults can be identified.

To solve the above issues, a frequency-dependent parameter identification method based on characteristic signal injection of half-bridge modular multilevel converter (HB-MMC) is proposed. First, by controlling the HB-MMC, the specific frequency signal is injected into the line, enabling the accurate calculation of the frequency-dependent parameter. Afterward, the residues and poles of the frequency-dependent parameter can be obtained by the vector fitting method. Finally, the performance of DSEBP is improved based on the available residues and poles. The key contributions are as follows.

1) The redundant submodules of HB-MMC are used for characteristic signal injection, which solves the difficulty in parameter identification caused by the operating frequency of 0 Hz in the DC grid.

2) The proposed parameter identification method offers a precise calculation of frequency-dependent parameters, thereby enhancing the accuracy of parameter identification.

3) The utilization of the available residues and poles associated with the frequency-dependent parameter contributes to the development of the improved DSEBP. This improved method is deemed to be more practical and reliable in its application.

II. Characteristic Signal Injection of HB-MMC

A. Topology and Control Principles of HB-MMC

The topology of HB-MMC is shown in Fig. 1. The configuration comprises six bridge arms, wherein the upper and lower bridge arms are subsumed into a phase unit. Each bridge arm consists of multiple half-bridge submodules (HBSMs), characterized by identical structural composition, connected in series with the bridge inductance denoted as L_larm. HBSM consists of two insulated gate bipolar transistors (IGBTs) T1 and T2, capacitor C₀, and two anti-shunt diodes D1 and D2.

Fig. 1 Topology of HB-MMC.

HBSM has three operating statuses and can output two voltage levels U_c and 0. In Fig. 1, i_dc is the DC current; U_dc is the DC voltage; $i_{p φ}$ and $i_{n φ}$ ( $φ = a, b, c$ ) are the three-phase currents in the upper and lower bridge arms, respectively; and $u_{g φ}$ and $i_{g φ}$ are the three-phase voltage and current on the AC side of the converter, respectively.

Figure 2 shows the traditional control strategy for the HB-MMC during normal system operation, where PI stands for proportional-integral; PLL stands for phase-locked loop; NLM stands for nearest level modulation; $θ$ is the phase angle of $u_{g φ}$ ; P is the active power; P_ref is the reference value of the active power; U_dcref is the reference value of the DC voltage; Q is the reactive power; Q_ref is the reference value of the reactive power; I_d is the d-axis current; I_q is the q-axis current; U_d is the d-axis voltage; U_q is the q-axis voltage; $u_{d i f f φ}$ is the differential mode voltage; u_dcn is the rated value of the DC voltage; $u_{r e f}$ is the voltage reference for the bridge arm; $u_{r e f p φ}$ is the voltage reference for the upper bridge arm; and $u_{r e f n φ}$ is the voltage reference for the lower bridge arm.

Fig. 2 Traditional control strategy for HB-MMC.

The HB-MMC employs a traditional double closed-loop vector control strategy, wherein half of the submodule capacitors are engaged to sustain the DC voltage. The outer-loop control strategy involves the regulation of active power and reactive power, while the inner-loop control strategy focuses on current control.

B. Principle of Characteristic Signal Generation

Without considering the redundant submodules of the converter station, the relationship among the series number of bridge arm submodules N, the DC voltage U_dc, and the rated capacitance voltage U_cN can be expressed as:

N = U_{d c} / U_{c N}

(1)

In the design of MMC, the operating condition with modulation factor $K_{m} = 1$ should be considered. At this time, the voltage variation range of the upper and lower bridge arms is 0-U_dc, and in consequence, the number of submodules inserted by each bridge arm is 0-N. In this case, the additional submodules in series outside the N submodules are the redundant submodules that are considered in the design, and the number is N_de. However, in the actual operation of MMC, the modulation factor K_m is less than 1, and the maximum number of submodules to be invested by each bridge arm N_op can be expressed as:

N_{o p} = \frac{1 + K_{m}}{2} N

(2)

According to the above analysis, when the modulation factor K_m<1, the bridge arm submodules will not be fully utilized. At any given time, at least $N - N_{o p}$ submodules are in the bypass state, which are called operational redundant submodules $Δ N_{o p}$ . Therefore, when MMC is in normal operation, the total number of redundant submodules N_re is expressed as:

N_{r e} = N_{d e} + Δ N_{o p}

(3)

In the field of engineering, the most commonly employed MMC modulation strategies are those based on NLM, which allows for the theoretical control of the voltage difference between the MMC output voltage and the modulation wave within ±U_dc/2. The reference value of DC voltage is defined as U_dcref, then the number of submodules to be incorporated into the upper and lower bridge arms at time t can be expressed as:

\{\begin{array}{l} N_{p φ} (t) = \frac{N_{o p}}{2} - r o u n d (\frac{u_{v φ}^{*}}{U_{c N}}) \\ N_{n φ} (t) = \frac{N_{o p}}{2} + r o u n d (\frac{u_{v φ}^{*}}{U_{c N}}) \end{array}

(4)

where round(x) is the function to take the integer closest to x; and $u_{v φ}^{*}$ is the reference sinusoidal modulated voltage.

In this paper, the injection of characteristic signals is achieved through the control of input and the removal of redundant submodules. In other words, N_op is employed to modify the number and state of equivalent HBSM inputs. For the same phase, if the input submodules are simultaneously increased or decreased for both the upper and lower bridge arms, the reference potential of the AC-side voltage will remain unchanged. Therefore, the AC side of the MMC will not be affected.

The additional control strategy of signal injection is shown in Fig. 3. The moment at which the characteristic signal is injected serves as the initial moment, and the number of additional submodules $N_{Δ φ} (t)$ that need to be invested by each bridge arm at moment t is calculated.

Fig. 3 Additional control strategy of signal injection.

In Fig. 3, f_a is the frequency of the characteristic signal; $1 / f_{a}$ is the time of one period of the characteristic signal; T_D is the initial moment when the characteristic signal is at a high level; T_S is the initial moment when the characteristic signal is at a low level; n is the intermediate variable which is related to the duration of characteristic signal; and $Δ N$ is the number of additional submodules that are invested by each bridge arm.

The square wave is injected as the characteristic signal. When the injected characteristic signal is at a high level, the upper and lower bridge arms of the MMC will input $Δ N$ additional submodules. When the injected characteristic signal is at a low level, the upper and lower bridge arms of the MMC will remove $Δ N$ additional submodules. Afterwards, the number of additional submodules $Δ N$ is added to the number of normal operating submodules to be put into operation during the process of generating injected signals.

C. Selection of Characteristic Signals

In actual projects, it is typical for the bridge arms of converters to contain redundant submodules. For instance, the Nan’ao three-terminal flexible DC project in China comprises $134 + 13$ (redundant) submodules in the bridge arm of the Sucheng converter station and $200 + 20$ (redundant) submodules in the bridge arm of the Jinniu and Qing’ao converter stations. The ±500kV Zhangbei flexible DC grid project in China comprises $264 + 20$ (redundant) submodules in the bridge arm. Therefore, when HB-MMC is in normal operation, the number of redundant submodules can be employed for characteristic signal injection.

However, the selection of the characteristic signal necessitates consideration of a variety of factors, including the submodule switching control period, the limitations of practical engineering, and the duration of characteristic signal injection.

The choice of the injected characteristic signal frequency should accurately reflect the alterations in both the line characteristic impedance and attenuation function. Additionally, it is crucial for the characteristic signal frequency to be less than the sampling frequency of the protection system, which is inherently constrained by the switching speed of the submodules. In practical applications of high-voltage direct current (HVDC) projects, the typical sampling frequency ranges from 10 to 50 kHz, with a control cycle of approximately 100 $μ s$ . For instance, considering a sampling frequency of 10 kHz, the protection can correctly distinguish frequencies up to 5000 Hz according to the Shannon sampling theorem.

The amplitude of the injected characteristic signal is constrained by the number of redundant submodules. Many factors should be considered, such as the performance of the voltage/current transformer, the potential impact on the power system, and the overcurrent tolerance capability of power electronics. Consequently, the amplitude of the injected characteristic signal is constrained by rigorous limitations. However, it is worth noting that a larger amplitude of the injected characteristic signal can enhance the effectiveness of signal extraction. For power grids operating at 35 kV and above, the permissible voltage deviation range is typically set at 10% of the rated voltage. Additionally, the redundant submodules in the bridge arm usually account for 10% to 20% of the total. In this study, the amplitude of the injected characteristic signal is chosen to be 5% of the rated DC voltage.

For the duration of characteristic signal injection, the accuracy of signal extraction should be considered. First, if the duration of the characteristic signal injection is less than one period, the accuracy of the signal extraction is reduced. Second, during the first few periods of characteristic signal injection, some small disturbances may occur due to the switching of the control strategy. It is possible to await the stabilization of the injected characteristic signal before the signal extraction. Therefore, the duration is set to be ten periods of the characteristic signal.

D. Extraction of Characteristic Signal

Various methods can be employed for the extraction of amplitude and phase information from characteristic signals, such as wavelet transform, fast Fourier transform (FFT), and Prony algorithm. It is important to note that while the wavelet transform is adept at capturing the time-varying characteristics of a signal, it is unable to provide accurate phase information. The FFT reduces the computational complexity by minimizing multiplication and addition operations. However, its accuracy is constrained by the factors such as the chosen data window and the presence of DC attenuation components. The Prony algorithm offers the ability to describe transient characteristics of a signal, encompassing amplitude, phase, and frequency information.

The Prony algorithm is adopted in this study to extract the amplitude and phase of characteristic signals based on the measurement data obtained from both ends of the transmission line. The mathematical model is expressed as:

y (t) = \sum_{i = 1}^{q} A_{i} e^{α_{i} t} c o s (2 π f_{i} t + θ_{i})

(5)

where $y (t)$ is the actual measured signal; A_i is the amplitude; $α_{i}$ is the attenuation factor; f_i is the frequency; θ_i is the phase; and q is number of orders. The calculation amount of the Prony algorithm will increase exponentially when the order is high. In order to address this issue, the measured data are pre-filtered using a 4^th-order Butterworth band-pass filter. This filtering process serves to simplify the computational workload and also eliminates any interference caused by noise.

According to the extraction length of the injected characteristic signal, the estimation accuracy and calculation complexity should be considered. It is recommended that the extracted data length be at least one period to avoid errors in the estimation result. However, it is not advisable to extend this length excessively, as this will increase the complexity of the calculation and consequently the computational speed. Therefore, the extraction length is set to be three periods of the characteristic signal.

III. Identification Method of Frequency-dependent Parameters in Flexible DC Grids

A. Calculation of Characteristic Impedance and Attenuation Function

The traditional method for identifying transmission line parameters in AC grids encounters challenges when applied to DC grids due to the operating frequency in the DC grid being 0 Hz. In this subsection, a novel method for identifying frequency-dependent parameters of flexible DC lines is proposed, which considers the distributed capacitance of the transmission line. This method effectively avoids the interference caused by high-frequency signals, enabling accurate identification of frequency-dependent parameters.

The mathematical model of the transmission line with distributed parameters can be expressed as:

[\begin{array}{l} U (0, s) \\ I (0, s) \end{array}] = [\begin{matrix} c o s h (γ (s) l) & Z_{c} (s) s i n h (γ (s) l) \\ \frac{s i n h (γ (s) l)}{Z_{c} (s)} & c o s h (γ (s) l) \end{matrix}] [\begin{array}{l} U (l, s) \\ I (l, s) \end{array}]

(6)

Z_{c} (s) = \sqrt[]{(R_{0} + s L_{0}) / (G_{0} + s C_{0})}

(7)

γ (s) = \sqrt[]{(R_{0} + s L_{0}) (G_{0} + s C_{0})}

(8)

where U(0,s) and I(0,s) are the voltage and current at the beginning of a transmission line, respectively; U(l,s) and I(l,s) are the voltage and current at the end of a transmission line, respectively; l is the length of a transmission line; Z_c(s) is the characteristic impedance; $γ (s)$ is the propagation function; and R₀, L₀, G₀, and C₀ are the resistance, inductance, conductance, and capacitance per unit length, respectively.

Equation (6) can be rewritten as:

[\begin{array}{l} U (l, s) \\ I (l, s) \end{array}] = [\begin{matrix} c o s h (γ (s) l) & - Z_{c} (s) s i n h (γ (s) l) \\ - \frac{s i n h (γ (s) l)}{Z_{c} (s)} & c o s h (γ (s) l) \end{matrix}] [\begin{array}{l} U (0, s) \\ I (0, s) \end{array}]

(9)

According to (6) and (9), we can obtain:

\begin{array}{l} U (0, s) - U (l, s) = (U (l, s) - U (0, s)) c o s h (γ (s) l) + \\ Z_{c} s i n h (γ (s) l) (I (l, s) + I (0, s)) \end{array}

(10)

\begin{array}{l} I (0, s) - I (l, s) = (I (l, s) - I (0, s)) c o s h (γ (s) l) + \\ \frac{U (l, s) + U (0, s)}{Z_{c} (s)} s i n h (γ (s) l) \end{array}

(11)

U (0, s) + U (l, s) \frac{I (l, s)}{I (0, s)} = (U (l, s) + U (0, s) \frac{I (l, s)}{I (0, s)}) c o s h (γ (s) l)

(12)

By analyzing (10)-(12), the characteristic impedance, propagation function, and attenuation function $A (s)$ at a certain frequency can be calculated as:

Z_{c} (s) = \sqrt[]{\frac{(U (0, s) + U (l, s)) (U (0, s) - U (l, s))}{(I (0, s) + I (l, s)) (I (0, s) - I (l, s))}}

(13)

γ (s) = \frac{1}{l} a r c o s h (\frac{U (0, s) I (0, s) + U (l, s) I (l, s)}{U (l, s) I (0, s) + U (0, s) I (l, s)})

(14)

A (s) = e^{- γ (s) l}

(15)

Since the transmission line and the earth will produce the high-frequency skin effect under the influence of alternating electromagnetic fields, the distribution parameters of the transmission line will change with frequency. As a result, the characteristic impedance and propagation constant become frequency-dependent parameters.

It is a network synthesis problem to chose a linear network that has the frequency characteristics of Z_c(s). Z_c(s) can be expressed by the RLC element in the form of the sum of partial fractions to express the characteristic impedance. In the complex frequency domain, a characteristic impedance function can be expressed as:

Z_{c} (s) = H \frac{(s + z_{1}) (s + z_{2}) \dots (s + z_{j})}{(s + p_{1}) (s + p_{2}) \dots (s + p_{j})}

(16)

where z_j and p_j are the zeros and the poles of the characteristic impedance, respectively; j is the number of the zeros and the poles; and $H = \sqrt[]{L_{0} / C_{0}}$ is the value of the characteristic impedance when the frequency is infinite. Equation (16) can be written as a sum of partial fractions:

Z_{c} (s) = k_{0} + \frac{k_{1}}{s + p_{1}} + \dots + \frac{k_{j}}{s + p_{j}}

(17)

where k_j is the residue; and $k_{0} = \underset{s \to \infty}{l i m} Z_{c} (s) = H$ is a constant. The constant k₀ directly corresponds to a fixed resistance R_s₀, while the remaining terms represent complex frequency-domain functions of the parallel resistance-capacitance loop. Then, the $i^{t h}$ circuit can be expressed as:

\frac{k_{i}}{s + p_{i}} \Rightarrow \frac{R_{s i} \frac{1}{s C_{s i}}}{R_{s i} + \frac{1}{s C_{s i}}} = \frac{\frac{1}{C_{s i}}}{s + \frac{1}{R_{s i} C_{s i}}}

(18)

where R_si and C_si ( $i = 1,2, \dots, j$ ) are the resistors and capacitors of each item in the component fraction, respectively.

\{\begin{array}{l} R_{s 0} = k_{0} = Z_{c} (s) |_{s \to \infty} \\ R_{s i} = \frac{1}{p_{i} C_{s i}} = \frac{k_{i}}{p_{i}} i = 1,2, \dots, j \\ C_{s i} = \frac{1}{k_{i}} i = 1,2, \dots, j \end{array}

(19)

The characteristic impedance can be represented as an equivalent network by using resistors and capacitors that can satisfy the boundary conditions of the characteristic impedance:

1) At low frequencies, the capacitor branch is approximately open-circuit, and the characteristic impedance is represented by a series of resistors.

2) At high frequencies, the capacitor branch is approximately short-circuit, and the characteristic impedance is R_s₀.

\{\begin{array}{l} Z_{c} (0) = R_{s 0} + R_{s 1} + R_{s 2} + \dots + R_{s j} \\ Z_{c} (s) |_{s \to \infty} = R_{s 0} \end{array}

(20)

From the above analysis, the characteristic impedance value presents a monotonically decreasing trend with the increase in frequency. Therefore, the characteristic impedance at different frequencies is calculated by injecting characteristic signals of different frequencies several times, and the characteristic impedance in the full frequency domain is fitted by the linear interpolation method. The analytical method for the propagation constant is the same as that for the characteristic impedance.

B. Pole and Residue Identification of Frequency-dependent Parameters by Vector Fitting Method

Through the calculation in the previous subsection, we can obtain the amplitude-frequency characteristic curve of characteristic impedance and attenuation function. For unknown coefficients, (16) is a nonlinear equation, which makes the problem more complex. The vector fitting method is considered an effective technique for fitting data. By locating poles to make them become known, (17) is linearized, and then the unknown quantity can be easily calculated by the least squares method.

The approximate formula for the rational function is:

f (s) \approx \sum_{i = 1}^{j} \frac{c_{i}}{s - a_{i}} + d + s h

(21)

where c_i and a_i are the residues and the poles of the rational function, respectively; and d and h are the parameters. The poles and residues in the context of the vector fitting method are either real or complex conjugates, while the parameters d and h are of real values.

In (21), a set of initial poles ${\bar{a}}_{i}$ is specified, and the function $f (s)$ is multiplied by an unknown function $σ (s)$ .

[\begin{matrix} σ (s) f (s) \\ σ (s) \end{matrix}] \approx [\begin{matrix} \sum_{i = 1}^{j} \frac{c_{i}}{s - {\bar{a}}_{i}} + d + s h \\ \sum_{i = 1}^{j} \frac{{\tilde{c}}_{i}}{s - {\bar{a}}_{i}} + 1 \end{matrix}]

(22)

where ${\tilde{c}}_{i}$ is the residue of the function $σ (s)$ ; and $c_{i}$ , ${\tilde{c}}_{i}$ , d, and h are unknown.

In (22), the rational approximation for $σ (s)$ shares the same poles as the approximation for $σ (s) f (s)$ . It is worth noting that the uncertainty in the solution for $σ (s)$ has been eliminated due to the approximation of $σ (s)$ being approximately equal to 1 at high frequencies. By multiplying the second row in (22) with $f (s)$ , we obtain the relationship as:

\sum_{i = 1}^{j} \frac{c_{i}}{s - {\bar{a}}_{i}} + d + s h \approx (\sum_{i = 1}^{j} \frac{{\tilde{c}}_{i}}{s - {\bar{a}}_{i}} + 1) f (s)

(23)

Formula (23) is linear, which can be rewritten as:

(\sum_{i = 1}^{j} \frac{c_{i}}{s - {\bar{a}}_{i}} + d + s h) - (\sum_{i = 1}^{j} \frac{{\tilde{c}}_{i}}{s - {\bar{a}}_{i}}) f (s) \approx f (s)

(24)

When (24) is rewritten for multiple frequency points, it leads to a linear problem:

A x = b

(25)

A = [\frac{1}{s - a_{1}} \frac{1}{s - a_{2}} \dots \frac{1}{s - a_{j}} 1 s \frac{- f (s)}{s - a_{1}} \frac{- f (s)}{s - a_{2}} \dots \frac{- f (s)}{s - a_{j}}]

(26)

x = {[c_{1} c_{2} \dots c_{j} d h {\tilde{c}}_{1} {\tilde{c}}_{2} \dots {\tilde{c}}_{j}]}^{T}

(27)

b = f (s)

(28)

Equation (25) can be solved as a least squares problem. In order to make the result converge, we need to take the resulting new poles as the initial poles and repeat the calculation.

IV. Improved DSEBP

The DSEBP utilizes the concept of model matching to identify faults through the comparison between measurements and estimated states [

8]. The key aspect of this method lies in the development of a dynamic model that accounts for the frequency-dependent characteristic of transmission line parameters. However, acquiring the frequency-dependent parameters of transmission lines proves to be challenging in actual projects. Therefore, the improved DSEBP is proposed, which establishes a new dynamic model using the poles and residues calculated in Section III to improve the engineering practicability. Meanwhile, considering that the ground-mode component passes through the earth and its loop parameters are affected by complex factors such as earth resistivity, only the transmission line-mode parameters are used for the subsequent fault identification.

The relationship among the forward wave F(s), the backward wave B(s), the voltage measurements V(s), and the current measurements I(s) can be expressed as:

\{\begin{array}{l} F_{w} (s) = V_{w} (s) + Z_{c} (s) I_{w} (s) \\ B_{w} (s) = V_{w} (s) - Z_{c} (s) I_{w} (s) \end{array}

(29)

\{\begin{array}{l} B_{k} (s) = A (s) F_{m} (s) \\ B_{m} (s) = A (s) F_{k} (s) \end{array}

(30)

V_{k} (s) - Z_{c} (s) I_{k} (s) = A (s) F_{m} (s)

(31)

where $w = k, m$ , which represent the names of both ends of the transmission line.

From (31), we can observe that the accurate characteristic impedance and attenuation function are the key to the subsequent calculations. The characteristic impedance amplitude exhibits a rapid decrease in the low-frequency range, followed by a relatively stable behavior in the high-frequency range. While the attenuation function varies very little in the frequency range considered in this study. Therefore, by making A(s) approximately equal to 1, (31) will change as:

V_{k} (s) - Z_{c} (s) I_{k} (s) = V_{m} (s) + Z_{c} (s) I_{m} (s)

(32)

To convert (32) from the frequency domain to the time domain, the recursive convolution theorem is employed. The improved dynamic model of the DC line in the time domain can be obtained as:

\begin{array}{l} V_{k} (t) = V_{m} (t) + I_{k z} (t) + I_{m z} (t) + I_{k} (t) R_{s 0} + I_{m} (t) R_{s 0} = V_{m} (t) + \\ g_{1} I_{k z} (t - Δ t) + c_{1} I_{k} (t) + q_{1} I_{k} (t - Δ t) + g_{1} I_{m z} (t - Δ t) + \\ c_{1} I_{m} (t) + q_{1} I_{m} (t - Δ t) + I_{k} (t) R_{s 0} + I_{m} (t) R_{s 0} \end{array}

(33)

where g₁, c₁, and q₁ are the constants that can be calculated by recursive convolution theorem; and $I_{k z} (t)$ and $I_{m z} (t)$ are the intermediate state variables.

Combined with the improved DC line model, the system measurement equation containing the measurements, state variables, and historical values is established.

The measurements can be divided into two categories: actual measurements and virtual measurements. Actual measurements include the voltage and current values $V_{k} (t)$ , $V_{m} (t)$ , $I_{k} (t)$ , and $I_{m} (t)$ . Virtual measurements represent a physical law that must be satisfied. Note that the value of the virtual measurement (zero) is known with certainty.

The state variables include the voltage state variables $V_{v, k} (t)$ and $V_{v, m} (t)$ , current state variable $I_{v, m} (t)$ , and intermediate state variables $I_{k z} (t)$ and $I_{m z} (t)$ . The historical values include the voltage and current values before time $Δ t$ $I_{k} (t - Δ t)$ , $I_{m} (t - Δ t)$ , $I_{k z} (t - Δ t)$ , and $I_{m z} (t - Δ t)$ . The system measurement equations can be established as:

\begin{array}{l} \{\begin{array}{l} I_{k} (t) = \frac{V_{k} (t) - V_{m} (t) - I_{k z} (t) - I_{m z} (t) - I_{m} (t) R_{s 0}}{R_{s 0}} \\ 0 = I_{k z} (t) - g_{1} I_{k z} (t - Δ t) - c_{1} I_{k} (t) - q_{1} I_{k} (t - Δ t) \\ 0 = I_{m z} (t) - g_{1} I_{m z} (t - Δ t) - c_{1} I_{m} (t) - q_{1} I_{m} (t - Δ t) \\ I_{m} (t) = I_{v, m} (t) \\ V_{k} (t) = V_{v, k} (t) \\ V_{m} (t) = V_{v, m} (t) \end{array} \Rightarrow \\ z (t) = f (x (t), t) + e_{r} (t) = H x (t) + C (t) + e_{r} (t) \end{array}

(34)

where z(t) is the column vector of measurements; x(t) is the column vector of state variables; $f (x (t), t)$ is the function related to the state variables; e_r(t) is the column vector of measurement error; H is a Jacobian matrix; and C(t) is the column vector of historical values.

The concept of dynamic state estimation is integrated into the protection relay by determining if the measurements can match the established model. The specific implementation method is weighted least squares.

m i n J = (z (t) - f (x {(t), t))}^{T} W (z (t) - f (x (t), t))

(35)

where $W = d i a g {1 / σ_{1}^{2}, 1 / σ_{2}^{2}, \dots, 1 / σ_{n u m}^{2}}$ is a weight matrix, and $σ_{i}$ is the standard deviation, and num is the number of measurements.

When line parameters are known, $f (x (t), t)$ is a linear function and the Jacobi matrix is a constant matrix. Therefore, the column vector of the optimal estimated state variables is obtained as:

\hat{x} (t) = (H^{T} {W H)}^{- 1} H^{T} W (z (t) - C (t))

(36)

The column vector of the estimated measurements $\hat{z} (t)$ and the sum of square of the normalized residuals $ζ$ can be calculated as (37) and (38). And $ζ$ is applied to identify whether an internal fault occurs.

\hat{z} (t) = f (\hat{x} (t), t) = H \hat{x} (t) + C (t)

(37)

ζ = \sum_{i = 1}^{n u m} {(\frac{{\hat{z}}_{i} (t) - z_{i} (t)}{σ_{i}})}^{2}

(38)

where $z_{i}$ is the i^th measurement; and ${\hat{z}}_{i}$ is the i^th estimated measurement.

The deviation between measurements and estimated measurements is defined as the residual r, which is expressed as:

r (t) = \hat{z} (t) - z (t)

(39)

If there are no internal faults, the measurements match the established model and satisfy the chi-square distribution. However, during internal faults, the measurements will no longer match the established model, and $ζ$ will not meet the chi-square distribution. Therefore, $ζ$ can be employed to determine whether the line is healthy. The protection criterion of improved DSEBP $T r i p (t)$ is:

T r i p (t) = \{\begin{array}{l} 1 \int_{t - T_{s e t}}^{t} S (τ) d τ \geq T_{d e l a y} & & S l o p e_k (t) < 0 \\ 0 o t h e r w i s e \end{array}

(40)

S (t) = \{\begin{array}{l} 1 ζ \geq χ_{γ}^{2} (F) \\ 0 ζ < χ_{γ}^{2} (F) \end{array}

(41)

S l o p e_k (t) = \frac{I_{n} (t) - I_{n} (t - 10)}{{\hat{I}}_{n} (t) - {\hat{I}}_{n} (t - 10)}

(42)

where T_set is the reset time; $χ_{γ}^{2} (F)$ is the threshold for $ζ$ ; T_delay is the trip delay; ${\hat{I}}_{n} (t)$ is the estimated value of $I_{n} (t)$ ; and Slope_k(t) is the ratio of the slope of $I_{n} (t)$ and ${\hat{I}}_{n} (t)$ from time $t - 10$ to time $t$ .

In this paper, $χ_{γ}^{2} (F)$ is set to be 6.63. When $ζ$ is greater than 6.63, it is known that there is a 99% probability that an internal fault occurs on the protected line by consulting the chi-square distribution table. Nevertheless, the aforementioned probability does not indicate that there is a 1% chance that the protection will malfunction. This is because the protection criterion contains a user-defined trip delay. When $ζ$ remains greater than the $χ_{γ}^{2} (F)$ for T_delay, the protection relay will issue a trip signal.

Meanwhile, Slope_k indicates the ratio of the slope of estimated measurements to the slope of actual measurements. When the external fault occurs, the estimated measurements are basically consistent with the actual measurements, and Slope_k is greater than 0. But when the internal fault occurs, there is a notable discrepancy between the estimated measurements and the actual measurements, and Slope_k is less than 0.

Therefore, by combining $ζ$ and Slope_k, the improved DSEBP can identify the faults correctly. Considering the reliability and sensitivity, the reset time T_set is set to be 1 ms, and the trip delay T_delay is set to be 5. The standard deviation of actual measurements and virtual measurements are assigned values of 0.02 p.u. and 0.001 p.u., respectively. The flow chart of the proposed improved DSEBP is shown in Fig. 4, which can be divided into four key steps.

Fig. 4 Flow chart of proposed method.

Step 1: determine whether the line parameters need to be calculated. If the protection method does not require line parameter identification, it can proceed directly to Step 3. However, if the protection method is used for the first time, the line parameters must be calculated.

Alternatively, the protection method needs to update the line parameters to improve the protection performance, and then online identification of the line parameters is also required.

The control strategy of HB-MMC is changed from the conventional control strategy to the characteristic signal injection control strategy. Afterwards, the characteristic signals are injected by controlling the addition and removal of additional submodules $Δ N$ . And the amplitude and phase of the characteristic signals are extracted by the Prony method.

Step 2: using the amplitude and phase of the characteristic signals at different frequencies, the corresponding characteristic impedance Z_c(s) and attenuation function A(s) can be calculated according to (13)-(15). Subsequently, the poles and residues of the frequency-dependent parameters can be calculated by the vector fitting method.

Step 3: the system measurement equation (34) can be established based on the characteristic impedance, attenuation function, measurements, state variables, and historical values. And the weighted least squares method is utilized to calculate the optimal estimated state variables $\hat{x} (t)$ , the estimated measurements $\hat{z} (t)$ , and the normalized residuals $ζ$ .

Step 4: by combining the $ζ$ and $S l o p e_k$ in the protection criterion, the improved DSEBP can identify the faults correctly.

V. Simulation and Analysis

The ±500 kV four-terminal flexible DC transmission system shown in Fig. 5 is constructed by PSCAD/EMTDC.

Fig. 5 Four-terminal flexible DC transmission system.

The characteristic signal is injected into the flexible DC system through MMC1. The green squares are the location of the protection devices. The protection sampling rate is set to be 10 kHz. The main parameters of the simulation model are shown in Table I.

TABLE I Main Parameters of Simulation Model

Model	Parameter	Value
Transmission line	Length of line I (MMC1-MMC2)	200 km
	Length of line II (MMC1-MMC3)	50 km
	Length of line III (MMC2-MMC4)	190 km
	Length of line IV (MMC3-MMC4)	220 km
HB-MMC	Number of units	200
	Bridge inductance	100 mH
	Submodule capacitance	10 mF
	Rated capacity	1500 MVA

A. Identification Accuracy Analysis of Line Parameters

1)　Length of Line I is 200 km

To verify the accuracy of the parameter identification method proposed in this paper, the line parameters are calculated for line I. The parameter identification results of the line characteristic impedance amplitudes and phase angles at different frequencies are shown in Table II and Table III.

TABLE II Parameter Identification Results of Line Characteristic Impedance Amplitudes

Frequency (Hz)	Accurate value (Ω)	Fitted value (Ω)	Relative error (%)
1	337.845	338.885	<1
10	225.456	225.459	<1
100	219.172	219.196	<1
1000	216.241	216.245	<1
4000	216.194	216.220	<1

TABLE III Parameter Identification Results of Line Characteristic Impedance Phase Angles

Frequency (Hz)	Accurate value (rad)	Fitted value (rad)	Relative error (%)
1	-0.507100	-0.504200	<1
10	-0.098400	-0.098300	<1
100	-0.025800	-0.025700	<1
1000	-0.003699	-0.003691	<1
4000	-0.000929	-0.000933	<1

It can be observed that the parameter identification method proposed in this paper exhibits excellent accuracy in the calculation of the line characteristic impedance amplitudes and phase angles across a wide range of frequencies, from low to high. The relative error can be controlled to within 1% for both the calculation of line characteristic impedance amplitudes and phase angles.

2)　Length of Line I is 1000 km

As shown in Tables II and III, the calculation of the line characteristic impedance amplitudes and phase angles exhibits high accuracy with a length of transmission line of 200 km. It is essential to ascertain whether the proposed parameter identification method is applicable to transmission lines of greater length.

From (7), it can be observed that the line characteristic impedance is only related to the line parameters. Therefore, the length of the transmission line does not affect the line characteristic impedance magnitudes and phase angles.

In order to verify the accuracy of the aforementioned analysis, the line parameters are identified with the length of line I set to be 1000 km. The parameter identification results of the line characteristic impedance amplitudes and phase angles at different frequencies are shown in Tables IV and V.

TABLE IV Parameter Identification Results of Line Characteristic Impedance Amplitudes when Length of Transmission Line is 1000 km

Frequency (Hz)	Accurate value (Ω)	Fitted value (Ω)	Relative error (%)
1	337.845	338.847	<1
10	225.456	225.416	<1
100	219.172	219.224	<1
1000	216.241	216.219	<1
4000	216.194	216.135	<1

TABLE V Parameter Identification Results of Line Characteristic Impedance Phase Angles when Length of Transmission Line is 1000 km

Frequency (Hz)	Accurate value (rad)	Fitted value (rad)	Relative error (%)
1	-0.507100	-0.504600	<1
10	-0.098400	-0.098100	<1
100	-0.025800	-0.025700	<1
1000	-0.003699	-0.003695	<1
4000	-0.000929	-0.000932	<1

As evidenced by Tables IV and V, with an increased length of transmission line, the proposed parameter identification method also exhibits excellent accuracy in the calculation of the line characteristic impedance amplitudes and phase angles across a wide range of frequencies as well. Furthermore, it can be verified that the length of the transmission line does not affect the line characteristic impedance magnitudes and phase angles.

B. Identification Accuracy Analysis of Poles and Residues

The identification results of poles and residues of the line characteristic impedance are shown in Table VI. There are some differences between the fitted values of poles and residues and the accurate values. This is because the fitted values of poles and residues only consider the line characteristic impedance under 5000 Hz. However, the accurate values of poles and residues in PSCAD consider the line characteristic impedance under 100000 Hz.

TABLE VI Identification Results of Poles and Residues of Line Characteristic Impedance

Category	Accurate value	Fitted value
Poles	-1.16143	-1.04139
	-6.11132	-4.91769
	-619.66040	-301.97479
Residues	693.90621	592.66130
	688.34288	728.90311
	2849.97961	2289.83266
Constant	216.58498	216.49197

The fitted results of the line characteristic impedance amplitudes and phase angles are shown in Fig. 6. The amplitude-frequency curve and phase angle-frequency curve are basically consistent with the accurate curves in PSCAD. The relative errors of the parameters with different frequency bands are less than 1%, which demonstrates that the proposed parameter identification method has good fitting accuracy.

Fig. 6 Fitted results of line characteristic impedance amplitudes and phase angles. (a) Amplitudes. (b) Phase angles.

C. Accuracy of Established Dynamic Model

1)　Length of Line I is 200 km

The simulation results under the normal operation condition of the power system with the length of transmission line of 200 km, as shown in Fig. 7, serve as a means to evaluate the accuracy of the established dynamic model. Both the estimated and measured current/voltage curves exhibit a high level of consistency, with a relative error of less than 1%. Therefore, the established dynamic model is accurate enough to be used in the improved DSEBP.

Fig. 7 Simulation results under normal operation condition of power system with length of transmission line of 200 km. (a) Voltage curves. (b) Current curves.

2)　Length of Line I is 1000 km

In Section IV, it is assumed that the propagation function is equal to 1 to reduce the computational complexity. From Fig. 7, it can be observed that the established dynamic model can still maintain a high accuracy after simplifying the propagation function for a transmission line with a length of 200 km. It is necessary to find out if this simplification can be valid for longer transmission lines.

From (8) and (15), it can be observed that the attenuation function is affected by the line parameters and the length of transmission line. At the same frequency, the attenuation function of a long transmission line is less than that of a short transmission line. Therefore, this simplification applied to long transmission lines may introduce more modeling errors.

The simulation results under normal operation condition of power system with the length of transmission line of 1000 km are shown in Fig. 8.

Fig. 8 Simulation results under normal operation condition of power system with length of transmission line of 1000 km. (a) Voltage curves. (b) Current curves.

As illustrated in Fig. 8, when the length of transmission line is 1000 km, there are discrepancies between the measured and estimated current curves due to the simplification. However, both the estimated and measured current/voltage curves exhibit a high level of consistency, with a relative error of less than 1.5%. Therefore, the simplification can also be applied to transmission lines with a length of 1000 km.

D. Performance of Improved DSEBP

1)　Influence of Fault Locations and Types

According to Fig. 5, the fault locations are designated as f1, f2, f3 and f4. Notably, f4 is situated at distance of 0% along line III, starting from the m side. The remaining three fault points (f1, f2, and f3) are situated at distances of 0%, 50%, and 100% along the line I, starting from the k side. To validate the efficacy of the improved DSEBP, diverse fault types are simulated, including positive pole-to-ground (P-PTG) faults, negative pole-to-ground (N-PTG) faults, and pole-to-pole (PTP) faults.

The simulation results of protection performance under different fault locations and types are shown in Fig. 9.

Fig. 9 Simulation results of protection performance under different fault locations and types. (a) f2, P-PTG. (b) f2, PTP. (c) f4, PTP.

As shown in Fig. 9(a) and (b), before an internal fault occurs, the estimated current curve is consistent with the measured current curve and $ζ$ is nearly zero. However, there are significant differences between the estimated and measured current curves during the internal fault. According to the protection criterion of the improved DSEBP, $ζ$ exceeds $χ_{γ}^{2} (F)$ , and Slope_k is -7.8 at $f 2$ under P-PTG faults and -7.1 at f2 under PTP faults, respectively. Therefore, the improved DSEBP correctly sends out the trip signal.

As shown in Fig. 9(c), the estimated current curve is consistent with the measured current curve before an external fault occurs. During the external fault, the estimated and measured current curves remain consistent. According to the protection criterion of the improved DSEBP, $ζ$ exceeds the threshold value, but Slope_k is 0.1 at f4 under PTP faults, which does not meet the protection criterion. Therefore, the improved DSEBP can correctly identify the external fault and will not operate incorrectly.

The performance of the improved DSEBP under different fault locations and types is shown in Table VII. It shows that the improved DSEBP can quickly identify the internal faults and will not operate incorrectly during an external fault.

TABLE VII Simulation Results Under Different Fault Locations and Types

Fault location	Fault type	Behavior of protection
Fault location	Fault type	k side	m side
f1	P-PTG	Trip at 2.1 ms	Trip at 2.7 ms
	N-PTG	Trip at 2.2 ms	Trip at 2.7 ms
	PTP	Trip at 2.4 ms	Trip at 1.6 ms
f2	P-PTG	Trip at 2.4 ms	Trip at 2.4 ms
	N-PTG	Trip at 2.3 ms	Trip at 2.4 ms
	PTP	Trip at 2.0 ms	Trip at 2.0 ms
f3	P-PTG	Trip at 2.6 ms	Trip at 2.0 ms
	N-PTG	Trip at 2.7 ms	Trip at 2.1 ms
	PTP	Trip at 1.6 ms	Trip at 2.3 ms
f4	P-PTG
	N-PTG
	PTP

2)　Influence of Fault Resistance

To test the ability of the improved DSEBP to withstand fault resistance, an internal P-PTG fault (f2) with a fault resistance of 300 Ω is set and the simulation results are shown in Fig. 10.

Fig. 10 Simulation results of protection performance under internal P-PTG fault (f2) with fault resistance of 300 Ω.

It shows that there are significant differences between the estimated and the measured current curves during the internal fault.

According to the protection criterion of the improved DSEBP and fault phase selection criterion, $ζ$ exceeds $χ_{γ}^{2} (F)$ , and Slope_k is -8.45. Therefore, the improved DSEBP can operate normally with a fault resistance of 300 $Ω$ .

3)　Influence of Noise Disturbance

The ability of the improved DSEBP to withstand noise disturbance is also tested. The 25 dB noise is added to the measurement data to verify the protection performance.

As shown in Table VIII, the improved DSEBP can operate normally under the 25 dB noise disturbance.

TABLE VIII Simulation Results Under Noise Disturbance

Fault location	Fault type	$ζ$	Slope_k	Fault identification
f1	P-PTG	37.6	-6.40	Internal fault
	N-PTG	58.3	-5.50	Internal fault
	PTP	85.9	-8.10	Internal fault
f2	P-PTG	31.4	-11.40	Internal fault
	N-PTG	37.6	-6.30	Internal fault
	PTP	52.6	-8.00	Internal fault
f3	P-PTG	42.8	-5.80	Internal fault
	N-PTG	44.1	-7.20	Internal fault
	PTP	39.8	-8.80	Internal fault
f4	P-PTG	16.6	0.60	External fault
	N-PTG	24.2	0.90	External fault
	PTP	32.8	0.06	External fault

4)　Influence of Synchronization Error

For double-ended protection methods, the synchronization error may affect the correct fault identification. A synchronization error of 1 ms is added to the PTP fault (f4), and the simulation result is given in Fig. 11.

Fig. 11 Simulation results of protection performance with synchronization error of 1 ms.

Due to the synchronization error, there is a certain deviation between the estimated and measured current curves during the external fault. However, the trend and amplitude of the current curve remain consistent. Meanwhile, Slope_k is 0.02, which does not satisfy the protection criterion. Therefore, the simulation results verify that the improved DSEBP can withstand the effect of the synchronization error.

E. Reliability Verification of Slope_k

The trip signal is activated only when the slopes of the measured and estimated currents are opposite. It is important to determine if the slopes of the measured and estimated currents are always opposite when an internal fault occurs.

From (34), it can be observed that the system measurement equation in the time domain is more complex and it contains numerous state variables. Therefore, we qualitatively analyze this problem from the frequency domain based on (32). The current on the k side of the line in (32) can be expressed as:

\begin{array}{l} I_{k} (s) = \frac{V_{k} (s) - V_{m} (s) - Z_{c} (s) I_{m} (s)}{Z_{c} (s)} = \\ \frac{V_{k} (s) - V_{m} (s)}{Z_{c} (s)} - I_{m} (s) \approx - I_{m} (s) \end{array}

(41)

Formula (41) is correct when the system is operating normally. If the transmission line voltage drop is ignored, then $I_{k} (s) \approx - I_{m} (s)$ . For the measured current I_k(s), its direction is the same before and after the internal fault. However, when an internal fault occurs, the physical laws of the transmission line will break down and (41) will not be able to correctly calculate the current at the k side. Two extreme scenarios of internal faults are taken for illustration.

1)　Scenario with Fault at f1

When the fault occurs at f1, the voltage value at the k side is almost zero and the current value at the k side is:

{\hat{I}}_{k, k} (s) = \frac{- V_{m} (s) + Z_{c} (s) I_{m} (s)}{Z_{c} (s)} = - \frac{V_{m} (s)}{Z_{c} (s)} + I_{m} (s)

(42)

The current value before the fault is $I_{k 0 -} (s) \approx - I_{m 0 -} (s)$ , where $I_{m 0 -}$ is the current value at the m side of the line before the fault occurs. If $S l o p e_k > 0$ is required, ${\hat{I}}_{k, k} (s)$ must be less than $- I_{m 0 -} (s)$ . As shown in (42), V_m(s) should be as large as possible and I_m(s) should be as small as possible. If the length of the transmission line increases, V_m(s) after the fault will become larger. If the fault resistance increases, I_m(s) after the fault will become smaller.

Therefore, the length of line I is set to be 1000 km, and the following two tests are used to verify the performance of the proposed improved DSEBP.

Test 1: the internal PTG fault occurs at f1, and the simulation results are shown in Fig. 12. As shown in Fig. 12, before an internal fault occurs, the estimated current curve is consistent with the measured current curve, and $ζ$ is nearly zero. However, there are significant differences between the estimated and measured current curves during the internal fault. And it is evident that the slopes between the measured current I_k(t) and the estimated current ${\hat{I}}_{k} (t)$ are opposite.

Fig. 12 Simulation results of protection performance in test 1.

According to the protection criterion of the improved DSEBP, $ζ$ exceeds $χ_{γ}^{2} (F)$ , and Slope_k is -4.6. Therefore, the improved DSEBP correctly sends out the trip signal.

Test 2: the PTG fault with 300 $Ω$ fault resistance occurs at f1, and the simulation results are shown in Fig. 13. As shown in Fig. 13, before an internal fault occurs, the estimated current curve is consistent with the measured current curve, and $ζ$ is nearly zero. However, there are differences between the estimated and measured current curves during the internal fault with high fault resistance. It is evident that the slopes between the measured current I_k(t) and the estimated current ${\hat{I}}_{k} (t)$ are opposite within a period of time after the fault occurs.

Fig. 13 Simulation results of protection performance in test 2.

According to the protection criterion of the improved DSEBP, $ζ$ exceeds $χ_{γ}^{2} (F)$ , and Slope_k is -2.4. Therefore, the improved DSEBP correctly sends out the trip signal.

2)　Scenario with Fault at f3

When the fault occurs at f3, the voltage value at the m side is almost zero and the current value at the m side is:

{\hat{I}}_{k, m} (s) = \frac{V_{k} (s) + Z_{c} (s) I_{m} (s)}{Z_{c} (s)} = \frac{V_{k} (s)}{Z_{c} (s)} + I_{m} (s)

(43)

The value of the current before the fault is $I_{k 0 -} (s) \approx - I_{m 0 -} (s)$ . It is evident that the slopes of ${\hat{I}}_{k, m} (s)$ and $- I_{m 0 -} (s)$ are opposite. Therefore, the slopes between the measured current I_k(t) and the estimated current ${\hat{I}}_{k} (t)$ are opposite.

In summary, the qualitative analysis and simulation results show that the slopes between the measured current I_k(t) and the estimated current ${\hat{I}}_{k} (t)$ are opposite after the internal fault occurs. The proposed improved DSEBP can be applied to transmission lines with a length of 200-1000 km, which can meet the requirements of the majority of practical projects.

VI. Conclusion

The operating frequency of the DC system during steady-state operation is 0 Hz. Obtaining accurate transmission line parameters utilized by the protection is a challenging endeavor. Therefore, a method to identify the frequency-dependent parameters of flexible DC grids based on the HB-MMC is proposed. And the improved DSEBP can operate normally based on the identified parameters. The conclusions are as follows.

1) The frequency-dependent parameter can be identified with high accuracy by injecting characteristic signals into the HB-MMC. In the context of protection methods that rely on transmission line parameters, the practicality is enhanced.

2) The available residues and poles of the frequency-dependent parameter can be obtained by the vector fitting method, which can be used in the improved DSEBP to effectively identify faults.

3) The improved DSEBP can operate within 3 ms during an internal fault. And it can operate normally under a fault resistance of 300 $Ω$ , exhibiting a high degree of sensitivity. In addition, it has high reliability to withstand the 25 dB noise disturbance and 1 ms synchronization error.

References

J. He, K. Chen, M. Li et al., “Review of protection and fault handling for a flexible DC grid,” Protection and Control of Modern Power Systems, vol. 5, p. 15, Dec. 2020. [Baidu Scholar]

L. Liu, Z. Liu, M. Popov et al., “A fast protection of multi-terminal HVDC system based on transient signal detection,” IEEE Transactions on Power Delivery, vol. 36, no. 1, pp. 43-51, Feb. 2021. [Baidu Scholar]

T. Lan, Y. Li, and X. Duan, “High fault-resistance tolerable traveling wave protection for multi-terminal VSC-HVDC,” IEEE Transactions on Power Delivery, vol. 36, no. 2, pp. 943-956, Apr. 2021. [Baidu Scholar]

C. Zhang, G. Song, and X. Dong, “A novel traveling wave protection method for DC transmission lines using current fitting,” IEEE Transactions on Power Delivery, vol. 35, no. 6, pp. 2980-2991, Dec. 2020. [Baidu Scholar]

H. Xiao, S. Ren, and Y. Li, “Novel disturbance blocking criterion for reliable current differential protection of LCC-HVDC lines,” IEEE Transactions on Power Delivery, vol. 36, no. 1, pp. 477-480, Feb. 2021. [Baidu Scholar]

K. Chen, J. He, M. Li et al., “A similarity comparison based pilot protection scheme for VSC-HVDC grids considering fault current limiting strategy,” Journal of Modern Power Systems and Clean Energy, vol. 11, no. 4, pp. 1305-1315, Jul. 2023. [Baidu Scholar]

B. Li, M. Lv, B. Li et al., “Research on an improved protection principle based on differential voltage traveling wave for VSC-HVDC transmission lines,” IEEE Transactions on Power Delivery, vol. 35, no. 5, pp. 2319-2328, Oct. 2020. [Baidu Scholar]

J. He, M. Nie, M. Li et al., “Dynamic state estimation based protection for flexible DC grid,” IEEE Transactions on Industrial Electronics, vol. 70, no. 3, pp. 3069-3079, Mar. 2023. [Baidu Scholar]

R. Schulze, P. Schegner, and R. Zivanovic, “Parameter identification of unsymmetrical transmission lines using fault records obtained from protective relays,” IEEE Transactions on Power Delivery, vol. 26, no. 2, pp. 1265-1272, Apr. 2011. [Baidu Scholar]

Z. Liang, Z. Yang, P. Li, et al., “Study parameters power transmission line frequency measurement system,” Power System Technology, vol. 25, no. 3, pp. 34-37, Mar. 2021. [Baidu Scholar]

X. Yang, G. Geng, Y. Wang et al., “Method for accurately measuring the power-frequency parameters of EHV/UHV transmission lines,” IET Generation, Transmission & Distribution, vol. 11, no. 7, pp. 1726-1734, May 2017. [Baidu Scholar]

L. Hofmann, “Series expansions for line series impedances considering different specific resistances, magnetic permeabilities, and dielectric permittivities of conductors, air, and ground,” IEEE Transactions on Power Delivery, vol. 18, no. 2, pp. 564-570, Apr. 2003. [Baidu Scholar]

M. Asprou and E. Kyriakides, “Identification and estimation of erroneous transmission line parameters using PMU measurements,” IEEE Transactions on Power Delivery, vol. 32, no. 6, pp. 2510-2519, Dec. 2017. [Baidu Scholar]

A. Xue, F. Xu, J. Xu et al., “Online pattern recognition and data correction of PMU data under GPS spoofing attack,” Journal of Modern Power Systems and Clean Energy, vol. 8, no. 6, pp. 1240-1249, Dec. 2020. [Baidu Scholar]

H. Hong, H. Kong, L. Gu et al., “Steady-state PMU data selection for parameter identification of transmission line considering the influence of measurement error,” IET Generation, Transmission & Distribution, vol. 16, no. 22, pp. 4549-4562, Nov. 2022. [Baidu Scholar]

C. S. Indulkar and K. Ramalingam, “Estimation of transmission line parameters from measurements,” International Journal of Electrical Power & Energy Systems, vol. 30, no. 5, pp. 337-342, Jun. 2008. [Baidu Scholar]

S. Kurokawa, J. Pissolato, M. C. Tavares et al., “A new procedure to derive transmission-line parameters: applications and restrictions,” IEEE Transactions on Power Delivery, vol. 21, no. 1, pp. 492-498, Jan. 2006. [Baidu Scholar]

Y. Wang and W. Xu, “Algorithms and field experiences for estimating transmission line parameters based on fault record data,” IET Generation, Transmission & Distribution, vol. 9, no. 13, pp. 1773-1781, Oct. 2015. [Baidu Scholar]

T. Wang, G. Song, and K. S. T. Hussain, “Adaptive single-pole auto-reclosing scheme for hybrid MMC-HVDC systems,” IEEE Transactions on Power Delivery, vol. 34, no. 6, pp. 2194-2203, Dec. 2019. [Baidu Scholar]

M. Li, K. Chen, J. He et al., “Impedance phase based active injection protection for radial VSC-HVDC grid,” IEEE Transactions on Power Delivery, vol. 38, no. 3, pp. 1621-1635, Jun. 2023. [Baidu Scholar]

Address:No.19 Chengxin Avenue, Jiangning District, Nanjing 211106, China

E-mail: mpce@alljournals.cn

Tel:86-25-81093060

Fax:86-25-81093040

Home

Introduction

Editorial Board

For Author

Call For Papers

APC

Sponsor & Publisher