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

2.1 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.

2.1.1 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]:

(1) θ ̇ i = ω ̃ i , i = 1,2 , … , n

(2) M i ω ̃ ̇ i = P i − ∑ j = 1 , j ≠ i n C i j ⁡ sin θ i j + D i j ⁡ cos θ i j − M i M T P COI ︸ = f i ( θ ) ,

where f _i (θ) is the acceleration power, θ _i = δ _i − δ ₀ 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. δ ₀ and ω ₀ are the center of angle and center of speed defined by:

δ 0 = 1 M T ∑ i = 1 n M i δ i ;

ω 0 = 1 M T ∑ i = 1 n M i ω i ;

M T = ∑ i = 1 n M i .

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.

The variables P _i and P _COI in (2) are given by:

P i = P m i − E i 2 G i i ;

P COI = ∑ i = 1 n P i − 2 ∑ i = 1 n − 1 ∑ j = i + 1 n D i j ⁡ cos θ i j .

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 C _ij = E _i E _j B _ij and D _ij = E _i E _j G _ij , 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.

2.1.3 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 c r , 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 V _cr = V _PE,max.

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:

1. f T ( θ ) ⋅ ( θ − θ s ) < 0 within the region of stability,

2. f T ( θ ) ⋅ ( θ − θ s ) = 0 at the crossing point and

3. f T ( θ ) ⋅ ( θ − θ s ) > 0 outside the region of stability.

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).

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.

2.2 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

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 Y ̄ 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 Y ̂ s and vector I ̂ 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 Y ̂ s in (9) is computed at the beginning of the simulation; the matrix product Y ̄ i T ⋅ Y i − 1 in (10) is pre-computed; the matrices Y ̂ 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.

Figure 4:

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

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.

3.1 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:

Figure 5:

Overview of the PEBS-based direct method workflow.

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 x d ′ 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 r r , 0 − 1 and the inverse of the modified admittance matrix Y r r , mod − 1 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 c r . 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.

3.2 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.

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 ₁, p ₂) and the interconnect partition.

Figure 6:

Illustration of communication aspects for a two-partition system showing the initialization step (t _init) and one simulation time step (t ₀ − t ₆).

3.3 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-by-step 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].

4.1 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 three-phase 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 time-domain 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 time-domain 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.

4.2 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.

4.2.1 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.

Figure 7:

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

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.

Figure 8:

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

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.

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.

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.

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.

4.3 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 SIME-based method, i.e. 8 times faster for the smallest network and 2.5 times faster for the largest.

4.4 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.