PMU State Estimation
To initiate the process, let us construct the PowerSystem composite 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}$, 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. This process of adding measurement devices will be carried out in the State Estimation Model section. Currently, we are only initializing the Measurement type at this stage:
device = measurement()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}_\text{re},\mathbf{V}_\text{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 $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 accuracy of the state estimation 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, where $k = 2|\bar{\mathcal{P}}|$ is the number of 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.
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, varianceMagnitude = 1e-5,
angle = -0.1, 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 [1], 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 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{V}_{i})}, w_{\Im(\bar{V}_{i})}].\]
addPmu!(system, device; label = "V₃, θ₃", bus = 3, magnitude = 0.9, varianceMagnitude = 1e-5,
angle = -0.2, 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 End 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, varianceMagnitude = 1,
angle = 0.4, 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 [1], 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:
\[ \begin{aligned} A = \cfrac{g_{ij} + g_{\text{s}ij}}{\tau_{ij}^2},\;\;\; B = \cfrac{b_{ij}+b_{\text{s}ij}} {\tau_{ij}^2},\;\;\; C = \cfrac{g_{ij}}{\tau_{ij}},\;\;\; D = \cfrac{b_{ij}}{\tau_{ij}}. \end{aligned}\]
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, varianceMagnitude = 1,
angle = -0.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 End 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, varianceMagnitude = 1e-5,
angle = -2.9, 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, varianceMagnitude = 1e-5,
angle = 2.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 = pmuWlsStateEstimation(system, device)Coefficient Matrix
Using the above-described equations, JuliaGrid forms the coefficient matrix $\mathbf{H} \in \mathbb{R}^{k \times 2n}$:
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}_\text{re},\mathbf{V}_\text{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.87, angle = -0.15)
addPmu!(system, device; from = 1, magnitude = 0.30, angle = -0.71)
addPmu!(system, device; from = 2, magnitude = 0.31, angle = -0.49)Following this, we recreate the WLS state estimation model:
analysis = pmuWlsStateEstimation(system, device)Upon inspection, it becomes evident that the precision matrix maintains a diagonal structure:
julia> 𝐖 = analysis.method.precision8×8 SparseArrays.SparseMatrixCSC{Float64, Int64} with 8 stored entries: 100000.0 ⋅ ⋅ ⋅ … ⋅ ⋅ ⋅ ⋅ 100000.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 1.00546e5 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 131177.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ … 2.09799e5 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 1.25032e5 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 3.37492e5
Mean Vector
To retrieve the vector $\mathbf z$, containing the means of Gaussian distributions for each measurement, users can utilize:
julia> 𝐳 = analysis.method.mean8-element Vector{Float64}: 1.0 0.0 0.8602308378043567 -0.13001117525203132 0.22750856279715245 -0.19555013130646098 0.2735231861691376 -0.14589402533305898
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.\]
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 17 stored entries: 1.0 ⋅ ⋅ ⋅ ⋅ ⋅ -0.0470534 1.0 ⋅ ⋅ ⋅ ⋅ -0.522908 -8.10654e-5 1.0 ⋅ ⋅ ⋅ ⋅ ⋅ -1.78244 1.0 ⋅ ⋅ ⋅ ⋅ 0.0154393 -0.00521047 1.0 ⋅ 0.0165651 -0.344281 -0.0110068 0.0036689 -0.522906 1.0julia> 𝐔 = analysis.method.factorization.U6×6 SparseArrays.SparseMatrixCSC{Float64, Int64} with 17 stored entries: 0.494535 -0.00889224 -0.487608 ⋅ ⋅ 0.00896556 ⋅ 0.479929 -0.000191972 ⋅ ⋅ -0.47321 ⋅ ⋅ 0.25474 -0.219243 0.00155842 -0.00162742 ⋅ ⋅ ⋅ 0.148187 -0.000633631 0.000653543 ⋅ ⋅ ⋅ ⋅ 0.561653 -0.430201 ⋅ ⋅ ⋅ ⋅ ⋅ 0.130906
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.9972846475719777 0.8728310555274105 0.89548122281777julia> 𝚯 = analysis.voltage.angle3-element Vector{Float64}: 0.0011070674719133831 -0.1504959040059231 -0.21044407465996606
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 [2, 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 = pmuWlsStateEstimation(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.Q8×8 SuiteSparse.SPQR.QRSparseQ{Float64, Int64}: -0.273505 0.206182 -0.0052897 … 0.388626 0.551877 -0.0193801 -0.0193562 -0.747569 -0.28055 0.184628 0.651606 0.664687 -0.0210875 0.154778 0.237136 -0.0124726 0.00940244 -0.560241 -6.07153e-18 -3.46945e-17 0.707159 -0.533091 -0.0129838 5.55112e-17 1.11022e-16 0.0 0.480656 -0.00211542 … -0.392041 -0.557228 0.0 0.0 -0.355843 0.576773 -0.406818 0.0 0.0 0.0 -0.509844 0.359315julia> 𝐑 = analysis.method.factorization.R6×6 SparseArrays.SparseMatrixCSC{Float64, Int64} with 21 stored entries: -1156.2 497.316 14.2388 -3.53502 586.583 -10.4591 ⋅ 659.702 -3.91115 -3.88362 -442.195 7.88459 ⋅ ⋅ -888.673 503.067 -2.68756 283.077 ⋅ ⋅ ⋅ -574.977 -2.97106 247.685 ⋅ ⋅ ⋅ ⋅ -384.983 3.24637 ⋅ ⋅ ⋅ ⋅ ⋅ -337.119
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.9972846475719777 0.8728310555274102 0.8954812228177701julia> 𝚯 = analysis.voltage.angle3-element Vector{Float64}: 0.0011070674719133396 -0.15049590400592314 -0.21044407465996617
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 [2, 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 [3].
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 = pmuWlsStateEstimation(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 [2, 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}, i \in \mathcal{M} \}.\]
Users can access this information using the variable:
julia> outlier.maxNormalizedResidual387.72710598252246
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 5"
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.maxNormalizedResidual1.289041523534379
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 [2, 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 [1, 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.9999997979267847 0.8761497035131027 0.8983275932436297julia> 𝚯 = analysis.voltage.angle3-element Vector{Float64}: 2.1386863183372956e-7 -0.151458927244325 -0.21100496414291084
Optimal PMU Placement
JuliaGrid utilizes the optimal PMU placement algorithm proposed in [4]. The optimal positioning of PMUs is framed as an integer linear programming problem, expressed as:
\[ \begin{aligned} \text{minimize}& \;\;\; \sum_{i=1}^n d_i\\ \text{subject\;to}& \;\;\; \mathbf A \mathbf d \ge \mathbf a. \end{aligned}\]
Here, the vector $\mathbf d = [d_1,\dots,d_n]^T$ serves as the optimization variable, where $d_i \in \mathbb{F} = \{0,1\}$ is the PMU placement or a binary decision variable associated with the bus $i \in \mathcal{N}$. The all-one vector $\mathbf a$ is of dimension $n$. The binary connectivity matrix $\mathbf A \in \mathbb{F}^{n \times n}$ can be directly derived from the bus nodal matrix $\mathbf Y$ by converting its entries into binary form [5].
Consequently, we obtain the binary vector $\mathbf d = [d_1,\dots,d_n]^T$, where $d_i = 1$, $i \in \mathcal{N}$, suggests that a PMU should be placed at bus $i$. The primary aim of PMU placement in the power system is to determine a minimal set of PMUs such that the entire system is observable without relying on traditional measurements [4]. Specifically, when we observe $d_i = 1$, it indicates that the PMU is installed at bus $i \in \mathcal{N}$ to measure bus voltage phasor and also to measure all current phasors across branches incident to bus $i$.
Determining the optimal PMU placement involves analyzing the created power system. For example:
using GLPK
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 = 2, 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)
placement = pmuPlacement(system, GLPK.Optimizer)The placement variable contains data regarding the optimal placement of measurements. It lists all buses $i \in \mathcal{N}$ that satisfy $d_i = 1$:
julia> placement.busOrderedCollections.OrderedDict{String, Int64} with 1 entry: "2" => 2
This PMU installed at bus 2 will measure the bus voltage phasor at the corresponding bus and all current phasors at the branches incident to bus 2 located at the from-bus or to-bus ends. These data are stored in the variables:
julia> placement.fromOrderedCollections.OrderedDict{String, Int64} with 1 entry: "3" => 3julia> placement.toOrderedCollections.OrderedDict{String, Int64} with 2 entries: "1" => 1 "2" => 2
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:
| Bus | Active | Reactive |
|---|---|---|
| Injections | $\mathbf{P} = [P_i]$ | $\mathbf{Q} = [Q_i]$ |
| Generator injections | $\mathbf{P}_{\text{p}} = [P_{\text{p}i}]$ | $\mathbf{Q}_{\text{p}} = [Q_{\text{p}i}]$ |
| Shunt elements | $\mathbf{P}_{\text{sh}} = [{P}_{\text{sh}i}]$ | $\mathbf{Q}_{\text{sh}} = [{Q}_{\text{sh}i}]$ |
| Branch | Active | Reactive |
|---|---|---|
| From-bus end flows | $\mathbf{P}_{\text{i}} = [P_{ij}]$ | $\mathbf{Q}_{\text{i}} = [Q_{ij}]$ |
| To-bus end flows | $\mathbf{P}_{\text{j}} = [P_{ji}]$ | $\mathbf{Q}_{\text{j}} = [Q_{ji}]$ |
| Shunt elements | $\mathbf{P}_{\text{s}} = [P_{\text{s}ij}]$ | $\mathbf{P}_{\text{s}} = [P_{\text{s}ij}]$ |
| Series elements | $\mathbf{P}_{\text{l}} = [P_{\text{l}ij}]$ | $\mathbf{Q}_{\text{l}} = [Q_{\text{l}ij}]$ |
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.421663162249556 -0.19530673833832304 -0.2213236577243508julia> 𝐐 = analysis.power.injection.reactive3-element Vector{Float64}: 0.36126037933172706 -0.4503322130044093 0.11258234719120122
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}: 0.921663162249556 -0.19530673833832304 0.2786763422756492julia> 𝐐ₚ = analysis.power.supply.reactive3-element Vector{Float64}: 0.36126037933172706 -0.1503322130044093 0.11258234719120122
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.0julia> 𝐐ₛₕ = 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.22750847062143964 0.19415469162811635 0.22294829657037818julia> 𝐐ᵢ = analysis.power.from.reactive3-element Vector{Float64}: 0.1955501454954226 0.16571023383630398 -0.1278285740859068
Similarly, the vectors of active and reactive power flows at the to-bus end are stored as $\mathbf{P}_{\text{j}} = [P_{ji}]$ and $\mathbf{Q}_{\text{j}} = [Q_{ji}]$, respectively, and can be retrieved using the following code:
julia> 𝐏ⱼ = analysis.power.to.active3-element Vector{Float64}: -0.22554403050865718 -0.19271100440004374 -0.2213236577243508julia> 𝐐ⱼ = analysis.power.to.reactive3-element Vector{Float64}: -0.17196970008834428 -0.15053393883015864 0.11258234719120203
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.0julia> 𝐐ₛ = analysis.power.charging.reactive3-element Vector{Float64}: -0.035352757976394165 -0.035352757976394165 -0.0314926153549798
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.0019644401127824167 0.00144368722807256 0.001624638846027477julia> 𝐐ₗ = analysis.power.series.reactive3-element Vector{Float64}: 0.05893320338347258 0.05052905298253964 0.016246388460274763
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:
| Bus | Magnitude | Angle |
|---|---|---|
| Injections | $\mathbf{I} = [I_i]$ | $\bm{\psi} = [\psi_i]$ |
| Branch | Magnitude | Angle |
|---|---|---|
| From-bus end flows | $\mathbf{I}_{\text{i}} = [I_{ij}]$ | $\bm{\psi}_{\text{i}} = [\psi_{ij}]$ |
| To-bus end flows | $\mathbf{I}_{\text{j}} = [I_{ji}]$ | $\bm{\psi}_{\text{j}} = [\psi_{ji}]$ |
| Series elements | $\mathbf{I}_{\text{l}} = [I_{\text{l}ij}]$ | $\bm{\psi}_{\text{l}} = [\psi_{\text{l}ij}]$ |
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.555255804719925 0.5602470534877696 0.27641608388046557julia> 𝛙 = analysis.current.injection.angle3-element Vector{Float64}: -0.7084001797703948 1.8285495088268773 -2.882032194998405
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.2999999999682708 0.25525663989764175 0.29332239932133297julia> 𝛙ᵢ = analysis.current.from.angle3-element Vector{Float64}: -0.7100000222767826 -0.7065199041774624 0.3691383932409053
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.3237186097855176 0.27910332071770105 0.27641608388046557julia> 𝛙ⱼ = analysis.current.to.angle3-element Vector{Float64}: 2.3387018720258586 2.3270001254272197 -2.882032194998405
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.3134039017611634 0.2686714748603358 0.28501217921586053julia> 𝛙ₗ = analysis.current.series.angle3-element Vector{Float64}: -0.7584140912993005 -0.7631714293387878 0.3157764502116658
References
[1] ISO-IEC-OIML-BIPM: "Guide to the expression of uncertainty in measurement," 1992.
[2] A. Abur and A. Exposito, Power System State Estimation: Theory and Implementation, Taylor & Francis, 2004.
[3] G. Korres, "A distributed multiarea state estimation," IEEE Trans. Power Syst., vol. 26, no. 1, pp. 73–84, Feb. 2011.
[4] B. Gou, "Optimal placement of PMUs by integer linear programming," IEEE Trans. Power Syst., vol. 23, no. 3, pp. 1525–1526, Aug. 2008.
[5] B. Xu and A. Abur, "Observability analysis and measurement placement for systems with PMUs," in Proc. IEEE PES PSCE, New York, NY, 2004, pp. 943-946 vol.2.