Genetic Algorithms: Solving Trade-Offs with MOGA

NUM_COMPONENTS = 4
NUM_OPTIONS = 3
 
COMPONENT_NAMES = [
    ["Plastic Frame", "Aluminum Frame", "Carbon Fiber Frame"],
    ["Slow Motor", "Fast Motor", "High-Thrust Motor"],
    ["1000mAh Battery", "2200mAh Battery", "5000mAh Battery"],
    ["1080p Camera", "4K Camera", "4K Camera w/ Gimbal"]
]

! pip install -q numpy matplotlib
 
import numpy as np
 
COSTS = np.array([
    [10, 30, 70],   # Frame: Plastic, Aluminum, Carbon Fiber
    [20, 40, 90],   # Motor: Brushless (Slow), Brushless (Fast), High-Thrust
    [15, 35, 80],   # Battery: 1000mAh, 2200mAh, 5000mAh
    [25, 60, 150]   # Camera: 1080p, 4K, 4K w/ Gimbal
])
 
ASSEMBLY_TIMES = np.array([
    [20, 15, 25],   # Frame
    [10, 15, 20],   # Motor
    [8, 12, 18],    # Battery
    [15, 20, 35]    # Camera
])
 
QUALITY_SCORES = np.array([
    [2, 5, 9],      # Frame
    [4, 6, 9],      # Motor
    [3, 6, 8],      # Battery
    [3, 7, 10]      # Camera
])

MAX_COST = 250
MIN_QUALITY = 20
 
def calculate_objectives(chromosome):
    """
    Calculates the two objectives (cost and time) for a given chromosome (drone design).
    Also checks for constraint violations.
    """
    # Chromosome is an array of indices, e.g., [0, 2, 1, 0]
    # np.arange(NUM_COMPONENTS) creates [0, 1, 2, 3] to select one option for each component
    cost = COSTS[np.arange(NUM_COMPONENTS), chromosome].sum()
    time = ASSEMBLY_TIMES[np.arange(NUM_COMPONENTS), chromosome].sum()
    quality = QUALITY_SCORES[np.arange(NUM_COMPONENTS), chromosome].sum()
 
    # Apply Penalty for constraint violation
    if cost > MAX_COST or quality < MIN_QUALITY:
        # Return a very poor fitness score for invalid solutions
        # For minimization problems, this is a very large number.
        return float('inf'), float('inf')
 
    return cost, time

def fast_non_dominated_sort(values):
    """
    Sorts the population into Pareto fronts.
    Returns a list of fronts, where each front is a list of indices.
    """
    population_size = len(values)
    domination_counts = np.zeros(population_size, dtype=int)
    dominated_solutions = [[] for _ in range(population_size)]
    fronts = [[]]
 
    for p_idx in range(population_size):
        for q_idx in range(p_idx + 1, population_size):
            p_val = values[p_idx]
            q_val = values[q_idx]
 
            # Check for domination
            if all(p_val <= q_val) and any(p_val < q_val): # p dominates q
                domination_counts[q_idx] += 1
                dominated_solutions[p_idx].append(q_idx)
            elif all(q_val <= p_val) and any(q_val < p_val): # q dominates p
                domination_counts[p_idx] += 1
                dominated_solutions[q_idx].append(p_idx)
 
    # Find the first front (all solutions with domination_count = 0)
    for i in range(population_size):
        if domination_counts[i] == 0:
            fronts[0].append(i)
 
    front_idx = 0
    while len(fronts[front_idx]) > 0:
        next_front = []
        for p_idx in fronts[front_idx]:
            for q_idx in dominated_solutions[p_idx]:
                domination_counts[q_idx] -= 1
                if domination_counts[q_idx] == 0:
                    next_front.append(q_idx)
        front_idx += 1
        fronts.append(next_front)
    
    return fronts[:-1] # Remove the last empty front

def calculate_crowding_distance(values, front):
    """
    Calculates the crowding distance for each solution in a front.
    """
    num_objectives = values.shape[1]
    front_size = len(front)
    distances = np.zeros(front_size)
 
    for m in range(num_objectives):
        # Sort the front by the current objective
        sorted_indices = sorted(range(front_size), key=lambda i: values[front[i], m])
        
        # Assign infinite distance to the boundary solutions
        distances[sorted_indices[0]] = float('inf')
        distances[sorted_indices[-1]] = float('inf')
 
        # Get the range of the objective values
        obj_min = values[front[sorted_indices[0]], m]
        obj_max = values[front[sorted_indices[-1]], m]
        
        if obj_max == obj_min:
            continue
 
        # Calculate distances for intermediate solutions
        for i in range(1, front_size - 1):
            distances[sorted_indices[i]] += (values[front[sorted_indices[i+1]], m] - values[front[sorted_indices[i-1]], m]) / (obj_max - obj_min)
            
    return distances

def create_initial_population(pop_size):
    return np.random.randint(0, NUM_OPTIONS, size=(pop_size, NUM_COMPONENTS))

def selection(population, fronts, crowding_distances):
    """Selects parents using binary tournament selection based on rank and crowding."""
    pop_size = len(population)
    selected_indices = []
    for _ in range(pop_size):
        # Pick two random individuals
        p1_idx, p2_idx = np.random.choice(pop_size, 2, replace=False)
        
        # Compare them based on Pareto rank and crowding distance
        if (fronts[p1_idx] < fronts[p2_idx]) or \
           (fronts[p1_idx] == fronts[p2_idx] and crowding_distances[p1_idx] > crowding_distances[p2_idx]):
            selected_indices.append(p1_idx)
        else:
            selected_indices.append(p2_idx)
    return population[selected_indices]

def crossover(parent1, parent2):
    point = np.random.randint(1, NUM_COMPONENTS)
    child1 = np.concatenate([parent1[:point], parent2[point:]])
    child2 = np.concatenate([parent2[:point], parent1[point:]])
    return child1, child2

def mutation(chromosome, mutation_rate):
    if np.random.rand() < mutation_rate:
        gene_to_mutate = np.random.randint(NUM_COMPONENTS)
        new_value = np.random.randint(NUM_OPTIONS)
        # Ensure the new value is different
        while new_value == chromosome[gene_to_mutate]:
            new_value = np.random.randint(NUM_OPTIONS)
        chromosome[gene_to_mutate] = new_value
    return chromosome

def run_genetic_algorithm(pop_size, generations, mutation_rate):
    """Main function to run the NSGA-II algorithm."""
    # Initialization
    population = create_initial_population(pop_size)
    initial_obj_values = np.array([calculate_objectives(ind) for ind in population])
    
    for gen in range(generations):
        # Create offspring
        offspring_population = []
        # Selection -> Crossover -> Mutation
        parents = population # For simplicity in NSGA-II, selection happens on the combined pop
        while len(offspring_population) < pop_size:
            p1, p2 = parents[np.random.choice(len(parents), 2)]
            c1, c2 = crossover(p1, p2)
            offspring_population.append(mutation(c1, mutation_rate))
            offspring_population.append(mutation(c2, mutation_rate))
        
        offspring_population = np.array(offspring_population)[:pop_size]
 
        # Combine parent and offspring populations
        combined_population = np.vstack([population, offspring_population])
        obj_values = np.array([calculate_objectives(ind) for ind in combined_population])
        
        # Filter out invalid solutions before sorting
        valid_indices = np.where(obj_values[:, 0] != float('inf'))[0]
        if len(valid_indices) == 0:
            print(f"Generation {gen+1}: No valid solutions found, re-initializing population. \n")
            population = create_initial_population(pop_size)
            continue
            
        combined_population = combined_population[valid_indices]
        obj_values = obj_values[valid_indices]
 
        # Perform NSGA-II sorting
        fronts = fast_non_dominated_sort(obj_values)
        
        # Create the next generation
        new_population = []
        front_num = 0
        while len(new_population) + len(fronts[front_num]) <= pop_size:
            new_population.extend(combined_population[fronts[front_num]])
            front_num += 1
 
        # Fill remaining spots with the most spread-out solutions from the next front
        if len(new_population) < pop_size:
            crowding_distances = calculate_crowding_distance(obj_values, fronts[front_num])
            sorted_by_crowding = sorted(zip(fronts[front_num], crowding_distances), key=lambda x: x[1], reverse=True)
            
            remaining = pop_size - len(new_population)
            for i in range(remaining):
                new_population.append(combined_population[sorted_by_crowding[i][0]])
        
        population = np.array(new_population)
        if (gen + 1) % 50 == 0:
            print(f"Generation {gen+1}/{generations} completed. \n")
            
    # Final population results
    final_obj_values = np.array([calculate_objectives(ind) for ind in population])
    final_fronts = fast_non_dominated_sort(final_obj_values)
    
    return population, final_obj_values, final_fronts, initial_obj_values

POPULATION_SIZE = 200
GENERATIONS = 300
MUTATION_RATE = 0.1
 
# --- Run and Visualize ---
final_pop, final_values, final_fronts, initial_values = run_genetic_algorithm(
    pop_size=POPULATION_SIZE,
    generations=GENERATIONS,
    mutation_rate=MUTATION_RATE
)

import matplotlib.pyplot as plt
import matplotlib.colors as mcolors
 
def plot_results(initial_obj, final_obj, fronts):
    """Plots the initial population and the final Pareto front."""
    plt.style.use('seaborn-v0_8-whitegrid')
    fig, ax = plt.subplots(figsize=(12, 8))
 
    # Plot initial random population
    valid_initial = initial_obj[initial_obj[:, 0] != float('inf')]
    ax.scatter(valid_initial[:, 0], valid_initial[:, 1], c='gray', alpha=0.5, label='Initial Random Designs')
 
    # Plot the final population, with colors for different fronts
    colors = list(mcolors.TABLEAU_COLORS.values())
    
    for i, front in enumerate(fronts):
        if not front: continue
        front_values = final_obj[front]
        label = f'Pareto Front (Rank {i+1})' if i == 0 else f'Front {i+1}'
        zorder = len(fronts) - i # Draw best front on top
        ax.scatter(front_values[:, 0], front_values[:, 1], c=colors[i % len(colors)], label=label, s=80, zorder=zorder, edgecolors='k', alpha=0.8)
 
    ax.set_title('Drone Manufacturing: Cost vs. Assembly Time Optimization', fontsize=18)
    ax.set_xlabel('Total Cost ($)', fontsize=14)
    ax.set_ylabel('Total Assembly Time (minutes)', fontsize=14)
    ax.legend()
    ax.grid(True)
    plt.tight_layout()
    plt.show()

plot_results(initial_values, final_values, final_fronts)

print("\n" + "="*50 + "\n")
print("      Optimal Drone Designs on the Pareto Front\n")
print("="*50 + "\n")
 
front_indices = final_fronts[0]
front_solutions = final_pop[front_indices]
front_values = final_values[front_indices]
 
sorted_indices = np.argsort(front_values[:, 0])
 
for i, idx in enumerate(sorted_indices):
    solution_chromosome = front_solutions[idx]
    cost, time = front_values[idx]
    
    print(f"\n--- Design #{i+1} (The Trade-Off Choice) --- \n")
    print(f"    - Total Cost: ${cost:.2f}\n")
    print(f"    - Assembly Time: {int(time)} minutes\n")
    print("    - Components:")
    print(f"        1. Frame:   {COMPONENT_NAMES[0][solution_chromosome[0]]}\n")
    print(f"        2. Motor:   {COMPONENT_NAMES[1][solution_chromosome[1]]}\n")
    print(f"        3. Battery: {COMPONENT_NAMES[2][solution_chromosome[2]]}\n")
    print(f"        4. Camera:  {COMPONENT_NAMES[3][solution_chromosome[3]]}\n")
    break