preparing LODF [skip ci]

CHANGELOG.rst

@@ -26,6 +26,8 @@ Change Log
 - [IMPROVED] remove some compilation warnings for clang
 - [IMPROVED] possibility to specify generator used as slack by its name when initializing
   from `pypowsybl`.
+- [IMPROVED] removing some warnings when grid2op is not installed
+  (it should not raise any warning as lightsim2grid does not require grid2op)
 
 [0.8.1] 2024-03-26
 --------------------

lightsim2grid/tests/test_ptdf.py

@@ -75,7 +75,6 @@ def test_ptdf_from_ls(self):
             assert solve_error <= self.tol, f"error for line {line_id}: {solve_error:.2e}MW"
         # test powerflow are correct
         res_ptdf2 = np.dot(PTDF2, self.dcSbus * self.gridmodel.get_sn_mva())
-        np.where(np.abs(res_ptdf2 - self.res_powerflow) >= 1e3)
         assert np.abs(res_ptdf2 - self.res_powerflow).max() <= self.tol, f"max error for powerflow: {np.abs(res_ptdf2 - self.res_powerflow).max():.2e}MW"
 
 

src/ChooseSolver.h

@@ -276,6 +276,20 @@ class ChooseSolver
             return res;
         }
 
+        RealMat get_lodf(const Eigen::SparseMatrix<cplx_type> & dcYbus,
+                         const IntVect & from_bus,
+                         const IntVect & to_bus){
+                if(_solver_type != SolverType::DC && 
+                   _solver_type != SolverType::KLUDC && 
+                   _solver_type != SolverType::NICSLUDC &&
+                   _solver_type != SolverType::CKTSODC){
+                throw std::runtime_error("ChooseSolver::get_lodf: cannot get ptdf for a solver that is not DC.");
+                }
+            auto p_solver = get_prt_solver("get_lodf", true);
+            const auto & res =  p_solver -> get_lodf(dcYbus, from_bus, to_bus);
+            return res;
+        }
+
         void tell_solver_control(const SolverControl & solver_control){
             auto p_solver = get_prt_solver("tell_solver_control", false);
             p_solver -> tell_solver_control(solver_control);

src/GridModel.cpp

@@ -920,6 +920,17 @@ RealMat GridModel::get_ptdf(){
     return _dc_solver.get_ptdf(Ybus_dc_);
 }
 
+RealMat GridModel::get_lodf(){
+    if(Ybus_dc_.size() == 0){
+        throw std::runtime_error("GridModel::get_lodf: Cannot get the ptdf without having first computed a DC powerflow.");
+    }
+    IntVect from_bus(powerlines_.nb() + trafos_.nb());
+    IntVect to_bus(powerlines_.nb() + trafos_.nb());
+    from_bus << powerlines_.get_bus_from(), trafos_.get_bus_from();
+    to_bus << powerlines_.get_bus_to(), trafos_.get_bus_to();
+    return _dc_solver.get_lodf(Ybus_dc_, from_bus, to_bus);
+}
+
 Eigen::SparseMatrix<real_type> GridModel::get_Bf(){
     if(Ybus_dc_.size() == 0){
         throw std::runtime_error("GridModel::get_Bf: Cannot get the Bf matrix without having first computed a DC powerflow.");

src/GridModel.h

@@ -303,6 +303,8 @@ class GridModel : public GenericContainer
                        real_type tol  // not used for DC
                        );
         RealMat get_ptdf();
+        RealMat get_lodf();
+
         Eigen::SparseMatrix<real_type> get_Bf();
 
         // ac powerflow

src/linear_solvers/CKTSOSolver.cpp

@@ -33,7 +33,7 @@ ErrorType CKTSOLinearSolver::reset(){
     return ErrorType::NoError;
 }
 
-ErrorType CKTSOLinearSolver::initialize(Eigen::SparseMatrix<real_type> & J){
+ErrorType CKTSOLinearSolver::initialize(const Eigen::SparseMatrix<real_type> & J){
     const long long *oparm = oparm_;
     int ret_ = CKTSO_CreateSolver(&solver_, &iparm_, &oparm);
     if (ret_ < 0){   // fail
@@ -89,10 +89,9 @@ ErrorType CKTSOLinearSolver::initialize(Eigen::SparseMatrix<real_type> & J){
     return err;
 }
 
-ErrorType CKTSOLinearSolver::solve(Eigen::SparseMatrix<real_type> & J, RealVect & b, bool doesnt_need_refactor){
+ErrorType CKTSOLinearSolver::solve(const 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;

src/linear_solvers/CKTSOSolver.h

@@ -61,8 +61,8 @@ class CKTSOLinearSolver
 
         // public api
         ErrorType reset();
-        ErrorType initialize(Eigen::SparseMatrix<real_type> & J);
-        ErrorType solve(Eigen::SparseMatrix<real_type> & J, RealVect & b, bool doesnt_need_refactor);
+        ErrorType initialize(const Eigen::SparseMatrix<real_type> & J);
+        ErrorType solve(const 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

@@ -21,34 +21,45 @@ ErrorType KLULinearSolver::reset(){
     return ErrorType::NoError;
 }
 
-ErrorType KLULinearSolver::initialize(Eigen::SparseMatrix<real_type>&  J){
+ErrorType KLULinearSolver::initialize(const Eigen::SparseMatrix<real_type>&  J){
     // default Eigen representation: column major, which is good for klu !
-    // J is const here, even if it's not said in klu_analyze
+    // J is const here, but `klu_analyze` signature expects arrays and not const arrays
+    // so I const_cast
     const auto n = J.cols();
     common_ = klu_common();
     ErrorType res = ErrorType::NoError; 
-    symbolic_ = klu_analyze(n, J.outerIndexPtr(), J.innerIndexPtr(), &common_);
+    symbolic_ = klu_analyze(n,
+                            const_cast<Eigen::SparseMatrix<real_type>::StorageIndex *>(J.outerIndexPtr()),
+                            const_cast<Eigen::SparseMatrix<real_type>::StorageIndex *>(J.innerIndexPtr()),
+                            &common_);
     if(common_.status != KLU_OK){
         res = ErrorType::SolverAnalyze; 
     }else{
-        numeric_ = klu_factor(J.outerIndexPtr(), J.innerIndexPtr(), J.valuePtr(), symbolic_, &common_);
+        numeric_ = klu_factor(const_cast<Eigen::SparseMatrix<real_type>::StorageIndex *>(J.outerIndexPtr()),
+                              const_cast<Eigen::SparseMatrix<real_type>::StorageIndex *>(J.innerIndexPtr()),
+                              const_cast<real_type*>(J.valuePtr()),
+                              symbolic_, &common_);
         if(common_.status != KLU_OK) res = ErrorType::SolverFactor; 
     }
     return res;
 }
 
-ErrorType KLULinearSolver::solve(Eigen::SparseMatrix<real_type>& J, RealVect & b, bool doesnt_need_refactor){
+ErrorType KLULinearSolver::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 klu_solver.analyze() before calling that)
-    // J is const even if it does not compile if said const
+    // J is const here, but `klu_refactor` signature expects arrays and not const arrays
+    // so I const_cast
     int ok;
     ErrorType err = ErrorType::NoError;
     bool stop = false;
     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
-        ok = klu_refactor(J.outerIndexPtr(), J.innerIndexPtr(), J.valuePtr(), symbolic_, numeric_, &common_);
+        ok = klu_refactor(const_cast<Eigen::SparseMatrix<real_type>::StorageIndex *>(J.outerIndexPtr()),
+                          const_cast<Eigen::SparseMatrix<real_type>::StorageIndex *>(J.innerIndexPtr()),
+                          const_cast<real_type*>(J.valuePtr()),
+                          symbolic_, numeric_, &common_);
         if (ok != 1) {
             // std::cout << "\t KLU: refactor error" << std::endl;
             err = ErrorType::SolverReFactor;

src/linear_solvers/KLUSolver.h

@@ -45,8 +45,8 @@ class KLULinearSolver
 
         // public api
         ErrorType reset();
-        ErrorType initialize(Eigen::SparseMatrix<real_type>& J);
-        ErrorType solve(Eigen::SparseMatrix<real_type>& J, RealVect & b, bool doesnt_need_refactor);
+        ErrorType initialize(const Eigen::SparseMatrix<real_type>& J);
+        ErrorType solve(const 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

@@ -35,7 +35,7 @@ ErrorType NICSLULinearSolver::reset(){
     return ErrorType::NoError;
 }
 
-ErrorType NICSLULinearSolver::initialize(Eigen::SparseMatrix<real_type> & J){
+ErrorType NICSLULinearSolver::initialize(const Eigen::SparseMatrix<real_type> & J){
     // default Eigen representation: column major, which is good for klu !
     const auto n = J.cols(); // should be equal to J_.nrows()
     int ret;
@@ -75,10 +75,9 @@ ErrorType NICSLULinearSolver::initialize(Eigen::SparseMatrix<real_type> & J){
     return err;
 }
 
-ErrorType NICSLULinearSolver::solve(Eigen::SparseMatrix<real_type> & J, RealVect & b, bool doesnt_need_refactor){
+ErrorType NICSLULinearSolver::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 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;

src/linear_solvers/NICSLUSolver.h

@@ -53,8 +53,8 @@ class NICSLULinearSolver
 
         // public api
         ErrorType reset();
-        ErrorType initialize(Eigen::SparseMatrix<real_type> & J);
-        ErrorType solve(Eigen::SparseMatrix<real_type> & J, RealVect & b, bool doesnt_need_refactor);
+        ErrorType initialize(const Eigen::SparseMatrix<real_type> & J);
+        ErrorType solve(const 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

@@ -26,7 +26,6 @@ ErrorType SparseLULinearSolver::initialize(const Eigen::SparseMatrix<real_type>
 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(!doesnt_need_refactor){

src/linear_solvers/SparseLUSolver.h

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

src/main.cpp

@@ -895,6 +895,7 @@ PYBIND11_MODULE(lightsim2grid_cpp, m)
         .def("get_solver_control", &GridModel::get_solver_control, "TODO")
         .def("compute_newton", &GridModel::ac_pf, DocGridModel::ac_pf.c_str())
         .def("get_ptdf", &GridModel::get_ptdf, DocGridModel::_internal_do_not_use.c_str()) // TODO PTDF
+        .def("get_lodf", &GridModel::get_lodf, DocGridModel::_internal_do_not_use.c_str()) // TODO PTDF
         .def("get_Bf", &GridModel::get_Bf, DocGridModel::_internal_do_not_use.c_str()) // TODO PTDF
 
          // apply action faster (optimized for grid2op representation)

src/powerflow_algorithm/BaseAlgo.h

@@ -124,7 +124,9 @@ class BaseAlgo : public BaseConstants
         virtual RealMat get_ptdf(const Eigen::SparseMatrix<cplx_type> & dcYbus){
             throw std::runtime_error("Impossible to get the PTDF matrix with this solver type.");
         }
-        virtual Eigen::SparseMatrix<real_type> get_lodf(){  // TODO interface is likely to change
+        virtual RealMat get_lodf(const Eigen::SparseMatrix<cplx_type> & dcYbus,
+                                 const IntVect & from_bus,
+                                 const IntVect & to_bus){  // TODO interface is likely to change
             throw std::runtime_error("Impossible to get the LODF matrix with this solver type.");
         }
         virtual Eigen::SparseMatrix<real_type> get_bsdf(){  // TODO interface is likely to change

src/powerflow_algorithm/BaseDCAlgo.h

@@ -39,9 +39,11 @@ class BaseDCAlgo: public BaseAlgo
                         real_type tol
                         );
 
-        virtual RealMat get_ptdf(const Eigen::SparseMatrix<cplx_type> & dcYbus);  // TODO PTDF
-        virtual Eigen::SparseMatrix<real_type> get_lodf();  // TODO PTDF
-        virtual Eigen::SparseMatrix<real_type> get_bsdf();  // TODO PTDF
+        virtual RealMat get_ptdf(const Eigen::SparseMatrix<cplx_type> & dcYbus);
+        virtual RealMat get_lodf(const Eigen::SparseMatrix<cplx_type> & dcYbus,
+                                 const IntVect & from_bus,
+                                 const IntVect & to_bus);  // TODO LODF
+        virtual Eigen::SparseMatrix<real_type> get_bsdf();  // TODO BSDF
 
     private:
         // no copy allowed

src/powerflow_algorithm/BaseDCAlgo.tpp

@@ -256,9 +256,9 @@ RealMat BaseDCAlgo<LinearSolver>::get_ptdf(const Eigen::SparseMatrix<cplx_type>
 
     // solve iteratively the linear systems (one per powerline)
     PTDF = RealMat::Zero(Bf_T_with_slack.cols(), Bf_T_with_slack.rows());  // rows and cols are "inverted" because the matrix Bf is transposed
-    for (int line_id=0; line_id < nb_pow_tr; ++line_id){
+    for (int row_id=0; row_id < nb_pow_tr; ++row_id){
         // build the rhs vector
-        for (typename Eigen::SparseMatrix<real_type>::InnerIterator it(Bf_T_with_slack, line_id); it; ++it)
+        for (typename Eigen::SparseMatrix<real_type>::InnerIterator it(Bf_T_with_slack, row_id); it; ++it)
         {
             const auto bus_id = it.row();
             if(mat_bus_id_(bus_id) == -1) continue;  // I don't add anything if it's the slack
@@ -269,7 +269,7 @@ RealMat BaseDCAlgo<LinearSolver>::get_ptdf(const Eigen::SparseMatrix<cplx_type>
         _linear_solver.solve(dcYbus_noslack_, rhs, true);  // I don't need to refactorize the matrix (hence the `true`)
 
         // assign results to the PTDF matrix
-        PTDF(line_id, ind_no_slack) = rhs;
+        PTDF(row_id, ind_no_slack) = rhs;
 
         // reset the rhs vector to 0.
         rhs.array() = 0.;
@@ -280,10 +280,12 @@ RealMat BaseDCAlgo<LinearSolver>::get_ptdf(const Eigen::SparseMatrix<cplx_type>
 }
 
 template<class LinearSolver>
-Eigen::SparseMatrix<real_type> BaseDCAlgo<LinearSolver>::get_lodf(){
-    // TODO
-    return dcYbus_noslack_;
-
+RealMat BaseDCAlgo<LinearSolver>::get_lodf(const Eigen::SparseMatrix<cplx_type> & dcYbus,
+                                                                  const IntVect & from_bus,
+                                                                  const IntVect & to_bus){
+    auto PTDF = get_ptdf(dcYbus);  // size n_bus x n_bus
+    RealMat LODF = RealMat::Zero(from_bus.size(), from_bus.rows());  // nb_line, nb_line
+    return LODF;
 }
 
 template<class LinearSolver>