Quantifying Bounds of Model Gap for Synchronous Generators
aa r X i v : . [ m a t h . O C ] F e b Quantifying Bounds of Model Gap for Synchronous Generators
Peng Wang,
Member, IEEE,
Shaobu Wang,
Senior Member, IEEE,
Renke Huang,
Member, IEEE,
Zhenyu Huang,
Fellow, IEEE
Abstract —In practice, uncertainties in parameters and modelstructures always cause a gap between a model and the corre-sponding physical entity. Hence, to evaluate the performance ofa model, the bounds of this gap must be assessed. In this paper,we propose a trajectory-sensitivity–based approach to quantifythe bounds of the gap. The trajectory sensitivity is expressedas a linear time-varying system. We thus first derive severalbounds for a general linear time-varying system in differentscenarios. The derived bounds are then applied to obtain boundsof the model gap for generator plant models with different typesof structural information, e.g., models of different orders. Casestudies are carried out to show the efficacy of the bounds throughsynchronous generator models on different accuracy levels.
Index Terms —linear time-varying systems; model validation;model structure difference; trajectory sensitivity.
I. I
NTRODUCTION
Integrity of plant models in power systems is crucial fora reliable and robust power grid, which is key to economi-cal electricity delivery to power consumers. Also, long-termplanning and short-term operational decisions depend heavilyon studies and investigations using plant models. Imprecisemodels will result in conclusions that deviate from reality,which jeopardizes system designs and affects decision-making.Model mismatch has been an important cause of a fewblackouts in North America, e.g., the 1996 Western Statesoutage and the 2011 blackout in San Diego. In the studyon these blackouts, simulations with imprecise models showthat systems are stable and well-damped, but these systemseventually collapse, which contradicts the simulation studywith mismatched models. Thus, investigating the mismatchbetween models and plants is important for a reliable, robust,and economical grid. Much effort has been devoted to suchinvestigation, to improve model quality.Model calibration has been a hot research topic in thelast decade. The basic idea is to inject measurement data,e.g., voltage magnitude and frequency/phase angle, into thepower plant terminal bus in dynamic simulation and comparethe simulated response with the actual measured data, as in[1], [2]. Such a method is called ”event playback” and hasbeen included in major software vendor tools, e.g., GE PSLF,Siemens PTI PSS/E, Powertech Labs TSAT, and PowerWorldSimulator [3]. Once model discrepancy is observed with event
Financial support for this work was provided by the Advanced ScientificComputing Research (ASCR) program of the U.S. Department of Energy(DOE) Office of Science and the Advanced Grid Modeling (AGM) programof the DOE Office of Electricity. Pacific Northwest National Laboratory isoperated by Battelle for the DOE under Contract DE-AC05-76RL01830.Peng Wang, Shaobu Wang, and Zhenyu Huang are with Pa-cific Northwest National Laboratory, 902 Battelle Boulevard, Rich-land, WA. 99354. Corresponding emails: [email protected] , [email protected] , [email protected] . playback or another method, model identification and calibra-tion are needed. Methods and tools to calibrate models aredeveloped using linear or nonlinear curve fitting techniques,either for power plants as in [4], [5] or for general plants,as in [6], in which the methods can be applied to powerplants. But it has been reported that multiple solutions mayexist for the same model performance with the methods andtools in [3], [5], [6], making it difficult to identify the trueparameters of the model. Simultaneous perturbation stochasticapproximation methods based on particle swarm optimizationapproaches are developed in [7], [8]. Such methods are reliableand sufficiently precise, but they are time-consuming, forthey require extensive iterations of simulations to obtain goodperformance. A probabilistic collocation method is used toevaluate the uncertainty in power systems with a few uncertainparameters in [9]. Feature-based searching algorithms areinvestigated in [10]–[12], which help identify the problematicparameters but lack calibration accuracy. A dynamic-state-estimation-based generator parameter identification algorithmis proposed in [13] with good results for the test case, but ithas been shown to be very sensitive to measurement noise.A rule-based approach is reported in [14], which classifiesparameters based on their effects on simulation responses. Butthis approach requires extensive trial and error and knowledgeof the physics of the system, which may not be acquired withenough precision. A method and a tool suite are developed tovalidate and calibrate models based on trajectory sensitivityanalysis and advanced ensemble Kalman filters in [15].The latest results show that most existing approaches canhandle parameter errors well but model structure errors arevery difficult to be calibrated. For example, if one uses asecond-order model to approximate a fourth-order model’sresponse. No matter how to tune the parameters of thesecond-order model, their response curves won’t match. Thismotivates us to study the bound of model structures (modelgap). Such quantification will enable us to know the rangeof the gap in advance and be prepared to take precautions toensure power systems run reliably, robustly, and economically.This is important for both long-term planning and short-termoperational decisions.Quantifications of the model gap are investigated based ontrajectory sensitivity in this paper. Two bounds for modelgaps are derived via analysis of trajectory sensitivity withrespect to imprecise structural information. Because impreciseinformation probably results from changes of parameters andunmodeled dynamics, the first bound quantifies model gapscaused by imprecise parameters and the second one by unmod-eled dynamics. Scenarios for model gaps caused by changesof parameters and unmodeled dynamics are both simulated toshow the efficacy of the derived bounds. Unlike existing resultsn using trajectory sensitivity for model gaps in power systemsin the literature, e.g., [16], explicit bounds of the trajectorysensitivity are obtained in this paper.The rest of this paper is organized in a generic-to-specificway. That is, we first present the results for trajectory sensi-tivity for general models and then apply those general resultsto quantify the model gaps of synchronous generators. Somepreliminary information is introduced in Section II. Severalbounds for linear time-varying (LTV) systems are derived inSection III. We apply the bounds to quantify the model gapsof synchronous generator models on different accuracy levelsin Section IV. In Section V, it comes to the conclusions.II. P RELIMINARY G ROUNDWORK
This section provides preliminary development on trajec-tory sensitivity, integral inequality, and other necessary back-ground.The trajectory sensitivity of a system of general ordinarydifferential equations (ODEs) is as follows.
Lemma 1. [17] Consider a system of ordinary differentialequations d x d t = f ( t, x, λ ) x ( t ) = x . (1) Suppose f is continuous and has continuous partial deriva-tives with respect to x and λ in the region R = { ( t, x, λ ) : | t − t | ≤ a, k x − x k ≤ b, k λ − λ k ≤ c } . Then the solutionof (1) is continuously differentiable with respect to t in theregion R = { ( t, λ ) : | t − t | ≤ h, k λ − λ k ≤ c } , where h = min( a, bM ) and M = max ( t,x,λ ) ∈ R | f | . Let x s be thesolution of (1) . Let A = ∂f ( t,x s ,λ ) ∂x and B = ∂f ( t,x s ,λ ) ∂λ be thepartial derivatives of f with respect to x and λ , respectively,at ( t, x s , λ ) . Let z = ∂x s ∂t , z = ∂x s ∂x , and z = ∂x s ∂λ be the derivatives, which are also called sensitivities, of thesolution x s with respect to t , x , and λ , respectively. Thenthe sensitivities obey the following ODEs: d z d t = Az , z ( t ) = − f ( t , x , λ ) . (2) d z d t = Az , z ( t ) = 1 . (3) d z d t = Az + B, z ( t ) = 0 . (4) Remark 1.
In Lemma 1, the variable t is limited to a finiterange | t − t | ≤ h. However, if f is uniformly bounded withrespect to t , we can extend R in Lemma 1 to a larger regionwith respect to t . Specifically, if f depends only on x and λ ,we can extend t to infinity in the region R . Remark 2.
Note that the state matrix A = ∂f ( t,x s ,λ ) ∂x in (2) , (3) , and (4) is time-varying and independent of the sensitivities z , z , and z . Therefore, (2) , (3) , and (4) are LTV systems. For a scalar integral inequality, Gr¨onwall’s inequality ap-plies, as follows.
Lemma 2 (Gr¨onwall’s Lemma) . [18] Let α ( t ) , β ( t ) , γ ( t ) bea real continuous function defined in [ a, b ] . If γ ( t ) ≤ α ( t ) + Z ta β ( s ) γ ( s )d s, ∀ t ∈ [ a, b ] and β ( t ) ≥ , ∀ t ∈ [ a, b ] , then γ ( t ) ≤ α ( t ) + Z ta α ( s ) β ( s ) e R ts β ( ρ )d ρ d s, ∀ t ∈ [ a, b ] . (5) If, in addition, the function α is non-decreasing, then γ ( t ) ≤ α ( t ) e R ta β ( s )d s , ∀ t ∈ [ a, b ] . An ODE can be equivalently converted to an integralequation, and then transformed into an integral inequality.Thus, Gr¨onwall’s Lemma can be used to compute bounds ofthe states of the trajectory sensitivity equations (2), (3), and(4).Matrix exponentials are important for linear systems. Wehave the following lemma on the norm of matrix exponentials.
Lemma 3. [19] Let A be a stable matrix, (i.e., the realparts of its eigenvalues are negative), and H be the symmetricand positive-definite matrix such that A T H + HA T = − I , λ max ( H ) represents the maximum eigenvalue of H and λ min ( H ) the minimum one. Then, k e At k ≤ βe − ct , where β = q λ max ( H ) λ min ( H ) and c = λ max ( H ) . From Lemma 3, the norm of the exponential of a stablematrix can be upper-bounded by a vanishing exponentialfunction. Such a result is useful to investigate the steady stateof a linear system.The following lemma may be helpful for investigation of alinear system with vanishing input.
Lemma 4. lim t →∞ Z t e − c ( t − τ ) γ ( τ )d τ = 0 , (6) where c is a positive constant, γ ( t ) ≥ and lim t →∞ γ ( t ) = 0 . The proof of Lemma 4 may be found in a calculus textbook.We re-prove it in the appendix of this paper.III. B
OUNDS OF L INEAR T IME -V ARYING S YSTEMS
Model gaps are caused by imprecise parameters and unmod-eled dynamics. Bounds of model gaps resulting from these twocauses can be captured by analyzing the bound of trajectorysensitivity of the model. Since the trajectory sensitivity of amodel is described by LTV systems, bounds of LTV systemsare derived in this section.Consider an LTV system d z d t = A ( t ) z + u ( t ) z (0) = 0 . (7)The following assumption is imposed on the state matrix ofthe LTV system (7): Assumption 1. A ( ∞ ) = lim t →∞ A ( t ) exists and is a constantmatrix, and the eigenvalues of A ( ∞ ) have negative real parts. emark 3. Assumption 1 is reasonable because practicalpower plants usually reach a steady state, which results inconstant-limit matrices in its trajectory sensitivity.
In this section, we consider two cases: 1) the input to anLTV system is known, and 2) the input to an LTV system isnot known but bounded.The first case can be used to quantify model gaps causedby imprecise parameters and the second one can be applied toobtaining model gaps arising from unmodeled dynamics.
A. Bounds of LTV Systems with Known Input
When the input of the LTV system (7) of which the statematrix has a constant limit as stated in Assumption 1 canbe accurately known, bounds of the state of the LTV systemcan be obtained by decomposing the LTV system into a lineartime-invariant (LTI) system of which the state matrix is A ( ∞ ) as in Assumption 1 and an additional term. The result is statedas follows. Theorem 1.
Consider the LTV system in (7) which satisfiesAssumption 1. Let z LT I ( t ) = e A ( ∞ ) t Z t e − A ( ∞ ) τ u ( τ )d τ, (8) let g ( t ) be an upper bound of Z t k e A ( ∞ )( t − τ ) ( A ( τ ) − A ( ∞ )) z LT I ( τ ) k d τ, and let h ( t, τ ) be an upper bound of k e A ( ∞ )( t − τ ) ( A ( τ ) − A ( ∞ )) k . Then, the bounds of the state of the LTV system (7) can beestimated in the following two ways:Bound 1) The distance between z and z LT I can be boundedas follows: k z ( t ) − z LT I ( t ) k ≤ g ( t ) + max ≤ τ ≤ t { g ( τ ) } (cid:16) e R t h ( t,s )d s − (cid:17) . (9) Thus, the lower bound of the i th entry of z is z LT I,i ( t ) − g ( t ) − max ≤ τ ≤ t { g ( τ ) } (cid:16) e R t h ( t,s )d s − (cid:17) and the upper bound is z LT I,i ( t ) + g ( t ) + max ≤ τ ≤ t { g ( τ ) } (cid:16) e R t h ( t,s )d s − (cid:17) , where z LT I,i is the i th entry of z LT I .Bound 2) The norm of z can be bounded as follows: k z ( t ) k ≤ k z LT I ( t ) k + max ≤ τ ≤ t {k z LT I ( τ ) k} ( e R t h ( t,s )d s − . (10) Thus, the lower bound of the i th entry of z is −k z LT I ( t ) k − max ≤ τ ≤ t {k z LT I ( τ ) k} ( e R t h ( t,s )d s − and the upper bound is k z LT I ( t ) k + max ≤ τ ≤ t {k z LT I ( τ ) k} ( e R t h ( t,s )d s − . Proof.
It can be obtained from (7) that d z d t = A ( t ) z + u ( t )= A ( ∞ ) z + u ( t ) + ( A ( t ) − A ( ∞ )) z. Let β ( t ) = u ( t ) + ( A ( t ) − A ( ∞ )) z. Then d z d t = A ( ∞ ) z + β ( t ) . Then z ( t ) = e A ( ∞ ) t Z t e − A ( ∞ ) τ β ( τ )d τ = e A ( ∞ ) t Z t e − A ( ∞ ) τ u ( τ )d τ + e A ( ∞ ) t Z t e − A ( ∞ ) τ ( A ( τ ) − A ( ∞ )) z ( τ )d τ = z LT I ( t ) + e A ( ∞ ) t Z t e − A ( ∞ ) τ ( A ( τ ) − A ( ∞ )) z ( τ )d τ. (11)Proof of Bound 1): We first derive the upper bound of thedistance between z and z LT I . From (11), it can be deducedthat z ( t ) − z LT I ( t )= e A ( ∞ ) t Z t e − A ( ∞ ) τ ( A ( τ ) − A ( ∞ )) z ( τ )d τ = e A ( ∞ ) t Z t e − A ( ∞ ) τ ( A ( τ ) − A ( ∞ ))( z ( τ ) − z LT I ( τ ))d τ + e A ( ∞ ) t Z t e − A ( ∞ ) τ ( A ( τ ) − A ( ∞ )) z LT I ( τ )d τ (12)Then k z ( t ) − z LT I ( t ) k≤ Z t k e A ( ∞ )( t − τ ) ( A ( τ ) − A ( ∞ )) z LT I ( τ ) k d τ + Z t k e A ( ∞ )( t − τ ) ( A ( τ ) − A ( ∞ )) kk ( z ( τ ) − z LT I ( τ )) k d τ ≤ g ( t ) + Z t h ( t, τ ) k ( z ( τ ) − z LT I ( τ )) k d τ. Applying Lemma 2, we find that k z − z LT I k≤ g ( t ) + Z t g ( τ ) h ( t, τ ) e R tτ h ( t,s )d s d τ ≤ g ( t ) + max ≤ τ ≤ t { g ( τ ) } Z t h ( t, τ ) e R tτ h ( t,s )d s d τ = g ( t ) + max ≤ τ ≤ t { g ( τ ) } e R t h ( t,s )d s Z t h ( t, τ ) e − R τ h ( t,s )d s d τ (13)n (13), Z t h ( t, τ ) e − R τ h ( t,s )d s d τ = Z t e − R τ h ( t,s )d s d( Z τ h ( t, s )d s )= − (cid:16) e − R τ h ( t,s )d s (cid:17) (cid:12)(cid:12)(cid:12) t =1 − e − R t h ( t,s )d s Therefore, k z ( t ) − z LT I ( t ) k ≤ g ( t ) + max ≤ τ ≤ t { g ( τ ) } ( e R t h ( t,s )d s − . Let z i and z LT I,i denote the i th entries of z and z LT I ,respectively. It can be deduced that | z i − z LT I,i | ≤ k z − z LT I k and thus z LT I,i − k z − z LT I k ≤ z i ≤ z LT I,i + k z − z LT I k . Therefore, the lower bound of z i is z LT I,i ( t ) − g ( t ) − max ≤ τ ≤ t { g ( τ ) } (cid:16) e R t h ( t,s )d s − (cid:17) and the upper bound is z LT I,i ( t ) + g ( t ) + max ≤ τ ≤ t { g ( τ ) } (cid:16) e R t h ( t,s )d s − (cid:17) . Proof of Bound 2): We now derive a bound of k z k . Sucha bound can be easily obtained from the derivation of Bound1). However, in this part, we derive a smaller bound of z withfewer steps than that derived from Bound 1).From (11), it can be deduced that k z ( t ) k≤k z LT I ( t ) k + Z t k e A ( ∞ )( t − τ ) ( A ( τ ) − A ( ∞ )) kk z ( τ ) k d τ ≤k z LT I ( t ) k + Z t h ( t, τ ) k z ( τ ) k d τ. From Lemma 2, k z ( t ) k≤k z LT I ( t ) k + Z t k z LT I ( τ ) k h ( t, τ ) e R tτ h ( t,s )d s d τ = k z LT I ( t ) k + e R t h ( t,s )d s Z t k z LT I ( τ ) k h ( t, τ ) e − R τ h ( t,s )d s d τ ≤k z LT I ( t ) k + e R t h ( t,s )d s max ≤ τ ≤ t {k z LT I ( τ )( t ) k}× Z t h ( t, τ ) e − R τ h ( t,s )d s d τ = k z LT I ( t ) k + max ≤ τ ≤ t {k z LT I ( τ ) k} (cid:16) e R t h ( t,s )d s − (cid:17) . Let z i denote the i th entry of z . From the fact that −k z k ≤ z i ≤ k z k , it can be determined that the lower bound of z i is −k z LT I ( t ) k − max ≤ τ ≤ t {k z LT I ( τ ) k} ( e R t h ( t,s )d s − and the upper bound is k z LT I ( t ) k + max ≤ τ ≤ t {k z LT I ( τ ) k} ( e R t h ( t,s )d s − . Remark 4.
A tighter bound can be obtained for Bound 1 viaa more accurate derivation of (13) as follows. k z − z LT I k≤ g ( t ) + Z t g ( τ ) h ( t, τ ) e R tτ h ( t,s )d s d τ ≤ g ( t ) + e R t h ( t,s )d s N − X i =0 max t i ≤ τ ≤ t i +1 { g ( τ ) }× Z t i +1 t i h ( t, τ ) e − R τ h ( t,s )d s d τ, where t = 1 , t N = t , and t i < t i +1 , i = 0 , , · · · , N − .Accordingly, k z ( t ) − z LT I ( t ) k≤ g ( t ) + N − X i =0 max t i ≤ τ ≤ t i +1 { g ( τ ) } e R t h ( t,s )d s × ( e − R ti h ( t,s )d s − e − R ti +10 h ( t,s )d s ) . (14) The bound in (14) is tighter than the bound in Bound 1, butfinding it requires more computation.
With the derived bounds, the behavior of z ( t ) can beapproximated by that of z LT I ( t ) when t is sufficiently large,as shown below. Corollary 1.
Consider an LTV system in (7) with Assump-tion 1. Also, assume that u ( t ) is bounded. Let z LT I be definedas in (8) . Then lim t →∞ ( z ( t ) − z LT I ( t )) = 0 . Proof.
Since u ( t ) is bounded and A ( ∞ ) is stable, z LT I is alsobounded. Let M LT I be a bound of k z LT I k . g ( t ) in Theorem 1can be selected as Z t k e A ( ∞ )( t − τ ) kk ( A ( τ ) − A ( ∞ )) kk z LT I ( τ ) k d τ and h ( t, τ ) as k e A ( ∞ )( t − τ ) kk ( A ( τ ) − A ( ∞ )) k . Then, g ( t ) = Z t h ( t, τ ) k z LT I ( τ ) k d τ ≤ M LT I Z t h ( t, τ )d τ. From Lemma 3, k e A ( ∞ )( t − τ ) k ≤ βe − c ( t − τ ) , where β and c are selected as in Lemma 3. Thus, Z t h ( t, τ )d τ ≤ β Z t e − c ( t − τ ) k ( A ( τ ) − A ( ∞ )) k d τ. Note that lim τ →∞ k ( A ( τ ) − A ( ∞ )) k = 0 . Then, fromLemma 4, lim t →∞ Z t e − c ( t − τ ) k ( A ( τ ) − A ( ∞ )) k d τ = 0 . It can be determined that lim t →∞ Z t h ( t, τ )d τ ≤ . lso, as h ( t, τ ) ≥ and thus R t h ( t, τ )d τ ≥ , lim t →∞ Z t h ( t, τ )d τ = 0 . So, lim t →∞ g ( t ) ≤ . As g ( t ) ≥ , lim t →∞ g ( t ) = 0 . Also, from thefact that lim t →∞ Z t h ( t, τ )d τ = 0 it can be shown that lim t →∞ e R t h ( t,s )d s − . From Bound 1), lim t →∞ ( z ( t ) − z LT I ( t )) = 0 . B. Bounds of LTV Systems with Unknown but Bounded Input
When the disturbance to the ODEs in (1) is constant,the bounds derived in Theorem 1 can be directly applied toquantifying the bound of the trajectory sensitivity of (1). In thissubsection, we thus focus on the case in which the disturbancevaries with time.Let δλ ( t ) be an arbitrary continuous function such that k δλ ( t ) k ∞ = 1 and let η be a small positive number. ηδλ ( t ) can represent an arbitrary continuous signal of which themagnitude is η . The time-varying disturbance of the parameter λ in (1) can be expressed as ηδλ ( t ) . Let λ ( η ) ( t ) = λ + ηδλ ( t ) .Then when η = 0 , λ (0) ( t ) = λ and there is no disturbance.Let x ( η ) be the solution of the following ODE: d x d t = f ( t, x, λ ( η ) ) x ( t ) = x . Then z ( η ) = d x ( η ) ( t )d η is the sensitivity of x ( η ) with respect to η . Then from Lemma 1, d z (0) d t = ∂f∂x z (0) + ∂f∂λ δλ (15)and has initial values z (0) ( t ) = 0 .If we know the mathematical form of δλ , then we canapply the bounds in Theorem 1 to obtain the bounds of thetrajectory sensitivity of (1) when the disturbance varies withtime. However, the mathematical form of δλ is unknown orcannot be accurately known in many cases. Therefore, a boundof LTV systems for bounded disturbances is necessary in thosecases. Such a bound is derived as follows. Theorem 2.
Consider an LTV system (7) with Assumption 1.Also, assume that u ( t ) is a vector of which every entry isa continuous function. Let M u = sup t ∈ [0 , ∞ ] max i | u i ( t ) | ,where u i is the i th entry of u . Let z LT I , g ( t ) , and h ( t, τ ) be defined as in Theorem 1 and a ( t ) be an upper bound of R t k e A ( ∞ )( t − τ ) k d τ . Then, the bound of the state of the LTVsystem can be estimated as follows: k z k ≤ K , ∞ M u (cid:18) a ( t ) + max ≤ τ ≤ t { a ( t ) } ( e R t h ( t,s )d s − (cid:19) , where K , ∞ is a constant such that k u ( t ) k ≤ K , ∞ k u ( t ) k ∞ , ∀ t . (It is easy to see that K , ∞ ≤ √ n , where n is the dimension of u . ) Thus, the lower bound ofevery entry of z is − K , ∞ M u (cid:18) a ( t ) + max ≤ τ ≤ t { a ( t ) } ( e R t h ( t,s )d s − (cid:19) and the upper bound is K , ∞ M u (cid:18) a ( t ) + max ≤ τ ≤ t { a ( t ) } ( e R t h ( t,s )d s − (cid:19) . Proof.
To obtain a bound of k z k with bounded input u ( t ) ,we can first obtain a bound of k z LT I k and then combine thebound of k z LT I k with Bound 2) in Theorem 1. Thus, we nextderive a bound of k z LT I k . From the definition of z LT I in (8), k z LT I ( t ) k ≤ Z t k e A ( ∞ )( t − τ ) kk u ( τ ) k d τ ≤ K , ∞ Z t k e A ( ∞ )( t − τ ) kk u ( τ ) k ∞ d τ ≤ K , ∞ M u Z t k e A ( ∞ )( t − τ ) k d τ ≤ K , ∞ M u a ( t ) . Therefore, k z k ≤ K , ∞ M u (cid:18) a ( t ) + max ≤ τ ≤ t { a ( t ) } ( e R t h ( t,s )d s − (cid:19) . The lower bound and upper bound of every entry of z are − K , ∞ M u (cid:18) a ( t ) + max ≤ τ ≤ t { a ( t ) } ( e R t h ( t,s )d s − (cid:19) and K , ∞ M u (cid:18) a ( t ) + max ≤ τ ≤ t { a ( t ) } ( e R t h ( t,s )d s − (cid:19) , respectively. Remark 5.
The bound in Theorem 2 is derived for all boundedinputs. It might be very large for some types of bounded inputs.
IV. C
ASE S TUDIES
In this section, the bounds obtained in Section III are appliedto quantifying a generator model gap caused by mismatch ofmechanical power via simulation. The bound in Theorem 1 isapplied to quantifying the gap caused by imprecise parameters,while the bound in Theorem 2 is applied to quantifying thegap caused by unmodeled dynamics. Both causes are regardedas disturbances to the model. For a generator, the gap of therotor angle δ is more important than angular frequency ω because it is related to the gap of the electrical power output ofthe generator. In this section, the generator response withoutany disturbance is referred to as nominal response, while thatwith disturbance as disturbed response. Assume a generator ABLE IP
ARAMETER V ALUES IN S IMULATION E XAMPLE
Variable Name Physical Meaning Value Unit M mechanical starting time 13 kWs/kVA f frequency rating 60 Hz x ′ d d-axis transient reactance 0.2 p.u. x ′ q q-axis transient reactance 0.4 p.u. x l leakage reactance 0.15 p.u. V voltage rating 1 p.u. e ′ d d-axis transient potential 0.1 p.u. e ′ q q-axis transient potential 0.9 p.u. r a armature resistance 0.0005 p.u. P m mechanical power 1 p.u. D damping coefficient 100 Time (second) δ (r ad i an s ) -2.4-2.2-2-1.8-1.6-1.4-1.2-1 Time (second) ω (r ad i an s / s e c ond ) Fig. 1. Response of generator in nominal case is represented by the second-order generator model shownbelow: d δ d t = Ω b ( ω − dω d t = 1 M ( P m − P e − D ( ω − P e = ( v q + r a i q ) i q + ( v d + r a i d ) i d Q e = ( v q + r a i q ) i d − ( v d + r a i d ) i q v q + r a i q − e ′ q + ( x ′ d − x l ) i d v d + r a i d − e ′ d − (cid:0) x ′ q − x l (cid:1) i q v d = V cos δv q = V sin δ, The physical meanings, values, and units of these parametersare collected in Table I and Ω b = 2 πf .By solving i d , i q as a function of v d , v q and substituting thelast two equations into it, we can express P e as a function of δ . We denote the function P e = P e ( δ ) . Let P e,δ = d P e d δ . Thegenerator model can be simplified to the following: d δ d t = Ω b ( ω − dω d t = 1 M ( P m − P e ( δ ) − D ( ω − . (16)The nominal response of the generator with parametervalues being collected in Table I is presented in Fig. 1. Theinitial values of δ and ω are selected to be δ = − π − . and ω = 0 . , respectively. The generator model is stable withthis pair of initial values.Let x = ( δ, ω ) T represent the state of the generator and f ( x, P m ) = (Ω b ( ω − , ( P m − P e − D ( ω − /M ) T bethe right-hand side of the generator model. Next, we onlysimulate the model gap caused by uncertainty of P m . If welet z = d x d λ = (cid:0) d δ d P m d ω d P m (cid:1) T , then d z d t = A ( t ) z + f P m ( x, P m ) z (0) = (0 0) T , (17)where A ( t ) = f x ( x ( t ) , P m ) = (cid:18) b − P e,δ M − DM (cid:19) , and f P m = ∂f∂P m = (cid:18) , M (cid:19) T . The trajectory sensitivity equations in (17) are LTV, and itsstate matrix A ( t ) is convergent as δ is convergent, which canbe observed in Fig. 1. When δ ( ∞ ) is substituted into A ( t ) , A ( ∞ ) is stable and the imaginary parts of its eigenvalues arenonzero. As a result, instead of the upper bound of matrix ex-ponentials from Lemma 3, the upper bound of k e A ( ∞ )( t − τ ) k isselected to be C A e − k A ( t − τ ) (1+ S A sin( ω A ( t − τ )+ φ A )) . Notethat the bound obtained from Lemma 3 will result in boundsof a model gap that is too large to make sense. The upperbound of k A ( τ ) − A ( ∞ ) k is selected to be C dA e − k dA ( τ ) (1 + S dA sin( ω dA τ + φ dA )) . The way A ( τ ) − A ( ∞ ) vanishes isdetermined by the way δ converges to zero in (16). Thegenerator model represented by (16) is nonlinear. It is thusdifficult to tell whether δ ( τ ) and A ( τ ) converge in a waysimilar to C dA e − k dA ( τ ) (1 + S dA sin( ω dA τ + φ dA )) . However,we can always find such an upper bound of k A ( τ ) − A ( ∞ ) k .Next, the bounds of LTV systems derived in Section III areapplied to (17) to estimate the model gap caused by impreciseparameters and unmodeled dynamics. It should be noted thatthe bounds of the trajectory sensitivity equations in (17) canonly be applied when the disturbance is small. A. Bounds with imprecise Parameters
The constant disturbance to P m is selected to be . p.u.The difference between the nominal response and the dis-turbed response is shown in Fig. 2, which also showsthe upper bound and the lower bound corresponding tothose in Bound 1) in Theorem 1. Note that the upperbound and the lower bound are expressed as z LT I,i ( t ) ± (cid:16) g ( t ) + max ≤ τ ≤ t { g ( τ ) } ( e R t h ( t,s )d s − (cid:17) and they are notsymmetric about the x-axis.The quantification of the model gap caused by constantdisturbance is quite tight for the phase angle δ . The obtainedbounds have the same steady-state value as the actual modelgap, and the upper bound is very close to the actual model gapat the beginning. In the middle of the response, when the timeis around . second, there is mismatch between the obtainedbounds and the actual response gap. That is because the upperbound estimates of k e A ( ∞ )( t − τ ) k and k A ( τ ) − A ( ∞ ) k that are ime (second) ∆ δ (r ad i an s ) -0.01-0.00500.0050.010.015 upper bound ∆ δ lower boundTime (second) ∆ ω (r ad i an s / s e c ond ) -0.01-0.00500.0050.01 upper bound ∆ ω lower bound Fig. 2. Gap for generator with constant disturbance
Time (second) ∆ δ (r ad i an s ) -0.0500.05 upper bound ∆ δ lower boundTime (second) ∆ ω (r ad i an s / s e c ond ) -0.0500.05 upper bound ∆ ω lower bound Fig. 3. Gap for generator with a sine wave disturbance used in the simulation are larger than the actual k e A ( ∞ )( t − τ ) k and k A ( τ ) − A ( ∞ ) k . The quantification of the model gapobtained here is quite accurate for the phase angle δ .Note that the bounds of the model gap are obtained for thestate vector of the system, which is (cid:0) δ ω (cid:1) T in this example.If one state plays a major role in the gap, the obtained boundsmay be tight for this state and very loose for other states. Thisis why the quantification is not accurate for ω , which can beobserved in the plot at the bottom of Fig. 2. However, becausethe phase angle δ determines the power output of the generatorand is more important than ω , the quantification of model gapsworks well enough for the generator model. B. Bounds with Unmodeled Dynamics
In this subsection, two time-varying disturbance signals areadopted, representing unmodeled dynamics. The first one is asine wave while the second one is the difference between thegenerator model in (16) and a more complex one.The magnitude of the sine wave is selected to be . p.u.and the frequency to be the resonant frequency of the system, ˙ x = A ( ∞ ) x . The disturbance signal at the resonant frequencyis amplified to the largest degree by the LTI system of whichthe state matrix is A ( ∞ ) . Combined with Corollary 1, whichstates that the behavior of the LTV system of sensitivity Time (second) ∆ δ (r ad i an s ) -0.03-0.02-0.0100.010.020.03 upper bound ∆ δ lower boundTime (second) ∆ ω (r ad i an s / s e c ond ) -0.03-0.02-0.0100.010.020.03 upper bound ∆ ω lower bound Fig. 4. Gap between second-order model and higher-order model equations in (17) satisfying Assumption 1 can be approximatedby the LTI system ˙ z = A ( ∞ ) z + f P m , a sine signal at thegenerator’s resonant frequency results in a larger model gapthan other signals with the same magnitude. The model gapcaused by this sine signal is presented in Figure 3, in whichthe upper bound and the lower bound are those in Theorem 2and are symmetric about the x-axis. The bound is close forthe disturbance of δ , which shows that the bounds of modelgaps are tight for this case.The second signal for unmodeled dynamics originates fromthe difference between (16) and a higher-order model. Thehigher-order model of the generator has the following defini-tions in addition to (16), T ⋆in = T order + 1 R ( ω ref − ω ) T in = T ⋆in , T min ≤ T ⋆in ≤ T max ,T max , T ⋆in > T max ,T min , T ⋆in < T min , ˙ t g = T in − t g T s , ˙ t g = (1 − T T c ) t g − t g T c , ˙ t g = (1 − T T )( t g + T T c t g ) − t g T ,P m = t g + T T ( t g + T T c t g ) , (18)where the physical meanings and values of the variables arepresented in Table II. Note that P m in (16) is regarded as aconstant, while it is dynamic in (18). P m in (18) is an outputof an LTI system with saturated input. The gap between thehigher-order model and the second-order model is presented inFig. 4, in which the upper bound and the lower bound are thosein Theorem 2 and are symmetric about the x-axis. It can beobserved that the gap between the second-order model in (16)and the higher-order model is within the bound developed inTheorem 2.It should be noted that the bounds used in this subsection aredeveloped for model gaps caused by all unmodeled dynamics. ABLE IIP
ARAMETER V ALUES IN H IGHER - ORDER M ODEL
Variable Name Physical Meaning Value Unit ω ref reference speed 1 p.u. R droop 0.02 p.u. T max maximum turbine output 1.2 p.u. T min minimum turbine output 0.3 p.u. T s governor time constant 0.1 s T c servo time constant 0.45 s T transient gain time constant 0 s T power fraction time constant 12 s T reheat time constant 50 s Thus, the bounds may be very loose for some specific distur-bance, such as that resulting from a higher-order model, aspresented in Fig. 4. But results for the sine signal selected forFig. 3 show that the developed bounds are tight for the classof bounded dynamics.V. C
ONCLUSIONS
Bounds of model gaps were quantified by estimating thebounds of trajectory sensitivity. First, bounds for responses ofgeneral linear time-varying systems were obtained with bothconstant input and dynamic input. Then, the bounds for lineartime-varying systems were applied to quantifying the modelgaps. Simulations with a generator model showed the efficacyof the quantification. Future work may include investigationof more general models and online quantification.A
PPENDIX P ROOF OF L EMMA γ ( t ) → , k γ k ∞ = max ≤ t< ∞ { γ ( t ) } is finite. Andfor any ǫ > , there exists T γ , such that for all τ > T γ , γ ( τ ) < c k γ k ∞ +1 ǫ ; there also exists T e , such that for all t e >T e , e − ct < c k γ k ∞ +1 ǫ . Let T = max { T γ , T e } . Then, for all t > T , Z t e − c ( t − τ ) γ ( τ )d τ = Z T e − c ( t − τ ) γ ( τ )d τ + Z tT e − c ( t − τ ) γ ( τ )d τ ≤ max ≤ τ ≤ T { γ ( τ ) } Z T e − c ( t − τ ) d τ + c k γ k ∞ + 1 ǫ Z tT e − c ( t − τ ) d τ = max ≤ τ ≤ T { γ ( τ ) } e − c ( t − T ) − e − ct c + c k γ k ∞ + 1 ǫ − e − c ( t − T ) c< max ≤ τ ≤ T { γ ( τ ) } e − c ( t − T )+ e − ct c + c k γ k ∞ + 1 ǫ − e − c ( t − T ) c< k γ k ∞ c c k γ k ∞ + 1 ǫ + c k γ k ∞ + 1 1 c ǫ = ǫ. As a result, lim t →∞ R t e − c ( t − τ ) γ ( τ )d τ = 0 . R EFERENCES[1] D. Kosterev, “Hydro turbine-governor model validation inPacific Northwest,”
IEEE Transactions on Power Systems ,vol. 19, no. 2, pp. 1144–1149, 2004. [Online]. Available:https://ieeexplore.ieee.org/document/1295026 [2] Z. Huang, P. Du, D. Kosterev, and S. Yang, “Generatordynamic model validation and parameter calibration using phasormeasurements at the point of connection,”
IEEE Transactions on PowerSystems , vol. 28, no. 2, pp. 1939–1949, 2013. [Online]. Available:https://ieeexplore.ieee.org/document/6487426[3] North American SynchroPhasor Initiative (NASPI), “Technicalworkshop - model verification tools,” 2016. [Online]. Available:https://naspi.org/node/528[4] P. Pourbeik, R. Rhinier, S. Hsu, B. L. Agrawal, and R. Bisbee,“Semiautomated model validation of power plant equipment usingonline measurements,”
IEEE Transactions on Energy Conversion , pp. 1–5, 2016. [Online]. Available:https://ieeexplore.ieee.org/abstract/document/7519886[8] C. Tsai, L. Changchien, I. Chen, C. Lin, W. Lee, C. Wu,and H. Lan, “Practical considerations to calibrate generator modelparameters using phasor measurements,”
IEEE Transactions on SmartGrid , vol. 8, no. 5, pp. 2228–2238, 2017. [Online]. Available:https://ieeexplore.ieee.org/document/7401126[9] J. R. Hockenberry and B. C. Lesieutre, “Evaluation of uncertainty indynamic simulations of power system models: The probabilisticcollocation method,”
IEEE Transactions on Power Systems ,vol. 19, no. 3, pp. 1483–1491, 2004. [Online]. Available:https://ieeexplore.ieee.org/document/1318685[10] A. R. Borden and B. C. Lesieutre, “Determining power system modalcontent of data motivated by normal forms,” in , 2012, pp. 1–6. [Online].Available: https://ieeexplore.ieee.org/document/6336387[11] A. R. Borden, B. C. Lesieutre, and J. Gronquist, “Power system modalanalysis tool developed for industry use,” in , 2013, pp. 1–6. [Online].Available: https://ieeexplore.ieee.org/document/6666956[12] A. R. Borden and B. C. Lesieutre, “Variable projection method forpower system modal identification,”
IEEE Transactions on PowerSystems , vol. 29, no. 6, pp. 2613–2620, 2014. [Online]. Available:https://ieeexplore.ieee.org/document/1318685[13] R. Huang, E. Farantatos, G. J. Cokkinides, and A. P. Meliopoulos,“Physical parameters identification of synchronous generators bya dynamic state estimator,” , pp. 1–5, 2013. [Online]. Available:https://ieeexplore.ieee.org/document/6672882[14] A. A. Hajnoroozi, F. Aminifar, and H. Ayoubzadeh, “Generating unitmodel validation and calibration through synchrophasor measurements,”
IEEE Transactions on Smart Grid , vol. 6, no. 1, pp. 441–449, 2015.[Online]. Available: https://ieeexplore.ieee.org/document/696561[15] R. Huang, R. Diao, Y. Li, J. J. Sanchez-Gasca, Z. Huang,B. Thomas, P. V. Etingov, S. Kincic, S. Wang, R. Fan, et al. ,“Calibrating parameters of power system stability models usingadvanced ensemble Kalman filter,”
IEEE Transactions on PowerSystems , vol. 33, no. 3, pp. 2895–2905, 2018. [Online]. Available:https://ieeexplore.ieee.org/document/8068956.[16] I. A. Hiskens and J. Alseddiqui, “Sensitivity, approximation, anduncertainty in power system dynamic simulation,”
IEEE Transactionson Power Systems , vol. 21, no. 4, pp. 1808–1820, 2006. [Online].Available: https://ieeexplore.ieee.org/document/1717585[17] T. Ding and C. Li,
A Course on Ordinary Differential Equations .Beijing: Higher Education Press, 2004.[18] S. S. Dragomir,
Some Gronwall type inequalities andapplications . Nova Science, 2003. [Online]. Available:https://papers.ssrn.com/sol3/papers.cfm?abstract id=3158353[19] G. Hu and M. Liu, “The weighted logarithmic matrix normand bounds of the matrix exponential,”