Optimization of Wind Turbine plant in Scipy

Question

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?

Is the simulation output deterministic? There are only 2 values in c. Why not compare the power output when c=0 and when c=1, then select the c with highest power output. — ferdy, Dec 02 '21 at 14:42
i need a combination of c=0 and c=1, I mean I need to switch off and switch on some wt — sadra, Dec 02 '21 at 14:49
If I will run `wf_model(x, y)` 10 times (without c yet), will I get the same output power? I am trying to install pywake but not done yet. I will try scipy, and hyperparameter optimizers. — ferdy, Dec 02 '21 at 15:07
if you want to delete C you must run wf_model(x, y) with x ,y not x_new,y_new. since x,y are for the original layout and x/y_new are used when we implement c and switch off some wind turbines. — sadra, Dec 02 '21 at 15:11
Yes I get that. By the way is the `import function` available? — ferdy, Dec 02 '21 at 15:14
thanks a lot for your help. here is the function: def funC(x,y,c): x_selected = x[c==1] y_selected = y[c==1] return (x_selected, y_selected) — sadra, Dec 02 '21 at 15:17
I run that code (without c) 10x and it returns `682.0407252944838` 10x too, so this function is deterministic. — ferdy, Dec 02 '21 at 16:09
what do you mean by deterministic? maybe i did not explain well the problem. i must turn off some wind turbines,absolutely if I can run all of them I can get the highest result. but I must to switch off some of them. now I want to know if I change the c and indeed change the number of on/off WT, which combination gives better results which is the highest power? — sadra, Dec 02 '21 at 16:37
Deterministic means if I run it 100 times or so, the answer is the same 100 times too. I run it once the answer is 682.0407252944838, then again the answer is the same 682.0407252944838. I just want to make sure that when c is introduced in the system the output power is only affected by varying the c. I already tried it with scipy and to deliver maximum power all turbines needs to run. Of course because our objective is the maximum power output. Will post the code later. — ferdy, Dec 02 '21 at 22:21
let me clarify this you said `now I want to know if I change the c and indeed change the number of on/off WT, which combination gives better results which is the highest power?` Can you give an example of what that is about? — ferdy, Dec 03 '21 at 11:05
yes exactly, indeed i want to complete the code in future and I want to add some new constraint to the code, but for the first step I want to know if I have to switch off some wt which combination gives better result? for sure if I can use all the turbines I can get the highest output, but as I said I want to add some limitations in future so I can not use all the turbines so I need to switch off some of them . as a initial guess you can consider any combinations for c values it just must be 0 or 1 — sadra, Dec 03 '21 at 11:37
All right, could you accept my answer for now, given two approaches, scipy and optuna? — ferdy, Dec 03 '21 at 13:06
yes, i think it is completely correct, just I have some questions regarding the code in Scipy: if I want to add other constraints what should I do? and also how can I implement penalty function to the minimization procedure? — sadra, Dec 03 '21 at 15:36
Good ref is [here](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html), also some examples [here](https://scipy-lectures.org/advanced/mathematical_optimization/). Some people here are also willing to help. I thought about changing the objective, instead of power you can use the price or amount of money it generates. For example objective can be the sale of producing the power less operating and maintenance cost, repairs etc. — ferdy, Dec 03 '21 at 16:49
sorry, I faced another question, I really appreciate you if you can also help me with this. I want to define a constraint for the optimization. sim_res.TI_eff gives the turbulence intensity of each wd and ws. the problem is that for constraint I need to define a function for it, but we define the code in the way that we can only return .aep so how can I define a function for TI? I mean how can I use the output of wt_simulation(c) to check TI_eff ? — sadra, Dec 04 '21 at 21:49
I guess we need to close this issue first, then you need to create another question. — ferdy, Dec 04 '21 at 23:42

ferdy · Accepted Answer · 2021-12-03T09:35:05.307

(1) SCIPY

Code

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

Output

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

(2) OPTUNA

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.

Code

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

Output

...

[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.

Optimization of Wind Turbine plant in Scipy

1 Answers1

(1) SCIPY

Code

Output

(2) OPTUNA

Code

Output