Stability Region Estimation Under Low Voltage Ride Through Constraints using Sum of Squares
Chetan Mishra, James S. Thorp, Virgilio A. Centeno, Anamitra Pal
IIEEE Copyright Notice © 2017 IEEE. Personal use of this material is permitted. Permission from IEEE must be obtained for all other uses, in any current or future media, including reprinting/republishing this material for advertising or promotional purposes, creating new collective works, for resale or redistribution to servers or lists, or reuse of any copyrighted component of this work in other works.
Accepted to be Published in: Proceedings of the North American Power Symposium, September 17-19, 2017, Morgantown, WV, USA tability Region Estimation Under Low Voltage RideThrough Constraints using Sum of Squares
Chetan Mishra, James S. Thorp, Virgilio A. Centeno
Bradley Dept. of Electrical and Computer Engineering,Virginia Polytechnic Institute and State University,Blacksburg-24061, Virginia, USAEmail: {chetan31, jsthorp, virgilio}@vt.edu
Anamitra Pal
School of Electrical, Computer, and Energy Engineering,Arizona State University,Tempe-85287, Arizona, USAEmail: [email protected]
Abstract — The increasing penetration of inverter basedrenewable generation (RG) in the form of solar photo-voltaic (PV)or wind has introduced numerous operational challenges anduncertainties. According to the standards [1], [2], these generatorsare made to trip offline if their operating requirements are notmet. In an RG-rich system, this might alter the system dynamicsand/or cause shifting of the equilibrium points to the extent that acascaded tripping scenario is manifested. The present workattempts at avoiding such scenarios by estimating the constrainedstability region (CSR) inside which the system must operate usingmaximal level set of a Lyapunov function estimated through sumof squares (SOS) technique. A time-independent conservativeapproximation of the LVRT constraint is initially derived for aclassical model of the power system. The proposed approach iseventually validated by evaluating the stability of a 3 machine testsystem with trip-able RG.
Keywords— Constrained stability region (CSR), direct methods,low voltage ride through (LVRT), sum of squares (SOS).
I. I
NTRODUCTION
Renewable generation (RG) poses numerous reliabilitychallenges to successful power system operation. Besides theuncertainty in output, which can be dealt with up to a certainextent with forecasting, it is the response of inverter protectionto system disturbances that has emerged as a major threat tosystem stability [3]. Traditionally inverter based generation wasnot expected to support the grid during disturbances and thuswas tripped offline when such a scenario arose. However, withincrease in RG penetration, this tripping could translate to losinglarge amounts of generation, more than what the system mightbe capable of sustaining, at the time of disturbance. A majordevelopment that occurred to safeguard against this problem wasthe development of inverter ride-through features which, as thename suggests, made the inverters stay online (ride through)during disturbances. The ride through standards would clearlydefine operating requirements (in terms of voltage andfrequency) for staying online [1], [2].Direct methods are one of the most popular methods used inindustry and academia alike for analyzing transient stability ofpower systems [4]. However, direct methods are traditionallyfocused on capturing the mode that resulted in loss ofsynchronism of synchronous generators . Furthermore, theirsuccess heavily depends on the numerous system relatedassumptions that are made [5]. Loss of synchronism ischaracterized by a system trajectory leaving the stability region (SR) of the stable equilibrium point (SEP) of interest. Anotherway a system can leave the SR is when it transitions to anothersystem having a different SEP. The latter can happen in a self-protecting system, such as a power system that is operating withconsiderable RG penetration as these generators are constrained to trip if the ride-through standards are not met during thedisturbance. While at first it may not occur to be a threat, thisconstrained tripping in an RG-rich system could initiate apotentially cascaded tripping scenario that can have disastrousconsequences. Making transient stability assessments for suchconstrained systems is also extremely difficult as along with lossof synchronism of synchronous generators, the direct methodsmust also estimate tripping of renewable generators.The only previous research on this topic [6] converted suchconstrained systems into unconstrained systems, paving the wayfor energy function based methods like boundary controllingunstable (BCU) equilibrium point method [7] to becomeapplicable. Also proposed was a time independent formulationof the low voltage ride through (LVRT) constraint. In this paper,we present an alternate formulation for solving this problem; onethat is based on a Lyapunov function approach [8]. We willdemonstrate the application of sum of squares (SOS) techniquesin estimating the CSR for LVRT constrained systems.The rest of the paper is structured as follows. The idea ofCSR along with the modeling of LVRT constraint for a classicalpower system model is presented in Section II. The formulationof estimation of constrained region of attraction (ROA) in theform of SOS optimization is proposed in Section III. Lastly,effectiveness of the proposed methodology is demonstrated inSection IV and validated against time-domain simulations.II. L OW V OLTAGE R IDE T HROUGH (LVRT) C
ONSTRAINT
A. Constrained systems and stability region
Equality constraints when present, force the systems toevolve over a manifold. An example is the power system whichis constrained to satisfy balance of nodal injections. Systemswith equality constraints are usually modeled using differentialalgebraic equations. It is also common to find systems havinginequality constraints. Such constraints can result from physicallimitations of the system under study, requirement to meet someoperating standard, or due to a preference of the system designer.An example is the condition that no RG should trip during a lowvoltage situation, which can be mathematically expressed as, ( )| > ( ) (1)
In (1), ( ) is the voltage at the point of interconnection forthe given inverter based generation and ( ) is the voltagereference value from the corresponding LVRT curve. Whendealing with constrained systems, it is important to discuss theidea of feasible region of state space which is a collection of allpoints that do not violate any of the constraints. In order to doso, consider a generalized autonomous constrained system withstate vector governed by the following equations: = ( ) (2) (3) (4)
It should be noted that (2)-(4) are set-up such that ( ) = 0 , ( ) ( ) = 0 i.e. if the systeminitializes on the manifold defined by ( ) = 0, then it can onlyevolve over that manifold. On the other hand, the inequalityconstraint, ( ) 0 , needs to be satisfied, but does notinfluence the system trajectories. This implies that theequilibrium points for the unconstrained system that lie insidethe feasibility region make up the equilibrium points of theconstrained system. Given a stable equilibrium point for thissystem, the CSR is defined as a collection of all those points instate space from which the emerging trajectories remain totallyinside the feasibility region and eventually converge to .We estimate the CSR using the largest invariant setcontaining which does not intersect the infeasible part of statespace. A Lyapunov function ( ) [8], is estimated such that itscorresponding level set { | ( ) 1} serves as the invariant setto be found. We are only concerned with the portion of the statespace = { | ( ) = 0} where the trajectories are constrainedto evolve. Assuming = 0 , the Lyapunov conditions to befulfilled are [9]:1. ( ) > 0 \{0} and (0) = 0 ( ) = ( ) < 0 \{0} { | ( ) 1} ( ) 0 { | ( ) 1} B. LVRT constraint modeling for the classical power system
LVRT curve is a time-dependent lower limit of the voltageat the renewable generator bus. Now, having time-independentconstraints is much simpler to deal with. We will be taking aconservative approach in this study by approximating LVRTcurve as a constant value equal to its maximum. While thisconstraint might appear to be very limiting, for systems withstatic load models the voltage recovery at the instance of faultclearing is almost instantaneous. Therefore, this will not be aconservative approximation for such systems. Mathematically,for the inverter based generation this can be written as, ( ) > ( ) (5) We will first use the network reduced model of the classicalpower system. In accordance with the common modelingpractice, inverter based generation will be modeled as anegative real load. Thus, inverter dynamics will not bemodeled. Since we will be using a sum of squares (SOS) based approach which requires polynomial systems, the first task is toconvert the system into a polynomial system through variabletransformation. We will use the transformation described inTable I [10], where superscript indicates the SEP of interestand refers to the synchronous machine.
TABLE I. V
ARIABLE T RANSFORMATION [10]
New Variable Function of Original States (0) =sin( )1 cos( )
The classical power system model for in a singlemachine reference frame for uniform damping is shown below. == cos : + sin + (cos + : sin = , == 1,2 … ( 1) Using the variable transformation described in Table I, thesystem becomes a constrained system with the origin being therelevant SEP as shown in (9)-(12). = +++++ (9) = (1 ) × (10) = × (11) (12) = 1,2 … 1
Since the network buses are reduced, the expression forvoltages at the inverter buses needs to be derived as an explicitfunction of states ( , ). The LVRT constraint also needs to bea polynomial function of the transformed states. Using thenetwork reduction process described in [11], we get (13)where is the injected current vector, is the internal EMFvector for generators (excluding inverter based generatorswhich are modeled as negative loads), is the voltage vectorfor network buses (to be reduced), and ’s are the componentsof the admittance matrix. Re-arranging (13), we get, [ ] × =[ ] [ ] | | cos +cos + (14)ubstituting the transformed states in (14), we get (15). [ ] × = [ ] × × 1 +1 + (15) In (15), is a constant matrix consisting of terms from ,and . Now, as the bus voltage magnitudes involve asquare root term they are not polynomial functions. Therefore,we use the LVRT constraint in terms of the square of busvoltage magnitudes which then become polynomial functionsof the transformed variables. Let inverter based generatorswith respective LVRT curves be present. Then, the LVRTconstraint for each will be written in terms of lower limit on thevoltage magnitude squared as shown below. , … = ( )0 [1, ] (16)
III. C
ONSTRAINED
ROA E
STIMATION U SING S UM OF S QUARES
Sum of squares (SOS) based SR estimation has shownpromising results in comparison to pre-existing Lyapunov basedtechniques for power systems [10]. Owing to the ease withwhich complex systems can be handled in the SOS approach, wewill be using this technique for estimating our CSR. Beforedoing so, we provide a brief background on SOS.
A. Relevant concepts in sum of squares (SOS)
A function ( ) is called a sum of squares if it can be writtenas ( ) = ( ) . For a polynomial ( ) , the SOS check isequivalent to finding a symmetric positive semidefinite matrix(gram matrix) such that ( ) = ( ) ( ) where ( ) areall the monomials of degree half the degree of ( ) [12]. Ithas been shown that belongs to an affine subspace ofsymmetric matrices and is written as = + (17) where ( ) ( ) = ( ) and are basis matrices fulfillingthe condition = 0 . Thus, SOS problem can be convertedto a linear matrix inequality (LMI) problem of finding . + 0
Now, Positivestellensatz (P-Satz) Theorem [13] states thatgiven polynomials { ( ), ( ) … . ( )} , { ( ), ( ) … . ( )} , { ( ), ( ) … . ( )} thefollowing are equivalent: ( ) 0 , ( ) 0 … . ( ) 0( ) 0, ( ) 0 … . ( ) 0( ) = 0, ( ) = 0 … . ( ) = 0= (18) Polynomials ( ), ( ) … . ( ) ,( ), ( ) … . ( ) ,( ), ( ) … . ( )+ + = 0 (19)
This theorem is a powerful tool for modeling the constraintsas SOS constraints. For the rest of the paper, we will be usingand not to represent the state vector.
B. Local estimate of Lyapunov function
To get an initial Lyapunov function and its correspondingROA estimate, we use an iteration of expanding D algorithm[14]. We start with a known polynomial ( ) and a semi-algebraic set = { | ( ) } and expand it such that , unknown function ( ) is decreasing < 0, > 0 ,and ( ) 0 . Here it should be kept in mind that since thesystem dynamics are constrained to the manifold ={ | ( ) = 0} , the constraints are only required to be fulfilledon it. In order to apply P-Satz theorem, we represent theconstraints as set emptiness conditions as shown below. max. Constr: ( ) cannot be : P-Satz: ( , ( )) +( ( )) + ( ( )) = 0
Constr: ( ) cannot be : P-Satz: ( , ( )) +( ( )) + ( ( )) = 0( ) =
Constr: ( ) cannot be < 0 : P-Satz: ( ( ), ( )) +( ( )) + ( ( )) = 0[1, ]
Due to the presence of term in the first twoconstraints, those are not semi-algebraic. This is converted to ( ) 0 where ( ) = , is a small number [14]. Themodified SOS constraint conditions are given by, max, , , , , , ,( ) +( )( ) +[1, ] The functions , & in formulation according to the P-Satz theorem are not unique. Also, there can be terms which area multiplication of multiple unknown functions. For instance, inthe above expression, besides an unknown , there are otherunknown multiplier functions ( and ). These cannot besolved through SOS programming as the problem cannot beconverted to a corresponding LMI. Therefore, the degrees needto be chosen carefully to make it possible to find a feasiblesolution [14]. For example, for a known function , if theonstraints are is SOS and is SOS, first constraint canonly be an SOS if its degree is even. So, if has an odd degree,degree of > degree of . As such, some of the rules wefollow when choosing the degrees for multiplier functions are:a. To limit the size of the problem, each term in theconstraint expression is chosen to have the same/similardegree with the overall degree of the known positiveterms degrees of other terms.b. Degrees of SOS multipliers are even by definition andso each constraint’s degree is even as well.c. Since is 0 at = 0 , does not have constant terms.For further reading on simplifications, please see [12]. C. Iterative ROA estimate using expanding interior algorithm
From Section III.B, we get an initial Lyapunov function .One way of estimating ROA using a given Lyapunov functionis through its level sets. Ref. [14] states that if be adomain containing equilibrium = 0 of the system = ( ) with ( ) being the corresponding Lyapunov function definedin , then any region = { | ( ) } is a positiveinvariant region contained in the equilibrium’s ROA. Thus, fora given Lyapunov function, its largest level set contained insidethe region gives an estimate of its ROA. Since our system alsohas inequality constraints, another condition to be fulfilled bythe level set is that it stays within the feasible region.We use the expanding interior algorithm [14] to convert theproblem of finding the biggest level set for an unknownLyapunov function into an SOS optimization problem. The ideabeing that a semi-algebraic set contained inside a given localestimate for ROA can be expanded. During the expansionprocess, the Lyapunov function as well as its associated ROAestimate are allowed to change as long as they totally contain .Now, this makes the problem heavily dependent on the shape ofchosen. An alternative to this dependency was proposed in[10] where an outer iteration loop was added to change intothe obtained Lyapunov level set at the end of the inner expansionprocess. The optimization problem with set emptiness ( ) constraints after incorporating all the previously discussedconditions can be written as: Given = { | ( ) } withradius and an unknown Lyapunov function with theassociated SR estimate given by { | ( ) 1 }, .. .( ) cannot be for : { ( )0, ( ) = 0, 0} = (20)Is contained inside the SR estimate ( ( )1 ) : { ( ) , ( ) = 0, ( ) 1, ( )1} = (21)Inside the SR estimate, ( ) strictly decreasesalong all trajectories: ( ) 1, ( ) 0, ( ) = 0, 0 = (22)Inside the SR estimate, ( ) is never < 0 : ( ) 1, ( ) 0, ( ) 0, ( ) =0 =[1, ] (23) Using P-Satz theorem as well as the definitions for Cones,Multiplicative Monoids and Identity generator, we can convertthese constraints into the following equivalent conditions: + + = 0 (24) + ( ) + ( 1)+ ( )( 1) ++ ( 1) = 0 (25) + (1 ) + + (1 ) ++ = 0 (26) + (1 ) (1 )+ + = 0 (27)
Here, we deliberately leave out the higher degree terms in allthe constraints (for example in the first constraint) to makeLMI formulation possible as discussed before. Since a productof SOS functions is also SOS, in the first and third constraints,assuming = = 1 and = 0, = × , = ×, = × , = × , = × , = × , =× , we will eliminate and . Similarly in the secondconstraint, assuming = 1, = × ( 1), = = 0 ,we can take ( 1) common from each term. Finally, for thefourth constraint we assume = = 0, = × and take common. Taken together, we get, = = (28) ( ) ( 1) = = (29) (1 ) = = (30) (1 ) + = = (31)
Here, it can be seen that some of the terms have 2 unknownfunctions multiplied (for example × in third constraint)which makes the problem unsolvable by semi-definiteprogramming. To fix this problem, is maximized for a givenfixed using the algorithm summarized below. Expanding Interior (fixed p, maximize ) i. Input , ( ) . Initialize i = 0, ( ) , ( ) ii. Substitute = ( ) , = ( ) in constraints in Eqn.28-31, line search over starting from ( ) till the SOSproblem is infeasible. Let the final obtained values be , , , . Set ( ) = iii. Substitute = , = , = in constraintsin Eqn 28-31, line search over starting from ( ) till the SOS problem is infeasible. Let the finalobtained values be , . Update ( ) =, ( ) = , ( ) = .iv. If ( ) ( ) > , = + 1 , Goto ii.v. Stop. As mentioned before, since the technique is dependent on thechoice of the function , a workaround is to update = at theend of each expanding interior loop [10]. This gives thefollowing overall algorithm. verall Algorithm i. Initialize = 0, ( ) = , ( ) = 0 ii. Run Expanding Interior with ( ) , ( ) as inputs toget .iii. If ( ) then ( ) = , ( ) = 1 , = + 1 ,Goto ii.iv. Stop IV. R
ESULTS AND D ISCUSSION
We will demonstrate the implementation as well as theeffectiveness of the proposed technique on a standard 3 machinesystem [15] with RG connected as shown in Fig. 1. A small PVmodeled as negative real load ( 0.4 + 0 ) p.u. is added to bus1 with the generator’s output at that bus reduced accordingly.IEEE 1547 standard LVRT curve is used for the PV which hasa maximum value of 0.85 p.u. (Fig. 1) Thus, the feasibilityregion is given by where is the voltage at bus 1.The post-fault system with the PV connected has the SEP at (0.3165,0.3451,0,0) resulting in the transformed statevariables given in Table II. The chosen degrees for the multiplierfunctions are given in Table III.
Figure 1a. Three machine system with PV connected at bus 1;Figure 1b. IEEE 1547 standard LVRT curve
The obtained estimated Lyapunov function for theconstrained system is given by, ( ) = 0.019 × 3.7 × 10 × × + 0.025 ×× 0.019 × × 0.018 × × +0.012 × × + 3.3 × 10 × + 2.4 ×10 × × 4.4 × 10 × × + 0.021 ×× + 4.8 × 10 × × + 1.0 × + 0.21 × × 0.14 × × 4.0 × 10 ×× + 0.65 × 0.018 × × 0.19 ×× + 0.21 × + 0.74 × 0.024 × ×+ 0.7 × 0.73 ×
TABLE II. P
ROPOSED V ARIABLE T RANSFORMATION
New Variable In terms of Original System State sin( 0.3165)1 cos( 0.3165)sin( 0.3451)1 cos( 0.3451)
TABLE III. M
ULTIPLIER F UNCTION D EGREES
Function Degree Function Degree , ,
22 00 , , The estimated SRs { | ( ) 1} in terms of original statesfor the given constrained system vs. the same system with noLVRT constraint are shown for = 2 in Fig. 2. It can benoticed that since the feasibility region is more binding in thehorizontal direction, the estimate hits that feasibility boundaryfirst and stops expanding in that direction.
Figure 2. Constrained vs Unconstrained ROA
In order to validate the estimated constrained ROA, voltageat bus 1 was plotted for post fault trajectories starting fromrandom points within it. As seen in Fig. 3 below, the voltagesfor all the trajectories never fall below 0.85 p.u. Also, theprojection of trajectories emerging from the SR boundary on theangle plane shows that they are stable as well (Fig. 4). Now, inorder to analyze the level of conservativeness of the proposedapproach, we obtain the true CSR through time domainsimulation using the maximum of LVRT curve. As discussedbefore, if the trajectory emerging from a point in state spaceintersects the infeasible region or does not return to the SEP, itcannot belong to the CSR. It can be observed from Fig. 5 thatwhile the estimate covers a large portion of the one obtainedthrough time domain simulation, since the SEP was closer to theright side of the feasibility boundary, when expanding the ROAestimate with SEP as the center, it stops as soon as it intersectsthe right side boundary. The proposed technique could yieldconservative results for systems having slower voltage recovery,s it could be difficult to maintain the immediate post faultsystem voltage above the maximum value of LVRT curve.However, systems with static load models can be dealt witheffectively using the proposed approach.
Figure 3. Voltage along trajectories emerging from the estimatedstability boundaryFigure 4. Projection of trajectories from the estimated stabilityboundary on the angle planeFigure 5. Estimated constrained ROA (brown) vs actual (blue)
V. C
ONCLUSIONS AND F UTURE W ORK
In this paper, estimation of stability region (SR) under thelow voltage ride through (LVRT) constraint is tackled using aLyapunov-based approach through SOS programming. It wasseen that SOS serves as a powerful tool to deal with the analysis of complex dynamical systems. A conservative formulation ofLVRT constraint followed by the derivation of SOS constraintswas proposed for a classical power system model. It was seenfrom the results obtained for a 3 machine test system that theproposed approach gave satisfactory estimates of the CSR.However, there is still scope for improvement. One majorchallenge is that due to computational limitations, it is extremelydifficult to use the proposed approach in large scale systems. Inthis regard, the idea of decomposition using a vector Lyapunovfunction (found to be applicable to large interconnected systemsin [16]) can be explored. With regards to the conservativeness inour approach, modeling of LVRT constraints similar to [6] couldbe tried. Lastly, since the final estimate depended heavily on therelative location of the SEP and the feasibility boundary, in theexpanding interior algorithm, shifting of the center of theexpanding region could be attempted.R
EFERENCES[1] “IEEE Draft Standard for Interconnecting Distributed Resources withElectric Power Systems - Amendment 1,”
IEEE P1547aD3 Dec. 2013
Ecofys
Direct Methods for Stability Analysis of Electric PowerSystems: Theoretical Foundation, BCU Methodologies, andApplications , 1 edition. Wiley, 2011.[5] L. F. C. Alberto and H. D. Chiang, “Towards development ofgeneralized energy functions for electric power systems,” in , 2012, pp. 1–6.[6] A. P. Sohn, A. L. Abrantes, L. F. C. Alberto, and H. D. Chiang,“Stability region of a wind power system under low-voltage ride-through constraint,” in , 2016, pp. 1–5.[7] H.-D. Chiang, F. F. Wu, and P. P. Varaiya, “A BCU method for directanalysis of power system transient stability,”
IEEE Trans. Power Syst. ,vol. 9, no. 3, pp. 1194–1208, Aug. 1994.[8] A. M. Lyapunov,
General Problem of the Stability Of Motion
Washington, DC: CRC Press, 1992.[9] A. Papachristodoulou and S. Prajna, “A tutorial on sum of squarestechniques for systems analysis,” in
Proceedings of the 2005, AmericanControl Conference, 2005. , 2005, pp. 2686–2700 vol. 4.[10] M. Anghel, F. Milano, and A. Papachristodoulou, “AlgorithmicConstruction of Lyapunov Functions for Power System StabilityAnalysis,”
IEEE Trans. Circuits Syst. Regul. Pap. , vol. 60, no. 9, pp.2533–2546, Sep. 2013.[11] M. A. Pai,
Energy Function Analysis for Power System Stability ,Softcover reprint of the original 1st ed. 1989 edition. Springer, 2013.[12] P. A. Parrilo, “Structured semidefinite programs and semialgebraicgeometry methods in robustness and optimization,” phd, CaliforniaInstitute of Technology, 2000.[13] G. Stengle, “A Nullstellensatz and a Positivstellensatz in SemialgebraicGeometry.,”
Math. Ann. , vol. 207, pp. 87–98, 1974.[14] Zachary William Jarvis-Wloszek, “Lyapunov Based Analysis andController Synthesis for Polynomial Systems using Sum-of-SquaresOptimization,” University of California, Berkley, 2013.[15] T. Athay, R. Podmore, and S. Virmani, “A Practical Method for theDirect Analysis of Transient Stability,”
IEEE Trans. Power Appar.Syst. , vol. PAS-98, no. 2, pp. 573–584, Mar. 1979.[16] J. Anderson and A. Papachristodoulou, “A Decomposition Techniquefor Nonlinear Dynamical System Analysis,”