From 91939fdc4637a2849a847d47bdf096b4e153989f Mon Sep 17 00:00:00 2001 From: DONNOT Benjamin Date: Thu, 14 Dec 2023 16:00:08 +0100 Subject: [PATCH] let's see what tests fail now --- CHANGELOG.rst | 2 + benchmarks/benchmark_solvers.py | 2 +- lightsim2grid/tests/test_GridModel.py | 5 +- lightsim2grid/tests/test_SameResPP.py | 4 +- lightsim2grid/tests/test_case118.py | 2 +- lightsim2grid/tests/test_issue_56.py | 2 +- lightsim2grid/tests/test_ptdf.py | 3 +- src/BaseNRSolver.h | 14 +++++- src/BaseNRSolver.tpp | 71 ++++++++++++++++----------- src/BaseNRSolverSingleSlack.h | 3 +- src/BaseNRSolverSingleSlack.tpp | 39 ++++++++------- src/DataGen.cpp | 4 +- src/DataGen.h | 2 + src/DataLine.h | 13 +++-- src/GridModel.cpp | 40 ++++++++++++--- src/GridModel.h | 10 +++- src/Utils.h | 2 +- 17 files changed, 144 insertions(+), 74 deletions(-) diff --git a/CHANGELOG.rst b/CHANGELOG.rst index 56c4a88..c880b94 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -30,6 +30,8 @@ Change Log - [FIXED] a bug in init from pypowsybl when some object were disconnected. It raises an error (because they are not connected to a bus): now this function properly handles these cases. +- [FIXED] a bug leading to not propagate correctly the "compute_results" flag when the + environment was copied (for example) - [ADDED] sets of methods to extract the main component of a grid and perform powerflow only on this one. - [ADDED] possibility to set / retrieve the names of each elements of the grid. diff --git a/benchmarks/benchmark_solvers.py b/benchmarks/benchmark_solvers.py index 7281e06..55dfa24 100644 --- a/benchmarks/benchmark_solvers.py +++ b/benchmarks/benchmark_solvers.py @@ -43,10 +43,10 @@ lightsim2grid.SolverType.SparseLU: "NR (SLU)", lightsim2grid.SolverType.KLU: "NR (KLU)", lightsim2grid.SolverType.NICSLU: "NR (NICSLU *)", + lightsim2grid.SolverType.CKTSO: "NR (CKTSO *)", lightsim2grid.SolverType.SparseLUSingleSlack: "NR single (SLU)", lightsim2grid.SolverType.KLUSingleSlack: "NR single (KLU)", lightsim2grid.SolverType.NICSLUSingleSlack: "NR single (NICSLU *)", - lightsim2grid.SolverType.CKTSO: "NR (CKTSO *)", lightsim2grid.SolverType.CKTSOSingleSlack: "NR single (CKTSO *)", lightsim2grid.SolverType.FDPF_XB_SparseLU: "FDPF XB (SLU)", lightsim2grid.SolverType.FDPF_BX_SparseLU: "FDPF BX (SLU)", diff --git a/lightsim2grid/tests/test_GridModel.py b/lightsim2grid/tests/test_GridModel.py index 0d6329e..7230fc7 100644 --- a/lightsim2grid/tests/test_GridModel.py +++ b/lightsim2grid/tests/test_GridModel.py @@ -103,7 +103,8 @@ def make_v0(self, net): return V0 def run_me_pf(self, V0): - return self.model.compute_newton(V0, self.max_it, self.tol) + res = self.model.ac_pf(V0, self.max_it, self.tol) + return res def run_ref_pf(self, net): with warnings.catch_warnings(): @@ -111,7 +112,7 @@ def run_ref_pf(self, net): pp.runpp(net, init="flat", lightsim2grid=False, - numba=True, + numba=False, distributed_slack=False) def do_i_skip(self, func_name): diff --git a/lightsim2grid/tests/test_SameResPP.py b/lightsim2grid/tests/test_SameResPP.py index 9e7e1f7..13ae6a8 100644 --- a/lightsim2grid/tests/test_SameResPP.py +++ b/lightsim2grid/tests/test_SameResPP.py @@ -131,7 +131,7 @@ def _aux_test(self, pn_net): V = backend._grid.ac_pf(v_tmp, 10, 1e-5) assert V.shape[0], "? lightsim diverge when initialized with pp final voltage ?" - backend._grid.tell_topo_changed() + backend._grid.tell_solver_need_reset() Y_pp = backend.init_pp_backend._grid._ppc["internal"]["Ybus"] Sbus = backend.init_pp_backend._grid._ppc["internal"]["Sbus"] @@ -218,7 +218,7 @@ def _aux_test(self, pn_net): backend._grid.deactivate_result_computation() Vdc = backend._grid.dc_pf(Vinit, max_iter, tol_this) backend._grid.reactivate_result_computation() - backend._grid.tell_topo_changed() + backend._grid.tell_solver_need_reset() Ydc_me = copy.deepcopy(backend._grid.get_dcYbus()) Sdc_me = copy.deepcopy(backend._grid.get_dcSbus()) assert np.max(np.abs(V_init_ref[pp_vect_converter] - Vdc[:nb_sub])) <= 100.*self.tol,\ diff --git a/lightsim2grid/tests/test_case118.py b/lightsim2grid/tests/test_case118.py index c4e4890..95b5d87 100644 --- a/lightsim2grid/tests/test_case118.py +++ b/lightsim2grid/tests/test_case118.py @@ -131,7 +131,7 @@ def test_neurips_track2(self): with warnings.catch_warnings(): warnings.filterwarnings("ignore") ls_grid = init(self.pp_net) - ls_grid.tell_topo_changed() + ls_grid.tell_solver_need_reset() V = np.ones(2 * self.nb_bus_total, dtype=np.complex_) V = ls_grid.ac_pf(V, self.max_it, self.tol) self.check_results(V[:self.nb_bus_total], ls_grid, self.pp_net) diff --git a/lightsim2grid/tests/test_issue_56.py b/lightsim2grid/tests/test_issue_56.py index 6b0a3a6..9aed2a0 100644 --- a/lightsim2grid/tests/test_issue_56.py +++ b/lightsim2grid/tests/test_issue_56.py @@ -49,7 +49,7 @@ def test_dc(self): grid_model.deactivate_powerline(l_id) else: grid_model.deactivate_trafo(l_id - nb_powerline) - grid_model.tell_topo_changed() + grid_model.tell_solver_need_reset() V = 1.0 * self.env.backend.V # np.ones(2 * self.env.n_sub, dtype=np.complex_) res = grid_model.dc_pf(V, 10, 1e-8) if len(res): diff --git a/lightsim2grid/tests/test_ptdf.py b/lightsim2grid/tests/test_ptdf.py index 8bfb6da..4a98e4a 100644 --- a/lightsim2grid/tests/test_ptdf.py +++ b/lightsim2grid/tests/test_ptdf.py @@ -36,7 +36,8 @@ def setUp(self) -> None: if solver_type not in self.gridmodel.available_solvers(): self.skipTest("Solver type not supported on this platform") self.gridmodel.change_solver(solver_type) - self.gridmodel.dc_pf(self.V_init, 1, 1e-8) + V = self.gridmodel.dc_pf(self.V_init, 1, 1e-8) + assert len(V), f"dc pf has diverged with error {self.gridmodel.get_dc_solver().get_error()}" self.dcYbus = 1.0 * self.gridmodel.get_dcYbus() self.dcSbus = 1.0 * self.gridmodel.get_dcSbus().real self.Bbus = 1.0 * self.dcYbus.real diff --git a/src/BaseNRSolver.h b/src/BaseNRSolver.h index 8b7c5fa..5a46123 100644 --- a/src/BaseNRSolver.h +++ b/src/BaseNRSolver.h @@ -136,8 +136,16 @@ class BaseNRSolver : public BaseSolver void fill_value_map(Eigen::Index slack_bus_id, const Eigen::VectorXi & pq, - const Eigen::VectorXi & pvpq); - + const Eigen::VectorXi & pvpq, + bool reset_J); + + void reset_if_needed(){ + if(_solver_control.need_reset_solver() || + _solver_control.has_dimension_changed() || + _solver_control.has_ybus_some_coeffs_zero()){ + reset(); + } + } protected: // used linear solver LinearSolver _linear_solver; @@ -151,6 +159,8 @@ class BaseNRSolver : public BaseSolver // to store the mapping from the element of J_ in dS_dVm_ and dS_dVa_ // it does not own any memory at all ! std::vector value_map_; + // std::vector col_map_; + // std::vector row_map_; // timers double timer_initialize_; diff --git a/src/BaseNRSolver.tpp b/src/BaseNRSolver.tpp index b099558..4dbbfb8 100644 --- a/src/BaseNRSolver.tpp +++ b/src/BaseNRSolver.tpp @@ -46,20 +46,14 @@ bool BaseNRSolver::compute_pf(const Eigen::SparseMatrix exc_ << "V (" << V.size() << ") and Ybus (" << Ybus.rows()<< ", " << Ybus.cols() << ")."; throw std::runtime_error(exc_.str()); } - reset_timer(); - // std::cout << "dist slack" << std::endl; - - if(_solver_control.need_reset_solver() || - _solver_control.has_dimension_changed()){ - reset(); - } - auto timer = CustTimer(); if(!is_linear_solver_valid()) { // err_ = ErrorType::NotInitError; return false; } - + reset_timer(); + reset_if_needed(); err_ = ErrorType::NoError; // reset the error if previous error happened + auto timer = 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 @@ -94,9 +88,11 @@ bool BaseNRSolver::compute_pf(const Eigen::SparseMatrix bool has_just_been_initialized = false; // to avoid a call to klu_refactor follow a call to klu_factor in the same loop // std::cout << "iter " << nr_iter_ << " dx(0): " << -F(0) << " dx(1): " << -F(1) << std::endl; // std::cout << "slack_absorbed " << slack_absorbed << std::endl; - BaseNRSolver::value_map_.clear(); // TODO smarter solver: only needed if ybus has changed - BaseNRSolver::dS_dVm_.resize(0,0); // TODO smarter solver: only needed if ybus has changed - BaseNRSolver::dS_dVa_.resize(0,0); // TODO smarter solver: only needed if ybus has changed + value_map_.clear(); // TODO smarter solver: only needed if ybus has changed + // BaseNRSolver::col_map_.clear(); // TODO smarter solver: only needed if ybus has changed + // BaseNRSolver::row_map_.clear(); // TODO smarter solver: only needed if ybus has changed + dS_dVm_.resize(0,0); // TODO smarter solver: only needed if ybus has changed + dS_dVa_.resize(0,0); // TODO smarter solver: only needed if ybus has changed // BaseNRSolver::dS_dVm_.setZero(); // TODO smarter solver: only needed if ybus has changed // BaseNRSolver::dS_dVa_.setZero(); // TODO smarter solver: only needed if ybus has changed while ((!converged) & (nr_iter_ < max_iter)){ @@ -181,6 +177,7 @@ void BaseNRSolver::reset(){ template void BaseNRSolver::_dSbus_dV(const Eigen::Ref > & Ybus, const Eigen::Ref & V){ + // std::cout << "Ybus.nonZeros(): " << Ybus.nonZeros() << std::endl; auto timer = CustTimer(); const auto size_dS = V.size(); const CplxVect Vnorm = V.array() / V.array().abs(); @@ -198,6 +195,12 @@ void BaseNRSolver::_dSbus_dV(const Eigen::Ref >::InnerIterator it(Ybus, col_id); it; ++it) @@ -298,14 +301,14 @@ void BaseNRSolver::_get_values_J(int & nb_obj_this_col, template void BaseNRSolver::fill_jacobian_matrix(const Eigen::SparseMatrix & Ybus, - const CplxVect & V, - Eigen::Index slack_bus_id, - const RealVect & slack_weights, - const Eigen::VectorXi & pq, - const Eigen::VectorXi & pvpq, - const std::vector & pq_inv, - const std::vector & pvpq_inv - ) + const CplxVect & V, + Eigen::Index slack_bus_id, + const RealVect & slack_weights, + const Eigen::VectorXi & pq, + const Eigen::VectorXi & pvpq, + const std::vector & pq_inv, + const std::vector & pvpq_inv + ) { /** Remember, J has the shape: @@ -346,7 +349,7 @@ void BaseNRSolver::fill_jacobian_matrix(const Eigen::SparseMatrix< #endif // __COUT_TIMES // first time i initialized the matrix, so i need to compute its sparsity pattern fill_jacobian_matrix_unkown_sparsity_pattern(Ybus, V, slack_bus_id, slack_weights, pq, pvpq, pq_inv, pvpq_inv); - fill_value_map(slack_bus_id, pq, pvpq); + fill_value_map(slack_bus_id, pq, pvpq, false); #ifdef __COUT_TIMES std::cout << "\t\t fill_jacobian_matrix_unkown_sparsity_pattern : " << timer2.duration() << std::endl; #endif // __COUT_TIMES @@ -356,7 +359,7 @@ void BaseNRSolver::fill_jacobian_matrix(const Eigen::SparseMatrix< #ifdef __COUT_TIMES auto timer3 = CustTimer(); #endif // - if (BaseNRSolver::value_map_.size() == 0) fill_value_map(slack_bus_id,pq, pvpq); + if (BaseNRSolver::value_map_.size() == 0) fill_value_map(slack_bus_id, pq, pvpq, true); fill_jacobian_matrix_kown_sparsity_pattern(slack_bus_id, pq, pvpq ); @@ -541,11 +544,15 @@ template void BaseNRSolver::fill_value_map( Eigen::Index slack_bus_id, const Eigen::VectorXi & pq, - const Eigen::VectorXi & pvpq + const Eigen::VectorXi & pvpq, + bool reset_J ) { const int n_pvpq = static_cast(pvpq.size()); - value_map_ = std::vector (J_.nonZeros()); + value_map_ = std::vector (); + value_map_.reserve(BaseNRSolver::J_.nonZeros()); + // col_map_ = std::vector (J_.nonZeros()); + // row_map_ = std::vector (J_.nonZeros()); const auto n_row = J_.cols(); unsigned int pos_el = 0; @@ -554,16 +561,18 @@ void BaseNRSolver::fill_value_map( { auto row_id = it.row(); const auto col_id = it.col() - 1; // it's equal to "col_" + if(reset_J) it.valueRef() = 0.; // "forget" previous J value in this setting + if(row_id==0){ // this is the row of the slack bus const Eigen::Index row_id_dS_dVx_r = slack_bus_id; // same for both matrices if(col_id < n_pvpq){ const int col_id_dS_dVa_r = pvpq[col_id]; - value_map_[pos_el] = &dS_dVa_.coeffRef(row_id_dS_dVx_r, col_id_dS_dVa_r); + value_map_.push_back(&dS_dVa_.coeffRef(row_id_dS_dVx_r, col_id_dS_dVa_r)); } else{ const int col_id_dS_dVm_r = pq[col_id - n_pvpq]; - value_map_[pos_el] = &dS_dVm_.coeffRef(row_id_dS_dVx_r, col_id_dS_dVm_r); + value_map_.push_back(&dS_dVm_.coeffRef(row_id_dS_dVx_r, col_id_dS_dVm_r)); } }else{ row_id -= 1; // "do not consider" the row for slack bus (handled above) @@ -573,7 +582,7 @@ void BaseNRSolver::fill_value_map( const int row_id_dS_dVa_r = pvpq[row_id]; const int col_id_dS_dVa_r = pvpq[col_id]; // this_el = dS_dVa_r.coeff(row_id_dS_dVa_r, col_id_dS_dVa_r); - value_map_[pos_el] = &dS_dVa_.coeffRef(row_id_dS_dVa_r, col_id_dS_dVa_r); + value_map_.push_back(&dS_dVa_.coeffRef(row_id_dS_dVa_r, col_id_dS_dVa_r)); // I don't need to perform these checks: if they failed, the element would not be in J_ in the first place // const int is_row_non_null = pq_inv[row_id_dS_dVa_r]; @@ -587,25 +596,27 @@ void BaseNRSolver::fill_value_map( const int row_id_dS_dVa_i = pq[row_id - n_pvpq]; const int col_id_dS_dVa_i = pvpq[col_id]; // this_el = dS_dVa_i.coeff(row_id_dS_dVa_i, col_id_dS_dVa_i); - value_map_[pos_el] = &dS_dVa_.coeffRef(row_id_dS_dVa_i, col_id_dS_dVa_i); + value_map_.push_back(&dS_dVa_.coeffRef(row_id_dS_dVa_i, col_id_dS_dVa_i)); }else if((col_id >= n_pvpq) && (row_id < n_pvpq)){ // this is the J12 part (dS_dVm_r) const int row_id_dS_dVm_r = pvpq[row_id]; const int col_id_dS_dVm_r = pq[col_id - n_pvpq]; // this_el = dS_dVm_r.coeff(row_id_dS_dVm_r, col_id_dS_dVm_r); - value_map_[pos_el] = &dS_dVm_.coeffRef(row_id_dS_dVm_r, col_id_dS_dVm_r); + value_map_.push_back(&dS_dVm_.coeffRef(row_id_dS_dVm_r, col_id_dS_dVm_r)); }else if((col_id >= n_pvpq) && (row_id >= n_pvpq)){ // this is the J22 part (dS_dVm_i) const int row_id_dS_dVm_i = pq[row_id - n_pvpq]; const int col_id_dS_dVm_i = pq[col_id - n_pvpq]; // this_el = dS_dVm_i.coeff(row_id_dS_dVm_i, col_id_dS_dVm_i); - value_map_[pos_el] = &dS_dVm_.coeffRef(row_id_dS_dVm_i, col_id_dS_dVm_i); + value_map_.push_back(&dS_dVm_.coeffRef(row_id_dS_dVm_i, col_id_dS_dVm_i)); } } // go to the next element ++pos_el; } } + dS_dVa_.makeCompressed(); + dS_dVm_.makeCompressed(); } template diff --git a/src/BaseNRSolverSingleSlack.h b/src/BaseNRSolverSingleSlack.h index 8cd8522..09f8042 100644 --- a/src/BaseNRSolverSingleSlack.h +++ b/src/BaseNRSolverSingleSlack.h @@ -58,7 +58,8 @@ class BaseNRSolverSingleSlack : public BaseNRSolver ); void fill_value_map(const Eigen::VectorXi & pq, - const Eigen::VectorXi & pvpq); + const Eigen::VectorXi & pvpq, + bool reset_J); }; diff --git a/src/BaseNRSolverSingleSlack.tpp b/src/BaseNRSolverSingleSlack.tpp index fd9e1ac..7e39ba2 100644 --- a/src/BaseNRSolverSingleSlack.tpp +++ b/src/BaseNRSolverSingleSlack.tpp @@ -40,19 +40,13 @@ bool BaseNRSolverSingleSlack::compute_pf(const Eigen::SparseMatrix exc_ << "V (" << V.size() << ") and Ybus (" << Ybus.rows()<<", "<::reset_timer(); - // std::cout << "singleslack" << std::endl; - - if(BaseNRSolver::_solver_control.need_reset_solver() || - BaseNRSolver::_solver_control.has_dimension_changed()){ - BaseNRSolver::reset(); - } - if(!BaseNRSolver::is_linear_solver_valid()){ return false; } - + BaseNRSolver::reset_timer(); + BaseNRSolver::reset_if_needed(); BaseNRSolver::err_ = ErrorType::NoError; // reset the error if previous error happened + auto timer = CustTimer(); // initialize once and for all the "inverse" of these vectors // Eigen::VectorXi my_pv = BaseNRSolver::retrieve_pv_with_slack(slack_ids, pv); @@ -84,6 +78,7 @@ bool BaseNRSolverSingleSlack::compute_pf(const Eigen::SparseMatrix BaseNRSolver::value_map_.clear(); // TODO smarter solver: only needed if ybus has changed or pq changed or pv changed BaseNRSolver::dS_dVm_.resize(0,0); // TODO smarter solver: only needed if ybus has changed or pq changed or pv changed BaseNRSolver::dS_dVa_.resize(0,0); // TODO smarter solver: only needed if ybus has changed or pq changed or pv changed + // BaseNRSolver::J_.setZero(); // TODO smarter solver: only needed if ybus has changed or pq changed or pv changed or ybus_some_coeffs_zero_ // BaseNRSolver::dS_dVm_.setZero(); // TODO smarter solver: only needed if ybus has changed // BaseNRSolver::dS_dVa_.setZero(); // TODO smarter solver: only needed if ybus has changed while ((!converged) & (BaseNRSolver::nr_iter_ < max_iter)){ @@ -196,7 +191,7 @@ void BaseNRSolverSingleSlack::fill_jacobian_matrix(const Eigen::Sp #endif // __COUT_TIMES // first time i initialized the matrix, so i need to compute its sparsity pattern fill_jacobian_matrix_unkown_sparsity_pattern(Ybus, V, pq, pvpq, pq_inv, pvpq_inv); - fill_value_map(pq, pvpq); + fill_value_map(pq, pvpq, false); // std::cout << "\t\tfill_jacobian_matrix_unkown_sparsity_pattern" << std::endl; #ifdef __COUT_TIMES std::cout << "\t\t fill_jacobian_matrix_unkown_sparsity_pattern : " << timer2.duration() << std::endl; @@ -209,7 +204,7 @@ void BaseNRSolverSingleSlack::fill_jacobian_matrix(const Eigen::Sp #endif // __COUT_TIMES if (BaseNRSolver::value_map_.size() == 0){ // std::cout << "\t\tfill_value_map called" << std::endl; - fill_value_map(pq, pvpq); + fill_value_map(pq, pvpq, true); } fill_jacobian_matrix_kown_sparsity_pattern(pq, pvpq); // std::cout << "\t\tfill_jacobian_matrix_kown_sparsity_pattern" << std::endl; @@ -264,9 +259,9 @@ void BaseNRSolverSingleSlack::fill_jacobian_matrix_unkown_sparsity if(BaseNRSolver::J_.cols() != size_j) { need_insert = true; - BaseNRSolver::J_ = Eigen::SparseMatrix(size_j,size_j); + BaseNRSolver::J_ = Eigen::SparseMatrix(size_j, size_j); // pre allocate a large enough matrix - BaseNRSolver::J_.reserve(2*(BaseNRSolver::dS_dVa_.nonZeros()+BaseNRSolver::dS_dVm_.nonZeros())); + BaseNRSolver::J_.reserve(2*(BaseNRSolver::dS_dVa_.nonZeros() + BaseNRSolver::dS_dVm_.nonZeros())); // from an experiment, outerIndexPtr is initialized, with the number of columns // innerIndexPtr and valuePtr are not. } @@ -359,11 +354,14 @@ it requires that J_ is initialized, in compressed mode. template void BaseNRSolverSingleSlack::fill_value_map( const Eigen::VectorXi & pq, - const Eigen::VectorXi & pvpq + const Eigen::VectorXi & pvpq, + bool reset_J ) { const int n_pvpq = static_cast(pvpq.size()); - BaseNRSolver::value_map_ = std::vector (BaseNRSolver::J_.nonZeros()); + BaseNRSolver::value_map_.clear(); + // std::cout << "BaseNRSolver::J_.nonZeros(): " << BaseNRSolver::J_.nonZeros() << std::endl; + BaseNRSolver::value_map_.reserve(BaseNRSolver::J_.nonZeros()); const int n_col = static_cast(BaseNRSolver::J_.cols()); unsigned int pos_el = 0; @@ -372,13 +370,14 @@ void BaseNRSolverSingleSlack::fill_value_map( { const int row_id = static_cast(it.row()); const int col_id = static_cast(it.col()); // it's equal to "col_" + if(reset_J) it.valueRef() = 0.; // "forget" previous J value in this setting // real_type & this_el = J_x_ptr[pos_el]; if((col_id < n_pvpq) && (row_id < n_pvpq)){ // this is the J11 part (dS_dVa_r) const int row_id_dS_dVa_r = pvpq[row_id]; const int col_id_dS_dVa_r = pvpq[col_id]; // this_el = dS_dVa_r.coeff(row_id_dS_dVa_r, col_id_dS_dVa_r); - BaseNRSolver::value_map_[pos_el] = &BaseNRSolver::dS_dVa_.coeffRef(row_id_dS_dVa_r, col_id_dS_dVa_r); + BaseNRSolver::value_map_.push_back(&BaseNRSolver::dS_dVa_.coeffRef(row_id_dS_dVa_r, col_id_dS_dVa_r)); // I don't need to perform these checks: if they failed, the element would not be in J_ in the first place // const int is_row_non_null = pq_inv[row_id_dS_dVa_r]; @@ -393,25 +392,27 @@ void BaseNRSolverSingleSlack::fill_value_map( const int row_id_dS_dVa_i = pq[row_id - n_pvpq]; const int col_id_dS_dVa_i = pvpq[col_id]; // this_el = dS_dVa_i.coeff(row_id_dS_dVa_i, col_id_dS_dVa_i); - BaseNRSolver::value_map_[pos_el] = &BaseNRSolver::dS_dVa_.coeffRef(row_id_dS_dVa_i, col_id_dS_dVa_i); + BaseNRSolver::value_map_.push_back(&BaseNRSolver::dS_dVa_.coeffRef(row_id_dS_dVa_i, col_id_dS_dVa_i)); }else if((col_id >= n_pvpq) && (row_id < n_pvpq)){ // this is the J12 part (dS_dVm_r) const int row_id_dS_dVm_r = pvpq[row_id]; const int col_id_dS_dVm_r = pq[col_id - n_pvpq]; // this_el = dS_dVm_r.coeff(row_id_dS_dVm_r, col_id_dS_dVm_r); - BaseNRSolver::value_map_[pos_el] = &BaseNRSolver::dS_dVm_.coeffRef(row_id_dS_dVm_r, col_id_dS_dVm_r); + BaseNRSolver::value_map_.push_back(&BaseNRSolver::dS_dVm_.coeffRef(row_id_dS_dVm_r, col_id_dS_dVm_r)); }else if((col_id >= n_pvpq) && (row_id >= n_pvpq)){ // this is the J22 part (dS_dVm_i) const int row_id_dS_dVm_i = pq[row_id - n_pvpq]; const int col_id_dS_dVm_i = pq[col_id - n_pvpq]; // this_el = dS_dVm_i.coeff(row_id_dS_dVm_i, col_id_dS_dVm_i); - BaseNRSolver::value_map_[pos_el] = &BaseNRSolver::dS_dVm_.coeffRef(row_id_dS_dVm_i, col_id_dS_dVm_i); + BaseNRSolver::value_map_.push_back(&BaseNRSolver::dS_dVm_.coeffRef(row_id_dS_dVm_i, col_id_dS_dVm_i)); } // go to the next element ++pos_el; } } + // BaseNRSolver::dS_dVa_.makeCompressed(); + // BaseNRSolver::dS_dVm_.makeCompressed(); } template diff --git a/src/DataGen.cpp b/src/DataGen.cpp index c504fc4..e92152e 100644 --- a/src/DataGen.cpp +++ b/src/DataGen.cpp @@ -333,6 +333,7 @@ void DataGen::set_vm(CplxVect & V, const std::vector & id_grid_to_solver) c Eigen::VectorXi DataGen::get_slack_bus_id() const{ std::vector tmp; + tmp.reserve(gen_slackbus_.size()); Eigen::VectorXi res; const auto nb_gen = nb(); for(int gen_id = 0; gen_id < nb_gen; ++gen_id){ @@ -342,7 +343,8 @@ Eigen::VectorXi DataGen::get_slack_bus_id() const{ if(!is_in_vect(my_bus, tmp)) tmp.push_back(my_bus); } } - res = Eigen::VectorXi::Map(&tmp[0], tmp.size()); // force the copy of the data apparently + if(tmp.empty()) throw std::runtime_error("DataGen::get_slack_bus_id: no generator are tagged slack bus for this grid."); + res = Eigen::VectorXi::Map(tmp.data(), tmp.size()); // force the copy of the data apparently return res; } diff --git a/src/DataGen.h b/src/DataGen.h index c143017..0f531cb 100644 --- a/src/DataGen.h +++ b/src/DataGen.h @@ -234,6 +234,7 @@ class DataGen: public DataGeneric void deactivate(int gen_id, SolverControl & solver_control) { if (status_[gen_id]){ solver_control.tell_recompute_sbus(); + solver_control.tell_pq_changed(); // bus might now be pq if(voltage_regulator_on_[gen_id]) solver_control.tell_v_changed(); if(!turnedoff_gen_pv_) solver_control.tell_pv_changed(); if(gen_slack_weight_[gen_id] != 0. || gen_slackbus_[gen_id]){ @@ -246,6 +247,7 @@ class DataGen: public DataGeneric void reactivate(int gen_id, SolverControl & solver_control) { if(!status_[gen_id]){ solver_control.tell_recompute_sbus(); + solver_control.tell_pq_changed(); // bus might now be pv if(voltage_regulator_on_[gen_id]) solver_control.tell_v_changed(); if(!turnedoff_gen_pv_) solver_control.tell_pv_changed(); if(gen_slack_weight_[gen_id] != 0. || gen_slackbus_[gen_id]){ diff --git a/src/DataLine.h b/src/DataLine.h index 10b4be0..49f07db 100644 --- a/src/DataLine.h +++ b/src/DataLine.h @@ -189,6 +189,7 @@ class DataLine : public DataGeneric virtual void get_graph(std::vector > & res) const; void deactivate(int powerline_id, SolverControl & solver_control) { + // std::cout << "line: deactivate called\n"; if(status_[powerline_id]){ solver_control.tell_recompute_ybus(); // but sparsity pattern do not change here (possibly one more coeff at 0.) @@ -199,12 +200,18 @@ class DataLine : public DataGeneric void reactivate(int powerline_id, SolverControl & solver_control) { if(!status_[powerline_id]){ solver_control.tell_recompute_ybus(); - solver_control.tell_ybus_change_sparsity_pattern(); // this might change + solver_control.tell_ybus_change_sparsity_pattern(); // sparsity pattern might change: a non zero coeff can pop up } _reactivate(powerline_id, status_); } - void change_bus_or(int powerline_id, int new_bus_id, SolverControl & solver_control, int nb_bus) {_change_bus(powerline_id, new_bus_id, bus_or_id_, solver_control, nb_bus);} - void change_bus_ex(int powerline_id, int new_bus_id, SolverControl & solver_control, int nb_bus) {_change_bus(powerline_id, new_bus_id, bus_ex_id_, solver_control, nb_bus);} + void change_bus_or(int powerline_id, int new_bus_id, SolverControl & solver_control, int nb_bus) { + // std::cout << "line: change_bus_or called\n"; + _change_bus(powerline_id, new_bus_id, bus_or_id_, solver_control, nb_bus); + } + void change_bus_ex(int powerline_id, int new_bus_id, SolverControl & solver_control, int nb_bus) { + // std::cout << "line: change_bus_or called\n"; + _change_bus(powerline_id, new_bus_id, bus_ex_id_, solver_control, nb_bus); + } int get_bus_or(int powerline_id) {return _get_bus(powerline_id, status_, bus_or_id_);} int get_bus_ex(int powerline_id) {return _get_bus(powerline_id, status_, bus_ex_id_);} virtual void fillYbus(std::vector > & res, diff --git a/src/GridModel.cpp b/src/GridModel.cpp index 7bdbea3..b717565 100644 --- a/src/GridModel.cpp +++ b/src/GridModel.cpp @@ -17,6 +17,7 @@ GridModel::GridModel(const GridModel & other) init_vm_pu_ = other.init_vm_pu_; sn_mva_ = other.sn_mva_; + compute_results_ = other.compute_results_; // copy the powersystem representation // 1. bus @@ -72,8 +73,11 @@ GridModel::GridModel(const GridModel & other) _solver.change_solver(other._solver.get_type()); _dc_solver.change_solver(other._dc_solver.get_type()); compute_results_ = other.compute_results_; + solver_control_.tell_all_changed(); _dc_solver.set_gridmodel(this); _solver.set_gridmodel(this); + _dc_solver.tell_solver_control(solver_control_); + _solver.tell_solver_control(solver_control_); } //pickle @@ -264,8 +268,13 @@ void GridModel::reset(bool reset_solver, bool reset_ac, bool reset_dc) dcSbus_ = CplxVect(); bus_pv_ = Eigen::VectorXi(); bus_pq_ = Eigen::VectorXi(); - slack_weights_ = RealVect(); solver_control_.tell_all_changed(); + + slack_bus_id_ac_me_ = Eigen::VectorXi(); // slack bus id, gridmodel number + slack_bus_id_ac_solver_ = Eigen::VectorXi(); // slack bus id, solver number + slack_bus_id_dc_me_ = Eigen::VectorXi(); + slack_bus_id_dc_solver_ = Eigen::VectorXi(); + slack_weights_ = RealVect(); // reset the solvers if (reset_solver){ @@ -274,6 +283,7 @@ void GridModel::reset(bool reset_solver, bool reset_ac, bool reset_dc) _solver.set_gridmodel(this); _dc_solver.set_gridmodel(this); } + } CplxVect GridModel::ac_pf(const CplxVect & Vinit, @@ -301,6 +311,7 @@ CplxVect GridModel::ac_pf(const CplxVect & Vinit, Ybus_ac_, id_me_to_ac_solver_, id_ac_solver_to_me_, + slack_bus_id_ac_me_, slack_bus_id_ac_solver_, is_ac, solver_control_); @@ -397,6 +408,7 @@ CplxVect GridModel::check_solution(const CplxVect & V_proposed, bool check_q_lim Ybus_ac_, id_me_to_ac_solver_, id_ac_solver_to_me_, + slack_bus_id_ac_me_, slack_bus_id_ac_solver_, is_ac, reset_solver); @@ -429,6 +441,7 @@ CplxVect GridModel::pre_process_solver(const CplxVect & Vinit, Eigen::SparseMatrix & Ybus, std::vector & id_me_to_solver, std::vector & id_solver_to_me, + Eigen::VectorXi & slack_bus_id_me, Eigen::VectorXi & slack_bus_id_solver, bool is_ac, const SolverControl & solver_control) @@ -451,7 +464,8 @@ CplxVect GridModel::pre_process_solver(const CplxVect & Vinit, if (solver_control.need_reset_solver() || solver_control.has_dimension_changed() || solver_control.has_slack_participate_changed()){ - slack_bus_id_solver = generators_.get_slack_bus_id(); + // std::cout << "slack_bus_id_solver\n"; + slack_bus_id_me = generators_.get_slack_bus_id(); // this is the slack bus ids with the gridmodel ordering, not the solver ordering. // conversion to solver ordering is done in init_slack_bus @@ -462,6 +476,7 @@ CplxVect GridModel::pre_process_solver(const CplxVect & Vinit, if (solver_control.need_reset_solver() || solver_control.ybus_change_sparsity_pattern() || solver_control.has_dimension_changed()){ + // std::cout << "init_Ybus;" << std::endl; init_Ybus(Ybus, id_me_to_solver, id_solver_to_me); // std::cout << "init_Ybus;" << std::endl; } @@ -469,12 +484,14 @@ CplxVect GridModel::pre_process_solver(const CplxVect & Vinit, solver_control.ybus_change_sparsity_pattern() || solver_control.has_dimension_changed() || solver_control.need_recompute_ybus()){ + // std::cout << "fillYbus;" << std::endl; 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 + // std::cout << "init_Sbus;" << std::endl; Sbus = CplxVect::Constant(id_solver_to_me.size(), 0.); // std::cout << "init_Sbus;" << std::endl; } @@ -483,7 +500,9 @@ CplxVect GridModel::pre_process_solver(const CplxVect & Vinit, solver_control.has_slack_participate_changed() || solver_control.has_pv_changed() || solver_control.has_pq_changed()) { - init_slack_bus(Sbus, id_me_to_solver, id_solver_to_me, slack_bus_id_solver); + // std::cout << "init_slack_bus;" << std::endl; + init_slack_bus(Sbus, id_me_to_solver, id_solver_to_me, slack_bus_id_me, slack_bus_id_solver); + // std::cout << "fillpv_pq;" << std::endl; fillpv_pq(id_me_to_solver, id_solver_to_me, slack_bus_id_solver, solver_control); // std::cout << "fillpv_pq;" << std::endl; } @@ -491,6 +510,7 @@ CplxVect GridModel::pre_process_solver(const CplxVect & Vinit, if (is_ac && (solver_control.need_reset_solver() || solver_control.has_dimension_changed() || solver_control.need_recompute_sbus())){ + // std::cout << "total_gen_per_bus_;" << std::endl; int nb_bus_total = static_cast(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.); @@ -504,6 +524,7 @@ CplxVect GridModel::pre_process_solver(const CplxVect & Vinit, solver_control.has_dimension_changed() || solver_control.has_slack_participate_changed() || solver_control.has_pq_changed()) { + // std::cout << "fillSbus_me;" << std::endl; fillSbus_me(Sbus, is_ac, id_me_to_solver); // std::cout << "fillSbus_me;" << std::endl; } @@ -566,6 +587,7 @@ void GridModel::process_results(bool conv, static_cast(Vinit.size())); } else { //powerflow diverge + // std::cout << "powerflow diverge" << std::endl; reset_results(); // TODO solver control ??? something to do here ? } @@ -598,11 +620,12 @@ void GridModel::init_Ybus(Eigen::SparseMatrix & Ybus, void GridModel::init_slack_bus(const CplxVect & Sbus, const std::vector& id_me_to_solver, const std::vector& id_solver_to_me, + const Eigen::VectorXi & slack_bus_id_me, Eigen::VectorXi & slack_bus_id_solver){ const int nb_bus = static_cast(id_solver_to_me.size()); // slack_bus_id_solver = Eigen::VectorXi::Zero(slack_bus_id_.size()); - // slack_bus_id_solver = Eigen::VectorXi::Constant(slack_bus_id_solver.size(), _deactivated_bus_id); + slack_bus_id_solver = Eigen::VectorXi::Constant(slack_bus_id_me.size(), _deactivated_bus_id); size_t i = 0; // std::cout << "slack_bus_id_solver 2: "; @@ -616,7 +639,7 @@ void GridModel::init_slack_bus(const CplxVect & Sbus, // std::cout << std::endl; // for(auto el: slack_bus_id_) { - for(auto el: slack_bus_id_solver) { + for(auto el: slack_bus_id_me) { auto tmp = id_me_to_solver[el]; if(tmp == _deactivated_bus_id){ std::ostringstream exc_; @@ -624,7 +647,7 @@ void GridModel::init_slack_bus(const CplxVect & Sbus, exc_ << " You can check element "; exc_ << el; exc_ << ": ["; - for(auto el2 : slack_bus_id_solver) exc_ << el2 << ", "; + for(auto el2 : slack_bus_id_me) exc_ << el2 << ", "; exc_ << "]."; throw std::out_of_range(exc_.str()); } @@ -772,6 +795,7 @@ void GridModel::compute_results(bool ac){ } void GridModel::reset_results(){ + // std::cout << "reset_results\n"; powerlines_.reset_results(); // TODO have a function to dispatch that to all type of elements shunts_.reset_results(); trafos_.reset_results(); @@ -812,8 +836,10 @@ CplxVect GridModel::dc_pf(const CplxVect & Vinit, Ybus_dc_, id_me_to_dc_solver_, id_dc_solver_to_me_, + slack_bus_id_dc_me_, slack_bus_id_dc_solver_, - is_ac, solver_control_); + is_ac, + solver_control_); // std::cout << "after pre proces\n"; // start the solver if(solver_control_.need_reset_solver() || diff --git a/src/GridModel.h b/src/GridModel.h index 2bf85b7..19717ba 100644 --- a/src/GridModel.h +++ b/src/GridModel.h @@ -74,6 +74,7 @@ class GridModel : public DataGeneric GridModel(): solver_control_(), + compute_results_(true), init_vm_pu_(1.04), sn_mva_(1.0){ _solver.change_solver(SolverType::SparseLU); @@ -633,6 +634,7 @@ class GridModel : public DataGeneric Eigen::SparseMatrix & Ybus, std::vector & id_me_to_solver, std::vector & id_solver_to_me, + Eigen::VectorXi & slack_bus_id_me, Eigen::VectorXi & slack_bus_id_solver, bool is_ac, const SolverControl & solver_control); @@ -647,7 +649,9 @@ class GridModel : public DataGeneric void init_slack_bus(const CplxVect & Sbus, const std::vector & id_me_to_solver, const std::vector& id_solver_to_me, - Eigen::VectorXi & slack_bus_id_solver); + const Eigen::VectorXi & slack_bus_id_me, + Eigen::VectorXi & slack_bus_id_solver + ); void fillYbus(Eigen::SparseMatrix & res, bool ac, const std::vector& id_me_to_solver); void fillSbus_me(CplxVect & res, bool ac, const std::vector& id_me_to_solver); void fillpv_pq(const std::vector& id_me_to_solver, @@ -802,7 +806,9 @@ class GridModel : public DataGeneric // 8. slack bus // std::vector slack_bus_id_; - Eigen::VectorXi slack_bus_id_ac_solver_; + Eigen::VectorXi slack_bus_id_ac_me_; // slack bus id, gridmodel number + Eigen::VectorXi slack_bus_id_ac_solver_; // slack bus id, solver number + Eigen::VectorXi slack_bus_id_dc_me_; Eigen::VectorXi slack_bus_id_dc_solver_; RealVect slack_weights_; diff --git a/src/Utils.h b/src/Utils.h index 406c125..eeaec05 100644 --- a/src/Utils.h +++ b/src/Utils.h @@ -96,7 +96,7 @@ class SolverControl void tell_none_changed(){ change_dimension_ = false; pv_changed_ = false; - pq_changed_ = true; + pq_changed_ = false; slack_participate_changed_ = false; need_reset_solver_ = false; need_recompute_sbus_ = false;