From 62ad9b1770d2318ec309b0d866cba19d6f7df721 Mon Sep 17 00:00:00 2001 From: lukelowry Date: Sun, 17 May 2026 14:41:26 -0500 Subject: [PATCH 1/3] REGCA skeleton --- GridKit/Model/PhasorDynamics/CMakeLists.txt | 1 + .../Model/PhasorDynamics/ComponentLibrary.hpp | 1 + .../PhasorDynamics/Converter/CMakeLists.txt | 6 + .../Converter/REGCA/CMakeLists.txt | 49 +++ .../PhasorDynamics/Converter/REGCA/Regca.cpp | 27 ++ .../PhasorDynamics/Converter/REGCA/Regca.hpp | 168 ++++++++++ .../Converter/REGCA/RegcaData.hpp | 77 +++++ .../REGCA/RegcaDependencyTracking.cpp | 27 ++ .../Converter/REGCA/RegcaEnzyme.cpp | 31 ++ .../Converter/REGCA/RegcaImpl.hpp | 303 ++++++++++++++++++ GridKit/Model/PhasorDynamics/INPUT_FORMAT.md | 34 +- .../Model/PhasorDynamics/SystemModelData.hpp | 3 + .../SystemModelDataJSONParser.hpp | 6 + .../Model/PhasorDynamics/SystemModelImpl.hpp | 57 ++++ tests/UnitTests/PhasorDynamics/CMakeLists.txt | 10 + .../PhasorDynamics/ConverterRegcaTests.hpp | 232 ++++++++++++++ .../SystemSingleComponentTests.hpp | 24 ++ .../PhasorDynamics/runConverterRegcaTests.cpp | 16 + .../runSystemSingleComponentTests.cpp | 1 + 19 files changed, 1072 insertions(+), 1 deletion(-) create mode 100644 GridKit/Model/PhasorDynamics/Converter/CMakeLists.txt create mode 100644 GridKit/Model/PhasorDynamics/Converter/REGCA/CMakeLists.txt create mode 100644 GridKit/Model/PhasorDynamics/Converter/REGCA/Regca.cpp create mode 100644 GridKit/Model/PhasorDynamics/Converter/REGCA/Regca.hpp create mode 100644 GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaData.hpp create mode 100644 GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaDependencyTracking.cpp create mode 100644 GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaEnzyme.cpp create mode 100644 GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaImpl.hpp create mode 100644 tests/UnitTests/PhasorDynamics/ConverterRegcaTests.hpp create mode 100644 tests/UnitTests/PhasorDynamics/runConverterRegcaTests.cpp diff --git a/GridKit/Model/PhasorDynamics/CMakeLists.txt b/GridKit/Model/PhasorDynamics/CMakeLists.txt index 90ddb2abc..73d80e510 100644 --- a/GridKit/Model/PhasorDynamics/CMakeLists.txt +++ b/GridKit/Model/PhasorDynamics/CMakeLists.txt @@ -33,6 +33,7 @@ add_subdirectory(Branch) add_subdirectory(Bus) add_subdirectory(BusFault) add_subdirectory(BusToSignalAdapter) +add_subdirectory(Converter) add_subdirectory(Exciter) add_subdirectory(Governor) add_subdirectory(Load) diff --git a/GridKit/Model/PhasorDynamics/ComponentLibrary.hpp b/GridKit/Model/PhasorDynamics/ComponentLibrary.hpp index e2ae075d2..8d989fede 100644 --- a/GridKit/Model/PhasorDynamics/ComponentLibrary.hpp +++ b/GridKit/Model/PhasorDynamics/ComponentLibrary.hpp @@ -5,6 +5,7 @@ #include #include #include +#include #include #include #include diff --git a/GridKit/Model/PhasorDynamics/Converter/CMakeLists.txt b/GridKit/Model/PhasorDynamics/Converter/CMakeLists.txt new file mode 100644 index 000000000..cafb2cc36 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Converter/CMakeLists.txt @@ -0,0 +1,6 @@ +# [[ +# Author(s): +# - Luke Lowery +# ]] + +add_subdirectory(REGCA) diff --git a/GridKit/Model/PhasorDynamics/Converter/REGCA/CMakeLists.txt b/GridKit/Model/PhasorDynamics/Converter/REGCA/CMakeLists.txt new file mode 100644 index 000000000..9e9b582d7 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Converter/REGCA/CMakeLists.txt @@ -0,0 +1,49 @@ +# [[ +# Author(s): +# - Luke Lowery +# ]] + +set(_install_headers + Regca.hpp + RegcaData.hpp) + +if(GRIDKIT_ENABLE_ENZYME) + gridkit_add_library(phasor_dynamics_converter_regca + SOURCES + RegcaEnzyme.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_converter_regca + SOURCES + Regca.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_converter_regca_dependency_tracking + SOURCES + RegcaDependencyTracking.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_converter_regca) +target_link_libraries(phasor_dynamics_components_dependency_tracking + INTERFACE GridKit::phasor_dynamics_converter_regca_dependency_tracking) diff --git a/GridKit/Model/PhasorDynamics/Converter/REGCA/Regca.cpp b/GridKit/Model/PhasorDynamics/Converter/REGCA/Regca.cpp new file mode 100644 index 000000000..b018c4bc3 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Converter/REGCA/Regca.cpp @@ -0,0 +1,27 @@ +/** + * @file Regca.cpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Non-Enzyme instantiation for the REGCA converter model. + */ + +#include "RegcaImpl.hpp" + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Converter + { + template + int Regca::evaluateJacobian() + { + Log::misc() << "Evaluate Jacobian for Regca..." << std::endl; + Log::misc() << "Jacobian evaluation is not implemented!" << std::endl; + return 0; + } + + template class Regca; + template class Regca; + } // namespace Converter + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Converter/REGCA/Regca.hpp b/GridKit/Model/PhasorDynamics/Converter/REGCA/Regca.hpp new file mode 100644 index 000000000..cbafdc7d8 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Converter/REGCA/Regca.hpp @@ -0,0 +1,168 @@ +/** + * @file Regca.hpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Declaration of the REGCA phasor-dynamics converter model. + */ + +#pragma once + +#include +#include +#include + +#include +#include +#include +#include + +namespace GridKit +{ + namespace PhasorDynamics + { + template + class BusBase; + + template + class SignalNode; + } // namespace PhasorDynamics +} // namespace GridKit + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Converter + { + /// Internal variables of a `Regca` + enum class RegcaInternalVariables : size_t + { + VM, ///< Filtered terminal voltage + IQ, ///< Reactive-current state + IP, ///< Active-current state + VT, ///< Terminal voltage magnitude + II, ///< Imaginary injected current + IQEXTRA, ///< HVRCM extra reactive current + IL, ///< LVPL upper-limit current curve + IR, ///< Real injected current + LP, ///< Active-current lower rate bound + UP, ///< Active-current upper rate bound + MAXIMUM, + }; + + /// External variables of a `Regca` + enum class RegcaExternalVariables : size_t + { + IPCMD, ///< Active-current command signal + IQCMD, ///< Reactive-current command signal + MAXIMUM, + }; + + template + class Regca : public Component + { + using Component::gridkit_component_id_; + using Component::alpha_; + using Component::f_; + using Component::h_; + using Component::J_; + using Component::J_cols_buffer_; + using Component::J_rows_buffer_; + using Component::J_vals_buffer_; + using Component::mva_system_base_; + using Component::nnz_; + using Component::residual_indices_; + using Component::size_; + using Component::tag_; + using Component::time_; + using Component::variable_indices_; + using Component::wb_; + using Component::y_; + using Component::yp_; + + public: + using RealT = typename Component::RealT; + using bus_type = BusBase; + using signal_type = SignalNode; + using model_data_type = RegcaData; + using MonitorT = Model::VariableMonitor; + + Regca(bus_type* bus); + Regca(bus_type* bus, const model_data_type& data); + ~Regca(); + + int setGridKitComponentID(IdxT) override final; + int allocate() override final; + int verify() const override final; + int initialize() override final; + int tagDifferentiable() 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( + ScalarT*, ScalarT*, ScalarT*, ScalarT*, ScalarT*); + __attribute__((always_inline)) inline int evaluateBusResidual( + ScalarT*, ScalarT*, ScalarT*, ScalarT*); + + private: + void initializeParameters(const model_data_type& data); + void initializeMonitor(); + + ScalarT& Vr() + { + return bus_->Vr(); + } + + ScalarT& Vi() + { + return bus_->Vi(); + } + + ScalarT& Ir() + { + return bus_->Ir(); + } + + ScalarT& Ii() + { + return bus_->Ii(); + } + + bus_type* bus_{nullptr}; + + RealT P0_{0}; + RealT Q0_{0}; + RealT Sconv_{0}; + RealT Tg_{0}; + RealT TM_{0}; + RealT Rqmax_{0}; + RealT Rqmin_{0}; + RealT Rpmax_{0}; + bool sL_{false}; + RealT IL1_{0}; + RealT VL0_{0}; + RealT VL1_{0}; + RealT VA0_{0}; + RealT VA1_{0}; + RealT Vhvmax_{0}; + IdxT bus_id_{0}; + + ComponentSignals signals_; + std::unique_ptr monitor_; + + std::vector ws_; + std::vector ws_indices_; + }; + } // namespace Converter + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaData.hpp b/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaData.hpp new file mode 100644 index 000000000..29f8530fe --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaData.hpp @@ -0,0 +1,77 @@ +/** + * @file RegcaData.hpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Modeling data for the REGCA converter model. + */ + +#pragma once + +#include + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Converter + { + /// Parameter keys for the REGCA converter model. + enum class RegcaParameters + { + P0, ///< Initial active power injection on system base + Q0, ///< Initial reactive power injection on system base + Sconv, ///< Converter/model power base + Tg, ///< Converter current-control lag time constant + TM, ///< Terminal voltage sensor time constant + Rqmax, ///< Reactive-current recovery positive rate limit + Rqmin, ///< Reactive-current recovery negative rate limit + Rpmax, ///< Active-current magnitude recovery rate limit + sL, ///< LVPL switch + IL1, ///< LVPL upper-current ceiling + VL0, ///< LVPL zero-crossing voltage + VL1, ///< LVPL upper breakpoint voltage + VA0, ///< LVACM lower breakpoint voltage + VA1, ///< LVACM upper breakpoint voltage + Vhvmax ///< Terminal-voltage ceiling for HV reactive management + }; + + /// Ports for the REGCA converter model. + enum class RegcaPorts + { + bus, ///< Terminal bus ID + ipcmd, ///< Optional active-current command signal ID + iqcmd ///< Optional reactive-current command signal ID + }; + + /// Variables available through the monitor interface. + enum class RegcaMonitorableVariables + { + ir, ///< Real current injection on converter base + ii, ///< Imaginary current injection on converter base + p, ///< Active power injection on system base + q, ///< Reactive power injection on system base + vt, ///< Terminal voltage magnitude + vm, ///< Filtered terminal voltage + ip, ///< Active-current state + iq, ///< Reactive-current state + iqextra, ///< HVRCM extra reactive current + il, ///< LVPL upper-limit current curve + lp, ///< Active-current lower rate bound + up ///< Active-current upper rate bound + }; + + template + struct RegcaData : public ComponentData + { + RegcaData() = default; + + using Parameters = RegcaParameters; + using Ports = RegcaPorts; + using MonitorableVariables = RegcaMonitorableVariables; + }; + } // namespace Converter + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaDependencyTracking.cpp b/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaDependencyTracking.cpp new file mode 100644 index 000000000..20d688688 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaDependencyTracking.cpp @@ -0,0 +1,27 @@ +/** + * @file RegcaDependencyTracking.cpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Dependency-tracking instantiations for the REGCA converter model. + */ + +#include "RegcaImpl.hpp" + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Converter + { + template + int Regca::evaluateJacobian() + { + Log::misc() << "Evaluate Jacobian for Regca..." << std::endl; + Log::misc() << "Jacobian evaluation is not implemented!" << std::endl; + return 0; + } + + template class Regca; + template class Regca; + } // namespace Converter + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaEnzyme.cpp b/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaEnzyme.cpp new file mode 100644 index 000000000..7ad8ee7fc --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaEnzyme.cpp @@ -0,0 +1,31 @@ +/** + * @file RegcaEnzyme.cpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Enzyme sparse Jacobian for the REGCA converter model. + */ + +#include + +#include "RegcaImpl.hpp" + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Converter + { + template + int Regca::evaluateJacobian() + { + Log::misc() << "Evaluate Jacobian for Regca..." << std::endl; + Log::misc() << "Jacobian evaluation is experimental!" << std::endl; + + J_.zeroMatrix(); + return 0; + } + + template class Regca; + template class Regca; + } // namespace Converter + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaImpl.hpp b/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaImpl.hpp new file mode 100644 index 000000000..6598ea24e --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaImpl.hpp @@ -0,0 +1,303 @@ +/** + * @file RegcaImpl.hpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Definition of the REGCA phasor-dynamics converter model. + */ + +#pragma once + +#include +#include + +#include +#include +#include +#include +#include +#include + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Converter + { + using Log = ::GridKit::Utilities::Logger; + + template + Regca::Regca(bus_type* bus) + : bus_(bus) + { + size_ = static_cast(RegcaInternalVariables::MAXIMUM); + } + + template + Regca::Regca(bus_type* bus, const model_data_type& data) + : bus_(bus), + monitor_(std::make_unique(data)) + { + initializeParameters(data); + initializeMonitor(); + size_ = static_cast(RegcaInternalVariables::MAXIMUM); + } + + template + Regca::~Regca() + { + } + + template + void Regca::initializeParameters(const model_data_type& data) + { + using Params = typename model_data_type::Parameters; + using Ports = typename model_data_type::Ports; + + auto loadReal = [&](auto key, RealT& target) + { + if (!data.parameters.contains(key)) + { + return; + } + + const auto& value = data.parameters.at(key); + if (const auto* real_value = std::get_if(&value)) + { + target = *real_value; + } + else if (const auto* index_value = std::get_if(&value)) + { + target = static_cast(*index_value); + } + }; + + auto loadBool = [&](auto key, bool& target) + { + if (!data.parameters.contains(key)) + { + return; + } + + const auto& value = data.parameters.at(key); + if (const auto* bool_value = std::get_if(&value)) + { + target = *bool_value; + } + else if (const auto* index_value = std::get_if(&value)) + { + target = (*index_value != 0); + } + }; + + loadReal(Params::P0, P0_); + loadReal(Params::Q0, Q0_); + loadReal(Params::Sconv, Sconv_); + loadReal(Params::Tg, Tg_); + loadReal(Params::TM, TM_); + loadReal(Params::Rqmax, Rqmax_); + loadReal(Params::Rqmin, Rqmin_); + loadReal(Params::Rpmax, Rpmax_); + loadBool(Params::sL, sL_); + loadReal(Params::IL1, IL1_); + loadReal(Params::VL0, VL0_); + loadReal(Params::VL1, VL1_); + loadReal(Params::VA0, VA0_); + loadReal(Params::VA1, VA1_); + loadReal(Params::Vhvmax, Vhvmax_); + + if (data.ports.contains(Ports::bus)) + { + bus_id_ = data.ports.at(Ports::bus); + } + } + + template + const Model::VariableMonitorBase* Regca::getMonitor() const + { + return monitor_.get(); + } + + template + void Regca::initializeMonitor() + { + using Variable = typename model_data_type::MonitorableVariables; + auto index = [](RegcaInternalVariables variable) + { + return static_cast(variable); + }; + + monitor_->set(Variable::ir, [this, index] + { return y_[index(RegcaInternalVariables::IR)]; }); + monitor_->set(Variable::ii, [this, index] + { return y_[index(RegcaInternalVariables::II)]; }); + monitor_->set(Variable::p, [] + { return ScalarT{0}; }); + monitor_->set(Variable::q, [] + { return ScalarT{0}; }); + monitor_->set(Variable::vt, [this, index] + { return y_[index(RegcaInternalVariables::VT)]; }); + monitor_->set(Variable::vm, [this, index] + { return y_[index(RegcaInternalVariables::VM)]; }); + monitor_->set(Variable::ip, [this, index] + { return y_[index(RegcaInternalVariables::IP)]; }); + monitor_->set(Variable::iq, [this, index] + { return y_[index(RegcaInternalVariables::IQ)]; }); + monitor_->set(Variable::iqextra, [this, index] + { return y_[index(RegcaInternalVariables::IQEXTRA)]; }); + monitor_->set(Variable::il, [this, index] + { return y_[index(RegcaInternalVariables::IL)]; }); + monitor_->set(Variable::lp, [this, index] + { return y_[index(RegcaInternalVariables::LP)]; }); + monitor_->set(Variable::up, [this, index] + { return y_[index(RegcaInternalVariables::UP)]; }); + } + + template + int Regca::setGridKitComponentID(IdxT component_id) + { + gridkit_component_id_ = component_id; + return 0; + } + + template + int Regca::allocate() + { + size_ = static_cast(RegcaInternalVariables::MAXIMUM); + auto size = static_cast(size_); + + f_.assign(size, ScalarT{0}); + y_.assign(size, ScalarT{0}); + yp_.assign(size, ScalarT{0}); + tag_.assign(size, false); + variable_indices_.resize(size); + residual_indices_.resize(size); + + wb_.assign(2, ScalarT{0}); + h_.assign(2, ScalarT{0}); + + auto signal_size = static_cast(RegcaExternalVariables::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); + } + + return 0; + } + + template + int Regca::verify() const + { + int ret = 0; + + if (bus_ == nullptr) + { + Log::error() << "Regca: bus pointer is null\n"; + ret += 1; + } + + if (signals_.template isAttached()) + { + if (!signals_.template isLinked()) + { + Log::error() << "Regca: ipcmd signal attached with no linked source\n"; + ret += 1; + } + } + + if (signals_.template isAttached()) + { + if (!signals_.template isLinked()) + { + Log::error() << "Regca: iqcmd signal attached with no linked source\n"; + ret += 1; + } + } + + return ret; + } + + template + int Regca::initialize() + { + std::fill(y_.begin(), y_.end(), ScalarT{0}); + std::fill(yp_.begin(), yp_.end(), ScalarT{0}); + return 0; + } + + template + int Regca::tagDifferentiable() + { + std::fill(tag_.begin(), tag_.end(), false); + tag_[static_cast(RegcaInternalVariables::VM)] = true; + tag_[static_cast(RegcaInternalVariables::IQ)] = true; + tag_[static_cast(RegcaInternalVariables::IP)] = true; + return 0; + } + + template + __attribute__((always_inline)) inline int Regca::evaluateInternalResidual( + [[maybe_unused]] ScalarT* y, + [[maybe_unused]] ScalarT* yp, + [[maybe_unused]] ScalarT* wb, + [[maybe_unused]] ScalarT* ws, + ScalarT* f) + { + for (IdxT i = 0; i < size_; ++i) + { + f[static_cast(i)] = ScalarT{0}; + } + + return 0; + } + + template + __attribute__((always_inline)) inline int Regca::evaluateBusResidual( + [[maybe_unused]] ScalarT* y, + [[maybe_unused]] ScalarT* yp, + [[maybe_unused]] ScalarT* wb, + ScalarT* h) + { + h[0] = ScalarT{0}; + h[1] = ScalarT{0}; + return 0; + } + + template + int Regca::evaluateResidual() + { + std::fill(ws_.begin(), ws_.end(), ScalarT{0}); + std::fill(ws_indices_.begin(), ws_indices_.end(), INVALID_INDEX); + + if (signals_.template isAttached()) + { + const auto index = static_cast(RegcaExternalVariables::IPCMD); + ws_[index] = signals_.template readExternalVariable(); + ws_indices_[index] = + signals_.template readExternalVariableIndex(); + } + + if (signals_.template isAttached()) + { + const auto index = static_cast(RegcaExternalVariables::IQCMD); + ws_[index] = signals_.template readExternalVariable(); + ws_indices_[index] = + signals_.template readExternalVariableIndex(); + } + + wb_[0] = Vr(); + wb_[1] = Vi(); + + evaluateInternalResidual(y_.data(), yp_.data(), wb_.data(), ws_.data(), f_.data()); + evaluateBusResidual(y_.data(), yp_.data(), wb_.data(), h_.data()); + + Ir() += h_[0]; + Ii() += h_[1]; + + return 0; + } + } // namespace Converter + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/INPUT_FORMAT.md b/GridKit/Model/PhasorDynamics/INPUT_FORMAT.md index 32df36cdc..48db096f8 100644 --- a/GridKit/Model/PhasorDynamics/INPUT_FORMAT.md +++ b/GridKit/Model/PhasorDynamics/INPUT_FORMAT.md @@ -144,16 +144,48 @@ are specified: `Genrou` | 6th order machine model | `bus`, `pmech`\*, `speed`\*, `efd`\* | `p0`, `q0`, `H`, `D`, `Ra`, `Tdop`, `Tdopp`, `Tqop`, `Tqopp`, `Xd`, `Xdp`, `Xdpp`, `Xq`, `Xqp`, `Xqpp`, `Xl`, `S10`, `S12`, `mva` | `ir`, `ii`, `p`, `q`, `delta`, `omega`, `speed` `Gensal` | 5th order salient-pole machine model | `bus`, `pmech`\*, `speed`\*, `efd`\* | `p0`, `q0`, `H`, `D`, `Ra`, `Tdop`, `Tdopp`, `Tqopp`, `Xd`, `Xdp`, `Xdpp`, `Xq`, `Xl`, `S10`, `S12`, `mva` | `ir`, `ii`, `p`, `q`, `delta`, `omega`, `speed`, `Eqp`, `psidp`, `psiqpp`, `psidpp`, `vd`, `vq`, `te`, `id`, `iq` `GenClassical` | the classical machine model | `bus`, `pmech`\*, `speed`\*, `efd`\* | `p0`, `q0`, `H`, `D`, `Ra`, `Xdp`, `mva` | `ir`, `ii`, `p`, `q`, `delta`, `omega` + `Regca` | WECC REGCA renewable generator/converter skeleton | `bus`, `ipcmd`\*, `iqcmd`\* | `P0`, `Q0`, `Sconv`, `Tg`, `TM`, `Rqmax`, `Rqmin`, `Rpmax`, `sL`, `IL1`, `VL0`, `VL1`, `VA0`, `VA1`, `Vhvmax` | `ir`, `ii`, `p`, `q`, `vt`, `vm`, `ip`, `iq`, `iqextra`, `il`, `lp`, `up` `Tgov1 ` | the TGOV1 governor model | `pmech`, `speed` | `R`, `T1`, `T2`, `T3`, `Pvmax`, `Pvmin`, `Dt` | `none` `Ieeet1` | the IEEET1 exciter model | `bus`, `speed`, `efd`, `vs`\* | `Tr`, `Ka`, `Ta`, `Ke`, `Te`, `Kf`, `Tf`, `Vrmin`, `Vrmax`, `E1`, `E2`, `Se1`, `Se2`, `Ispdlim` | `efd`, `ksat` `SexsPti` | the SEXS-PTI simplified exciter model | `bus`, `efd`, `vs`\* | `Ta`, `Tb`, `Te`, `K`, `Efdmax`, `Efdmin` | `efd` - `Ieeest` | the IEEEST stabilizer model | `input`, `output` | `A1`, `A2`, `A3`, `A4`, `A5`, `A6`, `T1`, `T2`, `T3`, `T4`, `T5`, `T6`, `Ks`, `Lsmin`, `Lsmax`, `Vcl`, `Vcu`, `Tdelay` | `vss` + `Ieeest` | the IEEEST stabilizer model | `input`, `output` | `A1`, `A2`, `A3`, `A4`, `A5`, `A6`, `T1`, `T2`, `T3`, `T4`, `T5`, `T6`, `Ks`, `Lsmin`, `Lsmax`, `Vcl`, `Vcu`, `Tdelay` | `vss` `BusFault` | simple impedance-based fault at a bus | `bus`, `status`\* | `state0`, `R`, `X` | `state`, `ir`, `ii` `BusToSignalAdapter` | signal adapter component for a bus | `bus`, `vr`, `vi`, `ir`, `ii` | | Ports marked with \* are optional and, if missing, will be assumed to be connected to a constant value. This list is subject to change. +For `Regca`, `ipcmd` and `iqcmd` are optional in the skeleton implementation +and default to constant zero command inputs when omitted. + +Compact `Regca` device example: + +```json +{ + "class": "Regca", + "ports": { "bus": 1, "ipcmd": 10, "iqcmd": 11 }, + "id": "CV1", + "params": { + "P0": 1.0, + "Q0": 0.0, + "Sconv": 100.0, + "Tg": 0.02, + "TM": 0.02, + "Rqmax": 999.0, + "Rqmin": -999.0, + "Rpmax": 999.0, + "sL": true, + "IL1": 1.1, + "VL0": 0.4, + "VL1": 0.9, + "VA0": 0.4, + "VA1": 0.9, + "Vhvmax": 1.2 + }, + "mon": ["ir", "ii"] +} +``` + For `Branch`, `tap` and `phase` are optional parameters. If omitted, `tap` defaults to `1.0` and `phase` defaults to `0.0` radians. Bus `bus1` is the tap diff --git a/GridKit/Model/PhasorDynamics/SystemModelData.hpp b/GridKit/Model/PhasorDynamics/SystemModelData.hpp index d4deef66e..7cb2765c5 100644 --- a/GridKit/Model/PhasorDynamics/SystemModelData.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModelData.hpp @@ -10,6 +10,7 @@ #include #include #include +#include #include #include #include @@ -40,6 +41,7 @@ namespace GridKit using BusDataT = BusData; using BusToSignalAdapterDataT = BusToSignalAdapterData; using BusFaultDataT = BusFaultData; + using RegcaDataT = Converter::RegcaData; using Tgov1DataT = Governor::Tgov1Data; using Ieeet1DataT = Exciter::Ieeet1Data; using SexsPtiDataT = Exciter::SexsPtiData; @@ -93,6 +95,7 @@ namespace GridKit std::vector adapter; ///< bus-to-signal adapters within the model std::vector branch; ///< Branches within the model std::vector bus_fault; ///< Bus faults within the model + std::vector regca; ///< REGCA converter instances within the model std::vector genrou; ///< GENROU instances within the model std::vector gensal; ///< GENSAL instances within the model std::vector genclassical; ///< Classical generator instances within the model diff --git a/GridKit/Model/PhasorDynamics/SystemModelDataJSONParser.hpp b/GridKit/Model/PhasorDynamics/SystemModelDataJSONParser.hpp index 00a6278c2..b215c0c86 100644 --- a/GridKit/Model/PhasorDynamics/SystemModelDataJSONParser.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModelDataJSONParser.hpp @@ -136,6 +136,12 @@ namespace GridKit raw_component.get_to(loadzip); sm.loadzip.push_back(loadzip); } + else if (kind == "Regca") + { + typename SystemModelData::RegcaDataT regca; + raw_component.get_to(regca); + sm.regca.push_back(regca); + } else if (kind == "Tgov1") { typename SystemModelData::Tgov1DataT gov; diff --git a/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp b/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp index 891c2c57d..c28e7864f 100644 --- a/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp @@ -46,6 +46,7 @@ namespace GridKit using namespace Governor; using namespace Exciter; using namespace Stabilizer; + using namespace Converter; // Set system model tolerances rel_tol_ = 1e-7; @@ -115,6 +116,62 @@ namespace GridKit addComponent(adapter); } + // Add REGCA converters + for (const auto& regcadata : data.regca) + { + IdxT bus_index = 0; + if (regcadata.ports.contains(RegcaPorts::bus)) + { + bus_index = regcadata.ports.at(RegcaPorts::bus); + } + + auto* regca = new Regca(getBus(bus_index), regcadata); + + if (regcadata.ports.contains(RegcaPorts::ipcmd)) + { + IdxT ipcmd = regcadata.ports.at(RegcaPorts::ipcmd); + constexpr auto IPCMD = RegcaExternalVariables::IPCMD; + regca->getSignals().template attachSignalNode(getSignal(ipcmd)); + } + + if (regcadata.ports.contains(RegcaPorts::iqcmd)) + { + IdxT iqcmd = regcadata.ports.at(RegcaPorts::iqcmd); + constexpr auto IQCMD = RegcaExternalVariables::IQCMD; + regca->getSignals().template attachSignalNode(getSignal(iqcmd)); + } + + if (regcadata.ports.contains(RegcaPorts::ibranchr)) + { + IdxT ibranchr = regcadata.ports.at(RegcaPorts::ibranchr); + constexpr auto IR = RegcaInternalVariables::IR; + regca->getSignals().template assignSignalNode(getSignal(ibranchr)); + } + + if (regcadata.ports.contains(RegcaPorts::ibranchi)) + { + IdxT ibranchi = regcadata.ports.at(RegcaPorts::ibranchi); + constexpr auto II = RegcaInternalVariables::II; + regca->getSignals().template assignSignalNode(getSignal(ibranchi)); + } + + if (regcadata.ports.contains(RegcaPorts::pbranch)) + { + IdxT pbranch = regcadata.ports.at(RegcaPorts::pbranch); + constexpr auto PBR = RegcaInternalVariables::PBR; + regca->getSignals().template assignSignalNode(getSignal(pbranch)); + } + + if (regcadata.ports.contains(RegcaPorts::qbranch)) + { + IdxT qbranch = regcadata.ports.at(RegcaPorts::qbranch); + constexpr auto QBR = RegcaInternalVariables::QBR; + regca->getSignals().template assignSignalNode(getSignal(qbranch)); + } + + addComponent(regca); + } + // Add branches for (const auto& branchdata : data.branch) { diff --git a/tests/UnitTests/PhasorDynamics/CMakeLists.txt b/tests/UnitTests/PhasorDynamics/CMakeLists.txt index bac53ec7b..ed5f171a3 100644 --- a/tests/UnitTests/PhasorDynamics/CMakeLists.txt +++ b/tests/UnitTests/PhasorDynamics/CMakeLists.txt @@ -94,6 +94,14 @@ target_link_libraries( GridKit::phasor_dynamics_systemmodel_dependency_tracking GridKit::testing) +add_executable(test_phasor_converter_regca runConverterRegcaTests.cpp) +target_link_libraries( + test_phasor_converter_regca + GridKit::definitions + GridKit::phasor_dynamics_components + GridKit::phasor_dynamics_components_dependency_tracking + GridKit::testing) + add_executable(test_phasor_stabilizer_ieeest runStabilizerIeeestTests.cpp) target_link_libraries( test_phasor_stabilizer_ieeest @@ -138,6 +146,7 @@ add_test(NAME PhasorDynamicsGovernorTgov1Test COMMAND test_phasor_governor_tgov1 add_test(NAME PhasorDynamicsExciterIeeet1Test COMMAND test_phasor_exciter_ieeet1) add_test(NAME PhasorDynamicsGensalTest COMMAND test_phasor_gensal) add_test(NAME PhasorDynamicsExciterSexsPtiTest COMMAND test_phasor_exciter_sexspti) +add_test(NAME PhasorDynamicsConverterRegcaTest COMMAND test_phasor_converter_regca) add_test(NAME PhasorDynamicsStabilizerIeeestTest COMMAND test_phasor_stabilizer_ieeest) add_test(NAME PhasorDynamicsGenClassicalTest COMMAND test_phasor_gen_classical) add_test(NAME PhasorDynamicsLoadTest COMMAND test_phasor_load) @@ -159,6 +168,7 @@ install( test_phasor_exciter_ieeet1 test_phasor_gensal test_phasor_exciter_sexspti + test_phasor_converter_regca test_phasor_stabilizer_ieeest test_phasor_gen_classical test_phasor_system diff --git a/tests/UnitTests/PhasorDynamics/ConverterRegcaTests.hpp b/tests/UnitTests/PhasorDynamics/ConverterRegcaTests.hpp new file mode 100644 index 000000000..f74da86e6 --- /dev/null +++ b/tests/UnitTests/PhasorDynamics/ConverterRegcaTests.hpp @@ -0,0 +1,232 @@ +#pragma once + +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace GridKit +{ + namespace Testing + { + template + class ConverterRegcaTests + { + public: + using RealT = typename PhasorDynamics::Component::RealT; + + ConverterRegcaTests() = default; + ~ConverterRegcaTests() = default; + + static constexpr ScalarT kTol = static_cast(1.0e-14); + + TestOutcome constructor() + { + TestStatus success = true; + + PhasorDynamics::Bus bus(1.0, 0.0); + + PhasorDynamics::Converter::Regca minimal(&bus); + success *= (minimal.size() == static_cast(PhasorDynamics::Converter::RegcaInternalVariables::MAXIMUM)); + success *= (minimal.getMonitor() == nullptr); + + auto data = makeTestData(); + PhasorDynamics::Converter::Regca from_data(&bus, data); + success *= (from_data.size() == static_cast(PhasorDynamics::Converter::RegcaInternalVariables::MAXIMUM)); + success *= (from_data.getMonitor() != nullptr); + + return success.report(__func__); + } + + TestOutcome lifecycle() + { + TestStatus success = true; + + PhasorDynamics::Bus bus(1.0, 0.0); + bus.allocate(); + bus.initialize(); + bus.evaluateResidual(); + + auto data = makeTestData(); + PhasorDynamics::Converter::Regca regca(&bus, data); + success *= (regca.allocate() == 0); + success *= (regca.initialize() == 0); + success *= (regca.tagDifferentiable() == 0); + success *= (regca.evaluateResidual() == 0); + success *= (regca.evaluateJacobian() == 0); + + const auto& y = regca.y(); + const auto& yp = regca.yp(); + const auto& f = regca.getResidual(); + for (size_t i = 0; i < static_cast(regca.size()); ++i) + { + success *= isEqual(y[i], static_cast(0.0), kTol); + success *= isEqual(yp[i], static_cast(0.0), kTol); + success *= isEqual(f[i], static_cast(0.0), kTol); + } + + using Vars = PhasorDynamics::Converter::RegcaInternalVariables; + for (size_t i = 0; i < static_cast(Vars::MAXIMUM); ++i) + { + const bool expected = (i == static_cast(Vars::VM)) + || (i == static_cast(Vars::IQ)) + || (i == static_cast(Vars::IP)); + success *= (regca.tag()[i] == expected); + } + + success *= isEqual(bus.Ir(), static_cast(0.0), kTol); + success *= isEqual(bus.Ii(), static_cast(0.0), kTol); + + return success.report(__func__); + } + + TestOutcome signalVerification() + { + TestStatus success = true; + + PhasorDynamics::Bus bus(1.0, 0.0); + PhasorDynamics::Converter::Regca regca(&bus, makeTestData()); + + PhasorDynamics::SignalNode ipcmd_node; + PhasorDynamics::SignalNode iqcmd_node; + ScalarT ipcmd_value{0.25}; + ScalarT iqcmd_value{-0.10}; + IdxT ipcmd_index = 21; + IdxT iqcmd_index = 22; + + regca.getSignals().template attachSignalNode(&ipcmd_node); + success *= (regca.verify() > 0); + + ipcmd_node.set(&ipcmd_value, &ipcmd_index); + success *= (regca.verify() == 0); + + regca.getSignals().template attachSignalNode(&iqcmd_node); + success *= (regca.verify() > 0); + + iqcmd_node.set(&iqcmd_value, &iqcmd_index); + success *= (regca.verify() == 0); + + return success.report(__func__); + } + + TestOutcome nullBusVerification() + { + TestStatus success = true; + + PhasorDynamics::Converter::Regca regca(nullptr); + success *= (regca.verify() > 0); + + return success.report(__func__); + } + + TestOutcome jsonParseAndSystemAssembly() + { + TestStatus success = true; + + std::istringstream input(R"json( +{ + "header": { + "format_version": 0, + "format_revision": 1, + "case_name": "REGCA skeleton", + "case_description": "REGCA parser smoke test", + "case_comments": "", + "freq_base": 60.0, + "va_base": 100000000.0 + }, + "buses": [ + { "number": 1, "class": "bus", "name": "Bus 1", "init": { "Vr": 1.0, "Vi": 0.0 }, "v_base": 1.0 } + ], + "devices": [ + { + "class": "Regca", + "ports": { "bus": 1 }, + "id": "CV1", + "params": { + "P0": 1.0, + "Q0": 0.0, + "Sconv": 100, + "Tg": 0.02, + "TM": 0.02, + "Rqmax": 999.0, + "Rqmin": -999.0, + "Rpmax": 999.0, + "sL": true, + "IL1": 1.1, + "VL0": 0.4, + "VL1": 0.9, + "VA0": 0.4, + "VA1": 0.9, + "Vhvmax": 1.2 + }, + "mon": ["ir", "ii"] + } + ] +} +)json"); + + auto data = PhasorDynamics::parseSystemModelData(input); + success *= (data.regca.size() == 1); + success *= (data.regca[0].device_class == "Regca"); + success *= (data.regca[0].ports.at(PhasorDynamics::Converter::RegcaPorts::bus) == 1); + success *= (std::get_if(&data.regca[0].parameters.at(PhasorDynamics::Converter::RegcaParameters::Sconv)) != nullptr); + success *= (std::get_if(&data.regca[0].parameters.at(PhasorDynamics::Converter::RegcaParameters::sL)) != nullptr); + + 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() == 12); + success *= isEqual(system.getResidual()[0], 0.0, static_cast(kTol)); + success *= isEqual(system.getResidual()[1], 0.0, static_cast(kTol)); + + return success.report(__func__); + } + + private: + auto makeTestData() -> PhasorDynamics::Converter::RegcaData + { + using Params = PhasorDynamics::Converter::RegcaParameters; + using Ports = PhasorDynamics::Converter::RegcaPorts; + using Mon = PhasorDynamics::Converter::RegcaMonitorableVariables; + + PhasorDynamics::Converter::RegcaData data; + data.device_class = "Regca"; + data.disambiguation_string = "regca_test"; + data.ports[Ports::bus] = 1; + data.monitored_variables.insert(Mon::ir); + data.monitored_variables.insert(Mon::ii); + + data.parameters[Params::P0] = static_cast(1.0); + data.parameters[Params::Q0] = static_cast(0.0); + data.parameters[Params::Sconv] = static_cast(100); + data.parameters[Params::Tg] = static_cast(0.02); + data.parameters[Params::TM] = static_cast(0.02); + data.parameters[Params::Rqmax] = static_cast(999.0); + data.parameters[Params::Rqmin] = static_cast(-999.0); + data.parameters[Params::Rpmax] = static_cast(999.0); + data.parameters[Params::sL] = true; + data.parameters[Params::IL1] = static_cast(1.1); + data.parameters[Params::VL0] = static_cast(0.4); + data.parameters[Params::VL1] = static_cast(0.9); + data.parameters[Params::VA0] = static_cast(0.4); + data.parameters[Params::VA1] = static_cast(0.9); + data.parameters[Params::Vhvmax] = static_cast(1.2); + + return data; + } + }; + } // namespace Testing +} // namespace GridKit diff --git a/tests/UnitTests/PhasorDynamics/SystemSingleComponentTests.hpp b/tests/UnitTests/PhasorDynamics/SystemSingleComponentTests.hpp index c32df9f58..790d7048a 100644 --- a/tests/UnitTests/PhasorDynamics/SystemSingleComponentTests.hpp +++ b/tests/UnitTests/PhasorDynamics/SystemSingleComponentTests.hpp @@ -166,6 +166,30 @@ namespace GridKit return success.report(__func__); } + TestOutcome regca() + { + TestStatus success = true; + + PhasorDynamics::SystemModel* system = new PhasorDynamics::SystemModel(); + + PhasorDynamics::BusInfinite bus; + system->addBus(&bus); + + PhasorDynamics::Converter::Regca regca(&bus); + system->addComponent(®ca); + + success *= system->allocate() == 0; + success *= system->initialize() == 0; + success *= system->evaluateResidual() == 0; + success *= system->evaluateJacobian() == 0; + success *= system->size() == regca.size(); + + delete system; + system = nullptr; + + return success.report(__func__); + } + TestOutcome genrou() { TestStatus success = true; diff --git a/tests/UnitTests/PhasorDynamics/runConverterRegcaTests.cpp b/tests/UnitTests/PhasorDynamics/runConverterRegcaTests.cpp new file mode 100644 index 000000000..ae8462396 --- /dev/null +++ b/tests/UnitTests/PhasorDynamics/runConverterRegcaTests.cpp @@ -0,0 +1,16 @@ +#include "ConverterRegcaTests.hpp" + +int main() +{ + GridKit::Testing::TestingResults result; + + GridKit::Testing::ConverterRegcaTests test; + + result += test.constructor(); + result += test.lifecycle(); + result += test.signalVerification(); + result += test.nullBusVerification(); + result += test.jsonParseAndSystemAssembly(); + + return result.summary(); +} diff --git a/tests/UnitTests/PhasorDynamics/runSystemSingleComponentTests.cpp b/tests/UnitTests/PhasorDynamics/runSystemSingleComponentTests.cpp index fc4f07e80..2eb8cecb8 100644 --- a/tests/UnitTests/PhasorDynamics/runSystemSingleComponentTests.cpp +++ b/tests/UnitTests/PhasorDynamics/runSystemSingleComponentTests.cpp @@ -14,6 +14,7 @@ int main() result += test.ieeet1(); result += test.load(); result += test.loadZIP(); + result += test.regca(); result += test.genrou(); result += test.genClassical(); result += test.tgov1(); From a04b0dd78ef13321b1ea9a7a6b78cab14fd6f65a Mon Sep 17 00:00:00 2001 From: lukelowry Date: Sun, 17 May 2026 17:21:56 -0500 Subject: [PATCH 2/3] REGCA Implementation and Tests --- GridKit/CommonMath.hpp | 30 +- .../PhasorDynamics/Converter/REECA/README.md | 8 +- .../PhasorDynamics/Converter/REGCA/README.md | 192 +++--- .../PhasorDynamics/Converter/REGCA/Regca.hpp | 31 +- .../Converter/REGCA/RegcaData.hpp | 16 +- .../Converter/REGCA/RegcaEnzyme.cpp | 79 +++ .../Converter/REGCA/RegcaImpl.hpp | 373 +++++++++-- GridKit/Model/PhasorDynamics/INPUT_FORMAT.md | 34 +- .../Math/SmoothnessIndicatorTests.hpp | 50 +- .../Math/runSmoothnessIndicatorTests.cpp | 3 +- .../PhasorDynamics/ConverterRegcaTests.hpp | 589 +++++++++++++++++- .../SystemSingleComponentTests.hpp | 31 +- .../PhasorDynamics/runConverterRegcaTests.cpp | 10 + 13 files changed, 1225 insertions(+), 221 deletions(-) diff --git a/GridKit/CommonMath.hpp b/GridKit/CommonMath.hpp index e4a7ec87f..8f9f8c5cf 100644 --- a/GridKit/CommonMath.hpp +++ b/GridKit/CommonMath.hpp @@ -162,7 +162,31 @@ namespace GridKit } /** - * @brief Smooth two-sided deadband function + * @brief Smooth Type 1 no-offset two-sided deadband function + * + * Smooth approximation to a deadband that returns zero inside the band and + * passes the input through unchanged outside the band. + * + * @tparam ScalarT - scalar data type + * @tparam RealT - Real data type (see GridKit::ScalarTraits::RealT) + * + * @param[in] x - Input signal + * @param[in] lower - Lower breakpoint + * @param[in] upper - Upper breakpoint + * @return Smooth no-offset deadbanded value + */ + template + __attribute__((always_inline)) inline ScalarT deadband1( + const ScalarT x, + const RealT lower, + const RealT upper) + { + assert(lower <= upper); + return x * (sigmoid(lower - x) + sigmoid(x - upper)); + } + + /** + * @brief Smooth Type 2 offset two-sided deadband function * * Smooth approximation to x - min(max(x, lower), upper), composed from the * smooth ramp function. @@ -173,10 +197,10 @@ namespace GridKit * @param[in] x - Input signal * @param[in] lower - Lower breakpoint * @param[in] upper - Upper breakpoint - * @return Smooth deadbanded value + * @return Smooth offset deadbanded value */ template - __attribute__((always_inline)) inline ScalarT deadband( + __attribute__((always_inline)) inline ScalarT deadband2( const ScalarT x, const RealT lower, const RealT upper) diff --git a/GridKit/Model/PhasorDynamics/Converter/REECA/README.md b/GridKit/Model/PhasorDynamics/Converter/REECA/README.md index 6e0d0b83e..834533492 100644 --- a/GridKit/Model/PhasorDynamics/Converter/REECA/README.md +++ b/GridKit/Model/PhasorDynamics/Converter/REECA/README.md @@ -161,7 +161,7 @@ Symbol | Units | Description | $V_T$ | [p.u.] | Terminal voltage magnitude | $V_{\mathrm{meas}}^{\mathrm{safe}}$ | [p.u.] | Safe filtered terminal voltage for divider blocks | Lower bounded by 0.01 $s_{\mathrm{dip}}$ | [binary] | Voltage-dip/overvoltage freeze indicator | 1 when outside voltage thresholds -$V_{\mathrm{err}}$ | [p.u.] | Deadbanded voltage error | Defined by CommonMath `deadband` +$V_{\mathrm{err}}$ | [p.u.] | Deadbanded voltage error | Defined by CommonMath `deadband2` $I_{\mathrm{qv}}$ | [p.u.] | Reactive-current injection candidate | Converter base $Q_{\mathrm{ref}}$ | [p.u.] | Selected reactive-power reference | From power-factor or external reactive-power command $e_Q$ | [p.u.] | Reactive-power control error | Limited $Q_{\mathrm{ref}}$ minus $Q_{\mathrm{gen}}$ @@ -262,7 +262,7 @@ The algebraic targets use CommonMath helper notation where applicable: 0 &= -V_T^2 + V_\mathrm r^2 + V_\mathrm i^2 \\ 0 &= -V_\mathrm{meas}^\mathrm{safe} + \max(V_\mathrm{meas}, 0.01) \\ 0 &= -s_\mathrm{dip} + \text{outside}(V_T, V_\mathrm{dip}, V_\mathrm{up}) \\ - 0 &= -V_\mathrm{err} + \text{deadband}(V_\mathrm{ref0} - V_\mathrm{meas}, D_\mathrm{bd1}, D_\mathrm{bd2}) \\ + 0 &= -V_\mathrm{err} + \text{deadband2}(V_\mathrm{ref0} - V_\mathrm{meas}, D_\mathrm{bd1}, D_\mathrm{bd2}) \\ 0 &= -I_\mathrm{qv} + \text{clamp}(K_\mathrm{qv} V_\mathrm{err}, I_\mathrm{qinj}^{\min}, I_\mathrm{qinj}^{\max}) \\ 0 &= -Q_\mathrm{ref} + s_\mathrm{pf} P_\mathrm{meas}\tan(\phi_\mathrm{pf}^\mathrm{ref}) @@ -290,7 +290,7 @@ The algebraic targets use CommonMath helper notation where applicable: The $V_T$, $I_{\mathrm{q}}^{\mathrm{circ}}$, and $I_{\mathrm{p}}^{\mathrm{circ}}$ variables use nonnegative branches of squared algebraic residuals. This preserves the $s_{PQ}=0$ Q-priority and $s_{PQ}=1$ P-priority current-circle behavior without explicit square roots; a consistent solution should satisfy the nonnegative branch and nonnegative radicands. -CommonMath defines the helper targets and smooth approximations for [min, max, clamp, deadband, and outside](../../../../CommonMath.md#derived-functions). +CommonMath defines the helper targets and smooth approximations for [min, max, clamp, deadband2, and outside](../../../../CommonMath.md#derived-functions). ## Initialization @@ -324,7 +324,7 @@ Then evaluate the upstream algebraic chain: \begin{aligned} V_{\mathrm{meas},0}^{\mathrm{safe}} &= \text{max}(V_{\mathrm{meas},0}, 0.01) \\ s_{\mathrm{dip},0} &= \text{outside}(V_{T,0}, V_{\mathrm{dip}}, V_{\mathrm{up}}) \\ - V_{\mathrm{err},0} &= \text{deadband}(V_{\mathrm{ref0}} - V_{\mathrm{meas},0}, D_{\mathrm{bd1}}, D_{\mathrm{bd2}}) \\ + V_{\mathrm{err},0} &= \text{deadband2}(V_{\mathrm{ref0}} - V_{\mathrm{meas},0}, D_{\mathrm{bd1}}, D_{\mathrm{bd2}}) \\ I_{\mathrm{qv},0} &= \text{clamp}(K_{\mathrm{qv}} V_{\mathrm{err},0}, I_{\mathrm{qinj}}^{\min}, I_{\mathrm{qinj}}^{\max}) \\ Q_{\mathrm{ref},0} &= s_{\mathrm{pf}} P_{\mathrm{meas},0}\tan(\phi_{\mathrm{pf},0}^{\mathrm{ref}}) + s_{\mathrm{pf}}^{\mathrm{off}} Q_{\mathrm{ext},0} \\ e_{Q,0} &= \text{clamp}(Q_{\mathrm{ref},0}, Q^{\min}, Q^{\max}) - Q_{\mathrm{gen},0} \\ diff --git a/GridKit/Model/PhasorDynamics/Converter/REGCA/README.md b/GridKit/Model/PhasorDynamics/Converter/REGCA/README.md index ea352c520..c2961c13f 100644 --- a/GridKit/Model/PhasorDynamics/Converter/REGCA/README.md +++ b/GridKit/Model/PhasorDynamics/Converter/REGCA/README.md @@ -7,6 +7,7 @@ current source at the network interface. Notes: - LVACM uses the unfiltered terminal voltage $V_T$; LVPL uses the filtered voltage $V_M$. - Internal currents are on converter base; bus injections are converted to system base in the network interface. +- Exported branch power measurements are on system base. - HVRCM is represented by internal algebraic current $I_{\mathrm{q}}^{\mathrm{extra}}$. ## Block Diagram @@ -23,9 +24,9 @@ Standard REGCA converter-interface model. Symbol | Units | Description | Typical Value | Note ---------------------------------|----------|-------------------------------------------------------|---------------|------ -$P_{\mathrm{0}}$ | [p.u.] | Initial active power injection | | On system base -$Q_{\mathrm{0}}$ | [p.u.] | Initial reactive power injection | | On system base -$S^{\mathrm{conv}}$ | [MVA] | Converter/model power base | TBD | +$P_0$ | [p.u.] | Initial active power injection on system base | TBD | JSON key: `P0`; fallback command source +$Q_0$ | [p.u.] | Initial reactive power injection on system base | TBD | JSON key: `Q0`; fallback command source +$S^{\text{base}}$ | [MVA] | REGCA model power base | TBD | JSON key: `mva` $T_{\mathrm{g}}$ | [sec] | Converter current-control lag time constant | TBD | $T_M$ | [sec] | Terminal voltage sensor time constant | TBD | Block name: `Tfltr` $R_{\mathrm{q}}^{\max}$ | [p.u./s] | Reactive-current recovery positive rate limit | TBD | Block name: `Iqrmax` @@ -45,7 +46,7 @@ Implementations should reject or report invalid parameter sets: ```math \begin{aligned} - S^{\mathrm{conv}} &> 0 & + S^{\text{base}} &> 0 & T_{\mathrm{g}} &> 0 & T_M &> 0 \\ R_{\mathrm{p}}^{\max} &> 0 & @@ -94,6 +95,8 @@ $I_L$ | [p.u.] | LVPL upper-limit current curve $I_{\mathrm{r}}$ | [p.u.] | Injected current, real component on network reference frame | Converter base $\ell_{\mathrm{p}}$ | [p.u./s] | Smooth active-current lower rate bound | Equivalent to diagram `Rdown` $u_{\mathrm{p}}$ | [p.u./s] | Smooth active-current upper rate bound | Effective `Rup`; includes LVPL anti-windup when $s_L=1$ +$P_{\mathrm{br}}$ | [p.u.] | Active-power output | System base +$Q_{\mathrm{br}}$ | [p.u.] | Reactive-power output | System base ### External Variables @@ -106,8 +109,12 @@ Symbol | Units | Description --------------------------------|--------|------------------------------------------------------------------|------ $V_{\mathrm{r}}$ | [p.u.] | Terminal voltage, real component on network reference frame | Owned by bus object $V_{\mathrm{i}}$ | [p.u.] | Terminal voltage, imaginary component on network reference frame | Owned by bus object -$I_{\mathrm{q}}^{\mathrm{cmd}}$ | [p.u.] | Reactive-current command | Converter base; owned by REEC, constant if no REEC is connected -$I_{\mathrm{p}}^{\mathrm{cmd}}$ | [p.u.] | Active-current command | Converter base; owned by REEC, constant if no REEC is connected +$I_{\mathrm{p}}^{\mathrm{cmd}}$ | [p.u.] | Active-current command in the terminal-voltage reference frame | Converter base; owned by REEC, constant if no REEC is connected +$I_{\mathrm{q}}^{\mathrm{cmd}}$ | [p.u.] | Reactive-current command in the terminal-voltage reference frame | Converter base; owned by REEC, constant if no REEC is connected + +When either `ipcmd` or `iqcmd` port is omitted, REGCA derives the missing +constant command from the optional `P0`/`Q0` pair. If `P0`/`Q0` are omitted, +the fallback commands are zero. ## Model Equations @@ -122,54 +129,39 @@ Define the pre-limit current derivatives: ### Differential Equations -The exact state equations are +The state equations use CommonMath helper notation: ```math \begin{aligned} 0 &= -T_M \dot V_M - V_M + V_T \\ 0 &= -\dot I_{\mathrm{q}} + \begin{cases} - \min(f_{\mathrm{q}}, R_{\mathrm{q}}^{\max}) & Q_{\mathrm{0}} > 0 \\ - \max(f_{\mathrm{q}}, R_{\mathrm{q}}^{\min}) & Q_{\mathrm{0}} \le 0 + \min(f_{\mathrm{q}}, R_{\mathrm{q}}^{\max}) & I_{\mathrm{q},0}^{\mathrm{cmd}} > 0 \\ + \max(f_{\mathrm{q}}, R_{\mathrm{q}}^{\min}) & I_{\mathrm{q},0}^{\mathrm{cmd}} \le 0 \end{cases} \\ 0 &= -\dot I_{\mathrm{p}} + \text{clamp}(f_{\mathrm{p}}, \ell_{\mathrm{p}}, u_{\mathrm{p}}) \end{aligned} ``` -The implemented smooth state equations are - -```math -\begin{aligned} - 0 &= -T_M \dot V_M - V_M + V_T \\ - 0 &= -\dot I_{\mathrm{q}} + - \begin{cases} - f_{\mathrm{q}} - \rho(f_{\mathrm{q}} - R_{\mathrm{q}}^{\max}) - & Q_{\mathrm{0}} > 0 \\ - f_{\mathrm{q}} + \rho(R_{\mathrm{q}}^{\min} - f_{\mathrm{q}}) - & Q_{\mathrm{0}} \le 0 - \end{cases} \\ - 0 &= -\dot I_{\mathrm{p}} + - \ell_{\mathrm{p}} - + \rho(f_{\mathrm{p}} - \ell_{\mathrm{p}}) - - \rho(f_{\mathrm{p}} - u_{\mathrm{p}}) -\end{aligned} -``` - -Here $\rho$ is GridKit's smooth ramp function. The $I_{\mathrm{q}}$ branch is -selected by initial reactive power $Q_{\mathrm{0}}$. The $I_{\mathrm{p}}$ -equation is the smooth clamp of $f_{\mathrm{p}}$ between the algebraic bounds -$\ell_{\mathrm{p}}$ and $u_{\mathrm{p}}$. +The $I_{\mathrm{q}}$ branch is selected by the initial reactive-current +command. CommonMath defines the helper targets and smooth approximations for +[min, max, and clamp](../../../../CommonMath.md#derived-functions). ### Algebraic Equations -The piecewise definitions in this section switch on continuous states, unlike -the $I_{\mathrm{q}}$ differential branch selected by initial conditions. The -exact algebraic targets are: +The exact algebraic targets are: ```math \begin{aligned} 0 &= -V_T^2 + V_{\mathrm{r}}^2 + V_{\mathrm{i}}^2 \\ - I_{\mathrm{i}} &= -I_{\mathrm{q}} + I_{\mathrm{q}}^{\mathrm{extra}} \\ + 0 &= -V_T I_{\mathrm{i}} + + I_{\mathrm{p}}V_{\mathrm{i}} + \begin{cases} + 0 & V_T \le V_{A0} \\ + \dfrac{V_T - V_{A0}}{V_{A1} - V_{A0}} & V_{A0} < V_T < V_{A1} \\ + 1 & V_T \ge V_{A1} + \end{cases} + - (I_{\mathrm{q}} - I_{\mathrm{q}}^{\mathrm{extra}})V_{\mathrm{r}} \\ 0 &= \begin{cases} I_{\mathrm{q}}^{\mathrm{extra}} @@ -177,45 +169,59 @@ exact algebraic targets are: V_T - V_{\mathrm{hv}}^{\max} & I_{\mathrm{q}}^{\mathrm{extra}} > 0 \end{cases} \\ - I_L &= I_{L1} - \begin{cases} - 0 & V_M \le V_{L0} \\ - \dfrac{V_M - V_{L0}}{V_{L1} - V_{L0}} & V_{L0} < V_M < V_{L1} \\ - 1 & V_M \ge V_{L1} - \end{cases} \\ - I_{\mathrm{r}} &= I_{\mathrm{p}} - \begin{cases} - 0 & V_T \le V_{A0} \\ - \dfrac{V_T - V_{A0}}{V_{A1} - V_{A0}} & V_{A0} < V_T < V_{A1} \\ - 1 & V_T \ge V_{A1} - \end{cases} \\ - \ell_{\mathrm{p}} &= + 0 &= -I_L + + I_{L1} + \begin{cases} + 0 & V_M \le V_{L0} \\ + \dfrac{V_M - V_{L0}}{V_{L1} - V_{L0}} & V_{L0} < V_M < V_{L1} \\ + 1 & V_M \ge V_{L1} + \end{cases} \\ + 0 &= -V_T I_{\mathrm{r}} + + I_{\mathrm{p}}V_{\mathrm{r}} + \begin{cases} + 0 & V_T \le V_{A0} \\ + \dfrac{V_T - V_{A0}}{V_{A1} - V_{A0}} & V_{A0} < V_T < V_{A1} \\ + 1 & V_T \ge V_{A1} + \end{cases} + + (I_{\mathrm{q}} - I_{\mathrm{q}}^{\mathrm{extra}})V_{\mathrm{i}} \\ + 0 &= -\ell_{\mathrm{p}} + \begin{cases} -R_{\mathrm{p}}^{\max} & I_{\mathrm{p}} \le 0 \\ -\infty & I_{\mathrm{p}} > 0 \end{cases} \\ - u_{\mathrm{p}} &= + 0 &= -u_{\mathrm{p}} + \begin{cases} R_{\mathrm{p}}^{\max} & I_{\mathrm{p}} \ge 0 \ \land\ (s_L = 0 \lor I_{\mathrm{p}} < I_L) \\ 0 & s_L = 1 \ \land\ I_{\mathrm{p}} \ge I_L \\ \infty & I_{\mathrm{p}} < 0 - \end{cases} + \end{cases} \\ + 0 &= -P_{\mathrm{br}} + + \dfrac{S^{\text{base}}}{S^{\text{sys}}} + \left(V_{\mathrm{r}} I_{\mathrm{r}} + V_{\mathrm{i}} I_{\mathrm{i}}\right) \\ + 0 &= -Q_{\mathrm{br}} + + \dfrac{S^{\text{base}}}{S^{\text{sys}}} + \left(V_{\mathrm{i}} I_{\mathrm{r}} - V_{\mathrm{r}} I_{\mathrm{i}}\right) \end{aligned} ``` -The implemented algebraic residuals use smooth $\text{linseg}$, -$\rho$, and $\sigma$ operators: +GridKit implements the algebraic residuals with smooth $\text{linseg}$, +$\rho$, and $\sigma$ operators; CommonMath defines the +[linear segment, ramp, and sigmoid helpers](../../../../CommonMath.md#derived-functions). +An equivalent residual form is: ```math \begin{aligned} 0 &= -V_T^2 + V_{\mathrm{r}}^2 + V_{\mathrm{i}}^2 \\ - 0 &= -I_{\mathrm{i}} - I_{\mathrm{q}} + I_{\mathrm{q}}^{\mathrm{extra}} \\ + 0 &= -V_T I_{\mathrm{i}} + + I_{\mathrm{p}}\text{linseg}(V_T;\ V_{A0},\ V_{A1},\ 1)V_{\mathrm{i}} + - (I_{\mathrm{q}} - I_{\mathrm{q}}^{\mathrm{extra}})V_{\mathrm{r}} \\ 0 &= -I_{\mathrm{q}}^{\mathrm{extra}} + \rho\!\left(I_{\mathrm{q}}^{\mathrm{extra}} - (V_{\mathrm{hv}}^{\max} - V_T)\right) \\ 0 &= -I_L + \text{linseg}(V_M;\ V_{L0},\ V_{L1},\ I_{L1}) \\ - 0 &= -I_{\mathrm{r}} - + I_{\mathrm{p}}\text{linseg}(V_T;\ V_{A0},\ V_{A1},\ 1) \\ + 0 &= -V_T I_{\mathrm{r}} + + I_{\mathrm{p}}\text{linseg}(V_T;\ V_{A0},\ V_{A1},\ 1)V_{\mathrm{r}} + + (I_{\mathrm{q}} - I_{\mathrm{q}}^{\mathrm{extra}})V_{\mathrm{i}} \\ 0 &= -\ell_{\mathrm{p}} - R_{\mathrm{p}}^{\max} - (M_{\mathrm{p}} - R_{\mathrm{p}}^{\max})\sigma(I_{\mathrm{p}}) \\ @@ -228,7 +234,13 @@ $\rho$, and $\sigma$ operators: M_{\mathrm{p}}(1-\sigma(I_{\mathrm{p}})) + R_{\mathrm{p}}^{\max}\sigma(I_{\mathrm{p}}) & s_L = 0 - \end{cases} + \end{cases} \\ + 0 &= -P_{\mathrm{br}} + + \dfrac{S^{\text{base}}}{S^{\text{sys}}} + \left(V_{\mathrm{r}} I_{\mathrm{r}} + V_{\mathrm{i}} I_{\mathrm{i}}\right) \\ + 0 &= -Q_{\mathrm{br}} + + \dfrac{S^{\text{base}}}{S^{\text{sys}}} + \left(V_{\mathrm{i}} I_{\mathrm{r}} - V_{\mathrm{r}} I_{\mathrm{i}}\right) \end{aligned} ``` @@ -241,8 +253,8 @@ REGCA currents: ```math \begin{aligned} - I_{\mathrm{r}}^{\mathrm{inj}} &:= I_{\mathrm{r}}\dfrac{S^{\mathrm{conv}}}{S^{\mathrm{sys}}} \\ - I_{\mathrm{i}}^{\mathrm{inj}} &:= I_{\mathrm{i}}\dfrac{S^{\mathrm{conv}}}{S^{\mathrm{sys}}} + I_{\mathrm{r}}^{\mathrm{inj}} &:= I_{\mathrm{r}}\dfrac{S^{\text{base}}}{S^{\text{sys}}} \\ + I_{\mathrm{i}}^{\mathrm{inj}} &:= I_{\mathrm{i}}\dfrac{S^{\text{base}}}{S^{\text{sys}}} \end{aligned} ``` @@ -256,52 +268,66 @@ steady-state initial values: ```math \begin{aligned} V_T &= \sqrt{V_\mathrm{r}^2 + V_\mathrm{i}^2} \\ - I_\mathrm{r0} &= \dfrac{P_0 V_\mathrm{r} + Q_0 V_\mathrm{i}}{V_T^2} - \dfrac{S^\mathrm{sys}}{S^\mathrm{conv}} \\ - I_\mathrm{i0} &= \dfrac{P_0 V_\mathrm{i} - Q_0 V_\mathrm{r}}{V_T^2} - \dfrac{S^\mathrm{sys}}{S^\mathrm{conv}} \\ V_{M0} &= V_T \\ I_{L0} &= \text{linseg}(V_T;\ V_{L0},\ V_{L1},\ I_{L1}) \\ - I_\mathrm{p0} &= \dfrac{I_\mathrm{r0}} - {\text{linseg}(V_T;\ V_{A0},\ V_{A1},\ 1)} \\ - \ell_\mathrm{p0} &= -R_\mathrm{p}^{\max} - - (M_\mathrm{p} - R_\mathrm{p}^{\max})\sigma(I_\mathrm{p0}) \\ - u_\mathrm{p0} &= + I_\mathrm{p0} &= + \begin{cases} + I_\mathrm{p}^\mathrm{cmd}(0) & \text{if the } I_\mathrm{p}^\mathrm{cmd} \text{ port is attached} \\ + \dfrac{P_0} + {V_T\text{linseg}(V_T;\ V_{A0},\ V_{A1},\ 1)} + \dfrac{S^{\text{sys}}}{S^{\text{base}}} + & \text{otherwise} + \end{cases} \\ + I_\mathrm{q0} &= \begin{cases} - M_\mathrm{p}(1-\sigma(I_\mathrm{p0})) - + R_\mathrm{p}^{\max}\sigma(I_\mathrm{p0}) - \sigma(I_{L0} - I_\mathrm{p0}) - & s_L = 1 \\ - M_\mathrm{p}(1-\sigma(I_\mathrm{p0})) - + R_\mathrm{p}^{\max}\sigma(I_\mathrm{p0}) - & s_L = 0 + I_\mathrm{q}^\mathrm{cmd}(0) & \text{if the } I_\mathrm{q}^\mathrm{cmd} \text{ port is attached} \\ + \dfrac{Q_0}{V_T}\dfrac{S^{\text{sys}}}{S^{\text{base}}} + & \text{otherwise} \end{cases} \\ + \ell_\mathrm{p0} &= -R_\mathrm{p}^{\max} + - (M_\mathrm{p} - R_\mathrm{p}^{\max})\sigma(I_\mathrm{p0}) \\ + u_\mathrm{p0} &= M_\mathrm{p}(1-\sigma(I_\mathrm{p0})) + + R_\mathrm{p}^{\max}\sigma(I_\mathrm{p0}) + \left(1 - s_L + s_L\sigma(I_{L0} - I_\mathrm{p0})\right) \\ I_\mathrm{q0}^\mathrm{extra} &= 0 \\ - I_\mathrm{q0} &= -I_\mathrm{i0} + I_\mathrm{q0}^\mathrm{extra} \\ - I_\mathrm{p0}^\mathrm{cmd} &= I_\mathrm{p0} \\ - I_\mathrm{q0}^\mathrm{cmd} &= I_\mathrm{q0} + I_\mathrm{i0} &= \dfrac{ + \text{linseg}(V_T;\ V_{A0},\ V_{A1},\ 1)I_{\mathrm{p0}}V_{\mathrm{i}} + - I_{\mathrm{q0}}V_{\mathrm{r}} + }{V_T} \\ + I_\mathrm{r0} &= \dfrac{ + \text{linseg}(V_T;\ V_{A0},\ V_{A1},\ 1)I_{\mathrm{p0}}V_{\mathrm{r}} + + I_{\mathrm{q0}}V_{\mathrm{i}} + }{V_T} \\ + P_{\mathrm{br},0} &= \dfrac{S^{\text{base}}}{S^{\text{sys}}} + \left(V_{\mathrm{r}} I_{\mathrm{r0}} + V_{\mathrm{i}} I_{\mathrm{i0}}\right) \\ + Q_{\mathrm{br},0} &= \dfrac{S^{\text{base}}}{S^{\text{sys}}} + \left(V_{\mathrm{i}} I_{\mathrm{r0}} - V_{\mathrm{r}} I_{\mathrm{i0}}\right) \end{aligned} ``` For normal power-flow starts, $V_T > V_{A1}$, so $\text{linseg}(V_T;\ V_{A0},\ V_{A1},\ 1) = 1$ and the $I_{\mathrm{p0}}$ formula is well defined. +When a command port is omitted, GridKit stores the resolved initial value as +the constant command used during residual evaluation. Initialization should verify: -- $V_T \le V_{\mathrm{hv}}^{\max}$. If $V_T \ge V_{\mathrm{hv}}^{\max}$, +- $V_T < V_{\mathrm{hv}}^{\max}$. If $V_T \ge V_{\mathrm{hv}}^{\max}$, $I_{\mathrm{q0}}^{\mathrm{extra}} = 0$ may not satisfy the HVRCM algebraic condition, and a nonzero value should be solved or the initialization rejected. - $\text{linseg}(V_T;\ V_{A0},\ V_{A1},\ 1) > 0$ when - $I_{\mathrm{r0}} \ne 0$. If the LVACM gain is zero, no finite - $I_{\mathrm{p0}}$ can reproduce nonzero initial active current. + the `ipcmd` port is omitted and $P_0 \ne 0$. If the LVACM gain is zero, + no finite $I_{\mathrm{p0}}$ can reproduce nonzero initial active power. All internal derivatives initialize to zero. ## Model Outputs Real and imaginary injected currents, $I_{\mathrm{r}}$ and -$I_{\mathrm{i}}$, are converter-base algebraic variables. System-base power -outputs use the bus-facing currents: +$I_{\mathrm{i}}$, are converter-base algebraic variables. Optional +`ibranchr`, `ibranchi`, `pbranch`, and `qbranch` ports export REGCA-owned +algebraic measurement states for downstream controllers such as REPCA. The +`pbranch` and `qbranch` signal ports export system-base power outputs: ```math \begin{aligned} P &= V_{\mathrm{r}} I_{\mathrm{r}}^{\mathrm{inj}} + V_{\mathrm{i}} I_{\mathrm{i}}^{\mathrm{inj}} \\ diff --git a/GridKit/Model/PhasorDynamics/Converter/REGCA/Regca.hpp b/GridKit/Model/PhasorDynamics/Converter/REGCA/Regca.hpp index cbafdc7d8..db10ae6c9 100644 --- a/GridKit/Model/PhasorDynamics/Converter/REGCA/Regca.hpp +++ b/GridKit/Model/PhasorDynamics/Converter/REGCA/Regca.hpp @@ -46,6 +46,8 @@ namespace GridKit IR, ///< Real injected current LP, ///< Active-current lower rate bound UP, ///< Active-current upper rate bound + PBR, ///< Active-power output on system base + QBR, ///< Reactive-power output on system base MAXIMUM, }; @@ -68,12 +70,12 @@ namespace GridKit using Component::J_cols_buffer_; using Component::J_rows_buffer_; using Component::J_vals_buffer_; - using Component::mva_system_base_; using Component::nnz_; using Component::residual_indices_; using Component::size_; using Component::tag_; using Component::time_; + using Component::va_system_base_; using Component::variable_indices_; using Component::wb_; using Component::y_; @@ -117,6 +119,20 @@ namespace GridKit private: void initializeParameters(const model_data_type& data); void initializeMonitor(); + void setDerivedParameters(); + + ScalarT toComponentBase(ScalarT value) const + { + return value * va_system_base_ / va_converter_base_; + } + + ScalarT toSystemBase(ScalarT value) const + { + return value / toComponentBase(static_cast(ONE)); + } + + ScalarT activeCurrentLowerRateBound(ScalarT ip) const; + ScalarT activeCurrentUpperRateBound(ScalarT ip, ScalarT il) const; ScalarT& Vr() { @@ -142,7 +158,7 @@ namespace GridKit RealT P0_{0}; RealT Q0_{0}; - RealT Sconv_{0}; + RealT mva_base_{0}; RealT Tg_{0}; RealT TM_{0}; RealT Rqmax_{0}; @@ -157,6 +173,17 @@ namespace GridKit RealT Vhvmax_{0}; IdxT bus_id_{0}; + IdxT parameter_error_count_{0}; + RealT Mp_{0}; + RealT va_converter_base_{0}; + RealT use_lvpl_{0}; + RealT bypass_lvpl_{1}; + RealT iq_use_upper_{0}; + RealT iq_use_lower_{1}; + + ScalarT ipcmd_set_{0}; + ScalarT iqcmd_set_{0}; + ComponentSignals signals_; std::unique_ptr monitor_; diff --git a/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaData.hpp b/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaData.hpp index 29f8530fe..49d23b861 100644 --- a/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaData.hpp +++ b/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaData.hpp @@ -17,9 +17,9 @@ namespace GridKit /// Parameter keys for the REGCA converter model. enum class RegcaParameters { - P0, ///< Initial active power injection on system base - Q0, ///< Initial reactive power injection on system base - Sconv, ///< Converter/model power base + P0, ///< Initial active power on system base + Q0, ///< Initial reactive power on system base + mva, ///< MVA base of the REGCA model Tg, ///< Converter current-control lag time constant TM, ///< Terminal voltage sensor time constant Rqmax, ///< Reactive-current recovery positive rate limit @@ -37,9 +37,13 @@ namespace GridKit /// Ports for the REGCA converter model. enum class RegcaPorts { - bus, ///< Terminal bus ID - ipcmd, ///< Optional active-current command signal ID - iqcmd ///< Optional reactive-current command signal ID + bus, ///< Terminal bus ID + ipcmd, ///< Optional active-current command signal ID + iqcmd, ///< Optional reactive-current command signal ID + ibranchr, ///< Optional real current measurement output signal ID + ibranchi, ///< Optional imaginary current measurement output signal ID + pbranch, ///< Optional active-power measurement output signal ID + qbranch ///< Optional reactive-power measurement output signal ID }; /// Variables available through the monitor interface. diff --git a/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaEnzyme.cpp b/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaEnzyme.cpp index 7ad8ee7fc..8bfb167d7 100644 --- a/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaEnzyme.cpp +++ b/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaEnzyme.cpp @@ -21,6 +21,85 @@ namespace GridKit Log::misc() << "Jacobian evaluation is experimental!" << std::endl; J_.zeroMatrix(); + if (J_rows_buffer_ == nullptr) + { + // Reserve space for a dense internal block. Enzyme's sparse storage + // keeps only structural nonzeros for each differentiated block. + J_rows_buffer_ = new IdxT[static_cast(size_) * static_cast(size_)]; + J_cols_buffer_ = new IdxT[static_cast(size_) * static_cast(size_)]; + J_vals_buffer_ = new RealT[static_cast(size_) * static_cast(size_)]; + } + + using RegcaT = GridKit::PhasorDynamics::Converter::Regca; + using Fn = GridKit::Enzyme::Sparse::MemberFunctions; + + GridKit::Enzyme::Sparse::DfDy::eval(this, + f_.size(), + y_.size(), + (this->getResidualIndices()).data(), + (this->getVariableIndices()).data(), + y_.data(), + yp_.data(), + wb_.data(), + ws_.data(), + alpha_, + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + J_); + + GridKit::Enzyme::Sparse::DfDwb::eval(this, + f_.size(), + static_cast(bus_->size()), + (this->getResidualIndices()).data(), + (bus_->getVariableIndices()).data(), + y_.data(), + yp_.data(), + wb_.data(), + ws_.data(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + J_); + + GridKit::Enzyme::Sparse::DfDws::eval(this, + f_.size(), + ws_.size(), + (this->getResidualIndices()).data(), + ws_indices_.data(), + y_.data(), + yp_.data(), + wb_.data(), + ws_.data(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + J_); + + GridKit::Enzyme::Sparse::DhDy::eval(this, + static_cast(bus_->size()), + y_.size(), + (bus_->getResidualIndices()).data(), + (this->getVariableIndices()).data(), + y_.data(), + yp_.data(), + wb_.data(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + J_); return 0; } diff --git a/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaImpl.hpp b/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaImpl.hpp index 6598ea24e..212d29788 100644 --- a/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaImpl.hpp @@ -7,6 +7,7 @@ #pragma once #include +#include #include #include @@ -46,16 +47,50 @@ namespace GridKit { } + template + void Regca::setDerivedParameters() + { + Mp_ = static_cast(100.0) * Rpmax_; + use_lvpl_ = ZERO; + bypass_lvpl_ = ONE; + if (sL_) + { + use_lvpl_ = ONE; + bypass_lvpl_ = ZERO; + } + iq_use_upper_ = ZERO; + iq_use_lower_ = ONE; + va_converter_base_ = mva_base_ * static_cast(1.0e6); + } + + template + ScalarT Regca::activeCurrentLowerRateBound(ScalarT ip) const + { + return -Rpmax_ - (Mp_ - Rpmax_) * Math::sigmoid(ip); + } + + template + ScalarT Regca::activeCurrentUpperRateBound(ScalarT ip, ScalarT il) const + { + const ScalarT sigma_ip = Math::sigmoid(ip); + return Mp_ * (ONE - sigma_ip) + + Rpmax_ * sigma_ip * (bypass_lvpl_ + use_lvpl_ * Math::sigmoid(il - ip)); + } + template void Regca::initializeParameters(const model_data_type& data) { using Params = typename model_data_type::Parameters; using Ports = typename model_data_type::Ports; - auto loadReal = [&](auto key, RealT& target) + parameter_error_count_ = 0; + + auto loadRequiredReal = [&](auto key, RealT& target, const char* name) { if (!data.parameters.contains(key)) { + Log::error() << "Regca: missing required parameter '" << name << "'\n"; + ++parameter_error_count_; return; } @@ -68,46 +103,84 @@ namespace GridKit { target = static_cast(*index_value); } + else + { + Log::error() << "Regca: parameter '" << name << "' must be numeric\n"; + ++parameter_error_count_; + } }; - auto loadBool = [&](auto key, bool& target) + auto loadOptionalReal = [&](auto key, RealT& target, const char* name) { if (!data.parameters.contains(key)) { return; } + const auto& value = data.parameters.at(key); + if (const auto* real_value = std::get_if(&value)) + { + target = *real_value; + } + else if (const auto* index_value = std::get_if(&value)) + { + target = static_cast(*index_value); + } + else + { + Log::error() << "Regca: parameter '" << name << "' must be numeric\n"; + ++parameter_error_count_; + } + }; + + auto loadRequiredSwitch = [&](auto key, bool& target, const char* name) + { + if (!data.parameters.contains(key)) + { + Log::error() << "Regca: missing required parameter '" << name << "'\n"; + ++parameter_error_count_; + return; + } + const auto& value = data.parameters.at(key); if (const auto* bool_value = std::get_if(&value)) { target = *bool_value; } - else if (const auto* index_value = std::get_if(&value)) + else if (const auto* index_value = std::get_if(&value); + index_value && (*index_value == 0 || *index_value == 1)) { - target = (*index_value != 0); + target = (*index_value == 1); + } + else + { + Log::error() << "Regca: parameter '" << name << "' must be bool or 0/1\n"; + ++parameter_error_count_; } }; - loadReal(Params::P0, P0_); - loadReal(Params::Q0, Q0_); - loadReal(Params::Sconv, Sconv_); - loadReal(Params::Tg, Tg_); - loadReal(Params::TM, TM_); - loadReal(Params::Rqmax, Rqmax_); - loadReal(Params::Rqmin, Rqmin_); - loadReal(Params::Rpmax, Rpmax_); - loadBool(Params::sL, sL_); - loadReal(Params::IL1, IL1_); - loadReal(Params::VL0, VL0_); - loadReal(Params::VL1, VL1_); - loadReal(Params::VA0, VA0_); - loadReal(Params::VA1, VA1_); - loadReal(Params::Vhvmax, Vhvmax_); + loadOptionalReal(Params::P0, P0_, "P0"); + loadOptionalReal(Params::Q0, Q0_, "Q0"); + loadRequiredReal(Params::mva, mva_base_, "mva"); + loadRequiredReal(Params::Tg, Tg_, "Tg"); + loadRequiredReal(Params::TM, TM_, "TM"); + loadRequiredReal(Params::Rqmax, Rqmax_, "Rqmax"); + loadRequiredReal(Params::Rqmin, Rqmin_, "Rqmin"); + loadRequiredReal(Params::Rpmax, Rpmax_, "Rpmax"); + loadRequiredSwitch(Params::sL, sL_, "sL"); + loadRequiredReal(Params::IL1, IL1_, "IL1"); + loadRequiredReal(Params::VL0, VL0_, "VL0"); + loadRequiredReal(Params::VL1, VL1_, "VL1"); + loadRequiredReal(Params::VA0, VA0_, "VA0"); + loadRequiredReal(Params::VA1, VA1_, "VA1"); + loadRequiredReal(Params::Vhvmax, Vhvmax_, "Vhvmax"); if (data.ports.contains(Ports::bus)) { bus_id_ = data.ports.at(Ports::bus); } + + setDerivedParameters(); } template @@ -129,10 +202,10 @@ namespace GridKit { return y_[index(RegcaInternalVariables::IR)]; }); monitor_->set(Variable::ii, [this, index] { return y_[index(RegcaInternalVariables::II)]; }); - monitor_->set(Variable::p, [] - { return ScalarT{0}; }); - monitor_->set(Variable::q, [] - { return ScalarT{0}; }); + monitor_->set(Variable::p, [this, index] + { return y_[index(RegcaInternalVariables::PBR)]; }); + monitor_->set(Variable::q, [this, index] + { return y_[index(RegcaInternalVariables::QBR)]; }); monitor_->set(Variable::vt, [this, index] { return y_[index(RegcaInternalVariables::VT)]; }); monitor_->set(Variable::vm, [this, index] @@ -184,13 +257,50 @@ namespace GridKit this->setResidualIndex(j, j); } + if (signals_.template isAssigned()) + { + signals_.template getSignalNode()->set( + &y_[static_cast(RegcaInternalVariables::IR)], + &(this->getVariableIndex(static_cast(RegcaInternalVariables::IR)))); + } + + if (signals_.template isAssigned()) + { + signals_.template getSignalNode()->set( + &y_[static_cast(RegcaInternalVariables::II)], + &(this->getVariableIndex(static_cast(RegcaInternalVariables::II)))); + } + + if (signals_.template isAssigned()) + { + signals_.template getSignalNode()->set( + &y_[static_cast(RegcaInternalVariables::PBR)], + &(this->getVariableIndex(static_cast(RegcaInternalVariables::PBR)))); + } + + if (signals_.template isAssigned()) + { + signals_.template getSignalNode()->set( + &y_[static_cast(RegcaInternalVariables::QBR)], + &(this->getVariableIndex(static_cast(RegcaInternalVariables::QBR)))); + } + return 0; } template int Regca::verify() const { - int ret = 0; + int ret = static_cast(parameter_error_count_); + + auto check = [&](bool condition, const char* message) + { + if (!condition) + { + Log::error() << "Regca: " << message << '\n'; + ret += 1; + } + }; if (bus_ == nullptr) { @@ -198,6 +308,16 @@ namespace GridKit ret += 1; } + check(mva_base_ > ZERO, "mva must be positive"); + check(Tg_ > ZERO, "Tg must be positive"); + check(TM_ > ZERO, "TM must be positive"); + check(Rpmax_ > ZERO, "Rpmax must be positive"); + check(Rqmin_ < ZERO && ZERO < Rqmax_, "Rqmin < 0 < Rqmax is required"); + check(IL1_ >= ZERO, "IL1 must be non-negative"); + check(ZERO <= VL0_ && VL0_ < VL1_, "VL0/VL1 must satisfy 0 <= VL0 < VL1"); + check(ZERO <= VA0_ && VA0_ < VA1_, "VA0/VA1 must satisfy 0 <= VA0 < VA1"); + check(Vhvmax_ > ZERO, "Vhvmax must be positive"); + if (signals_.template isAttached()) { if (!signals_.template isLinked()) @@ -222,8 +342,107 @@ namespace GridKit template int Regca::initialize() { - std::fill(y_.begin(), y_.end(), ScalarT{0}); - std::fill(yp_.begin(), yp_.end(), ScalarT{0}); + if (bus_ == nullptr) + { + Log::error() << "Regca: cannot initialize with null bus\n"; + return 1; + } + + if (parameter_error_count_ > 0 || mva_base_ <= ZERO || Tg_ <= ZERO + || TM_ <= ZERO || Rpmax_ <= ZERO + || !(Rqmin_ < ZERO && ZERO < Rqmax_) + || IL1_ < ZERO || !(ZERO <= VL0_ && VL0_ < VL1_) + || !(ZERO <= VA0_ && VA0_ < VA1_) || Vhvmax_ <= ZERO) + { + Log::error() << "Regca: cannot initialize with invalid parameters\n"; + return 1; + } + + const auto VM = static_cast(RegcaInternalVariables::VM); + const auto IQ = static_cast(RegcaInternalVariables::IQ); + const auto IP = static_cast(RegcaInternalVariables::IP); + const auto VT = static_cast(RegcaInternalVariables::VT); + const auto II = static_cast(RegcaInternalVariables::II); + const auto IQEXTRA = static_cast(RegcaInternalVariables::IQEXTRA); + const auto IL = static_cast(RegcaInternalVariables::IL); + const auto IR = static_cast(RegcaInternalVariables::IR); + const auto LP = static_cast(RegcaInternalVariables::LP); + const auto UP = static_cast(RegcaInternalVariables::UP); + const auto PBR = static_cast(RegcaInternalVariables::PBR); + const auto QBR = static_cast(RegcaInternalVariables::QBR); + + const ScalarT vr = Vr(); + const ScalarT vi = Vi(); + const ScalarT vt2 = vr * vr + vi * vi; + const ScalarT vt = std::sqrt(vt2); + + if (vt <= ZERO) + { + Log::error() << "Regca: terminal voltage magnitude must be positive at initialization\n"; + return 1; + } + + if (vt >= Vhvmax_) + { + Log::error() << "Regca: terminal voltage magnitude must be less than Vhvmax at initialization\n"; + return 1; + } + + const ScalarT lvacm = Math::linseg(vt, VA0_, VA1_, ONE); + ScalarT ipcmd0{ZERO}; + ScalarT iqcmd0{ZERO}; + + if (signals_.template isAttached()) + { + ipcmd0 = signals_.template readExternalVariable(); + } + else if (P0_ != ZERO) + { + if (vt <= VA0_ || lvacm <= ZERO) + { + Log::error() << "Regca: LVACM gain is zero with nonzero initial active power\n"; + return 1; + } + + ipcmd0 = toComponentBase(static_cast(P0_) / vt) / lvacm; + } + + if (signals_.template isAttached()) + { + iqcmd0 = signals_.template readExternalVariable(); + } + else if (Q0_ != ZERO) + { + iqcmd0 = toComponentBase(static_cast(Q0_) / vt); + } + + const ScalarT iqextra0{ZERO}; + const ScalarT qnet0 = iqcmd0 - iqextra0; + + y_[VM] = vt; + y_[IQ] = iqcmd0; + y_[IP] = ipcmd0; + y_[VT] = vt; + y_[II] = (lvacm * y_[IP] * vi - qnet0 * vr) / vt; + y_[IQEXTRA] = iqextra0; + y_[IL] = Math::linseg(vt, VL0_, VL1_, IL1_); + y_[IR] = (lvacm * y_[IP] * vr + qnet0 * vi) / vt; + y_[LP] = activeCurrentLowerRateBound(y_[IP]); + y_[UP] = activeCurrentUpperRateBound(y_[IP], y_[IL]); + y_[PBR] = toSystemBase(vr * y_[IR] + vi * y_[II]); + y_[QBR] = toSystemBase(vi * y_[IR] - vr * y_[II]); + + ipcmd_set_ = ipcmd0; + iqcmd_set_ = iqcmd0; + iq_use_upper_ = ZERO; + iq_use_lower_ = ONE; + if (static_cast(iqcmd0) > ZERO) + { + iq_use_upper_ = ONE; + iq_use_lower_ = ZERO; + } + + std::fill(yp_.begin(), yp_.end(), ZERO); return 0; } @@ -239,51 +458,115 @@ namespace GridKit template __attribute__((always_inline)) inline int Regca::evaluateInternalResidual( - [[maybe_unused]] ScalarT* y, - [[maybe_unused]] ScalarT* yp, - [[maybe_unused]] ScalarT* wb, - [[maybe_unused]] ScalarT* ws, - ScalarT* f) + ScalarT* y, + ScalarT* yp, + ScalarT* wb, + ScalarT* ws, + ScalarT* f) { - for (IdxT i = 0; i < size_; ++i) - { - f[static_cast(i)] = ScalarT{0}; - } + const auto VM = static_cast(RegcaInternalVariables::VM); + const auto IQ = static_cast(RegcaInternalVariables::IQ); + const auto IP = static_cast(RegcaInternalVariables::IP); + const auto VT = static_cast(RegcaInternalVariables::VT); + const auto II = static_cast(RegcaInternalVariables::II); + const auto IQEXTRA = static_cast(RegcaInternalVariables::IQEXTRA); + const auto IL = static_cast(RegcaInternalVariables::IL); + const auto IR = static_cast(RegcaInternalVariables::IR); + const auto LP = static_cast(RegcaInternalVariables::LP); + const auto UP = static_cast(RegcaInternalVariables::UP); + const auto PBR = static_cast(RegcaInternalVariables::PBR); + const auto QBR = static_cast(RegcaInternalVariables::QBR); + + const auto IPCMD = static_cast(RegcaExternalVariables::IPCMD); + const auto IQCMD = static_cast(RegcaExternalVariables::IQCMD); + + const ScalarT vm = y[VM]; + const ScalarT iq = y[IQ]; + const ScalarT ip = y[IP]; + const ScalarT vt = y[VT]; + const ScalarT ii = y[II]; + const ScalarT iqextra = y[IQEXTRA]; + const ScalarT il = y[IL]; + const ScalarT ir = y[IR]; + const ScalarT lp = y[LP]; + const ScalarT up = y[UP]; + const ScalarT pbr = y[PBR]; + const ScalarT qbr = y[QBR]; + + const ScalarT vm_dot = yp[VM]; + const ScalarT iq_dot = yp[IQ]; + const ScalarT ip_dot = yp[IP]; + + const ScalarT vr = wb[0]; + const ScalarT vi = wb[1]; + + const ScalarT ipcmd = ws[IPCMD]; + const ScalarT iqcmd = ws[IQCMD]; + + const ScalarT fq = (iqcmd - iq) / Tg_; + const ScalarT fp = (ipcmd - ip) / Tg_; + + const ScalarT iq_rate = + iq_use_upper_ * (fq - Math::ramp(fq - Rqmax_)) + + iq_use_lower_ * (fq + Math::ramp(Rqmin_ - fq)); + + const ScalarT ip_rate = lp + Math::ramp(fp - lp) - Math::ramp(fp - up); + + f[VM] = -vm_dot + (vt - vm) / TM_; + f[IQ] = -iq_dot + iq_rate; + f[IP] = -ip_dot + ip_rate; + const ScalarT lvacm = Math::linseg(vt, VA0_, VA1_, ONE); + const ScalarT qnet = iq - iqextra; + + f[VT] = -vt * vt + vr * vr + vi * vi; + f[II] = -vt * ii + lvacm * ip * vi - qnet * vr; + f[IQEXTRA] = -iqextra + Math::ramp(iqextra - (Vhvmax_ - vt)); + f[IL] = -il + Math::linseg(vm, VL0_, VL1_, IL1_); + f[IR] = -vt * ir + lvacm * ip * vr + qnet * vi; + f[LP] = -lp + activeCurrentLowerRateBound(ip); + f[UP] = -up + activeCurrentUpperRateBound(ip, il); + f[PBR] = -pbr + toSystemBase(vr * ir + vi * ii); + f[QBR] = -qbr + toSystemBase(vi * ir - vr * ii); return 0; } template __attribute__((always_inline)) inline int Regca::evaluateBusResidual( - [[maybe_unused]] ScalarT* y, + ScalarT* y, [[maybe_unused]] ScalarT* yp, [[maybe_unused]] ScalarT* wb, ScalarT* h) { - h[0] = ScalarT{0}; - h[1] = ScalarT{0}; + const auto II = static_cast(RegcaInternalVariables::II); + const auto IR = static_cast(RegcaInternalVariables::IR); + + h[0] = toSystemBase(y[IR]); + h[1] = toSystemBase(y[II]); return 0; } template int Regca::evaluateResidual() { - std::fill(ws_.begin(), ws_.end(), ScalarT{0}); + const auto IPCMD = static_cast(RegcaExternalVariables::IPCMD); + const auto IQCMD = static_cast(RegcaExternalVariables::IQCMD); + + ws_[IPCMD] = ipcmd_set_; + ws_[IQCMD] = iqcmd_set_; std::fill(ws_indices_.begin(), ws_indices_.end(), INVALID_INDEX); if (signals_.template isAttached()) { - const auto index = static_cast(RegcaExternalVariables::IPCMD); - ws_[index] = signals_.template readExternalVariable(); - ws_indices_[index] = + ws_[IPCMD] = signals_.template readExternalVariable(); + ws_indices_[IPCMD] = signals_.template readExternalVariableIndex(); } if (signals_.template isAttached()) { - const auto index = static_cast(RegcaExternalVariables::IQCMD); - ws_[index] = signals_.template readExternalVariable(); - ws_indices_[index] = + ws_[IQCMD] = signals_.template readExternalVariable(); + ws_indices_[IQCMD] = signals_.template readExternalVariableIndex(); } diff --git a/GridKit/Model/PhasorDynamics/INPUT_FORMAT.md b/GridKit/Model/PhasorDynamics/INPUT_FORMAT.md index 48db096f8..b33849140 100644 --- a/GridKit/Model/PhasorDynamics/INPUT_FORMAT.md +++ b/GridKit/Model/PhasorDynamics/INPUT_FORMAT.md @@ -144,7 +144,7 @@ are specified: `Genrou` | 6th order machine model | `bus`, `pmech`\*, `speed`\*, `efd`\* | `p0`, `q0`, `H`, `D`, `Ra`, `Tdop`, `Tdopp`, `Tqop`, `Tqopp`, `Xd`, `Xdp`, `Xdpp`, `Xq`, `Xqp`, `Xqpp`, `Xl`, `S10`, `S12`, `mva` | `ir`, `ii`, `p`, `q`, `delta`, `omega`, `speed` `Gensal` | 5th order salient-pole machine model | `bus`, `pmech`\*, `speed`\*, `efd`\* | `p0`, `q0`, `H`, `D`, `Ra`, `Tdop`, `Tdopp`, `Tqopp`, `Xd`, `Xdp`, `Xdpp`, `Xq`, `Xl`, `S10`, `S12`, `mva` | `ir`, `ii`, `p`, `q`, `delta`, `omega`, `speed`, `Eqp`, `psidp`, `psiqpp`, `psidpp`, `vd`, `vq`, `te`, `id`, `iq` `GenClassical` | the classical machine model | `bus`, `pmech`\*, `speed`\*, `efd`\* | `p0`, `q0`, `H`, `D`, `Ra`, `Xdp`, `mva` | `ir`, `ii`, `p`, `q`, `delta`, `omega` - `Regca` | WECC REGCA renewable generator/converter skeleton | `bus`, `ipcmd`\*, `iqcmd`\* | `P0`, `Q0`, `Sconv`, `Tg`, `TM`, `Rqmax`, `Rqmin`, `Rpmax`, `sL`, `IL1`, `VL0`, `VL1`, `VA0`, `VA1`, `Vhvmax` | `ir`, `ii`, `p`, `q`, `vt`, `vm`, `ip`, `iq`, `iqextra`, `il`, `lp`, `up` + `Regca` | WECC REGCA renewable generator/converter model | `bus`, `ipcmd`\*, `iqcmd`\*, `ibranchr`\*, `ibranchi`\*, `pbranch`\*, `qbranch`\* | `P0`, `Q0`, `mva`, `Tg`, `TM`, `Rqmax`, `Rqmin`, `Rpmax`, `sL`, `IL1`, `VL0`, `VL1`, `VA0`, `VA1`, `Vhvmax` | `ir`, `ii`, `p`, `q`, `vt`, `vm`, `ip`, `iq`, `iqextra`, `il`, `lp`, `up` `Tgov1 ` | the TGOV1 governor model | `pmech`, `speed` | `R`, `T1`, `T2`, `T3`, `Pvmax`, `Pvmin`, `Dt` | `none` `Ieeet1` | the IEEET1 exciter model | `bus`, `speed`, `efd`, `vs`\* | `Tr`, `Ka`, `Ta`, `Ke`, `Te`, `Kf`, `Tf`, `Vrmin`, `Vrmax`, `E1`, `E2`, `Se1`, `Se2`, `Ispdlim` | `efd`, `ksat` `SexsPti` | the SEXS-PTI simplified exciter model | `bus`, `efd`, `vs`\* | `Ta`, `Tb`, `Te`, `K`, `Efdmax`, `Efdmin` | `efd` @@ -155,38 +155,6 @@ are specified: Ports marked with \* are optional and, if missing, will be assumed to be connected to a constant value. This list is subject to change. -For `Regca`, `ipcmd` and `iqcmd` are optional in the skeleton implementation -and default to constant zero command inputs when omitted. - -Compact `Regca` device example: - -```json -{ - "class": "Regca", - "ports": { "bus": 1, "ipcmd": 10, "iqcmd": 11 }, - "id": "CV1", - "params": { - "P0": 1.0, - "Q0": 0.0, - "Sconv": 100.0, - "Tg": 0.02, - "TM": 0.02, - "Rqmax": 999.0, - "Rqmin": -999.0, - "Rpmax": 999.0, - "sL": true, - "IL1": 1.1, - "VL0": 0.4, - "VL1": 0.9, - "VA0": 0.4, - "VA1": 0.9, - "Vhvmax": 1.2 - }, - "mon": ["ir", "ii"] -} -``` - - For `Branch`, `tap` and `phase` are optional parameters. If omitted, `tap` defaults to `1.0` and `phase` defaults to `0.0` radians. Bus `bus1` is the tap side for off-nominal transformer branches. diff --git a/tests/UnitTests/Math/SmoothnessIndicatorTests.hpp b/tests/UnitTests/Math/SmoothnessIndicatorTests.hpp index 1ce2d8b1f..9c85c6388 100644 --- a/tests/UnitTests/Math/SmoothnessIndicatorTests.hpp +++ b/tests/UnitTests/Math/SmoothnessIndicatorTests.hpp @@ -55,7 +55,31 @@ namespace GridKit return success.report(__func__); } - TestOutcome deadband() + TestOutcome deadband1() + { + TestStatus success = true; + + const ScalarT lower = scalar(-0.05); + const ScalarT upper = scalar(0.10); + const ScalarT tol = scalar(kSmoothTolerance); + + // Type 1 (no-offset) deadband returns ~0 inside the band ... + const ScalarT inside_input = scalar(0.02); + success *= (std::abs(Math::deadband1(inside_input, lower, upper)) < tol * tol); + + // ... and passes the input through unchanged (no offset) outside the band. + const ScalarT far_above_input = scalar(4.0); + const ScalarT far_below_input = scalar(-4.0); + + success *= std::isfinite(Math::deadband1(far_above_input, lower, upper)); + success *= within(Math::deadband1(far_above_input, lower, upper), far_above_input, tol); + success *= std::isfinite(Math::deadband1(far_below_input, lower, upper)); + success *= within(Math::deadband1(far_below_input, lower, upper), far_below_input, tol); + + return success.report(__func__); + } + + TestOutcome deadband2() { TestStatus success = true; @@ -69,12 +93,12 @@ namespace GridKit const ScalarT expected_below = below_input - lower; const ScalarT expected_above = above_input - upper; - success *= within(Math::deadband(below_input, lower, upper), expected_below, tol); - success *= (std::abs(Math::deadband(inside_input, lower, upper)) < tol * tol); - success *= within(Math::deadband(above_input, lower, upper), expected_above, tol); + success *= within(Math::deadband2(below_input, lower, upper), expected_below, tol); + success *= (std::abs(Math::deadband2(inside_input, lower, upper)) < tol * tol); + success *= within(Math::deadband2(above_input, lower, upper), expected_above, tol); - const ScalarT lower_breakpoint = Math::deadband(lower, lower, upper); - const ScalarT upper_breakpoint = Math::deadband(upper, lower, upper); + const ScalarT lower_breakpoint = Math::deadband2(lower, lower, upper); + const ScalarT upper_breakpoint = Math::deadband2(upper, lower, upper); success *= (lower_breakpoint < scalar(0.0)); success *= (upper_breakpoint > scalar(0.0)); @@ -82,7 +106,7 @@ namespace GridKit success *= (std::abs(upper_breakpoint) < tol); const ScalarT x = scalar(-0.4); - success *= (std::abs(Math::deadband(x, lower, upper) + success *= (std::abs(Math::deadband2(x, lower, upper) - (x - Math::clamp(x, lower, upper))) < scalar(kRoundoffTolerance)); @@ -91,17 +115,17 @@ namespace GridKit const ScalarT expected_far_above = far_above_input - upper; const ScalarT expected_far_below = far_below_input - lower; - success *= std::isfinite(Math::deadband(far_above_input, lower, upper)); - success *= within(Math::deadband(far_above_input, lower, upper), expected_far_above, tol); - success *= std::isfinite(Math::deadband(far_below_input, lower, upper)); - success *= within(Math::deadband(far_below_input, lower, upper), expected_far_below, tol); + success *= std::isfinite(Math::deadband2(far_above_input, lower, upper)); + success *= within(Math::deadband2(far_above_input, lower, upper), expected_far_above, tol); + success *= std::isfinite(Math::deadband2(far_below_input, lower, upper)); + success *= within(Math::deadband2(far_below_input, lower, upper), expected_far_below, tol); const ScalarT point = scalar(0.25); const ScalarT above_point = scalar(0.75); const ScalarT below_point = scalar(-0.25); - success *= (std::abs(Math::deadband(above_point, point, point) - (above_point - point)) + success *= (std::abs(Math::deadband2(above_point, point, point) - (above_point - point)) < scalar(kRoundoffTolerance)); - success *= (std::abs(Math::deadband(below_point, point, point) - (below_point - point)) + success *= (std::abs(Math::deadband2(below_point, point, point) - (below_point - point)) < scalar(kRoundoffTolerance)); return success.report(__func__); diff --git a/tests/UnitTests/Math/runSmoothnessIndicatorTests.cpp b/tests/UnitTests/Math/runSmoothnessIndicatorTests.cpp index 233c36ddd..aecdc7ae8 100644 --- a/tests/UnitTests/Math/runSmoothnessIndicatorTests.cpp +++ b/tests/UnitTests/Math/runSmoothnessIndicatorTests.cpp @@ -7,7 +7,8 @@ int main() GridKit::Testing::SmoothnessIndicatorTests test; result += test.clamp(); - result += test.deadband(); + result += test.deadband1(); + result += test.deadband2(); result += test.limitIndicators(); result += test.slew(); result += test.linseg(); diff --git a/tests/UnitTests/PhasorDynamics/ConverterRegcaTests.hpp b/tests/UnitTests/PhasorDynamics/ConverterRegcaTests.hpp index f74da86e6..94fc10574 100644 --- a/tests/UnitTests/PhasorDynamics/ConverterRegcaTests.hpp +++ b/tests/UnitTests/PhasorDynamics/ConverterRegcaTests.hpp @@ -1,10 +1,12 @@ #pragma once +#include #include -#include #include #include +#include +#include #include #include #include @@ -14,6 +16,7 @@ #include #include #include +#include namespace GridKit { @@ -28,7 +31,7 @@ namespace GridKit ConverterRegcaTests() = default; ~ConverterRegcaTests() = default; - static constexpr ScalarT kTol = static_cast(1.0e-14); + static constexpr ScalarT kTol = static_cast(1.0e-10); TestOutcome constructor() { @@ -39,11 +42,42 @@ namespace GridKit PhasorDynamics::Converter::Regca minimal(&bus); success *= (minimal.size() == static_cast(PhasorDynamics::Converter::RegcaInternalVariables::MAXIMUM)); success *= (minimal.getMonitor() == nullptr); + success *= (minimal.verify() > 0); auto data = makeTestData(); PhasorDynamics::Converter::Regca from_data(&bus, data); success *= (from_data.size() == static_cast(PhasorDynamics::Converter::RegcaInternalVariables::MAXIMUM)); success *= (from_data.getMonitor() != nullptr); + success *= (from_data.verify() == 0); + + return success.report(__func__); + } + + TestOutcome parameterValidation() + { + TestStatus success = true; + + PhasorDynamics::Bus bus(1.0, 0.0); + + auto missing = makeTestData(); + missing.parameters.erase(Params::Tg); + PhasorDynamics::Converter::Regca missing_model(&bus, missing); + success *= (missing_model.verify() > 0); + + auto bad_switch = makeTestData(); + bad_switch.parameters[Params::sL] = static_cast(2); + PhasorDynamics::Converter::Regca bad_switch_model(&bus, bad_switch); + success *= (bad_switch_model.verify() > 0); + + success *= invalidParameterCase(bus, Params::mva, static_cast(0.0)); + success *= invalidParameterCase(bus, Params::Tg, static_cast(0.0)); + success *= invalidParameterCase(bus, Params::TM, static_cast(0.0)); + success *= invalidParameterCase(bus, Params::Rpmax, static_cast(0.0)); + success *= invalidParameterCase(bus, Params::Rqmin, static_cast(0.0)); + success *= invalidParameterCase(bus, Params::IL1, static_cast(-0.1)); + success *= invalidParameterCase(bus, Params::VL1, static_cast(0.3)); + success *= invalidParameterCase(bus, Params::VA1, static_cast(0.3)); + success *= invalidParameterCase(bus, Params::Vhvmax, static_cast(0.0)); return success.report(__func__); } @@ -65,14 +99,10 @@ namespace GridKit success *= (regca.evaluateResidual() == 0); success *= (regca.evaluateJacobian() == 0); - const auto& y = regca.y(); - const auto& yp = regca.yp(); - const auto& f = regca.getResidual(); for (size_t i = 0; i < static_cast(regca.size()); ++i) { - success *= isEqual(y[i], static_cast(0.0), kTol); - success *= isEqual(yp[i], static_cast(0.0), kTol); - success *= isEqual(f[i], static_cast(0.0), kTol); + success *= isEqual(regca.yp()[i], static_cast(0.0), kTol); + success *= isEqual(regca.getResidual()[i], static_cast(0.0), kTol); } using Vars = PhasorDynamics::Converter::RegcaInternalVariables; @@ -90,6 +120,167 @@ namespace GridKit return success.report(__func__); } + TestOutcome steadyStateInitializationGolden() + { + TestStatus success = true; + + PhasorDynamics::Bus bus(0.8, 0.6); + bus.allocate(); + bus.initialize(); + + auto data = makeGoldenTestData(true); + + PhasorDynamics::Converter::Regca regca(&bus, data); + + PhasorDynamics::SignalNode ipcmd_node; + PhasorDynamics::SignalNode iqcmd_node; + ScalarT ipcmd_value{0.92}; + ScalarT iqcmd_value{-0.44}; + IdxT ipcmd_index = 21; + IdxT iqcmd_index = 22; + ipcmd_node.set(&ipcmd_value, &ipcmd_index); + iqcmd_node.set(&iqcmd_value, &iqcmd_index); + regca.getSignals().template attachSignalNode(&ipcmd_node); + regca.getSignals().template attachSignalNode(&iqcmd_node); + + regca.allocate(); + success *= (regca.initialize() == 0); + + std::vector expected_y(static_cast(Vars::MAXIMUM), static_cast(0.0)); + expected_y[index(Vars::VM)] = static_cast(1.0); + expected_y[index(Vars::IQ)] = static_cast(-0.44); + expected_y[index(Vars::IP)] = static_cast(0.92); + expected_y[index(Vars::VT)] = static_cast(1.0); + expected_y[index(Vars::II)] = static_cast(0.904); + expected_y[index(Vars::IQEXTRA)] = static_cast(0.0); + expected_y[index(Vars::IL)] = static_cast(1.1); + expected_y[index(Vars::IR)] = static_cast(0.472); + expected_y[index(Vars::LP)] = static_cast(-70.0); + expected_y[index(Vars::UP)] = static_cast(0.7); + expected_y[index(Vars::PBR)] = static_cast(0.92); + expected_y[index(Vars::QBR)] = static_cast(-0.44); + + success *= vectorMatches(regca.y(), expected_y, "REGCA initialization state"); + for (size_t i = 0; i < static_cast(regca.size()); ++i) + { + success *= isEqual(regca.yp()[i], static_cast(0.0), kTol); + } + + return success.report(__func__); + } + + TestOutcome attachedSignalInitialization() + { + TestStatus success = true; + + PhasorDynamics::Bus bus(1.0, 0.0); + bus.allocate(); + bus.initialize(); + + auto data = makeTestData(); + PhasorDynamics::Converter::Regca regca(&bus, data); + + PhasorDynamics::SignalNode ipcmd_node; + PhasorDynamics::SignalNode iqcmd_node; + ScalarT ipcmd_value{0.75}; + ScalarT iqcmd_value{-0.20}; + IdxT ipcmd_index = 21; + IdxT iqcmd_index = 22; + + ipcmd_node.set(&ipcmd_value, &ipcmd_index); + iqcmd_node.set(&iqcmd_value, &iqcmd_index); + regca.getSignals().template attachSignalNode(&ipcmd_node); + regca.getSignals().template attachSignalNode(&iqcmd_node); + + regca.allocate(); + success *= (regca.initialize() == 0); + + success *= isEqual(ipcmd_value, static_cast(0.75), kTol); + success *= isEqual(iqcmd_value, static_cast(-0.20), kTol); + success *= isEqual(regca.y()[index(Vars::IP)], static_cast(0.75), kTol); + success *= isEqual(regca.y()[index(Vars::IQ)], static_cast(-0.20), kTol); + + return success.report(__func__); + } + + TestOutcome fallbackInit() + { + TestStatus success = true; + + PhasorDynamics::Bus bus(0.8, 0.6); + bus.allocate(); + bus.initialize(); + + auto data = makeTestData(); + data.parameters[Params::mva] = static_cast(50.0); + data.parameters[Params::IL1] = static_cast(2.0); + data.parameters[Params::P0] = static_cast(0.8); + data.parameters[Params::Q0] = static_cast(-0.1); + + PhasorDynamics::Converter::Regca regca(&bus, data); + regca.allocate(); + success *= (regca.initialize() == 0); + + success *= isEqual(regca.y()[index(Vars::IR)], static_cast(1.16), kTol); + success *= isEqual(regca.y()[index(Vars::II)], static_cast(1.12), kTol); + success *= isEqual(regca.y()[index(Vars::PBR)], static_cast(0.8), kTol); + success *= isEqual(regca.y()[index(Vars::QBR)], static_cast(-0.1), kTol); + success *= isEqual(regca.y()[index(Vars::IP)], static_cast(1.6), kTol); + success *= isEqual(regca.y()[index(Vars::IQ)], static_cast(-0.2), kTol); + + bus.evaluateResidual(); + success *= (regca.evaluateResidual() == 0); + + for (size_t i = 0; i < static_cast(regca.size()); ++i) + { + success *= isEqual(regca.getResidual()[i], static_cast(0.0), kTol); + } + success *= isEqual(bus.Ir(), static_cast(0.58), kTol); + success *= isEqual(bus.Ii(), static_cast(0.56), kTol); + + return success.report(__func__); + } + + TestOutcome invalidInitialization() + { + TestStatus success = true; + + auto data = makeTestData(); + + PhasorDynamics::Bus high_voltage_bus(1.3, 0.0); + high_voltage_bus.allocate(); + high_voltage_bus.initialize(); + PhasorDynamics::Converter::Regca high_voltage_regca(&high_voltage_bus, data); + high_voltage_regca.allocate(); + success *= (high_voltage_regca.initialize() > 0); + + auto low_voltage_fallback_data = data; + low_voltage_fallback_data.parameters[Params::P0] = static_cast(0.75); + PhasorDynamics::Bus low_voltage_fallback_bus(0.2, 0.0); + low_voltage_fallback_bus.allocate(); + low_voltage_fallback_bus.initialize(); + PhasorDynamics::Converter::Regca low_voltage_fallback_regca( + &low_voltage_fallback_bus, low_voltage_fallback_data); + low_voltage_fallback_regca.allocate(); + success *= (low_voltage_fallback_regca.initialize() > 0); + + PhasorDynamics::Bus low_voltage_attached_bus(0.2, 0.0); + low_voltage_attached_bus.allocate(); + low_voltage_attached_bus.initialize(); + PhasorDynamics::Converter::Regca low_voltage_attached_regca( + &low_voltage_attached_bus, low_voltage_fallback_data); + PhasorDynamics::SignalNode ipcmd_node; + ScalarT ipcmd_value{0.75}; + IdxT ipcmd_index = 21; + ipcmd_node.set(&ipcmd_value, &ipcmd_index); + low_voltage_attached_regca.getSignals().template attachSignalNode( + &ipcmd_node); + low_voltage_attached_regca.allocate(); + success *= (low_voltage_attached_regca.initialize() == 0); + + return success.report(__func__); + } + TestOutcome signalVerification() { TestStatus success = true; @@ -104,13 +295,13 @@ namespace GridKit IdxT ipcmd_index = 21; IdxT iqcmd_index = 22; - regca.getSignals().template attachSignalNode(&ipcmd_node); + regca.getSignals().template attachSignalNode(&ipcmd_node); success *= (regca.verify() > 0); ipcmd_node.set(&ipcmd_value, &ipcmd_index); success *= (regca.verify() == 0); - regca.getSignals().template attachSignalNode(&iqcmd_node); + regca.getSignals().template attachSignalNode(&iqcmd_node); success *= (regca.verify() > 0); iqcmd_node.set(&iqcmd_value, &iqcmd_index); @@ -123,12 +314,92 @@ namespace GridKit { TestStatus success = true; - PhasorDynamics::Converter::Regca regca(nullptr); + PhasorDynamics::Converter::Regca regca(nullptr, makeTestData()); success *= (regca.verify() > 0); return success.report(__func__); } + TestOutcome residualGoldenVectors() + { + TestStatus success = true; + + std::vector positive_q_lvpl(static_cast(Vars::MAXIMUM), static_cast(0.0)); + positive_q_lvpl[index(Vars::VM)] = static_cast(0.29); + positive_q_lvpl[index(Vars::IQ)] = static_cast(0.52); + positive_q_lvpl[index(Vars::IP)] = static_cast(0.21999997439919405); + positive_q_lvpl[index(Vars::VT)] = static_cast(0.0046000000000000485); + positive_q_lvpl[index(Vars::II)] = static_cast(0.2546); + positive_q_lvpl[index(Vars::IQEXTRA)] = static_cast(-0.03); + positive_q_lvpl[index(Vars::IL)] = static_cast(0.29199937917427254); + positive_q_lvpl[index(Vars::IR)] = static_cast(0.26); + positive_q_lvpl[index(Vars::LP)] = static_cast(-69.6); + positive_q_lvpl[index(Vars::UP)] = static_cast(-0.2999999999999803); + positive_q_lvpl[index(Vars::PBR)] = static_cast(0.0); + positive_q_lvpl[index(Vars::QBR)] = static_cast(0.0); + + std::vector nonpositive_q_no_lvpl(static_cast(Vars::MAXIMUM), static_cast(0.0)); + nonpositive_q_no_lvpl[index(Vars::VM)] = static_cast(0.29); + nonpositive_q_no_lvpl[index(Vars::IQ)] = static_cast(1.5200000000000002); + nonpositive_q_no_lvpl[index(Vars::IP)] = static_cast(0.21999997439919405); + nonpositive_q_no_lvpl[index(Vars::VT)] = static_cast(0.0046000000000000485); + nonpositive_q_no_lvpl[index(Vars::II)] = static_cast(0.2546); + nonpositive_q_no_lvpl[index(Vars::IQEXTRA)] = static_cast(-0.03); + nonpositive_q_no_lvpl[index(Vars::IL)] = static_cast(0.29199937917427254); + nonpositive_q_no_lvpl[index(Vars::IR)] = static_cast(0.26); + nonpositive_q_no_lvpl[index(Vars::LP)] = static_cast(-69.6); + nonpositive_q_no_lvpl[index(Vars::UP)] = static_cast(0.39999999999999997); + nonpositive_q_no_lvpl[index(Vars::PBR)] = static_cast(0.0); + nonpositive_q_no_lvpl[index(Vars::QBR)] = static_cast(0.0); + + success *= residualGoldenVectorCase(static_cast(0.1), true, positive_q_lvpl); + success *= residualGoldenVectorCase(static_cast(-0.1), false, nonpositive_q_no_lvpl); + + return success.report(__func__); + } + + TestOutcome busInjection() + { + TestStatus success = true; + + PhasorDynamics::Bus bus(0.8, 0.6); + bus.allocate(); + bus.initialize(); + + auto data = makeTestData(); + data.parameters[Params::mva] = static_cast(50.0); + + PhasorDynamics::Converter::Regca regca(&bus, data); + + PhasorDynamics::SignalNode ipcmd_node; + PhasorDynamics::SignalNode iqcmd_node; + PhasorDynamics::SignalNode pbranch_node; + PhasorDynamics::SignalNode qbranch_node; + ScalarT ipcmd_value{1.6}; + ScalarT iqcmd_value{-0.2}; + IdxT ipcmd_index = 21; + IdxT iqcmd_index = 22; + ipcmd_node.set(&ipcmd_value, &ipcmd_index); + iqcmd_node.set(&iqcmd_value, &iqcmd_index); + regca.getSignals().template attachSignalNode(&ipcmd_node); + regca.getSignals().template attachSignalNode(&iqcmd_node); + regca.getSignals().template assignSignalNode(&pbranch_node); + regca.getSignals().template assignSignalNode(&qbranch_node); + + regca.allocate(); + success *= (regca.initialize() == 0); + + bus.evaluateResidual(); + success *= (regca.evaluateResidual() == 0); + + success *= isEqual(bus.Ir(), static_cast(0.58), kTol); + success *= isEqual(bus.Ii(), static_cast(0.56), kTol); + success *= isEqual(pbranch_node.read(), static_cast(0.8), kTol); + success *= isEqual(qbranch_node.read(), static_cast(-0.1), kTol); + + return success.report(__func__); + } + TestOutcome jsonParseAndSystemAssembly() { TestStatus success = true; @@ -138,8 +409,8 @@ namespace GridKit "header": { "format_version": 0, "format_revision": 1, - "case_name": "REGCA skeleton", - "case_description": "REGCA parser smoke test", + "case_name": "REGCA full model", + "case_description": "REGCA parser behavior test", "case_comments": "", "freq_base": 60.0, "va_base": 100000000.0 @@ -153,9 +424,9 @@ namespace GridKit "ports": { "bus": 1 }, "id": "CV1", "params": { - "P0": 1.0, + "P0": 0.0, "Q0": 0.0, - "Sconv": 100, + "mva": 100, "Tg": 0.02, "TM": 0.02, "Rqmax": 999.0, @@ -169,18 +440,21 @@ namespace GridKit "VA1": 0.9, "Vhvmax": 1.2 }, - "mon": ["ir", "ii"] + "mon": ["ir", "ii", "p", "q"] } ] } )json"); - auto data = PhasorDynamics::parseSystemModelData(input); - success *= (data.regca.size() == 1); - success *= (data.regca[0].device_class == "Regca"); - success *= (data.regca[0].ports.at(PhasorDynamics::Converter::RegcaPorts::bus) == 1); - success *= (std::get_if(&data.regca[0].parameters.at(PhasorDynamics::Converter::RegcaParameters::Sconv)) != nullptr); - success *= (std::get_if(&data.regca[0].parameters.at(PhasorDynamics::Converter::RegcaParameters::sL)) != nullptr); + auto data = PhasorDynamics::parseSystemModelData(input); + success *= (data.regca.size() == 1); + const auto& regca_data = data.regca[0]; + success *= (regca_data.device_class == "Regca"); + success *= (regca_data.ports.at(PhasorDynamics::Converter::RegcaPorts::bus) == 1); + success *= (std::get_if(®ca_data.parameters.at(Params::P0)) != nullptr); + success *= (std::get_if(®ca_data.parameters.at(Params::Q0)) != nullptr); + success *= (std::get_if(®ca_data.parameters.at(Params::mva)) != nullptr); + success *= (std::get_if(®ca_data.parameters.at(Params::sL)) != nullptr); PhasorDynamics::SystemModel system(data); success *= (system.allocate() == 0); @@ -188,19 +462,50 @@ namespace GridKit success *= (system.tagDifferentiable() == 0); success *= (system.evaluateResidual() == 0); success *= (system.evaluateJacobian() == 0); - success *= (system.size() == 12); + success *= (system.size() == 14); success *= isEqual(system.getResidual()[0], 0.0, static_cast(kTol)); success *= isEqual(system.getResidual()[1], 0.0, static_cast(kTol)); + for (size_t i = 2; i < system.getResidual().size(); ++i) + { + success *= isEqual(system.getResidual()[i], 0.0, static_cast(kTol)); + } return success.report(__func__); } +#ifdef GRIDKIT_ENABLE_ENZYME + TestOutcome jacobian() + { + TestStatus success = true; + + auto dependency_tracking_jacobian = DependencyTrackingJacobian(); + auto enzyme_jacobian = EnzymeJacobian(); + + success *= (dependency_tracking_jacobian.size() == enzyme_jacobian.size()); + const auto nrows = std::min(dependency_tracking_jacobian.size(), enzyme_jacobian.size()); + for (size_t i = 0; i < nrows; ++i) + { + success *= isEqual(dependency_tracking_jacobian[i], enzyme_jacobian[i], static_cast(1.0e-8)); + } + + return success.report(__func__); + } +#endif + private: + using Params = PhasorDynamics::Converter::RegcaParameters; + using Vars = PhasorDynamics::Converter::RegcaInternalVariables; + using Ext = PhasorDynamics::Converter::RegcaExternalVariables; + + static size_t index(Vars variable) + { + return static_cast(variable); + } + auto makeTestData() -> PhasorDynamics::Converter::RegcaData { - using Params = PhasorDynamics::Converter::RegcaParameters; - using Ports = PhasorDynamics::Converter::RegcaPorts; - using Mon = PhasorDynamics::Converter::RegcaMonitorableVariables; + using Ports = PhasorDynamics::Converter::RegcaPorts; + using Mon = PhasorDynamics::Converter::RegcaMonitorableVariables; PhasorDynamics::Converter::RegcaData data; data.device_class = "Regca"; @@ -209,9 +514,7 @@ namespace GridKit data.monitored_variables.insert(Mon::ir); data.monitored_variables.insert(Mon::ii); - data.parameters[Params::P0] = static_cast(1.0); - data.parameters[Params::Q0] = static_cast(0.0); - data.parameters[Params::Sconv] = static_cast(100); + data.parameters[Params::mva] = static_cast(100); data.parameters[Params::Tg] = static_cast(0.02); data.parameters[Params::TM] = static_cast(0.02); data.parameters[Params::Rqmax] = static_cast(999.0); @@ -227,6 +530,234 @@ namespace GridKit return data; } + + auto makeGoldenTestData(bool use_lvpl) -> PhasorDynamics::Converter::RegcaData + { + auto data = makeTestData(); + data.parameters[Params::Tg] = static_cast(0.2); + data.parameters[Params::TM] = static_cast(0.4); + data.parameters[Params::Rqmax] = static_cast(0.5); + data.parameters[Params::Rqmin] = static_cast(-0.6); + data.parameters[Params::Rpmax] = static_cast(0.7); + data.parameters[Params::sL] = use_lvpl; + data.parameters[Params::IL1] = static_cast(1.1); + data.parameters[Params::VL0] = static_cast(0.4); + data.parameters[Params::VL1] = static_cast(0.9); + data.parameters[Params::VA0] = static_cast(0.4); + data.parameters[Params::VA1] = static_cast(0.9); + data.parameters[Params::Vhvmax] = static_cast(1.3); + return data; + } + + bool invalidParameterCase(PhasorDynamics::Bus& bus, Params param, RealT value) + { + auto data = makeTestData(); + data.parameters[param] = value; + PhasorDynamics::Converter::Regca model(&bus, data); + return model.verify() > 0; + } + + bool vectorMatches(const std::vector& actual, + const std::vector& expected, + const char* label) const + { + bool success = (actual.size() == expected.size()); + const auto n = std::min(actual.size(), expected.size()); + for (size_t i = 0; i < n; ++i) + { + if (!isEqual(actual[i], expected[i], kTol)) + { + std::cout << label << " mismatch at row " << i << ": " + << actual[i] << " != " << expected[i] << "\n"; + success = false; + } + } + return success; + } + + void setGoldenResidualState(PhasorDynamics::Converter::Regca& regca) + { + regca.y()[index(Vars::VM)] = static_cast(0.86); + regca.y()[index(Vars::IQ)] = static_cast(-0.2); + regca.y()[index(Vars::IP)] = static_cast(0.85); + regca.y()[index(Vars::VT)] = static_cast(0.98); + regca.y()[index(Vars::II)] = static_cast(0.18); + regca.y()[index(Vars::IQEXTRA)] = static_cast(0.03); + regca.y()[index(Vars::IL)] = static_cast(0.72); + regca.y()[index(Vars::IR)] = static_cast(0.5); + regca.y()[index(Vars::LP)] = static_cast(-0.4); + regca.y()[index(Vars::UP)] = static_cast(0.3); + regca.y()[index(Vars::PBR)] = static_cast(0.52); + regca.y()[index(Vars::QBR)] = static_cast(-0.046); + + regca.yp()[index(Vars::VM)] = static_cast(0.01); + regca.yp()[index(Vars::IQ)] = static_cast(-0.02); + regca.yp()[index(Vars::IP)] = static_cast(0.03); + } + + bool residualGoldenVectorCase(ScalarT initial_iqcmd, bool use_lvpl, const std::vector& expected) + { + bool success = true; + + PhasorDynamics::Bus bus(0.95, 0.25); + bus.allocate(); + bus.initialize(); + + auto data = makeGoldenTestData(use_lvpl); + PhasorDynamics::Converter::Regca regca(&bus, data); + + PhasorDynamics::SignalNode ipcmd_node; + PhasorDynamics::SignalNode iqcmd_node; + ScalarT ipcmd_value{0.9}; + ScalarT iqcmd_value{initial_iqcmd}; + IdxT ipcmd_index = 21; + IdxT iqcmd_index = 22; + ipcmd_node.set(&ipcmd_value, &ipcmd_index); + iqcmd_node.set(&iqcmd_value, &iqcmd_index); + regca.getSignals().template attachSignalNode(&ipcmd_node); + regca.getSignals().template attachSignalNode(&iqcmd_node); + + regca.allocate(); + regca.initialize(); + iqcmd_value = static_cast(0.1); + setGoldenResidualState(regca); + + bus.evaluateResidual(); + regca.evaluateResidual(); + + success *= vectorMatches(regca.getResidual(), expected, "REGCA residual"); + + return success; + } + +#ifdef GRIDKIT_ENABLE_ENZYME + void setJacobianState(PhasorDynamics::Converter::Regca& regca, + PhasorDynamics::Bus& bus) + { + bus.y()[0] = static_cast(0.95); + bus.y()[1] = static_cast(0.25); + + setGoldenResidualState(regca); + } + + void setJacobianStateDep( + PhasorDynamics::Converter::Regca& regca, + PhasorDynamics::Bus& bus) + { + bus.y()[0].setValue(0.95); + bus.y()[1].setValue(0.25); + + regca.y()[index(Vars::VM)].setValue(0.86); + regca.y()[index(Vars::IQ)].setValue(-0.2); + regca.y()[index(Vars::IP)].setValue(0.85); + regca.y()[index(Vars::VT)].setValue(0.98); + regca.y()[index(Vars::II)].setValue(0.18); + regca.y()[index(Vars::IQEXTRA)].setValue(0.03); + regca.y()[index(Vars::IL)].setValue(0.72); + regca.y()[index(Vars::IR)].setValue(0.5); + regca.y()[index(Vars::LP)].setValue(-0.4); + regca.y()[index(Vars::UP)].setValue(0.3); + regca.y()[index(Vars::PBR)].setValue(0.52); + regca.y()[index(Vars::QBR)].setValue(-0.046); + + regca.yp()[index(Vars::VM)].setValue(0.01); + regca.yp()[index(Vars::IQ)].setValue(-0.02); + regca.yp()[index(Vars::IP)].setValue(0.03); + } + + std::vector DependencyTrackingJacobian() + { + using DepVar = DependencyTracking::Variable; + + auto data = makeGoldenTestData(true); + + PhasorDynamics::Bus bus(DepVar{0.95}, DepVar{0.25}); + PhasorDynamics::Converter::Regca regca(&bus, data); + + PhasorDynamics::SignalNode ipcmd_node; + PhasorDynamics::SignalNode iqcmd_node; + DepVar ipcmd_value{0.9}; + DepVar iqcmd_value{0.1}; + IdxT ipcmd_index = static_cast(regca.size() + bus.size()); + IdxT iqcmd_index = static_cast(regca.size() + bus.size() + 1); + + bus.allocate(); + regca.allocate(); + bus.initialize(); + setJacobianStateDep(regca, bus); + + for (IdxT i = 0; i < regca.size(); ++i) + { + regca.y()[static_cast(i)].setVariableNumber(i); + regca.yp()[static_cast(i)].setVariableNumber(i); + } + for (IdxT i = 0; i < bus.size(); ++i) + { + bus.y()[static_cast(i)].setVariableNumber(i + regca.size()); + } + ipcmd_value.setVariableNumber(ipcmd_index); + iqcmd_value.setVariableNumber(iqcmd_index); + + ipcmd_node.set(&ipcmd_value, &ipcmd_index); + iqcmd_node.set(&iqcmd_value, &iqcmd_index); + regca.getSignals().template attachSignalNode(&ipcmd_node); + regca.getSignals().template attachSignalNode(&iqcmd_node); + + bus.evaluateResidual(); + regca.evaluateResidual(); + + std::vector dependencies( + static_cast(regca.size() + bus.size())); + for (IdxT i = 0; i < regca.size(); ++i) + { + dependencies[static_cast(i)] = regca.getResidual()[static_cast(i)].getDependencies(); + } + dependencies[static_cast(regca.size())] = bus.Ir().getDependencies(); + dependencies[static_cast(regca.size() + 1)] = bus.Ii().getDependencies(); + + return dependencies; + } + + std::vector EnzymeJacobian() + { + auto data = makeGoldenTestData(true); + + PhasorDynamics::Bus bus(0.95, 0.25); + PhasorDynamics::Converter::Regca regca(&bus, data); + + PhasorDynamics::SignalNode ipcmd_node; + PhasorDynamics::SignalNode iqcmd_node; + ScalarT ipcmd_value{0.9}; + ScalarT iqcmd_value{0.1}; + IdxT ipcmd_index = static_cast(regca.size() + bus.size()); + IdxT iqcmd_index = static_cast(regca.size() + bus.size() + 1); + + bus.allocate(); + regca.allocate(); + for (IdxT i = 0; i < bus.size(); ++i) + { + bus.setVariableIndex(i, i + regca.size()); + bus.setResidualIndex(i, i + regca.size()); + } + + bus.initialize(); + setJacobianState(regca, bus); + regca.updateTime(0.0, 1.0); + + ipcmd_node.set(&ipcmd_value, &ipcmd_index); + iqcmd_node.set(&iqcmd_value, &iqcmd_index); + regca.getSignals().template attachSignalNode(&ipcmd_node); + regca.getSignals().template attachSignalNode(&iqcmd_node); + + bus.evaluateResidual(); + regca.evaluateResidual(); + regca.evaluateJacobian(); + + auto model_jacobian = regca.getJacobian(); + model_jacobian.deduplicate(); + return MapFromCOO(model_jacobian); + } +#endif }; } // namespace Testing } // namespace GridKit diff --git a/tests/UnitTests/PhasorDynamics/SystemSingleComponentTests.hpp b/tests/UnitTests/PhasorDynamics/SystemSingleComponentTests.hpp index 790d7048a..4d67d39f5 100644 --- a/tests/UnitTests/PhasorDynamics/SystemSingleComponentTests.hpp +++ b/tests/UnitTests/PhasorDynamics/SystemSingleComponentTests.hpp @@ -172,10 +172,10 @@ namespace GridKit PhasorDynamics::SystemModel* system = new PhasorDynamics::SystemModel(); - PhasorDynamics::BusInfinite bus; + PhasorDynamics::BusInfinite bus(1.0, 0.0); system->addBus(&bus); - PhasorDynamics::Converter::Regca regca(&bus); + PhasorDynamics::Converter::Regca regca(&bus, makeRegcaData()); system->addComponent(®ca); success *= system->allocate() == 0; @@ -190,6 +190,33 @@ namespace GridKit return success.report(__func__); } + private: + auto makeRegcaData() -> PhasorDynamics::Converter::RegcaData + { + using Params = PhasorDynamics::Converter::RegcaParameters; + + PhasorDynamics::Converter::RegcaData data; + data.device_class = "Regca"; + data.disambiguation_string = "regca_test"; + data.parameters[Params::P0] = static_cast(1.0); + data.parameters[Params::Q0] = static_cast(0.0); + data.parameters[Params::mva] = static_cast(100.0); + data.parameters[Params::Tg] = static_cast(0.02); + data.parameters[Params::TM] = static_cast(0.02); + data.parameters[Params::Rqmax] = static_cast(999.0); + data.parameters[Params::Rqmin] = static_cast(-999.0); + data.parameters[Params::Rpmax] = static_cast(999.0); + data.parameters[Params::sL] = true; + data.parameters[Params::IL1] = static_cast(1.1); + data.parameters[Params::VL0] = static_cast(0.4); + data.parameters[Params::VL1] = static_cast(0.9); + data.parameters[Params::VA0] = static_cast(0.4); + data.parameters[Params::VA1] = static_cast(0.9); + data.parameters[Params::Vhvmax] = static_cast(1.2); + return data; + } + + public: TestOutcome genrou() { TestStatus success = true; diff --git a/tests/UnitTests/PhasorDynamics/runConverterRegcaTests.cpp b/tests/UnitTests/PhasorDynamics/runConverterRegcaTests.cpp index ae8462396..0d3daf3b4 100644 --- a/tests/UnitTests/PhasorDynamics/runConverterRegcaTests.cpp +++ b/tests/UnitTests/PhasorDynamics/runConverterRegcaTests.cpp @@ -7,10 +7,20 @@ int main() GridKit::Testing::ConverterRegcaTests test; result += test.constructor(); + result += test.parameterValidation(); result += test.lifecycle(); + result += test.steadyStateInitializationGolden(); + result += test.attachedSignalInitialization(); + result += test.fallbackInit(); + result += test.invalidInitialization(); result += test.signalVerification(); result += test.nullBusVerification(); + result += test.residualGoldenVectors(); + result += test.busInjection(); result += test.jsonParseAndSystemAssembly(); +#ifdef GRIDKIT_ENABLE_ENZYME + result += test.jacobian(); +#endif return result.summary(); } From 89aac1edecd5f7ea54db23c015c5c8e4a17bd81a Mon Sep 17 00:00:00 2001 From: lukelowry Date: Mon, 1 Jun 2026 23:57:57 -0500 Subject: [PATCH 3/3] Doc polish, Corrections and Update CHANGELOG --- CHANGELOG.md | 1 + .../Converter/REGCA/CMakeLists.txt | 69 +- .../PhasorDynamics/Converter/REGCA/README.md | 287 +++---- .../PhasorDynamics/Converter/REGCA/Regca.cpp | 4 +- .../PhasorDynamics/Converter/REGCA/Regca.hpp | 66 +- .../Converter/REGCA/RegcaData.hpp | 6 +- .../REGCA/RegcaDependencyTracking.cpp | 4 +- .../Converter/REGCA/RegcaEnzyme.cpp | 4 +- .../Converter/REGCA/RegcaImpl.hpp | 227 +++--- tests/UnitTests/PhasorDynamics/CMakeLists.txt | 4 +- .../PhasorDynamics/ConverterRegcaTests.hpp | 726 +++++++++++------- .../PhasorDynamics/runConverterRegcaTests.cpp | 15 +- 12 files changed, 757 insertions(+), 656 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 066786e09..6911fb4fd 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -62,6 +62,7 @@ - Added cmake-format hooks, including in pre-commit. - Added off-nominal tap ratio and phase shift support to the PhasorDynamics `Branch` model. - Added portable Vector class to GridKit +- Added `REGCA` converter model implementation for PhasorDynamics. ## v0.1 diff --git a/GridKit/Model/PhasorDynamics/Converter/REGCA/CMakeLists.txt b/GridKit/Model/PhasorDynamics/Converter/REGCA/CMakeLists.txt index 9e9b582d7..3ffe1fef7 100644 --- a/GridKit/Model/PhasorDynamics/Converter/REGCA/CMakeLists.txt +++ b/GridKit/Model/PhasorDynamics/Converter/REGCA/CMakeLists.txt @@ -3,47 +3,52 @@ # - Luke Lowery # ]] -set(_install_headers - Regca.hpp - RegcaData.hpp) +set(_install_headers Regca.hpp RegcaData.hpp) if(GRIDKIT_ENABLE_ENZYME) - gridkit_add_library(phasor_dynamics_converter_regca - SOURCES - RegcaEnzyme.cpp - HEADERS - ${_install_headers} - INCLUDE_DIRECTORIES - PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include + gridkit_add_library( + phasor_dynamics_converter_regca + SOURCES RegcaEnzyme.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 + PUBLIC + GridKit::phasor_dynamics_core + PUBLIC + GridKit::phasor_dynamics_signal + PRIVATE + ClangEnzymeFlags COMPILE_OPTIONS - PRIVATE -mllvm -enzyme-auto-sparsity=1 -fno-math-errno) + PRIVATE + -mllvm + -enzyme-auto-sparsity=1 + -fno-math-errno) else() - gridkit_add_library(phasor_dynamics_converter_regca - SOURCES - Regca.cpp - HEADERS - ${_install_headers} - INCLUDE_DIRECTORIES - PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include + gridkit_add_library( + phasor_dynamics_converter_regca + SOURCES Regca.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) + PUBLIC + GridKit::phasor_dynamics_core + PUBLIC + GridKit::phasor_dynamics_signal) endif() -gridkit_add_library(phasor_dynamics_converter_regca_dependency_tracking - SOURCES - RegcaDependencyTracking.cpp - INCLUDE_DIRECTORIES - PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include +gridkit_add_library( + phasor_dynamics_converter_regca_dependency_tracking + SOURCES RegcaDependencyTracking.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) + PUBLIC + GridKit::phasor_dynamics_core + PUBLIC + GridKit::phasor_dynamics_signal_dependency_tracking) -target_link_libraries(phasor_dynamics_components +target_link_libraries( + phasor_dynamics_components INTERFACE GridKit::phasor_dynamics_converter_regca) -target_link_libraries(phasor_dynamics_components_dependency_tracking +target_link_libraries( + phasor_dynamics_components_dependency_tracking INTERFACE GridKit::phasor_dynamics_converter_regca_dependency_tracking) diff --git a/GridKit/Model/PhasorDynamics/Converter/REGCA/README.md b/GridKit/Model/PhasorDynamics/Converter/REGCA/README.md index c2961c13f..66500962f 100644 --- a/GridKit/Model/PhasorDynamics/Converter/REGCA/README.md +++ b/GridKit/Model/PhasorDynamics/Converter/REGCA/README.md @@ -5,10 +5,10 @@ inverter-coupled resources. In GridKit it is represented as a controlled current source at the network interface. Notes: -- LVACM uses the unfiltered terminal voltage $V_T$; LVPL uses the filtered voltage $V_M$. -- Internal currents are on converter base; bus injections are converted to system base in the network interface. -- Exported branch power measurements are on system base. -- HVRCM is represented by internal algebraic current $I_{\mathrm{q}}^{\mathrm{extra}}$. +- Internal currents are on converter base. +- Bus injections and exported branch powers are on system base. +- LVACM uses $V_T$; LVPL uses $V_M$. +- HVRCM is represented by algebraic current $I_{\mathrm{q}}^{\mathrm{extra}}$. ## Block Diagram @@ -22,23 +22,26 @@ Standard REGCA converter-interface model. ## Model Parameters -Symbol | Units | Description | Typical Value | Note ----------------------------------|----------|-------------------------------------------------------|---------------|------ -$P_0$ | [p.u.] | Initial active power injection on system base | TBD | JSON key: `P0`; fallback command source -$Q_0$ | [p.u.] | Initial reactive power injection on system base | TBD | JSON key: `Q0`; fallback command source -$S^{\text{base}}$ | [MVA] | REGCA model power base | TBD | JSON key: `mva` -$T_{\mathrm{g}}$ | [sec] | Converter current-control lag time constant | TBD | -$T_M$ | [sec] | Terminal voltage sensor time constant | TBD | Block name: `Tfltr` -$R_{\mathrm{q}}^{\max}$ | [p.u./s] | Reactive-current recovery positive rate limit | TBD | Block name: `Iqrmax` -$R_{\mathrm{q}}^{\min}$ | [p.u./s] | Reactive-current recovery negative rate limit | TBD | Block name: `Iqrmin` -$R_{\mathrm{p}}^{\max}$ | [p.u./s] | Active-current magnitude recovery rate limit | TBD | Block name: `rrpwr` -$s_L$ | [binary] | LVPL switch | TBD | Block name: `LPVLSW` -$I_{L1}$ | [p.u.] | LVPL upper-current ceiling | TBD | Block name: `Lvpl1` -$V_{L0}$ | [p.u.] | LVPL zero-crossing voltage | TBD | Block name: `xerox` -$V_{L1}$ | [p.u.] | LVPL upper breakpoint voltage | TBD | Block name: `brkpt` -$V_{A0}$ | [p.u.] | LVACM lower breakpoint voltage | TBD | Block name: `Lvpnt0` -$V_{A1}$ | [p.u.] | LVACM upper breakpoint voltage | TBD | Block name: `Lvpnt1` -$V_{\mathrm{hv}}^{\max}$ | [p.u.] | Terminal-voltage ceiling for HV reactive management | TBD | Block name: `Vlim` +Symbol | Units | JSON | Description | Typical Value | Note +---------------------------------|----------|----------|-------------------------------------------------------|---------------|------ +$P_0$ | [p.u.] | `P0` | Initial active power injection on system base | 1.0 | Required initialization source +$Q_0$ | [p.u.] | `Q0` | Initial reactive power injection on system base | 0.0 | Required initialization source +$S^{\mathrm{base}}$ | [MVA] | `mva` | REGCA model power base | 100.0 | +$T_{\mathrm{g}}$ | [sec] | `Tg` | Converter current-control lag time constant | 0.02 | Block name: `Tg` +$T_M$ | [sec] | `TM` | Terminal voltage sensor time constant | 0.02 | Block name: `Tfltr` +$R_{\mathrm{q}}^{\max}$ | [p.u./s] | `Rqmax` | Reactive-current recovery positive rate limit | 999.0 | Block name: `Iqrmax`; GridKit requires $R_{\mathrm{q}}^{\max} > 0$ +$R_{\mathrm{q}}^{\min}$ | [p.u./s] | `Rqmin` | Reactive-current recovery negative rate limit | -999.0 | Block name: `Iqrmin`; GridKit requires $R_{\mathrm{q}}^{\min} < 0$ +$R_{\mathrm{p}}^{\max}$ | [p.u./s] | `Rpmax` | Active-current magnitude recovery rate limit | 999.0 | Block name: `rrpwr` +$s_L$ | [binary] | `sL` | LVPL switch | 1 | Block name: `LPVLSW` +$I_{L1}$ | [p.u.] | `IL1` | LVPL upper-current ceiling | 1.1 | Block name: `LVPL1` +$V_{L0}$ | [p.u.] | `VL0` | LVPL zero-crossing voltage | 0.4 | Block name: `zerox` +$V_{L1}$ | [p.u.] | `VL1` | LVPL upper breakpoint voltage | 0.9 | Block name: `brkpt` +$V_{A0}$ | [p.u.] | `VA0` | LVACM lower breakpoint voltage | 0.4 | Block name: `LVPnt0` +$V_{A1}$ | [p.u.] | `VA1` | LVACM upper breakpoint voltage | 0.9 | Block name: `LVPnt1` +$V_{\mathrm{hv}}^{\max}$ | [p.u.] | `Vhvmax` | Terminal-voltage ceiling for HV reactive management | 1.2 | Block name: `VLim` + +PowerWorld fields `Qmin`, `Khv`, and `Xe` are not parameters of this GridKit +REGCA implementation. ### Parameter Validation @@ -46,7 +49,7 @@ Implementations should reject or report invalid parameter sets: ```math \begin{aligned} - S^{\text{base}} &> 0 & + S^{\mathrm{base}} &> 0 & T_{\mathrm{g}} &> 0 & T_M &> 0 \\ R_{\mathrm{p}}^{\max} &> 0 & @@ -61,16 +64,14 @@ Implementations should reject or report invalid parameter sets: ### Model Derived Parameters -The smooth active-current bound equations use $M_{\mathrm{p}}$, a numerical -relaxation for inactive $\pm\infty$ rate bounds: +The active-current bounds use finite surrogate $M_{\mathrm{p}}$ for inactive +$\pm\infty$ limits: ```math M_{\mathrm{p}} = 100 R_{\mathrm{p}}^{\max} ``` -$M_{\mathrm{p}}$ is not a physical REGCA parameter; it should be large enough -that inactive bounds do not bind expected $f_{\mathrm{p}}$ values while staying -moderate enough to keep the smooth clamp well conditioned. +$M_{\mathrm{p}}$ is not a physical REGCA parameter. ## Model Variables @@ -112,134 +113,67 @@ $V_{\mathrm{i}}$ | [p.u.] | Terminal voltage, imaginary component $I_{\mathrm{p}}^{\mathrm{cmd}}$ | [p.u.] | Active-current command in the terminal-voltage reference frame | Converter base; owned by REEC, constant if no REEC is connected $I_{\mathrm{q}}^{\mathrm{cmd}}$ | [p.u.] | Reactive-current command in the terminal-voltage reference frame | Converter base; owned by REEC, constant if no REEC is connected -When either `ipcmd` or `iqcmd` port is omitted, REGCA derives the missing -constant command from the optional `P0`/`Q0` pair. If `P0`/`Q0` are omitted, -the fallback commands are zero. +REGCA initializes its current commands from the required `P0`/`Q0` power-flow +injection and writes those resolved commands to any attached `ipcmd`/`iqcmd` +ports. If no controller is connected, the resolved initialization commands are +held constant during residual evaluation. ## Model Equations -Define the pre-limit current derivatives: - -```math -\begin{aligned} - f_{\mathrm{q}} &= \dfrac{1}{T_{\mathrm{g}}}(I_{\mathrm{q}}^{\mathrm{cmd}} - I_{\mathrm{q}}) \\ - f_{\mathrm{p}} &= \dfrac{1}{T_{\mathrm{g}}}(I_{\mathrm{p}}^{\mathrm{cmd}} - I_{\mathrm{p}}) -\end{aligned} -``` - ### Differential Equations -The state equations use CommonMath helper notation: +The state equations use CommonMath helper notation. The $I_{\mathrm{q}}$ +limiter branch is selected by the initialized reactive-current command. ```math \begin{aligned} 0 &= -T_M \dot V_M - V_M + V_T \\ - 0 &= -\dot I_{\mathrm{q}} + + 0 &= -T_{\mathrm{g}}\dot I_{\mathrm{q}} + + (I_{\mathrm{q}}^{\mathrm{cmd}} - I_{\mathrm{q}}) + \begin{cases} - \min(f_{\mathrm{q}}, R_{\mathrm{q}}^{\max}) & I_{\mathrm{q},0}^{\mathrm{cmd}} > 0 \\ - \max(f_{\mathrm{q}}, R_{\mathrm{q}}^{\min}) & I_{\mathrm{q},0}^{\mathrm{cmd}} \le 0 + -\rho(I_{\mathrm{q}}^{\mathrm{cmd}} - I_{\mathrm{q}} - T_{\mathrm{g}}R_{\mathrm{q}}^{\max}) + & I_{\mathrm{q},0}^{\mathrm{cmd}} > 0 \\ + \rho(T_{\mathrm{g}}R_{\mathrm{q}}^{\min} - (I_{\mathrm{q}}^{\mathrm{cmd}} - I_{\mathrm{q}})) + & I_{\mathrm{q},0}^{\mathrm{cmd}} \le 0 \end{cases} \\ - 0 &= -\dot I_{\mathrm{p}} + \text{clamp}(f_{\mathrm{p}}, \ell_{\mathrm{p}}, u_{\mathrm{p}}) + 0 &= -T_{\mathrm{g}}\dot I_{\mathrm{p}} + + T_{\mathrm{g}}\ell_{\mathrm{p}} + + \rho(I_{\mathrm{p}}^{\mathrm{cmd}} - I_{\mathrm{p}} - T_{\mathrm{g}}\ell_{\mathrm{p}}) + - \rho(I_{\mathrm{p}}^{\mathrm{cmd}} - I_{\mathrm{p}} - T_{\mathrm{g}}u_{\mathrm{p}}) \end{aligned} ``` -The $I_{\mathrm{q}}$ branch is selected by the initial reactive-current -command. CommonMath defines the helper targets and smooth approximations for -[min, max, and clamp](../../../../CommonMath.md#derived-functions). +CommonMath defines the [linear segment, ramp, and sigmoid helpers](../../../../CommonMath.md#derived-functions). ### Algebraic Equations -The exact algebraic targets are: - -```math -\begin{aligned} - 0 &= -V_T^2 + V_{\mathrm{r}}^2 + V_{\mathrm{i}}^2 \\ - 0 &= -V_T I_{\mathrm{i}} - + I_{\mathrm{p}}V_{\mathrm{i}} - \begin{cases} - 0 & V_T \le V_{A0} \\ - \dfrac{V_T - V_{A0}}{V_{A1} - V_{A0}} & V_{A0} < V_T < V_{A1} \\ - 1 & V_T \ge V_{A1} - \end{cases} - - (I_{\mathrm{q}} - I_{\mathrm{q}}^{\mathrm{extra}})V_{\mathrm{r}} \\ - 0 &= - \begin{cases} - I_{\mathrm{q}}^{\mathrm{extra}} - & V_T < V_{\mathrm{hv}}^{\max} \\ - V_T - V_{\mathrm{hv}}^{\max} - & I_{\mathrm{q}}^{\mathrm{extra}} > 0 - \end{cases} \\ - 0 &= -I_L - + I_{L1} - \begin{cases} - 0 & V_M \le V_{L0} \\ - \dfrac{V_M - V_{L0}}{V_{L1} - V_{L0}} & V_{L0} < V_M < V_{L1} \\ - 1 & V_M \ge V_{L1} - \end{cases} \\ - 0 &= -V_T I_{\mathrm{r}} - + I_{\mathrm{p}}V_{\mathrm{r}} - \begin{cases} - 0 & V_T \le V_{A0} \\ - \dfrac{V_T - V_{A0}}{V_{A1} - V_{A0}} & V_{A0} < V_T < V_{A1} \\ - 1 & V_T \ge V_{A1} - \end{cases} - + (I_{\mathrm{q}} - I_{\mathrm{q}}^{\mathrm{extra}})V_{\mathrm{i}} \\ - 0 &= -\ell_{\mathrm{p}} + - \begin{cases} - -R_{\mathrm{p}}^{\max} & I_{\mathrm{p}} \le 0 \\ - -\infty & I_{\mathrm{p}} > 0 - \end{cases} \\ - 0 &= -u_{\mathrm{p}} + - \begin{cases} - R_{\mathrm{p}}^{\max} & I_{\mathrm{p}} \ge 0 \ \land\ (s_L = 0 \lor I_{\mathrm{p}} < I_L) \\ - 0 & s_L = 1 \ \land\ I_{\mathrm{p}} \ge I_L \\ - \infty & I_{\mathrm{p}} < 0 - \end{cases} \\ - 0 &= -P_{\mathrm{br}} - + \dfrac{S^{\text{base}}}{S^{\text{sys}}} - \left(V_{\mathrm{r}} I_{\mathrm{r}} + V_{\mathrm{i}} I_{\mathrm{i}}\right) \\ - 0 &= -Q_{\mathrm{br}} - + \dfrac{S^{\text{base}}}{S^{\text{sys}}} - \left(V_{\mathrm{i}} I_{\mathrm{r}} - V_{\mathrm{r}} I_{\mathrm{i}}\right) -\end{aligned} -``` - -GridKit implements the algebraic residuals with smooth $\text{linseg}$, -$\rho$, and $\sigma$ operators; CommonMath defines the -[linear segment, ramp, and sigmoid helpers](../../../../CommonMath.md#derived-functions). -An equivalent residual form is: +The algebraic equations are: ```math \begin{aligned} 0 &= -V_T^2 + V_{\mathrm{r}}^2 + V_{\mathrm{i}}^2 \\ 0 &= -V_T I_{\mathrm{i}} - + I_{\mathrm{p}}\text{linseg}(V_T;\ V_{A0},\ V_{A1},\ 1)V_{\mathrm{i}} - - (I_{\mathrm{q}} - I_{\mathrm{q}}^{\mathrm{extra}})V_{\mathrm{r}} \\ + - (I_{\mathrm{q}} - I_{\mathrm{q}}^{\mathrm{extra}})V_{\mathrm{r}} + + V_{\mathrm{i}}I_{\mathrm{p}}\text{linseg}(V_T;\ V_{A0},\ V_{A1},\ 1) \\ 0 &= -I_{\mathrm{q}}^{\mathrm{extra}} - + \rho\!\left(I_{\mathrm{q}}^{\mathrm{extra}} - (V_{\mathrm{hv}}^{\max} - V_T)\right) \\ + + \rho(V_T - V_{\mathrm{hv}}^{\max}) \\ 0 &= -I_L + \text{linseg}(V_M;\ V_{L0},\ V_{L1},\ I_{L1}) \\ 0 &= -V_T I_{\mathrm{r}} - + I_{\mathrm{p}}\text{linseg}(V_T;\ V_{A0},\ V_{A1},\ 1)V_{\mathrm{r}} - + (I_{\mathrm{q}} - I_{\mathrm{q}}^{\mathrm{extra}})V_{\mathrm{i}} \\ + + (I_{\mathrm{q}} - I_{\mathrm{q}}^{\mathrm{extra}})V_{\mathrm{i}} + + V_{\mathrm{r}}I_{\mathrm{p}}\text{linseg}(V_T;\ V_{A0},\ V_{A1},\ 1) \\ 0 &= -\ell_{\mathrm{p}} - R_{\mathrm{p}}^{\max} - (M_{\mathrm{p}} - R_{\mathrm{p}}^{\max})\sigma(I_{\mathrm{p}}) \\ 0 &= -u_{\mathrm{p}} - + \begin{cases} - M_{\mathrm{p}}(1-\sigma(I_{\mathrm{p}})) - + R_{\mathrm{p}}^{\max}\sigma(I_{\mathrm{p}}) - \sigma(I_L - I_{\mathrm{p}}) - & s_L = 1 \\ - M_{\mathrm{p}}(1-\sigma(I_{\mathrm{p}})) - + R_{\mathrm{p}}^{\max}\sigma(I_{\mathrm{p}}) - & s_L = 0 - \end{cases} \\ + + M_{\mathrm{p}}(1-\sigma(I_{\mathrm{p}})) + + R_{\mathrm{p}}^{\max}\sigma(I_{\mathrm{p}}) + \left(1 - s_L + s_L\sigma(I_L - I_{\mathrm{p}})\right) \\ 0 &= -P_{\mathrm{br}} - + \dfrac{S^{\text{base}}}{S^{\text{sys}}} + + \dfrac{S^{\mathrm{base}}}{S^{\mathrm{sys}}} \left(V_{\mathrm{r}} I_{\mathrm{r}} + V_{\mathrm{i}} I_{\mathrm{i}}\right) \\ 0 &= -Q_{\mathrm{br}} - + \dfrac{S^{\text{base}}}{S^{\text{sys}}} + + \dfrac{S^{\mathrm{base}}}{S^{\mathrm{sys}}} \left(V_{\mathrm{i}} I_{\mathrm{r}} - V_{\mathrm{r}} I_{\mathrm{i}}\right) \end{aligned} ``` @@ -253,8 +187,8 @@ REGCA currents: ```math \begin{aligned} - I_{\mathrm{r}}^{\mathrm{inj}} &:= I_{\mathrm{r}}\dfrac{S^{\text{base}}}{S^{\text{sys}}} \\ - I_{\mathrm{i}}^{\mathrm{inj}} &:= I_{\mathrm{i}}\dfrac{S^{\text{base}}}{S^{\text{sys}}} + I_{\mathrm{r}}^{\mathrm{inj}} &:= I_{\mathrm{r}}\dfrac{S^{\mathrm{base}}}{S^{\mathrm{sys}}} \\ + I_{\mathrm{i}}^{\mathrm{inj}} &:= I_{\mathrm{i}}\dfrac{S^{\mathrm{base}}}{S^{\mathrm{sys}}} \end{aligned} ``` @@ -267,71 +201,52 @@ steady-state initial values: ```math \begin{aligned} - V_T &= \sqrt{V_\mathrm{r}^2 + V_\mathrm{i}^2} \\ - V_{M0} &= V_T \\ - I_{L0} &= \text{linseg}(V_T;\ V_{L0},\ V_{L1},\ I_{L1}) \\ - I_\mathrm{p0} &= - \begin{cases} - I_\mathrm{p}^\mathrm{cmd}(0) & \text{if the } I_\mathrm{p}^\mathrm{cmd} \text{ port is attached} \\ - \dfrac{P_0} - {V_T\text{linseg}(V_T;\ V_{A0},\ V_{A1},\ 1)} - \dfrac{S^{\text{sys}}}{S^{\text{base}}} - & \text{otherwise} - \end{cases} \\ - I_\mathrm{q0} &= - \begin{cases} - I_\mathrm{q}^\mathrm{cmd}(0) & \text{if the } I_\mathrm{q}^\mathrm{cmd} \text{ port is attached} \\ - \dfrac{Q_0}{V_T}\dfrac{S^{\text{sys}}}{S^{\text{base}}} - & \text{otherwise} - \end{cases} \\ - \ell_\mathrm{p0} &= -R_\mathrm{p}^{\max} - - (M_\mathrm{p} - R_\mathrm{p}^{\max})\sigma(I_\mathrm{p0}) \\ - u_\mathrm{p0} &= M_\mathrm{p}(1-\sigma(I_\mathrm{p0})) - + R_\mathrm{p}^{\max}\sigma(I_\mathrm{p0}) - \left(1 - s_L + s_L\sigma(I_{L0} - I_\mathrm{p0})\right) \\ - I_\mathrm{q0}^\mathrm{extra} &= 0 \\ - I_\mathrm{i0} &= \dfrac{ - \text{linseg}(V_T;\ V_{A0},\ V_{A1},\ 1)I_{\mathrm{p0}}V_{\mathrm{i}} - - I_{\mathrm{q0}}V_{\mathrm{r}} - }{V_T} \\ - I_\mathrm{r0} &= \dfrac{ - \text{linseg}(V_T;\ V_{A0},\ V_{A1},\ 1)I_{\mathrm{p0}}V_{\mathrm{r}} - + I_{\mathrm{q0}}V_{\mathrm{i}} - }{V_T} \\ - P_{\mathrm{br},0} &= \dfrac{S^{\text{base}}}{S^{\text{sys}}} - \left(V_{\mathrm{r}} I_{\mathrm{r0}} + V_{\mathrm{i}} I_{\mathrm{i0}}\right) \\ - Q_{\mathrm{br},0} &= \dfrac{S^{\text{base}}}{S^{\text{sys}}} - \left(V_{\mathrm{i}} I_{\mathrm{r0}} - V_{\mathrm{r}} I_{\mathrm{i0}}\right) + V_T &= \sqrt{V_{\mathrm{r}}^2 + V_{\mathrm{i}}^2} \\ + V_{M,0} &= V_T \\ + I_{L,0} &= \text{linseg}(V_T;\ V_{L0},\ V_{L1},\ I_{L1}) \\ + I_{\mathrm{p},0} &= V_T^{-1}P_0 + \dfrac{S^{\mathrm{sys}}}{S^{\mathrm{base}}} + \left[\text{linseg}(V_T;\ V_{A0},\ V_{A1},\ 1)\right]^{-1} \\ + I_{\mathrm{q},0} &= V_T^{-1}Q_0\dfrac{S^{\mathrm{sys}}}{S^{\mathrm{base}}} \\ + \ell_{\mathrm{p},0} &= -R_{\mathrm{p}}^{\max} + - (M_{\mathrm{p}} - R_{\mathrm{p}}^{\max})\sigma(I_{\mathrm{p},0}) \\ + u_{\mathrm{p},0} &= M_{\mathrm{p}}(1-\sigma(I_{\mathrm{p},0})) + + R_{\mathrm{p}}^{\max}\sigma(I_{\mathrm{p},0}) + \left(1 - s_L + s_L\sigma(I_{L,0} - I_{\mathrm{p},0})\right) \\ + I_{\mathrm{q},0}^{\mathrm{extra}} &= 0 \\ + I_{\mathrm{i},0} &= V_T^{-1}\left[ + -I_{\mathrm{q},0}V_{\mathrm{r}} + + V_{\mathrm{i}}I_{\mathrm{p},0}\text{linseg}(V_T;\ V_{A0},\ V_{A1},\ 1) + \right] \\ + I_{\mathrm{r},0} &= V_T^{-1}\left[ + I_{\mathrm{q},0}V_{\mathrm{i}} + + V_{\mathrm{r}}I_{\mathrm{p},0}\text{linseg}(V_T;\ V_{A0},\ V_{A1},\ 1) + \right] \\ + P_{\mathrm{br},0} &= \dfrac{S^{\mathrm{base}}}{S^{\mathrm{sys}}} + \left(V_{\mathrm{r}} I_{\mathrm{r},0} + V_{\mathrm{i}} I_{\mathrm{i},0}\right) \\ + Q_{\mathrm{br},0} &= \dfrac{S^{\mathrm{base}}}{S^{\mathrm{sys}}} + \left(V_{\mathrm{i}} I_{\mathrm{r},0} - V_{\mathrm{r}} I_{\mathrm{i},0}\right) \end{aligned} ``` -For normal power-flow starts, $V_T > V_{A1}$, so -$\text{linseg}(V_T;\ V_{A0},\ V_{A1},\ 1) = 1$ and the -$I_{\mathrm{p0}}$ formula is well defined. -When a command port is omitted, GridKit stores the resolved initial value as -the constant command used during residual evaluation. - -Initialization should verify: -- $V_T < V_{\mathrm{hv}}^{\max}$. If $V_T \ge V_{\mathrm{hv}}^{\max}$, - $I_{\mathrm{q0}}^{\mathrm{extra}} = 0$ may not satisfy the HVRCM algebraic - condition, and a nonzero value should be solved or the initialization rejected. -- $\text{linseg}(V_T;\ V_{A0},\ V_{A1},\ 1) > 0$ when - the `ipcmd` port is omitted and $P_0 \ne 0$. If the LVACM gain is zero, - no finite $I_{\mathrm{p0}}$ can reproduce nonzero initial active power. - -All internal derivatives initialize to zero. +REGCA writes the resolved initial commands to attached `ipcmd` and `iqcmd` +ports. Initialization rejects nonpositive $V_T$, terminal voltage at or above +$V_{\mathrm{hv}}^{\max}$, and nonzero $P_0$ with zero LVACM gain. All internal +derivatives initialize to zero. ## Model Outputs -Real and imaginary injected currents, $I_{\mathrm{r}}$ and -$I_{\mathrm{i}}$, are converter-base algebraic variables. Optional -`ibranchr`, `ibranchi`, `pbranch`, and `qbranch` ports export REGCA-owned -algebraic measurement states for downstream controllers such as REPCA. The -`pbranch` and `qbranch` signal ports export system-base power outputs: -```math -\begin{aligned} - P &= V_{\mathrm{r}} I_{\mathrm{r}}^{\mathrm{inj}} + V_{\mathrm{i}} I_{\mathrm{i}}^{\mathrm{inj}} \\ - Q &= V_{\mathrm{i}} I_{\mathrm{r}}^{\mathrm{inj}} - V_{\mathrm{r}} I_{\mathrm{i}}^{\mathrm{inj}} -\end{aligned} -``` -Power outputs are positive leaving the converter and entering the bus. +Output | Units | Description | Note +----------------|----------|--------------------------------------|------ +`ir` | [p.u.] | Real current injection | Converter base; exported through `ibranchr` when assigned +`ii` | [p.u.] | Imaginary current injection | Converter base; exported through `ibranchi` when assigned +`p` | [p.u.] | Active-power output | System base; exported through `pbranch` when assigned +`q` | [p.u.] | Reactive-power output | System base; exported through `qbranch` when assigned +`vt` | [p.u.] | Terminal voltage magnitude | +`vm` | [p.u.] | Filtered terminal voltage | +`ip` | [p.u.] | Active-current state | Converter base +`iq` | [p.u.] | Reactive-current state | Converter base +`iqextra` | [p.u.] | HVRCM extra reactive current | Converter base +`il` | [p.u.] | LVPL upper-limit current curve | +`lp` | [p.u./s] | Active-current lower rate bound | +`up` | [p.u./s] | Active-current upper rate bound | diff --git a/GridKit/Model/PhasorDynamics/Converter/REGCA/Regca.cpp b/GridKit/Model/PhasorDynamics/Converter/REGCA/Regca.cpp index b018c4bc3..0f6a2aca6 100644 --- a/GridKit/Model/PhasorDynamics/Converter/REGCA/Regca.cpp +++ b/GridKit/Model/PhasorDynamics/Converter/REGCA/Regca.cpp @@ -12,8 +12,8 @@ namespace GridKit { namespace Converter { - template - int Regca::evaluateJacobian() + template + int Regca::evaluateJacobian() { Log::misc() << "Evaluate Jacobian for Regca..." << std::endl; Log::misc() << "Jacobian evaluation is not implemented!" << std::endl; diff --git a/GridKit/Model/PhasorDynamics/Converter/REGCA/Regca.hpp b/GridKit/Model/PhasorDynamics/Converter/REGCA/Regca.hpp index db10ae6c9..0a91d5b0c 100644 --- a/GridKit/Model/PhasorDynamics/Converter/REGCA/Regca.hpp +++ b/GridKit/Model/PhasorDynamics/Converter/REGCA/Regca.hpp @@ -19,10 +19,10 @@ namespace GridKit { namespace PhasorDynamics { - template + template class BusBase; - template + template class SignalNode; } // namespace PhasorDynamics } // namespace GridKit @@ -59,37 +59,39 @@ namespace GridKit MAXIMUM, }; - template - class Regca : public Component + template + class Regca : public Component { - using Component::gridkit_component_id_; - using Component::alpha_; - using Component::f_; - using Component::h_; - using Component::J_; - 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::time_; - using Component::va_system_base_; - using Component::variable_indices_; - using Component::wb_; - using Component::y_; - using Component::yp_; + using Component::gridkit_component_id_; + using Component::alpha_; + using Component::f_; + using Component::h_; + using Component::J_; + 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::time_; + using Component::va_system_base_; + using Component::variable_indices_; + using Component::wb_; + using Component::y_; + using Component::yp_; public: - using RealT = typename Component::RealT; - using bus_type = BusBase; - using signal_type = SignalNode; - using model_data_type = RegcaData; - using MonitorT = Model::VariableMonitor; - - Regca(bus_type* bus); - Regca(bus_type* bus, const model_data_type& data); + using ScalarT = scalar_type; + using IdxT = index_type; + using RealT = typename Component::RealT; + using BusT = BusBase; + using SignalT = SignalNode; + using ModelDataT = RegcaData; + using MonitorT = Model::VariableMonitor; + + Regca(BusT* bus); + Regca(BusT* bus, const ModelDataT& data); ~Regca(); int setGridKitComponentID(IdxT) override final; @@ -117,7 +119,7 @@ namespace GridKit ScalarT*, ScalarT*, ScalarT*, ScalarT*); private: - void initializeParameters(const model_data_type& data); + void initializeParameters(const ModelDataT& data); void initializeMonitor(); void setDerivedParameters(); @@ -154,7 +156,7 @@ namespace GridKit return bus_->Ii(); } - bus_type* bus_{nullptr}; + BusT* bus_{nullptr}; RealT P0_{0}; RealT Q0_{0}; diff --git a/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaData.hpp b/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaData.hpp index 49d23b861..998049475 100644 --- a/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaData.hpp +++ b/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaData.hpp @@ -63,9 +63,9 @@ namespace GridKit up ///< Active-current upper rate bound }; - template - struct RegcaData : public ComponentData + struct RegcaData : public ComponentData diff --git a/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaDependencyTracking.cpp b/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaDependencyTracking.cpp index 20d688688..bf470cc26 100644 --- a/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaDependencyTracking.cpp +++ b/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaDependencyTracking.cpp @@ -12,8 +12,8 @@ namespace GridKit { namespace Converter { - template - int Regca::evaluateJacobian() + template + int Regca::evaluateJacobian() { Log::misc() << "Evaluate Jacobian for Regca..." << std::endl; Log::misc() << "Jacobian evaluation is not implemented!" << std::endl; diff --git a/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaEnzyme.cpp b/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaEnzyme.cpp index 8bfb167d7..41b4a51b6 100644 --- a/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaEnzyme.cpp +++ b/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaEnzyme.cpp @@ -14,8 +14,8 @@ namespace GridKit { namespace Converter { - template - int Regca::evaluateJacobian() + template + int Regca::evaluateJacobian() { Log::misc() << "Evaluate Jacobian for Regca..." << std::endl; Log::misc() << "Jacobian evaluation is experimental!" << std::endl; diff --git a/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaImpl.hpp b/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaImpl.hpp index 212d29788..0a195e0c9 100644 --- a/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaImpl.hpp @@ -25,15 +25,15 @@ namespace GridKit { using Log = ::GridKit::Utilities::Logger; - template - Regca::Regca(bus_type* bus) + template + Regca::Regca(BusT* bus) : bus_(bus) { size_ = static_cast(RegcaInternalVariables::MAXIMUM); } - template - Regca::Regca(bus_type* bus, const model_data_type& data) + template + Regca::Regca(BusT* bus, const ModelDataT& data) : bus_(bus), monitor_(std::make_unique(data)) { @@ -42,13 +42,13 @@ namespace GridKit size_ = static_cast(RegcaInternalVariables::MAXIMUM); } - template - Regca::~Regca() + template + Regca::~Regca() { } - template - void Regca::setDerivedParameters() + template + void Regca::setDerivedParameters() { Mp_ = static_cast(100.0) * Rpmax_; use_lvpl_ = ZERO; @@ -63,29 +63,32 @@ namespace GridKit va_converter_base_ = mva_base_ * static_cast(1.0e6); } - template - ScalarT Regca::activeCurrentLowerRateBound(ScalarT ip) const + template + scalar_type Regca::activeCurrentLowerRateBound( + scalar_type ip) const { return -Rpmax_ - (Mp_ - Rpmax_) * Math::sigmoid(ip); } - template - ScalarT Regca::activeCurrentUpperRateBound(ScalarT ip, ScalarT il) const + template + scalar_type Regca::activeCurrentUpperRateBound( + scalar_type ip, + scalar_type il) const { const ScalarT sigma_ip = Math::sigmoid(ip); return Mp_ * (ONE - sigma_ip) + Rpmax_ * sigma_ip * (bypass_lvpl_ + use_lvpl_ * Math::sigmoid(il - ip)); } - template - void Regca::initializeParameters(const model_data_type& data) + template + void Regca::initializeParameters(const ModelDataT& data) { - using Params = typename model_data_type::Parameters; - using Ports = typename model_data_type::Ports; + using Params = typename ModelDataT::Parameters; + using Ports = typename ModelDataT::Ports; parameter_error_count_ = 0; - auto loadRequiredReal = [&](auto key, RealT& target, const char* name) + auto load_required_real = [&](auto key, RealT& target, const char* name) { if (!data.parameters.contains(key)) { @@ -110,30 +113,7 @@ namespace GridKit } }; - auto loadOptionalReal = [&](auto key, RealT& target, const char* name) - { - if (!data.parameters.contains(key)) - { - return; - } - - const auto& value = data.parameters.at(key); - if (const auto* real_value = std::get_if(&value)) - { - target = *real_value; - } - else if (const auto* index_value = std::get_if(&value)) - { - target = static_cast(*index_value); - } - else - { - Log::error() << "Regca: parameter '" << name << "' must be numeric\n"; - ++parameter_error_count_; - } - }; - - auto loadRequiredSwitch = [&](auto key, bool& target, const char* name) + auto load_required_switch = [&](auto key, bool& target, const char* name) { if (!data.parameters.contains(key)) { @@ -159,21 +139,21 @@ namespace GridKit } }; - loadOptionalReal(Params::P0, P0_, "P0"); - loadOptionalReal(Params::Q0, Q0_, "Q0"); - loadRequiredReal(Params::mva, mva_base_, "mva"); - loadRequiredReal(Params::Tg, Tg_, "Tg"); - loadRequiredReal(Params::TM, TM_, "TM"); - loadRequiredReal(Params::Rqmax, Rqmax_, "Rqmax"); - loadRequiredReal(Params::Rqmin, Rqmin_, "Rqmin"); - loadRequiredReal(Params::Rpmax, Rpmax_, "Rpmax"); - loadRequiredSwitch(Params::sL, sL_, "sL"); - loadRequiredReal(Params::IL1, IL1_, "IL1"); - loadRequiredReal(Params::VL0, VL0_, "VL0"); - loadRequiredReal(Params::VL1, VL1_, "VL1"); - loadRequiredReal(Params::VA0, VA0_, "VA0"); - loadRequiredReal(Params::VA1, VA1_, "VA1"); - loadRequiredReal(Params::Vhvmax, Vhvmax_, "Vhvmax"); + load_required_real(Params::P0, P0_, "P0"); + load_required_real(Params::Q0, Q0_, "Q0"); + load_required_real(Params::mva, mva_base_, "mva"); + load_required_real(Params::Tg, Tg_, "Tg"); + load_required_real(Params::TM, TM_, "TM"); + load_required_real(Params::Rqmax, Rqmax_, "Rqmax"); + load_required_real(Params::Rqmin, Rqmin_, "Rqmin"); + load_required_real(Params::Rpmax, Rpmax_, "Rpmax"); + load_required_switch(Params::sL, sL_, "sL"); + load_required_real(Params::IL1, IL1_, "IL1"); + load_required_real(Params::VL0, VL0_, "VL0"); + load_required_real(Params::VL1, VL1_, "VL1"); + load_required_real(Params::VA0, VA0_, "VA0"); + load_required_real(Params::VA1, VA1_, "VA1"); + load_required_real(Params::Vhvmax, Vhvmax_, "Vhvmax"); if (data.ports.contains(Ports::bus)) { @@ -183,16 +163,16 @@ namespace GridKit setDerivedParameters(); } - template - const Model::VariableMonitorBase* Regca::getMonitor() const + template + const Model::VariableMonitorBase* Regca::getMonitor() const { return monitor_.get(); } - template - void Regca::initializeMonitor() + template + void Regca::initializeMonitor() { - using Variable = typename model_data_type::MonitorableVariables; + using Variable = typename ModelDataT::MonitorableVariables; auto index = [](RegcaInternalVariables variable) { return static_cast(variable); @@ -224,15 +204,15 @@ namespace GridKit { return y_[index(RegcaInternalVariables::UP)]; }); } - template - int Regca::setGridKitComponentID(IdxT component_id) + template + int Regca::setGridKitComponentID(IdxT component_id) { gridkit_component_id_ = component_id; return 0; } - template - int Regca::allocate() + template + int Regca::allocate() { size_ = static_cast(RegcaInternalVariables::MAXIMUM); auto size = static_cast(size_); @@ -288,8 +268,8 @@ namespace GridKit return 0; } - template - int Regca::verify() const + template + int Regca::verify() const { int ret = static_cast(parameter_error_count_); @@ -339,8 +319,8 @@ namespace GridKit return ret; } - template - int Regca::initialize() + template + int Regca::initialize() { if (bus_ == nullptr) { @@ -371,67 +351,62 @@ namespace GridKit const auto PBR = static_cast(RegcaInternalVariables::PBR); const auto QBR = static_cast(RegcaInternalVariables::QBR); - const ScalarT vr = Vr(); - const ScalarT vi = Vi(); - const ScalarT vt2 = vr * vr + vi * vi; - const ScalarT vt = std::sqrt(vt2); + const ScalarT vr = Vr(); + const ScalarT vi = Vi(); + const ScalarT vt = std::sqrt(vr * vr + vi * vi); if (vt <= ZERO) { Log::error() << "Regca: terminal voltage magnitude must be positive at initialization\n"; return 1; } - if (vt >= Vhvmax_) { - Log::error() << "Regca: terminal voltage magnitude must be less than Vhvmax at initialization\n"; + Log::error() + << "Regca: terminal voltage magnitude must be below Vhvmax at initialization\n"; return 1; } + // REGCA owns the network terminal and establishes the initial + // converter operating point from the power-flow injection. Controller + // command ports may not have been initialized yet, so initialization + // resolves the commands from P0/Q0 and publishes them to attached + // ports below. P0/Q0 are given on the system base; toComponentBase + // converts them to the converter base the internal states use. const ScalarT lvacm = Math::linseg(vt, VA0_, VA1_, ONE); - ScalarT ipcmd0{ZERO}; - ScalarT iqcmd0{ZERO}; - if (signals_.template isAttached()) - { - ipcmd0 = signals_.template readExternalVariable(); - } - else if (P0_ != ZERO) + if (P0_ != ZERO && (vt <= VA0_ || lvacm <= ZERO) ) { - if (vt <= VA0_ || lvacm <= ZERO) - { - Log::error() << "Regca: LVACM gain is zero with nonzero initial active power\n"; - return 1; - } - - ipcmd0 = toComponentBase(static_cast(P0_) / vt) / lvacm; + Log::error() << "Regca: LVACM gain is zero with nonzero initial active power\n"; + return 1; } - if (signals_.template isAttached()) + ScalarT ipcmd0{ZERO}; + if (P0_ != ZERO) { - iqcmd0 = signals_.template readExternalVariable(); - } - else if (Q0_ != ZERO) - { - iqcmd0 = toComponentBase(static_cast(Q0_) / vt); + ipcmd0 = toComponentBase(static_cast(P0_) / vt) / lvacm; } + const ScalarT iqcmd0 = toComponentBase(static_cast(Q0_) / vt); const ScalarT iqextra0{ZERO}; const ScalarT qnet0 = iqcmd0 - iqextra0; y_[VM] = vt; - y_[IQ] = iqcmd0; - y_[IP] = ipcmd0; y_[VT] = vt; - y_[II] = (lvacm * y_[IP] * vi - qnet0 * vr) / vt; + y_[IP] = ipcmd0; + y_[IQ] = iqcmd0; y_[IQEXTRA] = iqextra0; y_[IL] = Math::linseg(vt, VL0_, VL1_, IL1_); - y_[IR] = (lvacm * y_[IP] * vr + qnet0 * vi) / vt; + y_[IR] = (qnet0 * vi + lvacm * ipcmd0 * vr) / vt; + y_[II] = (-qnet0 * vr + lvacm * ipcmd0 * vi) / vt; y_[LP] = activeCurrentLowerRateBound(y_[IP]); y_[UP] = activeCurrentUpperRateBound(y_[IP], y_[IL]); y_[PBR] = toSystemBase(vr * y_[IR] + vi * y_[II]); y_[QBR] = toSystemBase(vi * y_[IR] - vr * y_[II]); + // Retain the resolved commands as the constant source used during + // residual evaluation when no controller drives the command ports, and + // select the reactive-current rate-limit branch from the command sign. ipcmd_set_ = ipcmd0; iqcmd_set_ = iqcmd0; iq_use_upper_ = ZERO; @@ -442,12 +417,24 @@ namespace GridKit iq_use_lower_ = ZERO; } + // Seed attached command nodes with the steady-state values. Controller + // initialization can use these signal values, and unattached ports fall + // back to the constants stored above. + if (signals_.template isAttached()) + { + signals_.template writeExternalVariable(ipcmd0); + } + if (signals_.template isAttached()) + { + signals_.template writeExternalVariable(iqcmd0); + } + std::fill(yp_.begin(), yp_.end(), ZERO); return 0; } - template - int Regca::tagDifferentiable() + template + int Regca::tagDifferentiable() { std::fill(tag_.begin(), tag_.end(), false); tag_[static_cast(RegcaInternalVariables::VM)] = true; @@ -456,8 +443,9 @@ namespace GridKit return 0; } - template - __attribute__((always_inline)) inline int Regca::evaluateInternalResidual( + template + __attribute__((always_inline)) inline int + Regca::evaluateInternalResidual( ScalarT* y, ScalarT* yp, ScalarT* wb, @@ -503,26 +491,27 @@ namespace GridKit const ScalarT ipcmd = ws[IPCMD]; const ScalarT iqcmd = ws[IQCMD]; - const ScalarT fq = (iqcmd - iq) / Tg_; - const ScalarT fp = (ipcmd - ip) / Tg_; + const ScalarT iq_error = iqcmd - iq; + const ScalarT ip_error = ipcmd - ip; - const ScalarT iq_rate = - iq_use_upper_ * (fq - Math::ramp(fq - Rqmax_)) - + iq_use_lower_ * (fq + Math::ramp(Rqmin_ - fq)); + const ScalarT iq_step = + iq_use_upper_ * (iq_error - Math::ramp(iq_error - Tg_ * Rqmax_)) + + iq_use_lower_ * (iq_error + Math::ramp(Tg_ * Rqmin_ - iq_error)); - const ScalarT ip_rate = lp + Math::ramp(fp - lp) - Math::ramp(fp - up); + const ScalarT ip_step = + Tg_ * lp + Math::ramp(ip_error - Tg_ * lp) - Math::ramp(ip_error - Tg_ * up); - f[VM] = -vm_dot + (vt - vm) / TM_; - f[IQ] = -iq_dot + iq_rate; - f[IP] = -ip_dot + ip_rate; + f[VM] = -TM_ * vm_dot - vm + vt; + f[IQ] = -Tg_ * iq_dot + iq_step; + f[IP] = -Tg_ * ip_dot + ip_step; const ScalarT lvacm = Math::linseg(vt, VA0_, VA1_, ONE); const ScalarT qnet = iq - iqextra; f[VT] = -vt * vt + vr * vr + vi * vi; - f[II] = -vt * ii + lvacm * ip * vi - qnet * vr; - f[IQEXTRA] = -iqextra + Math::ramp(iqextra - (Vhvmax_ - vt)); + f[II] = -vt * ii - qnet * vr + lvacm * ip * vi; + f[IQEXTRA] = -iqextra + Math::ramp(vt - Vhvmax_); f[IL] = -il + Math::linseg(vm, VL0_, VL1_, IL1_); - f[IR] = -vt * ir + lvacm * ip * vr + qnet * vi; + f[IR] = -vt * ir + qnet * vi + lvacm * ip * vr; f[LP] = -lp + activeCurrentLowerRateBound(ip); f[UP] = -up + activeCurrentUpperRateBound(ip, il); f[PBR] = -pbr + toSystemBase(vr * ir + vi * ii); @@ -531,8 +520,8 @@ namespace GridKit return 0; } - template - __attribute__((always_inline)) inline int Regca::evaluateBusResidual( + template + __attribute__((always_inline)) inline int Regca::evaluateBusResidual( ScalarT* y, [[maybe_unused]] ScalarT* yp, [[maybe_unused]] ScalarT* wb, @@ -546,8 +535,8 @@ namespace GridKit return 0; } - template - int Regca::evaluateResidual() + template + int Regca::evaluateResidual() { const auto IPCMD = static_cast(RegcaExternalVariables::IPCMD); const auto IQCMD = static_cast(RegcaExternalVariables::IQCMD); diff --git a/tests/UnitTests/PhasorDynamics/CMakeLists.txt b/tests/UnitTests/PhasorDynamics/CMakeLists.txt index ed5f171a3..b09a2f450 100644 --- a/tests/UnitTests/PhasorDynamics/CMakeLists.txt +++ b/tests/UnitTests/PhasorDynamics/CMakeLists.txt @@ -98,8 +98,8 @@ add_executable(test_phasor_converter_regca runConverterRegcaTests.cpp) target_link_libraries( test_phasor_converter_regca GridKit::definitions - GridKit::phasor_dynamics_components - GridKit::phasor_dynamics_components_dependency_tracking + GridKit::phasor_dynamics_systemmodel + GridKit::phasor_dynamics_systemmodel_dependency_tracking GridKit::testing) add_executable(test_phasor_stabilizer_ieeest runStabilizerIeeestTests.cpp) diff --git a/tests/UnitTests/PhasorDynamics/ConverterRegcaTests.hpp b/tests/UnitTests/PhasorDynamics/ConverterRegcaTests.hpp index 94fc10574..7ef9e971e 100644 --- a/tests/UnitTests/PhasorDynamics/ConverterRegcaTests.hpp +++ b/tests/UnitTests/PhasorDynamics/ConverterRegcaTests.hpp @@ -1,13 +1,13 @@ #pragma once #include +#include #include #include #include #include #include -#include #include #include #include @@ -22,33 +22,34 @@ namespace GridKit { namespace Testing { - template + template class ConverterRegcaTests { public: - using RealT = typename PhasorDynamics::Component::RealT; + using ScalarT = scalar_type; + using IdxT = index_type; + using RealT = typename PhasorDynamics::Component::RealT; ConverterRegcaTests() = default; ~ConverterRegcaTests() = default; - static constexpr ScalarT kTol = static_cast(1.0e-10); + static constexpr ScalarT kTol = static_cast(1.0e-9); - TestOutcome constructor() + TestOutcome constructionAndValidation() { TestStatus success = true; PhasorDynamics::Bus bus(1.0, 0.0); PhasorDynamics::Converter::Regca minimal(&bus); - success *= (minimal.size() == static_cast(PhasorDynamics::Converter::RegcaInternalVariables::MAXIMUM)); + success *= (minimal.size() == static_cast(Vars::MAXIMUM)); success *= (minimal.getMonitor() == nullptr); success *= (minimal.verify() > 0); - auto data = makeTestData(); - PhasorDynamics::Converter::Regca from_data(&bus, data); - success *= (from_data.size() == static_cast(PhasorDynamics::Converter::RegcaInternalVariables::MAXIMUM)); - success *= (from_data.getMonitor() != nullptr); - success *= (from_data.verify() == 0); + PhasorDynamics::Converter::Regca configured(&bus, makeData()); + success *= (configured.size() == static_cast(Vars::MAXIMUM)); + success *= (configured.getMonitor() != nullptr); + success *= (configured.verify() == 0); return success.report(__func__); } @@ -59,12 +60,22 @@ namespace GridKit PhasorDynamics::Bus bus(1.0, 0.0); - auto missing = makeTestData(); + auto missing = makeData(); missing.parameters.erase(Params::Tg); PhasorDynamics::Converter::Regca missing_model(&bus, missing); success *= (missing_model.verify() > 0); - auto bad_switch = makeTestData(); + auto missing_p0 = makeData(); + missing_p0.parameters.erase(Params::P0); + PhasorDynamics::Converter::Regca missing_p0_model(&bus, missing_p0); + success *= (missing_p0_model.verify() > 0); + + auto missing_q0 = makeData(); + missing_q0.parameters.erase(Params::Q0); + PhasorDynamics::Converter::Regca missing_q0_model(&bus, missing_q0); + success *= (missing_q0_model.verify() > 0); + + auto bad_switch = makeData(); bad_switch.parameters[Params::sL] = static_cast(2); PhasorDynamics::Converter::Regca bad_switch_model(&bus, bad_switch); success *= (bad_switch_model.verify() > 0); @@ -82,201 +93,213 @@ namespace GridKit return success.report(__func__); } - TestOutcome lifecycle() + TestOutcome initializesFromPowerFlowAndPublishesSignals() { TestStatus success = true; - PhasorDynamics::Bus bus(1.0, 0.0); + const ScalarT vr{0.8}; + const ScalarT vi{0.6}; + const ScalarT p0{0.4}; + const ScalarT q0{-0.1}; + const RealT mva{50.0}; + + PhasorDynamics::Bus bus(vr, vi); bus.allocate(); bus.initialize(); - bus.evaluateResidual(); - auto data = makeTestData(); + auto data = makeData(); + data.parameters[Params::mva] = mva; + data.parameters[Params::P0] = static_cast(p0); + data.parameters[Params::Q0] = static_cast(q0); + PhasorDynamics::Converter::Regca regca(&bus, data); - success *= (regca.allocate() == 0); - success *= (regca.initialize() == 0); - success *= (regca.tagDifferentiable() == 0); - success *= (regca.evaluateResidual() == 0); - success *= (regca.evaluateJacobian() == 0); - for (size_t i = 0; i < static_cast(regca.size()); ++i) - { - success *= isEqual(regca.yp()[i], static_cast(0.0), kTol); - success *= isEqual(regca.getResidual()[i], static_cast(0.0), kTol); - } + ScalarT ipcmd_value{-1.0}; + ScalarT iqcmd_value{1.0}; + IdxT ipcmd_index = 20; + IdxT iqcmd_index = 21; - using Vars = PhasorDynamics::Converter::RegcaInternalVariables; - for (size_t i = 0; i < static_cast(Vars::MAXIMUM); ++i) - { - const bool expected = (i == static_cast(Vars::VM)) - || (i == static_cast(Vars::IQ)) - || (i == static_cast(Vars::IP)); - success *= (regca.tag()[i] == expected); - } + PhasorDynamics::SignalNode ipcmd_node; + PhasorDynamics::SignalNode iqcmd_node; + PhasorDynamics::SignalNode ir_node; + PhasorDynamics::SignalNode ii_node; + PhasorDynamics::SignalNode pbranch_node; + PhasorDynamics::SignalNode qbranch_node; + ipcmd_node.set(&ipcmd_value, &ipcmd_index); + iqcmd_node.set(&iqcmd_value, &iqcmd_index); - success *= isEqual(bus.Ir(), static_cast(0.0), kTol); - success *= isEqual(bus.Ii(), static_cast(0.0), kTol); + regca.getSignals().template attachSignalNode(&ipcmd_node); + regca.getSignals().template attachSignalNode(&iqcmd_node); + regca.getSignals().template assignSignalNode(&ir_node); + regca.getSignals().template assignSignalNode(&ii_node); + regca.getSignals().template assignSignalNode(&pbranch_node); + regca.getSignals().template assignSignalNode(&qbranch_node); + success *= (regca.allocate() == 0); + success *= (regca.verify() == 0); + success *= (regca.initialize() == 0); + success *= (regca.evaluateResidual() == 0); + + const ScalarT vt = std::sqrt(vr * vr + vi * vi); + const ScalarT lvacm = + Math::linseg(vt, static_cast(0.4), static_cast(0.9), ONE); + const ScalarT ipcmd0 = toComponentBase(p0 / vt, mva) / lvacm; + const ScalarT iqcmd0 = toComponentBase(q0 / vt, mva); + const ScalarT qnet0 = iqcmd0; + const ScalarT ir0 = (lvacm * ipcmd0 * vr + qnet0 * vi) / vt; + const ScalarT ii0 = (lvacm * ipcmd0 * vi - qnet0 * vr) / vt; + + success *= scalarMatches(regca.y()[index(Vars::VM)], vt, "VM"); + success *= scalarMatches(regca.y()[index(Vars::VT)], vt, "VT"); + success *= scalarMatches(regca.y()[index(Vars::IP)], ipcmd0, "IP"); + success *= scalarMatches(regca.y()[index(Vars::IQ)], iqcmd0, "IQ"); + success *= scalarMatches(regca.y()[index(Vars::IR)], ir0, "IR"); + success *= scalarMatches(regca.y()[index(Vars::II)], ii0, "II"); + success *= scalarMatches(regca.y()[index(Vars::PBR)], p0, "PBR"); + success *= scalarMatches(regca.y()[index(Vars::QBR)], q0, "QBR"); + success *= scalarMatches(ipcmd_node.read(), ipcmd0, "ipcmd signal"); + success *= scalarMatches(iqcmd_node.read(), iqcmd0, "iqcmd signal"); + success *= scalarMatches(ir_node.read(), ir0, "ir signal"); + success *= scalarMatches(ii_node.read(), ii0, "ii signal"); + success *= scalarMatches(pbranch_node.read(), p0, "pbranch signal"); + success *= scalarMatches(qbranch_node.read(), q0, "qbranch signal"); + + success *= allResidualsZero(regca); return success.report(__func__); } - TestOutcome steadyStateInitializationGolden() + TestOutcome unconnectedCommandsRemainConstant() { TestStatus success = true; - PhasorDynamics::Bus bus(0.8, 0.6); + auto data = makeData(); + data.parameters[Params::P0] = static_cast(0.6); + data.parameters[Params::Q0] = static_cast(0.2); + + PhasorDynamics::Bus bus(1.0, 0.0); bus.allocate(); bus.initialize(); - auto data = makeGoldenTestData(true); - PhasorDynamics::Converter::Regca regca(&bus, data); - PhasorDynamics::SignalNode ipcmd_node; - PhasorDynamics::SignalNode iqcmd_node; - ScalarT ipcmd_value{0.92}; - ScalarT iqcmd_value{-0.44}; - IdxT ipcmd_index = 21; - IdxT iqcmd_index = 22; - ipcmd_node.set(&ipcmd_value, &ipcmd_index); - iqcmd_node.set(&iqcmd_value, &iqcmd_index); - regca.getSignals().template attachSignalNode(&ipcmd_node); - regca.getSignals().template attachSignalNode(&iqcmd_node); - - regca.allocate(); + success *= (regca.allocate() == 0); + success *= (regca.verify() == 0); success *= (regca.initialize() == 0); + success *= (regca.evaluateResidual() == 0); - std::vector expected_y(static_cast(Vars::MAXIMUM), static_cast(0.0)); - expected_y[index(Vars::VM)] = static_cast(1.0); - expected_y[index(Vars::IQ)] = static_cast(-0.44); - expected_y[index(Vars::IP)] = static_cast(0.92); - expected_y[index(Vars::VT)] = static_cast(1.0); - expected_y[index(Vars::II)] = static_cast(0.904); - expected_y[index(Vars::IQEXTRA)] = static_cast(0.0); - expected_y[index(Vars::IL)] = static_cast(1.1); - expected_y[index(Vars::IR)] = static_cast(0.472); - expected_y[index(Vars::LP)] = static_cast(-70.0); - expected_y[index(Vars::UP)] = static_cast(0.7); - expected_y[index(Vars::PBR)] = static_cast(0.92); - expected_y[index(Vars::QBR)] = static_cast(-0.44); - - success *= vectorMatches(regca.y(), expected_y, "REGCA initialization state"); - for (size_t i = 0; i < static_cast(regca.size()); ++i) - { - success *= isEqual(regca.yp()[i], static_cast(0.0), kTol); - } + success *= isEqual(regca.y()[index(Vars::IP)], static_cast(0.6), kTol); + success *= isEqual(regca.y()[index(Vars::IQ)], static_cast(0.2), kTol); + success *= allResidualsZero(regca); return success.report(__func__); } - TestOutcome attachedSignalInitialization() + TestOutcome externalCommandsDriveRuntimeResidual() { TestStatus success = true; + auto data = makeData(); + data.parameters[Params::P0] = static_cast(0.5); + data.parameters[Params::Q0] = static_cast(-0.1); + PhasorDynamics::Bus bus(1.0, 0.0); bus.allocate(); bus.initialize(); - auto data = makeTestData(); - PhasorDynamics::Converter::Regca regca(&bus, data); + ScalarT ipcmd_value{0.0}; + ScalarT iqcmd_value{0.0}; + IdxT ipcmd_index = 22; + IdxT iqcmd_index = 23; PhasorDynamics::SignalNode ipcmd_node; PhasorDynamics::SignalNode iqcmd_node; - ScalarT ipcmd_value{0.75}; - ScalarT iqcmd_value{-0.20}; - IdxT ipcmd_index = 21; - IdxT iqcmd_index = 22; - ipcmd_node.set(&ipcmd_value, &ipcmd_index); iqcmd_node.set(&iqcmd_value, &iqcmd_index); + + PhasorDynamics::Converter::Regca regca(&bus, data); regca.getSignals().template attachSignalNode(&ipcmd_node); regca.getSignals().template attachSignalNode(&iqcmd_node); - regca.allocate(); + success *= (regca.allocate() == 0); + success *= (regca.verify() == 0); success *= (regca.initialize() == 0); - success *= isEqual(ipcmd_value, static_cast(0.75), kTol); - success *= isEqual(iqcmd_value, static_cast(-0.20), kTol); - success *= isEqual(regca.y()[index(Vars::IP)], static_cast(0.75), kTol); - success *= isEqual(regca.y()[index(Vars::IQ)], static_cast(-0.20), kTol); + ipcmd_value = static_cast(0.55); + iqcmd_value = static_cast(-0.05); + bus.evaluateResidual(); + success *= (regca.evaluateResidual() == 0); + + const ScalarT tg = static_cast(0.02); + const ScalarT ip_error = ipcmd_value - regca.y()[index(Vars::IP)]; + const ScalarT iq_error = iqcmd_value - regca.y()[index(Vars::IQ)]; + + const ScalarT expected_ip_residual = + tg * regca.y()[index(Vars::LP)] + + Math::ramp(ip_error - tg * regca.y()[index(Vars::LP)]) + - Math::ramp(ip_error - tg * regca.y()[index(Vars::UP)]); + const ScalarT expected_iq_residual = + iq_error + Math::ramp(tg * static_cast(-999.0) - iq_error); + + success *= isEqual(regca.getResidual()[index(Vars::IP)], expected_ip_residual, kTol); + success *= isEqual(regca.getResidual()[index(Vars::IQ)], expected_iq_residual, kTol); + success *= (regca.getResidual()[index(Vars::IP)] > ZERO); + success *= (regca.getResidual()[index(Vars::IQ)] > ZERO); return success.report(__func__); } - TestOutcome fallbackInit() + TestOutcome invalidInitialization() { TestStatus success = true; - PhasorDynamics::Bus bus(0.8, 0.6); - bus.allocate(); - bus.initialize(); + auto data = makeData(); - auto data = makeTestData(); - data.parameters[Params::mva] = static_cast(50.0); - data.parameters[Params::IL1] = static_cast(2.0); - data.parameters[Params::P0] = static_cast(0.8); - data.parameters[Params::Q0] = static_cast(-0.1); + { + PhasorDynamics::Bus bus(1.3, 0.0); + bus.allocate(); + bus.initialize(); - PhasorDynamics::Converter::Regca regca(&bus, data); - regca.allocate(); - success *= (regca.initialize() == 0); + PhasorDynamics::Converter::Regca regca(&bus, data); + success *= (regca.allocate() == 0); + success *= (regca.initialize() > 0); + } - success *= isEqual(regca.y()[index(Vars::IR)], static_cast(1.16), kTol); - success *= isEqual(regca.y()[index(Vars::II)], static_cast(1.12), kTol); - success *= isEqual(regca.y()[index(Vars::PBR)], static_cast(0.8), kTol); - success *= isEqual(regca.y()[index(Vars::QBR)], static_cast(-0.1), kTol); - success *= isEqual(regca.y()[index(Vars::IP)], static_cast(1.6), kTol); - success *= isEqual(regca.y()[index(Vars::IQ)], static_cast(-0.2), kTol); + { + PhasorDynamics::Bus bus(0.0, 0.0); + bus.allocate(); + bus.initialize(); - bus.evaluateResidual(); - success *= (regca.evaluateResidual() == 0); + PhasorDynamics::Converter::Regca regca(&bus, data); + success *= (regca.allocate() == 0); + success *= (regca.initialize() > 0); + } - for (size_t i = 0; i < static_cast(regca.size()); ++i) { - success *= isEqual(regca.getResidual()[i], static_cast(0.0), kTol); + auto low_voltage_data = data; + low_voltage_data.parameters[Params::P0] = static_cast(0.5); + + PhasorDynamics::Bus bus(0.2, 0.0); + bus.allocate(); + bus.initialize(); + + PhasorDynamics::Converter::Regca regca(&bus, low_voltage_data); + success *= (regca.allocate() == 0); + success *= (regca.initialize() > 0); } - success *= isEqual(bus.Ir(), static_cast(0.58), kTol); - success *= isEqual(bus.Ii(), static_cast(0.56), kTol); - return success.report(__func__); - } + { + auto low_voltage_data = data; + low_voltage_data.parameters[Params::P0] = static_cast(0.0); + low_voltage_data.parameters[Params::Q0] = static_cast(0.1); - TestOutcome invalidInitialization() - { - TestStatus success = true; + PhasorDynamics::Bus bus(0.2, 0.0); + bus.allocate(); + bus.initialize(); - auto data = makeTestData(); - - PhasorDynamics::Bus high_voltage_bus(1.3, 0.0); - high_voltage_bus.allocate(); - high_voltage_bus.initialize(); - PhasorDynamics::Converter::Regca high_voltage_regca(&high_voltage_bus, data); - high_voltage_regca.allocate(); - success *= (high_voltage_regca.initialize() > 0); - - auto low_voltage_fallback_data = data; - low_voltage_fallback_data.parameters[Params::P0] = static_cast(0.75); - PhasorDynamics::Bus low_voltage_fallback_bus(0.2, 0.0); - low_voltage_fallback_bus.allocate(); - low_voltage_fallback_bus.initialize(); - PhasorDynamics::Converter::Regca low_voltage_fallback_regca( - &low_voltage_fallback_bus, low_voltage_fallback_data); - low_voltage_fallback_regca.allocate(); - success *= (low_voltage_fallback_regca.initialize() > 0); - - PhasorDynamics::Bus low_voltage_attached_bus(0.2, 0.0); - low_voltage_attached_bus.allocate(); - low_voltage_attached_bus.initialize(); - PhasorDynamics::Converter::Regca low_voltage_attached_regca( - &low_voltage_attached_bus, low_voltage_fallback_data); - PhasorDynamics::SignalNode ipcmd_node; - ScalarT ipcmd_value{0.75}; - IdxT ipcmd_index = 21; - ipcmd_node.set(&ipcmd_value, &ipcmd_index); - low_voltage_attached_regca.getSignals().template attachSignalNode( - &ipcmd_node); - low_voltage_attached_regca.allocate(); - success *= (low_voltage_attached_regca.initialize() == 0); + PhasorDynamics::Converter::Regca regca(&bus, low_voltage_data); + success *= (regca.allocate() == 0); + success *= (regca.initialize() == 0); + } return success.report(__func__); } @@ -286,14 +309,14 @@ namespace GridKit TestStatus success = true; PhasorDynamics::Bus bus(1.0, 0.0); - PhasorDynamics::Converter::Regca regca(&bus, makeTestData()); + PhasorDynamics::Converter::Regca regca(&bus, makeData()); PhasorDynamics::SignalNode ipcmd_node; PhasorDynamics::SignalNode iqcmd_node; ScalarT ipcmd_value{0.25}; ScalarT iqcmd_value{-0.10}; - IdxT ipcmd_index = 21; - IdxT iqcmd_index = 22; + IdxT ipcmd_index = 24; + IdxT iqcmd_index = 25; regca.getSignals().template attachSignalNode(&ipcmd_node); success *= (regca.verify() > 0); @@ -314,88 +337,172 @@ namespace GridKit { TestStatus success = true; - PhasorDynamics::Converter::Regca regca(nullptr, makeTestData()); + PhasorDynamics::Converter::Regca regca(nullptr, makeData()); success *= (regca.verify() > 0); + success *= (regca.initialize() > 0); return success.report(__func__); } - TestOutcome residualGoldenVectors() + TestOutcome busInjectionUsesSystemBase() { TestStatus success = true; - std::vector positive_q_lvpl(static_cast(Vars::MAXIMUM), static_cast(0.0)); - positive_q_lvpl[index(Vars::VM)] = static_cast(0.29); - positive_q_lvpl[index(Vars::IQ)] = static_cast(0.52); - positive_q_lvpl[index(Vars::IP)] = static_cast(0.21999997439919405); - positive_q_lvpl[index(Vars::VT)] = static_cast(0.0046000000000000485); - positive_q_lvpl[index(Vars::II)] = static_cast(0.2546); - positive_q_lvpl[index(Vars::IQEXTRA)] = static_cast(-0.03); - positive_q_lvpl[index(Vars::IL)] = static_cast(0.29199937917427254); - positive_q_lvpl[index(Vars::IR)] = static_cast(0.26); - positive_q_lvpl[index(Vars::LP)] = static_cast(-69.6); - positive_q_lvpl[index(Vars::UP)] = static_cast(-0.2999999999999803); - positive_q_lvpl[index(Vars::PBR)] = static_cast(0.0); - positive_q_lvpl[index(Vars::QBR)] = static_cast(0.0); - - std::vector nonpositive_q_no_lvpl(static_cast(Vars::MAXIMUM), static_cast(0.0)); - nonpositive_q_no_lvpl[index(Vars::VM)] = static_cast(0.29); - nonpositive_q_no_lvpl[index(Vars::IQ)] = static_cast(1.5200000000000002); - nonpositive_q_no_lvpl[index(Vars::IP)] = static_cast(0.21999997439919405); - nonpositive_q_no_lvpl[index(Vars::VT)] = static_cast(0.0046000000000000485); - nonpositive_q_no_lvpl[index(Vars::II)] = static_cast(0.2546); - nonpositive_q_no_lvpl[index(Vars::IQEXTRA)] = static_cast(-0.03); - nonpositive_q_no_lvpl[index(Vars::IL)] = static_cast(0.29199937917427254); - nonpositive_q_no_lvpl[index(Vars::IR)] = static_cast(0.26); - nonpositive_q_no_lvpl[index(Vars::LP)] = static_cast(-69.6); - nonpositive_q_no_lvpl[index(Vars::UP)] = static_cast(0.39999999999999997); - nonpositive_q_no_lvpl[index(Vars::PBR)] = static_cast(0.0); - nonpositive_q_no_lvpl[index(Vars::QBR)] = static_cast(0.0); - - success *= residualGoldenVectorCase(static_cast(0.1), true, positive_q_lvpl); - success *= residualGoldenVectorCase(static_cast(-0.1), false, nonpositive_q_no_lvpl); + const ScalarT vr{0.8}; + const ScalarT vi{0.6}; + const ScalarT p0{0.4}; + const ScalarT q0{-0.1}; + const RealT mva{50.0}; + + PhasorDynamics::Bus bus(vr, vi); + bus.allocate(); + bus.initialize(); + + auto data = makeData(); + data.parameters[Params::mva] = mva; + data.parameters[Params::P0] = static_cast(p0); + data.parameters[Params::Q0] = static_cast(q0); + + PhasorDynamics::Converter::Regca regca(&bus, data); + success *= (regca.allocate() == 0); + success *= (regca.initialize() == 0); + + bus.evaluateResidual(); + success *= (regca.evaluateResidual() == 0); + + success *= isEqual(bus.Ir(), toSystemBase(regca.y()[index(Vars::IR)], mva), kTol); + success *= isEqual(bus.Ii(), toSystemBase(regca.y()[index(Vars::II)], mva), kTol); + success *= isEqual(bus.Ir(), static_cast(0.26), kTol); + success *= isEqual(bus.Ii(), static_cast(0.32), kTol); return success.report(__func__); } - TestOutcome busInjection() + TestOutcome residualEquations() { TestStatus success = true; - PhasorDynamics::Bus bus(0.8, 0.6); + PhasorDynamics::Bus bus(0.95, 0.25); bus.allocate(); bus.initialize(); - auto data = makeTestData(); - data.parameters[Params::mva] = static_cast(50.0); - - PhasorDynamics::Converter::Regca regca(&bus, data); + ScalarT ipcmd_value{0.9}; + ScalarT iqcmd_value{0.1}; + IdxT ipcmd_index = 26; + IdxT iqcmd_index = 27; PhasorDynamics::SignalNode ipcmd_node; PhasorDynamics::SignalNode iqcmd_node; - PhasorDynamics::SignalNode pbranch_node; - PhasorDynamics::SignalNode qbranch_node; - ScalarT ipcmd_value{1.6}; - ScalarT iqcmd_value{-0.2}; - IdxT ipcmd_index = 21; - IdxT iqcmd_index = 22; ipcmd_node.set(&ipcmd_value, &ipcmd_index); iqcmd_node.set(&iqcmd_value, &iqcmd_index); + + PhasorDynamics::Converter::Regca regca(&bus, makeDynamicData()); regca.getSignals().template attachSignalNode(&ipcmd_node); regca.getSignals().template attachSignalNode(&iqcmd_node); - regca.getSignals().template assignSignalNode(&pbranch_node); - regca.getSignals().template assignSignalNode(&qbranch_node); - regca.allocate(); + success *= (regca.allocate() == 0); success *= (regca.initialize() == 0); + setResidualState(regca); bus.evaluateResidual(); success *= (regca.evaluateResidual() == 0); - success *= isEqual(bus.Ir(), static_cast(0.58), kTol); - success *= isEqual(bus.Ii(), static_cast(0.56), kTol); - success *= isEqual(pbranch_node.read(), static_cast(0.8), kTol); - success *= isEqual(qbranch_node.read(), static_cast(-0.1), kTol); + const auto expected = expectedResidualForState(ipcmd_value, iqcmd_value); + success *= vectorMatches(regca.getResidual(), expected, "REGCA residual"); + + return success.report(__func__); + } + + TestOutcome highVoltageReactiveCurrentRoot() + { + TestStatus success = true; + + auto data = makeDynamicData(); + + PhasorDynamics::Bus bus(1.0, 0.0); + bus.allocate(); + bus.initialize(); + + PhasorDynamics::Converter::Regca regca(&bus, data); + success *= (regca.allocate() == 0); + success *= (regca.initialize() == 0); + + regca.y()[index(Vars::VT)] = static_cast(1.35); + regca.y()[index(Vars::IQEXTRA)] = + Math::ramp(regca.y()[index(Vars::VT)] - static_cast(1.3)); + + bus.evaluateResidual(); + success *= (regca.evaluateResidual() == 0); + success *= scalarMatches(regca.getResidual()[index(Vars::IQEXTRA)], + static_cast(0.0), + "HVRCM residual root"); + + return success.report(__func__); + } + + TestOutcome limiterBranchCoverage() + { + TestStatus success = true; + + { + auto data = makeDynamicData(); + data.parameters[Params::Q0] = static_cast(0.2); + + PhasorDynamics::Bus bus(1.0, 0.0); + bus.allocate(); + bus.initialize(); + + ScalarT iqcmd_value{0.2}; + IdxT iqcmd_index = 28; + PhasorDynamics::SignalNode iqcmd_node; + iqcmd_node.set(&iqcmd_value, &iqcmd_index); + + PhasorDynamics::Converter::Regca regca(&bus, data); + regca.getSignals().template attachSignalNode(&iqcmd_node); + + success *= (regca.allocate() == 0); + success *= (regca.verify() == 0); + success *= (regca.initialize() == 0); + + iqcmd_value = static_cast(0.4); + bus.evaluateResidual(); + success *= (regca.evaluateResidual() == 0); + + const ScalarT iq_error = iqcmd_value - regca.y()[index(Vars::IQ)]; + const ScalarT expected = + iq_error - Math::ramp(iq_error - static_cast(0.2) * static_cast(0.5)); + + success *= scalarMatches(regca.getResidual()[index(Vars::IQ)], + expected, + "positive IQ branch"); + } + + { + auto data = makeDynamicData(); + data.parameters[Params::sL] = false; + + PhasorDynamics::Bus bus(1.0, 0.0); + bus.allocate(); + bus.initialize(); + + PhasorDynamics::Converter::Regca regca(&bus, data); + success *= (regca.allocate() == 0); + success *= (regca.initialize() == 0); + + bus.evaluateResidual(); + success *= (regca.evaluateResidual() == 0); + + const ScalarT ip = regca.y()[index(Vars::IP)]; + const ScalarT sigma_ip = Math::sigmoid(ip); + const RealT rpmax = static_cast(0.7); + const RealT mp = static_cast(100.0) * rpmax; + const ScalarT expected = mp * (ONE - sigma_ip) + rpmax * sigma_ip; + + success *= scalarMatches(regca.y()[index(Vars::UP)], expected, "UP without LVPL"); + success *= scalarMatches(regca.getResidual()[index(Vars::UP)], + static_cast(0.0), + "UP residual without LVPL"); + } return success.report(__func__); } @@ -416,7 +523,13 @@ namespace GridKit "va_base": 100000000.0 }, "buses": [ - { "number": 1, "class": "bus", "name": "Bus 1", "init": { "Vr": 1.0, "Vi": 0.0 }, "v_base": 1.0 } + { + "number": 1, + "class": "bus", + "name": "Bus 1", + "init": { "Vr": 1.0, "Vi": 0.0 }, + "v_base": 1.0 + } ], "devices": [ { @@ -450,11 +563,16 @@ namespace GridKit success *= (data.regca.size() == 1); const auto& regca_data = data.regca[0]; success *= (regca_data.device_class == "Regca"); - success *= (regca_data.ports.at(PhasorDynamics::Converter::RegcaPorts::bus) == 1); - success *= (std::get_if(®ca_data.parameters.at(Params::P0)) != nullptr); - success *= (std::get_if(®ca_data.parameters.at(Params::Q0)) != nullptr); - success *= (std::get_if(®ca_data.parameters.at(Params::mva)) != nullptr); - success *= (std::get_if(®ca_data.parameters.at(Params::sL)) != nullptr); + success *= (regca_data.ports.at(PhasorDynamics::Converter::RegcaPorts::bus) + == 1); + success *= (std::get_if(®ca_data.parameters.at(Params::P0)) + != nullptr); + success *= (std::get_if(®ca_data.parameters.at(Params::Q0)) + != nullptr); + success *= (std::get_if(®ca_data.parameters.at(Params::mva)) + != nullptr); + success *= (std::get_if(®ca_data.parameters.at(Params::sL)) + != nullptr); PhasorDynamics::SystemModel system(data); success *= (system.allocate() == 0); @@ -478,14 +596,16 @@ namespace GridKit { TestStatus success = true; - auto dependency_tracking_jacobian = DependencyTrackingJacobian(); - auto enzyme_jacobian = EnzymeJacobian(); + auto dependency_tracking_jacobian = dependencyTrackingJacobian(); + auto enzyme_jacobian = enzymeJacobian(); success *= (dependency_tracking_jacobian.size() == enzyme_jacobian.size()); const auto nrows = std::min(dependency_tracking_jacobian.size(), enzyme_jacobian.size()); for (size_t i = 0; i < nrows; ++i) { - success *= isEqual(dependency_tracking_jacobian[i], enzyme_jacobian[i], static_cast(1.0e-8)); + success *= isEqual(dependency_tracking_jacobian[i], + enzyme_jacobian[i], + static_cast(1.0e-8)); } return success.report(__func__); @@ -502,7 +622,7 @@ namespace GridKit return static_cast(variable); } - auto makeTestData() -> PhasorDynamics::Converter::RegcaData + auto makeData() -> PhasorDynamics::Converter::RegcaData { using Ports = PhasorDynamics::Converter::RegcaPorts; using Mon = PhasorDynamics::Converter::RegcaMonitorableVariables; @@ -513,8 +633,12 @@ namespace GridKit data.ports[Ports::bus] = 1; data.monitored_variables.insert(Mon::ir); data.monitored_variables.insert(Mon::ii); + data.monitored_variables.insert(Mon::p); + data.monitored_variables.insert(Mon::q); - data.parameters[Params::mva] = static_cast(100); + data.parameters[Params::P0] = static_cast(0.0); + data.parameters[Params::Q0] = static_cast(0.0); + data.parameters[Params::mva] = static_cast(100.0); data.parameters[Params::Tg] = static_cast(0.02); data.parameters[Params::TM] = static_cast(0.02); data.parameters[Params::Rqmax] = static_cast(999.0); @@ -531,32 +655,76 @@ namespace GridKit return data; } - auto makeGoldenTestData(bool use_lvpl) -> PhasorDynamics::Converter::RegcaData + auto makeDynamicData() -> PhasorDynamics::Converter::RegcaData { - auto data = makeTestData(); + auto data = makeData(); + data.parameters[Params::P0] = static_cast(0.6); + data.parameters[Params::Q0] = static_cast(-0.2); data.parameters[Params::Tg] = static_cast(0.2); data.parameters[Params::TM] = static_cast(0.4); data.parameters[Params::Rqmax] = static_cast(0.5); data.parameters[Params::Rqmin] = static_cast(-0.6); data.parameters[Params::Rpmax] = static_cast(0.7); - data.parameters[Params::sL] = use_lvpl; data.parameters[Params::IL1] = static_cast(1.1); - data.parameters[Params::VL0] = static_cast(0.4); - data.parameters[Params::VL1] = static_cast(0.9); - data.parameters[Params::VA0] = static_cast(0.4); - data.parameters[Params::VA1] = static_cast(0.9); data.parameters[Params::Vhvmax] = static_cast(1.3); return data; } bool invalidParameterCase(PhasorDynamics::Bus& bus, Params param, RealT value) { - auto data = makeTestData(); + auto data = makeData(); data.parameters[param] = value; PhasorDynamics::Converter::Regca model(&bus, data); return model.verify() > 0; } + ScalarT toComponentBase(ScalarT value, RealT mva) const + { + return value * static_cast(100.0) / mva; + } + + ScalarT toSystemBase(ScalarT value, RealT mva) const + { + return value * mva / static_cast(100.0); + } + + ScalarT activeLower(ScalarT ip, RealT rpmax) const + { + const RealT mp = static_cast(100.0) * rpmax; + return -rpmax - (mp - rpmax) * Math::sigmoid(ip); + } + + ScalarT activeUpper(ScalarT ip, ScalarT il, RealT rpmax, bool use_lvpl) const + { + const RealT mp = static_cast(100.0) * rpmax; + const ScalarT sigma_ip = Math::sigmoid(ip); + ScalarT output = mp * (ONE - sigma_ip) + rpmax * sigma_ip; + if (use_lvpl) + { + output = mp * (ONE - sigma_ip) + rpmax * sigma_ip * Math::sigmoid(il - ip); + } + return output; + } + + bool allResidualsZero(PhasorDynamics::Converter::Regca& regca) const + { + bool success = true; + for (size_t i = 0; i < static_cast(regca.size()); ++i) + { + if (!isEqual(regca.getResidual()[i], static_cast(0.0), kTol)) + { + std::cout << "REGCA residual row " << i << " is " << regca.getResidual()[i] << "\n"; + success = false; + } + if (!isEqual(regca.yp()[i], static_cast(0.0), kTol)) + { + std::cout << "REGCA derivative row " << i << " is " << regca.yp()[i] << "\n"; + success = false; + } + } + return success; + } + bool vectorMatches(const std::vector& actual, const std::vector& expected, const char* label) const @@ -575,7 +743,18 @@ namespace GridKit return success; } - void setGoldenResidualState(PhasorDynamics::Converter::Regca& regca) + bool scalarMatches(ScalarT actual, ScalarT expected, const char* label) const + { + if (isEqual(actual, expected, kTol)) + { + return true; + } + + std::cout << label << " mismatch: " << actual << " != " << expected << "\n"; + return false; + } + + void setResidualState(PhasorDynamics::Converter::Regca& regca) { regca.y()[index(Vars::VM)] = static_cast(0.86); regca.y()[index(Vars::IQ)] = static_cast(-0.2); @@ -595,52 +774,60 @@ namespace GridKit regca.yp()[index(Vars::IP)] = static_cast(0.03); } - bool residualGoldenVectorCase(ScalarT initial_iqcmd, bool use_lvpl, const std::vector& expected) + std::vector expectedResidualForState(ScalarT ipcmd, ScalarT iqcmd) const { - bool success = true; - - PhasorDynamics::Bus bus(0.95, 0.25); - bus.allocate(); - bus.initialize(); - - auto data = makeGoldenTestData(use_lvpl); - PhasorDynamics::Converter::Regca regca(&bus, data); - - PhasorDynamics::SignalNode ipcmd_node; - PhasorDynamics::SignalNode iqcmd_node; - ScalarT ipcmd_value{0.9}; - ScalarT iqcmd_value{initial_iqcmd}; - IdxT ipcmd_index = 21; - IdxT iqcmd_index = 22; - ipcmd_node.set(&ipcmd_value, &ipcmd_index); - iqcmd_node.set(&iqcmd_value, &iqcmd_index); - regca.getSignals().template attachSignalNode(&ipcmd_node); - regca.getSignals().template attachSignalNode(&iqcmd_node); - - regca.allocate(); - regca.initialize(); - iqcmd_value = static_cast(0.1); - setGoldenResidualState(regca); - - bus.evaluateResidual(); - regca.evaluateResidual(); - - success *= vectorMatches(regca.getResidual(), expected, "REGCA residual"); - - return success; + constexpr RealT tg = static_cast(0.2); + constexpr RealT tm = static_cast(0.4); + constexpr RealT rqmin = static_cast(-0.6); + constexpr RealT rpmax = static_cast(0.7); + constexpr RealT vl0 = static_cast(0.4); + constexpr RealT vl1 = static_cast(0.9); + constexpr RealT va0 = static_cast(0.4); + constexpr RealT va1 = static_cast(0.9); + constexpr RealT il1 = static_cast(1.1); + constexpr ScalarT vr = static_cast(0.95); + constexpr ScalarT vi = static_cast(0.25); + constexpr ScalarT vm = static_cast(0.86); + constexpr ScalarT iq = static_cast(-0.2); + constexpr ScalarT ip = static_cast(0.85); + constexpr ScalarT vt = static_cast(0.98); + constexpr ScalarT ii = static_cast(0.18); + constexpr ScalarT iqext = static_cast(0.03); + constexpr ScalarT il = static_cast(0.72); + constexpr ScalarT ir = static_cast(0.5); + constexpr ScalarT lp = static_cast(-0.4); + constexpr ScalarT up = static_cast(0.3); + constexpr ScalarT pbr = static_cast(0.52); + constexpr ScalarT qbr = static_cast(-0.046); + constexpr ScalarT vmdot = static_cast(0.01); + constexpr ScalarT iqdot = static_cast(-0.02); + constexpr ScalarT ipdot = static_cast(0.03); + + const ScalarT iq_error = iqcmd - iq; + const ScalarT ip_error = ipcmd - ip; + const ScalarT lvacm = Math::linseg(vt, va0, va1, ONE); + const ScalarT qnet = iq - iqext; + + std::vector expected(static_cast(Vars::MAXIMUM), + static_cast(0.0)); + expected[index(Vars::VM)] = -tm * vmdot - vm + vt; + expected[index(Vars::IQ)] = -tg * iqdot + iq_error + Math::ramp(tg * rqmin - iq_error); + expected[index(Vars::IP)] = + -tg * ipdot + tg * lp + Math::ramp(ip_error - tg * lp) - Math::ramp(ip_error - tg * up); + expected[index(Vars::VT)] = -vt * vt + vr * vr + vi * vi; + expected[index(Vars::II)] = -vt * ii + lvacm * ip * vi - qnet * vr; + expected[index(Vars::IQEXTRA)] = -iqext + Math::ramp(vt - static_cast(1.3)); + expected[index(Vars::IL)] = -il + Math::linseg(vm, vl0, vl1, il1); + expected[index(Vars::IR)] = -vt * ir + lvacm * ip * vr + qnet * vi; + expected[index(Vars::LP)] = -lp + activeLower(ip, rpmax); + expected[index(Vars::UP)] = -up + activeUpper(ip, il, rpmax, true); + expected[index(Vars::PBR)] = -pbr + vr * ir + vi * ii; + expected[index(Vars::QBR)] = -qbr + vi * ir - vr * ii; + return expected; } #ifdef GRIDKIT_ENABLE_ENZYME - void setJacobianState(PhasorDynamics::Converter::Regca& regca, - PhasorDynamics::Bus& bus) - { - bus.y()[0] = static_cast(0.95); - bus.y()[1] = static_cast(0.25); - - setGoldenResidualState(regca); - } - - void setJacobianStateDep( + void setResidualStateDep( PhasorDynamics::Converter::Regca& regca, PhasorDynamics::Bus& bus) { @@ -665,11 +852,11 @@ namespace GridKit regca.yp()[index(Vars::IP)].setValue(0.03); } - std::vector DependencyTrackingJacobian() + std::vector dependencyTrackingJacobian() { using DepVar = DependencyTracking::Variable; - auto data = makeGoldenTestData(true); + auto data = makeDynamicData(); PhasorDynamics::Bus bus(DepVar{0.95}, DepVar{0.25}); PhasorDynamics::Converter::Regca regca(&bus, data); @@ -684,7 +871,7 @@ namespace GridKit bus.allocate(); regca.allocate(); bus.initialize(); - setJacobianStateDep(regca, bus); + setResidualStateDep(regca, bus); for (IdxT i = 0; i < regca.size(); ++i) { @@ -710,7 +897,8 @@ namespace GridKit static_cast(regca.size() + bus.size())); for (IdxT i = 0; i < regca.size(); ++i) { - dependencies[static_cast(i)] = regca.getResidual()[static_cast(i)].getDependencies(); + dependencies[static_cast(i)] = + regca.getResidual()[static_cast(i)].getDependencies(); } dependencies[static_cast(regca.size())] = bus.Ir().getDependencies(); dependencies[static_cast(regca.size() + 1)] = bus.Ii().getDependencies(); @@ -718,9 +906,9 @@ namespace GridKit return dependencies; } - std::vector EnzymeJacobian() + std::vector enzymeJacobian() { - auto data = makeGoldenTestData(true); + auto data = makeDynamicData(); PhasorDynamics::Bus bus(0.95, 0.25); PhasorDynamics::Converter::Regca regca(&bus, data); @@ -741,7 +929,7 @@ namespace GridKit } bus.initialize(); - setJacobianState(regca, bus); + setResidualState(regca); regca.updateTime(0.0, 1.0); ipcmd_node.set(&ipcmd_value, &ipcmd_index); diff --git a/tests/UnitTests/PhasorDynamics/runConverterRegcaTests.cpp b/tests/UnitTests/PhasorDynamics/runConverterRegcaTests.cpp index 0d3daf3b4..2e9c1c482 100644 --- a/tests/UnitTests/PhasorDynamics/runConverterRegcaTests.cpp +++ b/tests/UnitTests/PhasorDynamics/runConverterRegcaTests.cpp @@ -6,17 +6,18 @@ int main() GridKit::Testing::ConverterRegcaTests test; - result += test.constructor(); + result += test.constructionAndValidation(); result += test.parameterValidation(); - result += test.lifecycle(); - result += test.steadyStateInitializationGolden(); - result += test.attachedSignalInitialization(); - result += test.fallbackInit(); + result += test.initializesFromPowerFlowAndPublishesSignals(); + result += test.unconnectedCommandsRemainConstant(); + result += test.externalCommandsDriveRuntimeResidual(); result += test.invalidInitialization(); result += test.signalVerification(); result += test.nullBusVerification(); - result += test.residualGoldenVectors(); - result += test.busInjection(); + result += test.busInjectionUsesSystemBase(); + result += test.residualEquations(); + result += test.highVoltageReactiveCurrentRoot(); + result += test.limiterBranchCoverage(); result += test.jsonParseAndSystemAssembly(); #ifdef GRIDKIT_ENABLE_ENZYME result += test.jacobian();