feature/update main model for automatic tap changer

power_grid_model_c/power_grid_model/include/power_grid_model/calculation_parameters.hpp

@@ -126,9 +126,11 @@ struct TransformerTapRegulatorCalcParam {
     IntS status{};
 };
 
-struct TransformerTapRegulatorOptimizerOutput {
-    IntS tap_pos{na_IntS};
+struct TransformerTapPosition {
+    ID transformer_id{};
+    IntS tap_position{};
 };
+using TransformerTapPositionOutput = std::vector<TransformerTapPosition>;
 
 // from side, to side
 // in case of indices for math model, -1 means the branch is not connected to that side
@@ -269,7 +271,6 @@ template <symmetry_tag sym_type> struct SolverOutput {
     std::vector<ApplianceSolverOutput<sym>> source;
     std::vector<ApplianceSolverOutput<sym>> shunt;
     std::vector<ApplianceSolverOutput<sym>> load_gen;
-    std::vector<TransformerTapRegulatorOptimizerOutput> transformer_tap_regulator;
 };
 
 template <symmetry_tag sym_type> struct ShortCircuitSolverOutput {
@@ -325,13 +326,13 @@ static_assert(short_circuit_solver_output_type<ShortCircuitSolverOutput<symmetri
 static_assert(short_circuit_solver_output_type<ShortCircuitSolverOutput<asymmetric_t>>);
 
 struct OptimizerOutput {
-    std::vector<TransformerTapRegulatorOptimizerOutput> transformer_tap_regulator;
+    TransformerTapPositionOutput transformer_tap_positions;
 };
 
-template <solver_output_type T> struct MathOutput {
+template <typename T> struct MathOutput {
     using SolverOutputType = T;
 
-    std::vector<SolverOutputType> solver_output;
+    SolverOutputType solver_output;
     OptimizerOutput optimizer_output;
 };
 

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

@@ -30,6 +30,14 @@ struct Idx2D {
     friend constexpr bool operator==(Idx2D x, Idx2D y) = default;
 };
 
+struct Idx2DHash {
+    std::size_t operator()(const Idx2D& idx) const {
+        size_t const h1 = std::hash<Idx>{}(idx.group);
+        size_t const h2 = std::hash<Idx>{}(idx.pos);
+        return h1 ^ (h2 << 1);
+    }
+};
+
 struct symmetric_t {};
 struct asymmetric_t {};
 

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

@@ -139,6 +139,13 @@ class InvalidRegulatedObject : public PowerGridError {
     }
 };
 
+class DuplicativelyRegulatedObject : public PowerGridError {
+  public:
+    DuplicativelyRegulatedObject() {
+        append_msg("There are objects regulated by more than one regulator. Maximum one regulator is allowed.");
+    }
+};
+
 class AutomaticTapCalculationError : public PowerGridError {
   public:
     AutomaticTapCalculationError(ID id) {

power_grid_model_c/power_grid_model/include/power_grid_model/component/transformer_tap_regulator.hpp

@@ -56,13 +56,17 @@ class TransformerTapRegulator : public Regulator {
         return update_data;
     }
 
+    constexpr TransformerTapRegulatorOutput get_null_output() const {
+        return {.id = id(), .energized = 0, .tap_pos = na_IntS};
+    }
     TransformerTapRegulatorOutput get_output(IntS const& tap_pos) const {
         TransformerTapRegulatorOutput output{};
         output.id = id();
         output.energized = static_cast<IntS>(energized(true));
         output.tap_pos = tap_pos;
         return output;
     }
+    constexpr RegulatorShortCircuitOutput get_null_sc_output() const { return {.id = id(), .energized = 0}; }
 
     template <symmetry_tag sym> TransformerTapRegulatorCalcParam calc_param() const {
         TransformerTapRegulatorCalcParam param{};
@@ -77,7 +81,6 @@ class TransformerTapRegulator : public Regulator {
 
     // getter
     ControlSide control_side() const { return control_side_; }
-    RegulatorShortCircuitOutput get_null_sc_output() const { return {.id = id(), .energized = 0}; }
 
   private:
     // transformer tap regulator parameters

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

@@ -9,6 +9,8 @@
 
 #include "../all_components.hpp"
 
+#include <unordered_set>
+
 namespace power_grid_model::main_core {
 
 constexpr std::array<Branch3Side, 3> const branch3_sides = {Branch3Side::side_1, Branch3Side::side_2,
@@ -22,6 +24,8 @@ template <std::derived_from<Base> Component, class ComponentContainer, std::forw
 inline void add_component(MainModelState<ComponentContainer>& state, ForwardIterator begin, ForwardIterator end,
                           double system_frequency) {
     reserve_component<Component>(state, std::distance(begin, end));
+    // do sanity check on the transformer tap regulator
+    std::vector<Idx2D> regulated_objects;
     // loop to add component
     for (auto it = begin; it != end; ++it) {
         auto const& input = *it;
@@ -96,6 +100,7 @@ inline void add_component(MainModelState<ComponentContainer>& state, ForwardIter
             emplace_component<Component>(state, id, input);
         } else if constexpr (std::derived_from<Component, TransformerTapRegulator>) {
             Idx2D const regulated_object_idx = get_component_idx_by_id(state, input.regulated_object);
+            regulated_objects.push_back(regulated_object_idx);
 
             ID const regulated_terminal = [&input, &state, &regulated_object_idx] {
                 using enum ControlSide;
@@ -158,6 +163,13 @@ inline void add_component(MainModelState<ComponentContainer>& state, ForwardIter
             emplace_component<Component>(state, id, input, regulated_object_type, u_rated);
         }
     }
+    // Make sure that each regulated object has at most one regulator
+    const std::unordered_set<Idx2D, Idx2DHash> unique_regulated_objects(regulated_objects.begin(),
+                                                                        regulated_objects.end());
+    if (unique_regulated_objects.size() != regulated_objects.size()) {
+        // There are duplicates
+        throw DuplicativelyRegulatedObject{};
+    }
 }
 
 } // namespace power_grid_model::main_core

power_grid_model_c/power_grid_model/include/power_grid_model/main_core/output.hpp

@@ -335,19 +335,25 @@ inline auto output_result(Component const& fault, MainModelState<ComponentContai
 template <std::derived_from<TransformerTapRegulator> Component, class ComponentContainer,
           steady_state_solver_output_type SolverOutputType>
     requires model_component_state_c<MainModelState, ComponentContainer, Component>
-constexpr auto output_result(Component const& /* transformer_tap_regulator */,
+constexpr auto output_result(Component const& transformer_tap_regulator,
                              MainModelState<ComponentContainer> const& /* state */,
-                             std::vector<SolverOutputType> const& /* solver_output */, Idx const /* obj_seq */) {
-    // TODO: this function is not implemented
-    using sym = typename SolverOutputType::sym;
-    return typename TransformerTapRegulator::OutputType<sym>{};
+                             MathOutput<std::vector<SolverOutputType>> const& math_output, Idx const /* obj_seq */) {
+    if (auto const it = std::ranges::find_if(
+            math_output.optimizer_output.transformer_tap_positions,
+            [regulated_object = transformer_tap_regulator.regulated_object()](auto const& transformer_tap_pos) {
+                return transformer_tap_pos.transformer_id == regulated_object;
+            });
+        it != std::end(math_output.optimizer_output.transformer_tap_positions)) {
+        return transformer_tap_regulator.get_output(it->tap_position);
+    }
+    return transformer_tap_regulator.get_null_output();
 }
 template <std::derived_from<TransformerTapRegulator> Component, class ComponentContainer,
           short_circuit_solver_output_type SolverOutputType>
     requires model_component_state_c<MainModelState, ComponentContainer, Component>
-inline auto output_result(Component const& transformer_tap_regulator,
-                          MainModelState<ComponentContainer> const& /* state */,
-                          std::vector<SolverOutputType> const& /* solver_output */, Idx const /* obj_seq */) {
+constexpr auto
+output_result(Component const& transformer_tap_regulator, MainModelState<ComponentContainer> const& /* state */,
+              MathOutput<std::vector<SolverOutputType>> const& /* math_output */, Idx const /* obj_seq */) {
     return transformer_tap_regulator.get_null_sc_output();
 }
 
@@ -361,10 +367,10 @@ template <std::derived_from<Base> Component, class ComponentContainer, solver_ou
                      } -> std::convertible_to<std::iter_value_t<ResIt>>;
              }
 constexpr ResIt output_result(MainModelState<ComponentContainer> const& state,
-                              std::vector<SolverOutputType> const& solver_output, ResIt res_it) {
+                              MathOutput<std::vector<SolverOutputType>> const& math_output, ResIt res_it) {
     return detail::produce_output<Component, Idx2D>(
-        state, res_it, [&solver_output](Component const& component, Idx2D math_id) {
-            return output_result<Component>(component, solver_output, math_id);
+        state, res_it, [&math_output](Component const& component, Idx2D math_id) {
+            return output_result<Component>(component, math_output.solver_output, math_id);
         });
 }
 template <std::derived_from<Base> Component, class ComponentContainer, solver_output_type SolverOutputType,
@@ -377,10 +383,10 @@ template <std::derived_from<Base> Component, class ComponentContainer, solver_ou
                      } -> std::convertible_to<std::iter_value_t<ResIt>>;
              }
 constexpr ResIt output_result(MainModelState<ComponentContainer> const& state,
-                              std::vector<SolverOutputType> const& solver_output, ResIt res_it) {
+                              MathOutput<std::vector<SolverOutputType>> const& math_output, ResIt res_it) {
     return detail::produce_output<Component, Idx2D>(
-        state, res_it, [&state, &solver_output](Component const& component, Idx2D const math_id) {
-            return output_result<Component>(component, state, solver_output, math_id);
+        state, res_it, [&state, &math_output](Component const& component, Idx2D const math_id) {
+            return output_result<Component>(component, state, math_output.solver_output, math_id);
         });
 }
 template <std::derived_from<Base> Component, class ComponentContainer, solver_output_type SolverOutputType,
@@ -393,10 +399,10 @@ template <std::derived_from<Base> Component, class ComponentContainer, solver_ou
                      } -> std::convertible_to<std::iter_value_t<ResIt>>;
              }
 constexpr ResIt output_result(MainModelState<ComponentContainer> const& state,
-                              std::vector<SolverOutputType> const& solver_output, ResIt res_it) {
+                              MathOutput<std::vector<SolverOutputType>> const& math_output, ResIt res_it) {
     return detail::produce_output<Component, Idx>(
-        state, res_it, [&state, &solver_output](Component const& component, Idx const obj_seq) {
-            return output_result<Component, ComponentContainer>(component, state, solver_output, obj_seq);
+        state, res_it, [&state, &math_output](Component const& component, Idx const obj_seq) {
+            return output_result<Component, ComponentContainer>(component, state, math_output.solver_output, obj_seq);
         });
 }
 template <std::derived_from<Base> Component, class ComponentContainer, solver_output_type SolverOutputType,
@@ -409,10 +415,26 @@ template <std::derived_from<Base> Component, class ComponentContainer, solver_ou
                      } -> std::convertible_to<std::iter_value_t<ResIt>>;
              }
 constexpr ResIt output_result(MainModelState<ComponentContainer> const& state,
-                              std::vector<SolverOutputType> const& solver_output, ResIt res_it) {
+                              MathOutput<std::vector<SolverOutputType>> const& math_output, ResIt res_it) {
     return detail::produce_output<Component, Idx2DBranch3>(
-        state, res_it, [&solver_output](Component const& component, Idx2DBranch3 const& math_id) {
-            return output_result<Component>(component, solver_output, math_id);
+        state, res_it, [&math_output](Component const& component, Idx2DBranch3 const& math_id) {
+            return output_result<Component>(component, math_output.solver_output, math_id);
+        });
+}
+template <std::derived_from<Base> Component, class ComponentContainer, typename SolverOutputType,
+          std::forward_iterator ResIt>
+    requires model_component_state_c<MainModelState, ComponentContainer, Component> &&
+             requires(Component const& component, MainModelState<ComponentContainer> const& state,
+                      MathOutput<SolverOutputType> const& math_output, Idx const obj_seq) {
+                 {
+                     output_result<Component>(component, state, math_output, obj_seq)
+                     } -> std::convertible_to<std::iter_value_t<ResIt>>;
+             }
+constexpr ResIt output_result(MainModelState<ComponentContainer> const& state,
+                              MathOutput<SolverOutputType> const& math_output, ResIt res_it) {
+    return detail::produce_output<Component, Idx>(
+        state, res_it, [&state, &math_output](Component const& component, Idx const obj_seq) {
+            return output_result<Component, ComponentContainer>(component, state, math_output, obj_seq);
         });
 }
 
@@ -421,10 +443,10 @@ template <std::same_as<Appliance> Component, class ComponentContainer, solver_ou
           std::forward_iterator ResIt>
     requires model_component_state_c<MainModelState, ComponentContainer, Component>
 constexpr ResIt output_result(MainModelState<ComponentContainer> const& state,
-                              std::vector<SolverOutputType> const& solver_output, ResIt res_it) {
-    res_it = output_result<Source>(state, solver_output, res_it);
-    res_it = output_result<GenericLoadGen>(state, solver_output, res_it);
-    res_it = output_result<Shunt>(state, solver_output, res_it);
+                              MathOutput<std::vector<SolverOutputType>> const& math_output, ResIt res_it) {
+    res_it = output_result<Source>(state, math_output, res_it);
+    res_it = output_result<GenericLoadGen>(state, math_output, res_it);
+    res_it = output_result<Shunt>(state, math_output, res_it);
     return res_it;
 }
 

power_grid_model_c/power_grid_model/include/power_grid_model/main_model.hpp

@@ -678,13 +678,11 @@ class MainModelImpl<ExtraRetrievableTypes<ExtraRetrievableType...>, ComponentLis
     }
 
     template <symmetry_tag sym> auto calculate_power_flow(Options const& options) {
-        auto result_pf =
-            optimizer::get_optimizer<MainModelState, ConstDataset>(
-                options.optimizer_type, options.optimizer_strategy,
-                calculate_power_flow_<sym>(options.err_tol, options.max_iter),
-                [this](ConstDataset update_data) { this->update_component<permanent_update_t>(update_data); })
-                ->optimize(state_, options.calculation_method);
-        return MathOutput<SolverOutput<sym>>{.solver_output = std::move(result_pf), .optimizer_output = {}};
+        return optimizer::get_optimizer<MainModelState, ConstDataset>(
+                   options.optimizer_type, options.optimizer_strategy,
+                   calculate_power_flow_<sym>(options.err_tol, options.max_iter),
+                   [this](ConstDataset update_data) { this->update_component<permanent_update_t>(update_data); })
+            ->optimize(state_, options.calculation_method);
     }
 
     // Single load flow calculation, propagating the results to result_data
@@ -716,9 +714,10 @@ class MainModelImpl<ExtraRetrievableTypes<ExtraRetrievableType...>, ComponentLis
 
     // Single state estimation calculation, returning math output results
     template <symmetry_tag sym> auto calculate_state_estimation(Options const& options) {
-        return MathOutput<SolverOutput<sym>>{.solver_output = calculate_state_estimation_<sym>(
-                                                 options.err_tol, options.max_iter)(state_, options.calculation_method),
-                                             .optimizer_output = {}};
+        return MathOutput<std::vector<SolverOutput<sym>>>{
+            .solver_output =
+                calculate_state_estimation_<sym>(options.err_tol, options.max_iter)(state_, options.calculation_method),
+            .optimizer_output = {}};
     }
 
     // Single state estimation calculation, propagating the results to result_data
@@ -750,7 +749,7 @@ class MainModelImpl<ExtraRetrievableTypes<ExtraRetrievableType...>, ComponentLis
 
     // Single short circuit calculation, returning short circuit math output results
     template <symmetry_tag sym> auto calculate_short_circuit(Options const& options) {
-        return MathOutput<ShortCircuitSolverOutput<sym>>{
+        return MathOutput<std::vector<ShortCircuitSolverOutput<sym>>>{
             .solver_output = calculate_short_circuit_<sym>(options.short_circuit_voltage_scaling)(
                 state_, options.calculation_method),
             .optimizer_output = {}};
@@ -783,30 +782,29 @@ class MainModelImpl<ExtraRetrievableTypes<ExtraRetrievableType...>, ComponentLis
     }
 
     template <typename Component, typename MathOutputType, std::forward_iterator ResIt>
-        requires solver_output_type<typename MathOutputType::SolverOutputType>
+        requires solver_output_type<typename MathOutputType::SolverOutputType::value_type>
     ResIt output_result(MathOutputType const& math_output, ResIt res_it) const {
         assert(construction_complete_);
-        return main_core::output_result<Component, ComponentContainer>(state_, math_output.solver_output, res_it);
+        return main_core::output_result<Component, ComponentContainer>(state_, math_output, res_it);
     }
 
     template <solver_output_type SolverOutputType>
-    void output_result(MathOutput<SolverOutputType> const& math_output, Dataset const& result_data, Idx pos = 0) {
-        using OutputFunc = void (*)(MainModelImpl & x, MathOutput<SolverOutputType> const& math_output,
+    void output_result(MathOutput<std::vector<SolverOutputType>> const& math_output, Dataset const& result_data,
+                       Idx pos = 0) {
+        using OutputFunc = void (*)(MainModelImpl & x, MathOutput<std::vector<SolverOutputType>> const& math_output,
                                     MutableDataPointer const& data_ptr, Idx position);
 
-        static constexpr std::array<OutputFunc, n_types> get_result{[](MainModelImpl& model,
-                                                                       MathOutput<SolverOutputType> const& math_output_,
-                                                                       MutableDataPointer const& data_ptr,
-                                                                       Idx position) {
-            auto const begin =
-                data_ptr
-                    .get_iterators<
-                        std::conditional_t<steady_state_solver_output_type<SolverOutputType>,
+        static constexpr std::array<OutputFunc, n_types> get_result{
+            [](MainModelImpl& model, MathOutput<std::vector<SolverOutputType>> const& math_output_,
+               MutableDataPointer const& data_ptr, Idx position) {
+                auto const begin = data_ptr
+                                       .get_iterators<std::conditional_t<
+                                           steady_state_solver_output_type<SolverOutputType>,
                                            typename ComponentType::template OutputType<typename SolverOutputType::sym>,
                                            typename ComponentType::ShortCircuitOutputType>>(position)
-                    .first;
-            model.output_result<ComponentType>(math_output_, begin);
-        }...};
+                                       .first;
+                model.output_result<ComponentType>(math_output_, begin);
+            }...};
 
         Timer const t_output(calculation_info_, 3000, "Produce output");
         for (ComponentEntry const& entry : AllComponents::component_index_map) {