diff --git a/CHANGELOG.md b/CHANGELOG.md
index 59e9c884f..a7238d9ca 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 subsystem partitioning support for `PowerElectronics` models, including partition interfaces, and independent residual and Jacobian evaluation.
 
 ## v0.1
 
diff --git a/GridKit/Model/PowerElectronics/CMakeLists.txt b/GridKit/Model/PowerElectronics/CMakeLists.txt
index c985d69d6..368f89515 100644
--- a/GridKit/Model/PowerElectronics/CMakeLists.txt
+++ b/GridKit/Model/PowerElectronics/CMakeLists.txt
@@ -21,6 +21,7 @@ add_subdirectory(TransmissionLine)
 add_subdirectory(MicrogridLoad)
 add_subdirectory(MicrogridLine)
 add_subdirectory(MicrogridBusDQ)
+add_subdirectory(PartitionInterface)
 
 install(
   FILES CircuitComponent.hpp
@@ -29,4 +30,5 @@ install(
         SystemModelPowerElectronics.hpp
         NodeBase.hpp
         ExternalConnection.hpp
+        SubsystemModel.hpp
   DESTINATION include/GridKit/Model/PowerElectronics)
diff --git a/GridKit/Model/PowerElectronics/Capacitor/Capacitor.cpp b/GridKit/Model/PowerElectronics/Capacitor/Capacitor.cpp
index aab348679..c9c3d1150 100644
--- a/GridKit/Model/PowerElectronics/Capacitor/Capacitor.cpp
+++ b/GridKit/Model/PowerElectronics/Capacitor/Capacitor.cpp
@@ -52,6 +52,25 @@ namespace GridKit
     return 0;
   }
 
+  /**
+   * @brief Compute the absolute tolerance for each variable in the model
+   *
+   * @param rel_tol The relative tolerance which can be used to pick the
+   *        absolute tolerance.
+   * @tparam ScalarT Scalar data type
+   * @tparam IdxT Index data type
+   * @return int 0 if successful, non-zero otherwise.
+   *
+   * This represents a "noise" level close to zero for which pure relative
+   * error cannot be used.
+   */
+  template <class ScalarT, typename IdxT>
+  int Capacitor<ScalarT, IdxT>::setAbsoluteTolerance(RealT rel_tol)
+  {
+    abs_tol_.setToConst(static_cast<ScalarT>(rel_tol));
+    return 0;
+  }
+
   /**
    * @brief Evaluate the resisdual of the Capcitor
    *
@@ -117,6 +136,18 @@ namespace GridKit
     return 0;
   }
 
+  template <class ScalarT, typename IdxT>
+  bool Capacitor<ScalarT, IdxT>::isCloneable() const
+  {
+    return true;
+  }
+
+  template <class ScalarT, typename IdxT>
+  CircuitComponent<ScalarT, IdxT>* Capacitor<ScalarT, IdxT>::clone() const
+  {
+    return new Capacitor<ScalarT, IdxT>(*this);
+  }
+
   // Available template instantiations
   template class Capacitor<double, long int>;
   template class Capacitor<double, size_t>;
diff --git a/GridKit/Model/PowerElectronics/Capacitor/Capacitor.hpp b/GridKit/Model/PowerElectronics/Capacitor/Capacitor.hpp
index d7a131a9a..0b3353c47 100644
--- a/GridKit/Model/PowerElectronics/Capacitor/Capacitor.hpp
+++ b/GridKit/Model/PowerElectronics/Capacitor/Capacitor.hpp
@@ -29,6 +29,7 @@ namespace GridKit
     using CircuitComponent<ScalarT, IdxT>::y_int_;
     using CircuitComponent<ScalarT, IdxT>::yp_ext_;
     using CircuitComponent<ScalarT, IdxT>::yp_int_;
+    using CircuitComponent<ScalarT, IdxT>::abs_tol_;
     using CircuitComponent<ScalarT, IdxT>::tag_;
     using CircuitComponent<ScalarT, IdxT>::f_ext_;
     using CircuitComponent<ScalarT, IdxT>::f_int_;
@@ -50,15 +51,19 @@ namespace GridKit
 
     int initialize();
     int tagDifferentiable();
+    int setAbsoluteTolerance(RealT);
     int evaluateInternalResidual() final;
     int evaluateExternalResidual() final;
     int evaluateJacobian();
     int evaluateIntegrand();
 
-    int initializeAdjoint();
-    int evaluateAdjointResidual();
+    int  initializeAdjoint();
+    int  evaluateAdjointResidual();
     // int evaluateAdjointJacobian();
-    int evaluateAdjointIntegrand();
+    int  evaluateAdjointIntegrand();
+    bool isCloneable() const;
+
+    CircuitComponent<ScalarT, IdxT>* clone() const;
 
   private:
     RealT C_;
diff --git a/GridKit/Model/PowerElectronics/CircuitComponent.hpp b/GridKit/Model/PowerElectronics/CircuitComponent.hpp
index becbc4208..4dac238b1 100644
--- a/GridKit/Model/PowerElectronics/CircuitComponent.hpp
+++ b/GridKit/Model/PowerElectronics/CircuitComponent.hpp
@@ -27,6 +27,138 @@ namespace GridKit
 
     CircuitComponent() = default;
 
+    CircuitComponent(const CircuitComponent& other)
+      : n_extern_(other.n_extern_),
+        n_intern_(other.n_intern_),
+        extern_indices_(other.extern_indices_),
+        size_(other.size_),
+        nnz_(other.nnz_),
+        size_quad_(other.size_quad_),
+        size_opt_(other.size_opt_),
+        current_jac_size_(other.current_jac_size_),
+
+        // These pointers refer to storage supplied by a parent system.
+        // The copied component must be connected to its own storage later.
+        y_int_(nullptr),
+        yp_int_(nullptr),
+        f_int_(nullptr),
+
+        tag_(other.tag_),
+        time_(other.time_),
+        alpha_(other.alpha_),
+        max_steps_(other.max_steps_),
+        idc_(other.idc_),
+        allocated_(other.allocated_)
+    {
+      /*
+       * VectorT disables its normal copy constructor and copy-assignment
+       * operator. Use its provided copyFromExternal() operation to perform
+       * an independent copy of the vector data.
+       */
+      auto copyVector = [](VectorT& destination, const VectorT& source)
+      {
+        const IdxT source_size = source.getSize();
+
+        if (source_size == 0)
+        {
+          return;
+        }
+
+        destination.resize(source_size);
+        destination.copyFromExternal(source);
+      };
+
+      /*
+       * Deep-copy the local-to-global connection mapping.
+       */
+      if (other.connection_nodes_)
+      {
+        connection_nodes_ = std::make_unique<IdxT[]>(static_cast<size_t>(size_));
+
+        for (IdxT i = 0; i < size_; ++i)
+        {
+          connection_nodes_[static_cast<size_t>(i)] = other.connection_nodes_[static_cast<size_t>(i)];
+        }
+      }
+
+      /*
+       * Deep-copy the COO Jacobian row indices.
+       */
+      if (other.jacobian_coo_rows_)
+      {
+        jacobian_coo_rows_ = std::make_unique<IdxT[]>(static_cast<size_t>(nnz_));
+
+        for (IdxT i = 0; i < nnz_; ++i)
+        {
+          jacobian_coo_rows_[static_cast<size_t>(i)] = other.jacobian_coo_rows_[static_cast<size_t>(i)];
+        }
+      }
+
+      /*
+       * Deep-copy the COO Jacobian column indices.
+       */
+      if (other.jacobian_coo_cols_)
+      {
+        jacobian_coo_cols_ = std::make_unique<IdxT[]>(static_cast<size_t>(nnz_));
+
+        for (IdxT i = 0; i < nnz_; ++i)
+        {
+          jacobian_coo_cols_[static_cast<size_t>(i)] = other.jacobian_coo_cols_[static_cast<size_t>(i)];
+        }
+      }
+
+      /*
+       * Deep-copy the COO Jacobian values.
+       */
+      if (other.jacobian_coo_values_)
+      {
+        jacobian_coo_values_ = std::make_unique<RealT[]>(static_cast<size_t>(nnz_));
+
+        for (IdxT i = 0; i < nnz_; ++i)
+        {
+          jacobian_coo_values_[static_cast<size_t>(i)] = other.jacobian_coo_values_[static_cast<size_t>(i)];
+        }
+      }
+
+      if (size_ > 0)
+      {
+        y_ext_  = std::make_unique<const ScalarT*[]>(static_cast<size_t>(size_));
+        yp_ext_ = std::make_unique<const ScalarT*[]>(static_cast<size_t>(size_));
+        f_ext_  = std::make_unique<ScalarT*[]>(static_cast<size_t>(size_));
+
+        for (IdxT i = 0; i < size_; ++i)
+        {
+          y_ext_[i]  = nullptr;
+          yp_ext_[i] = nullptr;
+          f_ext_[i]  = nullptr;
+        }
+      }
+
+      // State, state derivative, residual, and absolute tolerance.
+      copyVector(y_, other.y_);
+      copyVector(yp_, other.yp_);
+      copyVector(f_, other.f_);
+      copyVector(abs_tol_, other.abs_tol_);
+      copyVector(g_, other.g_);
+      copyVector(yB_, other.yB_);
+      copyVector(ypB_, other.ypB_);
+      copyVector(fB_, other.fB_);
+      copyVector(gB_, other.gB_);
+      copyVector(param_, other.param_);
+      copyVector(param_up_, other.param_up_);
+      copyVector(param_lo_, other.param_lo_);
+    }
+
+    virtual CircuitComponent<ScalarT, IdxT>* clone() const
+    {
+      return nullptr;
+    }
+
+    virtual bool isCloneable() const
+    {
+      return false;
+    }
+
     /**
      * @note Cannot be marked final, since it is overriden to recurse in the system model.
      */
@@ -51,7 +183,7 @@ namespace GridKit
       return this->n_intern_;
     }
 
-    std::set<size_t> getExternIndices()
+    std::set<IdxT> getExternIndices()
     {
       return this->extern_indices_;
     }
@@ -69,7 +201,7 @@ namespace GridKit
     int setInternalConnectionNodes(size_t local_index, IdxT global_index)
     {
       assert(!extern_indices_.contains(static_cast<IdxT>(local_index)));
-      connection_nodes_[local_index] = global_index;
+      setConnectionNodes(local_index, global_index);
       return 0;
     }
 
@@ -95,6 +227,23 @@ namespace GridKit
       return 0;
     }
 
+    /**
+     * @brief Update the connection index for a variable.
+     *
+     * Changes only the connection index without modifying the variable's
+     * internal/external classification or its associated data pointers.
+     *
+     * @param local_index Index of the local variable.
+     * @param connection_index New connection index for the variable.
+     *
+     * @return int 0 if successful.
+     */
+    int setConnectionNodes(size_t local_index, IdxT connection_index)
+    {
+      connection_nodes_[local_index] = connection_index;
+      return 0;
+    }
+
     /**
      * @brief Given the location of value in the local vector map to global index
      *
@@ -132,9 +281,10 @@ namespace GridKit
       jacobian_coo_cols_   = std::make_unique<IdxT[]>(static_cast<size_t>(nnz_));
       jacobian_coo_values_ = std::make_unique<RealT[]>(static_cast<size_t>(nnz_));
 
-      y_ext_            = std::make_unique<const ScalarT*[]>(static_cast<size_t>(size_));
-      yp_ext_           = std::make_unique<const ScalarT*[]>(static_cast<size_t>(size_));
-      f_ext_            = std::make_unique<ScalarT*[]>(static_cast<size_t>(size_));
+      y_ext_  = std::make_unique<const ScalarT*[]>(static_cast<size_t>(size_));
+      yp_ext_ = std::make_unique<const ScalarT*[]>(static_cast<size_t>(size_));
+      f_ext_  = std::make_unique<ScalarT*[]>(static_cast<size_t>(size_));
+
       connection_nodes_ = std::make_unique<IdxT[]>(static_cast<size_t>(size_));
 
       if (!allocated_)
@@ -225,6 +375,51 @@ namespace GridKit
       f_int_ = internal_res;
     }
 
+    /**
+     * @brief Clear all internal and external data pointers.
+     *
+     * This disconnects the component from the data storage currently associated
+     * with its internal and external variables. The external pointer arrays
+     * themselves remain allocated so that the component can be connected to
+     * new data later.
+     *
+     * @return int 0 if successful.
+     */
+    int clearPointers()
+    {
+      // Clear internal data pointers.
+      y_int_  = nullptr;
+      yp_int_ = nullptr;
+      f_int_  = nullptr;
+
+      // Clear external data pointers, but keep the pointer arrays allocated.
+      if (y_ext_)
+      {
+        for (IdxT i = 0; i < size_; ++i)
+        {
+          y_ext_[i] = nullptr;
+        }
+      }
+
+      if (yp_ext_)
+      {
+        for (IdxT i = 0; i < size_; ++i)
+        {
+          yp_ext_[i] = nullptr;
+        }
+      }
+
+      if (f_ext_)
+      {
+        for (IdxT i = 0; i < size_; ++i)
+        {
+          f_ext_[i] = nullptr;
+        }
+      }
+
+      return 0;
+    }
+
   protected:
     /**
      * @brief Reset the Jacobian so it can be constructed. Helper method for \ref setJacValues().
@@ -439,6 +634,16 @@ namespace GridKit
       return idc_;
     }
 
+    /**
+     * @brief Check whether the component has already been allocated.
+     *
+     * @return true if allocate() has previously completed, false otherwise.
+     */
+    bool isAllocated() const
+    {
+      return allocated_;
+    }
+
   protected:
     /**
      * @brief Allocate state and residual storage owned by this component.
diff --git a/GridKit/Model/PowerElectronics/DistributedGenerator/DistributedGenerator.cpp b/GridKit/Model/PowerElectronics/DistributedGenerator/DistributedGenerator.cpp
index dfa6f1a98..bc12c6fa6 100644
--- a/GridKit/Model/PowerElectronics/DistributedGenerator/DistributedGenerator.cpp
+++ b/GridKit/Model/PowerElectronics/DistributedGenerator/DistributedGenerator.cpp
@@ -421,6 +421,18 @@ namespace GridKit
     return 0;
   }
 
+  template <class ScalarT, typename IdxT>
+  bool DistributedGenerator<ScalarT, IdxT>::isCloneable() const
+  {
+    return true;
+  }
+
+  template <class ScalarT, typename IdxT>
+  CircuitComponent<ScalarT, IdxT>* DistributedGenerator<ScalarT, IdxT>::clone() const
+  {
+    return new DistributedGenerator<ScalarT, IdxT>(*this);
+  }
+
   // Available template instantiations
   template class DistributedGenerator<double, long int>;
   template class DistributedGenerator<double, size_t>;
diff --git a/GridKit/Model/PowerElectronics/DistributedGenerator/DistributedGenerator.hpp b/GridKit/Model/PowerElectronics/DistributedGenerator/DistributedGenerator.hpp
index ab7e340f1..8a5219848 100644
--- a/GridKit/Model/PowerElectronics/DistributedGenerator/DistributedGenerator.hpp
+++ b/GridKit/Model/PowerElectronics/DistributedGenerator/DistributedGenerator.hpp
@@ -75,18 +75,21 @@ namespace GridKit
                          NodeT*                                      node_bus);
     virtual ~DistributedGenerator();
 
-    int initialize();
-    int allocate() final;
-    int tagDifferentiable();
-    int setAbsoluteTolerance(RealT);
-    int evaluateInternalResidual() final;
-    int evaluateExternalResidual() final;
-    int evaluateJacobian();
-    int evaluateIntegrand();
-    int initializeAdjoint();
-    int evaluateAdjointResidual();
+    int  initialize();
+    int  allocate() final;
+    int  tagDifferentiable();
+    int  setAbsoluteTolerance(RealT);
+    int  evaluateInternalResidual() final;
+    int  evaluateExternalResidual() final;
+    int  evaluateJacobian();
+    int  evaluateIntegrand();
+    int  initializeAdjoint();
+    int  evaluateAdjointResidual();
     // int evaluateAdjointJacobian();
-    int evaluateAdjointIntegrand();
+    int  evaluateAdjointIntegrand();
+    bool isCloneable() const;
+
+    CircuitComponent<ScalarT, IdxT>* clone() const;
 
   private:
     RealT wb_;
diff --git a/GridKit/Model/PowerElectronics/InductionMotor/InductionMotor.cpp b/GridKit/Model/PowerElectronics/InductionMotor/InductionMotor.cpp
index 0224e6ebe..a17456b65 100644
--- a/GridKit/Model/PowerElectronics/InductionMotor/InductionMotor.cpp
+++ b/GridKit/Model/PowerElectronics/InductionMotor/InductionMotor.cpp
@@ -63,6 +63,25 @@ namespace GridKit
     return 0;
   }
 
+  /**
+   * @brief Compute the absolute tolerance for each variable in the model
+   *
+   * @param rel_tol The relative tolerance which can be used to pick the
+   *        absolute tolerance.
+   * @tparam ScalarT Scalar data type
+   * @tparam IdxT Index data type
+   * @return int 0 if successful, non-zero otherwise.
+   *
+   * This represents a "noise" level close to zero for which pure relative
+   * error cannot be used.
+   */
+  template <class ScalarT, typename IdxT>
+  int InductionMotor<ScalarT, IdxT>::setAbsoluteTolerance(RealT rel_tol)
+  {
+    abs_tol_.setToConst(static_cast<ScalarT>(rel_tol));
+    return 0;
+  }
+
   /**
    * @brief Contributes to the resisdual
    *
@@ -129,6 +148,18 @@ namespace GridKit
     return 0;
   }
 
+  template <class ScalarT, typename IdxT>
+  bool InductionMotor<ScalarT, IdxT>::isCloneable() const
+  {
+    return true;
+  }
+
+  template <class ScalarT, typename IdxT>
+  CircuitComponent<ScalarT, IdxT>* InductionMotor<ScalarT, IdxT>::clone() const
+  {
+    return new InductionMotor<ScalarT, IdxT>(*this);
+  }
+
   // Available template instantiations
   template class InductionMotor<double, long int>;
   template class InductionMotor<double, size_t>;
diff --git a/GridKit/Model/PowerElectronics/InductionMotor/InductionMotor.hpp b/GridKit/Model/PowerElectronics/InductionMotor/InductionMotor.hpp
index 122cd4d30..6f1724670 100644
--- a/GridKit/Model/PowerElectronics/InductionMotor/InductionMotor.hpp
+++ b/GridKit/Model/PowerElectronics/InductionMotor/InductionMotor.hpp
@@ -29,6 +29,7 @@ namespace GridKit
     using CircuitComponent<ScalarT, IdxT>::y_int_;
     using CircuitComponent<ScalarT, IdxT>::yp_ext_;
     using CircuitComponent<ScalarT, IdxT>::yp_int_;
+    using CircuitComponent<ScalarT, IdxT>::abs_tol_;
     using CircuitComponent<ScalarT, IdxT>::tag_;
     using CircuitComponent<ScalarT, IdxT>::f_ext_;
     using CircuitComponent<ScalarT, IdxT>::f_int_;
@@ -50,15 +51,19 @@ namespace GridKit
 
     int initialize();
     int tagDifferentiable();
+    int setAbsoluteTolerance(RealT);
     int evaluateInternalResidual() final;
     int evaluateExternalResidual() final;
     int evaluateJacobian();
     int evaluateIntegrand();
 
-    int initializeAdjoint();
-    int evaluateAdjointResidual();
+    int  initializeAdjoint();
+    int  evaluateAdjointResidual();
     // int evaluateAdjointJacobian();
-    int evaluateAdjointIntegrand();
+    int  evaluateAdjointIntegrand();
+    bool isCloneable() const;
+
+    CircuitComponent<ScalarT, IdxT>* clone() const;
 
   private:
     RealT Lls_;
diff --git a/GridKit/Model/PowerElectronics/Inductor/Inductor.cpp b/GridKit/Model/PowerElectronics/Inductor/Inductor.cpp
index b192cebc5..f9e2f5d3d 100644
--- a/GridKit/Model/PowerElectronics/Inductor/Inductor.cpp
+++ b/GridKit/Model/PowerElectronics/Inductor/Inductor.cpp
@@ -147,6 +147,18 @@ namespace GridKit
     return 0;
   }
 
+  template <class ScalarT, typename IdxT>
+  bool Inductor<ScalarT, IdxT>::isCloneable() const
+  {
+    return true;
+  }
+
+  template <class ScalarT, typename IdxT>
+  CircuitComponent<ScalarT, IdxT>* Inductor<ScalarT, IdxT>::clone() const
+  {
+    return new Inductor<ScalarT, IdxT>(*this);
+  }
+
   // Available template instantiations
   template class Inductor<double, long int>;
   template class Inductor<double, size_t>;
diff --git a/GridKit/Model/PowerElectronics/Inductor/Inductor.hpp b/GridKit/Model/PowerElectronics/Inductor/Inductor.hpp
index b82ce1070..30d25db79 100644
--- a/GridKit/Model/PowerElectronics/Inductor/Inductor.hpp
+++ b/GridKit/Model/PowerElectronics/Inductor/Inductor.hpp
@@ -60,10 +60,13 @@ namespace GridKit
     int evaluateJacobian();
     int evaluateIntegrand();
 
-    int initializeAdjoint();
-    int evaluateAdjointResidual();
+    int  initializeAdjoint();
+    int  evaluateAdjointResidual();
     // int evaluateAdjointJacobian();
-    int evaluateAdjointIntegrand();
+    int  evaluateAdjointIntegrand();
+    bool isCloneable() const;
+
+    CircuitComponent<ScalarT, IdxT>* clone() const;
 
   private:
     RealT  L_;
diff --git a/GridKit/Model/PowerElectronics/LinearTransformer/LinearTransformer.cpp b/GridKit/Model/PowerElectronics/LinearTransformer/LinearTransformer.cpp
index fcfc4a02d..b1537c8ba 100644
--- a/GridKit/Model/PowerElectronics/LinearTransformer/LinearTransformer.cpp
+++ b/GridKit/Model/PowerElectronics/LinearTransformer/LinearTransformer.cpp
@@ -65,6 +65,25 @@ namespace GridKit
     return 0;
   }
 
+  /**
+   * @brief Compute the absolute tolerance for each variable in the model
+   *
+   * @param rel_tol The relative tolerance which can be used to pick the
+   *        absolute tolerance.
+   * @tparam ScalarT Scalar data type
+   * @tparam IdxT Index data type
+   * @return int 0 if successful, non-zero otherwise.
+   *
+   * This represents a "noise" level close to zero for which pure relative
+   * error cannot be used.
+   */
+  template <class ScalarT, typename IdxT>
+  int LinearTransformer<ScalarT, IdxT>::setAbsoluteTolerance(RealT rel_tol)
+  {
+    abs_tol_.setToConst(static_cast<ScalarT>(rel_tol));
+    return 0;
+  }
+
   /**
    * @brief Computes the component resisdual
    */
@@ -114,6 +133,18 @@ namespace GridKit
     return 0;
   }
 
+  template <class ScalarT, typename IdxT>
+  bool LinearTransformer<ScalarT, IdxT>::isCloneable() const
+  {
+    return true;
+  }
+
+  template <class ScalarT, typename IdxT>
+  CircuitComponent<ScalarT, IdxT>* LinearTransformer<ScalarT, IdxT>::clone() const
+  {
+    return new LinearTransformer<ScalarT, IdxT>(*this);
+  }
+
   // Available template instantiations
   template class LinearTransformer<double, long int>;
   template class LinearTransformer<double, size_t>;
diff --git a/GridKit/Model/PowerElectronics/LinearTransformer/LinearTransformer.hpp b/GridKit/Model/PowerElectronics/LinearTransformer/LinearTransformer.hpp
index 62dc84305..6b75dcc1d 100644
--- a/GridKit/Model/PowerElectronics/LinearTransformer/LinearTransformer.hpp
+++ b/GridKit/Model/PowerElectronics/LinearTransformer/LinearTransformer.hpp
@@ -29,6 +29,7 @@ namespace GridKit
     using CircuitComponent<ScalarT, IdxT>::y_int_;
     using CircuitComponent<ScalarT, IdxT>::yp_ext_;
     using CircuitComponent<ScalarT, IdxT>::yp_int_;
+    using CircuitComponent<ScalarT, IdxT>::abs_tol_;
     using CircuitComponent<ScalarT, IdxT>::tag_;
     using CircuitComponent<ScalarT, IdxT>::f_ext_;
     using CircuitComponent<ScalarT, IdxT>::f_int_;
@@ -50,15 +51,19 @@ namespace GridKit
 
     int initialize();
     int tagDifferentiable();
+    int setAbsoluteTolerance(RealT);
     int evaluateInternalResidual() final;
     int evaluateExternalResidual() final;
     int evaluateJacobian();
     int evaluateIntegrand();
 
-    int initializeAdjoint();
-    int evaluateAdjointResidual();
+    int  initializeAdjoint();
+    int  evaluateAdjointResidual();
     // int evaluateAdjointJacobian();
-    int evaluateAdjointIntegrand();
+    int  evaluateAdjointIntegrand();
+    bool isCloneable() const;
+
+    CircuitComponent<ScalarT, IdxT>* clone() const;
 
   private:
     RealT L0_;
diff --git a/GridKit/Model/PowerElectronics/MicrogridBusDQ/MicrogridBusDQ.cpp b/GridKit/Model/PowerElectronics/MicrogridBusDQ/MicrogridBusDQ.cpp
index ffc02dfb3..765286583 100644
--- a/GridKit/Model/PowerElectronics/MicrogridBusDQ/MicrogridBusDQ.cpp
+++ b/GridKit/Model/PowerElectronics/MicrogridBusDQ/MicrogridBusDQ.cpp
@@ -154,6 +154,18 @@ namespace GridKit
     return 0;
   }
 
+  template <class ScalarT, typename IdxT>
+  bool MicrogridBusDQ<ScalarT, IdxT>::isCloneable() const
+  {
+    return true;
+  }
+
+  template <class ScalarT, typename IdxT>
+  CircuitComponent<ScalarT, IdxT>* MicrogridBusDQ<ScalarT, IdxT>::clone() const
+  {
+    return new MicrogridBusDQ<ScalarT, IdxT>(*this);
+  }
+
   // Available template instantiations
   template class MicrogridBusDQ<double, long int>;
   template class MicrogridBusDQ<double, size_t>;
diff --git a/GridKit/Model/PowerElectronics/MicrogridBusDQ/MicrogridBusDQ.hpp b/GridKit/Model/PowerElectronics/MicrogridBusDQ/MicrogridBusDQ.hpp
index 111ca7be8..2314700a3 100644
--- a/GridKit/Model/PowerElectronics/MicrogridBusDQ/MicrogridBusDQ.hpp
+++ b/GridKit/Model/PowerElectronics/MicrogridBusDQ/MicrogridBusDQ.hpp
@@ -60,10 +60,13 @@ namespace GridKit
     int evaluateJacobian();
     int evaluateIntegrand();
 
-    int initializeAdjoint();
-    int evaluateAdjointResidual();
+    int  initializeAdjoint();
+    int  evaluateAdjointResidual();
     // int evaluateAdjointJacobian();
-    int evaluateAdjointIntegrand();
+    int  evaluateAdjointIntegrand();
+    bool isCloneable() const;
+
+    CircuitComponent<ScalarT, IdxT>* clone() const;
 
   private:
     RealT  RN_;
diff --git a/GridKit/Model/PowerElectronics/MicrogridLine/MicrogridLine.cpp b/GridKit/Model/PowerElectronics/MicrogridLine/MicrogridLine.cpp
index 65c3fb0a0..1292cd045 100644
--- a/GridKit/Model/PowerElectronics/MicrogridLine/MicrogridLine.cpp
+++ b/GridKit/Model/PowerElectronics/MicrogridLine/MicrogridLine.cpp
@@ -182,6 +182,18 @@ namespace GridKit
     return 0;
   }
 
+  template <class ScalarT, typename IdxT>
+  bool MicrogridLine<ScalarT, IdxT>::isCloneable() const
+  {
+    return true;
+  }
+
+  template <class ScalarT, typename IdxT>
+  CircuitComponent<ScalarT, IdxT>* MicrogridLine<ScalarT, IdxT>::clone() const
+  {
+    return new MicrogridLine<ScalarT, IdxT>(*this);
+  }
+
   // Available template instantiations
   template class MicrogridLine<double, long int>;
   template class MicrogridLine<double, size_t>;
diff --git a/GridKit/Model/PowerElectronics/MicrogridLine/MicrogridLine.hpp b/GridKit/Model/PowerElectronics/MicrogridLine/MicrogridLine.hpp
index 1354a9d9b..041562cbb 100644
--- a/GridKit/Model/PowerElectronics/MicrogridLine/MicrogridLine.hpp
+++ b/GridKit/Model/PowerElectronics/MicrogridLine/MicrogridLine.hpp
@@ -60,10 +60,13 @@ namespace GridKit
     int evaluateJacobian();
     int evaluateIntegrand();
 
-    int initializeAdjoint();
-    int evaluateAdjointResidual();
+    int  initializeAdjoint();
+    int  evaluateAdjointResidual();
     // int evaluateAdjointJacobian();
-    int evaluateAdjointIntegrand();
+    int  evaluateAdjointIntegrand();
+    bool isCloneable() const;
+
+    CircuitComponent<ScalarT, IdxT>* clone() const;
 
   private:
     RealT  R_;
diff --git a/GridKit/Model/PowerElectronics/MicrogridLoad/MicrogridLoad.cpp b/GridKit/Model/PowerElectronics/MicrogridLoad/MicrogridLoad.cpp
index a14a038f1..7739b39be 100644
--- a/GridKit/Model/PowerElectronics/MicrogridLoad/MicrogridLoad.cpp
+++ b/GridKit/Model/PowerElectronics/MicrogridLoad/MicrogridLoad.cpp
@@ -173,6 +173,18 @@ namespace GridKit
     return 0;
   }
 
+  template <class ScalarT, typename IdxT>
+  bool MicrogridLoad<ScalarT, IdxT>::isCloneable() const
+  {
+    return true;
+  }
+
+  template <class ScalarT, typename IdxT>
+  CircuitComponent<ScalarT, IdxT>* MicrogridLoad<ScalarT, IdxT>::clone() const
+  {
+    return new MicrogridLoad<ScalarT, IdxT>(*this);
+  }
+
   // Available template instantiations
   template class MicrogridLoad<double, long int>;
   template class MicrogridLoad<double, size_t>;
diff --git a/GridKit/Model/PowerElectronics/MicrogridLoad/MicrogridLoad.hpp b/GridKit/Model/PowerElectronics/MicrogridLoad/MicrogridLoad.hpp
index 931983d50..6189e941c 100644
--- a/GridKit/Model/PowerElectronics/MicrogridLoad/MicrogridLoad.hpp
+++ b/GridKit/Model/PowerElectronics/MicrogridLoad/MicrogridLoad.hpp
@@ -60,10 +60,13 @@ namespace GridKit
     int evaluateJacobian();
     int evaluateIntegrand();
 
-    int initializeAdjoint();
-    int evaluateAdjointResidual();
+    int  initializeAdjoint();
+    int  evaluateAdjointResidual();
     // int evaluateAdjointJacobian();
-    int evaluateAdjointIntegrand();
+    int  evaluateAdjointIntegrand();
+    bool isCloneable() const;
+
+    CircuitComponent<ScalarT, IdxT>* clone() const;
 
   private:
     RealT  R_;
diff --git a/GridKit/Model/PowerElectronics/NodeBase.hpp b/GridKit/Model/PowerElectronics/NodeBase.hpp
index 3c7f7c80b..11ddb1117 100644
--- a/GridKit/Model/PowerElectronics/NodeBase.hpp
+++ b/GridKit/Model/PowerElectronics/NodeBase.hpp
@@ -136,6 +136,23 @@ namespace GridKit
         };
       }
 
+      /**
+       * @brief Update the connection index for a variable.
+       *
+       * Changes only the connection index without modifying the variable's
+       * internal/external classification or its associated data pointers.
+       *
+       * @param local_index Index of the local variable.
+       * @param connection_index New connection index for the variable.
+       *
+       * @return int 0 if successful.
+       */
+      int setConnectionNodes(size_t local_index, IdxT connection_index)
+      {
+        connection_nodes_[local_index] = connection_index;
+        return 0;
+      }
+
       int allocate() override
       {
         size_t size = static_cast<size_t>(n_intern_ + n_extern_);
@@ -392,6 +409,16 @@ namespace GridKit
         return gB_;
       }
 
+      /**
+       * @brief Check whether the Node has already been allocated.
+       *
+       * @return true if allocate() has previously completed, false otherwise.
+       */
+      bool isAllocated() const
+      {
+        return allocated_;
+      }
+
     private:
       void allocateVectors(IdxT n)
       {
diff --git a/GridKit/Model/PowerElectronics/PartitionInterface/BusPartitionInterface.cpp b/GridKit/Model/PowerElectronics/PartitionInterface/BusPartitionInterface.cpp
new file mode 100644
index 000000000..98c5ac704
--- /dev/null
+++ b/GridKit/Model/PowerElectronics/PartitionInterface/BusPartitionInterface.cpp
@@ -0,0 +1,376 @@
+
+#include "BusPartitionInterface.hpp"
+
+#include <cassert>
+#include <cmath>
+#include <cstddef>
+
+#include <GridKit/Model/PowerElectronics/ExternalConnection.hpp>
+
+namespace GridKit
+{
+
+  /**
+   * @brief Construct a partition interface between a bus and a circuit component.
+   *
+   * The interface wraps a copy of a component located across a partition
+   * boundary and evaluates the contributions that the component makes to the
+   * connected bus. All interface variables are treated as external variables
+   * and have the same size as the wrapped component.
+   *
+   *
+   * @param bus Bus associated with the partition interface.
+   * @param component Copy of the component across the partition boundary.
+   * @param id Unique identifier for the interface.
+   */
+  template <class ScalarT, typename IdxT>
+  BusPartitionInterface<ScalarT, IdxT>::BusPartitionInterface(node_type* bus, component_type* component, IdxT id)
+    : component_(component->clone()),
+      bus_(bus)
+  {
+    size_     = component_->size();
+    n_intern_ = 0;
+    n_extern_ = static_cast<size_t>(component_->size());
+    idc_      = id;
+
+    // All variables of the bus interface are external to the interface.
+    for (IdxT i = 0; i < size_; i++)
+    {
+      extern_indices_.insert(i);
+    }
+
+    // Map each global bus connection index to its position in the bus.
+    std::unordered_map<IdxT, size_t> bus_connections;
+
+    for (size_t i = 0; i < static_cast<size_t>(bus_->size()); ++i)
+    {
+      bus_connections[bus_->getNodeConnection(i).idx_] = i;
+    }
+
+    // The interface Jacobian contains only rows corresponding to variables
+    // owned by the bus.
+    const IdxT* coo_rows = component_->jacobianCooRows();
+
+    nnz_ = 0;
+
+    for (IdxT k = 0; k < component_->nnz(); ++k)
+    {
+      const IdxT row_node = component_->getNodeConnection(static_cast<size_t>(coo_rows[k]));
+
+      if (bus_connections.contains(row_node))
+      {
+        ++nnz_;
+        jac_map_.push_back(k); // Keep track of the entries so they can be easily extracted
+      }
+    }
+  }
+
+  template <class ScalarT, typename IdxT>
+  BusPartitionInterface<ScalarT, IdxT>::~BusPartitionInterface()
+  {
+    delete component_;
+  }
+
+  /**
+   * @brief Allocate storage and initialize the interface mappings.
+   *
+   * Identifies the external variables of the wrapped component that are
+   * connected to the bus and builds the mappings used to transfer residual
+   * contributions between the component and the interface. Private storage is
+   * also allocated for the internal variables of the wrapped component.
+   *
+   * @return 0 on success and a nonzero value if the component does not contain
+   *         all variables required by the bus interface.
+   */
+  template <class ScalarT, typename IdxT>
+  int BusPartitionInterface<ScalarT, IdxT>::allocate()
+  {
+
+    CircuitComponent<ScalarT, IdxT>::allocate();
+
+    // Build a lookup of the global indices belonging to the bus.
+    std::unordered_map<IdxT, size_t> bus_connections;
+
+    for (size_t i = 0; i < static_cast<size_t>(bus_->size()); ++i)
+    {
+      bus_connections[bus_->getNodeConnection(i).idx_] = i;
+    }
+
+    const auto& external_indices = component_->getExternIndices();
+
+    is_external_.assign(static_cast<size_t>(size_), false);
+
+    bus_input_ports_.clear();
+    bus_output_ports_.clear();
+
+    size_t external_index = 0;
+
+    // Identify which component variables are external and which of those external
+    // variables are connected to the bus. bus_input_ports_ stores the corresponding
+    // variable position in this interface, while bus_output_ports_ stores its position
+    // in the wrapped component's external residual vector.
+    for (size_t i = 0; i < static_cast<size_t>(size_); ++i)
+    {
+      const IdxT connection_index = component_->getNodeConnection(i);
+
+      this->setConnectionNodes(i, connection_index);
+
+      if (!external_indices.contains(static_cast<IdxT>(i)))
+      {
+        continue;
+      }
+
+      is_external_[i] = true;
+
+      if (bus_connections.contains(connection_index))
+      {
+        bus_input_ports_.push_back(i);
+        bus_output_ports_.push_back(external_index);
+      }
+
+      ++external_index;
+    }
+
+    // A valid bus interface must contain every variable belonging to the bus.
+    // If fewer bus variables are found, the wrapped component does not provide
+    // the complete coupling required by this interface.
+    const size_t bus_size = static_cast<size_t>(bus_->size());
+
+    if (bus_input_ports_.size() != bus_size)
+    {
+      std::cerr << "ERROR: Invalid partition interface detected. "
+                << "Bus(ID=" << bus_->busID()
+                << "), Component(ID=" << component_->getIDcomponent()
+                << "). Expected " << bus_size
+                << " bus connections, but found "
+                << bus_input_ports_.size() << "."
+                << std::endl;
+
+      return 1;
+    }
+
+    // The wrapped component is evaluated independently by the interface.
+    // Allocate storage for its internal variables and residuals and redirect
+    // its internal pointers to this private storage.
+    const size_t internal_size = static_cast<size_t>(component_->getInternalSize());
+
+    y_ptr_  = std::make_unique<ScalarT[]>(internal_size);
+    yp_ptr_ = std::make_unique<ScalarT[]>(internal_size);
+    f_ptr_  = std::make_unique<ScalarT[]>(internal_size);
+
+    component_->setInternalPointer(y_ptr_.get());
+    component_->setInternalDerivativePointer(yp_ptr_.get());
+    component_->setInternalResidualPointer(f_ptr_.get());
+
+    return 0;
+  }
+
+  template <class ScalarT, typename IdxT>
+  int BusPartitionInterface<ScalarT, IdxT>::setAbsoluteTolerance(RealT rel_tol)
+  {
+    abs_tol_.setToConst(static_cast<ScalarT>(rel_tol));
+    return 0;
+  }
+
+  /**
+   * @brief Initialize the partition interface.
+   *
+   * @return 0 on success.
+   */
+  template <class ScalarT, typename IdxT>
+  int BusPartitionInterface<ScalarT, IdxT>::initialize()
+  {
+    return 0;
+  }
+
+  /**
+   * @brief Identify differential variables
+   */
+  template <class ScalarT, typename IdxT>
+  int BusPartitionInterface<ScalarT, IdxT>::tagDifferentiable()
+  {
+    return 0;
+  }
+
+  /**
+   * @brief Eval Internal Residual
+   */
+  template <class ScalarT, typename IdxT>
+  int BusPartitionInterface<ScalarT, IdxT>::evaluateInternalResidual()
+  {
+    return 0;
+  }
+
+  /**
+   * @brief Evaluate the wrapped component's contributions to the bus residual.
+   *
+   * The wrapped component is evaluated using state information supplied through
+   * the interface. Residual contributions associated with bus variables are
+   * extracted from the component's external residual and accumulated into the
+   * corresponding interface residual entries.
+   *
+   * @return 0 on success, or the error code returned by the wrapped component.
+   */
+  template <class ScalarT, typename IdxT>
+  int BusPartitionInterface<ScalarT, IdxT>::evaluateExternalResidual()
+  {
+
+    // Evaluate all external residual contributions produced by the wrapped
+    // component.
+    std::vector<ScalarT> component_ext_residual(static_cast<size_t>(component_->getExternSize()));
+
+    updateComponentPointers(component_ext_residual.data());
+
+    if (int err_code = component_->evaluateExternalResidual())
+    {
+      return err_code;
+    }
+
+    // Only contributions associated with the bus are accumulated
+    // into the interface residual below.
+    for (size_t i = 0; i < bus_input_ports_.size(); ++i)
+    {
+      *f_ext_[bus_input_ports_[i]] += component_ext_residual[bus_output_ports_[i]];
+    }
+
+    return 0;
+  }
+
+  /**
+   * @brief Evaluate the Jacobian contributions associated with the bus.
+   *
+   * Evaluates the Jacobian of the wrapped component and extracts the entries
+   * whose residual rows correspond to variables belonging to the connected bus.
+   * The selected entries are then used to assemble the interface Jacobian.
+   *
+   * @return 0 on success.
+   */
+  template <class ScalarT, typename IdxT>
+  int BusPartitionInterface<ScalarT, IdxT>::evaluateJacobian()
+  {
+
+    this->zeroJacMatrix();
+
+    // The Jacobian only requires the component state pointers.
+    //  Residual output is not needed, so external residual pointers
+    // are redirected to dummy storage.
+    updateComponentPointers(nullptr);
+
+    component_->evaluateJacobian();
+
+    const IdxT*  cooRows = component_->jacobianCooRows();
+    const IdxT*  cooCols = component_->jacobianCooCols();
+    const RealT* cooVals = component_->jacobianCooValues();
+
+    std::vector<IdxT>  rows;
+    std::vector<IdxT>  cols;
+    std::vector<RealT> vals;
+
+    rows.reserve(jac_map_.size());
+    cols.reserve(jac_map_.size());
+    vals.reserve(jac_map_.size());
+
+    // Extract only the Jacobian entries whose residual rows belong to the bus.
+    for (const IdxT index : jac_map_)
+    {
+      rows.push_back(cooRows[index]);
+      cols.push_back(cooCols[index]);
+      vals.push_back(cooVals[index]);
+    }
+
+    this->setJacValues(rows, cols, vals);
+
+    return 0;
+  }
+
+  /**
+   * @brief Update the wrapped component with state and output residual pointers.
+   *
+   * Reconstructs the state required to evaluate the wrapped component using
+   * data supplied through the interface. External component variables are
+   * connected directly to the interface data, while internal component
+   * variables are copied into the interface's private storage, which the wrapped
+   * component already keeps track of for its internal state.
+   *
+   * If residual storage is provided, external residual contributions are
+   * written to that storage. Otherwise, they are redirected to dummy storage.
+   *
+   * @pre The interface must be allocated and its external state pointers must
+   *      reference valid state and derivative data.
+   *
+   * @post The wrapped component is configured with the state, derivative, and
+   *       residual pointers required for evaluation.
+   *
+   * @param residual Storage for the component's external residual contributions,
+   *                 or nullptr when residual output is not required.
+   *
+   * @return 0 on success.
+   */
+  template <class ScalarT, typename IdxT>
+  int BusPartitionInterface<ScalarT, IdxT>::updateComponentPointers(ScalarT* residual)
+  {
+    size_t internal_index = 0;
+    size_t external_index = 0;
+
+    const auto& external_indices = component_->getExternIndices();
+
+    // Reconstruct the state expected by the wrapped component from the
+    // interface's external data. Route the residual to dummy variable if
+    // output data is not needed such as when evaluating Jacobians.
+    for (size_t i = 0; i < static_cast<size_t>(component_->size()); ++i)
+    {
+      if (is_external_[i])
+      {
+        ExternalConnection<ScalarT, IdxT> connection{
+            .y_   = y_ext_[i],
+            .yp_  = yp_ext_[i],
+            .f_   = residual ? &residual[external_index] : &dummy_residual_,
+            .idx_ = component_->getNodeConnection(i)};
+
+        component_->setExternalConnectionNodes(i, connection);
+
+        ++external_index;
+      }
+      else
+      {
+        y_ptr_[internal_index]  = *y_ext_[i];
+        yp_ptr_[internal_index] = *yp_ext_[i];
+
+        ++internal_index;
+      }
+    }
+
+    return 0;
+  }
+
+  template <class ScalarT, typename IdxT>
+  int BusPartitionInterface<ScalarT, IdxT>::evaluateIntegrand()
+  {
+    return 0;
+  }
+
+  template <class ScalarT, typename IdxT>
+  int BusPartitionInterface<ScalarT, IdxT>::initializeAdjoint()
+  {
+    return 0;
+  }
+
+  template <class ScalarT, typename IdxT>
+  int BusPartitionInterface<ScalarT, IdxT>::evaluateAdjointResidual()
+  {
+    return 0;
+  }
+
+  template <class ScalarT, typename IdxT>
+  int BusPartitionInterface<ScalarT, IdxT>::evaluateAdjointIntegrand()
+  {
+    return 0;
+  }
+
+  // Available template instantiations
+  template class BusPartitionInterface<double, long int>;
+  template class BusPartitionInterface<double, size_t>;
+  template class BusPartitionInterface<DependencyTracking::Variable, long int>;
+  template class BusPartitionInterface<DependencyTracking::Variable, size_t>;
+
+} // namespace GridKit
diff --git a/GridKit/Model/PowerElectronics/PartitionInterface/BusPartitionInterface.hpp b/GridKit/Model/PowerElectronics/PartitionInterface/BusPartitionInterface.hpp
new file mode 100644
index 000000000..0a61407e4
--- /dev/null
+++ b/GridKit/Model/PowerElectronics/PartitionInterface/BusPartitionInterface.hpp
@@ -0,0 +1,90 @@
+
+
+#pragma once
+
+#include <GridKit/Model/PowerElectronics/CircuitComponent.hpp>
+#include <GridKit/Model/PowerElectronics/NodeBase.hpp>
+#include <GridKit/Model/PowerElectronics/PartitionInterface/PartitionInterface.hpp>
+
+namespace GridKit
+{
+  /*!
+   * @brief Declaration of a passive BusPartitionInterface class.
+   *
+   */
+  template <class ScalarT, typename IdxT>
+  class BusPartitionInterface : public PartitionInterface<ScalarT, IdxT>
+  {
+
+    using component_type = CircuitComponent<ScalarT, IdxT>;
+    using node_type      = typename PowerElectronics::NodeBase<ScalarT, IdxT>;
+    using RealT          = typename CircuitComponent<ScalarT, IdxT>::RealT;
+
+    using CircuitComponent<ScalarT, IdxT>::size_;
+    using CircuitComponent<ScalarT, IdxT>::nnz_;
+    using CircuitComponent<ScalarT, IdxT>::time_;
+    using CircuitComponent<ScalarT, IdxT>::alpha_;
+    using CircuitComponent<ScalarT, IdxT>::y_ext_;
+    using CircuitComponent<ScalarT, IdxT>::y_int_;
+    using CircuitComponent<ScalarT, IdxT>::yp_ext_;
+    using CircuitComponent<ScalarT, IdxT>::yp_int_;
+    using CircuitComponent<ScalarT, IdxT>::tag_;
+    using CircuitComponent<ScalarT, IdxT>::abs_tol_;
+    using CircuitComponent<ScalarT, IdxT>::f_ext_;
+    using CircuitComponent<ScalarT, IdxT>::f_int_;
+    using CircuitComponent<ScalarT, IdxT>::g_;
+    using CircuitComponent<ScalarT, IdxT>::yB_;
+    using CircuitComponent<ScalarT, IdxT>::ypB_;
+    using CircuitComponent<ScalarT, IdxT>::fB_;
+    using CircuitComponent<ScalarT, IdxT>::gB_;
+    using CircuitComponent<ScalarT, IdxT>::param_;
+    using CircuitComponent<ScalarT, IdxT>::idc_;
+
+    using CircuitComponent<ScalarT, IdxT>::extern_indices_;
+    using CircuitComponent<ScalarT, IdxT>::n_extern_;
+    using CircuitComponent<ScalarT, IdxT>::n_intern_;
+    using CircuitComponent<ScalarT, IdxT>::connection_nodes_;
+
+  public:
+    BusPartitionInterface(node_type* bus, component_type* component, IdxT id);
+    virtual ~BusPartitionInterface();
+
+    int allocate();
+    int initialize();
+    int tagDifferentiable();
+    int setAbsoluteTolerance(RealT);
+    int evaluateInternalResidual() final;
+    int evaluateExternalResidual() final;
+    int evaluateJacobian();
+    int evaluateIntegrand();
+
+    int initializeAdjoint();
+    int evaluateAdjointResidual();
+    // int evaluateAdjointJacobian();
+    int evaluateAdjointIntegrand();
+
+  private:
+    int updateComponentPointers(ScalarT* residual);
+
+    // Component and bus associated with the partition interface.
+    component_type* component_;
+    node_type*      bus_;
+
+    // Storage for the wrapped component's internal variables.
+    std::unique_ptr<ScalarT[]> y_ptr_;
+    std::unique_ptr<ScalarT[]> yp_ptr_;
+    std::unique_ptr<ScalarT[]> f_ptr_;
+
+    std::vector<size_t> bus_input_ports_;
+    std::vector<size_t> bus_output_ports_;
+
+    // Dummy storage used when the component residual is not needed.
+    ScalarT dummy_residual_{};
+
+    // Maps interface Jacobian entries to the wrapped component Jacobian.
+    std::vector<IdxT> jac_map_;
+
+    // Identifies which wrapped component variables are external.
+    std::vector<bool> is_external_;
+  };
+} // namespace GridKit
diff --git a/GridKit/Model/PowerElectronics/PartitionInterface/CMakeLists.txt b/GridKit/Model/PowerElectronics/PartitionInterface/CMakeLists.txt
new file mode 100644
index 000000000..58b4fa884
--- /dev/null
+++ b/GridKit/Model/PowerElectronics/PartitionInterface/CMakeLists.txt
@@ -0,0 +1,5 @@
+gridkit_add_library(
+  power_elec_partition_interfaces
+  SOURCES BusPartitionInterface.cpp
+  HEADERS BusPartitionInterface.hpp PartitionInterface.hpp
+  LINK_LIBRARIES GridKit::dense_vector)
diff --git a/GridKit/Model/PowerElectronics/PartitionInterface/PartitionInterface.hpp b/GridKit/Model/PowerElectronics/PartitionInterface/PartitionInterface.hpp
new file mode 100644
index 000000000..9a2f0e072
--- /dev/null
+++ b/GridKit/Model/PowerElectronics/PartitionInterface/PartitionInterface.hpp
@@ -0,0 +1,24 @@
+
+
+#pragma once
+
+#include <GridKit/AutomaticDifferentiation/DependencyTracking/Variable.hpp>
+#include <GridKit/Model/Evaluator.hpp>
+#include <GridKit/Model/PowerElectronics/CircuitComponent.hpp>
+
+namespace GridKit
+{
+  /*!
+   * @brief Base class for partition interface components.
+   *
+   */
+  template <class ScalarT, typename IdxT>
+  class PartitionInterface : public CircuitComponent<ScalarT, IdxT>
+  {
+  public:
+    PartitionInterface() = default;
+
+    ~PartitionInterface() = default;
+  };
+
+} // namespace GridKit
diff --git a/GridKit/Model/PowerElectronics/Resistor/Resistor.cpp b/GridKit/Model/PowerElectronics/Resistor/Resistor.cpp
index 5c6c4a1cb..45efeab33 100644
--- a/GridKit/Model/PowerElectronics/Resistor/Resistor.cpp
+++ b/GridKit/Model/PowerElectronics/Resistor/Resistor.cpp
@@ -141,6 +141,18 @@ namespace GridKit
     return 0;
   }
 
+  template <class ScalarT, typename IdxT>
+  bool Resistor<ScalarT, IdxT>::isCloneable() const
+  {
+    return true;
+  }
+
+  template <class ScalarT, typename IdxT>
+  CircuitComponent<ScalarT, IdxT>* Resistor<ScalarT, IdxT>::clone() const
+  {
+    return new Resistor<ScalarT, IdxT>(*this);
+  }
+
   // Available template instantiations
   template class Resistor<double, long int>;
   template class Resistor<double, size_t>;
diff --git a/GridKit/Model/PowerElectronics/Resistor/Resistor.hpp b/GridKit/Model/PowerElectronics/Resistor/Resistor.hpp
index 540ab4d20..19c90a190 100644
--- a/GridKit/Model/PowerElectronics/Resistor/Resistor.hpp
+++ b/GridKit/Model/PowerElectronics/Resistor/Resistor.hpp
@@ -60,10 +60,13 @@ namespace GridKit
     int evaluateJacobian();
     int evaluateIntegrand();
 
-    int initializeAdjoint();
-    int evaluateAdjointResidual();
+    int  initializeAdjoint();
+    int  evaluateAdjointResidual();
     // int evaluateAdjointJacobian();
-    int evaluateAdjointIntegrand();
+    int  evaluateAdjointIntegrand();
+    bool isCloneable() const;
+
+    CircuitComponent<ScalarT, IdxT>* clone() const;
 
   private:
     RealT  R_;
diff --git a/GridKit/Model/PowerElectronics/SubsystemModel.hpp b/GridKit/Model/PowerElectronics/SubsystemModel.hpp
new file mode 100644
index 000000000..1d222b94e
--- /dev/null
+++ b/GridKit/Model/PowerElectronics/SubsystemModel.hpp
@@ -0,0 +1,833 @@
+
+
+#pragma once
+
+#include <cassert>
+#include <functional>
+#include <memory>
+#include <optional>
+#include <stdexcept>
+#include <tuple>
+#include <unordered_map>
+#include <vector>
+
+#include <GridKit/Constants.hpp>
+#include <GridKit/Model/PowerElectronics/Bus/MicrogridBus.hpp>
+#include <GridKit/Model/PowerElectronics/CircuitComponent.hpp>
+#include <GridKit/Model/PowerElectronics/CircuitNode.hpp>
+#include <GridKit/Model/PowerElectronics/PartitionInterface/PartitionInterface.hpp>
+#include <GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp>
+#include <GridKit/ScalarTraits.hpp>
+
+namespace GridKit
+{
+
+  /**
+   * @brief Represents a subset of a PowerElectronicsModel that can be evaluated
+   *        independently.
+   *
+   * A SubsystemModel contains a collection of existing GridKit components and
+   * nodes taken from a larger system. Variables owned by those components and
+   * nodes become internal variables of the subsystem. Variables needed by those
+   * components but owned outside the subsystem become external coupling
+   * variables.
+   *
+   * Components normally store connection indices in the global system indexing.
+   * During subsystem allocation, these indices are temporarily replaced with a
+   * contiguous local subsystem indexing so that the subsystem can be evaluated
+   * like an independent PowerElectronicsModel.
+   *
+   * External coupling values must be supplied before residual or Jacobian
+   * evaluation, either directly through the external-data vectors or through a
+   * forcing function.
+   *
+   * @tparam ScalarT Scalar type used by the model.
+   * @tparam IdxT Index type used for variable and connection indices.
+   */
+  template <class ScalarT, typename IdxT>
+  class SubsystemModel : public PowerElectronicsModel<ScalarT, IdxT>
+  {
+  public:
+    using ForcingData  = std::tuple<std::vector<ScalarT>, std::vector<ScalarT>>;
+    using TimeFunction = std::function<ForcingData(ScalarT)>;
+
+  protected:
+    using SystemModel    = PowerElectronicsModel<ScalarT, IdxT>;
+    using RealT          = typename CircuitComponent<ScalarT, IdxT>::RealT;
+    using CsrMatrixT     = typename CircuitComponent<ScalarT, IdxT>::CsrMatrixT;
+    using component_type = CircuitComponent<ScalarT, IdxT>;
+    using node_type      = PowerElectronics::NodeBase<ScalarT, IdxT>;
+    using interface_type = PartitionInterface<ScalarT, IdxT>;
+
+    using SystemModel::abs_tol_;
+    using SystemModel::allocated_;
+    using SystemModel::allocateVectors;
+    using SystemModel::alpha_;
+    using SystemModel::connection_nodes_;
+    using SystemModel::f_int_;
+    using SystemModel::n_extern_;
+    using SystemModel::n_intern_;
+    using SystemModel::nnz_;
+    using SystemModel::size_;
+    using SystemModel::tag_;
+    using SystemModel::time_;
+    using SystemModel::y_int_;
+    using SystemModel::yp_int_;
+
+    using SystemModel::components_;
+    using SystemModel::csr_jac_;
+    using SystemModel::jac_call_count_;
+    using SystemModel::map_to_csr_;
+    using SystemModel::neg1_;
+    using SystemModel::nodes_;
+    using SystemModel::use_jac_;
+
+  public:
+    /**
+     * @brief Default constructor for the system model
+     *
+     * @post System model parameters set as default
+     */
+
+    SubsystemModel()
+      : SystemModel(false)
+    {
+    }
+
+    /**
+     * @brief Default constructor for the system model
+     *
+     * @post System model parameters set as default
+     */
+
+    SubsystemModel(bool use_jac)
+      : SystemModel(use_jac)
+    {
+    }
+
+    ~SubsystemModel() override
+    {
+      for (auto* comp : interfaces_)
+      {
+        delete comp;
+      }
+      // SubsystemModel does not own components_/nodes_ — they belong to (and are
+      // deleted by) the parent system model.
+      components_.clear();
+      nodes_.clear();
+    }
+
+    /**
+     * @brief Allocate the subsystem for independent evaluation.
+     *
+     * Converts the selected components and nodes from global system indexing to
+     * local subsystem indexing, allocates storage for the subsystem's internal
+     * variables, and creates storage for coupling variables that lie outside the
+     * subsystem.
+     *
+     * Component pointers are then redirected to either the subsystem's internal
+     * vectors or its external coupling-data vectors. Finally, the subsystem
+     * Jacobian structure is assembled using only entries whose row and column
+     * belong to internal subsystem variables.
+     *
+     * @pre Components and nodes added to the subsystem must already be allocated
+     *      by the parent system.
+     *
+     * @post Component and node connection indices use local subsystem indexing,
+     *       and the subsystem is ready for residual and Jacobian evaluation.
+     *
+     * @return 0 on success, otherwise an error code returned during allocation.
+     */
+    int allocate() override
+    {
+      if (int err = hold())
+      {
+        return err;
+      }
+
+      n_intern_ = internal_map_.size();
+      n_extern_ = external_map_.size();
+      size_     = n_intern_;
+
+      // Allocate subsystem vectors.
+      y_ext_data_.resize(n_extern_);
+      yp_ext_data_.resize(n_extern_);
+      f_ext_data_.resize(n_extern_);
+
+      connection_nodes_ = std::make_unique<IdxT[]>(size_);
+      if (!allocated_)
+      {
+        allocateVectors(static_cast<IdxT>(size_), true);
+        abs_tol_.setToZero(memory::HOST);
+      }
+
+      tag_.resize(size_);
+
+      external_data_indices_.resize(n_extern_);
+
+      // Store the mapping from local subsystem indices back to their global system indices
+      for (const auto [global_idx, local_idx] : internal_map_)
+      {
+        this->setConnectionNodes(local_idx, global_idx);
+      }
+
+      // Store the global indices of all external coupling variables.
+      for (const auto [global_idx, local_idx] : external_map_)
+      {
+        external_data_indices_[local_idx - n_intern_] = global_idx;
+      }
+
+      { // Start node internal indexing after all component internals for proper KLU ordering
+        size_t node_internal_idx;
+        for (node_type* node : nodes_)
+        {
+          for (size_t i = 0; i < node->getInternalSize(); i++)
+          {
+            node_internal_idx = node->getNodeConnection(i).idx_;
+
+            ExternalConnection<ScalarT, IdxT> node_connection{
+                .y_   = y_int_ + node_internal_idx,
+                .yp_  = yp_int_ + node_internal_idx,
+                .f_   = f_int_ + node_internal_idx,
+                .idx_ = static_cast<IdxT>(node_internal_idx)};
+
+            node->setExternalConnectionNodes(i, node_connection);
+          }
+        }
+      }
+
+      {
+        // The offset for each component's internal variables in the system vector.
+        // They start at 0, and are stacked on top of each other.
+        size_t component_internal_idx = 0;
+        for (component_type* comp : components_)
+        {
+          // Update component internal pointers to their correct offsets
+          comp->setInternalPointer(&y_int_[component_internal_idx]);
+          comp->setInternalDerivativePointer(&yp_int_[component_internal_idx]);
+          comp->setInternalResidualPointer(&f_int_[component_internal_idx]);
+
+          component_internal_idx += comp->getInternalSize();
+
+          const auto& external_indices = comp->getExternIndices();
+
+          for (IdxT local_index : external_indices)
+          {
+            const IdxT connection_index = comp->getNodeConnection(local_index);
+
+            // A variable can be external from the component's point of view while still
+            // being owned by this subsystem. In that case, connect the component directly
+            // to the corresponding entry in the subsystem's internal vectors.
+            if (connection_index < n_intern_)
+            {
+              ExternalConnection<ScalarT, IdxT> connection{
+                  .y_   = y_int_ + connection_index,
+                  .yp_  = yp_int_ + connection_index,
+                  .f_   = f_int_ + connection_index,
+                  .idx_ = connection_index};
+
+              comp->setExternalConnectionNodes(local_index, connection);
+
+              continue;
+            }
+
+            // Otherwise the variable is owned outside this subsystem. Connect the
+            // component to the subsystem's external coupling-data storage instead.
+            const IdxT external_offset = connection_index - static_cast<IdxT>(n_intern_);
+
+            ExternalConnection<ScalarT, IdxT> connection{
+                .y_   = &y_ext_data_[external_offset],
+                .yp_  = &yp_ext_data_[external_offset],
+                .f_   = &f_ext_data_[external_offset],
+                .idx_ = connection_index};
+
+            comp->setExternalConnectionNodes(local_index, connection);
+          }
+        }
+      }
+
+      // Allocation always rebuilds the system Jacobian and its COO-to-CSR map.
+      delete csr_jac_;
+      csr_jac_ = nullptr;
+
+      delete[] map_to_csr_;
+      map_to_csr_ = nullptr;
+
+      // Evaluate component Jacobians to get sparsity
+      for (component_type* component : components_)
+      {
+        component->evaluateJacobian();
+      }
+
+      // Check whether a Jacobian entry belongs to the subsystem Jacobian.
+      // Only entries whose row and column are both internal subsystem
+      // variables are retained.
+      auto isValidEntry = [this](IdxT row, IdxT col)
+      {
+        if (row == neg1_ || col == neg1_)
+        {
+          return false;
+        }
+
+        const bool row_is_internal = row < this->getInternalSize();
+        const bool col_is_internal = col < this->getInternalSize();
+
+        return (row_is_internal && col_is_internal);
+      };
+
+      IdxT nnz_dup = 0;
+
+      for (const component_type* component : components_)
+      {
+        const IdxT* r   = component->jacobianCooRows();
+        const IdxT* c   = component->jacobianCooCols();
+        const IdxT  nnz = component->nnz();
+
+        for (IdxT i = 0; i < nnz; ++i)
+        {
+          const IdxT row = component->getNodeConnection(r[i]);
+          const IdxT col = component->getNodeConnection(c[i]);
+
+          if (isValidEntry(row, col))
+          {
+            ++nnz_dup;
+          }
+        }
+      }
+
+      // Allocate COO triplet arrays (we own these until we hand off to CsrMatrix)
+      IdxT*  rows_dup = new IdxT[nnz_dup];
+      IdxT*  cols_dup = new IdxT[nnz_dup];
+      RealT* vals_dup = new RealT[nnz_dup];
+
+      IdxT counter = 0;
+
+      for (const component_type* component : components_)
+      {
+        const IdxT*  r   = component->jacobianCooRows();
+        const IdxT*  c   = component->jacobianCooCols();
+        const RealT* v   = component->jacobianCooValues();
+        const IdxT   nnz = component->nnz();
+
+        for (IdxT i = 0; i < nnz; ++i)
+        {
+          const IdxT row = component->getNodeConnection(r[i]);
+          const IdxT col = component->getNodeConnection(c[i]);
+
+          if (!isValidEntry(row, col))
+          {
+            continue;
+          }
+
+          rows_dup[counter] = row;
+          cols_dup[counter] = col;
+          vals_dup[counter] = v[i];
+
+          ++counter;
+        }
+      }
+
+      // Build the system COO Jacobian
+      LinearAlgebra::CooMatrix<RealT, IdxT> jac(size_, size_, nnz_dup, &rows_dup, &cols_dup, &vals_dup);
+
+      // Populate CSR data with sort and deduplicate
+      IdxT* row_ptrs = jac.getCsrRowData();
+
+      // Deduplicated nnz
+      nnz_ = jac.getNnz();
+
+      // Allocate cols/vals with deduplicated nnz
+      IdxT*  cols = new IdxT[nnz_];
+      RealT* vals = new RealT[nnz_];
+
+      std::copy(jac.getColData(), jac.getColData() + nnz_, cols);
+      std::copy(jac.getValues(), jac.getValues() + nnz_, vals);
+
+      // Create the CSR Jacobian
+      csr_jac_ = new CsrMatrixT(size_, size_, nnz_, &row_ptrs, &cols, &vals);
+
+      const IdxT* map_to_sorted = jac.getMapToSorted();
+      const IdxT* map_to_dedup  = jac.getMapToDeduplicated();
+
+      // Build a mappping from original COO index to CSR index
+      map_to_csr_ = new IdxT[nnz_dup];
+      for (IdxT i = 0; i < nnz_dup; ++i)
+      {
+        map_to_csr_[map_to_sorted[i]] = map_to_dedup[i];
+      }
+
+      allocated_ = true;
+      return 0;
+    }
+
+    /**
+     * @brief Update the subsystem external state and derivative data.
+     *
+     * If a forcing function is provided, evaluate it at the current subsystem
+     * time and copy the returned coupling values into the external state vectors.
+     *
+     * @post y_ext_data_ and yp_ext_data_ contain the external coupling values
+     *       returned by the forcing function, if one is set.
+     *
+     * @throws std::runtime_error If the forcing function returns vectors whose
+     *         sizes do not match the subsystem external-data vectors.
+     *
+     * @return 0 on success.
+     */
+    int distributeExternalVectors()
+    {
+
+      if (forcing_function_)
+      {
+        const auto [y_forcing, yp_forcing] = (*forcing_function_)(time_);
+
+        if (y_forcing.size() != y_ext_data_.size() || yp_forcing.size() != yp_ext_data_.size())
+        {
+          throw std::runtime_error(
+              "SubsystemModel::distributeExternalVectors: forcing function "
+              "returned vectors with incorrect sizes.");
+        }
+
+        std::copy(y_forcing.begin(), y_forcing.end(), y_ext_data_.begin());
+        std::copy(yp_forcing.begin(), yp_forcing.end(), yp_ext_data_.begin());
+      }
+
+      return 0;
+    }
+
+    /**
+     * @brief Evaluate Residuals at each component then collect them
+     *
+     * @return int 0 if successful, positive if there's a recoverable error, negative if unrecoverable
+     */
+    int evaluateInternalResidual() override
+    {
+      if (int err_code = distributeExternalVectors())
+      {
+        return err_code;
+      }
+
+      return SystemModel::evaluateInternalResidual();
+    }
+
+    /**
+     * @brief Creates the system Jacobian representing \f$\alpha dF/dy' + dF/dy\f$
+     *
+     * Updates the CSR Jacobian values using the per-component mappings
+     * computed during allocate().
+     *
+     * @return int 0 if successful, positive if there's a recoverable error, negative if unrecoverable
+     */
+    int evaluateJacobian() override
+    {
+      if (int err_code = distributeExternalVectors())
+      {
+        return err_code;
+      }
+
+      return SystemModel::evaluateJacobian();
+    }
+
+    /**
+     * @brief Add a component to the subsystem.
+     *
+     * Rejected while connections are in the local-indexed state, since the
+     * component's stored connection indices would otherwise be interpreted
+     * inconsistently with the rest of the subsystem. Call release() first.
+     *
+     * @param[in] component Component to add.
+     */
+    void addComponent(component_type* component)
+    {
+      if (!component->isAllocated())
+      {
+        throw std::logic_error(
+            "SubsystemModel::addComponent: cannot add an unallocated component.");
+      }
+
+      if (connections_are_local_)
+      {
+        throw std::logic_error(
+            "SubsystemModel::addComponent: cannot add component while "
+            "in local-indexed state. Call release() first.");
+      }
+
+      SystemModel::addComponent(component);
+    }
+
+    /**
+     * @brief Add a node to the subsystem.
+     *
+     * Rejected while connections are in the local-indexed state, since the
+     * node's stored connection indices would otherwise be interpreted
+     * inconsistently with the rest of the subsystem. Call release() first.
+     *
+     * @param[in] node Node to add.
+     */
+    void addNode(node_type* node)
+    {
+      if (!node->isAllocated())
+      {
+        throw std::logic_error(
+            "SubsystemModel::addNode: cannot add an unallocated node.");
+      }
+
+      if (connections_are_local_)
+      {
+        throw std::logic_error(
+            "SubsystemModel::addNode: cannot add node while in "
+            "local-indexed state. Call release() first.");
+      }
+
+      SystemModel::addNode(node);
+    }
+
+    /**
+     * @brief Add a partition interface to the subsystem.
+     *
+     * Adds the partition interface to the subsystem's component list and keeps a
+     * separate reference to it in the interface list.
+     *
+     * @param component Pointer to the interface component to add.
+     */
+    void addInterface(interface_type* component)
+    {
+      addComponent(component);
+      interfaces_.push_back(component);
+    }
+
+    /**
+     * @brief Prepare the subsystem topology for local evaluation.
+     *
+     * Components and nodes initially contain the same global connection indices
+     * assigned by the parent system. This method determines which of those
+     * variables belong to the subsystem and which are external coupling variables,
+     * assigns each a local subsystem index, and replaces the stored global
+     * connection indices with those local indices.
+     *
+     * This local indexing is used while the subsystem is allocated and evaluated.
+     * release() restores the original global indices.
+     *
+     * @pre Component and node connections use global system indices.
+     *
+     * @post Component and node connections use local subsystem indices.
+     *
+     * @return int 0 if successful, positive if there's a recoverable error, negative if unrecoverable
+     */
+    int hold()
+    {
+      buildIndexMappings();
+
+      if (int err_code = mapGlobalToLocal())
+      {
+        return err_code;
+      }
+
+      return 0;
+    }
+
+    /**
+     * @brief Restore the subsystem topology to global system indexing.
+     *
+     * Reverses hold() by replacing local subsystem connection indices with the
+     * original global system indices. The local internal/external mappings are
+     * then cleared so that components or nodes can safely be added or removed.
+     *
+     * @pre Component and node connections may use local subsystem indices.
+     *
+     * @post Component and node connections use their original global indices,
+     *       the subsystem mappings are cleared, and the subsystem is marked
+     *       unallocated.
+     *
+     * @return int 0 if successful, positive if there's a recoverable error, negative if unrecoverable
+     */
+    int release()
+    {
+      if (int err_code = mapLocalToGlobal())
+      {
+        return err_code;
+      }
+
+      internal_map_.clear();
+      external_map_.clear();
+
+      allocated_ = false;
+
+      return 0;
+    }
+
+    const std::vector<IdxT>& getExternalDataIndices() const
+    {
+      return external_data_indices_;
+    }
+
+    std::vector<ScalarT>& getExternalDataY()
+    {
+      return y_ext_data_;
+    }
+
+    std::vector<ScalarT>& getExternalDataYP()
+    {
+      return yp_ext_data_;
+    }
+
+    std::vector<ScalarT>& getExternalDataF()
+    {
+      return f_ext_data_;
+    }
+
+    void setForcingFunction(TimeFunction function)
+    {
+      forcing_function_ = std::move(function);
+    }
+
+    const std::unordered_map<IdxT, IdxT>& getInternalMap() const
+    {
+      return internal_map_;
+    }
+
+    const std::unordered_map<IdxT, IdxT>& getExternalMap() const
+    {
+      return external_map_;
+    }
+
+  private:
+    /**
+     * @brief Restore local subsystem connection indices to global system indices.
+     *
+     * Reverses mapGlobalToLocal(). Internal subsystem indices are translated
+     * through the subsystem connection-node table, while external subsystem
+     * indices are translated through external_data_indices_.
+     *
+     * The component/node connectivity is unchanged; only the index representation
+     * is restored.
+     *
+     * @pre Component and node connections use local subsystem indices.
+     *
+     * @post Component and node connections use their original global indices.
+     *
+     * @return 0 on success.
+     */
+    int mapLocalToGlobal()
+    {
+      if (!connections_are_local_)
+      {
+        return 0;
+      }
+
+      for (component_type* component : components_)
+      {
+
+        for (IdxT i = 0; i < component->size(); i++)
+        {
+          const IdxT index = component->getNodeConnection(i);
+
+          if (index == neg1_)
+          {
+            continue;
+          }
+          else if (index < this->getInternalSize())
+          {
+            component->setConnectionNodes(i, this->getNodeConnection(index));
+          }
+          else
+          {
+            component->setConnectionNodes(i, external_data_indices_[index - this->getInternalSize()]);
+          }
+        }
+      }
+
+      for (node_type* node : nodes_)
+      {
+
+        for (IdxT i = 0; i < node->size(); i++)
+        {
+          const IdxT index = node->getNodeConnection(i).idx_;
+
+          if (index == neg1_)
+          {
+            continue;
+          }
+          node->setConnectionNodes(i, this->getNodeConnection(index));
+        }
+      }
+
+      connections_are_local_ = false;
+
+      return 0;
+    }
+
+    /**
+     * @brief Replace global connection indices with subsystem-local indices.
+     *
+     * Uses the mappings created by buildIndexMappings() to rewrite the connection
+     * indices stored by every component and node. Variables owned by the subsystem
+     * use internal_map_; variables owned outside the subsystem use external_map_.
+     *
+     * This changes only the indexing used to identify connections; it does not
+     * change the physical component/node connectivity.
+     *
+     * @pre internal_map_ and external_map_ have been constructed from global
+     *      connection indices.
+     *
+     * @post All valid component and node connections use local subsystem indices.
+     *
+     * @return 0 on success.
+     */
+    int mapGlobalToLocal()
+    {
+      if (connections_are_local_)
+      {
+        return 0;
+      }
+
+      for (component_type* component : components_)
+      {
+        for (IdxT i = 0; i < component->size(); ++i)
+        {
+          const IdxT index = component->getNodeConnection(i);
+
+          if (index == neg1_)
+          {
+            continue;
+          }
+
+          if (internal_map_.contains(index))
+          {
+            component->setConnectionNodes(i, internal_map_.at(index));
+          }
+          else
+          {
+            component->setConnectionNodes(i, external_map_.at(index));
+          }
+        }
+      }
+
+      for (node_type* node : nodes_)
+      {
+        for (IdxT i = 0; i < node->size(); ++i)
+        {
+          const IdxT index = node->getNodeConnection(i).idx_;
+
+          if (index == neg1_)
+          {
+            continue;
+          }
+
+          node->setConnectionNodes(i, internal_map_.at(index));
+        }
+      }
+
+      connections_are_local_ = true;
+
+      return 0;
+    }
+
+    /**
+     * @brief Build the global-to-local variable mappings for the subsystem.
+     *
+     * Examines the connection indices of all components and nodes in the
+     * subsystem and divides the referenced variables into two groups:
+     *
+     * - Internal variables are owned by a component or node in this subsystem.
+     * - External variables are required by a subsystem component but are owned
+     *   outside the subsystem.
+     *
+     * Internal variables receive local indices first. External coupling variables
+     * are then assigned indices immediately after the internal range. This gives
+     * every variable referenced by the subsystem a unique local index while
+     * preserving its original global index in the corresponding map.
+     *
+     * @pre Component and node connections use global system indices.
+     *
+     * @post internal_map_ and external_map_ contain the global-to-local mappings
+     *       needed to convert the subsystem topology to local indexing.
+     */
+    void buildIndexMappings()
+    {
+
+      if (connections_are_local_)
+      {
+        return;
+      }
+
+      internal_map_.clear();
+      external_map_.clear();
+
+      size_t component_internal_idx = 0;
+      // Pass 1: Add variables owned internally by subsystem components.
+      for (component_type* comp : components_)
+      {
+        const auto& extern_indices = comp->getExternIndices();
+
+        for (IdxT i = 0; i < comp->size(); i++)
+        {
+          const IdxT index = comp->getNodeConnection(i);
+
+          if (index != neg1_ && !extern_indices.contains(i))
+          {
+            internal_map_[index] = component_internal_idx++;
+          }
+        }
+      }
+
+      // Pass 2: Add variables owned by subsystem nodes.
+      for (node_type* node : nodes_)
+      {
+
+        for (IdxT i = 0; i < node->size(); i++)
+        {
+          const IdxT index = node->getNodeConnection(i).idx_;
+
+          if (index != neg1_)
+          {
+            internal_map_[index] = component_internal_idx++;
+          }
+        }
+      }
+
+      // Pass 3: Add component dependencies that are owned outside the subsystem.
+      for (component_type* comp : components_)
+      {
+        auto extern_indices = comp->getExternIndices();
+        for (IdxT j = 0; j < comp->size(); j++)
+        {
+          if (!extern_indices.contains(j))
+          {
+            continue;
+          }
+          const IdxT index = comp->getNodeConnection(j);
+
+          if (internal_map_.count(index) < 1 && external_map_.count(index) < 1 && index != neg1_)
+          {
+            external_map_[index] = component_internal_idx++;
+          }
+        }
+      }
+    }
+
+    // Global system index -> local subsystem index for subsystem internal variables
+    std::unordered_map<IdxT, IdxT> internal_map_;
+
+    // Global system index -> local subsystem index for subsystem external variables
+    std::unordered_map<IdxT, IdxT> external_map_;
+
+    // Global system index corresponding to each entry in the subsystem external state vectors
+    std::vector<IdxT> external_data_indices_;
+
+    // subsystem external State, derivative, and residual vectors.
+    std::vector<ScalarT> y_ext_data_;
+    std::vector<ScalarT> yp_ext_data_;
+    std::vector<ScalarT> f_ext_data_;
+
+    std::optional<TimeFunction> forcing_function_;
+
+    // Partition interfaces are owned by SubsystemModel.
+    // components_ contains non-owning pointers, including these interfaces.
+    std::vector<component_type*> interfaces_;
+
+    bool connections_are_local_{false};
+
+  }; // class SubsystemModel
+
+} // namespace GridKit
diff --git a/GridKit/Model/PowerElectronics/SynchronousMachine/SynchronousMachine.cpp b/GridKit/Model/PowerElectronics/SynchronousMachine/SynchronousMachine.cpp
index c3a647762..b10480f33 100644
--- a/GridKit/Model/PowerElectronics/SynchronousMachine/SynchronousMachine.cpp
+++ b/GridKit/Model/PowerElectronics/SynchronousMachine/SynchronousMachine.cpp
@@ -80,6 +80,25 @@ namespace GridKit
     return 0;
   }
 
+  /**
+   * @brief Compute the absolute tolerance for each variable in the model
+   *
+   * @param rel_tol The relative tolerance which can be used to pick the
+   *        absolute tolerance.
+   * @tparam ScalarT Scalar data type
+   * @tparam IdxT Index data type
+   * @return int 0 if successful, non-zero otherwise.
+   *
+   * This represents a "noise" level close to zero for which pure relative
+   * error cannot be used.
+   */
+  template <class ScalarT, typename IdxT>
+  int SynchronousMachine<ScalarT, IdxT>::setAbsoluteTolerance(RealT rel_tol)
+  {
+    abs_tol_.setToConst(static_cast<ScalarT>(rel_tol));
+    return 0;
+  }
+
   /**
    * @brief Compute the resisdual of the component.
    *
@@ -163,6 +182,18 @@ namespace GridKit
     return 0;
   }
 
+  template <class ScalarT, typename IdxT>
+  bool SynchronousMachine<ScalarT, IdxT>::isCloneable() const
+  {
+    return true;
+  }
+
+  template <class ScalarT, typename IdxT>
+  CircuitComponent<ScalarT, IdxT>* SynchronousMachine<ScalarT, IdxT>::clone() const
+  {
+    return new SynchronousMachine<ScalarT, IdxT>(*this);
+  }
+
   // Available template instantiations
   template class SynchronousMachine<double, long int>;
   template class SynchronousMachine<double, size_t>;
diff --git a/GridKit/Model/PowerElectronics/SynchronousMachine/SynchronousMachine.hpp b/GridKit/Model/PowerElectronics/SynchronousMachine/SynchronousMachine.hpp
index 80b68eed2..7cf040583 100644
--- a/GridKit/Model/PowerElectronics/SynchronousMachine/SynchronousMachine.hpp
+++ b/GridKit/Model/PowerElectronics/SynchronousMachine/SynchronousMachine.hpp
@@ -31,6 +31,7 @@ namespace GridKit
     using CircuitComponent<ScalarT, IdxT>::y_int_;
     using CircuitComponent<ScalarT, IdxT>::yp_ext_;
     using CircuitComponent<ScalarT, IdxT>::yp_int_;
+    using CircuitComponent<ScalarT, IdxT>::abs_tol_;
     using CircuitComponent<ScalarT, IdxT>::tag_;
     using CircuitComponent<ScalarT, IdxT>::f_ext_;
     using CircuitComponent<ScalarT, IdxT>::f_int_;
@@ -52,15 +53,19 @@ namespace GridKit
 
     int initialize();
     int tagDifferentiable();
+    int setAbsoluteTolerance(RealT);
     int evaluateInternalResidual() final;
     int evaluateExternalResidual() final;
     int evaluateJacobian();
     int evaluateIntegrand();
 
-    int initializeAdjoint();
-    int evaluateAdjointResidual();
+    int  initializeAdjoint();
+    int  evaluateAdjointResidual();
     // int evaluateAdjointJacobian();
-    int evaluateAdjointIntegrand();
+    int  evaluateAdjointIntegrand();
+    bool isCloneable() const;
+
+    CircuitComponent<ScalarT, IdxT>* clone() const;
 
   private:
     RealT                    Lls_;
diff --git a/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp b/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp
index 174e28568..3f1aee396 100644
--- a/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp
+++ b/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp
@@ -23,10 +23,12 @@ namespace GridKit
     using component_type = CircuitComponent<ScalarT, IdxT>;
     using node_type      = PowerElectronics::NodeBase<ScalarT, IdxT>;
 
+  protected:
     using Base::abs_tol_;
     using Base::allocated_;
     using Base::allocateVectors;
     using Base::alpha_;
+    using Base::connection_nodes_;
     using Base::f_ext_;
     using Base::f_int_;
     using Base::n_extern_;
@@ -41,17 +43,6 @@ namespace GridKit
     using Base::yp_int_;
 
   public:
-    /**
-     * @brief Default constructor for the system model
-     *
-     * @post System model parameters set as default
-     */
-    PowerElectronicsModel()
-    {
-      // By default don't use the jacobian
-      use_jac_ = false;
-    }
-
     /**
      * @brief Constructor for the system model
      *
@@ -117,7 +108,7 @@ namespace GridKit
      *
      * @return int 0 if successful, positive if there's a recoverable error, negative if unrecoverable
      */
-    int allocate() final
+    int allocate() override
     {
       size_t component_internal_size = 0;
       for (component_type* comp : components_)
@@ -293,7 +284,7 @@ namespace GridKit
       return Base::initialize();
     }
 
-    int tagDifferentiable() final
+    int tagDifferentiable() override
     {
       return 0;
     }
@@ -321,7 +312,7 @@ namespace GridKit
      *
      * @return int 0 if successful, positive if there's a recoverable error, negative if unrecoverable
      */
-    int evaluateInternalResidual() final
+    int evaluateInternalResidual() override
     {
       for (IdxT i = 0; i < size_; i++)
       {
@@ -362,7 +353,7 @@ namespace GridKit
      *
      * @return int 0 if successful, positive if there's a recoverable error, negative if unrecoverable
      */
-    int evaluateJacobian() final
+    int evaluateJacobian() override
     {
       // Zero out values
       RealT* vals = csr_jac_->getValues();
@@ -384,11 +375,18 @@ namespace GridKit
 
         for (IdxT i = 0; i < nnz; ++i)
         {
-          if (component->getNodeConnection(r[i]) != neg1_ && component->getNodeConnection(c[i]) != neg1_)
+          const IdxT row = component->getNodeConnection(r[i]);
+          const IdxT col = component->getNodeConnection(c[i]);
+
+          const bool is_internal_entry = row != neg1_ && col != neg1_ && row < n_intern_ && col < n_intern_;
+
+          if (!is_internal_entry)
           {
-            vals[map_to_csr_[counter]] += v[i];
-            ++counter;
+            continue;
           }
+
+          vals[map_to_csr_[counter]] += v[i];
+          ++counter;
         }
       }
 
@@ -468,7 +466,7 @@ namespace GridKit
       allocated_ = false;
     }
 
-  private:
+  protected:
     static constexpr IdxT neg1_ = INVALID_INDEX<IdxT>;
 
     std::vector<component_type*> components_;
diff --git a/GridKit/Model/PowerElectronics/TransmissionLine/TransmissionLine.cpp b/GridKit/Model/PowerElectronics/TransmissionLine/TransmissionLine.cpp
index 3adbbef58..c959c2b5c 100644
--- a/GridKit/Model/PowerElectronics/TransmissionLine/TransmissionLine.cpp
+++ b/GridKit/Model/PowerElectronics/TransmissionLine/TransmissionLine.cpp
@@ -63,6 +63,25 @@ namespace GridKit
     return 0;
   }
 
+  /**
+   * @brief Compute the absolute tolerance for each variable in the model
+   *
+   * @param rel_tol The relative tolerance which can be used to pick the
+   *        absolute tolerance.
+   * @tparam ScalarT Scalar data type
+   * @tparam IdxT Index data type
+   * @return int 0 if successful, non-zero otherwise.
+   *
+   * This represents a "noise" level close to zero for which pure relative
+   * error cannot be used.
+   */
+  template <class ScalarT, typename IdxT>
+  int TransmissionLine<ScalarT, IdxT>::setAbsoluteTolerance(RealT rel_tol)
+  {
+    abs_tol_.setToConst(static_cast<ScalarT>(rel_tol));
+    return 0;
+  }
+
   /**
    * @brief Evaluate residual of transmission line
    *
@@ -184,6 +203,18 @@ namespace GridKit
     return 0;
   }
 
+  template <class ScalarT, typename IdxT>
+  bool TransmissionLine<ScalarT, IdxT>::isCloneable() const
+  {
+    return true;
+  }
+
+  template <class ScalarT, typename IdxT>
+  CircuitComponent<ScalarT, IdxT>* TransmissionLine<ScalarT, IdxT>::clone() const
+  {
+    return new TransmissionLine<ScalarT, IdxT>(*this);
+  }
+
   // Available template instantiations
   template class TransmissionLine<double, long int>;
   template class TransmissionLine<double, size_t>;
diff --git a/GridKit/Model/PowerElectronics/TransmissionLine/TransmissionLine.hpp b/GridKit/Model/PowerElectronics/TransmissionLine/TransmissionLine.hpp
index 091984c89..c9743e9dd 100644
--- a/GridKit/Model/PowerElectronics/TransmissionLine/TransmissionLine.hpp
+++ b/GridKit/Model/PowerElectronics/TransmissionLine/TransmissionLine.hpp
@@ -33,6 +33,7 @@ namespace GridKit
     using CircuitComponent<ScalarT, IdxT>::y_int_;
     using CircuitComponent<ScalarT, IdxT>::yp_ext_;
     using CircuitComponent<ScalarT, IdxT>::yp_int_;
+    using CircuitComponent<ScalarT, IdxT>::abs_tol_;
     using CircuitComponent<ScalarT, IdxT>::tag_;
     using CircuitComponent<ScalarT, IdxT>::f_ext_;
     using CircuitComponent<ScalarT, IdxT>::f_int_;
@@ -54,15 +55,19 @@ namespace GridKit
 
     int initialize();
     int tagDifferentiable();
+    int setAbsoluteTolerance(RealT);
     int evaluateInternalResidual() final;
     int evaluateExternalResidual() final;
     int evaluateJacobian();
     int evaluateIntegrand();
 
-    int initializeAdjoint();
-    int evaluateAdjointResidual();
+    int  initializeAdjoint();
+    int  evaluateAdjointResidual();
     // int evaluateAdjointJacobian();
-    int evaluateAdjointIntegrand();
+    int  evaluateAdjointIntegrand();
+    bool isCloneable() const;
+
+    CircuitComponent<ScalarT, IdxT>* clone() const;
 
   private:
     RealT R_;
diff --git a/GridKit/Model/PowerElectronics/VoltageSource/VoltageSource.cpp b/GridKit/Model/PowerElectronics/VoltageSource/VoltageSource.cpp
index 9a4accd0b..a0d0d3946 100644
--- a/GridKit/Model/PowerElectronics/VoltageSource/VoltageSource.cpp
+++ b/GridKit/Model/PowerElectronics/VoltageSource/VoltageSource.cpp
@@ -140,6 +140,18 @@ namespace GridKit
     return 0;
   }
 
+  template <class ScalarT, typename IdxT>
+  bool VoltageSource<ScalarT, IdxT>::isCloneable() const
+  {
+    return true;
+  }
+
+  template <class ScalarT, typename IdxT>
+  CircuitComponent<ScalarT, IdxT>* VoltageSource<ScalarT, IdxT>::clone() const
+  {
+    return new VoltageSource<ScalarT, IdxT>(*this);
+  }
+
   // Available template instantiations
   template class VoltageSource<double, long int>;
   template class VoltageSource<double, size_t>;
diff --git a/GridKit/Model/PowerElectronics/VoltageSource/VoltageSource.hpp b/GridKit/Model/PowerElectronics/VoltageSource/VoltageSource.hpp
index 814fa5ec5..ab9125f40 100644
--- a/GridKit/Model/PowerElectronics/VoltageSource/VoltageSource.hpp
+++ b/GridKit/Model/PowerElectronics/VoltageSource/VoltageSource.hpp
@@ -60,10 +60,13 @@ namespace GridKit
     int evaluateJacobian() final;
     int evaluateIntegrand();
 
-    int initializeAdjoint();
-    int evaluateAdjointResidual();
+    int  initializeAdjoint();
+    int  evaluateAdjointResidual();
     // int evaluateAdjointJacobian();
-    int evaluateAdjointIntegrand();
+    int  evaluateAdjointIntegrand();
+    bool isCloneable() const;
+
+    CircuitComponent<ScalarT, IdxT>* clone() const;
 
   private:
     RealT  V_;
diff --git a/examples/PowerElectronics/CMakeLists.txt b/examples/PowerElectronics/CMakeLists.txt
index 2bf311e22..bdf85eb07 100644
--- a/examples/PowerElectronics/CMakeLists.txt
+++ b/examples/PowerElectronics/CMakeLists.txt
@@ -10,5 +10,6 @@ if(TARGET SUNDIALS::idas)
     add_subdirectory(RLCircuit)
     add_subdirectory(Microgrid)
     add_subdirectory(ScaleMicrogrid)
+    add_subdirectory(Partition)
   endif()
 endif()
diff --git a/examples/PowerElectronics/Common/JacTestHelper.hpp b/examples/PowerElectronics/Common/JacTestHelper.hpp
new file mode 100644
index 000000000..8cd7db7e2
--- /dev/null
+++ b/examples/PowerElectronics/Common/JacTestHelper.hpp
@@ -0,0 +1,282 @@
+#pragma once
+
+#include <cmath>
+#include <iostream>
+#include <optional>
+#include <unordered_map>
+#include <vector>
+
+#include <GridKit/LinearAlgebra/SparseMatrix/CsrMatrix.hpp>
+#include <GridKit/MemoryUtilities/MemoryUtils.hpp>
+#include <GridKit/Model/PowerElectronics/SubsystemModel.hpp>
+#include <GridKit/Testing/TestHelpers.hpp>
+
+#include "GridKit/Testing/Testing.hpp"
+
+namespace GridKit
+{
+  namespace Testing
+  {
+    /**
+     * @brief Verify that the subsystem Jacobian is an exact subset of the
+     * full-system Jacobian to ensure accuracy of subsystem Jacobians.
+     *
+     * Each subsystem Jacobian entry should match the corresponding entry in the
+     * full-system Jacobian after converting local subsystem indices back to global
+     * indices. The check verifies that:
+     *
+     * - every subsystem entry exists in the full-system Jacobian,
+     * - every full-system entry whose row and column belong to the subsystem
+     * exists in the subsystem Jacobian, and
+     * - matching entries agree within the specified tolerance.
+     *
+     * @note When the components in a subsystem are evaluated in a different order
+     *       than in the monolithic reference system, a larger tolerance may be
+     *       required to account for floating-point roundoff introduced by the
+     *       different summation order.
+     *
+     * @param full_jac Full-system Jacobian.
+     * @param sub_jac Subsystem Jacobian.
+     * @param subsystem Subsystem associated with sub_jac.
+     * @param tolerance Maximum allowed difference between matching entries.
+     *
+     * @return true if the subsystem Jacobian is an exact subset of the
+     * full-system Jacobian, false otherwise.
+     */
+
+    template <typename RealT, typename IdxT>
+    bool verifySubsystemJacobian(
+        GridKit::LinearAlgebra::CsrMatrix<RealT, IdxT>& full_jac,
+        GridKit::LinearAlgebra::CsrMatrix<RealT, IdxT>& sub_jac,
+        GridKit::SubsystemModel<RealT, IdxT>&           subsystem,
+        std::optional<RealT>                            tolerance = std::nullopt)
+    {
+      constexpr auto host      = memory::HOST;
+      // Full-system CSR data.
+      const IdxT*    full_rows = full_jac.getRowData(host);
+      const IdxT*    full_cols = full_jac.getColData(host);
+      const RealT*   full_vals = full_jac.getValues(host);
+
+      // Subsystem CSR data.
+      const IdxT*  sub_rows = sub_jac.getRowData(host);
+      const IdxT*  sub_cols = sub_jac.getColData(host);
+      const RealT* sub_vals = sub_jac.getValues(host);
+
+      const IdxT num_sub_rows = sub_jac.getNumRows();
+      const IdxT num_sub_cols = sub_jac.getNumColumns();
+
+      // The subsystem internal map stores global-to-local indices.
+      const auto& global_to_local = subsystem.getInternalMap();
+
+      // The internal map should contain one entry for every subsystem row.
+      if (global_to_local.size() != num_sub_rows)
+      {
+        std::cout << "Internal map size mismatch: map has "
+                  << global_to_local.size()
+                  << " entries, but subsystem Jacobian has "
+                  << num_sub_rows
+                  << " rows\n";
+
+        return false;
+      }
+
+      // The subsystem Jacobian is expected to be square.
+      if (num_sub_rows != num_sub_cols)
+      {
+        std::cout << "Subsystem Jacobian is not square: "
+                  << num_sub_rows
+                  << " rows and "
+                  << num_sub_cols
+                  << " columns\n";
+
+        return false;
+      }
+
+      // Build the reverse map from local subsystem indices to global indices.
+      // This is needed to compare subsystem rows and columns with the
+      // corresponding entries in the full-system Jacobian.
+      std::vector<IdxT> local_to_global(num_sub_rows);
+      std::vector<bool> local_index_found(num_sub_rows, false);
+
+      for (const auto& [global_index, local_index] : global_to_local)
+      {
+
+        // Make sure the global index is valid for the full-system Jacobian.
+        if (global_index >= full_jac.getNumRows())
+        {
+          std::cout << "Invalid global index "
+                    << global_index
+                    << " in subsystem internal map\n";
+
+          return false;
+        }
+
+        // Make sure the local index is valid for the subsystem Jacobian.
+        if (local_index >= num_sub_rows)
+        {
+          std::cout << "Invalid local index "
+                    << local_index
+                    << " mapped from global index "
+                    << global_index << '\n';
+
+          return false;
+        }
+
+        // Each local index should correspond to only one global index.
+        if (local_index_found[local_index])
+        {
+          std::cout << "Duplicate local index "
+                    << local_index
+                    << " in subsystem internal map\n";
+
+          return false;
+        }
+
+        local_to_global[local_index]   = global_index;
+        local_index_found[local_index] = true;
+      }
+
+      // Make sure every subsystem local index has a global index.
+      for (IdxT local_index = 0; local_index < num_sub_rows; ++local_index)
+      {
+        if (!local_index_found[local_index])
+        {
+          std::cout << "No global index maps to local index "
+                    << local_index << '\n';
+
+          return false;
+        }
+      }
+
+      bool matches = true;
+
+      // Compare each subsystem row with the corresponding full-system row.
+      for (IdxT local_row = 0; local_row < num_sub_rows; ++local_row)
+      {
+        const IdxT global_row = local_to_global[local_row];
+
+        const IdxT sub_begin  = sub_rows[local_row];
+        const IdxT sub_end    = sub_rows[local_row + 1];
+        const IdxT full_begin = full_rows[global_row];
+        const IdxT full_end   = full_rows[global_row + 1];
+
+        /*
+         * Store the current subsystem row using global column indices.
+         *
+         * The subsystem Jacobian uses local column indices, while the full-system
+         * Jacobian uses global column indices. Convert each local column to its
+         * corresponding global column so the two rows can be compared directly.
+         */
+        std::unordered_map<IdxT, RealT> sub_row_entries;
+
+        for (IdxT sub_index = sub_begin; sub_index < sub_end; ++sub_index)
+        {
+          const IdxT local_column = sub_cols[sub_index];
+
+          // Fail if the subsystem column index is outside the valid local range.
+          if (local_column >= num_sub_cols)
+          {
+            std::cout << "Invalid subsystem column index "
+                      << local_column
+                      << " in local row "
+                      << local_row << '\n';
+
+            matches = false;
+            continue;
+          }
+
+          const IdxT global_column = local_to_global[local_column];
+
+          // Fail if the subsystem row contains more than one entry for the same column.
+          if (sub_row_entries.find(global_column) != sub_row_entries.end())
+          {
+            std::cout << "Duplicate subsystem entry at ("
+                      << global_row << ", "
+                      << global_column << ")\n";
+
+            matches = false;
+            continue;
+          }
+
+          sub_row_entries[global_column] = sub_vals[sub_index];
+        }
+
+        // Compare full-system entries whose columns belong to the subsystem.
+        for (IdxT full_index = full_begin; full_index < full_end; ++full_index)
+        {
+          const IdxT global_column = full_cols[full_index];
+
+          // No need to check if column belongs outside the subsystem.
+          if (global_to_local.find(global_column) == global_to_local.end())
+          {
+            continue;
+          }
+
+          auto sub_entry = sub_row_entries.find(global_column);
+
+          // Then it must exist in the subsystem Jacobian; fail otherwise.
+          if (sub_entry == sub_row_entries.end())
+          {
+            std::cout << "Entry exists only in full Jacobian at ("
+                      << global_row << ", "
+                      << global_column << ")\n";
+
+            matches = false;
+            continue;
+          }
+
+          const RealT full_value = full_vals[full_index];
+          const RealT sub_value  = sub_entry->second;
+          const RealT difference = std::abs(full_value - sub_value);
+
+          // Different component evaluation orders can change the order of floating-point
+          // summation and introduce small roundoff differences. Use an appropriate
+          // tolerance when comparing against the monolithic reference.
+          // if tolerance is not supplied we simply use default machine precision provided
+          // by GridKit's Test::isEqual
+          auto isEqual = [&tolerance](RealT value, RealT reference)
+          {
+            if (tolerance)
+            {
+              return GridKit::Testing::isEqual(value, reference, *tolerance);
+            }
+
+            return GridKit::Testing::isEqual(value, reference);
+          };
+
+          // Then the values must agree, fail otherwise
+          if (!isEqual(sub_value, full_value))
+          {
+            std::cout << "Jacobian value mismatch at ("
+                      << global_row << ", "
+                      << global_column << "): "
+                      << "full = " << full_value
+                      << ", subsystem = " << sub_value
+                      << ", difference = " << difference << '\n';
+
+            matches = false;
+          }
+
+          // Remove the matched entry
+          sub_row_entries.erase(sub_entry);
+        }
+
+        // Any entries remaining must be missing from the
+        // full-system Jacobian, so this subsystem row contains incorrect entries, fail!
+        for (const auto& entry : sub_row_entries)
+        {
+          const IdxT global_column = entry.first;
+
+          std::cout << "Entry exists only in subsystem Jacobian at ("
+                    << global_row << ", "
+                    << global_column << ")\n";
+
+          matches = false;
+        }
+      }
+
+      return matches;
+    }
+
+  } // namespace Testing
+} // namespace GridKit
diff --git a/examples/PowerElectronics/Common/MicrogridNetwork.hpp b/examples/PowerElectronics/Common/MicrogridNetwork.hpp
new file mode 100644
index 000000000..84bc1b118
--- /dev/null
+++ b/examples/PowerElectronics/Common/MicrogridNetwork.hpp
@@ -0,0 +1,309 @@
+// MicrogridNetwork.hpp
+
+#pragma once
+
+#include <cassert>
+#include <cmath>
+#include <cstddef>
+#include <vector>
+
+#include <GridKit/Model/PowerElectronics/Bus/MicrogridBus.hpp>
+#include <GridKit/Model/PowerElectronics/Bus/SignalNode.hpp>
+#include <GridKit/Model/PowerElectronics/DistributedGenerator/DistributedGenerator.hpp>
+#include <GridKit/Model/PowerElectronics/MicrogridBusDQ/MicrogridBusDQ.hpp>
+#include <GridKit/Model/PowerElectronics/MicrogridLine/MicrogridLine.hpp>
+#include <GridKit/Model/PowerElectronics/MicrogridLoad/MicrogridLoad.hpp>
+#include <GridKit/Model/PowerElectronics/PartitionInterface/BusPartitionInterface.hpp>
+#include <GridKit/Model/PowerElectronics/SubsystemModel.hpp>
+#include <GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp>
+
+namespace GridKit
+{
+
+  using index_type = size_t;
+  using real_type  = double;
+
+  using SignalNode  = GridKit::PowerElectronics::SignalNode<real_type, index_type>;
+  using Bus         = GridKit::PowerElectronics::MicrogridBus<real_type, index_type>;
+  using BusDQ       = GridKit::MicrogridBusDQ<real_type, index_type>;
+  using DGGenerator = GridKit::DistributedGenerator<real_type, index_type>;
+  using Line        = GridKit::MicrogridLine<real_type, index_type>;
+  using Load        = GridKit::MicrogridLoad<real_type, index_type>;
+  using GenParams   = GridKit::DistributedGeneratorParameters<real_type, index_type>;
+
+  /*
+   * Contains components and nodes that make up the scaled microgrid network.
+   *
+   * The network contains 2 * N_size IBRs and stores the physical components
+   * that make up the scale microgrid.
+   */
+  struct ScaleMicrogridNetwork
+  {
+    size_t                    model_id_next;
+    size_t                    N_size;
+    SignalNode                dg_signal;
+    std::vector<Bus>          buses;
+    std::vector<BusDQ*>       busesDQ;
+    std::vector<DGGenerator*> generators;
+    std::vector<Line*>        lines;
+    std::vector<Load*>        loads;
+    std::vector<GenParams>    DGParam_list;
+
+    ScaleMicrogridNetwork(size_t n_size)
+      : model_id_next(0),
+        N_size(n_size),
+        buses(2 * n_size),
+        busesDQ(2 * n_size, nullptr),
+        generators(2 * n_size, nullptr),
+        lines(2 * n_size, nullptr),
+        loads(2 * n_size, nullptr),
+        DGParam_list(2 * n_size)
+    {
+    }
+  };
+
+  /**
+   * @brief Construct all components of a scaled microgrid network.
+   *
+   * Builds a microgrid containing @c 2*N_size generators and buses connected
+   * in a chain by transmission lines. Loads are connected to every other bus,
+   * and a virtual DQ bus is associated with each physical bus.
+   *
+   * The first two generators use the reference generator parameter set, while
+   * all remaining generators use the second parameter set. Line parameters
+   * alternate along the network, and the first load uses a different parameter
+   * set from the remaining loads.
+   *
+   * The created components are stored in @p network and assigned unique model
+   * identifiers using @c network.model_id_next.
+   *
+   * @param[in,out] network Network in which the microgrid components are
+   *                        constructed and stored.
+   *
+   * @pre @c network.N_size is greater than zero.
+   * @pre The component storage in @p network has been sized consistently with
+   *      @c network.N_size.
+   *
+   * @post @p network contains @c 2*N_size generators and virtual DQ buses.
+   * @post @p network contains @c 2*N_size-1 transmission lines connecting
+   *       consecutive buses.
+   * @post @p network contains @c N_size loads connected to every other bus.
+   * @post @c network.DGParam_list contains the parameters associated with each
+   *       generator.
+   * @post @c network.model_id_next is advanced for every component created.
+   *
+   * @note Components are dynamically allocated and their pointers are stored
+   *       in the corresponding network component vectors.
+   */
+  inline void buildScaleMicrogridNetwork(ScaleMicrogridNetwork& network)
+  {
+    size_t N_size = network.N_size;
+
+    assert(N_size > 0);
+    // Every Bus has the same virtual resistance. This is due to numerical stability as mentioned in the paper.
+    real_type RN = 1.0e4;
+
+    // DG Params Vector
+    // All DGs have the same set of parameters except for the first two.
+    GenParams DG_parms1;
+    DG_parms1.wb_  = 2.0 * M_PI * 50.0;
+    DG_parms1.wc_  = 31.41;
+    DG_parms1.mp_  = 9.4e-5;
+    DG_parms1.Vn_  = 380.0;
+    DG_parms1.nq_  = 1.3e-3;
+    DG_parms1.F_   = 0.75;
+    DG_parms1.Kiv_ = 420.0;
+    DG_parms1.Kpv_ = 0.1;
+    DG_parms1.Kic_ = 2.0e4;
+    DG_parms1.Kpc_ = 15.0;
+    DG_parms1.Cf_  = 5.0e-5;
+    DG_parms1.rLf_ = 0.1;
+    DG_parms1.Lf_  = 1.35e-3;
+    DG_parms1.rLc_ = 0.03;
+    DG_parms1.Lc_  = 0.35e-3;
+
+    GenParams DG_parms2;
+    DG_parms2.wb_  = 2.0 * M_PI * 50.0;
+    DG_parms2.wc_  = 31.41;
+    DG_parms2.mp_  = 12.5e-5;
+    DG_parms2.Vn_  = 380.0;
+    DG_parms2.nq_  = 1.5e-3;
+    DG_parms2.F_   = 0.75;
+    DG_parms2.Kiv_ = 390.0;
+    DG_parms2.Kpv_ = 0.05;
+    DG_parms2.Kic_ = 16.0e3;
+    DG_parms2.Kpc_ = 10.5;
+    DG_parms2.Cf_  = 50.0e-6;
+    DG_parms2.rLf_ = 0.1;
+    DG_parms2.Lf_  = 1.35e-3;
+    DG_parms2.rLc_ = 0.03;
+    DG_parms2.Lc_  = 0.35e-3;
+
+    network.DGParam_list.assign(2 * N_size, DG_parms2);
+
+    // First two generators use parameters 1
+    if (network.DGParam_list.size() >= 1)
+    {
+      network.DGParam_list[0] = DG_parms1;
+    }
+    if (network.DGParam_list.size() >= 2)
+    {
+      network.DGParam_list[1] = DG_parms1;
+    }
+
+    // line vector params
+    // Every odd line has the same parameters and every even line has the same parameters
+    real_type              rline1 = 0.23;
+    real_type              Lline1 = 0.1 / (2.0 * M_PI * 50.0);
+    real_type              rline2 = 0.35;
+    real_type              Lline2 = 0.58 / (2.0 * M_PI * 50.0);
+    std::vector<real_type> rline_list(2 * N_size - 1, 0.0);
+    std::vector<real_type> Lline_list(2 * N_size - 1, 0.0);
+    for (index_type i = 0; i < rline_list.size(); i++)
+    {
+      rline_list[i] = (i % 2) ? rline2 : rline1;
+      Lline_list[i] = (i % 2) ? Lline2 : Lline1;
+    }
+
+    // load parms
+    // Only the first load has the same paramaters.
+    real_type rload1 = 3.0;
+    real_type Lload1 = 2.0 / (2.0 * M_PI * 50.0);
+    real_type rload2 = 2.0;
+    real_type Lload2 = 1.0 / (2.0 * M_PI * 50.0);
+
+    std::vector<real_type> rload_list(N_size, rload2);
+    std::vector<real_type> Lload_list(N_size, Lload2);
+    if (rload_list.size() >= 1)
+    {
+      rload_list[0] = rload1;
+      Lload_list[0] = Lload1;
+    }
+
+    // Create the reference generator
+    auto* dg_ref = new DGGenerator(network.model_id_next++,
+                                   network.DGParam_list[0],
+                                   true,
+                                   &network.dg_signal,
+                                   &network.buses[0]);
+
+    network.generators[0] = dg_ref;
+
+    // Create the remaining generators.
+    for (index_type i = 1; i < 2 * N_size; i++)
+    {
+      auto* dg = new DGGenerator(network.model_id_next++,
+                                 network.DGParam_list[i],
+                                 false,
+                                 &network.dg_signal,
+                                 &network.buses[i]);
+
+      network.generators[i] = dg;
+    }
+
+    //  // Create transmission lines between consecutive buses.
+    for (index_type i = 0; i < 2 * N_size - 1; i++)
+    {
+      auto* line_model = new Line(network.model_id_next++,
+                                  rline_list[i],
+                                  Lline_list[i],
+                                  &network.dg_signal,
+                                  &network.buses[i],
+                                  &network.buses[i + 1]);
+
+      network.lines[i + 1] = line_model;
+    }
+
+    // Create loads on every other bus.
+    for (index_type i = 0; i < N_size; i++)
+    {
+      auto* load_model = new Load(network.model_id_next++,
+                                  rload_list[i],
+                                  Lload_list[i],
+                                  &network.dg_signal,
+                                  &network.buses[2 * i]);
+
+      network.loads[2 * i] = load_model;
+    }
+
+    // Create and Add all the microgrid Virtual DQ Buses
+    for (index_type i = 0; i < 2 * N_size; i++)
+    {
+      auto* virDQbus_model = new BusDQ(network.model_id_next++,
+                                       RN,
+                                       &network.buses[i]);
+
+      network.busesDQ[i] = virDQbus_model;
+    }
+  }
+
+  /**
+   * @brief Assemble a scaled microgrid network into a power electronics model.
+   *
+   * Adds the signal node, physical buses, generators, transmission lines,
+   * loads, and virtual DQ buses stored in @p network to @p sys_model.
+   *
+   * This function does not construct or allocate any network components. The
+   * physical network must already have been created by
+   * buildScaleMicrogridNetwork().
+   *
+   * @param[in] network Constructed scaled microgrid network whose components
+   *                    are added to the system model.
+   * @param[in,out] sys_model Power electronics model to which the network
+   *                          components and nodes are added.
+   *
+   * @pre @c network.N_size is greater than zero.
+   * @pre @p network has been constructed by buildScaleMicrogridNetwork().
+   * @pre All component and node pointers referenced by @p network are valid.
+   *
+   * @post The signal node and all physical buses in @p network have been added
+   *       to @p sys_model.
+   * @post All generators, transmission lines, loads, and virtual DQ buses in
+   *       @p network have been added to @p sys_model.
+   *
+   * @note This function only assembles the network into the system model. It
+   *       does not call PowerElectronicsModel::allocate().
+   */
+  inline void assembleSystem(ScaleMicrogridNetwork& network, GridKit::PowerElectronicsModel<real_type, index_type>& sys_model)
+  {
+    size_t N_size = network.N_size;
+
+    // Ensure minimum size requirement
+    assert(N_size > 0);
+
+    // Add all bus nodes
+    sys_model.addNode(&network.dg_signal);
+
+    for (size_t i = 0; i < 2 * N_size; i++)
+    {
+      sys_model.addNode(&network.buses[i]);
+    }
+
+    // Add all generators
+    for (index_type i = 0; i < 2 * N_size; i++)
+    {
+      sys_model.addComponent(network.generators[i]);
+    }
+
+    // Load all the Line components
+    for (index_type i = 1; i < 2 * N_size; i++)
+    {
+      sys_model.addComponent(network.lines[i]);
+    }
+
+    //  Load all the Load components
+    for (index_type i = 0; i < 2 * N_size; i++)
+    {
+      if (network.loads[i] != nullptr)
+      {
+        sys_model.addComponent(network.loads[i]);
+      }
+    }
+
+    // Add all the microgrid Virtual DQ Buses
+    for (index_type i = 0; i < 2 * N_size; i++)
+    {
+      sys_model.addComponent(network.busesDQ[i]);
+    }
+  }
+} // namespace GridKit
diff --git a/examples/PowerElectronics/Common/PartitionUtilities.hpp b/examples/PowerElectronics/Common/PartitionUtilities.hpp
new file mode 100644
index 000000000..d55119da4
--- /dev/null
+++ b/examples/PowerElectronics/Common/PartitionUtilities.hpp
@@ -0,0 +1,238 @@
+#pragma once
+
+#include <omp.h>
+
+#include <algorithm>
+#include <cassert>
+#include <cstddef>
+#include <vector>
+
+#include <GridKit/Model/PowerElectronics/MicrogridLine/MicrogridLine.hpp>
+#include <GridKit/Model/PowerElectronics/PartitionInterface/BusPartitionInterface.hpp>
+#include <GridKit/Model/PowerElectronics/SubsystemModel.hpp>
+
+#include "MicrogridNetwork.hpp"
+
+namespace GridKit
+{
+
+  /**
+   * @brief Partition a scaled microgrid network into subsystem models.
+   *
+   * Divides the network into contiguous groups of IBRs and creates one
+   * subsystem for each group. The IBRs are distributed as evenly as possible,
+   * with any remainder assigned to the first partitions.
+   *
+   * The reference signal node is added to the first subsystem. A partition
+   * interface is added at the right boundary of every subsystem except the
+   * last to represent its connection to the neighboring subsystem.
+   *
+   * @param[in,out] network Scaled microgrid network to partition.
+   * @param[out] subsystems Vector populated with the created subsystem models.
+   * @param[in] num_partitions Number of subsystems to create.
+   *
+   * @pre The network has been constructed and contains @c 2*N_size IBRs.
+   * @pre @p num_partitions is greater than zero and does not exceed the
+   *      number of IBRs in the network.
+   * @pre The generators, buses, virtual DQ buses, loads, and lines referenced
+   *      by the network have valid lifetimes for use by the created subsystems.
+   *
+   * @post @p subsystems contains exactly @p num_partitions subsystem models.
+   * @post Every IBR in the network belongs to exactly one subsystem.
+   * @post Every subsystem except the last has a partition interface at its
+   *       right boundary.
+   *
+   * @note The subsystem models are dynamically allocated. The caller is
+   *       responsible for releasing and deleting them.
+   */
+  template <class ScalarT, typename IdxT>
+  void partitionNetwork(
+      ScaleMicrogridNetwork&                                network,
+      std::vector<GridKit::SubsystemModel<ScalarT, IdxT>*>& subsystems,
+      size_t                                                num_partitions)
+  {
+    const size_t num_ibrs = 2 * network.N_size;
+
+    assert(num_partitions <= num_ibrs);
+
+    subsystems.resize(num_partitions);
+
+    IdxT q     = num_ibrs / num_partitions;
+    IdxT r     = num_ibrs % num_partitions;
+    IdxT index = 0;
+
+    for (IdxT j = 0; j < num_partitions; j++)
+    {
+      auto* partition =
+          new GridKit::SubsystemModel<ScalarT, IdxT>();
+
+      // Add the reference signal node to the first partition.
+      if (j == 0)
+      {
+        partition->addNode(&network.dg_signal);
+      }
+
+      IdxT part_size = q + (j < r ? 1 : 0);
+      IdxT end       = std::min(index + part_size, num_ibrs);
+
+      // Add all components belonging to this partition.
+      for (; index < end; ++index)
+      {
+        partition->addComponent(network.generators[index]);
+        partition->addComponent(network.busesDQ[index]);
+
+        if (network.loads[index] != nullptr)
+        {
+          partition->addComponent(network.loads[index]);
+        }
+
+        if (network.lines[index] != nullptr)
+        {
+          partition->addComponent(network.lines[index]);
+        }
+
+        partition->addNode(&network.buses[index]);
+      }
+
+      // Add the interface at the right boundary of the partition.
+      if (index < num_ibrs)
+      {
+
+        auto* busInterface = new GridKit::BusPartitionInterface<ScalarT, IdxT>(
+            &network.buses[index - 1],
+            network.lines[index],
+            network.model_id_next++);
+
+        busInterface->allocate();
+        partition->addInterface(busInterface);
+      }
+
+      subsystems[j] = partition;
+    }
+  }
+
+  /**
+   * @brief Evaluate subsystem residuals and reconstruct the global residual.
+   *
+   * Distributes the global state and state-derivative vectors to each
+   * subsystem using its internal and external index mappings. Each subsystem
+   * residual is then evaluated in parallel, and its internal residual entries
+   * are gathered into the global residual vector.
+   *
+   * @param[in] subsystems Subsystem models to evaluate.
+   * @param[in] y Global state vector.
+   * @param[in] yp Global state-derivative vector.
+   * @param[out] f Global residual vector reconstructed from the subsystem
+   *               residuals.
+   * @param[in] time Current simulation time.
+   * @param[in] alpha Jacobian scaling parameter associated with the current
+   *                  time-integration evaluation.
+   *
+   * @pre All subsystem models have been allocated.
+   * @pre @p y and @p yp contain all global entries
+   * @pre @p f is large enough to contain every global residual entry referenced
+   *      by the subsystems.
+   * @pre The subsystems have disjoint internal index sets so that parallel
+   *      writes to @p f do not overlap.
+   *
+   * @post Each subsystem residual has been evaluated at the supplied
+   *       @p time and @p alpha.
+   * @post @p f contains the reconstructed residual for all internal variables
+   *       represented by the subsystems.
+   */
+  template <class ScalarT, typename IdxT>
+  void evaluatePartitionResiduals(
+      const std::vector<GridKit::SubsystemModel<ScalarT, IdxT>*>& subsystems,
+      const std::vector<ScalarT>&                                 y,
+      const std::vector<ScalarT>&                                 yp,
+      std::vector<ScalarT>&                                       f,
+      ScalarT                                                     time,
+      ScalarT                                                     alpha)
+  {
+#ifdef _OPENMP
+#pragma omp parallel for schedule(guided)
+#endif
+    for (auto* partition : subsystems)
+    {
+      partition->updateTime(time, alpha);
+
+      for (size_t i = 0; i < partition->getExternSize(); i++)
+      {
+        partition->getExternalDataY()[i]  = y[partition->getExternalDataIndices()[i]];
+        partition->getExternalDataYP()[i] = yp[partition->getExternalDataIndices()[i]];
+      }
+
+      auto* partition_y  = partition->y().getData();
+      auto* partition_yp = partition->yp().getData();
+
+      for (size_t i = 0; i < partition->getInternalSize(); i++)
+      {
+        partition_y[i]  = y[partition->getNodeConnection(i)];
+        partition_yp[i] = yp[partition->getNodeConnection(i)];
+      }
+
+      partition->y().setDataUpdated();
+      partition->yp().setDataUpdated();
+
+      partition->evaluateResidual();
+
+      auto* residual = partition->getResidual().getData();
+
+      for (size_t i = 0; i < partition->getInternalSize(); i++)
+      {
+        f[partition->getNodeConnection(i)] = residual[i];
+      }
+    }
+  }
+
+  /**
+   * @brief Assemble the scaled microgrid from left to right.
+   *
+   * Adds the network components to @p sys_model in physical left-to-right
+   * order. For each bus location, the generator, virtual DQ bus, load, line,
+   * and bus associated with that location are added before moving to the next
+   * location in the network.
+   *
+   * This ordering is useful when comparing the monolithic system with
+   * partition-based evaluations that traverse the network from left to right.
+   *
+   * @param[in] network Constructed scaled microgrid network.
+   * @param[in,out] sys_model Power electronics model to populate.
+   *
+   * @pre @c network.N_size is greater than zero.
+   * @pre @p network has been constructed by buildScaleMicrogridNetwork().
+   * @pre All component and node pointers stored in @p network are valid.
+   *
+   * @post All nodes and components in @p network have been added to
+   *       @p sys_model in left-to-right network order.
+   */
+  template <class ScalarT, typename IdxT>
+  void assembleSystemLeftToRight(
+      ScaleMicrogridNetwork&                         network,
+      GridKit::PowerElectronicsModel<ScalarT, IdxT>& sys_model)
+  {
+    const size_t num_ibrs = 2 * network.N_size;
+
+    assert(network.N_size > 0);
+
+    sys_model.addNode(&network.dg_signal);
+
+    for (IdxT i = 0; i < num_ibrs; ++i)
+    {
+      sys_model.addComponent(network.generators[i]);
+      sys_model.addComponent(network.busesDQ[i]);
+
+      if (network.loads[i] != nullptr)
+      {
+        sys_model.addComponent(network.loads[i]);
+      }
+
+      if (network.lines[i] != nullptr)
+      {
+        sys_model.addComponent(network.lines[i]);
+      }
+
+      sys_model.addNode(&network.buses[i]);
+    }
+  }
+} // namespace GridKit
diff --git a/examples/PowerElectronics/Microgrid/CMakeLists.txt b/examples/PowerElectronics/Microgrid/CMakeLists.txt
index 6eae26be3..86dea6678 100644
--- a/examples/PowerElectronics/Microgrid/CMakeLists.txt
+++ b/examples/PowerElectronics/Microgrid/CMakeLists.txt
@@ -1,4 +1,9 @@
 add_executable(microgrid Microgrid.cpp)
+
+target_include_directories(
+  microgrid
+  PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/..)
+
 target_link_libraries(
   microgrid
   GridKit::power_elec_disgen
diff --git a/examples/PowerElectronics/Microgrid/Microgrid.cpp b/examples/PowerElectronics/Microgrid/Microgrid.cpp
index 74497fc41..5e844e91d 100644
--- a/examples/PowerElectronics/Microgrid/Microgrid.cpp
+++ b/examples/PowerElectronics/Microgrid/Microgrid.cpp
@@ -1,21 +1,13 @@
 #include <cmath>
-#include <filesystem>
-#include <fstream>
-#include <iomanip>
+#include <cstddef>
 #include <iostream>
-#include <memory>
-#include <numbers>
 
-#include <GridKit/Model/PowerElectronics/Bus/MicrogridBus.hpp>
-#include <GridKit/Model/PowerElectronics/Bus/SignalNode.hpp>
-#include <GridKit/Model/PowerElectronics/DistributedGenerator/DistributedGenerator.hpp>
-#include <GridKit/Model/PowerElectronics/MicrogridBusDQ/MicrogridBusDQ.hpp>
-#include <GridKit/Model/PowerElectronics/MicrogridLine/MicrogridLine.hpp>
-#include <GridKit/Model/PowerElectronics/MicrogridLoad/MicrogridLoad.hpp>
 #include <GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp>
 #include <GridKit/Solver/Dynamic/DynamicSolver.hpp>
 #include <GridKit/Solver/Dynamic/Ida.hpp>
 
+#include "Common/MicrogridNetwork.hpp"
+
 int main(int /* argc */, char const** /* argv */)
 {
   /// @todo Needs to be modified. Some components are small relative to others thus
@@ -28,139 +20,16 @@ int main(int /* argc */, char const** /* argv */)
   // Create model
   auto* sysmodel = new GridKit::PowerElectronicsModel<double, size_t>(use_jac);
 
-  // Modeled after the problem in the paper
-  double RN = 1.0e4;
-
-  // DG Params
-  static constexpr auto pi = std::numbers::pi_v<double>;
-
-  GridKit::DistributedGeneratorParameters<double, size_t> parms1;
-  parms1.wb_  = 2.0 * pi * 50.0;
-  parms1.wc_  = 31.41;
-  parms1.mp_  = 9.4e-5;
-  parms1.Vn_  = 380.0;
-  parms1.nq_  = 1.3e-3;
-  parms1.F_   = 0.75;
-  parms1.Kiv_ = 420.0;
-  parms1.Kpv_ = 0.1;
-  parms1.Kic_ = 2.0e4;
-  parms1.Kpc_ = 15.0;
-  parms1.Cf_  = 5.0e-5;
-  parms1.rLf_ = 0.1;
-  parms1.Lf_  = 1.35e-3;
-  parms1.rLc_ = 0.03;
-  parms1.Lc_  = 0.35e-3;
-
-  GridKit::DistributedGeneratorParameters<double, size_t> parms2;
-  // Parameters from MATLAB Microgrid code for first DG
-  parms2.wb_  = 2.0 * pi * 50.0;
-  parms2.wc_  = 31.41;
-  parms2.mp_  = 12.5e-5;
-  parms2.Vn_  = 380.0;
-  parms2.nq_  = 1.5e-3;
-  parms2.F_   = 0.75;
-  parms2.Kiv_ = 390.0;
-  parms2.Kpv_ = 0.05;
-  parms2.Kic_ = 16.0e3;
-  parms2.Kpc_ = 10.5;
-  parms2.Cf_  = 50.0e-6;
-  parms2.rLf_ = 0.1;
-  parms2.Lf_  = 1.35e-3;
-  parms2.rLc_ = 0.03;
-  parms2.Lc_  = 0.35e-3;
-
-  // Line params
-  double rline1 = 0.23;
-  double Lline1 = 0.1 / (2.0 * pi * 50.0);
-
-  double rline2 = 0.35;
-  double Lline2 = 0.58 / (2.0 * pi * 50.0);
-
-  double rline3 = 0.23;
-  double Lline3 = 0.1 / (2.0 * pi * 50.0);
-
-  // load parms
-  double rload1 = 3.0;
-  double Lload1 = 2.0 / (2.0 * pi * 50.0);
-
-  double rload2 = 2.0;
-  double Lload2 = 1.0 / (2.0 * pi * 50.0);
-
-  using SignalNode = GridKit::PowerElectronics::SignalNode<double, size_t>;
-  SignalNode dg_signal;
-
-  sysmodel->addNode(&dg_signal);
-
-  using Bus = GridKit::PowerElectronics::MicrogridBus<double, size_t>;
-  Bus bus1;
-  Bus bus2;
-  Bus bus3;
-  Bus bus4;
-
-  sysmodel->addNode(&bus1);
-  sysmodel->addNode(&bus2);
-  sysmodel->addNode(&bus3);
-  sysmodel->addNode(&bus4);
-
-  // dg 1
-  GridKit::DistributedGenerator<double, size_t>* dg1 = new GridKit::DistributedGenerator<double, size_t>(
-      0, parms1, true, &dg_signal, &bus1);
-  sysmodel->addComponent(dg1);
-
-  // dg 2
-  GridKit::DistributedGenerator<double, size_t>* dg2 = new GridKit::DistributedGenerator<double, size_t>(
-      1, parms1, false, &dg_signal, &bus2);
-  sysmodel->addComponent(dg2);
-
-  // dg 3
-  GridKit::DistributedGenerator<double, size_t>* dg3 = new GridKit::DistributedGenerator<double, size_t>(
-      2, parms2, false, &dg_signal, &bus3);
-  sysmodel->addComponent(dg3);
-
-  // dg 4
-  GridKit::DistributedGenerator<double, size_t>* dg4 = new GridKit::DistributedGenerator<double, size_t>(
-      3, parms2, false, &dg_signal, &bus4);
-  sysmodel->addComponent(dg4);
-
-  // Lines
-
-  // line 1
-  GridKit::MicrogridLine<double, size_t>* l1 = new GridKit::MicrogridLine<double, size_t>(
-      4, rline1, Lline1, &dg_signal, &bus1, &bus2);
-  sysmodel->addComponent(l1);
-
-  // line 2
-  GridKit::MicrogridLine<double, size_t>* l2 = new GridKit::MicrogridLine<double, size_t>(
-      5, rline2, Lline2, &dg_signal, &bus2, &bus3);
-  sysmodel->addComponent(l2);
-
-  // line 3
-  GridKit::MicrogridLine<double, size_t>* l3 = new GridKit::MicrogridLine<double, size_t>(
-      6, rline3, Lline3, &dg_signal, &bus3, &bus4);
-  sysmodel->addComponent(l3);
-
-  //  loads
-
-  // load 1
-  GridKit::MicrogridLoad<double, size_t>* load1 = new GridKit::MicrogridLoad<double, size_t>(7, rload1, Lload1, &dg_signal, &bus1);
-  sysmodel->addComponent(load1);
-
-  // load 2
-  GridKit::MicrogridLoad<double, size_t>* load2 = new GridKit::MicrogridLoad<double, size_t>(8, rload2, Lload2, &dg_signal, &bus3);
-  sysmodel->addComponent(load2);
-
-  // Virtual PQ Buses
-  GridKit::MicrogridBusDQ<double, size_t>* bus_para_1 = new GridKit::MicrogridBusDQ<double, size_t>(9, RN, &bus1);
-  sysmodel->addComponent(bus_para_1);
-
-  GridKit::MicrogridBusDQ<double, size_t>* bus_para_2 = new GridKit::MicrogridBusDQ<double, size_t>(10, RN, &bus2);
-  sysmodel->addComponent(bus_para_2);
+  // Build the four-generator microgrid network.
+  size_t                         N_size = 2;
+  GridKit::ScaleMicrogridNetwork network(N_size);
 
-  GridKit::MicrogridBusDQ<double, size_t>* bus_para_3 = new GridKit::MicrogridBusDQ<double, size_t>(11, RN, &bus3);
-  sysmodel->addComponent(bus_para_3);
+  GridKit::buildScaleMicrogridNetwork(network);
+  GridKit::assembleSystem(network, *sysmodel);
 
-  GridKit::MicrogridBusDQ<double, size_t>* bus_para_4 = new GridKit::MicrogridBusDQ<double, size_t>(12, RN, &bus4);
-  sysmodel->addComponent(bus_para_4);
+  // Generator parameters used to construct the initial conditions.
+  const auto& parms1 = network.DGParam_list[0];
+  const auto& parms2 = network.DGParam_list[2];
 
   sysmodel->allocate();
 
@@ -192,7 +61,7 @@ int main(int /* argc */, char const** /* argv */)
   }
 
   // since the intial P_com = 0
-  y[dg_signal.getNodeConnection(0).idx_] = parms1.wb_;
+  y[network.dg_signal.getNodeConnection(0).idx_] = parms1.wb_;
 
   sysmodel->y().setDataUpdated();
   sysmodel->yp().setDataUpdated();
diff --git a/examples/PowerElectronics/Partition/CMakeLists.txt b/examples/PowerElectronics/Partition/CMakeLists.txt
new file mode 100644
index 000000000..cd182dcd3
--- /dev/null
+++ b/examples/PowerElectronics/Partition/CMakeLists.txt
@@ -0,0 +1,39 @@
+find_package(OpenMP)
+
+add_executable(PartitionMicrogrid PartitionMicrogrid.cpp)
+target_include_directories(
+  PartitionMicrogrid
+  PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/..)
+
+target_link_libraries(
+  PartitionMicrogrid
+  PRIVATE GridKit::power_elec_disgen
+          GridKit::power_elec_microline
+          GridKit::power_elec_microload
+          GridKit::solvers_dyn
+          GridKit::power_elec_microbusdq
+          GridKit::power_elec_partition_interfaces)
+
+add_executable(PartitionScaleMicrogrid PartitionScaleMicrogrid.cpp)
+target_include_directories(
+  PartitionScaleMicrogrid
+  PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/..)
+
+target_link_libraries(
+  PartitionScaleMicrogrid
+  PRIVATE GridKit::power_elec_disgen
+          GridKit::power_elec_microline
+          GridKit::power_elec_microload
+          GridKit::solvers_dyn
+          GridKit::power_elec_microbusdq
+          GridKit::power_elec_partition_interfaces)
+
+if(OpenMP_CXX_FOUND)
+  target_link_libraries(PartitionScaleMicrogrid PRIVATE OpenMP::OpenMP_CXX)
+endif()
+
+add_test(NAME PartitionMicrogrid COMMAND PartitionMicrogrid)
+add_test(NAME PartitionScaleMicrogrid COMMAND PartitionScaleMicrogrid)
+
+install(TARGETS PartitionMicrogrid RUNTIME DESTINATION bin)
+install(TARGETS PartitionScaleMicrogrid RUNTIME DESTINATION bin)
diff --git a/examples/PowerElectronics/Partition/PartitionMicrogrid.cpp b/examples/PowerElectronics/Partition/PartitionMicrogrid.cpp
new file mode 100644
index 000000000..72e638ca0
--- /dev/null
+++ b/examples/PowerElectronics/Partition/PartitionMicrogrid.cpp
@@ -0,0 +1,295 @@
+#include <cmath>
+#include <cstddef>
+#include <cstdio>
+#include <cstdlib>
+#include <iostream>
+#include <limits>
+#include <vector>
+
+#include <GridKit/Model/PowerElectronics/Bus/MicrogridBus.hpp>
+#include <GridKit/Model/PowerElectronics/Bus/SignalNode.hpp>
+#include <GridKit/Model/PowerElectronics/DistributedGenerator/DistributedGenerator.hpp>
+#include <GridKit/Model/PowerElectronics/MicrogridBusDQ/MicrogridBusDQ.hpp>
+#include <GridKit/Model/PowerElectronics/MicrogridLine/MicrogridLine.hpp>
+#include <GridKit/Model/PowerElectronics/MicrogridLoad/MicrogridLoad.hpp>
+#include <GridKit/Model/PowerElectronics/PartitionInterface/BusPartitionInterface.hpp>
+#include <GridKit/Model/PowerElectronics/SubsystemModel.hpp>
+#include <GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp>
+#include <GridKit/Solver/Dynamic/DynamicSolver.hpp>
+#include <GridKit/Solver/Dynamic/Ida.hpp>
+
+#include "Common/JacTestHelper.hpp"
+#include "Common/MicrogridNetwork.hpp"
+#include "Common/PartitionUtilities.hpp"
+
+using Component = GridKit::CircuitComponent<double, size_t>;
+using Node      = GridKit::PowerElectronics::NodeBase<double, size_t>;
+using Subsystem = GridKit::SubsystemModel<double, size_t>;
+using System    = GridKit::PowerElectronicsModel<double, size_t>;
+
+std::vector<size_t> getComponentConnections(const std::vector<Component*>& components);
+std::vector<size_t> getNodeConnections(const std::vector<Node*>& nodes);
+
+/*
+ * Verify partitioned residual and Jacobian evaluation against the
+ *        monolithic microgrid model.
+ *
+ * The microgrid is manually divided into two subsystems. A partition interface
+ * is introduced at the boundary between bus 1 and the line connecting buses 1
+ * and 2. Each subsystem is evaluated independently, and its residuals are gathered
+ * and compared with the monolithic reference.
+ */
+int main()
+{
+
+  constexpr size_t num_network_sections = 2;
+  constexpr double time                 = 0.1;
+  constexpr double alpha                = 0.1;
+
+  bool use_jac = true;
+
+  // ---------------------------------------------------------------------------
+  // build the grid network and assemble the system model
+  // ---------------------------------------------------------------------------
+  GridKit::ScaleMicrogridNetwork network(num_network_sections);
+
+  GridKit::buildScaleMicrogridNetwork(network);
+
+  auto* system = new System(use_jac);
+
+  GridKit::assembleSystemLeftToRight(network, *system);
+  system->allocate();
+
+  std::vector<double> y(system->size());
+  std::vector<double> yp(system->size());
+
+  for (size_t i = 0; i < system->size(); ++i)
+  {
+    y[i]  = static_cast<double>(i + 1);
+    yp[i] = static_cast<double>(i + 1);
+  }
+
+  auto* system_y  = system->y().getData();
+  auto* system_yp = system->yp().getData();
+
+  for (size_t i = 0; i < system->size(); ++i)
+  {
+    system_y[i]  = y[i];
+    system_yp[i] = yp[i];
+  }
+
+  system->y().setDataUpdated();
+  system->yp().setDataUpdated();
+
+  system->updateTime(time, alpha);
+  system->evaluateResidual();
+  system->evaluateJacobian();
+
+  auto* system_jacobian = system->getCsrJacobian();
+  auto* system_residual = system->getResidual().getData();
+
+  //------------------------------------------------------------------------------
+  // Gather all global indices belonging partition 1 and 2 to test release() later
+  //------------------------------------------------------------------------------
+  std::vector<Component*> components = {
+      network.generators[0],
+      network.generators[1],
+      network.lines[1],
+      network.loads[0],
+      network.busesDQ[0],
+      network.busesDQ[1],
+      network.generators[2],
+      network.generators[3],
+      network.lines[2],
+      network.lines[3],
+      network.loads[2],
+      network.busesDQ[2],
+      network.busesDQ[3]};
+
+  std::vector<Node*> nodes = {
+      &network.dg_signal,
+      &network.buses[0],
+      &network.buses[1],
+      &network.buses[2],
+      &network.buses[3]};
+
+  const auto original_component_connections = getComponentConnections(components);
+  const auto original_node_connections      = getNodeConnections(nodes);
+
+  // --------------------------------------------------------------------------------
+  // Create 2 Partitions and Partition interfaces
+  // --------------------------------------------------------------------------------
+
+  auto* bus_interface = new GridKit::BusPartitionInterface<double, size_t>(
+      &network.buses[1],
+      network.lines[2],
+      14);
+
+  bus_interface->allocate();
+
+  auto* partition1 = new Subsystem();
+  auto* partition2 = new Subsystem();
+
+  // --------------------------------------------------------------------------------
+  // Manually add components, nodes and a bus partition interface to Partition 1
+  // --------------------------------------------------------------------------------
+
+  partition1->addNode(&network.dg_signal);
+  partition1->addComponent(network.generators[0]);
+  partition1->addComponent(network.busesDQ[0]);
+  partition1->addComponent(network.loads[0]);
+  partition1->addNode(&network.buses[0]);
+  partition1->addComponent(network.lines[1]);
+  partition1->addComponent(network.generators[1]);
+  partition1->addComponent(network.busesDQ[1]);
+  partition1->addInterface(bus_interface);
+  partition1->addNode(&network.buses[1]);
+
+  // ---------------------------------------------------------------------------
+  // Manually add components and nodes to Partition 2
+  // ---------------------------------------------------------------------------
+
+  partition2->addComponent(network.generators[2]);
+  partition2->addComponent(network.busesDQ[2]);
+  partition2->addComponent(network.loads[2]);
+  partition2->addComponent(network.lines[2]);
+  partition2->addNode(&network.buses[2]);
+  partition2->addComponent(network.generators[3]);
+  partition2->addComponent(network.busesDQ[3]);
+  partition2->addComponent(network.lines[3]);
+  partition2->addNode(&network.buses[3]);
+
+  std::vector<Subsystem*> partitions = {partition1, partition2};
+
+  for (auto* partition : partitions)
+  {
+    partition->allocate();
+  }
+
+  std::vector<double> partition_residual(system->size(), 0.0);
+
+  GridKit::evaluatePartitionResiduals(partitions, y, yp, partition_residual, time, alpha);
+
+  // ---------------------------------------------------------------------------
+  // Verify the subsystem Jacobians
+  // ---------------------------------------------------------------------------
+
+  bool jacobians_match = true;
+
+  for (auto* partition : partitions)
+  {
+    partition->evaluateJacobian();
+
+    auto* partition_jacobian = partition->getCsrJacobian();
+
+    jacobians_match = jacobians_match && GridKit::Testing::verifySubsystemJacobian(*system_jacobian, *partition_jacobian, *partition);
+  }
+
+  if (!jacobians_match)
+  {
+    std::cout << "ERROR: At least one subsystem Jacobian is incorrect!\n";
+    return 1;
+  }
+
+  // ---------------------------------------------------------------------------
+  // Gather and verify the subsystem residuals against monolithic residuals
+  // ---------------------------------------------------------------------------
+
+  double max_error = 0.0;
+
+  for (size_t i = 0; i < system->size(); ++i)
+  {
+    const double error = std::abs(system_residual[i] - partition_residual[i]) / std::abs(system_residual[i] + 1);
+
+    max_error = std::max(max_error, error);
+  }
+
+  const double machine_epsilon =
+      std::numeric_limits<double>::epsilon();
+
+  const bool residuals_match = max_error <= machine_epsilon;
+
+  std::cout << "\nPartition Microgrid Validation\n";
+  std::cout << "------------------------------\n";
+
+  std::cout << std::left
+            << std::setw(32) << "Maximum residual error:"
+            << std::setprecision(16)
+            << max_error << '\n';
+
+  std::cout << std::left
+            << std::setw(32) << "Machine epsilon:"
+            << machine_epsilon << '\n';
+
+  std::cout << std::left
+            << std::setw(32) << "Residuals matched:"
+            << (residuals_match ? "True" : "False")
+            << '\n';
+
+  // ---------------------------------------------------------------------------
+  // Verify subsystem release() from SubsystemModel
+  // ---------------------------------------------------------------------------
+
+  for (auto* partition : partitions)
+  {
+    partition->release();
+  }
+
+  const bool components_restored = getComponentConnections(components) == original_component_connections;
+  const bool nodes_restored      = getNodeConnections(nodes) == original_node_connections;
+
+  if (!components_restored || !nodes_restored)
+  {
+    std::cout << "ERROR: Subsystem release did not restore "
+                 "the original global connection indices!\n";
+
+    return 1;
+  }
+
+  // ---------------------------------------------------------------------------
+  // Clean up
+  // ---------------------------------------------------------------------------
+  delete system;
+
+  for (auto* partition : partitions)
+  {
+    delete partition;
+  }
+
+  return residuals_match ? 0 : 1;
+}
+
+/**
+ * @brief Collect the connection indices of a set of components.
+ */
+std::vector<size_t> getComponentConnections(const std::vector<Component*>& components)
+{
+  std::vector<size_t> connections;
+
+  for (const auto* component : components)
+  {
+    for (size_t i = 0; i < component->size(); ++i)
+    {
+      connections.push_back(component->getNodeConnection(i));
+    }
+  }
+
+  return connections;
+}
+
+/**
+ * @brief Collect the connection indices of a set of nodes.
+ */
+std::vector<size_t> getNodeConnections(const std::vector<Node*>& nodes)
+{
+  std::vector<size_t> connections;
+
+  for (auto* node : nodes)
+
+    for (size_t i = 0; i < node->size(); ++i)
+    {
+      connections.push_back(node->getNodeConnection(i).idx_);
+    }
+
+  return connections;
+}
diff --git a/examples/PowerElectronics/Partition/PartitionScaleMicrogrid.cpp b/examples/PowerElectronics/Partition/PartitionScaleMicrogrid.cpp
new file mode 100644
index 000000000..1834cd2e3
--- /dev/null
+++ b/examples/PowerElectronics/Partition/PartitionScaleMicrogrid.cpp
@@ -0,0 +1,428 @@
+#include <omp.h>
+
+#include <cassert>
+#include <chrono>
+#include <cmath>
+#include <cstddef>
+#include <format>
+#include <iostream>
+#include <string>
+
+#include <GridKit/Model/PowerElectronics/Bus/MicrogridBus.hpp>
+#include <GridKit/Model/PowerElectronics/Bus/SignalNode.hpp>
+#include <GridKit/Model/PowerElectronics/DistributedGenerator/DistributedGenerator.hpp>
+#include <GridKit/Model/PowerElectronics/MicrogridBusDQ/MicrogridBusDQ.hpp>
+#include <GridKit/Model/PowerElectronics/MicrogridLine/MicrogridLine.hpp>
+#include <GridKit/Model/PowerElectronics/MicrogridLoad/MicrogridLoad.hpp>
+#include <GridKit/Model/PowerElectronics/PartitionInterface/BusPartitionInterface.hpp>
+#include <GridKit/Model/PowerElectronics/SubsystemModel.hpp>
+#include <GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp>
+#include <GridKit/Solver/Dynamic/DynamicSolver.hpp>
+#include <GridKit/Solver/Dynamic/Ida.hpp>
+#include <GridKit/Testing/Testing.hpp>
+
+#include "Common/JacTestHelper.hpp"
+#include "Common/MicrogridNetwork.hpp"
+#include "Common/PartitionUtilities.hpp"
+
+using index_type = size_t;
+using real_type  = double;
+
+using SignalNode  = GridKit::PowerElectronics::SignalNode<real_type, index_type>;
+using Bus         = GridKit::PowerElectronics::MicrogridBus<real_type, index_type>;
+using BusDQ       = GridKit::MicrogridBusDQ<real_type, index_type>;
+using DGGenerator = GridKit::DistributedGenerator<real_type, index_type>;
+using Line        = GridKit::MicrogridLine<real_type, index_type>;
+using Load        = GridKit::MicrogridLoad<real_type, index_type>;
+
+/*
+ * Output data for each parallel function evaluation run
+ */
+struct RunResult
+{
+  bool        success = true;
+  index_type  num_partitions;
+  real_type   partition_eval_time;  // seconds
+  real_type   monolithic_eval_time; // seconds
+  real_type   speedup;              // monolithic / partition
+  real_type   max_error;
+  std::string jacobian_status;
+};
+
+/*
+ * Reference data from the monolithic system evaluation.
+ *
+ * The monolithic system is evaluated once and its states, residual, and
+ * evaluation time are stored here. Each partition configuration reuses this
+ * data for validation and performance comparison.
+ */
+struct MonolithicReference
+{
+  std::vector<real_type> y;
+  std::vector<real_type> yp;
+  std::vector<real_type> residual;
+  real_type              monolithic_eval_time;
+};
+
+/**
+ * @brief Evaluate the monolithic system and store reference data.
+ *
+ * Initializes the global state and state-derivative vectors with deterministic
+ * values, evaluates and times the monolithic residual, evaluates the full
+ * Jacobian, and stores the resulting data for later comparison with
+ * partitioned evaluations.
+ *
+ * The residual evaluation time is measured independently from the Jacobian
+ * evaluation so that the stored timing represents only the cost of the
+ * monolithic residual evaluation.
+ *
+ * @param[in,out] system Allocated monolithic power electronics system.
+ *
+ * @return MonolithicReference containing the state vector, state-derivative
+ *         vector, residual, and residual evaluation time.
+ *
+ * @pre @p system is non-null.
+ * @pre @p system has already been assembled and allocated.
+ *
+ * @post The state and state-derivative vectors of @p system contain the
+ *       deterministic reference values.
+ * @post The monolithic residual and Jacobian have been evaluated.
+ * @post The returned reference contains an independent copy of the
+ *       monolithic residual.
+ */
+MonolithicReference evaluateMonolithicSystem(GridKit::PowerElectronicsModel<real_type, index_type>* system)
+{
+  MonolithicReference reference;
+
+  // ---------------------------------------------------------------------------
+  // Initialize the reference state and state-derivative vectors
+  // ---------------------------------------------------------------------------
+  reference.y.resize(system->size());
+  reference.yp.resize(system->size());
+
+  for (size_t i = 0; i < system->size(); ++i)
+  {
+    reference.y[i]  = static_cast<real_type>(i + 1);
+    reference.yp[i] = static_cast<real_type>(i + 1);
+  }
+
+  // Copy the reference values into the vectors owned by the system model.
+  auto* system_y  = system->y().getData();
+  auto* system_yp = system->yp().getData();
+
+  for (size_t i = 0; i < system->size(); ++i)
+  {
+    system_y[i]  = reference.y[i];
+    system_yp[i] = reference.yp[i];
+  }
+
+  system->y().setDataUpdated();
+  system->yp().setDataUpdated();
+
+  // ---------------------------------------------------------------------------
+  // Evaluate and time the monolithic residual
+  // ---------------------------------------------------------------------------
+
+  auto start_time = std::chrono::high_resolution_clock::now();
+
+  system->updateTime(0.1, 0.1);
+  system->evaluateResidual();
+
+  auto end_time = std::chrono::high_resolution_clock::now();
+
+  reference.monolithic_eval_time = std::chrono::duration<real_type>(end_time - start_time).count();
+
+  // evaluate the monolithic Jacobian
+  system->evaluateJacobian();
+
+  // Store a copy of the residual return reference data
+  auto* residual = system->getResidual().getData();
+
+  reference.residual.assign(residual, residual + system->size());
+
+  return reference;
+}
+
+/**
+ * @brief Evaluate one partition configuration and validate it against the
+ *        monolithic reference system.
+ *
+ * Creates the requested subsystem decomposition, allocates each subsystem,
+ * evaluates the partitioned residual, verifies each subsystem Jacobian against
+ * the monolithic Jacobian, and compares the reconstructed global residual with
+ * the stored monolithic residual.
+ *
+ * The partitioned residual evaluation is timed independently and compared with
+ * the previously measured monolithic residual evaluation time to compute the
+ * resulting speedup.
+ *
+ * @param[in,out] network Scaled microgrid network to partition.
+ * @param[in] system Allocated monolithic system used as the reference.
+ * @param[in] reference Stored monolithic state, derivative, residual, and
+ *                      timing information.
+ * @param[in] num_partitions Number of subsystem partitions to create.
+ *
+ * @return Performance and validation results for the requested partition count.
+ *
+ * @pre @p system is non-null and has already been evaluated.
+ * @pre @p network has been constructed using buildScaleMicrogridNetwork().
+ * @pre @p num_partitions is greater than zero and does not exceed the number
+ *      of IBRs in the network.
+ * @pre @p reference contains state and derivative vectors consistent with the
+ *      size of @p system.
+ *
+ * @post All temporary subsystem models created by this function are released
+ *       and deleted before returning.
+ * @post The returned result reports the partition timing, speedup, residual
+ *       error, and subsystem Jacobian validation status.
+ */
+RunResult evaluatePartitioning(
+    GridKit::ScaleMicrogridNetwork&                        network,
+    GridKit::PowerElectronicsModel<real_type, index_type>* system,
+    const MonolithicReference&                             reference,
+    index_type                                             num_partitions)
+{
+
+  // ---------------------------------------------------------------------------
+  // Create and allocate subsystem partitions
+  // ---------------------------------------------------------------------------
+
+  std::vector<GridKit::SubsystemModel<real_type, index_type>*> subsystems;
+
+  GridKit::partitionNetwork(network, subsystems, num_partitions);
+
+  for (auto* partition : subsystems)
+  {
+    partition->allocate();
+  }
+
+  // Global residual reconstructed from the subsystem residuals.
+  std::vector<real_type> f(system->size(), 1.0);
+  // Elementwise error between the monolithic and reconstructed residuals.
+  std::vector<real_type> error(system->size(), 1.0);
+
+  // ---------------------------------------------------------------------------
+  // Evaluate and time the partitioned residual
+  // ---------------------------------------------------------------------------
+  auto start_time = std::chrono::high_resolution_clock::now();
+
+  GridKit::evaluatePartitionResiduals(subsystems, reference.y, reference.yp, f, 0.1, 0.1);
+
+  auto end_time = std::chrono::high_resolution_clock::now();
+
+  auto partition_eval_time = std::chrono::duration<real_type>(end_time - start_time);
+
+  // ---------------------------------------------------------------------------
+  // Verify subsystem Jacobians against the monolithic Jacobian
+  // ---------------------------------------------------------------------------
+
+  auto* system_jacobian = system->getCsrJacobian();
+
+  bool jacobian_match = true;
+
+  for (auto* partition : subsystems)
+  {
+    partition->evaluateJacobian();
+
+    jacobian_match = jacobian_match && GridKit::Testing::verifySubsystemJacobian(*system_jacobian, *partition->getCsrJacobian(), *partition);
+  }
+
+  // ---------------------------------------------------------------------------
+  // Compare the reconstructed and monolithic residuals
+  // ---------------------------------------------------------------------------
+  real_type max_error = 0.0;
+
+  for (size_t i = 0; i < system->size(); ++i)
+  {
+    error[i] = std::abs(f[i] - reference.residual[i]) / (reference.residual[i] + 1.0);
+
+    if (max_error < error[i])
+    {
+      max_error = error[i];
+    }
+  }
+
+  // ---------------------------------------------------------------------------
+  // Store performance and validation results
+  // ---------------------------------------------------------------------------
+  RunResult result;
+
+  result.num_partitions       = num_partitions;
+  result.partition_eval_time  = partition_eval_time.count();
+  result.monolithic_eval_time = reference.monolithic_eval_time;
+  result.speedup              = reference.monolithic_eval_time / partition_eval_time.count();
+  result.max_error            = max_error;
+  result.jacobian_status      = jacobian_match ? "Correct" : "Wrong";
+
+  // ---------------------------------------------------------------------------
+  // Check validation results
+  // ---------------------------------------------------------------------------
+  if (!jacobian_match)
+  {
+    std::cout << "ERROR: At least one subsystem Jacobian is incorrect!"
+              << std::endl;
+
+    result.success = false;
+  }
+
+  if (max_error > std::numeric_limits<real_type>::epsilon())
+  {
+    std::cout << "ERROR: Max Error too high!: " << max_error << std::endl;
+    result.success = false;
+  }
+
+  // ---------------------------------------------------------------------------
+  // Release and destroy temporary subsystem models
+  // ---------------------------------------------------------------------------
+
+  for (auto* partition : subsystems)
+  {
+    partition->release();
+    delete partition;
+  }
+
+  return result;
+}
+
+/**
+ * @brief Benchmark and validate partitioned residual evaluation of a large
+ *        scaled microgrid. It also confirms the correctness of the subsystem
+ *        Jacobian.
+ *
+ * Builds a scaled microgrid once, constructs a monolithic reference system,
+ * and compares several subsystem decompositions of the same physical network.
+ *
+ * For each requested partition count, the benchmark measures the partitioned
+ * residual evaluation time, computes its speedup relative to the monolithic
+ * evaluation, verifies the reconstructed residual, and checks every subsystem
+ * Jacobian against the monolithic Jacobian.
+ *
+ * The physical network and monolithic system are constructed only once.
+ * Individual subsystem models are created and destroyed separately for each
+ * partition count.
+ *
+ * @return 0 if all residual and Jacobian validation checks pass; otherwise 1.
+ */
+int main(int argc, char const* argv[])
+{
+  index_type N_size = 5000;
+
+  std::vector<index_type> num_partitions_list = {500};
+
+  /*
+   * If command-line arguments are provided, the first argument specifies
+   * N_size and all remaining arguments specify partition counts.
+   *
+   * Example:
+   *   ./PartitionScaleMicrogrid 5000 10 48 100 500
+   */
+  if (argc > 1)
+  {
+    try
+    {
+      N_size = static_cast<index_type>(std::stoull(argv[1]));
+
+      if (N_size < 1)
+      {
+        std::cerr << "ERROR: N_size must be at least 1.\n";
+        return 1;
+      }
+
+      // When N_size is supplied explicitly, at least one partition count
+      // must also be supplied.
+      if (argc < 3)
+      {
+        std::cerr << "ERROR: At least one partition count must be provided "
+                  << "when N_size is specified.\n";
+        return 1;
+      }
+
+      num_partitions_list.clear();
+
+      for (int i = 2; i < argc; ++i)
+      {
+        index_type num_partitions = static_cast<index_type>(std::stoull(argv[i]));
+
+        if (num_partitions < 1)
+        {
+          std::cerr << "ERROR: Number of partitions must be at least 1.\n";
+          return 1;
+        }
+
+        if (num_partitions > 2 * N_size)
+        {
+          std::cerr << "ERROR: Number of partitions ("
+                    << num_partitions
+                    << ") cannot exceed the number of IBRs ("
+                    << 2 * N_size
+                    << ").\n";
+          return 1;
+        }
+
+        num_partitions_list.push_back(num_partitions);
+      }
+    }
+    catch (const std::exception& e)
+    {
+      std::cerr
+          << "ERROR: Invalid command-line argument: "
+          << e.what()
+          << "\n";
+
+      return 1;
+    }
+  }
+
+  bool use_jac = true;
+
+  // Build the physical network once.
+  GridKit::ScaleMicrogridNetwork network(N_size);
+  GridKit::buildScaleMicrogridNetwork(network);
+
+  // Build, assemble and allocate the monolithic system once.
+  auto* system = new GridKit::PowerElectronicsModel<real_type, index_type>(use_jac);
+
+  GridKit::assembleSystemLeftToRight(network, *system);
+  system->allocate();
+
+  // Evaluate the monolithic reference once.
+  MonolithicReference reference = evaluateMonolithicSystem(system);
+
+  std::cout << std::format("{:<16}{:>16}{:>18}{:>12}{:>14}{:>16}\n",
+                           "num_partitions",
+                           "partition_time",
+                           "monolithic_time",
+                           "speedup",
+                           "error",
+                           "Jacobians");
+
+  std::cout << std::string(93, '-') << "\n";
+
+  // Only the partitioned system is rebuilt and evaluated for each partition count.
+  for (index_type p : num_partitions_list)
+  {
+    assert(p <= 2 * N_size);
+
+    // Takes in the network and partition it into p partitions, and performs parallel function eval
+    RunResult r = evaluatePartitioning(network, system, reference, p);
+
+    if (!r.success)
+    {
+      delete system;
+      return 1;
+    }
+
+    // Output the results from partition evaluation
+    std::cout << std::format("{:<16d}{:>14.4f} s{:>16.4f} s{:>11.2f}x{:>14.3e} {:>16s}\n",
+                             r.num_partitions,
+                             r.partition_eval_time,
+                             r.monolithic_eval_time,
+                             r.speedup,
+                             r.max_error,
+                             r.jacobian_status);
+  }
+
+  delete system;
+
+  return 0;
+}
diff --git a/examples/PowerElectronics/ScaleMicrogrid/CMakeLists.txt b/examples/PowerElectronics/ScaleMicrogrid/CMakeLists.txt
index 0334a939b..6bee549d3 100644
--- a/examples/PowerElectronics/ScaleMicrogrid/CMakeLists.txt
+++ b/examples/PowerElectronics/ScaleMicrogrid/CMakeLists.txt
@@ -1,5 +1,14 @@
 add_executable(scalemicrogrid ScaleMicrogrid.cpp)
 add_executable(scalemicrogridarbitrary ScaleMicrogridArbitrary.cpp)
+
+target_include_directories(
+  scalemicrogrid
+  PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/..)
+
+target_include_directories(
+  scalemicrogridarbitrary
+  PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/..)
+
 target_link_libraries(
   scalemicrogrid
   GridKit::power_elec_disgen
diff --git a/examples/PowerElectronics/ScaleMicrogrid/ScaleMicrogrid.cpp b/examples/PowerElectronics/ScaleMicrogrid/ScaleMicrogrid.cpp
index cc8dc28ee..1e97cd56a 100644
--- a/examples/PowerElectronics/ScaleMicrogrid/ScaleMicrogrid.cpp
+++ b/examples/PowerElectronics/ScaleMicrogrid/ScaleMicrogrid.cpp
@@ -1,22 +1,14 @@
 #include <cmath>
-#include <filesystem>
-#include <fstream>
-#include <iomanip>
 #include <iostream>
-#include <memory>
-#include <numbers>
-
-#include <GridKit/Model/PowerElectronics/Bus/MicrogridBus.hpp>
-#include <GridKit/Model/PowerElectronics/Bus/SignalNode.hpp>
-#include <GridKit/Model/PowerElectronics/DistributedGenerator/DistributedGenerator.hpp>
-#include <GridKit/Model/PowerElectronics/MicrogridBusDQ/MicrogridBusDQ.hpp>
-#include <GridKit/Model/PowerElectronics/MicrogridLine/MicrogridLine.hpp>
-#include <GridKit/Model/PowerElectronics/MicrogridLoad/MicrogridLoad.hpp>
+#include <vector>
+
 #include <GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp>
 #include <GridKit/Solver/Dynamic/DynamicSolver.hpp>
 #include <GridKit/Solver/Dynamic/Ida.hpp>
 #include <GridKit/Testing/Testing.hpp>
 
+#include "Common/MicrogridNetwork.hpp"
+
 using index_type = size_t;
 using real_type  = double;
 
@@ -94,143 +86,11 @@ int test(index_type Nsize, real_type error_tol, bool debug_output)
     std::cout << "Using default Nsize = 8.\n";
   }
 
-  // Modeled after the problem in the paper
-  // Every Bus has the same virtual resistance. This is due to the numerical stability as mentioned in the paper.
-  real_type RN = 1.0e4;
-
-  // DG Params Vector
-  // All DGs have the same set of parameters except for the first two.
-  static constexpr auto pi = std::numbers::pi_v<real_type>;
-
-  GridKit::DistributedGeneratorParameters<real_type, index_type> DG_parms1;
-  DG_parms1.wb_  = 2.0 * pi * 50.0;
-  DG_parms1.wc_  = 31.41;
-  DG_parms1.mp_  = 9.4e-5;
-  DG_parms1.Vn_  = 380.0;
-  DG_parms1.nq_  = 1.3e-3;
-  DG_parms1.F_   = 0.75;
-  DG_parms1.Kiv_ = 420.0;
-  DG_parms1.Kpv_ = 0.1;
-  DG_parms1.Kic_ = 2.0e4;
-  DG_parms1.Kpc_ = 15.0;
-  DG_parms1.Cf_  = 5.0e-5;
-  DG_parms1.rLf_ = 0.1;
-  DG_parms1.Lf_  = 1.35e-3;
-  DG_parms1.rLc_ = 0.03;
-  DG_parms1.Lc_  = 0.35e-3;
-
-  GridKit::DistributedGeneratorParameters<real_type, index_type> DG_parms2;
-  DG_parms2.wb_  = 2.0 * pi * 50.0;
-  DG_parms2.wc_  = 31.41;
-  DG_parms2.mp_  = 12.5e-5;
-  DG_parms2.Vn_  = 380.0;
-  DG_parms2.nq_  = 1.5e-3;
-  DG_parms2.F_   = 0.75;
-  DG_parms2.Kiv_ = 390.0;
-  DG_parms2.Kpv_ = 0.05;
-  DG_parms2.Kic_ = 16.0e3;
-  DG_parms2.Kpc_ = 10.5;
-  DG_parms2.Cf_  = 50.0e-6;
-  DG_parms2.rLf_ = 0.1;
-  DG_parms2.Lf_  = 1.35e-3;
-  DG_parms2.rLc_ = 0.03;
-  DG_parms2.Lc_  = 0.35e-3;
-
-  std::vector<GridKit::DistributedGeneratorParameters<real_type, index_type>> DGParams_list(2 * Nsize, DG_parms2);
-
-  DGParams_list[0] = DG_parms1;
-  DGParams_list[1] = DG_parms1;
-
-  // line vector params
-  // Every odd line has the same parameters and every even line has the same parameters
-  real_type              rline1 = 0.23;
-  real_type              Lline1 = 0.1 / (2.0 * pi * 50.0);
-  real_type              rline2 = 0.35;
-  real_type              Lline2 = 0.58 / (2.0 * pi * 50.0);
-  std::vector<real_type> rline_list(2 * Nsize - 1, 0.0);
-  std::vector<real_type> Lline_list(2 * Nsize - 1, 0.0);
-  for (index_type i = 0; i < rline_list.size(); i++)
-  {
-    rline_list[i] = (i % 2) ? rline2 : rline1;
-    Lline_list[i] = (i % 2) ? Lline2 : Lline1;
-  }
-
-  // load parms
-  // Only the first load has the same paramaters.
-  real_type rload1 = 3.0;
-  real_type Lload1 = 2.0 / (2.0 * pi * 50.0);
-  real_type rload2 = 2.0;
-  real_type Lload2 = 1.0 / (2.0 * pi * 50.0);
-
-  std::vector<real_type> rload_list(Nsize, rload2);
-  std::vector<real_type> Lload_list(Nsize, Lload2);
-  rload_list[0] = rload1;
-  Lload_list[0] = Lload1;
-
-  using SignalNode = GridKit::PowerElectronics::SignalNode<double, size_t>;
-  SignalNode dg_signal;
-  sys_model->addNode(&dg_signal);
-
-  using Bus                                     = GridKit::PowerElectronics::MicrogridBus<double, size_t>;
-  std::unique_ptr<std::unique_ptr<Bus>[]> buses = std::make_unique<std::unique_ptr<Bus>[]>(2 * Nsize);
-  for (size_t i = 0; i < 2 * Nsize; i++)
-  {
-    buses[i] = std::make_unique<Bus>();
-    sys_model->addNode(buses[i].get());
-  }
-
-  // Create the reference DG
-  auto* dg_ref = new DistributedGenerator<real_type, index_type>(0,
-                                                                 DGParams_list[0],
-                                                                 true,
-                                                                 &dg_signal,
-                                                                 buses[0].get());
-  sys_model->addComponent(dg_ref);
-
-  // Keep track of models and index location
-  index_type model_id = 1;
-  // Add all other DGs
-  for (index_type i = 1; i < 2 * Nsize; i++)
-  {
-    // current DG to add
-    auto* dg = new DistributedGenerator<real_type, index_type>(model_id++,
-                                                               DGParams_list[i],
-                                                               false,
-                                                               &dg_signal,
-                                                               buses[i].get());
-    sys_model->addComponent(dg);
-  }
-
-  // Load all the Line compoenents
-  for (index_type i = 0; i < 2 * Nsize - 1; i++)
-  {
-    // line
-    auto* line_model = new MicrogridLine<real_type, index_type>(model_id++,
-                                                                rline_list[i],
-                                                                Lline_list[i],
-                                                                &dg_signal,
-                                                                buses[i].get(),
-                                                                buses[i + 1].get());
-    sys_model->addComponent(line_model);
-  }
-
-  //  Load all the Load components
-  for (index_type i = 0; i < Nsize; i++)
-  {
-    auto* load_model = new MicrogridLoad<real_type, index_type>(model_id++,
-                                                                rload_list[i],
-                                                                Lload_list[i],
-                                                                &dg_signal,
-                                                                buses[2 * i].get());
-    sys_model->addComponent(load_model);
-  }
+  // Build and assemble the scaled microgrid network.
+  ScaleMicrogridNetwork network(Nsize);
 
-  // Add all the microgrid Virtual DQ Buses
-  for (index_type i = 0; i < 2 * Nsize; i++)
-  {
-    auto* virDQbus_model = new MicrogridBusDQ<real_type, index_type>(model_id++, RN, buses[i].get());
-    sys_model->addComponent(virDQbus_model);
-  }
+  buildScaleMicrogridNetwork(network);
+  assembleSystem(network, *sys_model);
 
   // allocate all the intial conditions
   sys_model->allocate();
@@ -253,13 +113,15 @@ int test(index_type Nsize, real_type error_tol, bool debug_output)
   // Create initial derivatives specifics generated in MATLAB
   for (index_type i = 0; i < 2 * Nsize; i++)
   {
-    yp[13 * i - 1 + 2] = DGParams_list[i].Vn_;
-    yp[13 * i - 1 + 4] = DGParams_list[i].Kpv_ * DGParams_list[i].Vn_;
-    yp[13 * i - 1 + 6] = (DGParams_list[i].Kpc_ * DGParams_list[i].Kpv_ * DGParams_list[i].Vn_) / DGParams_list[i].Lf_;
+    const auto& params = network.DGParam_list[i];
+
+    yp[13 * i - 1 + 2] = params.Vn_;
+    yp[13 * i - 1 + 4] = params.Kpv_ * params.Vn_;
+    yp[13 * i - 1 + 6] = (params.Kpc_ * params.Kpv_ * params.Vn_) / params.Lf_;
   }
 
   // since the intial P_com = 0, the set the intial vector to the reference frame
-  y[dg_signal.getNodeConnection(0).idx_] = DG_parms1.wb_;
+  y[network.dg_signal.getNodeConnection(0).idx_] = network.DGParam_list[0].wb_;
 
   sys_model->y().setDataUpdated();
   sys_model->yp().setDataUpdated();
diff --git a/examples/PowerElectronics/ScaleMicrogrid/ScaleMicrogridArbitrary.cpp b/examples/PowerElectronics/ScaleMicrogrid/ScaleMicrogridArbitrary.cpp
index a8ab6a4a1..c83ce2979 100644
--- a/examples/PowerElectronics/ScaleMicrogrid/ScaleMicrogridArbitrary.cpp
+++ b/examples/PowerElectronics/ScaleMicrogrid/ScaleMicrogridArbitrary.cpp
@@ -1,20 +1,13 @@
 #include <cmath>
-#include <filesystem>
-#include <fstream>
-#include <iomanip>
 #include <iostream>
 
-#include <GridKit/Model/PowerElectronics/Bus/MicrogridBus.hpp>
-#include <GridKit/Model/PowerElectronics/Bus/SignalNode.hpp>
-#include <GridKit/Model/PowerElectronics/DistributedGenerator/DistributedGenerator.hpp>
-#include <GridKit/Model/PowerElectronics/MicrogridBusDQ/MicrogridBusDQ.hpp>
-#include <GridKit/Model/PowerElectronics/MicrogridLine/MicrogridLine.hpp>
-#include <GridKit/Model/PowerElectronics/MicrogridLoad/MicrogridLoad.hpp>
 #include <GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp>
 #include <GridKit/Solver/Dynamic/DynamicSolver.hpp>
 #include <GridKit/Solver/Dynamic/Ida.hpp>
 #include <GridKit/Testing/Testing.hpp>
 
+#include "Common/MicrogridNetwork.hpp"
+
 using index_type = size_t;
 using real_type  = double;
 
@@ -79,149 +72,11 @@ int printMicrogridSystems(index_type N_size)
     return 1;
   }
 
-  // Modeled after the problem in the paper
-  // Every Bus has the same virtual resistance. This is due to numerical stability as mentioned in the paper.
-  real_type RN = 1.0e4;
-
-  // DG Params Vector
-  // All DGs have the same set of parameters except for the first two.
-  static constexpr auto pi = std::numbers::pi_v<real_type>;
-
-  GridKit::DistributedGeneratorParameters<real_type, index_type> DG_parms1;
-  DG_parms1.wb_  = 2.0 * pi * 50.0;
-  DG_parms1.wc_  = 31.41;
-  DG_parms1.mp_  = 9.4e-5;
-  DG_parms1.Vn_  = 380.0;
-  DG_parms1.nq_  = 1.3e-3;
-  DG_parms1.F_   = 0.75;
-  DG_parms1.Kiv_ = 420.0;
-  DG_parms1.Kpv_ = 0.1;
-  DG_parms1.Kic_ = 2.0e4;
-  DG_parms1.Kpc_ = 15.0;
-  DG_parms1.Cf_  = 5.0e-5;
-  DG_parms1.rLf_ = 0.1;
-  DG_parms1.Lf_  = 1.35e-3;
-  DG_parms1.rLc_ = 0.03;
-  DG_parms1.Lc_  = 0.35e-3;
-
-  GridKit::DistributedGeneratorParameters<real_type, index_type> DG_parms2;
-  DG_parms2.wb_  = 2.0 * pi * 50.0;
-  DG_parms2.wc_  = 31.41;
-  DG_parms2.mp_  = 12.5e-5;
-  DG_parms2.Vn_  = 380.0;
-  DG_parms2.nq_  = 1.5e-3;
-  DG_parms2.F_   = 0.75;
-  DG_parms2.Kiv_ = 390.0;
-  DG_parms2.Kpv_ = 0.05;
-  DG_parms2.Kic_ = 16.0e3;
-  DG_parms2.Kpc_ = 10.5;
-  DG_parms2.Cf_  = 50.0e-6;
-  DG_parms2.rLf_ = 0.1;
-  DG_parms2.Lf_  = 1.35e-3;
-  DG_parms2.rLc_ = 0.03;
-  DG_parms2.Lc_  = 0.35e-3;
-
-  std::vector<GridKit::DistributedGeneratorParameters<real_type, index_type>> DGParams_list(2 * N_size, DG_parms2);
-
-  // First two generators use parameters 1
-  if (DGParams_list.size() >= 1)
-    DGParams_list[0] = DG_parms1;
-  if (DGParams_list.size() >= 2)
-    DGParams_list[1] = DG_parms1;
-
-  // line vector params
-  // Every odd line has the same parameters and every even line has the same parameters
-  real_type              rline1 = 0.23;
-  real_type              Lline1 = 0.1 / (2.0 * pi * 50.0);
-  real_type              rline2 = 0.35;
-  real_type              Lline2 = 0.58 / (2.0 * pi * 50.0);
-  std::vector<real_type> rline_list(2 * N_size - 1, 0.0);
-  std::vector<real_type> Lline_list(2 * N_size - 1, 0.0);
-  for (index_type i = 0; i < rline_list.size(); i++)
-  {
-    rline_list[i] = (i % 2) ? rline2 : rline1;
-    Lline_list[i] = (i % 2) ? Lline2 : Lline1;
-  }
-
-  // load parms
-  // Only the first load has the same paramaters.
-  real_type rload1 = 3.0;
-  real_type Lload1 = 2.0 / (2.0 * pi * 50.0);
-  real_type rload2 = 2.0;
-  real_type Lload2 = 1.0 / (2.0 * pi * 50.0);
-
-  std::vector<real_type> rload_list(N_size, rload2);
-  std::vector<real_type> Lload_list(N_size, Lload2);
-  if (rload_list.size() >= 1)
-  {
-    rload_list[0] = rload1;
-    Lload_list[0] = Lload1;
-  }
-
-  using SignalNode = GridKit::PowerElectronics::SignalNode<double, size_t>;
-  SignalNode dg_signal;
-  sys_model.addNode(&dg_signal);
+  // Build and assemble the scaled microgrid network.
+  ScaleMicrogridNetwork network(N_size);
 
-  using Bus                                     = GridKit::PowerElectronics::MicrogridBus<double, size_t>;
-  std::unique_ptr<std::unique_ptr<Bus>[]> buses = std::make_unique<std::unique_ptr<Bus>[]>(2 * N_size);
-  for (size_t i = 0; i < 2 * N_size; i++)
-  {
-    buses[i] = std::make_unique<Bus>();
-    sys_model.addNode(buses[i].get());
-  }
-
-  // Create the reference DG
-  auto* dg_ref = new DistributedGenerator<real_type, index_type>(0,
-                                                                 DGParams_list[0],
-                                                                 true,
-                                                                 &dg_signal,
-                                                                 buses[0].get());
-  sys_model.addComponent(dg_ref);
-
-  // Keep track of models and index location
-  index_type model_id = 1;
-  // Add all other DGs
-  for (index_type i = 1; i < 2 * N_size; i++)
-  {
-    // current DG to add
-    auto* dg = new DistributedGenerator<real_type, index_type>(model_id++,
-                                                               DGParams_list[i],
-                                                               false,
-                                                               &dg_signal,
-                                                               buses[i].get());
-    sys_model.addComponent(dg);
-  }
-
-  // Load all the Line components
-  for (index_type i = 0; i < 2 * N_size - 1; i++)
-  {
-    // line
-    auto* line_model = new MicrogridLine<real_type, index_type>(model_id++,
-                                                                rline_list[i],
-                                                                Lline_list[i],
-                                                                &dg_signal,
-                                                                buses[i].get(),
-                                                                buses[i + 1].get());
-    sys_model.addComponent(line_model);
-  }
-
-  //  Load all the Load components
-  for (index_type i = 0; i < N_size; i++)
-  {
-    auto* load_model = new MicrogridLoad<real_type, index_type>(model_id++,
-                                                                rload_list[i],
-                                                                Lload_list[i],
-                                                                &dg_signal,
-                                                                buses[2 * i].get());
-    sys_model.addComponent(load_model);
-  }
-
-  // Add all the microgrid Virtual DQ Buses
-  for (index_type i = 0; i < 2 * N_size; i++)
-  {
-    auto* virDQbus_model = new MicrogridBusDQ<real_type, index_type>(model_id++, RN, buses[i].get());
-    sys_model.addComponent(virDQbus_model);
-  }
+  buildScaleMicrogridNetwork(network);
+  assembleSystem(network, sys_model);
 
   // allocate all the initial conditions
   sys_model.allocate();
@@ -236,16 +91,22 @@ int printMicrogridSystems(index_type N_size)
     yp[i] = 0.0;
   }
 
+  //-------------------------------------------------------------------
   // Create Initial derivatives specifics generated in MATLAB
+  //-------------------------------------------------------------------
   for (index_type i = 0; i < 2 * N_size; i++)
   {
-    yp[13 * i - 1 + 3] = DGParams_list[i].Vn_;
-    yp[13 * i - 1 + 5] = DGParams_list[i].Kpv_ * DGParams_list[i].Vn_;
-    yp[13 * i - 1 + 7] = (DGParams_list[i].Kpc_ * DGParams_list[i].Kpv_ * DGParams_list[i].Vn_) / DGParams_list[i].Lf_;
+    const auto& params = network.DGParam_list[i];
+
+    yp[13 * i - 1 + 3] = params.Vn_;
+    yp[13 * i - 1 + 5] = params.Kpv_ * params.Vn_;
+    yp[13 * i - 1 + 7] = (params.Kpc_ * params.Kpv_ * params.Vn_) / params.Lf_;
   }
 
+  //---------------------------------------------------------------------------
   // since the initial P_com = 0, set the initial vector to the reference frame
-  y[dg_signal.getNodeConnection(0).idx_] = DG_parms1.wb_;
+  //---------------------------------------------------------------------------
+  y[network.dg_signal.getNodeConnection(0).idx_] = network.DGParam_list[0].wb_;
 
   sys_model.y().setDataUpdated();
   sys_model.yp().setDataUpdated();
diff --git a/tests/UnitTests/PowerElectronics/CMakeLists.txt b/tests/UnitTests/PowerElectronics/CMakeLists.txt
index 74d50eb1e..88a0bc7f6 100644
--- a/tests/UnitTests/PowerElectronics/CMakeLists.txt
+++ b/tests/UnitTests/PowerElectronics/CMakeLists.txt
@@ -3,6 +3,26 @@ target_link_libraries(
   test_power_electronics_node
   PRIVATE GridKit::power_electronics_circuit_node GridKit::testing)
 
+add_executable(test_subsystem_model_with_hires runSubsystemModelWithHiresTest.cpp)
+target_link_libraries(
+  test_subsystem_model_with_hires
+  PRIVATE GridKit::power_elec_partition_interfaces 
+          GridKit::testing
+          GridKit::sparse_matrix)
+
+add_executable(test_power_electronics_component_clone runComponentCloneTests.cpp)
+target_link_libraries(
+  test_power_electronics_component_clone
+  PRIVATE GridKit::power_elec_disgen
+          GridKit::power_elec_microline
+          GridKit::power_elec_microload
+          GridKit::power_elec_microbusdq
+          GridKit::testing)
+
 add_test(NAME PowerElectronicsNodeTest COMMAND $<TARGET_FILE:test_power_electronics_node>)
+add_test(NAME SubsystemModelWithHires COMMAND $<TARGET_FILE:test_subsystem_model_with_hires>)
+add_test(NAME PowerElectronicsComponentCloneTest COMMAND $<TARGET_FILE:test_power_electronics_component_clone>)
 
 install(TARGETS test_power_electronics_node RUNTIME DESTINATION bin)
+install(TARGETS test_subsystem_model_with_hires RUNTIME DESTINATION bin)
+install(TARGETS test_power_electronics_component_clone RUNTIME DESTINATION bin)
diff --git a/tests/UnitTests/PowerElectronics/ComponentCloneTests.hpp b/tests/UnitTests/PowerElectronics/ComponentCloneTests.hpp
new file mode 100644
index 000000000..5124442d8
--- /dev/null
+++ b/tests/UnitTests/PowerElectronics/ComponentCloneTests.hpp
@@ -0,0 +1,263 @@
+
+#include <GridKit/Model/PowerElectronics/Bus/MicrogridBus.hpp>
+#include <GridKit/Model/PowerElectronics/Bus/SignalNode.hpp>
+#include <GridKit/Model/PowerElectronics/DistributedGenerator/DistributedGenerator.hpp>
+#include <GridKit/Model/PowerElectronics/MicrogridBusDQ/MicrogridBusDQ.hpp>
+#include <GridKit/Model/PowerElectronics/MicrogridLine/MicrogridLine.hpp>
+#include <GridKit/Model/PowerElectronics/MicrogridLoad/MicrogridLoad.hpp>
+#include <GridKit/Testing/TestHelpers.hpp>
+#include <GridKit/Testing/Testing.hpp>
+
+namespace GridKit
+{
+  template <typename ComponentT, typename ScalarT, typename IdxT>
+  bool verifyComponentClone(ComponentT& component)
+  {
+    using RealT = typename CircuitComponent<ScalarT, IdxT>::RealT;
+
+    bool success = true;
+
+    auto* clone = dynamic_cast<ComponentT*>(component.clone());
+
+    if (clone == nullptr)
+    {
+      return false;
+    }
+
+    /**************************************************************************
+     * Verify the Clone Initially Matches the Original
+     **************************************************************************/
+
+    success &= clone != &component;
+    success &= clone->size() == component.size();
+    success &= clone->nnz() == component.nnz();
+    success &= clone->getInternalSize() == component.getInternalSize();
+    success &= clone->getExternSize() == component.getExternSize();
+    success &= clone->getExternIndices() == component.getExternIndices();
+    success &= clone->getIDcomponent() == component.getIDcomponent();
+
+    /**************************************************************************
+     * Verify Connection Independence
+     **************************************************************************/
+
+    for (IdxT i = 0; i < component.size(); ++i)
+    {
+      const IdxT connection = component.getNodeConnection(i);
+
+      success &= clone->getNodeConnection(i) == connection;
+
+      clone->setConnectionNodes(i, connection + 1);
+
+      success &= component.getNodeConnection(i) == connection;
+      success &= clone->getNodeConnection(i) == connection + 1;
+
+      clone->setConnectionNodes(i, connection);
+    }
+
+    /**************************************************************************
+     * Verify State Independence
+     **************************************************************************/
+
+    auto checkVectorIndependence = [&success](auto& original, auto& copy)
+    {
+      success &= original.getSize() == copy.getSize();
+
+      if (original.getSize() == 0)
+      {
+        return;
+      }
+
+      auto* original_data = original.getData();
+      auto* copy_data     = copy.getData();
+
+      success &= original_data != copy_data;
+
+      const auto original_value = original_data[0];
+
+      copy_data[0] = original_value + 1.0;
+
+      success &= original_data[0] == original_value;
+      success &= copy_data[0] == original_value + 1.0;
+
+      copy_data[0] = original_value;
+    };
+
+    checkVectorIndependence(component.y(), clone->y());
+    checkVectorIndependence(component.yp(), clone->yp());
+    checkVectorIndependence(component.getResidual(), clone->getResidual());
+    checkVectorIndependence(component.absoluteTolerance(), clone->absoluteTolerance());
+    checkVectorIndependence(component.param(), clone->param());
+
+    /**************************************************************************
+     * Verify Jacobian Independence
+     **************************************************************************/
+
+    if (component.nnz() > 0)
+    {
+      auto* original_rows   = component.jacobianCooRows();
+      auto* original_cols   = component.jacobianCooCols();
+      auto* original_values = component.jacobianCooValues();
+
+      auto* clone_rows   = clone->jacobianCooRows();
+      auto* clone_cols   = clone->jacobianCooCols();
+      auto* clone_values = clone->jacobianCooValues();
+
+      success &= clone_rows != original_rows;
+      success &= clone_cols != original_cols;
+      success &= clone_values != original_values;
+
+      for (IdxT i = 0; i < component.nnz(); ++i)
+      {
+        success &= clone_rows[i] == original_rows[i];
+        success &= clone_cols[i] == original_cols[i];
+        success &= clone_values[i] == original_values[i];
+      }
+
+      const IdxT  original_row   = original_rows[0];
+      const IdxT  original_col   = original_cols[0];
+      const RealT original_value = original_values[0];
+
+      clone_rows[0]   = original_row + 1;
+      clone_cols[0]   = original_col + 1;
+      clone_values[0] = original_value + 1.0;
+
+      success &= original_rows[0] == original_row;
+      success &= original_cols[0] == original_col;
+      success &= original_values[0] == original_value;
+
+      clone_rows[0]   = original_row;
+      clone_cols[0]   = original_col;
+      clone_values[0] = original_value;
+    }
+
+    delete clone;
+
+    return success;
+  }
+
+  namespace Testing
+  {
+    template <typename ScalarT, typename IdxT>
+    class CircuitComponentCloneTests
+    {
+      using SignalNode          = PowerElectronics::SignalNode<ScalarT, IdxT>;
+      using Bus                 = PowerElectronics::MicrogridBus<ScalarT, IdxT>;
+      using BusDQ               = MicrogridBusDQ<ScalarT, IdxT>;
+      using Generator           = DistributedGenerator<ScalarT, IdxT>;
+      using GeneratorParameters = DistributedGeneratorParameters<ScalarT, IdxT>;
+      using Line                = MicrogridLine<ScalarT, IdxT>;
+      using Load                = MicrogridLoad<ScalarT, IdxT>;
+
+    public:
+      CircuitComponentCloneTests()
+      {
+        /**************************************************************************
+         * Construct Network Nodes
+         **************************************************************************/
+
+        signal_.allocate();
+
+        bus1_.allocate();
+        bus2_.allocate();
+
+        /**************************************************************************
+         * Distributed Generator Parameters
+         **************************************************************************/
+
+        generator_parameters_.wb_  = 2.0 * M_PI * 50.0;
+        generator_parameters_.wc_  = 31.41;
+        generator_parameters_.mp_  = 9.4e-5;
+        generator_parameters_.Vn_  = 380.0;
+        generator_parameters_.nq_  = 1.3e-3;
+        generator_parameters_.F_   = 0.75;
+        generator_parameters_.Kiv_ = 420.0;
+        generator_parameters_.Kpv_ = 0.1;
+        generator_parameters_.Kic_ = 2.0e4;
+        generator_parameters_.Kpc_ = 15.0;
+        generator_parameters_.Cf_  = 5.0e-5;
+        generator_parameters_.rLf_ = 0.1;
+        generator_parameters_.Lf_  = 1.35e-3;
+        generator_parameters_.rLc_ = 0.03;
+        generator_parameters_.Lc_  = 0.35e-3;
+
+        /**************************************************************************
+         * Construct Components
+         **************************************************************************/
+
+        generator_ = new Generator(1, generator_parameters_, true, &signal_, &bus1_);
+
+        line_ = new Line(2, 0.23, 0.1 / (2.0 * M_PI * 50.0), &signal_, &bus1_, &bus2_);
+
+        load_ = new Load(3, 3.0, 2.0 / (2.0 * M_PI * 50.0), &signal_, &bus1_);
+
+        bus_dq_ = new BusDQ(4, 1.0e4, &bus1_);
+
+        /**************************************************************************
+         * Allocate Components
+         **************************************************************************/
+
+        generator_->allocate();
+        line_->allocate();
+        load_->allocate();
+        bus_dq_->allocate();
+      }
+
+      ~CircuitComponentCloneTests()
+      {
+        delete generator_;
+        delete line_;
+        delete load_;
+        delete bus_dq_;
+      }
+
+      TestOutcome distributedGeneratorClone()
+      {
+        TestStatus success = true;
+
+        success *= verifyComponentClone<Generator, ScalarT, IdxT>(*generator_);
+
+        return success.report(__func__);
+      }
+
+      TestOutcome microgridLineClone()
+      {
+        TestStatus success = true;
+
+        success *= verifyComponentClone<Line, ScalarT, IdxT>(*line_);
+
+        return success.report(__func__);
+      }
+
+      TestOutcome microgridLoadClone()
+      {
+        TestStatus success = true;
+
+        success *= verifyComponentClone<Load, ScalarT, IdxT>(*load_);
+
+        return success.report(__func__);
+      }
+
+      TestOutcome microgridBusDQClone()
+      {
+        TestStatus success = true;
+
+        success *= verifyComponentClone<BusDQ, ScalarT, IdxT>(*bus_dq_);
+
+        return success.report(__func__);
+      }
+
+    private:
+      SignalNode signal_;
+
+      Bus bus1_;
+      Bus bus2_;
+
+      GeneratorParameters generator_parameters_;
+
+      Generator* generator_{nullptr};
+      Line*      line_{nullptr};
+      Load*      load_{nullptr};
+      BusDQ*     bus_dq_{nullptr};
+    };
+  } // namespace Testing
+} // namespace GridKit
diff --git a/tests/UnitTests/PowerElectronics/README.md b/tests/UnitTests/PowerElectronics/README.md
new file mode 100644
index 000000000..be5c30248
--- /dev/null
+++ b/tests/UnitTests/PowerElectronics/README.md
@@ -0,0 +1,139 @@
+## HIRES Partitioning Test Problem
+
+The HIRES test problem is a simple ODE system with eight variables. To
+demonstrate the partitioning machinery in GridKit, the system is divided into
+three components. Equations 1--3 belong to **Component 1**, equations 4--5
+belong to **Component 2**, which is modeled as a bus, and equations 6--8 belong
+to **Component 3**.
+
+> **Note:** HIRES is not a circuit problem. In this example, it is modeled to
+> resemble GridKit circuit components so that the existing partitioning
+> machinery can be used. The example also provides a simple test problem for
+> verifying the order of co-simulation methods.
+
+The full HIRES system is
+
+$$
+\begin{aligned}
+f_1 &= \frac{dy_1}{dt} +1.71y_1 -0.43y_2 -8.32y_3 -0.0007, \\
+f_2 &= \frac{dy_2}{dt} -1.71y_1 +8.75y_2, \\
+f_3 &= \frac{dy_3}{dt} +10.03y_3 -0.43y_4 -0.035y_5, \\
+f_4 &= \frac{dy_4}{dt} -8.32y_2 -1.71y_3 +1.12y_4, \\
+f_5 &= \frac{dy_5}{dt} +1.745y_5 -0.43y_6 -0.43y_7, \\
+f_6 &= \frac{dy_6}{dt} +280y_6y_8 -0.69y_4 -1.71y_5
+       +0.43y_6 -0.69y_7, \\
+f_7 &= \frac{dy_7}{dt} -280y_6y_8 +1.81y_7, \\
+f_8 &= \frac{dy_8}{dt} +280y_6y_8 -1.81y_7.
+\end{aligned}
+$$
+
+For the component representation, equations $f_4$ and $f_5$ are decomposed
+to expose the contributions from each component:
+
+$$
+\begin{aligned}
+f_4
+&= \frac{dy_4}{dt} -8.32y_2 -1.71y_3 +1.12y_4 \qquad
+&\longrightarrow
+\left(\frac{dy_4}{dt}+y_4\right)
++\left(0.1y_4-8.32y_2-1.71y_3\right)
++0.02y_4,
+\\[6pt]
+f_5
+&= \frac{dy_5}{dt}+1.745y_5-0.43y_6-0.43y_7 \qquad
+&\longrightarrow
+\left(\frac{dy_5}{dt}+y_5\right)
++0.7y_5
++\left[0.045y_5-0.43y_6-0.43y_7\right].
+\end{aligned}
+$$
+
+The decomposition does not change the HIRES system. It separates the terms in
+the bus equations (conviniently choosen to be equation 4 and 5) according to the component responsible for each contribution.
+
+### Component 1
+
+Component 1 has three internal equations:
+
+$$
+\begin{aligned}
+f_1 &= \frac{dy_1}{dt} +1.71y_1 -0.43y_2 -8.32y_3 -0.0007, \\
+f_2 &= \frac{dy_2}{dt} -1.71y_1 +8.75y_2, \\
+f_3 &= \frac{dy_3}{dt} +10.03y_3 -0.43y_4 -0.035y_5.
+\end{aligned}
+$$
+
+It also contributes the following terms to the bus equations as its external contribution:
+
+$$
+\begin{aligned}
+f_4^{(1)} &= 0.1y_4 -8.32y_2 -1.71y_3, \\
+f_5^{(1)} &= 0.7y_5.
+\end{aligned}
+$$
+
+### Component 2 (HiresBus)
+
+Component 2 represents the bus and owns the following contributions to
+equations 4 and 5:
+
+$$
+\begin{aligned}
+f_4^{(2)} &= \frac{dy_4}{dt} + y_4, \\
+f_5^{(2)} &= \frac{dy_5}{dt} + y_5.
+\end{aligned}
+$$
+
+The remaining terms in these equations are supplied by the components
+connected to the bus.
+
+### Component 3
+
+Component 3 has three internal equations:
+
+$$
+\begin{aligned}
+f_6 &= \frac{dy_6}{dt} -280y_6y_8 +0.69y_4 +1.71y_5
+       -0.43y_6 +0.69y_7, \\
+f_7 &= \frac{dy_7}{dt} +280y_6y_8 -1.81y_7, \\
+f_8 &= \frac{dy_8}{dt} -280y_6y_8 +1.81y_7.
+\end{aligned}
+$$
+
+It also contributes the following terms to the bus equations as its external contribution:
+
+$$
+\begin{aligned}
+f_4^{(3)} &= 0.02y_4, \\
+f_5^{(3)} &= 0.045y_5 -0.43y_6 -0.43y_7.
+\end{aligned}
+$$
+
+### Bus Residual Assembly
+
+The complete bus residuals are obtained by adding the contributions from
+Components 1, 2, and 3:
+
+$$
+f_4 = f_4^{(1)} + f_4^{(2)} + f_4^{(3)},
+$$
+
+$$
+f_5 = f_5^{(1)} + f_5^{(2)} + f_5^{(3)}.
+$$
+
+This reconstruction gives exactly the corresponding equations in the original
+HIRES system.
+
+### Partitioning
+
+The full HIRES system in component form looks like this:
+
+```text
+Component 1  --------  Component 2 (HiresBus)  --------  Component 3
+```
+
+The system is divided between **Component 2 (HiresBus)** and **Component 3**, and
+a `BusPartitionInterface` is added to the first partition. The resulting
+partition residuals are then evaluated independently and are then compared with the full-system residual to verify
+that the partitioned evaluation reproduces the original system.
diff --git a/tests/UnitTests/PowerElectronics/SubsystemModelWithHiresTest.hpp b/tests/UnitTests/PowerElectronics/SubsystemModelWithHiresTest.hpp
new file mode 100644
index 000000000..fa524bfc1
--- /dev/null
+++ b/tests/UnitTests/PowerElectronics/SubsystemModelWithHiresTest.hpp
@@ -0,0 +1,749 @@
+#include <algorithm>
+#include <array>
+#include <cmath>
+#include <cstddef>
+#include <limits>
+#include <vector>
+
+#include <GridKit/Model/PowerElectronics/Bus/MicrogridBus.hpp>
+#include <GridKit/Model/PowerElectronics/CircuitComponent.hpp>
+#include <GridKit/Model/PowerElectronics/NodeBase.hpp>
+#include <GridKit/Model/PowerElectronics/PartitionInterface/BusPartitionInterface.hpp>
+#include <GridKit/Model/PowerElectronics/SubsystemModel.hpp>
+#include <GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp>
+#include <GridKit/Testing/TestHelpers.hpp>
+#include <GridKit/Testing/Testing.hpp>
+
+namespace GridKit
+{
+
+  /*!
+   * @brief Hires Component 1 class.
+   *
+   */
+  template <class ScalarT, typename IdxT>
+  class HiresComponent1 : public CircuitComponent<ScalarT, IdxT>
+  {
+    using RealT = typename CircuitComponent<ScalarT, IdxT>::RealT;
+    using NodeT = typename PowerElectronics::NodeBase<ScalarT, IdxT>;
+
+    using CircuitComponent<ScalarT, IdxT>::size_;
+    using CircuitComponent<ScalarT, IdxT>::nnz_;
+    using CircuitComponent<ScalarT, IdxT>::alpha_;
+    using CircuitComponent<ScalarT, IdxT>::y_ext_;
+    using CircuitComponent<ScalarT, IdxT>::y_int_;
+    using CircuitComponent<ScalarT, IdxT>::yp_ext_;
+    using CircuitComponent<ScalarT, IdxT>::yp_int_;
+    using CircuitComponent<ScalarT, IdxT>::abs_tol_;
+    using CircuitComponent<ScalarT, IdxT>::f_ext_;
+    using CircuitComponent<ScalarT, IdxT>::f_int_;
+    using CircuitComponent<ScalarT, IdxT>::idc_;
+
+    using CircuitComponent<ScalarT, IdxT>::extern_indices_;
+    using CircuitComponent<ScalarT, IdxT>::n_extern_;
+    using CircuitComponent<ScalarT, IdxT>::n_intern_;
+
+  public:
+    HiresComponent1(NodeT* bus, IdxT id)
+      : node_ref_(bus)
+    {
+      size_           = 5;
+      n_intern_       = 3;
+      n_extern_       = 2;
+      extern_indices_ = {0, 1};
+      idc_            = id;
+      nnz_            = 12;
+    }
+
+    ~HiresComponent1()
+    {
+    }
+
+    int allocate() final
+    {
+      CircuitComponent<ScalarT, IdxT>::allocate();
+
+      this->setExternalConnectionNodes(0, node_ref_->getNodeConnection(0));
+      this->setExternalConnectionNodes(1, node_ref_->getNodeConnection(1));
+
+      return 0;
+    }
+
+    int initialize() final
+    {
+      return 0;
+    }
+
+    int tagDifferentiable() final
+    {
+      return 0;
+    }
+
+    int evaluateInternalResidual() final
+    {
+      // Internals
+      f_int_[0] = -yp_int_[0] - 1.71 * y_int_[0] + 0.43 * y_int_[1] + 8.32 * y_int_[2] + 0.0007;
+      f_int_[1] = -yp_int_[1] + 1.71 * y_int_[0] - 8.75 * y_int_[1];
+      f_int_[2] = -yp_int_[2] - 10.03 * y_int_[2] + 0.43 * *y_ext_[0] + 0.035 * *y_ext_[1];
+
+      return 0;
+    }
+
+    int evaluateExternalResidual()
+    {
+      // outputs
+      *f_ext_[0] += 8.32 * y_int_[1] + 1.71 * y_int_[2] - 0.1 * *y_ext_[0];
+      *f_ext_[1] += -0.7 * *y_ext_[1];
+
+      return 0;
+    }
+
+    int evaluateJacobian() final
+    {
+
+      this->zeroJacMatrix();
+
+      // Internal Jacobian Entries
+      std::vector<IdxT>  row = {2, 2, 2, 3, 3, 4, 4, 4};
+      std::vector<IdxT>  col = {2, 3, 4, 2, 3, 4, 0, 1};
+      std::vector<RealT> val = {-1.71 - alpha_, 0.43, 8.32, 1.71, -8.75 - alpha_, -10.03 - alpha_, 0.43, 0.035};
+
+      this->setJacValues(row, col, val);
+
+      // External Jacobian Entries
+      row = {0, 0, 0, 1};
+      col = {3, 4, 0, 1};
+      val = {8.32, 1.71, -0.1, -0.7};
+
+      this->setJacValues(row, col, val);
+
+      return 0;
+    }
+
+    int evaluateIntegrand() final
+    {
+      return 0;
+    }
+
+    int initializeAdjoint() final
+    {
+      return 0;
+    }
+
+    int evaluateAdjointResidual() final
+    {
+      return 0;
+    }
+
+    int evaluateAdjointIntegrand() final
+    {
+      return 0;
+    }
+
+    /**
+     * @brief Compute the absolute tolerance for each variable in the model
+     */
+    int setAbsoluteTolerance(RealT rel_tol) final
+    {
+      abs_tol_.setToConst(static_cast<ScalarT>(rel_tol));
+      return 0;
+    }
+
+    CircuitComponent<ScalarT, IdxT>* clone() const
+    {
+      return new HiresComponent1<ScalarT, IdxT>(*this);
+    }
+
+  private:
+    NodeT* node_ref_;
+  };
+
+  /*!
+   * @brief Hires Bus Component (Component 2).
+   *
+   */
+  template <class ScalarT, typename IdxT>
+  class HiresBus : public CircuitComponent<ScalarT, IdxT>
+  {
+
+    using RealT = typename CircuitComponent<ScalarT, IdxT>::RealT;
+    using NodeT = typename PowerElectronics::NodeBase<ScalarT, IdxT>;
+
+    using CircuitComponent<ScalarT, IdxT>::size_;
+    using CircuitComponent<ScalarT, IdxT>::nnz_;
+    using CircuitComponent<ScalarT, IdxT>::alpha_;
+    using CircuitComponent<ScalarT, IdxT>::y_ext_;
+    using CircuitComponent<ScalarT, IdxT>::y_int_;
+    using CircuitComponent<ScalarT, IdxT>::yp_ext_;
+    using CircuitComponent<ScalarT, IdxT>::yp_int_;
+    using CircuitComponent<ScalarT, IdxT>::abs_tol_;
+    using CircuitComponent<ScalarT, IdxT>::f_ext_;
+    using CircuitComponent<ScalarT, IdxT>::f_int_;
+    using CircuitComponent<ScalarT, IdxT>::idc_;
+
+    using CircuitComponent<ScalarT, IdxT>::extern_indices_;
+    using CircuitComponent<ScalarT, IdxT>::n_extern_;
+    using CircuitComponent<ScalarT, IdxT>::n_intern_;
+
+  public:
+    HiresBus(NodeT* bus, IdxT id)
+      : node_ref_(bus)
+    {
+      size_           = 2;
+      n_intern_       = 0;
+      n_extern_       = 2;
+      extern_indices_ = {0, 1};
+      idc_            = id;
+      nnz_            = 2;
+    }
+
+    ~HiresBus()
+    {
+    }
+
+    int allocate() final
+    {
+      CircuitComponent<ScalarT, IdxT>::allocate();
+
+      this->setExternalConnectionNodes(0, node_ref_->getNodeConnection(0));
+      this->setExternalConnectionNodes(1, node_ref_->getNodeConnection(1));
+
+      return 0;
+    }
+
+    int initialize() final
+    {
+      return 0;
+    }
+
+    int tagDifferentiable() final
+    {
+      return 0;
+    }
+
+    int evaluateInternalResidual() final
+    {
+      return 0;
+    }
+
+    int evaluateExternalResidual() final
+    {
+      *f_ext_[0] += -*yp_ext_[0] - *y_ext_[0];
+      *f_ext_[1] += -*yp_ext_[1] - *y_ext_[1];
+
+      return 0;
+    }
+
+    int evaluateJacobian() final
+    {
+      this->zeroJacMatrix();
+
+      std::vector<IdxT>  row = {0, 1};
+      std::vector<IdxT>  col = {0, 1};
+      std::vector<RealT> val = {-alpha_ - 1.0, -alpha_ - 1.0};
+
+      this->setJacValues(row, col, val);
+
+      return 0;
+    }
+
+    int evaluateIntegrand() final
+    {
+      return 0;
+    }
+
+    int initializeAdjoint() final
+    {
+      return 0;
+    }
+
+    int evaluateAdjointResidual() final
+    {
+      return 0;
+    }
+
+    int evaluateAdjointIntegrand() final
+    {
+      return 0;
+    }
+
+    /**
+     * @brief Compute the absolute tolerance for each variable in the model
+     */
+    int setAbsoluteTolerance(RealT rel_tol) final
+    {
+      abs_tol_.setToConst(static_cast<ScalarT>(rel_tol));
+      return 0;
+    }
+
+    CircuitComponent<ScalarT, IdxT>* clone() const
+    {
+      return new HiresBus<ScalarT, IdxT>(*this);
+    }
+
+  private:
+    NodeT* node_ref_;
+  };
+
+  /*!
+   * @brief Hires Component 3 class.
+   *
+   */
+  template <class ScalarT, typename IdxT>
+  class HiresComponent3 : public CircuitComponent<ScalarT, IdxT>
+  {
+    using RealT = typename CircuitComponent<ScalarT, IdxT>::RealT;
+    using NodeT = typename PowerElectronics::NodeBase<ScalarT, IdxT>;
+
+    using CircuitComponent<ScalarT, IdxT>::size_;
+    using CircuitComponent<ScalarT, IdxT>::nnz_;
+    using CircuitComponent<ScalarT, IdxT>::alpha_;
+    using CircuitComponent<ScalarT, IdxT>::y_ext_;
+    using CircuitComponent<ScalarT, IdxT>::y_int_;
+    using CircuitComponent<ScalarT, IdxT>::yp_ext_;
+    using CircuitComponent<ScalarT, IdxT>::yp_int_;
+    using CircuitComponent<ScalarT, IdxT>::abs_tol_;
+    using CircuitComponent<ScalarT, IdxT>::f_ext_;
+    using CircuitComponent<ScalarT, IdxT>::f_int_;
+    using CircuitComponent<ScalarT, IdxT>::idc_;
+
+    using CircuitComponent<ScalarT, IdxT>::extern_indices_;
+    using CircuitComponent<ScalarT, IdxT>::n_extern_;
+    using CircuitComponent<ScalarT, IdxT>::n_intern_;
+
+  public:
+    HiresComponent3(NodeT* bus, IdxT id)
+      : node_ref_(bus)
+    {
+      size_           = 5;
+      n_intern_       = 3;
+      n_extern_       = 2;
+      extern_indices_ = {0, 1};
+      idc_            = id;
+      nnz_            = 15;
+    }
+
+    ~HiresComponent3()
+    {
+    }
+
+    int allocate()
+    {
+      CircuitComponent<ScalarT, IdxT>::allocate();
+
+      this->setExternalConnectionNodes(0, node_ref_->getNodeConnection(0));
+      this->setExternalConnectionNodes(1, node_ref_->getNodeConnection(1));
+
+      return 0;
+    }
+
+    int initialize()
+    {
+      return 0;
+    }
+
+    int tagDifferentiable()
+    {
+      return 0;
+    }
+
+    int evaluateInternalResidual()
+    {
+
+      // Internals
+      f_int_[0] = -yp_int_[0] - 280 * y_int_[0] * y_int_[2] + 0.69 * *y_ext_[0] + 1.71 * *y_ext_[1] - 0.43 * y_int_[0] + 0.69 * y_int_[1];
+      f_int_[1] = -yp_int_[1] + 280 * y_int_[0] * y_int_[2] - 1.81 * y_int_[1];
+      f_int_[2] = -yp_int_[2] - 280 * y_int_[0] * y_int_[2] + 1.81 * y_int_[1];
+
+      return 0;
+    }
+
+    int evaluateExternalResidual()
+    {
+      // Externals
+      *f_ext_[0] += -0.02 * *y_ext_[0];
+      *f_ext_[1] += -0.045 * *y_ext_[1] + 0.43 * y_int_[0] + 0.43 * y_int_[1];
+
+      return 0;
+    }
+
+    int evaluateJacobian()
+    {
+      this->zeroJacMatrix();
+
+      // Internal Jacobian Entries [row 1]
+      std::vector<IdxT>  row = {2, 2, 2, 2, 2};
+      std::vector<IdxT>  col = {2, 3, 4, 0, 1};
+      std::vector<RealT> val = {-280 * y_int_[2] - 0.43 - alpha_, 0.69, -280 * y_int_[0], 0.69, 1.71};
+
+      this->setJacValues(row, col, val);
+
+      // Internal Jacobian Entries [row 2]
+      row = {3, 3, 3};
+      col = {2, 3, 4};
+      val = {280 * y_int_[2], -1.81 - alpha_, 280 * y_int_[0]};
+
+      this->setJacValues(row, col, val);
+
+      // Internal Jacobian Entries [row 3]
+      row = {4, 4, 4};
+      col = {2, 3, 4};
+      val = {-280 * y_int_[2], 1.81, -280 * y_int_[0] - alpha_};
+
+      this->setJacValues(row, col, val);
+
+      // External Jacobian Entries
+      row = {0, 1, 1, 1};
+      col = {0, 2, 3, 1};
+      val = {-0.02, 0.43, 0.43, -0.045};
+
+      this->setJacValues(row, col, val);
+
+      return 0;
+    }
+
+    int evaluateIntegrand()
+    {
+      return 0;
+    }
+
+    int initializeAdjoint()
+    {
+      return 0;
+    }
+
+    int evaluateAdjointResidual()
+    {
+      return 0;
+    }
+
+    int evaluateAdjointIntegrand()
+    {
+      return 0;
+    }
+
+    /**
+     * @brief Compute the absolute tolerance for each variable in the model
+     *
+     */
+    int setAbsoluteTolerance(RealT rel_tol)
+    {
+      abs_tol_.setToConst(static_cast<ScalarT>(rel_tol));
+      return 0;
+    }
+
+    CircuitComponent<ScalarT, IdxT>* clone() const
+    {
+      return new HiresComponent3<ScalarT, IdxT>(*this);
+    }
+
+  private:
+    NodeT* node_ref_;
+  };
+
+  namespace Testing
+  {
+
+    /**
+     * This example assembles the HIRES problem and partitions it to demonstrate
+     * that the partitioned system can reconstruct the full system.
+     *
+     * The HIRES problem consists of eight ODE equations and is divided into three
+     * components. HiresComponent1 contains equations 1, 2, and 3,
+     * HiresComponent3 contains equations 6, 7, and 8, and HiresBus contains
+     * equations 4 and 5. HiresBus plays a role similar to MicrogridBusDQ, where
+     * other components connected to the bus contribute terms to its equations.
+     *
+     * In the equations below, terms enclosed in parentheses ( ) are contributions
+     * from HiresComponent1, terms enclosed in square brackets [ ] are contributions
+     * from HiresComponent3, and the remaining terms in equations 4 and 5 belong
+     * to HiresBus.
+     *\f[
+     * f_1 = dy_1/dt + 1.71y_1 - 0.43y_2 - 8.32y_3 - 0.0007
+     *
+     * f_2 = dy_2/dt - 1.71y_1 + 8.75y_2
+     *
+     * f_3 = dy_3/dt + 10.03y_3 - 0.43y_4 - 0.035y_5
+     *
+     * f_4 = dy_4/dt + y_4 + (0.1y_4 - 8.32y_2 - 1.71y_3) + [0.02y_4]
+     *
+     * f_5 = dy_5/dt + y_5 + \left(0.7y_5 \right) + \left[0.045y_5 - 0.43y_6 - 0.43y_7 \right]
+     *
+     * f_6 = dy_6/dt - 280y_6y_8 + 0.69y_4 + 1.71y_5 - 0.43y_6 + 0.69y_7
+     *
+     * f_7 = dy_7/dt + 280y_6y_8 - 1.81y_7
+     *
+     * f_8 = dy_8/dt - 280y_6y_8 + 1.81y_7
+     *\f]
+     *
+     * The assembled system has the following structure:
+     *
+     * (HiresComponent1) -------- (HiresBus) -------- (HiresComponent3)
+     *
+     * The system is partitioned between HiresBus and HiresComponent3. A bus
+     * partition interface is introduced in the partition containing HiresBus to
+     * preserve the contribution of HiresComponent3 to the bus equations. The
+     * residuals evaluated independently by the partitions are then printed with
+     * the residual of the full system for eye-ball comparision.
+     *
+     * This example will also be useful later for testing the order of accuracy of co-simulation methods.
+     */
+    template <class ScalarT, typename IdxT>
+    class SubsystemModelWithHiresTests
+    {
+      using RealT     = typename CircuitComponent<ScalarT, IdxT>::RealT;
+      using Bus       = PowerElectronics::MicrogridBus<ScalarT, IdxT>;
+      using Subsystem = SubsystemModel<ScalarT, IdxT>;
+      using System    = PowerElectronicsModel<ScalarT, IdxT>;
+
+    public:
+      /**
+       * @brief Construct the HIRES reference system and subsystem partitions.
+       */
+      SubsystemModelWithHiresTests()
+        : system_(new System()),
+          partition1_(new Subsystem()),
+          partition2_(new Subsystem()),
+          comp1_(new HiresComponent1<ScalarT, IdxT>(&bus_, 1)),
+          bus1_(new HiresBus<ScalarT, IdxT>(&bus_, 2)),
+          comp3_(new HiresComponent3<ScalarT, IdxT>(&bus_, 3))
+      {
+
+        y_  = {1, 2, 3, 4, 5, 6, 7, 8};
+        yp_ = {1, 2, 3, 4, 5, 6, 7, 8};
+        // ---------------------------------------------------------------------
+        // Assemble and allocate the monolithic reference system
+        // ---------------------------------------------------------------------
+
+        system_->addComponent(comp1_);
+        system_->addComponent(comp3_);
+        system_->addComponent(bus1_);
+        system_->addNode(&bus_);
+
+        system_->allocate();
+
+        distributeVariables(y_, yp_);
+        system_->updateTime(1.0, 2.0);
+        system_->evaluateResidual();
+
+        // ---------------------------------------------------------------------
+        // Construct the partition interface
+        // ---------------------------------------------------------------------
+
+        auto* comp3_copy = new HiresComponent3<ScalarT, IdxT>(*comp3_);
+        bus_interface_   = new BusPartitionInterface<ScalarT, IdxT>(&bus_, comp3_copy, 4);
+
+        bus_interface_->allocate();
+
+        // ---------------------------------------------------------------------
+        // Assemble the subsystem partitions
+        // ---------------------------------------------------------------------
+
+        partition1_->addComponent(comp1_);
+        partition1_->addComponent(bus1_);
+        partition1_->addInterface(bus_interface_);
+        partition1_->addNode(&bus_);
+
+        partition2_->addComponent(comp3_);
+
+        partition1_->allocate();
+        partition2_->allocate();
+        partitions_ = {partition1_, partition2_};
+      }
+
+      ~SubsystemModelWithHiresTests()
+      {
+        delete partition1_;
+        delete partition2_;
+        delete system_;
+      }
+
+      void distributeVariables(const std::vector<ScalarT>& y, const std::vector<ScalarT>& yp)
+      {
+        auto* system_y  = system_->y().getData();
+        auto* system_yp = system_->yp().getData();
+
+        for (size_t i = 0; i < system_->size(); ++i)
+        {
+          system_y[i]  = y[i];
+          system_yp[i] = yp[i];
+        }
+
+        system_->y().setDataUpdated();
+        system_->yp().setDataUpdated();
+
+        for (auto* partition : partitions_)
+        {
+          for (size_t i = 0; i < partition->getExternSize(); ++i)
+          {
+            const auto global_index = partition->getExternalDataIndices()[i];
+
+            partition->getExternalDataY()[i]  = y[global_index];
+            partition->getExternalDataYP()[i] = yp[global_index];
+          }
+
+          auto* partition_y  = partition->y().getData();
+          auto* partition_yp = partition->yp().getData();
+
+          for (size_t i = 0; i < partition->getInternalSize(); ++i)
+          {
+            const auto global_index = partition->getNodeConnection(i);
+
+            partition_y[i]  = y[global_index];
+            partition_yp[i] = yp[global_index];
+          }
+
+          partition->y().setDataUpdated();
+          partition->yp().setDataUpdated();
+        }
+      }
+
+      /**
+       * @brief Verify that the partitioned HIRES residual matches the
+       *        monolithic residual.
+       */
+      TestOutcome residual()
+      {
+        TestStatus success = true;
+
+        // Distribute variables to all partitions
+        distributeVariables(y_, yp_);
+
+        std::vector<ScalarT> partition_residual(system_->size(), 0.0);
+
+        for (auto* partition : partitions_)
+        {
+          partition->updateTime(1.0, 2.0);
+          partition->evaluateResidual();
+
+          auto* residual = partition->getResidual().getData();
+
+          for (size_t i = 0; i < partition->getInternalSize(); ++i)
+          {
+            partition_residual[partition->getNodeConnection(i)] = residual[i];
+          }
+
+          partition->getResidual().setDataUpdated();
+        }
+
+        auto* reference_residual = system_->getResidual().getData();
+
+        RealT max_error = 0.0;
+
+        for (size_t i = 0; i < system_->size(); ++i)
+        {
+          double error = std::abs(partition_residual[i] - reference_residual[i]) / std::abs(reference_residual[i] + 1);
+          max_error    = std::max(max_error, error);
+        }
+
+        std::cout << "max error " << max_error << std::endl;
+
+        success *= max_error <= std::numeric_limits<RealT>::epsilon();
+
+        return success.report(__func__);
+      }
+
+      /**
+       * @brief Verify that each subsystem Jacobian matches the corresponding
+       *        entries of the monolithic HIRES Jacobian.
+       */
+      TestOutcome jacobian()
+      {
+        TestStatus success = true;
+
+        RealT alpha = 2.0;
+
+        distributeVariables(y_, yp_);
+
+        /*
+         * GridKit stores the HIRES variables in component assembly order:
+         *
+         *   System index:  0   1   2   3   4   5   6   7
+         *   HIRES index:   0   1   2   5   6   7   3   4
+         *   Variable:      y1  y2  y3  y6  y7  y8  y4  y5
+         */
+        const std::array<IdxT, 8> sysmodel_to_hires = {
+            0, 1, 2, 5, 6, 7, 3, 4};
+
+        const std::array<IdxT, 8> hires_to_sysmodel = {
+            0, 1, 2, 6, 7, 3, 4, 5};
+
+        const RealT y6 = y_[hires_to_sysmodel[5]];
+        const RealT y8 = y_[hires_to_sysmodel[7]];
+
+        std::array<std::array<RealT, 8>, 8> reference_jac =
+            {{{-alpha - 1.71, 0.43, 8.32, 0.0, 0.0, 0.0, 0.0, 0.0},
+              {1.71, -alpha - 8.75, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
+              {0.0, 0.0, -alpha - 10.03, 0.43, 0.035, 0.0, 0.0, 0.0},
+              {0.0, 8.32, 1.71, -alpha - 1.12, 0.0, 0.0, 0.0, 0.0},
+              {0.0, 0.0, 0.0, 0.0, -alpha - 1.745, 0.43, 0.43, 0.0},
+              {0.0, 0.0, 0.0, 0.69, 1.71, -alpha - 280.0 * y8 - 0.43, 0.69, -280.0 * y6},
+              {0.0, 0.0, 0.0, 0.0, 0.0, 280.0 * y8, -alpha - 1.81, 280.0 * y6},
+              {0.0, 0.0, 0.0, 0.0, 0.0, -280.0 * y8, 1.81, -alpha - 280.0 * y6}}};
+
+        for (auto* partition : partitions_)
+        {
+          partition->updateTime(1.0, alpha);
+          partition->evaluateJacobian();
+
+          auto* partition_jac = partition->getCsrJacobian();
+
+          const auto* row_ptr = partition_jac->getRowData();
+          const auto* cols    = partition_jac->getColData();
+          const auto* vals    = partition_jac->getValues();
+
+          const size_t n = partition->getInternalSize();
+
+          std::vector<std::vector<RealT>> dense_jac(n, std::vector<RealT>(n, 0.0));
+
+          // Convert the partition CSR Jacobian to a dense matrix.
+          for (size_t row = 0; row < n; ++row)
+          {
+            for (IdxT k = row_ptr[row]; k < row_ptr[row + 1]; ++k)
+            {
+              dense_jac[row][cols[k]] = vals[k];
+            }
+          }
+
+          // Compare with the corresponding entries of the full Jacobian.
+          for (size_t row = 0; row < n; ++row)
+          {
+            const IdxT global_row = partition->getNodeConnection(row);
+            const IdxT ref_row    = sysmodel_to_hires[global_row];
+
+            for (size_t col = 0; col < n; ++col)
+            {
+              const IdxT global_col = partition->getNodeConnection(col);
+              const IdxT ref_col    = sysmodel_to_hires[global_col];
+
+              success *= std::abs(dense_jac[row][col] - reference_jac[ref_row][ref_col]) <= std::numeric_limits<RealT>::epsilon();
+            }
+          }
+        }
+
+        return success.report(__func__);
+      }
+
+    private:
+      Bus bus_;
+
+      System* system_;
+
+      Subsystem* partition1_;
+      Subsystem* partition2_;
+
+      HiresComponent1<ScalarT, IdxT>* comp1_;
+      HiresBus<ScalarT, IdxT>*        bus1_;
+      HiresComponent3<ScalarT, IdxT>* comp3_;
+
+      BusPartitionInterface<ScalarT, IdxT>* bus_interface_;
+
+      std::vector<Subsystem*> partitions_;
+
+      std::vector<ScalarT> y_;
+      std::vector<ScalarT> yp_;
+    };
+
+  } // namespace Testing
+
+} // namespace GridKit
diff --git a/tests/UnitTests/PowerElectronics/runComponentCloneTests.cpp b/tests/UnitTests/PowerElectronics/runComponentCloneTests.cpp
new file mode 100644
index 000000000..dfb5fa97f
--- /dev/null
+++ b/tests/UnitTests/PowerElectronics/runComponentCloneTests.cpp
@@ -0,0 +1,15 @@
+#include "ComponentCloneTests.hpp"
+
+int main()
+{
+  GridKit::Testing::CircuitComponentCloneTests<double, size_t> tests;
+
+  GridKit::Testing::TestingResults result;
+
+  result += tests.distributedGeneratorClone();
+  result += tests.microgridLineClone();
+  result += tests.microgridLoadClone();
+  result += tests.microgridBusDQClone();
+
+  return result.summary();
+}
diff --git a/tests/UnitTests/PowerElectronics/runSubsystemModelWithHiresTest.cpp b/tests/UnitTests/PowerElectronics/runSubsystemModelWithHiresTest.cpp
new file mode 100644
index 000000000..2b8ef11b0
--- /dev/null
+++ b/tests/UnitTests/PowerElectronics/runSubsystemModelWithHiresTest.cpp
@@ -0,0 +1,12 @@
+#include "SubsystemModelWithHiresTest.hpp"
+
+int main()
+{
+  GridKit::Testing::TestingResults                               result;
+  GridKit::Testing::SubsystemModelWithHiresTests<double, size_t> test;
+
+  result += test.residual();
+  result += test.jacobian();
+
+  return result.summary();
+}