Add option for IDA consistent initial conditions mode

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
 

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 <class ScalarT, typename IdxT>
+    int Ida<ScalarT, IdxT>::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<int>(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 <class ScalarT, typename IdxT>
+    void Ida<ScalarT, IdxT>::setConsistentICType(IdaConsistentICType consistent_ic_type)
+    {
+      consistent_ic_type_ = consistent_ic_type;
+    }
+
     /**
      * @brief Set the maximum number of steps
      *

GridKit/Solver/Dynamic/Ida.hpp

@@ -40,6 +40,12 @@ namespace AnalysisManager
       std::string report() const;
     };
 
+    enum class IdaConsistentICType
+    {
+      Y,
+      YA_YDP
+    };
+
     template <class ScalarT, typename IdxT>
     class Ida : public DynamicSolver<ScalarT, IdxT>
     {
@@ -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};

application/PhasorDynamics/AnalysisUtilities.hpp

@@ -12,6 +12,7 @@
 #include <nlohmann/json.hpp>
 
 #include <GridKit/Model/PhasorDynamics/SystemModelData.hpp>
+#include <GridKit/Solver/Dynamic/Ida.hpp>
 #include <GridKit/Testing/TestHelpers.hpp>
 #include <GridKit/Utilities/Logger/Logger.hpp>
 
@@ -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<SystemEvent> events;
+      std::vector<SystemEvent>                       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<double>      error_tol;
+      std::vector<double>                            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<std::string>();
+        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"))
       {

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;

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

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)
 

tests/UnitTests/Solver/Dynamic/IdaTests.hpp

@@ -369,6 +369,74 @@ namespace GridKit
     private:
       RealT t_{};
     };
+
+    template <class ScalarT, typename IdxT>
+    class ConsistentICTypeEvaluator : public NullEvaluator<ScalarT, IdxT>
+    {
+    protected:
+      using NullEvaluator<ScalarT, IdxT>::allocated_;
+      using NullEvaluator<ScalarT, IdxT>::y_;
+      using NullEvaluator<ScalarT, IdxT>::yp_;
+      using NullEvaluator<ScalarT, IdxT>::abs_tol_;
+      using NullEvaluator<ScalarT, IdxT>::tag_;
+      using NullEvaluator<ScalarT, IdxT>::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<ScalarT>::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<RealT>::epsilon();
+
+        // IdaConsistentICType::YA_YDP with non-steady-state initial guess for the derivatives
+        {
+          Model::ConsistentICTypeEvaluator<ScalarT, IdxT> model(false);
+
+          Ida<ScalarT, IdxT> 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<ScalarT, IdxT> model(true);
+
+          Ida<ScalarT, IdxT> 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<ScalarT, IdxT> model;
+
+          Ida<ScalarT, IdxT> 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