diff --git a/CHANGELOG.rst b/CHANGELOG.rst index 025c310..4a4113e 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -32,11 +32,15 @@ Change Log these cases. - [FIXED] a bug leading to not propagate correctly the "compute_results" flag when the environment was copied (for example) +- [FIXED] a bug where copying a lightsim2grid `GridModel` did not fully copy it +- [FIXED] a bug in the "topo_vect" comprehension cpp side (sometimes some buses + might not be activated / deactivated correctly) - [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. - [ADDED] embed in the generator models the "non pv" behaviour. (TODO need to be able to change Q from python side) - [ADDED] computation of PTPF (Power Transfer Distribution Factor) is now possible +- [ADDED] (not tested) support for more than 2 busbars per substation - [IMPROVED] now performing the new grid2op `create_test_suite` - [IMPROVED] now lightsim2grid properly throw `BackendError` - [IMPROVED] clean ce cpp side by refactoring: making clearer the difference (linear) solver diff --git a/lightsim2grid/tests/test_solver_control.py b/lightsim2grid/tests/test_solver_control.py index aa1a441..4a1c413 100644 --- a/lightsim2grid/tests/test_solver_control.py +++ b/lightsim2grid/tests/test_solver_control.py @@ -25,6 +25,8 @@ import numpy as np import grid2op from grid2op.Action import CompleteAction +from grid2op.Parameters import Parameters + from lightsim2grid import LightSimBackend from lightsim2grid.solver import SolverType @@ -50,6 +52,41 @@ def setUp(self) -> None: self.tol_solver = 1e-8 # solver self.tol_equal = 1e-10 # for comparing with and without the "smarter solver" things, and make sure everything is really equal! + def test_update_topo_ac(self, runpf_fun="_run_ac_pf"): + """test when I disconnect a line alone at one end: it changes the size of the ybus / sbus vector AC""" + LINE_ID = 2 + dim_topo = type(self.env).dim_topo + mask_changed = np.zeros(dim_topo, dtype=bool) + mask_val = np.zeros(dim_topo, dtype=np.int32) + mask_changed[type(self.env).line_ex_pos_topo_vect[LINE_ID]] = True + mask_val[type(self.env).line_ex_pos_topo_vect[LINE_ID]] = 2 + self.gridmodel.update_topo(mask_changed, mask_val) + V = getattr(self, runpf_fun)(gridmodel=self.gridmodel) + assert len(V), "it should not have diverged here" + self.gridmodel.unset_changes() + + mask_changed = np.zeros(dim_topo, dtype=bool) + mask_val = np.zeros(dim_topo, dtype=np.int32) + mask_changed[type(self.env).line_or_pos_topo_vect[LINE_ID]] = True + mask_val[type(self.env).line_or_pos_topo_vect[LINE_ID]] = -1 + # mask_changed[type(self.env).line_ex_pos_topo_vect[LINE_ID]] = True + # mask_val[type(self.env).line_ex_pos_topo_vect[LINE_ID]] = -1 + self.gridmodel.update_topo(mask_changed, mask_val) + solver = self.gridmodel.get_dc_solver() if runpf_fun == "_run_dc_pf" else self.gridmodel.get_solver() + V1 = getattr(self, runpf_fun)(gridmodel=self.gridmodel) + assert len(V1), f"it should not have diverged here. Error : {solver.get_error()}" + + self.gridmodel.tell_solver_need_reset() + V2 = getattr(self, runpf_fun)(gridmodel=self.gridmodel) + assert len(V2), f"it should not have diverged here. Error : {solver.get_error()}" + assert np.allclose(V1, V2, rtol=self.tol_equal, atol=self.tol_equal) + + def test_update_topo_dc(self): + """test when I disconnect a line alone at one end: it changes the size of the ybus / sbus vector AC""" + if not self.need_dc: + self.skipTest("Useless to run DC") + self.test_update_topo_ac("_run_dc_pf") + def test_pf_run_dc(self): """test I have the same results if nothing is done with and without restarting from scratch when running dc powerflow""" if not self.need_dc: @@ -138,7 +175,7 @@ def aux_do_undo_ac(self, getattr(self, funname_undo)(gridmodel=self.gridmodel, el_id=el_id, el_val=el_val) V_reco = getattr(self, runpf_fun)(gridmodel=self.gridmodel) assert len(V_reco), f"error for el_id={el_id}: gridmodel should converge in {pf_mode}" - assert np.allclose(V_reco, V_init, rtol=self.tol_equal, atol=self.tol_equal), f"error for el_id={el_id}: do an action and then undo it should not have any impact in {pf_mode}" + assert np.allclose(V_reco, V_init, rtol=self.tol_equal, atol=self.tol_equal), f"error for el_id={el_id}: do an action and then undo it should not have any impact in {pf_mode}: max {np.abs(V_init - V_reco).max():.2e}" self.gridmodel.unset_changes() V_reco1 = getattr(self, runpf_fun)(gridmodel=self.gridmodel) assert len(V_reco1), f"error for el_id={el_id}: gridmodel should converge in {pf_mode}" diff --git a/src/GridModel.cpp b/src/GridModel.cpp index 5f4579a..ba83336 100644 --- a/src/GridModel.cpp +++ b/src/GridModel.cpp @@ -14,6 +14,7 @@ GridModel::GridModel(const GridModel & other) { reset(true, true, true); + max_nb_bus_per_sub_ = other.max_nb_bus_per_sub_; init_vm_pu_ = other.init_vm_pu_; sn_mva_ = other.sn_mva_; @@ -97,6 +98,22 @@ GridModel::StateRes GridModel::get_state() const auto res_storage = storages_.get_state(); auto res_dc_line = dc_lines_.get_state(); + std::vector load_pos_topo_vect(load_pos_topo_vect_.begin(), load_pos_topo_vect_.end()); + std::vector gen_pos_topo_vect(gen_pos_topo_vect_.begin(), gen_pos_topo_vect_.end()); + std::vector line_or_pos_topo_vect(line_or_pos_topo_vect_.begin(), line_or_pos_topo_vect_.end()); + std::vector line_ex_pos_topo_vect(line_ex_pos_topo_vect_.begin(), line_ex_pos_topo_vect_.end()); + std::vector trafo_hv_pos_topo_vect(trafo_hv_pos_topo_vect_.begin(), trafo_hv_pos_topo_vect_.end()); + std::vector trafo_lv_pos_topo_vect(trafo_lv_pos_topo_vect_.begin(), trafo_lv_pos_topo_vect_.end()); + std::vector storage_pos_topo_vect(storage_pos_topo_vect_.begin(), storage_pos_topo_vect_.end()); + + std::vector load_to_subid(load_to_subid_.begin(), load_to_subid_.end()); + std::vector gen_to_subid(gen_to_subid_.begin(), gen_to_subid_.end()); + std::vector line_or_to_subid(line_or_to_subid_.begin(), line_or_to_subid_.end()); + std::vector line_ex_to_subid(line_ex_to_subid_.begin(), line_ex_to_subid_.end()); + std::vector trafo_hv_to_subid(trafo_hv_to_subid_.begin(), trafo_hv_to_subid_.end()); + std::vector trafo_lv_to_subid(trafo_lv_to_subid_.begin(), trafo_lv_to_subid_.end()); + std::vector storage_to_subid(storage_to_subid_.begin(), storage_to_subid_.end()); + GridModel::StateRes res(version_major, version_medium, version_minor, @@ -112,7 +129,23 @@ GridModel::StateRes GridModel::get_state() const res_load, res_sgen, res_storage, - res_dc_line + res_dc_line, + n_sub_, + max_nb_bus_per_sub_, + load_pos_topo_vect, + gen_pos_topo_vect, + line_or_pos_topo_vect, + line_ex_pos_topo_vect, + trafo_hv_pos_topo_vect, + trafo_lv_pos_topo_vect, + storage_pos_topo_vect, + load_to_subid, + gen_to_subid, + line_or_to_subid, + line_ex_to_subid, + trafo_hv_to_subid, + trafo_lv_to_subid, + storage_to_subid ); return res; }; @@ -139,11 +172,11 @@ void GridModel::set_state(GridModel::StateRes & my_state) exc_ << "It is not possible. Please reinstall it."; throw std::runtime_error(exc_.str()); } - std::vector ls_to_pp_ = std::get<3>(my_state); + const std::vector & ls_to_pp = std::get<3>(my_state); init_vm_pu_ = std::get<4>(my_state); sn_mva_ = std::get<5>(my_state); - std::vector & bus_vn_kv = std::get<6>(my_state); - std::vector & bus_status = std::get<7>(my_state); + const std::vector & bus_vn_kv = std::get<6>(my_state); + const std::vector & bus_status = std::get<7>(my_state); // powerlines LineContainer::StateRes & state_lines = std::get<8>(my_state); @@ -162,13 +195,45 @@ void GridModel::set_state(GridModel::StateRes & my_state) // dc lines DCLineContainer::StateRes & state_dc_lines = std::get<15>(my_state); - // assign it to this instance + // grid2op specific + n_sub_ = std::get<16>(my_state); + max_nb_bus_per_sub_ = std::get<17>(my_state); + const std::vector & load_pos_topo_vect = std::get<18>(my_state); + const std::vector & gen_pos_topo_vect = std::get<19>(my_state); + const std::vector & line_or_pos_topo_vect = std::get<20>(my_state); + const std::vector & line_ex_pos_topo_vect = std::get<21>(my_state); + const std::vector & trafo_hv_pos_topo_vect = std::get<22>(my_state); + const std::vector & trafo_lv_pos_topo_vect = std::get<23>(my_state); + const std::vector & storage_pos_topo_vect = std::get<24>(my_state); + const std::vector & load_to_subid = std::get<25>(my_state); + const std::vector & gen_to_subid = std::get<26>(my_state); + const std::vector & line_or_to_subid = std::get<27>(my_state); + const std::vector & line_ex_to_subid = std::get<28>(my_state); + const std::vector & trafo_hv_to_subid = std::get<29>(my_state); + const std::vector & trafo_lv_to_subid = std::get<30>(my_state); + const std::vector & storage_to_subid = std::get<31>(my_state); + + load_pos_topo_vect_ = IntVectRowMaj::Map(load_pos_topo_vect.data(), load_pos_topo_vect.size()); + gen_pos_topo_vect_ = IntVectRowMaj::Map(gen_pos_topo_vect.data(), gen_pos_topo_vect.size()); + line_or_pos_topo_vect_ = IntVectRowMaj::Map(line_or_pos_topo_vect.data(), line_or_pos_topo_vect.size()); + line_ex_pos_topo_vect_ = IntVectRowMaj::Map(line_ex_pos_topo_vect.data(), line_ex_pos_topo_vect.size()); + trafo_hv_pos_topo_vect_ = IntVectRowMaj::Map(trafo_hv_pos_topo_vect.data(), trafo_hv_pos_topo_vect.size()); + trafo_lv_pos_topo_vect_ = IntVectRowMaj::Map(trafo_lv_pos_topo_vect.data(), trafo_lv_pos_topo_vect.size()); + storage_pos_topo_vect_ = IntVectRowMaj::Map(storage_pos_topo_vect.data(), storage_pos_topo_vect.size()); + load_to_subid_ = IntVectRowMaj::Map(load_to_subid.data(), load_to_subid.size()); + gen_to_subid_ = IntVectRowMaj::Map(gen_to_subid.data(), gen_to_subid.size()); + line_or_to_subid_ = IntVectRowMaj::Map(line_or_to_subid.data(), line_or_to_subid.size()); + line_ex_to_subid_ = IntVectRowMaj::Map(line_ex_to_subid.data(), line_ex_to_subid.size()); + trafo_hv_to_subid_ = IntVectRowMaj::Map(trafo_hv_to_subid.data(), trafo_hv_to_subid.size()); + trafo_lv_to_subid_ = IntVectRowMaj::Map(trafo_lv_to_subid.data(), trafo_lv_to_subid.size()); + storage_to_subid_ = IntVectRowMaj::Map(storage_to_subid.data(), storage_to_subid.size()); - set_ls_to_orig(IntVect::Map(&ls_to_pp_[0], ls_to_pp_.size())); // set also _orig_to_ls + // assign it to this instance + set_ls_to_orig(IntVect::Map(ls_to_pp.data(), ls_to_pp.size())); // set also _orig_to_ls // buses // 1. bus_vn_kv_ - bus_vn_kv_ = RealVect::Map(&bus_vn_kv[0], bus_vn_kv.size()); + bus_vn_kv_ = RealVect::Map(bus_vn_kv.data(), bus_vn_kv.size()); // 2. bus status bus_status_ = bus_status; @@ -240,7 +305,8 @@ void GridModel::init_bus(const RealVect & bus_vn_kv, int nb_line, int nb_trafo){ and initialize the Ybus_ matrix at the proper shape **/ - const int nb_bus = static_cast(bus_vn_kv.size()); + const int nb_bus = static_cast(bus_vn_kv.size()); // size of buses are checked in set_max_nb_bus_per_sub + bus_vn_kv_ = bus_vn_kv; // base_kv bus_status_ = std::vector(nb_bus, true); // by default everything is connected @@ -997,9 +1063,6 @@ void GridModel::update_storages_p(Eigen::Ref > has_changed, Eigen::Ref > new_values) { - const int nb_bus = static_cast(bus_status_.size()); - for(int i = 0; i < nb_bus; ++i) bus_status_[i] = false; - update_topo_generic(has_changed, new_values, load_pos_topo_vect_, load_to_subid_, &GridModel::reactivate_load, @@ -1045,6 +1108,19 @@ void GridModel::update_topo(Eigen::Ref(bus_status_.size()); + for(int i = 0; i < nb_bus; ++i) bus_status_[i] = false; + + powerlines_.update_bus_status(bus_status_); // TODO have a function to dispatch that to all type of elements + shunts_.update_bus_status(bus_status_); + trafos_.update_bus_status(bus_status_); + loads_.update_bus_status(bus_status_); + sgens_.update_bus_status(bus_status_); + storages_.update_bus_status(bus_status_); + generators_.update_bus_status(bus_status_); + dc_lines_.update_bus_status(bus_status_); } // for FDPF (implementation of the alg 2 method FDBX (FDXB will follow) // TODO FDPF diff --git a/src/GridModel.h b/src/GridModel.h index 123a711..7a0a104 100644 --- a/src/GridModel.h +++ b/src/GridModel.h @@ -45,6 +45,8 @@ class GridModel : public GenericContainer { public: + typedef Eigen::Array IntVectRowMaj; + typedef std::tuple< int, // version major int, // version medium @@ -69,14 +71,32 @@ class GridModel : public GenericContainer // storage units LoadContainer::StateRes, //dc lines - DCLineContainer::StateRes + DCLineContainer::StateRes, + // grid2op specific + int, // n_sub + int, // max_nb_bus_per_sub + std::vector, // load_pos_topo_vect_ + std::vector, + std::vector, + std::vector, + std::vector, + std::vector, + std::vector, + std::vector, // load_to_subid_ + std::vector, + std::vector, + std::vector, + std::vector, + std::vector, + std::vector > StateRes; GridModel(): solver_control_(), compute_results_(true), init_vm_pu_(1.04), - sn_mva_(1.0){ + sn_mva_(1.0), + max_nb_bus_per_sub_(2){ _solver.change_solver(SolverType::SparseLU); _dc_solver.change_solver(SolverType::DC); _solver.set_gridmodel(this); @@ -591,6 +611,20 @@ class GridModel : public GenericContainer { n_sub_ = n_sub; } + void set_max_nb_bus_per_sub(int max_nb_bus_per_sub) + { + max_nb_bus_per_sub_ = max_nb_bus_per_sub; + if(bus_vn_kv_.size() != n_sub_ * max_nb_bus_per_sub_){ + std::ostringstream exc_; + exc_ << "GridModel::set_max_nb_bus_per_sub: "; + exc_ << "your model counts "; + exc_ << bus_vn_kv_.size() << " buses according to `bus_vn_kv_` but "; + exc_ << n_sub_ * max_nb_bus_per_sub_ << " according to n_sub_ * max_nb_bus_per_sub_."; + exc_ << "Both should match: either reinit it with another call to `init_bus` or set properly the number of "; + exc_ << "substations / buses per substations with `set_n_sub` / `set_max_nb_bus_per_sub`"; + throw std::runtime_error(exc_.str()); + } + } void fillSbus_other(CplxVect & res, bool ac, const std::vector& id_me_to_solver){ fillSbus_me(res, ac, id_me_to_solver); @@ -706,27 +740,23 @@ class GridModel : public GenericContainer { for(int el_id = 0; el_id < vect_pos.rows(); ++el_id) { + int el_pos = vect_pos(el_id); + if(! has_changed(el_pos)) continue; int new_bus = new_values(el_pos); if(new_bus > 0){ // new bus is a real bus, so i need to make sure to have it turned on, and then change the bus int init_bus_me = vect_subid(el_id); - int new_bus_backend = new_bus == 1 ? init_bus_me : init_bus_me + n_sub_ ; + int new_bus_backend = new_bus == 1 ? init_bus_me : init_bus_me + n_sub_ * (max_nb_bus_per_sub_ - 1); bus_status_[new_bus_backend] = true; - if(has_changed(el_pos)) - { - (this->*fun_react)(el_id); // eg reactivate_load(load_id); - (this->*fun_change)(el_id, new_bus_backend); // eg change_bus_load(load_id, new_bus_backend); - } + (this->*fun_react)(el_id); // eg reactivate_load(load_id); + (this->*fun_change)(el_id, new_bus_backend); // eg change_bus_load(load_id, new_bus_backend); } else{ - if(has_changed(el_pos)) - { - // new bus is negative, we deactivate it - (this->*fun_deact)(el_id);// eg deactivate_load(load_id); - // bus_status_ is set to "false" in GridModel.update_topo - // and a bus is activated if (and only if) one element is connected to it. - // I must not set `bus_status_[new_bus_backend] = false;` in this case ! - } + // new bus is negative, we deactivate it + (this->*fun_deact)(el_id);// eg deactivate_load(load_id); + // bus_status_ is set to "false" in GridModel.update_topo + // and a bus is activated if (and only if) one element is connected to it. + // I must not set `bus_status_[new_bus_backend] = false;` in this case ! } } } @@ -824,21 +854,22 @@ class GridModel : public GenericContainer // specific grid2op int n_sub_; - Eigen::Array load_pos_topo_vect_; - Eigen::Array gen_pos_topo_vect_; - Eigen::Array line_or_pos_topo_vect_; - Eigen::Array line_ex_pos_topo_vect_; - Eigen::Array trafo_hv_pos_topo_vect_; - Eigen::Array trafo_lv_pos_topo_vect_; - Eigen::Array storage_pos_topo_vect_; - - Eigen::Array load_to_subid_; - Eigen::Array gen_to_subid_; - Eigen::Array line_or_to_subid_; - Eigen::Array line_ex_to_subid_; - Eigen::Array trafo_hv_to_subid_; - Eigen::Array trafo_lv_to_subid_; - Eigen::Array storage_to_subid_; + int max_nb_bus_per_sub_; + IntVectRowMaj load_pos_topo_vect_; + IntVectRowMaj gen_pos_topo_vect_; + IntVectRowMaj line_or_pos_topo_vect_; + IntVectRowMaj line_ex_pos_topo_vect_; + IntVectRowMaj trafo_hv_pos_topo_vect_; + IntVectRowMaj trafo_lv_pos_topo_vect_; + IntVectRowMaj storage_pos_topo_vect_; + + IntVectRowMaj load_to_subid_; + IntVectRowMaj gen_to_subid_; + IntVectRowMaj line_or_to_subid_; + IntVectRowMaj line_ex_to_subid_; + IntVectRowMaj trafo_hv_to_subid_; + IntVectRowMaj trafo_lv_to_subid_; + IntVectRowMaj storage_to_subid_; }; diff --git a/src/element_container/DCLineContainer.h b/src/element_container/DCLineContainer.h index 14abf21..5ee14dc 100644 --- a/src/element_container/DCLineContainer.h +++ b/src/element_container/DCLineContainer.h @@ -189,7 +189,11 @@ class DCLineContainer : public GenericContainer virtual void get_graph(std::vector > & res) const {}; virtual void disconnect_if_not_in_main_component(std::vector & busbar_in_main_component); virtual void nb_line_end(std::vector & res) const; - + virtual void update_bus_status(std::vector & bus_status) const { + from_gen_.update_bus_status(bus_status); + to_gen_.update_bus_status(bus_status); + } + real_type get_qmin_or(int dcline_id) {return from_gen_.get_qmin(dcline_id);} real_type get_qmax_or(int dcline_id) {return from_gen_.get_qmax(dcline_id);} real_type get_qmin_ex(int dcline_id) {return to_gen_.get_qmin(dcline_id);} diff --git a/src/element_container/GeneratorContainer.cpp b/src/element_container/GeneratorContainer.cpp index 254813f..e45767b 100644 --- a/src/element_container/GeneratorContainer.cpp +++ b/src/element_container/GeneratorContainer.cpp @@ -432,7 +432,7 @@ void GeneratorContainer::set_q(const RealVect & reactive_mismatch, void GeneratorContainer::update_slack_weights(Eigen::Ref > could_be_slack, - SolverControl & solver_control) + SolverControl & solver_control) { const int nb_gen = nb(); for(int gen_id = 0; gen_id < nb_gen; ++gen_id) diff --git a/src/element_container/GeneratorContainer.h b/src/element_container/GeneratorContainer.h index 1c46327..2dbdb4c 100644 --- a/src/element_container/GeneratorContainer.h +++ b/src/element_container/GeneratorContainer.h @@ -175,14 +175,15 @@ class GeneratorContainer: public GenericContainer **/ void add_slackbus(int gen_id, real_type weight, SolverControl & solver_control){ // TODO DEBUG MODE - if(weight <= 0.) throw std::runtime_error("GeneratorContainer::add_slackbus Cannot assign a negative weight to the slack bus."); + if(weight <= 0.) throw std::runtime_error("GeneratorContainer::add_slackbus Cannot assign a negative (<=0) weight to the slack bus."); if(!gen_slackbus_[gen_id]) solver_control.tell_slack_participate_changed(); gen_slackbus_[gen_id] = true; if(gen_slack_weight_[gen_id] != weight) solver_control.tell_slack_weight_changed(); gen_slack_weight_[gen_id] = weight; } void remove_slackbus(int gen_id, SolverControl & solver_control){ - if(!gen_slackbus_[gen_id]) solver_control.tell_slack_participate_changed(); + if(gen_slackbus_[gen_id]) solver_control.tell_slack_participate_changed(); + if(gen_slack_weight_[gen_id] != 0.) solver_control.tell_slack_weight_changed(); gen_slackbus_[gen_id] = false; gen_slack_weight_[gen_id] = 0.; } @@ -243,7 +244,7 @@ class GeneratorContainer: public GenericContainer 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(); + solver_control.tell_pv_changed(); if(gen_slack_weight_[gen_id] != 0. || gen_slackbus_[gen_id]){ solver_control.tell_slack_participate_changed(); solver_control.tell_slack_weight_changed(); @@ -256,7 +257,7 @@ class GeneratorContainer: public GenericContainer 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(); + solver_control.tell_pv_changed(); if(gen_slack_weight_[gen_id] != 0. || gen_slackbus_[gen_id]){ solver_control.tell_slack_participate_changed(); solver_control.tell_slack_weight_changed(); @@ -272,7 +273,14 @@ class GeneratorContainer: public GenericContainer 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); - + virtual void update_bus_status(std::vector & bus_status) const { + const int nb_gen = nb(); + for(int gen_id = 0; gen_id < nb_gen; ++gen_id) + { + if(!status_[gen_id]) continue; + bus_status[bus_id_[gen_id]] = true; + } + } real_type get_qmin(int gen_id) {return min_q_.coeff(gen_id);} real_type get_qmax(int gen_id) {return max_q_.coeff(gen_id);} void change_p(int gen_id, real_type new_p, SolverControl & solver_control); diff --git a/src/element_container/GenericContainer.h b/src/element_container/GenericContainer.h index c0e6da3..d48ef00 100644 --- a/src/element_container/GenericContainer.h +++ b/src/element_container/GenericContainer.h @@ -91,7 +91,8 @@ class GenericContainer : public BaseConstants const std::vector & id_grid_to_solver) const {}; virtual void get_q(std::vector& q_by_bus) {}; - + virtual void update_bus_status(std::vector & bus_status) const {}; + void set_p_slack(const RealVect& node_mismatch, const std::vector & id_grid_to_solver) {}; static const int _deactivated_bus_id; diff --git a/src/element_container/LineContainer.h b/src/element_container/LineContainer.h index fb78817..7e6acf4 100644 --- a/src/element_container/LineContainer.h +++ b/src/element_container/LineContainer.h @@ -180,6 +180,15 @@ class LineContainer : public GenericContainer virtual void disconnect_if_not_in_main_component(std::vector & busbar_in_main_component); virtual void nb_line_end(std::vector & res) const; virtual void get_graph(std::vector > & res) const; + virtual void update_bus_status(std::vector & bus_status) const { + const int nb_ = nb(); + for(int el_id = 0; el_id < nb_; ++el_id) + { + if(!status_[el_id]) continue; + bus_status[bus_or_id_[el_id]] = true; + bus_status[bus_ex_id_[el_id]] = true; + } + } void deactivate(int powerline_id, SolverControl & solver_control) { // std::cout << "line: deactivate called\n"; @@ -187,6 +196,7 @@ class LineContainer : public GenericContainer 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(); + solver_control.tell_dimension_changed(); // if the extremity of the line is alone on a bus, this can happen... } _deactivate(powerline_id, status_); } @@ -194,6 +204,7 @@ class LineContainer : public GenericContainer if(!status_[powerline_id]){ solver_control.tell_recompute_ybus(); solver_control.tell_ybus_change_sparsity_pattern(); // sparsity pattern might change: a non zero coeff can pop up + solver_control.tell_dimension_changed(); // if the extremity of the line is alone on a bus, this can happen... } _reactivate(powerline_id, status_); } diff --git a/src/element_container/LoadContainer.h b/src/element_container/LoadContainer.h index 3786d52..0483f3f 100644 --- a/src/element_container/LoadContainer.h +++ b/src/element_container/LoadContainer.h @@ -158,6 +158,14 @@ class LoadContainer : public GenericContainer virtual void disconnect_if_not_in_main_component(std::vector & busbar_in_main_component); virtual void fillSbus(CplxVect & Sbus, const std::vector & id_grid_to_solver, bool ac) const; + virtual void update_bus_status(std::vector & bus_status) const { + const int nb_ = nb(); + for(int el_id = 0; el_id < nb_; ++el_id) + { + if(!status_[el_id]) continue; + bus_status[bus_id_[el_id]] = true; + } + } void compute_results(const Eigen::Ref & Va, const Eigen::Ref & Vm, diff --git a/src/element_container/SGenContainer.h b/src/element_container/SGenContainer.h index ae05cdb..b50a269 100644 --- a/src/element_container/SGenContainer.h +++ b/src/element_container/SGenContainer.h @@ -179,6 +179,14 @@ class SGenContainer: public GenericContainer virtual void fillSbus(CplxVect & Sbus, const std::vector & id_grid_to_solver, bool ac) const ; virtual void gen_p_per_bus(std::vector & res) const; + virtual void update_bus_status(std::vector & bus_status) const { + const int nb_ = nb(); + for(int el_id = 0; el_id < nb_; ++el_id) + { + if(!status_[el_id]) continue; + bus_status[bus_id_[el_id]] = true; + } + } void compute_results(const Eigen::Ref & Va, const Eigen::Ref & Vm, diff --git a/src/element_container/ShuntContainer.h b/src/element_container/ShuntContainer.h index 4fe808c..8458acf 100644 --- a/src/element_container/ShuntContainer.h +++ b/src/element_container/ShuntContainer.h @@ -161,6 +161,14 @@ class ShuntContainer : public GenericContainer FDPFMethod xb_or_bx) const; virtual void fillYbus_spmat(Eigen::SparseMatrix & res, bool ac, const std::vector & id_grid_to_solver); virtual void fillSbus(CplxVect & Sbus, const std::vector & id_grid_to_solver, bool ac) const; // in DC i need that + virtual void update_bus_status(std::vector & bus_status) const { + const int nb_ = nb(); + for(int el_id = 0; el_id < nb_; ++el_id) + { + if(!status_[el_id]) continue; + bus_status[bus_id_[el_id]] = true; + } + } void compute_results(const Eigen::Ref & Va, const Eigen::Ref & Vm, diff --git a/src/element_container/TrafoContainer.h b/src/element_container/TrafoContainer.h index 6e15025..fbb07f8 100644 --- a/src/element_container/TrafoContainer.h +++ b/src/element_container/TrafoContainer.h @@ -167,7 +167,7 @@ class TrafoContainer : public GenericContainer } if(id >= nb()) { - throw std::range_error("Generator out of bound. Not enough transformers on the grid."); + throw std::range_error("Trafo out of bound. Not enough transformers on the grid."); } return TrafoInfo(*this, id); } @@ -178,6 +178,7 @@ class TrafoContainer : public GenericContainer 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(); + solver_control.tell_dimension_changed(); // if the extremity of the line is alone on a bus, this can happen... } _deactivate(trafo_id, status_); } @@ -185,6 +186,7 @@ class TrafoContainer : public GenericContainer if(!status_[trafo_id]){ solver_control.tell_recompute_ybus(); solver_control.tell_ybus_change_sparsity_pattern(); // this might change + solver_control.tell_dimension_changed(); // if the extremity of the line is alone on a bus, this can happen... } _reactivate(trafo_id, status_); } @@ -214,6 +216,15 @@ class TrafoContainer : public GenericContainer int nb_powerline, bool transpose) const; virtual void hack_Sbus_for_dc_phase_shifter(CplxVect & Sbus, bool ac, const std::vector & id_grid_to_solver); // needed for dc mode + virtual void update_bus_status(std::vector & bus_status) const { + const int nb_ = nb(); + for(int el_id = 0; el_id < nb_; ++el_id) + { + if(!status_[el_id]) continue; + bus_status[bus_hv_id_[el_id]] = true; + bus_status[bus_lv_id_[el_id]] = true; + } + } void compute_results(const Eigen::Ref & Va, const Eigen::Ref & Vm, diff --git a/src/main.cpp b/src/main.cpp index 6545d63..b798613 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -837,6 +837,7 @@ PYBIND11_MODULE(lightsim2grid_cpp, m) // auxiliary functions .def("set_n_sub", &GridModel::set_n_sub, DocGridModel::_internal_do_not_use.c_str()) + .def("set_max_nb_bus_per_sub", &GridModel::set_max_nb_bus_per_sub, DocGridModel::_internal_do_not_use.c_str()) .def("set_load_pos_topo_vect", &GridModel::set_load_pos_topo_vect, DocGridModel::_internal_do_not_use.c_str()) .def("set_gen_pos_topo_vect", &GridModel::set_gen_pos_topo_vect, DocGridModel::_internal_do_not_use.c_str()) .def("set_line_or_pos_topo_vect", &GridModel::set_line_or_pos_topo_vect, DocGridModel::_internal_do_not_use.c_str()) diff --git a/src/powerflow_algorithm/BaseDCAlgo.tpp b/src/powerflow_algorithm/BaseDCAlgo.tpp index 0f249cb..f56c9e7 100644 --- a/src/powerflow_algorithm/BaseDCAlgo.tpp +++ b/src/powerflow_algorithm/BaseDCAlgo.tpp @@ -11,15 +11,15 @@ // TODO SLACK !!! template bool BaseDCAlgo::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 + ) { // max_iter is ignored // tol is ignored @@ -27,7 +27,6 @@ bool BaseDCAlgo::compute_pf(const Eigen::SparseMatrix & // and for the slack bus both the magnitude and the angle are used. if(!is_linear_solver_valid()) { - // err_ = ErrorType::NotInitError; return false; } BaseAlgo::reset_timer(); @@ -35,11 +34,12 @@ bool BaseDCAlgo::compute_pf(const Eigen::SparseMatrix & auto timer = CustTimer(); if(_solver_control.need_reset_solver() || _solver_control.has_dimension_changed() || - _solver_control.has_ybus_some_coeffs_zero()){ + _solver_control.has_slack_participate_changed() || // the full "ybus without slack" has changed, everything needs to be recomputed_solver_control.ybus_change_sparsity_pattern() + _solver_control.ybus_change_sparsity_pattern() || + _solver_control.has_ybus_some_coeffs_zero() + ){ reset(); } - - if (_solver_control.ybus_change_sparsity_pattern()) need_factorize_ = true; sizeYbus_with_slack_ = static_cast(Ybus.rows()); @@ -47,34 +47,32 @@ bool BaseDCAlgo::compute_pf(const Eigen::SparseMatrix & auto timer_preproc = CustTimer(); #endif // __COUT_TIMES - const CplxVect & Sbus_tmp = Sbus; - // TODO SLACK (for now i put all slacks as PV, except the first one) // this should be handled in Sbus, because we know the amount of power absorbed by the slack // so we can compute it correctly ! - // std::cout << "\t\t\tretrieve_pv_with_slack \n"; - // std::cout << "slack_ids: "; - // for(auto el: slack_ids) std::cout << el << ", "; - // std::cout << std::endl; - my_pv_ = retrieve_pv_with_slack(slack_ids, pv); - // std::cout << "my_pv_: "; - // for(auto el: my_pv_) std::cout << el << ", "; - // std::cout << std::endl; - // const Eigen::VectorXi & my_pv = pv; - - // find the slack buses - // std::cout << "\t\t\textract_slack_bus_id \n"; - slack_buses_ids_solver_ = extract_slack_bus_id(my_pv_, pq, sizeYbus_with_slack_); - sizeYbus_without_slack_ = sizeYbus_with_slack_ - slack_buses_ids_solver_.size(); - - // corresp bus -> solverbus - // std::cout << "\t\t\tfill_mat_bus_id \n"; - fill_mat_bus_id(sizeYbus_with_slack_); + if(need_factorize_ || + _solver_control.has_pv_changed()) { + my_pv_ = retrieve_pv_with_slack(slack_ids, pv); + } + + if(need_factorize_ || + _solver_control.has_pv_changed() || + _solver_control.has_pq_changed()) { + // find the slack buses + slack_buses_ids_solver_ = extract_slack_bus_id(my_pv_, pq, sizeYbus_with_slack_); + sizeYbus_without_slack_ = sizeYbus_with_slack_ - slack_buses_ids_solver_.size(); + + // corresp bus -> solverbus + fill_mat_bus_id(sizeYbus_with_slack_); + } // remove the slack bus from Ybus - // and extract only real part - // std::cout << "\t\t\tfill_dcYbus_noslack \n"; - fill_dcYbus_noslack(sizeYbus_with_slack_, Ybus); + if(need_factorize_ || + _solver_control.need_recompute_ybus() || + _solver_control.ybus_change_sparsity_pattern() || + _solver_control.has_ybus_some_coeffs_zero()) { + fill_dcYbus_noslack(sizeYbus_with_slack_, Ybus); + } #ifdef __COUT_TIMES std::cout << "\t dc: preproc: " << 1000. * timer_preproc.duration() << "ms" << std::endl; @@ -84,7 +82,14 @@ bool BaseDCAlgo::compute_pf(const Eigen::SparseMatrix & #ifdef __COUT_TIMES auto timer_solve = CustTimer(); #endif // __COUT_TIMES + bool just_factorize = false; + if(_solver_control.ybus_change_sparsity_pattern() || + _solver_control.has_ybus_some_coeffs_zero()){ + need_factorize_ = true; + } + + // initialize the solver if needed if(need_factorize_){ ErrorType status_init = _linear_solver.initialize(dcYbus_noslack_); if(status_init != ErrorType::NoError){ @@ -96,16 +101,16 @@ bool BaseDCAlgo::compute_pf(const Eigen::SparseMatrix & } // remove the slack bus from Sbus - // std::cout << "\t\t\t dcSbus_noslack_ \n"; - dcSbus_noslack_ = RealVect::Constant(sizeYbus_without_slack_, my_zero_); - for (int k=0; k < sizeYbus_with_slack_; ++k){ - if(mat_bus_id_(k) == -1) continue; // I don't add anything to the slack bus - const int col_res = mat_bus_id_(k); - dcSbus_noslack_(col_res) = std::real(Sbus_tmp(k)); + if(need_factorize_ || _solver_control.need_recompute_sbus()){ + dcSbus_noslack_ = RealVect::Constant(sizeYbus_without_slack_, my_zero_); + for (int k=0; k < sizeYbus_with_slack_; ++k){ + if(mat_bus_id_(k) == -1) continue; // I don't add anything to the slack bus + const int col_res = mat_bus_id_(k); + dcSbus_noslack_(col_res) = std::real(Sbus(k)); + } } // solve for theta: Sbus = dcY . theta (make a copy to keep dcSbus_noslack_) - // std::cout << "\t\t\t Va_dc_without_slack \n"; RealVect Va_dc_without_slack = dcSbus_noslack_; ErrorType error = _linear_solver.solve(dcYbus_noslack_, Va_dc_without_slack, just_factorize); if(error != ErrorType::NoError){