Using the CMA-ES algorithm to optimize the position of wind turbines

Question

I have an implementation of an optimization algorithm called the Covariance Matrix Adaptation Evolutionary Strategy (CMA-ES). I am using this algorithm to optimize the position of wind turbines inside a wind farm, but if I place more than a certain number of turbines (equivalent to a high value of N at the beginning of the code) after around 50 iterations of the while loop, my RAM memory is filled. I have tried using deletion of certain variables after they are being used and using slots in my class declaration but to no noticeable effect.

Could you take a look at the code and tell me if there is something very memory consuming that I'm not seeing?

This is a stand-alone, coupled with the Rosenbrock function to optimize.

def frosenbrock(x):
    x1=np.delete(x,len(x)-1);
    x2=np.delete(x,0);
    f = 100*sum(((x1**2) - (x2))**2) + sum((x1-1)**2);
    return f;

def optimise_CMA(self,initial_guess):   
    #--------------INITIALIZATION----------------------
    N=20 #Number of variables, problem dimension
    #Container for coordinates called xmean
    xmean=[]                                  #Objective variables (coordinates)initial point
    for i in range(int(N)):
        xmean.append([random.uniform(0,1)]);
    xmean=np.asarray(xmean,dtype=np.float_)
    sigma=0.5                              #Coordinate wise standard deviation
    stopfitness = 1e-10                    #stop if fitness < stopfitness (minimization)
    stopeval =1e3*N**2                     #stop after this number of evaluations ()

    #Strategy parameter setting: Selection
    lambd=int(4.0+floor(3.0*log(N)))                #population size, offspring number
    mu=(lambd/2)                           #number of parents/points for recombination
    weights=np.array([])
    for j in range(mu):
        weights=np.append(weights,log(mu+0.5)-log(j+1))#muXone array for weighted recombination
    mu = floor(mu)
    weights=weights/sum(weights)
    mueff=sum(weights)**2/sum(weights**2)
    weights_array=weights.reshape([mu,1])  

    #Strategy parameter setting: Adaptation
    cc=(4+mueff/N)/(N+4 + 2*mueff/N)      #time constant for cumulation for C
    cs=(mueff+2)/(N+mueff+5)              #t-const for cumulation for sigma control 
    c1=2/((N+1.3)**2+mueff)               #learning rate for rank-one update of C
    cmu=min(1-c1,2*(mueff-2+1/mueff)/((N+2)**2+mueff)) #and for rank-mu update
    damps =1 +2*max(0,sqrt((mueff-1)/(N+1))-1)+cs  #damping for sigma, usually close to 1

    #Initialize dynamic (internal) strategy parameters and constants
    pc=np.zeros((N,1),dtype=np.float_)                      #evolution path for C  
    ps=np.zeros((N,1),dtype=np.float_)                      #evolution path for sigma 
    D=np.ones((N,1),dtype=np.float_)                        #diagonal D defines the scaling  
    B =np.eye(N,dtype=np.float_)                #B defines the coordinate system
    D_sqd=D**2           
    diag_D_sqd=(np.diag(D_sqd[:,0]))       #Generate diagonal matrix with D_sqd in the diagonal, the rest are zeros 
    trans_B=B.transpose()                 #Tranpose matrix of B 
    C=np.dot(np.dot(B,diag_D_sqd),trans_B)#Covariance matrix C
    nega_D=D**(-1)                        #Elementwise to power -1 
    diag_D_neg=(np.diag(nega_D[:,0]))     #Generate diagonal matrix with nega_D in diagonal, rest zeros
    invsqrtC =np.dot(np.dot(B,diag_D_neg),trans_B)#C^-1/2 
    eigeneval = 0                         #track update of B and D
    chiN=N**0.5*(1-1/(4*N)+1/(21*N**2))   #expectation of ||N(0,I)|| == norm(randn(N,1)) 

    #---------GENERATION LOOP---------------
    counteval=0
    arx=[]
    while counteval<stopeval:
        #Generate and evaluate lamba offspring
        itera=1
        arfitness=[]
        for l in range(lambd):
            offspring=[]                  #Create a container for the offspring 
            offspring=xmean+np.dot(sigma*B,(D*np.random.standard_normal((N,1)))) #m + sig * Normal(0,C)
            if itera==1:
                arx=offspring
            else:
                arx=np.hstack((arx,offspring))
            arfitness.append(frosenbrock(offspring))          #EVALUATE OBJ FUNCTION
            counteval=counteval+1
            itera=itera+1   
        #Sort by fitness and compute weighted mean into xmean
        ordered=[]
        ordered=sorted(enumerate(arfitness), key=lambda x: x[1])    #minimization, list with (Index,Fitness) elements
        xold=xmean
        best_off_indexes=[]
        for i in range(int(mu)):
            best_off_indexes.append(ordered[i][0]) #List of best indexes of mu offspring
        recomb=np.zeros([N,mu])
        cont=0
        for index in best_off_indexes:
            recomb[:,cont]=arx[:,index]
            cont=cont+1
        xmean=np.dot(recomb,weights_array)

        #Cumulation: update evolution paths
        ps=(1-cs)*ps+sqrt(cs*(2-cs)*mueff)*np.dot(invsqrtC,(xmean-xold))*(1/sigma)
        hsig=np.linalg.norm(ps)/sqrt(1-(1-cs)**(2*counteval/lambd))/chiN < 1.4 + 2/(N+1)
        if hsig==True:
            hsig=1
        else:
            hsig=0
        pc=(1-cc)*pc+hsig*sqrt(cc*(2-cc)*mueff)*(xmean-xold)/sigma

        #Adapt covariance matrix C
        artmp=(1/sigma)*(recomb-np.tile(xold,(1,mu)))
        weights_diag=np.diag(weights_array[:,0])
        C=(1-c1-cmu)*C+c1*(np.dot(pc,pc.transpose())+((1-hsig)*cc*(2-cc)*C))+cmu*np.dot(artmp,np.dot(weights_diag,artmp.transpose()))

        #Adapt step size sigma
        sigma=sigma*exp((cs/damps)*(np.linalg.norm(ps)/chiN - 1))

        #Decomposition of C into B*diag(D.^2)*B' (diagonalization)
        if counteval-eigeneval > lambd/(c1+cmu)/N/10:
            eigeneval=counteval                    #to achieve O(N^2)
            C=np.triu(C)+np.triu(C,1).transpose()  #enforce symmetry
            B=np.linalg.eig(C)[1]                  #eigen decomposition, B==normalized eigenvectors
            diag_D=np.linalg.eig(C)[0]
            diag_D=diag_D.reshape([N,1])
            D=diag_D**0.5                       #D is a vector of standard deviations now
            D_inv=D**(-1)
            D_inv=np.diag(D_inv[:,0])
            invsqrtC=np.dot(B,np.dot(D_inv,B.transpose()))
            print ordered[0][1];               # Uncomment to see convergence in the console
        #Break, if fitness is good enough or condition exceeds 1e14, better termination methods are advisable           
        del arfitness
        if ordered[0][1]<=stopfitness or max(D)>1e7*min(D):
            break
        else:
            del arx
            del ordered
        #Return best point at last iteration
    xmin=arx[:,best_off_indexes[0]]
    return xmin

Answer 1

1. First off, style. Your code is... not very well styled. You have many PEP8 violations. Here's a list of those violations.

   1. You have no space between operators/variable declarations. You should have one space between operators, like this: a + b, a == b, or var = 10.
   2. Secondly, your variable names need to be more descriptive. Names like cont, arx, or l are not good. Names should be not too short, not too long, and as descriptive as possible.
   3. You need two blank lines between functions.

2. I'm not quite sure on this, but I think you issue with memory stems from the fact that you're doing many intensive calculations. Python is not Mathematica or R, which means that it isn't built specifically to do calculations like this. You could use a mathematical library, like numpy, as well. See if you can shorten or simplify some of your expressions.

3. You're also doing a ton of looping, which can hang the program, especially with all those calculations. See where loops can be shortened.
4. Finally, it seems like you're allocating quite a bit of items to lists. If you can remove lists after they're needed, this may speed up your execution time.

Answer 2

import numpy as np
import math
def frosenbrock(x):

    x1 = x[0][0]
    x2 = x[1][0]
    val = 100. * (x2 - x1 ** 2.) ** 2. + (1 - x1) ** 2.
    return val;

def optimise_CMA(initial_guess):   
    #--------------INITIALIZATION----------------------
    N=2 #Number of variables, problem dimension
    #Container for coordinates called xmean
    xmean=initial_guess                           #Objective variables (coordinates)initial point
    if len(xmean) == 0:
        for i in range(int(N)):
            xmean.append([random.uniform(0,1)]);

    xmean=np.array(xmean,dtype=np.float_).reshape(N,1)
    sigma=0.5                              #Coordinate wise standard deviation
    stopfitness = 1e-10                    #stop if fitness < stopfitness (minimization)
#     stopeval =1e3*N**2                     #stop after this number of evaluations ()
    stopeval = 100                    #stop after this number of evaluations ()


    #Strategy parameter setting: Selection
    lambd=int(4.0+math.floor(3.0*np.log(N)))                #population size, offspring number
    mu=int(lambd/2)                           #number of parents/points for recombination
    weights=np.array([])
    for j in range(mu):
        weights=np.append(weights,np.log(mu+0.5)-np.log(j+1))#muXone array for weighted recombination
    mu = math.floor(mu)
    weights=weights/sum(weights)
    mueff=sum(weights)**2/sum(weights**2)
    weights_array=weights.reshape([mu,1])  

    #Strategy parameter setting: Adaptation
    cc=(4+mueff/N)/(N+4 + 2*mueff/N)      #time constant for cumulation for C
    cs=(mueff+2)/(N+mueff+5)              #t-const for cumulation for sigma control 
    c1=2/((N+1.3)**2+mueff)               #learning rate for rank-one update of C
    cmu=min(1-c1,2*(mueff-2+1/mueff)/((N+2)**2+mueff)) #and for rank-mu update
    damps =1 +2*max(0,np.sqrt((mueff-1)/(N+1))-1)+cs  #damping for sigma, usually close to 1

    #Initialize dynamic (internal) strategy parameters and constants
    pc=np.zeros((N,1),dtype=np.float_)                      #evolution path for C  
    ps=np.zeros((N,1),dtype=np.float_)                      #evolution path for sigma 
    D=np.ones((N,1),dtype=np.float_)                        #diagonal D defines the scaling  
    B =np.eye(N,dtype=np.float_)                #B defines the coordinate system
    D_sqd=D**2           
    diag_D_sqd=(np.diag(D_sqd[:,0]))       #Generate diagonal matrix with D_sqd in the diagonal, the rest are zeros 
    trans_B=B.transpose()                 #Tranpose matrix of B 
    C=np.dot(np.dot(B,diag_D_sqd),trans_B)#Covariance matrix C
    nega_D=D**(-1)                        #Elementwise to power -1 
    diag_D_neg=(np.diag(nega_D[:,0]))     #Generate diagonal matrix with nega_D in diagonal, rest zeros
    invsqrtC =np.dot(np.dot(B,diag_D_neg),trans_B)#C^-1/2 
    eigeneval = 0                         #track update of B and D
    chiN=N**0.5*(1-1/(4*N)+1/(21*N**2))   #expectation of ||N(0,I)|| == norm(randn(N,1)) 

    #---------GENERATION LOOP---------------
    counteval=0
    arx=[]
    best_fitness = 100000000
    while counteval<stopeval:
        #Generate and evaluate lamba offspring
        itera=1
        arfitness=[]
        for l in range(lambd):
            offspring=[]                  #Create a container for the offspring 
            offspring = xmean + sigma * np.dot(B,(np.multiply(D,np.random.standard_normal((N,1))))) #m + sig * Normal(0,C)
            if itera==1:
                arx=offspring
            else:
                arx=np.hstack((arx,offspring))
            arfitness.append(frosenbrock(offspring))          #EVALUATE OBJ FUNCTION
            counteval=counteval+1
            itera=itera+1   
        #Sort by fitness and compute weighted mean into xmean
        ordered=[]
        ordered=sorted(enumerate(arfitness), key=lambda x: x[1])    #minimization, list with (Index,Fitness) elements
        xold=xmean
        best_off_indexes=[]
        for i in range(int(mu)):
            best_off_indexes.append(ordered[i][0]) #List of best indexes of mu offspring
        recomb=np.zeros([N,mu])
        cont=0
        for index in best_off_indexes:
            recomb[:,cont]=arx[:,index]
            cont=cont+1
        xmean=np.dot(recomb,weights_array)

        #Cumulation: update evolution paths
        ps=(1-cs)*ps+np.sqrt(cs*(2-cs)*mueff)*np.dot(invsqrtC,(xmean-xold))*(1/sigma)
        hsig=np.linalg.norm(ps)/np.sqrt(1-(1-cs)**(2*counteval/lambd))/chiN < 1.4 + 2/(N+1)
        if hsig==True:
            hsig=1
        else:
            hsig=0
        pc=(1-cc)*pc+hsig*np.sqrt(cc*(2-cc)*mueff)*(xmean-xold)/sigma

        #Adapt covariance matrix C
        artmp=(1/sigma)*(recomb-np.tile(xold,(1,mu)))
        weights_diag=np.diag(weights_array[:,0])
        C=(1-c1-cmu)*C+c1*(np.dot(pc,pc.transpose())+((1-hsig)*cc*(2-cc)*C))+cmu*np.dot(artmp,np.dot(weights_diag,artmp.transpose()))

        #Adapt step size sigma
        sigma=sigma*np.exp((cs/damps)*(np.linalg.norm(ps)/chiN - 1))

        #Decomposition of C into B*diag(D.^2)*B' (diagonalization)
        if counteval-eigeneval > lambd/(c1+cmu)/N/10:
            eigeneval=counteval                    #to achieve O(N^2)
            C=np.triu(C)+np.triu(C,1).transpose()  #enforce symmetry
            B=np.linalg.eig(C)[1]                  #eigen decomposition, B==normalized eigenvectors
            diag_D=np.linalg.eig(C)[0]
            diag_D=diag_D.reshape([N,1])
            D=diag_D**0.5                       #D is a vector of standard deviations now
            D_inv=D**(-1)
            D_inv=np.diag(D_inv[:,0])
            invsqrtC=np.dot(B,np.dot(D_inv,B.transpose()))
#             print (ordered[0][1]);               # Uncomment to see convergence in the console
        #Break, if fitness is good enough or condition exceeds 1e14, better termination methods are advisable           
#         del arfitness
#         if ordered[0][1]<=stopfitness or max(D)>1e7*min(D):
#             print('bad convergence')
#             break
#         else:
#             del arx
#             del ordered
        #Return best point at last iteration


        if min(arfitness) < best_fitness:
            best_fitness = min(arfitness)
            best_conf = arx[:,best_off_indexes[0]]
    xmin=arx[:,best_off_indexes[0]]
    print(best_fitness,best_conf)
    return xmin

xmean must be a 2d array and offspring formulation was wrong. D .* randn(nGenes,1) involves Element wise multiplication.

offspring = xmean + sigma * np.dot(B,(np.multiply(D,np.random.standard_normal((N,1))))) #m + sig * Normal(0,C)