diff --git a/lightsim2grid/lightSimBackend.py b/lightsim2grid/lightSimBackend.py index a8c66c5..81732db 100644 --- a/lightsim2grid/lightSimBackend.py +++ b/lightsim2grid/lightSimBackend.py @@ -894,9 +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") + # print(f"{self.V[:14] = }") Vdc = self._grid.dc_pf(copy.deepcopy(self.V), self.max_it, self.tol) - print("after dc pf") + # print(f"{Vdc[:14] = }") 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()}") @@ -905,7 +905,6 @@ 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: diff --git a/setup.py b/setup.py index 99193f8..9a0f4f0 100644 --- a/setup.py +++ b/setup.py @@ -157,7 +157,8 @@ "src/BaseMultiplePowerflow.cpp", "src/Computers.cpp", "src/SecurityAnalysis.cpp", - "src/Solvers.cpp"] + "src/Solvers.cpp", + "src/Utils.cpp"] if KLU_SOLVER_AVAILABLE: src_files.append("src/KLUSolver.cpp") diff --git a/src/BaseNRSolver.tpp b/src/BaseNRSolver.tpp index 9c25925..b099558 100644 --- a/src/BaseNRSolver.tpp +++ b/src/BaseNRSolver.tpp @@ -13,15 +13,15 @@ template bool BaseNRSolver::compute_pf(const Eigen::SparseMatrix & 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 @@ -36,19 +36,28 @@ bool BaseNRSolver::compute_pf(const Eigen::SparseMatrix // TODO DEBUG MODE std::ostringstream exc_; exc_ << "BaseNRSolver::compute_pf: Size of the Sbus should be the same as the size of Ybus. Currently: "; - exc_ << "Sbus (" << Sbus.size() << ") and Ybus (" << Ybus.rows() << ", " << Ybus.rows() << ")."; + exc_ << "Sbus (" << Sbus.size() << ") and Ybus (" << Ybus.rows() << ", " << Ybus.cols() << ")."; throw std::runtime_error(exc_.str()); } if(V.size() != Ybus.rows() || V.size() != Ybus.cols() ){ // TODO DEBUG MODE std::ostringstream exc_; exc_ << "BaseNRSolver::compute_pf: Size of V (init voltages) should be the same as the size of Ybus. Currently: "; - exc_ << "V (" << V.size() << ") and Ybus (" << Ybus.rows()<< ", " << Ybus.rows() << ")."; + 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()) return false; + if(!is_linear_solver_valid()) { + // err_ = ErrorType::NotInitError; + return false; + } err_ = ErrorType::NoError; // reset the error if previous error happened @@ -85,6 +94,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 + // 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)){ nr_iter_++; fill_jacobian_matrix(Ybus, V_, slack_bus_id, slack_weights, pq, pvpq, pq_inv, pvpq_inv); @@ -146,6 +160,7 @@ bool BaseNRSolver::compute_pf(const Eigen::SparseMatrix << "\n\t timer_total_nr_: " << timer_total_nr_ << "\n\n"; #endif // __COUT_TIMES + _solver_control.tell_none_changed(); return res; } @@ -165,7 +180,7 @@ void BaseNRSolver::reset(){ template void BaseNRSolver::_dSbus_dV(const Eigen::Ref > & Ybus, - const Eigen::Ref & V){ + const Eigen::Ref & V){ auto timer = CustTimer(); const auto size_dS = V.size(); const CplxVect Vnorm = V.array() / V.array().abs(); @@ -331,6 +346,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); #ifdef __COUT_TIMES std::cout << "\t\t fill_jacobian_matrix_unkown_sparsity_pattern : " << timer2.duration() << std::endl; #endif // __COUT_TIMES @@ -339,7 +355,8 @@ void BaseNRSolver::fill_jacobian_matrix(const Eigen::SparseMatrix< // properly and faster (approx 3 times faster than the previous one) #ifdef __COUT_TIMES auto timer3 = CustTimer(); - #endif // __COUT_TIMES + #endif // + if (BaseNRSolver::value_map_.size() == 0) fill_value_map(slack_bus_id,pq, pvpq); fill_jacobian_matrix_kown_sparsity_pattern(slack_bus_id, pq, pvpq ); @@ -404,7 +421,7 @@ void BaseNRSolver::fill_jacobian_matrix_unkown_sparsity_pattern( // optim : if the matrix was already computed, i don't initialize it, i instead reuse as much as i can // i can do that because the matrix will ALWAYS have the same non zero coefficients. // in this if, i allocate it in a "large enough" place to avoid copy when first filling it - if(J_.cols() != size_j) J_ = Eigen::SparseMatrix(size_j,size_j); + if(J_.cols() != size_j) J_ = Eigen::SparseMatrix(size_j, size_j); std::vector > coeffs; // HERE FOR PERF OPTIM (3) coeffs.reserve(2*(dS_dVa_.nonZeros()+dS_dVm_.nonZeros()) + slack_weights.size()); // HERE FOR PERF OPTIM (3) @@ -512,7 +529,6 @@ void BaseNRSolver::fill_jacobian_matrix_unkown_sparsity_pattern( J_.setFromTriplets(coeffs.begin(), coeffs.end()); // HERE FOR PERF OPTIM (3) // std::cout << "end fill jacobian unknown " << std::endl; J_.makeCompressed(); - fill_value_map(slack_bus_id, pq, pvpq); // std::cout << "end fill_value_map" << std::endl; } diff --git a/src/BaseNRSolverSingleSlack.tpp b/src/BaseNRSolverSingleSlack.tpp index df0a6e8..fd9e1ac 100644 --- a/src/BaseNRSolverSingleSlack.tpp +++ b/src/BaseNRSolverSingleSlack.tpp @@ -11,15 +11,15 @@ template bool BaseNRSolverSingleSlack::compute_pf(const Eigen::SparseMatrix & Ybus, - CplxVect & V, - const CplxVect & Sbus, - const Eigen::VectorXi & slack_ids, - const RealVect & slack_weights, // unused here - 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, // unused here + 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 @@ -31,18 +31,26 @@ bool BaseNRSolverSingleSlack::compute_pf(const Eigen::SparseMatrix if(Sbus.size() != Ybus.rows() || Sbus.size() != Ybus.cols() ){ std::ostringstream exc_; exc_ << "BaseNRSolverSingleSlack::compute_pf: Size of the Sbus should be the same as the size of Ybus. Currently: "; - exc_ << "Sbus (" << Sbus.size() << ") and Ybus (" << Ybus.rows() << ", " << Ybus.rows() << ")."; + exc_ << "Sbus (" << Sbus.size() << ") and Ybus (" << Ybus.rows() << ", " << Ybus.cols() << ")."; throw std::runtime_error(exc_.str()); } if(V.size() != Ybus.rows() || V.size() != Ybus.cols() ){ std::ostringstream exc_; exc_ << "BaseNRSolverSingleSlack::compute_pf: Size of V (init voltages) should be the same as the size of Ybus. Currently: "; - exc_ << "V (" << V.size() << ") and Ybus (" << Ybus.rows()<<", "<::reset_timer(); - - if(!BaseNRSolver::is_linear_solver_valid()) return false; + // 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::err_ = ErrorType::NoError; // reset the error if previous error happened auto timer = CustTimer(); @@ -73,14 +81,20 @@ bool BaseNRSolverSingleSlack::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 const cplx_type m_i = BaseNRSolver::my_i; // otherwise it does not compile - + 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::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)){ BaseNRSolver::nr_iter_++; + // std::cout << "\tnr_iter_ " << BaseNRSolver::nr_iter_ << std::endl; fill_jacobian_matrix(Ybus, BaseNRSolver::V_, pq, pvpq, pq_inv, pvpq_inv); if(BaseNRSolver::need_factorize_){ BaseNRSolver::initialize(); if(BaseNRSolver::err_ != ErrorType::NoError){ // I got an error during the initialization of the linear system, i need to stop here + // std::cout << BaseNRSolver::err_ << std::endl; res = false; break; } @@ -95,6 +109,7 @@ bool BaseNRSolverSingleSlack::compute_pf(const Eigen::SparseMatrix has_just_been_initialized = false; if(BaseNRSolver::err_ != ErrorType::NoError){ // I got an error during the solving of the linear system, i need to stop here + // std::cout << BaseNRSolver::err_ << std::endl; res = false; break; } @@ -117,7 +132,11 @@ bool BaseNRSolverSingleSlack::compute_pf(const Eigen::SparseMatrix F = BaseNRSolver::_evaluate_Fx(Ybus, BaseNRSolver::V_, Sbus, my_pv, pq); bool tmp = F.allFinite(); - if(!tmp) break; // divergence due to Nans + if(!tmp){ + BaseNRSolver::err_ = ErrorType::InifiniteValue; + // std::cout << BaseNRSolver::err_ << std::endl; + break; // divergence due to Nans + } converged = BaseNRSolver::_check_for_convergence(F, tol); } if(!converged){ @@ -135,17 +154,18 @@ bool BaseNRSolverSingleSlack::compute_pf(const Eigen::SparseMatrix << "\n\t timer_total_nr_: " << BaseNRSolver::timer_total_nr_ << "\n\n"; #endif // __COUT_TIMES + BaseNRSolver::_solver_control.tell_none_changed(); return res; } template void BaseNRSolverSingleSlack::fill_jacobian_matrix(const Eigen::SparseMatrix & Ybus, - const CplxVect & V, - const Eigen::VectorXi & pq, - const Eigen::VectorXi & pvpq, - const std::vector & pq_inv, - const std::vector & pvpq_inv - ) + const CplxVect & V, + const Eigen::VectorXi & pq, + const Eigen::VectorXi & pvpq, + const std::vector & pq_inv, + const std::vector & pvpq_inv + ) { /** J has the shape @@ -165,7 +185,6 @@ void BaseNRSolverSingleSlack::fill_jacobian_matrix(const Eigen::Sp const int n_pvpq = static_cast(pvpq.size()); const int n_pq = static_cast(pq.size()); const int size_j = n_pvpq + n_pq; - // TODO to gain a bit more time below, try to compute directly, in _dSbus_dV(Ybus, V); // TODO the `dS_dVa_[pvpq, pvpq]` // TODO so that it's easier to retrieve in the next few lines ! @@ -177,6 +196,8 @@ 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); + // 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; #endif // __COUT_TIMES @@ -186,7 +207,12 @@ void BaseNRSolverSingleSlack::fill_jacobian_matrix(const Eigen::Sp #ifdef __COUT_TIMES auto timer3 = CustTimer(); #endif // __COUT_TIMES + if (BaseNRSolver::value_map_.size() == 0){ + // std::cout << "\t\tfill_value_map called" << std::endl; + fill_value_map(pq, pvpq); + } fill_jacobian_matrix_kown_sparsity_pattern(pq, pvpq); + // std::cout << "\t\tfill_jacobian_matrix_kown_sparsity_pattern" << std::endl; #ifdef __COUT_TIMES std::cout << "\t\t fill_jacobian_matrix_kown_sparsity_pattern : " << timer3.duration() << std::endl; #endif // __COUT_TIMES @@ -323,7 +349,6 @@ void BaseNRSolverSingleSlack::fill_jacobian_matrix_unkown_sparsity } // J_.setFromTriplets(coeffs.begin(), coeffs.end()); // HERE FOR PERF OPTIM (3) BaseNRSolver::J_.makeCompressed(); - fill_value_map(pq, pvpq); } /** @@ -340,9 +365,9 @@ void BaseNRSolverSingleSlack::fill_value_map( const int n_pvpq = static_cast(pvpq.size()); BaseNRSolver::value_map_ = std::vector (BaseNRSolver::J_.nonZeros()); - const int n_row = static_cast(BaseNRSolver::J_.cols()); + const int n_col = static_cast(BaseNRSolver::J_.cols()); unsigned int pos_el = 0; - for (int col_=0; col_ < n_row; ++col_){ + for (int col_=0; col_ < n_col; ++col_){ for (Eigen::SparseMatrix::InnerIterator it(BaseNRSolver::J_, col_); it; ++it) { const int row_id = static_cast(it.row()); diff --git a/src/BaseSolver.cpp b/src/BaseSolver.cpp index 61b9d49..624b523 100644 --- a/src/BaseSolver.cpp +++ b/src/BaseSolver.cpp @@ -122,6 +122,7 @@ bool BaseSolver::_check_for_convergence(const RealVect & F, { auto timer = CustTimer(); const auto norm_inf = F.lpNorm(); + // std::cout << "\t\tnorm_inf: " << norm_inf << std::endl; bool res = norm_inf < tol; timer_check_ += timer.duration(); return res; diff --git a/src/DCSolver.tpp b/src/DCSolver.tpp index fee67b2..87c2fb3 100644 --- a/src/DCSolver.tpp +++ b/src/DCSolver.tpp @@ -26,6 +26,15 @@ bool BaseDCSolver::compute_pf(const Eigen::SparseMatrix // V is used the following way: at pq buses it's completely ignored. For pv bus only the magnitude is used, // and for the slack bus both the magnitude and the angle are used. + if(!is_linear_solver_valid()) { + // err_ = ErrorType::NotInitError; + return false; + } + if(_solver_control.need_reset_solver() || + _solver_control.has_dimension_changed()){ + reset(); + } + auto timer = CustTimer(); BaseSolver::reset_timer(); sizeYbus_with_slack_ = static_cast(Ybus.rows()); @@ -135,7 +144,7 @@ bool BaseDCSolver::compute_pf(const Eigen::SparseMatrix #ifdef __COUT_TIMES std::cout << "\t dc postproc: " << 1000. * timer_postproc.duration() << "ms" << std::endl; #endif // __COUT_TIMES - + _solver_control.tell_none_changed(); timer_total_nr_ += timer.duration(); return true; } diff --git a/src/DataDCLine.h b/src/DataDCLine.h index 8f2576a..b9363d2 100644 --- a/src/DataDCLine.h +++ b/src/DataDCLine.h @@ -224,7 +224,7 @@ class DataDCLine : public DataGeneric virtual void fillpv(std::vector& bus_pv, std::vector & has_bus_been_added, - Eigen::VectorXi & slack_bus_id_solver, + const Eigen::VectorXi & slack_bus_id_solver, const std::vector & id_grid_to_solver) const { from_gen_.fillpv(bus_pv, has_bus_been_added, slack_bus_id_solver, id_grid_to_solver); to_gen_.fillpv(bus_pv, has_bus_been_added, slack_bus_id_solver, id_grid_to_solver); diff --git a/src/DataGen.cpp b/src/DataGen.cpp index bc60fd2..c504fc4 100644 --- a/src/DataGen.cpp +++ b/src/DataGen.cpp @@ -168,7 +168,7 @@ void DataGen::fillSbus(CplxVect & Sbus, const std::vector & id_grid_to_solv void DataGen::fillpv(std::vector & bus_pv, std::vector & has_bus_been_added, - Eigen::VectorXi & slack_bus_id_solver, + const Eigen::VectorXi & slack_bus_id_solver, const std::vector & id_grid_to_solver) const { const int nb_gen = nb(); @@ -216,7 +216,7 @@ void DataGen::reset_results(){ res_q_ = RealVect(); // in MVar res_v_ = RealVect(); // in kV res_theta_ = RealVect(); // in deg - bus_slack_weight_ = RealVect(); + // bus_slack_weight_ = RealVect(); } void DataGen::get_vm_for_dc(RealVect & Vm){ @@ -331,16 +331,18 @@ void DataGen::set_vm(CplxVect & V, const std::vector & id_grid_to_solver) c } } -std::vector DataGen::get_slack_bus_id() const{ - std::vector res; +Eigen::VectorXi DataGen::get_slack_bus_id() const{ + std::vector tmp; + Eigen::VectorXi res; const auto nb_gen = nb(); for(int gen_id = 0; gen_id < nb_gen; ++gen_id){ if(gen_slackbus_[gen_id]){ const auto my_bus = bus_id_(gen_id); // do not add twice the same "slack bus" - if(!is_in_vect(my_bus, res)) res.push_back(my_bus); + 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 return res; } @@ -349,7 +351,7 @@ void DataGen::set_p_slack(const RealVect& node_mismatch, { if(bus_slack_weight_.size() == 0){ // TODO DEBUG MODE: perform this check only in debug mode - throw std::runtime_error("Impossible to set the active value of generators for the slack bus"); + throw std::runtime_error("DataGen::set_p_slack: Impossible to set the active value of generators for the slack bus: no known slack (you should haved called DataGen::get_slack_weights first)"); } const auto nb_gen = nb(); for(int gen_id = 0; gen_id < nb_gen; ++gen_id){ diff --git a/src/DataGen.h b/src/DataGen.h index e2fb9cb..c143017 100644 --- a/src/DataGen.h +++ b/src/DataGen.h @@ -221,7 +221,7 @@ class DataGen: public DataGeneric **/ RealVect get_slack_weights(Eigen::Index nb_bus_solver, const std::vector & id_grid_to_solver); - std::vector get_slack_bus_id() const; + Eigen::VectorXi get_slack_bus_id() const; void set_p_slack(const RealVect& node_mismatch, const std::vector & id_grid_to_solver); // modification @@ -236,7 +236,7 @@ class DataGen: public DataGeneric solver_control.tell_recompute_sbus(); 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]){ + if(gen_slack_weight_[gen_id] != 0. || gen_slackbus_[gen_id]){ solver_control.tell_slack_participate_changed(); solver_control.tell_slack_weight_changed(); } @@ -248,14 +248,18 @@ class DataGen: public DataGeneric solver_control.tell_recompute_sbus(); 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]){ + if(gen_slack_weight_[gen_id] != 0. || gen_slackbus_[gen_id]){ solver_control.tell_slack_participate_changed(); solver_control.tell_slack_weight_changed(); } } _reactivate(gen_id, status_); } - void change_bus(int gen_id, int new_bus_id, SolverControl & solver_control, int nb_bus) {_change_bus(gen_id, new_bus_id, bus_id_, solver_control, nb_bus);} + void change_bus(int gen_id, int new_bus_id, SolverControl & solver_control, int nb_bus) { + if (new_bus_id != bus_id_[gen_id]){ + if (gen_slack_weight_[gen_id] != 0. || gen_slackbus_[gen_id]) solver_control.has_slack_participate_changed(); + } + _change_bus(gen_id, new_bus_id, bus_id_, solver_control, nb_bus);} int get_bus(int gen_id) {return _get_bus(gen_id, status_, bus_id_);} virtual void reconnect_connected_buses(std::vector & bus_status) const; virtual void disconnect_if_not_in_main_component(std::vector & busbar_in_main_component); @@ -269,7 +273,7 @@ class DataGen: public DataGeneric virtual void fillSbus(CplxVect & Sbus, const std::vector & id_grid_to_solver, bool ac) const; virtual void fillpv(std::vector& bus_pv, std::vector & has_bus_been_added, - Eigen::VectorXi & slack_bus_id_solver, + const Eigen::VectorXi & slack_bus_id_solver, const std::vector & id_grid_to_solver) const; void init_q_vector(int nb_bus, Eigen::VectorXi & total_gen_per_bus, diff --git a/src/DataGeneric.cpp b/src/DataGeneric.cpp index 08d30ed..3b7c330 100644 --- a/src/DataGeneric.cpp +++ b/src/DataGeneric.cpp @@ -86,8 +86,8 @@ void DataGeneric::_change_bus(int el_id, int new_bus_me_id, Eigen::VectorXi & el // TODO speed: sparsity pattern might not change if something is already there solver_control.tell_ybus_change_sparsity_pattern(); - solver_control.tell_recompute_sbus(); - solver_control.tell_recompute_ybus(); + solver_control.tell_recompute_sbus(); // if a bus changed for load / generator + solver_control.tell_recompute_ybus(); // if a bus changed for shunts / line / trafo } bus_me_id = new_bus_me_id; } diff --git a/src/DataGeneric.h b/src/DataGeneric.h index d844f4a..768cc7d 100644 --- a/src/DataGeneric.h +++ b/src/DataGeneric.h @@ -88,7 +88,7 @@ class DataGeneric : public BaseConstants virtual void fillSbus(CplxVect & Sbus, const std::vector & id_grid_to_solver, bool ac) const {}; virtual void fillpv(std::vector& bus_pv, std::vector & has_bus_been_added, - Eigen::VectorXi & slack_bus_id_solver, + const Eigen::VectorXi & slack_bus_id_solver, const std::vector & id_grid_to_solver) const {}; virtual void get_q(std::vector& q_by_bus) {}; diff --git a/src/DataLine.h b/src/DataLine.h index 2026153..10b4be0 100644 --- a/src/DataLine.h +++ b/src/DataLine.h @@ -192,6 +192,7 @@ class DataLine : public DataGeneric if(status_[powerline_id]){ solver_control.tell_recompute_ybus(); // but sparsity pattern do not change here (possibly one more coeff at 0.) + solver_control.tell_ybus_some_coeffs_zero(); } _deactivate(powerline_id, status_); } diff --git a/src/DataShunt.cpp b/src/DataShunt.cpp index 6813a7a..ba52bbd 100644 --- a/src/DataShunt.cpp +++ b/src/DataShunt.cpp @@ -170,7 +170,10 @@ void DataShunt::change_p(int shunt_id, real_type new_p, SolverControl & solver_c { bool my_status = status_.at(shunt_id); // and this check that load_id is not out of bound if(!my_status) throw std::runtime_error("Impossible to change the active value of a disconnected shunt"); - if(p_mw_(shunt_id) != new_p) solver_control.tell_recompute_sbus(); + if(p_mw_(shunt_id) != new_p){ + solver_control.tell_recompute_ybus(); + solver_control.tell_recompute_sbus(); // in dc mode sbus is modified + } p_mw_(shunt_id) = new_p; } @@ -179,7 +182,9 @@ void DataShunt::change_q(int shunt_id, real_type new_q, SolverControl & solver_c { bool my_status = status_.at(shunt_id); // and this check that load_id is not out of bound if(!my_status) throw std::runtime_error("Impossible to change the reactive value of a disconnected shunt"); - if(q_mvar_(shunt_id) != new_q) solver_control.tell_recompute_sbus(); + if(q_mvar_(shunt_id) != new_q){ + solver_control.tell_recompute_ybus(); + } q_mvar_(shunt_id) = new_q; } diff --git a/src/DataTrafo.h b/src/DataTrafo.h index 90e05e9..033a58a 100644 --- a/src/DataTrafo.h +++ b/src/DataTrafo.h @@ -178,6 +178,7 @@ class DataTrafo : public DataGeneric if(status_[trafo_id]){ solver_control.tell_recompute_ybus(); // but sparsity pattern do not change here (possibly one more coeff at 0.) + solver_control.tell_ybus_some_coeffs_zero(); } _deactivate(trafo_id, status_); } diff --git a/src/GridModel.cpp b/src/GridModel.cpp index f4768f2..7bdbea3 100644 --- a/src/GridModel.cpp +++ b/src/GridModel.cpp @@ -260,7 +260,7 @@ void GridModel::reset(bool reset_solver, bool reset_ac, bool reset_dc) Ybus_dc_ = Eigen::SparseMatrix(); } - Sbus_ = CplxVect(); + acSbus_ = CplxVect(); dcSbus_ = CplxVect(); bus_pv_ = Eigen::VectorXi(); bus_pq_ = Eigen::VectorXi(); @@ -291,9 +291,14 @@ CplxVect GridModel::ac_pf(const CplxVect & Vinit, bool conv = false; CplxVect res = CplxVect(); + reset_results(); // reset the results + // pre process the data to define a proper jacobian matrix, the proper voltage vector etc. bool is_ac = true; - CplxVect V = pre_process_solver(Vinit, Ybus_ac_, + // std::cout << "before pre process" << std::endl; + CplxVect V = pre_process_solver(Vinit, + acSbus_, + Ybus_ac_, id_me_to_ac_solver_, id_ac_solver_to_me_, slack_bus_id_ac_solver_, @@ -301,10 +306,20 @@ CplxVect GridModel::ac_pf(const CplxVect & Vinit, solver_control_); // start the solver - slack_weights_ = generators_.get_slack_weights(Ybus_ac_.rows(), id_me_to_ac_solver_); - conv = _solver.compute_pf(Ybus_ac_, V, Sbus_, slack_bus_id_ac_solver_, slack_weights_, bus_pv_, bus_pq_, max_iter, tol / sn_mva_); + // std::cout << "before get_slack_weights" << std::endl; + if(solver_control_.need_reset_solver() || + solver_control_.has_dimension_changed() || + solver_control_.has_slack_participate_changed() || + solver_control_.has_pv_changed() || + solver_control_.has_slack_weight_changed()){ + // std::cout << "get_slack_weights" << std::endl; + slack_weights_ = generators_.get_slack_weights(Ybus_ac_.rows(), id_me_to_ac_solver_); + } + // std::cout << "before compute_pf" << std::endl; + conv = _solver.compute_pf(Ybus_ac_, V, acSbus_, slack_bus_id_ac_solver_, slack_weights_, bus_pv_, bus_pq_, max_iter, tol / sn_mva_); - // store results (in ac mode) + // store results (in ac mode) + // std::cout << "before process_results" << std::endl; process_results(conv, res, Vinit, true, id_me_to_ac_solver_); // return the vector of complex voltage at each bus @@ -377,14 +392,18 @@ CplxVect GridModel::check_solution(const CplxVect & V_proposed, bool check_q_lim bool is_ac = true; SolverControl reset_solver; reset_solver.tell_none_changed(); // TODO reset solver - CplxVect V = pre_process_solver(V_proposed, Ybus_ac_, - id_me_to_ac_solver_, id_ac_solver_to_me_, slack_bus_id_ac_solver_, + CplxVect V = pre_process_solver(V_proposed, + acSbus_, + Ybus_ac_, + id_me_to_ac_solver_, + id_ac_solver_to_me_, + slack_bus_id_ac_solver_, is_ac, reset_solver); // compute the mismatch CplxVect tmp = Ybus_ac_ * V; // this is a vector tmp = tmp.array().conjugate(); // i take the conjugate - auto mis = V.array() * tmp.array() - Sbus_.array(); + auto mis = V.array() * tmp.array() - acSbus_.array(); // TODO ac or dc here // store results CplxVect res = _get_results_back_to_orig_nodes(mis, @@ -406,6 +425,7 @@ CplxVect GridModel::check_solution(const CplxVect & V_proposed, bool check_q_lim }; CplxVect GridModel::pre_process_solver(const CplxVect & Vinit, + CplxVect & Sbus, Eigen::SparseMatrix & Ybus, std::vector & id_me_to_solver, std::vector & id_solver_to_me, @@ -416,44 +436,56 @@ CplxVect GridModel::pre_process_solver(const CplxVect & Vinit, // 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); - if(solver_control.need_reset_solver()) _solver.reset(); + if(solver_control.need_reset_solver()){ + // std::cout << "_ac_solver.reset();" << std::endl; + _solver.reset(); + } } else { _dc_solver.tell_solver_control(solver_control); if(solver_control.need_reset_solver()){ - std::cout << "_dc_solver.reset();" << std::endl; + // std::cout << "_dc_solver.reset();" << std::endl; _dc_solver.reset(); } } if (solver_control.need_reset_solver() || + solver_control.has_dimension_changed() || solver_control.has_slack_participate_changed()){ - std::cout << "get_slack_bus_id;" << std::endl; - slack_bus_id_ = generators_.get_slack_bus_id(); + slack_bus_id_solver = 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 + + // std::cout << "slack_bus_id_solver 1: "; + // for(auto el: slack_bus_id_solver) std::cout << el << ", "; + // std::cout << std::endl; } 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; + // 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; + // 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; + // init Sbus + Sbus = CplxVect::Constant(id_solver_to_me.size(), 0.); + // std::cout << "init_Sbus;" << std::endl; } if (solver_control.need_reset_solver() || + solver_control.has_dimension_changed() || 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); fillpv_pq(id_me_to_solver, id_solver_to_me, slack_bus_id_solver, solver_control); - std::cout << "fillpv_pq;" << std::endl; + // std::cout << "fillpv_pq;" << std::endl; } if (is_ac && (solver_control.need_reset_solver() || @@ -465,14 +497,15 @@ CplxVect GridModel::pre_process_solver(const CplxVect & Vinit, 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_); - std::cout << "total_gen_per_bus_;" << std::endl; + // std::cout << "total_gen_per_bus_;" << std::endl; } if (solver_control.need_reset_solver() || + solver_control.has_dimension_changed() || 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; + fillSbus_me(Sbus, is_ac, id_me_to_solver); + // std::cout << "fillSbus_me;" << std::endl; } const int nb_bus_solver = static_cast(id_solver_to_me.size()); @@ -484,7 +517,10 @@ 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; + // std::cout << "pre_process_solver: V result: "< & Ybus, id_me_to_solver = std::vector(nb_bus_init, _deactivated_bus_id); // by default, if a bus is disconnected, then it has a -1 there id_solver_to_me = std::vector(); id_solver_to_me.reserve(nb_bus_init); - int bus_id_solver=0; + int bus_id_solver = 0; for(int bus_id_me=0; bus_id_me < nb_bus_init; ++bus_id_me){ if(bus_status_[bus_id_me]){ // bus is connected @@ -552,31 +589,56 @@ void GridModel::init_Ybus(Eigen::SparseMatrix & Ybus, ++bus_id_solver; } } - const int nb_bus = static_cast(id_solver_to_me.size()); + const int nb_bus_solver = static_cast(id_solver_to_me.size()); - Ybus = Eigen::SparseMatrix(nb_bus, nb_bus); - Ybus.reserve(nb_bus + 2*powerlines_.nb() + 2*trafos_.nb()); + Ybus = Eigen::SparseMatrix(nb_bus_solver, nb_bus_solver); + Ybus.reserve(nb_bus_solver + 2*powerlines_.nb() + 2*trafos_.nb()); } -void GridModel::init_Sbus(CplxVect & Sbus, - std::vector& id_me_to_solver, - std::vector& id_solver_to_me, - Eigen::VectorXi & slack_bus_id_solver){ +void GridModel::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 int nb_bus = static_cast(id_solver_to_me.size()); - Sbus = CplxVect::Constant(nb_bus, 0.); - slack_bus_id_solver = Eigen::VectorXi::Zero(slack_bus_id_.size()); + 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); size_t i = 0; - for(auto el: slack_bus_id_) { - slack_bus_id_solver(i) = id_me_to_solver[el]; + // std::cout << "slack_bus_id_solver 2: "; + // for(auto el: slack_bus_id_solver) std::cout << el << ", "; + // std::cout << std::endl; + // std::cout << "id_me_to_solver: "; + // for(auto el: id_me_to_solver) std::cout << el << ", "; + // std::cout << std::endl; + // std::cout << "id_solver_to_me: "; + // for(auto el: id_solver_to_me) std::cout << el << ", "; + // std::cout << std::endl; + + // for(auto el: slack_bus_id_) { + for(auto el: slack_bus_id_solver) { + auto tmp = id_me_to_solver[el]; + if(tmp == _deactivated_bus_id){ + std::ostringstream exc_; + exc_ << "GridModel::init_Sbus: One of the slack bus is disconnected."; + exc_ << " You can check element "; + exc_ << el; + exc_ << ": ["; + for(auto el2 : slack_bus_id_solver) exc_ << el2 << ", "; + exc_ << "]."; + throw std::out_of_range(exc_.str()); + } + slack_bus_id_solver(i) = tmp; ++i; } + // std::cout << "slack_bus_id_solver 3: "; + // for(auto el: slack_bus_id_solver) std::cout << el << ", "; + // std::cout << std::endl; if(is_in_vect(_deactivated_bus_id, slack_bus_id_solver)){ // TODO improve error message with the gen_id // TODO DEBUG MODE: only check that in debug mode - throw std::runtime_error("One of the slack bus is disconnected !"); + throw std::runtime_error("GridModel::init_Sbus: One of the slack bus is disconnected !"); } } void GridModel::fillYbus(Eigen::SparseMatrix & res, bool ac, const std::vector& id_me_to_solver){ @@ -586,6 +648,7 @@ void GridModel::fillYbus(Eigen::SparseMatrix & res, bool ac, const st **/ // init the Ybus matrix + res.setZero(); // it should not be needed but might not hurt too much either. std::vector > tripletList; tripletList.reserve(bus_vn_kv_.size() + 4*powerlines_.nb() + 4*trafos_.nb() + shunts_.nb()); powerlines_.fillYbus(tripletList, ac, id_me_to_solver, sn_mva_); // TODO have a function to dispatch that to all type of elements @@ -596,13 +659,14 @@ void GridModel::fillYbus(Eigen::SparseMatrix & res, bool ac, const st storages_.fillYbus(tripletList, ac, id_me_to_solver, sn_mva_); generators_.fillYbus(tripletList, ac, id_me_to_solver, sn_mva_); dc_lines_.fillYbus(tripletList, ac, id_me_to_solver, sn_mva_); - res.setFromTriplets(tripletList.begin(), tripletList.end()); + res.setFromTriplets(tripletList.begin(), tripletList.end()); // works because "The initial contents of *this is destroyed" res.makeCompressed(); } void GridModel::fillSbus_me(CplxVect & Sbus, bool ac, const std::vector& id_me_to_solver) { - // init the Sbus vector + // init the Sbus + Sbus.array() = 0.; // reset to 0. powerlines_.fillSbus(Sbus, id_me_to_solver, ac); // TODO have a function to dispatch that to all type of elements trafos_.fillSbus(Sbus, id_me_to_solver, ac); shunts_.fillSbus(Sbus, id_me_to_solver, ac); @@ -617,8 +681,8 @@ void GridModel::fillSbus_me(CplxVect & Sbus, bool ac, const std::vector& id } void GridModel::fillpv_pq(const std::vector& id_me_to_solver, - std::vector& id_solver_to_me, - Eigen::VectorXi & slack_bus_id_solver, + const std::vector& id_solver_to_me, + const Eigen::VectorXi & slack_bus_id_solver, const SolverControl & solver_control) { // Nothing to do if neither pv, nor pq nor the dimension of the problem has changed @@ -633,8 +697,6 @@ void GridModel::fillpv_pq(const std::vector& id_me_to_solver, bus_pv.reserve(nb_bus); std::vector has_bus_been_added(nb_bus, false); - // std::cout << "id_me_to_solver.size(): " << id_me_to_solver.size() << std::endl; - bus_pv_ = Eigen::VectorXi(); bus_pq_ = Eigen::VectorXi(); powerlines_.fillpv(bus_pv, has_bus_been_added, slack_bus_id_solver, id_me_to_solver); // TODO have a function to dispatch that to all type of elements @@ -652,8 +714,8 @@ void GridModel::fillpv_pq(const std::vector& id_me_to_solver, bus_pq.push_back(bus_id); has_bus_been_added[bus_id] = true; // don't add it a second time } - bus_pv_ = Eigen::Map(bus_pv.data(), bus_pv.size()); - bus_pq_ = Eigen::Map(bus_pq.data(), bus_pq.size()); + bus_pv_ = Eigen::VectorXi::Map(&bus_pv[0], bus_pv.size()); + bus_pq_ = Eigen::VectorXi::Map(&bus_pq[0], bus_pq.size()); } void GridModel::compute_results(bool ac){ // retrieve results from powerflow @@ -681,12 +743,12 @@ void GridModel::compute_results(bool ac){ //handle_slack_bus active power CplxVect mismatch; // power mismatch at each bus (SOLVER BUS !!!) - RealVect ractive_mismatch; // not used in dc mode (DO NOT ATTEMPT TO USE IT THERE) + RealVect reactive_mismatch; // not used in dc mode (DO NOT ATTEMPT TO USE IT THERE) RealVect active_mismatch; if(ac){ // In AC mode i am not forced to run through all the grid auto tmp = (Ybus_ac_ * V).conjugate(); - mismatch = V.array() * tmp.array() - Sbus_.array(); + mismatch = V.array() * tmp.array() - acSbus_.array(); active_mismatch = mismatch.real() * sn_mva_; } else{ active_mismatch = RealVect::Zero(V.size()); @@ -697,13 +759,15 @@ void GridModel::compute_results(bool ac){ const auto id_slack = slack_bus_id_dc_solver_(0); active_mismatch(id_slack) = -dcSbus_.real().sum() * sn_mva_; } + // for(auto el: active_mismatch) std::cout << el << ", "; + // std::cout << std::endl; generators_.set_p_slack(active_mismatch, id_me_to_solver); - if(ac) ractive_mismatch = mismatch.imag() * sn_mva_; + if(ac) reactive_mismatch = mismatch.imag() * sn_mva_; // mainly to initialize the Q value of the generators in dc (just fill it with 0.) - generators_.set_q(ractive_mismatch, id_me_to_solver, ac, + generators_.set_q(reactive_mismatch, id_me_to_solver, ac, total_gen_per_bus_, total_q_min_per_bus_, total_q_max_per_bus_); - dc_lines_.set_q(ractive_mismatch, id_me_to_solver, ac, + dc_lines_.set_q(reactive_mismatch, id_me_to_solver, ac, total_gen_per_bus_, total_q_min_per_bus_, total_q_max_per_bus_); } @@ -739,22 +803,39 @@ CplxVect GridModel::dc_pf(const CplxVect & Vinit, bool conv = false; CplxVect res = CplxVect(); + reset_results(); // reset the results + // pre process the data to define a proper jacobian matrix, the proper voltage vector etc. bool is_ac = false; - CplxVect V = pre_process_solver(Vinit, Ybus_dc_, - id_me_to_dc_solver_, id_dc_solver_to_me_, slack_bus_id_dc_solver_, + CplxVect V = pre_process_solver(Vinit, + dcSbus_, + 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"; + // std::cout << "after pre proces\n"; // start the solver - if(solver_control_.has_slack_participate_changed() || + if(solver_control_.need_reset_solver() || + solver_control_.has_dimension_changed() || + 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"; + solver_control_.has_slack_weight_changed()){ + // TODO smarter solver: this is done both in ac and in dc ! + // std::cout << "get_slack_weights" << std::endl; + slack_weights_ = generators_.get_slack_weights(Ybus_dc_.rows(), id_me_to_dc_solver_); + } + // std::cout << "V (init to dc pf)\n"; + // for(auto el: V) std::cout << el << ", "; + // std::cout << std::endl; + // std::cout << "dcSbus (init to dc pf)\n"; + // for(auto el: dcSbus_) std::cout << el << ", "; + // std::cout << std::endl; 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"; + // 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"; + // std::cout << "after compute_pf\n"; return res; } @@ -1047,7 +1128,9 @@ std::tuple GridModel::assign_slack_to_most_connected(){ generators_.remove_all_slackbus(); res_gen_id = generators_.assign_slack_bus(res_bus_id, gen_p_per_bus, solver_control_); std::get<1>(res) = res_gen_id; - slack_bus_id_ = std::vector(); + // slack_bus_id_ = std::vector(); + slack_bus_id_ac_solver_ = Eigen::VectorXi(); + slack_bus_id_dc_solver_ = Eigen::VectorXi(); slack_weights_ = RealVect(); return res; } diff --git a/src/GridModel.h b/src/GridModel.h index 0030c4d..2bf85b7 100644 --- a/src/GridModel.h +++ b/src/GridModel.h @@ -255,7 +255,7 @@ class GridModel : public DataGeneric unsigned int size_th = 6; if (my_state.size() != size_th) { - std::cout << "LightSim::GridModel state size " << my_state.size() << " instead of "<< size_th << std::endl; + // std::cout << "LightSim::GridModel state size " << my_state.size() << " instead of "<< size_th << std::endl; // TODO more explicit error message throw std::runtime_error("Invalid state when loading LightSim::GridModel"); } @@ -471,7 +471,7 @@ class GridModel : public DataGeneric return Ybus_dc_; // This is copied to python } Eigen::Ref get_Sbus() const{ - return Sbus_; + return acSbus_; } Eigen::Ref get_dcSbus() const{ return dcSbus_; @@ -629,23 +629,30 @@ class GridModel : public DataGeneric if you will perform a powerflow after it or not. (usually put ``true`` here). **/ CplxVect pre_process_solver(const CplxVect & Vinit, + CplxVect & Sbus, Eigen::SparseMatrix & Ybus, std::vector & id_me_to_solver, std::vector & id_solver_to_me, Eigen::VectorXi & slack_bus_id_solver, bool is_ac, const SolverControl & solver_control); + + // init the Ybus matrix (its size, it is filled up elsewhere) and also the + // converter from "my bus id" to the "solver bus id" (id_me_to_solver and id_solver_to_me) void init_Ybus(Eigen::SparseMatrix & Ybus, std::vector & id_me_to_solver, std::vector& id_solver_to_me); - void init_Sbus(CplxVect & Sbus, - std::vector & id_me_to_solver, - std::vector& id_solver_to_me, - Eigen::VectorXi & slack_bus_id_solver); + + // converts the slack_bus_id from gridmodel ordering into solver ordering + 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); 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, std::vector& id_solver_to_me, - Eigen::VectorXi & slack_bus_id_solver, + void fillpv_pq(const std::vector& id_me_to_solver, + const std::vector& id_solver_to_me, + const Eigen::VectorXi & slack_bus_id_solver, const SolverControl & solver_control); // results @@ -794,7 +801,7 @@ class GridModel : public DataGeneric DataDCLine dc_lines_; // 8. slack bus - std::vector slack_bus_id_; + // std::vector slack_bus_id_; Eigen::VectorXi slack_bus_id_ac_solver_; Eigen::VectorXi slack_bus_id_dc_solver_; RealVect slack_weights_; @@ -802,7 +809,7 @@ class GridModel : public DataGeneric // as matrix, for the solver Eigen::SparseMatrix Ybus_ac_; Eigen::SparseMatrix Ybus_dc_; - CplxVect Sbus_; + CplxVect acSbus_; CplxVect dcSbus_; Eigen::VectorXi bus_pv_; // id are the solver internal id and NOT the initial id Eigen::VectorXi bus_pq_; // id are the solver internal id and NOT the initial id diff --git a/src/Utils.cpp b/src/Utils.cpp new file mode 100644 index 0000000..835bbb8 --- /dev/null +++ b/src/Utils.cpp @@ -0,0 +1,49 @@ +// Copyright (c) 2020, RTE (https://www.rte-france.com) +// See AUTHORS.txt +// This Source Code Form is subject to the terms of the Mozilla Public License, version 2.0. +// If a copy of the Mozilla Public License, version 2.0 was not distributed with this file, +// 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. + +#include "Utils.h" + +std::ostream& operator<<(std::ostream& out, const ErrorType & error_type){ + switch (error_type) + { + case ErrorType::NoError: + out << "NoError"; + break; + case ErrorType::SingularMatrix: + out << "SingularMatrix"; + break; + case ErrorType::TooManyIterations: + out << "TooManyIterations"; + break; + case ErrorType::InifiniteValue: + out << "InifiniteValue"; + break; + case ErrorType::SolverAnalyze: + out << "SolverAnalyze"; + break; + case ErrorType::SolverFactor: + out << "SolverFactor"; + break; + case ErrorType::SolverReFactor: + out << "SolverReFactor"; + break; + case ErrorType::SolverSolve: + out << "SolverSolve"; + break; + case ErrorType::NotInitError: + out << "NotInitError"; + break; + case ErrorType::LicenseError: + out << "LicenseError"; + break; + default: + out << "unknown error (check utils.cpp)"; + break; + } + return out; +} diff --git a/src/Utils.h b/src/Utils.h index 0d56bc0..406c125 100644 --- a/src/Utils.h +++ b/src/Utils.h @@ -37,7 +37,17 @@ typedef Eigen::Matrix RealMat; typedef Eigen::Matrix CplxMat; // type of error in the different solvers -enum class ErrorType {NoError, SingularMatrix, TooManyIterations, InifiniteValue, SolverAnalyze, SolverFactor, SolverReFactor, SolverSolve, NotInitError, LicenseError}; +enum class ErrorType {NoError, + SingularMatrix, + TooManyIterations, + InifiniteValue, + SolverAnalyze, + SolverFactor, + SolverReFactor, + SolverSolve, + NotInitError, + LicenseError}; +std::ostream& operator<<(std::ostream& out, const ErrorType & error_type); // define some constant for compilation outside of "setup.py" #ifndef VERSION_MAJOR @@ -65,6 +75,7 @@ class SolverControl need_recompute_ybus_(true), v_changed_(true), slack_weight_changed_(true), + ybus_some_coeffs_zero_(true), ybus_change_sparsity_pattern_(true) {}; @@ -78,6 +89,7 @@ class SolverControl need_recompute_ybus_ = true; v_changed_ = true; slack_weight_changed_ = true; + ybus_some_coeffs_zero_ = true; ybus_change_sparsity_pattern_ = true; } @@ -91,6 +103,7 @@ class SolverControl need_recompute_ybus_ = false; v_changed_ = false; slack_weight_changed_ = false; + ybus_some_coeffs_zero_ = false; ybus_change_sparsity_pattern_ = false; } @@ -114,6 +127,10 @@ class SolverControl void tell_v_changed(){v_changed_ = true;} // at least one generator has changed its slack participation void tell_slack_weight_changed(){slack_weight_changed_ = true;} + // tells that some coeff of ybus might have been set to 0. + // (and ybus compressed again, so these coeffs are really completely hidden) + // might need to trigger some recomputation of some solvers (eg NR based ones) + void tell_ybus_some_coeffs_zero(){ybus_some_coeffs_zero_ = true;} bool has_dimension_changed() const {return change_dimension_;} bool has_pv_changed() const {return pv_changed_;} @@ -125,6 +142,7 @@ class SolverControl bool ybus_change_sparsity_pattern() const {return ybus_change_sparsity_pattern_;} bool has_slack_weight_changed() const {return slack_weight_changed_;} bool has_v_changed() const {return v_changed_;} + bool has_ybus_some_coeffs_zero() const {return ybus_some_coeffs_zero_;} protected: bool change_dimension_; @@ -136,6 +154,7 @@ class SolverControl bool need_recompute_ybus_; // some coeff of ybus changed, but not its sparsity pattern bool v_changed_; bool slack_weight_changed_; + bool ybus_some_coeffs_zero_; // tells that some coeff of ybus might have been set to 0. (and ybus compressed again, so these coeffs are really completely hidden) bool ybus_change_sparsity_pattern_; // sparsity pattern of ybus changed (and so are its coeff), or ybus change of dimension };