A two-stage parallel hybrid method for power system dynamic security assessment

: The present paper introduces a hybrid method based on both the Lyapunov’s direct approach using the Potential Energy Boundary Surface method and parallel time-domain simulation for fast power system dynamic securityassessment.Themaincontributioninthepresented methodisanapproachthatcombinesthedirectionalderiva-tive of the potential energy function and the maximum potential energy criterion in order to reduce the conserva-tiveness in the estimation of the stability region at the direct method stage. For detailed investigation of the critical contingencies,a parallel approach isapplied at the time-domain simulation stage. From the simulation results, the proposed method achieves a high level of accuracy in classifying network contingencies and shows great computational performance potential in light of the requirement for fast dynamic security assessment.


Introduction
Transient stability analysis is an important but challenging task for the reliable and secure operation of power systems.Simulation tools are frequently used for system analysis with a significant level of simplification [1].Accurate analysis is carried out using time-domain approaches, which rely on the step-by-step numerical integration of the non-linear differential algebraic equations describing the system.However, applying time-domain simulations for a complete system analysis -consisting of studying a number of disturbances -requires large computation time.This in turn imposes constraints on the number of cases and scenarios that can be studied during the analysis process.In view of changing network conditions such as continuous load variations and the current operation of grids with a large portion of generation from intermittent renewable energy sources, there are continuous deviations in the operating point of the network.Therefore, continuous analysis becomes a vital task to ensure system stability and reliability.
Alternative approaches using direct methods based on the equal area criterion or Lyapunov's method have been proposed for fast stability assessment: Direct methods achieve computational efficiency in the stability assessment process by using algorithmic techniques to compare the system energies during the fault period and at the initial post-fault state; thus, eliminating the time-consuming numerical integration of the system differential equations in the post-fault period.Ribbens-Pavella et al. [2] present direct methods, with emphasis on using the Lyapunovbased method, for studying dynamics of large-scale electric power systems.Lyapunov's direct method has further been used for specific applications such as fast stability analysis and efficient classification of contingencies [3,4].The direct methods based on the equal area criterion have also found application in the area of contingency filtering as presented in [5].A comparison of Lyapunov's direct method and equal area criterion method is given in [6] in an effort to assess the adaptability of the methods for real time operation.
The main drawback of direct approaches is that they do not derive the actual system trajectory during the analysis process.Hybrid approaches, combining the advantages of direct methods and time-domain simulations, are considered to provide improved speed and accuracy of transient stability analysis.Zhang et al. [7] describe an implementation of a hybrid method, known as Single Machine Equivalent (SIME).The method has been applied for different contingency filtering and assessment studies as described in [8,9].Recent studies have seen a renewed interest in application and improvement of the SIME assessment approach [10,11].The direct method used in the SIME approach is based on the equal area criterion.From the computation point of view, SIME applies time-domain simulations during the fault-on period and partly during the postfault period in order to rank the machines in the system.This implies that the computational time is increased by the application of numerical integration during the post-fault period.According to the assessment of the direct methods in [6], the Lyapunov-based method has a potential of reducing the computational time by limiting the explicit numerical integration only to the fault-on period.This results in faster contingency selection criterion and therefore speeds up the analysis process.
The present paper introduces a hybrid dynamic security assessment method based on coupled Lyapunov's direct approach using the Potential Energy Boundary Surface wie oben (PEBS) method and parallel time-domain simulations.The hybrid method takes advantage of the two approaches to achieve fast contingency screening and selection, and further detailed analysis of critical contingencies.The first stage of the analysis applies the PEBS-based method to derive the system state in terms of stability or instability for given contingencies and identifies critical contingencies.In the second stage, the critical contingencies are further analyzed by applying a time-domain simulation approach.The main contributions in the present paper are: (1) An approach that combines the directional derivative of the potential energy function and the maximum potential energy criterion in order to reduce the conservativeness in the estimation of the stability region at the direct method step.(2) A parallel approach applied at the time-domain simulation stage for detailed analysis of the identified critical contingencies.
The remainder of this article is organized as follows: Section 2 describes the formulation of the two stages of the hybrid method.The implementation of the components of the hybrid method and the workflow are described in Section 3. In Section 4, simulation results are illustrated to evaluate the accuracy and analysis process of the hybrid method.Finally, Section 5 gives a summary of the paper and discussion of future work.

Hybrid method
The hybrid dynamic security assessment method proposed in the present paper is a two-stage method based on both a Lyapunov-based direct approach and a parallel timedomain simulation approach.Figure 1 illustrates the stages of the proposed hybrid method.The first stage of the method is a screening step.This uses the PEBS-based Lyapunov's method to derive the state of the system in terms of stability or instability for each contingency, thus identifying the critical and noncritical contingencies.At this stage, numerical integration is applied to derive the system trajectory only during the fault-on period.In addition to establishing the state of system stability, it is important to identify, quantify and visualize the actual system behavior during the transient and post-fault period for the critical contingencies.Therefore, the second stage of the hybrid method uses a parallel time-domain simulation approach for further analysis of the identified critical contingencies using detailed modeling and more accurate step-by-step numerical integration.The following sections describe the formulation of the methods applied at each stage of the hybrid analysis method.

Formulation of the direct stability assessment approach
The direct method is applied at the first stage of the proposed hybrid method as shown in Figure 1.Stability analysis using direct methods involves two main steps.The first step computes the relationship between the pre-fault stable equilibrium point (SEP) and the system state at the time of fault clearing.This is executed by numerical integration of the system equations during the fault period to derive the fault-on trajectory until the fault clearing point.The second step directly determines, without numerical integration of the post-fault system equations, whether the initial state of the post-fault system lies inside the stability region of a desired stable equilibrium point.An overview of the direct methods for transient stability analysis in power systems is given in [12].The direct method implemented in the approach presented in this paper is based on the PEBS method.This method directly determines system stability based on a transient energy function defined for the post-fault system and a region of stability relative to the disturbance under consideration [13,14].The transient energy function is defined depending on the system model, whereas suitable region of stability is defined using a critical energy relative to the fault-on trajectory as described in the following subsections.

System model representation
The model of the system used in the presented direct method is based on the following assumptions: classical generator model represented by constant voltage behind transient reactance; constant input mechanical power; and constant impedance loads.Figure 2 illustrates the interface between a classical generator model and the network.The dynamics of the ith generator are represented by the equation motion in the center of inertia (COI) reference given by ( 1) and (2) [15]: where f i () is the acceleration power,  i =  i −  0 is the relative angle, and ωi =  i −  0 is the relative rotor speed in COI reference. i and  i are the rotor angle and speed deviation of the ith generator in synchronous reference frame, respectively. 0 and  0 are the center of angle and center of speed defined by: The variables P i and P COI in (2) are given by: M i is the moment of inertia; M T is the total inertia of the system; and P mi is the mechanical input power; The terms C ij and D ij are derived from the electrical power output and are given by , where E i and E j are constant voltages behind direct axis transient reactance of machine i and j.The terms B ij and G ij are the transfer susceptance and conductance of the (i, j)th element in the reduced admittance matrix of the system, respectively.During stability analysis, the system is divided into the following states: Pre-fault state -system state before the disturbance; Faulted state -system during the disturbance; Post-fault state -system after clearance of the disturbance.The equations defined in (1) and ( 2) are used to describe the system during the three states.

Transient energy function
The energy function used in direct methods for assessing transient stability was derived by [15] and is commonly referred to as the transient energy function.The general form of the transient energy function is given by The variable  s is the post-fault stable equilibrium point (SEP).According to [15], the four terms on the right in the transient energy function given in (3) can be defined with respect to their physical interpretations as follows: Term 1 is the total change in rotor kinetic energy; Term 2 is the change in potential energy of the rotor; Term 3 is the change in stored magnetic energy in the branch between generators i and j; Term 4 is the change in dissipated energy in the branch between generators i and j.The dissipated energy defined by Term 4 depends on the actual path from a stable point to the final fault clearing point.This term can only be evaluated if the system trajectory is known or approximated through specific assumptions.One way of approximating this term is using the linear trajectory approximation.In this approximation, a linear trajectory of the system is assumed between the initial stable equilibrium point and the final angle [15].The resulting expression for evaluating the change in dissipated energy is given by (4) [15]. where ) .

Derivation of transient stability region
An important step in assessing system stability using the direct approach is the estimation of the region of stability.The region of stability is defined by a set of initial states from which the system trajectory converges to a stable equilibrium point.The direct assessment of system stability involves analyzing the location of the last state of the system trajectory at the point of disturbance clearance with respect to the region of stability.For power system applications, the critical energy of the system is used as an approximate measure of the stability region.In other words, the critical energy defines the level of energy as the system trajectory exits the region of stability.The interior of the stability region is then defined by the inequality V(, ω) < V cr , where V(, ω) is an energy function, and V cr is the critical energy [16].
As described in Section 2.1.2,the transient energy includes kinetic energy due to the rotation of the generators, and potential energy component as a result of the rotor position, stored magnetic energy and energy dissipated in the elements of the system.The PEBS method as proposed in [15,17] takes advantage of the comparison of the system energy.According to [17], the critical value of the energy function is estimated by following the values of function V(, ω) and V PE () at every time step during the fault-on period until the maximum value of the potential energy V PE,max is reached.This point is considered to be the boundary crossing and gives an estimate of the critical energy expressed as However, the first maximum value of the potential energy does not always indicate the true maximum point in the full state space.Figure 3 illustrates such a scenario with more than one maximum point.Therefore, an additional criterion is necessary to prove the existence of the true maximum critical point of the system within a specified time of assessment.The approach presented in the present paper additionally applies the directional derivative of the potential energy function along the fault-on trajectory in the angle space to derive the true crossing point of the stability boundary [15].This is mathematically represented by the inner product f T () ⋅ ( −  s ), where f () is the acceleration power of the fault-on system in (2), referred to as the PEBS function in the rest of the present paper.The regions around the stable point are defined as follows: f T () ⋅ ( −  s ) < 0 within the region of stability, -f T () ⋅ ( −  s ) = 0 at the crossing point and f T () ⋅ ( −  s ) > 0 outside the region of stability.
In the presented combined approach, the potential energy and directional derivative of the potential energy function are monitored at every time step until a change in sign is observed in the PEBS function, i.e. from negative to positive; this point defines the true stability boundary crossing at which  =  * .Once the exit point is determined, the maximum of the potential energy profile from the stable equilibrium point to the crossing point gives an estimate of the critical energy V cr = V PE,max ().
Figure 3 illustrates the estimation of the critical energy using the presented approach.As depicted in the figure, the first maximum potential energy V PE,max1 represents a local maximum point of the potential energy.The additional test using the directional derivative identifies another maximum point V PE,max2 before the exit of the stability region.In this case, the maximum of the two detected points defines the global maximum potential energy, which gives the estimate of the critical energy; and thus defines the stability boundary.This method eliminates the computationally challenging requirement of defining the controlling unstable equilibrium point (controlling UEP) used for determining the critical energy during the assessment process.

Formulation of the parallel time-domain simulation
The second stage of the proposed hybrid method is used for detailed analysis of selected contingencies, which requires more detailed models and more accurate solution techniques.In this case, the time-domain simulation approach is applied for the further detailed analysis as shown in Figure 1.
In time domain simulations, the solution is based on the step-by-step numerical integration of the equations to compute the system trajectory over time.The computational speed and accuracy mainly depend on the model details and simulation step size, which in turn define the time scale of the system dynamics that can be captured during the simulation.The time-domain simulation approach used in the presented hybrid method is specific for electromechanical transients analysis, assuming fundamental frequency for the network system quantities.Generally, the power system model is represented by where x is a vector of dynamic state variable, y is a vector of algebraic variables, and u is a set of model parameters.These differential and algebraic equations represent the interaction between the different components.Further details about the models and solution methodology of the time-domain simulation used in the hybrid method are described in [18,19].The interface between the generators and network system is via the stator algebraic equations included in the network nodal equation given by YU = I, (6) where Y is the network nodal admittance matrix, U and I are voltages and current injections vectors, respectively.Analysis of the computation time of the individual stages in time-domain simulations shows that the solution of ( 6) takes up a large percentage of the total runtime [20].
In the present paper, a block bordered diagonal form (BBDF) [21] is used to reformulate the network nodal equation into a form that can be parallelized during the network solution step.This restructures the network in a way that the nodes are grouped into p + 1 sub-blocks; where p is the number of subnetworks and the (p + 1)th sub-block consists of the boundary nodes.The subnetworks are created from partitioning the grid.The restructured network nodal equation is given by where Y i are the elements of the original Y matrix within subnetwork i; Y s is the matrix formed by the boundary nodes in the (p + 1)th sub-block.The Ȳi elements define data of the branches connecting subnetwork i to the boundary nodes.
In the new formulation, the solution of the network nodal equation in ( 7) is rearranged into two steps; the preprocessing step and one step for every subnetwork of the matrix given in ( 8)- (11).
The preprocessing step is given in (8), where the matrix Ŷs and vector Îs are calculated in ( 9) and (10), respectively.
This step is used to compute the boundary node voltages U s required for the second step to separately solve for the node voltages U i in the ith subnetwork as given in (11).The solution of (11) for the subnetwork node voltages constitutes the parallelizable task in the network solution implemented in the parallel time-domain simulation.The solutions of the main equations ( 8) and ( 11) are based on the LU factorization of the respective admittance matrices.
The following steps are additionally taken to optimize specific computations in the preprocessing step: The solution of matrix Ŷs in ( 9) is computed at the beginning of the simulation; the matrix product 10) is precomputed; the matrices Ŷs and Y i are factorized by LU decomposition beforehand to speed up the calculation during the simulation loop.
Solving the block bordered diagonal form initially requires a sequential solution for a vector of boundary node voltages from (8).The complexity of the solution grows with more than linear complexity.Therefore an efficient partitioning strategy of the network is important in the formulation of the parallel solution.The partitioning should result in a balanced number of nodes in the subnetworks and a minimum number of boundary nodes in the (p + 1)th subnetwork (interconnect partition) for an optimal performance.As described in [22], a graph partitioning strategy is applied in the presented method and extended for application to dynamic simulations.The network partitioning is summarized in Figure 4.
The formulation in Figure 4 shows that the equations in the subnetworks can be solved in parallel.The equations formed by the interconnect partition have to be solved sequentially at every time step.The sequential part of the solution depends on the size of the interconnect partition.The number of nodes in the interconnect partition is optimized by minimizing the number of cut branches between partitions from the basic partitioning format.

Hybrid method implementation
The new hybrid method is implemented in Matlab and Julia [23].The input power system networks applied in the simulations are described using the Matpower case format [24].The present section describes the implementation of the PEBS direct method, the solution process of the time-domain simulation, and the coupled simulation workflow of the hybrid method.

Direct method implementation
The implementation of the direct method for selection of contingencies is based on the PEBS approach introduced in Section 2.1.3.The following steps summarize the process of the implemented direct approach for the selection of network contingencies as shown in Figure 5: Step 1: The initial step is to compute the power flow of the pre-fault network in order to determine the steady state values of the system and the initial values of the dynamic variables (, ω).The power flow computation is executed using the Matpower power flow algorithm.The time-derivatives of the state variables are set to zero at the initialization stage.
Step 2: The parameters of the model and the system are set up for the pre-fault, faulted, and post-fault states.The load impedance is connected between the load bus and the reference node.Internal generator nodes are created as additional nodes, and the internal generator reactance connected between the internal nodes and the generator terminal nodes.For the faulted system, the fault impedance simulation is added depending on the location.The admittance matrix is recomputed at each switching operation.
The nodal admittance matrix Y of the network is formulated for each state as described in [24].A reduced internal network matrix Y int is formed for each state by eliminating all nodes (r) from the nodal admittance matrix except the internal generator nodes (n) by Kron reduction [25].The resulting matrix is given as The network reduction in (12) is carried out at every occurrence of network modification and scales with more than quadratic complexity, dominated by the inversion of matrix Y rr .In the method presented in this paper, the formulation of the reduced network admittance matrix following a network modification is simplified using the compensation method based on the matrix inverse lemma described in [26].This simplifies the inverse computation using a relation between the inverse of the original network admittance matrix Y −1 rr,0 and the inverse of the modified admittance matrix Y −1 rr,mod given by where Δy defines the modifications in the parameters of the original network, and M is a connection matrix based on the type of modification [27], i.e. for branch oriented modifications a column vector is defined for each modified branch in M with +1 and −1 entries on the respective nodes and 0 elsewhere, whereas for a short circuit connection of a bus to ground, +1 is defined for the corresponding bus position in the M matrix.Formulation of the modified network reduced internal admittance matrix Y int,mod is thus given by Step 3: The post fault stable equilibrium point (SEP) is determined from the solution of f i () = 0.It is important to determine whether the equilibrium point of the system is stable or unstable after the disturbance is cleared.Therefore, the stability of the equilibrium point is tested using the linearization principle about the equilibrium point, and computing the eigenvalues of the system dynamics' state Jacobian matrix as follows: Linearizing the state equations in (1) and ( 2) results in a state matrix of the form given in (15).
At the equilibrium point, the system is steady state stable if all the eigenvalues  x of the characteristic equation have negative real parts [28], i.e.R{ x } < 0, for all x.The equilibrium point is unstable if R{ x } > 0, for any x.If the stable equilibrium point does not exist, the corresponding switching operation is considered to be an unclassifiable case for direct stability analysis.The stability of such contingencies is analyzed using time domain simulations.
Step 4: For the system with stable equilibrium points, the fault-on trajectory given by the power system model equations is computed by numerical integration.The numerical integration algorithms in the time-domain simulation approach are applied at this stage.The PEBS crossing point is determined by monitoring the inner product f T ( t ) ⋅ ( t −  s ) and the value of the potential energy V PE () at each time step.The value of  is obtained from the fault-on angle and admittance matrix used for the calculations is based on the post-fault system.The PEBS crossing is reached when f T ( t ) ⋅ ( t −  s ) = 0.The critical energy, given by V cr = V PE ( t ), defines an estimate of the region of stability for the post-fault system.
The state at each time step during the integration of the fault-on system is tested to determine whether it lies within the stability region given by V( t , ωt ) < V cr .Therefore, the value of the energy function V( t , ωt ) is computed at successive times using the fault-on  and ω obtained from the step-by-step integration of the fault on trajectory.The point at which the transient energy V( t , ωt ) is equal to the critical energy V cr defines the critical clearing time (CCT).The stability margin of a critical contingency is given by S m = V cr − V cl , where V cl is the total energy of the system at the clearing time.Figure 5 shows the direct method simulation process.

Implementation of the parallel time-domain simulation approach
The parallel time-domain simulation algorithm is implemented in Julia programming language, which offers good performance and great parallelization potential [23].The main steps of the parallel algorithm are summarized in Algorithm 1.The inputs to the algorithm include the network description casefile in Matpower format and the subnetworks formed after the network partitioning process.The first step in the simulation is the initialization step to determine the initial dynamic and algebraic values required to establish a quasi-steady state starting point for the simulation.In the presented implementation, steady state power flow calculations are used to define the initial operating conditions.Additional inputs for the parallel dynamic simulation are the dynamic model parameters and the events file that defines the modifications in the network.
It should be noted that the solutions in the dynamic simulations are based on the step-by-step numerical solution.Therefore the parallelization is limited to a single time step.Since the generators in the system are naturally decoupled from each other, the first parallelization step is the computation of the decoupled machine differential equations.For this, the fourth order Runge-Kutta method is used to obtain the node current injections in the different subsystems.The second parallel step is the computation of the BBDF formulated network equation.The solution of the network consists of the precomputation steps, which are mainly matrix construction steps required for the sequential solution of the interconnect partition equations, and the parallel solution of the subnetwork equations as described in Section 2.2.The sparse LU factorization solver using the UMFPACK library [29] is used for solving the linear network system.Algorithm 1. Parallel computation procedure.The parallel time-domain simulation algorithm is based on a single node parallelization approach.The main computations in the algorithm are memory bound tasks dealing with vector arithmetic, matrix multiplication, and solving of the network equation.In the Julia environment, such a parallelization problem is effectively handled using the multithreading construct [23]. Figure 6 illustrates the communication process during each time step of the simulation for an example with two partitions (p 1 , p 2 ) and the interconnect partition.

Workflow of the hybrid method
The approach for dynamic security assessment presented in this paper is implemented by coupling the two algorithms; i.e. running the PEBS-based direct method in the first stage and the parallel time-domain simulation in the second stage.In the first stage, the simulation time is set according to the protection settings of the network.The goal of the simulation at this stage is to identify contingencies that result in system instability before isolating the cause of the contingency from the network.If the disturbance results in a gain in transient energy greater than the critical energy, it is considered to be a critical disturbance and added to a list of contingencies requiring further analysis.Otherwise, the disturbance is considered to be noncritical with no further analysis required.
The unclassifiable cases-identified in Step 3 of the direct method process-and the critical contingencies together with the respective estimated critical clearing times are sent to the time-domain simulation module in the second stage for further analysis.This involves the step-bystep simulation of all the selected contingencies.The fault clearing time is set according to the critical clearing time obtained in the first stage for the contingency under consideration.However, since no critical clearing time is derived for the unclassifiable cases in the first stage, the initial fault clearing time for these cases is set equal to one simulation step size.The aim of the analysis at this stage is to derive the actual transient response of the system during disturbance and in the post-disturbance period.
Algorithm 2 summarizes the steps of the hybrid dynamic security assessment process.The two stages are considered to be complementary to each other.Whereas the first stage provides a fast classification of disturbances and provides the stability margin of the network for each contingency, the actual transient response of the system is analyzed in the second stage.The coupled analysis therefore provides a quick selection and detailed analysis of only the dangerous contingencies.
For interactive modelling and visualization of assessment results, the hybrid method is coupled to the eASiMOV framework [30].The eASiMOV software framework -comprising of energy system analysis, simulation, modeling, optimization and visualization components -provides an extendable framework for interactive modeling and analysis of energy systems.The hybrid method is integrated into the eASiMOV framework as a designated part of the ePowSim simulation module.Details of the modules in the framework are given in [30].
Algorithm 2. Overview of the simulation steps in the hybrid dynamic security assessment method.

Evaluation of hybrid method
This section describes the simulation tests performed to evaluate the accuracy and computational performance of the presented hybrid method.The evaluation mainly focuses on the direct method stage.This is carried out by comparing the results obtained using the direct method to the results from the validated time-domain simulation method in terms of the ability to categorize network contingencies and identification of the critical clearing time of disturbances.The time-domain simulation results are used as the benchmark estimates in the present paper.
In the first test case, the method is applied to the standard IEEE 9-bus test network to verify the accuracy of the estimated critical fault clearing time.In the second test case, the hybrid method is tested using the model of the Baden-Württemberg state (Germany).This test focuses on the analysis process in the proposed approach.Additional networks with varying complexity are used to evaluate the accuracy, reliability and computational performance of the proposed method.The following networks are used for the evaluations: Case68, Case118 and Case1354, which are standard IEEE test networks, and CaseBW, CaseDE, which are network representations of the transmission grids in the Baden-Württemberg (Germany) state and in Germany, respectively.In all cases, the tested systems are limited to balanced network configurations.

Test case I: IEEE 9-bus network
The standard IEEE 9-bus network consists of three generators, nine buses, nine branches and three loads.The parameters of the generators and network structure are derived from the IEEE 9-bus network described in [31].The parameters of the generator considered in the test case are given in Table 1 showing the inertia constant H, damping constant D, steady state reactance x d and transient reactance x ′ d .The per unit (pu) values are expressed on a base of 100 MVA.
For the contingency analysis in this section, a threephase short circuit fault is applied on the different buses in the network.The complete contingency is defined by the location of the fault and switching action applied to clear the fault.In this case, a total of 18 contingencies are tested as shown in Table 2.The simulation time in each case is set to 2.0 s, which represents the minimum operation time of the protection devices.Practical power system protection settings present much lower clearing times.In the present simulation case, a setting of 2.0 s is used for testing purposes to assess the accuracy of the method in comparison to timedomain simulations.The step-by-step numerical integration stage of the direct method is computed with a step size of 1 ms.The same step size is used for the time domain simulation stage.
The simulation results comparing the critical clearing time obtained using the direct method approach (DM) and the step-by-step time-domain simulation approach (TDS) are shown in Table 2.The table includes the simulated contingencies (in terms of the fault bus and line trip action), the critical energy at the PEBS crossing (V cr ) and the Critical Clearing Time (CCT).The dash (−) in the line trip column indicates cases with no line switching action to clear the fault; the fault is considered to be self-clearing in such a case.The contingencies are classified using the relationship between the critical clearing time (t cct ) and the set simulation time (t cl = 2 s).The classification is summarized as follows: Critical contingencies if t cct < 2 s; Non-critical contingencies if t cct ≥ 2 s (stable).
From the results shown in Table 2, the direct method categorizes all the simulated cases as critical contingencies.This is observed from the estimated critical clearing time (CCT-DM), which is less than the fault duration of 2.0 s in all cases.A similar categorization is derived by the timedomain simulation approach (CCT-TDS).The critical clearing times obtained using the direct method are in close agreement with those from the time domain simulations based on numerical integration.The notable differences in the critical clearing time estimates of the two approaches can be attributed to the conservative estimations from the direct method.
The simulation results with the IEEE 9-bus network model show that the critical energy and critical clearing time depend on the location of the fault and nature of fault clearing.

Test case II: real power system
The method is further tested using the Baden-Württemberg network representing a real power system.The modified network consists of 17 generators, 150 buses and 225 branches.The generators are connected through step-up transformers to the transmission network at voltage levels of 220 kV and 380 kV.The loads, which represent the aggregated loads in distribution subnetworks, are at 110 kV buses and connected via transformers to the transmission network.A network representing the transmission grid of the Baden-Württemberg state is described in [32].
The test network consists of 568 contingencies.For testing purposes, the simulation time is set to 2.5 s in the first stage for each test case.The direct method categorizes the tested network contingencies as follows: 28 unclassifiable contingencies; 409 noncritical contingencies; and 131 critical contingencies.

Analysis of a critical contingency
A critical contingency is considered in this case to illustrate the analysis process in the hybrid method.A detailed analysis of selected cases in all categories (i.e.infeasible, noncritical, and critical contingencies) is given in [33].For the critical contingency in this case, a fault is simulated near bus 95, and cleared by switching out the line between bus 95 and 90.The plot of the eigenvalues of the resulting state matrix at the post-fault equilibrium point is shown in Figure 7.In this case the state matrix has no eigenvalues with positive real parts.This implies that the post-fault equilibrium point is stable.
The next step in the process is to monitor the PEBS function, potential energy and transient energy function as shown in Figure 8.The zero crossing of the PEBS function occurs at approximately 0.58 s.The potential energy at this time gives an estimate of the critical energy V cr = 6.353 pu.
The simulated fault is classified as a critical contingency.The critical clearing time of the fault is estimated at the point when V total = 6.353 pu, resulting into t cr = 0.471 s. 131 faults are classified as critical contingencies.The critical contingency specified above is then further analyzed using time-domain simulations.Figure 9a shows the analysis of the generator rotational speed response in the second stage of the simulation method.Furthermore, the corresponding response in the interactive visualization module of the eASiMOV framework is shown in Figure 9.
To determine the critical clearing time using the time domain simulation, the fault clearing time is successively varied until diverging responses are observed in the post fault period.Figure 10a shows the speed response at the clearing time where instability first occurs.The estimated critical clearing time in this case is 0.482 s, with a relative difference of 2.282 % compared to the estimate from the direct method.Figure 10 shows the details of the critical network section in the visualization framework.Such analysis provides more insight regarding the actual state of the network during the dynamic security assessment process.

Accuracy and performance
The accuracy, reliability and computational performance of the proposed method are evaluated in comparison to the time-domain simulation method and the state-of-the-art SIME-base method [11].The platform for the simulations is an Intel(R) Core(TM) i7-8700 CPU @ 3.20 GHz, 32 GB memory system, running on 64-bit Microsoft Windows 10.The accuracy and reliability in contingency selection is assessed in reference to the classification by the step-by-step time-domain simulation approach.For a specific network, assessment in the time-domain simulation is based on the angular separation between the machines in the network.A point of instability is reached if the rotor angular separation of any machine from the rest of the system is greater than a given angular threshold.A common measure of the maximum angular separation of generators before loss of synchronism is 120 • [34].Four test cases of the IEEE 118-bus network with different fault clearing times are considered here, i.e.Case I -2.0 s, Case II -1.0 s, Case III -0.5 s, Case IV -0.2 s.
Table 3 shows the contingency selection results by the time-domain simulation approach and the accuracy of the proposed method with respect to the benchmark results of the time-domain simulation approach.The new method shows a high level of accuracy in classifying the network contingencies.In this specific case, the method identifies unstable contingencies with an accuracy of 100 % in all scenarios.This implies that all critical contingencies in the network are captured during the selection process.However, less accurancy is attained for the noncritical contingencies (see Table 3 direct method accuracy).This value represents the number of selected noncritical cases by the presented method compared to the actual stable (noncritical) cases identified by the time-domain simulation.In this case, the method shows a small percentage of conservativeness whereby noncritical contingencies are falsely classified as critical (i.e.false critical cases).
Furthermore, the accuracy, reliability and computational performance of the proposed assessment method is are compared to the state-of-the-art SIME-based method described in [11] using networks of different complexity as shown in Table 4. Accuracy is measured by the ability to correctly select noncritical contingencies, whereas  reliability is measured by ability to correct identify critical cases that require detailed analysis.In each case, the fault clearing time is set to 200 ms to measure the computational runtime.Table 4 shows the comparison of the two methods with respect to benchmark time-domain simulations.The proposed assessment method shows better accuracy and reliability than the SIME-based method in all the tested cases.In addition, the proposed method shows better computational performance per contingency than the SIMEbased method, i.e. 8 times faster for the smallest network and 2.5 times faster for the largest.

Discussion
The simulation results with the IEEE 9-bus network model show that the critical energy and critical clearing time depend on the location of the disturbance, nature of clearing the disturbance and system parameters.Therefore, continuous system analysis is necessary for defining the required network settings in presence of changing operating conditions.The test case using the Baden-Württemberg system shows the selection criterion under varying fault clearing times in the new hybrid method.A selected case is analyzed to illustrate the underlying procedure during the selection process.The hybrid method generally reduces the amount of time required for the time-domain simulation by limiting the detailed analysis to only potentially harmful disturbances.Furthermore, the results show that the critical cases decrease with shorter fault clearing times, which shows that the runtime for detailed analysis using time-domain simulations can further be reduced.The proposed hybrid approach can therefore be applied for continuous dynamic security assessment with changing network operating conditions to analyze the ability to maintain the N-1 criterion of the power system.

Conclusions
The present paper introduces a two-stage parallel hybrid method to address the need for fast dynamic security assessment.A direct approach based on Lyapunov's method is applied in the first stage for the selection of contingencies, and a quick estimation of the critical clearing time.
The selection stage takes advantage of the fast stability analysis capability of the PEBS method.The novelty in the presented method is an approach that combines the directional derivative of the potential energy function and the maximum potential energy criterion in order to reduce the conservativeness in the estimation of the stability region in the analysis process.From the evaluation results, it is shown that the presented method can accurately and reliably classify contingencies into critical and noncritical contingencies.The estimates of the critical clearing time in the direct method are in close agreement with the numerical integration approach in the time-domain simulations.In terms of computational performance, the presented method shows better performance than the SIME-based method, and therefore shows great potential in light of the requirement for fast dynamic security assessment.This is due to the fact that the direct method stage applied in the proposed approach limits the numerical integration to the faulton state, unlike the SIME hybrid method which continues to compute the trajectory partly beyond the clearing time.In a further step, a parallel time-domain simulation approach is applied to the selected critical contingencies for detailed analysis of the network response.The detailed analysis provides information that cannot be derived at the direct method stage regarding the actual system transient response.
Further work is, however, required to address the notable conservativeness identified at the selection stage.

Figure 1 :
Figure 1: Overview of the proposed two-stage hybrid method workflow for contingency screening and detailed analysis.

Figure 2 :
Figure 2: Illustration of the generator-network interface.The admittance matrix Y accounts for the generator transient reactances, in addition to the network and loads.

Figure 3 :
Figure 3: Combined approach for estimation of critical energy; monitoring the potential energy (V PE ) and directional derivative of the potential energy function (PEBS function).

Figure 4 :
Figure 4: Extension of network partitioning with an interconnect partition for a three-subsystem network.

Figure 5 :
Figure 5: Overview of the PEBS-based direct method workflow.

Figure 6 :
Figure 6: Illustration of communication aspects for a two-partition system showing the initialization step (t init ) and one simulation time step (t 0 − t 6 ).

stage 2 :
Inputs: import list of selected contingencies c define simulation duration = T for each contingency in c do while t < T do fault settings: bus = n; Δt = t cct ; line = b run parallel time-domain simulation compute detailed system response Return: State variables at each time step

Figure 7 :
Figure 7: Eigenvalues plot for stable equilibrium point of post-fault system with line 95-90 switched out.

Figure 8 :
Figure 8: Monitoring PEBS function, total energy and potential energy for a critical contingency.

Figure 9 :
Figure 9: Time-domain generator speed response and visualization of the response of system frequency in the eASiMOV framework for a fault cleared before the critical clearing point.Critical network section is shown with frequency scale ranging between 49.6 Hz and 50.4 Hz.(a) Generator rotational frequency response for fault clearing time t cl = 0.471 s.(b) Visualization of system frequency at t = 2.270 s.(c) Visualization of system frequency at t = 9.350 s.

Figure 10 :
Figure 10: Time-domain generator speed response and visualization of the response of system frequency in the eASiMOV framework for a fault cleared at a time beyond critical clearing time.Critical network section is shown with frequency scale ranging between 49.0 Hz and 51.0 Hz.(a) Generator rotational frequency response for fault clearing time t cl = 0.483 s(t cl > t cr ).(b) Visualization of system frequency at t = 2.230 s.(c) Visualization of system frequency at t = 10.00 s.

Table 1 :
Set of generator parameters.

Table 2 :
List of contingencies and corresponding estimates of critical clearing time (CCT).

Table 3 :
Accuracy of the proposed direct method for screening contingencies with respect to the time-domain simulation approach.

Table 4 :
Comparison of the proposed direct method to the SIME-based method.