DC State Estimation

To initiate the process, let us construct the PowerSystem 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:

monitoring = measurement(system)

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!(monitoring; 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!(monitoring; 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!(monitoring; 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!(
  monitoring; label = "V₁, θ₁", bus = 1, magnitude = 1.0, angle = 0,
  varianceMagnitude = 1e-5, 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 = dcStateEstimation(monitoring)

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!(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 [5, Sec. 3.2].

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

analysis = dcStateEstimation(monitoring, 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!(analysis)

Access to the factorized matrix is possible through:

julia> 𝐐 = analysis.method.factorization.Q4×4 SparseArrays.SPQR.QRSparseQ{Float64, Int64}
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

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 analysis, as discussed in [5, Ch. 6]. It is important to note that robustness often comes at the cost of increased computational complexity.

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

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

The LAV state estimator is then formulated as the solution to the following optimization problem:

\[ \begin{aligned} \text{minimize}& \;\;\; \sum_{i \in \mathcal M} |r_i|\\ \text{subject\;to}& \;\;\; z_i - h_i(\bm {\Theta}) = r_i, \;\;\; \forall i \in \mathcal M, \end{aligned}\]

where $r_i$ denotes the residual of the $i$-th measurement.

To explicitly handle absolute values, we introduce two nonnegative variables $u_i \ge 0$ and $v_i \ge 0$, referred to as positive and negative deviations. This allows the optimization problem to be rewritten as:

\[ \begin{aligned} \text{minimize}& \;\;\; \sum_{i \in \mathcal M} (u_i + v_i) \\ \text{subject\;to} & \;\;\; z_i - h_i(\bm {\Theta}) = u_i - v_i, \;\;\; \forall i \in \mathcal M \\ & \;\;\; u_i \geq 0, \; v_i \geq 0, \;\;\; \forall i \in \mathcal M. \end{aligned}\]

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

using Ipopt

analysis = dcLavStateEstimation(monitoring, Ipopt.Optimizer)

Then the user can solve the optimization problem by:

solve!(analysis)

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.05600000000000001
 -0.1115

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!(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.395
 -0.09500000000000008
 -1.3

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 into 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.395
 0.004999999999999921
 0.0

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.28
 1.115
 0.185

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.28
 -1.115
 -0.185