diff --git a/floris/simulation/turbine.py b/floris/simulation/turbine.py
index b1cef9858..b284b10cd 100644
--- a/floris/simulation/turbine.py
+++ b/floris/simulation/turbine.py
@@ -665,7 +665,9 @@ def __attrs_post_init__(self) -> None:
         )
         self.power_interp = interp1d(
             wind_speeds,
-            inner_power
+            inner_power,
+            bounds_error=False,
+            fill_value=0
         )
 
         """
diff --git a/floris/tools/floris_interface.py b/floris/tools/floris_interface.py
index 4e4cc352e..a09003ec6 100644
--- a/floris/tools/floris_interface.py
+++ b/floris/tools/floris_interface.py
@@ -595,6 +595,10 @@ def get_turbine_powers(self) -> NDArrayFloat:
                 "Can't run function `FlorisInterface.get_turbine_powers` without "
                 "first running `FlorisInterface.calculate_wake`."
             )
+        # Check for negative velocities, which could indicate bad model
+        # parameters or turbines very closely spaced.
+        if (self.turbine_effective_velocities < 0.).any():
+            self.logger.warning("Some rotor effective velocities are negative.")
 
         turbine_powers = power(
             ref_density_cp_ct=self.floris.farm.ref_density_cp_cts,
diff --git a/floris/tools/optimization/yaw_optimization/yaw_optimization_base.py b/floris/tools/optimization/yaw_optimization/yaw_optimization_base.py
index ddddff55a..6a50247d9 100644
--- a/floris/tools/optimization/yaw_optimization/yaw_optimization_base.py
+++ b/floris/tools/optimization/yaw_optimization/yaw_optimization_base.py
@@ -338,7 +338,9 @@ def _normalize_control_problem(self):
             / self._normalization_length
         )
 
-    def _calculate_farm_power(self, yaw_angles=None, wd_array=None, turbine_weights=None):
+    def _calculate_farm_power(self, yaw_angles=None, wd_array=None, turbine_weights=None,
+            heterogeneous_speed_multipliers=None
+        ):
         """
         Calculate the wind farm power production assuming the predefined
         probability distribution (self.unc_options/unc_pmf), with the
@@ -358,6 +360,9 @@ def _calculate_farm_power(self, yaw_angles=None, wd_array=None, turbine_weights=
             yaw_angles = self._yaw_angles_baseline_subset
         if turbine_weights is None:
             turbine_weights = self._turbine_weights_subset
+        if heterogeneous_speed_multipliers is not None:
+            fi_subset.floris.flow_field.\
+                heterogenous_inflow_config['speed_multipliers'] = heterogeneous_speed_multipliers
 
         # Ensure format [incompatible with _subset notation]
         yaw_angles = self._unpack_variable(yaw_angles, subset=True)
diff --git a/floris/tools/optimization/yaw_optimization/yaw_optimizer_scipy.py b/floris/tools/optimization/yaw_optimization/yaw_optimizer_scipy.py
index 81ffc32a3..66339e426 100644
--- a/floris/tools/optimization/yaw_optimization/yaw_optimizer_scipy.py
+++ b/floris/tools/optimization/yaw_optimization/yaw_optimizer_scipy.py
@@ -108,6 +108,16 @@ def optimize(self):
                 yaw_template = np.tile(yaw_template, (1, 1, 1))
                 turbine_weights = np.tile(turbine_weights, (1, 1, 1))
 
+                # Handle heterogeneous inflow, if there is one
+                if (hasattr(self.fi.floris.flow_field, 'heterogenous_inflow_config') and
+                    self.fi.floris.flow_field.heterogenous_inflow_config is not None):
+                    het_sm_orig = np.array(
+                        self.fi.floris.flow_field.heterogenous_inflow_config['speed_multipliers']
+                    )
+                    het_sm = het_sm_orig[nwdi,:].reshape(1,-1)
+                else:
+                    het_sm = None
+
                 # Define cost function
                 def cost(x):
                     x_full = np.array(yaw_template, copy=True)
@@ -116,7 +126,8 @@ def cost(x):
                         - 1.0 * self._calculate_farm_power(
                             yaw_angles=x_full,
                             wd_array=[wd],
-                            turbine_weights=turbine_weights
+                            turbine_weights=turbine_weights,
+                            heterogeneous_speed_multipliers=het_sm
                         )[0, 0] / J0
                     )
 
diff --git a/floris/tools/optimization/yaw_optimization/yaw_optimizer_sr.py b/floris/tools/optimization/yaw_optimization/yaw_optimizer_sr.py
index 801c59312..82d50ef08 100644
--- a/floris/tools/optimization/yaw_optimization/yaw_optimizer_sr.py
+++ b/floris/tools/optimization/yaw_optimization/yaw_optimizer_sr.py
@@ -141,10 +141,19 @@ def _calc_powers_with_memory(self, yaw_angles_subset, use_memory=True):
         if not np.all(idx):
             # Now calculate farm powers for conditions we haven't yet evaluated previously
             start_time = timerpc()
+            if (hasattr(self.fi.floris.flow_field, 'heterogenous_inflow_config') and
+                self.fi.floris.flow_field.heterogenous_inflow_config is not None):
+                het_sm_orig = np.array(
+                    self.fi.floris.flow_field.heterogenous_inflow_config['speed_multipliers']
+                )
+                het_sm = np.tile(het_sm_orig, (Ny, 1))[~idx, :]
+            else:
+                het_sm = None
             farm_powers[~idx, :] = self._calculate_farm_power(
                 wd_array=wd_array_subset[~idx],
                 turbine_weights=turbine_weights_subset[~idx, :, :],
                 yaw_angles=yaw_angles_subset[~idx, :, :],
+                heterogeneous_speed_multipliers=het_sm
             )
             self.time_spent_in_floris += (timerpc() - start_time)