From e75f402324f0b967260894571161142b19a8ae28 Mon Sep 17 00:00:00 2001 From: lukelowry Date: Thu, 3 Sep 2026 20:20:25 -0500 Subject: [PATCH 1/2] initial mu configuration configurability --- GridKit/CommonMath.hpp | 8 ++++- .../PhasorDynamics/Controller/REECB/Reecb.hpp | 5 --- .../Controller/REECB/ReecbImpl.hpp | 8 +++-- .../PhasorDynamics/AnalysisUtilities.hpp | 26 ++++++++++++-- .../PhasorDynamics/ContingencyAnalysis.cpp | 2 ++ .../PhasorDynamics/DynamicSimulation.cpp | 2 ++ application/PhasorDynamics/README.md | 1 + .../Math/SmoothnessIndicatorTests.hpp | 8 +++++ .../PhasorDynamics/ControllerReecbTests.hpp | 28 ++++++++------- .../PhasorDynamics/ConverterRegcaTests.hpp | 35 +++++++++++-------- 10 files changed, 84 insertions(+), 39 deletions(-) diff --git a/GridKit/CommonMath.hpp b/GridKit/CommonMath.hpp index 8f9f8c5cf..19209f94e 100644 --- a/GridKit/CommonMath.hpp +++ b/GridKit/CommonMath.hpp @@ -10,16 +10,22 @@ namespace GridKit { namespace Math { + template + inline constexpr RealT DEFAULT_MU = 240.0; + /** * @brief Smoothing scale shared by CommonMath primitives * * Used by @ref sigmoid, @ref ramp, and functions composed from them to set * the width of smooth transitions. * + * @warning Process-wide setting; configure before constructing models or + * launching workers. + * * @tparam RealT - real data type */ template - inline constexpr RealT MU = 240.0; + inline RealT MU = DEFAULT_MU; /** * @brief Scaled sigmoid activation function diff --git a/GridKit/Model/PhasorDynamics/Controller/REECB/Reecb.hpp b/GridKit/Model/PhasorDynamics/Controller/REECB/Reecb.hpp index 9188a0920..6e298bcfc 100644 --- a/GridKit/Model/PhasorDynamics/Controller/REECB/Reecb.hpp +++ b/GridKit/Model/PhasorDynamics/Controller/REECB/Reecb.hpp @@ -138,11 +138,6 @@ namespace GridKit ScalarT* f); private: - /// Hinge width chosen so closed-circle leakage stays below the - /// initialization tolerance [p.u. current squared]. - static constexpr RealT CURRENT_CIRCLE_KNEE = - INITIALIZATION_TOLERANCE / Math::MU; - struct InitialPoint; struct InitialCurrentLimit diff --git a/GridKit/Model/PhasorDynamics/Controller/REECB/ReecbImpl.hpp b/GridKit/Model/PhasorDynamics/Controller/REECB/ReecbImpl.hpp index 6d39b88ec..061f6a4f8 100644 --- a/GridKit/Model/PhasorDynamics/Controller/REECB/ReecbImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Controller/REECB/ReecbImpl.hpp @@ -1081,7 +1081,8 @@ namespace GridKit Reecb::sqrtramp(ValueT x) { const RealT root_width = ONE / Math::MU; - const RealT knee = CURRENT_CIRCLE_KNEE; + // Keep closed-circle leakage below the initialization tolerance. + const RealT knee = INITIALIZATION_TOLERANCE / Math::MU; const ValueT absolute = std::abs(x); const ValueT normalizer = absolute + knee; @@ -1117,12 +1118,13 @@ namespace GridKit } const RealT mu = Math::MU; + const RealT knee = INITIALIZATION_TOLERANCE / mu; const RealT scaled_y = HALF * mu * y; const RealT hinged = y / mu * (scaled_y + std::hypot(scaled_y, ONE)); const RealT square = hinged - - QUARTER * CURRENT_CIRCLE_KNEE - * (CURRENT_CIRCLE_KNEE / hinged); + - QUARTER * knee + * (knee / hinged); return std::max(ZERO, square); } diff --git a/application/PhasorDynamics/AnalysisUtilities.hpp b/application/PhasorDynamics/AnalysisUtilities.hpp index 867a73834..2e904a2ea 100644 --- a/application/PhasorDynamics/AnalysisUtilities.hpp +++ b/application/PhasorDynamics/AnalysisUtilities.hpp @@ -1,16 +1,19 @@ #pragma once +#include #include #include #include #include #include +#include #include #include #include #include +#include #include #include #include @@ -60,6 +63,8 @@ namespace GridKit double rel_tol; /// absolute tolerance for the solver double abs_tol; + /// CommonMath smoothing scale + double mu{Math::DEFAULT_MU}; /// fixed solver time step size, or 0 for adaptive stepping double dt_fixed; /// maximum number of solver time steps, or 0 for the IDA default @@ -82,6 +87,18 @@ namespace GridKit SystemModelData<> model_data; }; + /** + * @brief Apply study-wide CommonMath configuration + * + * Must be called before constructing a system model or launching study + * workers. + */ + template + inline void configureCommonMath(const StudyData& study) + { + Math::MU = static_cast(study.mu); + } + using json = ::nlohmann::json; using Log = ::GridKit::Utilities::Logger; @@ -99,8 +116,13 @@ namespace GridKit j.at("system_model_file").get_to(c.system_model_file); c.dt_monitor = j.value("dt_monitor", 0.0); j.at("tmax").get_to(c.tmax); - c.rel_tol = j.value("rel_tol", DEFAULT_SOLVER_REL_TOL); - c.abs_tol = j.value("abs_tol", DEFAULT_SOLVER_ABS_TOL); + c.rel_tol = j.value("rel_tol", DEFAULT_SOLVER_REL_TOL); + c.abs_tol = j.value("abs_tol", DEFAULT_SOLVER_ABS_TOL); + c.mu = j.value("mu", Math::DEFAULT_MU); + if (!std::isfinite(c.mu) || c.mu <= 0.0) + { + throw std::invalid_argument("\"mu\" must be a positive finite number"); + } c.dt_fixed = j.value("dt_fixed", 0.0); c.max_steps = j.value("max_steps", std::size_t{0}); c.consistent_ic_type = AnalysisManager::Sundials::IdaConsistentICType::YA_YDP; diff --git a/application/PhasorDynamics/ContingencyAnalysis.cpp b/application/PhasorDynamics/ContingencyAnalysis.cpp index 01d20821a..f5f8a55fe 100644 --- a/application/PhasorDynamics/ContingencyAnalysis.cpp +++ b/application/PhasorDynamics/ContingencyAnalysis.cpp @@ -154,6 +154,8 @@ int main(int argc, const char* argv[]) checkCommandLine(argc, "ContingencyAnalysis"); auto study_data = parseStudyData(argv[1]); + configureCommonMath(study_data); + const auto start = Clock::now(); auto faults = study_data.model_data.bus_fault; diff --git a/application/PhasorDynamics/DynamicSimulation.cpp b/application/PhasorDynamics/DynamicSimulation.cpp index 52c59e228..a376737c6 100644 --- a/application/PhasorDynamics/DynamicSimulation.cpp +++ b/application/PhasorDynamics/DynamicSimulation.cpp @@ -23,6 +23,8 @@ int main(int argc, const char* argv[]) checkCommandLine(argc, "DynamicSimulation"); auto study = parseStudyData(argv[1]); + configureCommonMath(study); + // Instantiate system SystemModel sys(study.model_data); sys.allocate(); diff --git a/application/PhasorDynamics/README.md b/application/PhasorDynamics/README.md index 969979fe0..f872fd923 100644 --- a/application/PhasorDynamics/README.md +++ b/application/PhasorDynamics/README.md @@ -9,6 +9,7 @@ `tmax` | A floating-point value for max time `rel_tol` | Relative solver tolerance (default: 1.0e-7) `abs_tol` | Absolute solver tolerance override (default: 1.0e-9) + `mu` | CommonMath smoothing scale; must be positive and finite (default: 240.0) `dt_fixed` | Fixed solver time step size, or 0 for adaptive stepping (default: 0) `max_steps` | Maximum number of solver time steps, 0 for the IDA default, or a negative number for unlimited steps (default: 0) `consistent_ic_type` | IDA consistent initial condition calculation type; one of { "y", "ya_ydp" } (default: "ya_ydp") diff --git a/tests/UnitTests/Math/SmoothnessIndicatorTests.hpp b/tests/UnitTests/Math/SmoothnessIndicatorTests.hpp index caaae8ebc..072d06b68 100644 --- a/tests/UnitTests/Math/SmoothnessIndicatorTests.hpp +++ b/tests/UnitTests/Math/SmoothnessIndicatorTests.hpp @@ -225,6 +225,14 @@ namespace GridKit success *= std::isfinite(Math::ramp(far_below)); success *= (Math::ramp(far_below) < roundoff); + // Verify ramp uses the configured runtime smoothing scale. + const RealT original_mu = Math::MU; + const RealT configured_mu = 500.0; + + Math::MU = configured_mu; + success *= within(Math::ramp(scalar(0.0)), std::log(scalar(2.0)) / scalar(configured_mu), roundoff); + Math::MU = original_mu; + success *= (smooth_clip > lower); success *= (smooth_clip < upper); success *= std::isfinite(Math::clamp(far_above, lower, upper)); diff --git a/tests/UnitTests/PhasorDynamics/ControllerReecbTests.hpp b/tests/UnitTests/PhasorDynamics/ControllerReecbTests.hpp index a4c5d34e9..a552b34ae 100644 --- a/tests/UnitTests/PhasorDynamics/ControllerReecbTests.hpp +++ b/tests/UnitTests/PhasorDynamics/ControllerReecbTests.hpp @@ -314,7 +314,7 @@ namespace GridKit success *= stateMatches(fixture.reecb, {{Vars::ILCAP, 2.0}}, "initial current-circle capacity", - kCircleTol); + circleTolerance()); const auto* initial_values = fixture.reecb.y().getData(); success *= scalarMatches(initial_values[index(Vars::IQMAX)], @@ -378,10 +378,10 @@ namespace GridKit {{Vars::PMEAS, 0.75}, {Vars::PORD, 1.5}}, "omitted component rating"); - success *= stateMatches(system_base.reecb, - {{Vars::ILCAP, 2.0}}, + success *= stateMatches(system_base.reecb, + {{Vars::ILCAP, 2.0}}, "omitted-rating current-circle capacity", - kCircleTol); + circleTolerance()); success *= allResidualsWithinInitTolerance(system_base.reecb); return success.report(__func__); @@ -710,7 +710,7 @@ namespace GridKit success *= stateMatches(boundary.reecb, {{Vars::ILCAP, test_case.capacity}}, test_case.label, - kCircleTol); + circleTolerance()); success *= allResidualsWithinInitTolerance(boundary.reecb); } @@ -784,7 +784,7 @@ namespace GridKit success *= stateMatches(exhausted.reecb, {{Vars::ILCAP, 0.0}}, "injection does not expand current circle", - kCircleTol); + circleTolerance()); success *= allResidualsWithinInitTolerance(exhausted.reecb); Log::setVerbosity(previous_verbosity); @@ -846,7 +846,7 @@ namespace GridKit RealT tolerance = kTol; if (variable == Vars::ILCAP) { - tolerance = kCircleTol; + tolerance = circleTolerance(); } success *= variableMatches(residuals[index(variable)], @@ -914,7 +914,7 @@ namespace GridKit success *= stateMatches(fixture.reecb, {{Vars::ILCAP, 2.0}}, "selector ILCAP", - kCircleTol); + circleTolerance()); // Exactly one reactive path carries the operating point. const auto* y = fixture.reecb.y().getData(); @@ -1545,7 +1545,7 @@ namespace GridKit success *= residualsMatch(fixture.reecb, {{Vars::ILCAP, circleLeg(square) - capacity_state}}, label, - kCircleTol); + circleTolerance()); } } @@ -1596,7 +1596,7 @@ namespace GridKit success *= residualsMatch(open.reecb, {{Vars::ILCAP, open_limit}}, "finite open current circle", - kCircleTol); + circleTolerance()); success *= allResidualsFinite(open.reecb); } @@ -1635,7 +1635,7 @@ namespace GridKit success *= residualsMatch(fixture.reecb, {{Vars::ILCAP, ideal_capacity}}, "off-axis capacity", - kCircleTol); + circleTolerance()); const RealT capacity = fixture.reecb.getResidual().getData()[index(Vars::ILCAP)]; if (!std::isfinite(capacity) || capacity < ZERO) @@ -1977,8 +1977,10 @@ namespace GridKit } /// MU-aware tolerance for ideal current-circle comparisons. - static constexpr RealT kCircleTol = - std::max(ONE / (Math::MU * Math::MU), kTol); + static RealT circleTolerance() + { + return std::max(ONE / (Math::MU * Math::MU), kTol); + } static constexpr size_t kBusVrColumn = index(Vars::MAXIMUM); diff --git a/tests/UnitTests/PhasorDynamics/ConverterRegcaTests.hpp b/tests/UnitTests/PhasorDynamics/ConverterRegcaTests.hpp index d2524f824..fd2ddf269 100644 --- a/tests/UnitTests/PhasorDynamics/ConverterRegcaTests.hpp +++ b/tests/UnitTests/PhasorDynamics/ConverterRegcaTests.hpp @@ -38,8 +38,11 @@ namespace GridKit static constexpr ScalarT kTol = static_cast(100.0) * std::numeric_limits::epsilon(); - static constexpr RealT kSmoothTol = - std::max(ONE / (Math::MU * Math::MU), kTol); + + static RealT smoothTolerance() + { + return std::max(ONE / (Math::MU * Math::MU), kTol); + } /// Construction, the monitor, and every verify() error class: missing /// and invalid parameters, a null bus, and an unlinked command port. @@ -187,8 +190,8 @@ namespace GridKit // The latched commands restore both displaced states at their ideal // interior first-order rates. const auto* f = latched.regca.getResidual().getData(); - success *= scalarMatches(f[index(Vars::IP)], 0.5, "latched active-current rate", kSmoothTol); - success *= scalarMatches(f[index(Vars::IQ)], 0.3, "latched reactive-current rate", kSmoothTol); + success *= scalarMatches(f[index(Vars::IP)], 0.5, "latched active-current rate", smoothTolerance()); + success *= scalarMatches(f[index(Vars::IQ)], 0.3, "latched reactive-current rate", smoothTolerance()); return success.report(__func__); } @@ -378,7 +381,7 @@ namespace GridKit const auto* f = residual.getData(); for (const auto& row : expected) { - success *= scalarMatches(f[index(row.variable)], row.value, row.name, kSmoothTol); + success *= scalarMatches(f[index(row.variable)], row.value, row.name, smoothTolerance()); } return success.report(__func__); @@ -577,7 +580,7 @@ namespace GridKit // Initialization above Vhvmax activates HVRCM and preserves Q0. { - const RealT terminal_voltage = kHvrcmVoltageLimit + kHvrcmOffset; + const RealT terminal_voltage = kHvrcmVoltageLimit + hvrcmOffset(); auto data = makeData(); data.parameters[Params::q0] = 0.1; @@ -589,8 +592,8 @@ namespace GridKit const auto* y = fixture.regca.y().getData(); const RealT extra_current = y[index(Vars::IQEXTRA)]; - success *= scalarMatches( - extra_current, kHvrcmGain * kHvrcmOffset, "IQEXTRA with the default Khv", kSmoothTol); + success *= scalarMatches( + extra_current, kHvrcmGain * hvrcmOffset(), "IQEXTRA with the default Khv", smoothTolerance()); success *= scalarMatches( y[index(Vars::IQ)] - y[index(Vars::IQEXTRA)], 0.1 / terminal_voltage, "IQ preserves Q0 after HVRCM compensation"); success *= scalarMatches(y[index(Vars::QBR)], 0.1, "QBR"); @@ -630,12 +633,12 @@ namespace GridKit return fixture.regca.getResidual().getData()[index(Vars::IQEXTRA)]; }; - const RealT below = residualAt(kHvrcmVoltageLimit - kHvrcmOffset, 0.0); + const RealT below = residualAt(kHvrcmVoltageLimit - hvrcmOffset(), 0.0); const RealT at = residualAt(kHvrcmVoltageLimit, 0.0); - const RealT above = residualAt(kHvrcmVoltageLimit + kHvrcmOffset, 0.0); + const RealT above = residualAt(kHvrcmVoltageLimit + hvrcmOffset(), 0.0); success *= scalarMatches(at, kHvrcmGain * hvrcmTransition(), "HVRCM residual at the threshold"); - success *= scalarMatches(above - below, kHvrcmGain * kHvrcmOffset, "HVRCM symmetric residual difference"); + success *= scalarMatches(above - below, kHvrcmGain * hvrcmOffset(), "HVRCM symmetric residual difference"); const RealT shifted = residualAt(kHvrcmVoltageLimit, 0.1); success *= scalarMatches(shifted - at, -0.1, "HVRCM extra-current residual shift"); @@ -832,11 +835,13 @@ namespace GridKit static constexpr RealT kHvrcmGain = 0.7; /// Offset whose smooth-ramp tail decays by one significand width. - static constexpr RealT kHvrcmOffset = - static_cast(std::numeric_limits::digits) - * std::numbers::ln2_v / Math::MU; + static RealT hvrcmOffset() + { + return static_cast(std::numeric_limits::digits) + * std::numbers::ln2_v / Math::MU; + } - static constexpr RealT hvrcmTransition() + static RealT hvrcmTransition() { return std::numbers::ln2_v / Math::MU; } From dd96f515b1b2cd9106339c75b8f51460d742b642 Mon Sep 17 00:00:00 2001 From: lukelowry Date: Sun, 6 Sep 2026 07:08:45 +0000 Subject: [PATCH 2/2] Apply pre-commit fixes --- tests/UnitTests/PhasorDynamics/ControllerReecbTests.hpp | 4 ++-- tests/UnitTests/PhasorDynamics/ConverterRegcaTests.hpp | 6 +++--- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/tests/UnitTests/PhasorDynamics/ControllerReecbTests.hpp b/tests/UnitTests/PhasorDynamics/ControllerReecbTests.hpp index a552b34ae..41dea041d 100644 --- a/tests/UnitTests/PhasorDynamics/ControllerReecbTests.hpp +++ b/tests/UnitTests/PhasorDynamics/ControllerReecbTests.hpp @@ -378,8 +378,8 @@ namespace GridKit {{Vars::PMEAS, 0.75}, {Vars::PORD, 1.5}}, "omitted component rating"); - success *= stateMatches(system_base.reecb, - {{Vars::ILCAP, 2.0}}, + success *= stateMatches(system_base.reecb, + {{Vars::ILCAP, 2.0}}, "omitted-rating current-circle capacity", circleTolerance()); success *= allResidualsWithinInitTolerance(system_base.reecb); diff --git a/tests/UnitTests/PhasorDynamics/ConverterRegcaTests.hpp b/tests/UnitTests/PhasorDynamics/ConverterRegcaTests.hpp index fd2ddf269..903725b1a 100644 --- a/tests/UnitTests/PhasorDynamics/ConverterRegcaTests.hpp +++ b/tests/UnitTests/PhasorDynamics/ConverterRegcaTests.hpp @@ -190,8 +190,8 @@ namespace GridKit // The latched commands restore both displaced states at their ideal // interior first-order rates. const auto* f = latched.regca.getResidual().getData(); - success *= scalarMatches(f[index(Vars::IP)], 0.5, "latched active-current rate", smoothTolerance()); - success *= scalarMatches(f[index(Vars::IQ)], 0.3, "latched reactive-current rate", smoothTolerance()); + success *= scalarMatches(f[index(Vars::IP)], 0.5, "latched active-current rate", smoothTolerance()); + success *= scalarMatches(f[index(Vars::IQ)], 0.3, "latched reactive-current rate", smoothTolerance()); return success.report(__func__); } @@ -592,7 +592,7 @@ namespace GridKit const auto* y = fixture.regca.y().getData(); const RealT extra_current = y[index(Vars::IQEXTRA)]; - success *= scalarMatches( + success *= scalarMatches( extra_current, kHvrcmGain * hvrcmOffset(), "IQEXTRA with the default Khv", smoothTolerance()); success *= scalarMatches( y[index(Vars::IQ)] - y[index(Vars::IQEXTRA)], 0.1 / terminal_voltage, "IQ preserves Q0 after HVRCM compensation");