基于罚函数的风电场优化

Question

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

使用第一个函数(newsite)，我定义了风力涡轮机的数量和布局。然后在下一个函数中，在导入x0(c=x0=初始猜测)之后，对于每10个风向(wd)的范围，我取每个范围的平均wd的c值。例如，对于wd:[0,10]平均值是5，我取了wd=5的c值，并将其用于范围[0,10]中的所有wd和每个风速(ws)。我必须指出的是，c是指风力涡轮机关闭或开启的值(c=0的意思是wt关闭)。然后，我根据operating定义了c，这意味着如果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优化。

1. 调用SciPy minimize如下所示：

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

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

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

在目标函数的开头，

1. 正在这样做：

对于范围内的项目(4)：j在范围内(10,370，10)：i在范围内(j-10，j)：citem=citem

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

1. 您的目标函数需要20-25秒才能在我强大的工作站上进行评估。即使您只有10-15个优化参数，从优化器获得任何答案也需要几天的时间。您有33,000+变量。不可能，

我不知道你为什么要这么做，为什么你要这样做。你应该重新考虑你的方法。