Merge pull request #917 from PowerGridModel/feature/refactor-measured-values-statistics

Refactor to DecomposedComplexRandVar

Author: mgovers

power_grid_model_c/power_grid_model/include/power_grid_model/calculation_parameters.hpp

@@ -88,15 +88,7 @@ template <symmetry_tag sym> using VoltageSensorCalcParam = UniformComplexRandVar
 // The value is the complex power
 //   * for appliances, it is always in injection direction
 //   * for branches, the direction is node -> branch
-template <symmetry_tag sym_type> struct PowerSensorCalcParam {
-    using sym = sym_type;
-
-    static constexpr bool symmetric{is_symmetric_v<sym>};
-
-    ComplexValue<sym> value{};
-    RealValue<sym> p_variance{}; // variance (sigma^2) of the error range of the active power, in p.u.
-    RealValue<sym> q_variance{}; // variance (sigma^2) of the error range of the reactive power, in p.u.
-};
+template <symmetry_tag sym> using PowerSensorCalcParam = DecomposedComplexRandVar<sym>;
 
 // current sensor calculation parameters for state estimation
 // The value is the complex current
@@ -107,9 +99,7 @@ template <symmetry_tag sym_type> struct CurrentSensorCalcParam {
     static constexpr bool symmetric{is_symmetric_v<sym>};
 
     AngleMeasurementType angle_measurement_type{};
-    ComplexValue<sym> value{};
-    RealValue<sym> i_real_variance{}; // variance (sigma^2) of the error range of real part of the current, in p.u.
-    RealValue<sym> i_imag_variance{}; // variance (sigma^2) of the error range of imaginary part of the current, in p.u.
+    DecomposedComplexRandVar<sym> measurement;
 };
 
 template <typename T>

power_grid_model_c/power_grid_model/include/power_grid_model/component/current_sensor.hpp

@@ -134,11 +134,8 @@ template <symmetry_tag current_sensor_symmetry_> class CurrentSensor : public Ge
     template <symmetry_tag sym_calc> CurrentSensorCalcParam<sym_calc> calc_decomposed_param() const {
         auto const i_polar = PolarComplexRandVar<current_sensor_symmetry>(
             {i_measured_, i_sigma_ * i_sigma_}, {i_angle_measured_, i_angle_sigma_ * i_angle_sigma_});
-        auto const i_decomposed = DecomposedComplexRandVar<sym_calc>(i_polar);
         return CurrentSensorCalcParam<sym_calc>{.angle_measurement_type = angle_measurement_type(),
-                                                .value = i_decomposed.value(),
-                                                .i_real_variance = i_decomposed.real_component.variance,
-                                                .i_imag_variance = i_decomposed.imag_component.variance};
+                                                .measurement = DecomposedComplexRandVar<sym_calc>(i_polar)};
     }
     CurrentSensorOutput<symmetric_t> get_sym_output(ComplexValue<symmetric_t> const& i) const final {
         return get_generic_output<symmetric_t>(i);

power_grid_model_c/power_grid_model/include/power_grid_model/component/power_sensor.hpp

@@ -131,34 +131,37 @@ template <symmetry_tag power_sensor_symmetry_> class PowerSensor : public Generi
     PowerSensorCalcParam<symmetric_t> sym_calc_param() const final {
         PowerSensorCalcParam<symmetric_t> calc_param{};
         if (is_normal(p_sigma_) && is_normal(q_sigma_)) {
-            calc_param.p_variance = mean_val(p_sigma_ * p_sigma_);
-            calc_param.q_variance = mean_val(q_sigma_ * q_sigma_);
+            calc_param.real_component.variance = mean_val(p_sigma_ * p_sigma_);
+            calc_param.imag_component.variance = mean_val(q_sigma_ * q_sigma_);
         } else {
             auto const variance = is_nan(p_sigma_) ? apparent_power_sigma_ * apparent_power_sigma_ / 2
                                                    : std::numeric_limits<double>::infinity();
-            calc_param.p_variance = variance;
-            calc_param.q_variance = variance;
+            calc_param.real_component.variance = variance;
+            calc_param.imag_component.variance = variance;
         }
-        assert(is_nan(calc_param.p_variance) == is_nan(calc_param.q_variance));
+        assert(is_nan(calc_param.real_component.variance) == is_nan(calc_param.imag_component.variance));
 
-        calc_param.value = mean_val(s_measured_);
+        auto const mean_s_measured = mean_val(s_measured_);
+        calc_param.real_component.value = real(mean_s_measured);
+        calc_param.imag_component.value = imag(mean_s_measured);
         return calc_param;
     }
     PowerSensorCalcParam<asymmetric_t> asym_calc_param() const final {
         PowerSensorCalcParam<asymmetric_t> calc_param{};
         if (is_normal(p_sigma_) && is_normal(q_sigma_)) {
-            calc_param.p_variance = RealValue<asymmetric_t>{p_sigma_ * p_sigma_};
-            calc_param.q_variance = RealValue<asymmetric_t>{q_sigma_ * q_sigma_};
+            calc_param.real_component.variance = RealValue<asymmetric_t>{p_sigma_ * p_sigma_};
+            calc_param.imag_component.variance = RealValue<asymmetric_t>{q_sigma_ * q_sigma_};
         } else {
             auto const variance =
                 RealValue<asymmetric_t>{is_nan(p_sigma_) ? apparent_power_sigma_ * apparent_power_sigma_ / 2
                                                          : std::numeric_limits<double>::infinity()};
-            calc_param.p_variance = variance;
-            calc_param.q_variance = variance;
+            calc_param.real_component.variance = variance;
+            calc_param.imag_component.variance = variance;
         }
-        assert(is_nan(calc_param.p_variance) == is_nan(calc_param.q_variance));
+        assert(is_nan(calc_param.real_component.variance) == is_nan(calc_param.imag_component.variance));
 
-        calc_param.value = piecewise_complex_value(s_measured_);
+        calc_param.real_component.value = RealValue<asymmetric_t>{real(s_measured_)};
+        calc_param.imag_component.value = RealValue<asymmetric_t>{imag(s_measured_)};
         return calc_param;
     }
     PowerSensorOutput<symmetric_t> get_sym_output(ComplexValue<symmetric_t> const& s) const final {

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

@@ -196,9 +196,10 @@ template <symmetry_tag sym_type> class IterativeLinearSESolver {
                         if (measured_value.has_shunt(obj)) {
                             // G += (-Ys)^H * (variance^-1) * (-Ys)
                             auto const& shunt_power = measured_value.shunt_power(obj);
-                            block.g() += dot(hermitian_transpose(param.shunt_param[obj]),
-                                             diagonal_inverse(shunt_power.p_variance + shunt_power.q_variance),
-                                             param.shunt_param[obj]);
+                            block.g() +=
+                                dot(hermitian_transpose(param.shunt_param[obj]),
+                                    diagonal_inverse(static_cast<IndependentComplexRandVar<sym>>(shunt_power).variance),
+                                    param.shunt_param[obj]);
                         }
                     }
                     // branch
@@ -214,7 +215,7 @@ template <symmetry_tag sym_type> class IterativeLinearSESolver {
                                 auto const& power = std::invoke(branch_power_[measured_side], measured_value, obj);
                                 block.g() +=
                                     dot(hermitian_transpose(param.branch_param[obj].value[measured_side * 2 + b0]),
-                                        diagonal_inverse(power.p_variance + power.q_variance),
+                                        diagonal_inverse(static_cast<IndependentComplexRandVar<sym>>(power).variance),
                                         param.branch_param[obj].value[measured_side * 2 + b1]);
                             }
                         }
@@ -229,8 +230,8 @@ template <symmetry_tag sym_type> class IterativeLinearSESolver {
                     if (row == col) {
                         // assign variance to diagonal of 3x3 tensor, for asym
                         auto const& injection = measured_value.bus_injection(row);
-                        block.r() = ComplexTensor<sym>{
-                            static_cast<ComplexValue<sym>>(-(injection.p_variance + injection.q_variance))};
+                        block.r() = ComplexTensor<sym>{static_cast<ComplexValue<sym>>(
+                            -(static_cast<IndependentComplexRandVar<sym>>(injection).variance))};
                     }
                 }
                 // injection measurement not exist
@@ -286,8 +287,10 @@ template <symmetry_tag sym_type> class IterativeLinearSESolver {
                     if (measured_value.has_shunt(obj)) {
                         PowerSensorCalcParam<sym> const& m = measured_value.shunt_power(obj);
                         // eta += (-Ys)^H * (variance^-1) * i_shunt
-                        rhs_block.eta() -= dot(hermitian_transpose(param.shunt_param[obj]),
-                                               diagonal_inverse(m.p_variance + m.q_variance), conj(m.value / u[bus]));
+                        rhs_block.eta() -=
+                            dot(hermitian_transpose(param.shunt_param[obj]),
+                                diagonal_inverse(static_cast<IndependentComplexRandVar<sym>>(m).variance),
+                                conj(m.value() / u[bus]));
                     }
                 }
                 // branch
@@ -307,14 +310,15 @@ template <symmetry_tag sym_type> class IterativeLinearSESolver {
                             // eta += Y{side, b}^H * (variance^-1) * i_branch_{f, t}
                             rhs_block.eta() +=
                                 dot(hermitian_transpose(param.branch_param[obj].value[measured_side * 2 + b]),
-                                    diagonal_inverse(m.p_variance + m.q_variance), conj(m.value / u[measured_bus]));
+                                    diagonal_inverse(static_cast<IndependentComplexRandVar<sym>>(m).variance),
+                                    conj(m.value() / u[measured_bus]));
                         }
                     }
                 }
             }
             // fill block with injection measurement, need to convert to current
             if (measured_value.has_bus_injection(bus)) {
-                rhs_block.tau() = conj(measured_value.bus_injection(bus).value / u[bus]);
+                rhs_block.tau() = conj(measured_value.bus_injection(bus).value() / u[bus]);
             }
         }
     }

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

@@ -332,8 +332,8 @@ template <symmetry_tag sym> class MeasuredValues {
 
         // combine valid appliance_injection_measurement and direct_injection_measurement
         // three scenarios; check if we have valid injection measurement
-        auto const uncertain_direct_injection =
-            is_inf(direct_injection_measurement.p_variance) || is_inf(direct_injection_measurement.q_variance);
+        auto const uncertain_direct_injection = is_inf(direct_injection_measurement.real_component.variance) ||
+                                                is_inf(direct_injection_measurement.imag_component.variance);
 
         bus_injection_[bus].idx_bus_injection = static_cast<Idx>(power_main_value_.size());
         if (n_unmeasured > 0) {
@@ -343,8 +343,8 @@ template <symmetry_tag sym> class MeasuredValues {
                 // only direct injection
                 power_main_value_.push_back(direct_injection_measurement);
             }
-        } else if (uncertain_direct_injection || any_zero(appliance_injection_measurement.p_variance) ||
-                   any_zero(appliance_injection_measurement.q_variance)) {
+        } else if (uncertain_direct_injection || any_zero(appliance_injection_measurement.real_component.variance) ||
+                   any_zero(appliance_injection_measurement.imag_component.variance)) {
             // only appliance injection if
             //    there is no direct injection measurement,
             //    or we have zero injection
@@ -370,13 +370,15 @@ template <symmetry_tag sym> class MeasuredValues {
         }
 
         auto const& appliance_measurement = extra_value_[appliance_idx];
-        if (is_inf(appliance_measurement.p_variance) || is_inf(appliance_measurement.q_variance)) {
+        if (is_inf(appliance_measurement.real_component.variance) ||
+            is_inf(appliance_measurement.imag_component.variance)) {
             ++n_unmeasured;
             return;
         }
-        measurements.value += appliance_measurement.value;
-        measurements.p_variance += appliance_measurement.p_variance;
-        measurements.q_variance += appliance_measurement.q_variance;
+        measurements.real_component.value += appliance_measurement.real_component.value;
+        measurements.imag_component.value += appliance_measurement.imag_component.value;
+        measurements.real_component.variance += appliance_measurement.real_component.variance;
+        measurements.imag_component.variance += appliance_measurement.imag_component.variance;
     }
 
     void process_branch_measurements(StateEstimationInput<sym> const& input) {
@@ -463,22 +465,25 @@ template <symmetry_tag sym> class MeasuredValues {
         for (auto pos : sensors) {
             auto const& measurement = data[pos];
 
-            accumulated_inverse_p_variance += RealValue<sym>{1.0} / measurement.p_variance;
-            accumulated_inverse_q_variance += RealValue<sym>{1.0} / measurement.q_variance;
+            accumulated_inverse_p_variance += RealValue<sym>{1.0} / measurement.real_component.variance;
+            accumulated_inverse_q_variance += RealValue<sym>{1.0} / measurement.imag_component.variance;
 
-            accumulated_p_value += real(measurement.value) / measurement.p_variance;
-            accumulated_q_value += imag(measurement.value) / measurement.q_variance;
+            accumulated_p_value += measurement.real_component.value / measurement.real_component.variance;
+            accumulated_q_value += measurement.imag_component.value / measurement.imag_component.variance;
         }
 
         if (is_normal(accumulated_inverse_p_variance) && is_normal(accumulated_inverse_q_variance)) {
-            return PowerSensorCalcParam<sym>{accumulated_p_value / accumulated_inverse_p_variance +
-                                                 1.0i * accumulated_q_value / accumulated_inverse_q_variance,
-                                             RealValue<sym>{1.0} / accumulated_inverse_p_variance,
-                                             RealValue<sym>{1.0} / accumulated_inverse_q_variance};
+            return PowerSensorCalcParam<sym>{
+                .real_component = {.value = accumulated_p_value / accumulated_inverse_p_variance,
+                                   .variance = RealValue<sym>{1.0} / accumulated_inverse_p_variance},
+                .imag_component = {.value = accumulated_q_value / accumulated_inverse_q_variance,
+                                   .variance = RealValue<sym>{1.0} / accumulated_inverse_q_variance}};
         }
-        return PowerSensorCalcParam<sym>{accumulated_p_value + 1.0i * accumulated_q_value,
-                                         RealValue<sym>{std::numeric_limits<double>::infinity()},
-                                         RealValue<sym>{std::numeric_limits<double>::infinity()}};
+        return PowerSensorCalcParam<sym>{
+            .real_component = {.value = accumulated_p_value,
+                               .variance = RealValue<sym>{std::numeric_limits<double>::infinity()}},
+            .imag_component = {.value = accumulated_q_value,
+                               .variance = RealValue<sym>{std::numeric_limits<double>::infinity()}}};
     }
 
     template <sensor_calc_param_type CalcParam, bool only_magnitude = false>
@@ -538,7 +543,7 @@ template <symmetry_tag sym> class MeasuredValues {
             unconstrained_min(x.variance);
         }
         for (auto const& x : power_main_value_) {
-            auto const variance = x.p_variance + x.q_variance;
+            auto const variance = x.real_component.variance + x.imag_component.variance;
             if constexpr (is_symmetric_v<sym>) {
                 unconstrained_min(variance);
             } else {
@@ -552,8 +557,8 @@ template <symmetry_tag sym> class MeasuredValues {
         auto const inv_norm_var = 1.0 / min_var;
         std::ranges::for_each(voltage_main_value_, [inv_norm_var](auto& x) { x.variance *= inv_norm_var; });
         std::ranges::for_each(power_main_value_, [inv_norm_var](auto& x) {
-            x.p_variance *= inv_norm_var;
-            x.q_variance *= inv_norm_var;
+            x.real_component.variance *= inv_norm_var;
+            x.imag_component.variance *= inv_norm_var;
         });
     }
 
@@ -563,17 +568,17 @@ template <symmetry_tag sym> class MeasuredValues {
                                                  FlowVector& source_flow) const {
         // calculate residual, divide, and assign to unmeasured (but connected) appliances
         ComplexValue<sym> const s_residual_per_appliance =
-            (s - bus_appliance_injection.value) / static_cast<double>(n_unmeasured);
+            (s - bus_appliance_injection.value()) / static_cast<double>(n_unmeasured);
         for (Idx const load_gen : load_gens) {
             if (has_load_gen(load_gen)) {
-                load_gen_flow[load_gen].s = load_gen_power(load_gen).value;
+                load_gen_flow[load_gen].s = load_gen_power(load_gen).value();
             } else if (idx_load_gen_power_[load_gen] == unmeasured) {
                 load_gen_flow[load_gen].s = s_residual_per_appliance;
             }
         }
         for (Idx const source : sources) {
             if (has_source(source)) {
-                source_flow[source].s = source_power(source).value;
+                source_flow[source].s = source_power(source).value();
             } else if (idx_source_power_[source] == unmeasured) {
                 source_flow[source].s = s_residual_per_appliance;
             }
@@ -586,13 +591,14 @@ template <symmetry_tag sym> class MeasuredValues {
                                              FlowVector& source_flow) const {
         // residual normalized by variance
         // mu = (sum[S_i] - S_cal) / sum[variance]
-        auto const delta = bus_appliance_injection.value - s;
-        ComplexValue<sym> const mu =
-            real(delta) / bus_appliance_injection.p_variance + 1.0i * imag(delta) / bus_appliance_injection.q_variance;
+        auto const delta = ComplexValue<sym>{bus_appliance_injection.value() - s};
+        ComplexValue<sym> const mu = real(delta) / bus_appliance_injection.real_component.variance +
+                                     1.0i * imag(delta) / bus_appliance_injection.imag_component.variance;
 
         // S_i = S_i_mea - var_i * mu
         auto const calculate_injection = [&mu](auto const& power) {
-            return power.value - (power.p_variance * real(mu) + 1.0i * power.q_variance * imag(mu));
+            return ComplexValue<sym>{power.value() - ((power.real_component.variance * real(mu)) +
+                                                      (1.0i * power.imag_component.variance * imag(mu)))};
         };
 
         for (Idx const load_gen : load_gens) {