Merge branch 'main' into drop-python-3.10

Author: mgovers

.clang-tidy

@@ -55,6 +55,7 @@ readability-*,
 -readability-function-cognitive-complexity,
 -readability-identifier-length,
 -readability-magic-numbers,
+-bugprone-unchecked-optional-access,
 '
 
 WarningsAsErrors: '*'

CMakeLists.txt

@@ -29,7 +29,11 @@ set(CMAKE_LIBRARY_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/bin)
 
 include(GNUInstallDirs)
 
-cmake_policy(SET CMP0167 OLD)  # libboost-headers packaged by conda does not come with BoostConfig.cmake
+if(CMAKE_VERSION VERSION_GREATER_EQUAL "3.30.0") 
+#  libboost-headers packaged by conda does not come with BoostConfig.cmake
+cmake_policy(SET CMP0167 OLD)  
+endif()
+
 find_package(Boost REQUIRED)
 find_package(Eigen3 CONFIG REQUIRED)
 find_package(nlohmann_json CONFIG REQUIRED)

CMakePresets.json

@@ -254,7 +254,10 @@
     },
     {
       "name": "ci-clang-release",
-      "inherits": "clang-release"
+      "inherits": "clang-release",
+      "cacheVariables": {
+        "CMAKE_INTERPROCEDURAL_OPTIMIZATION": "TRUE"
+      }
     },
     {
       "name": "ci-gcc-debug",
@@ -270,7 +273,13 @@
       "name": "ci-gcc-release",
       "inherits": [
         "gcc-release"
-      ]
+      ],
+      "cacheVariables": {
+        "CMAKE_INTERPROCEDURAL_OPTIMIZATION": "TRUE"
+      },
+      "environment": {
+        "ASAN_OPTIONS": "detect_odr_violation=0"
+      }
     },
     {
       "name": "ci-sonar",

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

@@ -91,4 +91,10 @@ using IntSVector = std::vector<IntS>;
 template <class T, class... Ts>
 concept is_in_list_c = (std::same_as<std::remove_const_t<T>, Ts> || ...);
 
+// functor to include all
+struct IncludeAll {
+    template <class... T> constexpr bool operator()(T&&... /*ignored*/) const { return true; }
+};
+constexpr IncludeAll include_all{};
+
 } // namespace power_grid_model

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

@@ -108,7 +108,7 @@ class SparseMatrixError : public PowerGridError {
                    std::string("If you get this error from state estimation, ") +
                    "it might mean the system is not fully observable, i.e. not enough measurements.\n" +
                    "It might also mean that you are running into a corner case where PGM cannot resolve yet." +
-                   "See https://github.com/PowerGridModel/power-grid-model/issues/853.");
+                   "See https://github.com/PowerGridModel/power-grid-model/issues/864.");
     }
 };
 

power_grid_model_c/power_grid_model/include/power_grid_model/common/three_phase_tensor.hpp

@@ -152,9 +152,18 @@ inline ComplexValue<asymmetric_t> piecewise_complex_value(ComplexValue<asymmetri
 inline double cabs(double x) { return std::abs(x); }
 inline double cabs(DoubleComplex const& x) { return std::sqrt(std::norm(x)); }
 inline double abs2(DoubleComplex const& x) { return std::norm(x); }
-template <column_vector_or_tensor DerivedA> inline auto cabs(Eigen::ArrayBase<DerivedA> const& m) {
+template <column_vector_or_tensor DerivedA>
+inline auto cabs(Eigen::ArrayBase<DerivedA> const& m)
+    requires(std::same_as<typename DerivedA::Scalar, DoubleComplex>)
+{
     return sqrt(abs2(m));
 }
+template <column_vector_or_tensor DerivedA>
+inline auto cabs(Eigen::ArrayBase<DerivedA> const& m)
+    requires(std::same_as<typename DerivedA::Scalar, double>)
+{
+    return m.abs();
+}
 
 // phase_shift(x) = e^{i arg(x)} = x / |x|
 inline DoubleComplex phase_shift(DoubleComplex const x) {

power_grid_model_c/power_grid_model/include/power_grid_model/main_model_impl.hpp

@@ -1014,8 +1014,6 @@ class MainModelImpl<ExtraRetrievableTypes<ExtraRetrievableType...>, ComponentLis
         return math_param_increment;
     }
 
-    static constexpr auto include_all = [](Idx) { return true; };
-
     /** This is a heavily templated member function because it operates on many different variables of many
      *different types, but the essence is ever the same: filling one member (vector) of the calculation calc_input
      *struct (soa) with the right calculation symmetric or asymmetric calculation parameters, in the same order as
@@ -1068,7 +1066,7 @@ class MainModelImpl<ExtraRetrievableTypes<ExtraRetrievableType...>, ComponentLis
      */
     template <calculation_input_type CalcStructOut, typename CalcParamOut,
               std::vector<CalcParamOut>(CalcStructOut::*comp_vect), class ComponentIn,
-              std::invocable<Idx> PredicateIn = decltype(include_all)>
+              std::invocable<Idx> PredicateIn = IncludeAll>
         requires std::convertible_to<std::invoke_result_t<PredicateIn, Idx>, bool>
     static void prepare_input(MainModelState const& state, std::vector<Idx2D> const& components,
                               std::vector<CalcStructOut>& calc_input, PredicateIn include = include_all) {
@@ -1087,7 +1085,7 @@ class MainModelImpl<ExtraRetrievableTypes<ExtraRetrievableType...>, ComponentLis
 
     template <calculation_input_type CalcStructOut, typename CalcParamOut,
               std::vector<CalcParamOut>(CalcStructOut::*comp_vect), class ComponentIn,
-              std::invocable<Idx> PredicateIn = decltype(include_all)>
+              std::invocable<Idx> PredicateIn = IncludeAll>
         requires std::convertible_to<std::invoke_result_t<PredicateIn, Idx>, bool>
     static void prepare_input(MainModelState const& state, std::vector<Idx2D> const& components,
                               std::vector<CalcStructOut>& calc_input,

power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/iterative_linear_se_solver.hpp

@@ -95,11 +95,14 @@ template <symmetry_tag sym_type> class IterativeLinearSESolver {
         // preprocess measured value
         sub_timer = Timer(calculation_info, 2221, "Pre-process measured value");
         MeasuredValues<sym> const measured_values{y_bus.shared_topology(), input};
-        necessary_observability_check(measured_values, y_bus.math_topology(), y_bus.y_bus_structure());
+        auto const observability_result =
+            necessary_observability_check(measured_values, y_bus.math_topology(), y_bus.y_bus_structure());
 
-        // prepare matrix, including pre-factorization
+        // prepare matrix
         sub_timer = Timer(calculation_info, 2222, "Prepare matrix, including pre-factorization");
         prepare_matrix(y_bus, measured_values);
+        // prefactorize
+        sparse_solver_.prefactorize(data_gain_, perm_, observability_result.use_perturbation());
 
         // initialize voltage with initial angle
         sub_timer = Timer(calculation_info, 2223, "Initialize voltages");
@@ -252,8 +255,6 @@ template <symmetry_tag sym_type> class IterativeLinearSESolver {
             Idx const data_idx_tranpose = y_bus.lu_transpose_entry()[data_idx_lu];
             data_gain_[data_idx_lu].qh() = hermitian_transpose(data_gain_[data_idx_tranpose].q());
         }
-        // prefactorize
-        sparse_solver_.prefactorize(data_gain_, perm_);
     }
 
     void prepare_rhs(YBus<sym> const& y_bus, MeasuredValues<sym> const& measured_value,

power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/measured_values.hpp

@@ -498,9 +498,13 @@ template <symmetry_tag sym> class MeasuredValues {
     }
 
     // process one object
-    static constexpr auto default_status_checker = [](auto x) -> bool { return x; };
+    struct DefaultStatusChecker {
+        template <class T> bool operator()(T x) const { return x; }
+    };
+
+    static constexpr DefaultStatusChecker default_status_checker{};
 
-    template <class TS, class StatusChecker = decltype(default_status_checker)>
+    template <class TS, class StatusChecker = DefaultStatusChecker>
     static Idx process_one_object(Idx const object, grouped_idx_vector_type auto const& sensors_per_object,
                                   std::vector<TS> const& object_status,
                                   std::vector<PowerSensorCalcParam<sym>> const& input_data,

power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/newton_raphson_se_solver.hpp

@@ -159,7 +159,8 @@ template <symmetry_tag sym_type> class NewtonRaphsonSESolver {
         // preprocess measured value
         sub_timer = Timer(calculation_info, 2221, "Pre-process measured value");
         MeasuredValues<sym> const measured_values{y_bus.shared_topology(), input};
-        necessary_observability_check(measured_values, y_bus.math_topology(), y_bus.y_bus_structure());
+        auto const observability_result =
+            necessary_observability_check(measured_values, y_bus.math_topology(), y_bus.y_bus_structure());
 
         // initialize voltage with initial angle
         sub_timer = Timer(calculation_info, 2223, "Initialize voltages");
@@ -175,7 +176,8 @@ template <symmetry_tag sym_type> class NewtonRaphsonSESolver {
             prepare_matrix_and_rhs(y_bus, measured_values, output.u);
             // solve with prefactorization
             sub_timer = Timer(calculation_info, 2225, "Solve sparse linear equation");
-            sparse_solver_.prefactorize_and_solve(data_gain_, perm_, delta_x_rhs_, delta_x_rhs_);
+            sparse_solver_.prefactorize_and_solve(data_gain_, perm_, delta_x_rhs_, delta_x_rhs_,
+                                                  observability_result.use_perturbation());
             sub_timer = Timer(calculation_info, 2226, "Iterate unknown");
             max_dev = iterate_unknown(output.u, measured_values);
         };
@@ -519,7 +521,8 @@ template <symmetry_tag sym_type> class NewtonRaphsonSESolver {
     /// eta(row) += w_k . (z_k - f_k(x))
     ///
     /// In case there is no angle measurement, the slack bus or arbitray bus measurement is considered to have a virtual
-    /// angle measurement of zero. w_theta = 1.0 by default for all measurements
+    /// angle measurement of zero. w_theta = w_k by default for all measurements
+    ///    angle_error = u_error / u_rated (1.0) = w_k
     ///
     /// @param block LHS(row, col), ie. LHS(row, row)
     /// @param rhs_block RHS(row)
@@ -545,10 +548,10 @@ template <symmetry_tag sym_type> class NewtonRaphsonSESolver {
         RealValue<sym> delta_theta{};
         if (measured_values.has_angle_measurement(bus)) {
             delta_theta = RealValue<sym>{arg(measured_values.voltage(bus))} - RealValue<sym>{x_[bus].theta()};
-            w_theta = RealTensor<sym>{1.0};
+            w_theta = RealTensor<sym>{w_v};
         } else if (bus == virtual_angle_measurement_bus && !measured_values.has_angle()) {
             delta_theta = arg(ComplexValue<sym>{1.0}) - RealValue<sym>{x_[bus].theta()};
-            w_theta = RealTensor<sym>{1.0};
+            w_theta = RealTensor<sym>{w_v};
         }
 
         block.g_P_theta() += w_theta;