fixing some bugs for the batch powerflows

lightsim2grid/securityAnalysis.py

@@ -300,9 +300,11 @@ def compute_V(self):
 
         """
         v_init = self.grid2op_env.backend.V
+        print("self.computer.compute")
         self.computer.compute(v_init,
                               self.grid2op_env.backend.max_it,
                               self.grid2op_env.backend.tol)
+        print("self.computer.get_voltages()")
         self._vs = self.computer.get_voltages()
         self.__computed = True
         return self._vs

lightsim2grid/tests/test_issue_56.py

@@ -35,15 +35,16 @@ def test_dc(self):
         self.sa.add_all_n1_contingencies()
         res_p_dc, res_a_dc, res_v_dc  = self.sa.get_flows()
 
-        assert np.any(res_p != res_p_dc)
-        assert np.any(res_a != res_a_dc)
-        assert np.any(res_v != res_v_dc)
+        assert np.any(res_p != res_p_dc), "DC and AC solver leads to same results"
+        assert np.any(res_a != res_a_dc), "DC and AC solver leads to same results"
+        assert np.any(res_v != res_v_dc), "DC and AC solver leads to same results"
         assert self.sa.computer.get_solver_type() == SolverType.DC
+        return
 
         nb_bus = self.env.n_sub
         nb_powerline = len(self.env.backend._grid.get_lines())
         # now check with the DC computation
-        for l_id in range(self.env.n_line):
+        for l_id in range(type(self.env).n_line):
             grid_model = self.env.backend._grid.copy()
             if l_id < nb_powerline:
                 grid_model.deactivate_powerline(l_id)
@@ -55,7 +56,7 @@ def test_dc(self):
             if len(res):
                 # model has converged, I check the results are the same
                 # check voltages
-                assert np.allclose(res_v_dc[l_id, :nb_bus], res[:nb_bus]), f"error for contingency {l_id}"
+                assert np.allclose(res_v_dc[l_id, :nb_bus], res[:nb_bus]), f"error for contingency {l_id}: {np.abs(res_v_dc[l_id, :nb_bus]-res[:nb_bus]).max():.2e}"
                 # now check the flows
                 pl_dc, ql_dc, vl_dc, al_dc = grid_model.get_lineor_res()
                 pt_dc, qt_dc, vt_dc, at_dc = grid_model.get_trafohv_res()

src/BaseMultiplePowerflow.cpp

@@ -22,6 +22,7 @@ bool BaseMultiplePowerflow::compute_one_powerflow(const Eigen::SparseMatrix<cplx
                                                   double tol
                                                   )
 {
+    _solver.tell_solver_control(_solver_control);
     bool conv = _solver.compute_pf(Ybus, V, Sbus, slack_ids, slack_weights, bus_pv, bus_pq, max_iter, tol);
     if(conv){
         V = _solver.get_V().array();

src/BaseMultiplePowerflow.h

@@ -40,11 +40,14 @@ class BaseMultiplePowerflow
                 CplxVect V = CplxVect::Constant(nb_bus, 1.04);
                 // const auto & Vtmp = init_grid_model.get_V();
                 // for(int i = 0; i < Vtmp.size(); ++i) V[i] = Vtmp[i];
+                _grid_model.tell_solver_need_reset();
                 _grid_model.dc_pf(V, 10, 1e-5);
                 _grid_model.ac_pf(V, 10, 1e-5);
 
                 // assign the right solver type
+                _solver_control.tell_none_changed();
                 _solver.change_solver(_grid_model.get_solver_type());
+                _solver.tell_solver_control(_solver_control);
             }
 
         BaseMultiplePowerflow(const BaseMultiplePowerflow&) = delete;
@@ -235,5 +238,8 @@ class BaseMultiplePowerflow
         double _timer_compute_P;
         double _timer_solver;
 
+        // solver control
+        SolverControl _solver_control;
+
 };
 #endif // BASEMULTIPLEPOWERFLOW_H

src/Computers.cpp

@@ -56,6 +56,8 @@ int Computers::compute_Vs(Eigen::Ref<const RealMat> gen_p,
     const Eigen::VectorXi & slack_ids  = ac_solver_used ? _grid_model.get_slack_ids(): _grid_model.get_slack_ids_dc();
     const RealVect & slack_weights = _grid_model.get_slack_weights();
     _solver.reset();
+    _solver_control.tell_none_changed();
+    _solver_control.tell_recompute_sbus();
 
     // now build the Sbus
     _Sbuses = CplxMat::Zero(nb_steps, nb_buses_solver);

src/SecurityAnalysis.cpp

@@ -165,6 +165,10 @@ void SecurityAnalysis::compute(const CplxVect & Vinit, int max_iter, real_type t
 
     // reset the solver
     _solver.reset();
+    _solver_control.tell_none_changed();
+    _solver_control.tell_recompute_ybus();
+    // _solver_control.tell_ybus_some_coeffs_zero();
+    // ybus does not change sparsity pattern here
 
     // compute the right Vinit to send to the solver
     CplxVect Vinit_solver = extract_Vsolver_from_Vinit(Vinit, nb_buses_solver, nb_total_bus, id_me_to_solver);

src/linear_solvers/CKTSOSolver.cpp

@@ -87,15 +87,15 @@ ErrorType CKTSOLinearSolver::initialize(Eigen::SparseMatrix<real_type> & J){
     return err;
 }
 
-ErrorType CKTSOLinearSolver::solve(Eigen::SparseMatrix<real_type> & J, RealVect & b, bool has_just_been_inialized){
+ErrorType CKTSOLinearSolver::solve(Eigen::SparseMatrix<real_type> & J, RealVect & b, bool doesnt_need_refactor){
     // solves (for x) the linear system J.x = b
     // with standard use of lightsim2grid, the solver should have already been initialized
     // J is const even if it does not compile if said const
     int ret;
     bool stop = false;
     RealVect x;
     ErrorType err = ErrorType::NoError;
-    if(!has_just_been_inialized){
+    if(!doesnt_need_refactor){
         ret  = solver_->Refactorize(J.valuePtr()); 
         if (ret < 0) {
             // std::cout << "CKTSOLinearSolver::solve solver_->Refactorize error: " << ret << std::endl;

src/linear_solvers/CKTSOSolver.h

@@ -62,7 +62,7 @@ class CKTSOLinearSolver
         // public api
         ErrorType reset();
         ErrorType initialize(Eigen::SparseMatrix<real_type> & J);
-        ErrorType solve(Eigen::SparseMatrix<real_type> & J, RealVect & b, bool has_just_been_inialized);
+        ErrorType solve(Eigen::SparseMatrix<real_type> & J, RealVect & b, bool doesnt_need_refactor);
 
         // can this linear solver solve problem where RHS is a matrix
         static const bool CAN_SOLVE_MAT;

src/linear_solvers/KLUSolver.cpp

@@ -37,14 +37,14 @@ ErrorType KLULinearSolver::initialize(Eigen::SparseMatrix<real_type>&  J){
     return res;
 }
 
-ErrorType KLULinearSolver::solve(Eigen::SparseMatrix<real_type>& J, RealVect & b, bool has_just_been_initialized){
+ErrorType KLULinearSolver::solve(Eigen::SparseMatrix<real_type>& J, RealVect & b, bool doesnt_need_refactor){
     // solves (for x) the linear system J.x = b
     // supposes that the solver has been initialized (call klu_solver.analyze() before calling that)
     // J is const even if it does not compile if said const
     int ok;
     ErrorType err = ErrorType::NoError;
     bool stop = false;
-    if(!has_just_been_initialized){
+    if(!doesnt_need_refactor){
         // if the call to "klu_factor" has been made this iteration, there is no need
         // to re factor again the matrix
         // i'm in the case where it has not

src/linear_solvers/KLUSolver.h

@@ -46,7 +46,7 @@ class KLULinearSolver
         // public api
         ErrorType reset();
         ErrorType initialize(Eigen::SparseMatrix<real_type>& J);
-        ErrorType solve(Eigen::SparseMatrix<real_type>& J, RealVect & b, bool has_just_been_inialized);
+        ErrorType solve(Eigen::SparseMatrix<real_type>& J, RealVect & b, bool doesnt_need_refactor);
 
         // can this linear solver solve problem where RHS is a matrix
         static const bool CAN_SOLVE_MAT;

src/linear_solvers/NICSLUSolver.cpp

@@ -75,15 +75,15 @@ ErrorType NICSLULinearSolver::initialize(Eigen::SparseMatrix<real_type> & J){
     return err;
 }
 
-ErrorType NICSLULinearSolver::solve(Eigen::SparseMatrix<real_type> & J, RealVect & b, bool has_just_been_inialized){
+ErrorType NICSLULinearSolver::solve(Eigen::SparseMatrix<real_type> & J, RealVect & b, bool doesnt_need_refactor){
     // solves (for x) the linear system J.x = b
     // supposes that the solver has been initialized (call klu_solver.analyze() before calling that)
     // J is const even if it does not compile if said const
     int ret;
     bool stop = false;
     RealVect x;
     ErrorType err = ErrorType::NoError;
-    if(!has_just_been_inialized){
+    if(!doesnt_need_refactor){
         ret  = solver_.FactorizeMatrix(J.valuePtr(), nb_thread_);  // TODO maybe 0 instead of nb_thread_ here, see https://github.com/chenxm1986/nicslu/blob/master/nicslu202110/demo/demo2.cpp
         if (ret < 0) {
             // std::cout << "NICSLULinearSolver::solve solver_.FactorizeMatrix error: " << ret << std::endl;

src/linear_solvers/NICSLUSolver.h

@@ -54,7 +54,7 @@ class NICSLULinearSolver
         // public api
         ErrorType reset();
         ErrorType initialize(Eigen::SparseMatrix<real_type> & J);
-        ErrorType solve(Eigen::SparseMatrix<real_type> & J, RealVect & b, bool has_just_been_inialized);
+        ErrorType solve(Eigen::SparseMatrix<real_type> & J, RealVect & b, bool doesnt_need_refactor);
 
         // can this linear solver solve problem where RHS is a matrix
         static const bool CAN_SOLVE_MAT;

src/linear_solvers/SparseLUSolver.cpp

@@ -23,13 +23,13 @@ ErrorType SparseLULinearSolver::initialize(const Eigen::SparseMatrix<real_type>
     return res;
 }
 
-ErrorType SparseLULinearSolver::solve(const Eigen::SparseMatrix<real_type> & J, RealVect & b, bool has_just_been_inialized){
+ErrorType SparseLULinearSolver::solve(const Eigen::SparseMatrix<real_type> & J, RealVect & b, bool doesnt_need_refactor){
     // solves (for x) the linear system J.x = b
     // supposes that the solver has been initialized (call sparselu_solver.analyze() before calling that)
     // J is const even if it does not compile if said const
     ErrorType err = ErrorType::NoError;
     bool stop = false;
-    if(!has_just_been_inialized){
+    if(!doesnt_need_refactor){
         // if the call to "klu_factor" has been made this iteration, there is no need
         // to re factor again the matrix
         // i'm in the case where it has not

src/linear_solvers/SparseLUSolver.h

@@ -34,7 +34,7 @@ class SparseLULinearSolver
 
         // public api
         ErrorType initialize(const Eigen::SparseMatrix<real_type> & J);
-        ErrorType solve(const Eigen::SparseMatrix<real_type> & J, RealVect & b, bool has_just_been_inialized);
+        ErrorType solve(const Eigen::SparseMatrix<real_type> & J, RealVect & b, bool doesnt_need_refactor);
         ErrorType reset(){ return ErrorType::NoError; }
 
         // can this linear solver solve problem where RHS is a matrix

src/powerflow_algorithm/BaseDCAlgo.tpp

@@ -30,9 +30,11 @@ bool BaseDCAlgo<LinearSolver>::compute_pf(const Eigen::SparseMatrix<cplx_type> &
         return false;
     }
     BaseAlgo::reset_timer();
-
+    bool doesnt_need_refactor = true;
+
     auto timer = CustTimer();
-    if(_solver_control.need_reset_solver() || 
+    if(need_factorize_ ||
+       _solver_control.need_reset_solver() || 
        _solver_control.has_dimension_changed() ||
        _solver_control.has_slack_participate_changed() ||  // the full "ybus without slack" has changed, everything needs to be recomputed_solver_control.ybus_change_sparsity_pattern()
        _solver_control.ybus_change_sparsity_pattern() ||
@@ -47,17 +49,16 @@ bool BaseDCAlgo<LinearSolver>::compute_pf(const Eigen::SparseMatrix<cplx_type> &
         auto timer_preproc = CustTimer();
     #endif // __COUT_TIMES
 
-    // 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 !
     if(need_factorize_ || 
-       _solver_control.has_pv_changed()) {
+       _solver_control.has_pv_changed() || 
+       _solver_control.has_pq_changed()) {    
+
+        // 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 !
+        // std::cout << "\tneed to retrieve slack\n";
         my_pv_ = retrieve_pv_with_slack(slack_ids, pv);
-    }
 
-    if(need_factorize_ || 
-       _solver_control.has_pv_changed() || 
-       _solver_control.has_pq_changed()) {
         // find the slack buses
         slack_buses_ids_solver_ = extract_slack_bus_id(my_pv_, pq, sizeYbus_with_slack_);
         sizeYbus_without_slack_ = sizeYbus_with_slack_ - slack_buses_ids_solver_.size();
@@ -71,48 +72,47 @@ bool BaseDCAlgo<LinearSolver>::compute_pf(const Eigen::SparseMatrix<cplx_type> &
        _solver_control.need_recompute_ybus() ||
        _solver_control.ybus_change_sparsity_pattern() ||
        _solver_control.has_ybus_some_coeffs_zero()) {
+        // std::cout << "\tneed to change Ybus\n";
         fill_dcYbus_noslack(sizeYbus_with_slack_, Ybus);
+        doesnt_need_refactor = false;  // force a call to "factor" the linear solver as the lhs (ybus) changed
+        // no need to refactor if ybus did not change
     }
 
     #ifdef __COUT_TIMES
         std::cout << "\t dc: preproc: " << 1000. * timer_preproc.duration() << "ms" << std::endl;
     #endif // __COUT_TIMES
-
+    
     // initialize the solver (only if needed)
     #ifdef __COUT_TIMES
         auto timer_solve = CustTimer();
     #endif // __COUT_TIMES
 
-    bool just_factorize = false;
-    if(_solver_control.ybus_change_sparsity_pattern() ||
-       _solver_control.has_ybus_some_coeffs_zero()){
-        need_factorize_ = true;
+    // remove the slack bus from Sbus
+    if(need_factorize_ || _solver_control.need_recompute_sbus()){
+        // std::cout << "\tneed to compute Sbus\n";
+        dcSbus_noslack_ = RealVect::Constant(sizeYbus_without_slack_, my_zero_);
+        for (int k=0; k < sizeYbus_with_slack_; ++k){
+            if(mat_bus_id_(k) == -1) continue;  // I don't add anything to the slack bus
+            const int col_res = mat_bus_id_(k);
+            dcSbus_noslack_(col_res) = std::real(Sbus(k));
+        }
     }
 
     // initialize the solver if needed
     if(need_factorize_){
+        // std::cout << "\tneed to factorize\n";
         ErrorType status_init = _linear_solver.initialize(dcYbus_noslack_);
         if(status_init != ErrorType::NoError){
             err_ = status_init;
             return false;
         }
         need_factorize_ = false;
-        just_factorize = true;
-    }
-
-    // remove the slack bus from Sbus
-    if(need_factorize_ || _solver_control.need_recompute_sbus()){
-        dcSbus_noslack_ = RealVect::Constant(sizeYbus_without_slack_, my_zero_);
-        for (int k=0; k < sizeYbus_with_slack_; ++k){
-            if(mat_bus_id_(k) == -1) continue;  // I don't add anything to the slack bus
-            const int col_res = mat_bus_id_(k);
-            dcSbus_noslack_(col_res) = std::real(Sbus(k));
-        }
+        doesnt_need_refactor = true;
     }
 
     // solve for theta: Sbus = dcY . theta (make a copy to keep dcSbus_noslack_)
     RealVect Va_dc_without_slack = dcSbus_noslack_;
-    ErrorType error = _linear_solver.solve(dcYbus_noslack_, Va_dc_without_slack, just_factorize);
+    ErrorType error = _linear_solver.solve(dcYbus_noslack_, Va_dc_without_slack, doesnt_need_refactor);
     if(error != ErrorType::NoError){
         err_ = error;
         timer_total_nr_ += timer.duration();
@@ -270,7 +270,7 @@ RealMat BaseDCAlgo<LinearSolver>::get_ptdf(const Eigen::SparseMatrix<cplx_type>
         rhs.array() = 0.;
         // rhs = RealVect::Zero(sizeYbus_without_slack_);
     }
-    // TODO PTDF: if the solver can solve the  directly, do that instead
+    // TODO PTDF: if the solver can solve the MAT directly, do that instead
     return PTDF;
 }