DC State Estimation

In this example, we monitor a 6-bus power system, shown in Figure 1, and estimate bus voltage angles using DC state estimation.

Figure 1: The 6-bus power system.

Info

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

We start by defining the units for voltage angles, which will be used throughout this example:

@voltage(pu, deg)

Next, we define the power system by adding buses, branches, and generators with cost functions:

system = powerSystem()

addBus!(system; label = "Bus 1", type = 3)
addBus!(system; label = "Bus 2", type = 1, active = 0.217)
addBus!(system; label = "Bus 3", type = 1, active = 0.478)
addBus!(system; label = "Bus 4", type = 2, active = 0.076)
addBus!(system; label = "Bus 5", type = 1, active = 0.112)
addBus!(system; label = "Bus 6", type = 2, active = 0.295)

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", active = 0.8, maxActive = 2.3)
addGenerator!(system; label = "Generator 2", bus = "Bus 4", active = 0.4, maxActive = 2.3)

cost!(system; generator = "Generator 1", active = 2, polynomial = [1100.0; 500.0; 150.0])
cost!(system; generator = "Generator 2", active = 2, polynomial = [1500.0; 700.0; 140.0])

After defining the power system data, we generate the DC model, including the key matrices and vectors used in the analysis:

dcModel!(system)

Measurement Model

Next, we define the measurements. The question is how to obtain measurement values. In this example, synthetic measurements are generated using DC optimal power flow results.

To start, we initialize the measurement variable:

monitoring = measurement(system)

DC Optimal Power Flow

To obtain bus voltage angles, we solve the DC optimal power flow. Using these values, we compute the active powers associated with buses and branches:

powerFlow = dcOptimalPowerFlow(system, Ipopt.Optimizer)
powerFlow!(powerFlow; power = true, verbose = 1)

EXIT: The optimal solution was found.

Active Power Injection Measurements

We obtain active power injection measurements from the DC optimal power flow analysis:

printBusData(powerFlow)

|---------------------------------------------------------------------|
| Bus Data                                                            |
|---------------------------------------------------------------------|
| Label | Voltage | Power Generation | Power Demand | Power Injection |
|       |         |                  |              |                 |
|   Bus |   Angle |           Active |       Active |          Active |
|       |   [deg] |             [pu] |         [pu] |            [pu] |
|-------|---------|------------------|--------------|-----------------|
| Bus 1 |  0.0000 |           0.7181 |       0.0000 |          0.7181 |
| Bus 2 | -1.5346 |           0.0000 |       0.2170 |         -0.2170 |
| Bus 3 | -4.1563 |           0.0000 |       0.4780 |         -0.4780 |
| Bus 4 | -1.1185 |           0.4599 |       0.0760 |          0.3839 |
| Bus 5 | -2.1399 |           0.0000 |       0.1120 |         -0.1120 |
| Bus 6 | -2.1948 |           0.0000 |       0.2950 |         -0.2950 |
|---------------------------------------------------------------------|

Next, we define these measurements:

@wattmeter(label = "Wattmeter ?")
for (label, idx) in system.bus.label
    Pᵢ = powerFlow.power.injection.active[idx]
    addWattmeter!(monitoring; bus = label, active = Pᵢ, variance = 1e-4, noise = true)
end

Setting noise = true adds white Gaussian noise with variance = 1e-4 to the exact values, generating the final measurement values.

Active Power Flow Measurements

Next, we include selected active power flow measurements using the DC optimal power flow results:

printBranchData(powerFlow)

|-----------------------------------------------------------------------|
| Branch Data                                                           |
|-----------------------------------------------------------------------|
|            Label             | From-Bus Power | To-Bus Power | Status |
|                              |                |              |        |
|   Branch | From-Bus | To-Bus |         Active |       Active |        |
|          |          |        |           [pu] |         [pu] |        |
|----------|----------|--------|----------------|--------------|--------|
| Branch 1 | Bus 1    | Bus 2  |         0.5357 |      -0.5357 |      1 |
| Branch 2 | Bus 2    | Bus 3  |         0.1989 |      -0.1989 |      1 |
| Branch 3 | Bus 3    | Bus 4  |        -0.2791 |       0.2791 |      1 |
| Branch 4 | Bus 4    | Bus 5  |         0.1049 |      -0.1049 |      1 |
| Branch 5 | Bus 5    | Bus 6  |         0.0239 |      -0.0239 |      1 |
| Branch 6 | Bus 1    | Bus 6  |         0.1824 |      -0.1824 |      1 |
| Branch 7 | Bus 2    | Bus 6  |         0.0886 |      -0.0886 |      1 |
| Branch 8 | Bus 5    | Bus 2  |        -0.0311 |       0.0311 |      1 |
|-----------------------------------------------------------------------|

We add two active power flow measurements:

addWattmeter!(monitoring; from = "Branch 1", active = powerFlow.power.from.active[1])
addWattmeter!(monitoring; from = "Branch 4", active = powerFlow.power.from.active[4])

Here, noise is not set, so the measurement values remain exact.

Active Power Measurements

Finally, we display the complete measurement set:

printWattmeterData(monitoring)

|-----------------------------------------------|
| Wattmeter Data                                |
|-----------------------------------------------|
|    Label    |          Active Power           |
|             |                                 |
|             | Measurement | Variance | Status |
|             |        [pu] |     [pu] |        |
|-------------|-------------|----------|--------|
| Wattmeter 1 |      0.7178 | 1.00e-04 |      1 |
| Wattmeter 2 |     -0.2017 | 1.00e-04 |      1 |
| Wattmeter 3 |     -0.4926 | 1.00e-04 |      1 |
| Wattmeter 4 |      0.3781 | 1.00e-04 |      1 |
| Wattmeter 5 |     -0.1272 | 1.00e-04 |      1 |
| Wattmeter 6 |     -0.2833 | 1.00e-04 |      1 |
| Wattmeter 7 |      0.5357 | 1.00e-04 |      1 |
| Wattmeter 8 |      0.1049 | 1.00e-04 |      1 |
|-----------------------------------------------|

Figure 2 shows the measurement configuration, including active power injection measurements at all buses and two active power flow measurements.

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

Base Case Analysis

After obtaining the measurements, we create the DC state estimation model:

analysis = dcStateEstimation(monitoring)

Next, we solve the model to determine the WLS estimator for bus voltage angles, then use the results to compute power values:

stateEstimation!(analysis; power = true, verbose = 1)

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

We can then inspect the estimated bus voltage angles and corresponding power values:

printBusData(analysis)

|---------------------------------------------------------------------|
| Bus Data                                                            |
|---------------------------------------------------------------------|
| Label | Voltage | Power Generation | Power Demand | Power Injection |
|       |         |                  |              |                 |
|   Bus |   Angle |           Active |       Active |          Active |
|       |   [deg] |             [pu] |         [pu] |            [pu] |
|-------|---------|------------------|--------------|-----------------|
| Bus 1 |  0.0000 |           0.7190 |       0.0000 |          0.7190 |
| Bus 2 | -1.5329 |           0.0161 |       0.2170 |         -0.2009 |
| Bus 3 | -4.2839 |          -0.0123 |       0.4780 |         -0.4903 |
| Bus 4 | -1.2193 |           0.4578 |       0.0760 |          0.3818 |
| Bus 5 | -2.1962 |          -0.0149 |       0.1120 |         -0.1269 |
| Bus 6 | -2.2133 |           0.0122 |       0.2950 |         -0.2828 |
|---------------------------------------------------------------------|

We can also inspect the measurement results:

printWattmeterData(analysis)

|---------------------------------------------------------------------|
| Wattmeter Data                                                      |
|---------------------------------------------------------------------|
|    Label    |                     Active Power                      |
|             |                                                       |
|             | Measurement | Variance | Estimate | Residual | Status |
|             |        [pu] |     [pu] |     [pu] |     [pu] |        |
|-------------|-------------|----------|----------|----------|--------|
| Wattmeter 1 |      0.7178 | 1.00e-04 |   0.7190 |  -0.0012 |      1 |
| Wattmeter 2 |     -0.2017 | 1.00e-04 |  -0.2009 |  -0.0008 |      1 |
| Wattmeter 3 |     -0.4926 | 1.00e-04 |  -0.4903 |  -0.0024 |      1 |
| Wattmeter 4 |      0.3781 | 1.00e-04 |   0.3818 |  -0.0037 |      1 |
| Wattmeter 5 |     -0.1272 | 1.00e-04 |  -0.1269 |  -0.0003 |      1 |
| Wattmeter 6 |     -0.2833 | 1.00e-04 |  -0.2828 |  -0.0005 |      1 |
| Wattmeter 7 |      0.5357 | 1.00e-04 |   0.5351 |   0.0006 |      1 |
| Wattmeter 8 |      0.1049 | 1.00e-04 |   0.1003 |   0.0046 |      1 |
|---------------------------------------------------------------------|

Modifying Measurement Data

We now modify the measurement values. Instead of recreating the measurement set and the DC state estimation model, we update both together:

updateWattmeter!(analysis; label = "Wattmeter 7", active = 1.1)
updateWattmeter!(analysis; label = "Wattmeter 8", active = 1.6)

Changing these measurement values introduces two outliers into the dataset, affecting the estimates.

Next, we solve the DC state estimation problem again to compute the updated estimate:

stateEstimation!(analysis; power = true, verbose = 1)

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

We can now inspect the bus data:

printBusData(analysis)

|---------------------------------------------------------------------|
| Bus Data                                                            |
|---------------------------------------------------------------------|
| Label | Voltage | Power Generation | Power Demand | Power Injection |
|       |         |                  |              |                 |
|   Bus |   Angle |           Active |       Active |          Active |
|       |   [deg] |             [pu] |         [pu] |            [pu] |
|-------|---------|------------------|--------------|-----------------|
| Bus 1 |  0.0000 |           0.8587 |       0.0000 |          0.8587 |
| Bus 2 | -1.7757 |          -0.2352 |       0.2170 |         -0.4522 |
| Bus 3 | -1.4966 |           0.1411 |       0.4780 |         -0.3369 |
| Bus 4 |  2.4015 |           0.9455 |       0.0760 |          0.8695 |
| Bus 5 | -2.5798 |          -0.3121 |       0.1120 |         -0.4241 |
| Bus 6 | -2.8745 |          -0.2200 |       0.2950 |         -0.5150 |
|---------------------------------------------------------------------|

With the modified values for Wattmeter 7 and Wattmeter 8, the estimates deviate more from the exact values computed by DC optimal power flow because the altered measurements no longer match their corresponding values.

We then compute the LAV estimator instead of the WLS estimator:

analysis = dcLavStateEstimation(monitoring, Ipopt.Optimizer)
stateEstimation!(analysis; power = true, verbose = 1)

EXIT: The optimal solution was found.

We can inspect the bus data:

printBusData(analysis)

|---------------------------------------------------------------------|
| Bus Data                                                            |
|---------------------------------------------------------------------|
| Label | Voltage | Power Generation | Power Demand | Power Injection |
|       |         |                  |              |                 |
|   Bus |   Angle |           Active |       Active |          Active |
|       |   [deg] |             [pu] |         [pu] |            [pu] |
|-------|---------|------------------|--------------|-----------------|
| Bus 1 |  0.0000 |           0.7267 |       0.0000 |          0.7267 |
| Bus 2 | -1.5490 |           0.0153 |       0.2170 |         -0.2017 |
| Bus 3 | -4.3389 |          -0.0146 |       0.4780 |         -0.4926 |
| Bus 4 | -1.2805 |           0.4541 |       0.0760 |          0.3781 |
| Bus 5 | -2.2270 |          -0.0152 |       0.1120 |         -0.1272 |
| Bus 6 | -2.2379 |           0.0117 |       0.2950 |         -0.2833 |
|---------------------------------------------------------------------|

The LAV estimates are closer to the exact DC optimal power flow values because LAV is more robust to outliers than WLS.

Modifying Measurement Set

We continue with the LAV state estimation model and set two measurements out-of-service:

updateWattmeter!(analysis; label = "Wattmeter 1", status = 0)
updateWattmeter!(analysis; label = "Wattmeter 5", status = 0)

We then recompute the LAV estimator and active power values:

stateEstimation!(analysis; power = true, verbose = 1)

EXIT: The optimal solution was found.

We can now inspect the bus data:

printBusData(analysis)

|---------------------------------------------------------------------|
| Bus Data                                                            |
|---------------------------------------------------------------------|
| Label | Voltage | Power Generation | Power Demand | Power Injection |
|       |         |                  |              |                 |
|   Bus |   Angle |           Active |       Active |          Active |
|       |   [deg] |             [pu] |         [pu] |            [pu] |
|-------|---------|------------------|--------------|-----------------|
| Bus 1 |  0.0000 |           1.5931 |       0.0000 |          1.5931 |
| Bus 2 | -3.1513 |           0.0153 |       0.2170 |         -0.2017 |
| Bus 3 | -7.2823 |          -0.0146 |       0.4780 |         -0.4926 |
| Bus 4 | -5.3317 |           0.4541 |       0.0760 |          0.3781 |
| Bus 5 | -7.2694 |          -0.8816 |       0.1120 |         -0.9936 |
| Bus 6 | -5.9327 |           0.0117 |       0.2950 |         -0.2833 |
|---------------------------------------------------------------------|

Although LAV is more robust than WLS when handling outliers, estimate accuracy still depends on factors such as outlier magnitude, outlier count, and meter placement within the power system. Removing two accurate measurements while keeping outliers in the system shows that even LAV cannot fully compensate for the loss of reliable data, leading to less accurate estimates.