DC State Estimation

Performing DC state estimation first requires the PowerSystem type that has been created with the DC model, alongside the Measurement type that retains measurement data. Next, 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.

To obtain bus voltage angles and solve the DC state estimation problem, users can use:

• solve!.

After solving the DC state estimation, JuliaGrid provides a function for computing powers:

• power!.

Alternatively, users can call the wrapper function to solve the WLS or LAV model and, optionally, compute powers:

• stateEstimation!.

Users can also access specialized functions for computing specific types of powers for individual buses, branches, or generators within the power system.

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, monitoring = ems()

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

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)

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

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

analysis = dcStateEstimation(monitoring)

analysis = dcStateEstimation(monitoring, LDLt)

To obtain the bus voltage angles, use the solve! function as shown:

solve!(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

Wrapper Function

JuliaGrid provides the stateEstimation! wrapper function for DC state estimation. It solves the WLS or LAV model, and it can optionally compute powers:

analysis = dcStateEstimation(monitoring)
stateEstimation!(analysis; verbose = 2)

Number of entries in the coefficient matrix: 10
Number of measurement functions:              4
Number of state variables:                    2

EXIT: The solution of the DC state estimation was found.

Print Results in the REPL

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

@voltage(pu, deg)
printBusData(analysis)

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

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

printBusData(analysis; label = "Bus 1", header = true)
printBusData(analysis; label = "Bus 2")
printBusData(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:

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

Alternative Formulations

The conventional approach to solving the WLS state estimation problem generally performs well. However, it is well known that, under certain conditions often encountered in real-world systems, this method may suffer from numerical instabilities. These issues can prevent the algorithm from reaching a satisfactory solution. In such scenarios, users may choose to apply an alternative formulation of the WLS estimator.

Orthogonal Method

One such alternative is the orthogonal method [6, Sec. 3.2], which offers increased numerical robustness, particularly in cases where measurement variances differ significantly. By specifying the Orthogonal argument in the dcStateEstimation function, users can solve the WLS problem via QR factorization applied to a rectangular matrix formed by multiplying the square root of the precision matrix with the coefficient matrix:

analysis = dcStateEstimation(monitoring, Orthogonal)
stateEstimation!(analysis)

Peters and Wilkinson Method

Another option is the Peters and Wilkinson method [6, Sec. 3.4], which applies LU factorization to the same rectangular matrix, constructed using the square root of the precision matrix and the coefficient matrix. This method can be selected by passing the PetersWilkinson argument to the dcStateEstimation function:

analysis = dcStateEstimation(monitoring, PetersWilkinson)
stateEstimation!(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 analysis procedures [6, Ch. 6]. However, it is important to note that achieving robustness typically involves increased computational complexity.

To obtain an LAV estimator, users need to use 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(monitoring, Ipopt.Optimizer)

Setup Initial Primal Values

In JuliaGrid, the assignment of initial primal values for optimization variables takes place when the solve! function is executed. These values are derived from the voltage angles stored in the PowerSystem type and are assigned to the corresponding voltage field within the DcStateEstimation type:

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

Users can customize these values, which will be used as the initial primal values when executing the solve! function. One practical approach is to obtain the WLS estimator and then apply the resulting solution as the starting point for state estimation:

wls = dcStateEstimation(monitoring)
stateEstimation!(wls)

setInitialPoint!(analysis, wls)

As a result, the initial primal values will now reflect the outcome of the WLS state estimation solution:

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

Solution

To solve the formulated LAV state estimation model, call the wrapper function:

stateEstimation!(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.08999999999999997
Bus 3: -0.08399999999999999

Measurement Set Update

Begin by creating the PowerSystem and Measurement types with the ems function. The DC model is then configured using the dcModel! function. After that, initialize the DcStateEstimation type through the dcStateEstimation function and solve the resulting state estimation problem:

system, monitoring = ems()

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

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

dcModel!(system)

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

analysis = dcStateEstimation(monitoring)
stateEstimation!(analysis)

Next, modify the existing Measurement type using add and update functions. Then, create the new DcStateEstimation type based on the modified system and solve the state estimation problem:

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

analysis = dcStateEstimation(monitoring)
stateEstimation!(analysis)

State Estimation Update

For advanced workflows, users can create the DcStateEstimation type once using dcStateEstimation or dcLavStateEstimation. They can then modify existing measurement devices without recreating the DcStateEstimation type.

This approach extends the previous workflow by also avoiding recreation of the DcStateEstimation object. This is particularly useful when JuliaGrid can reuse the gain matrix factorization.

Now revisit the defined PowerSystem, Measurement and DcStateEstimation types:

system, monitoring = ems()

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

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

dcModel!(system)

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

analysis = dcStateEstimation(monitoring)
stateEstimation!(analysis)

Next, modify the existing Measurement type as well as the DcStateEstimation type using add and update functions. Then, immediately solve the state estimation problem:

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

stateEstimation!(analysis)

Reusing Weighted Least-Squares Matrix Factorization

Continuing from the preceding example, consider modifications that adjust the measurement value of the Wattmeter 2. 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 call the stateEstimation! wrapper function:

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

stateEstimation!(analysis)

Power Analysis

After obtaining the solution from the DC state estimation, use the power! function to calculate powers related to buses and branches. For instance, consider the model used to obtain the DC state estimation solution:

system, monitoring = ems()

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

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)

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

analysis = dcStateEstimation(monitoring)
stateEstimation!(analysis)

Compute active powers using the following function:

power!(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 use 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, use:

@power(MW, pu)
printWattmeterData(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(analysis, io; style = false)
CSV.write("wattmeter.csv", CSV.File(take!(io); delim = "|"))

Active Power Injection

To calculate the active power injection associated with a specific bus, use:

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

Active Power Injection from Generators

To calculate the active power supply associated with a specific bus, use:

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

Active Power Flow

Similarly, to compute the active power flow at both the from-bus and to-bus ends of a specific branch, use:

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