I am doing a project(using pywake library which is user-defined lib) and I have written the following code:
enter code here
import numpy as np
from py_wake.examples.data.hornsrev1 import V80
from py_wake.examples.data.hornsrev1 import Hornsrev1Site # We work with the Horns Rev 1 site, which comes already set up with PyWake.
from py_wake.examples.data.hornsrev1 import wt_x, wt_y # The existing layout coordinates are also available from PyWake
from py_wake import BastankhahGaussian
import function
site = Hornsrev1Site()
x, y = site.initial_position.T
windTurbines = V80()
wf_model = BastankhahGaussian(site, windTurbines)
c=np.random.randint(0,2,size=len(x))
# c=np.random.randint(0,2,size=len(x)) # Switched Off/On WT
print(c)
x_new,y_new=function.funC(x, y, c)
print(wf_model)
# run wind farm simulation
sim_res = wf_model(x_new, y_new, # wind turbine positions
h=None, # wind turbine heights (defaults to the heights defined in windTurbines)
type=0, # Wind turbine types
wd=None, # Wind direction (defaults to site.default_wd (0,1,...,360 if not overriden))
ws=None, # Wind speed (defaults to site.default_ws (3,4,...,25m/s if not overriden))
)
print(sim_res.aep().sum())
just to make a bit clear the code, we input some data about a wind farm and then do the simulation and get the power output of wind turbine(print(sim_res.aep().sum()) ).
now I have defined a new variable(c) in which we have two values 0 and 1. if c=1 we say that that wind turbine is one otherwise it is off so power production will be decreased.
now by using the defined scripts I want to do an optimization in Scipy by using Penalty Function: I want to maximize the power production by changing the c values. I mean we want to switch off/on different wind turbines and see the output power production. I know that the optimized value is when all the parameters in c will be one but there are some limits so I need to use the penalty function. so would you please help me to how optimize the power production by using c and penalty?
"""
References:
https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html
https://github.com/DTUWindEnergy/PyWake
"""
import time
from py_wake.examples.data.hornsrev1 import V80
from py_wake.examples.data.hornsrev1 import Hornsrev1Site # We work with the Horns Rev 1 site, which comes already set up with PyWake.
from py_wake import BastankhahGaussian
from scipy.optimize import minimize
import numpy as np
def funC(x, y, c):
"""
Turns on/off the use of wind turbine depending on the value of c.
scipy generates c real values in the range [0, 1] as specified by the bounds including 0.2 etc.
If c is 0.5 or more turbine will be used otherwise turbine will not be used.
"""
x_selected = x[c >= 0.5]
y_selected = y[c >= 0.5]
return (x_selected, y_selected)
def wt_simulation(c):
"""
This is our objective function. It will return the aep=annual energy production in GWh.
We will maximize aep.
"""
site = Hornsrev1Site()
x, y = site.initial_position.T
windTurbines = V80()
wf_model = BastankhahGaussian(site, windTurbines)
x_new, y_new = funC(x, y, c)
# run wind farm simulation
sim_res = wf_model(
x_new, y_new, # wind turbine positions
h=None, # wind turbine heights (defaults to the heights defined in windTurbines)
type=0, # Wind turbine types
wd=None, # Wind direction (defaults to site.default_wd (0,1,...,360 if not overriden))
ws=None, # Wind speed (defaults to site.default_ws (3,4,...,25m/s if not overriden))
)
aep_output = sim_res.aep().sum() # we maximize aep
return -float(aep_output) # negate because of scipy minimize
def solve():
t0 = time.perf_counter()
wt = 80 # for V80
x0 = np.ones(wt) # initial value
bounds = [(0, 1) for _ in range(wt)]
res = minimize(wt_simulation, x0=x0, bounds=bounds)
print(f'success status: {res.success}')
print(f'aep: {-res.fun}') # negate to get the true maximum aep
print(f'c values: {res.x}\n')
print(f'elapse: {round(time.perf_counter() - t0)}s')
# start
solve()
success status: True
aep: 682.0407252944838
c values: [1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
1. 1. 1. 1. 1. 1. 1. 1.]
elapse: 274s
This is using optuna framework considering c as hyperparameter that we are going to optimize to achieve maximum aep (annual energy production). I am using SkoptSampler here from scikit-optimize. Optuna has some samplers that we can use. This will enforce that what scipy saw is also seen by other optimizer.
"""
Additional modules
pip install optuna
pip install scikit-optimize
"""
import time
from py_wake.examples.data.hornsrev1 import V80
from py_wake.examples.data.hornsrev1 import Hornsrev1Site # We work with the Horns Rev 1 site, which comes already set up with PyWake.
from py_wake import BastankhahGaussian
import numpy as np
import optuna
def funC(x, y, c):
"""
Turns on/off the use of wind turbine depending on the value of c.
optuna generates c integer values in the range [0, 1] as specified by the bounds.
If c is 1 turbine will be run otherwise turbine will not be run.
"""
c = np.array(c)
x_selected = x[c == 1]
y_selected = y[c == 1]
return (x_selected, y_selected)
def objective(trial):
"""
This is our objective function. It will return the aep=annual energy production in GWh.
We will maximize aep.
"""
site = Hornsrev1Site()
x, y = site.initial_position.T
windTurbines = V80()
wt = 80 # for v80
# We ask values of c from optuna.
c = []
for i in range(wt):
varname = f'c{i}'
minv, maxv, stepv = 0, 1, 1
c.append(trial.suggest_int(varname, minv, maxv, step=stepv))
wf_model = BastankhahGaussian(site, windTurbines)
x_new, y_new = funC(x, y, c)
# run wind farm simulation
sim_res = wf_model(
x_new, y_new, # wind turbine positions
h=None, # wind turbine heights (defaults to the heights defined in windTurbines)
type=0, # Wind turbine types
wd=None, # Wind direction (defaults to site.default_wd (0,1,...,360 if not overriden))
ws=None, # Wind speed (defaults to site.default_ws (3,4,...,25m/s if not overriden))
)
aep_output = sim_res.aep().sum()
return float(aep_output) # give feedback to optuna of how the c we asks has performed
def optuna_hpo():
t0 = time.perf_counter()
num_trials = 300
sampler = optuna.integration.SkoptSampler()
study = optuna.create_study(sampler=sampler, direction="maximize")
study.optimize(objective, n_trials=num_trials)
print(f"Best params: {study.best_params}")
print(f"Best value: {study.best_value}\n")
print(f'elapse: {round(time.perf_counter() - t0)}s')
# start
optuna_hpo()
...
[I 2021-12-03 16:48:07,935] Trial 299 finished with value: 379.6195135529371 and parameters: {'c0': 0, 'c1': 0, 'c2': 1, 'c3': 1, 'c4': 1, 'c5': 0, 'c6': 1, 'c7': 0, 'c8': 1, 'c9': 0, 'c10': 0, 'c11': 1, 'c12': 0, 'c13': 0, 'c14': 1, 'c15': 0, 'c16': 1, 'c17': 0, 'c18': 1, 'c19': 0, 'c20': 1, 'c21': 0, 'c22': 0, 'c23': 1, 'c24': 1, 'c25': 1, 'c26': 1, 'c27': 1, 'c28': 0, 'c29': 0, 'c30': 0, 'c31': 0, 'c32': 0, 'c33': 1, 'c34': 0, 'c35': 0, 'c36': 1, 'c37': 1, 'c38': 1, 'c39': 0, 'c40': 0, 'c41': 1, 'c42': 0, 'c43': 0, 'c44': 1, 'c45': 1, 'c46': 1, 'c47': 0, 'c48': 0, 'c49': 1, 'c50': 0, 'c51': 1, 'c52': 0, 'c53': 1, 'c54': 1, 'c55': 1, 'c56': 0, 'c57': 1, 'c58': 1, 'c59': 1, 'c60': 1, 'c61': 1, 'c62': 1, 'c63': 0, 'c64': 0, 'c65': 1, 'c66': 0, 'c67': 0, 'c68': 0, 'c69': 1, 'c70': 1, 'c71': 0, 'c72': 1, 'c73': 1, 'c74': 0, 'c75': 1, 'c76': 0, 'c77': 1, 'c78': 1, 'c79': 1}. Best is trial 110 with value: 682.0407252944838.
Best params: {'c0': 1, 'c1': 1, 'c2': 1, 'c3': 1, 'c4': 1, 'c5': 1, 'c6': 1, 'c7': 1, 'c8': 1, 'c9': 1, 'c10': 1, 'c11': 1, 'c12': 1, 'c13': 1, 'c14': 1, 'c15': 1, 'c16': 1, 'c17': 1, 'c18': 1, 'c19': 1, 'c20': 1, 'c21': 1, 'c22': 1, 'c23': 1, 'c24': 1, 'c25': 1, 'c26': 1, 'c27': 1, 'c28': 1, 'c29': 1, 'c30': 1, 'c31': 1, 'c32': 1, 'c33': 1, 'c34': 1, 'c35': 1, 'c36': 1, 'c37': 1, 'c38': 1, 'c39': 1, 'c40': 1, 'c41': 1, 'c42': 1, 'c43': 1, 'c44': 1, 'c45': 1, 'c46': 1, 'c47': 1, 'c48': 1, 'c49': 1, 'c50': 1, 'c51': 1, 'c52': 1, 'c53': 1, 'c54': 1, 'c55': 1, 'c56': 1, 'c57': 1, 'c58': 1, 'c59': 1, 'c60': 1, 'c61': 1, 'c62': 1, 'c63': 1, 'c64': 1, 'c65': 1, 'c66': 1, 'c67': 1, 'c68': 1, 'c69': 1, 'c70': 1, 'c71': 1, 'c72': 1, 'c73': 1, 'c74': 1, 'c75': 1, 'c76': 1, 'c77': 1, 'c78': 1, 'c79': 1}
Best value: 682.0407252944838
elapse: 4533s
It found the maximum aep of 682.0407252944838 as early as trial 110 out of 300 trials.
300 trials are completed after 4533s or 1.25hr. The best parameters are all 1 meaning all turbines must run to get maximum aep.