在 Scipy 中使用 Penalty 函数优化风电场

Question

在下面的代码中，我想使用惩罚函数优化风电场。

使用第一个函数(newsite)，我已经定义了风力涡轮机的数量和布局。然后在下一个函数中，在导入 x0(c=x0=initial guess) 之后，对于 10 个风向的每个范围 (wd)，我将 c 值作为每个范围的平均 wd。例如，对于wd:[0,10]，平均值为5，我将c 的值设为wd=5，并将其用于[0,10] 范围内的所有wd 和每个风速（ws）。我不得不提到c 是表明风力涡轮机关闭或开启的值（c=0 表示wt 关闭）。那么我根据c定义了operating，也就是说如果operating为0，c=0，则wt为off。

然后我定义了惩罚函数来优化功率输出。事实上，无论TI_eff>0.14，我都需要实现一个惩罚函数，所以这个函数必须从原始功率输出中减去。比如sim_res.TI_eff[1][2][3] > 0.14，那么我需要应用惩罚函数所以curr_func[1][2][3]=sim_res.Power[1][2][3]-10000*(sim_res.TI_eff[1][2][3]-0.14)**2。

问题是我运行了这段代码，但它没有给我任何结果，我等了很长时间，我认为它卡在了一个无法收敛的循环中。所以我想知道是什么问题？

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 py_wake.turbulence_models import GCLTurbulence
from py_wake.deflection_models.jimenez import JimenezWakeDeflection
from scipy.optimize import minimize
from py_wake.wind_turbines.power_ct_functions import PowerCtFunctionList, PowerCtTabular
import numpy as np

def newSite(x,y):
    xNew=np.array([x[0]+560*i for i in range(4)])
    yNew=np.array([y[0]+560*i for i in range(4)])    
    x_newsite=np.array([xNew[0],xNew[0],xNew[0],xNew[1]])
    y_newsite=np.array([yNew[0],yNew[1],yNew[2],yNew[0]])

    return (x_newsite,y_newsite)



def wt_simulation(c):
    
    c = c.reshape(4,360,23)
    site = Hornsrev1Site()
    x, y = site.initial_position.T
    x_newsite,y_newsite=newSite(x,y)
    windTurbines = V80()
    
    
    for item in range(4):
        for j in range(10,370,10):
            for i in range(j-10,j):
                c[item][i]=c[item][j-5]
          
    windTurbines.powerCtFunction = PowerCtFunctionList(
    key='operating',
    powerCtFunction_lst=[PowerCtTabular(ws=[0, 100], power=[0, 0], power_unit='w', ct=[0, 0]), # 0=No power and ct
                         windTurbines.powerCtFunction], # 1=Normal operation
    default_value=1)

    operating = np.ones((4,360,23)) # shape=(#wt,wd,ws)
    operating[c <= 0.5]=0
    
    wf_model = BastankhahGaussian(site, windTurbines,deflectionModel=JimenezWakeDeflection(),turbulenceModel=GCLTurbulence())

    # run wind farm simulation
    sim_res = wf_model(
        x_newsite, y_newsite, # wind turbine positions
        h=None, # wind turbine heights (defaults to the heights defined in windTurbines)
        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))
        operating=operating
   )
    
    curr_func=np.ones((4,360,23))
    for i in range(4):
        for l in range(360):
            for k in range(23):
                  if sim_res.TI_eff[i][l][k]-0.14 > 0 :
                      curr_func[i][l][k]=sim_res.Power[i][l][k]-10000*(sim_res.TI_eff[i][l][k]-0.14)**2
                  else:
                      curr_func[i][l][k]=sim_res.Power[i][l][k]
    
    return -float(np.sum(curr_func))  # negative because of scipy minimize
         
t0 = time.perf_counter()

        
def solve():
    
    wt =4  # for V80  
    wd=360
    ws=23
    x0 = np.ones((wt,wd,ws)).reshape(-1)  # initial value for c
    b=(0,1)
    bounds=np.full((wt,wd,ws,2),b).reshape(-1, 2)

    
    res = minimize(wt_simulation, x0=x0, bounds=bounds)

    return res
    
res=solve()

print(f'success status: {res.success}')
print(f'aep: {-res.fun}')  # negative to get the true maximum aep
print(f'c values: {res.x}\n')
print(f'elapse: {round(time.perf_counter() - t0)}s')  

sim_res=wt_simulation(res.x)

Answer 1

您的方法中有很多地方对我来说要么是错误的，要么是无法理解的。只是为了好玩，我试过你的代码。几点观察：

1. 您的参数集（优化变量）的形状为 (4, 360, 23)，即您正在查看 33,120 个参数。对于这么大的问题，没有任何非线性优化算法可以为您提供任何有意义的答案。曾经。但话又说回来，如果您的优化变量应该只假设 0/1 值，那么您不应该查看 SciPy 优化。

2. 像这样调用 SciPy 最小化：

   res = minimize(wt_simulation, x0=x0, bounds=bounds)

   将在 BFGS、L-BFGS-B 或 SLSQP 之间选择非线性优化器（根据https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html 的文档）

   这些算法是基于梯度的，由于您没有提供目标函数的梯度，SciPy 将以数值方式计算它们。当你有 33,000 个参数时，祝你好运。永远不会完成。

3. 在您的目标函数开始时，您正在这样做：

   对于范围 (4) 中的项目： 对于范围内的 j (10,370,10)： 对于范围内的 i (j-10,j)： c[项目][i]=c[项目][j-5]

   我不明白您为什么要这样做，但您正在覆盖来自优化器的 c 的输入值。

4. 在我强大的工作站上评估您的目标函数需要 20-25 秒。即使您只有 10-15 个优化参数，您也需要几天时间才能从优化器中得到任何答案。您有 33,000 多个变量。没办法。

我不知道您为什么要这样做，以及为什么要按照自己的方式进行操作。你应该重新考虑你的方法。