Scipy 中的优化

Posted 2023-03-12

【中文标题】Scipy 中的优化

【英文标题】：Optimization in Scipy

【发布时间】：2022-01-11 00:43:52

【问题描述】：

我想在下面的代码中添加一些约束，我想在其中使用 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 ，其中我有另一个必须优化的返回，所以我不知道如何返回不受优化影响的 TI。

【问题讨论】：

这个ti_eff只有在仿真后才可用。我们不能在模拟之前输入 ti_eff。我们的目标是最大 aep。在我看来，scipy 无法做到这一点。但是 optuna 超参数调谐器可以。我们的想法是我们首先模拟，然后获取所有涡轮机的 ti_eff，如果甚至有一个涡轮机的 ti_eff 为 0.2 及以上，我们告诉 optuna aep 为零，这样它将尝试找到所有涡轮机的最大 aep，其中 ti_eff低于 0.2。

如果您在函数评估中有模拟步骤，则需要将求解器限制为“无导数优化”（DFO）类。进行有限差分不是一个好主意。

@ferdy 所以你的意思是我必须使用 Optuna 来处理这种类型的约束？如果是，则在模拟后的每次迭代中，如果 TI 大于 0.2，则使该 wt=0 的 aep 因此该特定 wt 的 c 必须更改为 0，它将作为下一次迭代的初始猜测（ c(wt=i)=0)。但是 optuna 可以这样做吗？

@ErwinKalvelagen 你能解释一下吗？因为我没有太多的编程和优化经验，所以有些概念是新的

使用模拟你不会得到渐变。大多数非线性求解器依赖于梯度。如果不可用，则可以使用有限差分来近似它们。这假设函数评估很便宜。在模拟的情况下，情况并非如此。所以 DFO 求解器是可以使用的求解器类型：它们不需要梯度，也不做有限差分。一个很好的参考是：link.springer.com/article/10.1007%252Fs10898-012-9951-y。我不完全理解你的约束。首先制定一个数学模型有很大帮助。

【参考方案1】：

这是处理 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