Repository: NREL/floris

Issue:
[BUG] Missed dimension reduction in `YawOptimization` base class
[These lines](https://github.com/NREL/floris/blob/440549ce6a603eb726971ccec645d3f7e86e249f/floris/optimization/yaw_optimization/yaw_optimization_base.py#L251-L264) 
```python
        idx = (minimum_yaw_angle_subset > 0.0) | (maximum_yaw_angle_subset < 0.0)
        if np.any(idx):
            # Find bounds closest to 0.0 deg
            combined_bounds = np.concatenate(
                (
                    np.expand_dims(minimum_yaw_angle_subset, axis=3),
                    np.expand_dims(maximum_yaw_angle_subset, axis=3)
                ),
                axis=3
            )
            # Overwrite all values that are not allowed to be 0.0 with bound value closest to zero
            ids_closest = np.expand_dims(np.argmin(np.abs(combined_bounds), axis=3), axis=3)
            yaw_mb = np.squeeze(np.take_along_axis(combined_bounds, ids_closest, axis=3))
            yaw_angles_template_subset[idx] = yaw_mb[idx]
```
were missed in the v4 dimension reduction process #764 , and can be fixed by replacing with
```python
        idx = (minimum_yaw_angle_subset > 0.0) | (maximum_yaw_angle_subset < 0.0)
        if np.any(idx):
            # Find bounds closest to 0.0 deg
            combined_bounds = np.concatenate(
                (
                    np.expand_dims(minimum_yaw_angle_subset, axis=2),
                    np.expand_dims(maximum_yaw_angle_subset, axis=2)
                ),
                axis=2
            )
            # Overwrite all values that are not allowed to be 0.0 with bound value closest to zero
            ids_closest = np.expand_dims(np.argmin(np.abs(combined_bounds), axis=2), axis=2)
            yaw_mb = np.squeeze(np.take_along_axis(combined_bounds, ids_closest, axis=2), axis=2)
            yaw_angles_template_subset[idx] = yaw_mb[idx]
```

(that is, switching from `axis=3` to `axis=2` throughout, and also specifying the axis for the `np.squeeze` operation).

The reason this hasn't come up before is that this code only runs if `minimum_yaw_angle` is strictly greater 0 or `maximum_yaw_angle` is strictly less than 0, which is an unusual circumstance.

I will create a bugfix PR in due course that will include tests for these "unusual" minimum and maximum angles.

Hints from issue comments:
@misi9170: Closed by #1067 

Solution patch:
diff --git a/examples/examples_control_optimization/002_opt_yaw_single_ws_uncertain.py b/examples/examples_control_optimization/002_opt_yaw_single_ws_uncertain.py
index 4b9ceda1e..1f6c80881 100644
--- a/examples/examples_control_optimization/002_opt_yaw_single_ws_uncertain.py
+++ b/examples/examples_control_optimization/002_opt_yaw_single_ws_uncertain.py
@@ -76,7 +76,7 @@
         color="r",
         marker="x",
     )
-    ax.set_ylabel("Yaw Offset (deg")
+    ax.set_ylabel("Yaw Offset T{0:03d} (deg)".format(tindex))
     ax.legend()
     ax.grid(True)
 
diff --git a/floris/optimization/yaw_optimization/yaw_optimization_base.py b/floris/optimization/yaw_optimization/yaw_optimization_base.py
index 358ecc40f..5dab7fe30 100644
--- a/floris/optimization/yaw_optimization/yaw_optimization_base.py
+++ b/floris/optimization/yaw_optimization/yaw_optimization_base.py
@@ -260,14 +260,14 @@ def _reduce_control_problem(self):
             # Find bounds closest to 0.0 deg
             combined_bounds = np.concatenate(
                 (
-                    np.expand_dims(minimum_yaw_angle_subset, axis=3),
-                    np.expand_dims(maximum_yaw_angle_subset, axis=3)
+                    np.expand_dims(minimum_yaw_angle_subset, axis=2),
+                    np.expand_dims(maximum_yaw_angle_subset, axis=2)
                 ),
-                axis=3
+                axis=2
             )
             # Overwrite all values that are not allowed to be 0.0 with bound value closest to zero
-            ids_closest = np.expand_dims(np.argmin(np.abs(combined_bounds), axis=3), axis=3)
-            yaw_mb = np.squeeze(np.take_along_axis(combined_bounds, ids_closest, axis=3))
+            ids_closest = np.expand_dims(np.argmin(np.abs(combined_bounds), axis=2), axis=2)
+            yaw_mb = np.squeeze(np.take_along_axis(combined_bounds, ids_closest, axis=2), axis=2)
             yaw_angles_template_subset[idx] = yaw_mb[idx]
 
         # Save all subset variables to self


Test patch:
diff --git a/tests/yaw_optimization_integration_test.py b/tests/yaw_optimization_integration_test.py
new file mode 100644
index 000000000..a0c3011fc
--- /dev/null
+++ b/tests/yaw_optimization_integration_test.py
@@ -0,0 +1,66 @@
+import numpy as np
+import pandas as pd
+import pytest
+
+from floris import FlorisModel
+from floris.optimization.yaw_optimization.yaw_optimizer_sr import YawOptimizationSR
+
+
+DEBUG = False
+VELOCITY_MODEL = "gauss"
+DEFLECTION_MODEL = "gauss"
+
+def test_yaw_optimization_limits(sample_inputs_fixture):
+    """
+    The Serial Refine (SR) method optimizes yaw angles based on a sequential, iterative yaw
+    optimization scheme. This test compares the optimization results from the SR method for
+    a simple farm with a simple wind rose to stored baseline results.
+    """
+    sample_inputs_fixture.core["wake"]["model_strings"]["velocity_model"] = VELOCITY_MODEL
+    sample_inputs_fixture.core["wake"]["model_strings"]["deflection_model"] = DEFLECTION_MODEL
+
+    fmodel = FlorisModel(sample_inputs_fixture.core)
+    wd_array = np.arange(0.0, 360.0, 90.0)
+    ws_array = 8.0 * np.ones_like(wd_array)
+    ti_array = 0.1 * np.ones_like(wd_array)
+
+    D = 126.0 # Rotor diameter for the NREL 5 MW
+    fmodel.set(
+        layout_x=[0.0, 5 * D, 10 * D],
+        layout_y=[0.0, 0.0, 0.0],
+        wind_directions=wd_array,
+        wind_speeds=ws_array,
+        turbulence_intensities=ti_array,
+    )
+
+    # Asymmetric limits
+    yaw_opt = YawOptimizationSR(
+        fmodel,
+        minimum_yaw_angle=-10.0,
+        maximum_yaw_angle=20.0,
+    )
+    yaw_opt.optimize()
+
+    # Strictly positive limits
+    yaw_opt = YawOptimizationSR(
+        fmodel,
+        minimum_yaw_angle=5.0,
+        maximum_yaw_angle=20.0,
+    )
+    yaw_opt.optimize()
+
+    # Strictly negative limits
+    yaw_opt = YawOptimizationSR(
+        fmodel,
+        minimum_yaw_angle=-20.0,
+        maximum_yaw_angle=-5.0,
+    )
+    yaw_opt.optimize()
+
+    # Infeasible limits
+    with pytest.raises(ValueError):
+        yaw_opt = YawOptimizationSR(
+            fmodel,
+            minimum_yaw_angle=20.0,
+            maximum_yaw_angle=5.0,
+        )