diff --git a/CHANGELOG.md b/CHANGELOG.md index 59e9c884f..270ec7077 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -80,6 +80,7 @@ - Added `HYGOV` governor model implementation for PhasorDynamics. - Added `REPCA` controller model implementation for PhasorDynamics. - Added `REECB` electrical-control model implementation for PhasorDynamics. +- Added IDA option to choose the consistent initial condition calculation type. ## v0.1 diff --git a/GridKit/Solver/Dynamic/Ida.cpp b/GridKit/Solver/Dynamic/Ida.cpp index f2f2d5446..ffeedebb9 100644 --- a/GridKit/Solver/Dynamic/Ida.cpp +++ b/GridKit/Solver/Dynamic/Ida.cpp @@ -238,12 +238,8 @@ namespace AnalysisManager // Find a consistent set of initial conditions for DAE if (findConsistent) { - int initType = IDA_Y_INIT; - - if (tag_) - initType = IDA_YA_YDP_INIT; - - retval = IDACalcIC(solver_, initType, t0 + 0.1); + const int consistentICType = getIDAConsistentICType(); + retval = IDACalcIC(solver_, consistentICType, t0 + 0.1); checkOutput(retval, "IDACalcIC"); retval = IDAGetConsistentIC(solver_, yy_, yp_); @@ -307,6 +303,24 @@ namespace AnalysisManager model_->updateTime(t, 0.0); } + template + int Ida::getIDAConsistentICType() const + { + switch (consistent_ic_type_) + { + case IdaConsistentICType::Y: + return IDA_Y_INIT; + case IdaConsistentICType::YA_YDP: + return IDA_YA_YDP_INIT; + default: + GridKit::Utilities::Logger::error() + << "Invalid IDA consistent initial condition type " + << static_cast(consistent_ic_type_) + << ".\n"; + throw SundialsException(); + } + } + /** * @brief Run the IDA solver and optionally produce monitor output every `dt_monitor`. * @@ -1197,6 +1211,19 @@ namespace AnalysisManager backward_suppress_alg_ = suppress; } + /** + * @brief Set the IDA consistent-initial-condition calculation type + * + * @param consistent_ic_type IDA consistent initial condition type. + * @tparam ScalarT Scalar data type + * @tparam IdxT Index data type + */ + template + void Ida::setConsistentICType(IdaConsistentICType consistent_ic_type) + { + consistent_ic_type_ = consistent_ic_type; + } + /** * @brief Set the maximum number of steps * diff --git a/GridKit/Solver/Dynamic/Ida.hpp b/GridKit/Solver/Dynamic/Ida.hpp index bbe87f8df..bbda253ee 100644 --- a/GridKit/Solver/Dynamic/Ida.hpp +++ b/GridKit/Solver/Dynamic/Ida.hpp @@ -40,6 +40,12 @@ namespace AnalysisManager std::string report() const; }; + enum class IdaConsistentICType + { + Y, + YA_YDP + }; + template class Ida : public DynamicSolver { @@ -133,6 +139,7 @@ namespace AnalysisManager ScalarT abs_tol_override = 0); void setSuppressAlgebraicErrors(bool suppress); void setBackwardSuppressAlgebraicErrors(bool suppress); + void setConsistentICType(IdaConsistentICType consistent_ic_type); void setMaxSteps(IdxT maxSteps) override; void setBackwardMaxSteps(IdxT maxSteps); @@ -180,6 +187,7 @@ namespace AnalysisManager int getMonitorStepCount(RealT tf, RealT dt_monitor) const; RealT getMonitorTime(RealT tf, RealT dt_monitor, int step, int nsteps) const; + int getIDAConsistentICType() const; void updateModelState(RealT t); private: @@ -208,11 +216,12 @@ namespace AnalysisManager int backwardID_{}; - RealT time_step_{}; - RealT rel_tol_{DEFAULT_REL_TOL}; - RealT abs_tol_override_{}; - IdxT max_steps_{}; - bool suppress_alg_{false}; + RealT time_step_{}; + RealT rel_tol_{DEFAULT_REL_TOL}; + RealT abs_tol_override_{}; + IdxT max_steps_{}; + bool suppress_alg_{false}; + IdaConsistentICType consistent_ic_type_{IdaConsistentICType::YA_YDP}; RealT backward_time_step_{}; RealT backward_rel_tol_{DEFAULT_REL_TOL}; diff --git a/application/PhasorDynamics/AnalysisUtilities.hpp b/application/PhasorDynamics/AnalysisUtilities.hpp index bb30d2861..867a73834 100644 --- a/application/PhasorDynamics/AnalysisUtilities.hpp +++ b/application/PhasorDynamics/AnalysisUtilities.hpp @@ -12,6 +12,7 @@ #include #include +#include #include #include @@ -50,33 +51,35 @@ namespace GridKit struct StudyData { /// path to system model JSON file - fs::path system_model_file; + fs::path system_model_file; /// monitor output time step size, or 0 for no intermediate monitoring - double dt_monitor; + double dt_monitor; /// max time - double tmax; + double tmax; /// relative tolerance for the solver - double rel_tol; + double rel_tol; /// absolute tolerance for the solver - double abs_tol; + double abs_tol; /// fixed solver time step size, or 0 for adaptive stepping - double dt_fixed; + double dt_fixed; /// maximum number of solver time steps, or 0 for the IDA default - std::size_t max_steps; + std::size_t max_steps; + /// IDA consistent initial condition calculation type + AnalysisManager::Sundials::IdaConsistentICType consistent_ic_type; /// set of system events - std::vector events; + std::vector events; /// path to output file - fs::path output_file; + fs::path output_file; /// path to reference file for validation - fs::path reference_file; + fs::path reference_file; /// Error tolerance (between output file and reference file) - std::vector error_tol; + std::vector error_tol; /// Type of total error (relative or absolute) - Testing::ErrorType error_type; + Testing::ErrorType error_type; /// Smallest value at which to scale for relative error - double abs_err_threshold; + double abs_err_threshold; /// Instance of model data - SystemModelData<> model_data; + SystemModelData<> model_data; }; using json = ::nlohmann::json; @@ -96,10 +99,29 @@ 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.dt_fixed = j.value("dt_fixed", 0.0); - c.max_steps = j.value("max_steps", std::size_t{0}); + c.rel_tol = j.value("rel_tol", DEFAULT_SOLVER_REL_TOL); + c.abs_tol = j.value("abs_tol", DEFAULT_SOLVER_ABS_TOL); + 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; + if (j.contains("consistent_ic_type")) + { + const auto consistent_ic_type_str = j.at("consistent_ic_type").get(); + if (consistent_ic_type_str == "y") + { + c.consistent_ic_type = AnalysisManager::Sundials::IdaConsistentICType::Y; + } + else if (consistent_ic_type_str == "ya_ydp") + { + c.consistent_ic_type = AnalysisManager::Sundials::IdaConsistentICType::YA_YDP; + } + else + { + Log::error() << "Invalid IDA consistent initial condition type \"" + << consistent_ic_type_str << "\"; " + << "must be either \"y\" or \"ya_ydp\""; + } + } for (auto& raw_event : j.at("events")) { diff --git a/application/PhasorDynamics/ContingencyAnalysis.cpp b/application/PhasorDynamics/ContingencyAnalysis.cpp index 31415c6e2..01d20821a 100644 --- a/application/PhasorDynamics/ContingencyAnalysis.cpp +++ b/application/PhasorDynamics/ContingencyAnalysis.cpp @@ -38,6 +38,7 @@ TestStatus runStudy(StudyData study_data) ida.setTolerance(study_data.rel_tol, study_data.abs_tol); ida.setFixedStep(study_data.dt_fixed); ida.setMaxSteps(study_data.max_steps); + ida.setConsistentICType(study_data.consistent_ic_type); ida.configureSimulation(); using EventType = SystemEvent::Type; diff --git a/application/PhasorDynamics/DynamicSimulation.cpp b/application/PhasorDynamics/DynamicSimulation.cpp index c9f8f32ea..52c59e228 100644 --- a/application/PhasorDynamics/DynamicSimulation.cpp +++ b/application/PhasorDynamics/DynamicSimulation.cpp @@ -32,6 +32,7 @@ int main(int argc, const char* argv[]) ida.setTolerance(study.rel_tol, study.abs_tol); ida.setFixedStep(study.dt_fixed); ida.setMaxSteps(study.max_steps); + ida.setConsistentICType(study.consistent_ic_type); ida.configureSimulation(); // Start timer diff --git a/application/PhasorDynamics/README.md b/application/PhasorDynamics/README.md index 35ab99a5d..969979fe0 100644 --- a/application/PhasorDynamics/README.md +++ b/application/PhasorDynamics/README.md @@ -2,21 +2,22 @@ ## Root elements - Name | Value - ---------------------|------------------------------------------------------- - `system_model_file` | Path to the system model file[^1] - `dt_monitor` | Monitor output time interval for recorded simulation results (default: 0, no intermediate monitoring) - `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) - `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) - `events` | An array of event groups (see [Events](#events) below) - `output_file` | Path to output (CSV) file (optional) - `reference_file` | A string containing the name of the case (optional) - `error_type` | One of { "relative" (default), "absolute" } - `error_tolerance` | A floating-point value for highest allowable total error (default: 1.0e-4) - `abs_err_threshold` | A floating-point value for the smallest value at which to scale relative error (default: machine epsilon for double-precision) + Name | Value + ----------------------|------------------------------------------------------- + `system_model_file` | Path to the system model file[^1] + `dt_monitor` | Monitor output time interval for recorded simulation results (default: 0, no intermediate monitoring) + `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) + `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") + `events` | An array of event groups (see [Events](#events) below) + `output_file` | Path to output (CSV) file (optional) + `reference_file` | A string containing the name of the case (optional) + `error_type` | One of { "relative" (default), "absolute" } + `error_tolerance` | A floating-point value for highest allowable total error (default: 1.0e-4) + `abs_err_threshold` | A floating-point value for the smallest value at which to scale relative error (default: machine epsilon for double-precision) [^1]: See system model [case format](../../GridKit/Model/PhasorDynamics/INPUT_FORMAT.md) diff --git a/tests/UnitTests/Solver/Dynamic/IdaTests.hpp b/tests/UnitTests/Solver/Dynamic/IdaTests.hpp index 655f83afa..aea551aa5 100644 --- a/tests/UnitTests/Solver/Dynamic/IdaTests.hpp +++ b/tests/UnitTests/Solver/Dynamic/IdaTests.hpp @@ -369,6 +369,74 @@ namespace GridKit private: RealT t_{}; }; + + template + class ConsistentICTypeEvaluator : public NullEvaluator + { + protected: + using NullEvaluator::allocated_; + using NullEvaluator::y_; + using NullEvaluator::yp_; + using NullEvaluator::abs_tol_; + using NullEvaluator::tag_; + using NullEvaluator::f_; + + public: + ConsistentICTypeEvaluator() = default; + + explicit ConsistentICTypeEvaluator(bool steady_state) + : steady_state_(steady_state) + { + } + + int initialize() override + { + if (!allocated_) + { + this->allocate(); + } + + auto* y = y_.getData(); + auto* yp = yp_.getData(); + auto* abs_tol = abs_tol_.getData(); + auto* f = f_.getData(); + + y[0] = 0.0; + y[1] = 10.0; + yp[0] = steady_state_ ? 0.0 : 2.0; // Purposefully inconsistent guess for non-steady-state case. + yp[1] = 0.0; + tag_ = {true, false}; + abs_tol[0] = 0.0; + abs_tol[1] = 0.0; + f[0] = 0.0; + f[1] = 0.0; + y_.setDataUpdated(); + yp_.setDataUpdated(); + abs_tol_.setDataUpdated(); + f_.setDataUpdated(); + return 0; + } + + IdxT size() override + { + return 2; + } + + int evaluateResidual() override + { + auto* f = f_.getData(); + const auto* y = y_.getData(); + const auto* yp = yp_.getData(); + + f[0] = yp[0] + y[0] + y[1] - 1.0; + f[1] = y[1] - y[0]; + f_.setDataUpdated(); + return 0; + } + + private: + bool steady_state_{false}; + }; } // namespace Model namespace Testing @@ -503,6 +571,62 @@ namespace GridKit return success.report(__func__); } + + TestOutcome consistentICType() + { + TestStatus success = true; + + using RealT = typename ScalarTraits::RealT; + + // If the tolerances are too tight, a finite difference approximation to the Jacobian + // cannot be generated, and initialization will fail with bad initial guesses. + static constexpr auto tol = 100.0 * std::numeric_limits::epsilon(); + + // IdaConsistentICType::YA_YDP with non-steady-state initial guess for the derivatives + { + Model::ConsistentICTypeEvaluator model(false); + + Ida ida(&model); + ida.setConsistentICType(AnalysisManager::Sundials::IdaConsistentICType::YA_YDP); + ida.setTolerance(tol); + ida.configureSimulation(); + ida.initializeSimulation(0.0); + + success *= isEqual(model.yp().getData()[0], 1.0, tol); + success *= isEqual(model.yp().getData()[1], 0.0, tol); + success *= isEqual(model.y().getData()[0], 0.0, tol); + success *= isEqual(model.y().getData()[1], 0.0, tol); + } + + // IdaConsistentICType::Y with steady-state initial guess for the derivatives + { + Model::ConsistentICTypeEvaluator model(true); + + Ida ida(&model); + ida.setConsistentICType(AnalysisManager::Sundials::IdaConsistentICType::Y); + ida.setTolerance(tol); + ida.configureSimulation(); + ida.initializeSimulation(0.0); + + success *= isEqual(model.yp().getData()[0], 0.0, tol); + success *= isEqual(model.yp().getData()[1], 0.0, tol); + success *= isEqual(model.y().getData()[0], 0.5, tol); + success *= isEqual(model.y().getData()[1], 0.5, tol); + } + + { + Model::NullEvaluator model; + + Ida ida(&model); + ida.setConsistentICType(AnalysisManager::Sundials::IdaConsistentICType::Y); + ida.configureSimulation(); + ida.initializeSimulation(0.0); + + success *= isEqual(model.y().getData()[0], 0.0); + } + + return success.report(__func__); + } }; } // namespace Testing } // namespace GridKit diff --git a/tests/UnitTests/Solver/Dynamic/runIdaTests.cpp b/tests/UnitTests/Solver/Dynamic/runIdaTests.cpp index e4f25a296..6b127ce1a 100644 --- a/tests/UnitTests/Solver/Dynamic/runIdaTests.cpp +++ b/tests/UnitTests/Solver/Dynamic/runIdaTests.cpp @@ -13,6 +13,7 @@ int main() result += test.dtMonitorSuppressesEpsilonFinalStep(); result += test.fixedStep(); result += test.suppressAlgebraicErrors(); + result += test.consistentICType(); return result.summary(); }