Solution: 

The source code in python would look like below,

import random
import numpy as np

# Distance matrix representing distances between delivery points
distance_matrix = np.array([
   [0, 29, 20, 21],
   [29, 0, 15, 17],
   [20, 15, 0, 28],
   [21, 17, 28, 0]
])

# Number of delivery points
num_points = len(distance_matrix)

# Parameters for GA
population_size = 10
generations = 100
mutation_rate = 0.1

def create_route():
   """ Create a random route """
   route = list(range(num_points))
   random.shuffle(route)
   return route

def initial_population(pop_size):
   """ Generate initial population """
   return [create_route() for _ in range(pop_size)]

def route_distance(route):
   """ Calculate the total distance of the route """
   distance = 0
   for i in range(len(route)):
      from_city = route[i]
      to_city = route[(i + 1) % num_points]
      distance += distance_matrix[from_city][to_city]
   return distance

def fitness(route):
   """ Fitness function """
   return 1 / route_distance(route)

def selection(population, fitness_scores):
   """ Select pairs of individuals based on fitness """
   fitness_total = sum(fitness_scores)
   selection_probs = [f/fitness_total for f in fitness_scores]
   selected = np.random.choice(population, size=len(population), p=selection_probs)
   return selected

def crossover(parent1, parent2):
   """ Perform order crossover """
   start, end = sorted(random.sample(range(num_points), 2))
   child = [-1] * num_points
   child[start:end] = parent1[start:end]

   fill_values = [item for item in parent2 if item not in child]
   fill_idx = [i for i, x in enumerate(child) if x == -1]

   for i, idx in enumerate(fill_idx):
      child[idx] = fill_values[i]

   return child

def mutate(route):
   """ Perform swap mutation """
   if random.random() < mutation_rate:
      i, j = random.sample(range(num_points), 2)
      route[i], route[j] = route[j], route[i]
   return route

def evolve_population(population):
   """ Evolve population by selection, crossover and mutation """
   fitness_scores = [fitness(route) for route in population]
   selected_population = selection(population, fitness_scores)

   new_population = []
   for i in range(0, len(selected_population), 2):
      parent1, parent2 = selected_population[i], selected_population[i+1]
      child1 = crossover(parent1, parent2)
      child2 = crossover(parent2, parent1)
      new_population.extend([mutate(child1), mutate(child2)])

   return new_population

# GA process
population = initial_population(population_size)
for generation in range(generations):
   population = evolve_population(population)

# Find the best route in the final population
best_route = min(population, key=route_distance)
best_distance = route_distance(best_route)

print("Optimized Delivery Route:", best_route)
print("Total Distance:", best_distance)