Repository: NREL/floris

Issue:
Bug report: `WindRose.plot()` argument `legend_kwargs` are said to be passed to `legend` in the doc strings, but they are passed to `colorbar` instead
<!--
IMPORTANT NOTES

Thank you for taking the time to report a bug. If you aren't certain whether an issue
is a bug, please first open a Discussion. Before submitting, please reread your
description to ensure that other readers can reasonably understand the issue
you're facing and the impact on your workflow or results.

This form is written in GitHub's Markdown format. For a reference on this type
of syntax, see GitHub's documentation:
https://docs.github.com/en/get-started/writing-on-github/getting-started-with-writing-and-formatting-on-github/basic-writing-and-formatting-syntax

This template contains guidance for your submission within the < ! - -, - - > blocks.
These are comments in HTML syntax and will not appear in the submission.
Be sure to use the "Preview" feature on GitHub to ensure your submission is formatted as intended.

When including code snippets, please paste the text itself and wrap the code block with
ticks (see the other character on the tilde ~ key in a US keyboard) to format it as code.
For example, Python code should be wrapped in ticks like this:
```python
def a_func():
    return 1

a = 1
b = a_func()
print(a + b)
```
This is preferred over screen shots since it is searchable and others can copy/paste
the text to run it.
-->

# Bug report: `WindRose.plot()` argument `legend_kwargs` are said to be passed to `legend` in the doc strings, but they are passed to `colorbar` instead
<!--
A clear and concise description of the bug including what happened
and what you expected to happen.
-->
`WindRose.plot()` argument `legend_kwargs` are said to be passed to `legend` in the doc strings, but they are passed to `colorbar` instead. This makes usage non-intuitive.

## How to reproduce
<!--
Describe how another person with no context can recreate this issue.
It is typically very helpful to reduce the problem as much as possible
and share a minimum working example. See the note above on including
code as text rather than screenshots.
-->
Pass typical legend arguments to `WindRose.plot()`
```python

    import matplotlib.pyplot as plt

    wind_data = fmodel.wind_data
  
    time_series = TimeSeries(
        wind_directions=wind_data.wind_directions,
        wind_speeds=wind_data.wind_speeds,
        turbulence_intensities=wind_data.turbulence_intensities,
    )
          
    wind_rose = time_series.to_WindRose(wd_edges=np.arange(0, 360, 3.0), ws_edges=np.arange(2, 20, 2.0))
      
    fig, ax = plt.subplots(1, figsize=(4,4), subplot_kw={'projection': 'polar'})

    ax = wind_rose.plot(ax=ax, wd_step=10, legend_kwargs={"title": "Wind Speed (m/s)", "frameon": True})
```

## Relevant output
<!--
Include any output, plots, or other means of communication here to add context to the problem.
 -->
```text
  File "floris/floris/wind_data.py", line 782, in plot
    ax.figure.colorbar(sm_ws, ax=ax, **legend_kwargs)
  File "matplotlib/figure.py", line 1253, in colorbar
    cb = cbar.Colorbar(cax, mappable, **{
         ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
```

## Floris version
<!--
Share your floris version and how you installed it. You can print the floris version from
a Python REPL or script with these commands:
```python
import floris
print(floris.__version__)
```
 -->
4.2

## System Information
<!-- Add your information here. -->
 - OS: macOS 14.7.1
 - Python version: 3.11
 - Library versions
   - matplotlib: 3.9.2

Discussion:
@misi9170: Hi @jaredthomas68 , thanks for raising this. I believe the problem was introduced in #969. I think you can get around this for now by setting `legend_kwargs` to `{"label": Wind Speed (m/s)"}` (and dropping the `frameon` argument). But this is a user-interface change that we should have handled more thoroughly in a minor release, so I'll look into a more robust bugfix.@misi9170: Hi @jaredthomas68 , I've believe I've now addressed this in #1028 

Solution patch:
diff --git a/.github/workflows/python-publish.yml b/.github/workflows/python-publish.yml
index ead77afce..9281f86c0 100644
--- a/.github/workflows/python-publish.yml
+++ b/.github/workflows/python-publish.yml
@@ -20,11 +20,11 @@ jobs:
     - name: Install dependencies
       run: |
         python -m pip install --upgrade pip
-        pip install setuptools wheel twine
+        pip install build twine
     - name: Build and publish
       env:
         TWINE_USERNAME: ${{ secrets.PYPI_USERNAME }}
         TWINE_PASSWORD: ${{ secrets.PYPI_PASSWORD }}
       run: |
-        python setup.py sdist bdist_wheel
+        python -m build
         twine upload dist/*
diff --git a/docs/dev_guide.md b/docs/dev_guide.md
index 4a7c5dd76..a4a08e195 100644
--- a/docs/dev_guide.md
+++ b/docs/dev_guide.md
@@ -210,7 +210,7 @@ is located at `floris/.github/workflows/continuous-integration-workflow.yaml`.
 The online documentation is built with Jupyter Book which uses Sphinx
 as a framework. It is automatically built and hosted by GitHub, but it
 can also be compiled locally. Additional dependencies are required
-for the documentation, and they are listed in the `EXTRAS` of `setup.py`.
+for the documentation, and they are listed in the `project.optional-dependencies` of `pyproject.toml`.
 The commands to build the docs are given below. After successfully
 compiling, a file should be located at ``docs/_build/html/index.html``.
 This file can be opened in any browser.
@@ -246,7 +246,7 @@ Be sure to complete each step in the sequence as described.
     with a commit message such as "Update version to vN.M".
     The version number must be updated in the following two files:
     - [floris/README.md](https://github.com/NREL/floris/blob/main/README.md)
-    - [floris/floris/version.py](https://github.com/NREL/floris/blob/main/floris/version.py)
+    - [pyproject.toml](https://github.com/NREL/floris/blob/main/pyproject.toml)
     Note that a `.0` version number is left off meaning that valid versions
     are `v3`, `v3.1`, `v3.1.1`, etc.
 
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/__init__.py b/floris/__init__.py
index d97bb24eb..717ea3175 100644
--- a/floris/__init__.py
+++ b/floris/__init__.py
@@ -1,9 +1,9 @@
 
+from importlib.metadata import version
 from pathlib import Path
 
 
-with open(Path(__file__).parent / "version.py") as _version_file:
-    __version__ = _version_file.read().strip()
+__version__ = version("floris")
 
 
 from .floris_model import FlorisModel
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 1f92fa092..83610e993 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)
 
diff --git a/floris/version.py b/floris/version.py
deleted file mode 100644
index bf77d5496..000000000
--- a/floris/version.py
+++ /dev/null
@@ -1 +0,0 @@
-4.2
diff --git a/floris/wind_data.py b/floris/wind_data.py
index b470fe515..3d49dec4b 100644
--- a/floris/wind_data.py
+++ b/floris/wind_data.py
@@ -779,8 +779,12 @@ def plot(
                 )
 
         # Configure the plot
-        ax.figure.colorbar(sm_ws, ax=ax, **legend_kwargs)
-        ax.figure.tight_layout()
+        try:
+           ax.figure.colorbar(sm_ws, ax=ax, **legend_kwargs)
+           ax.figure.tight_layout()
+        except TypeError:
+           ax.legend(reversed(rects), ws_bins, **legend_kwargs)
+           ax.figure.get_children()[-1].remove() # Remove the empty colorbar
         ax.set_theta_direction(-1)
         ax.set_theta_offset(np.pi / 2.0)
         ax.set_theta_zero_location("N")
@@ -1822,8 +1826,12 @@ def plot(
                 )
 
         # Configure the plot
-        ax.figure.colorbar(sm_wv, ax=ax, **legend_kwargs)
-        ax.figure.tight_layout()
+        try:
+            ax.figure.colorbar(sm_wv, ax=ax, **legend_kwargs)
+            ax.figure.tight_layout()
+        except TypeError:
+            ax.legend(reversed(rects), var_bins, **legend_kwargs)
+            ax.figure.get_children()[-1].remove() # Remove the empty colorbar
         ax.set_theta_direction(-1)
         ax.set_theta_offset(np.pi / 2.0)
         ax.set_theta_zero_location("N")
diff --git a/pyproject.toml b/pyproject.toml
index 330c5a2d8..dea224866 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -2,6 +2,71 @@
 requires = ["setuptools >= 40.6.0", "wheel"]
 build-backend = "setuptools.build_meta"
 
+[project]
+name = "floris"
+version = "4.2"
+description = "A controls-oriented engineering wake model."
+readme = "README.md"
+requires-python = ">=3.8"
+authors = [
+    { name = "Rafael Mudafort", email = "rafael.mudafort@nrel.gov" },
+    { name = "Paul Fleming", email = "paul.fleming@nrel.gov" },
+    { name = "Michael (Misha) Sinner", email = "Michael.Sinner@nrel.gov" },
+    { name = "Eric Simley", email = "Eric.Simley@nrel.gov" },
+    { name = "Christopher Bay", email = "Christopher.Bay@nrel.gov" },
+]
+license = { file = "LICENSE.txt" }
+keywords = ["floris"]
+classifiers = [
+        "License :: OSI Approved :: BSD License",
+        "Programming Language :: Python",
+        "Programming Language :: Python :: 3",
+        "Programming Language :: Python :: 3.8",
+        "Programming Language :: Python :: 3.9",
+        "Programming Language :: Python :: 3.10",
+        "Programming Language :: Python :: Implementation :: CPython",
+        "Programming Language :: Python :: Implementation :: PyPy"
+]
+dependencies = [
+    "attrs",
+    "pyyaml~=6.0",
+    "numexpr~=2.0",
+    "numpy~=1.20",
+    "scipy~=1.1",
+    "matplotlib~=3.0",
+    "pandas~=2.0",
+    "shapely~=2.0",
+    "coloredlogs~=15.0",
+    "pathos~=0.3",
+]
+
+[project.optional-dependencies]
+docs = [
+    "jupyter-book",
+    "sphinx-book-theme",
+    "sphinx-autodoc-typehints",
+    "sphinxcontrib-autoyaml",
+    "sphinxcontrib.mermaid",
+]
+develop = [
+    "pytest",
+    "pre-commit",
+    "ruff",
+    "isort"
+]
+
+[tool.setuptools.packages.find]
+include = ["floris*"]
+
+[tool.setuptools.package-data]
+floris = [
+    "turbine_library/*.yaml",
+    "core/wake_velocity/turbopark_lookup_table.mat"
+]
+
+[project.urls]
+Homepage = "https://github.com/NREL/floris"
+Documentation = "https://nrel.github.io/floris/"
 
 [coverage.run]
 # Coverage.py configuration file
diff --git a/setup.py b/setup.py
deleted file mode 100644
index 7e3a11be4..000000000
--- a/setup.py
+++ /dev/null
@@ -1,89 +0,0 @@
-
-from pathlib import Path
-
-from setuptools import find_packages, setup
-
-
-# Package meta-data.
-NAME = "FLORIS"
-DESCRIPTION = "A controls-oriented engineering wake model."
-URL = "https://github.com/NREL/FLORIS"
-EMAIL = "rafael.mudafort@nrel.gov"
-AUTHOR = "NREL National Wind Technology Center"
-REQUIRES_PYTHON = ">=3.8.0"
-
-# What packages are required for this module to be executed?
-REQUIRED = [
-    # simulation
-    "attrs",
-    "pyyaml~=6.0",
-    "numexpr~=2.0",
-    "numpy~=1.20",
-    "scipy~=1.1",
-
-    # tools
-    "matplotlib~=3.0",
-    "pandas~=2.0",
-    "shapely~=2.0",
-
-    # utilities
-    "coloredlogs~=15.0",
-    "pathos~=0.3",
-]
-
-# What packages are optional?
-# To use:
-#   pip install -e ".[docs,develop]"    install both sets of extras in editable install
-#   pip install -e ".[develop]"         installs only developer packages in editable install
-#   pip install "floris[develop]"       installs developer packages in non-editable install
-EXTRAS = {
-    "docs": {
-        "jupyter-book",
-        "sphinx-book-theme",
-        "sphinx-autodoc-typehints",
-        "sphinxcontrib-autoyaml",
-        "sphinxcontrib.mermaid",
-    },
-    "develop": {
-        "pytest",
-        "pre-commit",
-        "ruff",
-        "isort",
-    },
-}
-
-ROOT = Path(__file__).parent
-with open(ROOT / "floris" / "version.py") as version_file:
-    VERSION = version_file.read().strip()
-
-setup(
-    name=NAME,
-    version=VERSION,
-    description=DESCRIPTION,
-    long_description=DESCRIPTION,
-    long_description_content_type="text/markdown",
-    author=AUTHOR,
-    author_email=EMAIL,
-    python_requires=REQUIRES_PYTHON,
-    url=URL,
-    packages=find_packages(exclude=["tests", "*.tests", "*.tests.*", "tests.*"]),
-    package_data={
-        'floris': ['turbine_library/*.yaml', 'core/wake_velocity/turbopark_lookup_table.mat']
-    },
-    install_requires=REQUIRED,
-    extras_require=EXTRAS,
-    include_package_data=True,
-    license_files = ('LICENSE.txt',),
-    classifiers=[
-        # Trove classifiers
-        # Full list: https://pypi.python.org/pypi?%3Aaction=list_classifiers
-        "License :: OSI Approved :: BSD License",
-        "Programming Language :: Python",
-        "Programming Language :: Python :: 3",
-        "Programming Language :: Python :: 3.8",
-        "Programming Language :: Python :: 3.9",
-        "Programming Language :: Python :: 3.10",
-        "Programming Language :: Python :: Implementation :: CPython",
-        "Programming Language :: Python :: Implementation :: PyPy"
-    ],
-)


Test patch:
diff --git a/tests/floris_model_integration_test.py b/tests/floris_model_integration_test.py
index 395ff03a1..88a242c56 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)