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

added opt course examples, and added gradient calculation of more than one...

added opt course examples, and added gradient calculation of more than one output when using the wrapper
parent c7002cde
Pipeline #23639 passed with stages
in 13 minutes and 34 seconds
......@@ -52,6 +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
......
This diff is collapsed.
This diff is collapsed.
......@@ -66,7 +66,7 @@ try:
import topfarm
except ModuleNotFoundError:
!pip install topfarm"""
if not name == 'loads':
if not name in ['loads', 'wake_steering_and_loads', 'layout_and_loads']:
nb.insert_code_cell(2, code)
nb.save(dst_path + name + ".ipynb")
......@@ -91,7 +91,8 @@ def check_notebooks(notebooks=None):
if __name__ == '__main__':
notebooks = ['constraints', 'cost_models', 'drivers', 'loads', 'problems', 'roads_and_cables']
notebooks = ['constraints', 'cost_models', 'drivers', 'loads', 'problems',
'roads_and_cables', 'wake_steering_and_loads', 'layout_and_loads']
notebooks.remove('roads_and_cables')
notebooks.remove('loads')
check_notebooks(notebooks)
......
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
......@@ -17,5 +17,7 @@ notebooks is `here on GitLab <https://gitlab.windenergy.dtu.dk/TOPFARM/TopFarm2/
examples_nblinks/drivers_nb
examples_nblinks/problems_nb
examples_nblinks/cost_models_nb
examples_nblinks/loads_nb
examples_nblinks/roads_and_cables_nb
\ No newline at end of file
examples_nblinks/roads_and_cables_nb
examples_nblinks/wake_steering_and_loads_nb
examples_nblinks/layout_and_loads_nb
{
"path": "../../notebooks/layout_and_loads.ipynb"
}
\ No newline at end of file
{
"path": "../../notebooks/loads.ipynb"
}
\ No newline at end of file
{
"path": "../../notebooks/wake_steering_and_loads.ipynb"
}
\ No newline at end of file
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
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
site = LillgrundSite()
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
# # Wind turbines and wind farm model definition
windTurbines = IEA34_130_1WT_Surrogate()
wfm = IEA37SimpleBastankhahGaussian(site, windTurbines, turbulenceModel=STF2017TurbulenceModel(addedTurbulenceSuperpositionModel=MaxSum()))
load_signals = ['del_blade_flap', 'del_blade_edge', 'del_tower_bottom_fa',
'del_tower_bottom_ss', 'del_tower_top_torsion']
wsp = np.asarray([10, 15])
wdir = np.arange(0,360,45)
n_wt = x.size
i = n_wt
k = wsp.size
l = wdir.size
yaw_zero = np.zeros((i, l, k))
maxiter = 10
yaw_min, yaw_max = - 40, 40
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
step = 1e-4
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_func(x, y):
simres = wfm(x, y, wd=wdir, ws=wsp)
aep = simres.aep().sum()
loads = simres.loads('OneWT')['LDEL'].values
return aep, loads
tol = 1e-8
ec = 1e-1
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]])
cost_comp = AEPMaxLoadCostModelComponent(input_keys=[('x', x_init),('y', y_init)],
n_wt = n_wt,
aep_load_function = aep_load_func,
# aep_load_gradient = aep_load_gradient,
max_loads = max_loads,
objective=True,
step={'x': step, 'y': step},
output_keys=[('AEP', 0), ('loads', np.zeros((s, i)))]
)
problem = TopFarmProblem(design_vars={'x': x_init, 'y': y_init},
constraints=[XYBoundaryConstraint(boundary),
SpacingConstraint(min_spacing)],
# post_constraints=[(ls, val * load_fact) for ls, val in loads_nom.items()],
cost_comp=cost_comp,
driver=EasyScipyOptimizeDriver(optimizer='SLSQP', maxiter=maxiter, tol=tol),
plot_comp=NoPlot(),
expected_cost=ec)
tic = time.time()
if 1:
cost, state, recorder = problem.optimize()
toc = time.time()
print('Optimization took: {:.0f}s'.format(toc-tic))
if 0:
with open(f'./check_partials_{int(toc)}_{ec}_{step}.txt', 'w') as fid:
partials = problem.check_partials(out_stream=fid,
compact_print=True,
show_only_incorrect=True,
step=step)
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
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
site = LillgrundSite()
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]
# # Wind turbines and wind farm model definition
windTurbines = IEA34_130_1WT_Surrogate() #additional_models=[SimpleYawModel()]
wfm = IEA37SimpleBastankhahGaussian(site, windTurbines,deflectionModel=JimenezWakeDeflection(), turbulenceModel=STF2017TurbulenceModel(addedTurbulenceSuperpositionModel=MaxSum()))
load_signals = ['del_blade_flap', 'del_blade_edge', 'del_tower_bottom_fa',
'del_tower_bottom_ss', 'del_tower_top_torsion']
wsp = np.asarray([10, 15])
wdir = np.asarray([90])
n_wt = x.size
i = n_wt
k = wsp.size
l = wdir.size
yaw_zero = np.zeros((i, l, k))
maxiter = 10
driver = EasyScipyOptimizeDriver(optimizer='SLSQP', maxiter=maxiter)
yaw_min, yaw_max = - 40, 40
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
step = 1e-2
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_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
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
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,
output_keys=[('AEP', 0), ('loads', np.zeros((s, i)))]
)
yaw_init = np.zeros((i, l, k))
yaw_30 = np.full_like(yaw_init, 30)
yaw_init_rand = np.random.rand(i, l, k)*80-40
tol = 1e-8
ec = 1e-4
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)
tic = time.time()
if 1:
cost, state, recorder = problem.optimize()
toc = time.time()
print('Optimization took: {:.0f}s'.format(toc-tic))
if 0:
with open(f'./check_partials_{int(toc)}.txt', 'w') as fid:
partials = problem.check_partials(out_stream=fid,
compact_print=True,
show_only_incorrect=False,
step=step)
pycodestyle --ignore=E501,W504 --exclude="*Colonel*" topfarm
PAUSE
\ No newline at end of file
python -m pytest --ignore=topfarm/cost_models/fuga/Colonel
PAUSE
\ No newline at end of file
......@@ -37,5 +37,6 @@ setup(name='topfarm',
'sphinx_rtd_theme', # docs theme
'scikit-learn', # load surrogate
'mock', # replace variables during tests
'tensorflow', # loads examples with surrogates
],
zip_safe=True)
......@@ -3,14 +3,15 @@ import numpy as np
import time
from _collections import defaultdict
from topfarm.constraint_components.post_constraint import PostConstraint
import warnings
class CostModelComponent(ExplicitComponent):
"""Wrapper for pure-Python cost functions"""
def __init__(self, input_keys, n_wt, cost_function, cost_gradient_function=None,
output_key="Cost", output_unit="", additional_input=[], additional_output=[], max_eval=None,
objective=True, income_model=False, output_val=0.0, input_units=[], **kwargs):
output_keys=["Cost"], output_unit="", additional_input=[], additional_output=[], max_eval=None,
objective=True, maximize=False, output_vals=[0.0], input_units=[], step={}, **kwargs):
"""Initialize wrapper for pure-Python cost function
Parameters
......@@ -23,7 +24,7 @@ class CostModelComponent(ExplicitComponent):
Function to evaluate cost
cost_gradient_function : function handle, optional
Function to evaluate gradient of the cost function
output_key : str
output_key : list
Name of output key
output_unit : str
Units of output of cost function
......@@ -39,31 +40,58 @@ class CostModelComponent(ExplicitComponent):
Maximum number of function evaluations
objective : boolean
Must be True for standalone CostModelComponents and the final component in a TopFarmGroup
income_model : boolean
maximize : boolean
If True: objective is maximised during optimization\n
If False: Objective is minimized during optimization
output_val : float or array_like
Format of output
input_units : list of str
Units of the respective input_keys
step : dict of {str : float}
Finite difference step size for each input key, e.g. {input_key : step_size, input_key2 : step_size2}
"""
if 'income_model' in kwargs:
warnings.warn("""income_model is deprecated; use keyword maximize instead""",
DeprecationWarning, stacklevel=2)
maximize = kwargs['income_model']
kwargs.pop('income_model')
if 'output_key' in kwargs:
warnings.warn("""output_key is deprecated; use keyword output_keys instead""",
DeprecationWarning, stacklevel=2)
self.output_key = kwargs['output_key']
# if self.output_key not in output_keys:
output_keys = self.output_key
kwargs.pop('output_key')
else:
self.output_key = output_keys[0]
if isinstance(self.output_key, tuple):
self.output_key = self.output_key[0]
if 'output_val' in kwargs:
warnings.warn("""output_val is deprecated; use keyword output_vals instead""",
DeprecationWarning, stacklevel=2)
output_vals = [kwargs['output_val']]
kwargs.pop('output_val')
super().__init__(**kwargs)
assert isinstance(n_wt, int), n_wt
self.input_keys = list(input_keys)
self.cost_function = cost_function
self.cost_gradient_function = cost_gradient_function
self.n_wt = n_wt
self.output_key = output_key
if not isinstance(output_keys, list):
output_keys = [output_keys]
self.output_keys = output_keys
self.out_keys_only = list([(o, o[0])[isinstance(o, tuple)] for o in self.output_keys])
self.output_unit = output_unit
self.additional_input = additional_input
self.additional_output = [((x, np.zeros(self.n_wt)), x)[isinstance(x, tuple)] for x in additional_output]
self.max_eval = max_eval or 1e100
self.objective = objective
if income_model:
if maximize:
self.cost_factor = -1.0
else:
self.cost_factor = 1.0
self.output_val = output_val
output_vals = output_vals[:len(output_keys)] + [output_vals[0]] * (len(output_keys) - len(output_vals)) # extend output_vals if it is not same length as output_keys
self.output_vals = output_vals
n_input = len(self.input_keys) + len(additional_input)
self.input_units = (input_units + [None] * n_input)[:n_input]
......@@ -72,7 +100,7 @@ class CostModelComponent(ExplicitComponent):
self.func_time_sum = 0
self.n_grad_eval = 0
self.grad_time_sum = 0
self.step = {}
self.step = step
def setup(self):
for i, u in zip(self.input_keys + self.additional_input, self.input_units):
......@@ -84,30 +112,39 @@ class CostModelComponent(ExplicitComponent):
if self.objective:
self.add_output('cost', val=0.0)
self.add_output('cost_comp_eval', val=0.0)
self.add_output(self.output_key, val=self.output_val)
for o, v in zip(self.output_keys, self.output_vals):
if isinstance(o, tuple) and len(o) == 2:
self.add_output(o[0], val=o[1])
else:
self.add_output(o, val=v)
for key, val in self.additional_output:
self.add_output(key, val=val)
input_keys = list([(i, i[0])[isinstance(i, tuple)] for i in self.input_keys])
self.inp_keys = input_keys + list([(i, i[0])[isinstance(i, tuple)] for i in self.additional_input])
self.input_keys = input_keys
if self.objective:
if self.cost_gradient_function:
self.declare_partials('cost', input_keys)
else:
# Finite difference all partials.
if self.step == {}:
self.declare_partials('cost', input_keys, method='fd')
else:
for i in input_keys:
self.declare_partials('cost', i, step=self.step.get('cost', {}).get(i, None), method='fd')
self.declare_partials(self.output_key, input_keys)
if self.cost_gradient_function:
if self.objective:
self.declare_partials('cost', input_keys, method='exact')
# else:
for o in self.out_keys_only:
self.declare_partials(o, input_keys, method='exact')
else:
if self.cost_gradient_function:
self.declare_partials(self.output_key, input_keys)
if self.step == {}:
if self.objective:
self.declare_partials('cost', input_keys, method='fd')
# else:
for o in self.out_keys_only:
self.declare_partials(o, input_keys, method='fd')
else:
# Finite difference all partials.
self.declare_partials(self.output_key, input_keys, method='fd')
for i in input_keys:
if self.objective:
self.declare_partials('cost', i, step=self.step[i], method='fd')
# else:
for o in self.out_keys_only:
# self.declare_partials(o, input_keys, method='fd')
self.declare_partials(o, i, step=self.step[i], method='fd')
@property
def counter(self):
......@@ -133,10 +170,13 @@ class CostModelComponent(ExplicitComponent):
outputs[k] = v
else:
c = self.cost_function(**{x: inputs[x] for x in self.inp_keys})
if not isinstance(c, list):
c = [c]
if self.objective:
outputs['cost'] = c * self.cost_factor
outputs['cost'] = c[0] * self.cost_factor
outputs['cost_comp_eval'] = self.counter
outputs[self.output_key] = c
for o, _c in zip(self.out_keys_only, c):
outputs[o] = _c
self.func_time_sum += time.time() - t
self.n_func_eval += 1
......@@ -149,37 +189,52 @@ class CostModelComponent(ExplicitComponent):
for k, dCostdk in zip(self.input_keys,
self.cost_gradient_function(**{x: inputs[x] for x in self.inp_keys})):
if dCostdk is not None:
if not isinstance(dCostdk, list):
dCostdk = [dCostdk]
if self.objective:
J['cost', k] = dCostdk * self.cost_factor
J[self.output_key, k] = dCostdk
J['cost', k] = dCostdk[0] * self.cost_factor
for o, _d in zip(self.out_keys_only, dCostdk):
J[o, k] = _d
self.grad_time_sum += time.time() - t
self.n_grad_eval += 1
class AEPCostModelComponent(CostModelComponent):
def __init__(self, input_keys, n_wt, cost_function, cost_gradient_function=None,
output_unit="", additional_input=[], additional_output=[], max_eval=None, **kwargs):
output_unit="", additional_input=[], additional_output=[], max_eval=None,
output_key="AEP", **kwargs):
CostModelComponent.__init__(self, input_keys, n_wt, cost_function,
cost_gradient_function=cost_gradient_function,
output_key="AEP", output_unit=output_unit,
output_key=output_key, output_unit=output_unit,
additional_input=additional_input, additional_output=additional_output,
max_eval=max_eval, income_model=True, **kwargs)
max_eval=max_eval, maximize=True, **kwargs)
class AEPMaxLoadCostModelComponent(AEPCostModelComponent, PostConstraint):
def __init__(self, input_keys, n_wt, aep_load_function, max_loads, **kwargs):
class AEPMaxLoadCostModelComponent(CostModelComponent, PostConstraint):
def __init__(self, input_keys, n_wt, aep_load_function, max_loads,
aep_load_gradient=None, output_keys=["AEP", 'loads'], step={},
maximize=True, **kwargs):
self.max_loads = max_loads
def cost_function(**kwargs):
aep, load, additional_output = (aep_load_function(**kwargs) + ({},))[:3]
additional_output['loads'] = load
return aep, additional_output
additional_output = kwargs.get('additional_output', []) + ['loads']
AEPCostModelComponent.__init__(self, input_keys, n_wt, cost_function=cost_function,
additional_output=additional_output, **kwargs)
aep, load = aep_load_function(**kwargs)
return [aep, load]
if aep_load_gradient:
def gradient_function(**kwargs):
d_aep, d_load = aep_load_gradient(**kwargs)
return [[d_aep, d_load]]
else:
gradient_function = None
# additional_output = kwargs.get('additional_output', []) + [('loads', max_loads)]
CostModelComponent.__init__(self, input_keys, n_wt, cost_function=cost_function,
cost_gradient_function=gradient_function,
output_keys=output_keys, step=step, maximize=maximize,
**kwargs)
PostConstraint.__init__(self, 'loads', upper=max_loads)
def setup(self):
AEPCostModelComponent.setup(self)
input_keys = list([(i, i[0])[isinstance(i, tuple)] for i in self.input_keys])
self.declare_partials('loads', input_keys, method='fd')
# def setup(self):
# AEPCostModelComponent.setup(self)
# input_keys = list([(i, i[0])[isinstance(i, tuple)] for i in self.input_keys])
# self.declare_partials('loads', input_keys, method='fd')
......@@ -119,6 +119,7 @@ def test_AEPMaxLoadCostModelComponent_as_penalty_multi_wt():
design_vars={'x': ([0, 1]), 'y': ([0, 0])},
cost_comp=AEPMaxLoadCostModelComponent(
input_keys='xy', n_wt=2,
output_keys=["AEP", ('loads', [3, 3])],
aep_load_function=lambda x, y: (-np.sin(np.hypot(x, y)).sum(), np.hypot(x, y)),
max_loads=[3, 3]),
constraints=[],
......
......@@ -15,7 +15,8 @@ def get_notebooks():
@pytest.mark.parametrize("notebook", get_notebooks())
def test_notebooks(notebook):
if os.path.basename(notebook.filename) in ['loads.ipynb', 'roads_and_cables.ipynb']:
if os.path.basename(notebook.filename) in ['loads.ipynb', 'roads_and_cables.ipynb',
'wake_steering_and_loads.ipynb', 'layout_and_loads.ipynb']:
pytest.xfail("Notebook, %s, has known issues" % notebook)
import matplotlib.pyplot as plt
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment