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

upd

parent 8960246f
Pipeline #11584 passed with stages
in 11 minutes and 33 seconds
%% Cell type:markdown id: tags:
# Loads
Topfarm supports load constraints by evaluating a trained surrogate model at each iteration. There are currently two different high-level load workflows implemented, namely the Dynamic Wake Meandering (DWM) and the Frandsen approach.
### DWM
In DWM the loads are calculated for all turbines for all flowcases and for all turbine interactions. This means that for each turbine it will look at the wakes coming from all other turbines (individually). The results are subsequently reduced to only one load per turbine by applying a soft max. This is done to create a differentiable expression for the max load instead of just taking the max load, which would not work well in a gradient-based optimization.
### Frandsen
In the Frandsen implementation, the turbulence is aggregated from different sectors with the Wöhler exponent as described in the IEC. Loads are then calculated based on the effective turbulence. Alternatively, if one would like to improve the fidelity by trading off some memory consumption the loads could be calculated based on the all wind direction sectors and subsequently aggregated with the Wöhler exponent.
## Load Surrogates
Topfarm can utilize surrogates trained with a range of different algorithms and softwares, namely scikit-learn, OpenTURNS, TensorFlow and the artificial neural networks code (wind2loads ANN) developed at DTU.
%% Cell type:markdown id: tags:
# Import dependencies
## Import dependencies
%% Cell type:code id: tags:
``` python
import numpy as np
import topfarm
from topfarm import TopFarmProblem
from topfarm.easy_drivers import EasyScipyOptimizeDriver
from topfarm.constraint_components.spacing import SpacingConstraint
from topfarm.constraint_components.boundary import XYBoundaryConstraint
from topfarm.plotting import XYPlotComp
```
%% Cell type:markdown id: tags:
# Install wind2loads
## Install wind2loads
In this exercisw we will use the wind2loads artificial neural networks code to predict the turbine loads of the DTU 10MW reference turbine.
%% Cell type:code id: tags:
``` python
%%capture
try:
!git clone https://gitlab.windenergy.dtu.dk/TOPFARM/workshop-december-2019.git
except:
pass
!pip install -e ./workshop-december-2019
!pip install ./workshop-december-2019/wind2loads
```
%% Cell type:code id: tags:
``` python
import workshop
import os
workshop.w2l_path = 'workshop-december-2019\workshop\wind2loads_ann'
import workshop.workflow as wf
```
%% Output
C:\ProgramData\Anaconda3\lib\site-packages\sklearn\base.py:253: UserWarning: Trying to unpickle estimator MinMaxScaler from version 0.21.3 when using version 0.20.3. This might lead to breaking code or invalid results. Use at your own risk.
UserWarning)
%% Cell type:markdown id: tags:
# Set up the problem
## Set up the problem
In this example we will import most of the code for setting up the workflow.
%% Cell type:code id: tags:
``` python
def get_problem(maxiter=10, allowable_load=1.0, tol=1e-2):
problem = TopFarmProblem(
design_vars={topfarm.x_key: wf.site.initial_position[:, 0],
topfarm.y_key: wf.site.initial_position[:, 1]},
cost_comp=wf.get_load_cost_comp(),
driver=EasyScipyOptimizeDriver(maxiter=maxiter, tol=tol),
constraints=[SpacingConstraint(wf.min_spacing),
XYBoundaryConstraint(wf.boundary),],
post_constraints=[(ls, allowable_load) for ls in wf.load_signals],
plot_comp=XYPlotComp(),
approx_totals={'step':10},
expected_cost=1e-6,)
cost, state = problem.evaluate()
max_load = {ls: np.array([problem[ls+'_abs'].max()]) for ls in wf.load_signals}
problem.cost_comp.analytical_group.lifetime_comp.options['max_load'] = max_load
return problem
```
%% Cell type:markdown id: tags:
# Run the optimization
## Run the optimization
In this example there is the option to change the maximum number of iterations, the tolerance of the optimization as well as the maximum allowable load. The maximum allowable load is set as a percentage of the loads calculated for the initial layout, meaning e.g. that allowable_load=1.0 is 100% of the initial loads. We will optimize on AEP of the wind farm.
%% Cell type:code id: tags:
``` python
problem = get_problem(5, 1.0)
cost, state, recorder = problem.optimize(disp=True)
```
%% Output
Iteration limit exceeded (Exit mode 9)
Current function value: -1285467307.7406654
Iterations: 6
Function evaluations: 6
Gradient evaluations: 6
Optimization FAILED.
Iteration limit exceeded
-----------------------------------
Optimized in 174.461s
Optimized in 130.683s
%% Cell type:markdown id: tags:
# Post Process
## Post Process
During the optimization all variable are recorded. After it is finished one can see all the recorded variables by typing `recorder.keys()` and plot them by writing `plt.plot([your variable],'.')`. We can use this procedure to examine how the load constraint has developed during the optimization:
%% Cell type:markdown id: tags:
First we show the relative loads meaning the ratio between the load at each iteration and the initial loads:
%% Cell type:code id: tags:
``` python
import matplotlib.pyplot as plt
plt.figure()
for ls in wf.load_signals:
plt.plot(np.max(recorder[ls], axis=1),'.',label=ls)
plt.grid()
plt.title('Load constraints')
plt.legend()
plt.ylabel('Load / Initial Load')
plt.xlabel('iteration')
plt.show()
```
%% Output
%% Cell type:markdown id: tags:
Secondly we can also plot the absolute magnitude of the loads as a function of the iterations:
%% Cell type:code id: tags:
``` python
plt.figure()
for ls in wf.load_signals:
ls += '_abs'
plt.plot(np.max(recorder[ls], axis=1),'.',label=ls)
plt.grid()
plt.title('Load constraints')
plt.legend()
plt.ylabel('Load [kNm]')
plt.xlabel('iteration')
plt.show()
```
%% Output
%% Cell type:code id: tags:
``` python
```
......
%% Cell type:markdown id: tags:
# Roads and Cables
%% Cell type:code id: tags:
``` python
%%capture
try:
!git clone https://gitlab.windenergy.dtu.dk/TOPFARM/workshop-december-2019.git
except:
pass
# !cd workshop-december-2019
!pip install -e ./workshop-december-2019
# !cd wind2loads
!pip install ./workshop-december-2019/wind2loads
```
%% Cell type:markdown id: tags:
In colab, use the "inline" backend
%% Cell type:code id: tags:
``` python
# non-updating, inline plots
%matplotlib inline
# ...or updating plots in new window
#%matplotlib qt
```
%% Cell type:markdown id: tags:
Let's import a few classes
%% Cell type:code id: tags:
``` python
import numpy as np
# from openmdao.api import view_model
# from py_wake.examples.data.iea37._iea37 import IEA37_WindTurbines, IEA37Site
from topfarm.cost_models.cost_model_wrappers import CostModelComponent
from topfarm import TopFarmGroup, TopFarmProblem
from topfarm.easy_drivers import EasyRandomSearchDriver, EasyScipyOptimizeDriver
from topfarm.drivers.random_search_driver import RandomizeTurbinePosition_Circle
# from topfarm.constraint_components.boundary import CircleBoundaryConstraint
from topfarm.constraint_components.spacing import SpacingConstraint
from topfarm.constraint_components.boundary import XYBoundaryConstraint
from topfarm.cost_models.electrical.simple_msp import ElNetLength, ElNetCost
import matplotlib.pylab as plt
# Convenient plot component for plotting the cables from simple_msp in jupyter notebooks
from workshop.cabling import XYCablePlotComp, get_site
```
%% Output
---------------------------------------------------------------------------
ModuleNotFoundError Traceback (most recent call last)
<ipython-input-3-50701ef3030d> in <module>
15
16 # Convenient plot component for plotting the cables from simple_msp in jupyter notebooks
---> 17 from workshop.cabling import XYCablePlotComp, get_site
ModuleNotFoundError: No module named 'workshop.cabling'
%% Cell type:markdown id: tags:
### Setting up the site to optimize
We will use the IEA-37 site, using the DTU 10MW reference turbine
%% Cell type:code id: tags:
``` python
from workshop.cabling import get_site
from py_wake.examples.data.DTU10MW_RWT import DTU10MW
site = get_site()
n_wt = len(site.initial_position)
windTurbines = DTU10MW()
Drotor_vector = [windTurbines.diameter()] * n_wt
power_rated_vector = [float(windTurbines.power(20)/1000)] * n_wt
hub_height_vector = [windTurbines.hub_height()] * n_wt
rated_rpm_array = 12. * np.ones([n_wt])
print('Number of turbines:', n_wt)
```
%% Cell type:markdown id: tags:
### Setting up the AEP calculator
- Using the Gaussian wake model from Bastankhah & Porté Agel
- Based on 16 wind direction to speed things up (not critical here because we will be using the RandomSearch algorithm)
%% Cell type:code id: tags:
``` python
from py_wake.wake_models.gaussian import IEA37SimpleBastankhahGaussian
from py_wake.aep_calculator import AEPCalculator
## We use the Gaussian wake model
wake_model = IEA37SimpleBastankhahGaussian(site, windTurbines)
AEPCalc = AEPCalculator(wake_model)
## The AEP is calculated using n_wd wind directions
n_wd = 16
wind_directions = np.linspace(0., 360., n_wd, endpoint=False)
def aep_func(x, y, **kwargs):
"""A simple function that takes as input the x,y position of the turbines and return the AEP"""
return AEPCalc.calculate_AEP(x_i=x, y_i=y, wd=wind_directions).sum(-1).sum(-1)*10**6
```
%% Cell type:markdown id: tags:
### Setting up the NREL IRR cost model
Based on the 2006 NREL report
%% Cell type:code id: tags:
``` python
from topfarm.cost_models.economic_models.turbine_cost import economic_evaluation as EE_NREL
def irr_nrel(aep, electrical_connection_cost, **kwargs):
return EE_NREL(Drotor_vector, power_rated_vector, hub_height_vector, aep, electrical_connection_cost).calculate_irr()
```
%% Cell type:markdown id: tags:
### Setting up the DTU IRR cost model
Based on Witold's recent work
%% Cell type:code id: tags:
``` python
from topfarm.cost_models.economic_models.dtu_wind_cm_main import economic_evaluation as EE_DTU
distance_from_shore = 10.0 # [km]
energy_price = 0.2 / 7.4 # [DKK/kWh] / [DKK/EUR] -> [EUR/kWh]
project_duration = 20 # [years]
water_depth_array = 20 * np.ones([n_wt])
Power_rated_array = np.array(power_rated_vector)/1.0E3 # [MW]
ee_dtu = EE_DTU(distance_from_shore, energy_price, project_duration)
def irr_dtu(aep, electrical_connection_cost, **kwargs):
ee_dtu.calculate_irr(
rated_rpm_array,
Drotor_vector,
Power_rated_array,
hub_height_vector,
water_depth_array,
aep,
electrical_connection_cost)
return ee_dtu.IRR
```
%% Cell type:markdown id: tags:
### Setting up the Topfarm problem
%% Cell type:code id: tags:
``` python
## Some user options
#@markdown Which IRR Cost model to use
IRR_COST = 'DTU' #@param ["DTU", "NREL"]
#@markdown Minimum spacing between the turbines
min_spacing = 2 #@param {type:"slider", min:2, max:10, step:1}
#@markdown Minimum spacing between the turbines
cable_cost_per_meter = 750. # {type:"slider", min:0, max:10000, step:1}
## Electrical grid cable components (Minimum spanning tree from Topfarm report 2010)
elnetlength = ElNetLength(n_wt=n_wt)
elnetcost = ElNetCost(n_wt=n_wt, output_key='electrical_connection_cost', cost_per_meter=cable_cost_per_meter)
# The Topfarm IRR cost model components
irr_dtu_comp = CostModelComponent(input_keys=['aep', ('electrical_connection_cost', 0.0)], n_wt=n_wt,
cost_function=irr_dtu, output_key="irr", output_unit="%", objective=True, income_model=True)
irr_nrel_comp = CostModelComponent(input_keys=['aep', ('electrical_connection_cost', 0.0)], n_wt=n_wt,
cost_function=irr_nrel, output_key="irr", output_unit="%", objective=True, income_model=True)
irr_cost_models = {'DTU': irr_dtu_comp, 'NREL': irr_nrel_comp}
## The Topfarm AEP component
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))
## Plotting component
plot_comp = XYCablePlotComp(memory=0, plot_improvements_only=False, plot_initial=False)
## The group containing all the components
group = TopFarmGroup([aep_comp, elnetlength, elnetcost, irr_cost_models[IRR_COST]])
problem = TopFarmProblem(
design_vars=dict(zip('xy', site.initial_position.T)),
cost_comp=group,
driver=EasyRandomSearchDriver(randomize_func=RandomizeTurbinePosition_Circle(), max_iter=1000),
constraints=[SpacingConstraint(min_spacing * windTurbines.diameter(0)),
XYBoundaryConstraint(site.boundary)],
expected_cost=1.0,
plot_comp=plot_comp)
cost, state, recorder = problem.optimize()
```
%% Cell type:markdown id: tags:
### Exercises
- Try to see what is the effect of increasing or decreasing the cost of the cable
- Change between IRR cost model. Ask Witold about the difference between DTU and NREL models
%% Cell type:code id: tags:
``` python
```
......
%% Cell type:markdown id: tags:
# Loads
Topfarm supports load constraints by evaluating a trained surrogate model at each iteration. There are currently two different high-level load workflows implemented, namely the Dynamic Wake Meandering (DWM) and the Frandsen approach.
### DWM
In DWM the loads are calculated for all turbines for all flowcases and for all turbine interactions. This means that for each turbine it will look at the wakes coming from all other turbines (individually). The results are subsequently reduced to only one load per turbine by applying a soft max. This is done to create a differentiable expression for the max load instead of just taking the max load, which would not work well in a gradient-based optimization.
### Frandsen
In the Frandsen implementation, the turbulence is aggregated from different sectors with the Wöhler exponent as described in the IEC. Loads are then calculated based on the effective turbulence. Alternatively, if one would like to improve the fidelity by trading off some memory consumption the loads could be calculated based on the all wind direction sectors and subsequently aggregated with the Wöhler exponent.
## Load Surrogates
Topfarm can utilize surrogates trained with a range of different algorithms and softwares, namely scikit-learn, OpenTURNS, TensorFlow and the artificial neural networks code (wind2loads ANN) developed at DTU.
%% Cell type:markdown id: tags:
[Try this yourself](https://colab.research.google.com/github/DTUWindEnergy/TopFarm2/blob/master/docs/notebooks/loads.ipynb) (requires google account)
%% Cell type:code id: tags:
``` python
%%capture
# Install Topfarm if needed
try:
import topfarm
except ModuleNotFoundError:
!pip install topfarm
```
%% Cell type:markdown id: tags:
# Import dependencies
## Import dependencies
%% Cell type:code id: tags:
``` python
import numpy as np
import topfarm
from topfarm import TopFarmProblem
from topfarm.easy_drivers import EasyScipyOptimizeDriver
from topfarm.constraint_components.spacing import SpacingConstraint
from topfarm.constraint_components.boundary import XYBoundaryConstraint
from topfarm.plotting import XYPlotComp
```
%% Cell type:markdown id: tags:
# Install wind2loads
## Install wind2loads
In this exercisw we will use the wind2loads artificial neural networks code to predict the turbine loads of the DTU 10MW reference turbine.
%% Cell type:code id: tags:
``` python
%%capture
try:
!git clone https://gitlab.windenergy.dtu.dk/TOPFARM/workshop-december-2019.git
except:
pass
!pip install -e ./workshop-december-2019
!pip install ./workshop-december-2019/wind2loads
```
%% Cell type:code id: tags:
``` python
import workshop
import os
workshop.w2l_path = 'workshop-december-2019\workshop\wind2loads_ann'
import workshop.workflow as wf
```
%% Output
C:\ProgramData\Anaconda3\lib\site-packages\sklearn\base.py:253: UserWarning: Trying to unpickle estimator MinMaxScaler from version 0.21.3 when using version 0.20.3. This might lead to breaking code or invalid results. Use at your own risk.
UserWarning)
%% Cell type:markdown id: tags:
# Set up the problem
## Set up the problem
In this example we will import most of the code for setting up the workflow.
%% Cell type:code id: tags:
``` python
def get_problem(maxiter=10, allowable_load=1.0, tol=1e-2):
problem = TopFarmProblem(
design_vars={topfarm.x_key: wf.site.initial_position[:, 0],
topfarm.y_key: wf.site.initial_position[:, 1]},
cost_comp=wf.get_load_cost_comp(),
driver=EasyScipyOptimizeDriver(maxiter=maxiter, tol=tol),
constraints=[SpacingConstraint(wf.min_spacing),
XYBoundaryConstraint(wf.boundary),],
post_constraints=[(ls, allowable_load) for ls in wf.load_signals],
plot_comp=XYPlotComp(),
approx_totals={'step':10},
expected_cost=1e-6,)
cost, state = problem.evaluate()
max_load = {ls: np.array([problem[ls+'_abs'].max()]) for ls in wf.load_signals}
problem.cost_comp.analytical_group.lifetime_comp.options['max_load'] = max_load
return problem
```
%% Cell type:markdown id: tags:
# Run the optimization
## Run the optimization
In this example there is the option to change the maximum number of iterations, the tolerance of the optimization as well as the maximum allowable load. The maximum allowable load is set as a percentage of the loads calculated for the initial layout, meaning e.g. that allowable_load=1.0 is 100% of the initial loads. We will optimize on AEP of the wind farm.
%% Cell type:code id: tags:
``` python
problem = get_problem(5, 1.0)
cost, state, recorder = problem.optimize(disp=True)
```
%% Output
Iteration limit exceeded (Exit mode 9)
Current function value: -1285467307.7406654
Iterations: 6
Function evaluations: 6
Gradient evaluations: 6
Optimization FAILED.
Iteration limit exceeded
-----------------------------------
Optimized in 174.461s
Optimized in 130.683s
%% Cell type:markdown id: tags:
# Post Process
## Post Process
During the optimization all variable are recorded. After it is finished one can see all the recorded variables by typing `recorder.keys()` and plot them by writing `plt.plot([your variable],'.')`. We can use this procedure to examine how the load constraint has developed during the optimization:
%% Cell type:markdown id: tags:
First we show the relative loads meaning the ratio between the load at each iteration and the initial loads:
%% Cell type:code id: tags:
``` python
import matplotlib.pyplot as plt
plt.figure()
for ls in wf.load_signals:
plt.plot(np.max(recorder[ls], axis=1),'.',label=ls)
plt.grid()
plt.title('Load constraints')
plt.legend()
plt.ylabel('Load / Initial Load')
plt.xlabel('iteration')
plt.show()
```
%% Output
%% Cell type:markdown id: tags:
Secondly we can also plot the absolute magnitude of the loads as a function of the iterations:
%% Cell type:code id: tags:
``` python
plt.figure()
for ls in wf.load_signals:
ls += '_abs'
plt.plot(np.max(recorder[ls], axis=1),'.',label=ls)
plt.grid()
plt.title('Load constraints')
plt.legend()
plt.ylabel('Load [kNm]')
plt.xlabel('iteration')
plt.show()
```
%% Output