[FEATURE] Provide a generic_branch

[petersalemink95]

Following the discussion in #642 it was concluded that the implementation of a generic branch could be an addition to Power Grid Model.

We can make a new component generic_branch, which models branches as a Pi model; i.e. the user can specify:

• k: transformer ratio (float)
• theta0 / theta1 / theta2: angle shift
• r0_series / r1_series / r2_series
• x0_series / x1_series / x2_series
• g0_shunt / g1_shunt / g2_shunt
• b0_shunt / b1_shunt / b2_shunt

Note: except for k, all attributes are given for zero, positive and negative sequence.

##steps to follow
In order to add a new component to Power Grid Model the following steps need to be taken:

• Check if the available input/update/output data suffices (check power-grid-model/power_grid_model_c/power_grid_model/include/power_grid_model/auxiliary/ input.hpp / update.hpp / output.hpp or in the documentation directly)
• If not, add the new data format to code_generation/data/attribute_classes/ input.json / update.json / output.json + run code_generation/code_gen.py
• Create a new component in a new power-grid-model/power_grid_model_c/power_grid_model/include/power_grid_model/component/[component].hpp file that at least inherits from Base, but in this case GenericBranch should inherit from Branch
• If necessary: add new enums or exceptions
• Create the necessary unit tests in power-grid-model/tests/cpp_unit_tests/test_[component].cpp
• Add the test_[component].cpp to power-grid-model/tests/cpp_unit_tests/CMakeLists.txt
• Add component to power_grid_model_c/power_grid_model/include/power_grid_model/all_components.hpp
• Not necessary for this component (If necessary update main_core/topology.hpp / input.hpp / output.hpp / update.hpp)
• Add component to code_generation/data/dataset_class_maps/dataset_definitions.json + re-run code_generation/code_gen.py
• Add validation test cases to tests/data
• Update input/update data validator for the new component: src/power_grid_model/validation/validation.py + add corresponding tests

Note: the order is recommended, but not necessary

[TonyXiang8787]

@petersalemink95

The user needs to be able to provide the attributes of r, x, g, b, k, theta for positive, negative and zero sequence. So the attribute list is actually 3 times as described.

[sudo-ac]

Hi. I have now created a PI model Transformer in my fork. However, when implementing the asymmetric part, I stuck to the formulation of the transformer (i.e.
r_grounding_from, x_grounding_from,.r_grounding_to, x_rounding_to). Would that also be ok? Can you have a look at the class in my fork?

[TonyXiang8787]

Hi @sudo-ac, @petersalemink95,

Now after re-think about the issue maybe the way that we define attribute as r/x_series and g/b_shunt is not very logical. The grounding impedance of the transformer will be appended in one of the legs in the pi model. Making the shunt part of the pi model not the same between from and to side. In the mathematical model the branch is described as:

$$ \begin{bmatrix} I_f \\ I_t \end{bmatrix} = \begin{bmatrix} Y_{ff} & Y_{ft} \\ Y_{tf} & Y_{tt} \end{bmatrix} \begin{bmatrix} U_f \\ U_t \end{bmatrix} $$

The above formula applies three time for positive, negative and zero sequence. Using for complex parameters $Y_{ff}$, $Y_{ft}$, $Y_{tf}$, $Y_{tt}$ we can take the tap ratio, phase shift, grounding, and everything into account.

Maybe we need to let user directly define those instead of r/x_series and g/b_shunt, which is a half-way solution.

[sudo-ac]

Hi Tony (@petersalemink95),

this means you prefer to define the input parameters in this way:

r0_series / r1_series / r2_series
x0_series / x1_series / x2_series
g0_shunt / g1_shunt / g2_shunt
b0_shunt / b1_shunt / b2_shunt

0 = zero sysytem
1 = positiv system
2 = negative system

[TonyXiang8787]

Hi Tony (@petersalemink95),

this means you prefer to define the input parameters in this way:

r0_series / r1_series / r2_series x0_series / x1_series / x2_series g0_shunt / g1_shunt / g2_shunt b0_shunt / b1_shunt / b2_shunt

0 = zero sysytem 1 = positiv system 2 = negative system

Hi @sudo-ac,

Actually, not after my rethinking. The current way (in the issue) we define series impedance and shunt admittance does not allow us to give different values for the two legs of the pi model, which is needed for the transformer modelling. My new idea would be asking user to provide attributes which are closer to the internal mathematical model, namely yff_1, yft_1, ytf_1, and ytt_1. Same for _2 and _0. But it is a lot of attributes. Note in the interface we do not support complex numbers, so it eventually becomes gff_1 and bff_1, etc. This becomes 8 attributes per sequence and 24 attributes in total.

@petersalemink95 what do you think?

[petersalemink95]

@TonyXiang8787 The idea of generic_branch was to make the life of the end user easier, by not having to do conversions of the data themselves if R/X/G/B are present. I think letting the end user calculate all ff/ft/tf/tt parameters does not fit the original purpose of the issue anymore...

What do you think @sudo-ac?

[TonyXiang8787]

But now we are running into issue that the transformer has not the same values in the two legs of the pi model. Especially for zero-sequence the two legs of the pi model are almost always different.

We can still use the R/X/G/B model, then we have to use like g1_shunt_from and g1_shunt_to.

[sudo-ac]

confussing, i just added it the following way:

 explicit GenericBranch(GenericBranchInput const& genericbranch_input, double u1_rated, double u2_rated)
        : Branch{genericbranch_input},
          sn_{genericbranch_input.sn},
          r1_{genericbranch_input.r1},
          x1_{genericbranch_input.x1},
          g1_{genericbranch_input.g1},
          b1_{genericbranch_input.b1},          
          ratio_{is_nan(genericbranch_input.ratio) ? 1.0 : genericbranch_input.ratio},
          shift_{is_nan(genericbranch_input.shift) ? 0.0 : fmod(genericbranch_input.shift, 2 * M_PI)},
          base_i_from_{base_power_3p / u1_rated / sqrt3},
          base_i_to_{base_power_3p / u2_rated / sqrt3},
          r0_{is_nan(genericbranch_input.r0) ? 0.0 : genericbranch_input.r0},
          x0_{is_nan(genericbranch_input.x0) ? 0.0 : genericbranch_input.x0},
          g0_{is_nan(genericbranch_input.g0) ? 0.0 : genericbranch_input.g0},
          b0_{is_nan(genericbranch_input.b0) ? 0.0 : genericbranch_input.b0},
          r2_{is_nan(genericbranch_input.r2) ? genericbranch_input.r0 : genericbranch_input.r2},
          x2_{is_nan(genericbranch_input.x2) ? genericbranch_input.x0 : genericbranch_input.x2},
          g2_{is_nan(genericbranch_input.g2) ? genericbranch_input.g0 : genericbranch_input.g2},
          b2_{is_nan(genericbranch_input.b2) ? genericbranch_input.b0 : genericbranch_input.b2} {

        double const base_y_to = base_i_to_ * base_i_to_ / base_power_1p;
        y1_series_ = 1.0 / (r1_ + 1.0i * x1_) / base_y_to;
        y1_shunt_ = (g1_ + 1.0i * b1_) / base_y_to;

        // Add new calculations for r0, x0, g0, b0 (zero-sequence)        
        DoubleComplex const z0_series = 1.0 / y1_series_ + 3.0 * ((r0_ + 1.0i * x0_) / ratio_ / ratio_);
        y0_series_ = 1.0 / z0_series / base_y_to;
        y0_shunt_ = (g0_ + 1.0i * b0_) / base_y_to;

        // Add new calculations for r2, x2, g2, b2 (negative-sequence)
        // Check for zero impedance before dividing
        if (r2_ == 0.0 && x2_ == 0.0) {
           y2_series_ = 0.0;  // Set to 0 or a very large value to avoid division by 0
        } else {    
            y2_series_ = 1.0 / (r2_ + 1.0i * x2_) / base_y_to;
        }    
        y2_shunt_ = (g2_ + 1.0i * b2_) / base_y_to;
    }

and:

  BranchCalcParam<asymmetric_t> asym_calc_param() const override {
        // Positive sequence (param1)
        auto const param1 = calc_param_y_sym(y1_series_, y1_shunt_, ratio_ * std::exp(1.0i * shift_));
        // Negative sequence (param2) using r2, x2, g2, b2
        auto const param2 = calc_param_y_sym(y2_series_, y2_shunt_, ratio_ * std::exp(1.0i * (-shift_)));
        // Zero sequence (param0) using r0, x0, g0, b0
        BranchCalcParam<symmetric_t> param0{};
        param0 = calc_param_y_sym(y0_series_, y0_shunt_, ratio_ * std::exp(1.0i * shift_));
        // Calculation of the symmetrical matrices
        ComplexTensor<asymmetric_t> const sym_matrix = get_sym_matrix();
        ComplexTensor<asymmetric_t> const sym_matrix_inv = get_sym_matrix_inv();
        // Initialize the asymmetric parameter (Yabc) based on the sequence parameters
        BranchCalcParam<asymmetric_t> param;
        for (size_t i = 0; i != 4; ++i) {
            // Yabc = A * Y012 * A^-1
            ComplexTensor<asymmetric_t> y012;
            y012 << param0.value[i], 0.0, 0.0, 0.0, param1.value[i], 0.0, 0.0, 0.0, param2.value[i];
            param.value[i] = dot(sym_matrix, y012, sym_matrix_inv);
        }

        return param;
    }

Would this confirm to your requirements?

[sudo-ac]

Oh I see, theta0 and theta2 is missing, also I have to rename shift and ratio....

[petersalemink95]

@sudo-ac We're also thinking hard about the best solution.
The asymmetric calculation is a bit harder than expected here.

To include you in our thinking:

• For transformers there is a deviation between the from and to side of b0 and g0 (zero sequence). Therefore the original proposed idea is not valid.
• Adding r_grounding_from, x_grounding_from, r_grounding_to, x_rounding_to as you proposed can be a solution for transformers, but for transformers only; making it not a generic_branch anymore.
• Adding y_ff, y_ft, y_tf and y_tt , as @TonyXiang8787 proposed, can be a solution, but this would put more work at the user for asymmetric calculations

At this moment, do you just want to do symmetric calculations? If that is the case we can also implement symmetric calculations first and have some more time to think about the best asymmetric implementation: throw a NotImplementedError for now.

We have a brainstorm with the team tomorrow about this topic. If you have any ideas on how you would like to have this implemented from a user perspective, please let us know!

[sudo-ac]

Hi Peter, thanx to your answer. I see the problem in the transformer class for asymetrical calculation for the different winding types. For the moment my implementation looks like this:

// SPDX-FileCopyrightText: Contributors to the Power Grid Model project <[email protected]>
//
// SPDX-License-Identifier: MPL-2.0

#pragma once

#include "branch.hpp"

#include "../auxiliary/input.hpp"
#include "../auxiliary/output.hpp"
#include "../auxiliary/update.hpp"
#include "../calculation_parameters.hpp"
#include "../common/common.hpp"
#include "../common/exception.hpp"
#include "../common/three_phase_tensor.hpp"

/*
    generic Branch: either a line N = 1 or a transformoer N = t*e^(j*theta1)
    paramaters should be given as r1, x1, ....

     -----| |-----------y1_series-------
          | |   |                 |
          | |   y1_shunt          y1_shunt
          | |   |                 |
          | |   |                 |
     -----| |--------------------------
          N = k * e^(j*theta1)

*/

namespace power_grid_model {

class GenericBranch final : public Branch {
  public:
    using InputType = GenericBranchInput;
    using UpdateType = BranchUpdate;
    static constexpr char const* name = "generic_branch";

    explicit GenericBranch(GenericBranchInput const& genericbranch_input, double u1_rated, double u2_rated)
        : Branch{genericbranch_input},
          sn_{genericbranch_input.sn},
          r1_{genericbranch_input.r1},
          x1_{genericbranch_input.x1},
          g1_{genericbranch_input.g1},
          b1_{genericbranch_input.b1},          
          k_{is_nan(genericbranch_input.k) ? 1.0 : genericbranch_input.k},          
          theta1_{is_nan(genericbranch_input.theta1) ? 0.0 : fmod(genericbranch_input.theta1, 2 * M_PI)},  
          theta0_{is_nan(genericbranch_input.theta0) ? theta1_ : fmod(genericbranch_input.theta0, 2 * M_PI)},
          theta2_{is_nan(genericbranch_input.theta2) ? -1.0*theta1_ : fmod(genericbranch_input.theta2, 2 * M_PI)},          
          base_i_from_{base_power_3p / u1_rated / sqrt3},
          base_i_to_{base_power_3p / u2_rated / sqrt3},
          r0_{is_nan(genericbranch_input.r0) ? 0.0 : genericbranch_input.r0},
          x0_{is_nan(genericbranch_input.x0) ? 0.0 : genericbranch_input.x0},
          g0_{is_nan(genericbranch_input.g0) ? 0.0 : genericbranch_input.g0},
          b0_{is_nan(genericbranch_input.b0) ? 0.0 : genericbranch_input.b0},
          r2_{is_nan(genericbranch_input.r2) ? genericbranch_input.r0 : genericbranch_input.r2},
          x2_{is_nan(genericbranch_input.x2) ? genericbranch_input.x0 : genericbranch_input.x2},
          g2_{is_nan(genericbranch_input.g2) ? genericbranch_input.g0 : genericbranch_input.g2},
          b2_{is_nan(genericbranch_input.b2) ? genericbranch_input.b0 : genericbranch_input.b2} {

        double const base_y_from = base_i_from_ / (u1_rated / sqrt3);
        double const base_y_to = base_i_to_ / (u2_rated / sqrt3);

        y1_series_ = 1.0 / (r1_ + 1.0i * x1_) / base_y_to;
        y1_shunt_ = (g1_ + 1.0i * b1_) / base_y_to;

        // Add new calculations for r0, x0, g0, b0 (zero-sequence)  
        DoubleComplex const z0 = (r0_ + 1i * x0_) * base_y_from;
        DoubleComplex const z2 = (r2_ + 1i * x2_) * base_y_to;
        
        DoubleComplex const z2_series = (r2_ + 1i * x2_) * base_y_to;
        DoubleComplex const z1_series = (r1_ + 1i * x1_) * base_y_to;
        DoubleComplex const z0_series = z1_series + 3.0 * (z0 + z2 / k_ / k_);

        
        y0_series_ = 1.0 / z0_series;
        y0_shunt_ = (g0_ + 1.0i * b0_) / base_y_to;
        y2_series_ = 1.0 / z2_series;
        y2_shunt_ = (g2_ + 1.0i * b2_) / base_y_to;
    }

    // override getter
    double base_i_from() const override { return base_i_from_; }
    double base_i_to() const override { return base_i_to_; }
    double loading(double max_s, double /*max_i*/) const override { return is_nan(sn_) ? 0.0 : (max_s / sn_); };
    double phase_shift() const override { return theta1_; }
    bool is_param_mutable() const override { return false; }

    template <symmetry_tag sym>
    BranchOutput<sym> get_output(BranchSolverOutput<sym> const& branch_solver_output) const {
        // result object
        BranchOutput<sym> output{};
        static_cast<BaseOutput&>(output) = base_output(true);
        // calculate result
        output.p_from = base_power<sym> * real(branch_solver_output.s_f);
        output.q_from = base_power<sym> * imag(branch_solver_output.s_f);
        output.i_from = base_i_from() * cabs(branch_solver_output.i_f);
        output.s_from = base_power<sym> * cabs(branch_solver_output.s_f);
        output.p_to = base_power<sym> * real(branch_solver_output.s_t);
        output.q_to = base_power<sym> * imag(branch_solver_output.s_t);
        output.i_to = base_i_to() * cabs(branch_solver_output.i_t);
        output.s_to = base_power<sym> * cabs(branch_solver_output.s_t);
        double const max_s = std::max(sum_val(output.s_from), sum_val(output.s_to));
        double const max_i = std::max(max_val(output.i_from), max_val(output.i_to));
        output.loading = loading(max_s, max_i);
        output.shift = theta1_;
        output.ratio = k_;

        output.r1 = r1_;
        output.x1 = x1_;
        output.b1 = b1_;
        output.g1 = g1_;

        return output;
    }

  private:
    double sn_;
    double r1_, x1_, g1_, b1_;
    double k_;
    double theta1_, theta0_, theta2_;
    double base_i_from_, base_i_to_;    
    double r0_, x0_, g0_, b0_;
    double r2_, x2_, g2_, b2_;

    DoubleComplex y0_series_;
    DoubleComplex y0_shunt_;
    DoubleComplex y1_series_;
    DoubleComplex y1_shunt_;
    DoubleComplex y2_series_;
    DoubleComplex y2_shunt_;


    BranchCalcParam<symmetric_t> sym_calc_param() const override {
        return calc_param_y_sym(y1_series_, y1_shunt_, k_ * std::exp(1.0i * theta1_));
    }

    BranchCalcParam<asymmetric_t> asym_calc_param() const override {
        // Positive sequence (param1)
        auto const param1 = calc_param_y_sym(y1_series_, y1_shunt_, k_ * std::exp(1.0i * theta1_));
        // Negative sequence (param2) using r2, x2, g2, b2
        auto const param2 = calc_param_y_sym(y2_series_, y2_shunt_, k_ * std::exp(1.0i * (theta2_)));
        // Zero sequence (param0) using r0, x0, g0, b0
        BranchCalcParam<symmetric_t> param0{};
        param0 = calc_param_y_sym(y0_series_, y0_shunt_, k_ * std::exp(1.0i * theta0_));
        // Calculation of the symmetrical matrices
        ComplexTensor<asymmetric_t> const sym_matrix = get_sym_matrix();
        ComplexTensor<asymmetric_t> const sym_matrix_inv = get_sym_matrix_inv();
        // Initialize the asymmetric parameter (Yabc) based on the sequence parameters
        BranchCalcParam<asymmetric_t> param;
        for (size_t i = 0; i != 4; ++i) {
            // Yabc = A * Y012 * A^-1
            ComplexTensor<asymmetric_t> y012;
            y012 << param0.value[i], 0.0, 0.0, 0.0, param1.value[i], 0.0, 0.0, 0.0, param2.value[i];
            param.value[i] = dot(sym_matrix, y012, sym_matrix_inv);
        }

        return param;
    }
};

} // namespace power_grid_model

I need only the symetrical calculation. So, let me know your ideas.

[sudo-ac]

maybe the user could give an optional type parameter for the asymetrical calculation?

[sudo-ac]

....by the way, the section in methode BranchOutput is just for our debugging purposes and should not be contributed....

[petersalemink95]

Thanks @sudo-ac. I'll come back to you after the brainstorm tomorrow!

[sudo-ac]

....by the way, the section in methode BranchOutput is just for debugging purposes and should not be

[sudo-ac]

...to give y_ff, y_ft, etc. as input params, i could be possible to give them in polar notation (Magnitude + Angle)

[sudo-ac]

Hi @petersalemink95, @TonyXiang8787
we discussed the problem, and I understand that there is no generic way to calculate BranchCalcParam<asymmetric_t>. Therefore, Tony's suggestion could provide a generic solution to this issue. We could add it to the input.json.

  {

            "name": "GenericBranchInput",
            "base": "BranchInput",
            "attributes": [
                {    
                    "data_type": "double",
                    "names": ["r1","x1","g1","b1"],
                    "description": "positive sequence parameters"
                },
                {
                    "data_type": "double",
                    "names": "k",
                    "description": "transformer ratio, default = 1.0"
                },
                {
                    "data_type": "double",
                    "names": ["theta"],
                    "description": "angle shift"
                },
                {
                    "data_type": "double",
                    "names": "sn",
                    "description": "rated power for calculation of loading (optional)"
                },
               {

                    "data_type": "double",
                    "names": ["yff_mag","yff_angle"],
                    "description": "complex yff in polar ....."
                },
                 {
                    "data_type": "double",
                    "names": ["yft_mag","yft_angle"],
                    "description": "complex yft in polar ....."
                },
                {
                    "data_type": "double",
                    "names": ["ytf_mag","ytf_angle"],
                    "description": "complex ytf in polar ....."
                },
                {
                    "data_type": "double",
                    "names": ["ytt_mag","ytt_angle"],
                    "description": "complex ytt in polar ....."
                },

The user can calculate these parameters on their own, for example, by using a Python script. You could perhaps provide a sample script. The GenericBranch returns these parameters in the function BranchCalcParam<asymmetric_t>.

[sudo-ac]

...or ytt... could be given in cartesian form ytt_re, ytt_img.....

[sudo-ac]

Oh, I see, so you would have to pass a 3x3 matrix for each of the parameters (yff...ytt)? Then of course the above would not be correct.

[petersalemink95]

Hi @sudo-ac thanks for thinking along!

We've had a brainstorm with the team and concluded that the asymmetric part is quite difficult and needs some more work.
Since you're only interested in doing symmetric calculations at the moment you can go ahead by only providing the symmetric calculation (and positive sequence r,x,g,b). For the asymmetric calculation you can throw a NotImplementedError. I hope this idea works for you as well!

I will make a new issue for the asymmetric calculation, including some of the options we discussed.

[sudo-ac]

Hi @petersalemink95,

Thank you for your answer. We are currently discussing a simple solution. I'll get back to you when we have formulated the idea.

[sudo-ac]

Hi @petersalemink95

from our perspective, our work on the generic_branch (#729) is "code complete". What remains is updating the documentation, and adding the ability to update the ratio and theta could be done in the future. My colleague will start working on #739 next week. I’ll let you know when we are done with that.

[petersalemink95]

Hi @sudo-ac,

That is great to hear! Once the PR is made we'll do the review.
Regarding #739: there is no decision yet on how to implement asymmetric calculations. For us it would be fine to go ahead by just merging the symmetrical generic_branch for now; if that would suit your needs.