Skip to content
GitLab
Commit 4a7886d4 authored by Mikkel Friis-Møller's avatar Mikkel Friis-Møller
Browse files

npv calc also moved to np financial

parent 4792703e
Pipeline #23991 passed with stage
in 13 minutes and 41 seconds
Expand all Hide whitespace changes
Inline Side-by-side
.gitlab-ci.yml
View file @ 4a7886d4
......@@ -52,7 +52,7 @@ pages: # "pages" is a job specifically for GitLab pages [1]
only: # only run for these branches
- master
- /^test_doc.*/
- opt_course
tags: # only runners with this tag can do the job [3]
- python
......
_notebooks/elements/layout_and_loads.ipynb
View file @ 4a7886d4
%% Cell type:markdown id: tags:
# Load Constrained Layout Optimization
%% Cell type:markdown id: tags:
## Install TopFarm and PyWake
%% Cell type:code id: tags:
``` python
%%capture
try:
import py_wake
except:
!pip install git+https://gitlab.windenergy.dtu.dk/TOPFARM/PyWake.git@setup_for_surrogate
!pip install git+https://gitlab.windenergy.dtu.dk/TOPFARM/PyWake.git
try:
import topfarm
except:
!pip install git+https://gitlab.windenergy.dtu.dk/TOPFARM/TopFarm2.git@opt_course
!pip install topfarm
```
%% Cell type:markdown id: tags:
## Import section
%% Cell type:code id: tags:
``` python
import numpy as np
from numpy import newaxis as na
import time
from topfarm.cost_models.cost_model_wrappers import AEPMaxLoadCostModelComponent
from topfarm.easy_drivers import EasyScipyOptimizeDriver
from topfarm import TopFarmProblem
from topfarm.plotting import NoPlot
from topfarm.constraint_components.boundary import XYBoundaryConstraint
from topfarm.constraint_components.spacing import SpacingConstraint
from py_wake.examples.data.lillgrund import LillgrundSite
from py_wake.deficit_models.gaussian import IEA37SimpleBastankhahGaussian, NiayifarGaussian
from py_wake.turbulence_models.stf import STF2017TurbulenceModel
from py_wake.examples.data.iea34_130rwt import IEA34_130_1WT_Surrogate
from py_wake.superposition_models import MaxSum
from py_wake.wind_turbines.power_ct_functions import SimpleYawModel
```
%% Cell type:markdown id: tags:
## Select site, turbines, wake model and additional models and set up PyWake objects
%% Cell type:code id: tags:
``` python
site = LillgrundSite()
windTurbines = IEA34_130_1WT_Surrogate()
wfm = IEA37SimpleBastankhahGaussian(site, windTurbines, turbulenceModel=STF2017TurbulenceModel(addedTurbulenceSuperpositionModel=MaxSum()))
```
%% Cell type:markdown id: tags:
## Choose flow cases
(this will determine the speed and accuracy of the simulation). In this example we will focus on only a few flow cases.
%% Cell type:code id: tags:
``` python
wsp = np.asarray([10, 15])
wdir = np.arange(0,360,45)
```
%% Cell type:markdown id: tags:
## Constrain loads
In this example we will calculate nominal loads and use this as a basis for the load constraint.
%% Cell type:code id: tags:
``` python
x, y = site.initial_position.T
#keeping only every second turbine as lillegrund turbines are approx. half the size of the iea 3.4MW
x = x[::2]
y = y[::2]
x_init = x
y_init = y
n_wt = x.size
i = n_wt
k = wsp.size
l = wdir.size
load_fact = 1.002
simulationResult = wfm(x,y,wd=wdir, ws=wsp)
nom_loads = simulationResult.loads('OneWT')['LDEL'].values
max_loads = nom_loads * load_fact
s = nom_loads.shape[0]
load_signals = ['del_blade_flap', 'del_blade_edge', 'del_tower_bottom_fa',
'del_tower_bottom_ss', 'del_tower_top_torsion']
```
%% Cell type:markdown id: tags:
## Configure the optimization
this includes e.g. selection of maximum number of iterations, convergence tolerance, optimizer algorithm and design variable boundaries
%% Cell type:code id: tags:
``` python
maxiter = 5
tol = 1e-8
driver = EasyScipyOptimizeDriver(optimizer='SLSQP', maxiter=maxiter, tol=tol)
ec = 1e-2
step = 1e-4
min_spacing = 260
xi, xa = x_init.min()-min_spacing, x_init.max()+min_spacing
yi, ya = y_init.min()-min_spacing, y_init.max()+min_spacing
boundary = np.asarray([[xi, ya], [xa, ya], [xa, yi], [xi, yi]])
```
%% Cell type:markdown id: tags:
## Setup cost function
%% Cell type:code id: tags:
``` python
def aep_load_func(x, y):
simres = wfm(x, y, wd=wdir, ws=wsp)
aep = simres.aep().sum()
loads = simres.loads('OneWT')['LDEL'].values
return aep, loads
```
%% Cell type:markdown id: tags:
## Setup gradient function
In this example we will rely on the automatic finite difference, so no need to specify gradients
%% Cell type:markdown id: tags:
## Wrap your pure python cost and gradient functions in a topfarm component
%% Cell type:code id: tags:
``` python
cost_comp = AEPMaxLoadCostModelComponent(input_keys=[('x', x_init),('y', y_init)],
n_wt = n_wt,
aep_load_function = aep_load_func,
max_loads = max_loads,
objective=True,
step={'x': step, 'y': step},
output_keys=[('AEP', 0), ('loads', np.zeros((s, i)))]
)
```
%% Cell type:markdown id: tags:
## Set up the TopFarm problem
%% Cell type:code id: tags:
``` python
problem = TopFarmProblem(design_vars={'x': x_init, 'y': y_init},
constraints=[XYBoundaryConstraint(boundary),
SpacingConstraint(min_spacing)],
cost_comp=cost_comp,
driver=driver,
plot_comp=NoPlot(),
expected_cost=ec)
```
%% Cell type:markdown id: tags:
## Optimize
%% Cell type:code id: tags:
``` python
cost, state, recorder = problem.optimize()
```
%% Output
Iteration limit exceeded (Exit mode 9)
Current function value: -36978.307022695444
Iterations: 6
Function evaluations: 6
Gradient evaluations: 6
Optimization FAILED.
Iteration limit exceeded
-----------------------------------
%% Cell type:markdown id: tags:
## Plot results
Try to run the commands below to watch the resulting wake map for different flow cases
%% Cell type:code id: tags:
``` python
# import matplotlib.pyplot as plt
# plt.plot(recorder['AEP'])
```
%% Output
[<matplotlib.lines.Line2D at 0x7f330e4166d0>]
%% Cell type:code id: tags:
``` python
# n_i = recorder['counter'].size
# loads = recorder['loads'].reshape((n_i,s,n_wt))
# wt = 0
# for n, ls in enumerate(load_signals):
# plt.plot(loads[:,n,wt])
# plt.title(ls+f' turbine {wt}')
# plt.plot([0, n_i], [max_loads[n, wt], max_loads[n, wt]])
# plt.show()
```
%% Output
%% Cell type:code id: tags:
``` python
```
......
_notebooks/elements/roads_and_cables.ipynb
View file @ 4a7886d4
This source diff could not be displayed because it is too large. You can view the blob instead.
_notebooks/elements/wake_steering_and_loads.ipynb
View file @ 4a7886d4
%% Cell type:markdown id: tags:
# Load Constrained Wake Steering Optimization
%% Cell type:markdown id: tags:
## Install TopFarm and PyWake
%% Cell type:code id: tags:
``` python
%%capture
try:
import py_wake
except:
!pip install git+https://gitlab.windenergy.dtu.dk/TOPFARM/PyWake.git@setup_for_surrogate
!pip install git+https://gitlab.windenergy.dtu.dk/TOPFARM/PyWake.git
try:
import topfarm
except:
!pip install git+https://gitlab.windenergy.dtu.dk/TOPFARM/TopFarm2.git@opt_course
!pip install topfarm
```
%% Cell type:markdown id: tags:
## Import section
%% Cell type:code id: tags:
``` python
import numpy as np
from numpy import newaxis as na
import time
from topfarm.cost_models.cost_model_wrappers import AEPMaxLoadCostModelComponent
from topfarm.easy_drivers import EasyScipyOptimizeDriver
from topfarm import TopFarmProblem
from topfarm.plotting import NoPlot
from py_wake.examples.data.lillgrund import LillgrundSite
from py_wake.deficit_models.gaussian import IEA37SimpleBastankhahGaussian, NiayifarGaussian
from py_wake.turbulence_models.stf import STF2017TurbulenceModel
from py_wake.examples.data.iea34_130rwt import IEA34_130_1WT_Surrogate
from py_wake.deflection_models.jimenez import JimenezWakeDeflection
from py_wake.superposition_models import MaxSum
from py_wake.wind_turbines.power_ct_functions import SimpleYawModel
```
%% Cell type:markdown id: tags:
## Select site, turbines, wake model and additional models and set up PyWake objects
%% Cell type:code id: tags:
``` python
site = LillgrundSite()
windTurbines = IEA34_130_1WT_Surrogate()
wfm = IEA37SimpleBastankhahGaussian(site, windTurbines,deflectionModel=JimenezWakeDeflection(), turbulenceModel=STF2017TurbulenceModel(addedTurbulenceSuperpositionModel=MaxSum()))
```
%% Cell type:markdown id: tags:
## Choose flow cases
(this will determine the speed and accuracy of the simulation). In this example we will focus on only a few flow cases.
%% Cell type:code id: tags:
``` python
wsp = np.asarray([10, 15])
wdir = np.asarray([90])
```
%% Cell type:markdown id: tags:
## Constrain loads
In this example we will calculate nominal loads (without yaw) and use this as a basis for the load constraint.
%% Cell type:code id: tags:
``` python
x, y = site.initial_position.T
#keeping only every second turbine as lillegrund turbines are approx. half the size of the iea 3.4MW
x = x[::2]
y = y[::2]
n_wt = x.size
i = n_wt
k = wsp.size
l = wdir.size
yaw_zero = np.zeros((i, l, k))
load_fact = 1.02
simulationResult = wfm(x,y,wd=wdir, ws=wsp, yaw=yaw_zero)
nom_loads = simulationResult.loads('OneWT')['LDEL'].values
max_loads = nom_loads * load_fact
```
%% Output
WARNING:tensorflow:5 out of the last 32 calls to <function Model.make_predict_function.<locals>.predict_function at 0x000001E42258FC10> triggered tf.function retracing. Tracing is expensive and the excessive number of tracings could be due to (1) creating @tf.function repeatedly in a loop, (2) passing tensors with different shapes, (3) passing Python objects instead of tensors. For (1), please define your @tf.function outside of the loop. For (2), @tf.function has experimental_relax_shapes=True option that relaxes argument shapes that can avoid unnecessary retracing. For (3), please refer to https://www.tensorflow.org/guide/function#controlling_retracing and https://www.tensorflow.org/api_docs/python/tf/function for more details.
WARNING:tensorflow:6 out of the last 34 calls to <function Model.make_predict_function.<locals>.predict_function at 0x000001E420D444C0> triggered tf.function retracing. Tracing is expensive and the excessive number of tracings could be due to (1) creating @tf.function repeatedly in a loop, (2) passing tensors with different shapes, (3) passing Python objects instead of tensors. For (1), please define your @tf.function outside of the loop. For (2), @tf.function has experimental_relax_shapes=True option that relaxes argument shapes that can avoid unnecessary retracing. For (3), please refer to https://www.tensorflow.org/guide/function#controlling_retracing and https://www.tensorflow.org/api_docs/python/tf/function for more details.
%% Cell type:markdown id: tags:
## Configure the optimization
this includes e.g. selection of maximum number of iterations, convergence tolerance, optimizer algorithm and design variable boundaries
%% Cell type:code id: tags:
``` python
maxiter = 5
driver = EasyScipyOptimizeDriver(optimizer='SLSQP', maxiter=maxiter)
yaw_min, yaw_max = - 40, 40
yaw_init = np.zeros((i, l, k))
tol = 1e-8
ec = 1e-4
step = 1e-2
```
%% Cell type:markdown id: tags:
## Setup cost function
%% Cell type:code id: tags:
``` python
def aep_load_func(yaw_ilk):
simres = wfm(x, y, wd=wdir, ws=wsp, yaw=yaw_ilk)
aep = simres.aep().sum()
loads = simres.loads('OneWT')['LDEL'].values
return aep, loads
```
%% Cell type:markdown id: tags:
## Setup gradient function
For some problems it is sufficient to rely on the automatic finite difference calculated by OpenMDAO or you can specify the explicit gradients from your model. In this case we don't have explicit gradients but the automatic finite difference is also inefficient so we do a manual population of the Jacobian
%% Cell type:code id: tags:
``` python
s = nom_loads.shape[0]
P_ilk = np.broadcast_to(simulationResult.P.values[na], (i, l, k))
lifetime = float(60 * 60 * 24 * 365 * 20)
f1zh = 10.0 ** 7.0
lifetime_on_f1zh = lifetime / f1zh
indices = np.arange(i * l * k).reshape((i, l, k))
def aep_load_gradient(yaw_ilk):
simres0 = wfm(x, y, wd=wdir, ws=wsp, yaw=yaw_ilk)
aep0 = simres0.aep()
DEL0 = simulationResult.loads('OneWT')['DEL'].values
LDEL0 = simulationResult.loads('OneWT')['LDEL'].values
d_aep_d_yaw = np.zeros(i*l*k)
d_load_d_yaw = np.zeros((s * i, i * l * k))
for n in range(n_wt):
yaw_step = yaw_ilk.copy()
yaw_step = yaw_step.reshape(i, l, k)
yaw_step[n, :, :] += step
simres_fd = wfm(x, y, wd=wdir, ws=wsp, yaw=yaw_step)
aep_fd = simres_fd.aep()
d_aep_d_yaw[n * l * k : (n + 1) * l * k] = (((aep_fd.values - aep0.values) / step).sum((0))).ravel()
DEL_fd = simres_fd.loads('OneWT')['DEL'].values
for _ls in range(s):
m = simulationResult.loads('OneWT').m.values[_ls]
for _wd in range(l):
for _ws in range(k):
DEL_fd_fc = DEL0.copy()
DEL_fd_fc[:, :, _wd, _ws] = DEL_fd[:, :, _wd, _ws]
DEL_fd_fc = DEL_fd_fc[_ls, :, :, :]
f = DEL_fd_fc.mean()
LDEL_fd = (((P_ilk * (DEL_fd_fc/f) ** m).sum((1, 2)) * lifetime_on_f1zh) ** (1/m))*f
d_load_d_yaw[n_wt * _ls : n_wt * (_ls + 1), indices[n, _wd, _ws]] = (LDEL_fd - LDEL0[_ls]) / step
return d_aep_d_yaw, d_load_d_yaw
```
%% Cell type:markdown id: tags:
## Wrap your pure python cost and gradient functions in a topfarm component
%% Cell type:code id: tags:
``` python
cost_comp = AEPMaxLoadCostModelComponent(input_keys=[('yaw_ilk', np.zeros((i, l, k)))],
n_wt = n_wt,
aep_load_function = aep_load_func,
aep_load_gradient = aep_load_gradient,
max_loads = max_loads,
objective=True,
income_model=True,
output_keys=[('AEP', 0), ('loads', np.zeros((s, i)))]
)
```
%% Cell type:markdown id: tags:
## Set up the TopFarm problem
%% Cell type:code id: tags:
``` python
problem = TopFarmProblem(design_vars={'yaw_ilk': (yaw_init, yaw_min, yaw_max)},
cost_comp=cost_comp,
driver=EasyScipyOptimizeDriver(optimizer='SLSQP', maxiter=maxiter, tol=tol),
plot_comp=NoPlot(),
expected_cost=ec)
```
%% Cell type:markdown id: tags:
## Optimize
%% Cell type:code id: tags:
``` python
cost, state, recorder = problem.optimize()
```
%% Output
Iteration limit reached (Exit mode 9)
Current function value: [-738670.14118039]
Iterations: 5
Function evaluations: 8
Gradient evaluations: 5
Optimization FAILED.
Iteration limit reached
-----------------------------------
%% Cell type:markdown id: tags:
## Plot results
Try to run the commands below to watch the resulting wake map for different flow cases
%% Cell type:code id: tags:
``` python
# simulationResult = wfm(x,y,wd=wdir[0], ws=wsp[0], yaw=state['yaw_ilk'][:,0,0])
# simulationResult.flow_map().plot_wake_map()
```
%% Output
<matplotlib.contour.QuadContourSet at 0x1e423f219d0>
%% Cell type:code id: tags:
``` python
```
......
_notebooks/make_notebooks.py
View file @ 4a7886d4
......@@ -93,7 +93,7 @@ if __name__ == '__main__':
notebooks = ['constraints', 'cost_models', 'drivers', 'loads', 'problems',
'roads_and_cables', 'wake_steering_and_loads', 'layout_and_loads']
notebooks.remove('roads_and_cables')
# notebooks.remove('roads_and_cables')
notebooks.remove('loads')
check_notebooks(notebooks)
make_doc_notebooks(notebooks)
......
docs/notebooks/layout_and_loads.ipynb
View file @ 4a7886d4
%% Cell type:markdown id: tags:
# Load Constrained Layout Optimization
%% Cell type:markdown id: tags:
[Try this yourself](https://colab.research.google.com/github/DTUWindEnergy/TopFarm2/blob/master/docs/notebooks/layout_and_loads.ipynb) (requires google account)
%% Cell type:markdown id: tags:
## Install TopFarm and PyWake
%% Cell type:code id: tags:
``` python
%%capture
try:
import py_wake
except:
!pip install git+https://gitlab.windenergy.dtu.dk/TOPFARM/PyWake.git@setup_for_surrogate
!pip install git+https://gitlab.windenergy.dtu.dk/TOPFARM/PyWake.git
try:
import topfarm
except:
!pip install git+https://gitlab.windenergy.dtu.dk/TOPFARM/TopFarm2.git@opt_course
!pip install topfarm
```
%% Cell type:markdown id: tags:
## Import section
%% Cell type:code id: tags:
``` python
import numpy as np
from numpy import newaxis as na
import time
from topfarm.cost_models.cost_model_wrappers import AEPMaxLoadCostModelComponent
from topfarm.easy_drivers import EasyScipyOptimizeDriver
from topfarm import TopFarmProblem
from topfarm.plotting import NoPlot
from topfarm.constraint_components.boundary import XYBoundaryConstraint
from topfarm.constraint_components.spacing import SpacingConstraint
from py_wake.examples.data.lillgrund import LillgrundSite
from py_wake.deficit_models.gaussian import IEA37SimpleBastankhahGaussian, NiayifarGaussian
from py_wake.turbulence_models.stf import STF2017TurbulenceModel
from py_wake.examples.data.iea34_130rwt import IEA34_130_1WT_Surrogate
from py_wake.superposition_models import MaxSum
from py_wake.wind_turbines.power_ct_functions import SimpleYawModel
```
%% Cell type:markdown id: tags:
## Select site, turbines, wake model and additional models and set up PyWake objects
%% Cell type:code id: tags:
``` python
site = LillgrundSite()
windTurbines = IEA34_130_1WT_Surrogate()
wfm = IEA37SimpleBastankhahGaussian(site, windTurbines, turbulenceModel=STF2017TurbulenceModel(addedTurbulenceSuperpositionModel=MaxSum()))
```
%% Cell type:markdown id: tags:
## Choose flow cases
(this will determine the speed and accuracy of the simulation). In this example we will focus on only a few flow cases.
%% Cell type:code id: tags:
``` python
wsp = np.asarray([10, 15])
wdir = np.arange(0,360,45)
```
%% Cell type:markdown id: tags:
## Constrain loads
In this example we will calculate nominal loads and use this as a basis for the load constraint.
%% Cell type:code id: tags:
``` python
x, y = site.initial_position.T
#keeping only every second turbine as lillegrund turbines are approx. half the size of the iea 3.4MW
x = x[::2]
y = y[::2]
x_init = x
y_init = y
n_wt = x.size
i = n_wt
k = wsp.size
l = wdir.size
load_fact = 1.002
simulationResult = wfm(x,y,wd=wdir, ws=wsp)
nom_loads = simulationResult.loads('OneWT')['LDEL'].values
max_loads = nom_loads * load_fact
s = nom_loads.shape[0]
load_signals = ['del_blade_flap', 'del_blade_edge', 'del_tower_bottom_fa',
'del_tower_bottom_ss', 'del_tower_top_torsion']
```
%% Cell type:markdown id: tags:
## Configure the optimization
this includes e.g. selection of maximum number of iterations, convergence tolerance, optimizer algorithm and design variable boundaries
%% Cell type:code id: tags:
``` python
maxiter = 5
tol = 1e-8
driver = EasyScipyOptimizeDriver(optimizer='SLSQP', maxiter=maxiter, tol=tol)
ec = 1e-2
step = 1e-4
min_spacing = 260
xi, xa = x_init.min()-min_spacing, x_init.max()+min_spacing
yi, ya = y_init.min()-min_spacing, y_init.max()+min_spacing
boundary = np.asarray([[xi, ya], [xa, ya], [xa, yi], [xi, yi]])
```
%% Cell type:markdown id: tags:
## Setup cost function
%% Cell type:code id: tags:
``` python
def aep_load_func(x, y):
simres = wfm(x, y, wd=wdir, ws=wsp)
aep = simres.aep().sum()
loads = simres.loads('OneWT')['LDEL'].values
return aep, loads
```
%% Cell type:markdown id: tags:
## Setup gradient function
In this example we will rely on the automatic finite difference, so no need to specify gradients
%% Cell type:markdown id: tags:
## Wrap your pure python cost and gradient functions in a topfarm component
%% Cell type:code id: tags:
``` python
cost_comp = AEPMaxLoadCostModelComponent(input_keys=[('x', x_init),('y', y_init)],
n_wt = n_wt,
aep_load_function = aep_load_func,
max_loads = max_loads,
objective=True,
step={'x': step, 'y': step},
output_keys=[('AEP', 0), ('loads', np.zeros((s, i)))]
)
```
%% Cell type:markdown id: tags:
## Set up the TopFarm problem
%% Cell type:code id: tags:
``` python
problem = TopFarmProblem(design_vars={'x': x_init, 'y': y_init},
constraints=[XYBoundaryConstraint(boundary),
SpacingConstraint(min_spacing)],
cost_comp=cost_comp,
driver=driver,
plot_comp=NoPlot(),
expected_cost=ec)
```
%% Cell type:markdown id: tags:
## Optimize
%% Cell type:code id: tags:
``` python
cost, state, recorder = problem.optimize()
```
%% Output
Iteration limit exceeded (Exit mode 9)
Current function value: -36978.307022695444
Iterations: 6
Function evaluations: 6
Gradient evaluations: 6
Optimization FAILED.
Iteration limit exceeded
-----------------------------------
%% Cell type:markdown id: tags:
## Plot results
Try to run the commands below to watch the resulting wake map for different flow cases
%% Cell type:code id: tags:
``` python
# import matplotlib.pyplot as plt
# plt.plot(recorder['AEP'])
```
%% Output
[<matplotlib.lines.Line2D at 0x7f330e4166d0>]
%% Cell type:code id: tags:
``` python
# n_i = recorder['counter'].size
# loads = recorder['loads'].reshape((n_i,s,n_wt))
# wt = 0
# for n, ls in enumerate(load_signals):
# plt.plot(loads[:,n,wt])
# plt.title(ls+f' turbine {wt}')
# plt.plot([0, n_i], [max_loads[n, wt], max_loads[n, wt]])
# plt.show()
```
%% Output
%% Cell type:code id: tags:
``` python
```
......
docs/notebooks/roads_and_cables.ipynb
View file @ 4a7886d4
This diff is collapsed.
docs/notebooks/wake_steering_and_loads.ipynb
View file @ 4a7886d4
%% Cell type:markdown id: tags:
# Load Constrained Wake Steering Optimization
%% Cell type:markdown id: tags:
[Try this yourself](https://colab.research.google.com/github/DTUWindEnergy/TopFarm2/blob/master/docs/notebooks/wake_steering_and_loads.ipynb) (requires google account)
%% Cell type:markdown id: tags:
## Install TopFarm and PyWake
%% Cell type:code id: tags:
``` python
%%capture
try:
import py_wake
except:
!pip install git+https://gitlab.windenergy.dtu.dk/TOPFARM/PyWake.git@setup_for_surrogate
!pip install git+https://gitlab.windenergy.dtu.dk/TOPFARM/PyWake.git
try:
import topfarm
except:
!pip install git+https://gitlab.windenergy.dtu.dk/TOPFARM/TopFarm2.git@opt_course
!pip install topfarm
```
%% Cell type:markdown id: tags:
## Import section
%% Cell type:code id: tags:
``` python
import numpy as np
from numpy import newaxis as na
import time
from topfarm.cost_models.cost_model_wrappers import AEPMaxLoadCostModelComponent
from topfarm.easy_drivers import EasyScipyOptimizeDriver
from topfarm import TopFarmProblem
from topfarm.plotting import NoPlot
from py_wake.examples.data.lillgrund import LillgrundSite
from py_wake.deficit_models.gaussian import IEA37SimpleBastankhahGaussian, NiayifarGaussian
from py_wake.turbulence_models.stf import STF2017TurbulenceModel
from py_wake.examples.data.iea34_130rwt import IEA34_130_1WT_Surrogate
from py_wake.deflection_models.jimenez import JimenezWakeDeflection
from py_wake.superposition_models import MaxSum
from py_wake.wind_turbines.power_ct_functions import SimpleYawModel
```
%% Cell type:markdown id: tags:
## Select site, turbines, wake model and additional models and set up PyWake objects
%% Cell type:code id: tags:
``` python
site = LillgrundSite()
windTurbines = IEA34_130_1WT_Surrogate()
wfm = IEA37SimpleBastankhahGaussian(site, windTurbines,deflectionModel=JimenezWakeDeflection(), turbulenceModel=STF2017TurbulenceModel(addedTurbulenceSuperpositionModel=MaxSum()))
```
%% Cell type:markdown id: tags:
## Choose flow cases
(this will determine the speed and accuracy of the simulation). In this example we will focus on only a few flow cases.
%% Cell type:code id: tags:
``` python
wsp = np.asarray([10, 15])
wdir = np.asarray([90])
```
%% Cell type:markdown id: tags:
## Constrain loads
In this example we will calculate nominal loads (without yaw) and use this as a basis for the load constraint.
%% Cell type:code id: tags:
``` python
x, y = site.initial_position.T
#keeping only every second turbine as lillegrund turbines are approx. half the size of the iea 3.4MW
x = x[::2]
y = y[::2]
n_wt = x.size
i = n_wt
k = wsp.size
l = wdir.size
yaw_zero = np.zeros((i, l, k))
load_fact = 1.02
simulationResult = wfm(x,y,wd=wdir, ws=wsp, yaw=yaw_zero)
nom_loads = simulationResult.loads('OneWT')['LDEL'].values
max_loads = nom_loads * load_fact
```
%% Output
WARNING:tensorflow:5 out of the last 32 calls to <function Model.make_predict_function.<locals>.predict_function at 0x000001E42258FC10> triggered tf.function retracing. Tracing is expensive and the excessive number of tracings could be due to (1) creating @tf.function repeatedly in a loop, (2) passing tensors with different shapes, (3) passing Python objects instead of tensors. For (1), please define your @tf.function outside of the loop. For (2), @tf.function has experimental_relax_shapes=True option that relaxes argument shapes that can avoid unnecessary retracing. For (3), please refer to https://www.tensorflow.org/guide/function#controlling_retracing and https://www.tensorflow.org/api_docs/python/tf/function for more details.
WARNING:tensorflow:6 out of the last 34 calls to <function Model.make_predict_function.<locals>.predict_function at 0x000001E420D444C0> triggered tf.function retracing. Tracing is expensive and the excessive number of tracings could be due to (1) creating @tf.function repeatedly in a loop, (2) passing tensors with different shapes, (3) passing Python objects instead of tensors. For (1), please define your @tf.function outside of the loop. For (2), @tf.function has experimental_relax_shapes=True option that relaxes argument shapes that can avoid unnecessary retracing. For (3), please refer to https://www.tensorflow.org/guide/function#controlling_retracing and https://www.tensorflow.org/api_docs/python/tf/function for more details.
%% Cell type:markdown id: tags:
## Configure the optimization
this includes e.g. selection of maximum number of iterations, convergence tolerance, optimizer algorithm and design variable boundaries
%% Cell type:code id: tags:
``` python
maxiter = 5
driver = EasyScipyOptimizeDriver(optimizer='SLSQP', maxiter=maxiter)
yaw_min, yaw_max = - 40, 40
yaw_init = np.zeros((i, l, k))
tol = 1e-8
ec = 1e-4
step = 1e-2
```
%% Cell type:markdown id: tags:
## Setup cost function
%% Cell type:code id: tags:
``` python
def aep_load_func(yaw_ilk):
simres = wfm(x, y, wd=wdir, ws=wsp, yaw=yaw_ilk)
aep = simres.aep().sum()
loads = simres.loads('OneWT')['LDEL'].values
return aep, loads
```
%% Cell type:markdown id: tags:
## Setup gradient function
For some problems it is sufficient to rely on the automatic finite difference calculated by OpenMDAO or you can specify the explicit gradients from your model. In this case we don't have explicit gradients but the automatic finite difference is also inefficient so we do a manual population of the Jacobian
%% Cell type:code id: tags:
``` python
s = nom_loads.shape[0]
P_ilk = np.broadcast_to(simulationResult.P.values[na], (i, l, k))
lifetime = float(60 * 60 * 24 * 365 * 20)
f1zh = 10.0 ** 7.0
lifetime_on_f1zh = lifetime / f1zh
indices = np.arange(i * l * k).reshape((i, l, k))
def aep_load_gradient(yaw_ilk):
simres0 = wfm(x, y, wd=wdir, ws=wsp, yaw=yaw_ilk)
aep0 = simres0.aep()
DEL0 = simulationResult.loads('OneWT')['DEL'].values
LDEL0 = simulationResult.loads('OneWT')['LDEL'].values
d_aep_d_yaw = np.zeros(i*l*k)
d_load_d_yaw = np.zeros((s * i, i * l * k))
for n in range(n_wt):
yaw_step = yaw_ilk.copy()
yaw_step = yaw_step.reshape(i, l, k)
yaw_step[n, :, :] += step
simres_fd = wfm(x, y, wd=wdir, ws=wsp, yaw=yaw_step)
aep_fd = simres_fd.aep()
d_aep_d_yaw[n * l * k : (n + 1) * l * k] = (((aep_fd.values - aep0.values) / step).sum((0))).ravel()
DEL_fd = simres_fd.loads('OneWT')['DEL'].values
for _ls in range(s):
m = simulationResult.loads('OneWT').m.values[_ls]
for _wd in range(l):
for _ws in range(k):
DEL_fd_fc = DEL0.copy()
DEL_fd_fc[:, :, _wd, _ws] = DEL_fd[:, :, _wd, _ws]
DEL_fd_fc = DEL_fd_fc[_ls, :, :, :]
f = DEL_fd_fc.mean()
LDEL_fd = (((P_ilk * (DEL_fd_fc/f) ** m).sum((1, 2)) * lifetime_on_f1zh) ** (1/m))*f
d_load_d_yaw[n_wt * _ls : n_wt * (_ls + 1), indices[n, _wd, _ws]] = (LDEL_fd - LDEL0[_ls]) / step
return d_aep_d_yaw, d_load_d_yaw
```
%% Cell type:markdown id: tags:
## Wrap your pure python cost and gradient functions in a topfarm component
%% Cell type:code id: tags:
``` python
cost_comp = AEPMaxLoadCostModelComponent(input_keys=[('yaw_ilk', np.zeros((i, l, k)))],
n_wt = n_wt,
aep_load_function = aep_load_func,
aep_load_gradient = aep_load_gradient,
max_loads = max_loads,
objective=True,
income_model=True,
output_keys=[('AEP', 0), ('loads', np.zeros((s, i)))]
)
```
%% Cell type:markdown id: tags:
## Set up the TopFarm problem
%% Cell type:code id: tags:
``` python
problem = TopFarmProblem(design_vars={'yaw_ilk': (yaw_init, yaw_min, yaw_max)},
cost_comp=cost_comp,
driver=EasyScipyOptimizeDriver(optimizer='SLSQP', maxiter=maxiter, tol=tol),
plot_comp=NoPlot(),
expected_cost=ec)
```
%% Cell type:markdown id: tags:
## Optimize
%% Cell type:code id: tags:
``` python
cost, state, recorder = problem.optimize()
```
%% Output
Iteration limit reached (Exit mode 9)
Current function value: [-738670.14118039]
Iterations: 5
Function evaluations: 8
Gradient evaluations: 5
Optimization FAILED.
Iteration limit reached
-----------------------------------
%% Cell type:markdown id: tags:
## Plot results
Try to run the commands below to watch the resulting wake map for different flow cases
%% Cell type:code id: tags:
``` python
# simulationResult = wfm(x,y,wd=wdir[0], ws=wsp[0], yaw=state['yaw_ilk'][:,0,0])
# simulationResult.flow_map().plot_wake_map()
```
%% Output
<matplotlib.contour.QuadContourSet at 0x1e423f219d0>
%% Cell type:code id: tags:
``` python
```
......
setup.py
View file @ 4a7886d4
......@@ -28,6 +28,7 @@ setup(name='topfarm',
install_requires=[
'matplotlib', # for plotting
'numpy', # for numerical calculations
'numpy-financial', # for irr calculations
'openmdao==2.6', # for optimization
'pytest', # for testing
'pytest-cov', # for calculating coverage
......
topfarm/cost_models/economic_models/dtu_wind_cm_main.py
View file @ 4a7886d4
......@@ -4,6 +4,7 @@ Implementation of Offshore Turbine Cost Model, ref. Soeren Oemann Lind, soli@dtu
Structure based on 2017 turbine_cost.py by Arvydas Berzonskis, Goldwind
"""
import numpy as np
import numpy_financial as npf
class economic_evaluation():
......@@ -45,7 +46,7 @@ class economic_evaluation():
self.calculate_cash_flow()
self.NPV = np.npv(self.discount_rate, self.CWF)
self.NPV = npf.npv(self.discount_rate, self.CWF)
return self.NPV
......@@ -74,7 +75,7 @@ class economic_evaluation():
self.calculate_cash_flow()
self.IRR = 100 * np.irr(self.CWF) # [%]
self.IRR = 100 * npf.irr(self.CWF) # [%]
return self.IRR
......@@ -89,13 +90,10 @@ class economic_evaluation():
self.CWF.append(-self.project_costs_sums["DEVEX"] - self.project_costs_sums["CAPEX"])
elif i == self.project_duration:
self.CWF.append(int(self.energy_price * sum(self.aep_vector) -
self.project_costs_sums["ABEX"] -
self.project_costs_sums["OPEX"]))
self.CWF.append(self.energy_price * sum(self.aep_vector) - self.project_costs_sums["ABEX"] - self.project_costs_sums["OPEX"])
else:
self.CWF.append(int(self.energy_price * sum(self.aep_vector) -
self.project_costs_sums["OPEX"]))
self.CWF.append(self.energy_price * sum(self.aep_vector) - self.project_costs_sums["OPEX"])
def calculate_expenditures(self, rated_rpm, rotor_diameter, rated_power,
hub_height, aep_array, water_depth, cabling_cost=None):
......
topfarm/cost_models/electrical/simple_msp.py
View file @ 4a7886d4
import topfarm
import openmdao.api as om
from topfarm.cost_models.utils.spanning_tree import spanning_tree
from topfarm.cost_models.utils.spanning_tree import mst
from topfarm.cost_models.cost_model_wrappers import CostModelComponent
from topfarm.plotting import XYPlotComp, NoPlot
from IPython.display import display
import numpy as np
import matplotlib.pylab as plt
......@@ -30,14 +31,14 @@ class ElNetLength(ElNetComp):
def compute(self, inputs, outputs):
x, y = inputs[topfarm.x_key], inputs[topfarm.y_key]
elnet_layout = spanning_tree(x, y)
elnet_layout = mst(x, y)
outputs['elnet_length'] = sum(list(elnet_layout.values()))
class PlotElNet(ElNetComp):
def compute(self, inputs, outputs):
x, y = inputs[topfarm.x_key], inputs[topfarm.y_key]
elnet_layout = spanning_tree(x, y)
elnet_layout = mst(x, y)
indices = np.array(list(elnet_layout.keys())).T
plt.plot(x[indices], y[indices], color='r')
plt.plot(x, y, 'o')
......@@ -50,7 +51,7 @@ class ElNetCost(CostModelComponent):
self.n_wt = n_wt
self.length_key = length_key
CostModelComponent.__init__(self, [(self.length_key, 0.0)], self.n_wt,
self.cost, self.grad, objective=False,
self.cost, self.grad, objective=False, input_units=['m'],
**kwargs)
def initialize(self):
......@@ -66,10 +67,30 @@ class ElNetCost(CostModelComponent):
class XYElPlotComp(XYPlotComp):
"""Plotting component for turbine locations"""
def plot_current_position(self, x, y):
elnet_layout = spanning_tree(x, y)
elnet_layout = mst(x, y)
indices = np.array(list(elnet_layout.keys())).T
plt.plot(x[indices], y[indices], color='r')
for c, x_, y_ in zip(self.colors, x, y):