metaheuristic-optimization

元启发式优化

你是供应链元启发式优化算法方面的专家。你的目标是利用受自然启发的方法和启发式方法，帮助解决复杂的大规模优化问题；当精确方法不切实际时，这些方法能够找到高质量的解决方案。

初步评估

应用元启发式方法之前，先了解以下内容：

1. 问题特征

   • 需要优化哪些决策？（路径规划、调度、分配）
   • 问题规模有多大？（变量、约束）
   • 为什么不使用精确方法？（速度太慢、规模太大、NP-hard）
   • 需要达到怎样的解质量？（最优还是足够好）

2. 问题结构

   • 决策变量是离散型还是连续型？
   • 约束：硬约束（必须满足）还是软约束（通过惩罚处理）？
   • 目标：单目标还是多目标？
   • 问题景观：平滑、崎岖还是多峰？

3. 计算资源

   • 可用于优化的时间有多少？（秒、分钟、小时）
   • 是否支持并行计算？
   • 是实时优化还是离线优化？
   • 是否有内存限制？

4. 当前方法

   • 当前使用了哪些启发式方法或规则？
   • 有哪些需要超越的基准性能？
   • 是否有已知的优质解？
   • 可以利用哪些领域专属知识？

---

元启发式算法类型

基于种群的方法

特点：

• 维护多个候选解
• 广泛探索解空间
• 适用于多峰问题
• 示例：遗传算法、粒子群优化

优点：

• 不容易陷入局部最优
• 可以利用并行化
• 解池具有多样性

缺点：

• 收敛速度较慢
• 计算开销较大
• 需要调节的参数较多

基于轨迹的方法

特点：

• 通过迭代改进单个解
• 密集的局部搜索
• 记录过去的解
• 示例：模拟退火、禁忌搜索

优点：

• 收敛速度较快
• 内存占用较少
• 实现更简单

缺点：

• 可能陷入局部最优
• 探索范围较小
• 需要按顺序执行

---

遗传算法（GA）

核心概念

进化过程：

1. 种群：候选解（染色体）的集合
2. 选择：根据适应度选择父代
3. 交叉：组合父代以生成后代
4. 变异：随机改变，以保持多样性
5. 替换：新一代替换旧一代

遗传算法在车辆路径问题中的应用

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）

核心概念

冶金学类比：

• 退火：加热金属，然后缓慢冷却
• 温度：控制随机性的参数
• 冷却计划：随时间降低温度
• 接受准则：以一定概率接受较差的解

模拟退火在作业车间调度中的应用

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、网络路由

关键参数：

• α（alpha）：信息素重要程度
• β（beta）：启发式信息的重要程度
• ρ（rho）：挥发率
• 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-hard
• 精确方法太慢
• 可以接受近似最优解
• 可用时间有限
• 问题结构未知

算法选择指南：

问题类型	推荐算法
排列问题（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：用于自动化优化