metaheuristic-optimization

メタヒューリスティック最適化

サプライチェーン向けのメタヒューリスティック最適化アルゴリズムの専門家として、自然に着想を得た手法やヒューリスティック手法を用いて、厳密解法では現実的に解くことが難しい複雑で大規模な最適化問題に対し、高品質な解を見つける支援をします。

初期評価

メタヒューリスティックを適用する前に、以下を把握します。

1. 問題の特性

   • どのような意思決定を最適化する必要がありますか？（ルーティング、スケジューリング、割り当て）
   • 問題の規模はどの程度ですか？（変数、制約）
   • 厳密解法を使えない理由は何ですか？（時間がかかりすぎる、規模が大きすぎる、NP 困難）
   • 必要な解の品質はどの程度ですか？（最適解か、十分に良い解か）

2. 問題の構造

   • 決定変数は離散ですか、それとも連続ですか？
   • 制約はハード制約（必ず満たす必要がある）ですか、それともソフト制約（ペナルティを課す）ですか？
   • 目的は単一ですか、それとも複数ですか？
   • 問題のランドスケープは滑らかですか、険しいですか、それとも多峰性ですか？

3. 計算リソース

   • 最適化に利用できる時間はどの程度ですか？（秒、分、時間）
   • 並列計算を利用できますか？
   • リアルタイム最適化ですか、それともオフライン最適化ですか？
   • メモリに制約はありますか？

4. 現在のアプローチ

   • 現在使用しているヒューリスティックやルールはありますか？
   • 上回るべきベースラインの性能はありますか？
   • 既知の優良な解はありますか？
   • 活用できるドメイン固有の知識はありますか？

---

メタヒューリスティックアルゴリズムの種類

集団ベースの手法

特徴：

• 複数の候補解を維持します
• 解空間を広く探索します
• 多峰性の問題に適しています
• 例：遺伝的アルゴリズム、粒子群最適化

利点：

• 局所最適解に陥りにくいです
• 並列化を活用できます
• 解の多様なプールを維持できます

欠点：

• 収束が遅くなります
• 計算オーバーヘッドが大きくなります
• 調整するパラメーターが増えます

軌道ベースの手法

特徴：

• 単一の解を反復的に改善します
• 集中的な局所探索を行います
• 過去の解を記憶します
• 例：シミュレーテッドアニーリング、タブー探索

利点：

• 収束が速いです
• メモリ使用量が少ないです
• 実装が簡単です

欠点：

• 局所最適解に陥ることがあります
• 探索の範囲が狭くなります
• 逐次実行になります

---

遺伝的アルゴリズム（GA）

基本概念

進化プロセス：

1. 集団：候補解（染色体）の集合
2. 選択：適応度に基づいて親を選びます
3. 交叉：親を組み合わせて子を作ります
4. 突然変異：多様性を維持するためにランダムな変更を加えます
5. 置換：新しい世代が古い世代に置き換わります

車両経路問題に対する GA

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()

---

シミュレーテッドアニーリング（SA）

基本概念

冶金学のアナロジー：

• 焼きなまし：金属を加熱した後、ゆっくりと冷却します
• 温度：ランダム性を制御するパラメーターです
• 冷却スケジュール：時間の経過とともに温度を下げます
• 受理基準：悪化する解も確率的に受け入れます

ジョブショップスケジューリングに対する SA

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()

---

タブー探索

基本概念

メモリベースの探索：

• タブーリスト：循環を避けるため、最近の移動を記憶します
• 願望基準：これまでに見つかった最良解であれば、タブー状態を解除します
• 集中化：有望な領域に重点を置きます
• 多様化：新しい領域を探索します

施設配置問題に対するタブー探索

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
    ts_flp.plot_solution(result['best_solution'])
    ts_flp.plot_convergence()

---

その他のメタヒューリスティック

蟻コロニー最適化（ACO）

概念：アリの採餌行動をシミュレートします

• アリは経路にフェロモンを残します
• 後続のアリは、フェロモンの濃い経路をたどります
• 時間の経過とともに、良い解が強化されます

用途：TSP、VRP、ネットワークルーティング

主なパラメーター：

• α（アルファ）：フェロモンの重要度
• β（ベータ）：ヒューリスティック情報の重要度
• ρ（ロー）：蒸発率
• Q：フェロモンの付与量

粒子群最適化（PSO）

概念：鳥の群れや魚の群泳をシミュレートします

• 粒子が探索空間内を移動します
• 個体最良解と大域最良解の影響を受けます
• 速度と位置を更新します

用途：連続最適化、パラメーター調整

更新規則：

v_new = w * v + c1 * r1 * (pbest - x) + c2 * r2 * (gbest - x)
x_new = x + v_new

可変近傍探索（VNS）

概念：近傍構造を体系的に変化させます

• 複数の近傍で局所探索を行います
• シェイクによって局所最適解から脱出します
• 構造化された探索を行います

用途：経路探索、スケジューリング、立地

---

ハイブリッドメタヒューリスティクス

遺伝的アルゴリズム＋局所探索

def hybrid_ga_local_search(problem, ga_params, ls_params):
    """
    Hybrid GA with local search

    GA for global exploration + Local search for exploitation
    """

    # Initialize GA population
    population = initialize_population(ga_params['pop_size'])

    for generation in range(ga_params['generations']):
        # GA operations
        offspring = []
        for _ in range(ga_params['offspring_size']):
            parents = tournament_selection(population, 2)
            child = crossover(parents[0], parents[1])
            child = mutate(child, ga_params['mutation_rate'])
            offspring.append(child)

        # Apply local search to best offspring
        offspring_fitness = [evaluate(ind) for ind in offspring]
        best_idx = np.argmin(offspring_fitness)

        # Local search improves best individual
        offspring[best_idx] = local_search(offspring[best_idx], ls_params)

        # Update population
        population = elitism_replacement(population, offspring,
                                        ga_params['elite_size'])

    return best_individual(population)

---

パラメーター調整

重要なパラメーター

集団ベース（GA、PSO）：

• 集団サイズ：ほとんどの問題で50～200
• 突然変異率：0.05～0.20
• 交叉率：0.7～0.9
• エリートサイズ：集団の5～20%

軌跡ベース（SA、TS）：

• 初期温度（SA）：予想コストの10～1000倍
• 冷却率（SA）：0.90～0.99
• 禁忌期間（TS）：7～20反復
• 近傍サイズ：問題に依存

自動調整

from scipy.optimize import differential_evolution

def tune_ga_parameters(problem_instances):
    """
    Automated parameter tuning using meta-optimization
    """

    def objective(params):
        pop_size, mutation_rate, crossover_rate = params

        # Run GA with these parameters on validation set
        avg_quality = 0
        for problem in problem_instances:
            result = run_ga(problem,
                          pop_size=int(pop_size),
                          mutation_rate=mutation_rate,
                          crossover_rate=crossover_rate)
            avg_quality += result['best_fitness']

        return avg_quality / len(problem_instances)

    # Bounds for parameters
    bounds = [
        (50, 200),    # population size
        (0.01, 0.3),  # mutation rate
        (0.6, 0.95)   # crossover rate
    ]

    # Meta-optimize
    result = differential_evolution(objective, bounds)

    optimal_params = {
        'pop_size': int(result.x[0]),
        'mutation_rate': result.x[1],
        'crossover_rate': result.x[2]
    }

    return optimal_params

---

ツールとライブラリ

Pythonライブラリ

メタヒューリスティックフレームワーク：

• pymoo：多目的最適化
• deap：分散進化アルゴリズム
• scipy.optimize：差分進化、盆地ホッピング
• optuna：ハイパーパラメーター最適化

問題固有：

• python-tsp：TSPソルバー
• or-tools (Google)：経路探索、スケジューリング
• jmetal：多目的メタヒューリスティクス

並列コンピューティング：

• joblib：関数呼び出しの並列化
• dask：分散コンピューティング
• ray：スケーラブルなコンピューティング

商用ソフトウェア

• FICO Xpress Mosel：メタヒューリスティックモジュール
• CPLEX with Heuristics：MIPヒューリスティクス
• LocalSolver：局所探索ソルバー

---

よくある課題と解決策

課題：早期収束

問題：

• 集団の多様性が急速に失われます
• 局所最適解に陥ります

解決策：

• 突然変異率を上げます
• 多様性を維持する仕組みを導入します
• 再スタート戦略を採用します
• パラメーター制御を適応的に行います
• ニッチング手法を用います

課題：収束の遅さ

問題：

• 良好な解を見つけるのに時間がかかりすぎます
• 活用が不十分です

解決策：

• 局所探索とハイブリッド化します
• 初期化を改善します（貪欲構築）
• 選択圧を強めます
• 近傍サイズを適応的に調整します
• 問題固有のオペレーターを使用します

課題：パラメーターへの感度

問題：

• パラメーターによって性能が大きく変動します
• 調整が困難です

解決策：

• パラメーター制御を適応的に行います
• パラメーター不要のメタヒューリスティクスを使用します
• 自動調整（メタ最適化）を行います
• 文献で確立されたデフォルト値を使用します

課題：改善の停滞

問題：

• アルゴリズムの改善が止まります
• 探索を使い果たしています

解決策：

• 摂動／再スタート機構を導入します
• 近傍構造を変更します（VNS）
• 多様化戦略を採用します
• マルチスタート方式を使用します
• 複数のアルゴリズムを組み合わせます

---

ベストプラクティス

アルゴリズムの選択

メタヒューリスティクスを選ぶ場合：

• 問題がNP困難です
• 厳密解法では時間がかかりすぎます
• 準最適解で許容できます
• 利用できる時間が限られています
• 問題の構造が不明です

アルゴリズム選択ガイド：

問題の種類	推奨アルゴリズム
順列（TSP、VRP）	OX交叉を用いるGA、ACO
二値（施設立地）	GA、タブーサーチ
連続	PSO、差分進化、SA
スケジューリング	タブーサーチ、SA、GA
多目的	NSGA-II、MOEA/D

実装のヒント

1. シンプルに始める：まず基本的なバージョンを実装します
2. 問題に関する知識を活用する：カスタムオペレーターを設計します
3. ベンチマークを行う：単純なヒューリスティクスと比較します
4. 進捗を記録する：最良解の推移をログに記録します
5. 複数回実行する：異なるシード値で10～30回実行します
6. 解を検証する：制約が満たされているか確認します

---

出力形式

メタヒューリスティクスレポートのテンプレート

エグゼクティブサマリー：

• 問題の説明
• 選択したアルゴリズムとその理由
• 最良解の品質
• 計算時間

問題の定式化：

• 決定変数
• 目的関数
• 制約条件
• 問題の規模

アルゴリズムの設定：

パラメーター	値	根拠
アルゴリズム	遺伝的アルゴリズム	順列問題に適しています
集団サイズ	100	多様性と速度のバランス
世代数	500	第400世代までに収束を確認
突然変異率	0.15	この種の問題で標準的な値
交叉	順序交叉（OX）	順列の有効性を維持

結果：

指標	値
最良目的値	$125,456
平均（10回の実行）	$128,234
標準偏差	$1,823
最短実行時間	45秒
下界からのギャップ	5.2%

解の品質：

• 収束プロット
• 解の可視化
• ベースラインとの比較
• 統計分析（複数回の実行）

推奨事項：

• 最良解を導入する
• 期待される改善
• さらなる最適化の機会

---

質問

さらに詳しい情報が必要な場合：

1. どのような問題を解決しようとしていますか？
2. 厳密最適化を実用的でなくしている要因は何ですか？
3. 問題の規模はどれくらいですか？（変数、制約条件）
4. 最適化にどのくらいの時間をかけられますか？
5. 解の品質は重要ですか、それとも「十分に良い」解で問題ありませんか？
6. 制約条件はありますか？ハード制約ですか、それともソフト制約ですか？
7. 活用できる分野固有の知識はありますか？
8. これは一度限りの最適化ですか、それとも繰り返し行う最適化ですか？

---

関連スキル

• optimization-modeling：厳密最適化手法向け
• multi-objective-optimization：パレート最適化向け
• constraint-programming：制約ベースの探索向け
• route-optimization：VRP アプリケーション向け
• production-scheduling：スケジューリング問題向け
• network-design：ネットワーク最適化向け
• prescriptive-analytics：意思決定の最適化向け
• supply-chain-automation：最適化の自動化向け