The physical grid measurement data reflect the real-time operating status of the grid under different working conditions. The measurement data of each section of the grid reflect the operating status at that moment. We do not consider the reasons for changes in the grid state (caused by cyber-attacks or general equipment failures); it can be described as a specific interval _Δt(t₁∼t_n). According to the acquisition sequence, all state data fragments S_p(t_i) consisting of the physical grid operating data link _Qp(Δt) are defined as follows:

where X_p(t_i) represents the h measured attributes obtained from the physical-side device X_p at time t_i, including the voltage, current, phase angle, active power, and reactive power; S_p(t_i) represents all the measurement data collected by m devices on the physical side at time t_i.

The transmission delay and data packet loss rate typically reflect the performance status of a cyber-network. When the control signal or status information is lost during the transmission of the data packet because it exceeds the allowable proportion, the control of the device has been lost due to a network attack (Davarikia and Barati, 2018; Wang Q. et al., 2019). Three indicators (delay rate, packet loss rate, and threat degree) are established to characterize the operating data link of the cyber network.

where n is the number of communication links transmitting data, P′ is the number of data packet losses for link k, P_k is the number of data packets sent by link k, and P_T is the threshold of the packet loss rate of link k.

where m is the number of devices that send information, T_l^send is the sending time of data packet l, and T_l^receive is the receiving time of data packet l.

where w_ij is the threat degree of alarm events a_i,j, w_i is the average value of the threat degrees of all alarm events in IP_i, and n_i is the number of all alarm events in IP_i.

The three performance indicators of the operating status of the cyber network are used to establish the cyber system operating data link Q_c(Δt) in the interval Δt(t₁∼t_n):

where Y_c(t_i) represents the R_dr, R_pr, and W_th obtained from the cyber-side device Y_c at time t_i; S_c(t_i) represents the status data obtained from k devices on the cyber side at time t_i.

We use topological mapping to couple and map the data links of the two heterogeneous networks to form a data link of the cyber–physical operating state. The grid can be divided into m areas according to the physical grid connectivity, and each area has n transmission lines. It is assumed that a line consists of k electrical components {X₁, X₂, …, X_k}, each line is connected to n communication devices {Y₁, Y₂, …, Y_n}, and each communication device has a unique IP address {IP₁, IP₂, …, IP_n} in the cyber network. We sequentially connect each electrical component number, line number, and connected area in the data chain to create an index table linking the <connected area number Area, line number Line, electrical component ID number, and information component IP address>. The cyber network operating data link Q_c and the physical power grid operating data link Q_p in the interval are compared using the index table, and the data are stored in the corresponding index.

The cyber network clock with a collection period of T is used, and we set the sampling time window to ε = [T−αT′,T], where _α is the window size parameter. The larger the value, the longer the collection period is. In the sampling time window _ɛ many identical state events may occur in the cyber–physical coupling state chain. Therefore, these repetitive events are filtered and compressed to form the cyber–physical operating state data link, which is expressed as follows:

There are no labels for the different state categories in the original cyber–physical operating state data link Q(ε). It is necessary to distinguish the different state categories using cluster analysis. In this paper, the two-step principal component analysis (PCA) clustering algorithm is proposed. The PCA algorithm is used to cluster, transform, and filter the correlated attributes to extract linear uncorrelated attributes (Jian et al., 2004). The two-step algorithm is used to cluster the attribute set; it reduces the computational complexity and provides high clustering accuracy (Dom et al., 2003; Northrup et al., 2004; Phelps et al., 2009). The algorithm steps are as follows:

Input: cyber–physical operating state data link Q(ε) = {x₁, x₂, …, x_n}.

Output: D = {x_i, C_i}, where C_i is the operating state of the clusters C = {C₁, C₂, …, C_k}.

Step 1: Feature selection for clustering. The PCA algorithm is used to map n attributes in the data link to m dimension (m < n). The correlated attributes are filtered using an orthogonal transformation to obtain m-dimensional new features, A = {A₁, A₂, …A_m}. The centralizing mean xi=xi-1n⁢∑i=1nxi is used to derive the covariance matrix _XX^T, whose eigenvalues and eigenvectors are obtained. The data link set Q’(ε) is obtained after dimensionality reduction.

Step 2: Calculate the number of category clusters in the operating state. After the traversal process, the clustering feature (CF tree) growth in the balanced iterative reducing and clustering using hierarchies (BIRCH) algorithm is applied to the data link set Q’(ε). The data points in the data set are evaluated one by one to collect all data points in the dense area while generating the CF tree. The log-likelihood distance d(C_s,C_t) = ζ_s + ζ_t−ζ_{< s,t >} between the two clusters is used to create many small subclusters. The Bayes information criterion (BIC) is used to calculate the number of possible division schemes for the state category.

Step 3: Determine the number of categories C_J in the Q’(ε). The agglomerative hierarchical clustering (AHC) method is used to merge the subclusters one by one, and the desired number of clusters is reached according to the R(k) between the two clusters.

where C_k and C_k+1 is a partition scheme with k or k+1 cluster numbers; d_min(C_k+1) and d_min(C_k) is the distance between the two smallest clusters in the scheme.

Step 4: Label the sample data in each operating state cluster. The data points in each cluster are determined; the data points x_i in the state data link set Q’(ε) are regarded as single-point clusters according to the clustering results C_J. The logarithm similarity between x_i and each cluster in C_J is determined. Given the distance d{{x_i}, C_J}, x_i is placed into the nearest cluster, and labels are generated for each operating state category C = {C₁, C₂,…, C_k}.

A coordinated cyber-attack event of the CPPS has a small probability and high risk. In the data link _Q(ɛ) normal operation data account for the largest proportion, whereas the proportion of attack data is relatively small, resulting in unbalanced data. Therefore, the ADASYN algorithm is used to deal with the imbalance of the operating state classes (Qu et al., 2018; Wang et al., 2020). Balanced data distribution is obtained by adaptive synthetic oversampling. Different minority samples are given different weights to generate different numbers of samples. The algorithm process is as follows:

Input:D = {x_i, C_i}, where x_i is the cyber–physical operating state data link Q(ε), C_i is the class label. α is the imbalance threshold, C_k is a minority class, and C_l is the majority class.

Output: Balanced data set D′.

Step 1: Calculate the class imbalance, where Imbalance=Lagre⁢num⁢(Cl)Small⁢num⁢(Ck). Calculate the number of samples to be synthesized based on the degree of imbalance G = (Lagrenum(C_l)−Smallnum(C_k))×β,β∈[0,1].

Step 2: Calculate the proportion of the majority class in the K-nearest neighbors (KNNs). r_i = Δ_i/K, where Δ_i is the number of samples of the majority class in the KNN.

Step 3: Calculate the majority class surrounding each minority sample.

Step 4: Calculate the number of samples that need to be generated for each minority sample C_k.

Step 5: Select a minority sample among k neighbors around each minority sample and synthesize using Eq. (10). Repeat the synthesis g_i times until the desired number of synthesized samples is obtained.

where s_i is the composite sample, _{X_i} is the i-th sample in the minority sample, X_i ∈ [0,1], X_zi is a randomly selected minority sample among the KNNs of X_i.

Repeat the synthesis until the desired number of synthesized samples in Eq. (5) has been obtained.

The purpose of attack detection is to minimize the harm to the power grid caused by the attack. The harm caused by misinterpreting an attack as a normal event is far greater than that caused by misinterpreting a normal event as an attack (Huang et al., 2018). We propose using the CS function to improve the GBDT (Sakhnovich, 2011; Liao et al., 2016). The CS loss function replaces the standard cost loss function to prevent attack event misclassification. The improved CS loss function is defined as follows:

where K is the class of all attacks, _{C_k} is the sample of the k-th attack, and _{p_k (x)} is the probability of the k-th attack, _{w_k} is the CS function, it can be divided into two costs, i.e., the missed detection cost w(-,+)⁢(1-p⁢(x)),p⁢(x)≥w-w++w- and the misdetection cost w(+,-)⁢p⁢(x),p⁢(x)<w-w++w-.

Coordinated cyber-attack detection is a multi-classification task. A total of K types of attacks are assumed. The sample x in the cyber–physical operating state set is obtained, and the CS-GBDT algorithm is used to determine which class the x sample belongs to. The specific steps of the algorithm are as follows:

Input: Balanced data set D=′{(x1,C1),(x2,C2),…,(xN,CN)}, loss function Loss (C_k, f_k(x)), and the number of classifiers M.

Output: A strong learner for attack classification F(x).

Step 1: Initialize f_k0(x) = 0, the number of categories classified k = 1,2,…K.

Step 2: Starting from t = 1 to t = M, there are M iterations in total, repeating steps 3 through 6, at last building M classifiers.

Step 3: The one-hot code for each class y_i is generated. We calculate the probability of sending the k-th attack sample p_k(x).

Step 4: Start from k = 1 to k = K, repeating steps 5 through 6, we generate K different CART classification trees f₁(x), f₂(x),…f_K(x).

Step5: Calculate the negative gradient of each class in the m class and obtain the negative gradient error of the i-th sample corresponding to category k in the t-th iteration:

where N is the number of sample data.

We use the estimated residual {(x₁, r_k1),…(x_N, r_kN))} as an input to calculate the leaf node area of the m-th decision tree:

where R_mkj is the leaf node region R_mj of the m-th tree. K is the number of categories.

Step 6: Update the classifier f_mk(x).

where J is the number of leaf nodes per tree.

Step 7: Build final classification tree with high accuracy used for attack detection.

We simulate the different fault states of the physical power grid caused by cyber-attacks on the IEEE39-bus system in the RT-LAB and OPNET co-simulation environment. We collect the DR, PR, and threat information at different times on the cyber side. The voltage, current, impedance, and other data are collected on the physical side. The 10 data sets are obtained. Each set contains 56 attributes, and the cumulative number of records is about 50,000, including five types of operating states in the CPPS, as follows:

(1) Normal operating state (S1): there is no network attack on the cyber side, and the power grid on the physical side is operating normally. (2) Distributed denial of service (DDOS) attack state (S2): the data in the communication system are blocked by a DDOS attack, affecting the normal operation of the power system, measurement acquisition, and control commands. (3) Data injection attack state (S3): malicious data injection into physical power grid disguised as a normal fault, resulting in the operator mistakenly assuming a short-circuit fault. (4) Protection device parameter tampering attack state (S4): the attacker tampers with the distance parameter of the protection device, causing a failure of the protection device to disconnect the fault area. (5) Fault operation state (S5): the physical power grid has a single-phase, two-phase, or short-circuit fault.

The data set 1 with 4,966 records is selected in the experiment. After implementing the two-step PCA clustering algorithm, the number of outliers is 89, and there are five operating states, as shown in Figure 2A. Cluster-3 (S2) has the largest number of records (1,325). The clustering superiority is 0.93, and the clustering importance is 0.85, accounting for 29.3% of the records. The smallest cluster is Cluster-4 (S5), with 59 records, accounting for 1.3% of all records. The clustering superiority is 0.97, and the clustering importance is 0.93. Clustering superiority is a measure of cluster separation (−1∼0.2 poor| 0.2∼0.5 medium| 0.5∼1 good), and clustering importance is a measure of cluster cohesion (0∼0.2 poor| 0.2∼0.6 medium| 0.6∼1 good) (Nair and Narendran, 1997).

According to the negative sequence current and zero sequence current amplitude of each cluster in the experiment, the curves of the three attack states and the fault state are obtained, as shown in Figures 2B–E. Cluster-2 is significantly different from the other four states, while Cluster-4 and Cluster-5 have high similarities. The reason is that Cluster-2 is an attack that causes network blocking and delay, which is significantly different from the other types of data tampering attacks. Cluster-3 and Cluster-5 are physical power grid failures caused by information tampering attacks. These states are similar to the changes occurring in the Cluster-4 power grid normal fault.

The adjusted Rand index (ARI) is used to measure the accuracy of the clustering results; ARI ∈ [−1,1], the closer the value is to 1, the better the clustering performance is. The index is calculated as follows:

where RI is the Rand coefficient, and E(RI) is the expected value of each class.

Four typical clustering algorithms are selected for performance comparison, i.e., K-means, density-based spatial clustering of applications with noise (DBSCAN), clustering using representatives (CURE), and BIRCH. In the experiment, the sample size of the test data set is randomly selected and ranges from 5 to 100% of the data set. The ARI values of the different algorithms are shown in Figure 2F. As the proportion of the test data set increases, the ARI increases significantly. The accuracy of the proposed two-step PCA method is 97% for a sample size of 100%, demonstrating the excellent performance of this method. The K-means algorithm has the lowest ARI values.

The number of samples in the operating state classes in 10 data sets before implementing the algorithm is shown in Figure 3A. The number of samples is imbalanced in the different operating states. The largest number of records (143,766) occurs in the S3 state, and the fewest number (3,080) is observed in the S5 state. The maximum class imbalance is 3.77. There are multiple minority and majority categories in the joint data set, showing multi-category imbalance.

The ADASYN algorithm is used to oversample the categories whose number is less than the threshold. We set the maximum imbalance threshold to 1.2. The results in the 10 data sets are shown in Figure 3B. The proportion of records in each dataset is close to 20%. The ADASYN algorithm uses local screening and sampling to reduce the influence of data imbalance on the false alarm rate of coordinated cyber-attack detection.

The balanced data set is divided into a training set (70% samples) and a test set (30% samples). The model loss parameters are set according to the improved CS loss function. There are 130 integrated base classifiers, and the depth of each independent tree (max_depth) is seven.

The receiver operating characteristic (ROC) curve obtained by classifying the test data set is shown in Figure 4A. The curves of the five categories are close to the (0,1) position, and the average area under the ROC curve (AUC) is 0.982. This result shows that the attack detection model has a low false alarm rate and high accuracy.

The precision-recall curve obtained by classifying the test data set is shown in Figure 4B. The precision-recall curve are all close to the (1,1) position, indicating that the attack detection model has high recall and accuracy, even when the ratio of positive and negative samples is large. Therefore, the proposed attack detection model has a high classification accuracy for unbalanced data.

The confusion matrix of the attack detection results is shown in Figure 4C. The detection accuracy for the DDOS blocking attack (S2) is 98%, that of the data injection attack (S3) is 96%, that of the protection device parameter tampering attack (S4) is 97%, that of the normal operation (S1) is 99%, and that of the fault operation (S5) is 98%. These results demonstrate that the proposed coordinated cyber-attack detection model accurately detects coordinated attack events on the network and distinguishes attack states from the fault operation state, with a maximum false-positive rate of only 4%.

Finally, the proposed model is compared with typical classification algorithms, including the KNN, Xgboost, Random Forest, Adaboost, and support vector machine (SVM). The overall accuracy, average recall, average precision, and average F1-score of the algorithms are shown in Figure 4D. The recall and precision of the CS-GBDT algorithm are higher than 97%. The algorithm performance is stable, and it provides better performance for detecting various attack events than comparable algorithms.