diff --git a/CHANGELOG.md b/CHANGELOG.md
index ca382a44..b8dfddd8 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -41,6 +41,7 @@
 - Added a header file defining common math functions (e.g., sigmoid) to be used throughout the code.
 - Added capability to print monitored variables in multiple formats, triggered from `Ida::runSimulation`.
 - Added IDA statistics object which can be accumulated over multiple simulations.
+- Minor performance improvements to residual evaluation in PowerElectronics module.
 
 ## v0.1
 
diff --git a/GridKit/Model/PowerElectronics/CircuitComponent.hpp b/GridKit/Model/PowerElectronics/CircuitComponent.hpp
index 9d7f53ee..55bcf3f1 100644
--- a/GridKit/Model/PowerElectronics/CircuitComponent.hpp
+++ b/GridKit/Model/PowerElectronics/CircuitComponent.hpp
@@ -21,10 +21,20 @@ namespace GridKit
     using RealT   = typename Model::Evaluator<ScalarT, IdxT>::RealT;
     using MatrixT = typename Model::Evaluator<ScalarT, IdxT>::MatrixT;
 
-    CircuitComponent()  = default;
-    ~CircuitComponent() = default;
+    CircuitComponent() = default;
 
-    void updateTime(RealT t, RealT a)
+    ~CircuitComponent()
+    {
+      if (connection_nodes_ != nullptr)
+      {
+        delete[] connection_nodes_;
+      }
+    };
+
+    /**
+     * @note Cannot be marked final, since it is overriden to recurse in the system model.
+     */
+    void updateTime(RealT t, RealT a) override
     {
       this->time_  = t;
       this->alpha_ = a;
@@ -59,6 +69,11 @@ namespace GridKit
      */
     int setExternalConnectionNodes(IdxT local_index, IdxT global_index)
     {
+      if (connection_nodes_ == nullptr)
+      {
+        connection_nodes_ = new IdxT[size_];
+      }
+
       connection_nodes_[local_index] = global_index;
       return 0;
     }
@@ -71,185 +86,205 @@ namespace GridKit
      * @param local_index index of local value in vector
      * @return IdxT Index of the same value in the global vector
      */
-    IdxT getNodeConnection(IdxT local_index)
+    IdxT getNodeConnection(IdxT local_index) const
     {
-      return connection_nodes_.at(local_index);
+      return connection_nodes_[local_index];
     }
 
   public:
-    virtual IdxT size()
+    IdxT size() final
     {
       return size_;
     }
 
-    virtual IdxT nnz()
+    IdxT size() const
+    {
+      return size_;
+    }
+
+    IdxT nnz() final
     {
       return nnz_;
     }
 
-    virtual IdxT sizeQuadrature()
+    IdxT nnz() const
+    {
+      return nnz_;
+    }
+
+    IdxT sizeQuadrature() final
+    {
+      return size_quad_;
+    }
+
+    IdxT sizeQuadrature() const
     {
       return size_quad_;
     }
 
-    virtual IdxT sizeParams()
+    IdxT sizeParams() final
+    {
+      return size_opt_;
+    }
+
+    IdxT sizeParams() const
     {
       return size_opt_;
     }
 
-    virtual void setTolerances(RealT& rel_tol, RealT& abs_tol) const
+    void setTolerances(RealT& rel_tol, RealT& abs_tol) const final
     {
       rel_tol = rel_tol_;
       abs_tol = abs_tol_;
     }
 
-    virtual void setMaxSteps(IdxT& msa) const
+    void setMaxSteps(IdxT& msa) const final
     {
       msa = max_steps_;
     }
 
-    std::vector<ScalarT>& y()
+    std::vector<ScalarT>& y() final
     {
       return y_;
     }
 
-    const std::vector<ScalarT>& y() const
+    const std::vector<ScalarT>& y() const final
     {
       return y_;
     }
 
-    std::vector<ScalarT>& yp()
+    std::vector<ScalarT>& yp() final
     {
       return yp_;
     }
 
-    const std::vector<ScalarT>& yp() const
+    const std::vector<ScalarT>& yp() const final
     {
       return yp_;
     }
 
-    std::vector<bool>& tag()
+    std::vector<bool>& tag() final
     {
       return tag_;
     }
 
-    const std::vector<bool>& tag() const
+    const std::vector<bool>& tag() const final
     {
       return tag_;
     }
 
-    std::vector<ScalarT>& yB()
+    std::vector<ScalarT>& yB() final
     {
       return yB_;
     }
 
-    const std::vector<ScalarT>& yB() const
+    const std::vector<ScalarT>& yB() const final
     {
       return yB_;
     }
 
-    std::vector<ScalarT>& ypB()
+    std::vector<ScalarT>& ypB() final
     {
       return ypB_;
     }
 
-    const std::vector<ScalarT>& ypB() const
+    const std::vector<ScalarT>& ypB() const final
     {
       return ypB_;
     }
 
-    std::vector<ScalarT>& param()
+    std::vector<ScalarT>& param() final
     {
       return param_;
     }
 
-    const std::vector<ScalarT>& param() const
+    const std::vector<ScalarT>& param() const final
     {
       return param_;
     }
 
-    std::vector<ScalarT>& param_up()
+    std::vector<ScalarT>& param_up() final
     {
       return param_up_;
     }
 
-    const std::vector<ScalarT>& param_up() const
+    const std::vector<ScalarT>& param_up() const final
     {
       return param_up_;
     }
 
-    std::vector<ScalarT>& param_lo()
+    std::vector<ScalarT>& param_lo() final
     {
       return param_lo_;
     }
 
-    const std::vector<ScalarT>& param_lo() const
+    const std::vector<ScalarT>& param_lo() const final
     {
       return param_lo_;
     }
 
-    std::vector<ScalarT>& getResidual()
+    std::vector<ScalarT>& getResidual() final
     {
       return f_;
     }
 
-    const std::vector<ScalarT>& getResidual() const
+    const std::vector<ScalarT>& getResidual() const final
     {
       return f_;
     }
 
-    MatrixT& getJacobian()
+    MatrixT& getJacobian() final
     {
       return jac_;
     }
 
-    const MatrixT& getJacobian() const
+    const MatrixT& getJacobian() const final
     {
       return jac_;
     }
 
-    std::vector<ScalarT>& getIntegrand()
+    std::vector<ScalarT>& getIntegrand() final
     {
       return g_;
     }
 
-    const std::vector<ScalarT>& getIntegrand() const
+    const std::vector<ScalarT>& getIntegrand() const final
     {
       return g_;
     }
 
-    std::vector<ScalarT>& getAdjointResidual()
+    std::vector<ScalarT>& getAdjointResidual() final
     {
       return fB_;
     }
 
-    const std::vector<ScalarT>& getAdjointResidual() const
+    const std::vector<ScalarT>& getAdjointResidual() const final
     {
       return fB_;
     }
 
-    std::vector<ScalarT>& getAdjointIntegrand()
+    std::vector<ScalarT>& getAdjointIntegrand() final
     {
       return gB_;
     }
 
-    const std::vector<ScalarT>& getAdjointIntegrand() const
+    const std::vector<ScalarT>& getAdjointIntegrand() const final
     {
       return gB_;
     }
 
     //@todo Fix ID naming
-    IdxT getIDcomponent()
+    IdxT getIDcomponent() const
     {
       return idc_;
     }
 
   protected:
-    size_t               n_extern_;
-    size_t               n_intern_;
-    std::set<IdxT>       extern_indices_;
+    size_t         n_extern_;
+    size_t         n_intern_;
+    std::set<IdxT> extern_indices_;
     ///@todo may want to replace the mapping of connection_nodes to Node objects instead of IdxT. Allows for container free setup
-    std::map<IdxT, IdxT> connection_nodes_;
+    IdxT*          connection_nodes_ = nullptr;
 
   protected:
     IdxT size_{0};
diff --git a/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp b/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp
index 736f9d00..1bffd0d8 100644
--- a/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp
+++ b/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp
@@ -133,7 +133,7 @@ namespace GridKit
      *
      * @return int
      */
-    int allocate()
+    int allocate() final
     {
       return 1;
     }
@@ -145,7 +145,7 @@ namespace GridKit
      * @return true if all components have jacobian
      * @return false otherwise
      */
-    bool hasJacobian()
+    bool hasJacobian() final
     {
       if (!this->use_jac_)
         return false;
@@ -192,7 +192,7 @@ namespace GridKit
      *
      * @return int 0 if successful, positive if there's a recoverable error, negative if unrecoverable
      */
-    int initialize()
+    int initialize() final
     {
       // Initialize components
       for (const auto& component : components_)
@@ -215,17 +215,21 @@ namespace GridKit
     {
       for (const auto& component : components_)
       {
-        for (IdxT j = 0; j < component->size(); ++j)
+        IdxT                  size = component->size();
+        std::vector<ScalarT>& y    = component->y();
+        std::vector<ScalarT>& yp   = component->yp();
+
+        for (IdxT j = 0; j < size; ++j)
         {
           if (component->getNodeConnection(j) != neg1_)
           {
-            component->y()[j]  = y_[component->getNodeConnection(j)];
-            component->yp()[j] = yp_[component->getNodeConnection(j)];
+            y[j]  = y_[component->getNodeConnection(j)];
+            yp[j] = yp_[component->getNodeConnection(j)];
           }
           else
           {
-            component->y()[j]  = 0.0;
-            component->yp()[j] = 0.0;
+            y[j]  = 0.0;
+            yp[j] = 0.0;
           }
         }
       }
@@ -242,7 +246,7 @@ namespace GridKit
      *
      * @return int 0 if successful, positive if there's a recoverable error, negative if unrecoverable
      */
-    int evaluateResidual()
+    int evaluateResidual() final
     {
       for (IdxT i = 0; i < this->f_.size(); i++)
       {
@@ -257,12 +261,15 @@ namespace GridKit
       {
         // TODO:check return type
         component->evaluateResidual();
-        for (IdxT j = 0; j < component->size(); ++j)
+
+        IdxT                        size     = component->size();
+        const std::vector<ScalarT>& residual = component->getResidual();
+        for (IdxT j = 0; j < size; ++j)
         {
           //@todo should do a different grounding check
           if (component->getNodeConnection(j) != neg1_)
           {
-            f_[component->getNodeConnection(j)] += component->getResidual()[j];
+            f_[component->getNodeConnection(j)] += residual[j];
           }
         }
       }
@@ -278,7 +285,7 @@ namespace GridKit
      *
      * @return int 0 if successful, positive if there's a recoverable error, negative if unrecoverable
      */
-    int evaluateJacobian()
+    int evaluateJacobian() final
     {
       jac_.zeroMatrix();
       distributeVectors();
@@ -320,7 +327,7 @@ namespace GridKit
     /**
      * @brief Evaluate integrands for the system quadratures.
      */
-    int evaluateIntegrand()
+    int evaluateIntegrand() final
     {
 
       return 0;
@@ -332,7 +339,7 @@ namespace GridKit
      * Updates variables and optimization parameters, then initializes
      * adjoints locally and copies them to the system adjoint vector.
      */
-    int initializeAdjoint()
+    int initializeAdjoint() final
     {
       return 0;
     }
@@ -342,7 +349,7 @@ namespace GridKit
      *
      *
      */
-    int evaluateAdjointResidual()
+    int evaluateAdjointResidual() final
     {
       return 0;
     }
@@ -352,7 +359,7 @@ namespace GridKit
      *
      *
      */
-    int evaluateAdjointIntegrand()
+    int evaluateAdjointIntegrand() final
     {
       return 0;
     }
@@ -363,7 +370,7 @@ namespace GridKit
      * @param t
      * @param a
      */
-    void updateTime(RealT t, RealT a)
+    void updateTime(RealT t, RealT a) final
     {
       for (const auto& component : components_)
       {