This code is provided as supplementary material of the lecture Machine Learning and Optimization in Communications (MLOC).
This code illustrates:
In [1]:
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import ticker
For illustration, we optimize Rosenbrock's function in 2 dimensions. The function is given by \begin{equation*} f(\boldsymbol{x}) = (1-x_2)^2 + 100(x_2-x_1^2)^2 \end{equation*} such that $x_1 \geq 0$ and $x_2 \geq 0$.
In [2]:
# Rosenbrocks function
def rosenbrock(x,y):
return (1-x)**2+100*(y-(x**2))**2
def fun(x,y):
value = rosenbrock(x,y)
# CONSTRAINT: x and y should be non-negative. Hence, add a penalty for negative numbers that becomes larger the more negative the number is
if x < 0:
value += 100*(-x)
if y < 0:
value += 100*(-y)
return value
# vector-version of the function
vrosenbrock = np.vectorize(rosenbrock)
vfun = np.vectorize(fun)
Plot the function as a heat map.
In [3]:
# plot map of Griewanks function
x = np.arange(-2.0, 5.1, 0.1)
y = np.arange(-2.0, 5.1, 0.1)
X, Y = np.meshgrid(x, y)
fZ = vrosenbrock(X,Y)
plt.figure(1,figsize=(10,9))
plt.rcParams.update({'font.size': 14})
temp = np.log(fZ)
temp[temp < -1] = 0
plt.contourf(X,Y,temp,levels=20)
plt.fill_betweenx([-2, 5],[-2,-2],[0,0], alpha=0.5, color='white',linewidth=0)
plt.fill_betweenx([-2, 0],[0,0],[5,5], alpha=0.5, color='white',linewidth=0)
#plt.axhspan(-5, 0, alpha=0.5, color='white')
plt.axis('scaled')
plt.xlabel("$x_1$")
plt.ylabel("$x_2$")
plt.show()
Helper function to generate a random initial population of $D$ dimensions $N_P$ elements. The elementsof the population are randomly distributed in the interval $[x_{\min},x_{\max}]$ (in every dimension).
In [4]:
def initial_population(D, NP, xmin, xmax):
v = np.random.rand(NP, D)*(xmax - xmin) + xmin
return v
Carry out differential evolution similar (not identical, slightly modified) to the scheme DE1 described in [1]
[1] R. Storn and K. Price, "Differential Evolution - A simple and efficient adaptive scheme for global optimization over continuous spaces", Technical Report TR-95-012, March 1995
In [5]:
#dimension
D = 2
# population
NP = 15*D
# twiddling parameter
F = 0.8
# cross-over probability
CR = 0.3
# maximum 1000 iterations
max_iter = 1000
# generate initial population
population = initial_population(D, NP, -2, 5)[:]
# compute initial cost
cost = vfun(population[:,0], population[:,1])
best_index = np.argmin(cost)
best_cost = cost[best_index]
iteration = 0
# keep track of population
save_population = []
while iteration < max_iter:
# loop over every element from the population
for k in range(NP):
# get 4 random elements
rp = np.random.permutation(NP)[0:4]
# remove ourselves from the list
rp = [j for j in rp if j != k]
# generate new candidate vector
v = population[rp[0],:] + F*( population[rp[1],:] - population[rp[2],:] )
# take vector from population
u = np.array(population[k,:])
# cross-over each coordinate with probability CR with entry from candidate vector v
idx = np.random.rand(D) < CR
# cross-over
u[idx] = v[idx]
new_cost = fun(u[0], u[1])
if new_cost < cost[k]:
# better cost? keep!
cost[k] = new_cost
population[k,:] = u
if new_cost < best_cost:
best_cost = new_cost
best_index = k
save_population.append(np.array(population[:]))
iteration += 1
if iteration % 100 == 0:
print('After iteration %d, best cost %1.4f (obtained for (%1.2f,%1.2f))' % (iteration, best_cost, population[best_index,0], population[best_index,1]))
In [7]:
plt.figure(1,figsize=(9,9))
plt.rcParams.update({'font.size': 14})
plt.contourf(X,Y,temp,levels=20)
index = 180
cost = vfun(save_population[index][:,0], save_population[index][:,1])
best_index = np.argmin(cost)
plt.fill_betweenx([-2, 5],[-2,-2],[0,0], alpha=0.5, color='white',linewidth=0)
plt.fill_betweenx([-2, 0],[0,0],[5,5], alpha=0.5, color='white',linewidth=0)
plt.scatter(save_population[index][:,0], save_population[index][:,1], c='w')
plt.scatter(save_population[index][best_index,0], save_population[index][best_index,1], c='r')
plt.xlim((-2,5))
plt.ylim((-2,5))
plt.xlabel("$x_1$")
plt.ylabel("$x_2$")
plt.savefig('DE_Rosenbrock.pdf',bbox_inches='tight')
Generate animation.
In [49]:
%matplotlib notebook
# Generate animation
from matplotlib import animation, rc
from matplotlib.animation import PillowWriter # Disable if you don't want to save any GIFs.
font = {'size' : 18}
plt.rc('font', **font)
fig, ax = plt.subplots(1, figsize=(10,10))
ax.set_xlim(( -2, 5))
ax.set_ylim(( -2, 5))
ax.axis('scaled')
written = False
def animate(i):
ax.clear()
ax.contourf(X,Y,temp,levels=20)
cost = vfun(save_population[i][:,0], save_population[i][:,1])
best_index = np.argmin(cost)
ax.fill_betweenx([-2, 5],[-2,-2],[0,0], alpha=0.5, color='white',linewidth=0)
ax.fill_betweenx([-2, 0],[0,0],[5,5], alpha=0.5, color='white',linewidth=0)
ax.scatter(save_population[i][:,0], save_population[i][:,1], c='w')
ax.scatter(save_population[i][best_index,0], save_population[i][best_index,1], c='r')
ax.set_xlabel(r'$x_1$',fontsize=18)
ax.set_ylabel(r'$x_2$',fontsize=18)
ax.set_xlim(( -2, 5))
ax.set_ylim(( -2, 5))
anim = animation.FuncAnimation(fig, animate, frames=300, interval=80, blit=False)
fig.show()
anim.save('differential_evolution_Rosenbrock.gif', writer=PillowWriter(fps=7))