{'fraction_non_alphanumeric': 0.03421595598349381, 'fraction_numerical': 0.013554562127464465, 'mean_word_length': 4.676154847104749, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Modal-Space Method for Online Power System Steady-State Stability Monitoring.

This paper proposes a novel approach to estimate the steady-state angle stability limit (SSASL) by using the nonlinear power system dynamic model in the modal space. Through two linear changes of coordinates and a simplification introduced by the steady-state condition, the nonlinear power system dynamic model is transformed into a number of single-machine-like power systems whose power-angle curves can be derived and used for estimating the SSASL. The proposed approach estimates the SSASL of angles at all machines and all buses without the need for manually specifying the scenario, i.e. setting sink and source areas, and also without the need for solving multiple nonlinear power flows. Case studies on 9-bus and 39-bus power systems demonstrate that the proposed approach is always able to capture the aperiodic instability in an online environment, showing promising performance in the online monitoring of the steady-state angle stability over the traditional power flow-based analysis.


I. INTRODUCTION
The principal cause of the Northeast Blackout in 2003 was a lack of situational awareness due to inadequate reliability tools [1]. The comprehensive and accurate assessment of reliability, i.e. dynamic security assessment (DSA), requires solving computational-extensive numerical integrations which make DSA viable only for offline planning or long-cycle online applications [2]. In practical power system planning and operations, there are usually several thousand contingencies and tens or hundreds of interface limits to be calculated. In addition, considering more diverse power flow patterns brought by a high penetration of intermittent energy resources, future DSA is expected to be executed more frequently with a shorter cycle or even approaching real-time. These trends motivate the research needs for fast stability assessment tools to increase the situational awareness and lower the risk of large-scale blackouts.
A promising way to address such a challenge is the steadystate stability analysis, which only focuses on the steady-state condition of power systems such that the stability assessment can be significantly simplified. The static stability of the base case (without any contingencies) is a necessary condition of system stability and it has been proved practically useful [2] [3], including estimating the stability margin as well as other purposes such as the forward reserve procurement. The basic idea is to identify the steady-state voltage stability limit (SSVSL) in a very speedy fashion by solving static power flow Corresponding author: Le Xie equations such that SSVSL can be calculated and monitored online. The online implementation of this type of approaches will guarantee the system operator's awareness of the situation where a margin becomes dangerously low, allowing preparation and execution of necessary control actions. Such an application seems to be more needed in future grids with a high penetration of intermittent renewable energy where system conditions may change more frequently and significantly.
However, there are three major challenges with existing methods for assessing SSVSL when dealing with future power grid with massive renewables: (i) selection of interfaces and scenarios, i.e. which interface is of concern and how to select the sink and the source [4] [5]; (ii) need for solving a series of power flows for each interface and for each scenario; and (iii) SSVSL is an optimistic stability analysis since smallsignal instabilities may still occur without violating SSVSL (to be shown in case study section of this paper). Even regardless of the computational burden involved in (ii), the selections in (i) are always conducted manually, where peoples experiences and understanding about historic behaviors of the power grid play an important role. Such manual selections can be either difficult or unreliable, especially when considering the power flow pattern and critical interfaces may change more significantly and more frequently with loading and operating conditions for future power grids with a high penetration of intermittent distributed energy resources. To this end, an automated, fast and less-optimistic method is proposed in this paper to estimate the steady-state angle stability limit (SSASL), which captures the Saddle-node bifurcation points along certain stressing directions defined by the associated linearized dynamic system. Although these stressing directions may be unrealistic, it is found by extensive numerical studies that the stability margin defined by the proposed approach along any realistic stressing direction would always become low when approaching the stability boundary.
The key contributions of this paper are summarized below. 1) We explore the nonlinear modal structure of nonlinear power system dynamic model and propose a direct analysis of SSASL leveraging the nonlinear modal structure. We name it the modal space method, or MS method.
2) The proposed MS method is free of the selection of interface and scenario, and is always slightly conservative compared with the small-signal analysis (SSA). 3) We propose an MS method-based online stability monitoring application, which is shown to be less-optimistic than the traditional power flow based SSVSL analysis. 4) We perform extensive numerical studies on the 9-bus and 39-bus power systems, compare with the traditional SSVSL analysis and the standard SSA, and demonstrate the potential of the proposed MS method.

II. A BRIEF REVIEW OF SELECTED STABILITY PROBLEMS
This section briefly reviews a few selected power system stability problems, i.e. TSA, SSA and static voltage stability analysis (VSA), to point out the scope of this work. The review is neither exhaustive nor comprehensive, and we refer interested readers to reference [6] for a comprehensive review of power system stability problems.
A. TSA, SSA and static VSA TSA of nonlinear dynamical power systems is usually considered to be the most realistic dynamic performance analysis, which is used for benchmarking other simplified stability analyses, e.g. small-signal and voltage stability analyses, and evaluating new control schemes, e.g. damping controls. TSA is often defined upon a set of nonlinear differential-algebraic equations (DAEs), whose accurate assessment by a numerical integration can be very computationally-expensive [7]. Thus, in practical applications, numerical integration based TSA is used for most utilities, while for some utilities it is used as the final check of very few scenarios screened out by fast but not very accurate analyses, e.g. a direct method [8]. With today's analytical and computing capability, it is still extremely difficult, if not impossible, to implement online TSA for analyzing all contingencies and all potential changes in loading and generation conditions. SSA is a common simplification of the stability analysis of nonlinear dynamical systems, which is only valid in a small neighborhood of the equilibrium point, i.e. the steady-state condition, since all nonlinearities with power system DAEs are ignored. To evaluate the small-signal stability of power systems under other steady-state conditions or when subject to possible changes of loading and generations, SSA needs to be re-executed, which is also considered to be computationalexpensive for the computing resources in today's control room.
Static voltage stability analysis (VSA) is a further simplification which ignores all dynamics, while only retaining Kirchhoff's current and voltage laws in the power network, resulting in a set of nonlinear algebraic equations known as power flow equations. Static VSA may imply SSA under some extreme and unrealistic assumptions, while in general there is no such an implication [9][10] and voltage stability is usually an optimistic estimate of small-signal stability as shown in Fig. 1. Most production-level solutions of power flow equations are iterative in nature, e.g. Newton-Raphson and fast decouple methods, although computationally demanding, which are computationally much cheaper than SSA and TSA.

B. Scope of this work
There are two mechanisms for small-signal instability over slow changes of system conditions or parameters [10], i.e. self-oscillation instability and aperiodic instability, which respectively correspond to Hopf-bifurcation and Saddle-node bifurcation in the field of nonlinear dynamics. This work only 17 Voltage stability boundary focuses on the aperiodic instability, i.e. Saddle-node bifurcation induced instability. The red curve in Fig. 2 represents a typical self-oscillation instability, where the real parts of two conjugate eigenvalues change from negative to positive at points C so heading for an instability without meeting each other. The self-oscillation instability is mostly caused by improper designs of automatic control [10], and it is highly related to the damping of oscillatory modes. Therefore, the damping ratio of poorly-damped modes deserve dedicated monitoring programs in operations and planning [11] to avoid self-oscillation instability. The aperiodic instability refers to the mechanism that two conjugate eigenvalues first collide on the real axis at point P co , split into two real eigenvalues and one of them becomes positive at point C ap leading to an aperiodic instability, as illustrated by black curves in Fig. 2. Note that the aperiodic instability represents a loss of synchronism among generators, i.e. steady-state angle instability [10]. This work aims at a novel online steady-state angle stability monitoring application, where the nonlinearity of system dynamics are retained and utilized to directly infer the steadystate angle stability limit (SSASL). The proposed stability analysis is expected to be more accurate, i.e. less optimistic, than static VSA, while being slightly optimistic only when self-oscillation instability occurs before an aperiodic instability, as shown in Fig. 1.

III. POWER SYSTEM NONLINEAR DYNAMIC MODEL IN MODAL SPACE
Consider a classic N -machine power system where δ sij is the steady-state angle difference between generators i and j, ω s , P mi , where all loads are represented by constant impedance and included in parameters C ij and D ij . For simplicity, (1) can be re-written as (3) The models in (1) or (3) are called the nonlinear dynamic model in angle-speed space in this paper, while the one to be derived in the following is called the nonlinear dynamic model in modal space. Eigen-analysis based on the linearized model of (3) is briefly reviewed below, whose eigenvector matrix will be adopted in deriving the nonlinear dynamic model in modal space. It should be emphasized that although a linear change of coordinate is adopted in our derivation, the nonlinearity of system dynamical model is fully retained without any approximation.
Step 1: Linearize (3) at the origin and obtain (4), where A is the Jacobian of f at the origin.
Step 2: Calculate the eigenvalues and eigenvectors of A and define the transformation in (5), where P is the matrix consisting of the right eigenvectors of matrix A, and y is the state vector in modal space, i.e. y = {y 1 , y 2 , · · · , y 2N }.
In the eigen-analysis, the linear approximation in (4) is utilized to study the dynamic behaviors of (3) when subject to small disturbances. As a comparison, the idea to be proposed below will fully maintain all nonlinearities with the original DEs in (3). With the linear change of coordinates in (5), the nonlinear DEs in (3) are transformed to (7), called the dynamic model in modal space. Note that the system in (7) is mathematically equivalent to (3) given the matrix P in (5) is invertible. All subsequent analyses are performed on the nonlinear DEs in (7). Remarks: Corresponding to the concept of the mode in the linear analysis (4)-(6), the two nonlinear DEs corresponding to a conjugate pair of eigenvalues and their dominated nonlinear dynamic behaviors are referred to as a nonlinear mode in this paper. Unlike linear analysis where any two modes and their dynamics are independent, two nonlinear modes are only linearly independent but still coupled nonlinearly, which is termed the nonlinear modal interaction in normal form analysis [13], i.e. dynamics initiated in one nonlinear mode may propagate to and affect the dynamics of another nonlinear mode through high-order nonlinear terms.
It has been proved in [14] that for any classic N -machine power system with a uniform damping, the nonlinear dynamics associated with the two real modes do not affect the nonlinear dynamics associated with complex modes. Therefore, the angle stability of systems in (7) are dominated by its first (2N − 2) nonlinear DEs of (7). Expanding (7), we can obtain (8) whose first (2N − 2) nonlinear DEs in y 1 , y 2 , · · · , y 2N −2 are completely decoupled from y 2N −1 or y 2N , i.e. the states corresponding to the two real modes.

IV. A MODAL SPACE METHOD TO ESTIMATE SSASL
In this section, we propose the modal space method to estimate the SSASL of a general classic N -machine power system using its nonlinear dynamic model in modal space, as shown in (8). Three versions of the proposed method are developed with different requirements on the system steady-state condition at limit points. The following starts with the introduction of the proposed MS method, continues with a discussion on their differences, and ends up with an introduction of an online stability monitoring application.

A. Simplifying nonlinear dynamic model in modal space
This subsection presents a simplification of (8) into a number of single-degree-of-freedom (SDOF) nonlinear dynamical systems by using the steady-state condition. Note that the first (2N − 2) nonlinear DEs of (8) are still coupled together through nonlinear modal interactions, which place a difficulty for further theoretical analysis of system dynamics when subject to large disturbances. Fortunately, since this work focuses on the steady-state behavior of the system, i.e. how the equilibrium changes nonlinearly with loading condition/generation dispatch, instead of dynamic behaviors [15] [16], then it is reasonable and helpful to make the following simplification such that (8) becomes (9): the steady-state condition of the system in (8) implies that when analyzing any nonlinear oscillatory mode i ∈ {1, 2, ..., N − 1}, the impact from dynamics of all other modes can be ignored, i.e. y 2j−1 = y 2j = 0 for all j = i.
B. Constructing real-valued DEs for each nonlinear mode This subsection applies a linear change of coordinates to (9) and constructs a real-valued nonlinear SDOF dynamical system for each nonlinear oscillatory mode. The linear change of coordinates is intentionally designed to make the resulting real-valued system similar to a single-machine system, but with a differential characteristic of nonlinearities.
Note that the two state variables and the two nonlinear DEs w.r.t. the nonlinear mode i in (9) are complex-valued. More specifically, y 2j−1 and y 2j and the two DEs in (9) are respectively conjugate to each other, which is true if the two eigenvectors corresponding to a conjugate pair of eigenvalues of matrix P in (5) are conjugate to each other. The following adopts the linear change of coordinates [15] as shown in (10) to transform (9) to a new set of coordinates, such that the two new state variables turn out to be real-valued and possess similar meanings to angle and speed.
The determination of the transformation in (10) is intuitively generalized from the investigation of a general singlemachine-infinite-bus (SMIB) power system [17][18] as shown in (11), where the system has an equilibrium at origin and β = P max ω s /2H for simplicity (note that P m = P max sin δ s and P e = P max sin(δ + δ s )). The linearization of (11) at the origin is obtained as (12). By eigen-analysis, eigenvalues of (11)'s Jacobian are λ 1 = j √ β cos δ s and λ 1 = −j √ β cos δ s , and their right eigenvectors are shown in (13). The linear change of coordinates in (14) can transform (12) to its modal space as shown in (15), which corresponds to the linearization of (9). Thus, the transformation in (10) is chosen according to the inverse of matrix P , as shown in (16), which can transform (15) back to (12).
Adopt notations in (17) by separating the real and imaginary parts of (9) and its eigenvalues, where a i , b i , c i and d i are real-valued numbers or functions, then substitute (10) into (9) and obtain (18), which is real-valued. Note that the realness of (18) only relies on the conjugate relationships between the two DEs in (9) and between λ 2i−1 and λ 2i , as shown in (17).
C. Identifying generalized power-angle curve for each nonlinear mode It has been found by extensive studies on different system models and under different operating conditions that the equations in (18) always follow a special form as shown in (19), i.e. the derivative of the generalized speed is a singlevariate function of the generalized angle and the derivative of the generalized angle is exactly the generalized speed. This paper assumes such a relationship without proof. We further define the "generalized power" by (20) for the nonlinear mode i. Note that the desired information about SSASL is completely maintained in (20), which will be explored in the next subsection.

D. Estimating SSASL for each nonlinear mode
In an SMIB system, the SSASL is determined by the two closest extreme points on the power-angle curve around the stable equilibrium point as shown in Fig. 3, which are corresponding to the largest power transfer from (to) the machine to (from) the infinite bus. Similarly, the SSASL of each nonlinear mode can be estimated by identifying the two closest extreme points on the generalized power-angle curve around the origin. Since the generalized power-angle curve of each nonlinear mode is a univariate function in generalized angle as shown in (20), a bisection based numerical search for these two extreme points, say δ gi,1 and δ gi,2 , is not very computationally expensive and is therefore adopted in this paper. An approximate analytical solution may be possible and deserves further investigations.
It can be seen that two SSASL points can be obtained for each nonlinear mode. Thus, for the system in (8) having (N − 1) nonlinear modes, we have (2N − 2) SSASL points in total. An SSASL actually represents a system steady state where the system is about to lose stability. This subsection introduces procedures to determine such a system steady state, i.e. voltage magnitudes and angles of all buses, from the estimated SSASL for each nonlinear mode.
Note that the generalized angle and speed, i.e. coordinates of the system in (19), are transformed from original angle and speed derivations, i.e. coordinates of the system in (1), over two linear changes of coordinates in (5) and (10). Therefore, each estimated SSASL point represented by the generalized angle and speed can be transformed back to original anglespeed coordinates by the inverse of these two linear transformations, then all other states can be determined accordingly. Denote the two limit points w.r.t. nonlinear mode i as δ gi,k with k = 1 or 2. Detailed procedures are summarized below.
Step 3: Substitute y into (5) and obtain x.
Step 5: Calculate the terminal current vector I t at generator buses by I t = Y r E, where Y r is the reduced admittance matrix including generator source impedance and all load impedance.
Step 6: Calculate the terminal voltage vector V t at generator by (21), where Z s = diag{Z s1 , Z s2 , . . . , Z sN } is the source impedance matrix.
Step 8: By now, voltage phasors at all buses have been obtained, as shown in (21) and (23). Other states, including generator power output, load power consumption can be calculated respectively by (24)-(26), while line currents and line flows can be calculated accordingly by Ohm's law (omitted here).
We name the above derivation of the system steady state at an estimated limit point as the first version of the MS method, or MS1, to be distinguished from two more versions named MS2 and MS3 to be introduced below. Note that MS implicitly assumes that when the system goes from the current steady state to the steady state at the limit point, (i) generators' internal bus voltage magnitudes maintain unchanged and (ii) loads are modeled as constant impedance. These two implicit assumptions are not very realistic since (i) generators always have automatic voltage regulations to maintain the terminal voltage magnitude to be around a given reference, and (ii) load power consumption is more often used to measure the system stress level, instead of load impedance. Therefore, it would be more reasonable to maintain generator terminal voltage magnitude and load power unchanged. The following will briefly introduce MS2, which guarantees generator terminal voltage magnitude unchanged, and MS3, which guarantees both generator terminal voltage magnitude and load power unchanged.
MS2 starts with steps 1-6 of MS1. Then, the magnitude of each entry of the obtained vector V t is overwritten by the corresponding generator terminal bus voltage magnitude under the base-case condition, while all angles of V t are intact. The updated generator terminal voltage vector is denoted as V (update) t . The terminal current is updated by (27) and nongenerator bus voltage is updated similarly to (23) but with V (update) t . Finally, other states can be updated accordingly.
MS3 also starts with steps 1-6 of MS1, continues with all additional steps of MS2, then updates the load impedance using the solved bus voltage and the base-case load power consumption (such an operation alters the admittance matrix in (22), the non-generator bus voltage in (23) and generator terminal current), and finally iterates the updates of generator terminal voltage and current and load impedance until a convergence is reached.

F. Remarks on MS method
The proposed MS method approach exploits the nonlinearities within the power system nonlinear dynamic model and directly determines SSASL without repetitively solving the nonlinear power flow problem and the linear eigen-analysis. All system states, i.e. a solved power flow condition, at each SSASL point can be determined such that the stability limit can be reflected in all system states, including all transmission line flows and all bus voltages. Therefore, this new method eliminates the need for manually selecting interfaces, since all lines have been considered. In addition, this new approach also eliminates the need for manually selecting scenarios for stressing the system, since each estimated SSASL point represents the largest stressing condition of a nonlinear mode and such stressing directions are uniquely and automatically determined in the proposed method as implied by the generalized power-angle curve (20). A few more important remarks are summarized below. 1) Number: Two SSASL points can be obtained for each nonlinear oscillatory mode. Therefore, for a classical N -machine power system having (N − 1) oscillatory modes, there are (2N − 2) SSASL points. 2) Physical meaning: Each SSASL point is an estimate of Saddle-node bifurcation point, which represents the largest steady-state angle separations without causing steady-state instability when stressing the system about that nonlinear mode, i.e. increasing the power transfer by generation re-dispatch between two groups of generators determined by the mode shape of that mode. 3) Limitation in handling load change: As discussed in 2), the proposed approach assumes a constant load model (constant impedance in MS1 and MS2, while constant power in MS3) and then estimates the generation dispatchability limit. However, when the stability limit against load changes is of concern, the proposed approach seems currently not applicable in that regard. One would have to resort to other approaches, e.g. Thevenin's equivalent approach for static voltage stability. 4) Conservativeness: It is worth mentioning that the estimated SSASL corresponds to the point P co in Fig. 2, instead of C ap . Therefore, along any stressing direction, P co is always reached before C ap , which makes our proposed method conservative. Numerical studies will show that the MW distance between P co and C ap is almost negligible. Thus, the conservativeness of the proposed method is usually small. 5) Comparison of MS1, MS2 and MS3: From MS1 to MS3, the computational complexity increases while the accuracy is expected to improve. A brief summary is shown in Table. I.

G. Applying MS method to online stability monitoring
This subsection presents a way to apply the proposed MS method to the online steady-state angle stability monitoring. For any given operating condition of a classic N -machine power system, the proposed MS method can give (2N − 2) steady-state conditions respectively representing (2N −2) limit points on the aperiodic stability boundary. Therefore, distances between the current operating condition to those steady-state conditions at limit points can be used as stability margins for the monitoring purpose. Although any states can be used to define the distance, the MW power is usually a practically meaningful measure and is adopted to define the stability margin, as shown in (28).
where P = {P e1 , P e2 , . . . , P eN } is the generator active power vector at current operating condition, P SSASL,k represents the generator active power vector at the condition w.r.t. the kth SSASL point, k ∈ {1, 2, . . . , 2N − 2}, || * || is the norm operation (L 2 -norm is used in this paper). Thus, the proposed MS method defines (2N − 2) MW margins for a given operating condition. If monitoring these MW margins when the system is increasingly stressed, we will have (2N − 2) MW margin curves. An instability may occur when any of those curves approaches zero.

V. CASE STUDIES
This section illustrates the steady-state angle stability monitoring application based on the proposed MS method in detail on the IEEE 9-bus power system [19], including a comparison to the reference results by static VSA and SSA. The 39-bus power system [20] is then tested to show the effectiveness of the proposed method in a larger power system. All computations in this section are conducted in Matlab on a PC with dual-core Intel i5 at 2.6GHz and 8GB memory.

A. Tests on IEEE 9-bus power system
The IEEE 9-bus power system, whose one-line diagram can be found in [19] and is omitted here, is selected mainly because it is the smallest multi-machine power system whose steady-state angle stability boundary can be visualized in a 2-D plot, say P e2 -P e3 space or δ 21 -δ 31 space. There are three machines, therefore, defining two electromechanical modes. Based on the power flow data from [19], dynamic data from [21] and using classic generator model, the two electromechanical modes are found to be 1.38Hz and 2.13Hz, where the 1.38Hz mode represents the oscillation between generator 1 and generators 2 and 3, denoted as mode 1, while the 2.13Hz mode is between generator 3 and generators 1 and 2, denoted as mode 2.
The first test is to illustrate the accuracy of limit points identified by the proposed MS method by comparing to reference stability boundaries numerically identified by VSA and SSA. These reference results are obtained by a ray-scanning scheme, whose steps are summarized below for identifying the static voltage stability boundary, while steps for identifying aperiodic (or small-signal) stability boundaries only has a different step 2 where aperiodic (or small-signal) stability is also checked in addition to power flow convergence.
Step 1: Given an operating condition with a generation dispatch as P = {P e1 , P e2 , P e3 }, select a direction in the 2-D P e2 -P e3 space, say n = (∆P e2 , ∆P e3 ) and ||n|| = ∆, where ∆ represents the step size and ∆ = 10 MW is used in this paper to initialize the step size for each ray.
Step 2: Solve power flow with modified generations P e2 = P e2 + ∆P e2 and P e3 = P e3 + ∆P e3 . Step 3: If power flow in step 2 converges, then repeat step 2 to check a farther point along the ray in the 2-D space.
Otherwise, go to step 4.
Step 5: Recover generations to the last converged power flow by P e2 = P e2 − ∆P e2 and P e3 = P e3 − ∆P e3 , and then record (P e2 , P e3 ) as an estimate of the stability boundary.
Step 6: Go back to step 1 to start over with another direction until all desired directions are searched. Numerically identified reference stability boundaries and the SSASL points estimated by MS1, MS2 and MS3 are shown in Fig. 4. Several observations can be obtained: (i) VS boundary is optimistic in general compared to SSS boundary; (ii) AS boundary mostly coincides with SSS boundary while is slightly optimistic when self-oscillation instability occurs before an aperiodic instability; (iii) In δ 21 -δ 31 space, SSASL points estimated by MS1, MS2 and MS3 are fairly close to each other on the AS boundary, showing that all three versions are able to give accurate angle limit; and (iv) In P e2 -P e3 space, SSASL points by MS2 and MS3 are fairly close to the AS boundary while those by MS1 are not quite accurate, though being conservative.
The first test shows how accurate the estimated SSASL points are for a given operating condition by comparing to reference stability boundaries. In fact, if the actual stressing direction does not point to any SSASL points, the system will not exit the stability region through one of these SSASL points estimated at the base-case condition. This is not a problem, since we can always re-calculate these SSASL points by continuously applying the proposed MS method over the change of system steady state. To show an online stability monitoring application based on the proposed MS method, in the second test, the system operating condition is intentionally stressed in a specific direction. A number of steady states, to be checked by the proposed MS method, are selected between the base case and the stability limit. Fig. 5 shows the estimated SSASL points over the stressing process where P e3 increases from 85MW to 414MW (all loads and P e2 are maintained unchanged) to cause an instability. Five 1 represents the stopping criterion which takes 0.1MW in this paper. steady states, from the base case to the limit, are checked by the proposed MS method and the resulting SSASL points change from red to yellow. Fig. 5 shows that one of the four SSASL points, i.e. the one in the first quadrant, arrests the system when it tries to exit the stability region. It is also observed that when the system gets close to the stability boundary, other SSASL points, than the one arresting the system on the boundary, may not be very accurate. This is fine as long as there is always an SSASL point that accurately arrests the system when it tries to exit the stability region. This has been found true by exhausting all stressing directions in P e2 -P e3 space with a small resolution of 2 degrees. The above visualizes the accuracy of the proposed method and its application in online steady-state angle stability monitoring. However, such a visualization might not be possible for large power systems. To this end, the MW margin defined in (28) can help measure and visualize the distances from an operating condition to SSASL points on the AS boundary. Fig. 6 shows these MW margins when applying the proposed approach to multiple steady states over the stressing process, where the MW margins corresponding to the SSASL point for which the system is approaching are decreasing to zero, while MW margins of other SSASL points are relatively sufficient, either increasing or staying at 200MW or above.

B. Tests on New England 39-bus power system
This subsection applies the proposed MS method to a 39-bus power system [20], and two scenarios are tested. In the first scenario, P e37 increases and P e30 decreases by the same MW, resulting in a stress in the local power transfer and causing a small-signal instability when the change of P e37 (or P e30 ) reaches 2143.04MW. Along such a stressing direction, the system loses AS and VS respectively when the change of P e37 reaches 2142.95MW and 2422.42MW. The MW margins over stressing are monitored by the proposed method and shown in  Fig. 7, where MW margins of most SSASL points are above 1000MW even when being close to the instability while only two or three SSASL points may encounter low MW margins over stressing. Note that the smallest MW margin may switch from/to among these two/three SSASL points, as pinpointed in the black circles, over the stressing process, which calls for the need for monitoring all these critical SSASL points to not to miss any potential risk. It is also worth mentioning that (i) aperiodic instability and small-signal instability are extremely close to each other, i.e. less than 0.1MW in the MW change of the stress, and (ii) right before the system loses it aperiodic stability, the MW margin by (28) using the voltage stability limit is as huge as 482.68MW, while the MW margins by the proposed MS method are 54.5MW, 39.79MW and 13.70MW respectively for MS1, MS2 and MS3. Similar phenomena can also be seen in scenario 2, where P e33 , P e34 , P e35 , P e36 and P e38 increase and P e30 , P e31 , P e32 and P e37 decrease, stressing the interface between two distantly located generator groups. The current implementation of the proposed approach relies on the Symbolic Math Toolbox in Matlab, which turns out to be inefficient. In addition, estimating SSASL points involves the calculation of all eigenvalues and eigenvectors of the linearized system. For a single operating condition of the 39bus system, it takes up to 30s. A more efficient implementation without symbolic derivations, along with the selection of critical modes for a partial eigen-analysis are our future work.

VI. CONCLUSION
This paper proposes a modal-space (MS) method for estimating the steady-state angle stability limit (SSASL) using the power system nonlinear dynamic model in modal space. The MS method can estimate the SSASL for all system steady states in a single run. A steady-state angle stability online monitoring application is developed based on the MS method and tested on the IEEE 9-bus system and New England 39-bus system. Numerical results show that the proposed MS method is always able to show a meaningful MW margin and arrest the system when it tries to exit the aperiodic stability region.
To further show the potential in online environment, dedicated algorithms without symbolic derivations will be developed to reduce the computation. Future work will also consider N − 1 contingency analysis and remedial control actions.