DC State Estimation

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

system = powerSystem()

addBus!(system; label = 1, type = 3, angle = 0.0)
addBus!(system; label = 2, type = 1, active = 0.1)
addBus!(system; label = 3, type = 1, active = 1.3)
addBranch!(system; label = 1, from = 1, to = 2, reactance = 0.2)
addBranch!(system; label = 2, from = 1, to = 3, reactance = 0.1)
addBranch!(system; label = 3, from = 2, to = 3, reactance = 0.3)
addGenerator!(system; label = 1, bus = 1, active = 3.2)

dcModel!(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}$. In typical scenarios, the DC state estimation model relies solely on active power measurements originating from the set of wattmeters $\mathcal{P}$. However, we provide the option for users to include measurements from the set of PMUs $\bar{\mathcal{P}}$. Specifically, we utilize only the PMUs installed at the buses $\bar{\mathcal{P}}_\text{b} \subset \bar{\mathcal{P}}$ that measure bus voltage angles. 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()

State Estimation Model

In accordance with the DC Model, the DC state estimation is derived through the linearization of the non-linear model. In this linearized model, all bus voltage magnitudes are assumed to be $V_i \approx 1$, $i \in \mathcal{N}$. Additionally, shunt elements and branch resistances are neglected. This simplification implies that the DC model disregards reactive powers and transmission losses, focusing solely on active powers. Consequently, the DC state estimation considers only bus voltage angles, represented as $\mathbf x \equiv \bm {\Theta}$, as the state variables. As a result, the total number of state variables is $n-1$, with one voltage angle corresponding to the slack bus.

Within the JuliaGrid framework for DC state estimation, the methodology encompasses both active power flow and injection measurements from the set $\mathcal{P}$, along with bus voltage angle measurements represented by the set $\bar{\mathcal{P}}_\text{b}$. These measurements contribute to the construction of a linear system of equations:

\[ \mathbf{z}=\mathbf{h}(\bm {\Theta})+\mathbf{u},\]

where $\mathbf{h}(\bm {\Theta})=$ $[h_1(\bm {\Theta})$, $\dots$, $h_k(\bm {\Theta})]^{{T}}$ is the vector of linear measurement functions, $\mathbf{z} = [z_1,\dots,z_k]^{\mathrm{T}}$ is the vector of measurement values, and $\mathbf{u} = [u_1,\dots,u_k]^{\mathrm{T}}$ is the vector of uncorrelated measurement errors, and this defines the vector of measurement variances $\mathbf{v} = [v_1,\dots,v_k]^{\mathrm{T}}$, where $k = |\mathcal{M}|$.

Therefore, the linear system of equations can be represented based on the specific devices from which measurements originate, whether wattmeters or PMUs:

\[ \begin{bmatrix} \mathbf{z}_\mathcal{P}\\[3pt] \mathbf{z}_{\bar{\mathcal{P}}_\text{b}} \end{bmatrix} = \begin{bmatrix} \mathbf{h}_\mathcal{P}(\bm {\Theta})\\[3pt] \mathbf{h}_{\bar{\mathcal{P}}_\text{b}}(\bm {\Theta}) \end{bmatrix} + \begin{bmatrix} \mathbf{u}_\mathcal{P}\\[3pt] \mathbf{u}_{\bar{\mathcal{P}}_\text{b}}. \end{bmatrix}\]

In summary, upon user definition of the measurement devices, each $i$-th measurement device is linked to the measurement function $h_i(\bm {\Theta})$, the corresponding measurement value $z_i$, and the measurement variance $v_i$.

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}(\bm {\Theta}) = [h_{P_{i}}(\bm {\Theta})].\]

For example:

addWattmeter!(system, device; label = "P₃", bus = 3, active = -1.30, variance = 1e-3)

Here, utilizing the DC Model, we derive the function defining the active power injection as follows:

\[ h_{P_{i}}(\bm {\Theta}) = B_{ii}\theta_i + \sum_{j \in \mathcal{N}_i \setminus i} {B}_{ij} \theta_j + P_{\text{tr}i} + P_{\text{sh}i},\]

where $\mathcal{N}_i \setminus i$ contains buses incident to bus $i$, excluding bus $i$, with the following coefficient expressions:

\[ \begin{aligned} \cfrac{\mathrm \partial{h_{P_{i}}(\bm {\Theta})}}{\mathrm \partial \theta_{i}} = B_{ii}, \;\;\; \cfrac{\mathrm \partial{{h_{P_{i}}}(\bm {\Theta})}}{\mathrm \partial \theta_{j}} = {B}_{ij}. \end{aligned}\]

From-Bus End 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}(\bm {\Theta}) = [h_{P_{ij}}(\bm {\Theta})].\]

For example:

addWattmeter!(system, device; label = "P₁₂", from = 1, active = 0.28, variance = 1e-4)

Here, the function describing active power flow at the from-bus end is defined as follows:

\[ h_{P_{ij}}(\bm {\Theta}) = \cfrac{1}{\tau_{ij} x_{ij}} (\theta_{i} -\theta_{j}-\phi_{ij}),\]

with the following coefficient expressions:

\[ \begin{aligned} \cfrac{\mathrm \partial{h_{P_{ij}}(\bm {\Theta})}}{\mathrm \partial \theta_{i}} = \cfrac{1}{\tau_{ij} x_{ij}}, \;\;\; \cfrac{\mathrm \partial{{h_{P_{ij}}}(\bm {\Theta})}}{\mathrm \partial \theta_{j}} = -\cfrac{1}{\tau_{ij} x_{ij}}. \end{aligned}\]

To-Bus End 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}(\bm {\Theta}) = [h_{P_{ji}}(\bm {\Theta})].\]

For example:

addWattmeter!(system, device; label = "P₂₁", to = 1, active = -0.28, variance = 1e-4)

Thus, the function describing active power flow at the to-bus end is defined as follows:

\[ h_{P_{ji}}(\bm {\Theta}) = -\cfrac{1}{\tau_{ij} x_{ij}} (\theta_{i} -\theta_{j}-\phi_{ij}),\]

with the following coefficient expressions:

\[ \cfrac{\mathrm \partial{h_{P_{ji}}(\bm {\Theta})}}{\mathrm \partial \theta_{i}} = -\cfrac{1}{\tau_{ij} x_{ij}}, \;\;\; \cfrac{\mathrm \partial{{h_{P_{ji}}}(\bm {\Theta})}}{\mathrm \partial \theta_{j}} = \cfrac{1}{\tau_{ij} x_{ij}}.\]

Bus Voltage Angle Measurements

If the user opts to include phasor measurements that measure bus voltage angle at bus $i \in \mathcal{N}$, denoted as $\theta_i \in \bar{\mathcal{P}}_\text{b}$, the user will specify the measurement values, variances, and measurement functions of vectors:

\[ \mathbf{z}_{\bar{\mathcal{P}}_\text{b}} = [z_{\theta_i}], \;\;\; \mathbf{v}_{\bar{\mathcal{P}}_\text{b}} = [v_{\theta_i}], \;\;\; \mathbf{h}_{\bar{\mathcal{P}}_\text{b}}(\bm {\Theta}) = [h_{\theta_{i}}(\bm {\Theta})].\]

For example:

addPmu!(system, device; label = "V₁, θ₁", bus = 1, magnitude = 1.0, varianceMagnitude = 1e-5,
angle = 0, varianceAngle = 1e-6)

Here, the function defining the bus voltage angle measurement is straightforward:

\[ h_{\theta_{i}}(\bm {\Theta}) = \theta_{i},\]

with the following coefficient expression:

\[ \cfrac{\mathrm \partial{{h_{\theta_i}(\bm {\Theta})}}}{\mathrm \partial \theta_{i}}=1.\]

Weighted Least-squares Estimation

The solution to the DC 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 \bm {\Theta} = \mathbf H^{T} \bm \Sigma^{-1} (\mathbf z - \mathbf{c}).\]

Here, $\mathbf z \in \mathbb {R}^{k}$ denotes the vector of measurement values, the vector $\mathbf c \in \mathbb {R}^{k}$ holds constant terms, $\mathbf {H} \in \mathbb {R}^{k \times (n-1)}$ represents the coefficient matrix, and $\bm \Sigma \in \mathbb {R}^{k \times k}$ is the measurement error covariance matrix, where the diagonal elements hold measurement variances.

The inclusion of the vector $\mathbf{c}$ is necessary due to the fact that measurement functions associated with active power measurements may include constant terms, especially when there are non-zero shift angles of transformers or shunt elements in the system consuming active powers, as evident from the provided measurement functions.

Implementation

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

analysis = dcWlsStateEstimation(system, device)

Coefficient Matrix

Using the above-described equations, JuliaGrid forms the coefficient matrix $\mathbf{H} \in \mathbb{R}^{k \times (n-1)}$:

julia> 𝐇 = analysis.method.coefficient4×3 SparseArrays.SparseMatrixCSC{Float64, Int64} with 8 stored entries:
 -10.0  -3.33333  13.3333
   5.0  -5.0        ⋅
  -5.0   5.0        ⋅
   1.0    ⋅         ⋅

Each row in the matrix corresponds to a specific measurement. The first $|\mathcal{P}|$ rows correspond to wattmeters, ordered as users add wattmeters, while the last $|{\bar{\mathcal{P}}_\text{b}}|$ rows correspond to PMUs, also in the order users add PMUs.

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.precision4×4 SparseArrays.SparseMatrixCSC{Float64, Int64} with 4 stored entries:
 1000.0       ⋅        ⋅    ⋅
     ⋅   10000.0       ⋅    ⋅
     ⋅        ⋅   10000.0   ⋅
     ⋅        ⋅        ⋅   1.0e6

Mean Vector

Users can access the vector $\mathbf z - \mathbf{c}$, which contains the means of Gaussian distributions describing each measurement, using the following command:

julia> 𝐳 = analysis.method.mean4-element Vector{Float64}:
 -1.3
  0.28
 -0.28
  0.0

In the context of the power system, where phase-shifting transformers and shunt elements consuming active powers are absent, and the slack angle has a zero value, the vector $\mathbf{c}= \mathbf{0}$. Consequently, the vector of means holds values that are equal to the measurement values.

Estimate of State Variables

Once the model is established, we solve the WLS equation to derive the estimate of bus voltage angles:

\[ \hat{\bm {\Theta}} = [\mathbf H^{T} \bm \Sigma^{-1} \mathbf H]^{-1} \mathbf H^{T} \bm \Sigma^{-1} (\mathbf z - \mathbf{c}).\]

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.\]

Access to the factorized gain matrix is available through:

julia> 𝐋 = analysis.method.factorization.L3×3 SparseArrays.SparseMatrixCSC{Float64, Int64} with 4 stored entries:
 1.0    ⋅         ⋅
  ⋅    1.0        ⋅
  ⋅   -0.217391  1.0
julia> 𝐔 = analysis.method.factorization.U3×3 SparseArrays.SparseMatrixCSC{Float64, Int64} with 4 stored entries:
 1.0   ⋅      ⋅
  ⋅   0.92  -0.08
  ⋅    ⋅     0.782609

Finally, the estimated bus voltage angles $\hat{\bm {\Theta}} = [\hat{\theta}_i]$, $i \in \mathcal{N}$, can be retrieved using the variable:

julia> 𝚯 = analysis.voltage.angle3-element Vector{Float64}:
  0.0
 -0.05600000000000001
 -0.11150000000000002

It is essential to note that the slack bus voltage angle is temporarily excluded from the gain matrix $\mathbf G$ during computation. It is important to emphasize that this internal handling does not alter the stored elements.

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 [1, Sec. 3.2].

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

analysis = dcWlsStateEstimation(system, device, Orthogonal)

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

\[ \mathbf H^{T} \mathbf W \mathbf H \bm {\Theta} = \mathbf H^{T} \mathbf W (\mathbf z - \mathbf{c}),\]

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 \bm {\Theta} = \left({\mathbf W^{1/2}} \mathbf H\right)^{T} {\mathbf W^{1/2}} (\mathbf z - \mathbf{c}).\]

Consequently, we have:

\[ \bar{\mathbf{H}}^{T} \bar{\mathbf{H}} \bm {\Theta} = \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 - \mathbf{c}).\]

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.Q4×4 SuiteSparse.SPQR.QRSparseQ{Float64, Int64}:
 1.0   0.0       0.0       0.0
 0.0  -0.707107  0.707107  0.0
 0.0   0.707107  0.707107  0.0
 0.0   0.0       0.0       1.0
julia> 𝐑 = analysis.method.factorization.R3×3 SparseArrays.SparseMatrixCSC{Float64, Int64} with 4 stored entries:
 421.637  -105.409  0.0
    ⋅      707.107   ⋅
    ⋅         ⋅      ⋅

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

julia> 𝚯 = analysis.voltage.angle3-element Vector{Float64}:
  0.0
 -0.05600000000000001
 -0.1115

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 [1, 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 [2].

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

addWattmeter!(system, device; label = "P₁", bus = 1, active = 13.1, variance = 1e-4)

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

analysis = dcWlsStateEstimation(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 [1, 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 {\bm {\Theta}}), \;\;\; 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.maxNormalizedResidual420.45601204468034

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 the variable:

julia> outlier.detecttrue

This indicates that the measurement labeled as:

julia> outlier.label"P₁"

is removed from the DC model and marked as out-of-service.

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.maxNormalizedResidual4.76837158203125e-7

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 [1, 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}(\bm {\Theta})+\mathbf{u}.\]

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}\bm {\Theta} - \mathbf{c} =\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:

\[ \bm {\Theta}_x \in \mathbb {R}_{\ge 0}^{n}; \;\;\; \bm {\Theta}_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:

\[ \bm {\Theta} = \bm {\Theta}_x - \bm {\Theta}_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}(\bm {\Theta}_x - \bm {\Theta}_y) + \mathbf {r}_x - \mathbf {r}_y = \mathbf{z} - \mathbf{c} \\ & \;\;\; \bm {\Theta}_x \succeq \mathbf 0, \; \bm {\Theta}_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 = dcLavStateEstimation(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{\bm {\Theta}} = \bm {\Theta}_x - \bm {\Theta}_y.\]

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

julia> 𝚯 = analysis.voltage.angle3-element Vector{Float64}:
  0.0
 -0.055999999999915506
 -0.1114999999999376

Power Analysis

After obtaining the solution from the DC state estimation, we can calculate powers related to buses and branches using the power! function:

power!(system, analysis)

Power Injections

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

julia> 𝐏 = analysis.power.injection.active3-element Vector{Float64}:
  1.3949999999989533
 -0.09499999999950393
 -1.2999999999994496

Generator Power Injections

We can determine the active power supplied by generators to the buses by summing the active power injections and the active power demanded by consumers at each bus:

\[ P_{\text{p}i} = P_i + P_{\text{d}i},\;\;\; i \in \mathcal{N}.\]

The vector of active power injected by generators to the buses, denoted by $\mathbf{P}_{\text{p}} = [P_{\text{p}i}]$, can be obtained using:

julia> 𝐏ₚ = analysis.power.supply.active3-element Vector{Float64}:
 1.3949999999989533
 0.00500000000049608
 5.504485756091526e-13

Power Flows

The resulting active power flows are stored as the vector $\mathbf{P}_{\text{i}} = [P_{ij}]$, which can be retrieved using:

julia> 𝐏ᵢ = analysis.power.from.active3-element Vector{Float64}:
 0.27999999999957753
 1.1149999999993758
 0.18500000000007363

Similarly, the resulting active power flows are stored as the vector $\mathbf{P}_{\text{j}} = [P_{ji}]$, which can be retrieved using:

julia> 𝐏ⱼ = analysis.power.to.active3-element Vector{Float64}:
 -0.27999999999957753
 -1.1149999999993758
 -0.18500000000007363

References

[1] A. Abur and A. Exposito, Power System State Estimation: Theory and Implementation, Taylor & Francis, 2004.

[2] G. Korres, "A distributed multiarea state estimation," IEEE Trans. Power Syst., vol. 26, no. 1, pp. 73–84, Feb. 2011.

[3] M. Cosovic, M. Delalic, D. Raca, and D. Vukobratovic, "Observability analysis for large-scale power systems using factor graphs," IEEE Trans. Power Syst., 36(5), 4791-4799.

[4] A. Monticelli, State Estimation in Electric Power Systems: A Generalized Approach, ser. Kluwer international series in engineering and computer science. Springer US, 1999.

[5] H. Horisberger, "Observability analysis for power systems with measurement deficiencies," IFAC Proceedings Volumes, vol. 18, no. 7, pp.51–58, 1985.

[6] N. M. Manousakis and G. N. Korres, "Observability analysis for power systems including conventional and phasor measurements," in Proc. MedPower 2010, Agia Napa, 2010, pp. 1-8.