Repository: NREL/floris

Issue:
[BUG] Yaw optimization ignores provided `power_setpoints`
# Yaw optimization ignores provided power setpoints
As raised by @MiguelMarante in #1016, `YawOptimizerSR` (as well as other possibly other yaw optimization routines) do not correctly handle other controls being set, in particular, `disable_turbines` (which is a wrapper for near-zero `power_setpoints`).

The issue is that (several) local copies of the `FlorisModel` is made within the yaw optimizer, and the `power_setpoints` aren't copied over correctly. 

## How to reproduce
The following code demonstrates the problem (thanks to @MiguelMarante !)
```python
import numpy as np

from floris import FlorisModel
from floris.optimization.yaw_optimization.yaw_optimizer_sr import YawOptimizationSR

# Create floris object
yaml_file = 'inputs/emgauss.yaml'
fmodel = FlorisModel(yaml_file)

# Set wind conditions
fmodel.set(wind_directions=[270.]*2, wind_speeds=[8.]*2, turbulence_intensities=[.06]*2)

# Disable turbines
fmodel.set_operation_model('mixed')
fmodel.set(disable_turbines=np.array([[False, True, False], [True, False, False]]))

# Run optimization and print results
yaw_opt = YawOptimizationSR(fmodel)
df_opt = yaw_opt.optimize()
print(df_opt.yaw_angles_opt)
```

## Relevant output
Output:
```
0    [25.0, 22.65625, 0.0]
1    [25.0, 22.65625, 0.0]
Name: yaw_angles_opt, dtype: object
```
The expected output is not totally clear; the disabled turbine's angle is not well defined (since it produces no thrust and no power, any yaw angle should be optimal). However, we might expect the middle turbine to have a yaw angle of 25 degrees for the second case.

## Floris version
v4.1, v4.2

Discussion:
@misi9170: Addressed in #1031 (previously #1027)

Solution patch:
diff --git a/examples/examples_control_optimization/008_optimize_yaw_with_disabled_turbines.py b/examples/examples_control_optimization/008_optimize_yaw_with_disabled_turbines.py
new file mode 100644
index 000000000..bd4d80b7b
--- /dev/null
+++ b/examples/examples_control_optimization/008_optimize_yaw_with_disabled_turbines.py
@@ -0,0 +1,48 @@
+"""Example: Optimizing yaw angles with disabled turbines
+
+This example demonstrates how to optimize yaw angles in FLORIS, when some turbines are disabled.
+The example optimization is run using both YawOptimizerSR and YawOptimizerGeometric, the two
+yaw optimizers that support disabling turbines.
+"""
+
+import numpy as np
+
+from floris import FlorisModel
+from floris.optimization.yaw_optimization.yaw_optimizer_geometric import YawOptimizationGeometric
+from floris.optimization.yaw_optimization.yaw_optimizer_sr import YawOptimizationSR
+
+
+# Load a 3-turbine model
+fmodel = FlorisModel("../inputs/gch.yaml")
+
+# Set wind conditions to be the same for two cases
+fmodel.set(wind_directions=[270.]*2, wind_speeds=[8.]*2, turbulence_intensities=[.06]*2)
+
+# First run the case where all turbines are active and print results
+yaw_opt = YawOptimizationSR(fmodel)
+df_opt = yaw_opt.optimize()
+print("Serial Refine optimized yaw angles (all turbines active) [deg]:\n", df_opt.yaw_angles_opt)
+
+yaw_opt = YawOptimizationGeometric(fmodel)
+df_opt = yaw_opt.optimize()
+print("\nGeometric optimized yaw angles (all turbines active) [deg]:\n", df_opt.yaw_angles_opt)
+
+# Disable turbines (different pattern for each of the two cases)
+# First case: disable the middle turbine
+# Second case: disable the front turbine
+fmodel.set_operation_model('mixed')
+fmodel.set(disable_turbines=np.array([[False, True, False], [True, False, False]]))
+
+# Rerun optimizations and print results
+yaw_opt = YawOptimizationSR(fmodel)
+df_opt = yaw_opt.optimize()
+print(
+    "\nSerial Refine optimized yaw angles (some turbines disabled) [deg]:\n",
+    df_opt.yaw_angles_opt
+)
+# Note that disabled turbines are assigned a zero yaw angle, but their yaw angle is arbitrary as it
+# does not affect the total power output.
+
+yaw_opt = YawOptimizationGeometric(fmodel)
+df_opt = yaw_opt.optimize()
+print("\nGeometric optimized yaw angles (some turbines disabled) [deg]:\n", df_opt.yaw_angles_opt)
diff --git a/floris/core/turbine/operation_models.py b/floris/core/turbine/operation_models.py
index a6c1ff160..066e26b04 100644
--- a/floris/core/turbine/operation_models.py
+++ b/floris/core/turbine/operation_models.py
@@ -382,22 +382,22 @@ def axial_induction(
 @define
 class MixedOperationTurbine(BaseOperationModel):
 
+    @staticmethod
     def power(
         yaw_angles: NDArrayFloat,
         power_setpoints: NDArrayFloat,
         **kwargs
     ):
-        # Yaw angles mask all yaw_angles not equal to zero
-        yaw_angles_mask = yaw_angles != 0.0
-        power_setpoints_mask = power_setpoints < POWER_SETPOINT_DEFAULT
-        neither_mask = np.logical_not(yaw_angles_mask) & np.logical_not(power_setpoints_mask)
-
-        if (power_setpoints_mask & yaw_angles_mask).any():
-            raise ValueError((
-                "Power setpoints and yaw angles are incompatible."
-                "If yaw_angles entry is nonzero, power_setpoints must be greater than"
-                " or equal to {0}.".format(POWER_SETPOINT_DEFAULT)
-            ))
+        (
+            yaw_angles,
+            power_setpoints,
+            yaw_angles_mask,
+            power_setpoints_mask,
+            neither_mask
+        ) = MixedOperationTurbine._handle_mixed_operation_setpoints(
+            yaw_angles=yaw_angles,
+            power_setpoints=power_setpoints
+        )
 
         powers = np.zeros_like(power_setpoints)
         powers[yaw_angles_mask] += CosineLossTurbine.power(
@@ -414,21 +414,22 @@ def power(
 
         return powers
 
+    @staticmethod
     def thrust_coefficient(
         yaw_angles: NDArrayFloat,
         power_setpoints: NDArrayFloat,
         **kwargs
     ):
-        yaw_angles_mask = yaw_angles != 0.0
-        power_setpoints_mask = power_setpoints < POWER_SETPOINT_DEFAULT
-        neither_mask = np.logical_not(yaw_angles_mask) & np.logical_not(power_setpoints_mask)
-
-        if (power_setpoints_mask & yaw_angles_mask).any():
-            raise ValueError((
-                "Power setpoints and yaw angles are incompatible."
-                "If yaw_angles entry is nonzero, power_setpoints must be greater than"
-                " or equal to {0}.".format(POWER_SETPOINT_DEFAULT)
-            ))
+        (
+            yaw_angles,
+            power_setpoints,
+            yaw_angles_mask,
+            power_setpoints_mask,
+            neither_mask
+        ) = MixedOperationTurbine._handle_mixed_operation_setpoints(
+            yaw_angles=yaw_angles,
+            power_setpoints=power_setpoints
+        )
 
         thrust_coefficients = np.zeros_like(power_setpoints)
         thrust_coefficients[yaw_angles_mask] += CosineLossTurbine.thrust_coefficient(
@@ -445,21 +446,22 @@ def thrust_coefficient(
 
         return thrust_coefficients
 
+    @staticmethod
     def axial_induction(
         yaw_angles: NDArrayFloat,
         power_setpoints: NDArrayFloat,
         **kwargs
     ):
-        yaw_angles_mask = yaw_angles != 0.0
-        power_setpoints_mask = power_setpoints < POWER_SETPOINT_DEFAULT
-        neither_mask = np.logical_not(yaw_angles_mask) & np.logical_not(power_setpoints_mask)
-
-        if (power_setpoints_mask & yaw_angles_mask).any():
-            raise ValueError((
-                "Power setpoints and yaw angles are incompatible."
-                "If yaw_angles entry is nonzero, power_setpoints must be greater than"
-                " or equal to {0}.".format(POWER_SETPOINT_DEFAULT)
-            ))
+        (
+            yaw_angles,
+            power_setpoints,
+            yaw_angles_mask,
+            power_setpoints_mask,
+            neither_mask
+        ) = MixedOperationTurbine._handle_mixed_operation_setpoints(
+            yaw_angles=yaw_angles,
+            power_setpoints=power_setpoints
+        )
 
         axial_inductions = np.zeros_like(power_setpoints)
         axial_inductions[yaw_angles_mask] += CosineLossTurbine.axial_induction(
@@ -476,6 +478,34 @@ def axial_induction(
 
         return axial_inductions
 
+    @staticmethod
+    def _handle_mixed_operation_setpoints(
+        yaw_angles: NDArrayFloat,
+        power_setpoints: NDArrayFloat,
+    ):
+        """
+        Check for incompatible yaw angles and power setpoints and raise an error if found.
+        Return masks and updated setpoints.
+        """
+        # If any turbines are disabled, set their yaw angles to zero
+        yaw_angles[power_setpoints <= POWER_SETPOINT_DISABLED] = 0.0
+
+        # Create masks for whether yaw angles and power setpoints are set
+        yaw_angles_mask = yaw_angles != 0.0
+        power_setpoints_mask = power_setpoints < POWER_SETPOINT_DEFAULT
+        neither_mask = np.logical_not(yaw_angles_mask) & np.logical_not(power_setpoints_mask)
+
+        # Check for incompatibility and raise error if found.
+        if (power_setpoints_mask & yaw_angles_mask).any():
+            raise ValueError((
+                "Power setpoints and yaw angles are incompatible."
+                "If yaw_angles entry is nonzero, power_setpoints must be greater than"
+                " or equal to {0}.".format(POWER_SETPOINT_DEFAULT)
+            ))
+
+        # Return updated setpoints as well as masks
+        return yaw_angles, power_setpoints, yaw_angles_mask, power_setpoints_mask, neither_mask
+
 @define
 class AWCTurbine(BaseOperationModel):
     """
diff --git a/floris/floris_model.py b/floris/floris_model.py
index 09a5aa5d0..af8255232 100644
--- a/floris/floris_model.py
+++ b/floris/floris_model.py
@@ -453,10 +453,7 @@ def set(
         # previous setting
         if not (_yaw_angles == 0).all():
             self.core.farm.set_yaw_angles(_yaw_angles)
-        if not (
-            (_power_setpoints == POWER_SETPOINT_DEFAULT)
-            | (_power_setpoints == POWER_SETPOINT_DISABLED)
-        ).all():
+        if not (_power_setpoints == POWER_SETPOINT_DEFAULT).all():
             self.core.farm.set_power_setpoints(_power_setpoints)
         if _awc_modes is not None:
             self.core.farm.set_awc_modes(_awc_modes)
diff --git a/floris/optimization/yaw_optimization/yaw_optimization_base.py b/floris/optimization/yaw_optimization/yaw_optimization_base.py
index 07a2f7e11..1949d132c 100644
--- a/floris/optimization/yaw_optimization/yaw_optimization_base.py
+++ b/floris/optimization/yaw_optimization/yaw_optimization_base.py
@@ -5,6 +5,7 @@
 import numpy as np
 import pandas as pd
 
+from floris.core.turbine.operation_models import POWER_SETPOINT_DISABLED
 from floris.logging_manager import LoggingManager
 
 from .yaw_optimization_tools import derive_downstream_turbines
@@ -99,7 +100,7 @@ def __init__(
         """
 
         # Save turbine object to self
-        self.fmodel = fmodel.copy()
+        self.fmodel = copy.deepcopy(fmodel)
         self.nturbs = len(self.fmodel.layout_x)
 
         # # Check floris options
@@ -131,6 +132,11 @@ def __init__(
         self.minimum_yaw_angle = self._unpack_variable(minimum_yaw_angle)
         self.maximum_yaw_angle = self._unpack_variable(maximum_yaw_angle)
 
+        # Limit yaw angles to zero for disabled turbines
+        active_turbines = fmodel.core.farm.power_setpoints > POWER_SETPOINT_DISABLED
+        self.minimum_yaw_angle[~active_turbines] = 0.0
+        self.maximum_yaw_angle[~active_turbines] = 0.0
+
         # Set initial condition for optimization
         if x0 is not None:
             self.x0 = self._unpack_variable(x0)
@@ -224,7 +230,7 @@ def _reduce_control_problem(self):
         self.turbs_to_opt = (self.maximum_yaw_angle - self.minimum_yaw_angle >= 0.001)
 
         # Initialize subset variables as full set
-        self.fmodel_subset = self.fmodel.copy()
+        self.fmodel_subset = copy.deepcopy(self.fmodel)
         n_findex_subset = copy.deepcopy(self.fmodel.core.flow_field.n_findex)
         minimum_yaw_angle_subset = copy.deepcopy(self.minimum_yaw_angle)
         maximum_yaw_angle_subset = copy.deepcopy(self.maximum_yaw_angle)
@@ -301,6 +307,7 @@ def _calculate_farm_power(
             ti_array=None,
             turbine_weights=None,
             heterogeneous_speed_multipliers=None,
+            power_setpoints=None,
         ):
         """
         Calculate the wind farm power production assuming the predefined
@@ -353,6 +360,7 @@ def _calculate_farm_power(
             wind_speeds=ws_array,
             turbulence_intensities=ti_array,
             yaw_angles=yaw_angles,
+            power_setpoints=power_setpoints,
         )
         fmodel_subset.run()
         turbine_power = fmodel_subset.get_turbine_powers()
diff --git a/floris/optimization/yaw_optimization/yaw_optimizer_geometric.py b/floris/optimization/yaw_optimization/yaw_optimizer_geometric.py
index ea68204b4..68c687512 100644
--- a/floris/optimization/yaw_optimization/yaw_optimizer_geometric.py
+++ b/floris/optimization/yaw_optimization/yaw_optimizer_geometric.py
@@ -1,6 +1,7 @@
 
 import numpy as np
 
+from floris.core.turbine.operation_models import POWER_SETPOINT_DISABLED
 from floris.utilities import rotate_coordinates_rel_west
 
 from .yaw_optimization_base import YawOptimization
@@ -46,10 +47,11 @@ def optimize(self):
         # Loop through every WD individually. WS ignored!
         wd_array = self.fmodel_subset.core.flow_field.wind_directions
 
+        active_turbines = self.fmodel_subset.core.farm.power_setpoints > POWER_SETPOINT_DISABLED
         for nwdi, wd in enumerate(wd_array):
-            self._yaw_angles_opt_subset[nwdi, :] = geometric_yaw(
-                self.fmodel_subset.layout_x,
-                self.fmodel_subset.layout_y,
+            self._yaw_angles_opt_subset[nwdi, active_turbines[nwdi]] = geometric_yaw(
+                self.fmodel_subset.layout_x[active_turbines[nwdi]],
+                self.fmodel_subset.layout_y[active_turbines[nwdi]],
                 wd,
                 self.fmodel.core.farm.turbine_definitions[0]["rotor_diameter"],
                 top_left_yaw_upper=self.maximum_yaw_angle[0, 0],
diff --git a/floris/optimization/yaw_optimization/yaw_optimizer_scipy.py b/floris/optimization/yaw_optimization/yaw_optimizer_scipy.py
index cdde87656..810144c50 100644
--- a/floris/optimization/yaw_optimization/yaw_optimizer_scipy.py
+++ b/floris/optimization/yaw_optimization/yaw_optimizer_scipy.py
@@ -30,6 +30,12 @@ def __init__(
         Instantiate YawOptimizationScipy object with a FlorisModel object
         and assign parameter values.
         """
+        valid_op_models = ["cosine-loss"]
+        if fmodel.get_operation_model() not in valid_op_models:
+            raise ValueError(
+                "YawOptimizationScipy is currently limited to the following operation models: "
+                + ", ".join(valid_op_models)
+            )
         if opt_options is None:
             # Default SciPy parameters
             opt_options = {
diff --git a/floris/optimization/yaw_optimization/yaw_optimizer_sr.py b/floris/optimization/yaw_optimization/yaw_optimizer_sr.py
index 2b5b7ad1b..127cca01b 100644
--- a/floris/optimization/yaw_optimization/yaw_optimizer_sr.py
+++ b/floris/optimization/yaw_optimization/yaw_optimizer_sr.py
@@ -93,6 +93,7 @@ def _calc_powers_with_memory(self, yaw_angles_subset, use_memory=True):
         wd_array_subset = self.fmodel_subset.core.flow_field.wind_directions
         ws_array_subset = self.fmodel_subset.core.flow_field.wind_speeds
         ti_array_subset = self.fmodel_subset.core.flow_field.turbulence_intensities
+        power_setpoints_subset = self.fmodel_subset.core.farm.power_setpoints
         turbine_weights_subset = self._turbine_weights_subset
 
         # Reformat yaw_angles_subset, if necessary
@@ -108,6 +109,7 @@ def _calc_powers_with_memory(self, yaw_angles_subset, use_memory=True):
             wd_array_subset = np.tile(wd_array_subset, Ny)
             ws_array_subset = np.tile(ws_array_subset, Ny)
             ti_array_subset = np.tile(ti_array_subset, Ny)
+            power_setpoints_subset = np.tile(power_setpoints_subset, (Ny, 1))
             turbine_weights_subset = np.tile(turbine_weights_subset, (Ny, 1))
 
         # Initialize empty matrix for floris farm power outputs
@@ -143,7 +145,8 @@ def _calc_powers_with_memory(self, yaw_angles_subset, use_memory=True):
                 ti_array=ti_array_subset[~idx],
                 turbine_weights=turbine_weights_subset[~idx, :],
                 yaw_angles=yaw_angles_subset[~idx, :],
-                heterogeneous_speed_multipliers=het_sm
+                heterogeneous_speed_multipliers=het_sm,
+                power_setpoints=power_setpoints_subset[~idx, :],
             )
             self.time_spent_in_floris += (timerpc() - start_time)
 


Test patch:
diff --git a/tests/floris_model_integration_test.py b/tests/floris_model_integration_test.py
index fa05d43d3..a09a53ac1 100644
--- a/tests/floris_model_integration_test.py
+++ b/tests/floris_model_integration_test.py
@@ -10,7 +10,7 @@
     TimeSeries,
     WindRose,
 )
-from floris.core.turbine.operation_models import POWER_SETPOINT_DEFAULT
+from floris.core.turbine.operation_models import POWER_SETPOINT_DEFAULT, POWER_SETPOINT_DISABLED
 
 
 TEST_DATA = Path(__file__).resolve().parent / "data"
@@ -44,7 +44,24 @@ def test_assign_setpoints():
 
     # power_setpoints and disable_turbines (disable_turbines overrides power_setpoints)
     fmodel.set(power_setpoints=[[1e6, 2e6]], disable_turbines=[[True, False]])
-    assert np.allclose(fmodel.core.farm.power_setpoints, np.array([[0.001, 2e6]]))
+    assert np.allclose(fmodel.core.farm.power_setpoints, np.array([[POWER_SETPOINT_DISABLED, 2e6]]))
+
+    # Setting sequentially is equivalent to setting together
+    fmodel.reset_operation()
+    fmodel.set(disable_turbines=[[True, False]])
+    fmodel.set(yaw_angles=[[0, 30]])
+    assert np.allclose(
+        fmodel.core.farm.power_setpoints,
+        np.array([[POWER_SETPOINT_DISABLED, POWER_SETPOINT_DEFAULT]])
+    )
+    assert np.allclose(fmodel.core.farm.yaw_angles, np.array([[0, 30]]))
+
+    fmodel.set(disable_turbines=[[True, False]], yaw_angles=[[0, 30]])
+    assert np.allclose(
+        fmodel.core.farm.power_setpoints,
+        np.array([[POWER_SETPOINT_DISABLED, POWER_SETPOINT_DEFAULT]])
+    )
+    assert np.allclose(fmodel.core.farm.yaw_angles, np.array([[0, 30]]))
 
 def test_set_run():
     """
diff --git a/tests/geometric_yaw_unit_test.py b/tests/geometric_yaw_unit_test.py
new file mode 100644
index 000000000..2fb0452bb
--- /dev/null
+++ b/tests/geometric_yaw_unit_test.py
@@ -0,0 +1,115 @@
+
+import numpy as np
+import pandas as pd
+
+from floris import FlorisModel
+from floris.optimization.yaw_optimization.yaw_optimizer_geometric import YawOptimizationGeometric
+
+
+DEBUG = False
+VELOCITY_MODEL = "gauss"
+DEFLECTION_MODEL = "gauss"
+
+# Inputs for basic yaw optimizations
+WIND_DIRECTIONS = [0.0, 90.0, 180.0, 270.0]
+WIND_SPEEDS = [8.0] * 4
+TURBULENCE_INTENSITIES = [0.06] * 4
+LAYOUT_X = [0.0, 600.0, 1200.0]
+LAYOUT_Y = [0.0, 0.0, 0.0]
+MAXIMUM_YAW_ANGLE = 25.0
+
+def test_basic_optimization(sample_inputs_fixture):
+    """
+    The Serial Refine (SR) method optimizes yaw angles based on a sequential, iterative yaw
+    optimization scheme. This test checks basic properties of the optimization result.
+    """
+    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)
+
+    fmodel.set(
+        layout_x=LAYOUT_X,
+        layout_y=LAYOUT_Y,
+        wind_directions=WIND_DIRECTIONS,
+        wind_speeds=WIND_SPEEDS,
+        turbulence_intensities=TURBULENCE_INTENSITIES
+    )
+    fmodel.set_operation_model("cosine-loss")
+
+    yaw_opt = YawOptimizationGeometric(
+        fmodel,
+        minimum_yaw_angle=0.0,
+        maximum_yaw_angle=MAXIMUM_YAW_ANGLE
+    )
+    df_opt = yaw_opt.optimize()
+
+    # Unaligned conditions
+    assert np.allclose(df_opt.loc[0, "yaw_angles_opt"], 0.0)
+    assert np.allclose(df_opt.loc[2, "yaw_angles_opt"], 0.0)
+
+    # Check aligned conditions
+    # Check maximum and minimum are respected
+    assert (df_opt.loc[1, "yaw_angles_opt"] <= MAXIMUM_YAW_ANGLE).all()
+    assert (df_opt.loc[3, "yaw_angles_opt"] <= MAXIMUM_YAW_ANGLE).all()
+    assert (df_opt.loc[1, "yaw_angles_opt"] >= 0.0).all()
+    assert (df_opt.loc[3, "yaw_angles_opt"] >= 0.0).all()
+
+    # Check 90.0 and 270.0 are symmetric
+    assert np.allclose(df_opt.loc[1, "yaw_angles_opt"], np.flip(df_opt.loc[3, "yaw_angles_opt"]))
+
+    # Check last turbine's angles are zero at 270.0
+    assert np.allclose(df_opt.loc[3, "yaw_angles_opt"][-1], 0.0)
+
+    # YawOptimizationGeometric does not compute farm powers
+
+def test_disabled_turbines(sample_inputs_fixture):
+    """
+    Tests SR when some turbines are disabled and checks that the results are equivalent to removing
+    those turbines from the wind farm. Need a tight layout to ensure that the front-to-back distance
+    is not too large.
+    """
+
+    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)
+
+    fmodel.set(
+        layout_x=LAYOUT_X,
+        layout_y=LAYOUT_Y,
+        wind_directions=WIND_DIRECTIONS,
+        wind_speeds=WIND_SPEEDS,
+        turbulence_intensities=TURBULENCE_INTENSITIES
+    )
+    fmodel.set_operation_model("mixed")
+
+    # Disable the middle turbine in all wind conditions, run optimization, and extract results
+    fmodel.set(disable_turbines=[[False, True, False]]*4)
+    yaw_opt = YawOptimizationGeometric(
+        fmodel,
+        minimum_yaw_angle=0.0,
+        maximum_yaw_angle=MAXIMUM_YAW_ANGLE
+    )
+    df_opt = yaw_opt.optimize()
+    yaw_angles_opt_disabled = df_opt.loc[3, "yaw_angles_opt"]
+
+    # Set up a new wind farm with the middle turbine removed
+    fmodel = FlorisModel(sample_inputs_fixture.core)
+    fmodel.set(
+        layout_x=np.array(LAYOUT_X)[[0, 2]],
+        layout_y=np.array(LAYOUT_Y)[[0, 2]],
+        wind_directions=WIND_DIRECTIONS,
+        wind_speeds=WIND_SPEEDS,
+        turbulence_intensities=TURBULENCE_INTENSITIES
+    )
+    fmodel.set_operation_model("cosine-loss")
+    yaw_opt = YawOptimizationGeometric(
+        fmodel,
+        minimum_yaw_angle=0.0,
+        maximum_yaw_angle=MAXIMUM_YAW_ANGLE
+    )
+    df_opt = yaw_opt.optimize()
+    yaw_angles_opt_removed = df_opt.loc[3, "yaw_angles_opt"]
+
+    assert np.allclose(yaw_angles_opt_disabled[[0, 2]], yaw_angles_opt_removed)
diff --git a/tests/serial_refine_unit_test.py b/tests/serial_refine_unit_test.py
new file mode 100644
index 000000000..cfda030a7
--- /dev/null
+++ b/tests/serial_refine_unit_test.py
@@ -0,0 +1,107 @@
+
+import numpy as np
+import pandas as pd
+
+from floris import FlorisModel
+from floris.optimization.yaw_optimization.yaw_optimizer_sr import YawOptimizationSR
+
+
+DEBUG = False
+VELOCITY_MODEL = "gauss"
+DEFLECTION_MODEL = "gauss"
+
+# Inputs for basic yaw optimizations
+WIND_DIRECTIONS = [0.0, 90.0, 180.0, 270.0]
+WIND_SPEEDS = [8.0] * 4
+TURBULENCE_INTENSITIES = [0.06] * 4
+LAYOUT_X = [0.0, 600.0, 1200.0]
+LAYOUT_Y = [0.0, 0.0, 0.0]
+MAXIMUM_YAW_ANGLE = 25.0
+
+def test_basic_optimization(sample_inputs_fixture):
+    """
+    The Serial Refine (SR) method optimizes yaw angles based on a sequential, iterative yaw
+    optimization scheme. This test checks basic properties of the optimization result.
+    """
+    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)
+
+    fmodel.set(
+        layout_x=LAYOUT_X,
+        layout_y=LAYOUT_Y,
+        wind_directions=WIND_DIRECTIONS,
+        wind_speeds=WIND_SPEEDS,
+        turbulence_intensities=TURBULENCE_INTENSITIES
+    )
+    fmodel.set_operation_model("cosine-loss")
+
+    yaw_opt = YawOptimizationSR(fmodel, minimum_yaw_angle=0.0, maximum_yaw_angle=MAXIMUM_YAW_ANGLE)
+    df_opt = yaw_opt.optimize()
+
+    # Unaligned conditions
+    assert np.allclose(df_opt.loc[0, "yaw_angles_opt"], 0.0)
+    assert np.allclose(df_opt.loc[2, "yaw_angles_opt"], 0.0)
+
+    # Check aligned conditions
+    # Check maximum and minimum are respected
+    assert (df_opt.loc[1, "yaw_angles_opt"] <= MAXIMUM_YAW_ANGLE).all()
+    assert (df_opt.loc[3, "yaw_angles_opt"] <= MAXIMUM_YAW_ANGLE).all()
+    assert (df_opt.loc[1, "yaw_angles_opt"] >= 0.0).all()
+    assert (df_opt.loc[3, "yaw_angles_opt"] >= 0.0).all()
+
+    # Check 90.0 and 270.0 are symmetric
+    assert np.allclose(df_opt.loc[1, "yaw_angles_opt"], np.flip(df_opt.loc[3, "yaw_angles_opt"]))
+
+    # Check last turbine's angles are zero at 270.0
+    assert np.allclose(df_opt.loc[3, "yaw_angles_opt"][-1], 0.0)
+
+    # Check that optimizer reports a power improvement
+    assert (df_opt["farm_power_opt"] >= df_opt["farm_power_baseline"]).all()
+
+def test_disabled_turbines(sample_inputs_fixture):
+    """
+    Tests SR when some turbines are disabled and checks that the results are equivalent to removing
+    those turbines from the wind farm. Need a tight layout to ensure that the front-to-back distance
+    is not too large.
+    """
+
+    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)
+
+    fmodel.set(
+        layout_x=LAYOUT_X,
+        layout_y=LAYOUT_Y,
+        wind_directions=WIND_DIRECTIONS,
+        wind_speeds=WIND_SPEEDS,
+        turbulence_intensities=TURBULENCE_INTENSITIES
+    )
+    fmodel.set_operation_model("mixed")
+
+    # Disable the middle turbine in all wind conditions, run optimization, and extract results
+    fmodel.set(disable_turbines=[[False, True, False]]*4)
+    yaw_opt = YawOptimizationSR(fmodel, minimum_yaw_angle=0.0, maximum_yaw_angle=MAXIMUM_YAW_ANGLE)
+    df_opt = yaw_opt.optimize()
+    yaw_angles_opt_disabled = df_opt.loc[3, "yaw_angles_opt"]
+    farm_power_opt_disabled = df_opt.loc[3, "farm_power_opt"]
+
+    # Set up a new wind farm with the middle turbine removed
+    fmodel = FlorisModel(sample_inputs_fixture.core)
+    fmodel.set(
+        layout_x=np.array(LAYOUT_X)[[0, 2]],
+        layout_y=np.array(LAYOUT_Y)[[0, 2]],
+        wind_directions=WIND_DIRECTIONS,
+        wind_speeds=WIND_SPEEDS,
+        turbulence_intensities=TURBULENCE_INTENSITIES
+    )
+    fmodel.set_operation_model("cosine-loss")
+    yaw_opt = YawOptimizationSR(fmodel, minimum_yaw_angle=0.0, maximum_yaw_angle=MAXIMUM_YAW_ANGLE)
+    df_opt = yaw_opt.optimize()
+    yaw_angles_opt_removed = df_opt.loc[3, "yaw_angles_opt"]
+    farm_power_opt_removed = df_opt.loc[3, "farm_power_opt"]
+
+    assert np.allclose(yaw_angles_opt_disabled[[0, 2]], yaw_angles_opt_removed)
+    assert np.allclose(farm_power_opt_disabled, farm_power_opt_removed)