Experimental OPF Implementation using MATPOWER Input

CHANGELOG.md

@@ -95,6 +95,8 @@
 - Added `GridKit_ENABLE_DEVELOPER_MODE` CMake option; in developer mode the Logger defaults to `EVERYTHING` verbosity.
 - Added `Logger::raiseVerbosity()`, which raises the verbosity without lowering a higher level.
 - Changed most of the core library code outputs to use the Logger instead of `std::cout`/`std::cerr`.
+- Added the `StateData` reader and writer for the state file format, with `applyState` and `extractState` for PhasorDynamics and a `state_file` option in phasor dynamics solver JSON files.
+- Added the `OptimalPowerFlow` model family and the `EconomicDispatch` application, with sparse Enzyme Jacobians and Hessians and an Ipopt interface, to initialize PhasorDynamics at an optimal operating point with the limits and costs of a MATPOWER case.
 
 ## v0.1
 

CMakeLists.txt

@@ -126,7 +126,7 @@ endif("${isSystemDir}" STREQUAL "-1")
 option(BUILD_SHARED_LIBS "Build shared libraries" ON)
 
 if(GRIDKIT_ENABLE_IPOPT)
-  include(FindIpopt)
+  find_package(Ipopt REQUIRED)
   enable_language(Fortran) # Needed for linking to HSL
 endif()
 
@@ -217,7 +217,10 @@ configure_package_config_file(
   INSTALL_DESTINATION share/cmake/gridkit
   NO_SET_AND_CHECK_MACRO NO_CHECK_REQUIRED_COMPONENTS_MACRO)
 
+# Make the dependency finder available beside the generated package config.
+configure_file(cmake/FindIpopt.cmake ${CMAKE_CURRENT_BINARY_DIR}/FindIpopt.cmake COPYONLY)
+
 # Install configuration file
 install(FILES ${CMAKE_CURRENT_BINARY_DIR}/GridKitConfig.cmake
               ${CMAKE_CURRENT_BINARY_DIR}/GridKitConfigVersion.cmake
-        DESTINATION share/cmake/gridkit)
+              ${CMAKE_CURRENT_BINARY_DIR}/FindIpopt.cmake DESTINATION share/cmake/gridkit)

GridKit/AutomaticDifferentiation/CMakeLists.txt

@@ -1,3 +1,4 @@
 install(FILES DependencyTracking/Variable.hpp DependencyTracking/VariableImplementation.hpp
               DependencyTracking/VariableOperators.hpp
         DESTINATION include/GridKit/AutomaticDifferentiation/DependencyTracking)
+install(FILES Enzyme/CooEntries.hpp DESTINATION include/GridKit/AutomaticDifferentiation/Enzyme)

GridKit/AutomaticDifferentiation/Enzyme/CooEntries.hpp

@@ -0,0 +1,38 @@
+/**
+ * @file CooEntries.hpp
+ * @brief Coordinate entries filled by Enzyme sparse stores.
+ */
+
+#pragma once
+
+#include <cstddef>
+#include <vector>
+
+namespace GridKit
+{
+  namespace Enzyme
+  {
+    namespace Sparse
+    {
+      /**
+       * @brief Coordinate entries in global indices
+       *
+       * Filled by `mapped_store` and `lower_store`. Assembly sums duplicate
+       * coordinates.
+       */
+      struct CooEntries
+      {
+        std::vector<size_t> rows;
+        std::vector<size_t> cols;
+        std::vector<double> values;
+
+        void clear()
+        {
+          rows.clear();
+          cols.clear();
+          values.clear();
+        }
+      };
+    } // namespace Sparse
+  } // namespace Enzyme
+} // namespace GridKit

GridKit/AutomaticDifferentiation/Enzyme/D2LDx2.hpp

@@ -0,0 +1,132 @@
+/**
+ * @file D2LDx2.hpp
+ * @brief Enzyme sparse Hessian of an optimization Lagrangian.
+ */
+
+#pragma once
+
+#include <array>
+#include <cstddef>
+#include <type_traits>
+
+#include <GridKit/AutomaticDifferentiation/Enzyme/EnzymeDefinitions.hpp>
+#include <GridKit/AutomaticDifferentiation/Enzyme/LowerSparseStorage.hpp>
+#include <GridKit/AutomaticDifferentiation/Enzyme/ModelWrappers.hpp>
+
+namespace GridKit
+{
+  namespace Enzyme
+  {
+    namespace Sparse
+    {
+      /**
+       * @brief Enzyme automatic differentiation Hessian evaluator: Lagrangian Hessian, d2L/dx2
+       *
+       * Lower triangle of the Hessian of
+       * \f$L = \sigma f(x) + \lambda^T g(x)\f$ by one forward sweep per local
+       * variable over the reverse-mode gradient of \f$L\f$. Only structural
+       * entries are stored, in global indices.
+       *
+       * @tparam ModelT - model type with `VARIABLE_SIZE` local variables and
+       * `CONSTRAINT_SIZE` local constraints
+       */
+      template <typename ModelT>
+      struct D2LDx2
+      {
+        using ScalarT = typename ModelT::ScalarT;
+        using IdxT    = typename ModelT::IdxT;
+        using RealT   = typename ModelT::RealT;
+
+        static_assert(std::is_same_v<ScalarT, double> && std::is_same_v<IdxT, size_t>,
+                      "D2LDx2 supports double values and size_t indices");
+
+        /**
+         * @param[in] model - Pointer to the model to be differentiated
+         * @param[in] n_var - Number of local variables
+         * @param[in] variable_indices - Global index of each local variable
+         * @param[in] x - Local variables
+         * @param[in] sigma - Objective factor
+         * @param[in] lambda - Local constraint multipliers
+         * @param[in,out] entries - Hessian entries
+         */
+        static void eval(const ModelT*  model,
+                         const size_t   n_var,
+                         const IdxT*    variable_indices,
+                         const ScalarT* x,
+                         const RealT    sigma,
+                         const RealT*   lambda,
+                         CooEntries*    entries)
+        {
+          std::array<ScalarT, ModelT::VARIABLE_SIZE> gradient{};
+          for (size_t var_i = 0; var_i < n_var; ++var_i)
+          {
+            // Sparse storage. @see LowerSparseStorage.hpp
+            ScalarT* seed      = __enzyme_todense<ScalarT*>((void*) ident_load<ScalarT, IdxT>,
+                                                       (void*) ident_store<ScalarT, IdxT>,
+                                                       var_i);
+            ScalarT* dgradient = __enzyme_todense<ScalarT*>((void*) mapped_load,
+                                                            (void*) lower_store,
+                                                            var_i,
+                                                            variable_indices,
+                                                            variable_indices,
+                                                            entries);
+
+            __enzyme_fwddiff<void>((void*) lagrangianGradient,
+                                   enzyme_const,
+                                   model,
+                                   enzyme_dup,
+                                   x,
+                                   seed,
+                                   enzyme_dupnoneed,
+                                   gradient.data(),
+                                   dgradient,
+                                   enzyme_const,
+                                   sigma,
+                                   enzyme_const,
+                                   lambda);
+          }
+        }
+
+      private:
+        __attribute__((always_inline)) inline static ScalarT lagrangian(const ModelT*  model,
+                                                                        const ScalarT* x,
+                                                                        const RealT    sigma,
+                                                                        const RealT*   lambda)
+        {
+          std::array<ScalarT, ModelT::CONSTRAINT_SIZE> g{};
+          ModelWrapper<ModelT, MemberFunctions::Constraints>::eval(model, x, g.data());
+
+          ScalarT value = sigma * ModelWrapper<ModelT, MemberFunctions::Objective>::eval(model, x);
+          for (size_t i = 0; i < ModelT::CONSTRAINT_SIZE; ++i)
+          {
+            value += lambda[i] * g[i];
+          }
+          return value;
+        }
+
+        /**
+         * @brief Add the gradient of the Lagrangian to `gradient`
+         *
+         * Only the tangent of `gradient` is used, so it is not reset.
+         */
+        __attribute__((always_inline)) inline static void lagrangianGradient(const ModelT*  model,
+                                                                             const ScalarT* x,
+                                                                             ScalarT*       gradient,
+                                                                             const RealT    sigma,
+                                                                             const RealT*   lambda)
+        {
+          __enzyme_autodiff<void>((void*) lagrangian,
+                                  enzyme_const,
+                                  model,
+                                  enzyme_dup,
+                                  x,
+                                  gradient,
+                                  enzyme_const,
+                                  sigma,
+                                  enzyme_const,
+                                  lambda);
+        }
+      };
+    } // namespace Sparse
+  } // namespace Enzyme
+} // namespace GridKit

GridKit/AutomaticDifferentiation/Enzyme/DfDx.hpp

@@ -0,0 +1,52 @@
+/**
+ * @file DfDx.hpp
+ * @brief Enzyme gradient of an optimization objective.
+ */
+
+#pragma once
+
+#include <algorithm>
+#include <cstddef>
+
+#include <GridKit/AutomaticDifferentiation/Enzyme/EnzymeDefinitions.hpp>
+#include <GridKit/AutomaticDifferentiation/Enzyme/ModelWrappers.hpp>
+
+namespace GridKit
+{
+  namespace Enzyme
+  {
+    namespace Sparse
+    {
+      /**
+       * @brief Enzyme automatic differentiation gradient evaluator: objective gradient, df/dx
+       *
+       * One reverse sweep.
+       *
+       * @tparam ModelT - model type with `VARIABLE_SIZE` local variables
+       */
+      template <typename ModelT>
+      struct DfDx
+      {
+        using ScalarT = typename ModelT::ScalarT;
+        using RealT   = typename ModelT::RealT;
+
+        /**
+         * @param[in] model - Pointer to the model to be differentiated
+         * @param[in] x - Local variables
+         * @param[out] gradient - Gradient with respect to the local variables
+         */
+        static void eval(const ModelT* model, const ScalarT* x, RealT* gradient)
+        {
+          std::fill_n(gradient, ModelT::VARIABLE_SIZE, 0.0);
+
+          __enzyme_autodiff<void>((void*) ModelWrapper<ModelT, MemberFunctions::Objective>::eval,
+                                  enzyme_const,
+                                  model,
+                                  enzyme_dup,
+                                  x,
+                                  gradient);
+        }
+      };
+    } // namespace Sparse
+  } // namespace Enzyme
+} // namespace GridKit

GridKit/AutomaticDifferentiation/Enzyme/DgDx.hpp

@@ -0,0 +1,82 @@
+/**
+ * @file DgDx.hpp
+ * @brief Enzyme sparse Jacobian of optimization constraints.
+ */
+
+#pragma once
+
+#include <array>
+#include <cstddef>
+#include <type_traits>
+
+#include <GridKit/AutomaticDifferentiation/Enzyme/EnzymeDefinitions.hpp>
+#include <GridKit/AutomaticDifferentiation/Enzyme/LowerSparseStorage.hpp>
+#include <GridKit/AutomaticDifferentiation/Enzyme/ModelWrappers.hpp>
+
+namespace GridKit
+{
+  namespace Enzyme
+  {
+    namespace Sparse
+    {
+      /**
+       * @brief Enzyme automatic differentiation Jacobian evaluator: constraint Jacobian, dg/dx
+       *
+       * One forward sweep per local variable. Only structural entries are
+       * stored, in global indices.
+       *
+       * @tparam ModelT - model type with `CONSTRAINT_SIZE` local constraints
+       */
+      template <typename ModelT>
+      struct DgDx
+      {
+        using ScalarT = typename ModelT::ScalarT;
+        using IdxT    = typename ModelT::IdxT;
+
+        static_assert(std::is_same_v<ScalarT, double> && std::is_same_v<IdxT, size_t>,
+                      "DgDx supports double values and size_t indices");
+
+        /**
+         * @param[in] model - Pointer to the model to be differentiated
+         * @param[in] n_var - Number of local variables
+         * @param[in] constraint_indices - Global row of each local constraint
+         * @param[in] variable_indices - Global column of each local variable
+         * @param[in] x - Local variables
+         * @param[in,out] entries - Jacobian entries
+         */
+        static void eval(const ModelT*  model,
+                         const size_t   n_var,
+                         const IdxT*    constraint_indices,
+                         const IdxT*    variable_indices,
+                         const ScalarT* x,
+                         CooEntries*    entries)
+        {
+          std::array<ScalarT, ModelT::CONSTRAINT_SIZE> g{};
+          for (size_t var_i = 0; var_i < n_var; ++var_i)
+          {
+            // Sparse storage. @see LowerSparseStorage.hpp
+            ScalarT* seed = __enzyme_todense<ScalarT*>((void*) ident_load<ScalarT, IdxT>,
+                                                       (void*) ident_store<ScalarT, IdxT>,
+                                                       var_i);
+            ScalarT* dg   = __enzyme_todense<ScalarT*>((void*) mapped_load,
+                                                     (void*) mapped_store,
+                                                     var_i,
+                                                     constraint_indices,
+                                                     variable_indices,
+                                                     entries);
+
+            __enzyme_fwddiff<void>((void*) ModelWrapper<ModelT, MemberFunctions::Constraints>::eval,
+                                   enzyme_const,
+                                   model,
+                                   enzyme_dup,
+                                   x,
+                                   seed,
+                                   enzyme_dupnoneed,
+                                   g.data(),
+                                   dg);
+          }
+        }
+      };
+    } // namespace Sparse
+  } // namespace Enzyme
+} // namespace GridKit

GridKit/AutomaticDifferentiation/Enzyme/EnzymeDefinitions.hpp

@@ -48,6 +48,15 @@ namespace GridKit
        */
       template <typename T, typename... model_type>
       extern T __enzyme_fwddiff(void*, model_type...) noexcept;
+
+      /**
+       * @brief Enzyme reverse-mode template for GridKit models
+       *
+       * @tparam T - return type
+       * @tparam model_type - model type
+       */
+      template <typename T, typename... model_type>
+      extern T __enzyme_autodiff(void*, model_type...) noexcept;
     } // namespace Sparse
   } // namespace Enzyme
 } // namespace GridKit