Unit test and example for Genrou class.

examples/DistributedGeneratorTest/DGTest.cpp

@@ -1,5 +1,5 @@
 
-
+#define _USE_MATH_DEFINES
 #include <cmath>
 #include <filesystem>
 #include <fstream>

examples/Grid3Bus/Grid3BusSys.cpp

@@ -9,6 +9,7 @@
  *
  */
 
+#define _USE_MATH_DEFINES
 #include <cmath>
 #include <filesystem>
 #include <fstream>

examples/Microgrid/Microgrid.cpp

@@ -1,5 +1,6 @@
 
 
+#define _USE_MATH_DEFINES
 #include <cmath>
 #include <filesystem>
 #include <fstream>

examples/PhasorDynamics/CMakeLists.txt

@@ -1 +1,2 @@
 add_subdirectory(Example1)
+add_subdirectory(Example2)

examples/PhasorDynamics/Example1/CMakeLists.txt

@@ -4,11 +4,7 @@ target_link_libraries(phasordynamics_example1
 	GRIDKIT::phasor_dynamics_bus_fault 
 	GRIDKIT::phasor_dynamics_branch
 	GRIDKIT::phasor_dynamics_genrou
-	SUNDIALS::sunlinsolklu
-	SUNDIALS::core
-	SUNDIALS::ida
-	SUNDIALS::idas
-	SUNDIALS::sunmatrixdense)
+	GRIDKIT::solvers_dyn)
 install(TARGETS phasordynamics_example1 RUNTIME DESTINATION bin)
 
-add_test(NAME GenrouTest COMMAND $<TARGET_FILE:phasordynamics_example1>)
+add_test(NAME GenrouTest1 COMMAND $<TARGET_FILE:phasordynamics_example1>)

examples/PhasorDynamics/Example1/example1.cpp

@@ -1,37 +1,33 @@
-#include <stdio.h>
-#define _USE_MATH_DEFINES
-#include <math.h>
-#include <time.h>
-
-// #include <sundials_core.h>
-#include <idas/idas.h>
-#include <nvector/nvector_serial.h>
-#include <sunlinsol/sunlinsol_dense.h>
-#include <sunlinsol/sunlinsol_klu.h>
-#include <sunmatrix/sunmatrix_sparse.h>
-
-#include "Example1_Powerworld_Reference.hpp"
-#include "Model/PhasorDynamics/Branch/Branch.cpp"
-#include "Model/PhasorDynamics/Branch/Branch.hpp"
-#include "Model/PhasorDynamics/Bus/Bus.cpp"
-#include "Model/PhasorDynamics/Bus/Bus.hpp"
-#include "Model/PhasorDynamics/Bus/BusInfinite.cpp"
-#include "Model/PhasorDynamics/Bus/BusInfinite.hpp"
-#include "Model/PhasorDynamics/BusFault/BusFault.hpp"
-#include "Model/PhasorDynamics/SynchronousMachine/GENROUwS/Genrou.hpp"
-#include "Model/PhasorDynamics/SystemModel.hpp"
-#include "Solver/Dynamic/Ida.cpp"
-#include "Solver/Dynamic/Ida.hpp"
+/**
+ * @file example1.cpp
+ * @author Adam Birchfield (abirchfield@tamu.edu)
+ * @author Slaven Peles (peless@ornl.gov)
+ * @brief Example running a 2-bus system
+ *
+ * Simulates a 2-bus system with Genrou 6th order generator model and
+ * compares results with data generated for the same system by Poweworld.
+ *
+ */
+#include "example1.hpp"
+
+#include <ctime>
+#include <iostream>
+
+#include <Model/PhasorDynamics/Branch/Branch.hpp>
+#include <Model/PhasorDynamics/Bus/Bus.hpp>
+#include <Model/PhasorDynamics/Bus/BusInfinite.hpp>
+#include <Model/PhasorDynamics/BusFault/BusFault.hpp>
+#include <Model/PhasorDynamics/SynchronousMachine/GENROUwS/Genrou.hpp>
+#include <Model/PhasorDynamics/SystemModel.hpp>
+#include <Solver/Dynamic/Ida.hpp>
 #include <Utilities/Testing.hpp>
 
-#define _CRT_SECURE_NO_WARNINGS
-
 int main()
 {
   using namespace GridKit::PhasorDynamics;
   using namespace AnalysisManager::Sundials;
 
-  printf("Example 1 version 2\n");
+  std::cout << "Example 1 version 2\n";
 
   /* Create model parts */
   SystemModel<double, size_t> sys;
@@ -115,7 +111,6 @@ int main()
           / (1.0 + std::abs(reference_solution[i / 48][2]));
       if (err > error_V)
         error_V = err;
-      // std::cout << "t = " << ti << ": Vr = " << Vr << ", Vi = " << Vi << ", dw = " << dw;
       std::cout << "GridKit: t = " << ti
                 << ", |V| = " << std::sqrt(Vr * Vr + Vi * Vi)
                 << ", w = " << (1.0 + dw) << "\n";
@@ -149,11 +144,8 @@ int main()
     }
     ++j;
     ++i;
-    // if (i > 500)
-    //   break;
   }
 
-  // std::cout << buffer.str();
   int status = 0;
   std::cout << "Max error in |V| = " << error_V << "\n";
   if (error_V > 2e-4)
@@ -164,5 +156,5 @@ int main()
 
   std::cout << "\n\nComplete in " << (stop - start) / CLOCKS_PER_SEC << " seconds\n";
 
-  return 0;
+  return status;
 }

examples/PhasorDynamics/Example1/Example1_Powerworld_Reference.hpp → examples/PhasorDynamics/Example1/example1.hpp

@@ -1,3 +1,10 @@
+/**
+ * @file example1.hpp
+ * @author Adam Birchfield (abirchfield@tamu.edu)
+ * @author Slaven Peles (peless@ornl.gov)
+ * @brief Reference solution for the 2-bus system obtained with Powerworld
+ *
+ */
 #include <vector>
 
 // Columns:

examples/PhasorDynamics/Example2/CMakeLists.txt

@@ -0,0 +1,11 @@
+add_executable(phasordynamics_example2 example2.cpp)
+target_link_libraries(phasordynamics_example2 
+	GRIDKIT::phasor_dynamics_bus 
+	GRIDKIT::phasor_dynamics_bus_fault 
+	GRIDKIT::phasor_dynamics_branch
+	GRIDKIT::phasor_dynamics_genrou
+	GRIDKIT::phasor_dynamics_load
+	GRIDKIT::solvers_dyn)
+install(TARGETS phasordynamics_example2 RUNTIME DESTINATION bin)
+
+add_test(NAME GenrouTest2 COMMAND $<TARGET_FILE:phasordynamics_example2>)

examples/PhasorDynamics/Example2/example2.cpp

@@ -0,0 +1,141 @@
+/**
+ * @file example2.cpp
+ * @author Adam Birchfield (abirchfield@tamu.edu)
+ * @author Slaven Peles (peless@ornl.gov)
+ * @brief Example running a 3-bus system
+ *
+ * Simulates a 3-bus system with two Genrou 6th order generator models and
+ * compares results with data generated for the same system by Poweworld.
+ *
+ */
+#include "example2.hpp"
+
+#include <cstdio>
+#include <ctime>
+#include <fstream>
+#include <vector>
+
+#include <Model/PhasorDynamics/Branch/Branch.hpp>
+#include <Model/PhasorDynamics/Bus/Bus.hpp>
+#include <Model/PhasorDynamics/Bus/BusInfinite.hpp>
+#include <Model/PhasorDynamics/BusFault/BusFault.hpp>
+#include <Model/PhasorDynamics/Load/Load.hpp>
+#include <Model/PhasorDynamics/SynchronousMachine/GENROUwS/Genrou.hpp>
+#include <Model/PhasorDynamics/SystemModel.hpp>
+#include <Solver/Dynamic/Ida.hpp>
+#include <Utilities/Testing.hpp>
+
+int main()
+{
+  using namespace GridKit::PhasorDynamics;
+  using namespace AnalysisManager::Sundials;
+  using scalar_type = double;
+  using real_type   = double;
+  using index_type  = size_t;
+
+  auto error_allowed = static_cast<real_type>(1e-4);
+
+  std::cout << "Example 2 version 1\n";
+
+  /* Create model parts */
+  SystemModel<scalar_type, index_type> sys;
+  BusInfinite<scalar_type, index_type> bus1(1.06, 0.0);
+  Bus<scalar_type, index_type>         bus2(1.0599558398065716, -0.009675621941024773);
+  Bus<scalar_type, index_type>         bus3(0.9610827543495831, -0.13122476630506485);
+  Branch<scalar_type, index_type>      branch12(&bus1, &bus2, 0.05, 0.21, 0, 0.1);
+  Branch<scalar_type, index_type>      branch13(&bus1, &bus3, 0.06, 0.15, 0, 0.12);
+  Branch<scalar_type, index_type>      branch23(&bus2, &bus3, 0.08, 0.27, 0, 0.45);
+  Genrou<scalar_type, index_type>      gen2(&bus2, 1, 0.5, -0.07588, 2.7, 0., 0., 7., .04, .05, .75, 1.9, 0.17, 0.15, 0.4, 0.35, 0.15, 0.14999, 0., 0.);
+  Genrou<scalar_type, index_type>      gen3(&bus3, 1, 0.25, 0.26587, 1.6, 0., 0., 7.5, .04, .05, .75, 2.3, 0.2, 0.18, 0.5, 0.5, 0.18, 0.15, 0., 0.);
+  Load<scalar_type, index_type>        load3(&bus3, 0.4447197839297772, 0.20330047265361242);
+  BusFault<scalar_type, index_type>    fault(&bus3, 0, 1e-5, 0);
+
+  /* Connect everything together */
+  sys.addBus(&bus1);
+  sys.addBus(&bus2);
+  sys.addBus(&bus3);
+  sys.addComponent(&branch12);
+  sys.addComponent(&branch13);
+  sys.addComponent(&branch23);
+  sys.addComponent(&fault);
+  sys.addComponent(&load3);
+  sys.addComponent(&gen2);
+  sys.addComponent(&gen3);
+  sys.allocate();
+
+  real_type dt = 1.0 / 4.0 / 60.0;
+
+  std::stringstream buffer;
+
+  /* Set up simulation */
+  Ida<scalar_type, index_type> ida(&sys);
+  ida.configureSimulation();
+
+  /* Run simulation */
+  scalar_type start = static_cast<scalar_type>(clock());
+  ida.initializeSimulation(0.0, false);
+  ida.runSimulationFixed(0.0, dt, 1.0, buffer);
+  fault.setStatus(1);
+  ida.initializeSimulation(1.0, false);
+  ida.runSimulationFixed(1.0, dt, 1.1, buffer);
+  fault.setStatus(0);
+  ida.initializeSimulation(1.1, false);
+  ida.runSimulationFixed(1.1, dt, 10.0, buffer);
+  double stop = static_cast<double>(clock());
+
+  /* Check worst-case error */
+  real_type worst_error      = 0;
+  real_type worst_error_time = 0;
+
+  const index_type stride = 94;
+  const index_type nt     = 2401;
+  scalar_type      results[stride];
+  buffer.seekg(0, std::ios::beg);
+
+  std::ostream  nullout(nullptr);
+  std::ostream& out = nullout;
+
+  // // Uncomment code below to print output to a file:
+  // std::ofstream fileout;
+  // fileout.open("example2_results.csv");
+  // std::ostream& out = fileout;
+
+  out << "Time,gen2speed,gen3speed,v2mag,v3mag\n";
+  out << 0. << "," << 1. << "," << 1. << "," << 1. << "," << 1. << "\n";
+
+  for (index_type i = 0; i < nt - 1; ++i)
+  {
+    for (index_type j = 0; j < stride; ++j)
+    {
+      buffer >> results[j];
+    }
+    real_type   t             = results[0];
+    scalar_type gen2speed     = 1 + results[7];
+    scalar_type gen2speed_ref = reference_solution[i + 1][1];
+    scalar_type gen3speed     = 1 + results[28];
+    scalar_type gen3speed_ref = reference_solution[i + 1][2];
+    scalar_type v2mag         = sqrt(results[2] * results[2] + results[3] * results[3]);
+    scalar_type v2mag_ref     = reference_solution[i + 1][4];
+    scalar_type v3mag         = sqrt(results[4] * results[4] + results[5] * results[5]);
+    scalar_type v3mag_ref     = reference_solution[i + 1][5];
+
+    out << t << "," << gen2speed << "," << gen3speed << "," << v2mag << "," << v3mag << "\n";
+
+    real_type err = std::max({std::abs(gen2speed - gen2speed_ref),
+                              std::abs(gen3speed - gen3speed_ref),
+                              std::abs(v2mag - v2mag_ref),
+                              std::abs(v3mag - v3mag_ref)});
+    if (err > worst_error)
+    {
+      worst_error      = err;
+      worst_error_time = t;
+    }
+  }
+  // fileout.close();
+
+  std::cout << "Worst error " << worst_error
+            << " at time t = " << worst_error_time << "\n";
+  std::cout << "\n\nComplete in " << (stop - start) / CLOCKS_PER_SEC << " seconds\n";
+
+  return worst_error < error_allowed ? 0 : 1;
+}