Repository: Grid2op/lightsim2grid (language: python)

Issue:
Fixed Decoupled Power Flow algorithm in LightSim2Grid
I was wondering if you looked at implementing the Fixed Decoupled Power Flow algorithm in LightSim2Grid? I think this method is faster as the factorization can be reused directly between iterations. The method is implemented in MATPOWER. The relevant parts of the code can be found at the following links:
1. [https://github.com/MATPOWER/matpower/blob/3232f1ec3df214b162c9073227aa23ce933fcb73/lib/runpf.m#L336](https://github.com/MATPOWER/matpower/blob/3232f1ec3df214b162c9073227aa23ce933fcb73/lib/runpf.m#L336)
2. [https://github.com/MATPOWER/matpower/blob/3232f1ec3df214b162c9073227aa23ce933fcb73/lib/runpf.m#L337](https://github.com/MATPOWER/matpower/blob/3232f1ec3df214b162c9073227aa23ce933fcb73/lib/runpf.m#L337)
3. [https://github.com/MATPOWER/matpower/blob/master/lib/fdpf.m](https://github.com/MATPOWER/matpower/blob/master/lib/fdpf.m)

The B matrices needed for this algorithm can be easily created from the Ybus matrix . The MATPOWER code for creating B matrices can be found in the following link. The two key lines are line 75 and 84. These matrices are formed from submatrices of the imaginary part of Ybus. 

[https://github.com/MATPOWER/matpower/blob/3232f1ec3df214b162c9073227aa23ce933fcb73/lib/makeB.m](https://github.com/MATPOWER/matpower/blob/3232f1ec3df214b162c9073227aa23ce933fcb73/lib/makeB.m)

Feature Request : Multi-Slack  (Non Distributed) 
Hello Mr Bonnot,

We are using PandaPower to perform our power flow calculations, we have discovered LightSim2grid via PandaPower, we are interested in using it because of it's good shown performance.

We have several sub-stations which means several slacks non distributed, a use case example is mv_oberrhein in PandaPower Network, we have remarked that LightSim2grid is not activated for this use case.  We would like to use it in our works in order to help us to improve the speed performance.

Is there any way to use LighSim2Grid on a multi-slack (Non Distributed) network ?, if yes, is it easy to handle ?, are you planning to add this feature in the future ?

Thank you so much in advance.

Regards




Discussion / hints:

Hello,

Thanks for your interest in this package.

Currently only single or distributed slack i'are supported. I will have a lookin how to add multi slacks. And depending on how complicated it is I might add the feature in the future.

Note however that this will take some time (best case scenario I can look end of September) as I cannot work on lightsim2grid at the moment.

Benjamin 
Hello Bonnot,

Hope you are doing well, did you have some time to check the feasibility of introducing Multi-Slack on Lightsim2grid ? can you share with us your update and plans? thank you in advance.

Regards
Hello,

I am on parental leave until beginning of December. And then I have a 2 weeks holiday for Christmas.

I cannot work on this issue before early 2023 at best. Any pull request will drastically reduce the time it takes for me to do it.

I think it should be possible to implement it rather easily as I'm pretty sure it does not require heavy development.

If I may ask, out of curiosity, why would you like multiple slack in a non distributed way? This model does not really make sense for me but I'm not a power system expert (but I'm really curious 😉) 

Best

Benjamin 
Hello,

Congratulation for the new baby 🎉, and thanks for your interest. In fact, it's for the power flow calculation of several sub-stations with each one having its own Slack, it will be very quick to compute PF for all sub-stations at once instead of each one separately, Pandapower actually handles this parallel computation of PF (power flow) over multiple sub-stations, this is much faster compared to the sequential computation of PF for those multiple sub-stations.

you can imagine the task as parallelization of newton-raphson computation over multiple grids.

Thanks 
Thanks 😊


Oh OK I see. In this case it makes perfect sense indeed. I'll see what I can do asap (but again probably not in 2022 :-/)

Best

Benjamin 
Hello,
I think i just need to change one single function in the c++ side so it should be easy. By any chance, do you have any file you are using (for the grid) and script (for the code) so that i can make experiments with it asap ?

If not i will actually spend more time to find a suitable grid and then create a code with pandapower to do what you want to do rather than to code the actual feature in lightsim2grid.

Thanks
Hello Benjamin,
Thank you for getting back to us. Since we are not allowed to upload our networks, the easiest would be for you to use the MV network already available in pandapower :

import pandapower.networks
import pandapower as pp
net = pandapower.networks.mv_oberrhein()
pp.runpp(net, lightsim2grid=False) # works
pp.runpp(net, lightsim2grid=True) # NotImplementedError: option 'lightsim2grid' is True and multiple ext_grids are found, but distributed_slack=False.

The runtime with this network won’t be representative since the network is really small but it might still be useful for you to experiment with. Don’t hesitate to tell us if it is not
Hello,

As a first step this sounds great. I'll get back to you once i've done it (in a separate branch, probably next week) with the binary you can install and you can let me know if that works for you. If it does i'll merge it to the right repo.

Does that sound ok for you ?
Hello,

That's great! thank you! 
yes, it's ok for us.


Hello,

I have not forget this feature. I was swamped by work and there were lots of different bugs in grid2op and lightsim2grid... 

I'll work on it as soon as I can
A quick update on the matter: things are a bit more complicated, especially if i need to initialize the grid with DC.

I will try to see if I get the same results as pandapower for `mv_oberrhein()` grid when I let pandapower compute the initial vector.
Hello,

Thanks for keeping me updated and for your involvement in the matter. 
Hello,

Do you have any updates to share, is there something I can do to help? can you give some details about the complexity and what needs to be modified?

Thank you in advance

Regards 
Hello,

I tried to run it with the "simple ideas" that I got, but unfortunately it did not work. I think I need to double check in the case of multiple (non distributed) slack how pandapower does it. Because clearly, removing some buses that would be neither pv nor pq (so the slack) is not enough.
There is probably another trick that I don't really "master" right now :-/

Fix patch:
diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml
index f361d22..cc48a5a 100644
--- a/.github/workflows/main.yml
+++ b/.github/workflows/main.yml
@@ -41,6 +41,11 @@ jobs:
             abi: cp311,
             version: '3.11',
           }
+          - {
+            name: cp312,
+            abi: cp312,
+            version: '3.12',
+          }
 
     steps:
 
@@ -73,21 +78,33 @@ jobs:
       - name: Install wheel
         run: pip3 install wheelhouse/*.whl --user
 
-      - name: Install GDB
-        run: yum install -y gdb
+      # - name: Install GDB
+      #   run: yum install -y gdb
 
       - name: Check package can be imported
         run: |
-          # python3 -c "import lightsim2grid"
-          # python3 -c "from lightsim2grid import *"
-          # python3 -c "from lightsim2grid.newtonpf import newtonpf"
-          gdb -batch -ex r --args python3 -c "from lightsim2grid import *" -ex backtrace
-
-      - name: Check LightSimBackend can be imported
-        run: |
+          python3 -c "import lightsim2grid"
+          python3 -c "from lightsim2grid import *"
+          python3 -c "from lightsim2grid.newtonpf import newtonpf"
+      
+      - name: Fix urllib3 (python 3.7)
+        if: matrix.python.name == 'cp37'
+        run:
+          pip install urllib3==1.26.6
+          # otherwise urllib3 v2.0 only supports OpenSSL 1.1.1+, currently the 'ssl' module is compiled with 'OpenSSL 1.0.2k-fips  26 Jan 2017'
+
+      - name: Check LightSimBackend can be imported (install grid2op)
+        run: | 
           python3 -m pip install grid2op
-          python3 -c "from lightsim2grid import LightSimBackend"
-          python3 -c "from lightsim2grid import LightSimBackend; import grid2op; env = grid2op.make('l2rpn_case14_sandbox', test=True, backend=LightSimBackend())"
+          python3 -m pip freeze
+
+      - name: Check LightSimBackend can be imported 1
+        run:
+          python3 -v -c "from lightsim2grid import LightSimBackend"
+
+      - name: Check LightSimBackend can be imported 2
+        run:
+          python3 -v -c "from lightsim2grid import LightSimBackend; import grid2op; env = grid2op.make('l2rpn_case14_sandbox', test=True, backend=LightSimBackend())"
 
       - name: Upload wheel
         uses: actions/upload-artifact@v3
@@ -129,6 +146,10 @@ jobs:
             name: cp311,
             version: '3.11',
           }
+          - {
+            name: cp312,
+            version: '3.12',
+          }
     env:
       RUNNER_OS: windows-2019
       PYTHON_VERSION: ${{ matrix.python.version }}
@@ -148,6 +169,7 @@ jobs:
       - name: Install Python dependencies
         run: |
           python -m pip install --upgrade pip
+          python -m pip install --upgrade setuptools
           python -m pip install -r requirements.txt
 
       - name: Compile SuiteSparse cmake
@@ -270,6 +292,10 @@ jobs:
             name: cp311,
             version: '3.11',
           }
+          - {
+            name: cp312,
+            version: '3.12',
+          }
     env:
       RUNNER_OS: macos-latest
       PYTHON_VERSION: ${{ matrix.python.version }}
@@ -293,15 +319,16 @@ jobs:
           echo "DEPLOY=$( [[ $GITHUB_EVENT_NAME == 'push' && $GITHUB_REF == 'refs/tags'* ]] && echo 'True' || echo 'False' )" >> $GITHUB_ENV
 
       - name: Compile with cibuildwheel
-        uses: pypa/cibuildwheel@v2.15.0
+        uses: pypa/cibuildwheel@v2.16.2
         env:
           CIBW_BUILD: "cp3${{env.PYTHON_SUBVERSION}}-*"  # see https://github.com/cvxpy/cvxpy/blob/master/.github/workflows/build.yml
           CIBW_ARCHS_LINUX: auto aarch64  # not used I believe
           CIBW_ARCHS_MACOS: x86_64 arm64
           CIBW_ENVIRONMENT: __O3_OPTIM=1
-          CIBW_BEFORE_BUILD: make clean && make  # and not CIBW_BEFORE_ALL !
+          CIBW_BEFORE_ALL: python -m pip install --upgrade setuptools wheel pybind11
+          CIBW_BEFORE_BUILD: make clean && make CC=clang CXX=clang++  # and not CIBW_BEFORE_ALL !
           CIBW_TEST_REQUIRES: grid2op pandapower
-          CIBW_TEST_SKIP: "*-macosx_arm64"  # to silence warning "While arm64 wheels can be built on x86_64, they cannot be tested."
+          CIBW_TEST_SKIP: "cp312-* *-macosx_arm64"  # to silence warning "While arm64 wheels can be built on x86_64, they cannot be tested."
           CIBW_TEST_COMMAND: >
             python -c "import lightsim2grid" && 
             python -c "from lightsim2grid import *" &&
diff --git a/.gitignore b/.gitignore
index 9b18bfa..fac4b54 100644
--- a/.gitignore
+++ b/.gitignore
@@ -246,8 +246,9 @@ benchmarks/benchmark.png
 benchmarks/benchmark2.png
 cktso
 _pandapower
-benchmarks/slacks.py
+_pypower
 _gridcal
+benchmarks/slacks.py
 *grid.json
 bug_pandapower.py
 bug_pp.py
@@ -263,3 +264,8 @@ nicslu
 **nicslu.lic
 test_milad.py
 test_issue_50.py
+reverse_pypow.py
+test_issue_58.py
+test_issue_58_old.py
+lightsim2grid/gridmodel/pypowsybl.ipynb
+lightsim2grid/tests/test.mat
diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index bda46b0..147c4ab 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -3,13 +3,14 @@ Change Log
 
 [TODO]
 --------
-- make an `init` function from pypowsybl
-- support DC powerline (as modeled in pandapower)
+- [refacto] have a structure in cpp for the buses
+- [refacto] have the id_grid_to_solver and id_solver_to_grid etc. directly in the solver and NOT in the gridmodel.
 - support 3w trafo (as modeled in pandapower)
 - improve speed by not performing internal checks 
   (keep check for boundaries and all for python API instead) [see `TODO DEBUG MODE` in c++ code]
 - improve speed
 - code parrallelism directly in the `Computer` and `SecurityAnalysisCPP` classes
+- a mode to do both `Computer` and `SecurityAnalysisCPP`
 - use the "multi slack hack" (see issue #50) for SecurityAnalysis or Computer for example
 - code `helm` powerflow method
 - possibility to read CGMES files
@@ -17,9 +18,27 @@ Change Log
 - interface with gridpack (to enforce q limits for example)
 - maybe have a look at suitesparse "sliplu" tools ?
 - easier building (get rid of the "make" part)
-- code NR with dense matrices
 
-[0.7.3] 2023-08-24
+[0.7.5] 2023-10-05
+--------------------
+- [FIXED] a bug in DC powerflow when asking for computation time: it was not reset to 0. when
+  multiple powerflows used the same solver
+- [FIXED] a bug in AC and DC powerflow when shunts had active values
+- [ADDED] possibility to initialize a powergrid based on pypowsybl 
+  see https://github.com/BDonnot/lightsim2grid/issues/53
+- [ADDED] some more algorithm to perform powerflow: Fast Decoupled Powerflow (in BX and XB variant)
+  see https://github.com/BDonnot/lightsim2grid/issues/63
+- [ADDED] build lightsim2grid for python 3.12
+- [ADDED] support for non distributed slack but multiple slack buses
+  see https://github.com/BDonnot/lightsim2grid/issues/50 (ONLY FOR AC powerflow)
+- [IMPROVED] now shipping `src` and `eigen` directory in the source of 
+  lightsim2grid to allow their installation if wheels are not provided.
+- [IMPROVED] in the underlying cpp GridModel powerlines can now have 2
+  different values for the `h` parameters (`h_or` and `h_ex`).
+- [IMPROVED] now lightsim2grid is able to load a pandapower network with non
+  contiguous non starting at 0 bus index
+
+[0.7.3/4] 2023-08-24
 --------------------
 - [FIXED] a bug where, when you disconnect a load (or gen), the next action cannot be performed
   if it modifies the load (or gen), because you "cannot change the value of a disconnected load (or gen)"
diff --git a/MANIFEST.in b/MANIFEST.in
new file mode 100644
index 0000000..e99357e
--- /dev/null
+++ b/MANIFEST.in
@@ -0,0 +1,18 @@
+graft eigen
+graft src
+prune eigen/doc
+prune eigen/unsupported
+prune eigen/test
+prune eigen/failtest
+prune eigen/demos
+prune eigen/scripts
+prune eigen/cmake
+prune eigen/ci
+prune eigen/bench
+prune eigen/.gitlab
+prune eigen/.git
+prune eigen/.gitignore
+prune eigen/debug
+recursive-exclude eigen *.git*
+global-exclude __pycache__
+global-exclude *.py[co]
\ No newline at end of file
diff --git a/Makefile b/Makefile
index 5ef3f93..563ea41 100644
--- a/Makefile
+++ b/Makefile
@@ -12,8 +12,10 @@
 ###################################
 default: all
 
-CC = gcc
-CXX = g++
+# CC = gcc
+# CXX = g++
+CC := $(if $(CC),$(CC),gcc)
+CXX := $(if $(CXX),$(CXX),g++)
 
 LIBPATH_REL = ./SuiteSparse
 LIBPATH = $(realpath $(LIBPATH_REL))
@@ -54,12 +56,12 @@ distclean: purge
 all: prelude
 
 prelude:
-	(cd $(LIBPATH)/SuiteSparse_config/ && CC=gcc make)
-	(cd $(LIBPATH)/CXSparse/Lib/ && CC=gcc make )
-	(cd $(LIBPATH)/AMD/Lib/ && Cc=gcc make)
-	(cd $(LIBPATH)/BTF/Lib/ && CC=gcc make)
-	(cd $(LIBPATH)/COLAMD/Lib/ && CC=gcc make)
-	(cd $(LIBPATH)/KLU/Lib/ && CC=gcc make)
+	(cd $(LIBPATH)/SuiteSparse_config/ && CC=$(CC) make)
+	(cd $(LIBPATH)/CXSparse/Lib/ && CC=$(CC) make )
+	(cd $(LIBPATH)/AMD/Lib/ && CC=$(CC) make)
+	(cd $(LIBPATH)/BTF/Lib/ && CC=$(CC) make)
+	(cd $(LIBPATH)/COLAMD/Lib/ && CC=$(CC) make)
+	(cd $(LIBPATH)/KLU/Lib/ && CC=$(CC) make)
 
 ##################################
 # this is the documentation
diff --git a/docs/conf.py b/docs/conf.py
index 8b64a4f..8438a0a 100644
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -22,7 +22,11 @@
 author = 'Benjamin DONNOT'
 
 # The full version, including alpha/beta/rc tags
+<<<<<<< HEAD
+release = "0.7.5.dev0"
+=======
 release = "0.7.4"
+>>>>>>> master
 version = '0.7'
 
 # -- General configuration ---------------------------------------------------
diff --git a/lightsim2grid/__init__.py b/lightsim2grid/__init__.py
index f7d6d89..145b6f3 100644
--- a/lightsim2grid/__init__.py
+++ b/lightsim2grid/__init__.py
@@ -5,11 +5,7 @@
 # you can obtain one at http://mozilla.org/MPL/2.0/.
 # SPDX-License-Identifier: MPL-2.0
 # This file is part of LightSim2grid, LightSim2grid implements a c++ backend targeting the Grid2Op platform.
-
-import faulthandler
-faulthandler.enable()
-
-__version__ = "0.7.4"
+__version__ = "0.7.5"
 
 __all__ = ["newtonpf", "SolverType", "ErrorType", "solver"]
 
@@ -20,9 +16,10 @@
 try:
     from lightsim2grid.lightSimBackend import LightSimBackend
     __all__.append("LightSimBackend")
-except ImportError:
+except ImportError as exc_:
     # grid2op is not installed, the Backend will not be available
     pass
+    print(f"LightSimBackend import error: {exc_}")
 
 try:
     from lightsim2grid.physical_law_checker import PhysicalLawChecker
@@ -30,6 +27,7 @@
 except ImportError as exc_:
     # grid2op is not installed, the PhysicalLawChecker will not be available
     pass
+    print(f"PhysicalLawChecker import error: {exc_}")
 
 try:
     from lightsim2grid.timeSerie import TimeSerie
@@ -38,6 +36,7 @@
 except ImportError as exc_:
     # grid2op is not installed, the TimeSeries module will not be available
     pass
+    print(f"TimeSerie import error: {exc_}")
 
 try:
     from lightsim2grid.securityAnalysis import SecurityAnalysis
@@ -46,3 +45,4 @@
 except ImportError as exc_:
     # grid2op is not installed, the SecurtiyAnalysis module will not be available
     pass
+    print(f"SecurityAnalysis import error: {exc_}")
diff --git a/lightsim2grid/gridmodel/_aux_add_dc_line.py b/lightsim2grid/gridmodel/_aux_add_dc_line.py
index 41a9bee..d3d8d8b 100644
--- a/lightsim2grid/gridmodel/_aux_add_dc_line.py
+++ b/lightsim2grid/gridmodel/_aux_add_dc_line.py
@@ -5,13 +5,15 @@
 # you can obtain one at http://mozilla.org/MPL/2.0/.
 # SPDX-License-Identifier: MPL-2.0
 # This file is part of LightSim2grid, LightSim2grid implements a c++ backend targeting the Grid2Op platform.
+
 import warnings
 
 import numpy as np
 from ._aux_add_sgen import SOME_KIND_OF_INF_FOR_PMIN_PMAX
+from ._pp_bus_to_ls_bus import pp_bus_to_ls
 
 
-def _aux_add_dc_line(model, pp_net):
+def _aux_add_dc_line(model, pp_net, pp_to_ls):
     """
     Add the transformers of the pp_net into the lightsim2grid "model"
 
@@ -29,9 +31,12 @@ def _aux_add_dc_line(model, pp_net):
         raise RuntimeError("Cannot handle 'parallel' dcline columns. Please duplicate the rows if that is the case. "
                            "Some pp_net.dcline[\"parallel\"] != 1 it is not handled by lightsim yet.")
 
-
-    branch_from_id = pp_net.dcline["from_bus"].values
-    branch_to_id = pp_net.dcline["to_bus"].values
+    if pp_net.dcline.shape[0] == 0:
+        # nothing to do if no dc line
+        return
+    
+    branch_from_id = pp_bus_to_ls(pp_net.dcline["from_bus"].values, pp_to_ls)
+    branch_to_id = pp_bus_to_ls(pp_net.dcline["to_bus"].values, pp_to_ls)
     p_mw = -pp_net.dcline["p_mw"].values
     if np.any(~np.isfinite(p_mw)):
         warnings.warn("Some non finite values are found for p_mw, they have been replaced by 0.")
diff --git a/lightsim2grid/gridmodel/_aux_add_gen.py b/lightsim2grid/gridmodel/_aux_add_gen.py
index 2590b4a..c864486 100644
--- a/lightsim2grid/gridmodel/_aux_add_gen.py
+++ b/lightsim2grid/gridmodel/_aux_add_gen.py
@@ -1,14 +1,16 @@
-# Copyright (c) 2020, RTE (https://www.rte-france.com)
+# Copyright (c) 2020-2023, RTE (https://www.rte-france.com)
 # See AUTHORS.txt
 # This Source Code Form is subject to the terms of the Mozilla Public License, version 2.0.
 # If a copy of the Mozilla Public License, version 2.0 was not distributed with this file,
 # you can obtain one at http://mozilla.org/MPL/2.0/.
 # SPDX-License-Identifier: MPL-2.0
 # This file is part of LightSim2grid, LightSim2grid implements a c++ backend targeting the Grid2Op platform.
+
 import numpy as np
+from ._pp_bus_to_ls_bus import pp_bus_to_ls
 
 
-def _aux_add_gen(model, pp_net):
+def _aux_add_gen(model, pp_net, pp_to_ls):
     """
     Add the generators of the pp_net into the lightsim2grid "model"
 
@@ -26,7 +28,7 @@ def _aux_add_gen(model, pp_net):
                           pp_net.gen["vm_pu"].values,
                           pp_net.gen["min_q_mvar"].values,
                           pp_net.gen["max_q_mvar"].values,
-                          pp_net.gen["bus"].values
+                          pp_bus_to_ls(pp_net.gen["bus"].values, pp_to_ls)
                           )
     for gen_id, is_connected in enumerate(pp_net.gen["in_service"].values):
         if not is_connected:
diff --git a/lightsim2grid/gridmodel/_aux_add_line.py b/lightsim2grid/gridmodel/_aux_add_line.py
index e63fbac..98accdf 100644
--- a/lightsim2grid/gridmodel/_aux_add_line.py
+++ b/lightsim2grid/gridmodel/_aux_add_line.py
@@ -1,15 +1,16 @@
-# Copyright (c) 2020, RTE (https://www.rte-france.com)
+# Copyright (c) 2020-2023, RTE (https://www.rte-france.com)
 # See AUTHORS.txt
 # This Source Code Form is subject to the terms of the Mozilla Public License, version 2.0.
 # If a copy of the Mozilla Public License, version 2.0 was not distributed with this file,
 # you can obtain one at http://mozilla.org/MPL/2.0/.
 # SPDX-License-Identifier: MPL-2.0
 # This file is part of LightSim2grid, LightSim2grid implements a c++ backend targeting the Grid2Op platform.
-import warnings
+
 import numpy as np
 
+from ._pp_bus_to_ls_bus import pp_bus_to_ls
 
-def _aux_add_line(converter, model, pp_net):
+def _aux_add_line(converter, model, pp_net, pp_to_ls=None):
     """
     Add the lines of the pp_net into the lightsim2grid "model"
 
@@ -34,14 +35,14 @@ def _aux_add_line(converter, model, pp_net):
             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"]
-        )
+            pp_net.bus.loc[pp_net.line["from_bus"]]["vn_kv"], 
+            pp_net.bus.loc[pp_net.line["to_bus"]]["vn_kv"], 
+            )
 
     ### add them to the grid
     model.init_powerlines(line_r, line_x, line_h,
-                          pp_net.line["from_bus"].values,
-                          pp_net.line["to_bus"].values
+                          pp_bus_to_ls(pp_net.line["from_bus"].values, pp_to_ls),
+                          pp_bus_to_ls(pp_net.line["to_bus"].values, pp_to_ls)
                           )
     for line_id, is_connected in enumerate(pp_net.line["in_service"].values):
         if not is_connected:
diff --git a/lightsim2grid/gridmodel/_aux_add_load.py b/lightsim2grid/gridmodel/_aux_add_load.py
index c30e3cf..14b3695 100644
--- a/lightsim2grid/gridmodel/_aux_add_load.py
+++ b/lightsim2grid/gridmodel/_aux_add_load.py
@@ -1,14 +1,16 @@
-# Copyright (c) 2020, RTE (https://www.rte-france.com)
+# Copyright (c) 2020-2023, RTE (https://www.rte-france.com)
 # See AUTHORS.txt
 # This Source Code Form is subject to the terms of the Mozilla Public License, version 2.0.
 # If a copy of the Mozilla Public License, version 2.0 was not distributed with this file,
 # you can obtain one at http://mozilla.org/MPL/2.0/.
 # SPDX-License-Identifier: MPL-2.0
 # This file is part of LightSim2grid, LightSim2grid implements a c++ backend targeting the Grid2Op platform.
+
 import numpy as np
 
+from ._pp_bus_to_ls_bus import pp_bus_to_ls
 
-def _aux_add_load(model, pp_net):
+def _aux_add_load(model, pp_net, pp_to_ls):
     """
     Add the load of the pp_net into the lightsim2grid "model"
 
@@ -25,9 +27,14 @@ def _aux_add_load(model, pp_net):
         raise RuntimeError("Cannot handle 'parallel' load columns. Please duplicate the rows if that is the case. "
                            "Some pp_net.load[\"parallel\"] != 1 it is not handled by lightsim yet.")
 
-    model.init_loads(pp_net.load["p_mw"].values,
-                     pp_net.load["q_mvar"].values,
-                     pp_net.load["bus"].values
+    ratio = 1.0
+    if "scaling" in pp_net.load:
+        ratio = pp_net.load["scaling"].values
+        ratio[~np.isfinite(ratio)] = 1.0
+        
+    model.init_loads(pp_net.load["p_mw"].values * ratio,
+                     pp_net.load["q_mvar"].values * ratio,
+                     pp_bus_to_ls(pp_net.load["bus"].values, pp_to_ls)
                      )
     for load_id, is_connected in enumerate(pp_net.load["in_service"].values):
         if not is_connected:
diff --git a/lightsim2grid/gridmodel/_aux_add_sgen.py b/lightsim2grid/gridmodel/_aux_add_sgen.py
index 448b658..c454728 100644
--- a/lightsim2grid/gridmodel/_aux_add_sgen.py
+++ b/lightsim2grid/gridmodel/_aux_add_sgen.py
@@ -1,17 +1,19 @@
-# Copyright (c) 2020, RTE (https://www.rte-france.com)
+# Copyright (c) 2020-2023, RTE (https://www.rte-france.com)
 # See AUTHORS.txt
 # This Source Code Form is subject to the terms of the Mozilla Public License, version 2.0.
 # If a copy of the Mozilla Public License, version 2.0 was not distributed with this file,
 # you can obtain one at http://mozilla.org/MPL/2.0/.
 # SPDX-License-Identifier: MPL-2.0
 # This file is part of LightSim2grid, LightSim2grid implements a c++ backend targeting the Grid2Op platform.
+
 import warnings
 import numpy as np
+from ._pp_bus_to_ls_bus import pp_bus_to_ls
 
 SOME_KIND_OF_INF_FOR_PMIN_PMAX = 99999.
 
 
-def _aux_add_sgen(model, pp_net):
+def _aux_add_sgen(model, pp_net, pp_to_ls):
     """
     Add the static generators  (=PQ generators) of the pp_net into the lightsim2grid "model"
 
@@ -64,13 +66,18 @@ def _aux_add_sgen(model, pp_net):
         warnings.warn("There were some Nan in the pp_net.sgen[\"max_q_mvar\"], they have been replaced by 0")
     max_q_mvar[~np.isfinite(max_q_mvar)] = 0.
 
-    model.init_sgens(pp_net.sgen["p_mw"].values,
-                     pp_net.sgen["q_mvar"].values,
+    ratio = 1.0
+    if "scaling" in pp_net.sgen:
+        ratio = pp_net.sgen["scaling"].values
+        ratio[~np.isfinite(ratio)] = 1.0
+        
+    model.init_sgens(pp_net.sgen["p_mw"].values * ratio,
+                     pp_net.sgen["q_mvar"].values * ratio,
                      min_p_mw,
                      max_p_mw,
                      min_q_mvar,
                      max_q_mvar,
-                     pp_net.sgen["bus"].values
+                     pp_bus_to_ls(pp_net.sgen["bus"].values, pp_to_ls)
                      )
     for sgen_id, is_connected in enumerate(pp_net.sgen["in_service"].values):
         if not is_connected:
diff --git a/lightsim2grid/gridmodel/_aux_add_shunt.py b/lightsim2grid/gridmodel/_aux_add_shunt.py
index 58d0db0..c1d3e52 100644
--- a/lightsim2grid/gridmodel/_aux_add_shunt.py
+++ b/lightsim2grid/gridmodel/_aux_add_shunt.py
@@ -5,10 +5,12 @@
 # you can obtain one at http://mozilla.org/MPL/2.0/.
 # SPDX-License-Identifier: MPL-2.0
 # This file is part of LightSim2grid, LightSim2grid implements a c++ backend targeting the Grid2Op platform.
+
 import numpy as np
+from ._pp_bus_to_ls_bus import pp_bus_to_ls
 
 
-def _aux_add_shunt(model, pp_net):
+def _aux_add_shunt(model, pp_net, pp_to_ls):
     """
     Add the shunts of the pp_net into the lightsim2grid "model"
 
@@ -25,9 +27,13 @@ def _aux_add_shunt(model, pp_net):
         raise RuntimeError("Cannot handle 'parallel' sgen columns. Please duplicate the rows if that is the case. "
                            "Some pp_net.sgen[\"parallel\"] != 1 it is not handled by lightsim yet.")
 
+    if pp_net.shunt.shape[0] == 0:
+        # nothing to do if no shunt
+        return
+    
     model.init_shunt(pp_net.shunt["p_mw"].values,
                      pp_net.shunt["q_mvar"].values,
-                     pp_net.shunt["bus"].values
+                     pp_bus_to_ls(pp_net.shunt["bus"].values, pp_to_ls)
                      )
     for sh_id, is_connected in enumerate(pp_net.shunt["in_service"].values):
         if not is_connected:
diff --git a/lightsim2grid/gridmodel/_aux_add_slack.py b/lightsim2grid/gridmodel/_aux_add_slack.py
index 8900c78..17b9c83 100644
--- a/lightsim2grid/gridmodel/_aux_add_slack.py
+++ b/lightsim2grid/gridmodel/_aux_add_slack.py
@@ -1,15 +1,17 @@
-# Copyright (c) 2020, RTE (https://www.rte-france.com)
+# Copyright (c) 2020-2023, RTE (https://www.rte-france.com)
 # See AUTHORS.txt
 # This Source Code Form is subject to the terms of the Mozilla Public License, version 2.0.
 # If a copy of the Mozilla Public License, version 2.0 was not distributed with this file,
 # you can obtain one at http://mozilla.org/MPL/2.0/.
 # SPDX-License-Identifier: MPL-2.0
 # This file is part of LightSim2grid, LightSim2grid implements a c++ backend targeting the Grid2Op platform.
+
 import warnings
 import numpy as np
+from ._pp_bus_to_ls_bus import pp_bus_to_ls
 
 
-def _aux_add_slack(model, pp_net):
+def _aux_add_slack(model, pp_net, pp_to_ls):
     """
     add the slack bus(es) to the lightsim2grid "model" based on the information in the pandapower network.
 
@@ -49,7 +51,7 @@ def _aux_add_slack(model, pp_net):
             
         # in this case the ext grid is not taken into account, i raise a warning if
         # there is one
-        slack_bus_ids = pp_net.ext_grid["bus"].values
+        slack_bus_ids = pp_bus_to_ls(pp_net.ext_grid["bus"].values, pp_to_ls)
         if pp_net.ext_grid.shape[0] >= 1:
             warnings.warn("LightSim will not consider the pandapower \"ext_grid\" as there "
                           "are already generators tagged as slack bus")
@@ -59,15 +61,15 @@ def _aux_add_slack(model, pp_net):
                       "I will attempt to add some from the ext_grid.")
         # first i try to see if a generator is connected to a slack bus
         # TODO SLACK: deactivate warnings (will be done at the end when more tests will be done)
-        slack_bus_ids = pp_net.ext_grid["bus"].values
+        slack_bus_ids = pp_bus_to_ls(pp_net.ext_grid["bus"].values, pp_to_ls)
         # if pp_net.ext_grid.shape[0] >= 2:
         #     warnings.warn("LightSim cannot handle multiple slack bus at the moment. Only the first "
         #                   "slack bus of pandapower will be used.")
 
-        if np.all(np.isin(slack_bus_ids, pp_net.gen["bus"].values)):
+        if np.all(np.isin(slack_bus_ids, pp_bus_to_ls(pp_net.gen["bus"].values, pp_to_ls))):
             # all slack buses have a generator connected to them
             # so i assume it was just a computation artifact, and assign these generators as slack buses
-            slack_gen_ids = np.isin(pp_net.gen["bus"].values, slack_bus_ids)  # id of generators connected to slack bus
+            slack_gen_ids = np.isin(pp_bus_to_ls(pp_net.gen["bus"].values, pp_to_ls), slack_bus_ids)  # id of generators connected to slack bus
             slack_gen_ids = np.where(slack_gen_ids)[0]  # keep only the id of the generators
             if "slack_weight" in pp_net.gen:
                 slack_coeff = pp_net.gen["slack_weight"].values[slack_gen_ids]
@@ -100,7 +102,7 @@ def _aux_add_slack(model, pp_net):
             vm_pu = 1.0 * pp_net.ext_grid["vm_pu"].values
             gen_p = np.concatenate((pp_net.gen["p_mw"].values, slack_contrib))
             gen_v = np.concatenate((pp_net.gen["vm_pu"].values, vm_pu))
-            gen_bus = np.concatenate((pp_net.gen["bus"].values, slack_bus_ids))
+            gen_bus = np.concatenate((pp_bus_to_ls(pp_net.gen["bus"].values, pp_to_ls), slack_bus_ids))
             gen_min_q = np.concatenate((pp_net.gen["min_q_mvar"].values, [-999999.]))
             gen_max_q = np.concatenate((pp_net.gen["max_q_mvar"].values, [+99999.]))
             model.init_generators(gen_p, gen_v, gen_min_q, gen_max_q, gen_bus)
diff --git a/lightsim2grid/gridmodel/_aux_add_storage.py b/lightsim2grid/gridmodel/_aux_add_storage.py
index 4ebf47b..57efe28 100644
--- a/lightsim2grid/gridmodel/_aux_add_storage.py
+++ b/lightsim2grid/gridmodel/_aux_add_storage.py
@@ -1,14 +1,16 @@
-# Copyright (c) 2020, RTE (https://www.rte-france.com)
+# Copyright (c) 2020-2023, RTE (https://www.rte-france.com)
 # See AUTHORS.txt
 # This Source Code Form is subject to the terms of the Mozilla Public License, version 2.0.
 # If a copy of the Mozilla Public License, version 2.0 was not distributed with this file,
 # you can obtain one at http://mozilla.org/MPL/2.0/.
 # SPDX-License-Identifier: MPL-2.0
 # This file is part of LightSim2grid, LightSim2grid implements a c++ backend targeting the Grid2Op platform.
+
 import numpy as np
+from ._pp_bus_to_ls_bus import pp_bus_to_ls
 
 
-def _aux_add_storage(model, pp_net):
+def _aux_add_storage(model, pp_net, pp_to_ls):
     """
     Add the storages of the pp_net into the lightsim2grid "model"
 
@@ -28,7 +30,7 @@ def _aux_add_storage(model, pp_net):
 
     model.init_storages(pp_net.storage["p_mw"].values,
                         pp_net.storage["q_mvar"].values,
-                        pp_net.storage["bus"].values
+                        pp_bus_to_ls(pp_net.storage["bus"].values, pp_to_ls)
                         )
     for stor_id, is_connected in enumerate(pp_net.storage["in_service"].values):
         if not is_connected:
diff --git a/lightsim2grid/gridmodel/_aux_add_trafo.py b/lightsim2grid/gridmodel/_aux_add_trafo.py
index 21a6753..f6aa299 100644
--- a/lightsim2grid/gridmodel/_aux_add_trafo.py
+++ b/lightsim2grid/gridmodel/_aux_add_trafo.py
@@ -1,15 +1,17 @@
-# Copyright (c) 2020, RTE (https://www.rte-france.com)
+# Copyright (c) 2020-2023, RTE (https://www.rte-france.com)
 # See AUTHORS.txt
 # This Source Code Form is subject to the terms of the Mozilla Public License, version 2.0.
 # If a copy of the Mozilla Public License, version 2.0 was not distributed with this file,
 # you can obtain one at http://mozilla.org/MPL/2.0/.
 # SPDX-License-Identifier: MPL-2.0
 # This file is part of LightSim2grid, LightSim2grid implements a c++ backend targeting the Grid2Op platform.
+
 import warnings
 import numpy as np
+from ._pp_bus_to_ls_bus import pp_bus_to_ls
 
 
-def _aux_add_trafo(converter, model, pp_net):
+def _aux_add_trafo(converter, model, pp_net, pp_to_ls):
     """
     Add the transformers of the pp_net into the lightsim2grid "model"
 
@@ -89,8 +91,8 @@ def _aux_add_trafo(converter, model, pp_net):
                      tap_pos,
                      shift_,
                      is_tap_hv_side,
-                     pp_net.trafo["hv_bus"].values,
-                     pp_net.trafo["lv_bus"].values)
+                     pp_bus_to_ls(pp_net.trafo["hv_bus"].values, pp_to_ls),
+                     pp_bus_to_ls(pp_net.trafo["lv_bus"].values, pp_to_ls))
 
     for tr_id, is_connected in enumerate(pp_net.trafo["in_service"].values):
         if not is_connected:
diff --git a/lightsim2grid/gridmodel/_aux_check_legit.py b/lightsim2grid/gridmodel/_aux_check_legit.py
index c7bf09b..d6da02a 100644
--- a/lightsim2grid/gridmodel/_aux_check_legit.py
+++ b/lightsim2grid/gridmodel/_aux_check_legit.py
@@ -1,10 +1,11 @@
-# Copyright (c) 2020, RTE (https://www.rte-france.com)
+# Copyright (c) 2020-2023, RTE (https://www.rte-france.com)
 # See AUTHORS.txt
 # This Source Code Form is subject to the terms of the Mozilla Public License, version 2.0.
 # If a copy of the Mozilla Public License, version 2.0 was not distributed with this file,
 # you can obtain one at http://mozilla.org/MPL/2.0/.
 # SPDX-License-Identifier: MPL-2.0
 # This file is part of LightSim2grid, LightSim2grid implements a c++ backend targeting the Grid2Op platform.
+
 import warnings
 import numpy as np
 
@@ -47,10 +48,10 @@ def _aux_check_legit(pp_net):
     #                        "in pandapower network")
 
     # bus indexes should start at 0 and be contiguous
-    if np.any(np.sort(pp_net.bus.index) != np.arange(pp_net.bus.shape[0])):
-        raise RuntimeError("In order to work, pandapower bus indexes should start at 0 and be contiguous. "
-                           "Make sure that `pp_net.bus.index` have this property. You can write a github "
-                           "issue if you want improvment on this regard.")
+    # if np.any(np.sort(pp_net.bus.index) != np.arange(pp_net.bus.shape[0])):
+    #     raise RuntimeError("In order to work, pandapower bus indexes should start at 0 and be contiguous. "
+    #                        "Make sure that `pp_net.bus.index` have this property. You can write a github "
+    #                        "issue if you want improvment on this regard.")
     
     if "_options" in pp_net and \
        "trafo_model" in pp_net["_options"] and \
diff --git a/lightsim2grid/gridmodel/_pp_bus_to_ls_bus.py b/lightsim2grid/gridmodel/_pp_bus_to_ls_bus.py
new file mode 100644
index 0000000..2f6ebd2
--- /dev/null
+++ b/lightsim2grid/gridmodel/_pp_bus_to_ls_bus.py
@@ -0,0 +1,16 @@
+# Copyright (c) 2023, RTE (https://www.rte-france.com)
+# See AUTHORS.txt
+# This Source Code Form is subject to the terms of the Mozilla Public License, version 2.0.
+# If a copy of the Mozilla Public License, version 2.0 was not distributed with this file,
+# you can obtain one at http://mozilla.org/MPL/2.0/.
+# SPDX-License-Identifier: MPL-2.0
+# This file is part of LightSim2grid, LightSim2grid implements a c++ backend targeting the Grid2Op platform.
+
+import numpy as np
+
+def pp_bus_to_ls(pp_bus_id, pp_to_ls_converter):
+    if pp_to_ls_converter is None:
+        res = pp_bus_id
+    else:
+        res = np.array([pp_to_ls_converter[pp_id] for pp_id in pp_bus_id])
+    return res
diff --git a/lightsim2grid/gridmodel/from_pypowsybl.py b/lightsim2grid/gridmodel/from_pypowsybl.py
new file mode 100644
index 0000000..3afc031
--- /dev/null
+++ b/lightsim2grid/gridmodel/from_pypowsybl.py
@@ -0,0 +1,160 @@
+# Copyright (c) 2023, RTE (https://www.rte-france.com)
+# See AUTHORS.txt
+# This Source Code Form is subject to the terms of the Mozilla Public License, version 2.0.
+# If a copy of the Mozilla Public License, version 2.0 was not distributed with this file,
+# you can obtain one at http://mozilla.org/MPL/2.0/.
+# SPDX-License-Identifier: MPL-2.0
+# This file is part of LightSim2grid, LightSim2grid implements a c++ backend targeting the Grid2Op platform.
+
+import numpy as np
+import pypowsybl as pypo
+from lightsim2grid_cpp import GridModel
+
+
+def init(net : pypo.network, gen_slack_id: int = None):
+    model = GridModel()
+    # for substation
+    # network.get_voltage_levels()["substation_id"]
+    # network.get_substations()
+    # network.get_busbar_sections()
+    
+    
+    # initialize and use converters
+    sn_mva_ = 100.  # TODO read from net
+    f_hz = 50.   # TODO read from net
+    
+    # assign unique id to the buses
+    bus_df = net.get_buses().copy()
+    bus_df["bus_id"] = np.arange(bus_df.shape[0])
+    model.set_sn_mva(sn_mva_)
+    model.set_init_vm_pu(1.06)
+    model.init_bus(net.get_voltage_levels().loc[bus_df["voltage_level_id"].values]["nominal_v"].values,
+                   0, 0  # unused
+                   )
+        
+    # do the generators
+    df_gen = net.get_generators()
+    model.init_generators(df_gen["target_p"].values,
+                          df_gen["target_v"].values / net.get_voltage_levels().loc[df_gen["voltage_level_id"].values]["nominal_v"].values,
+                          df_gen["min_q"].values,
+                          df_gen["max_q"].values,
+                          1 * bus_df.loc[df_gen["bus_id"].values]["bus_id"].values
+                          )
+    # TODO dist slack
+    if gen_slack_id is None:
+        model.add_gen_slackbus(0, 1.)
+    else:
+        model.add_gen_slackbus(gen_slack_id, 1.)
+    
+    # for loads
+    df_load = net.get_loads()
+    model.init_loads(df_load["p0"].values,
+                     df_load["q0"].values,
+                     1 * bus_df.loc[df_load["bus_id"].values]["bus_id"].values
+                     )
+    
+    # for lines
+    df_line = net.get_lines()
+    # TODO add g1 / b1 and g2 / b2 in lightsim2grid
+    # line_h = (1j*df_line["g1"].values + df_line["b1"].values + 1j*df_line["g2"].values + df_line["b2"].values)
+    # per unit
+    branch_from_kv = net.get_voltage_levels().loc[df_line["voltage_level1_id"].values]["nominal_v"].values
+    branch_to_kv = net.get_voltage_levels().loc[df_line["voltage_level2_id"].values]["nominal_v"].values
+    
+    # only valid for lines with same voltages at both side...
+    # branch_from_pu = branch_from_kv * branch_from_kv / sn_mva_
+    # line_r = df_line["r"].values / branch_from_pu
+    # line_x = df_line["x"].values / branch_from_pu  
+    # line_h_or = (1j*df_line["g1"].values + df_line["b1"].values) * branch_from_pu
+    # line_h_ex = (1j*df_line["g2"].values + df_line["b2"].values) * branch_from_pu
+    # real per unit conversion 
+    # see https://github.com/powsybl/pypowsybl/issues/642
+    # see https://github.com/powsybl/powsybl-core/blob/266442cbbd84f630acf786018618eaa3d496c6ba/ieee-cdf/ieee-cdf-converter/src/main/java/com/powsybl/ieeecdf/converter/IeeeCdfImporter.java#L347
+    # for right formula
+    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)
+    model.init_powerlines_full(line_r,
+                               line_x,
+                               line_h_or,
+                               line_h_ex,
+                               1 * bus_df.loc[df_line["bus1_id"].values]["bus_id"].values,
+                               1 * bus_df.loc[df_line["bus2_id"].values]["bus_id"].values
+                              )
+            
+    # for trafo
+    df_trafo = net.get_2_windings_transformers()
+    # TODO net.get_ratio_tap_changers()
+    # TODO net.get_phase_tap_changers()
+    shift_ = np.zeros(df_trafo.shape[0])
+    tap_pos = 1.0 * shift_
+    is_tap_hv_side = np.ones(df_trafo.shape[0], dtype=bool)  # TODO
+    
+    # per unit
+    trafo_from_kv = net.get_voltage_levels().loc[df_trafo["voltage_level1_id"].values]["nominal_v"].values
+    trafo_to_kv = net.get_voltage_levels().loc[df_trafo["voltage_level2_id"].values]["nominal_v"].values
+    trafo_to_pu = trafo_to_kv * trafo_to_kv / sn_mva_
+    # tap
+    tap_step_pct = (df_trafo["rated_u1"] / trafo_from_kv - 1.) * 100.
+    has_tap = tap_step_pct != 0.
+    tap_pos[has_tap] += 1
+    tap_step_pct[~has_tap] = 1.0  # or any other values...
+    model.init_trafo(df_trafo["r"].values / trafo_to_pu,
+                     df_trafo["x"].values / trafo_to_pu,
+                     2.*(1j*df_trafo["g"].values + df_trafo["b"].values) * trafo_to_pu,
+                     tap_step_pct,
+                     tap_pos,
+                     shift_,
+                     is_tap_hv_side,
+                     1 * bus_df.loc[df_trafo["bus1_id"].values]["bus_id"].values, # TODO do I need to change hv / lv
+                     1 * bus_df.loc[df_trafo["bus2_id"].values]["bus_id"].values)    
+    
+    # for shunt
+    df_shunt = net.get_shunt_compensators()
+    shunt_kv = net.get_voltage_levels().loc[df_shunt["voltage_level_id"].values]["nominal_v"].values
+    model.init_shunt(-df_shunt["g"].values * shunt_kv**2,
+                     -df_shunt["b"].values * shunt_kv**2,
+                     1 * bus_df.loc[df_shunt["bus_id"].values]["bus_id"].values
+                    )
+    
+    # for hvdc (TODO not tested yet)
+    df_dc = net.get_hvdc_lines()
+    df_sations = net.get_vsc_converter_stations()
+    bus_from_id = df_sations.loc[df_dc["converter_station1_id"].values]["bus_id"].values
+    bus_to_id = df_sations.loc[df_dc["converter_station2_id"].values]["bus_id"].values
+    loss_percent = np.zeros(df_dc.shape[0])  # TODO 
+    loss_mw = np.zeros(df_dc.shape[0])  # TODO
+    model.init_dclines(bus_df.loc[bus_from_id]["bus_id"].values,
+                       bus_df.loc[bus_to_id]["bus_id"].values,
+                       df_dc["target_p"].values,
+                       loss_percent,
+                       loss_mw,
+                       net.get_voltage_levels().loc[df_sations.loc[df_dc["converter_station1_id"].values]["voltage_level_id"].values]["nominal_v"].values,
+                       net.get_voltage_levels().loc[df_sations.loc[df_dc["converter_station2_id"].values]["voltage_level_id"].values]["nominal_v"].values,
+                       df_sations.loc[df_dc["converter_station1_id"].values]["min_q"].values,
+                       df_sations.loc[df_dc["converter_station1_id"].values]["max_q"].values,
+                       df_sations.loc[df_dc["converter_station2_id"].values]["min_q"].values,
+                       df_sations.loc[df_dc["converter_station2_id"].values]["max_q"].values
+                       )
+    
+    # storage units  (TODO not tested yet)
+    df_batt = net.get_batteries()
+    model.init_storages(df_batt["target_p"].values,
+                        df_batt["target_q"].values,
+                        1 * bus_df.loc[df_batt["bus_id"].values]["bus_id"].values
+                        )
+    
+    # TODO
+    # sgen
+    
+    # TODO checks
+    # no 3windings trafo and other exotic stuff
+    return model
diff --git a/lightsim2grid/gridmodel/initGridModel.py b/lightsim2grid/gridmodel/initGridModel.py
index e13f0f5..8a6840c 100644
--- a/lightsim2grid/gridmodel/initGridModel.py
+++ b/lightsim2grid/gridmodel/initGridModel.py
@@ -80,6 +80,11 @@ def init(pp_net):
     model.set_sn_mva(pp_net.sn_mva)
 
     tmp_bus_ind = np.argsort(pp_net.bus.index)
+    if np.any(np.sort(pp_net.bus.index) != np.arange(pp_net.bus.shape[0])):
+        model._ls_to_pp = 1 * pp_net.bus.index.values.astype(int)
+        pp_to_ls = {pp_bus: ls_bus for pp_bus, ls_bus in zip(pp_net.bus.index, tmp_bus_ind)}
+    else:
+        pp_to_ls = None
     model.init_bus(pp_net.bus.iloc[tmp_bus_ind]["vn_kv"].values,
                    pp_net.line.shape[0],
                    pp_net.trafo.shape[0])
@@ -87,33 +92,37 @@ def init(pp_net):
     # deactivate in lightsim the deactivated bus in pandapower
     for bus_id in range(pp_net.bus.shape[0]):
         if not pp_net.bus["in_service"].values[bus_id]:
-            model.deactivate_bus(bus_id)
+            if pp_to_ls is None:
+                pp_bus_id = bus_id
+            else:
+                pp_bus_id = pp_to_ls[bus_id]
+            model.deactivate_bus(pp_bus_id)
 
     # init the powerlines
-    _aux_add_line(converter, model, pp_net)
+    _aux_add_line(converter, model, pp_net, pp_to_ls)
 
     # init the shunts
-    _aux_add_shunt(model, pp_net)
+    _aux_add_shunt(model, pp_net, pp_to_ls)
 
     # handle the trafos
-    _aux_add_trafo(converter, model, pp_net)
+    _aux_add_trafo(converter, model, pp_net, pp_to_ls)
 
     # handle loads
-    _aux_add_load(model, pp_net)
+    _aux_add_load(model, pp_net, pp_to_ls)
 
     # handle static generators (PQ generator)
-    _aux_add_sgen(model, pp_net)
+    _aux_add_sgen(model, pp_net, pp_to_ls)
 
     # handle generators
-    _aux_add_gen(model, pp_net)
+    _aux_add_gen(model, pp_net, pp_to_ls)
 
     # handle storage units
-    _aux_add_storage(model, pp_net)
+    _aux_add_storage(model, pp_net, pp_to_ls)
 
     # handle dc line
-    _aux_add_dc_line(model, pp_net)
+    _aux_add_dc_line(model, pp_net, pp_to_ls)
 
     # deal with slack bus
-    _aux_add_slack(model, pp_net)
+    _aux_add_slack(model, pp_net, pp_to_ls)
 
     return model
diff --git a/lightsim2grid/solver/__init__.py b/lightsim2grid/solver/__init__.py
index a3f495b..6f4909e 100644
--- a/lightsim2grid/solver/__init__.py
+++ b/lightsim2grid/solver/__init__.py
@@ -13,7 +13,9 @@
            "GaussSeidelSynchSolver",
            "SparseLUSolver",
            "SparseLUSolverSingleSlack",
-           "DCSolver"]
+           "DCSolver",
+           "FDPF_XB_SparseLUSolver",
+           "FDPF_BX_SparseLUSolver"]
 
 from lightsim2grid_cpp import SolverType
 from lightsim2grid_cpp import ErrorType
@@ -24,14 +26,20 @@
 from lightsim2grid_cpp import SparseLUSolver  # SolverType.SparseLU
 from lightsim2grid_cpp import SparseLUSolverSingleSlack  # SolverType.SparseLUSingleSlack
 from lightsim2grid_cpp import DCSolver  # SolverType.DC
+from lightsim2grid_cpp import FDPF_XB_SparseLUSolver  # SolverType.FDPF_XB_SparseLU
+from lightsim2grid_cpp import FDPF_BX_SparseLUSolver  # SolverType.FDPF_BX_SparseLU
 
 try:
     from lightsim2grid_cpp import KLUSolver  # SolverType.KLU
     from lightsim2grid_cpp import KLUSolverSingleSlack  # SolverType.KLUSingleSlack
     from lightsim2grid_cpp import KLUDCSolver  # SolverType.KLUDC
+    from lightsim2grid_cpp import FDPF_XB_KLUSolver  # SolverType.FDPF_XB_KLU
+    from lightsim2grid_cpp import FDPF_BX_KLUSolver  # SolverType.FDPF_BX_KLU
     __all__.append("KLUSolver")
     __all__.append("KLUSolverSingleSlack")
     __all__.append("KLUDCSolver")
+    __all__.append("FDPF_XB_KLUSolver")
+    __all__.append("FDPF_BX_KLUSolver")
 except Exception as exc_:
     # KLU is not available
     pass
@@ -40,9 +48,13 @@
     from lightsim2grid_cpp import NICSLUSolver  # SolverType.NICSLU
     from lightsim2grid_cpp import NICSLUSolverSingleSlack  # SolverType.NICSLUSingleSlack
     from lightsim2grid_cpp import NICSLUDCSolver  # SolverType.NICSLUDC
+    from lightsim2grid_cpp import FDPF_XB_NICSLUSolver  # SolverType.FDPF_XB_NICSLU
+    from lightsim2grid_cpp import FDPF_BX_NICSLUSolver  # SolverType.FDPF_BX_NICSLU
     __all__.append("NICSLUSolver")
     __all__.append("NICSLUSolverSingleSlack")
     __all__.append("NICSLUDCSolver")
+    __all__.append("FDPF_XB_NICSLUSolver")
+    __all__.append("FDPF_BX_NICSLUSolver")
 except Exception as exc_:
     # NICSLU is not available
     pass
@@ -51,9 +63,13 @@
     from lightsim2grid_cpp import CKTSOSolver  # SolverType.CKTSO
     from lightsim2grid_cpp import CKTSOSolverSingleSlack  # SolverType.CKTSOSingleSlack
     from lightsim2grid_cpp import CKTSODCSolver  # SolverType.CKTSODC
+    from lightsim2grid_cpp import FDPF_XB_CKTSOSolver  # SolverType.FDPF_XB_CKTSO
+    from lightsim2grid_cpp import FDPF_BX_CKTSOSolver  # SolverType.FDPF_BX_CKTSO
     __all__.append("CKTSOSolver")
     __all__.append("CKTSOSolverSingleSlack")
     __all__.append("CKTSODCSolver")
+    __all__.append("FDPF_XB_CKTSOSolver")
+    __all__.append("FDPF_BX_CKTSOSolver")
 except Exception as exc_:
     # NICSLU is not available
     pass
diff --git a/setup.py b/setup.py
index 2025432..b28c51d 100644
--- a/setup.py
+++ b/setup.py
@@ -14,7 +14,7 @@
 from pybind11.setup_helpers import Pybind11Extension, build_ext
 
 
-__version__ = "0.7.4"
+__version__ = "0.7.5"
 KLU_SOLVER_AVAILABLE = False
 
 # Try to link against SuiteSparse (if available)
@@ -156,7 +156,8 @@
              "src/BaseSolver.cpp",
              "src/BaseMultiplePowerflow.cpp",
              "src/Computers.cpp",
-             "src/SecurityAnalysis.cpp"]
+             "src/SecurityAnalysis.cpp",
+             "src/Solvers.cpp"]
 
 if KLU_SOLVER_AVAILABLE:
     src_files.append("src/KLUSolver.cpp")
@@ -333,6 +334,7 @@
             "autodocsumm>=0.1.13",
             "grid2op>=1.6.4",
             "recommonmark",
+            "pypowsybl"
         ],
         "benchmark": [
             "tabulate",
@@ -349,7 +351,8 @@
             "grid2op>=1.6.4",
             "numba",
             "pandapower>=2.8.0",
-            "packaging"
+            "packaging", 
+            "pypowsybl"
         ]
     }
 }
diff --git a/src/BaseConstants.h b/src/BaseConstants.h
index c728fed..7d6f268 100644
--- a/src/BaseConstants.h
+++ b/src/BaseConstants.h
@@ -26,4 +26,8 @@ class BaseConstants
         static const real_type my_zero_;
 };
 
+enum class FDPFMethod {XB, BX};  // Different type of FDPF powerflow
+// FDPFMethod::XB => alg = 2 in pypower / pandapower
+// FDPFMethod::BX => alg = 3 in pypower / pandapower
+
 #endif // BASECONSTANTS_H
diff --git a/src/BaseFDPFSolver.h b/src/BaseFDPFSolver.h
new file mode 100644
index 0000000..4fa95e2
--- /dev/null
+++ b/src/BaseFDPFSolver.h
@@ -0,0 +1,221 @@
+// Copyright (c) 2023, RTE (https://www.rte-france.com)
+// See AUTHORS.txt
+// This Source Code Form is subject to the terms of the Mozilla Public License, version 2.0.
+// If a copy of the Mozilla Public License, version 2.0 was not distributed with this file,
+// you can obtain one at http://mozilla.org/MPL/2.0/.
+// SPDX-License-Identifier: MPL-2.0
+// This file is part of LightSim2grid, LightSim2grid implements a c++ backend targeting the Grid2Op platform.
+
+#ifndef BASEFDPFSOLVER_H
+#define BASEFDPFSOLVER_H
+
+#include "BaseSolver.h"
+
+/**
+Base class for Fast Decoupled Powerflow based solver
+**/
+template<class LinearSolver, FDPFMethod XB_BX>
+class BaseFDPFSolver : public BaseSolver
+{
+    public:
+        BaseFDPFSolver():BaseSolver(true), need_factorize_(true) {}
+
+        virtual
+        bool compute_pf(const Eigen::SparseMatrix<cplx_type> & Ybus,
+                        CplxVect & V,
+                        const CplxVect & Sbus,
+                        const Eigen::VectorXi & slack_ids,
+                        const RealVect & slack_weights,
+                        const Eigen::VectorXi & pv,
+                        const Eigen::VectorXi & pq,
+                        int max_iter,
+                        real_type tol
+                        ) ;  // requires a gridmodel !
+
+        // bool compute_pf(const Eigen::SparseMatrix<cplx_type> & Ybus,
+        //                 const Eigen::SparseMatrix<cplx_type> & Bp,
+        //                 const Eigen::SparseMatrix<cplx_type> & Bpp,
+        //                 CplxVect & V,
+        //                 const CplxVect & Sbus,
+        //                 const Eigen::VectorXi & slack_ids,
+        //                 const RealVect & slack_weights,
+        //                 const Eigen::VectorXi & pv,
+        //                 const Eigen::VectorXi & pq,
+        //                 int max_iter,
+        //                 real_type tol
+        //                 ) ;  // TODO add Bp and Bpp as argument for use in python directly !
+
+        virtual void reset()
+        {   
+            BaseSolver::reset();
+            // solution of the problem
+            Bp_ = Eigen::SparseMatrix<real_type> ();  // the B prime matrix (size n_pvpq)
+            Bpp_ = Eigen::SparseMatrix<real_type>();  // the B double prime matrix  (size n_pq)
+            p_ = RealVect();
+            q_ = RealVect();
+            need_factorize_ = true;
+
+            // reset linear solvers
+            ErrorType reset_status = _linear_solver_Bp.reset();
+            if(reset_status != ErrorType::NoError) err_ = reset_status;
+            reset_status = _linear_solver_Bpp.reset();
+            if((reset_status != ErrorType::NoError) && (err_ != ErrorType::NotInitError)) err_ = reset_status;
+        }
+
+        Eigen::SparseMatrix<real_type> debug_get_Bp_python() { return Bp_;}
+        Eigen::SparseMatrix<real_type> debug_get_Bpp_python() { return Bpp_;}
+
+    protected:
+        virtual void reset_timer(){
+            BaseSolver::reset_timer();
+        }
+
+        CplxVect evaluate_mismatch(const Eigen::SparseMatrix<cplx_type> &  Ybus,
+                                   const CplxVect & V,
+                                   const CplxVect & Sbus,
+                                   Eigen::Index slack_id,  // id of the ref slack bus
+                                   real_type slack_absorbed,
+                                   const RealVect & slack_weights)
+        {
+            CplxVect tmp = Ybus * V;  // this is a vector
+            tmp = tmp.array().conjugate();  // i take the conjugate
+            auto mis = V.array() * tmp.array() - Sbus.array() + slack_absorbed * slack_weights.array();
+            return mis;
+        }
+        
+        void fillBp_Bpp(Eigen::SparseMatrix<real_type> & Bp, Eigen::SparseMatrix<real_type> & Bpp) const;  // defined in Solvers.cpp !
+
+        virtual
+        void initialize(){
+            auto timer = CustTimer();
+            err_ = ErrorType::NoError; // reset error message
+            // init Bp solver
+            ErrorType init_status = _linear_solver_Bp.initialize(Bp_);
+            if(init_status != ErrorType::NoError){
+                _linear_solver_Bp.reset();
+                _linear_solver_Bpp.reset();
+                err_ = init_status;
+                need_factorize_ = true;
+                timer_initialize_ += timer.duration();
+                return;
+            }
+            // init Bpp solver (if Bp is sucesfull of course)
+            init_status = _linear_solver_Bpp.initialize(Bpp_);
+            if(init_status != ErrorType::NoError){
+                _linear_solver_Bp.reset();
+                _linear_solver_Bpp.reset();
+                err_ = init_status;
+                need_factorize_ = true;
+                timer_initialize_ += timer.duration();
+                return;
+            }
+
+            // everything went well
+            need_factorize_ = false;
+            timer_initialize_ += timer.duration();
+        }
+
+        virtual
+        void solve(LinearSolver& linear_solver,
+                   Eigen::SparseMatrix<real_type>& mat,
+                   RealVect & b,
+                   bool has_just_been_inialized){
+            auto timer = CustTimer();
+            const ErrorType solve_status = linear_solver.solve(mat, b, has_just_been_inialized);
+            if(solve_status != ErrorType::NoError){
+                // std::cout << "solve error: " << solve_status << std::endl;
+                err_ = solve_status;
+            }
+            timer_solve_ += timer.duration();
+        }
+
+        bool has_converged(const Eigen::Ref<const CplxVect > & tmp_va,
+                           const Eigen::SparseMatrix<cplx_type> & Ybus,
+                           const CplxVect & Sbus,
+                           Eigen::Index slack_bus_id,
+                           real_type & slack_absorbed,
+                           const RealVect & slack_weights,
+                           const Eigen::Ref<const Eigen::VectorXi> & pvpq,
+                           const Eigen::Ref<const Eigen::VectorXi> & pq,
+                           real_type tol)
+        {
+            /**
+            It is suppose to implement the python code bellow (so it updates p_, q_, v_, va_ and vm_):
+
+            .. code-block:: python
+
+                V = Vm * exp(1j * Va)
+
+                ## evalute mismatch
+                mis = (V * conj(Ybus * V) - Sbus) / Vm
+                P = mis[pvpq].real
+                Q = mis[pq].imag
+
+                ## check tolerance
+                normP = linalg.norm(P, Inf)
+                normQ = linalg.norm(Q, Inf)
+                if verbose > 1:
+                    sys.stdout.write('\n  Q  %3d   %10.3e   %10.3e' % (i, normP, normQ))
+                if normP < tol and normQ < tol:
+                    converged = 1
+                    if verbose:
+                        sys.stdout.write('\nFast-decoupled power flow converged in %d '
+                            'P-iterations and %d Q-iterations.\n' % (i, i))
+                    break
+            **/
+
+            // V = Vm * exp(1j * Va) : 
+            V_ = Vm_.array() * tmp_va.array();
+            Vm_ = V_.array().abs();
+            Va_ = V_.array().arg();
+            auto mis = evaluate_mismatch(Ybus, V_, Sbus, slack_bus_id, slack_absorbed, slack_weights);  // mis = (V * conj(Ybus * V) - Sbus) / Vm
+            mis.array() /= Vm_.array();  // mis = (V * conj(Ybus * V) - Sbus) / Vm (do not forget the / Vm !)
+
+            bool tmp = mis.allFinite();
+            if(!tmp){
+                err_ = ErrorType::InifiniteValue;
+                return false; // divergence due to Nans
+            }
+            p_ = mis(pvpq).real();  // P = mis[pvpq].real
+            q_ = mis(pq).imag();  // Q = mis[pq].imag
+            return _check_for_convergence(p_, q_, tol);
+        }
+
+        void fill_sparse_matrices(const Eigen::SparseMatrix<real_type> & grid_Bp,
+                                  const Eigen::SparseMatrix<real_type> & grid_Bpp,
+                                  const std::vector<int> & pvpq_inv,
+                                  const std::vector<int> & pq_inv,
+                                  Eigen::Index n_pvpq,
+                                  Eigen::Index n_pq);
+
+        void aux_fill_sparse_matrices(const Eigen::SparseMatrix<real_type> & grid_Bp_Bpp,
+                                      const std::vector<int> & ind_inv,
+                                      Eigen::Index mat_dim,
+                                      Eigen::SparseMatrix<real_type> & res);
+
+    protected:
+        // use 2 linear solvers
+        LinearSolver _linear_solver_Bp;
+        LinearSolver _linear_solver_Bpp;
+
+        // solution of the problem
+        Eigen::SparseMatrix<real_type> Bp_;  // the B prime matrix (size n_pvpq)
+        Eigen::SparseMatrix<real_type> Bpp_;  // the B double prime matrix  (size n_pq)
+        RealVect p_;  // (size n_pvpq)
+        RealVect q_;  // (size n_pq)
+        bool need_factorize_;
+
+        // timers
+        double timer_initialize_;
+        // double timer_dSbus_;
+        // double timer_fillJ_;
+
+    private:
+        // no copy allowed
+        BaseFDPFSolver( const BaseFDPFSolver & ) =delete ;
+        BaseFDPFSolver & operator=( const BaseFDPFSolver & ) =delete ;
+};
+
+#include "BaseFDPFSolver.tpp"
+
+#endif // BASEFDPFSOLVER_H
diff --git a/src/BaseFDPFSolver.tpp b/src/BaseFDPFSolver.tpp
new file mode 100644
index 0000000..080d191
--- /dev/null
+++ b/src/BaseFDPFSolver.tpp
@@ -0,0 +1,190 @@
+// Copyright (c) 2023, RTE (https://www.rte-france.com)
+// See AUTHORS.txt
+// This Source Code Form is subject to the terms of the Mozilla Public License, version 2.0.
+// If a copy of the Mozilla Public License, version 2.0 was not distributed with this file,
+// you can obtain one at http://mozilla.org/MPL/2.0/.
+// SPDX-License-Identifier: MPL-2.0
+// This file is part of LightSim2grid, LightSim2grid implements a c++ backend targeting the Grid2Op platform.
+
+// inspired from pypower https://github.com/rwl/PYPOWER/blob/master/pypower/fdpf.py
+
+template<class LinearSolver, FDPFMethod XB_BX>
+bool BaseFDPFSolver<LinearSolver, XB_BX>::compute_pf(const Eigen::SparseMatrix<cplx_type> & Ybus,
+                                                     CplxVect & V,
+                                                     const CplxVect & Sbus,
+                                                     const Eigen::VectorXi & slack_ids,
+                                                     const RealVect & slack_weights,
+                                                     const Eigen::VectorXi & pv,
+                                                     const Eigen::VectorXi & pq,
+                                                     int max_iter,
+                                                     real_type tol
+                                                     )
+{
+    /**
+    This method uses the newton raphson algorithm to compute voltage angles and magnitudes at each bus
+    of the system.
+    If the Ybus matrix changed, please uses the appropriate method to recomptue it!
+    **/
+    // TODO DEBUG MODE check what can be checked: no voltage at 0, Ybus is square, Sbus same size than V and
+    // TODO DEBUG MODE Ybus (nrow or ncol), pv and pq have value that are between 0 and nrow etc.
+    // TODO DEBUG MODE: check that all nodes is either pv, pq or slack
+    // TODO DEBUG MODE: check that the slack_weights sum to 1.0
+    if(Sbus.size() != Ybus.rows() || Sbus.size() != Ybus.cols() ){
+        // TODO DEBUG MODE
+        std::ostringstream exc_;
+        exc_ << "BaseFDPFSolver::compute_pf: Size of the Sbus should be the same as the size of Ybus. Currently: ";
+        exc_ << "Sbus  (" << Sbus.size() << ") and Ybus (" << Ybus.rows() << ", " << Ybus.rows() << ").";
+        throw std::runtime_error(exc_.str());
+    }
+    if(V.size() != Ybus.rows() || V.size() != Ybus.cols() ){
+        // TODO DEBUG MODE
+        std::ostringstream exc_;
+        exc_ << "BaseFDPFSolver::compute_pf: Size of V (init voltages) should be the same as the size of Ybus. Currently: ";
+        exc_ << "V  (" << V.size() << ") and Ybus (" << Ybus.rows()<< ", " << Ybus.rows() << ").";
+        throw std::runtime_error(exc_.str());
+    }
+    reset_timer();
+    auto timer = CustTimer();
+    if(!is_linear_solver_valid()) return false;
+    err_ = ErrorType::NoError;  // reset the error if previous error happened
+
+    Eigen::VectorXi my_pv = retrieve_pv_with_slack(slack_ids, pv);  // retrieve_pv_with_slack (not all), add_slack_to_pv (all)
+    real_type slack_absorbed = std::real(Sbus.sum());  // initial guess for slack_absorbed
+    const auto slack_bus_id = slack_ids(0);
+    
+    // initialize once and for all the "inverse" of these vectors
+    const auto n_pv = my_pv.size();
+    const auto n_pq = pq.size();
+    Eigen::VectorXi pvpq(n_pv + n_pq);
+    pvpq << my_pv, pq;  // pvpq = r_[pv, pq]
+
+    // fill the sparse matrix Bp and Bpp (depends on the method: BX or XB)
+    //TODO     if(need_factorize_) ???
+    Eigen::SparseMatrix<real_type> grid_Bp;
+    Eigen::SparseMatrix<real_type> grid_Bpp;
+    fillBp_Bpp(grid_Bp, grid_Bpp);
+    
+    // fill the solver matrices Bp and Bpp : 
+    // Bp_ = Bp[array([pvpq]).T, pvpq].tocsc() # splu requires a CSC matrix
+    // Bpp_ = Bpp[array([pq]).T, pq].tocsc()
+    const auto n_pvpq = pvpq.size();
+    std::vector<int> pvpq_inv(V.size(), -1);
+    for(int inv_id=0; inv_id < n_pvpq; ++inv_id) pvpq_inv[pvpq(inv_id)] = inv_id;
+    std::vector<int> pq_inv(V.size(), -1);
+    for(int inv_id=0; inv_id < n_pq; ++inv_id) pq_inv[pq(inv_id)] = inv_id;
+    fill_sparse_matrices(grid_Bp, grid_Bpp, pvpq_inv, pq_inv, n_pvpq, n_pq);
+
+    V_ = V;  // V = V0
+    Vm_ = V_.array().abs();   // Vm = abs(V)
+    Va_ = V_.array().arg();  // Va = angle(V)
+
+    // first check, if the problem is already solved, i stop there
+    // compute a first time the mismatch to initialize the slack bus
+    CplxVect mis = evaluate_mismatch(Ybus, V, Sbus, slack_bus_id, slack_absorbed, slack_weights);
+    mis.array() /= Vm_.array();  // mis = (V * conj(Ybus * V) - Sbus) / Vm
+    p_ = mis(pvpq).real();  // P = mis[pvpq].real
+    q_ = mis(pq).imag();  // Q = mis[pq].imag
+    
+    CplxVect tmp_va;
+    nr_iter_ = 0; //current step
+    bool converged = _check_for_convergence(p_, q_, tol);
+    bool res = true;  // have i converged or not
+    bool ls_bp_has_just_been_initialized = false;
+    bool ls_bpp_has_just_been_initialized = false;
+
+    while ((!converged) & (nr_iter_ < max_iter)){
+        nr_iter_++;
+
+        if(need_factorize_){
+            initialize();
+            if(err_ != ErrorType::NoError){
+                // I got an error during the initialization of the linear system, i need to stop here
+                res = false;
+                break;
+            }
+            ls_bp_has_just_been_initialized = true;
+            ls_bpp_has_just_been_initialized = true;
+        }
+
+        // do the P iteration (for Va)
+        RealVect x = p_;
+        solve(_linear_solver_Bp, Bp_, x, ls_bp_has_just_been_initialized);  //  dVa = -Bp_solver.solve(P)     
+        if(err_ != ErrorType::NoError){
+            // I got an error during the solving of the linear system, i need to stop here
+            res = false;
+            break;
+        }
+        Va_(pvpq) -= x;  // Va[pvpq] = Va[pvpq] + dVa
+        tmp_va.array() = (Va_.array().cos().template cast<cplx_type>() + my_i * Va_.array().sin().template cast<cplx_type>() );  // reused for Q iteration
+        if(has_converged(tmp_va, Ybus, Sbus, slack_bus_id, slack_absorbed, slack_weights, pvpq, pq, tol)){
+            converged = true;
+            break;
+        }
+        // do the Q iterations (for Vm)
+        solve(_linear_solver_Bpp, Bpp_, q_, ls_bpp_has_just_been_initialized);  //  dVm = -Bpp_solver.solve(Q)   
+        if(err_ != ErrorType::NoError){
+            // I got an error during the solving of the linear system, i need to stop here
+            res = false;
+            break;
+        }
+        Vm_(pq) -= q_;
+        converged = has_converged(tmp_va, Ybus, Sbus, slack_bus_id, slack_absorbed, slack_weights, pvpq, pq, tol);
+    }
+    if(!converged){
+        if (err_ == ErrorType::NoError) err_ = ErrorType::TooManyIterations;
+        res = false;
+    }
+    timer_total_nr_ += timer.duration();
+    #ifdef __COUT_TIMES
+        std::cout << "Computation time: " << "\n\t timer_initialize_: " << timer_initialize_
+                  << "\n\t timer_dSbus_ (called in _fillJ_): " << timer_dSbus_
+                  << "\n\t timer_fillJ_: " << timer_fillJ_
+                  << "\n\t timer_Fx_: " << timer_Fx_
+                  << "\n\t timer_check_: " << timer_check_
+                  << "\n\t timer_solve_: " << timer_solve_
+                  << "\n\t timer_total_nr_: " << timer_total_nr_
+                  << "\n\n";
+    #endif // __COUT_TIMES
+    return res;
+}
+
+template<class LinearSolver, FDPFMethod XB_BX>
+void BaseFDPFSolver<LinearSolver, XB_BX>::fill_sparse_matrices(const Eigen::SparseMatrix<real_type> & grid_Bp,
+                                                               const Eigen::SparseMatrix<real_type> & grid_Bpp,
+                                                               const std::vector<int> & pvpq_inv,
+                                                               const std::vector<int> & pq_inv,
+                                                               Eigen::Index n_pvpq,
+                                                               Eigen::Index n_pq)
+{
+  /**
+   Init Bp_ and Bpp_ such that:
+    // Bp_ = grid_Bp[array([pvpq]).T, pvpq].tocsc() # splu requires a CSC matrix
+    // Bpp_ = grid_Bpp[array([pq]).T, pq].tocsc()
+   **/
+  aux_fill_sparse_matrices(grid_Bp, pvpq_inv, n_pvpq, Bp_);
+  aux_fill_sparse_matrices(grid_Bpp, pq_inv, n_pq, Bpp_);
+}
+
+template<class LinearSolver, FDPFMethod XB_BX>
+void BaseFDPFSolver<LinearSolver, XB_BX>::aux_fill_sparse_matrices(const Eigen::SparseMatrix<real_type> & grid_Bp_Bpp,
+                                                                   const std::vector<int> & ind_inv,
+                                                                   Eigen::Index mat_dim,
+                                                                   Eigen::SparseMatrix<real_type> & res)
+{
+    // clear previous matrix
+    res = Eigen::SparseMatrix<real_type>(mat_dim, mat_dim);
+    // compute the coefficients
+    std::vector<Eigen::Triplet<real_type> > tripletList;
+    tripletList.reserve(grid_Bp_Bpp.nonZeros());
+    for (int k = 0; k < grid_Bp_Bpp.outerSize(); ++k){
+        for (Eigen::SparseMatrix<real_type>::InnerIterator it(grid_Bp_Bpp, k); it; ++it){
+            if ((ind_inv[it.row()] >= 0) && (ind_inv[it.col()] >= 0)){
+                auto tmp = Eigen::Triplet<real_type>(ind_inv[it.row()], ind_inv[it.col()], it.value());
+                tripletList.push_back(tmp);
+            }
+        }
+    }
+    // put them in the matrix
+    res.setFromTriplets(tripletList.begin(), tripletList.end());
+    res.makeCompressed();
+}
diff --git a/src/BaseMultiplePowerflow.h b/src/BaseMultiplePowerflow.h
index aa02282..34c6fae 100644
--- a/src/BaseMultiplePowerflow.h
+++ b/src/BaseMultiplePowerflow.h
@@ -53,8 +53,8 @@ class BaseMultiplePowerflow
         void change_solver(const SolverType & type){
             _solver.change_solver(type);
         }
-        std::vector<SolverType> available_solvers() {return _solver.available_solvers(); }
-        SolverType get_solver_type() {return _solver.get_type(); }
+        std::vector<SolverType> available_solvers() const {return _solver.available_solvers(); }
+        SolverType get_solver_type() const {return _solver.get_type(); }
 
         // utlities informations
         double amps_computation_time() const {return _timer_compute_A;}
diff --git a/src/BaseNRSolver.tpp b/src/BaseNRSolver.tpp
index accab85..9c25925 100644
--- a/src/BaseNRSolver.tpp
+++ b/src/BaseNRSolver.tpp
@@ -114,7 +114,7 @@ bool BaseNRSolver<LinearSolver>::compute_pf(const Eigen::SparseMatrix<cplx_type>
             Va_(pq) -= F.segment(n_pv + 1, n_pq);
             Vm_(pq) -= F.segment(n_pv + n_pq + 1, n_pq);
         }
-        slack_absorbed -= F(0); // by convention in fill_jacobian_matrix the slack bus is the last component
+        slack_absorbed -= F(0); // by convention in fill_jacobian_matrix the slack bus is the first component
 
         // std::cout << "iter " << nr_iter_ << " dx(0): " << -F(0) << " dx(1): " << -F(1) << std::endl;
         // std::cout << "slack_absorbed " << slack_absorbed << std::endl;
diff --git a/src/BaseNRSolverSingleSlack.tpp b/src/BaseNRSolverSingleSlack.tpp
index 65fc365..df0a6e8 100644
--- a/src/BaseNRSolverSingleSlack.tpp
+++ b/src/BaseNRSolverSingleSlack.tpp
@@ -30,13 +30,13 @@ bool BaseNRSolverSingleSlack<LinearSolver>::compute_pf(const Eigen::SparseMatrix
     // TODO Ybus (nrow or ncol), pv and pq have value that are between 0 and nrow etc.
     if(Sbus.size() != Ybus.rows() || Sbus.size() != Ybus.cols() ){
         std::ostringstream exc_;
-        exc_ << "BaseNRSolver::compute_pf: Size of the Sbus should be the same as the size of Ybus. Currently: ";
+        exc_ << "BaseNRSolverSingleSlack::compute_pf: Size of the Sbus should be the same as the size of Ybus. Currently: ";
         exc_ << "Sbus  (" << Sbus.size() << ") and Ybus (" << Ybus.rows() << ", " << Ybus.rows() << ").";
         throw std::runtime_error(exc_.str());
     }
     if(V.size() != Ybus.rows() || V.size() != Ybus.cols() ){
         std::ostringstream exc_;
-        exc_ << "BaseNRSolver::compute_pf: Size of V (init voltages) should be the same as the size of Ybus. Currently: ";
+        exc_ << "BaseNRSolverSingleSlack::compute_pf: Size of V (init voltages) should be the same as the size of Ybus. Currently: ";
         exc_ << "V  (" << V.size() << ") and Ybus (" << Ybus.rows()<<", "<<Ybus.rows() << ").";
         throw std::runtime_error(exc_.str());
     }
@@ -47,7 +47,8 @@ bool BaseNRSolverSingleSlack<LinearSolver>::compute_pf(const Eigen::SparseMatrix
     BaseNRSolver<LinearSolver>::err_ = ErrorType::NoError;  // reset the error if previous error happened
     auto timer = CustTimer();
     // initialize once and for all the "inverse" of these vectors
-    Eigen::VectorXi my_pv = BaseNRSolver<LinearSolver>::retrieve_pv_with_slack(slack_ids, pv);
+    // Eigen::VectorXi my_pv = BaseNRSolver<LinearSolver>::retrieve_pv_with_slack(slack_ids, pv);
+    Eigen::VectorXi my_pv = pv;
     // Eigen::VectorXi my_pv = pv; // BaseNRSolver<LinearSolver>::retrieve_pv_with_slack(slack_ids, pv);
 
     const int n_pv = static_cast<int>(my_pv.size());
diff --git a/src/BaseSolver.cpp b/src/BaseSolver.cpp
index 108d92f..05c837e 100644
--- a/src/BaseSolver.cpp
+++ b/src/BaseSolver.cpp
@@ -121,9 +121,21 @@ bool BaseSolver::_check_for_convergence(const RealVect & F,
     return res;
 }
 
-int BaseSolver::extract_slack_bus_id(const Eigen::VectorXi & pv,
-                                     const Eigen::VectorXi & pq,
-                                     unsigned int nb_bus)
+bool BaseSolver::_check_for_convergence(const RealVect & p,
+                                        const RealVect & q,
+                                        real_type tol)
+{
+    auto timer = CustTimer();
+    const auto norm_inf_p = p.lpNorm<Eigen::Infinity>();
+    const auto norm_inf_q = q.lpNorm<Eigen::Infinity>();
+    bool res =  (norm_inf_p  < tol) && (norm_inf_q < tol);
+    timer_check_ += timer.duration();
+    return res;
+}
+
+Eigen::VectorXi BaseSolver::extract_slack_bus_id(const Eigen::VectorXi & pv,
+                                                 const Eigen::VectorXi & pq,
+                                                 unsigned int nb_bus)
 {
     // pv: list of index of pv nodes
     // pq: list of index of pq nodes
@@ -131,7 +143,8 @@ int BaseSolver::extract_slack_bus_id(const Eigen::VectorXi & pv,
     // returns: res: the id of the unique slack bus (throw an error if no slack bus is found)
     // /!\ does not support multiple slack bus!!!
 
-    int res=-1;
+    Eigen::VectorXi res(nb_bus - pv.size() - pq.size());
+    Eigen::Index i_res = 0;
     bool found=false;
     // run through both pv and pq nodes and declare they are not slack bus
     std::vector<bool> tmp(nb_bus, true);
@@ -148,15 +161,18 @@ int BaseSolver::extract_slack_bus_id(const Eigen::VectorXi & pv,
     {
         if(tmp[k])
         {
-            if(!found){
-                res = k;
-                found = true;
-            }else{
-                throw std::runtime_error("BaseSolver::extract_slack_bus_id: multiple slack bus found on your grid !");
-            }
+            res[i_res] = k;
+            ++i_res;
+            // if(!found){
+            //     res = k;
+            //     found = true;
+            // }else{
+            //     throw std::runtime_error("BaseSolver::extract_slack_bus_id: multiple slack bus found on your grid !");
+            // }
         }
     }
-    if(res == -1){
+    // if(res == -1){
+    if(res.size() != i_res){
         throw std::runtime_error("BaseSolver::extract_slack_bus_id: No slack bus is found in your grid");
     }
     return res;
diff --git a/src/BaseSolver.h b/src/BaseSolver.h
index d31fddf..54d8710 100644
--- a/src/BaseSolver.h
+++ b/src/BaseSolver.h
@@ -28,6 +28,9 @@
 #include "CustTimer.h"
 #include "BaseConstants.h"
 
+class GridModel;
+
+
 /**
 This class represents a solver to compute powerflow.
 
@@ -51,6 +54,10 @@ class BaseSolver : public BaseConstants
 
         virtual ~BaseSolver(){}
 
+        void set_gridmodel(const GridModel * gridmodel){
+            _gridmodel = gridmodel;
+        }
+
         Eigen::Ref<const RealVect> get_Va() const{
             return Va_;
         }
@@ -126,6 +133,10 @@ class BaseSolver : public BaseConstants
         bool _check_for_convergence(const RealVect & F,
                                     real_type tol);
 
+        bool _check_for_convergence(const RealVect & p,
+                                    const RealVect & q,
+                                    real_type tol);
+
         void one_iter_all_at_once(CplxVect & tmp_Sbus,
                                   const Eigen::SparseMatrix<cplx_type> & Ybus,
                                   const Eigen::VectorXi & pv,
@@ -137,9 +148,9 @@ class BaseSolver : public BaseConstants
                       const Eigen::VectorXi & pq
                       );
 
-        int extract_slack_bus_id(const Eigen::VectorXi & pv,
-                                 const Eigen::VectorXi & pq,
-                                 unsigned int nb_bus);
+        Eigen::VectorXi extract_slack_bus_id(const Eigen::VectorXi & pv,
+                                             const Eigen::VectorXi & pq,
+                                             unsigned int nb_bus);
 
         /**
         When there are multiple slacks, add the other "slack buses" in the PV buses indexes
@@ -171,6 +182,14 @@ class BaseSolver : public BaseConstants
             return my_pv;
         }
         
+        // terribly inefficient way to know if an element is in a vector
+        bool isin(int k, const Eigen::VectorXi vect) const{
+            for(auto el : vect){
+                if(el == k) return true;
+            }
+            return false;
+        }
+
     protected:
         // solver initialization
         int n_;
@@ -190,11 +209,13 @@ class BaseSolver : public BaseConstants
         // 4: end of possible iterations (divergence because nr_iter_ >= max_iter)
 
         // timers
-         double timer_Fx_;
-         double timer_solve_;
-         double timer_check_;
-         double timer_total_nr_;
+        double timer_Fx_;
+        double timer_solve_;
+        double timer_check_;
+        double timer_total_nr_;
 
+        const GridModel * _gridmodel;  // does not have ownership so that's fine (pointer to the base gridmodel, can be used for some powerflow)
+        
     private:
         // no copy allowed
         BaseSolver( const BaseSolver & ) ;
diff --git a/src/ChooseSolver.h b/src/ChooseSolver.h
index 764bb4f..57d5e69 100644
--- a/src/ChooseSolver.h
+++ b/src/ChooseSolver.h
@@ -20,7 +20,12 @@
 enum class SolverType {SparseLU, KLU, GaussSeidel, DC, GaussSeidelSynch, NICSLU, 
                        SparseLUSingleSlack, KLUSingleSlack, NICSLUSingleSlack, 
                        KLUDC, NICSLUDC,
-                       CKTSO, CKTSOSingleSlack, CKTSODC};
+                       CKTSO, CKTSOSingleSlack, CKTSODC,
+                       FDPF_XB_SparseLU, FDPF_BX_SparseLU, // from 0.7.5
+                       FDPF_XB_KLU, FDPF_BX_KLU,  // from 0.7.5
+                       FDPF_XB_NICSLU, FDPF_BX_NICSLU,  // from 0.7.5
+                       FDPF_XB_CKTSO,  FDPF_BX_CKTSO  // from 0.7.5
+                       };
 
 // TODO define a template class instead of these weird stuff !!!
 // TODO export all methods from base class !
@@ -48,20 +53,28 @@ class ChooseSolver
             res.push_back(SolverType::DC);
             res.push_back(SolverType::GaussSeidelSynch);
             res.push_back(SolverType::SparseLUSingleSlack);
+            res.push_back(SolverType::FDPF_XB_SparseLU);
+            res.push_back(SolverType::FDPF_BX_SparseLU);
             #ifdef KLU_SOLVER_AVAILABLE
                 res.push_back(SolverType::KLU);
                 res.push_back(SolverType::KLUSingleSlack);
                 res.push_back(SolverType::KLUDC);
+                res.push_back(SolverType::FDPF_XB_KLU);
+                res.push_back(SolverType::FDPF_BX_KLU);
             #endif
             #ifdef NICSLU_SOLVER_AVAILABLE
                 res.push_back(SolverType::NICSLU);
                 res.push_back(SolverType::NICSLUSingleSlack);
                 res.push_back(SolverType::NICSLUDC);
+                res.push_back(SolverType::FDPF_XB_NICSLU);
+                res.push_back(SolverType::FDPF_BX_NICSLU);
             #endif
             #ifdef CKTSO_SOLVER_AVAILABLE
                 res.push_back(SolverType::CKTSO);
                 res.push_back(SolverType::CKTSOSingleSlack);
                 res.push_back(SolverType::CKTSODC);
+                res.push_back(SolverType::FDPF_XB_CKTSO);
+                res.push_back(SolverType::FDPF_BX_CKTSO);
             #endif
             return res;
         }
@@ -75,7 +88,42 @@ class ChooseSolver
             return res;
         }
         SolverType get_type() const {return _solver_type;}
-        
+
+        // TODO DO that for all solvers
+        FDPF_XB_SparseLUSolver & get_fdpf_xb_lu() { return _solver_fdpf_xb_lu;}
+        FDPF_BX_SparseLUSolver & get_fdpf_bx_lu() { return _solver_fdpf_bx_lu;}
+
+        void set_gridmodel(const GridModel * gridmodel){
+            _solver_lu.set_gridmodel(gridmodel);
+            _solver_lu_single.set_gridmodel(gridmodel);
+            _solver_gaussseidel.set_gridmodel(gridmodel);
+            _solver_gaussseidelsynch.set_gridmodel(gridmodel);
+            _solver_dc.set_gridmodel(gridmodel);
+            _solver_fdpf_xb_lu.set_gridmodel(gridmodel);
+            _solver_fdpf_bx_lu.set_gridmodel(gridmodel);
+            #ifdef KLU_SOLVER_AVAILABLE
+                _solver_klu.set_gridmodel(gridmodel);
+                _solver_klu_single.set_gridmodel(gridmodel);
+                _solver_klu_dc.set_gridmodel(gridmodel);
+                _solver_fdpf_xb_klu.set_gridmodel(gridmodel);
+                _solver_fdpf_bx_klu.set_gridmodel(gridmodel);
+            #endif  // KLU_SOLVER_AVAILABLE
+            #ifdef NICSLU_SOLVER_AVAILABLE
+                _solver_nicslu.set_gridmodel(gridmodel);
+                _solver_nicslu_single.set_gridmodel(gridmodel);
+                _solver_nicslu_dc.set_gridmodel(gridmodel);
+                _solver_fdpf_xb_nicslu.set_gridmodel(gridmodel);
+                _solver_fdpf_bx_nicslu.set_gridmodel(gridmodel);
+            #endif  // NICSLU_SOLVER_AVAILABLE
+            #ifdef CKTSO_SOLVER_AVAILABLE
+                _solver_cktso.set_gridmodel(gridmodel);
+                _solver_cktso_single.set_gridmodel(gridmodel);
+                _solver_cktso_dc.set_gridmodel(gridmodel);
+                _solver_fdpf_xb_cktso.set_gridmodel(gridmodel);
+                _solver_fdpf_bx_cktso.set_gridmodel(gridmodel);
+            #endif  // CKTSO_SOLVER_AVAILABLE
+        }
+
         bool ac_solver_used() const{
             auto p_solver = get_prt_solver("ac_solver_used", false);
             return p_solver->IS_AC;
@@ -86,7 +134,12 @@ class ChooseSolver
             if(type == _solver_type) return;
 
             #ifndef KLU_SOLVER_AVAILABLE
-                if((type == SolverType::KLU) || (type == SolverType::KLUDC) || (type == SolverType::KLUSingleSlack)){
+                if((type == SolverType::KLU) || 
+                   (type == SolverType::KLUDC) || 
+                   (type == SolverType::KLUSingleSlack) ||
+                   (type == SolverType::FDPF_XB_KLU) ||
+                   (type == SolverType::FDPF_BX_KLU)
+                   ){
                     std::string msg;
                     msg = "Impossible to change for a solver using KLU for linear algebra. Please compile lightsim2grid from source to benefit from this.";
                     throw std::runtime_error(msg);
@@ -94,7 +147,12 @@ class ChooseSolver
             #endif
 
             #ifndef NICSLU_SOLVER_AVAILABLE
-                if((type == SolverType::NICSLU) || (type == SolverType::NICSLUDC) || (type ==  SolverType::NICSLUSingleSlack)){
+                if((type == SolverType::NICSLU) || 
+                   (type == SolverType::NICSLUDC) || 
+                   (type ==  SolverType::NICSLUSingleSlack) || 
+                   (type ==  SolverType::FDPF_XB_NICSLU) ||
+                   (type ==  SolverType::FDPF_BX_NICSLU)
+                   ){
                     std::string msg;
                     msg = "Impossible to change for a solver using NICSLU for linear algebra. Please compile lightsim2grid from source to benefit from this.";
                     throw std::runtime_error(msg);
@@ -102,7 +160,12 @@ class ChooseSolver
             #endif
 
             #ifndef CKTSO_SOLVER_AVAILABLE
-                if((type == SolverType::CKTSO) || (type == SolverType::CKTSODC) || (type ==  SolverType::CKTSOSingleSlack)){
+                if((type == SolverType::CKTSO) || 
+                   (type == SolverType::CKTSODC) || 
+                   (type ==  SolverType::CKTSOSingleSlack) ||
+                   (type ==  SolverType::FDPF_XB_CKTSO) ||
+                   (type ==  SolverType::FDPF_BX_CKTSO)
+                   ){
                     std::string msg;
                     msg = "Impossible to change for a solver using CKTSO for linear algebra. Please compile lightsim2grid from source to benefit from this.";
                     throw std::runtime_error(msg);
@@ -178,16 +241,25 @@ class ChooseSolver
             else if(_solver_type == SolverType::CKTSOSingleSlack){
                 return _solver_cktso_single.get_J();}
             #endif // CKTSO_SOLVER_AVAILABLE
-            else if(_solver_type == SolverType::GaussSeidel){
-                throw std::runtime_error("ChooseSolver::get_J: There is not Jacobian matrix for the GaussSeidel powerflow.");}
             else if(_solver_type == SolverType::DC || 
                     _solver_type == SolverType::KLUDC || 
                     _solver_type == SolverType::NICSLUDC ||
                     _solver_type == SolverType::CKTSODC){
                 throw std::runtime_error("ChooseSolver::get_J: There is not Jacobian matrix for the DC powerflow.");}
+            else if(_solver_type == SolverType::FDPF_XB_SparseLU || 
+                    _solver_type == SolverType::FDPF_BX_SparseLU || 
+                    _solver_type == SolverType::FDPF_XB_KLU || 
+                    _solver_type == SolverType::FDPF_BX_KLU || 
+                    _solver_type == SolverType::FDPF_XB_NICSLU ||
+                    _solver_type == SolverType::FDPF_BX_NICSLU ||
+                    _solver_type == SolverType::FDPF_XB_CKTSO ||
+                    _solver_type == SolverType::FDPF_BX_CKTSO){
+                throw std::runtime_error("ChooseSolver::get_J: There is not Jacobian matrix for the FDPF powerflow.");}
             else if(_solver_type == SolverType::GaussSeidelSynch){
                 throw std::runtime_error("ChooseSolver::get_J: There is not Jacobian matrix for the GaussSeidelSynch powerflow.");}
-            else throw std::runtime_error("Unknown solver type encountered");
+            else if(_solver_type == SolverType::GaussSeidel){
+                throw std::runtime_error("ChooseSolver::get_J: There is not Jacobian matrix for the GaussSeidel powerflow.");}
+            else throw std::runtime_error("Unknown solver type encountered (get_J)");
         }
 
         /** apparently i cannot pass a const ref for a sparse matrix in python**/
@@ -224,39 +296,57 @@ class ChooseSolver
             
             #ifndef KLU_SOLVER_AVAILABLE
                 if(_solver_type == SolverType::KLU){
-                    std::string msg = "Impossible to use the KLU solver, it is not available on your platform.";
+                    std::string msg = "Impossible to use the KLU linear solver, your version of lightsim2grid has not been compiled to use it.";
                     throw std::runtime_error(msg);
                 } else if(_solver_type == SolverType::KLUSingleSlack){
-                    std::string msg = "Impossible to use the KLU solver, it is not available on your platform.";
+                    std::string msg = "Impossible to use the KLU linear solver, your version of lightsim2grid has not been compiled to use it.";
                     throw std::runtime_error(msg);
                 } else if(_solver_type == SolverType::KLUDC){
-                    std::string msg = "Impossible to use the KLU solver, it is not available on your platform.";
+                    std::string msg = "Impossible to use the KLU linear solver, your version of lightsim2grid has not been compiled to use it.";
+                    throw std::runtime_error(msg);
+                } else if(_solver_type == SolverType::FDPF_XB_KLU){
+                    std::string msg = "Impossible to use the KLU linear solver, your version of lightsim2grid has not been compiled to use it.";
+                    throw std::runtime_error(msg);
+                }else if(_solver_type == SolverType::FDPF_BX_KLU){
+                    std::string msg = "Impossible to use the KLU linear solver, your version of lightsim2grid has not been compiled to use it.";
                     throw std::runtime_error(msg);
                 }
             #endif  // KLU_SOLVER_AVAILABLE
 
             #ifndef NICSLU_SOLVER_AVAILABLE
                 if(_solver_type == SolverType::NICSLU){
-                    std::string msg = "Impossible to use the NICSLU solver, it is not available on your platform.";
+                    std::string msg = "Impossible to use the NICSLU linear solver, your version of lightsim2grid has not been compiled to use it.";
                     throw std::runtime_error(msg);
                 } else if(_solver_type == SolverType::NICSLUSingleSlack){
-                    std::string msg = "Impossible to use the NICSLU solver, it is not available on your platform.";
+                    std::string msg = "Impossible to use the NICSLU linear solver, your version of lightsim2grid has not been compiled to use it.";
                     throw std::runtime_error(msg);
                 } else if(_solver_type == SolverType::NICSLUDC){
-                    std::string msg = "Impossible to use the NICSLU solver, it is not available on your platform.";
+                    std::string msg = "Impossible to use the NICSLU linear solver, your version of lightsim2grid has not been compiled to use it.";
+                    throw std::runtime_error(msg);
+                } else if(_solver_type == SolverType::FDPF_XB_NICSLU){
+                    std::string msg = "Impossible to use the NICSLU linear solver, your version of lightsim2grid has not been compiled to use it.";
+                    throw std::runtime_error(msg);
+                } else if(_solver_type == SolverType::FDPF_BX_NICSLU){
+                    std::string msg = "Impossible to use the NICSLU linear solver, your version of lightsim2grid has not been compiled to use it.";
                     throw std::runtime_error(msg);
                 }
             #endif  // NICSLU_SOLVER_AVAILABLE
 
             #ifndef CKTSO_SOLVER_AVAILABLE
                 if(_solver_type == SolverType::CKTSO){
-                    std::string msg = "Impossible to use the CKTSO solver, it is not available on your platform.";
+                    std::string msg = "Impossible to use the CKTSO linear solver, your version of lightsim2grid has not been compiled to use it.";
                     throw std::runtime_error(msg);
                 } else if(_solver_type == SolverType::CKTSOSingleSlack){
-                    std::string msg = "Impossible to use the CKTSO solver, it is not available on your platform.";
+                    std::string msg = "Impossible to use the CKTSO linear solver, your version of lightsim2grid has not been compiled to use it.";
                     throw std::runtime_error(msg);
                 } else if(_solver_type == SolverType::CKTSODC){
-                    std::string msg = "Impossible to use the CKTSO solver, it is not available on your platform.";
+                    std::string msg = "Impossible to use the CKTSO linear solver, your version of lightsim2grid has not been compiled to use it.";
+                    throw std::runtime_error(msg);
+                } else if(_solver_type == SolverType::FDPF_XB_CKTSO){
+                    std::string msg = "Impossible to use the CKTSO linear solver, your version of lightsim2grid has not been compiled to use it.";
+                    throw std::runtime_error(msg);
+                } else if(_solver_type == SolverType::FDPF_BX_CKTSO){
+                    std::string msg = "Impossible to use the CKTSO linear solver, your version of lightsim2grid has not been compiled to use it.";
                     throw std::runtime_error(msg);
                 }
             #endif  // CKTSO_SOLVER_AVAILABLE
@@ -272,24 +362,32 @@ class ChooseSolver
             if(_solver_type == SolverType::SparseLU){res = &_solver_lu;}
             else if(_solver_type == SolverType::SparseLUSingleSlack){res = &_solver_lu_single;}
             else if(_solver_type == SolverType::DC){res = &_solver_dc;}
+            else if(_solver_type == SolverType::FDPF_XB_SparseLU){res = &_solver_fdpf_xb_lu;}
+            else if(_solver_type == SolverType::FDPF_BX_SparseLU){res = &_solver_fdpf_bx_lu;}
             #ifdef KLU_SOLVER_AVAILABLE
             else if(_solver_type == SolverType::KLU){res = & _solver_klu;}
             else if(_solver_type == SolverType::KLUSingleSlack){res = &_solver_klu_single;}
             else if(_solver_type == SolverType::KLUDC){res = &_solver_klu_dc;}
+            else if(_solver_type == SolverType::FDPF_XB_KLU){res = &_solver_fdpf_xb_klu;}
+            else if(_solver_type == SolverType::FDPF_BX_KLU){res = &_solver_fdpf_bx_klu;}
             #endif  // KLU_SOLVER_AVAILABLE
             #ifdef NICSLU_SOLVER_AVAILABLE
             else if(_solver_type == SolverType::NICSLU){res = &_solver_nicslu;}
             else if(_solver_type == SolverType::NICSLUSingleSlack){res = &_solver_nicslu_single;}
             else if(_solver_type == SolverType::NICSLUDC){res = &_solver_nicslu_dc;}
+            else if(_solver_type == SolverType::FDPF_XB_NICSLU){res = &_solver_fdpf_xb_nicslu;}
+            else if(_solver_type == SolverType::FDPF_BX_NICSLU){res = &_solver_fdpf_bx_nicslu;}
             #endif // NICSLU_SOLVER_AVAILABLE
             #ifdef CKTSO_SOLVER_AVAILABLE
             else if(_solver_type == SolverType::CKTSO){res = &_solver_cktso;}
             else if(_solver_type == SolverType::CKTSOSingleSlack){res = &_solver_cktso_single;}
             else if(_solver_type == SolverType::CKTSODC){res = &_solver_cktso_dc;}
+            else if(_solver_type == SolverType::FDPF_XB_CKTSO){res = &_solver_fdpf_xb_cktso;}
+            else if(_solver_type == SolverType::FDPF_BX_CKTSO){res = &_solver_fdpf_bx_cktso;}
             #endif // CKTSO_SOLVER_AVAILABLE
             else if(_solver_type == SolverType::GaussSeidel){res = &_solver_gaussseidel;}
             else if(_solver_type == SolverType::GaussSeidelSynch){res = &_solver_gaussseidelsynch;}
-            else throw std::runtime_error("Unknown solver type encountered");
+            else throw std::runtime_error("Unknown solver type encountered (ChooseSolver get_prt_solver const)");
             return res;
         }
         BaseSolver * get_prt_solver(const std::string & error_msg, bool check_right_solver_=true) {
@@ -298,24 +396,32 @@ class ChooseSolver
             if(_solver_type == SolverType::SparseLU){res = &_solver_lu;}
             else if(_solver_type == SolverType::SparseLUSingleSlack){res = &_solver_lu_single;}
             else if(_solver_type == SolverType::DC){res = &_solver_dc;}
+            else if(_solver_type == SolverType::FDPF_XB_SparseLU){res = &_solver_fdpf_xb_lu;}
+            else if(_solver_type == SolverType::FDPF_BX_SparseLU){res = &_solver_fdpf_bx_lu;}
             #ifdef KLU_SOLVER_AVAILABLE
             else if(_solver_type == SolverType::KLU){res = & _solver_klu;}
             else if(_solver_type == SolverType::KLUSingleSlack){res = &_solver_klu_single;}
             else if(_solver_type == SolverType::KLUDC){res = &_solver_klu_dc;}
+            else if(_solver_type == SolverType::FDPF_XB_KLU){res = &_solver_fdpf_xb_klu;}
+            else if(_solver_type == SolverType::FDPF_BX_KLU){res = &_solver_fdpf_bx_klu;}
             #endif  // KLU_SOLVER_AVAILABLE
             #ifdef NICSLU_SOLVER_AVAILABLE
             else if(_solver_type == SolverType::NICSLU){res = &_solver_nicslu;}
             else if(_solver_type == SolverType::NICSLUSingleSlack){res = &_solver_nicslu_single;}
             else if(_solver_type == SolverType::NICSLUDC){res = &_solver_nicslu_dc;}
+            else if(_solver_type == SolverType::FDPF_XB_NICSLU){res = &_solver_fdpf_xb_nicslu;}
+            else if(_solver_type == SolverType::FDPF_BX_NICSLU){res = &_solver_fdpf_bx_nicslu;}
             #endif // NICSLU_SOLVER_AVAILABLE
             #ifdef CKTSO_SOLVER_AVAILABLE
             else if(_solver_type == SolverType::CKTSO){res = &_solver_cktso;}
             else if(_solver_type == SolverType::CKTSOSingleSlack){res = &_solver_cktso_single;}
             else if(_solver_type == SolverType::CKTSODC){res = &_solver_cktso_dc;}
+            else if(_solver_type == SolverType::FDPF_XB_CKTSO){res = &_solver_fdpf_xb_cktso;}
+            else if(_solver_type == SolverType::FDPF_BX_CKTSO){res = &_solver_fdpf_bx_cktso;}
             #endif // CKTSO_SOLVER_AVAILABLE
             else if(_solver_type == SolverType::GaussSeidel){res = &_solver_gaussseidel;}
             else if(_solver_type == SolverType::GaussSeidelSynch){res = &_solver_gaussseidelsynch;}
-            else throw std::runtime_error("Unknown solver type encountered");
+            else throw std::runtime_error("Unknown solver type encountered (ChooseSolver get_prt_solver non const)");
             return res;
         }
     protected:
@@ -329,20 +435,28 @@ class ChooseSolver
         GaussSeidelSolver _solver_gaussseidel;
         GaussSeidelSynchSolver _solver_gaussseidelsynch;
         DCSolver _solver_dc;
+        FDPF_XB_SparseLUSolver _solver_fdpf_xb_lu;
+        FDPF_BX_SparseLUSolver _solver_fdpf_bx_lu;
         #ifdef KLU_SOLVER_AVAILABLE
             KLUSolver _solver_klu;
             KLUSolverSingleSlack _solver_klu_single;
             KLUDCSolver _solver_klu_dc;
+            FDPF_XB_KLUSolver _solver_fdpf_xb_klu;
+            FDPF_BX_KLUSolver _solver_fdpf_bx_klu;
         #endif  // KLU_SOLVER_AVAILABLE
         #ifdef NICSLU_SOLVER_AVAILABLE
             NICSLUSolver _solver_nicslu;
             NICSLUSolverSingleSlack _solver_nicslu_single;
             NICSLUDCSolver _solver_nicslu_dc;
+            FDPF_XB_NICSLUSolver _solver_fdpf_xb_nicslu;
+            FDPF_BX_NICSLUSolver _solver_fdpf_bx_nicslu;
         #endif  // NICSLU_SOLVER_AVAILABLE
         #ifdef CKTSO_SOLVER_AVAILABLE
             CKTSOSolver _solver_cktso;
             CKTSOSolverSingleSlack _solver_cktso_single;
             CKTSODCSolver _solver_cktso_dc;
+            FDPF_XB_CKTSOSolver _solver_fdpf_xb_cktso;
+            FDPF_BX_CKTSOSolver _solver_fdpf_bx_cktso;
         #endif  // CKTSO_SOLVER_AVAILABLE
 };
 
diff --git a/src/DCSolver.tpp b/src/DCSolver.tpp
index c46abcc..ff38cc1 100644
--- a/src/DCSolver.tpp
+++ b/src/DCSolver.tpp
@@ -27,7 +27,7 @@ bool BaseDCSolver<LinearSolver>::compute_pf(const Eigen::SparseMatrix<cplx_type>
     //   and for the slack bus both the magnitude and the angle are used.
 
     auto timer = CustTimer();
-
+    BaseSolver::reset_timer();
     const int nb_bus_solver = static_cast<int>(Ybus.rows());
 
     #ifdef __COUT_TIMES
@@ -40,25 +40,35 @@ bool BaseDCSolver<LinearSolver>::compute_pf(const Eigen::SparseMatrix<cplx_type>
     // TODO SLACK (for now i put all slacks as PV, except the first one)
     // this should be handled in Sbus, because we know the amount of power absorbed by the slack
     // so we can compute it correctly !
-    Eigen::VectorXi my_pv = retrieve_pv_with_slack(slack_ids, pv);
-
-    // find the slack bus
-    int slack_bus_id_solver = extract_slack_bus_id(my_pv, pq, nb_bus_solver);
-
+    const Eigen::VectorXi my_pv = retrieve_pv_with_slack(slack_ids, pv);
+    // const Eigen::VectorXi & my_pv = pv;
+
+    // find the slack buses
+    Eigen::VectorXi slack_bus_ids_solver = extract_slack_bus_id(my_pv, pq, nb_bus_solver);
+    // corresp bus -> solverbus
+    Eigen::VectorXi ybus_to_me = Eigen::VectorXi::Constant(nb_bus_solver, -1);
+    // Eigen::VectorXi me_to_ybus = Eigen::VectorXi::Constant(nb_bus_solver - slack_bus_ids_solver.size(), -1);
+    int solver_id = 0;
+    for (int ybus_id=0; ybus_id < nb_bus_solver; ++ybus_id){
+        if(isin(ybus_id, slack_bus_ids_solver)) continue;  // I don't add anything to the slack bus
+        ybus_to_me(ybus_id) = solver_id;
+        // me_to_ybus(solver_id) = ybus_id;
+        ++solver_id;
+    }
     // remove the slack bus from Ybus
     // and extract only real part
     // TODO see if "prune" might work here https://eigen.tuxfamily.org/dox/classEigen_1_1SparseMatrix.html#title29
     std::vector<Eigen::Triplet<real_type> > tripletList;
     tripletList.reserve(Ybus.nonZeros());
     for (int k=0; k < nb_bus_solver; ++k){
-        if(k == slack_bus_id_solver) continue;  // I don't add anything to the slack bus
+        if(ybus_to_me(k) == -1) continue;  // I don't add anything to the slack bus
         for (Eigen::SparseMatrix<cplx_type>::InnerIterator it(Ybus, k); it; ++it)
         {
-            int row_res = static_cast<int>(it.row());
-            if(row_res == slack_bus_id_solver) continue;
-            row_res = row_res > slack_bus_id_solver ? row_res - 1 : row_res;
-            int col_res = static_cast<int>(it.col());
-            col_res = col_res > slack_bus_id_solver ? col_res - 1 : col_res;
+            int row_res = static_cast<int>(it.row());  // TODO Eigen::Index here ?
+            if(ybus_to_me(row_res) == -1) continue;
+            row_res = ybus_to_me(row_res);
+            int col_res = static_cast<int>(it.col());  // should be k   // TODO Eigen::Index here ?
+            col_res = ybus_to_me(col_res);
             tripletList.push_back(Eigen::Triplet<real_type> (row_res, col_res, std::real(it.value())));
         }
     }
@@ -85,11 +95,10 @@ bool BaseDCSolver<LinearSolver>::compute_pf(const Eigen::SparseMatrix<cplx_type>
     }
 
     // remove the slack bus from Sbus
-    RealVect dcSbus = RealVect::Constant(nb_bus_solver - 1, my_zero_);
+    RealVect dcSbus = RealVect::Constant(nb_bus_solver - slack_bus_ids_solver.size(), my_zero_);
     for (int k=0; k < nb_bus_solver; ++k){
-        if(k == slack_bus_id_solver) continue;  // I don't add anything to the slack bus
-        int col_res = k;
-        col_res = col_res > slack_bus_id_solver ? col_res - 1 : col_res;
+        if(ybus_to_me(k) == -1) continue;  // I don't add anything to the slack bus
+        const int col_res = ybus_to_me(k);
         dcSbus(col_res) = std::real(Sbus_tmp(k));
     }
 
@@ -124,19 +133,19 @@ bool BaseDCSolver<LinearSolver>::compute_pf(const Eigen::SparseMatrix<cplx_type>
     // see the BaseSolver.cpp: _evaluate_Fx
     RealVect Va_dc = RealVect::Constant(nb_bus_solver, my_zero_);
     // fill Va from dc approx
-    for (int bus_id_with_slack=0; bus_id_with_slack < nb_bus_solver; ++bus_id_with_slack){
-        if(bus_id_with_slack == slack_bus_id_solver) continue;  // slack bus is handled elsewhere
-        int bus_id_without_slack = bus_id_with_slack > slack_bus_id_solver ? bus_id_with_slack - 1 : bus_id_with_slack;
-        Va_dc(bus_id_with_slack) = Va_dc_without_slack(bus_id_without_slack);
+    for (int ybus_id=0; ybus_id < nb_bus_solver; ++ybus_id){
+        if(ybus_to_me(ybus_id) == -1) continue;  // slack bus is handled elsewhere
+        const int bus_me = ybus_to_me(ybus_id);
+        Va_dc(ybus_id) = Va_dc_without_slack(bus_me);
     }
-    Va_dc.array() += std::arg(V(slack_bus_id_solver));
+    Va_dc.array() += std::arg(V(slack_bus_ids_solver(0)));
 
     // save the results
     Va_ = Va_dc;
 
     // add the Voltage setpoints of the generator
     Vm_ = V.array().abs();
-    Vm_(slack_bus_id_solver) = std::abs(V(slack_bus_id_solver));
+    Vm_(slack_bus_ids_solver) = V(slack_bus_ids_solver).array().abs();
 
     // now compute the resulting complex voltage
     V_ = (Va_.array().cos().template cast<cplx_type>() + my_i * Va_.array().sin().template cast<cplx_type>());
diff --git a/src/DataDCLine.h b/src/DataDCLine.h
index 7a84470..185be25 100644
--- a/src/DataDCLine.h
+++ b/src/DataDCLine.h
@@ -198,10 +198,10 @@ class DataDCLine : public DataGeneric
     }
 
     // solver stuff
-    virtual void fillSbus(CplxVect & Sbus, const std::vector<int> & id_grid_to_solver){
-        from_gen_.fillSbus(Sbus, id_grid_to_solver);   
-        to_gen_.fillSbus(Sbus, id_grid_to_solver);   
-    }
+    virtual void fillSbus(CplxVect & Sbus, const std::vector<int> & id_grid_to_solver, bool ac) const{
+        from_gen_.fillSbus(Sbus, id_grid_to_solver, ac);   
+        to_gen_.fillSbus(Sbus, id_grid_to_solver, ac);   
+    } 
 
     virtual void fillpv(std::vector<int>& bus_pv,
                         std::vector<bool> & has_bus_been_added,
diff --git a/src/DataGen.cpp b/src/DataGen.cpp
index 5f3d893..238a658 100644
--- a/src/DataGen.cpp
+++ b/src/DataGen.cpp
@@ -115,7 +115,7 @@ RealVect DataGen::get_slack_weights(Eigen::Index nb_bus_solver, const std::vecto
     return res;
 }
 
-void DataGen::fillSbus(CplxVect & Sbus, const std::vector<int> & id_grid_to_solver){
+void DataGen::fillSbus(CplxVect & Sbus, const std::vector<int> & id_grid_to_solver, bool ac) const {
     const int nb_gen = nb();
     int bus_id_me, bus_id_solver;
     real_type tmp;
diff --git a/src/DataGen.h b/src/DataGen.h
index 99f2ae9..237f118 100644
--- a/src/DataGen.h
+++ b/src/DataGen.h
@@ -183,7 +183,7 @@ class DataGen: public DataGeneric
     void change_p(int gen_id, real_type new_p, bool & need_reset);
     void change_v(int gen_id, real_type new_v_pu, bool & need_reset);
 
-    virtual void fillSbus(CplxVect & Sbus, const std::vector<int> & id_grid_to_solver);
+    virtual void fillSbus(CplxVect & Sbus, const std::vector<int> & id_grid_to_solver, bool ac) const;
     virtual void fillpv(std::vector<int>& bus_pv,
                         std::vector<bool> & has_bus_been_added,
                         Eigen::VectorXi & slack_bus_id_solver,
diff --git a/src/DataGeneric.h b/src/DataGeneric.h
index a2be43f..8d9bc60 100644
--- a/src/DataGeneric.h
+++ b/src/DataGeneric.h
@@ -70,9 +70,16 @@ class DataGeneric : public BaseConstants
         virtual void fillYbus(std::vector<Eigen::Triplet<cplx_type> > & res,
                               bool ac,
                               const std::vector<int> & id_grid_to_solver,
-                              real_type sn_mva) {};
+                              real_type sn_mva) const {};
+
+        virtual void fillBp_Bpp(std::vector<Eigen::Triplet<real_type> > & Bp,
+                                std::vector<Eigen::Triplet<real_type> > & Bpp,
+                                const std::vector<int> & id_grid_to_solver,
+                                real_type sn_mva,
+                                FDPFMethod xb_or_bx) const {};
+
         virtual void fillYbus(Eigen::SparseMatrix<cplx_type> & res, bool ac, const std::vector<int> & id_grid_to_solver) {};
-        virtual void fillSbus(CplxVect & Sbus, const std::vector<int> & id_grid_to_solver){};
+        virtual void fillSbus(CplxVect & Sbus, const std::vector<int> & id_grid_to_solver, bool ac) const {};
         virtual void fillpv(std::vector<int>& bus_pv,
                             std::vector<bool> & has_bus_been_added,
                             Eigen::VectorXi & slack_bus_id_solver,
diff --git a/src/DataLine.cpp b/src/DataLine.cpp
index ac5b4c1..ecf4f59 100644
--- a/src/DataLine.cpp
+++ b/src/DataLine.cpp
@@ -30,7 +30,38 @@ void DataLine::init(const RealVect & branch_r,
 
     bus_or_id_ = branch_from_id;
     bus_ex_id_ = branch_to_id;
-    powerlines_h_ = branch_h;
+    powerlines_h_or_ = 0.5 * branch_h;
+    powerlines_h_ex_ = 0.5 * branch_h;
+    powerlines_r_ = branch_r;
+    powerlines_x_ = branch_x;
+    status_ = std::vector<bool>(branch_r.size(), true); // by default everything is connected
+    _update_model_coeffs();
+}
+
+void DataLine::init(const RealVect & branch_r,
+                    const RealVect & branch_x,
+                    const CplxVect & branch_h_or,
+                    const CplxVect & branch_h_ex,
+                    const Eigen::VectorXi & branch_from_id,
+                    const Eigen::VectorXi & branch_to_id
+                    )
+{
+    /**
+    This method initialize the Ybus matrix from the branch matrix.
+    It has to be called once when the solver is initialized. Afterwards, a call to
+    "updateYbus" should be made for performance optimiaztion instead. //TODO
+    **/
+
+    // TODO check what can be checked: branch_* have same size, no id in branch_to_id that are
+    // TODO not in [0, .., buv_vn_kv.size()] etc.
+
+    //TODO consistency with trafo: have a converter methods to convert this value into pu, and store the pu
+    // in this method
+
+    bus_or_id_ = branch_from_id;
+    bus_ex_id_ = branch_to_id;
+    powerlines_h_or_ = branch_h_or;
+    powerlines_h_ex_ = branch_h_ex;
     powerlines_r_ = branch_r;
     powerlines_x_ = branch_x;
     status_ = std::vector<bool>(branch_r.size(), true); // by default everything is connected
@@ -41,11 +72,12 @@ DataLine::StateRes DataLine::get_state() const
 {
      std::vector<real_type> branch_r(powerlines_r_.begin(), powerlines_r_.end());
      std::vector<real_type> branch_x(powerlines_x_.begin(), powerlines_x_.end());
-     std::vector<cplx_type > branch_h(powerlines_h_.begin(), powerlines_h_.end());
+     std::vector<cplx_type > branch_hor(powerlines_h_or_.begin(), powerlines_h_or_.end());
+     std::vector<cplx_type > branch_hex(powerlines_h_ex_.begin(), powerlines_h_ex_.end());
      std::vector<int > branch_from_id(bus_or_id_.begin(), bus_or_id_.end());
      std::vector<int > branch_to_id(bus_ex_id_.begin(), bus_ex_id_.end());
      std::vector<bool> status = status_;
-     DataLine::StateRes res(branch_r, branch_x, branch_h, branch_from_id, branch_to_id, status);
+     DataLine::StateRes res(branch_r, branch_x, branch_hor, branch_hex, branch_from_id, branch_to_id, status);
      return res;
 }
 void DataLine::set_state(DataLine::StateRes & my_state)
@@ -54,16 +86,18 @@ void DataLine::set_state(DataLine::StateRes & my_state)
 
     std::vector<real_type> & branch_r = std::get<0>(my_state);
     std::vector<real_type> & branch_x = std::get<1>(my_state);
-    std::vector<cplx_type > & branch_h = std::get<2>(my_state);
-    std::vector<int> & branch_from_id = std::get<3>(my_state);
-    std::vector<int> & branch_to_id = std::get<4>(my_state);
-    std::vector<bool> & status = std::get<5>(my_state);
+    std::vector<cplx_type > & branch_h_or = std::get<2>(my_state);
+    std::vector<cplx_type > & branch_h_ex = std::get<3>(my_state);
+    std::vector<int> & branch_from_id = std::get<4>(my_state);
+    std::vector<int> & branch_to_id = std::get<5>(my_state);
+    std::vector<bool> & status = std::get<6>(my_state);
     // TODO check sizes
 
     // now assign the values
     powerlines_r_ = RealVect::Map(&branch_r[0], branch_r.size());
     powerlines_x_ = RealVect::Map(&branch_x[0], branch_x.size());
-    powerlines_h_ = CplxVect::Map(&branch_h[0], branch_h.size());
+    powerlines_h_or_ = CplxVect::Map(&branch_h_or[0], branch_h_or.size());
+    powerlines_h_ex_ = CplxVect::Map(&branch_h_ex[0], branch_h_ex.size());
 
     // input data
     bus_or_id_ = Eigen::VectorXi::Map(&branch_from_id[0], branch_from_id.size());
@@ -91,9 +125,10 @@ void DataLine::_update_model_coeffs()
         // for AC
         // see https://matpower.org/docs/MATPOWER-manual.pdf eq. 3.2
         const cplx_type ys = 1. / (powerlines_r_(i) + my_i * powerlines_x_(i));
-        const cplx_type h = my_i * powerlines_h_(i) * 0.5;
-        yac_ff_(i) = (ys + h);
-        yac_tt_(i) = (ys + h);
+        const cplx_type h_or = my_i * powerlines_h_or_(i);
+        const cplx_type h_ex = my_i * powerlines_h_ex_(i);
+        yac_ff_(i) = (ys + h_or);
+        yac_tt_(i) = (ys + h_ex);
         yac_tf_(i) = -ys;
         yac_ft_(i) = -ys;
 
@@ -116,7 +151,7 @@ void DataLine::fillYbus_spmat(Eigen::SparseMatrix<cplx_type> & res, bool ac, con
 void DataLine::fillYbus(std::vector<Eigen::Triplet<cplx_type> > & res,
                         bool ac,
                         const std::vector<int> & id_grid_to_solver,
-                        real_type sn_mva)
+                        real_type sn_mva) const
 {
     // fill the matrix
     //TODO template here instead of "if" for ac / dc
@@ -169,6 +204,90 @@ void DataLine::fillYbus(std::vector<Eigen::Triplet<cplx_type> > & res,
     }
 }
 
+void DataLine::fillBp_Bpp(std::vector<Eigen::Triplet<real_type> > & Bp,
+                          std::vector<Eigen::Triplet<real_type> > & Bpp,
+                          const std::vector<int> & id_grid_to_solver,
+                          real_type sn_mva,
+                          FDPFMethod xb_or_bx) const
+{
+
+    // For Bp
+    // temp_branch[:, BR_B] = zeros(nl)           ## zero out line charging shunts
+    // temp_branch[:, TAP] = ones(nl)             ## cancel out taps
+    // if alg == 2:                               ## if XB method
+    //    temp_branch[:, BR_R] = zeros(nl)       ## zero out line resistance
+
+    // For Bpp
+    // temp_branch[:, SHIFT] = zeros(nl)          ## zero out phase shifters
+    // if alg == 3:                               ## if BX method
+    //     temp_branch[:, BR_R] = zeros(nl)    ## zero out line resistance
+    const Eigen::Index nb_line = static_cast<int>(powerlines_r_.size());
+    real_type yft_bp, ytf_bp, yff_bp, ytt_bp;
+    real_type yft_bpp, ytf_bpp, yff_bpp, ytt_bpp;
+    //diagonal coefficients
+    for(Eigen::Index line_id=0; line_id < nb_line; ++line_id){
+        // i only add this if the powerline is connected
+        if(!status_[line_id]) continue;
+
+        // get the from / to bus id
+        int bus_or_id_me = bus_or_id_(line_id);
+        int bus_or_solver_id = id_grid_to_solver[bus_or_id_me];
+        if(bus_or_solver_id == _deactivated_bus_id){
+            std::ostringstream exc_;
+            exc_ << "DataLine::fillBp_Bpp: the line with id ";
+            exc_ << line_id;
+            exc_ << " is connected (or side) to a disconnected bus while being connected";
+            throw std::runtime_error(exc_.str());
+        }
+        int bus_ex_id_me = bus_ex_id_(line_id);
+        int bus_ex_solver_id = id_grid_to_solver[bus_ex_id_me];
+        if(bus_ex_solver_id == _deactivated_bus_id){
+            std::ostringstream exc_;
+            exc_ << "DataLine::fillBp_Bpp: the line with id ";
+            exc_ << line_id;
+            exc_ << " is connected (ex side) to a disconnected bus while being connected";
+            throw std::runtime_error(exc_.str());
+        }
+
+        // get the coefficients
+        cplx_type ys_bp, ys_bpp;
+        if(xb_or_bx==FDPFMethod::XB){
+            ys_bp = 1. / (0. + my_i * powerlines_x_(line_id));
+            ys_bpp = 1. / (powerlines_r_(line_id) + my_i * powerlines_x_(line_id));
+        }else if (xb_or_bx==FDPFMethod::BX){
+            ys_bp = 1. / (powerlines_r_(line_id) + my_i * powerlines_x_(line_id));
+            ys_bpp = 1. / (0. + my_i * powerlines_x_(line_id));
+        }else{
+            std::ostringstream exc_;
+            exc_ << "DataLine::fillBp_Bpp: unknown method for the FDPF powerflow for line id ";
+            exc_ << line_id;
+            throw std::runtime_error(exc_.str());            
+        }
+        const real_type ys_bp_r = std::imag(ys_bp); 
+        yff_bp = ys_bp_r;
+        ytt_bp = ys_bp_r;
+        yft_bp = -ys_bp_r;
+        ytf_bp = -ys_bp_r;
+        const real_type ys_bpp_r = std::imag(ys_bpp); 
+        yff_bpp = ys_bpp_r + std::imag(my_i * powerlines_h_or_(line_id));
+        ytt_bpp = ys_bpp_r + std::imag(my_i * powerlines_h_ex_(line_id));
+        yft_bpp = -ys_bpp_r;
+        ytf_bpp = -ys_bpp_r;
+
+        // and now add them
+        Bp.push_back(Eigen::Triplet<real_type> (bus_or_solver_id, bus_ex_solver_id, -yft_bp));
+        Bp.push_back(Eigen::Triplet<real_type> (bus_ex_solver_id, bus_or_solver_id, -ytf_bp));
+        Bp.push_back(Eigen::Triplet<real_type> (bus_or_solver_id, bus_or_solver_id, -yff_bp));
+        Bp.push_back(Eigen::Triplet<real_type> (bus_ex_solver_id, bus_ex_solver_id, -ytt_bp));
+
+        Bpp.push_back(Eigen::Triplet<real_type> (bus_or_solver_id, bus_ex_solver_id, -yft_bpp));
+        Bpp.push_back(Eigen::Triplet<real_type> (bus_ex_solver_id, bus_or_solver_id, -ytf_bpp));
+        Bpp.push_back(Eigen::Triplet<real_type> (bus_or_solver_id, bus_or_solver_id, -yff_bpp));
+        Bpp.push_back(Eigen::Triplet<real_type> (bus_ex_solver_id, bus_ex_solver_id, -ytt_bpp));
+
+    }
+}
+
 void DataLine::reset_results()
 {
     res_powerline_por_ = RealVect();  // in MW
@@ -181,7 +300,6 @@ void DataLine::reset_results()
     res_powerline_aex_ = RealVect();  // in kA
 }
 
-
 void DataLine::compute_results(const Eigen::Ref<const RealVect> & Va,
                                const Eigen::Ref<const RealVect> & Vm,
                                const Eigen::Ref<const CplxVect> & V,
diff --git a/src/DataLine.h b/src/DataLine.h
index 317ceee..dae3e8a 100644
--- a/src/DataLine.h
+++ b/src/DataLine.h
@@ -44,6 +44,8 @@ class DataLine : public DataGeneric
                 real_type r_pu;
                 real_type x_pu;
                 cplx_type h_pu;
+                cplx_type h_or_pu;
+                cplx_type h_ex_pu;
 
                 bool has_res;
                 real_type res_p_or_mw;
@@ -58,13 +60,15 @@ class DataLine : public DataGeneric
                 real_type res_theta_ex_deg;
 
                 LineInfo(const DataLine & r_data_line, int my_id):
-                id(-1),
+                id(my_id),
                 connected(false),
                 bus_or_id(-1),
                 bus_ex_id(-1),
                 r_pu(-1.0),
                 x_pu(-1.0),
                 h_pu(0., 0.),
+                h_or_pu(0., 0.),
+                h_ex_pu(0., 0.),
                 has_res(false),
                 res_p_or_mw(0.),
                 res_q_or_mvar(0.),
@@ -85,7 +89,9 @@ class DataLine : public DataGeneric
                         bus_ex_id = r_data_line.bus_ex_id_.coeff(my_id);
                         r_pu = r_data_line.powerlines_r_.coeff(my_id);
                         x_pu = r_data_line.powerlines_x_.coeff(my_id);
-                        h_pu = r_data_line.powerlines_h_.coeff(my_id);
+                        h_or_pu = r_data_line.powerlines_h_or_.coeff(my_id);
+                        h_ex_pu = r_data_line.powerlines_h_ex_.coeff(my_id);
+                        h_pu = (h_or_pu + h_ex_pu);
 
                         has_res = r_data_line.res_powerline_por_.size() > 0;
                         if(has_res)
@@ -114,6 +120,7 @@ class DataLine : public DataGeneric
                std::vector<real_type>, // branch_r
                std::vector<real_type>, // branch_x
                std::vector<cplx_type>, // branch_h
+               std::vector<cplx_type>, // branch_h
                std::vector<int>, // branch_from_id
                std::vector<int>, // branch_to_id
                std::vector<bool> // status_
@@ -128,6 +135,14 @@ class DataLine : public DataGeneric
               const Eigen::VectorXi & branch_to_id
               );
 
+    void init(const RealVect & branch_r,
+              const RealVect & branch_x,
+              const CplxVect & branch_h_or,
+              const CplxVect & branch_h_ex,
+              const Eigen::VectorXi & branch_from_id,
+              const Eigen::VectorXi & branch_to_id
+              );
+
     // pickle
     DataLine::StateRes get_state() const;
     void set_state(DataLine::StateRes & my_state );
@@ -173,7 +188,12 @@ class DataLine : public DataGeneric
                           bool ac,
                           const std::vector<int> & id_grid_to_solver,
                           real_type sn_mva
-                          );
+                          ) const;
+    virtual void fillBp_Bpp(std::vector<Eigen::Triplet<real_type> > & Bp,
+                            std::vector<Eigen::Triplet<real_type> > & Bpp,
+                            const std::vector<int> & id_grid_to_solver,
+                            real_type sn_mva,
+                            FDPFMethod xb_or_bx) const;
     virtual void fillYbus_spmat(Eigen::SparseMatrix<cplx_type> & res, bool ac, const std::vector<int> & id_grid_to_solver);
 
     void compute_results(const Eigen::Ref<const RealVect> & Va,
@@ -214,7 +234,8 @@ class DataLine : public DataGeneric
         // physical properties
         RealVect powerlines_r_;
         RealVect powerlines_x_;
-        CplxVect powerlines_h_;
+        CplxVect powerlines_h_or_;
+        CplxVect powerlines_h_ex_;
 
         // input data
         Eigen::VectorXi bus_or_id_;
diff --git a/src/DataLoad.cpp b/src/DataLoad.cpp
index 70e8f6a..003ca81 100644
--- a/src/DataLoad.cpp
+++ b/src/DataLoad.cpp
@@ -47,7 +47,7 @@ void DataLoad::set_state(DataLoad::StateRes & my_state )
 }
 
 
-void DataLoad::fillSbus(CplxVect & Sbus, const std::vector<int> & id_grid_to_solver){
+void DataLoad::fillSbus(CplxVect & Sbus, const std::vector<int> & id_grid_to_solver, bool ac) const {
     int nb_load = nb();
     int bus_id_me, bus_id_solver;
     cplx_type tmp;
diff --git a/src/DataLoad.h b/src/DataLoad.h
index 89fe216..78effd6 100644
--- a/src/DataLoad.h
+++ b/src/DataLoad.h
@@ -139,7 +139,7 @@ class DataLoad : public DataGeneric
     void change_p(int load_id, real_type new_p, bool & need_reset);
     void change_q(int load_id, real_type new_q, bool & need_reset);
 
-    virtual void fillSbus(CplxVect & Sbus, const std::vector<int> & id_grid_to_solver);
+    virtual void fillSbus(CplxVect & Sbus, const std::vector<int> & id_grid_to_solver, bool ac) const;
 
     void compute_results(const Eigen::Ref<const RealVect> & Va,
                          const Eigen::Ref<const RealVect> & Vm,
diff --git a/src/DataSGen.cpp b/src/DataSGen.cpp
index d6ebbd9..c0c6fec 100644
--- a/src/DataSGen.cpp
+++ b/src/DataSGen.cpp
@@ -82,7 +82,7 @@ void DataSGen::set_state(DataSGen::StateRes & my_state )
 }
 
 
-void DataSGen::fillSbus(CplxVect & Sbus, const std::vector<int> & id_grid_to_solver){
+void DataSGen::fillSbus(CplxVect & Sbus, const std::vector<int> & id_grid_to_solver, bool ac) const {
     const int nb_sgen = nb();
     int bus_id_me, bus_id_solver;
     cplx_type tmp;
diff --git a/src/DataSGen.h b/src/DataSGen.h
index 2a99517..41ec3ac 100644
--- a/src/DataSGen.h
+++ b/src/DataSGen.h
@@ -159,7 +159,7 @@ class DataSGen: public DataGeneric
     void change_p(int sgen_id, real_type new_p, bool & need_reset);
     void change_q(int sgen_id, real_type new_q, bool & need_reset);
 
-    virtual void fillSbus(CplxVect & Sbus, const std::vector<int> & id_grid_to_solver);
+    virtual void fillSbus(CplxVect & Sbus, const std::vector<int> & id_grid_to_solver, bool ac) const ;
 
     void compute_results(const Eigen::Ref<const RealVect> & Va,
                          const Eigen::Ref<const RealVect> & Vm,
diff --git a/src/DataShunt.cpp b/src/DataShunt.cpp
index fe685fe..4c697b2 100644
--- a/src/DataShunt.cpp
+++ b/src/DataShunt.cpp
@@ -49,17 +49,19 @@ void DataShunt::set_state(DataShunt::StateRes & my_state )
 void DataShunt::fillYbus(std::vector<Eigen::Triplet<cplx_type> > & res,
                          bool ac,
                          const std::vector<int> & id_grid_to_solver,
-                         real_type sn_mva){
-    const int nb_shunt = static_cast<int>(q_mvar_.size());
+                         real_type sn_mva) const
+{
+    if(!ac) return; // no shunt in DC
+
+    const Eigen::Index nb_shunt = static_cast<int>(q_mvar_.size());
     cplx_type tmp;
     int bus_id_me, bus_id_solver;
-    for(int shunt_id=0; shunt_id < nb_shunt; ++shunt_id){
+    for(Eigen::Index shunt_id=0; shunt_id < nb_shunt; ++shunt_id){
         // i don't do anything if the shunt is disconnected
         if(!status_[shunt_id]) continue;
 
         // assign diagonal coefficient
-        tmp = {p_mw_(shunt_id), my_zero_};
-        if(ac) tmp += my_i * q_mvar_(shunt_id);
+        tmp = {p_mw_(shunt_id), -q_mvar_(shunt_id)};
 
         bus_id_me = bus_id_(shunt_id);
         bus_id_solver = id_grid_to_solver[bus_id_me];
@@ -71,49 +73,61 @@ void DataShunt::fillYbus(std::vector<Eigen::Triplet<cplx_type> > & res,
             throw std::runtime_error(exc_.str());
         }
         if(sn_mva != 1.) tmp /= sn_mva;
-        res.push_back(Eigen::Triplet<cplx_type> (bus_id_solver, bus_id_solver, -tmp));
+        res.push_back(Eigen::Triplet<cplx_type> (bus_id_solver, bus_id_solver, tmp));
     }
 }
 
-void DataShunt::fillSbus(CplxVect & Sbus, const std::vector<int> & id_grid_to_solver)  // in DC i need that
+void DataShunt::fillBp_Bpp(std::vector<Eigen::Triplet<real_type> > & Bp,
+                           std::vector<Eigen::Triplet<real_type> > & Bpp,
+                           const std::vector<int> & id_grid_to_solver,
+                           real_type sn_mva,
+                           FDPFMethod xb_or_bx) const
 {
-    // if(ac) return;  // in AC I do not do that
-    // std::cout << " ok i use this function" << std::endl;
-    // - bus[:, GS] / baseMVA  # in pandapower
-    // yish=gish+jbish -> so g is the MW !
-    const int nb_shunt = static_cast<int>(q_mvar_.size());
+    const Eigen::Index nb_shunt = static_cast<int>(q_mvar_.size());
+    real_type tmp;
     int bus_id_me, bus_id_solver;
-    for(int shunt_id=0; shunt_id < nb_shunt; ++shunt_id){
+    for(Eigen::Index shunt_id=0; shunt_id < nb_shunt; ++shunt_id){
         // i don't do anything if the shunt is disconnected
         if(!status_[shunt_id]) continue;
+
         bus_id_me = bus_id_(shunt_id);
         bus_id_solver = id_grid_to_solver[bus_id_me];
         if(bus_id_solver == _deactivated_bus_id){
-            throw std::runtime_error("GridModel::fillSbus: A shunt is connected to a disconnected bus.");
+            std::ostringstream exc_;
+            exc_ << "DataShunt::fillBp_Bpp: the shunt with id ";
+            exc_ << shunt_id;
+            exc_ << " is connected to a disconnected bus while being connected";
+            throw std::runtime_error(exc_.str());
         }
-        Sbus.coeffRef(bus_id_solver) -= p_mw_(shunt_id);
+        // assign diagonal coefficient
+        tmp = q_mvar_(shunt_id);
+        if(sn_mva != 1.) tmp /= sn_mva;
+        Bpp.push_back(Eigen::Triplet<real_type> (bus_id_solver, bus_id_solver, tmp));  // -(-tmp) [-tmp for the "correct" value, but then for Bpp i have -(-tmp)]
     }
 }
 
-void DataShunt::fillYbus_spmat(Eigen::SparseMatrix<cplx_type> & res, bool ac, const std::vector<int> & id_grid_to_solver){
+void DataShunt::fillSbus(CplxVect & Sbus, const std::vector<int> & id_grid_to_solver, bool ac) const  // in DC i need that
+{
+    if(ac) return;  // in AC I do not do that
+    // std::cout << " ok i use this function" << std::endl;
+    // - bus[:, GS] / baseMVA  # in pandapower
+    // yish=gish+jbish -> so g is the MW !
     const int nb_shunt = static_cast<int>(q_mvar_.size());
-    //TODO this is no more used!!!! see the other fillYbus
-    cplx_type tmp;
     int bus_id_me, bus_id_solver;
     for(int shunt_id=0; shunt_id < nb_shunt; ++shunt_id){
         // i don't do anything if the shunt is disconnected
         if(!status_[shunt_id]) continue;
-
-        // assign diagonal coefficient
-        tmp = p_mw_(shunt_id) + my_i * q_mvar_(shunt_id);
         bus_id_me = bus_id_(shunt_id);
         bus_id_solver = id_grid_to_solver[bus_id_me];
         if(bus_id_solver == _deactivated_bus_id){
-            throw std::runtime_error("GridModel::fillYbusShunt: A shunt is connected to a disconnected bus.");
+            throw std::runtime_error("GridModel::fillSbus: A shunt is connected to a disconnected bus.");
         }
-        res.coeffRef(bus_id_solver, bus_id_solver) -= tmp;
+        Sbus.coeffRef(bus_id_solver) -= p_mw_(shunt_id);
     }
-    //TODO this is no more used!!!! see the other fillYbus
+}
+
+void DataShunt::fillYbus_spmat(Eigen::SparseMatrix<cplx_type> & res, bool ac, const std::vector<int> & id_grid_to_solver){
+    throw std::runtime_error("DataShunt::fillYbus_spmat: should not be used anymore !");
 }
 
 void DataShunt::compute_results(const Eigen::Ref<const RealVect> & Va,
diff --git a/src/DataShunt.h b/src/DataShunt.h
index 26f856c..9a74be9 100644
--- a/src/DataShunt.h
+++ b/src/DataShunt.h
@@ -135,9 +135,14 @@ class DataShunt : public DataGeneric
     virtual void fillYbus(std::vector<Eigen::Triplet<cplx_type> > & res,
                           bool ac,
                           const std::vector<int> & id_grid_to_solver,
-                          real_type sn_mva);
+                          real_type sn_mva) const;
+    virtual void fillBp_Bpp(std::vector<Eigen::Triplet<real_type> > & Bp,
+                            std::vector<Eigen::Triplet<real_type> > & Bpp,
+                            const std::vector<int> & id_grid_to_solver,
+                            real_type sn_mva,
+                            FDPFMethod xb_or_bx) const;
     virtual void fillYbus_spmat(Eigen::SparseMatrix<cplx_type> & res, bool ac, const std::vector<int> & id_grid_to_solver);
-    virtual void fillSbus(CplxVect & Sbus, const std::vector<int> & id_grid_to_solver);  // in DC i need that
+    virtual void fillSbus(CplxVect & Sbus, const std::vector<int> & id_grid_to_solver, bool ac) const;  // in DC i need that
 
     void compute_results(const Eigen::Ref<const RealVect> & Va,
                          const Eigen::Ref<const RealVect> & Vm,
diff --git a/src/DataTrafo.cpp b/src/DataTrafo.cpp
index ead36ce..ab53604 100644
--- a/src/DataTrafo.cpp
+++ b/src/DataTrafo.cpp
@@ -166,7 +166,7 @@ void DataTrafo::fillYbus_spmat(Eigen::SparseMatrix<cplx_type> & res, bool ac, co
 void DataTrafo::fillYbus(std::vector<Eigen::Triplet<cplx_type> > & res,
                          bool ac,
                          const std::vector<int> & id_grid_to_solver,
-                         real_type sn_mva)
+                         real_type sn_mva) const
 {
     //TODO merge that with fillYbusBranch!
     //TODO template here instead of "if" for ac / dc
@@ -230,7 +230,7 @@ void DataTrafo::hack_Sbus_for_dc_phase_shifter(CplxVect & Sbus, bool ac, const s
         bus_id_solver_lv = id_grid_to_solver[bus_id_me];
         if(bus_id_solver_lv == _deactivated_bus_id){
             std::ostringstream exc_;
-            exc_ << "DataTrafo::fillSbus: the trafo with id ";
+            exc_ << "DataTrafo::hack_Sbus_for_dc_phase_shifter: the trafo with id ";
             exc_ << trafo_id;
             exc_ << " is connected (lv side) to a disconnected bus while being connected";
             throw std::runtime_error(exc_.str());
@@ -239,7 +239,7 @@ void DataTrafo::hack_Sbus_for_dc_phase_shifter(CplxVect & Sbus, bool ac, const s
         bus_id_solver_hv = id_grid_to_solver[bus_id_me];
         if(bus_id_solver_hv == _deactivated_bus_id){
             std::ostringstream exc_;
-            exc_ << "DataTrafo::fillSbus: the trafo with id ";
+            exc_ << "DataTrafo::hack_Sbus_for_dc_phase_shifter: the trafo with id ";
             exc_ << trafo_id;
             exc_ << " is connected (hv side) to a disconnected bus while being connected";
             throw std::runtime_error(exc_.str());
@@ -351,3 +351,103 @@ void DataTrafo::reset_results(){
     res_v_lv_ = RealVect();  // in kV
     res_a_lv_ = RealVect();  // in kA
 }
+
+void DataTrafo::fillBp_Bpp(std::vector<Eigen::Triplet<real_type> > & Bp,
+                          std::vector<Eigen::Triplet<real_type> > & Bpp,
+                          const std::vector<int> & id_grid_to_solver,
+                          real_type sn_mva,
+                          FDPFMethod xb_or_bx) const
+{
+
+    // For Bp
+    // temp_branch[:, BR_B] = zeros(nl)           ## zero out line charging shunts
+    // temp_branch[:, TAP] = ones(nl)             ## cancel out taps
+    // if alg == 2:                               ## if XB method
+    //    temp_branch[:, BR_R] = zeros(nl)       ## zero out line resistance
+
+    // For Bpp
+    // temp_branch[:, SHIFT] = zeros(nl)          ## zero out phase shifters
+    // if alg == 3:                               ## if BX method
+    //     temp_branch[:, BR_R] = zeros(nl)    ## zero out line resistance
+    const Eigen::Index nb_trafo = static_cast<int>(nb());
+    real_type yft_bp, ytf_bp, yff_bp, ytt_bp;
+    real_type yft_bpp, ytf_bpp, yff_bpp, ytt_bpp;
+
+    //diagonal coefficients
+    for(Eigen::Index tr_id=0; tr_id < nb_trafo; ++tr_id){
+        // i only add this if the powerline is connected
+        if(!status_[tr_id]) continue;
+
+        // get the from / to bus id
+        int bus_or_id_me = bus_hv_id_(tr_id);
+        int bus_or_solver_id = id_grid_to_solver[bus_or_id_me];
+        if(bus_or_solver_id == _deactivated_bus_id){
+            std::ostringstream exc_;
+            exc_ << "DataTrafo::fillBp_Bpp: the trafo with id ";
+            exc_ << tr_id;
+            exc_ << " is connected (hv side) to a disconnected bus while being connected";
+            throw std::runtime_error(exc_.str());
+        }
+        int bus_ex_id_me = bus_lv_id_(tr_id);
+        int bus_ex_solver_id = id_grid_to_solver[bus_ex_id_me];
+        if(bus_ex_solver_id == _deactivated_bus_id){
+            std::ostringstream exc_;
+            exc_ << "DataTrafo::fillBp_Bpp: the trafo with id ";
+            exc_ << tr_id;
+            exc_ << " is connected (lv side) to a disconnected bus while being connected";
+            throw std::runtime_error(exc_.str());
+        }
+
+        // get the coefficients
+        // tau is needed for Bpp
+        double tau_bpp = is_tap_hv_side_[tr_id] ? ratio_(tr_id) : 1. / ratio_(tr_id);
+        // for Bp we need shift
+        const real_type theta_shift = shift_(tr_id);
+        cplx_type eitheta_shift_bp  = {my_one_, my_zero_};  // exp(j  * alpha)
+        cplx_type emitheta_shift_bp = {my_one_, my_zero_};  // exp(-j * alpha)
+        if(theta_shift != 0.)
+        {
+            const real_type cos_theta = std::cos(theta_shift);
+            const real_type sin_theta = std::sin(theta_shift);
+            eitheta_shift_bp = {cos_theta, sin_theta};
+            emitheta_shift_bp = {cos_theta, -sin_theta};
+        }
+        // depending on XB or BX we define the y differently
+        cplx_type ys_bp, ys_bpp;
+        if(xb_or_bx==FDPFMethod::XB){
+            ys_bp = 1. / (0. + my_i * x_(tr_id));
+            ys_bpp = 1. / (r_(tr_id) + my_i * x_(tr_id));
+        }else if (xb_or_bx==FDPFMethod::BX){
+            ys_bp = 1. / (r_(tr_id) + my_i * x_(tr_id));
+            ys_bpp = 1. / (0. + my_i * x_(tr_id));
+        }else{
+            std::ostringstream exc_;
+            exc_ << "DataTrafo::fillBp_Bpp: unknown method for the FDPF powerflow for line id ";
+            exc_ << tr_id;
+            throw std::runtime_error(exc_.str());            
+        }
+
+        const real_type ys_bp_r = std::imag(ys_bp); 
+        yff_bp = ys_bp_r;
+        ytt_bp = ys_bp_r;
+        ytf_bp = -std::imag(ys_bp * emitheta_shift_bp);
+        yft_bp = -std::imag(ys_bp * eitheta_shift_bp);
+        const real_type ys_bpp_r = std::imag(ys_bpp); 
+        yff_bpp = (ys_bpp_r + std::imag(0.5 * my_i * h_(tr_id))) / (tau_bpp * tau_bpp);
+        ytt_bpp = ys_bpp_r + std::imag(0.5 * my_i * h_(tr_id));
+        ytf_bpp = -ys_bpp_r / tau_bpp;
+        yft_bpp = -ys_bpp_r / tau_bpp;
+
+        // and now add them
+        Bp.push_back(Eigen::Triplet<real_type> (bus_or_solver_id, bus_ex_solver_id, -yft_bp));
+        Bp.push_back(Eigen::Triplet<real_type> (bus_ex_solver_id, bus_or_solver_id, -ytf_bp));
+        Bp.push_back(Eigen::Triplet<real_type> (bus_or_solver_id, bus_or_solver_id, -yff_bp));
+        Bp.push_back(Eigen::Triplet<real_type> (bus_ex_solver_id, bus_ex_solver_id, -ytt_bp));
+
+        Bpp.push_back(Eigen::Triplet<real_type> (bus_or_solver_id, bus_ex_solver_id, -yft_bpp));
+        Bpp.push_back(Eigen::Triplet<real_type> (bus_ex_solver_id, bus_or_solver_id, -ytf_bpp));
+        Bpp.push_back(Eigen::Triplet<real_type> (bus_or_solver_id, bus_or_solver_id, -yff_bpp));
+        Bpp.push_back(Eigen::Triplet<real_type> (bus_ex_solver_id, bus_ex_solver_id, -ytt_bpp));
+
+    }
+}
\ No newline at end of file
diff --git a/src/DataTrafo.h b/src/DataTrafo.h
index 3400b9d..84f29e6 100644
--- a/src/DataTrafo.h
+++ b/src/DataTrafo.h
@@ -179,7 +179,12 @@ class DataTrafo : public DataGeneric
     virtual void fillYbus(std::vector<Eigen::Triplet<cplx_type> > & res,
                           bool ac,
                           const std::vector<int> & id_grid_to_solver,
-                          real_type sn_mva);
+                          real_type sn_mva) const;
+    virtual void fillBp_Bpp(std::vector<Eigen::Triplet<real_type> > & Bp,
+                            std::vector<Eigen::Triplet<real_type> > & Bpp,
+                            const std::vector<int> & id_grid_to_solver,
+                            real_type sn_mva,
+                            FDPFMethod xb_or_bx) const;
     virtual void hack_Sbus_for_dc_phase_shifter(CplxVect & Sbus, bool ac, const std::vector<int> & id_grid_to_solver);  // needed for dc mode
 
     void compute_results(const Eigen::Ref<const RealVect> & Va,
diff --git a/src/GridModel.cpp b/src/GridModel.cpp
index 614fdd1..c10476b 100644
--- a/src/GridModel.cpp
+++ b/src/GridModel.cpp
@@ -7,11 +7,15 @@
 // This file is part of LightSim2grid, LightSim2grid implements a c++ backend targeting the Grid2Op platform.
 
 #include "GridModel.h"
+#include "ChooseSolver.h"  // to avoid circular references
+
 
 GridModel::GridModel(const GridModel & other)
 {
     reset(true, true, true);
 
+    _ls_to_pp = other._ls_to_pp;
+
     init_vm_pu_ = other.init_vm_pu_;
     sn_mva_ = other.sn_mva_;
 
@@ -68,12 +72,15 @@ GridModel::GridModel(const GridModel & other)
     _solver.change_solver(other._solver.get_type());
     _dc_solver.change_solver(other._dc_solver.get_type());
     compute_results_ = other.compute_results_;
+    _dc_solver.set_gridmodel(this);
+    _solver.set_gridmodel(this);
 }
 
 //pickle
 GridModel::StateRes GridModel::get_state() const
 {
     std::vector<real_type> bus_vn_kv(bus_vn_kv_.begin(), bus_vn_kv_.end());
+    std::vector<int> ls_to_pp(_ls_to_pp.begin(), _ls_to_pp.end());
     int version_major = VERSION_MAJOR;
     int version_medium = VERSION_MEDIUM;
     int version_minor = VERSION_MINOR;
@@ -89,6 +96,7 @@ GridModel::StateRes GridModel::get_state() const
     GridModel::StateRes res(version_major,
                             version_medium,
                             version_minor,
+                            ls_to_pp,
                             init_vm_pu_,
                             sn_mva_,
                             bus_vn_kv,
@@ -128,30 +136,32 @@ void GridModel::set_state(GridModel::StateRes & my_state)
         exc_ << "It is not possible. Please reinstall it.";
         throw std::runtime_error(exc_.str());
     }
-    init_vm_pu_ = std::get<3>(my_state);
-    sn_mva_ = std::get<4>(my_state);
-    std::vector<real_type> & bus_vn_kv = std::get<5>(my_state);
-    std::vector<bool> & bus_status = std::get<6>(my_state);
+    std::vector<int> ls_to_pp_ = std::get<3>(my_state);
+    init_vm_pu_ = std::get<4>(my_state);
+    sn_mva_ = std::get<5>(my_state);
+    std::vector<real_type> & bus_vn_kv = std::get<6>(my_state);
+    std::vector<bool> & bus_status = std::get<7>(my_state);
 
     // powerlines
-    DataLine::StateRes & state_lines = std::get<7>(my_state);
+    DataLine::StateRes & state_lines = std::get<8>(my_state);
     // shunts
-    DataShunt::StateRes & state_shunts = std::get<8>(my_state);
+    DataShunt::StateRes & state_shunts = std::get<9>(my_state);
     // trafos
-    DataTrafo::StateRes & state_trafos = std::get<9>(my_state);
+    DataTrafo::StateRes & state_trafos = std::get<10>(my_state);
     // generators
-    DataGen::StateRes & state_gens = std::get<10>(my_state);
+    DataGen::StateRes & state_gens = std::get<11>(my_state);
     // loads
-    DataLoad::StateRes & state_loads = std::get<11>(my_state);
+    DataLoad::StateRes & state_loads = std::get<12>(my_state);
     // static gen
-    DataSGen::StateRes & state_sgens= std::get<12>(my_state);
+    DataSGen::StateRes & state_sgens= std::get<13>(my_state);
     // storage units
-    DataLoad::StateRes & state_storages = std::get<13>(my_state);
+    DataLoad::StateRes & state_storages = std::get<14>(my_state);
     // dc lines
-    DataDCLine::StateRes & state_dc_lines = std::get<14>(my_state);
+    DataDCLine::StateRes & state_dc_lines = std::get<15>(my_state);
 
     // assign it to this instance
 
+    _ls_to_pp =IntVect::Map(&ls_to_pp_[0], ls_to_pp_.size());
     // buses
     // 1. bus_vn_kv_
     bus_vn_kv_ = RealVect::Map(&bus_vn_kv[0], bus_vn_kv.size());
@@ -217,6 +227,8 @@ void GridModel::reset(bool reset_solver, bool reset_ac, bool reset_dc)
     if (reset_solver){
         _solver.reset();
         _dc_solver.reset();
+        _solver.set_gridmodel(this);
+        _dc_solver.set_gridmodel(this);
     }
     // std::cout << "GridModel::reset called" << std::endl;
 }
@@ -517,14 +529,14 @@ void GridModel::fillYbus(Eigen::SparseMatrix<cplx_type> & res, bool ac, const st
 void GridModel::fillSbus_me(CplxVect & Sbus, bool ac, const std::vector<int>& id_me_to_solver)
 {
     // init the Sbus vector
-    powerlines_.fillSbus(Sbus, id_me_to_solver);
-    trafos_.fillSbus(Sbus, id_me_to_solver);
-    shunts_.fillSbus(Sbus, id_me_to_solver);
-    loads_.fillSbus(Sbus, id_me_to_solver);
-    sgens_.fillSbus(Sbus, id_me_to_solver);
-    storages_.fillSbus(Sbus, id_me_to_solver);
-    generators_.fillSbus(Sbus, id_me_to_solver);
-    dc_lines_.fillSbus(Sbus, id_me_to_solver);
+    powerlines_.fillSbus(Sbus, id_me_to_solver, ac);
+    trafos_.fillSbus(Sbus, id_me_to_solver, ac);
+    shunts_.fillSbus(Sbus, id_me_to_solver, ac);
+    loads_.fillSbus(Sbus, id_me_to_solver, ac);
+    sgens_.fillSbus(Sbus, id_me_to_solver, ac);
+    storages_.fillSbus(Sbus, id_me_to_solver, ac);
+    generators_.fillSbus(Sbus, id_me_to_solver, ac);
+    dc_lines_.fillSbus(Sbus, id_me_to_solver, ac);
     if (sn_mva_ != 1.0) Sbus /= sn_mva_;
     // in dc mode, this is used for the phase shifter, this should not be divided by sn_mva_ !
     trafos_.hack_Sbus_for_dc_phase_shifter(Sbus, ac, id_me_to_solver);
@@ -826,3 +838,34 @@ void GridModel::update_topo(Eigen::Ref<Eigen::Array<bool, Eigen::Dynamic, Eigen:
                         &GridModel::deactivate_trafo
                         );
 }
+
+// for FDPF (implementation of the alg 2 method FDBX (FDXB will follow)  // TODO FDPF
+void GridModel::fillBp_Bpp(Eigen::SparseMatrix<real_type> & Bp, 
+                           Eigen::SparseMatrix<real_type> & Bpp, 
+                           FDPFMethod xb_or_bx) const
+{
+    // clear the matrices
+    const int nb_bus_solver = static_cast<int>(id_ac_solver_to_me_.size());
+    Bp = Eigen::SparseMatrix<real_type>(nb_bus_solver, nb_bus_solver);
+    Bpp = Eigen::SparseMatrix<real_type>(nb_bus_solver, nb_bus_solver);
+
+    // init the Bp and Bpp matrices for Fast Decoupled Powerflow  (TODO FDPF: optim when it's NOT needed just like for Ybus)
+    std::vector<Eigen::Triplet<real_type> > tripletList_Bp;
+    std::vector<Eigen::Triplet<real_type> > tripletList_Bpp;
+    tripletList_Bp.reserve(bus_vn_kv_.size() + 4 * powerlines_.nb() + 4 * trafos_.nb() + shunts_.nb());
+    tripletList_Bpp.reserve(bus_vn_kv_.size() + 4 * powerlines_.nb() + 4 * trafos_.nb() + shunts_.nb());
+    // run through the grid and get the parameters to fill them
+    powerlines_.fillBp_Bpp(tripletList_Bp, tripletList_Bpp, id_me_to_ac_solver_, sn_mva_, xb_or_bx);
+    shunts_.fillBp_Bpp(tripletList_Bp, tripletList_Bpp, id_me_to_ac_solver_, sn_mva_, xb_or_bx);
+    trafos_.fillBp_Bpp(tripletList_Bp, tripletList_Bpp, id_me_to_ac_solver_, sn_mva_, xb_or_bx);
+    loads_.fillBp_Bpp(tripletList_Bp, tripletList_Bpp, id_me_to_ac_solver_, sn_mva_, xb_or_bx);
+    sgens_.fillBp_Bpp(tripletList_Bp, tripletList_Bpp, id_me_to_ac_solver_, sn_mva_, xb_or_bx);
+    storages_.fillBp_Bpp(tripletList_Bp, tripletList_Bpp, id_me_to_ac_solver_, sn_mva_, xb_or_bx);
+    generators_.fillBp_Bpp(tripletList_Bp, tripletList_Bpp, id_me_to_ac_solver_, sn_mva_, xb_or_bx);
+    dc_lines_.fillBp_Bpp(tripletList_Bp, tripletList_Bpp, id_me_to_ac_solver_, sn_mva_, xb_or_bx);
+    // now make the matrices effectively
+    Bp.setFromTriplets(tripletList_Bp.begin(), tripletList_Bp.end());
+    Bp.makeCompressed();
+    Bpp.setFromTriplets(tripletList_Bpp.begin(), tripletList_Bpp.end());
+    Bpp.makeCompressed();
+}
diff --git a/src/GridModel.h b/src/GridModel.h
index a731b6e..a32e092 100644
--- a/src/GridModel.h
+++ b/src/GridModel.h
@@ -36,18 +36,23 @@
 #include "DataSGen.h"
 #include "DataDCLine.h"
 
-
 // import newton raphson solvers using different linear algebra solvers
 #include "ChooseSolver.h"
+// class ChooseSolver;
+// enum class SolverType;
 
 //TODO implement a BFS check to make sure the Ymatrix is "connected" [one single component]
 class GridModel : public DataGeneric
 {
+    public:  // can be modified python side
+        IntVect _ls_to_pp;  // for converter from bus in lightsim2grid index to bus in pandapower index
+
     public:
         typedef std::tuple<
                 int, // version major
                 int, // version medium
                 int, // version minor
+                std::vector<int>, // ls_to_pp
                 real_type,  // init_vm_pu
                 real_type, //sn_mva
                 std::vector<real_type>,  // bus_vn_kv
@@ -72,6 +77,7 @@ class GridModel : public DataGeneric
 
         GridModel():need_reset_(true), topo_changed_(true), compute_results_(true), init_vm_pu_(1.04), sn_mva_(1.0){
             _dc_solver.change_solver(SolverType::DC);
+            _solver.set_gridmodel(this);
         }
         GridModel(const GridModel & other);
         GridModel copy() const{
@@ -132,6 +138,18 @@ class GridModel : public DataGeneric
                              ){
             powerlines_.init(branch_r, branch_x, branch_h, branch_from_id, branch_to_id);
         }
+        void init_powerlines_full(const RealVect & branch_r,
+                                  const RealVect & branch_x,
+                                  const CplxVect & branch_h_or,
+                                  const CplxVect & branch_h_ex,
+                                  const Eigen::VectorXi & branch_from_id,
+                                  const Eigen::VectorXi & branch_to_id
+                             ){
+            powerlines_.init(branch_r, branch_x, branch_h_or,
+                             branch_h_ex, branch_from_id, 
+                             branch_to_id);
+        }
+
         void init_shunt(const RealVect & shunt_p_mw,
                         const RealVect & shunt_q_mvar,
                         const Eigen::VectorXi & shunt_bus_id){
@@ -491,6 +509,24 @@ class GridModel : public DataGeneric
             fillSbus_me(res, ac, id_me_to_solver);
         }
 
+        //for FDPF
+        void fillBp_Bpp(Eigen::SparseMatrix<real_type> & Bp, 
+                        Eigen::SparseMatrix<real_type> & Bpp, 
+                        FDPFMethod xb_or_bx) const;
+
+        Eigen::SparseMatrix<real_type> debug_get_Bp_python(FDPFMethod xb_or_bx){
+            Eigen::SparseMatrix<real_type> Bp;
+            Eigen::SparseMatrix<real_type> Bpp;
+            fillBp_Bpp(Bp, Bpp, xb_or_bx);
+            return Bp;
+        }
+        Eigen::SparseMatrix<real_type> debug_get_Bpp_python(FDPFMethod xb_or_bx){
+            Eigen::SparseMatrix<real_type> Bp;
+            Eigen::SparseMatrix<real_type> Bpp;
+            fillBp_Bpp(Bp, Bpp, xb_or_bx);
+            return Bpp;
+        }
+
     protected:
     // add method to change topology, change ratio of transformers, change
 
@@ -607,6 +643,7 @@ class GridModel : public DataGeneric
     protected:
         // member of the grid
         // static const int _deactivated_bus_id;
+
         bool need_reset_;
         bool topo_changed_;
         bool compute_results_;
diff --git a/src/Solvers.cpp b/src/Solvers.cpp
new file mode 100644
index 0000000..a7ea0d7
--- /dev/null
+++ b/src/Solvers.cpp
@@ -0,0 +1,37 @@
+// Copyright (c) 2023, RTE (https://www.rte-france.com)
+// See AUTHORS.txt
+// This Source Code Form is subject to the terms of the Mozilla Public License, version 2.0.
+// If a copy of the Mozilla Public License, version 2.0 was not distributed with this file,
+// you can obtain one at http://mozilla.org/MPL/2.0/.
+// SPDX-License-Identifier: MPL-2.0
+// This file is part of LightSim2grid, LightSim2grid implements a c++ backend targeting the Grid2Op platform.
+
+#include "Solvers.h"
+#include "GridModel.h"
+
+// these functions use the _gridmodel that is a not a real type in the template class declaration.
+// this is why i need to define them here for every specialization.
+
+template<class LinearSolver, FDPFMethod XB_BX>
+void BaseFDPFSolver<LinearSolver, XB_BX>::fillBp_Bpp(Eigen::SparseMatrix<real_type> & Bp, Eigen::SparseMatrix<real_type> & Bpp) const
+{
+    _gridmodel->fillBp_Bpp(Bp, Bpp, XB_BX);
+}
+
+template void FDPF_XB_SparseLUSolver::fillBp_Bpp(Eigen::SparseMatrix<real_type> & Bp, Eigen::SparseMatrix<real_type> & Bpp) const;
+template void FDPF_BX_SparseLUSolver::fillBp_Bpp(Eigen::SparseMatrix<real_type> & Bp, Eigen::SparseMatrix<real_type> & Bpp) const;
+
+#ifdef KLU_SOLVER_AVAILABLE
+template void FDPF_XB_KLUSolver::fillBp_Bpp(Eigen::SparseMatrix<real_type> & Bp, Eigen::SparseMatrix<real_type> & Bpp) const;
+template void FDPF_BX_KLUSolver::fillBp_Bpp(Eigen::SparseMatrix<real_type> & Bp, Eigen::SparseMatrix<real_type> & Bpp) const;
+#endif  // KLU_SOLVER_AVAILABLE
+
+#ifdef NICSLU_SOLVER_AVAILABLE
+template void FDPF_XB_NICSLUSolver::fillBp_Bpp(Eigen::SparseMatrix<real_type> & Bp, Eigen::SparseMatrix<real_type> & Bpp) const
+template void FDPF_BX_NICSLUSolver::fillBp_Bpp(Eigen::SparseMatrix<real_type> & Bp, Eigen::SparseMatrix<real_type> & Bpp) const;
+#endif  // NICSLU_SOLVER_AVAILABLE
+
+#ifdef CKTSO_SOLVER_AVAILABLE
+template void FDPF_XB_CKTSOSolver::fillBp_Bpp(Eigen::SparseMatrix<real_type> & Bp, Eigen::SparseMatrix<real_type> & Bpp) const;
+template void FDPF_BX_CKTSOSolver::fillBp_Bpp(Eigen::SparseMatrix<real_type> & Bp, Eigen::SparseMatrix<real_type> & Bpp) const;
+#endif  // CKTSO_SOLVER_AVAILABLE
diff --git a/src/Solvers.h b/src/Solvers.h
index 350c254..4adaac8 100644
--- a/src/Solvers.h
+++ b/src/Solvers.h
@@ -1,4 +1,4 @@
-// Copyright (c) 2020, RTE (https://www.rte-france.com)
+// Copyright (c) 2020-2023, RTE (https://www.rte-france.com)
 // See AUTHORS.txt
 // This Source Code Form is subject to the terms of the Mozilla Public License, version 2.0.
 // If a copy of the Mozilla Public License, version 2.0 was not distributed with this file,
@@ -9,6 +9,7 @@
 #include "BaseNRSolver.h"
 #include "BaseNRSolverSingleSlack.h"
 #include "DCSolver.h"
+#include "BaseFDPFSolver.h"
 
 #include "SparseLUSolver.h"
 #include "KLUSolver.h"
@@ -21,6 +22,9 @@ typedef BaseNRSolver<SparseLULinearSolver> SparseLUSolver;
 typedef BaseNRSolverSingleSlack<SparseLULinearSolver> SparseLUSolverSingleSlack;
 /** Solver based on Newton Raphson, using the SparseLU decomposition of Eigen, only suitable for the DC approximation**/
 typedef BaseDCSolver<SparseLULinearSolver> DCSolver;
+/** Solver based on Fast Decoupled, using the SparseLU decomposition of Eigen**/
+typedef BaseFDPFSolver<SparseLULinearSolver, FDPFMethod::XB> FDPF_XB_SparseLUSolver;
+typedef BaseFDPFSolver<SparseLULinearSolver, FDPFMethod::BX> FDPF_BX_SparseLUSolver;
 
 #ifdef KLU_SOLVER_AVAILABLE
     /** Solver based on Newton Raphson, using the KLU linear solver**/
@@ -29,6 +33,9 @@ typedef BaseDCSolver<SparseLULinearSolver> DCSolver;
     typedef BaseNRSolverSingleSlack<KLULinearSolver> KLUSolverSingleSlack;
     /** Solver based on Newton Raphson, using the KLU linear solver, only suitable for the DC approximation**/
     typedef BaseDCSolver<KLULinearSolver> KLUDCSolver;
+    /** Solver based on Fast Decoupled, using the KLU linear solver**/
+    typedef BaseFDPFSolver<KLULinearSolver, FDPFMethod::XB> FDPF_XB_KLUSolver;
+    typedef BaseFDPFSolver<KLULinearSolver, FDPFMethod::BX> FDPF_BX_KLUSolver;
 #elif defined(_READ_THE_DOCS)
     // hack to display accurately the doc in read the doc even if the models are not compiled
     /** Solver based on Newton Raphson, using the KLU linear solver**/
@@ -37,6 +44,9 @@ typedef BaseDCSolver<SparseLULinearSolver> DCSolver;
     class KLUSolverSingleSlack : public SparseLUSolverSingleSlack {};
     /** Solver based on Newton Raphson, using the KLU linear solver, only suitable for the DC approximation**/
     class KLUDCSolver : public DCSolver {};
+    /** Solver based on Fast Decoupled, using the KLU linear solver**/
+    class FDPF_XB_KLUSolver : public FDPF_XB_SparseLUSolver {};
+    class FDPF_BX_KLUSolver : public FDPF_BX_SparseLUSolver {};
 #endif  // KLU_SOLVER_AVAILABLE
 
 #ifdef NICSLU_SOLVER_AVAILABLE
@@ -46,14 +56,20 @@ typedef BaseDCSolver<SparseLULinearSolver> DCSolver;
     typedef BaseNRSolverSingleSlack<NICSLULinearSolver> NICSLUSolverSingleSlack;
     /** Solver based on Newton Raphson, using the NICSLU linear solver (needs a specific license), only suitable for the DC approximation**/
     typedef BaseDCSolver<NICSLULinearSolver> NICSLUDCSolver;
+    /** Solver based on Fast Decoupled, using the NICSLU linear solver (needs a specific license)**/
+    typedef BaseFDPFSolver<NICSLULinearSolver, FDPFMethod::XB> FDPF_XB_NICSLUSolver;
+    typedef BaseFDPFSolver<NICSLULinearSolver, FDPFMethod::BX> FDPF_BX_NICSLUSolver;
 #elif defined(_READ_THE_DOCS)
     // hack to display accurately the doc in read the doc even if the models are not compiled
-    /** Solver based on Newton Raphson, using the KLU linear solver**/
+    /** Solver based on Newton Raphson, using the NICSLU linear solver (needs a specific license)**/
     class NICSLUSolver : public SparseLUSolver{};
-    /** Solver based on Newton Raphson, using the KLU linear solver, do not consider multiple slack bus**/
+    /** Solver based on Newton Raphson, using the NICSLU linear solver (needs a specific license), do not consider multiple slack bus**/
     class NICSLUSolverSingleSlack : public SparseLUSolverSingleSlack{};
-    /** Solver based on Newton Raphson, using the KLU linear solver, only suitable for the DC approximation**/
+    /** Solver based on Newton Raphson, using the NICSLU linear solver (needs a specific license), only suitable for the DC approximation**/
     class NICSLUDCSolver : public DCSolver{};
+    /** Solver based on Fast Decoupled, using the NICSLU linear solver (needs a specific license)**/
+    class FDPF_XB_NICSLUSolver : public FDPF_XB_SparseLUSolver {};
+    class FDPF_BX_NICSLUSolver : public FDPF_BX_SparseLUSolver {};
 #endif  // NICSLU_SOLVER_AVAILABLE
 
 #ifdef CKTSO_SOLVER_AVAILABLE
@@ -63,12 +79,18 @@ typedef BaseDCSolver<SparseLULinearSolver> DCSolver;
     typedef BaseNRSolverSingleSlack<CKTSOLinearSolver> CKTSOSolverSingleSlack;
     /** Solver based on Newton Raphson, using the CKTSO linear solver (needs a specific license), only suitable for the DC approximation**/
     typedef BaseDCSolver<CKTSOLinearSolver> CKTSODCSolver;
+    /** Solver based on Fast Decoupled, using the CKTSO linear solver (needs a specific license)**/
+    typedef BaseFDPFSolver<CKTSOLinearSolver, FDPFMethod::XB> FDPF_XB_CKTSOSolver;
+    typedef BaseFDPFSolver<CKTSOLinearSolver, FDPFMethod::BX> FDPF_BX_CKTSOSolver;
 #elif defined(_READ_THE_DOCS)
     // hack to display accurately the doc in read the doc even if the models are not compiled
-    /** Solver based on Newton Raphson, using the KLU linear solver**/
+    /** Solver based on Newton Raphson, using the CKTSO linear solver (needs a specific license)**/
     class CKTSOSolver : public SparseLUSolver{};
-    /** Solver based on Newton Raphson, using the KLU linear solver, do not consider multiple slack bus**/
+     /** Solver based on Newton Raphson, using the CKTSO linear solver (needs a specific license), do not consider multiple slack bus**/
     class CKTSOSolverSingleSlack : public SparseLUSolverSingleSlack{};
-    /** Solver based on Newton Raphson, using the KLU linear solver, only suitable for the DC approximation**/
+    /** Solver based on Newton Raphson, using the CKTSO linear solver (needs a specific license), only suitable for the DC approximation**/
     class CKTSODCSolver : public DCSolver{};
+    /** Solver based on Fast Decoupled, using the CKTSO linear solver (needs a specific license)**/
+    class FDPF_XB_CKTSOSolver : public FDPF_XB_SparseLUSolver {};
+    class FDPF_BX_CKTSOSolver : public FDPF_BX_SparseLUSolver {};
 #endif  // CKTSO_SOLVER_AVAILABLE
diff --git a/src/Utils.h b/src/Utils.h
index 319b40d..f8a6d02 100644
--- a/src/Utils.h
+++ b/src/Utils.h
@@ -29,6 +29,7 @@ typedef std::tuple<Eigen::Ref<const EigenPythonNumType>,
                    Eigen::Ref<const EigenPythonNumType>,
                    Eigen::Ref<const EigenPythonNumType>,
                    Eigen::Ref<const EigenPythonNumType> > tuple4d;
+typedef Eigen::Matrix<int, Eigen::Dynamic, 1> IntVect;
 typedef Eigen::Matrix<real_type, Eigen::Dynamic, 1> RealVect;
 typedef Eigen::Matrix<cplx_type, Eigen::Dynamic, 1> CplxVect;
 
diff --git a/src/help_fun_msg.cpp b/src/help_fun_msg.cpp
index fb4d857..6a3fe40 100644
--- a/src/help_fun_msg.cpp
+++ b/src/help_fun_msg.cpp
@@ -166,6 +166,40 @@ const std::string DocSolver::DCSolver =  R"mydelimiter(
 
 )mydelimiter";
 
+const std::string DocSolver::FDPF_XB_SparseLUSolver =  R"mydelimiter(
+    Default implementation of the Fast Decoupled Powerflow solver (XB version: "alg 2" / "fdxb"  in pypower / pandapower), it uses the default Eigen sparse lu decomposition for 
+    its underlying sparse matrix manipulation.
+
+    See :ref:`available-powerflow-solvers` for more information on how to use it.
+
+    .. note::
+
+        In the enum :attr:`lightsim2grid.solver.SolverType`, it is called `FDPF_SparseLU` 
+        
+        You can use it with:
+        
+        - `env_lightsim.backend.set_solver_type(lightsim2grid.solver.FDPF_SparseLU)` after creation
+        - `LightSimBackend(solver_type=lightsim2grid.solver.FDPF_SparseLU)` at creation time
+
+)mydelimiter";
+
+const std::string DocSolver::FDPF_BX_SparseLUSolver =  R"mydelimiter(
+    Default implementation of the Fast Decoupled Powerflow solver (XB version: "alg 3" / "fdbx"  in pypower / pandapower), it uses the default Eigen sparse lu decomposition for 
+    its underlying sparse matrix manipulation.
+
+    See :ref:`available-powerflow-solvers` for more information on how to use it.
+
+    .. note::
+
+        In the enum :attr:`lightsim2grid.solver.SolverType`, it is called `FDPF_SparseLU` 
+        
+        You can use it with:
+        
+        - `env_lightsim.backend.set_solver_type(lightsim2grid.solver.FDPF_SparseLU)` after creation
+        - `LightSimBackend(solver_type=lightsim2grid.solver.FDPF_SparseLU)` at creation time
+
+)mydelimiter";
+
 const std::string DocSolver::KLUSolver = R"mydelimiter(
     This classes implements the Newton Raphson algorithm, allowing for distributed slack and using the faster KLU solver available in the SuiteSparse library
     for the linear algebra (can be unavailable if you build lightsim2grid from source). It is usually faster than the :class:`lightsim2grid.solver.SparseLUSolver`.
@@ -229,6 +263,40 @@ const std::string DocSolver::KLUDCSolver = R"mydelimiter(
 
 )mydelimiter";
 
+const std::string DocSolver::FDPF_XB_KLUSolver =  R"mydelimiter(
+    Default implementation of the Fast Decoupled Powerflow solver (XB version: "alg 2" / "fdbx"  in pypower / pandapower), it uses the fast KLU library for 
+    its underlying sparse matrix manipulation.
+
+    See :ref:`available-powerflow-solvers` for more information on how to use it.
+
+    .. note::
+
+        In the enum :attr:`lightsim2grid.solver.SolverType`, it is called `FDPF_KLU` 
+        
+        You can use it with:
+        
+        - `env_lightsim.backend.set_solver_type(lightsim2grid.solver.FDPF_KLU)` after creation
+        - `LightSimBackend(solver_type=lightsim2grid.solver.FDPF_KLU)` at creation time
+
+)mydelimiter";
+
+const std::string DocSolver::FDPF_BX_KLUSolver =  R"mydelimiter(
+    Default implementation of the Fast Decoupled Powerflow solver (XB version: "alg 3" / "fdxb"  in pypower / pandapower), it uses the fast KLU library for 
+    its underlying sparse matrix manipulation.
+
+    See :ref:`available-powerflow-solvers` for more information on how to use it.
+
+    .. note::
+
+        In the enum :attr:`lightsim2grid.solver.SolverType`, it is called `FDPF_KLU` 
+        
+        You can use it with:
+        
+        - `env_lightsim.backend.set_solver_type(lightsim2grid.solver.FDPF_KLU)` after creation
+        - `LightSimBackend(solver_type=lightsim2grid.solver.FDPF_KLU)` at creation time
+
+)mydelimiter";
+
 const std::string DocSolver::NICSLUSolver = R"mydelimiter(
     This classes implements the Newton Raphson algorithm, allowing for distributed slack and using the faster NICSLU solver available in the NICSLU library
     for the linear algebra. It is usually faster than the :class:`lightsim2grid.solver.SparseLUSolver`. (requires a build from source)
@@ -244,6 +312,10 @@ const std::string DocSolver::NICSLUSolver = R"mydelimiter(
         - `env_lightsim.backend.set_solver_type(lightsim2grid.solver.NICSLU)` after creation
         - `LightSimBackend(solver_type=lightsim2grid.solver.NICSLU)` at creation time
 
+    .. warning::
+        
+        Use this solver requires a compilation of lightsim2grid from source (see readme) AND an appropriate license for nicslu.
+
     .. note::
 
         NICSLU is available at https://github.com/chenxm1986/nicslu
@@ -264,7 +336,11 @@ const std::string DocSolver::NICSLUSolverSingleSlack = R"mydelimiter(
         
         - `env_lightsim.backend.set_solver_type(lightsim2grid.solver.NICSLUSingleSlack)` after creation
         - `LightSimBackend(solver_type=lightsim2grid.solver.NICSLUSingleSlack)` at creation time    
-    
+
+    .. warning::
+        
+        Use this solver requires a compilation of lightsim2grid from source (see readme) AND an appropriate license for nicslu.
+
     .. note::
 
         NICSLU is available at https://github.com/chenxm1986/nicslu
@@ -292,12 +368,66 @@ const std::string DocSolver::NICSLUDCSolver = R"mydelimiter(
 
         Otherwise, it is used internally to find good starting point to intialize the real AC solver.
 
+    .. warning::
+        
+        Use this solver requires a compilation of lightsim2grid from source (see readme) AND an appropriate license for nicslu.
+
     .. note::
 
         NICSLU is available at https://github.com/chenxm1986/nicslu
  
 )mydelimiter";
 
+const std::string DocSolver::FDPF_XB_NICSLUSolver =  R"mydelimiter(
+    Default implementation of the Fast Decoupled Powerflow solver (XB version: "alg 2" / "fdxb"  in pypower / pandapower), it uses the fast NICSLU library for 
+    its underlying sparse matrix manipulation.
+
+    See :ref:`available-powerflow-solvers` for more information on how to use it.
+
+    .. note::
+
+        In the enum :attr:`lightsim2grid.solver.SolverType`, it is called `FDPF_NICSLU` 
+        
+        You can use it with:
+        
+        - `env_lightsim.backend.set_solver_type(lightsim2grid.solver.FDPF_NICSLU)` after creation
+        - `LightSimBackend(solver_type=lightsim2grid.solver.FDPF_NICSLU)` at creation time    
+
+    .. warning::
+        
+        Use this solver requires a compilation of lightsim2grid from source (see readme) AND an appropriate license for nicslu.
+
+    .. note::
+
+        NICSLU is available at https://github.com/chenxm1986/nicslu
+
+)mydelimiter";
+
+const std::string DocSolver::FDPF_BX_NICSLUSolver =  R"mydelimiter(
+    Default implementation of the Fast Decoupled Powerflow solver (XB version: "alg 3" / "fdbx"  in pypower / pandapower), it uses the fast NICSLU library for 
+    its underlying sparse matrix manipulation.
+
+    See :ref:`available-powerflow-solvers` for more information on how to use it.
+
+    .. note::
+
+        In the enum :attr:`lightsim2grid.solver.SolverType`, it is called `FDPF_NICSLU` 
+        
+        You can use it with:
+        
+        - `env_lightsim.backend.set_solver_type(lightsim2grid.solver.FDPF_NICSLU)` after creation
+        - `LightSimBackend(solver_type=lightsim2grid.solver.FDPF_NICSLU)` at creation time    
+
+    .. warning::
+        
+        Use this solver requires a compilation of lightsim2grid from source (see readme) AND an appropriate license for nicslu.
+
+    .. note::
+
+        NICSLU is available at https://github.com/chenxm1986/nicslu
+
+)mydelimiter";
+
 const std::string DocSolver::CKTSOSolver = R"mydelimiter(
     This classes implements the Newton Raphson algorithm, allowing for distributed slack and using the faster CKTSO solver available in the CKTSO library
     for the linear algebra (requires a build from source)
@@ -367,6 +497,56 @@ const std::string DocSolver::CKTSODCSolver = R"mydelimiter(
 
 )mydelimiter";
 
+const std::string DocSolver::FDPF_XB_CKTSOSolver =  R"mydelimiter(
+    Default implementation of the Fast Decoupled Powerflow solver (XB version: "alg 2" / "fdxb"  in pypower / pandapower), it uses the fast CKTSO library for 
+    its underlying sparse matrix manipulation.
+
+    See :ref:`available-powerflow-solvers` for more information on how to use it.
+
+    .. note::
+
+        In the enum :attr:`lightsim2grid.solver.SolverType`, it is called `FDPF_CKTSO` 
+        
+        You can use it with:
+        
+        - `env_lightsim.backend.set_solver_type(lightsim2grid.solver.FDPF_CKTSO)` after creation
+        - `LightSimBackend(solver_type=lightsim2grid.solver.FDPF_CKTSO)` at creation time    
+
+    .. warning::
+        
+        Use this solver requires a compilation of lightsim2grid from source (see readme) AND an appropriate license for cktso.
+
+    .. note::
+
+        CKTSO is available at https://github.com/chenxm1986/cktso
+
+)mydelimiter";
+
+const std::string DocSolver::FDPF_BX_CKTSOSolver =  R"mydelimiter(
+    Default implementation of the Fast Decoupled Powerflow solver (XB version: "alg 3" / "fdbx"  in pypower / pandapower), it uses the fast CKTSO library for 
+    its underlying sparse matrix manipulation.
+
+    See :ref:`available-powerflow-solvers` for more information on how to use it.
+
+    .. note::
+
+        In the enum :attr:`lightsim2grid.solver.SolverType`, it is called `FDPF_CKTSO` 
+        
+        You can use it with:
+        
+        - `env_lightsim.backend.set_solver_type(lightsim2grid.solver.FDPF_CKTSO)` after creation
+        - `LightSimBackend(solver_type=lightsim2grid.solver.FDPF_CKTSO)` at creation time    
+
+    .. warning::
+        
+        Use this solver requires a compilation of lightsim2grid from source (see readme) AND an appropriate license for cktso.
+
+    .. note::
+
+        CKTSO is available at https://github.com/chenxm1986/cktso
+
+)mydelimiter";
+
 const std::string DocSolver::GaussSeidelSolver = R"mydelimiter(
     Default implementation of the "Gauss Seidel" powerflow solver. We do not recommend to use it as the Newton Raphson based solvers
     are usually much (much) faster.
diff --git a/src/help_fun_msg.h b/src/help_fun_msg.h
index 113f1dd..17fa719 100644
--- a/src/help_fun_msg.h
+++ b/src/help_fun_msg.h
@@ -32,18 +32,26 @@ struct DocSolver
     static const std::string SparseLUSolver;
     static const std::string SparseLUSolverSingleSlack;
     static const std::string DCSolver;
+    static const std::string FDPF_XB_SparseLUSolver;
+    static const std::string FDPF_BX_SparseLUSolver;
 
     static const std::string KLUSolver;
     static const std::string KLUSolverSingleSlack;
     static const std::string KLUDCSolver;
+    static const std::string FDPF_XB_KLUSolver;
+    static const std::string FDPF_BX_KLUSolver;
 
     static const std::string NICSLUSolver;
     static const std::string NICSLUSolverSingleSlack;
     static const std::string NICSLUDCSolver;
+    static const std::string FDPF_XB_NICSLUSolver;
+    static const std::string FDPF_BX_NICSLUSolver;
 
     static const std::string CKTSOSolver;
     static const std::string CKTSOSolverSingleSlack;
     static const std::string CKTSODCSolver;
+    static const std::string FDPF_XB_CKTSOSolver;
+    static const std::string FDPF_BX_CKTSOSolver;
 
     static const std::string GaussSeidelSolver;
     static const std::string GaussSeidelSynchSolver;
diff --git a/src/main.cpp b/src/main.cpp
index 4b061ba..75e4e08 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -23,6 +23,12 @@ namespace py = pybind11;
 PYBIND11_MODULE(lightsim2grid_cpp, m)
 {
 
+    // solver method for FDPF
+    py::enum_<FDPFMethod>(m, "FDPFMethod", "This enum controls the type of method you can use for Fast Decoupled Powerflow (XB or BX)")
+        .value("XB", FDPFMethod::XB, "denotes the XB method")
+        .value("BX", FDPFMethod::BX, "denotes the BX method")
+        .export_values();
+
     // solvers
     py::enum_<SolverType>(m, "SolverType", "This enum controls the solver you want to use.")
         .value("GaussSeidel", SolverType::GaussSeidel, "denotes the :class:`lightsim2grid.solver.GaussSeidelSolver`")
@@ -39,6 +45,14 @@ PYBIND11_MODULE(lightsim2grid_cpp, m)
         .value("CKTSO", SolverType::CKTSO, "denotes the :class:`lightsim2grid.solver.CKTSOSolver`")
         .value("CKTSOSingleSlack", SolverType::CKTSOSingleSlack, "denotes the :class:`lightsim2grid.solver.CKTSOSolverSingleSlack`")
         .value("CKTSODC", SolverType::CKTSODC, "denotes the :class:`lightsim2grid.solver.CKTSODCSolver`")
+        .value("FDPF_XB_SparseLU", SolverType::FDPF_XB_SparseLU, "denotes the :class:`lightsim2grid.solver.FDPF_XB_SparseLUSolver`")
+        .value("FDPF_BX_SparseLU", SolverType::FDPF_BX_SparseLU, "denotes the :class:`lightsim2grid.solver.FDPF_BX_SparseLUSolver`")
+        .value("FDPF_XB_KLU", SolverType::FDPF_XB_KLU, "denotes the :class:`lightsim2grid.solver.FDPF_XB_KLUSolver`")
+        .value("FDPF_BX_KLU", SolverType::FDPF_BX_KLU, "denotes the :class:`lightsim2grid.solver.FDPF_BX_KLUSolver`")
+        .value("FDPF_XB_NICSLU", SolverType::FDPF_XB_NICSLU, "denotes the :class:`lightsim2grid.solver.FDPF_XB_NICSLUSolver`")
+        .value("FDPF_BX_NICSLU", SolverType::FDPF_BX_NICSLU, "denotes the :class:`lightsim2grid.solver.FDPF_BX_NICSLUSolver`")
+        .value("FDPF_XB_CKTSO", SolverType::FDPF_XB_CKTSO, "denotes the :class:`lightsim2grid.solver.FDPF_XB_CKTSOSolver`")
+        .value("FDPF_BX_CKTSO", SolverType::FDPF_BX_CKTSO, "denotes the :class:`lightsim2grid.solver.FDPF_BX_CKTSOSolver`")
         .export_values();
 
     py::enum_<ErrorType>(m, "ErrorType", "This enum controls the error encountered in the solver")
@@ -95,6 +109,36 @@ PYBIND11_MODULE(lightsim2grid_cpp, m)
         .def("get_timers", &DCSolver::get_timers, DocSolver::get_timers.c_str())  // returns the timers corresponding to times the solver spent in different part
         .def("solve", &DCSolver::compute_pf, py::call_guard<py::gil_scoped_release>(), DocSolver::compute_pf.c_str());  // perform the newton raphson optimization
 
+    py::class_<FDPF_XB_SparseLUSolver>(m, "FDPF_XB_SparseLUSolver", DocSolver::FDPF_XB_SparseLUSolver.c_str())
+        .def(py::init<>())
+        .def("get_Va", &FDPF_XB_SparseLUSolver::get_Va, DocSolver::get_Va.c_str())  // get the voltage angle vector (vector of double)
+        .def("get_Vm", &FDPF_XB_SparseLUSolver::get_Vm, DocSolver::get_Vm.c_str())  // get the voltage magnitude vector (vector of double)
+        .def("get_V", &FDPF_XB_SparseLUSolver::get_V, DocSolver::get_V.c_str()) 
+        .def("get_error", &FDPF_XB_SparseLUSolver::get_error, DocSolver::get_error.c_str())  // get the error message, see the definition of "err_" for more information
+        .def("get_nb_iter", &FDPF_XB_SparseLUSolver::get_nb_iter, DocSolver::get_nb_iter.c_str())  // return the number of iteration performed at the last optimization
+        .def("reset", &FDPF_XB_SparseLUSolver::reset, DocSolver::reset.c_str())  // reset the solver to its original state
+        .def("converged", &FDPF_XB_SparseLUSolver::converged, DocSolver::converged.c_str())  // whether the solver has converged
+        .def("compute_pf", &FDPF_XB_SparseLUSolver::compute_pf, py::call_guard<py::gil_scoped_release>(), DocSolver::compute_pf.c_str())  // compute the powerflow
+        .def("get_timers", &FDPF_XB_SparseLUSolver::get_timers, DocSolver::get_timers.c_str())  // returns the timers corresponding to times the solver spent in different part
+        .def("solve", &FDPF_XB_SparseLUSolver::compute_pf, py::call_guard<py::gil_scoped_release>(), DocSolver::compute_pf.c_str())  // perform the newton raphson optimization
+        .def("debug_get_Bp_python", &FDPF_XB_SparseLUSolver::debug_get_Bp_python, DocGridModel::_internal_do_not_use.c_str())  // perform the newton raphson optimization
+        .def("debug_get_Bpp_python", &FDPF_XB_SparseLUSolver::debug_get_Bpp_python, DocGridModel::_internal_do_not_use.c_str());  // perform the newton raphson optimization
+
+    py::class_<FDPF_BX_SparseLUSolver>(m, "FDPF_BX_SparseLUSolver", DocSolver::FDPF_BX_SparseLUSolver.c_str())
+        .def(py::init<>())
+        .def("get_Va", &FDPF_BX_SparseLUSolver::get_Va, DocSolver::get_Va.c_str())  // get the voltage angle vector (vector of double)
+        .def("get_Vm", &FDPF_BX_SparseLUSolver::get_Vm, DocSolver::get_Vm.c_str())  // get the voltage magnitude vector (vector of double)
+        .def("get_V", &FDPF_BX_SparseLUSolver::get_V, DocSolver::get_V.c_str()) 
+        .def("get_error", &FDPF_BX_SparseLUSolver::get_error, DocSolver::get_error.c_str())  // get the error message, see the definition of "err_" for more information
+        .def("get_nb_iter", &FDPF_BX_SparseLUSolver::get_nb_iter, DocSolver::get_nb_iter.c_str())  // return the number of iteration performed at the last optimization
+        .def("reset", &FDPF_BX_SparseLUSolver::reset, DocSolver::reset.c_str())  // reset the solver to its original state
+        .def("converged", &FDPF_BX_SparseLUSolver::converged, DocSolver::converged.c_str())  // whether the solver has converged
+        .def("compute_pf", &FDPF_BX_SparseLUSolver::compute_pf, py::call_guard<py::gil_scoped_release>(), DocSolver::compute_pf.c_str())  // compute the powerflow
+        .def("get_timers", &FDPF_BX_SparseLUSolver::get_timers, DocSolver::get_timers.c_str())  // returns the timers corresponding to times the solver spent in different part
+        .def("solve", &FDPF_BX_SparseLUSolver::compute_pf, py::call_guard<py::gil_scoped_release>(), DocSolver::compute_pf.c_str())  // perform the newton raphson optimization
+        .def("debug_get_Bp_python", &FDPF_BX_SparseLUSolver::debug_get_Bp_python, DocGridModel::_internal_do_not_use.c_str())  // perform the newton raphson optimization
+        .def("debug_get_Bpp_python", &FDPF_BX_SparseLUSolver::debug_get_Bpp_python, DocGridModel::_internal_do_not_use.c_str());  // perform the newton raphson optimization
+
     #if defined(KLU_SOLVER_AVAILABLE) || defined(_READ_THE_DOCS)
         py::class_<KLUSolver>(m, "KLUSolver", DocSolver::KLUSolver.c_str())
             .def(py::init<>())
@@ -136,6 +180,33 @@ PYBIND11_MODULE(lightsim2grid_cpp, m)
             .def("compute_pf", &KLUDCSolver::compute_pf, py::call_guard<py::gil_scoped_release>(), DocSolver::compute_pf.c_str())  // perform the newton raphson optimization
             .def("get_timers", &KLUDCSolver::get_timers, DocSolver::get_timers.c_str())  // returns the timers corresponding to times the solver spent in different part
             .def("solve", &KLUDCSolver::compute_pf, py::call_guard<py::gil_scoped_release>(), DocSolver::compute_pf.c_str());  // perform the newton raphson optimization
+        
+        py::class_<FDPF_XB_KLUSolver>(m, "FDPF_XB_KLUSolver", DocSolver::FDPF_XB_KLUSolver.c_str())
+            .def(py::init<>())
+            .def("get_Va", &FDPF_XB_KLUSolver::get_Va, DocSolver::get_Va.c_str())  // get the voltage angle vector (vector of double)
+            .def("get_Vm", &FDPF_XB_KLUSolver::get_Vm, DocSolver::get_Vm.c_str())  // get the voltage magnitude vector (vector of double)
+            .def("get_V", &FDPF_XB_KLUSolver::get_V, DocSolver::get_V.c_str()) 
+            .def("get_error", &FDPF_XB_KLUSolver::get_error, DocSolver::get_error.c_str())  // get the error message, see the definition of "err_" for more information
+            .def("get_nb_iter", &FDPF_XB_KLUSolver::get_nb_iter, DocSolver::get_nb_iter.c_str())  // return the number of iteration performed at the last optimization
+            .def("reset", &FDPF_XB_KLUSolver::reset, DocSolver::reset.c_str())  // reset the solver to its original state
+            .def("converged", &FDPF_XB_KLUSolver::converged, DocSolver::converged.c_str())  // whether the solver has converged
+            .def("compute_pf", &FDPF_XB_KLUSolver::compute_pf, py::call_guard<py::gil_scoped_release>(), DocSolver::compute_pf.c_str())  // compute the powerflow
+            .def("get_timers", &FDPF_XB_KLUSolver::get_timers, DocSolver::get_timers.c_str())  // returns the timers corresponding to times the solver spent in different part
+            .def("solve", &FDPF_XB_KLUSolver::compute_pf, py::call_guard<py::gil_scoped_release>(), DocSolver::compute_pf.c_str());  // perform the newton raphson optimization
+                
+        py::class_<FDPF_BX_KLUSolver>(m, "FDPF_BX_KLUSolver", DocSolver::FDPF_BX_KLUSolver.c_str())
+            .def(py::init<>())
+            .def("get_Va", &FDPF_BX_KLUSolver::get_Va, DocSolver::get_Va.c_str())  // get the voltage angle vector (vector of double)
+            .def("get_Vm", &FDPF_BX_KLUSolver::get_Vm, DocSolver::get_Vm.c_str())  // get the voltage magnitude vector (vector of double)
+            .def("get_V", &FDPF_BX_KLUSolver::get_V, DocSolver::get_V.c_str()) 
+            .def("get_error", &FDPF_BX_KLUSolver::get_error, DocSolver::get_error.c_str())  // get the error message, see the definition of "err_" for more information
+            .def("get_nb_iter", &FDPF_BX_KLUSolver::get_nb_iter, DocSolver::get_nb_iter.c_str())  // return the number of iteration performed at the last optimization
+            .def("reset", &FDPF_BX_KLUSolver::reset, DocSolver::reset.c_str())  // reset the solver to its original state
+            .def("converged", &FDPF_BX_KLUSolver::converged, DocSolver::converged.c_str())  // whether the solver has converged
+            .def("compute_pf", &FDPF_BX_KLUSolver::compute_pf, py::call_guard<py::gil_scoped_release>(), DocSolver::compute_pf.c_str())  // compute the powerflow
+            .def("get_timers", &FDPF_BX_KLUSolver::get_timers, DocSolver::get_timers.c_str())  // returns the timers corresponding to times the solver spent in different part
+            .def("solve", &FDPF_BX_KLUSolver::compute_pf, py::call_guard<py::gil_scoped_release>(), DocSolver::compute_pf.c_str());  // perform the newton raphson optimization
+                
     #endif  // KLU_SOLVER_AVAILABLE (or _READ_THE_DOCS)
 
     #if defined(NICSLU_SOLVER_AVAILABLE) || defined(_READ_THE_DOCS)
@@ -179,6 +250,33 @@ PYBIND11_MODULE(lightsim2grid_cpp, m)
             .def("compute_pf", &NICSLUDCSolver::compute_pf, py::call_guard<py::gil_scoped_release>(), DocSolver::compute_pf.c_str())  // perform the newton raphson optimization
             .def("get_timers", &NICSLUDCSolver::get_timers, DocSolver::get_timers.c_str())  // returns the timers corresponding to times the solver spent in different part
             .def("solve", &NICSLUDCSolver::compute_pf, py::call_guard<py::gil_scoped_release>(), DocSolver::compute_pf.c_str());  // perform the newton raphson optimization
+        
+        py::class_<FDPF_XB_NICSLUSolver>(m, "FDPF_XB_NICSLUSolver", DocSolver::FDPF_XB_NICSLUSolver.c_str())
+            .def(py::init<>())
+            .def("get_Va", &FDPF_XB_NICSLUSolver::get_Va, DocSolver::get_Va.c_str())  // get the voltage angle vector (vector of double)
+            .def("get_Vm", &FDPF_XB_NICSLUSolver::get_Vm, DocSolver::get_Vm.c_str())  // get the voltage magnitude vector (vector of double)
+            .def("get_V", &FDPF_XB_NICSLUSolver::get_V, DocSolver::get_V.c_str()) 
+            .def("get_error", &FDPF_XB_NICSLUSolver::get_error, DocSolver::get_error.c_str())  // get the error message, see the definition of "err_" for more information
+            .def("get_nb_iter", &FDPF_XB_NICSLUSolver::get_nb_iter, DocSolver::get_nb_iter.c_str())  // return the number of iteration performed at the last optimization
+            .def("reset", &FDPF_XB_NICSLUSolver::reset, DocSolver::reset.c_str())  // reset the solver to its original state
+            .def("converged", &FDPF_XB_NICSLUSolver::converged, DocSolver::converged.c_str())  // whether the solver has converged
+            .def("compute_pf", &FDPF_XB_NICSLUSolver::compute_pf, py::call_guard<py::gil_scoped_release>(), DocSolver::compute_pf.c_str())  // compute the powerflow
+            .def("get_timers", &FDPF_XB_NICSLUSolver::get_timers, DocSolver::get_timers.c_str())  // returns the timers corresponding to times the solver spent in different part
+            .def("solve", &FDPF_XB_NICSLUSolver::compute_pf, py::call_guard<py::gil_scoped_release>(), DocSolver::compute_pf.c_str());  // perform the newton raphson optimization
+        
+        py::class_<FDPF_BX_NICSLUSolver>(m, "FDPF_BX_NICSLUSolver", DocSolver::FDPF_BX_NICSLUSolver.c_str())
+            .def(py::init<>())
+            .def("get_Va", &FDPF_BX_NICSLUSolver::get_Va, DocSolver::get_Va.c_str())  // get the voltage angle vector (vector of double)
+            .def("get_Vm", &FDPF_BX_NICSLUSolver::get_Vm, DocSolver::get_Vm.c_str())  // get the voltage magnitude vector (vector of double)
+            .def("get_V", &FDPF_BX_NICSLUSolver::get_V, DocSolver::get_V.c_str()) 
+            .def("get_error", &FDPF_BX_NICSLUSolver::get_error, DocSolver::get_error.c_str())  // get the error message, see the definition of "err_" for more information
+            .def("get_nb_iter", &FDPF_BX_NICSLUSolver::get_nb_iter, DocSolver::get_nb_iter.c_str())  // return the number of iteration performed at the last optimization
+            .def("reset", &FDPF_BX_NICSLUSolver::reset, DocSolver::reset.c_str())  // reset the solver to its original state
+            .def("converged", &FDPF_BX_NICSLUSolver::converged, DocSolver::converged.c_str())  // whether the solver has converged
+            .def("compute_pf", &FDPF_BX_NICSLUSolver::compute_pf, py::call_guard<py::gil_scoped_release>(), DocSolver::compute_pf.c_str())  // compute the powerflow
+            .def("get_timers", &FDPF_BX_NICSLUSolver::get_timers, DocSolver::get_timers.c_str())  // returns the timers corresponding to times the solver spent in different part
+            .def("solve", &FDPF_BX_NICSLUSolver::compute_pf, py::call_guard<py::gil_scoped_release>(), DocSolver::compute_pf.c_str());  // perform the newton raphson optimization
+    
     #endif  // NICSLU_SOLVER_AVAILABLE (or _READ_THE_DOCS)
 
     #if defined(CKTSO_SOLVER_AVAILABLE) || defined(_READ_THE_DOCS)
@@ -222,6 +320,33 @@ PYBIND11_MODULE(lightsim2grid_cpp, m)
             .def("compute_pf", &CKTSODCSolver::compute_pf, py::call_guard<py::gil_scoped_release>(), DocSolver::compute_pf.c_str())  // perform the newton raphson optimization
             .def("get_timers", &CKTSODCSolver::get_timers, DocSolver::get_timers.c_str())  // returns the timers corresponding to times the solver spent in different part
             .def("solve", &CKTSODCSolver::compute_pf, py::call_guard<py::gil_scoped_release>(), DocSolver::compute_pf.c_str());  // perform the newton raphson optimization
+
+        py::class_<FDPF_XB_CKTSOSolver>(m, "FDPF_XB_CKTSOSolver", DocSolver::FDPF_XB_CKTSOSolver.c_str())
+            .def(py::init<>())
+            .def("get_Va", &FDPF_XB_CKTSOSolver::get_Va, DocSolver::get_Va.c_str())  // get the voltage angle vector (vector of double)
+            .def("get_Vm", &FDPF_XB_CKTSOSolver::get_Vm, DocSolver::get_Vm.c_str())  // get the voltage magnitude vector (vector of double)
+            .def("get_V", &FDPF_XB_CKTSOSolver::get_V, DocSolver::get_V.c_str()) 
+            .def("get_error", &FDPF_XB_CKTSOSolver::get_error, DocSolver::get_error.c_str())  // get the error message, see the definition of "err_" for more information
+            .def("get_nb_iter", &FDPF_XB_CKTSOSolver::get_nb_iter, DocSolver::get_nb_iter.c_str())  // return the number of iteration performed at the last optimization
+            .def("reset", &FDPF_XB_CKTSOSolver::reset, DocSolver::reset.c_str())  // reset the solver to its original state
+            .def("converged", &FDPF_XB_CKTSOSolver::converged, DocSolver::converged.c_str())  // whether the solver has converged
+            .def("compute_pf", &FDPF_XB_CKTSOSolver::compute_pf, py::call_guard<py::gil_scoped_release>(), DocSolver::compute_pf.c_str())  // compute the powerflow
+            .def("get_timers", &FDPF_XB_CKTSOSolver::get_timers, DocSolver::get_timers.c_str())  // returns the timers corresponding to times the solver spent in different part
+            .def("solve", &FDPF_XB_CKTSOSolver::compute_pf, py::call_guard<py::gil_scoped_release>(), DocSolver::compute_pf.c_str());  // perform the newton raphson optimization
+    
+        py::class_<FDPF_BX_CKTSOSolver>(m, "FDPF_BX_CKTSOSolver", DocSolver::FDPF_BX_CKTSOSolver.c_str())
+            .def(py::init<>())
+            .def("get_Va", &FDPF_BX_CKTSOSolver::get_Va, DocSolver::get_Va.c_str())  // get the voltage angle vector (vector of double)
+            .def("get_Vm", &FDPF_BX_CKTSOSolver::get_Vm, DocSolver::get_Vm.c_str())  // get the voltage magnitude vector (vector of double)
+            .def("get_V", &FDPF_BX_CKTSOSolver::get_V, DocSolver::get_V.c_str()) 
+            .def("get_error", &FDPF_BX_CKTSOSolver::get_error, DocSolver::get_error.c_str())  // get the error message, see the definition of "err_" for more information
+            .def("get_nb_iter", &FDPF_BX_CKTSOSolver::get_nb_iter, DocSolver::get_nb_iter.c_str())  // return the number of iteration performed at the last optimization
+            .def("reset", &FDPF_BX_CKTSOSolver::reset, DocSolver::reset.c_str())  // reset the solver to its original state
+            .def("converged", &FDPF_BX_CKTSOSolver::converged, DocSolver::converged.c_str())  // whether the solver has converged
+            .def("compute_pf", &FDPF_BX_CKTSOSolver::compute_pf, py::call_guard<py::gil_scoped_release>(), DocSolver::compute_pf.c_str())  // compute the powerflow
+            .def("get_timers", &FDPF_BX_CKTSOSolver::get_timers, DocSolver::get_timers.c_str())  // returns the timers corresponding to times the solver spent in different part
+            .def("solve", &FDPF_BX_CKTSOSolver::compute_pf, py::call_guard<py::gil_scoped_release>(), DocSolver::compute_pf.c_str());  // perform the newton raphson optimization
+    
     #endif  // CKTSO_SOLVER_AVAILABLE (or _READ_THE_DOCS)
 
     py::class_<GaussSeidelSolver>(m, "GaussSeidelSolver", DocSolver::GaussSeidelSolver.c_str())
@@ -266,7 +391,9 @@ PYBIND11_MODULE(lightsim2grid_cpp, m)
         .def("get_error", &ChooseSolver::get_error, DocSolver::get_V.c_str()) 
         .def("get_nb_iter", &ChooseSolver::get_nb_iter, DocSolver::get_nb_iter.c_str()) 
         .def("converged", &ChooseSolver::converged, DocSolver::converged.c_str()) 
-        .def("get_computation_time", &ChooseSolver::get_computation_time, DocSolver::get_computation_time.c_str());
+        .def("get_computation_time", &ChooseSolver::get_computation_time, DocSolver::get_computation_time.c_str())
+        .def("get_fdpf_xb_lu", &ChooseSolver::get_fdpf_xb_lu, py::return_value_policy::reference, DocGridModel::_internal_do_not_use.c_str())  // TODO this for all solver !
+        .def("get_fdpf_bx_lu", &ChooseSolver::get_fdpf_bx_lu, py::return_value_policy::reference, DocGridModel::_internal_do_not_use.c_str());
 
     // iterator for generators
     py::class_<DataGen>(m, "DataGen", DocIterator::DataGen.c_str())
@@ -402,7 +529,9 @@ PYBIND11_MODULE(lightsim2grid_cpp, m)
         .def_readonly("bus_ex_id", &DataLine::LineInfo::bus_ex_id, DocIterator::bus_ex_id.c_str())
         .def_readonly("r_pu", &DataLine::LineInfo::r_pu, DocIterator::r_pu.c_str())
         .def_readonly("x_pu", &DataLine::LineInfo::x_pu, DocIterator::x_pu.c_str())
-        .def_readonly("h_pu", &DataLine::LineInfo::h_pu, DocIterator::h_pu.c_str())
+        .def_readonly("h_pu", &DataLine::LineInfo::h_pu, DocIterator::x_pu.c_str())
+        .def_readonly("h_or_pu", &DataLine::LineInfo::h_or_pu, DocIterator::h_pu.c_str())
+        .def_readonly("h_ex_pu", &DataLine::LineInfo::h_ex_pu, DocIterator::h_pu.c_str())
         .def_readonly("has_res", &DataLine::LineInfo::has_res, DocIterator::has_res.c_str())
         .def_readonly("res_p_or_mw", &DataLine::LineInfo::res_p_or_mw, DocIterator::res_p_or_mw.c_str())
         .def_readonly("res_q_or_mvar", &DataLine::LineInfo::res_q_or_mvar, DocIterator::res_q_or_mvar.c_str())
@@ -457,6 +586,7 @@ PYBIND11_MODULE(lightsim2grid_cpp, m)
     py::class_<GridModel>(m, "GridModel", DocGridModel::GridModel.c_str())
         .def(py::init<>())
         .def("copy", &GridModel::copy)
+        .def_readwrite("_ls_to_pp", &GridModel::_ls_to_pp, "for converter from bus in lightsim2grid index to bus in pandapower index")
 
         // pickle
         .def(py::pickle(
@@ -499,6 +629,7 @@ PYBIND11_MODULE(lightsim2grid_cpp, m)
 
         // init its elements
         .def("init_powerlines", &GridModel::init_powerlines, DocGridModel::_internal_do_not_use.c_str())  // TODO code the possibility to add / remove a powerline after creation
+        .def("init_powerlines_full", &GridModel::init_powerlines_full, DocGridModel::_internal_do_not_use.c_str())  // TODO code the possibility to add / remove a powerline after creation
         .def("init_shunt", &GridModel::init_shunt, DocGridModel::_internal_do_not_use.c_str())  // same
         .def("init_trafo", &GridModel::init_trafo, DocGridModel::_internal_do_not_use.c_str())  // same 
         .def("init_generators", &GridModel::init_generators, DocGridModel::_internal_do_not_use.c_str())  // same
@@ -673,6 +804,10 @@ PYBIND11_MODULE(lightsim2grid_cpp, m)
         .def("set_trafo_hv_to_subid", &GridModel::set_trafo_hv_to_subid, DocGridModel::_internal_do_not_use.c_str())
         .def("set_trafo_lv_to_subid", &GridModel::set_trafo_lv_to_subid, DocGridModel::_internal_do_not_use.c_str())
         .def("set_storage_to_subid", &GridModel::set_storage_to_subid, DocGridModel::_internal_do_not_use.c_str())
+
+        // debug function (might disappear without further notice)
+        .def("debug_get_Bp_python", &GridModel::debug_get_Bp_python, DocGridModel::_internal_do_not_use.c_str())
+        .def("debug_get_Bpp_python", &GridModel::debug_get_Bpp_python, DocGridModel::_internal_do_not_use.c_str())
         ;
 
     py::class_<Computers>(m, "Computers", DocComputers::Computers.c_str())


Test patch:
diff --git a/lightsim2grid/tests/test_DoNothingACDC.py b/lightsim2grid/tests/test_DoNothingACDC.py
index ec412b7..86c5715 100644
--- a/lightsim2grid/tests/test_DoNothingACDC.py
+++ b/lightsim2grid/tests/test_DoNothingACDC.py
@@ -65,3 +65,7 @@ def test_dc(self):
     
     def test_ac(self):
         self._aux_test_do_nothing(is_dc=False)
+
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/lightsim2grid/tests/test_fdpf.py b/lightsim2grid/tests/test_fdpf.py
new file mode 100644
index 0000000..3dd8926
--- /dev/null
+++ b/lightsim2grid/tests/test_fdpf.py
@@ -0,0 +1,166 @@
+# Copyright (c) 2023, RTE (https://www.rte-france.com)
+# See AUTHORS.txt
+# This Source Code Form is subject to the terms of the Mozilla Public License, version 2.0.
+# If a copy of the Mozilla Public License, version 2.0 was not distributed with this file,
+# you can obtain one at http://mozilla.org/MPL/2.0/.
+# SPDX-License-Identifier: MPL-2.0
+# This file is part of LightSim2grid, LightSim2grid implements a c++ backend targeting the Grid2Op platform.
+
+import warnings
+import numpy as np
+from scipy.sparse import csc_matrix
+
+from lightsim2grid_cpp import FDPFMethod
+from lightsim2grid.gridmodel import init
+from lightsim2grid.solver import SolverType 
+
+import pandapower.networks as pn
+import pandapower as pp
+from pandapower.pypower.makeB import makeB
+from pandapower.pf.ppci_variables import _get_pf_variables_from_ppci
+from pandapower.pd2ppc import _pd2ppc
+from pandapower.auxiliary import _init_runpp_options
+
+import unittest
+
+
+
+class BaseFDPFTester:
+    def get_solving_method(self):
+        return FDPFMethod.XB
+    
+    def get_algo_pp(self):
+        # return "fdxb" if self.fdpf_meth == FDPFMethod.XB else "fdbx"  # in theory
+        return "fdbx" if self.fdpf_meth == FDPFMethod.XB else "fdxb"  # but.... https://github.com/e2nIEE/pandapower/issues/2142
+    
+    def get_network(self):
+        return pn.case14()
+        # self.net = pn.case118()
+        # self.net = pn.case300()
+    
+    def get_pp_options(self):
+        return dict(algorithm=self.get_algo_pp(),
+                    calculate_voltage_angles="auto",
+                    init="flat",
+                    max_iteration="auto",
+                    tolerance_mva=1e-8,
+                    trafo_model="t",
+                    trafo_loading="current",
+                    enforce_q_lims=False,
+                    check_connectivity=True,
+                    voltage_depend_loads=True, 
+                    consider_line_temperature=False,
+                    run_control=False,
+                    distributed_slack=False,
+                    tdpf=False,
+                    tdpf_delay_s=None)
+        
+    def setUp(self) -> None:
+        self.net = self.get_network()
+        with warnings.catch_warnings():
+            warnings.filterwarnings("ignore")
+            self.gridmodel = init(self.net)
+        self.tol = 1e-7
+        self.tol_solver = 1e-8
+        
+        # XB
+        self.fdpf_meth = self.get_solving_method()
+        self.alg = 2 if self.fdpf_meth == FDPFMethod.XB else 3
+        
+        # alg = 3  # BX  # TODO later
+        # self.fdpf_meth = FDPFMethod.BX  # TODO later
+        return super().setUp()
+    
+    def _aux_get_Bp_Bpp(self):
+        # init pp options
+        _init_runpp_options(self.net, **self.get_pp_options())
+        
+        # clear lookups
+        self.net._pd2ppc_lookups = {"bus": np.array([], dtype=np.int64), "ext_grid": np.array([], dtype=np.int64),
+                                    "gen": np.array([], dtype=np.int64), "branch": np.array([], dtype=np.int64)}
+
+        # convert pandapower net to ppc
+        ppc, ppci = _pd2ppc(self.net)
+        self.net["_ppc"] = ppc
+        baseMVA, bus, gen, branch, svc, tcsc, ref, pv, pq, *_, gbus, V0, ref_gens = _get_pf_variables_from_ppci(ppci)
+        pp_Bp, pp_Bpp = makeB(baseMVA, bus, np.real(branch), self.alg)
+        return pp_Bp, pp_Bpp, pv, pq, V0
+        
+    def test_bp_bpp(self):
+        """test that Bp and Bpp are correct (from the grid)"""
+        pp_Bp, pp_Bpp, *_= self._aux_get_Bp_Bpp()
+        
+        # now get the matrices 
+        self.gridmodel.ac_pf(1.04 * np.ones(self.net.bus.shape[0], dtype=complex), 10, 1e-7)  # need to init the "model bus id" id_grid_to_solver
+        ls_Bp = self.gridmodel.debug_get_Bp_python(self.fdpf_meth)
+        ls_Bpp = self.gridmodel.debug_get_Bpp_python(self.fdpf_meth)
+        
+        assert np.abs(pp_Bp - ls_Bp).max() <= self.tol, f"error in Bp: max {np.abs(pp_Bp - ls_Bp).max():.2e}"
+        assert np.abs(pp_Bpp - ls_Bpp).max() <= self.tol, f"error in Bpp: max {np.abs(pp_Bpp - ls_Bpp).max():.2e}"
+
+    def test_Bp_Bpp_solver(self):
+        """test that Bp and Bpp are correct (from the solver: where only some indices are used)"""
+        # retrieve the matrix in the solver
+        self.gridmodel.change_solver(SolverType.FDPF_XB_SparseLU if self.fdpf_meth == FDPFMethod.XB else SolverType.FDPF_BX_SparseLU)
+        V_ls = self.gridmodel.ac_pf(1.04 * np.ones(self.net.bus.shape[0], dtype=complex), 30, 1e0)  # to ensure it can "converge" to be able to retrieve Bp and Bpp
+        used_solver = self.gridmodel.get_solver().get_fdpf_xb_lu() if self.fdpf_meth == FDPFMethod.XB else self.gridmodel.get_solver().get_fdpf_bx_lu() 
+        ls_Bp = used_solver.debug_get_Bp_python()
+        ls_Bpp = used_solver.debug_get_Bpp_python()
+        # recreate the matrices from pandapower
+        grid_Bp, grid_Bpp, pv, pq, *_ = self._aux_get_Bp_Bpp()
+        pvpq = np.r_[pv, pq]
+        pp_Bp = grid_Bp[np.array([pvpq]).T, pvpq].tocsc()
+        pp_Bpp = grid_Bpp[np.array([pq]).T, pq].tocsc()
+        
+        # check they match
+        assert np.abs(pp_Bp - ls_Bp).max() <= self.tol, f"error in Bp: max {np.abs(pp_Bp - ls_Bp).max():.2f}"
+        assert np.abs(pp_Bpp - ls_Bpp).max() <= self.tol, f"error in Bpp: max {np.abs(pp_Bpp - ls_Bpp).max():.2f}"
+        
+    def test_solver(self):
+        self.gridmodel.change_solver(SolverType.FDPF_XB_SparseLU if self.fdpf_meth == FDPFMethod.XB else SolverType.FDPF_BX_SparseLU)
+        *_, V0 = self._aux_get_Bp_Bpp()
+        V_ls = self.gridmodel.ac_pf(1. * V0, 30, self.tol_solver * self.gridmodel.get_sn_mva())  # no division by sn_mva in pypower implementation TODO issue !
+        pp.runpp(self.net, algorithm=self.get_algo_pp(), tol=self.tol_solver, init="flat") 
+        assert V_ls.size > 0, f"lightsim2grid powerflow has diverged {self.gridmodel.get_solver().get_error()}"
+        nb_iter_ls = self.gridmodel.get_solver().get_nb_iter()
+        nb_iter_pp = self.net._ppc['iterations']
+        assert nb_iter_pp == nb_iter_ls, f"mismatch nb_iter {nb_iter_pp} vs {nb_iter_ls} maybe https://github.com/e2nIEE/pandapower/issues/2142 has been fixed ?"
+        V_pp = self.net.res_bus["vm_pu"].values * np.exp(1j * np.deg2rad(self.net.res_bus["va_degree"].values))
+        assert np.abs(V_ls - V_pp).max() <= self.tol, f"mismatch in V: {np.abs(V_ls - V_pp).max():.2e}"
+
+
+class FDPFTester_Case14_XB(BaseFDPFTester, unittest.TestCase):
+    pass
+
+
+class FDPFTester_Case14_BX(BaseFDPFTester, unittest.TestCase):
+    def get_solving_method(self):
+        return FDPFMethod.BX
+
+
+class FDPFTester_Case118_XB(BaseFDPFTester, unittest.TestCase):
+    def get_network(self):
+        return pn.case118()
+
+
+class FDPFTester_Case118_BX(BaseFDPFTester, unittest.TestCase):
+    def get_network(self):
+        return pn.case118()
+    def get_solving_method(self):
+        return FDPFMethod.BX
+
+
+class FDPFTester_Case300_XB(BaseFDPFTester, unittest.TestCase):
+    def get_network(self):
+        return pn.case300()
+
+
+class FDPFTester_Case300_BX(BaseFDPFTester, unittest.TestCase):
+    def get_network(self):
+        return pn.case300()
+    def get_solving_method(self):
+        return FDPFMethod.BX
+
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/lightsim2grid/tests/test_init_from_pypowsybl.py b/lightsim2grid/tests/test_init_from_pypowsybl.py
new file mode 100644
index 0000000..97528f5
--- /dev/null
+++ b/lightsim2grid/tests/test_init_from_pypowsybl.py
@@ -0,0 +1,299 @@
+# Copyright (c) 2023, RTE (https://www.rte-france.com)
+# See AUTHORS.txt
+# This Source Code Form is subject to the terms of the Mozilla Public License, version 2.0.
+# If a copy of the Mozilla Public License, version 2.0 was not distributed with this file,
+# you can obtain one at http://mozilla.org/MPL/2.0/.
+# SPDX-License-Identifier: MPL-2.0
+# This file is part of LightSim2grid, LightSim2grid implements a c++ backend targeting the Grid2Op platform.
+
+import pypowsybl as pp
+import pypowsybl.loadflow as lf
+import numpy as np
+import unittest
+import warnings
+
+from lightsim2grid.gridmodel.from_pypowsybl import init
+
+import pandapower.networks as pn
+import pandapower as pdp
+from lightsim2grid.gridmodel import init as init_from_pp
+
+
+class AuxInitFromPyPowSyBl:    
+    def get_pypo_grid(self):
+        return pp.network.create_ieee14()
+         
+    def get_gen_slack_id(self):
+        return 0
+    
+    def get_slackbus_id(self):
+        return 0
+    
+    def get_equiv_pdp_grid(self):
+        return pn.case14()
+    
+    def get_tol_eq(self):
+        return 1e-6
+    
+    def compare_pp(self):
+        """will this test suite compare pypowsybl and pandapower (cannot be used for ieee57 or ieee118)"""
+        return True
+    
+    def setUp(self) -> None:
+        self.network_ref = self.get_pypo_grid()
+        
+        # init equivalent pandapower grid (if any)
+        tmp = self.get_equiv_pdp_grid()
+        if tmp is not None:
+            self.pp_samecase = tmp
+            with warnings.catch_warnings():
+                warnings.filterwarnings("ignore")
+                self.ref_samecase = init_from_pp(self.pp_samecase)
+            self.can_pp = True
+        else:
+            # TODO will crash later if no equiv grid
+            self.can_pp = False
+            self.pp_samecase = None
+            self.ref_samecase = None
+            
+        # init lightsim2grid model
+        self.gridmodel = init(self.network_ref, gen_slack_id=self.get_gen_slack_id())
+        
+        # use some data
+        self.nb_bus_total = self.network_ref.get_buses().shape[0]
+        self.V_init_dc = np.ones(self.nb_bus_total, dtype=np.complex_)
+        self.V_init_ac = 1.04 * self.V_init_dc
+        self.tol = 1e-7  # for the solver
+        self.tol_eq = self.get_tol_eq()
+        return super().setUp()
+    
+    def test_basic(self):
+        """check that all elements are ok"""
+        assert len(self.gridmodel.get_lines()) == self.network_ref.get_lines().shape[0]    
+        assert len(self.gridmodel.get_loads()) == self.network_ref.get_loads().shape[0]    
+        assert len(self.gridmodel.get_generators()) == self.network_ref.get_generators().shape[0]    
+        assert len(self.gridmodel.get_trafos()) == self.network_ref.get_2_windings_transformers().shape[0]    
+        assert len(self.gridmodel.get_shunts()) == self.network_ref.get_shunt_compensators().shape[0]    
+    
+    def test_compare_pp(self):
+        """compare from the reference case14"""
+        if not self.can_pp:
+            self.skipTest("no equivalent pandapower grid")
+        if not self.compare_pp():
+            self.skipTest("It is expected that pandapower and pypowsybl gives different results") 
+            
+        # check the SBus and Ybus in DC (i need to run powerflow for that)      
+        v_ls = self.gridmodel.dc_pf(1.0 * self.V_init_dc, 2, self.tol)
+        v_ls_ref = self.ref_samecase.dc_pf(1.0 * self.V_init_dc, 2, self.tol)
+        slack_id = self.get_slackbus_id()
+        # (array([ 30,  31, 212, 218, 265]), array([265,  31, 212, 218,  30]))
+        
+        # for case 118
+        # bus_or = [64]
+        # bus_ex = [67]
+        # lines = [el.id for el in self.gridmodel.get_lines() if (el.bus_or_id in bus_or and el.bus_ex_id in bus_ex) or (el.bus_or_id in bus_ex and el.bus_ex_id in bus_or)]
+        # lines_ref = [el.id for el in self.ref_samecase.get_lines() if (el.bus_or_id in bus_or and el.bus_ex_id in bus_ex) or (el.bus_or_id in bus_ex and el.bus_ex_id in bus_or)]
+        # if not lines_ref:
+        #     lines_ref = [el.id for el in self.ref_samecase.get_trafos() if (el.bus_hv_id in bus_or and el.bus_lv_id in bus_ex) or (el.bus_hv_id in bus_ex and el.bus_lv_id in bus_or)]
+        # import pdb
+        # pdb.set_trace()
+        # # self.pp_samecase["_ppc"]["internal"]["Bbus"]
+        # self.pp_samecase["_ppc"]["internal"]["Ybus"][64,67]
+        # self.pp_samecase["_ppc"]["internal"]["bus"]
+
+        # for case 300
+        # np.where(tmp_.todense() >= 2000)
+        # (array([ 30,  31, 212, 218, 265]), array([265,  31, 212, 218,  30]))
+
+        # bus_or = [64]
+        # bus_ex = [67]
+        # lines = [el.id for el in self.gridmodel.get_lines() if (el.bus_or_id in bus_or and el.bus_ex_id in bus_ex) or (el.bus_or_id in bus_ex and el.bus_ex_id in bus_or)]
+        # lines_ref = [el.id for el in self.ref_samecase.get_lines() if (el.bus_or_id in bus_or and el.bus_ex_id in bus_ex) or (el.bus_or_id in bus_ex and el.bus_ex_id in bus_or)]
+        # if not lines_ref:
+        #     lines_ref = [el.id for el in self.ref_samecase.get_trafos() if (el.bus_hv_id in bus_or and el.bus_lv_id in bus_ex) or (el.bus_hv_id in bus_ex and el.bus_lv_id in bus_or)]
+        # # self.pp_samecase["_ppc"]["internal"]["Bbus"]
+        # self.pp_samecase["_ppc"]["internal"]["Ybus"][64,67]
+        # self.pp_samecase["_ppc"]["internal"]["bus"]
+                
+        assert np.abs(v_ls - v_ls_ref).max() <= self.tol_eq, "error for vresults for dc"
+        tmp_ = self.gridmodel.get_dcYbus() - self.ref_samecase.get_dcYbus()
+        assert np.abs(tmp_).max() <= self.tol_eq, "error for dcYbus"
+        # check Sbus without slack
+        if slack_id > 0:
+            assert np.abs(self.gridmodel.get_Sbus()[:slack_id] - self.ref_samecase.get_Sbus()[:slack_id]).max() <= self.tol_eq, "error for dc Sbus"
+        if slack_id != self.gridmodel.get_Sbus().shape[0] - 1:
+            assert np.abs(self.gridmodel.get_Sbus()[(slack_id+1):] - self.ref_samecase.get_Sbus()[(slack_id+1):]).max() <= self.tol_eq, "error for dc Sbus"
+
+        # same in AC
+        v_ls = self.gridmodel.ac_pf(self.V_init_ac, 2, self.tol)
+        v_ls_ref = self.ref_samecase.ac_pf(self.V_init_ac, 2, self.tol)
+        assert np.abs(self.gridmodel.get_Ybus() - self.ref_samecase.get_Ybus()).max() <= self.tol_eq, "error for Ybus"
+        # check Sbus without slack
+        if slack_id > 0:
+            assert np.abs(self.gridmodel.get_Sbus()[:slack_id] - self.ref_samecase.get_Sbus()[:slack_id]).max() <= self.tol_eq, "error for dc Sbus"
+        if slack_id != self.gridmodel.get_Sbus().shape[0] - 1:
+            assert np.abs(self.gridmodel.get_Sbus()[(slack_id+1):] - self.ref_samecase.get_Sbus()[(slack_id+1):]).max() <= self.tol_eq, "error for dc Sbus"
+
+    def test_dc_pf(self):
+        """test I get the same results as pandapower in dc"""
+        v_ls = self.gridmodel.dc_pf(self.V_init_dc, 2, self.tol)
+        if self.compare_pp():
+            v_ls_ref = self.ref_samecase.dc_pf(self.V_init_dc, 2, self.tol)
+            assert np.abs(v_ls - v_ls_ref).max() <= self.tol_eq, "error for vresults for dc"
+        lf.run_dc(self.network_ref)
+        if self.compare_pp():
+            pdp.rundcpp(self.pp_samecase)
+        
+        # for case 300
+        # np.where(np.abs(v_ang_ls - v_ang_pypo) > 3)  # => 173, 177
+        # bus_or = [173]
+        # bus_ex = [177]
+        # # lines = [el.id for el in self.gridmodel.get_lines() if (el.bus_or_id in bus_or and el.bus_ex_id in bus_ex) or (el.bus_or_id in bus_ex and el.bus_ex_id in bus_or)]
+        # lines = [el.id for el in self.gridmodel.get_lines() if (el.bus_or_id in bus_or or el.bus_or_id in bus_ex or el.bus_ex_id in bus_or or el.bus_ex_id in bus_ex)]
+        # trafos = [el.id for el in self.gridmodel.get_trafos() if (el.bus_hv_id in bus_or or el.bus_hv_id in bus_ex or el.bus_lv_id in bus_or or el.bus_lv_id in bus_ex)]
+        # shunts = [el.id for el in  self.gridmodel.get_shunts() if el.bus_id in bus_or or el.bus_id in bus_ex]
+        # trafo 77 => 204    2040    in cdf file
+        # trafo 85 => 
+        
+        # v_mag not really relevant in dc so i study only va
+        v_ang_pypo = self.network_ref.get_buses()["v_angle"].values
+        v_ang_ls = np.rad2deg(np.angle(v_ls))
+        if self.compare_pp():
+            v_ang_pp = self.pp_samecase.res_bus["va_degree"].values
+            assert np.abs(v_ang_ls - v_ang_pp).max() <= self.tol_eq, "error for va results for dc"
+        assert np.abs(v_ang_ls - v_ang_pypo).max() <= self.tol_eq
+        
+    def test_ac_pf(self):
+        # run the powerflows
+        v_ls = self.gridmodel.ac_pf(1.0 * self.V_init_ac, 10, self.tol)
+        if self.compare_pp():
+            v_ls_ref = self.ref_samecase.ac_pf(1.0 * self.V_init_ac, 10, self.tol)   
+            assert np.abs(v_ls - v_ls_ref).max() <= self.tol_eq, "error for vresults for ac"
+        
+        param = lf.Parameters(voltage_init_mode=pp._pypowsybl.VoltageInitMode.UNIFORM_VALUES,
+                              transformer_voltage_control_on=False,
+                              no_generator_reactive_limits=True,
+                              phase_shifter_regulation_on=False,
+                              simul_shunt=False,
+                              distributed_slack=False,
+                              provider_parameters={"slackBusSelectionMode": "NAME",
+                                                   "slackBusesIds": self.network_ref.get_buses().iloc[self.get_slackbus_id()].name}
+                              ) 
+        res_pypow = lf.run_ac(self.network_ref, parameters=param)
+        if self.compare_pp():
+            pdp.runpp(self.pp_samecase, init="flat")
+            
+        # check voltage angles
+        v_ang_pypo = self.network_ref.get_buses()["v_angle"].values
+        v_ang_ls = np.rad2deg(np.angle(v_ls))
+        if self.compare_pp():
+            v_ang_pp = self.pp_samecase.res_bus["va_degree"].values - self.pp_samecase.ext_grid["va_degree"].values
+            assert np.abs(v_ang_ls - v_ang_pp).max() <= self.tol_eq, "error for va results for ac"
+        if res_pypow[0].status == pp._pypowsybl.LoadFlowComponentStatus.CONVERGED:
+            assert np.abs(v_ang_ls - v_ang_pypo).max() <= self.tol_eq
+        
+        # check voltage magnitudes
+        bus_ref_kv = self.network_ref.get_voltage_levels().loc[self.network_ref.get_buses()["voltage_level_id"].values]["nominal_v"].values
+        v_mag_pypo = self.network_ref.get_buses()["v_mag"].values / bus_ref_kv
+        v_mag_ls = np.abs(v_ls)
+        if self.compare_pp():
+            v_mag_pp = self.pp_samecase.res_bus["vm_pu"].values
+            assert np.abs(v_mag_ls - v_mag_pp).max() <= self.tol_eq, "error for va results for dc"
+        if res_pypow[0].status == pp._pypowsybl.LoadFlowComponentStatus.CONVERGED:
+            assert np.abs(v_mag_ls - v_mag_pypo).max() <= self.tol_eq
+        
+        # check that pypow solution is "feasible" as seen by lightsim2grid
+        if res_pypow[0].status == pp._pypowsybl.LoadFlowComponentStatus.CONVERGED:
+            v_pypow = v_mag_pypo * np.exp(1j * np.deg2rad(v_ang_pypo))
+            v_pypow_ls = self.gridmodel.check_solution(v_pypow, False)
+            assert np.abs(v_pypow_ls).max() <= 10. * self.tol_eq
+     
+        
+class TestCase14FromPypo(AuxInitFromPyPowSyBl, unittest.TestCase):
+    pass
+        
+        
+class TestCase30FromPypo(AuxInitFromPyPowSyBl, unittest.TestCase):
+    # unittest.TestCase does not work because of https://github.com/powsybl/pypowsybl/issues/644
+    def get_pypo_grid(self):
+        return pp.network.create_ieee30()
+    
+    def get_equiv_pdp_grid(self):
+        return pn.case_ieee30()
+        
+        
+class TestCase57FromPypo(AuxInitFromPyPowSyBl, unittest.TestCase):
+    # does not appear to be the same grid !
+    def get_pypo_grid(self):
+        res = pp.network.create_ieee57()
+        # df = res.get_2_windings_transformers()[["rated_u1"]]
+        # df["rated_u1"] = 1.0
+        # res.update_2_windings_transformers(df)
+        
+        # df = res.get_lines()[["b1", "b2", "r"]]
+        # df["b1"] = 0.
+        # df["b2"] = 0.
+        # df["r"] = 0.
+        # res.update_lines(df)
+        return res
+    
+    def get_equiv_pdp_grid(self):
+        return pn.case57()
+    
+    def compare_pp(self):
+        return False
+    
+    def get_tol_eq(self):
+        return 1e-4  # otherwise vangle from pypowsybl and pandapower does not match
+                
+        
+class TestCase118FromPypo(AuxInitFromPyPowSyBl, unittest.TestCase):
+    # unittest.TestCase does not work because of https://github.com/powsybl/pypowsybl/issues/644
+    def get_pypo_grid(self):
+        return pp.network.create_ieee118()
+    
+    def get_equiv_pdp_grid(self):
+        return pn.case118()
+     
+    def get_gen_slack_id(self):
+        return 29
+    
+    def get_slackbus_id(self):
+        return 68
+    
+    def get_tol_eq(self):
+        # return 3e-3  # otherwise vangle from pypowsybl and pandapower does not match
+        return 3e-5  # otherwise vangle from pypowsybl and pandapower does not match
+    
+    def compare_pp(self):
+        return False
+                
+        
+class TestCase300FromPypo(AuxInitFromPyPowSyBl):
+    # does not work, probably grid not ordered the same
+    # need further investigation
+    def get_pypo_grid(self):
+        return pp.network.create_ieee300()
+    
+    def get_equiv_pdp_grid(self):
+        return pn.case300()
+    
+    def get_tol_eq(self):
+        # return 3e-3  # otherwise vangle from pypowsybl and pandapower does not match
+        return 3e-5  # otherwise vangle from pypowsybl and pandapower does not match
+    
+    def compare_pp(self):
+        return False
+     
+    def get_gen_slack_id(self):
+        # does not work with PP, probably bus not ordered the same
+        return 55
+    
+    def get_slackbus_id(self):
+        # does not work with PP, probably bus not ordered the same
+        return 257
+
+
+if __name__ == "__main__":
+    unittest.main()
+    
\ No newline at end of file
diff --git a/lightsim2grid/tests/test_multi_slack.py b/lightsim2grid/tests/test_multi_slack.py
new file mode 100644
index 0000000..6a98c65
--- /dev/null
+++ b/lightsim2grid/tests/test_multi_slack.py
@@ -0,0 +1,165 @@
+# Copyright (c) 2023, RTE (https://www.rte-france.com)
+# See AUTHORS.txt
+# This Source Code Form is subject to the terms of the Mozilla Public License, version 2.0.
+# If a copy of the Mozilla Public License, version 2.0 was not distributed with this file,
+# you can obtain one at http://mozilla.org/MPL/2.0/.
+# SPDX-License-Identifier: MPL-2.0
+# This file is part of LightSim2grid, LightSim2grid implements a c++ backend targeting the Grid2Op platform.
+
+import numpy as np
+import warnings
+
+import pandapower.networks
+import pandapower as pp
+
+from pandapower.pypower.makeBdc import makeBdc
+from pandapower.pypower.makeSbus import makeSbus
+from pandapower.pypower.idx_bus import GS
+from pandapower.pf.ppci_variables import _get_pf_variables_from_ppci
+from pandapower.pd2ppc import _pd2ppc
+from pandapower.auxiliary import _init_runpp_options
+
+import unittest
+
+from lightsim2grid.gridmodel import init
+from lightsim2grid.solver import SolverType
+
+
+class BaseMVOberrheinTester(unittest.TestCase):
+    def get_network(self):
+        res = pandapower.networks.mv_oberrhein()
+        res.switch["closed"] = True
+        res.trafo["shift_degree"] = 0.  # ignored in pandapower since `calculate_voltage_angles=False` 
+        # (otherwise powerflow diverge)
+        return res
+
+    def _aux_get_init_pp_data(self):
+        # init pp options
+        _init_runpp_options(self.net, **self.get_pp_options())
+        
+        # clear lookups
+        self.net._pd2ppc_lookups = {"bus": np.array([], dtype=np.int64), "ext_grid": np.array([], dtype=np.int64),
+                                    "gen": np.array([], dtype=np.int64), "branch": np.array([], dtype=np.int64)}
+
+        # convert pandapower net to ppc
+        ppc, ppci = _pd2ppc(self.net)
+        self.net["_ppc"] = ppc
+        baseMVA, bus, gen, branch, svc, tcsc, ref, pv, pq, *_, gbus, V0, ref_gens = _get_pf_variables_from_ppci(ppci)
+        return pv, pq, V0
+        
+    def get_pp_options(self):
+        return dict(algorithm="nr",
+                    calculate_voltage_angles=False,
+                    init="flat",
+                    max_iteration="auto",
+                    # max_iteration="auto",
+                    tolerance_mva=self.tol_solver,
+                    trafo_model="t",
+                    trafo_loading="current",
+                    enforce_q_lims=False,
+                    check_connectivity=True,
+                    voltage_depend_loads=True, 
+                    consider_line_temperature=False,
+                    run_control=False,
+                    distributed_slack=False,
+                    tdpf=False,
+                    tdpf_delay_s=None)
+        
+    def setUp(self) -> None:
+        self.net = self.get_network()
+        with warnings.catch_warnings():
+            warnings.filterwarnings("ignore")
+            self.gridmodel = init(self.net)
+        self.tol = 1e-7
+        self.tol_solver = 1e-8
+        
+    def test_Ybus(self):
+        # Ybus seen by pandapower
+        pp.runpp(self.net, **self.get_pp_options())
+        Ybus_pp_tmp = self.net["_ppc"]["internal"]["Ybus"]
+        id_bus_pp = self.net._pd2ppc_lookups["bus"][self.gridmodel._ls_to_pp]
+        assert (id_bus_pp == np.arange(id_bus_pp.size)).all()  # assert bus are sorted correctly and in the same order
+        # Ybus_pp = Ybus_pp_tmp[id_bus_pp, id_bus_pp]  # does not work in don't know why... anyway
+        Ybus_pp = Ybus_pp_tmp
+        *_, V0 = self._aux_get_init_pp_data()
+        # Ybus from lightsim2grid
+        V_ls = self.gridmodel.ac_pf(1. * V0, 30, self.tol_solver)
+        Ybus_ls = self.gridmodel.get_Ybus()
+        # np.where(np.abs(Ybus_pp - Ybus_ls).todense() >=100.)
+        assert np.abs((Ybus_pp - Ybus_ls)).max() <= self.tol, f"error in Ybus: max {np.abs((Ybus_pp - Ybus_ls).todense()).max():.2e}"
+        
+    def test_Sbus(self):
+        # Ybus seen by pandapower
+        pp.runpp(self.net, **self.get_pp_options())
+        bus_lookup = self.net._pd2ppc_lookups["bus"]  #[self.gridmodel._ls_to_pp]
+        Sbus_pp = self.net["_ppc"]["internal"]["Sbus"]
+        *_, V0 = self._aux_get_init_pp_data()
+        # Ybus from lightsim2grid
+        V_ls = self.gridmodel.ac_pf(1. * V0, 30, self.tol_solver)
+        Sbus_ls = self.gridmodel.get_Sbus()
+        slack_id = bus_lookup[self.net.ext_grid["bus"].values]
+        is_not_slack = np.ones(self.net.bus.shape[0], dtype=bool)
+        is_not_slack[slack_id] = False
+        # np.where((np.abs(Sbus_pp - Sbus_ls) >= 1e-3) & is_not_slack)
+        assert np.abs(Sbus_pp[is_not_slack] - Sbus_ls[is_not_slack]).max() <= self.tol, f"error in Sbus: max {np.abs((Sbus_pp - Sbus_ls)).max():.2e}"
+
+    def test_ac_solver(self):
+        *_, V0 = self._aux_get_init_pp_data()
+        self.gridmodel.change_solver(SolverType.SparseLUSingleSlack)
+        V_ls = self.gridmodel.ac_pf(1. * V0, 30, self.tol_solver)
+        pp.runpp(self.net, **self.get_pp_options()) 
+        assert V_ls.size > 0, f"lightsim2grid powerflow has diverged {self.gridmodel.get_solver().get_error()}"
+        nb_iter_ls = self.gridmodel.get_solver().get_nb_iter()
+        nb_iter_pp = self.net._ppc['iterations']
+        assert nb_iter_pp == nb_iter_ls, f"mismatch nb_iter {nb_iter_pp} vs {nb_iter_ls}"
+        V_pp = self.net.res_bus["vm_pu"].values * np.exp(1j * np.deg2rad(self.net.res_bus["va_degree"].values))
+        assert np.abs(V_ls - V_pp).max() <= self.tol, f"mismatch in V: {np.abs(V_ls - V_pp).max():.2e}"
+
+    def test_dc_Ybus(self):
+        # Ybus seen by pandapower
+        pp.rundcpp(self.net)
+        Ybus_pp = self.net["_ppc"]["internal"]["Bbus"]
+        *_, V0 = self._aux_get_init_pp_data()
+        # Ybus from lightsim2grid
+        V_ls = self.gridmodel.dc_pf(1. * V0, 30, self.tol_solver)
+        Ybus_ls = self.gridmodel.get_dcYbus()
+        # np.where(np.abs( (Ybus_pp - Ybus_ls).todense()) >= 1e-4)
+        assert np.abs( (Ybus_pp - Ybus_ls).todense()).max() <= self.tol, f"error in Ybus (dc): max {np.abs((Ybus_pp - Ybus_ls).todense()).max():.2e}"
+
+    def test_dc_Sbus(self):
+        self.skipTest("I don't really know what pandapower is doing")
+        # Ybus seen by pandapower
+        pp.rundcpp(self.net)
+        ppci = self.net._ppc
+        baseMVA = ppci['baseMVA']
+        bus = ppci['bus']
+        gen = ppci['gen']
+        branch = ppci['branch']
+        B, Bf, Pbusinj, Pfinj, Cft = makeBdc(bus, branch)
+        Sbus_pp = makeSbus(baseMVA, bus, gen) - Pbusinj - bus[:, GS] / baseMVA  # named Pbus in pandapower
+        self.net["_ppc"]["internal"]["Bbus"] * np.deg2rad(self.net.res_bus["va_degree"].values) - Sbus_pp
+        
+        *_, V0 = self._aux_get_init_pp_data()
+        # Ybus from lightsim2grid
+        V_ls = self.gridmodel.dc_pf(1. * V0, 30, self.tol_solver)
+        Sbus_ls = self.gridmodel.get_dcSbus()
+        slack_id = self.net.ext_grid["bus"].values
+        is_not_slack = np.ones(self.net.bus.shape[0], dtype=bool)
+        is_not_slack[slack_id] = False
+        # np.where((np.abs(Sbus_pp - Sbus_ls) >= 1e-3) & is_not_slack)
+        assert np.abs(Sbus_pp[is_not_slack] - Sbus_ls[is_not_slack]).max() <= self.tol, f"error in Sbus: max {np.abs((Sbus_pp - Sbus_ls)).max():.2e}"
+
+    def test_dc_solver(self):
+        self.skipTest("I don't really know what pandapower is doing")
+        *_, V0 = self._aux_get_init_pp_data()
+        V_ls = self.gridmodel.dc_pf(1. * V0, 30, self.tol_solver)
+        pp.rundcpp(self.net, tol=self.tol_solver, init="flat")
+        assert V_ls.size > 0, f"lightsim2grid powerflow has diverged {self.gridmodel.get_dc_solver().get_error()}"
+        V_pp = self.net.res_bus["vm_pu"].values * np.exp(1j * np.deg2rad(self.net.res_bus["va_degree"].values))
+        import pdb
+        pdb.set_trace()
+        assert np.abs(V_ls - V_pp).max() <= self.tol, f"mismatch in V (for DC): {np.abs(V_ls - V_pp).max():.2e}"
+        
+        
+if __name__ == "__main__":
+    unittest.main()
diff --git a/lightsim2grid/tests/test_nr_case300.py b/lightsim2grid/tests/test_nr_case300.py
new file mode 100644
index 0000000..cdec9f5
--- /dev/null
+++ b/lightsim2grid/tests/test_nr_case300.py
@@ -0,0 +1,166 @@
+# Copyright (c) 2023, RTE (https://www.rte-france.com)
+# See AUTHORS.txt
+# This Source Code Form is subject to the terms of the Mozilla Public License, version 2.0.
+# If a copy of the Mozilla Public License, version 2.0 was not distributed with this file,
+# you can obtain one at http://mozilla.org/MPL/2.0/.
+# SPDX-License-Identifier: MPL-2.0
+# This file is part of LightSim2grid, LightSim2grid implements a c++ backend targeting the Grid2Op platform.
+
+import warnings
+import numpy as np
+from scipy.sparse import csc_matrix
+
+from lightsim2grid.gridmodel import init
+from lightsim2grid.solver import SolverType 
+
+import pandapower.networks as pn
+import pandapower as pp
+from pandapower.pypower.makeBdc import makeBdc
+from pandapower.pypower.makeSbus import makeSbus
+from pandapower.pypower.idx_bus import GS
+from pandapower.pf.ppci_variables import _get_pf_variables_from_ppci
+from pandapower.pd2ppc import _pd2ppc
+from pandapower.auxiliary import _init_runpp_options
+
+import unittest
+
+
+
+class BaseCase300Tester(unittest.TestCase):
+    def get_network(self):
+        return pn.case300()
+
+    def _aux_get_init_pp_data(self):
+        # init pp options
+        _init_runpp_options(self.net, **self.get_pp_options())
+        
+        # clear lookups
+        self.net._pd2ppc_lookups = {"bus": np.array([], dtype=np.int64), "ext_grid": np.array([], dtype=np.int64),
+                                    "gen": np.array([], dtype=np.int64), "branch": np.array([], dtype=np.int64)}
+
+        # convert pandapower net to ppc
+        ppc, ppci = _pd2ppc(self.net)
+        self.net["_ppc"] = ppc
+        baseMVA, bus, gen, branch, svc, tcsc, ref, pv, pq, *_, gbus, V0, ref_gens = _get_pf_variables_from_ppci(ppci)
+        return pv, pq, V0
+        
+    def get_pp_options(self):
+        return dict(algorithm="nr",
+                    calculate_voltage_angles="auto",
+                    init="flat",
+                    max_iteration="auto",
+                    tolerance_mva=1e-8,
+                    trafo_model="t",
+                    trafo_loading="current",
+                    enforce_q_lims=False,
+                    check_connectivity=True,
+                    voltage_depend_loads=True, 
+                    consider_line_temperature=False,
+                    run_control=False,
+                    distributed_slack=False,
+                    tdpf=False,
+                    tdpf_delay_s=None)
+        
+    def setUp(self) -> None:
+        self.net = self.get_network()
+        with warnings.catch_warnings():
+            warnings.filterwarnings("ignore")
+            self.gridmodel = init(self.net)
+        self.tol = 1e-7
+        self.tol_solver = 1e-8
+        
+    def test_Ybus(self):
+        # Ybus seen by pandapower
+        pp.runpp(self.net)
+        Ybus_pp = self.net["_ppc"]["internal"]["Ybus"]
+        *_, V0 = self._aux_get_init_pp_data()
+        # Ybus from lightsim2grid
+        V_ls = self.gridmodel.ac_pf(1. * V0, 30, self.tol_solver)
+        Ybus_ls = self.gridmodel.get_Ybus()
+        assert np.abs( (Ybus_pp - Ybus_ls).todense()).max() <= self.tol, f"error in Ybus: max {np.abs((Ybus_pp - Ybus_ls).todense()).max():.2e}"
+        
+    def test_Sbus(self):
+        # Ybus seen by pandapower
+        pp.runpp(self.net)
+        Sbus_pp = self.net["_ppc"]["internal"]["Sbus"]
+        *_, V0 = self._aux_get_init_pp_data()
+        # Ybus from lightsim2grid
+        V_ls = self.gridmodel.ac_pf(1. * V0, 30, self.tol_solver)
+        Sbus_ls = self.gridmodel.get_Sbus()
+        slack_id = self.net.ext_grid["bus"].values
+        is_not_slack = np.ones(self.net.bus.shape[0], dtype=bool)
+        is_not_slack[slack_id] = False
+        # np.where((np.abs(Sbus_pp - Sbus_ls) >= 1e-3) & is_not_slack)
+        assert np.abs(Sbus_pp[is_not_slack] - Sbus_ls[is_not_slack]).max() <= self.tol, f"error in Sbus: max {np.abs((Sbus_pp - Sbus_ls)).max():.2e}"
+
+    def test_ac_solver(self):
+        *_, V0 = self._aux_get_init_pp_data()
+        V_ls = self.gridmodel.ac_pf(1. * V0, 30, self.tol_solver)
+        pp.runpp(self.net, algorithm="nr", tol=self.tol_solver, init="flat") 
+        assert V_ls.size > 0, f"lightsim2grid powerflow has diverged {self.gridmodel.get_solver().get_error()}"
+        nb_iter_ls = self.gridmodel.get_solver().get_nb_iter()
+        nb_iter_pp = self.net._ppc['iterations']
+        assert nb_iter_pp == nb_iter_ls, f"mismatch nb_iter {nb_iter_pp} vs {nb_iter_ls}"
+        V_pp = self.net.res_bus["vm_pu"].values * np.exp(1j * np.deg2rad(self.net.res_bus["va_degree"].values))
+        assert np.abs(V_ls - V_pp).max() <= self.tol, f"mismatch in V: {np.abs(V_ls - V_pp).max():.2e}"
+
+    def test_dc_Ybus(self):
+        # Ybus seen by pandapower
+        pp.rundcpp(self.net)
+        Ybus_pp = self.net["_ppc"]["internal"]["Bbus"]
+        *_, V0 = self._aux_get_init_pp_data()
+        # Ybus from lightsim2grid
+        V_ls = self.gridmodel.dc_pf(1. * V0, 30, self.tol_solver)
+        Ybus_ls = self.gridmodel.get_dcYbus()
+        # np.where(np.abs( (Ybus_pp - Ybus_ls).todense()) >= 1e-4)
+        assert np.abs( (Ybus_pp - Ybus_ls).todense()).max() <= self.tol, f"error in Ybus (dc): max {np.abs((Ybus_pp - Ybus_ls).todense()).max():.2e}"
+        
+    def test_Sbus(self):
+        # Ybus seen by pandapower
+        pp.runpp(self.net)
+        Sbus_pp = self.net["_ppc"]["internal"]["Sbus"]
+        *_, V0 = self._aux_get_init_pp_data()
+        
+        # Ybus from lightsim2grid
+        V_ls = self.gridmodel.ac_pf(1. * V0, 30, self.tol_solver)
+        Sbus_ls = self.gridmodel.get_Sbus()
+        slack_id = self.net.ext_grid["bus"].values
+        is_not_slack = np.ones(self.net.bus.shape[0], dtype=bool)
+        is_not_slack[slack_id] = False
+        # np.where((np.abs(Sbus_pp - Sbus_ls) >= 1e-3) & is_not_slack)
+        assert np.abs(Sbus_pp[is_not_slack] - Sbus_ls[is_not_slack]).max() <= self.tol, f"error in Sbus: max {np.abs((Sbus_pp - Sbus_ls)).max():.2e}"
+
+    def test_dc_Sbus(self):
+        self.skipTest("I don't really know what pandapower is doing")
+        # Ybus seen by pandapower
+        pp.rundcpp(self.net)
+        ppci = self.net._ppc
+        baseMVA = ppci['baseMVA']
+        bus = ppci['bus']
+        gen = ppci['gen']
+        branch = ppci['branch']
+        B, Bf, Pbusinj, Pfinj, Cft = makeBdc(bus, branch)
+        Sbus_pp = makeSbus(baseMVA, bus, gen) - Pbusinj - bus[:, GS] / baseMVA  # named Pbus in pandapower
+        self.net["_ppc"]["internal"]["Bbus"] * np.deg2rad(self.net.res_bus["va_degree"].values) - Sbus_pp
+        
+        *_, V0 = self._aux_get_init_pp_data()
+        # Ybus from lightsim2grid
+        V_ls = self.gridmodel.dc_pf(1. * V0, 30, self.tol_solver)
+        Sbus_ls = self.gridmodel.get_dcSbus()
+        slack_id = self.net.ext_grid["bus"].values
+        is_not_slack = np.ones(self.net.bus.shape[0], dtype=bool)
+        is_not_slack[slack_id] = False
+        # np.where((np.abs(Sbus_pp - Sbus_ls) >= 1e-3) & is_not_slack)
+        assert np.abs(Sbus_pp[is_not_slack] - Sbus_ls[is_not_slack]).max() <= self.tol, f"error in Sbus: max {np.abs((Sbus_pp - Sbus_ls)).max():.2e}"
+
+    def test_dc_solver(self):
+        *_, V0 = self._aux_get_init_pp_data()
+        V_ls = self.gridmodel.dc_pf(1. * V0, 30, self.tol_solver)
+        pp.rundcpp(self.net, tol=self.tol_solver, init="flat") 
+        assert V_ls.size > 0, f"lightsim2grid powerflow has diverged {self.gridmodel.get_solver().get_error()}"
+        V_pp = self.net.res_bus["vm_pu"].values * np.exp(1j * np.deg2rad(self.net.res_bus["va_degree"].values))
+        assert np.abs(V_ls - V_pp).max() <= self.tol, f"mismatch in V (for DC): {np.abs(V_ls - V_pp).max():.2e}"
+        
+        
+if __name__ == "__main__":
+    unittest.main()