Jacobian via dependency tracking for LoadTest.

Author: nkoukpaizan

src/Model/PhasorDynamics/Bus/BusInfinite.cpp

@@ -195,6 +195,7 @@ namespace GridKit
     // Available template instantiations
     template class BusInfinite<double, long int>;
     template class BusInfinite<double, size_t>;
+    template class BusInfinite<Sparse::Variable, long int>;
     template class BusInfinite<Sparse::Variable, size_t>;
 
   } // namespace PhasorDynamics

src/Model/PhasorDynamics/BusBase.hpp

@@ -2,8 +2,8 @@
 
 #include <vector>
 
-#include <Model/Evaluator.hpp>
 #include <LinearAlgebra/SparsityPattern/Variable.hpp>
+#include <Model/Evaluator.hpp>
 
 namespace GridKit
 {

src/Model/PhasorDynamics/Component.hpp

@@ -2,6 +2,7 @@
 
 #include <vector>
 
+#include <LinearAlgebra/SparsityPattern/Variable.hpp>
 #include <Model/Evaluator.hpp>
 
 namespace GridKit

src/Model/PhasorDynamics/Load/Load.cpp

@@ -173,6 +173,8 @@ namespace GridKit
     // Available template instantiations
     template class Load<double, long int>;
     template class Load<double, size_t>;
+    template class Load<Sparse::Variable, long int>;
+    template class Load<Sparse::Variable, size_t>;
 
   } // namespace PhasorDynamics
 } // namespace GridKit

src/Utilities/Testing.hpp

@@ -232,7 +232,6 @@ namespace GridKit
       for (const auto& pair_a : a)
       {
         auto it = b.find(pair_a.first);
-        std::cout << pair_a.second << ", " << it->second << std::endl;
         fail += !isEqual(pair_a.second, it->second);
       }
 

tests/UnitTests/PhasorDynamics/LoadTests.hpp

@@ -3,6 +3,7 @@
 #include <iomanip>
 #include <iostream>
 
+#include <LinearAlgebra/SparsityPattern/Variable.hpp>
 #include <Model/PhasorDynamics/Bus/Bus.hpp>
 #include <Model/PhasorDynamics/Bus/BusInfinite.hpp>
 #include <Model/PhasorDynamics/Load/Load.hpp>
@@ -65,6 +66,58 @@ namespace GridKit
 
         return success.report(__func__);
       }
+
+      TestOutcome jacobian()
+      {
+        TestStatus success = true;
+
+        real_type R{2.0}; ///< Load resistance
+        real_type X{4.0}; ///< Load reactance
+
+        Sparse::Variable Vr{10.0}; ///< Bus real voltage
+        Sparse::Variable Vi{20.0}; ///< Bus imaginary voltage
+                                   
+        Vr.setVariableNumber(0); ///< Independent variables: first
+        Vi.setVariableNumber(1); ///< Independent variables: second
+
+        PhasorDynamics::BusInfinite<Sparse::Variable, IdxT> bus(Vr, Vi);
+
+        PhasorDynamics::Load<Sparse::Variable, IdxT> load(&bus, R, X);
+        load.evaluateResidual(); ///< Computes the residual and the Jacobian values by tracking 
+                                 ///< the dependencies
+
+        std::vector<Sparse::Variable> residuals{bus.Ir(), bus.Ii()};
+        std::vector<Sparse::Variable::DependencyMap> ref = analyticalJacobian(R, X);
+        
+        /// Compare dependencies computed automatically to the ones computed analytically
+        for (size_t i = 0; i < residuals.size(); ++i)
+        {
+          Sparse::Variable res = residuals[i];
+          const Sparse::Variable::DependencyMap& dependencies = res.getDependencies();
+          success *= (GridKit::Testing::isEqual(dependencies, ref[i])); 
+        }
+
+        return success.report(__func__);
+      }
+    private:
+      std::vector<Sparse::Variable::DependencyMap> analyticalJacobian(const real_type R,
+                                                                      const real_type X)
+      {
+        const real_type b = -X / (R * R + X * X);
+        const real_type g = R / (R * R + X * X);
+
+        real_type dIr_dVr = -g;
+        real_type dIr_dVi = -b;
+
+        real_type dIi_dVr = b; 
+        real_type dIi_dVi = -g;
+                              
+        std::vector<Sparse::Variable::DependencyMap> dependencies(2);
+        dependencies[0] = {{0, dIr_dVr}, {1, dIr_dVi}};
+        dependencies[1] = {{0, dIi_dVr}, {1, dIi_dVi}};
+
+        return dependencies;
+      }
     };
 
   } // namespace Testing

tests/UnitTests/PhasorDynamics/runLoadTests.cpp

@@ -10,6 +10,7 @@ int main()
 
   result += test.constructor();
   result += test.residual();
+  result += test.jacobian();
 
   return result.summary();
 }