diff --git a/CHANGELOG.rst b/CHANGELOG.rst index 0861eed..7d8bacd 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -21,8 +21,12 @@ Change Log [0.8.1] 2024-xx-yy -------------------- - [FIXED] a bug with shunts when `nb_busbar_per_sub` >= 2 +- [FIXED] some bugs preventing backward compatibility - [IMPROVED] time measurments in python and c++ - [IMPROVED] now test lightsim2grid with oldest grid2op version +- [IMPROVED] speed, by accelerating the reading back of the data (now read only once and then + pointers are re used) +- [IMPROVED] c++ side avoid allocating memory (which allow to gain speed python side too) [0.8.0] 2024-03-18 -------------------- diff --git a/lightsim2grid/lightSimBackend.py b/lightsim2grid/lightSimBackend.py index 7458b7c..81b2ad4 100644 --- a/lightsim2grid/lightSimBackend.py +++ b/lightsim2grid/lightSimBackend.py @@ -141,6 +141,7 @@ def __init__(self, self._timer_postproc = 0. self._timer_solver = 0. self._timer_read_data_back = 0. + self._timer_test_read = 0. self.next_prod_p = None # this vector is updated with the action that will modify the environment # it is done to keep track of the redispatching @@ -221,6 +222,7 @@ def __init__(self, self._timer_preproc = 0. self._timer_solver = 0. self._timer_read_data_back = 0. + self._timer_test_read = 0. # hack for the storage unit: # in grid2op, for simplicity, I suppose that if a storage is alone on a busbar, and @@ -233,6 +235,16 @@ def __init__(self, # backend SHOULD not do these kind of stuff self._idx_hack_storage = [] + # speed optimization + self._lineor_res = None + self._lineex_res = None + self._load_res = None + self._gen_res = None + self._shunt_res = None + self._trafo_hv_res = None + self._trafo_lv_res = None + self._storage_res = None + def _aux_init_super(self, detailed_infos_for_cascading_failures, can_be_copied, @@ -267,7 +279,11 @@ def _aux_init_super(self, "you cannot set max_iter, tol nor solver_type arguments.") Backend.__init__(self, detailed_infos_for_cascading_failures=detailed_infos_for_cascading_failures) - + + if hasattr(type(self), "can_handle_more_than_2_busbar"): + # do not forget to propagate this if needed + self.can_handle_more_than_2_busbar() + def turnedoff_no_pv(self): self._turned_off_pv = False self._grid.turnedoff_no_pv() @@ -480,6 +496,7 @@ def load_grid(self, path=None, filename=None): else: raise BackendError(f"Impossible to initialize the backend with '{self._loader_method}'") self._grid.tell_solver_need_reset() + self._reset_res_pointers() # force the re reading of the accessors at next powerflow def _should_not_have_to_do_this(self, path=None, filename=None): # included in grid2op now ! @@ -871,6 +888,7 @@ def _aux_finish_setup_after_reading(self): self.__me_at_init = self._grid.copy() self.__init_topo_vect = np.ones(cls.dim_topo, dtype=dt_int) self.__init_topo_vect[:] = self.topo_vect + self.sh_bus[:] = 1 def assert_grid_correct_after_powerflow(self): """ @@ -991,6 +1009,8 @@ def apply_action(self, backendAction): self._grid.change_bus_shunt(sh_id, type(self).local_bus_to_global_int(new_bus, self.shunt_to_subid[sh_id])) else: self._grid.change_bus_shunt(sh_id, self.shunt_to_subid[sh_id] + (new_bus == 2) * type(self).n_sub) + # remember the topology not to need to read it back from the grid + self.sh_bus[shunt_bus.changed] = shunt_bus.values[shunt_bus.changed] for sh_id, new_p in shunt_p: self._grid.change_p_shunt(sh_id, new_p) for sh_id, new_q in shunt_q: @@ -1002,6 +1022,26 @@ def _handle_dist_slack(self): if self._dist_slack_non_renew: self._grid.update_slack_weights(type(self).gen_redispatchable) + def _fetch_grid_data(self): + beg_test = time.perf_counter() + if self._lineor_res is None: + self._lineor_res = self._grid.get_lineor_res() + if self._lineex_res is None: + self._lineex_res = self._grid.get_lineex_res() + if self._load_res is None: + self._load_res = self._grid.get_loads_res() + if self._gen_res is None: + self._gen_res = self._grid.get_gen_res() + if self._trafo_hv_res is None: + self._trafo_hv_res = self._grid.get_trafohv_res() + if self._trafo_lv_res is None: + self._trafo_lv_res = self._grid.get_trafolv_res() + if self._storage_res is None: + self._storage_res = self._grid.get_storages_res() + if self._shunt_res is None: + self._shunt_res = self._grid.get_shunts_res() + self._timer_test_read += time.perf_counter() - beg_test + def runpf(self, is_dc=False): my_exc_ = None res = False @@ -1026,9 +1066,7 @@ 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(f"{self.V[:14] = }") Vdc = self._grid.dc_pf(copy.deepcopy(self.V), self.max_it, self.tol) - # 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()}") @@ -1060,22 +1098,24 @@ def runpf(self, is_dc=False): beg_readback = time.perf_counter() self.V[:] = V + self._fetch_grid_data() + (self.p_or[:self.__nb_powerline], self.q_or[:self.__nb_powerline], self.v_or[:self.__nb_powerline], - self.a_or[:self.__nb_powerline]) = self._grid.get_lineor_res() + self.a_or[:self.__nb_powerline]) = self._lineor_res (self.p_or[self.__nb_powerline:], self.q_or[self.__nb_powerline:], self.v_or[self.__nb_powerline:], - self.a_or[self.__nb_powerline:]) = self._grid.get_trafohv_res() + self.a_or[self.__nb_powerline:]) = self._trafo_hv_res (self.p_ex[:self.__nb_powerline], self.q_ex[:self.__nb_powerline], self.v_ex[:self.__nb_powerline], - self.a_ex[:self.__nb_powerline]) = self._grid.get_lineex_res() + self.a_ex[:self.__nb_powerline]) = self._lineex_res (self.p_ex[self.__nb_powerline:], self.q_ex[self.__nb_powerline:], self.v_ex[self.__nb_powerline:], - self.a_ex[self.__nb_powerline:]) = self._grid.get_trafolv_res() + self.a_ex[self.__nb_powerline:]) = self._trafo_lv_res self.a_or *= 1000. # kA in lightsim, A expected in grid2op self.a_ex *= 1000. # kA in lightsim, A expected in grid2op @@ -1085,10 +1125,10 @@ def runpf(self, is_dc=False): self.a_ex[~np.isfinite(self.a_ex)] = 0. self.v_ex[~np.isfinite(self.v_ex)] = 0. - self.load_p[:], self.load_q[:], self.load_v[:] = self._grid.get_loads_res() - self.prod_p[:], self.prod_q[:], self.prod_v[:] = self._grid.get_gen_res() + self.load_p[:], self.load_q[:], self.load_v[:] = self._load_res + self.prod_p[:], self.prod_q[:], self.prod_v[:] = self._gen_res if self.__has_storage: - self.storage_p[:], self.storage_q[:], self.storage_v[:] = self._grid.get_storages_res() + self.storage_p[:], self.storage_q[:], self.storage_v[:] = self._storage_res self.storage_v[self.storage_v == -1.] = 0. # voltage is 0. for disconnected elements in grid2op self._timer_read_data_back += time.perf_counter() - beg_readback @@ -1168,7 +1208,18 @@ def _fill_nans(self): self.storage_v[:] = np.NaN self.storage_theta[:] = np.NaN self.V[:] = self._grid.get_init_vm_pu() # reset the V to its "original" value (see issue 30) - + self._reset_res_pointers() + + def _reset_res_pointers(self): + self._lineor_res = None + self._lineex_res = None + self._load_res = None + self._gen_res = None + self._trafo_hv_res = None + self._trafo_lv_res = None + self._storage_res = None + self._shunt_res = None + def __deepcopy__(self, memo): result = self.copy() memo[id(self)] = result @@ -1228,7 +1279,8 @@ def copy(self): "supported_grid_format", "max_it", "tol", "_turned_off_pv", "_dist_slack_non_renew", "_use_static_gen", "_loader_method", "_loader_kwargs", - "_stop_if_load_disco", "_stop_if_gen_disco" + "_stop_if_load_disco", "_stop_if_gen_disco", + "_timer_test_read" ] for attr_nm in li_regular_attr: if hasattr(self, attr_nm): @@ -1288,7 +1340,10 @@ def copy(self): res._backend_action_class = self._backend_action_class # this is const res.__init_topo_vect = self.__init_topo_vect res.available_solvers = self.available_solvers - + res._reset_res_pointers() + res._fetch_grid_data() + + # assign back "self" attributes self._grid = mygrid self.init_pp_backend = inippbackend self.__me_at_init = __me_at_init @@ -1336,10 +1391,7 @@ def storages_info(self): def shunt_info(self): return self.cst_1 * self.sh_p, self.cst_1 * self.sh_q, self.cst_1 * self.sh_v, self.sh_bus - def _set_shunt_info(self): - tick = time.perf_counter() - self.sh_p[:], self.sh_q[:], self.sh_v[:] = self._grid.get_shunts_res() - shunt_bus = np.array([self._grid.get_bus_shunt(i) for i in range(type(self).n_shunt)], dtype=dt_int) + def _compute_shunt_bus_with_compat(self, shunt_bus): cls = type(self) if hasattr(cls, "global_bus_to_local"): self.sh_bus[:] = cls.global_bus_to_local(shunt_bus, cls.shunt_to_subid) @@ -1353,8 +1405,30 @@ def _set_shunt_info(self): for i in range(n_busbar_per_sub): res[(i * cls.n_sub <= shunt_bus) & (shunt_bus < (i+1) * cls.n_sub)] = i + 1 res[shunt_bus == -1] = -1 + self.sh_bus[:] = res + + def _set_shunt_info(self): + tick = time.perf_counter() + self.sh_p[:], self.sh_q[:], self.sh_v[:] = self._shunt_res + # cls = type(self) + # shunt_bus = np.array([self._grid.get_bus_shunt(i) for i in range(cls.n_shunt)], dtype=dt_int) + # shunt_bus = self._grid.get_all_shunt_buses() + # if hasattr(cls, "global_bus_to_local"): + # self.sh_bus[:] = cls.global_bus_to_local(shunt_bus, cls.shunt_to_subid) + # else: + # res = (1 * shunt_bus).astype(dt_int) # make a copy + # if hasattr(cls, "n_busbar_per_sub"): + # n_busbar_per_sub = cls.n_busbar_per_sub + # else: + # # backward compat when this was not defined: + # n_busbar_per_sub = DEFAULT_N_BUSBAR_PER_SUB + # for i in range(n_busbar_per_sub): + # res[(i * cls.n_sub <= shunt_bus) & (shunt_bus < (i+1) * cls.n_sub)] = i + 1 + # res[shunt_bus == -1] = -1 + # self.sh_bus[:] = res self.sh_v[self.sh_v == -1.] = 0. # in grid2op disco element have voltage of 0. and -1. self._timer_read_data_back += time.perf_counter() - tick + # self._timer_test_read += time.perf_counter() - tick def _disconnect_line(self, id_): self.topo_vect[self.line_ex_pos_topo_vect[id_]] = -1 @@ -1373,11 +1447,13 @@ def reset(self, grid_path, grid_filename=None): self._grid.unset_changes() self._grid.change_solver(self.__current_solver_type) self._handle_turnedoff_pv() - self.topo_vect[:] = self.__init_topo_vect self.comp_time = 0. self.timer_gridmodel_xx_pf = 0. self._timer_postproc = 0. self._timer_preproc = 0. self._timer_solver = 0. self._timer_read_data_back = 0. + self._timer_test_read = 0. self._grid.tell_solver_need_reset() + self.sh_bus[:] = 1 # TODO self._compute_shunt_bus_with_compat(self._grid.get_all_shunt_buses()) + self.topo_vect[:] = self.__init_topo_vect # TODO diff --git a/src/GridModel.cpp b/src/GridModel.cpp index 67b4853..068e235 100644 --- a/src/GridModel.cpp +++ b/src/GridModel.cpp @@ -371,7 +371,7 @@ CplxVect GridModel::ac_pf(const CplxVect & Vinit, bool conv = false; CplxVect res = CplxVect(); - reset_results(); // reset the results + // reset_results(); // clear the results No need to do it, results are neceassirly set or reset in post process // pre process the data to define a proper jacobian matrix, the proper voltage vector etc. bool is_ac = true; @@ -901,7 +901,7 @@ CplxVect GridModel::dc_pf(const CplxVect & Vinit, bool conv = false; CplxVect res = CplxVect(); - reset_results(); // reset the results + // reset_results(); // clear the results No need to do it, results are neceassirly set or reset in post process // pre process the data to define a proper jacobian matrix, the proper voltage vector etc. bool is_ac = false; diff --git a/src/GridModel.h b/src/GridModel.h index 7a4b69c..f9a863b 100644 --- a/src/GridModel.h +++ b/src/GridModel.h @@ -494,6 +494,8 @@ class GridModel : public GenericContainer Eigen::Ref get_dclineor_theta() const {return dc_lines_.get_theta_or();} Eigen::Ref get_dclineex_theta() const {return dc_lines_.get_theta_ex();} + Eigen::Ref get_all_shunt_buses() const {return shunts_.get_buses();} + // get some internal information, be cerafull the ID of the buses might not be the same // TODO convert it back to this ID, that will make copies, but who really cares ? Eigen::SparseMatrix get_Ybus(){ diff --git a/src/element_container/GeneratorContainer.cpp b/src/element_container/GeneratorContainer.cpp index 7598301..973c6ba 100644 --- a/src/element_container/GeneratorContainer.cpp +++ b/src/element_container/GeneratorContainer.cpp @@ -54,6 +54,7 @@ void GeneratorContainer::init(const RealVect & generators_p, turnedoff_gen_pv_ = true; voltage_regulator_on_ = std::vector(generators_p.size(), true); q_mvar_ = RealVect::Zero(generators_p.size()); + reset_results(); } void GeneratorContainer::init_full(const RealVect & generators_p, @@ -222,10 +223,10 @@ void GeneratorContainer::compute_results(const Eigen::Ref & Va, } void GeneratorContainer::reset_results(){ - res_p_ = RealVect(); // in MW - res_q_ = RealVect(); // in MVar - res_v_ = RealVect(); // in kV - res_theta_ = RealVect(); // in deg + res_p_ = RealVect(nb()); // in MW + res_q_ = RealVect(nb()); // in MVar + res_v_ = RealVect(nb()); // in kV + res_theta_ = RealVect(nb()); // in deg // bus_slack_weight_ = RealVect(); } @@ -414,13 +415,23 @@ void GeneratorContainer::set_q(const RealVect & reactive_mismatch, const RealVect & total_q_max_per_bus) { const int nb_gen = nb(); - res_q_ = RealVect::Constant(nb_gen, 0.); - if(!ac) return; // do not consider Q values in dc mode + // res_q_ = RealVect::Constant(nb_gen, 0.); + if(!ac) + { + // do not consider Q values in dc mode + for(int gen_id = 0; gen_id < nb_gen; ++gen_id) res_q_(gen_id) = 0.; + return; + } + real_type eps_q = 1e-8; for(int gen_id = 0; gen_id < nb_gen; ++gen_id) { real_type real_q = 0.; - if(!status_[gen_id]) continue; // set at 0 for disconnected generators + if(!status_[gen_id]){ + // set at 0 for disconnected generators + res_q_(gen_id) = 0.; + continue; + } if (!voltage_regulator_on_[gen_id]) continue; // gen is purposedly not pv if ((!turnedoff_gen_pv_) && p_mw_(gen_id) == 0.) continue; // in this case turned off generators are not pv diff --git a/src/element_container/GenericContainer.cpp b/src/element_container/GenericContainer.cpp index cf6dc51..fc8ce14 100644 --- a/src/element_container/GenericContainer.cpp +++ b/src/element_container/GenericContainer.cpp @@ -14,7 +14,7 @@ const int GenericContainer::_deactivated_bus_id = -1; // TODO all functions bellow are generic ! Make a base class for that -void GenericContainer::_get_amps(RealVect & a, const RealVect & p, const RealVect & q, const RealVect & v){ +void GenericContainer::_get_amps(RealVect & a, const RealVect & p, const RealVect & q, const RealVect & v) const { const real_type _1_sqrt_3 = 1.0 / std::sqrt(3.); RealVect p2q2 = p.array() * p.array() + q.array() * q.array(); p2q2 = p2q2.array().cwiseSqrt(); @@ -96,7 +96,7 @@ void GenericContainer::_change_bus(int el_id, int new_bus_me_id, Eigen::VectorXi bus_me_id = new_bus_me_id; } -int GenericContainer::_get_bus(int el_id, const std::vector & status_, const Eigen::VectorXi & bus_id_) +int GenericContainer::_get_bus(int el_id, const std::vector & status_, const Eigen::VectorXi & bus_id_) const { int res; bool val = status_.at(el_id); // also check if the el_id is out of bound @@ -114,12 +114,14 @@ void GenericContainer::v_kv_from_vpu(const Eigen::Ref & Va, const Eigen::VectorXi & bus_me_id, const std::vector & id_grid_to_solver, const RealVect & bus_vn_kv, - RealVect & v) + RealVect & v) const { - v = RealVect::Constant(nb_element, -1.0); for(int el_id = 0; el_id < nb_element; ++el_id){ // if the element is disconnected, i leave it like that - if(!status[el_id]) continue; + if(!status[el_id]) { + v(el_id) = -1; + continue; + } int el_bus_me_id = bus_me_id(el_id); int bus_solver_id = id_grid_to_solver[el_bus_me_id]; if(bus_solver_id == _deactivated_bus_id){ @@ -135,7 +137,6 @@ void GenericContainer::v_kv_from_vpu(const Eigen::Ref & Va, } } - void GenericContainer::v_deg_from_va(const Eigen::Ref & Va, const Eigen::Ref & Vm, const std::vector & status, @@ -143,11 +144,14 @@ void GenericContainer::v_deg_from_va(const Eigen::Ref & Va, const Eigen::VectorXi & bus_me_id, const std::vector & id_grid_to_solver, const RealVect & bus_vn_kv, - RealVect & theta){ - theta = RealVect::Constant(nb_element, 0.0); + RealVect & theta) const +{ for(int el_id = 0; el_id < nb_element; ++el_id){ // if the element is disconnected, i leave it like that - if(!status[el_id]) continue; + if(!status[el_id]) { + theta(el_id) = -1; + continue; + } int el_bus_me_id = bus_me_id(el_id); int bus_solver_id = id_grid_to_solver[el_bus_me_id]; if(bus_solver_id == _deactivated_bus_id){ diff --git a/src/element_container/GenericContainer.h b/src/element_container/GenericContainer.h index b1fa679..98050a0 100644 --- a/src/element_container/GenericContainer.h +++ b/src/element_container/GenericContainer.h @@ -83,7 +83,9 @@ class GenericContainer : public BaseConstants int nb_line, bool transpose) const {}; - virtual void fillYbus(Eigen::SparseMatrix & res, bool ac, const std::vector & id_grid_to_solver) {}; + // no more used ! + virtual void fillYbus(Eigen::SparseMatrix & res, bool ac, const std::vector & id_grid_to_solver) const {}; + 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, @@ -120,12 +122,12 @@ class GenericContainer : public BaseConstants void _reactivate(int el_id, std::vector & status); void _deactivate(int el_id, std::vector & status); void _change_bus(int el_id, int new_bus_me_id, Eigen::VectorXi & el_bus_ids, SolverControl & solver_control, int nb_bus); - int _get_bus(int el_id, const std::vector & status_, const Eigen::VectorXi & bus_id_); + int _get_bus(int el_id, const std::vector & status_, const Eigen::VectorXi & bus_id_) const; /** compute the amps from the p, the q and the v (v should NOT be pair unit) **/ - void _get_amps(RealVect & a, const RealVect & p, const RealVect & q, const RealVect & v); + void _get_amps(RealVect & a, const RealVect & p, const RealVect & q, const RealVect & v) const; /** convert v from pu to v in kv (and assign it to the right element...) @@ -137,7 +139,7 @@ class GenericContainer : public BaseConstants const Eigen::VectorXi & bus_me_id, const std::vector & id_grid_to_solver, const RealVect & bus_vn_kv, - RealVect & v); + RealVect & v) const; /** @@ -150,13 +152,13 @@ class GenericContainer : public BaseConstants const Eigen::VectorXi & bus_me_id, const std::vector & id_grid_to_solver, const RealVect & bus_vn_kv, - RealVect & v); + RealVect & v) const; /** check the size of the elements **/ template - void check_size(const T & container, intType size, const std::string & container_name) + void check_size(const T & container, intType size, const std::string & container_name) const { if(static_cast(container.size()) != size) throw std::runtime_error(container_name + " do not have the proper size"); } diff --git a/src/element_container/LineContainer.cpp b/src/element_container/LineContainer.cpp index efeac7d..662a1b1 100644 --- a/src/element_container/LineContainer.cpp +++ b/src/element_container/LineContainer.cpp @@ -44,6 +44,7 @@ void LineContainer::init(const RealVect & branch_r, powerlines_x_ = branch_x; status_ = std::vector(branch_r.size(), true); // by default everything is connected _update_model_coeffs(); + reset_results(); } void LineContainer::init(const RealVect & branch_r, @@ -82,6 +83,7 @@ void LineContainer::init(const RealVect & branch_r, powerlines_x_ = branch_x; status_ = std::vector(branch_r.size(), true); // by default everything is connected _update_model_coeffs(); + reset_results(); } LineContainer::StateRes LineContainer::get_state() const @@ -98,7 +100,6 @@ LineContainer::StateRes LineContainer::get_state() const } void LineContainer::set_state(LineContainer::StateRes & my_state) { - reset_results(); names_ = std::get<0>(my_state); std::vector & branch_r = std::get<1>(my_state); std::vector & branch_x = std::get<2>(my_state); @@ -121,6 +122,7 @@ void LineContainer::set_state(LineContainer::StateRes & my_state) status_ = status; _update_model_coeffs(); + reset_results(); } void LineContainer::_update_model_coeffs() @@ -355,14 +357,16 @@ void LineContainer::fillBf_for_PTDF(std::vector > & Bf void LineContainer::reset_results() { - res_powerline_por_ = RealVect(); // in MW - res_powerline_qor_ = RealVect(); // in MVar - res_powerline_vor_ = RealVect(); // in kV - res_powerline_aor_ = RealVect(); // in kA - res_powerline_pex_ = RealVect(); // in MW - res_powerline_qex_ = RealVect(); // in MVar - res_powerline_vex_ = RealVect(); // in kV - res_powerline_aex_ = RealVect(); // in kA + res_powerline_por_ = RealVect(nb()); // in MW + res_powerline_qor_ = RealVect(nb()); // in MVar + res_powerline_vor_ = RealVect(nb()); // in kV + res_powerline_aor_ = RealVect(nb()); // in kA + res_powerline_pex_ = RealVect(nb()); // in MW + res_powerline_qex_ = RealVect(nb()); // in MVar + res_powerline_vex_ = RealVect(nb()); // in kV + res_powerline_aex_ = RealVect(nb()); // in kA + res_powerline_thetaor_ = RealVect(nb()); + res_powerline_thetaex_ = RealVect(nb()); } void LineContainer::compute_results(const Eigen::Ref & Va, @@ -375,19 +379,21 @@ void LineContainer::compute_results(const Eigen::Ref & Va, { // it needs to be initialized at 0. Eigen::Index nb_element = nb(); - res_powerline_por_ = RealVect::Constant(nb_element, my_zero_); // in MW - res_powerline_qor_ = RealVect::Constant(nb_element, my_zero_); // in MVar - res_powerline_vor_ = RealVect::Constant(nb_element, my_zero_); // in kV - res_powerline_aor_ = RealVect::Constant(nb_element, my_zero_); // in kA - res_powerline_pex_ = RealVect::Constant(nb_element, my_zero_); // in MW - res_powerline_qex_ = RealVect::Constant(nb_element, my_zero_); // in MVar - res_powerline_vex_ = RealVect::Constant(nb_element, my_zero_); // in kV - res_powerline_aex_ = RealVect::Constant(nb_element, my_zero_); // in kA - res_powerline_thetaor_ = RealVect::Constant(nb_element, my_zero_); // in kV - res_powerline_thetaex_ = RealVect::Constant(nb_element, my_zero_); // in kV for(Eigen::Index line_id = 0; line_id < nb_element; ++line_id){ // don't do anything if the element is disconnected - if(!status_[line_id]) continue; + if(!status_[line_id]) { + res_powerline_por_(line_id) = 0.0; // in MW + res_powerline_qor_(line_id) = 0.0; // in MVar + res_powerline_vor_(line_id) = 0.0; // in kV + res_powerline_aor_(line_id) = 0.0; // in kA + res_powerline_pex_(line_id) = 0.0; // in MW + res_powerline_qex_(line_id) = 0.0; // in MVar + res_powerline_vex_(line_id) = 0.0; // in kV + res_powerline_aex_(line_id) = 0.0; // in kA + res_powerline_thetaor_(line_id) = 0.0; // in kV + res_powerline_thetaex_(line_id) = 0.0; // in kV + continue; + } // connectivity int bus_or_id_me = bus_or_id_(line_id); @@ -447,6 +453,8 @@ void LineContainer::compute_results(const Eigen::Ref & Va, res_powerline_por_(line_id) = (std::real(ydc_ff_(line_id)) * Va(bus_or_solver_id) + std::real(ydc_ft_(line_id)) * Va(bus_ex_solver_id)) * sn_mva; res_powerline_pex_(line_id) = (std::real(ydc_tt_(line_id)) * Va(bus_ex_solver_id) + std::real(ydc_tf_(line_id)) * Va(bus_or_solver_id)) * sn_mva; + res_powerline_qor_(line_id) = 0.; + res_powerline_qex_(line_id) = 0.; // for the voltage (by hypothesis vm = 1) // res_powerline_vor_(line_id) = bus_vn_kv_or; // res_powerline_vex_(line_id) = bus_vn_kv_ex; diff --git a/src/element_container/LoadContainer.cpp b/src/element_container/LoadContainer.cpp index d9d0a1b..afd0718 100644 --- a/src/element_container/LoadContainer.cpp +++ b/src/element_container/LoadContainer.cpp @@ -8,6 +8,7 @@ #include "LoadContainer.h" #include +#include void LoadContainer::init(const RealVect & loads_p, const RealVect & loads_q, @@ -22,6 +23,7 @@ void LoadContainer::init(const RealVect & loads_p, q_mvar_ = loads_q; bus_id_ = loads_bus_id; status_ = std::vector(loads_p.size(), true); + reset_results(); } @@ -37,7 +39,6 @@ LoadContainer::StateRes LoadContainer::get_state() const void LoadContainer::set_state(LoadContainer::StateRes & my_state ) { - reset_results(); names_ = std::get<0>(my_state); std::vector & p_mw = std::get<1>(my_state); std::vector & q_mvar = std::get<2>(my_state); @@ -50,6 +51,7 @@ void LoadContainer::set_state(LoadContainer::StateRes & my_state ) q_mvar_ = RealVect::Map(&q_mvar[0], q_mvar.size()); bus_id_ = Eigen::VectorXi::Map(&bus_id[0], bus_id.size()); status_ = status; + reset_results(); } @@ -87,17 +89,23 @@ void LoadContainer::compute_results(const Eigen::Ref & Va, real_type sn_mva, bool ac) { - int nb_load = nb(); + const int nb_load = nb(); v_kv_from_vpu(Va, Vm, status_, nb_load, bus_id_, id_grid_to_solver, bus_vn_kv, res_v_); v_deg_from_va(Va, Vm, status_, nb_load, bus_id_, id_grid_to_solver, bus_vn_kv, res_theta_); res_p_ = p_mw_; - res_q_ = q_mvar_; + if(ac) res_q_ = q_mvar_; + else{ + // no q in DC mode + for(int load_id = 0; load_id < nb_load; ++load_id) res_q_(load_id) = 0.; + } } void LoadContainer::reset_results(){ - res_p_ = RealVect(); // in MW - res_q_ = RealVect(); // in MVar - res_v_ = RealVect(); // in kV + // std::cout << "Loads reset_results \n"; + res_p_ = RealVect(nb()); // in MW + res_q_ = RealVect(nb()); // in MVar + res_v_ = RealVect(nb()); // in kV + res_theta_ = RealVect(nb()); // in deg } void LoadContainer::change_p(int load_id, real_type new_p, SolverControl & solver_control) diff --git a/src/element_container/SGenContainer.cpp b/src/element_container/SGenContainer.cpp index 6300727..8690b74 100644 --- a/src/element_container/SGenContainer.cpp +++ b/src/element_container/SGenContainer.cpp @@ -34,6 +34,7 @@ void SGenContainer::init(const RealVect & sgen_p, q_max_mvar_ = sgen_qmax; bus_id_ = sgen_bus_id; status_ = std::vector(sgen_p.size(), true); + reset_results(); } SGenContainer::StateRes SGenContainer::get_state() const @@ -51,9 +52,7 @@ SGenContainer::StateRes SGenContainer::get_state() const } void SGenContainer::set_state(SGenContainer::StateRes & my_state ) -{ - reset_results(); - +{ names_ = std::get<0>(my_state); std::vector & p_mw = std::get<1>(my_state); std::vector & q_mvar = std::get<2>(my_state); @@ -82,6 +81,7 @@ void SGenContainer::set_state(SGenContainer::StateRes & my_state ) q_max_mvar_ = RealVect::Map(&q_max[0], size); bus_id_ = Eigen::VectorXi::Map(&bus_id[0], bus_id.size()); status_ = status; + reset_results(); } void SGenContainer::fillSbus(CplxVect & Sbus, const std::vector & id_grid_to_solver, bool ac) const { @@ -119,13 +119,17 @@ void SGenContainer::compute_results(const Eigen::Ref & Va, v_deg_from_va(Va, Vm, status_, nb_sgen, bus_id_, id_grid_to_solver, bus_vn_kv, res_theta_); res_p_ = p_mw_; if(ac) res_q_ = q_mvar_; - else res_q_ = RealVect::Zero(nb_sgen); + else{ + // no q in DC mode + for(int sgen_id = 0; sgen_id < nb_sgen; ++sgen_id) res_q_(sgen_id) = 0.; + } } void SGenContainer::reset_results(){ - res_p_ = RealVect(); // in MW - res_q_ = RealVect(); // in MVar - res_v_ = RealVect(); // in kV + res_p_ = RealVect(nb()); // in MW + res_q_ = RealVect(nb()); // in MVar + res_v_ = RealVect(nb()); // in kV + res_theta_ = RealVect(nb()); // in deg } void SGenContainer::change_p(int sgen_id, real_type new_p, SolverControl & solver_control) diff --git a/src/element_container/ShuntContainer.cpp b/src/element_container/ShuntContainer.cpp index 4b5687b..4bd38a5 100644 --- a/src/element_container/ShuntContainer.cpp +++ b/src/element_container/ShuntContainer.cpp @@ -23,6 +23,7 @@ void ShuntContainer::init(const RealVect & shunt_p_mw, q_mvar_ = shunt_q_mvar; bus_id_ = shunt_bus_id; status_ = std::vector(p_mw_.size(), true); // by default everything is connected + reset_results(); } ShuntContainer::StateRes ShuntContainer::get_state() const @@ -37,7 +38,6 @@ ShuntContainer::StateRes ShuntContainer::get_state() const void ShuntContainer::set_state(ShuntContainer::StateRes & my_state ) { - reset_results(); names_ = std::get<0>(my_state); std::vector & p_mw = std::get<1>(my_state); std::vector & q_mvar = std::get<2>(my_state); @@ -50,6 +50,7 @@ void ShuntContainer::set_state(ShuntContainer::StateRes & my_state ) q_mvar_ = RealVect::Map(&q_mvar[0], q_mvar.size()); bus_id_ = Eigen::VectorXi::Map(&bus_id[0], bus_id.size()); status_ = status; + reset_results(); } void ShuntContainer::fillYbus(std::vector > & res, @@ -147,10 +148,14 @@ void ShuntContainer::compute_results(const Eigen::Ref & Va, const int nb_shunt = static_cast(p_mw_.size()); v_kv_from_vpu(Va, Vm, status_, nb_shunt, bus_id_, id_grid_to_solver, bus_vn_kv, res_v_); v_deg_from_va(Va, Vm, status_, nb_shunt, bus_id_, id_grid_to_solver, bus_vn_kv, res_theta_); - res_p_ = RealVect::Constant(nb_shunt, my_zero_); - res_q_ = RealVect::Constant(nb_shunt, my_zero_); + // res_p_ = RealVect::Constant(nb_shunt, my_zero_); + // res_q_ = RealVect::Constant(nb_shunt, my_zero_); for(int shunt_id = 0; shunt_id < nb_shunt; ++shunt_id){ - if(!status_[shunt_id]) continue; + if(!status_[shunt_id]) { + res_p_(shunt_id) = my_zero_; + res_q_(shunt_id) = my_zero_; + continue; + } int bus_id_me = bus_id_(shunt_id); int bus_solver_id = id_grid_to_solver[bus_id_me]; if(bus_solver_id == _deactivated_bus_id){ @@ -163,13 +168,15 @@ void ShuntContainer::compute_results(const Eigen::Ref & Va, cplx_type s = E * I; res_p_(shunt_id) = std::real(s) * sn_mva; if(ac) res_q_(shunt_id) = std::imag(s) * sn_mva; + else res_q_(shunt_id) = my_zero_; } } void ShuntContainer::reset_results(){ - res_p_ = RealVect(); // in MW - res_q_ = RealVect(); // in MVar - res_v_ = RealVect(); // in kV + res_p_ = RealVect(nb()); // in MW + res_q_ = RealVect(nb()); // in MVar + res_v_ = RealVect(nb()); // in kV + res_theta_ = RealVect(nb()); // in deg } void ShuntContainer::change_p(int shunt_id, real_type new_p, SolverControl & solver_control) diff --git a/src/element_container/ShuntContainer.h b/src/element_container/ShuntContainer.h index 8458acf..40cf7cb 100644 --- a/src/element_container/ShuntContainer.h +++ b/src/element_container/ShuntContainer.h @@ -146,7 +146,9 @@ class ShuntContainer : public GenericContainer void change_bus(int shunt_id, int new_bus_id, SolverControl & solver_control, int nb_bus) {_change_bus(shunt_id, new_bus_id, bus_id_, solver_control, nb_bus);} void change_p(int shunt_id, real_type new_p, SolverControl & solver_control); void change_q(int shunt_id, real_type new_q, SolverControl & solver_control); - int get_bus(int shunt_id) {return _get_bus(shunt_id, status_, bus_id_);} + int get_bus(int shunt_id) const {return _get_bus(shunt_id, status_, bus_id_);} + Eigen::Ref get_buses() const {return 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); diff --git a/src/element_container/TrafoContainer.cpp b/src/element_container/TrafoContainer.cpp index f33e1cd..b38ce8a 100644 --- a/src/element_container/TrafoContainer.cpp +++ b/src/element_container/TrafoContainer.cpp @@ -52,6 +52,7 @@ void TrafoContainer::init(const RealVect & trafo_r, is_tap_hv_side_ = trafo_tap_hv; status_ = std::vector(trafo_r.size(), true); _update_model_coeffs(); + reset_results(); } @@ -74,8 +75,6 @@ TrafoContainer::StateRes TrafoContainer::get_state() const void TrafoContainer::set_state(TrafoContainer::StateRes & my_state) { - reset_results(); - names_ = std::get<0>(my_state); std::vector & branch_r = std::get<1>(my_state); std::vector & branch_x = std::get<2>(my_state); @@ -111,6 +110,7 @@ void TrafoContainer::set_state(TrafoContainer::StateRes & my_state) shift_ = RealVect::Map(&shift[0], size); is_tap_hv_side_ = is_tap_hv_side; _update_model_coeffs(); + reset_results(); } @@ -271,19 +271,21 @@ void TrafoContainer::compute_results(const Eigen::Ref & Va, { // it needs to be initialized at 0. const int nb_element = nb(); - res_p_hv_ = RealVect::Constant(nb_element, 0.0); // in MW - res_q_hv_ = RealVect::Constant(nb_element, 0.0); // in MVar - res_v_hv_ = RealVect::Constant(nb_element, 0.0); // in kV - res_a_hv_ = RealVect::Constant(nb_element, 0.0); // in kA - res_p_lv_ = RealVect::Constant(nb_element, 0.0); // in MW - res_q_lv_ = RealVect::Constant(nb_element, 0.0); // in MVar - res_v_lv_ = RealVect::Constant(nb_element, 0.0); // in kV - res_a_lv_ = RealVect::Constant(nb_element, 0.0); // in kA - res_theta_hv_ = RealVect::Constant(nb_element, 0.0); // in degree - res_theta_lv_ = RealVect::Constant(nb_element, 0.0); // in degree for(int trafo_id = 0; trafo_id < nb_element; ++trafo_id){ // don't do anything if the element is disconnected - if(!status_[trafo_id]) continue; + if(!status_[trafo_id]) { + res_p_hv_(trafo_id) = 0.0; // in MW + res_q_hv_(trafo_id) = 0.0; // in MVar + res_v_hv_(trafo_id) = 0.0; // in kV + res_a_hv_(trafo_id) = 0.0; // in kA + res_p_lv_(trafo_id) = 0.0; // in MW + res_q_lv_(trafo_id) = 0.0; // in MVar + res_v_lv_(trafo_id) = 0.0; // in kV + res_a_lv_(trafo_id) = 0.0; // in kA + res_theta_hv_(trafo_id) = 0.0; // in degree + res_theta_lv_(trafo_id) = 0.0; // in degree + continue; + } // connectivity int bus_hv_id_me = bus_hv_id_(trafo_id); @@ -341,7 +343,10 @@ void TrafoContainer::compute_results(const Eigen::Ref & Va, // result of the dc powerflow res_p_hv_(trafo_id) = (std::real(ydc_ff_(trafo_id)) * Va(bus_hv_solver_id) + std::real(ydc_ft_(trafo_id)) * Va(bus_lv_solver_id) - dc_x_tau_shift_(trafo_id) ) * sn_mva; res_p_lv_(trafo_id) = (std::real(ydc_tt_(trafo_id)) * Va(bus_lv_solver_id) + std::real(ydc_tf_(trafo_id)) * Va(bus_hv_solver_id) + dc_x_tau_shift_(trafo_id) ) * sn_mva; - + + res_q_hv_(trafo_id) = 0.; + res_q_lv_(trafo_id) = 0.; + // for voltages, because vm = 1. pu by hypothesis // res_v_hv_(trafo_id) = bus_vn_kv_hv; // res_v_lv_(trafo_id) = bus_vn_kv_lv; @@ -353,14 +358,16 @@ void TrafoContainer::compute_results(const Eigen::Ref & Va, } void TrafoContainer::reset_results(){ - res_p_hv_ = RealVect(); // in MW - res_q_hv_ = RealVect(); // in MVar - res_v_hv_ = RealVect(); // in kV - res_a_hv_ = RealVect(); // in kA - res_p_lv_ = RealVect(); // in MW - res_q_lv_ = RealVect(); // in MVar - res_v_lv_ = RealVect(); // in kV - res_a_lv_ = RealVect(); // in kA + res_p_hv_ = RealVect(nb()); // in MW + res_q_hv_ = RealVect(nb()); // in MVar + res_v_hv_ = RealVect(nb()); // in kV + res_a_hv_ = RealVect(nb()); // in kA + res_p_lv_ = RealVect(nb()); // in MW + res_q_lv_ = RealVect(nb()); // in MVar + res_v_lv_ = RealVect(nb()); // in kV + res_a_lv_ = RealVect(nb()); // in kA + res_theta_hv_ = RealVect(nb()); + res_theta_lv_ = RealVect(nb()); } diff --git a/src/main.cpp b/src/main.cpp index 49e8e88..a73fad7 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -817,6 +817,8 @@ PYBIND11_MODULE(lightsim2grid_cpp, m) .def("get_trafohv_theta", &GridModel::get_trafohv_theta, DocGridModel::_internal_do_not_use.c_str()) .def("get_trafolv_theta", &GridModel::get_trafolv_theta, DocGridModel::_internal_do_not_use.c_str()) + .def("get_all_shunt_buses", &GridModel::get_all_shunt_buses, DocGridModel::_internal_do_not_use.c_str()) + // do something with the grid // .def("init_Ybus", &DataModel::init_Ybus) // temporary .def("deactivate_result_computation", &GridModel::deactivate_result_computation, DocGridModel::deactivate_result_computation.c_str())