Turbine hub height optimization
Another type of design variable that can be specified is the turbines’ hub heights. In this case, the hub heights are converted into continuous variables so a gradient-based optimization with autograd can be performed. This approach is taken to speed up the optimization process.
In this example, it is necessary to specify a “GenericWindTurbine” since the hub heights need to change throughout the optimization and thus the power and thrust coefficient curve will change for each hub height studied.
Install TOPFARM if needed
[1]:
# Install TopFarm if needed
import importlib
if not importlib.util.find_spec("topfarm"):
%pip install git+https://gitlab.windenergy.dtu.dk/TOPFARM/TopFarm2.git
We now import the site, turbine positions and wake models necessary to run the optimization. The Lillgrund site will be used and its layout extracted from PyWake’s data.
[2]:
import numpy as np
import matplotlib.pyplot as plt
[3]:
from py_wake.examples.data.lillgrund import LillgrundSite, power_curve, ct_curve
from py_wake.wind_turbines._wind_turbines import WindTurbines, WindTurbine
from py_wake.wind_turbines.generic_wind_turbines import GenericWindTurbine
from py_wake.wind_turbines.power_ct_functions import PowerCtTabular
from topfarm.easy_drivers import EasyScipyOptimizeDriver
from py_wake.deficit_models.gaussian import BastankhahGaussian
from topfarm.cost_models.cost_model_wrappers import AEPCostModelComponent
from py_wake.site.shear import PowerShear
from py_wake.utils.gradients import autograd
from topfarm import TopFarmProblem
/builds/TOPFARM/TopFarm2/.pixi/envs/default/lib/python3.13/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
from .autonotebook import tqdm as notebook_tqdm
Now we specify the initial conditions and hub height boundaries. The turbines’ hub heights are staggered and assigned to every couple wind turbine in the farm.
[4]:
### Initial inputs
n_wt = 7
n_wd = 1
wd = [270]
lb = 82 # lower boundary constraint
ub = 160 # upper boundary constraint
hh = 80 # starting hub height condition
hg = 120 # second starting hub height condition
h = np.ones([n_wt]) * hh # hub height array
h_max = np.ones([n_wt]) * ub # baseline hub height
for i in range(n_wt):
if i % 2 == 0:
h[i] = hh
else:
h[i] = hg
print('wind farm hub heights:',h)
power = 2300
diameter = np.ones([n_wt]) * 93 # diameter [m]
# Site specification
h_ref = 110
alpha = 0.04
site = LillgrundSite(shear=PowerShear(h_ref, alpha))
x = np.linspace(0, 93 * 4 * n_wt, n_wt)
y = [0] * n_wt
wind farm hub heights: [ 80. 120. 80. 120. 80. 120. 80.]
Then we need to set up the GenericWindTurbine object and the WindFarmModel.
[5]:
nom_power_array = power * np.ones([n_wt]) # rated power array
class SWT23(WindTurbine): # Siemens 2.3 MW
def __init__(self, method='linear'):
"""
Parameters
----------
method : {'linear', 'pchip'}
linear(fast) or pchip(smooth and gradient friendly) interpolation
"""
WindTurbine.__init__(self, name='SWT23', diameter=93, hub_height=80,
powerCtFunction=PowerCtTabular(power_curve[:, 0], power_curve[:, 1], 'kw',
ct_curve[:, 1], method=method))
wind_turbines = WindTurbines(
names=['SWT23' for i in range(len(x))],
diameters = diameter,
hub_heights = h,
powerCtFunctions=[GenericWindTurbine(name='SWT23',
diameter = diameter[i],
hub_height = h[i],
power_norm = nom_power_array[i]).powerCtFunction for i in range(len(x))])
wf_model = BastankhahGaussian(site, wind_turbines)
/builds/TOPFARM/TopFarm2/.pixi/envs/default/lib/python3.13/site-packages/py_wake/deficit_models/gaussian.py:124: UserWarning: The BastankhahGaussian model is not representative of the setup used in the literature. For this, use py_wake.literature.gaussian_models.Bastankhah_PorteAgel_2014 instead
DeprecatedModel.__init__(self, 'py_wake.literature.gaussian_models.Bastankhah_PorteAgel_2014')
Lastly, we set up the CostModelComponent that is responsible for calculating the AEP and works as the objective function in the optimization.
[6]:
# AEP Calculation
class PyWakeAEPCostModelComponent(AEPCostModelComponent):
def __init__(self, windFarmModel, n_wt, wd=None, ws=None, max_eval=None, grad_method=autograd, n_cpu=1, **kwargs):
self.windFarmModel = windFarmModel
#objective function
def get_aep_func(h):
h_new = h[:n_wt]
simres = windFarmModel(x, y, h=h_new)
aep = simres.aep().sum()
return aep
#specifying the gradients
def daep_h(h):
return windFarmModel.aep_gradients(autograd, wrt_arg=['h'])(x, y, h)
AEPCostModelComponent.__init__(self,
input_keys=['h'],
n_wt=n_wt,
cost_function=get_aep_func,
cost_gradient_function=daep_h,
output_unit='GWh',
max_eval=max_eval, **kwargs)
cost_comp = PyWakeAEPCostModelComponent(windFarmModel=wf_model, n_wt=len(x), grad_method=autograd, n_cpu=1, wd=None, ws=None)
Lastly, we set up the TopFarmProblem along with some optimization parameters.
[7]:
# optimization specs and problem formulation
maxiter = 20
tol = 1e-6
problem = TopFarmProblem(design_vars= {'h':(h, lb, ub)},
cost_comp=cost_comp,
driver=EasyScipyOptimizeDriver(optimizer='SLSQP', maxiter=maxiter, tol=tol),
n_wt=n_wt,
expected_cost=0.001
)
[8]:
_,state,_=problem.optimize()
Optimization terminated successfully (Exit mode 0)
Current function value: -62751.3352971155
Iterations: 16
Function evaluations: 17
Gradient evaluations: 16
Optimization Complete
-----------------------------------
Now we plot the turbines in the XZ plane to visualize the final hub heights.
[9]:
from py_wake import XZGrid
h = np.around(state['h'])
print('final hub heights:',h)
#taking only the first row of turbines
x = x[:6]
y = y[:6]
h = h[:6]
sim_res_ref = wf_model(x, y, wd=[270])
sim_res_opt = wf_model(x, y, h=h, wd=[270])
plt.figure(figsize=(12,4))
sim_res_opt.flow_map(XZGrid(y=0, z=np.linspace(0.5, 1000))).plot_wake_map()
plt.ylabel('Height [m]')
plt.xlabel('x [m]')
final hub heights: [ 99. 160. 92. 160. 92. 160. 108.]
[9]:
Text(0.5, 0, 'x [m]')
[10]:
x_locs = [-100, 1000, 2000]
box = sim_res_opt.flow_box(x=x_locs, y=0, h=np.linspace(0, 1000))
[11]:
plt.plot(box.WS_eff.isel(ws=5).values[0, :, :, 0].T, box.h, label=x_locs)
leg = plt.legend(title='X Location (m)')
plt.xlabel('Speed (m/s)')
plt.ylabel('Height (m)')
[11]:
Text(0, 0.5, 'Height (m)')
