Problems

Try this yourself (requires google account)

In TOPFARM the problem is the collection/container of components and drivers. Optimizing a problem is executing a workflow - finding the best feasible solution under a series of constraints using a specific driver and initial conditions. To the problem is also attached a plotting routine so one can follow the optimization path as well as a recorder so that intermediate variable values are accessible after the optimization is finished.

The TopfarmProblem inherits its fundamental nature from the OpenMDAO problem, and is being adapted so it can connect the given workflow with the driver. For example, if the user specifies a boundary constraint and this is not supported by the driver, the Problem is equipped with a penalty component that will deter the driver giving unfeasible solutions. Or if your workflow does not have gradients for all components and a gradient based driver is specified, finite differencing is applied to obtain the gradients.

The TopFarmProblem is defined as a class, and it is comprised of the following elements:

def __init__(self, design_vars, cost_comp=None, driver=EasyScipyOptimizeDriver(),
             constraints=[], plot_comp=NoPlot(), record_id=None,
             expected_cost=1, ext_vars={}, post_constraints=[], approx_totals=False,
             recorder=None, additional_recorders=None,
             n_wt=0, grid_layout_comp=None):

It takes the design variables, cost model component, driver and constraints as the most important parameters. In addition, it has a feature to visualize the optimization with a plotting component, and the use of additional constraints (for example, fatigue loading) can be specified with the post_constraints option. The latest development also allows the user to have a regular wind farm layout instead of irregular and individual position of turbines during a layout optimization. To do this, use the grid_layout_comp.

Before working on the notebook exercises, make sure that you import or install TOPFARM on your newly created environment.

# Install TopFarm if needed
import importlib
if not importlib.util.find_spec("topfarm"):
    !pip install git+https://gitlab.windenergy.dtu.dk/TOPFARM/TopFarm2.git

Import supporting libraries in Python

import numpy as np
import matplotlib.pyplot as plt

Next we import and initialize several functions and classes from TOPFARM to set up the problem including:

• TopFarmProblem. Overall topfarm problem class to which the objectives, design variables, and constraints are added

• EasyScipyOptimizeDriver. A subclass of ScipyOptimizeDriver which is configured for the given workflow

• get_iea37_initial, get_iea37_constraints, get_iea37_cost. Functions to get the initial layout, the constraints and the cost function for the IEA task 37 benchmark example

• NoPlot, XYPlotComp. Plotting components to visualize the results

from topfarm import TopFarmProblem
from topfarm.easy_drivers import EasyScipyOptimizeDriver
from topfarm.examples.iea37 import get_iea37_initial, get_iea37_constraints, get_iea37_cost
from topfarm.plotting import NoPlot, XYPlotComp

Turbine location optimization - Example 1

This example optimizes the locations of the 9-turbine benchmark wind farm from IEA Task 37 using the provided initial locations and the EasyScipyOptimizeDriver. Details on the benchmark can be found in the following reference:

• Baker et al. (2019) “Best Practices for Wake Model and Optimization Algorithm Selection in Wind Farm Layout Optimization”. AIAA 2019.

First we set up the number of turbines and initial positions.

n_wt = 9          #valid number of turbines are: 9, 16, 36, 64
x, y = get_iea37_initial(n_wt).T

Now we import from PyWake the wake model, wind turbine and site objects to define our objective function, which is represented by the AEP (see https://topfarm.pages.windenergy.dtu.dk/PyWake/index.html for more information).

from py_wake.deficit_models.gaussian import IEA37SimpleBastankhahGaussian   #wake model
from py_wake.examples.data.iea37 import IEA37_WindTurbines, IEA37Site       #wind turbines and site used

site = IEA37Site(n_wt)
wind_turbines = IEA37_WindTurbines()

wake_model = IEA37SimpleBastankhahGaussian(site, wind_turbines)
wd = np.linspace(0., 360., 16, endpoint=False)                              #number of wind directions to study

#objective function
def aep_func(x,y,wd=wd):
    sim_res = wake_model(x,y, wd=wd)
    aep = sim_res.aep().sum()
    return aep

/builds/TOPFARM/TopFarm2/.pixi/envs/default/lib/python3.11/site-packages/py_wake/deficit_models/gaussian.py:277: UserWarning: The IEA37SimpleBastankhahGaussian model is not representative of the setup used in the literature. For this, use py_wake.literature.iea37_case_study1.IEA37CaseStudy1 instead
  DeprecatedModel.__init__(self, 'py_wake.literature.iea37_case_study1.IEA37CaseStudy1')

Then we set up the CostModelComponent that will be used for the optimization as well as the constraints of the problem.

The cost model component is defined as:

def __init__(self, input_keys, n_wt, cost_function, cost_gradient_function=None,
             output_keys=["Cost"], output_unit="", additional_input=[], additional_output=[], max_eval=None,
             objective=True, maximize=False, output_vals=[0.0], input_units=[], step={}, use_penalty=True, **kwargs):

where it takes the design variables as input keys and the objective function as the cost function. In addition, the gradient function can be specified when necessary. By default, the objective is set to True and the maximize to False, which treats the optimization as a minimization problem.

from topfarm.cost_models.cost_model_wrappers import CostModelComponent

constraint = get_iea37_constraints(n_wt)

aep_comp = CostModelComponent(input_keys=['x','y'],
                              n_wt=n_wt,
                              cost_function=aep_func,
                              output_keys=[('AEP', 0)],
                              output_unit="GWh",
                              objective=True,
                              maximize=True
                             )

Lastly, we set up the TopFarmProblem with the previously defined CostModelComponents and design variables x and y.

problem = TopFarmProblem(design_vars={'x': x, 'y': y},
                          n_wt=n_wt,
                          cost_comp=aep_comp,
                          constraints=constraint,
                          driver=EasyScipyOptimizeDriver(disp=False),
                          plot_comp=XYPlotComp()
                        )

INFO: checking out_of_order...
INFO:     out_of_order check complete (0.000496 sec).
INFO: checking system...
INFO:     system check complete (0.000020 sec).
INFO: checking solvers...
INFO:     solvers check complete (0.000196 sec).
INFO: checking dup_inputs...
INFO:     dup_inputs check complete (0.000035 sec).
INFO: checking missing_recorders...
INFO:     missing_recorders check complete (0.000005 sec).
INFO: checking unserializable_options...
INFO:     unserializable_options check complete (0.000175 sec).
INFO: checking comp_has_no_outputs...
INFO:     comp_has_no_outputs check complete (0.000044 sec).
INFO: checking auto_ivc_warnings...
INFO:     auto_ivc_warnings check complete (0.000004 sec).

_, state,_ = problem.optimize()

INFO: checking out_of_order...
INFO:     out_of_order check complete (0.000173 sec).
INFO: checking system...
INFO:     system check complete (0.000018 sec).
INFO: checking solvers...
INFO:     solvers check complete (0.000141 sec).
INFO: checking dup_inputs...
INFO:     dup_inputs check complete (0.000030 sec).
INFO: checking missing_recorders...
INFO:     missing_recorders check complete (0.000006 sec).
INFO: checking unserializable_options...
INFO:     unserializable_options check complete (0.000138 sec).
INFO: checking comp_has_no_outputs...
INFO:     comp_has_no_outputs check complete (0.000040 sec).
INFO: checking auto_ivc_warnings...
INFO:     auto_ivc_warnings check complete (0.000003 sec).

Ignoring fixed y limits to fulfill fixed data aspect with adjustable data limits.

Ignoring fixed y limits to fulfill fixed data aspect with adjustable data limits.

state

{'x': array([-4.94541050e+02, -6.29539540e+02, -7.48570820e+02,  1.09396139e+02,
        -1.72143211e-01, -1.09518679e+02,  7.49035949e+02,  6.30012909e+02,
         4.94529807e+02]),
 'y': array([-751.95023118, -146.37045879,  489.61739038, -629.91381474,
          -1.2565284 ,  630.35702783, -489.33957981,  143.80596575,
         751.95762515])}

Turbine location optimization - Example 2

Alternatively, we can use the PyWakeAEPCostModelComponent as the cost model for the optimization.

from py_wake.deficit_models.gaussian import IEA37SimpleBastankhahGaussian   #wake model
from py_wake.examples.data.iea37 import IEA37_WindTurbines, IEA37Site       #wind turbines and site used
from topfarm.cost_models.py_wake_wrapper import PyWakeAEPCostModelComponent #cost model

#creating a function to simulate the flow conditions of the IEA37 Site, specifying the desired number of turbines and wind directions
def get_iea37_cost(n_wt=9, n_wd=16):

    """Cost component that wraps the IEA 37 AEP calculator"""
    wd = np.linspace(0., 360., n_wd, endpoint=False)
    site = IEA37Site(n_wt)
    wind_turbines = IEA37_WindTurbines()
    wake_model = IEA37SimpleBastankhahGaussian(site, wind_turbines)

    return PyWakeAEPCostModelComponent(wake_model, n_wt, wd=wd)

Now we create a function that specifies the optimization problem, with the desired design variables, constraints, and driver.

def optimize_iea37_locs(n_wt, n_wd, driver, state=None):
    """
    Parameters
    ----------
    - n_wt: int
        Number of wind turbines

    - n_wd: int
        Number of wind directions to consider for the AEP

    - driver: TopfarmDriver instance
        The optimization algorithm to use

    - state: dict(x=[], y=[]) [default=None]
        State to start from the optimization

    Returns
    -------
    - state: The state after the optimization

    """

    initial = get_iea37_initial(n_wt)
    design_vars = dict(zip('xy', (initial[:, :2]).T))

    tf = TopFarmProblem(
            design_vars,
            get_iea37_cost(n_wt, n_wd=n_wd),
            constraints=get_iea37_constraints(n_wt),
            driver=driver,
            plot_comp=XYPlotComp())

    if not state:
        _, state = tf.evaluate()

    _, state, _ = tf.optimize(state)
    return state,

state = optimize_iea37_locs(9, 16, EasyScipyOptimizeDriver(disp=False))

INFO: checking out_of_order...
INFO:     out_of_order check complete (0.000191 sec).
INFO: checking system...
INFO:     system check complete (0.000018 sec).
INFO: checking solvers...
INFO:     solvers check complete (0.000164 sec).
INFO: checking dup_inputs...
INFO:     dup_inputs check complete (0.000027 sec).
INFO: checking missing_recorders...
INFO:     missing_recorders check complete (0.000006 sec).
INFO: checking unserializable_options...
INFO:     unserializable_options check complete (0.000137 sec).
INFO: checking comp_has_no_outputs...
INFO:     comp_has_no_outputs check complete (0.000036 sec).
INFO: checking auto_ivc_warnings...
INFO:     auto_ivc_warnings check complete (0.000004 sec).
INFO: checking out_of_order...
INFO:     out_of_order check complete (0.000155 sec).
INFO: checking system...
INFO:     system check complete (0.000016 sec).
INFO: checking solvers...
INFO:     solvers check complete (0.000140 sec).
INFO: checking dup_inputs...
INFO:     dup_inputs check complete (0.000027 sec).
INFO: checking missing_recorders...
INFO:     missing_recorders check complete (0.000006 sec).
INFO: checking unserializable_options...
INFO:     unserializable_options check complete (0.000142 sec).
INFO: checking comp_has_no_outputs...
INFO:     comp_has_no_outputs check complete (0.000036 sec).
INFO: checking auto_ivc_warnings...
INFO:     auto_ivc_warnings check complete (0.000004 sec).

/builds/TOPFARM/TopFarm2/.pixi/envs/default/lib/python3.11/site-packages/py_wake/deficit_models/gaussian.py:277: UserWarning: The IEA37SimpleBastankhahGaussian model is not representative of the setup used in the literature. For this, use py_wake.literature.iea37_case_study1.IEA37CaseStudy1 instead
  DeprecatedModel.__init__(self, 'py_wake.literature.iea37_case_study1.IEA37CaseStudy1')
Ignoring fixed y limits to fulfill fixed data aspect with adjustable data limits.

INFO: checking out_of_order...
INFO:     out_of_order check complete (0.000205 sec).
INFO: checking system...
INFO:     system check complete (0.000023 sec).
INFO: checking solvers...
INFO:     solvers check complete (0.000138 sec).
INFO: checking dup_inputs...
INFO:     dup_inputs check complete (0.000034 sec).
INFO: checking missing_recorders...
INFO:     missing_recorders check complete (0.000008 sec).
INFO: checking unserializable_options...
INFO:     unserializable_options check complete (0.000154 sec).
INFO: checking comp_has_no_outputs...
INFO:     comp_has_no_outputs check complete (0.000040 sec).
INFO: checking auto_ivc_warnings...
INFO:     auto_ivc_warnings check complete (0.000003 sec).

Ignoring fixed y limits to fulfill fixed data aspect with adjustable data limits.

We can see that the optimized AEP corresponds to 239.69 GWh and an improvement of 34.37% is found. However, these value is lower than the one reported in Baker et al. (2019) of 257.790 GWh.

Moreover, the layout does not match the figures given in Appendix A in the same reference. This is due to the fact that the SLSQP optimizer was attracted to a local minimum. To find the global optimum, more advanced optimization procedures should be used. This benchmark is discussed in more detail in the validation report linked in TOPFARM’s documentation.

To obtain the final turbine positions that correspond to the optimum AEP found, we need to call the state object.

state

({'x': array([-4.94541038e+02, -6.29538241e+02, -7.48569543e+02,  1.09398064e+02,
         -1.70674714e-01, -1.09519176e+02,  7.49037349e+02,  6.30015331e+02,
          4.94530553e+02]),
  'y': array([-751.9502388 , -146.37144376,  489.61742462, -629.91364594,
           -1.25798282,  630.35702136, -489.33930019,  143.80366339,
          751.9571348 ])},)