Optimization of wind farm using Penalty function in Scipy

Question

In the following code I want to optimize a wind farm using a penalty function.

Using the first function(newsite), I have defined the wind turbines numbers and layout. Then in the next function, after importing x0(c=x0=initial guess), for each range of 10 wind directions (wd) I took the c values for the mean wd of each range. For instance, for wd:[0,10] mean value is 5 and I took c values of wd=5 and put it for all wd in the range[0,10] and for each wind speed(ws). I have to mention that c is the value that shows that wind turbines are off or on(c=0 means wt is off). then I have defined operating according to the c, which means that if operating is 0,c=0 and that wt is off.

Then I defined the penalty function to optimize power output. indeed wherever TI_eff>0.14, I need to implement a penalty function so this function must be subtracted from the original power output. For instance, if sim_res.TI_eff[1][2][3] > 0.14, so I need to apply penalty function so curr_func[1][2][3]=sim_res.Power[1][2][3]-10000*(sim_res.TI_eff[1][2][3]-0.14)**2.

The problem is that I run this code but it did not give me any results and I waited for long hours, I think it was stuck in a loop that could not reach converge. so I want to know what is the problem?

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

There are a number of things in your approach that are either wrong or incomprehensible to me. Just for fun I have tried your code. A few observations:

1. Your set of parameters (optimization variables) has a shape of (4, 360, 23), i.e. you are looking at 33,120 parameters. There is no nonlinear optimization algorithm that is going to give you any meaningful answer to a problem that big. Ever. But then again, you shouldn't be looking at SciPy optimize if your optimization variables should only assume 0/1 values.

2. Calling SciPy minimize like this:

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

   Is going to select a nonlinear optimizer between BFGS, L-BFGS-B or SLSQP (according to the documentation at https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html)

   Those algorithms are gradient-based, and since you're not providing a gradient of your objective function SciPy is going to calculate them numerically. Good luck with that when you have 33,000 parameters. Never going to finish.

3. At the beginning of your objective function you are doing this:

   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]

   I don't understand why you're doing it but you are overriding the input values of c coming from the optimizer.

4. Your objective function takes 20-25 seconds to evaluate on my powerful workstation. Even if you had only 10-15 optimization parameters, it would take you several days to get any answer out of an optimizer. You have 33,000+ variables. No way.

I don't know why you are doing this and why you're doing it the way you're doing it. You should rethink your approach.