Optimization-based Settingless Algorithm Combining Protection and Fault Identification

TL;DR

This paper introduces a settingless, optimization-based protection algorithm for medium-voltage lines that quickly and accurately detects faults, identifies fault types, and estimates fault characteristics without relying on empirically set parameters.

Contribution

It proposes a novel settingless, convex optimization-based protection method that combines fault detection, location, and type identification in a unified framework.

Findings

Operates reliably with noisy measurements

Accurately estimates fault characteristics

Works quickly under various grid conditions

Abstract

The demand for faster protection algorithms is growing due to the increasingly faster dynamics in the system. The majority of existing algorithms require empirically selected set-points, which may reduce sensitivity to internal faults and cause security problems. This paper addresses these challenges by proposing a settingless time-domain unit protection algorithm for medium-voltage lines. The main idea of the algorithm is to identify which model of a protected line, i.e. healthy or with an internal fault, is more consistent with the input measurements. This is done by solving a number of small-scale convex optimization problems, which at the same time determine the characteristics of an internal fault that best fit the measurements. Thus, the proposed algorithm merges protection, fault location and fault type identification functionalities. The algorithm's performance is extensively tested on a grid model in MATLAB Simulink for different types of generation and grid operating conditions. The results demonstrate that the algorithm can operate quickly and reliably, and accurately estimate fault characteristics even in the presence of noisy measurements and uncertain line parameters.

Tables

Table 1. TABLE I: Parameters of test grid components

Protected Line				Grid
$L_L$, mH/km	$R_L$, $\Omega$/km	$K_L$, p.u.	$d$, km	$U_g$, kV	$f$, Hz
[1.3, 1.4]	[0.2, 0.42]	3	[10, 80]	36	50
Systems $S1$ and $S2$
$L_S$, mH	$R_S$, $\Omega$	$K_S$, p.u.	$E_1$, p.u.	$E_2$, p.u.	$\theta$, ∘
[46, 250]	[1.4, 19.4]	1.5	[0.9, 1.1]	[0.9, 1.1]	[-30, 30]

Table 2. TABLE II: Considered fault types and ranges of fault characteristics

Type	$R_a$, $\Omega$	$R_b$, $\Omega$	$R_c$, $\Omega$	$R_g$, $\Omega$	$T$, ms
$K_3$	$R_a=R_b=R_c\in$ [0, 100]			$\infty$	[5, 25]
$K_2$	$\left(R_a+R_b\right)\in$ [0, 100]		$\infty$	$\infty$	[5, 25]
$K_{2g}$	[0, 100]	[0, 100]	$\infty$	[0, 100]	[5, 25]
$K_1$	[0, 100]	$\infty$	$\infty$	0	[5, 25]

Table 3. TABLE III: Results for internal faults for base test cases (Section IV-D )

Group	Fault Type	Protection	Fault Characteristics
		% Faults Detected	Location Error, m		Resistance Error, $\Omega$
			max	mean	max	mean
2	$K_3$	100	6.66	0.78	9.13	0.009
	$K_2$	100	48.71	0.38	1.98	0.005
	$K_{2g}$	100	30.99	0.65	0.41	0.031
	$K_1$	100	6.57	0.54	0.09	0.006
3	$K_3$	100	4.27	0.59	0.29	0.002
	$K_2$	100	18.31	0.23	0.01	0.001
	$K_{2g}$	100	29.24	0.42	0.08	0.006
	$K_1$	100	2.72	0.29	0.01	0.001

Table 4. TABLE IV: Setup for test cases with converter-based generation (Section  IV-E )

Protected Line				$S1$
$L_L$, H/km	$R_L$, $\Omega$/km	$K_L$	$d$, km	$L_S$, H	$R_S$, $\Omega$	$K_S$	$E_1$	${\theta }_1$, ∘
0.0013	0.2	3	40	0.046	1.4	1.5	1	0
$S2$ (data for one wind turbine)
$P$, MW	$X_d$	$X_d^{\prime }$	$X_d^{\prime \prime }$	$X_q$	$X_q^{\prime \prime }$	$X_l$	$R_s$	$H$, s
2	1.305	0.296	0.252	0.474	0.243	0.18	0.006	0.62

Table 5. TABLE V: Results for internal faults for test cases with converter-based generation (Section  IV-E )

Group	Fault Type	Protection	Fault Characteristics
		% Faults Detected	Location Error, m		Resistance Error, $\Omega$
			max	mean	max	mean
2	$K_3$	100	30.83	5.79	1.45	0.036
	$K_2$	100	82.79	7.15	2.50	0.055
	$K_{2g}$	100	39.77	6.57	0.87	0.058
	$K_1$	100	21.43	6.69	0.78	0.021
3	$K_3$	100	17.16	4.45	0.94	0.024
	$K_2$	100	40.84	5.01	1.72	0.053
	$K_{2g}$	100	27.61	4.83	0.10	0.023
	$K_1$	100	14.55	5.69	0.08	0.011

Equations

$$u_1-u_2-\left(R\cdot i_1+L\cdot \frac{\partial i_1}{\partial t}\right)=0$$

$$u_2-u_1-\left(R\cdot i_2+L\cdot \frac{\partial i_2}{\partial t}\right)=0,$$

$$u_1-u_2-\alpha \cdot Z\cdot i_1+(1-\alpha )\cdot Z\cdot i_2=0$$

$$u_1-\alpha \cdot Z\cdot i_1-Z_F\cdot (i_1+i_2)=0,$$

$$Z_F:=\left[\begin{array}{ccc}{R}_a+R_g& R_g& R_g\\ R_g& R_b+R_g& R_g\\ R_g& R_g& R_c+R_g\end{array}\right].$$

$$S:=\left[\begin{array}{c}{U}_1-U_2-\left(R\cdot I_1+L\cdot I_1^{\prime }\right)\\ U_2-U_1-\left(R\cdot I_2+L\cdot I_2^{\prime }\right)\end{array}\right]$$

$$W(R_F,\alpha ):=\left[\begin{array}{c}{U}_1-U_2-\alpha \cdot Z\cdot I_1+(1-\alpha )\cdot Z\cdot I_2\\ U_1-\alpha \cdot Z\cdot I_1-Z_F\cdot (I_1+I_2)\end{array}\right],$$

$${\Delta }_1:=\frac{1}{6N}\sum _{i=1}^6\sum _{j=1}^N{S}_{i,j}^2,$$

$${\Delta }_2:=$$

$$0\le R_F\le R_F^{\max }$$

$$0\le \alpha \le 1,$$

$${\Delta }_2:=$$

$$0\le x\le x^{\max },$$

$${\mathcal{N}}_k$$

$${\mathcal{D}}_k$$

$${\mathcal{F}}_k$$

$${\Delta }_m:=$$

$$0\le R_F\le R_F^{\max }$$

$$0\le \alpha \le 1,$$

Full text

I Introduction

Deployment of distributed generation (DGs) and increasing power demand are pushing power systems to their limits and, as a consequence, are tightening the requirements for the fault clearing time [1]. In addition, converter-based DGs produce voltages and currents in case of a fault that may cause malfunctioning of conventional phasor-based relays. These problems necessitate the development and utilization of time-domain protection algorithms, which by their nature should operate well below one cycle of the power frequency and work reliably in grids with various types of DGs. Two main approaches have been proposed in the literature in this direction: traveling wave (TW) protection algorithms and protection algorithms based on differential equations for lumped circuits [2].

The first approach is based on detecting and analyzing traveling waves in a power grid with line models represented by circuits with distributed parameters. Since fault-generated TWs propagate close to the speed of light and contain valuable information about their source, TW protection algorithms can be very fast and sensitive to grid disturbances. Various methods have been proposed in the literature to detect a faulty line using TW theory. For instance, the authors in [3] propose a two-terminal TW protection that is based on the analysis of the arrival times of the first wavefronts of the disturbance-generated TWs. A differential protection based on equivalent TWs is presented in [4]. An extensive review of these and other TW protection algorithms can be found in [2] and [5], in which their advantages and disadvantages are highlighted. A common problem for a majority of the TW protection algorithms is that they may malfunction in case of close-in faults and faults with small inception angles.

The second approach is based on analyzing differential equations that describe a lumped-parameter model of a power grid. While it is naturally slower than the TW-based approach, it does not have the aforementioned problems of the TW protection algorithms. The authors in [6, 7, 8] propose protection algorithms based on incremental quantities, which are fault-generated voltages and currents. Particular combinations of the incremental quantities are compared to thresholds to classify internal and external faults. A summary of algorithms based on incremental quantities is presented in [2]. Taking it a step further, the authors in [9] propose a protection scheme that merges incremental quantity- and TW-based algorithms to improve the performance of time-domain protection. In [10], detailed guidelines for the calculation of set-points for this scheme are presented. Another protection principle that is based on dynamic state estimation is proposed in [11] and [12]. It utilizes measurements of terminal voltages and currents to estimate the goodness of fit of a healthy line model to the measurements. The goodness of fit obtained with the chi-square test is compared with a threshold to identify whether the protected line is healthy or not. While the proposed algorithm operates reliably in case of internal faults, transients due to certain external faults may lead to an erroneous tripping of a protected line. Therefore, a half-a-cycle delay is introduced before making a final tripping decision.

To the best of the authors’ knowledge, all proposed and conventional protection algorithms require set-points. The majority of these set-points are selected empirically, which may lead to reduced sensitivity to internal faults and security problems. The aim of this paper is to address this drawback by proposing an optimization-based settingless unit protection algorithm for medium-voltage power lines. This time-domain algorithm belongs to the group of approaches based on differential equations for lumped circuits. The main idea of the algorithm is to identify which model of a protected line, healthy or with an internal fault, is more consistent with the input measurements. This is done by solving a number of optimization problems, which at the same time determine the characteristics of an internal fault that best fit the measurements. The proposed algorithm has the following distinctive features and advantages over the existing algorithms:

• •

relies on instantaneous voltage and current measurements from both line ends,

• •

does not require selection of set-points, only parameters of the protected line have to be predefined,

• •

merges protection, fault location and fault type identification functionalities to unify these closely related functions and determine all key fault characteristics at once,

• •

is efficient and reliable since the optimization problems are small-scale and convex, hence they can be quickly solved to global optimality.

II Models of Protected Line

This work assumes that the protected line can be in two states: healthy or with an internal fault, which is represented by a generic model shown in Fig. 1. Both line ends are connected to the grid, a representation of which is omitted since the proposed algorithm does not require any grid information. A healthy line at any time instance can be described with the following differential equations in matrix form:

[TABLE]

where $R$ and $L$ are matrices of self and mutual resistances and inductances of phases, respectively, and $u_{1}$, $u_{2}$, $i_{1}$ and $i_{2}$ are vectors of three-phase instantaneous voltages and currents. Applying Kirchhoff’s laws to the line model with an internal fault yields the following equations:

[TABLE]

where $Z:=R+L\cdot\frac{\partial}{\partial t}$ and $\alpha$ is the fault location defined as relative distance from the left line terminal to the fault with respect to the line length. The matrix of fault resistances $Z_{F}$ is defined as follows:

[TABLE]

Equations (1)-(4) serve as the basis for the proposed algorithm.

III Proposed Algorithm

III-A Algorithm Overview

The proposed algorithm is shown in Fig. 2 and is intended to perform settingless unit protection of a medium-voltage line and identification of fault characteristics. Its inputs are:

• •

the parameters of the protected line (matrices $R$ and $L$ described in Section II), which are assumed to be accurately estimated using data-driven methods proposed in literature [13, 14];

• •

short observation windows containing samples of instantaneous three-phase voltages and currents from both ends of the line, which must be synchronized.

The algorithm consists of three main steps. The first step is to preprocess the input measurements to account for the missing data and obtain time derivatives of measured currents, which the algorithm requires since it relies on differential equations of a power line. At the second step, the algorithm evaluates how well two models of a power line, i.e. healthy and with an internal fault, fit the input measurements. To do this, a number of optimization problems are solved, which allows to reliably detect an internal fault and identify its characteristics as the ones that are the most consistent with the data. Based on the evaluation carried out at the second step, the last step chooses the model that best represents the current state of the protected line. This results in the identification of whether the line is faulty or not and, if yes, the most probable fault characteristics given the data. A detailed description of these steps is given below along with the hardware requirements for the implementation of the proposed algorithm.

III-B Preprocessing of Input Measurements

The proposed algorithm utilizes instantaneous measurements of voltages and currents and their derivatives for evaluation of healthy and faulty line models, which are based on differential equations. To obtain the correct evaluation, the windows of measurements have to be time aligned and their derivatives reliably estimated. This is ensured by the preprocessing step of the proposed algorithm.

At this step, the proposed algorithm initially checks if any samples of measurements from the other end of the line are lost during the data exchange via communication channels. If missing data is detected, the samples that correspond to the time stamps of this data are excluded from all the obtained windows of measurements. In this way the proposed algorithm ensures that all samples of measurements are time aligned. Note that by design of the second step of the algorithm, this procedure has no negative impact on the performance of the protection and fault identification functions provided that a small percentage of samples are lost. The methods for identification of lost data packages are covered in [15] and therefore, outside of the scope of this paper.

Next, the proposed algorithm determines the time derivatives of measured currents. This can be done by fitting successive sub-sets of a predefined number $l\geq 3$ of adjacent samples with a second-degree polynomial using linear least squares. The obtained polynomials are then differentiated at the time instances corresponding to the time stamps of the measurements’ samples. This method filters out high frequency noise while not impairing the quality of the derivative estimation.

Unfortunately, this method may incorrectly calculate the derivatives for several samples around the fault inception moment. To see this, consider a window of samples of the measured phase current with the fault occurring in the middle of the window, shown in Fig. 3a. At the fault inception moment, the true derivative of the current is undefined due to a sudden change of the line model. While the numerical estimate of this derivative always exists, it is incorrect whenever it is computed using the samples both before and after the inception moment. Therefore, there are $l-2$ samples for which the derivative estimate is erroneous. To enable reliable operation of the algorithm, they must be identified and removed from the windows of measurements and derivatives.

The identification of such samples is based on the numeric estimate of the second derivative of the current. Since it shows the rate of change of the first derivative, it has the maximum absolute value for the samples closest to the fault inception moment (see Fig. 3a). Hence, $l-2$ consecutive samples with the highest absolute values of the numerically estimated second derivative are removed from the observation windows. Note that these samples are removed even if the observation windows do not contain the fault inception moment as in Fig. 3b. However, it has virtually no negative impact on the protection and fault identification functions in the case with no fault provided that $l$ is small. It is worth noting that the wavelet transform can be alternatively used for detection of erroneous estimates of derivatives since it is well-known for finding sudden changes in the signals [16]. However, an empirical analysis indicated that the proposed method based on the second derivative has better performance than the wavelet transform for this task.

This step of the algorithm results in matrices $U_{1}$, $U_{2}$, $I_{1}$, $I_{2}$, $I^{\prime}_{1}$, $I^{\prime}_{2}\in\mathbb{R}^{3\times N}$ of samples of measured voltages and currents from both line terminals as well as the computed estimates of current derivatives. Here, $N$ denotes the total number of timestamps that remain after accounting for the missing data and errors in the derivative estimation. Each row of the matrices corresponds to a particular phase and each column corresponds to a particular time stamp. The obtained matrices are further utilized at the second step of the algorithm.

III-C Evaluation of Line Models

This is the main step of the algorithm as it is used to quantitatively evaluate how well the measurements fit two models: a healthy line and a line with an internal fault, described by equations (1)-(2) and (3)-(4), respectively. For the latter model this step also identifies the internal fault characteristics that best correspond to the input data. The proposed algorithm determines the fault location, type, values of fault resistances and inception moment. The main idea behind this step is to consider $(M+2)$ possible hypotheses that reflect the potential nature of the measurements in the observation window as shown in Fig. 4:

• •

all samples correspond to the normal operation of a protected line,

• •

all samples correspond to the time period after the inception of an internal fault,

• •

$M\geq 2$ possible cases of mixtures of samples before and after an internal fault in the observation window. In principle, $M$ should be equal to $N$ since an internal fault can occur at any moment in the window. However, to limit the number of considered cases, it is assumed that at each mixture case a fault has occurred within a particular interval of $R$ timestamps111Since $M$ must be a natural number, the values of $N$ and $R$ have to be selected accordingly., which thus gives $M=\frac{N}{R}$.

Clearly, one of these $M+2$ hypotheses must be true. Finding it requires quantitative assessment of the validity of all hypotheses. The validity of the $m$-th hypothesis is quantified through a mean squared error, denoted by $\Delta_{m}$, of the corresponding line model (healthy or faulty) or their mixture given the set of input measurements. Hence, the proposed algorithm constructs matrices $S$ and $W(R_{F},\alpha)\in\mathbb{R}^{6\times N}$ of mismatches of equations (1)-(2) and (3)-(4), respectively, defined as follows:

[TABLE]

[TABLE]

where $R_{F}\,:=\,\left[\,R_{a},\,R_{b},\,R_{c},\,R_{g}\,\right]^{T}$ is the vector of fault resistances shown in Fig 1. Note that while $S$ is a constant matrix because its entries are purely based on the measurement values, the elements of $W$ depend on the values of $R_{F}$ and $\alpha$, which are unknown. The procedure for computing $\Delta_{m}$ for all considered cases and estimating $R_{F}$ and $\alpha$ that fit the measurements the best is described below.

III-C1 Normal operation (case 1 in Fig. 4)

First, the proposed algorithm assesses the case that assumes that all samples in the observation window correspond to a healthy line. Hence, the value of $\Delta_{1}$ is computed as follows:

[TABLE]

The value of $\Delta_{1}$ serves as a measure of validity of this hypothesis given the input data.

III-C2 After-fault operation (case 2 in Fig. 4)

Next, the case that assumes that all samples in the observation window correspond to the time period after the internal fault inception is considered. For this case, the input data is fitted to the line model with a fault, described by (3)-(4). Since the fault characteristics are unknown, the proposed algorithm computes the mean squared error $\Delta_{2}$ of the faulty line model by determining the values of $R_{F}$ and $\alpha$ that minimize $\Delta_{2}$. This is done by solving the following optimization problem:

[TABLE]

where $R_{F}^{\max}$ is an upper bound placed on $R_{F}$ in order to consider only realistic values of fault resistances and improve numerical stability of the solution algorithm. Problem (9) can be rewritten in a standard form as

[TABLE]

where $x:=\left[\,R_{F},\,\alpha\,\right]$, $x^{\max}:=\left[\,R_{F}^{\max},\,1\,\right]$, $H\in\mathbb{R}^{5\times 5}$, $F\in\mathbb{R}^{5}$, and $d\in\mathbb{R}$. It follows from (7) and (9a) that $H$ is a positive-definite matrix, hence problem (10) is convex and can be solved reliably to global optimality by any commercial solver. The obtained value of $\Delta_{2}$ represents the smallest possible mean squared error of the faulty line model. Note that the optimal value of the variable vector, denoted by $x^{*}_{2}$, provides the fault characteristics that best fit the input data.

III-C3 Mixed operating conditions (cases 3 to M+2 in Fig. 4)

Next, $M$ cases are assessed in which it is assumed that the observation window includes measurements from before and after an internal fault. Let $1\leq k\leq M$ be the index of a given mixture case, thus $k=m-2$. Further, let $\mathcal{P}$ be the set of the indices of all timestamps in the observation window, i.e. $\mathcal{P}:=\left\{\,1,\,2,\,\dots,\,N\,\right\}$. As mentioned above, the proposed algorithm assumes that for each mixture case $k$ the fault has occurred within a particular time interval containing $R$ timestamps. Set $\mathcal{P}$ can be then divided into the following three disjoint sets:

• •

set $\mathcal{N}_{k}$ of the indices of timestamps corresponding to the assumed time interval before the fault;

• •

set $\mathcal{D}_{k}$ of the indices of timestamps corresponding to the assumed time interval of the fault inception;

• •

set $\mathcal{F}_{k}$ of the indices of timestamps corresponding to the assumed time interval after the fault.

This idea is illustrated in Fig. 4, in which set $\mathcal{N}_{k}$ corresponds to the green part of the observation window, set $\mathcal{D}_{k}$ to the red part, and set $\mathcal{F}_{k}$ to the blue part. These sets are defined for the $k$-th mixture case as follows:

[TABLE]

Naturally, it holds that $\mathcal{N}_{k}\cup\mathcal{D}_{k}\cup\mathcal{F}_{k}=\mathcal{P}$ for $\forall k$. Note that by design of (14) sets $D_{k}$ are selected in such a way that $\,\bigcup_{k=1}^{M}\mathcal{D}_{k}=\mathcal{P}$ and $\mathcal{D}_{i}\cap\mathcal{D}_{j}=\emptyset$ for all $i\neq j$. This ensures that for every possible fault inception moment in the observation window there exists one corresponding mixture case whose assumed fault inception interval contains this moment.

The mean squared error of each mixture case is computed by fitting the green part of the measurement windows to the healthy line model and the blue part to the faulty line model. The red part of each observation window is not utilized in fitting since there is an uncertainty which of the two line models this part corresponds to. Therefore, the following optimization problem is solved to obtain the value of $\Delta_{m}$:

[TABLE]

where $\eta=6(N-R)$. Similarly to (9), problem (16) can be rewritten in a standard form given by (10) and is also convex. Therefore, it can be solved fast and reliably to obtain $\Delta_{m}$, which is a measure of validity of the $k$-th mixture case (or $m$-th case among all $M+2$ cases). The optimal value of the variable vector $x^{*}_{m}:=[R_{F}^{*},\alpha^{*}]$ provides the fault characteristics that best fit the input data.

Note, the values of $N$ (length of the observation window) and $M$ (number of considered mixture cases) must be provided by the user. These values control the trade-off between the speed of the algorithm and the accuracy of estimating the fault characteristics. A smaller $N$ reduces the amount of data that has to be stored and transmitted over the communication channel, while a smaller $M$ decreases the number of mixture cases to be evaluated. On the other hand, higher values for $N$ and $M$ improve the accuracy of estimating $R_{F}$ and $\alpha$ at the expense of increasing the computational burden. The selection of optimal values for $N$ and $M$ is outside of the scope of this paper as it depends also on the capabilities of the protection devices and the communication channel.

III-D Comparison of Evaluated Cases

Once all $M+2$ cases have been processed, the algorithm proceeds to the selection of the case that best corresponds to the data in the measurement windows. This case is identified as the one that has the smallest mean squared error $\Delta_{m}$. If the model with no internal fault is determined to be the most likely one, the algorithm takes new measurement windows and repeats the procedure described above. Otherwise, it returns an estimate of the fault characteristics and sends a tripping signal to the corresponding circuit breaker. Hence, the proposed algorithm simultaneously performs protection, fault location and fault type detection for a power line. All these functions are crucial for stable and reliable operation of the grid. For instance, fast and accurate fault location performed by the proposed algorithm can reduce the grid restoration time. Fault type detection accomplished by identifying $R_{F}$ is important for power lines with single pole switching [17]. In addition, the values of $R_{F}$ can be utilized for collecting statistical records on faults.

It is important to note that for both cases $m=1$ and $m=3$, the faulty line model is not utilized, i.e. the mean squared errors are computed using only the healthy line model (see Fig. 4). The difference between these two cases is that all $N$ timestamps are used for computing $\Delta_{1}$, whereas $\Delta_{3}$ is obtained using only $N-R$ timestamps. Hence, limited computer precision and the presence of noisy measurements could result in $\Delta_{3}<\Delta_{1}$ even though in reality no internal fault has occurred within the observation window. To enable reliable operation of the proposed algorithm in such situations, an extra step is carried out whenever $\Delta_{3}$ and $\Delta_{1}$ have the smallest and the second smallest values among all $\Delta_{m}$. The idea is to compare the relative validity of case $m=1$ with respect to cases $m=3$ and $m=4$ in order to identify whether case $m=1$ is closer to $m=3$ (candidate for the true hypothesis) or case $m=4$ (false hypothesis). To do this, parameters $a_{1}:=\Delta_{1}/\Delta_{3}$ and $a_{2}:=\Delta_{4}/\Delta_{1}$ are computed and compared to each other. If $a_{1}<a_{2}$, the algorithm concludes that the line is healthy. Otherwise, the occurrence of an internal fault at the time interval corresponding to $\mathcal{D}_{1}$ is selected as the valid hypothesis.

A similar procedure is utilized for cases $m=2$ and $m=M+2$, since in both of these cases the mean squared error is computed using only the faulty line model.

III-E Hardware Requirements for Algorithm Implementation

The efficiency and reliability of the proposed algorithm depend on the quality of the input data. Hence, certain hardware requirements must be satisfied to ensure its correct operation. First, the input measurements need to be unbiased and have a low noise level. This requires the use of optical current and voltage transformers since conventional transformers might have significant deviations of the secondary signals from the primary ones [18]. Second, the proposed algorithm requires accurate estimation of time derivatives of the measured currents because it is based on analyzing differential equations of a power line model. To achieve this for an observation window of several milliseconds, a high sampling frequency in the range of dozens of kHz should be supported by a protection device. Third, since the algorithm relies on the exchange of measurements between the protection devices at two ends of the protected line, high-speed and high-bandwidth communication channels are required. Therefore, a fiber-optic channel or 5G communication is assumed to be utilized by the protection algorithm. It is worth noting that the aforementioned technologies are becoming increasingly available [2], which makes the proposed algorithm applicable in modern grids.

IV Simulation Results

The proposed algorithm was extensively tested for different types of generation, various levels of measurement noise, inaccuracies in assumed line parameters and different operating conditions: normal operation, internal and external faults.

IV-A Simulation Setup

The proposed algorithm was evaluated on the grid representation shown in Fig. 5, with the protected line modeled as described in Section II. The solidly grounded supplying systems $S1$ and $S2$ were modeled as voltage sources with magnitudes $E_{1}$ and $E_{2}$, angle difference $\theta$ and coupled $RL$ branches. The protection system of a power line consists of two relays $P1$ and $P2$ that are connected by a communication channel and run the proposed algorithm. Since $P1$ and $P2$ are identical, only the operation of $P1$ was analyzed in detail. The grid model and the proposed algorithm were implemented in Simulink and MATLAB, respectively. Table I shows the considered ranges of the parameters of the grid components. The parameters of the protected line include the line length $d$, the positive sequence resistance and inductance ${R}_{L}$ and ${L}_{L}$, respectively, and the ratio $K_{L}$ between the zero and positive sequence resistances and inductances. Similar parameters are defined for the models of the equivalent systems $S1$ and $S2$. Table I also provides the considered grid voltage level $U_{g}$ and its power frequency $f$. A fault in the grid was modeled as shown in Fig. 1. The considered fault types and corresponding ranges of fault characteristics are presented in Table II. Here, $T$ is the fault inception moment, $K_{3}$, $K_{2}$, $K_{2g}$ and $K_{1}$ denote three-phase, two-phase without ground, two-phase with ground and one-phase-to-ground fault types, respectively.

In the proposed algorithm, the sampling rate was set to 100 kHz and the derivatives of measured currents were obtained using subsets of $l=5$ adjacent samples. The upper bound on the elements of vector $x$ was set to $x^{\max}:=\left[\,\infty\,,\,\infty\,,\,\infty\,,\,\infty\,,\,1\,\right]^{T}$. This bound was utilized in optimization problems (9), (10) and (16), which were solved with Gurobi through its MATLAB interface [19]. The length of the observation window was set to 2 ms, which corresponds to 200 timestamps, and the number of considered mixture cases was set to $M=10$. The successive observation windows were assumed to be adjacent to each other, which minimizes the time between the fault inception and its detection.

IV-B Normal Grid Operation Conditions

Initially, the security of the proposed algorithm was evaluated on one hundred scenarios corresponding to normal operating conditions of the grid. For each scenario, the parameters of grid components were chosen randomly from the ranges provided in Table I and the grid model was simulated for 20 ms. Then, the measurements obtained from both line ends were fed to $P1$, where the proposed algorithm processed the measurement windows one by one. For all considered scenarios and all observation windows, the proposed algorithm correctly identified that the protected line was healthy for the entire simulation period.

IV-C External faults

Next, a number of tests were carried out to assess the security of the proposed algorithm in case of external faults. For this, one hundred sets of parameter values for the grid components were selected randomly from the ranges given in Table I. For each set, a large number of external fault instances were simulated. The diversity of these instances was ensured by simulating faults with all possible combinations of their characteristics within the limits defined in Table II with the following step sizes: $\Delta R_{a}=\Delta R_{b}=\Delta R_{c}=\Delta R_{g}=10\Omega$ and $\Delta T=0.5$ ms. Only external faults located at the buses incident to the protected line were examined as they are the most critical ones. These faults were assumed to be cleared by bus protection within 15 to 30 ms after their inception, which was chosen randomly for each simulation. The total simulation time was set to 70 ms. The measurements acquired after each simulation were analyzed by the proposed algorithm in the same way as described in the previous subsection. For all considered tests and observation windows, the proposed algorithm correctly identified that the protected line was healthy for the entire simulation time.

IV-D Internal faults

The dependability of the proposed algorithm for internal faults ($\alpha\!\in\![0,1]$) was assessed using one hundred sets of different values for the parameters of the grid components, which were selected randomly from Table I. For each set, internal faults were simulated with all possible combinations of their characteristics within the limits defined in Table II and the following step sizes: $\Delta R_{a}=\Delta R_{b}=\Delta R_{c}=\Delta R_{g}=10\Omega$, $\Delta\alpha=0.1$ p.u. and $\Delta T=0.5$ ms. The total simulation time of the grid model was set to 32 ms. The measurements acquired after each simulation were analyzed by the proposed algorithm in the same way as described above. For presentation clarity, the obtained results are classified into three groups:

• •

group 1: the observation windows that contained samples corresponding only to the pre-fault conditions;

• •

group 2: the observation windows that contained the inception moment of an internal fault;

• •

group 3: the observation windows that contained samples corresponding only to after-fault conditions.

IV-D1 Group 1

For all considered observation windows in this group, the proposed algorithm correctly identified that the protected line was healthy.

IV-D2 Group 2

The results for the second group are presented in Table III. As can be seen, the proposed algorithm correctly detected an internal fault for all considered fault instances. This table also presents the maximum and mean estimation errors of the fault location and fault resistances for each fault type. Note that the mean errors of fault resistances were calculated by averaging the maximum errors of estimating $R_{a}$, $R_{b}$, $R_{c}$ and $R_{g}$ over all tests.

IV-D3 Group 3

The results for the last group are also presented in Table III, with the same criteria for assessing the algorithm’s performance as for the second group. The algorithm correctly identified the state of the line as faulted for all considered faults. The maximum and mean estimation errors of fault location and resistances were significantly reduced for this group compared to the errors obtained for the second group. The reason for this is that the observation windows in group 3 contained more samples corresponding to after-fault conditions, which increased the estimation accuracy of the fault characteristics. Therefore, in practice the fault characteristics should be estimated using windows in this group.

The results obtained in Sections IV-B to IV-D demonstrate that the proposed algorithm can perform reliable line protection and reasonably accurate estimation of fault characteristics. The simulation results also indicate that the proposed algorithm correctly estimated the fault inception interval for all simulated instances of internal faults.

It is worth mentioning that successive observation windows can be made adjacent to each other only if the combined time of data collection and algorithm execution is below the window’s length of 2 ms. The data collection time comprises the delays of communication channel propagation and communication equipment. For the highest considered line length, it does not exceed 1 ms [3]. Hence, the proposed algorithm should be executed in under 1 ms. While the MATLAB implementation resulted in a slightly higher runtime, it could be significantly reduced by using a lower level programming language and dedicated hardware of a microprocessor relay.

IV-E Impact of converter-based generation sources

Converter-based generators may introduce a number of challenges for conventional protection algorithms [20]. Therefore, it is important to test the proposed algorithm on a grid model with this type of sources. Hence, source $S2$ in Fig. 5 was substituted by a model of a wind farm consisting of five identical Type 4 wind turbines. The parameters of grid components were set as given in Table IV, where $P$ is the nominal power, $H$ is the inertia constant, $R_{s}$ is the stator resistance and $X_{l}$ is the leakage reactance. The values of d- and q-axis reactances are also reported in the table. All parameters for which no units are provided are given in per unit.

To assess the security of the proposed algorithm for external faults and normal operation of the grid, a number of external fault instances were simulated and analyzed as described in Section IV-C. For all considered tests and observation windows, the proposed algorithm correctly identified that the protected line was healthy for the entire simulation time.

To assess the dependability of the proposed algorithm, a number of internal faults were simulated and analyzed as described in Section IV-D. For the observation windows that did not have samples corresponding to fault conditions (group 1), the proposed algorithm correctly identified the protected line as healthy. For the rest of observation windows (groups 2 and 3), the results are presented in Table V. The proposed algorithm correctly detected the internal fault for all considered fault instances. The maximum and mean estimation errors of fault location and resistances had similar values and followed the same trend as in Section IV-D. Note that the proposed algorithm also correctly estimated the fault inception interval for all simulated instances of internal faults.

IV-F Impact of measurement noise

Despite the high accuracy of optical instrumental transformers and advanced filtering methods utilized in modern protection relays, some noise can still be present in signals processed by protection algorithms. Therefore, the proposed algorithm was tested using scenarios described in Sections IV-B to IV-D and adding Gaussian noise to all measurements. Normal grid operation, external and internal faults were simulated for the same one hundred sets of the values of the grid component parameters, which were selected randomly from the ranges in Table I. Figure 6 presents the assessment of dependability and security of the proposed algorithm and estimation errors of fault location and resistances for different signal-to-noise ratios (SNRs) from 30 to 110 dB. Here, $N_{f}$ is the number of observation windows with samples corresponding to internal fault conditions, $N_{f}^{alg}$ is the number of times the proposed algorithm correctly identified an internal fault when analyzing these windows. Similarly, $N_{n}$ is the number of windows with samples corresponding only to normal grid conditions and/or external faults and $N_{n}^{alg}$ is the number of times the algorithm correctly identified the protected line state as healthy when analyzing these windows. The estimation errors in Fig. 6c and Fig. 6d are shown for observation windows of group 3 since these windows lead to a more reliable evaluation of fault characteristics as described above. Figure 6 shows the results only for $K_{3}$ faults; the other fault types led to similar results.

The results indicate that the dependability of the proposed algorithm was not affected by the considered increase in the noise level, whereas security was adequate only for SNR $\geq$ 60 dB. The fault location and estimation of fault resistances achieved reasonable accuracy approximately for the same range of SNRs. Therefore, the proposed algorithm was able to reliably perform all its functions for SNR $\geq$ 60 dB, which covers a realistic range of noise levels induced by modern optical instrumental transformers [21]. Note that the estimation error of the fault inception interval followed a similar trend to the one exhibited by the measure of security.

IV-G Impact of inaccuracies in line parameters

Finally, the robustness of the proposed algorithm to the inaccuracies in the input line parameters $L_{L}$ and $R_{L}$ was evaluated. This was done by simulating the same scenarios as described in Section IV-F but for fixed line parameters provided in Table IV and infinite SNR. The input values of $L_{L}$ and $R_{L}$ fed to the algorithm were changed by up to $\pm 20\%$ compared to the values utilized for the grid simulation. For all considered scenarios the algorithm was secure and dependable. Figure 7 presents the estimation errors of fault location and resistances for different deviations of line parameters utilized in the algorithm from the parameters employed for the simulations. The fault location and estimation of fault resistances achieved reasonable accuracy for $L_{L}$ and $R_{L}$ in the range of $\pm 10\%$, which is larger than the accuracy of available methods for line parameters’ estimation [13, 14]. Note that the proposed algorithm also correctly estimated the fault inception interval for all simulated internal fault instances.

V Conclusion

This paper presents a settingless optimization-based unit protection algorithm for medium-voltage power lines. The fact that no set-points need to be defined simplifies the algorithm’s application and reduces the risk of human errors. The algorithm combines protection, fault location and fault type identification functionalities, which allows for the unification of these closely related functions. The main idea of the proposed algorithm is to identify which model of a protected line, healthy or with an internal fault, best fits the input measurements. For that purpose, a number of convex optimization problems are solved, which allows to reliably detect an internal fault and estimate its characteristics that are most consistent with the data.

The proposed algorithm was extensively tested on a high number of scenarios. They included different levels of measurement noise, uncertainty of line parameters, and presence of converter-based generation in the grid. The results demonstrate that the algorithm was secure and dependable for all considered scenarios with realistic noise levels even without utilization of filtering methods. The algorithm also achieved accurate estimation of characteristics of internal faults for realistic levels of noise and errors in the assumed line parameters, even for high-impedance faults.

An important direction of future work is to make the proposed algorithm applicable for high voltage transmission lines, which require more detailed modeling of the lines.