A Highly Relevant Method for Incorporation of Shunt Connected FACTS Device into Multi-Machine Power System to Dampen Electromechanical Oscillations

Abstract

A number of techniques have been proposed to dampen the power system oscillations in the electric power systems. Flexible alternating current transmission system (FACTS) devices are becoming one of them. Among the FACTS family, the static synchronous compensator (STATCOM), a shunt connected FACTS device, has been widely used to provide smooth and rapid steady state, limit transient voltage, and improve the power system stability and performance by absorbing or injecting reactive power. However, the influence ability depends on its placement, control signal, and place of receiving-signal in the network. In order to satisfy these issues, this paper proposes a method for optimal setting and signal position of the STATCOM into the multi-machine power systems with the aim for damping the electromechanical oscillations. This method is developed from the energy approach based on Gramian matrices considering multiple tasks on the Lyapunov equation, in which the observability Gramian matrix is used to seek an optimal location for STATCOM placement. The another is the controllability one used to determine the best local input signal placement that is chosen as a feedback signal for the power oscillation damping (POD) of STATCOM. In addition, the Krylov-based model reduction method is introduced to shorten the calculation time. The proposed method has been verified on the IEEE 24-bus system by analyzing the small-signal stability to search several feasible placements, and then the transient stability is analyzed to compare and determine an optimal placement through testing various cases. The obtained result is also compared with other optimal method.

1. Introduction

In recent years, the demand of electrical power increases due to the social, political, and technological perspective. The distributed generators (DGs) and renewable energy generation have become an inevitable trend because of the economic, environmental, and critical factors and limiting available primary energy resources. When they are connected into the existing power system far away from the load centers, leading on an increase in the scale and complexity of the electric power systems, most power systems operate with requirement very close to their limits and always constrained by small-signal and transient stability conditions. If the inertial response and power damping capability are an insufficient quality when disturbances in the power system occur, low frequency oscillations may appear. These oscillations are due to the outcome of dynamical coactions between the generators within a system. Thus, it can lead to the generators outage, lines tripping, and regions blackout of the system.

In the past five decades, the power systems around the world have faced more blackouts with various severities due to human and technical faults. Power-technology.com has listed out the world’s top 10 worst blackouts, such as Indian on 30 and 31 July 2012; Paraguay-Brazil on 10 November 2009; Chenzhou in China in January and February 2008; Germany, France, Italy and Spain on November 2006; Indonesia on 18 August 2005; Greece on 12 July 2004; Switzerland on 28 September 2003; Eight U.S. states and Canada on 14 August 2003; Quebec and parts of the U.S. on March 1989; New York City on 13 July 1977; and USA on 9 November 1965, all of which have happened in the last 50 years. Among them, the Indian grid disturbance is a notable event, leaving 620 and 370 million people without electric power for several hours on two days, 30 and 31 July 2012, respectively. The least event is the power system blackout in South Australia on 28 September 2016. Figure 1 is an example of the instable inter-area event of recording a wide area measuring system (WAMS) on 9 February 2006 in Spain [1]. Clearly, observing from this figure, it can be shown significant changes of the dynamic behavior of the system due to the loss of generated power 1.2 GW. The frequency decreases immediately in the proximity of one of power plant outage. The decrease of frequency spreads over the whole system and finally reaches the Romania system with a time delay of about 2.0 s. This event tells us that the inter-area oscillations have a frequency range 0.21–0.25 Hz.

The low frequency electromechanical oscillations happen in the power systems due to the contingencies. These oscillations have a frequency range 0.1–2.0 Hz, classified into two groups based on the energy of oscillation: (i) occurring between generators in a small group, called local mode, having a frequency range 0.7–2.0 Hz; and (ii) occurring between generators in large groups, called inter-area oscillations, having a frequency range 0.1–0.8 Hz [2,3]. Such oscillations menace the stability point of the power system, thus damping these oscillations is necessary for the sure system operation. From the viewpoint of the system instability due to the limited value increment for current network structure, accordingly, most of the efforts for damping oscillations focus on using the controllers to improve the oscillation damping ability, such as Flexible alternating current transmission system (FACTS) devices, power system stabilizers (PSSs), energy storage systems (ESSs), and so on. The installation of PSS in the system is the traditional method. For this method, the adverse effects are the local-area mode oscillations damping, the variation in the voltage profile under serious disturbances, the leading power factor [4]. FACTS is widely used to enhance the operational flexible of the system under other operation conditions, and especially to damp out local- and inter-area mode oscillations [5,6]. This device has the ability to control the parameters of the power system or to modify the load flow, so that the power system may be improved stability and damped oscillations when the grid faults occurred. Thus, FACTS could be the best alternative to enhance the stability margin of the existing power systems and it provides against the inter-area oscillation damping better than PSSs [7]. For ESS device, from the point of view of the power system control, the ESSs are better than the FACTS and PSSs, because they cannot only regulate the reactive power exchange but also the active power with the power system where it is installed. The purpose of applying ESS technology in power is the same FACTS and PSS. However, in order to enhance the power system oscillation damping, ESSs are more flexible than PSSs [8] because of exchanging active power directly. For instance, Du et al. [9] have explained why and how the battery energy storage system (BESS) based on stabilizer can dampen the inter-area oscillations in a multi-machine power system. Chen et al. [10] introduced the complex torque coefficient method to analyze the use of the energy storage system ESS based on flywheel storage technology for damping the power system oscillations of the infinite-bus power system considering single-machine. For this paper, we only focus to dampen the power system oscillations using FACTS device. However, the benefits of FACTS depend greatly on their placement in network [2,7,11]. Consequently, the optimal location of FACTS devices, and especially static synchronous compensator (STATCOM) in the multi-machine, is an interesting research topic.

In the last decades, many researchers have proposed techniques for locating the location of FACT based on two groups: (i) analytical methods; and (ii) heuristic optimization approaches [12,13,14]. Methods based on heuristic algorithms include the genetic algorithm (GA) [15], the particle swarm optimization (PSO) technique [16], the cuckoo search algorithm [17], modified artificial bee colony algorithm (MABCA) [18], and simulate annealing (SA) [19]. These methods were used to increase system load-ability and the system security margin. Evolutionary algorithm (EA) [20] and bacterial swarming algorithm (BSA) [21] are used to minimize the real power loss in the transmission lines and to improve voltage profile at load buses. The bees algorithm (BA) is applied by [22] to maximize the available transfer capability of power agreement between sink areas and source. Contrasting to the heuristic methods, some authors proposed the analytical methods. For example, Sadikovic et al., Magaji et al. [23,24] have proposed the residue factor method to damp inter-area mode of oscillations. Kunnar et al. [25] have utilized the modal controllability index to dampen inter-area oscillations. Van Dai et al., Le et al. [2,11], have introduced the energy method to dampen the power system oscillations in IEEE 39-bus power system, and Zamora-Cárdenas et al. [26] proposed the multi-parameter trajectory sensitivity approach to determine the optimal placement of series FACTS devices with the objective of providing enough level of transient stability.

In addition, the optimal location of STATCOM highly depends on the supplemental controller, namely the power oscillation damping (POD) controller and its feedback signal position [27]. Leading in recent years, many of the optimization methods have been proposed to find out ways for answering the question of which POD controller and its feedback signal could result in the STATCOM having the considerable effect on the system. For example, genetic algorithm (GA) is proposed by Eshtehardiha et al., Panda et al. [28,29], but this method requires a very long run time to the large-scale systems. Safari et al. [30] introduced the Artificial Bee Colony (ABC) method, but the convergence is very slow. Abd-Elazim et al. [31] introduced the imperialist competitive algorithm (ICA) for optimal design of STATCOM parameters. In addition, local input signals as the line reactive and active power, line current, and bus voltage are all good selections for the feedback signal of POD controller of FACTS to dampen the power system oscillations [2,11]. The active power in the transmission line is considered as an effective input signal for POD controller [32]. Reference [33] concluded that the current or active power is not difference when using them as the feedback signal. For this paper, the active power in the transmission line is chosen as a feedback signal. However, in the previous proposed researches, the authors did not mention the location of input signal placement.

It can be observed that most employed techniques in the previous literature have several drawbacks: (i) There is no mention of the optimal location of the local input signal for the POD of STATCOM; (ii) Slow convergence in search stage, limit of local search ability, and the requested process time is very long when studying the large-scale systems, so they just focus on analyzing the small-scale power systems; (iii) The calculation of critical modes may be doubtful in the case of the larger-scale and complex power systems because they may not be unique. Furthermore, the calculation of these critical modes also depends on the local- or inter-area mode; (iv) The calculation of participation coefficients is only based on the state variables and neglects the input–output behaviors. Therefore, in order to overcome these drawbacks, this paper proposes a feasible method. This method is a conjunction between the energy method based on the observability and controllability Gramian matrices and the Krylov-based model reduction method to determine the optimal location of STATCOM and the local input signal position of POD controller with aim of improving the stability of the large-scale power systems, respectively.

The main contribution of this paper is to propose a feasible method:

(i)

To determine the optimal placement of STATCOM and the best local input signal position of its POD controller;

(ii)

To limit the time calculation when analyzing the complexity and large-scale power systems on the small-signal stability analysis.

The remainder of this paper is organized as follows: Section 2 analyzes the mathematical model of the STATCOM, power system, and incorporation the STATCOM in the power system. Section 3 explains the methods for designing POD of STATCOM and selecting the location of STATCOM-POD and location of the feedback. The details of the test cases are given in Section 4. Finally, the conclusion is outlined in Section 5 and the Arnoldi algorithm, QR decomposition, and parameters of the STATCOM are detailed in Appendix A.

2. Analysis of Mathematical Models

2.1. STATCOM Model

The STATCOM is one the most noteworthy series of the FACTS derives. The main function of STATCOM is to control a fast response of the transmission voltage, increase the power transfer on transmission line, and improve the transient stability of the power network by generating and/or absorbing reactive power after clearing fault. In this paper, the structure of a STATCOM consisting of a voltage source converter (VSC), a coupling transformer, a DC capacitor, and controller unit is used to study, as shown in Figure 2a. Its operational Voltage–Current characteristic is shown in Figure 2b.

The objective of this paper is to investigate the damping of power system oscillations when the grid fault occurs, STATCOM should operate over the rated maximum inductive and/or capacitive compensation range being independent of the AC system voltage. Accordingly, the controller is designed to control STATCOM, so that it always keeps full capability during most severe contingencies. The output current variable can be obtained by using various control approaches. Besides, this output current variable depends on the position and kind (active power, reactive power, current of transmission line, transmission angle, or generator rotor speed) of the feedback signal. For this purpose, the dynamic model is proposed for STATCOM [34]. This dynamic mode is compound of the voltage regulator with transient gain, and can be determined by the time constants T₁ through T₄ and integrator K. The steady-state gain is equal to the inverse of Droop, which allows the sharing of voltage of generators. The inputs to model are the monitored voltage V_bus, voltage reference on the connected STATCOM bus V_ref, the STATCOM internal voltage V_T, and the auxiliary signal. This auxiliary signal is generated from the power oscillation damping (POD) controller. This POD consists of an amplification block, a wash-out filter, a low-pass filter, and m_c stages of lead-lag blocks as shown in Figure 3, where T_m, T_w, T_lead, and T_lag are the measurement, washout, lead, and lag time constants, respectively. The input signal could be the line active and reactive power, current, bus voltage, generator angle speed, or phase angle.

The output current of STATCOM is varied according to firing angle control, which is adjusted through the voltage regulator with transient gain and integrator gain to main the STATCOM bus voltage at the desired reference value. By linearizing the STATCOM dynamic model with POD controller, we can be obtained as follows [35]:

$$\Delta {\dot{X}}_{\mathrm{Stat}}=A_{\mathrm{Stat}}\Delta X_{\mathrm{Stat}}+B_{\mathrm{Stat}}\Delta V_{\mathrm{Stat}},$$

(1)

where the matrices A_Stat and B_Stat depend on the time constants T₁ through T₄, gain K, the monitored voltage V_bus, the voltage reference V_ref, and the STATCOM internal voltage V_T.

2.2. Power System Model

The methodological method of the dynamic modeling of the multi-machine power system with m-machine and n-bus has been described in [36] applied to this case study. In this model, each synchronous generator is represented by two-axis flux decay dynamic and the Excitation system is used the IEEE type I to provide the generator [37]. The differential algebraic equation (DAE) model of the multimachine power system without FACTS can be expressed the following equation:

$$\{\begin{array}{l}x=f(x,y,u),x(0)=x_0;\\ 0=g(x,y,u),y(0)=y_0;\end{array}$$

(2)

in which x, y, and u are the state, algebraic, and input vectors, respectively and are defined as:

$$\{\begin{array}{l}x={\left[{\mathsf{\delta }}_i,{\mathsf{\omega }}_i,E_{qi}^{\text{'}},E_{di}^{\text{'}},E_{fdi},V_{Ri},R_{Fi}\right]}^{\mathrm{T}},\\ y={\left[V_j,{\mathrm{\theta }}_j,I_{di},I_{qi}\right]}^{\mathrm{T}},\\ u={\left[T_{Mi},V_{REFi}\right]}^{\mathrm{T}},\\ i=1,\mathrm{...},m;j=1,\mathrm{...},n;\end{array}$$

(3)

where

$\mathrm{T}$	the transpose operator;
${\mathsf{\delta }}_i$	the rotor angle of generator ith;
${\mathsf{\omega }}_i$	the speed of generator ith;
$V_j$	the voltage magnitude at bus jth;
${\mathrm{\theta }}_j$	the phase angle at bus jth;
$I_{di}$	the d-axis component of the current of generator ith;
$I_{qi}$	the q-axis component of the current of generator ith;
$V_{Ri}$	the input amplifier voltage the excitation of generator ith;
$R_{Fi}$	the stabilizer feedback variable of the excitation of generator ith;
$T_{Mi}$	the electrical power of generator ith;
$E_{qi}^{\text{'}}$	the q-axis component of the internal voltage of generator ith;
$E_{di}^{\text{'}}$	the d-axis component of the internal voltage of generator ith;
$E_{fdi}$	the d-axis component of the field voltages of generator ith; and
$V_{REFi}$	the reference voltage of generator ith.

According to the methodology presented in [36], the linearized model is considered and given the following symbolic form:

$$\left[\begin{array}{c}\Delta \dot{X}\\ 0\\ 0\end{array}\right]=\left[\begin{array}{ccc}{A}^{\text{'}}& B_1^{\text{'}}& B_2^{\text{'}}\\ C_1^{\text{'}}& D_{11}^{\text{'}}& D_{12}^{\text{'}}\\ C_2^{\text{'}}& D_{21}^{\text{'}}& D_{22}^{\text{'}}\end{array}\right]\left[\begin{array}{c}\Delta X\\ \Delta y_1\\ \Delta y_2\end{array}\right]+\left[\begin{array}{c}{E}^{\text{'}}\\ 0\\ 0\end{array}\right]\Delta U,$$

(4)

Equation (4) can be rewritten as:

$$\{\begin{array}{l}\Delta \dot{x}=A^{\text{'}}\Delta x+B^{\text{'}}\Delta y+E^{\text{'}}\Delta u,\\ 0=C^{\text{'}}\Delta x+D^{\text{'}}\Delta y,\end{array}$$

(5)

where A’, B’, C’, D’, and E’ are appropriate Jacobians of the system in Equation (2) evaluated at the operating point, and can identify as follows:

$$\{\begin{array}{l}{A}^{\text{'}}=\frac{\partial f}{\partial x},B^{\text{'}}=\left[B_1^{\text{'}}{B}_2^{\text{'}}\right]=\frac{\partial f}{\partial y},C^{\text{'}}=\left[\begin{array}{l}{C}_1^{\text{'}}\\ C_2^{\text{'}}\end{array}\right]=\frac{\partial g}{\partial x},\\ D^{\text{'}}=\left[\begin{array}{l}{D}_{11}^{\text{'}}{D}_{12}^{\text{'}}\\ D_{21}^{\text{'}}{D}_{22}^{\text{'}}\end{array}\right]=\frac{\partial g}{\partial y},E^{\text{'}}=\left[\begin{array}{l}{E}_1\\ 0\\ 0\end{array}\right]=\frac{\partial g}{\partial u},\Delta y=\left[\Delta y_1,\Delta y_1\right],\end{array}$$

(6)

and

$$\{\begin{array}{l}\Delta x={\left[\Delta {\mathsf{\delta }}_i,\Delta {\mathsf{\omega }}_i,\Delta E_{qi}^{\text{'}},\Delta E_{di}^{\text{'}},\Delta E_{fdi},\Delta V_{Ri},\Delta R_{Fi}\right]}^{\mathrm{T}},\\ \Delta y={\left[{\mathrm{\theta }}_1{V}_1\mathrm{...}{V}_m|{\mathrm{\theta }}_2\mathrm{...}{\mathrm{\theta }}_n{V}_{m+1}\mathrm{...}{V}_n\right]}^{\mathrm{T}},\\ \Delta u={\left[\Delta T_{Mi},\Delta V_{REFi}\right]}^{\mathrm{T}}.\end{array}$$

(7)

Therefore, this methodology utilized to analyze the small-signal stability programs using PSS/E (Power System Simulator for Engineering) and Matlab softwares.

2.3. The STATCOM-Power System Model

The STATCOM can be installed at the existing load bus or the extra bus that created between two the existing buses of the multimachine power system. As such, in this case study, a STATCOM is installed at the bus t in the n-bus power system, the injected power flow into the bus t are given the following equations [38]:

$$\{\begin{array}{l}{P}_t=P_{\mathrm{Stat}}+{\displaystyle \sum _{k=1}^n(V_t^2{G}_t-V_t{V}_k[G_{tk}\cos ({\mathrm{\theta }}_t-{\mathrm{\theta }}_k)+B_{tk}\sin ({\mathrm{\theta }}_t-{\mathrm{\theta }}_k)]),}\\ Q_t=Q_{\mathrm{Stat}}-{\displaystyle \sum _{k=1}^n(V_t^2{B}_t-V_t{V}_k[G_{tk}\cos ({\mathrm{\theta }}_t-{\mathrm{\theta }}_k)+B_{tk}\sin ({\mathrm{\theta }}_t-{\mathrm{\theta }}_k)],}\end{array}$$

(8)

where P_Stat and Q_Stat are, respectively, the active and reactive power at converter terminal at the bus t, that are:

$$\{\begin{array}{l}{P}_{\mathrm{Stat}}=(V_{\mathrm{Stat}}^2{G}_{\mathrm{Stat}}-V_t{V}_{\mathrm{Stat}}[G_{\mathrm{Stat}}\cos ({\mathrm{\theta }}_{\mathrm{Stat}}-{\mathrm{\theta }}_t)+B_{\mathrm{Stat}}\sin ({\mathrm{\theta }}_{\mathrm{Stat}}-{\mathrm{\theta }}_t)]),\\ Q_{\mathrm{Stat}}=-(V_{\mathrm{Stat}}^2{B}_{\mathrm{Stat}}-V_t{V}_{\mathrm{Stat}}[B_{\mathrm{Stat}}\cos ({\mathrm{\theta }}_{\mathrm{Stat}}-{\mathrm{\theta }}_t)+B_{\mathrm{Stat}}\sin ({\mathrm{\theta }}_{\mathrm{Stat}}-{\mathrm{\theta }}_t)])\end{array}$$

(9)

in which V_t, θ_t, V_k, θ_k, V_Stat and θ_Stat are the voltage magnitudes and phase angles at bus t, k, and converter terminal, respectively; and Y_t = G_t + jB_t, Y_tk = G_tk + jB_tk, and Y_Stat = G_Stat + jB_Stat are the values of self-admittance at bus t for the existing n-bus power system, between buses t and k, and converter terminal, respectively.

Linearizing Equation (8) and combining Equation (1), the linearized model of STATCOM can be derived as the following equation:

$$\{\begin{array}{l}\Delta {\dot{X}}_{\mathrm{Stat}}=A_{\mathrm{Stat}}\Delta X_{\mathrm{Stat}}+B_{\mathrm{Stat}}{\left[\Delta {\mathrm{\theta }}_t,\Delta V_t,\Delta {\mathrm{\theta }}_{\mathrm{Stat}},\Delta V_{\mathrm{Stat}}\right]}^{\mathrm{T}},\\ {\left[\Delta P_t,\Delta Q_t,\Delta P_{\mathrm{Stat}},\Delta Q_{\mathrm{Stat}}\right]}^{\mathrm{T}}=C_{\mathrm{Stat}}\Delta X_{\mathrm{Stat}}+D_{\mathrm{Stat}}{\left[\Delta {\mathrm{\theta }}_t,\Delta V_t,\Delta {\mathrm{\theta }}_{\mathrm{Stat}},\Delta V_{\mathrm{Stat}}\right]}^{\mathrm{T}},\end{array}$$

(10)

where A_Stat, B_Stat, C_Stat, and D_Stat are appropriate Jacobians of the system in Equation (8) evaluated at the operating point; the voltage magnitude V_Stat and phase angle θ_Stat at converter terminal are considered as state variables [38].

The methodology for incorporating FACTS devices into the n-bus power system is developed [35,36,38,39,40]. Based on this methodology, the DAE model of m-machine and n-bus power system with the STATCOM can be obtained by combining Equations (5) and (10):

$$\left[\begin{array}{c}\Delta \dot{X}\\ \Delta {\dot{X}}_{\mathrm{Stat}}\\ 0\\ 0\end{array}\right]=\left[\begin{array}{llll}{A}^{\text{'}}& A_{1\mathrm{Stat}}& B_1^{\text{'}}& B_2^{\text{'}}\\ A_{2\mathrm{Stat}}& A_{\mathrm{Stat}}& B_{1\mathrm{Stat}}& B_{2\mathrm{Stat}}\\ C_1^{\text{'}}& C_{1\mathrm{Stat}}& D_{11}^{\text{'}}& D_{12}^{\text{'}}\\ C_2^{\text{'}}& C_{2\mathrm{Stat}}& D_{1\mathrm{Stat}}^{\text{'}}& D_{2\mathrm{Stat}}^{\text{'}}\end{array}\right]\left[\begin{array}{c}\Delta X\\ \Delta {\dot{X}}_{\mathrm{Stat}}\\ \Delta y_1\\ \Delta y_2\end{array}\right]+\left[\begin{array}{c}{E}^{\text{'}}\\ 0\\ 0\\ 0\end{array}\right]\Delta U.$$

(11)

Equation (11) can be rewritten as:

$$\{\begin{array}{l}\Delta {\dot{x}}_{\mathrm{new}}=A_{\mathrm{new}}\Delta x_{\mathrm{new}}+B_{\mathrm{new}}\Delta y_{\mathrm{new}}+E_{\mathrm{new}}\Delta u_{\mathrm{new}},\\ 0=C_{\mathrm{new}}\Delta x_{\mathrm{new}}+D_{\mathrm{new}}\Delta y_{\mathrm{new}}\end{array}$$

(12)

in which, A_new, B_new, C_new, D_new, and E_new are appropriate Jacobians of the new system evaluated at the operating point, and can identify as follows:

$$\{\begin{array}{l}{A}_{new}=\left[\begin{array}{ll}{A}^{\text{'}}& A_{1\mathrm{Stat}}\\ A_{2\mathrm{Stat}}& A_{\mathrm{Stat}}\end{array}\right],B_{new}=\left[\begin{array}{ll}{B}_1^{\text{'}}& B_2^{\text{'}}\\ B_{1\mathrm{Stat}}^{\text{'}}& B_{1\mathrm{Stat}}^{\text{'}}\end{array}\right],C_{new}=\left[\begin{array}{ll}{C}_1^{\text{'}}& C_{1\mathrm{Stat}}\\ C_2^{\text{'}}& C_{2\mathrm{Stat}}\end{array}\right],\\ D_{new}=\left[\begin{array}{ll}{D}_{11}^{\text{'}}& D_{12}^{\text{'}}\\ D_{1\mathrm{Stat}}^{\text{'}}& D_{2\mathrm{Stat}}^{\text{'}}\end{array}\right],E_{new}={\left[E^{\text{'}},0,0,0\right]}^{\mathrm{T}},\Delta y_{\mathrm{new}}=\left[\Delta y_1,\Delta y_2\right]\end{array}$$

(13)

and

$$\{\begin{array}{l}\Delta x_{\mathrm{new}}={\left[{(\Delta {\mathsf{\delta }}_{\mathrm{new}})}_i,{(\Delta {\mathsf{\omega }}_{\mathrm{new}})}_i,{(\Delta E_{\mathrm{new}}^{\text{'}})}_{qi},{(\Delta E_{\mathrm{new}}^{\text{'}})}_{di},(\Delta E_{\mathrm{new}fdi}),(\Delta V_{\mathrm{new}Ri}),{(\Delta R_{\mathrm{new}})}_{Fi}\right]}^{\mathrm{T}},\\ \Delta y_{\mathrm{new}}={\left[{({\mathrm{\theta }}_{\mathrm{new}})}_1\left(V_{\mathrm{new}}{)}_1\mathrm{...}{(V_{\mathrm{new}})}_m|{({\mathrm{\theta }}_{\mathrm{new}})}_2\mathrm{...}{({\mathrm{\theta }}_{\mathrm{new}})}_n\right(V_{\mathrm{new}}{)}_{m+1}\mathrm{...}{(V_{\mathrm{new}})}_n\right]}^{\mathrm{T}},\\ \Delta u_{\mathrm{new}}={\left[{(\Delta T_{\mathrm{new}})}_{Mi},{(\Delta V_{\mathrm{new}})}_{REFi}\right]}^{\mathrm{T}}.\end{array}$$

(14)

The meaning of all state variables in Equation (14) is the same as that in Equation (3), but corresponding to the power system with STATCOM.

3. Methods Analysis

In this section, it elaborates on mathematical steps taken to formulate expressions for finding optimal STATCOM placement. The method based on two controllability and observability properties of the system, playing an important role in determining the location of placement and input signal of STATCOM, is explained. In order to evaluate these controllability and observability indexes, the linear dynamic system will be considered under the state-space form is necessary to detail conjunction between output and input signals, which contain more information about the internal states of system. Based on Equation (12), a linear time invariant system with FACTS except POD controller described in the state-space form as follows:

$$\{\begin{array}{l}\dot{x}(t)=Ax(t)+Bu(t),\\ y(t)=Cx(t)+Du(t),\end{array}$$

(15)

where

x(t)	the state vector of length equal to the number of states n;
y(t)	the output vector of length equal to the number of output m;
u(t)	the input vector of length equal to the number of input r;
A	the state matrix of size n × n;
B	the control matrix of size n × r;
C	the output matrix of size m × n; and
D	the feed-forward matrix of size m × r.

3.1. The Order Reduction Method

For a large-scale power system, the number of state variables could be big, taking a long time to process even when the modern computer systems are used. The small-signal stability analysis only needs to take into account the number of the important state variables, being normally smaller than the number of the state variables of the original linearized power system. In order to solve this problem, the order reduction solutions were studied in [41], in which the Krylov-based model reduction solution is selected and reusable to study, and can be summed up as below.

The transfer function of the system in Equation (15) can be described as follows:

$$G(s)=C{(sI-A)}^{-1}B+D,$$

(16)

and can be expanded under the Laurent series form around certain point s₀

$$G(s)={\mathrm{\mu }}_0(s_0)+{\mathrm{\mu }}_1(s_1)\frac{(s-s_0)}{1!}+\cdots +{\mathrm{\mu }}_k(s_k)\frac{{(s-s_k)}^k}{k!}+\cdots ,$$

(17)

where μ_{(0, 1, 2, …, k)} denotes the moments and the kth moment can be defined as:

$${{\mathrm{\mu }}_k(s_0)={(-1)}^k\frac{d^k}{d{s}^k}G(s)|}_{s=s_0}\in {\mathcal{R}}^{m\times r},k>0,$$

(18)

it should be noted that the Laurent series expansion can be done around different points.

Expanding the Laurent series in the neighborhood of the arbitrary complex number (s₀ ∈ $\infty$), the moments are:

$${\mathrm{\mu }}_0(s_0)=D+C{(s_{0}I-A)}^{-1}B,{\mathrm{\mu }}_k(s_0)=C{(s_{0}I-A)}^{-(k+1)}B,k>0.$$

(19)

If the moments are selected at zero (s₀ = 0), the resulting problem is known as Pade approximation, these moments are:

$${\mathrm{\mu }}_0(0)=D+C{A}^{-1}B,{\mathrm{\mu }}_k(0)=C{A}^{-(k+1)}B,k>0.$$

(20)

Expanding the Laurent series around infinity (s₀ = $\infty$), the moments are obtained by deriving from Equations (16) and (18) as:

$${\mathrm{\mu }}_0(\infty )=D,{\mathrm{\mu }}_k(\infty )=C{A}^{k-1}B,k>0,$$

(21)

are called Markov parameters, and the resulting problem is known as partial realization.

In order to solve the above-mentioned model reduction, the algorithm is proposed in [42], namely, the Arnoldi algorithm is reusable. Its procedure is shown in Appendices A.1 and A.2.

After applying the Krylov-based model reduction solution for the original system of large dimension (15), the reduced system obtained as:

$$\{\begin{array}{l}\tilde{\dot{x}}(t)=\tilde{A}\tilde{x}(t)+\tilde{B}u(t),\\ y(t)=\tilde{C}\tilde{x}(t)+\tilde{D}u(t),\\ {\mathrm{\mu }}_j={\tilde{\mathrm{\mu }}}_j,j=1,2,\cdots ,k\end{array}\mathrm{in}\mathrm{which}|\begin{array}{l}\tilde{A}=V_t^{\mathrm{T}}A{V}_t,\tilde{B}=V_t^{\mathrm{T}}B,\\ \tilde{C}=C{V}_t,\tilde{D}=D,\end{array}$$

(22)

As can be observed in Equation (22), the moments of the original and reduced systems equal to the kth term of the Laurent series expansion, and the reduction does not change the inputs and outputs of the original system.

3.2. The Conventional Method

For the purpose of comparison and recognition to the proposed method, an approach, namely the residue method is introduced in [6], again employed in this paper to determine the best siting of STATCOM by tuning parameters of the POD controller. The procedure is re-summed up to suit the objective of this study as follow:

The transfer function of the reduced system in Equation (22) can be described

$$\tilde{G}(s)=\frac{y(s)}{u(s)}=\tilde{C}{(sI-\tilde{A})}^{-1}\tilde{B}+\tilde{D}.$$

(23)

The transfer function in Equation (23) between the gth output and hth input can be expressed in terms of modes and residues as:

$${\tilde{G}}_{gh}(s)=\frac{y_g(s)}{u_h(s)}={\displaystyle \sum _{i=1}^N\frac{{\tilde{C}}_g{t}_i{\upsilon }_i{\tilde{B}}_g}{(s-{\lambda }_i)}}={\displaystyle \sum _{i=1}^n\frac{R_{igh}}{(s-{\lambda }_i)}},$$

(24)

where λ_i is the ith eigenvalue; t_i and υ_i denote the right and left eigenvectors associated with the ith eigenvalue λ_i, respectively; N is the total number of eigenvalues; R_igh is the residue associated with ith mode and it gives the measure of the sensitivity of mode to a feedback signal between the output gth and input hth; and this R_igh can also be expressed in terms of the controllability and observability of mode.

Considering a power system shown in Figure 4a, in which the system includes STATCOM unit and its POD controller. The control close-loop system shown in Figure 4b, where POD(s) is the transfer function of STATCOM-POD controller and $\tilde{G}(s)$ is the transfer function of the reduced system. As proposed in this study, POD is a transfer function, as shown in Figure 3, including an amplification block, a wash-out and low-pass filters, and m_c stages of lead-lag block, and can be given by:

$$POD(s)=K_p\underset{PO{D}_1(s)}{\underbrace{(\frac{1}{1+s{T}_m})(\frac{s{T}_w}{1+s{T}_w}){(\frac{1+s{T}_{lead}}{1+s{T}_{lag}})}^{m_c}}},$$

(25)

where K_p is the constant gain. T_m, T_w, T_lead, and T_lag are the measurement, washout, lead, and lag time constants, respectively.

When applying the feedback control, the eigenvalues of the reduced system $\tilde{G}(s)$ are changed, so that the eigenvalue sensitivity ith to feedback gain can be calculated as [6,43]:

$$\frac{\partial {\lambda }_i}{\partial K_p}=R_i\frac{\partial POD({\lambda }_i)}{\partial K_p}=R_{i}PO{D}_1({\lambda }_i),$$

(26)

where R_i denotes the open-loop residue corresponding to eigenvalue ith, λ_i denotes the mode that should be influenced by the STATCOM damping controller. At the initial operating point that is usually the open-loop system, the eigenvalue λ_i can be describes under small changes as follows:

$$\Delta {\lambda }_i=R_{i}POD({\lambda }_i)=R_i\Delta K_{p}PO{D}_1({\lambda }_i),$$

(27)

and this displacement can show in Figure 4c [6]. It can also be seen in this figure that the angle Φ_com is the compensation angle that needs to drive the displacement of eigenvalue. This angle depends on the time constants of the lead-lag function of the POD controller, and these time constants can be calculated by using the following equations:

$$\{\begin{array}{l}{T}_{lead}=\frac{1}{{\mathsf{\omega }}_i\sqrt{{\alpha }_c}},\\ T_{lag}={\alpha }_c{T}_{lead}\end{array}$$

(28)

in which

$${\alpha }_c=\frac{1-\sin \left(\frac{({180}^o-\arg (R_i))}{m_c}\right)}{1+\sin \left(\frac{({180}^o-\arg (R_i))}{m_c}\right)},$$

(29)

where ω_i is the frequency of the ith oscillation mode, arg(R_i) is the residue phase angle associated with the ith oscillation mode, and m_c is the number of compensation stages.

After tuning time constants of the lead-lag, the gain K_p is determined by root-locus method, and it can also calculate this K_p as a function of the desired eigenvalue location ${\lambda }_i^{des}$ based on Equation (27) as follows [6]:

$$K_p=\left|\frac{{\lambda }_i^{des}-{\lambda }_i}{R_{i}PO{D}_1({\lambda }_i)}\right|.$$

(30)

3.3. The Proposed Method

In this section, a method is developed from the Gramian matrix of controllability and observability proposed to determine the optimal location of the sensor (is device using to receive the real power flow from the line and send it to the POD) and STATCOM-POD placements. Two proposed algorithms shown schematically in Figure 5 for these developed applications. Below is the explanation.

3.3.1. Determine Location of Input Signal for POD

The generator speed deviations, bus voltage, real and reactive power flow in the transmission line, and line current are all good nomination in selecting the input signal for the STATCOM control loop [2,11]. The selection of location of these signals placing in the systems with the objective to dampen the power system oscillations is more difficult. In order to solve this problem, a method based on the controllability Gramian matrix of the reduced linearized system is proposed to determine an optimal location for sensor placement. This sensor can receive and send the signal to POD.

For linear time-invariant systems of Equation (22), the observability matrix can be defined as:

$$Co=\left[\tilde{B}\tilde{A}\tilde{B}{\tilde{A}}^2\tilde{B}\cdots {\tilde{A}}^{n-1}\tilde{B}\right]\in {\mathcal{R}}^{(nr)\times n},$$

(31)

and the observability Gramian matrix of the pair ($\tilde{C},\tilde{A}$) at time t is

$$W_{Co}(t)={\displaystyle \underset{0}{\overset{t}{\int }}{\mathrm{e}}^{\tilde{A}\tau }\tilde{B}{\tilde{B}}^{\mathrm{T}}{\mathrm{e}}^{{\tilde{A}}^{\mathrm{T}}\tau }d\tau }\in {\mathcal{R}}^{n\times n},$$

(32)

can be used to calculate the observability of a system. Therefore, there are two possible ways to check the controllability of a system based on these matrices. Firstly, the pair ($\tilde{B},\tilde{A}$) is observable if and only if the controllability matrix Co has full rank n, i.e., rank (Co) = n. Secondly, the pair ($\tilde{B},\tilde{A}$) is controllability if and only if the controllability Gramian matrix W_Co(t) is positive definite for $\forall$t > 0, and it satisfies the following differential Lyapunov matrix equation:

$$\{\begin{array}{l}-\frac{d}{dt}{W}_{Co}(t)+\tilde{A}{W}_{Co}(t)+W_{Co}(t){\tilde{A}}^{\mathrm{T}}=-\tilde{B}{\tilde{B}}^{\mathrm{T}},\\ W_{Co}(0)=0.\end{array}$$

(33)

If the system in Equation (22) is asymptotically stable around the origin, then t → $\infty$, the controllability Gramian matrix can be computed as

$$W_{Co}=\underset{x\to \infty }{\lim }{W}_{Co}(t)={\displaystyle \underset{0}{\overset{\infty }{\int }}{\mathrm{e}}^{\tilde{A}t}\tilde{B}{\tilde{B}}^{\mathrm{T}}{\mathrm{e}}^{{\tilde{A}}^{\mathrm{T}}t}dt},$$

(34)

and this matrix satisfying the following Lyapunov matrix equation:

$$W_{Co}\tilde{A}+W_{Co}{\tilde{A}}^{\mathrm{T}}+\tilde{B}{\tilde{B}}^T=0.$$

(35)

Obviously, Equation (35) can see that the observability Gramian matrix W_Codepends on the input signal matrix $\tilde{B}$ and state matrix $\tilde{A}$; thus, the input/control energy could be influenced due to these matrixes.

In the case the system is controllable but a transient disturbance occurs, setting up the initial condition x(0) = x₀ at time t = 0; it is required to force a final state x₁ at final time t₁, there exists some input u(t) with t ∈ (−$\infty$, 0). The amount of energy required to use the input signal for reaching the initial state x₀, this energy can be compute as follows [44]:

$$E_{Co}(x_0)=\frac{1}{2}{\Vert u(t)\Vert }^2={\displaystyle \underset{-\infty }{\overset{0}{\int }}{u}^{\mathrm{T}}(t)u(t)dt}=x_0^{\mathrm{T}}{\left({\displaystyle \underset{0}{\overset{\infty }{\int }}{\mathrm{e}}^{\tilde{A}t}\tilde{B}{\tilde{B}}^{\mathrm{T}}{\mathrm{e}}^{{\tilde{A}}^{\mathrm{T}}t}dt}\right)}^{-1}{x}_0=x_0^{\mathrm{T}}{\left(W_{Co}\right)}^{-1}{x}_0.$$

(36)

From Equation (36), the input energy is a quadratic form in the initial state, and it associates with the inverse controllability Gramian matrix. In order to minimize the input/control energy over the entire set, we have to minimize the inverse controllability Gramian matrix or equivalently to maximize controllability Gramian matrix via Trace. Therefore, the optimal location of the sensor can be selected based on the optimization criterion is used as follows:

$${(J_{Co})}_i=\underset{j=1,2,3,\cdots ,m}{\max }{(trace(W_{Co}))}_j.$$

(37)

Figure 5b illustrates the flow chart of calculation steps needed to find the results of optimal location of sensor placement based on the proposed analytical expressions.

3.3.2. Determine Location of STATCOM-POD

As known, the states of a system are internal variables; in a general situation, it cannot to directly measure them, but at the same time, the outputs can be measured quite easily. The observability is the capability of the system to allow the reconstruction of the internal state variables through the time from the measured outputs. Therefore, this property plays an important role in the optimal location analysis of FACTS. Based on this property, a method is developed from the observability Gramian matrix of the reduced linearized system is proposed to determine an optimal location for STATCOM-POD placement.

For linear time invariant systems of Equation (22), the observability matrix can be defined as:

$$Ob=\left[\begin{array}{c}\tilde{C}\\ \tilde{C}\tilde{A}\\ \tilde{C}{\tilde{A}}^2\\ ⋮\\ \tilde{C}{\tilde{A}}^{n-1}\end{array}\right]\in {\mathcal{R}}^{(nm)\times n},$$

(38)

and the observability Gramian matrix of the pair ($\tilde{C},\tilde{A}$) at time t is

$$W_{Ob}(t)={\displaystyle \underset{0}{\overset{t}{\int }}\exp ({\tilde{A}}^{\mathrm{T}}\tau ){\tilde{C}}^{\mathrm{T}}\tilde{C}\exp (\tilde{A}\tau )d\tau }\in {\mathcal{R}}^{n\times n},$$

(39)

can be used to calculate the observability of a system. Therefore, there are two possible ways to check the observability of a system based on these matrices. Firstly, the pair ($\tilde{C},\tilde{A}$) is observable if and only if the observable matrix Ob has full rank n, i.e., rank (Ob) = n. Secondly, the pair ($\tilde{C},\tilde{A}$) is observable if and only if the observability Gramian matrix W_Ob(t) is positive definite for $\forall$t > 0, and it satisfies the following differential Lyapunov matrix equation:

$$\{\begin{array}{l}-\frac{d}{dt}{W}_{Ob}(t)+{\tilde{A}}^{\mathrm{T}}{W}_{Ob}(t)+W_{Ob}(t)\tilde{A}=-{\tilde{C}}^T\tilde{C},\\ W_{Ob}(0)=0.\end{array}$$

(40)

If the system in Equation (22) is asymptotically stable around the origin, then t → $\infty$, the observability Gramian matrix can be calculated:

$$W_{Ob}=\underset{x\to \infty }{\lim }{W}_{Ob}(t)={\displaystyle \underset{0}{\overset{\infty }{\int }}{\mathrm{e}}^{{\tilde{A}}^{\mathrm{T}}t}{\tilde{C}}^{\mathrm{T}}\tilde{C}{\mathrm{e}}^{\tilde{A}t}dt},$$

(41)

and it can be satisfied the following Lyapunov matrix equation:

$${\tilde{A}}^{\mathrm{T}}{W}_{Ob}+W_{Ob}\tilde{A}+{\tilde{C}}^T\tilde{C}=0.$$

(42)

Obviously, Equation (35) can see that the observability Gramian matrix W_Ob depends only on the output matrix $\tilde{C}$ and state matrix $\tilde{A}$. Thus, the output energy could be influenced due to these matrices.

If state x₀ ∈ $ℛ$, it exists the outputs y(t): 0 ≤ t < t₁, and assuming that the inputs u(t) = 0, are known; this is true, the above system in Equation (22) is observable over the interval [0, t1]. The corresponding output can be computed as follows [45]:

$$y(t)=\tilde{C}{e}^{\tilde{A}t}{x}_0,$$

(43)

and the output energy at the system is decayed from initial condition x₀ to zero in the inexistence of inputs, can be obtained as [44]:

$$E_{Ob}(x_0)=\frac{1}{2}{\Vert y(t)\Vert }^2={\displaystyle \underset{0}{\overset{\infty }{\int }}{y}^{\mathrm{T}}(t)y(t)dt}=x_0^{\mathrm{T}}\left({\displaystyle \underset{0}{\overset{\infty }{\int }}{\mathrm{e}}^{{\tilde{A}}^{\mathrm{T}}t}{\tilde{C}}^{\mathrm{T}}\tilde{C}{\mathrm{e}}^{\tilde{A}t}dt}\right)x_0=x_0^{\mathrm{T}}\left(W_{Ob}\right)x_0.$$

(44)

From Equation (44), the output energy is a quadratic form in the initial state, and it associates with the observability Gramian matrix. In order to maximize output energy over the entire set, we have to maximize the observability Gramian matrix via Trace. Therefore, the optimal location of STATCOM-FOD placement can be selected based on the optimization criterion is used as follows:

$$J_{Ob}=\underset{i=1,2,3,\cdots ,n}{\max }{(trace(W_{Ob}))}_i.$$

(45)

Figure 5a illustrates the flow chart of calculation steps needed to find the results of optimal location of STATCOM-POD placement based on proposed analytical expressions.

4. Case Study

In this section, the methods described in Section 3 are applied to determine the location of STATCOM-POD. The IEEE 24-bus reliability test system is used to investigate the effects of the proposed method for damping the power system oscillations. This network has two voltage levels that are 138 kV and 230 kV, consisting of a total 38 lines and 24 buses, and is model as described in Section 2.2. The data of IEEE 24-bus reliability test system is taken from [46]. All components are established in the PSS/E and MATLAB environments. The result of load flow calculation is shown in Figure 6.

In order to validate the proposed method for the optimal location of STACOM-POD with purpose for damping the power system oscillation, different simulations are carried out on the region of 230 kV voltage level of the IEEE 24-bus reliability test system through the study cases that are discussed as follows:

4.1. Applying the Conventional Method for Optimizing STATCOM-POD Placement and Parameters

The POD controller structure of STATCOM in Figure 3 is considered in this study. In order to seek the best sitting for the STATCOM for damping the oscillatory mode, the STACOM is placed in different buses on the test network except the generator buses. The model residue values associated with critical mode are computed by using the transfer function between the STATCOM active power deviation and its inputs. The parameters for the STACOM are listed in Appendix A.3. The numerical results of the model residues for the transfer function by applying the conventional method introduced in Section 3.2 are shown in

Table 1. Location indices of STATCOM-POD.

Model Residue Values of the Transfer Function
STATCOM-POD Placed at the Bus	$\left|R_i\right|$
14	0.976
11	0.778
17	0.736
12	0.198
24	0.012
19	0.005
20	0.001

. The test system is supposed that the critical oscillatory mode to the uncontrolled system characterized with the eigenvalue and the relative damping ratio are $\lambda =-0.1278\pm j3.5659$ and $\xi =0.029$, respectively. According to Table 1, the bus number 14 is the most effective placement of the STATCOM since this bus has the largest model residue value. The transfer function for the POD controller of STATCOM can be obtained as:

$$POD(s)=K_p(\frac{1}{1+0.1s})(\frac{10s}{1+10s}){(\frac{1+0.1359s}{1+0.4355s})}^2.$$

(46)

The value K_p in Equation (46) is calculated by supposing the eigenvalue moves to the desired location $\lambda =-0.725\pm j3.621$ from the original location $\lambda =-0.1278\pm j3.5659$, and can be obtained as:

$$K_p=\left|\frac{{\lambda }_i^{des}-{\lambda }_i}{R_{i}PO{D}_1({\lambda }_i)}\right|=\mathrm{26.157.}$$

(47)

In order to check the efficiency and reliability of the POD controller for STATCOM to the system stability, a three-phase fault is introduced in the line 11–14 close to the bus 14 at 1.0 s and cleared at 1.2 s after by opening the faulted line. A direct comparison between the relative rotor angle swing response of generators 7 and 13 to a three-phase fault with and without POD plotted in Figure 7. It is observed from this figure that the oscillations are damped out after the placement of STATCOM with the POD.

4.2. Applying the Proposed Method for Optimizing STATCOM-POD Placement

The IEEE test network of 24-bus is considered as a larger tested network. This network created the state matrix A. this matrix has dimension [126 × 126] including three states as the transducer, washout, and lead/lag of a STACOM added in the system. The parameters for the STACOM and POD and the input signal for controller are introduced in Section 4.1.

As mentioned in Section 3.1, the order reduction method to eliminate the singular values that are smaller than 10⁻³ (the threshold value chosen in this study) is implemented. A direct comparison between the eigenvalue of the system state matrix before and after reducing order is given in Figure 8a. According to this figure, the state matrix of the system after reducing order has the dimension [25 × 25]. The new system is only 25 Hankel singular values that are smaller than 10⁻³ after implementing the order reduction, as shown in Figure 8b. In addition, the result shown in Figure 8c is a comparison of the frequency responses between original and reduced systems. There is no difference between two input–output in the bandwidth from 10^−∞ to 10⁶. While the bandwidth from 10⁶ to 10^∞, the frequency responses of the reduced system is flat, this means the order of the new system is smaller than that of original system. Therefore, the new system now has 25 state variables, which can fully meet effects of the system to calculate the matrices, and we can conclude that the system before and after reduction is equivalent in terms of robustness.

4.2.1. Determining the Location of the STATCOM-POD Placement

As mentioned in Section 3.3.2, the optimal location of the STATCOM-POD placement is determined based on the objective function in Equation (45) by using the proposed algorithm, and the flowchart shown in Figure 5a, the procedure is given below:

First.	Start.
Step 1.	Place the STATCOM-POD controller at the bus ith (i = 1, 2, 3, …, n).,
Step 2.	Run base case load flow,
Step 3.	Generate the active power perturbation signal in the transmission line jth. (j = 1, 2, 3, …, m), in which m is large enough. For each contingency of the active power perturbation j.
Step 4.	Compute the matrices A, B, C, and D corresponding to the placement ith and the active power perturbation signal jth.
Step 5.	Perform the order reduction based on the proposed method in Section 3.1.
Step 6.	Estimate the stable condition based on the state matrix ${\tilde{A}}_{ij}$ to exclude the unstable cases. If the condition is satisfied, Step 6 will be performed; conversely, Step 1 will be iterated.
Step 7.	Compute the observability Gramian matrix of the new system corresponding to the active power perturbation signal jth.
Step 8.	Iterate the Steps 3 to 7 for computing a set of the active power perturbation (j = 1, 2, 3, …, m).
Step 9.	Compute the energy for each placement corresponding to a set of the active power perturbation (j = 1, 2, 3, …, m).
Step 10.	Iterate Steps 1 to 7 to calculate all the STATCOM-POD placements (i = 1, 2, 3, …, n).
Step 11.	Compare the maximum total energy to evaluate the optimal placement for STATCOM-POD controller based on the objective function in Equation (45).
Finally.	End.

A set of the active power perturbation was selected from single-line outage cases based on the real power flow performance index (PI) introduced in [47] and successfully applied in the literature [2,11,25,48,49]. It can be obtained as:

$$\mathrm{PI}={\displaystyle \sum _{i=1}^{nl}\frac{{\upsilon }_i}{2\vartheta }{(\frac{P_{li}}{P_{li}^{rated}})}^{2\vartheta }},$$

(48)

where nl denotes total number of lines, $\vartheta$ denotes the exponent, ${\upsilon }_i$ denotes a the non-negative real weighting factor used to reflect the importance for line i, and $P_{li}$ and $P_{li}^{\max }$ are the real power and rated power flows of the line i, respectively. In this simulation, the value of ${\upsilon }_i$ and $\vartheta$ is selected equal to 1.0 and 2.0, respectively. From that, we have a set of the active power perturbation to consider as the contingency cases listed in

Table 2. Studied contingency cases.

Case Number	Contingency Case (Active Power Perturbation Signal in the Line between Buses)	Case Number	Contingency Case (Active Power Perturbation Signal in the Line between Buses)
(1)	11–13	(10)	16–17
(2)	11–14	(11)	16–19
(3)	12–13	(12)	17–18
(4)	12–23	(13)	17–22
(5)	13–23	(14)	18–21
(6)	14–16	(15)	19–20
(7)	15–16	(16)	20–23
(8)	15–21	(17)	21–22
(9)	15–24	-	-

.

The energy values based on Trace of observability Gramian matrix corresponding to the contingency cases as listed in Table 2, which were determined by using the proposed algorithm as shown in Figure 5a, are given in

Table 3. Energy values according to several feasible locations of STATCOM-POD controller.

Case Number	Contingency Case (From Table 2)	Energy Values Based on Trace of Observability Gramian Matrix
		STATCOM-POD Placed at the Bus
		11	12	14	17	19	20	24
(1)	11–13	82.58	62.90	89.57	294.03	153.84	112.79	157.13
(2)	11–14	494.38	277.40	427.34	680.79	503.95	382.72	373.93
(3)	12–13	170.41	97.55	136.50	229.47	149.17	114.52	144.60
(4)	12–23	251.29	223.05	221.48	356.44	289.13	311.12	210.69
(5)	13–23	234.69	177.58	221.68	404.43	261.51	268.73	225.35
(6)	14–16	618.80	218.62	691.30	619.27	434.24	360.01	333.03
(7)	15–16	299.01	235.69	387.10	642.85	452.45	307.92	366.29
(8)	15–21	469.97	387.73	498.84	818.31	579.79	419.27	609.25
(9)	15–24	308.28	205.19	272.09	447.37	361.44	298.43	397.58
(10)	16–17	491.35	305.15	573.04	852.46	656.52	463.65	490.95
(11)	16–19	428.88	321.11	548.22	1020.29	855.97	600.79	562.90
(12)	17–18	291.69	230.46	337.48	666.07	441.14	310.05	338.04
(13)	17–22	271.33	292.01	278.07	454.16	295.42	209.23	237.74
(14)	18–21	104.00	94.42	109.06	255.49	140.14	106.10	169.33
(15)	19–20	359.08	282.60	545.56	1071.40	760.24	523.01	611.26
(16)	20–23	300.10	208.97	529.71	1108.86	776.81	491.31	646.05
(17)	21–22	356.14	295.10	363.93	557.35	382.87	285.16	366.80
Total		5531.97	3915.53	6230.94	10479.04	7494.62	5564.79	6240.92

, where we only listed out several feasible locations (the cases have the high energy values). It can be seen in

Table 4. Scenarios on three-phase fault for test of transient stability.

Scenario Number	Fault Is nearly Bus	Start Time	Closed Time	Line Outage between Buses
(1)	14	1	1.2	14–11
(2)	15	1	1.1	15–21
(3)	17	1	1.1	17–22
(4)	12	1	1.2	12–13

that the bus 17 is the optimal location for STATCOM-POD controller since this location has the highest total energy value compared to other locations; in other words, the energy need to drive the observable state variables is high. Therefore, the bus 17 is the best location to install the controller.

In order to verify the above conclusion, the transient stability analysis of the system is proposed to compare the optimal location with other feasible ones under considering fourth three-phase fault scenarios as shown in Table 4.

• The scenario number 1: A solid three-phase fault occurs at 1 s on line 14–11, close to bus 14, and is cleared after 0.15 s by tripping the faulted line.
• The scenario number 2: A solid three-phase fault occurs at 1 s on line 15–21, close to bus 15, and is cleared after 0.15 s by tripping the faulted line.
• The scenario number 3: A solid three-phase fault occurs at 1 s on line 17–22, close to bus 17, and is cleared after 0.15 s by tripping the faulted line.
• The scenario number 4: A solid three-phase fault occurs at 1 s on line 12–13, close to bus 12, and is cleared after 0.15 s by tripping the faulted line.

The simulations have been performed for each placement of STACOM corresponding to the diminishing total energy value through the following cases:

Case 1. The simulation was done on scenario number 1, when the STATCOM_POD is placed at bus 17, the transient respond improves and the oscillations of the relative angle between generators 7 and 13 decreases significantly compared with that is not placed at bus 17, as can be observed from Figure 9.

Case 2. In Table 3, it can be seen that buses 17, 14, and 12 have maximum, third, and seventh highest total energy value, respectively. A simulation was done on scenario number 2, when STATCOM-POD is placed at these buses. Figure 10a plots the relative angle between generators 15 and 13; it can be that the transient response improves significantly with STATCOM-POD placed at bus 17. Another simulation was done, and the result is plotted in Figure 10b: the relative angle between generators 16 and 13 is also significantly improved when STATCOM-POD placed at bus 17.

Case 3. Simulation was done on scenario number 3 considering the placement of STATCOM-POD as done in Case 2. Figure 10a,b shows the plot of the relative angle between generators 18 and 22 compared to generator 13. In this case, the STATCOM-POD is placed at bus 17, damping out oscillations of rotor relative angle faster than if placed at buses 14 and 12.

It is also observed in Figure 10 and Figure 11 that the damping of rotor angle oscillations improves with STATCOM-POD placed at bus 14 and this placement provides better damping compared to bus 12.

Case 4. In order to make sure “the higher the total energy value, the quicker the damping oscillations”, the simulation in this case was done on scenario number 4, when STATCOM-POD placed at buses 17, 14, and 19. Generator 13 oscillating against generators 1 and 7 was plotted, as shown in Figure 9 and Figure 12. With a fault occurred, the damping ratio is low, so that the damping of oscillations of relative angle between generators 13 and 1 is very slow without STATCOM-POD, as shown in Figure 9.

As can be observed in Figure 12, the transient response improves and oscillations relative angle between generators 1 and 13 decreases significantly with STATCOM-POD is placed at the bus 17 compared with that is placed at buses 14 and 19. In addition, the rotor angle swing response of generators 13 and 1 damped out significantly with STATCOM-POD placed at bus 19 compared with that placed at bus 14. Thus, we could conclude that the location of the STATCOM-POD placement with the purpose for damping the power system oscillations is more reliable on the high observable energy value. Therefore, bus 17 is best location to install the STATCOM-POD controller.

4.2.2. Determining the Location of the Input Signal Placement for STATCOM-POD

After determining the a best location to install STATCOM-POD, the optimal location of feedback signal placement needs to consider for FACTS controllers, it is necessary condition to enhance the system stability.

As mentioned in Section 3.3.1, the optimal location of the signal placement for STATCOM-POD is determined based on the controllable function in Equation (37) by using the proposed algorithm, and the flowchart shown in Figure 5b, the procedure is given below:

First.	Start.
Step 1.	Place the STATCOM-POD controller at the optimal bus i (at bus 17 in this study).
Step 2.	Run base case load flow.
Step 3.	Generate the active power signal in the transmission line jth (j = 1, 2, 3, …, m).
Step 4.	Compute the matrices A, B, C, and D corresponding to the placement i and the active power signal jth.
Step 5.	Perform the order reduction based on the proposed method in Section 3.1.
Step 6.	Estimate the stable condition based on the state matrix ${\tilde{A}}_{ij}$ to exclude the unstable cases. If the condition is satisfied, this step will be performed; otherwise, Step 1 will be iterated.
Step 7.	Compute the controllability Gramian matrix of the new system corresponding to the active power perturbation signal jth.
Step 8.	Iterate Steps 3 to 7 for computing a set of the active power perturbation (j = 1, 2, 3, …, m) corresponding to the i placement.
Step 9.	Compute the energy for a set of the active power perturbation (j = 1, 2, 3, …, m).
Step 10.	Evaluate the maximum energy in order to determine the optimal location for the input signal placement based on the objective function in Equation (37).
Finally.	End.

The values based on Trace of controllability Gramian matrix, corresponding to a set of input signal, are considered on 17 lines of 230 kV level of IEEE 24-bus network, as given in

Table 5. Energy values according to a set of input signal location with STATCOM-POD placed at bus 17.

Case Number	Active Power Signal in the Line between Buses	Trace Value of Controllability Gramian Matrix	Index
(1)	14–16	622.40	0.49
(2)	19–16	1030.32	0.82
(3)	18–17	617.00	0.49
(4)	17–16	984.65	0.78
(5)	11–14	645.46	0.51
(6)	15–21	794.08	0.63
(7)	15–16	607.81	0.48
(8)	11–13	486.66	0.39
(9)	12–13	406.34	0.32
(10)	13–23	395.03	0.31
(11)	18–21	266.74	0.21
(12)	12–23	395.03	0.31
(13)	17–22	342.55	0.27
(14)	21–22	595.58	0.47
(15)	23–20	1261.78	1.00
(16)	15–24	566.31	0.45
(17)	19–20	1058.28	0.84

. These values were determined by using the proposed algorithm as shown in Figure 5a. It can be seen in Table 5 that the line between buses 23 and 20 is the optimal location of input signal placement for STATCOM-POD controller since this location has the highest value of Trace of controllability Gramian matrix compared to another locations, in other words, the energy need to drive the controllable state variables is small. Therefore, based on the result of the Trace of controllability Gramian matrix, the line between buses 23 and 20 is the best location to be used to install the sensor for receiving and sending the real power flow signal in the transmission line (has been taken as the feedback signal) to the POD of STATCOM.

In order to test the system transient stability ability with the received feedback signal from the line having the highest trace value for POD of STATCOM, a three-phase fault was applied on Scenario 3. Comparing between the optimal location (line 23–20) and another location (line 18–21) having minimum energy value based on the output reactive power of controller unit and the rotor angle swing response of between generators 13 and 7 is done. Figure 13 show the plot the transient response of the relative angle between generators 7 and 13 and the output reactive power of STATCOM-POD.

Observing Figure 13a, when the sensor (the device to be used to receive the real power flow from the line and then send it to the POD) is placed in line 23–20, the rotor angle swing response of between generators 13 and 7 is significantly improved compared with that in the line 18–21.

A direct comparison of the output reactive power of controller unit was done to observe the capable of meeting of controller unit. Figure 13b plots this comparison. Evidently, it can be observed from this figure that the sensor is placed in the line 23–20 (POD controller received the active power signal in the line 23–20); the output response of controller unit is quicker than that in the line 18–21. As a result, we could conclude that the location of the input signal placement for STATCOM-POD with the purpose to dampen the power system oscillations is a lot more reliant on the high controllable energy value. Therefore, line 23–20 is best location. This line active power flow has been taken as the feedback signal for supplement controller of STATCOM.

4.3. Comparative Simulation

After applying the proposed method based on the small-signal analysis and retesting one based on the transient analysis, the best location for STATCOM-POD placement is bus 17 and the best location for sensor placement is line 23–20 (POD controller received the active power signal in line 23–20). In addition, according to Table 3, the second-highest total energy value was obtained on the electric power supply system when STATCOM is placed at the bus 19. Vice versa, applying the conventional method as explained and applied in Section 3.1 and Section 4.1, respectively, the best location for STATCOM-POD placement is bus 14.

To compare the meeting capable of controller unit for selected location by two methods with main objective to validate further the effectiveness of the proposed method, a direct simulation was done on Scenario number 2, when STATCOM-POD is placed at three locations, buses 14, 17, and 19, based on:

(i)

The active power oscillations damping on the transmission line;

(ii)

The oscillations damping of the relative angle between generators 7 and 13; and

(iii)

The output reactive power response of controller unit.

Figure 14 show the plot of the transient response on Scenario number 2, in which Figure 14a shows the plot of the active power flow in line 1–2. The oscillations are dampened out after STATCOM-POD placement at bus 17 in about 6 s. With STACOM-POD placed at bus 19, the oscillations are dampened out in about 10 s, whereas when STACOM-POD is placed at bus 14, the oscillations are dampened out in about 18 s.

The generator 13 oscillates against generators 1, 2 and 7. With the occurred fault, the critical mode, which damping ratio is very low and is an inter-area mode, gets excited. Because this, the oscillations damping of the relative angle between generators 7 and 13 is slow without STACTCOM-POD, as shown from Figure 9. Figure 14b shows the plot of the oscillations of the relative angle between generators 1 and 13. It can be observed from this figure that the oscillations are damped out in about 6 s, 10 s, and 18 s with STATCOM-POD placed at buses 17, 19, and 14, respectively.

Figure 14c shows the plot of the dynamic response of STATCOM-POD unit corresponding to three locations. Observing this figure, injects to the gird with a reactive power range of −2.0 p.u. to 2.0 p.u. during the period of fault, this output reactive power flow reacted immediately, damped faster, and extremely fast returned to the compensative mode to reach new steady-state after clearing the fault when STATCOM-POD is placed at bus 17. Evidently, it can be observed through above simulation cases, the optimal STATCOM_POD placement and the optimal local input signal placement for POD have been sought by the proposed method provided the effectiveness and robustness of the performance under transient conditions better compared with the conventional method.

5. Conclusions

In this paper, a highly relevant stochastic method for seeking the optimal STATCOM-POD placement and the best location of the sensor (the device using to receive the real power flow from the line and send it to the POD) in the multimachine systems has been proposed with the purpose to dampen the rotor angle oscillations and enhance the transient stability. It is mainly focused on calculation the Gramian matrix of the controllability and observability of the linearized multimachine system based on small-signal analyses. The parameters of the controller obtained from the conventional method. On that basis, the proposed method aims to develop in two steps. Firstly, the optimal STATCOM-POD placement is determined by the algorithm based on the Gramian matrix of the observability on the small-signal stability analyses. Then, test cases of the transient stability are analyzed. The obtained result is no satisfactory since the POD of STACTCOM received the feedback local signal with placement that is not optimal location. In order to overcome this problem, the second step of the proposed method is developed to select the optimal location of the sensor based on the Gramian matrix of the controllability; it was also analyzed based on the small-signal stability and test cases of the transient stability.

The simulation results show that the rotor angle oscillations of generators are significantly dampened when the STATCOM-POD is placed at bus 17 and the sensor is placed in the line between buses number 23 and 20, which have the highest Trace of the Gramian matrix of the controllability and observability. Therefore, it could be concluded that the location of the STATCOM-POD placement with the purpose for damping the power system oscillations is more reliant on the high observable and controllable energy values.

The reduction method based on Krylov-based model reduction solution has also been introduced in this paper for reducing the calculation time of Gramian matrix when dealing with the large-scale power systems.