Implementation of pandapower net impedance element

[rbolgaryn]

Is your feature request related to a problem? Please describe.

Dear Benjamin, we just implemented a new feature in pandapower for unsymmetric (different for from and to bus) shunt admittances. Earlier, we relied on the column BR_B for the shunt admittance, where there was a complex value to represent B and G. Actually, the values in this column were:

b = 2 * net.f_hz * math.pi * line["c_nf_per_km"].values * 1e-9 * baseR * length_km * parallel
g = line["g_us_per_km"].values * 1e-6 * baseR * length_km * parallel
branch[f:t, BR_B] = b - g * 1j

You can notice that the b here is defined as b-g1j, instead of the correct way g + b1j, because at a later point it was multiplied by 1j in makeYbus, where it became correct g + b*1j.

We made it explicit now in the PR e2nIEE/pandapower#2347

This also allows us creating an impedance element with unequal shunt admittances, and also transformer elements with unequal leakage impedances.

In the case 118, there are some branches that connect different voltage levels and have positive shunt reactances B. When imported from ppc earlier, we used to represent them as transfromers. However, this would be incorrect because transformers are inductive elements and they only can have negative B. To rectify this, we are importing them as impedance elements now, because they can have positive B and still connect different voltage levels (https://github.com/pawellytaev/pandapower/blob/1f40fa0192e997d694c3fc620ac1e5ab83c47ce8/pandapower/converter/pypower/from_ppc.py#L380).

Due to the fact that the impedance elements are not converted from pandapower in init at the moment, we now have failing tests in the contingency function that relies on the lightsim2grid model.

Describe the solution you'd like

An addition of impedance element in the pandapower -> lightsim2grid init function, that supports unequal values for rft and rtf, xft and xtf, gf and dt, bf and bt. This would actually be a very versatile element because of this possibility of unequal impedances.

Describe alternatives you've considered

We can take the case118 out of the test coverage temporarily, to be able to have passing tests.

Additional context

Here relevant parts of the code in pandapower:

Calculating the impedance parameters for ppc:

https://github.com/pawellytaev/pandapower/blob/1f40fa0192e997d694c3fc620ac1e5ab83c47ce8/pandapower/build_branch.py#L614

Build Ybus matrix vectors while considering unequal from and to values:
https://github.com/pawellytaev/pandapower/blob/1f40fa0192e997d694c3fc620ac1e5ab83c47ce8/pandapower/pypower/makeYbus.py#L96

What do you think about this?

Best regards,
Roman

[BDonnot]

Hello,

TL;DR: please send me a piece of code as small as possible so I can check if the fix I have in mind works. It will not require too many changes in lightsim2grid because it can already handle different shunt admittance internally (just not in the pandapower converter at the moment).

Long detailed answer:

This also allows us creating an impedance element with unequal shunt admittances, and also transformer elements with unequal leakage impedances.

This is a nice addition and much needed for real powergrid cases.

Lightsim2grid can already handle this internally. I already use this when initializing the gridmodel from pypowsybl for example with:

    gridmodel.init_powerlines_full(line_r,
                               line_x,
                               line_h_or,
                               line_h_ex,
                               lor_bus,
                               lex_bus
                              )

In this case line_h_or are complex vector equal to (b_or + 1j * g_or)

Right now, when initializing a gridmodel from pandapower, what is done (assuming parrallel is 1) :

            pp_net.line["r_ohm_per_km"].values * pp_net.line["length_km"].values,
            pp_net.line["x_ohm_per_km"].values * pp_net.line["length_km"].values,
            pp_net.line["c_nf_per_km"].values * pp_net.line["length_km"].values,
            pp_net.line["g_us_per_km"].values * pp_net.line["length_km"].values,
            pp_net.bus.loc[pp_net.line["from_bus"]]["vn_kv"], 
            pp_net.bus.loc[pp_net.line["to_bus"]]["vn_kv"], 

Internally, lightsim2grid will compute the correct g (based on the c_nf_per_km here) but will discard the g.

Also I am not impacted by the "multiplying by 1j" part because I don't do it in lightsim2grid, so the right value is applied.

In the case 118, there are some branches that connect different voltage levels and have positive shunt reactances B.

You are totally correct, there are some weird things in the 118 especially with the voltages as you pointed out.

To rectify this, we are importing them as impedance elements now,

There are generic formula to pass from a "transformer" (with given ratio and different voltages at both its ends) to powerlines. For example, you can use this one:

    v1 = branch_from_kv
    v2 = branch_to_kv
    line_r = sn_mva *  df_line["r"].values / v1 / v2
    line_x = sn_mva *  df_line["x"].values / v1 / v2
    tmp_ = np.reciprocal(df_line["r"].values + 1j*df_line["x"].values)
    b1 = df_line["b1"].values * v1*v1/sn_mva + (v1-v2)*tmp_.imag*v1/sn_mva
    b2 = df_line["b2"].values * v2*v2/sn_mva + (v2-v1)*tmp_.imag*v2/sn_mva
    g1 = df_line["g1"].values * v1*v1/sn_mva + (v1-v2)*tmp_.real*v1/sn_mva
    g2 = df_line["g2"].values * v2*v2/sn_mva + (v2-v1)*tmp_.real*v2/sn_mva
    line_h_or = (b1 + 1j * g1)
    line_h_ex = (b2 + 1j * g2)

This is the formula used in lightsim2grid and it comes from industrial software "powsybl", for example you can find it use here:
https://github.com/powsybl/powsybl-core/blob/266442cbbd84f630acf786018618eaa3d496c6ba/ieee-cdf/ieee-cdf-converter/src/main/java/com/powsybl/ieeecdf/converter/IeeeCdfImporter.java#L347

Due to the fact that the impedance elements are not converted from pandapower in init at the moment, we now have failing tests in the contingency function that relies on the lightsim2grid model.

I totally understand the whys. The way the "init_from_pandapower" works at the moment will fail for the ieee118 because of the error you pointed out.

If you provide me with the failing tests, I can try to work on it.

It will be rather easy, nothing has to change deeply in lightsim2grid. Only this function needs to be modified:
https://github.com/BDonnot/lightsim2grid/blob/master/src/DataConverter.cpp#L126

The rest of lightsim2grid will behave normally, because as I said internally the cpp part can handle different h1 and h2.

A quick note, i'll leave for holiday for 3 weeks. Not sure i'll be able to fix it and release a new lightsim2grid version before leaving.

Best

Benjamin

[rbolgaryn]

Dear Benjamin,

Thank you for the quick response! It sounds like lightsim2grid can integrate this functionality!

You can find the failing tests here (I added one separately for case 118):

https://github.com/pawellytaev/pandapower/blob/5abbea2573a85e7fd3e2d21f3ea6c1aee5ffa3cf/pandapower/test/contingency/test_contingency.py#L136

Also a test for a grid with unequal transformer impedances (I assume it would maybe require a bit more adjustments):

https://github.com/pawellytaev/pandapower/blob/da4b5ae6acd42e75dd37ef20d4dcbd823fef48d3/pandapower/test/contingency/test_contingency.py#L156

The code for calculating the R, X, B, G for from/to buses is in the function _wye_delta:

https://github.com/pawellytaev/pandapower/blob/5abbea2573a85e7fd3e2d21f3ea6c1aee5ffa3cf/pandapower/build_branch.py#L345

https://github.com/pawellytaev/pandapower/blob/647fca74bf5b02eaebb0d0556de04cf4844e954c/pandapower/test/contingency/test_contingency.py#L53

https://github.com/pawellytaev/pandapower/blob/647fca74bf5b02eaebb0d0556de04cf4844e954c/pandapower/test/contingency/test_contingency.py#L102

As for the timeline, do not worry if you do not have a new release before your vacation, we can set the tests to xfail and reactivate them as soon as the new lightsim2grid release is out.

Best regards,
Roman

[BDonnot]

Hello,

About this issue, I noticed something "interesting".

There is an error which is not related to the internal modelling with the lines having different parameters at both sides.

I don't remember the previous behaviour, but apparently when a contingency leads to a divergence of the powerflow now lightsim2grid returns 0. for the voltage (previously it might have returned Nan, I don't remember).

And as far as I can tell, this is what makes the test for "min" (I did not check the "max" yet) failed: all "min voltages" are 0 because one contigency (at least) cannot be computed due to a divergence of the powerflow.

Moving from

    v_init = net._ppc["internal"]["V"]
    s.compute(v_init, net._options["max_iteration"], net._options["tolerance_mva"])
    v_res = s.get_voltages()
    s.compute_flows()
    kamps_all = s.get_flows()
    vm_pu = np.abs(v_res)

to this: (adding the last line)

    v_init = net._ppc["internal"]["V"]
    s.compute(v_init, net._options["max_iteration"], net._options["tolerance_mva"])
    v_res = s.get_voltages()
    s.compute_flows()
    kamps_all = s.get_flows()
    vm_pu = np.abs(v_res)
   vm_pu[np.abs(vm_pu) <= 1e-6] = np.nan  # voltages are 0. in case of divergence

fixed part of the problem.

Now i'm diving into the grid modelling :-)

[BDonnot]

Hi again :-)

I am looking at the failing test and I found some other difference of behaviour between lightsim2grid and pandapower that have an impact on the contingency analysis.

I made a "unittest" out of the the failing test above that you can find here:

def run_for_from_bus_loading(net, **kwargs):
    # copied from https://github.com/pawellytaev/pandapower/blob/5abbea2573a85e7fd3e2d21f3ea6c1aee5ffa3cf/pandapower/test/contingency/test_contingency.py#L331
    pp.runpp(net, **kwargs)
    print(net.res_bus.iloc[6]["vm_pu"])
    net.res_line["loading_percent"] = net.res_line.i_from_ka / net.line.max_i_ka * 100.
    if len(net.trafo) > 0:
        max_i_ka_limit = net.trafo.sn_mva.values / (net.trafo.vn_hv_kv.values * np.sqrt(3))
        net.res_trafo["loading_percent"] = net.res_trafo.i_hv_ka / max_i_ka_limit * 100.
        
        
class PandapowerIntegration(unittest.TestCase):
    def setUp(self) -> None:
        if version.parse(pp.__version__) < MIN_PP_VERSION:
            self.skipTest(f"Pandapower too old for this test: requires >={MIN_PP_VERSION}, found {pp.__version__}")
        return super().setUp()
    
    def test_case118(self, case=None):
        """test inspired from the test provided in the github issue https://github.com/BDonnot/lightsim2grid/issues/88#issuecomment-2265150641
        linked https://github.com/pawellytaev/pandapower/blob/da4b5ae6acd42e75dd37ef20d4dcbd823fef48d3/pandapower/test/contingency/test_contingency.py#L156"""
        if case is None:
            case = pp.networks.case118()
        net = copy.deepcopy(case)
        net2 = copy.deepcopy(net)
        nminus1_cases = {"line": {"index": net.line.index.values},
                         "trafo": {"index": net.trafo.index.values}}

        res = run_contingency(net2, nminus1_cases, contingency_evaluation_function=run_for_from_bus_loading)

        run_contingency_ls2g(net, nminus1_cases)
        for s in ("min", "max"):
            assert np.allclose(res["bus"][f"{s}_vm_pu"], net.res_bus[f"{s}_vm_pu"].values, atol=1e-9, rtol=0), f'for bus voltages for {s}, max error {np.abs(res["bus"][f"{s}_vm_pu"] - net.res_bus[f"{s}_vm_pu"].values).max():.2e}'
            assert np.allclose(np.nan_to_num(res["line"][f"{s}_loading_percent"]),
                            net.res_line[f"{s}_loading_percent"].values, atol=1e-6, rtol=0), s
            if len(net.trafo) > 0:
                assert np.allclose(np.nan_to_num(res["trafo"][f"{s}_loading_percent"]),
                                net.res_trafo[f"{s}_loading_percent"].values, atol=1e-6, rtol=0), s
           
    def test_case14(self):
        self.test_case118(pp.networks.case14())

And to simplify my analysis, I use the pandapower available on pypi (to make sure lightsim2grid and pandapower have the same interpration of the parameters). I also study the case14 as it is much easier to understand what is going on on small network.

The test fail with a difference of min voltages on bus 6. And after looking contingency per contingency, the results for this bus match perfectly (within tolerance) for all of them except for the penultimate transformer connecting bus 6 to bus 7 (where there is a generator).
If I simulate this transformer with pandapower, the powerflow converges.
But it does not in lightsim2grid. And it does not because the grid would be split in two: bus 1 to 6, 8 to 14 on one side and bus 7 alone on the other.

I suspect there is a similar issue for the case118 where pandapower considers only the "main component" and discard the rest of the grid whereas lightsim2grid consider it's a divergence (user has to manually turn off all elements on the "rest of the grid" for lightsim2grid to converge) - a less error prone behaviour in my opinion.

From there, I can propose different options

1. make available python side which contingencies are simulated by lightsim2grid and for those who are not you fail back in using pandapower for these ones
2. have an option in lightsim2grid to "discard element not in the main component" and run the powerflow in the contingency analysis
3. keep the results as is and remove these cases from the unit test

[BDonnot]

Final update on the topic (waiting for some inputs from your side for the contingency analysis)

I can confirm that when I remove the "diverging powerflow in lightsim2grid due to disconnected grid" then, for the pandapower version on pypi, I get the right behaviour between lightsim2grid and pandapower (with a tolerance of 1e-8)

This is the code I used:

class PandapowerIntegration(unittest.TestCase):
    def setUp(self) -> None:
        if version.parse(pp.__version__) < MIN_PP_VERSION:
            self.skipTest(f"Pandapower too old for this test: requires >={MIN_PP_VERSION}, found {pp.__version__}")
        return super().setUp()
    
    def _aux_test_case(self, case=None, nminus1_cases=None):
        """test inspired from the test provided in the github issue https://github.com/BDonnot/lightsim2grid/issues/88#issuecomment-2265150641
        linked https://github.com/pawellytaev/pandapower/blob/da4b5ae6acd42e75dd37ef20d4dcbd823fef48d3/pandapower/test/contingency/test_contingency.py#L156"""
        if case is None:
            case = pp.networks.case118()
        net = copy.deepcopy(case)
        net2 = copy.deepcopy(net)
        if nminus1_cases is None:
            nminus1_cases = {"line": {"index": net.line.index.values},
                             "trafo": {"index": net.trafo.index.values}}

        res = run_contingency(net2, nminus1_cases, contingency_evaluation_function=run_for_from_bus_loading)

        run_contingency_ls2g(net, nminus1_cases)
        for s in ("min", "max"):
            # import pdb
            # pdb.set_trace()
            # (np.abs(res["bus"][f"{s}_vm_pu"] - net.res_bus[f"{s}_vm_pu"].values) >= 1e-5).nonzero()
            assert np.allclose(res["bus"][f"{s}_vm_pu"], net.res_bus[f"{s}_vm_pu"].values, atol=1e-8, rtol=0), f'for bus voltages for {s}, max error {np.abs(res["bus"][f"{s}_vm_pu"] - net.res_bus[f"{s}_vm_pu"].values).max():.2e}'
            assert np.allclose(np.nan_to_num(res["line"][f"{s}_loading_percent"]),
                            net.res_line[f"{s}_loading_percent"].values, atol=1e-6, rtol=0), s
            if len(net.trafo) > 0:
                assert np.allclose(np.nan_to_num(res["trafo"][f"{s}_loading_percent"]),
                                net.res_trafo[f"{s}_loading_percent"].values, atol=1e-6, rtol=0), s
            
    def test_case118(self):
        net = pp.networks.case118()
        cont_lines = net.line.index.values
        cont_trafo = net.trafo.index.values
        
        # not the same behaviour in pp and ls for contingency of these lines / trafo
        lines_diff_behaviour = [6, 7, 103, 121, 163, 164, 170]
        trafos_diff_behaviour = [11, 12]
        #################
        # check that indeed the behaviour is different (ie that at least one bus is disconnected in pandapower)
        for el in lines_diff_behaviour:
            test_net = copy.deepcopy(net)
            test_net.line.loc[el, "in_service"] = False
            pp.runpp(test_net)
            assert (~np.isfinite(test_net.res_bus["vm_pu"])).any()
        for el in trafos_diff_behaviour:
            test_net = copy.deepcopy(net)
            test_net.trafo.loc[el, "in_service"] = False
            pp.runpp(test_net)
            assert (~np.isfinite(test_net.res_bus["vm_pu"])).any()
        # end of the "test of the integration test"
        #############
        
        # now the contingency analysis
        cont_lines = np.delete(cont_lines, lines_diff_behaviour)
        cont_trafo = np.delete(cont_trafo, trafos_diff_behaviour)
        nminus1_cases = {"line": {"index": cont_lines},
                         "trafo": {"index": cont_trafo
                                  }}
        self._aux_test_case(net, nminus1_cases)

[BDonnot]

Ok, now (see messages above) that the test pass with "legacy" pandapower, I installed pandapower from github, from the "develop" branch.

I then loaded the ieee118 and to my surprise, there are not "impedance" elements there...

import pandapower
print(pandapower.__version__)
case = pandapower.networks.case118()
print(case.impedance)

Give me the result:

3.0.0
Empty DataFrame
Columns: [name, from_bus, to_bus, rft_pu, xft_pu, rtf_pu, xtf_pu, sn_mva, in_service, gf_pu, gt_pu, bf_pu, bt_pu]
Index: []

By investigating however, I noticed that there has been some impacting changes in the function to compute the transfomer "model parameters" from the transformer input data (net.trafo) .

This lead to a difference in the Ybus admittance matrix and as such a difference in the resulting powerflows.

And this is why lightsim2grid gives wrong results, because the _wye_delta function does not give the same output as before for regular transformers. Is this expected ?