From c74201ab5778bdf1c7ee1f165cc5f2b567aa20ae Mon Sep 17 00:00:00 2001 From: Alexander Novotny Date: Mon, 5 May 2025 16:09:13 -0400 Subject: [PATCH 1/5] Keep sub-evaluators up to date after integration step Evaluators which contain "children" evaluators such as systems should be given the chance to keep their sub-evaluators up to date after events like an integrator step. With the `update()` method, an external source can notify the evaluator when its state has been updated so that it can then go and update its children. This is useful for callbacks, since callbacks will probably want to query individual components rather than systems. --- examples/PhasorDynamics/Example1/example1.cpp | 8 +- src/Model/Evaluator.hpp | 6 + src/Model/PhasorDynamics/SystemModel.hpp | 229 +++++++----------- src/Solver/Dynamic/Ida.cpp | 2 + 4 files changed, 96 insertions(+), 149 deletions(-) diff --git a/examples/PhasorDynamics/Example1/example1.cpp b/examples/PhasorDynamics/Example1/example1.cpp index 1c8508d1..f70c61c8 100644 --- a/examples/PhasorDynamics/Example1/example1.cpp +++ b/examples/PhasorDynamics/Example1/example1.cpp @@ -91,12 +91,10 @@ int main() // push it into output, which is updated outside the callback. auto output_cb = [&](double t) { - std::vector& yval = sys.y(); - output.push_back({.ti = t, - .Vr = yval[0], - .Vi = yval[1], - .dw = yval[3]}); + .Vr = bus1.Vr(), + .Vi = bus1.Vi(), + .dw = gen.y()[1]}); }; // The above lambda is equivalent to writing diff --git a/src/Model/Evaluator.hpp b/src/Model/Evaluator.hpp index b4a447dd..3f94153e 100644 --- a/src/Model/Evaluator.hpp +++ b/src/Model/Evaluator.hpp @@ -56,6 +56,12 @@ namespace GridKit virtual void setTolerances(real_type& rtol, real_type& atol) const = 0; virtual void setMaxSteps(IdxT& msa) const = 0; + /// @brief Let the `Evaluator` know that its state has been updated externally (for example by an integrator), + /// in case it needs to keep other things up to date based on its internal state. + virtual void update() + { + } + virtual std::vector& y() = 0; virtual const std::vector& y() const = 0; diff --git a/src/Model/PhasorDynamics/SystemModel.hpp b/src/Model/PhasorDynamics/SystemModel.hpp index 03347b22..24864357 100644 --- a/src/Model/PhasorDynamics/SystemModel.hpp +++ b/src/Model/PhasorDynamics/SystemModel.hpp @@ -343,42 +343,15 @@ namespace GridKit */ int evaluateIntegrand() { - // Update variables - IdxT varOffset = 0; - IdxT optOffset = 0; + updateChildren(); + for (const auto& bus : buses_) { - for (IdxT j = 0; j < bus->size(); ++j) - { - bus->y()[j] = y_[varOffset + j]; - bus->yp()[j] = yp_[varOffset + j]; - } - varOffset += bus->size(); - - for (IdxT j = 0; j < bus->sizeParams(); ++j) - { - bus->param()[j] = param_[optOffset + j]; - } - optOffset += bus->sizeParams(); - bus->evaluateIntegrand(); } for (const auto& component : components_) { - for (IdxT j = 0; j < component->size(); ++j) - { - component->y()[j] = y_[varOffset + j]; - component->yp()[j] = yp_[varOffset + j]; - } - varOffset += component->size(); - - for (IdxT j = 0; j < component->sizeParams(); ++j) - { - component->param()[j] = param_[optOffset + j]; - } - optOffset += component->sizeParams(); - component->evaluateIntegrand(); } @@ -413,45 +386,9 @@ namespace GridKit */ int initializeAdjoint() { - IdxT offset = 0; - IdxT optOffset = 0; - - // Update bus variables and optimization parameters - for (const auto& bus : buses_) - { - for (IdxT j = 0; j < bus->size(); ++j) - { - bus->y()[j] = y_[offset + j]; - bus->yp()[j] = yp_[offset + j]; - } - offset += bus->size(); + updateChildren(); - for (IdxT j = 0; j < bus->sizeParams(); ++j) - { - bus->param()[j] = param_[optOffset + j]; - } - optOffset += bus->sizeParams(); - } - - // Update component variables and optimization parameters - for (const auto& component : components_) - { - for (IdxT j = 0; j < component->size(); ++j) - { - component->y()[j] = y_[offset + j]; - component->yp()[j] = yp_[offset + j]; - } - offset += component->size(); - - for (IdxT j = 0; j < component->sizeParams(); ++j) - { - component->param()[j] = param_[optOffset + j]; - } - optOffset += component->sizeParams(); - } - - // Reset counter - offset = 0; + size_t offset = 0; // Initialize bus adjoints for (const auto& bus : buses_) @@ -491,45 +428,7 @@ namespace GridKit */ int evaluateAdjointResidual() { - IdxT varOffset = 0; - IdxT optOffset = 0; - - // Update variables in component models - for (const auto& bus : buses_) - { - for (IdxT j = 0; j < bus->size(); ++j) - { - bus->y()[j] = y_[varOffset + j]; - bus->yp()[j] = yp_[varOffset + j]; - bus->yB()[j] = yB_[varOffset + j]; - bus->ypB()[j] = ypB_[varOffset + j]; - } - varOffset += bus->size(); - - for (IdxT j = 0; j < bus->sizeParams(); ++j) - { - bus->param()[j] = param_[optOffset + j]; - } - optOffset += bus->sizeParams(); - } - - for (const auto& component : components_) - { - for (IdxT j = 0; j < component->size(); ++j) - { - component->y()[j] = y_[varOffset + j]; - component->yp()[j] = yp_[varOffset + j]; - component->yB()[j] = yB_[varOffset + j]; - component->ypB()[j] = ypB_[varOffset + j]; - } - varOffset += component->size(); - - for (IdxT j = 0; j < component->sizeParams(); ++j) - { - component->param()[j] = param_[optOffset + j]; - } - optOffset += component->sizeParams(); - } + updateChildrenAdjoint(); for (const auto& bus : buses_) { @@ -575,44 +474,7 @@ namespace GridKit */ int evaluateAdjointIntegrand() { - // First, update variables - IdxT varOffset = 0; - IdxT optOffset = 0; - for (const auto& bus : buses_) - { - for (IdxT j = 0; j < bus->size(); ++j) - { - bus->y()[j] = y_[varOffset + j]; - bus->yp()[j] = yp_[varOffset + j]; - bus->yB()[j] = yB_[varOffset + j]; - bus->ypB()[j] = ypB_[varOffset + j]; - } - varOffset += bus->size(); - - for (IdxT j = 0; j < bus->sizeParams(); ++j) - { - bus->param()[j] = param_[optOffset + j]; - } - optOffset += bus->sizeParams(); - } - - for (const auto& component : components_) - { - for (IdxT j = 0; j < component->size(); ++j) - { - component->y()[j] = y_[varOffset + j]; - component->yp()[j] = yp_[varOffset + j]; - component->yB()[j] = yB_[varOffset + j]; - component->ypB()[j] = ypB_[varOffset + j]; - } - varOffset += component->size(); - - for (IdxT j = 0; j < component->sizeParams(); ++j) - { - component->param()[j] = param_[optOffset + j]; - } - optOffset += component->sizeParams(); - } + updateChildrenAdjoint(); // Evaluate integrand and update global vector for (const auto& component : components_) @@ -648,10 +510,89 @@ namespace GridKit components_.push_back(component); } + void update() override + { + updateChildrenAdjoint(); + } + private: std::vector buses_; std::vector components_; + /// @brief Whenever the internal state of a system is updated, the system + /// needs to reflect those changes in its "children" aka buses and + /// components. `updateChildren` will only update `y`, `yp`, and + /// `param` - call when you know only the internal state of those + /// variables has been updated. + void updateChildren() + { + // Update variables + IdxT varOffset = 0; + IdxT optOffset = 0; + for (const auto& bus : buses_) + { + for (IdxT j = 0; j < bus->size(); ++j) + { + bus->y()[j] = y_[varOffset + j]; + bus->yp()[j] = yp_[varOffset + j]; + } + varOffset += bus->size(); + + for (IdxT j = 0; j < bus->sizeParams(); ++j) + { + bus->param()[j] = param_[optOffset + j]; + } + optOffset += bus->sizeParams(); + } + + for (const auto& component : components_) + { + for (IdxT j = 0; j < component->size(); ++j) + { + component->y()[j] = y_[varOffset + j]; + component->yp()[j] = yp_[varOffset + j]; + } + varOffset += component->size(); + + for (IdxT j = 0; j < component->sizeParams(); ++j) + { + component->param()[j] = param_[optOffset + j]; + } + optOffset += component->sizeParams(); + } + } + + /// @brief Whenever the internal state of a system is updated, the system + /// needs to reflect those changes in its "children" aka buses and + /// components. `updateChildrenAdjoint` will update all internal + /// variables of this system's children, including adjoint variables. + void updateChildrenAdjoint() + { + // Update all non-adjoint state + updateChildren(); + + IdxT varOffset = 0; + IdxT optOffset = 0; + for (const auto& bus : buses_) + { + for (IdxT j = 0; j < bus->size(); ++j) + { + bus->yB()[j] = yB_[varOffset + j]; + bus->ypB()[j] = ypB_[varOffset + j]; + } + varOffset += bus->size(); + } + + for (const auto& component : components_) + { + for (IdxT j = 0; j < component->size(); ++j) + { + component->yB()[j] = yB_[varOffset + j]; + component->ypB()[j] = ypB_[varOffset + j]; + } + varOffset += component->size(); + } + } }; // class SystemModel } // namespace PhasorDynamics diff --git a/src/Solver/Dynamic/Ida.cpp b/src/Solver/Dynamic/Ida.cpp index 7cce626d..3d4ba474 100644 --- a/src/Solver/Dynamic/Ida.cpp +++ b/src/Solver/Dynamic/Ida.cpp @@ -231,6 +231,7 @@ namespace AnalysisManager model_->updateTime(tret, 0.0); copyVec(yy_, model_->y()); copyVec(yp_, model_->yp()); + model_->update(); (*step_callback)(tret); } @@ -246,6 +247,7 @@ namespace AnalysisManager model_->updateTime(tf, 0.0); copyVec(yy_, model_->y()); copyVec(yp_, model_->yp()); + model_->update(); // std::cout << "\n"; return retval; From 4b6edba1b989da95e664258ec92b6412f9d97317 Mon Sep 17 00:00:00 2001 From: Alexander Novotny Date: Tue, 6 May 2025 12:09:03 -0400 Subject: [PATCH 2/5] Make sure to notify children in a system when updated Since other systems could be children in a system, it's important to notify them when they've been updated. --- src/Model/PhasorDynamics/SystemModel.hpp | 44 ++++++++++++++++++++---- 1 file changed, 38 insertions(+), 6 deletions(-) diff --git a/src/Model/PhasorDynamics/SystemModel.hpp b/src/Model/PhasorDynamics/SystemModel.hpp index 24864357..d30fcd45 100644 --- a/src/Model/PhasorDynamics/SystemModel.hpp +++ b/src/Model/PhasorDynamics/SystemModel.hpp @@ -343,7 +343,9 @@ namespace GridKit */ int evaluateIntegrand() { - updateChildren(); + // No need to notify children yet, since we're about to + // evaluate them + updateChildren(); for (const auto& bus : buses_) { @@ -386,6 +388,8 @@ namespace GridKit */ int initializeAdjoint() { + // No need to notify children yet, since we're about to + // initialize them updateChildren(); size_t offset = 0; @@ -428,7 +432,9 @@ namespace GridKit */ int evaluateAdjointResidual() { - updateChildrenAdjoint(); + // No need to notify children yet, since we're about to + // evaluate them + updateChildrenAdjoint(); for (const auto& bus : buses_) { @@ -474,7 +480,9 @@ namespace GridKit */ int evaluateAdjointIntegrand() { - updateChildrenAdjoint(); + // No need to notify children yet, since we're about to + // evaluate them + updateChildrenAdjoint(); // Evaluate integrand and update global vector for (const auto& component : components_) @@ -512,7 +520,7 @@ namespace GridKit void update() override { - updateChildrenAdjoint(); + updateChildrenAdjoint(); } private: @@ -524,6 +532,11 @@ namespace GridKit /// components. `updateChildren` will only update `y`, `yp`, and /// `param` - call when you know only the internal state of those /// variables has been updated. + /// @tparam NOTIFY_CHILDREN Whether or not buses and components should + /// be notified through their `update()` method that their state + /// has been updated. If you are planning to notify them otherwise, + /// set this to `false`. + template void updateChildren() { // Update variables @@ -543,6 +556,9 @@ namespace GridKit bus->param()[j] = param_[optOffset + j]; } optOffset += bus->sizeParams(); + + if constexpr (NOTIFY_CHILDREN) + bus->update(); } for (const auto& component : components_) @@ -559,6 +575,9 @@ namespace GridKit component->param()[j] = param_[optOffset + j]; } optOffset += component->sizeParams(); + + if constexpr (NOTIFY_CHILDREN) + component->update(); } } @@ -566,10 +585,17 @@ namespace GridKit /// needs to reflect those changes in its "children" aka buses and /// components. `updateChildrenAdjoint` will update all internal /// variables of this system's children, including adjoint variables. + /// @tparam NOTIFY_CHILDREN Whether or not buses and components should + /// be notified through their `update()` method that their state + /// has been updated. If you are planning to notify them otherwise, + /// set this to `false`. + template void updateChildrenAdjoint() { - // Update all non-adjoint state - updateChildren(); + // Update all non-adjoint state. + // No need to notify children that they've been updated yet, + // because we still have more updating to do. + updateChildren(); IdxT varOffset = 0; IdxT optOffset = 0; @@ -581,6 +607,9 @@ namespace GridKit bus->ypB()[j] = ypB_[varOffset + j]; } varOffset += bus->size(); + + if constexpr (NOTIFY_CHILDREN) + bus->update(); } for (const auto& component : components_) @@ -591,6 +620,9 @@ namespace GridKit component->ypB()[j] = ypB_[varOffset + j]; } varOffset += component->size(); + + if constexpr (NOTIFY_CHILDREN) + component->update(); } } }; // class SystemModel From 7ef877e8c87e98677f47b03f2d9565d3cb96abf6 Mon Sep 17 00:00:00 2001 From: Alexander Novotny Date: Tue, 6 May 2025 12:19:56 -0400 Subject: [PATCH 3/5] evaluateResidual() uses updateChildren() --- src/Model/PhasorDynamics/SystemModel.hpp | 33 +++--------------------- 1 file changed, 4 insertions(+), 29 deletions(-) diff --git a/src/Model/PhasorDynamics/SystemModel.hpp b/src/Model/PhasorDynamics/SystemModel.hpp index d30fcd45..e2d458b7 100644 --- a/src/Model/PhasorDynamics/SystemModel.hpp +++ b/src/Model/PhasorDynamics/SystemModel.hpp @@ -263,42 +263,17 @@ namespace GridKit */ int evaluateResidual() { - // Update variables - IdxT varOffset = 0; - IdxT optOffset = 0; + // No need to notify children yet, since we're about to + // evaluate them + updateChildren(); + for (const auto& bus : buses_) { - for (IdxT j = 0; j < bus->size(); ++j) - { - bus->y()[j] = y_[varOffset + j]; - bus->yp()[j] = yp_[varOffset + j]; - } - varOffset += bus->size(); - - for (IdxT j = 0; j < bus->sizeParams(); ++j) - { - bus->param()[j] = param_[optOffset + j]; - } - optOffset += bus->sizeParams(); - bus->evaluateResidual(); } for (const auto& component : components_) { - for (IdxT j = 0; j < component->size(); ++j) - { - component->y()[j] = y_[varOffset + j]; - component->yp()[j] = yp_[varOffset + j]; - } - varOffset += component->size(); - - for (IdxT j = 0; j < component->sizeParams(); ++j) - { - component->param()[j] = param_[optOffset + j]; - } - optOffset += component->sizeParams(); - component->evaluateResidual(); } From 0ac539107a118c92b69ea5f264aa709ed7f4a91a Mon Sep 17 00:00:00 2001 From: Alexander Novotny Date: Tue, 6 May 2025 14:25:43 -0400 Subject: [PATCH 4/5] Update example2 --- examples/PhasorDynamics/Example2/example2.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/examples/PhasorDynamics/Example2/example2.cpp b/examples/PhasorDynamics/Example2/example2.cpp index fd580dce..11634b5e 100644 --- a/examples/PhasorDynamics/Example2/example2.cpp +++ b/examples/PhasorDynamics/Example2/example2.cpp @@ -114,10 +114,10 @@ int main() std::vector& yval = sys.y(); output.push_back({.t = t, - .gen2speed = 1 + yval[5], - .gen3speed = 1 + yval[26], - .v2mag = sqrt(yval[0] * yval[0] + yval[1] * yval[1]), - .v3mag = sqrt(yval[2] * yval[2] + yval[3] * yval[3])}); + .gen2speed = 1 + gen2.y()[1], + .gen3speed = 1 + gen3.y()[1], + .v2mag = sqrt(bus2.Vr() * bus2.Vr() + bus2.Vi() * bus2.Vi()), + .v3mag = sqrt(bus3.Vr() * bus3.Vr() + bus3.Vi() * bus3.Vi())}); }; /* Set up simulation */ From 38bf189372fd3d9906cf93f9db3eeb1b13bf7d28 Mon Sep 17 00:00:00 2001 From: Alexander Novotny Date: Tue, 6 May 2025 14:26:59 -0400 Subject: [PATCH 5/5] Remove unused variable --- examples/PhasorDynamics/Example2/example2.cpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/examples/PhasorDynamics/Example2/example2.cpp b/examples/PhasorDynamics/Example2/example2.cpp index 11634b5e..fdab584f 100644 --- a/examples/PhasorDynamics/Example2/example2.cpp +++ b/examples/PhasorDynamics/Example2/example2.cpp @@ -111,8 +111,6 @@ int main() auto output_cb = [&](real_type t) { - std::vector& yval = sys.y(); - output.push_back({.t = t, .gen2speed = 1 + gen2.y()[1], .gen3speed = 1 + gen3.y()[1],