Skip to content

Commit

Permalink
implement and test computation of ptdf
Browse files Browse the repository at this point in the history
  • Loading branch information
BDonnot committed Dec 7, 2023
1 parent 4602b15 commit 14d62bd
Show file tree
Hide file tree
Showing 3 changed files with 129 additions and 39 deletions.
106 changes: 106 additions & 0 deletions lightsim2grid/tests/test_ptdf.py
Original file line number Diff line number Diff line change
@@ -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()
10 changes: 8 additions & 2 deletions src/DCSolver.h
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,12 @@ template<class LinearSolver>
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(){}

Expand Down Expand Up @@ -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<real_type> dcYbus_noslack_;
Eigen::VectorXi my_pv_;
Expand Down
52 changes: 15 additions & 37 deletions src/DCSolver.tpp
Original file line number Diff line number Diff line change
Expand Up @@ -28,7 +28,7 @@ bool BaseDCSolver<LinearSolver>::compute_pf(const Eigen::SparseMatrix<cplx_type>

auto timer = CustTimer();
BaseSolver::reset_timer();
nb_dcbus_solver_ = static_cast<int>(Ybus.rows());
sizeYbus_with_slack_ = static_cast<int>(Ybus.rows());

#ifdef __COUT_TIMES
auto timer_preproc = CustTimer();
Expand All @@ -43,14 +43,15 @@ bool BaseDCSolver<LinearSolver>::compute_pf(const Eigen::SparseMatrix<cplx_type>
// 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;
Expand All @@ -72,8 +73,8 @@ bool BaseDCSolver<LinearSolver>::compute_pf(const Eigen::SparseMatrix<cplx_type>
}

// 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));
Expand Down Expand Up @@ -108,9 +109,9 @@ bool BaseDCSolver<LinearSolver>::compute_pf(const Eigen::SparseMatrix<cplx_type>
// 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);
Expand Down Expand Up @@ -161,7 +162,7 @@ void BaseDCSolver<LinearSolver>::fill_dcYbus_noslack(int nb_bus_solver, const Ei
template<class LinearSolver>
template<typename ref_mat_type> // ref_mat_type should be `real_type` or `cplx_type`
void BaseDCSolver<LinearSolver>::remove_slack_buses(int nb_bus_solver, const Eigen::SparseMatrix<ref_mat_type> & ref_mat, Eigen::SparseMatrix<real_type> & res_mat){
res_mat = Eigen::SparseMatrix<real_type>(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<real_type>(sizeYbus_without_slack_, sizeYbus_without_slack_); // TODO dist slack: -1 or -mat_bus_id_.size() here ????
std::vector<Eigen::Triplet<real_type> > tripletList;
tripletList.reserve(ref_mat.nonZeros());
for (int k=0; k < nb_bus_solver; ++k){
Expand All @@ -185,7 +186,8 @@ void BaseDCSolver<LinearSolver>::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<real_type>();
my_pv_ = Eigen::VectorXi();
Expand All @@ -197,7 +199,7 @@ template<class LinearSolver>
RealMat BaseDCSolver<LinearSolver>::get_ptdf(const Eigen::SparseMatrix<cplx_type> & dcYbus){
Eigen::SparseMatrix<real_type> 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
Expand All @@ -218,8 +220,7 @@ RealMat BaseDCSolver<LinearSolver>::get_ptdf(const Eigen::SparseMatrix<cplx_type
const Eigen::VectorXi ind_no_slack = Eigen::VectorXi::Map(&ind_no_slack_[0], ind_no_slack_.size());

// solve iteratively the linear systems (one per powerline)
PTDF = RealMat(Bf_T_with_slack.cols(), Bf_T_with_slack.rows()); // rows and cols are "inverted" because the matrix Bf is transposed
rhs = RealVect::Zero(Bf_T_with_slack.cols() - slack_buses_ids_solver_.size()); // TODO dist slack: -1 or -mat_bus_id_.size() here ????
PTDF = RealMat::Zero(Bf_T_with_slack.cols(), Bf_T_with_slack.rows()); // rows and cols are "inverted" because the matrix Bf is transposed
for (int line_id=0; line_id < nb_pow_tr; ++line_id){
// build the rhs vector
for (typename Eigen::SparseMatrix<real_type>::InnerIterator it(Bf_T_with_slack, line_id); it; ++it)
Expand All @@ -229,38 +230,15 @@ RealMat BaseDCSolver<LinearSolver>::get_ptdf(const Eigen::SparseMatrix<cplx_type
const auto col_res = mat_bus_id_(bus_id);
rhs[col_res] = it.value();
}
if (line_id == 16){
std::cout << "line 16\n";
for(auto i = 0; i < nb_bus; ++i) std::cout << rhs[i] << ", ";
std::cout << std::endl;
}
if (line_id == 17){
std::cout << "line 17\n";
for(auto i = 0; i < nb_bus; ++i) std::cout << rhs[i] << ", ";
std::cout << std::endl;
}
if (line_id == 18){
std::cout << "line 18\n";
for(auto i = 0; i < nb_bus; ++i) std::cout << rhs[i] << ", ";
std::cout << std::endl;
}

// solve the linear system
_linear_solver.solve(dcYbus_noslack_, rhs, true); // I don't need to refactorize the matrix (hence the `true`)

if (line_id == 18){
std::cout << "res for 18\n";
RealVect res = dcYbus_noslack_ * rhs;
for(auto i = 0; i < nb_bus; ++i) std::cout << res[i] << ", ";
std::cout << std::endl;
}

// assign results to the PTDF matrix
PTDF(line_id, ind_no_slack) = rhs;

// reset the rhs vector to 0.
rhs.array() = 0.;
// rhs = RealVect::Zero(Bf_T_with_slack.cols() - slack_buses_ids_solver_.size());
// rhs = RealVect::Zero(sizeYbus_without_slack_);
}
// TODO PTDF: if the solver can solve the directly, do that instead
return PTDF;
Expand Down

0 comments on commit 14d62bd

Please sign in to comment.