From b465c2b26f2ddfe2011990ca11990d1b554ec723 Mon Sep 17 00:00:00 2001 From: DONNOT Benjamin Date: Thu, 7 Mar 2024 17:58:41 +0100 Subject: [PATCH] fixing bug in the benchmarking of grid size, improve solver control --- .circleci/config.yml | 15 ++- .readthedocs.yml | 11 +- CHANGELOG.rst | 1 + benchmarks/benchmark_grid_size.py | 113 +++++++++++++----- lightsim2grid/lightSimBackend.py | 2 + src/element_container/GeneratorContainer.cpp | 30 +++-- src/element_container/LoadContainer.cpp | 24 ++-- src/element_container/SGenContainer.cpp | 40 ++++--- src/element_container/ShuntContainer.cpp | 19 ++- src/powerflow_algorithm/BaseDCAlgo.tpp | 5 + src/powerflow_algorithm/BaseFDPFAlgo.tpp | 51 +++++--- src/powerflow_algorithm/BaseNRAlgo.tpp | 28 +++-- .../BaseNRSingleSlackAlgo.tpp | 45 ++++--- 13 files changed, 259 insertions(+), 125 deletions(-) diff --git a/.circleci/config.yml b/.circleci/config.yml index 6973ecc..4ea5c7c 100644 --- a/.circleci/config.yml +++ b/.circleci/config.yml @@ -46,6 +46,7 @@ jobs: steps: - checkout - run: apt-get update && apt-get install python3-pip -y + - run: python3 -m pip install --upgrade pip setuptools - run: python3 -m pip install virtualenv - run: python3 -m virtualenv venv_test - run: @@ -64,6 +65,7 @@ jobs: steps: - checkout - run: apt-get update && apt-get install python3-pip -y + - run: python3 -m pip install --upgrade pip setuptools - run: python3 -m pip install virtualenv - run: python3 -m virtualenv venv_test - run: @@ -82,6 +84,7 @@ jobs: steps: - checkout - run: apt-get update && apt-get install python3-pip -y + - run: python3 -m pip install --upgrade pip setuptools - run: python3 -m pip install virtualenv - run: python3 -m virtualenv venv_test - run: @@ -100,6 +103,7 @@ jobs: steps: - checkout - run: apt-get update && apt-get install python3-pip -y + - run: python3 -m pip install --upgrade pip setuptools - run: python3 -m pip install virtualenv - run: python3 -m virtualenv venv_test - run: @@ -118,6 +122,7 @@ jobs: steps: - checkout - run: apt-get update && apt-get install python3-pip python3-virtualenv -y + - run: python3 -m pip install --upgrade pip setuptools # - run: python3 -m pip install virtualenv - run: python3 -m virtualenv venv_test - run: @@ -148,6 +153,7 @@ jobs: steps: - checkout - run: apt-get update && apt-get install python3-pip -y + - run: python3 -m pip install --upgrade pip setuptools - run: python3 -m pip install virtualenv - run: python3 -m virtualenv venv_test - run: @@ -165,6 +171,7 @@ jobs: steps: - checkout - run: apt-get update && apt-get install python3-pip git -y + - run: python3 -m pip install --upgrade pip setuptools - run: command: | git submodule init @@ -184,6 +191,7 @@ jobs: steps: - checkout - run: apt-get update && apt-get install python3-pip git -y + - run: python3 -m pip install --upgrade pip setuptools - run: python3 -m pip install virtualenv - run: python3 -m virtualenv venv_test - run: @@ -202,6 +210,7 @@ jobs: steps: - checkout - run: apt-get update && apt-get install python3-pip git -y + - run: python3 -m pip install --upgrade pip setuptools - run: python3 -m pip install virtualenv - run: python3 -m virtualenv venv_test - run: @@ -220,6 +229,7 @@ jobs: steps: - checkout - run: apt-get update && apt-get install python3-pip git -y + - run: python3 -m pip install --upgrade pip setuptools - run: python3 -m pip install virtualenv - run: python3 -m virtualenv venv_test - run: @@ -238,6 +248,7 @@ jobs: steps: - checkout - run: apt-get update && apt-get install python3-pip git -y + - run: python3 -m pip install --upgrade pip setuptools - run: python3 -m pip install virtualenv - run: python3 -m virtualenv venv_test - run: @@ -268,9 +279,7 @@ jobs: size: medium # ("medium" "large" "xlarge" "2xlarge") steps: - checkout - # - run: - # name: "Install Python" - # command: choco install python --version=3.9 # use python 3.9 for windows test + - run: python3 -m pip install --upgrade pip setuptools - run: py -m pip install virtualenv - run: py -m virtualenv venv_test - run: diff --git a/.readthedocs.yml b/.readthedocs.yml index 97df330..ab6ff9c 100644 --- a/.readthedocs.yml +++ b/.readthedocs.yml @@ -1,4 +1,4 @@ -version: 2 +version: "2" submodules: include: @@ -6,8 +6,15 @@ submodules: - eigen recursive: true +build: + os: "ubuntu-22.04" + tools: + python: "3.10" + +sphinx: + configuration: docs/conf.py + python: - version: 3.8 install: - method: pip path: . diff --git a/CHANGELOG.rst b/CHANGELOG.rst index 666a2a4..2fa05e1 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -39,6 +39,7 @@ Change Log - [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) +- [FIXED] read the docs was broken - [ADDED] sets of methods to extract the main component of a grid and perform powerflow only on this one. - [ADDED] possibility to set / retrieve the names of each elements of the grid. diff --git a/benchmarks/benchmark_grid_size.py b/benchmarks/benchmark_grid_size.py index 09d627e..a2e2668 100644 --- a/benchmarks/benchmark_grid_size.py +++ b/benchmarks/benchmark_grid_size.py @@ -7,12 +7,15 @@ # This file is part of LightSim2grid, LightSim2grid a implements a c++ backend targeting the Grid2Op platform. import warnings +import copy import pandapower as pp -import numpy as np +import numpy as np +import hashlib from scipy.interpolate import interp1d import matplotlib.pyplot as plt from grid2op import make, Parameters from grid2op.Chronics import FromNPY +from grid2op.Backend import PandaPowerBackend from lightsim2grid import LightSimBackend, TimeSerie try: from lightsim2grid import ContingencyAnalysis @@ -33,6 +36,8 @@ VERBOSE = False MAKE_PLOT = False +WITH_PP = False +DEBUG = False case_names = [ "case14.json", @@ -69,7 +74,28 @@ def make_grid2op_env(pp_case, casse_name, load_p, load_q, gen_p, sgen_p): ) return env_lightsim -def get_loads_gens(load_p_init, load_q_init, gen_p_init, sgen_p_init): + +def make_grid2op_env_pp(pp_case, casse_name, load_p, load_q, gen_p, sgen_p): + param = Parameters.Parameters() + param.init_from_dict({"NO_OVERFLOW_DISCONNECTION": True}) + + with warnings.catch_warnings(): + warnings.filterwarnings("ignore") + env_pp = make("blank", + param=param, test=True, + backend=PandaPowerBackend(lightsim2grid=False), + chronics_class=FromNPY, + data_feeding_kwargs={"load_p": load_p, + "load_q": load_q, + "prod_p": gen_p + }, + grid_path=case_name, + _add_to_name=f"{case_name}", + ) + return env_pp + + +def get_loads_gens(load_p_init, load_q_init, gen_p_init, sgen_p_init, prng): # scale loads # use some French time series data for loads @@ -137,28 +163,31 @@ def get_loads_gens(load_p_init, load_q_init, gen_p_init, sgen_p_init): vals = np.array(vals) * coeffs["month"]["oct"] * coeffs["day"]["mon"] x_interp = 12 * np.arange(len(vals)) coeffs = interp1d(x=x_interp, y=vals, kind="cubic") - all_vals = coeffs(x_final) + all_vals = coeffs(x_final).reshape(-1, 1) + if DEBUG: + all_vals[:] = 1 # compute the "smooth" loads matrix - load_p_smooth = all_vals.reshape(-1, 1) * load_p_init.reshape(1, -1) - load_q_smooth = all_vals.reshape(-1, 1) * load_q_init.reshape(1, -1) + load_p_smooth = all_vals * load_p_init.reshape(1, -1) + load_q_smooth = all_vals * load_q_init.reshape(1, -1) # add a bit of noise to it to get the "final" loads matrix - load_p = load_p_smooth * np.random.lognormal(mean=0., sigma=0.003, size=load_p_smooth.shape) - load_q = load_q_smooth * np.random.lognormal(mean=0., sigma=0.003, size=load_q_smooth.shape) - + load_p = load_p_smooth * prng.lognormal(mean=0., sigma=0.003, size=load_p_smooth.shape) + load_q = load_q_smooth * prng.lognormal(mean=0., sigma=0.003, size=load_q_smooth.shape) + if DEBUG: + load_p[:] = load_p_smooth + load_q[:] = load_q_smooth + # scale generators accordingly gen_p = load_p.sum(axis=1).reshape(-1, 1) / load_p_init.sum() * gen_p_init.reshape(1, -1) sgen_p = load_p.sum(axis=1).reshape(-1, 1) / load_p_init.sum() * sgen_p_init.reshape(1, -1) - return load_p, load_q, gen_p, sgen_p if __name__ == "__main__": - np.random.seed(42) - + prng = np.random.default_rng(42) case_names_displayed = [get_env_name_displayed(el) for el in case_names] - g2op_times = [] + solver_preproc_solver_time = [] g2op_speeds = [] g2op_sizes = [] g2op_step_time = [] @@ -197,8 +226,14 @@ def get_loads_gens(load_p_init, load_q_init, gen_p_init, sgen_p_init): res_unit = "ms" # simulate the data - load_p, load_q, gen_p, sgen_p = get_loads_gens(load_p_init, load_q_init, gen_p_init, sgen_p_init) - + load_p, load_q, gen_p, sgen_p = get_loads_gens(load_p_init, load_q_init, gen_p_init, sgen_p_init, prng) + if DEBUG: + hash_fun = hashlib.blake2b(digest_size=16) + hash_fun.update(load_p.tobytes()) + hash_fun.update(load_q.tobytes()) + hash_fun.update(gen_p.tobytes()) + hash_fun.update(sgen_p.tobytes()) + print(hash_fun.hexdigest()) # create the grid2op env nb_ts = gen_p.shape[0] # add slack ! @@ -206,15 +241,33 @@ def get_loads_gens(load_p_init, load_q_init, gen_p_init, sgen_p_init): if "res_ext_grid" in case: slack_gens += np.tile(case.res_ext_grid["p_mw"].values.reshape(1,-1), (nb_ts, 1)) gen_p_g2op = np.concatenate((gen_p, slack_gens), axis=1) - # get the env + # get the env + if WITH_PP: + env_pp = make_grid2op_env_pp(case, + case_name, + load_p, + load_q, + gen_p_g2op, + sgen_p) + _ = env_pp.reset() + done = False + nb_step_pp = 0 + changed_sgen = case.sgen["in_service"].values + while not done: + # hack for static gen... + env_pp.backend._grid.sgen["p_mw"] = sgen_p[nb_step_pp, :] + obs, reward, done, info = env_pp.step(env_pp.action_space()) + nb_step_pp += 1 + if nb_step_pp != nb_ts: + warnings.warn("Divergence even with pandapower !") + print("Pandapower stops, lightsim starts") + env_lightsim = make_grid2op_env(case, case_name, load_p, load_q, gen_p_g2op, sgen_p) - - # Perform the computation using grid2op _ = env_lightsim.reset() done = False @@ -222,16 +275,20 @@ def get_loads_gens(load_p_init, load_q_init, gen_p_init, sgen_p_init): changed_sgen = case.sgen["in_service"].values while not done: # hack for static gen... - env_lightsim.backend._grid.update_sgens_p(changed_sgen, - sgen_p[nb_step, changed_sgen].astype(np.float32)) + changed_sgen = copy.deepcopy(case.sgen["in_service"].values) + this_sgen = sgen_p[nb_step, :].astype(np.float32) + # this_sgen = sgen_p_init[changed_sgen].astype(np.float32) + env_lightsim.backend._grid.update_sgens_p(changed_sgen, this_sgen) obs, reward, done, info = env_lightsim.step(env_lightsim.action_space()) nb_step += 1 # NB lightsim2grid does not handle "static gen" because I cannot set "p" in gen in grid2op # so results will vary between TimeSeries and grid2op ! - + # env_lightsim.backend._grid.tell_solver_need_reset() + # env_lightsim.backend._grid.dc_pf(env_lightsim.backend.V, 1, 1e-7) + # env_lightsim.backend._grid.get_bus_status() if nb_step != nb_ts: warnings.warn(f"only able to make {nb_step} (out of {nb_ts}) for {case_name} in grid2op. Results will not be availabe for grid2op step") - g2op_times.append(None) + solver_preproc_solver_time.append(None) g2op_speeds.append(None) g2op_step_time.append(None) ls_solver_time.append(None) @@ -239,8 +296,8 @@ def get_loads_gens(load_p_init, load_q_init, gen_p_init, sgen_p_init): g2op_sizes.append(env_lightsim.n_sub) else: total_time = env_lightsim.backend._timer_preproc + env_lightsim.backend._timer_solver # + env_lightsim.backend._timer_postproc - total_time = env_lightsim._time_step - g2op_times.append(total_time) + # total_time = env_lightsim._time_step + solver_preproc_solver_time.append(total_time) g2op_speeds.append(1.0 * nb_step / total_time) g2op_step_time.append(1.0 * env_lightsim._time_step / nb_step) ls_solver_time.append(env_lightsim.backend.comp_time) @@ -262,7 +319,7 @@ def get_loads_gens(load_p_init, load_q_init, gen_p_init, sgen_p_init): env_lightsim.backend.tol) time_serie._TimeSerie__computed = True a_or = time_serie.compute_A() - assert status, f"some powerflow diverge for {case_name}: {computer.nb_solved()} " + assert status, f"some powerflow diverge for Time Series for {case_name}: {computer.nb_solved()} " if VERBOSE: # print detailed results if needed @@ -311,9 +368,9 @@ def get_loads_gens(load_p_init, load_q_init, gen_p_init, sgen_p_init): for i, nm_ in enumerate(case_names_displayed): tab_g2op.append((nm_, ts_sizes[i], + 1000. * g2op_step_time[i] if g2op_step_time[i] else None, 1000. / g2op_speeds[i] if g2op_speeds[i] else None, g2op_speeds[i], - 1000. * g2op_step_time[i] if g2op_step_time[i] else None, 1000. * ls_gridmodel_time[i] / nb_step if ls_gridmodel_time[i] else None, 1000. * ls_solver_time[i] / nb_step if ls_solver_time[i] else None, )) @@ -321,9 +378,9 @@ def get_loads_gens(load_p_init, load_q_init, gen_p_init, sgen_p_init): res_use_with_grid2op_2 = tabulate(tab_g2op, headers=["grid", "size (nb bus)", - "step time (ms / pf)", - "speed (pf / s)", "avg step duration (ms)", + "time [DC + AC] (ms / pf)", + "speed (pf / s)", "time in 'gridmodel' (ms / pf)", "time in 'pf algo' (ms / pf)", ], @@ -362,7 +419,7 @@ def get_loads_gens(load_p_init, load_q_init, gen_p_init, sgen_p_init): if MAKE_PLOT: # make the plot summarizing all results - plt.plot(g2op_sizes, g2op_times, linestyle='solid', marker='+', markersize=8) + plt.plot(g2op_sizes, solver_preproc_solver_time, linestyle='solid', marker='+', markersize=8) plt.xlabel("Size (number of substation)") plt.ylabel("Time taken (s)") plt.title(f"Time to compute {g2op_sizes[0]} powerflows using Grid2Op.step (dc pf [init] + ac pf)") diff --git a/lightsim2grid/lightSimBackend.py b/lightsim2grid/lightSimBackend.py index bbecd58..f0d57c5 100644 --- a/lightsim2grid/lightSimBackend.py +++ b/lightsim2grid/lightSimBackend.py @@ -415,6 +415,7 @@ def load_grid(self, path=None, filename=None): self._load_grid_pypowsybl(path, filename) else: raise BackendError(f"Impossible to initialize the backend with '{self._loader_method}'") + self._grid.tell_solver_need_reset() def _should_not_have_to_do_this(self, path=None, filename=None): # included in grid2op now ! @@ -1230,3 +1231,4 @@ def reset(self, grid_path, grid_filename=None): self._timer_postproc = 0. self._timer_preproc = 0. self._timer_solver = 0. + self._grid.tell_solver_need_reset() diff --git a/src/element_container/GeneratorContainer.cpp b/src/element_container/GeneratorContainer.cpp index e45767b..b580c40 100644 --- a/src/element_container/GeneratorContainer.cpp +++ b/src/element_container/GeneratorContainer.cpp @@ -167,9 +167,9 @@ void GeneratorContainer::fillSbus(CplxVect & Sbus, const std::vector & id_g } void GeneratorContainer::fillpv(std::vector & bus_pv, - std::vector & has_bus_been_added, - const Eigen::VectorXi & slack_bus_id_solver, - const std::vector & id_grid_to_solver) const + std::vector & has_bus_been_added, + const Eigen::VectorXi & slack_bus_id_solver, + const std::vector & id_grid_to_solver) const { const int nb_gen = nb(); int bus_id_me, bus_id_solver; @@ -198,12 +198,12 @@ void GeneratorContainer::fillpv(std::vector & bus_pv, } void GeneratorContainer::compute_results(const Eigen::Ref & Va, - const Eigen::Ref & Vm, - const Eigen::Ref & V, - const std::vector & id_grid_to_solver, - const RealVect & bus_vn_kv, - real_type sn_mva, - bool ac) + const Eigen::Ref & Vm, + const Eigen::Ref & V, + const std::vector & id_grid_to_solver, + const RealVect & bus_vn_kv, + real_type sn_mva, + bool ac) { const int nb_gen = nb(); v_kv_from_vpu(Va, Vm, status_, nb_gen, bus_id_, id_grid_to_solver, bus_vn_kv, res_v_); @@ -257,8 +257,10 @@ void GeneratorContainer::change_p(int gen_id, real_type new_p, SolverControl & s solver_control.tell_pv_changed(); } } - if (p_mw_(gen_id) != new_p) solver_control.tell_recompute_sbus(); - p_mw_(gen_id) = new_p; + if (p_mw_(gen_id) != new_p){ + solver_control.tell_recompute_sbus(); + p_mw_(gen_id) = new_p; + } } void GeneratorContainer::change_q(int gen_id, real_type new_q, SolverControl & solver_control) @@ -275,8 +277,10 @@ void GeneratorContainer::change_q(int gen_id, real_type new_q, SolverControl & s } // TODO DEBUG MODE : raise an error if generator is regulating voltage, maybe ? // this would have not effect - if (q_mvar_(gen_id) != new_q) solver_control.tell_recompute_sbus(); - q_mvar_(gen_id) = new_q; + if (q_mvar_(gen_id) != new_q){ + solver_control.tell_recompute_sbus(); + q_mvar_(gen_id) = new_q; + } } void GeneratorContainer::change_v(int gen_id, real_type new_v_pu, SolverControl & solver_control) diff --git a/src/element_container/LoadContainer.cpp b/src/element_container/LoadContainer.cpp index 4a4de71..d7eb561 100644 --- a/src/element_container/LoadContainer.cpp +++ b/src/element_container/LoadContainer.cpp @@ -71,12 +71,12 @@ void LoadContainer::fillSbus(CplxVect & Sbus, const std::vector & id_grid_t } void LoadContainer::compute_results(const Eigen::Ref & Va, - const Eigen::Ref & Vm, - const Eigen::Ref & V, - const std::vector & id_grid_to_solver, - const RealVect & bus_vn_kv, - real_type sn_mva, - bool ac) + const Eigen::Ref & Vm, + const Eigen::Ref & V, + const std::vector & id_grid_to_solver, + const RealVect & bus_vn_kv, + real_type sn_mva, + bool ac) { int nb_load = nb(); v_kv_from_vpu(Va, Vm, status_, nb_load, bus_id_, id_grid_to_solver, bus_vn_kv, res_v_); @@ -102,8 +102,10 @@ void LoadContainer::change_p(int load_id, real_type new_p, SolverControl & solve exc_ << ")"; throw std::runtime_error(exc_.str()); } - if (p_mw_(load_id) != new_p) solver_control.tell_recompute_sbus(); - p_mw_(load_id) = new_p; + if (p_mw_(load_id) != new_p) { + solver_control.tell_recompute_sbus(); + p_mw_(load_id) = new_p; + } } void LoadContainer::change_q(int load_id, real_type new_q, SolverControl & solver_control) @@ -117,8 +119,10 @@ void LoadContainer::change_q(int load_id, real_type new_q, SolverControl & solve exc_ << ")"; throw std::runtime_error(exc_.str()); } - if (q_mvar_(load_id) != new_q) solver_control.tell_recompute_sbus(); - q_mvar_(load_id) = new_q; + if (q_mvar_(load_id) != new_q) { + solver_control.tell_recompute_sbus(); + q_mvar_(load_id) = new_q; + } } void LoadContainer::reconnect_connected_buses(std::vector & bus_status) const { diff --git a/src/element_container/SGenContainer.cpp b/src/element_container/SGenContainer.cpp index 48c1baa..6300727 100644 --- a/src/element_container/SGenContainer.cpp +++ b/src/element_container/SGenContainer.cpp @@ -7,14 +7,15 @@ // This file is part of LightSim2grid, LightSim2grid implements a c++ backend targeting the Grid2Op platform. #include "SGenContainer.h" +#include void SGenContainer::init(const RealVect & sgen_p, - const RealVect & sgen_q, - const RealVect & sgen_pmin, - const RealVect & sgen_pmax, - const RealVect & sgen_qmin, - const RealVect & sgen_qmax, - const Eigen::VectorXi & sgen_bus_id) + const RealVect & sgen_q, + const RealVect & sgen_pmin, + const RealVect & sgen_pmax, + const RealVect & sgen_qmin, + const RealVect & sgen_qmax, + const Eigen::VectorXi & sgen_bus_id) { int size = static_cast(sgen_p.size()); GenericContainer::check_size(sgen_p, size, "sgen_p"); @@ -100,19 +101,18 @@ void SGenContainer::fillSbus(CplxVect & Sbus, const std::vector & id_grid_t exc_ << " is connected to a disconnected bus while being connected to the grid."; throw std::runtime_error(exc_.str()); } - tmp = static_cast(p_mw_(sgen_id)); - tmp += my_i * q_mvar_(sgen_id); + tmp = {p_mw_(sgen_id), q_mvar_(sgen_id)}; Sbus.coeffRef(bus_id_solver) += tmp; } } void SGenContainer::compute_results(const Eigen::Ref & Va, - const Eigen::Ref & Vm, - const Eigen::Ref & V, - const std::vector & id_grid_to_solver, - const RealVect & bus_vn_kv, - real_type sn_mva, - bool ac) + const Eigen::Ref & Vm, + const Eigen::Ref & V, + const std::vector & id_grid_to_solver, + const RealVect & bus_vn_kv, + real_type sn_mva, + bool ac) { const int nb_sgen = nb(); v_kv_from_vpu(Va, Vm, status_, nb_sgen, bus_id_, id_grid_to_solver, bus_vn_kv, res_v_); @@ -139,8 +139,10 @@ void SGenContainer::change_p(int sgen_id, real_type new_p, SolverControl & solve exc_ << ")"; throw std::runtime_error(exc_.str()); } - if (p_mw_(sgen_id) != new_p) solver_control.tell_recompute_sbus(); - p_mw_(sgen_id) = new_p; + if (p_mw_(sgen_id) != new_p){ + solver_control.tell_recompute_sbus(); + p_mw_(sgen_id) = new_p; + } } void SGenContainer::change_q(int sgen_id, real_type new_q, SolverControl & solver_control) @@ -154,8 +156,10 @@ void SGenContainer::change_q(int sgen_id, real_type new_q, SolverControl & solve exc_ << ")"; throw std::runtime_error(exc_.str()); } - if (q_mvar_(sgen_id) != new_q) solver_control.tell_recompute_sbus(); - q_mvar_(sgen_id) = new_q; + if (q_mvar_(sgen_id) != new_q){ + solver_control.tell_recompute_sbus(); + q_mvar_(sgen_id) = new_q; + } } void SGenContainer::reconnect_connected_buses(std::vector & bus_status) const { diff --git a/src/element_container/ShuntContainer.cpp b/src/element_container/ShuntContainer.cpp index 5817834..76a6129 100644 --- a/src/element_container/ShuntContainer.cpp +++ b/src/element_container/ShuntContainer.cpp @@ -48,9 +48,9 @@ void ShuntContainer::set_state(ShuntContainer::StateRes & my_state ) } void ShuntContainer::fillYbus(std::vector > & res, - bool ac, - const std::vector & id_grid_to_solver, - real_type sn_mva) const + bool ac, + const std::vector & id_grid_to_solver, + real_type sn_mva) const { if(!ac) return; // no shunt in DC @@ -79,10 +79,10 @@ void ShuntContainer::fillYbus(std::vector > & res, } void ShuntContainer::fillBp_Bpp(std::vector > & Bp, - std::vector > & Bpp, - const std::vector & id_grid_to_solver, - real_type sn_mva, - FDPFMethod xb_or_bx) const + std::vector > & Bpp, + const std::vector & id_grid_to_solver, + real_type sn_mva, + FDPFMethod xb_or_bx) const { const Eigen::Index nb_shunt = static_cast(q_mvar_.size()); real_type tmp; @@ -174,9 +174,8 @@ void ShuntContainer::change_p(int shunt_id, real_type new_p, SolverControl & sol 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; } - p_mw_(shunt_id) = new_p; - } void ShuntContainer::change_q(int shunt_id, real_type new_q, SolverControl & solver_control) @@ -185,8 +184,8 @@ void ShuntContainer::change_q(int shunt_id, real_type new_q, SolverControl & sol 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_ybus(); + q_mvar_(shunt_id) = new_q; } - q_mvar_(shunt_id) = new_q; } void ShuntContainer::reconnect_connected_buses(std::vector & bus_status) const { diff --git a/src/powerflow_algorithm/BaseDCAlgo.tpp b/src/powerflow_algorithm/BaseDCAlgo.tpp index bb59466..7de831f 100644 --- a/src/powerflow_algorithm/BaseDCAlgo.tpp +++ b/src/powerflow_algorithm/BaseDCAlgo.tpp @@ -27,6 +27,7 @@ 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()) { + // std::cout << "!is_linear_solver_valid()\n"; return false; } BaseAlgo::reset_timer(); @@ -104,6 +105,7 @@ bool BaseDCAlgo::compute_pf(const Eigen::SparseMatrix & ErrorType status_init = _linear_solver.initialize(dcYbus_noslack_); if(status_init != ErrorType::NoError){ err_ = status_init; + // std::cout << "_linear_solver.initialize\n"; return false; } need_factorize_ = false; @@ -116,6 +118,7 @@ bool BaseDCAlgo::compute_pf(const Eigen::SparseMatrix & if(error != ErrorType::NoError){ err_ = error; timer_total_nr_ += timer.duration(); + // std::cout << "_linear_solver.solve\n"; return false; } @@ -128,8 +131,10 @@ bool BaseDCAlgo::compute_pf(const Eigen::SparseMatrix & Vm_ = RealVect(); Va_ = RealVect(); timer_total_nr_ += timer.duration(); + // std::cout << "_linear_solver.allFinite" << Va_dc_without_slack.array().allFinite() <<", " << Va_dc_without_slack.lpNorm() <<"\n"; return false; } + // std::cout << "\t " << Va_dc_without_slack.lpNorm() << "\n"; #ifdef __COUT_TIMES std::cout << "\t dc solve: " << 1000. * timer_solve.duration() << "ms" << std::endl; diff --git a/src/powerflow_algorithm/BaseFDPFAlgo.tpp b/src/powerflow_algorithm/BaseFDPFAlgo.tpp index 28a7c1b..1c063b3 100644 --- a/src/powerflow_algorithm/BaseFDPFAlgo.tpp +++ b/src/powerflow_algorithm/BaseFDPFAlgo.tpp @@ -67,22 +67,35 @@ bool BaseFDPFAlgo::compute_pf(const Eigen::SparseMatrix grid_Bp; - Eigen::SparseMatrix grid_Bpp; - fillBp_Bpp(grid_Bp, grid_Bpp); - - // fill the solver matrices Bp and Bpp : - // Bp_ = Bp[array([pvpq]).T, pvpq].tocsc() # splu requires a CSC matrix - // Bpp_ = Bpp[array([pq]).T, pq].tocsc() - const auto n_pvpq = pvpq.size(); - std::vector pvpq_inv(V.size(), -1); - for(int inv_id=0; inv_id < n_pvpq; ++inv_id) pvpq_inv[pvpq(inv_id)] = inv_id; - std::vector pq_inv(V.size(), -1); - for(int inv_id=0; inv_id < n_pq; ++inv_id) pq_inv[pq(inv_id)] = inv_id; - fill_sparse_matrices(grid_Bp, grid_Bpp, pvpq_inv, pq_inv, n_pvpq, n_pq); + if(need_factorize_ || + _solver_control.need_reset_solver() || + _solver_control.has_dimension_changed() || + _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() || + _solver_control.need_recompute_ybus() || + _solver_control.has_slack_participate_changed() || + _solver_control.has_pv_changed() || + _solver_control.has_pq_changed() + ) + { + // extract Bp and Bpp for the whole grid + Eigen::SparseMatrix grid_Bp; + Eigen::SparseMatrix grid_Bpp; + fillBp_Bpp(grid_Bp, grid_Bpp); + + // fill the solver matrices Bp_ and Bpp_ + // Bp_ = Bp[array([pvpq]).T, pvpq].tocsc() + // Bpp_ = Bpp[array([pq]).T, pq].tocsc() + std::vector pvpq_inv(V.size(), -1); + for(int inv_id=0; inv_id < n_pvpq; ++inv_id) pvpq_inv[pvpq(inv_id)] = inv_id; + std::vector pq_inv(V.size(), -1); + for(int inv_id=0; inv_id < n_pq; ++inv_id) pq_inv[pq(inv_id)] = inv_id; + fill_sparse_matrices(grid_Bp, grid_Bpp, pvpq_inv, pq_inv, n_pvpq, n_pq); + } V_ = V; // V = V0 Vm_ = V_.array().abs(); // Vm = abs(V) @@ -160,11 +173,11 @@ bool BaseFDPFAlgo::compute_pf(const Eigen::SparseMatrix void BaseFDPFAlgo::fill_sparse_matrices(const Eigen::SparseMatrix & grid_Bp, - const Eigen::SparseMatrix & grid_Bpp, - const std::vector & pvpq_inv, - const std::vector & pq_inv, - Eigen::Index n_pvpq, - Eigen::Index n_pq) + const Eigen::SparseMatrix & grid_Bpp, + const std::vector & pvpq_inv, + const std::vector & pq_inv, + Eigen::Index n_pvpq, + Eigen::Index n_pq) { /** Init Bp_ and Bpp_ such that: diff --git a/src/powerflow_algorithm/BaseNRAlgo.tpp b/src/powerflow_algorithm/BaseNRAlgo.tpp index 0977d7d..33ae20e 100644 --- a/src/powerflow_algorithm/BaseNRAlgo.tpp +++ b/src/powerflow_algorithm/BaseNRAlgo.tpp @@ -88,13 +88,27 @@ bool BaseNRAlgo::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; - value_map_.clear(); // TODO smarter solver: only needed if ybus has changed - // BaseNRAlgo::col_map_.clear(); // TODO smarter solver: only needed if ybus has changed - // BaseNRAlgo::row_map_.clear(); // TODO smarter solver: only needed if ybus has changed - dS_dVm_.resize(0,0); // TODO smarter solver: only needed if ybus has changed - dS_dVa_.resize(0,0); // TODO smarter solver: only needed if ybus has changed - // BaseNRAlgo::dS_dVm_.setZero(); // TODO smarter solver: only needed if ybus has changed - // BaseNRAlgo::dS_dVa_.setZero(); // TODO smarter solver: only needed if ybus has changed + if(need_factorize_ || + _solver_control.need_reset_solver() || + _solver_control.has_dimension_changed() || + _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() || + _solver_control.need_recompute_ybus() || + _solver_control.has_slack_participate_changed() || + _solver_control.has_pv_changed() || + _solver_control.has_pq_changed() + ) + { + value_map_.clear(); // TODO smarter solver: only needed if ybus has changed + // BaseNRAlgo::col_map_.clear(); // TODO smarter solver: only needed if ybus has changed + // BaseNRAlgo::row_map_.clear(); // TODO smarter solver: only needed if ybus has changed + dS_dVm_.resize(0,0); // TODO smarter solver: only needed if ybus has changed + dS_dVa_.resize(0,0); // TODO smarter solver: only needed if ybus has changed + // BaseNRAlgo::dS_dVm_.setZero(); // TODO smarter solver: only needed if ybus has changed + // BaseNRAlgo::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); diff --git a/src/powerflow_algorithm/BaseNRSingleSlackAlgo.tpp b/src/powerflow_algorithm/BaseNRSingleSlackAlgo.tpp index 8c0a429..612f963 100644 --- a/src/powerflow_algorithm/BaseNRSingleSlackAlgo.tpp +++ b/src/powerflow_algorithm/BaseNRSingleSlackAlgo.tpp @@ -11,15 +11,15 @@ template bool BaseNRSingleSlackAlgo::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 @@ -75,12 +75,27 @@ bool BaseNRSingleSlackAlgo::compute_pf(const Eigen::SparseMatrix::my_i; // otherwise it does not compile - BaseNRAlgo::value_map_.clear(); // TODO smarter solver: only needed if ybus has changed or pq changed or pv changed - BaseNRAlgo::dS_dVm_.resize(0,0); // TODO smarter solver: only needed if ybus has changed or pq changed or pv changed - BaseNRAlgo::dS_dVa_.resize(0,0); // TODO smarter solver: only needed if ybus has changed or pq changed or pv changed - // BaseNRAlgo::J_.setZero(); // TODO smarter solver: only needed if ybus has changed or pq changed or pv changed or ybus_some_coeffs_zero_ - // BaseNRAlgo::dS_dVm_.setZero(); // TODO smarter solver: only needed if ybus has changed - // BaseNRAlgo::dS_dVa_.setZero(); // TODO smarter solver: only needed if ybus has changed + if(BaseNRAlgo::need_factorize_ || + BaseNRAlgo::_solver_control.need_reset_solver() || + BaseNRAlgo::_solver_control.has_dimension_changed() || + BaseNRAlgo::_solver_control.has_slack_participate_changed() || // the full "ybus without slack" has changed, everything needs to be recomputed_solver_control.ybus_change_sparsity_pattern() + BaseNRAlgo::_solver_control.ybus_change_sparsity_pattern() || + BaseNRAlgo::_solver_control.has_ybus_some_coeffs_zero() || + BaseNRAlgo::_solver_control.need_recompute_ybus() || + // BaseNRAlgo:: _solver_control.has_slack_participate_changed() || + BaseNRAlgo::_solver_control.has_pv_changed() || + BaseNRAlgo::_solver_control.has_pq_changed() + ) + { + BaseNRAlgo::value_map_.clear(); // TODO smarter solver: only needed if ybus has changed + // BaseNRAlgo::col_map_.clear(); // TODO smarter solver: only needed if ybus has changed + // BaseNRAlgo::row_map_.clear(); // TODO smarter solver: only needed if ybus has changed + BaseNRAlgo::dS_dVm_.resize(0,0); // TODO smarter solver: only needed if ybus has changed + BaseNRAlgo::dS_dVa_.resize(0,0); // TODO smarter solver: only needed if ybus has changed + // BaseNRAlgo::dS_dVm_.setZero(); // TODO smarter solver: only needed if ybus has changed + // BaseNRAlgo::dS_dVa_.setZero(); // TODO smarter solver: only needed if ybus has changed + + } while ((!converged) & (BaseNRAlgo::nr_iter_ < max_iter)){ BaseNRAlgo::nr_iter_++; // std::cout << "\tnr_iter_ " << BaseNRAlgo::nr_iter_ << std::endl;