Min Xie Cyber-Physical Distributed Systems

3.3.1.2 Parameter Estimation of DHMM

The parameters (π, P, B) of DHMM for the network state and the time delays in (16) need to be estimated using the set of time delay observations o = {o1, o2, ⋯, oK} obtained in the K periods. To this aim, the Expectation Minimization (EM) is widely used for its efficiency in finding the maximum likelihood estimation (MLE) of the parameters of a statistical model given observations [632]. Therefore, the MLE of the parameters (π, P, B) based on the EM algorithm is:

(3.35)

where Pr(o| λ) is the likelihood function which needs to be maximized during the EM algorithm.

The EM algorithm aims to find optimal λ* for the DHMM which can most likely generate the set of the time delay observations o = {o1, o2, ⋯, oK}. Pr(o| λ) in (3.35) is actually an incomplete data likelihood function, i.e., does not employ the underlying data (network state q) in the optimization process, and can be rewritten as lg Pr (o| λ) due to the monotonicity of logarithmic function.

To solve the missing data problem, based on the MDEM‐based Baum‐Welch algorithm, the hidden network state q is considered, and a complete data log‐likelihood function lg Pr (o, q| λ) is achieved. Thus, the cost function used in the Baum‐Welch algorithm is given as

(3.36)

where λ′ is the optimal parameter estimation obtained in previous iteration and used to evaluate the expectation. λ is the new parameter estimation to be optimized for the purpose of increasing G.

As the underlying network state q is a discrete sequence which takes values from Q, G can be modified as

(3.37)

As the cost function in (3.37) is maximized only with respect to λ, i.e., Pr(o| λ′) can be neglected, the likelihood function can be rewritten as

(3.38)

If G(λ, λ′) ≥ G(λ′, λ′), Pr(o| λ) ≥  Pr (o| λ′) holds.

Proof: Based on [633][634], the following relationship holds

(3.39)

Thus, when G(λ, λ′) ≥ G(λ′, λ′), there exists log Pr {o| λ} − log  Pr {o| λ′} ≥ 0, i.e., Pr{o| λ} ≥ Pr {o| λ′}.

Theorem 3.1 illustrates the reason that the MDEM‐based Baum‐Welch algorithm can maximize the likelihood function of the DHMM for time delays. Each iteration in the MDEM‐based Baum‐Welch algorithm consists of two steps, which are given as follows

• E‐step: compute

(3.40)

• M‐step: find a local maximum λ* that maximizes G(λ, λ′), i.e.,

As shown by (3.40), the optimization problem of maximizing G(λ, λ′) are independently split into three problems of deriving maximum values for three terms: π*, P* and B*. These values are achieved by solving its derivative via three procedures, i.e., forward procedure, backward procedure and update procedure, given by the MDEM‐based Baum‐Welch algorithm [635]. The results of π*, P* and B* for the DHMM have already been given by [636].

These two steps, i.e., E‐step and M‐step, are repeated until a desired level of convergence, i.e., |λ* − λ′| < ε, is reached. The final result would converge to an extremum of the likelihood function.

However, as λ consists of π, P and B, it is very hard to compute the distance between λ* and λ′. Notice that

If

Proof: The proof can be referred to [638][639].

Therefore, when |λ* − λ′| < ε is reached, the MDEM‐based Baum‐Welch algorithm halts. The most likely hidden state sequence is computed and used in the time delay prediction.

3.3.1.3 Delay Prediction via DHMM

Denoting the prediction of τk as

The technique of quantifying

• Use λ* and

(3.41)

• Use λ* and

(3.42)

• Based on the uniform quantization technique and

(3.43) i.e.,

The technique of quantifying

• The first two steps are same as the steps presented in (3.41) and (3.42).

• Based on the K‐means clustering quantization technique and

(3.44) i.e.,

In the period K, the control center has S‐C time delay measurements

Finally, by using the Viterbi algorithm, the most likely network state sequence

In the EWMA method,

(3.45)

where α is the weight and n is the window length.

3.3.1.4 Smith Predictor Structure

The delayed control process can be effectively handled by the Smith predictor if the information on time delays is known [641][642]. In this chapter, the time delays of the network τk are measured by the timestamp technique of Section 2.3, and

Figure 3.8 The illustrative WAPS with the DHMM‐based Smith predictor.

As shown in Figure 3.8, P(s) denotes the transfer function of the real plant in Figure 3.1 and equals to GG(s)GT(s)GP(s).

(3.46)

The Smith predictor uses an internal model,

The estimated disturbance is fed back through a filter Fd(s) to improve the disturbance rejection performance. Based on [643], the Fd(s) is as follows:

(3.47)

where n is the number of the poles of

Defining the poles to be canceled as r1, ⋯, rn, the coefficients a1, ⋯, an satisfy:

(3.48)

Finally, the discrete PID controller is driven by a delay‐free signal

The application of the proposed DHMM‐based Smith predictor to the WAPS with communication delays is represented in Figure 3.9.

Figure 3.9 The flowchart of the DHMM‐based Smith Predictor.

3.3.2 Delay Predictions

3.3.2.1 Settings of DHMM

The interval (hL, hU] of time delays is divided into seven complete subintervals based on the dataset of time delays [647], i.e., M = 7 and O = {1, 2, 3, 4, 5, 6, 7} [648]. Given the parameters for the Ethernet in Table 3.1, the C‐A time delays are comprised within the interval (0.048,0.158] and the S‐C time delays are within the interval (0.037,0.126]. For the EWMA, α = 0.8 and n = 10 [649].

Applying the uniform quantization technique to the C‐A time delays, (3.27) is determined as follows:

(3.49)

Based on the MDEM‐based Baum‐Welch algorithm, the DHMM parameters estimated from one sample containing 100 C‐A time delay measurements, i.e., the length of most recent time delays is 100, are

(3.50)

(3.51)

(3.52)

where

are the estimation results derived from the uniform quantization technique.

In order to make comparisons, the estimation results

derived from the K‐means clustering quantization technique are:

(3.53)

(3.54)

(3.55)

where the seven centroids are: 0.072, 0.082, 0.086, 0.092, 0.102, 0.112 and 0.120.

3.3.2.2 Prediction Comparison

The results of the estimation processes show that

and

are different, because different quantization techniques produce different observations {o1, ⋯, o100}.

Figure 3.10 displays the observations for the periods 30 to 44. In the period 40, the delay measurement, i.e., 0.084, is comprised within the third subinterval (0.080,0.096] based on the uniform quantization technique, and, therefore, o40 = 3. Conversely, the same delay measurement, i.e., 0.084, is closer to the second cluster centroid, i.e., 0.082, for the K‐means quantization technique, and, therefore, o40 = 2.

Figure 3.10 The C‐A time delay observations derived from different quantization techniques.

Therefore, different quantization techniques may produce different observation values for the same delay measurement. Since the centroids in the K‐means quantization technique are statistically inferred from the distribution of time delays, the observations derived from the K‐means quantization technique can better track the actual time delays as compared to the uniform quantization technique as shown by Figure 3.10.

The tracking ability of the K‐means quantization technique results in more accurate predictions as compared to the uniform quantization technique. Based on the estimated

,

and the proposed time delay prediction method defined by (3.41) ‐ (3.44), the corresponding prediction results of the EWMA and DHMM method are shown in Figure 3.11.

Figure 3.11 The prediction results for C‐A time delays

(Source: Mo and Sansavini, 2021 [650]).

The mean square error (MSE) of the EWMA, Uniform‐based DHMM and K‐means‐based DHMM are 0.0636, 0.0563 and 0.0478, respectively. The prediction of EWMA lags behind the real time delay, because it is based on the moving average of actual time delays. The prediction of DHMM relies on the most likely communication network state and stays constant until the network state changes. Figure 3.11 shows that the K‐means‐based DHMM provides more accurate predictions than the uniform‐based DHMM because the observations of the K‐mean quantization technique are statistically inferred.

Figure 3.12 displays the most likely network states estimated by the Viterbi algorithm. For the uniform quantization technique, the most likely network state in period 40 is q40 = 2. Based on (3.41) and

, the predicted network state in period 41 is

. Based on (3.42) and

, the predicted observation in period 41 is

. Therefore, the predicted time delay in the period 41 determined by (3.43) is 0.0873s.

Figure 3.12 Most likely network states derived from the two different quantization techniques.

For the K‐means quantization technique, the most likely network state of period 40 is q40 = 2, and, thereby, the predicted network state of period 41 is

, based on (3.41) and

. Based on (3.42) and

, the predicted observation of period 41 is

, and, therefore, the predicted time delay in period 43 belongs to the third centroid, i.e., 0.0859s, based on (3.44).

Table 3.5 shows the statistical properties of the MSE of the EWMA and DHMM method for C‐A and S‐C delays, computed from 1000 samples each containing 100 delay measurements to ensure accuracy.

Table 3.5 Statistical properties of MSE of delay predictions.

The DHMM outperforms the EWMA in the delay prediction because of the smaller mean MSE. Compared with the EWMA, the uniform‐based DHMM reduces the MSE of C‐A and S‐C delay predictions by 3.47% and 3.72%. Compared with the EWMA, the K‐means‐based DHMM reduces the MSE of C‐A and S‐C delay predictions by 7.19% and 7.26%.

Indeed, the prediction of EWMA depends on recent delays while the prediction of DHMM depends on the current network state, which is the underlying process causing delays. Because the observations of K‐means‐based DHMM is able to capture the statistical distribution of time delays, the MSE of the K‐means‐based DHMM is smaller than the one of the uniform‐based DHMM.

3.3.3 Performance of Smith Predictor

3.3.3.1 Settings of Smith Predictor

At t = 35s, a positive load disturbance of 0.1 pu is added to the WAPS with communication delays described in Table 3.1 [651][652]. We study two different cases, i.e., the LFC via discrete PI controller derived from the delay margin method [653] and LFC via robust PID controller derived using genetic algorithms (GA) [654]. These two controllers are used by most real power systems [655][656][657][658][659][660][661][662][663][664].

The internal plant model

with droop characteristic is

(3.56)

where

has a pair of complex poles −1.302 ± 2.184i.