# Created: 2026-06-04
# (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.
"""Pluggable rotor-shaft strategy test (multi-mass torsional shaft).
The rotor mechanics are a pluggable strategy. The default :class:`SingleMass`
reproduces the classic single rigid-mass swing equation (its byte-equivalence is
gated by the unchanged positional baselines in
test_recur_sim/test_recur_est/test_run_ideal). A multi-mass (torsional) shaft
keeps the generator mass on the canonical ``delta``/``omega`` names -- so the
electromagnetic model, the network injection, the governor, the PSS and the
initialisation are untouched -- and adds further rotor masses
``delta_<name>``/``omega_<name>``.
This validates :class:`Shaft4Mass` (HP--IP--LP--GEN) and :class:`Shaft5Mass`
(HP--IP--LP--GEN--EXC):
1. Structural: each extra mass adds exactly two differential states and no
private algebraics; the generator mass stays 'delta'/'omega' as states 0/1.
2. Steady state: every mass is synchronous (omega_i -> 1 p.u.) and the static
shaft twist satisfies the torque balance -- each shaft section carries exactly
the cumulative mechanical power of the turbine masses upstream of it.
3. The disturbance excites the torsional dynamics (the inter-mass twist
oscillates), i.e. the elastic shaft is genuinely modelled.
See docs/algebraic_equations_design.md (shaft strategy).
"""
import numpy as np
from pydynamicestimator.config import config
from pydynamicestimator.run import run
from pydynamicestimator.devices.shaft import Shaft4Mass, Shaft5Mass
from pydynamicestimator.devices.synchronous import SynchronousSubtransientSP
from pydynamicestimator.tests.test_avr_algebraic import _COMMON, _machine
_T_DIST = 3.0 # disturbance time in the shaft_4mass fixture (sim_dist.txt)
[docs]
def test_shaft_structural():
"""A multi-mass shaft adds exactly two differential states per extra mass and
no private algebraics; the generator mass stays 'delta'/'omega' as states 0,1."""
single = SynchronousSubtransientSP()
assert single.states[:2] == ["delta", "omega"]
m4 = SynchronousSubtransientSP(shaft=Shaft4Mass())
# gen stays first; HP/IP/LP add delta_*/omega_* (3 masses * 2 states = +6).
assert m4.states[:2] == ["delta", "omega"]
assert m4.ns == single.ns + 6
for s in ("delta_hp", "omega_hp", "delta_ip", "omega_ip", "delta_lp", "omega_lp"):
assert s in m4.states
# The shaft declares per-mass / per-section / per-turbine parameters.
for p in ("H_hp", "H_ip", "H_lp", "K_hp_ip", "K_ip_lp", "K_lp_gen",
"F_hp", "F_ip", "F_lp"):
assert p in m4._params
m5 = SynchronousSubtransientSP(shaft=Shaft5Mass())
assert m5.ns == single.ns + 8 # +1 EXC mass over the 4-mass shaft
assert "delta_exc" in m5.states and "K_gen_exc" in m5._params
[docs]
def test_shaft_steady_state_torque_balance():
"""At the pre-disturbance steady state every mass is synchronous and the
static shaft twist carries the cumulative upstream mechanical power: a shaft
section between masses a (upstream) and b transmits sum(F_upstream)*pm, so
K_ab*(delta_a - delta_b) == sum(F_upstream)*pm."""
_, sim = run(config.updated(testsystemfile="shaft_4mass", **_COMMON))
m = _machine()
assert sim.n_priv == 0 # the shaft introduces only differential states
x = np.array(sim.x_full)
assert np.isfinite(x).all(), "non-finite trajectory (init or integration diverged)"
dt = (_COMMON["T_end"] - _COMMON["T_start"]) / (x.shape[1] - 1)
i_pre = int(round((_T_DIST - 0.5) / dt)) # a settled column before the disturbance
# Every mass synchronous at steady state (omega is ABSOLUTE p.u. speed).
for s in ("omega", "omega_hp", "omega_ip", "omega_lp"):
assert abs(x[getattr(m, s)[0]][i_pre] - 1.0) < 1e-6, f"{s} not synchronous"
pm = x[m.pm[0]][i_pre]
d_gen = x[m.delta[0]][i_pre]
d_hp = x[m.delta_hp[0]][i_pre]
d_ip = x[m.delta_ip[0]][i_pre]
d_lp = x[m.delta_lp[0]][i_pre]
K_hp_ip, K_ip_lp, K_lp_gen = m.K_hp_ip[0], m.K_ip_lp[0], m.K_lp_gen[0]
F_hp, F_ip, F_lp = m.F_hp[0], m.F_ip[0], m.F_lp[0]
# HP--IP carries F_hp*pm; IP--LP carries (F_hp+F_ip)*pm; LP--GEN carries pm.
assert abs(K_hp_ip * (d_hp - d_ip) - F_hp * pm) < 1e-4
assert abs(K_ip_lp * (d_ip - d_lp) - (F_hp + F_ip) * pm) < 1e-4
assert abs(K_lp_gen * (d_lp - d_gen) - (F_hp + F_ip + F_lp) * pm) < 1e-4
# Sanity: the turbine masses lead the generator (they push torque downstream).
assert d_hp > d_ip > d_lp > d_gen
[docs]
def test_shaft_torsional_excitation():
"""The disturbance excites the elastic shaft: the inter-mass twist is flat
before the disturbance (rigid-body steady state) and oscillates afterwards."""
_, sim = run(config.updated(testsystemfile="shaft_4mass", **_COMMON))
m = _machine()
x = np.array(sim.x_full)
twist = x[m.delta_hp[0]] - x[m.delta[0]] # HP-to-GEN shaft twist
dt = (_COMMON["T_end"] - _COMMON["T_start"]) / (x.shape[1] - 1)
i_dist = int(round(_T_DIST / dt))
pre = twist[: i_dist - 1]
post = twist[i_dist:]
assert np.std(pre) < 1e-8, "shaft twist not steady before the disturbance"
assert (post.max() - post.min()) > 1e-6, (
"no torsional oscillation excited -- the shaft is behaving rigidly"
)