Merge pull request #668 from PowerGridModel/feature/binary-search-tap…

…-changer

Feature / Binary search in Automatic Tap Changer

docs/examples/Transformer Examples.ipynb

@@ -40,7 +40,8 @@
     "import numpy as np\n",
     "import warnings\n",
     "\n",
-    "with warnings.catch_warnings(action=\"ignore\", category=DeprecationWarning):\n",
+    "with warnings.catch_warnings():\n",
+    "    warnings.simplefilter(\"ignore\", category=DeprecationWarning)\n",
     "    # suppress warning about pyarrow as future required dependency\n",
     "    import pandas as pd\n",
     "\n",
@@ -891,15 +892,160 @@
   },
   {
    "cell_type": "markdown",
-   "id": "7828a4e3",
+   "id": "9f1a6f30",
    "metadata": {},
    "source": [
-    "**NOTE:** the tap positions obtained using the `any_valid_tap` strategy may depend on the initial tap position of the transformers."
+    "You could also opt for fast_any_tap that takes advantage of binary search instead of linear search internally."
    ]
   },
   {
    "cell_type": "code",
    "execution_count": 12,
+   "id": "1fa08adb",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "------node result------\n"
+     ]
+    },
+    {
+     "data": {
+      "text/html": [
+       "<div>\n",
+       "<style scoped>\n",
+       "    .dataframe tbody tr th:only-of-type {\n",
+       "        vertical-align: middle;\n",
+       "    }\n",
+       "\n",
+       "    .dataframe tbody tr th {\n",
+       "        vertical-align: top;\n",
+       "    }\n",
+       "\n",
+       "    .dataframe thead th {\n",
+       "        text-align: right;\n",
+       "    }\n",
+       "</style>\n",
+       "<table border=\"1\" class=\"dataframe\">\n",
+       "  <thead>\n",
+       "    <tr style=\"text-align: right;\">\n",
+       "      <th></th>\n",
+       "      <th>id</th>\n",
+       "      <th>u</th>\n",
+       "    </tr>\n",
+       "  </thead>\n",
+       "  <tbody>\n",
+       "    <tr>\n",
+       "      <th>0</th>\n",
+       "      <td>2</td>\n",
+       "      <td>9999.994675</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>1</th>\n",
+       "      <td>4</td>\n",
+       "      <td>401.932237</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>2</th>\n",
+       "      <td>6</td>\n",
+       "      <td>346.203709</td>\n",
+       "    </tr>\n",
+       "  </tbody>\n",
+       "</table>\n",
+       "</div>"
+      ],
+      "text/plain": [
+       "   id            u\n",
+       "0   2  9999.994675\n",
+       "1   4   401.932237\n",
+       "2   6   346.203709"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    },
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "\n",
+      "------tap regulator result------\n",
+      "\n",
+      "----------fast_any_tap----------\n"
+     ]
+    },
+    {
+     "data": {
+      "text/html": [
+       "<div>\n",
+       "<style scoped>\n",
+       "    .dataframe tbody tr th:only-of-type {\n",
+       "        vertical-align: middle;\n",
+       "    }\n",
+       "\n",
+       "    .dataframe tbody tr th {\n",
+       "        vertical-align: top;\n",
+       "    }\n",
+       "\n",
+       "    .dataframe thead th {\n",
+       "        text-align: right;\n",
+       "    }\n",
+       "</style>\n",
+       "<table border=\"1\" class=\"dataframe\">\n",
+       "  <thead>\n",
+       "    <tr style=\"text-align: right;\">\n",
+       "      <th></th>\n",
+       "      <th>id</th>\n",
+       "      <th>energized</th>\n",
+       "      <th>tap_pos</th>\n",
+       "    </tr>\n",
+       "  </thead>\n",
+       "  <tbody>\n",
+       "    <tr>\n",
+       "      <th>0</th>\n",
+       "      <td>8</td>\n",
+       "      <td>1</td>\n",
+       "      <td>-1</td>\n",
+       "    </tr>\n",
+       "  </tbody>\n",
+       "</table>\n",
+       "</div>"
+      ],
+      "text/plain": [
+       "   id  energized  tap_pos\n",
+       "0   8          1       -1"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "# one-time power flow calculation with automatic tap changing\n",
+    "output_data6f = model5.calculate_power_flow(tap_changing_strategy=TapChangingStrategy.fast_any_tap)\n",
+    "\n",
+    "# the node at the control side of the transformer now has a voltage within the specified voltage band\n",
+    "print(\"------node result------\")\n",
+    "display(pd.DataFrame(output_data6f[ComponentType.node])[[\"id\", \"u\"]])\n",
+    "\n",
+    "print(\"\\n------tap regulator result------\")\n",
+    "print(\"\\n----------fast_any_tap----------\")\n",
+    "display(pd.DataFrame(output_data6f[ComponentType.transformer_tap_regulator]))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "7828a4e3",
+   "metadata": {},
+   "source": [
+    "**NOTE:** the tap positions obtained using the `any_valid_tap` and `fast_any_tap` strategy may depend on the initial tap position of the transformers."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 13,
    "id": "be276930",
    "metadata": {},
    "outputs": [
@@ -1094,7 +1240,7 @@
     "\n",
     "# power flow batch calculation with automatic tap changing\n",
     "output_data = model5.calculate_power_flow(\n",
-    "    update_data=update_data, tap_changing_strategy=TapChangingStrategy.any_valid_tap\n",
+    "    update_data=update_data, tap_changing_strategy=TapChangingStrategy.fast_any_tap\n",
     ")\n",
     "\n",
     "print(\"------node_4 batch result------\")\n",
@@ -1116,7 +1262,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 13,
+   "execution_count": 14,
    "id": "21b01a77",
    "metadata": {},
    "outputs": [
@@ -1259,7 +1405,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 14,
+   "execution_count": 15,
    "id": "c0314199",
    "metadata": {},
    "outputs": [
@@ -1412,7 +1558,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 15,
+   "execution_count": 16,
    "id": "ac091392",
    "metadata": {},
    "outputs": [],
@@ -1437,7 +1583,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 16,
+   "execution_count": 17,
    "id": "9b66d16e",
    "metadata": {},
    "outputs": [
@@ -1585,7 +1731,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.12.2"
+   "version": "3.10.14"
   },
   "vscode": {
    "interpreter": {

docs/user_manual/calculations.md

@@ -588,12 +588,13 @@ These {hoverxreftooltip}`user_manual/components:Transformer Tap Regulator`s try
 The $U_{\text{control}}$ may be compensated for the voltage drop during transport.
 Power flow calculations that take the behavior of these regulators into account may be toggled by providing one of the following strategies to the {py:meth}`tap_changing_strategy <power_grid_model.PowerGridModel.calculate_power_flow>` option.
 
-| Algorithm                                                              | Default  | Speed    | Algorithm call                                                                                              |
-| ---------------------------------------------------------------------- | -------- | -------- | ----------------------------------------------------------------------------------------------------------- |
-| No automatic tap changing (regular power flow)                         | &#10004; | &#10004; | {py:class}`TapChangingStrategy.disabled <power_grid_model.enum.TapChangingStrategy.disabled>`               |
-| Optimize tap positions for any value in the voltage band               |          | &#10004; | {py:class}`TapChangingStrategy.any_valid_tap <power_grid_model.enum.TapChangingStrategy.any_valid_tap>`     |
-| Optimize tap positions for lowest possible voltage in the voltage band |          |          | {py:class}`TapChangingStrategy.min_voltage_tap <power_grid_model.enum.TapChangingStrategy.min_voltage_tap>` |
-| Optimize tap positions for lowest possible voltage in the voltage band |          |          | {py:class}`TapChangingStrategy.max_voltage_tap <power_grid_model.enum.TapChangingStrategy.max_voltage_tap>` |
+| Algorithm                                                                   | Default  | Speed    | Algorithm call                                                                                              |
+| --------------------------------------------------------------------------- | -------- | -------- | ----------------------------------------------------------------------------------------------------------- |
+| No automatic tap changing (regular power flow)                              | &#10004; | &#10004; | {py:class}`TapChangingStrategy.disabled <power_grid_model.enum.TapChangingStrategy.disabled>`               |
+| Optimize tap positions for any value in the voltage band                    |          |          | {py:class}`TapChangingStrategy.any_valid_tap <power_grid_model.enum.TapChangingStrategy.any_valid_tap>`     |
+| Optimize tap positions for lowest possible voltage in the voltage band      |          |          | {py:class}`TapChangingStrategy.min_voltage_tap <power_grid_model.enum.TapChangingStrategy.min_voltage_tap>` |
+| Optimize tap positions for lowest possible voltage in the voltage band      |          |          | {py:class}`TapChangingStrategy.max_voltage_tap <power_grid_model.enum.TapChangingStrategy.max_voltage_tap>` |
+| Optimize tap positions for any value in the voltage band with binary search |          | &#10004; | {py:class}`TapChangingStrategy.fast_any_tap <power_grid_model.enum.TapChangingStrategy.fast_any_tap>`       |
 
 ##### Control logic for power flow with automatic tap changing
 
@@ -647,11 +648,21 @@ Hence, this assumption is reflected in the requirements mentioned in {hoverxreft
 
 Internally, to achieve an optimal regulated tap position, the control algorithm sets initial tap positions and exploits neighborhoods around local optima, depending on the strategy as follows.
 
-| strategy                                                                                                    | initial tap position | exploitation direction | description                                                                           |
-| ----------------------------------------------------------------------------------------------------------- | -------------------- | ---------------------- | ------------------------------------------------------------------------------------- |
-| {py:class}`TapChangingStrategy.any_valid_tap <power_grid_model.enum.TapChangingStrategy.any_valid_tap>`     | current tap position | no exploitation        | Find any tap position that gives a control side voltage within the `u_band`           |
-| {py:class}`TapChangingStrategy.min_voltage_tap <power_grid_model.enum.TapChangingStrategy.min_voltage_tap>` | `tap_max`            | step up                | Find the tap position that gives the lowest control side voltage within the `u_band`  |
-| {py:class}`TapChangingStrategy.max_voltage_tap <power_grid_model.enum.TapChangingStrategy.max_voltage_tap>` | `tap_min`            | step down              | Find the tap position that gives the highest control side voltage within the `u_band` |
+| strategy                                                                                                    | initial tap position | exploitation direction | search method | description                                                                           |
+| ----------------------------------------------------------------------------------------------------------- | -------------------- | ---------------------- | ------------- | ------------------------------------------------------------------------------------- |
+| {py:class}`TapChangingStrategy.any_valid_tap <power_grid_model.enum.TapChangingStrategy.any_valid_tap>`     | current tap position | no exploitation        | linear search | Find any tap position that gives a control side voltage within the `u_band`           |
+| {py:class}`TapChangingStrategy.min_voltage_tap <power_grid_model.enum.TapChangingStrategy.min_voltage_tap>` | `tap_max`            | step up                | binary search | Find the tap position that gives the lowest control side voltage within the `u_band`  |
+| {py:class}`TapChangingStrategy.max_voltage_tap <power_grid_model.enum.TapChangingStrategy.max_voltage_tap>` | `tap_min`            | step down              | binary search | Find the tap position that gives the highest control side voltage within the `u_band` |
+| {py:class}`TapChangingStrategy.fast_any_tap <power_grid_model.enum.TapChangingStrategy.fast_any_tap>`       | current tap position | no exploitation        | binary search | Find any tap position that gives a control side voltage within the `u_band`           |
+
+##### Search methods used for tap changing optimization
+
+Given the discrete nature of the finite tap ranges, we use the following search methods to find the next tap position along the exploitation direction.
+
+| Search method | Description                                                                            |
+| ------------- | -------------------------------------------------------------------------------------- |
+| linear search | Start with an initial guess and do a local search with step size 1 for each iteration step. |
+| binary search | Start with a large search region and reduce the search region by half for every iteration step. |
 
 
 ## Batch Calculations

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

@@ -25,7 +25,7 @@ using IntS = int8_t;
 // struct of indexing to sub modules
 struct Idx2D {
     Idx group; // sequence number of outer module/groups
-    Idx pos;   //  sequence number inside the group
+    Idx pos;   // sequence number inside the group
 
     friend constexpr bool operator==(Idx2D x, Idx2D y) = default;
 };

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

@@ -118,6 +118,12 @@ enum class OptimizerStrategy : IntS { // Conventions for optimization strategies
     global_maximum = 2,               // global_maximum = argmax{f(x) \in Range} for x in Domain
     local_minimum = 3,                // local_minimum = Any{argmin{f(x) \in Range}} for x in Domain
     local_maximum = 4,                // local_maximum = Any{argmax{f(x) \in Range}} for x in Domain
+    fast_any = 5,                     // fast_any = Any{f(x) \in Range} for x \in Domain, but faster
+};
+
+enum class SearchMethod : IntS { // Which type of tap search method for finite element optimization process
+    linear_search = 0,           // use linear_search method: one step per iteration
+    binary_search = 1,           // use binary search: half a tap range at a time
 };
 
 } // namespace power_grid_model

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

@@ -255,4 +255,13 @@ class UnreachableHit : public PowerGridError {
     }
 };
 
+class TapSearchStrategyIncompatibleError : public InvalidArguments {
+  public:
+    template <typename T1, typename T2>
+    TapSearchStrategyIncompatibleError(std::string const& method, const T1& value1, const T2& value2)
+        : InvalidArguments{
+              method, std::string{typeid(T1).name()} + " #" + detail::to_string(static_cast<IntS>(value1)) + " and " +
+                          std::string{typeid(T2).name()} + " #" + detail::to_string(static_cast<IntS>(value2))} {}
+};
+
 } // namespace power_grid_model

power_grid_model_c/power_grid_model/include/power_grid_model/main_model_fwd.hpp

@@ -21,7 +21,7 @@ struct MainModelOptions {
     CalculationSymmetry calculation_symmetry{CalculationSymmetry::symmetric};
     CalculationMethod calculation_method{CalculationMethod::default_method};
     OptimizerType optimizer_type{OptimizerType::no_optimization};
-    OptimizerStrategy optimizer_strategy{OptimizerStrategy::any};
+    OptimizerStrategy optimizer_strategy{OptimizerStrategy::fast_any};
 
     double err_tol{1e-8};
     Idx max_iter{20};

power_grid_model_c/power_grid_model/include/power_grid_model/main_model_impl.hpp

@@ -710,7 +710,6 @@ class MainModelImpl<ExtraRetrievableTypes<ExtraRetrievableType...>, ComponentLis
         return std::ranges::all_of(update_independent, [](bool const is_independent) { return is_independent; });
     }
 
-    // Single calculation
     template <calculation_type_tag calculation_type, symmetry_tag sym> auto calculate(Options const& options) {
         auto const calculator = [this, &options] {
             if constexpr (std::derived_from<calculation_type, power_flow_t>) {
@@ -726,10 +725,14 @@ class MainModelImpl<ExtraRetrievableTypes<ExtraRetrievableType...>, ComponentLis
             throw UnreachableHit{"MainModelImpl::calculate", "Unknown calculation type"};
         }();
 
+        SearchMethod const& search_method = options.optimizer_strategy == OptimizerStrategy::any
+                                                ? SearchMethod::linear_search
+                                                : SearchMethod::binary_search;
+
         return optimizer::get_optimizer<MainModelState, ConstDataset>(
                    options.optimizer_type, options.optimizer_strategy, calculator,
                    [this](ConstDataset update_data) { this->update_component<permanent_update_t>(update_data); },
-                   *meta_data_)
+                   *meta_data_, search_method)
             ->optimize(state_, options.calculation_method);
     }
 

power_grid_model_c/power_grid_model/include/power_grid_model/optimizer/optimizer.hpp

@@ -18,7 +18,7 @@ template <typename State, typename UpdateType, typename StateCalculator, typenam
     requires detail::state_calculator_c<StateCalculator, State> &&
              std::invocable<std::remove_cvref_t<StateUpdater>, UpdateType>
 constexpr auto get_optimizer(OptimizerType optimizer_type, OptimizerStrategy strategy, StateCalculator calculator,
-                             StateUpdater updater, meta_data::MetaData const& meta_data) {
+                             StateUpdater updater, meta_data::MetaData const& meta_data, SearchMethod search) {
     using enum OptimizerType;
     using namespace std::string_literals;
     using BaseOptimizer = detail::BaseOptimizer<StateCalculator, State>;
@@ -31,7 +31,7 @@ constexpr auto get_optimizer(OptimizerType optimizer_type, OptimizerStrategy str
                       std::invocable<std::remove_cvref_t<StateUpdater>, ConstDataset const&> &&
                       main_core::component_container_c<typename State::ComponentContainer, TransformerTapRegulator>) {
             return BaseOptimizer::template make_shared<TapPositionOptimizer<StateCalculator, StateUpdater, State>>(
-                std::move(calculator), std::move(updater), strategy, meta_data);
+                std::move(calculator), std::move(updater), strategy, meta_data, search);
         }
         [[fallthrough]];
     default: