Skip to content

Commit

Permalink
adding lodf
Browse files Browse the repository at this point in the history
  • Loading branch information
BDonnot committed Apr 9, 2024
1 parent 7342285 commit fdf1869
Show file tree
Hide file tree
Showing 10 changed files with 216 additions and 27 deletions.
5 changes: 5 additions & 0 deletions CHANGELOG.rst
Original file line number Diff line number Diff line change
Expand Up @@ -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`.
Expand Down
156 changes: 156 additions & 0 deletions lightsim2grid/tests/test_lodf.py
Original file line number Diff line number Diff line change
@@ -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()
8 changes: 4 additions & 4 deletions src/element_container/DCLineContainer.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -66,8 +66,8 @@ void DCLineContainer::init(const Eigen::VectorXi & branch_from_id,
void DCLineContainer::nb_line_end(std::vector<int> & 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);
Expand All @@ -81,8 +81,8 @@ void DCLineContainer::nb_line_end(std::vector<int> & res) const
void DCLineContainer::disconnect_if_not_in_main_component(std::vector<bool> & 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;
Expand Down
4 changes: 2 additions & 2 deletions src/element_container/DCLineContainer.h
Original file line number Diff line number Diff line change
Expand Up @@ -289,8 +289,6 @@ class DCLineContainer : public GenericContainer
}

const std::vector<bool>& 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();}
Expand All @@ -299,6 +297,8 @@ class DCLineContainer : public GenericContainer

Eigen::Ref<const RealVect> get_theta_or() const {return from_gen_.get_theta();}
Eigen::Ref<const RealVect> get_theta_ex() const {return to_gen_.get_theta();}
Eigen::Ref<const Eigen::VectorXi> get_bus_from() const {return from_gen_.get_bus_id();}
Eigen::Ref<const Eigen::VectorXi> get_bus_to() const {return to_gen_.get_bus_id();}

protected:
// it is modeled as 2 generators that are "linked" together
Expand Down
2 changes: 1 addition & 1 deletion src/element_container/GeneratorContainer.h
Original file line number Diff line number Diff line change
Expand Up @@ -328,7 +328,7 @@ class GeneratorContainer: public GenericContainer
Eigen::Ref<const RealVect> get_theta() const {return res_theta_;}

const std::vector<bool>& get_status() const {return status_;}
const Eigen::VectorXi & get_bus_id() const {return bus_id_;}
Eigen::Ref<const Eigen::VectorXi> get_bus_id() const {return bus_id_;}

void cout_v(){
for(const auto & el : vm_pu_){
Expand Down
2 changes: 1 addition & 1 deletion src/element_container/LoadContainer.h
Original file line number Diff line number Diff line change
Expand Up @@ -181,7 +181,7 @@ class LoadContainer : public GenericContainer

Eigen::Ref<const RealVect> get_theta() const {return res_theta_;}
const std::vector<bool>& get_status() const {return status_;}
const Eigen::VectorXi & get_bus_id() const {return bus_id_;}
Eigen::Ref<const Eigen::VectorXi> get_bus_id() const {return bus_id_;}

protected:
// physical properties
Expand Down
2 changes: 1 addition & 1 deletion src/element_container/SGenContainer.h
Original file line number Diff line number Diff line change
Expand Up @@ -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<const RealVect> get_theta() const {return res_theta_;}
const std::vector<bool>& get_status() const {return status_;}
const Eigen::VectorXi & get_bus_id() const {return bus_id_;}
Eigen::Ref<const Eigen::VectorXi> get_bus_id() const {return bus_id_;}

protected:
// physical properties
Expand Down
1 change: 1 addition & 0 deletions src/element_container/ShuntContainer.h
Original file line number Diff line number Diff line change
Expand Up @@ -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<const RealVect> get_theta() const {return res_theta_;}
Eigen::Ref<const Eigen::VectorXi> get_bus_id() const {return bus_id_;}
const std::vector<bool>& get_status() const {return status_;}

protected:
Expand Down
45 changes: 30 additions & 15 deletions src/main.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -453,9 +453,11 @@ PYBIND11_MODULE(lightsim2grid_cpp, m)
py::class_<GeneratorContainer>(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_<GeneratorContainer::GenInfo>(m, "GenInfo", DocIterator::GenInfo.c_str())
.def_readonly("id", &GeneratorContainer::GenInfo::id, DocIterator::id.c_str())
Expand All @@ -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_<SGenContainer::SGenInfo>(m, "SGenInfo", DocIterator::SGenInfo.c_str())
.def_readonly("id", &SGenContainer::SGenInfo::id, DocIterator::id.c_str())
Expand All @@ -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_<LoadContainer::LoadInfo>(m, "LoadInfo", DocIterator::LoadInfo.c_str())
.def_readonly("id", &LoadContainer::LoadInfo::id, DocIterator::id.c_str())
Expand All @@ -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_<ShuntContainer::ShuntInfo>(m, "ShuntInfo", DocIterator::ShuntInfo.c_str())
.def_readonly("id", &ShuntContainer::ShuntInfo::id, DocIterator::id.c_str())
Expand All @@ -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_<TrafoContainer::TrafoInfo>(m, "TrafoInfo", DocIterator::TrafoInfo.c_str())
.def_readonly("id", &TrafoContainer::TrafoInfo::id, DocIterator::id.c_str())
Expand Down Expand Up @@ -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_<LineContainer::LineInfo>(m, "LineInfo", DocIterator::LineInfo.c_str())
.def_readonly("id", &LineContainer::LineInfo::id, DocIterator::id.c_str())
Expand Down Expand Up @@ -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_<DCLineContainer::DCLineInfo>(m, "DCLineInfo", DocIterator::DCLineInfo.c_str())
.def_readonly("id", &DCLineContainer::DCLineInfo::id, DocIterator::id.c_str())
Expand Down
18 changes: 15 additions & 3 deletions src/powerflow_algorithm/BaseDCAlgo.tpp
Original file line number Diff line number Diff line change
Expand Up @@ -7,6 +7,8 @@
// This file is part of LightSim2grid, LightSim2grid implements a c++ backend targeting the Grid2Op platform.

// #include "DCSolver.h"
#include <limits> // for nans
#include <cmath> // for nans

// TODO SLACK !!!
template<class LinearSolver>
Expand Down Expand Up @@ -281,10 +283,20 @@ RealMat BaseDCAlgo<LinearSolver>::get_ptdf(const Eigen::SparseMatrix<cplx_type>

template<class LinearSolver>
RealMat BaseDCAlgo<LinearSolver>::get_lodf(const Eigen::SparseMatrix<cplx_type> & 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<real_type>::quiet_NaN();
}
}
return LODF;
}

Expand Down

0 comments on commit fdf1869

Please sign in to comment.