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, each measuring active and reactive power injection as well as bus voltage magnitude. Additionally, Meter 4 and Meter 5 measure active and reactive power flows along with current magnitudes.
Users can download a Julia script containing the scenarios from this section using the following link.
We define the power system, specify 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 angle of the slack bus (type = 3) remains fixed at the specified value.
AC State Estimation Wrapper Function
Throughout the simulations below, AC state estimation is run multiple times. To avoid repeatedly calling multiple JuliaGrid built-in functions, we define a wrapper function that performs the state estimation, allowing us to call a single function each time. This wrapper function computes estimate of bus voltage magnitudes and angles. Once the algorithm converges, it then calculates the powers at buses and branches:
function acStateEstimation!(system::PowerSystem, analysis::ACStateEstimation)
for iteration = 1:20
stopping = solve!(system, analysis)
if stopping < 1e-8
println("The algorithm converged in $iteration iterations.")
break
end
end
power!(system, analysis)
endDisplay Data Settings
Before running simulations, we configure the data display settings:
show = Dict("Shunt Power" => false, "Status" => false)Measurement Model
The first question is how to obtain measurement data, specifically values. These can either be predefined or generated artificially. This example explores different approaches for defining measurement data.
One way to obtain measurement data is by solving the power system to determine exact electrical quantities, which will serve as the source for measurements. More precisely, AC power flow analysis is used to compute voltages and powers, and these exact values are then used to generate measurement data.
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 2.
Figure 2: The 3-bus power system.
Data for loads and generators is added as follows:
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, 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)
endThe algorithm converged in 3 iterations.Bus Voltage Magnitude Measurements
First, let us examine the results obtained from AC power flow, as these voltage magnitude values will be used to generate measurements:
printBusData(system, 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 by creating all bus voltage magnitude measurements using a function that directly utilizes the results from AC power flow:
@voltmeter(label = "Meter ?")
addVoltmeter!(system, device, powerFlow; variance = 1e-4, noise = true)By setting noise = true, white Gaussian noise with variance = 1e-4 is added to the exact bus voltage magnitude values to generate measurement values:
printVoltmeterData(system, device)|-------------------------------------------|
| Voltmeter Data |
|-------------------------------------------|
| Label | Voltage Magnitude |
| | |
| | Measurement | Variance | Status |
| | [pu] | [pu] | |
|---------|-------------|----------|--------|
| Meter 1 | 0.9995 | 1.00e-04 | 1 |
| Meter 2 | 0.9429 | 1.00e-04 | 1 |
| Meter 3 | 0.9370 | 1.00e-04 | 1 |
|-------------------------------------------|Active and Reactive Power Measurements
Active and reactive power injection measurements can be created using the same method as for voltage magnitude measurements, where we first utilize the function power!. However, in this case, we take a different approach, adding measurements one by one to Meter 1, Meter 2, and Meter 3 using data from AC power flow:
@wattmeter(label = "Meter ?")
@varmeter(label = "Meter ?")
for (label, idx) in system.bus.label
Pᵢ, Qᵢ = injectionPower(system, powerFlow; label)
addWattmeter!(system, device; bus = label, active = Pᵢ, variance = 1e-3, noise = true)
addVarmeter!(system, device; bus = label, reactive = Qᵢ, variance = 1e-4, noise = true)
endNext, we define active and reactive power flow measurements at Branch 1 for both the from-bus and to-bus ends. Again, we use data from AC power flow, computing the powers and using exact values to generate measurements. The default setting noise = false remains unchanged, meaning exact values are used:
Pᵢⱼ, Qᵢⱼ = fromPower(system, powerFlow; label = "Branch 1")
addWattmeter!(system, device; label = "Meter 4", from = "Branch 1", active = Pᵢⱼ)
addVarmeter!(system, device; label = "Meter 4", from = "Branch 1", reactive = Qᵢⱼ)
Pⱼᵢ, Qⱼᵢ = toPower(system, powerFlow; label = "Branch 1")
addWattmeter!(system, device; label = "Meter 5", to = "Branch 1", active = Pⱼᵢ)
addVarmeter!(system, device; label = "Meter 5", to = "Branch 1", reactive = Qⱼᵢ)As a result, we obtain the set of active power measurements:
printWattmeterData(system, device)|-------------------------------------------|
| Wattmeter Data |
|-------------------------------------------|
| Label | Active Power |
| | |
| | Measurement | Variance | Status |
| | [pu] | [pu] | |
|---------|-------------|----------|--------|
| Meter 1 | 3.6529 | 1.00e-03 | 1 |
| Meter 2 | -1.0897 | 1.00e-03 | 1 |
| Meter 3 | -2.2774 | 1.00e-03 | 1 |
| Meter 4 | 2.2734 | 1.00e-04 | 1 |
| Meter 5 | -2.1610 | 1.00e-04 | 1 |
|-------------------------------------------|The set of reactive power measurements can be viewed as follows:
printVarmeterData(system, device)|-------------------------------------------|
| Varmeter Data |
|-------------------------------------------|
| Label | Reactive Power |
| | |
| | Measurement | Variance | Status |
| | [pu] | [pu] | |
|---------|-------------|----------|--------|
| Meter 1 | 0.7230 | 1.00e-04 | 1 |
| Meter 2 | -0.3113 | 1.00e-04 | 1 |
| Meter 3 | -0.1950 | 1.00e-04 | 1 |
| Meter 4 | 0.6706 | 1.00e-04 | 1 |
| Meter 5 | -0.5021 | 1.00e-04 | 1 |
|-------------------------------------------|Current Magnitude Measurements
Finally, current magnitude measurements need to be defined at Branch 1 for both the from-bus and to-bus ends. Here, we assume these values are known in advance and provide them directly to the functions:
@ammeter(statusFrom = 0, statusTo = 0)
addAmmeter!(system, device; label = "Meter 4", from = "Branch 1", magnitude = 1.36)
addAmmeter!(system, device; 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. However, even though they do not impact the computation, they still occupy positions in the matrices and vectors used for state estimation. If there is no plan to put these measurements in-service later, it is advisable to exclude them from the measurement set.
The obtained set of current magnitude measurements is:
printAmmeterData(system, device)|-------------------------------------------|
| 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 obtaining the measurements, we proceed with solving the AC state estimation by initializing the Gauss-Newton method to compute the WLS estimator:
analysis = gaussNewton(system, device)Next, we run the iterative algorithm to estimate bus voltages and use these results to compute power values:
acStateEstimation!(system, analysis)The algorithm converged in 4 iterations.This allows users to observe estimated bus voltages along with power values associated with buses:
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.0037 | 0.0000 | 3.6154 | 0.7269 | 0.0000 | 0.0000 | 3.6154 | 0.7269 |
| Bus 2 | 0.9359 | -0.0397 | -0.0155 | -0.0086 | 1.1000 | 0.3000 | -1.1155 | -0.3086 |
| Bus 3 | 0.9400 | -0.0582 | 0.0061 | 0.0065 | 2.3000 | 0.2000 | -2.2939 | -0.1935 |
|-------------------------------------------------------------------------------------------|Power flows at branches can also be examined:
printBranchData(system, 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.2745 | 0.6686 | -2.1629 | -0.5012 | 0.1116 | 0.1674 |
| Branch 2 | Bus 1 | Bus 2 | 1.3409 | 0.0584 | -1.2515 | -0.0047 | 0.0894 | 0.0536 |
| Branch 3 | Bus 2 | Bus 3 | 0.1360 | -0.3039 | -0.1309 | 0.3077 | 0.0051 | 0.0038 |
|-------------------------------------------------------------------------------------------|Additionally, users can retrieve results related to measurement devices. For instance, estimated values and corresponding residuals for wattmeters can be displayed:
printWattmeterData(system, device, analysis)|-----------------------------------------------------------------|
| Wattmeter Data |
|-----------------------------------------------------------------|
| Label | Active Power |
| | |
| | Measurement | Variance | Estimate | Residual | Status |
| | [pu] | [pu] | [pu] | [pu] | |
|---------|-------------|----------|----------|----------|--------|
| Meter 1 | 3.6529 | 1.00e-03 | 3.6154 | 0.0375 | 1 |
| Meter 2 | -1.0897 | 1.00e-03 | -1.1155 | 0.0258 | 1 |
| Meter 3 | -2.2774 | 1.00e-03 | -2.2939 | 0.0165 | 1 |
| Meter 4 | 2.2734 | 1.00e-04 | 2.2745 | -0.0011 | 1 |
| Meter 5 | -2.1610 | 1.00e-04 | -2.1629 | 0.0019 | 1 |
|-----------------------------------------------------------------|Modifying Measurement Data
Measurement values and variances will now be updated. Instead of recreating the measurement set and the Gauss-Newton method from the beginning, both models will be modified simultaneously:
updateVoltmeter!(system, device, analysis; label = "Meter 1", magnitude = 1.0, noise = false)
updateWattmeter!(system, device, analysis; label = "Meter 2", active = -1.1, variance = 1e-6)
updateVarmeter!(system, device, analysis; label = "Meter 3", variance = 1e-1)These updates demonstrate the flexibility of JuliaGrid in modifying measurements. For the voltmeter, we now generate a measurement without adding noise, for the wattmeter, we change both the value and variance, while the varmeter retains its previous value with only a variance adjustment.
Next, AC state estimation is run again to compute the new estimate without recreating the Gauss-Newton model. Additionally, this step initializes the iterative algorithm with a warm start, as the initial voltage magnitudes and angles correspond to the solution from the base case analysis:
acStateEstimation!(system, analysis)The algorithm converged in 3 iterations.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.0035 | 0.0000 | 3.6077 | 0.7265 | 0.0000 | 0.0000 | 3.6077 | 0.7265 |
| Bus 2 | 0.9361 | -0.0395 | -0.0000 | -0.0096 | 1.1000 | 0.3000 | -1.1000 | -0.3096 |
| Bus 3 | 0.9398 | -0.0582 | -0.0025 | 0.0074 | 2.3000 | 0.2000 | -2.3025 | -0.1926 |
|-------------------------------------------------------------------------------------------|Modifying Measurement Set
Now, we modify the measurement set by including current magnitude measurements:
updateAmmeter!(system, device, analysis; label = "Meter 4", status = 1)
updateAmmeter!(system, device, analysis; label = "Meter 5", status = 1)We then solve the AC state estimation again:
acStateEstimation!(system, analysis)The algorithm converged in 8 iterations.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.1373 | 0.0000 | 3.4577 | 0.6979 | 0.0000 | 0.0000 | 3.4577 | 0.6979 |
| Bus 2 | 1.0795 | -0.0290 | -0.0001 | -0.0258 | 1.1000 | 0.3000 | -1.1001 | -0.3258 |
| Bus 3 | 1.0834 | -0.0425 | 0.0894 | -0.0122 | 2.3000 | 0.2000 | -2.2106 | -0.2122 |
|-------------------------------------------------------------------------------------------|The bus voltage estimates appear suspicious, indicating the presence of bad data among the newly added ammeter measurements. To address this, we perform bad data processing:
outlier = residualTest!(system, device, analysis; threshold = 4.0)An outlier with a significantly high normalized residual is detected, specifically related to the current magnitude measurements in Meter 4:
julia> outlier.detecttruejulia> outlier.maxNormalizedResidual78.72321756657138
The bad data processing function automatically removes the detected outlier. Before repeating the AC state estimation, using a warm start is not advisable, as the previous state was obtained in the presence of bad data. Instead, it is useful to reset the initial point, for example, by using the values defined within the power system data:
setInitialPoint!(system, analysis)
acStateEstimation!(system, analysis)The algorithm converged in 5 iterations.Observing the bus-related data, we can confirm that the results now align with expectations, as the measurement set is free from bad data:
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.0018 | 0.0000 | 3.6095 | 0.7269 | 0.0000 | 0.0000 | 3.6095 | 0.7269 |
| Bus 2 | 0.9343 | -0.0396 | -0.0000 | -0.0094 | 1.1000 | 0.3000 | -1.1000 | -0.3094 |
| Bus 3 | 0.9380 | -0.0585 | -0.0034 | 0.0078 | 2.3000 | 0.2000 | -2.3034 | -0.1922 |
|-------------------------------------------------------------------------------------------|