diff --git a/CHANGELOG.rst b/CHANGELOG.rst index 9fbafbd..c117ac5 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -23,6 +23,11 @@ Change Log - [FIXED] CI was broken after migration to artifact v4, set it back to v3 (and make the names of the folder clearer) - [FIXED] CI when using latest pandapower version (2.14) which broke some previous tests +- [ADDED] the computation of the LODF (line outage distribution factor) in + lightsim2grid +- [ADDED] some convenience functions to retrieve in a vectorized way the + buses to which each elements of a given container is connected + (*eg* `gridmodel.get_lines().get_bus_from()`) - [IMPROVED] remove some compilation warnings for clang - [IMPROVED] possibility to specify generator used as slack by its name when initializing from `pypowsybl`. diff --git a/lightsim2grid/tests/test_lodf.py b/lightsim2grid/tests/test_lodf.py new file mode 100644 index 0000000..1571c08 --- /dev/null +++ b/lightsim2grid/tests/test_lodf.py @@ -0,0 +1,156 @@ +# Copyright (c) 2023, RTE (https://www.rte-france.com) +# See AUTHORS.txt +# This Source Code Form is subject to the terms of the Mozilla Public License, version 2.0. +# If a copy of the Mozilla Public License, version 2.0 was not distributed with this file, +# you can obtain one at http://mozilla.org/MPL/2.0/. +# SPDX-License-Identifier: MPL-2.0 +# This file is part of LightSim2grid, LightSim2grid implements a c++ backend targeting the Grid2Op platform. + +import unittest +import pandapower as pp +import pandapower.networks as pn +import numpy as np +import scipy +import warnings + +from pandapower.pypower.makeLODF import update_LODF_diag +from lightsim2grid.gridmodel import init +from lightsim2grid.solver import SolverType + +import pdb + +class TestLODFCase14SLU(unittest.TestCase): + def make_grid(self): + case14 = pn.case14() + return case14 + + def get_solver_type(self): + return SolverType.DC + + def setUp(self) -> None: + self.case = self.make_grid() + with warnings.catch_warnings(): + warnings.filterwarnings("ignore") + self.gridmodel = init(self.case) + self.V_init = 1. * self.gridmodel.get_bus_vn_kv() + solver_type = self.get_solver_type() + if solver_type not in self.gridmodel.available_solvers(): + self.skipTest("Solver type not supported on this platform") + self.gridmodel.change_solver(solver_type) + V = self.gridmodel.dc_pf(1. * self.V_init, 1, 1e-8) + assert len(V), f"dc pf has diverged with error {self.gridmodel.get_dc_solver().get_error()}" + self.dcYbus = 1.0 * self.gridmodel.get_dcYbus() + self.dcSbus = 1.0 * self.gridmodel.get_dcSbus().real + self.Bbus = 1.0 * self.dcYbus.real + self.res_powerflow = 1.0 * np.concatenate((self.gridmodel.get_lineor_res()[0], self.gridmodel.get_trafohv_res()[0])) + self.nb = self.case.bus.shape[0] + self.nbr = self.case.line.shape[0] + self.case.trafo.shape[0] + self.slack_bus = self.case.ext_grid.iloc[0]["bus"] + self.noref = np.arange(1, self.nb) ## use bus 1 for voltage angle reference + self.noslack = np.flatnonzero(np.arange(self.nb) != self.slack_bus) + self.tol = 1e-6 + return super().setUp() + + def test_from_PTDF(self): + """compare the ls implementation pypower (in pandapowerer) + implementation and see if they match""" + PTDF = 1.0 * self.gridmodel.get_ptdf() + + # pypower implementation + f_ = np.concatenate((1 * self.gridmodel.get_lines().get_bus_from(), 1 * self.gridmodel.get_trafos().get_bus_from())) + t_ = np.concatenate((1 * self.gridmodel.get_lines().get_bus_to(), 1 * self.gridmodel.get_trafos().get_bus_to())) + Cft = scipy.sparse.csr_matrix((np.r_[np.ones(self.nbr), -np.ones(self.nbr)], + (np.r_[f_, t_], np.r_[np.arange(self.nbr), np.arange(self.nbr)])), + (self.nb, self.nbr)) + + H = PTDF * Cft + h = np.diag(H, 0) + den = (np.ones((self.nbr, 1)) * h.T * -1 + 1.) + with np.errstate(divide='ignore', invalid='ignore'): + LODF_pypower = (H / den) + update_LODF_diag(LODF_pypower) + # end pypower implementation + + LODF_ls = 1.0 * self.gridmodel.get_lodf() + + for l_id in range(LODF_ls.shape[0]): + pypow = 1. * LODF_pypower[:,l_id] + ls = 1. * LODF_ls[:,l_id] + if np.all(np.isfinite(pypow)): + # if everything works on pyower, it should work on lightsim2grid too + assert np.all(np.isfinite(ls)), f"error for line id {l_id}" + if np.all(np.isfinite(ls)): + # and the opposite too + assert np.all(np.isfinite(pypow)), f"error for line id {l_id}" + if np.any(~np.isfinite(pypow)): + # if it does not work for a given line, then + # every other flow are nan + assert np.all(~np.isfinite(pypow)), f"error for line id {l_id}" + # every flow for ls are nans too + assert np.all(~np.isfinite(ls)), f"error for line id {l_id}" + if np.any(~np.isfinite(ls)): + # if it does not work for a given line, then + # every other flow are nan + assert np.all(~np.isfinite(ls)), f"error for line id {l_id}" + # every flow for pypower are nans too + assert np.all(~np.isfinite(pypow)), f"error for line id {l_id}" + if np.all(np.isfinite(pypow)): + assert np.abs(pypow - ls).max() <= 1e-6, f"error for line id {l_id}" + + def test_with_powerflow(self): + """compute powerflow with LODF and with "standard" powerflow + and see if it matches + """ + LODF_ls = 1.0 * self.gridmodel.get_lodf() + nb_powerlines = len(self.gridmodel.get_lines()) + for l_id in range(LODF_ls.shape[0]): + if l_id < nb_powerlines: + self.gridmodel.deactivate_powerline(l_id) + else: + self.gridmodel.deactivate_trafo(l_id - nb_powerlines) + + V = self.gridmodel.dc_pf(1. * self.V_init, 1, 1e-8) + if V.shape[0] > 0: + # it has converged + por_pow = np.concatenate((self.gridmodel.get_lineor_res()[0], + self.gridmodel.get_trafohv_res()[0])) + por_lodf = self.res_powerflow + LODF_ls[:, l_id] * self.res_powerflow[l_id] + assert np.abs(por_lodf - por_pow).max() <= 1e-6, f"error for line id {l_id}" + else: + # it has diverged + assert np.all(~np.isfinite(LODF_ls[:,l_id])), f"error for line id {l_id}" + + if l_id < nb_powerlines: + self.gridmodel.reactivate_powerline(l_id) + else: + self.gridmodel.reactivate_trafo(l_id - nb_powerlines) + +class TestLODFCase30SLU(TestLODFCase14SLU): + def make_grid(self): + res = pn.case30() + return res + + +class TestLODFCase118SLU(TestLODFCase14SLU): + def make_grid(self): + res = pn.case118() + return res + + +class TestLODFCase14KLU(TestLODFCase14SLU): + def get_solver_type(self): + return SolverType.KLUDC + + +class TestLODFCase30KLU(TestLODFCase30SLU): + def get_solver_type(self): + return SolverType.KLUDC + + +class TestLODFCase118KLU(TestLODFCase118SLU): + def get_solver_type(self): + return SolverType.KLUDC + + +if __name__ == "__main__": + unittest.main() diff --git a/src/element_container/DCLineContainer.cpp b/src/element_container/DCLineContainer.cpp index 7992598..6c00533 100644 --- a/src/element_container/DCLineContainer.cpp +++ b/src/element_container/DCLineContainer.cpp @@ -66,8 +66,8 @@ void DCLineContainer::init(const Eigen::VectorXi & branch_from_id, void DCLineContainer::nb_line_end(std::vector & res) const { const Eigen::Index nb = from_gen_.nb(); - const auto & bus_or_id = get_bus_id_or(); - const auto & bus_ex_id = get_bus_id_ex(); + const auto & bus_or_id = get_bus_from(); + const auto & bus_ex_id = get_bus_to(); for(Eigen::Index i = 0; i < nb; ++i){ if(!status_[i]) continue; auto bus_or = bus_or_id(i); @@ -81,8 +81,8 @@ void DCLineContainer::nb_line_end(std::vector & res) const void DCLineContainer::disconnect_if_not_in_main_component(std::vector & busbar_in_main_component) { const Eigen::Index nb = from_gen_.nb(); - const auto & bus_or_id = get_bus_id_or(); - const auto & bus_ex_id = get_bus_id_ex(); + const auto & bus_or_id = get_bus_from(); + const auto & bus_ex_id = get_bus_to(); SolverControl unused_solver_control; for(Eigen::Index i = 0; i < nb; ++i){ if(!status_[i]) continue; diff --git a/src/element_container/DCLineContainer.h b/src/element_container/DCLineContainer.h index 400c1a5..c244987 100644 --- a/src/element_container/DCLineContainer.h +++ b/src/element_container/DCLineContainer.h @@ -289,8 +289,6 @@ class DCLineContainer : public GenericContainer } const std::vector& get_status() const {return status_;} - const Eigen::VectorXi & get_bus_id_or() const {return from_gen_.get_bus_id();} - const Eigen::VectorXi & get_bus_id_ex() const {return to_gen_.get_bus_id();} tuple3d get_or_res() const {return from_gen_.get_res();} tuple3d get_ex_res() const {return to_gen_.get_res();} @@ -299,6 +297,8 @@ class DCLineContainer : public GenericContainer Eigen::Ref get_theta_or() const {return from_gen_.get_theta();} Eigen::Ref get_theta_ex() const {return to_gen_.get_theta();} + Eigen::Ref get_bus_from() const {return from_gen_.get_bus_id();} + Eigen::Ref get_bus_to() const {return to_gen_.get_bus_id();} protected: // it is modeled as 2 generators that are "linked" together diff --git a/src/element_container/GeneratorContainer.h b/src/element_container/GeneratorContainer.h index 8f5e0a0..299311c 100644 --- a/src/element_container/GeneratorContainer.h +++ b/src/element_container/GeneratorContainer.h @@ -328,7 +328,7 @@ class GeneratorContainer: public GenericContainer Eigen::Ref get_theta() const {return res_theta_;} const std::vector& get_status() const {return status_;} - const Eigen::VectorXi & get_bus_id() const {return bus_id_;} + Eigen::Ref get_bus_id() const {return bus_id_;} void cout_v(){ for(const auto & el : vm_pu_){ diff --git a/src/element_container/LoadContainer.h b/src/element_container/LoadContainer.h index 0dc86cd..9a6710b 100644 --- a/src/element_container/LoadContainer.h +++ b/src/element_container/LoadContainer.h @@ -181,7 +181,7 @@ class LoadContainer : public GenericContainer Eigen::Ref get_theta() const {return res_theta_;} const std::vector& get_status() const {return status_;} - const Eigen::VectorXi & get_bus_id() const {return bus_id_;} + Eigen::Ref get_bus_id() const {return bus_id_;} protected: // physical properties diff --git a/src/element_container/SGenContainer.h b/src/element_container/SGenContainer.h index 6a3f831..0540147 100644 --- a/src/element_container/SGenContainer.h +++ b/src/element_container/SGenContainer.h @@ -201,7 +201,7 @@ class SGenContainer: public GenericContainer tuple4d get_res_full() const {return tuple4d(res_p_, res_q_, res_v_, res_theta_);} Eigen::Ref get_theta() const {return res_theta_;} const std::vector& get_status() const {return status_;} - const Eigen::VectorXi & get_bus_id() const {return bus_id_;} + Eigen::Ref get_bus_id() const {return bus_id_;} protected: // physical properties diff --git a/src/element_container/ShuntContainer.h b/src/element_container/ShuntContainer.h index 59c4a8a..e0fd79c 100644 --- a/src/element_container/ShuntContainer.h +++ b/src/element_container/ShuntContainer.h @@ -185,6 +185,7 @@ class ShuntContainer : public GenericContainer tuple4d get_res_full() const {return tuple4d(res_p_, res_q_, res_v_, res_theta_);} Eigen::Ref get_theta() const {return res_theta_;} + Eigen::Ref get_bus_id() const {return bus_id_;} const std::vector& get_status() const {return status_;} protected: diff --git a/src/main.cpp b/src/main.cpp index e3eba52..b88e01e 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -453,9 +453,11 @@ PYBIND11_MODULE(lightsim2grid_cpp, m) py::class_(m, "GeneratorContainer", DocIterator::GeneratorContainer.c_str()) .def("__len__", [](const GeneratorContainer & data) { return data.nb(); }) .def("__getitem__", [](const GeneratorContainer & data, int k){return data[k]; } ) - .def("__iter__", [](const GeneratorContainer & data) { - return py::make_iterator(data.begin(), data.end()); - }, py::keep_alive<0, 1>()); /* Keep vector alive while iterator is used */ + .def("__iter__", [](const GeneratorContainer & data) { + return py::make_iterator(data.begin(), data.end()); + }, py::keep_alive<0, 1>()) /* Keep vector alive while iterator is used */ + .def("get_bus_id", &GeneratorContainer::get_bus_id, "TODO doc", py::keep_alive<0, 1>()) + ; py::class_(m, "GenInfo", DocIterator::GenInfo.c_str()) .def_readonly("id", &GeneratorContainer::GenInfo::id, DocIterator::id.c_str()) @@ -481,8 +483,10 @@ PYBIND11_MODULE(lightsim2grid_cpp, m) .def("__len__", [](const SGenContainer & data) { return data.nb(); }) .def("__getitem__", [](const SGenContainer & data, int k){return data[k]; } ) .def("__iter__", [](const SGenContainer & data) { - return py::make_iterator(data.begin(), data.end()); - }, py::keep_alive<0, 1>()); /* Keep vector alive while iterator is used */ + return py::make_iterator(data.begin(), data.end()); + }, py::keep_alive<0, 1>()) /* Keep vector alive while iterator is used */ + .def("get_bus_id", &SGenContainer::get_bus_id, "TODO doc", py::keep_alive<0, 1>()) + ; py::class_(m, "SGenInfo", DocIterator::SGenInfo.c_str()) .def_readonly("id", &SGenContainer::SGenInfo::id, DocIterator::id.c_str()) @@ -506,8 +510,10 @@ PYBIND11_MODULE(lightsim2grid_cpp, m) .def("__len__", [](const LoadContainer & data) { return data.nb(); }) .def("__getitem__", [](const LoadContainer & data, int k){return data[k]; } ) .def("__iter__", [](const LoadContainer & data) { - return py::make_iterator(data.begin(), data.end()); - }, py::keep_alive<0, 1>()); /* Keep vector alive while iterator is used */ + return py::make_iterator(data.begin(), data.end()); + }, py::keep_alive<0, 1>()) /* Keep vector alive while iterator is used */ + .def("get_bus_id", &LoadContainer::get_bus_id, "TODO doc", py::keep_alive<0, 1>()) + ; py::class_(m, "LoadInfo", DocIterator::LoadInfo.c_str()) .def_readonly("id", &LoadContainer::LoadInfo::id, DocIterator::id.c_str()) @@ -527,8 +533,10 @@ PYBIND11_MODULE(lightsim2grid_cpp, m) .def("__len__", [](const ShuntContainer & data) { return data.nb(); }) .def("__getitem__", [](const ShuntContainer & data, int k){return data[k]; } ) .def("__iter__", [](const ShuntContainer & data) { - return py::make_iterator(data.begin(), data.end()); - }, py::keep_alive<0, 1>()); /* Keep vector alive while iterator is used */ + return py::make_iterator(data.begin(), data.end()); + }, py::keep_alive<0, 1>()) /* Keep vector alive while iterator is used */ + .def("get_bus_id", &ShuntContainer::get_bus_id, "TODO doc", py::keep_alive<0, 1>()) + ; py::class_(m, "ShuntInfo", DocIterator::ShuntInfo.c_str()) .def_readonly("id", &ShuntContainer::ShuntInfo::id, DocIterator::id.c_str()) @@ -548,8 +556,11 @@ PYBIND11_MODULE(lightsim2grid_cpp, m) .def("__len__", [](const TrafoContainer & data) { return data.nb(); }) .def("__getitem__", [](const TrafoContainer & data, int k){return data[k]; } ) .def("__iter__", [](const TrafoContainer & data) { - return py::make_iterator(data.begin(), data.end()); - }, py::keep_alive<0, 1>()); /* Keep vector alive while iterator is used */ + return py::make_iterator(data.begin(), data.end()); + }, py::keep_alive<0, 1>()) /* Keep vector alive while iterator is used */ + .def("get_bus_from", &TrafoContainer::get_bus_from, "TODO doc", py::keep_alive<0, 1>()) + .def("get_bus_to", &TrafoContainer::get_bus_to, "TODO doc", py::keep_alive<0, 1>()) + ; py::class_(m, "TrafoInfo", DocIterator::TrafoInfo.c_str()) .def_readonly("id", &TrafoContainer::TrafoInfo::id, DocIterator::id.c_str()) @@ -580,8 +591,10 @@ PYBIND11_MODULE(lightsim2grid_cpp, m) .def("__len__", [](const LineContainer & data) { return data.nb(); }) .def("__getitem__", [](const LineContainer & data, int k){return data[k]; } ) .def("__iter__", [](const LineContainer & data) { - return py::make_iterator(data.begin(), data.end()); - }, py::keep_alive<0, 1>()); /* Keep vector alive while iterator is used */ + return py::make_iterator(data.begin(), data.end()); + }, py::keep_alive<0, 1>()) /* Keep vector alive while iterator is used */ + .def("get_bus_from", &LineContainer::get_bus_from, "TODO doc", py::keep_alive<0, 1>()) + .def("get_bus_to", &LineContainer::get_bus_to, "TODO doc", py::keep_alive<0, 1>()); py::class_(m, "LineInfo", DocIterator::LineInfo.c_str()) .def_readonly("id", &LineContainer::LineInfo::id, DocIterator::id.c_str()) @@ -611,8 +624,10 @@ PYBIND11_MODULE(lightsim2grid_cpp, m) .def("__len__", [](const DCLineContainer & data) { return data.nb(); }) .def("__getitem__", [](const DCLineContainer & data, int k){return data[k]; } ) .def("__iter__", [](const DCLineContainer & data) { - return py::make_iterator(data.begin(), data.end()); - }, py::keep_alive<0, 1>()); /* Keep vector alive while iterator is used */ + return py::make_iterator(data.begin(), data.end()); + }, py::keep_alive<0, 1>()) /* Keep vector alive while iterator is used */ + .def("get_bus_from", &DCLineContainer::get_bus_from) + .def("get_bus_to", &DCLineContainer::get_bus_to); py::class_(m, "DCLineInfo", DocIterator::DCLineInfo.c_str()) .def_readonly("id", &DCLineContainer::DCLineInfo::id, DocIterator::id.c_str()) diff --git a/src/powerflow_algorithm/BaseDCAlgo.tpp b/src/powerflow_algorithm/BaseDCAlgo.tpp index 507af77..3e59cfb 100644 --- a/src/powerflow_algorithm/BaseDCAlgo.tpp +++ b/src/powerflow_algorithm/BaseDCAlgo.tpp @@ -7,6 +7,8 @@ // This file is part of LightSim2grid, LightSim2grid implements a c++ backend targeting the Grid2Op platform. // #include "DCSolver.h" +#include // for nans +#include // for nans // TODO SLACK !!! template @@ -281,10 +283,20 @@ RealMat BaseDCAlgo::get_ptdf(const Eigen::SparseMatrix template RealMat BaseDCAlgo::get_lodf(const Eigen::SparseMatrix & dcYbus, - const IntVect & from_bus, - const IntVect & to_bus){ - auto PTDF = get_ptdf(dcYbus); // size n_bus x n_bus + const IntVect & from_bus, + const IntVect & to_bus){ + const RealMat PTDF = get_ptdf(dcYbus); // size n_line x n_bus RealMat LODF = RealMat::Zero(from_bus.size(), from_bus.rows()); // nb_line, nb_line + for(Eigen::Index line_id=0; line_id < from_bus.size(); ++line_id){ + LODF.col(line_id).array() = PTDF.col(from_bus(line_id)).array() - PTDF.col(to_bus(line_id)).array(); + const real_type diag_coeff = LODF(line_id, line_id); + if (diag_coeff != 1.){ + LODF.col(line_id).array() /= (1. - diag_coeff); + LODF(line_id, line_id) = -1.; + }else{ + LODF.col(line_id).array() = std::numeric_limits::quiet_NaN(); + } + } return LODF; }