AC State Estimation
In this example, we analyze a 3-bus power system, shown in Figure 1. The objective is to estimate bus voltage magnitudes and angles for a given measurement configuration.
Figure 1: The 3-bus power system with the given measurement configuration.
The measurement set consists of Meter 1, Meter 2, and Meter 3, which measure active and reactive power injections and bus voltage magnitudes. Meter 4 and Meter 5 measure active and reactive power flows and current magnitudes.
Users can download a Julia script containing the scenarios from this section using the following link.
We define the power system, add buses and branches, and assign the generator to the slack bus:
system = powerSystem()
addBus!(system; label = "Bus 1", magnitude = 1.01, angle = 0.0, type = 3)
addBus!(system; label = "Bus 2", magnitude = 0.92, angle = -0.04)
addBus!(system; label = "Bus 3", magnitude = 0.93, angle = -0.05)
@branch(reactance = 0.03)
addBranch!(system; label = "Branch 1", from = "Bus 1", to = "Bus 3", resistance = 0.02)
addBranch!(system; label = "Branch 2", from = "Bus 1", to = "Bus 2", resistance = 0.05)
addBranch!(system; label = "Branch 3", from = "Bus 2", to = "Bus 3", resistance = 0.04)
addGenerator!(system; label = "Generator 1", bus = "Bus 1")The magnitude and angle values define the initial point for the iterative AC state estimation algorithm, while the slack bus angle (type = 3) remains fixed.
Display Data Settings
Before running simulations, we set verbose from its default silent mode (0) to basic output (1):
@config(verbose = 1)For more detailed solver output, verbose can also be adjusted in the functions that solve specific analyses. Next, we configure the data display settings:
show = Dict("Shunt Power" => false, "Status" => false)Measurement Model
The first question is how to obtain measurement values. They can be predefined or generated artificially. This example explores both approaches.
One way to obtain measurement values is to solve the power system and use the resulting electrical quantities as measurement sources. Here, AC power flow computes voltages and powers, and these exact values are then used to generate measurements.
We begin by initializing the measurement variable:
monitoring = measurement(system)AC Power Flow
AC power flow analysis requires generator and load data. Figure 2 shows the system configuration.
Figure 2: The 3-bus power system.
We add the load and generator data:
updateBus!(system; label = "Bus 2", type = 1, active = 1.1, reactive = 0.3)
updateBus!(system; label = "Bus 3", type = 1, active = 2.3, reactive = 0.2)
updateGenerator!(system; label = "Generator 1", active = 3.3, reactive = 2.1)Next, we run AC power flow analysis to obtain bus voltages:
acModel!(system)
powerFlow = newtonRaphson(system)
powerFlow!(powerFlow)EXIT: The solution was found using the Newton-Raphson method in 3 iterations.Bus Voltage Magnitude Measurements
First, we inspect the AC power flow results because the voltage magnitudes will be used to generate measurements:
printBusData(powerFlow)|-----------------------------|
| Bus Data |
|-----------------------------|
| Label | Voltage |
| | |
| Bus | Magnitude | Angle |
| | [pu] | [rad] |
|-------|-----------|---------|
| Bus 1 | 1.0000 | 0.0000 |
| Bus 2 | 0.9324 | -0.0399 |
| Bus 3 | 0.9360 | -0.0586 |
|-----------------------------|One way to obtain measurements is to create all bus voltage magnitude measurements with a function that uses the AC power flow results directly:
@voltmeter(label = "Meter ?")
addVoltmeter!(monitoring, powerFlow; variance = 1e-4, noise = true)Setting noise = true adds white Gaussian noise with variance = 1e-4 to the exact bus voltage magnitudes to generate measurement values:
printVoltmeterData(monitoring)|-------------------------------------------|
| Voltmeter Data |
|-------------------------------------------|
| Label | Voltage Magnitude |
| | |
| | Measurement | Variance | Status |
| | [pu] | [pu] | |
|---------|-------------|----------|--------|
| Meter 1 | 0.9988 | 1.00e-04 | 1 |
| Meter 2 | 0.9158 | 1.00e-04 | 1 |
| Meter 3 | 0.9381 | 1.00e-04 | 1 |
|-------------------------------------------|Active and Reactive Power Measurements
Active and reactive power injection measurements can be created with the same method as voltage magnitude measurements, using power! first. Here, however, we add measurements one by one to Meter 1, Meter 2, and Meter 3 using AC power flow data:
@wattmeter(label = "Meter ?")
@varmeter(label = "Meter ?")
for (label, idx) in system.bus.label
Pᵢ, Qᵢ = injectionPower(powerFlow; label)
addWattmeter!(monitoring; bus = label, active = Pᵢ, variance = 1e-3, noise = true)
addVarmeter!(monitoring; bus = label, reactive = Qᵢ, variance = 1e-4, noise = true)
endNext, we define active and reactive power flow measurements at both ends of Branch 1. We again use AC power flow data, compute the powers, and use the exact values to generate measurements. The default setting noise = false remains unchanged, so exact values are used:
Pᵢⱼ, Qᵢⱼ = fromPower(powerFlow; label = "Branch 1")
addWattmeter!(monitoring; label = "Meter 4", from = "Branch 1", active = Pᵢⱼ)
addVarmeter!(monitoring; label = "Meter 4", from = "Branch 1", reactive = Qᵢⱼ)
Pⱼᵢ, Qⱼᵢ = toPower(powerFlow; label = "Branch 1")
addWattmeter!(monitoring; label = "Meter 5", to = "Branch 1", active = Pⱼᵢ)
addVarmeter!(monitoring; label = "Meter 5", to = "Branch 1", reactive = Qⱼᵢ)This gives the set of active power measurements:
printWattmeterData(monitoring)|-------------------------------------------|
| Wattmeter Data |
|-------------------------------------------|
| Label | Active Power |
| | |
| | Measurement | Variance | Status |
| | [pu] | [pu] | |
|---------|-------------|----------|--------|
| Meter 1 | 3.6565 | 1.00e-03 | 1 |
| Meter 2 | -1.1098 | 1.00e-03 | 1 |
| Meter 3 | -2.2504 | 1.00e-03 | 1 |
| Meter 4 | 2.2734 | 1.00e-04 | 1 |
| Meter 5 | -2.1610 | 1.00e-04 | 1 |
|-------------------------------------------|The reactive power measurements are:
printVarmeterData(monitoring)|-------------------------------------------|
| Varmeter Data |
|-------------------------------------------|
| Label | Reactive Power |
| | |
| | Measurement | Variance | Status |
| | [pu] | [pu] | |
|---------|-------------|----------|--------|
| Meter 1 | 0.7428 | 1.00e-04 | 1 |
| Meter 2 | -0.2844 | 1.00e-04 | 1 |
| Meter 3 | -0.1935 | 1.00e-04 | 1 |
| Meter 4 | 0.6706 | 1.00e-04 | 1 |
| Meter 5 | -0.5021 | 1.00e-04 | 1 |
|-------------------------------------------|Current Magnitude Measurements
Finally, we define current magnitude measurements at both ends of Branch 1. Here, we assume these values are known in advance and pass them directly to the functions:
@ammeter(statusFrom = 0, statusTo = 0)
addAmmeter!(monitoring; label = "Meter 4", from = "Branch 1", magnitude = 1.36)
addAmmeter!(monitoring; label = "Meter 5", to = "Branch 1", magnitude = 2.37)For current magnitude measurements, we set statusFrom = 0 and statusTo = 0, indicating that these measurements are out-of-service and do not influence state estimation. Although they do not affect the computation, they still occupy positions in the matrices and vectors used for state estimation. If they will not be put in-service later, they should be excluded from the measurement set.
The current magnitude measurements are:
printAmmeterData(monitoring)|-------------------------------------------|
| Ammeter Data |
|-------------------------------------------|
| Label | Current Magnitude |
| | |
| | Measurement | Variance | Status |
| | [pu] | [pu] | |
|---------|-------------|----------|--------|
| Meter 4 | 1.3600 | 1.00e-04 | 0 |
| Meter 5 | 2.3700 | 1.00e-04 | 0 |
|-------------------------------------------|Base Case Analysis
After collecting the measurements, we solve the AC state estimation problem with the Gauss-Newton method to estimate bus voltages. The resulting estimates are then used to compute powers:
analysis = gaussNewton(monitoring)
stateEstimation!(analysis; power = true, verbose = 2)Number of entries in the Jacobian: 63
Number of measurement functions: 15
Number of state variables: 5
Number of buses: 3
Number of branches: 3
-----------------------------------------------
Iteration Objective Value Maximum Increment
-----------------------------------------------
0 1.27888280e+04 1.82876175e-02
1 1.11196852e+01 2.00720981e-03
2 8.40332990e+00 1.56669885e-05
3 8.40331530e+00 9.32655037e-08
4 8.40331529e+00 6.85528103e-10
Measurement Maximum Value
Absolute Residual: 7 3.3066e-02
Objective Value: 11 1.7445e+00
EXIT: The solution was found using the Gauss-Newton method in 4 iterations.We can then inspect the estimated bus voltages and bus power values:
printBusData(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 | 0.9933 | 0.0000 | 3.6314 | 0.7296 | 0.0000 | 0.0000 | 3.6314 | 0.7296 |
| Bus 2 | 0.9242 | -0.0413 | -0.0429 | 0.0041 | 1.1000 | 0.3000 | -1.1429 | -0.2959 |
| Bus 3 | 0.9288 | -0.0594 | 0.0238 | -0.0029 | 2.3000 | 0.2000 | -2.2762 | -0.2029 |
|-------------------------------------------------------------------------------------------|We can also inspect the branch power flows:
printBranchData(analysis; show)|-------------------------------------------------------------------------------------------|
| Branch Data |
|-------------------------------------------------------------------------------------------|
| Label | From-Bus Power | To-Bus Power | Series Power |
| | | | |
| Branch | From-Bus | To-Bus | Active | Reactive | Active | Reactive | Active | Reactive |
| | | | [pu] | [pu] | [pu] | [pu] | [pu] | [pu] |
|----------|----------|--------|--------|----------|---------|----------|--------|----------|
| Branch 1 | Bus 1 | Bus 3 | 2.2750 | 0.6747 | -2.1609 | -0.5035 | 0.1141 | 0.1712 |
| Branch 2 | Bus 1 | Bus 2 | 1.3564 | 0.0549 | -1.2630 | 0.0011 | 0.0934 | 0.0560 |
| Branch 3 | Bus 2 | Bus 3 | 0.1201 | -0.2970 | -0.1153 | 0.3006 | 0.0048 | 0.0036 |
|-------------------------------------------------------------------------------------------|Users can also retrieve measurement-device results. For example, wattmeter estimates and residuals can be displayed:
printWattmeterData(analysis)|-----------------------------------------------------------------|
| Wattmeter Data |
|-----------------------------------------------------------------|
| Label | Active Power |
| | |
| | Measurement | Variance | Estimate | Residual | Status |
| | [pu] | [pu] | [pu] | [pu] | |
|---------|-------------|----------|----------|----------|--------|
| Meter 1 | 3.6565 | 1.00e-03 | 3.6314 | 0.0251 | 1 |
| Meter 2 | -1.1098 | 1.00e-03 | -1.1429 | 0.0331 | 1 |
| Meter 3 | -2.2504 | 1.00e-03 | -2.2762 | 0.0258 | 1 |
| Meter 4 | 2.2734 | 1.00e-04 | 2.2750 | -0.0016 | 1 |
| Meter 5 | -2.1610 | 1.00e-04 | -2.1609 | -0.0001 | 1 |
|-----------------------------------------------------------------|Modifying Measurement Data
We now update measurement values and variances. Instead of recreating the measurement set and the Gauss-Newton method, we modify both models together:
updateVoltmeter!(analysis; label = "Meter 1", magnitude = 1.0, noise = false)
updateWattmeter!(analysis; label = "Meter 2", active = -1.1, variance = 1e-6)
updateVarmeter!(analysis; label = "Meter 3", variance = 1e-1)These updates demonstrate JuliaGrid’s flexibility in modifying measurements. The voltmeter measurement is generated without noise, the wattmeter value and variance are changed, and the varmeter keeps its previous value with only a variance adjustment.
Next, we run AC state estimation again without recreating the Gauss-Newton model. This enables a warm start because the initial voltage magnitudes and angles come from the base case solution:
stateEstimation!(analysis; power = true)EXIT: The solution was found using the Gauss-Newton method in 3 iterations.We can now inspect the bus data:
printBusData(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 | 0.9939 | 0.0000 | 3.6089 | 0.7321 | 0.0000 | 0.0000 | 3.6089 | 0.7321 |
| Bus 2 | 0.9259 | -0.0406 | -0.0001 | 0.0101 | 1.1000 | 0.3000 | -1.1001 | -0.2899 |
| Bus 3 | 0.9293 | -0.0592 | 0.0005 | -0.0133 | 2.3000 | 0.2000 | -2.2995 | -0.2133 |
|-------------------------------------------------------------------------------------------|Modifying Measurement Set
We now modify the measurement set by including current magnitude measurements:
updateAmmeter!(analysis; label = "Meter 4", status = 1)
updateAmmeter!(analysis; label = "Meter 5", status = 1)We then solve the AC state estimation problem again:
stateEstimation!(analysis; power = true)EXIT: The solution was found using the Gauss-Newton method in 8 iterations.We can now inspect the bus data:
printBusData(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.1336 | 0.0000 | 3.4520 | 0.7016 | 0.0000 | 0.0000 | 3.4520 | 0.7016 |
| Bus 2 | 1.0758 | -0.0294 | -0.0001 | -0.0069 | 1.1000 | 0.3000 | -1.1001 | -0.3069 |
| Bus 3 | 1.0793 | -0.0425 | 0.0955 | -0.0343 | 2.3000 | 0.2000 | -2.2045 | -0.2343 |
|-------------------------------------------------------------------------------------------|The bus voltage estimates appear suspicious, indicating bad data among the newly added ammeter measurements. To identify it, we perform bad data analysis:
outlier = residualTest!(analysis; threshold = 4.0)The analysis detects an outlier with a high normalized residual, associated with the current magnitude measurement in Meter 4:
julia> outlier.detecttruejulia> outlier.maxNormalizedResidual79.34131378048356
The bad data analysis function automatically removes the detected outlier. Before repeating AC state estimation, a warm start is not advisable because the previous state was computed with bad data. Instead, we reset the initial point, for example, using the values defined in the power system data:
setInitialPoint!(analysis)
stateEstimation!(analysis; power = true)EXIT: The solution was found using the Gauss-Newton method in 4 iterations.The bus data confirm that the results now align with expectations because the measurement set is free from bad data:
printBusData(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 | 0.9979 | 0.0000 | 3.6048 | 0.7312 | 0.0000 | 0.0000 | 3.6048 | 0.7312 |
| Bus 2 | 0.9302 | -0.0402 | -0.0001 | 0.0096 | 1.1000 | 0.3000 | -1.1001 | -0.2904 |
| Bus 3 | 0.9336 | -0.0586 | 0.0025 | -0.0142 | 2.3000 | 0.2000 | -2.2975 | -0.2142 |
|-------------------------------------------------------------------------------------------|