Hybrid Non-Dominated Sorting Genetic Algorithm (Hybrid NSGA-II) in Python

import random as rn
import numpy as np
import matplotlib.pyplot as plt
import math

#_____________________________________________________________________________
def random_population(nv,n,lb,ub):
    # nv = numver of variables
    # n = number of random solutions
    # lb = lower bound
    # ub = upper bound
    pop=np.zeros((n, nv)) 
    for i in range(n):
        pop[i,:] = np.random.uniform(lb,ub)

    return pop
#_____________________________________________________________________________
def crossover(pop, crossover_rate):
    offspring = np.zeros((crossover_rate, pop.shape[1]))
    for i in range(int(crossover_rate/2)):
        r1=np.random.randint(0, pop.shape[0])
        r2 = np.random.randint(0, pop.shape[0])
        while r1 == r2:
            r1 = np.random.randint(0, pop.shape[0])
            r2 = np.random.randint(0, pop.shape[0])
        cutting_point = np.random.randint(1, pop.shape[1])
        offspring[2*i, 0:cutting_point] = pop[r1, 0:cutting_point]
        offspring[2*i, cutting_point:] = pop[r2, cutting_point:]
        offspring[2*i+1, 0:cutting_point] = pop[r2, 0:cutting_point]
        offspring[2*i+1, cutting_point:] = pop[r1, cutting_point:]

    return offspring
#_____________________________________________________________________________
def mutation(pop, mutation_rate):
    offspring = np.zeros((mutation_rate, pop.shape[1]))
    for i in range(int(mutation_rate/2)):
        r1=np.random.randint(0, pop.shape[0])
        r2 = np.random.randint(0, pop.shape[0])
        while r1 == r2:
            r1 = np.random.randint(0, pop.shape[0])
            r2 = np.random.randint(0, pop.shape[0])
        cutting_point = np.random.randint(0, pop.shape[1])
        offspring[2*i] = pop[r1]
        offspring[2*i,cutting_point] = pop[r2,cutting_point]
        offspring[2*i+1] = pop[r2]
        offspring[2*i+1, cutting_point] = pop[r1, cutting_point]

    return offspring
#_____________________________________________________________________________
def local_search(pop, n, step_size):
    # number of offspring chromosomes generated from the local search
    offspring = np.zeros((n, pop.shape[1]))
    for i in range(n):
    ...    
          offspring[i,:] = chromosome
    return offspring
#_____________________________________________________________________________
def evaluation(pop):
    fitness_values = np.zeros((pop.shape[0], 2)) # because of 2 objective functions
    for i, x in enumerate(pop):
     
        obj1 = 0
        for j in range(2):
            obj1 += - 10*math.exp(-0.2*math.sqrt((x[j])**2 + (x[j+1])**2))

        obj2 = 0
        for j in range(3):
            obj2 += (abs(x[j]))**0.8 + 5*math.sin((x[j])**3)

        fitness_values[i,0] = obj1 # objective 1
        fitness_values[i,1] = obj2 # objective 2

    return fitness_values
#_____________________________________________________________________________
def crowding_calculation(fitness_values):
    
    pop_size = len(fitness_values[:, 0])
    fitness_value_number = len(fitness_values[0, :])
    matrix_for_crowding = np.zeros((pop_size, fitness_value_number))
    ...
    for i in range(fitness_value_number):
        crowding_results = np.zeros(pop_size)
        ...
    return crowding_distance
#_____________________________________________________________________________
def remove_using_crowding(fitness_values, number_solutions_needed):
  
    pop_index = np.arange(fitness_values.shape[0])
    crowding_distance = crowding_calculation(fitness_values)
    selected_pop_index = np.zeros((number_solutions_needed))
    ...
    
    for i in range(number_solutions_needed):
        pop_size = pop_index.shape[0]
        solution_1 = rn.randint(0, pop_size - 1)
        ...
    
    return (selected_pop_index)
#_____________________________________________________________________________
def pareto_front_finding(fitness_values, pop_index):

    pop_size = fitness_values.shape[0]
    ...
    for i in range(pop_size):
        for j in range(pop_size):
    ...
    
    return pop_index[pareto_front]
#_____________________________________________________________________________
def selection(pop, fitness_values, pop_size):
    
    pop_index_0 = np.arange(pop.shape[0])
    pop_index = np.arange(pop.shape[0])
    pareto_front_index = []
    ...
   
    selected_pop = pop[pareto_front_index.astype(int)]

    return selected_pop
#_____________________________________________________________________________
# Parameters
nv = 3
lb = [-5, -5, -5]
ub = [5, 5, 5]
pop_size = 150
rate_crossover = 20
rate_mutation = 20
rate_local_search = 10
step_size = 0.1
pop = random_population(nv,pop_size,lb,ub)
#_____________________________________________________________________________
# Main loop of NSGA II

for i in range(200):
    offspring_from_crossover = crossover(pop,rate_crossover)
    offspring_from_mutation = mutation(pop,rate_mutation)
    ...
    fitness_values = evaluation(pop)
    pop = selection(pop,fitness_values, pop_size)
    print('iteration:', i)
#_____________________________________________________________________________
# Pareto front visualization
fitness_values = evaluation(pop)
index = np.arange(pop.shape[0]).astype(int)
pareto_front_index = pareto_front_finding(fitness_values, index)
pop = pop[pareto_front_index, :]
print("_________________")
print("Optimal solutions:")
print("       x1               x2                 x3")
print(pop) # show optimal solutions
fitness_values = fitness_values[pareto_front_index]
print("______________")
print("Fitness values:")
print("  objective 1    objective 2")
print("      |              |")
print(fitness_values)
plt.scatter(fitness_values[:, 0],fitness_values[:, 1])
plt.xlabel('Objective function 1')
plt.ylabel('Objective function 2')
plt.show()
.
.
.
Sorry! This is only a half of the Python code.