4___AITI#7853_in press Advances in Technology Innovation, vol. 6, no. 4, 2021, pp. 235-250 A Recursive Least-Squares Approach with Memorizing Factor for Deriving Dynamic Equivalents of Power Systems Ali Karami * Faculty of Engineering, University of Guilan, Rasht, Iran Received 06 June 2021; received in revised form 01 August 2021; accepted 02 August 2021 DOI: https://doi.org/10.46604/aiti.2021.7853 Abstract In this research, a two-stage identification-based approach is proposed to obtain a two-machine equivalent (TME) system of an interconnected power system for transient stability studies. To estimate the parameters of the equivalent system, a three-phase fault is applied near and/or at the bus of a local machine in the original multimachine system. The electrical parameters of the equivalent system are calculated in the first stage by equating the active and reactive powers of the local machine in both the original and the predefined equivalent systems. The mechanical parameters are estimated in the second stage by using a recursive least-squares estimation (RLSE) technique with a factor called “memorizing factor”. The approach is demonstrated on New England 10-machine 39-bus system, and its accuracy and efficiency are verified by computer simulation in MATLAB software. The results obtained from the TME system agree well with those obtained from the original multimachine system. Keywords: network reduction, transient stability, critical clearing time, system identification 1. Introduction In modern power systems, for economic and security reasons, individual power companies are connected together with tie lines to form an interconnected system. Sophisticated computer programs are usually used for analysis and control of large power systems. However, it is not necessary to model the whole power network when the aim is to investigate a specific phenomenon in a small part of the system. In addition, sometimes it is very difficult to model the entire system in detail due to the lack of complete system data and/or inaccuracies in the system data [1]. To reduce the computational requirement, it is highly desirable to model only a small part of the system called “internal system” with enough accuracy, and represent the rest part of the system known as “external system”, by an equivalent model. Power system equivalencing (or reduction) is generally divided into two main groups, namely static equivalencing and dynamic equivalencing. The static equivalents are employed in investigating the steady-state problems of power systems such as power flow and system planning studies. The dynamic equivalents are used in dynamical studies of power systems, including large and/or small disturbance stability problems as well as dynamic security assessment. This study is concerned with dynamic equivalencing of power systems for transient stability studies [2-3]. To reduce the computational burden of transient stability analysis, a dynamic equivalent is sometimes used for multimachine power systems [4]. However, a single-machine-infinite-bus (SMIB) system is usually employed for this purpose. An infinite bus is used to approximately model a large area in a power system [1]. An equivalent synchronous machine rather than an infinite bus was considered in the work of Li et al. [5] to include the external system effect on the internal system. This early work was the main motivation for the present study. However, the method proposed in this study differs significantly from the method of this early work. * Corresponding author. E-mail address: karami_s@guilan.ac.ir Tel.: +98-13-33690274; Fax: +98-13-33690271 Advances in Technology Innovation, vol. 6, no. 4, 2021, pp. 235-250 The main objective of this study is to derive a two-machine equivalent (TME) system of a multimachine power system, by looking into the system from the bus of a local machine. The desired parameters of the defined equivalent system can be separated into two groups. The first group consists of the electrical parameters and the second group consists of the mechanical parameters of the TME system. To estimate the parameters of the equivalent system, a systematic two-stage identification-based approach is proposed. The electrical parameters of the equivalent system are calculated in the first stage of the approach by using an innovative and simple method. The mechanical parameters of the equivalent system are estimated in the second stage of the approach by employing a recursive least-squares estimation (RLSE) technique with a factor called here “memorizing factor”. One advantage of the proposed methodology is that it is simple to understand. In addition, it is relatively easy to implement in a computer program. The effectiveness of the proposed methodology is demonstrated on the New England 10-machine 39-bus test system. The simulation results obtained are thoroughly described and discussed. The rest of the study is organized as follows. Section 2 presents the literature review. Section 3 represents a brief description of the RLSE algorithm. The multimachine power system model and problem statement are described in section 4. Section 5 presents the first stage of the proposed approach. The second stage of the proposed approach is described in section 6. The simulation results are presented in section 7. Section 8 discusses the performance of the approach, and finally section 9 is devoted to the concluding remarks. 2. Literature Review Research in power system reduction dates back to several decades ago; however, only a few references are mentioned here due to space limitation. The most well-known methods of static equivalencing are Ward equivalent [6] and radial equivalent independent (REI) equivalent [7]. The Ward equivalent is based on Gaussian elimination technique, and the REI equivalent is indeed an extension of the Ward equivalent. These methods separate a power system into three parts: internal, boundary, and external systems. The size of power system model is decreased by replacing the external system with small equivalents, whereas the internal system remains unchanged [1]. An extended Ward equivalent method has been used for power system contingency analysis and static security assessment in the work of Srivani et al. [8]. Fu et al. [9] have proposed a hybrid method based on artificial neural network (ANN) and Ward-type equivalency for fast and online voltage security assessment of power systems. The proposed approach possesses good properties of Ward equivalent method, and can update the parameters of the equivalent model for representing real-time topology change of the external system. Lee et al. [10] have presented a static network equivalent model for Korean power systems. The proposed equivalent model preserves the overall transmission network characteristics focusing on power flows among the areas in Korean power systems. The already existing methods of dynamic equivalencing are modal methods, coherency methods, and measurement or simulation-based methods [11]. Modal methods are based on the linearized state-space model of power systems. In the modal methods, the system equations are first linearized. Then, the eigenvalues, eigenvectors, and participation factors of the linearized equations are employed to eliminate the less relevant part of the system [12]. Paternina et al. [13] have proposed a strategy for developing dynamic equivalents in large-scale power system by mode preservation. Using the state-space representation of power systems, other reduction techniques have also been investigated in the literature. Rergis et al. [14] have proposed a general approach for constructing reduced-order models of power systems, which is based on a truncated balanced realization (TBR) linear reduction procedure. Chaniotis et al. [15] have described the use of Krylov subspace methods in the model reduction of power systems. Zhu et al. [16] have proposed a method for large-scale power system model reduction based on the extended Krylov subspace technique. Ritschel et al. [17] have reported a balanced truncation model reduction approach for reducing two commonly used nonlinear and dynamical models of power grids. 236 Advances in Technology Innovation, vol. 6, no. 4, 2021, pp. 235-250 In coherency methods, coherent groups of generators are first identified by using a suitable technique. Then, the generators in coherent groups are represented by an equivalent generator [18]. Unlike the modal methods, the coherency methods retain the physical models of the generators in an equivalent form [1, 11]. Ayon et al. [19] have presented a methodology to identify coherent areas based on sliding power spectral density algorithm. Miah [20] has proposed a coherency-based dynamic equivalent of a power system external area for transient stability assessment. Based on the use of coherency aggregation techniques and nonlinear black-box optimization, Matar [21] has developed a dynamic model reduction methodology. Chittora et al. [22] have proposed a coherency based dynamic reduction technique for reducing the size of the power system network, and have tested it on a network having high voltage direct current (HVDC) transmission line. In measurement or simulation-based methods, the external system response is either measured or simulated, and system identification or curve-fitting techniques are used to determine the model parameters [11]. In the system identification based methods, a structure is first defined for the equivalent system. The measured data are then utilized to estimate the parameters of the predefined equivalent system. These methods require no knowledge of the parameters and topology (structure) of the external system. Stankovic et al. [23] have proposed an ANN-based strategy for the identification of the reduced-order dynamic equivalents, which uses only measurements at points where the internal and external systems are interfaced. Chakrabortty et al. [24] have described an algorithm for constructing simplified interarea models of large power networks, using dynamic measurements available from phasor measurement units (PMUs). Ramirez et al. [25] have proposed a robust dynamic reduction approach to reduce the computational burden and time of transient stability studies. Rueda et al. [26] have introduced an identification approach of the dynamic equivalent parameters, from measured operational dynamic responses associated to different disturbances. In order to increase both model accuracy and simulation speed, Zhang et al. [27] have proposed a measurement-based dynamic equivalent. Tong et al. [28] have proposed an artificial neuro-fuzzy inference system (ANFIS)-based dynamic reduction approach, in which the wide-area measurements obtained by PMUs at the boundaries between the reduced system and the study system have been used to represent the external area. Jiang et al. [29] have proposed a measurement-based reduction approach using system identification techniques, in which instead of a single event, multiple grid events with different patterns have been considered for estimating the parameters of the equivalent loads. Almeida et al. [30] have reported the development of three software tools for the computation of dynamic equivalent models of power systems. The available techniques for power system dynamic reduction suffer from two major drawbacks: (1) they are usually time-consuming; (2) they are difficult to implement in a conventional transient stability program. For instance, as noted in the work of Ayon et al. [19], the eigenvalues and eigen-matrix of the linearized system equations are usually difficult and also time-consuming to compute in the modal based methods. 3. RLSE Algorithm with Memorizing Factor Assume that a dependent variable y can be evaluated by a linear function of K independent variables denoted by a row vector x = (x1, x2, …, xK) [31], then Eq. (1) can be obtained: 1 1 2 2 ( x ) ... x Φ K K y f x x xφ φ φ= = + + + = (1) where Φ = ���,��,… ,� � shows an unknown column vector of parameters to be estimated, and T denotes transpose. Assume further that S samples of measurements y(t) and x1(t), x2(t), …, xK(t) are available at time t = 1, 2, …, S. In the least-squares (LS) technique, the aim is to fit the S samples by the linear function of Eq. (1). By filling the samples into Eq. (1), a set of linear equations is obtained and can be written in matrix form as: Y X Φ= (2) 237 Advances in Technology Innovation, vol. 6, no. 4, 2021, pp. 235-250 where Y denotes an S × 1 vector and X shows an S K× matrix. Note that if S > K, that is for the case in which the number of samples is greater than the number of unknown parameters (i.e., overdetermined system of equations), there exists no exact solution for the K unknown parameters Φ. Notice that sometimes a time-varying process is encountered and it is required that the estimated parameters are updated when new observations (samples) are available. To improve the computational efficiency, RLSE technique can be employed in such cases. In RLSE algorithm, the previously estimated parameters are utilized to obtain new estimations for the unknown parameters. To reduce the effect of past data points on the current estimate, a forgetting factor (� ) is used in the LS criterion. The forgetting factor � is a number less than 1 (e.g., � = 0.99) and causes the recent measured data to affect the estimation more, i.e., past samples are discounted. The RLSE updating procedure with a forgetting factor � can be summarized as follows [31]: P ( ) x ( 1) G ( 1) [ x ( 1) P ( ) x ( 1)] T T f t t t t t tλ ++ = + + + (3) Φ( 1) Φ( ) G( 1)[ ( 1) x( 1)Φ( )]t t t y t t t+ = + + + − + (4) P ( 1) [P ( ) G ( 1) x ( 1)P ( )] / f t t t t t λ+ = − + + (5) where G(t+1) denotes an K × 1 vector and P(t) is an K × K matrix. The matrix P is usually initialized to a very large identity matrix since no information is available in the beginning of RLSE algorithm. The initial values of unknown parameters, that is Φ(0), can be obtained by LS method or can be simply assumed as a vector of zeroes. In some cases, however, a time-varying process may be encountered, in which the older samples convey more information about the “true” parameter vector Φ than the newly available data points. In such cases, more weights are to be given to the older samples so that the parameter values are weighted more by the older samples. To fix this problem, a factor called “memorizing factor” (��) is introduced in this study. The idea of a memorizing factor is very simple. It can be shown that if a number of more than 1 is set to the forgetting factor � used in Eqs. (3)-(5), then rather than the older data points, less weight is given to the newly available samples. In other words, if the forgetting factor � is set to a value just more than 1, then compared with the older samples, the recent samples are discounted. This forgetting factor of more than 1 is called here the memorizing factor ��. Notice that the greater the memorizing factor is, the less the influence of the recent data points on the estimation algorithm is. It is to be noted that the equations for RLSE method with a memorizing factor �� are almost same as Eqs. (3)-(5). The only difference is that the forgetting factor � is replaced by the memorizing factor �� in those equations. As will be seen later on in this study in section 7.1., the parameters of the TME system are estimated with much more accuracy if a memorizing factor (e.g., �� = 1.002) is used in the estimation algorithm. 4. Power System Model and Problem Statement Consider that there is a power system consisting of n generators. In the classical model, the ith generator is represented with internal electromotive forces (EMF) ��� = ��∠�� behind direct-axis transient reactance ��� � as shown in Fig. 1. In Fig. 1, it is assumed that a local machine i is connected to the rest of the power system. The aim is to derive a TME system for the entire power system, by looking into the system from the bus of this local machine. The loads are assumed to be represented as constant shunt admittances in the classical model of power system. The reduced (internal) admittance matrix Yred (or Yint) of the system is then obtained [32]. 238 Advances in Technology Innovation, vol. 6, no. 4, 2021, pp. 235-250 i i E δ∠ dixj ′iP iQ Fig. 1 The connection of a local machine to a multimachine system Letting ��� ≅ �� − �� and denoting the elements of matrix ���� as ��� = ��� + � ��, the active power and reactive power of the ith machine can be expressed as follows [32]: 2 1 ( cos sin ) n j i i i ii i j ij ij i j ij ijP E G E E G E E Bδ δ = ≠ = + +∑ (6) 2 1 ( sin cos ) n j i i i ii i j ij ij i j ij ijQ E B E E G E E Bδ δ = ≠ = − + −∑ (7) The dynamics of the synchronous machine in classical model are governed by the following differential equations [32]: ; with 1, 2, ..., i i s d i n dt δ ω ω= − = (8) 1 [ ( )] ; w ith 1, 2, ...,i mi i i i s i d P P D i n dt M ω ω ω= − − − = (9) where !�� is the input mechanical power (assumed constant), "� is the rotor inertia constant, and #� is the natural damping coefficient of the ith machine. Also, %& is the electrical synchronous speed. To simulate the behaviour of the actual system in the event of a large disturbance, it is necessary to numerically solve Eqs. (8) and (9) in both faulted period and post fault period. Having obtained the time variation of the rotor angles, the time variation of the ith machine’s active power (!�) and reactive power ('�) can be easily obtained by using Eqs. (6) and (7), respectively. As can be seen in Eqs. (6) and (7), both the angle of ith machine itself and the angle of remaining system machines are used in calculating !� and '�. It is also required to use the elements of the reduced matrix ����. It is important to mention here that the information of the whole system is implicitly included in the expressions for !� and '�. In other words, !� and '� are indeed global system variables, in which the effects of both the machine itself and the rest of the network are included. Suppose that a large disturbance occurs near and/or at the bus of local machine i. Furthermore, assume that it is only required to approximately evaluate the transient stability of this machine with respect to the rest of the system. A dynamic equivalent model of a multimachine system, viewed from the bus of this local machine, is derived for this purpose. Here, machine i represents the internal system (of interest), and the rest of the system represents the external system (to be reduced). Instead of an infinite bus, the aim is to replace the external system by an equivalent machine e in series with a lossless transmission line of reactance ��, as depicted in Fig. 2. Fig. 2 shows the defined TME system. The equivalent machine e shown in Fig. 2 is represented as a classical machine with an internal EMF ��� = ��∠��, a rotor inertia constant "�, and a natural damping coefficient #�. Direct-axis transient reactance of the equivalent machine e is included in the reactance of transmission line �� for the sake of simplicity. 239 Advances in Technology Innovation, vol. 6, no. 4, 2021, pp. 235-250 , iM iD iP iQ dixj ′ ejx , eM eD i i E δ∠ e e E δ∠ Fig. 2 The TME system of a power system To determine the parameters of the TME system, the original multimachine system is initially excited to obtain the time variations of several system variables. A three-phase fault of short duration (e.g., 50 ms) is considered at the terminal of the local machine for this purpose. The obtained time solution for the system variables is then saved. Notice that the two-stage approach is performed only in post fault period. For practical application of the proposed approach, the required time variations of the system variables can be obtained from real system events. The desired parameters of the TME system can be separated into two groups. The first group consists of three electrical parameters (i.e., ��, ��, and ��); whereas the second group includes "� and #�, representing the mechanical parameters of the equivalent machine e. It is to be noted that all parameters of the local machine i are assumed to be known. 5. The First Stage of the Proposed Approach To estimate the parameters of the TME system, a three-phase short circuit is applied near the terminal of machine i in the original multimachine system. In other words, a time-domain simulation is carried out in the original test system. From the simulation, the active and reactive powers of machine i (i.e., !� and '�), are obtained at each time-step in post fault period. The time variation of local machine angle (i.e., ��), is also obtained from the time-domain simulation. The electrical parameters of the TME system are then estimated as follows. The active and reactive powers of machine i in Fig. 2 can be calculated by using the following equations [1]: sin( )i ei i e edi E E P x x δ δ= − ′ + (10) [ cos( )]ii i e i e edi E Q E E x x δ δ= − − ′ + (11) The electrical parameters ��, ��, and �� can be estimated by equating the active and reactive powers of machine i, obtained from the original system, with the corresponding powers calculated from Eqs. (10) and (11). In other words, those electrical parameters are calculated in such a way that if they are substituted into Eqs. (10) and (11), the active and reactive powers of machine i from the TME system become equal to their corresponding powers from the original system. There are only two known variables (i.e., !� and '�) here, but three unknown parameters (i.e., ��, ��, and ��) need to be estimated. Therefore, many solutions can be found for the unknown parameters. To obtain suitable solutions, the rotor angle �� of machine e is assumed to be equal to the inertia weighted average (iwa) angle of all machines in the external system. The angle ��() is defined as: 1 1 n iwa j j j it M M δ δ = ≠ = ∑ (12) where "* = ∑ "� , �-�.� , and n is the number of generators. Replacing the assumed value for �� from Eq. (12) into Eqs. (10) and (11), one can find �� and �� by the following equations (after a little algebra): 240 Advances in Technology Innovation, vol. 6, no. 4, 2021, pp. 235-250 [ cos( ) sin( )] i i e i i e i i e E P E P Qδ δ δ δ = − + − (13) sin( )i ee i e di i E E x x P δ δ ′= − − (14) It should be emphasized that, at each time-step of integration, it is required to re-calculate the parameters ��, ��, and �� by repeatedly solving Eqs. (12)-(14). Then, the time solutions of ��, ��, and �� are saved. The values obtained for these parameters at the end of the first stage are considered their final estimated values. 6. The Second Stage of the Proposed Approach To estimate the values of "� and #�, it is assumed that the time variation of �� obtained from the first stage of the approach fits the following second-order swing equation: 2 2 e e e e me e d d M D P P dtdt δ δ + = − (15) where !�� represents the input mechanical power and !� denotes the electrical output power of machine e. Note that the defined TME system shown in Fig. 2 is purely reactive; therefore, the electrical output power of machine i is completely absorbed by the equivalent machine e. The equivalent machine e acts here as a synchronous motor, and this leads to the following equations: e iP P= − (16) me miP P= − (17) Substituting Eqs. (16) and (17) into Eq. (15), it yields: 2 2 e e e e i mi d d M D P P dtdt δ δ + = − (18) To estimate the values of "� and #�, the following approach is used. Since time variation of �� is known, time variations of the first and second time-derivatives of �� (i.e., � / � = 0��/02 and � 3 � = 0 ���/02 �) can be calculated by using suitable numerical differentiation techniques. The first and second time derivatives of eδ for the lth step in post fault period are denoted by �/��4 and � 3 ��4 , respectively. Substituting � / ��4 and � 3 ��4 into Eq. (18) gives: miieeee PlPlDlM −=+ )()()( ... δδ (19) where !��4 shows the electrical output power of machine i, for the lth step in post fault period. Notice that Eq. (19) represents a discrete and algebraic equation, which is also linear with respect to "� and #�. Also, Eq. (19) needs to be satisfied at each time-step in post fault. Therefore, at the end of the post fault, there will be a set of algebraic equations from which only two unknown parameters are to be estimated. A set of overdetermined linear algebraic equations needs to be solved here. In this study, the RLSE algorithm with a memorizing factor is employed to estimate the values of "� and #� (referring to section 3). Based on what mentioned so far, the main steps of the proposed method are summarized in the flowchart illustrated in Fig. 3. In Fig. 3, dt denotes the time-step used in the numerical integration, and T represents the total simulation time. 241 Advances in Technology Innovation, vol. 6, no. 4, 2021, pp. 235-250 Start Perform time-domain simulation to solve swing equations F irst sta g e o f th e a p p ro a c h Is t > T ? NO YES NO YES End S e c o n d sta g e o f th e a p p ro a c h Define a fault scenario in the power system Set t = t0 , where t0 is a few time after the fault clearing time Calculate by using Eq. (12), and then setiwaδ iwae δδ = E x c ite th e sy ste m b y a p p ly in g a fa u lt Obtain and by solving Eqs. (13) and (14)eE ex Set t = t + dt Save final estimated values of , , and eE eδ ex Obtain and for the post-fault period using numerical differentiation techniques eδɺ eδɺɺ Set t = t0 , where t0 is a few time after the fault clearing time Obtain and by employing an RLSE algorithmeM eD Set t = t + dt Is t > T ? Save final estimated values of and eM eD Read power system data and perform power flow analysis Fig. 3 The proposed two-stage equivalencing approach 7. Simulation Results The proposed two-stage approach is applied to the New England 10-machine 39-bus (or the IEEE 39-bus) test system. A single-line diagram of this system is shown in Fig. 4. The bus and line data of this system can be found in the work of Pai et al. [32]. The New England test system is called “the test system” hereafter, for the sake of simplicity. The time-domain simulations are carried out by using a computer program written by the author in MATLAB software. This program solves system differential equations by utilizing the fourth-order Runge-Kutta integration technique. However, power flow solutions for the test system are obtained by using the MATPOWER package [33]. A time-step of 0.001 s is used in all calculations. The classical model is used and a uniform damping (i.e., 567806 = #�/"�; with i = 1, 2, …, n) of 1.0 is assumed for all machines of the test system. Fig. 4 Single-line diagram of the New England 39-bus test system 242 Advances in Technology Innovation, vol. 6, no. 4, 2021, pp. 235-250 7.1. Main results of simulation To excite the test system, a three-phase self-clearing type fault of 0.05 s duration is considered at bus 35. The simulation is carried out for 10 s in order to capture full dynamic behaviour of the test system. The location of the fault is very close to generator 2, as shown in Fig. 4. For this fault, the internal system consists of generator 2 and the rest of the system represents the external system. Generator 2 is the local machine here (i.e., i = 2), and the aim is to obtain a TME system from the bus of this local machine. It is to be noted that the proposed approach is performed only for post fault period. The time variation of the local machine angle (rotor angle of generator 2), obtained from the time-domain simulation in the original test system, is plotted in Fig. 5 (blue curve). The time variation of the equivalent machine angle ��, obtained from the first-stage of the approach, is also illustrated in Fig. 5 (red curve). The variation of the local machine active power is plotted in Fig. 6. The variations of �� and �� are also illustrated in Figs. 7 and 8, respectively. Figs. 5, 7, and 8 clearly show that the values of ��, ��, and �� progressively settle at certain values as the time increases. The final estimates of several system variables including the rotor angle of both the local machine and the equivalent machine, along with some other system variables such as the inertia constant and damping coefficient of the local machine, are listed in Table 1. The values of "� and #� are estimated in the second stage of the approach. Before presenting the simulation results for "� and #�, it is worthwhile to investigate the effect of the inertia constant and damping coefficient of a machine on its dynamic behaviour (or on its rotor angle oscillations). For this purpose, consider the swing equations of a synchronous machine, e.g. Eq. (18). As can be seen in Eq. (18), only for the case in which the machine angle varies with time, its dynamic behaviour will be affected by its inertia constant and damping coefficient. This is due to the fact that when the machine angle remains constant, the first and second time derivatives of its angle will be equal to zero. For such a case, the machine inertia constant and damping coefficient are indeed omitted in the machine swing equations. Consequently, the machine inertia constant and damping coefficient will have no effect on the rotor angle oscillations. Instead, if the machine angle quickly varies with time, the machine dynamic behaviour is greatly affected by its inertia constant and damping coefficient. Table 1 Final estimated values for parameters of the TME system �� (rad) �� (p.u.) �� (p.u.) �� (rad) �� (p.u.) !�� (p.u.) "� (p.u.) #� (p.u.) "� (p.u.) #� (p.u.) 0.7597 0.9709 0.0051 1.1694 1.2255 6.3339 0.1607 0.1607 2.2152 0.9447 From the above discussion, it is concluded that in the beginning of post fault period, rotor angle of a synchronous machine conveys much information regarding its inertia constant and damping coefficient. Therefore, instead of a forgetting factor, a memorizing factor can be used in the RLSE algorithm of the second stage. As mentioned in section 3, a forgetting factor is a number less than 1; whereas a memorizing factor is a number more than 1. If a memorizing factor is used in the RLSE algorithm, then compared with the recent measurements, more weights are given to the older data points. Fig. 5 Variations of the local and equivalent machines angles Fig. 6 Active power variations of the local machine after an 0.05 s fault on the original system 243 Advances in Technology Innovation, vol. 6, no. 4, 2021, pp. 235-250 Fig. 7 Voltage magnitude variations of the equivalent machine Fig. 8 Variations of the transmission line reactance The time variations of "� and #� obtained by utilizing a fixed memorizing factor of 1.002 are illustrated (green curves) in Figs. 9 and 10, respectively. It can be easily seen in Figs. 9 and 10 that, except for the first few seconds, the values of "� and #� remain almost constant in post fault period. The time plots in Figs. 9 and 10 are consistent with the time variations of the local machine and equivalent machine angles shown in Fig. 5. As can be seen in Fig. 5, those angles progressively settle at certain values in post fault period. However, for the first few seconds in post fault period, considerable time variations in both rotor angles can be observed. The final estimated values of "� and #� are also given in Table 1. To show the advantage of memorizing factor, a fixed forgetting factor of 0.99 is also used in the second stage of the proposed approach, and then the time variations of "� and #� are obtained. The time variations obtained for "� and #� are plotted (red curves) in Figs. 11 and 12, respectively. As can be seen, now there exist big fluctuations in the values of "� and #�. As a matter of fact, it is very difficult or even impossible to estimate the values of "� and #�, in this case. This proves the superiority of the memorizing factor to the forgetting factor. The time variations of "� and #� obtained by considering three different values for the memorizing factor are also obtained and plotted in Figs. 9 and 10, respectively. The values chosen for the memorizing factor are 1.004, 1.01, and 1.05. In addition, the time variations of "� and #� obtained by employing LS technique are illustrated in Figs. 9 and 10, respectively. For the case of LS technique, the memorizing factor is assumed to be one (i.e., ��= 1.00). As can be seen in Figs. 9 and 10, depending on the value of the memorizing factor, both "� and #� vary differently over time. Furthermore, for each value of the memorizing factor, both "� and #� progressively settle at some different values in the post-fault period. The final estimated values of "� and #�, for various values of the memorizing factor, are listed in Table 2. Fig. 9 Inertia constant variations of the equivalent machine by using memorizing factor 244 Advances in Technology Innovation, vol. 6, no. 4, 2021, pp. 235-250 Fig. 10 Damping coefficient variations of the equivalent machine by using memorizing factor Fig. 11 Inertia constant variations of the equivalent machine by using forgetting factor Fig. 12 Damping coefficient variations of the equivalent machine by using forgetting factor Table 2 Comparison of parameters of the TME system �� "� (p.u) #� (p.u.) CCT (sec) 1.00 1.6071 1.0257 0.251 1.002 2.2152 0.9447 0.254 1.004 2.4208 0.9575 0.255 1.01 1.6417 2.0033 0.253 1.05 1.7766 2.2474 0.254 In an RLSE based estimation method, a system of overdetermined equations needs to be solved (referring to section 3). For a system of overdetermined equations, there exists no solution that satisfies all the equations exactly. Therefore, RLSE algorithm is used to estimate the unknown parameters in an optimal sense. In LS algorithm, to obtain solutions for the unknown 245 Advances in Technology Innovation, vol. 6, no. 4, 2021, pp. 235-250 parameters, each equation is weighted equally. However, in RLSE algorithm, based on the value chosen for the memorizing factor (or the forgetting factor), each equation is weighted differently. Therefore, for each value chosen for ��, the parameters "� and #� vary differently over time. The system of overdetermined equations obtained from Eq. (19) is solved in the second stage of the proposed approach. As can be seen in Eq. (19), when the first and second time derivatives of the local machine angle become zero, and at the same time, the right hand side of Eq. (19) becomes zero, then no new equation is added to the set of aforementioned overdetermined system of equations. As shown in Fig. 5, after a few seconds in the post-fault period, rotor angle of the equivalent machine remains unchanged. Therefore, after a few seconds, the first and second time derivatives of the equivalent machine angle become zero. In addition, as shown in Fig. 6, the active power of the local machine remains unchanged and reaches to its input mechanical power after a few seconds. Consequently, the right hand side of Eq. (19) becomes zero. In fact, the right hand side of Eq. (19) represents the dependent variable y used in Eq. (4). Moreover, the first and second time derivatives of the equivalent machine angle are indeed the independent variables denoted by the vector x in Eq. (4). Referring to Eq. (4), it can be concluded that when both y and x are equal to zero, the unknown vector of parameters Φ is no longer updated by RLSE algorithm. The vector Φ consists of the mechanical parameters "� and #� of the equivalent machine e. As mentioned in section 3, the greater the memorizing factor is, the less the influence of the recent data points (or equations) on the estimation algorithm is. However, if the memorizing factor is a value very greater than 1, then almost all equations in a system of overdetermined equations are greatly discounted. This means that even those past equations that are very important for estimating the unknowns are neglected. For such a case, the parameters cannot be estimated by RLSE algorithm. The curves shown in Figs. 9 and 10 confirm this fact. As seen, for a memorizing factor of greater than 1.05, RLSE algorithm may fail to correctly estimate the parameters "� and #�. The memorizing factor leading to the best (or optimal) estimation for the unknown parameters will depend on the system of overdetermined equations being considered. As stated before, the values of "� and #� cannot be estimated by using RLSE algorithm with a forgetting factor. For the case of forgetting factor, past samples (equations) are discounted. The less the forgetting factor is, the less the influence of past equations on the estimation algorithm is. The variations of "� and #� obtained by using a fixed forgetting factor of 0.986 are also illustrated (blue curves) in Figs. 10 and 11, respectively. As shown, compared with the case of a fixed forgetting factor of 0.99, there exist higher fluctuations in the values of "� and #� in this case. This is due to the fact that the less the forgetting factor, the less the effect of past equations in estimating the parameters "� and #�. It should be noted that only the values of "� and #� estimated in the second stage of the approach are dependent on the value of the memorizing factor. The parameters ��, ��, and �� are estimated in the first stage of the approach and their values do not depend on the value chosen for the memorizing factor. Therefore, for all cases considered in this sub-section, the values of ��, ��, and �� are the same as those given in Table 1. 7.2. Results for the derived TME system In this sub-section, the derived TME system is used to replicate the behaviour of the original unreduced system by considering a three-phase fault on both systems. A three-phase self-clearing type fault of 0.1 s duration is considered at bus 2 in the original test system. The same fault is also considered at bus 2 in the derived TME system. As it is well known, relative rotor angles rather than absolute rotor angles need to be used in evaluating the stability of a power system. To obtain the relative rotor angles, it is required to take a reference rotor angle. The relative rotor angle in the TME system is easily obtained as: ��� − �� , where �� and �� are the angles of local machine 2 and equivalent machine e, respectively. For the original test system, the inertia weighted average angle (��()) of machines in the external system is used as an approximate reference rotor angle for local machine 2. The angle ��() is defined earlier in Eq. (12). This relative rotor angle is denoted as: ��� − ��() . 246 Advances in Technology Innovation, vol. 6, no. 4, 2021, pp. 235-250 The time variation of the relative angle ��� − ��() is shown in Fig. 13 (solid blue curve). The time variation of the relative angle ��� − �� is also illustrated in Fig. 13 (dotted red color). As can be seen in Fig. 13, the relative rotor angle ��� − �� of the TME system closely approximates the relative rotor angle ��� − ��() of the original system. The time variations of the active and reactive powers of local machine 2 on both systems are shown in Figs. 14 and 15, respectively. It can again be observed in Figs. 14 and 15 that the active and reactive powers of the local machine in the TME system closely approximate the corresponding powers in the original system. To further investigate the effectiveness of the proposed method, a three-phase self-clearing type fault is considered at bus 2 on both the TME system and the original test system. The critical clearing time (CCT) corresponding to this fault on both systems is then obtained. Note that CCT is used as an index for transient stability analysis of power systems [34]. The CCT is found to be 0.254 s for the TME system, whereas the CCT of 0.228 s is obtained for the original system. The values of CCT on both systems are obtained via a trial-and-error approach and performing several time-domain simulations. It is obvious that the CCT obtained from the TME system is near to its value in the original system, and this confirms the correctness of the proposed approach. Fig. 13 Relative rotor angles of local machine after an 0.1 s fault on the original and TME systems Fig. 14 Active power output of local machine after an 0.1 s fault on the original and TME systems Fig. 15 Reactive power output of local machine after an 0.1 s fault on the original and TME systems The CCTs for other values of "� and #� mentioned in Table 2 are also obtained in the TME system. The obtained CCTs are listed in Table 2. As seen in this table, the CCTs obtained by using different values for the memorizing factor are almost equal to each other. This can be justified as follows. On the one hand, it is to be noted that the differences between the values of "� and #� are rather small, for various values of the memorizing factor. On the other hand, it can be noticed that the assumed fault is a severe three-phase short circuit at the terminal of the local machine. For this fault, the rest of the system, which is modelled here by an equivalent machine, has a little effect on the dynamic performance of the local machine. Therefore, the values of "� and #� will only have a little impact on the CCT value as well. 247 Advances in Technology Innovation, vol. 6, no. 4, 2021, pp. 235-250 8. Discussion The simulation results presented in section 7 shows that there exist some differences between the behaviour of the original system and that of the derived TME system. These differences may be associated with the model defined for the TME system. As shown in Fig. 2, the equivalent machine e is assumed to be connected to the local machine i through a lossless transmission line. As a modification, a resistance can be included in series with the transmission line in the defined TME system. It is expected that by this modification, the TME system can approximate the behaviour of the original system with more accuracy. In this case, however, all the necessary equations mentioned in sections 5 and 6 are to be changed appropriately. In this study, the parameters of the equivalent system are obtained by using simulation data from a single fault scenario. In addition, just one system operating condition is considered. In order to improve the performance of the reduced (equivalent) model, multiple faults at different locations with different system operating conditions can be considered for estimating the parameters of the equivalent model. Furthermore, the approach proposed in this study can be used in combination with a heuristic optimization-based method such as particle swarm optimization (PSO) algorithm [35], to enhance the estimation accuracy of the TME system parameters. The time variations of the required system variables are obtained by considering a three-phase fault on the original test system in this study. It is important to point out that the required time variations of system variables can be acquired from an actual system disturbance event or an actual field measurement on a multimachine power system. The first three steps of the proposed approach, shown by the dotted red color in Fig. 3, can be avoided if the measured (recorded) data are utilized in estimating the parameters of the TME system. Furthermore, if some of the required system data are unavailable, they can be estimated by employing a suitable state estimation technique. The use of local measurements available from PMUs is also an alternative way for estimating the required data [24]. Nevertheless, these issues have not been addressed in the present study. The TME system is derived for the classical model of power system; however, it can be used for more accurate transient stability analysis of the local machine as well. To this end, the detailed generator model with excitation system and other control systems can be applied to the local machine in the derived TME system. The TME system can then be employed for more accurate design of power system stabilizer (PSS) [1] or other controllers for the local machine. As described in section 7.1., the idea of memorizing factor is inspired by analyzing the behaviour of synchronous generators, following a three-phase fault on a multimachine power system. However, the memorizing factor based RLSE algorithm can, in general, be used for solving any set of over-determined algebraic equations in which the older equations are more important than the recent ones. As a matter of fact, the dynamic equivalencing approach proposed in this study can be considered one of the applications of the memorizing factor concept. As stated in section 7.1., for a given system of overdetermined equations, there exists an optimal value for the memorizing factor that leads to the best estimation of the unknown parameters. Future research can be focused on finding such an optimal value for the memorizing factor. 9. Conclusions In this study, a systematic two-stage identification-based approach was presented to tackle the challenging problem of multimachine power systems dynamic equivalencing. In particular, a TME system was derived for a multimachine power system, by looking into the system from the bus of a local system machine. First, the TME system was defined, and then the desired parameters of the predefined equivalent system were estimated by using a two-stage approach. The electrical parameters of the equivalent system were obtained by utilizing a simple technique in the first stage of the approach. The mechanical parameters of the equivalent system were obtained in the second stage of the approach by employing an innovative RLSE technique with a “memorizing factor”. The effectiveness of the proposed methodology was validated by using the New England 39-bus test system. 248 Advances in Technology Innovation, vol. 6, no. 4, 2021, pp. 235-250 Many simulation results were presented to show the correctness of the proposed approach. The advantages of using a memorizing factor rather than a forgetting factor were justified and also verified by many simulation results. Instead of the traditional SMIB system, the derived TME system can be employed for transient stability analysis as well as other dynamical studies in multimachine power systems. It can also be used for more accurate design of PSS or other controllers for the local machine. The memorizing factor based RLSE algorithm introduced in this study is completely new, and can be applied in any system identification-based problem in which the older samples are more important than the newly available data points. The concept of varying memorizing factor can be regarded as a future work. Conflicts of Interest The author declares no conflict of interest. References [1] J. Machowski, Z. Lubosny, J. Bialek, and J. R. Bumby, Power System Dynamics: Stability and Control, 3rd ed. Hoboken: John Wiley & Sons, 2020. [2] A. Karami and K. M. Galougahi, “Improvement in Power System Transient Stability by Using STATCOM and Neural Networks,” Electrical Engineering, vol. 101, no. 1, pp. 19-33, April 2019. [3] E. C. Machado and J. E. Pessanha, “Foundations on the Hamiltonian Energy Balance Method for Power System Transient Stability Analysis: Theory and Simulation,” Journal of Control, Automation, and Electrical Systems, vol. 31, no. 1, pp. 226-232, February 2020. [4] Y. F. Gao, J. Q. Wang, T. N. Xiao, and D. Z. Jiang, “Fast Emergency Control Strategy Calculation Based on Dynamic Equivalence and Integral Sensitivity,” Frontiers of Information Technology and Electronic Engineering, vol. 20, no. 8, pp. 1119-1132, August 2019. [5] X. Y. Li and O. P. Malik, “Estimation of Equivalent Models for Emergency State Control of Interconnected Power-Systems Based on Multistage Recursive Least-Squares Identification,” IEE Proceedings C (Generation, Transmission, and Distribution), vol. 140, no. 4, pp. 319-325, July 1993. [6] J. B. Ward, “Equivalent Circuits for Power-Flow Studies,” Transactions of the American Institute of Electrical Engineers, vol. 68, no. 1, pp. 373-382, July 1949. [7] P. Dimo, Nodal Analysis of Power Systems, Lancashire: Abacus Press, 1975. [8] J. Srivani and K. S. Swarup, “Power System Static Security Assessment and Evaluation Using External System Equivalents,” International Journal of Electrical Power and Energy Systems, vol. 30, no. 2, pp. 83-92, February 2008. [9] Y. Fu and T. S. Chung, “A Hybrid Artificial Neural Network (ANN) and Ward Equivalent Approach for On-Line Power System Voltage Security Assessment,” Electric Power System Research, vol. 53, no. 3, pp. 165-171, March 2000. [10] B. G. Lee, J. Lee, and S. Kim, “Development of a Static Equivalent Model for Korean Power Systems Using Power Transfer Distribution Factor-Based k-Means++ Algorithm,” Energies, vol. 13, no. 24, 6663, December 2020. [11] U. D. Annakkage, N. K. C. Nair, Y. Liang, A. M. Gole, V. Dinavahi, B. Gustavsen, et al., “Dynamic System Equivalents: A Survey of Available Techniques,” IEEE Transactions on Power Delivery, vol. 27, no. 1, pp. 411-420, January 2012. [12] J. H. Chow, Power System Coherency and Model Reduction, Springer Science and Business Media, 2013. [13] M. R. A. Paternina, J. M. Ramirez-Arredondo, J. D. Lara-Jiménez, and A. Zamora-Mendez, “Dynamic Equivalents by Modal Decomposition of Tie-Line Active Power Flows,” IEEE Transactions on Power Systems, vol. 32, no. 2, pp. 1304-1314, March 2017. [14] C. M. Rergis, R. J. Betancourt, and A. R. Messina, “Order Reduction of Power Systems by Modal Truncated Balanced Realization,” Electric Power Components and Systems, vol. 45, no. 2, pp. 147-158, 2017. [15] D. Chaniotis and M. A. Pai, “Model Reduction in Power Systems Using Krylov Subspace Methods,” IEEE Transactions on Power Systems, vol. 20, no. 2, pp. 888-894, May 2005. [16] Z. Zhu, G. Geng, and Q. Jiang, “Power System Dynamic Model Reduction Based on Extended Krylov Subspace Method,” IEEE Transactions on Power Systems, vol. 31, no. 6, pp. 4483-4494, November 2016. [17] T. K. Ritschel, F. Weiß, M. Baumann, and S. Grundel, “Nonlinear Model Reduction of Dynamical Power Grid Models Using Quadratization and Balanced Truncation,” Automatisierungstechnik, vol. 68, no. 12, pp. 1022-1034, November 2020. 249 Advances in Technology Innovation, vol. 6, no. 4, 2021, pp. 235-250 [18] R. Podmore, “Identification of Coherent Generators for Dynamic Equivalents,” IEEE Transactions on Power Apparatus and Systems, vol. 97, no. 4, pp. 1344-1354, July 1978. [19] J. J. Ayon, E. Barocio, I. R. Cabrera, and R. Betancourt, “Identification of Coherent Areas Using a Power Spectral Density Algorithm,” Electrical Engineering, vol. 100, no. 2, pp. 1009-1019, June 2018. [20] A. M. Miah, “Study of a Coherency-Based Simple Dynamic Equivalent for Transient Stability Assessment,” IET Generation, Transmission, and Distribution, vol. 5, no. 4, pp. 405-416, April 2011. [21] M. Matar, N. Fernandopulle, and A. Maria, “Dynamic Model Reduction of Large Power Systems Based on Coherency Aggregation Techniques and Black-Box Optimization,” International Conference on Power Systems Transients, June 2013, pp. 1-7. [22] S. Chittora and S. N. Singh, “Coherency Based Dynamic Equivalencing of Electric Power System,” 18th National Power Systems Conference, December 2014, pp. 1-6. [23] A. M. Stankovic, A. T. Saric, and M. Milosevic, “Identification of Nonparametric Dynamic Power System Equivalents with Artificial Neural Networks,” IEEE Transactions on Power Systems, vol. 18, no. 4, pp. 1478-1486, November 2003. [24] A. Chakrabortty and A. Salazar, “Building a Dynamic Electro-Mechanical Model for the Pacific AC Intertie Using Distributed Synchrophasor Measurements,” European Transactions on Electrical Power, vol. 21, no. 4, pp. 1657-1672, May 2011. [25] J. M. Ramirez, B. V. Hernández, and R. E. Correa, “Dynamic Equivalence by an Optimal Strategy,” Electric Power System Research, vol. 84, no. 1, pp. 58-64, March 2012. [26] J. L. Rueda, J. Cepeda, I. Erlich, D. Echeverría, and G. Argüello, “Heuristic Optimization Based Approach for Identification of Power System Dynamic Equivalents,” International Journal of Electrical Power and Energy Systems, vol. 64, pp. 185-193, January 2015. [27] X. Zhang, Y. Xue, S. You, Y. Liu, Z. Yuan, J. Chai, et al., “Measurement-Based Power System Dynamic Model Reductions,” North American Power Symposium, September 2017, pp. 1-6. [28] N. Tong, Z. Jiang, S. You, L. Zhu, X. Deng, Y. Xue, et al., “Dynamic Equivalence of Large-Scale Power Systems Based on Boundary Measurements,” American Control Conference, July 2020, pp. 3164-3169. [29] Z. Jiang, N. Tong, Y. Liu, Y. Xue, and A. G. Tarditi, “Enhanced Dynamic Equivalent Identification Method of Large-Scale Power Systems Using Multiple Events,” Electric Power Systems Research, vol. 189, 106569, December 2020. [30] A. B. Almeida, R. Reginatto, and R. J. G. C. D. Silva, “A Software Tool for the Determination of Dynamic Equivalents of Power Systems,” IREP Symposium Bulk Power System Dynamics and Control, August 2010, pp. 1-10. [31] K. J. Keesman, System Identification: An Introduction, London: Springer, 2011. [32] M. A. Pai, Energy Function Analysis for Power System Stability, Boston: Kluwer Academic Publishers, 1989. [33] R. D. Zimmerman, C. E. Murillo-Sánchez, and R. J. Thomas, “MATPOWER: Steady-State Operations, Planning, and Analysis Tools for Power Systems Research and Education,” IEEE Transactions on Power Systems, vol. 26, no. 1, pp. 12-19, February 2011. [34] A. Karami and S. Z. Esmaili, “Transient Stability Assessment of Power Systems Described with Detailed Models Using Neural Networks,” International Journal of Electrical Power and Energy Systems, vol. 45, no. 1, pp. 279-292, February 2013. [35] A. Slowik, Swarm Intelligence Algorithms: A Tutorial, Boca Raton: CRC Press, 2021. Copyright© by the authors. Licensee TAETI, Taiwan. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY-NC) license (https://creativecommons.org/licenses/by-nc/4.0/). 250