AC State Estimation
To perform nonlinear or AC state estimation, the initial requirement is to have the PowerSystem composite type configured with the AC model, along with the Measurement composite type storing measurement data. Next, we can develop either the weighted least-squares (WLS) model, utilizing the Gauss-Newton method, or the least absolute value (LAV) model. These models are encapsulated within the ACStateEstimation type:
For resolving the AC state estimation problem and obtaining bus voltage magnitudes and angles, utilize the following function:
After executing the function solve!, where the user employs the Gauss-Newton method, the user has the ability to check if the measurement set contains outliers throughout bad data analysis and remove those measurements using:
Moreover, before the creating ACStateEstimation type, users can initiate observability analysis to identify observable islands and restore observability by employing:
After obtaining the AC state estimation solution, JuliaGrid offers post-processing analysis functions for calculating powers and currents associated with buses and branches:
Furthermore, there are specialized functions dedicated to calculating specific types of powers related to particular buses and branches:
Likewise, there are specialized functions dedicated to calculating specific types of currents related to particular buses or branches:
Bus Type Modification
In AC state estimation, it is necessary to designate a slack bus, where the bus voltage angle is known. Therefore, when establishing the ACStateEstimation type, the initially assigned slack bus is evaluated and may be altered. If the designated slack bus (type = 3) lacks a connected in-service generator, it will be changed to a demand bus (type = 1). Conversely, the first generator bus (type = 2) with an active in-service generator linked to it will be reassigned as the new slack bus (type = 3).
Observability Analysis
To initiate the power system with measurements at specific locations, follow the provided example:
system = powerSystem()
device = measurement()
addBus!(system; label = "Bus 1", type = 3)
addBus!(system; label = "Bus 2", type = 1, active = 0.1, reactive = 0.01)
addBus!(system; label = "Bus 3", type = 2, active = 0.5, reactive = 0.01)
addBus!(system; label = "Bus 4", type = 1, active = 0.2, reactive = 0.02)
addBus!(system; label = "Bus 5", type = 1, active = 0.3, reactive = 0.03)
addBus!(system; label = "Bus 6", type = 1, active = 0.1, reactive = 0.01)
addBus!(system; label = "Bus 7", type = 1, active = 0.1, reactive = 0.01)
@branch(resistance = 0.02, conductance = 1e-4, susceptance = 0.002)
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.01)
addBranch!(system; label = "Branch 3", from = "Bus 2", to = "Bus 5", reactance = 0.02)
addBranch!(system; label = "Branch 4", from = "Bus 3", to = "Bus 4", reactance = 0.03)
addBranch!(system; label = "Branch 5", from = "Bus 5", to = "Bus 6", reactance = 0.05)
addBranch!(system; label = "Branch 6", from = "Bus 3", to = "Bus 5", reactance = 0.05)
addBranch!(system; label = "Branch 7", from = "Bus 6", to = "Bus 7", reactance = 0.05)
addGenerator!(system; label = "Generator 1", bus = "Bus 1", active = 4.2, reactive = 0.2)
addGenerator!(system; label = "Generator 2", bus = "Bus 3", active = 0.2, reactive = 0.1)
addWattmeter!(system, device; label = "Wattmeter 1", from = "Branch 1", active = 1.15)
addVarmeter!(system, device; label = "Varmeter 1", from = "Branch 1", reactive = -0.50)
addWattmeter!(system, device; label = "Wattmeter 2", from = "Branch 4", active = 0.20)
addVarmeter!(system, device; label = "Varmeter 2", from = "Branch 4", reactive = -0.02)
addWattmeter!(system, device; label = "Wattmeter 3", from = "Branch 5", active = -0.20)
addVarmeter!(system, device; label = "Varmeter 3", from = "Branch 5", reactive = 0.02)
addWattmeter!(system, device; label = "Wattmeter 4", bus = "Bus 2", active = -0.1)
addVarmeter!(system, device; label = "Varmeter 4", bus = "Bus 2", reactive = -0.01)
addWattmeter!(system, device; label = "Wattmeter 5", bus = "Bus 3", active = -0.30)
addVarmeter!(system, device; label = "Varmeter 5", bus = "Bus 3", reactive = 0.66)Attempting to solve this system immediately may not be possible because the gain matrix will be singular. To avoid this situation, users can perform observability analysis. JuliaGrid employs standard observability analysis performed on the linear decoupled measurement model. Active power measurements from wattmeters are utilized to estimate bus voltage angles, while reactive power measurements from varmeters are used to estimate bus voltage magnitudes. This necessitates that measurements of active and reactive power come in pairs.
We suggest that readers refer to the tutorial on Observability Analysis for insights into the implementation.
However, the initial step involves defining observable islands. JuliaGrid offers users two options for obtaining observable islands: flow observable islands or maximal observable islands. The selection depends on the power system's structure and available measurements. Identifying only flow observable islands reduces complexity in the island detection function but increases complexity in the restoration function.
Flow Observable Islands
Now, let us identify flow observable islands:
islands = islandTopologicalFlow(system, device)As a result, four flow observable islands are identified: Bus 1 and Bus 2 form the first island, Bus 3 and Bus 4 form the second island, Bus 5 and Bus 6 constitute the third island, while Bus 7 form the fourth island:
julia> islands.island4-element Vector{Vector{Int64}}: [1, 2] [3, 4] [5, 6] [7]
Maximal Observable Islands
Following that, we will instruct the user on obtaining maximal observable islands:
islands = islandTopological(system, device)The outcome reveals the identification of two maximal observable islands:
julia> islands.island2-element Vector{Vector{Int64}}: [1, 2, 3, 4, 5, 6] [7]
It is evident that upon comparing this result with the flow islands, the merging of the two injection measurements at Bus 2 and Bus 3 consolidated the first, second, and third flow observable islands into a single island.
Restore Observability
Before commencing the restoration of observability in the context of the linear decoupled measurement model and observability analysis, it is imperative to ensure that the system possesses one bus voltage magnitude measurement. This necessity arises from the fact that observable islands are identified based on wattmeters, where wattmeters are tasked with estimating voltage angles. Since one voltage angle is already known from the slack bus, the same principle should be applied to bus voltage magnitudes. Therefore, to address this requirement, we add:
addVoltmeter!(system, device; bus = "Bus 1", magnitude = 1.0)Subsequently, the user needs to establish a set of pseudo-measurements, where measurements must come in pairs as well. Let us create that set:
pseudo = measurement()
addWattmeter!(system, pseudo; label = "Pseudo-Wattmeter 1", bus = "Bus 1", active = 0.31)
addVarmeter!(system, pseudo; label = "Pseudo-Varmeter 1", bus = "Bus 1", reactive = -0.19)
addWattmeter!(system, pseudo; label = "Pseudo-Wattmeter 2", from = "Branch 7", active = 0.10)
addVarmeter!(system, pseudo; label = "Pseudo-Varmeter 2", from = "Branch 7", reactive = 0.01)The labels for specific pseudo-measurements must differ from those defined in the measurements stored in the device set. This is necessary because the next step involves adding pseudo-measurements to the device set.
Subsequently, the user can execute the restorationGram! function:
restorationGram!(system, device, pseudo, islands)This function attempts to restore observability using pseudo-measurements. As a result, the inclusion of measurements from Pseudo-Wattmeter 2 and Pseudo-Varmeter 2 facilitates observability restoration, and these measurements are subsequently added to the device variable:
julia> device.wattmeter.labelOrderedCollections.OrderedDict{String, Int64} with 6 entries: "Wattmeter 1" => 1 "Wattmeter 2" => 2 "Wattmeter 3" => 3 "Wattmeter 4" => 4 "Wattmeter 5" => 5 "Pseudo-Wattmeter 2" => 6julia> device.varmeter.labelOrderedCollections.OrderedDict{String, Int64} with 6 entries: "Varmeter 1" => 1 "Varmeter 2" => 2 "Varmeter 3" => 3 "Varmeter 4" => 4 "Varmeter 5" => 5 "Pseudo-Varmeter 2" => 6
Consequently, the power system becomes observable, allowing the user to proceed with forming the AC state estimation model and solving it. Ensuring the observability of the system does not guarantee obtaining accurate estimates of the state variables. Numerical ill-conditioning may adversely impact the state estimation algorithm. However, in most cases, efficient estimation becomes feasible when the system is observable [1].
Additionally, it is worth mentioning that restoration might encounter difficulties due to the default zero pivot threshold set at 1e-5. This threshold can be modified using the restorationGram! function.
During the restoration step, if users define bus phasor measurements, these measurements will be considered. Consequently, the system may achieve observability even if multiple islands persist.
Weighted Least-squares Estimator
To begin, we will define the PowerSystem and Measurement types:
system = powerSystem()
device = measurement()
addBus!(system; label = "Bus 1", type = 3)
addBus!(system; label = "Bus 2", type = 1, active = 0.1, reactive = 0.01)
addBus!(system; label = "Bus 3", type = 1, active = 2.5, reactive = 0.2)
@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)
addGenerator!(system; label = "Generator 1", bus = "Bus 1", active = 3.2, reactive = 0.3)
@voltmeter(label = "Voltmeter ? (!)")
addVoltmeter!(system, device; bus = "Bus 1", magnitude = 1.0, variance = 1e-2)
@ammeter(label = "Ammeter ? (!)")
addAmmeter!(system, device; from = "Branch 3", magnitude = 0.947, variance = 1e-1)
addAmmeter!(system, device; to = "Branch 2", magnitude = 1.674, variance = 1e-1)
@wattmeter(label = "Wattmeter ? (!)")
addWattmeter!(system, device; from = "Branch 1", active = 1.046, variance = 1e-3)
addWattmeter!(system, device; bus = "Bus 2", active = -0.1, variance = 2e-3)
@varmeter(label = "Varmeter ? (!)")
addVarmeter!(system, device; from = "Branch 1", reactive = 0.059, variance = 1e-4)
addVarmeter!(system, device; bus = "Bus 2", reactive = -0.01, variance = 1e-3)Next, to establish the AC state estimation model, we will utilize the gaussNewton function:
analysis = gaussNewton(system, device)Here, the user triggers LU factorization as the default method for solving the system of linear equations within each iteration of the Gauss-Newton method. However, the user also has the option to select alternative factorization methods such as LDLt or QR:
analysis = gaussNewton(system, device, LDLt)Setup Starting Voltages
The initial voltages for the Gauss-Newton method are determined based on the specified initial voltage magnitudes and angles within the buses of the PowerSystem type. These values are then forwarded to the ACStateEstimation during the execution of the gaussNewton function. Therefore, the starting voltages in this example are as follows:
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 modify these vectors according to their own requirements in order to adjust the starting voltages. For instance, users can conduct an initial AC power flow analysis and utilize the obtained solution as the starting voltages for AC state estimation:
powerFlow = newtonRaphson(system)
for iteration = 1:10
mismatch!(system, powerFlow)
solve!(system, powerFlow)
end
for i = 1:system.bus.number
analysis.voltage.magnitude[i] = powerFlow.voltage.magnitude[i]
analysis.voltage.angle[i] = powerFlow.voltage.angle[i]
endState Estimator
To conduct an iterative process using the Gauss-Newton method, it is essential to include the solve! function inside the iteration loop. For example:
for iteration = 1:20
solve!(system, analysis)
endOnce the state estimator is obtained, users can access the bus voltage magnitudes and angles using:
julia> print(system.bus.label, analysis.voltage.magnitude, analysis.voltage.angle)Bus 1: 0.9999492012048501, 0.0 Bus 2: 0.9763971958413307, -0.05196984653510395 Bus 3: 0.9562624707681898, -0.08054149367597473
Breaking the Iterative Process
The iterative process can be terminated using the solve! function. The following code demonstrates how to utilize this function to break out of the iteration loop:
analysis = gaussNewton(system, device)
for iteration = 1:20
stopping = solve!(system, analysis)
if stopping < 1e-8
break
end
endThe solve! function returns the maximum absolute values of the state variable increment, which are commonly used as a convergence criterion in the iterative Gauss-Newton algorithm.
We suggest that readers refer to the tutorial on AC State Estimation for insights into the implementation.
Inclusion of PMUs in Rectangular Coordinates
In the example above, our focus is solely on solving the AC state estimation using SCADA measurements. However, users have the option to also integrate PMUs into the AC state estimation, either in the rectangular or polar coordinate system.
The default approach is to include PMUs in the rectangular coordinate system:
@pmu(label = "PMU ? (!)")
addPmu!(system, device; to = "Branch 1", magnitude = 1.05, angle = 3.047)In the case of the rectangular system, inclusion resolves ill-conditioned problems arising in polar coordinates due to small values of current magnitudes. However, this approach's main disadvantage is related to measurement errors, as measurement errors correspond to polar coordinates. Therefore, the covariance matrix must be transformed from polar to rectangular coordinates [4]. As a result, measurement errors of a single PMU are correlated, and the covariance matrix does not have a diagonal form. Despite that, the measurement error covariance matrix is usually considered as a diagonal matrix, affecting the accuracy of the state estimation.
In the example above, we specifically include PMUs where measurement error correlations are disregarded. This is evident through the precision matrix, which maintains a diagonal form:
julia> analysis = gaussNewton(system, device);julia> analysis.method.precision9×9 SparseArrays.SparseMatrixCSC{Float64, Int64} with 9 stored entries: 100.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 10.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 10.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 1000.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 500.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 10000.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 1000.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 99908.6 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 90778.2
Lastly, we incorporate correlation into our model by adding a new PMU with the desired error correlation:
addPmu!(system, device; bus = "Bus 3", magnitude = 0.95, angle = -0.08, correlated = true)Now, we can observe the precision matrix that does not hold a diagonal form:
julia> analysis = gaussNewton(system, device);julia> analysis.method.precision11×11 SparseArrays.SparseMatrixCSC{Float64, Int64} with 13 stored entries: 100.0 ⋅ ⋅ ⋅ ⋅ … ⋅ ⋅ ⋅ ⋅ 10.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 10.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 1000.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 500.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ … ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 90778.2 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 100069.0 860.583 ⋅ ⋅ ⋅ ⋅ ⋅ … ⋅ 860.583 1.10734e5
Inclusion of PMUs in Polar Coordinates
The second approach involves incorporating these measurements into the polar coordinate system. For instance:
addPmu!(system, device; from = "Branch 1", magnitude = 1.048, angle = -0.057, polar = true)This inclusion of PMUs provides more accurate state estimates compared to rectangular inclusion but demands longer computing time. PMUs are handled in the same manner as SCADA measurements. However, this approach is susceptible to ill-conditioned problems arising in polar coordinates due to small values of current magnitudes [2, 3].
It is important to note that with each individual phasor measurement, we can set the coordinate system, providing flexibility to include some in polar and some in rectangular systems.
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!(system, device; label = "PMU 2 (Bus 3)", correlated = false)Subsequently, by specifying the Orthogonal argument in the gaussNewton function, JuliaGrid implements a more robust approach to obtain the WLS estimator, which proves particularly beneficial when substantial differences exist among measurement variances:
analysis = gaussNewton(system, device, Orthogonal)
for iteration = 1:20
stopping = solve!(system, analysis)
if stopping < 1e-8
break
end
endBad Data Processing
After acquiring the WLS solution using the Gauss-Newton method, users can conduct bad data analysis employing the largest normalized residual test. Continuing with our defined power system and measurement set, let us introduce a new measurement. Upon proceeding to find the solution for this updated state:
addWattmeter!(system, device; from = "Branch 2", active = 31.1)
analysis = gaussNewton(system, device)
for iteration = 1:20
stopping = solve!(system, analysis)
if stopping < 1e-8
break
end
endHere, we can observe the impact of the outlier on the solution:
julia> print(system.bus.label, analysis.voltage.magnitude, analysis.voltage.angle)Bus 1: 1.105634031595695, 0.0 Bus 2: 1.0813928031976263, -0.04784530109433603 Bus 3: 1.0117796176799636, -0.16528213297101343
Following the solution acquisition, we can verify the presence of erroneous data. Detection of such data is determined by the threshold keyword. If the largest normalized residual's value exceeds the threshold, the measurement will be identified as bad data and consequently removed from the AC state estimation model:
outlier = residualTest!(system, device, analysis; threshold = 4.0)Users can examine the data obtained from the bad data analysis:
julia> outlier.detecttruejulia> outlier.maxNormalizedResidual282.7197102603639julia> outlier.label"Wattmeter 3 (From Branch 2)"
Hence, upon detecting bad data, the detect variable will hold true. The maxNormalizedResidual variable retains the value of the largest normalized residual, while the label contains the label of the measurement identified as bad data. JuliaGrid will mark the respective measurement as out-of-service within the Measurement type.
After removing bad data, a new estimate can be computed without considering this specific measurement. The user has the option to either restart the gaussNewton function or proceed directly to the iteration loop. However, if the latter option is chosen, using voltages obtained with outlier presence as the starting point could significantly impede algorithm convergence. To avoid this undesirable outcome, the user should first establish a new starting point and commence the iteration procedure. For instance:
for i = 1:system.bus.number
analysis.voltage.magnitude[i] = system.bus.voltage.magnitude[i]
analysis.voltage.angle[i] = system.bus.voltage.angle[i]
end
for iteration = 1:20
stopping = solve!(system, analysis)
if stopping < 1e-8
break
end
endConsequently, we obtain a new solution devoid of the impact of the outlier measurement:
julia> print(system.bus.label, analysis.voltage.magnitude, analysis.voltage.angle)Bus 1: 0.9938191289507472, 0.0 Bus 2: 0.9702353516839954, -0.05230134628742669 Bus 3: 0.9501354580385988, -0.08084132445043153
We suggest that readers refer to the tutorial on Bad Data Processing for insights into the implementation.
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 processing 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 = acLavStateEstimation(system, device, Ipopt.Optimizer)Setup Starting Primal Values
In JuliaGrid, the assignment of starting primal values for optimization variables takes place when the solve! function is executed. Starting primal values are determined based on the voltage fields within the ACStateEstimation type. By default, these values are initially established using the the initial bus voltage magnitudes and angles from PowerSystem 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 starting primal values when executing the solve! function.
Solution
To solve the formulated LAV state estimation model, simply execute the following function:
solve!(system, 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.9997142894573017, 0.0 Bus 2: 0.9761238285490208, -0.05198283480797272 Bus 3: 0.9559672017526282, -0.08055409509688616
We suggest that readers refer to the tutorial on Least Absolute Value Estimation for insights into the implementation.
Measurement Set Update
After establishing the Measurement composite type using the measurement function, users gain the capability to incorporate new measurement devices or update existing ones.
Once updates are completed, users can seamlessly progress towards generating the ACStateEstimation type using the gaussNewton or acLavStateEstimation function. Ultimately, resolving the AC state estimation is achieved through the utilization of the solve! function:
system = powerSystem()
device = measurement() # <- Initializing a Measurement instance
addBus!(system; label = "Bus 1", type = 3)
addBus!(system; label = "Bus 2", type = 1, active = 0.1, reactive = 0.01)
addBus!(system; label = "Bus 3", type = 1, active = 2.5, reactive = 0.2)
@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)
addGenerator!(system; label = "Generator 1", bus = "Bus 1", active = 3.2, reactive = 0.3)
@voltmeter(label = "Voltmeter ? (!)", variance = 1e-3)
addVoltmeter!(system, device; bus = "Bus 1", magnitude = 1.0)
@wattmeter(label = "Wattmeter ? (!)", varianceBus = 1e-2)
addWattmeter!(system, device; from = "Branch 1", active = 1.046)
addWattmeter!(system, device; bus = "Bus 2", active = -0.1)
@varmeter(label = "Varmeter ? (!)", varianceFrom = 1e-3)
addVarmeter!(system, device; from = "Branch 1", reactive = 0.059)
addVarmeter!(system, device; bus = "Bus 2", reactive = -0.01)
@pmu(label = "PMU ? (!)")
addPmu!(system, device; bus = "Bus 2", magnitude = 0.976, angle = -0.052)
analysis = gaussNewton(system, device) # <- Creating ACStateEstimation for the defined model
for iteration = 1:20
stopping = solve!(system, analysis)
if stopping < 1e-8
break
end
end
addWattmeter!(system, device; from = "Branch 3", active = 0.924)
updateWattmeter!(system, device; label = "Wattmeter 2 (Bus 2)", variance = 1e-4)
addVarmeter!(system, device; to = "Branch 3", reactive = -0.044, variance = 1e-5)
updateVarmeter!(system, device; label = "Varmeter 2 (Bus 2)", reactive = -0.011)
updatePmu!(system, device; label = "PMU 1 (Bus 2)", polar = false)
analysis = gaussNewton(system, device) # <- Creating ACStateEstimation for the updated model
for iteration = 1:20
stopping = solve!(system, analysis)
if stopping < 1e-8
break
end
endThis method 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 ACStateEstimation composite type using gaussNewton or acLavStateEstimation just once. After this initial setup, users can seamlessly modify existing measurement devices without the need to recreate the ACStateEstimation type.
This advancement extends beyond the previous scenario where recreating the Measurement type was unnecessary, to now include the scenario where ACStateEstimation also does not need to be recreated.
The addition of new measurements after the creation of ACStateEstimation is not practical in terms of reusing this 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.
Weighted Least-squares Estimator
We can modify the prior example to achieve the same model without establishing ACStateEstimation twice:
system = powerSystem()
device = measurement() # <- Initializing a Measurement instance
addBus!(system; label = "Bus 1", type = 3)
addBus!(system; label = "Bus 2", type = 1, active = 0.1, reactive = 0.01)
addBus!(system; label = "Bus 3", type = 1, active = 2.5, reactive = 0.2)
@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)
addGenerator!(system; label = "Generator 1", bus = "Bus 1", active = 3.2, reactive = 0.3)
@voltmeter(label = "Voltmeter ? (!)", variance = 1e-3)
addVoltmeter!(system, device; bus = "Bus 1", magnitude = 1.0)
@wattmeter(label = "Wattmeter ? (!)", varianceBus = 1e-2)
addWattmeter!(system, device; from = "Branch 1", active = 1.046)
addWattmeter!(system, device; bus = "Bus 2", active = -0.1)
addWattmeter!(system, device; from = "Branch 3", active = 0.924, status = 0)
@varmeter(label = "Varmeter ? (!)", varianceFrom = 1e-3)
addVarmeter!(system, device; from = "Branch 1", reactive = 0.059)
addVarmeter!(system, device; bus = "Bus 2", reactive = -0.01)
addVarmeter!(system, device; to = "Branch 3", reactive = -0.044, variance = 1e-5, status = 0)
@pmu(label = "PMU ? (!)")
addPmu!(system, device; bus = "Bus 2", magnitude = 0.976, angle = -0.052)
analysis = gaussNewton(system, device) # <- Creating ACStateEstimation for the defined model
for iteration = 1:20
stopping = solve!(system, analysis)
if stopping < 1e-8
break
end
end
updateWattmeter!(system, device, analysis; label = "Wattmeter 3 (From Branch 3)", status = 1)
updateWattmeter!(system, device, analysis; label = "Wattmeter 2 (Bus 2)", variance = 1e-4)
updateVarmeter!(system, device, analysis; label = "Varmeter 3 (To Branch 3)", status = 1)
updateVarmeter!(system, device, analysis; label = "Varmeter 2 (Bus 2)", reactive = -0.011)
updatePmu!(system, device, analysis; label = "PMU 1 (Bus 2)", polar = false)
# <- No need for re-creation; we have already updated the existing ACStateEstimation instance
for iteration = 1:20
stopping = solve!(system, analysis)
if stopping < 1e-8
break
end
endThis method removes the need to restart and recreate both the Measurement and the ACStateEstimation from the beginning when implementing changes to the existing measurement set. Subsequently, JuliaGrid can leverage the symbolic factorizations of LU or LDLt, provided that the nonzero pattern of the gain matrix remains unchanged. This method avoids the need to initialize the AC state estimation model from scratch.
Least Absolute Value Estimator
The same methodology can be applied to the LAV method, thereby circumventing the need to construct an optimization model from scratch.
Power and Current Analysis
After obtaining the solution from the AC 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 AC state estimation solution:
system = powerSystem()
device = measurement()
addBus!(system; label = "Bus 1", type = 3, susceptance = 0.002)
addBus!(system; label = "Bus 2", type = 1, active = 0.1, reactive = 0.01)
addBus!(system; label = "Bus 3", type = 1, active = 2.5, reactive = 0.2)
@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)
addGenerator!(system; label = "Generator 1", bus = "Bus 1", active = 3.2, reactive = 0.3)
addWattmeter!(system, device; from = "Branch 1", active = 1.046, variance = 1e-2)
addWattmeter!(system, device; bus = "Bus 2", active = -0.1, variance = 1e-3)
addWattmeter!(system, device; from = "Branch 3", active = 0.924, variance = 1e-3)
addVarmeter!(system, device; from = "Branch 1", reactive = 0.059, variance = 1e-3)
addVarmeter!(system, device; bus = "Bus 2", reactive = -0.01, variance = 1e-2)
addVarmeter!(system, device; to = "Branch 3", reactive = -0.044, variance = 1e-3)
analysis = gaussNewton(system, device)
for iteration = 1:20
stopping = solve!(system, analysis)
if stopping < 1e-8
break
end
endWe can now utilize the provided functions to compute powers and currents:
power!(system, analysis)
current!(system, analysis)For instance, if we want to show the active power injections and the from-bus current angles, we can employ the following code:
julia> print(system.bus.label, analysis.power.injection.active)Bus 1: 2.6951392445386286 Bus 2: -0.09999730749216137 Bus 3: -2.499287793840303julia> print(system.branch.label, analysis.current.from.angle)Branch 1: -0.05634389548038201 Branch 2: -0.15413096811690216 Branch 3: -0.08788060713203448
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 AC State Estimation.
To compute specific quantities for particular components, rather than calculating powers or currents for all components, users can utilize one of the provided functions below.
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(system, analysis; label = "Bus 1")(2.6951392445386286, 0.31320449579404946)
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(system, analysis; label = "Bus 1")(2.6951392445386286, 0.31320449579404946)
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(system, analysis; label = "Bus 1")(0.0, -0.002006637952521858)
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(system, analysis; label = "Branch 2")(1.6491130561317453, 0.2562114962034615)julia> active, reactive = toPower(system, analysis; label = "Branch 2")(-1.5932872079562277, -0.15530958620622323)
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(system, analysis; label = "Branch 1")(9.800363372293309e-5, -0.039201453489173234)
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(system, analysis; label = "Branch 2")(0.055729791752061136, 0.13932447938015285)
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(system, analysis; label = "Bus 1")(2.708785622827295, -0.11569193435017591)
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(system, analysis; label = "Branch 2")(1.6661346611128796, -0.15413096811690216)julia> magnitude, angle = toCurrent(system, analysis; label = "Branch 2")(1.6709804057016453, 2.964170792091241)
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(system, analysis; label = "Branch 2")(1.6692781636393188, -0.16599467621735425)
References
[1] G. Korres, Observability analysis based on echelon form of a reduced dimensional Jacobian matrix, IEEE Trans. Power Syst., vol. 26, no. 4, pp. 2572–2573, 2011.
[2] G. N. Korres and N. M. Manousakis, State estimation and observability analysis for phasor measurement unit measured systems, IET Gener. Transm. Dis., vol. 6, no. 9, 2012.
[3] A. Gomez-Exposito, A. Abur, P. Rousseaux, A. de la Villa Jaen, and C. Gomez-Quiles, On the use of PMUs in power system state estimation, Proc. IEEE PSCC, 2011.
[4] M. Zhou, V. A. Centeno, J. S. Thorp, and A. G. Phadke, An alternative for including phasor measurements in state estimators, IEEE Trans. Power Syst., vol. 21, no. 4, 2006.
[5] A. Abur and A. Exposito, Power System State Estimation: Theory and Implementation, Taylor & Francis, 2004.