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.utils.gradients import autograd
from topfarm import TopFarmProblem
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
site = LillgrundSite()
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.11/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
)
INFO: checking out_of_order...
INFO: out_of_order check complete (0.000371 sec).
INFO: checking system...
INFO: system check complete (0.000011 sec).
INFO: checking solvers...
INFO: solvers check complete (0.000088 sec).
INFO: checking dup_inputs...
INFO: dup_inputs check complete (0.000023 sec).
INFO: checking missing_recorders...
INFO: missing_recorders check complete (0.000005 sec).
INFO: checking unserializable_options...
INFO: unserializable_options check complete (0.000112 sec).
INFO: checking comp_has_no_outputs...
INFO: comp_has_no_outputs check complete (0.000021 sec).
INFO: checking auto_ivc_warnings...
INFO: auto_ivc_warnings check complete (0.000004 sec).
[8]:
_,state,_=problem.optimize()
INFO: checking out_of_order...
INFO: out_of_order check complete (0.000123 sec).
INFO: checking system...
INFO: system check complete (0.000012 sec).
INFO: checking solvers...
INFO: solvers check complete (0.000087 sec).
INFO: checking dup_inputs...
INFO: dup_inputs check complete (0.000019 sec).
INFO: checking missing_recorders...
INFO: missing_recorders check complete (0.000005 sec).
INFO: checking unserializable_options...
INFO: unserializable_options check complete (0.000091 sec).
INFO: checking comp_has_no_outputs...
INFO: comp_has_no_outputs check complete (0.000021 sec).
INFO: checking auto_ivc_warnings...
INFO: auto_ivc_warnings check complete (0.000004 sec).
Optimization terminated successfully (Exit mode 0)
Current function value: -62621.32172858771
Iterations: 3
Function evaluations: 3
Gradient evaluations: 3
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)).plot_wake_map()
plt.ylabel('Height [m]')
plt.xlabel('x [m]')
final hub heights: [ 82. 160. 82. 160. 82. 160. 82.]
[9]:
Text(0.5, 0, 'x [m]')
