Skip to content

Commit

Permalink
rework review comments pt3
Browse files Browse the repository at this point in the history
  • Loading branch information
@leovsch
leovsch committed Feb 12, 2025
1 parent 0ccb30a commit 845b780
Show file tree
Hide file tree
Showing 4 changed files with 29 additions and 52 deletions.
2 changes: 1 addition & 1 deletion power_grid_model_c/power_grid_model/include/power_grid_model/common/matrix_utils.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,7 +11,7 @@ inline DoubleComplex average_of_diagonal_of_matrix(const ComplexTensor<asymmetri
}

inline DoubleComplex average_of_off_diagonal_of_matrix(const ComplexTensor<asymmetric_t> &matrix) {
return (matrix(0,2) + matrix(1,1) + matrix(2,0)) / 3.0;
return (matrix(0,1) + matrix(1,2) + matrix(1,0) + matrix(1,2) + matrix(2,0) + matrix(2,1)) / 6.0;
}

} // namespace power_grid_model
35 changes: 17 additions & 18 deletions power_grid_model_c/power_grid_model/include/power_grid_model/component/asym_line.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -31,8 +31,8 @@ class AsymLine : public Branch {
ComplexTensor<asymmetric_t> c_matrix = compute_c_matrix_from_input(asym_line_input);
ComplexTensor<asymmetric_t> z_series = compute_z_series_from_input(asym_line_input);

y_series = inv(z_series);
y_shunt = 2 * pi * system_frequency * c_matrix * 1.0i;
_y_series_abc = inv(z_series);
_y_shunt_abc = 2 * pi * system_frequency * c_matrix * 1.0i;
}

// override getter
Expand All @@ -45,8 +45,8 @@ class AsymLine : public Branch {
private:
double i_n_;
double base_i_;
ComplexTensor<asymmetric_t> y_series;
ComplexTensor<asymmetric_t> y_shunt;
ComplexTensor<asymmetric_t> _y_series_abc;
ComplexTensor<asymmetric_t> _y_shunt_abc;

ComplexTensor<asymmetric_t> compute_z_series_from_input(const power_grid_model::AsymLineInput& asym_line_input) {
ComplexTensor<asymmetric_t> z_series_abc;
Expand All @@ -58,34 +58,30 @@ class AsymLine : public Branch {
else {
ComplexTensor4 r_matrix = ComplexTensor4(asym_line_input.r_aa, asym_line_input.r_bb, asym_line_input.r_cc, asym_line_input.r_nn, asym_line_input.r_ba, asym_line_input.r_ca, asym_line_input.r_na, asym_line_input.r_cb, asym_line_input.r_nb, asym_line_input.r_nc);
ComplexTensor4 x_matrix = ComplexTensor4(asym_line_input.x_aa, asym_line_input.x_bb, asym_line_input.x_cc, asym_line_input.x_nn, asym_line_input.x_ba, asym_line_input.x_ca, asym_line_input.x_na, asym_line_input.x_cb, asym_line_input.x_nb, asym_line_input.x_nc);

ComplexTensor4 y = r_matrix + 1.0i * x_matrix;
z_series_abc = kron_reduction(y);
}
ComplexTensor<asymmetric_t> a_matrix = get_sym_matrix();
ComplexTensor<asymmetric_t> a_matrix_inv = get_sym_matrix_inv();
ComplexTensor<asymmetric_t> z_series = (a_matrix_inv.matrix() * z_series_abc.matrix() * a_matrix.matrix()).array();
return z_series;
return z_series_abc;
}

ComplexTensor<asymmetric_t> compute_c_matrix_from_input(const power_grid_model::AsymLineInput& asym_line_input) {
ComplexTensor<asymmetric_t> c_matrix;
if (!is_nan(asym_line_input.c0) && !is_nan(asym_line_input.c1)) {
c_matrix = ComplexTensor<asymmetric_t>(asym_line_input.c0 + asym_line_input.c1, -asym_line_input.c1);
c_matrix = ComplexTensor<asymmetric_t>{(2.0 * asym_line_input.c1 + asym_line_input.c0) / 3.0, (asym_line_input.c0 - asym_line_input.c1) / 3.0 };
}
else if (is_nan(asym_line_input.c_nn)) {
c_matrix = ComplexTensor<asymmetric_t>(asym_line_input.c_aa, asym_line_input.c_ba, asym_line_input.c_bb, asym_line_input.c_ca, asym_line_input.c_cb, asym_line_input.c_cc);
c_matrix = ComplexTensor<asymmetric_t>(asym_line_input.c_aa, asym_line_input.c_bb, asym_line_input.c_cc, asym_line_input.c_ba, asym_line_input.c_ca, asym_line_input.c_cb);
}
else {
ComplexTensor4 c_matrix_neutral = ComplexTensor4(asym_line_input.c_aa, asym_line_input.c_ba, asym_line_input.c_bb, asym_line_input.c_ca, asym_line_input.c_cb, asym_line_input.c_cc, asym_line_input.c_na, asym_line_input.c_nb, asym_line_input.c_nc, asym_line_input.c_nn);
ComplexTensor4 c_matrix_neutral = ComplexTensor4(asym_line_input.c_aa, asym_line_input.c_bb, asym_line_input.c_cc, asym_line_input.c_nn, asym_line_input.c_ba, asym_line_input.c_ca, asym_line_input.c_na, asym_line_input.c_cb, asym_line_input.c_nb, asym_line_input.c_nc);
c_matrix = kron_reduction(c_matrix_neutral);
}
return c_matrix;
}

BranchCalcParam<symmetric_t> sym_calc_param() const override {
DoubleComplex y1_series_ = average_of_diagonal_of_matrix(y_series) - average_of_off_diagonal_of_matrix(y_series);
DoubleComplex y1_shunt_ = average_of_diagonal_of_matrix(y_shunt) - average_of_off_diagonal_of_matrix(y_shunt);
DoubleComplex y1_series_ = average_of_diagonal_of_matrix(_y_series_abc) - average_of_off_diagonal_of_matrix(_y_series_abc);
DoubleComplex y1_shunt_ = average_of_diagonal_of_matrix(_y_shunt_abc) - average_of_off_diagonal_of_matrix(_y_shunt_abc);
return calc_param_y_sym(y1_series_, y1_shunt_, 1.0);
}

Expand All @@ -96,18 +92,21 @@ class AsymLine : public Branch {
// single connected
if (from_status() || to_status()) {
// branch_shunt = 0.5 * y_shunt + 1.0 / (1.0 / y_series + 2.0 / y_shunt);
ComplexTensor<asymmetric_t> branch_shunt = 0.5 * inv(y_shunt) + inv(inv(y_series) + 2.0 * inv(y_shunt));
ComplexTensor<asymmetric_t> branch_shunt = ComplexTensor<asymmetric_t>();
if ((cabs(_y_shunt_abc) >= numerical_tolerance).all()) {
branch_shunt = 0.5 * inv(_y_shunt_abc) + inv(inv(_y_series_abc) + 2.0 * inv(_y_shunt_abc));
}
// from or to connected
param.yff() = from_status() ? branch_shunt : ComplexTensor<asymmetric_t>();
param.ytt() = to_status() ? branch_shunt : ComplexTensor<asymmetric_t>();
}
}
// both connected
else {
param.ytt() = y_series + 0.5 * y_shunt;
param.ytt() = _y_series_abc + 0.5 * _y_shunt_abc;
param.yff() = param.ytt();
param.yft() = y_series;
param.ytf() = y_series;
param.yft() = _y_series_abc;
param.ytf() = _y_series_abc;
}
return param;
}
Expand Down
2 changes: 1 addition & 1 deletion power_grid_model_c/power_grid_model/include/power_grid_model/component/line_utils.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,7 @@ namespace power_grid_model {

inline ComplexTensor<asymmetric_t> kron_reduction(const ComplexTensor4& matrix_to_reduce) {
ComplexTensor4 Y = matrix_to_reduce;
ComplexTensor<asymmetric_t> Y_aa = ComplexTensor<asymmetric_t>(Y(0,0), Y(1,0), Y(1, 1), Y(2,0), Y(2,1), Y(2,2));
ComplexTensor<asymmetric_t> Y_aa = ComplexTensor<asymmetric_t>(Y(0,0), Y(1,1), Y(2, 2), Y(1,0), Y(2,0), Y(2,1));
ComplexValue<asymmetric_t> Y_ab(Y(0,3), Y(1,3), Y(2,3));
ComplexValue<asymmetric_t> Y_ba(Y(3,0), Y(3,1), Y(3,2));
DoubleComplex Y_bb_inv = 1.0 / Y(3,3);
Expand Down
42 changes: 10 additions & 32 deletions tests/cpp_unit_tests/test_asym_line.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -44,12 +44,11 @@ TEST_CASE("Test asym line") {
double const base_i = base_power_1p / (10.0e3 / sqrt3);
double const base_y = base_i * base_i / base_power_1p;
Branch& branch = asym_line;
ComplexTensor<asymmetric_t> const y_series = ComplexTensor<asymmetric_t>(0.66303418-0.34266364i, -0.02971114+0.03783535i, 0.04762194+0.00681293i, 0.04762194+0.00681293i, 2.48612768-0.46271628i, 0.05942228-0.07567069i, -0.02971114+0.03783535i, -0.09524388-0.01362585i, 2.48612768-0.46271628i);
ComplexTensor<asymmetric_t> const y_series = ComplexTensor<asymmetric_t>(1.87842984-0.42269873i, 1.87842984-0.42269873i, 1.87842984-0.42269873i, -0.62560863-0.00463073i, -0.57187623+0.12931409i, -0.62560863-0.00463073i);
ComplexTensor<asymmetric_t> const y_shunt = 2 * pi * system_frequency * ComplexTensor<asymmetric_t>{(2.0 * input.c1 + input.c0) / 3.0, (input.c0 - input.c1) / 3.0 } * 1.0i;

ComplexTensor<asymmetric_t> const y_shunt = 2 * pi * system_frequency * ComplexTensor<asymmetric_t>(input.c0 + input.c1, -input.c1) * 1.0i;

DoubleComplex const y1_series = (y_series(0,0) + y_series(1,1) + y_series(2,2)) / 3.0 - (y_series(0,2) + y_series(1,1) + y_series(2,0)) / 3.0;
DoubleComplex const y1_shunt = (y_shunt(0,0) + y_shunt(1,1) + y_shunt(2,2)) / 3.0 - (y_shunt(0,2) + y_shunt(1,1) + y_shunt(2,0)) / 3.0;
DoubleComplex const y1_series = (y_series(0,0) + y_series(1,1) + y_series(2,2)) / 3.0 - (y_series(0,1) + y_series(1,2) + y_series(1,0) + y_series(1,2) + y_series(2,0) + y_series(2,1)) / 6.0;
DoubleComplex const y1_shunt = (y_shunt(0,0) + y_shunt(1,1) + y_shunt(2,2)) / 3.0 - (y_shunt(0,1) + y_shunt(1,2) + y_shunt(1,0) + y_shunt(1,2) + y_shunt(2,0) + y_shunt(2,1)) / 6.0;

// symmetric
DoubleComplex const yff1 = y1_series + 0.5 * y1_shunt;
Expand All @@ -58,28 +57,7 @@ TEST_CASE("Test asym line") {

// asymmetric
ComplexTensor<asymmetric_t> ytt = y_series + 0.5 * y_shunt;
DoubleComplex const yff0 = ytt(0,0);
DoubleComplex const yft0 = -y_series(0,0);
DoubleComplex const ys0 = 0.5 * y_shunt(0,0) + 1.0 / (1.0 / y_series(0,0) + 2.0 / y_shunt(0,0));
ComplexTensor<asymmetric_t> const yffa{(2.0 * yff1 + yff0) / 3.0, (yff0 - yff1) / 3.0};
ComplexTensor<asymmetric_t> const yfta{(2.0 * yft1 + yft0) / 3.0, (yft0 - yft1) / 3.0};
ComplexTensor<asymmetric_t> const ysa{(2.0 * ys1 + ys0) / 3.0, (ys0 - ys1) / 3.0};

DoubleComplex const u1f = 1.0;
DoubleComplex const u1t = 0.9;
ComplexValue<asymmetric_t> const uaf{1.0};
ComplexValue<asymmetric_t> const uat{0.9};
DoubleComplex const i1f = (yff1 * u1f + yft1 * u1t) * base_i;
DoubleComplex const i1t = (yft1 * u1f + yff1 * u1t) * base_i;
DoubleComplex const s_f = conj(i1f) * u1f * 10e3 * sqrt3;
DoubleComplex const s_t = conj(i1t) * u1t * 10e3 * sqrt3;
double const loading = std::max(cabs(i1f), cabs(i1t)) / 200.0;

// Short circuit results
DoubleComplex const if_sc{1.0, 1.0};
DoubleComplex const it_sc{2.0, 2.0 * sqrt(3)};
ComplexValue<asymmetric_t> const if_sc_asym{1.0 + 1.0i};
ComplexValue<asymmetric_t> const it_sc_asym{2.0 + (2.0i * sqrt(3))};
ComplexTensor<asymmetric_t> branch_shunt = 0.5 * inv(y_shunt) + inv(inv(y_series) + 2.0 * inv(y_shunt));

CHECK(asym_line.math_model_type() == ComponentType::branch);

Expand Down Expand Up @@ -134,10 +112,10 @@ TEST_CASE("Test asym line") {
SUBCASE("Asymmetric parameters") {
// double connected
BranchCalcParam<asymmetric_t> param = asym_line.calc_param<asymmetric_t>();
CHECK((cabs(param.yff() - yffa) < numerical_tolerance).all());
CHECK((cabs(param.ytt() - yffa) < numerical_tolerance).all());
CHECK((cabs(param.ytf() - yfta) < numerical_tolerance).all());
CHECK((cabs(param.yft() - yfta) < numerical_tolerance).all());
CHECK((cabs(param.yff() - ytt) < numerical_tolerance).all());
CHECK((cabs(param.ytt() - ytt) < numerical_tolerance).all());
CHECK((cabs(param.ytf() - y_series) < numerical_tolerance).all());
CHECK((cabs(param.yft() - y_series) < numerical_tolerance).all());
// no source
param = branch.calc_param<asymmetric_t>(false);
CHECK((cabs(param.yff() - 0.0) < numerical_tolerance).all());
Expand All @@ -147,7 +125,7 @@ TEST_CASE("Test asym line") {
// from connected
CHECK(branch.set_status(na_IntS, false));
param = asym_line.calc_param<asymmetric_t>();
CHECK((cabs(param.yff() - ysa) < numerical_tolerance).all());
CHECK((cabs(param.yff() - branch_shunt) < numerical_tolerance).all()); // Fail
CHECK((cabs(param.ytt() - 0.0) < numerical_tolerance).all());
CHECK((cabs(param.ytf() - 0.0) < numerical_tolerance).all());
CHECK((cabs(param.yft() - 0.0) < numerical_tolerance).all());
Expand Down

0 comments on commit 845b780

Please sign in to comment.