1414#include < iostream>
1515
1616#include < Model/PhasorDynamics/Branch/Branch.hpp>
17+ #include < Model/PhasorDynamics/Branch/BranchData.hpp>
1718#include < Model/PhasorDynamics/Bus/Bus.hpp>
19+ #include < Model/PhasorDynamics/Bus/BusData.hpp>
1820#include < Model/PhasorDynamics/Bus/BusInfinite.hpp>
1921#include < Model/PhasorDynamics/BusFault/BusFault.hpp>
22+ #include < Model/PhasorDynamics/BusFault/BusFaultData.hpp>
2023#include < Model/PhasorDynamics/SynchronousMachine/GENROUwS/Genrou.hpp>
24+ #include < Model/PhasorDynamics/SynchronousMachine/GENROUwS/GenrouData.hpp>
2125#include < Model/PhasorDynamics/SystemModel.hpp>
2226#include < Solver/Dynamic/Ida.hpp>
2327#include < Utilities/Testing.hpp>
@@ -27,45 +31,79 @@ int main()
2731 using namespace GridKit ::PhasorDynamics;
2832 using namespace AnalysisManager ::Sundials;
2933
34+ using scalar_type = double ;
35+ using real_type = double ;
36+ using index_type = size_t ;
37+
3038 std::cout << " Example 1 version 2\n " ;
3139
32- /* Create model parts */
33- SystemModel<double , size_t > sys;
34- Bus<double , size_t > bus1 (0.9949877346411762 , 0.09999703952427966 );
35- BusInfinite<double , size_t > bus2 (1.0 , 0.0 );
36- Branch<double , size_t > branch (&bus1, &bus2, 0 , 0.1 , 0 , 0 );
37- BusFault<double , size_t > fault (&bus1, 0 , 1e-3 , 0 );
38-
39- Genrou<double , size_t > gen (&bus1,
40- 1 ,
41- 1 .,
42- 0.05013 ,
43- 3 .,
44- 0 .,
45- 0 .,
46- 7 .,
47- .04 ,
48- .05 ,
49- .75 ,
50- 2.1 ,
51- 0.2 ,
52- 0.18 ,
53- 0.5 ,
54- 0.5 ,
55- 0.18 ,
56- 0.15 ,
57- 0 .,
58- 0 .);
59-
60- /* Connect everything together */
40+ //
41+ // Create model data
42+ //
43+ BusData<real_type, index_type> bus_data_1;
44+ bus_data_1.Vr0 = 0.9949877346411762 ;
45+ bus_data_1.Vi0 = 0.09999703952427966 ;
46+
47+ BusData<real_type, index_type> bus_data_2;
48+ bus_data_2.Vr0 = 1.0 ;
49+ bus_data_2.Vi0 = 0.0 ;
50+
51+ BranchData<real_type, index_type> branch_data_1_2;
52+ branch_data_1_2.R = 0.0 ;
53+ branch_data_1_2.X = 0.1 ;
54+ branch_data_1_2.G = 0.0 ;
55+ branch_data_1_2.B = 0.0 ;
56+
57+ GenrouData<real_type, index_type> gen_data_1;
58+ gen_data_1.unit_id = 1 ;
59+ gen_data_1.p0 = 1 .;
60+ gen_data_1.q0 = 0.05013 ;
61+ gen_data_1.H = 3 .;
62+ gen_data_1.D = 0 .;
63+ gen_data_1.Ra = 0 .;
64+ gen_data_1.Tdop = 7 .;
65+ gen_data_1.Tdopp = .04 ;
66+ gen_data_1.Tqopp = .05 ;
67+ gen_data_1.Tqop = .75 ;
68+ gen_data_1.Xd = 2.1 ;
69+ gen_data_1.Xdp = 0.2 ;
70+ gen_data_1.Xdpp = 0.18 ;
71+ gen_data_1.Xq = 0.5 ;
72+ gen_data_1.Xqp = 0.5 ;
73+ gen_data_1.Xqpp = 0.18 ;
74+ gen_data_1.Xl = 0.15 ;
75+ gen_data_1.S10 = 0 .;
76+ gen_data_1.S12 = 0 .;
77+
78+ BusFaultData<real_type, index_type> bus_fault_data_1;
79+ bus_fault_data_1.R = 0.0 ;
80+ bus_fault_data_1.X = 1e-3 ;
81+ bus_fault_data_1.status = false ;
82+
83+ //
84+ // Instantiate model components
85+ //
86+
87+ Bus<scalar_type, size_t > bus1 (bus_data_1);
88+ BusInfinite<scalar_type, size_t > bus2 (bus_data_2);
89+ Branch<scalar_type, size_t > branch (&bus1, &bus2, branch_data_1_2);
90+ BusFault<scalar_type, size_t > fault (&bus1, bus_fault_data_1);
91+ Genrou<scalar_type, size_t > gen (&bus1, gen_data_1);
92+
93+ //
94+ // Create the 2-bus system
95+ //
96+
97+ SystemModel<scalar_type, size_t > sys;
6198 sys.addBus (&bus1);
6299 sys.addBus (&bus2);
63100 sys.addComponent (&branch);
64101 sys.addComponent (&fault);
65102 sys.addComponent (&gen);
66103 sys.allocate ();
67104
68- double dt = 1.0 / 4.0 / 60.0 ;
105+ // Set time step to 1/4 of a 60Hz cycle
106+ real_type dt = 1.0 / 4.0 / 60.0 ;
69107
70108 // A data structure to keep track of the data we want to
71109 // compare to the reference solution. Rather than keeping
@@ -77,7 +115,9 @@ int main()
77115 // (plain ol' data), which have some benefits in C++.
78116 struct OutputData
79117 {
80- double ti, Vr, Vi, dw;
118+ // Output variables are time, real and imaginary voltage and
119+ // frequency deviation
120+ real_type ti, Vr, Vi, dw;
81121 };
82122
83123 // A list of output for each time step.
@@ -93,50 +133,50 @@ int main()
93133 // reference to that variable inside the callback). We select
94134 // the subset of the output we're interested in recording and
95135 // push it into output, which is updated outside the callback.
96- auto output_cb = [&](double t)
136+ auto output_cb = [&](real_type t)
97137 {
98- std::vector<double >& yval = sys.y ();
138+ std::vector<scalar_type >& y_val = sys.y ();
99139
100- output.push_back (OutputData{t, yval [0 ], yval [1 ], yval [3 ]});
140+ output.push_back (OutputData{t, y_val [0 ], y_val [1 ], y_val [3 ]});
101141 };
102142
103143 // Set up simulation
104- Ida<double , size_t > ida (&sys);
144+ Ida<scalar_type , size_t > ida (&sys);
105145 ida.configureSimulation ();
106146
107147 // Run simulation - making sure to pass the callback to record output
108- double start = static_cast <double >(clock ());
148+ real_type start = static_cast <real_type >(clock ());
109149
110150 // Run for 1s
111151 ida.initializeSimulation (0.0 , false );
112152 int nout = static_cast <int >(std::round ((1.0 - 0.0 ) / dt));
113153 ida.runSimulation (1.0 , nout, output_cb);
114154
115155 // Introduce fault and run for the next 0.1s
116- fault.setStatus (1 );
156+ fault.setStatus (true );
117157 ida.initializeSimulation (1.0 , false );
118158 nout = static_cast <int >(std::round ((1.1 - 1.0 ) / dt));
119159 ida.runSimulation (1.1 , nout, output_cb);
120160
121161 // Clear the fault and run until t = 10s.
122- fault.setStatus (0 );
162+ fault.setStatus (false );
123163 ida.initializeSimulation (1.1 , false );
124164 nout = static_cast <int >(std::round ((10.0 - 1.1 ) / dt));
125165 ida.runSimulation (10.0 , nout, output_cb);
126- double stop = static_cast <double >(clock ());
166+ real_type stop = static_cast <real_type >(clock ());
127167
128- double error_V = 0.0 ; // error in |V|
129- double error_w = 0.0 ; // error in rotor speed
168+ real_type error_V = 0.0 ; // error in |V|
169+ real_type error_w = 0.0 ; // error in rotor speed
130170
131171 // Read through the simulation data stored in the buffer.
132172 // Since we captured by reference, output should be available
133173 // for us to read here, outside the callback.
134174 for (size_t i = 0 ; i < output.size (); i++)
135175 {
136- OutputData data = output[i];
137- std::vector<double >& ref_sol = Example1::reference_solution[i + 1 ];
176+ OutputData data = output[i];
177+ std::vector<real_type >& ref_sol = Example1::reference_solution[i + 1 ];
138178
139- double err =
179+ real_type err =
140180 std::abs (std::sqrt (data.Vr * data.Vr + data.Vi * data.Vi ) - ref_sol[2 ])
141181 / (1.0 + std::abs (ref_sol[2 ]));
142182 if (err > error_V)
@@ -160,8 +200,9 @@ int main()
160200 // std::cout << "\n";
161201 }
162202
163- double error_V_allowed = 2e-4 ;
164- double error_w_allowed = 1e-4 ;
203+ // Errors allowed for agreement with Powerworld results
204+ real_type error_V_allowed = 2e-4 ;
205+ real_type error_w_allowed = 1e-4 ;
165206
166207 // Tolerances based on Powerworld reference accuracy
167208 int status = 0 ;
0 commit comments