inverse update approach

docs/user_manual/calculations.md

@@ -435,8 +435,8 @@ pre-factorized, the computation cost of each iteration is much smaller than Newt
 method, where the Jacobian matrix needs to be constructed and factorized each time.
 
 ```{warning}
-Since the algorithm will assume angles to be zero if not given, this might result in not having a 
-crash due to an unobservable system, but succeeding with the calculations and giving faulty results.
+The algorithm will assume angles to be zero by default. This produces more correct outputs when the system is observable, but will
+prevent the calculation from raising an exception if it is unobservable, therefore giving faulty results.
 ```
 
 Algorithm call: `CalculationMethod.iterative_linear`

power_grid_model_c/power_grid_model/include/power_grid_model/auxiliary/dataset.hpp

@@ -69,7 +69,7 @@ template <bool is_const> class DataPointer {
         if (indptr_ != nullptr) {
             return indptr_[pos + 1] - indptr_[pos];
         }
-        return (Idx)elements_per_scenario_;
+        return elements_per_scenario_;
     }
 
     ptr_t<void> raw_ptr() const { return ptr_; }

power_grid_model_c/power_grid_model/include/power_grid_model/auxiliary/meta_data_gen.hpp

@@ -27,7 +27,7 @@ using DatasetMap = std::map<std::string, MetaComponent, std::less<>>;
 using AllDatasetMap = std::map<std::string, DatasetMap, std::less<>>;
 
 // template function to add meta data
-template <class CT> void add_meta_data(AllDatasetMap& meta) {
+template <class CT> inline void add_meta_data(AllDatasetMap& meta) {
     meta["input"].try_emplace(CT::name, MetaComponentImpl<typename CT::InputType>{}, CT::name);
     meta["update"].try_emplace(CT::name, MetaComponentImpl<typename CT::UpdateType>{}, CT::name);
     meta["sym_output"].try_emplace(CT::name, MetaComponentImpl<typename CT::template OutputType<true>>{}, CT::name);

power_grid_model_c/power_grid_model/include/power_grid_model/component/branch.hpp

@@ -148,13 +148,22 @@ class Branch : public Base {
     }
 
     // default update for branch, will be hidden for transformer
-    UpdateChange update(BranchUpdate const& update) {
-        assert(update.id == id());
-        bool const changed = set_status(update.from_status, update.to_status);
+    UpdateChange update(BranchUpdate const& update_data) {
+        assert(update_data.id == id());
+        bool const changed = set_status(update_data.from_status, update_data.to_status);
         // change branch connection will change both topo and param
         return {changed, changed};
     }
 
+    auto inverse(std::derived_from<BranchUpdate> auto update_data) const {
+        assert(update_data.id == id());
+
+        set_if_not_nan(update_data.from_status, static_cast<IntS>(from_status_));
+        set_if_not_nan(update_data.to_status, static_cast<IntS>(to_status_));
+
+        return update_data;
+    }
+
   protected:
     // calculate branch param based on symmetric component
     BranchCalcParam<true>

power_grid_model_c/power_grid_model/include/power_grid_model/component/branch3.hpp

@@ -160,13 +160,23 @@ class Branch3 : public Base {
     }
 
     // default update for branch3, will be hidden for three winding transformer
-    UpdateChange update(Branch3Update const& update) {
-        assert(update.id == id());
-        bool const changed = set_status(update.status_1, update.status_2, update.status_3);
+    UpdateChange update(Branch3Update const& update_data) {
+        assert(update_data.id == id());
+        bool const changed = set_status(update_data.status_1, update_data.status_2, update_data.status_3);
         // change in branch3 connection will change both topo and param
         return {changed, changed};
     }
 
+    auto inverse(std::derived_from<Branch3Update> auto update_data) const {
+        assert(update_data.id == id());
+
+        set_if_not_nan(update_data.status_1, static_cast<IntS>(status_1_));
+        set_if_not_nan(update_data.status_2, static_cast<IntS>(status_2_));
+        set_if_not_nan(update_data.status_3, static_cast<IntS>(status_3_));
+
+        return update_data;
+    }
+
   private:
     ID node_1_;
     ID node_2_;

power_grid_model_c/power_grid_model/include/power_grid_model/component/fault.hpp

@@ -111,6 +111,19 @@ class Fault final : public Base {
         return {false, false}; // topology and parameters do not change
     }
 
+    FaultUpdate inverse(FaultUpdate update_data) const {
+        assert(update_data.id == id());
+
+        set_if_not_nan(update_data.status, static_cast<IntS>(status_));
+        set_if_not_nan(update_data.fault_type, fault_type_);
+        set_if_not_nan(update_data.fault_phase, fault_phase_);
+        set_if_not_nan(update_data.fault_object, fault_object_);
+        set_if_not_nan(update_data.r_f, r_f_);
+        set_if_not_nan(update_data.x_f, x_f_);
+
+        return update_data;
+    }
+
     constexpr bool energized(bool is_connected_to_source) const override { return is_connected_to_source; }
 
     constexpr bool status() const { return status_; }

power_grid_model_c/power_grid_model/include/power_grid_model/component/load_gen.hpp

@@ -87,14 +87,26 @@ class LoadGen final : public std::conditional_t<is_gen, GenericGenerator, Generi
     }
 
     // update for load_gen
-    UpdateChange update(LoadGenUpdate<sym> const& update) {
-        assert(update.id == this->id());
-        this->set_status(update.status);
-        set_power(update.p_specified, update.q_specified);
+    UpdateChange update(LoadGenUpdate<sym> const& update_data) {
+        assert(update_data.id == this->id());
+        this->set_status(update_data.status);
+        set_power(update_data.p_specified, update_data.q_specified);
         // change load connection and/or value will not change topology or parameters
         return {false, false};
     }
 
+    LoadGenUpdate<sym> inverse(LoadGenUpdate<sym> update_data) const {
+        double const scalar = direction_ * base_power<sym>;
+
+        assert(update_data.id == this->id());
+
+        set_if_not_nan(update_data.status, static_cast<IntS>(this->status()));
+        set_if_not_nan(update_data.p_specified, real(s_specified_) * scalar);
+        set_if_not_nan(update_data.q_specified, imag(s_specified_) * scalar);
+
+        return update_data;
+    }
+
   private:
     ComplexValue<sym> s_specified_{}; // specified power injection
 

power_grid_model_c/power_grid_model/include/power_grid_model/component/node.hpp

@@ -29,6 +29,7 @@ class Node final : public Base {
 
     // update node, nothing happens here
     static constexpr UpdateChange update(BaseUpdate const& /* update_data */) { return {false, false}; }
+    static constexpr BaseUpdate inverse(BaseUpdate update_data) { return update_data; }
 
     // energized
     template <bool sym>

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

@@ -81,20 +81,32 @@ template <bool sym> class PowerSensor : public GenericPowerSensor {
         set_power(power_sensor_input.p_measured, power_sensor_input.q_measured);
     };
 
-    UpdateChange update(PowerSensorUpdate<sym> const& power_sensor_update) {
+    UpdateChange update(PowerSensorUpdate<sym> const& update_data) {
         constexpr double scalar = 1.0 / base_power<sym>;
 
-        set_power(power_sensor_update.p_measured, power_sensor_update.q_measured);
+        set_power(update_data.p_measured, update_data.q_measured);
 
-        if (!is_nan(power_sensor_update.power_sigma)) {
-            apparent_power_sigma_ = power_sensor_update.power_sigma * scalar;
+        if (!is_nan(update_data.power_sigma)) {
+            apparent_power_sigma_ = update_data.power_sigma * scalar;
         }
-        update_real_value<sym>(power_sensor_update.p_sigma, p_sigma_, scalar);
-        update_real_value<sym>(power_sensor_update.q_sigma, q_sigma_, scalar);
+        update_real_value<sym>(update_data.p_sigma, p_sigma_, scalar);
+        update_real_value<sym>(update_data.q_sigma, q_sigma_, scalar);
 
         return {false, false};
     }
 
+    PowerSensorUpdate<sym> inverse(PowerSensorUpdate<sym> update_data) const {
+        assert(update_data.id == this->id());
+
+        set_if_not_nan(update_data.p_measured, real(s_measured_) * base_power<sym>);
+        set_if_not_nan(update_data.q_measured, imag(s_measured_) * base_power<sym>);
+        set_if_not_nan(update_data.power_sigma, apparent_power_sigma_ * base_power<sym>);
+        set_if_not_nan(update_data.p_sigma, p_sigma_ * base_power<sym>);
+        set_if_not_nan(update_data.q_sigma, q_sigma_ * base_power<sym>);
+
+        return update_data;
+    }
+
   private:
     ComplexValue<sym> s_measured_{};
 

power_grid_model_c/power_grid_model/include/power_grid_model/component/shunt.hpp

@@ -47,16 +47,27 @@ class Shunt : public Appliance {
         }
     }
 
-    // update for shunt
-    UpdateChange update(ShuntUpdate const& update) {
-        assert(update.id == id());
-        bool changed = set_status(update.status);
-        changed = update_params(update) || changed;
+    UpdateChange update(ShuntUpdate const& update_data) {
+        assert(update_data.id == id());
+        bool changed = set_status(update_data.status);
+        changed = update_params(update_data) || changed;
 
         // change shunt connection will not change topology, but will change parameters
         return {false, changed};
     }
 
+    ShuntUpdate inverse(ShuntUpdate update_data) const {
+        assert(update_data.id == id());
+
+        set_if_not_nan(update_data.status, static_cast<IntS>(this->status()));
+        set_if_not_nan(update_data.g1, g1_);
+        set_if_not_nan(update_data.b1, b1_);
+        set_if_not_nan(update_data.g0, g0_);
+        set_if_not_nan(update_data.b0, b0_);
+
+        return update_data;
+    }
+
   private:
     double base_y_{nan};
 

power_grid_model_c/power_grid_model/include/power_grid_model/component/source.hpp

@@ -96,15 +96,25 @@ class Source : public Appliance {
     }
 
     // update for source
-    UpdateChange update(SourceUpdate const& update) {
-        assert(update.id == id());
-        bool const topo_changed = set_status(update.status);
-        bool const param_changed = set_u_ref(update.u_ref, update.u_ref_angle);
+    UpdateChange update(SourceUpdate const& update_data) {
+        assert(update_data.id == id());
+        bool const topo_changed = set_status(update_data.status);
+        bool const param_changed = set_u_ref(update_data.u_ref, update_data.u_ref_angle);
         // change source connection will change both topo and param
         // change u ref will change param
         return {topo_changed, param_changed || topo_changed};
     }
 
+    SourceUpdate inverse(SourceUpdate update_data) const {
+        assert(update_data.id == id());
+
+        set_if_not_nan(update_data.status, static_cast<IntS>(this->status()));
+        set_if_not_nan(update_data.u_ref, u_ref_);
+        set_if_not_nan(update_data.u_ref_angle, u_ref_angle_);
+
+        return update_data;
+    }
+
   private:
     double u_ref_;
     double u_ref_angle_;

power_grid_model_c/power_grid_model/include/power_grid_model/component/three_winding_transformer.hpp

@@ -117,13 +117,22 @@ class ThreeWindingTransformer : public Branch3 {
         return true;
     }
 
-    UpdateChange update(ThreeWindingTransformerUpdate const& update) {
-        assert(update.id == id());
-        bool const topo_changed = set_status(update.status_1, update.status_2, update.status_3);
-        bool const param_changed = set_tap(update.tap_pos) || topo_changed;
+    UpdateChange update(ThreeWindingTransformerUpdate const& update_data) {
+        assert(update_data.id == id());
+        bool const topo_changed = set_status(update_data.status_1, update_data.status_2, update_data.status_3);
+        bool const param_changed = set_tap(update_data.tap_pos) || topo_changed;
         return {topo_changed, param_changed};
     }
 
+    ThreeWindingTransformerUpdate inverse(ThreeWindingTransformerUpdate update_data) const {
+        assert(update_data.id == id());
+
+        update_data = Branch3::inverse(update_data);
+        set_if_not_nan(update_data.tap_pos, tap_pos_);
+
+        return update_data;
+    }
+
   private:
     // three winding transformer parameters
     double u1_;

power_grid_model_c/power_grid_model/include/power_grid_model/component/transformer.hpp

@@ -84,13 +84,22 @@ class Transformer : public Branch {
     }
 
     // update for transformer, hide default update for branch
-    UpdateChange update(TransformerUpdate const& update) {
-        assert(update.id == id());
-        bool const topo_changed = set_status(update.from_status, update.to_status);
-        bool const param_changed = set_tap(update.tap_pos) || topo_changed;
+    UpdateChange update(TransformerUpdate const& update_data) {
+        assert(update_data.id == id());
+        bool const topo_changed = set_status(update_data.from_status, update_data.to_status);
+        bool const param_changed = set_tap(update_data.tap_pos) || topo_changed;
         return {topo_changed, param_changed};
     }
 
+    TransformerUpdate inverse(TransformerUpdate update_data) const {
+        assert(update_data.id == id());
+
+        update_data = Branch::inverse(update_data);
+        set_if_not_nan(update_data.tap_pos, tap_pos_);
+
+        return update_data;
+    }
+
   private:
     // transformer parameter
     double u1_;

power_grid_model_c/power_grid_model/include/power_grid_model/component/voltage_sensor.hpp

@@ -79,18 +79,33 @@ template <bool sym> class VoltageSensor : public GenericVoltageSensor {
           u_measured_{voltage_sensor_input.u_measured / (u_rated_ * u_scale<sym>)},
           u_angle_measured_{voltage_sensor_input.u_angle_measured} {};
 
-    UpdateChange update(VoltageSensorUpdate<sym> const& voltage_sensor_update) {
+    UpdateChange update(VoltageSensorUpdate<sym> const& update_data) {
+        assert(update_data.id == this->id());
+
         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>);
+        update_real_value<sym>(update_data.u_measured, u_measured_, scalar);
+        update_real_value<sym>(update_data.u_angle_measured, u_angle_measured_, 1.0);
+
+        if (!is_nan(update_data.u_sigma)) {
+            u_sigma_ = update_data.u_sigma * scalar;
         }
 
         return {false, false};
     }
 
+    VoltageSensorUpdate<sym> inverse(VoltageSensorUpdate<sym> update_data) const {
+        assert(update_data.id == this->id());
+
+        double const scalar = u_rated_ * u_scale<sym>;
+
+        set_if_not_nan(update_data.u_measured, u_measured_ * scalar);
+        set_if_not_nan(update_data.u_angle_measured, u_angle_measured_);
+        set_if_not_nan(update_data.u_sigma, u_sigma_ * scalar);
+
+        return update_data;
+    }
+
   private:
     double u_rated_;
     double u_sigma_;

power_grid_model_c/power_grid_model/include/power_grid_model/container.hpp

@@ -186,24 +186,12 @@ class Container<RetrievableTypes<GettableTypes...>, StorageableTypes...> {
         cum_size_ = {accumulate_size_per_vector<GettableTypes>()...};
     };
 
-    // cache a Storagable item with index pos to restore to when restore_values() is called
-    template <class Storageable> void cache_item(Idx pos) {
-        auto const& value = get_raw<Storageable, Storageable>(pos);
-        auto& cached_vec = std::get<std::vector<std::pair<Idx, Storageable>>>(cached_reset_values_);
-
-        cached_vec.emplace_back(pos, value);
-    }
-
-    void restore_values() { (restore_values_impl<StorageableTypes>(), ...); }
-
   private:
     std::tuple<std::vector<StorageableTypes>...> vectors_;
     std::unordered_map<ID, Idx2D> map_;
     std::array<Idx, num_gettable> size_;
     std::array<std::array<Idx, num_storageable + 1>, num_gettable> cum_size_;
 
-    std::tuple<std::vector<std::pair<Idx, StorageableTypes>>...> cached_reset_values_; // indices + reset values
-
 #ifndef NDEBUG
     // set construction_complete is used for debug assertions only
     bool construction_complete_{false};
@@ -262,15 +250,6 @@ class Container<RetrievableTypes<GettableTypes...>, StorageableTypes...> {
         return res;
     }
 
-    template <class Storageable> void restore_values_impl() {
-        auto& cached_vec = std::get<std::vector<std::pair<Idx, Storageable>>>(cached_reset_values_);
-        for (auto it = cached_vec.crbegin(); it != cached_vec.crend(); ++it) {
-            auto const& cache = *it;
-            get_raw<Storageable, Storageable>(cache.first) = cache.second;
-        }
-        cached_vec.clear();
-    }
-
     // define iterator
     template <class Gettable>
     class Iterator : public boost::iterator_facade<Iterator<Gettable>, Gettable, boost::random_access_traversal_tag,

power_grid_model_c/power_grid_model/include/power_grid_model/main_core/input.hpp

@@ -17,8 +17,8 @@ namespace power_grid_model::main_core {
 // different selection based on component type
 template <std::derived_from<Base> Component, class ComponentContainer, std::forward_iterator ForwardIterator>
     requires model_component_state<MainModelState, ComponentContainer, Component>
-void add_component(MainModelState<ComponentContainer>& state, ForwardIterator begin, ForwardIterator end,
-                   double system_frequency) {
+inline void add_component(MainModelState<ComponentContainer>& state, ForwardIterator begin, ForwardIterator end,
+                          double system_frequency) {
     size_t const size = std::distance(begin, end);
     state.components.template reserve<Component>(size);
     // loop to add component