Applying a formula for generator redispatch to damp interarea oscillations using synchrophasors
tto appear in IEEE Transactions on Power Systems, accepted September 2015
Applying a formula for generator redispatch todamp interarea oscillations using synchrophasors
Sarai Mendoza–Armenta,
Member, IEEE
Ian Dobson,
Fellow, IEEE
Abstract —If an interarea oscillatory mode has insufficientdamping, generator redispatch can be used to improve itsdamping. We explain and apply a new analytic formula forthe modal sensitivity to rank the best pairs of generators toredispatch. The formula requires some dynamic power systemdata and we show how to obtain that data from synchrophasormeasurements. The application of the formula to damp interareamodes is explained and illustrated with interarea modes of theNew England 10-generator power system.
Index Terms —Power system dynamic stability, phasor mea-surement units, power system control.
I. I
NTRODUCTION
Power transmission systems have multiple electromechani-cal oscillatory modes in which power system areas can swingagainst each other. In large grids, these interarea oscillationstypically have low frequency in the range 0.1 to 1.0 Hz,and can appear for large or unusual power transfers. Poorlydamped or negatively damped oscillations can become morefrequent as power systems experience greater variability ofloading conditions and can lead to equipment damage, mal-function or blackouts. Practical rules for power system securityoften require sufficient damping of oscillatory modes [1], [2],such as damping ratio of at least 5%, and power transfers ontie lines are sometimes limited by oscillations [1], [3], [4], [5].There are several approaches to maintaining sufficientdamping of oscillatory modes, including limiting power trans-fers [6], installing closed loop controls [1], and the approachof this paper, which is to take operator actions such as redis-patching generation [6], [7], [8]. It is now feasible to monitormodal damping and frequency online from synchrophasors(also known as PMUs) [9], [10]. Suppose that a mode withinsufficient damping ratio is detected. Then what actionsshould be taken to restore the mode damping?This paper calculates the best generator pairs to redispatchto maintain the mode damping by combining synchrophasorand state estimator measurements with a new analytic formulafor the sensitivity of the mode eigenvalue with respect togenerator redispatch. This formula, previously thought to beunattainable, is derived with a combination of new and oldmethods in [11]. The length of the derivation (more than 8pages) precludes its presentation here. In this paper we state,explain, and demonstrate the application of the new formulaand show how the terms of the formula could be obtainedfrom power system measurements. In particular, we proposeusing synchrophasors to measure the terms of the formula that
The authors are with ECpE dept., Iowa State University, Ames IA USA;[email protected]. We gratefully acknowledge support in part from NSFgrant CPS-1135825 and Arend J. and Velma V. Sandbulte professorship funds. depend on dynamics and using the state estimator to measurethe terms of the formula that depend on statics. In applyinga first order sensitivity formula, we assume that the powersystem ambient or transient behavior is dominated by thelinearized dynamics associated with an asymptotically stableoperating equilibrium.Changes in generator dispatch change the oscillation damp-ing by exploiting nonlinearity of the power system: changingthe dispatch changes the operating equilibrium and hence thelinearization of the power system about that equilibrium thatdetermines the oscillatory modes and their damping. This openloop approach that applies an operator action after too littledamping is detected can be contrasted with an approach thatdesigns closed loop controls to damp the oscillations preven-tively. The closed loop control design chooses control gainsthat appear explicitly in the power system Jacobian, whereasgenerator redispatch changes the Jacobian indirectly by chang-ing the operating point at which the Jacobian is evaluated.We now review previous works using generator redispatch todamp oscillations; these have considered heuristics, brute forcecomputations, and formulas that are difficult to implementfrom measurements. These approaches have established thatgenerator redispatch can damp oscillatory modes.Fischer and Erlich pioneered heuristics for the redispatchin terms of the mode shapes for some simple grid structuresand for the European grid [8], [12]. Their heuristics seempromising for insights, elaborations and validation, especiallysince there has also been progress in determining the modeshape from measurements [13], [14], [15], [16].There are also previous approaches that require a dynamicgrid model. The effective generator redispatches can be deter-mined by repetitive computation of eigenvalues of a dynamicpower grid model to give numerical sensitivities [5], [6], [17],[18]. Also, there are exact computations of the sensitivity ofthe damping from a dynamic power grid model [7], [19], [20],[21] that are based on the eigenvalue sensitivity formula ∂λ∂p = wJ p vwv , (1)where w and v are left and right eigenvectors associated withthe eigenvalue λ and J p is the derivative of the Jacobian withrespect to the amount of generator redispatch p . The calcula-tion of J p involves the Hessian and the sensitivity of the oper-ating point to p . However, requiring a large scale power systemdynamic model poses difficulties. It is challenging to obtainvalidated models of generator dynamics over a wide area andparticularly difficult to determine dynamic load models thatwould be applicable online when poor modal damping arises.We think that a good way to solve the difficulties withPreprint c (cid:13) Personal use of this material is permitted. Permission from IEEE must be obtained for all other uses, in anycurrent or future media, including reprinting/republishing this material for advertising or promotional purposes, creating new collectiveworks, for resale or redistribution to servers or lists, or reuse of any copyrighted component of this work in other works. a r X i v : . [ c s . S Y ] S e p online large scale power system dynamic models is to combinemodels with synchrophasor measurements to get actionableinformation about mode damping. In particular, the dynamicinformation about the power system can be estimated fromsynchrophasors. However, formula (1) is not suitable for thispurpose since it is not feasible to estimate from measurementsthe left eigenvector w in (1) (or derivatives of eigenvectorsin other versions of (1)). This was a primary motivation fordeveloping our new formula in [11]. In particular, the formulashows that the first order effect of a generator redispatchlargely depends on the mode shape (the right eigenvector)and power flow quantities that can be measured online. Theassumed equivalent generator dynamics only appears as afactor common to all redispatches. Given a lightly dampedinterarea mode of a system, the method can rank all thepossible generators pairs. The rank is based on the size ofthe change of the damping ratio of the interarea mode.The main objective of ranking the generator pairs is toprovide advice to the operator of several effective generatorredispatches to damp the oscillations from which a correctiveaction can be selected. There are many economic goals andoperational constraints governing the final selection of anappropriate dispatch by the operator, and the integration ofthis decision with optimal power flow and the power marketsis left to future work.In this paper, the role in the method of measurements andthe formula are first illustrated in a simple 3-bus system. Thenthe general formula derived in [11] is presented. How staticand dynamic power quantities are obtained from measurementsis discussed. The complex denominator of the formula containsthe effect of the assumed equivalent generator dynamics, anda technique to obtain from measurements the phase of thiscomplex denominator is presented. The paper then explainshow to calculate the best generator pairs to damp a givenmode, and illustrates and verifies the calculation with interareamodes of the New England 10-generator system.II. F ORMULA FOR EIGENVALUE SENSITIVITY WITHRESPECT TO GENERATOR REDISPATCH
A. Special case of a 3-bus system
We first present the eigenvalue sensitivity formula for aspecial case in which the formula simplifies. Consider theinterarea mode of the simple 3-bus, 3-generator system shownin Fig. 1. The generator dynamics is described by the swingequation and the transmission lines are lossless. In this specialcase, the bus voltage magnitudes are assumed constant. At thebase case, the system is at an stable operating point and Fig. 1shows the line power flows and the interarea oscillating modepattern. The mode pattern shows that generator G1 is swingingagainst G3, and that G2 is not participating in the oscillation.The bus voltage phasor angles are δ , δ , δ . Let θ = δ − δ and θ = δ − δ be the voltage angle differences across line 1and line 2, and let p and p be the real power flows throughline 1 and line 2. The interarea mode has complex eigenvalue λ and complex right eigenvector x . Generator redispatch causeschanges in the angles across the lines dθ and dθ , and these in (cid:126) (cid:126) (cid:126) G2G1 G3 p p Fig. 1.
The gray lines joining the buses show the magnitude of thepower flow with the grayscale and the direction of the power flowwith the arrows. The red arrows at each bus show the oscillation modeshape; that is, the magnitude and direction of the complex entries ofthe right eigenvector x associated with the interarea mode. turn cause changes dλ in the eigenvalue. According to [11],the formula for the first-order change in λ with respect togenerator redispatch in a 3-bus system with constant voltagemagnitudes reduces to dλ = ( x (cid:48) θ ) p dθ + ( x (cid:48) θ ) p dθ α . (2)All quantities are in per unit unless otherwise indicated. Thecomplex denominator α = 0 . ∠ . ◦ seconds dependson the inertia and damping coefficients of the generators, theeigenvalue λ and its right eigenvector x . The right eigenvectorwritten with components corresponding to the bus angles is x = ( x δ , x δ , x δ ) t . The right eigenvector may also beexpressed in terms of the angles across the line as x (cid:48) =( x (cid:48) θ , x (cid:48) θ ) t . That is, x (cid:48) θ = x δ − x δ is the change in theeigenvector across line 1, and x (cid:48) θ = x δ − x δ is the changeof the eigenvector across line 2. It is these changes x (cid:48) θ and x (cid:48) θ in the right eigenvector across the lines that appear in (2).Formula (2) is in terms of static load flow quantities θ and p that are available from state estimation, and dynamic quantities λ and x that could be available from synchrophasor measure-ments. Table I shows the values at the base case. Formula (2)indicates which lines have suitable power flow and eigenvectorcomponents to affect oscillation damping. In particular, it iseffective for the redispatch to change the angles across thelines that have both changes in the mode shape across the lineand sufficient real power flow in the right direction. TABLE IQ
UANTITIES FROM MEASUREMENTS
Static quantities from Dynamic quantitiesstate estimator from synchrophasorsLine No. θ p [pu] x (cid:48) θ . ◦ . ◦ Let us define the complex coefficients of dθ and dθ in (2)as C θ and C θ respectively. Then dλ = C θ dθ + C θ dθ , where (3) C θ = − . − j . , C θ = 0 . − j . . Coefficient C θ has the largest real and imaginary components.It follows that dλ is more sensitive to redispatches donethrough line 1 than redispatches done through line 2. The Since x is the right eigenvector of a quadratic formulation of theeigenvalue problem [11], it contains the angle components, but not thefrequency components, of the conventional right eigenvector. That is,the conventional right eigenvector is ( x δ , x δ , x δ , x ω , x ω , x ω ) t =( x δ , x δ , x δ , λx δ , λx δ , λx δ ) t = ( x, λx ) t . damping ratio at the base case is given by ζ = − Re { λ }| λ | = − σ | σ + j ω | = − σ √ σ + ω . (4)Table II shows the results for the three possible generatorspairs in the system for a small redispatch of 0.01 pu. Asexpected from (3), the damping ratio has the largest changeswhen redispatch is done through line 1. This is the case ofgenerator pairs G1+,G2- and G1+,G3-; their changes in dζ are close. The damping ratio has almost no change whenredispatch is done only through line 2 with the pair G3+,G2-.Gi+,Gj- indicates the redispatch that increases power at gen-erator Gi and decreases power at generator Gj. TABLE IIG
ENERATOR PAIRS RANKED BY CHANGE IN ζ (%) ; REDISPATCH = 0.01Generator Pair dζ (%) B. General case of formula for the change in the eigenvalue
Consider a connected power grid that has m generators, n + m buses and (cid:96) lines. Buses n + 1 , . . . , n + m are theinternal buses of the generators. Assume AC power flowand lossless transmission lines. Every generator is modeledwith an equivalent second order equation (swing equation)and constant internal voltage magnitude. Loads are constantpower, but frequency dependence of the real power and voltagedependence of the reactive power can be accommodated[11]. Then the power system dynamics is described by aset of differential-algebraic equations with variables ( δ, V ) .The dynamic variables are δ n +1 , . . . , δ n + m , and the algebraicvariables are δ , . . . , δ n , V , . . . , V n . In this paper we use thequadratic formulation of the eigenvalue problem [11], [22],[23] with right eigenvector x and eigenvalue λ : ( M λ + Dλ + L ) x = 0 , (5)where L is part of the system Jacobian describedin [11] and M and D are diagonal matrices con-taining the generator equivalent inertias and dampings; M = diag { h /ω , h , /ω , . . . , h m /ω , , . . . , } and D = diag { d , d , . . . , d m , , . . . , } . The componentsof the right eigenvector or mode shape are x =( x δn +1 , . . . , x δn + m , x δ , . . . , x δn , x V , . . . , x Vn ) t .In other papers, the eigenstructure of differential-algebraicmodels with extended Jacobian J is often analyzed usingthe generalized eigenvalue problem Jv = λEv [24]. Thedifference between x and v is that x does not include thecomponents of v related to generator angular speeds.According to [11], the new formula for the sensitivity of anonresonant (algebraic multiplicity one) eigenvalue λ is dλ = 1 α (cid:32) (cid:96) (cid:88) i =1 ˜ C θ k dθ k + n (cid:88) i =1 ˜ C V i dV i (cid:33) . (6) The denominator of (6), which is the same for all redispatchesof a given mode, is α = 2 λx T M x + x T Dx, (7) dθ k is the change in angle across the line k and dV i is thechange in load i voltage magnitude due to the redispatch. ˜ C θ k = [( x (cid:48) θ k ) − ( x (cid:48) ν k ) ] p k + 2 x (cid:48) θ k x (cid:48) ν k q k , (8) ˜ C V i = (cid:96) (cid:88) k =1 | A ik | ( − C q k q k − C p k p k ) V − i − C Q i Q i V − i , (9) C q k = x (cid:48) ν k (cid:18) x (cid:48) ν k − x Vi V i (cid:19) − ( x (cid:48) θ k ) ,C p k = 2 x (cid:48) θ k (cid:18) x (cid:48) ν k − x Vi V i (cid:19) , C Q i = − (cid:18) x Vi V i (cid:19) ,p k = b k V i V j sin θ k , q k = − b k V i V j cos θ k , where θ k = (cid:26) δ i − δ j if bus i is sending end of line kδ j − δ i if bus i is receiving end of line kx (cid:48) θ k = (cid:40) x δi − x δj if bus i is sending end of line kx δj − x δi if bus i is receiving end of line kx (cid:48) ν k = x Vi V i + x Vj V j if line k joins load bus i to load bus jx Vi V i if line k joins load bus i to generator bus jp k is the real power flow in line k , q k is part of the reactivepower flow in line k , and | A ik | is the absolute value of the ik -component of the bus-line incidence matrix. Q i is the reactivepower demanded by load bus i .Formula (6) expresses the first-order change in the eigen-value dλ in terms of the changes dθ in angles across the linesand changes dV in load voltage magnitudes caused by thegenerator redispatch. dλ depends linearly on the active powerflow, part of the reactive power flow through every line ofthe network, and the reactive power demands at the loads, allevaluated at the base case. C. Generator modeling
The overall dynamics of each generator is described by anequivalent swing equation model with constant internal voltagemagnitude. For generator bus i , an “equivalent swing equa-tion” is the standard swing equation model with inertia anddamping coefficients h i and d i that produce the second ordermodel that best approximates the entire dynamics of the gener-ator i and its controls (this would be described as a second or-der dominant pole approximation in automatic controls). In ourapproach, there is no need to model the parameters h i and d i for each generator; it is only necessary to assume the existenceof a second order model that can approximate well enough thecontribution of the generator to the electromechanical mode.Indeed, since (6) only includes the combined generator dynam-ics as a common factor that is the same for all redispatches,we do not need to know the individual parameters of eachequivalent generator model. Other authors using synchropha-sor measurements to identify aggregated generator dynamics for studying oscillations also assume swing equation generatormodels but for their purposes identify the individual inertia anddamping coefficients [25], [26]. However, we have to considerthe generator modeling independently of previous studies sinceour generator redispatch application has novel and differentmodeling requirements as discussed in section VII.III. M EASUREMENTS
Formula (6) depends on power system quantities that couldbe observed from measurements from the state estimator andfrom synchrophasors.
A. Quantities obtained from state estimator: p, q, Q, dθ, dV
The state estimator can determine the active power p andpart of the reactive power flow q through every line of thenetwork, and the load reactive power injections Q . Load flowequations can be used to relate the generator redispatchesto changes in angles across the lines dθ and load voltagemagnitudes dV . B. Quantities obtained from synchrophasors: λ, x
Formula (6) depends on the dynamic quantities λ and x that satisfy the quadratic formulation of the eigenvalueproblem (5). For an electromechanical oscillation present ina system, synchrophasors can make online measurements [9],[10], [27], [28] of the damping and frequency of the eigenvalue λ associated with the oscillation. x is easy to obtain from aconventional right eigenvector. The right eigenvector of λ is inprinciple, and to some considerable extent in practice, avail-able from ambient or transient synchrophasor measurements[13], [14], [15]. It is conceivable that historical observationsor offline computations or general principles about the modeshape could be used to augment or interpolate the real-timeobservations, or that the real-time observations could be usedto verify a predicted mode shape. Thus some combination ofmeasurements and calculation could yield the mode shape x . C. Method for estimating phase of α from measurements The complex denominator α = 2 λx T M x + x T Dx isthe only part of formula (6) that depends on the generatorequivalent dynamic parameters in the inertia and dampingmatrices M and D . In estimating the change dλ in theeigenvalue from (6), the angle ∠ α controls the direction of dλ and the size of α controls the size of dλ . For each mode, α is the same for all generator redispatches. Therefore, inorder to rank the generator redispatches for a given mode,it is sufficient to estimate ∠ α from measurements.Formula (6) can be summarized as dλ = Numerator α , so ∠ α = ∠ Numerator − ∠ dλ. (10)To estimate ∠ α we propose to take advantage of the smallrandom load variations around the operating equilibrium.For such small random load variations, samples of dλ couldbe obtained from a series of synchrophasor estimates of λ .Samples of dθ and dV can be obtained from the load flow equations with simulated random load variations, so thatsamples of the numerator of (6) can be computed. It turns outthat both the dλ samples and numerator samples have preferredor principal directions in the complex plane. Analyzing the dλ and numerator samples with Principal Component Analysisgives a principal axis direction for dλ and a principal axisdirection for the numerator, and according to (10), the anglebetween these principal axis directions can be used to find ∠ α . In section VI-B, we illustrate and apply this calculationof ∠ α to interarea modes of the New England system. D. Measurement processing requirements
To be able to supply the dynamic quantities in formula(6), online synchrophasor monitoring of the oscillatory modeeigenvalue λ and mode shape x is needed. The modal eigen-value is used directly in (6) and, according to the suggestedapproach in section III-C, for estimating the phase of α . Thissubsection comments briefly about the likely measurementprocessing limitations. We would expect to be able to usestandard sampling rates and data concentration similar to thatalready in use for mode monitoring.Methods to estimate oscillatory modes and mode shapesfrom synchrophasor data are deployed and improving, andsome recent methods [14], [29], [30], [31], have used up to5 minutes of ambient or transient data for multiple parallelalgorithms to converge to consistent results. There is also con-sistency between results from methods for transient responseafter a disturbance and ambient methods. However, it does notmatter for our algorithm whether the dynamics is estimatedfrom transient or ambient responses.The mode shape estimation currently samples the modeshape at multiple spatial locations. We require the modeshape at other locations to be interpolated, perhaps guidedby previously observed mode shapes for specific modes. Ourapproach also requires static quantities to be estimated withstandard state estimation. Given state estimator convergence,the estimated state should be readily available within theseveral minute time scale already required for mode estima-tion. We do not anticipate significant computational delay inevaluating formula (6) and ranking the generators.IV. R ELATING REDISPATCH TO THE CHANGES dθ , dV , dλ Formula (6) expresses the eigenvalue change dλ in terms of dθ and dV . It remains to express dθ and dV in terms of theredispatch dP . Define the coefficient vectors of dθ and dV in(6) as C θ = 1 α ˜ C θ and C V = 1 α ˜ C V . (11)Then dλ = C θ · dθ + C V · dV = ( C θ , C V ) (cid:18) dθdV (cid:19) . (12)From [11] we know that the linear relationship betweenthe changes in angles dθ and changes in voltages dV withredispatch dP is given by (cid:18) dθdV (cid:19) = (cid:32) A T(cid:96) × ( n + m ) (cid:96) × n n × ( n + m ) I n × n (cid:33) L † (cid:18) dP (cid:19) , (13) where A is the network incidence matrix, and † indicatespseudoinverse. As an alternative to the linearized computationof ( dθ, dV ) from the redispatch dP in (13), one can simplyrecompute the AC load flow with the assumed change inredispatch dP . We can relate the change in an eigenvalue ofa mode to a generator redispatch dP by substituting (13) into(12): dλ = C P · dP = m (cid:88) i =1 C P i dP i . (14)The complex coefficient C P i gives the contribution of gener-ator i to dλ by increasing the real power of generator i by dP i . The redispatch is assumed to satisfy the active powerbalance constraint (cid:80) mi dP i = 0 . It is convenient to define ˜ C P = | α | C P . Then multiplying both sides of (14) by | α | gives | α | dλ = ˜ C P · dP = m (cid:88) i =1 ˜ C P i dP i . (15)It can be seen by considering (11)–(14) that ˜ C P can becalculated from ˜ C θ , ˜ C V and ∠ α . (Note, for example, that | α | C θ = e − ∠ α ˜ C θ .) Since | α | dλ is dλ multiplied by a realconstant, the complex number | α | dλ has the same angle as dλ and has magnitude proportional to dλ . It follows that wecan use | α | dλ to rank the generator redispatches.V. R ANKING GENERATORS PAIRS BY INCREASE INDAMPING RATIO DUE TO REDISPATCH
For a given oscillatory mode with insufficient damping ratio,the first step towards ranking the best generators to redispatchis to make the following measurements and calculations:1) The eigenvalue λ associated with the interarea oscil-lation and its mode shape x are estimated from syn-chrophasor measurements.2) p, q, Q are obtained from the state estimator.3) ∠ α is computed from measurements.4) Coefficients ˜ C θ k and ˜ C V i in (6) are computed from (8)and (9).5) ˜ C P is computed from ˜ C θ , ˜ C V , ∠ α .The most straightforward way of implementing generatorredispatch is with a pair Gi+,Gj- of generators; that is, in-creasing the real power of generator i by some amount dP i and increasing the real power of generator j by the negativeamount dP j = − dP i . We can use (14) to find the generatorpair Gi+,Gj- that gives the most favorable eigenvalue change dλ that best increases the damping ratio. For a redispatchof amount dP i with generator pair Gi+,Gj-, the resultingeigenvalue change satisfies | α | dλ ij = ( ˜ C P i − ˜ C P j ) dP i . (16)Also the opposite redispatch Gj-,Gi+ with the same generatorpair will yield | α | dλ ji = −| α | dλ ij . For every pair of gener-ators in the system we can compute | α | dλ ij and | α | dλ ji andthen we can compute the corresponding damping ratios ζ ij and ζ ji for every generator pair with (4). The damping ratiosfor all pairs { ζ , ζ , . . . , ζ m − ,m , ζ m,m − } are then ranked from the largest one to the smallest one. The highly rankedpairs indicate the generator pairs and redispatch directions thatwill increase the damping ratio of the mode the most.VI. D AMPING MODES OF THE N EW E NGLAND SYSTEM
This section presents the results for damping interareamodes of the New England 10-generator system [32] withgenerator redispatch. The equivalent parameters of the gen-erators are given in Appendix A and the rest of the data isprovided in [32]. All the numerical computation was done withthe software Mathematica. Table III shows the eigenvaluesof the 4 interarea modes at the base case. Since the systemhas 10 generators, there are 45 generators pairs. The methodpresented in section V was implemented for the interareamodes (this computation used a calculated value of α anda recalculated load flow to evaluate dθ and dV from dP ). TABLE IIII
NTERAREA MODE EIGENVALUES AT THE BASE CASE
Mode Eigenvalue λ [1/s] f [Hz] Damping Ratio ζ (%)1 -0.040336 + j3.4135 0.54327 1.181572 -0.018839 + j4.7631 0.75807 0.395513 -0.024903 + j5.4994 0.87526 0.452834 -0.055799 + j6.0159 0.95746 0.92748 A. Verifying the formula
Formula (6) was verified for the four interarea modes, buthere we present results only for mode 1. The eigenvalue λ was computed for a given redispatch using (6) and using theexact computation of the Jacobian at the new operating point.Table IV shows the exact and approximate λ for differentamounts of redispatch of generator pair G5+,G9-. Table IV confirms that (6) reproduces the first order varia-tion of λ with respect to the redispatch. TABLE IVE
IGENVALUE λ FOR GENERATOR REDISPATCH
G5+, G9-Redispatch Exact eigenvalue Approximate eigenvalue0.0 -0.0403355 + j3.4135 -0.0403355 + j3.41350.0005 -0.0403225 + j3.4131 -0.0403225 + j3.41310.001 -0.0403095 + j3.4127 -0.0403095 + j3.41270.01 -0.0400715 + j3.4059 -0.0400736 + j3.40590.02 -0.0397985 + j3.3981 -0.0398086 + j3.39830.03 -0.0395161 + j3.3899 -0.0395403 + j3.3905
B. Estimating the phase of α for random load variations Table V shows the exact ∠ α computed with (7) from themode shape and the equivalent generator parameters, and the ∠ α estimated from random load variations with the methodof section III-C. The base case generations of the generators G1 through G10 are . , . , . , . , . , . , . , . , . , . per unit respectively. As anexample of specifying the redispatch, a 0.03 redispatch of generator pairG5+,G9- changes the generation of G5 from 5.08 to 5.11 and the generationof G9 from 8.3 to 8.27. For each interarea mode, the set of random loads used in thecomputation were generated with the software Mathematica.The active power random load vector P r was sampled froma normal distribution of zero mean The components of thereactive power random load vector Q r were computed as Q ri = P i Q i P ri . (17) P i and Q i are the active and reactive power demanded byload i at the base case. Then the load flow solution for thevector of loads P + P r and Q + Q r were computed and dλ and the numerator of (6) were computed. 50 randomload scenarios were generated. Fig. 2 shows samples for λ after being trimmed by 30% to remove outliers and analyzedwith principal component analysis. This 2-dimensional λ data was trimmed with multidimensional trimming based onprojection depth [33]. Projection depth induces order for highdimensional data, which makes trimming straightforward. Theresults in Table V show that the method gives a very goodestimation of ∠ α . TABLE VE
STIMATED AND EXACT PHASES OF α FOR THE INTERAREA MODES
Mode Exact ∠ α Estimated ∠ α Exact ∠ α - Estimated ∠ α . ◦ ◦ − . ◦ . ◦ ◦ − . ◦ . ◦ ◦ . ◦ . ◦ ◦ − . ◦ C. Ranking generator pairs to damp mode 1
Fig. 3 shows the power flow p and oscillating mode patternof interarea mode 1 at the base case. The mode pattern showsthat generators G2 through G9 are oscillating against G10.The component of G10 is not very large compared withthe components of the other generators, but G10 is a largegenerator that represents an equivalent of the New York Stategrid. From Fig. 3 we can see that generator G5 participatesmost in the oscillation, G8 has a small participation in theoscillation, and G1 does not participate in the oscillation. Themethod of section V was applied to rank the generator pairsthat best increase the damping ratio for a small redispatch of0.01 pu, and the top 10 generator pairs are shown in Table VI.The top 9 generator pairs all include G5-; that is, decreasinggeneration at G5 with increasing generation elsewhere. Thechange in damping ratio for these pairs are of the same order ofmagnitude, so it is clear from Table VI that the largest changesin damping ratio are due to the generator pairs involving G5-. TABLE VIG
ENERATOR PAIRS RANKED BY CHANGE IN ζ (%) ; REDISPATCH = 0.01Generator Pair dζ (%) Generator Pair dζ (%) Load 12 has exceptionally high reactive power and was not varied. (cid:230) (cid:230) (cid:230)(cid:230) (cid:230)(cid:230)(cid:230) (cid:230)(cid:230) (cid:230) (cid:230)(cid:230) (cid:230)(cid:230)(cid:230) (cid:230)(cid:230) (cid:230) (cid:230)(cid:230) (cid:230) (cid:230)(cid:230) (cid:230)(cid:230) (cid:230) (cid:230)(cid:230) (cid:230)(cid:230) (cid:230)(cid:230)(cid:230)(cid:230) (cid:230) (cid:230)
PrincipalAxis (cid:45) (cid:45) (cid:45) (cid:45) (cid:230)(cid:230) (cid:230) (cid:230)(cid:230) (cid:230)(cid:230) (cid:230)(cid:230) (cid:230)(cid:230) (cid:230) (cid:230)(cid:230) (cid:230)(cid:230)(cid:230) (cid:230)(cid:230)(cid:230)(cid:230) (cid:230)(cid:230) (cid:230) (cid:230)(cid:230) (cid:230)(cid:230) (cid:230)(cid:230) (cid:230) (cid:230) (cid:230)(cid:230)(cid:230) (cid:230)
PrincipalAxis (cid:45) (cid:45) (cid:45) (cid:56) d Λ (cid:60) (cid:45) (cid:45) (cid:45) (cid:56) d Λ (cid:60) Fig. 2.
50 samples of the numerator of (6) and dλ after trimmingby 30%. Principal axes are computed and shown as lines. G5G7 G4 G3 G2G10G1G8G9G6
Fig. 3.
The gray lines joining the buses show the magnitude ofthe active power flow with the grayscale and the direction of theactive power flow with the arrows. The red arrows at each bus showthe pattern of the oscillation for mode 1; that is, the magnitude anddirection of the entries of the right eigenvector x δ . D. Ranking generator pairs to damp mode 2
Fig. 4 shows the power flow p and oscillating mode patternof interarea mode λ at the base case. The mode pattern showsthat generators G1-G3 and G6-G9 are oscillating against G5.Generators G5 and G9 participate most in the oscillation.The method of section V was applied to rank the generatorpairs that best increase the damping ratio for a small redispatch G5G7 G4G3 G2 G10G1G8G9G6
Fig. 4.
The gray lines joining the buses show the magnitude of theactive power flow with the grayscale and the direction of the activepower flow with the arrows. The red arrows at each bus show thepattern of the oscillation for mode 2. of 0.01 pu, and the top 10 generator pairs are shown inTable VII. The top 9 generator pairs all include G5+; thatis, increasing generation at G5 with decreasing generationelsewhere. Although the change in damping ratio is of thesame order of magnitude for the 10 pairs, the increase indamping ratio for pair 10 G4+,G9-, is roughly half of theincrease of damping ratio for pair 9 G5+,G4-.Now we analyze the most significant components of formula(6) to explain why G5 is playing a key role in damping mode 2.Table VII shows that the pairs with the largest change in
TABLE VIIG
ENERATOR PAIRS RANKED BY CHANGE IN ζ (%) ; REDISPATCH = 0.01Re { C V } Gen. Pair dζ (%) dλ Re { C θ } · dθ · dV damping ratio dζ are the ones with the largest increase in thedamping Re { dλ } , so we can focus on the real part of (12).Moreover, Table VII also shows that the changes Re { C θ } · dθ are larger than the changes Re { C V } · dV , so we focus onanalyzing the terms of (12) related to dθ :Re { dλ } = Re { C θ } · dθ + Re { C V } · dV ≈ Re { C θ } · dθ (18)Fig. 5 shows different quantities related to the 56 lines ofthe New England system in gray scale. Each generator isrepresented in the network by an internal bus and a terminalbus. There are two lines in series at the edge of the networkassociated with each generator. The line joining the internalbus to the terminal bus represents the generator transient G5 G9 (a) | Re { C θ }| for λ . G5 G9 (b) | dp | for redispatch G5+,G9-. G5 (c) | dθ | for G5+,G9-. G5 (d) | Re { C θ }· dθ | for λ & G5+,G9-.Fig. 5. Gray scale in the lines shows the components of the specified vector. reactance, and the other line lumps together the transformerand lines joining the generator terminal bus to the network.(G10 differs since it represents New York state.) Fig. 5(a)shows | Re { C θ }| , the absolute value of dθ ’s coefficient in (18).The lines that represents the transient reactance of G5 and thetransient reactance of G9 have large components of | Re { C θ }| .Fig. 5(b) shows | dp | , the absolute value of the change in powerflow in lines for the best ranked generator pair G5+,G9-. Asexpected, several lines have a significant change in powerflow, but Fig. 5(c) shows that only the line that representsthe transient reactance of G5 has a large change in angle | dθ | . This is due to the fact that transient reactance of G5is much larger than the reactance of any of the other 55lines of the system (see appendix A). So, for any possiblegenerator pair that involves G5, the line that represents thetransient reactance of G5 will always have the largest changein angle. This large change in angle, combined with a largecoefficient | Re { C θ }| , produces the dominant term of (18), asshown by Fig. 5(d). Thus G5 is the key generator to participatein redispatch for producing the largest changes in damping formode 2. Although Fig. 5(a) shows that there are other linesthat have large dθ ’s real coefficient, and Fig. 5(b) shows thatother lines have an important change in power | dp | , Fig. 5(c)shows that such lines do not have a large change in angle dθ ,and as a result their associated terms in (18) are not large, asseen in Fig. 5(d). This analysis shows a mechanism of howdamping by redispatch works by changing the angles acrosslines that have large coefficients C θ . E. Larger redispatches
Larger redispatches introduce some nonlinearity in thechange in the damping ratio. This nonlinearity arises fromthree sources: the change in the load flow, the change inthe eigenvalue dλ , and the change in the damping ratio.For a given size of redispatch of a generator pair, we cancompute the exact nonlinear change in the load flow dueto the redispatch, use (6) to linearly estimate the change inthe eigenvalue based on measurements, and then compute thenonlinear change in the damping ratio. This is a large signalapplication of (6), and the top 10 generator pairs for modes 1and 2 are shown in Table VIII for a redispatch of 0.4 pu. Forcomparison, Table VIII also shows the exact calculation thatuses the nonlinear computation of dλ . As shown in Fig. 6, theuse of (6) for a larger redispatch gives the same grouping ofthe top 9 effective generator pairs and a similar ranking andgrouping of the top 10 generator pairs as the exact calculation.An exception is that for mode 1, the 5th ranked generator pairG9+,G5- for a larger redispatch using (6) becomes the 9thranked generator pair for the exact calculation.Fig. 6 also compares the changes in the damping ratio ofthe top 10 generator pairs for a small redispatch using (6)to the larger redispatch using (6). The grouping of the top9 effective generator pairs is preserved and the approximateranking is preserved.Overall, some details of the rankings differ for similarlyeffective generator pairs, but since the ranking will be usedto provide a set of effective generators pairs from which anoperationally suitable pair can be selected for redispatch byoperators, the performance of the ranking is satisfactory. TABLE VIIIG
ENERATOR PAIRS RANKED BY CHANGE IN ζ (%) ; REDISPATCH = 0.4MODE 1 dζ (%) MODE 2 dζ (%) Gen. Pair Using (6) Exact Gen. Pair Using (6) Exact1 G6+,G5- 0.181 0.138 1 G5+,G9- 0.563 0.6002 G7+,G5- 0.177 0.134 2 G5+,G10- 0.558 0.6153 G3+,G5- 0.177 0.135 3 G5+,G1- 0.534 0.5954 G2+,G5- 0.171 0.128 4 G5+,G8- 0.506 0.5725 G9+,G5- 0.162 0.092 5 G5+,G2- 0.473 0.5456 G4+,G5- 0.160 0.122 6 G5+,G3- 0.453 0.5277 G8+,G5- 0.159 0.115 7 G5+,G7- 0.434 0.5068 G1+,G5- 0.155 0.112 8 G5+,G6- 0.432 0.5039 G10+,G5- 0.145 0.104 9 G5+,G4- 0.365 0.44410 G6+,G10- 0.030 0.034 10 G4+,G9- 0.139 0.120
VII. D
ISCUSSION OF GENERATOR MODELING
Since our approach depends in new ways on both mea-surements and an equivalent second-order swing equationgenerator model, our generator modeling requirements aredifferent than in other approaches to suppressing oscillations.Section II-C explains that the generator dynamics are approx-imated by a swing equation, but we do not need to determinethe parameters of the swing equation for any individual gen-erator. This section further discusses the generator modeling.As a general observation, in closed loop control of oscil-lations with power system stabilizers, which forms most ofthe literature on suppressing oscillations, the generator and itscontrols need to be modeled in sufficient detail. Indeed the d Ζ (cid:72) (cid:37) (cid:76) I n c r ea s e MODE 1 MODE 2Using formula (6) Using formula (6)Redispatch 0.01 pu Redispatch 0.01 puExact ExactRedispatch 0.4 pu Redispatch 0.4 pu
Fig. 6. Changes in damping ratio for Mode 1 and Mode 2 for the top10 generator pairs from Tables VI, VII, VIII rescaled to the same range toallow comparison. For each mode, left hand dots use formula (6) for a smallredispatch, middle dots use (6) for a larger redispatch, and the right handdots are the exact calculation for the larger redispatch. Clustering of similarlyeffective generator pairs is shown by close dots and changes in ranking appearas lines crossing. designed control gains directly affect entries of the Jacobianto damp the oscillatory mode. Our control is open loop andworks by the entirely different principle of exploiting systemnonlinearity by changing the operating point at which theJacobian is evaluated. This changes the focus from the linearparts of the model to the nonlinearities.Formula (6) computes the first order sensitivity of the oscil-latory mode eigenvalue to generator redispatch. Appendix Bproves that the first order eigenvalue sensitivity to redispatchdoes not depend on linear parts of the power system model. Inparticular, if the generator magnetic saturation and hysteresisare neglected, the higher order parts of the generator modelingare linear and the eigenvalue sensitivity only depends on thenonlinearity in the swing equation and any stator algebraicequations and does not depend on the linear higher order partof the generator dynamic modeling. This does suggest thatthe linear higher order generator dynamics can be omitted inderiving the formula.In applying formula (6), we do not use a model of thepower system dynamics. Instead we rely on measurementsof the system dynamics, particularly the eigenvalue and theright eigenvector of the mode and the phase of the complexscalar parameter α that combines together all of the generatordynamics. There are no model assumptions in these measuredquantities. That is, if part of the power system affects the oscil-latory dynamics, the effect will appear via the measurementsused by the formula.We did an initial test of the approximation involved ingenerator modeling with the interarea mode of the 3-generatormodel similar to Fig. 1 with a sixth-order generator model ateach bus. Formula (6) is applied to this detailed model bymeasuring the phase of α , and using the computed mode,mode shape, and load flow to estimate the change dλ inthe eigenvalue that would arise from small redispatches ingeneration. Then the detailed model is used to compute theexact change dλ in the eigenvalue. The comparison of theapproximate and exact dλ is shown in Table IX. Similarlyto the intended application of the formula to ranking redis-patches in a real power system, the power system modelwith sixth-order generators does not have the parameters ofthe equivalent second order generator models available. Thatis, the magnitude of α is not known, and only the phase of α is estimated, and so the formula predicts dλ to within a constant real multiplier. Therefore in Table IX we comparethe ratios of | dλ | for each redispatch to | dλ | for the redispatchof generators G1 and G2, as well as comparing the phases of dλ . The approximation of dλ in Table IX is close enough tobe acceptable for ranking of generator redispatches.While the generator modeling issues should be investigatedfurther in future work, and further analytic progress is notruled out, both theoretical considerations of the irrelevance oflinear parts of the generator model and an initial test indicatethat combining measurements of the dynamic quantities witha formula assuming a second order swing equation can beadequate for ranking generator redispatches. TABLE IXE
IGENVALUE CHANGES FOR REDISPATCH OF
PU IN GENERATORSYSTEM WITH SIXTH - ORDER GENERATOR MODELING
Generator Exact Approximate with formulapair ∠ dλ | dλ | ratio ∠ dλ | dλ | ratioG1+,G2- − . ◦ − . ◦ − . ◦ − . ◦ − . ◦ − . ◦ VIII. C
ONCLUSIONS
There has been success in monitoring interarea modaldamping with synchrophasor measurements [9], [10]; the nextstep is to leverage synchrophasor measurements to provideadvice to the operators to maintain a suitable modal dampingratio when the damping ratio is insufficient. Difficulties inaccomplishing this in the past include the lack of wide-areaonline dynamic models and standard formulas that dependon quantities that cannot be measured. In this paper, wecircumvent these difficulties by calculating the best generatorpairs to redispatch to maintain modal damping by combiningsynchrophasor and state estimator measurements with a newanalytic formula for the sensitivity of the mode eigenvalue withrespect to generator redispatch. The assumed equivalent gener-ator dynamics only appears as a complex factor /α commonto all redispatches and we propose a method of estimating thephase of α from ambient measurements. The new formula issomewhat complicated, and we explain and illustrate how itworks in 3 and 10 generator examples. Future work may welldiscover further insights and applications using the formula.In summary, we make substantial progress towards practicalapplication of a new formula to damp interarea oscillationsbased on measurable quantities.R EFERENCES[1] Cigr´e Task Force 07 of Advisory Group 01 of Study Committee 38,Analysis and control of power system oscillations, Paris, Dec. 1996.[2] G. Rogers, Power System Oscillations, Kluwer Academic, 2000.[3] IEEE Power system engineering committee, Eigenanalysis and fre-quency domain methods for system dynamic performance, IEEE Publi-cation 90TH0292-3-PWR, 1989.[4] IEEE PES Systems Oscillations Working Group, Inter-area oscillationsin power systems, IEEE Publication 95 TP 101, Oct. 1994.[5] C.Y. Chung, L. Wang, F. Howell, P. Kundur, Generation reschedulingmethods to improve power transfer capability constrained by small-signal stability, IEEE Trans. Power Syst., vol. 19, no. 1, pp. 524-530,Feb. 2004. [6] Z. Huang, N. Zhou, F.K. Tuffner, Y. Chen, D.J. Trudnowski, MANGO-Modal analysis for grid operation: A method for damping improvementthrough operating point adjustment, U.S. Dept. of Energy, Oct. 2010.[7] I. Dobson, F.L. Alvarado, C.L. DeMarco, P. Sauer, S. Greene, H. Eng-dahl, J. Zhang, Avoiding and suppressing oscillations, PSerc publication00-01, Dec. 1999.[8] A. Fischer, I. Erlich, Assessment of power system small signal stabilitybased on mode shape information, IREP Bulk Power System Dynamicsand Control V, Onomichi, Japan, Aug. 2001.[9] J.W. Pierre, D.J. Trudnowski, M.K. Donnelly, Initial results in electrome-chanical mode identification from ambient data, IEEE Trans. PowerSyst., vol. 12, no. 3, pp. 1245-1251, Aug. 1997.[10] R.W. Wies, J.W. Pierre, D.J. Trudnowski, Use of ARMA block process-ing for estimating stationary low-frequency electromechanical modes ofpower systems, IEEE Trans. Power Syst., vol. 18, no. 1, pp. 167-173,Feb. 2003.[11] S. Mendoza-Armenta, I. Dobson, A formula for damping interareaoscillations with generator redispatch, IREP Symposium - Bulk PowerSystem Dynamics and Control - IX Rethymnon, Greece, Aug. 2013.Also available online arXiv:1306.3590v2.[12] A. Fischer, I. Erlich, Impact of long-distance power transits on thedynamic security of large interconnected power systems, IEEE PortoPower Tech Conference, Porto, Portugal, Sept. 2001.[13] D.J. Trudnowski, Estimating electromechanical mode shape from syn-chrophasor measurements, IEEE Trans. Power Syst., vol. 23, no. 3, pp.1188-1195, Aug. 2008.[14] N.R. Chaudhuri, B. Chaudhuri, Damping and relative mode-shape esti-mation in near real-time through phasor approach, IEEE Trans. PowerSyst., vol. 26, no. 1, pp. 364-373, Feb. 2011.[15] L. Dosiek, N. Zhou, J.W. Pierre, Z. Huang, D.J. Trudnowski, Modeshape estimation algorithms under ambient conditions: A comparativereview, IEEE Trans. Power Syst., vol. 28, no. 2, pp. 779-787, May 2013.[16] E. Barocio, B. C. Pal, N. F. Thornhill, A. R. Messina, A dynamic modedecomposition framework for global power system oscillation analysis,IEEE Trans. Power Syst., vol. PP, issue: 99, 2014.[17] Z. Huang, N. Zhou, F. Tuffner, Y. Chen, D. Trudnowski, W. Mittelstadt,J. Hauer, J. Dagle, Improving small signal stability through operatingpoint adjustment, IEEE PES General Meeting, Minneapolis, MN, July2010.[18] R. Diao, Z. Huang, N. Zhou, Y. Chen, F. Tuffner, J. Fuller, S. Jin,J.E Dagle, Deriving optimal operational rules for mitigating inter-areaoscillations, Power Syst. Conf. & Exposition, Phoenix AZ, March 2011.[19] I. Dobson, F.L. Alvarado, C.L. DeMarco, Sensitivity of Hopf bifurca-tions to power system parameters, 31st Conference on Decision andControl, Tucson, Arizona, Dec. 1992.[20] H.K. Nam, Y.K. Kim, K.S. Shim, K.Y. Lee, A new eigen-sensitivitytheory of augmented matrix and its applications to power systemstability, IEEE Trans. Power Syst., vol. 15, pp. 363-369, Feb. 2000.[21] S. Wang, Q. Jiang, Y. Cao, WAMS-based monitoring and control ofHopf bifurcations in multi-machine power systems, Journal of ZhejiangUniversity Science A, vol. 9, No. 6, pp. 840-848, 2008.[22] B.E. Eliasson, D.J. Hill, Damping structure and sensitivity in the Nordelpower system, IEEE Trans. Power Syst., vol. 7, No. 1, Feb. 1992.[23] E. Mallada, A. Tang, Improving damping of power networks: powerscheduling and impedance adaptation, Conf. Decision and Control andEuropean Control Conf. (CDC-ECC), Orlando, FL, Dec. 2011.[24] T. Smed, Feasible eigenvalue sensitivity for large power systems, IEEETransactions on Power Systems, Vol. 8, No. 2, pp. 555-563, May 1993.[25] A. Chakrabortty, J. H. Chow, A. Salazar, A measurement-based frame-work for dynamic equivalencing of large power systems using WAMS,Innovative Smart Grid Technologies, Gaithersburg, MD, Jan. 2010.[26] N. Zhou, S. Lu, R. Singh, M. Elizondo, Calibration of reduced dynamicmodels of power systems using phasor measurement unit (PMU) data,North American Power Symp. (NAPS), Boston MA, Aug. 2011.[27] L. Vanfretti, J.H. Chow, Analysis of power system oscillations fordeveloping synchrophasor data applications, IREP Symposium BulkPower System Dynamics and Control VIII, Buzios, Brazil, Aug. 2010.[28] IEEE Task Force on Identification of Electromechanical Modes, Iden-tification of electromechanical modes in power systems, IEEE SpecialPublication TP462, June 2012.[29] G. Liu, J. Ning, Z. Tashman, V. Venkatasubramanian, P. Trachian,Oscillation monitoring system using synchrophasors, IEEE Power andEnergy Society General Meeting, San Diego CA, July 2012.[30] J. Ning, X. Pan, V. Venkatasubramanian, Oscillation modal analysisfrom ambient synchrophasor data using distributed frequency domainoptimization, IEEE Trans. Power Systems, vol. 28, no. 2, 2013, pp.1960-1968. [31] H. Khalilinia, L. Zhang, V. Venkatasubramanian, Fast frequency-domaindecomposition for ambient oscillation monitoring, IEEE Trans. PowerDelivery, vol. 30, no. 3, pp: 1631-1633, June 2015.[32] M. A. Pai, Energy function analysis for power system stability, Kluwer,Norwell, MA, 1989.[33] Y. Zuo, Multidimensional trimming based on projection depth, Annalsof Statistics, Vol. 34, No. 5, pp. 2211-2251, Oct. 2006. A PPENDIX
A: N EW E NGLAND GENERATOR DATAGen. Terminal Internal V InternalNo. Bus No. Bus No. h [s] d [s] x (cid:48) d Bus1 30 40 42.0 0.0267 0.031 1.05012 31 41 30.3 0.0161 0.0697 1.03883 32 42 35.8 0.0209 0.0531 1.04394 33 43 28.6 0.0243 0.0436 1.03485 34 44 26.0 0.0014 0.132 1.20986 35 45 34.8 0.0277 0.05 1.09417 36 46 26.4 0.0140 0.049 1.09448 37 47 24.3 0.0116 0.057 1.07059 38 48 34.5 0.0002 0.057 1.125210 39 49 500.0 0.3979 0.006 1.0317 V i is the internal constant voltage magnitude of generator i at the base case;i.e., at zero redispatch, and x (cid:48) d is the transient generator reactance. A PPENDIX
B: I
RRELEVANCE OF LINEAR MODELING
This appendix proves from (1) that the first order eigenvaluesensitivity to redispatch does not depend on linear parts of thepower system model. This result was first mentioned in [7,section 4.6]. It is convenient to use the extended differential-algebraic form of the power system equations [24] with statevector z , Jacobian ¯ J , and extended right and left eigenvectors ¯ v and ¯ w . In this notation, (1) becomes ∂λ∂p = ¯ w ¯ J p ¯ v ¯ w ¯ v (19)In this case, the parameter p parameterizes the generatorredispatch, and it appears linearly in the system equations.Therefore p does not appear explicitly in the Jacobian ¯ J , and ¯ J p = ∂ ¯ J∂p = (cid:88) k ∂ ¯ J∂z k z kp (20)where z p = ∂z∂p is the sensitivity of the operating point z tothe redispatch. Then substituting in (19) and writing it out incoordinates gives ∂λ∂p = (cid:88) i,j,k ¯ w i ∂ ¯ J ij ∂z k z kp ¯ v j (cid:88) i ¯ w i ¯ v i (21)It is clear from (21) that the linear parts of the model vanishin ∂ ¯ J ij ∂z k and that if a dynamic state is associated with alinear differential equation, the corresponding entry of ¯ w getsmultiplied by zero in the numerator of (21). Sarai Mendoza-Armenta (M 13) received the PhD in Physics from Institutode F´ısica y Matem´aticas, Universidad Michoacana, Mexico in 2013. She wasvisiting scholar at Iowa State University from March 2012 to March 2013.She was post-doctoral research associate in the Electrical and ComputerEngineering Department at Iowa State University.