# Created: 2024-12-01
# Last Modified: 2026-03-19
# (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.
#
# The code is based on the publication: Katanic, M., Lygeros, J., Hug, G.: Recursive
# dynamic state estimation for power systems with an incomplete nonlinear DAE model.
# IET Gener. Transm. Distrib. 18, 3657–3668 (2024). https://doi.org/10.1049/gtd2.13308
# The full paper version is available at: https://arxiv.org/abs/2305.10065v2
# See full metadata at: README.md
# For inquiries, contact: mkatanic@ethz.ch
from __future__ import annotations # Postponed type evaluation
from typing import TYPE_CHECKING, Tuple
from pydynamicestimator.devices.device import DeviceRect, sin, cos
from pydynamicestimator.devices.avr import AVR, IEEEDC1A
from pydynamicestimator.devices.governor import Governor, TGOV1
from pydynamicestimator.devices.pss import PSS
from pydynamicestimator.devices.shaft import Shaft, SingleMass
if TYPE_CHECKING:
from pydynamicestimator.system import Dae
import casadi as ca
import numpy as np
[docs]
class Synchronous(DeviceRect):
"""Metaclass for SG in rectangular coordinates with TGOV1 governor and pluggable AVR"""
def __init__(
self,
avr: AVR = None,
governor: Governor = None,
pss: PSS = None,
shaft: Shaft = None,
):
super().__init__()
# Store the shaft strategy (default to a single rigid mass -- the classic
# swing equation -- for backward compatibility). A multi-mass torsional
# shaft adds further rotor masses while keeping the generator mass on the
# canonical 'delta'/'omega' names.
self._shaft: Shaft = shaft or SingleMass()
# Store the AVR strategy (default to IEEEDC1A for backward compatibility)
self._avr: AVR = avr or IEEEDC1A()
# Store the governor strategy (default to TGOV1 for backward compatibility)
self._governor: Governor = governor or TGOV1()
# Store the PSS strategy. There is NO default: when None, the machine
# supplies no stabilizing signal to its AVR (pss_signal returns 0).
self._pss: PSS = pss
# --- Machine params ---
self._params.update(
{
"fn": 50,
"H": 30,
"R_s": 0.0,
"x_q": 0.2,
"x_d": 0.2,
"x_l": 0.1,
"D": 0.0,
"f": 0.01,
}
)
# --- Shaft params (from strategy; empty for the single-mass default) ---
self._params.update(self._shaft.params())
# --- Governor params (from strategy) ---
self._params.update(self._governor.params())
# --- AVR params (from strategy) ---
self._params.update(self._avr.params())
# --- PSS params (from strategy, if any) ---
if self._pss is not None:
self._params.update(self._pss.params())
self._descr.update(
{
"H": "inertia constant",
"D": "rotor damping",
"fn": "rated frequency",
"bus": "bus id",
"gen": "static generator id",
"R_s": "stator resistance",
"x_d": "reactance in d axis",
"x_q": "reactance in q axis",
"x_l": "leakage / stator-series reactance",
"f": "rotor friction coefficient",
}
)
# delta/omega (and any extra rotor masses) are described by the shaft.
self._descr.update(self._shaft.descriptions())
self._descr.update(self._governor.descriptions())
self._descr.update(self._avr.descriptions())
if self._pss is not None:
self._descr.update(self._pss.descriptions())
# params
# SG
self.fn = np.array([], dtype=float)
self.H = np.array([], dtype=float)
self.R_s = np.array([], dtype=float)
self.x_d = np.array([], dtype=float)
self.x_q = np.array([], dtype=float)
self.x_l = np.array([], dtype=float)
self.D = np.array([], dtype=float)
self.f = np.array([], dtype=float)
# Shaft param arrays (dynamic from strategy; none for the single-mass shaft)
for param_name in self._shaft.params():
setattr(self, param_name, np.array([], dtype=float))
# Governor param arrays (dynamic from strategy)
for param_name in self._governor.params():
setattr(self, param_name, np.array([], dtype=float))
# AVR param arrays (dynamic from strategy)
for param_name in self._avr.params():
setattr(self, param_name, np.array([], dtype=float))
# PSS param arrays (dynamic from strategy, if any)
if self._pss is not None:
for param_name in self._pss.params():
setattr(self, param_name, np.array([], dtype=float))
# --- Rotor/shaft states (from strategy) ---
# Registered FIRST so the default single-mass shaft yields the same
# leading [delta, omega, ...] ordering as before (preserves the
# positional baselines). Every shaft MUST expose the generator mass as
# the differential states 'delta' and 'omega'; a multi-mass shaft adds
# further rotor masses (delta_<name>/omega_<name>) after them.
self.ns = 0
shaft_states = self._shaft.states()
assert "delta" in shaft_states and "omega" in shaft_states, (
"Shaft must expose the generator mass as the states 'delta' and 'omega'"
)
self.ns += len(shaft_states)
self.states.extend(shaft_states)
self.units.extend(self._shaft.units())
for state_name in shaft_states:
setattr(self, state_name, np.array([], dtype=float))
# --- Governor states + private algebraics (from strategy) ---
# Registered BEFORE the AVR so the state-vector ordering stays
# [delta, omega, <governor>, <avr>, <machine flux>] -- identical to the
# previous hardcoded [delta, omega, psv, pm, ...] when the default TGOV1
# is used (preserves positional baseline comparisons). 'pm' is the
# mechanical-power coupling into the swing equation; it may be a state
# (TGOV1) or a private algebraic (read via self.var_sym("pm")).
gov_states = self._governor.states()
gov_algebs = self._governor.algebs()
assert "pm" in set(gov_states) | set(gov_algebs), (
"Governor must expose 'pm' as either a state or a private algebraic"
)
self.ns += len(gov_states)
self.states.extend(gov_states)
self.units.extend(self._governor.units())
for state_name in gov_states:
setattr(self, state_name, np.array([], dtype=float))
if gov_algebs:
self._algebs_int.extend(gov_algebs)
self._algebs_int_units.update(self._governor.algebs_units())
self._algebs_int_noise.update(self._governor.algebs_noise())
self._algebs_int_x0.update(self._governor.algebs_x0())
for algeb_name in gov_algebs:
setattr(self, algeb_name, np.array([], dtype=float))
# --- AVR states (from strategy) ---
avr_states = self._avr.states()
avr_algebs = self._avr.algebs()
# 'Efd' is the field-voltage coupling variable; it may be a differential
# state (pure-lag exciter, e.g. IEEEDC1A) or a device-private algebraic
# (direct-feedthrough exciter, e.g. AVRKundur). Either is fine --
# the machine reads it through self.var_sym("Efd").
assert "Efd" in set(avr_states) | set(avr_algebs), (
"AVR must expose 'Efd' as either a state or a private algebraic"
)
self.ns += len(avr_states)
self.states.extend(avr_states)
self.units.extend(self._avr.units())
for state_name in avr_states:
setattr(self, state_name, np.array([], dtype=float))
# --- AVR private algebraic variables (from strategy) ---
# Registered onto the host's device-private-algebraic machinery
# (_algebs_int + metadata), exactly like a machine-declared private
# (e.g. the SP-DAE stator currents). Empty for pure-lag exciters, so
# this is a no-op for every existing AVR.
if avr_algebs:
self._algebs_int.extend(avr_algebs)
self._algebs_int_units.update(self._avr.algebs_units())
self._algebs_int_noise.update(self._avr.algebs_noise())
self._algebs_int_x0.update(self._avr.algebs_x0())
for algeb_name in avr_algebs:
setattr(self, algeb_name, np.array([], dtype=float))
# --- PSS states + private algebraics (from strategy, if any) ---
# Registered after the AVR. 'Vs' is the stabilizing signal summed into
# the AVR error; it may be a state or a private algebraic and is read by
# the AVR through self.pss_signal(dae). Absent by default (no PSS), so
# this whole block is a no-op for every machine without one.
if self._pss is not None:
pss_states = self._pss.states()
pss_algebs = self._pss.algebs()
assert "Vs" in set(pss_states) | set(pss_algebs), (
"PSS must expose 'Vs' as either a state or a private algebraic"
)
self.ns += len(pss_states)
self.states.extend(pss_states)
self.units.extend(self._pss.units())
for state_name in pss_states:
setattr(self, state_name, np.array([], dtype=float))
if pss_algebs:
self._algebs_int.extend(pss_algebs)
self._algebs_int_units.update(self._pss.algebs_units())
self._algebs_int_noise.update(self._pss.algebs_noise())
self._algebs_int_x0.update(self._pss.algebs_x0())
for algeb_name in pss_algebs:
setattr(self, algeb_name, np.array([], dtype=float))
self._states_noise.update(self._shaft.states_noise())
self._states_noise.update(self._governor.states_noise())
self._states_noise.update(self._avr.states_noise())
if self._pss is not None:
self._states_noise.update(self._pss.states_noise())
self._states_init_error.update(self._shaft.states_init_error())
self._states_init_error.update(self._governor.states_init_error())
self._states_init_error.update(self._avr.states_init_error())
if self._pss is not None:
self._states_init_error.update(self._pss.states_init_error())
self._x0.update(self._shaft.x0())
self._x0.update(self._governor.x0())
self._x0.update(self._avr.x0())
if self._pss is not None:
self._x0.update(self._pss.x0())
# Set points (from strategies)
self._setpoints.update(self._governor.setpoints())
self._setpoints.update(self._avr.setpoints())
if self._pss is not None:
self._setpoints.update(self._pss.setpoints())
for sp_name in self._governor.setpoints():
setattr(self, sp_name, np.array([], dtype=float))
for sp_name in self._avr.setpoints():
setattr(self, sp_name, np.array([], dtype=float))
if self._pss is not None:
for sp_name in self._pss.setpoints():
setattr(self, sp_name, np.array([], dtype=float))
self.properties.update({"fplot": True})
[docs]
def gcall(self, dae: Dae, i_d: ca.SX, i_q: ca.SX) -> None:
# algebraic equations (current balance in rectangular coordinates) + scale the current back to the grid reference power
dae.g[self.vre] -= (
self.Sn
/ dae.Sb
* (i_d * sin(dae.x[self.delta]) + i_q * cos(dae.x[self.delta]))
)
dae.g[self.vim] -= (
self.Sn
/ dae.Sb
* (i_d * -cos(dae.x[self.delta]) + i_q * sin(dae.x[self.delta]))
)
[docs]
def _omega_ref(self, dae: Dae) -> ca.SX:
"""Per-machine reference-frame frequency: the nominal ``dae.omega_net``
(= 1 p.u.) by default, or the estimated per-bus reference frequency when
``dae.omega_ref_buses`` is present. Shared by the rotor seam (was
duplicated inside every electromagnetic method)."""
omega_ref_vec = ca.SX.ones(self.n, 1) * dae.omega_net
if getattr(dae, "omega_ref_buses", None) is not None:
omega_list = []
for k in range(self.n):
bus_label = str(self.bus[k])
bus_idx = dae.grid.idx_bus[bus_label]
omega_list.append(dae.omega_ref_buses[bus_idx])
omega_ref_vec = ca.vertcat(*omega_list)
return omega_ref_vec
[docs]
def rotor(self, dae: Dae, Pe: ca.SX) -> None:
"""Rotor mechanics seam: delegate the swing (rotor-motion) equations to
the pluggable shaft strategy. The default :class:`SingleMass` writes the
classic single-mass swing equation (reading the governor port
``var_sym("pm")``, the air-gap power ``Pe`` and the reference frequency);
a multi-mass (torsional) shaft adds further rotor masses while keeping the
generator mass on the ``delta`` / ``omega`` names, so the electromagnetic
models, the network injection and the controllers are untouched. See
docs/algebraic_equations_design.md (§9.6/§9.7)."""
self._shaft.fgcall(self, dae, Pe)
[docs]
def fgcall(self, dae: Dae) -> None:
"""Template orchestration shared by all machine models: the EM subclass
defines stator currents + flux dynamics and returns the air-gap power;
the base wires the rotor, the controller ports, and the network."""
i_d, i_q, Pe = self.electromagnetic(dae)
self.rotor(dae, Pe)
self.governor_fcall(dae)
self.avr_fcall(dae)
self.pss_fcall(dae)
self.gcall(dae, i_d, i_q)
[docs]
def electromagnetic(self, dae: Dae):
"""Subclass hook: define the stator currents (i_d, i_q), write the
flux/emf differential equations, and return ``(i_d, i_q, Pe)`` where
``Pe`` is the air-gap power/torque consumed by :meth:`rotor`."""
raise NotImplementedError
[docs]
def governor_fcall(self, dae: Dae) -> None:
"""Delegate the governor's equations to the strategy object. The strategy
writes its differential equations into ``dae.f`` and, if it declares
private algebraics (e.g. an algebraic 'pm'), their residuals into
``dae.g`` -- hence a single ``fgcall``."""
self._governor.fgcall(self, dae)
[docs]
def tgov1(self, dae: Dae) -> None:
"""Deprecated: use governor_fcall instead."""
self.governor_fcall(dae)
[docs]
def avr_fcall(self, dae: Dae) -> None:
"""Delegate the AVR's equations to the strategy object. The strategy
writes its differential equations into ``dae.f`` and, if it declares
private algebraics (e.g. an algebraic 'Efd'), their residuals into
``dae.g`` -- hence a single ``fgcall``."""
self._avr.fgcall(self, dae)
[docs]
def var_sym(self, dae: Dae, name: str) -> ca.SX:
"""Resolve a registered variable to its symbol, wherever it lives.
A coupling variable such as 'Efd' may be a differential state (in
``dae.x``) or a device-private algebraic (in ``dae.y``), depending on
the AVR strategy. This resolver lets the machine's electromagnetic
equations reference it without knowing which -- the keystone that makes
state-vs-algebraic exciters interchangeable.
"""
if name in self._algebs_int:
return dae.y[getattr(self, name)]
return dae.x[getattr(self, name)]
[docs]
def pss_fcall(self, dae: Dae) -> None:
"""Delegate the PSS's equations to the strategy object (no-op if there is
no PSS). The strategy writes its differential equations into ``dae.f`` and,
if it declares private algebraics (e.g. an algebraic 'Vs'), their residuals
into ``dae.g``."""
if self._pss is not None:
self._pss.fgcall(self, dae)
[docs]
def pss_signal(self, dae: Dae) -> ca.SX:
"""The supplementary stabilizing signal 'Vs' summed into the AVR voltage
error -- the host-mediated PSS->AVR coupling. Returns the 'Vs' symbol
(state or private algebraic, via :meth:`var_sym`) when a PSS is present,
else 0 (no PSS => no signal, byte-identical to the unstabilized AVR)."""
if self._pss is None:
return 0
return self.var_sym(dae, "Vs")
[docs]
def ieeedc1a(self, dae: Dae) -> None:
"""Deprecated: use avr_fcall instead."""
self.avr_fcall(dae)
[docs]
class SynchronousTransient(Synchronous):
r"""
Transient two-axis SG with TGOV1 governor and IEEEDC1A AVR
**Rotor Dynamics**
.. math::
\dot{\delta} &= 2 \pi f_n \Delta \omega \\
\Delta \dot{\omega} &= \frac{1}{2 H} \left( P_m - E_d I_d - E_q I_q + (X_q' - X_d') I_d I_q - D \Delta \omega - f (\Delta \omega + 1) \right)
**Electromagnetic Equations**
.. math::
\dot{E}_q &= \frac{1}{T_{d'}} \left( -E_q + E_f + (X_d - X_d') I_d \right)\\
\dot{E}_d &= \frac{1}{T_{q'}} \left( -E_d - (X_q - X_q') I_q \right)
**Excitation System Equations**
.. math::
\dot{E}_{\textup{fd}} &= \frac{1}{T_{E}} \left( -K_{E,i}E_{\textup{fd},i} + V_{R,i}\right)\\
\dot{R}_{f} &= \frac{1}{T_{F} } \left(- R_{f} + \frac{K_{F}}{T_{F}}E_{\textup{fd}} \right)\\
\dot{V}_{R} &= \frac{1}{T_{A}} \left(-V_{R} + K_{A}R_{f} - \frac{K_{A}K_{F}}{T_{F}}E_{\textup{fd}} + K_{A}(v_{\textup{ref}} - v_i) \right)
**Turbine-Governor System Equations**
.. math::
\dot{p}_{\textup{m}} &= \frac{1}{T_{ch}} \left( p_{\textup{sv}}- p_{\textup{m}} \right)\\
\dot{p}_{sv} &= \frac{1}{T_{sv}} \left( - \frac{\Delta \omega}{R_d} - p_{\textup{sv}} + p_{\textup{ref}} \right)
"""
def __init__(self, avr=None, governor=None, pss=None, shaft=None) -> None:
super().__init__(avr=avr, governor=governor, pss=pss, shaft=shaft)
self._type = "Synchronous_machine"
self._name = "Synchronous_machine_transient_model"
# States
self.ns += 2
self.states.extend(["e_dprim", "e_qprim"])
self.units.extend(["p.u.", "p.u."])
self.e_dprim = np.array([], dtype=float)
self.e_qprim = np.array([], dtype=float)
self._states_noise.update(
{
"e_dprim": 1,
"e_qprim": 1,
}
)
self._states_init_error.update({"e_dprim": 0.1, "e_qprim": 0.1})
self._x0.update(
{
"delta": 0.0,
"omega": 0.0,
"e_dprim": -0.4,
"e_qprim": 1,
"psv": 0.5,
"pm": 0.5,
"Efd": 2.5,
"Rf": 0.0,
"Vr": 2.5,
}
)
# Params
self._params.update(
{"x_dprim": 0.05, "x_qprim": 0.1, "T_dprim": 8.0, "T_qprim": 0.8}
)
self._descr.update(
{
"T_dprim": "d-axis transient time constant",
"T_qprim": "q-axis transient time constant",
"e_dprim": "d-axis voltage behind transient reactance",
"e_qprim": "q-axis voltage behind transient reactance",
"x_dprim": "d-axis transient reactance",
"x_qprim": "q-axis transient reactance",
}
)
# Parameters
self.x_dprim = np.array([], dtype=float)
self.x_qprim = np.array([], dtype=float)
self.T_dprim = np.array([], dtype=float)
self.T_qprim = np.array([], dtype=float)
self.properties.update(
{
"fgcall": True,
"finit": True,
"init_data": True,
"xy_index": True,
"save_data": True,
"qcall": True,
"gcall": True,
}
)
self._init_data()
[docs]
def electromagnetic(self, dae: Dae):
"""Two-axis (transient) electromagnetic model: stator currents, the
e'_d/e'_q dynamics, and the air-gap power. The swing equation and
orchestration live in the base (rotor / fgcall template)."""
i_d, i_q = self.input_current(dae)
# Air-gap power consumed by the rotor swing (Pe = E'_d I_d + E'_q I_q
# + (X'_q - X'_d) I_d I_q).
Pe = (
dae.x[self.e_dprim] * i_d
+ dae.x[self.e_qprim] * i_q
+ (self.x_qprim - self.x_dprim) * i_d * i_q
)
dae.f[self.e_qprim] = (
1
/ self.T_dprim
* (-dae.x[self.e_qprim] + self.var_sym(dae, "Efd") - (self.x_d - self.x_dprim) * i_d)
) # Eq
dae.f[self.e_dprim] = (
1 / self.T_qprim * (-dae.x[self.e_dprim] + (self.x_q - self.x_qprim) * i_q)
) # Ed
return i_d, i_q, Pe
[docs]
class SynchronousSubtransient(Synchronous):
r"""Subtransient Anderson Fouad SG with TGOV1 governor and IEEEDC1A AVR
The subtransient behavior of the synchronous generator is described by the following differential equations:
**Rotor Dynamics**
.. math::
\dot{\delta} &= 2 \pi f_n \Delta \omega \\
\Delta \dot{\omega} &= \frac{1}{2 H} \Big( P_m - E_{d}^{\prime\prime} I_d - E_{q}^{\prime\prime} I_q + (X_q^{\prime\prime} - X_d^{\prime\prime}) I_d I_q - D \Delta \omega - f (\Delta \omega + 1) \Big) \\
**Electromagnetic Equations**
.. math::
\dot{E}_q^{\prime} &= \frac{1}{T_{d}^{\prime}} \Big( -E_q + E_f + (X_d - X_d^{\prime}) I_d \Big) \\
\dot{E}_d^{\prime} &= \frac{1}{T_{q}^{\prime}} \Big( -E_d - (X_q - X_q^{\prime}) I_q \Big) \\
\dot{E}_{q}^{\prime\prime} &= \frac{1}{T_{d}^{\prime\prime}} \Big( E_q - E_{q}^{\prime\prime} + (X_d^{\prime} - X_d^{\prime\prime}) I_d \Big) \\
\dot{E}_{d}^{\prime\prime} &= \frac{1}{T_{q}^{\prime\prime}} \Big( E_d - E_{d}^{\prime\prime} - (X_q^{\prime} - X_q^{\prime\prime}) I_q \Big) \\
**Excitation System Equations**
.. math::
\dot{E}_{\textup{fd}} &= \frac{1}{T_{E}} \left( -K_{E,i}E_{\textup{fd},i} + V_{R,i}\right)\\
\dot{R}_{f} &= \frac{1}{T_{F} } \left(- R_{f} + \frac{K_{F}}{T_{F}}E_{\textup{fd}} \right)\\
\dot{V}_{R} &= \frac{1}{T_{A}} \left(-V_{R} + K_{A}R_{f} - \frac{K_{A}K_{F}}{T_{F}}E_{\textup{fd}} + K_{A}(v_{\textup{ref}} - v_i) \right)
**Turbine-Governor System Equations**
.. math::
\dot{p}_{\textup{m}} &= \frac{1}{T_{ch}} \left( p_{\textup{sv}}- p_{\textup{m}} \right)\\
\dot{p}_{sv} &= \frac{1}{T_{sv}} \left( - \frac{\Delta \omega}{R_d} - p_{\textup{sv}} + p_{\textup{ref}} \right)
"""
def __init__(self, avr=None, governor=None, pss=None, shaft=None) -> None:
super().__init__(avr=avr, governor=governor, pss=pss, shaft=shaft)
# private data
self._type = "Synchronous_machine"
self._name = "Synchronous_machine_subtransient_model"
# States
self.ns += 4
self.states.extend(["e_dprim", "e_qprim", "e_dsec", "e_qsec"])
self.units.extend(["p.u.", "p.u.", "p.u.", "p.u."])
self.e_dprim = np.array([], dtype=float)
self.e_qprim = np.array([], dtype=float)
self.e_dsec = np.array([], dtype=float)
self.e_qsec = np.array([], dtype=float)
self._states_noise.update(
{"e_dprim": 1, "e_qprim": 1, "e_dsec": 1, "e_qsec": 1}
)
self._states_init_error.update(
{"e_dprim": 0.1, "e_qprim": 0.1, "e_dsec": 0.1, "e_qsec": 0.1}
)
self._x0.update(
{
"delta": 0.1,
"omega": 0.0,
"e_dprim": 0.0,
"e_qprim": 1.0,
"psv": 0.5,
"pm": 0.5,
"Efd": 2.3,
"Rf": 0.0,
"Vr": 2.3,
"e_dsec": 0.0,
"e_qsec": 1.0,
}
)
# Params
self._params.update(
{
"x_dprim": 0.05,
"x_qprim": 0.1,
"T_dprim": 8.0,
"T_qprim": 0.8,
"x_dsec": 0.01,
"x_qsec": 0.01,
"T_dsec": 0.001,
"T_qsec": 0.001,
}
)
self._descr.update(
{
"T_dprim": "d-axis transient time constant",
"T_qprim": "q-axis transient time constant",
"x_dprim": "d-axis transient reactance",
"x_qprim": "q-axis transient reactance",
"e_dprim": "d-axis voltage behind transient reactance",
"e_qprim": "q-axis voltage behind transient reactance",
"e_dsec": "d-axis voltage behind subtransient reactance",
"e_qsec": "q-axis voltage behind subtransient reactance",
"T_dsec": "d-axis subtransient time constant",
"T_qsec": "q-axis subtransient time constant",
"x_dsec": "d-axis subtransient reactance",
"x_qsec": "q-axis subtransient reactance",
}
)
# Parameters
self.x_dprim = np.array([], dtype=float)
self.x_qprim = np.array([], dtype=float)
self.T_dprim = np.array([], dtype=float)
self.T_qprim = np.array([], dtype=float)
self.x_dsec = np.array([], dtype=float)
self.x_qsec = np.array([], dtype=float)
self.T_dsec = np.array([], dtype=float)
self.T_qsec = np.array([], dtype=float)
self.properties.update(
{
"fgcall": True,
"finit": True,
"init_data": True,
"xy_index": True,
"save_data": True,
"qcall": True,
}
)
self._init_data()
[docs]
def electromagnetic(self, dae: Dae):
"""Anderson-Fouad (subtransient) electromagnetic model: stator currents,
the e'/e'' dynamics, and the air-gap power. Swing/orchestration in base."""
i_d, i_q = self.input_current(dae)
# Air-gap power (subtransient emfs): Pe = E''_d I_d + E''_q I_q
# + (X''_q - X''_d) I_d I_q.
Pe = (
dae.x[self.e_dsec] * i_d
+ dae.x[self.e_qsec] * i_q
+ (self.x_qsec - self.x_dsec) * i_d * i_q
)
dae.f[self.e_qprim] = (
1
/ self.T_dprim
* (-dae.x[self.e_qprim] + self.var_sym(dae, "Efd") - (self.x_d - self.x_dprim) * i_d)
) # Eq
dae.f[self.e_dprim] = (
1 / self.T_qprim * (-dae.x[self.e_dprim] + (self.x_q - self.x_qprim) * i_q)
) # Ed
dae.f[self.e_qsec] = (
1
/ self.T_dsec
* (
dae.x[self.e_qprim]
- dae.x[self.e_qsec]
- (self.x_dprim - self.x_dsec) * i_d
)
)
dae.f[self.e_dsec] = (
1
/ self.T_qsec
* (
dae.x[self.e_dprim]
- dae.x[self.e_dsec]
+ (self.x_qprim - self.x_qsec) * i_q
)
)
return i_d, i_q, Pe
[docs]
class SynchronousSubtransientSP(Synchronous):
r"""Subtransient Sauer and Pai SG model with stator dynamics
with TGOV1 governor and IEEEDC1A AVR
The model includes the following equations for rotor dynamics, stator dynamics, and the excitation system:
**Rotor Dynamics**
.. math::
\dot{\delta} &= 2 \pi f_n \Delta \omega \\
\Delta \dot{\omega} &= \frac{1}{2 H} \Big( P_m - (\psi_d I_q - \psi_q I_d) - D \Delta \omega - f (\Delta \omega + 1) \Big)
**Electromagnetic Equations**
The stator dynamics include the following equations for the flux linkages in the d and q axes:
.. math::
\dot{E}_d' &= \frac{1}{T_q'} \Big( -E_d' + (X_q - X_q') (i_q - g_{q2} \Psi_{q2} - (1 - g_{q1}) i_q - g_{q2} E_d') \Big) \\
\dot{E}_q' &= \frac{1}{T_d'} \Big( -E_q' - (X_d - X_d') (i_d - g_{d2} \Psi_{d2} - (1 - g_{d1}) i_d + g_{d2} E_q') + E_f \Big)\\
\dot{\Psi}_{d2} &= \frac{1}{T_{d2}} \Big( -\Psi_{d2} + E_q' - (X_d' - X_l) i_d \Big) \\
\dot{\Psi}_{q2} &= \frac{1}{T_{q2}} \Big( -\Psi_{q2} - E_d' - (X_q' - X_l) i_q \Big)
**Flux Linkage Dynamics**
The following equations describe the stator flux linkage dynamics in the d and q axes:
.. math::
\dot{\Psi}_d &= 2 \pi f_n (R_s i_d + (1 + \Delta \omega) \Psi_q + v_d) \\
\dot{\Psi}_q &= 2 \pi f_n (R_s i_q - (1 + \Delta \omega) \Psi_d + v_q)
**Algebraic Equations**
The following algebraic equations govern the system:
.. math::
i_d &= \frac{1}{x_d''} \Big( -\psi_d + g_{d1} e_q' + (1 - g_{d1}) \psi_{d2} \Big)\\
i_q &= \frac{1}{x_q''} \Big( -\psi_q - g_{q1} e_d' + (1 - g_{q1}) \psi_{q2} \Big)\\
g_{d1} &= \frac{x_d'' - x_l}{x_d' - x_l}\\
g_{q1} &= \frac{x_q'' - x_l}{x_q' - x_l}\\
g_{d2} &= \frac{1 - g_{d1}}{x_d' - x_l}\\
g_{q2} &= \frac{1 - g_{q1}}{x_q' - x_l}
"""
def __init__(self, avr=None, governor=None, pss=None, shaft=None) -> None:
super().__init__(avr=avr, governor=governor, pss=pss, shaft=shaft)
# private data
self._type = "Synchronous_machine"
self._name = "Synchronous_machine_subtransient_model_Sauer_Pai"
# States
self.ns += 6
self.states.extend(["e_dprim", "e_qprim", "psid", "psiq", "psid2", "psiq2"])
self.units.extend(["p.u.", "p.u.", "p.u.", "p.u.", "p.u.", "p.u."])
self.e_dprim = np.array([], dtype=float)
self.e_qprim = np.array([], dtype=float)
self.psid = np.array([], dtype=float)
self.psiq = np.array([], dtype=float)
self.psid2 = np.array([], dtype=float)
self.psiq2 = np.array([], dtype=float)
self._states_noise.update(
{
"e_dprim": 1.0,
"e_qprim": 1.0,
"psid": 1.0,
"psiq": 1.0,
"psid2": 1.0,
"psiq2": 1.0,
}
)
self._states_init_error.update(
{
"e_dprim": 0.1,
"e_qprim": 0.1,
"psid": 0.1,
"psiq": 0.1,
"psid2": 0.1,
"psiq2": 0.1,
}
)
self._x0.update(
{
"delta": 0.5,
"omega": 0.0,
"e_dprim": 0.2,
"e_qprim": 1.0,
"psid": 1.0,
"psiq": -0.5,
"psid2": 1.0,
"psiq2": -0.5,
"psv": 0.5,
"pm": 0.5,
"Efd": 2.3,
"Rf": 0.0,
"Vr": 2.3,
}
)
# Params
self._params.update(
{
"gd1": 1.0,
"gq1": 1.0,
"gd2": 1.0,
"gq2": 1.0,
"x_l": 0.1,
"x_dprim": 0.05,
"x_qprim": 0.1,
"T_dprim": 8.0,
"T_qprim": 0.8,
"x_dsec": 0.01,
"x_qsec": 0.01,
"T_dsec": 0.001,
"T_qsec": 0.001,
}
)
self._descr.update(
{
"T_dprim": "d-axis transient time constant",
"T_qprim": "q-axis transient time constant",
"x_dprim": "d-axis transient reactance",
"x_qprim": "q-axis transient reactance",
"e_dprim": "d-axis voltage behind transient reactance",
"e_qprim": "q-axis voltage behind transient reactance",
"T_dsec": "d-axis subtransient time constant",
"T_qsec": "q-axis subtransient time constant",
"x_dsec": "d-axis subtransient reactance",
"x_qsec": "q-axis subtransient reactance",
"x_l": "leakage reactance",
"psid": "stator flux in d axis",
"psiq": "stator flux in q axis",
"psiq2": "subtransient stator flux in q axis",
"psid2": "subtransient stator flux in d axis",
}
)
# Parameters
self.x_l = np.array([], dtype=float)
self.gd1 = np.array([], dtype=float)
self.gq1 = np.array([], dtype=float)
self.gd2 = np.array([], dtype=float)
self.gq2 = np.array([], dtype=float)
self.x_dprim = np.array([], dtype=float)
self.x_qprim = np.array([], dtype=float)
self.T_dprim = np.array([], dtype=float)
self.T_qprim = np.array([], dtype=float)
self.x_dsec = np.array([], dtype=float)
self.x_qsec = np.array([], dtype=float)
self.T_dsec = np.array([], dtype=float)
self.T_qsec = np.array([], dtype=float)
self.properties.update(
{
"fgcall": True,
"finit": True,
"init_data": True,
"xy_index": True,
"save_data": True,
"qcall": True,
}
)
self._init_data()
[docs]
def sauer_pai(self, dae: Dae, i_d: ca.SX, i_q: ca.SX):
"""
Sauer and Pai model.
Parameters
----------
dae : DAE
i_d : casadi.SX
i_q : casadi.SX
Returns
-------
"""
vd = dae.y[self.vre] * np.sin(dae.x[self.delta]) + dae.y[self.vim] * -np.cos(
dae.x[self.delta]
)
vq = dae.y[self.vre] * np.cos(dae.x[self.delta]) + dae.y[self.vim] * np.sin(
dae.x[self.delta]
)
# Air-gap power/torque consumed by the rotor swing (Pe = psi_d I_q - psi_q I_d).
Pe = dae.x[self.psid] * i_q - dae.x[self.psiq] * i_d
dae.f[self.e_dprim] = (
1
/ self.T_qprim
* (
-dae.x[self.e_dprim]
+ (self.x_q - self.x_qprim)
* (
i_q
- self.gq2 * dae.x[self.psiq2]
- (1 - self.gq1) * i_q
- self.gq2 * dae.x[self.e_dprim]
)
)
)
dae.f[self.e_qprim] = (
1
/ self.T_dprim
* (
-dae.x[self.e_qprim]
- (self.x_d - self.x_dprim)
* (
i_d
- self.gd2 * dae.x[self.psid2]
- (1 - self.gd1) * i_d
+ self.gd2 * dae.x[self.e_qprim]
)
+ self.var_sym(dae, "Efd")
)
)
dae.f[self.psid2] = (
1
/ self.T_dsec
* (
-dae.x[self.psid2]
+ dae.x[self.e_qprim]
- (self.x_dprim - self.x_l) * i_d
)
)
dae.f[self.psiq2] = (
1
/ self.T_qsec
* (
-dae.x[self.psiq2]
- dae.x[self.e_dprim]
- (self.x_qprim - self.x_l) * i_q
)
)
omega_rh = dae.x[self.omega]
# # For the distributed mode we have assumtion of linearity
# mode = getattr(dae, "omega_mode", "nom")
# if mode == "dist":
# # Frequency along the internal stator reactances
# # Defined as average of rotor speed and bus frequency from FDF
# omega_rh = 0.5 * (dae.x[self.omega] + omega_ref_vec)
dae.f[self.psid] = (
2 * np.pi * dae.fn * (self.R_s * i_d + (omega_rh) * dae.x[self.psiq] + vd)
)
dae.f[self.psiq] = (
2 * np.pi * dae.fn * (self.R_s * i_q - (omega_rh) * dae.x[self.psid] + vq)
)
return Pe
[docs]
def _gd_coeffs(self) -> None:
"""Subtransient coupling coefficients used by both the stator-current
expressions and the flux ODEs."""
self.gd1 = (self.x_dsec - self.x_l) / (self.x_dprim - self.x_l)
self.gq1 = (self.x_qsec - self.x_l) / (self.x_qprim - self.x_l)
self.gd2 = (1 - self.gd1) / (self.x_dprim - self.x_l)
self.gq2 = (1 - self.gq1) / (self.x_qprim - self.x_l)
[docs]
def electromagnetic(self, dae: Dae):
"""Sauer-Pai (subtransient, stator dynamics) electromagnetic model.
The stator-current derivation is :meth:`input_current` (overridden by the
_DAE variant to expose i_d/i_q as private algebraics); the shared flux/
stator physics and the air-gap power live in :meth:`sauer_pai`."""
i_d, i_q = self.input_current(dae)
Pe = self.sauer_pai(dae, i_d, i_q)
return i_d, i_q, Pe
[docs]
class SynchronousSubtransientSP6(Synchronous):
r"""Subtransient Sauer and Pai SG 6th order model with neglected stator dynamics
(stator modeled with algebraic equations) and included TGOV1 governor and
IEEEDC1A AVR
The model includes the following equations:
**Rotor Dynamics**
.. math::
\dot{\delta} &= 2 \pi f_n \Delta \omega \\
\Delta \dot{\omega} &= \frac{1}{2 H} \Big( P_m - (\psi_d I_q - \psi_q I_d) - D \Delta \omega - f (\Delta \omega + 1) \Big)
**Electromagnetic Equations**
The stator dynamics include the following equations for the flux linkages in the d and q axes:
.. math::
\dot{E}_d' &= \frac{1}{T_q'} \Big( -E_d' + (X_q - X_q') (i_q - g_{q2} \Psi_{q2} - (1 - g_{q1}) i_q - g_{q2} E_d') \Big) \\
\dot{E}_q' &= \frac{1}{T_d'} \Big( -E_q' - (X_d - X_d') (i_d - g_{d2} \Psi_{d2} - (1 - g_{d1}) i_d + g_{d2} E_q') + E_f \Big)\\
\dot{\Psi}_{d2} &= \frac{1}{T_{d2}} \Big( -\Psi_{d2} + E_q' - (X_d' - X_l) i_d \Big) \\
\dot{\Psi}_{q2} &= \frac{1}{T_{q2}} \Big( -\Psi_{q2} - E_d' - (X_q' - X_l) i_q \Big)
**Flux Linkage Dynamics**
The following equations describe the stator flux linkage dynamics in the d and q axes:
**Algebraic Equations**
The following algebraic equations govern the system:
.. math::
0 &=-i_d +\frac{1}{x_d''} \Big( -\psi_d + g_{d1} e_q' + (1 - g_{d1}) \psi_{d2} \Big)\\
0&=-i_q +\frac{1}{x_q''} \Big( -\psi_q - g_{q1} e_d' + (1 - g_{q1}) \psi_{q2} \Big)\\
0&= R_s i_d + (1 + \Delta \omega) \Psi_q + v_d \\
0&= R_s i_q - (1 + \Delta \omega) \Psi_d + v_q \\
g_{d1} &= \frac{x_d'' - x_l}{x_d' - x_l}\\
g_{q1} &= \frac{x_q'' - x_l}{x_q' - x_l}\\
g_{d2} &= \frac{1 - g_{d1}}{x_d' - x_l}\\
g_{q2} &= \frac{1 - g_{q1}}{x_q' - x_l}
"""
def __init__(self, avr=None, governor=None, pss=None, shaft=None) -> None:
super().__init__(avr=avr, governor=governor, pss=pss, shaft=shaft)
# private data
self._type = "Synchronous_machine"
self._name = "Synchronous_machine_subtransient_model_Sauer_Pai_6th_order"
# States
self.ns += 4
self.states.extend(["e_dprim", "e_qprim", "psid2", "psiq2"])
self.units.extend(["p.u.", "p.u.", "p.u.", "p.u."])
self.e_dprim = np.array([], dtype=float)
self.e_qprim = np.array([], dtype=float)
self.psid2 = np.array([], dtype=float)
self.psiq2 = np.array([], dtype=float)
self._states_noise.update(
{
"e_dprim": 1.0,
"e_qprim": 1.0,
"psid2": 1.0,
"psiq2": 1.0,
}
)
self._states_init_error.update(
{
"e_dprim": 0.1,
"e_qprim": 0.1,
"psid2": 0.1,
"psiq2": 0.1,
}
)
# finit Newton seed. The original guess (omega=0.0, Efd==Vr==2.3) sits in
# a fragile basin: omega is ABSOLUTE p.u. speed (steady state = 1.0, not 0),
# and the equal Efd/Vr is the degenerate d-axis excitation/flux chain noted
# for SP_DAE below. macOS/Accelerate tolerated it, but Linux/OpenBLAS
# converged the per-device finit to a spurious root (delta~pi, omega~0,
# fluxes ~1e3) -> the eliminated SP6 reference diverged from SP6_DAE in CI.
# Seed in the correct basin: true synchronous speed, delta near its true
# value, and distinct Efd/Vr to break the symmetry (cf. SP_DAE remedy).
self._x0.update(
{
"delta": 0.3,
"omega": 1.0,
"e_dprim": 0.2,
"e_qprim": 1.0,
"psid2": 1.0,
"psiq2": -0.5,
"psv": 0.5,
"pm": 0.5,
"Efd": 2.3,
"Rf": 0.0,
"Vr": 2.0,
}
)
# Params
self._params.update(
{
"gd1": 1.0,
"gq1": 1.0,
"gd2": 1.0,
"gq2": 1.0,
"x_l": 0.1,
"x_dprim": 0.05,
"x_qprim": 0.1,
"T_dprim": 8.0,
"T_qprim": 0.8,
"x_dsec": 0.01,
"x_qsec": 0.01,
"T_dsec": 0.001,
"T_qsec": 0.001,
}
)
self._descr.update(
{
"T_dprim": "d-axis transient time constant",
"T_qprim": "q-axis transient time constant",
"x_dprim": "d-axis transient reactance",
"x_qprim": "q-axis transient reactance",
"e_dprim": "d-axis voltage behind transient reactance",
"e_qprim": "q-axis voltage behind transient reactance",
"T_dsec": "d-axis subtransient time constant",
"T_qsec": "q-axis subtransient time constant",
"x_dsec": "d-axis subtransient reactance",
"x_qsec": "q-axis subtransient reactance",
"x_l": "leakage reactance",
"psid2": "subtransient flux in d axis",
"psiq2": "subtransient flux in q axis",
}
)
# Parameters
self.x_l = np.array([], dtype=float)
self.gd1 = np.array([], dtype=float)
self.gq1 = np.array([], dtype=float)
self.gd2 = np.array([], dtype=float)
self.gq2 = np.array([], dtype=float)
self.x_dprim = np.array([], dtype=float)
self.x_qprim = np.array([], dtype=float)
self.T_dprim = np.array([], dtype=float)
self.T_qprim = np.array([], dtype=float)
self.x_dsec = np.array([], dtype=float)
self.x_qsec = np.array([], dtype=float)
self.T_dsec = np.array([], dtype=float)
self.T_qsec = np.array([], dtype=float)
self.properties.update(
{
"fgcall": True,
"finit": True,
"init_data": True,
"xy_index": True,
"save_data": True,
"qcall": True,
}
)
self._init_data()
[docs]
def sauer_pai_6(self, dae: Dae, i_d: ca.SX, i_q: ca.SX, psid: ca.SX, psiq: ca.SX):
r"""
Sauer and Pai model.
Parameters
----------
dae : DAE
i_d : casadi.SX
i_q : casadi.SX
psid : casadi.SX
psiq : casadi.SX
Returns
-------
Args:
psid ():
psiq ():
"""
# Air-gap power/torque consumed by the rotor swing (Pe = psi_d I_q - psi_q I_d).
Pe = psid * i_q - psiq * i_d
dae.f[self.e_dprim] = (
1
/ self.T_qprim
* (
-dae.x[self.e_dprim]
+ (self.x_q - self.x_qprim)
* (
i_q
- self.gq2 * dae.x[self.psiq2]
- (1 - self.gq1) * i_q
- self.gq2 * dae.x[self.e_dprim]
)
)
)
dae.f[self.e_qprim] = (
1
/ self.T_dprim
* (
-dae.x[self.e_qprim]
- (self.x_d - self.x_dprim)
* (
i_d
- self.gd2 * dae.x[self.psid2]
- (1 - self.gd1) * i_d
+ self.gd2 * dae.x[self.e_qprim]
)
+ self.var_sym(dae, "Efd")
)
)
dae.f[self.psid2] = (
1
/ self.T_dsec
* (
-dae.x[self.psid2]
+ dae.x[self.e_qprim]
- (self.x_dprim - self.x_l) * i_d
)
)
dae.f[self.psiq2] = (
1
/ self.T_qsec
* (
-dae.x[self.psiq2]
- dae.x[self.e_dprim]
- (self.x_qprim - self.x_l) * i_q
)
)
return Pe
[docs]
def stator_solve(self, dae: Dae):
"""Eliminate the four algebraic stator unknowns (i_d, i_q, psi_d, psi_q)
by a symbolic 4x4 solve per machine. Returns (i_d, i_q, psid, psiq).
The _DAE variant overrides this to expose them as private algebraics."""
self.gd1 = (self.x_dsec - self.x_l) / (self.x_dprim - self.x_l)
self.gq1 = (self.x_qsec - self.x_l) / (self.x_qprim - self.x_l)
self.gd2 = (1 - self.gd1) / (self.x_dprim - self.x_l)
self.gq2 = (1 - self.gq1) / (self.x_qprim - self.x_l)
i_d = ca.SX.sym("i_d", self.n)
i_q = ca.SX.sym("i_q", self.n)
psid = ca.SX.sym("psid", self.n)
psiq = ca.SX.sym("psiq", self.n)
for i in range(self.n):
# Symbolic variables for unknowns
algs = ca.SX.sym("algs", 4) # symbolic unknowns
# Define vd and vq in terms of symbolic variables (using symbolic dae.x and dae.y)
vd = dae.y[self.vre][i] * ca.sin(dae.x[self.delta][i]) + dae.y[self.vim][
i
] * -ca.cos(dae.x[self.delta][i])
vq = dae.y[self.vre][i] * ca.cos(dae.x[self.delta][i]) + dae.y[self.vim][
i
] * ca.sin(dae.x[self.delta][i])
# Define g as a symbolic vector in terms of algs, dae.x, and dae.y
g = ca.SX(4, 1)
g[0] = -algs[0] + (1 / self.x_dsec[i]) * (
-algs[2]
+ self.gd1[i] * dae.x[self.e_qprim][i]
+ (1 - self.gd1[i]) * dae.x[self.psid2][i]
)
g[1] = -algs[1] + (1 / self.x_qsec[i]) * (
-algs[3]
- self.gq1[i] * dae.x[self.e_dprim][i]
+ (1 - self.gq1[i]) * dae.x[self.psiq2][i]
)
# Algebraic stator (neglected stator transients): 0 = R_s*i + ...
# These rows previously carried a 2*pi*f_n factor inherited from the
# stator *dynamic* equations (where it is the per-unit-time
# coefficient d/dt = omega_base*(...)). For an algebraic constraint
# (= 0) that factor is mathematically redundant -- it does NOT change
# the solved currents/fluxes (ca.solve scales J and g together) --
# but it scaled this row by ~314 and, propagated through the symbolic
# elimination, blew the per-device finit Jacobian condition number up
# to ~1e11 (smax ~3.4e7), so the init Newton/LM could not converge.
# Dropping it makes SP6 well-posed (cond ~7e3, residual ~1e-14),
# matching the stator-dynamics SP model. See
# CHANGELOG_algebraic_equations.md. NOTE: also uses omega, not 1+omega.
g[2] = self.R_s[i] * algs[0] + (dae.x[self.omega][i]) * algs[3] + vd
g[3] = self.R_s[i] * algs[1] - (dae.x[self.omega][i]) * algs[2] + vq
# Calculate the Jacobian of g with respect to algs
J = ca.jacobian(g, algs)
g_eval = ca.substitute(g, algs, ca.DM.zeros(4))
# Solve J * x = -g symbolically
sol = ca.solve(J, -g_eval)
# Assign solutions to variables
i_d[i] = sol[0]
i_q[i] = sol[1]
psid[i] = sol[2]
psiq[i] = sol[3]
return i_d, i_q, psid, psiq
[docs]
def electromagnetic(self, dae: Dae):
"""Sauer-Pai 6th-order (algebraic stator) electromagnetic model. The four
stator unknowns come from :meth:`stator_solve` (overridden by the _DAE
variant to expose them as private algebraics); the shared flux physics and
the air-gap power live in :meth:`sauer_pai_6`."""
i_d, i_q, psid, psiq = self.stator_solve(dae)
Pe = self.sauer_pai_6(dae, i_d, i_q, psid, psiq)
return i_d, i_q, Pe
[docs]
class SynchronousSubtransientSP6DAE(SynchronousSubtransientSP6):
r"""Subtransient Sauer-Pai 6th-order machine expressed as an explicit DAE.
Identical physics to :class:`SynchronousSubtransientSP6`, but instead of
symbolically eliminating the four stator unknowns (``i_d, i_q, psi_d,
psi_q``) with a per-instance ``ca.solve``, the four stator equations are
declared as device-private algebraic variables/equations (``_algebs_int``)
and handed to the DAE solver. This validates the private-algebraic mechanism:
it must reproduce :class:`SynchronousSubtransientSP6` to integrator tolerance.
See ``docs/algebraic_equations_design.md`` and
``CHANGELOG_algebraic_equations.md``.
The stator block is *linear* in the unknowns, so the parent's ``ca.solve``
is the exact solution of ``g = 0`` -- the DAE integrator converges to the
same values, only via a numerical Newton solve per step instead of a
closed-form expression built once.
.. note::
:class:`SynchronousSubtransientSP6` initially appeared ill-posed at init,
but the root cause was a spurious ``2*pi*f_n`` factor on its algebraic
stator equations (a leftover from the stator *dynamic* equations) that
inflated the per-device ``finit`` Jacobian condition number to ~1e11.
With that factor removed (see the parent's ``fgcall``), SP6 is well-posed
and this DAE form reproduces it to ~9e-7 (see
``pydynamicestimator/tests/test_algebraic_parity.py``).
"""
def __init__(self, avr=None, governor=None, pss=None, shaft=None) -> None:
super().__init__(avr=avr, governor=governor, pss=pss, shaft=shaft)
self._name = "Synchronous_machine_subtransient_model_Sauer_Pai_6th_order_DAE"
# Stator currents/fluxes become private algebraic variables instead of
# being eliminated. id_alg/iq_alg <-> algs[0]/algs[1] and
# psid_alg/psiq_alg <-> algs[2]/algs[3] in the parent's fgcall.
# Extend/update (not assign) so these compose with any AVR-declared
# privates already registered by Synchronous.__init__.
machine_privs = ["id_alg", "iq_alg", "psid_alg", "psiq_alg"]
self._algebs_int.extend(machine_privs)
self._algebs_int_x0.update(
{
"id_alg": 0.5,
"iq_alg": 0.5,
"psid_alg": 1.0,
"psiq_alg": -0.5,
}
)
self._algebs_int_noise.update({name: 1.0 for name in machine_privs})
self._algebs_int_units.update({name: "p.u." for name in machine_privs})
for name in machine_privs:
setattr(self, name, np.array([], dtype=float))
self._descr.update(
{
"id_alg": "stator d-axis current (algebraic)",
"iq_alg": "stator q-axis current (algebraic)",
"psid_alg": "stator d-axis flux (algebraic)",
"psiq_alg": "stator q-axis flux (algebraic)",
}
)
[docs]
def stator_solve(self, dae: Dae):
"""The four stator unknowns are genuine device-private algebraic variables
in dae.y; write their defining equations and return them (the integrator
solves the 4x4 block instead of ca.solve). Inherits SP6's
:meth:`electromagnetic` and the base orchestration."""
self.gd1 = (self.x_dsec - self.x_l) / (self.x_dprim - self.x_l)
self.gq1 = (self.x_qsec - self.x_l) / (self.x_qprim - self.x_l)
self.gd2 = (1 - self.gd1) / (self.x_dprim - self.x_l)
self.gq2 = (1 - self.gq1) / (self.x_qprim - self.x_l)
i_d = dae.y[self.id_alg]
i_q = dae.y[self.iq_alg]
psid = dae.y[self.psid_alg]
psiq = dae.y[self.psiq_alg]
vd = dae.y[self.vre] * np.sin(dae.x[self.delta]) + dae.y[self.vim] * -np.cos(
dae.x[self.delta]
)
vq = dae.y[self.vre] * np.cos(dae.x[self.delta]) + dae.y[self.vim] * np.sin(
dae.x[self.delta]
)
# The four stator defining equations (cf. SynchronousSubtransientSP6
# g[0..3]), written into each private's own dae.g slot (the 4x4 block is
# solved by the integrator instead of ca.solve). NOTE: the flux equations
# use dae.x[omega], NOT (1 + omega). The 2*pi*f_n factor that the stator
# *dynamic* equations carry is dropped here: these are algebraic
# constraints (= 0), so the factor is mathematically redundant and only
# inflates the Jacobian condition number (see the SP6 fix and
# CHANGELOG_algebraic_equations.md).
dae.g[self.id_alg] = -i_d + (1 / self.x_dsec) * (
-psid + self.gd1 * dae.x[self.e_qprim] + (1 - self.gd1) * dae.x[self.psid2]
)
dae.g[self.iq_alg] = -i_q + (1 / self.x_qsec) * (
-psiq - self.gq1 * dae.x[self.e_dprim] + (1 - self.gq1) * dae.x[self.psiq2]
)
dae.g[self.psid_alg] = self.R_s * i_d + dae.x[self.omega] * psiq + vd
dae.g[self.psiq_alg] = self.R_s * i_q - dae.x[self.omega] * psid + vq
return i_d, i_q, psid, psiq
[docs]
class SynchronousSubtransientSP_DAE(SynchronousSubtransientSP):
r"""Subtransient Sauer-Pai machine (with stator dynamics) with the stator
currents exposed as explicit device-private algebraic variables.
Identical physics to :class:`SynchronousSubtransientSP`, but the stator
currents ``i_d, i_q`` -- normally inlined explicit expressions of the flux
states -- are declared as private algebraic variables (``_algebs_int``) with
their defining equations ``0 = -i + <expr>`` handed to the DAE solver. Must
reproduce :class:`SynchronousSubtransientSP` to integrator tolerance.
This is the working parity vehicle for the private-algebraic mechanism: the
parent model initialises robustly and the defining equations are linear with
a trivially non-singular Jacobian (index-1). See
``docs/algebraic_equations_design.md`` and
``CHANGELOG_algebraic_equations.md``.
"""
def __init__(self, avr=None, governor=None, pss=None, shaft=None) -> None:
super().__init__(avr=avr, governor=governor, pss=pss, shaft=shaft)
self._name = "Synchronous_machine_subtransient_model_Sauer_Pai_DAE"
# Extend/update (not assign) so these compose with any AVR-declared
# privates already registered by Synchronous.__init__.
machine_privs = ["id_alg", "iq_alg"]
self._algebs_int.extend(machine_privs)
self._algebs_int_x0.update({"id_alg": 0.0, "iq_alg": 0.6})
self._algebs_int_noise.update({name: 1.0 for name in machine_privs})
self._algebs_int_units.update({name: "p.u." for name in machine_privs})
# The stock SP _x0 sits near a singular point of the per-device finit
# Jacobian (the d-axis excitation/flux chain Efd-Vr-e_qprim-psid2-psid is
# degenerate there). The stock model tolerates it under plain Newton, but
# the larger DAE init system diverges. Break the symmetry with distinct
# guesses in the same solution basin.
self._x0.update(
{
"delta": 0.3,
"omega": 0.0,
"e_dprim": 0.25,
"e_qprim": 1.1,
"psid": 0.95,
"psiq": -0.45,
"psid2": 0.9,
"psiq2": -0.4,
"psv": 0.5,
"pm": 0.5,
"Efd": 2.5,
"Rf": 0.05,
"Vr": 2.4,
}
)
for name in machine_privs:
setattr(self, name, np.array([], dtype=float))
self._descr.update(
{
"id_alg": "stator d-axis current (algebraic)",
"iq_alg": "stator q-axis current (algebraic)",
}
)