Smart Grid State Estimation with PMUs Time Synchronization Errors ^†

Abstract

State Estimation (SE) is one of the essential tasks to monitor and control the smart power grid. This paper presents a method to estimate the state variables combining the measurement of power demand at each bus with the data collected from a limited number of Phasor Measurement Units (PMUs). Although PMU data are usually assumed to be perfectly synchronized with the Coordinated Universal Time (UTC), this work explicitly considers the presence of time-synchronization errors due, for instance, to the actual performance of GPS receivers and the limited stability of the internal oscillator. The proposed algorithm is a recursive Kalman filter which not only estimates the state variables of the power system, but also the frequency deviations causing clock offsets which eventually affect the timestamps of the measures returned by different PMUs. The proposed solution was tested and compared with alternative approaches using both synthetic data applied to the IEEE 123 bus distribution feeder and real-field data collected from a small-size medium-voltage (MV) distribution system located inside the EPFL campus in Lausanne. Results show the validity of the proposed method in terms of state estimation accuracy. In particular, when some synchronization errors are present, the proposed algorithm can estimate and compensate for them.

1. Introduction

Presently, the worldwide ambient and climate targets could be achieved only through the so-called energy transition where renewable energies will play a central role to reduce the CO${}_2$ emissions [1]. The current electrical grid should become more flexible and robust to host not only variable Renewable Energy Sources (RES)—such as solar and wind generators—but also an increasing number of time-variable and non-linear loads such as electric vehicles and heat pumps. This context requires advanced monitoring and control techniques to assure proper and effective behavior to both transmission and distribution grids. Among the different tasks, the State Estimation (SE) is one of the most important tools which allows knowledge of the state variable in each point of the grid to properly adopt necessary control actions [2]. The classic method for SE solution is based on a weighted least-square (WLS) algorithm and has been mainly formulated and used in transmission networks due to their stable topology and numerous measurement points [3] which guarantee the observability requirement. However, recently, SE gained interest also at the distribution level where the physical structure and the lack of measurement points with respect to the number of nodes make observability sometimes challenging [4]. Also, due to the development of smart grid framework system-level solutions for SE that can be applied to both transmission and distribution grids can be of great interest [5], so different state estimators have been proposed in the recent literature for both transmission as well as distribution grids [6,7]. Complementary to this aspect, a new measurement device called Phasor Measurement Unit (PMU) is becoming central to improve the SE accuracy. Initially, thanks to their capability to provide a grid snapshot through either the voltage or current phasors synchronized with the Universal Coordinated Time (UTC), PMUs have been proposed for applications on the wide-area monitoring systems (WAMS) in transmission networks. According to their accuracy, PMUs can operate in two different classes (i.e., M or P) for measurement or protection scopes with different reporting rates depending on the use and the grid frequency (e.g., 60 Hz in the USA or 50 Hz in Europe) [8]. These characteristics allow the collection of numerous synchrophasors per second which can be used to track the distribution grid dynamics caused by time-varying loads or generators (as for example PV generators, EV-charger or storage charging and discharging ) [9], or for power quality evaluation [10], or for protection and stability assessment [11] or for State estimation process [12,13,14]. In the recent literature, different works analyze the impact of measurement uncertainty on SE accuracy [15,16] highlighting that a low number of PMU can significantly increase the state variable estimation even in the presence of distributed generation [17]. For this purpose several SE methods have been proposed only based on synchrophasors as in [12] or combining the measurements provided by the traditional power devices with them [18,19,20].

It is important to highlight that despite their high accuracy with respect to the classical power meters, the PMUs are affected by three different sources of uncertainty, i.e., synchrophasor estimation algorithms, analog-to-digital converter (ADC), or synchronization process [21]. In this paper, we only focus on the time-synchronization process and its impact on state estimation. As reported in [22] the synchronization module is responsible for the time dissemination over the system and the synchronization of the internal clocks for acquisition and the time-tag of the performed measurements. Internally to PMU the clocks, can be disciplined by the time-synchronization source and widespread all over the architecture. To discipline the acquisition with the clock, the most typical technique is based on phase-locked-loop (PLL). This operation can cause frequency deviations producing not negligible time-skews. This phenomenon has been considered in the literature related to the SE problem. In particular in [23] the impact of the synchronization uncertainty using a static estimation algorithm has been studied. A possible solution based on a bilinear measurement model considering the grid state variable and synchronization parameters is presented in [24]. The algorithm consists of two parallel Kalman filters which solve at the same time the state and the synchronization inaccuracy. Based on similar model, results are also reported in [25]. The synchronization source can rely directly or indirectly from GPS receiver. When the source is external the information about the reference time can be distributed using network time methods such as net time protocol (NTP), precise time protocol (PTP) [26], white rabbit [27] or IRIG-B. This type of synchronization is particularly adopted in substation to synchronize PMU from different points respecting the accuracy of ±1 $\mathsf{\mu }$s required by the Power Profile Standard IEEE C37.238–2011 [28]. However also the GPS or the other timing system can be affected by errors, because the signals can be not reachable or corrupted at the source at the time signal receiver or within the PMU [29] and this represent a further synchronization inaccuracy.

Since the accuracy of the PMU affects considerably the accuracy of the SE, the aim of this work is to consider PMU measurement uncertainties coming from time synchronization and to propose a technique to simultaneously estimate the state variables and the synchronization parameters in order to possibly compensate the errors. With respect to this context, this paper aims to extend the work done in [30] and to present a real-time algorithm based on Kalman filter. The algorithm is based on the linear model proposed in [31], where the power flow manifold is approximated around any feasible working point and explicitly computed offline. Based on the work published in [18], the state estimator will use two sources of information: the power profiles of active and reactive power (coming from forecasted values or pseudo-measurements) and the PMUs. Also, as a byproduct of the model choice, our methodology seamlessly applies to transmission as well as distribution grids. The benefits of the present method is evaluated through simulations over two different grids and the results show that in the presence of synchronization error, the presented method is able to compensate the inaccuracy and perform better than other algorithms that do not consider synchronization error parameters in the model. Moreover, a test over a small real grid only monitored by PMUs highlights the capability for the algorithm to compensate the frequency deviation within the nodes.

The remainder of this paper is structured as follows. The end of this section is devoted to symbols and abbreviations. Section 2 presents the power grid modeling and the power flow linearization. Section 3 formalizes the model the two type of measurements used (i.e., power data and PMUs). Section 4 describes the state estimation and Section 5 a two-node example theoretical solution. In Section 6 simulation and results are reported and finally the conclusions are presented in Section 7.

2. Grid Model and Linearized Power Flow

As aforementioned, the electric grid is evolving toward the smart grid concept due to the deep deployment of instrumentation devices and the use of ICT for monitoring and control. With this respect the power grid can be considered to be a cyber-physical system consisting of three layers which from top to bottom are given by (i) the cybernetic layer, (ii) the metering layer and (iii) the physical layer. See Figure 1 for a graphical representation of a smart grid as a cyber-physical system. The cybernetic layer consists of two main actors: the data aggregator (DA) and the central processing unit (CPU). The DA gathers the measurements coming from the metering layer and sends it to the CPU which performs all the computations.

Among the different tasks performed by the CPU, in this paper we focus on the estimation of the power grid state at a certain time. To estimate the state of a grid, it is necessary that the grid is observable, and this is assured by the second layer through the metering infrastructure. In this specific case, we assume that at each bus we can reconstruct the power profile based on low-accurate power meters. Moreover, we consider that a small subset of buses is equipped with Phasor Measurement Units (PMUs) which can provide the bus voltage and current phasors synchronized with the UTC. For all the details about the measurements model we refer the reader to Section 3.

Finally, regarding the physical layer, under the hypothesis of sinusoidal steady-state condition, if we denote with $\mathcal{V}=\{1,\dots ,n\}$ the electrical buses and $\mathcal{E}$ the lines between connected buses, the grid can be modeled as a graph $\mathcal{G}(\mathcal{V},\mathcal{E})$. For each node $h\in \mathcal{V}$, it is assumed that:

• $u_h=v_h{e}^{j{\mathit{\theta }}_h}\in \mathbb{C}$ is the complex voltage at the bus h where $v_h,{\mathit{\theta }}_h\in \mathbb{R}$ are the modulus and phase of the complex phasor, respectively;
• $i_h\in \mathbb{C}$ is the complex current injected at the bus h;
• $s_h=p_h+j{q}_h$ is the apparent power absorbed by the bus where $p_h,q_h\in \mathbb{R}$ are the active and the reactive power, respectively.

In addition, if we consider two nodes, h and k, the element $y_{hk}$ is the admittance of the electric line $(h,k)$, while $y_h^{sh}$ is the shunt admittance (admittance to ground) at bus h. Because of that, the admittance matrix can then be written as:

$${\left[Y\right]}_{hk}=\left\{\begin{array}{c}\begin{array}{cc}{y}_h^{sh}+\sum _{\ell \ne h}{y}_{h\ell },\hfill & \mathrm{if}k=h;\\ -y_{hk},\hfill & \mathrm{otherwise};\end{array}\hfill \end{array}\right.$$

(1)

Defining the voltage, current and apparent power vectors as $\mathbf{u}={[u_1,\dots ,u_n]}^T$, $\mathbf{i}={[i_1,\dots ,i_n]}^T$, $\mathbf{s}={[s_1,\dots ,s_n]}^T$, respectively, we have the Kirchhoff’s law and the nodal power balance as

$$\mathbf{i}=Y\mathbf{u},\mathbf{s}=\mathrm{diag}\left(\mathbf{u}\right)\overline{\mathbf{i}}$$

(2)

where $\overline{(·)}$ denotes the complex conjugate operator and $\mathrm{diag}(·)$ denotes the diagonal matrix with $ii$-th diagonal element equal to the i-th element of its vector argument. From the previous equations one gets

$$\mathbf{s}=\mathrm{diag}\left(\mathbf{u}\right)\overline{Y\mathbf{u}},$$

(3)

that equals the power flow equations to be solved. At this point, based on the main result of [31] it is possible to linearize the non-linear power flow Equations (3) around any feasible point in the power flow manifold.

The grid state vector can be defined as $\mathit{\xi }:={[{\mathbf{v}}^T,{\mathit{\theta }}^T,{\mathbf{p}}^T,{\mathbf{q}}^T]}^T$, being $\mathbf{v},\mathit{\theta },\mathbf{p},\mathbf{q}\in {\mathbb{R}}^n$ the bus values of voltage magnitude and phase ($\mathbf{v},\mathit{\theta }$) and active and reactive power ($\mathbf{p},\mathbf{q}$), respectively. Then, by expressing the complex Equations (3) in rectangular coordinates, it is possible to rewrite them in implicit form as $F\left(\mathit{\xi }\right)=0$, $F:{\mathbb{R}}^{4n}\mapsto {\mathbb{R}}^{2n}$, and, in turn, implicitly define (Lemma 1 of [31]) the power flow manifold

$$\mathcal{M}:=\left\{\mathit{\xi }\right|F\left(\mathit{\xi }\right)=0\}.$$

(4)

[Proposition 1 of [31]] Let define ${\mathit{\xi }}^{*}\in \mathcal{M}$, i.e., ${\mathit{\xi }}^{*}=\left\{{[{\left({\mathbf{v}}^{*}\right)}^T,{\left({\mathit{\theta }}^{*}\right)}^T,{\left({\mathbf{p}}^{*}\right)}^T,{\left({\mathbf{q}}^{*}\right)}^T]}^T\right|F\left({\mathit{\xi }}^{*}\right)=0\}.$ Then, the linear manifold tangent to $\mathcal{M}$ in ${\mathbf{x}}^{*}$ is given by

$$A_{{\mathit{\xi }}^{*}}(\mathit{\xi }-{\mathit{\xi }}^{*})=0,$$

(5)

where

$$A_{{\mathit{\xi }}^{*}}=\left[\underset{A_{{\mathbf{u}}^{*}}}{\underbrace{\left(\langle \mathrm{diag}\overline{Y{\mathbf{u}}^{*}}\rangle +\langle \mathrm{diag}{\mathbf{u}}^{*}\rangle N\langle Y\rangle \right)R\left({\mathbf{u}}^{*}\right)}}-I\right],$$

${\mathbf{u}}^{*}:={\mathbf{v}}^{*}{e}^{j{\mathit{\theta }}^{*}}$, I is the identity matrix of suitable size and

$$\begin{array}{cc}\hfill & N:=\left[\begin{array}{cc}I& 0\\ 0& -I\end{array}\right],\langle A\rangle =\left[\begin{array}{cc}\mathrm{Re}A& -\mathrm{Im}A\\ \mathrm{Im}A& \mathrm{Re}A\end{array}\right],\hfill \\ \hfill & R\left(\mathbf{u}\right):=\left[\begin{array}{cc}\mathrm{diag}(cos\mathit{\theta })& -\mathrm{diag}\left(\mathbf{v}\right)\mathrm{diag}(sin\mathit{\theta })\\ \mathrm{diag}(sin\mathit{\theta })& \mathrm{diag}\left(\mathbf{v}\right)\mathrm{diag}(cos\mathit{\theta })\end{array}\right].\hfill \end{array}$$

The best linear approximation corresponding to the plane tangent to $\mathcal{M}$ at ${\mathit{\xi }}^{*}$, of the power manifold $\mathcal{M}$ at the feasible point ${\mathit{\xi }}^{*}$ is stated by Proposition 2. It is relevant to note that assuming $A_{{\mathbf{u}}^{*}}$ invertible (this is not restrictive due to the presence of node shunt admittance), the voltage deviations in polar coordinates $\delta \mathbf{v}:=\mathbf{v}-{\mathbf{v}}^{*}$ and $\delta \mathit{\theta }:=\mathit{\theta }-{\mathit{\theta }}^{*}$ can be expressed as a linear function in rectangular coordinates of the power deviations $\delta \mathbf{p}:=\mathbf{p}-{\mathbf{p}}^{*}$ and $\delta \mathbf{q}:=\mathbf{q}-{\mathbf{q}}^{*}$

$$\left[\begin{array}{c}\delta \mathbf{v}\\ \delta \mathit{\theta }\end{array}\right]=A_{{\mathbf{u}}^{*}}^{-1}\left[\begin{array}{c}\delta \mathbf{p}\\ \delta \mathbf{q}\end{array}\right].$$

(6)

It is worth noting that Equation (4) holds for the electrical quantities (i.e., voltages and power injection) independently from the considered grid model buses (e.g., PQ, PV or slack bus). Indeed, the linear formulation in [31] works for any working point ${\mathit{\xi }}^{*}\in \mathcal{M}$ and, as shown in [31], generalizes over previously proposed linear approximation presented in the literature (such as the Linear Coupled power flow model [32], the DC power flow model [33] and the rectangular DC power flow model [34]). All the mathematical details are reported in [31].

3. Measurement Models

This section presents the measurement models used as input for the proposed state estimator. Since the grid observability is a requirement to find a possible solution in the state estimation problem, the first source of information is given by the active and reactive power at every bus, that will be used in a prior stage of the Kalman filter. Since these data are not everywhere available, at least as direct measurements, they are usually replaced by historical or forecasted data coming from aggregated values while considering a not negligible uncertainty contribution. As commonly indicated in the literature, we will also refer to this information as pseudo-measurements. This prior estimate will be refined in the update stage of the Kalman filter by exploiting high-accurate information coming from PMUs placed in (possibly only few) strategic points of the grid [18].

3.1. Power Demand

Currently, the electrical grid (mainly at distribution level) lacks measurement points to guarantee the necessary observability requirement for monitoring and control actions. For this reason, it is common to use information based on historical time-series or forecasted data of active and reactive power demands at each bus extrapolated from the aggregation of smart meters data [35]. By assuming node $h\in \mathcal{V}$ and time $t\in {\mathbb{Z}}_{+}$ it is possible to write

$$\left[\begin{array}{c}{\tilde{p}}_h\left(t\right)\\ {\tilde{q}}_h\left(t\right)\end{array}\right]=\left[\begin{array}{c}{p}_h\\ q_h\end{array}\right]+\left[\begin{array}{c}{w}_h^p\left(t\right)\\ w_h^q\left(t\right)\end{array}\right],\left[\begin{array}{c}{w}_h^p\left(t\right)\\ w_h^q\left(t\right)\end{array}\right]\sim \mathcal{N}(0,{\Sigma }^w)$$

(7)

where

$${\Sigma }^w=\left[\begin{array}{cc}{\sigma }_p^2{\left|p_h\right|}^2& \eta {\sigma }_p{\sigma }_q|p_h\left|\right|q_h|\\ \eta {\sigma }_p{\sigma }_q|p_h\left|\right|q_h|& {\sigma }_q^2{\left|q_h\right|}^2\end{array}\right]$$

being $p_h,q_h$ the nominal values of the active and reactive power respectively, ${\sigma }_p={\sigma }_q\approx$ 30–50% based on the results in [35] and $\eta \in [0,1]$. Also, we assume $\mathbb{E}\left[w_k^p\left(t\right)w_h^p\left(t\right)\right]=\mathbb{E}\left[w_k^q\left(t\right)w_h^q\left(t\right)\right]=\mathbb{E}\left[w_k^p\left(t\right)w_h^q\left(t\right)\right]=0$ according to [18,35].

3.2. Phasor Measurement Units

The second source of information is given by synchrophasor measurements provided by PMUs deployed in a limited set of nodes of the grid. As aforementioned, these devices are commonly considered able to perfectly synchronize the voltage or current phasor, with a negligible uncertainty synchronization contribution, by means of a GPS module for example. However, despite this, PMUs can be affected by poor synchronization due for example to temporary occlusion of satellites [36] or cyber-attacks [37]. In addition to this, within successive synchronization instants with the GPS module, usually providing 1pps (pulse-per-second) synchronization signal, PMUs exploit an internal oscillator as a reference clock which, in turn, might cause additional synchronization error. Hence, depending on the type of GPS module and oscillator, different synchronization error, directly proportional to the cost of these devices, can be achieved.

Ultimately, the measurements at bus h at time $t\in {\mathbb{Z}}_{+}$ can be denoted as

$$\begin{array}{cc}\hfill & {\tilde{v}}_h\left(t\right)=v_h\left(t\right)+w_h^v\left(t\right),w_h^v\left(t\right)\sim \mathcal{N}(0,{\sigma }_{\mathrm{pmu},v}^2|v_h{|}^2),\hfill \\ \hfill & {\tilde{\theta }}_h\left(t\right)={\theta }_h\left(t\right)+w_h^{\theta }\left(t\right)+d_h\left(t\right),w_h^{\theta }\left(t\right)\sim \mathcal{N}(0,{\sigma }_{\mathrm{pmu},\theta }^2),\hfill \end{array}$$

(8)

where we set ${\sigma }_{\mathrm{pmu},v}=0.1\%$, ${\sigma }_{\mathrm{pmu},\mathit{\theta }}={10}^{-3}\left[\mathrm{rad}\right]$, and assume uncorrelated measurement noise within the same node and across different nodes, i.e., $\mathbb{E}\left[w_k^v\left(t\right)w_k^{\theta }\left(t\right)\right]=\mathbb{E}\left[w_k^v\left(t\right)w_h^v\left(t\right)\right]=\mathbb{E}\left[w_k^{\theta }\left(t\right)w_h^{\theta }\left(t\right)\right]=0$. Indeed, $w^{\theta }$ is mainly due to sampling jitter and synchronization error while $w^v$ to the instrumentation amplitude noise. Differently from standard Gaussian additive models, the additional term $d_h\left(t\right)$ in (8) indicates the uncertainty contribution with respect to the UTC reference. It is known from [21] that this component mainly affects the voltage phase angle measurements, so we assume the synchronization error is limited to the voltage state variable. At this point, it is worth noting how the linear model (6) conveniently expresses the voltage state variables in polar coordinates as functions of the power state variables in rectangular coordinates. This differs from for previously proposed linear models usually expressing both voltages and powers (or, more often, currents) in rectangular coordinates. Because of this, the synchronization error should be projected from polar to rectangular coordinates thus introducing a coupling term between the real and imaginary part as well as an additional source of uncertainty (due to the approximate projection). Conversely, Equation (6) lets us avoid this.

The model of clock delay within successive synchronization instants ($kT,(k+1)T$), can be expressed as:

$$d_h\left(t\right)={\beta }_h+{\alpha }_h\frac{T}{M-1}t,t\in \{0,\dots ,M-1\},$$

(9)

where T is the GPS synchronization period and M is the number of PMU measurements collected within two successive synchronization instants. The parameters ${\beta }_h,{\alpha }_h\in \mathbb{R}$ are an offset term due to GPS inaccuracy, and the clock skew due to the fact that the internal clock of the PMU, in general, does not oscillate at the reference frequency, respectively.

Observe that the measurements are characterized by three different time-scales. First, power demand information coming from pseudo-measurements are usually available from historical data or in the form of 1-day a-head predictions. Second, the PMU measurements live on a faster time-scale, $t,t+1,\dots$, depending on the PMU reporting rate. Finally, as GPS provides a sync signal every T[s], the PMU internal clock re-synchronize with the universal reference at $kT,(k+1)T,\dots$. To distinguish the behavior when PMU measurements are collected with respect to instants when the GPS re-synchronizes, it is possible to re-define discrete-time instants in a unique time reference given by:

$$\tau (k,t)=kT+\frac{T}{M-1}t,k\in {\mathbb{Z}}_{+},t\in \{0,\dots ,M-1\},$$

(10)

referring to the k-th resync instant and to the t-th measurements within $[kT,(k+1\left)T\right)$, see Figure 2. For simplification in the following $x(k,t)$ will be used instead $x\left(\tau \right(k,t\left)\right)$

4. State Estimation

Literature presents different methods to solve the state estimation problem, which can be classified as static or dynamic. The static algorithms, mainly based on WLS algorithm, are only based on data measured at a given time. Conversely, the dynamic state estimators rely on a system model describing state evolution and estimate the state in an iterative way by applying a prediction and an update step. Depending on how the system model is defined, many different kinds of dynamic state estimators have been proposed for power systems, e.g., linear Kalman filters (KFs) [12,38], Extended Kalman Filters (EKFs) [39], or robust Unscented Kalman filters (UKF) [40]. In our case, by taking advantage of the dynamic linear model (6), the measurements and their statistical information, we present a real-time Kalman-based state estimator which, conversely to standard procedures, explicitly considers the effect of the synchronization error in the phasor measurement units. The proposed method is regarded as a Bayesian inference process, where active and reactive power measurements are used to establish a prior for the Bayesian model while the estimation improvement relies on PMUs information.

4.1. State-Space Model

Since the proposed estimator is based on a Kalman filter [41], a suitable state-space model is required. Conventionally, in the power system literature, since the grid state are the voltage phasors in rectangular coordinates, those are also chosen as filter state. Conversely, in this work, the state variables are the power demand deviations ( $\delta \mathbf{p},\delta \mathbf{q}\in {\mathbb{R}}^n$) and the synchronization error parameters ($\mathit{\alpha },\mathit{\beta }\in {\mathbb{R}}^m$). This choice is not restrictive and comes naturally from (6) and (8). It is relevant to remember that the synchronization error parameters equal the number $0\le m\le n$ of PMUs deployed in the grid. Now, since the incremental linear model (6) is defined regarding a predefined operating point, we assume ${\mathbf{p}}^{*}$ and ${\mathbf{q}}^{*}$ are nominal demands to which correspond ${\mathbf{u}}^{*}=({\mathbf{v}}^{*},{\mathit{\theta }}^{*})$. Then, state and output equations of the Kalman filter at $\tau (k,t)$, $k\in {\mathbb{Z}}_{+},t\in \{0,\dots ,M-1\}$, are defined by:

$$\begin{array}{ccc}& & \mathbf{x}(k,t+1)=\mathbf{x}(k,t)+{\mathbf{w}}^x(k,t)\hfill \end{array}$$

(11)

$$\begin{array}{ccc}& & \mathbf{y}(k,t)=H\mathbf{x}(k,t)+{\mathbf{w}}^y(k,t)\hfill \end{array}$$

(12)

where

$$\begin{array}{ccc}& & \mathbf{x}(k,t)={\left[\delta \mathbf{p}{(k,t)}^T\delta \mathbf{q}{(k,t)}^T\mathit{\alpha }{(k,t)}^T\mathit{\beta }{(k,t)}^T\right]}^T,\hfill \\ & & \mathbf{x}(k,0)\sim \mathcal{N}(0,{\Sigma }_0),\hfill \\ & & {\Sigma }_0=\left[\begin{array}{cccc}{\sigma }_p^2\mathrm{diag}\left(\right|{\mathbf{p}}^{*}{\left|\right)}^2& {\Sigma }_0^{pq}& 0& 0\\ {\Sigma }_0^{qp}& {\sigma }_q^2\mathrm{diag}\left(\right|{\mathbf{q}}^{*}{\left|\right)}^2& 0& 0\\ 0& 0& {\sigma }_{\alpha }^{2}I& 0\\ 0& 0& 0& {\sigma }_{\beta }^{2}I\end{array}\right],\hfill \\ & & {\Sigma }_0^{pq}={\Sigma }_0^{qp}={\sigma }_p{\sigma }_q\mathrm{diag}\left(\right|{\mathbf{p}}^{*}\left|\right)\mathrm{diag}\left(\right|{\mathbf{q}}^{*}\left|\right),\hfill \\ & & {\mathbf{w}}^x(k,t)\sim \mathcal{N}(0,W),\hfill \\ & & \mathbf{y}(k,t)=\left[\begin{array}{c}\delta \tilde{\mathbf{v}}(k,t)\\ \delta \tilde{\mathit{\theta }}(k,t)\end{array}\right]:=\left[\begin{array}{c}\tilde{\mathbf{v}}(k,t)\\ \tilde{\mathit{\theta }}(k,t)\end{array}\right]-\left[\begin{array}{c}{\mathbf{v}}^{*}\\ {\mathit{\theta }}^{*}\end{array}\right],\hfill \\ & & \tilde{\mathbf{v}}(k,t)={[{\tilde{v}}_1(k,t)\dots {\tilde{v}}_m(k,t)]}^T,\hfill \\ & & \tilde{\mathit{\theta }}(k,t)={[{\tilde{\theta }}_1(k,t)\dots {\tilde{\theta }}_m(k,t)]}^T,\hfill \\ & & H=\left[\begin{array}{cc}{A}_{{\mathbf{u}}^{*}}^{-1}& \begin{array}{cc}0& 0\\ t\frac{T}{M-1}I& I\end{array}\end{array}\right],\hfill \\ & & {\mathbf{w}}^y(k,t):={\left[{\mathbf{w}}^v{(k,t)}^T{\mathbf{w}}^{\theta }{(k,t)}^T\right]}^T\sim \mathcal{N}(0,R),\hfill \\ & & R=\left[\begin{array}{cc}{\sigma }_{\mathrm{pmu},v}^2\mathrm{diag}\left(\right|{\mathbf{v}}^{*}{\left|\right)}^2& 0\\ 0& {\sigma }_{\mathrm{pmu},\theta }^{2}I\end{array}\right].\hfill \end{array}$$

The estimated process noise covariance information (e.g., from data) can be embedded in the W matrix. In the following, given the small re-synchronization period, we usually consider $W=0$ and outline an interesting closed-form analysis in Section 5. Since the above model is incremental with respect to the nominal value ${\mathit{\xi }}^{*}$ and since the synchronization error parameters can assume both positive and negative values, the state is reasonably initialized as a zero mean Gaussian random variable.

As aforementioned, it is worth observing that due to the linear relation (6) between buses power $\mathbf{p},\mathbf{q}$ expressed in rectangular coordinates and voltage $\mathbf{v},\mathit{\theta }$ expressed in polar coordinates, the synchronization error enters linearly in the output model (12) without any further approximation. This is opposed to standard approaches in the literature where, to deal with linear models, the grid state is expressed in rectangular coordinates, i.e., real and imaginary parts of the voltages. In this case, to resort to linear output models, the synchronization error must either be assumed or approximated as purely imaginary [18,23], under the additional assumption of small voltage angles differences (which is particularly true in distribution grids [42]). Also, even in the case when no synchronization error is considered, i.e., $d\left(t\right)=0$, observe that phasorial measurements are practically collected in polar coordinates. Hence, by formulating the output model with the same representation, we do not introduce any further manipulation of the data, i.e., projection from polar to rectangular coordinates, which, in turn, requires re-computation of the measurements correlation.

4.2. Synchronization-Aware State Estimator

Based on Equations (11) and (12) the model is linear and it is possible to build a Kalman filter [41] to simultaneously estimate grid state and synchronization error parameters. Algorithm 1 describes the proposed Synchronization-aware State Estimator, also referred to as SASE (code available at [43]).

Algorithm 1 SASE

Require: ${\Sigma }_0$, R, H. Initialize $\Sigma \left(0\right)={\Sigma }_0$. 1: for$t\in \{0,\dots ,M-1\}$do {//Offline} 2:  Compute and store $$\begin{array}{c}L(t+1)=(\Sigma \left(t\right)+W)H^T{(H(\Sigma \left(t\right)+W)H^T+R)}^{-1}\\ \Sigma (t+1)=(I-L(t+1\left)H\right)(\Sigma (t)+W)\hfill \end{array}$$ 3: end for 4: for$k\in {\mathbb{Z}}_{+}$do {//Online} 5:   Initialize $\widehat{\mathbf{x}}(k,0)=0$ 6:    for$t\in \{0,\dots ,M-1\}$do 7:     $\widehat{\mathbf{x}}(k,t+1)=\widehat{\mathbf{x}}(k,t)+L(t+1)(\mathbf{y}(k,t+1)-H\widehat{\mathbf{x}}(k,t))$ 8:   end for 9: end for

Observe that to run Algorithm 1 values for ${\mathbf{p}}^{*}$, ${\mathbf{q}}^{*}$, ${\mathbf{v}}^{*}$ and ${\mathit{\theta }}^{*}$ are required to compute $\mathbf{y}$ and H according to Equations (11) and (12). Also, it requires ${\Sigma }_0$ and R. By leveraging the information coming from the available power demand time-series, ${\mathbf{p}}^{*}$ and ${\mathbf{q}}^{*}$ are computed as one-day a-head predictions. Then, by means of a single full AC power flow computation it is possible to compute the corresponding values for ${\mathbf{v}}^{*}$ and ${\mathit{\theta }}^{*}$. Please note that since $\Sigma \left(t\right)$ does not depend on the measurements, its evolution can be computed offline and stored for $t\in \{0,\dots ,M-1\}$ (line 2) thus alleviating the computational burden. Then, at each time instant, standard Kalman equations can be applied (lines 6÷8) Finally, thanks to Equation (6), the estimated voltages are equal to

$$\left[\begin{array}{c}\widehat{\mathbf{v}}\\ \widehat{\mathit{\theta }}\end{array}\right]=\left[\begin{array}{c}{\mathbf{v}}^{*}\\ {\mathit{\theta }}^{*}\end{array}\right]+A_{{\mathbf{u}}^{*}}^{-1}\left[\begin{array}{c}\delta \widehat{\mathbf{p}}\\ \delta \widehat{\mathbf{q}}\end{array}\right],$$

(13)

where, by partitioning $\Sigma$ as ${\Sigma }_0$, the covariance is given by

$${\Sigma }^{\mathbf{u}}:=\left[\begin{array}{cc}{\Sigma }^v& {\Sigma }^{v\theta }\\ {\Sigma }^{\theta v}& {\Sigma }^{\theta }\end{array}\right]=A_{{\mathbf{u}}^{*}}^{-1}\left[\begin{array}{cc}{\Sigma }^p& {\Sigma }^{pq}\\ {\Sigma }^{qp}& {\Sigma }^q\end{array}\right]A_{{\mathbf{u}}^{*}}^{-T}.$$

As a side note, observe that both the model (11)–(12) and Algorithm 1 are outlined for $t\in \{0,\dots ,M-1\}$ for a given $k\in {\mathbb{Z}}_{+}$. As suggested by Figure 2, since at $\tau (k,0)$ the PMUs re-synchronizes with the GPS, the filter is re-initialized to reset $\mathit{\alpha }$ and $\mathit{\beta }$ and allows the computation of a new estimate (line 5 in Algorithm 1). Similarly, newly available ${\mathbf{p}}^{*}$, ${\mathbf{q}}^{*}$ and new data can be used to recompute the model and W, respectively. However, it is not mandatory to reinitialize the filter at every time window. For instance, in Section 6.2 we test the proposed SASE against the SoA algorithm proposed in [12] on measurements data collected from a 5-node distribution feeder located inside the EPFL campus in Lousanne. Equations (19)–(21) outline the state matrices used in that specific use case.

5. Two-Nodes Case: A Closed-Form Solution

Consider a two-nodes network consisting of one load connected to one generator (the PCC, $v_{\mathrm{pcc}}=1$, ${\theta }_{\mathrm{pcc}}=0$) through a purely inductive line with susceptance $b=-1$ [p.u.], in the absence of shunt admittance. For the sake of the analysis, we assume the load is absorbing only active power p while $q=0$. In this case, the flat profile is a particular solution which can be chosen as linearization point. Thus, by leveraging the linear model (5) one has $p=\theta$ and $v=1$ being v and $\theta$ the voltage magnitude and the phase at the load, respectively. Notice that since it is assumed $q=0$, the voltage magnitude is fixed and equal to 1 thus only $p=\theta$ is of interest. Now, assume to collect, within two successive GPS synchronization instants, M phase measurements of the form (8) which, as in (12), can be expressed as

$$\mathbf{y}=\left[\begin{array}{c}{y}_1\\ ⋮\\ y_{M-1}\end{array}\right]=\left[\begin{array}{ccc}1& 1& 0\\ 1& 1& \frac{T}{M-1}\\ ⋮& ⋮& \\ 1& 1& T\end{array}\right]\left[\begin{array}{c}\theta \\ \beta \\ \alpha \end{array}\right]+\left[\begin{array}{c}{w}_1\\ ⋮\\ w_{M-1}\end{array}\right]=C\mathbf{x}+{\mathbf{w}}^y,$$

(14)

with ${\mathbf{w}}^y\sim \mathcal{N}(0,R)$, $R={\sigma }_{\mathrm{pmu},\theta }^{2}I$, ${\mathbf{x}}_0\sim \mathcal{N}(0,{\Sigma }_0)$, ${\Sigma }_0=\mathrm{diag}({\sigma }_{\theta }^2,{\sigma }_{\beta }^2,{\sigma }_{\alpha }^2)$. Furthermore, for the sake of the analysis, let us assume absence of process noise, i.e., $w^x=0,W=0$ which, in the case of reasonably stable power demands within T[s], represents an acceptable first order approximation. Then, the posterior variance matrix in information form reads as

$$\Sigma ={\left({\Sigma }_0^{-1}+C^T{R}^{-1}C\right)}^{-1}$$

(15)

and, thanks to (14), after some tedious but straightforward algebraic manipulations, it is possible to compute $\Sigma =\Sigma ({\sigma }_{\mathrm{pmu},\theta },{\sigma }_{\theta },{\sigma }_{\beta },{\sigma }_{\alpha },M,T)$ in closed form (reported in [44] for space reasons). Interestingly, it can be seen that, in the limit of the product $MT$, it holds that

$$\underset{MT\to \infty }{lim}\Sigma =\left[\begin{array}{ccc}\frac{{\sigma }_{\theta }^2{\sigma }_{\beta }^2}{{\sigma }_{\theta }^2+{\sigma }_{\beta }^2}& -\frac{{\sigma }_{\theta }^2{\sigma }_{\beta }^2}{{\sigma }_{\theta }^2+{\sigma }_{\beta }^2}& 0\\ -\frac{{\sigma }_{\theta }^2{\sigma }_{\beta }^2}{{\sigma }_{\theta }^2+{\sigma }_{\beta }^2}& \frac{{\sigma }_{\theta }^2{\sigma }_{\beta }^2}{{\sigma }_{\theta }^2+{\sigma }_{\beta }^2}& 0\\ 0& 0& 0\end{array}\right].$$

and, in particular, as shown in [44], that

$${[\Sigma ]}_{22}={\sigma }_{22}\left({\sigma }_{\theta }{\sigma }_{\beta }\right),{[\Sigma ]}_{33}={\sigma }_{33}\left(\frac{1}{M{T}^2}\right).$$

(16)

Hence, while for growing M or T, ${\sigma }_{33}\to 0$ meaning that the uncertainty on the skew parameter goes to zero and, consequently, the parameter is perfectly estimated, residual uncertainty remains on both $\theta$ and $\beta$ for which ${\sigma }_{11},{\sigma }_{22}\nrightarrow 0$. As can be seen from the output matrix C, this is because $\theta$ and $\beta$ are linearly dependent. Nonetheless, similarly to ${\sigma }_{22}$, even ${\sigma }_{11}$ and ${\sigma }_{12}$ are functions of the product ${\sigma }_{\theta }{\sigma }_{\beta }$. Thus, ${\sigma }_{11},{\sigma }_{12},{\sigma }_{22}\to 0$ for ${\sigma }_{\theta }{\sigma }_{\beta }\to 0$, meaning that if ${\sigma }_{\beta }=0$ then $\theta$ is perfectly estimated and vice versa. As highlighted later in the simulation section, this suggests that the different performance between the proposed SASE and what will be referred to as Ground Truth (GT) is majorly due to this linear dependence.

6. Simulation and Results

This section presents the test for the proposed SASE algorithm on two different data sets: (i) synthetic data generated from the standard IEEE 123 nodes test-bed [45]; (ii) field data collected from the smart grid located inside the EPFL campus, Switzerland [46]. For power flow computations, the Matlab Matpower package [47] has been used. Finally, if not differently specified,

Table 1. Parameters used in the simulations.

Parameter	Value [Units]	Ref.
T	1 [s]–1 pps gps resync signal	[24,50]
M	{20, 25, 30, 50, 60} [samples]	[8]
${\sigma }_p,{\sigma }_q$	$50\%$	[35]
${\sigma }_{\mathrm{pmu},v},{\sigma }_{\mathrm{pmu},\theta }$	$0.1\%$, ${10}^{-3}$ [rad]	[8]
${\sigma }_{\alpha }$	${10}^{-2}$ [rad]	[48]
${\sigma }_{\beta }$	$2\times {10}^{-4}$ [rad]	[49]
${\mathbf{v}}^{*},{\theta }^{*},{\mathbf{p}}^{*},{\mathbf{q}}^{*}$	power flow nominal solution

summarizes the value for the different considered parameters. Some observations are in order. First, regarding the PMU reporting rate, since it depends on the grid frequency, a set of values for M is considered. Second, given the relatively small sync period $T=1$ [s], synchronization error parameters $\beta ,\alpha$ are assumed constant within the interval $[kT,(k+1\left)T\right)$. Third, for ${\sigma }_{\alpha }$ it has been assumed that PMUs are equipped with quartz-crystal oscillator characterized by a clock stability ≈$10÷30$ ppm [48]. also, for ${\sigma }_{\beta }$ a 50 Hz frequency signal with synchronization ≈$0.5÷$ 1 $\mathsf{\mu }$s accurate has been assumed [49]. This is also in compliance with the accuracy required for the network time distribution in power system [28]. Finally, given the small re-synchronization period of 1[s], reasonable also according to the last research in PMUs synchronization methods [27], no process noise, $W=0$, has been considered.

6.1. Synthetic Data Set—IEEE 123 Nodes Grid

The proposed SASE algorithm has been tested on synthetic data generated from the standard IEEE 123 nodes test-bed [45] whose topology is reported in Figure 3 and compared with:

• an online iterative version of the Bayesian Linear State Estimation algorithm presented in [18] (denoted as BLSE) assuming no synchronization error in the measurements;
• a Ground Truth (denoted as GT) strategy assuming perfect knowledge and compensation of the synchronization error.

The estimation performance is evaluated in terms of Average of Root Mean Square Error, defined as

$$\widehat{\mathrm{ARMSE}}(\mathbf{a},\widehat{\mathbf{a}},t)=\sqrt{\frac{1}{N}\sum _{i=1}^N\frac{1}{n}\sum _{h=1}^n{|a_h^i\left(t\right)-{\widehat{a}}_h^i\left(t\right)|}^2}$$

(17)

where $h=1,\dots ,n$ indicates the h-th node and $i=1,\dots ,N$ the i-th Monte Carlo run, and ($a_h^i,{\widehat{a}}_h^i$) the true and estimated values, respectively. The matrix $\Sigma \left(t\right)$ can be used to compute the theoretical ARMSE as

$$\mathrm{ARMSE}\left(t\right)=\sqrt{\frac{1}{n}\mathrm{Tr}(\Sigma \left(t\right))}.$$

(18)

and we expect $\widehat{\mathrm{ARMSE}}\left(t\right)\approx \mathrm{ARMSE}\left(t\right)$ for large N. After $N=500$ Monte Carlo runs the behaviors of ARMSE (both empirical (17) and theoretical (18)) for the voltage phasor vectors ($\mathbf{u},\widehat{\mathbf{u}}$) and the synchronization error parameters, as functions of the number of PMUs deployed in the grid, are shown in Figure 4a,b, respectively. For PMU deployment the greedy approach presented in [18] has been used. The PMU measurements have been processed with $M=30$, right before a new synchronization instant occurs.

Figure 4a shows that SASE performs quite well with respect to the ground truth (GT) while in the case of BLSE, where the measurement model does not consider any synchronization error, the performance deteriorates and the improvement with respect to the increasing number of PMU is about 30%. On the other hand, in the case of the proposed estimator, the performance is improved of about 60% with only one PMU deployed. Looking at Figure 4b, it is possible to note that the performance of the synchronization error parameters are almost insensitive to the number of PMUs since PMUs are not correlated. Moreover, in both Figure 4a,b the theoretical and empirical curves have same values suggesting the validity of the adopted linear approximation.

This leads to the conclusions that the adopted linearized model could be also effectively used for other optimization problems, e.g., optimal PMU placement or parameter sensitivity analysis also for large power grids with a relevant benefit in terms of computational burden with respect to the Monte Carlo simulations. It is worth highlighting that these results hold for the values in Table 1.

To emphasize the analysis of Section 5, in Figure 5 the performance for a fixed number of PMUs as a function of the number of collected PMU measurements M are shown. Also, in this case, Figure 5a reports the valuable performance of the proposed SASE where the behavior decreases when the collected measurement increase. Conversely, the BLSE results almost constant because does not consider the effect of the delay. These considerations are also confirmed by Figure 5b. Observe that the SASE can compensate the estimated skew, which improves for increasing M.

As stressed in Section 5, this is an intrinsic modeling problem due to the fact that offset and power demand happens to be linearly dependent. In addition, Figure 5b reports the values ${\sigma }_{22}$, ${\sigma }_{33}$ in (16) as a function of M computed for the two-nodes network using the parameters value of Table 1. Observe how the theoretical values corresponding to the two-nodes case turn out to be extremely close to those obtained from the real network. This fact is interesting mainly for two reasons: (i) it supports the claim that in the limit for M (or T), the proposed estimator perfectly reconstructs the skews while residual error remains in the offsets; (ii) from the closed-form expressions for ${\sigma }_{22}$ and ${\sigma }_{33}$, it is possible to retrieve, at least approximately, the value of the parameters needed to obtain a desired level of estimation accuracy.

6.2. Real-World Data Set—EPFL Smartgrid

In this section, we present the result on the test performed using the real data, available for research purpose here [51] and collected from the 20 kV 3-phase 6 nodes smart grid installed in Lausanne within the framework of the NanoTera S${}^3$-Grid project and located inside the EPFL campus [12,46]. Figure 6 shows a graphical representation of the grid topology, where from bus 2 to 5 are connected after the transformers loads and PV systems, while bus 6 has zero-injection. For a complete overview of the network parameters we refer to the original publications [13] and [46].

We recall that the grid is characterized by a line topology with nodes 1 to 5 monitored with PMUs (measuring current and voltage at 50Hz) and node 6 (the last along the line) is not monitored. Also, to estimate the measurements characteristics and noise variance values, we resort to the description reported in Section 4.2.1 of [13] where the variances are computed from the datasheets of the PMU devices. Before analyzing the simulation results and to comment the conclusions coming from the application of SASE to real-monitored data, we would like to stress our choice (already introduced in Section 4.2) regarding the state model during synchronization instants $\tau (k,0)$, $k\in \mathbb{Z}$. Assume the state is $\mathbf{x}={[\delta {\mathbf{p}}^T,\delta {\mathbf{q}}^T,{\mathit{\alpha }}^T,{\mathit{\beta }}^T]}^T$ with evolution

$$\mathbf{x}(k,t+1)=F_{\tau (k,t)}\mathbf{x}(k,t)+\mathbf{w}(k,t).$$

(19)

To address the drift in the measurements, the state matrix is

$$\left\{\begin{array}{cc}\hfill & F_{\tau (k,t)}=I,k\in \mathbb{Z},t\in \{1,\dots ,M-1\};\hfill \\ \hfill & F_{\tau (k,0)}=\left[\begin{array}{cccc}I& 0& 0& 0\\ 0& I& 0& 0\\ 0& 0& I& 0\\ 0& 0& I·T& I\end{array}\right],k\in \mathbb{Z}.\hfill \end{array}\right.$$

(20)

The last row-block of F acts as an integrator for $\beta$ and is used to set their mean values at $\tau (k,0)$ to

$$\beta (k,0)=\beta (k-1,M-1)+\alpha (k-1,M-1)T.$$

(21)

We now turn to the test where a small subset of data consisting of a time window of 6 s collected on 17 November 2014, starting at 10:03:20 A.M. To better characterize the test, we compare, in terms of estimation and prediction, the proposed SASE with the Kalman-based estimator in [12] which use PMU measurements only. We assume to have at disposal measurements from nodes 1, 2 and 3 to perform the estimation while we use nodes 4 and 5 for validation, i.e., we do not use their measurements during estimation. The zero-injection node 6 is considered to be a virtual measurement for the algorithm in [12] and it is eliminated for the SASE.

Figure 7 and Figure 8 show the evolution of estimates and predictions at node 3 and 5, respectively, in the time interval $[2,6]$ s. The first two seconds of simulation have been cut out in order to let the estimator in [12] to properly compute the covariance matrix Q and converge to its steady state. From Figure 7 it is possible to see that the estimator in [12] nicely follows the measurements. The SASE, in harmony with its static state-space model (since we assumed $W=0$), captures the average demand. However, the confidence interval returned by the estimator is in perfect accordance with the measurement values. Conversely, Figure 8 highlights the first difference between the two algorithms. Indeed, due to lack of (prior) information, in the prediction task, the algorithm in [12] does not provide any useful value (not even reported due to the high difference in the scaling factor). This analysis is supported by

Table 2. Total vector error (TVE) for estimation and prediction of SASE and the algorithm [12] in the time interval $[2,6]$ s, using 3 PMUs for estimation.

Alg. [12]—Est.	SASE—Est.	Alg. [12]—Pred.	SASE—Pred.
$2.1\times {10}^{-3}$	$1.4\times {10}^{-2}$	$6\times {10}^5$	$8.5\times {10}^{-3}$

reporting the values of TVE of estimates and prediction regarding the corresponding measurements.

Taking into account the formulation of two algorithms and the results we can make the following comments. The SASE estimator using the state vector based on power demand results more adaptable even in the case of pre-existent infrastructure or in the case of unknown or very noise measurements (like using forecasting or historical data) because use the power demand. On the contrary, the algorithm in [12] is based on the assumption that every node is instrumented with PMUs so the grid is observable. However, this condition is more possible in a future scenario than at present grid. No information regarding the posterior covariance characterizing the estimates is presented in [12]. Conversely, we show how the SASE nicely provides accurate estimates characterized by meaningful confidence intervals. Finally, considering the phase angle rotation due to the deviation between the actual grid frequency and the nominal one, SASE automatically accounts for any linear trend existing in the phase angle measurements. So, in conclusion, even if the SASE algorithm was originally motivated to compensate synchronization errors, it can be effectively used also to compensate linear frequency deviations at different nodes.

7. Conclusions and Future Directions

This paper presents a real-time state estimator based on Kalman filter which uses information from historical/forecasted power measurements and a limited number of PMUs. Despite their high accuracy, also PMUs can suffer from inaccuracy given by internal or external synchronization errors. By testing the proposed algorithm on both synthetic and real-field grids, a comparison with two state-of-the-art SE methods shows how the presence of synchronization errors can easily mislead the estimator if the measurement model does not properly account for it. In particular, in the case of IEEE 123-bus grid the accuracy of the estimated node voltages increases of about 30% using the proposed SASE with respect to the BLSE algorithm and with only one PMU the SASE algorithm improves its performance of about 60%. Moreover, the theoretical and empirical curves are overlapped to indicate the validity of the used linear approximation. The validity of the proposed approach is also evaluated in a small real grid, where SASE is compared with the algorithm proposed in [12] where the benefit of SASE in total vector error (TVE) computation differ with one order of magnitude in estimation stage. As future directions we foresee the extension of this methodology to the detection and mitigation of GPS spoofing attacks and time-synchronization attacks as well as more dynamic scenarios by considering fast changing loads.