DC State Estimation

To perform the DC state estimation, we first need to have the PowerSystem type that has been created with the DC model, alongside the Measurement type that retains measurement data. Subsequently, we can formulate either the weighted least-squares (WLS) or the least absolute value (LAV) DC state estimation model encapsulated within the type DCStateEstimation using:

• dcStateEstimation,
• dcLavStateEstimation.

For resolving the DC state estimation problem and obtaining bus voltage angles, utilize the following function:

• solve!.

After executing the function solve!, where the user employs the WLS method, the user has the ability to check if the measurement set contains outliers throughout bad data analysis and remove those measurements using:

• residualTest!.

Moreover, before creating the DCStateEstimation type, users can initiate observability analysis to identify observable islands and restore observability by employing:

• islandTopologicalFlow,
• islandTopological,
• restorationGram!.

For more detailed information, users can refer to the Observability Analysis section within AC state estimation documentation. It is worth noting that when the system becomes observable within the AC model, it will also be observable within the DC state estimation model.

After obtaining the solution for DC state estimation, JuliaGrid offers a post-processing analysis function to compute active powers associated with buses and branches:

• power!.

Additionally, specialized functions are available for calculating specific types of powers for individual buses or branches.

Bus Type Modification

Just like in the Bus Type Modification section, when establishing the DCStateEstimation type, the initially assigned slack bus is evaluated and may be altered. If the designated slack bus (type = 3) lacks a connected in-service generator, it will be changed to a demand bus (type = 1). Conversely, the first generator bus (type = 2) with an active in-service generator linked to it will be reassigned as the new slack bus (type = 3).

Weighted Least-Squares Estimator

To solve the DC state estimation and derive WLS estimates using JuliaGrid, the process initiates by defining PowerSystem and Measurement types. Here is an illustrative example:

system = powerSystem()
device = measurement()

addBus!(system; label = "Bus 1", type = 3)
addBus!(system; label = "Bus 2", type = 1, active = 0.2)
addBus!(system; label = "Bus 3", type = 1, active = 0.4)

addBranch!(system; label = "Branch 1", from = "Bus 1", to = "Bus 2", reactance = 0.5)
addBranch!(system; label = "Branch 2", from = "Bus 1", to = "Bus 3", reactance = 0.2)
addBranch!(system; label = "Branch 3", from = "Bus 2", to = "Bus 3", reactance = 0.3)

addGenerator!(system; label = "Generator 1", bus = "Bus 1", active = 3.2)

@wattmeter(label = "Wattmeter ?")
addWattmeter!(system, device; bus = "Bus 1", active = 0.6, variance = 1e-3)
addWattmeter!(system, device; bus = "Bus 3", active = -0.4, variance = 1e-2)
addWattmeter!(system, device; from = "Branch 1", active = 0.18, variance = 1e-4)
addWattmeter!(system, device; to = "Branch 2", active = -0.42, variance = 1e-4)

The dcStateEstimation function serves to establish the DC state estimation problem:

analysis = dcStateEstimation(system, device)

analysis = dcStateEstimation(system, device, LDLt)

To obtain the bus voltage angles, the solve! function can be invoked as shown:

solve!(system, analysis)

Upon obtaining the solution, access the bus voltage angles using:

julia> print(system.bus.label, analysis.voltage.angle)Bus 1: 0.0
Bus 2: -0.09000000000000001
Bus 3: -0.08399999999999999

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

Specifically, by specifying the Orthogonal argument in the dcStateEstimation function, JuliaGrid implements a more robust approach to obtain the WLS estimator, which proves particularly beneficial when substantial differences exist among measurement variances:

analysis = dcStateEstimation(system, device, Orthogonal)
solve!(system, analysis)

Print Results in the REPL

Users have the option to print the results in the REPL using any units that have been configured, such as:

@voltage(pu, deg, V)
printBusData(system, analysis)

|-----------------|
| Bus Data        |
|-----------------|
| Label | Voltage |
|       |         |
|   Bus |   Angle |
|       |   [deg] |
|-------|---------|
| Bus 1 |  0.0000 |
| Bus 2 | -5.1566 |
| Bus 3 | -4.8128 |
|-----------------|

Next, users can easily customize the print results for specific buses, for example:

printBusData(system, analysis; label = "Bus 1", header = true)
printBusData(system, analysis; label = "Bus 2")
printBusData(system, analysis; label = "Bus 3", footer = true)

Save Results to a File

Users can also redirect print output to a file. For example, data can be saved in a text file as follows:

open("bus.txt", "w") do file
    printBusData(system, analysis, file)
end

Bad Data Processing

After acquiring the WLS solution using the solve! function, users can conduct bad data analysis employing the largest normalized residual test. Continuing with our defined power system and measurement set, let us introduce a new wattmeter. Upon proceeding to find the solution for this updated state:

addWattmeter!(system, device; from = "Branch 2", active = 4.1, variance = 1e-4)

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

Following the solution acquisition, we can verify the presence of erroneous data. Detection of such data is determined by the threshold keyword. If the largest normalized residual's value exceeds the threshold, the measurement will be identified as bad data and consequently removed from the DC state estimation model:

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

Users can examine the data obtained from the bad data analysis:

julia> outlier.detecttrue
julia> outlier.maxNormalizedResidual267.65576812169894
julia> outlier.label"Wattmeter 5"

Hence, upon detecting bad data, the detect variable will hold true. The maxNormalizedResidual variable retains the value of the largest normalized residual, while the label contains the label of the measurement identified as bad data. JuliaGrid will mark the respective measurements as out-of-service within the Measurement type.

Moreover, JuliaGrid will adjust the coefficient matrix and mean vector within the DCStateEstimation type based on measurements now designated as out-of-service. To optimize the algorithm's efficiency, JuliaGrid resets non-zero elements to zero in the coefficient matrix and mean vector, effectively removing the impact of the corresponding measurement on the solution:

julia> analysis.method.mean5-element Vector{Float64}:
  0.6
 -0.4
  0.18
 -0.42
  0.0
julia> analysis.method.coefficient5×3 SparseArrays.SparseMatrixCSC{Float64, Int64} with 12 stored entries:
  7.0  -2.0      -5.0
 -5.0  -3.33333   8.33333
  2.0  -2.0        ⋅
 -5.0    ⋅        5.0
  0.0    ⋅        0.0

Hence, after removing bad data, a new estimate can be computed without considering this specific measurement:

solve!(system, analysis)

Least Absolute Value Estimator

The LAV method presents an alternative estimation technique known for its increased robustness compared to WLS. While the WLS method relies on specific assumptions regarding measurement errors, robust estimators like LAV are designed to maintain unbiasedness even in the presence of various types of measurement errors and outliers. This characteristic often eliminates the need for extensive bad data processing procedures [5, Ch. 6]. However, it is important to note that achieving robustness typically involves increased computational complexity.

To obtain an LAV estimator, users need to employ one of the solvers listed in the JuMP documentation. In many common scenarios, the Ipopt solver proves sufficient to obtain a solution:

using Ipopt

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

Setup Starting Primal Values

In JuliaGrid, the assignment of starting primal values for optimization variables takes place when the solve! function is executed. Starting primal values are determined based on the voltage fields within the DCStateEstimation type. By default, these values are initially established using the initial bus voltage angles from PowerSystem type:

julia> print(system.bus.label, analysis.voltage.angle)Bus 1: 0.0
Bus 2: 0.0
Bus 3: 0.0

Users have the flexibility to customize these values according to their requirements, and they will be utilized as the starting primal values when executing the solve! function.

Solution

To solve the formulated LAV state estimation model, simply execute the following function:

solve!(system, analysis)

Upon obtaining the solution, access the bus voltage angles using:

julia> print(system.bus.label, analysis.voltage.angle)Bus 1: 0.0
Bus 2: -0.08999999999956487
Bus 3: -0.08399999999987212

Measurement Set Update

After establishing the Measurement type using the measurement function, users gain the capability to incorporate new measurement devices or update existing ones.

Once updates are completed, users can seamlessly progress towards generating the DCStateEstimation type using the dcStateEstimation or dcLavStateEstimation function. Ultimately, resolving the DC state estimation is achieved through the utilization of the solve! function:

system = powerSystem()
device = measurement() # <- Initialize the Measurement instance

addBus!(system; label = "Bus 1", type = 3)
addBus!(system; label = "Bus 2", type = 1, active = 0.1)

addBranch!(system; label = "Branch 1", from = "Bus 1", to = "Bus 2", reactance = 0.5)

addGenerator!(system; label = "Generator 1", bus = "Bus 1", active = 0.1)

@wattmeter(label = "Wattmeter ?")
addWattmeter!(system, device; bus = "Bus 2", active = -0.11, variance = 1e-3)
addWattmeter!(system, device; from = "Branch 1", active = 0.09, variance = 1e-4)

analysis = dcStateEstimation(system, device) # <- Build DCStateEstimation for the model
solve!(system, analysis)

addWattmeter!(system, device; to = "Branch 1", active = -0.12, variance = 1e-4)
updateWattmeter!(system, device; label = "Wattmeter 1", status = 0)
updateWattmeter!(system, device; label = "Wattmeter 2", active = 0.1, noise = false)

analysis = dcStateEstimation(system, device) # <- Build DCStateEstimation for new model
solve!(system, analysis)

State Estimation Update

An advanced methodology involves users establishing the DCStateEstimation type using dcStateEstimation or dcLavStateEstimation just once. After this initial setup, users can seamlessly modify existing measurement devices without the need to recreate the DCStateEstimation type.

This advancement extends beyond the previous scenario where recreating the Measurement type was unnecessary, to now include the scenario where DCStateEstimation also does not need to be recreated. Such efficiency can be particularly advantageous in cases where JuliaGrid can reuse gain matrix factorization.

We can modify the prior example to achieve the same model without establishing DCStateEstimation twice:

system = powerSystem()
device = measurement() # <- Initialize the Measurement instance

addBus!(system; label = "Bus 1", type = 3)
addBus!(system; label = "Bus 2", type = 1, active = 0.1)

addBranch!(system; label = "Branch 1", from = "Bus 1", to = "Bus 2", reactance = 0.5)

addGenerator!(system; label = "Generator 1", bus = "Bus 1", active = 0.1)

@wattmeter(label = "Wattmeter ?")
addWattmeter!(system, device; bus = "Bus 2", active = -0.11, variance = 1e-3)
addWattmeter!(system, device; from = "Branch 1", active = 0.09, variance = 1e-4)
addWattmeter!(system, device; to = "Branch 1", active = -0.12, variance = 1e-4, status = 0)

analysis = dcStateEstimation(system, device) # <- Build DCStateEstimation for the model
solve!(system, analysis)

updateWattmeter!(system, device; label = "Wattmeter 1", status = 0)
updateWattmeter!(system, device; label = "Wattmeter 2", active = 0.1, noise = false)
updateWattmeter!(system, device; label = "Wattmeter 3", status = 1)

# <- No need for re-build; we have already updated the existing DCStateEstimation instance
solve!(system, analysis)

Reusing Weighted Least-Squares Matrix Factorization

Drawing from the preceding example, our focus now shifts to finding a solution involving modifications that entail adjusting the measurement value of the Wattmeter 2. It is important to note that these adjustments do not impact the variance or status of the measurement device, which can affect the gain matrix. To resolve this updated system, users can simply execute the solve! function:

updateWattmeter!(system, device, analysis; label = "Wattmeter 2", active = 0.091)

solve!(system, analysis)

Power Analysis

After obtaining the solution from the DC state estimation, calculating powers related to buses and branches is facilitated by using the power! function. For instance, let us consider the model for which we obtained the DC state estimation solution:

system = powerSystem()
device = measurement()

addBus!(system; label = "Bus 1", type = 3, conductance = 1e-3)
addBus!(system; label = "Bus 2", type = 1, active = 0.2)
addBus!(system; label = "Bus 3", type = 1, active = 0.4)

addBranch!(system; label = "Branch 1", from = "Bus 1", to = "Bus 2", reactance = 0.5)
addBranch!(system; label = "Branch 2", from = "Bus 1", to = "Bus 3", reactance = 0.2)
addBranch!(system; label = "Branch 3", from = "Bus 2", to = "Bus 3", reactance = 0.3)

addGenerator!(system; label = "Generator 1", bus = "Bus 1", active = 3.2)

addWattmeter!(system, device; bus = "Bus 1", active = 0.6, variance = 1e-3)
addWattmeter!(system, device; bus = "Bus 3", active = -0.4, variance = 1e-2)
addWattmeter!(system, device; from = "Branch 1", active = 0.18, variance = 1e-4)
addWattmeter!(system, device; to = "Branch 2", active = -0.42, variance = 1e-4)

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

We can compute active powers using the following function:

power!(system, analysis)

For example, active power injections corresponding to buses are:

julia> print(system.bus.label, analysis.power.injection.active)Bus 1: 0.6008333333333332
Bus 2: -0.19983333333333342
Bus 3: -0.4

Print Results in the REPL

Users can utilize any of the print functions outlined in the Print API related to the DC analysis. For example, to print state estimation data related to wattmeters, we can use:

@power(MW, pu, pu)
printWattmeterData(system, device, analysis)

|---------------------------------------------------------------|
| Wattmeter Data                                                |
|---------------------------------------------------------------|
| Label |                     Active Power                      |
|       |                                                       |
|       | Measurement | Variance | Estimate | Residual | Status |
|       |        [MW] |     [MW] |     [MW] |     [MW] |        |
|-------|-------------|----------|----------|----------|--------|
| 1     |     60.0000 | 1.00e-01 |  60.0833 |  -0.0833 |      1 |
| 2     |    -40.0000 | 1.00e+00 | -40.0000 |   0.0000 |      1 |
| 3     |     18.0000 | 1.00e-02 |  17.9917 |   0.0083 |      1 |
| 4     |    -42.0000 | 1.00e-02 | -41.9917 |  -0.0083 |      1 |
|---------------------------------------------------------------|

Save Results to a CSV File

For CSV output, users should first generate a simple table with style = false, and then save it to a CSV file:

using CSV

io = IOBuffer()
printWattmeterData(system, device, analysis, io; style = false)
CSV.write("bus.csv", CSV.File(take!(io); delim = "|"))

Active Power Injection

To calculate active power injection associated with a specific bus, the function can be used:

julia> active = injectionPower(system, analysis; label = "Bus 1")0.6008333333333332

Active Power Injection from Generators

To calculate active power injection from the generators at a specific bus, the function can be used:

julia> active = supplyPower(system, analysis; label = "Bus 1")0.6008333333333332

Active Power Flow

Similarly, we can compute the active power flow at both the from-bus and to-bus ends of the specific branch by utilizing the provided functions below:

julia> active = fromPower(system, analysis; label = "Branch 1")0.17991666666666667
julia> active = toPower(system, analysis; label = "Branch 1")-0.17991666666666667