A stochastic approach to estimate distribution grid state with confidence regions
Rasmus L. Olsen Sina Hassani Troels Pedersen Jakob Gulddahl Rasmussen Hans-Peter Schwefel

	



1]Rasmus L. Olsencor1fn1
rlo@es.aau.dk
1]Sina Hassanifn2
sinah@es.aau.dk
1]Troels Pedersenfn3
troels@es.aau.dk
2]Jakob Gulddahl Rasmussenfn4
jgr@math.aau.dk
1,3]Hans-Peter Schwefelfn5
schwefel@griddata.eu

[cor1]Corresponding author
[1]organization=Department of Electronic systems, Aalborg University, addressline=Fredrik Bajersvej 7A, city=Aalborg, postcode=9220, country=Denmark
[2]organization=Department of Mathematical Sciences, Aalborg University, addressline=Skjernvej 4A, city=Aalborg, postcode=9220, country=Denmark
[3]organization=GridData GmbH, city=Anger,  country=Germany







 

Widely available measurement equipment in electrical distribution grids, such as power-quality measurement devices, substation meters, or customer smart meters do not provide phasor measurements due to the lack of high resolution time synchronisation. Instead such measurement devices allow to obtain magnitudes of voltages and currents and the local phase angle between those. In addition, these measurements are subject to measurement errors of up to few percent of the measurand. In order to utilize such measurements for grid monitoring, this paper presents and assesses a stochastic grid calculation approach that allows to derive confidence regions for the resulting current and voltage phasors. Two different metering models are introduced: a PMU model, which is used to validate theoretical properties of the estimator, and an Electric Meter model for which a Gaussian approximation is introduced. The estimator results are compared for the two meter models and case study results for a real Danish distribution grid are presented.




  Grid, state estimation, modelling







§ INTRODUCTION

Due to the penetration of renewable energy sources into the distribution grid and due to additional loads and production of energy, e.g. from electrification of transport and building heatings, the reliability and efficiency of electricity distribution grids is challenged, and planning and operations of these grids can no longer be performed based on assumptions on loads and generation units. Therefore, grid management, grid optimization, and grid planning using knowledge of the current and historic operational state of the grid is increasingly important <cit.>. 
Cable loads and voltage levels that exceed specific limits, must be detected and addressed by proper means, starting by enabling monitoring of the complete distribution grid. 
Furthermore in order to optimize voltage quality and grid efficiency, communication-based hierarchical control systems have been proposed and shown to be efficient <cit.>. The aforementioned functions of the distribution grid show best performance when the electrical state of the grid is known in all grid locations - while it is infeasible and costly to measure all these grid locations.  The necessary data which shows the condition of the grid and is used for control purposes are called the state variables of the system or in short, states of the system.

The large scale of the distribution grid and the presence of communication failure restrict the presence of necessary data for the algorithms. In addition, the impact of the uncertainty of the measured data needs to be considered. Therefore, a process in order to estimate the necessary grid state from the measured data and quantify the impact of measurement errors and missing measurements is necessary. This process is called state estimation in the literature <cit.>. The measured or forecast data used for the estimation process are called the inputs of the system.The state estimation process has for decades been widely used for the high voltage transmission system. These state estimation methods are now starting to become relevant in distribution system state estimation (DSSE), i.e. at medium and low-voltage level. While there are many similarities, the are also some major challenges <cit.> when applying DSSE. The key difference is within the observation possibilities, i.e. type and accuracy of measurements, where in the high voltage level, so-called PMUs are able to provide highly accurate synchronized phasor measurements for a large number of measurement location. In contrast to that, measurement devices in the distribution grid entail power quality measurement devices on primary substations, electrical measurement devices on secondary substations, smart meters and inverters at customer connections; none of the latter use accurate clock syncronisation to millisecond level and therefore they only provide magnitudes and local phase angles of RMS voltages and currents. Furthermore, the measurement errors caused by the device type and frequently also from the use of current transformers leads to a higher degree of uncertainty in the measurements; which must be accounted for in the DSSE applied in the MV or LV grid. This uncertainty has strong implications on applications and use of data driven decisions, e.g. planning of new cables and loss reduction <cit.>, since not knowing the accuracy of the data may easily lead to wrong decisions taken by the DSO, e.g. replacing the wrong cable or not detecting inefficient grids or faulty grids. Adding confidence regions to data will aid decision makers to know how much trust to put into application results and act accordingly.  

High unobservability of the distribution system is the second challenge of the DSSE, i.e., often the amount of measured data is not enough to estimate the states of the system <cit.>. This problem is due to the huge number of connection points of the distribution system compared to the transmission system. The ratio of imaginary and real parts of the line impedance (x/r) is low in the distribution grid, making the linear state estimation techniques (DC state estimation) impractical <cit.>. 


Finally, the use of heterogeneous communication networks such as celullar, power-line communication, and low-bandwith technologies such as LoRa causes bandwidth limitations and latency and losses in collection of measurement data. 

According to the literature most of the DSSE methods rely on the weighted least square (WLS) approach <cit.>. The WLS method is used in different situations which have different types of states and measurements. It should be stated, that besides the measurement data, the algorithms use forecast data to avoid unobservability and achieve better performance. Other algorithms are used for the DSSE such as Kalman filter, interval state estimation, gradient-based method, and semi-definite programming <cit.>. 

A wide variety of the states of the DSSE algorithm are used in the literature. Many papers consider the voltage phasors of the nodes as the state variables of the estimation technique. These states are used in both polar and rectangular forms. The polar form of the voltage phasors is used in <cit.>, while the rectangular form of the voltage phasors is used in <cit.>. Another set of state variables in the DSSE methods are current phasors of the branches. These methods are the most popular methods in the literature <cit.>. The current phasors are used in polar <cit.> and rectangular forms <cit.>. It is worth mentioning that in <cit.>, the slack bus voltage magnitude is included in the states, and in <cit.>, the slack bus voltage phasor in rectangular form is included in the states.

The inputs of the DSSE algorithms vary widely according to the method and the availability of the data. According to the literature the data are different combinations of active and reactive power of the loads and generators in the distribution grid, node voltage magnitudes, branch current magnitudes, and load current magnitudes <cit.>. The availability of the phasor measurement units (PMUs) makes it in principle possible to have the phasor measurements, while their deployment in distribution grids in practice is non-existent due to the tight requirements to time synchronization. Nevertheless, the phasor measurements are considered as the inputs in some DSSE algorithms <cit.>. 


On the other hand, the increasing deployment of Digital Electrical Measurement devices in transformer stations, and of Smart Meters and Smart Inverters at customer connections provides an increasing set of measurement locations in distribution grids, while these devices cannot determine absolute phase angles, hence are limited to determine magnitudes of voltages and currents and relative phase angles between the voltage and the current at the measurement point. The general poor clock synchronization of smart meters has previously made SE and DSSE at the low voltage level challenging. The impact of clock offsets on loss calculation applications has been analyzed in <cit.>: the authors show that the choice of the specific algorithm for loss calculation has a drastic impact on the sensitivity of the calculated grid losses due to clock deviation errors. The paper furthermore addresses errors in the loss calculation caused by duration of measurement interval and by measurement errors; however, the considered loss calculation algorithms assume that all required input variables have been measured, so no DSSE approach is used in that paper.


In this paper, a DSSE algorithm is proposed based on the available measurements of the distribution systems. The DSSE algorithm is presented to fulfill the main objectives of the state estimation process, i.e. the estimation and calculation of confidence regions of the unknown parameters of the system. The DSSE algorithm estimates the node voltage phasors and branch current phasors of the distribution grid. The inputs of the system are the available measured node voltages and the branch currents at a subset of grid locations and these measurements are subject to measurement errors. The setup of the framework enables us to assess different metering models, ranging from PMUs to measurement devices with low synchronisation precision. 






Our contribution in this paper is a methodology to assess confidence regions for model based estimates of voltages and currents in electrical grids. We present a stochastic model framework for carrying out the estimation and derive confidence regions analytically, and subsequently show the use in a real world scenario from a Danish DSO. 

Section <ref> describes the system model and the model of the measurement devices. Section <ref> derives the maximum likelihood estimator and the confidence ellipses for the voltage and current phasors. Section <ref> introduces the assessment approach and Section <ref> applies that approach to a real-life scenario of a distribution grid based on actual measurements.
Finally, Section <ref> summarizes the results and presents an outlook to future work.


















 











 

 

 




§ STOCHASTIC MODELS OF ERRORS IN REALISTIC ELECTRICAL MEASUREMENT SCENARIOS






 §.§ System Model
 







Consider a graph  with nodes and edges representing possible measurement points and power lines, respectively. Here, either an MV or an LV grid is considered.  In the MV  case,  nodes represent all MV busbars and sleeves including the  MV side busbar of the HV:MV transformer. For LV grids, nodes represent this first busbar is the LV side of the secondary transformer, the junction boxes (JBs), sleeves and customer connection boxes (CCBs).








 



The grid state  is given by a  voltage phasor for each node and a current phasor for each edge of the graph. This amounts to a number n of voltage phasors x_1…,x_n and a number m of current phasors  x_n+1… x_n+m.  The kth phasor is  of the form








    x_k = 
        
         u_k e^jθ_k,    k=1,…,n,(voltage)
     
        i_k e^j(θ_k+ϕ_k),    k=n+1,…,n+m,(current),

where  u_k, i_k, θ_k and ϕ_k denote respectively the "true" values of voltage magnitude, current magnitude, voltage phase and local phase angle between voltage and current. Collectively, these N=n+m phasors constitute the state vector x∈ℂ^N.[Including both current and voltage phasors of the grid as the state vector facilitates inclusion of noise for all measurements. Furthermore, it enables a generic representation of the equation systems for a variable number of observed voltages or currents. 
]  
For simplicity, we focus on a single-phase representation of the grid.
We remark that generalizations to three phase unbalanced grids are possible; the number of phasors in the state vector will triple in that case. Note however, that another type of measurement error would need to be addressed for the practical application to three phase-grids, namely the wrong phase assignment at measurement devices. The latter is out of scope for this paper.


The state vector fulfills a number of linear constraints given by the grid. Kirchhoff's equations constrain the current phasors and further linear equations link the voltage phasors at two ends of a line with the current phasor on that line and the line impedance. These linear constraints can be summarized in a single matrix C by the equation 

    Cx = 0,

see <cit.>[In busses with injected power, the injected current is summed up with the load current; therefore, Kirchoff's equation holds in this situation. In addition, we aim to estimate all the currents and voltages hence, zero-injected busses are not Kron-reduced in this work.]. For generalizations of the contributing equations to 3-phase grids while remaining  in the space covered with linear equations, see <cit.> for instance. 

In most situations it is impossible or impractical to measure the whole state vector directly. In fact,  only noisy measurements for functions of a subset of x can be achieved.   For MV grids, those typically include measurements at the primary and secondary transformers and on a subset of MV lines. For LV grids measurements  typically include the secondary substation, a subset of the CCBs, or possibly in special cases also intermediate JBs. In this paper, we consider measurement devices that collect RMS measurements on voltage magnitude, current magnitude and local phase angle between voltage and current. We refer to these measurement devices as meters for brevity.  Measurement data is collected by wired or wireless communication. Furthermore, the data is  time-stamped and sampled in a synchronized manner so we focus on one particular time-interval across the grid, where clock synchronization errors are neglected here (see <cit.> for impact of inaccurate clocks). 

We assume that measurements can be obtained as in total K real-valued entities from meters that are placed in the grid. 
The concatenation of all measurements from all meters can then be represented by an K-dimensional real vector y = (y_1,…,y_K). The measurement vector y is modeled as a stochastic vector depending on a subset of the state variables, expressed as  

    y = f(Dx,ϵ) ∈ℝ^K,

with ϵ being a noise vector and the measurement matrix D ∈{0,1}^K× N is defined so that the complex vector Dx contains only the phasors for which measurements are obtained. 














 










The estimation problem at hand is to recover the full state vector x  from vector  y, which contains the accumulated measurements. At its disposal, the estimator has knowledge of the grid topology specified by the matrix C (the linear grid equations mapped from grid structure and grid impedance), the measurement matrix D along with information on the metering function f and  distribution of measurement errors.





 §.§ Meter Models: PMU and Electrical Meter
 
In order to exemplify the measurement model description from the previous section and in preparation for the use-case scenario later  in the paper, we now introduce two different meter models for comparison: (1) A PMU measures the two phasors of voltage and current at on specific end of a line; therefore it determines four real-valued measurands, namely the cartesian coordinates of the voltage and current phasor at this measurement location. As an error model, we assume independent complex 2-dimensional Gaussian noise to the Cartesian phasor representations of voltage and current respectively. This model is relevant to validate the theorems presented in the next section, and is used for validation of the confidence regions in the use-case section.


(2) An Electrical Meter (EM) measures voltage and current magnitudes and furthermore the local phase angle between voltage and current phasor. Consequently it determines three measurands per measurement location. As measurement error, we here assume independent Gaussian noise on the magnitudes of voltage and currents and on the angle, ϕ, between voltages and currents, reflecting real world low-cost meter measurement methodologies. 


As the metering models provide different measurands as output (collected in the real-valued vector y), a data preparation stage is used to translate the vector y into a unified input z to the estimator, see Figure <ref>. This data preparation stage is specific for each meter model and introduced in the next subsection. 




 §.§ Data Preparation for Estimation

The target is to derive (approximate) phasor representations for the current and voltage phasors as approximations for the true phasors in the complex vector Dx ∈ℂ^N from the measurements in the vector y∈ℝ^M.  To this end,  we define a vector function g = (g_1,…,g_K) where each entry gives the complex-valued form of a voltage or current phasor derived from the measurements

For the PMU measurement model, the function g is trivial to set up as the meter model already delivers the rectangular coordinates of the desired phasors. Therefore, assembling the complex phasors from the two coordinates is sufficient. 

For the EM measuring model, the voltage and current phasors need to be approximated from the measured magnitudes and local phase angle. , i.e.

    z_k = g_k( y_k) =
    
        (u_k + ϵ_u,k) e^j(θ_k + ϵ_θ,k)   (voltage)
    
        (i_k + ϵ_i,k)e^j(θ_k + ϕ_k + ϵ_ϕ,k + ϵ_θ,k)   (current)

Here,  ϵ_u,k,ϵ_i,k,ϵ_θ,k and ϵ_ϕ,k denote the respective measurement errors on magnitudes or phase angles (which are already added in the metering model itself). Measurement errors at meter k are regarded as independent random variables with zero mean (unbiased errors). Furthermore, measurement errors from different meters are considered independent. 


Since the EM metering model does not determine the voltage angle, extra assumptions are needed to define z_k. In particular,  in this case we introduce a pseudo-measurement by assuming the voltage phase to be zero. This is a good approximation since the voltage phase θ_k is usually small in single-phase representations of distribution grids, when the reference point (e.g. the substation) is set to phase-angle zero.  Thereby, we obtain a pseudo measurement of the form (<ref>) with θ_k + ϵ_θ,k= 0. 

Finally, we obtain a vector z = (z_1,…,z_K) derived from measurements and influenced by the introduction of (pseudo-) measurements (here for the voltage phase angle).



 §.§ Complex Gaussian Approximation of measurement errors provided by EM model
 
The probability distribution of z is unknown in general.  For the EM model, even if the distribution of the different ϵ_.k are known, the distribution of z_k is unavailable due to the non-linear transform g_k in (<ref>).  
 For mathematical tractability, however, we approximate the distribution of z_k by a complex normal distribution defined such that it gives the same first and second moments as z_k, i.e 

    z_k ∼ CN(μ_k,Σ_1,k,Σ_2,k)

with mean μ_k, variance Σ_1,k = Cov(z_k,z_k) and pseudo variance Σ_2,k = Cov(z_k,z_k^*), with ^* being the complex conjugate.
The assumption of independent errors allows us to state the approximation for the distribution of  the vector z as 

    z ∼ CN(μ,Σ_1,Σ_2)

where the mean vector μ contains μ_1,…, μ_K and the covariance matrix Σ_1 is a  diagonal matrix with entries Σ_1, k, k = 1,…,K. The pseudo-covariance matrix Σ_2  is also diagonal and defined similarly. 





 

The mean, covariance and pseudo-covariance for z_k can be computed for a given distribution of measurement errors by use of Proposition <ref> given in appendix. We obtain  the  mean as
 
    μ_k = 𝔼[z_k] =
         
         u_k e^jθ_k c_ϵ_θ,k(1)   (voltage)
    
         i_k e^j(θ_k+ϕ_k)c_ϵ_θ,k(1)c_ϵ_ϕ,k(1)   (current),

 the variances as

    Σ_1,k
    
    
    
         =
    
    
          
         (1-|c_ϵ_θ,k(1)^2|)u_k^2 + σ_u,k^2    (voltage)
    
         (1-|c_ϵ_θ,k(1)^2c_ϵ_ϕ,k(1)^2|)i_k^2 + σ_i,k^2   (current)

and pseudo-variances as

    Σ_2,k
           =
     
     
     
    
    
        e^2jθ_k((u_k^2+σ_u_k^2)c_ϵ_θ,k(2) - u_k^2 c_ϵ_θ,k(1)^2)   (voltage)
    
        e^2j(θ_k+ϕ_k)((i_k^2+σ_i_k^2)c_ϵ_θ,k(2)c_ϵ_ϕ,k(2)
         - i_k^2 c_ϵ_θ,k(1)^2c_ϵ_ϕ,k(1)^2)    (current),

 
 where c_ϵ(t) = 𝔼[e^jtϵ] denotes the characteristic function of a random variable ϵ.

The quality of the complex normal distribution approximation depends on the true distribution of z. For example, in the EM model, we assume that the measurement error of the measured magnitudes  and of the measured angle are  normally distributed. In such case for large values of the variance of the angle, the true distribution will have a banana shape, and the resulting approximation will be poor, while for small values of the variance of the angle, the approximation is good.

In Figure <ref> we illustrate how the two models, a) the PMU measurement model and b) the EM measurement model, are implemented. In essence, we need to do more than just creating a measurement vector y, but also need preparation and calculation of the co-variance matrices that are required for the estimator, which we will describe in the subsequent section. For a real setting, we will receive the measurement vector y and will have to calculate the covariance matrices based on the measurements and based on properties of the measurement device, see Section <ref>. 

















































































































 
























  
  







§ MAXIMUM LIKELIHOOD ESTIMATION OF GRID STATE


























We now turn to deriving a maximum likelihood estimator for the entire grid state vector x based on the phasor vector z representing derived observations of a smaller number of directly observable state variables Dx. The methods for obtaining the maximum likelihood estimator and its properties are an extension of the methods used in <cit.>, but with a more general model and with additional results. In this section, the measurement model is given by

    z ∼ CN(Dx,Σ_1,Σ_2).

To do this, it is assumed that the linear constraint (<ref>) along with the noise covariance and pseudo-covariance matrices are available.

For generality,  we derive the estimator under slightly generalized assumptions. Specifically,  we generalize the constraint (<ref>) to

    Cx = c,

and arbitrary complex values in D, c, C, Σ_1 and Σ_2, provided that the matrix 

    P = Σ_1^*-Σ_2^HΣ_1^-1Σ_2

is positive semi-definite where superscript H denotes complex conjugate transpose. Furthermore, we denote W=Σ_2 (Σ_1^*)^-1 as a convenient short notation, since both W and P appears in various places in the below results. 

For brevity, we formulate the estimator using complex augmented vectors and matrices. Given a general complex vector v or matrices B_1 and B_2, their augmented versions are marked by an overbar and defined as

    v̅ = [   v; v^* ],     (B_1,B_2) = [   B_1   B_2; B_2^* B_1^* ].

In augmented notation the model (<ref>) is written as

    z̅∼ CN(D̅x̅, Σ̅),

where D̅ = (D,O), Σ̅= (Σ_1, Σ_2), and O denotes an appropriately sized matrix of zeros.

We estimate x by maximum likelihood estimation, or equivalently least squares estimation, by minimising

    (z̅ - D̅x̅)^H Σ̅^-1 (z̅ - D̅x̅),

subject to the constraints c=Cx. This is solved by minimising the Lagrange function

    ℒ(x,λ) = (z̅ - D̅x̅)^H Σ̅^-1 (z̅ - D̅x̅) + (2λ^H (Cx-c)).

This leads to the following result:

Assume (<ref>) and (<ref>). Then the maximum likelihood estimate is a solution to the equation 

    [ g̅; c̅ ]
    = 
    [   G̅ C̅^H;   C̅    O ][ x̅; λ̅ ],

where G̅ = (G_1,G_2), C̅ = (C,O),
G_1 = 2D^H(P^*)^-1D,
G_2 = -2D^H(P^*)^-1 W D^*,
and
g = 2 D^H (P^*)^-1 (z-Wz^*).



Taking the derivative of the Lagrange function (<ref>) with respect to x^*, and simplifying this, we obtain

    ∂ℒ(x,λ)/∂ x^* = -g + G_1 x + G_2 x^* + C^Hλ.

Equating this to zero and using augmented matrices and vectors, we get

    g̅ = G̅x̅ + C̅^Hλ̅.

Combining this with the constraint equation (<ref>) on augmented form, i.e. c̅ = C̅x̅, we obtain (<ref>).


Throughout this paper we will assume that the matrix 

    A=
    [   G̅ C̅^H;   C̅    O ]

appearing in (<ref>) is invertible. If this is not the case, then the solution to (<ref>) is not unique. Even if it is not unique, the solution may still provided valuable information, e.g. some entries in the estimate may still be unique, but we will leave this issue for future study. When A is invertible,
we can formulate the inverse of A as a block matrix with blocks of the same size as A, i.e.
    A^-1 = [ F̅_11 F̅_12; F̅_21 F̅_22 ],

and we can get a closed form solution for the maximum likelihood estimate of x. 


Assume (<ref>) and (<ref>). If A given by (<ref>) is invertible, the maximum likelihood estimate is given by

    x̂̅̂ = F̅_11g̅ + F̅_12c̅




This follows directly from Proposition <ref> by multiplying (<ref>) by (<ref>).

In principle closed form expressions can be obtained for F̅_11, F̅_12, etc. by using block inversion on A, but this requires G̅ to be invertible, and a necessary requirement for this is that D has full column rank. This is rarely fulfilled for the application in the present paper, and instead A can be inverted numerically, and the submatrices can then be extracted from A^-1. Note that A is a sparse matrix, which makes the numerical inversion faster. When A is invertible, the following proposition gives the distribution of the maximum likelihood estimate. 


Assume (<ref>) and (<ref>). If A is invertible, then

    x̂̅̂∼ CN(x̅, 2F̅_11).




First notice that g̅ = J̅z̅, where 

    J̅ = 
    (2D^H(P^*)^-1,-2D^H(P^*)^-1W)

and thus from Proposition <ref>, we get that x̂̅̂ = F̅_11J̅z̅ + F̅_12c̅ is a widely linear transformation of z̅. Hence from the distributional assumptions on z, i.e. (<ref>), we obtain that

    x̂̅̂∼ CN(F̅_11J̅D̅x̅ + F̅_12c̅, F̅_11J̅Σ̅J̅^H F̅_11).

To simplify the expectation, observe that expanding the identity A^-1A=I into its blocks, we obtain

    F̅_11G̅ + F̅_12C̅ = I   and  C̅F̅_11 = O.

The first equation in (<ref>) together with the observation that J̅D̅ = G̅ implies that 

    𝔼x̂̅̂ = F_11J̅D̅x̅ + F̅_12c̅ = (F_11G̅ + F̅_12C̅)x̅ = x̅.

To simplify the augmented covariance matrix, firstly observe that combining both equations in (<ref>), we obtain that




    F̅_11G̅F̅_11 = F̅_11.

Secondly, expanding and simplifying, we get that

    J̅Σ̅J̅^H = 2G̅.

From (<ref>) and (<ref>), it immediately follows that the covariance matrix of x̂̅̂ is given by 2F̅_11.


We note that Proposition <ref> implies that the maximum likelihood estimator is unbiased. Furthermore, it shows that the estimator is efficient, i.e. has a minimal covariance matrix, according to the following proposition, which gives the Cramer-Rao lower bound.


Assume (<ref>) and (<ref>). If A is invertible, then

    Cov(x̃̅̃) ≽ 2F̅_11,

where ≽ means that Cov(x̃̅̃) - 2F̅_11 is positive semi-definite, and x̃̅̃ is an arbitrary unbiased estimator of x̅.



By Theorem 1 in <cit.>, the Cramer-Rao lower bound in the complex valued and constrained case is given by

    Cov(x̃̅̃) ≽ U(U^H I(x̅) U)^-1 U^H,

where I(x̅) is the Fisher information matrix, and U is a 2N×(2N-2Q) matrix fulfilling 

    U^T U = I   and  (∂/∂x̅ h(x̅)) U = O

and h(x̅)=0 is the constraint. 

Since h(x̅)=C̅x̅ - c̅, the latter part of (<ref>) becomes C̅U=O, and thus the column space of the matrix U is the orthogonal complement of the row space of C̅, since A is invertible and therefore rank(C̅)=2Q implying rank(U)=2N-2Q, i.e. U has full rank. Now observe that since the 2N× 2N matrix F̅_11 fulfills the second equation in (<ref>), F̅_11 has the same column space, provided rank(F̅_11)=2N-2Q. We verify the rank: Firstly, by the second equation in (<ref>), we get that rank(F̅_11) ≤ 2N-2Q. Secondly, by the first equation in (<ref>), we get that

    2N ≤rank(F̅_11G̅) + rank(F̅_12C̅) 
    ≤rank(F̅_11) + rank(C̅)

so rank(F̅_11) ≥ 2N-2Q. Since U and F̅_11 thus have the same column space, we can use derivations similar to those immediately following (24) in <cit.> to get that 

    Cov(x̃̅̃) ≽F̅_11(F̅_11^H I(x̅) F̅_11)^+ F̅_11^H,

where we have first exchanged the inverse with a Moore-Penrose pseudoinverse shown as a ^+. 

By similar derivations as is used in (<ref>) for the first order derivatives of the Lagrange function, the second order derivative of the log likelihood function can be obtained, and from this we can obtain the Fisher information 

    I(x̅) = E(-∂/∂x̅∂/∂x̅^Hlog l) = 1/2G̅.

Inserting this into (<ref>) and using (<ref>), we get that 

    Cov(x̃̅̃) ≽2̅F_11(F̅_11^H G̅F̅_11)^+ F̅_11^H = 2F_11F̅_11^+ F̅_11^H = 2F̅_11.



We can obtain confidence intervals and ellipses from Proposition <ref> in the following way. Let F_1 and F_2 denote blocks in F̅_11=(F_1,F_2).
Then we get by Proposition <ref> that the covariance matrix for x̂_i (formulated as a two-dimensional real vector instead of a complex number) is given by

    [  ((F_1)_ii + (F_2)_ii) (-(F_1)_ii + (F_2)_ii);  ((F_1)_ii + (F_2)_ii)  ((F_1)_ii - (F_2)_ii) ].

From this, it follows that a confidence interval with confidence level 1-α for the real and imaginary parts of x_i are given by

    (x̂_i) ± z_α/2√(((F_1)_ii + (F_2)_ii)),
       (x̂_i) ± z_α/2√(((F_1)_ii - (F_2)_ii)).

Furthermore it follows that a confidence ellipse for x_i is given by the ellipse with 

    
  * center x̂_i, 
    
  * angle tan^-1(v_2/v_1),
    
  * major axis √(e_1χ_1-α,2^2), and    
    
  * minor axis √(e_2χ_1-α,2^2),

where e_1 is the largest eigenvalue of (<ref>), e_2 is the smallest eigenvalue, and (v_1,v_2) is the eigenvector corresponding to the largest eigenvalue.



§ ASSESSMENT APPROACH



Figure <ref> shows the details of how the estimator is applied and assessed.  We use a real-life distribution grid with smart meter measurements that determine the load scenario in a load flow calculation to determine the 'true' voltage and current phasors everywhere in the grid. This load-flow calculation together with the structural data and the load data is called the Reference Grid, see Figure <ref>. The load-flow calculation in the  RGM is based in a standard Newton-Raphson algorithm to solve the power flow equations, see e.g. <cit.>.



























































The true state vector x is then used by the Data Generator to generate measurements according to the PMU or EM model, shown in Figure <ref>. The preprocessed data z obtained from the measurements is then fed into the estimator from Section <ref>. 

The estimator provides the estimated state vector x̂ together with the confidence ellipses. This output is then compared with the true state vector x. The actual metrics for comparison will be provided in the next section together with the results.

The data generator requires as input the standard deviations of the voltage magnitude, of the current magnitude, of the phase angle, and of the voltage angle. The standard deviation of the phase angle, ϕ, between voltage and current phasors has been stated in the literature as σ_ϕ=10^-2 rad <cit.>, so this value is used in the case study later.

 The standard deviation of the voltage angle, σ_θ, is in our assessment approach obtained by calculating the grid scenario for a set of realistic load conditions and calculating the empirical standard deviation of the samples of θ over all grid locations and over all load conditions. By using actual measurement data as realistic load conditions we get σ_Θ= 0.003 rad later in the case study. 

The errors on voltage magnitude and current magnitude are specified by relative errors, ρ_U and ρ_I. For voltages, the relative error is defined in reference to the nominal voltage, which is U_nominal=400V later in the case study. For currents, the relative error is in relation to the true current magnitude at each measurement location, i.e. the calculation of the standard deviation of the measurement error for current magnitudes will depend on the actual current magnitude obtained from the true state vector x in the assessment approach in Figure <ref>.  The standard deviations σ of error for voltage magnitudes and current magnitudes are chosen in a way, that β=99% of the erroneous samples are within a fraction ρ of the nominal value (for voltages), respectively the true value μ (for currents), i.e.

    Pr(X ∈ (μ(1-ρ),μ(1+ρ)))=β.

For the magnitudes of voltages and currents later in the case study, we use respectively, ρ_U=1% and ρ_I=3%, <cit.>. 
Using the (1-β)/2=0.5% quantiles of the normal distribution, we 
get σ_U and σ_I from: 



    σ_U r_0 = μ_Uρ_U,    σ_I r_0 = μ_Iρ_I

with r_0=norminv((1+β)/2,0,1).

Note that in the general case, different values of ρ_U and ρ_I can be used for different measurement locations, e.g. reflecting different types of measurement devices, so these relative errors could be vectors. In the case study of the next section, we assume the use of identical measurement devices at the different locations, so the scalars, ρ_U and ρ_I introduced above are used at every measurement location. 
´



§ CASE-STUDY INTRODUCTION


To validate the methods we apply the assessment approach on a case study using a real-life Danish distribution grid <cit.>.



 §.§ Grid Scenario

 The grid covers a typical small sized Danish town, supplying energy to several households, a school, a church, and some local industry. The grid consists of a secondary 10kV:400V substation with 98 connected customers with varying load characteristics. The information about grid structure and cable types was obtained through automatic processing from the Geographic Information System at the DSO <cit.> and an automated mapping of cable types to cable parameters was realized based on available cable data sheets. Using interfaces to the automatically derived digital twin of the low-voltage grid <cit.>, the C matrix was derived, which contains the system equation sets in Equation (<ref>). 

Figure <ref> shows the grid topology as a graph; since there are no cross-connections between the feeders, a tree topology results. The root of the tree is the low-voltage substation busbar of the secondary transformer. Intermediate nodes represent junction boxes. The leaves of the tree are the Customer Connection Boxes (CCBs), which represent the handover points from the grid operator to the customer and are measured by Smart Meters in the given grid. Although the Smart Meters measure voltages and currents per-phase in this 3-phase grid, we here want to avoid to deal with the issue of wrong phase assignment and therefore we use a single-phase representation of this 3-phase grid and we use the appropriate formulas to do a phase aggregation of the measurements.























 §.§ Measurement locations

The base case for this paper is the actual real-life measurement scenario in which meters measure the voltages at the CCBs as well as the currents into these CCBs. In the PMU model (used for validation of the theory and as a comparison case) the meters measure the actual voltage and current phasors. In reality, the meters follow the EM model, see Section <ref>, i.e. they measure voltage magnitude, current magnitude and local phase angles between currents and voltages. The intermediate nodes in the grid topology, Junction Boxes (JB), are not measured in the measurement scenario (and neither in the real-life grid). The LV side of the secondary substation is measured in the real-life grid, but these measurements are only used for the reference grid model.


















The smart meter technology deployed in the case study grid allows to collect data in cycles of 6 hrs, with a time resolution of 15 minute intervals. The magnitudes of the voltage and currents are measured, as well as the local phase angle ϕ. We assume for simplicity that all measurements have successfully been aligned in time, while typically some clock inaccuracies up to 9 seconds are allowed, <cit.>, which would lead to additional errors in the measurand, see <cit.>. 




















    































	

		






		














These assumed true values of voltage and current phasors, shown for chosen example locations in the topology in Table <ref>, are later used as comparison base for the estimator and they are also input to the measurement models. The relative measurement errors, ρ_U and ρ_I, for voltages and currents are derived from the measurement device type (e.g. measurement device class) as described in the previous section. The standard deviations of voltage magnitude error and current magnitude error are then calculated according to Equation (<ref>). 
These values as well as a reduced set of true voltages and currents are used as input for the Data Generators, which create a noisy variant of the 'true' voltages and current phasors from the RGM according to the selected meter model (PMU or EM), see Section <ref>. 



































§ CASE STUDY RESULTS

We now apply the assessment approach to the grid and measurement scenario as described in the previous sections. Purposes of the case study are: (1) to validate the calculations and their implementation; (2) to investigate to what extent the approximations via pseudo-measurements (Sect. <ref>) and by the 'banana-shape' approximation (Sect. <ref>) for the EM model affect the accuracy of the confidence regions of the estimator; (3) to gain insights into the behaviour of confidence regions for the model estimation results for a practically relevant example case.  

Purpose (1) is achieved via applying the assessment approach using the PMU data generator to derive a so-called hit-rate metric, which is introduced in the next subsection. Reason to use the PMU Data Generator is that it exactly fulfils the theoretical assumptions of the estimator. 
Purpose (2) is addressed by an assessment of the hit-rate metric for the EM Data Generator, and its comparison to the PMU results. Purpose (3) is addressed for the EM Data Generator as practical measurement deployments in the use-case and in almost all of today's distribution grids follow the EM model.

Two metrics are used for the validation and assessment: 1) the average hit-rate, where the hit-rate is defined by the fraction of times that the estimated confidence ellipse from the estimator includes the true value, averaged over all grid nodes (for voltages) respectively grid lines (for currents):

    HR_k:= 1/R∑_r=1^R I(x_k,r⊆ ce_k,r)

The hit-rate for node k, HR_k is averaged over R repetitions, which we then average for all K measurement nodes to a system average hit-rate HR.


2) In order to quantify the convergence, we also look at deviation of the hit-rate estimator as follows: 

    Dev.HR_k=UB{HR_k}-LB{HR_k}

where Dev.HR_k expresses the deviation between the 95% upper (UB) and lower (LB) bounds of the hit-rate estimate for node k over R repetitions. Then, for a simpler handle of the metric, we average over all nodes to get the average hit rate deviation:


    Dev.HR:=1/KΣ_k^K Dev.HR_k.





 §.§ Validation - Model impact on hit-rates for the PMU Data Generator




The introduced grid estimation approach provides not only the estimates of voltage and current phasors but also confidence ellipses for the obtained phasors.  In order to validate the calculation of the confidence ellipses, the  grid scenario from the use-case is fixed and ground truth is obtained,  and a large number of R=50000 repetitions of sets of erroneous preprocessed observations are generated by the PMU data generator, and the estimation approach is calculated for each of these R preprocessed observation vectors. 

For the PMU Data generator, the hit-rate (HR) should stochastically converge to the value of α=95%, with α chosen arbitrarily to reflect a 95% confidence in the estimate, when the error model follows a 2D-Gaussian model. 




The left half of the result in Table <ref> confirms the convergence of the hit-rates to α=95% for the PMU Data Generator. The first row use a baseline assumption on the standard deviations of errors on voltages, currents, and local phase angle, while subsequent rows increase or decrease one of the error standard deviations by a factor of 10, keeping the other two on their baseline value. These standard deviations are used in the Data Generator and the derived covariance matrices Σ_1 and Σ_2 are channeled as input into the estimator; therefore the estimator is made aware of the change of the standard deviations in the Data Generators.










 §.§ Assessment of hit-rate when using the EM Data Generator

When the assessment uses the EM data generator, there are two approximations that deviate from the theoretical and accurate case: (1) the Pseudo-measurement for the voltage angle, which is taken in the data preparation in the EM data generator (Sect. <ref>); (2) the assumption in the estimator of Gaussian noise, which is only approximately true, see Sect. <ref>. These two approximations will influence the hit-rates of the confidence regions.

The right half of Table <ref> shows the impact of these two approximations, which depends on parameters of the measurement errors. In addition the grid scenario will influence the impact of the pseudo-measurements, so here the results have to be interpreted in the context of the realistic distribution grid of the use-case. 

The base case of measurement error standard devations, which is derived from realistic measurement device parameters, is specified in the top numerical row. For this base case, the hit-rates for voltage and current confidence regions are very close to the theoretical value of 95%, the voltage hit rate is 1% lower while the current hit-rate is slightly higher than 95%. When decreasing the standard deviations of the voltage or of the current, see third and second row from the bottom, respectively, the hit-rates stay at the base level approximately or even increase slightly. When decreasing the standard deviation of the local phase angle error between voltage and currents, both hit-rates decrease slightly - the latter is suspected to be an impact of the pseudo-measurements for voltage angles, which are suspected to take over more strongly when the standard deviation of the local phase angle error decreases.

When increasing the standard deviations by a factor of 10 individually (numerical rows 2 to 5 in table), hit-rates for voltages and currents both decrease, however to a different extent. Note that a factor of 10 on the standard deviations corresponds approximately to measurement errors of 
10%-30% (or even much more for current magnitudes, as the actual smaller values of the current magnitudes further out in the grid are not calculated in here when using a constant standard deviation), which is typically beyond practical relevant measurement deployments, so this is an extreme case to demonstrate the impact.

When increasing the standard deviation of the voltage magnitude error by a factor of 10, hit-rates for current estimator condidence regions are only marginally affected, while the hit-rates for the voltage estimator drop by almost 7%. When increasing the standard deviation of the current magnitude errors, we see the inverse situation: Voltage estimator hit-rates are only marginally reduced while the impact on current estimator hit-rate becomes noticable, while still being less than 1% reduced.
The increase of the standard deviation of the error on the local phase angle has the strongest impact: hit-rate of voltage phasors drop by almost 20% and hit-rates of current-phasors drop by a bit more than 3% compared to the base-case.

In summary, the realistic base case leads to hit-rates very close to the theoretical 95% while extreme scenarios of huge measurement errors may reduce these hit-rates strongly. The practical applicability of the novel approach to EM measurements in realistic settings is confirmed by the analysis.



 §.§ Voltage and current estimates under the EM model

In the subsequent analysis the EM data generator is used with the relative errors according to the base-case parameters (so from now on, standard deviations of errors on the current magnitude are in fact proportional to the true current magnitude of the currents at that measurement location, in contrast to the analysis in the previous section).

Figures <ref> to <ref> show the estimated voltage or current phasors (visualized by a square) in one repetition in comparison to the true value (diamond) and the estimated confidence ellipse. The first figures show the estimates at customer connection boxes, which also provide measurements (visualized by plus symbol) as input to the estimation approach. For the example results, we chose meters with IDs {70, 109, 134}, see Figure <ref> for grid location. 

The measurements of voltages in Figure <ref> are actually all lying on the x-axis, since the EM metering model does not allow to measure absolute phase angles - instead the phase angle θ of the voltage phasor is assumed to be 0 in the measurement. Figure <ref> shows the current estimates with confidence ellipses for the same meters.



















Figures <ref> and <ref> show a similar result for intermediate junction boxes, with ID's {42,44,106}, which do not provide any measurements to the estimation approach.




























 §.§ Analysis of noise on confidence regions

In grid monitoring we are ultimately interested in estimation of voltage and current magnitudes, and possibly on phase angles. In order to assess how these estimation results are impacted by different measurement error magnitudes, we define the 95% confidence range for voltage magnitudes as the difference of the largest and smallest voltage magnitude on the confidence ellipse Δ C, indexed either for voltage (U) or current (I): 

    Δ C_U =max_ConfEllipse |U| - min_ConfEllipse |U|.

We assess in the following, how these confidence ranges are impacted by changing error magnitudes.



  §.§.§ Change of error standard deviation on measurements of voltage magnitude


We first change the standard deviation, σ_U of the normally distributed error on the voltage magnitude. Figure <ref> shows that the confidence range for the estimated voltage magnitude for all 7 considered grid locations changes almost linearly with σ_u within the considered range of the voltage measurement error standard deviation between 1V and 4V. Furthermore, the difference of confidence ranges between different grid locations is small - and reduced further for larger σ_U. Although being measured by Smart Meters as input to the estimator, Customer Connection Boxes (dashed) show a slightly larger confidence range for the voltage as compared to junction boxes (non-measured, dotted) and the substation itself (also not measured for the estimation, solid).



A very different behavior is obtained for the confidence range of the current magnitudes, shown in Figure <ref>: Varying the standard deviation of the voltage measurements in the range between 1V and 4V only has negligible impact on the confidence ranges of the estimated currents, the shown curves are almost horizontal lines. In contrast to the voltage magnitudes, the confidence ranges of the current magnitudes however depends strongly on the grid location (here the line feeding into the substation bus bar, JB or CCB).
Furthermore, the different CCBs show different confidence ranges for the current magnitudes, partly resulting from the quite different true current phasors, see Table <ref>.




























  §.§.§ Change of error standard deviation on measurement of current magnitude


For this set of experiments, we move away from a fixed relative error on the current magnitude measurements across the measured CCBs, and instead set a fixed standard deviation, σ_I of the error on current magnitude across all measured grid lines. Figure  <ref> shows the impact on the confidence range of the voltage magnitude estimates, while Figure <ref> shows the impact on the confidence range of the estimated current magnitudes. Compared to varying the error standard deviation on voltage magnitudes in the previous subsection, the impact is more widely varying depending on inspected grid location. 









  §.§.§ Change of error standard deviation on measurement of local phase angles

When varying the error magnitude on the local phase angle measurement (here a range of σ_ϕ up to approximately double of its value in the base scenario is applied), results are not visibly impacted for the parameter ranges. Within the specific range of σ_ϕ∈[0,0.018] rad confidence ranges for voltage and current estimates do not change within four digits accuracy.

 




































We also investigated the confidence range of voltage and current magnitudes when visualizing over the grid topology graph. However, as the confidence ranges for voltage magnitudes show only little differences across different grid locations, and the confidence ranges for current magnitudes are largely influenced by the true current magnitude, the results are not interesting to be shown here. However, this will change when changing the measurement scenario, e.g. assuming some CCBs as not being measured, or adding the measurement at the substation. The latter studies have been started and will be published in subsequent work.












































§ CONCLUSION AND OUTLOOK

 

In this paper we presented a methodology to derive and assess confidence regions of estimated voltage and current phasors in an electrical distribution grid. Typical measurement devices in low and medium voltage grids are not clock synchronized to high precision level and therefore only allow to measure local phase angles. Therefore, the paper introduced the EM metering model, which is subsequently compared to a PMU metering model. While the PMU metering model fulfills the theoretical assumptions for the estimator, two approximations need to be introduced for the EM metering model: (1) constant zero-valued pseudo measurements on the phase angle of voltages; (2) approximations of the measurement error by 2-D Gaussian noise. 
Furthermore, not all desirable grid locations may be observed by measurement devices in distribution grids. 

In a realistic use-case of a single-phase representation of a real distribution grid, we showed that the estimation and related confidence ellipses show the theoretically expected stochastic convergence for the PMU model, while there are some deviation from the specified confidence level for a more realistic EM measurement model. However, these deviations are only significant in case of extremely large measurement errors.  A confidence range was used to asses the impact of increasing measurement errors on measured voltage magnitudes, current magnitudes and phase angle measurements, respectively, and example results from the realistic use-case have been reported.

Further reduced sets of measurement locations, e.g. in which some customers do not have smart meters, or some smart meters only measure voltages (which are less privacy sensitive), and also additional measurement device locations at the secondary substation will be investigated in future case studies. In this context, the generation of extended types of pseudo-measurements and the derivation of properties of the measurement error of such pseudo-measurements will also need to be studied. In addition, a characterization and inclusion of errors in the structural information and cable properties of the grid will be interesting. Finally, the extension of the use-case study to use a 3-phase plus neutral model will be interesting; the presented mathematics hold as the grid equations stay linear <cit.>, while additional complexity results from the wrong-phase assignment cases of measurements of voltages and currents.

Recent technological developments with 5G and beyond, and the installment of these communication interfaces in modern smart meters, enables a higher accurate clock synchronization among the meters. For example <cit.> discusses the requirements and use of Precision Time Protocols (PTP) to achieve clock synchronization as low as ±1.5μsec. for smart grid applications, and <cit.> concludes in 5G end-to-end systems clock synchronization can within a 5G-TSN network be kept within boundaries of 900ns. 5G enabled smart meters are, however, at this stage not widely deployed, but the industrial trend moves towards using such networks in next generation smart meters, hereby enabling applications with high requirement for clock synchronization such as DSSE. The advancement of technology may lead to the partial availability of voltage phasor measurements. Use-cases of heterogeneous measurement types, in the terminology of this paper, some measurement locations with PMU meters, others with EM meters, will therefore be interesting for future studies.



 



















































 elsarticle-num 
  





§ 



Let Z=(z+ϵ_z) e^j(ν+ϵ_ν), z,ν∈ℝ, and ϵ_z and ϵ_ν be independent random variables with 𝔼[ϵ_z] = 𝔼[ϵ_ν] = 0, Var(ϵ_z) = σ_z^2 and Var(ϵ_ν) = σ_ν^2. Then the mean, variance and pseudo-variance of Z are given by

    𝔼[Z] =  z e^jν c_ϵ_ν(1),    Var(Z) = (1-|c_ϵ_ν(1)^2|) z^2 + σ_z^2

and 

    PVar(Z) = e^2jν((z^2+σ_z^2)c_ϵ_ν(2) - z^2 c_ϵ_ν(1)^2),

where c_ϵ_ν(t) = 𝔼[e^jtϵ_ν] is the characteristic function of ϵ_ν. 




Using that ϵ_z and ϵ_ν are independent and 𝔼[ϵ_z]=0, we get that

    𝔼[Z] = 𝔼[z+ϵ_z]𝔼[e^j(ν+ϵ_ν)] = z e^jν c_ϵ_ν(1),


    𝔼[|Z|^2] = 𝔼[(z+ϵ_z)^2] = z^2 + σ_z^2

and

    𝔼[Z^2] = 𝔼[(z+ϵ_z)^2]𝔼[e^2j(ν+ϵ_ν)] = (z^2+σ_z^2) e^2jν c_ϵ_ν(2).

The formulas for the variance and pseudo-variance in the general case follows by inserting the expected values into Var(Z) = 𝔼[|Z|^2] - |𝔼[Z]|^2 and Var(Z) = 𝔼[Z^2] - 𝔼[Z]^2.  

 Let Z defined as in Proposition <ref> and let ϵ_ν be normally distributed. Then

    𝔼[Z] =  z e^jν e^-σ_ν^2/2,    Var(Z) = (1-e^-σ_ν^2) z^2 + σ_z^2

and 

    PVar(Z) = e^2jν((z^2+σ_z^2)e^-2σ_ν^2 - z^2 e^-σ_ν^2).



Use the characteristic function for the normal distribution in Proposition <ref>.


















