A Framework for Analyzing System Loadability With Multiple VSCs Using a Hybrid Model

Due to the increasing number of power electronics-based devices in modern power systems, there are concerns about the impact of higher-frequency interactions on system stability. The typically-used algebraic representation of networks fails to capture these newly emergent issues, but dynamically modelling the entire network with differential equations will lead to a prohibitively long simulation time. This paper proposes a computationally-efficient analysis framework to determine the system loadability in power electronics-rich, large networks with the consideration of multiple types of stabilities. The framework developed in this paper identifies the critical network elements based on the eigenvalue sensitivity in order to exploit a hybrid modelling approach. Within the hybrid model, only the critical network portion is modelled dynamically with the rest of the network represented by algebraic equations. Bifurcation theory is used to simultaneously analyze both small-disturbance rotor angle and small-disturbance voltage stability. The results obtained show that the high-frequency interactions are accurately captured with the proposed framework. Additionally, applying this methodology results in a small dimension for the system matrix and a reduction in the computational burden. Two test networks, a three-bus system and the IEEE 39-bus system, are used to illustrate and verify the analysis framework.

Abstract-Due to the increasing number of power electronicsbased devices in modern power systems, there are concerns about the impact of higher-frequency interactions on system stability.The typically-used algebraic representation of networks fails to capture these newly emergent issues, but dynamically modelling the entire network with differential equations will lead to a prohibitively long simulation time.This paper proposes a computationally-efficient analysis framework to determine the system loadability in power electronics-rich, large networks with the consideration of multiple types of stabilities.The framework developed in this paper identifies the critical network elements based on the eigenvalue sensitivity in order to exploit a hybrid modelling approach.Within the hybrid model, only the critical network portion is modelled dynamically with the rest of the network represented by algebraic equations.Bifurcation theory is used to simultaneously analyze both smalldisturbance rotor angle and small-disturbance voltage stability.The results obtained show that the high-frequency interactions are accurately captured with the proposed framework.Additionally, applying this methodology results in a small dimension for the system matrix and a reduction in the computational burden.Two test networks, a three-bus system and the IEEE 39-bus system, are used to illustrate and verify the analysis framework.

Index Terms-Bifurcation
with the state variables g(•) : R n x × R n y × Nonlinear vector field associated R + → R n y with the algebraic variables F UTURE power systems will contain increased numbers of power-electronic devices to integrate large quantities of renewable generation.The interactions between these power electronics components and the AC network can result in numerous issues, an example of which is sub-synchronous control interactions (SSCI) [1].These issues and interactions have driven both the small-disturbance angular stability and voltage stability oscillations to a higher frequency range and blurred the boundary between different types of stabilities.This is attracting increasing attention by the research community [1].The work in [2] indicates that the system may encounter small-disturbance rotor angle stability issues before encountering the small-disturbance voltage limit when considering system loadability -particularly when multiple voltage-source converters (VSCs) are introduced into the system.This has been shown to be true by the real power system blackout in Pakistan in 2006 when small-disturbance rotor angle and voltage instability occurred simultaneously [3].Therefore, it is increasingly important to consider both smalldisturbance voltage stability and small-disturbance rotor angle stability when analyzing future power systems.
Previous stability analysis of large-scale power systems has typically focused on low-frequency electromechanical oscillations [4].The higher frequency network electrical dynamics (series and shunt inductance and capacitance, LC) and their 0885-8950 © 2023 IEEE.Personal use is permitted, but republication/redistribution requires IEEE permission.
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
ordinary differential equations (ODEs) are typically neglected since their eigenvalues occur at higher frequencies than are studied.The conventional modeling of AC networks with power flow equations does not include the LC dynamics of the AC network and hence fails to capture the newly emergent power system issues related to power electronic-based devices.It has been noted in [5] that the results from steady-state network modeling can result in errors in system stability determination when compared with dynamic modeling of AC networks using ODEs.Following this idea, the work in [6] has analyzed system stability when the AC network is modeled by sets of differential equations.However, dynamic modeling of the full system with ODEs would result in a prohibitively large number of state variables to be included.This will ultimately lead to a prohibitively long simulation time.It is impractical to model all network dynamics of real large power systems due to the modelling complexity [7].Additionally, the necessity of including the dynamic representation of all power elements remains unclear [8].Hence, there is a need to analyze interactions between VSCs in large-scale systems, but with a reduced computational burden.Previous research employs the network dynamic equivalent model to decrease the computational burden during simulations.The equivalent model is built by dividing the entire network into two sections: the study and external area.The power components and the network in the study area are represented by a detailed mathematical model, while the external area, i.e., the rest of the network, is represented by a simplified model.Hence, the critical characteristics within the study area are maintained, whereas the overall system's order and the network complexity is reduced.The coherency-based technique is employed to identify coherent machines which can be aggregated.Numerous coherency identification approaches have been developed in previous research, including the weak-link method [9], [10], [11], the two-time scale approach [12], [13], [14], [15], linear time simulation [10], [11], [16], and modal analysis [17].The coherent generators identified are then aggregated based on either classical aggregation [10] or adaptive aggregation algorithms [18] to compute the simplified model for the external area.Alternatively, the reduced-order model for the network's external area can be developed using the Ward-type equivalent, based on the reduction of the system admittance matrix.This type of simplified network model is typically used in power flow analysis [19], [20], [21], [22].Apart from network dynamic equivalent models, model order reduction (MOR) is another effective tool to reduce the dimension of the system.A typical application of MOR is in designing power system stabilizers to damp inter-area oscillations in power systems.The main goal of a MOR methodology is to replace the original high-order system model with a reduced-dimensional model, while keeping a good approximation of the dynamic response of the original.MOR techniques have been studied extensively in the literature based on Balanced-truncation [23], [24], [25], [26] and Krylov subspace methodologies [23], [27], [28], [29].However, the simplified network models built based on MOR and the aforementioned equivalent circuit techniques are valid only for one specific operating point.Different loading conditions and network parameter configurations will have a significant impact on the simplified network model constructed.This paper analyses the system equilibrium trajectory as the system loadings continuously vary.Consequently, neither the model reduction technique nor the dynamic equivalent model are suitable for the studies in this paper.
The hybrid model proposed in [7] provides a potential solution for capturing the higher-frequency oscillations caused by the VSCs, when analyzing the small-disturbance stability of large power systems.Unlike the other model reduction techniques mentioned above, hybrid model is suitable for system loadability assessment, in which the system loading is continuously changing.In the hybrid model, the network components in the vicinity of the power converter are modeled with ODEs, with the rest of the network modeled using algebraic equations.This results in a reduction in the dimension of the system Jacobian matrix.Similar modeling ideas have been applied in other areas.For example, the work in [30] divided a system into dynamic and static parts to analyze the dynamics of a VSC-HVDC link.However, the hybrid model adopted in the previous research was designed to preserve the converter-related critical dynamics, with instabilities caused by non-converter equipment neglected.The impact of dynamically modelling different portions of the network is discussed in [8], however, there is still a lack of a methodology that provides a clear and scientific definition of the network portion that should be modelled with ODEs in the hybrid model.
Bifurcation is used to identify quantitative changes in the topological structure of a system [31].Structural change in power systems is commonly related to a change in system stability.Hence, bifurcation theory provides a natural platform for studying power system stability by considering various types of stabilities.Bifurcation theory has been employed to analyze system small-disturbance voltage stability in previous research [32], [33], [34], [35], [36], [37], [38].The work in [37] clearly shows that the system voltage limit at the tip of the PV curve can be represented as a saddle-node bifurcation (SNB), which is a typical type of bifurcation.Bifurcation theory is further adopted in [38] to investigate system voltage collapse issues, in which the system voltage limit is formed as an optimization problem.In addition to studies considering voltage stability, bifurcation has been applied to analyze the small-disturbance rotor angle stability of systems [39], [40], [41], [42], [43].Within those works, the small-disturbance rotor angle instabilities are identified as Hopf-bifurcations (HBs).The results from [39] demonstrate that the system can become unstable after encountering a HB.
Although previous research has investigated system stability with multiple types of stabilities considered, it has typically analyzed each classification of instability individually [44].Unlike conventional analysis techniques, bifurcation theory has the capability to perform system stability assessments with the consideration of multiple instability mechanisms concurrently.With the application of bifurcation theory, small-disturbance rotor angles and small-disturbance voltage instabilities are represented as different types of bifurcation points along the system operation trajectory.Therefore, multiple types of instabilities can be simultaneously analyzed and identified along the system loadability curve using bifurcation theory.
The main contributions of the paper are as follows: r An analysis framework is developed, based on bifurcation theory, to identify system loadability limits when assessing both small-disturbance rotor angle and small-disturbance voltage stability.In combination with hybrid modelling, the developed analysis framework determines which critical components must be modelled with ODEs, with the rest of network represented by conventional algebraic equations.The developed analysis methodology uses a considerably smaller system matrix and has a smaller computational burden compared with the full dynamic model.
r The formulation of the system loadability limit is updated with the consideration of small-disturbance stability and small-disturbance voltage stability concurrently.The necessity of considering both types of stability when calculating system loadability is illustrated by the comparison of different models.
r The instabilities caused by the interactions between VSCs and the AC network, and the interactions between multiple VSCs are identified using the developed analysis framework.Specifically, VSC-related instabilities previously missed, due to the use of traditional network modelling, are now revealed by the developed methodology.

II. BACKGROUND
This section will introduce the bifurcation analysis theories and concepts used within this research.

A. Bifurcation Theory
Power systems are commonly modeled by sets of nonlinear differential-algebraic equations (DAEs), as shown in (1).
In (1), x is a vector of system state variables, for example, the generator rotor angle; y is a vector of system algebraic variables, such as generator voltage magnitudes; μ is the bifurcation parameter which drives the system from one equilibrium point to another.In this study, system loading is selected as a bifurcation parameter.
In regions in which g −1 is not singular, the original DAEs system can be simplified as ODEs, as shown in (2) [45].
(2) Thus, the system's small-disturbance dynamic characteristics can also be indicated by a reduced system matrix as in (3).
In (3), f x is the partial derivation of function f with respect to x, f y represents the partial derivation of function f with respect to y, g x stands for the partial derivation of function g with respect to x , g y is the partial derivation of function g with respect to y.If g −1  y is singular at certain points along the system equilibrium curve, the system will encounter a singularity-induced bifurcation (SIB).However, this paper does not consider SIBs further as few instabilities in real power systems are caused by SIBs.A detailed investigation into SIBs can be found in [46].
The system equilibrium point can be determined by equating the left-hand side of (2) to zero, as in (4).
Suppose (x 0 , μ 0 ) is the system equilibrium point, (4) specifies the location of the equilibrium point x 0 as a function of μ 0 .With μ 0 varying, x 0 changes correspondingly in state space.When the eigenvalues of the system Jacobian matrix have negative real parts, the equilibrium point from ( 4) is asymptotically stable.At the critical point μ c , the system eigenvalues will pass the imaginary axis and the system loses its structural stability.Thus, μ c is then defined as the bifurcation value.
Within this paper, the loading factor is considered as the bifurcation parameter.All loads in the network are varied with respect to this loading factor (i.e., a scaling factor of nominal loading), with each load modelled with constant power factor.This bifurcation configuration represents the system's worst loading scenario and is consistent with past research into system loadability assessments [35], [47], [48].

B. Types of Bifurcation
In nonlinear dynamic systems, bifurcation theory is used to identify the quantitative change in system behaviors with respect to small variations in system parameters.Two types of bifurcation are closely related to power system stability-Saddle node bifurcation (SNB) and Hopf bifurcation (HB) [49].
1) A SNB occurs when the change of bifurcation parameter causes one stable equilibrium and another unstable equilibrium point to merge and the system equilibrium point disappears.At a SNB, the system Jacobian matrix becomes singular and a real mode becomes zero [50].2) A HB occurs when a stable system equilibrium point disappears and emerges as a periodic orbit or vice versa.At a HB, a pair of complex conjugate eigenvalues pass the imaginary axis [50].

III. HYBRID NETWORK MODELING
The hybrid modeling ideas developed in [7] are built upon in this paper.Within the hybrid model, the entire AC network is separated into two parts, namely the dynamic and static parts.Different from the previous work in [7] which only dynamically models the components in the vicinity of the VSC, this research uses the participation factor analysis to identify the key components which require dynamic representation.The details of the proposed analysis framework, including critical component identification, are included in Section IV.This section introduces the proposed hybrid modelling method along with the models of the key power system components.

A. Synchronous Generator
In this work, synchronous generators (SGs) are modeled using the 4th order representation developed in [51] combined with IEEE-Type1 excitation systems [52].The model can be described by a set of DAEs in a dq rotating frame as in (5).
In ( 5), x sg i is the vector of state variables of the ith SG; u sgdi,qi represent the SG termination voltage in dq frame; i sgdi,qi are the currents injected from the SG in dq-frame.

B. Voltage Source Converters
For the VSC modelling, this study uses an ideal lossless average model recommended by CIGRE [53].The VSC is modelled on the dq-axis, with angles provided from the phase-lock loop (PLL).This paper implements the typical vector control structure with inner current control.With the vector control, the VSC can control its active power output or DC voltage by regulating the d-axis current component.The reactive power output or the AC voltage at the point of coupling (PCC) can be manipulated by regulating the q-axis current components of the VSC [54], [55].
The voltage modulation function in typical VSCs in effect provides a feed-forward compensation for all but the largest variations in DC link voltage.Thus unless the VSC is trying to regulate the DC voltage, for the types of study under consideration in this paper, the power transfers in the VSC are naturally independent from the dynamics of the DC-side network, as has been demonstrated in [56].The power output of the VSC is associated with its control schemes and the characteristics of the AC-grid adjacent to the VSC.Furthermore, the research in this paper focuses on interactions between VSCs that are integrated into the same AC network and the impacts of these interactions on the AC-side dynamics of the system.It is assumed that the VSC in this paper is connected to a stiff DC network, which is valid for the following illustrative scenarios: (1) a VSC connected to a battery storage system, and (2) a short VSC-HVDC link.Under such circumstances, the DC-side dynamics of the VSC are assumed to be well controlled by the remote converter station.A comparative study of the dynamic behavior of the system with different DC network representations for a VSC-HVDC link has been included in the Appendix.It is evident that the dynamic response of the VSC-HVDC half-link model, in which one VSC is connected to a constant DC source, shows a good agreement with the results produced using the VSC-HVDC full-link model, in which the entire HVDC model is included.Consequently, the dynamics of the DC network are not included in this paper, and the DC-side network is represented by a constant DC voltage source.This is consistent with past research in [57], [58], [59] when investigating the interactions between VSCs that are coupled through the AC network.Further investigation of the impact of interactions on the system stability with the consideration of DC-side dynamics can be found in [60], [61], [62].
In past studies into interactions between VSCs, a wide range of outer controllers have been considered, including [60], [63], [64].Although P − Q control schemes are typically implemented in the VSC, the results from [57] showed that VSCs with P − V ac are more prone to cause adverse interactions which can potentially deteriorate the system stability.Additionally, this paper focuses on system stability within the AC portion of the system.Consequently, VSCs in this paper adopt feedback P − V ac control in the outer loop.The VSC model can also be represented by a set of DAEs, as shown in (6).
In ( 6), x vsc i is the vector of state variables that represents the ith VSC; u vscdi,qi are the voltages at PCC in dq-frame; i vscdi,qi are the currents flowing out of the VSC in the dq-frame.

C. DQ − dq Transformation
Each SG and VSC is operating under its own local dq-frame, based on either its rotor angle, or the angle provided by a phase-locked loop (PLL).When multiple machines are included in the network, the current injected from each device must be transformed from the local dq frame to a global DQ frame.The transformation is shown in (7).
In (7), α is the reference angle of the global (DQ) frame and β is the reference angle of the local (dq) frame.

D. Static Network Modeling
The static area of the network is modeled by sets of algebraic equations based on a system admittance matrix, as shown in (8).
In (8), Y is the admittance matrix of the static area of the network, v busd,q is the vector of busbar voltages in the static region in the DQ frame; i busd,q is the vector of injection currents to the buses in the static area in the DQ frame.The busbars considered in (8) include the buses within the static region as well as the boundary buses that separate the dynamic and static regions.The current injections considered in (8) include the currents from SGs and VSCs in the static area, and the currents injected from lines that connect dynamic and static areas of the network.The load is considered as a self-admittance at the connected bus and is added to the admittance matrix.The value of the admittance from loads is determined as in (9).
In ( 9), p i and q i are the active and reactive power loads at bus i; p i0 and q i0 are the initial active and reactive power loads.The step size of each iteration in the continuation method represents the increase in the load with respect to the changes in the bifurcation parameter.k in ( 9) is a scaling factor to make Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.the step size smaller and hence improve the result accuracy.A smaller value of k will result in a more accurate result but with a longer simulation time.With the consideration of the simulation time, a value of 0.1 for k is selected in this paper.

E. Static Network Integration
The completed integrated hybrid model is shown in Fig. 1.The upper portion of Fig. 1 depicts the integration of the network's static area.The output (current) of the SGs and VSCs represented in ( 5) and ( 6) match the injection currents in the left-hand side of (8).The right side of (8) provides the voltages which are the input to the VSCs and SGs in ( 5) and (6).
In general, the current injections from the VSCs and SGs are fed into the network, and the network then feeds back the voltages.The overall state-space model of the static part of the grid can be represented as (10).
In (10), i interDi,Qi represents a vector of the currents through the interconnected lines that connect the static and the dynamic areas of the network, in the DQ frame; u bouDi,Qi is a vector of the voltages at the boundaries of the static and the dynamic network portions, in the DQ frame; x sta is the vector of the state variables of the power components included in the static area of the network (for example, the state variables of the SGs in the static network area).

F. Dynamic Network Modeling
The transmission lines in the dynamic part of the network are modeled using a π model [65].The dynamics of the line and the busbars are represented as DAEs, which are shown in (11) and ( 12) respectively. ˙ In ( 11) and ( 12), x line i is the vector of currents passing through the transmission lines; x bus i represents the vector of busbar voltages; u busDi,Qi is the busbar voltage in the DQ frame; i busDi,Qi is the total sum of current injections at the i th bus in the DQ frame.
The load within the dynamic area is considered as a seriesconnected RL impedance.The dynamics of the load employ similar DAEs to those utilized in (11).The value of the impedance that represents the load can be calculated using (13).
In ( 13), r loadi and x loadi are the resistance and reactance values that represent the load at the ith bus.

G. Dynamic Network Integration
The integration of the dynamic areas of a network is shown in the bottom part of Fig. 1.Generally speaking, the dynamics of a busbar work as a junction that interconnects other AC components, including AC transmission lines, loads, VSCs and SGs in the dynamic network areas.The input of the busbar DAEs is the current injections from other components, and its output is the calculated voltages.With the busbar differential equations, all the components in the AC-side of the network are integrated together to form a closed loop.Therefore, the completed statespace model for the dynamic area can be described by (14).
In ( 14), x dyn is the vector of the state variables of the power components included in the dynamic area.The overall hybrid model shown in Fig. 1 is completed by connecting the static area modelled by (10) and the dynamic area modelled by (14).

IV. PROPOSED METHODOLOGY
This paper proposes an analytical framework to assess loadability limits in systems with multiple VSCs considering different instability mechanisms.The method proposed in this paper adopts a hybrid modelling scheme that maintains the accuracy of results within an acceptable error range, whilst significantly reducing the computational burden.The continuation method is adopted to fully track the system loadability trajectory with consideration of multiple bifurcations (i.e., multiple instability types).Further details about the application of the continuation method for loadability study can be found in [66].
The developed approach can be divided into four stages, which are and detailed below: 1) A full dynamic model is created in which the entire network is modeled with ODEs.The critical mode is determined from the full dynamic system matrix based on eigenvalue sensitivity.
2) The key state variables are identified based on participation factor analysis of the critical mode.The key power components (related to key state variables) and their adjacent networks are selected to be modeled dynamically.
3) The hybrid model, which includes dynamic and static areas, is created based on the hybrid modelling approach proposed in Section III. 4) The system loadability curve is calculated by the application of the continuation method using the developed hybrid model.The system loadability assessment is performed using bifurcation analysis.

A. Determination of Critical Modes
The first stage within the proposed method is to determine the critical eigenvalue that first crosses the imaginary axis based on eigenvalue sensitivity.In general, the trajectory of an eigenvalue can be considered as a function of uniform acceleration movement, with the first-and second-order eigenvalue sensitivity mapped to velocity and acceleration concepts from kinematics [67], [68].The eigenvalue with the shortest time, i.e., the eigenvalue expected to reach the imaginary axis first, will be considered as the critical mode.In this study, eigenvalue sensitivity will be directly calculated from an augmented system matrix as the power system is modeled by sets of DAEs rather than ODEs.
1) First-Order Eigenvalue Sensitivity With Respect to a Bifurcation Parameter: To determine the eigenvalue sensitivity, a power system represented by (1) can be rearranged into its state-space form (15).
For a given eigenvalue λ i , the corresponding left and right eigenvectors are w i and v i , therefore ( 16) can be stated.
In ( 16), ν yi and w T yi are the expanded parts of the corresponding eigenvector since an augmented matrix is used.
If the system augmented matrix is Ã = f x f y g x g x , the corresponding left and right eigenvectors of the augmented matrix are w i = w i w yi and ν i = ν i ν yi .Eq. ( 16) can then be simplified as (17) [69].
In (17), Ĩ is an identity matrix in which the first n-diagonal positions have values of 1 and all other positions have values of 0.
The first-order eigenvalue sensitivity with respect to the bifurcation parameter μ is then determined as (18).
2) Second-Order Eigenvalue Sensitivity With Respect to Bifurcation Parameter: It is also possible to extend the sensitivity analysis to determine second-order effects.Differentiating (17) with respect to the bifurcation parameter μ gives (19) [68].
In (19), 19) is singular, the determination of eigenvector sensitivity is not uniquely guaranteed.Additional conditions, where the right eigenvector sensitivity vector is orthogonal to the left eigenvector and the left eigenvector sensitivity vector is orthogonal to the right eigenvector, are considered.Hence two additional equations are added, as shown in (20).
Unique solutions to the problem of left/right eigenvector sensitivity can be found by solving both (19) and (20).Differentiating (18) with respect to the bifurcation parameter μ gives the second-order eigenvalue sensitivity shown in (21) [68].
3) Determination of Key Eigenvalue: The critical condition to maintain system stability within this research is to keep the real part of all eigenvalues negative, i.e., (σ<0), which is also adopted in [2], [67].Consequently, the eigenvalue that first passes through the imaginary axis is considered as the critical mode.Therefore, the real part of eigenvalue sensitivities developed in the previous section, i.e., ∂σ i ∂μ and ∂ 2 σ i ∂μ , will be adopted in this section to identify the critical mode.
From a kinematics point of view, for a given uniform acceleration a and an initial speed V 0 , the time t required to finish a given path S can be expressed as (22).
Similarly, if a system is initially stable, for any eigenvalue moving to the right-hand side, it can be assumed that the movement has uniform acceleration.This assumption has been adopted in previous research [2], [67] when identifying the most critical mode that first passes through the imaginary axis.Although, in a small number of cases, the movements of some eigenvalues are less consistent, the assumption of motion with constant acceleration represents the typical movements of the majority of eigenvalues.Based on more than 10 4 numerical simulations completed in this research, the above assumption is correct for over 90% of cases.The kinematics concepts in Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.(22) can then be represented by eigenvalue sensitivity, with parameters mapping shown in (23).
Hence, ( 22) could be expressed as (24) with respect to eigenvalue sensitivity.
The expected time i.e., the change in the bifurcation parameter μ c_i required to reach the imaginary axis, is therefore determined as (25) [70].
If the system is initially stable, the eigenvalue which has the smallest value of μ c_i is considered as the critical mode.
For initially unstable systems, σ i > 0, the eigenvalue would now be expected to move toward the left-hand side of the complex plane with μ increasing.Hence, the eigenvalue that has the largest of μ c_i will be considered as the critical mode.

B. Identification of Key Power Components
After the critical eigenvalue is identified, the second stage of the proposed methodology is to identify the key power components that require dynamic modeling.Participation factor analysis of the critical mode is used to determine this.A large participation factor value reveals a significant level of state variable participation in the corresponding mode.The definition of the participation factor is shown in (26) [4].
In (26), γ ki represents the relative participation of the kth state variable in the ith mode; v ki is the kth row and ith column of the right eigenvector; w ik is the ith row and kth column of the left eigenvector.
With the critical mode determined as described in the previous section, the participation of each state variable with respect to the critical mode can be calculated.The state variables can then be ranked in descending order of participation factor value.The top state variables with high participation factors are the key state variables.Determining a threshold of what constitutes high participation is a subject that requires knowledge of the system being considered.In this work, all state variable with participation values greater than α times the highest value are considered key state variables, where α is a scalar factor which in this paper is 0.1.
The power components represented by the identified key state variables are considered as key components.Dynamic network representation is used for all portions of the network that include the key components, as well as for the transmission lines and buses immediately adjacent to the key components.For example, if a bus is identified as a key component, then that bus, as well as all the lines and devices connected to that bus, will require dynamic model representation.

C. Construct the Hybrid Model
Once the dynamic network area is identified, the rest of the network will be modeled by algebraic equations and the hybrid model can be created based on the modeling method described in Section III.

D. Complete the Loadability Assessment
The system loadability curve will then be calculated by applying the continuation method to the developed hybrid model [66].The continuation method is selected in this research as it can provide the most detailed information about the system trajectory with respect to multiple instability mechanisms.Other time-efficient techniques, such as optimization-based analysis frameworks as in [35], could also be adopted if desired.Once the system loadability curve is obtained, the loadability assessment can be performed using bifurcation analysis.
It should be noted that the proposed methodology builds the hybrid model from the information provided by the full-dynamic model.This is a typical step in various approaches used to simplify a detailed system model.Take the dynamic reduction program DYNRED, developed under the Electric Power Research Institute's (EPRI) project RP2447, as an example [71].The network reduction program consists of the following three steps to form a reduced system model: (1) Obtain information from the detailed system model; (2) Identify the power components and network area of interest; (3) The network area of interest is modelled with a detailed model, and the rest of the system is represented by a simplified model.The developed analysis framework follows a similar network simplification process to build the hybrid model.Although the proposed approach requires a full-dynamic model at the initial operating point in which the entire system is modelled with ODEs, the proposed framework still reduces the computational burden, as it computes the entire system equilibrium trajectory when performing the system loadability assessment.The system loadability curve within this study typically contains 10 3 ∼10 4 operating points.Therefore, compared to the full-dynamic model, the developed hybrid model with its lower dimensional matrix is a computationally efficient tool for analyzing system loadability trajectory.This is further illustrated and supported by the numerical results in the following section.

V. NUMERICAL RESULTS AND DISCUSSION
The previous sections developed an analysis framework to efficiently assess the system loadability limit considering different types of stabilities.This section applies the developed framework to two test systems, a three-bus system and the IEEE 39 bus system.The three-bus system is used to illustrate the effectiveness of the proposed modeling approach.Following this, the analysis framework is applied to a large system, a modified IEEE 39 bus network, to assess the system loadability limit considering interactions between VSCs and the network.All the simulations are performed on a Lenovo ThinkStation P520c with Intel Xeon W-2123 3.60 GHz CPU and 32 GB RAM.

A. Case-1: Three-Bus System
This section utilizes a modified three-bus system, adopted from [37], to illustrate the proposed modeling approach.The modified three-bus system is illustrated in Fig. 2. Compared with [37], the SG at bus 1 is replaced with a VSC operating in P-V ac mode [72].The various colors utilized in Fig. 2 distinguish the different modelling representations used.Details of the different models will be explained in the following sections.
To illustrate the effectiveness of the proposed modeling approach, three different models are created.The models used are a steady-state model, a full-dynamic model and a hybrid model.Within the steady-state model, all transmission lines are represented by algebraic equations with their dynamics ignored.The full dynamic model has a dynamic representation of every network component with ODEs.In the hybrid model, the key network components are modeled with ODEs and the rest of the network is represented with algebraic equations.Firstly, the number of state variables in the three models is calculated to illustrate the difference in the expected computational burden imposed by the models.This will be further considered by comparing the simulation times taken when the different models are used.Then system equilibrium trajectories of the different models are calculated to show that the same operating points are produced, i.e., the same trajectories are produced by the different models.Finally, a system stability assessment, i.e., loadability calculation, along the equilibrium trajectory will be completed using bifurcation theory.Unlike the equilibrium trajectories, the system loadability assessments from the three models may vary, if different bifurcations (i.e., different instability phenomena) are captured by the different models.
In this section, three different models are created.SS-Model is a steady-state model in which both lines 1-2 and 2-3 are modeled by algebraic equations.Hy-Model represents a hybrid model.In Hy-Model, the dynamics of line 1-2 are modeled with ODEs and the rest of the network has a static representation with algebraic equations.Referring to Fig. 2, the part of the network represented dynamically in Hy-Model, i.e., line 1-2, is colored red.It should be noted that the three-bus test system is used to illustrate the effectiveness of the hybrid modeling approach on a small (and easy to understand) system, rather than the effectiveness of the entire assessment framework.The hybrid model in Case-1, therefore, is created deliberately with manual selection of the dynamic network region using the framework, rather than developed by the proposed assessment framework using the set threshold to automatically identify key network components.The dynamic area of the network in Hy-Model is surrounded by a dashed line in Fig. 2. Finally, Dyn-Model is a full dynamic model, in which both lines 1-2 and 2-3 are dynamically modeled with ODEs.
Table I summarizes the overall number of state variables in each model.It reveals that SS-Model consists of 12 state variables, the smallest of the three models.When the entire network is modeled with ODEs (Dyn-Model), the state space model has 18 state variables.This number is reduced to 16 when the hybrid model (Hy-Model) is used.The reduced number of state variables in Hy-Model indicates a lower-dimension system Jacobian matrix and hence an improved computational efficiency compared with Dyn-Model.
The simulation time taken by the three different models to perform a loadability assessment is summarized in Table II.Table II reveals that the steady-state model requires the shortest simulation time, which is 3 s.Dyn-Model takes the longest computation time, which is 5 s.The required time is reduced to 10 s when the hybrid model is used and, hence, this hybrid model results in a reduction in computational burden.In combination with results in Table II, this shows that a model with a larger number of state variables comes with a longer simulation time and, hence, imposes a heavier computational burden.It should be noted that in Case-1, a small sized power system is used.Hence, the difference in time taken by models is not significant.In Case-2 (shown later), for a large system, the dynamic modeling region in the hybrid model is a small proportion of the entire network, and therefore, the difference between the simulation times taken by models is greater.
The system equilibrium trajectories of the different models when the impedance of Line 1-2 is 0.5 p.u. are shown in Fig. 3.The results are obtained using the continuation method.The horizontal axis is the loading factor of the system which is also the bifurcation parameter in this paper, and the vertical axis represents the voltage magnitude at Bus 2. As is shown by Fig. 3, the system trajectories using different models are visually identical.Considering the trajectory obtained from the  full dynamic model as the standard result, the root-mean-square deviations of the curve derived from the steady-state model and the hybrid model are 4.57 × 10 −5 and 5.78×10 −5 respectively.The acceptable error tolerance used in this paper is based on the typical tolerances used in the continuation method, typically around 10 -2 [35].Such small errors from the steady-state model (4.57× 10 −5 ) and the hybrid model (5.78 × 10 −5 ) could easily be due to the application of the continuation method process itself.This illustrates that the different models produce the same operating points as the loading increases.Based on the bifurcation analysis, multiple instability mechanisms along the system trajectory will be identified as SNB and HB.The system maximum loading points identified as SNBs for the three models are, therefore, essentially identical.The loadability assessment along system trajectories, including SNBs and HBs, will be investigated further in the following sections.The case considered will illustrate the differences in the loadability assessments when using different models.
In order to confirm the assumption in Section III-B that the DC-side dynamics may be neglected, the loadability assessments shown in Table III have been conducted using both the HVDC full-link and the half-link models.The results have demonstrated that the bifurcation locations produced by the VSC-HVDC full-link and half-link models are identical up to five decimal places.This illustrates the effectiveness of the halflink model in terms of producing accurate loadability results.All three system models, including SS-Model, Hy-model and Dyn-Model, in the following sections will adopt the HVDC half-link model.
Table III summarizes the loadability assessment along the system trajectories.The results from Dyn-Model reveal that the system is initially unstable due to a HB and becomes stable as the load increases.As the load is further increased, the system finally encounters a SNB.The previous discussion about system equilibrium trajectories using different models suggests that the SNB values from the three models would be the same.This is seen in Table III, where the SNB values from the three models are all equal to 8.17 p.u. to three significant figures.Table III also details the HB values obtained.It can be seen that SS-Model fails to capture the occurrence of a HB along the system trajectory.The other two models have similar results for the HB location.Hy-Model shows that the instability will end at 8.10 p.u., and this is similar to the HB value from Dyn-Model, which is 8.14 p.u. Considering the assessment from Dyn-Model as the standard result, the relative error of assessment of the HB from Hy-Model is 4.91 × 10 −3 , which is smaller than the acceptable error range of 10 −2 .
Fig. 4 shows the participation factor analysis of the HBs for the different models.The steady-state model is not included in Fig. 4, because it fails to capture the HBs.It can be seen that the state variables of line 1-2 (i 12x , i 12y ) and state variables belonging to the VSC (e d , e q , x P ) participate most in the critical modes from both models.The SS-Model (not shown in Fig. 4) does not include the dynamic representation of the critical line 1-2.This is why it fails to capture the HBs, resulting in an inaccurate stability assessment.
According to the results obtained from Dyn-Model, the state variables e d , e q participate most with participation values of 0.200 and 0.212 respectively.The state variable x P (representing the active power control of the VSC) and the state variables of line 1-2 (i 12x , i 12y ) also participate strongly in the critical mode, with values of 0.125, 0.089 and 0.091 respectively.Hy-Model produces similar results and it can be established that the fulldynamic and hybrid models share a similar participation factor value distribution and, hence, both accurately assess the system stability along the loadability trajectories.
This case clearly illustrates that conventional modelling of the AC network with algebraic equations fails to capture certain types of instabilities along the system trajectory.The steady-state model, therefore, may not always be the appropriate model for loadability assessment.The full dynamic model, on the other hand, always calculates the system loadability accurately but requires a high-dimensional system matrix, which leads to increased computational burdens.The hybrid model combines the advantages of those two models.It can assess the system loadability accurately but with a reduced computational burden.Hence, the application of the hybrid model is recommended for Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.Fig. 5. Modified IEEE-39 bus system for Case-2-Scenario-1.(key elements colored in red and the dynamic area surrounded with a dashed line).system loadability analysis.Additionally, it is also shown that a system may become oscillatory unstable (identified as HB) before reaching the static maximum load point (identified as SNB).Hence, there is a need to consider both small-disturbance voltage and small-disturbance rotor angle stability simultaneously when assessing system loadability limits.

B. Case-2: 39-Bus System
This section uses a modified IEEE 39-bus system to demonstrate the effectiveness of the proposed methodology.The scenarios considered represent HBs caused by 1) network dynamics, 2) interaction between VSC and network and 3) interaction between multiple VSCs.A similar analysis process to Case-1 (the three-bus system discussed in the previous section) is used to compare the results from the different models.
The system loadability results are compared using the three different models created, i.e., SS-Model, Hy-Model, and Dyn-Model.It should be noted that for this case, the hybrid model is determined using the developed assessment framework (and not with manual selection of dynamic network areas as was used for Case 1).With the developed assessment framework, the interactions between the VSCs and the network can be captured as bifurcations along the system loadability trajectory with a reduced computation burden.
1) Scenario-1: One VSC: A modified IEEE 39-bus system is used in Scenario-1 and is shown in Fig. 5. Fig. 5 shows that the original SG at bus 32 is replaced by a VSC which operates in P-V ac mode.The dashed green line which encapsulates the red network elements in Fig. 5 distinguishes the dynamic region in the hybrid model.The details of the hybrid model will be explained in the following section.This study scenario will illustrate the effectiveness of the developed methodology in a large power system.For brevity, the participation of each state variable in critical mode is not shown.According to the participation factor analysis, the critical mode is dominated by the state variables belonging to G31; the VSC; Buses 5,6,7,8,10,11 and 13; the load at Bus 7, 31; and Lines 4-5, 5-8, 5-6, 6-7, 7-8, 6-31, 6-11, 10-11 and 10-13.Those AC components are considered as key components and colored red.Following the proposed assessment framework, the dynamic area of the network includes key network components and their adjacent components.Referring to Fig. 5, the dynamic area of the network is surrounded by a green dashed line.According to the developed framework, the network elements within the dynamic area are modeled with ODEs and the rest of the network is represented by algebraic equations.Table IV summarizes the number of state variables in the different models.The full dynamic model (Dyn-Model) consists of 277 state variables, which is roughly twice the number in the hybrid model The simulation time taken to perform a loadability assessment with the different models is shown in Table V. Dyn-Model consists of the largest number of state variables and, hence, requires the longest computation time (684 s).The simulation time is significantly reduced, to 246 s, when the developed hybrid model is used.This further highlights the impact that reducing the number of state variables has on the computational burden.
Table VI summarizes the bifurcations along the system loadability curves for Scenario-1 using the three models.In Scenario-1, the system is stable at the beginning and all eigenvalues have negative real parts.As the load increases, the system will show oscillatory instability, which is represented as HBs.Eventually, the system reaches a SNB as the loading further increases.Table VI shows that the SNB values are the same from the different Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.Fig. 6.Modified IEEE-39 bus system for Case-2-Scenario-2.(key elements colored in red and the dynamic area surrounded with a dashed line).models, at 2.14 p.u.However, the HB value from the SS-Model is significantly different from the results using the other two models.As is illustrated in Table VI, Hy-Model and Dyn-Model have similar HB values of 1.67 p.u. and 1.68 p.u. respectively.It can also be seen that the SS-Model has a different HB value at 2.11 p.u.. Taking the results from Dyn-Model as the standard assessment, the relative error of the HB value from Hy-Model is 5.95×10 −3 , which is smaller than the acceptable error range of 10 −2 .However, the relative error for SS-Model is 0.26, and this would result in a very inaccurate stability assessment.This error is due to the static (algebraic) representation of the key network components (i.e., the network around the VSC and G31) in the SS-Model.
From these results, it can be concluded that the proposed hybrid model can precisely assess the system loadability with a reduced number of state variables.This scenario again shows the necessity of considering multiple types of stabilities simultaneously when assessing the system loadability, as a system may become unstable via a small-disturbance angle instability before meeting with a small-disturbance voltage instability.
2) Scenario-2: Interactions Between VSC and Network: In Scenario-2 of Case-2, one VSC is connected to Bus 38, replacing the original G38.This scenario is included to demonstrate that the developed analysis framework accurately captures the instability due to interactions between the VSCs and the network along the system loadability curve.The modified IEEE-39 bus system is shown in Fig. 6.The previous work in [73] indicates that under a weak grid condition, the interactions between the VSC and the AC network will be significant and may further degrade the system stability.The network impedances in the vicinity of the VSC (i.e., impedances of transmission Lines 26-27, 26-28, 26-29, 28-29 and 28-38) are increased by a factor of 1.89; the loads at Buses 26, 27, 28 and 29 are increased by a factor of 1.8; and the power output of G38 is multiplied by a factor of 1.2 to exacerbate the instability due to the interaction between the VSC and the network.
According to the participation factor analysis, the critical mode represents the interaction between VSC and the adjacent network.Correspondingly, the state variables belonging to the VSC and its adjacent network have high participation in the critical mode.The key elements identified by the proposed framework are the VSC and Lines 26-27, 26-28, 26-29 and 28-38.In Fig. 6, the key network components are highlighted with red.The dynamic modelling area of the network in the hybrid model is surrounded by a dashed line.Hy-Model consists of 103 state variables compared to 71 for SS-Model and 277 for Dyn-Model.
Table VII shows the bifurcations along the system loadability trajectory from the three models.The system in Scenario-2 is initially unstable and becomes stable as the load increases.Similar to the results of SNBs in the previous scenarios, all three models produce the same SNB value, which is 1.15 p.u.The results from Dyn-Model presented in Table VII reveal that the system becomes stable when the system loading increases to 0.80 p.u. Hy-Model also has the same HB value.The SS-Model fails to capture HBs along the system loadability curve due to the missing key network dynamics.Scenario-2 again illustrates the effectiveness of the proposed methodology in a large system.Additionally, it can be concluded from this scenario that the proposed framework can be used to reveal the instability caused by the interactions between the VSC and the network.
3) Scenario-3: Interactions Between Multiple VSCs: In Scenario-3 of Case-2, two VSCs are integrated into the IEEE 39-bus system.As is shown in Fig. 7, the modified system, two VSCs are connected at Buses 30 and 37 replacing G30 and G37.To exacerbate the interactions between the VSCs (in a similar way to Scenario-2), the impedances around the VSCs are increased in this scenario.Specifically, the impedances of Lines 2-30, 2-25 and 25-37 are multiplied by a factor of 4, the power output of VSC 30 and 37 is multiplied by a factor of 2.5 and 1.2, and the power output of G38 is decreased by a factor of 0.5.
The participation factor analysis indicates that the state variables belonging to the VSCs and the network around the VSC have high participation in the critical mode.The key network components identified are the VSCs at Buses 30 and 37, Lines 2-30 and 25-37, Bus 2, and the load at Bus 25 -all highlighted in red in Fig. 7. Hy-Model consists of 96 state variables compared to 74 for SS-Model and 280 for Dyn-Model.The bifurcations along the system loadability trajectory from different models are shown in Table VIII.Again, the SNB values from the three models are identical, at 2.08 p.u.The HB values calculated from Hy-Model and Dyn-Model are similar, at 1.23 p.u. and 1.22 p.u. respectively.However, the HB value from SS-Model is 1.51 p.u, which is significantly different from the results from the other two models.
The analytical results in this scenario are further verified using time-domain simulations.The time-domain simulation results using the different models for a value of λ = 1.3 p.u. have been produced and are shown in Fig. 8. Based on the analytical results shown in Table VIII, the system is unstable when the Dyn-model and Hy-model are used, and stable when the SS-model is adopted -as the HB is incorrectly located at a higher loading.A small disturbance, in which the power reference of G32 increases from 1.0 p.u. to 1.1 p.u., is injected into the system at a time of t = 10 s.The rotor speed of G32 (ω G32 ) is plotted for the three different models.The results from the Dyn-and Hy-models show that the system becomes unstable after the small disturbance as illustrated in Fig. 8(a) and (b).However, if the SS-model is used, the rotor speed of G32 quickly stabilizes after the change in the power reference (Fig. 8(c)).It can be concluded that the time-domain simulation results verify the previous analytical results.Therefore, this demonstrates that the instability due to the interaction between VSCs could occur along the system loadability curve, and that this phenomenon could be captured by the developed hybrid model.However, the conventional steady-state model (i.e., the SS-model) fails to capture this instability due to the incorrect HB location.This could lead to significant errors in loadability assessment, concluding that the stability limits are higher than in practice.

VI. CONCLUSION
This paper proposes an assessment framework for analyzing system loadability efficiently with the consideration of both small-disturbance voltage stability and small-disturbance stability.The proposed method is developed based on bifurcation theory and combines several techniques, including eigenvalue sensitivity and participation factor analysis.With those techniques and combined with hybrid modelling, the developed analysis framework determines the key network components that need to be modelled dynamically, whilst the rest of the network is represented by conventional algebraic equations.Three models, i.e., a steady-state model, a hybrid model and a full dynamic model, have been developed for different scenarios and applied to two test systems to show the effectiveness of the proposed assessment framework.
The results have shown that instabilities caused by different mechanisms including -poor-damped network oscillations, interactions between a VSC and the network and the interactions between VSCs -may all remain undetected if a conventional steady-state model is used.It has been demonstrated that the proposed analysis framework accurately captures these instabilities but with a reduced-dimensional matrix and a shorter simulation time when compared with the full dynamic model.The proposed analysis framework could, therefore, benefit power system operators in identifying the system loadability limit precisely with a lower computational burden.Additionally, it has been clearly demonstrated that it is necessary to consider both small-disturbance angle and voltage stability when calculating the system loadability.
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

APPENDIX
Consider a typical two-terminal VSC-HVDC system connecting two AC systems as is illustrated in Fig. 9(a) [53], [56].Such a system consists of two VSCs, two AC systems and a DC system.The external AC networks are commonly represented by a voltage source (for VSC 1 ) or a synchronous generator (for VSC 2 ) in series with the grid equivalent impedance.In a typical configuration of the VSC-HVDC link, the power flows from VSC 1 to VSC 2 [56].In such a VSC-HVDC system, it can be assumed that VSC 1 operates with the V dc − Q outer controller and VSC 2 operates with the P − V ac outer controller.This appendix compares the dynamic response of power system with a VSC-HVDC link using different DC-side network representations.Two different test network models employed in this appendix are shown in Fig. 9, which are (a) a full-link model, where detailed mathematical model of the entire VSC-HVDC link is included and (b) a half-link model, where the DC-side dynamics of VSC 2 are neglected and simply represented by a constant DC voltage source.
To simulate the scenario in which VSC 2 is connected to a stiff DC network, all control parameters of both VSCs, particularly the DC voltage outer loop of VSC 1 , are tuned properly so that the DC voltage in the full-link model is well controlled by the VSC 1 station.It should be noted that the same system parameter setting shown in Table IX is used in both models to ensure the comparability of results.All the EMT simulations are performed using DIgSILENT PowerFactory 2022.
The EMT simulation results obtained when using different HVDC system representations are illustrated in Fig. 10.A small disturbance, in which the power reference of VSC 2 increases from 0.30 p.u. to 0.33 p.u., is injected into the system at time t = 5 s.Fig. 10 shows (a) the power output of VSC 2 at the PCC and (b) Phase A currents flowing out of the VSC for both models.The results using the VSC-HVDC full-link and half-link models are represented by a solid orange line and a dashed black line, respectively.
It is clear from Fig. 10(a) that the VSC output power from both models increases from 300 MW to 330 MW after the change in the reference power setting.Fig. 10(a) has clearly illustrated that the power output of VSC 2 calculated from the half-link model accurately overlaps with that produced using the full-link model, except for oscillations which last a few milliseconds after the disturbance is injected.Within the full-link model, the DC voltage is controlled by VSC 1 at the other end of the HVDC link, while in the half-link model, the VSC is directly connected to a constant DC voltage source.Hence, the full-link model simulates the power system transient behavior slightly better than the half-link model.Although the difference between initial dynamics using two different models exists, it is not critical in this study as the paper focuses on system behaviors in a quasi-steady-state timeframe.Specifically, the system stability information at steady state is of more interests rather than the transient behavior after the disturbance is injected.This can be confirmed by the identical loadability assessment using HVDC full-link and half-link models, as included in Section V-A.The relative difference between the peak values of the two models with respect to power output from the full-link model is 0.15%.This clearly demonstrates a good agreement between the results when the two models are used.Additionally, the high degree of accuracy of results from the half-link model can be further seen from the almost perfectly overlapping phase A current curves in Fig. 10(b).Taking the curve produced from the full-link model as a baseline, the root-mean-square error (RMSE) of the current output using the half-link model 2.53×10 −8 .It is evident from Fig. 10 that the half-link model can thus be used to provide system dynamic responses that are suitably accurate compared to those obtained from modelling the full VSC-HVDC link when analyzing AC-side stability.Therefore, these EMT simulation results confirm the assumption in Section III-B that the DC-side dynamics of the VSC can be neglected when the VSC is connected to a stiff DC network.Based on this, the half-link model is used in this study to analyze the impact of interactions between VSCs on the network AC-side dynamics.
To further demonstrate the validity of the assumption in Section III-B, the dynamic response of the rotor angular speed of G32 (ω G32 ) in Fig. 8(a) is reproduced with the VSC-HVDC represented by a full-link model and a half-link model and is shown in Fig. 11.The results using the full-link and the half-link model are represented by a solid orange line and a dashed black line, respectively.As is shown in Fig. 11, the results using those two models are visually overlapping.Considering the results obtained from the full-link as the standard result, the RMSE of the curve obtained using the half-link model is 2.09×10 −7 .Such a small difference could easily be due to the computational process itself.These time-domain simulation results have further demonstrated the suitability of the HVDC half-link model for the loadability studies discussed within this paper.

Fig. 2 .
Fig. 2. Modified three-bus test system.(key elements colored in red and the dynamic area surrounded with a dashed curve).

Fig. 4 .
Fig. 4. Participation of state variables in critical mode.

Fig. 7 .
Fig. 7. Modified IEEE-39 bus system for Case2-Scenario-3.(key elements colored in red and the dynamic area surrounded with a dashed line).
A Framework for Analyzing System Loadability With Multiple VSCs Using a Hybrid Model Youhong Chen , Student Member, IEEE, Robin Preece , Senior Member, IEEE, and Mike Barnes , Senior Member, IEEE

, hybrid model, loadability, interactions, VSCs.
i Vector of current u ∈ R n v Vector of voltage B. Function

TABLE I STATE
VARIABLES FOR THE DIFFERENT MODELS IN CASE-1

TABLE IV STATE
VARIABLES IN THE DIFFERENT MODELS IN CASE-2-SCENARIO-1

TABLE V SIMULATION
TIME IN DIFFERENT MODELS IN CASE-2

TABLE VII BIFURCATIONS
ALONG THE SYSTEM LOADABILITY CURVES FROM DIFFERENT MODELS IN CASE-2 SCENARIO-2