A Novel Switching Algorithm to the Preferred Clock Skew Estimator Applicable for the PTP Case in the Fractional Gaussian Noise Environment

,


Introduction
Te PTP, also known as IEEE1588v2, is intended to estimate the time (ofset or phase) and frequency synchronization between each clock in the network [1]. Te clock skew and the ofset are estimated based on messages with time stamps that the Master and Slave clocks exchange. Te asymmetry between the Forward and the Reverse paths is challenging in this estimating task. Te frst asymmetry is in the fxed delay of the Forward and Reverse paths. Since the messages can pass through a diferent number of switches or routers in each path, the fxed delay may be diferent. Terefore, assuming symmetry in the fxed delay between the paths may cause an estimation error [2][3][4]. Te random delay, commonly known as PDV, contains the second asymmetry. Te asymmetry, in this case, is brought either by actual background trafc or, occasionally, by cyberattacks on the Forward or Reverse path. It makes sense to think that high background trafc will boost PDV. Higher PDV may, however, reduce the frequency or time synchronization accuracy that the PTP obtains [5][6][7]. Te last asymmetry is in the Hurst exponent parameters of the Forward and Reverse paths that may have diferent values. When we have a higher Hurst exponent parameter in one of the paths (Forward or Reverse), it also raises the PDV. Recently, we proposed three clock skew estimators (a TWD clock skew estimator [5], an OWD clock skew estimator for the Forward path [6], and an OWD clock skew estimator for the Reverse path [6]) which have better clock skew performances from the MSE point of view compared to [8][9][10]. But the question of which of these three new clock skew estimators ( [5,6]) we should use under which PDV condition was still unanswered. In [7], we presented a switching algorithm that switches to the preferred clock skew estimator among our three novel clock skew estimators, the TWD clock skew estimator [5], the OWD clock skew estimator for the Forward path [6], and the OWD clock skew estimator for the Reverse path [6]. Te algorithm in [7] switched to the clock skew estimator with the best performance from the MSE perspective. However, the algorithm in [7] only applies to the Gaussian case. In this paper, we propose a novel switching algorithm (switching between the TWD clock skew algorithm from [5], the OWD clock skew algorithm for the Forward path from [6], and between the OWD clock skew algorithm for the Reverse path from [6]) applicable to a practical system where the PDV is characterized by an fGn process where the Hurst exponent parameter is in the range of 0.5 ≤ H < 1 and where after a small set of PTP measurements, the switching algorithm should be able to switch efectively to the preferred clock skew estimator with the best performance in the MSE perspective. Note that the Gaussian case in [7] is a particular case of the fGn case where H � 0.5. For 0.5 < H < 1, the corresponding fGn is of long-range dependency (LRD) [11]. According to [12], long-term persistence in natural processes was frst discovered by the authors in [13]. Te basics of long-range dependent time series can be found in [14], while in [15], we have the relationship of fGn and trafc modeling as well as the relationship of generalized fractional Gaussian noise (gfGn) and trafc modeling. Obviously, the Hurst exponent parameter has to be given or estimated when applying the fGn model with 0.5 < H < 1. In the PTP case where the PDV process is characterized as an fGn model with 0.5 ≤ H < 1, the Hurst exponent has to be estimated since it is not known by the system designer. It should be pointed out that (1) we do not have on hand the PDV samples. We can only measure indirectly each PDV sample together with a time-dependent ofset (please refer to the discussion section). (2) Te Hurst exponent parameter estimation task has to be carried out in real time. Tus, the estimation of the Hurst exponent parameter should be obtained as fast as possible, using only a small set of measurements such as 400 − 1000 measurements. After this period of time, the estimation of the Hurst exponent parameter should be efectively used in the newly proposed switching algorithm.
In [12], three common stochastic tools, the climacogram, i.e., variance of the time averaged process over the averaging time scale, the autocovariance function, and the power spectrum were compared to each other to assess each one's advantages and disadvantages in stochastic modeling and statistical inference. It should be pointed out that the climacogram estimator was introduced as a concept by the authors in [16] and mathematically defned (including estimation bias) by the authors in [17]. According to [12], the classical climacogram estimator has a small error as well as an expected value always positive, well-behaved, and close to its mode (most probable value), all of which are important advantages in stochastic model building. In contrast, the power spectrum and the autocovariance do not have some of these properties. Tus, according to the authors in [12], we should start with the climacogram rather than the power spectrum or the autocovariance for the Hurst exponent parameter estimation task. But the classical climacogram estimator [12] is not applicable in the PTP case where each PDV sample can only be measured indirectly together with a time-dependent ofset (please refer to the discussion section). As a matter of fact, every estimation tool for the Hurst exponent parameter that uses a set of "k" PTP measurements out of "n" PTP measurements and subtracts the average of "n" PTP measurements does not work properly due to the time-varying ofset that is added to each PDV sample in the PTP measurements (please refer to the discussion section regarding the time-varying ofset).
In this paper, we derive an algorithm for estimating the Hurst exponent parameters for the Forward and Reverse paths. After we have the estimated Hurst exponent parameters, the switching algorithm can choose the preferred clock skew estimator with the best clock skew performance from the MSE perspective between the three clock skew estimators: (1) the TWD clock skew estimator [5], (2) the OWD clock skew estimator for the Forward path [6], and (3) the OWD clock skew estimator for the Reverse path [6]. In [7], we presented two algorithms to detect the unexpected load or cyber-attack in one of the paths (Forward or Reverse) by tracking an increase in the PDV. Tose detection algorithms [7] apply only to the Gaussian case. However, in the fGn case with 0.5 < H < 1, those detection algorithms cannot alarm properly. Estimating the Hurst exponent parameters for the Forward and Reverse paths must be carried out before starting to estimate the PDV variances for each path. After estimating the Hurst exponent parameters for the Forward and Reverse paths, we can estimate the PDV variances for each path and derive the overload, unexpected load, and cyber-attack algorithms for the fGn case. It should be pointed out that the Gaussian case is a particular case of the fGn process where the Hurst exponent parameter is equal to 0.5. Tus, our new proposed algorithms answer on a wider range for Hurst exponent parameters.
In our previous works ( [5][6][7]), we assumed that the clock skew is a constant parameter. Tis assumption is consistent with the literature [2][3][4][8][9][10]. However, in practical systems, malfunctions in the Master or Slave clock, such as a high temperature, environmental conditions, or physical phenomena, may cause drift in the clock skew over time [18,19]. In those cases, the system of Master and Slave may get out of synchronization. Terefore, it would be useful for the system designer to have an alarm indicating any drift in the clock skew associated with the system. Te summary of this paper's contribution is as follows: (1) estimating the Hurst exponent parameters for the Forward and Reverse paths based on time stamps related to the PTP case for the frst time in the literature. (2) Supplying a switching algorithm that prefers the best clock skew estimator between those obtained in [5,6] from the MSE perspective for the fGn case, where the Hurst exponent parameter is in the range of 0.5 ≤ H < 1. (3) Deriving new algorithms for the unexpected load and the cyber-attack detection mechanism with alarms for the fGn case (where 2 Mathematical Problems in Engineering 0.5 ≤ H < 1) that include the case in [7]. (4) A detection algorithm that supplies an alarm when detecting a clock skew that deviates with time.

System Description
Te IEEE1588v2 is a synchronization protocol based on a Master-Slave architecture. Te Master and the Slave exchange messages with time stamps for frequency and time synchronization. For a detailed system description, please refer to our previous works ( [5][6][7]). Let us recall Figure 1 from [5], where based on [2][3][4], we may write the following: Te TWD clock skew estimator described in [5] is given by where T l,j (i) for l � 1, 2, 3, 4 is Te OWD clock skew estimator for the Forward path is given by [6] and the OWD clock skew estimator for the Reverse path is presented by [6] In the following, we assume the fGn case, where the PDV is a zero mean process with the Hurst exponent parameter (H) in the range of 0.5 ≤ H < 1. Te PDV variance for j � m is and when j ≠ m, the PDV variance is given by [20,21] where n � 1, 2 for p � F, R, H F and H R are the Hurst exponent parameters for the Forward and Reverse, respectively, and E[.] denotes the expectation operator on (.). Based on (1) and (2), we may write the following: where for n � 1, 2, we have

The Switching Algorithm for the fGn Case
Te proposed switching algorithm needs the estimation of the amount of Sync message periods (J) to receive the system's requirement from the MSE perspective. But in order to be able to estimate the amount of required Sync message periods (J) correctly, we need to frst estimate the Hurst exponent parameters of the Forward and Reverse paths. After this is carried out, we have to estimate the PDV variances of the Forward and Reverse paths. Namely, the switching algorithm has three preliminary stages: (A) estimating the Hurst exponent parameters of the Forward and Reverse paths, (B) estimating the PDV variances of the Forward and Reverse paths, and (C) estimating the amount of Sync message periods (J) to receive the system's requirement from the MSE perspective.

Estimating the Hurst Exponent Parameters.
Te algorithm estimates the Hurst exponent parameters of the Forward and Reverse paths separately. It should be pointed out that based on (10), the random variables Ω 1,j (i) and Ω 2,j (i) are a subtraction between two fGn variables. Terefore, estimating the Hurst exponent parameter of Ω n,j (i) is not an easy task and cannot be carried out as proposed in [22]. Te authors in [22] estimated the Hurst exponent parameter only of a single fGn variable (and not the subtraction of two fGn variables). In this work, as already mentioned earlier, we assume that E[ω n [j]] � 0 for ∀j. Terefore, by using (10), the expectation of Ω n,j (i) for n � 1, 2 can be written as follows: Based on (10), the variance of Ω n,j (i) may be written as follows: For j � m and i � k, we may write (12) with the help of (7) and (8) for i ≥ 1 as follows: Based on (9) and (13), we may write the variances of Ω 1,j (i) and Ω 2,j (i) as follows: and Our Hurst exponent parameters estimation algorithm (estimating H F and H R ) is based on the following two cases of (12): case A and case B. In case A, j � m and i � k � 1. Based on (14) and (15), the variances of Ω 1,j (1) and Ω 2,j (1) for case A can be written as follows: and In case B, j � m, i � 1, and k � 2. Based on (7), (8) and (12), we can write the following for case B: Based on (9) and (18), we have the following for case B: and In the following, we present the three steps involved in the Hurst exponent parameter estimation algorithm, which are based on each other.
Step 1: (9) and (16)-(20), we may defne , and HH R 1 [j] as follows: where α is the estimated clock skew. Te estimated parameters 〈HH p 0 [j]〉 and 〈HH p 1 [j]〉 are given by and where β H is a small positive parameter. Please note that we set where J o is the initial amount of Sync message periods that have to be chosen by the user. In our simulation, we used for J o the values 400 and 500.
Step 2: in [22], the authors derived an algorithm for estimating the Hurst exponent parameter where ω was assumed to be an fGn noise. In (31)-(33), we recall the algorithm for estimating the Hurst exponent parameter based on [22] for a single fGn noise (ω) as follows: where ρ(k) is defned by and c(k) is Please note that when k � 0, we get c(0) � σ 2 ω . Based on (10), we have Ω 1,j (i) and Ω 2,j (i). Tose parameters are the outcome from subtracting two fGn variables from diferent Sync message periods. Terefore, in our work, we can not get a direct expression for the estimated Hurst exponent parameter via [22]. But, the main idea in [22] of dividing two cases (such as (29)) and getting from this outcome the estimated Hurst exponent parameter is still applied in this study. As mentioned above, our Hurst exponent parameter estimation algorithm is based on case A and case B. Now, by dividing these two cases, based on (16), (17) and (19)-(26), we can write the following: where H p [j] is estimated via (31), and L[j] is expressed by where β is a small positive valued parameter. Based on (16) and (17), we may approximate σ 2 ω 1 [j] and σ 2 ω 2 [j] as follows: and where are the initial estimated values. Based on (14) and (15), σ 2 ω 1 and σ 2 ω 2 are given by and where g 1 (i) and g 2 (i) are as follows: while for j > 1 and the OWD case, J[j] is derived via where T syn is denoted as the Sync message period, and E[e 2 ] is denoted as the system's requirement from the MSE perspective. 〈σ 2 〉 is a multiplication of the estimated PDV variance for the Reverse path with a compensation factor denoted as "s" and given by equation. Tus, we have 〈σ 2 Tis compensation factor compensates for the asymmetry between the Hurst exponent parameters for the Forward and Reverse paths. If there is a symmetry between the Hurst exponent parameters for the Forward and Reverse paths, then we have s � 1. Te compensation factor is defned by where Please note that we can also multiply the estimated PDV variance for the Forward path 〈σ 2 ω 1 [j]〉 by the inverse of the compensation factor (s). In this case, instead of (40) (for j > 1 and the TWD case or for j � 1), J[j] is derived via where 〈σ 2 ω 1 [j]〉 � 1/s〈σ 2 ω 1 [j]〉. C and D are given by [5] Mathematical Problems in Engineering where , and the function fG * H (.) is given by Please note that the parameter a is associated with the generalized fractional Gaussian noise (gfGn) case [21]. Tus, in this paper, we set a � 1 in order to have the fGn scenario. Please note that for a � 1 and H F � H R � 0.5, the Gaussian case is obtained for the Reverse and Forward paths. In addition, it should be mentioned out that in [5], the Hurst exponent for the Forward and Reverse paths were not estimated but assumed to be known and to be equal to each other. Namely, the Hurst exponent for the Forward path was assumed to be equal to the Hurst exponent associated with the Reverse path. Now, after we estimated (1) the Hurst exponent parameters for the Forward and Reverse paths, (2) the PDV variances of the Forward and Reverse paths, and (3) the amount of Sync message periods required to obtain the system's requirement from the MSE respective, we are ready to present the switching algorithm to the preferred clock skew estimator among the TWD clock skew estimator [5], the OWD clock estimator for the Forward path [6], and the OWD clock estimator for the Reverse path [6]. Our new proposed switching algorithm to the preferred clock skew estimator from the MSE perspective for the fGn case is based on [7] and given by where Z fGn , Z fGn , and σ 2 fGn are based on [6] and given by According to [5], A and B are as follows: where J � J[j] and Z[j] is defned as follows: where a is a small positive parameter. For the switching algorithm (46), we apply equation (42) for setting the value for s. Please note that unlike in [7], the switching conditions to the preferred clock skew estimator for the fGn case are based on the estimated Hurst exponent parameters for the Forward and Reverse paths, estimated in (32). In addition, the estimated PDV variances for the Reverse and Forward paths applied in (51) are a function of the Hurst exponent parameters related to the Reverse and Forward paths, respectively, and derived via (34)-(39). Since we deal with the fGn case, we use in (46) the parameter s (42), Z fGn (47), Z fGn (48), and σ 2 fGn (49) which makes again (46) very diferent from the obtained switching conditions obtained in [7] for the Gaussian case.

The Detection of the Cyber-Attack and the Unexpected Load
Te unexpected load and cyber-attack triggers alarms proposed in [7] are only applicable for the Gaussian case. In this section, we derive those alarms for the fGn case for 0.5 ≤ H < 1 where the Gaussian case is a special case with H � 0.5. Te indication of an unexpected load that is not defned as a cyber-attack will be denoted in the following as a "Load trigger [j]." Let us recall the unexpected load algorithm from [7]: and Te algorithm will alarm us that we have a load in the Forward or Reverse path based on a positive predefned parameter b. When the PDV variances of both paths are almost identical, the algorithm will alarm us that we have a load in both paths based on a positive predefned parameter d.
Te algorithm will alarm us that we have a suspect for a cyber-attack along one of the paths based on a positive predefned parameter c. Please note that we suspect a cyberattack in the Forward path when the Cyber − attack trigger[j] � 0.01 and cyber-attack in the Reverse path when Cyber − attack trigger[j] � 100. Te algorithm in (53) is diferent from [7], since Z[j] (51) as explained earlier is estimated with the help of the estimated Hurst exponent parameters for the Forward and Reverse paths.

The Detection of the Clock Skew Deviation
In our previous works ( [5][6][7]), we assumed that the clock skew (α) is constant, as is also assumed in other papers ( [2][3][4][8][9][10]). However, in a practical system, the clock skew may drift, as is explained earlier in this work. In this case, we defne the clock skew α[j] as follows: where f S and f M are the Master and Slave clock frequencies, respectively. Please note that for f M [j] � f M and f S [j] � f S ∀j, we have a constant clock skew. Since we assume that our PDV is a zero mean process, we may write based on (9) that where E[T 1,j− 1 (1)] � T syn and T syn as explained earlier is the Sync message period defned by the Master. Now, based on (55), when α[j] is constant, the expectation of T 2,j− 1 (1) (the Sync message period that the Slave measures) will remain approximately constant. However, when α[j] changes over time, the expectation of T 2,j− 1 (1) will drift. In the following, the expectation of T 2,j− 1 (1) is denoted as T 2,j (1). Te update equation for T 2,j (1) is where β T is a small positive valued parameter. Please note that when j � 1, is the time stamp after one Sync message measured by the Slave where is the initial expectation of T 2,j− 1 (1) and may be written as follows: Te difculty of tracking a drift in α[j] based on T 2,j (1) is that the PDV may change (increase or decrease) over time. On the one hand, when α[j] is changing over time and we have a constant PDV, we expect T 2,j (1) to change accordingly. On the other hand, when the PDV is changing and we have a constant α[j], we also expect T 2,j (1) to change. Tus, when T 2,j (1) changes, we can not determine for sure that there is a clock skew drift. In order to track clock skew drifts via T 2,j (1) due to drifts in α[j] without the PDV interference, we apply T 2,j (1) to indicate for any clock skew drifts caused by α[j] as follows: where sign(.) is the sign function and T 2,0 (1) � 0. Please note that the sign function in (58) compares the current value of T 2,j (1) with T 2,0 (1). If the current value of T 2,j (1) is above the initial expectation, it will increase T 2,j (1); otherwise, it will decrease T 2,j (1). In the following, we denote the indication of a clock skew deviation alarm as "Clock Skew Deviation Trigger [j]." Te algorithm for the clock skew deviation detection is given by where y is a positive predefned parameter. Please note that we suspect a clock skew deviation in the system when the Clock Skew Deviation Trigger[j] � 5.

Simulation Results
Tis section tests the algorithms that were derived in this paper.   (41)). We use J 0 Sync messages and the TWD clock skew estimator (3) for the initial stage. Te switching algorithm (46) begins to prefer the clock skew estimator with the lowest MSE only after the initial stage. Figures 2-6 show the performance of the TWD clock skew estimator (3), the OWD clock skew estimator for the Forward path (5), the OWD clock skew estimator for the Reverse path (6), and the performance of the fnal clock skew estimator chosen from (46) for the fGn case. Te switching algorithm (46) prefers the clock skew estimator, which achieves the system's requirement (MSE) with the lowest Sync message periods. According to Figures 2 and 3, the fnal clock skew estimator chosen by (46) achieves the lowest MSE. In case 1, the heavy load is in the Reverse path, and the higher Hurst exponent parameter is also in the Reverse path. In case 2, the heavy load is in the Forward path, and the higher Hurst exponent parameter is also in the Forward path. It can be seen according to Figures 2 and 3 that the algorithm prefers the OWD clock skew estimator for the Forward path (case 1) or the OWD clock skew estimator for the Reverse path (case 2). Please note that for j < 2000, the PDV variances are the same on both paths. But, since the Hurst exponent parameters for the Forward and Reverse paths are diferent, the path with the lower Hurst exponent parameter is chosen. Tus, the switching algorithm chooses the OWD clock skew estimator over the TWD one. According to Figure 4, the fnal clock skew estimator preferred by (46) achieves the lowest MSE for the unexpected load in both paths (case 3). According to Figures 5 and 6, the fnal clock skew estimator chosen by (46) achieves the lowest MSE. In case 4, the heavy load is in the Reverse path. However, the higher Hurst exponent parameter is in the Forward path. In case 5, the heavy load is in the Forward path. However, the higher Hurst exponent parameter is in the Reverse path (case 5). It can be seen according to Figures 5 and 6 that for j < 2200, the algorithm prefers the OWD clock skew estimator for the Reverse path (case 4) or the OWD clock skew estimator for the Forward path (case 5). Please note that for j < 2000, the PDV variances are the same on both paths. But since the Hurst exponent parameters for the Forward and Reverse paths are diferent, the path with the lower Hurst exponent parameter is chosen. Tus, the switching algorithm chooses the OWD clock skew estimator over the TWD one. For 2200 < j < 3000, the algorithm prefers the TWD clock skew estimator (case 4 and case 5) since in one path, we have an unexpected load, and in the second path, we have the higher Hurst exponent parameter. For j > 3000, the algorithm prefers the OWD clock skew estimator with the path with the higher Hurst exponent parameter (the OWD for the Forward path in case 4 and the OWD for the Reverse path in case 5) since we have in the second path a heavy load compared to the other path.

Te Detection of the Cyber-Attack and the Unexpected
Load. Te algorithms (52) and (53) supply the alarms for the unexpected load or cyber-attack in one of the paths. Figures 7-10 show the unexpected load and cyber-attack output triggers and Z[j] (51) by which we determine the output triggers. In Figure 7, we set H F � H R � 0.5, meaning we have the Gaussian case, where the heavy load is in the Forward path. According to Figure 7, the algorithms (52) and (53) detect the unexpected load and the cyber-attack in the right path and output their alarms accordingly. Note that this corresponds to the case in [7]. In Figures 8-10, we have the fGn case for 0.5 ≤ H < 1. Figure 8 shows that we have a heavy load in both paths. According to Figure 8, the algorithm (52) detects the unexpected load and outputs its alarm accordingly. In Figures 9 and 10, we have a diferent Hurst exponent parameter in each path. In Figure 9, we set H F � 0.8 and H R � 0.6 and a heavy load in the Forward path, and in Figure 10, we set H F � 0.6 and H R � 0.8 and a heavy load in the Reverse path. According to Figures 9 and 10, the algorithms (52) and (53) detect the unexpected load or the cyberattack in the right path and output their alarms accordingly.

Te Detection of the Clock Skew Deviation.
Te algorithm (59) supplies the alarm for clock skew deviation. Figures 11-16 show the trigger clock skew deviation output, the T 2,j (1) (58) by which we determine the output trigger. In order to track the actual α[j], we apply α � [j] a normalized clock skew defned as follows:

Mathematical Problems in Engineering
In Figures 11-16, we have three cases: (1) a constant clock skew value, (2) an increase in the clock skew over time, and (3) a decrease in the clock skew over time. In Figures 11-13, we test the algorithm (59) in those three cases without an increase in the PDV. Figures 14-16 test the algorithm (59) in those three cases with an increasing heavy load in the Forward path. According to Figures 11  and 14, the algorithm (59) detects that there is no clock skew deviation and outputs its alarm accordingly. According to Figures 12, 13, 15, and 16, the algorithm detects the clock skew deviation and outputs its alarm accordingly, even with an increasing or constant heavy load in the Forward path.

. Discussion
As it was already stated earlier in this paper, we can measure indirectly only the PDV sample in the PTP case together with a time-dependent ofset. In order to see this, we rewrite (1) and (2) as follows: Please note that we have on hand only the PTP timestamps. In other words, we have on hand only t n [.] for n � 1, 2, 3, 4. Since we do not know the exact value for α, we substitute for α some estimation. Let us consider for a moment the TWD case, thus α is defned as α. Terefore, the estimation error is Now, by using substituting (63) into (61) and (62), we have   6) words, those ofsets increase or decrease with time depending on the sign of e. Terefore, based on (64) and (65), we can understand why the classical climacogram algorithm [12] cannot be applied here for e ≠ 0. Please note that in our approach in estimating the Hurst exponent, we look at the diference between PTP timestamps (T n,j (1) and T n,j (2), where n � 1, 2, 3, 4). Tus, we do not have to deal in our algorithm with the problem discussed from above of having an error that increases or decreases with time.
In subsection 3.1 of Section 3, we derived the estimation for the Hurst

Conclusion
Tis paper proposed a switching algorithm to the preferred clock skew estimator for the PTP case in a practical system (fGn case where 0.5 ≤ H < 1). Te switching algorithm chooses the clock skew estimator (TWD, OWD for the Forward path, or OWD for the Reverse path) with the best performance from the MSE perspective. Te switching algorithm uses an algorithm that estimates the PDV variances and Hurst exponent parameters associated with the Forward and Reverse paths. In this work, it is the frst time where the estimation of the Hurst exponent parameters is carried out for the Forward and Reverse paths based on time stamps related to the PTP case. Simulation results have confrmed that our proposed switching algorithm may work well in a practical network where asymmetry exists in the Hurst exponents, PDV variances, and fxed delays associated with the Forward and Reverse paths. In addition, we proposed in this paper also two detection algorithms (one for unexpected loads and the other for cyberattacks) that detect unexpected loads or cyber-attack in one of the paths or unexpected loads in both paths for the fGn case where 0.5 ≤ H < 1. Tose two algorithms are based on the estimation of the PDV variances and Hurst exponent parameters associated with the Forward and Reverse paths and output trigger alarms when unexpected loads or a cyber-attack is detected. Simulation results have confrmed the efectiveness of those algorithms. Tis paper also proposed an algorithm that detects clock skew deviation over time, which occurs in practical Master or Slave clocks. Tus, the system designer may have an additional tool for diagnosing the system.

Data Availability
All the data supporting the results are already included within the manuscript.

Conflicts of Interest
Te authors declare that they have no conficts of interest.