Wind Farm AEP Optimization in SciPy – Horns Rev 1 Site

import os
os.chdir('YOURPATH')    # Set working directory
os. getcwd()

import py_wake
import scipy
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 py_wake.examples.data.hornsrev1 import wt_x, wt_y
from scipy.optimize import minimize
import optuna
import matplotlib.pyplot as plt
import numpy as np

site = Hornsrev1Site()
wt = V80()

site.plot_wd_distribution(n_wd=12);

plt.figure()
wt.plot(wt_x, wt_y)
plt.xlabel('x [m]')
plt.ylabel('y [m]')

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: 122s

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()

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: 3173s