diff --git a/CHANGELOG.md b/CHANGELOG.md index ded783307..ec2882196 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -76,6 +76,7 @@ - Added EMT model and operator documentation. - Added `REGCA` converter model implementation for PhasorDynamics. - Remove unnecessary data copying while evaluating `PowerElectronics` models, speeding up large simulations by up to 3x +- Added `HYGOV` governor model implementation for PhasorDynamics. ## v0.1 diff --git a/GridKit/Model/PhasorDynamics/ComponentLibrary.hpp b/GridKit/Model/PhasorDynamics/ComponentLibrary.hpp index 0e110fd67..51c2d78e5 100644 --- a/GridKit/Model/PhasorDynamics/ComponentLibrary.hpp +++ b/GridKit/Model/PhasorDynamics/ComponentLibrary.hpp @@ -9,6 +9,7 @@ #include #include #include +#include #include #include #include diff --git a/GridKit/Model/PhasorDynamics/Governor/CMakeLists.txt b/GridKit/Model/PhasorDynamics/Governor/CMakeLists.txt index 7c7269784..fbe71c740 100644 --- a/GridKit/Model/PhasorDynamics/Governor/CMakeLists.txt +++ b/GridKit/Model/PhasorDynamics/Governor/CMakeLists.txt @@ -4,3 +4,4 @@ # ]] add_subdirectory(Tgov1) +add_subdirectory(HYGOV) diff --git a/GridKit/Model/PhasorDynamics/Governor/HYGOV/CMakeLists.txt b/GridKit/Model/PhasorDynamics/Governor/HYGOV/CMakeLists.txt new file mode 100644 index 000000000..1719101b0 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Governor/HYGOV/CMakeLists.txt @@ -0,0 +1,54 @@ +# [[ +# Author(s): +# - Luke Lowery +# ]] + +set(_install_headers Hygov.hpp HygovData.hpp) + +if(GRIDKIT_ENABLE_ENZYME) + gridkit_add_library( + phasor_dynamics_governor_hygov + SOURCES HygovEnzyme.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_governor_hygov + SOURCES Hygov.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_governor_hygov_dependency_tracking + SOURCES HygovDependencyTracking.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_governor_hygov) +target_link_libraries( + phasor_dynamics_components_dependency_tracking + INTERFACE GridKit::phasor_dynamics_governor_hygov_dependency_tracking) diff --git a/GridKit/Model/PhasorDynamics/Governor/HYGOV/Hygov.cpp b/GridKit/Model/PhasorDynamics/Governor/HYGOV/Hygov.cpp new file mode 100644 index 000000000..bdb7e15c0 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Governor/HYGOV/Hygov.cpp @@ -0,0 +1,27 @@ +/** + * @file Hygov.cpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Non-Enzyme instantiation for the HYGOV governor model. + */ + +#include "HygovImpl.hpp" + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Governor + { + template + int Hygov::evaluateJacobian() + { + Log::misc() << "Evaluate Jacobian for Hygov..." << std::endl; + Log::misc() << "Jacobian evaluation is not implemented!" << std::endl; + return 0; + } + + template class Hygov; + template class Hygov; + } // namespace Governor + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Governor/HYGOV/Hygov.hpp b/GridKit/Model/PhasorDynamics/Governor/HYGOV/Hygov.hpp new file mode 100644 index 000000000..9844dffd3 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Governor/HYGOV/Hygov.hpp @@ -0,0 +1,199 @@ +/** + * @file Hygov.hpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Declaration of the HYGOV governor model. + */ + +#pragma once + +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +namespace GridKit +{ + namespace PhasorDynamics + { + template + class SignalNode; + + namespace Governor + { + /// Internal variables of a `Hygov`. + enum class HygovInternalVariables : size_t + { + XN, ///< \f$x_n\f$ Speed lead-lag denominator state + XF, ///< \f$x_f\f$ Governor error filter output on component base + C, ///< \f$c\f$ Desired-gate position on component base + G, ///< \f$g\f$ Gate position on component base + Q, ///< \f$q\f$ Turbine flow on component base + OMEGADB, ///< \f$\omega_{\mathrm{db}}\f$ Deadbanded speed deviation + EF, ///< \f$e_f\f$ Governor error on component base + FC, ///< \f$f_c\f$ Desired-gate derivative target + RC, ///< \f$r_c\f$ Rate-limited desired-gate derivative target + PGV, ///< \f$P_{\mathrm{GV}}\f$ Gate-to-power curve output on component base + H, ///< \f$H\f$ Turbine head on component base + PMECH, ///< \f$P_{\mathrm{m}}\f$ Mechanical-power output on system base + MAXIMUM, + }; + + /// External variables of a `Hygov`. + enum class HygovExternalVariables : size_t + { + OMEGA, ///< \f$\omega\f$ Machine speed deviation + PREF, ///< \f$P^{\mathrm{ref}}\f$ Active-power/load reference on system base + PAUX, ///< \f$P^{\mathrm{aux}}\f$ Auxiliary power input on system base + MAXIMUM, + }; + + /** + * @brief Hydro turbine-governor model with temporary droop, gate servo, + * and a nonlinear single-penstock turbine. + * + * @tparam scalar_type Plain real or differentiable scalar type. + * @tparam index_type Integer index type. + */ + template + class Hygov : 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::va_system_base_; + 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 SignalT = SignalNode; + using ModelDataT = HygovData; + using MonitorT = Model::VariableMonitor; + using InternalVariablesT = HygovInternalVariables; + using ExternalVariablesT = HygovExternalVariables; + + Hygov(); + explicit Hygov(const ModelDataT& data); + ~Hygov(); + + 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(); + + /// Evaluate the nonlinear gate-to-power curve as a fixed sum of + /// smooth linear segments. + __attribute__((always_inline)) inline ScalarT gatePower(ScalarT gate) const; + + /// Steady component-base mechanical power at a gate and dam head. + RealT initialMechanicalPower(RealT gate, RealT Hdam) const; + + /// Bisect a bracketed initialization residual to machine rounding. + template + static RealT bisectInitialRoot(RealT a, + RealT b, + RealT fa, + RealT fb, + FuncT residual); + + /// Solve the gate at the configured dam head. + RealT solveInitialGate(RealT pmech) const; + + /// Solve the dam head that reproduces mechanical power at Gv5. + RealT solveInitialDamHead(RealT pmech) const; + + ScalarT toComponentBase(ScalarT value) const; + ScalarT toSystemBase(ScalarT value) const; + + static constexpr RealT TIME_CONSTANT_MINIMUM = static_cast(1.0e-3); + + /// Accepted seed distance beyond the achievable-power range edge. + static constexpr RealT INITIALIZATION_TOLERANCE = + static_cast(100.0) * std::numeric_limits::epsilon(); + + RealT Rperm_{static_cast(0.04)}; + RealT Rtemp_{static_cast(0.3)}; + RealT Tr_{static_cast(5.0)}; + RealT Tf_{static_cast(0.05)}; + RealT Tg_{static_cast(0.5)}; + RealT Velm_{static_cast(0.2)}; + RealT Gmax_{ONE}; + RealT Gmin_{ZERO}; + RealT Tw_{ONE}; + RealT At_{static_cast(1.2)}; + RealT Dturb_{static_cast(0.5)}; + RealT Qnl_{static_cast(0.05)}; + RealT Tn_{ZERO}; + RealT Tnp_{ZERO}; + RealT db1_{ZERO}; + RealT db2_{ZERO}; + RealT Hdam_{ONE}; + std::array Gv_{}; + std::array Pgv_{}; + + RealT va_component_base_{ZERO}; + RealT leadlag_gain_{ZERO}; + + IdxT parameter_error_count_{0}; + + RealT Gmin_response_{Gmin_}; + RealT Gmax_response_{Gmax_}; + RealT Hdam_eff_{Hdam_}; + ScalarT pref_set_{0}; + ScalarT paux_set_{0}; + + ComponentSignals signals_; + std::unique_ptr monitor_; + + std::vector ws_; + std::vector ws_indices_; + }; + } // namespace Governor + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovData.hpp b/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovData.hpp new file mode 100644 index 000000000..733d39135 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovData.hpp @@ -0,0 +1,113 @@ +/** + * @file HygovData.hpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Modeling data for the HYGOV governor model. + */ + +#pragma once + +#include + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Governor + { + /// Parameter keys for the HYGOV governor model. + enum class HygovParameters + { + Trate, ///< \f$T^\mathrm{rate}\f$ Turbine-rating power base + Rperm, ///< \f$R_{\mathrm{perm}}\f$ Permanent droop + Rtemp, ///< \f$R_{\mathrm{temp}}\f$ Temporary droop + Tr, ///< \f$T_r\f$ Temporary-droop reset time constant + Tf, ///< \f$T_f\f$ Governor error filter time constant + Tg, ///< \f$T_g\f$ Gate servo time constant + Velm, ///< \f$V_{\mathrm{elm}}\f$ Maximum desired-gate velocity magnitude + Gmax, ///< \f$G^{\max}\f$ Configured upper gate response limit + Gmin, ///< \f$G^{\min}\f$ Configured lower gate response limit + Tw, ///< \f$T_w\f$ Water inertia time constant + At, ///< \f$A_t\f$ Turbine gain + Dturb, ///< \f$D_{\mathrm{turb}}\f$ Turbine damping coefficient + Qnl, ///< \f$q_{\mathrm{NL}}\f$ No-load flow at nominal head + Tn, ///< \f$T_n\f$ Speed lead-lag numerator time constant + Tnp, ///< \f$T_{\mathrm{np}}\f$ Speed lead-lag denominator time constant + db1, ///< \f$D_{\omega}\f$ Type 1 speed deadband threshold + db2, ///< \f$D_2\f$ Unsupported mechanical backlash. Nonzero values warn and are ignored + Hdam, ///< \f$H_{\mathrm{dam}}\f$ Head available at dam + Gv0, ///< \f$G_V^{(0)}\f$ Gate point 0 + Gv1, ///< \f$G_V^{(1)}\f$ Gate point 1 + Gv2, ///< \f$G_V^{(2)}\f$ Gate point 2 + Gv3, ///< \f$G_V^{(3)}\f$ Gate point 3 + Gv4, ///< \f$G_V^{(4)}\f$ Gate point 4 + Gv5, ///< \f$G_V^{(5)}\f$ Gate point 5 + Pgv0, ///< \f$P_{\mathrm{GV}}^{(0)}\f$ Power point 0 + Pgv1, ///< \f$P_{\mathrm{GV}}^{(1)}\f$ Power point 1 + Pgv2, ///< \f$P_{\mathrm{GV}}^{(2)}\f$ Power point 2 + Pgv3, ///< \f$P_{\mathrm{GV}}^{(3)}\f$ Power point 3 + Pgv4, ///< \f$P_{\mathrm{GV}}^{(4)}\f$ Power point 4 + Pgv5 ///< \f$P_{\mathrm{GV}}^{(5)}\f$ Power point 5 + }; + + /// Buses for the HYGOV governor model. + enum class HygovBuses : size_t + { + SIZE + }; + + /// Signal inputs for the HYGOV governor model. + enum class HygovSignalInputs : size_t + { + speed, ///< Optional machine speed-deviation signal ID + pref, ///< Optional active-power/load reference signal ID + paux, ///< Optional auxiliary power input signal ID + SIZE + }; + + /// Signal outputs for the HYGOV governor model. + enum class HygovSignalOutputs : size_t + { + pmech, ///< Required mechanical-power output signal ID + SIZE + }; + + /// Variables available through the monitor interface. + enum class HygovMonitorableVariables + { + pmech, ///< Mechanical power output on system base + filter, ///< Governor error filter output on component base + desiredgate, ///< Desired-gate position on component base + gate, ///< Gate position on component base + flow, ///< Turbine flow on component base + head ///< Turbine head on component base + }; + + /** + * @brief Model data for HYGOV: parameters, optional input signals, the + * required mechanical-power output, and monitored variables. + * + * @tparam real_type Real parameter value type. + * @tparam index_type Integer index type. + * + * @see Hygov + */ + template + struct HygovData : public ComponentData + { + HygovData() = default; + + using Parameters = HygovParameters; + using Buses = HygovBuses; + using SignalInputs = HygovSignalInputs; + using SignalOutputs = HygovSignalOutputs; + using MonitorableVariables = HygovMonitorableVariables; + }; + } // namespace Governor + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovDependencyTracking.cpp b/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovDependencyTracking.cpp new file mode 100644 index 000000000..760bac957 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovDependencyTracking.cpp @@ -0,0 +1,27 @@ +/** + * @file HygovDependencyTracking.cpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Dependency-tracking instantiations for the HYGOV governor model. + */ + +#include "HygovImpl.hpp" + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Governor + { + template + int Hygov::evaluateJacobian() + { + Log::misc() << "Evaluate Jacobian for Hygov..." << std::endl; + Log::misc() << "Jacobian evaluation is not implemented!" << std::endl; + return 0; + } + + template class Hygov; + template class Hygov; + } // namespace Governor + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovEnzyme.cpp b/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovEnzyme.cpp new file mode 100644 index 000000000..174f23959 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovEnzyme.cpp @@ -0,0 +1,90 @@ +/** + * @file HygovEnzyme.cpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Enzyme sparse Jacobian for the HYGOV governor model. + */ + +#include + +#include "HygovImpl.hpp" + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Governor + { + template + int Hygov::evaluateJacobian() + { + Log::misc() << "Evaluate Jacobian for Hygov..." << std::endl; + Log::misc() << "Jacobian evaluation is experimental!" << std::endl; + + if (J_rows_buffer_ == nullptr) + { + auto size = static_cast(size_); + auto signal_size = static_cast(ws_.size()); + auto buffer_size = 2 * size * 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::Governor::Hygov; + 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::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 Hygov; + template class Hygov; + } // namespace Governor + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovImpl.hpp b/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovImpl.hpp new file mode 100644 index 000000000..9761ca9c6 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovImpl.hpp @@ -0,0 +1,1089 @@ +/** + * @file HygovImpl.hpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Definition of the HYGOV governor model. + */ + +#pragma once + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Governor + { + using Log = ::GridKit::Utilities::Logger; + + /** + * @brief Construct a HYGOV governor without parameters + * + * The model is sized with the documented parameter defaults but without + * a monitor or assigned mechanical-power output. + */ + template + Hygov::Hygov() + { + size_ = static_cast(HygovInternalVariables::MAXIMUM); + } + + /** + * @brief Construct a HYGOV governor from model data + * + * @param[in] data Parameters and monitored-variable selections. + */ + template + Hygov::Hygov(const ModelDataT& data) + : monitor_(std::make_unique(data)) + { + initializeParameters(data); + initializeMonitor(); + size_ = static_cast(HygovInternalVariables::MAXIMUM); + } + + template + Hygov::~Hygov() + { + } + + /** + * @brief Set the component ID + * + * @param[in] component_id Identifier assigned by the system model. + * @return int 0 on success. + */ + template + int Hygov::setGridKitComponentID(IdxT component_id) + { + gridkit_component_id_ = component_id; + return 0; + } + + /** + * @brief Allocate the model vectors and wire the mechanical-power output + * + * Sizes the state, residual, and signal-interface buffers, initializes + * the identity index maps, and points the assigned `pmech` node at the + * internal state it publishes. That node aliases HYGOV storage from + * here on, which is how initialize() reads the machine value. + * HYGOV attaches to no bus, so the bus-interface buffer stays empty. + * Repeated calls reuse the allocated vectors. + * + * @return int 0 on success. + */ + template + int Hygov::allocate() + { + const auto PMECH = static_cast(HygovInternalVariables::PMECH); + + if (!allocated_) + { + this->allocateVectors(size_); + } + auto size = static_cast(size_); + + tag_.assign(size, false); + variable_indices_.resize(size); + residual_indices_.resize(size); + + wb_.clear(); + + const auto signal_size = static_cast(HygovExternalVariables::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); + } + + if (signals_.template isAssigned()) + { + auto* y = y_.getData(); + signals_.template getSignalNode()->set( + &y[PMECH], + &(this->getVariableIndex(static_cast(PMECH)))); + } + + allocated_ = true; + return 0; + } + + /** + * @brief Validate the HYGOV configuration + * + * Checks parameter-loading errors, static parameter relationships, the + * power-base domain, gate-curve shape, gate-limit domain, the assigned + * mechanical-power output, and attached external signals. Mechanical- + * power feasibility is operating-point dependent and is checked by + * initialize(). + * + * @return int Number of configuration errors; zero when valid. + */ + template + int Hygov::verify() const + { + int ret = static_cast(parameter_error_count_); + + auto check = [&](bool condition, const char* message) + { + if (!condition) + { + Log::error() << "Hygov: " << message << '\n'; + ret += 1; + } + }; + + RealT component_power_base = va_component_base_; + const bool component_base_is_omitted = + !(component_power_base > ZERO); + if (component_base_is_omitted) + { + component_power_base = va_system_base_; + } + + const bool valid_component_base = std::isfinite(component_power_base) + && component_power_base > ZERO; + const bool valid_system_base = std::isfinite(va_system_base_) + && va_system_base_ > ZERO; + check(valid_component_base, + "component power base must be finite and positive"); + check(valid_system_base, + "system power base must be finite and positive"); + if (valid_component_base && valid_system_base) + { + const RealT system_to_component = va_system_base_ / component_power_base; + const RealT component_to_system = component_power_base / va_system_base_; + const bool valid_base_ratios = std::isfinite(system_to_component) + && system_to_component > ZERO + && std::isfinite(component_to_system) + && component_to_system > ZERO; + check(valid_base_ratios, + "system/component power-base conversion ratios must be finite and positive"); + } + + check(Rtemp_ > ZERO, "Rtemp must be nonzero"); + check(Tn_ >= ZERO, "Tn must be non-negative"); + check(Velm_ >= ZERO, "Velm must be non-negative"); + check(Gmin_ < Gmax_, "Gmin must be less than Gmax"); + check(At_ > ZERO, "At must be positive"); + check(Dturb_ >= ZERO, "Dturb must be non-negative"); + check(db1_ >= ZERO, "db1 must be non-negative"); + check(Hdam_ > ZERO, "Hdam must be positive"); + + bool curve_shape_is_valid = true; + for (size_t i = 1; i < Gv_.size(); ++i) + { + const bool gate_points_increase = Gv_[i - 1] < Gv_[i]; + const bool power_points_increase = Pgv_[i - 1] <= Pgv_[i]; + + check(gate_points_increase, "Gv points must be strictly increasing"); + check(power_points_increase, "Pgv points must be non-decreasing"); + + if (!gate_points_increase || !power_points_increase) + { + curve_shape_is_valid = false; + } + } + const bool minimum_gate_is_valid = Gv_[0] <= Gmin_; + const bool maximum_gate_is_valid = Gmax_ <= Gv_[5]; + check(minimum_gate_is_valid, "Gmin must be at or above the first Gv point"); + check(maximum_gate_is_valid, "Gmax must be at or below the last Gv point"); + + const bool can_check_power_range = curve_shape_is_valid + && Gmin_ < Gmax_ + && minimum_gate_is_valid + && maximum_gate_is_valid + && At_ > ZERO + && Hdam_ > ZERO; + if (can_check_power_range) + { + // A rise no wider than the tolerance that pins a seed to a range + // edge leaves the gate undetermined by the mechanical power. + const RealT minimum_power = initialMechanicalPower(Gv_[0], Hdam_); + const RealT maximum_power = initialMechanicalPower(Gv_[5], Hdam_); + const RealT power_range = maximum_power - minimum_power; + const bool finite_power_range = std::isfinite(minimum_power) + && std::isfinite(maximum_power) + && std::isfinite(power_range); + check(finite_power_range, + "mechanical-power range must be finite"); + if (finite_power_range) + { + check(power_range > INITIALIZATION_TOLERANCE, + "mechanical power must rise across [Gv0, Gv5]"); + } + } + + check(signals_.template isAssigned(), + "pmech output signal must be assigned"); + + // An attached port must resolve to readable 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() << "Hygov: " << name << " signal attached with no linked source\n"; + ret += 1; + } + }; + + check_attached_signal.template operator()("speed"); + check_attached_signal.template operator()("pref"); + check_attached_signal.template operator()("paux"); + + return ret; + } + + /** + * @brief Initialize HYGOV from the mechanical-power port + * + * Reads the assigned system-base `pmech` node and the attached speed + * and auxiliary-power inputs, solves the component-base steady state + * that preserves the given value, and publishes the resolved load reference + * to an attached `pref` signal. + * + * @pre allocate() has completed. + * @pre The machine model has initialized the assigned `pmech` node. + * + * @post On success the state zeros every residual row at machine + * rounding; a value clipped to the achievable-power range edge + * leaves a mechanical-power residual up to the initialization + * tolerance. + * @post On failure state, effective response limits, effective Hdam, + * and signal storage are unchanged. + * + * @return int 0 on success; nonzero when the configuration or initial + * values are invalid, the initial speed deviation is + * nonzero, or the initial mechanical power cannot be + * reproduced. + */ + template + int Hygov::initialize() + { + const auto XN = static_cast(HygovInternalVariables::XN); + const auto XF = static_cast(HygovInternalVariables::XF); + const auto C = static_cast(HygovInternalVariables::C); + const auto G = static_cast(HygovInternalVariables::G); + const auto Q = static_cast(HygovInternalVariables::Q); + const auto OMEGADB = static_cast(HygovInternalVariables::OMEGADB); + const auto EF = static_cast(HygovInternalVariables::EF); + const auto FC = static_cast(HygovInternalVariables::FC); + const auto RC = static_cast(HygovInternalVariables::RC); + const auto PGV = static_cast(HygovInternalVariables::PGV); + const auto H = static_cast(HygovInternalVariables::H); + const auto PMECH = static_cast(HygovInternalVariables::PMECH); + + bool ret = verify() == 0; + if (!ret) + { + Log::error() << "Hygov: cannot initialize with invalid configuration\n"; + return 1; + } + + ret = va_component_base_ > ZERO; + if (!ret) + { + va_component_base_ = va_system_base_; + } + + auto* y = y_.getData(); + + // The assigned pmech node aliases this entry after allocate(). Its + // system-base value remains untouched throughout initialization. + const ScalarT pmech0_system = y[PMECH]; + + ScalarT omega0{ZERO}; + if (signals_.template isAttached()) + { + omega0 = signals_.template readExternalVariable(); + } + + ScalarT paux0_system{ZERO}; + if (signals_.template isAttached()) + { + paux0_system = signals_.template readExternalVariable(); + } + + auto is_finite = [](ScalarT value) + { + return std::isfinite(static_cast(value)); + }; + ret = is_finite(pmech0_system) + && is_finite(omega0) + && is_finite(paux0_system); + if (!ret) + { + Log::error() << "Hygov: initial pmech, speed, and paux values must be finite\n"; + return 1; + } + + const ScalarT pmech0 = toComponentBase(pmech0_system); + const ScalarT paux0 = toComponentBase(paux0_system); + ret = is_finite(pmech0) + && is_finite(paux0); + if (!ret) + { + Log::error() << "Hygov: initial power-base conversions must be finite\n"; + return 1; + } + + // Synchronous machines provide an exactly zero speed deviation. A + // moving machine would need a multi-root gate search, which this + // model does not support. The speed was verified finite above, so + // this is an exact comparison by intent rather than a tolerance test. + const RealT speed0 = static_cast(omega0); + + ret = speed0 == ZERO; + if (!ret) + { + Log::error() << "Hygov: initialization requires zero speed deviation\n"; + return 1; + } + + const RealT pmech0_value = static_cast(pmech0); + const RealT maximum_power = initialMechanicalPower(Gv_[5], Hdam_); + RealT Hdam0 = Hdam_; + RealT gate0 = Gv_[5]; + + if (pmech0_value > maximum_power) + { + Hdam0 = solveInitialDamHead(pmech0_value); + if (!std::isfinite(Hdam0)) + { + Log::error() << "Hygov: no finite Hdam reproduces the initial mechanical power\n"; + return 1; + } + } + else + { + gate0 = solveInitialGate(pmech0_value); + if (!std::isfinite(gate0)) + { + Log::error() << "Hygov: initial mechanical power is below the first Gv endpoint\n"; + return 1; + } + } + + const RealT Gmin_response = std::min(Gmin_, gate0); + const RealT Gmax_response = std::max(Gmax_, gate0); + + const ScalarT h0 = static_cast(Hdam0); + const ScalarT pgv0 = gatePower(static_cast(gate0)); + const ScalarT q0 = std::sqrt(Hdam0) * pgv0; + const ScalarT omegadb0 = Math::deadband1(omega0, -db1_, db1_); + const ScalarT xn0 = omegadb0; + const ScalarT yomega0 = xn0 + leadlag_gain_ * (omegadb0 - xn0); + const ScalarT pref0 = toSystemBase(yomega0 + Rperm_ * gate0 - paux0); + + ret = is_finite(h0) + && is_finite(pgv0) + && is_finite(q0) + && is_finite(omegadb0) + && is_finite(xn0) + && is_finite(yomega0) + && is_finite(pref0); + if (!ret) + { + Log::error() << "Hygov: initialization produced a nonfinite value\n"; + return 1; + } + + y[XN] = xn0; + y[XF] = ZERO; + y[C] = gate0; + y[G] = gate0; + y[Q] = q0; + y[OMEGADB] = omegadb0; + y[EF] = ZERO; + y[FC] = ZERO; + y[RC] = ZERO; + y[PGV] = pgv0; + y[H] = h0; + + Gmin_response_ = Gmin_response; + Gmax_response_ = Gmax_response; + Hdam_eff_ = Hdam0; + pref_set_ = pref0; + paux_set_ = paux0_system; + + if (signals_.template isAttached()) + { + signals_.template writeExternalVariable(pref_set_); + } + + if (Hdam_eff_ > Hdam_) + { + Log::warning() << "Hygov: effective Hdam raised to match initial mechanical power\n"; + } + if (gate0 < Gmin_ || gate0 > Gmax_) + { + Log::warning() << "Hygov: initial gate is outside [Gmin, Gmax]; " + "response limits are adjusted to include the initialized value\n"; + } + + y_.setDataUpdated(); + yp_.setToConst(static_cast(ZERO)); + return 0; + } + + /** + * @brief Identify the differential variables + * + * The speed lead-lag state, the governor error filter, the desired + * gate, the gate servo, and the turbine flow carry derivatives; every + * other internal variable is algebraic. + * + * @return int 0 on success. + */ + template + int Hygov::tagDifferentiable() + { + const auto XN = static_cast(HygovInternalVariables::XN); + const auto XF = static_cast(HygovInternalVariables::XF); + const auto C = static_cast(HygovInternalVariables::C); + const auto G = static_cast(HygovInternalVariables::G); + const auto Q = static_cast(HygovInternalVariables::Q); + + std::fill(tag_.begin(), tag_.end(), false); + tag_[XN] = true; + tag_[XF] = true; + tag_[C] = true; + tag_[G] = true; + tag_[Q] = true; + return 0; + } + + /** + * @brief Compute the absolute tolerance for each variable in the model + * + * All HYGOV variables are per-unit speeds, gates, flows, heads, and + * powers of the same order, so their absolute and relative tolerance + * have the same value. + * + * @param[in] rel_tol Solver relative tolerance. + * @return int 0 on success. + */ + template + int Hygov::setAbsoluteTolerance(RealT rel_tol) + { + abs_tol_.setToConst(static_cast(rel_tol)); + return 0; + } + + /** + * @brief Residuals of system equations + * + * Refreshes the signal interface buffers and evaluates the internal + * residual. HYGOV attaches to no bus, so there is no bus interface to + * refresh. An unattached reference or auxiliary port falls back to the + * value latched by initialize(); an unattached speed port reads zero + * deviation. + * + * @return int 0 on success. + */ + template + int Hygov::evaluateResidual() + { + const auto OMEGA = static_cast(HygovExternalVariables::OMEGA); + const auto PREF = static_cast(HygovExternalVariables::PREF); + const auto PAUX = static_cast(HygovExternalVariables::PAUX); + + ws_[OMEGA] = ZERO; + ws_[PREF] = pref_set_; + ws_[PAUX] = paux_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_[PREF] = signals_.template readExternalVariable(); + ws_indices_[PREF] = + signals_.template readExternalVariableIndex(); + } + if (signals_.template isAttached()) + { + ws_[PAUX] = signals_.template readExternalVariable(); + ws_indices_[PAUX] = + signals_.template readExternalVariableIndex(); + } + + 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* Hygov::getMonitor() const + { + return monitor_.get(); + } + + /** + * @brief Internal residual + * + * Evaluates the five governor states and the seven algebraic rows + * documented in the model README. The body is kept free of branches + * and loops so that sparse automatic differentiation resolves a fixed + * structure; the gate curve enters as a fixed sum of smooth linear + * segments. + * + * @param[in] y Internal variables. + * @param[in] yp Internal variable derivatives. + * @param[in] wb Bus voltage components; unused, HYGOV attaches to no bus. + * @param[in] ws External signal values on system base. + * @param[out] f Internal residuals. + * @return int 0 on success. + */ + template + __attribute__((always_inline)) inline int + Hygov::evaluateInternalResidual( + const ScalarT* y, + const ScalarT* yp, + [[maybe_unused]] const ScalarT* wb, + const ScalarT* ws, + ScalarT* f) + { + const auto XN = static_cast(HygovInternalVariables::XN); + const auto XF = static_cast(HygovInternalVariables::XF); + const auto C = static_cast(HygovInternalVariables::C); + const auto G = static_cast(HygovInternalVariables::G); + const auto Q = static_cast(HygovInternalVariables::Q); + const auto OMEGADB = static_cast(HygovInternalVariables::OMEGADB); + const auto EF = static_cast(HygovInternalVariables::EF); + const auto FC = static_cast(HygovInternalVariables::FC); + const auto RC = static_cast(HygovInternalVariables::RC); + const auto PGV = static_cast(HygovInternalVariables::PGV); + const auto H = static_cast(HygovInternalVariables::H); + const auto PMECH = static_cast(HygovInternalVariables::PMECH); + + const auto OMEGA = static_cast(HygovExternalVariables::OMEGA); + const auto PREF = static_cast(HygovExternalVariables::PREF); + const auto PAUX = static_cast(HygovExternalVariables::PAUX); + + const ScalarT xn = y[XN]; + const ScalarT xf = y[XF]; + const ScalarT c = y[C]; + const ScalarT g = y[G]; + const ScalarT q = y[Q]; + const ScalarT omegadb = y[OMEGADB]; + const ScalarT ef = y[EF]; + const ScalarT fc = y[FC]; + const ScalarT rc = y[RC]; + const ScalarT pgv = y[PGV]; + const ScalarT head = y[H]; + const ScalarT pmech = y[PMECH]; + + const ScalarT xn_dot = yp[XN]; + const ScalarT xf_dot = yp[XF]; + const ScalarT c_dot = yp[C]; + const ScalarT g_dot = yp[G]; + const ScalarT q_dot = yp[Q]; + + const ScalarT omega = ws[OMEGA]; + const ScalarT pref = ws[PREF]; + const ScalarT paux = ws[PAUX]; + + const ScalarT yomega = xn + leadlag_gain_ * (omegadb - xn); + + f[XN] = -xn_dot + (omegadb - xn) / Tnp_; + f[XF] = -xf_dot + (ef - xf) / Tf_; + f[C] = -c_dot + Math::antiwindup(c, rc, Gmin_response_, Gmax_response_); + f[G] = -g_dot + (c - g) / Tg_; + f[Q] = -q_dot + (Hdam_eff_ - head) / Tw_; + f[OMEGADB] = -omegadb + Math::deadband1(omega, -db1_, db1_); + f[EF] = -ef + toComponentBase(pref + paux) - yomega - Rperm_ * c; + f[FC] = -Rtemp_ * fc + xf / Tr_ + (ef - xf) / Tf_; + f[RC] = -rc + Math::clamp(fc, -Velm_, Velm_); + f[PGV] = -pgv + gatePower(g); + f[H] = -q * q + head * pgv * pgv; + f[PMECH] = -toComponentBase(pmech) + At_ * head * (q - Qnl_) - Dturb_ * omega * g; + + return 0; + } + + // + // Private methods + // + + /** + * @brief Read the parameters out of the model data + * + * Every omitted parameter keeps the default documented in the model + * README. A non-numeric or nonfinite value is counted and reported by + * verify() rather than throwing. Integer JSON values are accepted for + * real parameters. All-zero `Gv` and `Pgv` source points select the + * identity gate curve. + * + * @param[in] data Parameters and monitored-variable selections. + */ + template + void Hygov::initializeParameters(const ModelDataT& data) + { + using Params = typename ModelDataT::Parameters; + + parameter_error_count_ = 0; + + auto load_real = [&](auto key, RealT& target, const char* name) -> bool + { + if (!data.parameters.contains(key)) + { + return false; + } + + 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() << "Hygov: parameter '" << name << "' must be numeric\n"; + ++parameter_error_count_; + return false; + } + + const bool ret = std::isfinite(parsed_value); + if (!ret) + { + Log::error() << "Hygov: parameter '" << name << "' must be finite\n"; + ++parameter_error_count_; + return false; + } + + target = parsed_value; + return true; + }; + + bool ret = load_real(Params::Trate, va_component_base_, "Trate"); + if (ret) + { + ret = va_component_base_ > ZERO; + if (!ret) + { + Log::error() << "Hygov: Trate must be positive when provided\n"; + ++parameter_error_count_; + } + va_component_base_ *= static_cast(1.0e6); + } + load_real(Params::Rperm, Rperm_, "Rperm"); + load_real(Params::Rtemp, Rtemp_, "Rtemp"); + load_real(Params::Tr, Tr_, "Tr"); + load_real(Params::Tf, Tf_, "Tf"); + load_real(Params::Tg, Tg_, "Tg"); + load_real(Params::Velm, Velm_, "Velm"); + load_real(Params::Gmax, Gmax_, "Gmax"); + load_real(Params::Gmin, Gmin_, "Gmin"); + load_real(Params::Tw, Tw_, "Tw"); + load_real(Params::At, At_, "At"); + load_real(Params::Dturb, Dturb_, "Dturb"); + load_real(Params::Qnl, Qnl_, "Qnl"); + load_real(Params::Tn, Tn_, "Tn"); + load_real(Params::Tnp, Tnp_, "Tnp"); + load_real(Params::db1, db1_, "db1"); + if (load_real(Params::db2, db2_, "db2") && db2_ != ZERO) + { + Log::warning() << "Hygov: nonzero db2 requests mechanical backlash, " + "but backlash is not implemented and db2 is ignored\n"; + } + load_real(Params::Hdam, Hdam_, "Hdam"); + load_real(Params::Gv0, Gv_[0], "Gv0"); + load_real(Params::Gv1, Gv_[1], "Gv1"); + load_real(Params::Gv2, Gv_[2], "Gv2"); + load_real(Params::Gv3, Gv_[3], "Gv3"); + load_real(Params::Gv4, Gv_[4], "Gv4"); + load_real(Params::Gv5, Gv_[5], "Gv5"); + load_real(Params::Pgv0, Pgv_[0], "Pgv0"); + load_real(Params::Pgv1, Pgv_[1], "Pgv1"); + load_real(Params::Pgv2, Pgv_[2], "Pgv2"); + load_real(Params::Pgv3, Pgv_[3], "Pgv3"); + load_real(Params::Pgv4, Pgv_[4], "Pgv4"); + load_real(Params::Pgv5, Pgv_[5], "Pgv5"); + + setDerivedParameters(); + } + + /** + * @brief Bind the monitorable variables to their internal states + * + * The mechanical-power output is published on the system base and the + * remaining outputs on the component base, as documented in the model + * README. + */ + template + void Hygov::initializeMonitor() + { + using Variable = typename ModelDataT::MonitorableVariables; + + monitor_->set(Variable::pmech, [this] + { return y_.getData()[static_cast(HygovInternalVariables::PMECH)]; }); + monitor_->set(Variable::filter, [this] + { return y_.getData()[static_cast(HygovInternalVariables::XF)]; }); + monitor_->set(Variable::desiredgate, [this] + { return y_.getData()[static_cast(HygovInternalVariables::C)]; }); + monitor_->set(Variable::gate, [this] + { return y_.getData()[static_cast(HygovInternalVariables::G)]; }); + monitor_->set(Variable::flow, [this] + { return y_.getData()[static_cast(HygovInternalVariables::Q)]; }); + monitor_->set(Variable::head, [this] + { return y_.getData()[static_cast(HygovInternalVariables::H)]; }); + } + + /** + * @brief Resolve the parameter-derived constants + * + * Resolves the default gate curve, floors each governor time constant, + * derives the speed lead-lag gain, and initializes the effective + * response limits and head. + */ + template + void Hygov::setDerivedParameters() + { + // Model data uses an all-exact-zero curve to mean "no curve supplied", + // so this is an exact comparison by intent rather than a tolerance + // test. Any nonzero point selects the given curve. + auto is_nonzero = [](RealT value) + { return value != ZERO; }; + + const bool curve_supplied = + std::any_of(Gv_.begin(), Gv_.end(), is_nonzero) + || std::any_of(Pgv_.begin(), Pgv_.end(), is_nonzero); + if (!curve_supplied) + { + Gv_ = {ZERO, + static_cast(0.2), + static_cast(0.4), + static_cast(0.6), + static_cast(0.8), + ONE}; + Pgv_ = Gv_; + } + + // 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() << "Hygov: " << name << " must be non-negative\n"; + ++parameter_error_count_; + } + }; + + check_non_negative(Tr_, "Tr"); + check_non_negative(Tf_, "Tf"); + check_non_negative(Tg_, "Tg"); + check_non_negative(Tw_, "Tw"); + check_non_negative(Tnp_, "Tnp"); + + if (Tr_ < TIME_CONSTANT_MINIMUM || Tf_ < TIME_CONSTANT_MINIMUM + || Tg_ < TIME_CONSTANT_MINIMUM || Tw_ < TIME_CONSTANT_MINIMUM + || Tnp_ < TIME_CONSTANT_MINIMUM) + { + Log::warning() << "Hygov: Tr, Tf, Tg, Tw, and Tnp below " + << TIME_CONSTANT_MINIMUM + << " s are raised to preserve Hessenberg form\n"; + } + + // HYGOV residuals solve explicitly for the state derivatives to preserve + // Hessenberg form. A zero time constant would instead require an implicit + // residual formulation, so enforce a strictly positive lower bound. + Tr_ = std::max(Tr_, TIME_CONSTANT_MINIMUM); + Tf_ = std::max(Tf_, TIME_CONSTANT_MINIMUM); + Tg_ = std::max(Tg_, TIME_CONSTANT_MINIMUM); + Tw_ = std::max(Tw_, TIME_CONSTANT_MINIMUM); + Tnp_ = std::max(Tnp_, TIME_CONSTANT_MINIMUM); + + leadlag_gain_ = Tn_ / Tnp_; + Gmin_response_ = Gmin_; + Gmax_response_ = Gmax_; + Hdam_eff_ = Hdam_; + } + + /** + * @brief Evaluate the nonlinear gate-to-power curve + * + * Sums the five smooth CommonMath linear segments spanned by the + * `Gv`/`Pgv` points, so the same fixed expression serves the residual + * and both scalar instantiations. + * + * @param[in] gate Gate position. + * @return Turbine power at nominal head. + */ + template + __attribute__((always_inline)) inline scalar_type + Hygov::gatePower(scalar_type gate) const + { + ScalarT retval = Pgv_[0] + + Math::linseg(gate, Gv_[0], Gv_[1], Pgv_[1] - Pgv_[0]) + + Math::linseg(gate, Gv_[1], Gv_[2], Pgv_[2] - Pgv_[1]) + + Math::linseg(gate, Gv_[2], Gv_[3], Pgv_[3] - Pgv_[2]) + + Math::linseg(gate, Gv_[3], Gv_[4], Pgv_[4] - Pgv_[3]) + + Math::linseg(gate, Gv_[4], Gv_[5], Pgv_[5] - Pgv_[4]); + + return retval; + } + + /** + * @brief Steady component-base mechanical power at a gate and dam head + * + * At the steady state the head equals the given dam head and the flow + * follows the gate curve, so the PGV, H, and PMECH rows collapse to + * @f[ + * P_{\mathrm{m}}(g,H_{\mathrm{dam}}) + * = A_t H_{\mathrm{dam}} + * \left(\sqrt{H_{\mathrm{dam}}}\,N_{\mathrm{GV}}(g) + * - q_{\mathrm{NL}}\right). + * @f] + * The expression is composed exactly as those rows compose it, so a + * gate solved against it zeros the implemented residual at machine + * rounding. + * + * @param[in] gate Gate position. + * @param[in] Hdam Dam head. + * @return Steady mechanical power on the component base. + */ + template + typename Hygov::RealT + Hygov::initialMechanicalPower(RealT gate, + RealT Hdam) const + { + const RealT pgv = static_cast(gatePower(static_cast(gate))); + const RealT q = std::sqrt(Hdam) * pgv; + return At_ * Hdam * (q - Qnl_); + } + + /** + * @brief Bisect a bracketed initialization residual + * + * Each iteration replaces one endpoint with a representable midpoint + * strictly inside the interval. A finite floating-point interval has a + * finite number of representable values, so the loop terminates when no + * interior midpoint remains. + * + * @tparam FuncT Residual callable. + * @param[in] a Lower endpoint. + * @param[in] b Upper endpoint. + * @param[in] fa Residual at the lower endpoint. + * @param[in] fb Residual at the upper endpoint. + * @param[in] residual Residual callable. + * @pre The endpoints are finite, @f$a < b@f$, and @f$f(a) \le 0 \le f(b)@f$. + * @return The endpoint with the smaller residual magnitude, or a quiet + * NaN if an interior residual is NaN. + */ + template + template + typename Hygov::RealT + Hygov::bisectInitialRoot(RealT a, + RealT b, + RealT fa, + RealT fb, + FuncT residual) + { + while (true) + { + const RealT mid = std::midpoint(a, b); + if (mid <= a || b <= mid) + { + break; + } + + const RealT fmid = residual(mid); + if (std::isnan(fmid)) + { + return std::numeric_limits::quiet_NaN(); + } + if (fmid <= ZERO) + { + a = mid; + fa = fmid; + } + else + { + b = mid; + fb = fmid; + } + } + + if (std::abs(fa) <= std::abs(fb)) + { + return a; + } + return b; + } + + /** + * @brief Solve the steady gate position for a given mechanical power + * + * Initialization requires a zero speed deviation and verify() requires + * the steady power to rise across [Gv0, Gv5], so the full gate-curve + * endpoint residuals decide feasibility and bisection converges to a + * root of the nondecreasing steady-power curve at the configured dam + * head. + * + * @pre verify() reports no errors. + * + * @param[in] pmech Mechanical power on the component base. + * @return The gate position, or a quiet NaN when no gate inside + * [Gv0, Gv5] reproduces the value within the initialization + * tolerance. + * + * @warning This function contains conditional branching and may be used + * during initialization, but not during residual evaluation. + */ + template + typename Hygov::RealT + Hygov::solveInitialGate(RealT pmech) const + { + const auto residual = [this, pmech](RealT gate) + { return initialMechanicalPower(gate, Hdam_) - pmech; }; + + RealT a = Gv_[0]; + RealT b = Gv_[5]; + RealT fa = residual(a); + RealT fb = residual(b); + + // A value just below the achievable range pins to the first gate + // point when it is within the initialization tolerance of the edge. + if (fa > ZERO) + { + if (fa <= INITIALIZATION_TOLERANCE) + { + return a; + } + return std::numeric_limits::quiet_NaN(); + } + if (fb < ZERO) + { + return std::numeric_limits::quiet_NaN(); + } + + return bisectInitialRoot(a, b, fa, fb, residual); + } + + /** + * @brief Solve the effective dam head for a high initial power + * + * Starting from the configured dam head, brackets a higher value whose + * steady mechanical power at Gv5 reaches the initial value, then + * bisects to machine rounding. + * + * @pre verify() reports no errors. + * @pre The initial mechanical power exceeds the last gate-curve endpoint. + * + * @param[in] pmech Mechanical power on the component base. + * @return Effective dam head, or a quiet NaN when no finite value is found. + * + * @warning This function contains conditional branching and may be used + * during initialization, but not during residual evaluation. + */ + template + typename Hygov::RealT + Hygov::solveInitialDamHead(RealT pmech) const + { + const RealT nan = std::numeric_limits::quiet_NaN(); + const auto residual = [this, pmech](RealT Hdam) + { return initialMechanicalPower(Gv_[5], Hdam) - pmech; }; + + RealT a = Hdam_; + RealT b = Hdam_; + RealT fa = residual(a); + RealT fb = fa; + + if (!std::isfinite(fa) || !(fa < ZERO) ) + { + return nan; + } + + const RealT maximum_head = std::numeric_limits::max(); + while (fb < ZERO) + { + a = b; + fa = fb; + + if (b >= maximum_head) + { + return nan; + } + if (b > maximum_head / TWO) + { + b = maximum_head; + } + else + { + b *= TWO; + } + + fb = residual(b); + if (std::isnan(fb)) + { + return nan; + } + } + + return bisectInitialRoot(a, b, fa, fb, residual); + } + + /** + * @brief Convert a system-base power to HYGOV component base + * + * @param[in] value Quantity on the system base. + * @return The same quantity on the component base. + */ + template + scalar_type Hygov::toComponentBase(scalar_type value) const + { + return value * va_system_base_ / va_component_base_; + } + + /** + * @brief Convert a component-base power to the system base + * + * @param[in] value Quantity on the component base. + * @return The same quantity on the system base. + */ + template + scalar_type Hygov::toSystemBase(scalar_type value) const + { + return value / toComponentBase(static_cast(ONE)); + } + + } // namespace Governor + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Governor/HYGOV/README.md b/GridKit/Model/PhasorDynamics/Governor/HYGOV/README.md new file mode 100644 index 000000000..31330168e --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Governor/HYGOV/README.md @@ -0,0 +1,414 @@ +# **Hydro Turbine-Governor Model (HYGOV)** + +HYGOV is a hydro turbine-governor model with temporary droop, a gate servo, and +a nonlinear single-penstock turbine. + +## Notes + +- HYGOVD `dbL`/`dbH`, `db2` backlash, and Kaplan blade-servo fields are not + modeled. The `db2` JSON field is accepted for source-format compatibility. + A nonzero value logs a warning and is ignored. + +## Block Diagram + +![HYGOV governor block diagram](../../../../../docs/Figures/PhasorDynamics/HYGOV/diagram.png) + +Figure 1: HYGOV governor model. Figure courtesy of the +[PowerWorld HYGOV model reference](https://www.powerworld.com/WebHelp/Content/TransientModels_HTML/Governor%20HYGOV%20and%20HYGOVD.htm). + +## Model Parameters + +Symbol | Units | JSON | Description | Typical Value | Note +------------------------|----------|---------------|------------------------------------------|---------------|------ +$T^\mathrm{rate}$ | [MW] | `Trate` | Turbine-rating power base | 100.0 | System power base when omitted +$R_{\mathrm{perm}}$ | [p.u.] | `Rperm` | Permanent droop | 0.04 | Source label: `R` +$R_{\mathrm{temp}}$ | [p.u.] | `Rtemp` | Temporary droop | 0.3 | Source label: `r` +$T_r$ | [sec] | `Tr` | Temporary-droop reset time constant | 5.0 | +$T_f$ | [sec] | `Tf` | Governor error filter time constant | 0.05 | +$T_g$ | [sec] | `Tg` | Gate servo time constant | 0.5 | +$V_{\mathrm{elm}}$ | [p.u./s] | `Velm` | Maximum desired-gate velocity magnitude | 0.2 | +$G^{\max}$ | [p.u.] | `Gmax` | Configured upper gate response limit | 1.0 | +$G^{\min}$ | [p.u.] | `Gmin` | Configured lower gate response limit | 0.0 | +$T_w$ | [sec] | `Tw` | Water inertia time constant | 1.0 | +$A_t$ | [p.u.] | `At` | Turbine gain | 1.2 | +$D_{\mathrm{turb}}$ | [p.u.] | `Dturb` | Turbine damping coefficient | 0.5 | +$q_{\mathrm{NL}}$ | [p.u.] | `Qnl` | No-load flow at nominal head | 0.05 | +$T_n$ | [sec] | `Tn` | Speed lead-lag numerator time constant | 0.0 | +$T_{\mathrm{np}}$ | [sec] | `Tnp` | Speed lead-lag denominator time constant | 0.0 | +$D_{\omega}$ | [p.u.] | `db1` | Type 1 speed deadband threshold | 0.0 | +$D_2$ | [p.u.] | `db2` | Unsupported mechanical backlash deadband | 0.0 | Nonzero values warn and are ignored +$H_{\mathrm{dam}}$ | [p.u.] | `Hdam` | Configured dam head | 1.0 | Lower bound on effective head +$G_V^{(k)}$ | [p.u.] | `Gv0`-`Gv5` | Gate point $k$ of the gain curve | 0.0 | $k=0,\ldots,5$ +$P_{\mathrm{GV}}^{(k)}$ | [p.u.] | `Pgv0`-`Pgv5` | Power point $k$ of the gain curve | 0.0 | $k=0,\ldots,5$ + +Every parameter is optional. Real-valued parameters accept real or integer +JSON values. All-zero `Gv` and `Pgv` source points select the identity curve. + +### Parameter Validation + +Real-valued parameters, `Known` initial values, power bases, and base-conversion +ratios must be finite. The bases and ratios must also be positive. Invalid +HYGOV parameter sets are rejected by the following checks: + +```math +\begin{aligned} + T^\mathrm{rate} &> 0 \quad \text{when provided} \\ + T_r, T_f, T_g, T_w, T_{\mathrm{np}} + &\ge 0 \\ + R_{\mathrm{temp}} + &\ne 0 \\ + T_n + &\ge 0 \\ + V_{\mathrm{elm}} + &\ge 0 \\ + G^{\min} + &< G^{\max} \\ + A_t + &> 0 \\ + D_{\mathrm{turb}} + &\ge 0 \\ + D_{\omega} + &\ge 0 \\ + H_{\mathrm{dam}} + &> 0 \\ + G_V^{(k)} + &< G_V^{(k+1)} + \quad k\in\{0,\ldots,4\} \\ + P_{\mathrm{GV}}^{(k)} + &\le P_{\mathrm{GV}}^{(k+1)} + \quad k\in\{0,\ldots,4\} \\ + G_V^{(0)} \le G^{\min} + &< G^{\max} \le G_V^{(5)} \\ + P_{\mathrm{m}}(G_V^{(5)}) - P_{\mathrm{m}}(G_V^{(0)}) + &> \epsilon_{\mathrm{init}} +\end{aligned} +``` + +The final condition uses the steady mechanical power and tolerance defined +under [Internal Initialization](#internal-initialization). + +### 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,f,g,w,\mathrm{np}\} \\ + k_{\mathrm{base}} + &= \dfrac{S^\mathrm{sys}}{T^\mathrm{rate}} \\ + k_n + &= \dfrac{T_n}{T_{\mathrm{np}}} \\ + N_{\mathrm{GV}}(x) + &= + P_{\mathrm{GV}}^{(0)} + + \sum_{k\in\{0,\ldots,4\}} + \text{linseg}\!\left( + x;\, + G_V^{(k)},\, + G_V^{(k+1)},\, + P_{\mathrm{GV}}^{(k+1)} - P_{\mathrm{GV}}^{(k)} + \right) +\end{aligned} +``` + +Multiplying by $k_\mathrm{base}$ converts system base to component base. + +CommonMath defines the [`linseg`](../../../../CommonMath.md#linseg) helper +used by $N_{\mathrm{GV}}$. + +## Model Ports + +Name | Port | Init | Description +--------|--------|---------|------ +`speed` | Input | Known | Machine speed deviation +`pref` | Input | Unknown | Active-power/load reference +`paux` | Input | Known | Auxiliary power input +`pmech` | Output | Known | Mechanical power 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 +`pmech` output must be assigned. The signal inputs are optional. Unattached +`speed` and `paux` inputs default to zero. + +## Model Variables + +### Internal Variables + +#### Differential + +Symbol | Units | Description | Note +------------------------|--------|-------------------------------------|------ +$x_n$ | [p.u.] | Speed lead-lag denominator state | Not circled in Fig. 1. Realizes the `Tn`/`Tnp` block +$x_f$ | [p.u.] | Governor error filter output | State 1 in Fig. 1 +$c$ | [p.u.] | Desired-gate position | State 2 in Fig. 1 +$g$ | [p.u.] | Gate position | State 3 in Fig. 1 +$q$ | [p.u.] | Turbine flow | State 4 in Fig. 1 + +#### Algebraic + +Symbol | Units | Description | Note +------------------------|----------|---------------------------------------------|------ +$\omega_{\mathrm{db}}$ | [p.u.] | Type 1 deadbanded speed deviation | +$e_f$ | [p.u.] | Governor error into the filter | Reference path less conditioned speed and permanent-droop feedback +$f_c$ | [p.u./s] | Desired-gate derivative target | Before rate and position limits +$r_c$ | [p.u./s] | Rate-limited desired-gate derivative target | Limited by $\pm V_{\mathrm{elm}}$ +$P_{\mathrm{GV}}$ | [p.u.] | Nonlinear gate-to-power curve output | $N_{\mathrm{GV}}(g)$ +$H$ | [p.u.] | Turbine head | Implicit water-column head +$P_{\mathrm{m}}$ | [p.u.] | Mechanical power to generator | System base + +### External Variables + +#### Differential + +None. + +#### Algebraic + +Symbol | Units | Init | Description | Note +------------------|--------|---------|-----------------------------|------ +$\omega$ | [p.u.] | Known | Machine speed deviation | Optional signal port `speed`. Defaults to zero +$P^\mathrm{ref}$ | [p.u.] | Unknown | Active-power/load reference | Optional signal port `pref`, system base +$P^\mathrm{aux}$ | [p.u.] | Known | Auxiliary power input | Optional signal port `paux`, system base, defaults to zero + +## Model Equations + +### Differential Equations + +The effective desired-gate response limits +$G_{\mathrm{resp}}^{\min}$ and $G_{\mathrm{resp}}^{\max}$ and the effective +dam head $H_{\mathrm{dam}}^{\mathrm{eff}}$ are resolved during initialization. + +```math +\begin{aligned} + 0 &= + -\dot{x}_n + + \dfrac{1}{T_{\mathrm{np}}} + \left(\omega_{\mathrm{db}} - x_n\right) \\ + 0 &= + -\dot{x}_f + + \dfrac{1}{T_f} + \left(e_f - x_f\right) \\ + 0 &= + -\dot{c} + + \text{antiwindup} + \left(c, r_c;\, G_{\mathrm{resp}}^{\min}, + G_{\mathrm{resp}}^{\max}\right) \\ + 0 &= + -\dot{g} + + \dfrac{1}{T_g} + \left(c - g\right) \\ + 0 &= + -\dot{q} + + \dfrac{1}{T_w} + \left(H_{\mathrm{dam}}^{\mathrm{eff}} - H\right) +\end{aligned} +``` + +CommonMath defines the [`antiwindup`](../../../../CommonMath.md#antiwindup) +target and smooth approximation. + +### Algebraic Equations + +```math +\begin{aligned} + 0 &= + -\omega_{\mathrm{db}} + + \text{deadband1} + \left(\omega;\, -D_{\omega}, D_{\omega}\right) \\ + 0 &= + -e_f + + k_{\mathrm{base}}\left(P^\mathrm{ref} + P^\mathrm{aux}\right) + - x_n + - k_n\left(\omega_{\mathrm{db}} - x_n\right) + - R_{\mathrm{perm}}c \\ + 0 &= + -R_{\mathrm{temp}}f_c + + \dfrac{x_f}{T_r} + + \dfrac{e_f - x_f}{T_f} \\ + 0 &= + -r_c + + \text{clamp} + \left(f_c;\, -V_{\mathrm{elm}}, V_{\mathrm{elm}}\right) \\ + 0 &= + -P_{\mathrm{GV}} + + N_{\mathrm{GV}}(g) \\ + 0 &= + -q^2 + + H P_{\mathrm{GV}}^2 \\ + 0 &= + -k_{\mathrm{base}}P_{\mathrm{m}} + + A_t H\left(q - q_{\mathrm{NL}}\right) + - D_{\mathrm{turb}}\omega g +\end{aligned} +``` + +CommonMath defines helper targets and smooth approximations for +[deadband1 and clamp](../../../../CommonMath.md#derived-functions). + +## Initialization + +### Input Initialization + +```math +\begin{aligned} + \omega + &\leftarrow \text{machine speed deviation} \\ + P_{\mathrm{m}} + &\leftarrow \text{machine mechanical power on system base} \\ + P^\mathrm{aux} + &\leftarrow \text{auxiliary power input on system base} +\end{aligned} +``` + +Initialization never replaces the system-base value held in $P_{\mathrm{m}}$. + +### Internal Initialization + +Initialization requires an exactly zero speed deviation, $\omega = 0$. +Restart initialization of a moving machine is not supported. All internal +derivatives are set to zero. + +Initialization first solves the gate at the configured dam head over the full +$[G_V^{(0)},G_V^{(5)}]$ gate curve. If that gate lies outside the configured +$[G^{\min},G^{\max}]$ interval, the corresponding response limit is expanded +to include it. The configured parameters are unchanged. + +If the required mechanical power exceeds the value at $G_V^{(5)}$, the gate is +pinned there and an effective dam head +$H_{\mathrm{dam}}^{\mathrm{eff}} \ge H_{\mathrm{dam}}$ is raised to reproduce +the operating point. Both searches use the same smooth $N_{\mathrm{GV}}$ curve +as the residual. No upper limit is applied to the head adjustment. The +effective values remain the response limits and water-column setpoint during +simulation. + +```math +\begin{aligned} + H + &\leftarrow H_{\mathrm{dam}}^{\mathrm{eff}} \\ + g + &\leftarrow \text{gate in } [G_V^{(0)},G_V^{(5)}] \text{ satisfying} \\ + &\qquad k_{\mathrm{base}}P_{\mathrm{m}} + = A_t H\left(\sqrt{H}\,N_{\mathrm{GV}}(g) - q_{\mathrm{NL}}\right) \\ + G_{\mathrm{resp}}^{\min} + &\leftarrow \min\!\left(G^{\min},g\right) \\ + G_{\mathrm{resp}}^{\max} + &\leftarrow \max\!\left(G^{\max},g\right) \\ + P_{\mathrm{GV}} + &\leftarrow N_{\mathrm{GV}}(g) \\ + q + &\leftarrow \sqrt{H}\,P_{\mathrm{GV}} \\ + c + &\leftarrow g \\ + \omega_{\mathrm{db}} + &\leftarrow \text{deadband1}\!\left(\omega;\, -D_{\omega}, D_{\omega}\right) \\ + x_n + &\leftarrow \omega_{\mathrm{db}} \\ + x_f + &\leftarrow 0 \\ + e_f + &\leftarrow 0 \\ + f_c + &\leftarrow 0 \\ + r_c + &\leftarrow 0 +\end{aligned} +``` + +A value within $\epsilon_{\mathrm{init}} = 100\,\epsilon_{\mathrm{mach}}$ +below the $G_V^{(0)}$ endpoint initializes at $G_V^{(0)}$ with a +mechanical-power residual up to $\epsilon_{\mathrm{init}}$. A lower value, or +a high-side value without a finite effective dam head, is rejected. All other +accepted values initialize with every residual at machine rounding. + +Every check resolves before state, the effective response limits, the effective +dam head, or signals are written, so a rejected initialization leaves them +unchanged. + +### Output Initialization + +```math +\begin{aligned} + P^\mathrm{ref} + &\leftarrow + \dfrac{1}{k_{\mathrm{base}}} + \left[ + e_f + - k_{\mathrm{base}}P^\mathrm{aux} + + x_n + + k_n\left(\omega_{\mathrm{db}} - x_n\right) + + R_{\mathrm{perm}}c + \right] +\end{aligned} +``` + +## Monitorable Outputs + +Output | Units | Description | Note +---------------|--------|------------------------------|------ +`pmech` | [p.u.] | Mechanical-power output | $P_{\mathrm{m}}$ (system base) +`filter` | [p.u.] | Governor error filter output | $x_f$ (component base) +`desiredgate` | [p.u.] | Desired-gate position | $c$ (component base) +`gate` | [p.u.] | Gate position | $g$ (component base) +`flow` | [p.u.] | Turbine flow | $q$ (component base) +`head` | [p.u.] | Turbine head | $H$ (component base) + +## Testing + +- `validation()` checks construction, monitor creation, parameter validation, + signal configuration, and minimum time-constant handling. +- `initializationAndSignals()` checks initialization, base conversion, + signal publication, monitor output, and unattached-reference latching. +- `initializationDomain()` checks effective-limit and effective-head + initialization, rejection atomicity, and initialization boundaries. +- `initializationExactness()` checks that initialized steady residuals rest + at machine rounding across the gate curve. +- `residualEquations()` checks every model residual against a fixed + numerical answer key. +- `governorControl()` checks the speed deadband, the desired-gate velocity + limit, and the gate-position anti-windup. +- `turbineDynamics()` checks the gate-power curve, the water column, turbine + damping, and initialization through the nonlinear curve. +- `jacobian()` compares the dependency-tracking and Enzyme Jacobians across + the gate curve when Enzyme support is enabled. + +## Appendix A: Backlash + +Input $u$, output $y$, half-play $b$, with $|u - y| \le b$. + +```math +\begin{aligned} + \dot{y} + &= + \begin{cases} + \dot{u} & |u - y| = b \text{ and } \dot{u}\left(u - y\right) > 0 \\ + 0 & \text{otherwise} + \end{cases} +\end{aligned} +``` + +which can be written in terms of our smooth functions as + +```math +\begin{aligned} + 0 &= + -\dot{y} + + \text{ramp}(\dot{u})\,\text{above}(u - y;\, b) + - \text{ramp}(-\dot{u})\,\text{below}(u - y;\, -b) +\end{aligned} +``` + +CommonMath defines the [`ramp`](GridKit/CommonMath.md#-ramp), +[`above`](GridKit/CommonMath.md#above), and +[`below`](GridKit/CommonMath.md#below) targets and smooth approximations. This is deferred until we permit non Hessenberg forms. Once permitted we should define: + +```math +\begin{aligned} + \text{backlash}(u,\dot{u},y;b) &= + \text{ramp}(\dot{u})\,\text{above}(u - y;\, b) + - \text{ramp}(-\dot{u})\,\text{below}(u - y;\, -b) +\end{aligned} +``` diff --git a/GridKit/Model/PhasorDynamics/Governor/README.md b/GridKit/Model/PhasorDynamics/Governor/README.md index fc46e5eb3..6d3983adb 100644 --- a/GridKit/Model/PhasorDynamics/Governor/README.md +++ b/GridKit/Model/PhasorDynamics/Governor/README.md @@ -9,5 +9,6 @@ A governor models the control system that regulates the output power of a machin There are a few standard Governor models - Turbine Governor (See [TGOV1](Tgov1/README.md)) +- Hydro Turbine Governor (See [HYGOV](HYGOV/README.md)) - IEEE Type G1 Turbine Governor (See [IEEEG1](IEEEG1/README.md)) - General Governor (See [GGOV1](GGOV1/README.md)) diff --git a/GridKit/Model/PhasorDynamics/INPUT_FORMAT.md b/GridKit/Model/PhasorDynamics/INPUT_FORMAT.md index 0ebbcd7ba..7a8fa915f 100644 --- a/GridKit/Model/PhasorDynamics/INPUT_FORMAT.md +++ b/GridKit/Model/PhasorDynamics/INPUT_FORMAT.md @@ -153,6 +153,7 @@ are specified: [GenClassical](SynchronousMachine/GenClassical/README.md) | the classical machine model [Regca](Converter/REGCA/README.md) | WECC REGCA renewable generator/converter model [Tgov1](Governor/Tgov1/README.md) | the TGOV1 governor model + [Hygov](Governor/HYGOV/README.md) | the HYGOV hydro turbine-governor model [Ieeet1](Exciter/IEEET1/README.md) | the IEEET1 exciter model [Esdc1a](Exciter/ESDC1A/README.md) | the ESDC1A exciter model [SexsPti](Exciter/SEXS-PTI/README.md) | the SEXS-PTI simplified exciter model diff --git a/GridKit/Model/PhasorDynamics/SystemModelData.hpp b/GridKit/Model/PhasorDynamics/SystemModelData.hpp index d743ab63e..e82fc965e 100644 --- a/GridKit/Model/PhasorDynamics/SystemModelData.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModelData.hpp @@ -14,6 +14,7 @@ #include #include #include +#include #include #include #include @@ -46,6 +47,7 @@ namespace GridKit using RegcaDataT = Converter::RegcaData; using Tgov1DataT = Governor::Tgov1Data; using Esdc1aDataT = Exciter::Esdc1aData; + using HygovDataT = Governor::HygovData; using Ieeet1DataT = Exciter::Ieeet1Data; using SexsPtiDataT = Exciter::SexsPtiData; using IeeestDataT = Stabilizer::IeeestData; @@ -107,6 +109,7 @@ namespace GridKit 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 hygov; ///< HYGOV governors 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 dc9a7e74c..328053288 100644 --- a/GridKit/Model/PhasorDynamics/SystemModelDataJSONParser.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModelDataJSONParser.hpp @@ -147,6 +147,12 @@ namespace GridKit raw_component.get_to(gov); sm.gov.push_back(gov); } + else if (kind == "Hygov") + { + typename SystemModelData::HygovDataT hygov; + raw_component.get_to(hygov); + sm.hygov.push_back(hygov); + } else if (kind == "Ieeet1") { typename SystemModelData::Ieeet1DataT exciter; diff --git a/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp b/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp index d4c126cd7..ceaa8ed6d 100644 --- a/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp @@ -315,6 +315,42 @@ namespace GridKit addComponent(gov); } + // Add HYGOV governors + for (const auto& hygovdata : data.hygov) + { + auto* hygov = new Hygov(hygovdata); + + if (hygovdata.signal_inputs.contains(HygovSignalInputs::speed)) + { + IdxT speed = hygovdata.signal_inputs.at(HygovSignalInputs::speed); + constexpr auto OMEGA = HygovExternalVariables::OMEGA; + hygov->getSignals().template attachSignalNode(getSignal(speed)); + } + + if (hygovdata.signal_inputs.contains(HygovSignalInputs::pref)) + { + IdxT pref = hygovdata.signal_inputs.at(HygovSignalInputs::pref); + constexpr auto PREF = HygovExternalVariables::PREF; + hygov->getSignals().template attachSignalNode(getSignal(pref)); + } + + if (hygovdata.signal_inputs.contains(HygovSignalInputs::paux)) + { + IdxT paux = hygovdata.signal_inputs.at(HygovSignalInputs::paux); + constexpr auto PAUX = HygovExternalVariables::PAUX; + hygov->getSignals().template attachSignalNode(getSignal(paux)); + } + + if (hygovdata.signal_outputs.contains(HygovSignalOutputs::pmech)) + { + IdxT pmech = hygovdata.signal_outputs.at(HygovSignalOutputs::pmech); + constexpr auto PMECH = HygovInternalVariables::PMECH; + hygov->getSignals().template assignSignalNode(getSignal(pmech)); + } + + addComponent(hygov); + } + for (const auto& excitedata : data.exciter) { IdxT bus_index = 0; @@ -714,20 +750,22 @@ namespace GridKit template int SystemModel::initialize() { + int status = 0; + for (const auto& bus : buses_) { - bus->initialize(); + status += bus->initialize(); } for (const auto& component : components_) { - component->initialize(); + status += component->initialize(); } y_.setDataUpdated(); yp_.setDataUpdated(); - return 0; + return status; } /** diff --git a/docs/Figures/PhasorDynamics/HYGOV/diagram.png b/docs/Figures/PhasorDynamics/HYGOV/diagram.png new file mode 100644 index 000000000..fc7d7916d Binary files /dev/null and b/docs/Figures/PhasorDynamics/HYGOV/diagram.png differ diff --git a/docs/GridKit/Model/PhasorDynamics/Governor/HYGOV/README.md b/docs/GridKit/Model/PhasorDynamics/Governor/HYGOV/README.md new file mode 100644 index 000000000..0e4dd411a --- /dev/null +++ b/docs/GridKit/Model/PhasorDynamics/Governor/HYGOV/README.md @@ -0,0 +1,6 @@ +# HYGOV + +```{include} ../../../../../../GridKit/Model/PhasorDynamics/Governor/HYGOV/README.md +:start-line: 1 +:relative-images: +``` diff --git a/docs/GridKit/Model/PhasorDynamics/Governor/README.md b/docs/GridKit/Model/PhasorDynamics/Governor/README.md index fe3d5b5b7..45b32ae26 100644 --- a/docs/GridKit/Model/PhasorDynamics/Governor/README.md +++ b/docs/GridKit/Model/PhasorDynamics/Governor/README.md @@ -6,6 +6,7 @@ :hidden: TGOV1 +HYGOV IEEEG1 GGOV1 ``` diff --git a/tests/UnitTests/PhasorDynamics/CMakeLists.txt b/tests/UnitTests/PhasorDynamics/CMakeLists.txt index 87950ac62..122e925ee 100644 --- a/tests/UnitTests/PhasorDynamics/CMakeLists.txt +++ b/tests/UnitTests/PhasorDynamics/CMakeLists.txt @@ -88,6 +88,14 @@ target_link_libraries( GridKit::phasor_dynamics_bus_dependency_tracking GridKit::testing) +add_executable(test_phasor_governor_hygov runGovernorHygovTests.cpp) +target_link_libraries( + test_phasor_governor_hygov + GridKit::definitions + GridKit::phasor_dynamics_governor_hygov + GridKit::phasor_dynamics_governor_hygov_dependency_tracking + GridKit::testing) + add_executable(test_phasor_exciter_ieeet1 runExciterIeeet1Tests.cpp) target_link_libraries( test_phasor_exciter_ieeet1 @@ -173,6 +181,7 @@ add_test(NAME PhasorDynamicsBusToSignalAdapterTest COMMAND test_phasor_bustosign 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 PhasorDynamicsGovernorHygovTest COMMAND test_phasor_governor_hygov) 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) @@ -198,6 +207,7 @@ install( test_phasor_loadzip test_phasor_genrou test_phasor_governor_tgov1 + test_phasor_governor_hygov test_phasor_exciter_ieeet1 test_phasor_exciter_esdc1a test_phasor_gensal diff --git a/tests/UnitTests/PhasorDynamics/ComponentConnectionTests.hpp b/tests/UnitTests/PhasorDynamics/ComponentConnectionTests.hpp index b973b830e..3b50a9f0b 100644 --- a/tests/UnitTests/PhasorDynamics/ComponentConnectionTests.hpp +++ b/tests/UnitTests/PhasorDynamics/ComponentConnectionTests.hpp @@ -6,6 +6,8 @@ #include #include #include +#include +#include #include #include #include @@ -80,6 +82,55 @@ namespace GridKit return success.report(__func__); } + + /// GENROU initializes first and writes the mechanical power it needs to + /// the shared node. HYGOV then initializes around that value and must + /// leave it unchanged at a steady state. The speed port is left + /// unattached, so HYGOV reads the exactly zero deviation its + /// initialization requires. + TestOutcome genrouHygov() + { + using MachineExternal = PhasorDynamics::GenrouExternalVariables; + using GovernorInternal = PhasorDynamics::Governor::HygovInternalVariables; + using GovernorParams = PhasorDynamics::Governor::HygovParameters; + + TestStatus success = true; + + PhasorDynamics::SystemModel system; + PhasorDynamics::BusInfinite bus( + static_cast(1.0), + static_cast(0.0)); + PhasorDynamics::SignalNode pmech; + PhasorDynamics::Genrou machine(&bus); + + PhasorDynamics::Governor::HygovData governor_data; + governor_data.parameters[GovernorParams::Tnp] = static_cast(1.0); + + PhasorDynamics::Governor::Hygov governor(governor_data); + + machine.getSignals().template attachSignalNode(&pmech); + governor.getSignals().template assignSignalNode(&pmech); + + system.addBus(&bus); + system.addComponent(&machine); + system.addComponent(&governor); + + success *= system.allocate() == 0; + success *= pmech.linked(); + success *= system.initialize() == 0; + success *= system.evaluateResidual() == 0; + + // At zero machine power the required mechanical power is zero. + success *= isEqual(pmech.read(), static_cast(0.0), kTol); + + const auto* residual = governor.getResidual().getData(); + for (IdxT row = 0; row < governor.size(); ++row) + { + success *= isEqual(residual[row], static_cast(0.0), kTol); + } + + return success.report(__func__); + } }; } // namespace Testing diff --git a/tests/UnitTests/PhasorDynamics/ConverterRegcaTests.hpp b/tests/UnitTests/PhasorDynamics/ConverterRegcaTests.hpp index 57b477298..427648e1d 100644 --- a/tests/UnitTests/PhasorDynamics/ConverterRegcaTests.hpp +++ b/tests/UnitTests/PhasorDynamics/ConverterRegcaTests.hpp @@ -719,7 +719,9 @@ namespace GridKit for (size_t i = 0; i < nrows; ++i) { - if (!isEqual(dependency_tracking_jacobian[i], enzyme_jacobian[i])) + if (!isEqual(dependency_tracking_jacobian[i], + enzyme_jacobian[i], + kTol)) { std::cout << "Jacobian row " << i << " mismatch between dependency tracking and Enzyme" @@ -999,7 +1001,7 @@ namespace GridKit bool scalarMatches(ScalarT actual, ScalarT expected, const char* label) const { - if (isEqual(actual, expected)) + if (isEqual(actual, expected, kTol)) { return true; } diff --git a/tests/UnitTests/PhasorDynamics/GovernorHygovTests.hpp b/tests/UnitTests/PhasorDynamics/GovernorHygovTests.hpp new file mode 100644 index 000000000..333ff02d9 --- /dev/null +++ b/tests/UnitTests/PhasorDynamics/GovernorHygovTests.hpp @@ -0,0 +1,1642 @@ +#pragma once + +#include +#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 GovernorHygovTests + { + public: + using ScalarT = scalar_type; + using IdxT = index_type; + using RealT = typename PhasorDynamics::Component::RealT; + + GovernorHygovTests() = default; + ~GovernorHygovTests() = default; + + static constexpr RealT kTol = + static_cast(100.0) * std::numeric_limits::epsilon(); + + /// Construction and every verify() error class, including parameter + /// types and finiteness, parameter relationships, power bases, curve + /// shape, gate-limit domain, the required pmech assignment, and signal + /// linkage, plus differentiability tagging. + TestOutcome validation() + { + TestStatus success = true; + + PhasorDynamics::Governor::Hygov empty; + success *= (empty.size() == static_cast(Internal::MAXIMUM)); + success *= (empty.getMonitor() == nullptr); + + Fixture configured(makeData()); + success *= (configured.hygov.size() == static_cast(Internal::MAXIMUM)); + success *= (configured.hygov.getMonitor() != nullptr); + success *= (configured.hygov.verify() == 0); + + noteExpectedLogs("Testing HYGOV defaults and invalid configurations. " + "Logged errors, time-constant warnings, and an unsupported " + "backlash warning are expected."); + + Fixture minimal(makeMinimalData()); + success *= (minimal.hygov.verify() == 0); + success *= defaultsMatchDocumentedValues(); + + success *= (empty.verify() > 0); + + const RealT nan = std::numeric_limits::quiet_NaN(); + const RealT infinity = std::numeric_limits::infinity(); + + const std::array real_parameters{{ + Params::Trate, + Params::Rperm, + Params::Rtemp, + Params::Tr, + Params::Tf, + Params::Tg, + Params::Velm, + Params::Gmax, + Params::Gmin, + Params::Tw, + Params::At, + Params::Dturb, + Params::Qnl, + Params::Tn, + Params::Tnp, + Params::db1, + Params::db2, + Params::Hdam, + Params::Gv0, + Params::Gv1, + Params::Gv2, + Params::Gv3, + Params::Gv4, + Params::Gv5, + Params::Pgv0, + Params::Pgv1, + Params::Pgv2, + Params::Pgv3, + Params::Pgv4, + Params::Pgv5, + }}; + const std::array nonfinite_values{{nan, infinity, -infinity}}; + + for (const Params parameter : real_parameters) + { + for (const RealT value : nonfinite_values) + { + Fixture invalid_fixture(makeData(), {{parameter, value}}); + success *= (invalid_fixture.hygov.verify() > 0); + } + } + + // The pmech output is required, so a model without an assigned node + // is rejected even when every parameter is valid. + PhasorDynamics::Governor::Hygov unassigned(makeData()); + success *= (unassigned.verify() > 0); + + const std::array, 19> invalid_parameter_values{{ + {Params::Trate, 0.0}, + {Params::Trate, -1.0}, + {Params::Rtemp, 0.0}, + {Params::Tr, -0.1}, + {Params::Tf, -0.1}, + {Params::Tg, -0.1}, + {Params::Tw, -0.1}, + {Params::Tn, -0.1}, + {Params::Tnp, -0.1}, + {Params::Velm, -0.1}, + {Params::Gmin, 1.1}, + {Params::At, 0.0}, + {Params::Dturb, -0.1}, + {Params::db1, -0.1}, + {Params::Hdam, 0.0}, + {Params::Gv2, 0.1}, + {Params::Pgv2, 0.1}, + {Params::Gmin, -0.05}, + {Params::Gmax, 1.05}, + }}; + + for (const auto& [parameter, value] : invalid_parameter_values) + { + Fixture invalid_fixture(makeData(), {{parameter, value}}); + success *= (invalid_fixture.hygov.verify() > 0); + } + + // A curve with no rise cannot yield a unique gate. + Fixture flat_curve(makeData(), + {{Params::Pgv1, 0.0}, + {Params::Pgv2, 0.0}, + {Params::Pgv3, 0.0}, + {Params::Pgv4, 0.0}, + {Params::Pgv5, 0.0}}); + success *= (flat_curve.hygov.verify() > 0); + + // A curve that rises only outside the configured response limits is + // valid because initialization may expand those limits. + Fixture flat_configured_range( + makeData(), + {{Params::Gmin, 0.0}, + {Params::Gmax, 0.2}, + {Params::Pgv0, 0.5}, + {Params::Pgv1, 0.5}, + {Params::Pgv2, 0.5}, + {Params::Pgv3, 0.5}, + {Params::Pgv4, 0.5}, + {Params::Pgv5, 1.0}}); + success *= (flat_configured_range.hygov.verify() == 0); + + // A requested backlash is accepted, warns, and remains inactive. + Fixture backlash(makeData(), {{Params::db2, 0.5}}); + success *= (backlash.hygov.verify() == 0); + + // Integer JSON values are accepted for real parameters; booleans are + // not numeric. + auto integer_real = makeData(); + integer_real.parameters[Params::Tw] = static_cast(2); + Fixture integer_model(integer_real); + success *= (integer_model.hygov.verify() == 0); + + auto bad_numeric_type = makeData(); + bad_numeric_type.parameters[Params::Trate] = true; + Fixture bad_numeric_model(bad_numeric_type); + success *= (bad_numeric_model.hygov.verify() > 0); + + Fixture overflowing_component_base( + makeData(), + {{Params::Trate, std::numeric_limits::max()}}); + success *= (overflowing_component_base.hygov.verify() > 0); + + Fixture overflowing_base_ratio( + makeData(), + {{Params::Trate, std::numeric_limits::min()}}); + success *= (overflowing_base_ratio.hygov.verify() > 0); + + const std::array invalid_system_bases{{ + 0.0, + -1.0, + nan, + infinity, + -infinity, + std::numeric_limits::min(), + }}; + + for (const RealT system_base : invalid_system_bases) + { + Fixture invalid_base(makeData(), {}, system_base); + success *= (invalid_base.hygov.verify() > 0); + } + + success *= unlinkedSignalRejected(); + success *= unlinkedSignalRejected(); + success *= unlinkedSignalRejected(); + + // All five zero time constants use the documented numerical floor and + // still admit a consistent steady-state initialization. + Fixture floors(makeData(), + {{Params::Tr, 0.0}, + {Params::Tf, 0.0}, + {Params::Tg, 0.0}, + {Params::Tw, 0.0}, + {Params::Tnp, 0.0}}); + success *= floors.initialize(0.4); + success *= (floors.evaluate() == 0); + success *= allResidualsZero(floors.hygov); + + return success.report(__func__); + } + + /// A nonidentity power-base initialization with every port attached. + /// The machine-provided pmech value must remain unchanged while HYGOV + /// initializes and publishes its resolved load reference. + TestOutcome initializationAndSignals() + { + TestStatus success = true; + + Fixture fixture(makeData(), {{Params::Trate, 50.0}}); + fixture.attachAllInputs(); + fixture.input(External::PAUX) = 0.02; + fixture.input(External::PREF) = 99.0; // stale value the publication must replace + success *= fixture.initialize(0.4); + success *= (fixture.hygov.tagDifferentiable() == 0); + success *= (fixture.evaluate() == 0); + + const auto* y = fixture.hygov.y().getData(); + success *= scalarMatches(y[static_cast(Internal::XF)], 0.0, "XF at rest"); + success *= scalarMatches(y[static_cast(Internal::C)], + 0.9000000000001573, + "C on component base"); + success *= scalarMatches(y[static_cast(Internal::G)], + 0.9000000000001573, + "G on component base"); + success *= scalarMatches(y[static_cast(Internal::Q)], 0.9, "Q on component base"); + success *= scalarMatches(y[static_cast(Internal::PGV)], + 0.9, + "PGV on component base"); + success *= scalarMatches(y[static_cast(Internal::H)], 1.0, "H at the dam head"); + success *= scalarMatches(fixture.pmech(), 0.4, "preserved pmech value"); + + success *= scalarMatches(fixture.input(External::OMEGA), 0.0, "preserved omega input"); + success *= scalarMatches(fixture.input(External::PREF), 0.0025, "published pref"); + success *= scalarMatches(fixture.input(External::PAUX), 0.02, "preserved paux input"); + + // Verify the six documented outputs through the public monitor controller. + RealT time = 0.0; + Model::VariableMonitorController monitor(time); + monitor.addMonitor(fixture.hygov.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,Hygov_hygov_test_pmech,Hygov_hygov_test_filter," + "Hygov_hygov_test_desiredgate,Hygov_hygov_test_gate," + "Hygov_hygov_test_flow,Hygov_hygov_test_head"); + const auto monitored = Tokenizer(monitor_values, ',')(); + if (monitored.size() == 7) + { + success *= scalarMatches(monitored[1], 0.4, "monitored pmech"); + success *= scalarMatches(monitored[2], 0.0, "monitored filter"); + success *= scalarMatches(monitored[3], 0.9000000000001573, "monitored desiredgate"); + success *= scalarMatches(monitored[4], 0.9000000000001573, "monitored gate"); + success *= scalarMatches(monitored[5], 0.9, "monitored flow"); + success *= scalarMatches(monitored[6], 1.0, "monitored head"); + } + else + { + std::cout << "HYGOV monitor emitted " << monitored.size() + << " values instead of 7\n"; + success = false; + } + + // The five governor states carry derivatives; the rest is algebraic. + for (size_t i = 0; i < static_cast(fixture.hygov.size()); ++i) + { + const bool differential = i <= static_cast(Internal::Q); + if (fixture.hygov.tag()[i] != differential) + { + std::cout << "HYGOV differentiability tag " << i << " mismatch\n"; + success = false; + } + } + + success *= allResidualsZero(fixture.hygov); + + // A system-base reference step lands on the governor error scaled by + // the base ratio. + fixture.input(External::PREF) = 0.1025; // the published 0.0025 plus a 0.1 step + success *= (fixture.evaluate() == 0); + success *= residualsMatch(fixture.hygov, + {{Internal::EF, 0.2}}, + "reference step on the component base"); + + // Unattached ports fall back to the references latched by + // initialize(), so the same steady state holds without a controller. + Fixture fallback(makeData(), {{Params::Trate, 50.0}}); + success *= fallback.initialize(0.4); + success *= (fallback.evaluate() == 0); + success *= allResidualsZero(fallback.hygov); + + return success.report(__func__); + } + + /// Mechanical-power, gate-limit, speed-deviation, and finite-input + /// initialization domains. Response limits and high-power dam head + /// are adjusted when needed; every rejected initialization is atomic. + TestOutcome initializationDomain() + { + TestStatus success = true; + + noteExpectedLogs("Testing HYGOV initialization boundaries. " + "Logged errors, response-limit warnings, and dam-head warnings " + "are expected."); + + success *= initializationRejectedAtomically( + makeResidualData(), + -0.3, + {{External::OMEGA, 0.0}, {External::PREF, 77.0}, {External::PAUX, 0.02}}, + "mechanical power below the gate curve"); + + const auto no_finite_head = withParameters( + makeData(), + {{Params::Pgv0, -1.0}, + {Params::Pgv1, -0.8}, + {Params::Pgv2, -0.6}, + {Params::Pgv3, -0.4}, + {Params::Pgv4, -0.2}, + {Params::Pgv5, 0.0}}); + success *= initializationRejectedAtomically( + no_finite_head, + 0.0, + {{External::OMEGA, 0.0}, {External::PREF, 77.0}, {External::PAUX, 0.02}}, + "no finite effective Hdam"); + + // 4.5 MW on the system base is 2.5 pu on a 1.8 MW turbine base. + Fixture effective_fixture( + makeData(), + {{Params::Trate, 1.8}, + {Params::At, 1.25}, + {Params::Qnl, 0.07}, + {Params::Gmax, 0.5}}); + effective_fixture.attachAllInputs(); + success *= effective_fixture.initialize(0.045); + success *= stateMatches( + effective_fixture.hygov, + {{Internal::C, 1.0}, + {Internal::G, 1.0}, + {Internal::Q, 1.2812656647316965}, + {Internal::PGV, 0.9971118867476669}, + {Internal::H, 1.6511654364800423}}, + "effective dam head"); + success *= scalarMatches(effective_fixture.pmech(), 0.045, "preserved pmech value"); + success *= scalarMatches(effective_fixture.input(External::PREF), + 0.0009, + "published pref"); + success *= (effective_fixture.evaluate() == 0); + success *= allResidualsZero(effective_fixture.hygov); + + struct ResponseLimitCase + { + const char* label; + Params limit_parameter; + RealT limit; + RealT rate; + }; + + const std::array response_limit_cases{{ + {"expanded upper response limit", Params::Gmax, 0.5, 0.1}, + {"expanded lower response limit", Params::Gmin, 0.7, -0.1}, + }}; + + for (const auto& test_case : response_limit_cases) + { + Fixture fixture(makeResidualData(), + {{test_case.limit_parameter, test_case.limit}}); + success *= fixture.initialize(0.4); + const RealT gate = static_cast( + fixture.hygov.y().getData()[static_cast(Internal::C)]); + const bool gate_is_outside = test_case.rate > 0.0 + ? gate > test_case.limit + : gate < test_case.limit; + if (!gate_is_outside) + { + std::cout << test_case.label << " did not initialize outside the configured limit\n"; + success = false; + } + success *= stateMatches(fixture.hygov, + {{Internal::G, gate}, {Internal::H, 1.2}}, + test_case.label); + success *= (fixture.evaluate() == 0); + success *= allResidualsZero(fixture.hygov); + + // The effective response bound admits an outward rate between the + // configured limit and initialized gate. + setState(fixture.hygov, + {{Internal::C, 0.5 * (test_case.limit + gate)}, + {Internal::RC, test_case.rate}}); + setDerivative(fixture.hygov, {{Internal::C, 0.0}}); + success *= (fixture.evaluate() == 0); + const RealT response_rate = static_cast( + fixture.hygov.getResidual().getData()[static_cast(Internal::C)]); + const bool rate_is_admitted = test_case.rate > 0.0 + ? response_rate > 0.9 * test_case.rate + : response_rate < 0.9 * test_case.rate; + if (!rate_is_admitted) + { + std::cout << test_case.label << " did not admit the outward desired-gate rate\n"; + success = false; + } + } + + // A failed retry preserves the effective head and response bounds from + // the prior success. + const auto effective_y = copyVector(effective_fixture.hygov.y()); + const auto effective_yp = copyVector(effective_fixture.hygov.yp()); + effective_fixture.input(External::OMEGA) = 0.03; + success *= (effective_fixture.hygov.initialize() != 0); + success *= vectorUnchanged(effective_fixture.hygov.y(), effective_y, "state"); + success *= vectorUnchanged(effective_fixture.hygov.yp(), effective_yp, "derivative"); + success *= scalarMatches(effective_fixture.input(External::PREF), + 0.0009, + "preserved pref"); + effective_fixture.input(External::OMEGA) = 0.0; + success *= (effective_fixture.evaluate() == 0); + success *= allResidualsZero(effective_fixture.hygov); + + setState(effective_fixture.hygov, + {{Internal::C, 0.75}, {Internal::RC, 0.2}}); + setDerivative(effective_fixture.hygov, {{Internal::C, 0.0}}); + success *= (effective_fixture.evaluate() == 0); + const RealT preserved_rate = static_cast( + effective_fixture.hygov.getResidual().getData()[static_cast(Internal::C)]); + if (!(preserved_rate > 0.19)) + { + std::cout << "failed initialization did not preserve effective Gmax\n"; + success = false; + } + + // A later feasible initialization starts again from configured limits + // and Hdam. + effective_fixture.setPmech(0.009); + success *= (effective_fixture.hygov.initialize() == 0); + success *= stateMatches(effective_fixture.hygov, + {{Internal::H, 1.0}}, + "configured dam head after reinitialization"); + success *= (effective_fixture.evaluate() == 0); + success *= allResidualsZero(effective_fixture.hygov); + + setState(effective_fixture.hygov, + {{Internal::C, 0.75}, {Internal::RC, 0.2}}); + setDerivative(effective_fixture.hygov, {{Internal::C, 0.0}}); + success *= (effective_fixture.evaluate() == 0); + success *= scalarMatches( + static_cast(effective_fixture.hygov.getResidual().getData()[static_cast(Internal::C)]), + 0.0, + "configured Gmax after reinitialization"); + + // Initialization supports only a zero speed deviation; a moving + // machine would need a multi-root gate search. + success *= initializationRejectedAtomically(makeResidualData(), + 0.4, + {{External::OMEGA, 0.03}, + {External::PREF, 77.0}, + {External::PAUX, 0.02}}, + "nonzero initial speed deviation"); + + // An invalid configuration is rejected before any state is written. + Fixture invalid_fixture(makeResidualData(), {{Params::Rtemp, 0.0}}); + invalid_fixture.attachAllInputs(); + success *= (invalid_fixture.hygov.allocate() == 0); + poisonState(invalid_fixture, 0.4); + const auto invalid_y = copyVector(invalid_fixture.hygov.y()); + const auto invalid_yp = copyVector(invalid_fixture.hygov.yp()); + if (invalid_fixture.hygov.initialize() == 0) + { + std::cout << "Expected initialization rejection: invalid configuration\n"; + success = false; + } + success *= vectorUnchanged(invalid_fixture.hygov.y(), invalid_y, "state"); + success *= vectorUnchanged(invalid_fixture.hygov.yp(), invalid_yp, "derivative"); + + // Zero mechanical power lands on an in-range root and initializes at rest. + Fixture zero_power_fixture(makeData()); + success *= zero_power_fixture.initialize(0.0); + success *= stateMatches( + zero_power_fixture.hygov, + {{Internal::C, 0.09999999999984271}, {Internal::G, 0.09999999999984271}}, + "zero mechanical power"); + success *= (zero_power_fixture.evaluate() == 0); + success *= allResidualsZero(zero_power_fixture.hygov); + + // The smooth identity curve leaves a ln(2)/MU knee at each end, so + // makeData()'s achievable component-base power range is + // [knee - 0.1, 0.9 - knee]. + const RealT knee = std::log(static_cast(2.0)) / Math::MU; + const RealT p_max = static_cast(0.9) - knee; + const RealT p_min = knee - static_cast(0.1); + + Fixture lower_edge(makeData()); + success *= lower_edge.initialize(p_min - 0.5 * kTol); + success *= stateMatches(lower_edge.hygov, + {{Internal::C, 0.0}, {Internal::G, 0.0}}, + "half the tolerance below the achievable minimum"); + success *= scalarMatches(lower_edge.pmech(), + p_min - 0.5 * kTol, + "clipped pmech value"); + success *= (lower_edge.evaluate() == 0); + success *= allResidualsZero(lower_edge.hygov); + + Fixture effective_edge(makeData()); + success *= effective_edge.initialize(p_max + 0.5 * kTol); + success *= stateMatches(effective_edge.hygov, + {{Internal::C, 1.0}, {Internal::G, 1.0}}, + "half the tolerance beyond the achievable maximum"); + const RealT effective_edge_head = static_cast( + effective_edge.hygov.y().getData()[static_cast(Internal::H)]); + if (!(effective_edge_head > 1.0)) + { + std::cout << "effective head was not raised above configured Hdam\n"; + success = false; + } + success *= (effective_edge.evaluate() == 0); + success *= allResidualsZero(effective_edge.hygov); + + success *= initializationRejectedAtomically( + makeData(), + p_min - 2.0 * kTol, + {{External::OMEGA, 0.0}, {External::PREF, 77.0}, {External::PAUX, 0.02}}, + "twice the tolerance below the achievable minimum"); + + const RealT nan = std::numeric_limits::quiet_NaN(); + const RealT infinity = std::numeric_limits::infinity(); + + // A non-finite input is rejected atomically, NaN included: the + // exact-preservation check states what a tolerance comparison of a + // NaN input never could. + const std::array nonfinite_inputs{{nan, infinity, -infinity}}; + + for (const RealT value : nonfinite_inputs) + { + success *= initializationRejectedAtomically(makeData(), + 0.4, + {{External::OMEGA, value}, + {External::PREF, 77.0}, + {External::PAUX, 0.02}}, + "non-finite speed input"); + success *= initializationRejectedAtomically(makeData(), + 0.4, + {{External::OMEGA, 0.0}, + {External::PREF, 77.0}, + {External::PAUX, value}}, + "non-finite auxiliary-power input"); + + // A non-finite seed lands in the aliased pmech state itself, so the + // poisoned-state comparison cannot express its preservation. The + // inputs still must survive untouched. + Fixture pmech_fixture(makeData()); + pmech_fixture.attachAllInputs(); + pmech_fixture.input(External::PREF) = 77.0; + pmech_fixture.input(External::PAUX) = 0.02; + success *= pmech_fixture.prepare(value); + success *= (pmech_fixture.hygov.initialize() != 0); + success *= scalarPreserved( + static_cast(pmech_fixture.input(External::PREF)), + 77.0, + "external input", + static_cast(External::PREF)); + success *= scalarPreserved(static_cast(pmech_fixture.input(External::PAUX)), + 0.02, + "external input", + static_cast(External::PAUX)); + } + + return success.report(__func__); + } + + /// Initialization solves the smooth gate curve the residual evaluates, + /// so every steady residual rests at machine rounding even where the + /// smoothing bends the curve away from its piecewise-linear points. + TestOutcome initializationExactness() + { + TestStatus success = true; + + // Values landing mid-segment and within the smoothing knee of every + // interior curve breakpoint, where a piecewise-linear inversion + // misses the implemented curve by up to O(1e-3). The gate literal + // proves where each value lands. + struct ExactnessCase + { + const char* label; + RealT pmech; + RealT gate; + }; + + const std::array exactness_cases{{ + {"gate inside the Gv1 knee", 0.0556, 0.1982318164100278}, + {"gate inside the Gv2 knee", 0.2509, 0.4003865335541374}, + {"gate mid-segment", 0.4, 0.5719050089028755}, + {"gate inside the Gv3 knee", 0.4244, 0.6007061471851347}, + {"gate inside the Gv4 knee", 0.5617, 0.8006094230811988}, + }}; + + for (const auto& seed : exactness_cases) + { + Fixture fixture(makeResidualData()); + success *= fixture.initialize(seed.pmech); + success *= stateMatches(fixture.hygov, {{Internal::G, seed.gate}}, seed.label); + success *= (fixture.evaluate() == 0); + success *= allResidualsZero(fixture.hygov); + } + + return success.report(__func__); + } + + /// A fixed numerical answer key for all 12 HYGOV residual rows. The + /// expected values are literals, not a second implementation of HYGOV. + TestOutcome residualEquations() + { + TestStatus success = true; + + Fixture fixture(makeResidualData()); + fixture.attachAllInputs(); + success *= fixture.initialize(0.4); + setAnswerKeyInputs(fixture); + setAnswerKeyState(fixture.hygov); + success *= (fixture.evaluate() == 0); + + const std::array(Internal::MAXIMUM)> expected{{ + {Internal::XN, -0.07785714285714286}, + {Internal::XF, -0.7300000000000001}, + {Internal::C, 0.06}, + {Internal::G, 0.1233333333333334}, + {Internal::Q, 0.011538461538461414}, + {Internal::OMEGADB, 0.0033514666467982894}, + {Internal::EF, 0.5863}, + {Internal::FC, -0.7405000000000002}, + {Internal::RC, 0.029996890386450745}, + {Internal::PGV, -0.04600000003160343}, + {Internal::H, -0.033299999999999885}, + {Internal::PMECH, -0.012679999999999934}, + }}; + + success *= (static_cast(fixture.hygov.getResidual().getSize()) == expected.size()); + success *= residualsMatch(fixture.hygov, expected, "answer key"); + + return success.report(__func__); + } + + /// Speed deadband, desired-gate velocity limiting, and gate-position + /// anti-windup. + TestOutcome governorControl() + { + TestStatus success = true; + const auto data = makeResidualData(); + + // Exercise both sides and the interior of the type-1 +/-0.01 deadband. + const std::array deadband_cases{{ + {"speed deadband below the band", + {{External::OMEGA, -0.05}}, + {{Internal::OMEGADB, 0.0}}, + {}, + {{Internal::OMEGADB, -0.049996641662021946}}}, + {"speed deadband inside the band", + {{External::OMEGA, 0.004}}, + {{Internal::OMEGADB, 0.0}}, + {}, + {{Internal::OMEGADB, 0.0009004582873718001}}}, + {"speed deadband above the band", + {{External::OMEGA, 0.05}}, + {{Internal::OMEGADB, 0.0}}, + {}, + {{Internal::OMEGADB, 0.049996641662021946}}}, + }}; + success *= runResidualCases(data, 0.4, deadband_cases); + + const std::array gate_velocity_cases{{ + {"gate velocity below the rate limit", + {}, + {{Internal::FC, -0.6}, {Internal::RC, 0.0}}, + {}, + {{Internal::RC, -0.15}}}, + {"gate velocity inside the rate limit", + {}, + {{Internal::FC, 0.05}, {Internal::RC, 0.0}}, + {}, + {{Internal::RC, 0.04999999999984272}}}, + {"gate velocity above the rate limit", + {}, + {{Internal::FC, 0.6}, {Internal::RC, 0.0}}, + {}, + {{Internal::RC, 0.15000000000000002}}}, + }}; + success *= runResidualCases(data, 0.4, gate_velocity_cases); + + const std::array gate_antiwindup_cases{{ + {"Gmax blocks an outward desired-gate rate", + {}, + {{Internal::C, 1.2}, {Internal::RC, 0.2}}, + {{Internal::C, 0.0}}, + {{Internal::C, 0.0}}}, + {"Gmin blocks an outward desired-gate rate", + {}, + {{Internal::C, -0.2}, {Internal::RC, -0.2}}, + {{Internal::C, 0.0}}, + {{Internal::C, 0.0}}}, + {"Gmax admits a restoring desired-gate rate", + {}, + {{Internal::C, 1.2}, {Internal::RC, -0.2}}, + {{Internal::C, 0.0}}, + {{Internal::C, -0.2}}}, + {"Gmin admits a restoring desired-gate rate", + {}, + {{Internal::C, -0.2}, {Internal::RC, 0.2}}, + {{Internal::C, 0.0}}, + {{Internal::C, 0.2}}}, + }}; + success *= runResidualCases(data, 0.4, gate_antiwindup_cases); + + // At alpha = 1, a blocked desired-gate row has derivative coefficient + // -1 and no RC dependence, independently of either Jacobian backend. + { + using DepVar = DependencyTracking::Variable; + + Fixture blocked(data); + blocked.attachAllInputs(); + success *= blocked.initialize(0.4); + setState(blocked.hygov, {{Internal::C, 1.2}, {Internal::RC, 0.2}}); + setDerivative(blocked.hygov, {{Internal::C, 0.0}}); + numberVariables(blocked); + success *= (blocked.evaluate() == 0); + + const auto& dependencies = + blocked.hygov.getResidual().getData()[static_cast(Internal::C)].getDependencies(); + const DepVar::DependencyMap expected{{ + {static_cast(Internal::C), -1.0}, + {static_cast(Internal::RC), 0.0}, + }}; + success *= isEqual(dependencies, expected, kTol); + } + + return success.report(__func__); + } + + /// Gate-power, water-column, damping, and curve-inversion behavior, + /// including a flat segment. + TestOutcome turbineDynamics() + { + TestStatus success = true; + const auto data = makeResidualData(); + + const std::array gate_power_cases{{ + {"gate-power curve segment 1", + {}, + {{Internal::G, 0.1}, {Internal::PGV, 0.0}}, + {}, + {{Internal::PGV, 0.07500000000021236}}}, + {"gate-power curve segment 2", + {}, + {{Internal::G, 0.3}, {Internal::PGV, 0.0}}, + {}, + {{Internal::PGV, 0.28500000000007075}}}, + {"gate-power curve segment 3", + {}, + {{Internal::G, 0.5}, {Internal::PGV, 0.0}}, + {}, + {{Internal::PGV, 0.5399999999999371}}}, + {"gate-power curve segment 4", + {}, + {{Internal::G, 0.7}, {Internal::PGV, 0.0}}, + {}, + {{Internal::PGV, 0.7549999999999292}}}, + {"gate-power curve segment 5", + {}, + {{Internal::G, 0.9}, {Internal::PGV, 0.0}}, + {}, + {{Internal::PGV, 0.9249999999998506}}}, + }}; + success *= runResidualCases(data, 0.4, gate_power_cases); + + // A head away from the dam head drives the flow and head rows, and + // turbine damping scales with speed deviation and gate. + const std::array turbine_cases{{ + {"water column", + {}, + {{Internal::Q, 0.61}, {Internal::H, 0.9}, {Internal::PGV, 0.55}}, + {{Internal::Q, 0.05}}, + {{Internal::Q, 0.18076923076923068}, {Internal::H, -0.09984999999999994}}}, + {"turbine damping", + {{External::OMEGA, 0.05}}, + {{Internal::G, 0.6}, + {Internal::Q, 0.7}, + {Internal::H, 1.1}, + {Internal::PMECH, 0.5}}, + {}, + {{Internal::PMECH, -0.2677999999999999}}}, + }}; + success *= runResidualCases(data, 0.4, turbine_cases); + + Fixture curve_fixture(data); + curve_fixture.attachAllInputs(); + success *= curve_fixture.initialize(0.33761676); + success *= stateMatches( + curve_fixture.hygov, + {{Internal::C, 0.5000001394783365}, {Internal::G, 0.5000001394783365}}, + "nonidentity curve inversion"); + success *= scalarMatches(curve_fixture.input(External::PREF), + 0.015000004184348527, + "nonidentity-curve published pref"); + success *= scalarMatches(curve_fixture.pmech(), 0.33761676, "preserved pmech value"); + success *= (curve_fixture.evaluate() == 0); + success *= allResidualsZero(curve_fixture.hygov); + + // A flat source-curve segment must initialize to a gate on that segment. + // makeData() uses equal power bases, At = Hdam = 1, and Qnl = 0.1, + // so a 0.5 plateau maps to pmech = 0.4 without encoding Math::MU. + const RealT flat_gate_minimum = static_cast(0.4); + const RealT flat_gate_maximum = static_cast(0.6); + const RealT plateau_power = static_cast(0.5); + const RealT plateau_pmech = static_cast(0.4); + Fixture flat_fixture(makeData(), + {{Params::Pgv2, plateau_power}, + {Params::Pgv3, plateau_power}}); + success *= flat_fixture.initialize(plateau_pmech); + const RealT flat_gate = + flat_fixture.hygov.y().getData()[static_cast(Internal::G)]; + if (flat_gate < flat_gate_minimum || flat_gate > flat_gate_maximum) + { + std::cout << "flat-segment plateau gate " + << std::setprecision(std::numeric_limits::max_digits10) + << flat_gate << " is outside [" << flat_gate_minimum + << ", " << flat_gate_maximum << "]\n"; + success = false; + } + success *= (flat_fixture.evaluate() == 0); + success *= allResidualsZero(flat_fixture.hygov); + + return success.report(__func__); + } + +#ifdef GRIDKIT_ENABLE_ENZYME + /// Every Enzyme CSR row must match dependency tracking at gates inside + /// each curve segment and at each breakpoint, and both paths must + /// carry the PGV row's gate dependence. + TestOutcome jacobian() + { + TestStatus success = true; + + const auto data = makeResidualData(); + const std::array gate_points{{0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9}}; + + for (const RealT gate : gate_points) + { + const auto dependency_jacobian = dependencyTrackingJacobian(data, gate, success); + const auto enzyme_jacobian = enzymeJacobian(data, gate, 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 << "HYGOV Jacobian row " << row << " at gate " << gate + << " mismatch between dependency tracking and Enzyme\n"; + success = false; + } + } + + // Guard the required PGV/G dependency even if both paths agree. + success *= jacobianContains( + dependency_jacobian, Internal::PGV, Internal::G, "dependency-tracking"); + success *= jacobianContains(enzyme_jacobian, Internal::PGV, Internal::G, "Enzyme"); + } + + return success.report(__func__); + } +#endif + + private: + using Params = PhasorDynamics::Governor::HygovParameters; + using Internal = PhasorDynamics::Governor::HygovInternalVariables; + using External = PhasorDynamics::Governor::HygovExternalVariables; + using Mon = PhasorDynamics::Governor::HygovMonitorableVariables; + using Data = PhasorDynamics::Governor::HygovData; + using HygovT = PhasorDynamics::Governor::Hygov; + + using InternalRow = std::pair; + using InternalRows = std::vector; + using ExternalRow = std::pair; + using ExternalRows = std::vector; + + /// Failure-report names for the internal rows, ordered as `Internal`. + static constexpr std::array(Internal::MAXIMUM)> kRowNames{ + {"XN", "XF", "C", "G", "Q", "OMEGADB", "EF", "FC", "RC", "PGV", "H", "PMECH"}}; + + struct ResidualCase + { + const char* label; + ExternalRows inputs; + InternalRows state; + InternalRows derivative; + InternalRows expected; + }; + + static Data withParameters(Data data, + std::initializer_list> overrides) + { + for (const auto& [parameter, value] : overrides) + { + data.parameters[parameter] = value; + } + return data; + } + + /// Owns the HYGOV model, the assigned mechanical-power node, and the + /// attached input nodes. Signal storage is declared before the model so + /// every referenced node outlives HYGOV. 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 pmech_node_; + + public: + explicit Fixture(const Data& data, + std::initializer_list> overrides = {}, + RealT system_va_base = 100.0e6) + : hygov(withParameters(data, overrides)) + { + hygov.setSystemBase(60.0, system_va_base); + hygov.getSignals().template assignSignalNode(&pmech_node_); + } + + Fixture(const Fixture&) = delete; + Fixture& operator=(const Fixture&) = delete; + + void attachAllInputs(RealT initial_value = 0.0) + { + const IdxT external_index_base = hygov.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 = hygov.getSignals(); + signals.template attachSignalNode( + &input_nodes_[static_cast(External::OMEGA)]); + signals.template attachSignalNode( + &input_nodes_[static_cast(External::PREF)]); + signals.template attachSignalNode( + &input_nodes_[static_cast(External::PAUX)]); + } + + /// Set the assigned mechanical-power node on the system base. + void setPmech(RealT pmech) + { + pmech_node_.init(static_cast(pmech)); + } + + /// Everything HYGOV initialization requires: allocation, + /// verification, and a machine-provided mechanical-power value. + bool prepare(RealT pmech) + { + const bool success = (hygov.allocate() == 0) && (hygov.verify() == 0); + if (!success) + { + std::cout << "HYGOV fixture preparation failed\n"; + return false; + } + + setPmech(pmech); + return true; + } + + bool initialize(RealT pmech) + { + if (!prepare(pmech)) + { + return false; + } + if (hygov.initialize() != 0) + { + std::cout << "HYGOV initialization failed\n"; + return false; + } + return true; + } + + int evaluate() + { + return hygov.evaluateResidual(); + } + + T pmech() const + { + return pmech_node_.read(); + } + + T& input(External port) + { + return input_values_[static_cast(port)]; + } + + IdxT inputIndex(External port) const + { + return input_indices_[static_cast(port)]; + } + + PhasorDynamics::Governor::Hygov hygov; + }; + + Data makeMinimalData() const + { + Data data; + data.device_class = "Hygov"; + data.disambiguation_string = "hygov_test"; + data.monitored_variables.insert(Mon::pmech); + data.monitored_variables.insert(Mon::filter); + data.monitored_variables.insert(Mon::desiredgate); + data.monitored_variables.insert(Mon::gate); + data.monitored_variables.insert(Mon::flow); + data.monitored_variables.insert(Mon::head); + return data; + } + + Data makeExplicitDefaultData() const + { + // These are the documented defaults. The all-zero source curve + // selects the identity curve, spelled out here point by point. + return withParameters(makeMinimalData(), + {{Params::Rperm, 0.04}, + {Params::Rtemp, 0.3}, + {Params::Tr, 5.0}, + {Params::Tf, 0.05}, + {Params::Tg, 0.5}, + {Params::Velm, 0.2}, + {Params::Gmax, 1.0}, + {Params::Gmin, 0.0}, + {Params::Tw, 1.0}, + {Params::At, 1.2}, + {Params::Dturb, 0.5}, + {Params::Qnl, 0.05}, + {Params::Tn, 0.0}, + {Params::Tnp, 0.0}, + {Params::db1, 0.0}, + {Params::db2, 0.0}, + {Params::Hdam, 1.0}, + {Params::Gv0, 0.0}, + {Params::Gv1, 0.2}, + {Params::Gv2, 0.4}, + {Params::Gv3, 0.6}, + {Params::Gv4, 0.8}, + {Params::Gv5, 1.0}, + {Params::Pgv0, 0.0}, + {Params::Pgv1, 0.2}, + {Params::Pgv2, 0.4}, + {Params::Pgv3, 0.6}, + {Params::Pgv4, 0.8}, + {Params::Pgv5, 1.0}}); + } + + Data makeData() const + { + // The documented typical values with the floored time constants + // raised above the floor, so routine fixtures log no warnings. + return withParameters(makeMinimalData(), + {{Params::Trate, 100.0}, + {Params::Rperm, 0.05}, + {Params::Rtemp, 0.4}, + {Params::Tr, 5.0}, + {Params::Tf, 0.2}, + {Params::Tg, 0.5}, + {Params::Velm, 0.5}, + {Params::Gmax, 1.0}, + {Params::Gmin, 0.0}, + {Params::Tw, 1.0}, + {Params::At, 1.0}, + {Params::Dturb, 0.0}, + {Params::Qnl, 0.1}, + {Params::Tn, 0.0}, + {Params::Tnp, 1.0}, + {Params::db1, 0.0}, + {Params::db2, 0.0}, + {Params::Hdam, 1.0}, + {Params::Gv0, 0.0}, + {Params::Gv1, 0.2}, + {Params::Gv2, 0.4}, + {Params::Gv3, 0.6}, + {Params::Gv4, 0.8}, + {Params::Gv5, 1.0}, + {Params::Pgv0, 0.0}, + {Params::Pgv1, 0.2}, + {Params::Pgv2, 0.4}, + {Params::Pgv3, 0.6}, + {Params::Pgv4, 0.8}, + {Params::Pgv5, 1.0}}); + } + + Data makeResidualData() const + { + return withParameters(makeData(), + {{Params::Trate, 50.0}, + {Params::Rperm, 0.06}, + {Params::Rtemp, 0.4}, + {Params::Tr, 4.0}, + {Params::Tf, 0.2}, + {Params::Tg, 0.6}, + {Params::Velm, 0.15}, + {Params::Gmax, 0.95}, + {Params::Gmin, 0.05}, + {Params::Tw, 1.3}, + {Params::At, 1.1}, + {Params::Dturb, 0.6}, + {Params::Qnl, 0.08}, + {Params::Tn, 0.7}, + {Params::Tnp, 1.4}, + {Params::db1, 0.01}, + {Params::Hdam, 1.2}, + {Params::Pgv1, 0.15}, + {Params::Pgv2, 0.42}, + {Params::Pgv3, 0.66}, + {Params::Pgv4, 0.85}}); + } + + template + void setAnswerKeyInputs(Fixture& fixture) const + { + fixture.input(External::OMEGA) = static_cast(0.02); + fixture.input(External::PREF) = static_cast(0.31); + fixture.input(External::PAUX) = static_cast(0.07); + } + + /// The rich state shared by the residual answer key and the Jacobian + /// comparison. Every row is distinct so a swapped index cannot pass. + template + void setAnswerKeyState(PhasorDynamics::Governor::Hygov& hygov) const + { + setState(hygov, + {{Internal::XN, 0.11}, + {Internal::XF, 0.23}, + {Internal::C, 0.52}, + {Internal::G, 0.47}, + {Internal::Q, 0.61}, + {Internal::OMEGADB, 0.015}, + {Internal::EF, 0.08}, + {Internal::FC, 0.12}, + {Internal::RC, 0.09}, + {Internal::PGV, 0.55}, + {Internal::H, 1.12}, + {Internal::PMECH, 0.33}}); + setDerivative(hygov, + {{Internal::XN, 0.01}, + {Internal::XF, -0.02}, + {Internal::C, 0.03}, + {Internal::G, -0.04}, + {Internal::Q, 0.05}}); + } + + /// Omitting every optional 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(), {}, 200.0e6); + Fixture explicit_defaults(makeExplicitDefaultData(), {}, 200.0e6); + implicit_defaults.attachAllInputs(); + explicit_defaults.attachAllInputs(); + + bool success = implicit_defaults.initialize(0.3) + && explicit_defaults.initialize(0.3); + if (!success) + { + std::cout << "HYGOV 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.hygov.y(), + copyVector(explicit_defaults.hygov.y()), + "documented-default state")) + { + success = false; + } + if (!vectorUnchanged(implicit_defaults.hygov.yp(), + copyVector(explicit_defaults.hygov.yp()), + "documented-default derivative")) + { + success = false; + } + if (!vectorUnchanged(implicit_defaults.hygov.getResidual(), + copyVector(explicit_defaults.hygov.getResidual()), + "documented-default residual")) + { + success = false; + } + + setAnswerKeyInputs(implicit_defaults); + setAnswerKeyInputs(explicit_defaults); + setAnswerKeyState(implicit_defaults.hygov); + setAnswerKeyState(explicit_defaults.hygov); + if (implicit_defaults.evaluate() != 0) + { + success = false; + } + if (explicit_defaults.evaluate() != 0) + { + success = false; + } + if (!vectorUnchanged(implicit_defaults.hygov.getResidual(), + copyVector(explicit_defaults.hygov.getResidual()), + "documented-default dynamic residual")) + { + success = false; + } + return success; + } + + template + bool unlinkedSignalRejected() const + { + PhasorDynamics::SignalNode unlinked_node; + Fixture fixture(makeData()); + fixture.hygov.getSignals().template attachSignalNode(&unlinked_node); + return fixture.hygov.verify() > 0; + } + + template + std::vector copyVector(const VectorT& vector) const + { + const auto* values = vector.getData(); + return std::vector(values, + values + static_cast(vector.getSize())); + } + + 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 << "HYGOV " << what << " row " << row + << " changed mismatch: " << actual << " != " << expected << "\n"; + } + return ret; + } + + /// Fill the state and derivative with a recognizable ramp, then restore + /// the aliased pmech entry, so any write by a rejected initialization + /// is visible. + void poisonState(Fixture& fixture, RealT pmech) const + { + auto* y = fixture.hygov.y().getData(); + auto* yp = fixture.hygov.yp().getData(); + for (size_t i = 0; i < static_cast(fixture.hygov.y().getSize()); ++i) + { + y[i] = 0.125 + 0.01 * static_cast(i); + yp[i] = -0.25 - 0.01 * static_cast(i); + } + fixture.setPmech(pmech); + fixture.hygov.y().setDataUpdated(); + fixture.hygov.yp().setDataUpdated(); + } + + /// Initialization must fail and leave the poisoned state, the seeded + /// pmech value, and every supplied input untouched. + bool initializationRejectedAtomically(const Data& data, + RealT pmech, + const ExternalRows& inputs, + const char* label) const + { + Fixture fixture(data); + fixture.attachAllInputs(); + for (const auto& [port, value] : inputs) + { + fixture.input(port) = static_cast(value); + } + if (!fixture.prepare(pmech)) + { + return false; + } + + poisonState(fixture, pmech); + const auto y_before = copyVector(fixture.hygov.y()); + const auto yp_before = copyVector(fixture.hygov.yp()); + + bool success = true; + if (fixture.hygov.initialize() == 0) + { + std::cout << "Expected initialization rejection: " << label << "\n"; + success = false; + } + + if (!scalarMatches(fixture.pmech(), pmech, "rejected pmech 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.hygov.y(), y_before, "state")) + { + success = false; + } + if (!vectorUnchanged(fixture.hygov.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::Governor::Hygov& hygov, + const InternalRows& rows) const + { + auto* y = hygov.y().getData(); + for (const auto& [variable, value] : rows) + { + y[static_cast(variable)] = static_cast(value); + } + hygov.y().setDataUpdated(); + } + + template + void setDerivative(PhasorDynamics::Governor::Hygov& hygov, + const InternalRows& rows) const + { + auto* yp = hygov.yp().getData(); + for (const auto& [variable, value] : rows) + { + yp[static_cast(variable)] = static_cast(value); + } + hygov.yp().setDataUpdated(); + } + + /// Evaluate each scenario on a fresh fixture to prevent state leakage. + template + bool runResidualCases(const Data& data, + RealT pmech, + const std::array& cases) const + { + bool success = true; + for (const auto& test_case : cases) + { + Fixture fixture(data); + fixture.attachAllInputs(); + if (!fixture.initialize(pmech)) + { + success = false; + continue; + } + for (const auto& [port, value] : test_case.inputs) + { + fixture.input(port) = static_cast(value); + } + setState(fixture.hygov, test_case.state); + setDerivative(fixture.hygov, test_case.derivative); + if (fixture.evaluate() != 0) + { + success = false; + } + if (!residualsMatch(fixture.hygov, test_case.expected, test_case.label)) + { + success = false; + } + } + return success; + } + + /// Compare one named row and report mismatches consistently. + bool rowMatches(RealT actual, + RealT expected, + const char* what, + size_t row, + const char* context) const + { + if (isEqual(actual, expected, kTol)) + { + return true; + } + std::cout << "HYGOV " << what << " row "; + if (row < kRowNames.size()) + { + std::cout << kRowNames[row]; + } + else + { + std::cout << row; + } + std::cout << ' ' << context << " mismatch: " + << std::setprecision(std::numeric_limits::max_digits10) + << actual << " != " << expected << '\n'; + return false; + } + + 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 HygovT& hygov, + const InternalRows& rows, + const char* context = "") const + { + return rowsMatch(hygov.getResidual(), rows, "residual", context); + } + + template + bool residualsMatch(const HygovT& hygov, + const std::array& rows, + const char* context = "") const + { + return rowsMatch(hygov.getResidual(), rows, "residual", context); + } + + bool stateMatches(const HygovT& hygov, + const InternalRows& rows, + const char* context = "") const + { + return rowsMatch(hygov.y(), rows, "state", context); + } + + /// The model sits at a steady state: every residual and every + /// derivative is zero. + bool allResidualsZero(const HygovT& hygov) const + { + bool success = true; + const auto* f = hygov.getResidual().getData(); + const auto* yp = hygov.yp().getData(); + for (size_t row = 0; row < static_cast(hygov.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(std::numeric_limits::max_digits10) + << 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.hygov.y().getData(); + auto* yp = fixture.hygov.yp().getData(); + + const auto model_size = static_cast(fixture.hygov.size()); + for (size_t i = 0; i < model_size; ++i) + { + y[i].setVariableNumber(i); + yp[i].setVariableNumber(i); + } + for (External port : {External::OMEGA, External::PREF, External::PAUX}) + { + fixture.input(port).setVariableNumber(fixture.inputIndex(port)); + } + + fixture.hygov.y().setDataUpdated(); + fixture.hygov.yp().setDataUpdated(); + } + +#ifdef GRIDKIT_ENABLE_ENZYME + template + bool jacobianContains(const JacobianRowsT& rows, + Internal row_variable, + Internal column_variable, + const char* what) const + { + const auto row = static_cast(row_variable); + const auto column = static_cast(column_variable); + if (row < rows.size() && rows[row].count(column) == 1) + { + return true; + } + std::cout << "HYGOV " << what << " Jacobian row " << row + << " is missing column " << column << "\n"; + return false; + } + + std::vector dependencyTrackingJacobian( + const Data& data, + RealT gate, + TestStatus& success) const + { + using DepVar = DependencyTracking::Variable; + + Fixture fixture(data); + fixture.attachAllInputs(); + success *= fixture.initialize(0.4); + setAnswerKeyInputs(fixture); + setAnswerKeyState(fixture.hygov); + setState(fixture.hygov, {{Internal::G, gate}}); + numberVariables(fixture); + success *= (fixture.evaluate() == 0); + + const auto model_size = static_cast(fixture.hygov.size()); + std::vector rows(model_size); + const auto* f = fixture.hygov.getResidual().getData(); + for (size_t i = 0; i < model_size; ++i) + { + rows[i] = f[i].getDependencies(); + } + return rows; + } + + std::vector enzymeJacobian( + const Data& data, + RealT gate, + TestStatus& success) const + { + Fixture fixture(data); + fixture.attachAllInputs(); + success *= fixture.initialize(0.4); + setAnswerKeyInputs(fixture); + setAnswerKeyState(fixture.hygov); + setState(fixture.hygov, {{Internal::G, gate}}); + fixture.hygov.updateTime(0.0, 1.0); + success *= (fixture.evaluate() == 0); + success *= (fixture.hygov.evaluateJacobian() == 0); + success *= (fixture.hygov.constructCsr() == 0); + return MapFromCsr(fixture.hygov.getCsrJacobian()); + } +#endif + }; + } // namespace Testing +} // namespace GridKit diff --git a/tests/UnitTests/PhasorDynamics/SystemSingleComponentTests.hpp b/tests/UnitTests/PhasorDynamics/SystemSingleComponentTests.hpp index dae9bbb57..de98a440f 100644 --- a/tests/UnitTests/PhasorDynamics/SystemSingleComponentTests.hpp +++ b/tests/UnitTests/PhasorDynamics/SystemSingleComponentTests.hpp @@ -321,6 +321,58 @@ namespace GridKit return success.report(__func__); } + /// HYGOV through the production path: model data to system + /// construction to signal wiring. The declared pmech signal starts at + /// zero mechanical power, an admissible operating point. + TestOutcome hygov() + { + using Data = PhasorDynamics::Governor::HygovData; + using Outputs = typename Data::SignalOutputs; + using Params = typename Data::Parameters; + using Vars = PhasorDynamics::Governor::HygovInternalVariables; + + constexpr IdxT pmech_id = static_cast(2); + + TestStatus success = true; + + PhasorDynamics::SystemModelData data; + data.freq_base = 60.0; + data.va_base = 100.0e6; + data.signal.resize(1); + data.signal[0].signal_id = pmech_id; + data.signal[0].name = "Mechanical Power"; + + Data hygov_data; + hygov_data.device_class = "Hygov"; + hygov_data.disambiguation_string = "hygov_system"; + hygov_data.parameters[Params::Trate] = static_cast(100.0); + hygov_data.parameters[Params::Tnp] = static_cast(1.0); + hygov_data.signal_outputs[Outputs::pmech] = pmech_id; + data.hygov.push_back(hygov_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* pmech = system.getSignal(pmech_id); + success *= pmech->linked(); + success *= pmech->getVariableIndex() == static_cast(Vars::PMECH); + + auto missing_output_data = data; + missing_output_data.hygov[0].signal_outputs.clear(); + + PhasorDynamics::SystemModel missing_output_system(missing_output_data); + std::cout << "Testing expected HYGOV missing-output configuration error.\n"; + success *= missing_output_system.verify() > 0; + + return success.report(__func__); + } + private: auto makeRegcaData() -> PhasorDynamics::Converter::RegcaData { diff --git a/tests/UnitTests/PhasorDynamics/runComponentConnectionTests.cpp b/tests/UnitTests/PhasorDynamics/runComponentConnectionTests.cpp index b9e9253d1..127b5107b 100644 --- a/tests/UnitTests/PhasorDynamics/runComponentConnectionTests.cpp +++ b/tests/UnitTests/PhasorDynamics/runComponentConnectionTests.cpp @@ -8,6 +8,7 @@ int main() GridKit::Testing::ComponentConnectionTests test; result += test.genrouEsdc1a(); + result += test.genrouHygov(); return result.summary(); } diff --git a/tests/UnitTests/PhasorDynamics/runGovernorHygovTests.cpp b/tests/UnitTests/PhasorDynamics/runGovernorHygovTests.cpp new file mode 100644 index 000000000..3f69b4b90 --- /dev/null +++ b/tests/UnitTests/PhasorDynamics/runGovernorHygovTests.cpp @@ -0,0 +1,21 @@ +#include "GovernorHygovTests.hpp" + +int main() +{ + GridKit::Testing::TestingResults result; + + GridKit::Testing::GovernorHygovTests test; + + result += test.validation(); + result += test.initializationAndSignals(); + result += test.initializationDomain(); + result += test.initializationExactness(); + result += test.residualEquations(); + result += test.governorControl(); + result += test.turbineDynamics(); +#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 ebf13c6fb..abb274a0f 100644 --- a/tests/UnitTests/PhasorDynamics/runSystemSingleComponentTests.cpp +++ b/tests/UnitTests/PhasorDynamics/runSystemSingleComponentTests.cpp @@ -19,6 +19,7 @@ int main() result += test.genrou(); result += test.genClassical(); result += test.tgov1(); + result += test.hygov(); // @todo The following components are not tested here because they require non-trivial constructors // PhasorDynamics::Exciter::SexsPti diff --git a/tests/UnitTests/Utilities/CaseFormatTests.hpp b/tests/UnitTests/Utilities/CaseFormatTests.hpp index ef87d6bec..8416d2e1b 100644 --- a/tests/UnitTests/Utilities/CaseFormatTests.hpp +++ b/tests/UnitTests/Utilities/CaseFormatTests.hpp @@ -9,6 +9,7 @@ #include #include #include +#include #include #include #include @@ -189,6 +190,7 @@ namespace GridKit using BusData = BusData; using BusType = typename BusData::BusType; using Esdc1aData = Exciter::Esdc1aData; + using HygovData = Governor::HygovData; const char data[] = R"({ @@ -213,13 +215,20 @@ namespace GridKit { "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"} + { "signal_id": 6, "name": "Under-excitation Limiter"}, + { "signal_id": 7, "name": "Hydro Mechanical Power"}, + { "signal_id": 8, "name": "Governor Load Reference"}, + { "signal_id": 9, "name": "Governor Auxiliary Power"} ], "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": {"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": "Hygov", "ports": {"speed": 1, "pmech": 7, "pref": 8, "paux": 9}, "id": "DV6", "params": {"Trate": 80.0, "Rperm": 0.05, "Rtemp": 0.35, "Tr": 5.0, "Tf": 0.05, "Tg": 0.5, + "Velm": 0.2, "Gmax": 0.98, "Gmin": 0.02, "Tw": 1.2, "At": 1.1, "Dturb": 0.4, "Qnl": 0.08, "Tn": 0.7, "Tnp": 1.4, "db1": 0.01, "db2": 0.02, "Hdam": 1.05, + "Gv0": 0.0, "Gv1": 0.2, "Gv2": 0.4, "Gv3": 0.6, "Gv4": 0.8, "Gv5": 1.0, + "Pgv0": 0.0, "Pgv1": 0.15, "Pgv2": 0.42, "Pgv3": 0.66, "Pgv4": 0.85, "Pgv5": 1.0}, "mon": ["pmech", "filter", "desiredgate", "gate", "flow", "head"]}, { "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} } @@ -244,6 +253,7 @@ namespace GridKit success *= result.genrou.size() == 1; success *= result.gov.size() == 1; success *= result.esdc1a.size() == 1; + success *= result.hygov.size() == 1; success *= result.loadz.size() == 0; success *= result.exciter.size() == 1; success *= result.sexspti.size() == 1; @@ -276,6 +286,12 @@ namespace GridKit success *= result.signal[4].name == "Stabilizer Signal"; success *= result.signal[5].signal_id == 6; success *= result.signal[5].name == "Under-excitation Limiter"; + success *= result.signal[6].signal_id == 7; + success *= result.signal[6].name == "Hydro Mechanical Power"; + success *= result.signal[7].signal_id == 8; + success *= result.signal[7].name == "Governor Load Reference"; + success *= result.signal[8].signal_id == 9; + success *= result.signal[8].name == "Governor Auxiliary Power"; success *= std::get(result.branch[0].parameters[BranchParameters::R]) == 0.0; success *= std::get(result.branch[0].parameters[BranchParameters::X]) == 0.1; @@ -357,6 +373,49 @@ namespace GridKit success *= result.esdc1a[0].monitored_variables.contains(Esdc1aData::MonitorableVariables::se); success *= result.esdc1a[0].monitored_variables.contains(Esdc1aData::MonitorableVariables::vfe); + success *= std::get(result.hygov[0].parameters[HygovData::Parameters::Trate]) == 80.0; + success *= std::get(result.hygov[0].parameters[HygovData::Parameters::Rperm]) == 0.05; + success *= std::get(result.hygov[0].parameters[HygovData::Parameters::Rtemp]) == 0.35; + success *= std::get(result.hygov[0].parameters[HygovData::Parameters::Tr]) == 5.0; + success *= std::get(result.hygov[0].parameters[HygovData::Parameters::Tf]) == 0.05; + success *= std::get(result.hygov[0].parameters[HygovData::Parameters::Tg]) == 0.5; + success *= std::get(result.hygov[0].parameters[HygovData::Parameters::Velm]) == 0.2; + success *= std::get(result.hygov[0].parameters[HygovData::Parameters::Gmax]) == 0.98; + success *= std::get(result.hygov[0].parameters[HygovData::Parameters::Gmin]) == 0.02; + success *= std::get(result.hygov[0].parameters[HygovData::Parameters::Tw]) == 1.2; + success *= std::get(result.hygov[0].parameters[HygovData::Parameters::At]) == 1.1; + success *= std::get(result.hygov[0].parameters[HygovData::Parameters::Dturb]) == 0.4; + success *= std::get(result.hygov[0].parameters[HygovData::Parameters::Qnl]) == 0.08; + success *= std::get(result.hygov[0].parameters[HygovData::Parameters::Tn]) == 0.7; + success *= std::get(result.hygov[0].parameters[HygovData::Parameters::Tnp]) == 1.4; + success *= std::get(result.hygov[0].parameters[HygovData::Parameters::db1]) == 0.01; + success *= std::get(result.hygov[0].parameters[HygovData::Parameters::db2]) == 0.02; + success *= std::get(result.hygov[0].parameters[HygovData::Parameters::Hdam]) == 1.05; + success *= std::get(result.hygov[0].parameters[HygovData::Parameters::Gv0]) == 0.0; + success *= std::get(result.hygov[0].parameters[HygovData::Parameters::Gv1]) == 0.2; + success *= std::get(result.hygov[0].parameters[HygovData::Parameters::Gv2]) == 0.4; + success *= std::get(result.hygov[0].parameters[HygovData::Parameters::Gv3]) == 0.6; + success *= std::get(result.hygov[0].parameters[HygovData::Parameters::Gv4]) == 0.8; + success *= std::get(result.hygov[0].parameters[HygovData::Parameters::Gv5]) == 1.0; + success *= std::get(result.hygov[0].parameters[HygovData::Parameters::Pgv0]) == 0.0; + success *= std::get(result.hygov[0].parameters[HygovData::Parameters::Pgv1]) == 0.15; + success *= std::get(result.hygov[0].parameters[HygovData::Parameters::Pgv2]) == 0.42; + success *= std::get(result.hygov[0].parameters[HygovData::Parameters::Pgv3]) == 0.66; + success *= std::get(result.hygov[0].parameters[HygovData::Parameters::Pgv4]) == 0.85; + success *= std::get(result.hygov[0].parameters[HygovData::Parameters::Pgv5]) == 1.0; + success *= result.hygov[0].signal_inputs[HygovData::SignalInputs::speed] == 1; + success *= result.hygov[0].signal_inputs[HygovData::SignalInputs::pref] == 8; + success *= result.hygov[0].signal_inputs[HygovData::SignalInputs::paux] == 9; + success *= result.hygov[0].signal_outputs[HygovData::SignalOutputs::pmech] == 7; + success *= result.hygov[0].disambiguation_string == "DV6"; + success *= result.hygov[0].monitored_variables.contains(HygovData::MonitorableVariables::pmech); + success *= result.hygov[0].monitored_variables.contains(HygovData::MonitorableVariables::filter); + success *= result.hygov[0].monitored_variables.contains( + HygovData::MonitorableVariables::desiredgate); + success *= result.hygov[0].monitored_variables.contains(HygovData::MonitorableVariables::gate); + success *= result.hygov[0].monitored_variables.contains(HygovData::MonitorableVariables::flow); + success *= result.hygov[0].monitored_variables.contains(HygovData::MonitorableVariables::head); + 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;