PMU State Estimation
To perform linear state estimation solely based on PMU data, the initial requirement is to have the PowerSystem composite type configured with the AC model, along with the Measurement composite type storing measurement data. Subsequently, we can formulate either the weighted least-squares (WLS) or the least absolute value (LAV) PMU state estimation model encapsulated within the type PMUStateEstimation using:
For resolving the PMU state estimation problem and obtaining bus voltage magnitudes and angles, utilize the following function:
After executing the function solve!, where the user employs the WLS 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 PMUStateEstimation type, users can initiate an optimal PMU placement algorithm to determine the minimal set of PMUs required for an observable system:
After obtaining the PMU 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:
Optimal PMU Placement
Let us define the PowerSystem composite type and perform the AC power flow analysis solely for generating data to artificially create measurement values:
system = powerSystem()
addBus!(system; label = "Bus 1", type = 3, active = 0.5)
addBus!(system; label = "Bus 2", type = 1, reactive = 0.05)
addBus!(system; label = "Bus 3", type = 1, active = 0.5)
@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 2", reactance = 0.01)
addBranch!(system; label = "Branch 3", from = "Bus 2", to = "Bus 3", reactance = 0.04)
@generator(reactive = 0.1)
addGenerator!(system; label = "Generator 1", bus = "Bus 1", active = 3.2)
addGenerator!(system; label = "Generator 2", bus = "Bus 2", active = 2.1)
analysis = newtonRaphson(system)
for iteration = 1:10
mismatch!(system, analysis)
solve!(system, analysis)
endOptimal Solution
Upon defining the PowerSystem composite type, JuliaGrid provides the possibility to determine the minimal number of PMUs required for system observability using the pmuPlacement function:
using GLPK
placement = pmuPlacement(system, GLPK.Optimizer)The placement variable contains data regarding the optimal placement of measurements. In this instance, installing a PMU at Bus 2 renders the system observable:
julia> placement.busOrderedCollections.OrderedDict{String, Int64} with 1 entry: "Bus 2" => 2
This PMU installed at Bus 2 will measure the bus voltage phasor at the corresponding bus and all current phasors at the branches incident to Bus 2 located at the from-bus or to-bus ends:
julia> placement.fromOrderedCollections.OrderedDict{String, Int64} with 1 entry: "Branch 3" => 3julia> placement.toOrderedCollections.OrderedDict{String, Int64} with 2 entries: "Branch 1" => 1 "Branch 2" => 2
Measurement Data
Utilizing PMU placement and AC power flow data, which serves as the source for measurement values in this scenario, we can construct the Measurement composite type as follows:
device = measurement()
@pmu(label = "PMU ? (!)")
for (bus, k) in placement.bus
Vᵢ, θᵢ = analysis.voltage.magnitude[k], analysis.voltage.angle[k]
addPmu!(system, device; bus = bus, magnitude = Vᵢ, angle = θᵢ)
end
for branch in keys(placement.from)
Iᵢⱼ, ψᵢⱼ = fromCurrent(system, analysis; label = branch)
addPmu!(system, device; from = branch, magnitude = Iᵢⱼ, angle = ψᵢⱼ)
end
for branch in keys(placement.to)
Iⱼᵢ, ψⱼᵢ = toCurrent(system, analysis; label = branch)
addPmu!(system, device; to = branch, magnitude = Iⱼᵢ, angle = ψⱼᵢ)
endFor example, we can observe the obtained set of measurement values:
julia> print(device.pmu.label, device.pmu.magnitude.mean, device.pmu.angle.mean)PMU 1 (Bus 2): 1.021562852206714, 0.014634612107908078 PMU 2 (From Branch 3): 0.4952989219070908, 0.07710792487077446 PMU 3 (To Branch 1): 0.47443528099140664, -0.5459877069301735 PMU 4 (To Branch 2): 1.1722759871967534, 0.16223095483453315
Weighted Least-squares Estimator
Let us continue with the previous example, where we defined the PowerSystem and Measurement types. To establish the PMU state estimation model, we will use the pmuWlsStateEstimation function:
analysis = pmuWlsStateEstimation(system, device)Here, the user triggers LU factorization as the default method for solving the PMU state estimation problem. However, the user also has the option to select alternative factorization methods such as LDLt or QR:
analysis = pmuWlsStateEstimation(system, device, QR)To obtain the bus voltage magnitudes and angles, the solve! function can be invoked as shown:
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.9999999999998835, -9.992662462743827e-14 Bus 2: 1.021562852206596, 0.01463461210781196 Bus 3: 1.0122941004122825, -0.005105099383630845
We recommend that readers refer to the tutorial on PMU State Estimation for insights into the implementation.
Correlated Measurement Errors
In the above approach, we assume that measurement errors from a single PMU are uncorrelated. This assumption leads to the covariance matrix and its inverse matrix (i.e., precision matrix) maintaining a diagonal form:
julia> analysis.method.precision8×8 SparseArrays.SparseMatrixCSC{Float64, Int64} with 8 stored entries: 99999.1 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 95823.9 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 1.0045e5 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 4.00322e5 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 126413.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 2.30398e5 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 99033.1 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 73288.7
While this approach is suitable for many scenarios, linear PMU state estimation relies on transforming from polar to rectangular coordinate systems. Consequently, measurement errors from a single PMU become correlated due to this transformation. This correlation results in the covariance matrix, and hence the precision matrix, no longer maintaining a diagonal form but instead becoming a block diagonal matrix.
To accommodate this, users have the option to consider correlation when adding each PMU to the Measurement type. For instance, let us add a new PMU while considering correlation:
addPmu!(system, device; bus = "Bus 3", magnitude = 1.01, angle = -0.005, correlated = true)Following this, we recreate the WLS state estimation model:
analysis = pmuWlsStateEstimation(system, device)Upon inspection, it becomes evident that the precision matrix no longer maintains a diagonal structure:
julia> analysis.method.precision10×10 SparseArrays.SparseMatrixCSC{Float64, Int64} with 12 stored entries: 99999.1 ⋅ ⋅ ⋅ … ⋅ ⋅ ⋅ ⋅ 95823.9 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 1.0045e5 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 4.00322e5 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ … ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 73288.7 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 100000.0 -9.85181 ⋅ ⋅ ⋅ ⋅ ⋅ -9.85181 98029.7
Subsequently, we can address this new scenario and observe the solution:
julia> solve!(system, analysis)julia> print(system.bus.label, analysis.voltage.magnitude, analysis.voltage.angle)Bus 1: 0.9988543535981543, 5.992024114227854e-5 Bus 2: 1.0204177444279554, 0.014710234392911753 Bus 3: 1.0111477069326744, -0.00505223619159087
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 [1, 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 5 (Bus 3)", correlated = false)Subsequently, by specifying the Orthogonal argument in the pmuWlsStateEstimation function, JuliaGrid implements a more robust approach to obtain the WLS estimator, which proves particularly beneficial when substantial differences exist among measurement variances:
analysis = pmuWlsStateEstimation(system, device, Orthogonal)
solve!(system, analysis)Bad Data Processing
After acquiring the WLS solution using the solve! function, 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 phasor measurement. Upon proceeding to find the solution for this updated state:
addPmu!(system, device; bus = "Bus 3", magnitude = 3.2, angle = 0.0, noise = false)
analysis = pmuWlsStateEstimation(system, device)
solve!(system, analysis)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 PMU 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.maxNormalizedResidual565.215511905541julia> outlier.label"PMU 6 (Bus 3)"
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 phasor measurement as out-of-service within the Measurement type.
Moreover, JuliaGrid will adjust the coefficient matrix and mean vector within the PMUStateEstimation type based on measurements now designated as out-of-service. To optimize the algorithm's efficiency, JuliaGrid resets non-zero elements to zero in the coefficient matrix and mean vector, effectively removing the impact of the corresponding measurement on the solution:
julia> analysis.method.mean12-element Vector{Float64}: 1.0214534591451636 0.014949642440409396 0.4938272187286028 0.038153637901242377 0.40545942634957655 -0.24635642762807566 1.1568832944388845 0.18934633137812892 1.009987375026302 -0.005049978958359636 0.0 0.0julia> analysis.method.coefficient12×6 SparseArrays.SparseMatrixCSC{Float64, Int64} with 30 stored entries: ⋅ 1.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 1.0 ⋅ ⋅ 10.0 -10.0 ⋅ 19.98 -20.0 ⋅ -19.98 20.0 ⋅ 10.0 -10.0 -6.89655 6.8966 ⋅ -17.2414 17.2214 ⋅ 17.2414 -17.2214 ⋅ -6.89655 6.8966 ⋅ -40.0 40.0001 ⋅ -20.0 19.98 ⋅ 20.0 -19.98 ⋅ -40.0 40.0001 ⋅ ⋅ ⋅ 1.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 1.0 ⋅ ⋅ 0.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.0
After removing bad data, a new estimate can be computed without considering this specific phasor measurement:
solve!(system, analysis)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 [1, 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 = pmuLavStateEstimation(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 PMUStateEstimation 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. It is important to note that JuliaGrid utilizes the provided data to set starting primal values in the rectangular coordinate system.
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.9999987227348247, -4.0017727041788876e-7 Bus 2: 1.0215625907440575, 0.014634738933430869 Bus 3: 1.0122928555607549, -0.0051070145522880755
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 PMUStateEstimation type using the pmuWlsStateEstimation or pmuLavStateEstimation function. Ultimately, resolving the PMU 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 2", to = "Bus 3", reactance = 0.03)
addGenerator!(system; label = "Generator 1", bus = "Bus 1", active = 3.2, reactive = 0.3)
@pmu(label = "PMU ?")
addPmu!(system, device; bus = "Bus 1", magnitude = 1.0, angle = 0.0)
addPmu!(system, device; bus = "Bus 2", magnitude = 0.98, angle = -0.023)
addPmu!(system, device; from = "Branch 2", magnitude = 0.5, angle = -0.05)
analysis = pmuWlsStateEstimation(system, device) # <- Creating PMUStateEstimation for model
solve!(system, analysis)
addPmu!(system, device; to = "Branch 2", magnitude = 0.5, angle = 3.1)
updatePmu!(system, device; label = "PMU 1", varianceMagnitude = 1e-8)
updatePmu!(system, device; label = "PMU 3", statusMagnitude = 0, statusAngle = 0)
analysis = pmuWlsStateEstimation(system, device) #<- Creating PMUStateEstimation for new model
solve!(system, analysis)This 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 PMUStateEstimation composite type using pmuWlsStateEstimation or pmuLavStateEstimation just once. After this initial setup, users can seamlessly modify existing measurement devices without the need to recreate the PMUStateEstimation type.
This advancement extends beyond the previous scenario where recreating the Measurement type was unnecessary, to now include the scenario where PMUStateEstimation also does not need to be recreated.
The addition of new measurements after the creation of PMUStateEstimation is not practical in terms of reusing the PMUStateEstimation 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 PMUStateEstimation 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 2", to = "Bus 3", reactance = 0.03)
addGenerator!(system; label = "Generator 1", bus = "Bus 1", active = 3.2, reactive = 0.3)
@pmu(label = "PMU ?")
addPmu!(system, device; bus = "Bus 1", magnitude = 1.0, angle = 0.0)
addPmu!(system, device; bus = "Bus 2", magnitude = 0.98, angle = -0.023)
addPmu!(system, device; from = "Branch 2", magnitude = 0.5, angle = -0.05)
addPmu!(system, device; to = "Branch 2", magnitude = 0.5, angle = 3.1, statusAngle = 0)
analysis = pmuWlsStateEstimation(system, device) # <- Creating PMUStateEstimation for model
solve!(system, analysis)
updatePmu!(system, device, analysis; label = "PMU 1", varianceMagnitude = 1e-8)
updatePmu!(system, device, analysis; label = "PMU 3", statusMagnitude = 0, statusAngle = 0)
updatePmu!(system, device, analysis; label = "PMU 4", statusAngle = 1)
# <- No need for re-creation; we have already updated the existing PMUStateEstimation instance
solve!(system, analysis)This method removes the need to restart and recreate both the Measurement and the PMUStateEstimation from the beginning when implementing changes to the existing measurement set. Next, JuliaGrid can reuse symbolic factorizations of LU or LDLt, as long as the nonzero pattern of the gain matrix remains consistent.
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 PMU 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 PMU 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)
addPmu!(system, device; bus = "Bus 1", magnitude = 1.0, angle = 0.0)
addPmu!(system, device; bus = "Bus 2", magnitude = 0.97, angle = -0.051)
addPmu!(system, device; from = "Branch 2", magnitude = 1.66, angle = -0.15)
addPmu!(system, device; to = "Branch 2", magnitude = 1.67, angle = 2.96)
analysis = pmuWlsStateEstimation(system, device)
solve!(system, analysis)We 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 magnitudes, we can employ the following code:
julia> print(system.bus.label, analysis.power.injection.active)Bus 1: 2.7128789695735005 Bus 2: -0.21564110130496367 Bus 3: -2.402882441531764julia> print(system.branch.label, analysis.current.from.magnitude)Branch 1: 1.084804224496155 Branch 2: 1.6624663614086073 Branch 3: 0.8658336251368861
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 PMU 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.712878969573501, 0.43288460463162526)
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.712878969573501, 0.43288460463162526)
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.0020004242557670198)
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.6429108157161207, 0.2553913316369414)julia> active, reactive = toPower(system, analysis; label = "Branch 2")(-1.587330137890149, -0.1549833506017014)
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.705560639417549e-5, -0.03882224255767019)
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.05548491704701527, 0.13871229261753817)
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.746907492685579, -0.1581787316058082)
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.6624663614086073, -0.15416239142972624)julia> magnitude, angle = toCurrent(system, analysis; label = "Branch 2")(1.667309344330871, 2.964124299788863)
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.6656067520128401, -0.16603367260017668)
References
[1] A. Abur and A. Exposito, Power System State Estimation: Theory and Implementation, Taylor & Francis, 2004.