AC 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 measurement devices $\mathcal M$ into the graph $\mathcal G$. The AC state estimation includes a set of voltmeters $\mathcal V$, ammeters $\mathcal I$, wattmeters $\mathcal P$, varmeters $\mathcal Q$, and PMUs $\bar{\mathcal P}$, with PMUs being able to integrate into AC state estimation in either rectangular coordinates or polar coordinates. To start, we initialize the Measurement type:
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
In accordance with the AC Model, the AC state estimation treats bus voltages as state variables, which we denoted by $\mathbf x \equiv [\bm \Theta, \mathbf V]^T$. The state vector encompasses two components:
- $\bm \Theta \in \mathbb{R}^{n-1}$, representing bus voltage angles,
- $\mathbf V \in \mathbb{R}^n$, representing bus voltage magnitudes.
Consequently, the total number of state variables is $s = 2n-1$, accounting for the fact that the voltage angle for the slack bus is known.
Within the JuliaGrid framework for AC state estimation, the methodology encompasses bus voltage magnitudes, branch current magnitudes, active powers, reactive powers, and phasor measurements. These measurements contribute to the construction of a nonlinear 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 nonlinear measurement functions, where $k$ is the number of measurements, $\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. It is worth noting that the number of equations in the system is equal to $k = |\mathcal V \cup \mathcal I \cup \mathcal P \cup \mathcal Q| + 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$. These covariances exist only if PMUs are observed in rectangular coordinates and correlation is required.
Hence, the nonlinear system of equations is structured according to the specific devices:
\[ \begin{bmatrix} \mathbf z_\mathcal V\\[3pt] \mathbf z_\mathcal I\\[3pt] \mathbf z_\mathcal P\\[3pt] \mathbf z_\mathcal Q\\[3pt] \mathbf z_{\bar{\mathcal P}} \end{bmatrix} = \begin{bmatrix} \mathbf h_\mathcal V(\mathbf x)\\[3pt] \mathbf h_\mathcal I(\mathbf x)\\[3pt] \mathbf h_\mathcal P(\mathbf x)\\[3pt] \mathbf h_\mathcal Q(\mathbf x)\\[3pt] \mathbf h_{\bar{\mathcal P}}(\mathbf x) \end{bmatrix} + \begin{bmatrix} \mathbf u_\mathcal V\\[3pt] \mathbf u_\mathcal I\\[3pt] \mathbf u_\mathcal P\\[3pt] \mathbf u_\mathcal Q\\[3pt] \mathbf u_{\bar{\mathcal P}} \end{bmatrix}\]
Please note that each error vector, denoted as $\mathbf{u}_i$, where $i \in \{\mathcal V, \mathcal I, \mathcal P, \mathcal Q\}$, is associated with the variance vector $\mathbf{v}_i$. However, for PMUs, the error vector $\mathbf{u}_{\bar{\mathcal P}}$, along with its variance vector $\mathbf{v}_{\bar{\mathcal P}}$, can also be associated with the covariance vector $\mathbf{w}_{\bar{\mathcal P}}$.
In summary, upon user definition of the measurement devices, each $i$-th legacy measurement device is linked to the measurement function $h_i(\mathbf x)$, the corresponding measurement value $z_i$, and the measurement variance $v_i$. Meanwhile, 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}$.
Typically, the AC state estimator is obtained using the Gauss-Newton method or its variation, which involves constructing the Jacobian matrix. Therefore, in addition to the aforementioned elements, we also need Jacobian expressions corresponding to the measurement functions, which are also provided below.
Bus Voltage Magnitude Measurements
When introducing a voltmeter $V_i \in \mathcal V$ at bus $i \in \mathcal N$, users specify the measurement value, variance, and measurement function of vectors:
\[ \mathbf{z}_\mathcal{V} = [z_{V_i}], \;\;\; \mathbf{v}_\mathcal{V} = [v_{V_i}], \;\;\; \mathbf{h}_\mathcal{V}(\mathbf x) = [h_{V_i}(\mathbf x)].\]
For example:
addVoltmeter!(system, device; label = "V₁", bus = 1, magnitude = 1.0, variance = 1e-3)Here, the bus voltage magnitude measurement function is simply defined as:
\[ h_{V_i}(\mathbf x) = V_i,\]
with the following Jacobian expression:
\[ \cfrac{\mathrm \partial{{h_{V_i}(\mathbf x)}}} {\mathrm \partial V_i} = 1.\]
From-Bus Current Magnitude Measurements
When introducing an ammeter at branch $(i,j) \in \mathcal E$, it can be placed at the from-bus end, denoted as $I_{ij} \in \mathcal I$, specifying the measurement value, variance, and measurement function of vectors:
\[ \mathbf{z}_\mathcal{I} = [z_{I_{ij}}], \;\;\; \mathbf{v}_\mathcal{I} = [v_{I_{ij}}], \;\;\; \mathbf{h}_\mathcal{I}(\mathbf x) = [h_{I_{ij}}(\mathbf x)].\]
For example:
addAmmeter!(system, device; label = "I₁₂", from = 1, magnitude = 0.3, variance = 1e-2)Here, following the guidelines outlined in the AC Model, the function defining the current magnitude at the from-bus end is expressed as:
\[ h_{I_{ij}}(\mathbf x) = \sqrt{A_{I_{ij}}V_i^2 + B_{I_{ij}}V_j^2 - 2[C_{I_{ij}} \cos(\theta_{ij} - \phi_{ij}) - D_{I_{ij}}\sin(\theta_{ij} - \phi_{ij})]V_iV_j},\]
where:
\[ \begin{gathered} A_{I_{ij}} = \cfrac{(g_{ij} + g_{\mathrm{s}i})^2 + (b_{ij} + b_{\mathrm{s}i})^2}{\tau_{ij}^4}, \;\;\; B_{I_{ij}} = \cfrac{g_{ij}^2 + b_{ij}^2}{\tau_{ij}^2} \\ C_{I_{ij}} = \cfrac{g_{ij}(g_{ij} + g_{\mathrm{s}i}) + b_{ij}(b_{ij} + b_{\mathrm{s}i})}{\tau_{ij}^3}, \;\;\; D_{I_{ij}} = \cfrac{g_{ij}b_{\mathrm{s}i} - b_{ij}g_{\mathrm{s}i}}{\tau_{ij}^3}. \end{gathered}\]
Jacobian expressions corresponding to the measurement function $h_{I_{ij}}(\mathbf x)$ are defined as follows:
\[ \begin{aligned} \cfrac{\mathrm \partial{h_{I_{ij}}(\mathbf x)}}{\mathrm \partial \theta_i} &=- \cfrac{\mathrm \partial{h_{I_{ij}}(\mathbf x)}}{\mathrm \partial \theta_j} = \cfrac{ [C_{I_{ij}}\sin(\theta_{ij} - \phi_{ij}) + D_{I_{ij}}\cos(\theta_{ij} - \phi_{ij})]V_i V_j}{h_{I_{ij}}(\mathbf x)} \\ \cfrac{\mathrm \partial{h_{I_{ij}}(\mathbf x)}}{\mathrm \partial V_i} &= \cfrac{A_{I_{ij}}V_i - [C_{I_{ij}}\cos(\theta_{ij} - \phi_{ij}) - D_{I_{ij}}\sin(\theta_{ij} - \phi_{ij})]V_j}{h_{I_{ij}}(\mathbf x)} \\ \cfrac{\mathrm \partial{h_{I_{ij}}(\mathbf x)}}{\mathrm \partial V_j} &= \cfrac{B_{I_{ij}}V_j - [C_{I_{ij}}\cos(\theta_{ij} - \phi_{ij}) - D_{I_{ij}}\sin(\theta_{ij} - \phi_{ij})]V_i}{h_{I_{ij}}(\mathbf x)}. \end{aligned}\]
To-Bus Current Magnitude Measurements
In addition to the scenario where we add ammeters at the from-bus end, an ammeter can also be positioned at the to-bus end, denoted as $I_{ji} \in \mathcal I$, specifying the measurement value, variance, and measurement function of vectors:
\[ \mathbf{z}_\mathcal{I} = [z_{I_{ji}}], \;\;\; \mathbf{v}_\mathcal{I} = [v_{I_{ji}}], \;\;\; \mathbf{h}_\mathcal{I}(\mathbf x) = [h_{I_{ji}}(\mathbf x)].\]
For example:
addAmmeter!(system, device; label = "I₂₁", to = 1, magnitude = 0.3, variance = 1e-3)Now, the measurement function is as follows:
\[ h_{I_{ji}}(\mathbf x) = \sqrt{A_{I_{ji}}V_i^2 + B_{I_{ji}}V_j^2 - 2[C_{I_{ji}} \cos(\theta_{ij} - \phi_{ij}) + D_{I_{ji}}\sin(\theta_{ij} - \phi_{ij})]V_iV_j},\]
where:
\[ \begin{gathered} A_{I_{ji}} = \cfrac{g_{ij}^2 + b_{ij}^2}{\tau_{ij}^2}, \;\;\; B_{I_{ji}} = (g_{ij} + g_{\mathrm{s}i})^2 + (b_{ij} + b_{\mathrm{s}i})^2 \\ C_{I_{ji}} = \cfrac{g_{ij}(g_{ij} + g_{\mathrm{s}i}) + b_{ij}(b_{ij} + b_{\mathrm{s}i})}{\tau_{ij}}, \;\;\; D_{I_{ji}} = \cfrac{g_{ij}b_{\mathrm{s}i} - b_{ij}g_{\mathrm{s}i}}{\tau_{ij}}. \end{gathered}\]
Jacobian expressions corresponding to the measurement function $h_{I_{ji}}(\mathbf x)$ are defined as follows:
\[ \begin{aligned} \cfrac{\mathrm \partial{h_{I_{ji}}(\mathbf x)}}{\mathrm \partial \theta_i} &=- \cfrac{\mathrm \partial{h_{I_{ji}}(\mathbf x)}}{\mathrm \partial \theta_j} = \cfrac{[C_{I_{ji}}\sin(\theta_{ij} - \phi_{ij}) - D_{I_{ji}}\cos(\theta_{ij}- \phi_{ij})]V_i V_j}{h_{I_{ji}}(\mathbf x)} \\ \cfrac{\mathrm \partial{h_{I_{ji}}(\mathbf x)}}{\mathrm \partial V_i} &= \cfrac{A_{I_{ji}}V_i - [C_{I_{ji}}\cos(\theta_{ij} - \phi_{ij}) + D_{I_{ji}}\sin(\theta_{ij} - \phi_{ij})]V_j}{h_{I_{ji}}(\mathbf x)} \\ \cfrac{\mathrm \partial{h_{I_{ji}}(\mathbf x)}}{\mathrm \partial V_j} &= \cfrac{B_{I_{ji}}V_j - [C_{I_{ji}}\cos(\theta_{ij} - \phi_{ij}) + D_{I_{ji}}\sin(\theta_{ij} - \phi_{ij})]V_i}{h_{I_{ji}}(\mathbf x)} . \end{aligned}\]
Active Power Injection Measurements
When adding a wattmeter $P_i \in \mathcal P$ at bus $i \in \mathcal N$, users specify that the wattmeter measures active power injection and define measurement value, variance, and measurement function of vectors:
\[ \mathbf{z}_\mathcal{P} = [z_{P_i}], \;\;\; \mathbf{v}_\mathcal{P} = [v_{P_i}], \;\;\; \mathbf{h}_\mathcal{P}(\mathbf x) = [h_{P_i}(\mathbf x)].\]
For example:
addWattmeter!(system, device; label = "P₃", bus = 3, active = -0.5, variance = 1e-3)Here, utilizing the AC Model, we derive the function defining the active power injection as follows:
\[ h_{P_i}(\mathbf x) = V_i\sum\limits_{j \in \mathcal{N}_i} (G_{ij}\cos\theta_{ij} + B_{ij}\sin\theta_{ij}) V_j,\]
where $\mathcal{N}_i$ contains buses incident to bus $i$, including bus $i$, with the following Jacobian expressions:
\[ \begin{aligned} \cfrac{\mathrm \partial{h_{P_i}(\mathbf x)}}{\mathrm \partial \theta_i} &= V_i\sum_{j \in \mathcal{N}_i} (-G_{ij}\sin\theta_{ij} + B_{ij}\cos\theta_{ij})V_j - B_{ii}V_i^2\\ \cfrac{\mathrm \partial{h_{P_i}(\mathbf x)}}{\mathrm \partial \theta_j} &= (G_{ij}\sin\theta_{ij} - B_{ij}\cos\theta_{ij})V_iV_j \\ \cfrac{\mathrm \partial{h_{P_i}(\mathbf x)}}{\mathrm \partial V_i} &= \sum_{j \in \mathcal{N}_i} (G_{ij}\cos\theta_{ij} + B_{ij}\sin\theta_{ij})V_j + G_{ii}V_i\\ \cfrac{\mathrm \partial{h_{P_i}(\mathbf x)}}{\mathrm \partial V_j} &= (G_{ij}\cos\theta_{ij} + B_{ij}\sin\theta_{ij})V_i. \end{aligned}\]
From-Bus Active Power Flow Measurements
Additionally, when introducing a wattmeter at branch $(i,j) \in \mathcal E$, users specify that the wattmeter measures active power flow. It can be positioned at the from-bus end, denoted as $P_{ij} \in \mathcal P$, specifying the measurement value, variance, and measurement function of vectors:
\[ \mathbf{z}_\mathcal{P} = [z_{P_{ij}}], \;\;\; \mathbf{v}_\mathcal{P} = [v_{P_{ij}}], \;\;\; \mathbf{h}_\mathcal{P}(\mathbf x) = [h_{P_{ij}}(\mathbf x)].\]
For example:
addWattmeter!(system, device; label = "P₁₂", from = 1, active = 0.2, variance = 1e-4)Here, the function describing active power flow at the from-bus end is defined as follows:
\[ h_{P_{ij}}(\mathbf x) = \cfrac{g_{ij} + g_{\mathrm{s}i}}{\tau_{ij}^2} V_i^2 - \cfrac{1}{\tau_{ij}} \left[g_{ij}\cos(\theta_{ij} - \phi_{ij}) + b_{ij}\sin(\theta_{ij} - \phi_{ij})\right]V_iV_j,\]
with the following Jacobian expressions:
\[ \begin{aligned} \cfrac{\mathrm \partial{h_{P_{ij}}(\mathbf x)}}{\mathrm \partial \theta_i} &=- \cfrac{\mathrm \partial{h_{P_{ij}}(\mathbf x)}}{\mathrm \partial \theta_j} = \cfrac{1}{\tau_{ij}} \left[g_{ij}\sin(\theta_{ij} - \phi_{ij}) - b_{ij}\cos(\theta_{ij} - \phi_{ij})\right] V_iV_j \\ \cfrac{\mathrm \partial{h_{P_{ij}}(\mathbf x)}}{\mathrm \partial V_{i}} &= 2\cfrac{g_{ij} + g_{\mathrm{s}i}}{\tau_{ij}^2} V_i - \cfrac{1}{\tau_{ij}} \left[g_{ij}\cos(\theta_{ij} - \phi_{ij}) + b_{ij}\sin(\theta_{ij} - \phi_{ij})\right] V_j \\ \cfrac{\mathrm \partial{h_{P_{ij}}(\mathbf x)}}{\mathrm \partial V_j} &= - \cfrac{1}{\tau_{ij}} \left[g_{ij}\cos(\theta_{ij} - \phi_{ij}) + b_{ij}\sin(\theta_{ij} - \phi_{ij})\right] V_i. \end{aligned}\]
To-Bus Active Power Flow Measurements
Similarly, a wattmeter can be placed at the to-bus end, denoted as $P_{ji} \in \mathcal{P}$, specifying the measurement value, variance, and measurement function of vectors:
\[ \mathbf{z}_\mathcal{P} = [z_{P_{ji}}], \;\;\; \mathbf{v}_\mathcal{P} = [v_{P_{ji}}], \;\;\; \mathbf{h}_\mathcal{P}(\mathbf x) = [h_{P_{ji}}(\mathbf x)].\]
For example:
addWattmeter!(system, device; label = "P₂₁", to = 1, active = -0.2, variance = 1e-4)Thus, the function describing active power flow at the to-bus end is defined as follows:
\[ h_{P_{ji}}(\mathbf x) = (g_{ij} + g_{\mathrm{s}i}) V_j^2 - \cfrac{1}{\tau_{ij}} \left[g_{ij}\cos(\theta_{ij} - \phi_{ij}) - b_{ij}\sin(\theta_{ij} - \phi_{ij})\right] V_iV_j,\]
with the following Jacobian expressions:
\[ \begin{aligned} \cfrac{\mathrm \partial{h_{P_{ji}}(\mathbf x)}}{\mathrm \partial \theta_i} &= - \cfrac{\mathrm \partial{h_{P_{ji}}(\mathbf x)}}{\mathrm \partial \theta_j} = \cfrac{1}{\tau_{ij}} \left[g_{ij}\sin(\theta_{ij} - \phi_{ij}) + b_{ij}\cos(\theta_{ij} - \phi_{ij})\right] V_iV_j \\ \cfrac{\mathrm \partial{h_{P_{ji}}(\mathbf x)}}{\mathrm \partial V_{i}} &= - \cfrac{1}{\tau_{ij}} \left[g_{ij}\cos(\theta_{ij} - \phi_{ij}) - b_{ij}\sin(\theta_{ij} - \phi_{ij})\right] V_j \\ \cfrac{\mathrm \partial{h_{P_{ji}}(\mathbf x)}}{\mathrm \partial V_{j}} &= 2(g_{ij} + g_{\mathrm{s}i}) V_j - \cfrac{1}{\tau_{ij}} \left[g_{ij}\cos(\theta_{ij} - \phi_{ij}) - b_{ij}\sin(\theta_{ij} - \phi_{ij})\right] V_i. \end{aligned}\]
Reactive Power Injection Measurements
When adding a varmeter $Q_i \in \mathcal Q$ at bus $i \in \mathcal N$, users specify that the varmeter measures reactive power injection and define the measurement value, variance, and measurement function of vectors:
\[ \mathbf{z}_\mathcal{Q} = [z_{Q_i}], \;\;\; \mathbf{v}_\mathcal{Q} = [v_{Q_i}], \;\;\; \mathbf{h}_\mathcal{Q}(\mathbf x) = [h_{Q_i}(\mathbf x)].\]
For example:
addVarmeter!(system, device; label = "Q₃", bus = 3, reactive = 0, variance = 1e-3)Here, utilizing the AC Model, we derive the function defining the reactive power injection as follows:
\[ h_{Q_i}(\mathbf x) = V_i\sum\limits_{j \in \mathcal{N}_i} (G_{ij}\sin\theta_{ij} - B_{ij}\cos\theta_{ij})V_j,\]
where $\mathcal{N}_i$ contains buses incident to bus $i$, including bus $i$, with the following Jacobian expressions:
\[ \begin{aligned} \cfrac{\mathrm \partial{h_{Q_i}(\mathbf x)}}{\mathrm \partial \theta_i} &= V_i\sum_{j \in \mathcal{N}_i} (G_{ij}\cos\theta_{ij} + B_{ij}\sin\theta_{ij})V_j - G_{ii}V_i^2\\ \cfrac{\mathrm \partial{h_{Q_i}(\mathbf x)}}{\mathrm \partial \theta_j} &= -(G_{ij}\cos\theta_{ij} + B_{ij}\sin\theta_{ij})V_iV_j \\ \cfrac{\mathrm \partial{h_{Q_i}(\mathbf x)}}{\mathrm \partial V_i} &= \sum_{j \in \mathcal{N}_i} (G_{ij}\sin\theta_{ij} - B_{ij}\cos\theta_{ij})V_j - B_{ii}V_i\\ \cfrac{\mathrm \partial{h_{Q_i}(\mathbf x)}}{\mathrm \partial V_j} &= (G_{ij}\sin\theta_{ij} - B_{ij}\cos\theta_{ij})V_i. \end{aligned}\]
From-Bus Reactive Power Flow Measurements
Additionally, when introducing a varmeter at branch $(i,j) \in \mathcal E$, users specify that the varmeter measures reactive power flow. It can be positioned at the from-bus end, denoted as $Q_{ij} \in \mathcal Q$, with its measurement value, variance, and measurement function included in vectors:
\[ \mathbf{z}_\mathcal{Q} = [z_{Q_{ij}}], \;\;\; \mathbf{v}_\mathcal{Q} = [v_{Q_{ij}}], \;\;\; \mathbf{h}_\mathcal{Q}(\mathbf x) = [h_{Q_{ij}}(\mathbf x)].\]
For example:
addVarmeter!(system, device; label = "Q₁₂", from = 1, reactive = 0.2, variance = 1e-4)Here, the function describing reactive power flow at the from-bus end is defined as follows:
\[ h_{Q_{ij}}(\mathbf x) = -\cfrac{b_{ij} + b_{\mathrm{s}i}}{\tau_{ij}^2} V_i^2 - \cfrac{1}{\tau_{ij}} \left[g_{ij}\sin(\theta_{ij} - \phi_{ij}) - b_{ij}\cos(\theta_{ij} - \phi_{ij})\right] V_iV_j,\]
with the following Jacobian expressions:
\[ \begin{aligned} \cfrac{\mathrm \partial{h_{Q_{ij}}(\mathbf x)}}{\mathrm \partial \theta_i} &= - \cfrac{\mathrm \partial{h_{Q_{ij}}(\mathbf x)}}{\mathrm \partial \theta_j} = - \cfrac{1}{\tau_{ij}} \left[g_{ij}\cos(\theta_{ij} - \phi_{ij}) + b_{ij}\sin(\theta_{ij} - \phi_{ij})\right] V_iV_j \\ \cfrac{\mathrm \partial{h_{Q_{ij}}(\mathbf x)}}{\mathrm \partial V_i} &= - 2\cfrac{b_{ij} + b_{\mathrm{s}i}}{\tau_{ij}^2} V_i - \cfrac{1}{\tau_{ij}} \left[g_{ij}\sin(\theta_{ij} - \phi_{ij}) - b_{ij}\cos(\theta_{ij} - \phi_{ij})\right] V_j\\ \cfrac{\mathrm \partial{h_{Q_{ij}}(\mathbf x)}}{\mathrm \partial V_j} &= - \cfrac{1}{\tau_{ij}} \left[g_{ij}\sin(\theta_{ij} - \phi_{ij}) - b_{ij}\cos(\theta_{ij} - \phi_{ij})\right] V_i. \end{aligned}\]
To-Bus Reactive Power Flow Measurements
Similarly, a varmeter can be placed at the to-bus end, denoted as $Q_{ji} \in \mathcal Q$, with its own measurement value, variance, and measurement function included in vectors:
\[ \mathbf{z}_\mathcal{Q} = [z_{Q_{ji}}], \;\;\; \mathbf{v}_\mathcal{Q} = [v_{Q_{ji}}], \;\;\; \mathbf{h}_\mathcal{Q}(\mathbf x) = [h_{Q_{ji}}(\mathbf x)].\]
For example:
addVarmeter!(system, device; label = "Q₂₁", to = 1, reactive = -0.2, variance = 1e-4)Thus, the function describing reactive power flow at the to-bus end is defined as follows:
\[ h_{Q_{ji}}(\mathbf x) = -(b_{ij} + b_{\mathrm{s}i}) V_j^2 + \cfrac{1}{\tau_{ij}} \left[g_{ij}\sin(\theta_{ij} - \phi_{ij}) + b_{ij}\cos(\theta_{ij} - \phi_{ij})\right] V_iV_j,\]
with the following Jacobian expressions:
\[ \begin{aligned} \cfrac{\mathrm \partial{h_{Q_{ji}}(\mathbf x)}}{\mathrm \partial \theta_i} &= - \cfrac{\mathrm \partial{h_{Q_{ji}}(\mathbf x)}}{\mathrm \partial \theta_j} = \cfrac{1}{\tau_{ij}} \left[g_{ij}\cos(\theta_{ij} - \phi_{ij}) - b_{ij}\sin(\theta_{ij} - \phi_{ij})\right] V_iV_j \\ \cfrac{\mathrm \partial{h_{Q_{ji}}(\mathbf x)}}{\mathrm \partial V_i} &= \cfrac{1}{\tau_{ij}} \left[g_{ij}\sin(\theta_{ij} - \phi_{ij}) + b_{ij}\cos(\theta_{ij} - \phi_{ij})\right] V_j\\ \cfrac{\mathrm \partial{h_{Q_{ji}}(\mathbf x)}}{\mathrm \partial V_j} &= -2(b_{ij} + b_{\mathrm{s}i}) V_{j} + \cfrac{1}{\tau_{ij}} \left[g_{ij}\sin(\theta_{ij} - \phi_{ij}) + b_{ij}\cos(\theta_{ij} - \phi_{ij})\right] V_i. \end{aligned}\]
Rectangular Bus Voltage Phasor Measurements
When a PMU $(V_i, \theta_i) \in \bar{\mathcal P}$ is introduced at bus $i \in \mathcal N$, it will be incorporated into the AC state estimation model using rectangular coordinates by default. It will define the measurement values, variances, and measurement functions of vectors:
\[ \mathbf{z}_{\bar{\mathcal P}} = [z_{\Re(\bar{V}_i)}, z_{\Im(\bar{V}_i)}], \;\;\; \mathbf{v}_{\bar{\mathcal P}} = [v_{\Re(\bar{V}_i)}, v_{\Im(\bar{V}_i)}], \;\;\; \mathbf{h}_{\bar{\mathcal P}}(\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) = V_i \cos \theta_i\\ h_{\Im(\bar{V}_i)}(\mathbf x) = V_i \sin \theta_i. \end{aligned}\]
Jacobian expressions corresponding to the measurement function $h_{\Re(\bar{V}_i)}(\mathbf x)$ are defined as follows:
\[ \begin{aligned} \cfrac{\mathrm \partial{h_{\Re(\bar{V}_i)}(\mathbf x)}}{\mathrm \partial \theta_i} = -V_i \sin \theta_i, \;\;\; \cfrac{\mathrm \partial{h_{\Re(\bar{V}_i)}(\mathbf x)}}{\mathrm \partial V_i} = \cos \theta_i, \end{aligned}\]
while Jacobian expressions corresponding to the measurement function $h_{\Im(\bar{V}_i)}(\mathbf x)$ are:
\[ \begin{aligned} \cfrac{\mathrm \partial{h_{\Im(\bar{V}_i)}(\mathbf x)}}{\mathrm \partial \theta_i} = V_i \cos \theta_i,\;\;\; \cfrac{\mathrm \partial{h_{\Im(\bar{V}_i)}(\mathbf x)}}{\mathrm \partial V_i} = \sin \theta_i. \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}_{\bar{\mathcal P}} = [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).\]
Polar Bus Voltage Phasor Measurements
If the user chooses to include phasor measurement $(V_i, \theta_i) \in \bar{\mathcal P}$ in polar coordinates in the AC state estimation model, the user will specify the measurement values, variances, and measurement functions of vectors:
\[ \mathbf{z}_{\bar{\mathcal P}} = [z_{V_i}, z_{\theta_i}], \;\;\; \mathbf{v}_{\bar{\mathcal P}} = [v_{V_i}, v_{\theta_i}], \;\;\; \mathbf{h}_{\bar{\mathcal P}}(\mathbf x) = [h_{V_{i}}(\mathbf x), h_{\theta_{i}}(\mathbf x)].\]
For example:
addPmu!(system, device; label = "V₁, θ₁", bus = 1, magnitude = 1.0, angle = 0,
varianceMagnitude = 1e-5, varianceAngle = 1e-6, polar = true)Here, the functions defining the bus voltage phasor measurement are straightforward:
\[ \begin{aligned} h_{V_i}(\mathbf x) = V_i\\ h_{\theta_i}(\mathbf x) = \theta_i, \end{aligned}\]
with the following Jacobian expressions:
\[ \begin{aligned} \cfrac{\mathrm \partial{{h_{V_i}(\mathbf x)}}}{\mathrm \partial V_i} = 1, \;\;\; \cfrac{\mathrm \partial{{h_{\theta_i}(\mathbf x)}}}{\mathrm \partial \theta_i} = 1. \end{aligned}\]
Rectangular From-Bus Current Phasor Measurements
When introducing a PMU at branch $(i,j) \in \mathcal E$, it can be placed at the from-bus end, denoted as $(I_{ij}, \psi_{ij}) \in \bar{\mathcal P}$, and it will be integrated into the AC state estimation model using rectangular coordinates by default. Incorporating current phasor measurements in the polar coordinate system is highly susceptible to ill-conditioned problems, especially when dealing with small values of current magnitudes. This is the reason why we typically include PMUs in the rectangular coordinate system by default.
Therefore, here we specify the measurement values, variances, and measurement functions of vectors:
\[ \mathbf{z}_{\bar{\mathcal P}} = [z_{\Re(\bar{I}_{ij})}, z_{\Im(\bar{I}_{ij})}], \;\;\; \mathbf{v}_{\bar{\mathcal P}} = [v_{\Re(\bar{I}_{ij})}, v_{\Im(\bar{I}_{ij})}], \;\;\; \mathbf{h}_{\bar{\mathcal P}}(\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-5, 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_{\psi_{ij}} \cos\theta_i - B_{\psi_{ij}} \sin\theta_i)V_i - [C_{\psi_{ij}} \cos(\theta_j + \phi_{ij}) - D_{\psi_{ij}} \sin(\theta_j + \phi_{ij})]V_j, \\ h_{\Im(\bar{I}_{ij})}(\mathbf x) &= (A_{\psi_{ij}} \sin\theta_i + B_{\psi_{ij}} \cos\theta_i)V_i - [C_{\psi_{ij}} \sin(\theta_j + \phi_{ij}) + D_{\psi_{ij}} \cos(\theta_j + \phi_{ij})]V_j. \end{aligned}\]
Jacobian expressions corresponding to the measurement function $h_{\Re(\bar{I}_{ij})}(\mathbf x)$ are defined as follows:
\[ \begin{aligned} \cfrac{\mathrm \partial{h_{\Re(\bar{I}_{ij})}(\mathbf x)}}{\mathrm \partial \theta_i} &= -(A_{\psi_{ij}} \sin\theta_{i} + B_{\psi_{ij}} \cos\theta_i)V_i\\ \cfrac{\mathrm \partial{h_{\Re(\bar{I}_{ij})}(\mathbf x)}}{\mathrm \partial \theta_j} &= [C_{\psi_{ij}} \sin(\theta_j + \phi_{ij}) + D_{\psi_{ij}} \cos(\theta_j + \phi_{ij})] V_j \\ \cfrac{\mathrm \partial{h_{\Re(\bar{I}_{ij})}(\mathbf x)}}{\mathrm \partial V_i} &= A_{\psi_{ij}} \cos\theta_{i} - B_{\psi_{ij}} \sin\theta_i\\ \cfrac{\mathrm \partial{h_{\Re(\bar{I}_{ij})}(\mathbf x)}}{\mathrm \partial V_j} &= -C_{\psi_{ij}} \cos(\theta_j + \phi_{ij}) + D_{\psi_{ij}} \sin(\theta_j + \phi_{ij}). \end{aligned}\]
while Jacobian expressions corresponding to the measurement function $h_{\Im(\bar{I}_{ij})}(\mathbf x)$ are:
\[ \begin{aligned} \cfrac{\mathrm \partial{h_{\Im(\bar{I}_{ij})}(\mathbf x)}}{\mathrm \partial \theta_i} &= (A_{\psi_{ij}} \cos \theta_{i} - B_{\psi_{ij}} \sin \theta_i)V_i\\ \cfrac{\mathrm \partial{h_{\Im(\bar{I}_{ij})}(\mathbf x)}}{\mathrm \partial \theta_j} &= [-C_{\psi_{ij}} \cos(\theta_j + \phi_{ij}) + D_{\psi_{ij}} \sin(\theta_j + \phi_{ij})] V_j \\ \cfrac{\mathrm \partial{h_{\Im(\bar{I}_{ij})}(\mathbf x)}}{\mathrm \partial V_i} &= A_{\psi_{ij}} \sin \theta_{i} + B_{\psi_{ij}} \cos\theta_i\\ \cfrac{\mathrm \partial{h_{\Im(\bar{I}_{ij})}(\mathbf x)}}{\mathrm \partial V_j} &= -C_{\psi_{ij}} \sin(\theta_j + \phi_{ij}) - D_{\psi_{ij}} \cos(\theta_j + \phi_{ij}). \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}_{\bar{\mathcal{P}}} = [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-4, 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).\]
Polar From-Bus Current Phasor Measurements
If the user chooses to include phasor measurement $(I_{ij}, \psi_{ij}) \in \bar{\mathcal P}$ in polar coordinates in the AC state estimation model, the user will specify the measurement values, variances, and measurement functions of vectors:
\[ \mathbf{z}_{\bar{\mathcal P}} = [z_{I_{ij}}, z_{\psi_{ij}}], \;\;\; \mathbf{v}_{\bar{\mathcal P}} = [v_{I_{ij}}, v_{\psi_{ij}}], \;\;\; \mathbf{h}_{\bar{\mathcal P}}(\mathbf x) = [h_{I_{ij}}(\mathbf x), h_{\psi_{ij}}(\mathbf x)].\]
For example:
addPmu!(system, device; label = "I₁₂, ψ₁₂", from = 1, magnitude = 0.3, angle = -0.7,
varianceMagnitude = 1e-5, varianceAngle = 1e-4, polar = true)Here, the function associated with the branch current magnitude at the from-bus end remains identical to the one provided in From-Bus Current Magnitude Measurements. However, the function defining the branch current angle measurement is expressed as:
\[ h_{\psi_{ij}}(\mathbf x) = \mathrm{atan} \Bigg[ \cfrac{(A_{\psi_{ij}} \sin\theta_i + B_{\psi_{ij}} \cos\theta_i)V_i - [C_{\psi_{ij}} \sin(\theta_j + \phi_{ij}) + D_{\psi_{ij}}\cos(\theta_j + \phi_{ij})]V_j} {(A_{\psi_{ij}} \cos\theta_i - B_{\psi_{ij}} \sin\theta_i)V_i - [C_{\psi_{ij}} \cos(\theta_j + \phi_{ij}) - D_{\psi_{ij}} \sin(\theta_j + \phi_{ij})]V_j} \Bigg],\]
where:
\[ A_{\psi_{ij}} = \cfrac{g_{ij} + g_{\mathrm{s}i}}{\tau_{ij}^2}, \;\;\; B_{\psi_{ij}} = \cfrac{b_{ij}+b_{\mathrm{s}i}}{\tau_{ij}^2}, \;\;\; C_{\psi_{ij}} = \cfrac{g_{ij}}{\tau_{ij}}, \;\;\; D_{\psi_{ij}} = \cfrac{b_{ij}}{\tau_{ij}}.\]
Jacobian expressions associated with the branch current magnitude function $h_{I_{ij}}(\mathbf x)$ remains identical to the one provided in From-Bus Current Magnitude Measurements. Further, Jacobian expressions corresponding to the measurement function $h_{\psi_{ij}}(\mathbf x)$ are defined as follows:
\[ \begin{aligned} \cfrac{\mathrm \partial{h_{\psi_{ij}}(\mathbf x)}}{\mathrm \partial \theta_i} &= \cfrac{A_{I_{ij}} V_i^2- [C_{I_{ij}} \cos(\theta_{ij}- \phi_{ij}) - D_{I_{ij}} \sin (\theta_{ij} - \phi_{ij})]V_iV_j}{h_{{I}_{ij}}^2(\mathbf x)} \\ \cfrac{\mathrm \partial{h_{\psi_{ij}}(\mathbf x)}}{\mathrm \partial \theta_j} &= \cfrac{B_{I_{ij}} V_j^2 - [C_{I_{ij}} \cos (\theta_{ij} - \phi_{ij}) - D_{I_{ij}} \sin(\theta_{ij}- \phi_{ij})]V_iV_j}{h_{{I}_{ij}}^2(\mathbf x)} \\ \cfrac{\mathrm \partial{h_{\psi_{ij}}(\mathbf x)}}{\mathrm \partial V_i} &= - \cfrac{[C_{I_{ij}} \sin (\theta_{ij} - \phi_{ij}) + D_{I_{ij}} \cos(\theta_{ij}- \phi_{ij})]V_j }{h_{{I}_{ij}}^2(\mathbf x)}\\ \cfrac{\mathrm \partial{h_{\psi_{ij}}(\mathbf x)}}{\mathrm \partial V_j} &= \cfrac{[C_{I_{ij}} \sin (\theta_{ij} - \phi_{ij}) + D_{I_{ij}} \cos(\theta_{ij}- \phi_{ij})]V_i }{h_{{I}_{ij}}^2(\mathbf x)}. \end{aligned}\]
Rectangular To-Bus Current Phasor Measurements
When introducing a PMU at branch $(i,j) \in \mathcal E$, it can be placed at the to-bus end, denoted as $(I_{ji}, \psi_{ji}) \in \bar{\mathcal P}$, and it will be integrated into the AC state estimation model using rectangular coordinates by default. The user will specify the measurement values, variances, and measurement functions of vectors:
\[ \mathbf{z}_{\bar{\mathcal P}} = [z_{\Re(\bar{I}_{ji})}, z_{\Im(\bar{I}_{ji})}], \;\;\; \mathbf{v}_{\bar{\mathcal P}} = [v_{\Re(\bar{I}_{ji})}, v_{\Im(\bar{I}_{ji})}], \;\;\; \mathbf{h}_{\bar{\mathcal P}}(\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) &= (A_{\psi_{ji}} \cos\theta_j - B_{\psi_{ji}} \sin\theta_j)V_j - [C_{\psi_{ji}} \cos(\theta_i - \phi_{ij}) - D_{\psi_{ji}} \sin (\theta_i + \phi_{ij})]V_i\\ h_{\Im(\bar{I}_{ji})}(\mathbf x) &= (A_{\psi_{ji}} \sin\theta_j + B_{\psi_{ji}} \cos\theta_j)V_j - [C_{\psi_{ji}} \sin(\theta_i + \phi_{ij}) + D_{\psi_{ji}} \cos (\theta_i + \phi_{ij})]V_i. \end{aligned}\]
Jacobian expressions corresponding to the measurement function $h_{\Re(\bar{I}_{ji})}(\mathbf x)$ are defined as follows:
\[ \begin{aligned} \cfrac{\mathrm \partial{h_{\Re(\bar{I}_{ji})}(\mathbf x)}}{\mathrm \partial \theta_i} &= [C_{\psi_{ji}} \sin (\theta_{i} - \phi_{ij}) + D_{\psi_{ji}} \cos (\theta_i - \phi_{ij})]V_i\\ \cfrac{\mathrm \partial{h_{\Re(\bar{I}_{ji})}(\mathbf x)}}{\mathrm \partial \theta_j} &= -(A_{\psi_{ji}} \sin\theta_j + B_{\psi_{ji}} \cos \theta_j ) V_j \\ \cfrac{\mathrm \partial{h_{\Re(\bar{I}_{ji})}(\mathbf x)}}{\mathrm \partial V_i} &= - C_{\psi_{ji}} \cos (\theta_i - \phi_{ij}) + D_{\psi_{ji}} \sin(\theta_i - \phi_{ij})\\ \cfrac{\mathrm \partial{h_{\Re(\bar{I}_{ji})}(\mathbf x)}}{\mathrm \partial V_j} &= A_{\psi_{ji}} \cos \theta_j - B_{\psi_{ji}} \sin \theta_j. \end{aligned}\]
while Jacobian expressions corresponding to the measurement function $h_{\Im(\bar{I}_{ji})}(\mathbf x)$ are:
\[ \begin{aligned} \cfrac{\mathrm \partial{h_{\Im(\bar{I}_{ji})}(\mathbf x)}}{\mathrm \partial \theta_i} &= [-C_{\psi_{ji}} \cos (\theta_{i} - \phi_{ij}) + D_{\psi_{ji}} \sin (\theta_i - \phi_{ij})]V_i\\ \cfrac{\mathrm \partial{h_{\Im(\bar{I}_{ji})}(\mathbf x)}}{\mathrm \partial \theta_j} &= (A_{\psi_{ji}} \cos\theta_j - B_{\psi_{ji}} \sin \theta_j ) V_j \\ \cfrac{\mathrm \partial{h_{\Im(\bar{I}_{ji})}(\mathbf x)}}{\mathrm \partial V_i} &= - C_{\psi_{ji}} \sin (\theta_i - \phi_{ij}) - D_{\psi_{ji}} \cos(\theta_i - \phi_{ij})\\ \cfrac{\mathrm \partial{h_{\Im(\bar{I}_{ji})}(\mathbf x)}}{\mathrm \partial V_j} &= A_{\psi_{ji}} \sin \theta_j + B_{\psi_{ji}} \cos \theta_j. \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}_{\bar{\mathcal{P}}} = [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).\]
Polar To-Bus Current Phasor Measurements
If the user chooses to include phasor measurement $(I_{ji}, \psi_{ji}) \in \bar{\mathcal P}$ in polar coordinates in the AC state estimation model, the user will specify the measurement values, variances, and measurement functions of vectors:
\[ \mathbf{z}_{\bar{\mathcal P}} = [z_{I_{ji}}, z_{\psi_{ji}}], \;\;\; \mathbf{v}_{\bar{\mathcal P}} = [v_{I_{ji}}, v_{\psi_{ji}}], \;\;\; \mathbf{h}_{\bar{\mathcal P}}(\mathbf x) = [h_{I_{ji}}(\mathbf x), h_{\psi_{ji}}(\mathbf x)].\]
For example:
addPmu!(system, device; label = "I₂₁, ψ₂₁", to = 1, magnitude = 0.3, angle = 2.3,
varianceMagnitude = 1e-2, varianceAngle = 1e-3, polar = true)Here, the function associated with the branch current magnitude at the to-bus end remains identical to the one provided in To-Bus Current Magnitude Measurements. However, the function defining the branch current angle measurement is expressed as:
\[ h_{\psi_{ji}}(\mathbf {x}) = \mathrm{atan}\Bigg[ \cfrac{(A_{\psi_{ji}} \sin\theta_j + B_{\psi_{ji}} \cos\theta_j)V_j - [C_{\psi_{ji}} \sin(\theta_i - \phi_{ij}) + D_{\psi_{ji}}\cos(\theta_i - \phi_{ij})]V_i} {(A_{\psi_{ji}} \cos\theta_j - B_{\psi_{ji}} \sin\theta_j)V_j - [C_{\psi_{ji}} \cos(\theta_i - \phi_{ij}) - D_{\psi_{ji}} \sin(\theta_i - \phi_{ij})]V_i} \Bigg],\]
where:
\[ A_{\psi_{ji}} = g_{ij} + g_{\mathrm{s}i}, \;\;\; B_{\psi_{ji}} = b_{ij} + b_{\mathrm{s}i}, \;\;\; C_{\psi_{ji}} = \cfrac{g_{ij}}{\tau_{ij}}, \;\;\; D_{\psi_{ji}} = \cfrac{b_{ij}}{\tau_{ij}}.\]
Jacobian expressions associated with the branch current magnitude function $h_{I_{ji}}(\mathbf x)$ remains identical to the one provided in To-Bus Current Magnitude Measurements. Further, Jacobian expressions corresponding to the measurement function $h_{\psi_{ji}}(\mathbf x)$ are defined as follows:
\[ \begin{aligned} \cfrac{\mathrm \partial{h_{\psi_{ji}}(\mathbf x)}}{\mathrm \partial \theta_i} &= \cfrac{A_{I_{ji}} V_i^2- [C_{I_{ji}} \cos(\theta_{ij}- \phi_{ij}) + D_{I_{ji}} \sin (\theta_{ij} - \phi_{ij}) ]V_iV_j}{h_{{I}_{ji}}^2(\mathbf x)} \\ \cfrac{\mathrm \partial{h_{\psi_{ji}}(\mathbf x)}}{\mathrm \partial \theta_j} &= \cfrac{B_{I_{ji}} V_j^2 - [C_{I_{ji}} \cos (\theta_{ij} - \phi_{ij}) + D_{I_{ji}} \sin(\theta_{ij}- \phi_{ij})]V_iV_j}{h_{{I}_{ji}}^2(\mathbf x)} \\ \cfrac{\mathrm \partial{h_{\psi_{ji}}(\mathbf x)}}{\mathrm \partial V_i} &= -\cfrac{[C_{I_{ji}} \sin (\theta_{ij} - \phi_{ij}) - D_{I_{ji}} \cos(\theta_{ij}- \phi_{ij})]V_j }{h_{{I}_{ji}}^2(\mathbf x)}\\ \cfrac{\mathrm \partial{h_{\psi_{ji}}(\mathbf x)}}{\mathrm \partial V_j} &= \cfrac{[C_{I_{ji}} \sin (\theta_{ij} - \phi_{ij}) - D_{I_{ji}} \cos(\theta_{ij}- \phi_{ij})]V_i }{h_{{I}_{ji}}^2(\mathbf x)}. \end{aligned}\]
Weighted Least-Squares Estimation
Given the available set of measurements $\mathcal M$, the weighted least-squares (WLS) estimator $\hat{\mathbf x}$ can be found using the Gauss-Newton method:
\[\begin{gathered} \Big[\mathbf J (\mathbf x^{(\nu)})^T \bm \Sigma^{-1} \mathbf J (\mathbf x^{(\nu)})\Big] \mathbf \Delta \mathbf x^{(\nu)} = \mathbf J (\mathbf x^{(\nu)})^T \bm \Sigma^{-1} \mathbf r (\mathbf x^{(\nu)}) \\ \mathbf x^{(\nu+1)} = \mathbf x^{(\nu)} + \mathbf \Delta \mathbf x^{(\nu)}, \end{gathered}\]
where $\nu = \{0,1,2,\dots\}$ is the iteration index, $\mathbf \Delta \mathbf x \in \mathbb {R}^s$ is the vector of increments of the state variables, $\mathbf J (\mathbf x)\in \mathbb {R}^{k \times s}$ is the Jacobian matrix of measurement functions $\mathbf h (\mathbf x)$ at $\mathbf x = \mathbf x^{(\nu)}$, $\bm \Sigma \in \mathbb {R}^{k \times k}$ is a measurement error covariance matrix, and $\mathbf r (\mathbf x) = \mathbf{z} - \mathbf h (\mathbf x)$ is the vector of residuals [13, Ch. 10]. It is worth noting that assuming uncorrelated measurement errors leads to a diagonal covariance matrix $\bm \Sigma$ corresponding to measurement variances. However, when incorporating PMUs in a rectangular coordinate system and aiming to observe error correlation, this matrix loses its diagonal form.
The nonlinear or AC state estimation represents a non-convex problem arising from the nonlinear measurement functions [22]. Due to the fact that the values of state variables usually fluctuate in narrow boundaries, the model represents the mildly nonlinear problem, where solutions are in a reasonable-sized neighborhood which enables the use of the Gauss-Newton method. The Gauss-Newton method can produce different rates of convergence, which can be anywhere from linear to quadratic [23, Sec. 9.2]. The convergence rate in regard to power system state estimation depends on the topology and measurements, and if parameters are consistent (e.g., free bad data measurement set), the method shows near quadratic convergence rate [13, Sec. 11.2].
Initialization
Let us begin by setting up a new set of measurements for the defined power system:
device = measurement()
@wattmeter(label = "Watmeter ?")
addWattmeter!(system, device; bus = 3, active = -0.5, variance = 1e-3)
addWattmeter!(system, device; from = 1, active = 0.2, variance = 1e-4)
@varmeter(label = "Varmeter ?")
addVarmeter!(system, device; bus = 2, reactive = -0.3, variance = 1e-3)
addVarmeter!(system, device; from = 1, reactive = 0.2, variance = 1e-4)
@pmu(label = "PMU ?")
addPmu!(system, device; bus = 1, magnitude = 1.0, angle = 0, polar = true)
addPmu!(system, device; bus = 3, magnitude = 0.9, angle = -0.2)To compute the voltage magnitudes and angles of buses using the Gauss-Newton method in JuliaGrid, we need to first execute the acModel! function to set up the system. Then, initialize the Gauss-Newton method using the gaussNewton function. The following code snippet demonstrates this process:
acModel!(system)
analysis = gaussNewton(system, device)Initially, the gaussNewton function constructs the mean vector holding measurement values:
julia> 𝐳 = analysis.method.mean8-element Vector{Float64}: -0.5 0.2 -0.3 0.2 1.0 0.0 0.8820599200571175 -0.1788023977155551
Additionally, it forms the precision or weighting matrix denoted as $\mathbf W = \bm \Sigma^{-1}$. We can access these values using the following command:
julia> 𝐖 = analysis.method.precision8×8 SparseArrays.SparseMatrixCSC{Float64, Int64} with 8 stored entries: 1000.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 10000.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 1000.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 10000.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 1.0e8 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 1.0e8 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 1.00756e8 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 1.22324e8
Finally, using initial bus voltage magnitudes and angles from the PowerSystem type, the function creates the initial vector $\mathbf{x}^{(0)}$ of bus voltage magnitudes $\mathbf{V}^{(0)}$ and angles $\bm{\Theta}^{(0)}$ for the Gauss-Newton method:
julia> 𝐕⁽⁰⁾ = analysis.voltage.magnitude3-element Vector{Float64}: 1.0 1.0 1.0julia> 𝚯⁽⁰⁾ = analysis.voltage.angle3-element Vector{Float64}: 0.0 0.0 0.0
Here, we utilize a "flat start" approach in our method. It is important to keep in mind that when dealing with initial conditions in this manner, the Gauss-Newton method may encounter difficulties.
Iterative Process
To apply the Gauss-Newton method, JuliaGrid provides the solve! function. This function is utilized iteratively until a stopping criterion is met, as demonstrated in the following code snippet:
for iteration = 1:20
stopping = solve!(system, analysis)
if stopping < 1e-8
break
end
endThe function solve! calculates the vector of residuals at each iteration using the equation:
\[ \mathbf r (\mathbf x^{(\nu)}) = \mathbf{z} - \mathbf h (\mathbf x^{(\nu)}).\]
The resulting vector from these calculations is stored in the residual variable of the ACStateEstimation type and can be accessed through the following line of code:
julia> 𝐫 = analysis.method.residual8-element Vector{Float64}: 0.0035064869296839163 -0.00195862748866385 0.01806748801610575 0.00552270919867498 -6.954572575601503e-7 0.0 7.228497772571174e-7 -3.090376066161582e-7
The order of the residual vector follows a specific pattern. If all device types exist, the first $|\mathcal V|$ elements correspond to voltmeters, followed by $|\mathcal I|$ elements corresponding to ammeters. Then we have $|\mathcal P|$ elements for wattmeters and $|\mathcal Q|$ elements for varmeters. Finally, we have $2|\bar{\mathcal P}|$ elements for PMUs. The order of these elements within specific devices follows the same order as they appear in the input data defined by the Measurement type.
At the same time, the function forms the Jacobian matrix $\mathbf J (\mathbf x^{(\nu)})$ and calculates the gain matrix $\mathbf G (\mathbf{x}^{(\nu)})$ using:
\[ \mathbf G (\mathbf x^{(\nu)}) = \mathbf J (\mathbf x^{(\nu)})^{T} \bm \Sigma^{-1} \mathbf J (\mathbf x^{(\nu)})\]
The Jacobian matrix and factorized gain matrix are stored in the ACStateEstimation type and can be accessed after each iteration:
julia> 𝐉 = analysis.method.jacobian8×6 SparseArrays.SparseMatrixCSC{Float64, Int64} with 26 stored entries: -1.25177 -3.8648 5.11657 -0.291196 -0.738569 -0.0772033 1.45034 -1.45034 ⋅ 0.257452 0.167366 ⋅ 0.242722 -0.374494 0.131773 -1.43737 5.39071 -4.3514 0.146465 -0.146465 ⋅ 1.8393 -1.65731 ⋅ ⋅ ⋅ ⋅ 1.0 ⋅ ⋅ 1.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.178802 ⋅ ⋅ 0.980067 ⋅ ⋅ 0.882059 ⋅ ⋅ -0.198669julia> 𝐋 = analysis.method.factorization.L6×6 SparseArrays.SparseMatrixCSC{Float64, Int64} with 16 stored entries: 1.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 1.0 ⋅ ⋅ ⋅ ⋅ ⋅ -0.00033994 1.0 ⋅ ⋅ ⋅ ⋅ -8.34403e-5 -4.44128e-5 1.0 ⋅ ⋅ ⋅ 0.0119829 -0.0222009 -0.306728 1.0 ⋅ ⋅ 3.20528e-5 -0.0372364 6.04031e-5 -0.000476186 1.0julia> 𝐔 = analysis.method.factorization.U6×6 SparseArrays.SparseMatrixCSC{Float64, Int64} with 16 stored entries: 1.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.57047 -0.311314 -0.0748163 0.0131228 0.0302766 ⋅ ⋅ 0.96267 -4.18609e-5 -2.55572e-5 -0.0369733 ⋅ ⋅ ⋅ 0.999491 -0.000374433 6.36002e-5 ⋅ ⋅ ⋅ ⋅ 0.468897 -0.192589 ⋅ ⋅ ⋅ ⋅ ⋅ 0.962374
Then finally, the function computes the vector of state variable increments using the equation:
\[ \mathbf \Delta \mathbf x^{(\nu)} = \mathbf G (\mathbf x^{(\nu)})^{-1} \mathbf J (\mathbf x^{(\nu)})^{T} \bm \Sigma^{-1} \mathbf r (\mathbf x^{(\nu)})\]
By default, JuliaGrid uses LU factorization as the primary method for factorizing the gain matrix $\mathbf{G} = \mathbf{L}\mathbf{U}$, aiming to compute the increments. Nevertheless, users have the flexibility to opt for QR or LDLt factorization as an alternative method.
Increment values are stored in the ACStateEstimation type and can be accessed after each iteration:
julia> 𝚫𝐱 = analysis.method.increment6-element Vector{Float64}: 0.0 -2.2536116143951154e-9 -1.1813208586384162e-14 -1.5569564317005995e-13 9.603994185165705e-10 1.126924248894049e-13
Here again, the JuliaGrid implementation of the AC state estimation follows a specific order to store the increment vector and Jacobian matrix. The vector of increments first contains $n$ increments of bus voltage angles $\mathbf \Delta \bm{\Theta}$, followed by $n$ increments of bus voltage magnitudes $\mathbf \Delta \mathbf{V}$. This order also corresponds to the columns of the Jacobian matrix, while the order of rows of the Jacobian is defined according to the order of the residual vector.
Note that the increment vector and Jacobian matrix hold the slack bus with a known voltage angle. An element of the increment vector and a column of the Jacobian matrix are not deleted, and the presence on the slack bus is handled internally by JuliaGrid, which is evident from the factorization of the gain matrix.
The function solve! adds the computed increment term to the previous solution to obtain a new solution:
\[ \mathbf {x}^{(\nu + 1)} = \mathbf {x}^{(\nu)} + \mathbf \Delta \mathbf {x}^{(\nu)}.\]
Therefore, the bus voltage magnitudes $\mathbf{V} = [V_i]$ and angles $\bm{\Theta} = [\theta_i]$ are stored in the following vectors:
julia> 𝐕 = analysis.voltage.magnitude3-element Vector{Float64}: 1.000000695457102 0.875116305093976 0.8999992301629248julia> 𝚯 = analysis.voltage.angle3-element Vector{Float64}: 0.0 -0.13396608670042887 -0.19999982303391817
Finally, the solve! function provides the maximum absolute value of state variable increments, typically employed as the termination criterion for the iteration loop. Specifically, if it falls below a predefined stopping criterion $\epsilon$, the algorithm converges:
\[ \max \{|\Delta x_i|,\; \forall i \} < \epsilon.\]
Jacobian Matrix
As a reminder, the Jacobian matrix consists of $n$ columns representing bus voltage angles $\bm \Theta$, followed by $n$ columns representing bus voltage magnitudes $\mathbf V$. The arrangement of rows is structured such that the first $|\mathcal V|$ rows correspond to voltmeters, followed by $|\mathcal I|$ rows corresponding to ammeters. Then, we have $|\mathcal P|$ rows for wattmeters and $|\mathcal Q|$ rows for varmeters. Finally, there are $2|\bar{\mathcal P}|$ rows for PMUs. The elements are computed based on the provided Jacobian expressions.
Precision Matrix
Let us revisit the precision matrix $\mathbf W$. In the previous example, we introduced a PMU in rectangular coordinates without considering correlations between measurement errors. Now, let us update that PMU to include correlation between measurement errors:
updatePmu!(system, device, analysis; label = "PMU 2", correlated = true)Subsequently, we can examine the updated precision matrix $\mathbf W$:
julia> 𝐖 = analysis.method.precision8×8 SparseArrays.SparseMatrixCSC{Float64, Int64} with 10 stored entries: 1000.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 10000.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 1000.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 10000.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 1.0e8 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 1.0e8 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 1.00926e8 4.56725e6 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 4.56725e6 1.22531e8
Observing the precision matrix, we notice that it loses its diagonal form due to the inclusion of measurement covariances in the 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 method [5, Sec. 3.2].
This approach is suitable when measurement errors are uncorrelated, and the precision matrix remains diagonal. Therefore, as a preliminary step, we need to eliminate the correlation, as we did previously:
updatePmu!(system, device; label = "PMU 2", correlated = false)To address ill-conditioned situations arising from significant differences in measurement variances, users can now employ the orthogonal factorization approach:
analysis = gaussNewton(system, device, Orthogonal)
for iteration = 1:20
stopping = solve!(system, analysis)
if stopping < 1e-8
break
end
endTo explain the method, we begin with the WLS equation:
\[ \Big[\mathbf J (\mathbf x^{(\nu)})^T \mathbf W \mathbf J (\mathbf x^{(\nu)})\Big] \mathbf \Delta \mathbf x^{(\nu)} = \mathbf J (\mathbf x^{(\nu)})^T \mathbf W \mathbf r (\mathbf x^{(\nu)})\]
where $\mathbf W = \bm \Sigma^{-1}$. Subsequently, we can write:
\[ \left[{\mathbf W^{1/2}} \mathbf J (\mathbf x^{(\nu)})\right]^T {\mathbf W^{1/2}} \mathbf J (\mathbf x^{(\nu)}) \Delta \mathbf x^{(\nu)} = \left[{\mathbf W^{1/2}} \mathbf J (\mathbf x^{(\nu)})\right]^T {\mathbf W^{1/2}} \mathbf r (\mathbf x^{(\nu)}).\]
Consequently, we have:
\[ \bar{\mathbf J} (\mathbf x^{(\nu)})^T \bar{\mathbf J}(\mathbf x^{(\nu)}) \Delta \mathbf x^{(\nu)} = \bar{\mathbf J}(\mathbf x^{(\nu)})^{T} \bar{\mathbf r} (\mathbf x^{(\nu)}),\]
where:
\[ \bar{\mathbf J}(\mathbf x^{(\nu)}) = {\mathbf W^{1/2}} \mathbf J (\mathbf x^{(\nu)}), \;\;\; \bar{\mathbf r} (\mathbf x^{(\nu)}) = {\mathbf W^{1/2}} \mathbf r (\mathbf x^{(\nu)}).\]
Therefore, within each iteration of the Gauss-Newton method, JuliaGrid conducts QR factorization on the rectangular matrix:
\[ \bar{\mathbf J}(\mathbf x^{(\nu)}) = {\mathbf W^{1/2}} \mathbf J (\mathbf x^{(\nu)}) = \mathbf Q \mathbf R.\]
Access to the factorized matrix is possible through:
julia> 𝐐 = analysis.method.factorization.Q8×8 SparseArrays.SPQR.QRSparseQ{Float64, Int64}julia> 𝐑 = analysis.method.factorization.R6×6 SparseArrays.SparseMatrixCSC{Float64, Int64} with 15 stored entries: -9920.64 381.086 1.99825 0.169277 0.309314 ⋅ ⋅ 10073.8 0.115791 0.616705 -2.33458 ⋅ ⋅ ⋅ -190.584 24.9998 -4.38277 ⋅ ⋅ ⋅ ⋅ -10001.8 3.74677 ⋅ ⋅ ⋅ ⋅ ⋅ -239.402 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
To obtain the solution, JuliaGrid avoids explicitly forming the orthogonal matrix $\mathbf Q$. Once the algorithm converges, estimates of bus voltage magnitudes $\hat{\mathbf V} = [\hat{V}_i]$ and angles $\hat{\bm \Theta} = [\hat{\theta}_i]$, where $i \in \mathcal N$ can be accessed using variables:
julia> 𝐕 = analysis.voltage.magnitude3-element Vector{Float64}: 1.000000695457102 0.875116305093976 0.8999992301629248julia> 𝚯 = analysis.voltage.angle3-element Vector{Float64}: 0.0 -0.13396608670042887 -0.19999982303391817
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:
addWattmeter!(system, device; bus = 3, active = 5.1, variance = 1e-3)Subsequently, we will construct the AC state estimation model and solve it:
analysis = gaussNewton(system, device)
for iteration = 1:20
stopping = solve!(system, analysis)
if stopping < 1e-8
break
end
endNow, 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 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 J (\hat {\mathbf x}) [\mathbf J (\hat {\mathbf x})^T \bm \Sigma^{-1} \mathbf J (\hat {\mathbf x})]^{-1} \mathbf J (\hat {\mathbf x})^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.maxNormalizedResidual148.26877093119242
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"Watmeter 3"
is removed from the measurement set and marked as out-of-service.
Subsequently, we can solve the system again, but this time without the removed measurement:
analysis = gaussNewton(system, device)
for iteration = 1:20
stopping = solve!(system, analysis)
if stopping < 1e-8
break
end
endFollowing that, we check for outliers once more:
outlier = residualTest!(system, device, analysis; threshold = 4.0)To examine the value:
julia> outlier.maxNormalizedResidual0.8182344209553056
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.
This section outlines the method as described in [5, Sec. 6.5]. Hence, we consider the system of nonlinear equations:
\[ \mathbf{z}=\mathbf{h}(\mathbf x)+\mathbf{u}.\]
Subsequently, the LAV state estimator is derived as the solution to the optimization problem:
\[ \begin{aligned} \text{minimize}& \;\;\; \sum_{i \in \mathcal M} |r_i|\\ \text{subject\;to}& \;\;\; z_i - h_i(\mathbf x) = r_i, \;\;\; \forall i \in \mathcal M, \end{aligned}\]
where $r_i$ represents the $i$-th measurement residual. Let $\eta_i$ be defined in a manner that ensures:
\[ |r_i| \leq \eta_i,\]
and replace the above inequality with two equalities using two non-negative slack variables $q_i$ and $w_i$:
\[ \begin{aligned} r_i - q_i &= -\eta_i \\ r_i + w_i &= \eta_i. \end{aligned}\]
Let us now define two additional non-negative variables $\overline{r}_i$ and $\underline{r}_i$, which satisfy the following relationships:
\[ \begin{aligned} r_i &= \overline{r}_i - \underline{r}_i\\ \overline{r}_i &= \cfrac{1}{2} q_i \\ \underline{r}_i &= \cfrac{1}{2} w_i. \end{aligned}\]
Then, the above two equalities become:
\[ \begin{aligned} r_i - 2\overline{r}_i &= -2 \eta_i \\ r_i + 2 \underline{r}_i &= 2 \eta_i, \end{aligned}\]
that is:
\[ \begin{aligned} \overline{r}_i + \underline{r}_i = \eta_i. \end{aligned}\]
Next, we define a vector of state variables according to two additional nonnegative vectors $\overline{\mathbf x} \in \mathbb {R}_{\ge 0}^s$ and $\underline{\mathbf x} \in \mathbb {R}_{\ge 0}^s$, which satisfy the following relationships:
\[ \mathbf x = \overline{\mathbf x} - \underline{\mathbf x}\]
Hence, the optimization problem can be written:
\[ \begin{aligned} \text{minimize}& \;\;\; \sum_{i \in \mathcal M} (\overline{r}_i + \underline{r}_i) \\ \text{subject\;to} & \;\;\; h_i(\overline{\mathbf x} - \underline{\mathbf x}) + \overline{r}_i - \underline{r}_i = z_i, \;\;\; \forall i \in \mathcal M \\ & \;\;\; \overline{r}_i \geq 0, \; \underline{r}_i \geq 0, \;\;\; \forall i \in \mathcal M \\ & \;\;\; \overline{\mathbf x} \succeq \mathbf 0, \; \underline{\mathbf x} \succeq \mathbf 0. \end{aligned}\]
To form the above optimization problem, the user can call the following function:
using Ipopt
analysis = acLavStateEstimation(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 non-negative variables, while the state estimator is obtained by:
\[ \hat{\mathbf x} = \overline{\mathbf x} - \underline{\mathbf x}.\]
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}: 1.0000000064101187 0.8716991294249965 0.8915402891580015julia> 𝚯 = analysis.voltage.angle3-element Vector{Float64}: 0.0 -0.133006380063307 -0.19921333228739968
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} = [Q_{\mathrm{s}ij}]$ |
| Branch | Series elements | $\mathbf P_\mathrm{l} = [P_{\mathrm{l}ij}]$ | $\mathbf Q_\mathrm{l} = [Q_{\mathrm{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.4569858239661292 0.048580814062877356 -0.5000000002429149julia> 𝐐 = analysis.power.injection.reactive3-element Vector{Float64}: 0.35278912125044276 -0.3000000005296262 -0.016001985804536115
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_{\mathrm{p}i} &= P_i + P_{\mathrm{d}i}\\ Q_{\mathrm{p}i} &= Q_i + Q_{\mathrm{d}i}. \end{aligned}\]
The active and reactive power injections from the generators at each bus are stored as vectors, denoted by $\mathbf{P}_\mathrm{p} = [P_{\mathrm{p}i}]$ and $\mathbf{Q}_\mathrm{p} = [Q_{\mathrm{p}i}]$, which can be obtained using:
julia> 𝐏ₚ = analysis.power.supply.active3-element Vector{Float64}: 0.9569858239661292 0.048580814062877356 -2.429149104088424e-10julia> 𝐐ₚ = analysis.power.supply.reactive3-element Vector{Float64}: 0.35278912125044276 -5.296262317600053e-10 -0.016001985804536115
Power at Bus Shunt Elements
Active and reactive powers associated with the shunt elements at each bus are represented by the vectors $\mathbf{P}_\mathrm{sh} = [{P}_{\mathrm{sh}i}]$ and $\mathbf{Q}_\mathrm{sh} = [{Q}_{\mathrm{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}_\mathrm{i} = [P_{ij}]$ and $\mathbf{Q}_\mathrm{i} = [Q_{ij}],$ respectively, and can be retrieved using the following commands:
julia> 𝐏ᵢ = analysis.power.from.active3-element Vector{Float64}: 0.20000000151719766 0.25698582244893153 0.24681281558781917julia> 𝐐ᵢ = analysis.power.from.reactive3-element Vector{Float64}: 0.20000000006001498 0.1527891211904276 -0.11784281253607717
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.19823200152494175 -0.2550678666049177 -0.24493213363799723julia> 𝐐ⱼ = analysis.power.to.reactive3-element Vector{Float64}: -0.18215718799354885 -0.12155754865018834 0.10555556284565207
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}_\mathrm{s} = [P_{\mathrm{s}ij}]$ and $\mathbf{Q}_\mathrm{s} = [Q_{\mathrm{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.035197187701210685 -0.03589688200024341 -0.031094069188644595
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}_\mathrm{l} = [P_{\mathrm{l}ij}]$ and $\mathbf{Q}_\mathrm{l} = [Q_{\mathrm{l}ij}]$. We can retrieve these values using the following code:
julia> 𝐏ₗ = analysis.power.series.active3-element Vector{Float64}: 0.0017679999922558948 0.0019179558440137988 0.0018806819498218604julia> 𝐐ₗ = analysis.power.series.reactive3-element Vector{Float64}: 0.05303999976767691 0.0671284545404826 0.018806819498218604
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}]$ |
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.577318112573757 0.3486386852476296 0.5611142921313774julia> 𝛙 = analysis.current.injection.angle3-element Vector{Float64}: -0.6574276811383212 1.2772475358372901 2.91038626973281
Current Flows
To obtain the vectors of magnitudes $\mathbf{I}_\mathrm{i} = [I_{ij}]$ and angles $\bm{\psi}_\mathrm{i} = [\psi_{ij}]$ for the resulting complex current flows, you can use the following commands:
julia> 𝐈ᵢ = analysis.current.from.magnitude3-element Vector{Float64}: 0.2828427117768214 0.2989752955478835 0.31375765527495864julia> 𝛙ᵢ = analysis.current.from.angle3-element Vector{Float64}: -0.7853981597544916 -0.5363973366973052 0.3124457478414518
Similarly, we can obtain the vectors of magnitudes $\mathbf{I}_mathrm{j} = [I_{ji}]$ and angles $\bm{\psi}_\mathrm{j} = [\psi_{ji}]$ of the resulting complex current flows using the following code:
julia> 𝐈ⱼ = analysis.current.to.magnitude3-element Vector{Float64}: 0.3088403202038889 0.3169261057594029 0.2991553732467297julia> 𝛙ⱼ = analysis.current.to.angle3-element Vector{Float64}: 2.2654218416675307 2.497651247230784 -2.933899355980511
Current at Branch Series Elements
To obtain the vectors of magnitudes $\mathbf{I}_\mathrm{l} = [I_{\mathrm{l}ij}]$ and angles $\bm{\psi}_\mathrm{l} = [\psi_{\mathrm{l}ij}]$ of the resulting complex current flows, one can use the following code:
julia> 𝐈ₗ = analysis.current.series.magnitude3-element Vector{Float64}: 0.29732137429521416 0.3096736866456196 0.30664979616998445julia> 𝛙ₗ = analysis.current.series.angle3-element Vector{Float64}: -0.8329812636144827 -0.591939510081054 0.2611180954469363