From b6fe2d5f6ad89b9e3e9bbe6977bbf62c7e2e415f Mon Sep 17 00:00:00 2001 From: lukelowry Date: Mon, 5 Oct 2026 16:13:13 -0500 Subject: [PATCH] cleanup/simplify initial PR --- CHANGELOG.md | 1 + GridKit/Model/PhasorDynamics/Bus/Bus.hpp | 5 + .../Model/PhasorDynamics/Bus/BusEnzyme.cpp | 11 +- GridKit/Model/PhasorDynamics/Bus/BusImpl.hpp | 32 ++++- GridKit/Model/PhasorDynamics/Bus/README.md | 6 +- GridKit/Model/PhasorDynamics/BusBase.hpp | 7 ++ .../TenGen/Genrou/TenGenGenrou.cpp | 7 +- tests/UnitTests/PhasorDynamics/BusTests.hpp | 114 ++++++++++++++++++ tests/UnitTests/PhasorDynamics/CMakeLists.txt | 6 +- .../UnitTests/PhasorDynamics/runBusTests.cpp | 4 + 10 files changed, 176 insertions(+), 17 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index a63cd6b078..65d190ac74 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,6 +4,7 @@ - Added `BusSignalVoltageOut` bus model with voltage signal outlets and current signal inlets. - Added `BusSignalVoltageIn` bus model with voltage signal inlets and current signal outlets. +- Added `Bus::setFault` to apply or clear a fault to ground directly at a bus. ## v0.2 diff --git a/GridKit/Model/PhasorDynamics/Bus/Bus.hpp b/GridKit/Model/PhasorDynamics/Bus/Bus.hpp index cc2d9f495d..b6548d2024 100644 --- a/GridKit/Model/PhasorDynamics/Bus/Bus.hpp +++ b/GridKit/Model/PhasorDynamics/Bus/Bus.hpp @@ -47,6 +47,7 @@ namespace GridKit virtual ~Bus(); virtual int setBusID(IdxT) override final; + virtual int setFault(bool status, RealT R, RealT X) override final; virtual int allocate() override final; virtual int tagDifferentiable() override final; virtual int setAbsoluteTolerance(RealT rel_tol) override final; @@ -160,6 +161,10 @@ namespace GridKit private: ScalarT Vr0_{0.0}; ScalarT Vi0_{0.0}; + + /* Fault admittance */ + RealT fault_g_{0.0}; + RealT fault_b_{0.0}; }; } // namespace PhasorDynamics diff --git a/GridKit/Model/PhasorDynamics/Bus/BusEnzyme.cpp b/GridKit/Model/PhasorDynamics/Bus/BusEnzyme.cpp index 25230d60c4..8de69cf3a5 100644 --- a/GridKit/Model/PhasorDynamics/Bus/BusEnzyme.cpp +++ b/GridKit/Model/PhasorDynamics/Bus/BusEnzyme.cpp @@ -13,7 +13,7 @@ namespace GridKit /** * @brief Jacobian evaluation experimental. * - * This sets values to 0, and these remain unchanged. It is needed to get + * This sets values to the fault admittance, zero when cleared. It is needed to get * the indices into the list of entries that will later be deduplicated. * Contributions to bus Jacobians from other components are stored in those components. * @@ -36,14 +36,15 @@ namespace GridKit J_cols_buffer_[1] = variable_indices_.at(1); J_cols_buffer_[2] = variable_indices_.at(0); J_cols_buffer_[3] = variable_indices_.at(1); - J_vals_buffer_[0] = 0.0; - J_vals_buffer_[1] = 0.0; - J_vals_buffer_[2] = 0.0; - J_vals_buffer_[3] = 0.0; nnz_ = 4; this->constructCoo(); } + + J_vals_buffer_[0] = -fault_g_; + J_vals_buffer_[1] = fault_b_; + J_vals_buffer_[2] = -fault_b_; + J_vals_buffer_[3] = -fault_g_; return 0; } diff --git a/GridKit/Model/PhasorDynamics/Bus/BusImpl.hpp b/GridKit/Model/PhasorDynamics/Bus/BusImpl.hpp index 3b431c0abc..74cc89ba7d 100644 --- a/GridKit/Model/PhasorDynamics/Bus/BusImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Bus/BusImpl.hpp @@ -136,6 +136,32 @@ namespace GridKit return 0; } + /** + * @brief Apply or clear a fault to ground + * + * @param[in] status - true applies the fault, false clears it + * @param[in] R - fault resistance [p.u.] + * @param[in] X - fault reactance [p.u.] + */ + template + int Bus::setFault(bool status, RealT R, RealT X) + { + if (status && !(ZERO < R * R + X * X)) + { + Log::error() << "Bus: fault impedance R + jX must be nonzero\n"; + return 1; + } + + fault_g_ = 0.0; + fault_b_ = 0.0; + if (status) + { + fault_g_ = R / (X * X + R * R); + fault_b_ = -X / (X * X + R * R); + } + return 0; + } + /*! * @brief Bus variables are algebraic. */ @@ -193,7 +219,7 @@ namespace GridKit } /*! - * @brief PQ bus does not compute residuals, so here we just reset residual values. + * @brief Reset the current balance to the fault current. * * @warning This implementation assumes bus residuals are always evaluated * _before_ component model residuals. @@ -204,8 +230,8 @@ namespace GridKit { auto* f = f_.getData(); - f[0] = 0.0; - f[1] = 0.0; + f[0] = -(fault_g_ * Vr() - fault_b_ * Vi()); + f[1] = -(fault_b_ * Vr() + fault_g_ * Vi()); f_.setDataUpdated(); return 0; } diff --git a/GridKit/Model/PhasorDynamics/Bus/README.md b/GridKit/Model/PhasorDynamics/Bus/README.md index f012aa5c4c..d89165683f 100644 --- a/GridKit/Model/PhasorDynamics/Bus/README.md +++ b/GridKit/Model/PhasorDynamics/Bus/README.md @@ -69,12 +69,12 @@ None. #### Algebraic -Let $\mathcal{D}$ denote the set of components connected to the bus. +Let $\mathcal{D}$ denote the set of components connected to the bus, and let $G_f$, $B_f$ be the conductance and susceptance of a fault applied to the Bus, respectively, with $G_f + jB_f = 1/(R_f + jX_f)$, both zero when no fault is applied. ```math \begin{aligned} -0 &= \sum_{d \in \mathcal{D}} I_{r,d} \\ -0 &= \sum_{d \in \mathcal{D}} I_{i,d} +0 &= -(G_f V_r - B_f V_i) + \sum_{d \in \mathcal{D}} I_{r,d} \\ +0 &= -(B_f V_r + G_f V_i) + \sum_{d \in \mathcal{D}} I_{i,d} \end{aligned} ``` diff --git a/GridKit/Model/PhasorDynamics/BusBase.hpp b/GridKit/Model/PhasorDynamics/BusBase.hpp index 7fed07b6cb..989d4280d5 100644 --- a/GridKit/Model/PhasorDynamics/BusBase.hpp +++ b/GridKit/Model/PhasorDynamics/BusBase.hpp @@ -239,6 +239,13 @@ namespace GridKit virtual int setBusID(IdxT) = 0; + /// Apply (status true) or clear a fault to ground with impedance R + jX. + virtual int setFault(bool /* status */, RealT /* R */, RealT /* X */) + { + Log::error() << "Faults are not supported by this bus type\n"; + return 1; + } + virtual const IdxT busID() const { return bus_id_; diff --git a/tests/IntegrationTests/PhasorDynamics/TenGen/Genrou/TenGenGenrou.cpp b/tests/IntegrationTests/PhasorDynamics/TenGen/Genrou/TenGenGenrou.cpp index 4e0d6797b2..5bb6627b1a 100644 --- a/tests/IntegrationTests/PhasorDynamics/TenGen/Genrou/TenGenGenrou.cpp +++ b/tests/IntegrationTests/PhasorDynamics/TenGen/Genrou/TenGenGenrou.cpp @@ -71,8 +71,6 @@ int main() Genrou gen9(&bus9, 0.5, -0.09662372, 3., 0., 0., 7., .04, .05, .75, 2.1, 0.2, 0.18, 0.5, 0.5, 0.18, 0.15, 0., 0.); Genrou gen10(&bus10, 0.5, -0.09932297, 3., 0., 0., 7., .04, .05, .75, 2.1, 0.2, 0.18, 0.5, 0.5, 0.18, 0.15, 0., 0.); - BusFault fault(&bus10, 0, 1e-5, 0); - /* Connect everything together */ SystemModel sys; @@ -104,7 +102,6 @@ int main() sys.addComponent(&gen8); sys.addComponent(&gen9); sys.addComponent(&gen10); - sys.addComponent(&fault); sys.allocate(); real_type dt = 1.0 / 4.0 / 60.0; @@ -171,7 +168,7 @@ int main() } // Introduce fault to ground and run for 0.1s - fault.setStatus(1); + bus10.setFault(true, 0.0, 1e-5); ida.initializeSimulation(1.0); ida.runSimulation(1.1, dt, output_cb); @@ -183,7 +180,7 @@ int main() success *= isEqual(gen10.y().getData()[omega_index], omega_ref, 5e-5); // Clear fault and run until t = 10s. - fault.setStatus(0); + bus10.setFault(false, 0.0, 1e-5); ida.initializeSimulation(1.1); ida.runSimulation(10.0, dt, output_cb); real_type stop = static_cast(clock()); diff --git a/tests/UnitTests/PhasorDynamics/BusTests.hpp b/tests/UnitTests/PhasorDynamics/BusTests.hpp index a1ebb75936..7ae466059f 100644 --- a/tests/UnitTests/PhasorDynamics/BusTests.hpp +++ b/tests/UnitTests/PhasorDynamics/BusTests.hpp @@ -1,6 +1,10 @@ +#include #include #include +#include +#include +#include #include #include #include @@ -13,6 +17,9 @@ namespace GridKit template class BusTests { + private: + using RealT = typename PhasorDynamics::Bus::RealT; + public: BusTests() = default; ~BusTests() = default; @@ -96,6 +103,113 @@ namespace GridKit return success.report(__func__); } + + /// Fault current applied and cleared through setFault + TestOutcome fault() + { + TestStatus success = true; + + ScalarT Vr{1.0}; + ScalarT Vi{2.0}; + RealT R{1.0}; + RealT X{2.0}; + + PhasorDynamics::Bus bus(Vr, Vi); + bus.allocate(); + bus.initialize(); + + // Fault current I = -V / (R + jX) + const std::complex current = -std::complex(Vr, Vi) / std::complex(R, X); + + success *= bus.setFault(true, R, X) == 0; + bus.evaluateResidual(); + success *= isEqual(bus.Ir(), current.real()); + success *= isEqual(bus.Ii(), current.imag()); + + success *= bus.setFault(false, R, X) == 0; + bus.evaluateResidual(); + success *= isEqual(bus.Ir(), 0.0); + success *= isEqual(bus.Ii(), 0.0); + + // A zero fault impedance is rejected + success *= bus.setFault(true, 0.0, 0.0) != 0; + + // An infinite bus cannot be faulted + PhasorDynamics::BusInfinite bus_inf; + success *= bus_inf.setFault(true, R, X) != 0; + + return success.report(__func__); + } + +#ifdef GRIDKIT_ENABLE_ENZYME + /// Fault Jacobian matches DependencyTracking and keeps its pattern when cleared + TestOutcome jacobian() + { + TestStatus success = true; + + RealT R{1.0}; + RealT X{2.0}; + + // Jacobian via DependencyTracking, which numbers y entries 2 * index + DependencyTracking::Variable Vr{1.0}; + DependencyTracking::Variable Vi{2.0}; + + PhasorDynamics::Bus dependency_tracking_bus(Vr, Vi); + dependency_tracking_bus.allocate(); + dependency_tracking_bus.initialize(); + dependency_tracking_bus.setFault(true, R, X); + dependency_tracking_bus.evaluateResidual(); + + const auto* f = dependency_tracking_bus.getResidual().getData(); + const auto size = static_cast(dependency_tracking_bus.size()); + + std::vector dependency_tracking_jacobian(size); + for (size_t row = 0; row < size; ++row) + { + for (const auto& [number, value] : f[row].getDependencies()) + { + dependency_tracking_jacobian[row][number / 2] = value; + } + } + + // Jacobian from the bus + PhasorDynamics::Bus bus(1.0, 2.0); + bus.allocate(); + bus.initialize(); + bus.setFault(true, R, X); + bus.evaluateResidual(); + bus.evaluateJacobian(); + + auto* jacobian = bus.getCooJacobian(); + const IdxT* rows = jacobian->getRowData(); + const IdxT* cols = jacobian->getColData(); + const RealT* vals = jacobian->getValues(); + + std::vector bus_jacobian(size); + for (IdxT i = 0; i < jacobian->getNnz(); ++i) + { + bus_jacobian[static_cast(rows[i])][static_cast(cols[i])] = vals[i]; + } + + for (size_t row = 0; row < size; ++row) + { + success *= isEqual(dependency_tracking_jacobian[row], bus_jacobian[row]); + } + + // Clearing the fault keeps the same entries with zero values + const IdxT nnz = jacobian->getNnz(); + bus.setFault(false, R, X); + bus.evaluateResidual(); + bus.evaluateJacobian(); + success *= bus.nnz() == nnz; + for (IdxT i = 0; i < nnz; ++i) + { + success *= isEqual(vals[i], 0.0); + } + + return success.report(__func__); + } +#endif }; } // namespace Testing diff --git a/tests/UnitTests/PhasorDynamics/CMakeLists.txt b/tests/UnitTests/PhasorDynamics/CMakeLists.txt index 173be7e552..08f8bc0cd1 100644 --- a/tests/UnitTests/PhasorDynamics/CMakeLists.txt +++ b/tests/UnitTests/PhasorDynamics/CMakeLists.txt @@ -5,7 +5,11 @@ add_executable(test_phasor_bus runBusTests.cpp) target_link_libraries( - test_phasor_bus GridKit::phasor_dynamics_bus GridKit::testing) + test_phasor_bus + GridKit::definitions + GridKit::phasor_dynamics_bus + GridKit::phasor_dynamics_bus_dependency_tracking + GridKit::testing) add_executable(test_phasor_bus_fault runBusFaultTests.cpp) target_link_libraries( diff --git a/tests/UnitTests/PhasorDynamics/runBusTests.cpp b/tests/UnitTests/PhasorDynamics/runBusTests.cpp index c7063318d6..8c086f275c 100644 --- a/tests/UnitTests/PhasorDynamics/runBusTests.cpp +++ b/tests/UnitTests/PhasorDynamics/runBusTests.cpp @@ -10,6 +10,10 @@ int main() result += test.constructor(); result += test.residual(); + result += test.fault(); +#ifdef GRIDKIT_ENABLE_ENZYME + result += test.jacobian(); +#endif return result.summary(); }