PMU State Estimation

To perform linear state estimation solely based on PMU data, the initial requirement is to have the PowerSystem type configured with the AC model, along with the Measurement type storing measurement data. Next, formulate either the weighted least-squares (WLS) or the least absolute value (LAV) PMU state estimation model encapsulated within the type PmuStateEstimation using:

To obtain bus voltages and solve the PMU state estimation problem, users can use:

solve!.

After solving the PMU state estimation, JuliaGrid provides functions for computing powers and currents:

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

stateEstimation!.

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

Phasor Measurements

Define the PowerSystem type and perform the AC power flow analysis solely for generating data to artificially create measurement values:

system = powerSystem()

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

@branch(resistance = 0.02, conductance = 1e-4, susceptance = 0.04)
addBranch!(system; label = "Branch 1", from = "Bus 1", to = "Bus 2", reactance = 0.05)
addBranch!(system; label = "Branch 2", from = "Bus 1", to = "Bus 2", reactance = 0.01)
addBranch!(system; label = "Branch 3", from = "Bus 2", to = "Bus 3", reactance = 0.04)

@generator(reactive = 0.1)
addGenerator!(system; label = "Generator 1", bus = "Bus 1", active = 3.2)
addGenerator!(system; label = "Generator 2", bus = "Bus 2", active = 2.1)

pf = newtonRaphson(system)
powerFlow!(pf)

Optimal PMU Placement

After defining the PowerSystem type, the next step is to define the Measurement type:

monitoring = measurement(system)

JuliaGrid allows users to determine the minimum number of PMUs needed for observability while also generating phasor measurements based on results from AC power flow through the pmuPlacement! function:

using HiGHS

@pmu(label = "PMU ? (!)")
placement = pmuPlacement!(monitoring, pf, HiGHS.Optimizer)

Note that users can also generate phasor measurements using results from AC optimal power flow.

The placement variable contains data regarding the optimal placement of measurements. Here, installing a PMU at Bus 2 renders the system observable:

julia> keys(placement.bus)KeySet for a OrderedCollections.OrderedDict{String, Int64} with 1 entry. Keys:
  "Bus 2"

This PMU installed at Bus 2 will measure the bus voltage phasor at the corresponding bus and all current phasors at the branches incident to Bus 2 located at the from-bus or to-bus ends:

julia> keys(placement.from)KeySet for a OrderedCollections.OrderedDict{String, Int64} with 1 entry. Keys:
  "Branch 3"
julia> keys(placement.to)KeySet for a OrderedCollections.OrderedDict{String, Int64} with 2 entries. Keys:
  "Branch 1"
  "Branch 2"

Finally, observe the obtained set of measurement values:

julia> print(monitoring.pmu.label, monitoring.pmu.magnitude.mean, monitoring.pmu.angle.mean)PMU 1 (Bus 2): 1.0215628522075866, 0.014634612107772262
PMU 2 (From Branch 3): 0.49529892193563957, 0.07710792493557296
PMU 3 (To Branch 1): 0.4744352810029824, -0.545987706954647
PMU 4 (To Branch 2): 1.172275987225155, 0.16223095481110425

Weighted Least-Squares Estimator

Continue with the previous example, where the PowerSystem and Measurement types were defined. To establish the PMU state estimation model, use the pmuStateEstimation function:

analysis = pmuStateEstimation(monitoring)

Tip

Here, the user triggers LU factorization as the default method for solving the PMU state estimation problem. The available factorization methods are LL, LDLt, LU, KLU, and QR:

analysis = pmuStateEstimation(monitoring, QR)

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

solve!(analysis)

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

julia> print(system.bus.label, analysis.voltage.magnitude, analysis.voltage.angle)Bus 1: 1.000000000000723, 1.4496499478317466e-13
Bus 2: 1.021562852208312, 0.014634612107903407
Bus 3: 1.0122941004148518, -0.005105099385202981

Info

For implementation details, see the tutorial on PMU State Estimation.

Wrapper Function

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

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

Number of entries in the coefficient matrix: 26
Number of measurement functions:              8
Number of state variables:                    6

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

Correlated Measurement Errors

The approach above assumes that measurement errors from a single PMU are uncorrelated. This assumption leads to the covariance matrix and its inverse matrix (i.e., precision matrix) remaining diagonal:

julia> analysis.method.precision8×8 SparseArrays.SparseMatrixCSC{Float64, Int64} with 8 stored entries:
 9.99991e7   ⋅          ⋅        …   ⋅          ⋅          ⋅
  ⋅         9.58239e7   ⋅            ⋅          ⋅          ⋅
  ⋅          ⋅         1.0045e8      ⋅          ⋅          ⋅
  ⋅          ⋅          ⋅            ⋅          ⋅          ⋅
  ⋅          ⋅          ⋅            ⋅          ⋅          ⋅
  ⋅          ⋅          ⋅        …  2.30398e8   ⋅          ⋅
  ⋅          ⋅          ⋅            ⋅         9.90331e7   ⋅
  ⋅          ⋅          ⋅            ⋅          ⋅         7.32887e7

While this approach is suitable for many scenarios, linear PMU state estimation relies on transforming from polar to rectangular coordinate systems. Consequently, measurement errors from a single PMU become correlated due to this transformation. This correlation results in the covariance matrix, and hence the precision matrix, no longer remaining diagonal but instead becoming a block diagonal matrix.

To accommodate this, users can consider correlation when adding each PMU to the Measurement type. For instance, add a new PMU while considering correlation:

addPmu!(monitoring; bus = "Bus 3", magnitude = 1.01, angle = -0.005, correlated = true)

Following this, recreate the WLS state estimation model:

analysis = pmuStateEstimation(monitoring)

Upon inspection, it becomes evident that the precision matrix is no longer diagonal:

julia> analysis.method.precision10×10 SparseArrays.SparseMatrixCSC{Float64, Int64} with 12 stored entries:
 9.99991e7   ⋅          ⋅        …   ⋅              ⋅          ⋅
  ⋅         9.58239e7   ⋅            ⋅              ⋅          ⋅
  ⋅          ⋅         1.0045e8      ⋅              ⋅          ⋅
  ⋅          ⋅          ⋅            ⋅              ⋅          ⋅
  ⋅          ⋅          ⋅            ⋅              ⋅          ⋅
  ⋅          ⋅          ⋅        …   ⋅              ⋅          ⋅
  ⋅          ⋅          ⋅            ⋅              ⋅          ⋅
  ⋅          ⋅          ⋅           7.32887e7       ⋅          ⋅
  ⋅          ⋅          ⋅            ⋅             1.0e8  -9851.81
  ⋅          ⋅          ⋅            ⋅         -9851.81       9.80297e7

Next, address this new scenario and observe the solution:

julia> stateEstimation!(analysis)
julia> print(system.bus.label, analysis.voltage.magnitude, analysis.voltage.angle)Bus 1: 0.9988543535974481, 5.9920242188089627e-5
Bus 2: 1.0204177444281375, 0.014710234393810051
Bus 3: 1.0111477069336943, -0.005052236192381124

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       |
|       |                     |
|   Bus | Magnitude |   Angle |
|       |      [pu] |   [deg] |
|-------|-----------|---------|
| Bus 1 |    0.9989 |  0.0034 |
| Bus 2 |    1.0204 |  0.8428 |
| Bus 3 |    1.0111 | -0.2895 |
|-----------------------------|

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

Tip

JuliaGrid also provides functions to print or save state estimation results, such as estimated values and residuals. For more details, users can consult the Power and Current Analysis section of this manual.

Alternative Formulations

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 scenarios, users may choose to apply an alternative formulation of the WLS estimator.

These alternative methods are applicable when measurement errors are uncorrelated and the precision matrix is diagonal. Therefore, as a preliminary step, the correlation must be eliminated, as shown previously:

updatePmu!(monitoring; label = "PMU 5 (Bus 3)", correlated = false)

Orthogonal Method

One such alternative is the orthogonal method [6, Sec. 3.2], which offers increased numerical robustness, especially when measurement variances differ significantly. This method solves the WLS problem using QR factorization applied to a rectangular matrix formed by multiplying the square root of the precision matrix with the coefficient matrix. To enable this method, specify the Orthogonal argument in the pmuStateEstimation function:

analysis = pmuStateEstimation(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 pmuStateEstimation function:

analysis = pmuStateEstimation(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 = pmuLavStateEstimation(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 magnitudes and angles stored in the PowerSystem type and are assigned to the corresponding voltage field within the PmuStateEstimation type:

julia> print(system.bus.label, analysis.voltage.magnitude, analysis.voltage.angle)Bus 1: 1.0, 0.0
Bus 2: 1.0, 0.0
Bus 3: 1.0, 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 perform an AC power flow analysis and then apply the resulting solution as the starting point for state estimation:

pf = newtonRaphson(system)
powerFlow!(pf)

setInitialPoint!(analysis, pf)

As a result, the initial primal values will now reflect the outcome of the power flow solution:

julia> print(system.bus.label, analysis.voltage.magnitude, analysis.voltage.angle)Bus 1: 1.0, 0.0
Bus 2: 1.0215628522075866, 0.014634612107772262
Bus 3: 1.0122941004141301, -0.005105099385349472

Solution

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

stateEstimation!(analysis)

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

julia> print(system.bus.label, analysis.voltage.magnitude, analysis.voltage.angle)Bus 1: 0.9999996847875078, 1.358655429134381e-7
Bus 2: 1.0215625389870358, 0.01463474976908131
Bus 3: 1.012293784667062, -0.005104966819778244

Info

Readers can refer to the Least Absolute Value Estimation tutorial for implementation insights.

Measurement Set Update

Begin by creating the PowerSystem and Measurement types with the ems function. The AC model is then configured using the acModel! function. After that, initialize the PmuStateEstimation type through the pmuStateEstimation function and solve the resulting state estimation problem:

system, monitoring = ems()

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

@branch(resistance = 0.02, conductance = 1e-4, susceptance = 0.04)
addBranch!(system; label = "Branch 1", from = "Bus 1", to = "Bus 2", reactance = 0.05)
addBranch!(system; label = "Branch 2", from = "Bus 2", to = "Bus 3", reactance = 0.03)

acModel!(system)

addPmu!(monitoring; label = "PMU 1", bus = "Bus 1", magnitude = 1.0, angle = 0.0)
addPmu!(monitoring; label = "PMU 2", bus = "Bus 2", magnitude = 0.98, angle = -0.023)
addPmu!(monitoring; label = "PMU 3", from = "Branch 2", magnitude = 0.5, angle = -0.05)

analysis = pmuStateEstimation(monitoring)
stateEstimation!(analysis)

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

addPmu!(monitoring; label = "PMU 4", to = "Branch 2", magnitude = 0.5, angle = 3)
updatePmu!(monitoring; label = "PMU 1", varianceMagnitude = 1e-8)
updatePmu!(monitoring; label = "PMU 3", status = 0)

analysis = pmuStateEstimation(monitoring)
stateEstimation!(analysis)

Info

This concept removes the need to restart and recreate the Measurement type from the beginning when implementing changes to the existing measurement set.

State Estimation Update

For advanced workflows, users can create the PmuStateEstimation type once using pmuStateEstimation or pmuLavStateEstimation. They can then modify existing measurement devices without recreating the PmuStateEstimation type.

This approach extends the previous workflow by also avoiding recreation of the PmuStateEstimation object.

Tip

The addition of new measurements after the creation of PmuStateEstimation is not practical in terms of reusing the PmuStateEstimation type. Instead, users should create a final set of measurements and then use update functions to manage devices, either setting them in-service or out-of-service throughout the process.

Now revisit the defined PowerSystem, Measurement and PmuStateEstimation types:

system, monitoring = ems()

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

@branch(resistance = 0.02, conductance = 1e-4, susceptance = 0.04)
addBranch!(system; label = "Branch 1", from = "Bus 1", to = "Bus 2", reactance = 0.05)
addBranch!(system; label = "Branch 2", from = "Bus 2", to = "Bus 3", reactance = 0.03)

acModel!(system)

addPmu!(monitoring; label = "PMU 1", bus = "Bus 1", magnitude = 1.0, angle = 0.0)
addPmu!(monitoring; label = "PMU 2", bus = "Bus 2", magnitude = 0.98, angle = -0.023)
addPmu!(monitoring; label = "PMU 3", from = "Branch 2", magnitude = 0.5, angle = -0.05)
addPmu!(monitoring; label = "PMU 4", to = "Branch 2", magnitude = 0.5, angle = 3, status = 0)

analysis = pmuStateEstimation(monitoring)
stateEstimation!(analysis)

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

updatePmu!(analysis; label = "PMU 1", varianceMagnitude = 1e-8)
updatePmu!(analysis; label = "PMU 3", status = 0)
updatePmu!(analysis; label = "PMU 4", status = 1)

stateEstimation!(analysis)

Info

This concept removes the need to rebuild both the Measurement and the PmuStateEstimation from the beginning when implementing changes to the existing measurement set. In the scenario of using the WLS model, JuliaGrid can reuse symbolic factorizations of LL, LDLt, LU, and KLU, provided that the nonzero pattern of the gain matrix remains unchanged.

Power and Current Analysis

After obtaining the solution from the PMU state estimation, calculate various electrical quantities related to buses and branches using the power! and current! functions. For instance, consider the model used to obtain the PMU state estimation solution:

system, monitoring = ems()

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

@branch(resistance = 0.02, conductance = 1e-4, susceptance = 0.04)
addBranch!(system; label = "Branch 1", from = "Bus 1", to = "Bus 2", reactance = 0.05)
addBranch!(system; label = "Branch 2", from = "Bus 1", to = "Bus 3", reactance = 0.05)
addBranch!(system; label = "Branch 3", from = "Bus 2", to = "Bus 3", reactance = 0.03)

addPmu!(monitoring; bus = "Bus 1", magnitude = 1.0, angle = 0.0)
addPmu!(monitoring; bus = "Bus 2", magnitude = 0.97, angle = -0.051)
addPmu!(monitoring; from = "Branch 2", magnitude = 1.66, angle = -0.15)
addPmu!(monitoring; to = "Branch 2", magnitude = 1.67, angle = 2.96)

analysis = pmuStateEstimation(monitoring)
stateEstimation!(analysis)

Now use the provided functions to compute powers and currents:

power!(analysis)
current!(analysis)

For instance, to show the active power injections and the from-bus current magnitudes, use the following code:

julia> print(system.bus.label, analysis.power.injection.active)Bus 1: 2.712878969573828
Bus 2: -0.21564110130564348
Bus 3: -2.40288244153142
julia> print(system.branch.label, analysis.current.from.magnitude)Branch 1: 1.0848042244964458
Branch 2: 1.6624663614086042
Branch 3: 0.8658336251363701

Info

To better understand the powers and currents associated with buses and branches that are calculated by the power! and current! functions, see the tutorials on PMU State Estimation.

Print Results in the REPL

Users can use any of the print functions outlined in the Print API. For example, to print state estimation data related to PMUs, use:

@voltage(pu, deg)
show = Dict("Voltage Angle" => false, "Current Angle" => false)
printPmuData(analysis; show)

|-------------------------------------------------------|
| PMU Data                                              |
|-------------------------------------------------------|
|                   Voltage Magnitude                   |
|                                                       |
| Measurement | Variance | Estimate | Residual | Status |
|        [pu] |     [pu] |     [pu] |     [pu] |        |
|-------------|----------|----------|----------|--------|
|      1.0000 | 1.00e-08 |   1.0001 |  -0.0001 |      1 |
|      0.9700 | 1.00e-08 |   0.9700 |   0.0000 |      1 |
|-------------------------------------------------------|
|-------------------------------------------------------|
| PMU Data                                              |
|-------------------------------------------------------|
|                   Current Magnitude                   |
|                                                       |
| Measurement | Variance | Estimate | Residual | Status |
|        [pu] |     [pu] |     [pu] |     [pu] |        |
|-------------|----------|----------|----------|--------|
|      1.6600 | 1.00e-08 |   1.6625 |  -0.0025 |      1 |
|      1.6700 | 1.00e-08 |   1.6673 |   0.0027 |      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()
printPmuData(analysis, io; style = false)
CSV.write("pmu.csv", CSV.File(take!(io); delim = "|"))

Active and Reactive Power Injection

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

julia> active, reactive = injectionPower(analysis; label = "Bus 1")(2.712878969573828, 0.432884604630847)

Active and Reactive Power Injection from Generators

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

julia> active, reactive = supplyPower(analysis; label = "Bus 1")(2.712878969573828, 0.432884604630847)

Active and Reactive Power at Shunt Element

To calculate the active and reactive power associated with the shunt element at a specific bus, use:

julia> active, reactive = shuntPower(analysis; label = "Bus 1")(0.0, -0.0020004242557668814)

Active and Reactive Power Flow

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

julia> active, reactive = fromPower(analysis; label = "Branch 2")(1.6429108157160512, 0.2553913316369938)
julia> active, reactive = toPower(analysis; label = "Branch 2")(-1.5873301378900804, -0.1549833506017493)

Active and Reactive Power at Charging Admittances

To calculate the active and reactive power associated with the branch charging admittances of a particular branch, use:

julia> active, reactive = chargingPower(analysis; label = "Branch 1")(9.705560639417204e-5, -0.038822242557668814)

Active powers indicate active losses within the branch's charging admittances. Moreover, charging admittances inject reactive power into the power system due to their capacitive nature, as denoted by a negative sign.

Active and Reactive Power at Series Impedance

To calculate the active and reactive power across the series impedance of the branch, use:

julia> active, reactive = seriesPower(analysis; label = "Branch 2")(0.055484917047015136, 0.13871229261753784)

The active power also considers active losses originating from the series resistance of the branch, while the reactive power represents reactive losses resulting from the impedance's inductive characteristics.

Current Injection

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

julia> magnitude, angle = injectionCurrent(analysis; label = "Bus 1")(2.7469074926858754, -0.1581787316054701)

Current Flow

To compute the current flow at both the from-bus and to-bus ends of a specific branch, use:

julia> magnitude, angle = fromCurrent(analysis; label = "Branch 2")(1.6624663614086042, -0.15416239142972438)
julia> magnitude, angle = toCurrent(analysis; label = "Branch 2")(1.6673093443308697, 2.964124299788867)

Current Through Series Impedance

To calculate the current passing through the series impedance of the branch in the direction from the from-bus end to the to-bus end, use the following function:

julia> magnitude, angle = seriesCurrent(analysis; label = "Branch 2")(1.6656067520128384, -0.16603367260017451)