DC State Estimation
To perform the DC state estimation, we first need to have the PowerSystem composite type that has been created with the DC model, alongside the Measurement composite type that retains measurement data. Subsequently, we can formulate either the weighted least-squares (WLS) or the least absolute value (LAV) DC state estimation model encapsulated within the abstract type DCStateEstimation using:
For resolving the DC state estimation problem and obtaining bus voltage 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 DCStateEstimation type, users can initiate observability analysis to identify observable islands and restore observability by employing:
For more detailed information, users can refer to the Observability Analysis section within AC state estimation documentation. It is worth noting that when the system becomes observable within the AC model, it will also be observable within the DC state estimation model.
After obtaining the solution for DC state estimation, JuliaGrid offers a post-processing analysis function to compute active powers associated with buses and branches:
Additionally, there are specialized functions dedicated to calculating specific types of active powers related to particular buses or branches:
Bus Type Modification
Just like in the Bus Type Modification section, when establishing the DCStateEstimation 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).
Weighted Least-squares Estimator
To solve the DC state estimation and derive WLS estimates using JuliaGrid, the process initiates by defining the composite types PowerSystem and Measurement. Here is an illustrative example:
system = powerSystem()
device = measurement()
addBus!(system; label = "Bus 1", type = 3)
addBus!(system; label = "Bus 2", type = 1, active = 0.2)
addBus!(system; label = "Bus 3", type = 1, active = 0.4)
addBranch!(system; label = "Branch 1", from = "Bus 1", to = "Bus 2", reactance = 0.5)
addBranch!(system; label = "Branch 2", from = "Bus 1", to = "Bus 3", reactance = 0.2)
addBranch!(system; label = "Branch 3", from = "Bus 2", to = "Bus 3", reactance = 0.3)
addGenerator!(system; label = "Generator 1", bus = "Bus 1", active = 3.2)
@wattmeter(label = "Wattmeter ?")
addWattmeter!(system, device; bus = "Bus 1", active = 0.6, variance = 1e-3)
addWattmeter!(system, device; bus = "Bus 3", active = -0.4, variance = 1e-2)
addWattmeter!(system, device; from = "Branch 1", active = 0.18, variance = 1e-4)
addWattmeter!(system, device; to = "Branch 2", active = -0.42, variance = 1e-4)The dcWlsStateEstimation function serves to establish the DC state estimation problem:
analysis = dcWlsStateEstimation(system, device)Here, the user triggers LU factorization as the default method for solving the DC state estimation problem. However, the user also has the option to select alternative factorization methods such as LDLt or QR:
analysis = dcWlsStateEstimation(system, device, LDLt)To obtain the bus voltage angles, the solve! function can be invoked as shown:
solve!(system, analysis)Upon obtaining the solution, access the bus voltage angles using:
julia> print(system.bus.label, analysis.voltage.angle)Bus 1: 0.0 Bus 2: -0.09000000000000001 Bus 3: -0.08399999999999999
We recommend that readers refer to the tutorial on DC State Estimation for insights into the implementation.
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 [2, Sec. 3.2].
Specifically, by specifying the Orthogonal argument in the dcWlsStateEstimation function, JuliaGrid implements a more robust approach to obtain the WLS estimator, which proves particularly beneficial when substantial differences exist among measurement variances:
analysis = dcWlsStateEstimation(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 wattmeter. Upon proceeding to find the solution for this updated state:
addWattmeter!(system, device; from = "Branch 2", active = 4.1, variance = 1e-4)
analysis = dcWlsStateEstimation(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 DC 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.maxNormalizedResidual267.65576812169894julia> outlier.label"Wattmeter 5"
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 measurements as out-of-service within the Measurement type.
Moreover, JuliaGrid will adjust the coefficient matrix and mean vector within the DCStateEstimation 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.mean5-element Vector{Float64}: 0.6 -0.4 0.18 -0.42 0.0julia> analysis.method.coefficient5×3 SparseArrays.SparseMatrixCSC{Float64, Int64} with 12 stored entries: 7.0 -2.0 -5.0 -5.0 -3.33333 8.33333 2.0 -2.0 ⋅ -5.0 ⋅ 5.0 0.0 ⋅ 0.0
Hence, after removing bad data, a new estimate can be computed without considering this specific 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 [2, 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 = dcLavStateEstimation(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 DCStateEstimation type. By default, these values are initially established using the the initial bus voltage angles from PowerSystem type:
julia> print(system.bus.label, analysis.voltage.angle)Bus 1: 0.0 Bus 2: 0.0 Bus 3: 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 angles using:
julia> print(system.bus.label, analysis.voltage.angle)Bus 1: 0.0 Bus 2: -0.08999999999956482 Bus 3: -0.08399999999987212
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 DCStateEstimation type using the dcWlsStateEstimation or dcLavStateEstimation function. Ultimately, resolving the DC 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)
addBranch!(system; label = "Branch 1", from = "Bus 1", to = "Bus 2", reactance = 0.5)
addGenerator!(system; label = "Generator 1", bus = "Bus 1", active = 0.1)
@wattmeter(label = "Wattmeter ?")
addWattmeter!(system, device; bus = "Bus 2", active = -0.11, variance = 1e-3)
addWattmeter!(system, device; from = "Branch 1", active = 0.09, variance = 1e-4)
analysis = dcWlsStateEstimation(system, device) # <- Creating DCStateEstimation for the model
solve!(system, analysis)
addWattmeter!(system, device; to = "Branch 1", active = -0.12, variance = 1e-4)
updateWattmeter!(system, device; label = "Wattmeter 1", status = 0)
updateWattmeter!(system, device; label = "Wattmeter 2", active = 0.1, noise = false)
analysis = dcWlsStateEstimation(system, device) # <- Creating DCStateEstimation 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 DCStateEstimation composite type using dcWlsStateEstimation or dcLavStateEstimation just once. After this initial setup, users can seamlessly modify existing measurement devices without the need to recreate the DCStateEstimation type.
This advancement extends beyond the previous scenario where recreating the Measurement type was unnecessary, to now include the scenario where DCStateEstimation also does not need to be recreated. Such efficiency can be particularly advantageous in cases where JuliaGrid can reuse gain matrix factorization.
The addition of new measurements after the creation of DCStateEstimation 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 DCStateEstimation 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)
addBranch!(system; label = "Branch 1", from = "Bus 1", to = "Bus 2", reactance = 0.5)
addGenerator!(system; label = "Generator 1", bus = "Bus 1", active = 0.1)
@wattmeter(label = "Wattmeter ?")
addWattmeter!(system, device; bus = "Bus 2", active = -0.11, variance = 1e-3)
addWattmeter!(system, device; from = "Branch 1", active = 0.09, variance = 1e-4)
addWattmeter!(system, device; to = "Branch 1", active = -0.12, variance = 1e-4, status = 0)
analysis = dcWlsStateEstimation(system, device) # <- Creating DCStateEstimation for the model
solve!(system, analysis)
updateWattmeter!(system, device; label = "Wattmeter 1", status = 0)
updateWattmeter!(system, device; label = "Wattmeter 2", active = 0.1, noise = false)
updateWattmeter!(system, device; label = "Wattmeter 3", status = 1)
# <- No need for re-creation; we have already updated the existing DCStateEstimation instance
solve!(system, analysis)This method removes the need to restart and recreate both the Measurement and the DCStateEstimation 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.
Reusing Weighted Least-squares Matrix Factorization
Drawing from the preceding example, our focus now shifts to finding a solution involving modifications that entail adjusting the measurement value of the Wattmeter 2. It is important to note that these adjustments do not impact the variance or status of the measurement device, which can affect the gain matrix. To resolve this updated system, users can simply execute the solve! function:
updateWattmeter!(system, device, analysis; label = "Wattmeter 2", active = 0.091)
solve!(system, analysis)In this scenario, JuliaGrid will recognize instances where the user has not modified parameters that impact the gain matrix. Consequently, JuliaGrid will leverage the previously performed gain matrix factorization, resulting in a significantly faster solution compared to recomputing the factorization.
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 Analysis
After obtaining the solution from the DC state estimation, calculating powers related to buses and branches is facilitated by using the power! function. For instance, let us consider the model for which we obtained the DC state estimation solution:
system = powerSystem()
device = measurement()
addBus!(system; label = "Bus 1", type = 3, conductance = 1e-3)
addBus!(system; label = "Bus 2", type = 1, active = 0.2)
addBus!(system; label = "Bus 3", type = 1, active = 0.4)
addBranch!(system; label = "Branch 1", from = "Bus 1", to = "Bus 2", reactance = 0.5)
addBranch!(system; label = "Branch 2", from = "Bus 1", to = "Bus 3", reactance = 0.2)
addBranch!(system; label = "Branch 3", from = "Bus 2", to = "Bus 3", reactance = 0.3)
addGenerator!(system; label = "Generator 1", bus = "Bus 1", active = 3.2)
addWattmeter!(system, device; bus = "Bus 1", active = 0.6, variance = 1e-3)
addWattmeter!(system, device; bus = "Bus 3", active = -0.4, variance = 1e-2)
addWattmeter!(system, device; from = "Branch 1", active = 0.18, variance = 1e-4)
addWattmeter!(system, device; to = "Branch 2", active = -0.42, variance = 1e-4)
analysis = dcWlsStateEstimation(system, device)
solve!(system, analysis)We can compute active powers using the following function:
power!(system, analysis)For example, active power injections corresponding to buses are:
julia> print(system.bus.label, analysis.power.injection.active)Bus 1: 0.6008333333333332 Bus 2: -0.19983333333333342 Bus 3: -0.4
To better understand the powers associated with buses, and branches that are calculated by the power! function, we suggest referring to the tutorials on.
To calculate specific quantities related to particular buses or branches, rather than computing values for all buses and branches, users can utilize one of the provided functions below.
Active Power Injection
To calculate active power injection associated with a specific bus, the function can be used:
julia> active = injectionPower(system, analysis; label = "Bus 1")0.6008333333333333
Active Power Injection from Generators
To calculate active power injection from the generators at a specific bus, the function can be used:
julia> active = supplyPower(system, analysis; label = "Bus 1")0.6008333333333333
Active Power Flow
Similarly, we can compute the active power flow at both the from-bus and to-bus ends of the specific branch by utilizing the provided functions below:
julia> active = fromPower(system, analysis; label = "Branch 1")0.17991666666666667julia> active = toPower(system, analysis; label = "Branch 1")-0.17991666666666667
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] A. Abur and A. Exposito, Power System State Estimation: Theory and Implementation, Taylor & Francis, 2004.