From a40afcad1b2404e58ede60bdae4eaa6ba9bec6ab Mon Sep 17 00:00:00 2001 From: Nicholson Koukpaizan Date: Sun, 11 Oct 2026 10:20:23 -0400 Subject: [PATCH 1/5] Add hierarchical PhasorDynamics system model with a unit test. --- CHANGELOG.md | 1 + GridKit/Model/PhasorDynamics/INPUT_FORMAT.md | 37 +- GridKit/Model/PhasorDynamics/SystemModel.hpp | 38 +- .../Model/PhasorDynamics/SystemModelData.hpp | 10 + .../SystemModelDataJSONParser.hpp | 50 ++ .../SystemModelDependencyTracking.cpp | 8 + .../PhasorDynamics/SystemModelEnzyme.cpp | 216 +++++---- .../Model/PhasorDynamics/SystemModelImpl.hpp | 397 +++++++++++++--- tests/UnitTests/PhasorDynamics/CMakeLists.txt | 16 + .../PhasorDynamics/SubsystemFlat.json | 116 +++++ .../PhasorDynamics/SubsystemHierarchical.json | 142 ++++++ .../PhasorDynamics/SubsystemModelTests.hpp | 443 ++++++++++++++++++ .../PhasorDynamics/runSubsystemModelTests.cpp | 17 + 13 files changed, 1341 insertions(+), 150 deletions(-) create mode 100644 tests/UnitTests/PhasorDynamics/SubsystemFlat.json create mode 100644 tests/UnitTests/PhasorDynamics/SubsystemHierarchical.json create mode 100644 tests/UnitTests/PhasorDynamics/SubsystemModelTests.hpp create mode 100644 tests/UnitTests/PhasorDynamics/runSubsystemModelTests.cpp diff --git a/CHANGELOG.md b/CHANGELOG.md index 0ffbbeba06..8e041ecb65 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -10,6 +10,7 @@ - Added `enable` and `disable` solver events that put a branch in or out of service. - `verify()` returns `Model::ConfigurationChecks` listing every configuration problem instead of an error count. - Added `Model::ParameterReader` +- Added hierarchical PhasorDynamics system models. ## v0.2 diff --git a/GridKit/Model/PhasorDynamics/INPUT_FORMAT.md b/GridKit/Model/PhasorDynamics/INPUT_FORMAT.md index 8ac4231c8b..7e556daefc 100644 --- a/GridKit/Model/PhasorDynamics/INPUT_FORMAT.md +++ b/GridKit/Model/PhasorDynamics/INPUT_FORMAT.md @@ -35,9 +35,9 @@ also be encoded as [MessagePack](https://msgpack.org). ## Format The root object contains `header`, `params`, `buses`, and `devices`, and may -also contain `signals` and `monitors`. `header` contains information about -the case, `params` contains system-wide parameters, `buses` is an array of -electrical buses, and `devices` is an array of system devices. +also contain `signals`, `monitors` and `subsystems`. `header` contains +information about the case, `params` contains system-wide parameters, `buses` +is an array of electrical buses, and `devices` is an array of system devices. ### Header @@ -161,6 +161,37 @@ are specified: [Ieeest](Stabilizer/IEEEST/README.md) | the IEEEST stabilizer model [ConstantSignalSource](SignalSource/README.md) | Constant complex signal source +### Subsystems + +Contained in the optional `subsystems` key is an array of objects, each of +which describes a hierarchical subsystem with its own `buses` and `devices`, +and optional `signals` and `subsystems`. Subsystems use the root's `params`. +Each subsystem also requires: + +Key | Description +-- | -- +`id` | Subsystem ID (string), unique across the whole hierarchy +`level` | Hierarchy level (integer): the root system is level 0 (implied, not given in the file), its subsystems level 1, their subsystems level 2, and so on + +The `level` must match the depth of the subsystem. It is redundant in this +hierarchical layout, but makes the hierarchy explicit, so that a flat layout +can be supported later. +That flat layout will also need a `subsystem` tag (the ID of the containing +subsystem) on buses, devices and signals; this tag is not supported yet. + +The input must already be partitioned: + +- A system with subsystems contains only buses and subsystems: its + `devices` array must be empty. +- Devices may use buses of their own system and of its ancestors, but not + of sibling subsystems. +- Devices may only use signals of their own system. +- Bus numbers and signal IDs must be unique across the whole hierarchy. + +Not all of these are checked yet: duplicate subsystem IDs, duplicate signal +IDs across systems, and bus numbers repeated in sibling subsystems are not +detected. + ## Example File for a 2-Bus System ```json diff --git a/GridKit/Model/PhasorDynamics/SystemModel.hpp b/GridKit/Model/PhasorDynamics/SystemModel.hpp index 0748b9be02..79deff2acc 100644 --- a/GridKit/Model/PhasorDynamics/SystemModel.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModel.hpp @@ -3,6 +3,7 @@ #include #include #include +#include #include #include @@ -36,6 +37,14 @@ namespace GridKit * Model::Evaluator for the system model. This is still work in * progress and code is not optimized. * + * Subsystems own their buses and devices and may connect to the parent's + * buses (external buses). A system with subsystems contains only buses + * and subsystems. Only the top-level system allocates storage. + * + * evaluateJacobian() on the top-level system assembles the system Jacobian. + * On a subsystem, it assembles the contributions of the subsystem's buses + * and devices in global indices, of size offset() + size(). + * * @todo Address thread safety for the system model methods. * */ @@ -74,6 +83,9 @@ namespace GridKit SystemModel(const SystemModelData& data); + SystemModel(const SystemModelData& data, + const std::vector& external_buses); + virtual ~SystemModel(); int setGridKitComponentID(IdxT component_id) override; @@ -100,6 +112,7 @@ namespace GridKit void addBus(BusT* bus); void addComponent(ComponentT* component); + void addSubsystem(SystemModel* subsystem); void setSystemBase(RealT freq_system_base, RealT va_system_base); @@ -108,17 +121,34 @@ namespace GridKit ComponentT* getComponent(IdxT gridkit_component_id); Branch* getBranch(const std::string& id); + IdxT numSubsystems() const; + SystemModel* getSubsystem(IdxT index); + IdxT offset() const; + + const std::vector& getExternalBuses() const; + private: - std::vector buses_; - SignalNodeSetT signal_nodes_; - std::vector components_; + BusT* findBus(IdxT bus_id); + IdxT computeSize(); + void bindChildren(IdxT global_offset); + void addMonitorsTo(MonitorT& monitor) const; + void evaluateCooJacobians(); - std::map gridkit_bus_indices_; ///< Map between gridkit_bus_id and bus_id + std::vector buses_; + SignalNodeSetT signal_nodes_; + std::vector components_; + std::vector subsystems_; /// Map between branch id and gridkit_component_id /// @todo remove this map and use branch IDs to access individual branches. std::map gridkit_branch_indices_; + std::map gridkit_bus_indices_; ///< Map between gridkit_bus_id and bus_id + std::map external_buses_; ///< Parent's buses by bus_id, not owned + std::vector ordered_external_buses_; ///< External buses in construction order, not owned + + IdxT offset_{0}; ///< Global index of this system's first variable + bool is_subsystem_{false}; ///< Allocated and bound by a parent system bool owns_components_{false}; /// Variable monitor diff --git a/GridKit/Model/PhasorDynamics/SystemModelData.hpp b/GridKit/Model/PhasorDynamics/SystemModelData.hpp index a197f082c6..49a6e9dc01 100644 --- a/GridKit/Model/PhasorDynamics/SystemModelData.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModelData.hpp @@ -124,6 +124,16 @@ namespace GridKit /// Monitor sink specs std::vector monitor_sink; + + /// Subsystem ID, unique in the hierarchy; empty for the top system + std::string id; + + /// Hierarchy level: 0 for the top system, parent level + 1 for a subsystem + IdxT level{0}; + + /// Hierarchical subsystems; devices may use buses of this system and + /// its ancestors. System bases must match this system's. + std::vector subsystems; }; /** diff --git a/GridKit/Model/PhasorDynamics/SystemModelDataJSONParser.hpp b/GridKit/Model/PhasorDynamics/SystemModelDataJSONParser.hpp index 25d2287ff5..736853dfed 100644 --- a/GridKit/Model/PhasorDynamics/SystemModelDataJSONParser.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModelDataJSONParser.hpp @@ -2,6 +2,7 @@ #include #include +#include #include @@ -19,6 +20,9 @@ namespace GridKit using Log = ::GridKit::Utilities::Logger; using MonitorFormat = ::GridKit::Model::VariableMonitorFormat; + template + void parseSystemContents(const json& j, SystemModelData& sm); + /// JSON parser function implementation for the `SystemModelData` type /// /// See the `README.md` in `GridKit/Model/PhasorDynamics` for more information @@ -77,6 +81,13 @@ namespace GridKit } } + parseSystemContents(j, sm); + } + + /// Parse buses, signals, devices and subsystems of a system or subsystem + template + void parseSystemContents(const json& j, SystemModelData& sm) + { /// Gets all electrical buses j.at("buses").get_to(sm.bus); @@ -206,6 +217,45 @@ namespace GridKit throw std::runtime_error("JSON parser failed"); } } + + /// Gets hierarchical subsystems, which use this system's bases + if (j.contains("subsystems")) + { + if (!j.at("subsystems").empty() && !j.at("devices").empty()) + { + Log::error() << "\n\tA system with subsystems may only contain buses and subsystems; " + << "move its devices into a subsystem.\n"; + throw std::runtime_error("JSON parser failed"); + } + + for (auto& raw_subsystem : j.at("subsystems")) + { + SystemModelData subsystem; + subsystem.freq_base = sm.freq_base; + subsystem.va_base = sm.va_base; + subsystem.level = sm.level + 1; + + if (!raw_subsystem.contains("id") || !raw_subsystem.at("id").is_string()) + { + Log::error() << "\n\tEach subsystem requires a string \"id\".\n"; + throw std::runtime_error("JSON parser failed"); + } + raw_subsystem.at("id").get_to(subsystem.id); + + // The level is redundant here; it allows a flat layout later. + if (!raw_subsystem.contains("level") + || !raw_subsystem.at("level").is_number_integer() + || raw_subsystem.at("level").get() != static_cast(subsystem.level)) + { + Log::error() << "\n\tSubsystem \"" << subsystem.id << "\" requires \"level\": " + << subsystem.level << " (parent level + 1).\n"; + throw std::runtime_error("JSON parser failed"); + } + + parseSystemContents(raw_subsystem, subsystem); + sm.subsystems.push_back(std::move(subsystem)); + } + } } } // namespace PhasorDynamics } // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/SystemModelDependencyTracking.cpp b/GridKit/Model/PhasorDynamics/SystemModelDependencyTracking.cpp index f6105911eb..e4d3430550 100644 --- a/GridKit/Model/PhasorDynamics/SystemModelDependencyTracking.cpp +++ b/GridKit/Model/PhasorDynamics/SystemModelDependencyTracking.cpp @@ -25,6 +25,14 @@ namespace GridKit template int SystemModel::evaluateJacobian() { + // The top-level system builds the Jacobian from its whole residual. + // @todo Subsystem Jacobian from the subsystem's own residual. + if (is_subsystem_) + { + Log::error() << "Subsystem Jacobians are only available with Enzyme\n"; + return -1; + } + this->constructCsr(); // Log::misc() << "System DependencyTracking Jacobian\n"; diff --git a/GridKit/Model/PhasorDynamics/SystemModelEnzyme.cpp b/GridKit/Model/PhasorDynamics/SystemModelEnzyme.cpp index 5f847533e6..fe752d73cb 100644 --- a/GridKit/Model/PhasorDynamics/SystemModelEnzyme.cpp +++ b/GridKit/Model/PhasorDynamics/SystemModelEnzyme.cpp @@ -24,7 +24,12 @@ namespace GridKit has_jacobian = has_jacobian && bus->hasJacobian(); } - if (!has_jacobian) + for (const auto& subsystem : subsystems_) + { + has_jacobian = has_jacobian && subsystem->hasJacobian(); + } + + if (!has_jacobian && !is_subsystem_) { Log::warning() << "GridKit was built with Enzyme, but some models " "don't have Jacobians available. " @@ -35,16 +40,11 @@ namespace GridKit } /** - * @brief Evaluate system Jacobian using component-level Jacobians. - * - * - Evaluate component-level Jacobians (stored as CooMatrixT). - * - If not already constructed, construct system-level Jacobian. - * This used a system-level CooMatrix to deduplicate and sort the sparsity pattern. - * - If already constructed, reset Jacobian values to zero and accumulate contributions. - * + * @brief Evaluate the COO Jacobians of buses and devices, including + * subsystems' */ template - int SystemModel::evaluateJacobian() + void SystemModel::evaluateCooJacobians() { // Initialize bus Jacobians for (const auto& bus : buses_) @@ -58,36 +58,73 @@ namespace GridKit component->evaluateJacobian(); } - // Build or update system CSR Jacobian + for (const auto& subsystem : subsystems_) + { + subsystem->evaluateCooJacobians(); + } + } + + /** + * @brief Evaluate system Jacobian using component-level Jacobians. + * + * - Evaluate component-level Jacobians (stored as CooMatrixT), including subsystems'. + * - If not already constructed, construct system-level Jacobian. + * This used a system-level CooMatrix to deduplicate and sort the sparsity pattern. + * - If already constructed, reset Jacobian values to zero and accumulate contributions. + * + * A subsystem assembles its Jacobian in global indices. Its external buses + * precede its own variables, so the Jacobian size is offset() + size(). + * Entries on external buses come from its devices only, not from the + * external buses. + * + * @pre The top-level system is allocated. + * + * @return 0 on success + */ + template + int SystemModel::evaluateJacobian() + { + evaluateCooJacobians(); + + // Build or update system CSR Jacobian. All passes must visit systems in + // the same order. if (csr_jac_ == nullptr) { // Count the number of non-zeros IdxT nnz_dup = 0; - for (const auto& component : components_) - { - const auto component_jacobian = component->getCooJacobian(); - - if (component_jacobian != nullptr) - { - nnz_dup += component_jacobian->getNnz(); - } - else - { - Log::warning() << "A component has returned a nullptr Jacobian.\n"; - } - } - for (const auto& bus : buses_) + std::vector systems{this}; + for (size_t k = 0; k < systems.size(); ++k) { - auto bus_jacobian = bus->getCooJacobian(); + SystemModel* system = systems[k]; + systems.insert(systems.end(), system->subsystems_.begin(), system->subsystems_.end()); - if (bus_jacobian != nullptr) + for (const auto& component : system->components_) { - nnz_dup += bus_jacobian->getNnz(); + const auto component_jacobian = component->getCooJacobian(); + + if (component_jacobian != nullptr) + { + nnz_dup += component_jacobian->getNnz(); + } + else + { + Log::warning() << "A component has returned a nullptr Jacobian.\n"; + } } - else + + for (const auto& bus : system->buses_) { - Log::warning() << "A bus has returned a nullptr Jacobian.\n"; + auto bus_jacobian = bus->getCooJacobian(); + + if (bus_jacobian != nullptr) + { + nnz_dup += bus_jacobian->getNnz(); + } + else + { + Log::warning() << "A bus has returned a nullptr Jacobian.\n"; + } } } @@ -97,54 +134,63 @@ namespace GridKit RealT* vals_dup = new RealT[static_cast(nnz_dup)]; IdxT counter = 0; - for (const auto& component : components_) + + systems = {this}; + for (size_t k = 0; k < systems.size(); ++k) { - auto component_jacobian = component->getCooJacobian(); + SystemModel* system = systems[k]; + systems.insert(systems.end(), system->subsystems_.begin(), system->subsystems_.end()); - if (component_jacobian != nullptr) + for (const auto& component : system->components_) { - const IdxT* rows = component_jacobian->getRowData(); - const IdxT* columns = component_jacobian->getColData(); - const RealT* values = component_jacobian->getValues(); - for (IdxT i = 0; i < component_jacobian->getNnz(); ++i) + auto component_jacobian = component->getCooJacobian(); + + if (component_jacobian != nullptr) { - rows_dup[counter] = rows[i]; - cols_dup[counter] = columns[i]; - vals_dup[counter] = values[i]; - counter++; + const IdxT* rows = component_jacobian->getRowData(); + const IdxT* columns = component_jacobian->getColData(); + const RealT* values = component_jacobian->getValues(); + for (IdxT i = 0; i < component_jacobian->getNnz(); ++i) + { + rows_dup[counter] = rows[i]; + cols_dup[counter] = columns[i]; + vals_dup[counter] = values[i]; + counter++; + } + } + else + { + Log::warning() << "A component has returned a nullptr Jacobian.\n"; } } - else - { - Log::warning() << "A component has returned a nullptr Jacobian.\n"; - } - } - - for (const auto& bus : buses_) - { - auto bus_jacobian = bus->getCooJacobian(); - if (bus_jacobian != nullptr) + for (const auto& bus : system->buses_) { - const IdxT* rows = bus_jacobian->getRowData(); - const IdxT* columns = bus_jacobian->getColData(); - const RealT* values = bus_jacobian->getValues(); - for (IdxT i = 0; i < bus_jacobian->getNnz(); ++i) + auto bus_jacobian = bus->getCooJacobian(); + + if (bus_jacobian != nullptr) { - rows_dup[counter] = rows[i]; - cols_dup[counter] = columns[i]; - vals_dup[counter] = values[i]; - counter++; + const IdxT* rows = bus_jacobian->getRowData(); + const IdxT* columns = bus_jacobian->getColData(); + const RealT* values = bus_jacobian->getValues(); + for (IdxT i = 0; i < bus_jacobian->getNnz(); ++i) + { + rows_dup[counter] = rows[i]; + cols_dup[counter] = columns[i]; + vals_dup[counter] = values[i]; + counter++; + } + } + else + { + Log::warning() << "A bus has returned a nullptr Jacobian.\n"; } - } - else - { - Log::warning() << "A bus has returned a nullptr Jacobian.\n"; } } // Build the system COO Jacobian - CooMatrixT jac(size_, size_, nnz_dup, &rows_dup, &cols_dup, &vals_dup); + const IdxT n = offset_ + size_; + CooMatrixT jac(n, n, nnz_dup, &rows_dup, &cols_dup, &vals_dup); // Populate CSR data with sort and deduplicate IdxT* row_ptrs = jac.getCsrRowData(); @@ -160,7 +206,7 @@ namespace GridKit std::copy(jac.getValues(), jac.getValues() + nnz_, vals); // Create the CSR Jacobian - csr_jac_ = new CsrMatrixT(size_, size_, nnz_, &row_ptrs, &cols, &vals); + csr_jac_ = new CsrMatrixT(n, n, nnz_, &row_ptrs, &cols, &vals); const IdxT* map_to_sorted = jac.getMapToSorted(); const IdxT* map_to_dedup = jac.getMapToDeduplicated(); @@ -183,32 +229,40 @@ namespace GridKit // Update CSR values from component and bus Jacobians IdxT counter = 0; - for (const auto& component : components_) + + std::vector systems{this}; + for (size_t k = 0; k < systems.size(); ++k) { - auto component_jacobian = component->getCooJacobian(); + SystemModel* system = systems[k]; + systems.insert(systems.end(), system->subsystems_.begin(), system->subsystems_.end()); - if (component_jacobian != nullptr) + for (const auto& component : system->components_) { - const RealT* values = component_jacobian->getValues(); - for (IdxT i = 0; i < component_jacobian->getNnz(); ++i) + auto component_jacobian = component->getCooJacobian(); + + if (component_jacobian != nullptr) { - vals[map_to_csr_[counter]] += values[i]; - counter++; + const RealT* values = component_jacobian->getValues(); + for (IdxT i = 0; i < component_jacobian->getNnz(); ++i) + { + vals[map_to_csr_[counter]] += values[i]; + counter++; + } } } - } - - for (const auto& bus : buses_) - { - auto bus_jacobian = bus->getCooJacobian(); - if (bus_jacobian != nullptr) + for (const auto& bus : system->buses_) { - const RealT* values = bus_jacobian->getValues(); - for (IdxT i = 0; i < bus_jacobian->getNnz(); ++i) + auto bus_jacobian = bus->getCooJacobian(); + + if (bus_jacobian != nullptr) { - vals[map_to_csr_[counter]] += values[i]; - counter++; + const RealT* values = bus_jacobian->getValues(); + for (IdxT i = 0; i < bus_jacobian->getNnz(); ++i) + { + vals[map_to_csr_[counter]] += values[i]; + counter++; + } } } } diff --git a/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp b/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp index e5db0b4b43..aaed310e2d 100644 --- a/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp @@ -35,11 +35,25 @@ namespace GridKit * @pre SystemModelData contains consistent connectivity information * and physically meaningful model parameters. * - * @post All component models in SystemModelData are created, and - * correctly connected into the system model. + * @post All component models in SystemModelData, including its + * subsystems, are created, and correctly connected into the system model. */ template SystemModel::SystemModel(const SystemModelData& data) + : SystemModel(data, {}) + { + } + + /** + * @brief Construct a subsystem connected to its parent's buses (external buses) + * + * @param[in] data - Subsystem data with the buses it owns + * @param[in] external_buses - Parent's buses that subsystem devices may + * use + */ + template + SystemModel::SystemModel(const SystemModelData& data, + const std::vector& external_buses) : monitor_(std::make_unique(time_)) { using namespace Governor; @@ -66,6 +80,18 @@ namespace GridKit addBus(bus); } + for (auto* bus : external_buses) + { + if (gridkit_bus_indices_.contains(bus->busID()) + || !external_buses_.emplace(bus->busID(), bus).second) + { + Log::error() << "External bus " << bus->busID() + << " is already in the subsystem\n"; + throw std::runtime_error("SystemModel construction failed"); + } + ordered_external_buses_.push_back(bus); + } + // Add REGCA converters for (const auto& regcadata : data.regca) { @@ -297,6 +323,20 @@ namespace GridKit addComponent(source); } + // Add subsystems; they may use this system's buses and its external + // buses. Added after devices, so devices above cannot use subsystem buses. + if (!data.subsystems.empty()) + { + std::vector visible_buses(external_buses); + visible_buses.insert(visible_buses.end(), buses_.begin(), buses_.end()); + for (const auto& subsystem_data : data.subsystems) + { + auto subsystem = std::make_unique(subsystem_data, visible_buses); + addSubsystem(subsystem.get()); + subsystem.release(); + } + } + for (const auto& sink : data.monitor_sink) { monitor_->addSink(sink); @@ -314,6 +354,12 @@ namespace GridKit { if (owns_components_) { + // Subsystem devices may hold this system's buses; delete them first. + for (auto subsystem : subsystems_) + { + delete subsystem; + } + for (auto component : components_) { delete component; @@ -329,8 +375,7 @@ namespace GridKit /** * @brief Set component ID * - * @note Should default to 0. Nested system models are not currently - * supported. + * @note Should default to 0. */ template int SystemModel::setGridKitComponentID(IdxT component_id) @@ -347,8 +392,7 @@ namespace GridKit * * @pre Buses and components with nonzero size are unallocated or already * bound to external system storage. - * - * @note System model composition is flat; nested systems are not supported. + * @pre This system is not a subsystem. * * @throws std::runtime_error if storage allocation, child binding, model * verification, or initialization for sparse Jacobian discovery fails. @@ -356,17 +400,13 @@ namespace GridKit template int SystemModel::allocate() { - size_ = 0; - - for (const auto& bus : buses_) + if (is_subsystem_) { - size_ += bus->size(); + Log::error() << "A subsystem is allocated by its top-level system\n"; + throw std::runtime_error("SystemModel allocation failed"); } - for (const auto& component : components_) - { - size_ += component->size(); - } + computeSize(); // Allocate global vectors if (!allocated_) @@ -391,15 +431,96 @@ namespace GridKit throw std::runtime_error("SystemModel allocation failed"); } + bindChildren(0); + + // Verify component configuration + if (!this->verifyAndLog()) + { + throw std::runtime_error("SystemModel allocation failed"); + } + + // Sparse-pattern discovery requires an initialized operating point. A failed + // initialization aborts allocation before residual/Jacobian evaluation or + // monitor startup. + // @todo Replace with a sparsity analysis that sets the NNZ and allocates + // the Jacobian without needing the Jacobian values. + if (hasJacobian()) + { + const int status = initialize(); + if (status != 0) + { + Log::error() << "System model initialization failed with status " + << status << '\n'; + throw std::runtime_error("SystemModel allocation failed"); + } + evaluateResidual(); + evaluateJacobian(); + } + + initializeMonitor(); + startMonitor(); + + allocated_ = true; + return 0; + } + + /** + * @brief Set and return the size of this system, including subsystems + */ + template + typename SystemModel::IdxT + SystemModel::computeSize() + { + size_ = 0; + + for (const auto& bus : buses_) + { + size_ += bus->size(); + } + + for (const auto& component : components_) + { + size_ += component->size(); + } + + for (const auto& subsystem : subsystems_) + { + size_ += subsystem->computeSize(); + } + + return size_; + } + + /** + * @brief Bind buses, components and subsystems to this system's vectors + * + * @param[in] global_offset - Global index of this system's first variable + */ + template + void SystemModel::bindChildren(IdxT global_offset) + { + offset_ = global_offset; + + // New indices invalidate a subsystem's Jacobian. + if (is_subsystem_) + { + delete csr_jac_; + csr_jac_ = nullptr; + + delete[] map_to_csr_; + map_to_csr_ = nullptr; + + nnz_ = 0; + } + tag_.resize(size_); variable_indices_.resize(size_); residual_indices_.resize(size_); - // Default variable and residual index mapping to local index for (IdxT j = 0; j < size_; ++j) { - this->setVariableIndex(j, j); - this->setResidualIndex(j, j); + this->setVariableIndex(j, global_offset + j); + this->setResidualIndex(j, global_offset + j); } IdxT offset = 0; @@ -421,8 +542,8 @@ namespace GridKit for (IdxT j = 0; j < bus->size(); ++j) { - bus->setVariableIndex(j, offset + j); - bus->setResidualIndex(j, offset + j); + bus->setVariableIndex(j, global_offset + offset + j); + bus->setResidualIndex(j, global_offset + offset + j); } offset += bus->size(); @@ -445,48 +566,32 @@ namespace GridKit for (IdxT j = 0; j < component->size(); ++j) { - component->setVariableIndex(j, offset + j); - component->setResidualIndex(j, offset + j); + component->setVariableIndex(j, global_offset + offset + j); + component->setResidualIndex(j, global_offset + offset + j); } offset += component->size(); } - if (offset != size_) - { - Log::error() << "Bound vector sizes do not match the system size\n"; - throw std::runtime_error("SystemModel allocation failed"); - } - - // Verify component configuration - if (!this->verifyAndLog()) + for (const auto& subsystem : subsystems_) { - throw std::runtime_error("SystemModel allocation failed"); - } - - // Sparse-pattern discovery requires an initialized operating point. A failed - // initialization aborts allocation before residual/Jacobian evaluation or - // monitor startup. - // @todo Replace with a sparsity analysis that sets the NNZ and allocates - // the Jacobian without needing the Jacobian values. - if (hasJacobian()) - { - const int status = initialize(); - if (status != 0) + const int bind_status = subsystem->bind(y_, yp_, f_, abs_tol_, offset); + if (bind_status != 0) { - Log::error() << "System model initialization failed with status " - << status << '\n'; + Log::error() << "Failed to bind subsystem vectors to system storage\n"; throw std::runtime_error("SystemModel allocation failed"); } - evaluateResidual(); - evaluateJacobian(); - } - initializeMonitor(); - startMonitor(); + subsystem->bindChildren(global_offset + offset); - allocated_ = true; - return 0; + offset += subsystem->size(); + } + + if (offset != size_) + { + Log::error() << "Bound vector sizes do not match the system size\n"; + throw std::runtime_error("SystemModel allocation failed"); + } } /** @@ -510,6 +615,15 @@ namespace GridKit } } + for (const auto& subsystem : subsystems_) + { + const auto subsystem_checks = subsystem->verify(); + for (const auto& error : subsystem_checks.errors()) + { + checks.fail(std::format("subsystem {}: {}", subsystem->getGridKitComponentID(), error)); + } + } + return checks; } @@ -546,13 +660,22 @@ namespace GridKit status += component->initialize(); } + // Subsystem devices need the external buses initialized above. + for (const auto& subsystem : subsystems_) + { + status += subsystem->initialize(); + } + y_.setDataUpdated(); yp_.setDataUpdated(); - // For DependencyTracking::Variable, set variable numbers + // For DependencyTracking::Variable, set variable numbers at the top level if constexpr (std::is_same_v) { - this->initializeDependencyTrackingVariableNumbers(); + if (!is_subsystem_) + { + this->initializeDependencyTrackingVariableNumbers(); + } } return status; @@ -563,13 +686,22 @@ namespace GridKit */ template void SystemModel::initializeMonitor() + { + addMonitorsTo(*monitor_); + } + + /** + * @brief Add monitors from buses, components and subsystems to `monitor` + */ + template + void SystemModel::addMonitorsTo(MonitorT& monitor) const { for (const auto* bus : buses_) { auto* mon = bus->getMonitor(); if (mon && !mon->empty()) { - monitor_->addMonitor(mon); + monitor.addMonitor(mon); } } @@ -578,9 +710,14 @@ namespace GridKit auto* mon = component->getMonitor(); if (mon && !mon->empty()) { - monitor_->addMonitor(mon); + monitor.addMonitor(mon); } } + + for (const auto* subsystem : subsystems_) + { + subsystem->addMonitorsTo(monitor); + } } template @@ -617,13 +754,13 @@ namespace GridKit template int SystemModel::tagDifferentiable() { - // Set initial values for global solution vectors + // Set initial values for global solution vectors; tag_ starts at offset_ for (const auto& bus : buses_) { bus->tagDifferentiable(); for (IdxT j = 0; j < bus->size(); ++j) { - tag_[bus->getVariableIndex(j)] = bus->tag()[j]; + tag_[bus->getVariableIndex(j) - offset_] = bus->tag()[j]; } } @@ -632,7 +769,16 @@ namespace GridKit component->tagDifferentiable(); for (IdxT j = 0; j < component->size(); ++j) { - tag_[component->getVariableIndex(j)] = component->tag()[j]; + tag_[component->getVariableIndex(j) - offset_] = component->tag()[j]; + } + } + + for (const auto& subsystem : subsystems_) + { + subsystem->tagDifferentiable(); + for (IdxT j = 0; j < subsystem->size(); ++j) + { + tag_[subsystem->getVariableIndex(j) - offset_] = subsystem->tag()[j]; } } @@ -662,6 +808,11 @@ namespace GridKit component->setAbsoluteTolerance(rel_tol); } + for (const auto& subsystem : subsystems_) + { + subsystem->setAbsoluteTolerance(rel_tol); + } + abs_tol_.setDataUpdated(); return 0; @@ -679,6 +830,7 @@ namespace GridKit * from components. Buses assign their residual values, while components * add to those values by in-place addition. This is why (for now) bus * residuals need to be computed first. + * For the same reason, external buses are evaluated before subsystems. */ template int SystemModel::evaluateResidual() @@ -693,6 +845,11 @@ namespace GridKit component->evaluateResidual(); } + for (const auto& subsystem : subsystems_) + { + subsystem->evaluateResidual(); + } + f_.setDataUpdated(); return 0; @@ -711,6 +868,11 @@ namespace GridKit { component->updateTime(t, a); } + + for (const auto& subsystem : subsystems_) + { + subsystem->updateTime(t, a); + } } /** @@ -733,6 +895,8 @@ namespace GridKit * * Add component at the end of the components array and set GridKit's component ID * + * @pre This system has no subsystems. + * * @todo: No integer user-facing component_id for now, but we could map GridKit's * component ID to the disambiguation_string * @@ -740,6 +904,12 @@ namespace GridKit template void SystemModel::addComponent(ComponentT* component) { + if (!subsystems_.empty()) + { + Log::error() << "A system with subsystems cannot contain devices\n"; + throw std::runtime_error("SystemModel::addComponent failed"); + } + IdxT gridkit_component_id = static_cast(components_.size()); component->setGridKitComponentID(gridkit_component_id); component->setSystemBase(this->freq_system_base_, @@ -748,6 +918,39 @@ namespace GridKit allocated_ = false; } + /** + * @brief Add subsystem; deleted with this system if it owns its components + * + * @pre The subsystem is not allocated and not part of another system, and + * has the same system bases. This system has no devices. + */ + template + void SystemModel::addSubsystem(SystemModel* subsystem) + { + if (subsystem == nullptr || subsystem == this || subsystem->is_subsystem_) + { + Log::error() << "Invalid subsystem or subsystem already in a system\n"; + throw std::runtime_error("SystemModel::addSubsystem failed"); + } + + if (!components_.empty()) + { + Log::error() << "A system with devices cannot contain subsystems\n"; + throw std::runtime_error("SystemModel::addSubsystem failed"); + } + + if (subsystem->freq_system_base_ != this->freq_system_base_ + || subsystem->va_system_base_ != this->va_system_base_) + { + Log::error() << "Subsystem bases differ from the parent's\n"; + throw std::runtime_error("SystemModel::addSubsystem failed"); + } + + subsystem->is_subsystem_ = true; + subsystems_.push_back(subsystem); + allocated_ = false; + } + /** * @brief Set system bases and propagate them to existing components. * @@ -764,20 +967,90 @@ namespace GridKit { component->setSystemBase(freq_system_base, va_system_base); } + + for (auto* subsystem : subsystems_) + { + subsystem->setSystemBase(freq_system_base, va_system_base); + } } /** - * @brief Return pointer to a bus + * @brief Return pointer to a bus of this system, its subsystems, or an + * external bus * */ template SystemModel::BusT* SystemModel::getBus(IdxT bus_id) { - // Should fail if user-provided bus_id is incorrect - IdxT gridkit_bus_id = gridkit_bus_indices_.at(bus_id); - assert((buses_[gridkit_bus_id])->busID() == bus_id); - return buses_[gridkit_bus_id]; + BusT* bus = findBus(bus_id); + if (bus == nullptr) + { + // Should fail if user-provided bus_id is incorrect + return external_buses_.at(bus_id); + } + return bus; + } + + /** + * @brief Return pointer to a bus of this system or its subsystems, or + * nullptr if there is none + */ + template + SystemModel::BusT* + SystemModel::findBus(IdxT bus_id) + { + const auto it = gridkit_bus_indices_.find(bus_id); + if (it != gridkit_bus_indices_.end()) + { + IdxT gridkit_bus_id = it->second; + assert((buses_[gridkit_bus_id])->busID() == bus_id); + return buses_[gridkit_bus_id]; + } + + for (auto* subsystem : subsystems_) + { + BusT* bus = subsystem->findBus(bus_id); + if (bus != nullptr) + { + return bus; + } + } + return nullptr; + } + + template + typename SystemModel::IdxT + SystemModel::numSubsystems() const + { + return static_cast(subsystems_.size()); + } + + template + SystemModel* + SystemModel::getSubsystem(IdxT index) + { + return subsystems_.at(index); + } + + /** + * @brief Global index of this system's first variable and residual + */ + template + typename SystemModel::IdxT + SystemModel::offset() const + { + return offset_; + } + + /** + * @brief External buses of this subsystem, in construction order + */ + template + const std::vector::BusT*>& + SystemModel::getExternalBuses() const + { + return ordered_external_buses_; } /** diff --git a/tests/UnitTests/PhasorDynamics/CMakeLists.txt b/tests/UnitTests/PhasorDynamics/CMakeLists.txt index b430e75988..ebc662a498 100644 --- a/tests/UnitTests/PhasorDynamics/CMakeLists.txt +++ b/tests/UnitTests/PhasorDynamics/CMakeLists.txt @@ -194,6 +194,20 @@ target_link_libraries( configure_file(${CMAKE_CURRENT_SOURCE_DIR}/ThreeBusBasicBad.json ${CMAKE_CURRENT_BINARY_DIR}/ThreeBusBasicBad.json COPYONLY) +add_executable(test_phasor_subsystem_model runSubsystemModelTests.cpp) +target_link_libraries( + test_phasor_subsystem_model + GridKit::phasor_dynamics_systemmodel + GridKit::phasor_dynamics_systemmodel_dependency_tracking + GridKit::testing) +target_include_directories( + test_phasor_subsystem_model + PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/nlohmann-json/single_include) +configure_file(${CMAKE_CURRENT_SOURCE_DIR}/SubsystemFlat.json + ${CMAKE_CURRENT_BINARY_DIR}/SubsystemFlat.json COPYONLY) +configure_file(${CMAKE_CURRENT_SOURCE_DIR}/SubsystemHierarchical.json + ${CMAKE_CURRENT_BINARY_DIR}/SubsystemHierarchical.json COPYONLY) + add_executable(test_phasor_system_single_component runSystemSingleComponentTests.cpp) target_link_libraries( test_phasor_system_single_component @@ -232,6 +246,7 @@ add_test( NAME PhasorDynamicsSystemTest COMMAND test_phasor_system WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}) +add_test(NAME PhasorDynamicsSubsystemModelTest COMMAND test_phasor_subsystem_model) add_test(NAME PhasorDynamicsSystemSingleComponentTest COMMAND test_phasor_system_single_component) add_test(NAME PhasorDynamicsComponentConnectionTest COMMAND test_phasor_component_connection) @@ -256,6 +271,7 @@ install( test_phasor_stabilizer_ieeest test_phasor_gen_classical test_phasor_system + test_phasor_subsystem_model test_phasor_system_single_component test_phasor_component_connection RUNTIME DESTINATION bin) diff --git a/tests/UnitTests/PhasorDynamics/SubsystemFlat.json b/tests/UnitTests/PhasorDynamics/SubsystemFlat.json new file mode 100644 index 0000000000..9011724f0f --- /dev/null +++ b/tests/UnitTests/PhasorDynamics/SubsystemFlat.json @@ -0,0 +1,116 @@ +{ + "header": { + "case_name": "SubsystemFlat", + "case_description": "Flat system of the subsystem model tests", + "case_comments": "None" + }, + "params": { + "freq_base": 60.0, + "va_base": 100000000.0 + }, + "buses": [ + { + "number": 1, + "class": "Bus", + "name": "1", + "init": { + "Vr": 1.01, + "Vi": 0.1 + }, + "params": { + "kv": 138.0 + } + }, + { + "number": 2, + "class": "Bus", + "name": "2", + "init": { + "Vr": 1.02, + "Vi": 0.2 + }, + "params": { + "kv": 138.0 + } + }, + { + "number": 3, + "class": "Bus", + "name": "3", + "init": { + "Vr": 1.03, + "Vi": 0.3 + }, + "params": { + "kv": 138.0 + } + }, + { + "number": 4, + "class": "Bus", + "name": "4", + "init": { + "Vr": 1.04, + "Vi": 0.4 + }, + "params": { + "kv": 138.0 + } + } + ], + "devices": [ + { + "class": "Branch", + "ports": { + "bus1": 1, + "bus2": 2 + }, + "id": "1", + "params": { + "R": 0.03, + "X": 0.1, + "G": 0.0, + "B": 0.04 + } + }, + { + "class": "Branch", + "ports": { + "bus1": 1, + "bus2": 3 + }, + "id": "2", + "params": { + "R": 0.04, + "X": 0.1, + "G": 0.0, + "B": 0.06 + } + }, + { + "class": "Branch", + "ports": { + "bus1": 3, + "bus2": 4 + }, + "id": "3", + "params": { + "R": 0.07, + "X": 0.3, + "G": 0.0, + "B": 0.08 + } + }, + { + "class": "LoadZ", + "ports": { + "bus": 4 + }, + "id": "1", + "params": { + "R": 5.0, + "X": 0.8 + } + } + ] +} diff --git a/tests/UnitTests/PhasorDynamics/SubsystemHierarchical.json b/tests/UnitTests/PhasorDynamics/SubsystemHierarchical.json new file mode 100644 index 0000000000..4e5801dd3c --- /dev/null +++ b/tests/UnitTests/PhasorDynamics/SubsystemHierarchical.json @@ -0,0 +1,142 @@ +{ + "header": { + "case_name": "SubsystemHierarchical", + "case_description": "SubsystemFlat split into hierarchical subsystems", + "case_comments": "None" + }, + "params": { + "freq_base": 60.0, + "va_base": 100000000.0 + }, + "buses": [ + { + "number": 1, + "class": "Bus", + "name": "1", + "init": { + "Vr": 1.01, + "Vi": 0.1 + }, + "params": { + "kv": 138.0 + } + } + ], + "devices": [], + "subsystems": [ + { + "id": "A", + "level": 1, + "buses": [ + { + "number": 2, + "class": "Bus", + "name": "2", + "init": { + "Vr": 1.02, + "Vi": 0.2 + }, + "params": { + "kv": 138.0 + } + } + ], + "devices": [ + { + "class": "Branch", + "ports": { + "bus1": 1, + "bus2": 2 + }, + "id": "1", + "params": { + "R": 0.03, + "X": 0.1, + "G": 0.0, + "B": 0.04 + } + } + ] + }, + { + "id": "B", + "level": 1, + "buses": [ + { + "number": 3, + "class": "Bus", + "name": "3", + "init": { + "Vr": 1.03, + "Vi": 0.3 + }, + "params": { + "kv": 138.0 + } + } + ], + "devices": [], + "subsystems": [ + { + "id": "C", + "level": 2, + "buses": [ + { + "number": 4, + "class": "Bus", + "name": "4", + "init": { + "Vr": 1.04, + "Vi": 0.4 + }, + "params": { + "kv": 138.0 + } + } + ], + "devices": [ + { + "class": "Branch", + "ports": { + "bus1": 1, + "bus2": 3 + }, + "id": "2", + "params": { + "R": 0.04, + "X": 0.1, + "G": 0.0, + "B": 0.06 + } + }, + { + "class": "Branch", + "ports": { + "bus1": 3, + "bus2": 4 + }, + "id": "3", + "params": { + "R": 0.07, + "X": 0.3, + "G": 0.0, + "B": 0.08 + } + }, + { + "class": "LoadZ", + "ports": { + "bus": 4 + }, + "id": "1", + "params": { + "R": 5.0, + "X": 0.8 + } + } + ] + } + ] + } + ] +} diff --git a/tests/UnitTests/PhasorDynamics/SubsystemModelTests.hpp b/tests/UnitTests/PhasorDynamics/SubsystemModelTests.hpp new file mode 100644 index 0000000000..87894afb19 --- /dev/null +++ b/tests/UnitTests/PhasorDynamics/SubsystemModelTests.hpp @@ -0,0 +1,443 @@ +#pragma once + +#include +#include +#include +#include +#include + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace GridKit +{ + namespace Testing + { + using json = nlohmann::json; + + template + class SubsystemModelTests + { + private: + using ComponentT = PhasorDynamics::Component; + using RealT = typename ComponentT::RealT; + using DataT = PhasorDynamics::SystemModelData; + using SystemT = PhasorDynamics::SystemModel; + using BusT = typename SystemT::BusT; + + public: + SubsystemModelTests() = default; + ~SubsystemModelTests() = default; + + /** + * @brief Compare a flat system with the same system split into + * hierarchical subsystems built programmatically. + */ + TestOutcome flatVsHierarchical() + { + TestStatus success = true; + + // Data of each system of the hierarchical input without its + // subsystems: top, A, B, C + const json input = readJson("SubsystemHierarchical.json"); + const json& b_input = input["subsystems"][1]; + std::vector data; + for (json system : {input, input["subsystems"][0], b_input, b_input["subsystems"][0]}) + { + system.erase("subsystems"); + system["header"] = input["header"]; + system["params"] = input["params"]; + data.push_back(system.template get()); + } + + SystemT flat(readJson("SubsystemFlat.json").template get()); + + // The top system and B own (and delete) their subsystems. B contains + // only bus 3 and C, so C gets branch 2 to bus 1. + SystemT top(data[0]); + auto* a = new SystemT(data[1], {top.getBus(1)}); + auto* b = new SystemT(data[2], {top.getBus(1)}); + auto* c = new SystemT(data[3], {top.getBus(1), b->getBus(3)}); + b->addSubsystem(c); + top.addSubsystem(a); + top.addSubsystem(b); + + compareSystems(flat, top, success); + + return success.report(__func__); + } + + /** + * @brief Build the hierarchical system of flatVsHierarchical() from + * hierarchical input data. + */ + TestOutcome hierarchicalInput() + { + TestStatus success = true; + + const DataT top_data = readJson("SubsystemHierarchical.json").template get(); + + success *= top_data.subsystems.size() == 2; + success *= top_data.subsystems[1].subsystems.size() == 1; + success *= top_data.subsystems[1].subsystems[0].va_base == top_data.va_base; + success *= top_data.level == 0; + success *= top_data.subsystems[1].id == "B"; + success *= top_data.subsystems[1].level == 1; + success *= top_data.subsystems[1].subsystems[0].id == "C"; + success *= top_data.subsystems[1].subsystems[0].level == 2; + + SystemT flat(readJson("SubsystemFlat.json").template get()); + SystemT top(top_data); + + compareSystems(flat, top, success); + + return success.report(__func__); + } + + /** + * @brief Invalid hierarchies in the input data fail at construction. + */ + TestOutcome invalidHierarchy() + { + TestStatus success = true; + + const json hierarchical = readJson("SubsystemHierarchical.json"); + json input; + + // Device uses a bus of a sibling subsystem + input = hierarchical; + input["subsystems"][0]["devices"][0]["ports"]["bus2"] = 3; + success *= constructionThrows(input); + + // Device uses a signal of its parent + const json signal_r = {{"signal_id", 1}, {"name", "sr"}}; + const json signal_i = {{"signal_id", 2}, {"name", "si"}}; + const json source = {{"class", "ConstantSignalSource"}, + {"ports", {{"sr", 1}, {"si", 2}}}, + {"id", "1"}, + {"params", json::object()}}; + input = hierarchical; + input["signals"] = json::array({signal_r, signal_i}); + input["subsystems"][0]["devices"].push_back(source); + success *= constructionThrows(input); + + // Devices in a system with subsystems + input = hierarchical; + json load = input["subsystems"][1]["subsystems"][0]["devices"][2]; + load["ports"]["bus"] = 1; + input["devices"].push_back(load); + success *= constructionThrows(input); + + // Programmatic hierarchy: devices and subsystems in one system + { + SystemT with_device(readJson("SubsystemFlat.json").template get()); + SystemT with_subsystem; + SystemT subsystem; + SystemT other; + with_subsystem.addSubsystem(&subsystem); + + bool add_subsystem_throws = false; + try + { + with_device.addSubsystem(&other); + } + catch (const std::exception&) + { + add_subsystem_throws = true; + } + success *= add_subsystem_throws; + + bool add_component_throws = false; + try + { + with_subsystem.addComponent(with_device.getComponent(0)); + } + catch (const std::exception&) + { + add_component_throws = true; + } + success *= add_component_throws; + } + + // Subsystem without an ID + input = hierarchical; + input["subsystems"][0].erase("id"); + success *= constructionThrows(input); + + // Subsystem without a level + input = hierarchical; + input["subsystems"][0].erase("level"); + success *= constructionThrows(input); + + // Subsystem level not matching the hierarchy (grandchild at level 1) + input = hierarchical; + input["subsystems"][1]["subsystems"][0]["level"] = 1; + success *= constructionThrows(input); + + return success.report(__func__); + } + + /** + * @brief Subsystems must have their parent's system bases, which + * setSystemBase() passes down the hierarchy. + */ + TestOutcome systemBases() + { + TestStatus success = true; + + // Programmatically built systems do not own their subsystems. + SystemT top; + SystemT a; + SystemT b; + top.addSubsystem(&a); + + b.setSystemBase(50.0, 100.0e6); + try + { + a.addSubsystem(&b); + success *= false; + } + catch (const std::exception&) + { + } + + // A gets the new bases from top, so B now fits. + top.setSystemBase(50.0, 100.0e6); + try + { + a.addSubsystem(&b); + } + catch (const std::exception&) + { + success *= false; + } + + return success.report(__func__); + } + +#ifdef GRIDKIT_ENABLE_ENZYME + /** + * @brief Compare the hierarchical Jacobian with the flat one and with + * DependencyTracking, and the subsystem Jacobians with the top-level + * one, before and after a fault changes the Jacobian values. + */ + TestOutcome jacobian() + { + using DependencyMap = DependencyTracking::Variable::DependencyMap; + + TestStatus success = true; + + SystemT flat(readJson("SubsystemFlat.json").template get()); + SystemT top(readJson("SubsystemHierarchical.json").template get()); + success *= top.numSubsystems() == 2 && top.getSubsystem(1)->numSubsystems() == 1; + if (top.numSubsystems() != 2 || top.getSubsystem(1)->numSubsystems() != 1) + { + return success.report(__func__); + } + auto* a = top.getSubsystem(0); + auto* b = top.getSubsystem(1); + auto* c = b->getSubsystem(0); + + PhasorDynamics::SystemModel top_dt( + readJson("SubsystemHierarchical.json").template get()); + + for (auto* system : {&flat, &top}) + { + success *= system->allocate() == 0; + success *= system->initialize() == 0; + } + success *= top_dt.allocate() == 0; + success *= top_dt.initialize() == 0; + + const IdxT n = top.size(); + + // Duplicate entries may be summed in a different order. + const RealT tol = 10 * std::numeric_limits::epsilon(); + + for (const bool fault : {false, true}) + { + if (fault) + { + success *= flat.getBus(4)->setFault(true, 0.0, 0.01) == 0; + success *= top.getBus(4)->setFault(true, 0.0, 0.01) == 0; + success *= top_dt.getBus(4)->setFault(true, 0.0, 0.01) == 0; + } + for (auto* system : {&flat, &top}) + { + success *= system->evaluateResidual() == 0; + success *= system->evaluateJacobian() == 0; + } + success *= top_dt.evaluateResidual() == 0; + success *= top_dt.evaluateJacobian() == 0; + if (flat.getCsrJacobian() == nullptr || top.getCsrJacobian() == nullptr + || top_dt.getCsrJacobian() == nullptr) + { + success *= false; + return success.report(__func__); + } + + // Both systems have the same ordering, see compareSystems(). + const auto jacobian_flat = MapFromCsr(flat.getCsrJacobian()); + const auto jacobian_top = MapFromCsr(top.getCsrJacobian()); + const auto jacobian_dt = MapFromCsr(top_dt.getCsrJacobian()); + for (IdxT i = 0; i < n; ++i) + { + success *= isEqual(jacobian_top[i], jacobian_flat[i]); + success *= isEqual(jacobian_top[i], jacobian_dt[i], tol); + } + + for (auto* subsystem : {a, b, c}) + { + success *= subsystem->evaluateJacobian() == 0; + if (subsystem->getCsrJacobian() == nullptr) + { + success *= false; + return success.report(__func__); + } + + // Subsystem Jacobian in global indices. Coupling with external + // variables is not checked yet. + const auto jacobian = MapFromCsr(subsystem->getCsrJacobian()); + const IdxT begin = subsystem->offset(); + const IdxT end = begin + subsystem->size(); + success *= static_cast(jacobian.size()) == end; + if (static_cast(jacobian.size()) != end) + { + return success.report(__func__); + } + + for (IdxT i = begin; i < end; ++i) + { + DependencyMap expected_row; + DependencyMap actual_row; + for (const auto& [column, value] : jacobian_top[i]) + { + if (begin <= column && column < end) + { + expected_row[column] = value; + } + } + for (const auto& [column, value] : jacobian[i]) + { + if (begin <= column && column < end) + { + actual_row[column] = value; + } + } + success *= isEqual(actual_row, expected_row, tol); + } + } + } + + return success.report(__func__); + } +#endif + + private: + /// Read a JSON file copied next to the test executable + static json readJson(const std::string& file_name) + { + std::ifstream file(file_name); + return json::parse(file); + } + + /// Whether constructing a system from the input throws + static bool constructionThrows(const json& input) + { + try + { + const DataT data = input.template get(); + SystemT system(data); + } + catch (const std::exception&) + { + return true; + } + return false; + } + + /** + * @brief Compare top = [bus 1 | A[bus 2] | B[bus 3 | C[bus 4, load]]] + * with the flat system, before and after a fault on bus 4. + */ + static void compareSystems(SystemT& flat, SystemT& top, TestStatus& success) + { + for (auto* system : {&flat, &top}) + { + success *= system->allocate() == 0; + success *= system->initialize() == 0; + } + + success *= top.numSubsystems() == 2; + if (top.numSubsystems() != 2) + { + return; + } + auto* a = top.getSubsystem(0); + auto* b = top.getSubsystem(1); + success *= b->numSubsystems() == 1; + if (b->numSubsystems() != 1) + { + return; + } + auto* c = b->getSubsystem(0); + + // Branches have no variables, so both systems order variables as + // buses 1-4, then the load. + success *= a->offset() == top.getBus(1)->size(); + success *= a->size() == a->getBus(2)->size(); + success *= b->offset() == a->offset() + a->size(); + success *= b->size() == b->getBus(3)->size() + c->size(); + success *= c->offset() == b->offset() + b->getBus(3)->size(); + success *= c->size() == c->getBus(4)->size() + c->getComponent(2)->size(); + success *= top.size() == top.getBus(1)->size() + a->size() + b->size(); + success *= top.size() == flat.size(); + if (top.size() != flat.size()) + { + return; + } + + // External buses + BusT* bus1 = top.getBus(1); + BusT* bus3 = b->getBus(3); + success *= a->getExternalBuses() == std::vector{bus1}; + success *= b->getExternalBuses() == std::vector{bus1}; + success *= c->getExternalBuses() == std::vector{bus1, bus3}; + + // The top system finds subsystem buses. + success *= top.getBus(3) == bus3; + success *= top.getBus(4) == c->getBus(4); + + for (const bool fault : {false, true}) + { + if (fault) + { + success *= flat.getBus(4)->setFault(true, 0.0, 0.01) == 0; + success *= top.getBus(4)->setFault(true, 0.0, 0.01) == 0; + } + + for (auto* system : {&flat, &top}) + { + success *= system->evaluateResidual() == 0; + } + const auto* f_flat = flat.getResidual().getData(); + const auto* f_top = top.getResidual().getData(); + for (IdxT i = 0; i < flat.size(); ++i) + { + success *= isEqual(f_top[i], f_flat[i]); + } + } + } + }; + } // namespace Testing +} // namespace GridKit diff --git a/tests/UnitTests/PhasorDynamics/runSubsystemModelTests.cpp b/tests/UnitTests/PhasorDynamics/runSubsystemModelTests.cpp new file mode 100644 index 0000000000..400149408b --- /dev/null +++ b/tests/UnitTests/PhasorDynamics/runSubsystemModelTests.cpp @@ -0,0 +1,17 @@ +#include "SubsystemModelTests.hpp" + +int main() +{ + GridKit::Testing::TestingResults result; + GridKit::Testing::SubsystemModelTests test; + + result += test.flatVsHierarchical(); + result += test.hierarchicalInput(); + result += test.invalidHierarchy(); + result += test.systemBases(); +#ifdef GRIDKIT_ENABLE_ENZYME + result += test.jacobian(); +#endif + + return result.summary(); +} From 87561b88fcd95f12934a44abe3bf40fc7ff1b2fe Mon Sep 17 00:00:00 2001 From: Nicholson Koukpaizan Date: Wed, 14 Oct 2026 10:50:10 -0400 Subject: [PATCH 2/5] Revert CHANGELOG update for now. --- CHANGELOG.md | 1 - 1 file changed, 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 8e041ecb65..0ffbbeba06 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -10,7 +10,6 @@ - Added `enable` and `disable` solver events that put a branch in or out of service. - `verify()` returns `Model::ConfigurationChecks` listing every configuration problem instead of an error count. - Added `Model::ParameterReader` -- Added hierarchical PhasorDynamics system models. ## v0.2 From 262152b81a5b410eed226aeb8c27a7241a2a2e9d Mon Sep 17 00:00:00 2001 From: Nicholson Koukpaizan Date: Wed, 14 Oct 2026 11:50:12 -0400 Subject: [PATCH 3/5] Call SystemModel default constructor, so the destructor gets called if exceptions are thrown. --- GridKit/Model/PhasorDynamics/SystemModelImpl.hpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp b/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp index aaed310e2d..98976c9e8b 100644 --- a/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp @@ -54,7 +54,7 @@ namespace GridKit template SystemModel::SystemModel(const SystemModelData& data, const std::vector& external_buses) - : monitor_(std::make_unique(time_)) + : SystemModel() { using namespace Governor; using namespace Exciter; @@ -63,6 +63,7 @@ namespace GridKit using namespace Converter; owns_components_ = true; + monitor_ = std::make_unique(time_); // Store parsed system bases before constructing data-driven components. this->setSystemBase(data.freq_base, data.va_base); From 5ffaba8a7a27e56a4f5008feb782ea1dbd1103de Mon Sep 17 00:00:00 2001 From: Nicholson Koukpaizan Date: Wed, 14 Oct 2026 12:26:30 -0400 Subject: [PATCH 4/5] Add components before connecting ports so they are freed if construction throws. --- .../Model/PhasorDynamics/SystemModelImpl.hpp | 48 +++++++++---------- 1 file changed, 24 insertions(+), 24 deletions(-) diff --git a/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp b/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp index 98976c9e8b..8732201596 100644 --- a/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp @@ -103,9 +103,9 @@ namespace GridKit } auto* regca = new Regca(getBus(bus_index), regcadata); - regca->getPorts().connect(regcadata, signal_nodes_); - addComponent(regca); + + regca->getPorts().connect(regcadata, signal_nodes_); } // Add branches @@ -167,8 +167,8 @@ namespace GridKit bus_index = gendata.buses.at(GenrouBuses::bus); } auto* gen = new Genrou(getBus(bus_index), gendata); - gen->getPorts().connect(gendata, signal_nodes_); addComponent(gen); + gen->getPorts().connect(gendata, signal_nodes_); } // Add GENSAL generators @@ -181,8 +181,8 @@ namespace GridKit } auto* gen = new Gensal(getBus(bus_index), gendata); - gen->getPorts().connect(gendata, signal_nodes_); addComponent(gen); + gen->getPorts().connect(gendata, signal_nodes_); } // Add classical generators @@ -194,8 +194,8 @@ namespace GridKit bus_index = gendata.buses.at(GenClassicalBuses::bus); } auto* gen = new GenClassical(getBus(bus_index), gendata); - gen->getPorts().connect(gendata, signal_nodes_); addComponent(gen); + gen->getPorts().connect(gendata, signal_nodes_); } // Add REECB after its current-command and feedback producers because @@ -209,36 +209,36 @@ namespace GridKit } auto* reecb = new Reecb(bus, reecbdata); - reecb->getPorts().connect(reecbdata, signal_nodes_); - addComponent(reecb); + + reecb->getPorts().connect(reecbdata, signal_nodes_); } // Add Tgov1 governors for (const auto& govdata : data.gov) { auto* gov = new Tgov1(govdata); - gov->getPorts().connect(govdata, signal_nodes_); - addComponent(gov); + + gov->getPorts().connect(govdata, signal_nodes_); } // Add GASTPTI governors for (const auto& gastptidata : data.gastpti) { auto* gastpti = new GastPti(gastptidata); - gastpti->getPorts().connect(gastptidata, signal_nodes_); - addComponent(gastpti); + + gastpti->getPorts().connect(gastptidata, signal_nodes_); } // Add HYGOV governors for (const auto& hygovdata : data.hygov) { auto* hygov = new Hygov(hygovdata); - hygov->getPorts().connect(hygovdata, signal_nodes_); - addComponent(hygov); + + hygov->getPorts().connect(hygovdata, signal_nodes_); } // Add IEEEST stabilizers before exciters that consume their output during @@ -246,8 +246,8 @@ namespace GridKit for (const auto& stabdata : data.stabilizer) { auto* stabilizer = new Ieeest(stabdata); - stabilizer->getPorts().connect(stabdata, signal_nodes_); addComponent(stabilizer); + stabilizer->getPorts().connect(stabdata, signal_nodes_); } for (const auto& excitedata : data.exciter) @@ -259,9 +259,9 @@ namespace GridKit } auto* exciter = new Ieeet1(getBus(bus_index), excitedata); - exciter->getPorts().connect(excitedata, signal_nodes_); - addComponent(exciter); + + exciter->getPorts().connect(excitedata, signal_nodes_); } for (const auto& excitedata : data.esdc1a) @@ -273,9 +273,9 @@ namespace GridKit } auto* exciter = new Esdc1a(bus, excitedata); - exciter->getPorts().connect(excitedata, signal_nodes_); - addComponent(exciter); + + exciter->getPorts().connect(excitedata, signal_nodes_); } for (const auto& excitedata : data.sexspti) @@ -287,9 +287,9 @@ namespace GridKit } auto* exciter = new SexsPti(getBus(bus_index), excitedata); - exciter->getPorts().connect(excitedata, signal_nodes_); - addComponent(exciter); + + exciter->getPorts().connect(excitedata, signal_nodes_); } // Add REPCA plant controllers after the signal producers they read at @@ -303,25 +303,25 @@ namespace GridKit } auto* repca = new Repca(bus, repcadata); - repca->getPorts().connect(repcadata, signal_nodes_); - addComponent(repca); + + repca->getPorts().connect(repcadata, signal_nodes_); } // Add constant signal sources for (const auto& srcdata : data.constant_source) { auto* source = new ConstantSignalSource(srcdata); - source->getPorts().connect(srcdata, signal_nodes_); addComponent(source); + source->getPorts().connect(srcdata, signal_nodes_); } // Add constant signal sources for (const auto& srcdata : data.function_source) { auto* source = new FunctionSignalSource(srcdata); - source->getPorts().connect(srcdata, signal_nodes_); addComponent(source); + source->getPorts().connect(srcdata, signal_nodes_); } // Add subsystems; they may use this system's buses and its external From f967e963782d988c9cf1642702f281a96693eed7 Mon Sep 17 00:00:00 2001 From: Nicholson Koukpaizan Date: Wed, 14 Oct 2026 18:30:39 -0400 Subject: [PATCH 5/5] Set subsystem IDs in addSubsystem so verify() errors name the right subsystem. --- GridKit/Model/PhasorDynamics/SystemModelImpl.hpp | 1 + 1 file changed, 1 insertion(+) diff --git a/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp b/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp index 8732201596..3dffd672e1 100644 --- a/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp @@ -948,6 +948,7 @@ namespace GridKit } subsystem->is_subsystem_ = true; + subsystem->setGridKitComponentID(static_cast(subsystems_.size())); subsystems_.push_back(subsystem); allocated_ = false; }