State Estimation

islandTopologicalFlow(monitoring::Measurement)

The function utilizes a topological approach to detect flow-observable islands, resulting in the formation of disconnected and loop-free subgraphs. It is assumed that active and reactive power measurements are paired, indicating a standard observability analysis. In this analysis, islands formed by active power measurements correspond to those formed by reactive power measurements.

Arguments

To define flow-observable islands, this function needs the Measurement type.

Returns

The function returns an instance of the Island type.

Example

system, monitoring = ems("case14.h5", "monitoring.h5")

statusWattmeter!(monitoring; inservice = 15)
monitoring.varmeter.reactive.status = copy(monitoring.wattmeter.active.status)

islands = islandTopologicalFlow(monitoring)

islandTopological(monitoring::Measurement)

The function employs a topological method to identify maximal-observable islands. Specifically, it employs active power measurements to pinpoint flow-observable islands. Subsequently, these islands are merged based on the available injection measurements.

It is assumed that active and reactive power measurements are paired, indicating a standard observability analysis. In this analysis, islands formed by active power measurements correspond to those formed by reactive power measurements.

Arguments

To define maximal-observable islands, this function needs the Measurement type.

Returns

The function returns an instance of the Island type.

Example

system, monitoring = ems("case14.h5", "monitoring.h5")

statusWattmeter!(monitoring; inservice = 15)
monitoring.varmeter.reactive.status = copy(monitoring.wattmeter.active.status)

islands = islandTopological(monitoring)

restorationGram!(monitoring::Measurement, pseudo::Measurement, islands::Island; threshold)

Upon identifying the Island, the function incorporates measurements from the available pseudo-measurements in the pseudo variable into the monitoring variable to reinstate observability. This method relies on reduced coefficient matrices and the Gram matrix.

It is important to note that the monitoring labels in the monitoring and pseudo variables must be different to enable the function to successfully incorporate measurements from pseudo into the monitoring set of measurements.

Keyword

The keyword threshold defines the zero pivot threshold value, with a default value of 1e-5. More precisely, all computed pivots less than this value will be treated as zero pivots.

Updates

The function updates the monitoring variable of the Measurement type.

Example

system, monitoring, pseudo = ems("case14.h5", "monitoring.h5", "pseudo.h5")

statusWattmeter!(monitoring; inservice = 10)
islands = islandTopological(monitoring)

restorationGram!(monitoring, pseudomonitoring, islands)

pmuPlacement(system::PowerSystem, optimizer; bridge, name, verbose)

The function determines the optimal placement of PMUs through integer linear programming. It identifies the minimum set of PMUs required to ensure observability and a unique state estimator.

The function accepts a PowerSystem type as input to establish the framework for finding the optimal PMU placement. If the ac field within the PowerSystem type is not yet created, the function automatically initiates an update process.

Additionally, the optimizer argument is a crucial component for formulating and solving the optimization problem. Typically, using the HiGHS or GLPK solver is sufficient. For more detailed information, please refer to the JuMP documenatation.

Keywords

The function accepts the following keywords:

• bridge: Controls the bridging mechanism (default: false).
• name: Handles the creation of string names (default: false).
• verbose: Controls the output display, ranging from silent mode (0) to detailed output (3).

Returns

The function returns an instance of the PmuPlacement type.

Note that if a PMU is understood as a device that measures the bus voltage phasor and all branch current phasors incident to the bus, users only need the results stored in the bus variable. However, if a PMU is considered to measure individual phasor, then all required phasor measurements can be found in the bus, from, and to variables.

Example

using HiGHS, Ipopt

system, monitoring = ems("case14.h5")

analysis = acOptimalPowerFlow(system, Ipopt.Optimizer)
powerFlow!(analysis; current = true)

placement = pmuPlacement(system, HiGHS.Optimizer)

@pmu(label = "PMU ?: !")
for (bus, i) in placement.bus
    Vi, θi = analysis.voltage.magnitude[i], analysis.voltage.angle[i]
    addPmu!(monitoring; bus = bus, magnitude = Vi, angle = θi)
end
for branch in keys(placement.from)
    Iij, ψij = fromCurrent(analysis; label = branch)
    addPmu!(monitoring; from = branch, magnitude = Iij, angle = ψij)
end
for branch in keys(placement.to)
    Iji, ψji = toCurrent(analysis; label = branch)
    addPmu!(monitoring; to = branch, magnitude = Iji, angle = ψji)
end

pmuPlacement!(monitoring::Measurement, analysis::AC, optimizer;
    varianceMagnitudeBus, varianceAngleBus,
    varianceMagnitudeFrom, varianceAngleFrom,
    varianceMagnitudeTo, varianceAngleTo,
    noise, correlated, polar,
    bridge, name, verbose)

The function finds the optimal PMU placement by executing pmuPlacement function. Then, based on the results from the AC type, it generates phasor measurements and integrates them into the Measurement type. If current values are missing in the AC type, the function calculates the associated currents required to form measurement values.

Keywords

PMUs at the buses can be configured using:

• varianceMagnitudeBus (pu or V): Variance of bus voltage magnitude measurements.
• varianceAngleBus (rad or deg): Variance of bus voltage angle measurements.

PMUs at the from-bus ends of the branches can be configured using:

• varianceMagnitudeFrom (pu or A): Variance of current magnitude measurements.
• varianceAngleFrom (rad or deg): Variance of current angle measurements.

PMUs at the to-bus ends of the branches can be configured using:

• varianceMagnitudeTo (pu or A): Variance of current magnitude measurements.
• varianceAngleTo (rad or deg): Variance of current angle measurements.

Settings for generating measurements include:

• noise: Defines the method for generating the measurement means:
  • noise = true: adds white Gaussian noise to the phasor values, using the defined variances,
  • noise = false: uses the exact phasor values without adding noise.

Settings for handling phasor measurements include:

• correlated: Specifies error correlation for PMUs for algorithms utilizing rectangular coordinates:
  • correlated = true: considers correlated errors,
  • correlated = false: disregards correlations between errors.
• polar: Chooses the coordinate system for including phasor measurements in AC state estimation:
  • polar = true: adopts the polar coordinate system,
  • polar = false: adopts the rectangular coordinate system.

Settings for the optimization solver include:

• bridge: Controls the bridging mechanism (default: false).
• name: Handles the creation of string names (default: false).
• verbose: Controls the output display, ranging from silent mode (0) to detailed output (3).

Updates

The function updates the pmu field of the Measurement type.

Returns

The function returns an instance of the PmuPlacement type.

Example

using HiGHS, Ipopt

system, monitoring = ems("case14.h5")

analysis = acOptimalPowerFlow(system, Ipopt.Optimizer)
powerFlow!(analysis; current = true)

pmuPlacement!(monitoring, analysis, HiGHS.Optimizer)

gaussNewton(monitoring::Measurement, [method = LU])

The function sets up the Gauss-Newton method to solve the nonlinear or AC state estimation model, where the vector of state variables is given in polar coordinates. The Gauss-Newton method throughout iterations provided WLS estimator.

Arguments

This function requires the Measurement type to establish the WLS state estimation framework.

Moreover, the presence of the method parameter is not mandatory. To address the WLS state estimation method, users can opt to utilize factorization techniques to decompose the gain matrix, such as LU, QR, or LDLt especially when the gain matrix is symmetric. Opting for the Orthogonal method is advisable for a more robust solution in scenarios involving ill-conditioned data, particularly when substantial variations in variances are present.

If the user does not provide the method, the default method for solving the estimation model will be LU factorization.

Updates

If the AC model has not been created, the function will automatically trigger an update of the ac field within the PowerSystem type.

Returns

The function returns an instance of the AcStateEstimation type.

Examples

Set up the AC state estimation model to be solved using the default LU factorization:

system, monitoring = ems("case14.h5", "monitoring.h5")

analysis = gaussNewton(monitoring)

Set up the AC state estimation model to be solved using the orthogonal method:

system, monitoring = ems("case14.h5", "monitoring.h5")

analysis = gaussNewton(monitoring, Orthogonal)

acLavStateEstimation(monitoring::Measurement, optimizer;
    iteration, tolerance, bridge, name, magnitude, angle, positive, negative, verbose)

The function sets up the LAV method to solve the nonlinear or AC state estimation model, where the vector of state variables is given in polar coordinates.

Arguments

This function requires the Measurement type to establish the LAV state estimation model. The LAV method offers increased robustness compared to WLS, ensuring unbiasedness even in the presence of various measurement errors and outliers.

Users can employ the LAV method to find an estimator by choosing one of the available optimization solvers. Typically, Ipopt suffices for most scenarios.

Keywords

The function accepts the following keywords:

• iteration: Specifies the maximum number of iterations.
• tolerance: Specifies the allowed deviation from the optimal solution.
• bridge: Controls the bridging mechanism (default: false).
• name: Handles the creation of string names (default: true).
• verbose: Controls the output display, ranging from the default silent mode (0) to detailed output (3).

Additionally, users can modify variable names used for printing and writing through the keywords magnitude, angle, positive, and negative. For instance, users can choose magnitude = "V", angle = "θ", positive = "u", and negative = "v" to display equations in a more readable format.

Updates

If the AC model has not been created, the function will automatically trigger an update of the ac field within the PowerSystem composite type.

Returns

The function returns an instance of the AcStateEstimation type.

Example

using Ipopt

system, monitoring = ems("case14.h5", "monitoring.h5")

analysis = acLavStateEstimation(monitoring, Ipopt.Optimizer)

increment!(analysis::AcStateEstimation)

By solving the normal equation, this function computes the bus voltage magnitude and angle increments.

Updates

The function updates the residual, jacobian, and factorisation variables within the AcStateEstimation type. Using these results, it then computes and updates the increment variable within the same type. It should be used during the Gauss-Newton iteration loop before invoking the solve! function.

Returns

The function returns the maximum absolute increment value, which can be used to terminate the iteration loop of the Gauss-Newton method applied to solve the AC state estimation problem.

Example

system, monitoring = ems("case14.h5", "monitoring.h5")

analysis = gaussNewton(monitoring)
increment!(analysis)

solve!(analysis::AcStateEstimation)

By computing the bus voltage magnitudes and angles, the function solves the AC state estimation. Note that if the Gauss-Newton method is employed to obtain the WLS estimator, this function simply updates the state variables using the obtained increments.

Updates

The resulting bus voltage magnitudes and angles are stored in the voltage field of the AcStateEstimation type.

Examples

Solving the AC state estimation model and obtaining the WLS estimator:

system, monitoring = ems("case14.h5", "monitoring.h5")

analysis = gaussNewton(monitoring)
for iteration = 1:20
    stopping = increment!(analysis)
    if stopping < 1e-8
        break
    end
    solve!(analysis)
end

Solving the AC state estimation model and obtaining the LAV estimator:

using Ipopt

system, monitoring = ems("case14.h5", "monitoring.h5")

analysis = acLavStateEstimation(monitoring, Ipopt.Optimizer; verbose = 1)
solve!(analysis)

setInitialPoint!(analysis::AcStateEstimation)

The function sets the initial point of the AC state estimation to the values from the PowerSystem type.

Updates

The function modifies the voltage field of the AcStateEstimation type.

Example:

system, monitoring = ems("case14.h5", "monitoring.h5")

analysis = gaussNewton(monitoring)
stateEstimation!(analysis)

residualTest!(analysis; threshold = 1.0)

setInitialPoint!(analysis)
stateEstimation!(analysis)

setInitialPoint!(target::AcStateEstimation, source::Analysis)

The function initializes the AC state estimation based on results from the Analysis type, whether from an AC or DC analysis.

The function assigns the bus voltage magnitudes and angles in the target argument, using data from the source argument. This allows users to initialize AC state estimation as needed.

If source comes from a DC analysis, only the bus voltage angles are assigned in the target argument, while the bus voltage magnitudes remain unchanged.

Example

system, monitoring = ems("case14.h5", "monitoring.h5")

source = newtonRaphson(system)
powerFlow!(source)

target = gaussNewton(monitoring)

setInitialPoint!(target, source)
stateEstimation!(target)

stateEstimation!(analysis::AcStateEstimation;
    iteration, tolerance, power, current, verbose)

The function serves as a wrapper for solving AC state estimation and includes the functions:

• solve!,
• power!,
• current!.

Additionally, for the WLS model, it includes:

• increment!.

It computes bus voltage magnitudes and angles using the WLS or LAV model with the option to compute powers and currents.

Keywords

Users can use the following keywords:

• iteration: Specifies the maximum number of iterations (default for WLS model: 40).
• tolerance: Defines the tolerance for the iteration stopping criterion (default for WLS model: 1e-8).
• power: Enables the computation of powers (default: false).
• current: Enables the computation of currents (default: false).
• verbose: Controls the output display, ranging from the default silent mode (0) to detailed output (3).

For the WLS model, tolerance refers to the step size tolerance in the stopping criterion, whereas for the LAV model, it defines the allowed deviation from the optimal solution. If iteration and tolerance are not specified for the LAV model, the optimization solver settings are used.

Examples

Use the wrapper function to obtain the WLS estimator:

system, monitoring = ems("case14.h5", "monitoring.h5")

analysis = gaussNewton(monitoring)
stateEstimation!(analysis; tolerance = 1e-10, current = true, verbose = 3)

Use the wrapper function to obtain the LAV estimator:

using Ipopt

system, monitoring = ems("case14.h5", "monitoring.h5")

analysis = acLavStateEstimation(monitoring, Ipopt.Optimizer)
stateEstimation!(analysis; iteration = 30, power = true, verbose = 1)

pmuStateEstimation(monitoring::Measurement, [method = LU])

The function establishes the linear WLS model for state estimation with PMUs only. In this model, the vector of state variables contains bus voltages, given in rectangular coordinates.

Arguments

This function requires the Measurement type to establish the WLS state estimation model.

Moreover, the presence of the method parameter is not mandatory. To address the WLS state estimation method, users can opt to utilize factorization techniques to decompose the gain matrix, such as LU, QR, or LDLt especially when the gain matrix is symmetric. Opting for the Orthogonal method is advisable for a more robust solution in scenarios involving ill-conditioned data, particularly when substantial variations in variances are present.

If the user does not provide the method, the default method for solving the estimation model will be LU factorization.

Updates

If the AC model has not been created, the function will automatically trigger an update of the ac field within the PowerSystem composite type.

Returns

The function returns an instance of the PmuStateEstimation type.

Examples

Set up the PMU state estimation model to be solved using the default LU factorization:

system, monitoring = ems("case14.h5", "monitoring.h5")

analysis = pmuStateEstimation(monitoring)

Set up the PMU state estimation model to be solved using the orthogonal method:

system, monitoring = ems("case14.h5", "monitoring.h5")

analysis = pmuStateEstimation(monitoring, Orthogonal)

pmuLavStateEstimation(monitoring::Measurement, optimizer;
    iteration, tolerance, bridge, name, real, imag, positive, negative, verbose)

The function establishes the LAV model for state estimation with PMUs only. In this model, the vector of state variables contains bus voltages, given in rectangular coordinates.

Arguments

This function requires the Measurement type to establish the LAV state estimation model. The LAV method offers increased robustness compared to WLS, ensuring unbiasedness even in the presence of various measurement errors and outliers.

Users can employ the LAV method to find an estimator by choosing one of the available optimization solvers. Typically, Ipopt suffices for most scenarios.

Keywords

The function accepts the following keywords:

• iteration: Specifies the maximum number of iterations.
• tolerance: Specifies the allowed deviation from the optimal solution.
• bridge: Controls the bridging mechanism (default: false).
• name: Handles the creation of string names (default: true).
• verbose: Controls the output display, ranging from the default silent mode (0) to detailed output (3).

Additionally, users can modify variable names used for printing and writing through the keywords real, imag, positive, and negative. For instance, users can choose real = "Vr", imag = "Vi", positive = "u", and negative = "v" to display equations in a more readable format.

Updates

If the AC model has not been created, the function will automatically trigger an update of the ac field within the PowerSystem composite type.

Returns

The function returns an instance of the PmuStateEstimation type.

Example

using Ipopt

system, monitoring = ems("case14.h5", "monitoring.h5")

analysis = pmuLavStateEstimation(monitoring, Ipopt.Optimizer)

solve!(analysis::PmuStateEstimation)

By computing the bus voltage magnitudes and angles, the function solves the PMU state estimation model.

Updates

The resulting bus voltage magnitudes and angles are stored in the voltage field of the PmuStateEstimation type.

Examples

Solving the PMU state estimation model and obtaining the WLS estimator:

system, monitoring = ems("case14.h5", "monitoring.h5")

analysis = pmuStateEstimation(monitoring)
solve!(analysis)

Solving the PMU state estimation model and obtaining the LAV estimator:

using Ipopt

system, monitoring = ems("case14.h5", "monitoring.h5")

analysis = pmuLavStateEstimation(monitoring, Ipopt.Optimizer; verbose = 1)
solve!(analysis)

stateEstimation!(analysis::PmuStateEstimation; iteration, tolerance, power, verbose)

The function serves as a wrapper for solving PMU state estimation and includes the functions:

• solve!,
• power!,
• current!.

It computes bus voltage magnitudes and angles using the WLS or LAV model with the option to compute powers and currents.

Keywords

Users can use the following keywords:

• iteration: Specifies the maximum number of iterations for the LAV model.
• tolerance: Specifies the allowed deviation from the optimal solution for the LAV model.
• power: Enables the computation of powers (default: false).
• current: Enables the computation of currents (default: false).
• verbose: Controls the output display, ranging from the default silent mode (0) to detailed output (3).

If iteration and tolerance are not specified, the optimization solver settings are used.

Examples

Use the wrapper function to obtain the WLS estimator:

system, monitoring = ems("case14.h5", "monitoring.h5")

analysis = pmuStateEstimation(monitoring)
stateEstimation!(analysis; power = true, verbose = 3)

Use the wrapper function to obtain the LAV estimator:

using Ipopt

system, monitoring = ems("case14.h5", "monitoring.h5")

analysis = pmuLavStateEstimation(monitoring, Ipopt.Optimizer)
stateEstimation!(analysis; iteration = 30, tolerance = 1e-6, verbose = 3)

dcStateEstimation(monitoring::Measurement, [method = LU])

The function establishes the WLS model for DC state estimation, where the vector of state variables contains only bus voltage angles.

Arguments

This function requires the Measurement type to establish the WLS state estimation model.

Moreover, the presence of the method parameter is not mandatory. To address the WLS state estimation method, users can opt to utilize factorization techniques to decompose the gain matrix, such as LU, QR, or LDLt especially when the gain matrix is symmetric. Opting for the Orthogonal method is advisable for a more robust solution in scenarios involving ill-conditioned data, particularly when substantial variations in variances are present.

If the user does not provide the method, the default method for solving the estimation model will be LU factorization.

Updates

If the DC model was not created, the function will automatically initiate an update of the dc field within the PowerSystem composite type. Additionally, if the slack bus lacks an in-service generator, JuliaGrid considers it a mistake and defines a new slack bus as the first generator bus with an in-service generator in the bus type list.

Returns

The function returns an instance of the DcStateEstimation type.

Examples

Set up the DC state estimation model to be solved using the default LU factorization:

system, monitoring = ems("case14.h5", "monitoring.h5")

analysis = dcStateEstimation(monitoring)

Set up the DC state estimation model to be solved using the orthogonal method:

system, monitoring = ems("case14.h5", "monitoring.h5")

analysis = dcStateEstimation(monitoring, Orthogonal)

dcLavStateEstimation(monitoring::Measurement, optimizer;
    iteration, tolerance, bridge, name, angle, positive, negative, verbose)

The function establishes the LAV model for DC state estimation, where the vector of state variables contains only bus voltage angles.

Arguments

This function requires the Measurement type to establish the LAV state estimation model. The LAV method offers increased robustness compared to WLS, ensuring unbiasedness even in the presence of various measurement errors and outliers.

Users can employ the LAV method to find an estimator by choosing one of the available optimization solvers. Typically, Ipopt suffices for most scenarios.

Keywords

The function accepts the following keywords:

• iteration: Specifies the maximum number of iterations.
• tolerance: Specifies the allowed deviation from the optimal solution.
• bridge: Controls the bridging mechanism (default: false).
• name: Handles the creation of string names (default: true).
• verbose: Controls the output display, ranging from the default silent mode (0) to detailed output (3).

Additionally, users can modify variable names used for printing and writing through the keywords angle, positive, and negative. For instance, users can choose angle = "θ", positive = "u", and negative = "v" to display equations in a more readable format.

Updates

If the DC model was not created, the function will automatically initiate an update of the dc field within the PowerSystem composite type. Additionally, if the slack bus lacks an in-service generator, JuliaGrid considers it a mistake and defines a new slack bus as the first generator bus with an in-service generator in the bus type list.

Returns

The function returns an instance of the DcStateEstimation type.

Example

using Ipopt

system, monitoring = ems("case14.h5", "monitoring.h5")

analysis = dcLavStateEstimation(monitoring, Ipopt.Optimizer)

solve!(analysis::DcStateEstimation)

By computing the bus voltage angles, the function solves the DC state estimation model.

Updates

The resulting bus voltage angles are stored in the voltage field of the DcStateEstimation type.

Examples

Solving the DC state estimation model and obtaining the WLS estimator:

system, monitoring = ems("case14.h5", "monitoring.h5")

analysis = dcStateEstimation(monitoring)
solve!(analysis)

Solving the DC state estimation model and obtaining the LAV estimator:

using Ipopt

system, monitoring = ems("case14.h5", "monitoring.h5")

analysis = dcLavStateEstimation(monitoring, Ipopt.Optimizer; verbose = 1)
solve!(analysis)

stateEstimation!(analysis::DcStateEstimation; iteration, tolerance, power, verbose)

The function serves as a wrapper for solving DC state estimation and includes the functions:

• solve!,
• power!.

It computes bus voltage angles using the WLS or LAV model with the option to compute powers.

Keywords

Users can use the following keywords:

• iteration: Specifies the maximum number of iterations for the LAV model.
• tolerance: Specifies the allowed deviation from the optimal solution for the LAV model.
• power: Enables the computation of powers (default: false).
• verbose: Controls the output display, ranging from the default silent mode (0) to detailed output (3).

If iteration and tolerance are not specified, the optimization solver settings are used.

Examples

Use the wrapper function to obtain the WLS estimator:

system, monitoring = ems("case14.h5", "monitoring.h5")

analysis = dcStateEstimation(monitoring)
stateEstimation!(analysis; power = true, verbose = 3)

Use the wrapper function to obtain the LAV estimator:

using Ipopt

system, monitoring = ems("case14.h5", "monitoring.h5")

analysis = dcLavStateEstimation(monitoring, Ipopt.Optimizer)
stateEstimation!(analysis; iteration = 30, tolerance = 1e-6, verbose = 1)

chiTest(analysis::StateEstimation; confidence)

The function performs a Chi-squared bad data detection test. This test can be applied after obtaining WLS estimator.

Arguments

The abstract type StateEstimation can have the following subtypes:

• AcStateEstimation: Conducts bad data analysis within AC state estimation.
• PmuStateEstimation: Conducts bad data analysis within PMU state estimation.
• DcStateEstimation: Conducts bad data analysis within DC state estimation.

Keyword

The keyword confidence specifies the detection confidence level, with a default value of 0.95.

Returns

The function returns an instance of the ChiTest type.

Example

system, monitoring = ems("case14.h5", "monitoring.h5")

analysis = gaussNewton(monitoring)
stateEstimation!(analysis)

bad = chiTest(analysis; confidence = 0.96)

residualTest!(analysis::StateEstimation; threshold)

The function conducts bad data detection and identification using the largest normalized residual test, subsequently removing measurement outliers from the measurement set. It can be executed after obtaining WLS estimator.

Arguments

The abstract type StateEstimation can have the following subtypes:

• AcStateEstimation: Conducts bad data analysis within AC state estimation.
• PmuStateEstimation: Conducts bad data analysis within PMU state estimation.
• DcStateEstimation: Conducts bad data analysis within DC state estimation.

Keyword

The keyword threshold establishes the identification threshold. If the largest normalized residual surpasses this threshold, the measurement is flagged as bad data. The default threshold value is set to threshold = 3.0.

Updates

If bad data is detected, the function flags the corresponding measurement within the Measurement type as out-of-service.

Moreover, for DcStateEstimation and PmuStateEstimation types, the function removes the corresponding measurement from the coefficient matrix and mean vector. This facilitates direct progress to the function that solves the state estimation problem.

Returns

The function returns an instance of the ResidualTest type.

Example

system, monitoring = ems("case14.h5", "monitoring.h5")

analysis = dcStateEstimation(monitoring)
stateEstimation!(analysis)

outlier = residualTest!(analysis; threshold = 4.0)
stateEstimation!(analysis)