Using pivot perturbation for ill-conditioned state estimation

.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)

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/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/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;

power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/observability.hpp

@@ -32,16 +32,25 @@ std::tuple<Idx, Idx> count_voltage_sensors(const Idx n_bus, const MeasuredValues
 // lower triangle part is always zero
 // for diagonal part, it will be one if there is bus injection
 // for upper triangle part, it will be one if there is branch flow sensor and the branch is fully connected
+// return a pair of
+//      a vector of flow sensor count
+//      a boolean indicating if the system is possibly ill-conditioned
 template <symmetry_tag sym>
-std::vector<int8_t> count_flow_sensors(MeasuredValues<sym> const& measured_values, MathModelTopology const& topo,
-                                       YBusStructure const& y_bus_structure) {
+std::pair<std::vector<int8_t>, bool> count_flow_sensors(MeasuredValues<sym> const& measured_values,
+                                                        MathModelTopology const& topo,
+                                                        YBusStructure const& y_bus_structure) {
     Idx const n_bus{topo.n_bus()};
-    std::vector<int8_t> flow_sensors(y_bus_structure.row_indptr.back(), 0); // initialize all to zero
+    std::pair<std::vector<int8_t>, bool> result{};
+    auto& [flow_sensors, possibly_ill_conditioned] = result;
+    flow_sensors.resize(y_bus_structure.row_indptr.back(), 0); // initialize all to zero
+    possibly_ill_conditioned = false;
     for (Idx row = 0; row != n_bus; ++row) {
+        bool has_at_least_one_sensor{false};
         // lower triangle is ignored and kept as zero
         // diagonal for bus injection measurement
         if (measured_values.has_bus_injection(row)) {
             flow_sensors[y_bus_structure.bus_entry[row]] = 1;
+            has_at_least_one_sensor = true;
         }
         // upper triangle for branch flow measurement
         for (Idx ybus_index = y_bus_structure.bus_entry[row] + 1; ybus_index != y_bus_structure.row_indptr[row + 1];
@@ -57,13 +66,22 @@ std::vector<int8_t> count_flow_sensors(MeasuredValues<sym> const& measured_value
                     if ((measured_values.has_branch_from(branch) || measured_values.has_branch_to(branch)) &&
                         topo.branch_bus_idx[branch][0] != -1 && topo.branch_bus_idx[branch][1] != -1) {
                         flow_sensors[ybus_index] = 1;
+                        has_at_least_one_sensor = true;
                         break;
                     }
                 }
             }
         }
+        // check voltage sensor
+        if (measured_values.has_voltage(row) && measured_values.has_angle_measurement(row)) {
+            has_at_least_one_sensor = true;
+        }
+        // the system could be ill-conditioned if there is no flow sensor for one bus, except the last bus
+        if (!has_at_least_one_sensor && row != n_bus - 1) {
+            possibly_ill_conditioned = true;
+        }
     }
-    return flow_sensors;
+    return result;
 }
 
 // re-organize the flow sensor for radial grid without phasor measurement
@@ -103,17 +121,27 @@ inline void assign_injection_sensor_radial(YBusStructure const& y_bus_structure,
 
 } // namespace detail
 
+struct ObservabilityResult {
+    bool is_sufficiently_observable{false};
+    bool is_possibly_ill_conditioned{false};
+    constexpr bool use_perturbation() const { return is_possibly_ill_conditioned && is_sufficiently_observable; }
+};
+
 template <symmetry_tag sym>
-inline void necessary_observability_check(MeasuredValues<sym> const& measured_values, MathModelTopology const& topo,
-                                          YBusStructure const& y_bus_structure) {
+inline ObservabilityResult necessary_observability_check(MeasuredValues<sym> const& measured_values,
+                                                         MathModelTopology const& topo,
+                                                         YBusStructure const& y_bus_structure) {
     Idx const n_bus{topo.n_bus()};
+    ObservabilityResult result{};
 
     auto const [n_voltage_sensor, n_voltage_phasor_sensor] = detail::count_voltage_sensors(n_bus, measured_values);
     if (n_voltage_sensor < 1) {
         throw NotObservableError{"No voltage sensor found!\n"};
     }
 
-    std::vector<int8_t> flow_sensors = detail::count_flow_sensors(measured_values, topo, y_bus_structure);
+    std::vector<int8_t> flow_sensors;
+    std::tie(flow_sensors, result.is_possibly_ill_conditioned) =
+        detail::count_flow_sensors(measured_values, topo, y_bus_structure);
     // count flow sensors, note we manually specify the intial value type to avoid overflow
     Idx const n_flow_sensor = std::reduce(flow_sensors.cbegin(), flow_sensors.cend(), Idx{}, std::plus<Idx>{});
 
@@ -130,12 +158,16 @@ inline void necessary_observability_check(MeasuredValues<sym> const& measured_va
     if (topo.is_radial && n_voltage_phasor_sensor == 0) {
         detail::assign_injection_sensor_radial(y_bus_structure, flow_sensors);
         // count flow sensors again
-        Idx const n_flow_sensor_new = std::reduce(flow_sensors.cbegin(), flow_sensors.cend(), Idx{}, std::plus<Idx>{});
-        if (n_flow_sensor_new < n_bus - 1) {
+        if (Idx const n_flow_sensor_new =
+                std::reduce(flow_sensors.cbegin(), flow_sensors.cend(), Idx{}, std::plus<Idx>{});
+            n_flow_sensor_new < n_bus - 1) {
             throw NotObservableError{"The number of power sensors appears sufficient, but they are not independent "
                                      "enough. The system is still not observable.\n"};
         }
+        result.is_sufficiently_observable = true;
     }
+
+    return result;
 }
 
 } // namespace power_grid_model::math_solver