On-Line Method for Determining Power System Transient Stability

Electrical power systems are continually experiencing disturbances. Of primary concern to the power industry are those disturbances classified as event disturbances, which include generator outages, short-circuits caused by lightning, sudden large load changes, or a combination of such events. It is of critical importance to be able to predict the ability of the power system to reach an acceptable steady-state operating environment after an event disturbance. This is the area of concern in system transient stability analysis.

The present invention relates to methods for analyzing power systems and, more particularly, to an on-line method for determining power system transient stability.

An event disturbance usually leads to a change in the configuration of the power system. Load disturbances are the small random fluctuations in load demands. The system usually remains unchanged after load disturbances.

In a power system the transient stability of concern is the synchronization of the phase and rotational velocity of the generators.

Power system operation is divided into three operating phases or operational configurations: pre-fault, fault-on, and post-fault. Transient stability is the study of whether the post-fault system will converge (tend) toward an acceptable steady state as time increases.

Consider the following possible scenario. An event disturbance, such as lightning, causes a short circuit in a power system. The surge in current demand caused by the short circuit causes a protective circuit breaker to trip, disconnecting the short circuit from the system. The sudden change in the load upon the generators causes at least one of them to overspeed, increasing the voltage on the system. If the post-fault condition becomes stable, the overspeeding generator will be gradually brought down to the same speed as the other generators in the system. If the post-fault conditioning will not stabilize, however, other circuit breakers could open in a chain reaction. In one worst case scenario, a significant amount of the power system generators may go off-line to protect themselves. The subsequent reduction in the service to utility customers, and the time required to reconnect the various generators and bring them back on-line, could be very costly and inconvenient. At present there is no convenient and reliable way of predicting when such an unstable condition is imminent.

Transient stability analysis is routinely performed by a utility company during system planning. The industry standard for transient stability varies from region to region, but in general requires the ability of the system to withstand severe event disturbances, such as a three-phase fault condition in close proximity to a generator with a stuck breaker.

Until recently transient stability analysis has been performed by utilities exclusively by means of the numerical integration of nonlinear differential equations describing the fault-on condition and the post-fault condition. Computer programs that can handle 2000 buses and 300 generators with detailed models are available, but may take several hours of computer time to run. Since the time interval of interest in transient stability analysis is between a few milliseconds and a few seconds, simulation for on-line security analysis which require a fast and reliable method has not been feasible, until this invention.

Mathematically the problem of transient stability in accordance with this invention can be defined by three sets of differential equations, one each for the pre-fault, fault-on, and the post-fault conditions:

where t.sub.f is the time of fault occurrence and t.sub.p the time of fault clearance.

Since the power system varies with time a value can be calculated that is indicative of the state of a system at any given moment, which in turn can be used to determine if the system will be stable at a future point in time.

The state or condition of a system is represented by a set of numbers. These values, combined with input functions, are introduced into the equations describing the dynamics of the power system. The solution to the equations provides a prediction of the future state and output of the system. The state or condition of a system is expressed as a function of time.

All the states of a system are said to exist in state space. The path the state or condition of a system follows as time progresses is known as a "trajectory".

One of the most important determinations is whether the post-fault trajectory of the system will approach a stable equilibrium point. This will determine transient stability. Around each stable equilibrium point (SEP) is a region of stability. Any point within this stability region will approach the SEP with respect to time.

At the time a fault occurs, the system is at an SEP. While the fault is occurring, the system begins a fault trajectory out of the region of stability. If the fault is cleared before the trajectory leaves the stability region, then the post-fault condition will be stable. The time between the fault occurrence and the time when the post-fault trajectory leaves the stability region is defined by t.sub.c, the critical time.

For a power system having N generators, the classical model of the dynamics of the i.sup.th generator is governed by the equations: ##EQU1## Equations 2-a and 2-b represent rotor mechanics; .delta..sub.i, w.sub.i are the respective angle and speed deviations (from a synchronous reference) of the ith rotor; M.sub.i and D.sub.i are respective inertia and damping constants of the ith generation; P.sub.i.sup.m is the constant mechanical input power; and P.sub.i.sup.e is determined by the electrical network. The i.sup.th generator is configured as a voltage source of constant magnitude E.sub.i and phase .sigma..sub.i driving (through a transient reactance) an electric network consisting of transmission lines and loads, configured as constant impedances. Viewed as an N port driven by the voltage sources E.sub.i <.delta..sub.i, this network is described by the reduced complex admittance matrix Y with coefficients Y.sup.i <.phi..sub.ij (see FIG. 1). The real, positive diagonal terms G.sub.ij of Y are separated as indicated. Equation 2-c gives the real power delivered to the electric network by the ith generator. In this N generator system, node N+1 is an infinite bus, with E.sub.n+1 =0 and .delta..sub.n+i =0. Thus the voltage at node n+1 serves as a synchronous reference.

Equations 2a, 2b, 2c represent the post-fault condition. Suppose there is a set of (i elements) stable equilibrium points for the system described by 2a, 2b, 2c that is pairs of .delta..sub.i, w.sub.i, where w.sub.i is zero. Then, if it is assumed for purposes of this description that the transfer conductances are zero, a version of equations 2a, 2b, 2c can be derived where the stability region can be estimated: ##EQU2##

Using Lyapunov theory, a direct method of determining the transient stability of the post-fault condition can be deduced.

Let (.delta..sub.s, w.sub.s) be a stable post-fault equilibrium. Define the energy function V(.delta., w) of the state: ##EQU3##

At the post-fault equilibrium, V(.delta..sub.s, w.sub.s)=0. Differentiating V(.delta..sub.s,w.sub.s) along the trajectories of equation 3 gives ##EQU4## This inequality shows that V qualifies for use as a Lyapunov function since V(.delta.(t),w(t)) must decline along the post-fault trajectory.

Lyapunov theory then leads to the following approach to finding if the post-fault condition is converging toward stability:

Let v>0 be a constant such that on the boundary of the stability region one has V(.delta.,w)>v. Then the connected component containing (.delta..sub.s,w.sub.s) of the set R(v)={.delta.,w)<v} is inside the stability region.

Three observations can be made. First, Theorem 1 is not an immediate consequence of general Lyapunov theory, since V(.delta.,w) does not become unbounded as the absolute value of .delta. goes to infinity.

Second, Theorem 1 provides an estimate of the stability region. To get the best estimate (the largest R(v)), one must choose v as large as possible. The largest value V.sub.cr is V.sub.cr =min{V(.delta.,w)l(.delta.,w) that are elements of the set of stability points on the boundary of the stability region). It can be deduced that V.sub.cr =min UEP {V(.delta..sub.e,w.sub.e)} where the minimization is over all the unstable equilibrium points (UEPs) (.delta..sub.e,w.sub.e) that are elements of the set of stability points on the stability region boundary. The UEP that achieves this minimum is called the critical UEP.

Third, R(V.sub.cr) is an underestimate of the stability region; see FIG. 2.

Theorem 1 suggests a stability criterion, according to the following three step procedure:

Step 1: Calculate the UEP lying on the boundary of the stability region. Calculate V.sub.cr =min UEP {v(.delta..sub.e,w.sub.e)}

Step 2: Evaluate v along the faulted trajectory (.delta.(.delta.(.tau.), w(.tau.). If this v is less than V.sub.cr then .tau. is less than .tau..sub.c (critical time ) and the system is stable.

Step 3: The value .tau..sub.cr of .tau. for which .delta.(.tau.), w(.tau.)=V.sub.cr provides an estimate of the true critical clearing time .tau..sub.c.

This procedure is known as the critical UEP technique.

The main computational difficulty in implementing this procedure is the calculation of the UEP in step 1. Observe that (.delta..sub.e,w.sub.e) is an equilibrium of the aforementioned equations: ##EQU5## if and only if w.sub.s =0 and ##EQU6##

These are essentially the real power flow equations and are difficult to solve.

Moreover as observed already, R(V.sub.cr) is an underestimate of the stability region. Put differently, the boundary of the stability region is poorly approximated by the surface of constant energy {V(.delta.,w)=V.sub.cr }. For the same reason .tau..sub.cr is an underestimate of the true critical clearing time .tau..sub.c.

A much improved estimate of .tau..sub.c can be obtained. It can be shown that the boundary of the stability regions is composed of the stable manifolds of the UEP that lie on the boundary. For instance in FIG. 2, the stability boundary is the union of the stability regions of UEP 4 and UEP 3. Also UEP 3 is the critical UEP and so V.sub.cr <V(UEP i) for i=3,4. Now consider the faulted trajectory (.delta.(.tau.), w(.tau.)) as shown. It leaves the stability region by exiting the boundary through the region of stability of UEP 2. Hence, an estimate of the critical time is given by .tau..sub.co, where:

Since V.sub.co >V.sub.cr, this estimate will be superior to .tau..sub.cr, obtained previously. The equilibrium UEP 2 is called the controlling UEP corresponding to the faulted trajectory. Note that the controlling UEP depends on the faulted trajectory, but that the critical UEP does not.

A two step controlling UEP procedure can be outlined:

Step 1: From the faulted trajectory .delta.(.tau.),w(.tau.)) estimate the controlling UEP (.delta..sub.co,w.sub.co). Calculate V.sub.co = V(.delta..sub.co, w.sub.co).

Step 2: The value .tau..sub.co again by

gives an estimate of the critical clearing time .tau..sub.c.

Once the critical clearing time is established, a system can be controlled. Appropriate measures can be used to bring generators back on-line more rapidly, or to shutdown certain generators of the system if post-fault stability will not be reached.

A method for determining power system transient stability is presented for controlling power systems on-line. This method is based on the direct (Lyapunov theory) approach controlling UEP technique.

The fundamental problem with the Lyapunov technique is the difficulty of calculating the controlling Unstable Equilibrium Point, which is a computationally intensive problem, even for a computer. The Lyapunov technique, therefore, cannot be used on-line. The invention comprises an on-line approach utilizing a modified Potential Energy Boundary Surface determination representing a stability region. This region is used to predict the critical clearing time for various on-line faults in a power system. The invention eliminates the difficulty of having to find the controlling UEP of the original system (3). The method of the invention determines the controlling of a reduced system and relates the controlling UEP of the reduced system to the controlling UEP of the original system.

The results of the inventive method were in close agreement with calculated conventional off-line results utilizing the cumbersome numerical integration of the differential equations describing the power system. The conventional technique is not practical, but has been used herein as a check upon the accuracy and precision of the invention.

It should be understood that, while commercial high power generating and distributing systems are described primarily herein, the invention is not limited to such industrial systems, but may also find application in any other electrical power systems such as telecommunications, computer networks, distributed robotic controls and the like.

A complete understanding of the present invention may be obtained by reference to the accompanying figures, when taken in conjunction with the detailed description thereof and in which:

FIG. 1 is a schematic diagram of the classical model of an electric power network;

FIG. 1a shows a fault state diagram;

FIG. 2 is a sample state space diagram illustrating the concepts of faulted trajectory, controlling UEP and critical UEP;

FIG. 2a depicts a stability region diagram;

FIGS. 3a and 3b show post-fault system state space diagram for a one machine, infinite bus model;

FIGS. 4a and 4b show a one machine, infinite bus schematic diagram and potential energy graph;

FIG. 5 shows a diagram of a two machine, infinite bus state space example;

FIG. 6 shows a diagram of three mathematical steps to determine the relationship between the stability boundaries of the original and the reduced system;

FIG. 7 illustrates a diagram of the relationship between the PEBS and the stability boundary of the original system;

FIG. 8 illustrates a diagram of the relationship between the stability boundary of the original system and the Potential Energy Boundary Surface;

FIGS. 9a and 9b illustrate diagrams of the fact that the connected constant energy surface differential V(.delta..sub.o) can approximate part of the stability boundary, but that connected constant energy surface differential V(.delta..sub.1,w.sub.1) may not; and

FIGS. 10a and 10b illustrate diagrams of when the connected constant energy surface differential V(.delta..sub.1,w.sub.1) is used to approximate the stability boundary.

Generally speaking, the invention features a new method of stability analysis for power systems.

Referring now to FIG. 1a, when a power system has had an event disturbance it can be said to experience a move from a Pre-Fault State 10 (a stable steady state) to a Fault-On State 14 coming to a Post-Fault State (a stable steady state) 18 during a period of time 22. Protective relays are placed strategically throughout the power system to detect faults and to trigger the opening of circuit breakers to isolate the fault. The fault occurs and the power system is then in the faulted condition 14 before it is cleared by the protective system operation. The transient stability problem is the study of the asymptotic stability of the post-fault system 18, with the initial condition given by the faulted system 14 response.

As depicted in FIG. 1a, the dynamic system 32 can be described by the set of differential equations: 12, 16 and 20 respectively for each state. x(t) is the vector of the state variables of the system at time t 22. At some time t.sub.f 26, the system undergoes a fault (structural change). This is represented by the change in the dynamics governing the system behavior from f.sub.I 24 to f.sub.F 26. Prior to the fault, the system is in a pre-fault 10 condition f.sub.I 24.

The fault duration is confined to the time interval {t.sub.F,t.sub.P }. During this interval, the system is governed by the fault-on 14 dynamics f.sub.F. At time t.sub.P 28, the fault is cleared and the system is henceforth governed by the post-fault 18 dynamics f 30.

During pre-fault 10, the system is at a known stable equilibrium point, x.sub.I reducing the analysis to the following equations:

At the end of the fault period the state of the system is

which represents the state of the fault-on 14 system 16 at time t 22 if it is in state x at time 0.

Post-fault 18 system 20 has a stable equilibrium point x.sub.s. Then f(x.sup.s)=0, and x.sup.s is stable. This identifies the fundamental problem of transient stability as follows: Starting in the post-fault 18 initial state x.sub.p, will the post-fault 18 system converge to the equilibrium point x.sub.s ?

Here .phi.(x,t) is the state of the post-fault 18 system 20 at time t if its state at time 0 is x. In other words, transient stability analysis essentially determines whether the fault-on trajectory at clearing time is lying inside the stability region (domain of attraction) of a desired stable equilibrium point (SEP) of its post-fault 18 system. Transient stability analysis does not estimate the whole stability boundary of the post-fault 18 system. Instead, only the relevant part of the stability boundary towards which the fault-on trajectory is heading is estimated.

Until recently transient stability analysis was performed only off-line as a function of power system reliability planning. The method used was a time-domain numerical simulation. This required computer programs which handled large numbers of buses and generators in detailed simulation models for a given fault and examined the simulated system response to these faults. Typically, ten or more simulations were required, with each simulation requiring 15 minutes of CPU time. As the requirements of reliability planning became more stringent, the move to on-line security analysis became more attractive as a critical guarantor of power system reliability. The time-domain simulations used off-line were inadequate, as the time interval of interest for transient stability analysis ranges from a few milliseconds to a few seconds.

One set of alternative approaches, known in the art as "direct methods" were proposed based on the Lyapunov theory. Critical to an understanding of these methods is stability boundary theory and the definition of the classical model.

Referring now to FIG. 2a, the following definitions are needed to discuss stability boundaries of nonlinear autonomous dynamical systems described by x=f(x) 34 where the vector field f is C' and it maps R.sup.N into R.sup.N :

(a) The point x.sub.i (where i=1,2, . . . ) 36, 38, 42 is equilibrium point (EP) if f(x.sub.i)=0.

(b) The derivation of f at x.sub.i is called the Jacobian matrix at x and is denoted by J(x.sub.i).

(c) An equilibrium point is hyperbolic if the Jacobian has no eigenvalues with zero real part.

(d) If all the eigenvalues of the Jacobian have a negative real part then it is a stable equilibrium point (SEP) 42.

(e) An EP is an unstable equilibrium point (UEP) if not all of its eigenvalues have a negative real part 38.

(f) If all the eigenvalues of an EP have a positive real part it is called a source.

(g) The type of the EP is defined to be the number of eigenvalues with a positive real part. Thus if an EP has exactly one eigenvalue with positive real part it is called a type-one equilibrium point.

(h) For an SEP x.sub.s, there is a region in the state space from which trajectories converge to x.sub.s. This is called the stability region 39.

The following assumptions are necessary when considering the stability boundary of nonlinear autonomous dynamic systems:

(a) That a stable equilibrium point (SEP) exists.

(b) All the equilibrium points (EP) on the stability boundary are hyperbolic.

(c) The intersection of the stable manifold (W.sup.s (x.sub.i)) and the unstable manifold (W.sup.u (x.sub.j)) satisfies the transversality condition for all the EP's x.sub.i,x.sub.j, on the stability boundary.

(d) There exists an energy function (C') such that:

i. V(.PHI.(x,t)).ltoreq.0 at x.epsilon.E;

ii. If x is not an equilibrium point, then the set (t .epsilon.: V(.PHI.(x,t))=0) has the measure 0 in R;

iii. V (.PHI.(x,t)) is bounded implies .PHI.(x,t) is bounded.

Conditions (i) and (ii) above imply that every trajectory must either go to infinity or converge to one of the equilibrium points 41.

The two fundamental theorems, well known in the art and based on the above assumptions, concerning the topological properties of stability regions, show that the stability boundary is composed of the union of all stable manifolds of the UEP's on the stability boundary.

Theorem 1 (necessary and sufficient condition for UEP to be on the stability boundary) states x.sub.i is an unstable equilibrium point (UEP) on the stability boundary .differential.A(x.sub.s) of a stable equilibrium point (SEP) x.sub.s if and only if W.sup.u (x.sub.i).andgate.A(x.sub.s).noteq..PHI..

Theorem 2 (characterization of the stability boundary) states: let x.sub.i, where i=1,2, . . . , be the unstable equilibrium points on the stability boundary .differential.A(x.sub.s) of a stable equilibrium point x.sub.s. Then .differential.A(x.sub.s)=x.sub.i .epsilon.E.andgate..differential.AW.sub.s (x.sub.i).

In addition, the energy function and the equilibrium points theorem proves that:

(a) The point at which the energy function achieves the minimum is the equilibrium point x.sub.i itself.

(b) On the stability boundary of x.sub.s, the point at which the energy function achieves the minimum and typically is a type-one equilibrium point.

(c) If it is assumed that the stability region A(x.sub.s) of x.sub.s is bounded, then on the stability boundary .differential.A(x.sub.s), the point at which the energy function achieves a local maximum must be a source.

Considerable progress has been made in the application of direct methods to power system transient stability analysis. As aforementioned, power system transient stability is concerned with a system's ability to obtain an acceptable steady-state operating condition following a large disturbance, called a fault. The mathematical model of a power system for transient stability analysis is a set of nonlinear differential equations. Mathematically, the transient stability analysis is to determine the stability of the post-fault system trajectory.

Consider the state space of the post-fault system, as shown in FIG. 3a. Let x.sub.s be the stable equilibrium point (SEP), that is an acceptable steady state operating condition. If the initial state of the post-fault system x, (i.e., the state at which the fault is cleared) lies inside the stability boundary, then the post-fault system trajectory will be stable. To determine whether the post fault system trajectory is stable, the so-called closest unstable equilibrium point (UEP) method compares the value of the energy function V(.) at every UEP (x.sub.1, x.sub.2, . . . ), finds the lowest value V1 and defines a region, the set {x: V(x)<V1}. This region is a subset of the stability region, as shown in FIG. 3a. Clearly, for a fault on a trajectory moving towards x.sub.2, the region determined above yields a result that is too conservative to be useful.

Another approach, called the controlling UEP method, suggests the use of the region employing the value of the energy function at x.sub.2 (i.e., the set {(x: V(x)<V2=V(x.sub.2)}) as a local approximation of the stability boundary. The problem with this method is that there is no simple way of determining the controlling UEP.

The inventive method of this invention is intended to circumvent the need to determine the controlling UEP of the original system. For a given fault-on trajectory, the method finds a "local" approximation of the stability boundary as depicted in FIG. 3b. The process of finding this local approximation is associated with the determination of the stability boundary of a lower dimensional system and is computationally rather efficient. The derivation is described below.

Consider the one machine infinite bus system in FIG. 4a. The system is described by the following equations: ##EQU7## where x=.delta., the phase angle of the machine with respect to the infinite bus. The equilibrium points are at: ##EQU8## Let x.sub.s be the SEP of interest. Consider the energy function for (1-1): ##EQU9## The energy function can be divided into "kinetic energy" and "potential energy" functions: ##EQU10##

At the equilibrium points of systems of equation (1-1), notice that the derivative of x with respect to time equals zero, and V.sub.p (x) is a local maximum of V.sub.p with respect to x. The potential energy V.sub.p (.) as a function of x is shown in FIG. 4b. The stability boundary of the system is shown in FIG. 3a. Note that:

a) The intersection of the stability region of x.sub.s and the set {(x,x): x=0, x element of set R} is P:={x,x): x.sub.2 <x<x.sub.1, x=0}

b) The boundary of this one dimensional region P in the x space consists of two points (0 dimensional ) x.sub.1 and x.sub.2.

c) These two points are characterized as being the local maxima of the potential energy function V.sub.p (.).

The PEBS method is based on the above facts and it proceeds as follows:

1) From the potential energy function V.sub.p (.), find the local maxima x.sub.1 and x.sub.2.

2) For those faults on trajectories (x(t),x(t)), whose x(t) component moves towards x.sub.1, the set {(x,x): V (x,x)=V.sub.p (x.sub.1), is chosen as the local approximation of the stability boundary. For those faults on trajectories (x(t),x(t)), whose x(t) component moves towards x.sub.2, the set {(x,x): V(x,x)=V.sub.p (x).sub.2 } is chosen as the local approximation of the stability boundary. It should be pointed out that for the one machine infinite bus system (i), the potential energy boundary is determined solely by V.sub.p (x), and (ii) in this case. The PEBS method is a method that determines just the controlling UEP.

The extension of the above ideas to multi-machine cases is non-trivial. Consider, for example a system consisting of two machines and an infinite bus. In this case, x= (.delta..sub.1,.delta..sub.2) and x=(.delta..sub.1, .delta..sub.2). The potential energy V.sub.p (.) as a function of x is shown in FIG. 5. Again the UEP's x.sub.1,x.sub.2,x.sub.3, . . . are the local maxima of the potential energy function V.sub.p (.) along a certain direction. For example, x.sub.1 is the local maximum point of V.sub.p (.) along the direction indication represented by the solid line through x.sub.1 in FIG. 5. A Potential Energy Boundary Surface (PEBS) is formed by joining these V.sub.p maxima through the "ridge" of the potential energy level surface, shown as a dotted line in FIG. 5. The PEBS constructed this way is clearly orthogonal to the equi-potential curves V.sub.p (.). Furthermore, along the direction orthogonal to the PEBS, potential energy V.sub.p (.) achieves a local maximum at the PEBS.

In the report by N. Kakimoto, Y. Ohsawa and M. Hayashi, "Transient Stability Analysis of Electric Power System via Lure-Type Lyapunov Function, Part I and II", Trans. IEE of Japan, Vol-98, pp. 516, 1978, the authors suggested a two-step procedure as the basis of the PEBS method:

Step 1) Calculate the fault on trajectory (x(t),x(t)). Suppose x(t) crosses PEBS at the value V.

Step 2) Use the surface {(x,x): V(x,x)<V} as a local approximation of the stability boundary. The derivation of the PEBS method for the multi-machine case, as indicated above, is based on heuristic arguments. Theoretical justification for the PEBS procedure, except for the one machine infinite bus case, is lacking. It is not clear for the multi-machine case whether the PEBS is the intersection of the stability boundary of the multi-machine system with the subspace {(x,x): x=0}. Moreover, because of the heuristic way PEBS is derived, it is difficult to see whether the method provides a good local approximation of the stability boundary, or under what conditions it does.

PEBS can be viewed as the stability boundary of the associated gradient system: ##EQU11## of the original system, which is the classical model of power systems consisting of N generators: ##EQU12## where x,y is an element of R.sup.n, and M,D are positive diagonal matrices. The motivation behind this approach is the fact that the potential energy V.sub.p (.) is the energy function of the gradient system (1-3). Moreover, note that all the equilibrium points of (1-4) lie on the subspace {(x,y): y=0} and at these points V(x,y)=V.sub.p (x).

To determine the relationship between the stability boundary of system (1-4) and the stability boundary (PEBS) of system (1-3), equation (1-4) is transmuted step by step into (1-3), as shown in FIG. 6. The motivation behind the transmutation in FIG. 6 is that it is desired to compare and study the relationship between the stability boundaries of each pair of systems. Equation (1-4) is rewritten as: ##EQU13## and equation (1-3) is rewritten as: ##EQU14##

If it were possible to eliminate the term (m,x) in (1-5), it would be possible to simply compare: ##EQU15## To eliminate the term (M), first consider an intermediate step, where the coefficient matrix is small, i.e., ##EQU16## In other words, determining the relationship between the stability boundaries of (1-3) and (1-4) is done through the following three steps.

Step 1: determine the relationship between the stability boundaries of the systems ##EQU17## where D.sup.-1, the inverse of matrix D, is also a positive diagonal matrix.

Step 2: determine the relationship between the stability boundaries of the systems ##EQU18## where .epsilon. is a small positive number and I denotes the identity matrix.

Step 3: determine the relationship between the stability boundaries of the systems ##EQU19## Thus, the relationship between the stability boundaries of (1-3) and (1-4) is established by connecting the relationship between the stability boundaries of (1-9a) and (1-9b), (1-10a) and (1-10b), (1-11) and (1-12).

Instead of studying (1-9a) and (1-9b) in step 1, generalize the problem and determine the relationship between the stability boundaries of the following generalized gradient systems:

where: M.sub.1, M.sub.2 are positive diagonal matrices, and f(x) is a gradient vector.

The approaches that are employed to establish the relationship between the stability boundaries in steps 1 and 2 are quite similar. First, derive a complete characterization of the stability boundaries of these systems. Second, conduct a qualitative analysis of the change in the stability boundaries when the vector fields of these systems are perturbed. In step 3, use the singular perturbation technique to show the relationship between the stability boundaries of (1-11) and (1-12).

Consider a power system consisting of N generators. Let the loads be modeled as constant impedances. Under the assumption that the transfer conductances are zero, the dynamics of the i.sup.th generator can be represented by the following equations: ##EQU20## where: node n+1 serves as the reference node, i.e., E.sub.n-1 =1 and .delta..sub.N+1 =0. The damping constant D.sub.i is assumed to be positive.

Power system transient stability analysis is concerned with the system's ability to reach an acceptable steady state operating condition following a large disturbance. Suppose a power system is operated at a stable equilibrium point x.sub.o, when a large disturbance (a fault) occurs. The system is then in the fault-on condition (it is a fault-on system) for brief duration before the fault is cleared by protective relay system operation, The system then enters the post-fault condition (e.g., the faulted line may be disconnected). If the system state at the fault clearing time is in the stability region of a stable equilibrium point of the post-fault condition, then the system is transiently stable.

The equations describing the pre-fault, fault-on, and post-fault systems all have the same form as (2-1) except Y.sub.ij 's are different due to the change in network topology.

Let (.delta..sub.s,w.sub.s) be a stable equilibrium point of (2-1). Define a function V(.delta.,w) as follows: ##EQU21## is the kinetic energy and V.sub.p (.) is the potential energy. Now, ##EQU22## and (2-1a), (2-1b) can be expressed as ##EQU23## Assume that zero is a regular value of the partial derivative with respect to .delta. of V.sub.p (.delta.). Note that M.sub.i >0, D.sub.i >0; the negative of the partial derivative, .differential.V.sub.p (.delta.) with respect to .delta. is a bounded vector field with bounded Jacobian; and the negative partial derivative, .differential.V.sub.p (.delta.) with respect to .delta. satisfies: ##EQU24## since it is a periodic vector in .delta.. The dynamic system (2-4) is, therefore, of the same form as the dynamic system

with (x,y)=(.delta.,w). Using the notation d.sub.p (M,D) to denote system (2-4), and the notation .differential.A.sub.p (M,D) to denote the stability boundary of a stable equilibrium point of D.sub.p (M, D).

Differentiating the function V(.delta., w) along the trajectories of system (2-1) or system (2-4) gives ##EQU25## The function V(.delta.,w) is therefore an energy function for the system .differential.p(M,D).

Since all of the equilibrium points of system (2-4) lie on the subspace

the energy function V(.delta.,w) at the equilibrium point (.delta..sub..epsilon.,w.sub..epsilon.) is of the following form:

This observation is motivation to study the stability region in the .delta. subspace {(.delta.,w): w=0} instead of in the whole state space. In other words, instead of studying the original system (2-4) consider the following gradient system: ##EQU26##

Define a dynamic system: ##EQU27## where: D is a positive diagonal matrix. The system d.sub.p (D) becomes the gradient system (2-6) when the matrix D is an identity matrix. Note that d.sub.p (D) is of the same form as the generalized gradient system d(D) defined previously as:

Let dA(.delta..sub.s) denote the stability boundary of a stable equilibrium .delta..sub.s of system (2-6) (or d.sub.p (I)). Next, define the PEBS: The stability boundary differential dA(.delta..sub.s) of (2-6) is the PEBS.

Below, it is shown that the PEBS defined above is consistent with the geometrical construction suggested by Kakimoto:

The PEBS intersects the level surface {.delta.: V.sub.p (.delta.)=c} orthogonally.

Proof:

It can be shown that: ##EQU28## where .delta. are the equilibrium points on the stability boundary differential A(.delta.). Since the gradient of V.sub.p (.delta.) is a vector perpendicular to the level surface {.delta.: V.sub.p (.delta.)=c} in the increasing direction, it follows that the vector field is perpendicular to the level surface {.delta.: V.sub.p (.delta.)=c}. Thus Theorem 2-1 is true.

Theorem 2-1 implies that in the direction orthogonal to the PEBS, the potential energy function V.sub.p (.) achieves a local maximum at the PEBS. Kakimoto, cited hereinabove, extends this fact and suggests an approximating scheme to find the PEBS. The goal of this scheme is to find the point at which the directional derivative of V.sub.p along the direction of the fault on trajectory achieves a local maximum. It is possible to go even further in suggesting a working definition of PEBS: consider that all rays emanating from x.sub.s, and the local maxima of V.sub.p (.) on these rays are joined together as the "PEBS".

The intuition behind the use of the potential energy maxima is that they represent the maximum energy (along the fault on trajectory) that can be converted into potential energy. Hence, if the fault is cleared before this point, all the kinetic energy can be converted into potential energy before the fault-on trajectory exits the stability region.

Essentially, the PEBS method proposed by Kakimoto for transient stability analysis provides a local approximation of the stability boundary of the original system (2-4). The application of the PEBS as a local approximation of the stability boundary proceeds in two steps.

Step 1: From the direction of the fault-on trajectory (.delta.(t),w(t)), detect the point .delta..sub.* at which the projected trajectory (.delta.(t)) crosses the PEBS. Let the value of V.sub.p (.) at .delta..sub.* be v. That is, v=V.sub.p (.delta..sub.*).

Step 2: Use the connected constant energy surface of the set {(.delta.,w): V(.delta.,w)=V.sub.p (.delta..sub.*)} containing the stable equilibrium point (.delta..sub.s,0) as the local approximation of the stability boundary differential, dA(.delta..sub.s, 0).

Next, evaluate the local approximation of the stability boundary that the PEBS method provides. For this purpose, compare the stability boundary of the original system d.sub.p (M,D) and the stability boundary of the gradient system d.sub.p (I) (i.e., PEBS). Recall that the stability boundary of (2-4) or (2-6) is the union of the stable manifolds of the unstable equilibrium points on the stability boundary. For this reason, first examine the relationship between the equilibrium points on the stability boundaries of d.sub.p (D) and d.sub.p (M,D).

The following theorem reveals the relationship between the equilibrium points of the gradient system d.sub.p (D) and the original system d.sub.p (D,M).

(i) (.delta.) is an equilibrium point of d.sub.p (D) if and only if (.delta.,0) is an equilibrium point of d.sub.p (M,D).

(ii) (.delta.) is a stable equilibrium point of d.sub.p (D) if and only if (.delta..sub.s,0) is a stable equilibrium point of d.sub.p (M,D).

(iii) N.sub.s (J(d.sub.p (D))at .delta.)=N.sub.x (J(d.sub.p (M, D)) at .delta.,0) and N.sub.c (J(d.sub.p (D))at .delta.)=N.sub.s (J(d.sub.p (M,D)) at .delta.,0).

The relationship between the equilibrium points on the stability boundary of the original system (2-6) and the equilibrium points on the PEBS is shown below in theorem 2-3.

If the dynamical systems d.sub.p (M,D) and d.sub.p (D) satisfy assumptions C.sub.1 and C.sub.2 then

(i) the equilibrium points (.delta..sub.i) on the PEBS correspond to the equilibrium points (.delta..sub.i,0) on the differential A.sub.p (M,D) and ##EQU29## where the following are assumed: Assumption C.sub.1 : The intersections of the stable and unstable manifolds of the equilibrium points on the stability boundary satisfy the transversality condition.

Assumption C.sub.2 : The number of equilibrium points on the stability boundary is finite.

This theorem asserts that, if the one parameter transversality condition is satisfied, and the number of equilibrium points on the stability boundary is finite, then the equilibrium points on the PEBS correspond to the equilibrium points on the stability boundary of the original system (2-4). Moreover, the PEBS is the union of the stable manifolds of the equilibrium points, say (.delta..sub.i), i=1,2 . . . n and the stability boundary of the original system is the union of the stable manifolds of the equilibrium points (.delta..sub.i 0),i=1,2 . . . n.

It is still uncertain whether the stable manifold in the PEBS and the intersection between the subspace {(.delta.,w): w=0} and the stable manifold of the corresponding equilibrium point in the stability boundary of the original system (2-6) are identical (see FIG. 7). Therefore, unlike the one machine infinite bus case, it cannot be claimed that the PEBS is the intersection of the stability boundary with the subspace {(.delta.,w): w=0}.

Now, consider under what condition the PEBS method gives an approximation of the relevant local stability boundary of the original system (2-4). Recall that the point at which the fault-on trajectory leaves the stability boundary of the post-fault system is called the exit point. Note that the stability boundary of system (2-4) or system (2-6) is the union of the stable manifolds of the equilibrium points on the stability boundary. What is meant by the relevant local stability boundary is the stable manifold on which the exit point of the fault-on trajectory lies.

Given a fault-on trajectory (.delta.(t),w(t)), let .delta..sub.t, w.sub.t denote the exit point on the stability boundary of the original system (2-4), and let .delta..sub.p denote the exit point of the projected trajectory .delta.(t), which is the projection of the fault-on trajectory in the .delta. space of the reduced system (2-6) (see FIG. 8). Let dV(x) denote the connected constant energy surface of V(.) containing the stable equilibrium point (.delta..sub.s,0) passing through the point x (i.e., dV(x):=the connected component of the set {(.delta.,w): V(.delta.,w)=V(x)} containing the stable equilibrium point (.delta.,0)). The idea of the PEBS method is to use the set differential V(.delta..sub.i,0) to approximate the relevant local stability boundary of the original system (2-4).

To see how good the set dV(.delta..sub.i,0) is in approximating the relevant local stability boundary of the original system (2-4), proceed with the following theorem. Use the notation V.sub.c (x) to denote the connected component of the set {(.delta.,w): V(.delta.,w)<V(x)} containing the stable equilibrium point (.delta..sub.s,0).

Let (.delta.,0) be an equilibrium point on the stability boundary differential, dA(.delta..sub.s,0) of system (2-4). Then

(i) the connected constant energy surface differential, dV(.delta.,0) intersects with the stable manifold W.sub.s (.delta.,0) only at the point (.delta.,0); moreover, the set V.sub.c (.delta.,0) does not contain any point which belongs to the stable manifold W.sub.s (.delta.,0).

(ii) Suppose (.delta..sub.1,w.sub.1) is an element of set {W.sub.s (.delta.,0), (.delta..sub.1,W.sub.1) not equal (.delta.,0)} and the equilibrium point (.delta.,j.sub.0) is of type k. Then the dimension of the set which is the intersection between the connected constant energy surface differential, dV(.delta..sub.1,w.sub.1) and the stable manifold W.sub.s (.delta.,0) is 2n-k-1; moreover, the set V.sub.c (.delta..sub.1,w.sub.1) does contain points which belong to the stable manifold W.sub.s (.delta.,0).

Part (i) of theorem 2-4 implies that for any fault on trajectory x.sub.f(t) starting from a point p with p part of the set A(.delta..sub.s,0) and V(p)<V(.delta.,0), if the exit point of this trajectory lies on the set W.sub.s (.delta.,0), then the fault on trajectory must pass through the connected constant energy differential, dV(.delta.,0) before it passes through the stable manifold W.sub.s (.delta.,0), (thus leaving the stability boundary d(.delta..sub.s,0). This suggests that the connected constant energy surface differential (.delta.,0) can be used to approximate part of the stability boundary dA(.delta..sub.s,0) (i.e., the stable manifold W.sub.s (.delta.,0) part). This is the essence of the controlling UEP method. (See FIG. 10a.)

Part (iii) of theorem 2-4 indicates that the fault-on trajectory x.sub.f(t) may pass through the connected constant energy surface dV(.delta..sub.1,w.sub.1) after it passes through the stable manifold W.sub.s (.delta.,0). It is therefore unwise to approximate the relevant part of stability boundary (i.e., the stable manifold W.sub.s (.delta.,0)) by the connected constant energy surface dV(.delta..sub.1,w.sub.1). (See FIG. 10b.) It is clear that part of the connected constant energy surface d(.delta..sub.1,w.sub.1) may lie outside the stability boundary dA(.delta..sub.s,0).

However, for the purpose of power system transient stability analysis, if the fault-on trajectory x.sub.f(t) passes through the surface differential dV(.delta..sub.1,w.sub.1) before it passes through the stable manifold W.sub.s (.delta.,0), the PEBS method still provides a good approximation of the relevant local stability boundary for this fault-on trajectory x.sub.f(t). (See FIG. 10a.) On the other hand, the PEBS method fails when the fault-on trajectory passes through the surface differential dV(.delta..sub.1,w.sub.1) after it passes through the stable manifold W.sub.s (.delta.,0). Note that in this case, the fault-on trajectory x.sub.f(t), after passing through the stable manifold W.sub.s (.delta.,0) becomes unstable for the post-fault system, but is still classified as stable by the PEBS method as long as this fault-on trajectory does not pass through the surface, differential, dV(.delta..sub.1,w.sub.1). (See FIG. 10b.)

Therefore, based on the above arguments, it follows that if the conditions below are satisfied:

(i) (.delta..sub.1) is on the stable manifold of (.delta.) implies (.delta..sub.i, w.sub.i) is on the stable manifold of (.delta.,0); and

(ii) the fault-on trajectory (.delta.(t),w(t)) passes through the constant energy surface differential, dV(.delta..sub.i,0) before it passes through the stable manifold W.sub.s (.delta., 0);

then the PEBS method gives a good approximation of the relevant local stability boundary of the original system for the given fault-on trajectory (.delta.(t),w(t)).

It is difficult to check whether condition (ii) above is satisfied. This is mainly because finding the exit point (.delta..sub.t,w.sub.t) for the fault-on trajectory (.delta.(t),w(t)) is difficult. However, a computational scheme to detect the exit point (.delta..sub.p) of the projected fault-on trajectory was proposed by Kakimoto, cited hereinabove. Theorem 2-4 provides a solution for this difficulty and suggests that the condition (ii) is always satisfied if the PEBS method is modified as follows:

Step 1': From the fault-on trajectory (.delta.(t),w(t)) detect the point .delta..sub.* at which the projected trajectory .delta.(t) crosses the PEBS. Let the value of V.sub.p (.) at .delta..sub.* be v (i.e., v=v.sub.p (.delta.)).

Step 2': Find the equilibrium point of system (2-6) whose stable manifold contains the point .delta..sub.*, say .delta.. Let the value of V.sub.p (.) at .delta. be v (i.e., v=V.sub.p (.delta.)).

Step 3': Use the connected constant energy surface of the set {(.delta.,w): V(.delta.,w)=V.sub.p (.delta.)} containing the stable equilibrium point (.delta..sub.s,0) as the local approximation of the stability boundary .differential.A(.delta..sub.s, 0).

Based on previous work, step 2' can be implemented by the following procedure:

Step (i): Use the point .delta..sub.* as the initial condition and integrate the post-fault gradient system (2-6) until the local minimum of V.sub.p (.delta.) with respect to .delta. is achieved, say .delta..sub.o.

Step (ii): Find the equilibrium point of the post-fault reduced system using .delta..sub.o as the initial guess (i.e., solve the equation derivative V.sub.p (.delta.) with respect to .delta.=0). The solution is the required one, .delta..

Table 1 below shows the results of a test on a 118-bus power system that compares the critical clearing time prediction of different approaches to transient stability analysis.

The controlling UEP method of the invention has been applied to a 118-bus system with 20 generators. Faults were created, including faults on both generator and load buses. Table 1 shows the critical clearing time of six faulted systems using three different methods: numerical integration technique, potential energy boundary surface method and the controlling UEP method of the present invention. This table demonstrates that the inventive method offers fairly accurate direct analysis of transient stability for these faults. Note that the PEBS method gives an over estimate of the critical clearing time when a fault occurs at bus #100, a severe fault and gives a serious over estimate when a fault; occurs at bus #68, a mild fault.

Through the application of the inventive method system, transient stability is assured. The inventive method provides the capability to predetermine a time to close a circuit breakers opened when the event disturbance occurred, thus clearing the fault-on situation and returning the system to a post-fault, steady state operating environment.

Since other modifications and changes varied to fit particular operating requirements and environments will be apparent to those skilled in the art, the invention is not considered limited to the example chosen for purposes of disclosure, and covers all changes and modifications which do not constitute departures from the true spirit and scope of this invention.

For example, while commercial high power generating and distributing systems are described primarily herein, the invention is not limited to such industrial systems, but may also find application in any other electrical power systems such as telecommunications, computer networks, distributed robotic controls and the like.

Having thus described the invention what is desired to be protected by Letters Patent is presented by the appended claims.