SciPy Yaw optimization classes: improved efficiency, automated downstream turbine exclusion and turbine-dependent weighing terms

[Bartdoekemeijer]

This pull request is ready to be merged.

Feature or improvement description
Improvements to the SciPy yaw optimization classes for improved efficiency, automated downstream turbine exclusion and turbine-dependent weighing terms in the objective definition.

Specifically, the novel additions are:

• The variable bnds is now expanded to remove control variables in the optimization problem whenever a turbine's specified lower bound equals its upper bound, within a certain margin (0.001 deg). Removing these turbines from the optimization speeds up the optimization, since the gradients of the objective function w.r.t. these variables need not be calculated and the dimensionality of the problem is smaller. This also circumvents certain numerical issues that I was experiencing in the SciPy optimization with recent versions of SciPy and numpy. Practically, the usage of bnds stays exactly the same as before -- it is just handled in a smarter way behind the scenes.
• Additionally, bnds now allows one to restrict a certain turbine to a nonzero misalignment angle, yet exclude it from the optimization problem. For example, one can opt to fix an upstream turbine to 20 degrees and only optimize the yaw angle for the second turbine in the row, simply by bnds[0] = [20.0, 20.0]. Previously, in several parts in the code, the yaw angle of such a turbine would fall back to 0.0 instead of its assigned value.
• I have added an input to the various SciPy yaw optimization classes called turbine_weights. During each objective function call, the turbine power values from floris are multiplied with their corresponding weights as defined in this array. By setting certain turbines' weights to 0, the user can exclude turbines from the objective/cost function. For example, one may optimize a three-turbine array while modeling a larger wind farm. With this new input variable, one can simply assign zero weights to all turbines except the three-turbine array of interest. Yet, the other turbines are still part of each function call's simulation, thereby including their wake effects. turbine_weights defaults to an array with length equal to the number of turbines, and all values being 1.0.
• I have added an input variable to the SciPy yaw optimization classes called exclude_downstream_turbines. This is a boolean which, when defined as True, will exclude all turbines that have a wake which does not impinge any other turbine from the optimization problem. Hence, this makes it easy for the user to immediately remove any turbines that do not have any effect on the cost function, rather than having the optimization algorithm figure this out on its own. Here are some figures of the simplified wake model used to derive whether a turbine is waked or not, where the green dots indicate most-downstream turbines:
  
  
  
  
  Additionally, this is often more accurate than letting the optimization figure out certain turbines do not affect the cost function. Namely, those turbines' yaw angles do not necessarily converge to zero but instead often converge to a nonzero small value.

Related issue, if one exists

This is a more commonly requested feature within the wind farm controls group at NREL. This also fixes issue #207 .

Impacted areas of the software

The SciPy optimization module.

Additional supporting information

@ejsimley has volunteered as a reviewer for this PR.

Test results, if applicable

Here is a simple test script based on the optimize_yaw_wind_rose.py example:

import os

import floris.tools as wfct
from floris.tools.optimization.scipy.yaw_wind_rose import YawOptimizationWindRose

# Instantiate a FLORIS object with 4 turbines
file_dir = os.path.dirname(os.path.abspath(__file__))
fi = wfct.floris_interface.FlorisInterface(
    os.path.join(file_dir, "../../../example_input.json")
)
fi.reinitialize_flow_field(
    layout_array=[[0., 600., 1200., 1800.], [0., 0., 0, 0.]]
)

# Instantiate the Optimization object
yaw_opt = YawOptimizationWindRose(
    fi=fi,
    wd=[90., 180., 270., 45.],
    ws=[8., 8., 8., 8.],
    minimum_yaw_angle=-25.,
    maximum_yaw_angle=25.,
    bnds=[[0., 25.], [2., 15.], [5., 7.], [3., 3.]],
    opt_options={"maxiter": 100,"disp": True,"iprint": 2,"ftol": 1e-7},
    turbine_weights=[0., 1., 1., 1.],
    exclude_downstream_turbines=True,
    verbose=False
)

# Perform optimization
df_opt = yaw_opt.optimize()

print(df_opt)

which outputs:

    ws     wd     power_opt                                  turbine_power_opt                                         yaw_angles
0  8.0   90.0  3.046841e+06  [546299.7970654803, 469420.7167990642, 339859....                 [0.0, 2.0, 7.000000000000001, 3.0]
1  8.0  180.0  6.764128e+06  [1695368.645547269, 1693542.8021105803, 168395...                               [0.0, 2.0, 5.0, 3.0]
2  8.0  270.0  3.783548e+06  [1409447.771337562, 734251.9139146324, 830846....  [25.0, 14.999999999999995, 6.999999999999999, ...
3  8.0   45.0  6.764128e+06  [1695368.645547269, 1693542.8021105803, 168395...                               [0.0, 2.0, 5.0, 3.0]

Which seems like the right solutions, namely:

• The optimal yaw angles never exceed their bounds.
• The power production of the first turbine is weighed as 0.0 in turbine_weights. For the wd=90.0 situation, turbine 0 is most downstream, and turbine 1 is directly impacting exclusively turbine 0. Hence, the yaw angle of turbine 1 is also as small as possible (2.0, its lower bound). The yaw angle of turbine 2 is as high as possible (7.0, upper bound) to steer the wake as far as possible away from turbine 1. The yaw angle of the last turbine is fixed at 3.0 through the prespecified bounds. This turbine is therefore excluded from the optimization, and 3.0 is indeed its found optimal value.
• For wd=180. and wd=45., all turbines operate in freestream. Therefore, all turbines have a minimal yaw angle, according to their bounds.
• For wd=270. degrees, we experience the traditional aligned optimization problem. Turbines 0, 1 and 2 all have their yaw angle equal to their upper bound. The power of turbine 0 is ignored in the total cost, but this does not significantly impact the optimization: it would only push for more wake steering in turbine 0.

I have also tested the parallelized/MPI optimization code on Eagle successfully.

[ejsimley]

Hey Bart, overall I think this is really nice, thanks for working on this. I like how you included the option to weight individual turbine powers instead of just having a binary option to include or not include turbines in the objective function. Also great that you modified to speed up the optimizations when bounds are set to the same value. I left several comments where I think some changes are needed, and also a bigger comment about whether we should use some existing functionality to determine downstream turbines.

[Bartdoekemeijer]

Hi @ejsimley. Thank you very much for the detailed comments throughout the PR! I processed all your comments and hopefully resolved all of them. Please take a second look whenever is convenient for you. Thanks!

[ejsimley]

I don't think this needs to be a warning. For example, I imagine it would be very common for the baseline values to be zero, but the initial guess for optimal offsets to be some other value.

[ejsimley]

Same thing here, where setting the yaw angles for the uncontrolled turbines to their baseline values makes more sense to me. Especially since only a single set of x0 values can be specified for all wind directions. Depending on the wind direction, different turbines will be downstream (assuming downstream turbines are excluded), so it seems safer to set those yaw angles to their baseline values.

[ejsimley]

flagging this as another place where, to me, baseline yaw angles makes more sense than x0.

[ejsimley]

Hey @Bartdoekemeijer, the updates are really nice and resolve most of the issues I could think of. I left a few minor comments, though, mostly about whether x0 or x_baseline should be used for the uncontrolled turbines.

[Bartdoekemeijer]

Hey @ejsimley. Thank you for these suggestions! I think we are really converging now and the code is definitely a lot clearer and organized after addressing your comments. I have made the desired changes. The biggest change is that I introduced an additional input to the yaw optimization classes called yaw_angles_baseline, to allow the user to define the baseline yaw angles used to calculate the initial and baseline power productions. If this variable is unspecified, then it falls back to the yaw angles inherent to the floris object. If any of those yaw angles are nonzero, an INFO statement will be printed to inform the user of this.

Please review them whenever is convenient for you. Thank you again!

[ejsimley]

OK, I think adding the yaw_angles_baseline variable was a nice way to solve some of the ambiguity! And thanks for the expanded documentation on the simple wake model. There's just one scenario I want to make sure I understand correctly. If we want to exclude some turbines from yaw optimization by setting their lower and upper bnds to the same value, the actual yaw offset used will be the corresponding value in yaw_angles_baseline, even if it differs from the bnds value, right? So if we set bnds[i][0] = bnds[i][1] = 10 but yaw_angles_baseline[i] = 0, it will use the value of 0 instead of 10. Do you think it makes more sense to set the "optimal_yaw_angles" to the values in bnds for turbines where the bounds are equal? Otherwise, I would say just make it clear in the docstrings which value will get used in cases where the the lower and upper bounds are the same but the value in yaw_angles_baseline is different.

[Bartdoekemeijer]

Yes, indeed, the values for bnds are ignored if they are equal and instead self.yaw_angles_baseline is used for those turbines. Again rethinking this because you're rightfully pointing out that this is confusing. To clarify it for myself: summarizing what should be happening:

• If turbines are completely downstream and exclude_downstream_turbines=True, then we should set those turbines' yaw angles to 0.0 or the value closest to 0.0, bounds-permitting.
• If turbines have their upper and lower yaw bound equal to one another, then that's their optimal yaw angle by definition.
• The remaining yaw angles are part of self.turbs_to_opt and should be a result of the optimization.

This means self.yaw_angles_baseline does not appear at all in the optimization of the yaw angles. I actually do like this a lot better. Everything is then defined by the boundary conditions. Also, we should in this case not initialize self.x0=self.yaw_angles_baseline if x0=None is provided by the user. Namely, this initial condition can easily fall out of the specified bounds. Instead, x0 should represent a fair guess of the optimal solution and therefore I have specified it as:
x0=0.0 wherever it falls within the turbine yaw bounds, and x0=np.mean(self.bnds[ti]) wherever 0.0 is not included in the allowable yaw range for turbine ti.

[ejsimley]

That makes a lot of sense to me Bart! Regarding your first bullet point, I just remember you brought up the point that a user might want to force some turbines to have non-zero yaw angles to represent a yaw error bias or some other fault scenario. But maybe that was more for the baseline yaw offsets, not the optimized offsets. Either way, it is fine with me to set the downstream turbines to 0 yaw.

[Bartdoekemeijer]

Yes! I do want to facilitate the option to fix downstream turbines (or any turbines for that matter) to a non-zero yaw angle. By setting those turbines' yaw angles to a value closest to 0.0 bounds-permitting, then it will also meet any equality constraints specified by the user. Namely, if the lower bound and upper bound are the same value, then the "value closest to 0.0" is the lower bound (or the upper bound, for that matter). I also want to have this option for evaluating the baseline case but that should be separate from the actual optimization, which I will also facilitate. Finally, x0 also is completely separate from all of this. The only requirement for x0 is that it meets the bnds specified by the user.

[Bartdoekemeijer]

@ejsimley I think the latest commits should do it. Let me know what you think. Thanks!

[ejsimley]

Hey Bart, I think your latest changes address all the cases we could think of! It's possible there might be some use cases that aren't handled by the code now, but those can be addressed later if they come up. I left a couple comments where documentation should be cleaned up (one is a suggestion). After you address those, fine with me to merge!

[Bartdoekemeijer]

Hi @ejsimley. I just updated the docstrings. I feel like I keep forgetting small things -- thanks for reviewing everything with such eye for detail! Hopefully I didn't miss anything else now and we can get this merged into the develop branch.

[ejsimley]

Thanks @Bartdoekemeijer! The documentation is very clear now. Appreciate all your work on this! I know I'll be using these new features a lot.