PMU State Estimation

This example examines a 6-bus power system, illustrated in Figure 1. The goal is to estimate bus voltage magnitudes and angles using only phasor measurements.

Figure 1: The 6-bus power system.

Info

Users can download a Julia script containing the scenarios from this section using the following link.

The power system is defined by specifying buses and branches, with the generator assigned to the slack bus:

system = powerSystem()

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

@branch(resistance = 0.02)
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.23)
addBranch!(system; label = "Branch 3", from = "Bus 3", to = "Bus 4", reactance = 0.19)
addBranch!(system; label = "Branch 4", from = "Bus 4", to = "Bus 5", reactance = 0.17)
addBranch!(system; label = "Branch 5", from = "Bus 5", to = "Bus 6", reactance = 0.04)
addBranch!(system; label = "Branch 6", from = "Bus 1", to = "Bus 6", reactance = 0.21)
addBranch!(system; label = "Branch 7", from = "Bus 2", to = "Bus 6", reactance = 0.13)
addBranch!(system; label = "Branch 8", from = "Bus 5", to = "Bus 2", reactance = 0.34)

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

Display Data Settings

Before running simulations, we configure the data display settings:

show = Dict("Shunt Power" => false, "Status" => false, "Series Power" => false)

Optimal PMU Placement

Next, PMUs need to be assigned to the power system shown in Figure 1. The placement is determined using an optimal PMU placement strategy that ensures observability with the minimal number of PMUs:

placement = pmuPlacement(system, HiGHS.Optimizer; print = false)

This provides the bus configuration where PMUs should be installed:

julia> placement.busOrderedCollections.OrderedDict{String, Int64} with 2 entries:
  "Bus 2" => 2
  "Bus 5" => 5

PMUs installed at these buses will measure bus voltage phasors and all currents in branches connected to those buses. Specifically, current phasors will be measured in the following branches:

julia> placement.fromOrderedCollections.OrderedDict{String, Int64} with 4 entries:
  "Branch 2" => 2
  "Branch 7" => 7
  "Branch 5" => 5
  "Branch 8" => 8
julia> placement.toOrderedCollections.OrderedDict{String, Int64} with 3 entries:
  "Branch 1" => 1
  "Branch 8" => 8
  "Branch 4" => 4

If users choose to generate phasor measurement values using optimal power flow or power flow analysis, the integers within bus and branch labels indicate positions in vectors where these values are stored.

Hence, Figure 2 illustrates the phasor measurement configuration, which includes bus voltage and branch current phasor measurements.

Figure 2: The 6-bus power system with phasor measurement configuration.

Finally, phasor measurements need to be defined. The question is how to obtain measurement values. In this example, AC power flow results will be used to generate synthetic measurements.

Measurement Model

To begin, let us initialize the measurement variable:

device = measurement()

AC Power Flow

AC power flow analysis requires generator and load data. The system configuration is shown in Figure 3.

Figure 3: The 6-bus power system with generators and loads.

Data for loads and generators is added as follows:

updateBus!(system; label = "Bus 2", type = 1, active = 0.217, reactive = 0.127)
updateBus!(system; label = "Bus 3", type = 1, active = 0.478, reactive = -0.039)
updateBus!(system; label = "Bus 4", type = 2, active = 0.076, reactive = 0.016)
updateBus!(system; label = "Bus 5", type = 1, active = 0.112, reactive = 0.075)
updateBus!(system; label = "Bus 6", type = 1, active = 0.295, reactive = 0.166)

updateGenerator!(system; label = "Generator 1", active = 2.324, reactive = -0.169)
addGenerator!(system; label = "Generator 2", bus = "Bus 4", active = 0.412, reactive = 0.234)

Next, AC power flow analysis is performed to obtain bus voltages:

acModel!(system)

powerFlow = newtonRaphson(system)
for iteration = 1:20
    stopping = mismatch!(system, powerFlow)
    if all(stopping .< 1e-8)
        println("The algorithm converged in $(iteration - 1) iterations.")
        break
    end
    solve!(system, powerFlow)
end

The algorithm converged in 3 iterations.

Bus Voltage Phasor Measurements

To obtain bus voltage phasor measurements, the exact values and optimal PMU placement data are used. By setting noise = true, white Gaussian noise with a default variance of 1e-8 is added to the magnitude and angle values:

@pmu(label = "!")
for (bus, idx) in placement.bus
    Vᵢ, θᵢ = powerFlow.voltage.magnitude[idx], powerFlow.voltage.angle[idx]
    addPmu!(system, device; bus = bus, magnitude = Vᵢ, angle = θᵢ, noise = true)
end

Branch Current Phasor Measurements

To add branch current phasor measurements, the current magnitudes and angles are first computed. These values are then used to form measurements, where the exact values are used as the noise keyword is ignored:

for branch in keys(placement.from)
    Iᵢⱼ, ψᵢⱼ = fromCurrent(system, powerFlow; label = branch)
    addPmu!(system, device; from = branch, magnitude = Iᵢⱼ, angle = ψᵢⱼ)
end
for branch in keys(placement.to)
    Iⱼᵢ, ψⱼᵢ = toCurrent(system, powerFlow; label = branch)
    addPmu!(system, device; to = branch, magnitude = Iⱼᵢ, angle = ψⱼᵢ)
end

Current phasor measurements can also be generated in the same way as voltage phasors by invoking the current! function after AC state estimation has converged.

Note that the formation of bus voltage and branch current measurements can be performed by calling the pmuPlacement! function, as demonstrated in the PMU State Estimation manual.

Phasor Measurements

Finally, the complete set of phasor measurements is observed, as illustrated in Figure 2:

printPmuData(system, device; width = Dict("Label" => 15))

|-----------------------------------------------------------------------------------|
| PMU Data                                                                          |
|-----------------------------------------------------------------------------------|
|     Label     |        Voltage Magnitude        |          Voltage Angle          |
|               |                                 |                                 |
|               | Measurement | Variance | Status | Measurement | Variance | Status |
|               |        [pu] |     [pu] |        |       [rad] |    [rad] |        |
|---------------|-------------|----------|--------|-------------|----------|--------|
| Bus 2         |      0.9818 | 1.00e-08 |      1 |     -0.0267 | 1.00e-08 |      1 |
| Bus 5         |      0.9782 | 1.00e-08 |      1 |     -0.0408 | 1.00e-08 |      1 |
|-----------------------------------------------------------------------------------|
|-----------------------------------------------------------------------------------|
| PMU Data                                                                          |
|-----------------------------------------------------------------------------------|
|     Label     |        Current Magnitude        |          Current Angle          |
|               |                                 |                                 |
|               | Measurement | Variance | Status | Measurement | Variance | Status |
|               |        [pu] |     [pu] |        |       [rad] |    [rad] |        |
|---------------|-------------|----------|--------|-------------|----------|--------|
| From Branch 2 |      0.2282 | 1.00e-08 |      1 |      0.1829 | 1.00e-08 |      1 |
| From Branch 7 |      0.1058 | 1.00e-08 |      1 |     -0.3102 | 1.00e-08 |      1 |
| From Branch 5 |      0.0508 | 1.00e-08 |      1 |     -1.5712 | 1.00e-08 |      1 |
| From Branch 8 |      0.0414 | 1.00e-08 |      1 |      2.8992 | 1.00e-08 |      1 |
| To Branch 1   |      0.5984 | 1.00e-08 |      1 |      2.9124 | 1.00e-08 |      1 |
| To Branch 8   |      0.0414 | 1.00e-08 |      1 |     -0.2424 | 1.00e-08 |      1 |
| To Branch 4   |      0.1413 | 1.00e-08 |      1 |      2.0979 | 1.00e-08 |      1 |
|-----------------------------------------------------------------------------------|

Base Case Analysis

Once the measurements are obtained, the state estimation model is created:

analysis = pmuStateEstimation(system, device)

Next, the model is solved to obtain the WLS estimator for bus voltages, and the results are used to compute powers:

solve!(system, analysis)
power!(system, analysis)

This enables users to observe the estimated bus voltages along with the corresponding power values:

printBusData(system, analysis; show)

|------------------------------------------------------------------------------------------|
| Bus Data                                                                                 |
|------------------------------------------------------------------------------------------|
| Label |       Voltage       | Power Generation  |   Power Demand    |  Power Injection   |
|       |                     |                   |                   |                    |
|   Bus | Magnitude |   Angle | Active | Reactive | Active | Reactive |  Active | Reactive |
|       |      [pu] |   [rad] |   [pu] |     [pu] |   [pu] |     [pu] |    [pu] |     [pu] |
|-------|-----------|---------|--------|----------|--------|----------|---------|----------|
| Bus 1 |    1.0000 |  0.0001 | 0.7772 |   0.2349 | 0.0000 |   0.0000 |  0.7772 |   0.2349 |
| Bus 2 |    0.9819 | -0.0269 | 0.0000 |  -0.0000 | 0.2170 |   0.1270 | -0.2170 |  -0.1270 |
| Bus 3 |    0.9897 | -0.0797 | 0.0000 |  -0.0000 | 0.4780 |  -0.0390 | -0.4780 |   0.0390 |
| Bus 4 |    1.0000 | -0.0302 | 0.4120 |   0.1687 | 0.0760 |   0.0160 |  0.3360 |   0.1527 |
| Bus 5 |    0.9781 | -0.0407 | 0.0000 |  -0.0000 | 0.1120 |   0.0750 | -0.1120 |  -0.0750 |
| Bus 6 |    0.9761 | -0.0398 | 0.0000 |  -0.0000 | 0.2950 |   0.1660 | -0.2950 |  -0.1660 |
|------------------------------------------------------------------------------------------|

Users can also compare these results with those obtained from AC power flow:

power!(system, powerFlow)
printBusData(system, powerFlow; show)

|------------------------------------------------------------------------------------------|
| Bus Data                                                                                 |
|------------------------------------------------------------------------------------------|
| Label |       Voltage       | Power Generation  |   Power Demand    |  Power Injection   |
|       |                     |                   |                   |                    |
|   Bus | Magnitude |   Angle | Active | Reactive | Active | Reactive |  Active | Reactive |
|       |      [pu] |   [rad] |   [pu] |     [pu] |   [pu] |     [pu] |    [pu] |     [pu] |
|-------|-----------|---------|--------|----------|--------|----------|---------|----------|
| Bus 1 |    1.0000 |  0.0000 | 0.7773 |   0.2349 | 0.0000 |   0.0000 |  0.7773 |   0.2349 |
| Bus 2 |    0.9819 | -0.0269 | 0.0000 |   0.0000 | 0.2170 |   0.1270 | -0.2170 |  -0.1270 |
| Bus 3 |    0.9898 | -0.0798 | 0.0000 |   0.0000 | 0.4780 |  -0.0390 | -0.4780 |   0.0390 |
| Bus 4 |    1.0000 | -0.0302 | 0.4120 |   0.1687 | 0.0760 |   0.0160 |  0.3360 |   0.1527 |
| Bus 5 |    0.9782 | -0.0408 | 0.0000 |   0.0000 | 0.1120 |   0.0750 | -0.1120 |  -0.0750 |
| Bus 6 |    0.9761 | -0.0398 | 0.0000 |   0.0000 | 0.2950 |   0.1660 | -0.2950 |  -0.1660 |
|------------------------------------------------------------------------------------------|

Additionally, estimated power flows at branches can be examined:

printBranchData(system, analysis; show)

|------------------------------------------------------------------------|
| Branch Data                                                            |
|------------------------------------------------------------------------|
|            Label             |   From-Bus Power   |    To-Bus Power    |
|                              |                    |                    |
|   Branch | From-Bus | To-Bus |  Active | Reactive |  Active | Reactive |
|          |          |        |    [pu] |     [pu] |    [pu] |     [pu] |
|----------|----------|--------|---------|----------|---------|----------|
| Branch 1 | Bus 1    | Bus 2  |  0.5827 |   0.1360 | -0.5756 |  -0.1181 |
| Branch 2 | Bus 2    | Bus 3  |  0.2192 |  -0.0467 | -0.2181 |   0.0586 |
| Branch 3 | Bus 3    | Bus 4  | -0.2599 |  -0.0197 |  0.2613 |   0.0328 |
| Branch 4 | Bus 4    | Bus 5  |  0.0747 |   0.1199 | -0.0743 |  -0.1165 |
| Branch 5 | Bus 5    | Bus 6  |  0.0020 |   0.0496 | -0.0019 |  -0.0495 |
| Branch 6 | Bus 1    | Bus 6  |  0.1945 |   0.0989 | -0.1935 |  -0.0889 |
| Branch 7 | Bus 2    | Bus 6  |  0.0997 |   0.0290 | -0.0995 |  -0.0276 |
| Branch 8 | Bus 5    | Bus 2  | -0.0397 |  -0.0081 |  0.0397 |   0.0087 |
|------------------------------------------------------------------------|

Modifying Measurement Data

Measurement values and variances are now updated. Instead of recreating the measurement set and the PMU state estimation model from the beginning, both are modified simultaneously:

updatePmu!(system, device, analysis; label = "From Branch 8", magnitude = 1.1)
updatePmu!(system, device, analysis; label = "From Branch 2", angle = 0.2, noise = true)

These updates demonstrate the flexibility of JuliaGrid in modifying measurements. For the phasor measurement at From Branch 8, only the voltage magnitude is changed, while the angle measurement remains the same. For From Branch 2, only the angle value is updated by adding white Gaussian noise, while the magnitude remains unchanged.

Next, the PMU state estimation is solved again to compute the new estimate:

solve!(system, analysis)
power!(system, analysis)

Bus-related data can now be examined:

printBusData(system, analysis; show)

|-------------------------------------------------------------------------------------------|
| Bus Data                                                                                  |
|-------------------------------------------------------------------------------------------|
| Label |       Voltage       |  Power Generation  |   Power Demand    |  Power Injection   |
|       |                     |                    |                   |                    |
|   Bus | Magnitude |   Angle |  Active | Reactive | Active | Reactive |  Active | Reactive |
|       |      [pu] |   [rad] |    [pu] |     [pu] |   [pu] |     [pu] |    [pu] |     [pu] |
|-------|-----------|---------|---------|----------|--------|----------|---------|----------|
| Bus 1 |    1.0024 |  0.0184 |  0.9477 |   0.2634 | 0.0000 |   0.0000 |  0.9477 |   0.2634 |
| Bus 2 |    0.9838 | -0.0081 |  0.3811 |  -0.0038 | 0.2170 |   0.1270 |  0.1641 |  -0.1308 |
| Bus 3 |    0.9915 | -0.0609 |  0.1923 |  -0.0093 | 0.4780 |  -0.0390 | -0.2857 |   0.0297 |
| Bus 4 |    0.9985 | -0.0486 |  0.2200 |   0.1640 | 0.0760 |   0.0160 |  0.1440 |   0.1480 |
| Bus 5 |    0.9768 | -0.0595 | -0.1058 |   0.0308 | 0.1120 |   0.0750 | -0.2178 |  -0.0442 |
| Bus 6 |    0.9735 | -0.0581 | -0.4421 |  -0.0083 | 0.2950 |   0.1660 | -0.7371 |  -0.1743 |
|-------------------------------------------------------------------------------------------|

With the updated measurement values, the estimated results deviate more significantly from the exact values obtained through AC power flow, as the modified measurements no longer align with them.

Modifying Measurement Set

Excluding phasor measurements from the set, when obtained using optimal placement, should be done with caution, as it can easily render the system unobservable. In this example, two measurements will be taken out-of-service, and two additional measurements will be immediately included to maintain observability:

updatePmu!(system, device; label = "From Branch 2", status = 0)
updatePmu!(system, device; label = "From Branch 8", status = 0)

addPmu!(system, device; to = "Branch 2", magnitude = 0.2282, angle = -2.9587)
addPmu!(system, device; to = "Branch 8", magnitude = 0.0414, angle = -0.2424)

Since new measurements are being added, analysis is not passed to these functions. Directly modifying the existing PMU state estimation model is not possible in this case. To achieve this, users should define the new measurements beforehand with status = 0 and then activate them by setting status = 1.

Next, the PMU state estimation model is created and solved:

analysis = pmuStateEstimation(system, device)
solve!(system, analysis)
power!(system, analysis)

Bus-related data can now be examined:

printBusData(system, analysis; show)

|-------------------------------------------------------------------------------------------|
| Bus Data                                                                                  |
|-------------------------------------------------------------------------------------------|
| Label |       Voltage       |  Power Generation  |   Power Demand    |  Power Injection   |
|       |                     |                    |                   |                    |
|   Bus | Magnitude |   Angle |  Active | Reactive | Active | Reactive |  Active | Reactive |
|       |      [pu] |   [rad] |    [pu] |     [pu] |   [pu] |     [pu] |    [pu] |     [pu] |
|-------|-----------|---------|---------|----------|--------|----------|---------|----------|
| Bus 1 |    1.0000 |  0.0001 |  0.7772 |   0.2349 | 0.0000 |   0.0000 |  0.7772 |   0.2349 |
| Bus 2 |    0.9819 | -0.0269 |  0.0000 |  -0.0000 | 0.2170 |   0.1270 | -0.2170 |  -0.1270 |
| Bus 3 |    0.9897 | -0.0797 |  0.0000 |  -0.0000 | 0.4780 |  -0.0390 | -0.4780 |   0.0390 |
| Bus 4 |    1.0000 | -0.0302 |  0.4120 |   0.1687 | 0.0760 |   0.0160 |  0.3360 |   0.1527 |
| Bus 5 |    0.9781 | -0.0407 |  0.0000 |  -0.0000 | 0.1120 |   0.0750 | -0.1120 |  -0.0750 |
| Bus 6 |    0.9761 | -0.0398 | -0.0000 |  -0.0000 | 0.2950 |   0.1660 | -0.2950 |  -0.1660 |
|-------------------------------------------------------------------------------------------|

By taking certain measurements out-of-service, the estimation was affected. Adding more precise measurements while maintaining observability led to an estimation that more accurately reflects the exact power system state.