diff --git a/CHANGELOG.md b/CHANGELOG.md index 5c7d6bcff..ded783307 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -52,6 +52,7 @@ - Added component model developer checklist to a README file. - Added `IEEEST` Stabilizer Model - Added `SEXS-PTI` Exciter Model +- Added `ESDC1A` Exciter Model - Added `GENSAL` Machine Model - Added 200 Bus Synthetic Illinois Case - Added node objects to `PowerElectronics` module & updated all examples to make use of them. diff --git a/GridKit/Model/PhasorDynamics/ComponentLibrary.hpp b/GridKit/Model/PhasorDynamics/ComponentLibrary.hpp index 6ac2375f3..0e110fd67 100644 --- a/GridKit/Model/PhasorDynamics/ComponentLibrary.hpp +++ b/GridKit/Model/PhasorDynamics/ComponentLibrary.hpp @@ -6,6 +6,7 @@ #include #include #include +#include #include #include #include diff --git a/GridKit/Model/PhasorDynamics/Exciter/CMakeLists.txt b/GridKit/Model/PhasorDynamics/Exciter/CMakeLists.txt index 1120e20c2..7e5e4b9e3 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/CMakeLists.txt +++ b/GridKit/Model/PhasorDynamics/Exciter/CMakeLists.txt @@ -3,5 +3,6 @@ # - Luke Lowery # ]] +add_subdirectory(ESDC1A) add_subdirectory(IEEET1) add_subdirectory(SEXS-PTI) diff --git a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/CMakeLists.txt b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/CMakeLists.txt new file mode 100644 index 000000000..599f7fcb2 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/CMakeLists.txt @@ -0,0 +1,54 @@ +# [[ +# Author(s): +# - Luke Lowery +# ]] + +set(_install_headers Esdc1a.hpp Esdc1aData.hpp) + +if(GRIDKIT_ENABLE_ENZYME) + gridkit_add_library( + phasor_dynamics_exciter_esdc1a + SOURCES Esdc1aEnzyme.cpp + HEADERS ${_install_headers} + INCLUDE_DIRECTORIES PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include + LINK_LIBRARIES + PUBLIC + GridKit::phasor_dynamics_core + PUBLIC + GridKit::phasor_dynamics_signal + PRIVATE + ClangEnzymeFlags + COMPILE_OPTIONS + PRIVATE + -mllvm + -enzyme-auto-sparsity=1 + -fno-math-errno) +else() + gridkit_add_library( + phasor_dynamics_exciter_esdc1a + SOURCES Esdc1a.cpp + HEADERS ${_install_headers} + INCLUDE_DIRECTORIES PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include + LINK_LIBRARIES + PUBLIC + GridKit::phasor_dynamics_core + PUBLIC + GridKit::phasor_dynamics_signal) +endif() + +gridkit_add_library( + phasor_dynamics_exciter_esdc1a_dependency_tracking + SOURCES Esdc1aDependencyTracking.cpp + INCLUDE_DIRECTORIES PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include + LINK_LIBRARIES + PUBLIC + GridKit::phasor_dynamics_core + PUBLIC + GridKit::phasor_dynamics_signal_dependency_tracking) + +target_link_libraries( + phasor_dynamics_components + INTERFACE GridKit::phasor_dynamics_exciter_esdc1a) +target_link_libraries( + phasor_dynamics_components_dependency_tracking + INTERFACE GridKit::phasor_dynamics_exciter_esdc1a_dependency_tracking) diff --git a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1a.cpp b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1a.cpp new file mode 100644 index 000000000..0dfe1f692 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1a.cpp @@ -0,0 +1,30 @@ +/** + * @file Esdc1a.cpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Non-Enzyme instantiation for the ESDC1A exciter model. + */ + +#include "Esdc1aImpl.hpp" + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Exciter + { + /** + * @brief Report that a separate Jacobian is unavailable in the plain build. + */ + template + int Esdc1a::evaluateJacobian() + { + Log::misc() << "Evaluate Jacobian for Esdc1a..." << std::endl; + Log::misc() << "Jacobian evaluation is not implemented!" << std::endl; + return 0; + } + + template class Esdc1a; + template class Esdc1a; + } // namespace Exciter + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1a.hpp b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1a.hpp new file mode 100644 index 000000000..dd42f55ac --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1a.hpp @@ -0,0 +1,186 @@ +/** + * @file Esdc1a.hpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Declaration of the ESDC1A exciter model. + */ + +#pragma once + +#include +#include +#include + +#include +#include +#include +#include + +namespace GridKit +{ + namespace PhasorDynamics + { + template + class BusBase; + + template + class SignalNode; + + namespace Exciter + { + /// Internal variables of an `Esdc1a`. + enum class Esdc1aInternalVariables : size_t + { + EFDP, ///< \f$E_{\mathrm{fd}}'\f$ Differential exciter field-voltage state [p.u.] + VC, ///< \f$V_C\f$ Differential filtered terminal-voltage magnitude [p.u.] + VR, ///< \f$V_R\f$ Differential voltage-regulator output [p.u.] + VF, ///< \f$V_F\f$ Differential stabilizing feedback state [p.u.] + XLL, ///< \f$x_{\mathrm{LL}}\f$ Differential input lead-lag denominator state [p.u.] + EV, ///< \f$e_V\f$ Algebraic voltage-error summing output [p.u.] + VLL, ///< \f$V_{\mathrm{LL}}\f$ Algebraic input lead-lag output [p.u.] + VHV, ///< \f$V_{\mathrm{HV}}\f$ Algebraic high-value gate output [p.u.] + SE, ///< \f$S_E\f$ Algebraic exciter saturation coefficient [p.u.] + VFE, ///< \f$V_{\mathrm{FE}}\f$ Algebraic exciter feedback drive [p.u.] + EFD, ///< \f$E_{\mathrm{fd}}\f$ Algebraic field-voltage output [p.u.] + MAXIMUM, ///< Number of ESDC1A internal variables + }; + + /// External signal variables read or initialized by an `Esdc1a`. + enum class Esdc1aExternalVariables : size_t + { + OMEGA, ///< \f$\omega\f$ Known machine speed deviation [p.u.] + VREF, ///< \f$V_{\mathrm{ref}}\f$ Unknown voltage-control reference [p.u.] + VS, ///< \f$V_S\f$ Known stabilizer input signal [p.u.] + VUEL, ///< \f$V_{\mathrm{UEL}}\f$ Known under-excitation limiter input [p.u.] + MAXIMUM, ///< Number of ESDC1A external signal variables + }; + + /** + * @brief IEEE DC1A excitation-system model (ESDC1A). + * + * @tparam scalar_type Plain real or differentiable scalar type. + * @tparam index_type Integer index type. + */ + template + class Esdc1a : public Component + { + using Component::abs_tol_; + using Component::allocated_; + using Component::alpha_; + using Component::f_; + using Component::gridkit_component_id_; + using Component::J_cols_buffer_; + using Component::J_rows_buffer_; + using Component::J_vals_buffer_; + using Component::nnz_; + using Component::residual_indices_; + using Component::size_; + using Component::tag_; + using Component::variable_indices_; + using Component::wb_; + using Component::y_; + using Component::yp_; + + public: + using ScalarT = scalar_type; + using IdxT = index_type; + using RealT = typename Component::RealT; + using BusT = BusBase; + using SignalT = SignalNode; + using ModelDataT = Esdc1aData; + using MonitorT = Model::VariableMonitor; + using InternalVariablesT = Esdc1aInternalVariables; + using ExternalVariablesT = Esdc1aExternalVariables; + + Esdc1a(BusT* bus); + Esdc1a(BusT* bus, const ModelDataT& data); + ~Esdc1a(); + + int setGridKitComponentID(IdxT component_id) override final; + int allocate() override final; + int verify() const override final; + int initialize() override final; + int tagDifferentiable() override final; + int setAbsoluteTolerance(RealT rel_tol) override final; + int evaluateResidual() override final; + int evaluateJacobian() override final; + + auto getSignals() + -> ComponentSignals& + { + return signals_; + } + + const Model::VariableMonitorBase* getMonitor() const override; + + __attribute__((always_inline)) inline int evaluateInternalResidual( + const ScalarT* y, + const ScalarT* yp, + const ScalarT* wb, + const ScalarT* ws, + ScalarT* f); + + private: + void initializeParameters(const ModelDataT& data); + void initializeMonitor(); + void setDerivedParameters(); + + static __attribute__((always_inline)) inline ScalarT awmin( + ScalarT x, + ScalarT f, + RealT xmin); + + /// Recover the input that the smooth CommonMath ramp maps to a + /// requested strictly positive output. + RealT inverseRamp(RealT ramp_output) const; + + ScalarT& Vr(); + ScalarT& Vi(); + + static constexpr RealT TIME_CONSTANT_MINIMUM = static_cast(1.0e-3); + + BusT* bus_{nullptr}; + + RealT Tr_{ZERO}; + RealT Ka_{static_cast(40.0)}; + RealT Ta_{static_cast(0.1)}; + RealT Tb_{ZERO}; + RealT Tc_{ZERO}; + RealT Vrmax_{ONE}; + RealT Vrmin_{static_cast(-1.0)}; + RealT Ke_{static_cast(0.1)}; + RealT Te_{static_cast(0.5)}; + RealT Kf_{static_cast(0.05)}; + RealT Tf1_{static_cast(0.7)}; + bool Spdmlt_{false}; + RealT E1_{static_cast(2.8)}; + RealT Se1_{static_cast(0.08)}; + RealT E2_{static_cast(3.7)}; + RealT Se2_{static_cast(0.33)}; + IdxT UEL_{0}; + bool exclim_{true}; + RealT spd_on_{0}; + RealT uel_on_{0}; + RealT lim_on_{1}; + RealT SA_{0}; + RealT SB_{0}; + + IdxT parameter_error_count_{0}; + + ScalarT omega_set_{0}; + ScalarT vref_set_{0}; + ScalarT vs_set_{0}; + ScalarT vuel_set_{0}; + + ComponentSignals + signals_; + std::unique_ptr monitor_; + + std::vector ws_; + std::vector ws_indices_; + }; + } // namespace Exciter + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aData.hpp b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aData.hpp new file mode 100644 index 000000000..0dfc3bc42 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aData.hpp @@ -0,0 +1,104 @@ +/** + * @file Esdc1aData.hpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Modeling data for the ESDC1A exciter model. + */ + +#pragma once + +#include + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Exciter + { + /// Parameter keys for the ESDC1A exciter model. Every parameter is + /// optional and retains its documented default when omitted. + enum class Esdc1aParameters + { + Tr, ///< \f$T_R\f$ Voltage transducer time constant [sec] + Ka, ///< \f$K_A\f$ Voltage-regulator gain [p.u.] + Ta, ///< \f$T_A\f$ Voltage-regulator time constant [sec] + Tb, ///< \f$T_B\f$ Input lead-lag denominator time constant [sec] + Tc, ///< \f$T_C\f$ Input lead-lag numerator time constant [sec] + Vrmax, ///< \f$V_R^{\max}\f$ Maximum voltage-regulator output [p.u.] + Vrmin, ///< \f$V_R^{\min}\f$ Minimum voltage-regulator output [p.u.] + Ke, ///< \f$K_E\f$ Exciter constant [p.u.] + Te, ///< \f$T_E\f$ Exciter time constant [sec] + Kf, ///< \f$K_F\f$ Stabilizing feedback gain [p.u.] + Tf1, ///< \f$T_{F1}\f$ Stabilizing feedback time constant [sec] + Spdmlt, ///< \f$s_{\mathrm{spd}}\f$ Field-voltage speed-multiplier flag [boolean] + E1, ///< \f$E_1\f$ First saturation voltage point [p.u.] + Se1, ///< \f$S_E(E_1)\f$ Saturation coefficient at \f$E_1\f$ [p.u.] + E2, ///< \f$E_2\f$ Second saturation voltage point [p.u.] + Se2, ///< \f$S_E(E_2)\f$ Saturation coefficient at \f$E_2\f$ [p.u.] + UEL, ///< \f$I_{\mathrm{UEL}}\f$ UEL input-routing selector [integer] + exclim ///< \f$s_{\mathrm{lim}}\f$ Exciter field-voltage-state lower-limit flag [boolean] + }; + + /// Buses for the ESDC1A exciter model. + enum class Esdc1aBuses : size_t + { + bus, ///< \f$V_{\mathrm{r}},V_{\mathrm{i}}\f$ Required Known terminal-bus voltage [p.u.] + SIZE ///< Number of ESDC1A bus ports + }; + + /// Signal inputs for the ESDC1A exciter model. + enum class Esdc1aSignalInputs : size_t + { + speed, ///< \f$\omega\f$ Known machine speed-deviation input [p.u.]; required when \f$s_{\mathrm{spd}}=1\f$ + vref, ///< \f$V_{\mathrm{ref}}\f$ Optional Unknown voltage-reference input [p.u.] + vs, ///< \f$V_S\f$ Optional Known stabilizer input [p.u.] + vuel, ///< \f$V_{\mathrm{UEL}}\f$ Optional Known UEL input [p.u.] + SIZE ///< Number of ESDC1A input-signal ports + }; + + /// Signal outputs for the ESDC1A exciter model. + enum class Esdc1aSignalOutputs : size_t + { + efd, ///< \f$E_{\mathrm{fd}}\f$ Required Known field-voltage output [p.u.] + SIZE ///< Number of ESDC1A output-signal ports + }; + + /// Variables available through the monitor interface. + enum class Esdc1aMonitorableVariables + { + efd, ///< \f$E_{\mathrm{fd}}\f$ Field-voltage output [p.u.] + vc, ///< \f$V_C\f$ Filtered terminal-voltage magnitude [p.u.] + vr, ///< \f$V_R\f$ Voltage-regulator output [p.u.] + vf, ///< \f$V_F\f$ Stabilizing feedback state [p.u.] + se, ///< \f$S_E\f$ Exciter saturation coefficient [p.u.] + vfe ///< \f$V_{\mathrm{FE}}\f$ Exciter feedback drive [p.u.] + }; + + /** + * @brief Model data for ESDC1A parameters, terminal bus, signal ports, + * and monitored variables. + * + * @tparam real_type Real parameter value type. + * @tparam index_type Integer index type. + * + * @see Esdc1a + */ + template + struct Esdc1aData : public ComponentData + { + Esdc1aData() = default; + + using Parameters = Esdc1aParameters; + using Buses = Esdc1aBuses; + using SignalInputs = Esdc1aSignalInputs; + using SignalOutputs = Esdc1aSignalOutputs; + using MonitorableVariables = Esdc1aMonitorableVariables; + }; + } // namespace Exciter + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aDependencyTracking.cpp b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aDependencyTracking.cpp new file mode 100644 index 000000000..dced64f5b --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aDependencyTracking.cpp @@ -0,0 +1,33 @@ +/** + * @file Esdc1aDependencyTracking.cpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Dependency-tracking instantiations for the ESDC1A exciter model. + */ + +#include "Esdc1aImpl.hpp" + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Exciter + { + /** + * @brief Report that dependency tracking does not assemble a separate Jacobian. + * + * Dependency tracking exposes the Jacobian structure through the + * residual rather than a separately assembled matrix. + */ + template + int Esdc1a::evaluateJacobian() + { + Log::misc() << "Evaluate Jacobian for Esdc1a..." << std::endl; + Log::misc() << "Jacobian evaluation is not implemented!" << std::endl; + return 0; + } + + template class Esdc1a; + template class Esdc1a; + } // namespace Exciter + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aEnzyme.cpp b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aEnzyme.cpp new file mode 100644 index 000000000..812ea4031 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aEnzyme.cpp @@ -0,0 +1,116 @@ +/** + * @file Esdc1aEnzyme.cpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Enzyme sparse Jacobian for the ESDC1A exciter model. + */ + +#include + +#include "Esdc1aImpl.hpp" + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Exciter + { + /** + * @brief Assemble the sparse ESDC1A component Jacobian with Enzyme. + * + * Differentiates the internal residual with respect to internal states, + * state derivatives, terminal-bus variables, and linked signal values, + * then assembles the resulting entries in COO form. + * + * @pre allocate() has completed. + * @pre evaluateResidual() has refreshed the interface buffers at the + * current state. + * @pre Solver alpha and global variable and residual indices are set. + */ + template + int Esdc1a::evaluateJacobian() + { + Log::misc() << "Evaluate Jacobian for Esdc1a..." << std::endl; + Log::misc() << "Jacobian evaluation is experimental!" << std::endl; + + if (J_rows_buffer_ == nullptr) + { + auto size = static_cast(size_); + auto bus_size = static_cast(bus_->size()); + auto signal_size = ws_.size(); + auto buffer_size = 2 * size * size + size * bus_size + size * signal_size; + J_rows_buffer_ = new IdxT[buffer_size]; + J_cols_buffer_ = new IdxT[buffer_size]; + J_vals_buffer_ = new RealT[buffer_size]; + } + + using ModelT = GridKit::PhasorDynamics::Exciter::Esdc1a; + using Fn = GridKit::Enzyme::Sparse::MemberFunctions; + + nnz_ = 0; + + GridKit::Enzyme::Sparse::DfDy::eval(this, + static_cast(f_.getSize()), + static_cast(y_.getSize()), + (this->getResidualIndices()).data(), + (this->getVariableIndices()).data(), + y_.getData(), + yp_.getData(), + wb_.data(), + ws_.data(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); + + GridKit::Enzyme::Sparse::DfDyp::eval(this, + static_cast(f_.getSize()), + static_cast(y_.getSize()), + (this->getResidualIndices()).data(), + (this->getVariableIndices()).data(), + y_.getData(), + yp_.getData(), + wb_.data(), + ws_.data(), + alpha_, + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); + + GridKit::Enzyme::Sparse::DfDwb::eval(this, + static_cast(f_.getSize()), + static_cast(bus_->size()), + (this->getResidualIndices()).data(), + (bus_->getVariableIndices()).data(), + y_.getData(), + yp_.getData(), + wb_.data(), + ws_.data(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); + + GridKit::Enzyme::Sparse::DfDws::eval(this, + static_cast(f_.getSize()), + ws_.size(), + (this->getResidualIndices()).data(), + ws_indices_.data(), + y_.getData(), + yp_.getData(), + wb_.data(), + ws_.data(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); + this->constructCoo(); + + return 0; + } + + template class Esdc1a; + template class Esdc1a; + } // namespace Exciter + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aImpl.hpp b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aImpl.hpp new file mode 100644 index 000000000..11906d3d9 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aImpl.hpp @@ -0,0 +1,893 @@ +/** + * @file Esdc1aImpl.hpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Definition of the ESDC1A exciter model. + */ + +#pragma once + +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Exciter + { + /// Logger used for ESDC1A diagnostics. + using Log = ::GridKit::Utilities::Logger; + + /** + * @brief Construct an ESDC1A exciter without parameters + * + * The model is sized but left unconfigured. Every parameter keeps its + * documented default, and no monitor is created. verify() rejects the + * model until an `efd` output node is assigned. + * + * @param[in] bus Terminal bus the exciter measures. + */ + template + Esdc1a::Esdc1a(BusT* bus) + : bus_(bus) + { + size_ = static_cast(Esdc1aInternalVariables::MAXIMUM); + setDerivedParameters(); + } + + /** + * @brief Construct an ESDC1A exciter from model data + * + * @param[in] bus Terminal bus the exciter measures. + * @param[in] data Parameters and monitored-variable selections. + */ + template + Esdc1a::Esdc1a(BusT* bus, const ModelDataT& data) + : bus_(bus), + monitor_(std::make_unique(data)) + { + initializeParameters(data); + initializeMonitor(); + size_ = static_cast(Esdc1aInternalVariables::MAXIMUM); + } + + /** + * @brief Destroy the ESDC1A exciter. + */ + template + Esdc1a::~Esdc1a() + { + } + + /** + * @brief Set the component ID + * + * @param[in] component_id Identifier assigned by the system model. + * @return Zero on success. + */ + template + int Esdc1a::setGridKitComponentID(IdxT component_id) + { + gridkit_component_id_ = component_id; + return 0; + } + + /** + * @brief Allocate the model vectors and wire the field-voltage output + * + * Sizes the state, residual, bus-interface, and signal-interface + * buffers, seeds the identity index maps, and points an assigned `efd` + * node at the internal field-voltage state. That node aliases ESDC1A + * storage from here on, which is how initialize() reads the seed a + * machine model wrote. Repeated calls reuse the allocated vectors. + * + * @return Zero on success. + */ + template + int Esdc1a::allocate() + { + const auto EFD = static_cast(Esdc1aInternalVariables::EFD); + + if (!allocated_) + { + this->allocateVectors(size_); + } + auto size = static_cast(size_); + + tag_.assign(size, false); + variable_indices_.resize(size); + residual_indices_.resize(size); + + wb_.assign(2, ScalarT{0}); + + const auto signal_size = static_cast(Esdc1aExternalVariables::MAXIMUM); + ws_.assign(signal_size, ScalarT{0}); + ws_indices_.assign(signal_size, INVALID_INDEX); + + for (IdxT j = 0; j < size_; ++j) + { + this->setVariableIndex(j, j); + this->setResidualIndex(j, j); + } + + auto* y = y_.getData(); + + if (signals_.template isAssigned()) + { + signals_.template getSignalNode()->set( + &y[EFD], + &(this->getVariableIndex(static_cast(EFD)))); + } + + allocated_ = true; + return 0; + } + + /** + * @brief Validate the ESDC1A configuration + * + * Checks parameter-loading errors, static parameter relationships, + * terminal-bus association, the required field-voltage output, and + * attached external signals. Seed feasibility is operating-point + * dependent and is checked by initialize(). + * + * @return Number of configuration errors; zero when valid. + */ + template + int Esdc1a::verify() const + { + int ret = static_cast(parameter_error_count_); + + auto check = [&](bool condition, const char* message) + { + if (!condition) + { + Log::error() << "Esdc1a: " << message << '\n'; + ret += 1; + } + }; + + if (bus_ == nullptr) + { + Log::error() << "Esdc1a: bus pointer is null\n"; + ret += 1; + } + + check(Ka_ > ZERO, "Ka must be positive"); + check(Tc_ >= ZERO, "Tc must be non-negative"); + check(Vrmin_ <= Vrmax_, "Vrmin must be less than or equal to Vrmax"); + check(UEL_ >= static_cast(0) && UEL_ <= static_cast(3), + "UEL must be 0, 1, 2, or 3"); + + // Model data uses an exact zero to mean "saturation bypassed", so + // this is an exact comparison by intent rather than a tolerance test. + const bool saturation_disabled = + Se1_ == ZERO && Se2_ == ZERO; + + if (!saturation_disabled) + { + check(E1_ > ZERO, "E1 must be positive when saturation is enabled"); + check(E2_ > ZERO, "E2 must be positive when saturation is enabled"); + check(Se1_ > ZERO, "Se1 must be positive when saturation is enabled"); + check(Se2_ > ZERO, "Se2 must be positive when saturation is enabled"); + + const bool saturation_points_are_ordered = + (E2_ > E1_ && Se2_ > Se1_) + || (E2_ < E1_ && Se2_ < Se1_); + check(saturation_points_are_ordered, + "E1/E2 and Se1/Se2 must be ordered consistently"); + } + + if (!signals_.template isAssigned()) + { + Log::error() << "Esdc1a: required efd output signal is not assigned\n"; + ret += 1; + } + + if (Spdmlt_ && !signals_.template isAttached()) + { + Log::error() << "Esdc1a: speed signal is required when Spdmlt is enabled\n"; + ret += 1; + } + + // An attached port must resolve to writable signal storage. The + // enumerator is a template argument, so each port names itself once. + auto check_attached_signal = + [&](const char* name) + { + if (signals_.template isAttached() + && !signals_.template isLinked()) + { + Log::error() << "Esdc1a: " << name << " signal attached with no linked source\n"; + ret += 1; + } + }; + + check_attached_signal.template operator()("speed"); + check_attached_signal.template operator()("vref"); + check_attached_signal.template operator()("vs"); + check_attached_signal.template operator()("vuel"); + + return ret; + } + + /** + * @brief Initialize ESDC1A from the field-voltage output + * + * Resolves the steady internal state and voltage reference while + * preserving the seeded `efd`, latches attached Known inputs, and + * publishes the reference to an attached `vref` signal. + * + * @return Zero on success; nonzero when the configuration or operating point is rejected. + */ + template + int Esdc1a::initialize() + { + const auto EFDP = static_cast(Esdc1aInternalVariables::EFDP); + const auto VC = static_cast(Esdc1aInternalVariables::VC); + const auto VR = static_cast(Esdc1aInternalVariables::VR); + const auto VF = static_cast(Esdc1aInternalVariables::VF); + const auto XLL = static_cast(Esdc1aInternalVariables::XLL); + const auto EV = static_cast(Esdc1aInternalVariables::EV); + const auto VLL = static_cast(Esdc1aInternalVariables::VLL); + const auto VHV = static_cast(Esdc1aInternalVariables::VHV); + const auto SE = static_cast(Esdc1aInternalVariables::SE); + const auto VFE = static_cast(Esdc1aInternalVariables::VFE); + const auto EFD = static_cast(Esdc1aInternalVariables::EFD); + + bool ret = verify() == 0; + if (!ret) + { + Log::error() << "Esdc1a: cannot initialize with invalid configuration\n"; + return 1; + } + + auto* y = y_.getData(); + + // The assigned efd node aliases this entry after allocate(). Its + // seeded value remains untouched throughout initialization. + const ScalarT efd0 = y[EFD]; + + ScalarT omega0{ZERO}; + if (signals_.template isAttached()) + { + omega0 = signals_.template readExternalVariable(); + } + + ScalarT vs0{ZERO}; + if (signals_.template isAttached()) + { + vs0 = signals_.template readExternalVariable(); + } + + ScalarT vuel0{ZERO}; + if (signals_.template isAttached()) + { + vuel0 = signals_.template readExternalVariable(); + } + + const ScalarT vc0 = std::sqrt(Vr() * Vr() + Vi() * Vi()); + + ret = std::isfinite(static_cast(efd0)) + && std::isfinite(static_cast(vc0)) + && vc0 > ZERO; + if (!ret) + { + Log::error() << "Esdc1a: initial bus-voltage magnitude must be finite and positive, " + "and the field-voltage seed must be finite\n"; + return 1; + } + + ret = std::isfinite(static_cast(omega0)) + && std::isfinite(static_cast(vs0)) + && std::isfinite(static_cast(vuel0)); + if (!ret) + { + Log::error() << "Esdc1a: initial speed, stabilizer, and UEL inputs must be finite\n"; + return 1; + } + + const ScalarT speed_multiplier = ONE + spd_on_ * omega0; + + ret = std::isfinite(static_cast(speed_multiplier)) + && speed_multiplier > ZERO; + if (!ret) + { + Log::error() << "Esdc1a: speed multiplier must be finite and positive at initialization\n"; + return 1; + } + + const ScalarT efdp0 = efd0 / speed_multiplier; + + ret = !exclim_ || efdp0 >= ZERO; + if (!ret) + { + Log::error() << "Esdc1a: initial Efd' is below its enabled zero limit\n"; + return 1; + } + + const ScalarT se0 = SB_ * Math::qramp(efdp0 - SA_); + const ScalarT vfe0 = (Ke_ + se0) * efdp0; + const ScalarT vr0 = vfe0; + const ScalarT vhv0 = vr0 / Ka_; + + ret = vr0 >= Vrmin_ && vr0 <= Vrmax_; + if (!ret) + { + Log::error() << "Esdc1a: initialized VR is outside limits\n"; + return 1; + } + + // An inactive high-value gate is seeded with the gate input, so the + // residual reproduces VHV through the same smooth maximum. UEL modes + // below 2 route VUEL through the gate; setDerivedParameters() derives + // the uel_on_ blend mask from the same threshold. + ScalarT vll0 = vhv0; + if (UEL_ < static_cast(2)) + { + const RealT gate_margin0 = static_cast(vhv0 - vuel0); + + ret = gate_margin0 > ZERO; + if (!ret) + { + Log::error() << "Esdc1a: smooth high-value gate is active at initialization\n"; + return 1; + } + vll0 = vuel0 + inverseRamp(gate_margin0); + } + + const ScalarT vf0 = ScalarT{ZERO}; + const ScalarT ev0 = vll0; + const ScalarT xll0 = ev0; + const ScalarT vref0 = ev0 + vc0 + vf0 - vs0 - uel_on_ * vuel0; + + y[EFDP] = efdp0; + y[VC] = vc0; + y[VR] = vr0; + y[VF] = vf0; + y[XLL] = xll0; + y[EV] = ev0; + y[VLL] = vll0; + y[VHV] = vhv0; + y[SE] = se0; + y[VFE] = vfe0; + y[EFD] = efd0; + + omega_set_ = omega0; + vref_set_ = vref0; + vs_set_ = vs0; + vuel_set_ = vuel0; + + if (signals_.template isAttached()) + { + signals_.template writeExternalVariable(vref_set_); + } + + y_.setDataUpdated(); + yp_.setToConst(static_cast(ZERO)); + return 0; + } + + /** + * @brief Identify the differential variables + * + * The field-voltage state, the voltage transducer, the regulator, the + * stabilizing feedback, and the lead-lag state carry derivatives; + * every other internal variable is algebraic. + * + * @return Zero on success. + */ + template + int Esdc1a::tagDifferentiable() + { + const auto EFDP = static_cast(Esdc1aInternalVariables::EFDP); + const auto VC = static_cast(Esdc1aInternalVariables::VC); + const auto VR = static_cast(Esdc1aInternalVariables::VR); + const auto VF = static_cast(Esdc1aInternalVariables::VF); + const auto XLL = static_cast(Esdc1aInternalVariables::XLL); + + std::fill(tag_.begin(), tag_.end(), false); + tag_[EFDP] = true; + tag_[VC] = true; + tag_[VR] = true; + tag_[VF] = true; + tag_[XLL] = true; + return 0; + } + + /** + * @brief Compute the absolute tolerance for each variable in the model + * + * Every internal variable receives @p rel_tol as its absolute + * tolerance. + * + * @param[in] rel_tol Solver relative tolerance. + * @return Zero on success. + */ + template + int Esdc1a::setAbsoluteTolerance(RealT rel_tol) + { + abs_tol_.setToConst(static_cast(rel_tol)); + return 0; + } + + /** + * @brief Residuals of system equations + * + * Refreshes the bus and signal interface buffers and evaluates the + * internal residual. ESDC1A injects no current, so there is no bus + * residual. An unattached input port falls back to the value latched + * by initialize(). + * + * @return Zero on success. + */ + template + int Esdc1a::evaluateResidual() + { + const auto OMEGA = static_cast(Esdc1aExternalVariables::OMEGA); + const auto VREF = static_cast(Esdc1aExternalVariables::VREF); + const auto VS = static_cast(Esdc1aExternalVariables::VS); + const auto VUEL = static_cast(Esdc1aExternalVariables::VUEL); + + ws_[OMEGA] = omega_set_; + ws_[VREF] = vref_set_; + ws_[VS] = vs_set_; + ws_[VUEL] = vuel_set_; + std::fill(ws_indices_.begin(), ws_indices_.end(), INVALID_INDEX); + + if (signals_.template isAttached()) + { + ws_[OMEGA] = signals_.template readExternalVariable(); + ws_indices_[OMEGA] = + signals_.template readExternalVariableIndex(); + } + if (signals_.template isAttached()) + { + ws_[VREF] = signals_.template readExternalVariable(); + ws_indices_[VREF] = + signals_.template readExternalVariableIndex(); + } + if (signals_.template isAttached()) + { + ws_[VS] = signals_.template readExternalVariable(); + ws_indices_[VS] = + signals_.template readExternalVariableIndex(); + } + if (signals_.template isAttached()) + { + ws_[VUEL] = signals_.template readExternalVariable(); + ws_indices_[VUEL] = + signals_.template readExternalVariableIndex(); + } + + wb_[0] = Vr(); + wb_[1] = Vi(); + + const auto* y = y_.getData(); + const auto* yp = yp_.getData(); + auto* f = f_.getData(); + + evaluateInternalResidual(y, yp, wb_.data(), ws_.data(), f); + f_.setDataUpdated(); + return 0; + } + + /** + * @brief Access the monitor + * + * @return Monitor for this model, or nullptr when the model was + * constructed without data. + */ + template + const Model::VariableMonitorBase* Esdc1a::getMonitor() const + { + return monitor_.get(); + } + + /** + * @brief Evaluate the ESDC1A internal residual. + * + * Evaluates the five exciter states and the six algebraic equations + * documented in the model README. The body is kept free of branches + * and loops so sparse automatic differentiation resolves a fixed + * structure; the three selector decisions enter as multiplicative + * masks set by setDerivedParameters(). + * + * @param[in] y Internal variables in Esdc1aInternalVariables order. + * @param[in] yp Internal derivatives in the same enum order. + * @param[in] wb Terminal-bus \f$(V_{\mathrm{r}},V_{\mathrm{i}})\f$ + * voltage components. + * @param[in] ws Signal values in Esdc1aExternalVariables order. + * @param[out] f Residuals in Esdc1aInternalVariables order. + * @return Zero on success. + */ + template + __attribute__((always_inline)) inline int + Esdc1a::evaluateInternalResidual( + const ScalarT* y, + const ScalarT* yp, + const ScalarT* wb, + const ScalarT* ws, + ScalarT* f) + { + const auto EFDP = static_cast(Esdc1aInternalVariables::EFDP); + const auto VC = static_cast(Esdc1aInternalVariables::VC); + const auto VR = static_cast(Esdc1aInternalVariables::VR); + const auto VF = static_cast(Esdc1aInternalVariables::VF); + const auto XLL = static_cast(Esdc1aInternalVariables::XLL); + const auto EV = static_cast(Esdc1aInternalVariables::EV); + const auto VLL = static_cast(Esdc1aInternalVariables::VLL); + const auto VHV = static_cast(Esdc1aInternalVariables::VHV); + const auto SE = static_cast(Esdc1aInternalVariables::SE); + const auto VFE = static_cast(Esdc1aInternalVariables::VFE); + const auto EFD = static_cast(Esdc1aInternalVariables::EFD); + + const auto OMEGA = static_cast(Esdc1aExternalVariables::OMEGA); + const auto VREF = static_cast(Esdc1aExternalVariables::VREF); + const auto VS = static_cast(Esdc1aExternalVariables::VS); + const auto VUEL = static_cast(Esdc1aExternalVariables::VUEL); + + const ScalarT efdp = y[EFDP]; + const ScalarT vc = y[VC]; + const ScalarT vr = y[VR]; + const ScalarT vf = y[VF]; + const ScalarT xll = y[XLL]; + const ScalarT ev = y[EV]; + const ScalarT vll = y[VLL]; + const ScalarT vhv = y[VHV]; + const ScalarT se = y[SE]; + const ScalarT vfe = y[VFE]; + const ScalarT efd = y[EFD]; + + const ScalarT efdp_dot = yp[EFDP]; + const ScalarT vc_dot = yp[VC]; + const ScalarT vr_dot = yp[VR]; + const ScalarT vf_dot = yp[VF]; + const ScalarT xll_dot = yp[XLL]; + + const ScalarT omega = ws[OMEGA]; + const ScalarT vref = ws[VREF]; + const ScalarT vs = ws[VS]; + const ScalarT vuel = ws[VUEL]; + + const ScalarT ec = std::sqrt(wb[0] * wb[0] + wb[1] * wb[1]); + const ScalarT ev_target = vref + vs + uel_on_ * vuel - vc - vf; + const ScalarT vfe_target = (Ke_ + se) * efdp; + const ScalarT efdp_rate = (vr - vfe) / Te_; + const ScalarT limited_efdp_rate = awmin(efdp, efdp_rate, ZERO); + + f[EFDP] = -efdp_dot + (ONE - lim_on_) * efdp_rate + + lim_on_ * limited_efdp_rate; + f[VC] = -vc_dot + (ec - vc) / Tr_; + f[VR] = -vr_dot + Math::antiwindup(vr, -vr + Ka_ * vhv, Vrmin_, Vrmax_) / Ta_; + f[VF] = -vf_dot + (-vf + Kf_ * (vr - vfe) / Te_) / Tf1_; + f[XLL] = -xll_dot + (ev - xll) / Tb_; + f[EV] = -ev + ev_target; + f[VLL] = -vll + xll + (Tc_ / Tb_) * (ev - xll); + f[VHV] = -vhv + uel_on_ * vll + + (ONE - uel_on_) * Math::max(vll, vuel); + f[SE] = -se + SB_ * Math::qramp(efdp - SA_); + f[VFE] = -vfe + vfe_target; + f[EFD] = -efd + (ONE + spd_on_ * omega) * efdp; + + return 0; + } + + // + // Private methods + // + + /** + * @brief Smooth anti-windup derivative above a fixed lower bound + * + * Passes the unconstrained rate above the bound, admits restoring + * motion from below it, and smoothly blocks outward motion. + * + * @param[in] x State limited from below. + * @param[in] f Unconstrained derivative of @p x. + * @param[in] xmin Fixed lower bound on @p x. + * @return Anti-windup-limited derivative. + */ + template + __attribute__((always_inline)) inline scalar_type + Esdc1a::awmin( + const ScalarT x, + const ScalarT f, + const RealT xmin) + { + const ScalarT above = Math::above(x, xmin); + + return (above + (ONE - above) * Math::sigmoid(f)) * f; + } + + /** + * @brief Read the parameters out of the model data + * + * No parameter is required; every parameter keeps the default + * documented in the model README when omitted. A non-numeric value, a + * switch outside \f$\{0,1\}\f$, or a non-integer selector is counted and + * reported by verify() rather than throwing. Integer JSON values are + * accepted for real parameters. + * + * @param[in] data Parameters and monitored-variable selections. + */ + template + void Esdc1a::initializeParameters(const ModelDataT& data) + { + using Params = typename ModelDataT::Parameters; + + parameter_error_count_ = 0; + + auto load_real = [&](auto key, RealT& target, const char* name) + { + if (!data.parameters.contains(key)) + { + return; + } + + const auto& value = data.parameters.at(key); + RealT parsed_value{}; + if (const auto* real_value = std::get_if(&value)) + { + parsed_value = *real_value; + } + else if (const auto* index_value = std::get_if(&value)) + { + parsed_value = static_cast(*index_value); + } + else + { + Log::error() << "Esdc1a: parameter '" << name << "' must be numeric\n"; + ++parameter_error_count_; + return; + } + + const bool ret = std::isfinite(parsed_value); + if (!ret) + { + Log::error() << "Esdc1a: parameter '" << name << "' must be finite\n"; + ++parameter_error_count_; + return; + } + + target = parsed_value; + }; + + auto load_switch = [&](auto key, bool& target, const char* name) + { + if (!data.parameters.contains(key)) + { + return; + } + + const auto& value = data.parameters.at(key); + if (const auto* bool_value = std::get_if(&value)) + { + target = *bool_value; + } + else + { + Log::error() << "Esdc1a: parameter '" << name << "' must be boolean\n"; + ++parameter_error_count_; + } + }; + + auto load_selector = [&](auto key, IdxT& target, const char* name) + { + if (!data.parameters.contains(key)) + { + return; + } + + const auto& value = data.parameters.at(key); + if (const auto* index_value = std::get_if(&value)) + { + target = *index_value; + } + else + { + Log::error() << "Esdc1a: parameter '" << name << "' must be an integer selector\n"; + ++parameter_error_count_; + } + }; + + load_real(Params::Tr, Tr_, "Tr"); + load_real(Params::Ka, Ka_, "Ka"); + load_real(Params::Ta, Ta_, "Ta"); + load_real(Params::Tb, Tb_, "Tb"); + load_real(Params::Tc, Tc_, "Tc"); + load_real(Params::Vrmax, Vrmax_, "Vrmax"); + load_real(Params::Vrmin, Vrmin_, "Vrmin"); + load_real(Params::Ke, Ke_, "Ke"); + load_real(Params::Te, Te_, "Te"); + load_real(Params::Kf, Kf_, "Kf"); + load_real(Params::Tf1, Tf1_, "Tf1"); + load_switch(Params::Spdmlt, Spdmlt_, "Spdmlt"); + load_real(Params::E1, E1_, "E1"); + load_real(Params::Se1, Se1_, "Se1"); + load_real(Params::E2, E2_, "E2"); + load_real(Params::Se2, Se2_, "Se2"); + load_selector(Params::UEL, UEL_, "UEL"); + load_switch(Params::exclim, exclim_, "exclim"); + setDerivedParameters(); + } + + /** + * @brief Bind the monitorable variables to their internal states + * + * Binds configured monitor keys to their corresponding internal + * variables. + */ + template + void Esdc1a::initializeMonitor() + { + using Variable = typename ModelDataT::MonitorableVariables; + + monitor_->set(Variable::efd, [this] + { return y_.getData()[static_cast(Esdc1aInternalVariables::EFD)]; }); + monitor_->set(Variable::vc, [this] + { return y_.getData()[static_cast(Esdc1aInternalVariables::VC)]; }); + monitor_->set(Variable::vr, [this] + { return y_.getData()[static_cast(Esdc1aInternalVariables::VR)]; }); + monitor_->set(Variable::vf, [this] + { return y_.getData()[static_cast(Esdc1aInternalVariables::VF)]; }); + monitor_->set(Variable::se, [this] + { return y_.getData()[static_cast(Esdc1aInternalVariables::SE)]; }); + monitor_->set(Variable::vfe, [this] + { return y_.getData()[static_cast(Esdc1aInternalVariables::VFE)]; }); + } + + /** + * @brief Resolve the parameter-derived constants and selector masks + * + * Raises the transducer, regulator, lead-lag, exciter, and feedback + * lags to the well-posedness floor, fits the quadratic saturation + * curve, and turns the three selectors into multiplicative masks. The + * masks let the residual select signal routing without + * parameter-dependent control flow, which keeps its structure fixed for + * sparse automatic differentiation. + */ + template + void Esdc1a::setDerivedParameters() + { + // The lags are raised to the floor in place, so a negative value is + // rejected here while the value as read is still available. verify() + // reports the count. + auto check_non_negative = [&](RealT value, const char* name) + { + if (value < ZERO) + { + Log::error() << "Esdc1a: " << name << " must be non-negative\n"; + ++parameter_error_count_; + } + }; + + check_non_negative(Tr_, "Tr"); + check_non_negative(Ta_, "Ta"); + check_non_negative(Tb_, "Tb"); + check_non_negative(Te_, "Te"); + check_non_negative(Tf1_, "Tf1"); + + if (Tr_ < TIME_CONSTANT_MINIMUM || Ta_ < TIME_CONSTANT_MINIMUM + || Tb_ < TIME_CONSTANT_MINIMUM || Te_ < TIME_CONSTANT_MINIMUM + || Tf1_ < TIME_CONSTANT_MINIMUM) + { + Log::warning() << "Esdc1a: Tr, Ta, Tb, Te, and Tf1 below " + << TIME_CONSTANT_MINIMUM + << " s are raised to that floor to keep the exciter lags well posed\n"; + } + + Tr_ = std::max(Tr_, TIME_CONSTANT_MINIMUM); + Ta_ = std::max(Ta_, TIME_CONSTANT_MINIMUM); + Tb_ = std::max(Tb_, TIME_CONSTANT_MINIMUM); + Te_ = std::max(Te_, TIME_CONSTANT_MINIMUM); + Tf1_ = std::max(Tf1_, TIME_CONSTANT_MINIMUM); + + spd_on_ = ZERO; + if (Spdmlt_) + { + spd_on_ = ONE; + } + + uel_on_ = ZERO; + if (UEL_ >= static_cast(2)) + { + uel_on_ = ONE; + } + + lim_on_ = ZERO; + if (exclim_) + { + lim_on_ = ONE; + } + + // A disabled or inconsistent saturation curve keeps the zero fit so + // the coefficients stay finite; verify() reports inconsistent data. + // The disabled test matches the verify() predicate exactly. + const bool saturation_disabled = + Se1_ == ZERO && Se2_ == ZERO; + const bool saturation_points_are_ordered = + (E2_ > E1_ && Se2_ > Se1_) + || (E2_ < E1_ && Se2_ < Se1_); + const bool saturation_consistent = + E1_ > ZERO && E2_ > ZERO + && Se1_ > ZERO && Se2_ > ZERO + && saturation_points_are_ordered; + if (saturation_disabled || !saturation_consistent) + { + SA_ = ZERO; + SB_ = ZERO; + return; + } + + const RealT C = std::sqrt(Se2_ / Se1_); + SA_ = (C * E1_ - E2_) / (C - ONE); + SB_ = Se1_ / ((E1_ - SA_) * (E1_ - SA_)); + } + + /** + * @brief Invert the smooth CommonMath ramp + * + * Initialization seeds the inactive high-value gate with the gate + * *input*, so the residual reproduces the requested output through the + * same smooth ramp it evaluates. + * + * For large positive values, the ramp is effectively equal to the input, so the + * inverse is effectively the output. In that regime this function returns `ramp_output` directly. + * This branching is numerically more robust. + * + * @param[in] ramp_output Strictly positive requested ramp output. + * @return The input the smooth ramp maps to the requested output. + * + * @pre @p ramp_output is finite and strictly positive. + * @warning This function contains conditional branching and may be used + * during initialization, but not during residual or Jacobian + * evaluation. + */ + template + typename Esdc1a::RealT + Esdc1a::inverseRamp(RealT ramp_output) const + { + static constexpr RealT SOFTPLUS_WIDTH = static_cast(50.0); + + const RealT scaled_output = Math::MU * ramp_output; + if (scaled_output > SOFTPLUS_WIDTH) + { + return ramp_output; + } + return std::log(std::expm1(scaled_output)) / Math::MU; + } + + /** + * @brief Access the terminal-bus real voltage component. + * + * @return Reference to the bus variable. + */ + template + scalar_type& Esdc1a::Vr() + { + return bus_->Vr(); + } + + /** + * @brief Access the terminal-bus imaginary voltage component. + * + * @return Reference to the bus variable. + */ + template + scalar_type& Esdc1a::Vi() + { + return bus_->Vi(); + } + + } // namespace Exciter + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/README.md b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/README.md new file mode 100644 index 000000000..aed2f0d58 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/README.md @@ -0,0 +1,414 @@ +# **IEEE DC1A Excitation System Model (ESDC1A)** + +ESDC1A is an IEEE DC1A excitation-system model with a voltage transducer, +input lead-lag compensation, a limited voltage regulator, exciter feedback and +saturation, under-excitation limiter routing, and an optional speed multiplier. + +## Notes + +- Internal voltage signals are on component base. +- The source diagram labels the optional multiplier input as `Speed`; GridKit + uses machine speed deviation, so the enabled multiplier is $1+\omega$. +- The UEL selector routes $V_{\mathrm{UEL}}$ either through the high-value gate + or through the voltage-error summing junction. + +## Block Diagram + +![ESDC1A exciter block diagram](../../../../../docs/Figures/PhasorDynamics/ESDC1A/diagram.png) + +Figure 1: ESDC1A exciter model. Figure courtesy of the +[PowerWorld ESDC1A model reference](https://www.powerworld.com/WebHelp/Content/TransientModels_HTML/Exciter%20ESDC1A.htm). + +## Model Parameters + +Symbol | Units | JSON | Description | Typical Value +------------------------------------|-----------|-----------|-------------------------------------------------|-------------- +$T_R$ | [sec] | `Tr` | Voltage transducer time constant | 0.0 +$K_A$ | [p.u.] | `Ka` | Voltage-regulator gain | 40.0 +$T_A$ | [sec] | `Ta` | Voltage-regulator time constant | 0.1 +$T_B$ | [sec] | `Tb` | Input lead-lag denominator time constant | 0.0 +$T_C$ | [sec] | `Tc` | Input lead-lag numerator time constant | 0.0 +$V_R^{\max}$ | [p.u.] | `Vrmax` | Maximum voltage-regulator output | 1.0 +$V_R^{\min}$ | [p.u.] | `Vrmin` | Minimum voltage-regulator output | -1.0 +$K_E$ | [p.u.] | `Ke` | Exciter constant | 0.1 +$T_E$ | [sec] | `Te` | Exciter time constant | 0.5 +$K_F$ | [p.u.] | `Kf` | Stabilizing feedback gain | 0.05 +$T_{F1}$ | [sec] | `Tf1` | Stabilizing feedback time constant | 0.7 +$s_{\mathrm{spd}}$ | [boolean] | `Spdmlt` | Field-voltage speed-multiplier flag | `false` +$E_1$ | [p.u.] | `E1` | First saturation voltage point | 2.8 +$S_E(E_1)$ | [p.u.] | `Se1` | Saturation coefficient at $E_1$ | 0.08 +$E_2$ | [p.u.] | `E2` | Second saturation voltage point | 3.7 +$S_E(E_2)$ | [p.u.] | `Se2` | Saturation coefficient at $E_2$ | 0.33 +$I_{\mathrm{UEL}}$ | [integer] | `UEL` | Under-excitation limiter input-routing selector | 0 +$s_{\mathrm{lim}}$ | [boolean] | `exclim` | Exciter field-voltage-state lower-limit flag | `true` + +Every parameter is optional. +All real-valued parameters must be finite. `Spdmlt` and `exclim` must be +JSON booleans, and `UEL` must be a JSON integer. + +### Parameter Validation + +Invalid ESDC1A parameter sets are rejected by the following checks: + +```math +\begin{aligned} + K_A + &> 0 \\ + T_R, T_A, T_B, T_C, T_E, T_{F1} + &\ge 0 \\ + V_R^{\min} + &\le V_R^{\max} \\ + s_{\mathrm{spd}}, s_{\mathrm{lim}} + &\in \{0,1\} \\ + I_{\mathrm{UEL}} + &\in \{0,1,2,3\} +\end{aligned} +``` + +The saturation points are either disabled together, + +```math +S_E(E_1) = S_E(E_2) = 0, +``` + +or define a valid two-point quadratic fit: + +```math +\begin{aligned} + E_1, E_2, S_E(E_1), S_E(E_2) &> 0 \\ + \left(E_2-E_1\right) + \left[S_E(E_2)-S_E(E_1)\right] &> 0 +\end{aligned} +``` + +### Model Derived Parameters + +Let $\epsilon_T = 10^{-3}\ \mathrm{s}$. A time constant below $\epsilon_T$ is +raised to that floor in place, so every equation below uses the raised value: + +```math +\begin{aligned} + T_x + &\leftarrow \max\!\left(T_x,\epsilon_T\right), + \quad x\in\{R,A,B,E,F1\} \\ + s_{\mathrm{UEL}} + &= + \begin{cases} + 1 & I_{\mathrm{UEL}} \ge 2 \\ + 0 & I_{\mathrm{UEL}} < 2 + \end{cases} +\end{aligned} +``` + +When saturation is disabled, $S_A = 0$ and $S_B = 0$. Otherwise, + +```math +\begin{aligned} + C &= \sqrt{\dfrac{S_E(E_2)}{S_E(E_1)}} \\ + S_A &= \dfrac{C E_1 - E_2}{C - 1} \\ + S_B &= \dfrac{S_E(E_1)}{(E_1 - S_A)^2} +\end{aligned} +``` + +## Model Ports + +Name | Port | Init | Description +--------|--------|---------|------ +`bus` | Bus | Known | Terminal bus voltage +`speed` | Input | Known | Machine speed deviation +`vref` | Input | Unknown | Voltage-control reference +`vs` | Input | Known | Stabilizer input signal +`vuel` | Input | Known | Under-excitation limiter input +`efd` | Output | Known | Field-voltage output + +`Known` ports hold their initial values before `initialize()` and are preserved +by it. `Unknown` inputs are resolved during initialization and written to +attached signal storage, or retained as constant inputs when unattached. The +`efd` output must be assigned. The `speed` input is required when +$s_{\mathrm{spd}} = 1$; every other signal input is optional. Unattached `speed`, +`vs`, and `vuel` inputs default to zero. + +## Model Variables + +### Internal Variables + +#### Differential + +Symbol | Units | Description | Note +------------------------------------|--------|--------------------------------------------------|------ +$E_{\mathrm{fd}}'$ | [p.u.] | Exciter field-voltage state | State 1 in Fig. 1; lower bounded at zero when $s_{\mathrm{lim}} = 1$; before the optional speed multiplier +$V_C$ | [p.u.] | Filtered terminal-voltage magnitude | State 2 in Fig. 1 +$V_R$ | [p.u.] | Voltage-regulator output | State 3 in Fig. 1 +$V_F$ | [p.u.] | Stabilizing feedback state | State 4 in Fig. 1 +$x_{\mathrm{LL}}$ | [p.u.] | Input lead-lag denominator state | State 5 in Fig. 1 + +#### Algebraic + +Symbol | Units | Description | Note +------------------------------------|--------|-----------------------------------|------ +$e_V$ | [p.u.] | Voltage-error summing output | +$V_{\mathrm{LL}}$ | [p.u.] | Input lead-lag output | +$V_{\mathrm{HV}}$ | [p.u.] | High-value gate output | +$S_E$ | [p.u.] | Exciter saturation coefficient | Evaluated at $E_{\mathrm{fd}}'$ +$V_{\mathrm{FE}}$ | [p.u.] | Exciter feedback drive | +$E_{\mathrm{fd}}$ | [p.u.] | Field-voltage output | Published through `efd` + +### External Variables + +#### Differential + +None. + +#### Algebraic + +Symbol | Units | Init | Description | Note +------------------------------------|--------|---------|----------------------------------------|------ +$V_{\mathrm{r}}$ | [p.u.] | Known | Terminal voltage, real component | Bus input +$V_{\mathrm{i}}$ | [p.u.] | Known | Terminal voltage, imaginary component | Bus input +$\omega$ | [p.u.] | Known | Machine speed deviation | Signal port `speed` +$V_{\mathrm{ref}}$ | [p.u.] | Unknown | Voltage-control reference | Signal port `vref` +$V_S$ | [p.u.] | Known | Stabilizer input signal | Signal port `vs` +$V_{\mathrm{UEL}}$ | [p.u.] | Known | Under-excitation limiter input | Signal port `vuel` + +## Model Equations + +Define the pre-limit exciter field-voltage rate: + +```math +f_E = \dfrac{V_R-V_{\mathrm{FE}}}{T_E}. +``` + +### Differential Equations + +```math +\begin{aligned} + 0 &= + -\dot{E}_{\mathrm{fd}}' + + \left(1-s_{\mathrm{lim}}\right)f_E + + s_{\mathrm{lim}}\, + \text{awmin}\left(E_{\mathrm{fd}}',f_E;0\right) \\ + 0 &= + -\dot{V}_C + + \dfrac{1}{T_R} + \left( + \sqrt{V_{\mathrm{r}}^2+V_{\mathrm{i}}^2} + - V_C + \right) \\ + 0 &= + -\dot{V}_R + + \dfrac{1}{T_A} + \text{antiwindup} + \left( + V_R,\, + -V_R + K_A V_{\mathrm{HV}};\, + V_R^{\min}, V_R^{\max} + \right) \\ + 0 &= + -\dot{V}_F + + \dfrac{1}{T_{F1}} + \left[ + -V_F + + \dfrac{K_F}{T_E} + \left(V_R - V_{\mathrm{FE}}\right) + \right] \\ + 0 &= + -\dot{x}_{\mathrm{LL}} + + \dfrac{1}{T_B} + \left(e_V - x_{\mathrm{LL}}\right) +\end{aligned} +``` + +The field-voltage-state limiter uses the fixed-lower-bound anti-windup rule +of [Appendix A](#appendix-a-awmin). + +### Algebraic Equations + +```math +\begin{aligned} + 0 &= + -e_V + + V_{\mathrm{ref}} + + V_S + + s_{\mathrm{UEL}}V_{\mathrm{UEL}} + - V_C + - V_F \\ + 0 &= + -V_{\mathrm{LL}} + + x_{\mathrm{LL}} + + \dfrac{T_C}{T_B} + \left(e_V - x_{\mathrm{LL}}\right) \\ + 0 &= + -V_{\mathrm{HV}} + + \begin{cases} + \text{max}\left(V_{\mathrm{LL}}, V_{\mathrm{UEL}}\right) + & s_{\mathrm{UEL}} = 0 \\ + V_{\mathrm{LL}} + & s_{\mathrm{UEL}} = 1 + \end{cases} \\ + 0 &= + -S_E + + S_B q\left(E_{\mathrm{fd}}' - S_A\right) \\ + 0 &= + -V_{\mathrm{FE}} + + \left(K_E + S_E\right)E_{\mathrm{fd}}' \\ + 0 &= + -E_{\mathrm{fd}} + + \left(1 + s_{\mathrm{spd}}\omega\right)E_{\mathrm{fd}}' +\end{aligned} +``` + +CommonMath defines helper targets and smooth approximations for +[max](../../../../CommonMath.md#derived-functions), the [ramp](../../../../CommonMath.md#primitives) +$\rho$, and the [quadratic ramp](../../../../CommonMath.md#primitives) $q$. + +## Initialization + +### Input Initialization + +```math +\begin{aligned} + V_{\mathrm{r}}, V_{\mathrm{i}} + &\leftarrow \text{terminal-bus voltage} \\ + E_{\mathrm{fd}} + &\leftarrow \text{machine field voltage} \\ + \omega + &\leftarrow \text{machine speed deviation or }0 \\ + V_S + &\leftarrow \text{stabilizer signal or }0 \\ + V_{\mathrm{UEL}} + &\leftarrow \text{under-excitation limiter input or }0 +\end{aligned} +``` + +Initialization never replaces the seeded value held in $E_{\mathrm{fd}}$. + +### Internal Initialization + +All internal derivatives are set to zero. The steady-state residuals are then +resolved in dependency order. The smooth high-value gate requires its input to +be recovered through the inverse CommonMath +[ramp](../../../../CommonMath.md#primitives) $\rho^{-1}$ when the UEL input is +routed through the gate: + +```math +\begin{aligned} + V_C + &\leftarrow \sqrt{V_{\mathrm{r}}^2+V_{\mathrm{i}}^2} \\ + E_{\mathrm{fd}}' + &\leftarrow + \dfrac{E_{\mathrm{fd}}}{1 + s_{\mathrm{spd}}\omega} \\ + S_E + &\leftarrow S_B q\left(E_{\mathrm{fd}}' - S_A\right) \\ + V_{\mathrm{FE}} + &\leftarrow \left(K_E + S_E\right)E_{\mathrm{fd}}' \\ + V_R + &\leftarrow V_{\mathrm{FE}} \\ + V_{\mathrm{HV}} + &\leftarrow \dfrac{V_R}{K_A} \\ + V_{\mathrm{LL}} + &\leftarrow + \begin{cases} + V_{\mathrm{UEL}} + + \rho^{-1} + \left(V_{\mathrm{HV}}-V_{\mathrm{UEL}}\right) + & s_{\mathrm{UEL}} = 0 \\ + V_{\mathrm{HV}} + & s_{\mathrm{UEL}} = 1 + \end{cases} \\ + V_F + &\leftarrow 0 \\ + e_V + &\leftarrow V_{\mathrm{LL}} \\ + x_{\mathrm{LL}} + &\leftarrow e_V +\end{aligned} +``` + +Initialization rejects a non-finite or zero bus-voltage magnitude, a +non-finite field-voltage seed, non-finite Known signal inputs, a nonpositive +speed multiplier $1 + s_{\mathrm{spd}}\omega$, $E_{\mathrm{fd}}'<0$ while +$s_{\mathrm{lim}}=1$, $V_R$ outside +$[V_R^{\min},V_R^{\max}]$, and high-value-gate active starts with +$s_{\mathrm{UEL}} = 0$ and +$V_{\mathrm{HV}}\le V_{\mathrm{UEL}}$. + +Every check resolves before any storage is written, so a rejected +initialization leaves state, the `efd` seed, and external signals unchanged. + +### Output Initialization + +```math +\begin{aligned} + V_{\mathrm{ref}} + &\leftarrow + e_V + + V_C + + V_F + - V_S + - s_{\mathrm{UEL}}V_{\mathrm{UEL}} +\end{aligned} +``` + +ESDC1A writes the resolved voltage-control reference to an attached `vref` +signal input. If no controller is connected, that value is used as a constant +reference input. + +## Monitorable Outputs + +Output | Units | Description | Note +----------------|--------|-------------------------------------|------ +`efd` | [p.u.] | Field-voltage output | $E_{\mathrm{fd}}$ +`vc` | [p.u.] | Filtered terminal-voltage magnitude | $V_C$ +`vr` | [p.u.] | Voltage-regulator output | $V_R$ +`vf` | [p.u.] | Stabilizing feedback state | $V_F$ +`se` | [p.u.] | Exciter saturation coefficient | $S_E$ +`vfe` | [p.u.] | Exciter feedback drive | $V_{\mathrm{FE}}$ + +## Testing + +- `validation()` checks construction, documented defaults, parameter types + and domains, signal configuration, and minimum time-constant handling. +- `initializationAndSignals()` checks steady initialization, selector + combinations, signal publication and latching, monitor output, and + differentiability tags. +- `initializationDomain()` checks rejected and accepted field-voltage, + terminal-voltage, Known-input, speed-multiplier, regulator-limit, and + high-value-gate operating points. +- `residualEquations()` checks every model residual against a fixed + numerical answer key. +- `voltageRegulation()` checks the transducer, summing junction, lead-lag, + stabilizing feedback, and regulator anti-windup behavior. +- `excitationLimits()` checks high-value-gate routing, saturation, + field-voltage-state limiting, and the optional speed multiplier. +- `jacobian()` compares the dependency-tracking and Enzyme Jacobians when + Enzyme support is enabled. + +## Appendix A: `awmin` + +The exact anti-windup rule at a fixed lower bound $\ell$ is + +```math +\text{awmin}(x,f;\ell) = + \begin{cases} + f & x > \ell \\ + \text{max}(f,0) & x \le \ell + \end{cases} +``` + +Above the bound the unconstrained derivative passes. At or below the bound, +outward motion is blocked and restoring motion is admitted. + +The model evaluates this rule with the following smooth approximation: + +```math +\text{awmin}(x,f;\ell) + \approx + \left[ + \sigma(f) + + \left(1-\sigma(f)\right)\text{above}(x;\ell) + \right]f. +``` + +CommonMath defines the [`above`](../../../../CommonMath.md#derived-functions) +and [`sigmoid`](../../../../CommonMath.md#primitives) targets and smooth +approximations. diff --git a/GridKit/Model/PhasorDynamics/Exciter/README.md b/GridKit/Model/PhasorDynamics/Exciter/README.md index 5391918f0..69e1851d5 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/README.md +++ b/GridKit/Model/PhasorDynamics/Exciter/README.md @@ -1,7 +1,7 @@ # **Exciter Models** > [!NOTE] -> IEEET1 and SEXS-PTI exciters are currently implemented. +> IEEET1, SEXS-PTI, and ESDC1A exciters are currently implemented. ## Introduction @@ -15,6 +15,7 @@ There are a few standard Exciter models - ESAC6A Excitation Model (See [ESAC6A](ESAC6A/README.md)) - IEEE Type 1 Excitation Model (See [IEEET1](IEEET1/README.md)) - IEEE DC1 Excitation Model (See [EXDC1](EXDC1/README.md)) +- IEEE DC1A Excitation Model (See [ESDC1A](ESDC1A/README.md)) - ESDC2A Excitation Model (See [ESDC2A](ESDC2A/README.md)) - EXAC1 Excitation Model (See [EXAC1](EXAC1/README.md)) - IEEE ST4B Excitation Model (See [ESST4B](ESST4B/README.md)) diff --git a/GridKit/Model/PhasorDynamics/INPUT_FORMAT.md b/GridKit/Model/PhasorDynamics/INPUT_FORMAT.md index 69df8017f..84c03ab16 100644 --- a/GridKit/Model/PhasorDynamics/INPUT_FORMAT.md +++ b/GridKit/Model/PhasorDynamics/INPUT_FORMAT.md @@ -147,6 +147,7 @@ are specified: `Tgov1 ` | the TGOV1 governor model | `pmech`, `speed` | `Trate`, `R`, `T1`, `T2`, `T3`, `Pvmax`, `Pvmin`, `Dt` | `none` `Regca` | WECC REGCA renewable generator/converter model | `bus`, `ipcmd`\*, `iqcmd`\*, `ibranchr`\*, `ibranchi`\*, `pbranch`\*, `qbranch`\* | `p0`, `q0`, `mva`, `Tg`, `TM`, `Rqmax`, `Rqmin`, `Rpmax`, `sL`, `IL1`, `VL0`, `VL1`, `VA0`, `VA1`, `Vhvmax`, `Qmin`, `Khv`, `Xe` | `ir`, `ii`, `p`, `q` `Ieeet1` | the IEEET1 exciter model | `bus`, `speed`, `efd`, `vs`\* | `Tr`, `Ka`, `Ta`, `Ke`, `Te`, `Kf`, `Tf`, `Vrmin`, `Vrmax`, `E1`, `E2`, `Se1`, `Se2`, `Ispdlim` | `efd`, `ksat` + `Esdc1a` | the ESDC1A exciter model | `bus`, `efd`, `speed`\*, `vref`\*, `vs`\*, `vuel`\* | `Tr`, `Ka`, `Ta`, `Tb`, `Tc`, `Vrmax`, `Vrmin`, `Ke`, `Te`, `Kf`, `Tf1`, `Spdmlt`, `E1`, `Se1`, `E2`, `Se2`, `UEL`, `exclim` | `efd`, `vc`, `vr`, `vf`, `se`, `vfe` `SexsPti` | the SEXS-PTI simplified exciter model | `bus`, `efd`, `vs`\* | `Ta`, `Tb`, `Te`, `K`, `Efdmax`, `Efdmin` | `efd` `Ieeest` | the IEEEST stabilizer model | `input`, `output` | `A1`, `A2`, `A3`, `A4`, `A5`, `A6`, `T1`, `T2`, `T3`, `T4`, `T5`, `T6`, `Ks`, `Lsmin`, `Lsmax`, `Vcl`, `Vcu`, `Tdelay` | `vss` `BusFault` | simple impedance-based fault at a bus | `bus`, `status`\* | `state0`, `R`, `X` | `state`, `ir`, `ii` diff --git a/GridKit/Model/PhasorDynamics/SystemModelData.hpp b/GridKit/Model/PhasorDynamics/SystemModelData.hpp index 2ff9adea8..699d750e5 100644 --- a/GridKit/Model/PhasorDynamics/SystemModelData.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModelData.hpp @@ -11,6 +11,7 @@ #include #include #include +#include #include #include #include @@ -44,6 +45,7 @@ namespace GridKit using BusFaultDataT = BusFaultData; using RegcaDataT = Converter::RegcaData; using Tgov1DataT = Governor::Tgov1Data; + using Esdc1aDataT = Exciter::Esdc1aData; using Ieeet1DataT = Exciter::Ieeet1Data; using SexsPtiDataT = Exciter::SexsPtiData; using IeeestDataT = Stabilizer::IeeestData; @@ -104,6 +106,7 @@ namespace GridKit std::vector loadz; ///< LoadZ instances within the model std::vector loadzip; ///< LoadZIP instances within the model std::vector gov; ///< Governors within the model + std::vector esdc1a; ///< ESDC1A exciters within the model std::vector exciter; ///< Exciters within the model std::vector sexspti; ///< SEXS-PTI exciters within the model std::vector stabilizer; ///< Stabilizers within the model diff --git a/GridKit/Model/PhasorDynamics/SystemModelDataJSONParser.hpp b/GridKit/Model/PhasorDynamics/SystemModelDataJSONParser.hpp index 12492f995..8c564e07c 100644 --- a/GridKit/Model/PhasorDynamics/SystemModelDataJSONParser.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModelDataJSONParser.hpp @@ -151,6 +151,12 @@ namespace GridKit raw_component.get_to(exciter); sm.exciter.push_back(exciter); } + else if (kind == "Esdc1a") + { + typename SystemModelData::Esdc1aDataT exciter; + raw_component.get_to(exciter); + sm.esdc1a.push_back(exciter); + } else if (kind == "SexsPti") { typename SystemModelData::SexsPtiDataT exciter; diff --git a/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp b/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp index 0f22cacb5..d4c126cd7 100644 --- a/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp @@ -349,6 +349,54 @@ namespace GridKit addComponent(exciter); } + for (const auto& excitedata : data.esdc1a) + { + BusT* bus = nullptr; + if (excitedata.buses.contains(Esdc1aBuses::bus)) + { + bus = getBus(excitedata.buses.at(Esdc1aBuses::bus)); + } + + auto* exciter = new Esdc1a(bus, excitedata); + + if (excitedata.signal_inputs.contains(Esdc1aSignalInputs::speed)) + { + IdxT speed = excitedata.signal_inputs.at(Esdc1aSignalInputs::speed); + constexpr auto OMEGA = Esdc1aExternalVariables::OMEGA; + exciter->getSignals().template attachSignalNode(getSignal(speed)); + } + + if (excitedata.signal_inputs.contains(Esdc1aSignalInputs::vref)) + { + IdxT vref = excitedata.signal_inputs.at(Esdc1aSignalInputs::vref); + constexpr auto VREF = Esdc1aExternalVariables::VREF; + exciter->getSignals().template attachSignalNode(getSignal(vref)); + } + + if (excitedata.signal_outputs.contains(Esdc1aSignalOutputs::efd)) + { + IdxT efd = excitedata.signal_outputs.at(Esdc1aSignalOutputs::efd); + constexpr auto EFD = Esdc1aInternalVariables::EFD; + exciter->getSignals().template assignSignalNode(getSignal(efd)); + } + + if (excitedata.signal_inputs.contains(Esdc1aSignalInputs::vs)) + { + IdxT vs = excitedata.signal_inputs.at(Esdc1aSignalInputs::vs); + constexpr auto VS = Esdc1aExternalVariables::VS; + exciter->getSignals().template attachSignalNode(getSignal(vs)); + } + + if (excitedata.signal_inputs.contains(Esdc1aSignalInputs::vuel)) + { + IdxT vuel = excitedata.signal_inputs.at(Esdc1aSignalInputs::vuel); + constexpr auto VUEL = Esdc1aExternalVariables::VUEL; + exciter->getSignals().template attachSignalNode(getSignal(vuel)); + } + + addComponent(exciter); + } + for (const auto& excitedata : data.sexspti) { IdxT bus_index = 0; diff --git a/docs/Figures/PhasorDynamics/ESDC1A/diagram.png b/docs/Figures/PhasorDynamics/ESDC1A/diagram.png new file mode 100644 index 000000000..c77a100a6 Binary files /dev/null and b/docs/Figures/PhasorDynamics/ESDC1A/diagram.png differ diff --git a/docs/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/README.md b/docs/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/README.md new file mode 100644 index 000000000..aa9a4c776 --- /dev/null +++ b/docs/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/README.md @@ -0,0 +1,6 @@ +# ESDC1A + +```{include} ../../../../../../GridKit/Model/PhasorDynamics/Exciter/ESDC1A/README.md +:start-line: 1 +:relative-images: +``` diff --git a/docs/GridKit/Model/PhasorDynamics/Exciter/README.md b/docs/GridKit/Model/PhasorDynamics/Exciter/README.md index 2980fd53c..143bad900 100644 --- a/docs/GridKit/Model/PhasorDynamics/Exciter/README.md +++ b/docs/GridKit/Model/PhasorDynamics/Exciter/README.md @@ -8,6 +8,7 @@ ESAC6A IEEET1 EXDC1 +ESDC1A ESDC2A EXAC1 ESST4B diff --git a/tests/UnitTests/PhasorDynamics/CMakeLists.txt b/tests/UnitTests/PhasorDynamics/CMakeLists.txt index d643637a4..87950ac62 100644 --- a/tests/UnitTests/PhasorDynamics/CMakeLists.txt +++ b/tests/UnitTests/PhasorDynamics/CMakeLists.txt @@ -96,6 +96,16 @@ target_link_libraries( GridKit::phasor_dynamics_components_dependency_tracking GridKit::testing) +add_executable(test_phasor_exciter_esdc1a runExciterEsdc1aTests.cpp) +target_link_libraries( + test_phasor_exciter_esdc1a + GridKit::definitions + GridKit::phasor_dynamics_exciter_esdc1a + GridKit::phasor_dynamics_exciter_esdc1a_dependency_tracking + GridKit::phasor_dynamics_bus + GridKit::phasor_dynamics_bus_dependency_tracking + GridKit::testing) + add_executable(test_phasor_exciter_sexspti runExciterSexsPtiTests.cpp) target_link_libraries( test_phasor_exciter_sexspti @@ -150,6 +160,13 @@ target_link_libraries( GridKit::phasor_dynamics_systemmodel_dependency_tracking GridKit::testing) +add_executable(test_phasor_component_connection runComponentConnectionTests.cpp) +target_link_libraries( + test_phasor_component_connection + GridKit::definitions + GridKit::phasor_dynamics_systemmodel + GridKit::testing) + add_test(NAME PhasorDynamicsBusTest COMMAND test_phasor_bus) add_test(NAME PhasorDynamicsBusFaultTest COMMAND test_phasor_bus_fault) add_test(NAME PhasorDynamicsBusToSignalAdapterTest COMMAND test_phasor_bustosignaladapter) @@ -157,6 +174,7 @@ add_test(NAME PhasorDynamicsBranchTest COMMAND test_phasor_branch) add_test(NAME PhasorDynamicsGenrouTest COMMAND test_phasor_genrou) add_test(NAME PhasorDynamicsGovernorTgov1Test COMMAND test_phasor_governor_tgov1) add_test(NAME PhasorDynamicsExciterIeeet1Test COMMAND test_phasor_exciter_ieeet1) +add_test(NAME PhasorDynamicsExciterEsdc1aTest COMMAND test_phasor_exciter_esdc1a) add_test(NAME PhasorDynamicsGensalTest COMMAND test_phasor_gensal) add_test(NAME PhasorDynamicsExciterSexsPtiTest COMMAND test_phasor_exciter_sexspti) add_test(NAME PhasorDynamicsConverterRegcaTest COMMAND test_phasor_converter_regca) @@ -169,6 +187,7 @@ add_test( COMMAND test_phasor_system WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}) add_test(NAME PhasorDynamicsSystemSingleComponentTest COMMAND test_phasor_system_single_component) +add_test(NAME PhasorDynamicsComponentConnectionTest COMMAND test_phasor_component_connection) install( TARGETS test_phasor_bus @@ -180,6 +199,7 @@ install( test_phasor_genrou test_phasor_governor_tgov1 test_phasor_exciter_ieeet1 + test_phasor_exciter_esdc1a test_phasor_gensal test_phasor_exciter_sexspti test_phasor_converter_regca @@ -187,4 +207,5 @@ install( test_phasor_gen_classical test_phasor_system test_phasor_system_single_component + test_phasor_component_connection RUNTIME DESTINATION bin) diff --git a/tests/UnitTests/PhasorDynamics/ComponentConnectionTests.hpp b/tests/UnitTests/PhasorDynamics/ComponentConnectionTests.hpp new file mode 100644 index 000000000..b973b830e --- /dev/null +++ b/tests/UnitTests/PhasorDynamics/ComponentConnectionTests.hpp @@ -0,0 +1,86 @@ +#pragma once + +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace GridKit +{ + namespace Testing + { + /// Connection tests for pairs of components that share a signal node. + /// Each case checks that the node links both components and that they + /// agree on its value after initialization. Solver-driven cases live + /// in @ref PDIntegrationTests. + template + class ComponentConnectionTests + { + public: + using ScalarT = scalar_type; + using IdxT = index_type; + using RealT = typename PhasorDynamics::Component::RealT; + + ComponentConnectionTests() = default; + ~ComponentConnectionTests() = default; + + // The tolerance only absorbs floating-point roundoff. + static constexpr RealT kTol = std::numeric_limits::epsilon(); + + /// GENROU initializes first and writes the field voltage it needs to + /// the shared node. ESDC1A then initializes around that value and + /// must leave it unchanged at a steady state. + TestOutcome genrouEsdc1a() + { + using MachineExternal = PhasorDynamics::GenrouExternalVariables; + using ExciterInternal = PhasorDynamics::Exciter::Esdc1aInternalVariables; + using ExciterParams = PhasorDynamics::Exciter::Esdc1aParameters; + + TestStatus success = true; + + PhasorDynamics::SystemModel system; + PhasorDynamics::BusInfinite bus( + static_cast(1.0), + static_cast(0.0)); + PhasorDynamics::SignalNode efd; + PhasorDynamics::Genrou machine(&bus); + + PhasorDynamics::Exciter::Esdc1aData exciter_data; + exciter_data.parameters[ExciterParams::Tr] = static_cast(0.02); + exciter_data.parameters[ExciterParams::Tb] = static_cast(0.5); + + PhasorDynamics::Exciter::Esdc1a exciter(&bus, exciter_data); + + machine.getSignals().template attachSignalNode(&efd); + exciter.getSignals().template assignSignalNode(&efd); + + system.addBus(&bus); + system.addComponent(&machine); + system.addComponent(&exciter); + + success *= system.allocate() == 0; + success *= efd.linked(); + success *= system.initialize() == 0; + success *= system.evaluateResidual() == 0; + + // At zero power the required field voltage equals the terminal voltage. + success *= isEqual(efd.read(), static_cast(1.0), kTol); + + const auto* residual = exciter.getResidual().getData(); + for (IdxT row = 0; row < exciter.size(); ++row) + { + success *= isEqual(residual[row], static_cast(0.0), kTol); + } + + return success.report(__func__); + } + }; + + } // namespace Testing +} // namespace GridKit diff --git a/tests/UnitTests/PhasorDynamics/ExciterEsdc1aTests.hpp b/tests/UnitTests/PhasorDynamics/ExciterEsdc1aTests.hpp new file mode 100644 index 000000000..881cd48ce --- /dev/null +++ b/tests/UnitTests/PhasorDynamics/ExciterEsdc1aTests.hpp @@ -0,0 +1,1523 @@ +#pragma once + +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace GridKit +{ + namespace Testing + { + using Log = ::GridKit::Utilities::Logger; + + template + class ExciterEsdc1aTests + { + public: + using ScalarT = scalar_type; + using IdxT = index_type; + using RealT = typename PhasorDynamics::Component::RealT; + + ExciterEsdc1aTests() = default; + ~ExciterEsdc1aTests() = default; + + // Rows that divide by a time constant amplify by 1/T; the + // worst here is VR, which recovers Ka*(Vr/Ka) through 1/Ta at 31 eps. + static constexpr RealT kTol = + static_cast(100.0) * std::numeric_limits::epsilon(); + + /// Construction and every verify() error class, including parameter + /// types, parameter relationships, bus ownership, and signal linkage. + TestOutcome validation() + { + TestStatus success = true; + + noteExpectedLogs("Testing ESDC1A defaults and invalid configurations. " + "Logged errors and time-constant warnings are expected."); + + PhasorDynamics::Bus bus(1.0, 0.0); + + PhasorDynamics::Exciter::Esdc1a empty(&bus); + success *= (empty.size() == static_cast(Internal::MAXIMUM)); + success *= (empty.getMonitor() == nullptr); + success *= (empty.verify() > 0); + + Fixture configured(makeData()); + success *= (configured.esdc1a.size() == static_cast(Internal::MAXIMUM)); + success *= (configured.esdc1a.getMonitor() != nullptr); + success *= configured.prepare(1.2); + + // ESDC1A has no required parameters: an empty parameter set is the + // documented-default model. + Fixture minimal(makeMinimalData()); + success *= minimal.prepare(1.2); + success *= defaultsMatchDocumentedValues(); + + // A model without the required efd output assignment is rejected. + PhasorDynamics::Exciter::Esdc1a unassigned(&bus, makeData()); + success *= (unassigned.verify() > 0); + + success *= invalidParameterCase(Params::Ka, 0.0); + success *= invalidParameterCase(Params::Ta, -0.1); + success *= invalidParameterCase(Params::Te, -0.1); + success *= invalidParameterCase(Params::Tc, -0.1); + success *= invalidParameterCase(Params::Tr, -0.1); + success *= invalidParameterCase(Params::Tb, -0.1); + success *= invalidParameterCase(Params::Tf1, -0.1); + success *= invalidParameterCase(Params::Vrmin, 2.0); + success *= invalidParameterCase(Params::UEL, static_cast(4)); + success *= invalidParameterCase(Params::UEL, static_cast(2.0)); + success *= invalidParameterCase(Params::UEL, static_cast(2.5)); + success *= invalidParameterCase(Params::UEL, true); + success *= invalidParameterCase(Params::Se1, 0.0); + success *= invalidParameterCase(Params::E2, 2.8); + success *= invalidParameterCase(Params::Se2, 0.08); + success *= invalidParameterCase(Params::E1, -1.0); + + for (const Params parameter : {Params::Tr, + Params::Ka, + Params::Ta, + Params::Tb, + Params::Tc, + Params::Vrmax, + Params::Vrmin, + Params::Ke, + Params::Te, + Params::Kf, + Params::Tf1, + Params::E1, + Params::Se1, + Params::E2, + Params::Se2}) + { + success *= invalidParameterCase(parameter, std::numeric_limits::quiet_NaN()); + success *= invalidParameterCase(parameter, std::numeric_limits::infinity()); + } + + // Saturation voltage and coefficient pairs must move in the same + // direction; either enumeration direction is otherwise valid. + auto reversed_saturation = makeData(); + reversed_saturation.parameters[Params::E1] = 3.7; + reversed_saturation.parameters[Params::Se1] = 0.33; + reversed_saturation.parameters[Params::E2] = 2.8; + reversed_saturation.parameters[Params::Se2] = 0.08; + Fixture reversed_saturation_fixture(reversed_saturation); + success *= (reversed_saturation_fixture.esdc1a.verify() == 0); + + auto crossed_ascending = makeData(); + crossed_ascending.parameters[Params::Se1] = 0.33; + crossed_ascending.parameters[Params::Se2] = 0.08; + Fixture crossed_ascending_fixture(crossed_ascending); + success *= (crossed_ascending_fixture.esdc1a.verify() > 0); + + auto crossed_descending = makeData(); + crossed_descending.parameters[Params::E1] = 3.7; + crossed_descending.parameters[Params::E2] = 2.8; + Fixture crossed_descending_fixture(crossed_descending); + success *= (crossed_descending_fixture.esdc1a.verify() > 0); + + // Integer JSON values are accepted for real parameters; booleans are + // not numeric. + auto integer_real = makeData(); + integer_real.parameters[Params::Ka] = static_cast(40); + Fixture integer_real_fixture(integer_real); + success *= (integer_real_fixture.esdc1a.verify() == 0); + success *= invalidParameterCase(Params::Ka, true); + + // Binary selectors accept JSON booleans only. + auto boolean_switches = makeData(); + boolean_switches.parameters[Params::Spdmlt] = true; + boolean_switches.parameters[Params::exclim] = false; + Fixture boolean_switch_fixture(boolean_switches); + boolean_switch_fixture.attachAllInputs(); + success *= (boolean_switch_fixture.esdc1a.verify() == 0); + + for (const Params flag : {Params::Spdmlt, Params::exclim}) + { + success *= invalidParameterCase(flag, static_cast(0)); + success *= invalidParameterCase(flag, static_cast(1)); + success *= invalidParameterCase(flag, static_cast(2)); + success *= invalidParameterCase(flag, static_cast(0.0)); + success *= invalidParameterCase(flag, static_cast(0.5)); + success *= invalidParameterCase(flag, static_cast(1.0)); + } + + // The enabled speed multiplier requires an attached speed input. + auto speed_required = makeData(); + speed_required.parameters[Params::Spdmlt] = true; + Fixture speed_required_fixture(speed_required); + success *= (speed_required_fixture.esdc1a.verify() > 0); + + PhasorDynamics::SignalNode busless_efd_node; + PhasorDynamics::Exciter::Esdc1a busless(nullptr, makeData()); + busless.getSignals().template assignSignalNode(&busless_efd_node); + success *= (busless.verify() > 0); + + success *= unlinkedSignalRejected(); + success *= unlinkedSignalRejected(); + success *= unlinkedSignalRejected(); + success *= unlinkedSignalRejected(); + + // All five floored time constants at zero use the documented + // numerical floor and still admit a consistent steady-state + // initialization. + auto zero_time = makeData(); + zero_time.parameters[Params::Tr] = 0.0; + zero_time.parameters[Params::Ta] = 0.0; + zero_time.parameters[Params::Tb] = 0.0; + zero_time.parameters[Params::Te] = 0.0; + zero_time.parameters[Params::Tf1] = 0.0; + + Fixture floored(zero_time); + success *= floored.initialize(1.2); + success *= (floored.evaluate() == 0); + success *= allResidualsZero(floored.esdc1a); + + return success.report(__func__); + } + + /// Initialization with every port attached: the seeded field voltage + /// is preserved, the resolved voltage reference is published, the + /// known inputs are read but never overwritten, and every selector + /// combination reaches a zero-derivative, zero-residual state. + TestOutcome initializationAndSignals() + { + TestStatus success = true; + + Fixture fixture(makeData()); + fixture.attachAllInputs(99.0); + fixture.input(External::OMEGA) = 0.02; + fixture.input(External::VS) = 0.03; + fixture.input(External::VUEL) = -0.4; + success *= fixture.initialize(1.2); + success *= (fixture.esdc1a.tagDifferentiable() == 0); + success *= (fixture.evaluate() == 0); + + const auto* y = fixture.esdc1a.y().getData(); + success *= scalarMatches(y[static_cast(Internal::EFDP)], 1.2, "EFDP"); + success *= scalarMatches(y[static_cast(Internal::VC)], 1.0, "VC"); + success *= scalarMatches(y[static_cast(Internal::VR)], 0.12, "VR"); + success *= scalarMatches(y[static_cast(Internal::VF)], 0.0, "VF"); + success *= scalarMatches(y[static_cast(Internal::EV)], 0.003, "EV gate input"); + success *= scalarMatches(y[static_cast(Internal::VHV)], 0.003, "VHV"); + success *= scalarMatches(y[static_cast(Internal::SE)], 0.0, "SE"); + success *= scalarMatches(y[static_cast(Internal::VFE)], 0.12, "VFE"); + success *= scalarMatches(fixture.efd(), 1.2, "seeded efd"); + + success *= scalarMatches(fixture.input(External::VREF), 0.973, "published vref"); + success *= scalarMatches(fixture.input(External::OMEGA), 0.02, "preserved speed input"); + success *= scalarMatches(fixture.input(External::VS), 0.03, "preserved vs input"); + success *= scalarMatches(fixture.input(External::VUEL), -0.4, "preserved vuel input"); + + RealT time = 0.0; + Model::VariableMonitorController monitor(time); + monitor.addMonitor(fixture.esdc1a.getMonitor()); + std::stringstream monitor_output; + monitor.addSink({Model::VariableMonitorFormat::CSV}, monitor_output); + monitor.start(); + monitor.print(); + monitor.stop(); + + std::string monitor_header; + std::string monitor_values; + std::getline(monitor_output, monitor_header); + std::getline(monitor_output, monitor_values); + success *= (monitor_header == "t,Esdc1a_esdc1a_test_efd,Esdc1a_esdc1a_test_vc," + "Esdc1a_esdc1a_test_vr,Esdc1a_esdc1a_test_vf," + "Esdc1a_esdc1a_test_se,Esdc1a_esdc1a_test_vfe"); + const auto monitored = Tokenizer(monitor_values, ',')(); + if (monitored.size() == 7) + { + success *= scalarMatches(monitored[1], 1.2, "monitored efd"); + success *= scalarMatches(monitored[2], 1.0, "monitored vc"); + success *= scalarMatches(monitored[3], 0.12, "monitored vr"); + success *= scalarMatches(monitored[4], 0.0, "monitored vf"); + success *= scalarMatches(monitored[5], 0.0, "monitored se"); + success *= scalarMatches(monitored[6], 0.12, "monitored vfe"); + } + else + { + std::cout << "ESDC1A monitor emitted " << monitored.size() + << " values instead of 7\n"; + success = false; + } + + for (size_t i = 0; i < static_cast(fixture.esdc1a.size()); ++i) + { + const bool expected = i <= static_cast(Internal::XLL); + if (fixture.esdc1a.tag()[i] != expected) + { + std::cout << "ESDC1A differentiability tag " << i << " mismatch\n"; + success = false; + } + } + success *= allResidualsZero(fixture.esdc1a); + + // With no attached inputs the latched values act as constant + // references. + Fixture unattached(makeData()); + success *= unattached.initialize(1.2); + success *= (unattached.evaluate() == 0); + success *= allResidualsZero(unattached.esdc1a); + + // Every selector combination must preserve the seeded field voltage + // and produce a zero-derivative, zero-residual state. + for (IdxT uel = 0; uel < 4; ++uel) + { + for (const bool speed_flag : {false, true}) + { + for (const bool limit_flag : {false, true}) + { + auto scenario_data = makeData(); + scenario_data.parameters[Params::UEL] = uel; + scenario_data.parameters[Params::Spdmlt] = speed_flag; + scenario_data.parameters[Params::exclim] = limit_flag; + + Fixture scenario(scenario_data); + scenario.attachAllInputs(); + scenario.input(External::OMEGA) = 0.02; + scenario.input(External::VS) = 0.03; + scenario.input(External::VUEL) = -0.4; + if (!scenario.initialize(1.2)) + { + std::cout << "ESDC1A initialization scenario failed: UEL=" << uel + << ", Spdmlt=" << speed_flag + << ", exclim=" << limit_flag << '\n'; + success = false; + continue; + } + + success *= (scenario.evaluate() == 0); + success *= allResidualsZero(scenario.esdc1a); + success *= scalarMatches(scenario.efd(), 1.2, "scenario efd preservation"); + } + } + } + + return success.report(__func__); + } + + /// The inadmissible initialization points: every rejection is atomic, + /// and the admissible operating points next to them still initialize + /// to zero residuals. + TestOutcome initializationDomain() + { + TestStatus success = true; + + noteExpectedLogs("Testing inadmissible ESDC1A initialization points. " + "Logged errors are expected."); + + // An enabled speed multiplier must remain finite and strictly + // positive. Other initialization limits are relaxed so these cases + // isolate that domain. + auto speed_data = makeData(); + speed_data.parameters[Params::Spdmlt] = true; + speed_data.parameters[Params::exclim] = false; + speed_data.parameters[Params::Vrmin] = -100.0; + speed_data.parameters[Params::Vrmax] = 100.0; + success *= initializationRejectedAtomically(speed_data, + 1.2, + {{External::OMEGA, -1.0}, + {External::VREF, 77.0}, + {External::VS, 77.0}, + {External::VUEL, -77.0}}, + "zero speed multiplier"); + success *= initializationRejectedAtomically(speed_data, + 1.2, + {{External::OMEGA, -1.1}, + {External::VREF, 77.0}, + {External::VS, 77.0}, + {External::VUEL, -77.0}}, + "negative speed multiplier"); + + // The enabled exciter lower limit rejects a negative field-voltage + // state before initialization writes any storage. + success *= initializationRejectedAtomically(makeData(), + -0.2, + {{External::OMEGA, 0.0}, + {External::VREF, 77.0}, + {External::VS, 77.0}, + {External::VUEL, -0.5}}, + "field-voltage state below zero limit"); + + // The seeded field voltage maps to a regulator output above Vrmax. + auto limit_data = makeData(); + limit_data.parameters[Params::Vrmax] = 0.05; + limit_data.parameters[Params::Vrmin] = -0.05; + success *= initializationRejectedAtomically(limit_data, + 1.2, + {{External::OMEGA, 0.0}, + {External::VREF, 77.0}, + {External::VS, 77.0}, + {External::VUEL, -77.0}}, + "regulator output outside limits"); + + // A UEL input above the gate operating point holds the high-value + // gate active, which the smooth gate cannot represent at rest. + success *= initializationRejectedAtomically(makeData(), + 1.2, + {{External::OMEGA, 0.0}, + {External::VREF, 77.0}, + {External::VS, 77.0}, + {External::VUEL, 0.5}}, + "active high-value gate"); + + // A non-finite field-voltage seed is rejected before any signal is + // published. + Fixture nonfinite(makeData()); + nonfinite.attachAllInputs(77.0); + success *= nonfinite.prepare(std::numeric_limits::quiet_NaN()); + success *= (nonfinite.esdc1a.initialize() != 0); + success *= scalarMatches(nonfinite.input(External::VREF), 77.0, "rejected vref preservation"); + + success *= initializationRejectedAtomically( + speed_data, + 1.2, + {{External::OMEGA, std::numeric_limits::infinity()}, + {External::VREF, 77.0}, + {External::VS, 0.0}, + {External::VUEL, -77.0}}, + "non-finite speed input"); + + success *= initializationRejectedAtomically( + makeData(), + 1.2, + {{External::OMEGA, 0.0}, + {External::VREF, 77.0}, + {External::VS, std::numeric_limits::infinity()}, + {External::VUEL, -0.5}}, + "non-finite stabilizer input"); + + success *= initializationRejectedAtomically( + makeData(), + 1.2, + {{External::OMEGA, 0.0}, + {External::VREF, 77.0}, + {External::VS, 0.0}, + {External::VUEL, -std::numeric_limits::infinity()}}, + "non-finite UEL input"); + + success *= initializationRejectedAtomically( + makeData(), + 1.2, + {{External::OMEGA, 0.0}, + {External::VREF, 77.0}, + {External::VS, 0.0}, + {External::VUEL, -0.5}}, + "zero terminal voltage", + 0.0, + 0.0); + + // An invalid configuration is rejected before any state is written. + auto invalid_data = makeData(); + invalid_data.parameters[Params::Ka] = 0.0; + Fixture invalid_fixture(invalid_data); + invalid_fixture.attachAllInputs(); + success *= (invalid_fixture.esdc1a.allocate() == 0); + poisonState(invalid_fixture, 1.2); + const auto invalid_y = copyVector(invalid_fixture.esdc1a.y()); + const auto invalid_yp = copyVector(invalid_fixture.esdc1a.yp()); + if (invalid_fixture.esdc1a.initialize() == 0) + { + std::cout << "Expected initialization rejection: invalid configuration\n"; + success = false; + } + success *= vectorUnchanged(invalid_fixture.esdc1a.y(), invalid_y, "state"); + success *= vectorUnchanged(invalid_fixture.esdc1a.yp(), invalid_yp, "derivative"); + + // The zero boundary is admissible, and disabling the limiter admits + // the same negative seed rejected above. + Fixture zero_boundary(makeData()); + zero_boundary.attachAllInputs(); + zero_boundary.input(External::VUEL) = -0.5; + success *= zero_boundary.initialize(0.0); + success *= (zero_boundary.evaluate() == 0); + success *= allResidualsZero(zero_boundary.esdc1a); + + auto unlimited_data = makeData(); + unlimited_data.parameters[Params::exclim] = false; + Fixture unlimited(unlimited_data); + unlimited.attachAllInputs(); + unlimited.input(External::VUEL) = -0.5; + success *= unlimited.initialize(-0.2); + success *= (unlimited.evaluate() == 0); + success *= allResidualsZero(unlimited.esdc1a); + + // A depressed speed input rescales the seed without rejection. + auto admissible_speed = makeData(); + admissible_speed.parameters[Params::Spdmlt] = true; + Fixture speed_fixture(admissible_speed); + speed_fixture.attachAllInputs(); + speed_fixture.input(External::OMEGA) = -0.5; + success *= speed_fixture.initialize(1.2); + success *= (speed_fixture.evaluate() == 0); + success *= allResidualsZero(speed_fixture.esdc1a); + success *= scalarMatches( + speed_fixture.esdc1a.y().getData()[static_cast(Internal::EFDP)], + 2.4, + "rescaled EFDP"); + + // The gate stays representable arbitrarily close to its operating + // point: the ramp inverse seeds a 0.003 margin exactly. + Fixture near_gate(makeData()); + near_gate.attachAllInputs(); + success *= near_gate.initialize(1.2); + success *= (near_gate.evaluate() == 0); + success *= allResidualsZero(near_gate.esdc1a); + + // Summing-junction routing removes the gate constraint entirely. + auto junction_data = makeData(); + junction_data.parameters[Params::UEL] = static_cast(2); + Fixture junction(junction_data); + junction.attachAllInputs(); + junction.input(External::VUEL) = 0.7; + success *= junction.initialize(1.2); + success *= (junction.evaluate() == 0); + success *= allResidualsZero(junction.esdc1a); + + return success.report(__func__); + } + + /// A fixed numerical answer key for all 11 ESDC1A equations. The + /// expected values are literals, not a second implementation of ESDC1A. + TestOutcome residualEquations() + { + TestStatus success = true; + + Fixture fixture(makeResidualData(), kStateVr, kStateVi); + fixture.attachAllInputs(); + success *= fixture.initialize(1.2); + setAnswerKeyInputs(fixture); + setAnswerKeyState(fixture.esdc1a); + success *= (fixture.evaluate() == 0); + + // Values are pinned after an independent one-time evaluation of the + // documented equations at setAnswerKeyState()/setAnswerKeyInputs(). + const std::array(Internal::MAXIMUM)> expected{{ + {Internal::EFDP, 0.04000000000000004}, + {Internal::VC, 0.19442890089805262}, + {Internal::VR, 0.13666666666666663}, + {Internal::VF, -0.022222222222222213}, + {Internal::XLL, 0.024999999999999994}, + {Internal::EV, -0.27999999999999986}, + {Internal::VLL, -0.007500000000000031}, + {Internal::VHV, 0.31535073999664665}, + {Internal::SE, -0.07541019662496842}, + {Internal::VFE, 0.1600000000000001}, + {Internal::EFD, 0.9100000000000001}, + }}; + + success *= (static_cast(fixture.esdc1a.getResidual().getSize()) == expected.size()); + success *= residualsMatch(fixture.esdc1a, expected); + + return success.report(__func__); + } + + /// The transducer, summing junction, lead-lag, stabilizing feedback, + /// and regulator anti-windup behavior at driven states with literal + /// expectations. + TestOutcome voltageRegulation() + { + TestStatus success = true; + + Fixture fixture(makeData()); + fixture.attachAllInputs(); + success *= fixture.initialize(1.2); + + // Transducer: the sensed voltage relaxes toward the bus magnitude. + setState(fixture.esdc1a, {{Internal::VC, 1.1}}); + setDerivative(fixture.esdc1a, {{Internal::VC, 0.2}}); + success *= (fixture.evaluate() == 0); + success *= residualsMatch(fixture.esdc1a, {{Internal::VC, -5.2}}, "voltage transducer"); + + // The field-voltage state and the stabilizing feedback share the + // (VR - VFE) drive. + setState(fixture.esdc1a, {{Internal::VR, 0.6}, {Internal::VFE, 0.2}, {Internal::VF, 0.1}}); + setDerivative(fixture.esdc1a, {{Internal::EFDP, 0.1}, {Internal::VF, 0.05}}); + success *= (fixture.evaluate() == 0); + success *= residualsMatch(fixture.esdc1a, + {{Internal::EFDP, 0.7}, {Internal::VF, -0.13571428571428573}}, + "field-voltage and feedback drive"); + + // Summing junction: UEL < 2 excludes the UEL input from the error. + Fixture summing(makeData()); + summing.attachAllInputs(); + success *= summing.initialize(1.2); + summing.input(External::VREF) = 1.1; + summing.input(External::VS) = 0.05; + summing.input(External::VUEL) = 0.2; + setState(summing.esdc1a, {{Internal::VC, 0.9}, {Internal::VF, 0.02}, {Internal::EV, 0.1}}); + success *= (summing.evaluate() == 0); + success *= residualsMatch(summing.esdc1a, {{Internal::EV, 0.13}}, "summing junction"); + + // UEL >= 2 routes the UEL input through the summing junction and + // turns the high-value gate into a lead-lag passthrough. + auto junction_data = makeData(); + junction_data.parameters[Params::UEL] = static_cast(2); + Fixture junction(junction_data); + junction.attachAllInputs(); + success *= junction.initialize(1.2); + junction.input(External::VREF) = 1.1; + junction.input(External::VS) = 0.05; + junction.input(External::VUEL) = 0.2; + setState(junction.esdc1a, + {{Internal::VC, 0.9}, + {Internal::VF, 0.02}, + {Internal::EV, 0.1}, + {Internal::VLL, 0.5}, + {Internal::VHV, 0.2}}); + success *= (junction.evaluate() == 0); + success *= residualsMatch(junction.esdc1a, + {{Internal::EV, 0.33}, {Internal::VHV, 0.3}}, + "summing-junction UEL routing"); + + // An active lead-lag pair advances the error and relaxes its state. + auto lead_lag_data = makeData(); + lead_lag_data.parameters[Params::Tc] = 0.2; + Fixture lead_lag(lead_lag_data); + lead_lag.attachAllInputs(); + success *= lead_lag.initialize(1.2); + setState(lead_lag.esdc1a, {{Internal::XLL, 0.4}, {Internal::EV, 0.7}, {Internal::VLL, 0.5}}); + setDerivative(lead_lag.esdc1a, {{Internal::XLL, 0.0}}); + success *= (lead_lag.evaluate() == 0); + success *= residualsMatch(lead_lag.esdc1a, + {{Internal::XLL, 0.6}, {Internal::VLL, 0.02}}, + "lead-lag"); + + // The regulator anti-windup blocks outward rates at both limits and + // admits restoring rates. + struct AntiWindupCase + { + const char* label; + RealT vr; + RealT vhv; + RealT expected; + }; + + const std::array antiwindup_cases{{ + {"Vrmax blocks an outward regulator rate", 1.5, 0.05, 0.0}, + {"Vrmax admits a restoring regulator rate", 1.5, 0.025, -5.0}, + {"Vrmin blocks an outward regulator rate", -1.5, -0.05, 0.0}, + {"Vrmin admits a restoring regulator rate", -1.5, -0.025, 5.0}, + }}; + + for (const auto& test_case : antiwindup_cases) + { + setState(fixture.esdc1a, {{Internal::VR, test_case.vr}, {Internal::VHV, test_case.vhv}}); + setDerivative(fixture.esdc1a, {{Internal::VR, 0.0}}); + success *= (fixture.evaluate() == 0); + success *= residualsMatch(fixture.esdc1a, + {{Internal::VR, test_case.expected}}, + test_case.label); + } + + return success.report(__func__); + } + + /// High-value gate selection, quadratic saturation, field-voltage-state + /// limiting, and the speed multiplier at driven states with literal + /// expectations. + TestOutcome excitationLimits() + { + TestStatus success = true; + + // The gate passes the larger of VLL and VUEL when UEL < 2; the + // smooth maximum keeps two-sided sensitivity at the tie. + struct GateCase + { + const char* label; + RealT vuel; + RealT expected; + }; + + const std::array gate_cases{{ + {"gate selects the lead-lag branch", -0.5, 0.3}, + {"gate selects the UEL branch", 0.8, 0.6000000000000001}, + {"gate tie point", 0.5, 0.30288811325233306}, + }}; + + Fixture gate(makeData()); + gate.attachAllInputs(); + success *= gate.initialize(1.2); + for (const auto& test_case : gate_cases) + { + gate.input(External::VUEL) = test_case.vuel; + setState(gate.esdc1a, {{Internal::VLL, 0.5}, {Internal::VHV, 0.2}}); + success *= (gate.evaluate() == 0); + success *= residualsMatch(gate.esdc1a, {{Internal::VHV, test_case.expected}}, test_case.label); + } + + // Both valid point orderings produce the same quadratic curve at the + // supplied points and on either side of the fitted knee. + struct SaturationOrderCase + { + const char* label; + RealT e1; + RealT se1; + RealT e2; + RealT se2; + }; + + const std::array saturation_order_cases{{ + {"ascending saturation points", 2.4, 0.1, 3.2, 0.5}, + {"descending saturation points", 3.2, 0.5, 2.4, 0.1}, + }}; + + for (const auto& test_case : saturation_order_cases) + { + auto data = makeResidualData(); + data.parameters[Params::E1] = test_case.e1; + data.parameters[Params::Se1] = test_case.se1; + data.parameters[Params::E2] = test_case.e2; + data.parameters[Params::Se2] = test_case.se2; + Fixture saturation(data); + saturation.attachAllInputs(); + success *= saturation.initialize(1.2); + + setState(saturation.esdc1a, {{Internal::EFDP, 2.4}, {Internal::SE, 0.0}}); + success *= (saturation.evaluate() == 0); + success *= residualsMatch(saturation.esdc1a, + {{Internal::SE, 0.1}}, + test_case.label); + + setState(saturation.esdc1a, {{Internal::EFDP, 3.2}}); + success *= (saturation.evaluate() == 0); + success *= residualsMatch(saturation.esdc1a, + {{Internal::SE, 0.5}}, + test_case.label); + + setState(saturation.esdc1a, {{Internal::EFDP, 2.0}, {Internal::SE, 0.05}}); + success *= (saturation.evaluate() == 0); + success *= residualsMatch(saturation.esdc1a, + {{Internal::SE, -0.035410196624968436}}, + test_case.label); + + setState(saturation.esdc1a, {{Internal::EFDP, 1.0}}); + success *= (saturation.evaluate() == 0); + success *= residualsMatch(saturation.esdc1a, + {{Internal::SE, -0.05}}, + test_case.label); + } + + auto disabled_data = makeResidualData(); + disabled_data.parameters[Params::Se1] = 0.0; + disabled_data.parameters[Params::Se2] = 0.0; + Fixture disabled(disabled_data); + disabled.attachAllInputs(); + success *= disabled.initialize(1.2); + setState(disabled.esdc1a, {{Internal::EFDP, 2.0}, {Internal::SE, 0.05}}); + success *= (disabled.evaluate() == 0); + success *= residualsMatch(disabled.esdc1a, {{Internal::SE, -0.05}}, "saturation disabled"); + + // The field-voltage-state lower limit blocks outward motion, admits + // restoring motion, and preserves the CommonMath transition at zero. + struct FieldLimitCase + { + const char* label; + bool enabled; + RealT efdp; + RealT vr; + RealT expected; + }; + + const std::array field_limit_cases{{ + {"below zero blocks an outward field-voltage rate", true, -0.2, -0.1, 0.0}, + {"zero retains the smooth lower-limit transition", true, 0.0, -0.1, -0.1}, + {"above zero admits a downward field-voltage rate", true, 0.2, -0.1, -0.2}, + {"below zero admits a restoring field-voltage rate", true, -0.2, 0.1, 0.2}, + {"disabled lower limit admits an outward rate", false, -0.2, -0.1, -0.2}, + }}; + + for (const auto& test_case : field_limit_cases) + { + auto data = makeData(); + data.parameters[Params::exclim] = test_case.enabled; + Fixture limit(data); + limit.attachAllInputs(); + success *= limit.initialize(1.2); + setState(limit.esdc1a, + {{Internal::EFDP, test_case.efdp}, + {Internal::VR, test_case.vr}, + {Internal::VFE, 0.0}}); + setDerivative(limit.esdc1a, {{Internal::EFDP, 0.0}}); + success *= (limit.evaluate() == 0); + success *= residualsMatch(limit.esdc1a, + {{Internal::EFDP, test_case.expected}}, + test_case.label); + } + + // The lower-limit selector does not alter the algebraic exciter + // feedback drive. + auto feedback_data = makeData(); + feedback_data.parameters[Params::Ke] = -0.2; + feedback_data.parameters[Params::Se1] = 0.0; + feedback_data.parameters[Params::Se2] = 0.0; + for (const bool enabled : {false, true}) + { + auto data = feedback_data; + data.parameters[Params::exclim] = enabled; + Fixture feedback(data); + feedback.attachAllInputs(); + feedback.input(External::VUEL) = -0.5; + success *= feedback.initialize(1.2); + setState(feedback.esdc1a, + {{Internal::EFDP, 1.0}, {Internal::SE, 0.0}, {Internal::VFE, 0.0}}); + success *= (feedback.evaluate() == 0); + success *= residualsMatch(feedback.esdc1a, + {{Internal::VFE, -0.2}}, + enabled ? "feedback with lower limit enabled" + : "feedback with lower limit disabled"); + } + + // At the lower-limit transition, pin the assembled alpha = 1 + // field-voltage-state row independently of either Jacobian backend. + { + using DepVar = DependencyTracking::Variable; + + Fixture transition(makeData()); + transition.attachAllInputs(); + success *= transition.initialize(1.2); + setState(transition.esdc1a, + {{Internal::EFDP, 0.0}, {Internal::VR, -0.1}, {Internal::VFE, 0.0}}); + setDerivative(transition.esdc1a, {{Internal::EFDP, 0.0}}); + numberVariables(transition); + success *= (transition.evaluate() == 0); + + const auto& dependencies = + transition.esdc1a.getResidual().getData()[static_cast(Internal::EFDP)].getDependencies(); + const DepVar::DependencyMap expected{{ + {static_cast(Internal::EFDP), -13.0}, + {static_cast(Internal::VR), 1.0}, + {static_cast(Internal::VFE), -1.0}, + }}; + success *= isEqual(dependencies, expected, kTol); + } + + // The speed multiplier scales the published field voltage only when + // enabled. + for (const auto& [enabled, expected] : std::array, 2>{{ + {false, 0.0}, + {true, 0.06000000000000005}, + }}) + { + auto data = makeData(); + data.parameters[Params::Spdmlt] = enabled; + Fixture speed(data); + speed.attachAllInputs(); + success *= speed.initialize(1.2); + speed.input(External::OMEGA) = 0.05; + setState(speed.esdc1a, {{Internal::EFDP, 1.2}, {Internal::EFD, 1.2}}); + success *= (speed.evaluate() == 0); + success *= residualsMatch(speed.esdc1a, + {{Internal::EFD, expected}}, + enabled ? "speed multiplier enabled" + : "speed multiplier disabled"); + } + + return success.report(__func__); + } + +#ifdef GRIDKIT_ENABLE_ENZYME + /// A single rich state and all four external inputs drive both + /// sensitivity paths; every Enzyme CSR row must match dependency + /// tracking. + TestOutcome jacobian() + { + TestStatus success = true; + + const auto data = makeResidualData(); + + const auto dependency_jacobian = dependencyTrackingJacobian(data, success); + const auto enzyme_jacobian = enzymeJacobian(data, success); + + success *= (dependency_jacobian.size() == enzyme_jacobian.size()); + const auto rows = std::min(dependency_jacobian.size(), enzyme_jacobian.size()); + for (size_t row = 0; row < rows; ++row) + { + if (!isEqual(dependency_jacobian[row], enzyme_jacobian[row], kTol)) + { + std::cout << "ESDC1A Jacobian row " << row + << " mismatch between dependency tracking and Enzyme\n"; + success = false; + } + } + + return success.report(__func__); + } +#endif + + private: + using Esdc1aT = PhasorDynamics::Exciter::Esdc1a; + using Data = typename Esdc1aT::ModelDataT; + using Params = typename Data::Parameters; + using Mon = typename Data::MonitorableVariables; + using Internal = typename Esdc1aT::InternalVariablesT; + using External = typename Esdc1aT::ExternalVariablesT; + + using InternalRow = std::pair; + using InternalRows = std::vector; + using ExternalRow = std::pair; + using ExternalRows = std::vector; + + /// Owns the terminal bus, ESDC1A, the assigned field-voltage node, and + /// the attached input nodes. Signal storage is declared before the + /// model so every referenced node outlives ESDC1A. Copying would + /// invalidate the model and signal-node pointers. + template + class Fixture + { + private: + std::array(External::MAXIMUM)> input_values_{}; + std::array(External::MAXIMUM)> input_indices_{}; + std::array, + static_cast(External::MAXIMUM)> + input_nodes_{}; + + PhasorDynamics::SignalNode efd_node_; + + public: + explicit Fixture(const Data& data, RealT vr = 0.8, RealT vi = 0.6) + : bus(static_cast(vr), static_cast(vi)), + esdc1a(&bus, data) + { + esdc1a.getSignals().template assignSignalNode(&efd_node_); + } + + Fixture(const Fixture&) = delete; + Fixture& operator=(const Fixture&) = delete; + + /// Attach fixture-owned storage to every external input. + void attachAllInputs(RealT initial_value = 0.0) + { + const IdxT external_index_base = esdc1a.size() + bus.size(); + + for (size_t port = 0; port < input_values_.size(); ++port) + { + input_values_[port] = static_cast(initial_value); + input_indices_[port] = external_index_base + static_cast(port); + input_nodes_[port].set(&input_values_[port], &input_indices_[port]); + } + + auto& signals = esdc1a.getSignals(); + signals.template attachSignalNode( + &input_nodes_[static_cast(External::OMEGA)]); + signals.template attachSignalNode( + &input_nodes_[static_cast(External::VREF)]); + signals.template attachSignalNode( + &input_nodes_[static_cast(External::VS)]); + signals.template attachSignalNode( + &input_nodes_[static_cast(External::VUEL)]); + } + + /// Seed the assigned field-voltage node. + void seedEfd(RealT efd) + { + efd_node_.init(static_cast(efd)); + } + + /// Everything ESDC1A initialization requires: allocation, + /// verification, an initialized terminal bus, and a seeded + /// field-voltage node. + bool prepare(RealT efd) + { + const bool success = (bus.allocate() == 0) && (esdc1a.allocate() == 0) + && (esdc1a.verify() == 0) && (bus.initialize() == 0); + if (!success) + { + std::cout << "ESDC1A fixture preparation failed\n"; + return false; + } + + seedEfd(efd); + return true; + } + + /// prepare() plus successful ESDC1A initialization. + bool initialize(RealT efd) + { + if (!prepare(efd)) + { + return false; + } + if (esdc1a.initialize() != 0) + { + std::cout << "ESDC1A initialization failed\n"; + return false; + } + return true; + } + + int evaluate() + { + return esdc1a.evaluateResidual(); + } + + T efd() const + { + return efd_node_.read(); + } + + T& input(External port) + { + return input_values_[static_cast(port)]; + } + + IdxT inputIndex(External port) const + { + return input_indices_[static_cast(port)]; + } + + PhasorDynamics::Bus bus; + PhasorDynamics::Exciter::Esdc1a esdc1a; + }; + + static constexpr RealT kStateVr = 0.9; + static constexpr RealT kStateVi = 0.4; + + Data makeMinimalData() const + { + Data data; + data.device_class = "Esdc1a"; + data.disambiguation_string = "esdc1a_test"; + data.monitored_variables.insert(Mon::efd); + data.monitored_variables.insert(Mon::vc); + data.monitored_variables.insert(Mon::vr); + data.monitored_variables.insert(Mon::vf); + data.monitored_variables.insert(Mon::se); + data.monitored_variables.insert(Mon::vfe); + return data; + } + + Data makeExplicitDefaultData() const + { + auto data = makeMinimalData(); + + // These are the documented defaults. + data.parameters[Params::Tr] = 0.0; + data.parameters[Params::Ka] = 40.0; + data.parameters[Params::Ta] = 0.1; + data.parameters[Params::Tb] = 0.0; + data.parameters[Params::Tc] = 0.0; + data.parameters[Params::Vrmax] = 1.0; + data.parameters[Params::Vrmin] = -1.0; + data.parameters[Params::Ke] = 0.1; + data.parameters[Params::Te] = 0.5; + data.parameters[Params::Kf] = 0.05; + data.parameters[Params::Tf1] = 0.7; + data.parameters[Params::Spdmlt] = false; + data.parameters[Params::E1] = 2.8; + data.parameters[Params::Se1] = 0.08; + data.parameters[Params::E2] = 3.7; + data.parameters[Params::Se2] = 0.33; + data.parameters[Params::UEL] = static_cast(0); + data.parameters[Params::exclim] = true; + return data; + } + + Data makeData() const + { + auto data = makeMinimalData(); + + // The documented typical values with the floored time constants + // raised above the floor, so routine fixtures log no warnings. + data.parameters[Params::Tr] = 0.02; + data.parameters[Params::Ka] = 40.0; + data.parameters[Params::Ta] = 0.1; + data.parameters[Params::Tb] = 0.5; + data.parameters[Params::Tc] = 0.0; + data.parameters[Params::Vrmax] = 1.0; + data.parameters[Params::Vrmin] = -1.0; + data.parameters[Params::Ke] = 0.1; + data.parameters[Params::Te] = 0.5; + data.parameters[Params::Kf] = 0.05; + data.parameters[Params::Tf1] = 0.7; + data.parameters[Params::Spdmlt] = false; + data.parameters[Params::E1] = 2.8; + data.parameters[Params::Se1] = 0.08; + data.parameters[Params::E2] = 3.7; + data.parameters[Params::Se2] = 0.33; + data.parameters[Params::UEL] = static_cast(0); + data.parameters[Params::exclim] = true; + return data; + } + + Data makeResidualData() const + { + auto data = makeData(); + + // Dynamic-response parameters: every gain, lag, and saturation + // coefficient is nontrivial, the lead-lag is active, and the speed + // multiplier is enabled. + data.parameters[Params::Tr] = 0.2; + data.parameters[Params::Ka] = 25.0; + data.parameters[Params::Ta] = 0.3; + data.parameters[Params::Tb] = 0.8; + data.parameters[Params::Tc] = 0.3; + data.parameters[Params::Vrmax] = 2.0; + data.parameters[Params::Vrmin] = -2.0; + data.parameters[Params::Ke] = 0.2; + data.parameters[Params::Te] = 0.6; + data.parameters[Params::Kf] = 0.08; + data.parameters[Params::Tf1] = 0.9; + data.parameters[Params::Spdmlt] = true; + data.parameters[Params::E1] = 2.4; + data.parameters[Params::Se1] = 0.1; + data.parameters[Params::E2] = 3.2; + data.parameters[Params::Se2] = 0.5; + return data; + } + + /// The external inputs the residual answer key is evaluated against. + template + void setAnswerKeyInputs(Fixture& fixture) const + { + fixture.input(External::OMEGA) = 0.03; + fixture.input(External::VREF) = 1.05; + fixture.input(External::VS) = 0.04; + fixture.input(External::VUEL) = 0.334; + } + + /// The rich state shared by the residual answer key and the Jacobian + /// comparison. Every row is distinct so a swapped index cannot pass, + /// and VLL sits close enough to VUEL that the smooth gate keeps + /// two-sided sensitivity. + template + void setAnswerKeyState(PhasorDynamics::Exciter::Esdc1a& esdc1a) const + { + setState(esdc1a, + {{Internal::EFDP, 2.00}, + {Internal::VC, 0.95}, + {Internal::VR, 0.45}, + {Internal::VF, 0.06}, + {Internal::XLL, 0.30}, + {Internal::EV, 0.36}, + {Internal::VLL, 0.33}, + {Internal::VHV, 0.02}, + {Internal::SE, 0.09}, + {Internal::VFE, 0.42}, + {Internal::EFD, 1.15}}); + setDerivative(esdc1a, + {{Internal::EFDP, 0.01}, + {Internal::VC, -0.02}, + {Internal::VR, 0.03}, + {Internal::VF, -0.04}, + {Internal::XLL, 0.05}}); + } + + /// Omitting every parameter must give exactly the model built from the + /// defaults the README documents, at rest and under load. + bool defaultsMatchDocumentedValues() const + { + Fixture implicit_defaults(makeMinimalData(), kStateVr, kStateVi); + Fixture explicit_defaults(makeExplicitDefaultData(), kStateVr, kStateVi); + implicit_defaults.attachAllInputs(); + explicit_defaults.attachAllInputs(); + + bool success = implicit_defaults.initialize(1.2) + && explicit_defaults.initialize(1.2); + if (!success) + { + std::cout << "ESDC1A documented-default comparison failed to initialize\n"; + return false; + } + + if (implicit_defaults.evaluate() != 0) + { + success = false; + } + if (explicit_defaults.evaluate() != 0) + { + success = false; + } + if (!vectorUnchanged(implicit_defaults.esdc1a.y(), + copyVector(explicit_defaults.esdc1a.y()), + "documented-default state")) + { + success = false; + } + if (!vectorUnchanged(implicit_defaults.esdc1a.yp(), + copyVector(explicit_defaults.esdc1a.yp()), + "documented-default derivative")) + { + success = false; + } + if (!vectorUnchanged(implicit_defaults.esdc1a.getResidual(), + copyVector(explicit_defaults.esdc1a.getResidual()), + "documented-default residual")) + { + success = false; + } + + setAnswerKeyInputs(implicit_defaults); + setAnswerKeyInputs(explicit_defaults); + setAnswerKeyState(implicit_defaults.esdc1a); + setAnswerKeyState(explicit_defaults.esdc1a); + if (implicit_defaults.evaluate() != 0) + { + success = false; + } + if (explicit_defaults.evaluate() != 0) + { + success = false; + } + if (!vectorUnchanged(implicit_defaults.esdc1a.getResidual(), + copyVector(explicit_defaults.esdc1a.getResidual()), + "documented-default dynamic residual")) + { + success = false; + } + return success; + } + + template + bool invalidParameterCase(Params parameter, ValueT value) const + { + auto data = makeData(); + data.parameters[parameter] = value; + Fixture fixture(data); + return fixture.esdc1a.verify() > 0; + } + + template + bool unlinkedSignalRejected() const + { + PhasorDynamics::SignalNode unlinked_node; + Fixture fixture(makeData()); + fixture.esdc1a.getSignals().template attachSignalNode(&unlinked_node); + return fixture.esdc1a.verify() > 0; + } + + template + std::vector copyVector(const VectorT& vector) const + { + const auto* values = vector.getData(); + return std::vector(values, + values + static_cast(vector.getSize())); + } + + /// Every row of a vector still holds its snapshot value. + template + bool vectorUnchanged(const VectorT& vector, + const std::vector& snapshot, + const char* what) const + { + bool success = true; + const auto* values = vector.getData(); + for (size_t i = 0; i < snapshot.size(); ++i) + { + if (!rowMatches(static_cast(values[i]), snapshot[i], what, i, "changed")) + { + success = false; + } + } + return success; + } + + /// An initialization input retains exactly the value supplied by its + /// owner, including signed infinities and NaN. + bool scalarPreserved(RealT actual, + RealT expected, + const char* what, + size_t row) const + { + bool ret = actual == expected; + if (std::isnan(expected)) + { + ret = std::isnan(actual); + } + if (!ret) + { + std::cout << "ESDC1A " << what << " row " << row + << " changed mismatch: " << actual << " != " << expected << "\n"; + } + return ret; + } + + /// Fill the state and derivative with a recognizable ramp, then re-seed + /// the aliased efd entry, so any write by a rejected initialization + /// is visible. + void poisonState(Fixture& fixture, RealT efd) const + { + auto* y = fixture.esdc1a.y().getData(); + auto* yp = fixture.esdc1a.yp().getData(); + for (size_t i = 0; i < static_cast(fixture.esdc1a.y().getSize()); ++i) + { + y[i] = 0.125 + 0.01 * static_cast(i); + yp[i] = -0.25 - 0.01 * static_cast(i); + } + fixture.seedEfd(efd); + fixture.esdc1a.y().setDataUpdated(); + fixture.esdc1a.yp().setDataUpdated(); + } + + bool initializationRejectedAtomically(const Data& data, + RealT efd_seed, + const ExternalRows& inputs, + const char* label, + RealT vr = 0.8, + RealT vi = 0.6) const + { + Fixture fixture(data, vr, vi); + fixture.attachAllInputs(); + for (const auto& [port, value] : inputs) + { + fixture.input(port) = value; + } + if (!fixture.prepare(efd_seed)) + { + return false; + } + + poisonState(fixture, efd_seed); + const auto y_before = copyVector(fixture.esdc1a.y()); + const auto yp_before = copyVector(fixture.esdc1a.yp()); + + bool success = true; + if (fixture.esdc1a.initialize() == 0) + { + std::cout << "Expected initialization rejection: " << label << "\n"; + success = false; + } + + if (!scalarMatches(fixture.efd(), efd_seed, "rejected efd preservation")) + { + success = false; + } + for (const auto& [port, value] : inputs) + { + if (!scalarPreserved(static_cast(fixture.input(port)), + value, + "external input", + static_cast(port))) + { + success = false; + } + } + if (!vectorUnchanged(fixture.esdc1a.y(), y_before, "state")) + { + success = false; + } + if (!vectorUnchanged(fixture.esdc1a.yp(), yp_before, "derivative")) + { + success = false; + } + return success; + } + + /// Write state rows and publish the update, folding in the + /// setDataUpdated() that a hand-written write block has to remember. + template + void setState(PhasorDynamics::Exciter::Esdc1a& esdc1a, + const InternalRows& rows) const + { + auto* y = esdc1a.y().getData(); + for (const auto& [variable, value] : rows) + { + y[static_cast(variable)] = static_cast(value); + } + esdc1a.y().setDataUpdated(); + } + + /// setState() for the derivative vector. + template + void setDerivative(PhasorDynamics::Exciter::Esdc1a& esdc1a, + const InternalRows& rows) const + { + auto* yp = esdc1a.yp().getData(); + for (const auto& [variable, value] : rows) + { + yp[static_cast(variable)] = static_cast(value); + } + esdc1a.yp().setDataUpdated(); + } + + /// Compare one vector row against its expected value. Every row check + /// in this suite reports through here, so failures share one format. + /// Rows are named by their canonical internal-variable enumeration. + static bool rowMatches(RealT actual, + RealT expected, + const char* what, + size_t row, + const char* context) + { + if (isEqual(actual, expected, kTol)) + { + return true; + } + std::cout << "ESDC1A " << what << " row " << row << ' ' << context + << " mismatch: " << std::setprecision(16) << actual + << " != " << expected << '\n'; + return false; + } + + /// Check selected rows of a model vector against expected values. + template + bool rowsMatch(const VectorT& vector, + const RowsT& rows, + const char* what, + const char* context) const + { + bool success = true; + const auto* values = vector.getData(); + for (const auto& [variable, expected] : rows) + { + const auto row = static_cast(variable); + if (!rowMatches(static_cast(values[row]), expected, what, row, context)) + { + success = false; + } + } + return success; + } + + bool residualsMatch(const Esdc1aT& esdc1a, + const InternalRows& rows, + const char* context = "") const + { + return rowsMatch(esdc1a.getResidual(), rows, "residual", context); + } + + template + bool residualsMatch(const Esdc1aT& esdc1a, + const std::array& rows, + const char* context = "") const + { + return rowsMatch(esdc1a.getResidual(), rows, "residual", context); + } + + /// The model sits at a steady state: every residual and every + /// derivative is zero. + bool allResidualsZero(const Esdc1aT& esdc1a) const + { + bool success = true; + const auto* f = esdc1a.getResidual().getData(); + const auto* yp = esdc1a.yp().getData(); + for (size_t row = 0; row < static_cast(esdc1a.getResidual().getSize()); ++row) + { + if (!rowMatches(static_cast(f[row]), 0.0, "residual", row, "at rest")) + { + success = false; + } + if (!rowMatches(static_cast(yp[row]), 0.0, "derivative", row, "at rest")) + { + success = false; + } + } + return success; + } + + bool scalarMatches(ScalarT actual, + ScalarT expected, + const char* label, + ScalarT tolerance = kTol) const + { + if (isEqual(actual, expected, tolerance)) + { + return true; + } + std::cout << label << " mismatch: " << std::setprecision(16) << actual + << " != " << expected << "\n"; + return false; + } + + void noteExpectedLogs(const char* message) const + { + const auto previous_verbosity = Log::verbosity(); + Log::setVerbosity(Log::Verbosity::EVERYTHING); + Log::misc() << message << "\n"; + Log::setVerbosity(previous_verbosity); + } + + void numberVariables(Fixture& fixture) const + { + auto* y = fixture.esdc1a.y().getData(); + auto* yp = fixture.esdc1a.yp().getData(); + auto* bus_y = fixture.bus.y().getData(); + + const auto model_size = static_cast(fixture.esdc1a.size()); + for (size_t i = 0; i < model_size; ++i) + { + y[i].setVariableNumber(i); + yp[i].setVariableNumber(i); + } + for (size_t i = 0; i < static_cast(fixture.bus.size()); ++i) + { + bus_y[i].setVariableNumber(model_size + i); + } + for (External port : {External::OMEGA, External::VREF, External::VS, External::VUEL}) + { + fixture.input(port).setVariableNumber(fixture.inputIndex(port)); + } + + fixture.esdc1a.y().setDataUpdated(); + fixture.esdc1a.yp().setDataUpdated(); + fixture.bus.y().setDataUpdated(); + } + +#ifdef GRIDKIT_ENABLE_ENZYME + std::vector dependencyTrackingJacobian( + const Data& data, + TestStatus& success) const + { + using DepVar = DependencyTracking::Variable; + + Fixture fixture(data, kStateVr, kStateVi); + fixture.attachAllInputs(); + success *= fixture.initialize(1.2); + setAnswerKeyInputs(fixture); + setAnswerKeyState(fixture.esdc1a); + numberVariables(fixture); + success *= (fixture.evaluate() == 0); + + const auto model_size = static_cast(fixture.esdc1a.size()); + std::vector rows(model_size); + const auto* f = fixture.esdc1a.getResidual().getData(); + for (size_t i = 0; i < model_size; ++i) + { + rows[i] = f[i].getDependencies(); + } + return rows; + } + + std::vector enzymeJacobian( + const Data& data, + TestStatus& success) const + { + Fixture fixture(data, kStateVr, kStateVi); + fixture.attachAllInputs(); + success *= fixture.initialize(1.2); + + for (IdxT i = 0; i < fixture.bus.size(); ++i) + { + fixture.bus.setVariableIndex(i, fixture.esdc1a.size() + i); + } + + setAnswerKeyInputs(fixture); + setAnswerKeyState(fixture.esdc1a); + fixture.esdc1a.updateTime(0.0, 1.0); + success *= (fixture.evaluate() == 0); + success *= (fixture.esdc1a.evaluateJacobian() == 0); + success *= (fixture.esdc1a.constructCsr() == 0); + return MapFromCsr(fixture.esdc1a.getCsrJacobian()); + } +#endif + }; + } // namespace Testing +} // namespace GridKit diff --git a/tests/UnitTests/PhasorDynamics/SystemSingleComponentTests.hpp b/tests/UnitTests/PhasorDynamics/SystemSingleComponentTests.hpp index 34e5ffa3e..dae9bbb57 100644 --- a/tests/UnitTests/PhasorDynamics/SystemSingleComponentTests.hpp +++ b/tests/UnitTests/PhasorDynamics/SystemSingleComponentTests.hpp @@ -118,6 +118,67 @@ namespace GridKit return success.report(__func__); } + /// ESDC1A through the production path: model data to system + /// construction to required bus and field-voltage signal wiring. UEL + /// mode 2 makes zero field voltage an admissible standalone state. + TestOutcome esdc1a() + { + using Data = PhasorDynamics::Exciter::Esdc1aData; + using Buses = typename Data::Buses; + using Outputs = typename Data::SignalOutputs; + using Params = typename Data::Parameters; + using Vars = PhasorDynamics::Exciter::Esdc1aInternalVariables; + + constexpr IdxT bus_id = static_cast(1); + constexpr IdxT efd_id = static_cast(1); + + TestStatus success = true; + + PhasorDynamics::SystemModelData data; + data.bus.resize(1); + data.bus[0].bus_id = bus_id; + data.bus[0].bus_type = PhasorDynamics::BusData::BusType::SLACK; + data.bus[0].Vr0 = static_cast(1.0); + data.bus[0].Vi0 = static_cast(0.0); + + data.signal.resize(1); + data.signal[0].signal_id = efd_id; + data.signal[0].name = "Field Voltage"; + + Data esdc1a_data; + esdc1a_data.device_class = "Esdc1a"; + esdc1a_data.disambiguation_string = "esdc1a_system"; + esdc1a_data.buses[Buses::bus] = bus_id; + esdc1a_data.parameters[Params::Tr] = static_cast(0.02); + esdc1a_data.parameters[Params::Tb] = static_cast(0.5); + esdc1a_data.parameters[Params::UEL] = static_cast(2); + esdc1a_data.signal_outputs[Outputs::efd] = efd_id; + data.esdc1a.push_back(esdc1a_data); + + PhasorDynamics::SystemModel system(data); + + success *= system.allocate() == 0; + success *= system.initialize() == 0; + success *= system.tagDifferentiable() == 0; + success *= system.evaluateResidual() == 0; + success *= system.evaluateJacobian() == 0; + success *= system.size() == static_cast(Vars::MAXIMUM); + + auto* efd = system.getSignal(efd_id); + success *= efd->linked(); + success *= efd->getVariableIndex() == static_cast(Vars::EFD); + + auto missing_bus_data = data; + missing_bus_data.bus[0].bus_id = static_cast(0); + missing_bus_data.esdc1a[0].buses.clear(); + + PhasorDynamics::SystemModel missing_bus_system(missing_bus_data); + std::cout << "Testing expected ESDC1A missing-bus configuration error.\n"; + success *= missing_bus_system.verify() > 0; + + return success.report(__func__); + } + TestOutcome load() { TestStatus success = true; diff --git a/tests/UnitTests/PhasorDynamics/runComponentConnectionTests.cpp b/tests/UnitTests/PhasorDynamics/runComponentConnectionTests.cpp new file mode 100644 index 000000000..b9e9253d1 --- /dev/null +++ b/tests/UnitTests/PhasorDynamics/runComponentConnectionTests.cpp @@ -0,0 +1,13 @@ +#include + +#include "ComponentConnectionTests.hpp" + +int main() +{ + GridKit::Testing::TestingResults result; + GridKit::Testing::ComponentConnectionTests test; + + result += test.genrouEsdc1a(); + + return result.summary(); +} diff --git a/tests/UnitTests/PhasorDynamics/runExciterEsdc1aTests.cpp b/tests/UnitTests/PhasorDynamics/runExciterEsdc1aTests.cpp new file mode 100644 index 000000000..e3630861c --- /dev/null +++ b/tests/UnitTests/PhasorDynamics/runExciterEsdc1aTests.cpp @@ -0,0 +1,20 @@ +#include "ExciterEsdc1aTests.hpp" + +int main() +{ + GridKit::Testing::TestingResults result; + + GridKit::Testing::ExciterEsdc1aTests test; + + result += test.validation(); + result += test.initializationAndSignals(); + result += test.initializationDomain(); + result += test.residualEquations(); + result += test.voltageRegulation(); + result += test.excitationLimits(); +#ifdef GRIDKIT_ENABLE_ENZYME + result += test.jacobian(); +#endif + + return result.summary(); +} diff --git a/tests/UnitTests/PhasorDynamics/runSystemSingleComponentTests.cpp b/tests/UnitTests/PhasorDynamics/runSystemSingleComponentTests.cpp index 2eb8cecb8..ebf13c6fb 100644 --- a/tests/UnitTests/PhasorDynamics/runSystemSingleComponentTests.cpp +++ b/tests/UnitTests/PhasorDynamics/runSystemSingleComponentTests.cpp @@ -12,6 +12,7 @@ int main() result += test.bus(); result += test.busFault(); result += test.ieeet1(); + result += test.esdc1a(); result += test.load(); result += test.loadZIP(); result += test.regca(); diff --git a/tests/UnitTests/Utilities/CaseFormatTests.hpp b/tests/UnitTests/Utilities/CaseFormatTests.hpp index aa9e18ff0..c939b49a6 100644 --- a/tests/UnitTests/Utilities/CaseFormatTests.hpp +++ b/tests/UnitTests/Utilities/CaseFormatTests.hpp @@ -8,6 +8,7 @@ #include #include #include +#include #include #include #include @@ -183,8 +184,9 @@ namespace GridKit TestOutcome signalParse() { using namespace GridKit::PhasorDynamics; - using BusData = BusData; - using BusType = typename BusData::BusType; + using BusData = BusData; + using BusType = typename BusData::BusType; + using Esdc1aData = Exciter::Esdc1aData; const char data[] = R"({ @@ -204,12 +206,16 @@ namespace GridKit "signals": [ { "signal_id": 1, "name": "Machine Speed Deviation"}, { "signal_id": 2, "name": "Mechanical Power"}, - { "signal_id": 3, "name": "Excitation Field"} + { "signal_id": 3, "name": "Excitation Field"}, + { "signal_id": 4, "name": "Voltage Reference"}, + { "signal_id": 5, "name": "Stabilizer Signal"}, + { "signal_id": 6, "name": "Under-excitation Limiter"} ], "devices": [ { "class": "Branch", "ports": {"bus1":1, "bus2":2}, "id": "BR1", "params": {"R":0.0, "X":0.1, "G":0.0, "B":0.0, "tap":1.05, "phase":0.1} }, { "class": "Genrou", "ports": {"bus":1, "speed": 1, "pmech":2, "efd":3}, "id": "DV1", "params": {"p0":1.0, "q0":0.05013, "H":3.0, "D":0.0, "Ra":0.0, "Tdop":7.0, "Tdopp":0.04, "Tqopp":0.05, "Tqop":0.75, "Xd":2.1, "Xdp":0.2, "Xdpp":0.18, "Xq":0.5, "Xqp": 0.0, "Xqpp":0.18, "Xl":0.15, "S10":0.0, "S12":0.0}, "mon": ["delta", "omega"] }, { "class": "Tgov1", "ports": {"bus":1, "speed": 1, "pmech":2}, "id": "DV2", "params": {"R":0.05, "T1":0.5,"T2":2.5, "T3":7.5, "Pvmax":0.0, "Pvmin":1.0, "Dt":0.0}}, + { "class": "Esdc1a", "ports": {"bus":1, "speed":1, "vref":4, "vs":5, "vuel":6, "efd":3}, "id": "DV5", "params": {"Tr":0.0, "Ka":40.0, "Ta":0.1, "Tb":0.0, "Tc":0.0, "Vrmax":1.0, "Vrmin":-1.0, "Ke":0.1, "Te":0.5, "Kf":0.05, "Tf1":0.7, "Spdmlt":false, "E1":2.8, "Se1":0.08, "E2":3.7, "Se2":0.33, "UEL":0, "exclim":true}, "mon": ["efd", "vc", "vr", "vf", "se", "vfe"] }, { "class": "Ieeet1", "ports": {"bus":1, "speed": 1, "efd":3}, "id": "DV3", "params": {"Tr":0.0, "Ka":50.0, "Ta":0.04, "Ke":-0.06, "Te":0.6, "Kf":0.09, "Tf":1.46, "Vrmin":-1.0, "Vrmax":1.0, "E1":2.8, "E2":3.373, "Se1":0.04, "Se2":0.33, "Ispdlim":0.0}}, { "class": "SexsPti", "ports": {"bus":1, "efd":3}, "id": "DV4", "params": {"Ta":0.1, "Tb":0.5, "Te":0.8, "K":10.0, "Efdmax":5.0, "Efdmin":-5.0}}, { "class": "BusFault", "ports": {"bus":1}, "id": "1", "params": {"state0": false, "R":0.0, "X":1e-3} } @@ -233,6 +239,7 @@ namespace GridKit success *= result.bus_fault.size() == 1; success *= result.genrou.size() == 1; success *= result.gov.size() == 1; + success *= result.esdc1a.size() == 1; success *= result.loadz.size() == 0; success *= result.exciter.size() == 1; success *= result.sexspti.size() == 1; @@ -259,6 +266,12 @@ namespace GridKit success *= result.signal[1].name == "Mechanical Power"; success *= result.signal[2].signal_id == 3; success *= result.signal[2].name == "Excitation Field"; + success *= result.signal[3].signal_id == 4; + success *= result.signal[3].name == "Voltage Reference"; + success *= result.signal[4].signal_id == 5; + success *= result.signal[4].name == "Stabilizer Signal"; + success *= result.signal[5].signal_id == 6; + success *= result.signal[5].name == "Under-excitation Limiter"; success *= std::get(result.branch[0].parameters[BranchParameters::R]) == 0.0; success *= std::get(result.branch[0].parameters[BranchParameters::X]) == 0.1; @@ -309,6 +322,38 @@ namespace GridKit success *= result.gov[0].signal_outputs[Governor::Tgov1SignalOutputs::pmech] == 2; success *= result.gov[0].disambiguation_string == "DV2"; + success *= std::get(result.esdc1a[0].parameters[Esdc1aData::Parameters::Tr]) == 0.0; + success *= std::get(result.esdc1a[0].parameters[Esdc1aData::Parameters::Ka]) == 40.0; + success *= std::get(result.esdc1a[0].parameters[Esdc1aData::Parameters::Ta]) == 0.1; + success *= std::get(result.esdc1a[0].parameters[Esdc1aData::Parameters::Tb]) == 0.0; + success *= std::get(result.esdc1a[0].parameters[Esdc1aData::Parameters::Tc]) == 0.0; + success *= std::get(result.esdc1a[0].parameters[Esdc1aData::Parameters::Vrmax]) == 1.0; + success *= std::get(result.esdc1a[0].parameters[Esdc1aData::Parameters::Vrmin]) == -1.0; + success *= std::get(result.esdc1a[0].parameters[Esdc1aData::Parameters::Ke]) == 0.1; + success *= std::get(result.esdc1a[0].parameters[Esdc1aData::Parameters::Te]) == 0.5; + success *= std::get(result.esdc1a[0].parameters[Esdc1aData::Parameters::Kf]) == 0.05; + success *= std::get(result.esdc1a[0].parameters[Esdc1aData::Parameters::Tf1]) == 0.7; + success *= !std::get(result.esdc1a[0].parameters[Esdc1aData::Parameters::Spdmlt]); + success *= std::get(result.esdc1a[0].parameters[Esdc1aData::Parameters::E1]) == 2.8; + success *= std::get(result.esdc1a[0].parameters[Esdc1aData::Parameters::Se1]) == 0.08; + success *= std::get(result.esdc1a[0].parameters[Esdc1aData::Parameters::E2]) == 3.7; + success *= std::get(result.esdc1a[0].parameters[Esdc1aData::Parameters::Se2]) == 0.33; + success *= std::get(result.esdc1a[0].parameters[Esdc1aData::Parameters::UEL]) == 0; + success *= std::get(result.esdc1a[0].parameters[Esdc1aData::Parameters::exclim]); + success *= result.esdc1a[0].buses[Esdc1aData::Buses::bus] == 1; + success *= result.esdc1a[0].signal_inputs[Esdc1aData::SignalInputs::speed] == 1; + success *= result.esdc1a[0].signal_inputs[Esdc1aData::SignalInputs::vref] == 4; + success *= result.esdc1a[0].signal_inputs[Esdc1aData::SignalInputs::vs] == 5; + success *= result.esdc1a[0].signal_inputs[Esdc1aData::SignalInputs::vuel] == 6; + success *= result.esdc1a[0].signal_outputs[Esdc1aData::SignalOutputs::efd] == 3; + success *= result.esdc1a[0].disambiguation_string == "DV5"; + success *= result.esdc1a[0].monitored_variables.contains(Esdc1aData::MonitorableVariables::efd); + success *= result.esdc1a[0].monitored_variables.contains(Esdc1aData::MonitorableVariables::vc); + success *= result.esdc1a[0].monitored_variables.contains(Esdc1aData::MonitorableVariables::vr); + success *= result.esdc1a[0].monitored_variables.contains(Esdc1aData::MonitorableVariables::vf); + success *= result.esdc1a[0].monitored_variables.contains(Esdc1aData::MonitorableVariables::se); + success *= result.esdc1a[0].monitored_variables.contains(Esdc1aData::MonitorableVariables::vfe); + success *= std::get(result.exciter[0].parameters[Exciter::Ieeet1Parameters::Tr]) == 0.0; success *= std::get(result.exciter[0].parameters[Exciter::Ieeet1Parameters::Ka]) == 50.0; success *= std::get(result.exciter[0].parameters[Exciter::Ieeet1Parameters::Ta]) == 0.04;