fix an issue with the shunts

CHANGELOG.rst

@@ -18,7 +18,12 @@ Change Log
 - maybe have a look at suitesparse "sliplu" tools ?
 - easier building (get rid of the "make" part)
 
-[0.8.0] 2023-03-18
+[0.8.1] 2024-xx-yy
+--------------------
+- [FIXED] a bug with shunts when `nb_busbar_per_sub` >= 2
+- [IMPROVED] time measurments in python and c++
+
+[0.8.0] 2024-03-18
 --------------------
 - [BREAKING] now able to retrieve `dcSbus` with a dedicated method (and not with the old `get_Sbus`).
   If you previously used `gridmodel.get_Sbus()` to retrieve the Sbus used for DC powerflow, please use

benchmarks/benchmark_grid_size.py

@@ -46,7 +46,7 @@
               "case300.json",
               "case1354pegase.json",
               "case1888rte.json",
-              "GBnetwork.json",  # 2224 buses
+            #   "GBnetwork.json",  # 2224 buses
               "case2848rte.json",
               "case2869pegase.json",
               "case3120sp.json",

docs/conf.py

@@ -22,7 +22,7 @@
 author = 'Benjamin DONNOT'
 
 # The full version, including alpha/beta/rc tags
-release = "0.8.0"
+release = "0.8.1.dev0"
 version = '0.8'
 
 # -- General configuration ---------------------------------------------------

lightsim2grid/__init__.py

@@ -5,7 +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.
-__version__ = "0.8.0"
+__version__ = "0.8.1.dev0"
 
 __all__ = ["newtonpf", "SolverType", "ErrorType", "solver"]
 

lightsim2grid/lightSimBackend.py

@@ -130,6 +130,7 @@ def __init__(self,
         self._timer_preproc = 0.
         self._timer_postproc = 0.
         self._timer_solver = 0.
+        self._timer_read_data_back = 0.
 
         self.next_prod_p = None  # this vector is updated with the action that will modify the environment
         # it is done to keep track of the redispatching
@@ -209,6 +210,7 @@ def __init__(self,
         self._timer_postproc = 0.
         self._timer_preproc = 0.
         self._timer_solver = 0.
+        self._timer_read_data_back = 0.
 
         # hack for the storage unit:
         # in grid2op, for simplicity, I suppose that if a storage is alone on a busbar, and
@@ -265,6 +267,7 @@ def turnedoff_pv(self):
         self._grid.turnedoff_pv()
 
     def _fill_theta(self):
+        tick = time.perf_counter()
         # line_or_theta = np.empty(self.n_line)
         self.line_or_theta[:self.__nb_powerline] = self._grid.get_lineor_theta()
         self.line_or_theta[self.__nb_powerline:] = self._grid.get_trafohv_theta()
@@ -280,6 +283,7 @@ def _fill_theta(self):
 
         if self.__has_storage:
             self.storage_theta[:] = self._grid.get_storage_theta()
+        self._timer_read_data_back += time.perf_counter() - tick
 
     def get_theta(self):
         """
@@ -1043,7 +1047,8 @@ def runpf(self, is_dc=False):
 
                 self.timer_gridmodel_xx_pf += self._grid.timer_last_ac_pf
                 # timer_gridmodel_xx_pf takes all the time within the gridmodel "ac_pf"
-
+
+            beg_readback = time.perf_counter()
             self.V[:] = V
             (self.p_or[:self.__nb_powerline],
              self.q_or[:self.__nb_powerline],
@@ -1075,6 +1080,7 @@ def runpf(self, is_dc=False):
             if self.__has_storage:
                 self.storage_p[:], self.storage_q[:], self.storage_v[:] = self._grid.get_storages_res()
             self.storage_v[self.storage_v == -1.] = 0.  # voltage is 0. for disconnected elements in grid2op
+            self._timer_read_data_back +=  time.perf_counter() - beg_readback
 
             self.next_prod_p[:] = self.prod_p
 
@@ -1203,6 +1209,7 @@ def copy(self):
                            "_big_topo_to_obj", "dim_topo",
                            "_idx_hack_storage",
                            "_timer_preproc", "_timer_postproc", "_timer_solver",
+                           "_timer_read_data_back",
                            "supported_grid_format", 
                            "max_it", "tol", "_turned_off_pv", "_dist_slack_non_renew",
                            "_use_static_gen", "_loader_method", "_loader_kwargs",
@@ -1315,14 +1322,24 @@ def shunt_info(self):
         return self.cst_1 * self.sh_p, self.cst_1 * self.sh_q, self.cst_1 * self.sh_v, self.sh_bus
 
     def _set_shunt_info(self):
+        tick = time.perf_counter()
         self.sh_p[:], self.sh_q[:], self.sh_v[:]  = self._grid.get_shunts_res()
         shunt_bus = np.array([self._grid.get_bus_shunt(i) for i in range(type(self).n_shunt)], dtype=dt_int)
-        res_bus = np.ones(shunt_bus.shape[0], dtype=dt_int)  # by default all shunts are on bus one
-        res_bus[shunt_bus >= self.__nb_bus_before] = 2  # except the one that are connected to bus 2
-        res_bus[shunt_bus == -1] = -1  # or the one that are disconnected
-        self.sh_bus[:] = res_bus
+        cls = type(self)
+        if hasattr(cls, "global_bus_to_local"):
+            self.sh_bus[:] = cls.global_bus_to_local(shunt_bus, cls.shunt_to_subid)
+        else:
+            res = (1 * shunt_bus).astype(dt_int)  # make a copy
+            for i in range(cls.n_busbar_per_sub):
+                res[(i * cls.n_sub <= shunt_bus) & (shunt_bus < (i+1) * cls.n_sub)] = i + 1
+            res[shunt_bus == -1] = -1
+        # res_bus = np.ones(shunt_bus.shape[0], dtype=dt_int)  # by default all shunts are on bus one
+        # res_bus[shunt_bus >= self.__nb_bus_before] = 2  # except the one that are connected to bus 2
+        # res_bus[shunt_bus == -1] = -1  # or the one that are disconnected
+        # self.sh_bus[:] = res_bus
         self.sh_v[self.sh_v == -1.] = 0.  # in grid2op disco element have voltage of 0. and -1.
-
+        self._timer_read_data_back += time.perf_counter() - tick
+
     def _disconnect_line(self, id_):
         self.topo_vect[self.line_ex_pos_topo_vect[id_]] = -1
         self.topo_vect[self.line_or_pos_topo_vect[id_]] = -1
@@ -1346,4 +1363,5 @@ def reset(self, grid_path, grid_filename=None):
         self._timer_postproc = 0.
         self._timer_preproc = 0.
         self._timer_solver = 0.
+        self._timer_read_data_back = 0.
         self._grid.tell_solver_need_reset()

lightsim2grid/tests/test_n_busbar_per_sub.py

@@ -14,6 +14,7 @@
 
 from lightsim2grid import LightSimBackend
 
+
 class TestLightSimBackend_3busbars(unittest.TestCase):
     def get_nb_bus(self):
         return 3
@@ -147,7 +148,7 @@ def test_move_gen(self):
             self.env.backend.apply_action(bk_act)
             global_bus = sub_id + (new_bus -1) * cls.n_sub 
             if new_bus >= 1:
-                assert self.env.backend._grid.get_generators()[el_id].bus_id == global_bus, f"error for new_bus {new_bus}: {self.env.backend._grid.get_loads()[el_id].bus_id} vs {global_bus}"
+                assert self.env.backend._grid.get_generators()[el_id].bus_id == global_bus, f"error for new_bus {new_bus}: {self.env.backend._grid.get_generators()[el_id].bus_id} vs {global_bus}"
                 if line_or_id is not None:
                     assert self.env.backend._grid.get_lines()[line_or_id].bus_or_id == global_bus
                 else:
@@ -179,7 +180,7 @@ def test_move_storage(self):
             self.env.backend.apply_action(bk_act)
             global_bus = sub_id + (new_bus -1) * cls.n_sub 
             if new_bus >= 1:
-                assert self.env.backend._grid.get_storages()[el_id].bus_id == global_bus, f"error for new_bus {new_bus}: {self.env.backend._grid.get_loads()[el_id].bus_id} vs {global_bus}"
+                assert self.env.backend._grid.get_storages()[el_id].bus_id == global_bus, f"error for new_bus {new_bus}: {self.env.backend._grid.get_sotrages()[el_id].bus_id} vs {global_bus}"
                 if line_or_id is not None:
                     assert self.env.backend._grid.get_lines()[line_or_id].bus_or_id == global_bus
                 else:
@@ -243,9 +244,12 @@ def test_move_shunt(self):
             bk_act = self.env._backend_action_class()
             bk_act += act
             self.env.backend.apply_action(bk_act)
+            self.env.backend._set_shunt_info()
+            sh_p, sh_q, sh_v, sh_bus = self.env.backend.shunt_info()
+            assert sh_bus[el_id] == new_bus
             global_bus = sub_id + (new_bus -1) * cls.n_sub 
             if new_bus >= 1:
-                assert self.env.backend._grid.get_shunts()[el_id].bus_id == global_bus, f"error for new_bus {new_bus}: {self.env.backend._grid.get_loads()[el_id].bus_id} vs {global_bus}"
+                assert self.env.backend._grid.get_shunts()[el_id].bus_id == global_bus, f"error for new_bus {new_bus}: {self.env.backend._grid.get_shunts()[el_id].bus_id} vs {global_bus}"
                 if line_or_id is not None:
                     assert self.env.backend._grid.get_lines()[line_or_id].bus_or_id == global_bus
                 else:

setup.py

@@ -14,7 +14,7 @@
 from pybind11.setup_helpers import Pybind11Extension, build_ext
 
 
-__version__ = "0.8.0"
+__version__ = "0.8.1.dev0"
 KLU_SOLVER_AVAILABLE = False
 
 # Try to link against SuiteSparse (if available)

src/ChooseSolver.h

@@ -286,13 +286,26 @@ class ChooseSolver
             Eigen::SparseMatrix<real_type> res = get_J();
             return res;
         }
+
         double get_computation_time() const
         {
             auto p_solver = get_prt_solver("get_computation_time", true);
             const auto & res =  p_solver -> get_timers();
             return std::get<3>(res);
         }
 
+        std::tuple<double, double, double, double> get_timers() const
+        {
+            auto p_solver = get_prt_solver("get_timers", true);
+            return p_solver -> get_timers();
+        }
+
+        TimerJacType get_timers_jacobian() const
+        {
+            const BaseAlgo * p_solver = get_prt_solver("get_timers_jacobian", true);
+            return p_solver -> get_timers_jacobian();
+        }
+
         ErrorType get_error() const{
             auto p_solver = get_prt_solver("get_error", true);
             return p_solver -> get_error();

src/main.cpp

@@ -393,6 +393,8 @@ PYBIND11_MODULE(lightsim2grid_cpp, m)
         .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_timers", &ChooseSolver::get_timers, "TODO")
+        .def("get_timers_jacobian", &ChooseSolver::get_timers_jacobian, "TODO")
         .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());
 

src/powerflow_algorithm/BaseAlgo.h

@@ -30,6 +30,9 @@
 
 class GridModel;
 
+typedef std::tuple<double, double, double, double, 
+                   double, double, double, double, 
+                   double> TimerJacType;
 
 /**
 This class represents a algorithm to compute powerflow.
@@ -85,6 +88,22 @@ class BaseAlgo : public BaseConstants
               timer_Fx_, timer_solve_, timer_check_, timer_total_nr_);
             return res;
         }
+
+        virtual TimerJacType  get_timers_jacobian() const
+        {
+            TimerJacType res = {
+                timer_Fx_,
+                timer_solve_,
+                -1.,  // not available for non NR solver, so I put -1
+                timer_check_,
+                -1.,  // not available for non NR solver, so I put -1
+                -1.,  // not available for non NR solver, so I put -1
+                -1.,  // not available for non NR solver, so I put -1
+                -1.,  // not available for non NR solver, so I put -1
+                timer_total_nr_
+            };
+            return res;
+        }
 
         virtual
         bool compute_pf(const Eigen::SparseMatrix<cplx_type> & Ybus,

src/powerflow_algorithm/BaseNRAlgo.h

@@ -18,7 +18,14 @@ template<class LinearSolver>
 class BaseNRAlgo : public BaseAlgo
 {
     public:
-        BaseNRAlgo():BaseAlgo(true), need_factorize_(true), timer_initialize_(0.), timer_dSbus_(0.), timer_fillJ_(0.) {}
+        BaseNRAlgo():
+            BaseAlgo(true),
+            need_factorize_(true),
+            timer_initialize_(0.),
+            timer_dSbus_(0.),
+            timer_fillJ_(0.),
+            timer_Va_Vm_(0.),
+            timer_pre_proc_(0.){}
 
         virtual
         Eigen::Ref<const Eigen::SparseMatrix<real_type> > get_J() const {
@@ -32,11 +39,18 @@ class BaseNRAlgo : public BaseAlgo
         }
 
         virtual
-        std::tuple<double, double, double, double, double, double, double> get_timers_jacobian()
+        TimerJacType get_timers_jacobian() const
         {
             // TODO refacto that, and change the order
-            auto res = std::tuple<double, double, double, double, double, double, double>(
-              timer_Fx_, timer_solve_, timer_initialize_, timer_check_, timer_dSbus_, timer_fillJ_, timer_total_nr_);
+            auto res = TimerJacType(timer_Fx_,
+                                    timer_solve_,
+                                    timer_initialize_,
+                                    timer_check_,
+                                    timer_dSbus_,
+                                    timer_fillJ_,
+                                    timer_Va_Vm_,
+                                    timer_pre_proc_,
+                                    timer_total_nr_);
             return res;
         }
 
@@ -59,6 +73,8 @@ class BaseNRAlgo : public BaseAlgo
             BaseAlgo::reset_timer();
             timer_dSbus_ = 0.;
             timer_fillJ_ = 0.;
+            timer_Va_Vm_ = 0.;
+            timer_pre_proc_ = 0.;
             timer_initialize_ = 0.;
         }
         virtual
@@ -171,6 +187,8 @@ class BaseNRAlgo : public BaseAlgo
         double timer_initialize_;
         double timer_dSbus_;
         double timer_fillJ_;
+        double timer_Va_Vm_;
+        double timer_pre_proc_;
 
 
     Eigen::SparseMatrix<real_type>

src/powerflow_algorithm/BaseNRAlgo.tpp

@@ -13,15 +13,15 @@
 
 template<class LinearSolver>
 bool BaseNRAlgo<LinearSolver>::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
-                                            )
+                                          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
@@ -54,6 +54,7 @@ bool BaseNRAlgo<LinearSolver>::compute_pf(const Eigen::SparseMatrix<cplx_type> &
     reset_if_needed();
     err_ = ErrorType::NoError;  // reset the error if previous error happened
     auto timer = CustTimer();
+    auto timer_pre_proc = CustTimer();
 
     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
@@ -77,6 +78,7 @@ bool BaseNRAlgo<LinearSolver>::compute_pf(const Eigen::SparseMatrix<cplx_type> &
     V_ = V;
     Vm_ = V_.array().abs();  // update Vm and Va again in case
     Va_ = V_.array().arg();  // we "wrapped around" with a negative Vm
+    timer_pre_proc_ += timer_pre_proc.duration();
 
     // first check, if the problem is already solved, i stop there
     // compute a first time the mismatch to initialize the slack bus
@@ -132,20 +134,29 @@ bool BaseNRAlgo<LinearSolver>::compute_pf(const Eigen::SparseMatrix<cplx_type> &
         }
         // const auto dx = -F;  // removed for speed optimization (-= used below)
 
+        auto timer_va_vm = CustTimer();
+        slack_absorbed -= F(0); // by convention in fill_jacobian_matrix the slack bus is the first component
         // update voltage (this should be done consistently with "_evaluate_Fx")
         if (n_pv > 0) Va_(my_pv) -= F.segment(1, n_pv);
         if (n_pq > 0){
             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 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;
         // TODO change here for not having to cast all the time ... maybe
         V_ = Vm_.array() * (Va_.array().cos().template cast<cplx_type>() + my_i * Va_.array().sin().template cast<cplx_type>() );
-        Vm_ = V_.array().abs();  // update Vm and Va again in case
-        Va_ = V_.array().arg();  // we wrapped around with a negative Vm TODO more efficient way maybe ?
+        // V_ = Vm_.array() *  (my_i * Va_.array().template cast<cplx_type>()).exp() ;
+        if(Vm_.minCoeff() < 0.)
+        {
+            // update Vm and Va again in case
+            // we wrapped around with a negative Vm TODO more efficient way maybe ?
+            Vm_ = V_.array().abs();  
+            Va_ = V_.array().arg();  
+        }
+        timer_Va_Vm_ += timer_va_vm.duration();
 
         F = _evaluate_Fx(Ybus, V_, Sbus, slack_bus_id, slack_absorbed, slack_weights, my_pv, pq);
         bool tmp = F.allFinite();
@@ -345,10 +356,9 @@ void BaseNRAlgo<LinearSolver>::fill_jacobian_matrix(const Eigen::SparseMatrix<cp
     `slack` is the representation of the equation connecting together the slack buses (represented by slack_weights)
     the remaining pq components are all 0.
     **/
-
-    auto timer = CustTimer();
     _dSbus_dV(Ybus, V);
 
+    auto timer = CustTimer();
     const auto n_pvpq = pvpq.size();
     const auto n_pq = pq.size();
     const auto size_j = n_pvpq + n_pq + 1;   // +1 because i add the slack bus