PMU State Estimation

To initiate the process, let us construct the PowerSystem type and formulate the AC model:

system = powerSystem()

addBus!(system; label = 1, type = 3, active = 0.5)
addBus!(system; label = 2, type = 1, reactive = 0.3)
addBus!(system; label = 3, type = 1, active = 0.5)

@branch(resistance = 0.02, susceptance = 0.04)
addBranch!(system; label = 1, from = 1, to = 2, reactance = 0.6)
addBranch!(system; label = 2, from = 1, to = 3, reactance = 0.7)
addBranch!(system; label = 3, from = 2, to = 3, reactance = 0.2)

addGenerator!(system; label = 1, bus = 1, active = 3.2, reactive = 0.2)

acModel!(system)

To review, we can conceptualize the bus/branch model as the graph denoted by $\mathcal{G} = (\mathcal{N}, \mathcal{E})$, where we have the set of buses $\mathcal{N} = \{1, \dots, n\}$, and the set of branches $\mathcal{E} \subseteq \mathcal{N} \times \mathcal{N}$ within the power system:

julia> 𝒩 = collect(keys(system.bus.label))3-element Vector{String}:
 "1"
 "2"
 "3"
julia> ℰ = [𝒩[system.branch.layout.from] 𝒩[system.branch.layout.to]]3×2 Matrix{String}:
 "1"  "2"
 "1"  "3"
 "2"  "3"

Following that, we will introduce the Measurement type and incorporate a set of PMUs $\mathcal{M} \equiv \bar{\mathcal{P}}$ into the graph $\mathcal{G}$, that capture both bus voltage and branch current phasors. To construct the linear PMU state estimation model, we represent the vector of state variables, as well as phasor measurements, in the rectangular coordinate system. Thus, we initialize the Measurement type:

device = measurement()

Notation

Here, when referring to a vector $\mathbf a$, we use the notation $\mathbf a = [a_i]$ or $\mathbf a = [a_{ij}]$, where $a_i$ represents the element related with bus $i \in \mathcal N$ or measurement $i \in \mathcal M$, while $a_{ij}$ denotes the element related with branch $(i,j) \in \mathcal E$.

State Estimation Model

Initially, PMUs output phasor measurements in polar coordinates. However, these measurements can be interpreted in rectangular coordinates, where the real and imaginary parts of bus voltages and branch current phasors serve as measurements. Additionally, to obtain the linear system of equations, we observe a vector of state variables in rectangular coordinates $\mathbf x \equiv[\mathbf V_\mathrm{re}, \mathbf V_\mathrm{im}]$:

$\mathbf V_\mathrm{re} =\big[\Re(\bar{V}_1),\dots,\Re(\bar{V}_n)\big]^T$, representing the real parts of complex bus voltages,
$\mathbf V_\mathrm{im} =\big[\Im(\bar{V}_1),\dots,\Im(\bar{V}_n)\big]^T$, representing the imaginary parts of complex bus voltages.

Consequently, the total number of state variables is $s = 2n$. It is worth noting that in this approach to state estimation, we do not require the slack bus.

The primary drawback of this method stems from measurement errors, which are associated with polar coordinates. Consequently, the covariance matrix must be transformed from polar to rectangular coordinates. As a result, errors from a single PMU are correlated, leading to a non-diagonal covariance matrix. Despite this, the covariance matrix is commonly treated as diagonal, impacting the state estimation accuracy in such scenarios.

Hence, the model includes real and imaginary parts of bus voltage and current phasor measurements from the set $\mathcal M$, contributing to the formulation of a linear system of equations:

\[ \mathbf z=\mathbf h(\mathbf x) + \mathbf u.\]

Here, $\mathbf h(\mathbf x)= [h_1(\mathbf x)$, $\dots$, $h_k(\mathbf x)]^T$ represents the vector of linear measurement functions, $\mathbf z = [z_1, \dots, z_k]^T$ denotes the vector of measurement values, and $\mathbf u = [u_1, \dots, u_k]^T$ represents the vector of measurement errors, where $k = 2|\bar{\mathcal P}|$.

These errors are assumed to follow a Gaussian distribution with a zero mean and covariance matrix $\bm \Sigma$. The diagonal elements of $\bm \Sigma$ correspond to the measurement variances $\mathbf v = [v_1, \dots, v_k]^T$, while the off-diagonal elements represent the covariances between the measurement errors $\mathbf w = [w_1, \dots, w_k]^T$.

In summary, upon defining the PMU, each $i$-th PMU is associated with two measurement functions $h_{2i-1}(\mathbf x)$, $h_{2i}(\mathbf x)$, along with their respective measurement values $z_{2i-1}$, $z_{2i}$, as well as their variances $v_{2i-1}$, $v_{2i}$, and possibly covariances $w_{2i-1}$, $w_{2i}$.

Bus Voltage Phasor Measurements

When a PMU $(V_i, \theta_i) \in \bar{\mathcal P}$ is introduced at bus $i \in \mathcal N$ in this type of state estimation, users specify the measurement values, variances, and measurement functions of vectors as follows:

\[ \mathbf z = [z_{\Re(\bar{V}_i)}, z_{\Im(\bar{V}_i)}], \;\;\; \mathbf v = [v_{\Re(\bar{V}_i)}, v_{\Im(\bar{V}_i)}], \;\;\; \mathbf h(\mathbf x) = [h_{\Re(\bar{V}_i)}(\mathbf x), h_{\Im(\bar{V}_i)}(\mathbf x)].\]

For example:

addPmu!(system, device; label = "V₂, θ₂", bus = 2, magnitude = 0.9, angle = -0.1,
varianceMagnitude = 1e-5, varianceAngle = 1e-5)

Here, measurement values are obtained according to:

\[ \begin{aligned} z_{\Re(\bar{V}_i)} = z_{V_i} \cos z_{\theta_i}\\ z_{\Im(\bar{V}_i)} = z_{V_i} \sin z_{\theta_i}. \end{aligned}\]

Utilizing the classical theory of propagation of uncertainty [21], the variances can be calculated as follows:

\[ \begin{aligned} v_{\Re(\bar{V}_i)} &= v_{V_i} \left[ \cfrac{\mathrm \partial} {\mathrm \partial z_{V_i}} (z_{V_i} \cos z_{\theta_i}) \right]^2 + v_{\theta_i} \left[ \cfrac{\mathrm \partial} {\mathrm \partial z_{\theta_i}} (z_{V_i} \cos z_{\theta_i})\right]^2 = v_{V_i} (\cos z_{\theta_i})^2 + v_{\theta_i} (z_{V_i} \sin z_{\theta_i})^2\\ v_{\Im(\bar{V}_i)} &= v_{V_i} \left[ \cfrac{\mathrm \partial} {\mathrm \partial z_{V_i}} (z_{V_i} \sin z_{\theta_i}) \right]^2 + v_{\theta_i} \left[ \cfrac{\mathrm \partial} {\mathrm \partial z_{\theta_i}} (z_{V_i} \sin z_{\theta_i})\right]^2 = v_{V_i} (\sin z_{\theta_i})^2 + v_{\theta_i} (z_{V_i} \cos z_{\theta_i})^2. \end{aligned}\]

Lastly, the functions defining the bus voltage phasor measurement are:

\[ \begin{aligned} h_{\Re(\bar{V}_i)}(\mathbf x) &= \Re(\bar{V}_i)\\ h_{\Im(\bar{V}_i)}(\mathbf x) &= \Im(\bar{V}_i). \end{aligned}\]

The coefficient expressions for measurement functions are as follows:

\[ \cfrac{\mathrm \partial{h_{\Re(\bar{V}_i)}(\mathbf x)}}{\mathrm \partial \Re(\bar{V}_i)} = 1, \;\;\; \cfrac{\mathrm \partial{h_{\Im(\bar{V}_i)}(\mathbf x)}}{\mathrm \partial \Im(\bar{V}_i)} = 1.\;\;\;\]

In the previous example, the user neglected the covariances between the real and imaginary parts of the measurement. However, if desired, the user can also include them in the state estimation model by specifying the covariances of the vector:

\[ \mathbf w = [w_{\Re(\bar{V}_i)}, w_{\Im(\bar{V}_i)}].\]

addPmu!(system, device; label = "V₃, θ₃", bus = 3, magnitude = 0.9, angle = -0.2,
varianceMagnitude = 1e-5, varianceAngle = 1e-5, correlated = true)

Then, the covariances are obtained as follows:

\[ w_{\Re(\bar{V}_{i})} = w_{\Im(\bar{V}_{i})} = v_{V_i} \cfrac{\mathrm \partial} {\mathrm \partial z_{V_i}} (z_{V_i} \cos z_{\theta_i}) \cfrac{\mathrm \partial} {\mathrm \partial z_{V_i}} (z_{V_i} \sin z_{\theta_i}) + v_{\theta_i} \cfrac{\mathrm \partial} {\mathrm \partial z_{\theta_i}} (z_{V_i} \cos z_{\theta_i}) \cfrac{\mathrm \partial} {\mathrm \partial z_{\theta_i}} (z_{V_i} \sin z_{\theta_i}),\]

which results in the solution:

\[ w_{\Re(\bar{V}_{i})} = w_{\Im(\bar{V}_{i})} = \cos z_{\theta_i} \sin z_{\theta_i}(v_{V_i} - v_{\theta_i} z_{V_i}^2).\]

From-Bus Current Phasor Measurements

If the user chooses to include phasor measurement $(I_{ij}, \psi_{ij}) \in \bar{\mathcal P}$ in the state estimation model, the user will specify the measurement values, variances, and measurement functions of vectors:

\[ \mathbf z = [z_{\Re(\bar{I}_{ij})}, z_{\Im(\bar{I}_{ij})}], \;\;\; \mathbf v = [v_{\Re(\bar{I}_{ij})}, v_{\Im(\bar{I}_{ij})}], \;\;\; \mathbf h(\mathbf x) = [h_{\Re(\bar{I}_{ij})}(\mathbf x), h_{\Im(\bar{I}_{ij})}(\mathbf x)].\]

For example:

addPmu!(system, device; label = "I₂₃, ψ₂₃", from = 3, magnitude = 0.3, angle = 0.4,
varianceMagnitude = 1e-3, varianceAngle = 1e-4)

Here, measurement values are obtained according to:

\[ \begin{aligned} z_{\Re(\bar{I}_{ij})} = z_{I_{ij}} \cos z_{\psi_{ij}}\\ z_{\Im(\bar{I}_{ij})} = z_{I_{ij}} \sin z_{\psi_{ij}}. \end{aligned}\]

Utilizing the classical theory of propagation of uncertainty [21], the variances can be calculated as follows:

\[ \begin{aligned} v_{\Re(\bar{I}_{ij})} & = v_{I_{ij}} (\cos z_{\psi_{ij}})^2 + v_{\psi_{ij}} (z_{I_{ij}} \sin z_{\psi_{ij}})^2 \\ v_{\Im(\bar{I}_{ij})} &= v_{I_{ij}} (\sin z_{\psi_{ij}})^2 + v_{\psi_{ij}} (z_{I_{ij}} \cos z_{\psi_{ij}})^2. \end{aligned}\]

The functions defining the current phasor measurement at the from-bus end are:

\[ \begin{aligned} h_{\Re(\bar{I}_{ij})}(\mathbf x) &= A\Re(\bar{V}_i) - B\Im(\bar{V}_i) - \left(C\cos\phi_{ij} - D\sin\phi_{ij}\right) \Re(\bar{V}_j) + \left(C\sin\phi_{ij} + D\cos\phi_{ij} \right) \Im(\bar{V}_j) \\ h_{\Im(\bar{I}_{ij})}(\mathbf x) &= B\Re(\bar{V}_i) + A\Im(\bar{V}_i) - \left(C\sin \phi_{ij} + D\cos\phi_{ij}\right) \Re(\bar{V}_j) - \left(C\cos\phi_{ij} - D\sin\phi_{ij} \right)\Im(\bar{V}_j), \end{aligned}\]

where:

\[ A = \cfrac{g_{ij} + g_{\mathrm{s}ij}}{\tau_{ij}^2},\;\;\; B = \cfrac{b_{ij}+b_{\mathrm{s}ij}} {\tau_{ij}^2},\;\;\; C = \cfrac{g_{ij}}{\tau_{ij}},\;\;\; D = \cfrac{b_{ij}}{\tau_{ij}}.\]

The coefficient expressions for measurement functions are as follows:

\[ \begin{aligned} \cfrac{\mathrm \partial{h_{\Re(\bar{I}_{ij})}(\mathbf x)}}{\mathrm \partial \Re(\bar{V}_i)} &= \cfrac{\mathrm \partial{h_{\Im(\bar{I}_{ij})}(\mathbf x)}}{\mathrm \partial \Im(\bar{V}_i)} = A \\ \cfrac{\mathrm \partial{h_{\Re(\bar{I}_{ij})}(\mathbf x)}} {\mathrm \partial \Re(\bar{V}_j)} &= \cfrac{\mathrm \partial{h_{\Im(\bar{I}_{ij})}(\mathbf x)}} {\mathrm \partial \Im(\bar{V}_j)} = - \left(C\cos\phi_{ij} - D\sin\phi_{ij}\right)\\ \cfrac{\mathrm \partial{h_{\Re(\bar{I}_{ij})}(\mathbf x)}}{\mathrm \partial \Im(\bar{V}_i)} &=- \cfrac{\mathrm \partial{h_{\Im(\bar{I}_{ij})}(\mathbf x)}}{\mathrm \partial \Re(\bar{V}_i)} = -B \\ \cfrac{\mathrm \partial{h_{\Re(\bar{I}_{ij})}(\mathbf x)}}{\mathrm \partial \Im(\bar{V}_j)} &= - \cfrac{\mathrm \partial{h_{\Im(\bar{I}_{ij})}(\mathbf x)}}{\mathrm \partial\Re(\bar{V}_j)} = \left(C\sin\phi_{ij} + D\cos\phi_{ij} \right). \end{aligned}\]

In the previous example, the user neglects the covariances between the real and imaginary parts of the measurement. However, if desired, the user can also include them in the state estimation model by specifying the covariances of the vector:

\[ \mathbf w = [w_{\Re(\bar{I}_{ij})}, w_{\Im(\bar{I}_{ij})}].\]

addPmu!(system, device; label = "I₁₃, ψ₁₃", from = 2, magnitude = 0.3, angle = -0.5,
varianceMagnitude = 1e-5, varianceAngle = 1e-5, correlated = true)

Then, the covariances are obtained as follows:

\[ w_{\Re(\bar{I}_{ij})} = w_{\Im(\bar{I}_{ij})} = \sin z_{\psi_{ij}} \cos z_{\psi_{ij}}(v_{I_{ij}} - v_{\psi_{ij}} z_{I_{ij}}^2).\]

To-Bus Current Phasor Measurements

If the user chooses to include phasor measurement $(I_{ji}, \psi_{ji}) \in \bar{\mathcal P}$ in the state estimation model, the user will specify the measurement values, variances, and measurement functions of vectors:

\[ \mathbf z = [z_{\Re(\bar{I}_{ji})}, z_{\Im(\bar{I}_{ji})}], \;\;\; \mathbf v = [v_{\Re(\bar{I}_{ji})}, v_{\Im(\bar{I}_{ji})}], \;\;\; \mathbf h(\mathbf x) = [h_{\Re(\bar{I}_{ji})}(\mathbf x), h_{\Im(\bar{I}_{ji})}(\mathbf x)].\]

For example:

addPmu!(system, device; label = "I₃₂, ψ₃₂", to = 3, magnitude = 0.3, angle = -2.9,
varianceMagnitude = 1e-5, varianceAngle = 1e-5)

Here, measurement values are obtained according to:

\[ \begin{aligned} z_{\Re(\bar{I}_{ji})} = z_{I_{ji}} \cos z_{\psi_{ji}}\\ z_{\Im(\bar{I}_{ji})} = z_{I_{ji}} \sin z_{\psi_{ji}}. \end{aligned}\]

The variances can be calculated as follows:

\[ \begin{aligned} v_{\Re(\bar{I}_{ji})} &= v_{I_{ji}} (\cos z_{\psi_{ji}})^2 + v_{\psi_{ji}} (z_{I_{ji}} \sin z_{\psi_{ji}})^2 \\ v_{\Im(\bar{I}_{ji})} &= v_{I_{ji}} (\sin z_{\psi_{ji}})^2 + v_{\psi_{ji}} (z_{I_{ji}} \cos z_{\psi_{ji}})^2. \end{aligned}\]

The functions defining the current phasor measurement at the to-bus end are:

\[ \begin{aligned} h_{\Re(\bar{I}_{ji})}(\mathbf x) &= \tau_{ij}^2 A \Re(\bar{V}_j) - \tau_{ij}^2 B \Im(\bar{V}_j) - \left(C \cos\phi_{ij} + D \sin \phi_{ij}\right) \Re(\bar{V}_i) - \left( C\sin \phi_{ij} - D\cos \phi_{ij} \right) \Im(\bar{V}_i)\\ h_{\Im(\bar{I}_{ji})}(\mathbf x) &= \tau_{ij}^2 B \Re(\bar{V}_j) + \tau_{ij}^2 A \Im(\bar{V}_j) + \left(C \sin \phi_{ij} - D \cos\phi_{ij} \right) \Re(\bar{V}_i) - \left(C\cos \phi_{ij} + D\sin \phi_{ij}\right) \Im(\bar{V}_i). \end{aligned}\]

The coefficient expressions for measurement functions are as follows:

\[ \begin{aligned} \cfrac{\mathrm \partial{h_{\Re(\bar{I}_{ji})}(\mathbf x)}}{\mathrm \partial \Re(\bar{V}_i)} &= \cfrac{\mathrm \partial{h_{\Im(\bar{I}_{ji})}(\mathbf x)}}{\mathrm \partial \Im(\bar{V}_i)} = - \left(C \cos\phi_{ij} + D \sin \phi_{ij}\right)\\ \cfrac{\mathrm \partial{h_{\Re(\bar{I}_{ji})}(\mathbf x)}} {\mathrm \partial \Re(\bar{V}_j)} &= \cfrac{\mathrm \partial{h_{\Im(\bar{I}_{ji})}(\mathbf x)}} {\mathrm \partial \Im(\bar{V}_j)} = \tau_{ij}^2A\\ \cfrac{\mathrm \partial{h_{\Re(\bar{I}_{ji})}(\mathbf x)}}{\mathrm \partial \Im(\bar{V}_i)} &= - \cfrac{\mathrm \partial{h_{\Im(\bar{I}_{ji})}(\mathbf x)}}{\mathrm \partial \Re(\bar{V}_i)} = -\left(C\sin \phi_{ij} - D\cos \phi_{ij} \right) \\ \cfrac{\mathrm \partial{h_{\Re(\bar{I}_{ji})}(\mathbf x)}}{\mathrm \partial \Im(\bar{V}_j)} &= - \cfrac{\mathrm \partial{h_{\Im(\bar{I}_{ji})}(\mathbf x)}}{\mathrm \partial \Re(\bar{V}_j)} = -\tau_{ij}^2B. \end{aligned}\]

As before, we are neglecting the covariances between the real and imaginary parts of the measurement. If desired, we can include them in the state estimation model by specifying the covariances of the vector:

\[ \mathbf w = [w_{\Re(\bar{I}_{ji})}, w_{\Im(\bar{I}_{ji})}].\]

addPmu!(system, device; label = "I₃₁, ψ₃₁", to = 2, magnitude = 0.3, angle = 2.5,
varianceMagnitude = 1e-5, varianceAngle = 1e-5, correlated = true)

Then, the covariances are obtained as follows:

\[ w_{\Re(\bar{I}_{ji})} = w_{\Im(\bar{I}_{ji})} = \sin z_{\psi_{ji}} \cos z_{\psi_{ji}}(v_{I_{ji}} - v_{\psi_{ji}} z_{I_{ji}}^2).\]

Weighted Least-Squares Estimation

The solution to the PMU state estimation problem is determined by solving the linear weighted least-squares (WLS) problem, represented by the following formula:

\[ \mathbf H^T \bm \Sigma^{-1} \mathbf H \mathbf x = \mathbf H^T \bm \Sigma^{-1} \mathbf z.\]

Here, $\mathbf z \in \mathbb{R}^k$ denotes the vector of measurement values, $\mathbf H \in \mathbb {R}^{k \times 2n}$ represents the coefficient matrix, and $\bm \Sigma \in \mathbb {R}^{k \times k}$ is the measurement error covariance matrix.

Implementation

JuliaGrid initiates the PMU state estimation framework by setting up the WLS model, as illustrated in the following:

analysis = pmuStateEstimation(system, device)

Coefficient Matrix

Using the above-described equations, JuliaGrid forms the coefficient matrix $\mathbf H$:

julia> 𝐇 = analysis.method.coefficient12×6 SparseArrays.SparseMatrixCSC{Float64, Int64} with 36 stored entries:
   ⋅         1.0        ⋅          ⋅          ⋅         ⋅
   ⋅          ⋅         ⋅          ⋅         1.0        ⋅
   ⋅          ⋅        1.0         ⋅          ⋅         ⋅
   ⋅          ⋅         ⋅          ⋅          ⋅        1.0
   ⋅         0.49505  -0.49505     ⋅         4.9305   -4.9505
   ⋅        -4.9305    4.9505      ⋅         0.49505  -0.49505
  0.040783    ⋅       -0.040783   1.40741     ⋅       -1.42741
 -1.40741     ⋅        1.42741    0.040783    ⋅       -0.040783
   ⋅        -0.49505   0.49505     ⋅        -4.9505    4.9305
   ⋅         4.9505   -4.9305      ⋅        -0.49505   0.49505
 -0.040783    ⋅        0.040783  -1.42741     ⋅        1.40741
  1.42741     ⋅       -1.40741   -0.040783    ⋅        0.040783

In this matrix, each row corresponds to a specific measurement in the rectangular coordinate system. Therefore, the $i$-th PMU is associated with the $2i - 1$ index of the row, representing the real part of the phasor measurement, while the $2i$ row corresponds to the imaginary part of the phasor measurement. Columns are ordered based on how the state variables are defined $\mathbf x \equiv[\mathbf{V}_\mathrm{re},\mathbf{V}_\mathrm{im}]$.

Precision Matrix

JuliaGrid opts not to retain the covariance matrix $\bm \Sigma$ but rather stores its inverse, the precision or weighting matrix denoted as $\mathbf W = \bm \Sigma^{-1}$. The order of these values corresponds to the description provided for the coefficient matrix. Users can access these values using the following command:

julia> 𝐖 = analysis.method.precision12×12 SparseArrays.SparseMatrixCSC{Float64, Int64} with 18 stored entries:
 1.0019e5   ⋅             ⋅         …   ⋅          ⋅          ⋅
  ⋅        1.23169e5      ⋅             ⋅          ⋅          ⋅
  ⋅         ⋅            1.00926e5      ⋅          ⋅          ⋅
  ⋅         ⋅         4567.25           ⋅          ⋅          ⋅
  ⋅         ⋅             ⋅             ⋅          ⋅          ⋅
  ⋅         ⋅             ⋅         …   ⋅          ⋅          ⋅
  ⋅         ⋅             ⋅             ⋅          ⋅          ⋅
  ⋅         ⋅             ⋅             ⋅          ⋅          ⋅
  ⋅         ⋅             ⋅             ⋅          ⋅          ⋅
  ⋅         ⋅             ⋅            7.03786e5   ⋅          ⋅
  ⋅         ⋅             ⋅         …   ⋅         4.62149e5  4.84789e5
  ⋅         ⋅             ⋅             ⋅         4.84789e5  7.48963e5

The precision matrix does not maintain a diagonal form, indicating that correlations between the real and imaginary parts of the phasor measurements are included in the model. To ignore these correlations, simply omit the correlated keyword within the function that adds a PMU. For example:

device = measurement()

@pmu(label = "PMU ?", noise = false)
addPmu!(system, device; bus = 1, magnitude = 1.0, angle = 0.0)
addPmu!(system, device; bus = 2, magnitude = 0.8745, angle = -0.1529)
addPmu!(system, device; from = 1, magnitude = 0.3033, angle = -0.7136)
addPmu!(system, device; from = 2, magnitude = 0.3142, angle = -0.4950)
addPmu!(system, device; to = 3, magnitude = 0.2809, angle = -2.8954)

Following this, we recreate the WLS state estimation model:

analysis = pmuStateEstimation(system, device)

Upon inspection, it becomes evident that the precision matrix maintains a diagonal structure:

julia> 𝐖 = analysis.method.precision10×10 SparseArrays.SparseMatrixCSC{Float64, Int64} with 10 stored entries:
 1.0e8   ⋅      ⋅          ⋅         …   ⋅          ⋅          ⋅
  ⋅     1.0e8   ⋅          ⋅             ⋅          ⋅          ⋅
  ⋅      ⋅     1.00549e8   ⋅             ⋅          ⋅          ⋅
  ⋅      ⋅      ⋅         1.29835e8      ⋅          ⋅          ⋅
  ⋅      ⋅      ⋅          ⋅             ⋅          ⋅          ⋅
  ⋅      ⋅      ⋅          ⋅         …   ⋅          ⋅          ⋅
  ⋅      ⋅      ⋅          ⋅             ⋅          ⋅          ⋅
  ⋅      ⋅      ⋅          ⋅            3.31017e8   ⋅          ⋅
  ⋅      ⋅      ⋅          ⋅             ⋅         1.05788e8   ⋅
  ⋅      ⋅      ⋅          ⋅             ⋅          ⋅         7.48423e8

Mean Vector

To retrieve the vector $\mathbf z$, containing the means of Gaussian distributions for each measurement, users can utilize:

julia> 𝐳 = analysis.method.mean10-element Vector{Float64}:
  1.0
  0.0
  0.8642976896322709
 -0.1331906667012325
  0.22929794379665394
 -0.19852794002514235
  0.2764861686305354
 -0.14925494482933274
 -0.272430120445671
 -0.06845903500603255

These values represent measurement values in the rectangular coordinate system as described earlier.

Estimate of State Variables

Next, the WLS equation is solved to obtain the estimate of state variables:

\[ \hat{\mathbf x} = [\mathbf H^T \bm \Sigma^{-1} \mathbf H]^{-1} \mathbf H^T \bm \Sigma^{-1} \mathbf z.\]

This process is executed using the solve! function:

solve!(system, analysis)

The initial step involves the LU factorization of the gain matrix:

\[ \mathbf G = \mathbf H^T \bm \Sigma^{-1} \mathbf H = \mathbf L \mathbf U.\]

Tip

By default, JuliaGrid utilizes LU factorization as the primary method to factorize the gain matrix. However, users maintain the flexibility to opt for alternative factorization methods such as LDLt or QR.

Access to the factorized gain matrix is available through:

julia> 𝐋 = analysis.method.factorization.L6×6 SparseArrays.SparseMatrixCSC{Float64, Int64} with 21 stored entries:
  1.0            ⋅            ⋅          ⋅           ⋅         ⋅
 -0.0271958     1.0           ⋅          ⋅           ⋅         ⋅
 -0.0317928    -0.989721     1.0         ⋅           ⋅         ⋅
 -0.0203296    -0.00368218   0.213744   1.0          ⋅         ⋅
  0.000776847  -0.354335     0.226762  -0.0889487   1.0        ⋅
  0.0024918     0.375511    -0.275104  -0.0528311  -0.987063  1.0
julia> 𝐔 = analysis.method.factorization.U6×6 SparseArrays.SparseMatrixCSC{Float64, Int64} with 21 stored entries:
 0.510148  -0.2204    -0.257315    -0.00612351   0.00149083    0.0045228
  ⋅         0.457753  -0.452449    -6.26466e-5  -0.0384085     0.038498
  ⋅          ⋅         0.00471087   3.7474e-5    0.000253295  -0.000290642
  ⋅          ⋅          ⋅           0.519168    -0.294218     -0.165281
  ⋅          ⋅          ⋅            ⋅           0.30568      -0.285374
  ⋅          ⋅          ⋅            ⋅            ⋅            0.0224408

Finally, JuliaGrid obtains the solution in the rectangular coordinate system and then transforms these solutions into the standard form given in the polar coordinate system.

The estimated bus voltage magnitudes $\hat{\mathbf V} = [\hat{V}_i]$ and angles $\hat{\bm {\Theta}} = [\hat{\theta}_i]$, $i \in \mathcal N$, can be retrieved using the variables:

julia> 𝐕 = analysis.voltage.magnitude3-element Vector{Float64}:
 0.9999944798898596
 0.8745066044755129
 0.8963779829803771
julia> 𝚯 = analysis.voltage.angle3-element Vector{Float64}:
  5.93448287474239e-6
 -0.15290401874549062
 -0.2137671552173047

Info

It is essential to note that the slack bus does not exist in the case of the PMU state estimation model.

Alternative Formulation

The resolution of the WLS state estimation problem using the conventional method typically progresses smoothly. However, it is widely acknowledged that in certain situations common to real-world systems, this method can be vulnerable to numerical instabilities. Such conditions might impede the algorithm from converging to a satisfactory solution. In such cases, users may opt for an alternative formulation of the WLS state estimation, namely, employing an approach called orthogonal factorization [5, Sec. 3.2]. This approach is suitable when measurement errors are uncorrelated, and the precision matrix remains diagonal.

To address ill-conditioned situations arising from significant differences in measurement variances, users can employ an alternative approach:

analysis = pmuStateEstimation(system, device, Orthogonal)

To explain the method, we begin with the WLS equation:

\[ \mathbf H^T \mathbf W \mathbf H \hat{\mathbf x} = \mathbf H^T \mathbf W \mathbf z,\]

where $\mathbf W = \bm \Sigma^{-1}$. Subsequently, we can write:

\[ \left({\mathbf W^{1/2}} \mathbf H\right)^T {\mathbf W^{1/2}} \mathbf H \hat{\mathbf x} = \left({\mathbf W^{1/2}} \mathbf H\right)^T {\mathbf W^{1/2}} \mathbf z.\]

Consequently, we have:

\[ \bar{\mathbf H}^T \bar{\mathbf H} \hat{\mathbf x} = \bar{\mathbf H}^T \bar{\mathbf z},\]

where:

\[ \bar{\mathbf H} = {\mathbf W^{1/2}} \mathbf H; \;\;\; \bar{\mathbf z} = {\mathbf W^{1/2}} \mathbf z.\]

At this point, QR factorization is performed on the rectangular matrix:

\[ \bar{\mathbf H} = {\mathbf W^{1/2}} \mathbf H = \mathbf Q \mathbf R.\]

Executing this procedure involves the solve! function:

solve!(system, analysis)

Access to the factorized matrix is possible through:

julia> 𝐐 = analysis.method.factorization.Q10×10 SparseArrays.SPQR.QRSparseQ{Float64, Int64}
julia> 𝐑 = analysis.method.factorization.R6×6 SparseArrays.SparseMatrixCSC{Float64, Int64} with 21 stored entries:
 -36315.5  15689.5           435.909     -106.126   18317.3         -321.961
       ⋅       1.37109e5     -18.7643  -11504.4        -1.35521e5  11531.2
       ⋅        ⋅         -28150.7      15950.1       -54.5971      8965.19
       ⋅        ⋅               ⋅      -54528.5      -206.486      50907.3
       ⋅        ⋅               ⋅            ⋅     -13898.4           66.1322
       ⋅        ⋅               ⋅            ⋅           ⋅         14366.9

To obtain the solution, JuliaGrid avoids materializing the orthogonal matrix $\mathbf Q$ and proceeds to solve the system, resulting in the estimate of bus voltage magnitudes $\hat{\mathbf V} = [\hat{V}_i]$ and angles $\hat{\bm \Theta} = [\hat{\theta}_i]$, where $i \in \mathcal N$:

julia> 𝐕 = analysis.voltage.magnitude3-element Vector{Float64}:
 0.9999944798898643
 0.8745066044755183
 0.8963779829803826
julia> 𝚯 = analysis.voltage.angle3-element Vector{Float64}:
  5.934482874810098e-6
 -0.15290401874549017
 -0.21376715521730386

Bad Data Processing

Besides the state estimation algorithm, one of the essential state estimation routines is the bad data processing, whose main task is to detect and identify measurement errors, and eliminate them if possible. This is usually done by processing the measurement residuals [5, Ch. 5], and typically, the largest normalized residual test is used to identify bad data. The largest normalized residual test is performed after we obtained the solution of the state estimation in the repetitive process of identifying and eliminating bad data measurements one after another [24].

To illustrate this process, let us introduce a new measurement that contains an obvious outlier:

addPmu!(system, device; bus = 3, magnitude = 2.5, angle = 0.1)

Subsequently, we will construct the WLS state estimation model and solve it:

analysis = pmuStateEstimation(system, device)
solve!(system, analysis)

Now, the bad data processing can be executed:

outlier = residualTest!(system, device, analysis; threshold = 4.0)

In this step, we employ the largest normalized residual test, guided by the analysis outlined in [5, Sec. 5.7]. To be more precise, we compute all measurement residuals in the rectangular coordinate system based on the obtained estimate of state variables:

\[ r_{i} = z_i - h_i(\hat {\mathbf x}), \;\;\; i \in \mathcal{M}.\]

The normalized residuals for all measurements are computed as follows:

\[ \bar{r}_{i} = \cfrac{|r_i|}{\sqrt{C_{ii}}} = \cfrac{|r_i|}{\sqrt{S_{ii}\Sigma_{ii}}}, \;\;\; i \in \mathcal{M}.\]

In this equation, we denote the diagonal entries of the residual covariance matrix $\mathbf C \in \mathbb{R}^{k \times k}$ as $C_{ii} = S_{ii}\Sigma_{ii}$, where $S_{ii}$ is the diagonal entry of the residual sensitivity matrix $\mathbf S$ representing the sensitivity of the measurement residuals to the measurement errors. For this specific configuration, the relationship is expressed as:

\[ \mathbf C = \mathbf S \bm \Sigma = \bm \Sigma - \mathbf H [\mathbf H^T \bm \Sigma^{-1} \mathbf H]^{-1} \mathbf H^T.\]

It is important to note that only the diagonal entries of $\mathbf C$ are required. To obtain the inverse, the JuliaGrid package utilizes a computationally efficient sparse inverse method, retrieving only the necessary elements of the inverse.

The subsequent step involves selecting the largest normalized residual, and the $j$-th measurement is then suspected as bad data and potentially removed from the measurement set $\mathcal{M}$:

\[ \bar{r}_{j} = \text{max} \{\bar{r}_{i}, \forall i \in \mathcal{M} \}.\]

Users can access this information using the variable:

julia> outlier.maxNormalizedResidual12860.120552736831

If the largest normalized residual, denoted as $\bar{r}_j$, satisfies the inequality:

\[ \bar{r}_j \ge \epsilon,\]

the corresponding measurement is identified as bad data and subsequently removed. In this example, the bad data identification threshold is set to $\epsilon = 4$. Users can verify the satisfaction of this inequality by inspecting:

julia> outlier.detecttrue

This indicates that the measurement labeled as:

julia> outlier.label"PMU 6"

is removed from the PMU model and marked as out-of-service. Specifically, either the real or imaginary part of the corresponding measurement is identified as the outlier. Consequently, both parts of the measurement are removed from the PMU state estimation model.

Subsequently, we can immediately solve the system again, but this time without the removed measurement:

solve!(system, analysis)

Following that, we check for outliers once more:

outlier = residualTest!(system, device, analysis; threshold = 4.0)

To examine the value:

julia> outlier.maxNormalizedResidual0.40906113630721075

As this value is now less than the threshold $\epsilon = 4$, the measurement is not removed, or there are no outliers. This can also be verified by observing the bad data flag:

julia> outlier.detectfalse

Least Absolute Value Estimation

The least absolute value (LAV) method provides an alternative estimation approach that is considered more robust in comparison to the WLS method. The WLS state estimation problem relies on specific assumptions about measurement errors, whereas robust estimators aim to remain unbiased even in the presence of various types of measurement errors and outliers. This characteristic eliminates the need for bad data processing, as discussed in [5, Ch. 6]. It is important to note that robustness often comes at the cost of increased computational complexity.

It can be demonstrated that the problem can be expressed as a linear programming problem. This section outlines the method as described in [5, Sec. 6.5]. To revisit, we consider the system of linear equations:

\[ \mathbf{z}=\mathbf{h}(\mathbf x)+\mathbf{u}+\mathbf{w}.\]

Subsequently, the LAV state estimator is derived as the solution to the optimization problem:

\[ \begin{aligned} \text{minimize}& \;\;\; \mathbf a^T |\mathbf r|\\ \text{subject\;to}& \;\;\; \mathbf{z} - \mathbf{H} \mathbf x =\mathbf r. \end{aligned}\]

Here, $\mathbf a \in \mathbb {R}^{k}$ is the vector with all entries equal to one, and $\mathbf r$ represents the vector of measurement residuals. Let $\bm \eta$ be defined in a manner that ensures:

\[ |\mathbf r| \preceq \bm \eta,\]

and replace the above inequality with two equalities using the introduction of two non-negative slack variables $\mathbf q \in \mathbb {R}_{\ge 0}^{k}$ and $\mathbf w \in \mathbb {R}_{\ge 0}^{k}$:

\[ \begin{aligned} \mathbf r - \mathbf q &= -\bm \eta \\ \mathbf r + \mathbf w &= \bm \eta. \end{aligned}\]

Let us now define four additional non-negative variables:

\[ \mathbf x_x \in \mathbb {R}_{\ge 0}^{n}, \;\;\; \mathbf x_y \in \mathbb {R}_{\ge 0}^{n}, \;\;\; \mathbf {r}_x \in \mathbb {R}_{\ge 0}^{k}, \;\;\; \mathbf {r}_y \in \mathbb {R}_{\ge 0}^{k},\]

where:

\[ \mathbf x = \mathbf x_x - \mathbf x_y, \;\;\; \mathbf r = \mathbf {r}_x - \mathbf {r}_y\\ \mathbf {r}_x = \cfrac{1}{2} \mathbf q, \;\;\; \mathbf {r}_y = \cfrac{1}{2} \mathbf w.\]

Then, the above two equalities become:

\[ \begin{aligned} \mathbf r - 2\mathbf {r}_x &= -2\bm \eta \\ \mathbf r + 2 \mathbf {r}_y &= 2\bm \eta, \end{aligned}\]

that is:

\[ \begin{aligned} \mathbf {r}_x + \mathbf {r}_y = \bm \eta, \;\;\; \mathbf r = \mathbf {r}_x - \mathbf {r}_y. \end{aligned}\]

Hence, the optimization problem can be written:

\[ \begin{aligned} \text{minimize}& \;\;\; \mathbf a^T (\mathbf {r}_x + \mathbf {r}_y)\\ \text{subject\;to}& \;\;\; \mathbf{H}(\mathbf x_x - \mathbf x_y) + \mathbf {r}_x - \mathbf {r}_y = \mathbf{z} \\ & \;\;\; \mathbf x_x \succeq \mathbf 0, \; \mathbf x_y \succeq \mathbf 0 \\ & \;\;\; \mathbf {r}_x \succeq \mathbf 0, \; \mathbf {r}_y \succeq \mathbf 0. \end{aligned}\]

To form the above optimization problem, the user can call the following function:

using Ipopt

analysis = pmuLavStateEstimation(system, device, Ipopt.Optimizer)

Then the user can solve the optimization problem by:

solve!(system, analysis)

As a result, we obtain optimal values for the four additional non-negative variables, while the state estimator is obtained by:

\[ \hat{\mathbf x} = \mathbf x_x - \mathbf x_y.\]

Users can retrieve the estimated bus voltage magnitudes $\hat{\mathbf V} = [\hat{V}_i]$ and angles $\hat{\bm {\Theta}} = [\hat{\theta}_i]$, $i \in \mathcal{N}$, using:

julia> 𝐕 = analysis.voltage.magnitude3-element Vector{Float64}:
 0.9999999253259375
 0.8745000982594181
 0.8963710121324888
julia> 𝚯 = analysis.voltage.angle3-element Vector{Float64}:
  1.4260877988717208e-9
 -0.15289998429150575
 -0.21376512970247452

Power Analysis

Once the computation of voltage magnitudes and angles at each bus is completed, various electrical quantities can be determined. JuliaGrid offers the power! function, which enables the calculation of powers associated with buses and branches. Here is an example code snippet demonstrating its usage:

power!(system, analysis)

The function stores the computed powers in the rectangular coordinate system. It calculates the following powers related to buses and branches:

Type	Power	Active	Reactive
Bus	Injections	$\mathbf P = [P_i]$	$\mathbf Q = [Q_i]$
Bus	Generator injections	$\mathbf P_\mathrm{p} = [P_{\mathrm{p}i}]$	$\mathbf Q_\mathrm{p} = [Q_{\mathrm{p}i}]$
Bus	Shunt elements	$\mathbf P_\mathrm{sh} = [{P}_{\mathrm{sh}i}]$	$\mathbf Q_\mathrm{sh} = [{Q}_{\mathrm{sh}i}]$
Branch	From-bus end flows	$\mathbf P_\mathrm{i} = [P_{ij}]$	$\mathbf Q_\mathrm{i} = [Q_{ij}]$
Branch	To-bus end flows	$\mathbf P_\mathrm{j} = [P_{ji}]$	$\mathbf Q_\mathrm{j} = [Q_{ji}]$
Branch	Shunt elements	$\mathbf P_\mathrm{s} = [P_{\mathrm{s}ij}]$	$\mathbf Q_\mathrm{s} = [P_{\mathrm{s}ij}]$
Branch	Series elements	$\mathbf P_\mathrm{l} = [P_{\mathrm{l}ij}]$	$\mathbf Q_\mathrm{l} = [Q_{\mathrm{l}ij}]$

Info

For a clear comprehension of the equations, symbols presented in this section, as well as for a better grasp of power directions, please refer to the Unified Branch Model.

Power Injections

Active and reactive power injections are stored as the vectors $\mathbf{P} = [P_i]$ and $\mathbf{Q} = [Q_i]$, respectively, and can be retrieved using the following commands:

julia> 𝐏 = analysis.power.injection.active3-element Vector{Float64}:
  0.5057581630622112
  3.394630646283425e-5
 -0.5000097006767218
julia> 𝐐 = analysis.power.injection.reactive3-element Vector{Float64}:
  0.3478155281168221
 -0.30002942033952235
 -1.1464065208938168e-5

Generator Power Injections

We can calculate the active and reactive power injections supplied by generators at each bus $i \in \mathcal{N}$ by summing the active and reactive power injections and the active and reactive power demanded by consumers at each bus:

\[ \begin{aligned} P_{\text{p}i} &= P_i + P_{\text{d}i}\\ Q_{\text{p}i} &= Q_i + Q_{\text{d}i}. \end{aligned}\]

The active and reactive power injections from the generators at each bus are stored as vectors, denoted by $\mathbf{P}_{\text{p}} = [P_{\text{p}i}]$ and $\mathbf{Q}_{\text{p}} = [Q_{\text{p}i}]$, which can be obtained using:

julia> 𝐏ₚ = analysis.power.supply.active3-element Vector{Float64}:
  1.0057581630622112
  3.394630646283425e-5
 -9.70067672179109e-6
julia> 𝐐ₚ = analysis.power.supply.reactive3-element Vector{Float64}:
  0.3478155281168221
 -2.9420339522356898e-5
 -1.1464065208938168e-5

Power at Bus Shunt Elements

Active and reactive powers associated with the shunt elements at each bus are represented by the vectors $\mathbf{P}_{\text{sh}} = [{P}_{\text{sh}i}]$ and $\mathbf{Q}_{\text{sh}} = [{Q}_{\text{sh}i}]$. To retrieve these powers in JuliaGrid, use the following commands:

julia> 𝐏ₛₕ = analysis.power.shunt.active3-element Vector{Float64}:
 0.0
 0.0
 0.0
julia> 𝐐ₛₕ = analysis.power.shunt.reactive3-element Vector{Float64}:
 -0.0
 -0.0
 -0.0

Power Flows

The resulting active and reactive power flows at each from-bus end are stored as the vectors $\mathbf{P}_{\text{i}} = [P_{ij}]$ and $\mathbf{Q}_{\text{i}} = [Q_{ij}],$ respectively, and can be retrieved using the following commands:

julia> 𝐏ᵢ = analysis.power.from.active3-element Vector{Float64}:
 0.22926869632849395
 0.27648946673371727
 0.22729627058758717
julia> 𝐐ᵢ = analysis.power.from.reactive3-element Vector{Float64}:
  0.19852792401746913
  0.14928760409935296
 -0.12639765229298713

The vectors of active and reactive power flows at the to-bus end are stored as $\mathbf{P}_\mathrm{j} = [P_{ji}]$ and $\mathbf{Q}_\mathrm{j} = [Q_{ji}]$, respectively, and can be retrieved using the following code:

julia> 𝐏ⱼ = analysis.power.to.active3-element Vector{Float64}:
 -0.2272623242811244
 -0.2743873720776717
 -0.22562232859905024
julia> 𝐐ⱼ = analysis.power.to.reactive3-element Vector{Float64}:
 -0.1736317680465352
 -0.11178390797862477
  0.11177244391341566

Power at Branch Shunt Elements

Active and reactive powers associated with the branch shunt elements at each branch are represented by the vectors $\mathbf{P}_{\text{s}} = [P_{\text{s}ij}]$ and $\mathbf{Q}_{\text{s}} = [Q_{\text{s}ij}]$. We can retrieve these values using the following code:

julia> 𝐏ₛ = analysis.power.charging.active3-element Vector{Float64}:
 0.0
 0.0
 0.0
julia> 𝐐ₛ = analysis.power.charging.reactive3-element Vector{Float64}:
 -0.03529500545015225
 -0.03606961684086606
 -0.03136462826494309

Power at Branch Series Elements

Active and reactive powers associated with the branch series element at each branch are represented by the vectors $\mathbf{P}_{\text{l}} = [P_{\text{l}ij}]$ and $\mathbf{Q}_{\text{l}} = [Q_{\text{l}ij}]$. We can retrieve these values using the following code:

julia> 𝐏ₗ = analysis.power.series.active3-element Vector{Float64}:
 0.0020063720473695416
 0.002102094656045543
 0.0016739419885371514
julia> 𝐐ₗ = analysis.power.series.reactive3-element Vector{Float64}:
 0.06019116142108629
 0.07357331296159413
 0.01673941988537151

Current Analysis

JuliaGrid offers the current! function, which enables the calculation of currents associated with buses and branches. Here is an example code snippet demonstrating its usage:

current!(system, analysis)

The function stores the computed currents in the polar coordinate system. It calculates the following currents related to buses and branches:

Type	Current	Magnitude	Angle
Bus	Injections	$\mathbf I = [I_i]$	$\bm \psi = [\psi_i]$
Branch	From-bus end flows	$\mathbf I_\mathrm{i} = [I_{ij}]$	$\bm \psi_\mathrm{i} = [\psi_{ij}]$
Branch	To-bus end flows	$\mathbf I_\mathrm{j} = [I_{ji}]$	$\bm \psi_\mathrm{j} = [\psi_{ji}]$
Branch	Series elements	$\mathbf I_\mathrm{l} = [I_{\mathrm{l}ij}]$	$\bm \psi_\mathrm{l} = [\psi_{\mathrm{l}ij}]$

Info

For a clear comprehension of the equations, symbols presented in this section, as well as for a better grasp of power directions, please refer to the Unified Branch Model.

Current Injections

In JuliaGrid, complex current injections are stored in the vector of magnitudes denoted as $\mathbf{I} = [I_i]$ and the vector of angles represented as $\bm{\psi} = [\psi_i]$. You can retrieve them using the following commands:

julia> 𝐈 = analysis.current.injection.magnitude3-element Vector{Float64}:
 0.6138135037392803
 0.3430867793578193
 0.557815563020728
julia> 𝛙 = analysis.current.injection.angle3-element Vector{Float64}:
 -0.6024307111825417
  1.4177831992447079
  2.927804596201733

Current Flows

To obtain the vectors of magnitudes $\mathbf{I}_{\text{i}} = [I_{ij}]$ and angles $\bm{\psi}_{\text{i}} = [\psi_{ij}]$ for the resulting complex current flows, you can use the following commands:

julia> 𝐈ᵢ = analysis.current.from.magnitude3-element Vector{Float64}:
 0.3032779013833973
 0.3142184410534215
 0.2974005806144884
julia> 𝛙ᵢ = analysis.current.from.angle3-element Vector{Float64}:
 -0.7136630829982435
 -0.4950864748299688
  0.35460840029890817

Similarly, we can obtain the vectors of magnitudes $\mathbf{I}_{\text{j}} = [I_{ji}]$ and angles $\bm{\psi}_{\text{j}} = [\psi_{ji}]$ of the resulting complex current flows using the following code:

julia> 𝐈ⱼ = analysis.current.to.magnitude3-element Vector{Float64}:
 0.32704429810286606
 0.33053698238603463
 0.2808999987073061
julia> 𝛙ⱼ = analysis.current.to.angle3-element Vector{Float64}:
  2.3362822907000393
  2.5409628708152208
 -2.8953999935529944

Current at Branch Series Elements

To obtain the vectors of magnitudes $\mathbf{I}_{\text{l}} = [I_{\text{l}ij}]$ and angles $\bm{\psi}_{\text{l}} = [\psi_{\text{l}ij}]$ of the resulting complex current flows, one can use the following code:

julia> 𝐈ₗ = analysis.current.series.magnitude3-element Vector{Float64}:
 0.3167311199874069
 0.3241986008641575
 0.2893045098626317
julia> 𝛙ₗ = analysis.current.series.angle3-element Vector{Float64}:
 -0.7614169263002375
 -0.5493964117829945
  0.30174835164585967