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

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])

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)    #minimization, list with (Index,Fitness) elements
xold=xmean
best_off_indexes=[]
for i in range(int(mu)):
best_off_indexes.append(ordered[i]) #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

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()))

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)                  #eigen decomposition, B==normalized eigenvectors
diag_D=np.linalg.eig(C)
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;               # 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<=stopfitness or max(D)>1e7*min(D):
break
else:
del arx
del ordered
#Return best point at last iteration
xmin=arx[:,best_off_indexes]
return xmin

• Could you provide test code that invokes the above in a way that you expect reproduces the memory issue? – Reinderien Dec 20 '14 at 4:32

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.
import numpy as np
import math
def frosenbrock(x):

x1 = x
x2 = x
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])

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)    #minimization, list with (Index,Fitness) elements
xold=xmean
best_off_indexes=[]
for i in range(int(mu)):
best_off_indexes.append(ordered[i]) #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

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()))

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)                  #eigen decomposition, B==normalized eigenvectors
diag_D=np.linalg.eig(C)
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);               # 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<=stopfitness or max(D)>1e7*min(D):
#             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]
xmin=arx[:,best_off_indexes]
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)