- Email: [email protected]
Practical stability assessment of distributed synchronous generators under variations in the system equilibrium conditions
Electrical Power and Energy Systems 55 (2014) 275–284
Contents lists available at ScienceDirect
Electrical Power and Energy Systems journal homepage: www.elsevier.com/locate/ijepes
Practical stability assessment of distributed synchronous generators under variations in the system equilibrium conditions Roman Kuiava a,⇑, Rodrigo A. Ramos b, Hemanshu R. Pota c, Luis F.C. Alberto b a
Federal University of Parana (UFPR), Department of Electrical Engineering, Rua Cel. F. H. dos Santos, 100, Jardim das Americas, 81531-980 Curitiba, Brazil University of Sao Paulo (USP), Engineering School of Sao Carlos (EESC), Department of Electrical Engineering, Av. Trabalhador Saocarlense, 400, Centro, 13566-590 Sao Carlos, Brazil c University of New South Wales at Australian Defence Force Academy (UNSW@ADFA), School of Information Technology and Electrical Engineering, PO Box 7916, Canberra BC, ACT 2610, Australia b
a r t i c l e
i n f o
Article history: Received 16 October 2012 Received in revised form 11 September 2013 Accepted 24 September 2013
Keywords: Distributed synchronous generator Practical stability Fast varying loads Switched affine system Linear matrix inequality
a b s t r a c t This paper proposes a method to assess the practical stability of power distribution systems with synchronous generators subject to changes in the system equilibrium conditions due to fast varying loads. The concept of practical stability deals with two known state-space regions X1 (which contains all the initial conditions reflecting the perturbations at which the system is subject during its operation) and X2 (which represents the operating security region of the power distribution system) satisfying X1 X2. The practical stability problem and the focus of this paper is to determine under which conditions the system trajectories will be confined into a security region of operation for a certain time interval of interest, as the equilibrium point of the model changes. This study was carried out using a mathematical model of the distribution system with synchronous generators in the form of a switched affine system. This proposed model is capable of describing the system behavior over a certain period within which changes on the equilibrium conditions of the system can occur. Sufficient conditions for the power distribution system with synchronous generators described as a switched affine system to be practically stable with respect to its operating security region X2 are given in the form of matrix inequalities constraints. The results, obtained for the model of a cogeneration plant of 10 MW added to a distribution network constituted by a feeder and six buses, show that the less stringent properties of the concept of practical stability can be very well-suited to the security analysis of power systems subjected to frequent variations in the load level. Ó 2013 Elsevier Ltd. All rights reserved.
1. Introduction The interest in distributed generation (DG) has increased considerably around the world over the years, driven mainly by market deregulation, continuous demand rise, technological developments and concerns on environment impact. In Brazil, for example, the strong economical incentives for ethanol production are likely to boost the number of cogeneration power plants connected to both subtransmission and distribution grids. Located mainly in southeast and center-west areas of Brazil, this sector of the national industry had already considerable experience in the usage of steam turbines, which have been extensively employed within its internal production [1–3]. As a result, several synchronous generators driven by steam turbines have been connected to the system using a byproduct of the sugar cane production to
⇑ Corresponding author. E-mail addresses: [email protected] (R. Kuiava), [email protected] (R.A. Ramos), [email protected] (H.R. Pota), [email protected] (L.F.C. Alberto). 0142-0615/$ - see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.ijepes.2013.09.013
feed the boilers, which constitutes a very sustainable form of electrical energy production. In view of these facts, within the new environment provided by DG, the networks may present new and unforeseen problems that must be appropriately dealt with. Recent papers, for example, present research studies considering distinct scenarios and technical factors that can affect the dynamic performance of power distribution systems with distributed synchronous generators in terms of transient stability [1,4,5] and small-signal stability in both balanced [6–10] and unbalanced operating conditions [11,12]. By means of modal and eigenvalues analysis, papers [8,7,11] also show that the presence of sustained electromechanical oscillations in the response of these systems can have a significant impact on the power quality across the grid, which means that this impact has to be assessed before granting permission to connect a generator to the grid. This paper extends the security analysis mentioned in the previous paragraph by assessing the behavior of the system as the load levels vary intermittently over time, i.e., fast varying loads (FVLs) type, using the concept of practical stability (which has a for-
276
R. Kuiava et al. / Electrical Power and Energy Systems 55 (2014) 275–284
mal definition that will be shown later in the paper). FVLs in distribution systems can cause undesirable disturbances or an increase in energy losses, as discussed in [13]. Different types of FVLs can be found, for instance, in industrial [14,15] and medicine fields [13]. From a practical viewpoint, the power distribution system must operates within a security region, which is determined from the allowable range of the state variables (rotor angle, rotor speed, generator terminal voltage, and others), even when facing fast variations in its load levels. Hence, in a realistic scenario, the problem of interest and the focus of this paper is to determine under which conditions the system trajectories will be confined into a security region of operation for a certain time interval of interest, as the equilibrium point of the model changes, and this feature can be captured adequately by the concept of practical stability. Generally speaking, the concept of practical stability [16–19] deals with two sets X1 Rn and X2 Rn satisfying X1 X2, which are specified for the initial state and the system trajectories, respectively. These two sets are specified in terms of physical limitations of the system variables and they must contain all the equilibrium points in which the system is expected to operate within a certain period of time I. For our purposes, the set X2 is a representation of the operating security region of the power distribution system, while set X1 involves all the initial conditions reflecting the perturbations at which the system is subject during its operation. The sets X1 and X2 can be determined from a practical knowledge of the system operation and/or via numerical simulations analysis. The practical stability requires that if the initial state is in X1, then the system state should stay in X2 for all t 2 I [17]. Hence, unlike the classical stability definition which is based on the existence of X1 for any X2, here both of the sets X1 and X2 are fixed and so they do not vary. Stability with respect to a set, rather than a particular point, is then the basis of the practical stability concept. A mathematical model of the distribution system with synchronous generators in the form of a continuous-time switched affine system is proposed in this paper to be used for practical stability studies. This proposed model is capable of describing the system behavior over a certain period within which changes on the equilibrium conditions of the system can occur. Sufficient conditions for the power distribution system described as a switched affine system to be practically stable with respect to a known operating security region X2 are given in the form of matrix inequalities constraints. The paper is organized as follows: Section 2 presents a traditional arrangement of distributed system with synchronous generation, which serves as motivation for the application of the proposed stability assessment; Section 3 presents the problem formulation and a motivation example; Section 4 presents the proposed modeling approach used to describe the original power distribution system model by means of a switched affine system; Section 5 presents the fundamentals of the practical stability concept; Section 6 presents the proposed systematic method to assess the practical stability; Section 7 presents some tests and results and Section 8 presents the conclusions and final remarks.
2. The study system A typical arrangement of a distributed system with synchronous generation, usually viewed in cogeneration plants in the sugarcane industry of Brazil, is shown in Fig. 1. This network consists of a 132-kV, 60-Hz subtransmission system, which feeds a 33-kV distribution system through one 132/33-kV, D/Yg transformer
Fig. 1. Diagram of the study system.
(T1). The cogeneration plant (G) is represented by a gas-turbine driven synchronous generator injecting 10 into the grid. A local load (L) is connected to the bus 3. This local load is represented by a conventional static load model with constant impedance, which is calculated from the values of the active and reactive power demand level of the bus 3 and the nominal amplitude voltage of this same bus. We also assume this local load L as being a simplified representation of several different loads coupled to the bus 3 and each of them can be active or not during a period of operation of the distribution power system. Hence, the operating conditions of the system are supposed to be divided according to which load it supplies. The detailed description of the mathematical model of this system is presented in Section 7 of tests and results. In this point of the paper, it is only important to point out that, for dynamic stability studies, the power distribution system with synchronous generator under investigation (or even more general power systems) can be mathematically described by a model in the state-space form given by
_ xðtÞ ¼ f ðxðtÞ; kÞ;
ð1Þ
where xðtÞ 2 Rn and k 2 Rm are, respectively, the vector with the state variables (such as, rotor angle and rotor speed) and the vector with the parameters that are dependent of the operating conditions, such as the system loads [20], and f : Rn Rm ! Rn is a nonlinear vector function. 3. The problem formulation Consider again the power distribution system shown in Fig. 1 and a mathematical model of it written in the state-space form (1). An equilibrium point xe of (1) (such that f(xe, k) = 0) is calculated by using the power-flow outputs for a specific operating condition (or load level), i.e., for a fixed value of k. Let us suppose that an operating condition k = k1 (associated to a local load L = L1) corresponds to the equilibrium point xe1 of (1) and k = k2 (in which the local load is changed to L = L2) to the equilibrium xe2 . Hence, for a total of N load conditions, we have N equilibrium points given by xe1 ; xe2 ; . . . ; xeN 2 Rn . These operating conditions can, obviously, vary among them along a period of time. The frequency of occurrence of these load changes can be forecasted in advance and this information can be used to estimate the minimum amount of time in which the system stays in a particular operating condition (and this information will be explored later in the paper by means of the concept of dwell-time for switched systems). Typical loads that demand large amounts of varying power from the network include, for example, electric arc furnaces [15]. Obviously, this type of load may require a more accurate model [15] (different from the static load model with constant impedance adopted in the paper) to analyze its impacts on system dynamics, such as, a hidden Markov-based model [21] or an adaptive model [22], but this level of model complexity is not taken under consideration in this paper. Other common fast
R. Kuiava et al. / Electrical Power and Energy Systems 55 (2014) 275–284
varying loads are rolling mills and spot welders in industrial area. An overview about different types of fast varying loads in distribution systems is presented in [13]. Let us also assume that linearized models of (1) in the vicinity of these equilibrium points are supposed to be proper descriptions of the system with respect to the load variations of interest. These linearized models (Ri, . . . , RN) are given in the state-space form
Ri : Dx_ i ðtÞ ¼ Ai Dxi ðtÞ;
ð2Þ
where Ai 2 Rnn is the ith state matrix, Dxi ðtÞ ¼ xðtÞ xei and i 2 S, where S = 1, 2, . . . , N. So, Ri correspond to the linear approximation of (1) in the vicinity of the equilibrium point xei . Let us consider a period of time I = [t0, tf), where t0 is the initial time and tf is a finite constant, and the sets X1 Rn and X2 Rn , X1 X2. The set X1 includes the initial conditions created by a number of perturbations to which the system may be subject during its operation and the set X2 represents the security region of the system operation. Now, let us say that the load conditions change from Li to Lj (which means that Ri changes to Rj, for some i, j 2 S, i – j, in the time instants t1, t2, . . . , tk, . . . 2 I(t0 < t1 < t2<
Fig. 2. Load variations and the system description as time goes by.
277
frequency oscillation is smaller for the case where faster load variations are considered. Furthermore, for this same case it is possible to notice that the peak values of the oscillations tend to increase as time goes by, and this fact is an indicative that the system trajectory will not be always into the security region of operation. From this motivation example, this paper provides some sufficient conditions guaranteeing that the system trajectories (with the system subject to load variations, as illustrated in Fig. 2) will be confined into the security region of operation X2 for all the period of time I. This is done by representing the set of linearized models (2) as a switched affine system, as presented in next section. 4. Modeling the power distribution system as a switched affine system Let us consider again a power distribution system with synchronous generators modeled in the state-space form (1) and a set of N linearized models Ri, . . . , RN given in the state-space form (2) obtained from the linearization of (1) with respect to N different equilibrium points of this system. Let us also consider a period of time I = [t0, tf). We assume that, for all t 2 [tk, tk+1), where k = 0, 1, 2, . . . and t0, t1, . . . , tk, . . . 2 I, the system behavior is described by the linearized model Rik ; ik 2 S, and at tk+1, the system description changes from Rik to Rikþ1 ; ikþ1 2 S; ikþ1 – ik . As discussed in the previous section, this means that, at tk+1, the load conditions change from Lik to Likþ1 . Before presenting the proposed system representation in the form of a switched affine system, it is important to discuss the main reason of using a switched affine system for stability analysis instead of the set of linearized models (2). This is done in next subsection. 4.1. The main difficulty of analyzing the system stability via linearized models From the Lyapunov stability theory, the linearized models
R1, . . . , RN are stable if there exists a set of positive definite matrices P1, . . . , PN such that ATi Pi þ Pi Ai 0; i ¼ 1; . . . ; N [24]. These conditions guarantee local stability of the equilibrium points xe1 ; . . . ; xeN in the sense of Lyapunov stability theory, but they do not guarantee that the system trajectories will be confined into a certain security region of operation X2 over the period of time I in the occurrence of changes on the equilibrium conditions. One of the main difficulties here of taking account the transitions of the equilibrium operating conditions to the stability study, is the fact that the state vector also changes during an operating condition variation (from Dxi(t) to Dxj(t)), once that the coordinate axis of each linearized model has a different origin. To better clarify this point, consider a hypothetical second order nonlinear dynamical system, where the state vector is x(t) = [u(t) v(t)]T. Let us assume this system has two equilibrium points, xe1 ¼ ½ue1 v e1 T and xe2 ¼ ½ue2 v e2 T . Hence, the state vectors of the linearized models calculated in the vicinity of these equilibrium points are given by Dx1(t) = [Du1(t) Dv1(t)]T and Dx2(t) = [Du2(t) Dv2(t)]T, being Du1 ðtÞ ¼ uðtÞ ue1 , Du2 ðtÞ ¼ uðtÞ ue2 , Dv 1 ðtÞ ¼ v ðtÞ v e1 and Dv 2 ðtÞ ¼ v ðtÞ v e2 . Fig. 4(a) illustrates the coordinate axis (u v) of this nonlinear system, where its origin is at (0, 0), as well as, the coordinate axes (Du1 Dv1) and (Du2 Dv2) of the linearized models with respect to the equilibriums xe1 and xe2 , where their origins are at xe1 and xe2 , respectively. Notice that for all the time interval [t0, t1) the system trajectory is obtained from the linearized model R1. At t = t1, the operating conditions of the system is changed, which means that from t1 and on the system trajectory is obtained from R2. The coordinate axis was also changed (from (Du1 Dv1) to (Du2 Dv2)) at
278
R. Kuiava et al. / Electrical Power and Energy Systems 55 (2014) 275–284
Fig. 3. System trajectory (rotor speed, in Hertz) by considering that the time ellapsing two consecutive load variations is no smaller than 1 s (dashed line) and 5 s (full line).
Fig. 4. (a) System trajectory given by linearized approximations with different coordinate axes and (b) system trajectory given by linearized approximations with a same coordinate axis.
the time instant t1, and this fact increases the complexity of analyzing the system stability. Based on the above discussion, next section provides a possible solution for the pointed problem. 4.2. The proposed solution
ations (or to consecutive equilibrium conditions variations). As previously mentioned, the frequency of occurrence of these load changes can be forecasted in advance and this information is adopted to estimate the dwell-time T. For that, we model the set of affine systems R1 ; . . . ; RN as a switched affine system, as presented in next subsection. 4.3. The proposed model for practical stability studies
The difficulty presented in previous section is overcome in this paper by, firstly, maintaining the coordinate axis of the nonlinear system (1) to the linearized models (2), which is done by transforming each of R1, . . . , RN in an affine system. This transformation process is given by
The proposed mathematical model of the power distribution system with synchronous generators for practical stability studies is given by the switched affine system
Ri : Dx_ i ðtÞ ¼ Ai Dxi ðtÞ )
_ xðtÞ ¼ ArðtÞ xðtÞ þ brðtÞ ;
_ x_ ei ¼ Ai ðxðtÞ xei Þ ) xðtÞ _ xðtÞ ¼ Ai xðtÞ Ai xei )
ð3Þ
b i : xðtÞ _ R ¼ Ai xðtÞ þ bi ; where bi ¼ Ai xei . Notice that the equilibrium point of Ri is xei , while the equilibrium of Ri is the origin Dxi(t) = 0. Let us now consider again the hypothetical second order dynamical system presented in last subsection and its trajectory for a certain initial condition, as illustrated now in Fig. 4(b). In this case, for all the time interval [t0, t1) the system trajectory is obtained from the affine system R1 . At t = t1, the operating conditions of the system is changed, which means that from t1 and on the system trajectory is obtained from R2 . Different from the situation illustrated in Fig. 4(a), the coordinate axis was not changed at t = t1 (in fact, the affine systems R1 and R2 are refered to the same coordinate axis (u v) of the nonlinear system). In Problem 1 formulated in Section 3, we focused on the conditions in which the system trajectory will be confined into the security region of operation X2 for all t 2 I. These conditions are evaluated in this paper in terms of a minimum time T (that will be called dwell-time) ellapsed between two consecutive load vari-
xðt0 Þ ¼ x0 ;
ð4Þ
where r: I ? S is called switching signal. Notice that for all the time in which r(t) = i, i 2 S, the trajectory of (4) is in fact the trajectory of the affine system Ri , whose equilibrium point is xei . Also, the switching signal r(t) is given by
rðtÞ ¼ ik 2 S; 8t 2 ½tk ; tkþ1 Þ; tkþ1 tk P T;
ð5Þ
where t1, t2, . . . , tk, tk+1, . . . 2 I are the instants of time in which the system equilibrium conditions change and T is the minimum time ellapsed between two consecutive load variations. Next subsection validates this proposed model for some numerical simulations. 4.4. Validation of the proposed model In this subsection, the approximation of the power distribution system with synchronous generators model in the form of the switched affine system (4) with r(t) as as defined by (5) is validated by means of some numerical simulations. For that, consider again the study system described in Section 2. A switched affine system was obtained by considering three load conditions L1, L2 and L3, whose numerical values are shown in Section 7. Hence, for each of these load conditions, an affine system was obtained
R. Kuiava et al. / Electrical Power and Energy Systems 55 (2014) 275–284
by following (3), which means that the switched affine system is constituted by three affine systems R1 , R2 and R3 . Fig. 5(a) and (b) shows the system trajectory (rotor speed, in Hertz) obtained from the nonlinear model (1) (gray full line) and the system trajectory obtained from the approximated model in the form of a switched affine system (black dashed line). These figures were obtained by considering the time ellapsing between two consecutive load variations is no smaller than 1 s (Fig. 5(a)) and 5 s (Fig. 5(b)), which means that T in (5) is equal to 1 s for the first case and 5 s for the second case. As it can be seen in Fig. 5(a) and (b), the switched affine system describes adequately (at least, for the operating conditions under consideration, which are detailed in Section 7 of tests and results) the system trajectory over a period of time where changes on the equilibrium conditions occur due to load variations. Discrepancies between the trajectories of the nonlinear and affine system models can be found mainly if the operating conditions (equilibrium points) are far from each other. Hence, from now and on we consider the power distribution system with synchronous generators modeled in the form of the switched affine system (4) with r(t) given by (5). Next section provides some sufficient conditions of practical stability for this switched affine system in terms of a minimum time T that must be ellapse between two consecutive load variations. 5. Sufficient conditions of practical stability of switched affine systems
279
Fig. 6. Illustrative example of the practical stability concept.
It is important to emphasize again that the set X2 is assumed to be the representation of a region within which the safe operation of the power system is guaranteed, while the set X1 includes the initial conditions created by a number of perturbations to which the system may be subject during its operation (these perturbations are usually known in advance as the results of the process of contingency screening). Besides, the dwell-time T of r(t) is a practical information that gives an idea on how often changes on the equilibrium conditions occur in the system as time goes by. In this paper, the sets X1 and X2 are assumed to be given in the form of the ellipsoids
X1 ¼
N \
X1i ;
ð6Þ
i¼1
Consider a power distribution system described mathematically as a switched affine system in the form (4) with r(t) given by (5). In the sequence, the practical stability concept is formalized and Fig. 6 illustrates this concept for a second order dynamical system with the sets X1 and X2 given in the form of ellipsoids. Definition 1. The switched affine system (4) with r(t) given by (5) is considered practically stable with respect to the sets X1 Rn and X2 Rn (X1 X2) in the time interval I = [t0, tf), if x(t0) 2 X1 implies x(t) 2 X2, for all t 2 I.
X2 ¼
N [
X2i ;
ð7Þ
i¼1
where
X1i ¼ fx 2 Rn : ðx xei Þ0 Pi ðx xei Þ þ di 6 ag;
ð8Þ
X2i ¼ fx 2 Rn : ðx xei Þ0 Pi ðx xei Þ þ di 6 1g;
ð9Þ
Fig. 5. System trajectory (rotor speed, in Hertz) of the nonlinear model (1) (gray full line) and the system trajectory of the approximated model in the form of a switched affine system (black dashed line): (a) T = 1 s and (b) T = 5 s.
280
R. Kuiava et al. / Electrical Power and Energy Systems 55 (2014) 275–284
being a a positive scalar satisfying a < 1, Pi 2 Rnn a positive definite matrix, di a positive scalar satisfying di < a; xei ¼ A1 i bi the equilibrium point of Ri and i = 1, . . . , N. Next theorem provides sufficient conditions for practical stability of the switched affine system (4) with r(t) given by (5), considering the sets X1 and X2 given in the form of (6) and (7), respectively. This result was previously proposed by the authors in [25,26], so it is presented here without proof. Theorem 1. Let the sets X1 and X2 be given in the form of (6) and (7), respectively. If there exist a scalar TD > 0 and a number l such that the inequalities
1
pffiffiffiffiffiffiffiffiffi 1=a;
N ðt0 ;tf Þ
0
eAi T D Pj eAi T D lPi
ð11Þ 0
eAi T D Pj Dxeij 0 Dxeij Pj Dxeij þ dj ldi
0;
6. The proposed method for practical stability analysis of power distribution systems The proposed method for practical stability analysis of power distribution systems includes the computation of the set X2 in the form of the ellipsoid (7) from the knowledge of the allowable operating range of the system. This allowable operating range is given by
i ¼ 1; . . . ; ng;
ð13Þ
xi
where xi and are, respectively, the lower and upper bounds of the state variable xi. The set X can be written as
X ¼ X i :¼ fx 2 Rn : aik ðx xei Þ 6 1;
k ¼ 1; . . . ; ne g;
Lemma 1. Consider the state-space region Xi defined as (14). The condition x 2 Xi can be written as
k ¼ 1; . . . ; ne :
ð17Þ
Proof. The matrix inequality (17) can be written as
0 a0ik 2
P i
0 : di 1
ð18Þ
Once the constrain (18) is satisfied, the S-Procedure guarantees that for all x such that
x xe i 1
0
Pi
x xe i 1
0
0 di 1
0 a0ik 2
x xe i 1
0;
ð19Þ
ð15Þ
x xe i 1
0:
ð20Þ
Now, notice that (19) and (21) are the matrix form of (16) and (15), respectively. This completes the proof. h The step-by-step procedure proposed below systematizes the use of Theorem 1 to assess the practical stability of power distribution systems with synchronous generators modeled as a switched affine system in the form (4) with r(t) given by (5). This procedure includes the computation of the set X2 X with its size as large as possible. Step 1: Initialization – (i) compute the equilibrium point xei of Ri ; i ¼ 1; . . . ; N; (ii) calculate the matrix Ai and the vector bi ¼ Ai xei , i = 1, . . . , N; (iii) from a known value of TD, compute the maximum number of switchings N ðt0 ; tf Þ that can occur within the time interval of interest I = [t0, tf); and (iv) from the allowable operating range of the system given by (13), compute the equivalent set Xi, i = 1, . . . , N, as defined by (14). Step 2: Specifying the initial conditions – (i) specify a value for a such that a < 1; (ii) compute the upper bound of the inequality (10); and (iii) specify a value for l such that (10) is satisfied.1 Step 3: Checking the practical stability – (i) build the computational representation of the matrix variable Pi, the scalar variables di and ki, i = 1, . . . , N. Solve the LMI optimization problem (OP) presented below. If a feasible solution is found, stop2; otherwise, return to Step 2 and specify a smaller value of a and/or a different value of l.
OP :
minimize
N X ki ; i¼1
ð14Þ
by a proper choice of the row vectors aik 2 Rn , k = 1, . . . , ne, for some i 2 S, where ne is the number of edges of Xi. Hence, the set X can be written equivalently in N different forms, i.e., X ¼ X 1 ¼ X 2 ¼ ¼ X N . These N different forms will be useful to compute set X2 exploring LMIs. The procedure to compute the equivalent forms of (13) is discussed in Section 7 of tests and results. Next Lemma provides sufficient conditions to ensure that the set X2i ; i 2 S, as defined by (9) is bounded by the state-space region X written in the equivalent form Xi.
2 ðx xei Þ0 a0ik aik ðx xei Þ P 0;
a0ik 0; di þ 1
it is true that
Proof. The proof can be found in [25,26]. h In the above Theorem, the parameter N ðt0 ; t f Þ is the maximum number of switching times that can occur on the time interval I and, obviously, it depends on the lenght of I and the value of TD. ðt f t 0 Þ ðt t Þ If is an integer number, then N ðt ; t Þ ¼ fT D 0 1; otherwise, 0 f TD ðt f t 0 Þ N ðt 0 ; t f Þ ¼ int T D . The use of Theorem 1 to assess the practical stability of power distribution systems with synchronous generators modeled as a switched affine system is described in next section.
Pi
is satisfied for all k = 1, . . . , ne.
ð12Þ
X :¼ f½x1 xn 0 2 Rn : xi 6 xi 6 xi ;
ð16Þ
Therefore, the condition x 2 X2i X i is guaranteed if the inequality
#
are satisfied for all i, j 2 S, i – j, where Dxeij ¼ xei xej , and N ðt 0 ; t f Þ is the maximum number of switchings that can occur on the time interval I, then the switched affine system (4) is practically stable with respect to X1 and X2 in the time interval I, for every r(t) satisfying (5) with T = TD.
ðx xei Þ0 Pi ðx xei Þ þ di 1 6 0:
ð10Þ
A0i P i þ Pi Ai 0; "
Let the set X2i be defined as (9). If x 2 X2i , then
subject to ki > 0;
di < a;
ð21Þ
0 Pi ki I;
ð22Þ
A0i Pi þ P i Ai 0;
ð23Þ
Pi
a0ik 0; di þ 1
ð24Þ
1 Several numerical tests in different examples have shown that proper values of l are close to the upper bound of (10). 2 Although the LMI problem OP must be solved with a pre-specified value of a, it is important to emphasize that the size of the set X1 can be manipulated by means of an iteractive process. For example, its size can be enlarged as possible by solving the LMI problem OP repeatedly with increasing values of a until a terminating condition (the unfeasibility of the LMI problem) has been met.
281
R. Kuiava et al. / Electrical Power and Energy Systems 55 (2014) 275–284
"
0
eAi T D Pj eAi T D lPi
i; j ¼ 1; . . . ; N;
0
eAi T D Pj Dxeij þ dj ldi
Dx0eij Pj Dxeij
i – j;
# 0;
ð25Þ
k ¼ 1; . . . ; ne :
If a feasible solution is found in this optimization problem, then the trajectories of the power distribution system with synchronous generators (described in the form of the switched affine system (4)) will be confined in the security region of operation X2 X for all the time interval I, ever so the initial condition is in X1 and the time ellapsed between two consecutive load variations is no smaller than T = TD. It is important to notice that the size of X2 is being enlarged as possible by means of the minimization of the highest eigenvalue of Pi, i = 1, . . . , N, via the inequality on the right of (22). In addition, the specification of l and a are only needed in order to transform the inequality (25) in an LMI. Next section illustrates the application of the above step-by-step procedure to study the practical stability of the system presented in Section 2. 7. Tests and results
Let us consider an affine representation in the vicinity of the equilibrium points xe1 ; xe2 and xe3 , as calculated by (3). Let us also consider a time interval I = [0, 65). Hence, the switched affine system adopted to study the practical stability of the distributed system with synchronous generator under consideration is described by
_ xðtÞ ¼ ArðtÞ xðtÞ þ brðtÞ ;
xð0Þ ¼ x0 ;
ð28Þ
where r: I = [0, 65) ? S = [1, 2, 3], being
A1 ¼
A2 ¼
A3 ¼
0
376:9911
0:2444
0:0667
0
376:9911
0:2413
0:0667
0
376:9911
0:2358
0:0667
;
b1 ¼
;
b2 ¼
;
b3 ¼
376:9911 0:0927
376:9911 0:0944
376:9911 0:0862
7.1. Tests and results for a second order model
:
0:9917 6 w
6 1:0083g:
Consider the power distribution system with synchronous generator presented in Section 2. In this subsection, we consider a second order model (classical model) for describing the study system, where the state variables are the generator rotor angle and angular speed. The state-space model is:
ð26Þ
ð27Þ
where d(t) and x(t) are the generator rotor angle and the angular speed, respectively, so the state vector is x(t) = [d(t) x(t)]0 . The parameters Eh0q ; G11 ; G12 and iB012 are function of the operating conditions, so k ¼ E0q G11 G12 B12 , being E0q the quadrature axis transient internal voltage of the generator, while G11, G12 and B12 are elements of the reduced bus admittance matrix. This matrix reflects the topological characteristics of the distribution network, including the resistance and reactance of the line 2–3, the reactances of the transformers and the load (L) impedance. Finally, ws is the synchronous speed, H is the inertia constant, Pm is the mechanical power (turbine output) and jV1j is the magnitude of the voltage at bus 1. Detailed information regarding the equations of the presented model and their respective parameters can be obtained in [23]. The numerical values of the other model parameters are: ws = 376.99 rad/s, H = 1.5 s, D = 0.2 and Pm = 0.1 p.u. We considered three different load levels:
L1 ¼ PL1 þ jQ L1 ¼ 3:2 þ j1:0 MVA; L2 ¼ PL2 þ jQ L2 ¼ 5:2 þ j2:0 MVA;
ð29Þ
The ranges of d and w in (29) are expressed in rad and p.u., respectively. In order to calculate the set X2 via the optimization problem OP, we have to write the set (29) in the form of Xi, i = 1, 2, 3, as defined by (14). This was done following the ideas presented in [27]. At first, notice that (29) is equivalent to a polytope in R2 whose vertices are defined by H :¼ {c1, c2, c3, c4}, where
c1 ¼
_ ðtÞ ¼ Pm E02q G11 þ E0 q jV 1 jðG12 cosðdðtÞÞ þ B12 sinðdðtÞÞÞ 2Hx DðxðtÞ ws Þ;
;
It was adopted a TD = 6 s. Notice that t0 = 0 s, tf = 65 s and TD = 6 s, which means that N ðt0 ; tf Þ ¼ 10. We consider the following allowable operating range of the system:
X :¼ f½d w0 2 R2 : 0:2618 6 d 6 0:4363;
_ dðtÞ ¼ ws wðtÞ ws ;
;
c3 ¼
0:2618 0:9917 0:4363 0:9917
;
;
0:2618 c2 ¼ ; 1:0083 0:4363 c4 ¼ : 1:0083
A vertex representation of X is defined as the convex hull of c1, c2, c3 and c4. Equivalently, we can define this vertex form of X as (14), for i = 1, 2, 3. We present in details the construction of the set X1. The sets X2 and X3 can be obtained in a similar manner. In order to calculate the row vector a1k 2 R2 , notice that each inequality in (14), i = 1, defines an hyperplane fx : a1k ðx xe1 Þ 6 1g for which
c1 ; c2 2 fx : a11 ðx xe1 Þ ¼ 1g;
c2 ; c3 2 fx : a12 ðx xe1 Þ ¼ 1g;
c3 ; c4 2 fx : a13 ðx xe1 Þ ¼ 1g;
c4 ; c1 2 fx : a14 ðx xe1 Þ ¼ 1g:
For a given set H, the row vector a1k, k = 1, . . . , 4, can be determined by solving the following set of linear systems:
a11 ðc1 xe1 Þ ¼ 1 a11 ðc2 xe1 Þ ¼ 1 a13 ðc3 xe1 Þ ¼ 1 a13 ðc4 xe1 Þ ¼ 1
; ;
a12 ðc2 xe1 Þ ¼ 1
; a12 ðc3 xe1 Þ ¼ 1 a14 ðc1 xe1 Þ ¼ 1; a14 ðc4 xe1 Þ ¼ 1;
yielding the following row vectors
L3 ¼ PL3 þ jQ L3 ¼ 7:2 þ j3:2 MVA:
a11 ¼ ½ 3:033 0 ;
The step-by-step procedure presented in the last section was used to assess the dynamic behavior of the study power distribution system in terms of practical stability. For the Step 1, we have that the system equilibrium points associated to the load levels L = L1, L = L2 and L = L3 are, respectively, xe1 ¼ ½0:1065 10 ; xe2 ¼ ½0:1152 10 and xe3 ¼ ½0:0828 10 .
This completes the Step 1. For the Step 2, it was specified an initial value of a equal to 0.5. As previously commented, numerical tests have shown that proper values of l are close to its upper , where l is the upper bound of bound. So, we specify l ¼ 0:99l (14).
a13 ¼ ½ 2:714 0 ;
a12 ¼ ½ 0 120:5 ; a14 ¼ ½ 0 120:5 :
282
R. Kuiava et al. / Electrical Power and Energy Systems 55 (2014) 275–284
For the Step 3, the LMI optimization problem OP was solved repeatedly with increasing values of a (from its initial value specified in Step 2) until the terminating condition given by the unfeasibility of the LMI problem. For all the process, we adopted l ¼ 0:99l . The final value of a was 0.67. The final solution of the LMI optimization problem OP was: k1 = 9093.4, k2 = 9475.1, k3 = 9093.6, d1 = 0.596, d2 = 0.594 and d3 = 0.596. In addition:
; 9093:5 5:729 0:951 : P3 ¼ 9093:6
P1 ¼
5:896
0:012
P2 ¼
6:087
1:572
9475:0
;
This completes the application of the proposed procedure. So, the power distribution system under consideration is practically stable with respect to the calculated sets X1 and X2 within the time interval I = [0, 65) for every load variation with a dwell-time T = 6 s. The calculated sets X1 and X2 are shown in Fig. 7. In addition, this figure shows the allowable operating range X and the system trajectory for an initial condition x0 2 X1 and for a switching signal with a dwell-time T = 6 s. It can be seen that the trajectory remains confined in X2 for all t 2 I. Also, the set X2 represents a realistic operating region of a power system, where the rotor speed can vary between, approximately, 59.6 Hz and 60.4 Hz, depending on the value of the rotor angle. 7.2. Tests and results for a fifth order model Consider again the power distribution system with synchronous generator presented in Section 2. In this subsection, we consider a fourth order model for describing the synchronous generator. Also, we consider that the generator is equipped with a first order model of an automatic voltage regulator (AVR). The state-space model is [20]:
_ dðtÞ ¼ ws wðtÞ ws ;
ð30Þ
_ ðtÞ ¼ Pref Pe ðtÞ DðxðtÞ ws Þ; 2Hx
ð31Þ
s0q0 E_ 0d ðtÞ ¼ X q X 0q ðII ðtÞ cosðdðtÞÞ IR ðtÞ sinðdðtÞÞÞ E0d ðtÞ;
ð32Þ
where d(t) and x(t) are the generator rotor angle and angular speed, respectively; E0d ðtÞ and E0q ðtÞ are the direct and quadrature axis transient internal voltage of the generator, respectively, and Efd(t) is voltage applied to the field circuit.iTThe state vector x(t) is given h by xðtÞ ¼ dðtÞ xðtÞ E0d ðtÞ E0q ðtÞ Efd ðtÞ . Detailed information regarding this system model can be obtained in [20]. In Eqs. (30)–(34), IR(t) and II(t) are the real and imaginary components, respectively, of the output generator current; Pe(t) is the generated electrical power and jVtj is the absolute value of generator terminal voltage. These variables IR(t), II(t), Pe(t) and jVt(t)j are calculated by the equations
IR ðtÞ ¼ ya V R ðtÞ þ yb V I ðtÞ;
ð35Þ
II ðtÞ ¼ ya V I ðtÞ yb V R ðtÞ;
ð36Þ
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi V 2R ðtÞ þ V 2I ðtÞ;
ð37Þ
jV t ðtÞj ¼
Pe ðtÞ ¼ E0d ðtÞ½II ðtÞ sinðdðtÞÞ þ IR ðtÞ cosðdðtÞÞ þ E0q ðtÞ½II ðtÞ cosðdðtÞÞ IR ðtÞ sinðdðtÞÞ þ ½II ðtÞ sinðdðtÞÞ þ IR ðtÞ cosðdðtÞÞ½II ðtÞ cosðdðtÞÞ IR ðtÞ sinðdðtÞÞ X 0q X 0d ;
ð38Þ
being VR(t) and VI(t) the real and imaginary components of the generator terminal voltage (bus 4); ya and yb the parameters reflecting the topological characteristics of the distribution network, including the resistance (R) and reactance (X) of the line 2–3, the transformers (T1 and T2) reactances and the load (L) impedance (see [28] for more details). Once the parameters ya and yb are function of the operating conditions, we have k = [ya yb]T. The equations for VR(t) and VI(t) are given by
V R ðtÞ V I ðtÞ
¼
N11 ðtÞ N12 ðtÞ N21 ðtÞ N22 ðtÞ
1 "
E0q ðtÞ E0d ðtÞ
# ;
ð39Þ
where
s0d0 E_ 0q ðtÞ ¼ Efd ðtÞ X d X 0d ½II ðtÞ sinðdðtÞÞ þ IR ðtÞ cosðdðtÞÞ E0q ðtÞ; ð33Þ T a E_ fd ðtÞ ¼ Efd ðtÞ þ K a ðV ref jV t ðtÞjÞ;
ð34Þ
N11 ðtÞ ¼ X 0d ½yb ðtÞ sinðdðtÞÞ ya cosðdðtÞÞ þ sinðdðtÞÞ; N12 ðtÞ ¼ X 0d ½ya ðtÞ sinðdðtÞÞ þ yb ðtÞ cosðdðtÞÞ cosðdðtÞÞ; N21 ðtÞ ¼ X 0q ½yb ðtÞ cosðdðtÞÞ þ ya ðtÞ sinðdðtÞÞ cosðdðtÞÞ; N22 ðtÞ ¼ X 0q ½ya ðtÞ cosðdðtÞÞ yb ðtÞ sinðdðtÞÞ sinðdðtÞÞ: By applying some algebraic manipulations to the set of Eqs. (30)–(39) (replacing (39) into Eqs. (35)–(37) and the resulting equations together with (38) into Eqs. (30)–(34)), the power distribution system with synchronous generator can be mathematically
Fig. 7. The sets X1 and X2 representing a realistic operating region of the power distribution system. Also, system trajectory for an initial condition x0 2 X1 and a dwell-time T = 6 s.
R. Kuiava et al. / Electrical Power and Energy Systems 55 (2014) 275–284
283
Fig. 8. Sets X1 and X2 and the system trajectory for an initial condition x0 2 X1 and dwell-time T = 6 s.
described by a model in the state-space form given in the compact form (1), where the vectors xðtÞ 2 R5 and k 2 R2 where previously defined and f : R5 R2 ! R5 is a nonlinear vector function. The numerical values of the model parameters are: ws = 376.99 rad/s, H = 1.5 s, D ¼ 0:2; X 0d ¼ 0:1408 p:u:, X 0q ¼ 0:434 p:u:, Xq = 0.741 p.u., Xd = 0.759 p.u., s0d0 ¼ 4:75 s, s0q0 ¼ 1:50, K a ¼ 30, Ta = 0.01, and Pref = 0.1 p.u. The three load levels L1, L2 and L3 used in last subsection are also adopted here. The system equilibrium point associated to the load level L = L1 is xe1 ¼ ½0:101 1 0:029 1:021 1:0360 . We consider the following allowable operating range of the system:
nh
d w E0d E0q Efd
i0
2 R5 : 0:2618 6 d 6 0:4363; 0:9917 6 w
6 1:0083; 0 6 E0d 6 1:5; 0 6 E0d 6 1:5; 0 6 Efd 6 1:5 : ð40Þ
X :¼
This operating range was defined by numerical simulations in different load conditions but, obviously, this set should be obtained preferable from a practical knowledge of the system. In order to calculate the set X2 via the optimization problem OP, we have to write the set (40) in the form of Xi, i = 1, 2, 3, as defined by (14). Due to space limitations, the explanation of this process is simplified in this paper. At first, notice that (40) is equivalent to a polytope in R5 whose vertices are defined by H :¼ fc1 ; c2 ; . . . ; c25 g ¼ f½h1 h2 h5 0 2 R5 : hi 2 ½xi ; xi , i ¼ 1; . . . ; 5g, where xi and xi are, respectively, the lower and the upper bound of the ith state variable, as defined in (40). A vertex representation of X is defined as the convex hull of c1 ; c2 ; . . . ; c25 . Equivalently, we can define this vertex form of X as (14), for i = 1, 2, 3, by following the same procedure adopted for the example in Section 7.1. It was specified an initial value of a equal , where l is the upper bound of (14). to 0.5. We specify l ¼ 0:99l The LMI optimization problem OP was solved repeatedly with increasing values of a (from its initial value) until the terminating condition given by the unfeasibility of the LMI problem. The final value of a was 0.6. Fig. 8 shows the sets X1 and X2 by considering that E0d , E0q and E0fd are in its equilibrium conditions. This plot corresponds to the cut of the actual estimate of the sets X1 and X2 in the hyperplane defined by the system states d and x. Also, this figure shows the system trajectory for an initial condition x0 2 X1 and for a switching signal with a dwell-time T = 7 s. It can be seen that the trajectory remains confined in X2 for all t 2 I. 8. Conclusion This paper presented a method to assess the practical stability of power distribution systems with synchronous generators sub-
ject to changes in the system equilibrium conditions due to load variations. Our results provide means to evaluate the conditions in which the system trajectories will be confined within a security operating region, as illustrated in the example of Section 7. It is also important to point out that, although the numerical results presented in Section 7 were provided for a load variation sequence with regular time interval of 5 s between two consecutive variations, our theoretical results are valid for any load variation sequence such that (5) is satisfied. So, irregular time intervals between two consecutive load variations (but, satisfying (5)) will certainly provide practical stable trajectories. A significant advance is achieved with this proposed method, given that it enables the assessment of stability not only of the equilibrium conditions (which is usually carried out by linearized analysis), but also of the transition among the different equilibria that result from load variations in power systems with distributed generation. The results show that, using the proposed method, it is possible to guarantee a stable operation for a single DG operating connected to a feeder with varying loads, which is consistent with the usual practice of utilities to allow only one generator to connect to an individual feeder. However, in the context of microgrids, there might be a number of generators connected to the same islanded system, so the extension of these results to a case with more than one generator is among the future directions of this work. In addition, the application of the proposed method for distributed synchronous generators with power factor control is also of interest for a future work.
References [1] de Abreu RVL, Marques FAS, Morn J, Freitas W, da Silva LCP. Impact of distributed synchronous generators on the dynamic performance of electrical power distribution systems. In: Proc of the IEEE PES trans and distr conf and exp. Latin America; 2004. [2] Moura FAM, Camacho JR, Chaves MLR, Guimarpes GC. Dynamics analysis of two synchronous machines interconnected with a distribution network in ATP-EMTP. In: Proc of the 15th int conf on int syst appl to power syst. IEEE; 2009. [3] Tekiner-Mogulkoc H, Coit DW, Felder FA. Electric power system generation expansion plans considering the impact of smart grid technologies. Int J Electr Power Energy Syst 2012;42(1):229–39. http://dx.doi.org/10.1016/ j.ijepes.2012.04.014. [4] Freitas W, Vieira JCM, Morelato A, Silva LCP, Costa VF, Lemos FAB. Comparative analysis between synchronous and induction machines for distributed generation applications. IEEE Trans Power Syst 2006;21(1):301–11. http:// dx.doi.org/10.1109/TPWRS.2005.860931. [5] Khani D, Yazdankhah AS, Kojabadi HM. Impacts of distributed generations on power system transient and voltage stability. Int J Electr Power Energy Syst 2012;43(1):488–500. http://dx.doi.org/10.1016/j.ijepes.2012.06.007.
284
R. Kuiava et al. / Electrical Power and Energy Systems 55 (2014) 275–284
[6] Edwards FV, Dudgeon JW, MacDonald JR, Leithead WE. Dynamic of distribution networks with distributed generation. In: Proc of the IEEE PES general meeting. IEEE; 2000. [7] Kuiava R, Ramos RA, de Oliveira RV, Bretas NG. An analysis of the potential impacts of electromechanical oscillations on the stability and power quality of distributed generation systems. In: Proc of the IEEE PES general meeting. IEEE; 2008. [8] Salim RH, Kuiava R, Ramos RA, Bretas NG. Impact of power factor regulation on small-signal stability of power distribution systems with distributed synchronous generators. Eur Trans Electr Power 2011;21(7):1923–40. http:// dx.doi.org/10.1002/etep.504. [9] Nguyen TB, Pai MA. A sensitivity-based approach for studying stability impact of distributed generation. Int J Electr Power Energy Syst 2008;30(8):442–6. http://dx.doi.org/10.1016/j.ijepes.2008.03.001. [10] Noroozian R, Abedi M, Gharehpetian GB, Hosseini S. Distributed resources and dc distribution system combination for high power quality. Int J Electr Power Energy Syst 2010;32(7):769–81. http://dx.doi.org/10.1016/ j.ijepes.2010.01.013. [11] Salim RH, Ramos RA, Oleskovicz M. Power quality of distributed generation systems as affected by electromechanical oscillations – definitions and possible solutions. IET Gener Transm Distrib 2011;5(11):1114–23. http:// dx.doi.org/10.1049/iet-gtd.2011.0277. [12] Salim RH, Ramos RA, Bretas NG. Analysis of the small signal dynamic performance of synchronous generators under unbalanced operating conditions. In: Proc of the IEEE PES general meeting. IEEE; 2010. [13] Bohorquez VB. Fast varying loads. In: Proc of the 9th int conf electr power quality and utilisation. IEEE; 2007. [14] Arseneau R. The performance of demand meters under varying load conditions. IEEE Trans Power Deliv 1993;8(4):1708–11. http://dx.doi.org/ 10.1109/61.248277. [15] Boulet B. Modeling and control of an electric arc furnace. In: Proc of the 2003 American control conference. IEEE; 2003. [16] Zhai G, Mitchel AN. Generalized practical stability analysis of discontinuous dynamical systems. In: Proceedings of the 42nd IEEE conf on dec and contr; 2003.
[17] Zhai G, Michel AN. Generalized practical stability analysis of discontinuous dynamical systems. Int J Appl Math Comp Sci 2004;14(1):5–12. [18] Xu X, Zhai G. Practical stability and stabilization of hybrid and switched systems. IEEE Trans Autom Contr 2005;50(11):1867–903. [19] Xu X, Zhai G, He S. Stabilizability and practical stabilizability of continuoustime switched systems: a unified view. In: Proceedings of the 2007 IEEE American contr. conf; 2007. [20] Kuiava R. Projeto de controladores para o amortecimento de oscilaç oes em sistemas elétricos com geraç ao distribuida. PhD thesis, EESC/USP, Sao Carlos, Brazil [In Portuguese]. [21] Esfahani MT, Vahidi B. A new stochastic model of electric arc furnace based on hidden markov model: a study of its effects on the power system. IEEE Trans Power Deliv 2012;27(4):1893–901. http://dx.doi.org/10.1109/ TPWRD.2012.2206408. [22] Zheng T, Makram EB. An adaptive arc furnace model. IEEE Trans Power Deliv 2000;15(3):931–9. http://dx.doi.org/10.1109/61.871355. [23] Kundur P. Power system stability and control. McGraw-Hill; 1994. [24] Boyd S, Ghaoui LE, Feron E, Balakrishnam V. Linear matrix inequalities in system and control theory. Soc Indust Appl Math 1994. [25] Kuiava R, Ramos RA, Pota HR, Alberto LFC. Practical stability of continuoustime switched systems without a common equilibria and governed by a timedependent switching signal. In: Proc of the IEEE int conf on contr and automat. IEEE; 2011. [26] Kuiava R, Ramos RA, Pota HR, Alberto LFC. Practical stability of switched systems without a common equilibria and governed by a time-dependent switching signal. Eur J Control 2013;19(3):206–13. http://dx.doi.org/10.1016/ j.ejcon.2012.11.001. [27] Rohr ER, Pereira LFA, Coutinho DF. Robustness analysis of nonlinear systems subject to state feedback linearization. Controle Autom 2009;20(4):482–9. [28] Kuiava R, Ramos RA, Pota HR. A new method to design robust power oscillation dampers for distributed synchronous generation systems. J Dyn Syst Meas Control 2013;135(3):031011–031011–10. http://dx.doi.org/10.1115/ 1.4023225.
Contents lists available at ScienceDirect
Electrical Power and Energy Systems journal homepage: www.elsevier.com/locate/ijepes
Practical stability assessment of distributed synchronous generators under variations in the system equilibrium conditions Roman Kuiava a,⇑, Rodrigo A. Ramos b, Hemanshu R. Pota c, Luis F.C. Alberto b a
Federal University of Parana (UFPR), Department of Electrical Engineering, Rua Cel. F. H. dos Santos, 100, Jardim das Americas, 81531-980 Curitiba, Brazil University of Sao Paulo (USP), Engineering School of Sao Carlos (EESC), Department of Electrical Engineering, Av. Trabalhador Saocarlense, 400, Centro, 13566-590 Sao Carlos, Brazil c University of New South Wales at Australian Defence Force Academy (UNSW@ADFA), School of Information Technology and Electrical Engineering, PO Box 7916, Canberra BC, ACT 2610, Australia b
a r t i c l e
i n f o
Article history: Received 16 October 2012 Received in revised form 11 September 2013 Accepted 24 September 2013
Keywords: Distributed synchronous generator Practical stability Fast varying loads Switched affine system Linear matrix inequality
a b s t r a c t This paper proposes a method to assess the practical stability of power distribution systems with synchronous generators subject to changes in the system equilibrium conditions due to fast varying loads. The concept of practical stability deals with two known state-space regions X1 (which contains all the initial conditions reflecting the perturbations at which the system is subject during its operation) and X2 (which represents the operating security region of the power distribution system) satisfying X1 X2. The practical stability problem and the focus of this paper is to determine under which conditions the system trajectories will be confined into a security region of operation for a certain time interval of interest, as the equilibrium point of the model changes. This study was carried out using a mathematical model of the distribution system with synchronous generators in the form of a switched affine system. This proposed model is capable of describing the system behavior over a certain period within which changes on the equilibrium conditions of the system can occur. Sufficient conditions for the power distribution system with synchronous generators described as a switched affine system to be practically stable with respect to its operating security region X2 are given in the form of matrix inequalities constraints. The results, obtained for the model of a cogeneration plant of 10 MW added to a distribution network constituted by a feeder and six buses, show that the less stringent properties of the concept of practical stability can be very well-suited to the security analysis of power systems subjected to frequent variations in the load level. Ó 2013 Elsevier Ltd. All rights reserved.
1. Introduction The interest in distributed generation (DG) has increased considerably around the world over the years, driven mainly by market deregulation, continuous demand rise, technological developments and concerns on environment impact. In Brazil, for example, the strong economical incentives for ethanol production are likely to boost the number of cogeneration power plants connected to both subtransmission and distribution grids. Located mainly in southeast and center-west areas of Brazil, this sector of the national industry had already considerable experience in the usage of steam turbines, which have been extensively employed within its internal production [1–3]. As a result, several synchronous generators driven by steam turbines have been connected to the system using a byproduct of the sugar cane production to
⇑ Corresponding author. E-mail addresses: [email protected] (R. Kuiava), [email protected] (R.A. Ramos), [email protected] (H.R. Pota), [email protected] (L.F.C. Alberto). 0142-0615/$ - see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.ijepes.2013.09.013
feed the boilers, which constitutes a very sustainable form of electrical energy production. In view of these facts, within the new environment provided by DG, the networks may present new and unforeseen problems that must be appropriately dealt with. Recent papers, for example, present research studies considering distinct scenarios and technical factors that can affect the dynamic performance of power distribution systems with distributed synchronous generators in terms of transient stability [1,4,5] and small-signal stability in both balanced [6–10] and unbalanced operating conditions [11,12]. By means of modal and eigenvalues analysis, papers [8,7,11] also show that the presence of sustained electromechanical oscillations in the response of these systems can have a significant impact on the power quality across the grid, which means that this impact has to be assessed before granting permission to connect a generator to the grid. This paper extends the security analysis mentioned in the previous paragraph by assessing the behavior of the system as the load levels vary intermittently over time, i.e., fast varying loads (FVLs) type, using the concept of practical stability (which has a for-
276
R. Kuiava et al. / Electrical Power and Energy Systems 55 (2014) 275–284
mal definition that will be shown later in the paper). FVLs in distribution systems can cause undesirable disturbances or an increase in energy losses, as discussed in [13]. Different types of FVLs can be found, for instance, in industrial [14,15] and medicine fields [13]. From a practical viewpoint, the power distribution system must operates within a security region, which is determined from the allowable range of the state variables (rotor angle, rotor speed, generator terminal voltage, and others), even when facing fast variations in its load levels. Hence, in a realistic scenario, the problem of interest and the focus of this paper is to determine under which conditions the system trajectories will be confined into a security region of operation for a certain time interval of interest, as the equilibrium point of the model changes, and this feature can be captured adequately by the concept of practical stability. Generally speaking, the concept of practical stability [16–19] deals with two sets X1 Rn and X2 Rn satisfying X1 X2, which are specified for the initial state and the system trajectories, respectively. These two sets are specified in terms of physical limitations of the system variables and they must contain all the equilibrium points in which the system is expected to operate within a certain period of time I. For our purposes, the set X2 is a representation of the operating security region of the power distribution system, while set X1 involves all the initial conditions reflecting the perturbations at which the system is subject during its operation. The sets X1 and X2 can be determined from a practical knowledge of the system operation and/or via numerical simulations analysis. The practical stability requires that if the initial state is in X1, then the system state should stay in X2 for all t 2 I [17]. Hence, unlike the classical stability definition which is based on the existence of X1 for any X2, here both of the sets X1 and X2 are fixed and so they do not vary. Stability with respect to a set, rather than a particular point, is then the basis of the practical stability concept. A mathematical model of the distribution system with synchronous generators in the form of a continuous-time switched affine system is proposed in this paper to be used for practical stability studies. This proposed model is capable of describing the system behavior over a certain period within which changes on the equilibrium conditions of the system can occur. Sufficient conditions for the power distribution system described as a switched affine system to be practically stable with respect to a known operating security region X2 are given in the form of matrix inequalities constraints. The paper is organized as follows: Section 2 presents a traditional arrangement of distributed system with synchronous generation, which serves as motivation for the application of the proposed stability assessment; Section 3 presents the problem formulation and a motivation example; Section 4 presents the proposed modeling approach used to describe the original power distribution system model by means of a switched affine system; Section 5 presents the fundamentals of the practical stability concept; Section 6 presents the proposed systematic method to assess the practical stability; Section 7 presents some tests and results and Section 8 presents the conclusions and final remarks.
2. The study system A typical arrangement of a distributed system with synchronous generation, usually viewed in cogeneration plants in the sugarcane industry of Brazil, is shown in Fig. 1. This network consists of a 132-kV, 60-Hz subtransmission system, which feeds a 33-kV distribution system through one 132/33-kV, D/Yg transformer
Fig. 1. Diagram of the study system.
(T1). The cogeneration plant (G) is represented by a gas-turbine driven synchronous generator injecting 10 into the grid. A local load (L) is connected to the bus 3. This local load is represented by a conventional static load model with constant impedance, which is calculated from the values of the active and reactive power demand level of the bus 3 and the nominal amplitude voltage of this same bus. We also assume this local load L as being a simplified representation of several different loads coupled to the bus 3 and each of them can be active or not during a period of operation of the distribution power system. Hence, the operating conditions of the system are supposed to be divided according to which load it supplies. The detailed description of the mathematical model of this system is presented in Section 7 of tests and results. In this point of the paper, it is only important to point out that, for dynamic stability studies, the power distribution system with synchronous generator under investigation (or even more general power systems) can be mathematically described by a model in the state-space form given by
_ xðtÞ ¼ f ðxðtÞ; kÞ;
ð1Þ
where xðtÞ 2 Rn and k 2 Rm are, respectively, the vector with the state variables (such as, rotor angle and rotor speed) and the vector with the parameters that are dependent of the operating conditions, such as the system loads [20], and f : Rn Rm ! Rn is a nonlinear vector function. 3. The problem formulation Consider again the power distribution system shown in Fig. 1 and a mathematical model of it written in the state-space form (1). An equilibrium point xe of (1) (such that f(xe, k) = 0) is calculated by using the power-flow outputs for a specific operating condition (or load level), i.e., for a fixed value of k. Let us suppose that an operating condition k = k1 (associated to a local load L = L1) corresponds to the equilibrium point xe1 of (1) and k = k2 (in which the local load is changed to L = L2) to the equilibrium xe2 . Hence, for a total of N load conditions, we have N equilibrium points given by xe1 ; xe2 ; . . . ; xeN 2 Rn . These operating conditions can, obviously, vary among them along a period of time. The frequency of occurrence of these load changes can be forecasted in advance and this information can be used to estimate the minimum amount of time in which the system stays in a particular operating condition (and this information will be explored later in the paper by means of the concept of dwell-time for switched systems). Typical loads that demand large amounts of varying power from the network include, for example, electric arc furnaces [15]. Obviously, this type of load may require a more accurate model [15] (different from the static load model with constant impedance adopted in the paper) to analyze its impacts on system dynamics, such as, a hidden Markov-based model [21] or an adaptive model [22], but this level of model complexity is not taken under consideration in this paper. Other common fast
R. Kuiava et al. / Electrical Power and Energy Systems 55 (2014) 275–284
varying loads are rolling mills and spot welders in industrial area. An overview about different types of fast varying loads in distribution systems is presented in [13]. Let us also assume that linearized models of (1) in the vicinity of these equilibrium points are supposed to be proper descriptions of the system with respect to the load variations of interest. These linearized models (Ri, . . . , RN) are given in the state-space form
Ri : Dx_ i ðtÞ ¼ Ai Dxi ðtÞ;
ð2Þ
where Ai 2 Rnn is the ith state matrix, Dxi ðtÞ ¼ xðtÞ xei and i 2 S, where S = 1, 2, . . . , N. So, Ri correspond to the linear approximation of (1) in the vicinity of the equilibrium point xei . Let us consider a period of time I = [t0, tf), where t0 is the initial time and tf is a finite constant, and the sets X1 Rn and X2 Rn , X1 X2. The set X1 includes the initial conditions created by a number of perturbations to which the system may be subject during its operation and the set X2 represents the security region of the system operation. Now, let us say that the load conditions change from Li to Lj (which means that Ri changes to Rj, for some i, j 2 S, i – j, in the time instants t1, t2, . . . , tk, . . . 2 I(t0 < t1 < t2<
Fig. 2. Load variations and the system description as time goes by.
277
frequency oscillation is smaller for the case where faster load variations are considered. Furthermore, for this same case it is possible to notice that the peak values of the oscillations tend to increase as time goes by, and this fact is an indicative that the system trajectory will not be always into the security region of operation. From this motivation example, this paper provides some sufficient conditions guaranteeing that the system trajectories (with the system subject to load variations, as illustrated in Fig. 2) will be confined into the security region of operation X2 for all the period of time I. This is done by representing the set of linearized models (2) as a switched affine system, as presented in next section. 4. Modeling the power distribution system as a switched affine system Let us consider again a power distribution system with synchronous generators modeled in the state-space form (1) and a set of N linearized models Ri, . . . , RN given in the state-space form (2) obtained from the linearization of (1) with respect to N different equilibrium points of this system. Let us also consider a period of time I = [t0, tf). We assume that, for all t 2 [tk, tk+1), where k = 0, 1, 2, . . . and t0, t1, . . . , tk, . . . 2 I, the system behavior is described by the linearized model Rik ; ik 2 S, and at tk+1, the system description changes from Rik to Rikþ1 ; ikþ1 2 S; ikþ1 – ik . As discussed in the previous section, this means that, at tk+1, the load conditions change from Lik to Likþ1 . Before presenting the proposed system representation in the form of a switched affine system, it is important to discuss the main reason of using a switched affine system for stability analysis instead of the set of linearized models (2). This is done in next subsection. 4.1. The main difficulty of analyzing the system stability via linearized models From the Lyapunov stability theory, the linearized models
R1, . . . , RN are stable if there exists a set of positive definite matrices P1, . . . , PN such that ATi Pi þ Pi Ai 0; i ¼ 1; . . . ; N [24]. These conditions guarantee local stability of the equilibrium points xe1 ; . . . ; xeN in the sense of Lyapunov stability theory, but they do not guarantee that the system trajectories will be confined into a certain security region of operation X2 over the period of time I in the occurrence of changes on the equilibrium conditions. One of the main difficulties here of taking account the transitions of the equilibrium operating conditions to the stability study, is the fact that the state vector also changes during an operating condition variation (from Dxi(t) to Dxj(t)), once that the coordinate axis of each linearized model has a different origin. To better clarify this point, consider a hypothetical second order nonlinear dynamical system, where the state vector is x(t) = [u(t) v(t)]T. Let us assume this system has two equilibrium points, xe1 ¼ ½ue1 v e1 T and xe2 ¼ ½ue2 v e2 T . Hence, the state vectors of the linearized models calculated in the vicinity of these equilibrium points are given by Dx1(t) = [Du1(t) Dv1(t)]T and Dx2(t) = [Du2(t) Dv2(t)]T, being Du1 ðtÞ ¼ uðtÞ ue1 , Du2 ðtÞ ¼ uðtÞ ue2 , Dv 1 ðtÞ ¼ v ðtÞ v e1 and Dv 2 ðtÞ ¼ v ðtÞ v e2 . Fig. 4(a) illustrates the coordinate axis (u v) of this nonlinear system, where its origin is at (0, 0), as well as, the coordinate axes (Du1 Dv1) and (Du2 Dv2) of the linearized models with respect to the equilibriums xe1 and xe2 , where their origins are at xe1 and xe2 , respectively. Notice that for all the time interval [t0, t1) the system trajectory is obtained from the linearized model R1. At t = t1, the operating conditions of the system is changed, which means that from t1 and on the system trajectory is obtained from R2. The coordinate axis was also changed (from (Du1 Dv1) to (Du2 Dv2)) at
278
R. Kuiava et al. / Electrical Power and Energy Systems 55 (2014) 275–284
Fig. 3. System trajectory (rotor speed, in Hertz) by considering that the time ellapsing two consecutive load variations is no smaller than 1 s (dashed line) and 5 s (full line).
Fig. 4. (a) System trajectory given by linearized approximations with different coordinate axes and (b) system trajectory given by linearized approximations with a same coordinate axis.
the time instant t1, and this fact increases the complexity of analyzing the system stability. Based on the above discussion, next section provides a possible solution for the pointed problem. 4.2. The proposed solution
ations (or to consecutive equilibrium conditions variations). As previously mentioned, the frequency of occurrence of these load changes can be forecasted in advance and this information is adopted to estimate the dwell-time T. For that, we model the set of affine systems R1 ; . . . ; RN as a switched affine system, as presented in next subsection. 4.3. The proposed model for practical stability studies
The difficulty presented in previous section is overcome in this paper by, firstly, maintaining the coordinate axis of the nonlinear system (1) to the linearized models (2), which is done by transforming each of R1, . . . , RN in an affine system. This transformation process is given by
The proposed mathematical model of the power distribution system with synchronous generators for practical stability studies is given by the switched affine system
Ri : Dx_ i ðtÞ ¼ Ai Dxi ðtÞ )
_ xðtÞ ¼ ArðtÞ xðtÞ þ brðtÞ ;
_ x_ ei ¼ Ai ðxðtÞ xei Þ ) xðtÞ _ xðtÞ ¼ Ai xðtÞ Ai xei )
ð3Þ
b i : xðtÞ _ R ¼ Ai xðtÞ þ bi ; where bi ¼ Ai xei . Notice that the equilibrium point of Ri is xei , while the equilibrium of Ri is the origin Dxi(t) = 0. Let us now consider again the hypothetical second order dynamical system presented in last subsection and its trajectory for a certain initial condition, as illustrated now in Fig. 4(b). In this case, for all the time interval [t0, t1) the system trajectory is obtained from the affine system R1 . At t = t1, the operating conditions of the system is changed, which means that from t1 and on the system trajectory is obtained from R2 . Different from the situation illustrated in Fig. 4(a), the coordinate axis was not changed at t = t1 (in fact, the affine systems R1 and R2 are refered to the same coordinate axis (u v) of the nonlinear system). In Problem 1 formulated in Section 3, we focused on the conditions in which the system trajectory will be confined into the security region of operation X2 for all t 2 I. These conditions are evaluated in this paper in terms of a minimum time T (that will be called dwell-time) ellapsed between two consecutive load vari-
xðt0 Þ ¼ x0 ;
ð4Þ
where r: I ? S is called switching signal. Notice that for all the time in which r(t) = i, i 2 S, the trajectory of (4) is in fact the trajectory of the affine system Ri , whose equilibrium point is xei . Also, the switching signal r(t) is given by
rðtÞ ¼ ik 2 S; 8t 2 ½tk ; tkþ1 Þ; tkþ1 tk P T;
ð5Þ
where t1, t2, . . . , tk, tk+1, . . . 2 I are the instants of time in which the system equilibrium conditions change and T is the minimum time ellapsed between two consecutive load variations. Next subsection validates this proposed model for some numerical simulations. 4.4. Validation of the proposed model In this subsection, the approximation of the power distribution system with synchronous generators model in the form of the switched affine system (4) with r(t) as as defined by (5) is validated by means of some numerical simulations. For that, consider again the study system described in Section 2. A switched affine system was obtained by considering three load conditions L1, L2 and L3, whose numerical values are shown in Section 7. Hence, for each of these load conditions, an affine system was obtained
R. Kuiava et al. / Electrical Power and Energy Systems 55 (2014) 275–284
by following (3), which means that the switched affine system is constituted by three affine systems R1 , R2 and R3 . Fig. 5(a) and (b) shows the system trajectory (rotor speed, in Hertz) obtained from the nonlinear model (1) (gray full line) and the system trajectory obtained from the approximated model in the form of a switched affine system (black dashed line). These figures were obtained by considering the time ellapsing between two consecutive load variations is no smaller than 1 s (Fig. 5(a)) and 5 s (Fig. 5(b)), which means that T in (5) is equal to 1 s for the first case and 5 s for the second case. As it can be seen in Fig. 5(a) and (b), the switched affine system describes adequately (at least, for the operating conditions under consideration, which are detailed in Section 7 of tests and results) the system trajectory over a period of time where changes on the equilibrium conditions occur due to load variations. Discrepancies between the trajectories of the nonlinear and affine system models can be found mainly if the operating conditions (equilibrium points) are far from each other. Hence, from now and on we consider the power distribution system with synchronous generators modeled in the form of the switched affine system (4) with r(t) given by (5). Next section provides some sufficient conditions of practical stability for this switched affine system in terms of a minimum time T that must be ellapse between two consecutive load variations. 5. Sufficient conditions of practical stability of switched affine systems
279
Fig. 6. Illustrative example of the practical stability concept.
It is important to emphasize again that the set X2 is assumed to be the representation of a region within which the safe operation of the power system is guaranteed, while the set X1 includes the initial conditions created by a number of perturbations to which the system may be subject during its operation (these perturbations are usually known in advance as the results of the process of contingency screening). Besides, the dwell-time T of r(t) is a practical information that gives an idea on how often changes on the equilibrium conditions occur in the system as time goes by. In this paper, the sets X1 and X2 are assumed to be given in the form of the ellipsoids
X1 ¼
N \
X1i ;
ð6Þ
i¼1
Consider a power distribution system described mathematically as a switched affine system in the form (4) with r(t) given by (5). In the sequence, the practical stability concept is formalized and Fig. 6 illustrates this concept for a second order dynamical system with the sets X1 and X2 given in the form of ellipsoids. Definition 1. The switched affine system (4) with r(t) given by (5) is considered practically stable with respect to the sets X1 Rn and X2 Rn (X1 X2) in the time interval I = [t0, tf), if x(t0) 2 X1 implies x(t) 2 X2, for all t 2 I.
X2 ¼
N [
X2i ;
ð7Þ
i¼1
where
X1i ¼ fx 2 Rn : ðx xei Þ0 Pi ðx xei Þ þ di 6 ag;
ð8Þ
X2i ¼ fx 2 Rn : ðx xei Þ0 Pi ðx xei Þ þ di 6 1g;
ð9Þ
Fig. 5. System trajectory (rotor speed, in Hertz) of the nonlinear model (1) (gray full line) and the system trajectory of the approximated model in the form of a switched affine system (black dashed line): (a) T = 1 s and (b) T = 5 s.
280
R. Kuiava et al. / Electrical Power and Energy Systems 55 (2014) 275–284
being a a positive scalar satisfying a < 1, Pi 2 Rnn a positive definite matrix, di a positive scalar satisfying di < a; xei ¼ A1 i bi the equilibrium point of Ri and i = 1, . . . , N. Next theorem provides sufficient conditions for practical stability of the switched affine system (4) with r(t) given by (5), considering the sets X1 and X2 given in the form of (6) and (7), respectively. This result was previously proposed by the authors in [25,26], so it is presented here without proof. Theorem 1. Let the sets X1 and X2 be given in the form of (6) and (7), respectively. If there exist a scalar TD > 0 and a number l such that the inequalities
1
pffiffiffiffiffiffiffiffiffi 1=a;
N ðt0 ;tf Þ
0
eAi T D Pj eAi T D lPi
ð11Þ 0
eAi T D Pj Dxeij 0 Dxeij Pj Dxeij þ dj ldi
0;
6. The proposed method for practical stability analysis of power distribution systems The proposed method for practical stability analysis of power distribution systems includes the computation of the set X2 in the form of the ellipsoid (7) from the knowledge of the allowable operating range of the system. This allowable operating range is given by
i ¼ 1; . . . ; ng;
ð13Þ
xi
where xi and are, respectively, the lower and upper bounds of the state variable xi. The set X can be written as
X ¼ X i :¼ fx 2 Rn : aik ðx xei Þ 6 1;
k ¼ 1; . . . ; ne g;
Lemma 1. Consider the state-space region Xi defined as (14). The condition x 2 Xi can be written as
k ¼ 1; . . . ; ne :
ð17Þ
Proof. The matrix inequality (17) can be written as
0 a0ik 2
P i
0 : di 1
ð18Þ
Once the constrain (18) is satisfied, the S-Procedure guarantees that for all x such that
x xe i 1
0
Pi
x xe i 1
0
0 di 1
0 a0ik 2
x xe i 1
0;
ð19Þ
ð15Þ
x xe i 1
0:
ð20Þ
Now, notice that (19) and (21) are the matrix form of (16) and (15), respectively. This completes the proof. h The step-by-step procedure proposed below systematizes the use of Theorem 1 to assess the practical stability of power distribution systems with synchronous generators modeled as a switched affine system in the form (4) with r(t) given by (5). This procedure includes the computation of the set X2 X with its size as large as possible. Step 1: Initialization – (i) compute the equilibrium point xei of Ri ; i ¼ 1; . . . ; N; (ii) calculate the matrix Ai and the vector bi ¼ Ai xei , i = 1, . . . , N; (iii) from a known value of TD, compute the maximum number of switchings N ðt0 ; tf Þ that can occur within the time interval of interest I = [t0, tf); and (iv) from the allowable operating range of the system given by (13), compute the equivalent set Xi, i = 1, . . . , N, as defined by (14). Step 2: Specifying the initial conditions – (i) specify a value for a such that a < 1; (ii) compute the upper bound of the inequality (10); and (iii) specify a value for l such that (10) is satisfied.1 Step 3: Checking the practical stability – (i) build the computational representation of the matrix variable Pi, the scalar variables di and ki, i = 1, . . . , N. Solve the LMI optimization problem (OP) presented below. If a feasible solution is found, stop2; otherwise, return to Step 2 and specify a smaller value of a and/or a different value of l.
OP :
minimize
N X ki ; i¼1
ð14Þ
by a proper choice of the row vectors aik 2 Rn , k = 1, . . . , ne, for some i 2 S, where ne is the number of edges of Xi. Hence, the set X can be written equivalently in N different forms, i.e., X ¼ X 1 ¼ X 2 ¼ ¼ X N . These N different forms will be useful to compute set X2 exploring LMIs. The procedure to compute the equivalent forms of (13) is discussed in Section 7 of tests and results. Next Lemma provides sufficient conditions to ensure that the set X2i ; i 2 S, as defined by (9) is bounded by the state-space region X written in the equivalent form Xi.
2 ðx xei Þ0 a0ik aik ðx xei Þ P 0;
a0ik 0; di þ 1
it is true that
Proof. The proof can be found in [25,26]. h In the above Theorem, the parameter N ðt0 ; t f Þ is the maximum number of switching times that can occur on the time interval I and, obviously, it depends on the lenght of I and the value of TD. ðt f t 0 Þ ðt t Þ If is an integer number, then N ðt ; t Þ ¼ fT D 0 1; otherwise, 0 f TD ðt f t 0 Þ N ðt 0 ; t f Þ ¼ int T D . The use of Theorem 1 to assess the practical stability of power distribution systems with synchronous generators modeled as a switched affine system is described in next section.
Pi
is satisfied for all k = 1, . . . , ne.
ð12Þ
X :¼ f½x1 xn 0 2 Rn : xi 6 xi 6 xi ;
ð16Þ
Therefore, the condition x 2 X2i X i is guaranteed if the inequality
#
are satisfied for all i, j 2 S, i – j, where Dxeij ¼ xei xej , and N ðt 0 ; t f Þ is the maximum number of switchings that can occur on the time interval I, then the switched affine system (4) is practically stable with respect to X1 and X2 in the time interval I, for every r(t) satisfying (5) with T = TD.
ðx xei Þ0 Pi ðx xei Þ þ di 1 6 0:
ð10Þ
A0i P i þ Pi Ai 0; "
Let the set X2i be defined as (9). If x 2 X2i , then
subject to ki > 0;
di < a;
ð21Þ
0 Pi ki I;
ð22Þ
A0i Pi þ P i Ai 0;
ð23Þ
Pi
a0ik 0; di þ 1
ð24Þ
1 Several numerical tests in different examples have shown that proper values of l are close to the upper bound of (10). 2 Although the LMI problem OP must be solved with a pre-specified value of a, it is important to emphasize that the size of the set X1 can be manipulated by means of an iteractive process. For example, its size can be enlarged as possible by solving the LMI problem OP repeatedly with increasing values of a until a terminating condition (the unfeasibility of the LMI problem) has been met.
281
R. Kuiava et al. / Electrical Power and Energy Systems 55 (2014) 275–284
"
0
eAi T D Pj eAi T D lPi
i; j ¼ 1; . . . ; N;
0
eAi T D Pj Dxeij þ dj ldi
Dx0eij Pj Dxeij
i – j;
# 0;
ð25Þ
k ¼ 1; . . . ; ne :
If a feasible solution is found in this optimization problem, then the trajectories of the power distribution system with synchronous generators (described in the form of the switched affine system (4)) will be confined in the security region of operation X2 X for all the time interval I, ever so the initial condition is in X1 and the time ellapsed between two consecutive load variations is no smaller than T = TD. It is important to notice that the size of X2 is being enlarged as possible by means of the minimization of the highest eigenvalue of Pi, i = 1, . . . , N, via the inequality on the right of (22). In addition, the specification of l and a are only needed in order to transform the inequality (25) in an LMI. Next section illustrates the application of the above step-by-step procedure to study the practical stability of the system presented in Section 2. 7. Tests and results
Let us consider an affine representation in the vicinity of the equilibrium points xe1 ; xe2 and xe3 , as calculated by (3). Let us also consider a time interval I = [0, 65). Hence, the switched affine system adopted to study the practical stability of the distributed system with synchronous generator under consideration is described by
_ xðtÞ ¼ ArðtÞ xðtÞ þ brðtÞ ;
xð0Þ ¼ x0 ;
ð28Þ
where r: I = [0, 65) ? S = [1, 2, 3], being
A1 ¼
A2 ¼
A3 ¼
0
376:9911
0:2444
0:0667
0
376:9911
0:2413
0:0667
0
376:9911
0:2358
0:0667
;
b1 ¼
;
b2 ¼
;
b3 ¼
376:9911 0:0927
376:9911 0:0944
376:9911 0:0862
7.1. Tests and results for a second order model
:
0:9917 6 w
6 1:0083g:
Consider the power distribution system with synchronous generator presented in Section 2. In this subsection, we consider a second order model (classical model) for describing the study system, where the state variables are the generator rotor angle and angular speed. The state-space model is:
ð26Þ
ð27Þ
where d(t) and x(t) are the generator rotor angle and the angular speed, respectively, so the state vector is x(t) = [d(t) x(t)]0 . The parameters Eh0q ; G11 ; G12 and iB012 are function of the operating conditions, so k ¼ E0q G11 G12 B12 , being E0q the quadrature axis transient internal voltage of the generator, while G11, G12 and B12 are elements of the reduced bus admittance matrix. This matrix reflects the topological characteristics of the distribution network, including the resistance and reactance of the line 2–3, the reactances of the transformers and the load (L) impedance. Finally, ws is the synchronous speed, H is the inertia constant, Pm is the mechanical power (turbine output) and jV1j is the magnitude of the voltage at bus 1. Detailed information regarding the equations of the presented model and their respective parameters can be obtained in [23]. The numerical values of the other model parameters are: ws = 376.99 rad/s, H = 1.5 s, D = 0.2 and Pm = 0.1 p.u. We considered three different load levels:
L1 ¼ PL1 þ jQ L1 ¼ 3:2 þ j1:0 MVA; L2 ¼ PL2 þ jQ L2 ¼ 5:2 þ j2:0 MVA;
ð29Þ
The ranges of d and w in (29) are expressed in rad and p.u., respectively. In order to calculate the set X2 via the optimization problem OP, we have to write the set (29) in the form of Xi, i = 1, 2, 3, as defined by (14). This was done following the ideas presented in [27]. At first, notice that (29) is equivalent to a polytope in R2 whose vertices are defined by H :¼ {c1, c2, c3, c4}, where
c1 ¼
_ ðtÞ ¼ Pm E02q G11 þ E0 q jV 1 jðG12 cosðdðtÞÞ þ B12 sinðdðtÞÞÞ 2Hx DðxðtÞ ws Þ;
;
It was adopted a TD = 6 s. Notice that t0 = 0 s, tf = 65 s and TD = 6 s, which means that N ðt0 ; tf Þ ¼ 10. We consider the following allowable operating range of the system:
X :¼ f½d w0 2 R2 : 0:2618 6 d 6 0:4363;
_ dðtÞ ¼ ws wðtÞ ws ;
;
c3 ¼
0:2618 0:9917 0:4363 0:9917
;
;
0:2618 c2 ¼ ; 1:0083 0:4363 c4 ¼ : 1:0083
A vertex representation of X is defined as the convex hull of c1, c2, c3 and c4. Equivalently, we can define this vertex form of X as (14), for i = 1, 2, 3. We present in details the construction of the set X1. The sets X2 and X3 can be obtained in a similar manner. In order to calculate the row vector a1k 2 R2 , notice that each inequality in (14), i = 1, defines an hyperplane fx : a1k ðx xe1 Þ 6 1g for which
c1 ; c2 2 fx : a11 ðx xe1 Þ ¼ 1g;
c2 ; c3 2 fx : a12 ðx xe1 Þ ¼ 1g;
c3 ; c4 2 fx : a13 ðx xe1 Þ ¼ 1g;
c4 ; c1 2 fx : a14 ðx xe1 Þ ¼ 1g:
For a given set H, the row vector a1k, k = 1, . . . , 4, can be determined by solving the following set of linear systems:
a11 ðc1 xe1 Þ ¼ 1 a11 ðc2 xe1 Þ ¼ 1 a13 ðc3 xe1 Þ ¼ 1 a13 ðc4 xe1 Þ ¼ 1
; ;
a12 ðc2 xe1 Þ ¼ 1
; a12 ðc3 xe1 Þ ¼ 1 a14 ðc1 xe1 Þ ¼ 1; a14 ðc4 xe1 Þ ¼ 1;
yielding the following row vectors
L3 ¼ PL3 þ jQ L3 ¼ 7:2 þ j3:2 MVA:
a11 ¼ ½ 3:033 0 ;
The step-by-step procedure presented in the last section was used to assess the dynamic behavior of the study power distribution system in terms of practical stability. For the Step 1, we have that the system equilibrium points associated to the load levels L = L1, L = L2 and L = L3 are, respectively, xe1 ¼ ½0:1065 10 ; xe2 ¼ ½0:1152 10 and xe3 ¼ ½0:0828 10 .
This completes the Step 1. For the Step 2, it was specified an initial value of a equal to 0.5. As previously commented, numerical tests have shown that proper values of l are close to its upper , where l is the upper bound of bound. So, we specify l ¼ 0:99l (14).
a13 ¼ ½ 2:714 0 ;
a12 ¼ ½ 0 120:5 ; a14 ¼ ½ 0 120:5 :
282
R. Kuiava et al. / Electrical Power and Energy Systems 55 (2014) 275–284
For the Step 3, the LMI optimization problem OP was solved repeatedly with increasing values of a (from its initial value specified in Step 2) until the terminating condition given by the unfeasibility of the LMI problem. For all the process, we adopted l ¼ 0:99l . The final value of a was 0.67. The final solution of the LMI optimization problem OP was: k1 = 9093.4, k2 = 9475.1, k3 = 9093.6, d1 = 0.596, d2 = 0.594 and d3 = 0.596. In addition:
; 9093:5 5:729 0:951 : P3 ¼ 9093:6
P1 ¼
5:896
0:012
P2 ¼
6:087
1:572
9475:0
;
This completes the application of the proposed procedure. So, the power distribution system under consideration is practically stable with respect to the calculated sets X1 and X2 within the time interval I = [0, 65) for every load variation with a dwell-time T = 6 s. The calculated sets X1 and X2 are shown in Fig. 7. In addition, this figure shows the allowable operating range X and the system trajectory for an initial condition x0 2 X1 and for a switching signal with a dwell-time T = 6 s. It can be seen that the trajectory remains confined in X2 for all t 2 I. Also, the set X2 represents a realistic operating region of a power system, where the rotor speed can vary between, approximately, 59.6 Hz and 60.4 Hz, depending on the value of the rotor angle. 7.2. Tests and results for a fifth order model Consider again the power distribution system with synchronous generator presented in Section 2. In this subsection, we consider a fourth order model for describing the synchronous generator. Also, we consider that the generator is equipped with a first order model of an automatic voltage regulator (AVR). The state-space model is [20]:
_ dðtÞ ¼ ws wðtÞ ws ;
ð30Þ
_ ðtÞ ¼ Pref Pe ðtÞ DðxðtÞ ws Þ; 2Hx
ð31Þ
s0q0 E_ 0d ðtÞ ¼ X q X 0q ðII ðtÞ cosðdðtÞÞ IR ðtÞ sinðdðtÞÞÞ E0d ðtÞ;
ð32Þ
where d(t) and x(t) are the generator rotor angle and angular speed, respectively; E0d ðtÞ and E0q ðtÞ are the direct and quadrature axis transient internal voltage of the generator, respectively, and Efd(t) is voltage applied to the field circuit.iTThe state vector x(t) is given h by xðtÞ ¼ dðtÞ xðtÞ E0d ðtÞ E0q ðtÞ Efd ðtÞ . Detailed information regarding this system model can be obtained in [20]. In Eqs. (30)–(34), IR(t) and II(t) are the real and imaginary components, respectively, of the output generator current; Pe(t) is the generated electrical power and jVtj is the absolute value of generator terminal voltage. These variables IR(t), II(t), Pe(t) and jVt(t)j are calculated by the equations
IR ðtÞ ¼ ya V R ðtÞ þ yb V I ðtÞ;
ð35Þ
II ðtÞ ¼ ya V I ðtÞ yb V R ðtÞ;
ð36Þ
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi V 2R ðtÞ þ V 2I ðtÞ;
ð37Þ
jV t ðtÞj ¼
Pe ðtÞ ¼ E0d ðtÞ½II ðtÞ sinðdðtÞÞ þ IR ðtÞ cosðdðtÞÞ þ E0q ðtÞ½II ðtÞ cosðdðtÞÞ IR ðtÞ sinðdðtÞÞ þ ½II ðtÞ sinðdðtÞÞ þ IR ðtÞ cosðdðtÞÞ½II ðtÞ cosðdðtÞÞ IR ðtÞ sinðdðtÞÞ X 0q X 0d ;
ð38Þ
being VR(t) and VI(t) the real and imaginary components of the generator terminal voltage (bus 4); ya and yb the parameters reflecting the topological characteristics of the distribution network, including the resistance (R) and reactance (X) of the line 2–3, the transformers (T1 and T2) reactances and the load (L) impedance (see [28] for more details). Once the parameters ya and yb are function of the operating conditions, we have k = [ya yb]T. The equations for VR(t) and VI(t) are given by
V R ðtÞ V I ðtÞ
¼
N11 ðtÞ N12 ðtÞ N21 ðtÞ N22 ðtÞ
1 "
E0q ðtÞ E0d ðtÞ
# ;
ð39Þ
where
s0d0 E_ 0q ðtÞ ¼ Efd ðtÞ X d X 0d ½II ðtÞ sinðdðtÞÞ þ IR ðtÞ cosðdðtÞÞ E0q ðtÞ; ð33Þ T a E_ fd ðtÞ ¼ Efd ðtÞ þ K a ðV ref jV t ðtÞjÞ;
ð34Þ
N11 ðtÞ ¼ X 0d ½yb ðtÞ sinðdðtÞÞ ya cosðdðtÞÞ þ sinðdðtÞÞ; N12 ðtÞ ¼ X 0d ½ya ðtÞ sinðdðtÞÞ þ yb ðtÞ cosðdðtÞÞ cosðdðtÞÞ; N21 ðtÞ ¼ X 0q ½yb ðtÞ cosðdðtÞÞ þ ya ðtÞ sinðdðtÞÞ cosðdðtÞÞ; N22 ðtÞ ¼ X 0q ½ya ðtÞ cosðdðtÞÞ yb ðtÞ sinðdðtÞÞ sinðdðtÞÞ: By applying some algebraic manipulations to the set of Eqs. (30)–(39) (replacing (39) into Eqs. (35)–(37) and the resulting equations together with (38) into Eqs. (30)–(34)), the power distribution system with synchronous generator can be mathematically
Fig. 7. The sets X1 and X2 representing a realistic operating region of the power distribution system. Also, system trajectory for an initial condition x0 2 X1 and a dwell-time T = 6 s.
R. Kuiava et al. / Electrical Power and Energy Systems 55 (2014) 275–284
283
Fig. 8. Sets X1 and X2 and the system trajectory for an initial condition x0 2 X1 and dwell-time T = 6 s.
described by a model in the state-space form given in the compact form (1), where the vectors xðtÞ 2 R5 and k 2 R2 where previously defined and f : R5 R2 ! R5 is a nonlinear vector function. The numerical values of the model parameters are: ws = 376.99 rad/s, H = 1.5 s, D ¼ 0:2; X 0d ¼ 0:1408 p:u:, X 0q ¼ 0:434 p:u:, Xq = 0.741 p.u., Xd = 0.759 p.u., s0d0 ¼ 4:75 s, s0q0 ¼ 1:50, K a ¼ 30, Ta = 0.01, and Pref = 0.1 p.u. The three load levels L1, L2 and L3 used in last subsection are also adopted here. The system equilibrium point associated to the load level L = L1 is xe1 ¼ ½0:101 1 0:029 1:021 1:0360 . We consider the following allowable operating range of the system:
nh
d w E0d E0q Efd
i0
2 R5 : 0:2618 6 d 6 0:4363; 0:9917 6 w
6 1:0083; 0 6 E0d 6 1:5; 0 6 E0d 6 1:5; 0 6 Efd 6 1:5 : ð40Þ
X :¼
This operating range was defined by numerical simulations in different load conditions but, obviously, this set should be obtained preferable from a practical knowledge of the system. In order to calculate the set X2 via the optimization problem OP, we have to write the set (40) in the form of Xi, i = 1, 2, 3, as defined by (14). Due to space limitations, the explanation of this process is simplified in this paper. At first, notice that (40) is equivalent to a polytope in R5 whose vertices are defined by H :¼ fc1 ; c2 ; . . . ; c25 g ¼ f½h1 h2 h5 0 2 R5 : hi 2 ½xi ; xi , i ¼ 1; . . . ; 5g, where xi and xi are, respectively, the lower and the upper bound of the ith state variable, as defined in (40). A vertex representation of X is defined as the convex hull of c1 ; c2 ; . . . ; c25 . Equivalently, we can define this vertex form of X as (14), for i = 1, 2, 3, by following the same procedure adopted for the example in Section 7.1. It was specified an initial value of a equal , where l is the upper bound of (14). to 0.5. We specify l ¼ 0:99l The LMI optimization problem OP was solved repeatedly with increasing values of a (from its initial value) until the terminating condition given by the unfeasibility of the LMI problem. The final value of a was 0.6. Fig. 8 shows the sets X1 and X2 by considering that E0d , E0q and E0fd are in its equilibrium conditions. This plot corresponds to the cut of the actual estimate of the sets X1 and X2 in the hyperplane defined by the system states d and x. Also, this figure shows the system trajectory for an initial condition x0 2 X1 and for a switching signal with a dwell-time T = 7 s. It can be seen that the trajectory remains confined in X2 for all t 2 I. 8. Conclusion This paper presented a method to assess the practical stability of power distribution systems with synchronous generators sub-
ject to changes in the system equilibrium conditions due to load variations. Our results provide means to evaluate the conditions in which the system trajectories will be confined within a security operating region, as illustrated in the example of Section 7. It is also important to point out that, although the numerical results presented in Section 7 were provided for a load variation sequence with regular time interval of 5 s between two consecutive variations, our theoretical results are valid for any load variation sequence such that (5) is satisfied. So, irregular time intervals between two consecutive load variations (but, satisfying (5)) will certainly provide practical stable trajectories. A significant advance is achieved with this proposed method, given that it enables the assessment of stability not only of the equilibrium conditions (which is usually carried out by linearized analysis), but also of the transition among the different equilibria that result from load variations in power systems with distributed generation. The results show that, using the proposed method, it is possible to guarantee a stable operation for a single DG operating connected to a feeder with varying loads, which is consistent with the usual practice of utilities to allow only one generator to connect to an individual feeder. However, in the context of microgrids, there might be a number of generators connected to the same islanded system, so the extension of these results to a case with more than one generator is among the future directions of this work. In addition, the application of the proposed method for distributed synchronous generators with power factor control is also of interest for a future work.
References [1] de Abreu RVL, Marques FAS, Morn J, Freitas W, da Silva LCP. Impact of distributed synchronous generators on the dynamic performance of electrical power distribution systems. In: Proc of the IEEE PES trans and distr conf and exp. Latin America; 2004. [2] Moura FAM, Camacho JR, Chaves MLR, Guimarpes GC. Dynamics analysis of two synchronous machines interconnected with a distribution network in ATP-EMTP. In: Proc of the 15th int conf on int syst appl to power syst. IEEE; 2009. [3] Tekiner-Mogulkoc H, Coit DW, Felder FA. Electric power system generation expansion plans considering the impact of smart grid technologies. Int J Electr Power Energy Syst 2012;42(1):229–39. http://dx.doi.org/10.1016/ j.ijepes.2012.04.014. [4] Freitas W, Vieira JCM, Morelato A, Silva LCP, Costa VF, Lemos FAB. Comparative analysis between synchronous and induction machines for distributed generation applications. IEEE Trans Power Syst 2006;21(1):301–11. http:// dx.doi.org/10.1109/TPWRS.2005.860931. [5] Khani D, Yazdankhah AS, Kojabadi HM. Impacts of distributed generations on power system transient and voltage stability. Int J Electr Power Energy Syst 2012;43(1):488–500. http://dx.doi.org/10.1016/j.ijepes.2012.06.007.
284
R. Kuiava et al. / Electrical Power and Energy Systems 55 (2014) 275–284
[6] Edwards FV, Dudgeon JW, MacDonald JR, Leithead WE. Dynamic of distribution networks with distributed generation. In: Proc of the IEEE PES general meeting. IEEE; 2000. [7] Kuiava R, Ramos RA, de Oliveira RV, Bretas NG. An analysis of the potential impacts of electromechanical oscillations on the stability and power quality of distributed generation systems. In: Proc of the IEEE PES general meeting. IEEE; 2008. [8] Salim RH, Kuiava R, Ramos RA, Bretas NG. Impact of power factor regulation on small-signal stability of power distribution systems with distributed synchronous generators. Eur Trans Electr Power 2011;21(7):1923–40. http:// dx.doi.org/10.1002/etep.504. [9] Nguyen TB, Pai MA. A sensitivity-based approach for studying stability impact of distributed generation. Int J Electr Power Energy Syst 2008;30(8):442–6. http://dx.doi.org/10.1016/j.ijepes.2008.03.001. [10] Noroozian R, Abedi M, Gharehpetian GB, Hosseini S. Distributed resources and dc distribution system combination for high power quality. Int J Electr Power Energy Syst 2010;32(7):769–81. http://dx.doi.org/10.1016/ j.ijepes.2010.01.013. [11] Salim RH, Ramos RA, Oleskovicz M. Power quality of distributed generation systems as affected by electromechanical oscillations – definitions and possible solutions. IET Gener Transm Distrib 2011;5(11):1114–23. http:// dx.doi.org/10.1049/iet-gtd.2011.0277. [12] Salim RH, Ramos RA, Bretas NG. Analysis of the small signal dynamic performance of synchronous generators under unbalanced operating conditions. In: Proc of the IEEE PES general meeting. IEEE; 2010. [13] Bohorquez VB. Fast varying loads. In: Proc of the 9th int conf electr power quality and utilisation. IEEE; 2007. [14] Arseneau R. The performance of demand meters under varying load conditions. IEEE Trans Power Deliv 1993;8(4):1708–11. http://dx.doi.org/ 10.1109/61.248277. [15] Boulet B. Modeling and control of an electric arc furnace. In: Proc of the 2003 American control conference. IEEE; 2003. [16] Zhai G, Mitchel AN. Generalized practical stability analysis of discontinuous dynamical systems. In: Proceedings of the 42nd IEEE conf on dec and contr; 2003.
[17] Zhai G, Michel AN. Generalized practical stability analysis of discontinuous dynamical systems. Int J Appl Math Comp Sci 2004;14(1):5–12. [18] Xu X, Zhai G. Practical stability and stabilization of hybrid and switched systems. IEEE Trans Autom Contr 2005;50(11):1867–903. [19] Xu X, Zhai G, He S. Stabilizability and practical stabilizability of continuoustime switched systems: a unified view. In: Proceedings of the 2007 IEEE American contr. conf; 2007. [20] Kuiava R. Projeto de controladores para o amortecimento de oscilaç oes em sistemas elétricos com geraç ao distribuida. PhD thesis, EESC/USP, Sao Carlos, Brazil [In Portuguese]. [21] Esfahani MT, Vahidi B. A new stochastic model of electric arc furnace based on hidden markov model: a study of its effects on the power system. IEEE Trans Power Deliv 2012;27(4):1893–901. http://dx.doi.org/10.1109/ TPWRD.2012.2206408. [22] Zheng T, Makram EB. An adaptive arc furnace model. IEEE Trans Power Deliv 2000;15(3):931–9. http://dx.doi.org/10.1109/61.871355. [23] Kundur P. Power system stability and control. McGraw-Hill; 1994. [24] Boyd S, Ghaoui LE, Feron E, Balakrishnam V. Linear matrix inequalities in system and control theory. Soc Indust Appl Math 1994. [25] Kuiava R, Ramos RA, Pota HR, Alberto LFC. Practical stability of continuoustime switched systems without a common equilibria and governed by a timedependent switching signal. In: Proc of the IEEE int conf on contr and automat. IEEE; 2011. [26] Kuiava R, Ramos RA, Pota HR, Alberto LFC. Practical stability of switched systems without a common equilibria and governed by a time-dependent switching signal. Eur J Control 2013;19(3):206–13. http://dx.doi.org/10.1016/ j.ejcon.2012.11.001. [27] Rohr ER, Pereira LFA, Coutinho DF. Robustness analysis of nonlinear systems subject to state feedback linearization. Controle Autom 2009;20(4):482–9. [28] Kuiava R, Ramos RA, Pota HR. A new method to design robust power oscillation dampers for distributed synchronous generation systems. J Dyn Syst Meas Control 2013;135(3):031011–031011–10. http://dx.doi.org/10.1115/ 1.4023225.