Skip to content
Original file line number Diff line number Diff line change
Expand Up @@ -80,17 +80,17 @@ class IterativeCurrentPFSolver : public IterativePFSolver<sym_type, IterativeCur

static constexpr auto is_iterative = true;

IterativeCurrentPFSolver(YBus<sym> const& y_bus, std::shared_ptr<MathModelTopology const> const& topo_ptr)
IterativeCurrentPFSolver(YBus<sym> const& y_bus, MathModelTopology const& topo_ptr)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this is no longer a pointer

Suggested change
IterativeCurrentPFSolver(YBus<sym> const& y_bus, MathModelTopology const& topo_ptr)
IterativeCurrentPFSolver(YBus<sym> const& y_bus, MathModelTopology const& topo)

: IterativePFSolver<sym, IterativeCurrentPFSolver>{y_bus, topo_ptr},
rhs_u_(y_bus.size()),
sparse_solver_{y_bus.shared_indptr_lu(), y_bus.shared_indices_lu(), y_bus.shared_diag_lu()} {}
sparse_solver_{y_bus.row_indptr_lu(), y_bus.col_indices_lu(), y_bus.lu_diag()} {}

// Add source admittance to Y bus and set variable for prepared y bus to true
void initialize_derived_solver(YBus<sym> const& y_bus, PowerFlowInput<sym> const& input,
SolverOutput<sym>& output) {
make_flat_start(input, output.u);

auto const& sources_per_bus = *this->sources_per_bus_;
auto const& sources_per_bus = this->sources_per_bus_.get();
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
Expand Down Expand Up @@ -119,14 +119,14 @@ class IterativeCurrentPFSolver : public IterativePFSolver<sym_type, IterativeCur
// Prepare matrix calculates injected current, i.e., RHS of solver for each iteration.
void prepare_matrix_and_rhs(YBus<sym> const& y_bus, PowerFlowInput<sym> const& input,
ComplexValueVector<sym> const& u) {
std::vector<LoadGenType> const& load_gen_type = *this->load_gen_type_;
std::vector<LoadGenType> const& load_gen_type = this->load_gen_type_.get();

// set rhs to zero for iteration start
std::fill(rhs_u_.begin(), rhs_u_.end(), ComplexValue<sym>{0.0});

// loop buses: i
for (auto const& [bus_number, load_gens, sources] :
enumerated_zip_sequence(*this->load_gens_per_bus_, *this->sources_per_bus_)) {
enumerated_zip_sequence(this->load_gens_per_bus_.get(), this->sources_per_bus_.get())) {
add_loads(load_gens, bus_number, input, load_gen_type, u);
add_sources(sources, bus_number, y_bus, input);
}
Expand Down Expand Up @@ -197,11 +197,11 @@ class IterativeCurrentPFSolver : public IterativePFSolver<sym_type, IterativeCur
}

void make_flat_start(PowerFlowInput<sym> const& input, ComplexValueVector<sym>& output_u) {
std::vector<double> const& phase_shift = *this->phase_shift_;
std::vector<double> const& phase_shift = this->phase_shift_.get();
// average u_ref of all sources
DoubleComplex const u_ref = [&]() {
DoubleComplex sum_u_ref = 0.0;
for (auto const& [bus, sources] : enumerated_zip_sequence(*this->sources_per_bus_)) {
for (auto const& [bus, sources] : enumerated_zip_sequence(this->sources_per_bus_.get())) {
for (Idx const source : sources) {
sum_u_ref += input.source[source] * std::exp(1.0i * -phase_shift[bus]); // offset phase shift
}
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -75,12 +75,12 @@ template <symmetry_tag sym_type> class IterativeLinearSESolver {
static constexpr Idx bsr_block_size_ = is_symmetric_v<sym> ? 2 : 6;

public:
IterativeLinearSESolver(YBus<sym> const& y_bus, std::shared_ptr<MathModelTopology const> topo_ptr)
IterativeLinearSESolver(YBus<sym> const& y_bus, MathModelTopology const& topo_ptr)
: n_bus_{y_bus.size()},
math_topo_{std::move(topo_ptr)},
math_topo_{topo_ptr},
data_gain_(y_bus.nnz_lu()),
x_rhs_(y_bus.size()),
sparse_solver_{y_bus.shared_indptr_lu(), y_bus.shared_indices_lu(), y_bus.shared_diag_lu()},
sparse_solver_{y_bus.row_indptr_lu(), y_bus.col_indices_lu(), y_bus.lu_diag()},
perm_(y_bus.size()) {}

SolverOutput<sym> run_state_estimation(YBus<sym> const& y_bus, StateEstimationInput<sym> const& input,
Expand Down Expand Up @@ -111,7 +111,7 @@ template <symmetry_tag sym_type> class IterativeLinearSESolver {
sub_timer = Timer{log, LogEvent::initialize_voltages}; // TODO(mgovers): make scoped subtimers
RealValue<sym> const mean_angle_shift = measured_values.mean_angle_shift();
for (Idx bus = 0; bus != n_bus_; ++bus) {
output.u[bus] = exp(1.0i * (mean_angle_shift + math_topo_->phase_shift[bus]));
output.u[bus] = exp(1.0i * (mean_angle_shift + math_topo_.get().phase_shift[bus]));
}

// loop to iterate
Expand Down Expand Up @@ -155,7 +155,7 @@ template <symmetry_tag sym_type> class IterativeLinearSESolver {

Idx n_bus_;
// shared topo data
std::shared_ptr<MathModelTopology const> math_topo_;
std::reference_wrapper<MathModelTopology const> math_topo_;

// data for gain matrix
std::vector<ILSEGainBlock<sym>> data_gain_;
Expand Down Expand Up @@ -359,7 +359,7 @@ template <symmetry_tag sym_type> class IterativeLinearSESolver {
if (has_angle) {
return 1.0;
}
auto const& voltage = x_rhs_[math_topo_->slack_bus].u();
auto const& voltage = x_rhs_[math_topo_.get().slack_bus].u();
auto const& voltage_a = [&voltage]() -> auto const& {
if constexpr (is_symmetric_v<sym>) {
return voltage;
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -81,22 +81,22 @@ template <symmetry_tag sym, typename DerivedSolver> class IterativePFSolver {
}

void calculate_result(YBus<sym> const& y_bus, PowerFlowInput<sym> const& input, SolverOutput<sym>& output) {
detail::calculate_pf_result(y_bus, input, *sources_per_bus_, *load_gens_per_bus_, output,
[this](Idx i) { return (*load_gen_type_)[i]; });
detail::calculate_pf_result(y_bus, input, sources_per_bus_.get(), load_gens_per_bus_.get(), output,
[this](Idx i) { return (load_gen_type_.get())[i]; });
}

private:
Idx n_bus_;
std::shared_ptr<DoubleVector const> phase_shift_;
std::shared_ptr<SparseGroupedIdxVector const> load_gens_per_bus_;
std::shared_ptr<DenseGroupedIdxVector const> sources_per_bus_;
std::shared_ptr<std::vector<LoadGenType> const> load_gen_type_;
IterativePFSolver(YBus<sym> const& y_bus, std::shared_ptr<MathModelTopology const> const& topo_ptr)
std::reference_wrapper<DoubleVector const> phase_shift_;
std::reference_wrapper<SparseGroupedIdxVector const> load_gens_per_bus_;
std::reference_wrapper<DenseGroupedIdxVector const> sources_per_bus_;
std::reference_wrapper<std::vector<LoadGenType> const> load_gen_type_;
IterativePFSolver(YBus<sym> const& y_bus, MathModelTopology const& topo_ptr)
: n_bus_{y_bus.size()},
phase_shift_{topo_ptr, &topo_ptr->phase_shift},
load_gens_per_bus_{topo_ptr, &topo_ptr->load_gens_per_bus},
sources_per_bus_{topo_ptr, &topo_ptr->sources_per_bus},
load_gen_type_{topo_ptr, &topo_ptr->load_gen_type} {}
phase_shift_{std::cref(topo_ptr.phase_shift)},
load_gens_per_bus_{std::cref(topo_ptr.load_gens_per_bus)},
sources_per_bus_{std::cref(topo_ptr.sources_per_bus)},
load_gen_type_{std::cref(topo_ptr.load_gen_type)} {}
};

} // namespace power_grid_model::math_solver
Original file line number Diff line number Diff line change
Expand Up @@ -54,12 +54,12 @@ template <symmetry_tag sym_type> class LinearPFSolver {

static constexpr auto is_iterative = false;

LinearPFSolver(YBus<sym> const& y_bus, std::shared_ptr<MathModelTopology const> const& topo_ptr)
LinearPFSolver(YBus<sym> const& y_bus, MathModelTopology const& topo_ptr)
: n_bus_{y_bus.size()},
load_gens_per_bus_{topo_ptr, &topo_ptr->load_gens_per_bus},
sources_per_bus_{topo_ptr, &topo_ptr->sources_per_bus},
load_gens_per_bus_{std::cref(topo_ptr.load_gens_per_bus)},
sources_per_bus_{std::cref(topo_ptr.sources_per_bus)},
mat_data_(y_bus.nnz_lu()),
sparse_solver_{y_bus.shared_indptr_lu(), y_bus.shared_indices_lu(), y_bus.shared_diag_lu()},
sparse_solver_{y_bus.row_indptr_lu(), y_bus.col_indices_lu(), y_bus.lu_diag()},
perm_(n_bus_) {}

SolverOutput<sym> run_power_flow(YBus<sym> const& y_bus, PowerFlowInput<sym> const& input, Logger& log) {
Expand Down Expand Up @@ -92,20 +92,21 @@ template <symmetry_tag sym_type> class LinearPFSolver {
private:
Idx n_bus_;
// shared topo data
std::shared_ptr<SparseGroupedIdxVector const> load_gens_per_bus_;
std::shared_ptr<DenseGroupedIdxVector const> sources_per_bus_;
std::reference_wrapper<SparseGroupedIdxVector const> load_gens_per_bus_;
std::reference_wrapper<DenseGroupedIdxVector const> sources_per_bus_;
// sparse linear equation
ComplexTensorVector<sym> mat_data_;
// sparse solver
SparseSolverType sparse_solver_;
BlockPermArray perm_;

void prepare_matrix_and_rhs(YBus<sym> const& y_bus, PowerFlowInput<sym> const& input, SolverOutput<sym>& output) {
detail::prepare_linear_matrix_and_rhs(y_bus, input, *load_gens_per_bus_, *sources_per_bus_, output, mat_data_);
detail::prepare_linear_matrix_and_rhs(y_bus, input, load_gens_per_bus_.get(), sources_per_bus_.get(), output,
mat_data_);
}

void calculate_result(YBus<sym> const& y_bus, PowerFlowInput<sym> const& input, SolverOutput<sym>& output) {
detail::calculate_pf_result(y_bus, input, *sources_per_bus_, *load_gens_per_bus_, output,
detail::calculate_pf_result(y_bus, input, sources_per_bus_.get(), load_gens_per_bus_.get(), output,
[](Idx /*i*/) { return LoadGenType::const_y; });
}
};
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -87,7 +87,7 @@ template <symmetry_tag sym> class MathSolver : public MathSolverBase<sym> {
// construct model if needed
if (!iec60909_sc_solver_.has_value()) {
Timer const timer{log, LogEvent::create_math_solver};
iec60909_sc_solver_.emplace(y_bus, topo_ptr_);
iec60909_sc_solver_.emplace(y_bus, *topo_ptr_);
}

// call calculation
Expand Down Expand Up @@ -121,7 +121,7 @@ template <symmetry_tag sym> class MathSolver : public MathSolverBase<sym> {
Logger& log, YBus<sym> const& y_bus) {
if (!newton_raphson_pf_solver_.has_value()) {
Timer const timer{log, LogEvent::create_math_solver};
newton_raphson_pf_solver_.emplace(y_bus, topo_ptr_);
newton_raphson_pf_solver_.emplace(y_bus, *topo_ptr_);
}
return newton_raphson_pf_solver_.value().run_power_flow(y_bus, input, err_tol, max_iter, log);
}
Expand All @@ -130,7 +130,7 @@ template <symmetry_tag sym> class MathSolver : public MathSolverBase<sym> {
Logger& log, YBus<sym> const& y_bus) {
if (!linear_pf_solver_.has_value()) {
Timer const timer{log, LogEvent::create_math_solver};
linear_pf_solver_.emplace(y_bus, topo_ptr_);
linear_pf_solver_.emplace(y_bus, *topo_ptr_);
}
return linear_pf_solver_.value().run_power_flow(y_bus, input, log);
}
Expand All @@ -139,7 +139,7 @@ template <symmetry_tag sym> class MathSolver : public MathSolverBase<sym> {
Logger& log, YBus<sym> const& y_bus) {
if (!iterative_current_pf_solver_.has_value()) {
Timer const timer{log, LogEvent::create_math_solver};
iterative_current_pf_solver_.emplace(y_bus, topo_ptr_);
iterative_current_pf_solver_.emplace(y_bus, *topo_ptr_);
}
return iterative_current_pf_solver_.value().run_power_flow(y_bus, input, err_tol, max_iter, log);
}
Expand All @@ -154,7 +154,7 @@ template <symmetry_tag sym> class MathSolver : public MathSolverBase<sym> {
// construct model if needed
if (!iterative_linear_se_solver_.has_value()) {
Timer const timer{log, LogEvent::create_math_solver};
iterative_linear_se_solver_.emplace(y_bus, topo_ptr_);
iterative_linear_se_solver_.emplace(y_bus, *topo_ptr_);
}

// call calculation
Expand All @@ -166,7 +166,7 @@ template <symmetry_tag sym> class MathSolver : public MathSolverBase<sym> {
// construct model if needed
if (!newton_raphson_se_solver_.has_value()) {
Timer const timer{log, LogEvent::create_math_solver};
newton_raphson_se_solver_.emplace(y_bus, topo_ptr_);
newton_raphson_se_solver_.emplace(y_bus, *topo_ptr_);
}

// call calculation
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -36,7 +36,7 @@ template <symmetry_tag sym> class MeasuredValues {
public:
// construct
MeasuredValues(std::shared_ptr<MathModelTopology const> topo, StateEstimationInput<sym> const& input)
: math_topology_{std::move(topo)},
: math_topology_{*topo},
bus_appliance_injection_(math_topology().n_bus()),
idx_voltage_(math_topology().n_bus()),
bus_injection_(math_topology().n_bus()),
Expand Down Expand Up @@ -112,12 +112,12 @@ template <symmetry_tag sym> class MeasuredValues {

LoadGenSourceFlow calculate_load_gen_source(ComplexValueVector<sym> const& u,
ComplexValueVector<sym> const& s) const {
std::vector<ApplianceSolverOutput<sym>> load_gen_flow(math_topology_->n_load_gen());
std::vector<ApplianceSolverOutput<sym>> source_flow(math_topology_->n_source());
std::vector<ApplianceSolverOutput<sym>> load_gen_flow(math_topology_.n_load_gen());
std::vector<ApplianceSolverOutput<sym>> source_flow(math_topology_.n_source());

// loop all buses
for (auto const& [bus, load_gens, sources] :
enumerated_zip_sequence(math_topology_->load_gens_per_bus, math_topology_->sources_per_bus)) {
enumerated_zip_sequence(math_topology_.load_gens_per_bus, math_topology_.sources_per_bus)) {
// under-determined or exactly determined
if (bus_injection_[bus].n_unmeasured_appliances > 0) {
calculate_non_over_determined_injection(bus_injection_[bus].n_unmeasured_appliances, load_gens, sources,
Expand Down Expand Up @@ -173,7 +173,7 @@ template <symmetry_tag sym> class MeasuredValues {

private:
// cache topology
std::shared_ptr<MathModelTopology const> math_topology_;
MathModelTopology const& math_topology_;
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

member types should not be const& as it prevents it from being copyable and movable. Instead, please use either a raw pointer const* or a std::reference_wrapper


// flat arrays of all the relevant measurement for the main calculation
// branch/shunt flow, bus voltage, injection flow
Expand Down Expand Up @@ -214,7 +214,7 @@ template <symmetry_tag sym> class MeasuredValues {
// the lowest bus index with a voltage measurement
Idx first_voltage_measurement_{};

constexpr MathModelTopology const& math_topology() const { return *math_topology_; }
constexpr MathModelTopology const& math_topology() const { return math_topology_; }

void process_bus_related_measurements(StateEstimationInput<sym> const& input) {
/*
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -211,12 +211,12 @@ class NewtonRaphsonPFSolver : public IterativePFSolver<sym_type, NewtonRaphsonPF

static constexpr auto is_iterative = true;

NewtonRaphsonPFSolver(YBus<sym> const& y_bus, std::shared_ptr<MathModelTopology const> const& topo_ptr)
NewtonRaphsonPFSolver(YBus<sym> const& y_bus, MathModelTopology const& topo_ptr)
: IterativePFSolver<sym, NewtonRaphsonPFSolver>{y_bus, topo_ptr},
data_jac_(y_bus.nnz_lu()),
x_(y_bus.size()),
del_x_pq_(y_bus.size()),
sparse_solver_{y_bus.shared_indptr_lu(), y_bus.shared_indices_lu(), y_bus.shared_diag_lu()},
sparse_solver_{y_bus.row_indptr_lu(), y_bus.col_indices_lu(), y_bus.lu_diag()},
perm_(y_bus.size()) {}

// Initilize the unknown variable in polar form
Expand All @@ -225,13 +225,12 @@ class NewtonRaphsonPFSolver : public IterativePFSolver<sym_type, NewtonRaphsonPF
using LinearSparseSolverType = SparseLUSolver<ComplexTensor<sym>, ComplexValue<sym>, ComplexValue<sym>>;

ComplexTensorVector<sym> linear_mat_data(y_bus.nnz_lu());
LinearSparseSolverType linear_sparse_solver{y_bus.shared_indptr_lu(), y_bus.shared_indices_lu(),
y_bus.shared_diag_lu()};
LinearSparseSolverType linear_sparse_solver{y_bus.row_indptr_lu(), y_bus.col_indices_lu(), y_bus.lu_diag()};
typename LinearSparseSolverType::BlockPermArray linear_perm(y_bus.size());

detail::copy_y_bus<sym>(y_bus, linear_mat_data);
detail::prepare_linear_matrix_and_rhs(y_bus, input, *this->load_gens_per_bus_, *this->sources_per_bus_, output,
linear_mat_data);
detail::prepare_linear_matrix_and_rhs(y_bus, input, this->load_gens_per_bus_.get(),
this->sources_per_bus_.get(), output, linear_mat_data);
linear_sparse_solver.prefactorize_and_solve(linear_mat_data, linear_perm, output.u, output.u);

// get magnitude and angle of start voltage
Expand All @@ -244,13 +243,13 @@ class NewtonRaphsonPFSolver : public IterativePFSolver<sym_type, NewtonRaphsonPF
// Calculate the Jacobian and deviation
void prepare_matrix_and_rhs(YBus<sym> const& y_bus, PowerFlowInput<sym> const& input,
ComplexValueVector<sym> const& u) {
std::vector<LoadGenType> const& load_gen_type = *this->load_gen_type_;
std::vector<LoadGenType> const& load_gen_type = this->load_gen_type_.get();
IdxVector const& bus_entry = y_bus.lu_diag();

prepare_matrix_and_rhs_from_network_perspective(y_bus, u, bus_entry);

for (auto const& [bus_number, load_gens, sources] :
enumerated_zip_sequence(*this->load_gens_per_bus_, *this->sources_per_bus_)) {
enumerated_zip_sequence(this->load_gens_per_bus_.get(), this->sources_per_bus_.get())) {
Idx const diagonal_position = bus_entry[bus_number];
add_loads(load_gens, bus_number, diagonal_position, input, load_gen_type);
add_sources(sources, bus_number, diagonal_position, y_bus, input, u);
Expand Down
Loading
Loading