Fix zero-sequence model of Yyn transformer

code_generation/data/attribute_classes/input.json

@@ -223,7 +223,7 @@
                 },
                 {
                     "data_type": "double",
-                    "names": ["uk", "pk", "i0", "p0"],
+                    "names": ["uk", "pk", "i0", "p0", "i0_zero_sequence", "p0_zero_sequence"],
                     "description": "short circuit and open testing parameters"
                 },
                 {

power_grid_model_c/power_grid_model/include/power_grid_model/auxiliary/input.hpp

@@ -195,6 +195,8 @@ struct TransformerInput {
     double pk{nan};  // short circuit and open testing parameters
     double i0{nan};  // short circuit and open testing parameters
     double p0{nan};  // short circuit and open testing parameters
+    double i0_zero_sequence{nan};  // short circuit and open testing parameters
+    double p0_zero_sequence{nan};  // short circuit and open testing parameters
     WindingType winding_from{static_cast<WindingType>(na_IntS)};  // winding type at each side
     WindingType winding_to{static_cast<WindingType>(na_IntS)};  // winding type at each side
     IntS clock{na_IntS};  // clock number

power_grid_model_c/power_grid_model/include/power_grid_model/auxiliary/meta_gen/input.hpp

@@ -188,7 +188,7 @@ struct get_attributes_list<LinkInput> {
 
 template<>
 struct get_attributes_list<TransformerInput> {
-    static constexpr std::array<MetaAttribute, 29> value{
+    static constexpr std::array<MetaAttribute, 31> value{
             // all attributes including base class
 
             meta_data_gen::get_meta_attribute<&TransformerInput::id>(offsetof(TransformerInput, id), "id"),
@@ -203,6 +203,8 @@ struct get_attributes_list<TransformerInput> {
             meta_data_gen::get_meta_attribute<&TransformerInput::pk>(offsetof(TransformerInput, pk), "pk"),
             meta_data_gen::get_meta_attribute<&TransformerInput::i0>(offsetof(TransformerInput, i0), "i0"),
             meta_data_gen::get_meta_attribute<&TransformerInput::p0>(offsetof(TransformerInput, p0), "p0"),
+            meta_data_gen::get_meta_attribute<&TransformerInput::i0_zero_sequence>(offsetof(TransformerInput, i0_zero_sequence), "i0_zero_sequence"),
+            meta_data_gen::get_meta_attribute<&TransformerInput::p0_zero_sequence>(offsetof(TransformerInput, p0_zero_sequence), "p0_zero_sequence"),
             meta_data_gen::get_meta_attribute<&TransformerInput::winding_from>(offsetof(TransformerInput, winding_from), "winding_from"),
             meta_data_gen::get_meta_attribute<&TransformerInput::winding_to>(offsetof(TransformerInput, winding_to), "winding_to"),
             meta_data_gen::get_meta_attribute<&TransformerInput::clock>(offsetof(TransformerInput, clock), "clock"),

power_grid_model_c/power_grid_model/include/power_grid_model/component/transformer.hpp

@@ -18,6 +18,14 @@
 namespace power_grid_model {
 
 class Transformer : public Branch {
+  private:
+    struct TransformerParams {
+        DoubleComplex y_series{};
+        DoubleComplex y_shunt{};
+        DoubleComplex y0_shunt{};
+        double k{1.0};
+    };
+
   public:
     using InputType = TransformerInput;
     using UpdateType = TransformerUpdate;
@@ -33,6 +41,10 @@
           pk_{transformer_input.pk},
           i0_{transformer_input.i0},
           p0_{transformer_input.p0},
+          i0_zero_sequence_{is_nan(transformer_input.i0_zero_sequence) ? i0_ : transformer_input.i0_zero_sequence},
+          p0_zero_sequence_{is_nan(transformer_input.p0_zero_sequence)
+                                ? p0_ + pk_ * (i0_zero_sequence_ * i0_zero_sequence_ - i0_ * i0_)
+                                : transformer_input.p0_zero_sequence},
           winding_from_{transformer_input.winding_from},
           winding_to_{transformer_input.winding_to},
           clock_{transformer_input.clock},
@@ -120,6 +132,8 @@
     double pk_;
     double i0_;
     double p0_;
+    double i0_zero_sequence_;
+    double p0_zero_sequence_;
     WindingType winding_from_;
     WindingType winding_to_;
     IntS clock_;
@@ -156,7 +170,7 @@
     }
 
     // calculate transformer parameter
-    std::tuple<DoubleComplex, DoubleComplex, double> transformer_params() const {
+    TransformerParams transformer_params() const {
         double const base_y_to = base_i_to_ * base_i_to_ / base_power_1p;
         // off nominal tap ratio
         auto const [u1, u2] = [this]() {
@@ -192,29 +206,39 @@
         z_series.imag(uk_sign * (z_series_imag_squared > 0.0 ? std::sqrt(z_series_imag_squared) : 0.0));
         // y series
         y_series = (1.0 / z_series) / base_y_to;
+
         // shunt
         DoubleComplex y_shunt;
-        // Y = I0_2 / (U2/sqrt3) = i0 * (S / sqrt3 / U2) / (U2/sqrt3) = i0 * S * / U2 / U2
+        // Y = I0 / (U2/sqrt3) = i0 * (S / sqrt3 / U2) / (U2/sqrt3) = i0 * S * / U2 / U2
         double const y_shunt_abs = i0_ * sn_ / u2 / u2;
         // G = P0 / (U2^2)
         y_shunt.real(p0_ / u2 / u2);
-
         auto const y_shunt_imag_squared = y_shunt_abs * y_shunt_abs - y_shunt.real() * y_shunt.real();
         y_shunt.imag(y_shunt_imag_squared > 0.0 ? -std::sqrt(y_shunt_imag_squared) : 0.0);
-
-        // y shunt
         y_shunt = y_shunt / base_y_to;
+
+        // shunt zero sequence
+        DoubleComplex y0_shunt;
+        // Y0 = I0_0 / (U2/sqrt3) = i0_zero_sequence_ * (S / sqrt3 / U2) / (U2/sqrt3) = i0_zero_sequence_ * S / U2 /
+        // U2
+        double const y0_shunt_abs = i0_zero_sequence_ * sn_ / u2 / u2;
+        // G0 = P0_0 / (U2^2)
+        y0_shunt.real(p0_zero_sequence_ / u2 / u2);
+        auto const y0_shunt_imag_squared = y0_shunt_abs * y0_shunt_abs - y0_shunt.real() * y0_shunt.real();
+        y0_shunt.imag(y0_shunt_imag_squared > 0.0 ? -std::sqrt(y0_shunt_imag_squared) : 0.0);
+        y0_shunt = y0_shunt / base_y_to;
+
         // return
-        return std::make_tuple(y_series, y_shunt, k);
+        return TransformerParams{.y_series = y_series, .y_shunt = y_shunt, .y0_shunt = y0_shunt, .k = k};
     }
 
     // branch param
     BranchCalcParam<symmetric_t> sym_calc_param() const final {
-        auto const [y_series, y_shunt, k] = transformer_params();
+        auto const [y_series, y_shunt, y0_shunt, k] = transformer_params();
         return calc_param_y_sym(y_series, y_shunt, k * std::exp(1.0i * (clock_ * deg_30)));
     }
     BranchCalcParam<asymmetric_t> asym_calc_param() const final {
-        auto const [y_series, y_shunt, k] = transformer_params();
+        auto const [y_series, y_shunt, y0_shunt, k] = transformer_params();
         // positive sequence
         auto const param1 = calc_param_y_sym(y_series, y_shunt, k * std::exp(1.0i * (clock_ * deg_30)));
         // negative sequence
@@ -230,19 +254,37 @@
             }
             DoubleComplex const z0_series = 1.0 / y_series + 3.0 * (z_grounding_to_ + z_grounding_from_ / k / k);
             DoubleComplex const y0_series = 1.0 / z0_series;
-            param0 = calc_param_y_sym(y0_series, y_shunt, k * std::exp(1.0i * phase_shift_0));
+            param0 = calc_param_y_sym(y0_series, y0_shunt, k * std::exp(1.0i * phase_shift_0));
         }
-        // YNd
-        if (winding_from_ == WindingType::wye_n && winding_to_ == WindingType::delta && from_status()) {
-            DoubleComplex const z0_series = 1.0 / y_series + 3.0 * z_grounding_from_ / k / k;
-            DoubleComplex const y0_series = 1.0 / z0_series;
-            param0.yff() = (y0_series + y_shunt) / k / k;
+        // YN*
+        else if (winding_from_ == WindingType::wye_n && from_status()) {
+            // ground path always possible via magnetization branch
+            DoubleComplex y0 = y0_shunt;
+            if (winding_to_ == WindingType::delta) {
+                // additional path via zk
+                y0 += y_series;
+            }
+            if (y0 != DoubleComplex{0.0, 0.0}) {
+                // avoid division by zero
+                DoubleComplex const z0 = 1.0 / y0 + 3.0 * z_grounding_from_ / k / k;
+                y0 = 1.0 / z0;
+                param0.yff() = y0 / k / k;
+            }
         }
-        // Dyn
-        if (winding_from_ == WindingType::delta && winding_to_ == WindingType::wye_n && to_status()) {
-            DoubleComplex const z0_series = 1.0 / y_series + 3.0 * z_grounding_to_;
-            DoubleComplex const y0_series = 1.0 / z0_series;
-            param0.ytt() = (y0_series + y_shunt);
+        // *yn
+        else if (winding_to_ == WindingType::wye_n && to_status()) {
+            // ground path always possible via magnetization branch
+            DoubleComplex y0 = y0_shunt;
+            if (winding_from_ == WindingType::delta) {
+                // additional path via zk
+                y0 += y_series;
+            }
+            if (y0 != DoubleComplex{0.0, 0.0}) {
+                // avoid division by zero
+                DoubleComplex const z0 = 1.0 / y0 + 3.0 * z_grounding_to_;
+                y0 = 1.0 / z0;
+                param0.ytt() = y0;
+            }
         }
         // ZN*
         // Zero sequence impedance of zigzag winding is approximately 10% of positive sequence impedance

power_grid_model_c/power_grid_model_c/include/power_grid_model_c/dataset_definitions.h

@@ -122,6 +122,8 @@ PGM_API extern PGM_MetaAttribute const* const PGM_def_input_transformer_uk;
 PGM_API extern PGM_MetaAttribute const* const PGM_def_input_transformer_pk;
 PGM_API extern PGM_MetaAttribute const* const PGM_def_input_transformer_i0;
 PGM_API extern PGM_MetaAttribute const* const PGM_def_input_transformer_p0;
+PGM_API extern PGM_MetaAttribute const* const PGM_def_input_transformer_i0_zero_sequence;
+PGM_API extern PGM_MetaAttribute const* const PGM_def_input_transformer_p0_zero_sequence;
 PGM_API extern PGM_MetaAttribute const* const PGM_def_input_transformer_winding_from;
 PGM_API extern PGM_MetaAttribute const* const PGM_def_input_transformer_winding_to;
 PGM_API extern PGM_MetaAttribute const* const PGM_def_input_transformer_clock;

power_grid_model_c/power_grid_model_c/src/dataset_definitions.cpp

@@ -111,6 +111,8 @@ PGM_MetaAttribute const* const PGM_def_input_transformer_uk = PGM_meta_get_attri
 PGM_MetaAttribute const* const PGM_def_input_transformer_pk = PGM_meta_get_attribute_by_name(nullptr, "input", "transformer", "pk");
 PGM_MetaAttribute const* const PGM_def_input_transformer_i0 = PGM_meta_get_attribute_by_name(nullptr, "input", "transformer", "i0");
 PGM_MetaAttribute const* const PGM_def_input_transformer_p0 = PGM_meta_get_attribute_by_name(nullptr, "input", "transformer", "p0");
+PGM_MetaAttribute const* const PGM_def_input_transformer_i0_zero_sequence = PGM_meta_get_attribute_by_name(nullptr, "input", "transformer", "i0_zero_sequence");
+PGM_MetaAttribute const* const PGM_def_input_transformer_p0_zero_sequence = PGM_meta_get_attribute_by_name(nullptr, "input", "transformer", "p0_zero_sequence");
 PGM_MetaAttribute const* const PGM_def_input_transformer_winding_from = PGM_meta_get_attribute_by_name(nullptr, "input", "transformer", "winding_from");
 PGM_MetaAttribute const* const PGM_def_input_transformer_winding_to = PGM_meta_get_attribute_by_name(nullptr, "input", "transformer", "winding_to");
 PGM_MetaAttribute const* const PGM_def_input_transformer_clock = PGM_meta_get_attribute_by_name(nullptr, "input", "transformer", "clock");

tests/cpp_unit_tests/test_transformer.cpp

@@ -418,6 +418,7 @@ TEST_CASE("Test Transfomer - Test grounding - Dyn11") {
                                  .pk = 100e3,
                                  .i0 = 0.015,
                                  .p0 = 30.0e4,
+                                 .i0_zero_sequence = 1.0,
                                  .winding_from = WindingType::delta,
                                  .winding_to = WindingType::wye_n,
                                  .clock = 11,
@@ -462,6 +463,18 @@ TEST_CASE("Test Transfomer - Test grounding - Dyn11") {
     }
     DoubleComplex const y_1_shunt = (y_shunt_real + 1i * y_shunt_imag) / base_y_to;
 
+    double const p0_zero_sequence =
+        input.p0 + input.pk * (input.i0_zero_sequence * input.i0_zero_sequence - input.i0 * input.i0);
+    double const y0_shunt_abs = input.i0_zero_sequence * input.sn / input.u2 / input.u2;
+    double const y0_shunt_real = p0_zero_sequence / input.u2 / input.u2;
+    double y0_shunt_imag;
+    if (y0_shunt_real > y0_shunt_abs) {
+        y0_shunt_imag = 0.0;
+    } else {
+        y0_shunt_imag = -std::sqrt(y0_shunt_abs * y0_shunt_abs - y0_shunt_real * y0_shunt_real);
+    }
+    DoubleComplex const y_0_shunt = (y0_shunt_real + 1i * y0_shunt_imag) / base_y_to;
+
     DoubleComplex const tap_ratio_1 = k * std::exp(1.0i * (deg_30 * input.clock));
 
     DoubleComplex const y_1_tt = (1.0 / z_1_series) + 0.5 * y_1_shunt;
@@ -485,7 +498,7 @@ TEST_CASE("Test Transfomer - Test grounding - Dyn11") {
     DoubleComplex const y_0_ff = low_admittance;
     DoubleComplex const y_0_ft = 0.0;
     DoubleComplex const y_0_tf = 0.0;
-    DoubleComplex const y_0_tt = (1.0 / (z_1_series + 3.0 * z_grounding_to)) + y_1_shunt;
+    DoubleComplex const y_0_tt = (1.0 / (1.0 / (1.0 / z_1_series + y_0_shunt) + 3.0 * z_grounding_to));
 
     // Sequence admittances -> phase addmitance
     ComplexTensor<asymmetric_t> y_ff_diagonal;
@@ -707,7 +720,7 @@ TEST_CASE("Test Transformer - Dyn11 - tap_max and tap_min flipped") {
     DoubleComplex const y_0_ff = low_admittance;
     DoubleComplex const y_0_ft = 0.0;
     DoubleComplex const y_0_tf = 0.0;
-    DoubleComplex const y_0_tt = (1.0 / (z_1_series + 3.0 * z_grounding_to)) + y_1_shunt;
+    DoubleComplex const y_0_tt = (1.0 / (1.0 / (1.0 / z_1_series + y_1_shunt) + 3.0 * z_grounding_to));
 
     // Sequence admittances -> phase addmitance
     ComplexTensor<asymmetric_t> y_ff_diagonal;

tests/data/power_flow/zero_sequence_yyn_transformer/asym_output.json

@@ -0,0 +1,21 @@
+{
+  "version": "1.0",
+  "type": "asym_output",
+  "is_batch": false,
+  "attributes": {},
+  "data": {
+    "node": [
+      {"id": 1, "energized": 1, "u_pu": [1, 1, 1], "u": [5773.5026918962576, 5773.5026918962576, 5773.5026918962576], "u_angle": [-2.6439876230271684e-44, -2.0943951023931953, 2.0943951023931953], "p": [881329.20767572254, 1037364.8940881646, 1081305.8982361141], "q": [226203.41305523555, 85377.466355255776, 290921.30803535768]},
+      {"id": 2, "energized": 1, "u_pu": [1.0141430605715689, 0.96656989153700335, 1.0110998975479204], "u": [5855.1576901778626, 5580.4938706947623, 5837.5879802689487], "u_angle": [0.001939888439272864, -2.1247815189679589, 2.0324414694676975], "p": [-799999.99999999849, -1000000.0000000008, -1199999.9999999972], "q": [7.4587508688682828e-10, 1.7304176386161853e-09, 5.1725865757671641e-09]}
+    ],
+    "transformer": [
+      {"id": 4, "energized": 1, "loading": 0.30705253328631577, "p_from": [881329.20767572254, 1037364.8940881646, 1081305.8982361141], "q_from": [226203.41305523555, 85377.466355255761, 290921.30803535768], "i_from": [157.59846030784894, 180.28437829560764, 193.9477496411763], "s_from": [909895.13482607121, 1040872.3433965337, 1119757.8546405528], "p_to": [-799999.99999999849, -1000000.0000000008, -1199999.9999999972], "q_to": [7.4606055156303146e-10, 1.7076212935884444e-09, 5.1782885887484745e-09], "i_to": [136.63167455626581, 179.19560941574917, 205.56435364331952], "s_to": [799999.99999999849, 1000000.0000000008, 1199999.9999999972]}
+    ],
+    "asym_load": [
+      {"id": 6, "energized": 1, "p": [799999.99999999988, 1000000, 1200000], "q": [0, 0, 0], "i": [136.63167455626603, 179.19560941574903, 205.56435364331998], "s": [799999.99999999988, 1000000, 1200000], "pf": [1, 1, 1]}
+    ],
+    "source": [
+      {"id": 7, "energized": 1, "p": [881329.20767572254, 1037364.8940881646, 1081305.8982361141], "q": [226203.41305523555, 85377.466355255761, 290921.30803535774], "i": [157.59846030784894, 180.28437829560764, 193.9477496411763], "s": [909895.13482607121, 1040872.3433965337, 1119757.8546405528], "pf": [0.96860525344405857, 0.99663027908213619, 0.96566047181979187]}
+    ]
+  }
+}

tests/data/power_flow/zero_sequence_yyn_transformer/asym_output.json.license

@@ -0,0 +1,3 @@
+SPDX-FileCopyrightText: Contributors to the Power Grid Model project <powergridmodel@lfenergy.org>
+
+SPDX-License-Identifier: MPL-2.0

tests/data/power_flow/zero_sequence_yyn_transformer/input.json

@@ -0,0 +1,21 @@
+{
+  "version": "1.0",
+  "type": "input",
+  "is_batch": false,
+  "attributes": {},
+  "data": {
+    "node": [
+      {"id": 1, "u_rated": 10000},
+      {"id": 2, "u_rated": 10000}
+    ],
+    "transformer": [
+      {"id": 4, "from_node": 1, "to_node": 2, "from_status": 1, "to_status": 1, "u1": 10000, "u2": 10000, "sn": 10000000, "uk": 0.10000000000000001, "pk": 0, "i0": 0.050000000000000003, "p0": 0, "i0_zero_sequence": 1, "winding_from": 0, "winding_to": 1, "clock": 12, "tap_side": 0, "tap_pos": 0, "tap_min": 0, "tap_max": 0, "tap_nom": 0, "tap_size": 0}
+    ],
+    "asym_load": [
+      {"id": 6, "node": 2, "status": 1, "type": 0, "p_specified": [800000, 1000000, 1200000], "q_specified": [0, 0, 0]}
+    ],
+    "source": [
+      {"id": 7, "node": 1, "status": 1, "u_ref": 1, "sk": 1.0000000000000001e+50, "rx_ratio": 0, "z01_ratio": 1}
+    ]
+  }
+}

tests/data/power_flow/zero_sequence_yyn_transformer/input.json.license

@@ -0,0 +1,3 @@
+SPDX-FileCopyrightText: Contributors to the Power Grid Model project <powergridmodel@lfenergy.org>
+
+SPDX-License-Identifier: MPL-2.0

tests/data/power_flow/zero_sequence_yyn_transformer/params.json

@@ -0,0 +1,7 @@
+{
+  "calculation_method": [
+    "newton_raphson"
+  ],
+  "rtol": 1e-07,
+  "atol": 1e-07
+}

tests/data/power_flow/zero_sequence_yyn_transformer/params.json.license

@@ -0,0 +1,3 @@
+SPDX-FileCopyrightText: Contributors to the Power Grid Model project <powergridmodel@lfenergy.org>
+
+SPDX-License-Identifier: MPL-2.0