Scipy中风力涡轮机设备的优化

Question

我正在做一个项目（使用 pywake 库，它是用户定义的库），我编写了以下代码：

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

为了让代码清楚一点，我们输入了一些风电场的数据，然后进行仿真，得到风力发电机的功率输出(print(sim_res.aep().sum()) )。

现在我定义了一个新变量（c），其中我们有两个值 0 和 1。如果 c=1，我们说该风力涡轮机是一个，否则它已关闭，因此发电量将减少。

现在通过使用定义的脚本，我想通过使用惩罚函数在 Scipy 中进行优化：我想通过更改 c 值来最大化发电量。我的意思是我们想关闭/打开不同的风力涡轮机并查看输出功率。我知道优化值是当 c 中的所有参数都为一个但有一些限制时，我需要使用惩罚函数。那么请您帮助我如何通过使用 c 和惩罚来优化发电量？

Answer 1

(1) SCIPY

代码

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

(2) OPTUNA

这是使用optuna framework 考虑c 作为超参数，我们将对其进行优化以实现最大aep（年发电量）。我在这里使用来自scikit-optimize 的 SkoptSampler。 Optuna 有一些我们可以使用的采样器。这将强制 scipy 看到的内容也被其他优化器看到。

代码

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

早在300次试验中的第110次就发现最大aep为682.0407252944838。

300 次试验在 4533 秒或 1.25 小时后完成。最佳参数均为 1，这意味着所有涡轮机都必须运行才能获得最大 aep。