# An analytical critical clearing time for parametric analysis of transient stability in power systems

###### Abstract

An analytic approximation for the critical clearing time (CCT) metric is derived from direct methods for power system stability. The formula has been designed to incorporate as many features of transient stability analysis as possible such as different fault locations and different post-fault network states. The purpose of this metric is to analyse trends in stability (in terms of CCT) of power systems under the variation of a system parameter. We demonstrate the performance of this metric to measure stability trends on an aggregated power network, the so-called two machine infinite bus network, by varying load parameters in the full bus admittance matrix using numerical continuation. Our metric is compared to two other expressions for the CCT which incorporate additional non-linearities present in the model.

## I Introduction

The complex dynamics of electric power systems have long been the subject of intense research particularly in the area of stability. Effective stability metrics provide control inputs and assist the system operator to ensure that a power system maintains synchrony after the network suffers a fault, i.e. that it exhibits transient stability. A traditional transient stability metric for short circuit faults on a power network is the so-called critical clearing time (CCT) [1, 2]. The CCT provides an upper bound on the duration of a short circuit on a power network before it is removed - ‘cleared’ - by the action of protection mechanisms to isolate the faulted circuit such that the system will regain synchronisation once the fault is cleared. In general, the CCT is a useful metric for power system design; by allowing the severity of different situations and the effectiveness of different interventions (generation re-dispatches, control modifications or network reinforcements) to be compared.

Currently, there are practical developments in power systems that promise to radically change power system dynamic behaviour. For example, the gradual substitution of power generated from large, synchronous machines by asynchronous machines or power fed via power electronic interfaces (e.g. wind farms, solar PV and HVDC interconnections to other systems), in addition to the changing nature of electrical loads [3]. Previous work in the literature [4] has investigated the effect of changing loads on system stability by repeating fault studies for different loading levels. As a consequence, there is value in articulating metrics that exploit theoretical, if simplified, descriptions of the system which can provide a deep understanding of the impact of a wide range of features of the network from parametric investigations. This can inform efforts to design strategies to mitigate possible instabilities in the system.

In the recent literature, alternative methodologies have been used to study stability when modelling a power system using the so-called swing equations [1, 5]. These include synchronisation [6], non-linear dynamics [7], bifurcation theory [8], passivity-based methods [9] and the computation of basins of attraction [10]. Direct methods [11] cast the swing equations in an energetic framework to provide a critical energy boundary for the whole system during a fault. Despite the difficulty of including non-negligible transfer conductances in direct methods [12], their advantages include the possible estimation of an analytical stability boundary and relatively quick computation. Also, they require no need for further simplifications of a power system beyond the swing equation model and they can be applied to any system that can be parametrised. The system operator can use this analytical stability metric for initial safety checks and to assess the stability margins of the system once a fault has been cleared.

One of the drawbacks of the direct methods is that it is difficult to predict when the system energy will cross the critical energy boundary because of the non-linear nature of the system dynamics. So-called fault trajectory sensitivity techniques [13, 14, 15] have been proposed to consider the effect of parameters on stability by linearising about the trajectory of a fault in state space with respect to a given parameter. Furthermore, a method for computing a so-called “direct CCT” has been proposed [16] which is based on linearising the power system model about a specific fault trajectory with respect to the system energy itself. An estimate of the CCT is then found by extrapolation. However, to our knowledge, an analytic CCT metric is only available for induction generators [17] and there is no analytic estimate of the CCT developed for a network of synchronous generators.

The aim of this paper is to propose a new analytic expression of the CCT. This estimate is derived by recasting the energetic metric used in the direct methods in terms of a metric in time by simplifying the energy functions and the dynamics during a fault. As is true for direct methods in general, our metric can serve as a lower bound to the true CCT for lossless power systems (or for power networks with small transfer conductances [18]) and can be applied to systems suffering a large fault at any location on a network. However, the purpose of our metric is to capture trends in stability as network parameters are varied and as such, the investigations in this paper are limited to aggregated or clustered power networks. (See [2, Chapter 14] for aggregation techniques.)

In general, a power network’s topology changes from its original structure when a fault is cleared. This is generally due to some switching action that isolates the region of the network that suffers the fault. Choosing the best strategy to quickly identify a need for and carry out this action is a crucial step in maintaining the stability and synchronisation of a power system. There is some uncertainty regarding the success and speed of protection actions, and, as a consequence, power flows may need to be restricted and more expensive, or higher carbon, power sources utilised. We argue that choices both in operation and the design of the system and its control can be facilitated by parameter investigations of power system stability models such as the swing equation and applying quick but effective stability metrics to illustrate the effect of a given parameter value change. A rigorous study of the strategies available to the system operator could be provided in-part by the continuous variation of model parameters, which could possibly uncover optimal parameter values to maximise stability at the design stage or on-line. The analytic CCT metric derived in this paper is able to capture sensitivities in stability of a given fault as a network parameter is varied. In particular, this paper studies the effect of a load parameter on the stability of a given fault in an aggregated network.

The rest of this paper is organised as follows: In Section II we formulate a CCT estimate denoted (where the subscript ‘’ signifies ‘Hamiltonian’) using direct methods and introduce the aggregate network used to conduct our investigations, the two-machine infinite bus (TMIB) network. By considering polynomial approximations of we derive an analytic CCT metric denoted (where the subscript ‘’ signifies ‘analytic’) in Section III. A parametric investigation of the effect on stability of different loadings on an aggregated network given a particular fault is presented in Section IV and finally, conclusions are drawn in Section V together with suggestions for future work.

## Ii Fault analysis using energy functions

### Ii-a Model description

We consider the classic swing equation model [2, 1] to describe the stability effects of transient faults on a power system with synchronous generation. The generators are modelled as voltage sources behind reactances and the loads on the network are of constant impedance. In general, generators have small losses due to damping [19] so without loss of generality we assume zero damping for generators. This model can be written as a set of coupled one-dimensional ordinary differential equations (ODEs), which describe the dynamics of the rotor angles of each synchronous generator in a network by considering Newton’s second law of dynamics. In vector form the equation is

(1) |

where

and

The vectors and are the generator rotor angles and angular speeds respectively, and the elements of the vector function are

where is a lumped parameter, (where is the grid frequency: in Europe), is the inertia constant, is the mechanical power input and is the electrical power output.

The loads on the power system are assumed to be constant impedance loads such that Kron reduction [20] can be applied to the network. Therefore, the swing equations describe the dynamics of a reduced network comprising of constant voltage sources connected through a network of impedances [2]. The total power consumed by conductive loads at generator is given by

(2) |

where is the internal voltage of generator ( assumed constant), is the conductance between generators and and is the shunt conductance at bus . The total electric power leaving generator is

(3) |

where is the maximum active power flow between generators and and is the susceptance of the network connection between node and node . The admittances are the elements of the (symmetric) reduced bus admittance matrix . Kron reduction is fundamentally a matrix operation permitted by applying Kirchoff’s laws to a power network and constructing from a larger bus admittance matrix where . The bus matrix is a block matrix which contains the full topology and load distribution (including the synchronous reactance) of a power network with synchronous generators.

### Ii-B Fault analysis

#### Ii-B1 Stability analysis of transient faults

The objective of transient fault analysis is to investigate whether a system will remain stable once a fault has been cleared and, ideally, no further action from the system operator would be required. We assume, without loss of generality, that the moment a power system suffers a short-circuit is at time and the fault is cleared at time . These two points in time define three distinct regimes in order to analyse the dynamics of a fault on a power system. These are (pre-fault), (fault-on) and (post-fault).

The fault analysis method in [19, Chapter 2] (recently summarised in [21]) for power networks with constant impedance loads, is employed in this paper. Each regime has a different bus matrix (and therefore reduced admittance matrix ) which will change the values of the parameters , and in the vector function (1). Therefore, three separate sets of equations of the form (1) are required to model the power system for all time given by

(5) |

where the labels ‘’, ‘’ and ‘’ refer to the parameter values for the system in regimes , and respectively. Pre-fault, a power system is assumed to be balanced and therefore we assume that (5) is located at a stable (‘s’) equilibrium point

for where . The dynamics for are given by

(6) |

during the fault and

(7) |

after the fault. From these expressions we can define the CCT, denoted , formally as the maximum value of such that in the post-fault trajectory (7) there is one full swing of the rotor angles before some pairs of rotors angles begin to diverge [19]; this is also known as first swing stability and is generally found algorithmically using power systems software packages.

#### Ii-B2 A CCT approximation using energetic methods

In general, a conservative metric for the local stability of systems of the form (1) can be found by constructing a suitable Lyapunov function. Direct methods use so-called energy functions [5], which can also serve as Lyapunov functions, to measure the global stability of such systems. A stability boundary is constructed in terms of a critical system energy in the post-fault regime and a power system is classified as unstable when the total system energy surpasses this critical energy.

The total system energy can be measured when a power system is modelled as a Hamiltonian system. However, the power consumed by the loads is a path-dependent quantity [5] and cannot be modelled exactly by a conservative system. The survey paper [22] collects numerous attempts that have been used to approximate this term so that an appropriate Hamiltonian system can be used. The most accepted technique [2, p. 231] models the power consumed by the loads as a constant term given by

(8) |

where the point is a stable stationary point in the post-fault regime which solves with . The dynamics of a power system with assumption (8) employed can be written as

(9) |

where terms depending on the conductive parts of loads are isolated to obtain the vector function . This function has a similar structure as in (1) where

and the elements of the vector are given by

with

The Hamiltonian function

(10) |

quantifies the post-fault system energy for a system of the form (9) and is the sum of the kinetic energy and the potential energy for a power system with generators where

(11) |

and

(12) |

An approximation of the CCT, denoted can be found by integrating the dynamics of the system during a fault until the system energy reaches the critical boundary (which will be computed later). More specifically, such an estimate can be obtained by observing the first instance that the Hamiltonian

(13) |

for where, in general, the energy difference

(14) |

is positive for a suitably chosen post-fault network. Note that, the power system during the fault is not modelled as a Hamiltonian. The CCT approximation is much faster to compute than the traditional CCT because the dynamics of the post-fault system (7) do not need to be computed.

The so-called closest UEP (unstable equilibrium point) method [11] is used to find the critical system energy in this work because, although it is the most conservative method (compared to the controlling UEP method or the potential energy boundary surface method [11]) it can be applied to any power system without considering the specific fault that a system suffers. In the presence of large linear loads in the network, the use of direct methods might lead to overestimates of the actual stability boundary [23] however, the intention is to study the effect of stability trends, so the closest UEP serves as an adequate method to capture the system energy for initial parametric studies.

### Ii-C An aggregate network

In order to study trends in stability under parametric variations, the dynamics of each generator in a network can be grouped into synchronous regions according to the electrical distance between individual generators. Previous studies [24, 25, 26] have used an aggregate power network model to study the global dynamics of a power system. Typically, the models presented in these references study the GB power network using a three bus network with the inertia of one machine at least two order of magnitudes larger than the other two. These models lend themselves well to be studied using a so-called two machine infinite bus (TMIB) system. This network structure has been previously studied in [27, 11, 28, 29]. The ODE in the form (9) for this system is given by

(17) | ||||

where we have employed assumption (8) to get a conservative system. After Kron reduction there are three interconnected buses in the network: two buses connected to synchronous generators and an infinite bus (bus 3). The infinite bus models the dynamics of a large section of a network as a generator with infinite inertia and constant internal voltage . As such, is a constant and without loss of generality we can set and use it as a reference point for the other two rotor angles.

## Iii An analytic stability metric

An analytic stability metric, denoted as , which is purely a function of network parameters, is presented and derived here. The metric is formulated by considering (13) and approximating both the Hamiltonian (10) and the fault trajectory (6) in this equation as polynomial functions of rotor angles and time respectively. During a fault, it is assumed that the governor control systems for the mechanical input power in each generator are not able to act quickly enough to change the parameter value during (or immediately after) the fault; therefore throughout this analysis. In addition, the dynamics of the rotor angles during the fault are approximated as constant but different accelerations. Although some of these approximations may seem cumbersome, and will detract significantly from the true dynamics of the system, they are valid in the limit as the CCT tends to zero. Therefore, we assume that for small values of the CCT these approximations can be assumed to capture the dynamics of a power system modelled as a Hamiltonian system. In addition, we remind the reader that this metric is designed to provide an instant illustration of the stability of a power system under the variation of a chosen parameter. In the next section, we demonstrate how our analytic metric can be used to find approximate regions for values of a load parameter in an aggregate network that improves stability for a given fault on the network.

### Metric formulation

An analytical expression to approximate the CCT for three-phase to ground faults close to a given bus (on a balanced system such that it can be modelled by means of a single phase equivalent) can be found by adapting the energetic framework for the CCT presented in Section II-B2. The expression (13) is altered such that we solve

(19) |

where is a polynomial function such that

(20) |

with initial condition

(21) |

and .

In order to construct the function , we first approximate the Hamiltonian function as a polynomial function of the rotor angles, denoted . The kinetic term from the post-fault Hamiltonian function is removed by also modelling the dynamics during the fault as a Hamiltonian system. In general, there is no stable stationary point available during the fault so the power consumed by the conductive loads is approximated as a constant such that the dynamics can be written in the form (9) to give

(22) |

for and . The Hamiltonian during the fault is given by

(23) |

In accordance with conservative systems, is a constant in time and this property is used to recast the expression for in (13) by considering the trivial relation

(24) |

resulting in the succinct expression

(25) | ||||

which has no dependence on rotor speeds. The values of the parameters , , and are found from the fault-on and post-fault reduced admittance matrices and the internal voltages.

A candidate function for can be found by replacing the cosine terms in (25) with the function

for small . This substitution gives

(26) | ||||

where the constant is found by applying the initial condition (21), i.e. to (26). Therefore,

(27) | ||||

where and (26) can be re-written as

(28) | ||||

where the constant is written explicitly.

In order to make (28) an explicit polynomial function of time, the fault trajectory must also be written as a polynomial function of time. In general, the dynamics during a fault are non-trivial [30] but in order for (19) to be analytically solvable for time, the rotor angle dynamics in (28) must have the form

(29) |

where the initial condition holds for all . An appropriate value for the acceleration can be found by assuming that for small CCTs the rotor dynamics can be modelled as a constant acceleration equal to the initial rotor acceleration at . This is given by

(30) |

for short fault times. From equation (30) the rotor accelerations . By substituting expressions (29) for the rotor angles into (28), the function

(31) | ||||

is a quadratic in where . Now (19) can be written as

(32) |

where the coefficients

are functions of the power network parameters. The solution of (32) and thus the expression for our analytic CCT is given by

(33) |

The smallest real value of is taken for a given set of parameters. A purely imaginary value for is produced if the discriminant or if and . In the case where and two positive roots are produced, otherwise there is one real root to (33). However, in general the parameter is positive because the total electrical load of a network reduces during a fault and so it is reasonable to assume that .

## Iv Parametric stability analysis

### Iv-a Implementation details

The stability of a power network is not only dependent on the type or duration of a fault but also on the choice of system parameters. Optimal regions of parameter space that increase the stability of a power system can be identified by the variation of system parameters. Here, we investigate values for a load on a TMIB network which improve system stability for a given fault, by comparing the metrics and outlined in Sections II and III against the true CCT . The energy difference is also compared against the temporal stability metrics.

The variation of a load parameter in a power network will change one of the elements in the full bus admittance matrix , but due to Kron reduction all the parameters in the reduced admittance matrix in each regime will change. Therefore, for each incremental change in the parameter we consider, a new fault study is required to find the CCT. In each fault study, the mechanical input powers for each generator are found by performing a pre-fault power flow for the system. The power flow is conducted using the same (small) rotor angles found in [19] for each incremental change in the parameter value to ensure that the network is initially in a stable state.

It is relatively quick to conduct a single fault study to find the true CCT . However, for each incremental change in a parameter, a new fault analysis is required to find the CCT. An analytic CCT has been developed to study trends in stability of power systems, which can be found instantly once the relevant parameter values for the model in (5) have been collected. For a given fault, system parameters can be varied continuously using an analytic CCT and this can provide an initial picture of stability that informs more detailed analysis.

All the stability metrics introduced in this paper, except for the true CCT , are dependent on a critical energy boundary which is dependent on the location of the closest UEP in this work. The position of the closest UEP will change under the variation of the loads and there are techniques developed in the literature to find these quickly [31]. However, we choose to use numerical continuation (previously applied to power systems in [32]) to illustrate interesting features of the closest UEP under the variation of loads. The stationary points of a TMIB system, modelled by the ODEs in (17), are located using the continuation software AUTO^{1}^{1}1http://indy.cs.concordia.ca/auto/ as a load parameter is varied. There is no rigorous proof provided in this paper that all the possible unstable equilibria on the stability boundary of a stable equilibrium point can be found from the solution branches from numerical continuation. However, no other solutions were found for this system when performing an exhaustive search over state space using the root finding algorithm from the Scipy^{2}^{2}2http://docs.scipy.org/doc/ library in the Python^{3}^{3}3www.python.org programming language. Therefore, without further analysis, it is assumed that only in a TMIB system can all the necessary stationary points be found using this continuation method.

The stationary points for each value of the continuation parameters in Fig. 4a and Fig. 6a are obtained by the following method: A stable stationary point denoted by that solves the post fault equation (where and ) is found using the root finding algorithm , where . This point belongs to the lower branch (blue squares) of the bifurcation diagram in Figs 4a and 6a. The other (unstable) equilibria on the boundary of the stability region of the stable equilibrium point which satisfy are found by numerical continuation of the element from the reduced admittance matrix . Once the continuation branches are found, the stationary points at the value of found in are recorded. The local stability of the stationary points obtained are found by computing the eigenvalues of the Jacobian matrix for the system (17). The stability of the solution branches in Figs 4a and 6a are stated in terms of the number of eigenvalues with real part greater than zero which can be found in the figure caption of Fig 4.

In this study a 9 bus, 3 generator power network found in [19] is used, a schematic of this network is provided in Fig. 3. All parameter values for this network are taken from [19]. This network is adapted into a TMIB network by changing one of the generators, which has an inertia an order of magnitude larger than the other two, into an infinite bus. The specific fault we consider is a three-phase to ground fault close to bus 7 on the line connecting buses 5 and 7. The post fault network is identical to the pre-fault network except the line connecting buses 5 and 7 is switched out.

### Iv-B Results

The size and nature of loads on actual power systems can vary over time and, in respect of the susceptive part, can be modified by the addition of reactive compensation. As a consequence, the conductive and susceptive parts of load C (denoted and respectively) of the network in Fig. 3 are investigated by varying one part while maintaining the other constant at its original value. In Figs. 4a and 6a the domains of the parameters and respectively are constrained by two conditions: the energy margin and that the synchronous machines are operating as generators in the pre-fault power flow, i.e. and . There is an additional constraint in Fig. 4 where only positive values of conductance are explored.

The critical energy change for the system (black line), the CCT estimate (blue line) and the analytic CCT (red line) are plotted as functions of the continuation parameters in the lower panels of Figs. 4 and 6. In addition, the true CCT (green line) is plotted using a simple algorithm that uses a binary search to find the maximum duration which the fault can be left on-line such that the rotor angles have one full swing together before they diverge. There are two different scales to facilitate observing the functions in the lower panels of Figs. 4 and 6: the energy change should be read using the right-hand y-axis labels and the three time metrics should be read using the left-hand y-axis labels as indicated in the figures.

In Fig. 4 the dependence of the system stability, for the fault we consider, on the conductance is studied as the susceptance is held constant. In Fig. 4a there is a discontinuity in the closest UEP (thick line) at (vertical dotted line) which is located between two pairs of fold points at and . (See [33] for an explanation of fold points). In Fig. 4b there is a discontinuous change in the gradient of which coincides with the discontinuity at , but the maximum point for at does not coincide with the other discontinuity in the closest UEP, nor any other points of significance in Fig. 4a. The analytic CCT is observed to approximate the CCT estimate very well as the load parameter is varied. However, the energetic techniques used to find and have resulted in non-conservative estimates of the true CCT . This feature is a manifestation of the original issue with direct methods which concerns the dissipative term (2) at each bus in the reduced network. Direct methods can be used for conservative stability assessments where the transfer conductances in the reduced network matrix are assumed to be small or zero [11], therefore the CCT estimate and the analytic CCT are strict lower bounds of the true CCT for networks with zero transfer conductances. However, even for a network with lossless lines, Kron reduction invokes complications in which a shunt load conductance in the full bus admittance matrix will increase the absolute values of in the reduced admittance matrix [22, 12].

In Fig. 7, the susceptance is varied in a network with small load conductances and it is observed that and are lower bounds to when compared to the results in Fig. 6. Despite whether the analytic CCT is an over or underestimate of the true CCT, it performs well as an indicator of the expected increase or decrease of the CCT as is changed. The greatest CCT as measured by all metrics for the fault we have considered is produced at . (This behaviour was found for all possible faults on the network under the variation of one of the loads , or within an order of magnitude of its nominal value.) In general, a lower mechanical input power from each generator is required for lower network loadings and therefore the acceleration of the generator rotor angles is roughly proportional to the mechanical input power, assuming that the load of the network decreases during a fault. Therefore, there is more time for the rotors to reach a critical value where they begin to diverge.

In Fig. 5 the change to the energy boundary (plotted in the manifold where ) due to the discontinuity in the location of the closest UEP is illustrated in Fig. 4a. In this manifold, the energy boundary is plotted as the level set (black line) for conductance values , and where the discontinuity occurs at . In each sub-figure of Fig. 5 the stationary points of (17) are plotted using the same marker style found in the legend of Fig. 4a and the level set is observed to intersect the closest UEP.

In Fig. 6 the dependence of the system stability, for the fault we consider, on the susceptance is studied as the conductance is held constant. In Fig. 6a there is a discontinuity in the closest UEP (thick line) at which is located between two fold points at and . In Fig. 6b there is a discontinuous change in the gradient of , and at the discontinuity in the closest UEP. The maximum of occurs at the highest value of susceptance plotted, which shows that the energy margin is not the best metric to quantify stability. The analytic CCT is very close to the CCT estimate as is varied and are, again, overestimates due to the presence of non-negligible transfer conductances. The maximum points of and , both at are very close to the closest UEP discontinuity, however the maximum point of the true CCT is lower, at . Despite this, the change in the true CCT as the susceptance is varied is well captured by the CCT approximations, except in the region where the gradients are of different signs.

Optimum susceptance | ||
---|---|---|

s | N/A | |

s | s |

Optimum susceptance | ||
---|---|---|

s | N/A | |

s | s |

It is observed that the system stability can benefit by setting the susceptance of load C to the optimum susceptance as measured by the analytic CCT because these susceptance values can be evaluated without the need for numerical integration. From Table II, the true CCT using the parameter values as stated in [19] is . A network operating at the optimum value of susceptance would give a true CCT of and this is an increase of . However, the value of the true CCT at the optimum value of susceptance as measured by the analytic CCT is which is a smaller but significant increase of . The advantage of using the optimum values of susceptance as measured by the analytic CCT is that an improved susceptance value is known as soon as the relevant network parameters have been collected.

The continuation of the susceptive part of load B (with the other loads at original values) is considered for the same fault at bus and the results are given in Table I. The results for load are not included in the tables because there was no maximum point for CCT found as the susceptive part of load A was varied and the trends in stability were similar to the lower panel of Fig 4. The largest CCT was found for which is the lowest value of susceptance for which the mechanical input powers of the synchronous generators and are both greater than zero.

## V Discussion

In this paper we have presented a new analytic CCT metric designed to be able to capture trends in the true CCT as a system parameter is varied. Specifically the effects of the conductive and susceptive parts of a load parameter on the network were considered, given a fault on the network. The analytic CCT metric was formulated by taking a polynomial approximation of the CCT estimate developed from direct methods and it is found that despite the simplicity of the formula, the analytic CCT is a good approximation to for short times. Given the difficulty for all direct methods to incorporate power networks with non-negligible transfer conductances, it was expected that the two approximating metrics and would not perform well as good estimates to the true CCT . However, the two approximating metrics performed much better as indicators of trends in stability as a load parameter was varied. In addition, the results in this paper were generated for an aggregate power network, the two-machine infinite-bus, where studies on aggregated networks [26, 24] are generally conducted to analyse global trends of a network with a much larger number of generators, buses and loads. The approximating CCT metrics are valid in principle for power systems with or without an infinite bus and numerical results will be extended to alternative networks without an infinite bus in future work.

Direct methods were chosen to formulate the CCT approximations because of their ability to construct a well-defined stability boundary in terms of a critical energy . However, this stability boundary is dependent on finding a critical UEP of the system that can be used to approximate the energy when a power system is modelled as a Hamiltonian system. In this paper, the closest UEP was chosen to compute the energy boundary because it is valid for any fault on a power system. A more accurate method to quantify the system energy uses the controlling UEP [34], which is dependent on the fault a network suffers and can be found using the ‘boundary of stability region-based controlling unstable equilibrium point’ (BCU) method [11, 35]; the limitations of which are discussed in [27]. For the TMIB network it was possible to find all UEPs on the stability boundary of an appropriate stable equilibrium point of the system under the variation of a load parameter using numerical continuation methods. There was no attempt to formally prove that all UEPs are captured, but an exhaustive algorithmic search was used to confirm this. However, the scalability of this approach is limited because the number of equilibria increases as the system size increases [36] and it is increasingly difficult to identify the critical UEP and this is a significant area of research in itself [31, 37]. In addition, the positions of the equilibria have to be found for each incremental change of a load parameter. This is another reason the analysis in this paper was limited to study trends on an aggregated power network with a small number of generators.

The most general use for the analytic CCT proposed in this paper is to capture stability trends under the variation of a network or generator parameter and we have specifically studied the effect of varying loads. A more specific suggested use for the metric is to locate regions of parameter space that will improve the system stability in terms of CCT. For the fault studied in our analysis of the TMIB network, we found that the optimum value of CCT for the conductive part of load C is for non-negative conductance values and that the CCT decreases as increases. One interpretation of this result is that a high penetration of linear power sources, local to the point of power consumption has the effect of increasing system stability, but research in [4] suggests that the effect on system dynamics is dependent on the specific technologies in the generator. More promising results were found for the variation of the susceptive part of a load under constant conductance. The variation of susceptive loads can represent, for example, a network owner’s installation of reactive compensation, a measure that is known to contribute not only to voltage regulation but also transient stability [38]. An optimal value of susceptance that maximises the CCT was identified in all three temporal metrics, however the optimal value of susceptance that maximises was different to the value that maximises both and . It was found that true CCT would be improved by if the optimum susceptance loading for load C as measured by the analytic CCT was used instead of the original value of susceptance for this load. Given the quick assessment provided by an analytic CCT, an exhaustive study of all three-phase to ground faults on a network with their respective post-fault clearing strategies can be performed under a continuous range of loading conditions without having to do any formal fault study. Future research would then be required to test mathematical optimisation methods designed to find optimal loading distributions that improve the stability of the power system as measured by the CCT.

Despite its drawbacks, our analytic stability metric has the potential to inform optimal fault management strategies to improve system stability through parameteric investigation. Its key advantage is that it can be computed instantly once all the system parameter values for the pre-fault, fault-on and post-fault systems have been collected. This feature of stability metrics could be of use due to system dynamics becoming more unpredictable from the changing nature of loads [3] and generation [39, 40, 41] under the constraint of limited power flow through transmission lines. Particularly, investigations on the effect of low inertia on power system stability [25, 42] can potentially benefit and this is an area of future work. Furthermore, optimisation techniques could be applied to analytic metrics to find regions in parameter space that increase power system stability in terms of CCT.

## References

- [1] P. S. Kundur, Power System Stability. McGraw-Hill, 1994.
- [2] J. Machowski, J. W. Bialek, and J. R. Bumby, Power System Dynamics: Stability and Control. Wiley, 2008.
- [3] K. Yamashita, S. Djokic, J. Matevosyan, F. Resende, L. Korunovic, Z. Dong, and J. Milanovic, “Modelling and aggregation of loads in flexible power networks - Scope and status of the work of CIGRE WG C4. 605,” in Power Plants and Power Systems Control, vol. 8, no. 1, 2012, pp. 405–410.
- [4] J. Slootweg and W. Kling, “Impacts of distributed generation on power system transient stability,” IEEE Power Engineering Society Summer Meeting, vol. 2, pp. 862–867, 2002.
- [5] M. A. Pai, Energy Function Analysis For Power System Stability. Kluwer, 1989.
- [6] F. Dörfler and F. Bullo, “Synchronization and transient stability in power networks and non-uniform Kuramoto oscillators,” SIAM J. Control Optim., vol. 50, no. 3, pp. 1616–1642, 2012.
- [7] Y. Susuki, I. Mezić, and T. Hikihara, “Coherent swing instability of power grids,” Journal of Nonlinear Science, vol. 21, no. 3, pp. 403–439, Feb. 2011.
- [8] V. Ajjarapu and B. Lee, “Bifurcation theory and its application to nonlinear dynamical phenomena in an electrical power system,” IEEE Transactions on Power Systems, vol. 7, no. 1, pp. 424–431, 1992.
- [9] M. Galaz, R. Ortega, A. S. Bazanella, and A. M. Stankovic, “An energy-shaping approach to the design of excitation control of synchronous generators,” Automatica, vol. 39, pp. 111–119, 2003.
- [10] Y. Hasegawa and Y. Ueda, “Global basin structure of attraction of two degrees of freedom swing equation system,” International Journal of Bifurcation and Chaos, vol. 9, no. 8, pp. 1549–1569, 1999.
- [11] H.-D. Chiang, Direct Methods for Stability Analysis of Electric Power Systems: Theoretical Foundation, BCU Methodologies, and Applications. Wiley, 2011.
- [12] V. R. Sastry, “Validity of neglecting transfer conductances in transient-stability studies,” Proceedings of the Institution of Electrical Engineers, vol. 120, no. 12, p. 1539, 1973.
- [13] M. J. Laufenberg and M. A. Pai, “A new approach to dynamic security assessment using trajectory sensitivites,” in 20th International Conference on Power Industry Computer Applications, 1997, pp. 272–277.
- [14] A. A. Fouad and S. E. Stanton, “Transient stability of a multi-machine power system Part I: Investigation of system trajectories,” IEEE Transactions on Power Apparatus and Systems, vol. 100, no. 7, pp. 3408–3416, 1981.
- [15] I. A. Hiskens and M. A. Pai, “Trajectory sensitivity analysis of hybrid systems,” IEEE Transactions on Circuits and Systems I: Fundamental Theory and Applications, vol. 47, no. 2, pp. 204–220, 2000.
- [16] T. B. Nguyen and M. Pai, “Dynamic security-constrained rescheduling of power systems using trajectory sensitivities,” IEEE Transactions on Power Systems, vol. 18, no. 2, pp. 848–854, 2003.
- [17] A. Grilo, A. Mota, L. Mota, and W. Freitas, “An analytical method for analysis of large-disturbance stability of induction generators,” IEEE Transactions on Power Systems, vol. 22, no. 4, pp. 1861–1869, Nov 2007.
- [18] H.-D. Chiang and C.-C. Chu, “Theoretical foundation of the BCU method for direct stability analysis of network-reduction power system models with small transfer conductances,” IEEE Transactions on Circuits and Systems I: Fundamental Theory and Applications, vol. 42, no. 5, pp. 252–265, 1995.
- [19] P. M. Anderson and A. A. Fouad, Power System Control and Stability, 2nd ed. IEEE Press, 2002.
- [20] F. Dörfler and F. Bullo, “Kron reduction of graphs with applications to electrical networks,” IEEE Transactions on Circuits and Systems I: Regular Papers, vol. 60, no. 1, pp. 150–163, Jan. 2013.
- [21] A. Gajduk, M. Todorovski, and L. Kocarev, “Stability of power grids: An overview,” The European Physical Journal Special Topics, vol. 223, no. 12, pp. 2387–2409, Jun. 2014.
- [22] M. Ribbens-Pavella and F. Evans, “Direct methods for studying dynamics of large-scale electric power systems: A survey,” Automatica, vol. 21, no. 1, pp. 1–21, Jan. 1985.
- [23] T. Athay, R. Podmore, and S. Virmani, “A practical method for the direct analysis of transient stability,” IEEE Transactions on Power Apparatus and Systems, vol. 98, no. 2, pp. 573–584, 1979.
- [24] F. M. Hughes, O. Anaya-Lara, N. Jenkins, and G. Strbac, “A power system stabilizer for DFIG-based wind generation,” IEEE Transactions on Power Systems, vol. 21, no. 2, pp. 763–772, May 2006.
- [25] A. Ulbig, T. S. Borsche, and G. Andersson, “Impact of low rotational inertia on power system stability and operation,” in Proceedings of the 19th IFAC World Congress, 2014.
- [26] K. Johnstone, R. M. Tumilty, K. R. W. Bell, and C. D. Booth, “Transient stability assessment of the GB transmission system with high penetrations of wind power,” 13th International Workshop on Large-scale Integration of Wind Power into Power Systems as well as on Transmission Networks for Offshore Wind Power Plants, 2014.
- [27] A. Llamas, J. De La Ree Lopez, L. Mili, A. G. Phadke, and J. S. Thorp, “Clarifications of the BCU method for transient stability analysis,” IEEE Transactions on Power Systems, vol. 10, no. 1, pp. 210–219, 1995.
- [28] M. Anghel, F. Milano, and A. Papachristodoulou, “Algorithmic construction of Lyapunov functions for power system stability analysis,” Circuits and Systems I: Regular Papers, IEEE Transactions on, vol. 60, no. 9, pp. 2533–2546, 2013.
- [29] Y. Guo, D. J. Hill, and Y. Wang, “Nonlinear decentralized control of large-scale power systems,” Automatica, vol. 36, no. 9, pp. 1275–1289, 2000.
- [30] C. Chia-Chi, “Towards a theory of multi-swing transient instability problems in electric power systems,” IEICE Transactions on Fundamentals of Electronics, Communications and Computer Sciences, vol. E88-A, no. 10, pp. 2692–2695, 2005.
- [31] C. W. Liu and J. S. Thorp, “A novel method to compute the closest unstable equilibrium point for transient stability region estimate in power systems,” IEEE Transactions on Circuits and Systems I: Fundamental Theory and Applications, vol. 44, no. 7, pp. 630–635, 1997.
- [32] L. Chen, Y. Min, F. Xu, and K.-P. Wang, “A continuation-based method to compute the relevant unstable equilibrium points for power system transient stability analysis,” IEEE Transactions on Power Systems, vol. 24, no. 1, pp. 165–172, Feb 2009.
- [33] S. Strogatz, Nonlinear Dynamics and Chaos: With Applications to Physics, Biology, Chemistry and Engineering. Westview Press, 2001.
- [34] H.-D. Chiang, “A theory-based controlling UEP method for direct analysis of power system transient stability,” in IEEE International Symposium on Circuits and Systems, vol. 3, May 1989, pp. 1980–1983.
- [35] H.-D. Chiang, F. F. Wu, and P. P. Varaiya, “A BCU method for direct analysis of power system transient stability,” IEEE Transactions on Power Systems, vol. 9, no. 3, pp. 1194–1208, 1994.
- [36] C. J. Tavora and O. J. M. Smith, “Equilibrium analysis of power systems,” Transactions on Power Appartatus and Systems, vol. 91, no. 3, pp. 1131–1137, 1972.
- [37] J. Lee, “Dynamic gradient approaches to compute the closest unstable equilibrium point for stability region estimate and their computational limitations,” IEEE Transactions on Automatic Control, vol. 48, no. 2, pp. 321–324, Feb. 2003.
- [38] B. Delfino, G. Denegri, M. Invernizzi, and P. Pinceti, “Estimating first swing stability of synchronous machines as affected by saturation controls,” IEEE Transactions on Energy Conversion, vol. 3, no. 3, pp. 636–646, 1988.
- [39] H. Urdal, R. Ierna, J. Zhu, C. Ivanov, A. Dahresobh, and D. Rostom, “System strength considerations in a converter dominated power system,” IET Renewable Power Generation, 2014.
- [40] T. B. Nguyen, M. A. Pai, and E. Muljadi, “Impact of wind power plants on voltage and transient stability of power systems,” in IEEE Energy2030, no. November, 2008.
- [41] D. Gautam, V. Vittal, and T. Harbour, “Impact of increased penetration of DFIG-based wind turbine generators on transient and small signal stability of power systems,” IEEE Transactions on Power Systems, vol. 24, no. 3, pp. 1426–1434, 2009.
- [42] P. Tielens and D. V. Hertem, “Grid inertia and frequency control in power systems with high penetration of renewables,” in Young Researchers Symposium in Electrical Power Engineering, no. 2, 2012, pp. 1–6.