Commit cba70b20 authored by Mikkel Friis-Møller's avatar Mikkel Friis-Møller
Browse files

added bathymetry.ipynb

parent bb02eddc
Pipeline #25359 passed with stages
in 14 minutes and 11 seconds
%% Cell type:markdown id: tags:
# Optimization with bathymetry and max water depth constraint
%% Cell type:markdown id: tags:
## Install packages if running in Colab
%% Cell type:code id: tags:
``` python
try:
RunningInCOLAB = 'google.colab' in str(get_ipython())
except NameError:
RunningInCOLAB = False
```
%% Cell type:code id: tags:
``` python
%%capture
if RunningInCOLAB:
!pip install git+https://gitlab.windenergy.dtu.dk/TOPFARM/PyWake.git@grd_files
!pip install git+https://gitlab.windenergy.dtu.dk/TOPFARM/TopFarm2.git@bathymetry
!pip install git+https://gitlab.windenergy.dtu.dk/TOPFARM/PyWake.git
!pip install git+https://gitlab.windenergy.dtu.dk/TOPFARM/TopFarm2.git
!pip install scipy==1.6.3 # constraint is not continuous which trips vers. 1.4.1 which presently is the default version
import os
os.kill(os.getpid(), 9)
```
%% Cell type:markdown id: tags:
## Import section
%% Cell type:code id: tags:
``` python
import numpy as np
import matplotlib.pyplot as plt
import time
from topfarm.cost_models.cost_model_wrappers import CostModelComponent
from topfarm.easy_drivers import EasyScipyOptimizeDriver
from topfarm import TopFarmProblem
from topfarm.plotting import NoPlot, XYPlotComp
from topfarm.constraint_components.boundary import XYBoundaryConstraint
from topfarm.constraint_components.spacing import SpacingConstraint
from topfarm.examples.data.parque_ficticio_offshore import ParqueFicticioOffshore
from py_wake.deficit_models.gaussian import IEA37SimpleBastankhahGaussian
from py_wake.examples.data.iea37._iea37 import IEA37_WindTurbines
```
%% Cell type:markdown id: tags:
## Set up site and optimization problem
%% Cell type:code id: tags:
``` python
site = ParqueFicticioOffshore()
site.bounds = 'ignore'
x_init, y_init = site.initial_position[:,0], site.initial_position[:,1]
boundary = site.boundary
# # # Wind turbines and wind farm model definition
windTurbines = IEA37_WindTurbines()
wfm = IEA37SimpleBastankhahGaussian(site, windTurbines)
wsp = np.asarray([10, 15])
wdir = np.arange(0,360,45)
maximum_water_depth = -52
n_wt = x_init.size
maxiter = 10
def aep_func(x, y, **kwargs):
simres = wfm(x, y, wd=wdir, ws=wsp)
aep = simres.aep().values.sum()
water_depth = np.diag(wfm.site.ds.interp(x=x, y=y)['water_depth'])
return [aep, water_depth]
tol = 1e-8
ec = 1e-2
min_spacing = 260
cost_comp = CostModelComponent(input_keys=[('x', x_init),('y', y_init)],
n_wt=n_wt,
cost_function=aep_func,
objective=True,
maximize=True,
output_keys=[('AEP', 0), ('water_depth', np.zeros(n_wt))]
)
problem = TopFarmProblem(design_vars={'x': x_init, 'y': y_init},
constraints=[XYBoundaryConstraint(boundary),
SpacingConstraint(min_spacing)],
post_constraints=[('water_depth', {'lower': maximum_water_depth})],
cost_comp=cost_comp,
driver=EasyScipyOptimizeDriver(optimizer='SLSQP', maxiter=maxiter, tol=tol),
plot_comp=XYPlotComp(),
expected_cost=ec)
```
%% Cell type:markdown id: tags:
## Optimize
%% Cell type:code id: tags:
``` python
tic = time.time()
cost, state, recorder = problem.optimize()
toc = time.time()
print('Optimization took: {:.0f}s'.format(toc-tic))
```
%% Output
Iteration limit reached (Exit mode 9)
Current function value: [-20412.33040188]
Iterations: 10
Function evaluations: 10
Gradient evaluations: 10
Optimization FAILED.
Iteration limit reached
-----------------------------------
Optimization took: 20s
%% Cell type:markdown id: tags:
## Check the max water depth
%% Cell type:code id: tags:
``` python
# plt.plot(recorder['water_depth'].min((1)))
# plt.plot([0,recorder['water_depth'].shape[0]],[maximum_water_depth, maximum_water_depth])
# plt.xlabel('Iteration')
# plt.ylabel('Max depth [m]')
```
%% Output
%% Cell type:markdown id: tags:
## Check the initial- and optimized layout wrt. the water depth boundary
%% Cell type:code id: tags:
``` python
# values = site.ds.water_depth.values
# x = site.ds.x.values
# y = site.ds.y.values
# levels = np.arange(int(values.min()), int(values.max()))
# max_wd_index = int(np.argwhere(levels==maximum_water_depth))
# Y, X = np.meshgrid(y, x)
# cs = plt.contour(x, y , values.T, levels)
# plt.close()
# lines = []
# for line in cs.collections[max_wd_index].get_paths():
# lines.append(line.vertices)
# fig2, ax2 = plt.subplots(1)
# site.ds.water_depth.plot(ax=ax2, levels=100)
# for line in lines:
# ax2.plot(line[:, 0], line[:,1], 'r')
# problem.model.plot_comp.plot_initial2current(x_init, y_init, state['x'], state['y'])
# ax2.set_title(f'Max Water Depth Boundary: {maximum_water_depth} m')
```
%% Output
Text(0.5, 1.0, 'Max Water Depth Boundary: -52 m')
%% Cell type:code id: tags:
``` python
```
......
......@@ -61,10 +61,11 @@ def make_doc_notebooks(notebooks):
t = '[Try this yourself](https://colab.research.google.com/github/DTUWindEnergy/TopFarm2/blob/master/docs/notebooks/%s.ipynb) (requires google account)'
nb.insert_markdown_cell(1, t % name)
code = """%%capture
# Install Topfarm if needed
import importlib
if not importlib.util.find_spec("topfarm"):
!pip install topfarm
# Install Topfarm if needed"""
"""
if not name in ['loads', 'wake_steering_and_loads', 'layout_and_loads']:
nb.insert_code_cell(2, code)
nb.save(dst_path + name + ".ipynb")
......
%% Cell type:markdown id: tags:
# Optimization with bathymetry and max water depth constraint
%% Cell type:markdown id: tags:
[Try this yourself](https://colab.research.google.com/github/DTUWindEnergy/TopFarm2/blob/master/docs/notebooks/bathymetry.ipynb) (requires google account)
%% Cell type:code id: tags:
``` python
%%capture\nimport importlib\nif not importlib.util.find_spec("topfarm"):\n !pip install topfarm\n# Install Topfarm if needed
%%capture\n# Install Topfarm if needed\nimport importlib\nif not importlib.util.find_spec("topfarm"):\n !pip install topfarm\n
```
%% Cell type:markdown id: tags:
## Install packages if running in Colab
%% Cell type:code id: tags:
``` python
try:
RunningInCOLAB = 'google.colab' in str(get_ipython())
except NameError:
RunningInCOLAB = False
```
%% Cell type:code id: tags:
``` python
%%capture
if RunningInCOLAB:
!pip install git+https://gitlab.windenergy.dtu.dk/TOPFARM/PyWake.git@grd_files
!pip install git+https://gitlab.windenergy.dtu.dk/TOPFARM/TopFarm2.git@bathymetry
!pip install git+https://gitlab.windenergy.dtu.dk/TOPFARM/PyWake.git
!pip install git+https://gitlab.windenergy.dtu.dk/TOPFARM/TopFarm2.git
!pip install scipy==1.6.3 # constraint is not continuous which trips vers. 1.4.1 which presently is the default version
import os
os.kill(os.getpid(), 9)
```
%% Cell type:markdown id: tags:
## Import section
%% Cell type:code id: tags:
``` python
import numpy as np
import matplotlib.pyplot as plt
import time
from topfarm.cost_models.cost_model_wrappers import CostModelComponent
from topfarm.easy_drivers import EasyScipyOptimizeDriver
from topfarm import TopFarmProblem
from topfarm.plotting import NoPlot, XYPlotComp
from topfarm.constraint_components.boundary import XYBoundaryConstraint
from topfarm.constraint_components.spacing import SpacingConstraint
from topfarm.examples.data.parque_ficticio_offshore import ParqueFicticioOffshore
from py_wake.deficit_models.gaussian import IEA37SimpleBastankhahGaussian
from py_wake.examples.data.iea37._iea37 import IEA37_WindTurbines
```
%% Cell type:markdown id: tags:
## Set up site and optimization problem
%% Cell type:code id: tags:
``` python
site = ParqueFicticioOffshore()
site.bounds = 'ignore'
x_init, y_init = site.initial_position[:,0], site.initial_position[:,1]
boundary = site.boundary
# # # Wind turbines and wind farm model definition
windTurbines = IEA37_WindTurbines()
wfm = IEA37SimpleBastankhahGaussian(site, windTurbines)
wsp = np.asarray([10, 15])
wdir = np.arange(0,360,45)
maximum_water_depth = -52
n_wt = x_init.size
maxiter = 10
def aep_func(x, y, **kwargs):
simres = wfm(x, y, wd=wdir, ws=wsp)
aep = simres.aep().values.sum()
water_depth = np.diag(wfm.site.ds.interp(x=x, y=y)['water_depth'])
return [aep, water_depth]
tol = 1e-8
ec = 1e-2
min_spacing = 260
cost_comp = CostModelComponent(input_keys=[('x', x_init),('y', y_init)],
n_wt=n_wt,
cost_function=aep_func,
objective=True,
maximize=True,
output_keys=[('AEP', 0), ('water_depth', np.zeros(n_wt))]
)
problem = TopFarmProblem(design_vars={'x': x_init, 'y': y_init},
constraints=[XYBoundaryConstraint(boundary),
SpacingConstraint(min_spacing)],
post_constraints=[('water_depth', {'lower': maximum_water_depth})],
cost_comp=cost_comp,
driver=EasyScipyOptimizeDriver(optimizer='SLSQP', maxiter=maxiter, tol=tol),
plot_comp=XYPlotComp(),
expected_cost=ec)
```
%% Cell type:markdown id: tags:
## Optimize
%% Cell type:code id: tags:
``` python
tic = time.time()
cost, state, recorder = problem.optimize()
toc = time.time()
print('Optimization took: {:.0f}s'.format(toc-tic))
```
%% Output
Iteration limit reached (Exit mode 9)
Current function value: [-20412.33040188]
Iterations: 10
Function evaluations: 10
Gradient evaluations: 10
Optimization FAILED.
Iteration limit reached
-----------------------------------
Optimization took: 20s
%% Cell type:markdown id: tags:
## Check the max water depth
%% Cell type:code id: tags:
``` python
# plt.plot(recorder['water_depth'].min((1)))
# plt.plot([0,recorder['water_depth'].shape[0]],[maximum_water_depth, maximum_water_depth])
# plt.xlabel('Iteration')
# plt.ylabel('Max depth [m]')
```
%% Output
%% Cell type:markdown id: tags:
## Check the initial- and optimized layout wrt. the water depth boundary
%% Cell type:code id: tags:
``` python
# values = site.ds.water_depth.values
# x = site.ds.x.values
# y = site.ds.y.values
# levels = np.arange(int(values.min()), int(values.max()))
# max_wd_index = int(np.argwhere(levels==maximum_water_depth))
# Y, X = np.meshgrid(y, x)
# cs = plt.contour(x, y , values.T, levels)
# plt.close()
# lines = []
# for line in cs.collections[max_wd_index].get_paths():
# lines.append(line.vertices)
# fig2, ax2 = plt.subplots(1)
# site.ds.water_depth.plot(ax=ax2, levels=100)
# for line in lines:
# ax2.plot(line[:, 0], line[:,1], 'r')
# problem.model.plot_comp.plot_initial2current(x_init, y_init, state['x'], state['y'])
# ax2.set_title(f'Max Water Depth Boundary: {maximum_water_depth} m')
```
%% Output
Text(0.5, 1.0, 'Max Water Depth Boundary: -52 m')
%% Cell type:code id: tags:
``` python
```
......
%% Cell type:markdown id: tags:
# Constraints
%% Cell type:markdown id: tags:
[Try this yourself](https://colab.research.google.com/github/DTUWindEnergy/TopFarm2/blob/master/docs/notebooks/constraints.ipynb) (requires google account)
%% Cell type:code id: tags:
``` python
%%capture\nimport importlib\nif not importlib.util.find_spec("topfarm"):\n !pip install topfarm\n# Install Topfarm if needed
%%capture\n# Install Topfarm if needed\nimport importlib\nif not importlib.util.find_spec("topfarm"):\n !pip install topfarm\n
```
%% Cell type:markdown id: tags:
Constraints are the second key element of a an optimization problem formulation. They ensure that the optimization results conforms to feasible / realistic solutions. There are three types of constraints in optimization:
* Variable bounds - upper and lower boundary values for design variables
* inequality constraints - constraint function values must be less or more than a given threshold
* equality constraints - constraint function must be exactly equal to a value (not as commonly used)
This notebook walks through a process to set up typical constraints in Topfarm for wind farm design problems including:
* boundary constraints - these tell Topfarm within a region where the perimeter of the site is that turbines must be sited within and also any exclusion zones that must be avoided
* spacing constraints - these tell Topfarm the minimum allowable inter-turbine spacing in the farm
%% Cell type:markdown id: tags:
**First we import supporting libraries in Python numpy and matplotlib**
%% Cell type:code id: tags:
``` python
# Import numpy and matplotlib files
import numpy as np
%matplotlib inline
import matplotlib.pyplot as plt
```
%% Cell type:markdown id: tags:
**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**
* **XYBoundaryConstraint - for a boundary specified as a series of connected perimeter vertices**
* **CircleBoundaryConstraint - for a circular boundary with a central location and a radius**
* **SpacingConstraint - for the inter-turbine spacing distance constraints**
* **CostModelComponent - a generic class for setting up a problem objective function**
**We also import a helper function to plot**
**For documentation on Topfarm see:**
https://topfarm.pages.windenergy.dtu.dk/TopFarm2/user_guide.html#
%% Cell type:code id: tags:
``` python
# Import topfarm problem, plotting support, constraint classes and generic cost model component
from topfarm import TopFarmProblem
from topfarm.plotting import XYPlotComp
from topfarm.constraint_components.boundary import XYBoundaryConstraint, CircleBoundaryConstraint
from topfarm.constraint_components.spacing import SpacingConstraint
from topfarm.cost_models.cost_model_wrappers import CostModelComponent
```
%% Cell type:markdown id: tags:
## Boundary constraints
**Next we are going to demonstrate the use of the XYBoundaryConstraint to set up site boundaries of a variety of types including square, rectangle and an arbitrary polygon. Additionally, a "convex hull" example is provided which is a commonly used boundary type for wind farm design optimization problems.**
**(From wikipedia) In mathematics, the convex hull or convex envelope or convex closure of a set X of points in the Euclidean plane or in a Euclidean space (or, more generally, in an affine space over the reals) is the smallest convex set that contains X. For instance, when X is a bounded subset of the plane, the convex hull may be visualized as the shape enclosed by a rubber band stretched around X.**
%% Cell type:code id: tags:
``` python
# set up a "boundary" arry with arbitrary points for use in the example
boundary = np.array([(0, 0), (1, 1), (3, 0), (3, 2), (0, 2)])
# set up dummy design variables and cost model component.
# This example includes 2 turbines (n_wt=2) located at x,y=0.5,0.5 and 1.5,1.5 respectively
x = [0.5,1.5]
y = [.5,1.5]
dummy_cost = CostModelComponent(input_keys=[],
n_wt=2,
cost_function=lambda : 1)
# We introduce a simple plotting function so we can quickly plot different types of site boundaries
def plot_boundary(name, constraint_comp):
tf = TopFarmProblem(
design_vars={'x':x, 'y':y}, # setting up our two turbines as design variables
cost_comp=dummy_cost, # using dummy cost model
constraints=[constraint_comp], # constraint set up for the boundary type provided
plot_comp=XYPlotComp()) # support plotting function
plt.figure()
plt.title(name)
tf.plot_comp.plot_constraints() # plot constraints is a helper function in topfarm to plot constraints
plt.plot(boundary[:,0], boundary[:,1],'.r', label='Boundary points') # plot the boundary points
plt.axis('equal')
plt.legend() # add the legend
```
%% Cell type:markdown id: tags:
**Now that we have set up our dummy problem, we can illustrate how different boundary types can be created from our boundary vertices**
**First we show a convex hull type as described above. Note that for the convex hull, all boundary points are contained within a convex perimeter but one of the boundary points on the interior is not used.**
%% Cell type:code id: tags:
``` python
# Now we use our boundary defined above to plot different types of boundary constraints on the problem
plot_boundary('convex_hull', XYBoundaryConstraint(boundary, 'convex_hull'))
```
%% Output
%% Cell type:markdown id: tags:
**Next we show a square type of boundary. In this case the maximum distance between the x and y elements of the vertices is used to establish the perimeter.**
%% Cell type:code id: tags:
``` python
plot_boundary('square', XYBoundaryConstraint(boundary, 'square'))
```
%% Output
%% Cell type:markdown id: tags:
**Now a rectangle boundary. Here we use the maximum distance on both x and y axes of the boundary coordinates to establish the perimeter.**
%% Cell type:code id: tags:
``` python
plot_boundary('rectangle', XYBoundaryConstraint(boundary, 'rectangle'))
```
%% Output
%% Cell type:markdown id: tags:
**Now a polygon boundary. This connects all the points in sequence. Note that this results in a nonconvex boundary. Nonconvex functions in optimization problems introduce complexity that can be challenging to handle and often require more sophisticated algorithms and higher computational expense.**
%% Cell type:code id: tags:
``` python
plot_boundary('polygon', XYBoundaryConstraint(boundary, 'polygon'))
```
%% Output
%% Cell type:markdown id: tags:
**Finally a circular boundary. For this we have to specify the midpoint of the circle and the radius.**
%% Cell type:code id: tags:
``` python
plot_boundary('Circle',CircleBoundaryConstraint((1.5,1),1))
```
%% Output
%% Cell type:markdown id: tags:
**You can also quickly plot all boundary types using some quick python tricks**
%% Cell type:code id: tags:
``` python
# plot all boundary types using our dummy problem and boundaries
for boundary_type in ['convex_hull','square','rectangle','polygon']:
plot_boundary(boundary_type, XYBoundaryConstraint(boundary, boundary_type))
```
%% Output
%% Cell type:markdown id: tags:
### Exercise!!
**Now play around with a new set of boundary vertices and construct different perimeters to explore the functionality. See if you can make even more complex polygon shapes.**
%% Cell type:code id: tags:
``` python
# Make your own set of vertices - anything you would like!
boundary = np.array([(0, 0), (2, 1), (4, 7), (1, 1), (0, 2)])
# Then see what types of perimeters they generate
for boundary_type in ['convex_hull','square','rectangle','polygon']:
plot_boundary(boundary_type, XYBoundaryConstraint(boundary, boundary_type))
```
%% Output
%% Cell type:markdown id: tags:
## Spacing constraints
The next most common constraint in a wind farm design optimization problem is on the allowable inter-turbine spacing in the farm. Losses due to wakes should ensure that turbines do spread out but a minimum constraint can also help to ensure that turbines do not get placed too close together.
The following provides a simple example of implementation of a minimum spacing constraint.
%% Cell type:markdown id: tags:
**Implement a minimum spacing constraint example**
%% Cell type:code id: tags:
``` python
# set up dummy design variables and cost model component.
# This example includes 2 turbines (n_wt=2) located at x,y=0.5,0.5 and 1.5,1.5 respectively
x = [0.5,1.5]
y = [.5,1.5]
dummy_cost = CostModelComponent(input_keys=[],
n_wt=2,
cost_function=lambda : 1)
# a function to plot a spacing constraint for a Topfarm problem
def plot_spacing(name, constraint_comp):
tf = TopFarmProblem(
design_vars={'x':x, 'y':y}, # setting up our two turbines as design variables
cost_comp=dummy_cost, # using dummy cost model
constraints=[constraint_comp], # constraint set up for the boundary type provided
plot_comp=XYPlotComp()) # support plotting function
tf.evaluate()
plt.figure()
plt.title(name)
tf.plot_comp.plot_constraints() # plot constraints is a helper function in topfarm to plot constraints
plt.plot(x,y,'.b', label='Wind turbines') # plot the turbine locations
plt.axis('equal')
plt.legend() # add the legend
plt.ylim([0,2])
plot_spacing('spacing', SpacingConstraint(1))
```
%% Output
%% Cell type:markdown id: tags:
### Exercise!!
**Play around with the spacing constraint size**
%% Cell type:code id: tags:
``` python
plot_spacing('spacing', SpacingConstraint(3))
```
%% Output
%% Cell type:code id: tags:
``` python
```
......
%% Cell type:markdown id: tags:
## Cost Models
Topfarm now comes with two built-in cost models. Additional user-defined cost models can easily be integrated as well. In this example, the ability to switch from using each of the two models is demonstrated.
%% Cell type:markdown id: tags:
[Try this yourself](https://colab.research.google.com/github/DTUWindEnergy/TopFarm2/blob/master/docs/notebooks/cost_models.ipynb) (requires google account)
%% Cell type:code id: tags:
``` python
%%capture\nimport importlib\nif not importlib.util.find_spec("topfarm"):\n !pip install topfarm\n# Install Topfarm if needed
%%capture\n# Install Topfarm if needed\nimport importlib\nif not importlib.util.find_spec("topfarm"):\n !pip install topfarm\n
```
%% Cell type:markdown id: tags:
### Cost Model 1: DTU implementation of the NREL Cost and Scaling Model
The first cost model in Topfarm is a python implementation of the National Renewable Energy Laboratory (NREL) Cost and Scaling Model. The report on which the model is based can be found here:
https://www.nrel.gov/docs/fy07osti/40566.pdf
The model was developed from the early to mid-2000s as part of the Wind Partnership for Advanced Component Technology (WindPACT) which explored innovative turbine design (at that time) as well as innovations on the balance of plant and operations. Several detailed design studies on the turbine and plant design can cost were made. For publications associated with the WindPACT program, see: [WindPACT publication list](http://nrel-primo.hosted.exlibrisgroup.com/primo_library/libweb/action/search.do;jsessionid=00F1EA4B14428BED000D0D1E7C0E2C46?fn=search&ct=search&initialSearch=true&mode=Basic&tab=default_tab&indx=1&dum=true&srt=rank&vid=Pubs&frbg=&vl%28freeText0%29=windpact&scp.scps=scope%3A%28PUBS%29%2Cscope%3A%28NREL_INTERNAL%29&vl%28870446075UI1%29=all_items)
From the WindPACT studies, the NREL cost and scaling model was developed as a set of curve-fits to the underlying detailed design data and includes:
* Turbine component masses and costs
* Balance of system costs
* Operational expenditures
* Financing and other costs
Over time, changes in turbine and plant technology have rendered the NREL cost and scaling model obselete, but it is still useful as a publicly available, full levelized cost of energy (LCOE) model for wind energy.
%% Cell type:markdown id: tags:
**Import Topfarm models to set up an LCOE workflow including the cost model**
%% Cell type:code id: tags:
``` python
# Import numerical python
import numpy as np
# Import pywake models including the IEA Wind Task 37 case study site, the Gaussian wake model and the AEP calculator
from py_wake.examples.data.iea37._iea37 import IEA37_WindTurbines, IEA37Site
from py_wake.deficit_models.gaussian import IEA37SimpleBastankhahGaussian
# Import Topfarm implementation of NREL Cost and Scaling model
from topfarm.cost_models.economic_models.turbine_cost import economic_evaluation as ee_1
# Import Topfarm constraints for site boundary and spacing
from topfarm.drivers.random_search_driver import RandomizeTurbinePosition_Circle
from topfarm.constraint_components.boundary import CircleBoundaryConstraint
from topfarm.constraint_components.spacing import SpacingConstraint
# Import Topfarm support classes for setting up problem and workflow
from topfarm.cost_models.cost_model_wrappers import CostModelComponent, AEPCostModelComponent
from topfarm.cost_models.py_wake_wrapper import PyWakeAEPCostModelComponent
from topfarm import TopFarmGroup, TopFarmProblem
from topfarm.plotting import XYPlotComp, NoPlot
# Import Topfarm implementation of Random Search or Scipy drivers
from topfarm.easy_drivers import EasyRandomSearchDriver
from topfarm.easy_drivers import EasyScipyOptimizeDriver
from topfarm.easy_drivers import EasySimpleGADriver
```
%% Cell type:markdown id: tags:
**Set up plotting capability**
%% Cell type:code id: tags:
``` python
try:
import matplotlib.pyplot as plt
plt.gcf()
plot_comp = XYPlotComp()
plot = True
except RuntimeError:
plot_comp = NoPlot()
plot = False
```
%% Output
%% Cell type:markdown id: tags:
**Set up IEA Wind Task 37 case study site with 16 turbines.**
%% Cell type:code id: tags:
``` python
# site set up
n_wt = 16 # number of wind turbines
site = IEA37Site(n_wt) # site is the IEA Wind Task 37 site with a circle boundary
windTurbines = IEA37_WindTurbines() # wind turbines are the IEA Wind Task 37 3.4 MW reference turbine
wake_model = IEA37SimpleBastankhahGaussian(site, windTurbines) # select the Gaussian wake model
# vectors for turbine properties: diameter, rated power and hub height
# these are inputs to the cost model
Drotor_vector = [windTurbines.diameter()] * n_wt
power_rated_vector = [float(windTurbines.power(20)/1000)] * n_wt
hub_height_vector = [windTurbines.hub_height()] * n_wt
```
%% Cell type:markdown id: tags:
**Set up functions for the AEP and cost calculations. Here we are using the internal rate of return (IRR) as our financial metric of interest.**
%% Cell type:code id: tags:
``` python
# function for calculating aep as a function of x,y positions of the wind turbiens
def aep_func(x, y, **kwargs):
return wake_model(x, y).aep().sum(['wd','ws']).values*10**6
# function for calculating overall internal rate of return (IRR)
def irr_func(aep, **kwargs):
my_irr = ee_1(Drotor_vector, power_rated_vector, hub_height_vector, aep).calculate_irr()
print(my_irr)
return my_irr
```
%% Cell type:markdown id: tags:
**Now set up a problem to run an optimization using IRR as the objective function. Note that the turbines are fixed so the main driver changing the IRR will be the AEP as the turbine positions change.**
%% Cell type:code id: tags:
``` python
# create an openmdao component for aep and irr to add to the problem
aep_comp = CostModelComponent(input_keys=['x','y'], n_wt=n_wt, cost_function=aep_func, output_key="aep", output_unit="GWh", objective=False, output_val=np.zeros(n_wt))
irr_comp = CostModelComponent(input_keys=['aep'], n_wt=n_wt, cost_function=irr_func, output_key="irr", output_unit="%", objective=True, income_model=True)
# create a group for the aep and irr components that links their common input/output (aep)
irr_group = TopFarmGroup([aep_comp, irr_comp])
# add the group to an optimization problem and specify the design variables (turbine positions),
# cost function (irr_group), driver (random search), and constraints (circular boundary and spacing)
problem = TopFarmProblem(
design_vars=dict(zip('xy', site.initial_position.T)),
cost_comp=irr_group,
driver=EasyRandomSearchDriver(randomize_func=RandomizeTurbinePosition_Circle(), max_iter=50),
#driver=EasyScipyOptimizeDriver(optimizer='COBYLA', maxiter=200, tol=1e-6, disp=False),
#driver=EasySimpleGADriver(max_gen=100, pop_size=5, Pm=None, Pc=.5, elitism=True, bits={}),
constraints=[SpacingConstraint(200),
CircleBoundaryConstraint([0, 0], 1300.1)],
plot_comp=plot_comp)
# assign data from optimizationn to a set of accessible variables and run the optimization
cost, state, recorder = problem.optimize()
```
%% Output
59.15317035889765
59.15317035889765
59.15317035889765
59.14310441797548
59.15317035889765
59.15317035889765
59.15317035889765
59.15317035889765
59.15317035889765
58.29211938020884
59.15317035889765
59.15317035889765
58.306611078551015
59.15317035889765
58.65177988998911
59.15317035889765
59.15317035889765
59.15317035889765
59.15317035889765
59.15317035889765
58.59176040816389
59.15317035889765
59.31620110565699
59.31620110565699
59.31620110565699
59.31620110565699
59.31620110565699
58.66292698793363
59.31620110565699
59.31620110565699
59.31620110565699
59.31620110565699
59.31620110565699
59.31620110565699
59.31620110565699
59.303257188319215
59.31620110565699
59.31620110565699
57.47298599965627
59.31620110565699
59.31620110565699
59.31620110565699
59.31620110565699
59.31620110565699
57.6380941604278
59.31620110565699
59.31620110565699
59.31620110565699
59.31620110565699
58.43798216