Update branch

power_grid_model_c/power_grid_model/include/power_grid_model/common/exception.hpp

@@ -57,7 +57,7 @@ class InvalidBranch3 : public PowerGridError {
 class InvalidTransformerClock : public PowerGridError {
   public:
     InvalidTransformerClock(ID id, IntS clock) {
-        append_msg("Invalid clock for transformer " + std::to_string(id) + ", clock  " + std::to_string(clock) + '\n');
+        append_msg("Invalid clock for transformer " + std::to_string(id) + ", clock " + std::to_string(clock) + '\n');
     }
 };
 
@@ -137,13 +137,6 @@ class InvalidCalculationMethod : public CalculationError {
     InvalidCalculationMethod() : CalculationError("The calculation method is invalid for this calculation!") {}
 };
 
-class UnknownAttributeName : public PowerGridError {
-  public:
-    explicit UnknownAttributeName(std::string const& attr_name) {
-        append_msg("Unknown attribute name!" + attr_name + "\n");
-    }
-};
-
 class InvalidShortCircuitType : public PowerGridError {
   public:
     explicit InvalidShortCircuitType(FaultType short_circuit_type) {
@@ -168,7 +161,7 @@ class InvalidShortCircuitPhases : public PowerGridError {
 class InvalidShortCircuitPhaseOrType : public PowerGridError {
   public:
     InvalidShortCircuitPhaseOrType() {
-        append_msg("During one calculation the short circuit types phases should be similar for all faults \n");
+        append_msg("During one calculation the short circuit types phases should be similar for all faults\n");
     }
 };
 

power_grid_model_c/power_grid_model/include/power_grid_model/main_model.hpp

@@ -448,6 +448,8 @@ class MainModelImpl<ExtraRetrievableTypes<ExtraRetrievableType...>, ComponentLis
             calculation_fn(*this, {}, ignore_output);
         } catch (const SparseMatrixError&) {
             // missing entries are provided in the update data
+        } catch (const NotObservableError&) {
+            // missing entries are provided in the update data
         }
 
         // error messages

power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/common_solver_functions.hpp

@@ -23,6 +23,30 @@ inline void add_sources(IdxRange const& sources, Idx /* bus_number */, YBus<sym>
     }
 }
 
+template <symmetry_tag sym>
+inline void add_linear_loads(boost::iterator_range<IdxCount> const& load_gens_per_bus, Idx /* bus_number */,
+                             PowerFlowInput<sym> const& input, ComplexTensor<sym>& diagonal_element) {
+    for (auto load_number : load_gens_per_bus) {
+        // YBus_diag += -conj(S_base)
+        add_diag(diagonal_element, -conj(input.s_injection[load_number]));
+    }
+}
+
+template <symmetry_tag sym>
+inline void prepare_linear_matrix_and_rhs(YBus<sym> const& y_bus, PowerFlowInput<sym> const& input,
+                                          grouped_idx_vector_type auto const& load_gens_per_bus,
+                                          grouped_idx_vector_type auto const& sources_per_bus, MathOutput<sym>& output,
+                                          ComplexTensorVector<sym>& mat_data) {
+    IdxVector const& bus_entry = y_bus.lu_diag();
+    for (auto const& [bus_number, load_gens, sources] : enumerated_zip_sequence(load_gens_per_bus, sources_per_bus)) {
+        Idx const diagonal_position = bus_entry[bus_number];
+        auto& diagonal_element = mat_data[diagonal_position];
+        auto& u_bus = output.u[bus_number];
+        add_linear_loads(load_gens, bus_number, input, diagonal_element);
+        add_sources(sources, bus_number, y_bus, input.source, diagonal_element, u_bus);
+    }
+}
+
 template <symmetry_tag sym> inline 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) {

power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/iterative_current_pf_solver.hpp

@@ -73,16 +73,18 @@ namespace iterative_current_pf {
 template <symmetry_tag sym>
 class IterativeCurrentPFSolver : public IterativePFSolver<sym, IterativeCurrentPFSolver<sym>> {
   public:
-    using BlockPermArray =
-        typename SparseLUSolver<ComplexTensor<sym>, ComplexValue<sym>, ComplexValue<sym>>::BlockPermArray;
+    using SparseSolverType = SparseLUSolver<ComplexTensor<sym>, ComplexValue<sym>, ComplexValue<sym>>;
+    using BlockPermArray = typename SparseSolverType::BlockPermArray;
 
     IterativeCurrentPFSolver(YBus<sym> const& y_bus, std::shared_ptr<MathModelTopology const> const& topo_ptr)
         : 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()} {}
 
     // 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 */) {
+    void initialize_derived_solver(YBus<sym> const& y_bus, PowerFlowInput<sym> const& input, MathOutput<sym>& output) {
+        make_flat_start(input, output.u);
+
         auto const& sources_per_bus = *this->sources_per_bus_;
         IdxVector const& bus_entry = y_bus.lu_diag();
         // if Y bus is not up to date
@@ -149,7 +151,7 @@ class IterativeCurrentPFSolver : public IterativePFSolver<sym, IterativeCurrentP
     ComplexValueVector<sym> rhs_u_;
     std::shared_ptr<ComplexTensorVector<sym> const> mat_data_;
     // sparse solver
-    SparseLUSolver<ComplexTensor<sym>, ComplexValue<sym>, ComplexValue<sym>> sparse_solver_;
+    SparseSolverType sparse_solver_;
     std::shared_ptr<BlockPermArray const> perm_;
     bool parameters_changed_ = true;
 
@@ -187,6 +189,26 @@ class IterativeCurrentPFSolver : public IterativePFSolver<sym, IterativeCurrentP
                                       ComplexValue<sym>{input.source[source_number]});
         }
     }
+
+    void make_flat_start(PowerFlowInput<sym> const& input, ComplexValueVector<sym>& output_u) {
+        std::vector<double> const& phase_shift = *this->phase_shift_;
+        // 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 (Idx const source : sources) {
+                    sum_u_ref += input.source[source] * std::exp(1.0i * -phase_shift[bus]); // offset phase shift
+                }
+            }
+            return sum_u_ref / static_cast<double>(input.source.size());
+        }();
+
+        // assign u_ref as flat start
+        for (Idx i = 0; i != this->n_bus_; ++i) {
+            // consider phase shift
+            output_u[i] = ComplexValue<sym>{u_ref * std::exp(1.0i * phase_shift[i])};
+        }
+    }
 };
 
 template class IterativeCurrentPFSolver<symmetric_t>;

power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/iterative_pf_solver.hpp

@@ -28,7 +28,6 @@ template <symmetry_tag sym, typename DerivedSolver> class IterativePFSolver {
                                    Idx max_iter, CalculationInfo& calculation_info) {
         // get derived reference for derived solver class
         auto derived_solver = static_cast<DerivedSolver&>(*this);
-        std::vector<double> const& phase_shift = *phase_shift_;
 
         // prepare
         MathOutput<sym> output;
@@ -40,25 +39,8 @@ template <symmetry_tag sym, typename DerivedSolver> class IterativePFSolver {
         // initialize
         {
             Timer const sub_timer{calculation_info, 2221, "Initialize calculation"};
-            // average u_ref of all sources
-            DoubleComplex const u_ref = [&]() {
-                DoubleComplex sum_u_ref = 0.0;
-                for (auto const& [bus, sources] : enumerated_zip_sequence(*sources_per_bus_)) {
-                    for (Idx const source : sources) {
-                        sum_u_ref += input.source[source] * std::exp(1.0i * -phase_shift[bus]); // offset phase shift
-                    }
-                }
-                return sum_u_ref / (double)input.source.size();
-            }();
-
-            // assign u_ref as flat start
-            for (Idx i = 0; i != n_bus_; ++i) {
-                // consider phase shift
-                output.u[i] = ComplexValue<sym>{u_ref * std::exp(1.0i * phase_shift[i])};
-            }
-
             // Further initialization specific to the derived solver
-            derived_solver.initialize_derived_solver(y_bus, output);
+            derived_solver.initialize_derived_solver(y_bus, input, output);
         }
 
         // start calculation
@@ -94,7 +76,7 @@ template <symmetry_tag sym, typename DerivedSolver> class IterativePFSolver {
         main_timer.stop();
 
         auto const key = Timer::make_key(2226, "Max number of iterations");
-        calculation_info[key] = std::max(calculation_info[key], (double)num_iter);
+        calculation_info[key] = std::max(calculation_info[key], static_cast<double>(num_iter));
 
         return output;
     }

power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/linear_pf_solver.hpp

@@ -45,6 +45,10 @@ namespace linear_pf {
 template <symmetry_tag sym> class LinearPFSolver {
 
   public:
+    using SparseSolverType = SparseLUSolver<ComplexTensor<sym>, ComplexValue<sym>, ComplexValue<sym>>;
+    using BlockPermArray =
+        typename SparseLUSolver<ComplexTensor<sym>, ComplexValue<sym>, ComplexValue<sym>>::BlockPermArray;
+
     LinearPFSolver(YBus<sym> const& y_bus, std::shared_ptr<MathModelTopology const> const& topo_ptr)
         : n_bus_{y_bus.size()},
           load_gens_per_bus_{topo_ptr, &topo_ptr->load_gens_per_bus},
@@ -87,29 +91,11 @@ template <symmetry_tag sym> class LinearPFSolver {
     // sparse linear equation
     ComplexTensorVector<sym> mat_data_;
     // sparse solver
-    SparseLUSolver<ComplexTensor<sym>, ComplexValue<sym>, ComplexValue<sym>> sparse_solver_;
-    typename SparseLUSolver<ComplexTensor<sym>, ComplexValue<sym>, ComplexValue<sym>>::BlockPermArray perm_;
+    SparseSolverType sparse_solver_;
+    BlockPermArray perm_;
 
     void prepare_matrix_and_rhs(YBus<sym> const& y_bus, PowerFlowInput<sym> const& input, MathOutput<sym>& output) {
-        using detail::add_sources;
-
-        IdxVector const& bus_entry = y_bus.lu_diag();
-        for (auto const& [bus_number, load_gens, sources] :
-             enumerated_zip_sequence(*load_gens_per_bus_, *sources_per_bus_)) {
-            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_gens, bus_number, input, diagonal_element);
-            add_sources(sources, bus_number, y_bus, input.source, diagonal_element, u_bus);
-        }
-    }
-
-    static void add_loads(boost::iterator_range<IdxCount> const& load_gens_per_bus, Idx /* bus_number */,
-                          PowerFlowInput<sym> const& input, ComplexTensor<sym>& diagonal_element) {
-        for (auto load_number : load_gens_per_bus) {
-            // YBus_diag += -conj(S_base)
-            add_diag(diagonal_element, -conj(input.s_injection[load_number]));
-        }
+        detail::prepare_linear_matrix_and_rhs(y_bus, input, *load_gens_per_bus_, *sources_per_bus_, output, mat_data_);
     }
 
     void calculate_result(YBus<sym> const& y_bus, PowerFlowInput<sym> const& input, MathOutput<sym>& output) {

power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/newton_raphson_pf_solver.hpp

@@ -202,6 +202,10 @@ template <symmetry_tag sym> class PFJacBlock : public Block<double, sym, true, 2
 // solver
 template <symmetry_tag sym> class NewtonRaphsonPFSolver : public IterativePFSolver<sym, NewtonRaphsonPFSolver<sym>> {
   public:
+    using SparseSolverType = SparseLUSolver<PFJacBlock<sym>, ComplexPower<sym>, PolarPhasor<sym>>;
+    using BlockPermArray =
+        typename SparseLUSolver<PFJacBlock<sym>, ComplexPower<sym>, PolarPhasor<sym>>::BlockPermArray;
+
     NewtonRaphsonPFSolver(YBus<sym> const& y_bus, std::shared_ptr<MathModelTopology const> const& topo_ptr)
         : IterativePFSolver<sym, NewtonRaphsonPFSolver>{y_bus, topo_ptr},
           data_jac_(y_bus.nnz_lu()),
@@ -211,7 +215,19 @@ template <symmetry_tag sym> class NewtonRaphsonPFSolver : public IterativePFSolv
           perm_(y_bus.size()) {}
 
     // Initilize the unknown variable in polar form
-    void initialize_derived_solver(YBus<sym> const& /* y_bus */, MathOutput<sym> const& output) {
+    void initialize_derived_solver(YBus<sym> const& y_bus, PowerFlowInput<sym> const& input, MathOutput<sym>& output) {
+        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()};
+        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);
+        linear_sparse_solver.prefactorize_and_solve(linear_mat_data, linear_perm, output.u, output.u);
+
         // get magnitude and angle of start voltage
         for (Idx i = 0; i != this->n_bus_; ++i) {
             x_[i].v() = cabs(output.u[i]);
@@ -269,9 +285,10 @@ template <symmetry_tag sym> class NewtonRaphsonPFSolver : public IterativePFSolv
     // 2. power unbalance: p/q_specified - p/q_calculated
     // 3. unknown iterative
     std::vector<ComplexPower<sym>> del_x_pq_;
-    SparseLUSolver<PFJacBlock<sym>, ComplexPower<sym>, PolarPhasor<sym>> sparse_solver_;
+
+    SparseSolverType sparse_solver_;
     // permutation array
-    typename SparseLUSolver<PFJacBlock<sym>, ComplexPower<sym>, PolarPhasor<sym>>::BlockPermArray perm_;
+    BlockPermArray perm_;
 
     /// @brief power_flow_ij = (ui @* conj(uj))  .* conj(yij)
     /// Hij = diag(Vi) * ( Gij .* sin(theta_ij) - Bij .* cos(theta_ij) ) * diag(Vj)

power_grid_model_c/power_grid_model/include/power_grid_model/optimizer/tap_position_optimizer.hpp

@@ -101,7 +101,7 @@ inline auto get_edge_weights(TransformerGraph const& graph) -> WeightedTrafoList
 //       i. Infinity(INT_MAX), if tap side of the transformer is disconnected.
 //          The transformer regulation should be ignored
 //       ii.Rank of the vertex at the tap side of the transformer, if tap side of the transformer is connected
-inline auto rank_transformers(WeightedTrafoList const& w_trafo_list) -> std::vector<Idx2D> { return {}; }
+inline auto rank_transformers(WeightedTrafoList const& /*w_trafo_list*/) -> std::vector<Idx2D> { return {}; }
 
 template <main_core::main_model_state_c State> inline auto rank_transformers(State const& state) -> std::vector<Idx2D> {
     return rank_transformers(get_edge_weights(build_transformer_graph(state)));