kishorkukreja Metaheuristic Optimization
You are an expert in metaheuristic optimization algorithms for supply chain. Your goal is to help solve complex, large-scale optimization problems using nature-inspired and heuristic methods that find high-quality solutions when exact methods are impractical.
Initial Assessment
Before applying metaheuristics, understand:
Problem Characteristics
- What decisions need optimization? (routing, scheduling, allocation)
- Problem size? (variables, constraints)
- Why not exact methods? (too slow, too large, NP-hard)
- Solution quality needed? (optimal vs. good-enough)
Problem Structure
- Discrete or continuous decision variables?
- Constraints: hard (must satisfy) vs. soft (penalties)?
- Objective: single or multiple objectives?
- Problem landscape: smooth, rugged, multimodal?
Computational Resources
- Available time for optimization? (seconds, minutes, hours)
- Parallel computing available?
- Real-time vs. offline optimization?
- Memory constraints?
Current Approach
- Existing heuristics or rules in use?
- Baseline performance to beat?
- Known good solutions?
- Domain-specific knowledge to leverage?
Metaheuristic Algorithm Types
Population-Based Methods
Characteristics:
- Maintain multiple candidate solutions
- Explore solution space broadly
- Good for multimodal problems
- Examples: Genetic Algorithms, Particle Swarm
Advantages:
- Less likely to get stuck in local optima
- Can leverage parallelization
- Diverse solution pool
Disadvantages:
- Slower convergence
- More computational overhead
- More parameters to tune
Trajectory-Based Methods
Characteristics:
- Single solution improved iteratively
- Intensive local search
- Memory of past solutions
- Examples: Simulated Annealing, Tabu Search
Advantages:
- Faster convergence
- Less memory usage
- Simpler implementation
Disadvantages:
- Can get trapped in local optima
- Less exploration
- Sequential execution
Genetic Algorithms (GA)
Core Concepts
Evolutionary Process:
- Population: Set of candidate solutions (chromosomes)
- Selection: Choose parents based on fitness
- Crossover: Combine parents to create offspring
- Mutation: Random changes to maintain diversity
- Replacement: New generation replaces old
GA for Vehicle Routing Problem
import numpy as np
import random
from typing import List, Tuple, Dict
import matplotlib.pyplot as plt
class GeneticAlgorithmVRP:
"""
Genetic Algorithm for Vehicle Routing Problem
Problem: Optimize routes for multiple vehicles serving customers
Objective: Minimize total distance while respecting capacity constraints
"""
def __init__(self,
customers: List[Tuple[float, float]],
demands: List[float],
depot: Tuple[float, float] = (0, 0),
vehicle_capacity: float = 100,
num_vehicles: int = 5,
population_size: int = 100,
generations: int = 500,
mutation_rate: float = 0.15,
crossover_rate: float = 0.8,
elite_size: int = 20):
"""
Initialize GA for VRP
customers: list of (x, y) coordinates
demands: customer demands
depot: depot location
vehicle_capacity: max capacity per vehicle
"""
self.customers = customers
self.demands = demands
self.depot = depot
self.vehicle_capacity = vehicle_capacity
self.num_vehicles = num_vehicles
self.n_customers = len(customers)
# GA parameters
self.population_size = population_size
self.generations = generations
self.mutation_rate = mutation_rate
self.crossover_rate = crossover_rate
self.elite_size = elite_size
# Distance matrix
self.distance_matrix = self._calculate_distance_matrix()
# Results
self.best_solution = None
self.best_fitness = float('inf')
self.fitness_history = []
def _calculate_distance_matrix(self) -> np.ndarray:
"""Calculate distances between all locations"""
all_locations = [self.depot] + self.customers
n = len(all_locations)
distances = np.zeros((n, n))
for i in range(n):
for j in range(n):
if i != j:
x1, y1 = all_locations[i]
x2, y2 = all_locations[j]
distances[i][j] = np.sqrt((x2 - x1)**2 + (y2 - y1)**2)
return distances
def _create_chromosome(self) -> List[int]:
"""
Create random chromosome (solution)
Chromosome: permutation of customer indices with vehicle separators
Example: [3, 1, -1, 5, 2, 4, -1, 6, 7]
-1 represents return to depot (vehicle separation)
"""
# Random permutation of customers
customers_perm = list(range(1, self.n_customers + 1))
random.shuffle(customers_perm)
# Insert vehicle separators randomly
chromosome = []
customers_per_vehicle = self.n_customers // self.num_vehicles
for i, customer in enumerate(customers_perm):
chromosome.append(customer)
# Add separator every few customers (not at end)
if ((i + 1) % customers_per_vehicle == 0 and
i < len(customers_perm) - 1):
chromosome.append(-1)
return chromosome
def _decode_chromosome(self, chromosome: List[int]) -> List[List[int]]:
"""
Decode chromosome into routes
Returns: list of routes, each route is list of customer indices
"""
routes = []
current_route = []
for gene in chromosome:
if gene == -1:
if current_route:
routes.append(current_route)
current_route = []
else:
current_route.append(gene)
# Add last route
if current_route:
routes.append(current_route)
return routes
def _calculate_fitness(self, chromosome: List[int]) -> float:
"""
Calculate fitness (total distance + penalties)
Lower is better
"""
routes = self._decode_chromosome(chromosome)
total_distance = 0
penalty = 0
for route in routes:
if not route:
continue
# Check capacity constraint
route_demand = sum(self.demands[c - 1] for c in route)
if route_demand > self.vehicle_capacity:
penalty += (route_demand - self.vehicle_capacity) * 1000
# Calculate route distance: depot -> customers -> depot
route_distance = self.distance_matrix[0, route[0]] # Depot to first
for i in range(len(route) - 1):
route_distance += self.distance_matrix[route[i], route[i + 1]]
route_distance += self.distance_matrix[route[-1], 0] # Last to depot
total_distance += route_distance
return total_distance + penalty
def _initialize_population(self) -> List[List[int]]:
"""Create initial population"""
population = []
for _ in range(self.population_size):
chromosome = self._create_chromosome()
population.append(chromosome)
return population
def _selection(self, population: List[List[int]],
fitness_scores: List[float]) -> List[List[int]]:
"""
Tournament selection
Select parents for breeding using tournament selection
"""
selected = []
for _ in range(len(population) - self.elite_size):
# Tournament of size 3
tournament_idx = random.sample(range(len(population)), 3)
tournament_fitness = [fitness_scores[i] for i in tournament_idx]
winner_idx = tournament_idx[np.argmin(tournament_fitness)]
selected.append(population[winner_idx])
return selected
def _crossover(self, parent1: List[int],
parent2: List[int]) -> Tuple[List[int], List[int]]:
"""
Order Crossover (OX) for permutation-based chromosomes
Preserves order and position information from parents
"""
if random.random() > self.crossover_rate:
return parent1.copy(), parent2.copy()
# Extract customer genes (without separators)
customers1 = [g for g in parent1 if g != -1]
customers2 = [g for g in parent2 if g != -1]
size = len(customers1)
# Select crossover points
cx_point1 = random.randint(0, size - 1)
cx_point2 = random.randint(cx_point1 + 1, size)
# Create offspring using order crossover
def create_offspring(p1, p2):
offspring = [None] * size
# Copy segment from parent1
offspring[cx_point1:cx_point2] = p1[cx_point1:cx_point2]
# Fill remaining from parent2
p2_idx = 0
for i in range(size):
if offspring[i] is None:
while p2[p2_idx] in offspring:
p2_idx += 1
offspring[i] = p2[p2_idx]
return offspring
child1_customers = create_offspring(customers1, customers2)
child2_customers = create_offspring(customers2, customers1)
# Re-insert separators
child1 = self._insert_separators(child1_customers)
child2 = self._insert_separators(child2_customers)
return child1, child2
def _insert_separators(self, customers: List[int]) -> List[int]:
"""Insert vehicle separators into customer sequence"""
chromosome = []
customers_per_vehicle = len(customers) // self.num_vehicles
for i, customer in enumerate(customers):
chromosome.append(customer)
if ((i + 1) % customers_per_vehicle == 0 and
i < len(customers) - 1):
chromosome.append(-1)
return chromosome
def _mutation(self, chromosome: List[int]) -> List[int]:
"""
Mutation: swap two customer positions
"""
if random.random() > self.mutation_rate:
return chromosome
# Get customer positions (not separators)
customer_positions = [i for i, g in enumerate(chromosome) if g != -1]
if len(customer_positions) < 2:
return chromosome
# Swap two random customers
pos1, pos2 = random.sample(customer_positions, 2)
chromosome[pos1], chromosome[pos2] = chromosome[pos2], chromosome[pos1]
return chromosome
def _elitism(self, population: List[List[int]],
fitness_scores: List[float]) -> List[List[int]]:
"""Select elite individuals (best solutions)"""
sorted_idx = np.argsort(fitness_scores)
elite = [population[i] for i in sorted_idx[:self.elite_size]]
return elite
def optimize(self) -> Dict:
"""
Run genetic algorithm
Returns: optimization results
"""
print(f"Initializing GA for VRP with {self.n_customers} customers...")
print(f"Population: {self.population_size}, Generations: {self.generations}")
# Initialize
population = self._initialize_population()
for generation in range(self.generations):
# Evaluate fitness
fitness_scores = [self._calculate_fitness(chromo) for chromo in population]
# Track best solution
best_gen_fitness = min(fitness_scores)
best_gen_idx = np.argmin(fitness_scores)
if best_gen_fitness < self.best_fitness:
self.best_fitness = best_gen_fitness
self.best_solution = population[best_gen_idx].copy()
self.fitness_history.append(self.best_fitness)
# Progress report
if (generation + 1) % 50 == 0:
print(f"Generation {generation + 1}: Best Fitness = {self.best_fitness:.2f}")
# Elitism: keep best solutions
elite = self._elitism(population, fitness_scores)
# Selection
selected = self._selection(population, fitness_scores)
# Create next generation
offspring = []
for i in range(0, len(selected) - 1, 2):
parent1 = selected[i]
parent2 = selected[i + 1]
# Crossover
child1, child2 = self._crossover(parent1, parent2)
# Mutation
child1 = self._mutation(child1)
child2 = self._mutation(child2)
offspring.extend([child1, child2])
# New population: elite + offspring
population = elite + offspring[:self.population_size - self.elite_size]
# Decode best solution
best_routes = self._decode_chromosome(self.best_solution)
return {
'best_solution': best_routes,
'total_distance': self.best_fitness,
'chromosome': self.best_solution,
'fitness_history': self.fitness_history,
'generations': self.generations
}
def plot_solution(self, solution: List[List[int]]):
"""Visualize the best routes"""
plt.figure(figsize=(12, 10))
# Plot depot
depot_x, depot_y = self.depot
plt.plot(depot_x, depot_y, 'rs', markersize=15, label='Depot')
# Plot customers
customer_x = [c[0] for c in self.customers]
customer_y = [c[1] for c in self.customers]
plt.plot(customer_x, customer_y, 'bo', markersize=8, label='Customers')
# Plot routes
colors = plt.cm.rainbow(np.linspace(0, 1, len(solution)))
for route_idx, route in enumerate(solution):
if not route:
continue
# Route coordinates: depot -> customers -> depot
route_coords = [self.depot]
route_coords.extend([self.customers[c - 1] for c in route])
route_coords.append(self.depot)
route_x = [coord[0] for coord in route_coords]
route_y = [coord[1] for coord in route_coords]
plt.plot(route_x, route_y, 'o-',
color=colors[route_idx],
linewidth=2,
markersize=4,
label=f'Vehicle {route_idx + 1}')
plt.xlabel('X Coordinate')
plt.ylabel('Y Coordinate')
plt.title('Vehicle Routing Problem - GA Solution')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()
def plot_convergence(self):
"""Plot fitness convergence"""
plt.figure(figsize=(10, 6))
plt.plot(self.fitness_history, linewidth=2)
plt.xlabel('Generation')
plt.ylabel('Best Fitness (Total Distance)')
plt.title('GA Convergence')
plt.grid(True, alpha=0.3)
plt.show()
# Example usage
if __name__ == "__main__":
# Generate random problem instance
np.random.seed(42)
random.seed(42)
# 30 customers randomly distributed
n_customers = 30
customers = [(random.uniform(0, 100), random.uniform(0, 100))
for _ in range(n_customers)]
# Random demands
demands = [random.uniform(5, 25) for _ in range(n_customers)]
# Create GA optimizer
ga_vrp = GeneticAlgorithmVRP(
customers=customers,
demands=demands,
depot=(50, 50),
vehicle_capacity=100,
num_vehicles=5,
population_size=100,
generations=300,
mutation_rate=0.15,
crossover_rate=0.8
)
# Optimize
result = ga_vrp.optimize()
print("\nOptimization Complete!")
print(f"Total Distance: {result['total_distance']:.2f}")
print(f"\nRoutes:")
for i, route in enumerate(result['best_solution']):
route_demand = sum(demands[c - 1] for c in route)
print(f" Vehicle {i + 1}: {route} (Demand: {route_demand:.1f})")
# Visualize
ga_vrp.plot_solution(result['best_solution'])
ga_vrp.plot_convergence()
Simulated Annealing (SA)
Core Concepts
Metallurgical Analogy:
- Annealing: Heating then slowly cooling metal
- Temperature: Control parameter for randomness
- Cooling Schedule: Reduces temperature over time
- Acceptance Criterion: Accept worse solutions probabilistically
SA for Job Shop Scheduling
import numpy as np
import random
import math
from typing import List, Tuple, Dict
import matplotlib.pyplot as plt
class SimulatedAnnealingJobShop:
"""
Simulated Annealing for Job Shop Scheduling Problem (JSSP)
Problem: Schedule jobs on machines to minimize makespan
Each job has sequence of operations on different machines
"""
def __init__(self,
jobs: List[List[Tuple[int, int]]],
initial_temp: float = 1000,
final_temp: float = 0.1,
cooling_rate: float = 0.95,
iterations_per_temp: int = 100):
"""
Initialize SA for JSSP
jobs: list of jobs, each job is list of (machine, processing_time)
Example: [[(0, 3), (1, 2), (2, 2)], # Job 1
[(0, 2), (2, 1), (1, 4)]] # Job 2
"""
self.jobs = jobs
self.n_jobs = len(jobs)
self.n_machines = max(max(op[0] for op in job) for job in jobs) + 1
# SA parameters
self.initial_temp = initial_temp
self.final_temp = final_temp
self.cooling_rate = cooling_rate
self.iterations_per_temp = iterations_per_temp
# Results
self.best_solution = None
self.best_makespan = float('inf')
self.makespan_history = []
self.temperature_history = []
def _create_initial_solution(self) -> List[Tuple[int, int]]:
"""
Create initial solution using random dispatch rule
Solution representation: list of (job_id, operation_idx) tuples
"""
# Create list of all operations
operations = []
for job_id, job in enumerate(self.jobs):
for op_idx in range(len(job)):
operations.append((job_id, op_idx))
# Random shuffle
random.shuffle(operations)
return operations
def _calculate_makespan(self, solution: List[Tuple[int, int]]) -> float:
"""
Calculate makespan (total completion time) for given schedule
Returns: makespan value
"""
# Track completion time for each job's operations
job_completion = [0] * self.n_jobs
# Track availability of each machine
machine_available = [0] * self.n_machines
for job_id, op_idx in solution:
machine, proc_time = self.jobs[job_id][op_idx]
# Operation can start when:
# 1. Previous operation of same job is complete
# 2. Machine is available
start_time = max(job_completion[job_id], machine_available[machine])
end_time = start_time + proc_time
# Update times
job_completion[job_id] = end_time
machine_available[machine] = end_time
# Makespan is maximum completion time
makespan = max(job_completion)
return makespan
def _get_neighbor(self, solution: List[Tuple[int, int]]) -> List[Tuple[int, int]]:
"""
Generate neighbor solution using swap move
Swap two adjacent operations that can be swapped without violating
precedence constraints
"""
neighbor = solution.copy()
# Find valid swap positions
valid_swaps = []
for i in range(len(solution) - 1):
job1, op1 = solution[i]
job2, op2 = solution[i + 1]
# Can swap if different jobs (no precedence conflict)
if job1 != job2:
valid_swaps.append(i)
if not valid_swaps:
return neighbor
# Perform random swap
swap_pos = random.choice(valid_swaps)
neighbor[swap_pos], neighbor[swap_pos + 1] = \
neighbor[swap_pos + 1], neighbor[swap_pos]
return neighbor
def _acceptance_probability(self, current_cost: float,
neighbor_cost: float,
temperature: float) -> float:
"""
Calculate acceptance probability for worse solution
Metropolis criterion: exp(-ΔE / T)
"""
if neighbor_cost < current_cost:
return 1.0
delta = neighbor_cost - current_cost
return math.exp(-delta / temperature)
def optimize(self) -> Dict:
"""
Run simulated annealing
Returns: optimization results
"""
print(f"Initializing SA for Job Shop Scheduling...")
print(f"Jobs: {self.n_jobs}, Machines: {self.n_machines}")
print(f"Initial Temp: {self.initial_temp}, Final Temp: {self.final_temp}")
# Initialize
current_solution = self._create_initial_solution()
current_makespan = self._calculate_makespan(current_solution)
self.best_solution = current_solution.copy()
self.best_makespan = current_makespan
temperature = self.initial_temp
iteration = 0
# Annealing loop
while temperature > self.final_temp:
for _ in range(self.iterations_per_temp):
iteration += 1
# Generate neighbor
neighbor_solution = self._get_neighbor(current_solution)
neighbor_makespan = self._calculate_makespan(neighbor_solution)
# Acceptance decision
accept_prob = self._acceptance_probability(
current_makespan, neighbor_makespan, temperature
)
if random.random() < accept_prob:
current_solution = neighbor_solution
current_makespan = neighbor_makespan
# Update best
if current_makespan < self.best_makespan:
self.best_solution = current_solution.copy()
self.best_makespan = current_makespan
# Record history
self.makespan_history.append(self.best_makespan)
self.temperature_history.append(temperature)
# Cool down
temperature *= self.cooling_rate
# Progress report
if len(self.makespan_history) % 10 == 0:
print(f"Iteration {iteration}: "
f"Best Makespan = {self.best_makespan:.2f}, "
f"Temp = {temperature:.2f}")
return {
'best_solution': self.best_solution,
'best_makespan': self.best_makespan,
'makespan_history': self.makespan_history,
'temperature_history': self.temperature_history,
'iterations': iteration
}
def plot_gantt_chart(self, solution: List[Tuple[int, int]]):
"""
Visualize schedule as Gantt chart
"""
# Calculate schedule times
job_completion = [0] * self.n_jobs
machine_available = [0] * self.n_machines
schedule = []
for job_id, op_idx in solution:
machine, proc_time = self.jobs[job_id][op_idx]
start_time = max(job_completion[job_id], machine_available[machine])
end_time = start_time + proc_time
schedule.append({
'job': job_id,
'machine': machine,
'start': start_time,
'duration': proc_time
})
job_completion[job_id] = end_time
machine_available[machine] = end_time
# Plot Gantt chart
fig, ax = plt.subplots(figsize=(14, 8))
colors = plt.cm.Set3(np.linspace(0, 1, self.n_jobs))
for task in schedule:
ax.barh(task['machine'],
task['duration'],
left=task['start'],
height=0.6,
color=colors[task['job']],
edgecolor='black',
linewidth=1.5)
# Add job label
ax.text(task['start'] + task['duration'] / 2,
task['machine'],
f"J{task['job']}",
ha='center',
va='center',
fontsize=10,
fontweight='bold')
ax.set_xlabel('Time', fontsize=12)
ax.set_ylabel('Machine', fontsize=12)
ax.set_title(f'Job Shop Schedule - Makespan: {self.best_makespan:.2f}',
fontsize=14)
ax.set_yticks(range(self.n_machines))
ax.set_yticklabels([f'M{i}' for i in range(self.n_machines)])
ax.grid(True, axis='x', alpha=0.3)
# Legend
legend_elements = [plt.Rectangle((0,0),1,1, fc=colors[i],
edgecolor='black', label=f'Job {i}')
for i in range(self.n_jobs)]
ax.legend(handles=legend_elements, loc='upper right')
plt.tight_layout()
plt.show()
def plot_convergence(self):
"""Plot makespan and temperature convergence"""
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 8))
# Makespan convergence
ax1.plot(self.makespan_history, linewidth=2, color='blue')
ax1.set_xlabel('Cooling Iteration')
ax1.set_ylabel('Best Makespan')
ax1.set_title('SA Convergence')
ax1.grid(True, alpha=0.3)
# Temperature schedule
ax2.plot(self.temperature_history, linewidth=2, color='red')
ax2.set_xlabel('Cooling Iteration')
ax2.set_ylabel('Temperature')
ax2.set_title('Cooling Schedule')
ax2.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
# Example usage
if __name__ == "__main__":
random.seed(42)
np.random.seed(42)
# Example: 5 jobs, 3 machines
# Each job: list of (machine, processing_time)
jobs = [
[(0, 3), (1, 2), (2, 2)], # Job 0
[(0, 2), (2, 1), (1, 4)], # Job 1
[(1, 4), (2, 3)], # Job 2
[(0, 3), (1, 1), (2, 2)], # Job 3
[(2, 2), (0, 3), (1, 1)] # Job 4
]
# Create SA optimizer
sa_jssp = SimulatedAnnealingJobShop(
jobs=jobs,
initial_temp=1000,
final_temp=0.1,
cooling_rate=0.95,
iterations_per_temp=100
)
# Optimize
result = sa_jssp.optimize()
print("\nOptimization Complete!")
print(f"Best Makespan: {result['best_makespan']:.2f}")
# Visualize
sa_jssp.plot_gantt_chart(result['best_solution'])
sa_jssp.plot_convergence()
Tabu Search
Core Concepts
Memory-Based Search:
- Tabu List: Remember recent moves to avoid cycling
- Aspiration Criterion: Override tabu if solution is best found
- Intensification: Focus on promising regions
- Diversification: Explore new regions
Tabu Search for Facility Location
import numpy as np
import random
from typing import List, Set, Tuple, Dict
import matplotlib.pyplot as plt
class TabuSearchFacilityLocation:
"""
Tabu Search for Capacitated Facility Location Problem
Problem: Select facilities to open and assign customers to minimize
total cost (fixed + transportation)
"""
def __init__(self,
customer_locations: List[Tuple[float, float]],
facility_locations: List[Tuple[float, float]],
customer_demands: List[float],
facility_capacities: List[float],
fixed_costs: List[float],
tabu_tenure: int = 10,
max_iterations: int = 500,
intensification_freq: int = 50):
"""
Initialize Tabu Search for CFLP
customer_locations: list of (x, y) for customers
facility_locations: list of (x, y) for facilities
customer_demands: demand for each customer
facility_capacities: capacity of each facility
fixed_costs: fixed cost to open each facility
tabu_tenure: how long moves stay tabu
"""
self.customers = customer_locations
self.facilities = facility_locations
self.demands = customer_demands
self.capacities = facility_capacities
self.fixed_costs = fixed_costs
self.n_customers = len(customer_locations)
self.n_facilities = len(facility_locations)
# Tabu Search parameters
self.tabu_tenure = tabu_tenure
self.max_iterations = max_iterations
self.intensification_freq = intensification_freq
# Distance matrix
self.distance_matrix = self._calculate_distances()
# Results
self.best_solution = None
self.best_cost = float('inf')
self.cost_history = []
def _calculate_distances(self) -> np.ndarray:
"""Calculate distances between customers and facilities"""
distances = np.zeros((self.n_customers, self.n_facilities))
for i, cust in enumerate(self.customers):
for j, fac in enumerate(self.facilities):
dist = np.sqrt((cust[0] - fac[0])**2 + (cust[1] - fac[1])**2)
distances[i, j] = dist
return distances
def _create_initial_solution(self) -> Dict:
"""
Create initial solution
Solution: {
'open_facilities': set of facility indices,
'assignments': list of facility assigned to each customer
}
"""
# Greedy construction: open cheapest facilities until demand met
sorted_facilities = sorted(range(self.n_facilities),
key=lambda f: self.fixed_costs[f] /
self.capacities[f])
open_facilities = set()
total_capacity = 0
total_demand = sum(self.demands)
for fac in sorted_facilities:
if total_capacity >= total_demand:
break
open_facilities.add(fac)
total_capacity += self.capacities[fac]
# Assign customers to nearest open facility
assignments = self._assign_customers(open_facilities)
return {
'open_facilities': open_facilities,
'assignments': assignments
}
def _assign_customers(self, open_facilities: Set[int]) -> List[int]:
"""
Assign customers to open facilities (nearest with capacity)
Returns: list where assignments[i] = facility for customer i
"""
assignments = [-1] * self.n_customers
facility_load = {f: 0 for f in open_facilities}
# Sort customers by minimum distance to any open facility
customer_order = sorted(range(self.n_customers),
key=lambda c: min(self.distance_matrix[c, f]
for f in open_facilities))
for cust in customer_order:
# Find nearest facility with available capacity
sorted_facs = sorted(open_facilities,
key=lambda f: self.distance_matrix[cust, f])
for fac in sorted_facs:
if facility_load[fac] + self.demands[cust] <= self.capacities[fac]:
assignments[cust] = fac
facility_load[fac] += self.demands[cust]
break
# If no capacity, assign to nearest (will be penalized)
if assignments[cust] == -1:
assignments[cust] = sorted_facs[0]
facility_load[sorted_facs[0]] += self.demands[cust]
return assignments
def _calculate_cost(self, solution: Dict) -> float:
"""
Calculate total cost for solution
Cost = fixed costs + transportation costs + capacity penalty
"""
open_facilities = solution['open_facilities']
assignments = solution['assignments']
# Fixed costs
fixed_cost = sum(self.fixed_costs[f] for f in open_facilities)
# Transportation costs
transport_cost = sum(
self.distance_matrix[c, assignments[c]]
for c in range(self.n_customers)
)
# Capacity penalty
facility_load = {f: 0 for f in open_facilities}
for c, f in enumerate(assignments):
facility_load[f] = facility_load.get(f, 0) + self.demands[c]
penalty = 0
for f in open_facilities:
if facility_load[f] > self.capacities[f]:
penalty += (facility_load[f] - self.capacities[f]) * 1000
return fixed_cost + transport_cost + penalty
def _get_neighborhood(self, solution: Dict,
tabu_list: List) -> List[Dict]:
"""
Generate neighborhood solutions
Two types of moves:
1. Open/close a facility
2. Reassign a customer
"""
neighbors = []
open_facilities = solution['open_facilities'].copy()
# Move 1: Open/close facility
for fac in range(self.n_facilities):
if fac in open_facilities:
# Try closing
if len(open_facilities) > 1: # Keep at least one open
new_open = open_facilities - {fac}
if ('close', fac) not in tabu_list:
new_assignments = self._assign_customers(new_open)
neighbors.append({
'open_facilities': new_open,
'assignments': new_assignments,
'move': ('close', fac)
})
else:
# Try opening
new_open = open_facilities | {fac}
if ('open', fac) not in tabu_list:
new_assignments = self._assign_customers(new_open)
neighbors.append({
'open_facilities': new_open,
'assignments': new_assignments,
'move': ('open', fac)
})
return neighbors
def optimize(self) -> Dict:
"""
Run Tabu Search
Returns: optimization results
"""
print(f"Initializing Tabu Search for Facility Location...")
print(f"Customers: {self.n_customers}, Facilities: {self.n_facilities}")
print(f"Max Iterations: {self.max_iterations}")
# Initialize
current_solution = self._create_initial_solution()
current_cost = self._calculate_cost(current_solution)
self.best_solution = current_solution.copy()
self.best_cost = current_cost
tabu_list = []
for iteration in range(self.max_iterations):
# Generate neighbors
neighbors = self._get_neighborhood(current_solution, tabu_list)
if not neighbors:
print("No valid neighbors found!")
break
# Evaluate neighbors
neighbor_costs = [self._calculate_cost(n) for n in neighbors]
# Select best non-tabu neighbor (or best if aspiration met)
best_neighbor_idx = np.argmin(neighbor_costs)
best_neighbor = neighbors[best_neighbor_idx]
best_neighbor_cost = neighbor_costs[best_neighbor_idx]
# Aspiration criterion: accept if best overall
if best_neighbor_cost < self.best_cost:
current_solution = best_neighbor
current_cost = best_neighbor_cost
self.best_solution = current_solution.copy()
self.best_cost = best_neighbor_cost
print(f"Iteration {iteration + 1}: "
f"New best solution found! Cost = {self.best_cost:.2f}")
else:
# Accept best neighbor
current_solution = best_neighbor
current_cost = best_neighbor_cost
# Update tabu list
if 'move' in best_neighbor:
tabu_list.append(best_neighbor['move'])
if len(tabu_list) > self.tabu_tenure:
tabu_list.pop(0)
self.cost_history.append(self.best_cost)
# Progress report
if (iteration + 1) % 50 == 0:
print(f"Iteration {iteration + 1}: "
f"Best Cost = {self.best_cost:.2f}, "
f"Current Cost = {current_cost:.2f}")
return {
'best_solution': self.best_solution,
'best_cost': self.best_cost,
'cost_history': self.cost_history,
'open_facilities': self.best_solution['open_facilities'],
'assignments': self.best_solution['assignments']
}
def plot_solution(self, solution: Dict):
"""Visualize the facility location solution"""
plt.figure(figsize=(12, 10))
open_facilities = solution['open_facilities']
assignments = solution['assignments']
# Plot customers
cust_x = [c[0] for c in self.customers]
cust_y = [c[1] for c in self.customers]
plt.scatter(cust_x, cust_y, c='blue', s=100,
marker='o', label='Customers', zorder=3)
# Plot facilities
colors = plt.cm.Set1(np.linspace(0, 1, len(open_facilities)))
for i, fac_idx in enumerate(open_facilities):
fac = self.facilities[fac_idx]
plt.scatter(fac[0], fac[1], c=[colors[i]], s=400,
marker='s', edgecolors='black', linewidth=2,
label=f'Facility {fac_idx}', zorder=4)
# Draw assignments
for cust_idx, assigned_fac in enumerate(assignments):
if assigned_fac == fac_idx:
cust = self.customers[cust_idx]
plt.plot([cust[0], fac[0]], [cust[1], fac[1]],
color=colors[i], alpha=0.3, linewidth=1, zorder=1)
plt.xlabel('X Coordinate')
plt.ylabel('Y Coordinate')
plt.title(f'Facility Location - Total Cost: {self.best_cost:.2f}')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()
def plot_convergence(self):
"""Plot cost convergence"""
plt.figure(figsize=(10, 6))
plt.plot(self.cost_history, linewidth=2)
plt.xlabel('Iteration')
plt.ylabel('Best Cost')
plt.title('Tabu Search Convergence')
plt.grid(True, alpha=0.3)
plt.show()
# Example usage
if __name__ == "__main__":
random.seed(42)
np.random.seed(42)
# Generate problem instance
n_customers = 40
n_facilities = 8
customers = [(random.uniform(0, 100), random.uniform(0, 100))
for _ in range(n_customers)]
facilities = [(random.uniform(0, 100), random.uniform(0, 100))
for _ in range(n_facilities)]
demands = [random.uniform(10, 30) for _ in range(n_customers)]
capacities = [random.uniform(150, 250) for _ in range(n_facilities)]
fixed_costs = [random.uniform(500, 1500) for _ in range(n_facilities)]
# Create Tabu Search optimizer
ts_flp = TabuSearchFacilityLocation(
customer_locations=customers,
facility_locations=facilities,
customer_demands=demands,
facility_capacities=capacities,
fixed_costs=fixed_costs,
tabu_tenure=15,
max_iterations=300
)
# Optimize
result = ts_flp.optimize()
print("\nOptimization Complete!")
print(f"Total Cost: {result['best_cost']:.2f}")
print(f"Open Facilities: {result['open_facilities']}")
# Visualize
…(truncated)