12. State Estimation#
State estimation (SE) determines the most likely operating state of a power system from a set of noisy, redundant measurements. While power flow computes the exact solution given precise inputs, state estimation works with real-world data that is inherently imprecise. It processes measurements such as bus voltage magnitudes, power injections, and branch flows, and produces the best estimate of all bus voltage magnitudes and angles.
ANDES provides two estimation algorithms. The weighted least squares (WLS) method minimizes the weighted sum of squared measurement residuals and is optimal under Gaussian noise. The least absolute value (LAV) method minimizes the weighted sum of absolute residuals and is robust to gross measurement errors. Redundancy in measurements provides both noise rejection and the ability to detect bad data through statistical tests.
This tutorial covers running state estimation in ANDES, constructing custom measurement sets, interpreting results, performing bad data detection with the chi-squared test, comparing WLS and LAV under bad data conditions, and handling multi-island systems.
# Reduce logging verbosity for PDF builds
import os
if os.environ.get('SPHINX_BUILD_PDF'):
import andes
_orig_config_logger = andes.config_logger
def _quiet_logger(stream_level=20, *args, **kwargs):
stream_level = max(stream_level, 30)
return _orig_config_logger(stream_level, *args, **kwargs)
andes.config_logger = _quiet_logger
12.1. Setup#
We begin by importing ANDES, NumPy, and the measurement framework. State estimation requires a converged power flow as the starting point, so we load a case and run power flow first.
%matplotlib inline
import numpy as np
import andes
from andes.se import Measurements
andes.config_logger(stream_level=20)
ss = andes.load(andes.get_case('ieee14/ieee14.json'))
ss.PFlow.run()
Working directory: "/home/docs/checkouts/readthedocs.org/user_builds/andes/checkouts/stable/docs/source/tutorials"
> Loaded generated Python code in "/home/docs/.andes/pycode".
Parsing input file "/home/docs/checkouts/readthedocs.org/user_builds/andes/envs/stable/lib/python3.11/site-packages/andes/cases/ieee14/ieee14.json"...
Input file parsed in 0.0021 seconds.
Connectivity check completed in 0.0004 seconds.
-> System connectivity check results:
No islanded bus detected.
System is interconnected.
Each island has a slack bus correctly defined and enabled.
System internal structure set up in 0.0211 seconds.
-> Power flow calculation
Numba: Off
Sparse solver: KLU
Solution method: NR method
Power flow initialized in 0.0036 seconds.
0: |F(x)| = 0.5605182134
1: |F(x)| = 0.006202200332
2: |F(x)| = 5.819382827e-06
3: |F(x)| = 6.957087684e-12
Converged in 4 iterations in 0.0022 seconds.
Report saved to "ieee14_out.txt" in 0.0012 seconds.
True
12.2. Running State Estimation#
The simplest way to run state estimation is to call ss.SE.run() with no arguments. ANDES will automatically generate a default set of measurements from the power flow solution, add Gaussian noise, and estimate the system state using the WLS algorithm.
ss.SE.run()
-> State Estimation
Method: WLS
Tolerance: 0.0001
Max iterations: 100
No measurements provided; using default set (all bus voltages + all bus injections, 42 measurements).
Adding angle reference at bus 1 for island with 14 buses.
WLS converged in 2 iterations, J = 16.2091
SE converged in 2 iterations in 0.0030 seconds.
-> SE Report
Converged: True
Iterations: 2
Objective J: 16.2091
Measurements: 43
States: 28
Redundancy: 1.54x
Max |dV|: 0.010248 pu
Max |da|: 0.0102194 rad
True
The SE report shows convergence status, the number of iterations, the objective function value \(J\), the number of measurements and states, and the measurement redundancy ratio. The estimation error (Max |dV| and Max |da|) compares the estimated state against the power flow solution.
After convergence, the estimated bus voltage magnitudes and angles are available through the v_est and a_est properties.
print(f"Converged: {ss.SE.converged}")
print(f"Estimated voltages (pu): {np.round(ss.SE.v_est, 4)}")
print(f"Estimated angles (rad): {np.round(ss.SE.a_est, 4)}")
Converged: True
Estimated voltages (pu): [1.0342 1.0314 1.0067 1.0103 1.016 1.0268 1.0219 1.0264 1.0199 1.0138
1.0171 1.0103 1.0086 1.0061]
Estimated angles (rad): [ 0. -0.0306 -0.0622 -0.0775 -0.0675 -0.1181 -0.0899 -0.0332 -0.1334
-0.1396 -0.1303 -0.135 -0.1405 -0.1677]
12.3. Custom Measurements#
In practice, not every bus and branch is instrumented. The Measurements class allows you to construct a custom measurement set that reflects the actual placement of meters in the system. Each measurement is characterized by its type, location, and standard deviation (sigma), which represents the expected measurement noise.
Four types of measurements are supported:
Method |
Measurement type |
Description |
|---|---|---|
|
Voltage magnitude |
RTU or PMU voltage measurements |
|
Voltage angle |
PMU angle measurements |
|
P and Q injection |
Bus power injection measurements |
|
P and Q flow |
Branch power flow measurements (from-end) |
The standard deviation sigma controls the weight of each measurement in the WLS objective. A smaller sigma means the estimator trusts that measurement more. Typical values are 0.005-0.02 for voltage measurements and 0.02-0.05 for power measurements.
m = Measurements(ss)
# Voltage meters at all buses (high accuracy)
m.add_bus_voltage(sigma=0.01)
# Power injection meters at all buses
m.add_bus_injection(sigma_p=0.02, sigma_q=0.03)
# Generate noisy measurement values from the PFlow solution
m.generate_from_pflow(seed=42)
print(f"Total measurements: {m.nm}")
print(f"State dimension: {2 * ss.Bus.n}")
print(f"Redundancy: {m.nm / (2 * ss.Bus.n):.1f}x")
Total measurements: 42
State dimension: 28
Redundancy: 1.5x
The redundancy ratio (measurements divided by states) indicates how much extra information is available. A ratio above 1.0 means the system is overdetermined, which is necessary for state estimation to work. Higher redundancy provides better noise rejection and more robust bad data detection. A redundancy of 2-3x is typical for real systems.
Pass the custom measurements to SE.run() to use them instead of the default set.
ss.SE.run(measurements=m)
-> State Estimation
Method: WLS
Tolerance: 0.0001
Max iterations: 100
Adding angle reference at bus 1 for island with 14 buses.
WLS converged in 2 iterations, J = 13.027
SE converged in 2 iterations in 0.0021 seconds.
-> SE Report
Converged: True
Iterations: 2
Objective J: 13.027
Measurements: 43
States: 28
Redundancy: 1.54x
Max |dV|: 0.0059106 pu
Max |da|: 0.00930009 rad
True
12.3.1. Selective Meter Placement#
You can also place measurements on specific buses or lines rather than the entire system. This is useful for modeling realistic SCADA configurations where only key substations have full instrumentation. The bus_idx and line_idx parameters accept a list of device indices.
m2 = Measurements(ss)
# Voltage at selected buses only
m2.add_bus_voltage(bus_idx=[1, 2, 3, 6, 8], sigma=0.005)
# PMU angle measurement at the slack bus
m2.add_bus_angle(bus_idx=[1], sigma=0.01)
# Injection meters at all buses
m2.add_bus_injection(sigma_p=0.02, sigma_q=0.03)
# Flow meters on selected lines
m2.add_line_flow(line_idx=['Line_1', 'Line_2', 'Line_3', 'Line_7', 'Line_8'],
sigma_p=0.02, sigma_q=0.03)
m2.generate_from_pflow(seed=42)
print(f"Total measurements: {m2.nm}")
ss.SE.run(measurements=m2)
Total measurements: 44
-> State Estimation
Method: WLS
Tolerance: 0.0001
Max iterations: 100
WLS converged in 2 iterations, J = 10.0405
SE converged in 2 iterations in 0.0028 seconds.
-> SE Report
Converged: True
Iterations: 2
Objective J: 10.0405
Measurements: 44
States: 28
Redundancy: 1.57x
Max |dV|: 0.0082004 pu
Max |da|: 0.0132314 rad
True
12.3.2. Low-Level add() Method#
All convenience wrappers delegate to the add() method, which accepts a model name and variable name. This method is useful when you need fine-grained control over individual measurements.
m3 = Measurements(ss)
# Equivalent to m3.add_bus_voltage(bus_idx=[1, 2], sigma=0.005)
m3.add('Bus', 'v', idx=[1, 2], sigma=0.005)
# Equivalent to m3.add_bus_angle(bus_idx=[1], sigma=0.01)
m3.add('Bus', 'a', idx=[1], sigma=0.01)
print(f"Measurements added: {m3.nm}")
Measurements added: 3
12.4. Zero-Noise Verification#
A useful sanity check is to verify that state estimation recovers the exact power flow solution when given perfect (noise-free) measurements. In this case, the measurement vector is set to \(z = h(x_{\text{true}})\) without any added noise. The WLS estimator should converge to the true state to within numerical precision.
This test validates the correctness of the measurement functions \(h(x)\) and the Jacobian \(H(x)\).
from andes.se import StaticEvaluator
m_exact = Measurements(ss)
m_exact.add_bus_voltage(sigma=0.01)
m_exact.add_bus_angle(sigma=0.01)
m_exact.add_bus_injection(sigma_p=0.02, sigma_q=0.03)
m_exact.add_line_flow(sigma_p=0.02, sigma_q=0.03)
# Set z = h(x_true) exactly (no noise)
m_exact.finalize()
ev = StaticEvaluator(ss, m_exact)
theta_true = np.array(ss.Bus.a.v, dtype=float)
Vm_true = np.array(ss.Bus.v.v, dtype=float)
m_exact.z = ev.h(theta_true, Vm_true)
ss.SE.run(measurements=m_exact)
v_err = np.max(np.abs(ss.SE.v_est - Vm_true))
a_err = np.max(np.abs(ss.SE.a_est - theta_true))
print(f"Max voltage error: {v_err:.2e} pu")
print(f"Max angle error: {a_err:.2e} rad")
Max voltage error: 0.00e+00 pu
Max angle error: 0.00e+00 rad
-> State Estimation
Method: WLS
Tolerance: 0.0001
Max iterations: 100
WLS converged in 1 iterations, J = 0
SE converged in 1 iterations in 0.0022 seconds.
-> SE Report
Converged: True
Iterations: 1
Objective J: 0
Measurements: 96
States: 28
Redundancy: 3.43x
Max |dV|: 0 pu
Max |da|: 0 rad
The errors are at the level of machine precision (\(\sim 10^{-11}\)), confirming that the measurement model and solver are implemented correctly.
12.5. Bad Data Detection#
The chi-squared (\(\chi^2\)) test is a standard method for detecting the presence of bad data in the measurement set. Under the assumption that measurement errors are Gaussian, the weighted sum of squared residuals \(J(\hat{x}) = [z - h(\hat{x})]^T W [z - h(\hat{x})]\) follows a chi-squared distribution with \(m - n\) degrees of freedom, where \(m\) is the number of measurements and \(n\) is the number of states.
If \(J\) exceeds the chi-squared threshold at the desired confidence level, the test rejects the null hypothesis that all measurements are normal, indicating possible bad data.
# Run SE with normal measurements
m_normal = Measurements(ss)
m_normal.add_bus_voltage(sigma=0.01)
m_normal.add_bus_injection(sigma_p=0.02, sigma_q=0.03)
m_normal.generate_from_pflow(seed=42)
ss.SE.run(measurements=m_normal)
passed, J, threshold, dof = ss.SE.chi_squared_test()
print(f"Chi-squared test: {'PASSED' if passed else 'FAILED'}")
print(f" J = {J:.4f}, threshold = {threshold:.4f}, dof = {dof}")
Chi-squared test: PASSED
J = 13.0270, threshold = 24.9958, dof = 15
-> State Estimation
Method: WLS
Tolerance: 0.0001
Max iterations: 100
Adding angle reference at bus 1 for island with 14 buses.
WLS converged in 2 iterations, J = 13.027
SE converged in 2 iterations in 0.0021 seconds.
-> SE Report
Converged: True
Iterations: 2
Objective J: 13.027
Measurements: 43
States: 28
Redundancy: 1.54x
Max |dV|: 0.0059106 pu
Max |da|: 0.00930009 rad
Now we inject a gross error into one measurement (a 10-sigma outlier) and observe the effect on the chi-squared test.
# Inject a gross error
m_bad = Measurements(ss)
m_bad.add_bus_voltage(sigma=0.01)
m_bad.add_bus_injection(sigma_p=0.02, sigma_q=0.03)
m_bad.generate_from_pflow(seed=42)
# Corrupt one voltage measurement with a 10-sigma outlier
m_bad.z[0] += 10 * m_bad.sigma[0]
ss.SE.run(measurements=m_bad)
passed, J, threshold, dof = ss.SE.chi_squared_test()
print(f"Chi-squared test: {'PASSED' if passed else 'FAILED'}")
print(f" J = {J:.4f}, threshold = {threshold:.4f}, dof = {dof}")
Chi-squared test: FAILED
J = 109.0202, threshold = 24.9958, dof = 15
-> State Estimation
Method: WLS
Tolerance: 0.0001
Max iterations: 100
Adding angle reference at bus 1 for island with 14 buses.
WLS converged in 3 iterations, J = 109.02
SE converged in 3 iterations in 0.0033 seconds.
-> SE Report
Converged: True
Iterations: 3
Objective J: 109.02
Measurements: 43
States: 28
Redundancy: 1.54x
Max |dV|: 0.0148127 pu
Max |da|: 0.0046609 rad
The chi-squared test correctly detects the presence of bad data. In practice, after detection, further analysis (such as the largest normalized residual test) would be used to identify and remove the specific bad measurement.
12.6. Robust Estimation with LAV#
The weighted least squares estimator is optimal under Gaussian noise but is sensitive to gross measurement errors. The least absolute value (LAV) estimator provides a robust alternative by minimizing the weighted sum of absolute residuals rather than squared residuals:
This L1-norm formulation assigns less influence to measurements with large residuals compared to the L2-norm used by WLS. The LAV problem is solved via iteratively reweighted least squares (IRLS), where each iteration solves a WLS subproblem with weights inversely proportional to the current residual magnitudes. Because IRLS converges more slowly than Gauss-Newton, a higher iteration limit is generally recommended.
The lav algorithm function can be passed to SE.run() through the algorithm parameter. The following example compares WLS and LAV on the corrupted measurement set from the previous section.
from andes.se import lav
# True state from the power flow solution
Vm_true = np.array(ss.Bus.v.v, dtype=float)
theta_true = np.array(ss.Bus.a.v, dtype=float)
# WLS errors (from the run in the previous section on m_bad)
wls_v_err = np.max(np.abs(ss.SE.v_est - Vm_true))
wls_a_err = np.max(np.abs(ss.SE.a_est - theta_true))
# Run LAV on the same corrupted measurements
ss.SE.run(measurements=m_bad, algorithm=lav)
lav_v_err = np.max(np.abs(ss.SE.v_est - Vm_true))
lav_a_err = np.max(np.abs(ss.SE.a_est - theta_true))
print(f"{'Metric':<20} {'WLS':>12} {'LAV':>12}")
print(f"{'':-<44}")
print(f"{'Max |dV| (pu)':<20} {wls_v_err:>12.6f} {lav_v_err:>12.6f}")
print(f"{'Max |da| (rad)':<20} {wls_a_err:>12.6f} {lav_a_err:>12.6f}")
Metric WLS LAV
--------------------------------------------
Max |dV| (pu) 0.014813 0.004423
Max |da| (rad) 0.004661 0.005232
-> State Estimation
Method: LAV
Tolerance: 0.0001
Max iterations: 100
LAV converged in 34 iterations, J = 27.1812
SE converged in 34 iterations in 0.0223 seconds.
-> SE Report
Converged: True
Iterations: 34
Objective J: 27.1812
Measurements: 43
States: 28
Redundancy: 1.54x
Max |dV|: 0.00442327 pu
Max |da|: 0.00523227 rad
The LAV estimator produces smaller estimation errors than WLS in the presence of gross measurement errors. The L1-norm formulation effectively downweights the corrupted measurement, preventing it from biasing the state estimate. This robustness comes at the cost of slower convergence, as the IRLS procedure typically requires more iterations than the standard Gauss-Newton method.
Note that the chi-squared test is defined for the WLS objective function and is not applicable after an LAV run. ANDES will issue a warning if chi_squared_test() is called following LAV estimation.
12.7. Multi-Island Systems#
In systems with multiple electrical islands (disconnected components), each island requires its own voltage angle reference because absolute angles are not observable from power measurements alone. ANDES handles this automatically by detecting islands through connectivity analysis and adding a pseudo-measurement for the angle reference at the slack bus of each island.
The Kundur two-area system provides a convenient example when configured as two separate islands.
ss_island = andes.load(andes.get_case('kundur/kundur_islands.xlsx'),
default_config=True)
ss_island.PFlow.run()
print(f"Number of islands: {len(ss_island.Bus.island_sets)}")
Number of islands: 2
Working directory: "/home/docs/checkouts/readthedocs.org/user_builds/andes/checkouts/stable/docs/source/tutorials"
> Reloaded generated Python code of module "pycode".
Parsing input file "/home/docs/checkouts/readthedocs.org/user_builds/andes/envs/stable/lib/python3.11/site-packages/andes/cases/kundur/kundur_islands.xlsx"...
Input file parsed in 0.0322 seconds.
Connectivity check completed in 0.0001 seconds.
-> System connectivity check results:
No islanded bus detected.
System contains 2 island(s).
Each island has a slack bus correctly defined and enabled.
System internal structure set up in 0.0199 seconds.
-> Power flow calculation
Numba: Off
Sparse solver: KLU
Solution method: NR method
Power flow initialized in 0.0031 seconds.
0: |F(x)| = 14.9282832
1: |F(x)| = 3.670045871
2: |F(x)| = 0.1542941496
3: |F(x)| = 0.0003701017904
4: |F(x)| = 2.966633161e-09
Converged in 5 iterations in 0.0025 seconds.
Report saved to "kundur_islands_out.txt" in 0.0007 seconds.
m_island = Measurements(ss_island)
m_island.add_bus_voltage(sigma=0.01)
m_island.add_bus_injection(sigma_p=0.02, sigma_q=0.03)
m_island.generate_from_pflow(seed=42)
ss_island.SE.run(measurements=m_island)
-> State Estimation
Method: WLS
Tolerance: 0.0001
Max iterations: 100
Adding angle reference at bus 1 for island with 5 buses.
Adding angle reference at bus 3 for island with 5 buses.
WLS converged in 2 iterations, J = 8.07042
SE converged in 2 iterations in 0.0022 seconds.
-> SE Report
Converged: True
Iterations: 2
Objective J: 8.07042
Measurements: 32
States: 20
Redundancy: 1.60x
Max |dV|: 0.00799312 pu
Max |da|: 0.00501073 rad
True
The log shows that ANDES added angle references for both islands automatically. Without these references, the WLS gain matrix would be singular and the estimator would fail.
12.8. Configuration#
The SE routine provides several configuration options accessible through ss.SE.config.
ss.SE.config
OrderedDict([('linsolve', 0),
('tol', 0.0001),
('max_iter', 100),
('flat_start', 0),
('report', 1)])
Option |
Default |
Description |
|---|---|---|
|
1e-4 |
Convergence tolerance for the SE iteration |
|
100 |
Maximum number of iterations |
|
0 |
Use flat start (1.0 pu, 0 rad) instead of PFlow solution |
|
1 |
Print result summary after estimation |
The flat start option is useful for testing robustness of convergence when the power flow solution is not available as an initial guess.
# Run with flat start
ss.SE.config.flat_start = 1
ss.SE.run(measurements=m_normal)
# Reset to default
ss.SE.config.flat_start = 0
-> State Estimation
Method: WLS
Tolerance: 0.0001
Max iterations: 100
WLS converged in 3 iterations, J = 13.027
SE converged in 3 iterations in 0.0026 seconds.
-> SE Report
Converged: True
Iterations: 3
Objective J: 13.027
Measurements: 43
States: 28
Redundancy: 1.54x
Max |dV|: 0.00591063 pu
Max |da|: 0.00930006 rad
12.9. Cleanup#
!andes misc -C
"/home/docs/checkouts/readthedocs.org/user_builds/andes/checkouts/stable/docs/source/tutorials/kundur_islands_out.txt" removed.
"/home/docs/checkouts/readthedocs.org/user_builds/andes/checkouts/stable/docs/source/tutorials/ieee14_out.txt" removed.
12.10. Next Steps#
Eigenvalue Analysis - Small-signal stability assessment
Contingency Analysis - N-1 contingency screening