python calls updated, now need to solve the segfault [skip ci]

lightsim2grid/lightSimBackend.py

@@ -735,7 +735,7 @@ def _aux_finish_setup_after_reading(self):
             self.storage_theta = np.full(cls.n_storage, dtype=dt_float, fill_value=np.NaN).reshape(-1)
 
         self._count_object_per_bus()
-        self._grid.tell_topo_changed()
+        self._grid.tell_solver_need_reset()
         self.__me_at_init = self._grid.copy()
         self.__init_topo_vect = np.ones(cls.dim_topo, dtype=dt_int)
         self.__init_topo_vect[:] = self.topo_vect
@@ -894,7 +894,9 @@ def runpf(self, is_dc=False):
                     self._grid.deactivate_result_computation()
                     # if I init with dc values, it should depends on previous state
                     self.V[:] = self._grid.get_init_vm_pu()  # see issue 30
+                    print("before dc pf")
                     Vdc = self._grid.dc_pf(copy.deepcopy(self.V), self.max_it, self.tol)
+                    print("after dc pf")
                     self._grid.reactivate_result_computation()
                     if Vdc.shape[0] == 0:
                         raise BackendError(f"Divergence of DC powerflow (non connected grid) at the initialization of AC powerflow. Detailed error: {self._grid.get_dc_solver().get_error()}")
@@ -903,6 +905,7 @@ def runpf(self, is_dc=False):
                     V_init = copy.deepcopy(self.V)
                 tick = time.perf_counter()
                 self._timer_preproc += tick - beg_preproc
+                print("before ac pf")
                 V = self._grid.ac_pf(V_init, self.max_it, self.tol)
                 self._timer_solver += time.perf_counter() - tick
                 if V.shape[0] == 0:
@@ -968,11 +971,11 @@ def runpf(self, is_dc=False):
             if (self.line_or_theta >= 1e6).any() or (self.line_ex_theta >= 1e6).any():
                 raise BackendError(f"Some theta are above 1e6 which should not be happening !")
             res = True
-            self._grid.unset_topo_changed()
+            self._grid.unset_changes()
             self._timer_postproc += time.perf_counter() - beg_postroc
         except Exception as exc_:
             # of the powerflow has not converged, results are Nan
-            self._grid.tell_topo_changed()
+            self._grid.unset_changes()
             self._fill_nans()
             res = False
             my_exc_ = exc_
@@ -1190,7 +1193,7 @@ def get_current_solver_type(self):
     def reset(self, grid_path, grid_filename=None):
         self._fill_nans()
         self._grid = self.__me_at_init.copy()
-        self._grid.tell_topo_changed()
+        self._grid.unset_changes()
         self._grid.change_solver(self.__current_solver_type)
         self._handle_turnedoff_pv()
         self.topo_vect[:] = self.__init_topo_vect

src/BaseSolver.cpp

@@ -23,6 +23,7 @@ void BaseSolver::reset(){
     err_ = ErrorType::NotInitError; //error message:
 
     _solver_control = SolverControl();
+    _solver_control.tell_all_changed();
 }
 
 

src/ChooseSolver.cpp

@@ -7,3 +7,80 @@
 // This file is part of LightSim2grid, LightSim2grid implements a c++ backend targeting the Grid2Op platform.
 
 #include "ChooseSolver.h"
+
+std::ostream& operator<<(std::ostream& out, const SolverType& solver_type)
+{
+    switch (solver_type)
+    {
+    case SolverType::SparseLU:
+        out << "SparseLU";
+        break;
+    case SolverType::KLU:
+        out << "KLU";
+        break;
+    case SolverType::GaussSeidel:
+        out << "GaussSeidel";
+        break;
+    case SolverType::DC:
+        out << "DC";
+        break;
+    case SolverType::GaussSeidelSynch:
+        out << "GaussSeidelSynch";
+        break;
+    case SolverType::NICSLU:
+        out << "NICSLU";
+        break;
+    case SolverType::SparseLUSingleSlack:
+        out << "SparseLUSingleSlack";
+        break;
+    case SolverType::KLUSingleSlack:
+        out << "KLUSingleSlack";
+        break;
+    case SolverType::NICSLUSingleSlack:
+        out << "NICSLUSingleSlack";
+        break;
+    case SolverType::KLUDC:
+        out << "KLUDC";
+        break;
+    case SolverType::NICSLUDC:
+        out << "NICSLUDC";
+        break;
+    case SolverType::CKTSO:
+        out << "CKTSO";
+        break;
+    case SolverType::CKTSOSingleSlack:
+        out << "CKTSOSingleSlack";
+        break;
+    case SolverType::CKTSODC:
+        out << "CKTSODC";
+        break;
+    case SolverType::FDPF_XB_SparseLU:
+        out << "FDPF_XB_SparseLU";
+        break;
+    case SolverType::FDPF_BX_SparseLU:
+        out << "FDPF_BX_SparseLU";
+        break;
+    case SolverType::FDPF_XB_KLU:
+        out << "FDPF_XB_KLU";
+        break;
+    case SolverType::FDPF_BX_KLU:
+        out << "FDPF_BX_KLU";
+        break;
+    case SolverType::FDPF_XB_NICSLU:
+        out << "FDPF_XB_NICSLU";
+        break;
+    case SolverType::FDPF_BX_NICSLU:
+        out << "FDPF_BX_NICSLU";
+        break;
+    case SolverType::FDPF_XB_CKTSO:
+        out << "FDPF_XB_CKTSO";
+        break;
+    case SolverType::FDPF_BX_CKTSO:
+        out << "FDPF_BX_CKTSO";
+        break;
+    default:
+        out << "(unknown)";
+        break;
+    }
+   return out;
+}

src/ChooseSolver.h

@@ -27,6 +27,8 @@ enum class SolverType {SparseLU, KLU, GaussSeidel, DC, GaussSeidelSynch, NICSLU,
                        FDPF_XB_CKTSO,  FDPF_BX_CKTSO  // from 0.7.5
                        };
 
+
+std::ostream& operator<<(std::ostream& out, const SolverType& solver_type);
 // TODO define a template class instead of these weird stuff !!!
 // TODO export all methods from base class !
 
@@ -278,7 +280,7 @@ class ChooseSolver
         }
 
         void tell_solver_control(const SolverControl & solver_control){
-            auto p_solver = get_prt_solver("get_ptdf", true);
+            auto p_solver = get_prt_solver("tell_solver_control", false);
             p_solver -> tell_solver_control(solver_control);
         }
 
@@ -312,7 +314,17 @@ class ChooseSolver
     private:
         void check_right_solver(const std::string & error_msg) const
         {
-            if(_solver_type != _type_used_for_nr) throw std::runtime_error("ChooseSolver: Solver mismatch when calling '"+error_msg+"': current solver is not the last solver used to perform a powerflow");
+            if(_solver_type != _type_used_for_nr){
+                std::ostringstream exc_;
+                exc_ << "ChooseSolver: Solver mismatch when calling '";
+                exc_ << error_msg;
+                exc_ << ": current solver (";
+                exc_ << _solver_type;
+                exc_ << ") is not the one used to perform a powerflow (";
+                exc_ << _type_used_for_nr;
+                exc_ << ").";
+                throw std::runtime_error(exc_.str());
+            }
 
             #ifndef KLU_SOLVER_AVAILABLE
                 if(_solver_type == SolverType::KLU){

src/DataDCLine.cpp

@@ -87,11 +87,9 @@ void DataDCLine::disconnect_if_not_in_main_component(std::vector<bool> & busbar_
         auto bus_or = bus_or_id(i);
         auto bus_ex = bus_ex_id(i);
         if(!busbar_in_main_component[bus_or]) {
-            bool tmp = false;
             from_gen_.deactivate(i, unused_solver_control);
         }
         if(!busbar_in_main_component[bus_ex]) {
-            bool tmp = false;
             to_gen_.deactivate(i, unused_solver_control);
         }
         // if(!busbar_in_main_component[bus_or] || !busbar_in_main_component[bus_ex]){

src/GridModel.cpp

@@ -414,25 +414,67 @@ CplxVect GridModel::pre_process_solver(const CplxVect & Vinit,
                                        const SolverControl & solver_control)
 {
     // TODO get rid of the "is_ac" argument: this info is available in the _solver already
-    if(is_ac) _solver.tell_solver_control(solver_control);
-    else _dc_solver.tell_solver_control(solver_control);
-
-    if (solver_control.has_slack_participate_changed()) slack_bus_id_ = generators_.get_slack_bus_id();
-    if (solver_control.ybus_change_sparsity_pattern() || solver_control.has_dimension_changed()) init_Ybus(Ybus, id_me_to_solver, id_solver_to_me);
-    if (solver_control.ybus_change_sparsity_pattern() || solver_control.has_dimension_changed() || solver_control.need_recompute_ybus()) fillYbus(Ybus, is_ac, id_me_to_solver);
-    if (solver_control.has_dimension_changed()) init_Sbus(Sbus_, id_me_to_solver, id_solver_to_me, slack_bus_id_solver);
-    if (solver_control.has_slack_participate_changed() || solver_control.has_pv_changed() || solver_control.has_pq_changed()) fillpv_pq(id_me_to_solver, id_solver_to_me, slack_bus_id_solver, solver_control);
+    if(is_ac){
+        _solver.tell_solver_control(solver_control);
+        if(solver_control.need_reset_solver()) _solver.reset();
+    } else {
+        _dc_solver.tell_solver_control(solver_control);
+        if(solver_control.need_reset_solver()){
+            std::cout << "_dc_solver.reset();" << std::endl;
+            _dc_solver.reset();
+        }
+    }
+
+    if (solver_control.need_reset_solver() || 
+        solver_control.has_slack_participate_changed()){
+            std::cout << "get_slack_bus_id;" << std::endl;
+            slack_bus_id_ = generators_.get_slack_bus_id();
+        }
+    if (solver_control.need_reset_solver() ||
+        solver_control.ybus_change_sparsity_pattern() || 
+        solver_control.has_dimension_changed()){
+            init_Ybus(Ybus, id_me_to_solver, id_solver_to_me);
+            std::cout << "init_Ybus;" << std::endl;
+        }
+    if (solver_control.need_reset_solver() ||
+        solver_control.ybus_change_sparsity_pattern() || 
+        solver_control.has_dimension_changed() || 
+        solver_control.need_recompute_ybus()){
+            fillYbus(Ybus, is_ac, id_me_to_solver);
+            std::cout << "fillYbus;" << std::endl;
+        }
+    if (solver_control.need_reset_solver() || 
+        solver_control.has_dimension_changed()) {
+            init_Sbus(Sbus_, id_me_to_solver, id_solver_to_me, slack_bus_id_solver);
+            std::cout << "init_Sbus;" << std::endl;
+        }
+    if (solver_control.need_reset_solver() || 
+        solver_control.has_slack_participate_changed() || 
+        solver_control.has_pv_changed() || 
+        solver_control.has_pq_changed()) {
+            fillpv_pq(id_me_to_solver, id_solver_to_me, slack_bus_id_solver, solver_control);
+            std::cout << "fillpv_pq;" << std::endl;
+        }
 
-    if (solver_control.has_dimension_changed() || solver_control.need_recompute_sbus() && is_ac){
+    if (is_ac && (solver_control.need_reset_solver() || 
+                  solver_control.has_dimension_changed() || 
+                  solver_control.need_recompute_sbus())){
         int nb_bus_total = static_cast<int>(bus_vn_kv_.size());
         total_q_min_per_bus_ = RealVect::Constant(nb_bus_total, 0.);
         total_q_max_per_bus_ = RealVect::Constant(nb_bus_total, 0.);
         total_gen_per_bus_ = Eigen::VectorXi::Constant(nb_bus_total, 0);
         generators_.init_q_vector(nb_bus_total, total_gen_per_bus_, total_q_min_per_bus_, total_q_max_per_bus_);
         dc_lines_.init_q_vector(nb_bus_total, total_gen_per_bus_, total_q_min_per_bus_, total_q_max_per_bus_);
-        fillSbus_me(Sbus_, is_ac, id_me_to_solver);
+        std::cout << "total_gen_per_bus_;" << std::endl;
     }
 
+    if (solver_control.need_reset_solver() || 
+        solver_control.has_slack_participate_changed() || 
+        solver_control.has_pq_changed()) {
+            fillSbus_me(Sbus_, is_ac, id_me_to_solver);
+            std::cout << "fillSbus_me;" << std::endl;
+        }
+
     const int nb_bus_solver = static_cast<int>(id_solver_to_me.size());
     CplxVect V = CplxVect::Constant(nb_bus_solver, init_vm_pu_);
     for(int bus_solver_id = 0; bus_solver_id < nb_bus_solver; ++bus_solver_id){
@@ -442,6 +484,7 @@ CplxVect GridModel::pre_process_solver(const CplxVect & Vinit,
     }
     generators_.set_vm(V, id_me_to_solver);
     dc_lines_.set_vm(V, id_me_to_solver);
+    std::cout << nb_bus_solver << std::endl;
     return V;
 }
 
@@ -701,15 +744,17 @@ CplxVect GridModel::dc_pf(const CplxVect & Vinit,
     CplxVect V = pre_process_solver(Vinit, Ybus_dc_,
                                     id_me_to_dc_solver_, id_dc_solver_to_me_, slack_bus_id_dc_solver_,
                                     is_ac, solver_control_);
-
+    std::cout << "after pre proces\n";
     // start the solver
     if(solver_control_.has_slack_participate_changed() || 
        solver_control_.has_pv_changed() || 
        solver_control_.has_slack_weight_changed()) slack_weights_ = generators_.get_slack_weights(Ybus_dc_.rows(), id_me_to_dc_solver_);
+    std::cout << "slack_weights\n";
     conv = _dc_solver.compute_pf(Ybus_dc_, V, dcSbus_, slack_bus_id_dc_solver_, slack_weights_, bus_pv_, bus_pq_, max_iter, tol);
-
+    std::cout << "after compute_pf\n";
     // store results (fase -> because I am in dc mode)
     process_results(conv, res, Vinit, false, id_me_to_dc_solver_);
+    std::cout << "after compute_pf\n";
     return res;
 }
 

src/GridModel.h

@@ -76,8 +76,11 @@ class GridModel : public DataGeneric
           solver_control_(),
           init_vm_pu_(1.04),
           sn_mva_(1.0){
+            _solver.change_solver(SolverType::SparseLU);
             _dc_solver.change_solver(SolverType::DC);
             _solver.set_gridmodel(this);
+            _dc_solver.set_gridmodel(this);
+            solver_control_.tell_all_changed();
         }
         GridModel(const GridModel & other);
         GridModel copy() const{
@@ -166,7 +169,6 @@ class GridModel : public DataGeneric
                         const RealVect & trafo_x,
                         const CplxVect & trafo_b,
                         const RealVect & trafo_tap_step_pct,
-            //                        const RealVect & trafo_tap_step_degree,
                         const RealVect & trafo_tap_pos,
                         const RealVect & trafo_shift_degree,
                         const std::vector<bool> & trafo_tap_hv,  // is tap on high voltage (true) or low voltate