diff --git a/lightsim2grid/tests/test_ptdf.py b/lightsim2grid/tests/test_ptdf.py new file mode 100644 index 0000000..92ad48c --- /dev/null +++ b/lightsim2grid/tests/test_ptdf.py @@ -0,0 +1,106 @@ +# 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 warnings +from scipy.sparse.linalg import spsolve + +from lightsim2grid.gridmodel import init +from lightsim2grid.solver import SolverType + +import pdb + +class TestCase14SLU(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() + self.gridmodel.change_solver(self.get_solver_type()) + self.gridmodel.dc_pf(self.V_init, 1, 1e-8) + 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_pp(self): + # test the right computation of Bf matrix (PTDF derived from python, might be slower) + Bf = 1.0 * self.gridmodel.get_Bf() + PTDF = np.zeros((self.nbr, self.nb)) + PTDF[:, self.noslack] = spsolve(self.Bbus[np.ix_(self.noslack, self.noref)].T, Bf[:, self.noref].toarray().T).T + # test the solver works correctly + tmp_mat = self.Bbus[np.ix_(self.noslack, self.noref)].T.todense() + for line_id in range(self.nbr): + solve_error = np.abs(np.dot(tmp_mat, PTDF[line_id, self.noslack]).T - Bf[line_id, self.noref].toarray().T).max() + assert solve_error <= self.tol, f"error for line {line_id}: {solve_error:.2e}MW" + # test powerflow are correct + res_ptdf = np.dot(PTDF, self.dcSbus * self.gridmodel.get_sn_mva()) + assert np.abs(res_ptdf - self.res_powerflow).max() <= self.tol, f"max error for powerflow: {np.abs(res_ptdf - self.res_powerflow).max():.2e}MW" + + def test_ptdf_from_ls(self): + # now test the right computation of the PTDF + Bf = 1.0 * self.gridmodel.get_Bf() + PTDF2 = 1.0 * self.gridmodel.get_ptdf() + # test the solver works correctly + tmp_mat = self.Bbus[np.ix_(self.noslack, self.noref)].T.todense() + for line_id in range(self.nbr): + solve_error = np.abs(np.dot(tmp_mat, PTDF2[line_id, self.noslack]).T - Bf[line_id, self.noref].toarray().T).max() + assert solve_error <= self.tol, f"error for line {line_id}: {solve_error:.2e}MW" + # test powerflow are correct + res_ptdf2 = np.dot(PTDF2, self.dcSbus * self.gridmodel.get_sn_mva()) + np.where(np.abs(res_ptdf2 - self.res_powerflow) >= 1e3) + assert np.abs(res_ptdf2 - self.res_powerflow).max() <= self.tol, f"max error for powerflow: {np.abs(res_ptdf2 - self.res_powerflow).max():.2e}MW" + + +class TestCase30SLU(TestCase14SLU): + def make_grid(self): + res = pn.case30() + return res + + +class TestCase118SLU(TestCase14SLU): + def make_grid(self): + res = pn.case118() + return res + + +class TestCase14KLU(TestCase14SLU): + def get_solver_type(self): + return SolverType.KLUDC + + +class TestCase30KLU(TestCase30SLU): + def get_solver_type(self): + return SolverType.KLUDC + + +class TestCase118KLU(TestCase118SLU): + def get_solver_type(self): + return SolverType.KLUDC + + +if __name__ == "__main__": + unittest.main() diff --git a/src/DCSolver.h b/src/DCSolver.h index 38c110d..eeb14a2 100644 --- a/src/DCSolver.h +++ b/src/DCSolver.h @@ -15,7 +15,12 @@ template class BaseDCSolver: public BaseSolver { public: - BaseDCSolver():BaseSolver(false), _linear_solver(), need_factorize_(true), nb_dcbus_solver_(0){}; + BaseDCSolver(): + BaseSolver(false), + _linear_solver(), + need_factorize_(true), + sizeYbus_with_slack_(0), + sizeYbus_without_slack_(0){}; ~BaseDCSolver(){} @@ -56,7 +61,8 @@ class BaseDCSolver: public BaseSolver bool need_factorize_; // save this not to recompute them when not needed - int nb_dcbus_solver_; + int sizeYbus_with_slack_; + int sizeYbus_without_slack_; RealVect dcSbus_noslack_; Eigen::SparseMatrix dcYbus_noslack_; Eigen::VectorXi my_pv_; diff --git a/src/DCSolver.tpp b/src/DCSolver.tpp index 1e65858..fee67b2 100644 --- a/src/DCSolver.tpp +++ b/src/DCSolver.tpp @@ -28,7 +28,7 @@ bool BaseDCSolver::compute_pf(const Eigen::SparseMatrix auto timer = CustTimer(); BaseSolver::reset_timer(); - nb_dcbus_solver_ = static_cast(Ybus.rows()); + sizeYbus_with_slack_ = static_cast(Ybus.rows()); #ifdef __COUT_TIMES auto timer_preproc = CustTimer(); @@ -43,14 +43,15 @@ bool BaseDCSolver::compute_pf(const Eigen::SparseMatrix // const Eigen::VectorXi & my_pv = pv; // find the slack buses - slack_buses_ids_solver_ = extract_slack_bus_id(my_pv_, pq, nb_dcbus_solver_); + slack_buses_ids_solver_ = extract_slack_bus_id(my_pv_, pq, sizeYbus_with_slack_); + sizeYbus_without_slack_ = sizeYbus_with_slack_ - slack_buses_ids_solver_.size(); // corresp bus -> solverbus - fill_mat_bus_id(nb_dcbus_solver_); + fill_mat_bus_id(sizeYbus_with_slack_); // remove the slack bus from Ybus // and extract only real part - fill_dcYbus_noslack(nb_dcbus_solver_, Ybus); + fill_dcYbus_noslack(sizeYbus_with_slack_, Ybus); #ifdef __COUT_TIMES std::cout << "\t dc: preproc: " << 1000. * timer_preproc.duration() << "ms" << std::endl; @@ -72,8 +73,8 @@ bool BaseDCSolver::compute_pf(const Eigen::SparseMatrix } // remove the slack bus from Sbus - dcSbus_noslack_ = RealVect::Constant(nb_dcbus_solver_ - slack_buses_ids_solver_.size(), my_zero_); - for (int k=0; k < nb_dcbus_solver_; ++k){ + dcSbus_noslack_ = RealVect::Constant(sizeYbus_without_slack_, my_zero_); + for (int k=0; k < sizeYbus_with_slack_; ++k){ if(mat_bus_id_(k) == -1) continue; // I don't add anything to the slack bus const int col_res = mat_bus_id_(k); dcSbus_noslack_(col_res) = std::real(Sbus_tmp(k)); @@ -108,9 +109,9 @@ bool BaseDCSolver::compute_pf(const Eigen::SparseMatrix // retrieve back the results in the proper shape (add back the slack bus) // TODO have a better way for this, for example using `.segment(0,npv)` // see the BaseSolver.cpp: _evaluate_Fx - RealVect Va_dc = RealVect::Constant(nb_dcbus_solver_, my_zero_); + RealVect Va_dc = RealVect::Constant(sizeYbus_with_slack_, my_zero_); // fill Va from dc approx - for (int ybus_id=0; ybus_id < nb_dcbus_solver_; ++ybus_id){ + for (int ybus_id=0; ybus_id < sizeYbus_with_slack_; ++ybus_id){ if(mat_bus_id_(ybus_id) == -1) continue; // slack bus is handled elsewhere const int bus_me = mat_bus_id_(ybus_id); Va_dc(ybus_id) = Va_dc_without_slack(bus_me); @@ -161,7 +162,7 @@ void BaseDCSolver::fill_dcYbus_noslack(int nb_bus_solver, const Ei template template // ref_mat_type should be `real_type` or `cplx_type` void BaseDCSolver::remove_slack_buses(int nb_bus_solver, const Eigen::SparseMatrix & ref_mat, Eigen::SparseMatrix & res_mat){ - res_mat = Eigen::SparseMatrix(nb_bus_solver - slack_buses_ids_solver_.size(), nb_bus_solver - slack_buses_ids_solver_.size()); // TODO dist slack: -1 or -mat_bus_id_.size() here ???? + res_mat = Eigen::SparseMatrix(sizeYbus_without_slack_, sizeYbus_without_slack_); // TODO dist slack: -1 or -mat_bus_id_.size() here ???? std::vector > tripletList; tripletList.reserve(ref_mat.nonZeros()); for (int k=0; k < nb_bus_solver; ++k){ @@ -185,7 +186,8 @@ void BaseDCSolver::reset(){ BaseSolver::reset(); _linear_solver.reset(); need_factorize_ = true; - nb_dcbus_solver_ = 0; + sizeYbus_with_slack_ = 0; + sizeYbus_without_slack_ = 0; dcSbus_noslack_ = RealVect(); dcYbus_noslack_ = Eigen::SparseMatrix(); my_pv_ = Eigen::VectorXi(); @@ -197,7 +199,7 @@ template RealMat BaseDCSolver::get_ptdf(const Eigen::SparseMatrix & dcYbus){ Eigen::SparseMatrix Bf_T_with_slack; RealMat PTDF; - RealVect rhs; + RealVect rhs = RealVect::Zero(sizeYbus_without_slack_); // TODO dist slack: -1 or -mat_bus_id_.size() here ???? // TODO PTDF: sparse matrix ? // TODO PTDF: distributed slack // TODO PTDF: check that the solver has converged @@ -218,8 +220,7 @@ RealMat BaseDCSolver::get_ptdf(const Eigen::SparseMatrix::InnerIterator it(Bf_T_with_slack, line_id); it; ++it) @@ -229,38 +230,15 @@ RealMat BaseDCSolver::get_ptdf(const Eigen::SparseMatrix