Application of the Trajectory Sensitivity Theory to Small Signal Stability Analysis

Abstract

The security assessment of power systems represents one of the principal studies that must be carried out in energy control centers. In this context, small-signal stability analysis is very important to determine the corresponding control strategies to improve security under stressed operating conditions of power systems. This chapter details a practical approach for assessing the stability of power system’s equilibrium points in real time based on the concept of trajectory sensitivity theory. This approach provides complementary information to that given by selective modal analysis: it determines how the state variables linked with the critical eigenvalues are affected by the system’s parameters and also determines the way of judging how the system’s parameters affect the oscillatory behavior of a power system. The WSCC 9-bus and a 190-buses equivalent system of the Mexican power system are used to demonstrate the generality of the approach as well as how its application in energy management systems is suitable for power system operation and control.

1. Introduction

Small-signal stability (SSS) is the ability of a power system to maintain synchronism when subjected to small disturbances such as small load and/or generation changes [1]. The analysis of SSS consists of assessing the stability of an equilibrium point (EP), as well as determining the most influential state variables in the stability of the operating point. For small enough disturbances, the system behavior can be studied via the theory of linear systems around an equilibrium point [2]. The stability of an EP is assessed by eigenvalue analysis (eigenanalysis) according to the Lyapunov criterion [3], which states that an EP will be stable in the small-signal sense, if all system eigenvalues of the system matrix are located on the left side of the complex plane. On the contrary, the EP will be unstable if at least one eigenvalue is located on the right side of the imaginary axis. In this context, the resulting dominant eigenvalue from the eigen-analysis is called the critical eigenvalue, and its association to the state variables is investigated by selective modal analysis (SMA) [4, 5]. Based on the participation factors analysis (PFA), the SMA provides those state variables having the highest influence in the EP stability by means of their coupling to the critical eigenvalue [6, 7]. Since the system eigenvalues are directly related to its dynamic performance, different forms of instability in a power system can be studied by means of well-defined structures of eigenvalues which are called bifurcations.

The theory of bifurcations is a powerful mathematical tool based on eigen-analysis to assess the stability of EPs in nonlinear systems [8]. This theory consists of searching for specific eigenvalue structures associated to different instability forms that appear on power systems [9]. One of the most common local bifurcations that can appear in the power system operation is the Hopf bifurcation (HB), which occurs when the system matrix contains a pair of purely imaginary eigenvalues causing undamped oscillatory behavior [10, 11, 12, 13]. Any parameter variation in the system may result in complicated behavior until the system stability changes. This point, where the stability changes, is defined as a bifurcation point. In this chapter the system loads are changed in order to analyze the stability of the EPs. The maximum load in a specific direction that a power system can provide before the appearance of a bifurcation point establishes the loading limit in that direction. Thus, the loading limit is directly associated with the stability margin of the system. The critical eigenvalue of an EP is used as an index of the stability margin. After a load change, an eigen-analysis permits us to assess the stability margin closeness. A small margin indicates closeness to a bifurcation point (instability).

SSS analysis is very important to determine the corresponding control strategies to improve security under stressed operating conditions of power systems. Control strategies employed in electric power systems are usually tested by means of an assessment of the stability improvement. Thus, the influence of parameters and components in the EP stability provides an insight to achieve the best control. In this context, the participation factors let us know the highest association between state variables and the critical eigenvalue dominating the EP stability [6, 7]. In this way, PFA allows the selection via the associated states to the critical eigenvalue of those components that will provide the best control in EP stability. Although the PFA selects the most sensitive states in the EP stability, it is not possible to identify in a direct form the most influential parameters, e.g., those most sensitive loads influencing the stability. In order to achieve this, the PFA must be combined with other methods providing parameter sensitivity features. In [14, 15] the authors combined the PFA and modal controllability of the weak damping oscillatory modes to obtain an optimal location of static VAR compensators (SVCs). In [16] the authors presented an approach to examine the effect of loads in the system stability by using participation factors and mode shape analysis. In Ref. [7], both a voltage stability analysis and a low-frequency oscillation analysis were performed by using the SSS analysis in the EPs. As the system faces increased loading conditions, bifurcation points appear, and the participating generators are identified by means of the most associated states obtained from the SMA.

This chapter presents an alternative method based on the trajectory sensitivity (TS) theory [17, 18] to investigate the stability of the EPs by using a time-domain simulation. In this approach which will be referred to as selective modal analysis and trajectory sensitivities (SMA-TS), the TS were computed with respect to the load parameters; therefore, besides the stability assessment and the participating states, this approach also has the ability to identify those most sensitive load parameters influencing the critical states. The stability assessment consists of just examining the TS oscillations. However, the stability analysis of the EP is never perturbed. The TS oscillations required for the stability assessment are produced by means of the initial condition values selected for the sensitivity variables. The SMA-TS approach uses an index of sensitivity quantification, which facilitates identifying the influence of the system loads around the HB points. The index allows to rank the power system loads in order to predict the most critical loading directions toward a B point. Its application is suitable to monitoring in real time the power system operation and improve the security. The results of the study cases of two real systems, i.e., WSCC 9-bus and an equivalent of the Mexican power system of 190-buses, are used to show the performance and applicability of the proposed approach.

2. Trajectory sensitivity theory

An electric power system can be represented analytically by a set of differential-algebraic equations (DAEs), as given by (1), where $x$ is a n-dimensional vector of dynamic state variables with initial conditions $x\left(t_0\right)=x_0$ , $y$ is a m-dimensional vector of instantaneous state (algebraic) variables, (usually the real and imaginary parts or the magnitudes and phase angles of the complex node voltages) with initial conditions $y\left(t_0\right)=y_0$ , and $\beta$ is a set of time-invariant parameters of the system. The dynamics of the equipment, e.g., generators and controls, is explicitly modeled by the set of differential equations through the function $f\left(\cdot \right)$ . The set of algebraic equations $0=g\left(\cdot \right)$ represents the stator algebraic equations and mismatch power flow equations at each node:

$$\begin{array}{l}\dot{x}=f\left(x,y,\beta \right)f:{\mathfrak{R}}^{n+m+p}\to {\mathfrak{R}}^n\\ 0=g\left(x,y,\beta \right)g:{\mathfrak{R}}^{n+m+p}\to {\mathfrak{R}}^m\\ x\in X\subset {\mathfrak{R}}^{n}y\in Y\subset {\mathfrak{R}}^m\beta \in \mathit{\beta }\subset {\mathfrak{R}}^p.\end{array}$$ E1

2.1. Analytical formulation

Let ${\beta }_0$ be the nominal values of $\beta$ , and assume that the nominal set of DAEs $\dot{x}=f\left(x,y,{\beta }_0\right)$ , $0=g\left(x,y,{\beta }_0\right)$ has a unique nominal trajectory solution $x\left(t,x_0,y_0,{\beta }_0\right)$ and $y\left(t,x_0,y_0,{\beta }_0\right)$ over $t\in \left[t_0,t_{\mathit{end}}\right]$ , where $t_0$ and $t_{\mathit{end}}$ are the initial and final times, respectively, of the study time period. Thus, for all $\beta$ sufficiently close to ${\beta }_0$ , the set of DAEs (1) has a unique perturbed trajectory solution $x\left(t,x_0,y_0,\beta \right)$ and $y\left(t,x_0,y_0,\beta \right)$ over $t\in \left[t_0,t_{\mathit{end}}\right]$ that is close to the nominal trajectory solution. This perturbed solution is given by (2) and (3) [19, 20]:

$$x\left(\cdot \right)=x_0+{\int }_{t_0}^{t_{\mathit{end}}}f\left(x\left(\cdot \right),y\left(\cdot \right),\beta \right)\mathit{ds}$$ E2

$$0=g\left(x\left(\cdot \right),y\left(\cdot \right),\beta \right).$$ E3

The sensitivities of the dynamic and algebraic state vectors with respect to a chosen system’s parameter, $x_{\beta }=\partial x\left(\cdot \right)/\partial \beta$ and $y_{\beta }=\partial y\left(\cdot \right)/\partial \beta$ , at a time $t$ along the trajectory are obtained from (4) and (5), which in turn are obtained from the partial derivative of (2) and (3) with respect to $\beta$ :

$$\frac{\partial x\left(\cdot \right)}{\partial \beta }={\int }_{t_0}^{t_{\mathit{end}}}\left(\frac{\partial f\left(\cdot \right)}{\partial x}\frac{\partial x}{\partial \beta }+\frac{\partial f\left(\cdot \right)}{\partial y}\frac{\partial y}{\partial \beta }+\frac{\partial f\left(\cdot \right)}{\partial \beta }\right)\mathit{ds}$$ E4

$$0=\frac{\partial g\left(\cdot \right)}{\partial x}\frac{\partial x}{\partial \beta }+\frac{\partial g\left(\cdot \right)}{\partial y}\frac{\partial y}{\partial \beta }+\frac{\partial g\left(\cdot \right)}{\partial \beta }.$$ E5

Lastly, the smooth evolution of the sensitivities along the trajectory (6) and (7) is obtained by differentiating (4) and (5) with respect to $t$ :

$${\dot{x}}_{\beta }=\frac{\partial f\left(\cdot \right)}{\partial x}\frac{\partial x}{\partial \beta }+\frac{\partial f\left(\cdot \right)}{\partial y}\frac{\partial y}{\partial \beta }+\frac{\partial f\left(\cdot \right)}{\partial \beta }\equiv f_x{x}_{\beta }+f_y{y}_{\beta }+f_{\beta };x_{\beta }\left(t_0\right)=0$$ E6

$$0=\frac{\partial g\left(\cdot \right)}{\partial x}\frac{\partial x}{\partial \beta }+\frac{\partial g\left(\cdot \right)}{\partial y}\frac{\partial y}{\partial \beta }+\frac{\partial g\left(\cdot \right)}{\partial \beta }\equiv g_x{x}_{\beta }+g_y{y}_{\beta }+g_{\beta };y_{\beta }\left(t_0\right)=0$$ E7

where $f_x$ , $f_y$ , $f_{\beta }$ , $g_x$ , $g_y$ , and $g_{\beta }$ are time-varying matrices computed along the system trajectories.

3. Trajectory sensitivity analysis

3.1. Sensitivity discretization

TS computation is obtained by means of the sequential solution of the nonlinear DAE system (1) and the linear time-varying DAE system (6) and (7). By applying the trapezoidal rule to algebraize the differential equations, both DAE systems are converted into the following systems of algebraic difference Eqs. (8)–(11).

$$F_1\left(\cdot \right)=x^{k+1}-x^k-\frac{\mathit{\Delta t}}{2}\left(f^{k+1}+f^k\right)=0$$ E8

$$F_2\left(\cdot \right)=g^{k+1}=0$$ E9

$$F_3\left(\cdot \right)=x_{\beta }^{k+1}-x_{\beta }^k-\frac{\mathit{\Delta t}}{2}\left(f_x^{k+1}{x}_{\beta }^{k+1}+f_y^{k+1}{y}_{\beta }^{k+1}+f_{\beta }^{k+1}+f_x^k{x}_{\beta }^k+f_y^k{y}_{\beta }^k+f_{\beta }^k\right)=0$$ E10

$$F_4\left(\cdot \right)=g_x^{k+1}{x}_{\beta }^{k+1}+g_y^{k+1}{y}_{\beta }^{k+1}+g_{\beta }^{k+1}=0$$ E11

where $\mathit{\Delta t}$ is the integration time step and the superscript $k$ is an index for the time instant $t_k$ at which variables and functions are evaluated, e.g., $x^k=x\left(t_k\right)$ and $f^k=f\left(x^k,y^k\right)$ .

3.2. Linear sensitivity computation

Once algebraized both DAE systems (8)–(11), the Newton-Raphson (NR) algorithm is used to provide an approximate solution of the algebraized nonlinear systems (8) and (9). In this case, the resulting linearized system $J^i\mathrm{\Delta }{X}^i=-F{\left(X\right)}^i$ , whose representation in expanded form is given in (12), provides the approximate solution $\mathrm{\Delta }{X}^i$ = ${\left[\mathrm{\Delta }{x}^k\mathrm{\Delta }{y}^k\right]}^T$ , which is updated at each NR $i$ th iteration, i.e., ${\left[\begin{array}{ll}{x}^{k+1}=x^k+\mathrm{\Delta }{x}^k& y^{k+1}=y^k+\mathrm{\Delta }{y}^k\end{array}\right]}^T$ , until a selected convergence criterion is satisfied. Note also that $J$ is the Jacobian matrix resulting of the linearization around an EP.

The initial guess ${\left[\begin{array}{ll}{x}_0^{k+1}=x^k& y_0^{k+1}=y^k\end{array}\right]}^T$ is used to start the NR algorithm from given values ${\left[\begin{array}{ll}{x}^k& y^k\end{array}\right]}^T$ :

$${\displaystyle \underset{J^i}{\underbrace{{\left[\begin{array}{ll}I-\frac{\mathit{\Delta t}}{2}{f}_x^{k+1}& -\frac{\mathit{\Delta t}}{2}{f}_y^{k+1}\\ g_x^{k+1}& g_y^{k+1}\end{array}\right]}^i}}}{\displaystyle \underset{\mathrm{\Delta }{X}^i}{\underbrace{{\left[\begin{array}{l}\mathrm{\Delta }{x}^k\\ \mathrm{\Delta }{y}^k\end{array}\right]}^i}}}=-{\displaystyle \underset{F{\left(\cdot \right)}^i}{\underbrace{{\left[\begin{array}{l}{F}_1\left(\cdot \right)\\ F_2\left(\cdot \right)\end{array}\right]}^i}}}.$$ E12

Once the states have been computed for a new time step, the TS can be calculated. To do this, the linear time-varying systems (10) and (11) rearranged as indicated in (13) should be evaluated with the recently computed states and solved for the TS ${\left[\begin{array}{ll}{x}_{\beta }^{k+1}& y_{\beta }^{k+1}\end{array}\right]}^T$ :

$${\displaystyle \underset{J}{\underbrace{\left[\begin{array}{ll}I-\frac{\mathit{\Delta t}}{2}{f}_x^{k+1}& -\frac{\mathit{\Delta t}}{2}{f}_y^{k+1}\\ g_x^{k+1}& g_y^{k+1}\end{array}\right]}}}{\displaystyle \underset{S}{\underbrace{\left[\begin{array}{l}{x}_{\beta }^{k+1}\\ y_{\beta }^{k+1}\end{array}\right]}}}={\displaystyle \underset{B}{\underbrace{\left[\begin{array}{l}{x}_{\beta }^k+\frac{\mathit{\Delta t}}{2}\left(f_x^k{x}_{\beta }^k+f_y^k{y}_{\beta }^k+f_{\beta }^k+f_{\beta }^{k+1}\right)\\ -g_{\beta }^{k+1}\end{array}\right]}}}.$$ E13

It is very important to note that at each time step $\mathrm{\Delta }t$ , the solution of (13) for the TS calculation $x_{\beta }^{k+1}$ and $y_{\beta }^{k+1}$ uniquely requires a forward/backward substitution, which represents a very small computational burden [17]. This is because the coefficient matrix on the left side of (13) corresponds to the Jacobian matrix already factored used in the final NR iteration to solve (12).

3.3. Multiparameter sensitivity

It is clear that the solution of the linear time-varying system (13) uses the same Jacobian matrix already factored for solving the TS calculation with respect to any $\beta$ parameter. Taking advantage of this observation, the solution approach described in Section 3.2 can be directly extended to provide an efficient linear computation of multiparameter trajectory sensitivities, as shown in (14). In this case, it is only necessary to carry out $\mathit{Np}$ forward/backward substitutions to calculate all TS vectors with respect to $\mathit{Np}$ parameters of the system, i.e., $x_{\mathit{\beta i}}^{k+1}$ and $y_{\mathit{\beta i}}^{k+1}$ $\forall i=1,\dots ,\mathit{Np}$ , where $S$ is a $\left(n+m\right)\times N_P$ sensitivity matrix whose $\mathit{Np}$ columns are the $\left(n+m\right)$ -dimension TS vectors with respect to $\mathit{Np}$ parameters:

$${\displaystyle \underset{J}{\underbrace{\left[\begin{array}{ll}I-\frac{\mathit{\Delta t}}{2}{f}_x^{k+1}& -\frac{\mathit{\Delta t}}{2}{f}_y^{k+1}\\ g_x^{k+1}& g_y^{k+1}\end{array}\right]}}}\underset{S}{{\displaystyle \underbrace{\left[\begin{array}{llll}{S}_{\beta 1}& S_{\beta 2}& \cdots & S_{\mathit{\beta Np}}\end{array}\right]}}}={\displaystyle \underset{B}{\underbrace{\left[\begin{array}{llll}{B}_{\beta 1}& B_{\beta 2}& \cdots & B_{\mathit{\beta Np}}\end{array}\right]}}}.$$ E14

3.4. Sensitivity quantification

As numerically demonstrated in [19, 20, 21, 22, 23], when the system approaches an unstable operation condition, the sensitivities of state trajectories have more rapid changes in magnitudes and larger excursions than the state trajectories. Hence, as the system approaches its stability boundary, the TS approach infinity [19]. In this sense, it is possible to associate the sensitivity information with the stability level of the system for a particular system parameter. For this purpose, an Euclidian norm of the trajectory sensitivity vector, referred to as a sensitivity norm, is proposed in [23] as a measure of proximity to instability. This sensitivity norm also permits the computation of the critical parameters whose variations steer the system much faster to an unstable operation condition. In this chapter, such a sensitivity norm is used to provide a time-varying index of proximity to oscillatory instability for a $n_g$ -generator system defined at each integration time step $\mathrm{\Delta }t$ , as shown in (15), which by computing sensitivities of rotor angle and speed trajectories with respect to active power loads measures the load power’s effect on the system’s small-signal stability:

$${\mathit{SN}}_{\rho }\left(t_k\right)=\sqrt{\sum _{m=1}^{n_g}\left({\left(\frac{\partial {\delta }_m\left(t_k\right)}{\partial {\beta }_{\rho }}-\frac{\partial {\delta }_j\left(t_k\right)}{\partial {\beta }_{\rho }}\right)}^2+{\left(\frac{\partial {\omega }_m\left(t_k\right)}{\partial {\beta }_{\rho }}\right)}^2\right)}\forall \rho =1,\cdots ,\mathit{Np},$$ E15

where $j$ denotes the reference generator.

In this case, the time evolution of sensitivities is necessary to quantify the loads’ influence on the possible occurrence of an oscillatory instability, where the highest values of the sensitivity norms ${\mathit{SN}}_{\rho }$ indicate the most sensitive loads for the EP’s stability. Thus, it should be noted that the critical load powers are those with the largest values of sensitivity index within the integration period, not the largest final value.

Note that this sensitivity norm has been successfully applied before for developing suitable approaches to improve the transient stability of power systems, e.g., for the estimation of the critical clearing time of a faulted system [23], the best possible location of FACTS controllers for transient stability enhancement [24, 25, 26], and the thyristor-controlled series compensator (TCSC) control design to enhance transient stability [27].

3.5. Sensitivity initial conditions

The analysis of a stationary operating point based on this approach is performed by keeping the corresponding EP constant during the whole simulation. Thus, during this simulation period with no perturbation, which is considered from the instant of time $t_0$ and the infinity $t_{\infty }$ , the TS also remain into the same stationary behavior established at nonzero constant values $x_{\beta }\left(t_{\infty }\right)\ne 0$ and $y_{\beta }\left(t_{\infty }\right)\ne 0$ . However, as our approach is based on catching the oscillatory behavior of trajectory sensitivities, the values $x_{\beta }\left(t_0\right)=y_{\beta }\left(t_0\right)=0$ are arbitrarily used as the initial conditions of the sensitivity variables. Then, such an initial condition perturbation starts an oscillatory behavior on the TS, whose transient period evolves from the perturbed initial condition $x_{\beta }\left(t_0\right)=y_{\beta }\left(t_0\right)=0$ to the stationary final condition $\left(x_{\beta }\left(t_{\infty }\right),y_{\beta }\left(t_{\infty }\right)\right)$ . Such a transient behavior provides the necessary TS information to assess the influence of the states and parameters (load powers) in the stability of the analyzed EP. It is important to note that because of the EP stay undisturbed during the TS transient simulation, $f_x$ , $f_y$ , $f_{\beta }$ and $g_x$ , $g_y$ , $g_{\beta }$ are considered time-invariant matrices in the linear sensitivity model (13), where only the TS ${\left[\begin{array}{ll}{x}_{\beta }^{k+1}& y_{\beta }^{k+1}\end{array}\right]}^T$ are time-varying. Thus, a reduced computational burden is required for the TS simulation.

3.6. Approach for small-signal stability with trajectory sensitivities

The application of the SMA-TS approach to assess the effect of a set of system’s parameters on the stability of the equilibrium points is summarized as follows:

Step 1. For an arbitrary set of fixed parameters $\beta$ , the system’s equilibrium is computed by solving the set of nonlinear algebraic equations (1) for $x$ and $y$ considering $\dot{x}=0$ . The NR algorithm is applied to obtain this solution given by the values $x_e$ and $y_e$ that satisfy $0=f\left(x_e,y_e,\beta \right)$ and $0=g\left(x_e,y_e,\beta \right)$ .

Step 2. Based on the Schur and the implicit function theorems [28, 29], compute the reduced Jacobian matrix $J_R={\left.\left(f_x-f_y{\mathrm{g}}_y^{-1}{\mathrm{g}}_x\right)\right|}_{\left(x_e,y_e,\beta \right)}$ that has the same dynamic and algebraic properties of the system’s Jacobian matrix.

Step 3. Determine the critical eigenvalues of $J_R$ , and perform a SMA to identify the associate state variables.

Step 4. Compute the sensitivities of associated state variables with respect to the selected system’s parameters (load powers) at the equilibrium point, $x_{\beta }^{t\to \infty }$ and $y_{\beta }^{t\to \infty }$ , by solving (16). The integration process is started with initial conditions $x_{\beta }\left(t_0\right)=y_{\beta }\left(t_0\right)=0$ for the parameter sensitivities, while the state and algebraic variables are set at their equilibrium values during the solution process. The time evolution of sensitivities is computed under the assumption that a very small perturbation is carried out in the system such that the state and algebraic variables are infinitesimally perturbed; their values, therefore, can be considered constant during the computation of the sensitivity index. These assumptions permit us to directly relate our proposal to the theory of trajectory sensitivities:

$${\displaystyle \underset{J}{\underbrace{{\left.\left[\begin{array}{ll}I-\frac{\mathit{\Delta t}}{2}{f}_x& -\frac{\mathit{\Delta t}}{2}{f}_y\\ g_x& g_y\end{array}\right]\right|}_{\left(x_e,y_e,\beta \right)}}}}{\displaystyle \underset{S}{\underbrace{\left[\begin{array}{l}{x}_{\beta }^{k+1}\\ y_{\beta }^{k+1}\end{array}\right]}}}={\displaystyle \underset{B}{\underbrace{\left[\begin{array}{l}{x}_{\beta }^k+\frac{\mathit{\Delta t}}{2}\left(f_x{x}_{\beta }^k+f_y{y}_{\beta }^k+2f_{\beta }\right)\\ -g_{\beta }\end{array}\right]}}}$$ E16

Step 5. Quantify the interaction between the system parameters and the associated state variables by using the sensitivity index (15). Since this index is a function of the sensitivities of those state variables directly associated with the oscillatory modes and the critical eigenvalues, it can be used to quantify the effect of the $i$ th parameter on these variables. In this case, the highest values of the sensitivity norms indicate the most sensitive parameters. Furthermore, the sensitivity index value increases as the system is approaching an oscillatory stability problem.

4. Small-signal stability analysis

This section presents the analysis of the TS applied on assessing the SSS analysis. The effectiveness of the SMA-TS approach is numerically tested by analyzing the WSCC 9-buses and 3-generator systems [7] and a reduced equivalent system corresponding to the Mexican power system consisting of 190-buses with 46 generators. For the purpose of the studies presented in this section, the system generators are represented by means of the two-axis model with a simple fast exciter loop containing max/min ceiling limits. For this case the system loads are represented by means of the constant power load model; however, the generality of the method allows to consider any load model.

4.1. Modal analysis WSCC system

In this subsection the conventional modal analysis is employed in order to investigate the small-signal stability of the WSCC power system, whose diagram is given in Figure 1. Table 1 presents the modal analysis for different EPs as the active power at bus 5 is increasing.

Figure 1.

WSCC 9-buses, 3-generators.

$P_{L5}\left(\mathrm{pu}\right)$	${\mathit{\lambda }}_{\mathit{\text{crit}}}$	Associated states	PF
4.2	−0.5085 ± 7.3001i	${\omega }_2,{\delta }_2,E_{d2}^{\prime },E_{q1}^{\prime },E_{\mathit{fd}1}$	1.0, 0.99, 0.18, 0.18, 0.16
4.3	−0.5395 ± 6.8512i	${\omega }_2,{\delta }_2,E_{d2}^{\prime },E_{q1}^{\prime },E_{\mathit{fd}1}$	1.0, 0.99, 0.28, 0.39, 0.35
4.4	−0.0305 ± 6.1462i	${\delta }_2,{\omega }_2,E_{d2}^{\prime },E_{q1}^{\prime },E_{\mathit{fd}1}$	1.0, 0.99, 0.46, 0.99, 0.83
4.5	0.7064 ± 5.8935i	$E_{q1}^{\prime },{\delta }_2,{\omega }_2,E_{\mathit{fd}1},E_{d2}^{\prime }$	1.0, 0.85, 0.85, 0.79, 0.38
4.6	1.5118 ± 5.7190i	$E_{q1}^{\prime },{\delta }_2,{\omega }_2,E_{\mathit{fd}1},E_{d2}^{\prime }$	1.0, 0.81, 0.80, 0.76, 0.32
4.7	2.6677 ± 5.5020i	$E_{q1}^{\prime },{\delta }_2,{\omega }_2,E_{\mathit{fd}1},E_{d2}^{\prime }$	1.0, 0.79, 0.78, 0.74, 0.28
4.8	5.3798 ± 4.5835i	$E_{q1}^{\prime },{\delta }_2,{\omega }_2,E_{\mathit{fd}1},E_{d2}^{\prime }$	1.0, 0.77, 0.75, 0.76, 0.26

Table 1.

Selective modal analysis.

SSS and SMA were performed to investigate the different operating points corresponding to the different levels of loading. Eigen-analysis revealed the SSS of EPs, whereas participation factors were used to identify the most associated states to the critical eigenvalue at each EP, as given in Table 1. The first column represents the load changes at bus 5. The second column presents the critical eigenvalue for the EP. The third column presents the most associated states to the critical mode (eigenvalue), obtained by selecting the highest magnitudes of participation factors, which are given in the fourth column in the table.

As the load embedded at bus 5 increased, the stability of the new EP decreased with respect to the previous one. The power system oscillatory instability called the HB was detected when the load changed from 4.4 to 4.5 p.u. The SMA around the HB revealed generator 2 as the most participative in the unstable EP. It is important to observe in Table 1 that generator 2 is the most associated with the critical eigenvalues for all analyzed EPs.

5. Trajectory sensitivity analysis: WSCC system

In order to test the proposed method based on TS to assess the EP stability, operating points before and after the HB point were investigated [30]. As the active power increases at bus 5, as was performed in [7], the system proximity to the bifurcation points is assessed by using the analysis of TS. The rotor angle and speed sensitivities with respect to the load active powers were traced and observed through the time. The TS oscillations provide qualitative information used to investigate the proximity to bifurcation points. Such sensitivity oscillations agree with the critical eigenvalues obtained by the SMA at each EP, as reported in Table 1. Thus, the oscillatory behavior in the TS indicates that the EP is around the HB point and the critical eigenvalue of the EP is complex. The proximity to the HB point is qualitatively assessed by observing the damping of the TS oscillation, i.e., if the TS oscillation is positively damped, the system is operating before the HB (stable EP); however, if the TS oscillation is undamped, the system operates after the HB (unstable EP).

Multiparameter sensitivity is used to assess the load influence around the HB. Such an assessment requires only one simulation to identify the critical loading direction on approaching HB. The TS computation in all cases started from second one onward.

5.1. Stability around the Hopf bifurcation

In order to qualitatively analyze the oscillatory behavior of the EPs around the HB, Figures 2–4 show the TS of the dynamic variables (generator states) with respect to the load embedded at bus 5. The oscillation waveforms and their peak values allow to assess and to determine the EP stability as well as its most associated states because of the loading increase. Figure 2 shows the TS with respect to $P_5=4.3$ pu where the sensitivity oscillations are damped. This agrees with the corresponding critical eigenvalue ${\lambda }_{\mathit{\text{crit}}}=-0.5395\pm 6.8512i$ which indicates that the system is not too close to the HB point. However, Figure 3 shows sensitivity oscillations with a very small damping when $P_5=4.4$ pu, where the critical eigenvalue ${\lambda }_{\mathit{\text{crit}}}=-0.0305\pm 6.1462i$ is very close to the imaginary axis in the complex plane and hence close to a HB point. This means that a very small variation in the load parameter could steer the system to operate in an unstable EP. Figure 4 shows the TS behavior when $P_5=4.41$ pu, which corresponds to an unstable EP after a HB point. From Figures 2–4, it is clear that the highest peaks of the TS oscillations in all cases correspond to generator 2. Therefore, such a generator is the most influential in the EP stability according to the most associated states reported in Table 1. However, this correspondence between associated states and TS is not always kept as will be shown in the upcoming sections.

Figure 2.

Parameter sensitivities with respect to $P_{L5}=4.3$ pu, ${\lambda }_{\mathit{\text{crit}}}=-0.5395\pm 6.8512i,$ [30].

Figure 3.

Parameter sensitivities with respect to $P_{L5}=4.4$ pu, ${\lambda }_{\mathit{\text{crit}}}=-0.0305\pm 6.1462i,$ [30].

Figure 4.

Parameter sensitivities with respect to $P_{L5}=4.41$ pu, ${\lambda }_{\mathit{\text{crit}}}=0.0462\pm 6.1105i,$ [30].

It should be noted from the figures that the load demand at bus 5 increases as the peaks of TS are higher. This qualitative information indicates that such a load variation has a major impact in the load angle and speed of generator 2 (red lines) than the state variables associated with generator 3 (blue lines). Note that this line of reasoning also applies for unstable EPs, i.e., after passing the HB point, as can be observed in Figure 4. Thus, a short simulation in time is enough for this approach to determine the EP stability and their most influencing generators.

5.2. Simulation efficiency

In order to observe the effect of the time-step integration $\mathrm{\Delta }t$ in the computational burden, Figure 5 shows the rotor angle sensitivity of generator 2 with respect to the active power embedded at bus 5. Such a trajectory was computed using three different time-step values. The figure shows sensitivity oscillations considering a period of 14 seconds.

Figure 5.

Effect of the integration time step on the evolution of parameter sensitivities, [30].

It is important to observe that the big difference between $\mathrm{\Delta }t=0.001$ sec and $\mathrm{\Delta }t=0.01$ sec results in negligible differences in the resulting trajectory sensitivities. For the two cases, the TS calculation required to compute 14000 and 1400 forward/backward substitutions, respectively. Such a difference is equivalent to reduce in 90% the number of sensitivity solutions. For the case with $\mathrm{\Delta }t=0.1$ sec, 140 forward/backward substitutions were required, which represents a reduction of 99% of such solutions required to assess the transient. Although for this case, a small difference between trajectory sensitivities results, this is still negligible. This TS-based method can assess at the same time the effect of $N_P$ parameters in the EP stability, whereas in the method of eigenvalues and modal analysis, it is not possible.

5.3. Most sensitive loads to Hopf bifurcation

The participation factors provide the change of critical eigenvalues to the change of the states of critical machines $\left(\partial \lambda /\partial x\right)$ , which establishes the SMA. However, this analysis does not provide information about the critical parameter (critical load in this case) influencing the EP stability but only the critical generator and its most participative states. Thus, the modal analysis does not allow the direct identification of the most sensitive loading directions in the stability of the EPs around a HB point. In practice, load increments not have a unique loading direction as in the previous study, which results in a valid consideration only for academic interest. All loads are then constantly varying in $N_P$ directions, so that it is very important to have a general tool to assess the stability of the EPs, especially operating under stressed conditions of loading. In this context, besides identifying the critical generators (states), the SMA-TS approach based on TS identifies the loading directions which are most sensitive to oscillatory instabilities.

Multiparameter analysis of TS allows computing trajectory sensitivities with respect to $N_P$ parameters in a power system [17] at the same time, as explained in Section 3.3, and can be used to find out the influence of the different loading directions on the SSS around a HB point. For this purpose, Figure 6 shows the sensitivity norm $\mathit{SN}$ through the time with respect to the three embedded loads in the system, where the highest peaks indicate the maximum influence of the corresponding load in the oscillatory behavior. The load demand corresponds to the base case as provided in [7].

Figure 6.

Loads’ effect on the EP’s stability (WSCC system), [30].

The oscillation of ${\mathit{SN}}_{\rho }$ shows that the load embedded at bus 8 $\left(P_{L8}\right)$ is the most sensitive for the EP. The load at bus 6 $\left(P_{L6}\right)$ is the next most sensitive and finally the load $P_{L5}$ . In order to validate the information provided by the ${\mathit{SN}}_{\rho }$ in Figure 6, one parametric study at a time was carried out for the active powers $P_{L8}$ and $P_{L6}$ as performed in Table 1 for the load $P_{L5}$ . The results are reported in Table 2 as follows: the first column indicates the load nodes, and columns 2 and 3 present the measured powers in the base case and the active power increment from the base case to the HB point to each loading direction, respectively. The fourth column provides the active power magnitude at which a HB point occurs, and lastly, column 5 presents the obtained values of ${\mathit{SN}}_{\rho }$ . It is important to observe in Table 2 that the smallest load increment matches the highest SN value and vice versa, i.e., the highest SN indicates a major change in the EP stability with respect to the corresponding load variation, which determines a shorter way to the HB point. Thus, the smallest increment $\mathrm{\Delta }{P}_{L8}=299$ MW proves that the highest SN value indicates this load takes the system to the HB faster than $\mathrm{\Delta }{P}_{L6}=314$ MW and this in turn faster than $\mathrm{\Delta }{P}_{L5}=316$ MW. This agrees with the ${\mathit{SN}}_{\rho }$ reported in Figure 6.

Node	Pbase (MW)	ΔPHB (MW)	PHB (MW)	SN
8	100.0	299	399	1.723
6	90.0	314	404	0.997
5	125.0	316	441	0.739

Table 2.

Sensitivity norm and Hopf bifurcation (WSCC system).

6. Trajectory sensitivity analysis: the Mexican system

In this section, the study consisted of computing the TS norm for 91 loads embedded in a reduced equivalent of the Mexican energy system, which consists of 190 nodes and 46 generators. The transmission components are divided into 180 transmission lines and 83 power transformers. Lastly, the system contains 26 capacitive compensators in shunt connection. The unifilar diagram of the power system is shown in Figure 7.

Figure 7.

The Mexican power system.

In order to assess the effect of the system loads on the system’s dynamic performance, the sensitivity norms with respect to 91 loads were computed. Figure 8 shows the effect of these sensitivities on a critically stable EP. Note that the active power demanded by loads connected at buses from 150 to 152 is the most sensitive in the EP stability. Therefore, according to the reasoning used into the previous section, the loading increase in such directions will steer in a faster way the system to a HB than the rest of the system loads. It must be pointed out that the computation of the sensitivity norms for the 91 system loads were carried out by using one sole time-domain simulation, which corresponds to solving 91 sensitivity DAE systems, with each one consisting of 702 equations and variables. Thus, the assessment of the 91 loads is equivalent to solving 63882 equations in the same simulation at the same time. However, considering the linearity of the sensitivity systems, the same time-invariant Jacobian matrix is used during the whole time-domain simulation, which considerably reduces the computational burden.

Figure 8.

Loads’ effect on the equilibrium point stability (Mexican system), [30].

Once the critical loads have been identified from Figure 8, it is possible to know the most affected generators by the most sensitive loads. Figure 9 shows the TS with respect to the active power at bus 152, which resulted as the most sensitive in the sensitivity norm assessment. The damped oscillations in the TS indicate that the EP is stable and the operation point is not at a HB, which agrees with the corresponding critical eigenvalue of the EP $\lambda =-0.0501\pm 7.8518i$ . It must be observed from Figure 9 that the highest rotor angle sensitivities $\partial {\delta }_{32}/\partial P_{L152}$ and $\partial {\delta }_{33}/\partial P_{L152}$ have identified generators 32 and 33 as the most influenced by the active power embedded at bus 152.

Figure 9.

Evolution of parameter sensitivities with respect to $P_{L152}$ , ${\lambda }_{\mathit{\text{crit}}}=-0.05\pm 7.85i$ , [30].

In order to validate the load ranking influence via the sensitivity norm, Table 3 shows how the increments in the most sensitive loading directions influence the SSS, as well as the proximity to the HB point. Column 1 (Node) indicates the most sensitive loads resulted from the TS analysis shown in Figure 8. Column 2 $\left({\lambda }_{\mathit{\text{crit}}}\right)$ provides the critical eigenvalue for the new EP resulting from such an increment. Columns 3 and 4 show the measured value of active power in the analyzed base case $\left(P_{\mathit{\text{base}}}\right)$ corresponding to Figure 9 and the increment in the specified loading direction $\left(\mathrm{\Delta }P\left(60\mathrm{MW}\right)\right)$ , respectively. Lastly, in columns 5 and 6, $\left(\mathrm{\Delta }{P}_{\mathit{HB}}\right)$ and $\left(P_{\mathit{HB}}\right)$ indicate the increased amount and the value of the active power where the system crosses a HB point by following the corresponding loading directions.

Node	λcrit	Pbase (MW)	ΔP60 (60 MW)	ΔPHB (MW)	PHB (MW)
152	0.021 ± 5.02i	172.64	232.64	53.0	225.64
150	0.018 ± 5.03i	188.24	248.64	54.0	242.24
151	0.015 ± 5.03i	18.72	78.72	54.0	72.72
147	0.012 ± 5.05i	104.00	164.00	56.0	160.00
153	−0.011 ± 8.83i	78.00	138.00	63.0	141.00
145	−0.015 ± 5.13i	83.20	143.20	72.0	155.20
120	−0.050 ± 7.85i	308.88	368.88	189.0	497.88

Table 3.

Loads’ sensitivity norm to Hopf bifurcation in the Mexican system.

It is important to outline that the load effect in the EP stability is not only dependent on the magnitude but also on the topologic location of loads. For example, the power demand embedded at bus 120 is 17 times larger than the load at bus 151; however, the load at 151 resulted in being more sensitive than the load embedded at bus 120, as shown in Table 3, column 3. It must be observed that the most sensitive loads (loads 152–147) provided a major change in the critical eigenvalue and thus in the SSS. The same increment in the most sensitive loading directions (buses 152–147) led the system to oscillatory instability due to a HB point, whereas with the increment in the least sensitive loading directions, the system remained stable. Then, the stability margins in the most sensitive loading directions become more reduced; therefore, according to the sensitivity ranking in Table 3, as the most sensitive loads were increased, the appearance of the HB was found faster as can be observed in column 5. Once more the SMA-TS approach has been successfully proved by determining that the most sensitive loads indicate the shorter ways toward the small-signal instability of the electric power systems.

7. Conclusions

In this chapter an alternative approach for monitoring the Hopf bifurcations along variations in multidimensional loading directions by using a time-domain method is presented, which is based on trajectory sensitivities. This approach, SMA-TS, is general and flexible, i.e., the size of the power systems, as well as the complexity of their mathematical modeling, does not represent any restriction. SMA-TS allows to identify the critical loading directions that steer the system to Hopf bifurcation points. Such an approach was tested in the 9-buses, 3-generators system as well as in 190-buses, 46-generators system. Regardless of the number of sensitivity parameters and system dimensions, SMA-TS requires only one simulation. Such a method keeps constant the Jacobian matrix of the system, requiring only one evaluation and factorization during the whole simulation. The computational effort then consists of performing just one forward/backward substitution at each time step. Furthermore, the approach can handle a very large integration step to drastically reduce the computational effort. Lastly, its application is suitable for real-time monitoring and security assessment in energy management systems.

Acknowledgments

The authors thank the University of Guanajuato for providing the financial support with PFCE 2018 resource for this publication.

Conflict of interest

All the authors of this chapter declare to have no any conflict of interest related with any person, company, institution, etc.