Open and run in Colab (interactive) Edit on Gitlab

Optimization with TOPFARM

This section describes two optimization examples: a layout optimization with AEP and power production optimization with a de-ratable wind turbine.

Install PyWake if needed

# Install PyWake if needed
    import py_wake
except ModuleNotFoundError:
    !pip install git+

In the following sections the examples are optimized using TOPFARM, an open source Python package developed by DTU Wind Energy to help with wind-farm optimizations. It has a lot of nice built-in features and wrappers.

Install TOPFARM if needed

# Install TopFarm if needed
    import topfarm
except ImportError:
    !pip install topfarm --user

Example 1 - Optimize AEP wrt. wind turbine position (x,y)

In this example we optimize the AEP of the IEAWind Task 37 Case Study 1.

As TOPFARM already contains a cost model component, PyWakeAEPCostModelComponent, for this kind of problem, setting up the problem is really simple.

import numpy as np
import matplotlib.pyplot as plt
import multiprocessing

#setting up pywake models
from py_wake.deficit_models.gaussian import IEA37SimpleBastankhahGaussian, BastankhahGaussian
from import IEA37Site, IEA37WindTurbines
from import Hornsrev1Site, V80
from py_wake.utils.gradients import autograd, fd, cs
from py_wake.utils.plotting import setup_plot

#setting up topfarm models
from topfarm._topfarm import TopFarmProblem
from topfarm.cost_models.py_wake_wrapper import PyWakeAEPCostModelComponent
from topfarm.constraint_components.boundary import CircleBoundaryConstraint, XYBoundaryConstraint
from topfarm.constraint_components.spacing import SpacingConstraint
from topfarm.easy_drivers import EasyScipyOptimizeDriver

#wind farm model for the IEA 37 site
def IEA37_wfm(n_wt, n_wd):
    site = IEA37Site(n_wt)
    site.default_wd = np.linspace(0,360,n_wd, endpoint=False)
    wt = IEA37WindTurbines()
    return IEA37SimpleBastankhahGaussian(site, wt)

#wind farm model for the Hornsrev1 site
Hornsrev1_wfm = BastankhahGaussian(Hornsrev1Site(), V80())

#function to create a topfarm problem, following the elements of OpenMDAO architecture
def get_topfarmProblem_xy(wfm, grad_method, maxiter, n_cpu):
    x, y =
    boundary_constr = [XYBoundaryConstraint(np.array([x, y]).T),
                       CircleBoundaryConstraint(center=[0, 0], radius=np.round(np.hypot(x, y).max()))][int(isinstance(, IEA37Site))]

    return TopFarmProblem(design_vars={'x': x, 'y': y},
                          cost_comp=PyWakeAEPCostModelComponent(windFarmModel=wfm, n_wt=len(x),
                                                                grad_method=grad_method, n_cpu=n_cpu,
                                       SpacingConstraint(min_spacing=2 * wfm.windTurbines.diameter())])
#we create a function to optimize the problem and plot the results in terms of AEP and simulation time
def optimize_and_plot(wfm, maxiter, skip_fd=False):
    for method, n_cpu in [(fd,1), (autograd,1), (autograd, multiprocessing.cpu_count())][int(skip_fd):]:
        tf = get_topfarmProblem_xy(wfm,method,maxiter,n_cpu)
        cost, state, recorder = tf.optimize(disp=True)
        t,aep = [recorder[v] for v in ['timestamp','AEP']]
        plt.plot(t-t[0],aep, label=f'{method.__name__}, {n_cpu}CPU(s)')
    n_wt, n_wd, n_ws = len(, len(, len(

    setup_plot(ylabel='AEP [GWh]', xlabel='Time [s]',title = f'{n_wt} wind turbines, {n_wd} wind directions, {n_ws} wind speeds')
optimize_and_plot(IEA37_wfm(16, n_wd=16), maxiter=3)

Pre-computed result of the AEP during 100 iterations of optimization of the IEA task 37 case study 1 (16 wind turbines, 12 wind directions and one wind speed) plotted as a function of time.

Autograd is seen to be faster than finite difference and for this relatively small problem, 1 CPU is faster than 32CPUs.


optimize_and_plot(Hornsrev1_wfm, 100, skip_fd=True)

Precomputed result of the AEP during 100 iterations of optimization of the Hornsrev 1 wind farm (80 wind turbines, 360 wind directions and 23 wind speed) plotted as a function of time.

In this case the optimization with 32 CPUs is around 6 times faster than the optimization with 1 CPU. image4.png

Example 2 - Optimize WS, TI, Power and custom functions

To optimize some output, y, with respect to some input, x, you simply need to setup a function, y = f(x).

In the examle below, we will use a wind turbine that can be de-rated.

De-ratable wind turbine

The relation between power and CT of the de-ratable wind turbine is obtained from 1D momentum theory.

import autograd.numpy as np
import matplotlib.pyplot as plt

from py_wake.wind_turbines._wind_turbines import WindTurbine
from py_wake.wind_turbines.power_ct_functions import PowerCtFunction
from py_wake.utils.model_utils import fix_shape

def power_ct(ws, derating, run_only):
    derating = fix_shape(derating, ws)
    cp = 16 / 27 * (1 - derating)
    power = np.maximum(0, 1 / 2 * 1.225 * 50**2 * np.pi * cp * ws ** 3)

    # solve cp = 4 * a * (1 - a)**2 for a
    y = 27.0 / 16.0 * cp
    a = 2.0 / 3.0 * (1 - np.cos(np.arctan2(2 * np.sqrt(y * (1.0 - y)), 1 - 2 * y) / 3.0))
    ct = 4 * a * (1 - a)
    return [power, ct][run_only]

powerCtFunction = PowerCtFunction(input_keys=['ws', 'derating'], power_ct_func=power_ct, power_unit='w')
wt = WindTurbine(name="MyWT", diameter=100, hub_height=100, powerCtFunction=powerCtFunction)

The power and CT curves as a function wind speed is plotted below for 0, 5 and 10% derating.

ax1 = plt.gca()
ax2 = plt.twinx(ax1)
ws = np.linspace(3, 20)
for derating in [0, .05, .1]:
    ct = wt.ct(ws, derating=derating)
    ax1.plot(ws, wt.power(ws, derating=derating) / 1e6, label='%d%% derating, ct=%.2f' % (derating * 100, ct[0]))
    ax2.plot(ws, ct,'--')

ax1.legend(loc='lower right')
ax1.set_xlabel('Wind speed [m/s]')
ax1.set_ylabel('Power [MW]')
ax2.set_ylim([0, 1])
(0.0, 1.0)

Maximize mean power by optimizing de-rating factor and hub height

In this example we will maximize the mean power by optimizing the individual wind turbine hub height and derating factors.

First we setup the WindFarmModel and the function to maximize, mean_power, which takes the hub height and derating factors as input:

from py_wake.deficit_models.gaussian import IEA37SimpleBastankhahGaussian
from import UniformSite
from py_wake.utils.gradients import autograd

from topfarm._topfarm import TopFarmProblem
from topfarm.cost_models.cost_model_wrappers import CostModelComponent
from topfarm.easy_drivers import EasyScipyOptimizeDriver

n_wt = 5
wfm = IEA37SimpleBastankhahGaussian(site=UniformSite(p_wd=[1],ti=.1), windTurbines=wt)
wt_x = np.arange(n_wt) * 4 * wt.diameter()
wt_y = np.zeros_like(wt_x)

def mean_power(zhub, derating):
    power_ilk = wfm.calc_wt_interaction(x_i=wt_x, y_i=wt_y, h_i=zhub, wd=[270], ws=[10], derating=derating)[2]
    return np.mean(power_ilk)

Setup the gradient function with respect to both arguments.

Again the PyWake autograd method will, under the hood, calculate the gradients with respect to both inputs in one go.

dmean_power_dzhub_derating = autograd(mean_power,argnum=[0,1])

Initialize zhub and derating. The values are choosen to avoid zero gradients (e.g. if all wt has same hub height)

zhub = np.arange(n_wt)+100 # 100,101,102,...
print ('Mean power: %f MW'%(mean_power(zhub,derating)/1e6))
print ('Gradients:',dmean_power_dzhub_derating(zhub,derating))
Mean power: 1.429855 MW
Gradients: [array([-132.37798792,  -20.72030943,    9.44975613,   41.70595145,
        101.94258977]), array([   5539.55413631,  170225.9858809 ,  132322.76792935,
         39204.83341772, -234801.3362418 ])]

Next step is to setup the CostModelComponent and TopFarmProblem

cost_comp=CostModelComponent(input_keys=['zhub', 'derating'],
                              maximize=True # because we want to maximize the mean power

tf = TopFarmProblem(design_vars={
                        # (initial_values, lower limit, upper limit)
                        'zhub': ([100]*n_wt, 80, 130),
                        'derating': ([0] * n_wt, 0, .9)},
                    expected_cost=1000, # expected cost impacts the size of the moves performed by the optimizer
INFO: checking out_of_order
INFO: checking system
INFO: checking solvers
INFO: checking dup_inputs
INFO: checking missing_recorders
INFO: checking unserializable_options
INFO: checking comp_has_no_outputs
INFO: checking auto_ivc_warnings

As seen above the gradient of the mean power wrt. hub height is zero when all wind turbines have the same height. This means that the optimizer “thinks” that the solution cannot be improved. We therefore need to initialize the optimization with slightly different hub heights.

Furthermore, the derating must be above zero to avoid inequality constraint failure.

#perform the optimization
cost, state, recorder = tf.optimize(state={'zhub':np.arange(n_wt)+100,   # 100,101,102,...
                                           'derating':[0.1]*n_wt         # 10% initial derating
print ()
print ('Optimized mean power %f MW'% (cost/1e6))
print ('Optimized state', state)
INFO: checking out_of_order
INFO: checking system
INFO: checking solvers
INFO: checking dup_inputs
INFO: checking missing_recorders
INFO: checking unserializable_options
INFO: checking comp_has_no_outputs
INFO: checking auto_ivc_warnings
Optimization terminated successfully    (Exit mode 0)
            Current function value: -1716.1873043787978
            Iterations: 48
            Function evaluations: 92
            Gradient evaluations: 44
Optimization Complete

Optimized mean power -1.716187 MW
Optimized state {'zhub': array([ 80., 130.,  80.,  80., 130.]), 'derating': array([8.69167702e-02, 6.26083692e-02, 1.99902808e-01, 4.53089707e-02,

Plotting the results

from py_wake import XZGrid
derating = state['derating']
h = state['zhub']

sim_res_ref = wfm(wt_x, wt_y, wd=[270], ws=[10], derating=[0] * n_wt)
sim_res_opt = wfm(wt_x, wt_y, h=h, wd=[270], ws=[10], derating=derating)
for x_, d in zip(wt_x, derating):
    plt.text(x_ + 50, -80, "%d%%" % np.round(d * 100), fontsize=10)
plt.ylabel('Height [m]')
plt.xlabel('x [m]')

for i, (sim_res, l) in enumerate([(sim_res_ref, 'Baseline'), (sim_res_opt, 'Optimized')]): + i * .4, sim_res.Power.squeeze() * 1e-6, width=.4,
            label='%s, mean power= %.2fMW' % (l, sim_res.Power.mean() * 1e-6))
plt.ylabel('Power [MW]')
plt.xlabel('Wind turbine')