fix the handling of phase shifters when init from pandapower

Signed-off-by: DONNOT Benjamin <benjamin.donnot@rte-france.com>

Author: BDonnot

CHANGELOG.rst

@@ -29,6 +29,12 @@ TODO: in `main.cpp` check the returned policy of pybind11 and also the `py::call
 TODO: a cpp class that is able to compute (DC powerflow) ContingencyAnalysis and TimeSeries using PTDF and LODF
 TODO: integration test with pandapower (see `pandapower/contingency/contingency.py` and import `lightsim2grid_installed` and check it's True)
 
+[0.12.1]  2026-xx-yy
+---------------------
+- [FIXED] phase shift transformers are now properly modeled
+  for both pandapower (new in this version) and pypowsybl (already
+  the case in previous version)
+
 [0.12.0] 2026-01-06
 --------------------
 - [BREAKING] for better consistency, and following pypowsybl convention, trafo and lines "side"

lightsim2grid/gridmodel/from_pypowsybl/_from_pypowsybl.py

@@ -361,7 +361,7 @@ def init(net : pypo.network.Network,
     ratio_tap_changer = net_pu.get_ratio_tap_changers()
 
     if 'alpha' in df_trafo_pu:
-        shift_ = np.rad2deg(df_trafo_pu['alpha'].values)  # given in radian by pypowsybl
+        shift_ = -np.rad2deg(df_trafo_pu['alpha'].values)  # given in radian by pypowsybl
     else:
         if net.get_phase_tap_changers().shape[0] > 0:
             raise RuntimeError("Phase tap changer are not handled by the pypowsybl converter "

lightsim2grid/tests/test_ACPF.py

@@ -60,7 +60,7 @@ def create_transformers(n_pdp):
     return r, x, g, b
 
 
-class MyTestCase(unittest.TestCase):
+class TestACPF(unittest.TestCase):
     def setUp(self) -> None:
         self.tol = 3e-5  # results are equal if they match up to tol
         self.raise_error = True
@@ -69,6 +69,31 @@ def test_case14(self):
         case = pn.case14()
         self._aux_test(case)
 
+    def test_case14_with_phaseshift(self):
+        case = pn.case14()        
+        hv_bus=0
+        lv_bus=2
+        pp.create_transformer_from_parameters(case,
+                                              hv_bus=hv_bus,
+                                              lv_bus=lv_bus,
+                                            #   sn_mva=1184.0,   # case RTE
+                                              sn_mva=9900.0,
+                                              vn_hv_kv=case.bus.iloc[hv_bus]["vn_kv"],
+                                              vn_lv_kv=case.bus.iloc[lv_bus]["vn_kv"],
+                                            #   i0_percent=-0.05152,   # case RTE
+                                              i0_percent=0.0,
+                                            #   vk_percent=0.404445,  # case RTE
+                                              vk_percent=2070.288000,
+                                            #   vkr_percent=0.049728, # case RTE
+                                              vkr_percent=0.0,
+                                              shift_degree=-10.0,  # case RTE
+                                            #   shift_degree=0.,
+                                              pfe_kw=0.
+                                              )
+        # self.tol = 2e-3
+        case.name = "case14_2"
+        self._aux_test(case)
+        
     def test_case39(self):
         case = pn.case39()
         self.tol = 3.1e-5  # results are equal if they match up to tol

lightsim2grid/tests/test_DCPF.py

@@ -43,11 +43,10 @@ def setUp(self) -> None:
     def test_case14(self):
         case = pn.case14()
         case.name = "case14"
-        # self.tol = 2e-3
         self._aux_test(case)
 
     def test_case14_with_phaseshift(self):
-        self.skipTest("pypowsybl and pandapower does not align on DC computation and phase shift... Lightsim2grid choses pypowsybl")
+        # self.skipTest("pandapower might do something weird in DC with phase shift")
         case = pn.case14()        
         hv_bus=0
         lv_bus=2
@@ -84,43 +83,43 @@ def test_case118(self):
         self.tol = 4e-5
         self._aux_test(case)
 
-    # def test_case1888rte(self):
-    #     # does not work probably with None in converters
-    #     case = pn.case1888rte()
-    #     case.name = "case1888rte"
-    #     self.tol_kcl = 2e-3
-    #     self._aux_test(case)
+    def test_case1888rte(self):
+        # self.skipTest("pandapower might do something weird in DC with phase shift")
+        case = pn.case1888rte()
+        case.name = "case1888rte"
+        self.tol_kcl = 2e-3
+        self._aux_test(case)
 
     def test_case300(self):
         case = pn.case300()
         self._aux_test(case)
 
     # def test_case2848rte(self):
-    #     # does not work probably with None in converters
+    #     # too large for CI
     #     case = pn.case2848rte()
     #     case.name = "case2848rte"
     #     self.tol = 4e-5
     #     self.tol_kcl = 7e-4
     #     self._aux_test(case)
 
     # def test_case6470rte(self):
-    #     # does not work probably with None in converters
+    # #     # too large for CI
     #     case = pn.case6470rte()
     #     case.name = "case6470rte"
     #     self.tol = 2e-4
     #     self.tol_kcl = 7e-4
     #     self._aux_test(case)
 
     # def test_case6495rte(self):
-    #     # does not work probably with None in converters
+    #     # too large for CI
     #     case = pn.case6495rte()
     #     case.name = "case6495rte"
     #     self.tol = 2e-4
     #     self.tol_kcl = 5e-3
     #     self._aux_test(case)
 
     # def test_case6515rte(self):
-    #     # does not work probably with None in converters
+    #     # too large for CI
     #     case = pn.case6515rte()
     #     case.name = "case6515rte"
     #     self.tol = 2e-4

lightsim2grid/tests/test_GridModel_pypowsybl.py

@@ -8,11 +8,8 @@
 
 import unittest
 import numpy as np
-import pandas as pd
-from packaging import version as version_packaging 
 
 import pypowsybl.network as pp_network
-import pypowsybl as pp
 import pypowsybl.loadflow as pp_lf
 
 from lightsim2grid.gridmodel import init_from_pypowsybl
@@ -535,6 +532,8 @@ def test_change_trafo_shift(self):
         # check it has an impact on the load flow
         Vfinal = self._run_both_pf(self.net_ref)
         with self.assertRaises(AssertionError):
+            # this checks the results are not the same
+            # when pypowsybl grid is not updated
             self.check_res(Vfinal, self.net_ref)
 
         self.skipTest("Hard to change the alpha of phase tap changer in pypowsbl")

lightsim2grid/tests/test_init_from_pypowsybl.py

@@ -112,7 +112,7 @@ def test_basic(self):
         assert len(self.gridmodel.get_shunts()) == self.network_ref.get_shunt_compensators().shape[0]    
 
     def test_compare_pp(self):
-        """compare from the reference case14"""
+        """test I get the same results than pypowsybl, in AC"""
         if not self.can_pp:
             self.skipTest("no equivalent pandapower grid")
         if not self.compare_pp():
@@ -158,7 +158,7 @@ def test_compare_pp(self):
             assert max_ <= self.tol_eq, f"error for dc Sbus : {max_:.2e}"
 
     def test_dc_pf(self):
-        """test I get the same results as pandapower in dc"""
+        """test I get the same results as pypowsybl in dc"""
         # if CURRENT_PYPOW_VERSION >= MIN_PYPO_DC_NOT_WORKING:
             # self.skipTest("Test not correct: pypowsybl change DC approx, see https://github.com/powsybl/pypowsybl/issues/1127")
         v_ls = self.gridmodel.dc_pf(self.V_init_dc, 2, self.tol)
@@ -175,7 +175,27 @@ def test_dc_pf(self):
         # 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))
-        # np.where(np.abs(v_ang_ls[reorder] - v_ang_pypo) >= 9.)
+        # for case300
+        # self.gridmodel.get_2_windings_transformers()[85]
+        # shift of 0.19896753472735357 radian
+        # 161 (bus1_id) and 424 (bus2_id)
+        # DC
+        # -9.859724243884386
+        # 9.859724243884614
+        # AC
+        # 83.57031489368106
+        # -83.56243958793598
+        
+        # self.network_ref.get_2_windings_transformers(all_attributes=True).iloc[85]
+        # alpha = 11.4 deg
+        # bus1_id VL196_0, bus2_id VL196_1
+        # self.network_ref.get_buses().loc[["VL196_0", "VL196_1"]]
+        # DC
+        # p1                         88.74429
+        # p2                        -88.74429
+        # AC
+        # p1                          83.570175
+        # p2                           -83.5623
         if self.compare_pp():
             v_ang_pp = self.pp_samecase.res_bus["va_degree"].values
             assert np.abs(v_ang_ls[reorder] - v_ang_pp).max() <= self.tol_eq, f"error for va results for dc: {np.abs(v_ang_ls[reorder] - v_ang_pp).max():.2e}"

src/GridModel.hpp

@@ -608,6 +608,14 @@ class GridModel final : public GenericContainer
         void change_ratio_trafo(int trafo_id, real_type new_ratio){
             trafos_.change_ratio(trafo_id, new_ratio, solver_control_);
         }
+
+        /**
+         * The shift is in radian (not degree !)
+         * 
+         * It is the shift on the "side 1" (regardless of the value of "is_tap_hv_side").
+         * If the tap is on the other side, the user has the reponsibility to
+         * take the opposite (ie -0.1 instead of +0.1)
+         */
         void change_shift_trafo(int trafo_id, real_type new_shift){
             trafos_.change_shift(trafo_id, new_shift, solver_control_);
         }

src/element_container/TrafoContainer.cpp

@@ -136,7 +136,7 @@ void TrafoContainer::_update_model_coeffs_one_el(int el_id)
     real_type theta_shift = shift_(el_id);
     if(!is_tap_hv_side_[el_id]){
         tau = my_one_ / tau;
-        theta_shift = -theta_shift;
+        // theta_shift = -theta_shift;
     }
     cplx_type eitheta_shift  = {my_one_, my_zero_};  // exp(j  * alpha)
     cplx_type emitheta_shift = {my_one_, my_zero_};  // exp(-j * alpha)
@@ -162,9 +162,9 @@ void TrafoContainer::_update_model_coeffs_one_el(int el_id)
     ydc_21_(el_id) = -tmp;
     ydc_12_(el_id) = -tmp;
 
-    dc_x_tau_shift_(el_id) = std::real(tmp) * theta_shift;
-    // if(!is_tap_hv_side_[el_id]) dc_x_tau_shift_(el_id) = std::real(tmp) * theta_shift;
-    // else dc_x_tau_shift_(el_id) = -std::real(tmp) * theta_shift;
+    // dc_x_tau_shift_(el_id) = std::real(tmp) * theta_shift;
+    if(!is_tap_hv_side_[el_id]) dc_x_tau_shift_(el_id) = -std::real(tmp) * theta_shift;
+    else dc_x_tau_shift_(el_id) = +std::real(tmp) * theta_shift;
 }
 
 void TrafoContainer::hack_Sbus_for_dc_phase_shifter(
@@ -222,8 +222,8 @@ void TrafoContainer::hack_Sbus_for_dc_phase_shifter(
             exc_ << " is connected (side 1) to a disconnected bus while being connected";
             throw std::runtime_error(exc_.str());
         }
-        Sbus.coeffRef(bus_id_solver_hv.cast_int()) += dc_x_tau_shift_[trafo_id];
-        Sbus.coeffRef(bus_id_solver_lv.cast_int()) -= dc_x_tau_shift_[trafo_id];
+        Sbus.coeffRef(bus_id_solver_hv.cast_int()) -= dc_x_tau_shift_[trafo_id];
+        Sbus.coeffRef(bus_id_solver_lv.cast_int()) += dc_x_tau_shift_[trafo_id];
     }
 }
 

src/element_container/TrafoContainer.hpp

@@ -150,12 +150,19 @@ class TrafoContainer : public TwoSidesContainer_rxh_A<OneSideContainer_ForBranch
                 }
         }
 
+        /**
+         * The shift is in radian (not degree !)
+         * 
+         * It is the shift on the "side 1" (regardless of the value of "is_tap_hv_side").
+         * If the tap is on the other side, the user has the reponsibility to
+         * take the opposite (ie -0.1 instead of +0.1)
+         */
         void change_shift(
             int el_id,
-            real_type new_shift,
+            real_type new_shift_rad,
             SolverControl & solver_control){
-                if(std::abs(shift_(el_id) - new_shift) >_tol_equal_float){
-                    shift_(el_id) = new_shift;
+                if(std::abs(shift_(el_id) - new_shift_rad) >_tol_equal_float){
+                    shift_(el_id) = new_shift_rad;
                     // TODO speed: only some part needs to be recomputed
                     _update_model_coeffs_one_el(el_id); 
                     solver_control.tell_recompute_ybus();

src/main.cpp

@@ -919,7 +919,7 @@ between 0 and `n_sub_ * max_nb_bus_per_sub_`
         .def("get_bus1_trafo", &GridModel::get_bus1_trafo, DocGridModel::_internal_do_not_use.c_str(), py::return_value_policy::reference)
         .def("get_bus2_trafo", &GridModel::get_bus2_trafo, DocGridModel::_internal_do_not_use.c_str(), py::return_value_policy::reference)
         .def("change_ratio_trafo", &GridModel::change_ratio_trafo, "TODO")
-        .def("change_shift_trafo", &GridModel::change_shift_trafo, "TODO")
+        .def("change_shift_trafo", &GridModel::change_shift_trafo, "TODO Change the phase shift ratio for a given transformer. It should be expressed in rad (not in deg) and on side1 (and not side2, as opposed to pypowsybl. If alpha is 10 in pypowsybl, you should convert it to rad (*eg* np.deg2rad(10.)), and invert it, so -np.deg2rad(10.). )")
 
         .def("deactivate_load", &GridModel::deactivate_load, DocGridModel::_internal_do_not_use.c_str())
         .def("reactivate_load", &GridModel::reactivate_load, DocGridModel::_internal_do_not_use.c_str())