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. Subsequently, we can 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 the wrapper function:

solve!.

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

Alternatively, instead of using functions responsible for solving state estimation and computing powers and currents, users can use the wrapper function:

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

Let us 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 Placment

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. In this instance, 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, we can 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

Let us continue with the previous example, where we defined the PowerSystem and Measurement types. To establish the PMU state estimation model, we will 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. However, the user also has the option to select alternative factorization methods such as LDLt or QR:

analysis = pmuStateEstimation(monitoring, QR)

To obtain the bus voltage magnitudes and angles, the solve! function can be invoked 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.00000000000074, 1.324710038975461e-13
Bus 2: 1.0215628522083289, 0.014634612107890924
Bus 3: 1.012294100414869, -0.005105099385215231

Info

We recommend that readers refer to the tutorial on PMU State Estimation for insights into the implementation.

Wrapper Function

JuliaGrid provides a wrapper function for PMU state estimation analysis and also supports the computation of powers and currents using the stateEstimation! function:

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

In the above approach, we assume that measurement errors from a single PMU are uncorrelated. This assumption leads to the covariance matrix and its inverse matrix (i.e., precision matrix) maintaining a diagonal form:

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 maintaining a diagonal form but instead becoming a block diagonal matrix.

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

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

Following this, we recreate the WLS state estimation model:

analysis = pmuStateEstimation(monitoring)

Upon inspection, it becomes evident that the precision matrix no longer maintains a diagonal structure:

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

Subsequently, we can 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.9988543535974568, 5.9920242181893276e-5
Bus 2: 1.020417744428146, 0.014710234393803864
Bus 3: 1.011147706933703, -0.0050522361923871885

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

This approach is suitable when measurement errors are uncorrelated, and the precision matrix remains diagonal. Therefore, as a preliminary step, we need to eliminate the correlation, as we did previously:

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

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

analysis = pmuStateEstimation(monitoring, Orthogonal)
stateEstimation!(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)
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 easily 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 as follows:

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

Tip

We also provide 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.

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 [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 = 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 have the flexibility to customize these values according to their requirements, and they will be utilized 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, simply execute the following 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

We begin by creating the PowerSystem and Measurement types with the ems function. The AC model is then configured using acModel! function. After that, we 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, we modify the existing Measurement type using add and update functions. Then, we 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

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

This advancement extends beyond the previous scenario where recreating the Measurement type was unnecessary, to now include the scenario where PmuStateEstimation also does not need to be recreated.

Tip

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

Let us now revisit our 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, we modify the existing Measurement type as well as the PmuStateEstimation type using add and update functions. We then immediately proceed to 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 employing the WLS model, JuliaGrid can reuse the symbolic factorizations of LU or LDLt, provided that the nonzero pattern of the gain matrix remains unchanged.

Power and Current Analysis

After obtaining the solution from the PMU state estimation, we can calculate various electrical quantities related to buses and branches using the power! and current! functions. For instance, let us consider the model for which we obtained 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)

We can now utilize the provided functions to compute powers and currents:

power!(analysis)
current!(analysis)

For instance, if we want to show the active power injections and the from-bus current magnitudes, we can employ the following code:

julia> print(system.bus.label, analysis.power.injection.active)Bus 1: 2.7128789695735662
Bus 2: -0.2156411013050359
Bus 3: -2.4028824415317627
julia> print(system.branch.label, analysis.current.from.magnitude)Branch 1: 1.0848042244961928
Branch 2: 1.6624663614086033
Branch 3: 0.8658336251367414

Info

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

Print Results in the REPL

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

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

|---------------------------------------------------------------|
| PMU Data                                                      |
|---------------------------------------------------------------|
| Label |                   Voltage Magnitude                   |
|       |                                                       |
|       | Measurement | Variance | Estimate | Residual | Status |
|       |        [pu] |     [pu] |     [pu] |     [pu] |        |
|-------|-------------|----------|----------|----------|--------|
| 1     |      1.0000 | 1.00e-08 |   1.0001 |  -0.0001 |      1 |
| 2     |      0.9700 | 1.00e-08 |   0.9700 |   0.0000 |      1 |
|---------------------------------------------------------------|
|---------------------------------------------------------------|
| PMU Data                                                      |
|---------------------------------------------------------------|
| Label |                   Current Magnitude                   |
|       |                                                       |
|       | Measurement | Variance | Estimate | Residual | Status |
|       |        [pu] |     [pu] |     [pu] |     [pu] |        |
|-------|-------------|----------|----------|----------|--------|
| 3     |      1.6600 | 1.00e-08 |   1.6625 |  -0.0025 |      1 |
| 4     |      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("bus.csv", CSV.File(take!(io); delim = "|"))

Active and Reactive Power Injection

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

julia> active, reactive = injectionPower(analysis; label = "Bus 1")(2.7128789695735662, 0.4328846046307191)

Active and Reactive Power Injection from Generators

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

julia> active, reactive = supplyPower(analysis; label = "Bus 1")(2.7128789695735662, 0.4328846046307191)

Active and Reactive Power at Shunt Element

To calculate the active and reactive power associated with shunt element at a specific bus, the function can be used:

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

Active and Reactive Power Flow

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

julia> active, reactive = fromPower(analysis; label = "Branch 2")(1.6429108157160377, 0.25539133163697453)
julia> active, reactive = toPower(analysis; label = "Branch 2")(-1.5873301378900675, -0.1549833506017298)

Active and Reactive Power at Charging Admittances

To calculate the active and reactive power linked with branch charging admittances of the particular branch, the function can be used:

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

Active powers indicate active losses within the branch's charging admittances. Moreover, charging admittances injected reactive powers 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, the function can be used:

julia> active, reactive = seriesPower(analysis; label = "Branch 2")(0.05548491704701507, 0.13871229261753767)

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, the function can be used:

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

Current Flow

We can compute the current flow at both the from-bus and to-bus ends of the specific branch by utilizing the provided functions below:

julia> magnitude, angle = fromCurrent(analysis; label = "Branch 2")(1.6624663614086033, -0.15416239142972446)
julia> magnitude, angle = toCurrent(analysis; label = "Branch 2")(1.6673093443308689, 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, we can use the following function:

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