Scipy 中的风力发电机组优化
Optimization of Wind Turbine plant in Scipy
我正在做一个项目（使用用户定义的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 值。我的意思是我们想要切换 off/on 不同的风力涡轮机并查看输出功率。我知道优化值是当 c 中的所有参数都为一个时，但有一些限制，所以我需要使用惩罚函数。那么请您帮助我如何通过使用 c 和 penalty 来优化发电量？
(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 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 有一些我们可以使用的采样器。这将强制 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。
我正在做一个项目（使用用户定义的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 值。我的意思是我们想要切换 off/on 不同的风力涡轮机并查看输出功率。我知道优化值是当 c 中的所有参数都为一个时，但有一些限制，所以我需要使用惩罚函数。那么请您帮助我如何通过使用 c 和 penalty 来优化发电量？
(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 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 有一些我们可以使用的采样器。这将强制 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。