Source code for pydynamicestimator.tests.test_avr_algebraic

# Created: 2026-06-02
# (c) Copyright 2024 ETH Zurich, Milos Katanic
# https://doi.org/10.5905/ethz-1007-842
#
# Licensed under the GNU General Public License v3.0;
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at:
#
#     https://www.gnu.org/licenses/gpl-3.0.en.html
#
# This software is distributed "AS IS", WITHOUT WARRANTY OF ANY KIND,
# express or implied. See the License for specific language governing
# permissions and limitations under the License.

"""Controller-side parity test for the device-private algebraic-equation mechanism.

The generalized AVR strategy protocol lets an exciter expose its field voltage
``Efd`` as *either* a differential state (pure-lag exciter) *or* a device-private
algebraic variable (direct-feedthrough / lead-lag exciter). This test exercises
the latter via ``AVRKundur`` -- the Kundur transducer+lead-lag AVR with
``Efd`` as the lead's algebraic output -- and validates it two ways:

1. The recovered ``Efd`` (in ``dae.y``) exactly satisfies its defining equation
   along the trajectory (internal consistency of the private-algebraic solve).
2. The lead-lag (algebraic ``Efd``) reproduces ``AVRKundur_Filter`` -- the same
   exciter with a parasitic output filter that fakes ``Efd`` into a state -- in
   the ``Tfd -> 0`` singular-perturbation limit. The ``avr_filter`` fixture uses
   a small ``Tfd``, so the two systems must agree to ``O(Tfd)``.

See docs/algebraic_equations_design.md and CHANGELOG_algebraic_equations.md.
"""

from pathlib import Path

import numpy as np

from pydynamicestimator.config import config
from pydynamicestimator.run import run

FIXTURE_ROOT = Path(__file__).parent / "fixtures"

_COMMON = dict(
    system_root=FIXTURE_ROOT,
    fn=50,
    Sb=100,
    ts=0.005,
    te=0.02,
    T_start=0.0,
    T_end=6.0,
    int_scheme="backward",
    int_scheme_sim="idas",
    init_error_diff=1,
    init_error_alg=True,
    plot=False,
    plot_voltage=False,
    plot_diff=False,
    proc_noise_alg=1e-3,
    proc_noise_diff=1e-4,
    filter="iekf",
    log_level="ERROR",
    incl_lim=False,
    line_dyn=False,
    skip_est=True,
    print_power_flow=False,
)


[docs] def _machine(which="sim"): """The synchronous-machine device from the most recent run. ``run`` resets ``system`` (clear_module) on entry, so this must be called immediately after a run and before the next one. ``which`` selects the simulation or estimation device list. """ from pydynamicestimator import system as _sys devices = _sys.device_list_sim if which == "sim" else _sys.device_list_est return next( d for d in devices if getattr(d, "_name", "").startswith("Synchronous") )
[docs] def test_avr_leadlag_constraint_holds(): """The lead-lag's algebraic field voltage must satisfy its defining equation along the whole trajectory: the recovered Efd (in y_full) equals Vl*(1 - TA/TB) + (TA/TB)*KA*(Vf_ref - Vtr) evaluated on the states.""" _, sim = run(config.updated(testsystemfile="avr_leadlag", **_COMMON)) m = _machine() assert sim.n_priv == 1 # Efd only x = np.array(sim.x_full) y = np.array(sim.y_full) efd = y[m.Efd[0]] # algebraic, recovered from the DAE solve Vl = x[m.Vl[0]] Vtr = x[m.Vtr[0]] TA, TB, KA, Vfr = m.TA[0], m.TB[0], m.KA[0], m.Vf_ref[0] expr = Vl * (1 - TA / TB) + (TA / TB) * KA * (Vfr - Vtr) assert np.abs(efd - expr).max() < 1e-6, ( f"Efd violates its lead-lag defining equation by " f"{np.abs(efd - expr).max():.3e} (expected < 1e-6)." )
[docs] def test_avr_leadlag_matches_filtered_limit(): """AVRKundur (Efd algebraic) must reproduce AVRKundur_Filter (Efd a state behind a fast parasitic filter) in the Tfd -> 0 limit. The avr_filter fixture uses Tfd = 1e-3, so the machine response must agree to O(Tfd).""" _, sim_f = run(config.updated(testsystemfile="avr_filter", **_COMMON)) mf = _machine() xf = np.array(sim_f.x_full) efd_f = xf[mf.Efd[0]] # Efd is a STATE in the filtered model delta_f = xf[mf.delta[0]] omega_f = xf[mf.omega[0]] _, sim_l = run(config.updated(testsystemfile="avr_leadlag", **_COMMON)) ml = _machine() xl = np.array(sim_l.x_full) yl = np.array(sim_l.y_full) efd_l = yl[ml.Efd[0]] # Efd is ALGEBRAIC in the lead-lag model delta_l = xl[ml.delta[0]] omega_l = xl[ml.omega[0]] # Structural difference: removing the filter drops one differential state # (Efd) and adds one private algebraic. assert sim_f.n_priv == 0 and sim_l.n_priv == 1 assert sim_l.nx == sim_f.nx - 1 assert sim_l.ny == sim_f.ny + 1 n = min(xf.shape[1], xl.shape[1]) assert np.isfinite(efd_l).all() # Machine dynamics (slow states) are essentially identical... assert np.abs(delta_f[:n] - delta_l[:n]).max() < 1e-3 assert np.abs(omega_f[:n] - omega_l[:n]).max() < 1e-4 # ...and the field voltage agrees to the filter-time-constant order (Tfd=1e-3). assert np.abs(efd_f[:n] - efd_l[:n]).max() < 5e-3, ( f"Algebraic Efd differs from the Tfd->0 filtered Efd by " f"{np.abs(efd_f[:n] - efd_l[:n]).max():.3e} (expected O(Tfd), < 5e-3)." )
[docs] def test_avr_leadlag_estimation_runs_and_tracks(): """An exciter whose field voltage Efd is a device-private ALGEBRAIC variable runs through the full IEKF estimator (qcall noise overlay, init-from-simulation seeding of the private, save/restore) and tracks the simulated truth as well as the eliminated/state-Efd models. Verifies the controller-side private rides the same estimation plumbing as the machine-side SP-DAE privates.""" cfg = dict(_COMMON) cfg["skip_est"] = False est, sim = run(config.updated(testsystemfile="avr_leadlag", **cfg)) m = _machine("est") assert est.n_priv == 1 # algebraic Efd xe = np.array(est.x_full) ye = np.array(est.y_full) xs = np.array(sim.x_full) assert np.isfinite(xe).all() and np.isfinite(ye).all() assert np.isfinite(ye[m.Efd[0]]).all() # the algebraic Efd is recovered n = min(xe.shape[1], xs.shape[1]) track_err = np.abs(xe[:, :n] - xs[:, :n]) # Same accuracy band as the eliminated/state-Efd estimators. assert track_err.mean() < 0.05, f"mean tracking error {track_err.mean():.4f}" assert track_err.max() < 0.5, f"max tracking error {track_err.max():.4f}"