PMU-based voltage stability prediction using least square support vector machine with online learning

Electric Power Systems Research 160 (2018) 234–242

Contents lists available at ScienceDirect

Electric Power Systems Research journal homepage: www.elsevier.com/locate/epsr

PMU-based voltage stability prediction using least square support vector machine with online learning Hao Yang a , Wen Zhang a,∗ , Jian Chen a , Liang Wang b a b

Key Laboratory of Power System Intelligent Dispatch and Control of the Ministry of Education (Shandong University), Jingshi Road, Jinan 250061, China State Grid Shandong Electric Power Company, Jinger Road, Jinan 250002, China

a r t i c l e

i n f o

Article history: Received 23 August 2017 Received in revised form 4 February 2018 Accepted 26 February 2018 Keywords: Short-term voltage stability Composite load model PMU measurements Trajectory extrapolation Unstable equilibrium point

a b s t r a c t Composite load model comprising an impedance-current-power (ZIP) load and an induction motor (IM) is widely used for dynamic voltage stability analysis. With the load model, this paper proposes a new PMU-based method to predict short-term voltage instability (STVIS). Using synchronous measurements, the method firstly performs online contingency analysis to obtain a set of three-element look-up tables comprising the presumed contingency, the corresponding post-fault stable and unstable equilibrium point (SEP and UEP) for each IM load individually. Next, when a fault really occurs, the time-series of IM slip is computed by an Euler algorithm based on local PMU measurements, and then a new time-series prediction method is proposed for rolling prediction of IM slip trajectory by introducing least square support vector machine (LSSVM) with online learning. Finally, from the view of IM stability mechanism, the STVIS status can be detected in advance by monitoring that the predicted slip trajectory reaches the IM’s UEP in the look-up table. The effectiveness of the proposed method is verified on the New England 39-bus system. © 2018 Elsevier B.V. All rights reserved.

1. Introduction Short-term voltage stability (STVS) refers to the ability of power system to keep steady voltages at all buses after a large disturbance in a transient time frame [1]. Unlike the rotor angle stability focusing on generators’ synchronism, STVS closely relates to the stability of dynamic loads. The short-term voltage instability (STVIS) is mainly caused by the IM dynamics tending to restore the power consumption beyond the capability provided by the post-fault system, resulting in IM stalling. During the stalling process, much reactive power is absorbed, deteriorating the voltage levels. Further, this in turn can make more IMs installing, and consequently, STVIS occurs. In modern power system, as the penetration of IMs (e.g. air conditioners and industrial motors) increases, it is crucial to predict the STVIS status and timely activate emergency control measures. For the steady-state voltage stability analysis, there are some mature methods, such as continuation method [2], singular value decomposition [3], and Thevenin equivalent based method [4]. However, these approaches are derived from power flow equations

∗ Corresponding author at: Electrical Engineer School, Shandong University, Jingshi Road 17923, Jinan 250061, China. E-mail address: [email protected] (W. Zhang). https://doi.org/10.1016/j.epsr.2018.02.018 0378-7796/© 2018 Elsevier B.V. All rights reserved.

based on the static models of system, and hence cannot be applied for STVS analysis. Many studies on STVS analysis focus on analyzing the IM dynamics after a large disturbance [5–8]. In Ref. [9], by using the transient P–V curves, the load characteristics of IM in the P–V plane is introduced and then an analytical method for STVS is put forward in a single-load infinite-bus system. In Ref. [10], to prevent IM from losing stability, the critical clearance time (CCT) of voltage sag is analytically derived, also in a single-load infinite-bus system. Due to the same time frame of STVS and rotor angle stability, the impacts of generator out-of-step on IM stability are studied in Ref. [11]. Inspired by the energy function method on rotor angle stability analysis, Refs. [12,13] used it to analyze voltage stability. However, the method may have difficulties for STVS analysis when applied to real systems with many IM loads. So far, time domain simulation (TDS) is still an effective way for large-bulk systems. Although the above methods play important roles in addressing STVS problems, they are still difficult for fast real-time detection of STVIS. In practice, the STVIS is detected by industrial experiential criteria in which only the duration of voltage sag of any bus under a predefined threshold exceeds a time-period, the STVIS status is detected [14]. But, the experiential criteria may lack reliability and explanation for different systems. Given the increasing placement of phasor measurement units (PMUs) in system, some machine-learning techniques have been adopted for real-time stability assessment, such as neural network

H. Yang et al. / Electric Power Systems Research 160 (2018) 234–242

(NN), support vector machine (SVM), decision tree (DT), time series shapelet, and intelligent system [15–21]. Based on these techniques, the stability assessment is treated as a binary classification problem, and it mainly includes two stages: offline training and online application. The classifier is obtained after offline training stage where a large stability/instability database are needed. At the online stage, with post-fault dynamic features (e.g. voltage trajectories) inputting into the classifier, the stability status can be detected early. These methods rely on massive TDSs requiring detailed system models and their accurate parameters. However, for complex system, the components’ parameters and models may be uncertain, and the generated stability/instability database based on TDS may be insufficient for training stage. Both may result in an unreliable classifier, and give unreliable conclusions. Further, these methods focus on the data themselves rather than the intrinsic mechanism of STVS, which may also lead to an unreliable result. Except the above machine-learning methods, Lyapunov exponent method is an effective tool to diagnose the system chaotic behavior [22]. In Ref. [22], by measuring the post-fault time series of voltage with PMUs, the maximum Lyapunov exponent is computed and then the STVIS status can be detected by the Lyapunov exponent. Currently, for STVS assessment, the key problem is that the appropriate dynamic load model should be adopted. For engineering and research experience, the composite load model is increasingly recommended and widely adopted in network operators for dynamic stability analysis, especially STVS analysis [23,24]. And, thanks to the wide use of PMUs, the accuracy of parameter identification of composite load is improved greatly by using measurement-based methods [25–28]. In this paper, with the composite load, we propose a PMU-based method to predict STVIS for multi-bus power systems. The method is divided into two stages. The first stage use synchronous measurements in operation center to perform online presumed contingency analysis to create a look-up table. In second stage, with local measurement, the post-fault STVIS status is detected in advance based on time-series prediction and the look-up table. The main contributions are as follows. 1) The method is a semi-parameter-free approach where only the power flow data and the parameters of composite load are required, which greatly reduces the dependence on the dynamic modeling and parameters. 2) By introducing LSSVM with online sequential learning, a new time-series prediction method is proposed to predict the IM slip trajectory in a rolling way. 3) With the predicted slip trajectory, the STVIS status can be detected in advance based on the instability mechanism of IM, rather than the experience criteria. 4) The method has applicability under symmetrical/asymmetrical fault conditions. The rest of this paper is organized as follows. In Section 2, the composite load model is briefly reviewed, and then the method to calculate the post-fault SEP and UEP of IM is presented. The STVS assessment based on the SEP and UEP is also introduced here. Section 3 provides a new time-series prediction method based on online learning LSSVM, and the scheme of PMU-based STVIS prediction method is presented. Simulation results are illustrated in Section 4, followed by conclusions in Section 5.

235

Fig. 1. (a) Composite load structure. (b) Circuit model of composite load.

The static ZIP load is formulated as follows:

 

2

 

2

Pzip = P0 pz V/V0

Qzip = Q0 qz V/V0













+ pi V/V0 + pp

+ qi V/V0 + qp

(1)

where V0 and V are the rated voltage magnitude and operating voltage magnitude, respectively; P0 and Q0 are real and reactive power respectively, before disturbance; both parts of (pz , pi , and pp ) and (qz , qi , and qp ) are the coefficients for the real and reactive power, respectively, which the sum of each part is 1. For IM, its equations are described as: dE˙  1 ˙ = −jsE˙  −  [E˙  − j(X − X  )I] T dt

(3)

1 ds = (Tm − Te ) 2Tj dt

(4)

V˙ = E˙  + (rs + jX  )I˙



 ∗





2.2. Look-up table of SEP and UEP of composite load The SEPs and UEPs of post-fault system play important roles in analyzing the stability based on the stability theory of dynamic system. For the STVS analysis, to obtain the post-fault IM’s SEP and UEP is of great concern. Here, the method for solving the post-fault SEP and UEP is as follows. In Eqs. (3) and (4), by setting the differential term to 0, we get:



1  ˙ [E˙ − j(X − X  )I] T

(6)





0 = T0 A(1 − s)2 + B (1 − s) + C − Re E˙  I˙

2.1. Composite load model The structure and circuit diagram of the composite load are shown in Fig. 1 where rs is the stator resistance, xs is the stator reactance, rr is the rotor resistance, xr is the rotor reactance, xm is magnetizing reactance and s is the IM slip.

(5)

where Te = Re E˙  I˙ , Tm = T0 A(1 − s)2 + B (1 − s) + C , T  = (xr + xm ) /rr , X = xs + xm , X  = xs + xr xm / (xr + xm ) .In the above equations, E˙  refers to the transient EMF of IM; I˙ is the stator current phasor; Tj is the inertia constant; V˙ is the voltage phasor. Tm and Te are the mechanical and electrical torque, respectively; T0 is the initial mechanical torque and A, B and C are the torque coefficients that has A + B + C = 1. Usually, the values of B and C are zeros. Re denotes the function to obtain the real part of phasor. Based on Eqs. (1)–(5), the model of composite load can be formulated by a set of differential-algebraic equations (DAEs), where the STVS can be evaluated by analyzing the dynamics of the composite load. For more details about the composite load, please refer to Refs. [25,26]. The measurement-based methods for estimating the parameters of composite load can be seen in Refs. [25–28]. In this study, we mainly focus on the application of the composite load for STVS analysis.

0 = jsE˙  +

2. Short-term voltage stability assessment

(2)

 ∗

.

(7)

To solve SEP and UEP of IM slip in Eq. (7), the expressions of E˙  and I˙ are first obtained by combining Eqs. (5) and (6): E˙  = V˙ I˙ = V˙

j(X − X  ) (rs − sT  X  ) + j(sT  rs + X)





j(X − X  ) 1 1−   rs + jX (rs − sT X  ) + j(sT  rs + X)

(8)

 .

(9)

236

H. Yang et al. / Electric Power Systems Research 160 (2018) 234–242

Then, substituting Eqs. (8) and (9) into Eq. (7), one has:



T0 [A(1 − s)2 + B(1 − s) + C] = V 2 Re



j(X − X  ) 1 ×  (rs − sT X  ) + j(sT  rs + X) rs +jX 

where V is the post-fault steady-state voltage magnitude; * denotes the conjugation of the phasor. By solving Eq. (10), the SEP and UEP of IM slip are obtained. To obtain the SEP and UEP using Eq. (10), the post-fault steady voltage magnitude of load bus is needed. Therefore, the power flow of post-fault system should be solved ahead. Considering that when a fault occurs, it will spend several seconds calculating the post-fault power flow for a bulk system, there is not enough time for real-time STVS analysis owing to the time frame (say several seconds) of STVS problem. Based on this, a look-up table is built beforehand by online presumed contingency analysis. This table is a three-element set consisting of the presumed contingency, the corresponding post-fault SEP and UPE for each composite load. By synchronous measurements, the online presumed contingency analysis can be performed considering the real-time operation of system, and hence the look-up table is obtained and updated online. By the look-up table, the SEP and UEP of each composite load can be sought out timely for post-fault STVS prediction when a contingency has really occurred. 2.3. STVS assessment based on SEP and UEP STVS can be assessed by analyzing whether the post-fault trajectory of IM slip will converge to its post-fault SEP or diverge through the UEP when starting from the initial state after fault-clearance. Typical slip trajectories and the electrical-mechanical torque versus rotor slip are shown in Fig. 2. In the figure, Apre , Apost and Bpost are the pre-fault SEP, post-fault SEP and UEP corresponding to the IM slips spre-sep , ssep and suep , respectively. tc1 , tc2 and tc3 are the different fault clearing times with tc1 < tc2 < tc3 . L1 , L2 and L3 are three loci of IM slip with fault clearing times tc1 , tc2 and tc3 , respectively. TE -Pre and TE -Post are steady electrical torque for the pre-fault and postfault system, respectively. TE -Sag is electrical torque during fault. TM is the mechanical torque. The STVS status depends on the fault clearing time and the IM net torque (electrical torque minus mechanical torque). Locus L1 denotes a stable case in which the fault is cleared at an earlier time tc1 . The slip trajectory is attracted by the post-fault SEP ssep at point Apost . Owing to voltage restoration in the post-fault period, the net torque will become positive and the IM slip will converge to ssep . Locus L2 is for an unstable case where the fault is cleared at time tc2 . The variations of L2 are determined by the net torque affected by

Fig. 2. Electrical-mechanical torque versus rotor slip curves.



1−

j (X − X  )  (rs − sT X  ) + j (sT  rs + X)

∗ (10)

the dynamics of the whole system. Once the IM slip passes through suep , it will monotonically increase to be close to standstill, and the STVIS occurs. That is because the net torque will be negative after the slip passes through suep even if the voltage restores to its postfault steady value. Locus L3 is another unstable case. At fault clearing time tc3 , the IM slip has crossed suep . Then L3 shows a monotonically increasing property since the net torque is always negative. Based on the above analysis, once the IM slip crosses the postfault suep , it will keep on increasing and IM stalls, and the STVS loses stability. Hence, the assessment method of STVS for multi-bus power system can be described as: 1) when the IM slips of all composite loads restore to their postfault SEPs, the post-fault system is STVS. 2) when the IM slip of one composite load reaches the post-fault UEP, the post-fault system is STVIS. For real-time STVIS prediction, the key problem is to predict whether the IM slip trajectory will pass through the post-fault UEP or not. To obtained the future IM slip trajectory, a fast trajectory extrapolation method based on LSSVM with online learning is presented next section. 3. PMU-based short-term voltage instability prediction In this section, the PMU-based STVIS prediction method is presented. The unwanted random noise and errors from PMU measurements need to be filtered out to obtain ‘clean’ measurements by some filtering and smoothing techniques [29,30] before performing the process. These filtering techniques are not addressed here. Instead, we focus on the STVIS prediction method using the ‘clean’ measurements. The sketch of the proposed STVIS prediction method is given in Fig. 3. As shown in the figure, the procedure in the left-hand dotted box is applied to create a look-up table that is introduced in Section 2.2. To obtain the look-up table, the online synchronous measurements and the parameters of composite load are required first, and then the contingency analysis is performed to get the SEP and UEP for each IM load using Eq. (10). The procedure in the left is carried

Fig. 3. Sketch of the proposed methodology.

H. Yang et al. / Electric Power Systems Research 160 (2018) 234–242

out in the operation center continuously, say for every 5-min. After the look-up table is created, it is transmitted to each load station. In the right box, the procedure is performed in each load station when a fault really occurs. For the right procedure, we first compute the real-time time-series of IM slip using a forward Euler algorithm based on local measurements, and then a new time-series prediction of IM slip by using sequential learning LSSVM is proposed. The last module in Fig. 3 is to detect STVIS based on the predicted IM slip trajectory and the look-up table in each load station individually. The following subsections give the important steps of the timeseries computation method, the time-series prediction method and the specific flowchart of the methodology. 3.1. PMU-based real-time IM slip computation 3.1.1. Slip computation considering symmetrical fault The slip is an internal variable that cannot be measured directly by PMU. To compute slip, first, Te is calculated by, Te = Ptol − Pzip −

(Ptol − Pzip )2 + (Qtol − Qzip )2 V2

(11)

rs

where Ptol , Qtol are respective the bus real and reactive powers, and the third part in the right-hand side is the stator loss of IM. Then, substituting Eq. (11) into Eq. (4), the slip s(tn+1 ) at sample time tn+1 can be calculated using a forward Euler method as: s(tn+1 ) = s(tn ) +



+



1 2 × A(1 − s(tn )) + B (1 − s(tn )) + C − Ptol (tn+1 ) + Pzip (tn+1 ) 2Tj

2 

Ptol (tn+1 ) − Pzip (tn+1 )

2 

V 2 (tn+1 )

rs

Although the computation speed of LSSVM is improved, LSSVM is still a batch learning algorithm. For the batch learning algorithm, whenever a new data is arrived, it will use the past data together with the new arriving data to do a retraining. In every retraining phase, the LSSVM needs to recalculate the inversion of a larger matrix, which is time consuming. To overcome the deficiency, an online learning LSSVM method with recursive process is proposed for prediction of IM slip trajectory, which is given as follows. 3.2.1. Reconstruction of time-series data Based on Takens theorem [33], any scalar time series can be reconstructed to be a phase space where more useful information about the original dynamic system can be explored. Suppose that we have a time series {x1 , x2 , . . ., xn }, which here is the time series of IM slip, the phase space can be reconstructed as follows: xi = (xi , xi−1, . . ., xi−(d−1) ) i = n, n − 1, . . ., 1 + (d − 1)

t

where tn+1 , tn are two consecutive sample times. t = tn+1 − tn . The initial s(t0 ) is spre-sep at the pre-fault steady condition. Ptol (tn+1 ), Qtol (tn+1 ) and V(tn+1 ) are measured by PMU. With measurements, Pzip (tn+1 ), Qzip (tn+1 ) can be derived by Eqs. (1) and (2). By using Eq. (12), the time series of IM slip is obtained online. 3.1.2. Slip computation considering asymmetrical faults The symmetrical component method is used to analyze asymmetrical fault condition. As for the stability analysis, the negativeand zero-sequence components throughout the system are always not of interest. That is because the average electrical torque acting on rotor shaft produced by negative sequence components is approximately zero. And since the zero sequence components do not couple with rotor windings, the electrical torque generated by them is also zero. In stability studies, we pay most attention to the quantities of the positive sequence network. But note that the quantities of positive sequence currents and voltages are affected by the negative- and zero-sequence components. Therefore, under asymmetrical fault case, the IM slip computation only needs to consider the positive-sequence components. Based on the per-phase voltage and current phasors with PMU measurements, the relevant positive- sequence components can be derived by using symmetrical component method. Then, by substituting the positive-sequence components into Eq. (12), the time series of IM slip considering asymmetrical faults can be obtained.

⎡

xd+1

⎛⎡

⎤

⎢x ⎥ ⎜⎢ ⎢ d+2 ⎥ ⎜⎢ ⎢ ⎥ = f ⎜⎢ ⎢ . ⎥ ⎜⎢ ⎣ .. ⎦ ⎝⎣

x1

x2

···

x2

x3

···

.. .

.. .

..

xn−d

xn−d+1

···

.

xd

⎤⎞

⎥⎟ ⎥⎟ · .. ⎥⎟ . ⎦⎠

xd+1 ⎥⎟

(14)

xn−1

Using time-series reconstruction, the time-series is converted to a preferable input–output data form, and then by identifying the mapping function f, the future time-series trajectory in the prediction window can be extrapolated. Owing to f is not straightforward, here it is estimated by the proposed online learning LSSVM. In Eq. (14), the embedding dimension d should be determined first. In this study, d is determined by using CAO algorithm [34], which does not contain any subjective parameters and not strongly depends on the number of available data points. 3.2.2. LSSVM with online learning  m For m training data points (xi , yi ) |xi ∈ Rd , yi ∈ R , x is the i=1 i Here, x is the i-th reconinput vector and yi is the output one. i  structed vector xi , xi−1, · · ·, xi−(d−1) , and yi equals to xi+1 . The LSSVM regression can be formulated as an optimization problem of the following function:

⎧ ⎪ ⎨

1 T C! 2 ei w w+ 2 2 m

min J(w, e) =

w,b,e

⎪ ⎩ s.t.y = wT ϕ(x ) + b + e i=1 m i i i i = 1, · · ·, m

(15)

where ϕ(·) is a nonlinear function mapping the input space into a higher dimensional space; w is a weight vector; bm is the bias term; ei is the error variable; C is regularization parameter. Equivalently, Eq. (15) can be reformulated using the Lagrange function as: L(w, bm , e, ˛) = J(w, e) −

m !

˛i {wT ϕ(xi ) + bm + ei − yi }

(16)

i=1

3.2. Trajectory extrapolation with online learning LSSVM SVM is a tool applying the structural risk minimization principle for solving the problems characterized by small samples, nonlinearity, high dimension and local minima [31]. The solution procedure of the standard SVM needs to solve a quadratic program (QP) with inequality constraints, which is time-consuming. Because of this, a reformulation to the standard SVM, called LSSVM, is introduced in Ref. [32]. The LSSVM is obtained by solving a set of linear equations.

(13)

where d is the embedding dimension, xi is the remodeled vector. Note xi is a time-series data with d dimension. As with Takens theorem, there has a smooth map f: Rd → R that,

xn (12)

+ Qtol (tn+1 ) − Qzip (tn+1 )

237

where ˛i is the Lagrange multiplier. Based on the optimality conditions that the partial derivatives of L with respect to w, bm , e, ˛ are all equal to zero, and further, the following linear equation is given after reduction.

"

0

1Te

1e

Km + C −1 I

"

bm ˛m



# =

0 ym

$ (17)

238

H. Yang et al. / Electric Power Systems Research 160 (2018) 234–242

where 1e = [1,1,. . .,1]T with m dimension; I is the identity matrix; ˛m = [˛1 , ˛2 ,. . ., ˛m ]T ; ym = [y1 , y2 ,. . ., ym ]T ; Km is the kernel matrix and its elements satisfy Km (i, j) = k(xi , xj ), i, j = 1,2,. . .,m. k(xi , xj ) denotes the kernel function, and the typical kernel functions are linear kernel, polynomial kernel, and radial basis function (RBF) kernel. Here, the RBF is adopted. RBF is expressed by k(xi , xj ) = exp{−xj − xi 22 / 2 } where  is a constant parameter. Assuming Hm = Km + C−1 I, ˛m and bm in Eq. (17) are solved as:

⎧ 1Te H−1 ⎪ m ym ⎪ ⎨ bm = T −1 1e Hm 1e

(18)

⎪ 1e 1Te H−1 m ym ⎪ ] ⎩ ˛m = H−1 m [ym − T −1

where Gm+1 = [k(x1 , xm+1 ), . . ., k(xm , xm+1 )]T ; gm+1 = k(xm+1 , xm+1 ) + 1/C. Based on the block matrix inversion formula [35], the inverse of Hm+1 can be written as:

"

H−1 m+1

=

H−1 m

0

0

0



T

+

T −1 [GTm+1 H−1 m , −1] [Gm+1 Hm , −1]

gm+1 − GTm+1 H−1 m Gm+1

.

(24)

can be online calculated using the recursive forTherefore, H −1 m+1

mula in Eq. (24). Then, substituting H −1 and ym+1 into Eq. (18), m+1 ˛m+1 and bm+1 are obtained, and then f can be updated fast by putting ˛m+1 and bm+1 into Eq. (19). With the updated f, the future trajectory is predicted by Eq. (20).

1e Hm 1e

⎡

⎤

k(x1 , x1 ) + 1/C · · · k(xm , x1 ) . .. . ⎦. .. .. where Hm = ⎣ . k(x1 , xm ) · · · k(xm , xm ) + 1/C After obtaining ˛m and bm , the LSSVM-based regression function f can be given as: y = f (x) =

m !

˛i k(xi , x) + bm .

Using the mapping f, the trajectory extrapolation process with l prediction window is given by,

⎧ d ⎪ ⎪ ⎪ % &' ( ⎪ ⎪ xˆ n+1 = f (xn−d+1 , xn−d+2 , · · ·, xn ) ⎪ ⎪ ⎪ ⎪ ⎪ d ⎪ ⎪ ⎪ % &' ( ⎨ xˆ n+2 = f (xn−d+2 , · · ·, xn , xˆ n+1 )

(20)

.. . d

%

&'

(

xˆ n+l = f (xn−d+l , · · ·, xn , xˆ n+1 , · · ·, xˆ n+l−1 )

where x and xˆ are the observed value and predicted value. By using Eq. (20), the future time-series trajectory with l time steps can be predicted. The main computing burden for obtaining f is the inversion of Hm in Eq. (18). For online time series prediction, much data is collected online and the training set is increased with time. To updating the f, a fast-recursive algorithm is described as follows. When the (m + 1)-th time data arrives, Eq. (17) is rewritten as:

"

"

0

1Te

1e

Km+1 + C −1 I

bm+1



˛m+1

# =

$

0

(21)

ym+1

where ˛m+1 = [˛1 , ˛2 ,. . ., ˛m , ˛m+1 ]T ; ym+1 = [y1 , y2 ,. . ., ym , ym+1 ]T ; Hm+1 = Km+1 + C−1 I has:

⎡

k(xm+1 , x1 )

⎢ ⎢ ⎢ ⎣

.. .

Hm

k(x1 , xm+1 )

···

k(xm , xm+1 )

k(xm+1 , xm ) k(xm+1 , xm+1 ) + 1/C

⎤ ⎥ ⎥ ⎥ ⎦

(22)

Rewriting Eq. (22), we get:

" Hm+1 =

Hm

Gm+1

GTm+1

gm+1

As mentioned in Section 2.3, the criterion for STVIS detection can be derived from the IM dynamics. Hence, the criterion is written as follows: j

j

Istvs =

suep − sˆj (t) j

j

suep − ssep

≤0

j ∈ ˛cm

(25)

(19)

i=1

⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩

3.3. STVIS detection criterion

 (23)

j

j

where ˛cm is the set of composite loads; ssep and suep are the post-fault SEP and UEP of the j-th IM; sˆj (t) is the predicted IM slip j based on the proposed extrapolation method; Istvs is the index of j

STVS. For any composite load, if Istvs ≤ 0, the system is declared as status of STVIS. j Based on Istvs index, the weak-bus of STVS is defined as follows. For transient stable case, all IMs keep stable after the transient j time-scale and the weak-bus is the load bus whose Istvs has a minimum value during the transient time-scale (the IM has a maximum fluctuation of rotor slip). For unstable case, the weak-bus refers j to the load bus whose Istvs first satisfies the criterion in Eq. (23). Note in this paper, the concept of weak-bus based on IM dynamics differs from the existing weak-bus identification based on static approaches [36,37] using steady model. The static method based weak-bus refers to the load bus has the minimum load margin at one operating state. Fig. 4 gives the specific flowchart of the STVIS prediction method. 4. Case studies The New England 39-bus system [38] shown in Fig. 4 is used to demonstrate the validation of the proposed method. The dynamic simulation was simulated in the commercial software DIgSILENT/Power factory [39] and the real-time STVIS prediction was realized in the Matlab platform. The DIgSILENT and Matlab were in a real-time communication environment, which were executed on a 32-bit PC with 3.4-GHz CPU and 4-GB RAM. The sampling rate for the phasor output streaming from a commercially-available PMU is typically in the range of 0.5–2 samples per cycle [22]. That is the sampling time is in the range of 0.01 s–0.04 s for 50-HZ system. Here, the sampling time is set to be 0.01 s from the TDS by DIgSILENT. Refs. [19] and [40] also adopt the 0.01 s sampling time. 4.1. Online learning LSSVM settings For the LSSVM application, the first task is to determine the tuning parameters (regularization parameter C and the constant parameter  in RBF kernel). Here, the initial training data is collected by using the post-fault sampling data of 300 ms time period (e.g. the initial training time window is 300 ms), and the 4-fold cross validation is used to determine the LSSVM parameters (C and ). And once the parameters C and  are optimized, they keep

H. Yang et al. / Electric Power Systems Research 160 (2018) 234–242

239

Fig. 5. New England 39-bus test system. Table 1 Post-fault equilibrium points for tripping line 4–14.

Fig. 4. Flowchart of the proposed PMU-based STVIS prediction method.

unchanged, and then the online rolling prediction of IM slip trajectory is implemented. The initial prediction time window is set to 100 ms. With new incoming data, the training data is increased, and in this study, the ratio between the training time window and prediction time window is kept as 3:1 until the prediction time window reaches to 200 ms, and then the prediction time window is fixed to the maximum 200 ms. Note the selection of the ratio and the maximum prediction time window is a tradeoff between the prediction accuracy and the early prediction time advantage. The prediction accuracy for the IM slip trajectory is denoted by the mean absolute error (MAE), which is defined as: 1 !W MAE = W p=1 p

)

*

sj − sˆj 1 !l | | p l j=1 sj p

p

(26)

p

where sj and sˆj are the actual and predicted IM slip, respectively; W is the total number of prediction windows and l is the prediction length for every prediction window. 4.2. Simulation results and analysis In Fig. 5, the heavy load region is in the dashed box. All the loads in this area are composite loads whose parameters are identified at different operating time. The proposed method is conducted on this region after a large disturbance. Here, at the current operation state, the load proportion for the composite loads is: Z–5%, I–5%, P–10%, and IM–80%. The detailed results and analysis are given as follows. 4.2.1. Symmetrical fault scenarios A 3-phase fault is applied at the middle of the line 4–14 at time 1 s, and the fault is cleared by opening the line at different times: 1.100 s, 1.106 s, 1.112 s and 1.115 s. Fig. 6 gives the voltages and IM slips in the box area.

Bus no.

SEP

UEP

Bus no.

SEP

UEP

bus3 bus4 bus15 bus16 bus18 bus21

0.0071 0.0090 0.0086 0.0079 0.0086 0.0068

0.2053 0.1620 0.1683 0.1829 0.1686 0.2125

bus24 bus25 bus26 bus27 bus28

0.0074 0.0063 0.0067 0.0070 0.0066

0.1962 0.2291 0.2184 0.2080 0.2214

As seen in Fig. 6(a), the system is short-term voltage stable for fault clearing time (tcl ) = 1.100 s. All the curves of IM slip return to the post-fault SEPs. At tcl = 1.106 s, Fig. 6(b) shows the IM slip at bus 4 increases continuously and this leads to STVIS at bus 4. Except for the bus 4, all bus voltages are transiently stable. In Fig. 6(c), with tcl = 1.112 s, both buses 4 and 15 lost voltage stability due to stalling of IMs at these two buses. When tcl is 1.115s, it can be seen in Fig. 6(d) that all IMs are stalled and all buses are transiently unstable. For this situation, the STVS and rotor angle stability are both transiently unstable. And owing to the out-of-step of generators, it causes voltages to produce chaotic behaviors. The post-fault SEPs and UEPs of all the IMs for tripping line 4–14 are shown in Table 1. Based on the proposed method, the unstable status can early be predicted. When tcl = 1.106 s, only bus 4 voltage is transiently unstable and the corresponding IM slip passes through the post-fault UEP, 0.1620. This unstable status can be detected at time 3.76 s before the IM slip goes across 0.1620. When tcl = 1.112s, the IM slips at bus 4 and bus 15 will pass through their post-fault SEPs (0.1620 and 0.1683, respectively) and STVIS occurs. The unstable status of bus 4 and 15 are early detected at time 2.00 s and 3.02 s, respectively. At tcl = 1.115 s, all IMs stall and the postfault system loses STVS. In actual, to determine the first unstable load bus (weak bus) is more important for operators to take emergency control actions to prevent cascading outages. By using the proposed method, the STVIS detection time for the weak load bus is given in Fig. 6. The performance of the proposed method is given in Table 2 where the weak bus (WB), the corresponding UEP, the STVIS detection time (DT), and the MAE of trajectory extrapolation are listed. Furthermore, as compared with the Lyapunov exponent based STVIS detection method in Ref. [22], the validity of the proposed method is further demonstrated. Table 2 shows that the proposed trajectory extrapolation has a good fitting accuracy. The STVIS can

240

H. Yang et al. / Electric Power Systems Research 160 (2018) 234–242 Table 3 Performance of the STVIS prediction method. Faulty line

WB

UEP

DT(s)

DT(s) [22]

MAE (%)

Status

21–22 16–19 5–6 6–7 10–11 6–11 2–3 23–24 28–29 16–21

bus21 – bus4 bus4 bus4 bus4 bus3 bus21 bus28 bus21

0.1779 – 0.1627 0.1658 0.1649 0.1642 0.1555 0.2052 0.1825 0.2197

1.61 0 1.87 – 1.96 2.03 1.72 2.17 1.96 1.83

2.38 2.64 2.76 – 2.86 3.00 2.64 2.73 1.98 2.54

0.83 – 1.51 0.91 1.23 1.49 1.42 2.10 2.29 1.44

Unstable Unstable Unstable Stable Unstable Unstable Unstable Unstable Unstable Unstable

Table 4 Prediction results for asymmetrical fault at line 4–14. Fault duration

Fault type

WB

UEP

DT(s)

DT(s) [22]

MAE (%)

400 ms 300 ms 200 ms

i ii iii

bus4

0.1620

2.01 1.68 1.69

3.02 2.86 3.08

0.52 0.61 0.79

the system goes into chaotic state, the STVIS has taken place. Thus, the Lyapunov exponent based method cannot capture the essential mechanism of STVIS and the detection time is not fast enough for making emergency control. To further illustrate the validity of the proposed method, fault scanning is performed on the transmission lines. Here, the STVIS prediction results considering the 3-phase fault at the middle of the ten maximum load lines and 120 ms fault duration are shown in Table 3. For fault line 6–7, the post-fault system is transiently stable. With the fault on line 16–19, the post-fault system has no equilibrium point. That is because when line 16–19 is tripped, generators 4 and 5 are separated from the whole power system, which results in much active power deficiency so as to make the postfault power flow unsolvable. This scenario can earlier be affirmed in the stage of online contingency analysis before the fault really occurs. For other unstable cases, it is again verified that the proposed method can detect the unstable status in advance and the trajectory extrapolation using the proposed online learning LSSVM has a good accuracy.

Fig. 6. Simulation results. (a)–(d) Fault clearing time at 1.100 s, 1.106 s, 1.112 s and 1.115 s, respectively.

Table 2 Performance of the prediction results for fault at line 4–14. Fault duration

WB

UEP

DT (s)

DT(s) [22]

MAE (%)

Status

100 ms 106 ms 112 ms 115 ms

bus 4

0.1620

– 3.76 2.00 1.88

– n/a n/a 3.18

2.06 0.41 0.54 0.66

Stable Unstable Unstable Unstable

be detected early for the three unstable cases. However, the detection time based on the Lyapunov exponent method can be given only for the most serious case with 115 ms fault duration, and the detection time is 3.18s. In fact, the Lyapunov exponent method for detecting STVIS is only suitable for the scenario that the system enters chaotic behaviors. And the chaotic behaviors of voltages in system are mainly driven by rotor angle instability instead of the instability of dynamic loads. In Fig. 6(d), it is clearly seen that before

4.2.2. Asymmetrical fault scenarios For the same line 4–14, three asymmetrical fault types are applied, which are: (i) single phase to ground fault, (ii) phase to phase fault, and (iii) phase to phase to ground fault. In order to simulate the STVIS situation, the fault durations are 400, 300 and 200 ms, respectively. These fault durations may appear when the backup protection is considered. The voltages of bus connecting the composite loads are shown in Fig. 7, where the STVIS detection time based on the proposed method in this study is shown by the dark dotted line while the detection time based on Lyapunov exponent method is denoted by the red dotted one. It is clearly seen in Fig. 7 that the Lyapunov exponent based method for detecting STVIS works when the post-fault system enters chaotic state. The performance of the proposed method for these three cases is given in Table 4. In Table 4, the high accuracy of trajectory extrapolation can be seen by using the proposed LSSVM regression and the weak load bus can early be detected before the post-fault system enters chaos. Therefore, the proposed method also has a robust capability to predict STVIS status under different fault types. 4.2.3. STVIS prediction with different load composition Considering the load composition has great influence on STVS, the performance of the proposed method under different composition of loads is illustrated in Table 5. The different load proportion in

H. Yang et al. / Electric Power Systems Research 160 (2018) 234–242

241

composite load constituting IM and constant-power load has more severe impact on STVS than the load consisting of IM and voltagedependent load. That is because the voltage-dependent load will reduce power demand when its’ terminal voltage drops, and the voltage drop deterioration is alleviated. And the effects considering the different proportions of IM are also given in the last three rows of Table 5. It shows that the more IMs in system, the more severe STVS when system suffers large disturbance. 4.3. Response speed of the online STVS prediction The response speed of the proposed method is fast enough for real-time application. Based on presumed contingency analysis, the SEPs and UEPs of the composite loads for the post-fault system are obtained in advance at each operating state. When the fault really occurs, by collecting the initial training data in the 300 ms time window in a sampling rate of 10 ms, the determination of the embedding dimension in the phase space construction and the parameter tuning of LSSVM only take about 70 ms. The time spent in each prediction window in the rolling prediction process only needs about 10 ms. Note that since the STVIS prediction is performed at each bus station with composite load individually, there is no time delay for communication, which is more suitable for real-time post-fault STVIS prediction. 4.4. Time delay in synchronous measurements As illustrated before, the procedure in the left box of Fig. 3 is carried out to create a look-up table based on presumed contingency analysis. The procedure is run continuously (say for every 5-min) under normal operation condition based on the remote synchronous PMU measurements. Considering the remote synchronous measurements, the latency including the time delay of phasor measurement, communication delay and phasor data concentrator (PDC) delay is ranged in a few hundred milli-seconds. The latency is 167–217 ms in Ref. [41]. Compared with the running cycle in this procedure, this time delay can be ignored in the operation center. For the procedure in the right-hand box of Fig. 3, we only need the real-time local measurements and hence the delay can also be ignored.

Fig. 7. (a) Voltges for single phase to ground fault at tcl = 1.40 s, (b) voltges for phase to phase fault at tcl = 1.30 s, (c) Voltges for phase to phase to ground fault at tcl = 1.20 s. Table 5 Prediction results considering different load composition. Composite load IM

Z

I

P

80% 80% 80% 70% 60% 50%

20% 0 0 0 0 0

0 20% 0 0 0 0

0 0 20% 30% 40% 50%

WB

UEP

DT(s)

DT(s) [22]

MAE (%)

bus4 bus4 bus4 bus4 bus4 bus4

0.1625 0.1621 0.1616 0.1888 0.2240 0.2721

2.09 1.77 1.57 1.61 2.03 2.32

– 3.07 2.47 2.29 2.44 2.43

0.37 0.58 0.16 0.65 0.47 0.56

the composite load model is set in the table. Here, we consider the 3-phase fault at the middle of line 4–14 and a 200 ms fault duration. From the table, it is clearly seen that the proposed method can effectively detect the STVIS status in advance compared to the detection time with the Lyapunov exponent based method, and the proposed trajectory extrapolation method has a high accuracy. From the first three load compositions in the table, it shows that the

5. Conclusion Considering the composite load, this paper proposes a PMUbased method to predict the impending short-term voltage instability (STVIS) based on the mechanism of short-term voltage stability (STVS). First, online presumed contingency analysis using synchronous measurements in operation center is performed to build a look-up table consisting of the contingency, the corresponding SEP and UEP of the induction motor (IM) for each composite load. Second, with local measurements, the time series of IM slip is online computed by using a forward Euler method. Then a new trajectory extrapolation method is proposed based on online sequential learning LSSVM to predict the future time-series of IM slip in a rolling way. Finally, based on IM stability mechanism, the STVIS status can be early detected by monitoring whether the predicted IM slip reaches the IM’s UEP in the looked-up table. The weak load buses of STVIS can be picked out earlier so as to help make emergency control accordingly. Application of the method on the New England 39-bus system shows that the extrapolated trajectory method presents a high fitting accuracy and the response speed of the proposed method is fast enough for online use. The results show that the proposed method is more effective in STVIS identification and the early detection time of STVIS than the Lyapunov-exponent

242

H. Yang et al. / Electric Power Systems Research 160 (2018) 234–242

based method. The STVIS status can be detected effectively under both symmetrical and asymmetrical faults. Acknowledgment This work was supported by the National Natural Science Foundation of China (No. 51177093). References [1] P. Kundur, J. Paserba, V. Ajjarapu, G. Andersson, A. Bose, C. Canizares, et al., Definition and classification of power system stability IEEE/CIGRE joint task force on stability terms and definitions, IEEE Trans. Power Syst. 19 (3) (2004) 1387–1401. [2] V. Ajjarapu, C. Christy, The continuation power flow: a tool for steady state voltage stability analysis, IEEE Trans. Power Syst. 7 (1) (1992) 416–423. [3] P. Sauer, M.A. Pai, Power system steady-state stability and the load-flow Jacobian, IEEE Trans. Power Syst. 5 (4) (1990) 1374–1383. [4] P. Kessel, H. Glavitsch, Estimating the voltage stability of a power system, IEEE Trans. Power Deliv. 1 (3) (1986) 346–354. [5] Y. Sekine, H. Ohtsuki, Cascaded voltage collapse, IEEE Trans. Power Syst. 5 (1) (1990) 250–256. [6] J.A. Diaz de Leon II, C.W. Taylor, Understanding and solving short-term voltage stability problems, Proc. IEEE PES Summer Meet (2002) 21–25. [7] J.W. Shaffer, Air conditioner response to transmission faults, IEEE Trans. Power Syst. 12 (2) (1997) 614–621. [8] V. Stewart, E.H. Camm, Modeling of stalled motor loads for power system short-term voltage stability analysis, Proc. IEEE PES General Meeting (2005) 1887–1892. [9] K. Kawabe, K. Tanaka, Analytical method for short-term voltage stability using the stability boundary in the P–V plane, IEEE Trans. Power Syst. 29 (6) (2014) 3041–3047. [10] Z.J. Wang, X.Y. Wang, C.Y. Chung, An analytical method for calculating critical voltage sag clearance time of induction motors, IEEE Trans. Power Deliv. 27 (4) (2012) 2412–2414. [11] E.G. Potamianakis, C.D. Vournas, Short-term voltage instability: effects on synchronous and induction machines, IEEE Trans. Power Syst. 21 (2) (2006) 791–798. [12] I.A. Hiskens, D.J. Hill, Energy functions, transient stability, and voltage behavior in power systems with nonlinear loads, IEEE Trans. Power Syst. 4 (2) (1989) 1525–1533. [13] K.L. Praprost, K.A. Loparo, An energy function method for deter-mining voltage collapse during a power system transient, IEEE Trans. Circuits Syst. 41 (10) (1994) 635–651. [14] D.J. Shoup, J.J. Paserba, C.W. Taylor, A survey of current practices for transient voltage dip/sag criteria related to power system stability, Proc. IEEE PES Power Systems Conf Expo (2004) 1140–1147. [15] Z. Qin, J. Davidson, A.A. Fouad, Application of artificial neural networks in power system security and vulnerability assessment, IEEE Trans. Power Syst. 9 (1) (1994) 525–532. [16] L.S. Moulin, A.P.A. da Silva, M.A. El-Sharkawi, R.J. Marks II, Support vector machines for transient stability analysis of large-scale power systems, IEEE Trans. Power Syst. 19 (2) (2004) 818–825. [17] F.R. Gomez, A.D. Rajapakse, U.D. Annakkage, I.T. Fernando, Support vector machine-based algorithm for post-fault transient stability status prediction using synchronized measurements, IEEE Trans. Power Syst. 26 (3) (2011) 1474–1483. [18] K. Sun, S. Likhate, V. Vittal, V. Sharma Kolluri, S. Mandal, An online dynamic security assessment scheme using phasor measurements and decision trees, IEEE Trans. Power Syst. 22 (4) (2007) 1935–1943.

[19] L. Zhu, C. Lu, Y. Sun, Time series shapelet classification based online short-term voltage stability assessment, IEEE Trans. Power Syst. 31 (2) (2016) 1430–1439. [20] Y. Xu, R. Zhang, J. Zhao, Z.Y. Dong, D. Wang, H. Yang, K.P. Wong, Assessing short-term voltage stability of electric power systems by a hierarchical intelligent system, IEEE Trans. Neural Netw. Learn. Syst. 27 (8) (2016) 1686–1696. [21] Y. Xu, Z.Y. Dong, J.H. Zhao, P. Zhang, K.P. Wong, A reliable intelligent system for real-time dynamic security assessment of power systems, IEEE Trans. Power Syst. 27 (3) (2012) 1253–1263. [22] S. Dasgupta, M. Paramasivam, U. Vaidya, V. Ajjarapu, Real-time monitoring of short-term voltage stability using PMU data, IEEE Trans. Power Syst. 28 (4) (2013) 3702–3711. [23] L. Perira, D. Kosterev, P. Mackin, D. Davies, J. Undrill, An interim dynamic induction motor model for stability studies in the WSCC, IEEE Trans. Power Syst. 17 (4) (2002) 1108–1115. [24] X. Liang, W. Xu, C. Chung, W. Freitas, K. Xiong, Dynamic load models for industrial facilities, IEEE Trans. Power Syst. 27 (1) (2012) 69–80. [25] R. He, J. Ma, D.J. Hill, Composite load modeling via measurement approach, IEEE Trans. Power Syst. 21 (2) (2006) 663–672. [26] J. Ma, D. Han, R. He, Z.Y. Dong, D.J. Hill, Reducing identified parameters of measurement-based composite load model, IEEE Trans. Power Syst. 23 (1) (2008) 76–83. [27] S. Son, S.-H. Lee, D.-H. Choi, K.-B. Song, J.-D. Park, Y.-H. Kwon, K. Hur, J.-W. Park, Improvement of composite load modeling based on parameter sensitivity and dependency analyses, IEEE Trans. Power Syst. 29 (1) (2014) 242–250. [28] P. Ju, F. Wu, Z.Y. Shao, X.P. Zhang, H.J. Fu, P.F. Zhang, N.Q. He, J.D. Han, Composite load models based on field measurements and their applications in dynamic analysis, IET Gen. Transm. Distrib. 1 (5) (2007) 724–730. [29] L. Fan, Y. Wehbe, Extended Kalman filtering based real-time dynamic state and parameter estimation using PMU data, Electr. Power Syst. Res. 103 (2013) 168–177. [30] J.E. Tate, T. Cutsem, Extracting steady state values from phasor measurement unit data using FIR and median filters, Proc. IEEE Power Syst. Conf. Expo (2009) 1–8. [31] C.J.C. Burges, A Tutorial on Support Vector Machines for Pattern Recognition, Kluwer, Norwell, MA, 1998. [32] J.A.K. Suykens, J. Vandewalle, Least squares support vector machine classifiers, Neural Process Lett. 9 (3) (1999) 293–300. [33] F. Takens, Detecting strange attractors in turbulence, Lecture Notes Math. 898 (1) (1981) 366–381. [34] L. Cao, Practical method for determining the minimum embedding dimension of a scalar time series, Physica D 110 (1–2) (1997) 43–50. [35] A.H. Sayed, Fundamentals of Adaptive Filtering, Wiley, New York, 2003. [36] B. Milosevic, M. Begovic, Voltage-stability protection and control using a wide-area network of phasor measurements, IEEE Trans. Power Syst. 18 (1) (2003) 121–127. [37] Y. Wang, I.R. Pordanjani, W. Li, W. Xu, T. Chen, E. Vaahedi, J. Gurney, Voltage stability monitoring based on the concept of coupled single-port circuit, IEEE Trans. Power Syst. 26 (4) (2011) 2154–2163. [38] M.A. Pai, Energy Function Analysis for Power System Stability, Kluwer, Norwell, 1989. [39] Power Factory Manual, DIgSILENT Power Factory Version 15.0. [40] S. Wang, W. Gao, J. Wang, J. Lin, Synchronized sampling technology-based compensation for network effects in WAMS communication, IEEE Trans. Smart Grid 3 (2) (2012) 837–845. [41] J.C. Cepeda, J.L. Rueda, D.G. Colomé, I. Erlich, Data-mining-based approach for predicting the power system post-contingency dynamic vulnerability status, Int. Trans. Electr. Energy Syst. 25 (10) (2015) 2515–2546.