Skip to content

Commit b54744f

Browse files
nkoukpaizanalexander-novo
authored and
Wiktoria Zielinska
committed
Verify Branch and Load Jacobians computed with Sparse::Variable (#89)
--------- Co-authored-by: Alexander Novotny <[email protected]>
1 parent 50233c4 commit b54744f

File tree

12 files changed

+213
-0
lines changed

12 files changed

+213
-0
lines changed

‎src/LinearAlgebra/SparsityPattern/Variable.hpp‎

Lines changed: 15 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -14,6 +14,8 @@
1414
#include <string>
1515
#include <vector>
1616

17+
#include <ScalarTraits.hpp>
18+
1719
namespace GridKit
1820
{
1921
namespace Sparse
@@ -326,5 +328,18 @@ namespace GridKit
326328
} // namespace Sparse
327329
} // namespace GridKit
328330

331+
namespace GridKit
332+
{
333+
template <>
334+
class ScalarTraits<Sparse::Variable>
335+
{
336+
public:
337+
typedef double real_type;
338+
typedef double norm_type;
339+
typedef double scalar_type;
340+
};
341+
342+
} // namespace GridKit
343+
329344
#include "VariableImplementation.hpp"
330345
#include "VariableOperators.hpp"

‎src/Model/PhasorDynamics/Branch/Branch.cpp‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -221,6 +221,8 @@ namespace GridKit
221221
// Available template instantiations
222222
template class Branch<double, long int>;
223223
template class Branch<double, size_t>;
224+
template class Branch<Sparse::Variable, long int>;
225+
template class Branch<Sparse::Variable, size_t>;
224226

225227
} // namespace PhasorDynamics
226228
} // namespace GridKit

‎src/Model/PhasorDynamics/Bus/Bus.cpp‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -215,6 +215,8 @@ namespace GridKit
215215
// Available template instantiations
216216
template class Bus<double, long int>;
217217
template class Bus<double, size_t>;
218+
template class Bus<Sparse::Variable, long int>;
219+
template class Bus<Sparse::Variable, size_t>;
218220

219221
} // namespace PhasorDynamics
220222
} // namespace GridKit

‎src/Model/PhasorDynamics/Bus/BusInfinite.cpp‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -195,6 +195,8 @@ namespace GridKit
195195
// Available template instantiations
196196
template class BusInfinite<double, long int>;
197197
template class BusInfinite<double, size_t>;
198+
template class BusInfinite<Sparse::Variable, long int>;
199+
template class BusInfinite<Sparse::Variable, size_t>;
198200

199201
} // namespace PhasorDynamics
200202
} // namespace GridKit

‎src/Model/PhasorDynamics/BusBase.hpp‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -2,6 +2,7 @@
22

33
#include <vector>
44

5+
#include <LinearAlgebra/SparsityPattern/Variable.hpp>
56
#include <Model/Evaluator.hpp>
67

78
namespace GridKit

‎src/Model/PhasorDynamics/Component.hpp‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -2,6 +2,7 @@
22

33
#include <vector>
44

5+
#include <LinearAlgebra/SparsityPattern/Variable.hpp>
56
#include <Model/Evaluator.hpp>
67

78
namespace GridKit

‎src/Model/PhasorDynamics/Load/Load.cpp‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -173,6 +173,8 @@ namespace GridKit
173173
// Available template instantiations
174174
template class Load<double, long int>;
175175
template class Load<double, size_t>;
176+
template class Load<Sparse::Variable, long int>;
177+
template class Load<Sparse::Variable, size_t>;
176178

177179
} // namespace PhasorDynamics
178180
} // namespace GridKit

‎src/Utilities/Testing.hpp‎

Lines changed: 53 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -10,6 +10,8 @@
1010
#include <cmath>
1111
#include <iostream>
1212
#include <limits>
13+
#include <map>
14+
#include <vector>
1315

1416
#include <PowerSystemData.hpp>
1517

@@ -220,6 +222,57 @@ namespace GridKit
220222
return fail == 0;
221223
}
222224

225+
/**
226+
* @brief Equatlity comparison between maps with a tolerance for the scalar value
227+
*
228+
* @tparam IdxT
229+
* @tparam RealT
230+
*
231+
* @param[in] a - first map to compare
232+
* @param[in] b - second map to compare
233+
* @param[in] tol - tolerance
234+
* @return bool - true if the maps are equal; false otherwise
235+
*/
236+
template <typename IdxT = size_t, typename RealT = double>
237+
inline bool isEqual(std::map<IdxT, RealT> a,
238+
std::map<IdxT, RealT> b,
239+
const RealT tol = std::numeric_limits<RealT>::epsilon())
240+
{
241+
int fail = 0;
242+
243+
if (a.size() != b.size())
244+
{
245+
fail++;
246+
errs() << "Containers do not have the same size! "
247+
<< "a.size() = " << a.size() << ", and "
248+
<< "b.size() = " << b.size() << "\n";
249+
}
250+
251+
for (const auto& pair_a : a)
252+
{
253+
auto it = b.find(pair_a.first);
254+
if (it != b.end())
255+
{
256+
if (!isEqual(pair_a.second, it->second, tol))
257+
{
258+
fail++;
259+
errs() << "Mismatching map values! "
260+
<< "a.first = " << pair_a.first << ", "
261+
<< "a.second = " << pair_a.second << ", and "
262+
<< "b.second = " << it->second << "\n";
263+
}
264+
}
265+
else
266+
{
267+
fail++;
268+
errs() << "Entry not found in the second container! "
269+
<< "a.first = " << pair_a.first << "\n";
270+
}
271+
}
272+
273+
return fail == 0;
274+
}
275+
223276
} // namespace Testing
224277

225278
} // namespace GridKit

‎tests/UnitTests/PhasorDynamics/BranchTests.hpp‎

Lines changed: 79 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,7 @@
11
#include <iomanip>
22
#include <iostream>
33

4+
#include <LinearAlgebra/SparsityPattern/Variable.hpp>
45
#include <Model/PhasorDynamics/Branch/Branch.hpp>
56
#include <Model/PhasorDynamics/Bus/Bus.hpp>
67
#include <Model/PhasorDynamics/Bus/BusInfinite.hpp>
@@ -77,6 +78,46 @@ namespace GridKit
7778
return success.report(__func__);
7879
}
7980

81+
TestOutcome jacobian()
82+
{
83+
TestStatus success = true;
84+
85+
real_type R{2.0}; ///< Branch series resistance
86+
real_type X{4.0}; ///< Branch series reactance
87+
real_type G{0.2}; ///< Branch shunt conductance
88+
real_type B{1.2}; ///< Branch shunt charging
89+
90+
Sparse::Variable Vr1{10.0}; ///< Bus-1 real voltage
91+
Sparse::Variable Vi1{20.0}; ///< Bus-1 imaginary voltage
92+
Sparse::Variable Vr2{30.0}; ///< Bus-2 real voltage
93+
Sparse::Variable Vi2{40.0}; ///< Bus-2 imaginary voltage
94+
95+
Vr1.setVariableNumber(0); ///< Independent variables: first
96+
Vi1.setVariableNumber(1); ///< Independent variables: second
97+
Vr2.setVariableNumber(2); ///< Independent variables: third
98+
Vi2.setVariableNumber(3); ///< Independent variables: fourth
99+
100+
PhasorDynamics::BusInfinite<Sparse::Variable, IdxT> bus1(Vr1, Vi1);
101+
PhasorDynamics::BusInfinite<Sparse::Variable, IdxT> bus2(Vr2, Vi2);
102+
103+
PhasorDynamics::Branch<Sparse::Variable, IdxT> branch(&bus1, &bus2, R, X, G, B);
104+
branch.evaluateResidual(); ///< Computes the residual and the Jacobian values by tracking
105+
///< the dependencies
106+
107+
std::vector<Sparse::Variable> residuals{bus1.Ir(), bus1.Ii(), bus2.Ir(), bus2.Ii()};
108+
std::vector<Sparse::Variable::DependencyMap> ref = analyticalJacobian(R, X, G, B);
109+
110+
/// Compare dependencies computed automatically to the ones computed analytically
111+
for (size_t i = 0; i < residuals.size(); ++i)
112+
{
113+
Sparse::Variable res = residuals[i];
114+
const Sparse::Variable::DependencyMap& dependencies = res.getDependencies();
115+
success *= (GridKit::Testing::isEqual(dependencies, ref[i]));
116+
}
117+
118+
return success.report(__func__);
119+
}
120+
80121
TestOutcome accessors()
81122
{
82123
TestStatus success = true;
@@ -139,6 +180,44 @@ namespace GridKit
139180

140181
return success.report(__func__);
141182
}
183+
184+
private:
185+
std::vector<Sparse::Variable::DependencyMap> analyticalJacobian(const real_type R,
186+
const real_type X,
187+
const real_type G,
188+
const real_type B)
189+
{
190+
const real_type b = -X / (R * R + X * X);
191+
const real_type g = R / (R * R + X * X);
192+
193+
real_type dIr1_dVr1 = -(g + 0.5 * G);
194+
real_type dIr1_dVi1 = (b + 0.5 * B);
195+
real_type dIr1_dVr2 = g;
196+
real_type dIr1_dVi2 = -b;
197+
198+
real_type dIi1_dVr1 = -(b + 0.5 * B);
199+
real_type dIi1_dVi1 = -(g + 0.5 * G);
200+
real_type dIi1_dVr2 = b;
201+
real_type dIi1_dVi2 = g;
202+
203+
real_type dIr2_dVr1 = g;
204+
real_type dIr2_dVi1 = -b;
205+
real_type dIr2_dVr2 = -(g + 0.5 * G);
206+
real_type dIr2_dVi2 = (b + 0.5 * B);
207+
208+
real_type dIi2_dVr1 = b;
209+
real_type dIi2_dVi1 = g;
210+
real_type dIi2_dVr2 = -(b + 0.5 * B);
211+
real_type dIi2_dVi2 = -(g + 0.5 * G);
212+
213+
std::vector<Sparse::Variable::DependencyMap> dependencies(4);
214+
dependencies[0] = {{0, dIr1_dVr1}, {1, dIr1_dVi1}, {2, dIr1_dVr2}, {3, dIr1_dVi2}};
215+
dependencies[1] = {{0, dIi1_dVr1}, {1, dIi1_dVi1}, {2, dIi1_dVr2}, {3, dIi1_dVi2}};
216+
dependencies[2] = {{0, dIr2_dVr1}, {1, dIr2_dVi1}, {2, dIr2_dVr2}, {3, dIr2_dVi2}};
217+
dependencies[3] = {{0, dIi2_dVr1}, {1, dIi2_dVi1}, {2, dIi2_dVr2}, {3, dIi2_dVi2}};
218+
219+
return dependencies;
220+
}
142221
}; // class BranchTest
143222

144223
} // namespace Testing

‎tests/UnitTests/PhasorDynamics/LoadTests.hpp‎

Lines changed: 54 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -3,6 +3,7 @@
33
#include <iomanip>
44
#include <iostream>
55

6+
#include <LinearAlgebra/SparsityPattern/Variable.hpp>
67
#include <Model/PhasorDynamics/Bus/Bus.hpp>
78
#include <Model/PhasorDynamics/Bus/BusInfinite.hpp>
89
#include <Model/PhasorDynamics/Load/Load.hpp>
@@ -65,6 +66,59 @@ namespace GridKit
6566

6667
return success.report(__func__);
6768
}
69+
70+
TestOutcome jacobian()
71+
{
72+
TestStatus success = true;
73+
74+
real_type R{2.0}; ///< Load resistance
75+
real_type X{4.0}; ///< Load reactance
76+
77+
Sparse::Variable Vr{10.0}; ///< Bus real voltage
78+
Sparse::Variable Vi{20.0}; ///< Bus imaginary voltage
79+
80+
Vr.setVariableNumber(0); ///< Independent variables: first
81+
Vi.setVariableNumber(1); ///< Independent variables: second
82+
83+
PhasorDynamics::BusInfinite<Sparse::Variable, IdxT> bus(Vr, Vi);
84+
85+
PhasorDynamics::Load<Sparse::Variable, IdxT> load(&bus, R, X);
86+
load.evaluateResidual(); ///< Computes the residual and the Jacobian values by tracking
87+
///< the dependencies
88+
89+
std::vector<Sparse::Variable> residuals{bus.Ir(), bus.Ii()};
90+
std::vector<Sparse::Variable::DependencyMap> ref = analyticalJacobian(R, X);
91+
92+
/// Compare dependencies computed automatically to the ones computed analytically
93+
for (size_t i = 0; i < residuals.size(); ++i)
94+
{
95+
Sparse::Variable res = residuals[i];
96+
const Sparse::Variable::DependencyMap& dependencies = res.getDependencies();
97+
success *= (GridKit::Testing::isEqual(dependencies, ref[i]));
98+
}
99+
100+
return success.report(__func__);
101+
}
102+
103+
private:
104+
std::vector<Sparse::Variable::DependencyMap> analyticalJacobian(const real_type R,
105+
const real_type X)
106+
{
107+
const real_type b = -X / (R * R + X * X);
108+
const real_type g = R / (R * R + X * X);
109+
110+
real_type dIr_dVr = -g;
111+
real_type dIr_dVi = -b;
112+
113+
real_type dIi_dVr = b;
114+
real_type dIi_dVi = -g;
115+
116+
std::vector<Sparse::Variable::DependencyMap> dependencies(2);
117+
dependencies[0] = {{0, dIr_dVr}, {1, dIr_dVi}};
118+
dependencies[1] = {{0, dIi_dVr}, {1, dIi_dVi}};
119+
120+
return dependencies;
121+
}
68122
};
69123

70124
} // namespace Testing

0 commit comments

Comments
 (0)