Modeling of Time-Delayed Distributed Cyber-Physical Power Systems for Small-Signal Stability Analysis

Compared with a power system under centralized control, a power system under fully distributed control has stronger interdependence between its physical side and cyber side. Such a system is a typical cyber-physical system (CPS). The time delay problem is one of the most critical issues in the application of distributed control for power systems. First, a distributed information flow model is established considering multiple sources of time delay, especially the communication delay between neighboring agents. Second, combining the dynamics of physical power systems, a general model of time-delayed CPS under fully distributed control is proposed, especially for inverter-based microgrids. Third, the stability of the time-delayed CPS is analyzed based on the proposed model with a critical eigenvalue tracking algorithm. A real-world microgrid with DGs under fully distributed frequency control is tested for small-signal stability analysis. The findings show that the stability of the distributed CPS with multiple sources of time delay depends on not only the absolute value of time delays but also the coordination of time delays and that the delay merging property in centralized mode is not relevant in a distributed CPS. The accuracy of the model and the obtained stability margin are verified by Simulink simulations.


I. INTRODUCTION
T HE ASTOUNDING pace of the construction of smart grids has transformed the traditional power system into a typical cyber-physical system (CPS). With the application of advanced and high-performance information and communication technology (ICT) in modern power systems, the operational efficiency of the cyber-physical power system (CPPS) has undergone a qualitative leap [1]. In recent years, due to ubiquitous measurement units, control systems and advanced information systems such as the energy management system (EMS) and the wide area measurement system (WAMS), the awareness of modern power systems has become increasingly comprehensive, and the control philosophy is more efficient than ever.
Moreover, the safe and reliable operation of power systems is highly dependent on the cyber side of such systems. Disturbances or contingencies on the cyber side of an information system have significant impacts on the physical side. In particular, time delays are among the most critical issues in information systems. The time delay problem remains one of the main obstacles to the revolution of WAMSs and the industrial implementation of wide area control systems (WACSs) [2]. The time delay of feedback signals for a power system stabilizer (PSS) may threaten the stability of the low-frequency oscillation of a large power system. Field tests performed to measure the time delay of feedback channels in WAMSs in a real system, the Guizhou Power Grid in China, showed that time delays might change remarkably when routing changes or the communication load increases [3]. Therefore, the impact of a time delay on the stability of the most complex system, a CPPS, is a matter of great concern.
With the development of renewable energy, traditional single-entity power systems have tended to become multiagent and diversified CPPSs [4]. The distributed control mode has received extensive attention in the field of power systems due to its plug-and-play, scalability, fast response and privacy protection features, which distinguish it from the centralized control mode and the decentralized control mode [5]. Moreover, with the increasing penetration of renewable generation into power systems, the challenges of power system operation and regulation have gradually become increasingly prominent. Given the lack of frequency and voltage regulation capabilities of power systems with high-penetration renewable energy, especially low-inertia microgrids [6], conventional operation and control methods are unable to cope with the active control problem in high-penetration renewable generation [7], [8]. Compared with traditional centralized control, distributed control provides a more flexible control mode for a power system while making the power system highly dependent on the reliability and immediacy of communication on the cyber side. A pioneering distributed control architecture for microgirds is practically implemented in a utility-scale project, IIT-Bronzeville Microgrid Cluster in [9]. Thus, due to the strong coupling and interdependence of the cyber side and physical side, a power system under distributed control is a distributed cyber-physical power system (DCPPS). The time delay problem in a DCPPS is especially critical and is nonnegligible under distributed control [10].
The existing studies of stability analysis of time-delayed systems [11]- [17] mainly focus on power systems in centralized mode, such as a power system with a PSS [12]. The stability of a time-delayed energy storage system with centralized and decentralized control levels is analyzed in [15]. A general model and formulation of a large delayed CPPS with PSS controllers is proposed in [16], which verifies delay merging in centralized mode.
However, the stability of time-delayed power systems under distributed control, which has higher time delay requirements and more complex information flows than centralized control or decentralized control, has not yet been studied, especially from a CPS perspective.
Although pioneering studies [11]- [16] have analyzed the stability of time-delayed power systems under centralized or local control, in regard to CPPSs under fully distributed control, some challenges remain. First, in a DCPPS with multiple sources of time delay, both the local measurement and control execution delay and the communication delay affecting neighboring agents need to be considered on the cyber side, which makes the distributed information flow complex. Second, unlike the existing studies that are concerned with the effect of time delays on the convergence problem of distributed algorithms [10], [18]- [20], the dynamic characteristics of both the physical side and cyber side should be considered in the stability analysis of a time-delayed DCPPS. Third, because of the coupled multiple sources of time delay between neighboring agents, the properties in centralized mode may no longer apply to a DCPPS. The impact of multisource delays on the stability margin of a DCPPS must be determined.
To address the above challenges, based on the proposed distributed information flow model, a general delay differentialalgebraic equation (DDAE) model for a DCPPS is proposed. The stability of a real-world islanded microgrid with distributed generators (DGs) under fully distributed frequency control is analyzed in this article based on the proposed model and critical eigenvalue tracking algorithm. To the best of the authors' knowledge, the contributions of this article are as follows: 1) A general cyber-physical model for CPPSs under fully distributed control considering multiple sources of time delay is proposed in this article. The proposed model is formulated as a set of DDAEs combining the distributed information flow on the cyber side and the power system on the physical side.
2) The DDAE-based DCPPS model is transformed into a set of delay differential equations (DDEs) by eliminating algebraic variables and merging cyber-physical state variables. Then, a critical eigenvalue tracking method is proposed to solve the transcendental characteristic equation of the DDEs with multiple periodic delay terms. Accordingly, the accurate delay margin of the DCPPS in the marginally stable state can be obtained, and the result is verified by Simulink simulations.
3) A stable region for small-signal stability analysis of a CPPS under fully distributed control considering multiple sources of time delay can be obtained based on the proposed model and method. The findings indicate that the stability of a CPPS under fully distributed control has different characteristics from those of a CPPS under centralized or local control. The sensitivities of the critical eigenvalue with respect to time delays and system parameters are also derived for system guidance and controller design.
The remainder of this article is organized as follows. Section II introduces the mathematical model of the distributed CPPS. Section III presents a stability analysis of the timedelayed CPPS under fully distributed control, which is verified by case studies in Section IV. Finally, Section V concludes the paper.

A. Modeling of Distributed Information Flow
A centralized information-energy flow model based on a one-direction closed-loop information flow model is introduced in our previous work [21], [22]. However, since each neighbor in distributed control needs local measurement information and must iteratively communicate with neighbors, the information flow of a DCPPS is more complicated than that of other systems.
An illuminating general structure of online algorithms for optimization problems is presented in [23], and the interactions between physical and cyber networks are clearly described. To study the impacts of multiple sources of time delay on the stability of DCPPSs, according to the distributed control framework of an agent in a DCPPS, as shown in Fig. 1, we consider a distributed information flow to involve three activities in the form of DDAEs: 1) measurement, 2) communication and iteration, and 3) execution.
According to Fig. 1, physical variables in a physical power system include physical state variables and physical algebraic variables, and these variables can be used jointly to determine the energy flow of the physical system. Cyber variables include cyber state variables in iterative calculation tasks and other algebraic variables. The measurement variables z c i of agent i from the physical side to the cyber side are considered cyber variables. Conversely, the control variables u p i from the cyber side to the physical side are defined as physical variables. 1) Measurement: In a measurement task of distributed information flow, agent i first measures the required local physical state variables through its measurement unit. Additionally, information is received from neighbors N i . Since the information from neighbors mainly affects iterative calculations, we consider only the local measurement variables in the measurement task.
Considering the time delay in the local measurement task, the mapping from local physical state variables on the physical side to measurement variables on the cyber side of agent i is established as follows: where superscript p indicates the variables on the physical side, superscript c indicates the variables on the cyber side, x p i ∈ R n i ×1 and y p i ∈ R l i ×1 denote the vectors of the physical state variables and physical algebraic variables of agent i, respectively. We assume that for each agent n i = n, for i = 1, 2, . . . , N, and l i = l, for i = 1, 2, . . . , N, where N denotes the number of agents, n and l respectively refer to the dimension of the vector of the physical state variable and physical algebraic variable corresponding to an agent. z c i ∈ R (n+l)×1 is the vector of the cyber measurement variables. τ m,i refers to the local measurement time delay of agent i. M x,i ∈ R n×n and M y,i ∈ R l×l are the physical state variable measurement matrix and the physical algebraic variable measurement matrix, which are diagonal matrices that connect the physical side and the cyber side of agent i. When there are no measurement errors, diag(M x,i , M y,i ) = I.
2) Communication and Iteration: The communication and iteration task involves the specific implementation of a distributed algorithm. Each agent uses the information from its neighbors and the local measurement variables to perform iterative calculations.
A distributed communication network topology is expressed as G = (v, e, A), where v = {v 1 , . . . , v n c } is a set of nodes in the communication network, n c denotes the number of agents, e ⊆ v × v is a set of edges, and A = [A ij ] N×N is the weighted adjacency matrix of the communication network.
All neighboring nodes directly connected to node v i through communication links can be represented by set N i (including node v i itself) as follows: where j ∈ N i refers to a set of nodes adjacent to node v i . Considering the dynamic characteristics of the fully distributed algorithm, the iterative calculation of the fully distributed controller of agent i is expressed as the following continuous model (CT) in the form of a differential equation: where x c i ∈ R c i ×1 is the vector of the cyber state variables of agent i, x c j∈N i denotes the vector of the cyber state variables of the neighbors of agent i, and h i (·) is the adjustment function of the distributed algorithm of agent i based on the information obtained from node i and its neighbors. Since each agent uses the same distributed control strategy, we assume that the dimension of the vector of cyber variables is c i = c for i = 1, 2, . . . , N. Thus, we have x c i ∈ R c×1 . To consider the time delay, the communication and iteration task is represented as the following continuous model by combining equations (1) and (3): where τ m,i denotes the local measurement time delay, τ c,ij represents the communication time delay of link (i, j) ∈ e, and τ c,ii = 0 for i = 1, 2, . . . , N.
3) Execution: The execution task involves the interface between the control signal on the cyber side and the control variables on the physical side. The control delay affects the stability of the distributed system, which should be considered. Thus, the algebraic equation for the execution task is represented as follows: where τ e,i is the execution time delay, u p i ∈ R c×1 is the vector of the physical control variables of agent i, and C i ∈ R c×c is a diagonal matrix similar to M x,i and M y,i ; this matrix can be used for sensitivity analysis of control variables.

B. Modeling of a Dynamic Power System
To analyze the stability of a distributed CPPS, the electromechanical model of the power system on the physical side is represented as a set of differential-algebraic equations.
1) Synchronous Generator Model: According to [24], we introduce the following second-order swing equation model to present the differential equations of generator i in a dynamic power system: where δ i is the rotor angular position of generator i in rad, ω i denotes the angular speed of generator i, ω N denotes the rated angular speed of the system, M i is the inertia constant, and D i refers to the damping coefficient. P m,i and P e,i represent the shaft power input and electrical power output of generator i, respectively.

2) DG Model of an Inverter-Based Microgrid:
The dynamics of DGs under the P-ω control strategy in an inverter-based islanded microgrid can be modeled as the following set of differential-algebraic equations [25]: where ω * is the angular speed of the generator at a phase angle reference bus in the islanded microgrid; F pi denotes the cutoff frequency of the first-order low-pass filter used in power measurements for DG i; K pi is the frequency droop gain; and P m G,i , P G,i and P * G,i denote the measured, real and rated active power generation of DG i, respectively.
Then, by eliminating P m G,i in (7), the differential equations of the DGs under the P-ω control strategy can be expressed as follows: where P G,i = P inj i + P d,i , P d,i is the load demand at bus i and P inj i denotes the active power injection at bus i, which can be calculated from the power flow equations. We note that this article chooses a distributed frequency control framework as a typical example for the test, so Q-V characteristics are not considered; this approach does not affect the derivation given later in this article.

3) General Model of a Dynamic Power System:
The differential equation model of a power system is mainly related to the rotor motion of synchronous generators, which can be expressed in the following abstract form: where x p ∈ R N·n×1 , y p ∈ R N·l×1 and u p ∈ R N·c×1 denote the vectors of the physical state variables, physical algebraic variables and physical control variables of the power system, respectively, and f p (·) is the abstract form of the swing equation of the power system. In (9) where g p (·) is the abstract form of the power flow equation. The original DDAEs need to be converted into state space form based on the DDEs in small-signal stability analysis, and the corresponding delayed algebraic variables need to be eliminated. The power flow equation is memoryless and fully determined by the state and algebraic variables at the current time, as shown in Fig. 2. To combine the cyber side and the physical side of a distributed CPPS, it is necessary to establish corresponding DDAEs according to different delay conditions. Considering the time delay in the measurement and execution tasks, the DDAE-based model is defined as follows:

C. Modeling of a Distributed CPPS
Considering the multiple sources of time delay for measurement, neighboring communication and execution tasks, the above distributed CPPS model is rewritten in the following DDAE form by combining the distributed information flow model on the cyber side and the dynamic power system model on the physical side: Notably, unlike the model for the centralized control mode in [16], the power system in distributed control mode needs to consider not only the local measurement delay and execution delay but also the communication delays between different agents; thus, we have established the abovementioned distributed information flow model.
As stated in the general model of a dynamic power system, the power flow equation is memoryless. Thus, the above model is based on the index-1 Hessenberg DDAEs; the corresponding algebraic equations do not simultaneously depend on ordinary variables and time-delayed variables and can thus be further transformed into DDEs.

A. Characteristic Equation
In this section, we focus on stability analysis of a timedelayed DCPPS. Compared with the traditional time-delayed power system model, the cyber-side model of a distributed CPPS is constructed with additional specifications, especially in relation to the distributed information flow model for distributed control. The distributed CPPS model provides the foundation for stability analysis of time-delayed DCPPSs.
Small-signal stability analysis of time-delayed DCPPSs at the equilibrium point can give the necessary conditions for assessing the stability of time-delayed DCPPSs.
At the equilibrium point, linearizing the DDAEs of the distributed CPPS leads to where Since the algebraic equations considering different time delays on the physical side do not change the basic function, the Jacobian matrices for these equations are the same. Thus, By combining the state variables of the cyber side and the physical side, we can reconstruct a vector of cyber-physical state variables: Therefore, the separate differential-algebraic equations on the physical side (13) and the cyber side (14) can be rewritten based on the reconstructed cyber-physical state variables as follows: where A (ij) ∈ R N·c×N·c is the augmented matrix of A c ij ∈ R c×c , which satisfies In this case, A c ij is an element (i, j) of A (ij) considering the time delay between the communication with agent i and agent j.
i and C i , respectively. Note: Unlike the formulation of a CPPS in centralized mode, a communication delay in a system under fully distributed control is global and affects neighboring agents, so the corresponding matrix needs to be augmented at the whole system level.
Transforming the linearized DDAE model (16) into DDEs yields For explicit expression, the above DDEs are written in the following form: where Then, the characteristic equation of the distributed CPPS with time delays is expressed as follows: It is noteworthy that the above characteristic equation is a transcendental characteristic equation containing multiple sources of time delay (i.e., multiple time delay terms: e −sτ m,i , e −sτ c,ij and e −sτ e,i , which represent the measurement, communication and execution time delays, respectively). Since these terms are periodic, the transcendental characteristic equation cannot be solved directly.
In the next subsection, we focus on the algorithm used to solve the transcendental characteristic equation.

B. Critical Eigenvalue Tracking Algorithm
As mentioned above, unlike the characteristic equation of ordinary differential equations (ODEs), the characteristic equation of the DDEs, such as equation (21), is a transcendental equation with periodic exponential terms; such an equation has infinitely many solutions and cannot be solved directly. In addition, unlike the characteristic equation for a power system with a single-time delay local controller, which can be directly solved by the function transformation method (e.g., Padé approximants in [26]), the small-signal stability of a DCPPS is affected by multiple sources of time delay, which make the transcendental equation comparatively complicated.
Accurately determining the delay margin of a time-delayed DCPPS is the major concern of this article. For this purpose, it is essential to use the exact equivalent model of a time-delayed DCPPS in the marginally stable state. When the rightmost eigenvalue λ r reaches the imaginary axis as the delay increases, the system is in the marginally stable state, and the real part of the rightmost eigenvalue Re(λ r ) = 0 [15]. Notably, when Re(λ r ) = 0 is satisfied, s = σ + jw in the characteristic equation can be replaced by jw d in the marginally stable state, and the characteristic equation in (21) can be rewritten as follows: where τ mar m,i , τ mar e,i and τ mar c,ij denote the delay margins of the measurement delay of agent i, control delay of agent i and communication delay of link (i, j) ∈ e, respectively.
However, due to the existence of multiple unknown delay margins, the above characteristic equations cannot be solved directly. Thus, we define the following delay direction vector d to represent the direction of the multiple sources of time delay at a certain ratio: In the direction based on multiple sources of time delay for a certain ratio, the characteristic equation in (21) can be replaced by: where θ ∈ R is a scalar that reflects the magnitude of the delay direction vectord. The periodicity of the characteristic equation depends on the least-common multiple of θd m , θd e and θd c . It can be observed from equation (24) that the delay margin depends on the magnitude of the delay θ , which reflects the absolute value of the synthetic time delay in a certain direction, and the coefficient of the direction vector, which reflects the relative values of the time delays. Based on the above characteristic equation, a critical eigenvalue tracking algorithm is proposed, as illustrated in Algorithm 1.
Remark 1: The purpose of Algorithm 1 is to solve θ = θ d when the real part of the rightmost eigenvalue Re(λ r ) = 0, which indicates that the system is in a marginally stable state. Therefore, the delay margin can be accurately obtained.
The proposed algorithm considers the relationship between multiple sources of time delay in a DCPPS. θ reduces the dimensionality of the multiple sources of time delay. The transcendental delay terms in (21) are replaced by constant terms at a certain θ . Moreover, by setting the delay interval θ according to the least-common multiple of different delay ratios, the first marginally stable state caused by the increase in delay can be quickly located. According to Algorithm 1, when the absolute value of Re(λ r ) is less than the set threshold ε, the algorithm converges. Then, the corresponding delay margin vector can be obtained by Algorithm 1 Critical Eigenvalue Tracking Algorithm 1: Initialization: Determine the direction vector of multiple sources of time delayd. Define θ min = 0, θ max = min 2π/d ,θ = θ min , and θ = θ min . Set the convergence tolerance ε ≥ 0 and k = 1. Set the number of intervals Nwith interval θ = (θ max − θ min )/N 2: while k < N do 3: Setτ m,i =θd m,i , τ c,i =θd c,i , τ e,i =θd e,i 4: Solve the characteristic equation (21) at the currentθ , obtain the rightmost real part of eigenvalueλ r 5: if Re{λ r }≺ 0 then 6:θ ←θ + θ , k ← k + 1 7: else 8: θ max ←θ, break // the system is unstable at the currentθ 9: end if 10: end while 11: Loop 12: Set τ m,i = θd m,i , τ c,i = θd c,i , τ e,i = θd e,i 13: At the current θ , solve the characteristic equation (21) and obtain the rightmost real part of eigenvalue λ r 14: if Re{λ r }≺ ±ε then // the stability margin is sufficiently accurate 15: Based on the proposed algorithm, the delay margin of a DCPPS under different ratios of multiple sources of time delay can be efficiently obtained.

C. Sensitivity Analysis
To study the stability of a DCPPS, it is critical to derive the corresponding sensitivity, which provides the foundation for the optimal design of the controllers of the DCPPS.
First, to analyze the impacts of multiple sources of time delay on the delay margin of a time-delayed DCPPS, the sensitivity of eigenvalue λ with respect to the measurement time delay τ m,k , k = 1, . . . , N c is formulated as follows: The sensitivity of eigenvalue λ with respect to the communication time delay τ kl , (k, l) ∈ e and the execution time delay τ e,j , j = 1, . . . , N c are expressed as follows: and ∂λ ∂τ e,k = −λe −λτ e,k u H T e,k v u H (I + )v Remark 2: If the real part of the sensitivity of the rightmost eigenvalue λ r with respect to a time delay τ has Re(∂λ r /∂τ ) < 0, it means that a slight increase in τ at the current state is beneficial to the stability of a time-delayed DCPPS. Otherwise, such an increase is not conducive to the stability performance. A similar relationship also applies to the sensitivities of λ with respect to parameters.
From the above sensitivity expressions (26)- (29), it can be seen that the stability performance of a time-delayed DCPPS depends on the corresponding time delay, the characteristic root, the characteristic vector of the DCPPS and the coefficient matrices of the corresponding time delay and other time delays in formula (27).
Then, to analyze the impact of the parameters on the eigenvalues, especially the critical rightmost eigenvalue, the sensitivity of eigenvalue λ with respect to parameter p is expressed as follows: The derivation of the sensitivity is provided in Appendix A. Similarly, it can be observed from the above sensitivity expression (30) that the sensitivity highly depends on the coefficient matrix corresponding to the parameter in the characteristic equation.

IV. CASE STUDY
In this section, a 10.5 kV AC islanded microgrid system from a real-world 3-feeder 62-bus distribution network in Xiamen, China, is selected as a DCPPS and tested based on a fully distributed power dispatch method for frequency recovery.

A. System Under Fully Distributed Frequency Control
Fully distributed frequency control can provide a fast dynamic response and rapid frequency recovery. A fully distributed frequency control method that includes a subgradient term and consensus term can even minimize the generation cost by achieving the equal incremental principle during frequency recovery [27]. The detailed cyber-side model of the CPPS under fully distributed frequency control is explained in Appendix B.

B. Test Introduction
The 10.5 kV islanded microgrid with 4 DGs, including 2 synchronous generators (SGs), a virtual synchronous generator (VSG) based on a storage system and a PV under P-ω control, is shown in Fig. 3. The PV under P-ω control is based on the inverter-based microgrid model; F p is set to 60r ad/s, and K p is set to 0.1, in accordance with [25]. The phase angle at bus 2, δ 2 , is set as 0 • . The black solid line is the transmission line on the physical side, and the gray dashed lines denote the communication links between agents. The tolerance ε in Algorithm 1 is set as 10 −12 , which makes the obtained  stability margin very accurate. Because the communication networks in power systems are mostly duplex communication systems, for the convenience of analysis, we assume that the communication delays of the same link are equal, that is, τ c,ij = τ c,ji , (i, j) ∈ e. More test data can be found in [28].

C. Stability Analysis of the Time-Delayed System
Based on the proposed stability analysis method, the timedelayed margin of the tested islanded microgrid can be calculated.
1) One-Dimensional Time Delay: First, we reduce the multiple sources of time delay to a one-dimensional time delay τ ; that is, we assume that all delays in the system are equal and that the system only has a stability margin τ mar .
According to the proposed model and the stability analysis method, the stability margin of the tested time-delayed system is To verify the accuracy of the calculated stability margin, a DDAE model of the DCPPS is developed in MATLAB/Simulink. Fig. 4 presents the simulation results for the rotor angular speed of generator 1 under a small disturbance in three different time-delayed scenarios. The simulation of the system starts at the equilibrium point. The small disturbance is set as a 1% increase in the load at bus 1 and is added at 10 s in the simulation. Fig. 4(a) shows that when the time delay τ = 0.99τ mar is slightly shorter than the theoretical stability margin τ mar , the rotor angular velocity of generator 1 will gradually return to the stable equilibrium point with an increasing number of iterations. When the time delay is equal to the theoretical calculation value, that is, τ = 165.2ms, the system is in a marginally stable state, which means that the real part of the rightmost eigenvalue of the distributed CPPS is on the imaginary axis in this scenario, as shown in Fig. 4(b). When the time delay τ = 1.01τ mar is slightly longer than the calculated margin, the system is unstable. The accuracy of the stability margin calculated by the proposed method is verified by simulations.
We compared the calculation results of our method with those of the Padé approximant method [26] in the frequency domain and the linear matrix inequality (LMI) method based on the Lyapunov functional in the time domain [29]. The simulations were performed on a laptop with an Intel Core i5-6200U CPU and 16 GB RAM. The LMI was constructed as a semidefinite programming (SDP) problem and solved using MOSEK [30] (ver. 9.2) in MATLAB R2019b. The Padé approximant method and our method were coded and solved in MATLAB R2019b. The comparison results are shown in Table I. According to [26], the sixth-order Padé approximant is a common choice in numerical simulations. Table I shows that the delay margin calculated from the LMI approach in the form of an SDP problem based on the Lyapunov functional has a certain degree of conservatism (155 ms versus the accurate delay margin of 165.2 ms verified by Simulink simulations). Additionally, Table I shows that the algorithm presented in this article is more efficient than the other methods in the delay margin computation for the DCPPS in the marginally stable state.
2) Multiple Sources of Time Delay: To compare the relationship between different kinds of time delays, we assume that the same kinds of delays are equal for different agents. Thus, the measurement, communication and execution time delays in the tested system can be reduced to τ m , τ c and τ e , respectively.
First, to simply analyze the relationship between the local delay and communication delay with neighbors, we assume that the local delay involves only measurement delay. Fig. 5 illustrates the stable region considering the local measurement delay and communication delay calculated by  the proposed algorithm. The filled green region denotes the stable region, and the boundary is the stability margin in the marginally stable state. The stability of four representative points in Fig. 5 is analyzed based on the simulations in Simulink, as shown in Fig. 6.
According to Fig. 5, Fig. 6, and Fig. 7, the following key findings can be obtained for the CPPS under fully distributed control with time delays: 1) By comparing points A and B in Fig. 5, it can be observed that when measurement and communication delays exist, a shorter communication delay does not necessarily contribute to the stability of the system.
2) As shown by Area E in Fig. 5 and the simulation results for Area E in Fig. 7, when the measurement delay is less than 36.7 ms, the system is stable regardless of the communication delay. The result means that when the system meets the measurement time delay constraints, stability can be achieved, even if there is a very long communication delay among neighbors, such as 10 s.
3) Due to the nonlinear feature of the transcendental characteristic equation, when the root trajectory of the Fig. 7. Simulation of the rotor angular speed of SG1 at two points relative to Area E. dominant root is tangential to the imaginary axis, the dominant root will change at a certain point, which will cause the stable margin to jump at that point, as shown at points C and D.
Remark 3: It is noteworthy that since the consensus term in the fully distributed control algorithm (45) is mainly used for economy (to meet increment rate equality criteria) instead of stability, when the measurement delays in the subgradient term are less than a certain value, the stability performance of the time-delayed DCPPS can be guaranteed, as shown in Area E in Fig. 5.
Considering the tested fully distributed frequency control algorithm, the communication delay exists in the consensus term of the increment rate (ICR) between generators, which is mainly for economic optimization rather than improved stability. Therefore, it is reasonable to state that a delay on a certain small scale is beneficial to stability, as shown at points A and B in Fig. 5.
For a generalized analysis, multiple sources of time delay, including measurement delay τ m , communication delay τ c and execution delay τ e , are simultaneously considered in the stability analysis of the CPPS under fully distributed control. The direction vector of this three-dimensional time delay is defined in spherical coordinates: where θ is the polar angle and ϕ is the azimuthal angle in the spherical coordinate system.
The stability margins under various directions of multiple sources of time delay are presented in Fig. 8. The figure clearly shows that the stability margin of the system is different for different ratios of the multiple sources of time delay. Unlike the delay merging property of centralized control architecture [16], in this case, the sums of all three types of delay margins are not equal. Thus, all delays for each agent under distributed control cannot be merged. The primary reason for the noncumulative property of the DCPPS is that interactions with other neighbors are involved in the iterative calculation.
For some specific direction vectors such as (60 • , 60 • ) and (70 • , 70 • ), the stability performance of the tested DCPPS is significantly worse than that for other direction vectors. Notably, the system has a longer delay tolerance for all three types of delay and for multiple sources of time delay in the direction (45 • , 50 • ) than in the direction (40 • , 40 • ). 3

) Heterogeneous Communication Delays:
In this subsection, we consider that the delays of all communication links are not identical and are normally distributed with variance σ 2 = 0.1, where the means of the heterogeneous communication delays depend on the transmission rate and the length of the corresponding communication links in [28]. The stability margins of the tested DCPPS considering heterogeneous communication delays under different transmission rates of communication links are shown in Table II. The transmission rate is reflected by the ratio of communication delay τ c,12 to measurement delay τ m,1 . Here, we assume that the measurement delay and communication delay of the system are equal.
It can be seen from Table II that when the communication delay has the greatest proportion, the stability margin of the measurement delay is approximately 20 ms. Based on this relationship, we can determine that when the measurement delay is shorter than 19.3 ms, the system can still ensure stability despite large communication delays. From this table, it can be observed that the stability region of the system is very small when τ c,12 /τ m,1 = 0.399, which is similar to point D in Fig. 5. This effect is also due to the change in the dominant eigenvalue. This phenomenon indicates that the ratio of local delay to communication delay should be considered when designing a DCPPS.

D. Analysis of the Delay Merging Property
In centralized mode, the measurement delay and control execution delay in one control loop can be merged into a synthetic delay; in this case, the sum of the measurement delay and control delay is fixed at the stability margin, and τ mar i = τ mar m,i +τ mar c,i . It is worth exploring whether this property applies to a DCPPS.  We assume that the communication delays in the tested systems are the same. Table III shows the stability margins of DG 1 under different delay ratios. The table indicates that the measurement delay and control execution delay in the DCPPS cannot be aggregated into a synthetic delay, even when there is no communication delay, which means the delay merging property in centralized mode [16] does not apply to a DCPPS. The results for the above delay also verify the specificity of the CPPS under fully distributed control.
In summary, the stability of a DCPPS with multiple sources of time delay depends on not only the absolute values of the time delays but also the relative values and coordination of the time delays. To achieve improved stability performance, multiple sources of time delay in a DCPPS must be simultaneously considered. In addition, the delay merging property in centralized mode cannot be applied to a time-delayed CPPS under fully distributed control.

E. Sensitivity Analysis
The sensitivities of the rightmost eigenvalue λ r of the DCPPS in the marginally stable state with respect to the measurement delay τ m , the communication delay τ m and the control execution delay τ e in different directions based on the spherical coordinates in (33) are computed based on formulas (26), (28) and (29), as shown in Table IV. To study the effects of different sources of time delay on stability, we assume that delays of the same type are equal. Table IV shows that the real part of the sensitivities of the rightmost eigenvalue λ r with respect to local delays, such as measurement delay and control execution delay, is positive in any direction, which means that an increase in these local delays makes the rightmost eigenvalue tend to fall in the right half plane; this change disrupts the stability of the DCPPS in the marginally stable state. It is noteworthy that in some  Fig. 5. In the direction (90 • , 30 • ) based on spherical coordinates, there is no delay in the control execution tasks, and the delay vector is located in the lower right part of the line CD in Fig. 5 on the plane formed by the communication delay and measurement delay. In this case, when the communication delay τ c increases, the DCPPS tends to be stable. However, when the delay vector is located in the upper left part of line CD in Fig. 5, such as in the direction (90 • , 60 • ), increasing the communication delay will make the system unstable.

V. CONCLUSION
Considering cyber-physical interdependence, a general model of a time-delayed CPPS under fully distributed control is proposed in the form of DDAEs based on a distributed information flow model and dynamic model of power systems. The small-signal stability of a DCPPS with multiple sources of time delay is analyzed based on the critical eigenvalue tracking algorithm. An islanded microgrid from a real-world distribution network in China under fully distributed frequency control is tested in the case study. The accuracy of the obtained stability margin is verified by a simulation of the tested case in Simulink. Case studies illustrate that with multiple sources of time delay, the stability of a time-delayed DCPPS depends on not only the absolute value of the time delays but also the relative value and coordination of time delays, which means that decreased delays do not necessarily have a positive impact on stability. The theoretical work presented in this article can guide the optimization of delay threshold setting and the planning of measurement schemes, control units and communication links in DCPPSs in the future. In addition, the stability of time-varying DCPPSs is worthy of further study.

A. Sensitivities With Respect to Time Delays
It follows from (21) that ⎛ where (λ, v) denotes the eigenpair of (19) and v is the right eigenvector. For convenience, equation (34) can be rewritten in the following abstract form: Similarly, we have where u refers to the left eigenvector of (19) and u H denotes the Hermitian transpose of u. u can be easily obtained by solving the characteristic equation that includes the Hermitian transpose ofÃ.
To analyze the sensitivity of the eigenvalue λ with respect to various time delays, the partial derivative of the left-hand side of equation (34) with respect to τ m,k , k = 1, . . . , N c can be expressed as follows: By combining the partial derivative of the right-hand side of equation (34) Thus, the eigenvalue sensitivity with respect to the measurement time delay τ m,k , k = 1, . . . , N c is given as shown in formula (26), which completes the derivation of formula (26). Similarly, the sensitivities (28) and (29) can be obtained.

B. Sensitivities With Respect to System Parameters
By taking the partial derivative of (34) with respect to parameter p and multiplying by u H , we have Thus, the sensitivity of eigenvalue λ with respect to parameter p in formulas (30) and (31) can be determined.

APPENDIX B
The local measured physical state variable of fully distributed frequency control refers to ω i , i = 1, . . . , N c in this example. According to the measurement task of the distributed information flow and considering the local measurement time delay, the cyber measurement variable of agent i is defined as follows: where ω i denotes the rotor angular speed of generator i and M x.i = I when there is no failure or disturbance in the measurement task. Moreover, to combine the cyber side and physical side and consider the communication time delay among neighboring agents, the fully distributed frequency control algorithm with the communication time delay is written as the following differential equation: where λ i denotes the ICR of the cost of generator i, λ i = 2a i P m,i + b i , P m,i refers to the shaft power input of generator i, and ω N is the synchronous speed. μ ij ∈ U is the communication coefficient of the subgradient algorithm, and α is the step size for frequency recovery. Element μ ij in communication matrix U is defined by the following equation: where a i and a j are the numbers of agents connected to agent i and agent j, respectively. According to [18], the defined matrix U is a doubly stochastic matrix that can ensure the convergence of the consensus algorithm. By combining (42) and (43), a model of the differential equation of the distributed algorithm that considers both measurement and communication time delays is obtained as follows:λ The cyber control variable λ i in distributed frequency control affects the physical control variable P m,i . Therefore, the execution task considering the time delay in this example is expressed as follows: Similarly, C i = I when there is no failure or disturbance in the execution task.