Cascading Failures Analysis Considering Extreme Virus Propagation of Cyber-Physical Systems in Smart Grids

1School of Electrical Engineering and Electronic Information, Xihua University, China 2Key Laboratory of Fluid and Power Machinery, Ministry of Education, Xihua University, Chengdu 610039, China 3School of Electrical Engineering, Southwest Jiaotong University, China 4Department of Energy, Politecnico di Torino, Italy 5Research Group on Natural Computing, Department of Computer Science and Artificial Intelligence, University of Seville, Spain


Introduction
The smart grid, as a modern electrical network (EN) infrastructure, can enhance the efficiency, reliability, and security of traditional ENs based on the advancement of cyber-physical systems [1][2][3]. In a smart grid, the monitoring, control, and management of the EN depend closely on the smart information and communication (cyber) network [4][5][6], which works such that the EN ensures not only its own secure operation but also reliable operation of the entire communication network. Meanwhile, when the EN fails (especially, through cascading failure), fault cross-propagation between the electrical and communication networks (ECNs), called interdependent network, occurs, which increases the complexity of fault propagation owing to interactions between the ECNs. For example, the Italian blackout of 2003 was triggered by effects of the ECN [7]. Therefore, exploring the propagation mechanism of interactive cascading failures [8] in an interdependent network has been receiving increasing attention.
To date, the connection between the different coupling modes between ECNs and the robustness/vulnerability of interdependent networks has been investigated widely [9][10][11][12]. Studies have demonstrated that the different types of links between ECNs greatly impact the robustness of the network. For example, [10] reveals the double-network link allocation 2 Complexity strategy is superior to single-network link allocation strategy. Therefore, reasonably allocating the interconnecting links between ECNs is vital for improving the robustness of interdependent networks. Accordingly, a few models (e.g., Petri nets) have been introduced to reveal the mechanism of interactions leading to catastrophic blackouts [13]. However, these works have been done merely from the perspective of the structure of the coupling of the ECNs.
Meanwhile, the physical and operational characteristics considering the interactions have been focused on as well. In [14], the impact of communication network vulnerability on power system operation was assessed considering both latency and communication interruptions. The data exchange model is introduced for modeling cascading failures in interdependent networks [15]. Reference [16] proposes a simulation platform to analyze the ECN vulnerability by considering the control strategy of power balance. Similarly, other simulations [17,18] have been proposed to analyze the fault mechanism considering interactions between ECNs. Although these studies have included the interactions between the communication network and EN, their focus is only to study how to set up the simulation platform by considering both the communication network and the EN. Moreover, the operational characteristics are only studied from a steady state point of view.
In addition, as communication networks become increasingly smart and as smart grids are increasingly accessed using the Internet owing to advancement of the energy Internet [19, our study, we improved the dynamic power flow method to redistribute loads and adjust unbalanced power in the network during fault propagation by introducing the primary frequency regulation and the equations of rotors of generators. Notably, we only consider high-voltage transmission networks as the study objects.
For the virus propagation, virus spread models [30,31] with time delay have been developed based on the complex network theory (CNT) from the perspective of the topological structure of the communication network. Generally, the communication network mainly has two topical topological structures: scale-free and small-world networks. This paper focuses on investigating the impacts of two types of networks on electrical networks during virus propagation. In other words, we analyze which type of communication network is more vulnerable from the network-wise perspective, i.e., once viruses are propagated in a communication network, which type of communication network can cause more damage to the electrical network. Meanwhile, we further analyze that in a communication network, the vulnerability of communication nodes with different degree is revealed by investigating the number of fault branches and blackout level of coupled electrical networks.
In addition, it should be noted that the most modern malware, surely most malware is used in known attacks to power grids and industrial controls, limits its own effectiveness by prematurely destroying/disabling nodes and is not self-replicating. Furthermore, currently in a real-world communication network vertexes will not be homogeneous; thus they will not support most of self-propagating malicious code. However, in order to consider a low-probability but high consequence scenario, in this paper, we assume a random constant spread malicious code (called "virus" in this paper) with the following features to investigate the impact of extreme case of self-propagating virus on the power system from the network-wise perspective: (1) the virus can block the communication between infected vertexes and the control center; (2) the virus can self-propagate among homogeneous vertexes; (3) a few infectious vertexes can be cured with the probability owing to the strengthening of security measures; (4) the differences of security level of each node are not considered.
The remainder of this paper is organized as follows. Section 2 describes the interactions between the ECNs in the coupling relationships and topological structures. In Section 3, the virus propagation models with time delay and information exchange model in the communication network are introduced. The dynamic power flow method and the overload mechanism are established in Section 4. The cascading failures model considering the interactions and the corresponding simulation analysis are described in Sections 5 and 6, respectively. In Section 7, we further discuss the contribution of this paper and the external validity of the modeling. Finally, conclusions are given with possible future work in Section 8.  is considered for exchanging information packets with the control center [15].

. . Topological Structures of ECN
Electrical and Communication Network as Graphs. For simplifying analysis of the topological structures, we abstracted the ECN as graphs. The EN can be considered as a complex network with nodes and links. The buses, including generators, loads, and substations, can be viewed as nodes while transmission lines can be viewed as branches; therefore, the electrical network is represented as the graph G = (N, B). The adjacent matrix = ( ) × is employed to define G as follows: where = represents that there is a branch between nodes and .
Similarly, the optical fibers and SDHs of the communication network can be considered as edges and vertexes, respectively; therefore, the communication network is represented as the graph G = (V, E). The adjacent matrix = ( V ) × is employed to define G as follows: where V = represents that there is an edge between vertices and V . Topological Structure Analysis. We analyze the structural characteristics from the perspective of the CNT. Existing literature indicate that ENs have small-world networks [33][34][35], which demonstrates that ENs have a relatively small average shortest path but a very large cluster coefficient. Thus, the smallworld electrical networks reveal that if a node (or branch) in the network fails, the adjacent and even nonadjacent nodes (or branches) could fail, leading to cascading failures. Meanwhile, ENs has scale-free characteristics, as determined by analyzing changes in the network structure and function when one or more nodes (or branches) are removed from the network, which shows the networks are highly vulnerable under deliberate attacks but robust under random attacks [28]. However, fault propagation mechanism of ENs studied from the perspective of pure topological structure is not comprehensive and should more focus on the physical and operational features.
In communication networks, generally, there are two types of topological networks: scale-free networks and smallworld networks [15,36,37]. Communication networks with scale-free structures contain a few nodes with high degree, and they can be considered center nodes. Compared to the small-world networks, the distributions of degree of which are more uniform, scale-free networks have higher communication efficiency but are more vulnerable to deliberate attacks.
From the perspective of pure topological structures, compared to ENs, fault (or virus) propagation in communication networks is largely determined by its topological structure. That is, a fault node (or branch) only causes neighboring nodes to fail. Therefore, we employ the CNT to develop the virus propagation models (VPMs) based on the topological structures.
In summary, the EN and communication network in smart grids have two essential differences in terms of the interactions of cascading failures.
Features . From the perspective of time scales, the power flow is transmitted based on continuous time, while the communication flow is transmitted based on discrete time.
Features . From the perspective of topological structures, fault propagation in the communication networks depends more on the network structures. Compared with the communication networks, fault propagation in ENs depends more on physical and operational modes because ENs comply with operational rules, for example, Kirchhoff 's law.

Virus Propagation and Information Exchange Models in Communication Networks
Before analyzing VPMs, we introduce the following topological concepts: The degree of is the number of neighboring vertexes connected to , as expressed by Degree distribution p(k) is the distribution function of the degrees. That is, when a vertex is randomly selected from the network, the probability that its degree is equal to k is p(k).

. . Virus Propagation Model Based on CNT. According to
Feature , virus propagation in the communication network depends on the network structure; therefore, we employ the SI [38][39][40] and SIR [41,42] models based on CNT to simulate virus propagation. In the SI model, the vertexes of the communication network are divided into two groups: susceptible set S and infectious set I. In S, the probability that a susceptible vertex contracts the virus from the infectious vertexes is . Meanwhile, because the virus spends some time in destroying the functions of susceptible vertexes (e.g., tampering with data or instructions), the susceptible vertexes take time to get infected. Therefore, we introduce time delay (virus infection cycle) to develop the SI model as follows: On the basis of the SI model, the SIR model considers that a few infectious vertexes can be cured with the probability owing to the strengthening of related antivirus measures (e.g., formatting operation). Thus, the infectious vertexes can obtain immunity in a certain virus removal cycle. Therefore, the vertexes add a group called removed set R. The SIR model with time delay is given as Because an infectious vertex only transmits the virus to its neighboring vertexes, the topological structure of the network greatly influences virus propagation. When the communication network is a small-world network, which can be regarded as a uniform network owing to the relatively uniform distribution of degree [37], the degree of is approximately equal to ⟨ ⟩, and the SIR model can be presented as follows: When the communication network is a scale-free network, the vertexes have different damage levels from the perspective of virus propagation because the distribution of degree follows a power law. That is, the greater the of , the more serious it is for the to spread or contract the virus to more vertexes. The SIR model can be expressed as (7) in terms of the vertex degree [43,44].
Generally, in the SI and SIR models, virus propagation is faster in scale-free networks owing to the power law distribution. In addition, by comparing SI and SIR models, once virus propagation occurs in the communication network, we can investigate whether the related antivirus measures with time delay can play an important role to prevent the fault from spreading across the EN.

. . Information Exchange Model in Communication Network.
In the smart grid, the communication vertexes send operational data (parameters) associated with branches to the control center and receive commands from the control center step-by-step through the communication network in the form of information packets. At every step, the same information packets can be received and sent only by each communication vertex. Before constructing the information exchange model, three simplifications are made as follows.
(1) The communication blocks of vertexes (or edges) are not considered in process of the information transfer when the vertexes work orderly. That is, the capacity of the vertexes is adequate to exchange/handle the information packets.
(2) Because the vertexes send and receive information packets at intervals of 0.833 ms [45], the time required for information exchange between the vertexes and control center can be ignored because it is very small compared to the time required for fault propagation in the EN.
(3) Because we focus on the interactions between the ECNs, the methods of gathering and dealing with the information (such as the measuring units, transmission channels and protocols, encryption and decryption algorithms, etc.) are not considered.
Based on the above simplifications, the information exchange model is constructed based on the structure of the communication network.

Communication Rules between Target Vertex and Control
Center. The vertexes abide by the rule of first-in-first-out to send out information packets to avoid exchange of the information packets to be in the same edge. At first, the target vertex produces information packets. Then, the information packets are sent to all its neighbor vertexes. If the control center is one of the neighbor vertexes, the information transfer ends; otherwise, all neighbor vertexes, acting as target vertexes, continue to send the information packets to their corresponding neighbor vertexes until the information packets are sent to the control center. In the above process, the information packets are transmitted successfully between the target vertex and control center if a path exists between them in the communication network.

Dynamic Power Flow and Overload Mechanism Models in Electrical Networks
. . Dynamic Power Flow in Electrical Networks. To redistribute power flow during disturbances, we employ primary frequency regulation [46,47] and rotor equation [48] to model the dynamic power flow method.
System Frequency Characteristics. We employ primary frequency regulation to adjust the power flow. The characteristics of load and generation frequency are given by (8) and (9), respectively.
and are calculated as follows: Unbalanced Power Redistribution. To redistribute the unbalanced power P un due to disturbances of the system, we first calculate the change in system frequency by using the primary frequency regulation and the rotor equation.
In (12), P un is calculated as follows: When P un < 0, we consider the characteristics of load and the generation frequency to adjust the system frequency: When P un > 0, we consider only the generation frequency characteristic: By using (8)- (15), the changes in every generator and load can be expressed as follows: Then, we employ the P-Q power flow to calculate the power flows of each bus ( = = ) as follows: Dynamic Power Flow Method. When branches fault during fault propagation in the network, the dynamic power flow can be calculated in Algorithm 1.
. . Overload Mechanism of Electrical Networks. In this paper, cascading failures in the EN are analyzed from the perspective of the overload mechanism. When one or more lines are cut off, the other lines are overloaded owing to the redistribution of power flow in the EN [27,28]. When a branch is overloaded, the larger the power flow over the branch, the shorter is the operational time for which the branch is permitted to continue working [15]. As most of other studies [15,49], in this paper, we assume that during the fault propagation, the control center tries to maintain the secure operation of the EN and lower the load shedding amount, thus the control strategy includes which branches to trip, how to adjust the generators output, as well as to shed which load of how many percentages, etc. Therefore, for some of the branches, the tripping command has to come from the control center. Of course, for some of the faulted lines, the tripping signal should be issued by a local protection unit. However, to simplify the process, we simply assume that the tripping command comes from the control center. Thus, under thus simplification, when the corresponding communication vertexes send information about overloading to the control center via the communication network, the control center must quickly send trip commands to the target vertexes. We employ the inverse-time overcorrect protection scheme [15,49] to calculate the overloaded operational time.
If the data exchange between the target vertex and control center to trip the branch B j is completed within t Bj , the control is successful; otherwise, the control is unsuccessful.

Interactive Cascading Failure Model
A diagram of cascading failures considering the interactions between the ECN is shown in Figure 3. During the cascading Complexity 7 Input: Electrical network information Output: Power flows of branches, frequency violation Step 1: Unbalanced power: Employ Equation (13) to calculate the unbalanced power of the network.
Step 2: System frequency characteristic: IF < 0, employ Equation (14) to calculate and go to the next step; ELSE, IF > 0, employ Equation (15) Meanwhile, the virus will cause branches to trip or lead to an outage directly or indirectly because the infectious vertexes lose the function of information exchange. Accordingly, there are four types of fault branches.
Type . The branch is forced to trip because the corresponding communications get infected, also called forced outage branches B .
Type . The branch is tripped properly because the control center successfully sends commands to the corresponding vertexes based on the received overload information and control strategy with , also called overload tripping branches B .
Type . The branch is damaged irreparably owing to control failures via the communication networks, leading to overload operational time exceeding t Bj , also called irreparable fault branches B .
Type . The branch undergoes forced outage owing to network splitting, also called network splitting branches B . Based on the above analysis, cascading failure considering the interactions between the ECNs is modeled in Algorithm 2.

Case Study
The proposed model was applied to the IEEE 118-bus system and the French grid [50]. The small-world and scale-free networks were adopted to represent the respective communication networks. The computational work was performed in MATLAB running on a laptop. The laptop (Compaq, v3646TU) was equipped with an Intel5 Core6 2 Duo CPU T7250@2.00 GHz, 2.00 GB RAM, and 64-bit Windows 7 operating system. Step 6: Network operational status: Calculate the power flow over (j = 1, 2, . . . , M B ). IF is overloaded, calculate t Bj of by Equation (19), and add to B . Communication network: Step 8: Infectious vertexes detection: IF the infection time of the candidate vertex (k = 1, 2, . . . , M I ) in I is equal to t, add to I and delete it from I .
Step 10: Removal of vertex detection: IF the removal time of the candidate vertex (w = 1, 2, . . . , M R ) in R is equal to t, add to R, and delete it from R . Step 7: Virus propagation: The susceptible vertex ( = 1, 2, . . . , M S ) contracts the virus with probability according to Equation (6). IF gets infected, delete from S, add to I , and label its infectious time t + 1 .
Step 9: Vertex immunization: The infected vertex (l = 1, 2, . . . , M I ) is immunized with the probability according to Equation (6). IF obtains immunity, delete from I, add to R , and label its immunity time t + 2 . Step11: t = t + Δt; END WHILE.  Table 1. Figure 4 shows that the load shedding changes with the passage of time based on the different topological structures of the communication networks and the SI model. Owing to space limits, the SIR-model-based load shedding is not given herein. In Figure 4 Table 1. Therefore, the coupling of EN with the scale-free communication network is affected more severely as the propagation time increases because the connectivity between vertexes often depends on a few hub vertexes (i.e., high-degree vertexes), and the exchange of information packets becomes difficult, leading to network paralysis once a few hub vertexes are infected.
A comparison of the SI and the SIR models shows that different VPMs have very small impacts on fault propagation  However, the number of fault branches with irreparable faults can be reduced, especially in the case of coupling with the scale-free communication network. This is because the immune vertexes treated with the antivirus recover their function of data exchange, and a few overloaded branches can receive trip commands from the control center in a timely manner, thus avoiding irreparable faults.
Furthermore, we analyze the interactive cascading failures by selecting different initial infectious vertexes. We used the SI model as an example. Because the degree distributions of the small-world communication network are known, we take the scale-free network as the basis to select the highdegree (vertexes 115 and 116) and small-degree (vertexes 4 and 8) vertexes as the initial infectious vertexes. Figures 5  and 6 show the total and real-time load shedding changes with the passage of time for different virus propagation times 1 = 2 s, 5 s, and 8 s. The initial vertexes have small impacts on fault propagation in the EN owing to the known degree distribution. However, in case of the coupling of the EN with the scale-free communication network, the initial vertexes greatly impact fault propagation. Compared to the smallworld communication network, when the initial vertexes are high-degree vertexes in the scale-free communication network, the propagation time is obviously shorter, which demonstrates the hub can rapidly spread the virus, leading to rapid collapse of the EN. By contrast, the propagation time is longer when the initial vertexes are the low-degree vertexes, and when the load shedding peaks, as shown in Figure 6, the interactive cascading failures continue to spread, which indicates virus propagation times are longer than fault propagation times. That is, when fault propagation has stopped, virus propagation continues.
. . French Grid. A real French grid with 1951 nodes and 2956 branches was employed to simulate the interactive model. Owing to computational complexity, we only choose the high-degree vertexes as initial infectious vertexes considering the topological structures of the communication network with the parameters Δt = 0.01 s, 1 = 2 = 2s, = = 0.3, = 7, and = 1.5. Figure 7 shows that the total and real-time load shedding of system changes with the passage of time. Compared to the small-world communication network, the propagation times are longer in the case of EN coupled with scale-free communication network, but the load shedding peaks at approximately 8.3 s under both topological structures.
Furthermore, we investigated the numbers of the three types of fault branches (M IF , M OT , and M FO ) at different moments, as shown in Figures 8(a) and 8(b). Between 6 s and 7 s, the fault propagation is at its height, which indicates that the numbers of infectious vertexes and forced outage branches B are the highest, leading to rapid collapse of the EN and surging load loss. In addition, in the fault propagation process in the EN, the fault branches with irreparable faults are not found when the EN is coupled with the small-world communication network. By contrast, there are many fault branches of this type at different moments in the case of EN coupled with the scale-free communication network. Therefore, when the communication network is scale-free, it is more vulnerable which cannot effectively resist the virus propagation leading to the more severe damage to the EN.
The conclusions obtained from these two cases are summed up in Table 2. In practice, when ENs are faced with hacker attacks, because the attackers find it relatively difficult to obtain complete information about the communication networks, such as topological structures, their attacks are random to some extent, which means scale-free communication networks are more appropriate for the ENs [10,30,31]. However, when ENs are faced with the threat of a cyber virus, the virus must be cleared promptly. Once the virus spreads, regardless of whether the initial infectious vertexes are selected randomly or deliberately, the infection will lead to a severe damage in the scale-free communication networks. Therefore, due to the propagation features of the scalefree communication networks, the simulation results show that software engineers should strengthen more the software

Discussion
In this paper, we extend the state of the art for the study of the integrated communication network and electrical network. Compared with other literature, the main contributions of our paper are as follows: We propose an interdependent fault propagation model which holistically considers the extreme virus propagation in the communication network to reveal the vulnerability of electrical network coupled with different communication network structures at the first time.
In the fault propagation model, to better reproduce the ex-post behavior of the electrical networks, we extended the dynamic power flow by including the primary frequency regulation and the equations of rotors of generators.
To solve the issue of different time frames in the interdependent system, we adopt two time windows, i.e., the virus infection cycle of nodes and tripping time of overloaded branches during fault propagation to analyze the fault mechanism of both electrical branches and communication nodes along time.
It should be noted that even though the electrical network and communication network are both presented as graphs to conveniently describe their interdependent topological relationship in this paper, the modeling approach captured most of the relative features of the two networks.
For the electrical network, besides the commonly considered steady state physical and operational rules, we also adopt the rotor equation and system frequency to consider simple system dynamics in order to present the interactions between   the electrical network and the communication network. As for the communication networks, we assume that the security level of each node is the same in terms of the infected rate. By contrast, in reality, the probability of the communication nodes got infected may vary for different nodes. However, the assumption made in our paper does not change the essence of the analysis and results in terms of evaluating which topology of communication network would have higher impacts on the electrical network during cyber-attacks.

Conclusions
The cyber-physical security of power systems is attracting increasing attention, especially after more and more evidences show that failures or attacks happening in the cyber system can greatly destroy the secure operation of power systems and bring tremendous consequences. To investigate the possible consequences, we propose an approximate interactive model to study cascading failures in ENs caused by virus in communication networks via two types of propagation. Our simulation on a standard study case, i.e., IEEE 118-bus system, and a realistic network, i.e., French grid, shows that the structure of the communication network has decisive impacts on the ECN in terms of the propagation time of cascading failures, loading shedding, number of faulted branches, etc. However, due to the simplification of the communication network and the virus propagation mechanism, the model can still be refined. In addition, the analysis is only focused on the overload of the system which may limit the results to part of behaviors of the EN.
Owing to the complexity of the propagation mechanism of interactive cascading failures, future work in this field will focus on considering more factors, such as data transmission delay, to simulate interactive cascading failures. Meanwhile, we also can investigate the impacts of differences of virus infection of nodes on interdependent fault propagation for electrical and communication networks. In addition, we can also analyze other aspects of the integrated CPS system, such as reliability, resilience, etc., under the virus propagation, to provide other dimensions for understanding the CPS. Constants : The branch between N i and N r V : The edge between V m and V v P Wrc : Power rating of generator c P Rrq : Power rating of load q K Wc : cth generator unit power regulation K W : Equivalent generator unit power regulation K Xq : qth load frequency regulation K X : Equivalent load frequency regulation K E : System unit power regulation T J : Equivalent Inertia time constant T Jc : Inertia time constant of generator c k Vm : Degree of vertex m p(k): Degree distribution : Proportional coefficient of inverse-time overcurrent protection : Power coefficient of inverse-time overcurrent protection : Current limit of branch B j .
Variables t: Time / clock t Bj : Overloaded operational time of branch B j : Simulation time ⋅(t): Value / set of a variable at time t Δt: Time step ⋅ k : Value / set of vertex(es) with k degrees P un : System unbalanced power Δf: Frequency offset ΔP W : Changes of power of all generators ΔP Wc : Changes of power of generator c ΔP X : Changes of power of all loads ΔP Xq : Changes of power of load q P Xq : Power of load q I Bj : CurrentoverbranchB j P i : Injection active power of node i Q i : Injection reactive power of node i B iu : Equivalent susceptance between nodes i and u G iu : Equivalent conductance between nodes i and u : Voltage phase angle difference between nodes i and u ] : Voltageofnodei Δ : Changes of angular acceleration of equivalent generator : The

Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
All authors declare that they have no conflicts of interest.