Repository: NREL/floris

Issue:
Bug in AWC operation model

# Bug in AWC operation model

There appears to be a bug in the 'awc' operational model in floris concerning the computation of power and thrust of actuated turbines (turbines set to 'helix' as opposed to 'baseline').

The error stems from the use of any() here: https://github.com/NREL/floris/blob/c50389cfba596e29dcbdfa7670db591b2b65da3e/floris/core/turbine/operation_models.py#L540 
Instead of only reducing the power of 'helix' turbines, if any turbine is labelled as 'helix' in awc_modes, the power of all turbines will be reduced. 


## How to reproduce
Consider the following minimal example: 
```
import numpy as np 
from floris import (
    FlorisModel,
    WindRose,
)
cases = ["Baseline","AWC",]
fmodels = {}
fmodels['Baseline'] = FlorisModel("../floris/examples/inputs/emgauss_helix.yaml")
fmodels['AWC'] = FlorisModel("../floris/examples/inputs/emgauss_helix.yaml")
awc_array = ['helix','baseline','baseline']
baseline_array = ['baseline','baseline','baseline']
awc_amplitudes = [4,4,4]

fmodels['Baseline'].set_operation_model("awc")
fmodels['Baseline'].set(
    awc_modes=np.array([baseline_array]),
    awc_amplitudes=np.array([awc_amplitudes])
    )

fmodels['AWC'].set_operation_model("awc")
fmodels['AWC'].set(
    awc_modes=np.array([awc_array]),
    awc_amplitudes=np.array([awc_amplitudes])
    )

turbine_powers = {}
for case in cases:
    floris_model = fmodels[case]
    floris_model.run_no_wake()
    turbine_powers[case] = floris_model.get_turbine_powers()/1000
    print(case,": ",turbine_powers[case])
```
which results it: 
`Baseline :  [[6301.1399872 6301.1399872 6301.1399872]]`
`AWC :  [[5908.03837044 5908.03837044 5908.03837044]]`
wheras the correct answer should be: 
`Baseline :  [[6301.1399872 6301.1399872 6301.1399872]]`
`AWC :  [[5908.03837044 6301.1399872  6301.1399872 ]]`


Additionally, this will cause an error when performing a wind sweep or estimating AEP. If 'helix' occurs for any turbine at any wind condition, then all turbines at all wind conditions see a reduction in power. Consider the following minimal example with two wind conditions:  
```
for caseiter, case in enumerate(cases):
    fmodel = fmodels[case]
    fmodel.set(
    wind_directions=np.array([270.0, 270.0]),
    wind_speeds=[5.0, 12.0],
    turbulence_intensities=np.array([0.06, 0.06]),
    )
    baseline_setting =  np.array([['baseline', 'baseline', 'baseline'],['baseline', 'baseline', 'baseline']])
    awc_setting =  np.array([['helix', 'baseline', 'baseline'],['baseline', 'baseline', 'baseline']])
    awc_amps = np.array([awc_amplitudes]*2)
    if 'Baseline' in case:
        fmodel.set(awc_modes=baseline_setting,awc_amplitudes=awc_amps)
    else:
        fmodel.set(awc_modes=awc_setting,awc_amplitudes=awc_amps)
    fmodel.run()

turbine_powers = {}
for case in cases:
    floris_model = fmodels[case]
    floris_model.run_no_wake()
    turbine_powers[case] = floris_model.get_turbine_powers()/1000
    print(case,": ",turbine_powers[case])
```
This will actually raise an incorrect warning and abort because of the message starting here: 
https://github.com/NREL/floris/blob/c50389cfba596e29dcbdfa7670db591b2b65da3e/floris/core/turbine/operation_models.py#L545
Again, this is a results of the bug. 'helix' should only be applied for WS=5, which is below rated, but because of the .any() call, there is both a 'helix' turbine and a wind speed above rated (WS=12) so the warning will appear. Even if you comment out this warning, the code produces:
`Baseline :  [[ 1375.2462449  1375.2462449  1375.2462449]
 [15000.        15000.        15000.       ]]`
`AWC :  [[ 1292.7915507  1292.7915507  1292.7915507]
 [13999.8581506 13999.8581506 13999.8581506]]`
whereas the correct result should be: 
`Baseline :  [[ 1375.2462449  1375.2462449  1375.2462449]
 [15000.        15000.        15000.       ]]`
`AWC :  [[ 1292.7915507  1375.2462449  1375.2462449]
 [15000.        15000.        15000.       ]] `
 
One way around this issue is to assume that awc_amplitudes applies even for 'baseline' turbines. However, this is confusing and error prone, because if all turbines are 'baseline' then awc_amplitudes doesn't matter but if one turbine at one wind condition is 'helix' then awc_amplitude does matter for all 'baseline' turbines at all wind conditions. Moreover, just relying on awc_amplitudes will not fix the issue with the warning mentioned above. A general fix is on my branch of floris here:  https://github.com/gyalla/floris/tree/awc_fix, and I'm happy to open a pull request.

## Floris version
Originally found with version 4.2.2 installed with conda

Hints from issue comments:
@paulf81: @jfrederik-nrel  would you mind to take a look at this one?@jfrederik-nrel: Hi @gyalla,. Thank you for bringing this bug to our attention. I pulled your branch and the solution looks good in general. However, I did find a bug in your code when you run multiple cases at once, i.e., when `awc_modes` has different dimensions than `awc_amplitudes` (as is the case in examples_control_types example 4). 

I fixed this by adding the following two lines in `operation_models.py` after line 555, where you define `mask`:
`if mask.shape[0] != base_powers.shape[0]:`
`     mask = np.tile(mask, (base_powers.shape[0],1))`

I verified that this fixes this issue, but I think this is probably not the most elegant solution. 

@misi9170, @paulf81, I'm pretty sure FLORIS has a way of handling inputs of different dimensions in `fmodel.set()` in other functions, does it not? Would either one of you have a better suggestion of how to deal with this in a more general fashion?@misi9170: HI @gyalla , agreed with @jfrederik-nrel --- thanks for identifying this issue and suggesting a bugfix. Yes, please go ahead and open up a pull request to get your fix in; we may make some modifications to your branch once the PR is open, if that's ok. 

@jfrederik-nrel , interesting. Generally speaking, no, we expect the control setpoints should always be of dimension `(n_findex x n_turbines)`, and I hadn't noticed that the example you mentioned doesn't follow that. It seems that what happens here is that [these lines](https://github.com/NREL/floris/blob/c50389cfba596e29dcbdfa7670db591b2b65da3e/floris/core/farm.py#L260-L264) expand the `(1 x n_turbines)` along the findices. But things would get ambiguous (and I suppose would break) if `awc_modes` was, say, size `(m x n_turbines)` with `m` not 1 or `n_findex`.

Maybe that expansion is fine and expected, but based on your comment, it looks like that doesn't always work, and perhaps we should simply enforce that `awc_modes` (and all other inputs) be passed with dimensions `(n_findex x n_turbines)`@jfrederik-nrel: Ah okay, I somehow thought that this was more common practice. If not, then I think the bug is actually in example 4, and not in @gyalla's suggested fix, and we should probably fix that example (which should be very easy). I suggest we do that as/while we merge in Gopal's pull request.@misi9170: Agreed, that sounds good!@gyalla: Great, thanks for the quick response @jfrederik-nrel and @misi9170. I've push a fixed for example number 4 on my branch as well and will open a PR to address both of these issues now. Of course, feel free to make any further modifications to the code as you see fit. I like the idea of enforcing the dimensions of input variables!@misi9170: Closed by #1081 

Solution patch:
diff --git a/examples/examples_control_types/004_helix_active_wake_mixing.py b/examples/examples_control_types/004_helix_active_wake_mixing.py
index 7738c079c..0fa1a4e9e 100644
--- a/examples/examples_control_types/004_helix_active_wake_mixing.py
+++ b/examples/examples_control_types/004_helix_active_wake_mixing.py
@@ -97,6 +97,7 @@
     np.linspace(0, 5, N),
     np.zeros(N), np.zeros(N), np.zeros(N)
     ]).reshape(4, N).T
+awc_modes = np.tile(awc_modes, (N,1)) # Repeat over N findices
 
 # Reset FlorisModel for different helix amplitudes
 fmodel.set(
diff --git a/floris/core/turbine/operation_models.py b/floris/core/turbine/operation_models.py
index 066e26b04..544086747 100644
--- a/floris/core/turbine/operation_models.py
+++ b/floris/core/turbine/operation_models.py
@@ -519,6 +519,9 @@ class AWCTurbine(BaseOperationModel):
     the other operation models.
     """
 
+    def AWC_model(a, b, c, base_values, awc_amplitudes):
+            return base_values * (1 - (b + c*base_values)*awc_amplitudes**a)
+
     def power(
         power_thrust_table: dict,
         velocities: NDArrayFloat,
@@ -537,29 +540,34 @@ def power(
             cubature_weights=cubature_weights
         )
 
-        if (awc_modes == 'helix').any():
-            if np.any(np.isclose(
-                base_powers/1000,
-                np.max(power_thrust_table['power'])
-                )):
-                raise UserWarning(
-                    'The selected wind speed is above or near rated wind speed. '
-                    '`AWCTurbine` operation model is not designed '
-                    'or verified for above-rated conditions.'
-                    )
-            return base_powers * (1 - (
-                power_thrust_table['helix_power_b']
-                + power_thrust_table['helix_power_c']*base_powers
-                )
-                *awc_amplitudes**power_thrust_table['helix_a']
-            ) # TODO: Should probably add max function here
-        if (awc_modes == 'baseline').any():
-            return base_powers
-        else:
+        valid_entries = ['helix', 'baseline']
+        if not np.all(np.isin(awc_modes, valid_entries)):
             raise UserWarning(
                 'Active wake mixing strategies other than the `helix` strategy '
                 'have not yet been implemented in FLORIS. Returning baseline power.'
-                )
+            )
+
+        # Create a copy of the base power to modify according to different AWC modes
+        powers = base_powers.copy()
+
+        helix_mask = (awc_modes == 'helix')
+        if np.any(np.isclose(base_powers[helix_mask]/1000,np.max(power_thrust_table['power']))):
+            raise UserWarning(
+                'The selected wind speed is above or near rated wind speed. '
+                '`AWCTurbine` operation model is not designed '
+                'or verified for above-rated conditions.'
+            )
+
+        awc_powers = AWCTurbine.AWC_model(
+            power_thrust_table['helix_a'],
+            power_thrust_table['helix_power_b'],
+            power_thrust_table['helix_power_c'],
+            base_powers[helix_mask],
+            awc_amplitudes[helix_mask]
+        )
+        powers[helix_mask] = awc_powers
+
+        return powers
 
 
     def thrust_coefficient(
@@ -577,20 +585,22 @@ def thrust_coefficient(
             average_method=average_method,
             cubature_weights=cubature_weights
         )
-        if (awc_modes == 'helix').any():
-            return base_thrust_coefficients * (1 - (
-                power_thrust_table['helix_thrust_b']
-                + power_thrust_table['helix_thrust_c']*base_thrust_coefficients
-                )
-                *awc_amplitudes**power_thrust_table['helix_a']
-            )
-        if (awc_modes == 'baseline').any():
-            return base_thrust_coefficients
-        else:
-            raise UserWarning(
-                'Active wake mixing strategies other than the `helix` strategy '
-                'have not yet been implemented in FLORIS. Returning baseline power.'
-                )
+
+        # Create a copy of the base thrust coefficients to modify according to different AWC modes
+        thrust_coefficients = base_thrust_coefficients.copy()
+
+        helix_mask = (awc_modes == 'helix')
+
+        awc_thrust_coefficients = AWCTurbine.AWC_model(
+            power_thrust_table['helix_a'],
+            power_thrust_table['helix_thrust_b'],
+            power_thrust_table['helix_thrust_c'],
+            base_thrust_coefficients[helix_mask],
+            awc_amplitudes[helix_mask]
+        )
+        thrust_coefficients[helix_mask] = awc_thrust_coefficients
+
+        return thrust_coefficients
 
     def axial_induction(
         power_thrust_table: dict,


Test patch:
diff --git a/tests/turbine_operation_models_unit_test.py b/tests/turbine_operation_models_unit_test.py
index b50aab54b..591c3e01c 100644
--- a/tests/turbine_operation_models_unit_test.py
+++ b/tests/turbine_operation_models_unit_test.py
@@ -568,6 +568,36 @@ def test_AWCTurbine():
     assert test_ai < base_ai
     assert test_ai > 0
 
+    # Test a mixture of "baseline" and "helix" modes
+    n_turbines = 3
+    awc_modes_test = np.array(
+        [
+            ["baseline", "baseline", "baseline"],
+            ["baseline", "helix", "helix"]
+        ]
+    )
+    test_power = AWCTurbine.power(
+        power_thrust_table=turbine_data["power_thrust_table"],
+        velocities=wind_speed * np.ones((2, n_turbines, 3, 3)), # 2 findices, 3 turbines, 3x3 grid
+        air_density=turbine_data["power_thrust_table"]["ref_air_density"],
+        awc_modes=awc_modes_test,
+        awc_amplitudes=2*np.ones((2, n_turbines)),
+    )
+    assert np.allclose(test_power[0,:], base_power)
+    assert np.allclose(test_power[1,0], base_power)
+    assert (test_power[1,1:] < base_power).all()
+
+    test_Ct = AWCTurbine.thrust_coefficient(
+        power_thrust_table=turbine_data["power_thrust_table"],
+        velocities=wind_speed * np.ones((2, n_turbines, 3, 3)), # 2 findices, 3 turbines, 3x3 grid
+        air_density=turbine_data["power_thrust_table"]["ref_air_density"],
+        awc_modes=awc_modes_test,
+        awc_amplitudes=2*np.ones((2, n_turbines)),
+    )
+    assert np.allclose(test_Ct[0,:], base_Ct)
+    assert np.allclose(test_Ct[1,0], base_Ct)
+    assert (test_Ct[1,1:] < base_Ct).all()
+
 def test_PeakShavingTurbine():
 
     n_turbines = 1