diff --git a/examples/PhasorDynamics/Example1/CMakeLists.txt b/examples/PhasorDynamics/Example1/CMakeLists.txt
index e79e6103b..d0e498a5e 100644
--- a/examples/PhasorDynamics/Example1/CMakeLists.txt
+++ b/examples/PhasorDynamics/Example1/CMakeLists.txt
@@ -4,6 +4,7 @@ target_link_libraries(phasordynamics_example1
 	GRIDKIT::phasor_dynamics_bus_fault 
 	GRIDKIT::phasor_dynamics_branch
 	GRIDKIT::phasor_dynamics_genrou
+	GRIDKIT::phasor_dynamics_load
 	GRIDKIT::solvers_dyn)
 install(TARGETS phasordynamics_example1 RUNTIME DESTINATION bin)
 
diff --git a/examples/PhasorDynamics/Example1/example1.cpp b/examples/PhasorDynamics/Example1/example1.cpp
index ac42fc236..04c78c5dd 100644
--- a/examples/PhasorDynamics/Example1/example1.cpp
+++ b/examples/PhasorDynamics/Example1/example1.cpp
@@ -23,6 +23,7 @@
 #include <Model/PhasorDynamics/SynchronousMachine/GENROUwS/Genrou.hpp>
 #include <Model/PhasorDynamics/SynchronousMachine/GENROUwS/GenrouData.hpp>
 #include <Model/PhasorDynamics/SystemModel.hpp>
+#include <Model/PhasorDynamics/SystemModelData.hpp>
 #include <Solver/Dynamic/Ida.hpp>
 #include <Utilities/Testing.hpp>
 
@@ -40,68 +41,72 @@ int main()
   //
   // Create model data
   //
-  BusData<real_type, index_type> bus_data_1;
-  bus_data_1.Vr0 = 0.9949877346411762;
-  bus_data_1.Vi0 = 0.09999703952427966;
-
-  BusData<real_type, index_type> bus_data_2;
-  bus_data_2.Vr0 = 1.0;
-  bus_data_2.Vi0 = 0.0;
-
-  BranchData<real_type, index_type> branch_data_1_2;
-  branch_data_1_2.R = 0.0;
-  branch_data_1_2.X = 0.1;
-  branch_data_1_2.G = 0.0;
-  branch_data_1_2.B = 0.0;
-
-  GenrouData<real_type, index_type> gen_data_1;
-  gen_data_1.unit_id = 1;
-  gen_data_1.p0      = 1.;
-  gen_data_1.q0      = 0.05013;
-  gen_data_1.H       = 3.;
-  gen_data_1.D       = 0.;
-  gen_data_1.Ra      = 0.;
-  gen_data_1.Tdop    = 7.;
-  gen_data_1.Tdopp   = .04;
-  gen_data_1.Tqopp   = .05;
-  gen_data_1.Tqop    = .75;
-  gen_data_1.Xd      = 2.1;
-  gen_data_1.Xdp     = 0.2;
-  gen_data_1.Xdpp    = 0.18;
-  gen_data_1.Xq      = 0.5;
-  gen_data_1.Xqp     = 0.5;
-  gen_data_1.Xqpp    = 0.18;
-  gen_data_1.Xl      = 0.15;
-  gen_data_1.S10     = 0.;
-  gen_data_1.S12     = 0.;
-
-  BusFaultData<real_type, index_type> bus_fault_data_1;
-  bus_fault_data_1.R      = 0.0;
-  bus_fault_data_1.X      = 1e-3;
-  bus_fault_data_1.status = false;
 
-  //
-  // Instantiate model components
-  //
-
-  Bus<scalar_type, size_t>         bus1(bus_data_1);
-  BusInfinite<scalar_type, size_t> bus2(bus_data_2);
-  Branch<scalar_type, size_t>      branch(&bus1, &bus2, branch_data_1_2);
-  BusFault<scalar_type, size_t>    fault(&bus1, bus_fault_data_1);
-  Genrou<scalar_type, size_t>      gen(&bus1, gen_data_1);
+  SystemModelData<scalar_type, index_type> data;
+
+  // Set bus data
+  data.bus.resize(2);
+
+  data.bus[0].bus_id   = 0;
+  data.bus[0].bus_type = BusData<scalar_type, index_type>::DEFAULT;
+  data.bus[0].Vr0      = 0.9949877346411762;
+  data.bus[0].Vi0      = 0.09999703952427966;
+
+  data.bus[1].bus_id   = 1;
+  data.bus[1].bus_type = BusData<scalar_type, index_type>::SLACK;
+  data.bus[1].Vr0      = 1.0;
+  data.bus[1].Vi0      = 0.0;
+
+  // Set branch data
+  data.branch.resize(1);
+
+  data.branch[0].bus1_id = data.bus[0].bus_id;
+  data.branch[0].bus2_id = data.bus[1].bus_id;
+  data.branch[0].R       = 0.0;
+  data.branch[0].X       = 0.1;
+  data.branch[0].G       = 0.0;
+  data.branch[0].B       = 0.0;
+
+  // Set generator data
+  data.genrou.resize(1);
+
+  data.genrou[0].unit_id = 1;
+  data.genrou[0].p0      = 1.;
+  data.genrou[0].q0      = 0.05013;
+  data.genrou[0].H       = 3.;
+  data.genrou[0].D       = 0.;
+  data.genrou[0].Ra      = 0.;
+  data.genrou[0].Tdop    = 7.;
+  data.genrou[0].Tdopp   = .04;
+  data.genrou[0].Tqopp   = .05;
+  data.genrou[0].Tqop    = .75;
+  data.genrou[0].Xd      = 2.1;
+  data.genrou[0].Xdp     = 0.2;
+  data.genrou[0].Xdpp    = 0.18;
+  data.genrou[0].Xq      = 0.5;
+  data.genrou[0].Xqp     = 0.5;
+  data.genrou[0].Xqpp    = 0.18;
+  data.genrou[0].Xl      = 0.15;
+  data.genrou[0].S10     = 0.;
+  data.genrou[0].S12     = 0.;
+
+  // Add faults
+  data.bus_fault.resize(1);
+
+  data.bus_fault[0].R      = 0.0;
+  data.bus_fault[0].X      = 1e-3;
+  data.bus_fault[0].status = false;
 
   //
-  // Create the 2-bus system
+  // Instantiate system model
   //
 
-  SystemModel<scalar_type, size_t> sys;
-  sys.addBus(&bus1);
-  sys.addBus(&bus2);
-  sys.addComponent(&branch);
-  sys.addComponent(&fault);
-  sys.addComponent(&gen);
+  SystemModel<scalar_type, index_type> sys(data);
   sys.allocate();
 
+  // Get access to the fault
+  auto* fault = sys.getBusFault(0);
+
   // Set time step to 1/4 of a 60Hz cycle
   real_type dt = 1.0 / 4.0 / 60.0;
 
@@ -153,13 +158,13 @@ int main()
   ida.runSimulation(1.0, nout, output_cb);
 
   // Introduce fault and run for the next 0.1s
-  fault.setStatus(true);
+  fault->setStatus(true);
   ida.initializeSimulation(1.0, false);
   nout = static_cast<int>(std::round((1.1 - 1.0) / dt));
   ida.runSimulation(1.1, nout, output_cb);
 
   // Clear the fault and run until t = 10s.
-  fault.setStatus(false);
+  fault->setStatus(false);
   ida.initializeSimulation(1.1, false);
   nout = static_cast<int>(std::round((10.0 - 1.1) / dt));
   ida.runSimulation(10.0, nout, output_cb);
diff --git a/examples/PhasorDynamics/Example2/example2.cpp b/examples/PhasorDynamics/Example2/example2.cpp
index 8a9cbc15f..3b7eeb68d 100644
--- a/examples/PhasorDynamics/Example2/example2.cpp
+++ b/examples/PhasorDynamics/Example2/example2.cpp
@@ -27,6 +27,7 @@
 #include <Model/PhasorDynamics/SynchronousMachine/GENROUwS/Genrou.hpp>
 #include <Model/PhasorDynamics/SynchronousMachine/GENROUwS/GenrouData.hpp>
 #include <Model/PhasorDynamics/SystemModel.hpp>
+#include <Model/PhasorDynamics/SystemModelData.hpp>
 #include <Solver/Dynamic/Ida.hpp>
 #include <Utilities/Testing.hpp>
 
@@ -88,131 +89,132 @@ int main()
   std::cout << "Example 2 version 1\n";
 
   //
-  // Create (load) model data
+  // Create model data
   //
 
+  SystemModelData<scalar_type, index_type> data;
+
+  // Set bus data
+  data.bus.resize(3);
+
+  // Bus 0
+  data.bus[0].bus_id   = 0;
+  data.bus[0].bus_type = BusData<scalar_type, index_type>::SLACK;
+  data.bus[0].Vr0      = 1.06;
+  data.bus[0].Vi0      = 0.0;
+
   // Bus 1
-  BusData<real_type, index_type> bus_data_1;
-  bus_data_1.Vr0 = 1.06;
-  bus_data_1.Vi0 = 0.0;
+  data.bus[1].bus_id   = 1;
+  data.bus[1].bus_type = BusData<scalar_type, index_type>::DEFAULT;
+  data.bus[1].Vr0      = 1.0599558398065716;
+  data.bus[1].Vi0      = -0.009675621941024773;
 
   // Bus 2
-  BusData<real_type, index_type> bus_data_2;
-  bus_data_2.Vr0 = 1.0599558398065716;
-  bus_data_2.Vi0 = -0.009675621941024773;
-
-  // Bus 3
-  BusData<real_type, index_type> bus_data_3;
-  bus_data_3.Vr0 = 0.9610827543495831;
-  bus_data_3.Vi0 = -0.13122476630506485;
+  data.bus[2].bus_id   = 2;
+  data.bus[2].bus_type = BusData<scalar_type, index_type>::DEFAULT;
+  data.bus[2].Vr0      = 0.9610827543495831;
+  data.bus[2].Vi0      = -0.13122476630506485;
+
+  // Set branch data
+  data.branch.resize(3);
+
+  // Branch 0-1
+  data.branch[0].bus1_id = data.bus[0].bus_id;
+  data.branch[0].bus2_id = data.bus[1].bus_id;
+  data.branch[0].R       = 0.05;
+  data.branch[0].X       = 0.21;
+  data.branch[0].G       = 0;
+  data.branch[0].B       = 0.1;
+
+  // Branch 0-2
+  data.branch[1].bus1_id = data.bus[0].bus_id;
+  data.branch[1].bus2_id = data.bus[2].bus_id;
+  data.branch[1].R       = 0.06;
+  data.branch[1].X       = 0.15;
+  data.branch[1].G       = 0;
+  data.branch[1].B       = 0.12;
 
   // Branch 1-2
-  BranchData<real_type, index_type> branch_data_1_2;
-  branch_data_1_2.R = 0.05;
-  branch_data_1_2.X = 0.21;
-  branch_data_1_2.G = 0;
-  branch_data_1_2.B = 0.1;
-
-  // Branch 1-3
-  BranchData<real_type, index_type> branch_data_1_3;
-  branch_data_1_3.R = 0.06;
-  branch_data_1_3.X = 0.15;
-  branch_data_1_3.G = 0;
-  branch_data_1_3.B = 0.12;
-
-  // Branch 2-3
-  BranchData<real_type, index_type> branch_data_2_3;
-  branch_data_2_3.R = 0.08;
-  branch_data_2_3.X = 0.27;
-  branch_data_2_3.G = 0;
-  branch_data_2_3.B = 0.45;
+  data.branch[2].bus1_id = data.bus[1].bus_id;
+  data.branch[2].bus2_id = data.bus[2].bus_id;
+  data.branch[2].R       = 0.08;
+  data.branch[2].X       = 0.27;
+  data.branch[2].G       = 0;
+  data.branch[2].B       = 0.45;
+
+  // Set generator data
+  data.genrou.resize(2);
+
+  // Generator on bus 1
+  data.genrou[0].bus_id  = 1;
+  data.genrou[0].unit_id = 0;
+  data.genrou[0].p0      = 0.5;
+  data.genrou[0].q0      = -0.07588;
+  data.genrou[0].H       = 2.7;
+  data.genrou[0].D       = 0.;
+  data.genrou[0].Ra      = 0.;
+  data.genrou[0].Tdop    = 7.;
+  data.genrou[0].Tdopp   = .04;
+  data.genrou[0].Tqopp   = .05;
+  data.genrou[0].Tqop    = .75;
+  data.genrou[0].Xd      = 1.9;
+  data.genrou[0].Xdp     = 0.17;
+  data.genrou[0].Xdpp    = 0.15;
+  data.genrou[0].Xq      = 0.4;
+  data.genrou[0].Xqp     = 0.35;
+  data.genrou[0].Xqpp    = 0.15;
+  data.genrou[0].Xl      = 0.14999;
+  data.genrou[0].S10     = 0.;
+  data.genrou[0].S12     = 0.;
 
   // Generator on bus 2
-  GenrouData<real_type, index_type> gen_data_2;
-  gen_data_2.unit_id = 1;
-  gen_data_2.p0      = 0.5;
-  gen_data_2.q0      = -0.07588;
-  gen_data_2.H       = 2.7;
-  gen_data_2.D       = 0.;
-  gen_data_2.Ra      = 0.;
-  gen_data_2.Tdop    = 7.;
-  gen_data_2.Tdopp   = .04;
-  gen_data_2.Tqopp   = .05;
-  gen_data_2.Tqop    = .75;
-  gen_data_2.Xd      = 1.9;
-  gen_data_2.Xdp     = 0.17;
-  gen_data_2.Xdpp    = 0.15;
-  gen_data_2.Xq      = 0.4;
-  gen_data_2.Xqp     = 0.35;
-  gen_data_2.Xqpp    = 0.15;
-  gen_data_2.Xl      = 0.14999;
-  gen_data_2.S10     = 0.;
-  gen_data_2.S12     = 0.;
-
-  // Generator on bus 3
-  GenrouData<real_type, index_type> gen_data_3;
-  gen_data_3.unit_id = 1;
-  gen_data_3.p0      = 0.25;
-  gen_data_3.q0      = 0.26587;
-  gen_data_3.H       = 1.6;
-  gen_data_3.D       = 0.;
-  gen_data_3.Ra      = 0.;
-  gen_data_3.Tdop    = 7.5;
-  gen_data_3.Tdopp   = .04;
-  gen_data_3.Tqopp   = .05;
-  gen_data_3.Tqop    = .75;
-  gen_data_3.Xd      = 2.3;
-  gen_data_3.Xdp     = 0.2;
-  gen_data_3.Xdpp    = 0.18;
-  gen_data_3.Xq      = 0.5;
-  gen_data_3.Xqp     = 0.5;
-  gen_data_3.Xqpp    = 0.18;
-  gen_data_3.Xl      = 0.15;
-  gen_data_3.S10     = 0.;
-  gen_data_3.S12     = 0.;
-
-  // Load on bus 3
-  LoadData<real_type, index_type> load_data_3;
-  load_data_3.R      = 0.4447197839297772;
-  load_data_3.X      = 0.20330047265361242;
-  load_data_3.bus_id = 3;
-
-  BusFaultData<real_type, index_type> bus_fault_data_3;
-  bus_fault_data_3.R      = 0.0;
-  bus_fault_data_3.X      = 1e-5;
-  bus_fault_data_3.status = false;
+  data.genrou[1].bus_id  = 2;
+  data.genrou[1].unit_id = 1;
+  data.genrou[1].p0      = 0.25;
+  data.genrou[1].q0      = 0.26587;
+  data.genrou[1].H       = 1.6;
+  data.genrou[1].D       = 0.;
+  data.genrou[1].Ra      = 0.;
+  data.genrou[1].Tdop    = 7.5;
+  data.genrou[1].Tdopp   = .04;
+  data.genrou[1].Tqopp   = .05;
+  data.genrou[1].Tqop    = .75;
+  data.genrou[1].Xd      = 2.3;
+  data.genrou[1].Xdp     = 0.2;
+  data.genrou[1].Xdpp    = 0.18;
+  data.genrou[1].Xq      = 0.5;
+  data.genrou[1].Xqp     = 0.5;
+  data.genrou[1].Xqpp    = 0.18;
+  data.genrou[1].Xl      = 0.15;
+  data.genrou[1].S10     = 0.;
+  data.genrou[1].S12     = 0.;
+
+  // Set load data
+  data.load.resize(1);
+
+  // Load on bus 2
+  data.load[0].bus_id = 2;
+  data.load[0].R      = 0.4447197839297772;
+  data.load[0].X      = 0.20330047265361242;
+
+  // Set fault data
+  data.bus_fault.resize(1);
+
+  data.bus_fault[0].bus_id = 2;
+  data.bus_fault[0].R      = 0.0;
+  data.bus_fault[0].X      = 1e-5;
+  data.bus_fault[0].status = false;
 
   //
-  // Instantiate model components
+  // Instantiate system
   //
 
-  BusInfinite<scalar_type, index_type> bus1(bus_data_1);
-  Bus<scalar_type, index_type>         bus2(bus_data_2);
-  Bus<scalar_type, index_type>         bus3(bus_data_3);
-
-  Branch<scalar_type, index_type> branch12(&bus1, &bus2, branch_data_1_2);
-  Branch<scalar_type, index_type> branch13(&bus1, &bus3, branch_data_1_3);
-  Branch<scalar_type, index_type> branch23(&bus2, &bus3, branch_data_2_3);
-
-  Genrou<scalar_type, index_type>   gen2(&bus2, gen_data_2);
-  Genrou<scalar_type, index_type>   gen3(&bus3, gen_data_3);
-  Load<scalar_type, index_type>     load3(&bus3, load_data_3);
-  BusFault<scalar_type, index_type> fault(&bus3, bus_fault_data_3);
-
-  /* Connect everything together */
-  SystemModel<scalar_type, index_type> sys;
-  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);
+  SystemModel<scalar_type, index_type> sys(data);
   sys.allocate();
 
+  // Get access to fault 0
+  auto* fault = sys.getBusFault(0);
+
   real_type dt = 1.0 / 4.0 / 60.0;
 
   std::vector<OutputData> output;
@@ -241,13 +243,13 @@ int main()
   ida.runSimulation(1.0, nout, output_cb);
 
   // Introduce fault to ground and run for 0.1s
-  fault.setStatus(true);
+  fault->setStatus(true);
   ida.initializeSimulation(1.0, false);
   nout = static_cast<int>(std::round((1.1 - 1.0) / dt));
   ida.runSimulation(1.1, nout, output_cb);
 
   // Clear fault and run until t = 10s.
-  fault.setStatus(false);
+  fault->setStatus(false);
   ida.initializeSimulation(1.1, false);
   nout = static_cast<int>(std::round((10.0 - 1.1) / dt));
   ida.runSimulation(10.0, nout, output_cb);
diff --git a/src/Model/PhasorDynamics/Branch/Branch.cpp b/src/Model/PhasorDynamics/Branch/Branch.cpp
index 5132a7a8d..effc8bd5c 100644
--- a/src/Model/PhasorDynamics/Branch/Branch.cpp
+++ b/src/Model/PhasorDynamics/Branch/Branch.cpp
@@ -73,7 +73,7 @@ namespace GridKit
     }
 
     template <class ScalarT, typename IdxT>
-    Branch<ScalarT, IdxT>::Branch(bus_type* bus1, bus_type* bus2, model_data_type& data)
+    Branch<ScalarT, IdxT>::Branch(bus_type* bus1, bus_type* bus2, const model_data_type& data)
       : bus1_(bus1),
         bus2_(bus2),
         R_(data.R),
diff --git a/src/Model/PhasorDynamics/Branch/Branch.hpp b/src/Model/PhasorDynamics/Branch/Branch.hpp
index cf3cba0e1..caefa491f 100644
--- a/src/Model/PhasorDynamics/Branch/Branch.hpp
+++ b/src/Model/PhasorDynamics/Branch/Branch.hpp
@@ -53,7 +53,7 @@ namespace GridKit
     public:
       Branch(bus_type* bus1, bus_type* bus2);
       Branch(bus_type* bus1, bus_type* bus2, real_type R, real_type X, real_type G, real_type B);
-      Branch(bus_type* bus1, bus_type* bus2, model_data_type& data);
+      Branch(bus_type* bus1, bus_type* bus2, const model_data_type& data);
       virtual ~Branch();
 
       virtual int allocate() override;
diff --git a/src/Model/PhasorDynamics/Bus/Bus.cpp b/src/Model/PhasorDynamics/Bus/Bus.cpp
index d8a1a3481..03473c8dd 100644
--- a/src/Model/PhasorDynamics/Bus/Bus.cpp
+++ b/src/Model/PhasorDynamics/Bus/Bus.cpp
@@ -61,7 +61,7 @@ namespace GridKit
      * @param[in] data - structure with bus data
      */
     template <class ScalarT, typename IdxT>
-    Bus<ScalarT, IdxT>::Bus(DataT& data)
+    Bus<ScalarT, IdxT>::Bus(const DataT& data)
       : BusBase<ScalarT, IdxT>(data.bus_id),
         Vr0_(data.Vr0),
         Vi0_(data.Vi0)
diff --git a/src/Model/PhasorDynamics/Bus/Bus.hpp b/src/Model/PhasorDynamics/Bus/Bus.hpp
index 57828ef01..11dda1be8 100644
--- a/src/Model/PhasorDynamics/Bus/Bus.hpp
+++ b/src/Model/PhasorDynamics/Bus/Bus.hpp
@@ -43,7 +43,7 @@ namespace GridKit
 
       Bus();
       Bus(ScalarT Vr, ScalarT Vi);
-      Bus(DataT& data);
+      Bus(const DataT& data);
       virtual ~Bus();
 
       virtual int allocate() override;
@@ -58,7 +58,7 @@ namespace GridKit
 
       virtual int BusType() const override
       {
-        return BusBase<ScalarT, IdxT>::BusType::DEFAULT;
+        return BusData<real_type, IdxT>::BusType::DEFAULT;
       }
 
       virtual ScalarT& Vr() override
diff --git a/src/Model/PhasorDynamics/Bus/BusData.hpp b/src/Model/PhasorDynamics/Bus/BusData.hpp
index abd8ab595..115deb5dd 100644
--- a/src/Model/PhasorDynamics/Bus/BusData.hpp
+++ b/src/Model/PhasorDynamics/Bus/BusData.hpp
@@ -23,10 +23,18 @@ namespace GridKit
     template <typename RealT, typename IdxT>
     struct BusData
     {
+      enum BusType
+      {
+        INVALID = 0,
+        DEFAULT,
+        SLACK
+      };
+
       RealT Vr0{0.0}; ///< Initial value for real bus voltage
       RealT Vi0{0.0}; ///< Initial value for imaginary bus voltage
 
-      IdxT bus_id{0}; ///< Unique ID of bus 1
+      IdxT    bus_id{0}; ///< Unique ID of bus 1
+      BusType bus_type{INVALID};
     };
   } // namespace PhasorDynamics
 } // namespace GridKit
diff --git a/src/Model/PhasorDynamics/Bus/BusFactory.hpp b/src/Model/PhasorDynamics/Bus/BusFactory.hpp
new file mode 100644
index 000000000..34eae481e
--- /dev/null
+++ b/src/Model/PhasorDynamics/Bus/BusFactory.hpp
@@ -0,0 +1,41 @@
+
+#pragma once
+
+#include <Model/PhasorDynamics/Bus/Bus.hpp>
+#include <Model/PhasorDynamics/Bus/BusData.hpp>
+#include <Model/PhasorDynamics/Bus/BusInfinite.hpp>
+
+namespace GridKit
+{
+  namespace PhasorDynamics
+  {
+    template <typename ScalarT = double, typename IdxT = int>
+    class BusFactory
+    {
+    public:
+      using real_type = typename Model::Evaluator<ScalarT, IdxT>::real_type;
+      using BusData   = GridKit::PhasorDynamics::BusData<real_type, IdxT>;
+
+      BusFactory() = delete;
+
+      static BusBase<ScalarT, IdxT>* create(const BusData& data)
+      {
+        BusBase<ScalarT, IdxT>* bus = nullptr;
+
+        switch (data.bus_type)
+        {
+        case BusData::DEFAULT:
+          bus = new Bus<ScalarT, IdxT>(data);
+          break;
+        case BusData::SLACK:
+          bus = new BusInfinite<ScalarT, IdxT>(data);
+          break;
+        default:
+          // Throw exception
+          std::cout << "Bus type " << data.bus_type << " unrecognized.\n";
+        }
+        return bus;
+      }
+    };
+  } // namespace PhasorDynamics
+} // namespace GridKit
diff --git a/src/Model/PhasorDynamics/Bus/BusInfinite.cpp b/src/Model/PhasorDynamics/Bus/BusInfinite.cpp
index 6f4378d65..f1e96a339 100644
--- a/src/Model/PhasorDynamics/Bus/BusInfinite.cpp
+++ b/src/Model/PhasorDynamics/Bus/BusInfinite.cpp
@@ -66,7 +66,7 @@ namespace GridKit
      * @param[in] data - structure with bus data
      */
     template <class ScalarT, typename IdxT>
-    BusInfinite<ScalarT, IdxT>::BusInfinite(DataT& data)
+    BusInfinite<ScalarT, IdxT>::BusInfinite(const DataT& data)
       : BusBase<ScalarT, IdxT>(data.bus_id),
         Vr_(data.Vr0),
         Vi_(data.Vi0)
diff --git a/src/Model/PhasorDynamics/Bus/BusInfinite.hpp b/src/Model/PhasorDynamics/Bus/BusInfinite.hpp
index 0657ec374..b7247401f 100644
--- a/src/Model/PhasorDynamics/Bus/BusInfinite.hpp
+++ b/src/Model/PhasorDynamics/Bus/BusInfinite.hpp
@@ -41,7 +41,7 @@ namespace GridKit
 
       BusInfinite();
       BusInfinite(ScalarT Vr, ScalarT Vi);
-      BusInfinite(DataT& data);
+      BusInfinite(const DataT& data);
       virtual ~BusInfinite();
 
       virtual int allocate() override;
@@ -56,7 +56,7 @@ namespace GridKit
 
       virtual int BusType() const override
       {
-        return BusBase<ScalarT, IdxT>::BusType::SLACK;
+        return BusData<real_type, IdxT>::BusType::SLACK;
       }
 
       virtual ScalarT& Vr() override
diff --git a/src/Model/PhasorDynamics/BusBase.hpp b/src/Model/PhasorDynamics/BusBase.hpp
index 332aaaed1..322e6f095 100644
--- a/src/Model/PhasorDynamics/BusBase.hpp
+++ b/src/Model/PhasorDynamics/BusBase.hpp
@@ -19,12 +19,6 @@ namespace GridKit
     public:
       using real_type = typename Model::Evaluator<ScalarT, IdxT>::real_type;
 
-      enum BusType
-      {
-        DEFAULT = 1,
-        SLACK
-      };
-
       BusBase()
         : size_(0),
           size_quad_(0),
diff --git a/src/Model/PhasorDynamics/BusFault/BusFault.cpp b/src/Model/PhasorDynamics/BusFault/BusFault.cpp
index 727997721..597f20cbf 100644
--- a/src/Model/PhasorDynamics/BusFault/BusFault.cpp
+++ b/src/Model/PhasorDynamics/BusFault/BusFault.cpp
@@ -64,7 +64,7 @@ namespace GridKit
      * @param bus2 - pointer to bus-2
      */
     template <class ScalarT, typename IdxT>
-    BusFault<ScalarT, IdxT>::BusFault(bus_type* bus, DataT& data)
+    BusFault<ScalarT, IdxT>::BusFault(bus_type* bus, const DataT& data)
       : bus_(bus),
         R_(data.R),
         X_(data.X),
diff --git a/src/Model/PhasorDynamics/BusFault/BusFault.hpp b/src/Model/PhasorDynamics/BusFault/BusFault.hpp
index 8963a7b30..72be0717f 100644
--- a/src/Model/PhasorDynamics/BusFault/BusFault.hpp
+++ b/src/Model/PhasorDynamics/BusFault/BusFault.hpp
@@ -43,7 +43,7 @@ namespace GridKit
     public:
       BusFault(bus_type* bus);
       BusFault(bus_type* bus, real_type R, real_type X, int status);
-      BusFault(bus_type* bus, DataT& data);
+      BusFault(bus_type* bus, const DataT& data);
       ~BusFault() = default;
 
       int allocate() override;
diff --git a/src/Model/PhasorDynamics/INPUT_FORMAT.md b/src/Model/PhasorDynamics/INPUT_FORMAT.md
index 149048b6c..0fdfdcfb4 100644
--- a/src/Model/PhasorDynamics/INPUT_FORMAT.md
+++ b/src/Model/PhasorDynamics/INPUT_FORMAT.md
@@ -49,7 +49,7 @@ Contained in the `header` key is an object with the following items:
   `format_version`   | Non-negative integer indicating the format version
   `format_revision`  | Non-negative integer indicating the format revision
   `case_name`        | A string containing the name of the case
-  `case_datetime`    | Optional string in the ISO 8601 format indicating a datetime associated with the case. Including a time is optional, but if one is specified, it is recommended that a UTC offset be included.
+  `case_date_time`   | Optional string in the ISO 8601 format indicating a datetime associated with the case. Including a time is optional, but if one is specified, it is recommended that a UTC offset be included.
   `case_description` | A string with more specific description of what is modeled in the case
   `case_comments`    | A string with additional notes as needed
   `freq_base`        | A floating point value indicating the system frequency base in hertz (Hz). This is commonly 60 Hz
diff --git a/src/Model/PhasorDynamics/Load/Load.cpp b/src/Model/PhasorDynamics/Load/Load.cpp
index bbc187904..8f130cbae 100644
--- a/src/Model/PhasorDynamics/Load/Load.cpp
+++ b/src/Model/PhasorDynamics/Load/Load.cpp
@@ -1,29 +1,29 @@
-
-#include "LoadImpl.hpp"
-
-namespace GridKit
-{
-  namespace PhasorDynamics
-  {
-    /**
-     * @brief Jacobian evaluation not implemented
-     *
-     * @tparam ScalarT - scalar data type
-     * @tparam IdxT    - matrix index data type
-     * @return int - error code, 0 = success
-     */
-    template <class ScalarT, typename IdxT>
-    int Load<ScalarT, IdxT>::evaluateJacobian()
-    {
-      std::cout << "Evaluate Jacobian for Load..." << std::endl;
-      std::cout << "Jacobian evaluation is not implemented!" << std::endl;
-
-      return 0;
-    }
-
-    // Available template instantiations
-    template class Load<double, long int>;
-    template class Load<double, size_t>;
-
-  } // namespace PhasorDynamics
-} // namespace GridKit
+
+#include "LoadImpl.hpp"
+
+namespace GridKit
+{
+  namespace PhasorDynamics
+  {
+    /**
+     * @brief Jacobian evaluation not implemented
+     *
+     * @tparam ScalarT - scalar data type
+     * @tparam IdxT    - matrix index data type
+     * @return int - error code, 0 = success
+     */
+    template <class ScalarT, typename IdxT>
+    int Load<ScalarT, IdxT>::evaluateJacobian()
+    {
+      std::cout << "Evaluate Jacobian for Load..." << std::endl;
+      std::cout << "Jacobian evaluation is not implemented!" << std::endl;
+
+      return 0;
+    }
+
+    // Available template instantiations
+    template class Load<double, long int>;
+    template class Load<double, size_t>;
+
+  } // namespace PhasorDynamics
+} // namespace GridKit
diff --git a/src/Model/PhasorDynamics/Load/Load.hpp b/src/Model/PhasorDynamics/Load/Load.hpp
index ee20f428b..928898c9a 100644
--- a/src/Model/PhasorDynamics/Load/Load.hpp
+++ b/src/Model/PhasorDynamics/Load/Load.hpp
@@ -47,7 +47,7 @@ namespace GridKit
       Load(bus_type* bus);
       Load(bus_type* bus, real_type R, real_type X);
       Load(bus_type* bus, IdxT component_id);
-      Load(bus_type* bus, model_data_type& data);
+      Load(bus_type* bus, const model_data_type& data);
       virtual ~Load();
 
       virtual int allocate() override;
diff --git a/src/Model/PhasorDynamics/Load/LoadImpl.hpp b/src/Model/PhasorDynamics/Load/LoadImpl.hpp
index e8cfbd2ba..59e46668e 100644
--- a/src/Model/PhasorDynamics/Load/LoadImpl.hpp
+++ b/src/Model/PhasorDynamics/Load/LoadImpl.hpp
@@ -38,8 +38,8 @@ namespace GridKit
     }
 
     template <class ScalarT, typename IdxT>
-    Load<ScalarT, IdxT>::Load(bus_type*        bus,
-                              model_data_type& data)
+    Load<ScalarT, IdxT>::Load(bus_type*              bus,
+                              const model_data_type& data)
       : bus_(bus),
         R_(data.R),
         X_(data.X)
diff --git a/src/Model/PhasorDynamics/SynchronousMachine/GENROUwS/Genrou.cpp b/src/Model/PhasorDynamics/SynchronousMachine/GENROUwS/Genrou.cpp
index 55895d50b..31b7982db 100644
--- a/src/Model/PhasorDynamics/SynchronousMachine/GENROUwS/Genrou.cpp
+++ b/src/Model/PhasorDynamics/SynchronousMachine/GENROUwS/Genrou.cpp
@@ -119,7 +119,7 @@ namespace GridKit
      *
      */
     template <class ScalarT, typename IdxT>
-    Genrou<ScalarT, IdxT>::Genrou(bus_type* bus, model_data_type& data)
+    Genrou<ScalarT, IdxT>::Genrou(bus_type* bus, const model_data_type& data)
       : bus_(bus),
         busID_(0),
         unit_id_(data.unit_id),
diff --git a/src/Model/PhasorDynamics/SynchronousMachine/GENROUwS/Genrou.hpp b/src/Model/PhasorDynamics/SynchronousMachine/GENROUwS/Genrou.hpp
index 283840840..305d47966 100644
--- a/src/Model/PhasorDynamics/SynchronousMachine/GENROUwS/Genrou.hpp
+++ b/src/Model/PhasorDynamics/SynchronousMachine/GENROUwS/Genrou.hpp
@@ -52,7 +52,7 @@ namespace GridKit
 
     public:
       Genrou(bus_type* bus, IdxT unit_id);
-      Genrou(bus_type* bus, model_data_type& data);
+      Genrou(bus_type* bus, const model_data_type& data);
       Genrou(bus_type* bus,
              IdxT      unit_id,
              ScalarT   p0,
diff --git a/src/Model/PhasorDynamics/SystemModel.hpp b/src/Model/PhasorDynamics/SystemModel.hpp
index 03347b224..848e87f14 100644
--- a/src/Model/PhasorDynamics/SystemModel.hpp
+++ b/src/Model/PhasorDynamics/SystemModel.hpp
@@ -6,8 +6,16 @@
 
 #include <Model/PhasorDynamics/BusBase.hpp>
 #include <Model/PhasorDynamics/Component.hpp>
+#include <Model/PhasorDynamics/SystemModelData.hpp>
 #include <ScalarTraits.hpp>
 
+// Temporary
+#include <Model/PhasorDynamics/Branch/Branch.hpp>
+#include <Model/PhasorDynamics/Bus/BusFactory.hpp>
+#include <Model/PhasorDynamics/BusFault/BusFault.hpp>
+#include <Model/PhasorDynamics/Load/Load.hpp>
+#include <Model/PhasorDynamics/SynchronousMachine/GENROUwS/Genrou.hpp>
+
 namespace GridKit
 {
   namespace PhasorDynamics
@@ -63,11 +71,82 @@ namespace GridKit
         this->max_steps_ = 2000;
       }
 
+      /**
+       * @brief Construct a new System Model object
+       *
+       * @param[in] data - Data structure with complete system data
+       *
+       * @pre SystemModelData contains consistent connectivity information
+       * and physically meaningful model parameters.
+       *
+       * @post All component models in SystemModelData are created, and
+       * correctly connected into the system model.
+       */
+      SystemModel(SystemModelData<real_type, IdxT>& data)
+      {
+        // Set system model tolerances
+        rel_tol_         = 1e-7;
+        abs_tol_         = 1e-9;
+        this->max_steps_ = 2000;
+
+        owns_components_ = true;
+
+        // Add electrical buses
+        for (const auto& busdata : data.bus)
+        {
+          BusBase<ScalarT, IdxT>* bus = BusFactory<ScalarT, IdxT>::create(busdata);
+          addBus(bus);
+        }
+
+        // Add branches
+        for (const auto& branchdata : data.branch)
+        {
+          auto* branch = new Branch<ScalarT, IdxT>(getBus(branchdata.bus1_id), getBus(branchdata.bus2_id), branchdata);
+          addComponent(branch);
+        }
+
+        // Add loads
+        for (const auto& loaddata : data.load)
+        {
+          auto* load = new Load<ScalarT, IdxT>(getBus(loaddata.bus_id), loaddata);
+          addComponent(load);
+        }
+
+        // Add generators
+        for (const auto& gendata : data.genrou)
+        {
+          auto* gen = new Genrou<ScalarT, IdxT>(getBus(gendata.bus_id), gendata);
+          addComponent(gen);
+        }
+
+        // Add faults
+        for (const auto& faultdata : data.bus_fault)
+        {
+          auto* fault = new BusFault<ScalarT, IdxT>(getBus(faultdata.bus_id), faultdata);
+          addFault(fault);
+        }
+      }
+
       /**
        * @brief Destructor for the system model
+       *
+       * If the SystemModel owns the components, it needs to delete them upon
+       * destructor call.
        */
       virtual ~SystemModel()
       {
+        if (owns_components_)
+        {
+          for (auto component : components_)
+          {
+            delete component;
+          }
+
+          for (auto bus : buses_)
+          {
+            delete bus;
+          }
+        }
       }
 
       /**
@@ -648,9 +727,39 @@ namespace GridKit
         components_.push_back(component);
       }
 
+      void addFault(component_type* component)
+      {
+        components_.push_back(component);
+        faults_.push_back(component);
+      }
+
+      bus_type* getBus(IdxT busid)
+      {
+        // Need to implement mapping of bus IDs to buses in the system model
+        assert((buses_[busid])->BusID() == busid);
+        return buses_[busid];
+      }
+
+      /**
+       * @brief Return pointer to a bus fault model
+       *
+       * This function is used to provide easier access to setting and
+       * clearing faults from the SystemModel interface.
+       *
+       * @warning This is a hack to get access to bus faults in examples.
+       * A more comprehensive solution is needed.
+       */
+      BusFault<ScalarT, IdxT>* getBusFault(IdxT fault_id)
+      {
+        return dynamic_cast<BusFault<ScalarT, IdxT>*>(faults_[fault_id]);
+      }
+
     private:
       std::vector<bus_type*>       buses_;
       std::vector<component_type*> components_;
+      std::vector<component_type*> faults_;
+
+      bool owns_components_{false};
 
     }; // class SystemModel
 
diff --git a/src/Model/PhasorDynamics/SystemModelData.hpp b/src/Model/PhasorDynamics/SystemModelData.hpp
new file mode 100644
index 000000000..28696d0af
--- /dev/null
+++ b/src/Model/PhasorDynamics/SystemModelData.hpp
@@ -0,0 +1,42 @@
+#pragma once
+
+#include <Model/PhasorDynamics/Branch/BranchData.hpp>
+#include <Model/PhasorDynamics/Bus/BusData.hpp>
+#include <Model/PhasorDynamics/BusFault/BusFaultData.hpp>
+#include <Model/PhasorDynamics/Load/LoadData.hpp>
+#include <Model/PhasorDynamics/SynchronousMachine/GENROUwS/GenrouData.hpp>
+
+namespace GridKit
+{
+  namespace PhasorDynamics
+  {
+    template <typename RealT = double, typename IdxT = size_t>
+    struct SystemModelData
+    {
+      using BranchDataT   = BranchData<RealT, IdxT>;
+      using BusDataT      = BusData<RealT, IdxT>;
+      using BusFaultDataT = BusFaultData<RealT, IdxT>;
+      using GenrouDataT   = GenrouData<RealT, IdxT>;
+      using LoadDataT     = LoadData<RealT, IdxT>;
+
+      // Header info
+      unsigned short format_version;
+      unsigned short format_revision;
+      std::string    case_name;
+      std::string    case_date_time;
+      std::string    case_description;
+      std::string    case_comments;
+      RealT          freq_base;
+      RealT          va_base;
+
+      // Bus data
+      std::vector<BusDataT> bus;
+
+      // Component data
+      std::vector<BranchDataT>   branch;
+      std::vector<BusFaultDataT> bus_fault;
+      std::vector<GenrouDataT>   genrou;
+      std::vector<LoadDataT>     load;
+    };
+  } // namespace PhasorDynamics
+} // namespace GridKit
diff --git a/tests/UnitTests/PhasorDynamics/CMakeLists.txt b/tests/UnitTests/PhasorDynamics/CMakeLists.txt
index c8b0a5368..f461ce686 100644
--- a/tests/UnitTests/PhasorDynamics/CMakeLists.txt
+++ b/tests/UnitTests/PhasorDynamics/CMakeLists.txt
@@ -28,6 +28,8 @@ target_link_libraries(test_phasor_gen_classical GRIDKIT::phasor_dynamics_gen_cla
 add_executable(test_phasor_system runSystemTests.cpp)
 target_link_libraries(test_phasor_system GRIDKIT::phasor_dynamics_load
                                          GRIDKIT::phasor_dynamics_branch
+                                         GRIDKIT::phasor_dynamics_genrou
+                                         GRIDKIT::phasor_dynamics_bus_fault
                                          GRIDKIT::phasor_dynamics_bus)
 
 add_test(NAME PhasorDynamicsBusTest COMMAND $<TARGET_FILE:test_phasor_bus>)
diff --git a/tests/UnitTests/PhasorDynamics/SystemTests.hpp b/tests/UnitTests/PhasorDynamics/SystemTests.hpp
index ea1cf7f29..dc8ecf81e 100644
--- a/tests/UnitTests/PhasorDynamics/SystemTests.hpp
+++ b/tests/UnitTests/PhasorDynamics/SystemTests.hpp
@@ -7,6 +7,7 @@
 #include <Model/PhasorDynamics/Bus/Bus.hpp>
 #include <Model/PhasorDynamics/Bus/BusInfinite.hpp>
 #include <Model/PhasorDynamics/SystemModel.hpp>
+#include <Model/PhasorDynamics/SystemModelData.hpp>
 #include <Utilities/TestHelpers.hpp>
 #include <Utilities/Testing.hpp>
 
@@ -29,21 +30,72 @@ namespace GridKit
       {
         TestStatus success = true;
 
-        // ScalarT Vr{1.0};
-        // ScalarT Vi{2.0};
-
         PhasorDynamics::SystemModel<ScalarT, IdxT>* system = nullptr;
 
         // Create an empty system
         system = new PhasorDynamics::SystemModel<ScalarT, IdxT>();
 
-        success *= (system != nullptr);
-
-        if (system)
+        if (system == nullptr)
         {
-          delete system;
+          std::cout << "Default constructor failed!\n";
+          success = false;
+          return success.report(__func__);
         }
 
+        delete system;
+        system = nullptr;
+
+        PhasorDynamics::SystemModelData<ScalarT, IdxT> data;
+
+        // Set bus data
+        data.bus.resize(2);
+
+        // Bus 0
+        data.bus[0].bus_id   = 0;
+        data.bus[0].bus_type = PhasorDynamics::BusData<ScalarT, IdxT>::SLACK;
+        data.bus[0].Vr0      = 10.0;
+        data.bus[0].Vi0      = 20.0;
+
+        // Bus 1
+        data.bus[1].bus_id   = 1;
+        data.bus[1].bus_type = PhasorDynamics::BusData<ScalarT, IdxT>::SLACK;
+        data.bus[1].Vr0      = 30.0;
+        data.bus[1].Vi0      = 40.0;
+
+        // Set branch data
+        data.branch.resize(1);
+
+        // Branch 0-1
+        data.branch[0].bus1_id = data.bus[0].bus_id;
+        data.branch[0].bus2_id = data.bus[1].bus_id;
+        data.branch[0].R       = 2.0;
+        data.branch[0].X       = 4.0;
+        data.branch[0].G       = 0.2;
+        data.branch[0].B       = 1.2;
+
+        // Create an empty system model
+        system = new PhasorDynamics::SystemModel<ScalarT, IdxT>(data);
+        system->allocate();
+        system->initialize();
+        system->evaluateResidual();
+
+        // Answer keys
+        const ScalarT Ir0{17.0};  ///< Solution: real current entering bus-0
+        const ScalarT Ii0{-10.0}; ///< Solution: imaginary current entering bus-0
+        const ScalarT Ir1{15.0};  ///< Solution: real current entering bus-1
+        const ScalarT Ii1{-20.0}; ///< Solution: imaginary current entering bus-1
+
+        auto* bus0 = system->getBus(0);
+        auto* bus1 = system->getBus(1);
+
+        success *= isEqual(bus0->Ir(), Ir0);
+        success *= isEqual(bus0->Ii(), Ii0);
+        success *= isEqual(bus1->Ir(), Ir1);
+        success *= isEqual(bus1->Ii(), Ii1);
+
+        delete system;
+        system = nullptr;
+
         return success.report(__func__);
       }