Experimental OPF Implementation

CHANGELOG.md

@@ -10,6 +10,7 @@
 - Added `enable` and `disable` solver events that put a branch in or out of service.
 - `verify()` returns `Model::ConfigurationChecks` listing every configuration problem instead of an error count.
 - Added `Model::ParameterReader`
+- Added the `OptimalPowerFlow` model family, the `SteadyState::IpoptInterface`, and the `EconomicDispatch` example.
 - Removed legacy PhasorDynamics component constructors in favor of constructors taking model data.
 - Added hierarchical PhasorDynamics system models.
 - Migrated PhasorDynamics models with real and Boolean parameters to `Model::ParameterReader`.

GridKit/AutomaticDifferentiation/Enzyme/D2LDx2.hpp

@@ -0,0 +1,131 @@
+/**
+ * @file D2LDx2.hpp
+ * @brief Enzyme sparse Hessian of an optimization Lagrangian.
+ */
+
+#pragma once
+
+#include <array>
+#include <cstddef>
+
+#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 sparse Hessian of \f$L = \sigma f + \lambda^T g\f$, by one
+       * forward sweep per local variable over the reverse gradient of \f$L\f$
+       *
+       * @tparam ModelT - optimization component model
+       */
+      template <typename ModelT>
+      struct D2LDx2
+      {
+        /**
+         * @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 - Multipliers of the local constraint contributions
+         * @param[out] rows - Row indices
+         * @param[out] cols - Column indices
+         * @param[out] vals - Values
+         * @param[out] nnz - Number of nonzeros
+         */
+        template <typename ScalarT, typename IdxT, typename RealT>
+        static void eval(const ModelT*  model,
+                         const size_t   n_var,
+                         const IdxT*    variable_indices,
+                         const ScalarT* x,
+                         const RealT    sigma,
+                         const RealT*   lambda,
+                         IdxT*          rows,
+                         IdxT*          cols,
+                         ScalarT*       vals,
+                         IdxT&          nnz)
+        {
+          std::array<ScalarT, ObjectiveT::NUM_INPUTS> 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*) sparse_load<ScalarT, IdxT>,
+                                                            (void*) sparse_store<ScalarT, IdxT>,
+                                                            var_i,
+                                                            1.0, // value scaling
+                                                            variable_indices,
+                                                            variable_indices,
+                                                            rows,
+                                                            cols,
+                                                            vals,
+                                                            &nnz);
+
+            __enzyme_fwddiff<void>((void*) lagrangianGradient<ScalarT, RealT>,
+                                   enzyme_const,
+                                   model,
+                                   enzyme_dup,
+                                   x,
+                                   seed,
+                                   enzyme_dupnoneed,
+                                   gradient.data(),
+                                   dgradient,
+                                   enzyme_const,
+                                   sigma,
+                                   enzyme_const,
+                                   lambda);
+          }
+        }
+
+      private:
+        using ObjectiveT   = ModelWrapper<ModelT, MemberFunctions::Objective>;
+        using ConstraintsT = ModelWrapper<ModelT, MemberFunctions::Constraints>;
+
+        template <typename ScalarT, typename RealT>
+        [[gnu::always_inline]] inline static ScalarT lagrangian(const ModelT*  model,
+                                                                const ScalarT* x,
+                                                                const RealT    sigma,
+                                                                const RealT*   lambda)
+        {
+          std::array<ScalarT, ConstraintsT::NUM_OUTPUTS> g{};
+          ConstraintsT::eval(model, x, g.data());
+
+          ScalarT value = sigma * ObjectiveT::eval(model, x);
+          for (size_t i = 0; i < ConstraintsT::NUM_OUTPUTS; ++i)
+          {
+            value += lambda[i] * g[i];
+          }
+          return value;
+        }
+
+        /// Add the Lagrangian gradient to `gradient`, whose tangent alone is used
+        template <typename ScalarT, typename RealT>
+        [[gnu::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<ScalarT, RealT>,
+                                  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,50 @@
+/**
+ * @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 objective gradient, df/dx, by one reverse sweep
+       *
+       * @tparam ModelT - optimization component model
+       */
+      template <typename ModelT>
+      struct DfDx
+      {
+        /**
+         * @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
+         */
+        template <typename ScalarT, typename RealT>
+        static void eval(const ModelT* model, const ScalarT* x, RealT* gradient)
+        {
+          using WrapperT = ModelWrapper<ModelT, MemberFunctions::Objective>;
+
+          std::fill_n(gradient, WrapperT::NUM_INPUTS, 0.0);
+
+          __enzyme_autodiff<void>((void*) WrapperT::template eval<ScalarT>,
+                                  enzyme_const,
+                                  model,
+                                  enzyme_dup,
+                                  x,
+                                  gradient);
+        }
+      };
+    } // namespace Sparse
+  } // namespace Enzyme
+} // namespace GridKit

GridKit/AutomaticDifferentiation/Enzyme/DgDx.hpp

@@ -0,0 +1,86 @@
+/**
+ * @file DgDx.hpp
+ * @brief Enzyme sparse Jacobian of optimization constraints.
+ */
+
+#pragma once
+
+#include <array>
+#include <cstddef>
+
+#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 sparse constraint Jacobian, dg/dx, by one forward sweep
+       * per local variable
+       *
+       * @tparam ModelT - optimization component model
+       */
+      template <typename ModelT>
+      struct DgDx
+      {
+        /**
+         * @param[in] model - Pointer to the model to be differentiated
+         * @param[in] n_var - Number of local variables
+         * @param[in] row_indices - Global row of each local constraint contribution
+         * @param[in] col_indices - Global column of each local variable
+         * @param[in] x - Local variables
+         * @param[out] rows - Row indices
+         * @param[out] cols - Column indices
+         * @param[out] vals - Values
+         * @param[out] nnz - Number of nonzeros
+         */
+        template <typename ScalarT, typename IdxT>
+        static void eval(const ModelT*  model,
+                         const size_t   n_var,
+                         const IdxT*    row_indices,
+                         const IdxT*    col_indices,
+                         const ScalarT* x,
+                         IdxT*          rows,
+                         IdxT*          cols,
+                         ScalarT*       vals,
+                         IdxT&          nnz)
+        {
+          using WrapperT = ModelWrapper<ModelT, MemberFunctions::Constraints>;
+
+          std::array<ScalarT, WrapperT::NUM_OUTPUTS> 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*) sparse_load<ScalarT, IdxT>,
+                                                     (void*) sparse_store<ScalarT, IdxT>,
+                                                     var_i,
+                                                     1.0, // value scaling
+                                                     row_indices,
+                                                     col_indices,
+                                                     rows,
+                                                     cols,
+                                                     vals,
+                                                     &nnz);
+
+            __enzyme_fwddiff<void>((void*) WrapperT::template eval<ScalarT>,
+                                   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

GridKit/AutomaticDifferentiation/Enzyme/LowerSparseStorage.hpp

@@ -245,13 +245,15 @@ namespace GridKit
       }
 
       /**
-       * @brief Enzyme sparse load
+       * @brief Enzyme sparse load, with the arguments of `sparse_store`
+       *
+       * Returns zero, so an accumulation stores its increment.
        *
        * @tparam ScalarT - scalar data type
        * @tparam IdxT - matrix index data type
        */
       template <typename ScalarT, typename IdxT>
-      __attribute__((always_inline)) static ScalarT sparse_load(IdxT, IdxT, IdxT*, IdxT*, ScalarT*)
+      [[gnu::always_inline]] static ScalarT sparse_load(IdxT, IdxT, ScalarT, const IdxT*, const IdxT*, IdxT*, IdxT*, ScalarT*, IdxT*)
       {
         return 0.0;
       }

GridKit/AutomaticDifferentiation/Enzyme/ModelWrappers.hpp

@@ -6,9 +6,11 @@
 
 #pragma once
 
+#include <cstddef>
 #include <vector>
 
 #include <GridKit/Model/EvaluatorTraits.hpp>
+#include <GridKit/Utilities/Enum.hpp>
 
 namespace GridKit
 {
@@ -28,7 +30,9 @@ namespace GridKit
         BusResidual11, //< Special case for branches that are connected to two buses
         BusResidual12, //< Special case for branches that are connected to two buses
         BusResidual21, //< Special case for branches that are connected to two buses
-        BusResidual22  //< Special case for branches that are connected to two buses
+        BusResidual22, //< Special case for branches that are connected to two buses
+        Objective,     //< Optimization objective
+        Constraints    //< Optimization internal constraints and contributions at the ports
       };
 
       /**
@@ -227,6 +231,58 @@ namespace GridKit
           model->evaluateBusResidual22(y, yp, wb, h);
         }
       };
+
+      /**
+       * @brief Wrapper partial template specialization for an optimization objective
+       */
+      template <typename ModelT>
+      struct ModelWrapper<ModelT, MemberFunctions::Objective>
+      {
+        /// Number of internal variables
+        static constexpr std::size_t NUM_INTERNAL_INPUTS = Utilities::enum_size<typename ModelT::InternalVariablesT>();
+
+        /// Number of local variables
+        static constexpr std::size_t NUM_INPUTS = NUM_INTERNAL_INPUTS + Utilities::enum_size<typename ModelT::ExternalVariablesT>();
+
+        /**
+         * @param[in] model - Pointer to the model to be differentiated
+         * @param[in] x - Local variables
+         * @return Objective contribution
+         */
+        template <typename ScalarT>
+        [[gnu::always_inline]] inline static ScalarT eval(const ModelT* model, const ScalarT* x)
+        {
+          return model->objective(x, x + NUM_INTERNAL_INPUTS);
+        }
+      };
+
+      /**
+       * @brief Wrapper partial template specialization for optimization constraints
+       */
+      template <typename ModelT>
+      struct ModelWrapper<ModelT, MemberFunctions::Constraints>
+      {
+        /// Number of internal variables
+        static constexpr std::size_t NUM_INTERNAL_INPUTS = Utilities::enum_size<typename ModelT::InternalVariablesT>();
+
+        /// Number of internal constraints
+        static constexpr std::size_t NUM_INTERNAL_OUTPUTS = Utilities::enum_size<typename ModelT::InternalConstraintsT>();
+
+        /// Number of local constraint contributions
+        static constexpr std::size_t NUM_OUTPUTS = NUM_INTERNAL_OUTPUTS + Utilities::enum_size<typename ModelT::ExternalConstraintsT>();
+
+        /**
+         * @param[in] model - Pointer to the model to be differentiated
+         * @param[in] x - Local variables
+         * @param[out] g - Local constraint contributions
+         */
+        template <typename ScalarT>
+        [[gnu::always_inline]] inline static void eval(const ModelT* model, const ScalarT* x, ScalarT* g)
+        {
+          model->internalConstraints(x, x + NUM_INTERNAL_INPUTS, g);
+          model->externalConstraints(x, x + NUM_INTERNAL_INPUTS, g + NUM_INTERNAL_OUTPUTS);
+        }
+      };
     } // namespace Sparse
   } // namespace Enzyme
 } // namespace GridKit