Merge pull request #131 from alliander-opensource/feature/three-phase-nan

Change update for three phase values with NaN's

Author: TonyXiang8787

include/power_grid_model/component/load_gen.hpp

@@ -78,15 +78,12 @@ class LoadGen final : public std::conditional_t<is_gen, GenericGenerator, Generi
     }
 
     void set_power(RealValue<sym> const& new_p_specified, RealValue<sym> const& new_q_specified) {
-        RealValue<sym> ps, qs;
-        if (is_nan(new_p_specified))
-            ps = real(s_specified_);
-        else
-            ps = direction_ / base_power<sym> * new_p_specified;
-        if (is_nan(new_q_specified))
-            qs = imag(s_specified_);
-        else
-            qs = direction_ / base_power<sym> * new_q_specified;
+        double const scalar = direction_ / base_power<sym>;
+        RealValue<sym> ps = real(s_specified_);
+        RealValue<sym> qs = imag(s_specified_);
+        update_real_value<sym>(new_p_specified, ps, scalar);
+        update_real_value<sym>(new_q_specified, qs, scalar);
+
         s_specified_ = ps + 1.0i * qs;
     }
 

include/power_grid_model/component/power_sensor.hpp

@@ -77,10 +77,13 @@ class PowerSensor : public GenericPowerSensor {
           {};
 
     UpdateChange update(PowerSensorUpdate<sym> const& power_sensor_update) {
-        if (!is_nan(power_sensor_update.p_measured) && !is_nan(power_sensor_update.q_measured)) {
-            s_measured_ = ((power_sensor_update.p_measured + 1i * power_sensor_update.q_measured) /
-                           base_power<sym>)*convert_direction();
-        }
+        double const scalar = convert_direction() / base_power<sym>;
+        RealValue<sym> ps = real(s_measured_);
+        RealValue<sym> qs = imag(s_measured_);
+        update_real_value<sym>(power_sensor_update.p_measured, ps, scalar);
+        update_real_value<sym>(power_sensor_update.q_measured, qs, scalar);
+        s_measured_ = ps + 1.0i * qs;
+
         if (!is_nan(power_sensor_update.power_sigma)) {
             power_sigma_ = power_sensor_update.power_sigma / base_power<sym>;
         }

include/power_grid_model/component/voltage_sensor.hpp

@@ -70,15 +70,14 @@ class VoltageSensor : public GenericVoltageSensor {
           u_angle_measured_{voltage_sensor_input.u_angle_measured} {};
 
     UpdateChange update(VoltageSensorUpdate<sym> const& voltage_sensor_update) {
-        if (!is_nan(voltage_sensor_update.u_measured)) {
-            u_measured_ = voltage_sensor_update.u_measured / (u_rated_ * u_scale<sym>);
-        }
+        double const scalar = 1 / (u_rated_ * u_scale<sym>);
+        update_real_value<sym>(voltage_sensor_update.u_measured, u_measured_, scalar);
+        update_real_value<sym>(voltage_sensor_update.u_angle_measured, u_angle_measured_, 1.0);
+
         if (!is_nan(voltage_sensor_update.u_sigma)) {
             u_sigma_ = voltage_sensor_update.u_sigma / (u_rated_ * u_scale<sym>);
         }
-        if (!is_nan(voltage_sensor_update.u_angle_measured)) {
-            u_angle_measured_ = voltage_sensor_update.u_angle_measured;
-        }
+
         return {false, false};
     }
 

include/power_grid_model/three_phase_tensor.hpp

@@ -321,6 +321,33 @@ inline bool is_nan(IntS x) {
     return x == na_IntS;
 }
 
+/* update real values
+
+   RealValue is only updated when the update value is not nan
+
+   symmetric:  update 1.0 with nan -> 1.0
+               update 1.0 with 2.0 -> 2.0
+
+   asymmetric: update [1.0, nan, nan] with [nan, nan, 2.0] -> [1.0, nan, 2.0]
+
+   The function assumes that the current value is normalized and new value should be normalized with scalar
+*/
+template <bool sym, class Proxy>
+void update_real_value(RealValue<sym> const& new_value, Proxy&& current_value, double scalar) {
+    if constexpr (sym) {
+        if (!is_nan(new_value)) {
+            current_value = scalar * new_value;
+        }
+    }
+    else {
+        for (size_t i = 0; i != 3; ++i) {
+            if (!is_nan(new_value(i))) {
+                current_value(i) = scalar * new_value(i);
+            }
+        }
+    }
+}
+
 // symmetric component matrix
 inline ComplexTensor<false> get_sym_matrix() {
     ComplexTensor<false> m;

src/power_grid_model/validation/utils.py

@@ -109,12 +109,19 @@ def update_component_data(component: str, input_data: np.ndarray, update_data: n
             mask = ~np.isnan(update_data[field])
         else:
             mask = np.not_equal(update_data[field], nan)
+
         if mask.ndim == 2:
-            mask = np.any(mask, axis=1)
-        data = update_data[["id", field]][mask]
-        idx = np.where(input_data["id"] == np.reshape(data["id"], (-1, 1)))
-        if isinstance(idx, tuple):
-            input_data[field][idx[1]] = data[field]
+            for phase in range(mask.shape[1]):
+                # find indexers of to-be-updated object
+                sub_mask = mask[:, phase]
+                idx = get_indexer(input_data["id"], update_data["id"][sub_mask])
+                # update
+                input_data[field][idx, phase] = update_data[field][sub_mask, phase]
+        else:
+            # find indexers of to-be-updated object
+            idx = get_indexer(input_data["id"], update_data["id"][mask])
+            # update
+            input_data[field][idx] = update_data[field][mask]
 
 
 def errors_to_string(
@@ -158,3 +165,25 @@ def nan_type(component: str, field: str, data_type="input"):
     It silently returns float('nan') if data_type/component/field can't be found.
     """
     return power_grid_meta_data.get(data_type, {}).get(component, {}).get("nans", {}).get(field, float("nan"))
+
+
+def get_indexer(input_ids: np.ndarray, update_ids: np.ndarray) -> np.ndarray:
+    """
+    Given array of ids from input and update dataset.
+    Find the posision of each id in the update dataset in the context of input dataset.
+    This is needed to update values in the dataset by id lookup.
+    Internally this is done by sorting the input ids, then using binary search lookup.
+
+    Args:
+        input_ids: array of ids in the input dataset
+        update_ids: array of ids in the update dataset
+
+    Returns:
+        np.ndarray: array of positions of the ids from update dataset in the input dataset
+            the following should hold
+            input_ids[result] == update_ids
+    """
+    permutation_sort = np.argsort(input_ids)  # complexity O(N_input * logN_input)
+    return permutation_sort[
+        np.searchsorted(input_ids, update_ids, sorter=permutation_sort)
+    ]  # complexity O(N_update * logN_input)

tests/cpp_unit_tests/test_load_gen.cpp

@@ -245,6 +245,40 @@ TEST_CASE("Test load generator") {
         CHECK(asym_result.q(1) == doctest::Approx(1e5));
     }
 
+    SUBCASE("Test set_power - sym") {
+        // update with nan, nothing happens
+        sym_gen_pq.set_power(RealValue<true>{nan}, RealValue<true>{nan});
+        ComplexValue<true> s_1 = sym_gen_pq.calc_param<true>();
+        CHECK(real(s_1) == 3.0);
+        CHECK(imag(s_1) == 3.0);
+
+        // update with values, s changes
+        sym_gen_pq.set_power(RealValue<true>{4.0e6}, RealValue<true>{5.0e6});
+        ComplexValue<true> s_2 = sym_gen_pq.calc_param<true>();
+        CHECK(real(s_2) == 4.0);
+        CHECK(imag(s_2) == 5.0);
+    }
+
+    SUBCASE("Test set_power - asym") {
+        // update with {nan, nan, nan}, nothing happens
+        asym_load_pq.set_power(RealValue<false>{nan}, RealValue<false>{nan});
+        ComplexValue<false> s_1 = asym_load_pq.calc_param<false>();
+        for (size_t i = 0; i != 3; i++) {
+            CHECK(real(s_1(i)) == -3.0);
+            CHECK(imag(s_1(i)) == -3.0);
+        }
+
+        // update some with nan, some with values
+        asym_load_pq.set_power(RealValue<false>{2.0e6, nan, 3.0e6}, RealValue<false>{nan, 4.0e6, nan});
+        ComplexValue<false> s_2 = asym_load_pq.calc_param<false>();
+        CHECK(real(s_2(0)) == -6.0);
+        CHECK(real(s_2(1)) == -3.0);  // not updated
+        CHECK(real(s_2(2)) == -9.0);
+        CHECK(imag(s_2(0)) == -3.0);  // not updated
+        CHECK(imag(s_2(1)) == -12.0);
+        CHECK(imag(s_2(2)) == -3.0);  // not updated
+    }
+
     SUBCASE("Test no source") {
         auto const s = sym_gen_pq.calc_param<false>(false);
         CHECK(real(s)(0) == doctest::Approx(0.0));