BUG: Non-physical power flow solution returned when simulating MV/LV network

[YassineAbdelouadoud]

Describe the bug
When simulating a MV/LV network (MV feeder to which is connected a MV/LV transformer, a LV line and a constant power load), it seems that a non-physical power flow solution can be returned for specific combination of inputs.

To Reproduce
Minimal Working Example to understand the problem

import pandas as pd
from roseau.load_flow import (ElectricalNetwork, Line, Bus, Transformer, LineParameters, LineType, ConductorType,
                              InsulatorType, Q_, PowerLoad, CurrentLoad, TransformerParameters)

results = []

line_parameters = LineParameters.from_geometry(
    "T_AL_70",
    line_type=LineType.TWISTED,
    conductor_type=ConductorType.AL,
    insulator_type=InsulatorType.PEX,
    section=70,  # mm²
    section_neutral=54,  # mm²
    height=10,  # m
    external_diameter=Q_(3.9, "cm"),
    max_current=150.
)

transformer_parameters = TransformerParameters.from_catalogue(id='FT_Standard_Standard_250kVA')
transformer_parameters.max_power = transformer_parameters.sn * 1.1

for lv_load in tqdm(np.arange(55., 70., 1.)):

    en = ElectricalNetwork.from_catalogue(name="MVFeeder210", load_point_name="Winter")
    mv_lv_bus_id = 'MVLV02676'

    lv_side_bus = Bus(id=mv_lv_bus_id + '_LV',
                      phases="abcn",
                      min_voltage=230. * 0.915,
                      max_voltage=230. * 1.085)
    # create transformer
    Transformer(id=mv_lv_bus_id + '_Transformer',
                bus1=en.buses[mv_lv_bus_id],
                bus2=lv_side_bus,
                parameters=transformer_parameters,
                phases1='abc',
                phases2='abcn')

    lv_bus_0 = Bus(id=mv_lv_bus_id + '_LV_0',
                   phases="abcn",
                   min_voltage=230. * 0.915,
                   max_voltage=230. * 1.085)

    new_line = Line(id=mv_lv_bus_id + '_line',
                    bus1=lv_side_bus,
                    bus2=lv_bus_0,
                    parameters=line_parameters,
                    phases='abcn',
                    ground=en.grounds['ground'],
                    length=0.1)
    load = PowerLoad(id=mv_lv_bus_id + '_load', bus=lv_bus_0,
                       powers=[complex(lv_load * 1.e3, 0.3 * lv_load), 0., 0.])
  
    en.solve_load_flow(auth=auth)
    v_mag = en.res_buses_voltages['voltage'].abs()
    lv_v_mag = v_mag[v_mag < 1000.]

    results.append(
        {
            'load': lv_load,
            'transfo_load1': en.res_transformers['power1'].abs().sum() * 1.e-3,
            'min_voltage': lv_v_mag.min(),
            'max_lv_voltage': lv_v_mag.max()
        }
    )

res_df = pd.DataFrame.from_records(results)
res_df['trafo_power_ratio'] = res_df['transfo_load1'] / res_df['load']
print(res_df)
    load  transfo_load1  min_voltage  max_lv_voltage  trafo_power_ratio
0   55.0      79.510371   200.891534      236.545121           1.445643
1   56.0      81.041662   200.293016      236.703442           1.447173
2   57.0      82.584240   199.690045      236.862680           1.448846
3   58.0      84.138325   199.082522      237.022856           1.450661
4   59.0     492.998626    32.109242      303.108455           8.355909
5   60.0     491.491145    32.755285      302.807896           8.191519
6   61.0      88.871999   197.231602      237.509223           1.456918
7   62.0     488.440114    34.062084      302.200464           7.878066
8   63.0      92.089861   195.972916      237.838563           1.461744
9   64.0      93.718248   195.335794      238.004836           1.464348
10  65.0      95.360005   194.693316      238.172212           1.467077
11  66.0      97.015451   194.045349      238.340720           1.469931
12  67.0      98.684915   193.391752      238.510388           1.472909
13  68.0     100.368746   192.732379      238.681246           1.476011
14  69.0     102.067307   192.067078      238.853326           1.479236

The same simulation, but using CurrentLoad instead of PowerLoad yields coherent results :

load = CurrentLoad(id=mv_lv_bus_id + '_LV_0_load', bus=lv_bus_0,
                           currents=[complex(lv_load * 1.e3, 0.3 * lv_load) / 230., 0., 0.] )

    load  transfo_load1  min_voltage  max_lv_voltage  trafo_power_ratio
0   55.0      74.132574   203.418299      239.287685           1.347865
1   56.0      75.272632   202.956715      239.479149           1.344154
2   57.0      76.412936   202.495167      239.670675           1.340578
3   58.0      77.553474   202.033657      239.862261           1.337129
4   59.0      78.694233   201.572184      240.053908           1.333801
5   60.0      79.835203   201.110749      240.245616           1.330587
6   61.0      80.976373   200.649352      240.437384           1.327482
7   62.0      82.117732   200.187994      240.629213           1.324480
8   63.0      83.259271   199.726673      240.821102           1.321576
9   64.0      84.400982   199.265392      241.013051           1.318765
10  65.0      85.542856   198.804149      241.205061           1.316044
11  66.0      86.684885   198.342946      241.397130           1.313407
12  67.0      87.827062   197.881783      241.589258           1.310852
13  68.0      88.969380   197.420659      241.781447           1.308373
14  69.0      90.111832   196.959575      241.973694           1.305969

The same simulation as the first one, but replacing the MV feeder by a voltage source also yields coherent results :

    load  transfo_load1  min_voltage  max_lv_voltage  trafo_power_ratio
0   55.0      79.415539   203.461498      238.598885           1.443919
1   56.0      80.934474   202.876465      238.756005           1.445258
2   57.0      82.464217   202.287252      238.914001           1.446741
3   58.0      84.004967   201.693773      239.072891           1.448362
4   59.0      85.556934   201.095934      239.232694           1.450118
5   60.0      87.120337   200.493641      239.393431           1.452006
6   61.0      88.695406   199.886795      239.555123           1.454023
7   62.0      90.282379   199.275294      239.717790           1.456167
8   63.0      91.881507   198.659032      239.881455           1.458437
9   64.0      93.493051   198.037898      240.046142           1.460829
10  65.0      95.117283   197.411778      240.211874           1.463343
11  66.0      96.754490   196.780553      240.378677           1.465977
12  67.0      98.404968   196.144098      240.546577           1.468731
13  68.0     100.069028   195.502286      240.715600           1.471603
14  69.0     101.746997   194.854980      240.885776           1.474594

Expected behavior
In the first simulation, we can observe that the lowest voltage magnitude drops from 199 V to 32 V when increasing the load from 58 to 59 kW, while the voltage drop from 57 kW to 58 kW was only 0.6 V. The same happens for 60 kW and 62 kW. This voltage drop causes huge current (and thus losses) in the line, which translates into overloading of the transformer.

I expect the lowest voltage to be between 197 and 199 V for the 59 and 60 kW.

Version details
roseau-load-flow 0.6.0

System Information
------------------
python               3.10.8.final.0
os                   linux
os_name              posix
machine              x86_64
Installed Dependencies
----------------------
pandas               2.1.0
numpy                1.23.5
geopandas            0.13.2
shapely              2.0.1
regex                2022.10.31
pint                 0.22
requests             2.31.0

[Saelyos]

I've just tried with current develop branch, and the issue seems resolved : I get the same results, except for the lines 4, 5, 7, that are now coherents :

    load  transfo_load1  min_voltage  max_lv_voltage  trafo_power_ratio
0   55.0      79.510371   200.891534      236.545121           1.445643
1   56.0      81.041662   200.293016      236.703442           1.447173
2   57.0      82.584240   199.690045      236.862680           1.448846
3   58.0      84.138325   199.082522      237.022856           1.450661
4   59.0      85.704150   198.470346      237.183990           1.452613
5   60.0      87.281957   197.853409      237.346105           1.454699
6   61.0      88.871999   197.231602      237.509223           1.456918
7   62.0      90.474541   196.604810      237.673367           1.459267
8   63.0      92.089861   195.972916      237.838563           1.461744
9   64.0      93.718248   195.335794      238.004836           1.464348
10  65.0      95.360005   194.693316      238.172212           1.467077
11  66.0      97.015451   194.045349      238.340720           1.469931
12  67.0      98.684915   193.391752      238.510388           1.472909
13  68.0     100.368746   192.732379      238.681246           1.476011
14  69.0     102.067307   192.067078      238.853326           1.479236

[Saelyos]

I could reproduce the bug, it was coming from an incorrect initial potentials propagation. The good results are now always returned, with sometimes even less iterations, so thanks for the report !

[YassineAbdelouadoud]

Thanks @Saelyos, this is very good news as this issue was affecting most of the MV/LV simulations I was trying to make. Will this fix be available in the next release ?

[benoit9126]

Yes! :-)