Skip to content

Commit

Permalink
Merge pull request #390 from PowerGridModel/feature/split-physics-and…
Browse files Browse the repository at this point in the history 
…-math

Refactor solvers
  • Loading branch information
@mgovers
mgovers committed Oct 10, 2023
2 parents 712af45 + 8a8e580 commit c402a71
Show file tree
Hide file tree
Showing 7 changed files with 446 additions and 387 deletions.
75 changes: 75 additions & 0 deletions ...model_c/power_grid_model/include/power_grid_model/math_solver/common_solver_functions.hpp
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,75 @@
// SPDX-FileCopyrightText: 2022 Contributors to the Power Grid Model project <dynamic.grid.calculation@alliander.com>
//
// SPDX-License-Identifier: MPL-2.0

#pragma once
#ifndef POWER_GRID_MODEL_COMMON_SOLVER_FUNCTIONS_HPP
#define POWER_GRID_MODEL_COMMON_SOLVER_FUNCTIONS_HPP

#include "y_bus.hpp"

#include "../calculation_parameters.hpp"

namespace power_grid_model::common_solver_functions {

template <bool sym>
void add_sources(IdxVector const& source_bus_indptr, Idx const& bus_number, YBus<sym> const& y_bus,
ComplexVector const& u_source_vector, ComplexTensor<sym>& diagonal_element, ComplexValue<sym>& u_bus) {
for (Idx source_number = source_bus_indptr[bus_number]; source_number != source_bus_indptr[bus_number + 1];
++source_number) {
ComplexTensor<sym> const y_source = y_bus.math_model_param().source_param[source_number];
diagonal_element += y_source; // add y_source to the diagonal of Ybus
u_bus += dot(y_source, ComplexValue<sym>{u_source_vector[source_number]}); // rhs += Y_source * U_source
}
}

template <bool sym> void copy_y_bus(YBus<sym> const& y_bus, ComplexTensorVector<sym>& mat_data) {
ComplexTensorVector<sym> const& ydata = y_bus.admittance();
std::transform(y_bus.map_lu_y_bus().cbegin(), y_bus.map_lu_y_bus().cend(), mat_data.begin(), [&](Idx k) {
if (k == -1) {
return ComplexTensor<sym>{};
}
return ydata[k];
});
}

template <bool sym>
void calculate_source_result(Idx const& bus_number, YBus<sym> const& y_bus, PowerFlowInput<sym> const& input,
MathOutput<sym>& output, IdxVector const& source_bus_indptr) {
for (Idx source = source_bus_indptr[bus_number]; source != source_bus_indptr[bus_number + 1]; ++source) {
ComplexValue<sym> const u_ref{input.source[source]};
ComplexTensor<sym> const y_ref = y_bus.math_model_param().source_param[source];
output.source[source].i = dot(y_ref, u_ref - output.u[bus_number]);
output.source[source].s = output.u[bus_number] * conj(output.source[source].i);
}
}

template <bool sym, class LoadGenFunc>
void calculate_load_gen_result(Idx const& bus_number, PowerFlowInput<sym> const& input, MathOutput<sym>& output,
IdxVector const& load_gen_bus_indptr, LoadGenFunc const& load_gen_func) {
for (Idx load_gen = load_gen_bus_indptr[bus_number]; load_gen != load_gen_bus_indptr[bus_number + 1]; ++load_gen) {
switch (LoadGenType const type = load_gen_func(load_gen); type) {
using enum LoadGenType;

case const_pq:
// always same power
output.load_gen[load_gen].s = input.s_injection[load_gen];
break;
case const_y:
// power is quadratic relation to voltage
output.load_gen[load_gen].s = input.s_injection[load_gen] * abs2(output.u[bus_number]);
break;
case const_i:
// power is linear relation to voltage
output.load_gen[load_gen].s = input.s_injection[load_gen] * cabs(output.u[bus_number]);
break;
default:
throw MissingCaseForEnumError("Power injection", type);
}
output.load_gen[load_gen].i = conj(output.load_gen[load_gen].s / output.u[bus_number]);
}
}

} // namespace power_grid_model::common_solver_functions

#endif
80 changes: 41 additions & 39 deletions ...l_c/power_grid_model/include/power_grid_model/math_solver/iterative_current_pf_solver.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -55,6 +55,7 @@ Iterative Power Flow
*/

#include "block_matrix.hpp"
#include "common_solver_functions.hpp"
#include "iterative_pf_solver.hpp"
#include "sparse_lu_solver.hpp"
#include "y_bus.hpp"
Expand Down Expand Up @@ -85,19 +86,12 @@ template <bool sym> class IterativeCurrentPFSolver : public IterativePFSolver<sy
// Add source admittance to Y bus and set variable for prepared y bus to true
void initialize_derived_solver(YBus<sym> const& y_bus, MathOutput<sym> const& /* output */) {
IdxVector const& source_bus_indptr = *this->source_bus_indptr_;
ComplexTensorVector<sym> const& ydata = y_bus.admittance();
IdxVector const& bus_entry = y_bus.lu_diag();
// if Y bus is not up to date
// re-build matrix and prefactorize Build y bus data with source admittance
if (y_data_ptr_ != &y_bus.admittance()) {
ComplexTensorVector<sym> mat_data(y_bus.nnz_lu());
// copy y bus data
std::transform(y_bus.map_lu_y_bus().cbegin(), y_bus.map_lu_y_bus().cend(), mat_data.begin(), [&](Idx k) {
if (k == -1) {
return ComplexTensor<sym>{};
}
return ydata[k];
});
common_solver_functions::copy_y_bus<sym>(y_bus, mat_data);

// loop bus
for (Idx bus_number = 0; bus_number != this->n_bus_; ++bus_number) {
Expand Down Expand Up @@ -132,37 +126,8 @@ template <bool sym> class IterativeCurrentPFSolver : public IterativePFSolver<sy

// loop buses: i
for (Idx bus_number = 0; bus_number != this->n_bus_; ++bus_number) {
// loop loads/generation: j
for (Idx load_number = load_gen_bus_indptr[bus_number]; load_number != load_gen_bus_indptr[bus_number + 1];
++load_number) {
// load type
LoadGenType const type = load_gen_type[load_number];
switch (type) {
using enum LoadGenType;

case const_pq:
// I_inj_i = conj(S_inj_j/U_i) for constant PQ type
rhs_u_[bus_number] += conj(input.s_injection[load_number] / u[bus_number]);
break;
case const_y:
// I_inj_i = conj((S_inj_j * abs(U_i)^2) / U_i) = conj((S_inj_j) * U_i for const impedance type
rhs_u_[bus_number] += conj(input.s_injection[load_number]) * u[bus_number];
break;
case const_i:
// I_inj_i = conj(S_inj_j*abs(U_i)/U_i) for const current type
rhs_u_[bus_number] += conj(input.s_injection[load_number] * cabs(u[bus_number]) / u[bus_number]);
break;
default:
throw MissingCaseForEnumError("Injection current calculation", type);
}
}
// loop sources: j
for (Idx source_number = source_bus_indptr[bus_number]; source_number != source_bus_indptr[bus_number + 1];
++source_number) {
// I_inj_i += Y_source_j * U_ref_j
rhs_u_[bus_number] += dot(y_bus.math_model_param().source_param[source_number],
ComplexValue<sym>{input.source[source_number]});
}
add_loads(bus_number, input, load_gen_bus_indptr, load_gen_type, u);
add_sources(bus_number, y_bus, input, source_bus_indptr);
}
}

Expand Down Expand Up @@ -192,6 +157,43 @@ template <bool sym> class IterativeCurrentPFSolver : public IterativePFSolver<sy
// sparse solver
SparseLUSolver<ComplexTensor<sym>, ComplexValue<sym>, ComplexValue<sym>> sparse_solver_;
std::shared_ptr<BlockPermArray const> perm_;

void add_loads(Idx const& bus_number, PowerFlowInput<sym> const& input, IdxVector const& load_gen_bus_indptr,
std::vector<LoadGenType> const& load_gen_type, ComplexValueVector<sym> const& u) {
for (Idx load_number = load_gen_bus_indptr[bus_number]; load_number != load_gen_bus_indptr[bus_number + 1];
++load_number) {
// load type
LoadGenType const type = load_gen_type[load_number];
switch (type) {
using enum LoadGenType;

case const_pq:
// I_inj_i = conj(S_inj_j/U_i) for constant PQ type
rhs_u_[bus_number] += conj(input.s_injection[load_number] / u[bus_number]);
break;
case const_y:
// I_inj_i = conj((S_inj_j * abs(U_i)^2) / U_i) = conj((S_inj_j) * U_i for const impedance type
rhs_u_[bus_number] += conj(input.s_injection[load_number]) * u[bus_number];
break;
case const_i:
// I_inj_i = conj(S_inj_j*abs(U_i)/U_i) for const current type
rhs_u_[bus_number] += conj(input.s_injection[load_number] * cabs(u[bus_number]) / u[bus_number]);
break;
default:
throw MissingCaseForEnumError("Injection current calculation", type);
}
}
}

void add_sources(Idx const& bus_number, YBus<sym> const& y_bus, PowerFlowInput<sym> const& input,
IdxVector const& source_bus_indptr) {
for (Idx source_number = source_bus_indptr[bus_number]; source_number != source_bus_indptr[bus_number + 1];
++source_number) {
// I_inj_i += Y_source_j * U_ref_j
rhs_u_[bus_number] += dot(y_bus.math_model_param().source_param[source_number],
ComplexValue<sym>{input.source[source_number]});
}
}
};

template class IterativeCurrentPFSolver<true>;
Expand Down
41 changes: 6 additions & 35 deletions ...rid_model_c/power_grid_model/include/power_grid_model/math_solver/iterative_pf_solver.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,7 @@
*/

// Check if all includes needed
#include "common_solver_functions.hpp"
#include "y_bus.hpp"

#include "../calculation_parameters.hpp"
Expand Down Expand Up @@ -112,41 +113,11 @@ template <bool sym, typename DerivedSolver> class IterativePFSolver {
output.load_gen.resize(load_gen_bus_indptr_->back());
output.bus_injection.resize(n_bus_);

// loop all bus
for (Idx bus = 0; bus != n_bus_; ++bus) {
// source
for (Idx source = (*source_bus_indptr_)[bus]; source != (*source_bus_indptr_)[bus + 1]; ++source) {
ComplexValue<sym> const u_ref{input.source[source]};
ComplexTensor<sym> const y_ref = y_bus.math_model_param().source_param[source];
output.source[source].i = dot(y_ref, u_ref - output.u[bus]);
output.source[source].s = output.u[bus] * conj(output.source[source].i);
}

// load_gen
for (Idx load_gen = (*load_gen_bus_indptr_)[bus]; load_gen != (*load_gen_bus_indptr_)[bus + 1];
++load_gen) {
LoadGenType const type = (*load_gen_type_)[load_gen];
switch (type) {
using enum LoadGenType;

case const_pq:
// always same power
output.load_gen[load_gen].s = input.s_injection[load_gen];
break;
case const_y:
// power is quadratic relation to voltage
output.load_gen[load_gen].s =
input.s_injection[load_gen] * cabs(output.u[bus]) * cabs(output.u[bus]);
break;
case const_i:
// power is linear relation to voltage
output.load_gen[load_gen].s = input.s_injection[load_gen] * cabs(output.u[bus]);
break;
default:
throw MissingCaseForEnumError("Power injection", type);
}
output.load_gen[load_gen].i = conj(output.load_gen[load_gen].s / output.u[bus]);
}
for (Idx bus_number = 0; bus_number != n_bus_; ++bus_number) {
common_solver_functions::calculate_source_result<sym>(bus_number, y_bus, input, output,
*source_bus_indptr_);
common_solver_functions::calculate_load_gen_result<sym>(bus_number, input, output, *load_gen_bus_indptr_,
[this](Idx i) { return (*load_gen_type_)[i]; });
}
output.bus_injection = y_bus.calculate_injection(output.u);
}
Expand Down
82 changes: 32 additions & 50 deletions ...r_grid_model_c/power_grid_model/include/power_grid_model/math_solver/linear_pf_solver.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,6 +30,7 @@ if there are sources
*/

#include "common_solver_functions.hpp"
#include "sparse_lu_solver.hpp"
#include "y_bus.hpp"

Expand Down Expand Up @@ -57,9 +58,6 @@ template <bool sym> class LinearPFSolver {

MathOutput<sym> run_power_flow(YBus<sym> const& y_bus, PowerFlowInput<sym> const& input,
CalculationInfo& calculation_info) {
// getter
ComplexTensorVector<sym> const& ydata = y_bus.admittance();
IdxVector const& bus_entry = y_bus.lu_diag();
// output
MathOutput<sym> output;
output.u.resize(n_bus_);
Expand All @@ -68,36 +66,8 @@ template <bool sym> class LinearPFSolver {

// prepare matrix
Timer sub_timer(calculation_info, 2221, "Prepare matrix");

// copy y bus data
std::transform(y_bus.map_lu_y_bus().cbegin(), y_bus.map_lu_y_bus().cend(), mat_data_.begin(), [&](Idx k) {
if (k == -1) {
return ComplexTensor<sym>{};
}
return ydata[k];
});

// loop to all loads and sources, j as load number
IdxVector const& load_gen_bus_idxptr = *load_gen_bus_indptr_;
IdxVector const& source_bus_indptr = *source_bus_indptr_;
for (Idx bus_number = 0; bus_number != n_bus_; ++bus_number) {
Idx const data_sequence = bus_entry[bus_number];
// loop loads
for (Idx load_number = load_gen_bus_idxptr[bus_number]; load_number != load_gen_bus_idxptr[bus_number + 1];
++load_number) {
// YBus_diag += -conj(S_base)
add_diag(mat_data_[data_sequence], -conj(input.s_injection[load_number]));
}
// loop sources
for (Idx source_number = source_bus_indptr[bus_number]; source_number != source_bus_indptr[bus_number + 1];
++source_number) {
// YBus_diag += Y_source
mat_data_[data_sequence] += y_bus.math_model_param().source_param[source_number];
// rhs += Y_source_j * U_ref_j
output.u[bus_number] += dot(y_bus.math_model_param().source_param[source_number],
ComplexValue<sym>{input.source[source_number]});
}
}
common_solver_functions::copy_y_bus<sym>(y_bus, mat_data_);
prepare_matrix_and_rhs(y_bus, input, output);

// solve
// u vector will have I_injection for slack bus for now
Expand All @@ -123,6 +93,30 @@ template <bool sym> class LinearPFSolver {
SparseLUSolver<ComplexTensor<sym>, ComplexValue<sym>, ComplexValue<sym>> sparse_solver_;
typename SparseLUSolver<ComplexTensor<sym>, ComplexValue<sym>, ComplexValue<sym>>::BlockPermArray perm_;

void prepare_matrix_and_rhs(YBus<sym> const& y_bus, PowerFlowInput<sym> const& input, MathOutput<sym>& output) {
// getter
IdxVector const& bus_entry = y_bus.lu_diag();
IdxVector const& load_gen_bus_idxptr = *load_gen_bus_indptr_;
IdxVector const& source_bus_indptr = *source_bus_indptr_;
for (Idx bus_number = 0; bus_number != n_bus_; ++bus_number) {
Idx const diagonal_position = bus_entry[bus_number];
auto& diagonal_element = mat_data_[diagonal_position];
auto& u_bus = output.u[bus_number];
add_loads(load_gen_bus_idxptr, bus_number, input, diagonal_element);
common_solver_functions::add_sources<sym>(source_bus_indptr, bus_number, y_bus, input.source,
diagonal_element, u_bus);
}
}

static void add_loads(IdxVector const& load_gen_bus_idxptr, Idx const& bus_number, PowerFlowInput<sym> const& input,
ComplexTensor<sym>& diagonal_element) {
for (Idx load_number = load_gen_bus_idxptr[bus_number]; load_number != load_gen_bus_idxptr[bus_number + 1];
++load_number) {
// YBus_diag += -conj(S_base)
add_diag(diagonal_element, -conj(input.s_injection[load_number]));
}
}

void calculate_result(YBus<sym> const& y_bus, PowerFlowInput<sym> const& input, MathOutput<sym>& output) {
// call y bus
output.branch = y_bus.template calculate_branch_flow<BranchMathOutput<sym>>(output.u);
Expand All @@ -133,23 +127,11 @@ template <bool sym> class LinearPFSolver {
output.load_gen.resize(load_gen_bus_indptr_->back());
output.bus_injection.resize(n_bus_);

// loop all bus
for (Idx bus = 0; bus != n_bus_; ++bus) {
// source
for (Idx source = (*source_bus_indptr_)[bus]; source != (*source_bus_indptr_)[bus + 1]; ++source) {
ComplexValue<sym> const u_ref{input.source[source]};
ComplexTensor<sym> const y_ref = y_bus.math_model_param().source_param[source];
output.source[source].i = dot(y_ref, u_ref - output.u[bus]);
output.source[source].s = output.u[bus] * conj(output.source[source].i);
}

// load_gen
for (Idx load_gen = (*load_gen_bus_indptr_)[bus]; load_gen != (*load_gen_bus_indptr_)[bus + 1];
++load_gen) {
// power is always quadratic relation to voltage for linear pf
output.load_gen[load_gen].s = input.s_injection[load_gen] * abs2(output.u[bus]);
output.load_gen[load_gen].i = conj(output.load_gen[load_gen].s / output.u[bus]);
}
for (Idx bus_number = 0; bus_number != n_bus_; ++bus_number) {
common_solver_functions::calculate_source_result<sym>(bus_number, y_bus, input, output,
*source_bus_indptr_);
common_solver_functions::calculate_load_gen_result<sym>(bus_number, input, output, *load_gen_bus_indptr_,
[](Idx /*i*/) { return LoadGenType::const_y; });
}
output.bus_injection = y_bus.calculate_injection(output.u);
}
Expand Down

0 comments on commit c402a71

Please sign in to comment.