State Estimation in Power Distribution Systems Based on Ensemble Kalman Filtering
State estimation in power distribution systems is a key component for increased reliability and optimal system performance. Well understood in transmission systems, state estimation is now an area of active research in distribution networks. While se…
Authors: C. Carquex, C. Rosenberg, K. Bhattacharya
1 State Estimation in Po wer Distrib ution Sys tems Based on Ensemble Kalman Filtering C ˆ ome Carquex, Catherine Rosenberg and Kankar Bhattach arya Abstract —State estimation in power distribution systems is a k ey component fo r increased reliability and optimal system perfo rmance. W ell understood in tra nsmission systems, state estimation is now an area of activ e resear ch in di st ribut i on networks. While sev eral sn apshot-based approaches have been used to solve this problem, few solutions ha ve been propose d in a dynamic framewo rk. In this paper , a Past-A ware State Estimation (P ASE) method is p roposed for distribution systems th at takes prev ious estimates into account to improv e th e accuracy of the current one, using an Ensembl e Kalman Fi lter . Fewer phasor measureme nts unit s (PMU) are needed to achieve the same estimation error tar get than snapshot-based methods. Contrary to current methods, the proposed solution does not embed power flow equ ations into the estimator . A th eoretical fo rmulation i s presented to compute a p riori t he advantag es of th e p roposed method vis-a-vis the state-of-th e-art. The proposed approach is validated considerin g the 33-bus distribution system and using power co nsumption traces from real h ou seholds. Index T erms —Distribution system, distribution system state estimation, ensemble Kalman fi lter , ph asor measurem ent un it, state estimation I . I N T R O D U C T I O N T RADITION ALL Y, electric power d istribution systems have been designed an d operated as passive systems to meet the cu stomers’ dema nd. Howev er , with tr ansformatio n of the grid to a smart grid , the reliability an d operatio nal challenges o f distribution systems have incre a sed. An o p erator will n eed to m anage the distribution system mo re closely in the futur e, re quiring impr oved visibility of its states [1] which will inv olve real-time monitoring [2]. Indeed, m ost so lutions to smart g rid r elated challenges at the distribution level assume a knowledge of the states of the system, and th erefore e ssentially rely on Distrib ution System State Estimation (DSSE), which is a key fu nction of supervisory control that some utilities h av e already began rolling -out [3]. The state of a power system can b e comp letely defined from the k n owledge of all bus voltage m agnitudes and an gles at time t [4]; ty pically , state estimation is carried ou t based on measuremen ts of variables such as the voltage m agnitudes and angles, av ailable from Phasor Measurement Units ( PMUs). State estimation of power systems is a well und erstood problem at the tran smission le vel a nd is traditiona lly solved using a sn apshot-based weighted least square (WLS) meth od which relies o n hig h quality measur ement data from PMUs [4]. Howe ver , transmission systems g enerally h av e a lim ited number o f buses and are equ ip ped with many me a surement C. Ca rquex , C. Rosenbe rg and K. Bha ttacha rya are wit h the Department of Electric al and Computer E nginee ring, Unive rsity of W at erloo, W aterloo, ON N2L 3G1, Canada, e-mail: { cacarq ue, cath, kankar } @uw aterlo o.ca devices since it is im portant to precisely monitor and contro l the system at all times. On the other hand, distribution systems comprise a large num ber of buses with little to n o measure- ments av ailable. While se veral re cent stud ies h av e focused on developing lo w-cost, easy to dep loy PMUs [5], [6], it is not practical to install PMUs at every distribution bus. I f PMUs were to be placed at selected buses only , there would be infinitely many solution s to the DSSE pro b lem. In ord er to reduce the number of p ossible solutions, p seu do-measu rements can be used [7], which are load forecasts compu ted ah e ad of time to aid DSSE in finding a “go od” solutio n. T ypic a lly , a pseudo- m easuremen t at a gi ven load b us co mprises an estimate of the expected active and rea cti ve power co nsumption s at the bus. Loa d f orecasting at the distribution le vel is difficult, hence pseudo-m easurements are usually of poo r quality . These fundam ental difference s, and the need for affordab le solutio n s, mean that new state estimation a p proach es are needed fo r distribution systems. Many studies have extended the WLS ap proach from tr ans- mission to distribution systems. A revie w of literature on the different state estimation techn iques and their applicatio n to DSSE problems is presented in [8]. One of th e first applications of the snap shot a p proach to the DSSE pr o blem was r eported in [9], where a p robab ilistic fo rmulation b ased on pseudo- m easuremen ts was used. In [10], th e power-flo w equa- tions were lin earized and a compu tationally friendly solution method was pro posed. The authors also showed that PMUs are needed for accu rate state estimatio n . Compressed sensing theory was used for state estimation with sparse measureme n ts in [1 1], while [12] used line-curr ent m a gnitudes and angles. Finally a sem i- definite pr ogramm in g a pproach was used to solve the DSSE problem in [13]. Sev eral resear c h ers have used Kalm an filters in state es- timation problem s fo r tr ansmission systems [14]. Howe ver , in distribution systems, the p oor quality of the pseudo - measuremen ts renders such m ethods ineffectiv e. Theref ore, very fe w Kalman filtering b ased methods have been de veloped for DSSE and none improve over the WLS. Huan g e t al . compare d the extended Kalman filter to the u nscented Kalman filter in [15]. From the r eported results it was n oted that th e re was no visible impr ovement in perform ance of the Kalman filter based metho ds over WLS. In [16] the impact of choice of th e model and measuremen t covariance matr ix on the perfor mance of the extended Kalman filter was examined. It was noted fro m the results that th e proposed filtering ap proach did not result in any perfo rmance impr ovement. The above discussed Kalman filter based app roaches ap ply the m ethods directly from the tran smission to distribution systems. Th e 2 problem of po or quality of pseud o-measure m ents is alleviated by assumin g that me a su rements are av ailable at every bus in real-time or quasi-re al-time, usually from synchr o nized smart- meters, which is not rea listic. In this paper, a p ast-aware method for DSSE, nam ed P ASE (Past-A ware State Estimation) , where the estimate a t time t depend s o n anterio r estimates a n d based on the En semble Kalman Filter (E nKF) [1 7] is presen te d . Applyin g the EnKF to this pr oblem is non -trivial, since m e asurements from sou rces with different time-scales must be m erged. Contrary to WLS and other ap proache s using different variations of the Kalm an filter , P ASE does not em bed the power flow equation s into the estimator, making it a versatile technique. Instead it relies on an external p ower -flow solver , which is left to the cho ice of the operator . In a sn apshot-based context where the state a t time t is computed in depend ently of the estimates at tim e s an terior to t , the WLS ob jectiv e function pr ovid es the best p erforma nce possible (exclud in g ill-condition ed cases) [1 8]. Such an esti- mator is referred to as the State of the Art ( SoA) in this pape r , for the p urpose of comparison. Specifically , the c ontributions of the work are threefold : 1) A maiden attempt is made to apply EnKF to a d istribution system spar sely mo nitored by PMUs f o r state estimatio n ; 2) An analytical fr amew ork is dev eloped to ev aluate the perfor mance of P ASE; 3) The theo retical results ar e validated via extensive simulation s on a 33 - bus distribution system and using p ower consump tio n tr aces f rom r eal hou seholds. The perfor mances of the proposed P ASE ap p roach and WLS are compare d and eng ineering insigh ts ar e provided to u nderstand the im pact of each decision variable on the perf o rmance of P ASE, as well as the trade- offs to m ake. Based on the above discussions, the main message of this work is th at P ASE is the first tech nique to imp rove upo n the SoA. It does so significantly when the elapsed time betwee n two consecu ti ve state com putations is small (less than 1 5 minu tes), i.e., less PMUs are needed to ach iev e the same estimation error . The rest of the paper is o rganized as follows. The b a ck- groun d and assum p tions ar e presented in Section I I. The SoA method is p resented in Section III and the pro posed P ASE solution in Section IV. Th e validation results are reported in Section V. Finally , the con clusions are d rawn in Section VI. I I . S Y S T E M A N D A S S U M P T I O N S The assumptions ar e stated in this section . A three-p hase balanced distribution system unde r nor m al operations is con - sidered. Th e DSSE p roblem is solved by the loc al distribution company ( LDC) using an app ropriate compu ta tio nal platfo rm. The following inform ation is need e d to implem ent DSSE, both with the SoA m ethod and th e p roposed P ASE method. Computational t imescale : a new state estimate is com puted ev ery ∆ T . T ypically in transmission systems, a tim e -step of 5 to 15 min is c o nsidered. In d istribution systems smaller time- steps are nee d ed because of higher load volatility , which can arise for example with high pen etration o f r e new ables. T he value o f ∆ T has an imp act on the co m putational burden. In this work time-steps f rom 6 secon d s to 15 minu tes are considered . Altog ether, the choice of an a p propr iate timescale for DSSE problems is still an open question. T opology : the distribution system has a radial topolog y and is d efined by a set of buses I o f car dinality | I | as well as a set of branch es B of co nstant an d known impe d ances, con necting the buses. The substation tra nsformer is modeled as a referen ce voltage sou rce of ma g nitude V 0 . Measurements: the subset S ⊆ I of buses are eq uipped with PMUs that m o nitor every ∆ T b oth th e bus voltage m ag- nitudes ( V s ) and bus angles ( δ s ). The measureme nts repo rted by the PMUs are a ssum ed to be unb iased an d the variance of the error of the reading s is known. These assumption s are common ly made in state estimation works [4]. A br oadban d commun ication infrastru c ture is available to transmit the m e a- surements with low latency and h igh reliability . The PMUs are placed in the distribution system according to a giv en mapping S . Pseudo-measurements: th ese are forecasts that “measure” both active and reactiv e powers. They are available fo r each bus i in I . Forecasts are ma d e at perio dic intervals ∆ T ′ , typically on ce a d ay fo r the n ext day (day -ahead fore c a st). At the time of compu tation, th e most recent forecast is used. Clearly , forecasts and PMU mea su rements are on completely different time-scales ( ∆ T ′ ≫ ∆ T ), hence the non - triviality of the E n KF . Forecasts ar e mad e b ased on historical data. Previous estimation work b ased on Kalman filters assumed real-time consumptio n data. This strong requirement is relaxed with forecasts. Data requirements: both the SoA and P ASE approach es require a forecasting method as well as samp le power con - sumption traces (active an d reactiv e) from the system at the lev el of each distribution transfo rmer, from which th e forecasting m ethod can be calibrated . Using th e d a ta, erro r parameters can be ob tained offline. L et e i ( t ) b e the forecast error at bus i and time t (for activ e p ower , f o r exam ple); e i ( t ) is assumed to b e a stationary rando m proc e ss. More over, forecasts are assumed to be unb iased (E [ e i ( t )] = 0 ) an d the variance of the errors (E [ e i ( t ) 2 ] ) to b e known. The estimation of the variance of the for ecast erro rs comes fr om the acquire d data. The assump tion o f an un biased fo recast is a strong hypoth esis, althou gh it is alm ost always used b y r esearchers [10]. The pro posed P ASE metho d n eeds two additiona l info r- mation that can be d eriv ed from the sam e sample d ata: a load ev olution model (wh ic h will be discu ssed in Section IV -A) and the forecast erro r correlation coefficient, ev aluated between two (compu ta tio n) time-steps at a given bus (i.e. , E [ e i ( t ) e i ( t − ∆ T )] ). Giv en th at the data samples ar e ne e ded for both meth ods, not mu ch work is inv olved to der i ve the se additional quantities from it. Finally , th e load f orecast error s ar e assumed to be un corre- lated b e twe en buses, an assumption often made in the literature [10]. System state: it is r epresented by state vectors; different (equiv alent) state representation s may be used depe n ding on their ease of use in the problem formulation . For examp le, y [ t ] = [ V [ t ] T , δ [ t ] T ] T 3 is a possible state vector rep resentation, where V [ t ] is the vector of voltage magn itudes at each bus, an d δ [ t ] the vec- tor of voltage ang les. Another way is to define x [ t ] = [ P [ t ] T , Q [ t ] T ] T where P [ t ] a n d Q [ t ] den o tes the vectors o f activ e a nd re activ e p ower injections a t each bus, resp e cti vely . Note that the power-flo w equatio ns link the state-vector s x and y . A third way , used in theo retical formulation s, is w [ t ] = [ v 1 [ t ] , . . . , v | I | [ t ]] T where v i [ t ] is the voltage ph asor at bus i , time t ; this can also similarly be related to oth er representatio ns. Limitations: In this work, unbalan ced system, distributed generation and b iased measurem ents ar e not consider e d an d are left for f uture studies. I I I . S TA T E - O F - T H E - A RT D S S E M E T H O D The SoA metho d [ 4] used to solve the DSSE pro blem is a sna p shot app roach an d uses a nonlinear WLS objective function . Given the system characterized by th e sets I , B , S and the mapp in g S , th e system state, at a given time, is estimated using an overdeterm ined set of equations. In the following, th e tim e depen dency of the variables is dr opped for better readab ility . T h e variables to be determ ined are the 2 | I | state variables. Each measure m ent add s o ne constraint. There are either 2 o r 4 measu rements per bus (active/reacti ve power forecast, voltage magnitude, and angle), depen ding on whether th ere is a PMU at the bus. Th e num ber of con straints is given b y M = 2 | I | + 2 | S | . Th e PMU measurem e nts and the fo recasts are stored in a vector z of length M , an d are relate d to the system state as per th e following mode l: z = f ( y ) + η wher e f is the f u nction that maps the state vector to the measur e ment vector, an d η is th e vector co ntaining the noise term o f each measurem ent. For example, f ( y ) = [ V ( y ) T , δ ( y ) T , P ( y ) T , Q ( y ) T ] T where V ( y ) an d δ ( y ) are the vectors, r espectiv ely , co ntaining the voltage mag nitude an d angle measurem ents at the buses with PMUs and P ( y ) and Q ( y ) are vectors o f active and reactive power fo r ecasts of size | I | , respectively . Assuming that th e measur ement errors are uncorr elated an d ha ve zero mean, the covariance m atrix Σ of the erro r vector η is written as, Σ = d iag ( σ 2 1 , ..., σ 2 M ) , where σ 2 m is the variance of the m th measuremen t. The ob jectiv e function to b e minimized at eac h time-step is g iven below: J ( y ) = ( z − f ( y )) T Σ − 1 ( z − f ( y )) (1) Sev eral meth ods exist to minimize the objec ti ve functio n, the simplest be in g to iteratively linearize f and solve the r esulting objective u sing the normal eq uations. I V . P RO P O S E D M E T H O D : P A S E T o solve the DSSE pr oblem, P ASE, an En KF-based method, is proposed . Kalm an filters are sequential filtering m ethods. Each iteration is a two steps pro cess: 1 ) the system state is integrated in time using an ev olution mo del, defining an ( a priori ) state estimate. 2) A vailable m e asurements ( in cluding pseudo- m easuremen ts) are used to correct th e estimate and define th e up dated state. T he term assimilation is used to refer to the seco nd step. The lo ad ev olution mo del used in this app roach is presented in Section IV -A. The idea behin d the p roposed ap proach is simple: the additional informatio n provided by the loa d ev olution model and th e pr eviously estimated states are u sed to alle viate th e poo r q uality o f pseudo- m easuremen ts. A. Lo ad Evolution Mod el For ea ch distribution tran sformer bus, an evolution model for the aggregate load is needed, b oth fo r the active an d reactive power co n sumption s. Specifically , the lo ad variation between two (comp utation) time-steps is consider ed: let L p i ( t ) and L q i ( t ) de n ote the instantaneous active and reacti ve ag g re- gated power respectively , at bus i and time t . It is assumed that L p i ( t ) and L q i ( t ) are stationar y ran dom processes. The load variation (aka load ev olution mod el) for active and reactive powers are defined as th e stationary ran dom processes L p i ( t ) − L p i ( t − ∆ T ) and L q i ( t ) − L q i ( t − ∆ T ) r espectively , characterized by th eir p robability d ensity fu n ctions (pdf ). The m e an of th e pro cesses is ze r o and the variance of the processes can be comp uted from the pd f b o th for active and reactive powers at b us i , denoted ( σ p i ) 2 and ( σ q i ) 2 , respecti vely . Such an ev olution mo d el is simple and fits with in th e EnKF framework. The pd f can be derived empirica lly , for examp le, from th e existing req uired samp le traces, d iscu ssed in Section II as will be explained later . Clearly a given lo ad ev olution model is valid only for systems with similar load com p ositions, and will v ary for different geographical ar eas. B. E nsemble Kalman F ilter The tradition al Kalman filter main tains a covariance m atrix associated with the state estima te . Th e EnKF does not use such a matr ix an d r epresents the system state pd f using a set of state vectors called e nsemble. Such ensemb le at time- step k (i.e., time k ∆ T ) is named X k . The covariance matrix is rep laced by the em pirical covariance computed fro m the ensemble. Th e estimate d system state is simp ly the mean of the ensemble co lumns. The size of the ensemble , L , will impact perf ormance. A small ensem b le size will yield faster compu ta tio ns. Howev er the c ovariance estimate fr om the ensemble will be less accurate. Th erefore there is a trad e - off between com putational speed and accu racy and a typical choice is a size o f L = 50 0 or 1 000 . T he covariance estimato r cov ( A, B ) of two ensembles A, B is defined as [17]: cov ( A, B ) = 1 N − 1 ( A − E [ A ])( B − E [ B ]) T (2) where E [ A ] is the mean of the column vecto rs contained in ensemble A. For cov ( A, A ) th e shorter sy ntax cov ( A ) is u sed. Each iteration of the EnKF (c orrespon ding to a computation of the state vector at time-step k ) follows the pr ocedur e detailed in Algorith m 1, each steps of the algo rithm are discussed next. C. In itial E nsemble The state vector x = [ P T , Q T ] T (of size 2 | I | ) is used . It is chosen giv en that the load e volution model d escribed in Section IV - A is defin ed in terms of in jected power . The pdf of the state vector x is represented by an ensem ble of size L : X 0 = [ x 0 1 , . . . , x 0 L ] , X 0 is a 2 | I | × L ma tr ix containing the ensemble m embers. The initial en semble is b uilt by ch oosing a “best-guess” estimate x 0 of th e state vector, to 4 Algorithm 1 Estimation of the state at time-step k Input: X k − 1 , m easurements and p seudo-m easurements at time-step k . 1: Com pute X k p : integrate the ensemble in time (Eq. 3) 2: Com pute X k u : assimilate pseudo-measure m ents (E q. 11) 3: Com pute X k a : assimilate PMU me a surements ( Eq. 13) 4: X k ← X k u Output: Estimated state ˜ x k = E [ X k ] for time-step k . which perturbatio ns are added to represent the error statistics of the in itial gu ess. The erro r distribution chosen for the initial ensemble is discussed in Section V. D. E nsemble Inte gration The EnKF is considere d at time-step k . Th e prior ensemble X k p is ob ta in ed by individually integrating fo r ward in tim e each vector of the ensemble X k − 1 , which was co m puted at the previous time-step. The integration is such that: X k p = X k − 1 + [ n 1 , . . . , n L ] (3) where n l ( l = 1 , . . . , L ) a re column vectors o f size 2 | I | containing th e stochastic noise wh ich accounts fo r th e un- certainties of the load ev olution mo del. Based o n the lo ad ev olution m odel defined in Section IV -A, two variance values ( σ p i ) 2 and ( σ q i ) 2 are associated to each bus i ( i = 1 , . . . , | I | ), respectively for the active and reactiv e powers. The ir values depend o n th e emp irical pdf der i ved. Each n i,l and n | I | + i,l ( i = 1 , . . . , | I | ) is respectively d rawn fro m a distrib ution which represents the empirical pd f o f the load evolution model. Note that the EnKF can accep t a ny load ev olution mo del. E. A ssimilation of Pseudo-Measurements The assimilation of measurements and pseu do- measuremen ts correspo nd to the u pdate step of the Kalman filter , d escribed at the beginn ing of Section IV. An assumption in Kalman filtering is that th e measur e m ent error is white Gaussian no ise. Sinc e p seudo-me asurements are for ecasts a n d do n o t depend on the state of the sy stem , they do not satisfy th is requir ement; instead the fo r ecast error is correlated in time. This pro blem, which is recurrent in Kalman-ba sed k inematic GPS app lications has b een solved previously , and a summar y of the different existing techniqu es can b e found in [19]. The solu tio n chosen in this paper is the time-differencing appr oach describ ed in [20] to remove time- correlated er ror in the pseudo-m easurements. T his me th od was selected for two rea so ns: 1) it doe s n ot requir e any reinterpr e tation of the Kalman equation s and 2) it do es no t introdu c e any laten cy . T o remove the c o rrelations, the f ollowing pro cess is used. Let th e transition matr ix Ψ o f the time-corr elated err or be defined as: Ψ = diag ( ψ p 1 , . . . , ψ p | I | , ψ q 1 , . . . , ψ q | I | ) (4) where ψ p i and ψ q i ( i = 1 , . . . , | I | ) are th e forecast error correlation coefficients at bus i , respectively for active and reactive po wers, introduc e d in Section II; Ψ is d iagonal since the forecast errors b etween b uses are assumed to be uncorre- lated. Q is define d as th e m odel noise cov ariance matrix, and is gi ven as: Q = diag (( σ p 1 ) 2 , . . . , ( σ p | I | ) 2 , ( σ q | I | +1 ) 2 , . . . , ( σ q 2 | I | ) 2 ) (5) R is th e covariance m atrix of the forecast er ror, of size 2 | I | × 2 | I | . R is d iagonal sinc e the fo r ecast error s are a ssum ed not correlated across buses, and is gi ven as: R = diag (( σ f p 1 ) 2 , . . . , ( σ f p | I | ) 2 , ( σ f q | I | +1 ) 2 , . . . , ( σ f q 2 | I | ) 2 ) (6) where σ f p i and σ f q i are th e standa r d d eviations of the forec a st error at bus i , respecti vely for the acti ve and reactiv e powers. The pseud o mea su rements are containe d in a vector d of size 2 | I | . An ensemble D of L perturb ed observations is defined such that D = [ d 1 , . . . , d L ] with each d l = d + ǫ l ( l = 1 , . . . , L ), where ǫ l is a vector drawn from a distribution which models the p seudo-m easurement noise. Before establishing th e update step, intermediary matrices ar e de fin ed next, which will be reused for th e th eoretical deriv ations. H ∗ = H − Ψ H , C = QH T Ψ T , D ∗ = D − Ψ D (7) R ∗ = ( R − Ψ R Ψ T ) + Ψ H QH T Ψ T (8) The updated observation and measurement matrices ( H ∗ and D ∗ , respectively , ar e compu ted in (7). The upd ated measure- ment error matrix R ∗ is comp uted in (8); Ψ is used to r emove the time correlation of the for ecast error between two time- steps. The mo del no ise matrix Q is ne eded to ensure that th e noise in troduce d by the ev olution step is retain ed. Indeed such noise doe s not have any time co rrelation com ponent. In this context, the observation matrix H is the identity ma tr ix ( in Section IV -G the ob servation matrix will n ot be the same) . T he update e q uations for the assimilation of pseudo -measureme n ts are gi ven as: E = H ∗ cov ( X k p ) H ∗ T + R ∗ + H ∗ C + C T H ∗ T (9) K = ( cov ( X k p ) H ∗ T + C ) E − 1 (10) X k u = X k p + K ( D ∗ − H ∗ X k p ) (11) F . Assimilation of PMU Measur ements Similar to the pseud o-measur ements, th e m easurements coming fr o m the PMUs are co ntained in a vector z of size 2 | S | . An ensemb le Z of L p erturbed o bservation vectors is computed su c h that Z = [ z 1 , . . . , z L ] , with each z l = z + ξ l ( l = 1 , . . . , L ) , where ξ k is a vecto r drawn fro m a distribution which models the mea su rement noise. The measure m ents f rom the PMUs can be related to the state vector using a fun ction h , such that z l = h ( x l ) + γ k , where γ k is an error vector . T he fu nction h ( · ) takes as input th e system state an d return s a vector co n taining the measurem ents that would hav e been observed considerin g that particu lar system state. Given that x contain s the active and reactive powers injected at each bus, h ( · ) is the p ower -flow solution ; the EnKF does not ne e d to know the analytical expre ssion of h ( · ) . I t is the solutio n gi ven b y the LDC’ s po wer-flo w solver , for example. This m akes the EnKF indep endent o f the way 5 power -flows are comp uted. The cost of such indepe n dence is co mputation al: one n eed to compute L p ower-flo ws at each time-step. Since h ( · ) is n o n-linear, th e measureme n ts cannot be obtain e d directly from the state using a simple multiplication by an observation m a trix. Instead, h ( x ) needs to b e compu ted explicitly . A tempo rary augmented state b x and augmen te d en semble b X k u are u sed to per form the assimilation, where: b x l = [ x l T , h T ( x l )] T , b X k u = [ b x 1 , . . . , b x L ] (12) The updated ensemble X k a is then c o mputed: X k a = X k u + K ( Z − b H b X k u ) (13) K = cov ( X k u , b H b X k u )[ cov ( b H b X k u ) + cov ( Z )] − 1 (14) where b H is a selection matrix u sed to select th e rows of the state vector co rrespond ing to the d esired measurements. G. Th eor etical E stimate of P erformance In this section, a method to com pute a theoretical estimate of the perfor mance and th e improvemen t achieved by the propo sed P ASE m e thod is d ev eloped. It is based on [ 10], where the authors proposed a tec hnique for estimatin g a priori the perfo rmances of the WLS estimator . Their work is extended in th is paper to fit the EnKF and compu te the relative gain between the two. The deriv ation is perfo rmed under the following assumptions, also made in [ 10]. The state vecto r is repr esented by w = [ v 1 , . . . , v | I | ] T . The fore c ast v ariance, the for ecast error time- correlation and load evolution m odel variance are assumed to be constant and id entical f or a cti ve and reactive powers. They ar e den oted resp e cti vely ( σ f i ) 2 , ψ f i and ( σ d i ) 2 . At each bus i , the apparent power magn itude | S f i | is used to re present the lo ad f orecast. In the analy sis fr amew ork, the shape of the load evolution m odel is not ne e d to be known, the value of the v ariance is suf ficient. The theor etical co mputation s are performed by using a linear Kalman filter . The covariance matrices are mad e time- in variant in order to obtain a steady-state for mulation of the filter [21]. From th is formulatio n , the covariance matrix of the system state can be compu ted and used to appr oximate the performan c e of the non-linear E nKF . The perform ance o f WLS can also b e com p uted since it can be seen as a Kalman filter th at is reset for each new estimation. T o ev aluate the perfor mance of the two state estimators over a p eriod of time T , the average ro ot mean square error of the voltage estimate (ARMSEV) is used as metric: ARMSEV = v u u t 1 T | I | T X t =0 | I | X i =1 E [ | ˆ v i [ t ] − v i [ t ] | 2 ] (15) where v i is the true voltage a t bus i and ˆ v i the estimated one. A lin ear versio n of th e power-flo w e q uations is used ; it is the fir st iteration of backward-forward sweep . A vectorized formu latio n is obtained by using a distrib ution load flow ( D L F) matrix, denoted by M , as describ ed in [22]. Th e relationship between th e injected power at each bus (represen ted by the vector s = [ s 1 , . . . , s | I | ] T , with s i the in jected power at bus i ) and the state vector is gi ven as: w = [ V 0 , . . . , V 0 ] + 1 V 0 M × s (16) where s is the conjug ate o f s . Se veral m atrices used by the Kalman equation s are defined. The load ev olution noise covariance matrix Q expr essed in terms of th e ap p arent power , and the forecast error cov ariance matrix R S are computed as follows: Q = diag (( σ d 1 ) 2 , . . . , ( σ d | I | ) 2 ) (17) R S = d iag (( σ f 1 ) 2 , . . . , ( σ f | I | ) 2 ) (18) The PMU measurem e nt er ror covariance m atrix is a p proxi- mated by assuming that th e variance of the voltage error whe n projected on to the real and imaginary axes is the sam e and equal to σ 2 P M U V 2 0 , where σ 2 P M U is the relative variance of the PMU measureme n ts such that R P M U = 2 σ 2 P M U V 2 0 × I | S | , where I | S | is the | S | × | S | identity matrix. The steady state covariance matrix of the state vector is computed by iterating th e Kalman eq uations. The covariance matrix is den oted by Σ ( · ) a . Th e iteratio n numb er is in dicated in the p arenthesis ( · ) . Such ma trix will converge to a steady state covariance matr ix Σ ( ss ) a . For each iteration , two other matrices are used to track the c ovariance matrix during intermediar y steps: Σ ( · ) p and Σ ( · ) u . They rep resent r espectiv ely the cov ariance matrix of the prio r state a n d the state after assimilation of PMU measuremen ts. At iteration 0 , the prior cov ariance matrix of the state is co mputed suc h th at: Σ (0) p = M × R S × M H (19) where ( · ) H indicates the Hermitian transpose (tran spose co n- jugate operato r). The up dated covariance matrix obtained after the assimilation of th e PMU m easurements is then c omputed : Σ (0) a = Σ (0) p − K H Σ (0) p (20) K = Σ (0) p H T ( H Σ (0) p H T + R ) − 1 (21) where H is th e observation matr ix for PMU measur ements. It is a selection matrix that relates state v ariables to the mea- surement vecto r . One can estimate th e ARMSEV p e rforman ce of WLS b ased on Σ (0) a : ARMSEV WLS = q 1 | I | trace (Σ (0) a ) . Any iteration it ( i t 6 = 0 ) is perfo rmed in 3 steps: first the prior covariance matrix Σ ( it ) p is computed based on the previous iteration, then the covariance matrix is upd ated using the PMU measurem ent c ov ariance matrix, and finally the pseudo- m easuremen ts are assimilated. Th e first two steps are such that (where H is the same matrix as in (20)): Σ ( it ) p = Σ ( it − 1) a + M × Q × M H (22) Σ ( it ) a = Σ ( it ) p − K H Σ ( it ) p (23) K = Σ ( it ) p H T ( H Σ ( it ) p H T + R P M U ) − 1 (24) The third step differs b ecause of the pseudo -measurem ent error time cor relation (see Section IV -E). The same time- differentiation method is used. The same updated matric e s a r e computed acco rding to ( 7)- (8) with only a few differences. Now H is th e in verse DLF matrix, m apping the state vector 6 to th e injected power ( H = M − 1 ). The fore cast error time correlation matrix is such that Ψ = diag ( ψ f 1 , . . . , ψ f | I | ) . Fi- nally , the matr ix R used in (8) is such th a t R = R S . The update equations thus b ecome: Σ ( it ) a = Σ ( it ) u − (Σ ( it ) u ( H ∗ ) T + C ) × K T (25) K = [Σ ( it ) u ∗ ( H ∗ ) T + C ] × [ H ∗ Σ ( it ) u ( H ∗ ) T + R ∗ + H ∗ C + C T ( H ∗ ) T ] − 1 (26) Once the steady state is reach ed after a f ew itera tio ns, the theoretical perf ormance of the EnKF can be comp uted. Th e ARMSEV error is su ch that: ARMSEV EnKF = q 1 | I | trace (Σ ( ss ) a ) . The relativ e ga in is ex- pressed as: G a in = ARMSEV WLS − ARMSEV EnKF ARMSEV WLS . V . V A L I DAT I O N A N D R E S U LT S The impr ovement in perf ormance achieved by the prop o sed P ASE meth od over WLS is evaluated by c o nsidering a 33- bus test distrib ution feeder [ 23] under nor mal o p erations. The WLS estimation prob lem is modeled in GAMS environment and solved using the M INOS solver . Attention has been paid to avoid poten tial numerica l issues. The ensemble size is set to L = 500 an d the p ower flow solutions obta in ed from h ( · ) are comp uted using the b ackward/forward sweep method [2 2]. The system is simulated over a p eriod of 24 h o urs. For the theoretical estimatio n, 5 0 iter a tions ( ( ss ) = 5 0 ) a re en o ugh to compute the steady-state of the state cov ariance matrix. A. Lo ad Evolution Model The (bus) load ev olution model was developed using a fine- grained energy co nsumption d ataset from Ontario, Canada. The dataset used to build the model is describ ed in [2 4] and comprises instantaneous acti ve power co nsumption data from 20 hom es, collected ov er e ig ht month s, with a resolution o f 6 seconds. Th e dataset is split random ly into two sub sets, on e for der iving the ch a r acterization ( training set), an d one f or the validation process ( testing set). No d istinction is made b etween the size of the houses nor the time of the year . T he resulting dataset is a collection of a few thousan ds of traces. Alth ough 20 homes may seem to be a limited sample size, considering the daily power trac es indepen dently allows to hav e a large number of unique profiles. Moreover th e 20 households cover a wide range o f living are a sizes and energy co n sumption patterns which increases the trac e diversity . Let n be th e n u mber of househ olds con nected to a bus. Using th e train ing set, empirical distributions fo r load chan ges were constru c ted for different values o f tim e-steps ∆ T an d aggregation levels n . A Lap lace d istribution described by a scale p arameter b (and variance σ 2 = 2 b 2 ) was found to be a goo d fit. The mean value is set to zer o sin c e as many positive and negative load chan ges ar e expected. This implies that th e tran sition mod e l is the identity , while its u ncertainty is char a cterized by the Laplace distribution. The influence of ∆ T and n o n the distribution variance is illu strated in Fig. 1. The variance e ssentially describes the load v ariation over time, a small value imp lying little variations. It is n o ted that as n increases and ∆ T shrinks, the v alue of σ 2 diminishes. Load evolution model variance: 2b 2 (kW) 2 10 20 30 50 100 200 Aggregation level (number of houses) 1200 600 300 96 48 24 12 6 Time-step (s) 20 40 60 80 100 120 Variance (kW) 2 Fig. 1. Influence of aggre gation le ve l and time step on the scale parameter b . The fit of the Lapla ce distributio n is visually good for all the time-steps and aggregat ion lev els considered here. It was assumed th at load chan ges are unco rrelated be twe en buses; which c a n be verified to hold true from the dataset, fo r any value of n and ∆ T up to 30 minu tes. The values o f σ 2 are d eriv ed emp irically as a func tio n of n and ∆ T . They are used to co mpute the ev olution step of the EnKF . Since no reactiv e p ower consumption da ta set was a vail- able, a similar model is assumed for reactive power c h anges. Howe ver , activ e a nd re activ e power con sumption chan g es are assumed to be ind ependen t, which is a common assumption in DSSE liter ature. Th e pr oposed method is generic and can be applied to any dataset from across th e globe. B. T est Distrib ution System The 33- bus test feeder data include s a cti ve and reactive power loads at each bus; bus-1 is the substation tran sformer bus, with V 0 set to 12.66 k V . The number of houses n i aggre- gated at a b us i is selected such that n i = n 11 P 33 bus i /P 33 bus 11 where n 11 = 1 0 houses and P 33 bus i is the static 33-bus active power lo ad at bus i . The co rrespond ing distribution trans- former traces are generated f rom the secon d half o f the dataset, by summing the desired nu mber of profiles, picked rand omly . Each trace is then scaled so th at th e mean of the pr ofile matches the load values. The values given b y th e empiric a l function in Section V - A are scaled accordin gly . Because no dataset for reactive po wer consump tion is av ailable, activ e and reactive po wer pro files are g enerated indepen dently from the same dataset. C. Measurement Mod el The simulation mo dels used fo r measurements ar e d escribed in this section . PMU: the PMU measurement er ror is simulated as an additive white Gau ssian noise o f n ominal variance σ 2 P M U , fo r both voltage magn itudes and ang les. The r eadings ˜ V s and ˜ δ s provided by the PMU at each bus s ( s ∈ S ⊆ I ) have an error variance suc h that E [˜ a 2 ] = σ 2 P M U ∗ ˜ a 2 , where ˜ a ind icates either the voltage m a gnitude or angle. The measurement errors are ind ependen t across buses, and the voltage magnitud e error indepen dent of the angle error . The PMU resolu tion is set to 1% ( σ P M U = 0 . 01 ); the PMU p lacement m ap S is determin ed u sing a g reedy method [10], i.e., PMUs are seq uentially added at th e loca tion th at provides th e mo st 7 improvement (with 32 load buses, a maximum of 3 2 PMUs). The placement o f PMU s is b eyond the scope of this work; many r esearchers have a ddressed this issue, see f or example [25]. Pseudo-measurements : the forecasts P f i and Q f i are taken as the mean value of the load p rofile generated at each distribution transfor mer i , as in [1 0]. Th ey are constan t over the simu la te d period . Using the training set, the no m inal standard deviation of th e fore c ast was ev aluated and set to σ 0 = 30% , for both a c tive an d r eactiv e powers, irre sp ectiv e of the aggr egation level. Therefor e f or each bus i , σ f p i = σ 0 P f i and σ f q i = σ 0 Q f i (6). The constant app arent power f orecast | S f i | is such that | S f i | = | P f i + j Q f i | . Fina lly ea ch σ f i (18) is compu ted as σ f i = σ 0 | S f i | . Pseudo -measureme n ts w ith a Gaussian distribution are used as “b est-guess” initial en semble (Section IV -C). Error time-correlatio n : ψ p i and ψ q i are ev aluated as fol- lows: since the sam e d ata is used for g enerating the active and re a cti ve p ower pr ofiles, ψ p i and ψ q i are equ al. They are ev aluated on the tr aining set. Given a n ag gregation le vel an d a tim e-step length, load profiles are built. The autocorrelation function R e i of the d ifference between the p rofile and its mean (represen tin g the forecast er ror) is co mputed. The value of ψ p i and ψ q i is gi ven by R e i (∆ T ) . D. V alid ation The theore tica l an d simulation results are presented in Figs. 2a-2c, ob tained b y a veraging the results of se veral real- izations. A r ealization is defined as th e observed per f ormance of both th e WLS an d P ASE on the 33- bus system. For ea c h realization, new lo ad profiles are gener ated based on the testing set, while the oth er p arameters stay the same. The perfor mance of the WLS and P ASE are plo tted alongside with the theor etical ones in Fig. 2a, where a time-step of 6 seconds has been used . WLS has been stud ied in [1 0] using synthetic da ta. Similar trend s are observed here with real data. Note that since WLS is snap shot-based , the size of the tim e- step d oes n ot matter . For P ASE, the theoretical results are close to the actual p erforman ce observed in simulation a s the number of PMUs introd uced in the system increases, wh ich validates the theore tical app roach. Similar trend s a r e ob served for different time-steps. T he actual g ain br ought abou t by P ASE is compa r ed with the theoretical one in Fig. 2b for a time step of 6 secon ds. Finally the influence of the time-step on th e gain is comp ared in Fig. 2c fo r two PMU con figuration s ( 5 PMUs and 20 PMUs). Theory and simu lation follow the same trend. The gap between th e ory an d simulation is relativ e to that observed in Fig. 2b. For 5 PMUs, a separation between the curves is observed. E. Co mparison Between WLS a nd Pr oposed P A SE Metho d The results pr e sented in Fig. 2a illustrate th e improvements achieved by the pro posed P ASE m ethod. Clearly , usin g a lo ad ev olution m odel improves the perf ormance of the estimator; giv en an arbitra ry target err or of 0.00 4 p. u., WLS requ ires more than 1 0 PMUs wh ile P ASE on ly 4. Even wh en each bus of th e distribution system is mon itored b y a PMU, the propo sed P ASE m ethod still bring s ab out an improvement of more than 40% wh en using a time- step of 6 second s. As illustrated in Fig. 2 c, hig her gains are ob tained for smaller time-steps. Ind e ed, for larger time-steps, the load ha s more chances o f changin g by a large magnitu d e between two estimates and thus has less inertia. E ven fo r large time- step (e.g., 10 min s) there is a gain of abo ut 15%. In pr actice, the granularity of the time-step d epends on the a vailable computatio nal speed . The smallest time-step considered in this work is 6 seconds an d represen ts a lower -bound on what was tried o u t. In compa rison, the DSSE p roblem was solved at each step in un der 1 secon d. F . En gineering I nsights In practice, the LDC will need to make trade-o ffs in th e choice o f the f ollowing p arameters: n umber of PMUs, their accuracy an d the time scale. The influence of PMU accuracy on the theore tical gain ach ieved by P ASE is shown in Fig. 3a, the three par ameters considered are depicte d in the plo t. The maximum gain is attained for a PMU error v ariance of abou t 1%. Clearly as the PMU measure m ent standard deviation decreases (i.e . , the PMU becom es m ore an d more acc u rate) the g ain ach iev ed by P ASE decreases since the load ev olution model is no t as usefu l in such c ircumstances. Similarly , as the standard deviation of the PMU increa ses, the gain de c r eases, since the load ev olution model h as to co mpensate fo r both poor f orecast accur acy and po or PMU measur ement accuracy . This fig u re a lso illustrates the r ole o f th e time- step, th e gain achieved by the filtering techn ique de c reasing as the time-step increases, underlining the limits of the load e volution mod e l. The trade- o ff between the three pa rameters conside red is illustrated by Fig. 3b: two PMU ac c uracies are used to draw the plots. A n arbitra r y target error is fixed and the m inimum number of PMUs required is determ ined as a functio n of the time-step. Clearly , the time-step has little in fluence on a very accurate PMU. Ho wev er , th e mor e accurate the PMU, the more costly it will be. With th e same numb er o f PMUs placed in the sy stem (4 ), choo sing a PMU ten time s less accurate will provide the same per f ormance given that a time-step small enoug h (6 seconds) is chosen. V I . C O N C L U S I O N S A novel P ASE m e thod for DSSE and its analysis framework were presented . The P ASE metho d perf o rms the fusion of measuremen ts and pseudo- measuremen ts and requ ires fewer PMUs tha n WLS to ach ieve the sam e estimatio n err or , for time-steps u nder 1 5 minutes. Engin e ering insig h ts were pre- sented highlig hting the major trade-offs in th e ch oice of decision variables f or the LDC. Using a smaller time-step allows the LDC to r elax the require m ents on th e PMU qu ality and their nu mber . Th ere are sev eral re maining ch allenges, such as the influence o f distributed gen eration and its mod eling as well as the impact of an unbalance d system on P ASE. R E F E R E N C E S [1] S . Paudyal, C. A. Canizar es, and K. Bhatta charya, “Optimal Operation of Distrib ution Feeder s in Smart Grids, ” IEEE T ransactions on Industrial Electr onics , 2011. [2] O . Ardakani an, S. K esha v , and C. Rosenber g, “Re al-T ime Distrib uted Control for Smart Electric V ehicle Char gers: From a Static to a Dynamic Study, ” IEEE T ra ns. Smart Grid , 2014. 8 0 5 10 15 20 25 30 35 Number of PMUs placed 0 0.002 0.004 0.006 0.008 0.01 0.012 ARMSEV (p.u.) Time-step: 6s WLS theoretical PASE theoretica l WLS simulation PASE simulation (a) 0 5 10 15 20 25 30 35 Number of PMUs placed 5 10 15 20 25 30 35 40 45 50 55 Performance (ARMSEV) gain (%) Time-step: 6s Theoretica l Simulation (b) 10 0 10 1 10 2 10 3 Time-step value (s) 0 5 10 15 20 25 30 35 40 45 50 55 Mean gain (%) 5 and 20 PMUs Theoretical Simulation (5 PMUs) Simulation (20 PMUs) (c) Fig. 2. (a) ARMSEV va lue function of the number of PMUs. The performanc e of the proposed P ASE method is compared with W LS. The theoretica l results are also compared against simulation results. A lower valu e means better performance. (b) Comparison of the gain from using P ASE ov er WLS on ARMSEV depending on the number of PMUs. The theoret ical results are compared to the observed gain in simulation. (c) Influence of the time-step on the mean performance gain. The theoretic al results are compared to the observed gain in s imulati on. The time-step axis has a logarit hmic scale. T he error bars represent the vari ance. 10 -2 10 -1 10 0 10 1 PMU standard deviation (%) 5 10 15 20 25 30 35 40 45 50 55 Gain (%) 6s, 48s and 10 mins time-step 5 PMUs, 6s 20 PMUs, 6s 5 PMUs, 48s 20 PMUs, 48s 5 PMUs, 10min 20 PMUs, 10min 0 100 200 300 400 500 600 Time-step considered (s) 0 5 10 15 20 25 30 Number of PMUs required ARMSEV target = 2.25e-03 (p.u.) 0.1% accuracy 1.0% accuracy Fig. 3. (a) Influence of PMU qual ity on the (theore tical ) gain achie ved by P ASE ov er WLS. The PMU accurac y is characteriz ed by its measurement error standard devi ation, a lower valu e means a m ore accurate PMU. A logarithmi c scale is used for the v ariance axis. (b) Minimum number of PMUs required to achie ve an avera ge tar get error of 2 . 25 e − 3 p.u., function of time-step and for two PMU accuracies (computed using the theoret ical formulation). [3] D . Atanacko vic and V . Dabic, “Deployment of real-time state estimator and load fl ow in BC Hydro DMS - chall enges and opportunit ies, ” in 2013 IEE E P ower Energy Societ y General Meetin g , 2013. [4] A . Mont icelli , State Estimation in Electric P ower Systems . Springer , 1999. [5] A . von Meier , D. Culler , A. McEachern, and R. Argha ndeh, “Micro - synchropha sors for distribut ion systems, ” in ISGT, IE EE PES , 2014. [6] R. N. Rodrigues, J. K. Zatta, P . C. C. Vi eira, and L. C. M. Schlichting, “ A lo w cost prototype of a Phasor Measurement Unit using Digital S ignal Processor, ” in IEEE B iennial Congre ss of Arg entina , Jun. 2016. [7] C. A. Fanti n, M. R. C. Castill o, B. E. B. d. Carv alho, and J. B. A. London, “Using pseudo and virtua l measure ments in distributi on system state estimation, ” in IEEE PE S T&D-LA , 2014. [8] A . Primadianto and C. N. Lu, “ A Revi e w on Distribut ion System State Estimation, ” IEE E T ran sactions on P ower Systems , 2016. [9] A . K. Ghosh, D. L . Lubkeman, M. J. Downey , and R. H. Jones, “Distrib ution circuit state estimat ion using a probabilisti c approach, ” IEEE T ransact ions on P ower Systems , 1997. [10] L . Schenato, G. Barchi, D. Macii, R. Arghand eh, K. Poolla, and A. V . Meier , “Bayesian linear s tate estimation using smart meters and PMUs measurement s in distribut ion grids, ” in SmartGridComm , 2014. [11] S. Alam, B. Nataraj an, and A. Pahw a, “Distributio n Grid State Estima- tion from Compressed Measurements, ” IEE E T rans. Smart Grid , 2014. [12] H. W ang and N. N. Schulz, “ A re vised branch curre nt-based distributio n system state estimation algorithm and meter placement impac t, ” IEEE T ransactions on P ower Systems , vol. 19, no. 1, pp. 207–213, Feb . 2004 . [13] C. Klauber and H. Zhu, “Distri buti on system state estimation using semidefinit e programming, ” in NAPS , 2015. [14] M. B. D. C. Filho and J. C. S. d. S ouza, “Forec asting-Aid ed State Estimation - Part I: Panorama, ” IEEE T ra ns. on P ower Systems , 2009. [15] S. -C. Huang, C.-N. Lu, and Y .-L . Lo, “Ev al uation of AMI and SCAD A Data Synergy for Distribut ion Feeder Modeling, ” IEEE T ran s. Smart Grid , 2015. [16] S. Sarri, M. Paolone, R. Cherkaoui , A. Borghetti , F . Napolitano, and C. A. Nucci, “State estimati on of acti v e distrib ution netwo rks: com- parison between WLS and iterated Kalman-filter algorithm in tegra ting PMUs, ” in Innovative Smart Grid T echn ologi es (ISGT E ur ope) , 2012. [17] G. Evensen, “The Ensemble Kalman Filter: theoret ical formulation and practi cal implementation, ” Ocean Dynamics , vol. 53, Nov . 2003. [18] R. Singh, B. C. Pal, and R. A. Jabr , “Choice of estimator for distribution system state estimation, ” IET Gener . T ransm. Distrib . , 2009. [19] K. W ang, Y . Li, and C. Rizos, “Prac tical Approaches to Kalman Fil- tering with Time-Corre lated Measurement Errors, ” IEEE T rans. Aer osp. Electr on. Syst. , 2012. [20] M. G. Petov ello , K. OKeefe, G. Lachapelle , and M. E. Cannon, “Conside ration of time-correl ated errors in a Kal man filter applic able to GNSS, ” Jo urnal of Geodesy , 2009. [21] B. Anderson and J. Moore, Optimal Fi ltering . Prentice-Hal l, 1979. [22] J. -H. T eng, “ A direct approach for distributi on system load flo w solu- tions, ” IEEE T ra nsactions on P ower Deliv ery , J ul. 2003. [23] M. E. Baran and F . F . W u, “Network reconfigurati on in distributi on systems for loss reduction and load balancin g, ” IEEE T rans. on P ower Del. , 1989. [24] O. Arda kanian, S. Kesha v , and C. Rosenberg, “Marko vian Models for Home E lectri city Consumption, ” in ACM SIGCOMM , 2011. [25] R. Singh, B. C. Pal, and R. B. V inte r , “Measurement Place ment in Distrib ution System State Estimation, ” IEEE T rans. on P ower Syst. , 2009.
Original Paper
Loading high-quality paper...
Comments & Academic Discussion
Loading comments...
Leave a Comment