Scipy 中的优化
Optimization in Scipy
我想在下面的代码中添加一些约束，我想在其中使用 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()
现在我想添加一个约束，其中湍流强度：sim_res_TI_eff 每个风力涡轮机 (wt) 对于每个风速 (was) 和每个风向 (wd) 必须低于某个值（例如 0.2）。我必须补充一点，例如 sim_res.TI_eff.sel(wt=1) 给出了每个 wd 的 TI，并且是风力涡轮机 #1。问题是我需要使用函数 wt_simulation，其中我有另一个必须优化的 return，所以我不知道如何 return TI 不受优化的影响。
这是处理 ti_eff 的一种方法。
"""
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
TIEFF_THRESHOLD = 0.2
def ok_tieff_constraint(tieff):
"""
Returns True if tieff is below threshold otherwise False.
"""
if np.any(tieff >= TIEFF_THRESHOLD):
return False
return True
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.
"""
islogging = True
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))
)
if islogging:
print(sim_res)
aep_output = float(sim_res.aep().sum()) # we maximize aep
# Constraint on ti_eff, if constraint is violated we set aep to zero.
if not ok_tieff_constraint(sim_res.TI_eff):
aep_output = 0
return -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()
输出：
...
<xarray.SimulationResult>
Dimensions: (wt: 80, wd: 360, ws: 23)
Coordinates:
* wt (wt) int32 0 1 2 3 4 5 6 7 8 ... 72 73 74 75 76 77 78 79
* wd (wd) int32 0 1 2 3 4 5 6 7 ... 353 354 355 356 357 358 359
* ws (ws) int32 3 4 5 6 7 8 9 10 11 ... 18 19 20 21 22 23 24 25
x (wt) float64 4.24e+05 4.24e+05 ... 4.294e+05 4.295e+05
y (wt) float64 6.151e+06 6.151e+06 ... 6.148e+06 6.148e+06
h (wt) float64 70.0 70.0 70.0 70.0 ... 70.0 70.0 70.0 70.0
type (wt) int32 0 0 0 0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 0 0 0
Data variables: (12/15)
WS_eff (wt, wd, ws) float64 3.0 4.0 5.0 6.0 ... 22.87 23.88 24.89
TI_eff (wt, wd, ws) float64 0.1 0.1 0.1 0.1 ... 0.1 0.1 0.1 0.1
Power (wt, wd, ws) float64 0.0 6.66e+04 1.54e+05 ... 2e+06 2e+06
CT (wt, wd, ws) float64 0.0 0.818 0.806 ... 0.06084 0.05377
WS (ws) int32 3 4 5 6 7 8 9 10 11 ... 18 19 20 21 22 23 24 25
WD (wd) int32 0 1 2 3 4 5 6 7 ... 353 354 355 356 357 358 359
... ...
Weibull_A (wd) float64 9.177 9.177 9.177 9.177 ... 9.177 9.177 9.177
Weibull_k (wd) float64 2.393 2.393 2.393 2.393 ... 2.393 2.393 2.393
Sector_frequency (wd) float64 0.001199 0.001199 ... 0.001199 0.001199
P (wd, ws) float64 6.147e-05 8.559e-05 ... 2.193e-08
tilt (wt, wd, ws) float64 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0
yaw (wt, wd, ws) float64 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0
Attributes:
wd_bin_size: 1
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: 273s
我想在下面的代码中添加一些约束，我想在其中使用 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()
现在我想添加一个约束，其中湍流强度：sim_res_TI_eff 每个风力涡轮机 (wt) 对于每个风速 (was) 和每个风向 (wd) 必须低于某个值（例如 0.2）。我必须补充一点，例如 sim_res.TI_eff.sel(wt=1) 给出了每个 wd 的 TI，并且是风力涡轮机 #1。问题是我需要使用函数 wt_simulation，其中我有另一个必须优化的 return，所以我不知道如何 return TI 不受优化的影响。
这是处理 ti_eff 的一种方法。
"""
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
TIEFF_THRESHOLD = 0.2
def ok_tieff_constraint(tieff):
"""
Returns True if tieff is below threshold otherwise False.
"""
if np.any(tieff >= TIEFF_THRESHOLD):
return False
return True
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.
"""
islogging = True
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))
)
if islogging:
print(sim_res)
aep_output = float(sim_res.aep().sum()) # we maximize aep
# Constraint on ti_eff, if constraint is violated we set aep to zero.
if not ok_tieff_constraint(sim_res.TI_eff):
aep_output = 0
return -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()
输出：
...
<xarray.SimulationResult>
Dimensions: (wt: 80, wd: 360, ws: 23)
Coordinates:
* wt (wt) int32 0 1 2 3 4 5 6 7 8 ... 72 73 74 75 76 77 78 79
* wd (wd) int32 0 1 2 3 4 5 6 7 ... 353 354 355 356 357 358 359
* ws (ws) int32 3 4 5 6 7 8 9 10 11 ... 18 19 20 21 22 23 24 25
x (wt) float64 4.24e+05 4.24e+05 ... 4.294e+05 4.295e+05
y (wt) float64 6.151e+06 6.151e+06 ... 6.148e+06 6.148e+06
h (wt) float64 70.0 70.0 70.0 70.0 ... 70.0 70.0 70.0 70.0
type (wt) int32 0 0 0 0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 0 0 0
Data variables: (12/15)
WS_eff (wt, wd, ws) float64 3.0 4.0 5.0 6.0 ... 22.87 23.88 24.89
TI_eff (wt, wd, ws) float64 0.1 0.1 0.1 0.1 ... 0.1 0.1 0.1 0.1
Power (wt, wd, ws) float64 0.0 6.66e+04 1.54e+05 ... 2e+06 2e+06
CT (wt, wd, ws) float64 0.0 0.818 0.806 ... 0.06084 0.05377
WS (ws) int32 3 4 5 6 7 8 9 10 11 ... 18 19 20 21 22 23 24 25
WD (wd) int32 0 1 2 3 4 5 6 7 ... 353 354 355 356 357 358 359
... ...
Weibull_A (wd) float64 9.177 9.177 9.177 9.177 ... 9.177 9.177 9.177
Weibull_k (wd) float64 2.393 2.393 2.393 2.393 ... 2.393 2.393 2.393
Sector_frequency (wd) float64 0.001199 0.001199 ... 0.001199 0.001199
P (wd, ws) float64 6.147e-05 8.559e-05 ... 2.193e-08
tilt (wt, wd, ws) float64 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0
yaw (wt, wd, ws) float64 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0
Attributes:
wd_bin_size: 1
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: 273s