Study on the algorithm of fault recording analysis combining its time-domain waveforms with phase-domain trajectories

The untimely handling of faults in a power system has a negative impact on its operation and even the national economy, and this requires coordination in the functions of protective relaying as well as supervisory & control devices, where digital fault recorders are used to record fault waveforms of electrical physical quantities. The fault recording of a simulated current is taken as the research object in this article, and it is transformed from the time-domain waveform into a phase-domain trajectory, which is used to analyze fault feature parameters and then reformulate the waveform. The original waveform of the current will be substituted by the reformulated one with fault features to realize functions in the power system. The algorithm of reformulating fault recording, the correlativity of the reformulated waveform and its original one, and errors produced in the research process are researched. The high correlation coefficient between the reformulated waveform and its original one shows that the algorithm studied in the article offers a simple and convenient option for fault recording analysis.

in an anticlockwise and clockwise direction and rotated in an appropriate angle, thus various components are separated from the fundamental component in real time and accurately 12,13 .
There are fault detection methods where effective artificial intelligence algorithms are highly applied to fault detection 14 .For example, a novel transformer model based on a deep convolutional neural network (CNN) is proposed, where 1-D deep CNNs are used to extract features in power systems 15 .In a fault detection algorithm based on a random forest (RF), the magnitude of maximum angular difference between positive and zero sequence component of the current at distributed generation buses is computed and then fed to a trained RF classifier for detecting fault conditions 16 .And a data-driven fault detection method with an ensemble classifier is used in modern distribution system 17 .
An extreme learning machine (ELM) algorithm for spontaneous fault detection and fault classification is proposed, which generates ten different types of normal and faulty data by two different transmission lines and normalizes the data to be used as a classifier for training 18 .The fault type detection and localization are achieved according to a histogram-based gradient boost (HGB) algorithm considering hosting capacity amendment in active distribution network 19 .Discrete wavelet transformation (DWT) and back-propagation neural network (BPNN) are also applied to fault diagnosis, where the high-frequency components in fault currents are extracted by DWT and the first-order coefficients detecting faults are studied to construct a BPNN-based decision algorithm 20 .
Besides, Kalman filter 21 , Bayes estimation 22 , petri nets 23 , cluster analysis 24 , support vector machine 25 , fuzzy set theory 26 , and rough set theory 27 etc. are also widely applied to the area of fault diagnosis in such way as multisource information fusion.However, a large amount of data in these studies are demanded and the computational processes are cumbersome.This increases the time cost and complexity of fault analysis and limits the depth and breadth of its application in power grids.
It is still not deep enough for all the present strategies of fault analysis to mine the characteristics of electrical physical quantities themselves, and the relevant analyses in terms of the phase domain are not involved.In the algorithm studied in this article artificial intelligent models are not trained, and the research priority is focused on the phase-domain trajectories of fault recording, which is innovated to acquire fault features in a power system.
The research results in the article are listed in "Results", including the expression of the reformulated waveform from its original one, their correlation coefficient, relative error, and the absolute error curves and the maximum absolute errors post both offline and online analysis.The algorithm of calculating fault features in fault recording and reformulating its waveforms, the accuracy requirement and errors are discussed in detail in "Methods".The significance of the research is listed in "Conclusion".

Results
A current of the fault recording pre and post a three-phase short circuit in an infinite power supply system is shown as follows: where ω is the angular frequency of the current.
A polar coordinate system is created by function "polarplot(theta, rho)" in Matlab, where "theta" and "rho" are polar angle θ and polar diameter ρ.Substituting θ for ωt and ρ for i in Eq. ( 1), where the system frequency is 50 Hz (ωt = 100π), there is: Since ωt is the phase angle of current i(t), and θ = ωt, Eq. ( 2) is referred to a phase-domain expression for Eq.(1), which is known as the time-domain expression for the current.Their figures are shown in Fig. 1a,b respectively.The process of forming the trajectory in Fig. 1b from the waveform of the short-circuit current in Fig. 1a is demonstrated by Supplementary Video 1 in Supplementary Information at the end of this article.
Comparing Fig. 1a,b, when the dc component decays, the apple-shaped trajectory becomes smaller and the balloon-shaped trajectory becomes larger, and both of them approach the blue dash-line circle in the shortcircuit steady state.These two trajectories are overlapped and stabilize at the position of this circle, as shown by the overlapped red solid-line circle and blue dashed-line circle in Fig. 1b.
The phase-domain trajectory of the dc component in the current is the green dashed-line spiral in Fig. 1b starting from the short-circuit instant and decays around the origin, and one circle is corresponded to one power frequency period.Its polar diameter decays to zero in the short-circuit steady state, and the decaying speed depends on its time constant.
It is challenging to distinguish whether the dc component in the short-circuit current still decays in the waveform at t = 0.2 s from Fig. 1a; however, from Fig. 1b at this instant the red solid-line trajectory is still not completely overlapped with the blue dashed-line trajectory, i.e., the current is still in the decaying process.
Due to the complexity of faults themselves and fault recording, it is tough to analyze them both quickly and accurately, and a methodological choice depends on the analysis purpose.When the fault recording is used to find failure causes, determine fault types and identify system operation modes, which are contributed to developing measures to improve the safe and stable operation of a power system, the requirement for real-time analysis is not high, while the requirement for the accuracy of fault feature parameters is high.For this case offline analysis is an approach.
(1) i(t) = sin ωt + 27.6 • , 0 ≤ t < 0.04s 4.36sin ω(t − 0.04) − 43.7 • + 3.476e −(t−0.04)/0.05, t ≥ 0.04s In particular, if there are normal steady states, transient states, short-circuit steady states, and states post the removal of failure in fault recording, periodic and non-periodic components are separated from it by the characteristics of circular and spiral phase-domain trajectories, and other components such as harmonics are considered further.
Fault recording is also included in microcomputer relay protection and automation, and it is mainly for the purpose of identifying fault types and determining fault distances in order to decide whether to send a tripping signal to a tripping circuit.In this case the real-time requirement for fault recording analysis is relatively high, while the requirement for the accuracy of fault feature parameters is not as high as that in offline analysis.Therefore, online analysis is suitable in this case.
Whether offline or online analysis is applied to fault recording, there inevitably are errors between a reformulated time-domain waveform and its original one.Therefore, only the reformulated waveform which meets certain accuracy requirement is the analysis result.
In the offline analysis of the short-circuit current in Fig. 1a, a correlation coefficient (c xy ) or relative error (ε) is used to reduce the absolute errors (e) in the reformulated waveform.When the accuracy requirement is that the correlation coefficient between the reformulated waveform and its original one is greater than or equal to c M = 1-10 −10 (0.999 999 999 999 9), the following reformulated waveform is obtained: Between the waveforms in Eqs.The curve of the absolute errors is shown in Fig. 2a.
In the online analysis only the waveform of the first 6 periods in Fig. 1a is used because the action time of regular quick protection is up to 0.12s 28 .When calculating the fault feature parameters, there is no accuracy requirement for the reformulated waveform since the rapidity is the primary consideration in protective relaying.
The online analysis of the sinusoidal current in the normal steady state and the adjustment process of its reformulated waveform are completely the same as that in offline analysis.The reformulated time-domain waveform of the short-circuit current in the transient state is: Between the waveforms in Eqs. ( 4

Methods
The relationship between the time-domain waveform of the current and its phase-domain trajectory The time-domain waveforms of the current in Eq. ( 1) and its phase-domain trajectories in Eq. ( 2) in the time interval of 0-0.12 s, 0-0.06 s, 0-0.08 s, and 1.96-2 s are shown in Fig. 3a-e.The fault recording of the current is digitally sampled from the instant when its phase angle is equal to 27.6°, and this instant is taken as the 0-th sample point, which is also described as a reference point, and there are: t = 0 s and θ = 0°.The phase angle of the reference point is equal to the initial phase angle of the current, which is φ 0 = 27.6°.
The coordinates of the points (a-k) in the waveforms shown in Fig. 3a,d are corresponded with those in the trajectories shown in Fig. 3b,c,e, which are listed in Table 1 respectively.
Point o: the reference point, initial phase angle of the current in the normal steady state, φ 0 = 27.6°;Point a: the 1st crest instant of the current in the normal steady state in the waveforms, the position of the 1st black circle in the trajectories; Point b: the short-circuit instant in the waveforms, the start of the maximum apple-shaped curve, i.e., the short-circuit reference point in the trajectories; Point c: the 1st crest instant of the short-circuit current in the waveforms, the position of the maximum polar diameter of the maximum apple-shaped curve in the trajectories; Point d: the 1st trough instant of the short-circuit current in the waveforms, the position of the maximum polar diameter of the minimum balloon-shaped curve in the trajectories; point e: the 2nd crest instant of the short-circuit current in the waveforms, the position of the maximum polar diameter of the 2nd largest apple-shaped curve in the trajectories; Point f: the 2nd trough instant of the short-circuit current in the waveforms, the position of the maximum polar diameter of the 2nd smallest balloon shaped curve in the trajectories; Point g: the initial instant of the short-circuit current in the short-circuit steady state in the waveforms, the start of the blue dashed-line circle in the trajectories; Point h: the 1st crest instant of the short-circuit current in the short-circuit steady state in the waveforms, the position of the diameter of the blue dashed-line circle in the trajectories, which is through the origin; Point i: the 1st trough instant of the short-circuit current in the short-circuit steady state in the waveforms, it is overlapped with point h in the trajectories; Point j: the initial instant of the decaying dc component in the waveforms, the start terminal of the green dashed-line spiral in the trajectories; Point k: one of the crest instants of the short-circuit current in the short-circuit steady state in the waveforms when the dc component decays to end, the position of the diameter of the red solid-line circle in the trajectories, which is through the origin.
In Fig. 3b the black solid-line circle containing point a and b is the trajectory of the sinusoidal current in the normal steady state.Point b is the short-circuit instant (the 0-th sample point or the reference point), and also the start of the apple-shaped curve.Since the decaying dc component is greater than zero, the 1st crest located above the time axis is corresponded to the position of the maximum polar diameter of the maximum appleshaped curve, and the 1st trough located below the time axis-the minimum balloon-shaped curve.When the dc component decays, the "apple" becomes gradually smaller and the "balloon"-larger, and both of them approach the circle of the sinusoidal current in the short-circuit steady state and stabilize at the position of this circle, as shown by the overlapped red solid-line circle and blue dashed-line circle in Fig. 3e.The 1st largest apple and the 1st smallest balloon are shown in Fig. 3b, and the 2nd larger apple and the 2nd smaller balloon-in Fig. 3c.
The phase-domain trajectory of the sinusoidal ac component in the current in the short-circuit steady state is a circle too, which starts from the short-circuit instant (point g) and is also the trajectory of the sinusoidal current in the short-circuit steady state, i.e., the blue dashed-line circles shown in Fig. 3b,c,e.The phase-domain trajectory of the decaying dc component in the current is the spiral, which is the green dashed-line spiral in Fig. 3b,c.It starts from the short-circuit instant (point j) and decays around the origin, and rotates one circle per power frequency period.Its polar diameter decays to zero in the short-circuit steady state, and its decaying speed depends on the time constant of the exponential function.
After the phase-domain trajectory of a short-circuit current (Eq.( 1)) is formed by an appropriate way, which is the expression in Eq. ( 2), the following short-circuit fault feature parameters are accessible from analyzing the trajectory with its time-domain waveform: the amplitude and initial phase angle of the sinusoidal current in the normal steady state, the short-circuit initial phase angle (the initial phase angle of the sinusoidal ac component at the short-circuit instant), the value of the current at the short-circuit instant, the initial value and time constant of the decaying dc component, the amplitude of the sinusoidal current in the short-circuit steady state.The original waveform of the current is thus reformulated by these features.

Calculating the time-domain waveform from the circular phase-domain trajectory of the sinusoidal current in the normal steady state
The phase-domain trajectory of a sinusoid is circular, and the diameter of the circle is equal to its amplitude.The relationship between the position (θ m ) of the circle (the polar angle of the diameter through the origin) and the initial phase angle (φ 0 ) of the sinusoid is: It is uncomplicated to obtain the amplitude and initial phase angle of the sinusoid from the diameter of the circular trajectory and its position.An especial attention should be paid to that when the polar diameter of this diameter is positive, it is more convenient to convert its polar angle to the value between 0° and 360°.If the polar diameter is negative, it should be firstly reflected through the origin to the position where its polar diameter is positive, and then its polar angle is converted to the convenient value between 0° and 360°.
For instance, in the trajectory (Fig. 3b) of the sinusoidal current in the normal steady state, point a is the position of the diameter (through the origin) of the trajectory circle, its polar diameter and polar angle are 1 and 62.1° respectively, which means that the position of the circle is θ m = 62.1°.Therefore, the amplitude of the sinusoid is 1 and its calculated initial phase angle is obtained by Eq. ( 5), which is: The polar angle and polar diameter of point o (the reference point) and b (the short-circuit reference point) are found from the trajectory circle in Fig. 3b, which is also list in Table 1.Supposing the duration of the sinusoidal current in the normal steady state is 0.04 s, and its frequency is 50 Hz, the expression of the raw reformulated waveform of this current is: Comparing Eq. ( 6) with the original function in Eq. ( 1), the calculated initial phase angle (φ 0C = 27.9°)from the trajectory circle is not equal to the actual angle (φ 0 = 27.6°), and this results in errors between the raw reformulated waveform and the original one.The absolute error curve is shown in Fig. 4a, where the maximum absolute error (e max = 0.005 236) is marked. (5) Table 1.The comparison of the coordinates of the points (a-k) in the waveforms (shown in Fig. 3a,d) with those in trajectories (shown in Fig. 3b,c,e).In the table "Theta" expressed by "angle 1 / angle 2" is the polar angle, where "angle 1" is the phase angle of a point based on the reference point, "angle 2" is the value which "angle 1" is converted to between 0° and 360°."Position θ m " is the polar angle between 0° and 360°, and when the polar diameter of a point is negative, θ m is the polar angle of the positive polar diameter to which the point is reflected through the origin.The correlation analysis between the reformulated waveform and its original one in the normal steady state After the features such as the amplitude and initial phase angle of the sinusoid from the circular trajectory are obtained, the similarity between the reformulated waveform and the original one is required to be measured by certain indexes.Obviously, it is not appropriate to use the absolute error curve in Fig. 4a.The similarity or interdependency between the two signals of x(t) and y(t) is measured by a correlation coefficient (c xy ) or a relative error (ε).When y(t) is used to approximate to x(t) in the time interval of t 1 -t 2 , the correlation coefficient c xy is calculated by the following equation 29 : And the relative error ε is: The actual initial angle (φ 0 ) of the original sinusoidal waveform is unknown pre its determination.In order to obtain the calculated initial angle (φ 0C ) within the allowable error range, a sample window (Δφ) (its definition and calculation are found below) is taken as an adjustment range of the phase angle, and φ 0 is located in the interval which is centered on φ 0C : Taking this interval as an initial one and φ 0C as an initial value, the initial angle which meets the requirement for an accuracy will be found in the interval by a dichotomy method.The requirement is expressed by the maximum correlation coefficient (c M ) between the reformulated waveform and its original one, and when the coefficient calculated by Eq. ( 7) meets the condition of c xy ≥ c M , the adjustment process is ended.The middle value of the present interval is exactly the calculation result of the initial angle.
The calculated initial phase angle from the trajectory circle in Fig. 3b is φ 0C = 27.9°, and it is the initial value.Taking 0.9° (the value of the sample window) as the adjustment range, the actual initial phase angle (27.6°) is located in the interval of [27.0°, 28.8°].The adjustment algorithm is shown in Fig. 5a, the accuracy requirement is that the correlation coefficient c xy is not less than c M = 1-10 −10 (0.999 999 999 9), and the data in the adjustment process are listed in Table 2.The program list of adjusting process of the sinusoidal current in the normal steady state, which is corresponding to the block chart in Fig. 5a, is referred to Supplementary Information 1 in Supplementary Information at the end of this article.
When n = 9, c xy has exceeded c M and the adjustment result is φ 0C = 27.599414 06°, which is smaller than the actual one (27.6°)by 0.000 585 938°.The reformulated sinusoidal function meeting the accuracy requirement is changed from Eq. ( 6) to the following one: The curve of the absolute errors between the reformulated sinusoidal waveform and its original one is shown in Fig. 4b, where the maximum absolute error e max = − 1.023 × 10 −5 is marked.

Calculation of the sinusoidal waveform from the circular trajectory of the sinusoidal current in the short-circuit steady state and the correlation analysis between the reformulated waveform and its original one
Following the same method as analyzing the sinusoidal current in the normal steady state, the amplitude and the calculated short-circuit initial phase angle of the sinusoidal waveform (Fig. 3d) in the time interval 1.96-2 s in the short-circuit steady state are obtained from the polar diameter of 4.36 and the polar angle of 354 14° (or 134°) of point k in the trajectory in Fig. 3e (or Table 1), and they are 4.36 and φ 0kC = − 44° respectively.Taking point o as the reference point, the expression of the raw reformulated function of the waveform in this time interval is: There is an error between the calculated initial phase angle (φ 0kC = − 44°) and the actual one (φ 0k = − 43.7°).The curve of the absolute errors between the initial reformulated waveform and the original one caused by the error is shown in Fig. 4c, where the maximum absolute error e max = − 0.022 83 is marked.
Taking − 44° as the initial value and 0.9° as the adjustment range, the actual angle of − 43.7° is located in the initial interval of [− 44.9°, − 43.1°].Following the adjustment algorithm in Fig. 5a, the accuracy requirement is taken as the same one as in the normal steady state (the correlation coefficient is not less than c M = 1-10 −10 ), the data in the adjustment process are listed in Table 3.
Table 2.The data in the adjustment process of the initial phase angle of the sinusoidal waveform in the normal steady state, where the relative errors (ε) are calculated by Eq. ( 8), the maximum absolute errors (e max ) are the maximum among the absolute errors between the reformulated waveform and its original one.The program list of adjusting process of the sinusoidal ac component in the short-circuit current, which iscorresponding to the block chart in Fig. 5, is referred to Supplementary Information 1 in Supplementary Information at the end of this article.
When n = 9, the correlation coefficient between the reformulated function and its original one has exceeded c M .The adjustment result is φ 0kC = − 43.699 414 06° and is larger than actual angle (− 43.7°) by 0.000 585 938°.The reformulated sinusoidal function meeting the accuracy requirement is changed from Eq. ( 10) to the following one: The curve of the absolute errors between the reformulated sinusoidal waveform and the original one is shown in Fig. 4d, where the maximum absolute error e max = 4.459 × 10 −5 is marked.

Calculation of the decaying dc component from the time-domain waveform of the short-circuit current and the correlation analysis between the reformulated waveform and its original one
Equation ( 11) is the sinusoidal current in the short-circuit steady state, from it the expression of the sinusoidal ac component in the time interval of 0.04-0.44s in the transient process is directly written to be: Subtracting the sinusoidal ac component (Eq.( 12)) from the short-circuit current (Eq.( 1)), the raw expression of the decaying dc component is: The waveform is plotted in Fig. 6a,b.The initial phase angle is − 43.699 414 06° from the reformulated sinusoidal ac component (see Eq. ( 12)), and the actual one is − 43.7° (see Eq. ( 1)).The difference causes the errors in the expression in Eq. ( 13), and the waveform in Fig. 6a,b is not accurate dc component.The errors are: Such small errors cannot be obviously observed in Fig. 6a,b.Expressing the initial value and time constant of the dc component as M and T a , its expression of the waveform is: Taking any two points in the waveform: point m(t m , i m ) and point n(t n , i n ) (see Fig. 6a), and substituting their coordinates into Eq.( 14), there are: Solving this set, the time constant (T a ) and initial value (M) are obtained as follows: 14)  i ap (t) = Me −(t−0.04)/T a , 0.04s ≤ t < 0.44s www.nature.com/scientificreports/Substituting T a and M back into Eq.( 14), the expression of the raw reformulated waveform of the dc component is then obtained, and there inevitable are errors in it.Taking the waveform of a segment from the time intervals of the waveform as an adjustment object, calculating its time constant (T a ) and initial value (M) by Eq. ( 15) using the coordinates of the two boundary points, and an expression of the dc component is formed.When the accuracy requirement between the reformulated waveform and the original one (see Eq. ( 13)) is met, the adjustment process is completed.
Taking the initial instant where the dc component starts to decay as a start terminal (point m(t m , i m )) and the instant where it decays to just less than or equal to 1% of the initial value (M)-an end one (point n(t n , i n )), the waveform of this segment is described as an initial adjustment object, and the duration from the start terminal to the end one-an initial time interval [t m , t n ], as shown in Fig. 6b, there are: The value of correlation coefficient c M between the reformulated waveform and the original one is still taken as an accuracy requirement.The time constant (T a ) and decaying initial value (M) of the initial adjustment object are calculated by Eq. ( 15) using the coordinates of the two boundary points (m(t m , i m ) and n(t n , i n )), and they are substituted into Eq.( 14) to calculate correlation coefficient c xy between the reformulated waveform and the original one.When the condition of c xy ≥ c M is satisfied, the reformulated dc component meets the accuracy requirement, and the adjustment process of the waveform is completed.When the condition of c xy ≥ c M is not satisfied, the waveform of the initial adjustment object is dichotomized into two segments, the first one is taken as a new adjustment object.After obtaining its T a and M by Eq. ( 15), another waveform of the dc component is reformulated to calculate the correlation coefficient to determine whether its accuracy meets the requirement.
The adjustment process is detailed in the block chart shown in Fig. 5b, where after T a and M of each segment are calculated, whether the accuracy of the reformulated waveform meets the requirement, the initial adjustment object is reformulated using these T a and M. If the accuracy of the reformulated waveform meets the requirement, there is no need to continue the subsequent adjustment process, this reformulated waveform is the final result.
In the adjustment process of the dc component, if the duration of any segment is equal to one power frequency period, regardless of the accuracy of the reformulated waveform, there is no need to continue the subsequent adjustment in this time interval.After the adjustment process is completed, if all time intervals are equal to or less than one power frequency period and all the reformulated waveforms do not meet the accuracy requirement, the values of T a and M of the segment with the largest correlation coefficient are taken to reformulate the dc component, and it is also the final result of the adjustment.
The program list of adjusting process of the decaying dc component in the short-circuit current, which is corresponding to the block chart in Figure 5b, is referred to Supplementary Information 2 in Supplementary Information at the end of this article.
In the waveform of the dc component in Fig. 6b, the short-circuit instant is taken as point m, which is: The value of 3.475 967 764 is actually the initial value of the dc component, which is M = 3.475 967 764.When the dc component decays to less than 1% of the initial value (M), the decaying process is nearly completed, therefore, the point where the dc component is less than or equal to 1%M is taken as point n, which is: The time constant and initial value are calculated by Eq. ( 15) using the two coordinates (t m , i m ) and (t n , i n ) to be: Substituting them into Eq.( 14), the expression of the reformulated dc component is: Taking the time interval of 0.04-0.27035 s as the initial time interval of the adjustment, the data in the adjustment process are listed in Table 4, where the accuracy requirement is: the correlation coefficient (c xy ) is not less than c M = 1-10 −10 (0.999 999 999 9).The time constant and decay initial value of the dc component meeting the accuracy requirement are:   16) and ( 17) are the expressions of the dc component pre-and post-adjustment, and the curves of the absolute errors between either of them and the expression in Eq. ( 13) are shown in Fig. 4e,f, where the maximum absolute errors are also marked.

The online analysis of fault recording
The online analysis of the sinusoidal current in the normal steady state and the adjustment process of the reformulated waveform are completely the same as that in offline analysis.Due to the lack of the sinusoidal current in the short-circuit steady state in quick protection, it is trickier to separate the sinusoidal ac component from the short circuit current before analyzing the dc component than the above-mentioned offline analysis.

Estimating the decaying dc component and sinusoidal ac component in the short circuit current in 40 ms
The regular expression of the short-circuit current in Fig. 3a is: The decrease of the "apple trajectory" and increase of the "balloon trajectory" as well as their polar angles are utilized to calculate the decaying dc component in the current, as shown in Fig. 3c.Point c and e are the positions of the maximum polar diameters of the largest and second largest apples, and point d and f are the positions of the maximum polar diameters of the smallest and second smallest balloons.The corresponded to them same points are also shown in Fig. 3a, and the polar diameters and polar angles of these points are listed in Table 1.Substituting the coordinates of these points into Eq.( 18), there are: Since the time interval between point c and e and that between point d and f are all one power frequency period, the sinusoidal ac components in the first two equations are equal, and the ones in the second two equations are also equal.The above equations are simplified to be: Solving Eq. ( 19), there are: Substituting the two values of T a and M into Eq.( 18), the expression of the decaying dc component is: Subtracting this dc component from the short-circuit current (see Eq. ( 18)), the waveform of the sinusoidal ac component is obtained and shown in Fig. 7a, where there are only 6 periods post the short circuit.Taking (17)  i ′ ap (t) = 3.475967764e −(t−0.04)/0.05000030722627764, 0.04s ≤ t < 0.44s A k sin[100π(0.0473− 0.04) + ϕ 0k ] + Me −(0.0473−0.04)/T a = 7.36 A k sin[100π(0.0673− 0.04) + ϕ 0k ] + Me −(0.0673−0.04)/T a = 6.37A k sin[100π(0.0573− 0.04) + ϕ 0k ] + Me −(0.0573−0.04)/T a = −1.897A k sin[100π(0.0773− 0.04) + ϕ 0k ] + Me −(0.0773−0.04)/T a = −2.708 the instant at t = 0 as a reference point, the phase-domain trajectory of the first half period (0.04-0.05 s) of the waveform is plotted in Fig. 7b.Following the same way as that in the offline analysis, from the polar angle of 854.1° (or 134.1°) and the polar diameter of 4.354 of the trajectory circle shown in Fig. 7b, the expression of the sinusoidal ac component is immediately obtained: Since the dc component is calculated by two periods (0.04-0.08 s), at least 40 ms are needed to obtain the decaying dc component and the sinusoidal ac component of the current in this kind of online analysis.

Estimating the sinusoidal ac component in the short circuit current in 20 ms
The decaying dc component in the short-circuit current impacts on only the current values but not the phase angles of the sinusoidal ac component, which is known from Fig. 1b and Fig. 3b,c.To improve the rapidity of protection, the sinusoidal ac component can be estimated in one period by the phase-domain trajectory of the short-circuit current.
As shown in Fig. 3b, the phase angles of point c and d approach the phase angle of point h.The average of the positions of point c (131.4°) and d (131°) is taken as the position of the trajectory circle corresponded to waveform of the sinusoidal ac component: The initial phase angle of the sinusoidal ac component is calculated by Eq. ( 5), and it is: The amplitude of the sinusoidal ac component is estimated by the average of the absolute values of the polar diameters of point c (7.36) and d (− 1.897), which is: Then the sinusoidal ac component is estimated as: Obviously, the expression in Eq. ( 22) is more inaccurate than that in Eqs. ( 1) and (21).

Errors in reformulated waveforms
From the relationship between Eqs. (1) and (2), when the phase-domain trajectory of a sampled waveform in fault recording is plotted, the polar diameter of each sample point is the instantaneous value of the waveform, and the polar angle of each sample point-its phase angle.In order to guarantee the correspondence between a waveform and its trajectory, it is inevitable to determine the 0-th sample point (SP) (the 0-th SP or SP0), where there is t = 0 and θ = 0°.This sample point is referred to as a reference one, and its phase angle is equal to the initial phase angle of the waveform.
The reference point can be any one of sample points, and it is usually the first one following a zero-crossing instant from negative to positive (or the reverse).Once the phase angle of the reference point, which is also SP0, is fixed, the phase angles of following SP1, SP2, … are all calculated from the initial phase angle (φ 0 ) and the sample period (T s ).Obviously, only one reference point is needed in a waveform, otherwise there appeared confusion during its processing, and correct results will not be obtained.
The phase angle corresponded to the time interval between any two sample points is defined as a phase window.For example, the phase angle corresponded to one power frequency period (360° or 2π radians) is described as a power frequency window, the phase angle corresponded to a half power frequency period (180° or π radians)-a half power frequency window, and the phase angle corresponded to one sample period in one power frequency period-a sample window.
Expressing the system frequency, sample frequency, and sample window as f 0 , f s , and Δφ, the number of the sample points in a power frequency period is f s / f 0 , and the value of Δφ is calculated by the following equation: A power frequency period is fixed, which is 360° or 2π radians, however, the size of a sample window varies with the changes of the system frequency and the sample frequency (see Eq. ( 23)).When the sample frequency remains unchanged, the size of the sample window is needed to be adjusted with the change of the system frequency.
The zero-crossing instant from negative to positive of a current in Fig. 8 is t (0) ʹ.Current i −1 of the previous point (SP(t −1 , i −1 )) is negative, current i 0 of subsequent point (SP(t 0 , i 0 )) is positive, and t 0 − t −1 is the sample period (T s ).SP(t 0 , i 0 ) is the first positive sample point following the zero-crossing instant, and it is taken as a reference point (the 0-th SP or SP0).
Connecting SP(t −1 , i −1 ) with SP(t 0 , i 0 ) in Fig. 8.The figure demonstrates how to determine the phase angle of a reference point using the linear interpolation method, and to form a green straight line, its intersection (t (0) ) with the time axis is approximated to the actual zero-crossing instant (t (0) ′).When sample frequency f s is much www.nature.com/scientificreports/greater than system frequency f 0 , instant t (0) is approximated to instant t (0) ′.The value of t (0) is calculated by a linear interpolation method, which is: The calculated value (φ 0C ) of the actual initial phase angle (φ 0 ) is calculated by: The calculated phase angle of the n-th sample point (SP(t n , i n )) φ nC is calculated by the equation: Since the calculated zero-crossing instant (t (0) ) by Eq. ( 24) is not the actual one (t (0) ʹ), the error results in errors in φ 0C and φ nC calculated by Eqs.(25) and (26).These errors are smaller when sample frequency f s is much greater than system frequency f 0 .
As shown in Fig. 8, the error (e φ ) in the calculated initial phase angle (φ 0C ) is: Its percentage relative to the sample window is calculated by the following equation: When sample frequency f s = 20 000 Hz, actual zero-crossing instant t (0) ʹ = 0.004 674 559 726 073 8 s, and calculated zero-crossing instant t (0) = 0.004 674 597 260 738 6 s, substituting them into Eq.( 27), the error in the phase angles of the sample points is: The larger the correlation coefficient between the reformulated waveform and its original one, the smaller the error in the reformulated waveform.The maximum correlation coefficient is 1, which means the original waveform can be completely restored from the reformulated one.However, due to the interference in the digitally sample waveforms in practice, the exact restoration of the original waveform is not contributed to the analysis of fault recording.The accuracy requirement is decided by the application scenarios of the results of fault recording analysis ("Supplementary information").

Conclusion
The research results show that the accuracy of reformulated fault recording is related to the features of the fault recording itself, the analysis purposes and the accuracy requirement.The algorithm studied in this article is not complex, where there are not time-consuming calculations, and the accuracy from the simulation data is relatively high.It is more reliable and understandable than other algorithms and easy to implement and spread in practice.
At present, the process analyzing and studying fault recording in a phase domain is still in an initial stage.In order to establish a theoretical basis, the simulated waveform in the article is used to achieve the algorithm.There is still a large gap between a simulated waveform and a practical one, and it is necessary to study them further in our subsequent researches.Let us make great efforts to work and look forward to it together.
) and (1) the correlation coefficient and relative error are: c xy = 0.999980063830419, ε = 3.98719417112892 × 10 −5 .The maximum absolute error in the first two periods (40 ms) is e max = − 0.03722.The curve of the absolute errors in the time interval of 0.04-0.16s is shown in Fig. 2b.

Figure 2 .
Figure 2. The absolute error curves of the reformulated waveforms of the fault recording.The figure in (a) is the errors in the reformulated waveform in offline analysis, and the figure in (b)-in online analysis.The maximum absolute errors are all marked in the two figures.

Figure 3 .
Figure 3.The time-domain waveforms of the current expressed by Eq. (1) and its phase-domain trajectories expressed by Eq. (2).The waveform in the time interval of 0-0.12 s is shown in (a), its trajectories in the time interval of 0-0.06 s and 0-0.08 s are shown in (b) and (c), the waveform and its trajectory in the time interval of 1.96-2 s are shown in (d) and (e).The black solid lines indicate the current in the normal steady state, and the red solid lines-in the transient state.The blue dash lines indicate the sinusoidal ac component, and the green dash lines-the decaying dc component.

Figure 4 .
Figure 4.The figures in (a) and (b) are the absolute error curves of the reformulated waveforms of the sinusoidal current in the normal steady state pre and post the adjustment.The figures in (c) and (d) are the absolute error curves of the reformulated waveforms of the sinusoidal current in the short-circuit steady state pre and post the adjustment.The figures in (e) and (f) are the absolute error curves of the reformulated waveforms of the decaying dc component in the short-circuit current.

( 8 )Figure 5 .
Figure 5.The block chart for the adjustment algorithm of the initial phase angle of a sinusoidal current is (a), and (b) is the block chart of the adjustment algorithm of the decaying dc component in a short-circuit current.

Figure 6 .
Figure 6.The waveform of the raw expression of the decaying dc component in the short-circuit current and the calculation of its time constant (T a ) and initial value (M).In (a) point m and n are any two points on the waveform.In point m is the one where i m = M, and point n is the one where i n ≤ 1%M.

( 19 )Figure 7 .
Figure 7.The sinusoidal ac component obtained after the decaying dc component is subtracted from the shortcircuit current, where (a) is the time-domain waveform of 6 periods post the short circuit and (b) is the phasedomain trajectory of the first half period of the waveform in (a).

Figure 8 .
Figure8.The figure demonstrates how to determine the phase angle of a reference point using the linear interpolation method, and how to estimate the calculated error in the phase angle of a sample point.The red curves are the sampled waveform of a current, and the green straight lines are formed by connecting the two sample points of the previous and subsequent the zero-crossing instant.

Table 4 .
The data in the adjustment process of the decaying dc component in the short-circuit current.