Facility Location Problem (FLP)

Problem Classification

1. Uncapacitated Facility Location Problem (UFLP)

Description:

• Decide which facilities to open (no capacity limits)
• Assign customers to open facilities
• Minimize total fixed + transportation costs

Characteristics:

• Each facility can serve unlimited demand
• Most basic FLP variant
• NP-hard but solvable for moderate instances

Applications:

• Initial network design
• Strategic long-term planning
• High-level site selection

2. Capacitated Facility Location Problem (CFLP)

Description:

• Facilities have capacity constraints
• More realistic than UFLP
• Each customer can be served by multiple facilities

Characteristics:

• Capacity constraints make problem harder
• May require more facilities than UFLP
• More complex solution methods needed

Applications:

• Warehouse network design
• Production facility location
• Service center placement

3. p-Median Problem

Description:

• Locate exactly p facilities
• Minimize total (weighted) distance from customers to facilities
• No fixed costs, no capacity constraints

Characteristics:

• Number of facilities predetermined
• Focus on minimizing transportation distance
• Classic location science problem

Applications:

• Emergency service location
• Retail store placement
• School/hospital location

4. p-Center Problem

Description:

• Locate p facilities
• Minimize the maximum distance from any customer to nearest facility
• Minimax objective (equity-focused)

Characteristics:

• Ensures no customer too far from service
• Different objective than p-median
• Good for emergency services

Applications:

• Ambulance station location
• Fire station placement
• Emergency response planning

Mathematical Formulations

Uncapacitated Facility Location Problem (UFLP)

Sets:

• I = {1, ..., m}: Set of potential facility locations
• J = {1, ..., n}: Set of customers

Parameters:

• f_i: Fixed cost to open facility at location i
• c_{ij}: Cost to serve customer j from facility i (transportation cost)
• d_j: Demand of customer j

Decision Variables:

• y_i ∈ {0,1}: 1 if facility i is opened, 0 otherwise
• x_{ij} ∈ [0,1]: Fraction of customer j's demand served by facility i

Objective Function:

Minimize: Σ_i f_i * y_i + Σ_i Σ_j c_{ij} * d_j * x_{ij}
          \_____________/   \___________________________/
          Fixed costs       Transportation costs

Constraints:

1. Demand satisfaction: Each customer fully served
   Σ_i x_{ij} = 1,  ∀j ∈ J

2. Facility opening: Can only serve from open facilities
   x_{ij} ≤ y_i,  ∀i ∈ I, ∀j ∈ J

3. Binary facility decisions:
   y_i ∈ {0,1},  ∀i ∈ I

4. Assignment variables:
   x_{ij} ≥ 0,  ∀i ∈ I, ∀j ∈ J

Complexity: NP-hard

Capacitated Facility Location Problem (CFLP)

Additional Parameters:

• Q_i: Capacity of facility i

Modified Constraints:

1. Demand satisfaction:
   Σ_i x_{ij} = 1,  ∀j ∈ J

2. Capacity constraints:
   Σ_j d_j * x_{ij} ≤ Q_i * y_i,  ∀i ∈ I

3. Opening constraints:
   x_{ij} ≤ y_i,  ∀i ∈ I, ∀j ∈ J

4. Binary variables:
   y_i ∈ {0,1},  ∀i ∈ I
   x_{ij} ≥ 0,  ∀i ∈ I, ∀j ∈ J

p-Median Problem

Objective:

Minimize: Σ_i Σ_j c_{ij} * x_{ij}

Constraints:

1. Each customer assigned to exactly one facility:
   Σ_i x_{ij} = 1,  ∀j ∈ J

2. Exactly p facilities opened:
   Σ_i y_i = p

3. Assignment only to open facilities:
   x_{ij} ≤ y_i,  ∀i ∈ I, ∀j ∈ J

4. Binary variables:
   y_i ∈ {0,1},  ∀i ∈ I
   x_{ij} ∈ {0,1},  ∀i ∈ I, ∀j ∈ J

Exact Solution Methods

1. MIP Formulation with PuLP (UFLP)

from pulp import *
import numpy as np

def solve_uflp(fixed_costs, transport_costs, demands):
    """
    Solve Uncapacitated Facility Location Problem

    Args:
        fixed_costs: list of fixed costs for each facility
        transport_costs: 2D array [facilities x customers] of unit transport costs
        demands: list of customer demands

    Returns:
        dict with optimal solution
    """
    m = len(fixed_costs)  # Number of potential facilities
    n = len(demands)       # Number of customers

    # Create problem
    prob = LpProblem("UFLP", LpMinimize)

    # Decision variables
    # y[i] = 1 if facility i is opened
    y = LpVariable.dicts("facility", range(m), cat='Binary')

    # x[i,j] = fraction of customer j's demand served by facility i
    x = LpVariable.dicts("assign",
                         [(i, j) for i in range(m) for j in range(n)],
                         lowBound=0, upBound=1, cat='Continuous')

    # Objective: Minimize total cost (fixed + transportation)
    prob += (
        lpSum([fixed_costs[i] * y[i] for i in range(m)]) +
        lpSum([transport_costs[i][j] * demands[j] * x[i,j]
               for i in range(m) for j in range(n)]),
        "Total_Cost"
    )

    # Constraints

    # 1. Each customer must be fully served
    for j in range(n):
        prob += (
            lpSum([x[i,j] for i in range(m)]) == 1,
            f"Demand_Customer_{j}"
        )

    # 2. Can only serve from open facilities
    for i in range(m):
        for j in range(n):
            prob += (
                x[i,j] <= y[i],
                f"Open_Facility_{i}_{j}"
            )

    # Solve
    import time
    start_time = time.time()
    prob.solve(PULP_CBC_CMD(msg=1, timeLimit=300))
    solve_time = time.time() - start_time

    # Extract solution
    if LpStatus[prob.status] in ['Optimal', 'Feasible']:
        open_facilities = [i for i in range(m) if y[i].varValue > 0.5]

        assignments = {}
        for j in range(n):
            assignments[j] = []
            for i in range(m):
                if x[i,j].varValue > 0.01:  # Threshold for numerical issues
                    assignments[j].append((i, x[i,j].varValue))

        # Calculate cost breakdown
        fixed_cost_total = sum(fixed_costs[i] for i in open_facilities)
        transport_cost_total = sum(
            transport_costs[i][j] * demands[j] * x[i,j].varValue
            for i in range(m) for j in range(n)
        )

        return {
            'status': LpStatus[prob.status],
            'total_cost': value(prob.objective),
            'fixed_cost': fixed_cost_total,
            'transport_cost': transport_cost_total,
            'open_facilities': open_facilities,
            'num_facilities': len(open_facilities),
            'assignments': assignments,
            'solve_time': solve_time
        }
    else:
        return {
            'status': LpStatus[prob.status],
            'total_cost': None,
            'open_facilities': [],
            'solve_time': solve_time
        }


# Example usage
if __name__ == "__main__":
    # Problem data: 5 potential facilities, 10 customers
    np.random.seed(42)

    # Fixed costs to open each facility
    fixed_costs = [5000, 4500, 6000, 5500, 4800]

    # Transportation costs (facility x customer)
    # Lower cost = closer proximity
    transport_costs = np.array([
        [10, 15, 8, 20, 12, 18, 22, 14, 16, 11],
        [18, 12, 16, 14, 10, 15, 20, 25, 13, 19],
        [14, 20, 18, 10, 16, 12, 15, 18, 22, 14],
        [22, 18, 14, 16, 20, 10, 12, 16, 14, 18],
        [16, 14, 20, 18, 14, 16, 10, 12, 15, 20]
    ])

    # Customer demands
    demands = [100, 150, 80, 120, 90, 110, 130, 95, 105, 125]

    result = solve_uflp(fixed_costs, transport_costs, demands)

    print(f"\n{'='*70}")
    print(f"UNCAPACITATED FACILITY LOCATION PROBLEM - SOLUTION")
    print(f"{'='*70}")
    print(f"Status: {result['status']}")
    print(f"Total Cost: ${result['total_cost']:,.2f}")
    print(f"  Fixed Costs: ${result['fixed_cost']:,.2f}")
    print(f"  Transport Costs: ${result['transport_cost']:,.2f}")
    print(f"\nFacilities Opened: {result['num_facilities']}")
    print(f"Facility IDs: {result['open_facilities']}")
    print(f"\nSolve Time: {result['solve_time']:.2f} seconds")

    print(f"\nCustomer Assignments:")
    for customer_id, assignment in result['assignments'].items():
        print(f"  Customer {customer_id} (demand={demands[customer_id]}):")
        for facility_id, fraction in assignment:
            print(f"    → Facility {facility_id}: {fraction*100:.1f}%")

2. Capacitated Facility Location Problem (CFLP)

def solve_cflp(fixed_costs, transport_costs, demands, capacities):
    """
    Solve Capacitated Facility Location Problem

    Args:
        fixed_costs: list of fixed costs for each facility
        transport_costs: 2D array [facilities x customers]
        demands: list of customer demands
        capacities: list of facility capacities

    Returns:
        dict with optimal solution
    """
    m = len(fixed_costs)
    n = len(demands)

    # Create problem
    prob = LpProblem("CFLP", LpMinimize)

    # Decision variables
    y = LpVariable.dicts("facility", range(m), cat='Binary')
    x = LpVariable.dicts("assign",
                         [(i, j) for i in range(m) for j in range(n)],
                         lowBound=0, cat='Continuous')

    # Objective
    prob += (
        lpSum([fixed_costs[i] * y[i] for i in range(m)]) +
        lpSum([transport_costs[i][j] * demands[j] * x[i,j]
               for i in range(m) for j in range(n)]),
        "Total_Cost"
    )

    # Constraints

    # 1. Each customer fully served
    for j in range(n):
        prob += (
            lpSum([x[i,j] for i in range(m)]) == 1,
            f"Demand_{j}"
        )

    # 2. Capacity constraints at each facility
    for i in range(m):
        prob += (
            lpSum([demands[j] * x[i,j] for j in range(n)]) <= capacities[i] * y[i],
            f"Capacity_{i}"
        )

    # 3. Can only serve from open facilities
    for i in range(m):
        for j in range(n):
            prob += (
                x[i,j] <= y[i],
                f"Open_{i}_{j}"
            )

    # Solve
    start_time = time.time()
    prob.solve(PULP_CBC_CMD(msg=1, timeLimit=300))
    solve_time = time.time() - start_time

    # Extract solution
    if LpStatus[prob.status] in ['Optimal', 'Feasible']:
        open_facilities = [i for i in range(m) if y[i].varValue > 0.5]

        # Calculate utilization for each open facility
        utilization = {}
        for i in open_facilities:
            used_capacity = sum(demands[j] * x[i,j].varValue for j in range(n))
            utilization[i] = (used_capacity / capacities[i]) * 100

        assignments = {}
        for j in range(n):
            assignments[j] = []
            for i in range(m):
                if x[i,j].varValue > 0.01:
                    assignments[j].append((i, x[i,j].varValue))

        return {
            'status': LpStatus[prob.status],
            'total_cost': value(prob.objective),
            'open_facilities': open_facilities,
            'num_facilities': len(open_facilities),
            'utilization': utilization,
            'assignments': assignments,
            'solve_time': solve_time
        }
    else:
        return {
            'status': LpStatus[prob.status],
            'solve_time': solve_time
        }


# Example usage
fixed_costs = [8000, 7500, 9000, 8500, 7800]
transport_costs = np.array([
    [10, 15, 8, 20, 12, 18, 22, 14, 16, 11],
    [18, 12, 16, 14, 10, 15, 20, 25, 13, 19],
    [14, 20, 18, 10, 16, 12, 15, 18, 22, 14],
    [22, 18, 14, 16, 20, 10, 12, 16, 14, 18],
    [16, 14, 20, 18, 14, 16, 10, 12, 15, 20]
])
demands = [100, 150, 80, 120, 90, 110, 130, 95, 105, 125]
capacities = [400, 350, 450, 380, 420]  # Facility capacities

result = solve_cflp(fixed_costs, transport_costs, demands, capacities)

print(f"\n{'='*70}")
print(f"CAPACITATED FACILITY LOCATION PROBLEM - SOLUTION")
print(f"{'='*70}")
print(f"Status: {result['status']}")
print(f"Total Cost: ${result['total_cost']:,.2f}")
print(f"Facilities Opened: {result['num_facilities']}")

print(f"\nFacility Utilization:")
for facility_id in result['open_facilities']:
    print(f"  Facility {facility_id}: {result['utilization'][facility_id]:.1f}% "
          f"(capacity={capacities[facility_id]})")

3. p-Median Problem

def solve_p_median(distances, demands, p):
    """
    Solve p-Median Problem

    Locate exactly p facilities to minimize total weighted distance

    Args:
        distances: 2D array of distances [facilities x customers]
        demands: customer demands (weights)
        p: number of facilities to open

    Returns:
        optimal solution
    """
    m = len(distances)     # Potential facility locations
    n = len(demands)       # Customers

    prob = LpProblem("p_Median", LpMinimize)

    # Decision variables
    y = LpVariable.dicts("facility", range(m), cat='Binary')
    x = LpVariable.dicts("assign",
                         [(i, j) for i in range(m) for j in range(n)],
                         cat='Binary')

    # Objective: Minimize total weighted distance
    prob += (
        lpSum([distances[i][j] * demands[j] * x[i,j]
               for i in range(m) for j in range(n)]),
        "Total_Weighted_Distance"
    )

    # Constraints

    # 1. Each customer assigned to exactly one facility
    for j in range(n):
        prob += (
            lpSum([x[i,j] for i in range(m)]) == 1,
            f"Assign_{j}"
        )

    # 2. Exactly p facilities opened
    prob += (
        lpSum([y[i] for i in range(m)]) == p,
        "p_Facilities"
    )

    # 3. Assignment only to open facilities
    for i in range(m):
        for j in range(n):
            prob += (
                x[i,j] <= y[i],
                f"Open_{i}_{j}"
            )

    # Solve
    start_time = time.time()
    prob.solve(PULP_CBC_CMD(msg=0))
    solve_time = time.time() - start_time

    if LpStatus[prob.status] in ['Optimal', 'Feasible']:
        open_facilities = [i for i in range(m) if y[i].varValue > 0.5]

        assignments = {}
        for j in range(n):
            for i in range(m):
                if x[i,j].varValue > 0.5:
                    assignments[j] = i
                    break

        return {
            'status': LpStatus[prob.status],
            'total_distance': value(prob.objective),
            'open_facilities': open_facilities,
            'assignments': assignments,
            'solve_time': solve_time
        }
    else:
        return {'status': LpStatus[prob.status]}


# Example: Locate 3 facilities among 8 candidates to serve 12 customers
np.random.seed(42)

# Generate random coordinates
facility_coords = np.random.rand(8, 2) * 100
customer_coords = np.random.rand(12, 2) * 100

# Calculate Euclidean distance matrix
distances = np.zeros((8, 12))
for i in range(8):
    for j in range(12):
        distances[i][j] = np.linalg.norm(facility_coords[i] - customer_coords[j])

demands = [100, 150, 80, 120, 90, 110, 130, 95, 105, 125, 115, 140]
p = 3  # Open exactly 3 facilities

result = solve_p_median(distances, demands, p)

print(f"\n{'='*70}")
print(f"p-MEDIAN PROBLEM - SOLUTION")
print(f"{'='*70}")
print(f"Status: {result['status']}")
print(f"Total Weighted Distance: {result['total_distance']:,.2f}")
print(f"Facilities Opened (p={p}): {result['open_facilities']}")
print(f"\nCustomer Assignments:")
for customer_id, facility_id in result['assignments'].items():
    dist = distances[facility_id][customer_id]
    print(f"  Customer {customer_id} → Facility {facility_id} "
          f"(distance={dist:.2f}, demand={demands[customer_id]})")

Greedy Heuristics

1. Greedy Add Algorithm (p-Median)

def greedy_add_p_median(distances, demands, p):
    """
    Greedy heuristic for p-Median Problem

    Iteratively add facility that gives maximum cost reduction

    Time complexity: O(p * m * n)
    Quality: Typically within 10-20% of optimal

    Args:
        distances: distance matrix [facilities x customers]
        demands: customer demands
        p: number of facilities to open

    Returns:
        heuristic solution
    """
    m, n = distances.shape

    open_facilities = []
    assignments = {}

    # Initialize: assign all customers to closest potential facility
    for j in range(n):
        assignments[j] = np.argmin(distances[:, j])

    # Greedily add p facilities
    for _ in range(p):
        best_facility = None
        best_improvement = 0

        for i in range(m):
            if i in open_facilities:
                continue

            # Calculate improvement if we add facility i
            improvement = 0
            for j in range(n):
                current_dist = distances[assignments[j]][j]
                new_dist = distances[i][j]

                if new_dist < current_dist:
                    improvement += demands[j] * (current_dist - new_dist)

            if improvement > best_improvement:
                best_improvement = improvement
                best_facility = i

        if best_facility is not None:
            open_facilities.append(best_facility)

            # Update assignments
            for j in range(n):
                current_dist = distances[assignments[j]][j]
                new_dist = distances[best_facility][j]

                if new_dist < current_dist:
                    assignments[j] = best_facility

    # Calculate total cost
    total_cost = sum(
        distances[assignments[j]][j] * demands[j]
        for j in range(n)
    )

    return {
        'open_facilities': open_facilities,
        'assignments': assignments,
        'total_distance': total_cost,
        'method': 'Greedy Add'
    }

2. Greedy Drop Algorithm

def greedy_drop_p_median(distances, demands, p):
    """
    Greedy drop heuristic for p-Median

    Start with all facilities, iteratively close facility
    with least impact until p remain

    Args:
        distances: distance matrix
        demands: customer demands
        p: number of facilities to keep

    Returns:
        heuristic solution
    """
    m, n = distances.shape

    # Start with all facilities open
    open_facilities = set(range(m))

    # Close facilities until only p remain
    while len(open_facilities) > p:
        best_facility_to_close = None
        min_cost_increase = float('inf')

        for i in open_facilities:
            # Test closing facility i
            test_facilities = open_facilities - {i}

            # Calculate cost with this facility closed
            cost = 0
            for j in range(n):
                min_dist = min(distances[f][j] for f in test_facilities)
                cost += min_dist * demands[j]

            if cost < min_cost_increase:
                min_cost_increase = cost
                best_facility_to_close = i

        open_facilities.remove(best_facility_to_close)

    # Final assignments
    assignments = {}
    for j in range(n):
        best_facility = min(open_facilities,
                          key=lambda f: distances[f][j])
        assignments[j] = best_facility

    total_cost = sum(
        distances[assignments[j]][j] * demands[j]
        for j in range(n)
    )

    return {
        'open_facilities': list(open_facilities),
        'assignments': assignments,
        'total_distance': total_cost,
        'method': 'Greedy Drop'
    }

3. Savings-Based Heuristic (UFLP)

def savings_heuristic_uflp(fixed_costs, transport_costs, demands):
    """
    Savings-based heuristic for UFLP

    Compute savings from opening each facility and open
    facilities in order of savings

    Args:
        fixed_costs: fixed facility costs
        transport_costs: transportation costs
        demands: customer demands

    Returns:
        heuristic solution
    """
    m = len(fixed_costs)
    n = len(demands)

    # Calculate savings for each facility
    # Savings = reduction in transport cost - fixed cost
    savings = []

    for i in range(m):
        # Calculate transport cost if facility i serves customers
        transport_saving = 0

        # Assume without facility i, customers served by most expensive option
        for j in range(n):
            # Cost to serve customer j from facility i
            cost_from_i = transport_costs[i][j] * demands[j]

            # Cost to serve from next best facility
            other_costs = [transport_costs[k][j] * demands[j]
                          for k in range(m) if k != i]

            if other_costs:
                min_other_cost = min(other_costs)
                transport_saving += max(0, min_other_cost - cost_from_i)

        net_saving = transport_saving - fixed_costs[i]
        savings.append((net_saving, i))

    # Sort by savings (descending)
    savings.sort(reverse=True)

    # Open facilities with positive savings
    open_facilities = []
    for saving, facility_id in savings:
        if saving > 0:
            open_facilities.append(facility_id)

    # If no positive savings, open at least one facility
    if not open_facilities:
        # Open facility with lowest total cost for all customers
        best_facility = None
        min_cost = float('inf')

        for i in range(m):
            cost = fixed_costs[i] + sum(
                transport_costs[i][j] * demands[j] for j in range(n)
            )
            if cost < min_cost:
                min_cost = cost
                best_facility = i

        open_facilities = [best_facility]

    # Assign customers to nearest open facility
    assignments = {}
    for j in range(n):
        best_facility = min(open_facilities,
                          key=lambda i: transport_costs[i][j])
        assignments[j] = best_facility

    # Calculate total cost
    fixed_cost_total = sum(fixed_costs[i] for i in open_facilities)
    transport_cost_total = sum(
        transport_costs[assignments[j]][j] * demands[j]
        for j in range(n)
    )
    total_cost = fixed_cost_total + transport_cost_total

    return {
        'open_facilities': open_facilities,
        'assignments': assignments,
        'total_cost': total_cost,
        'fixed_cost': fixed_cost_total,
        'transport_cost': transport_cost_total,
        'method': 'Savings Heuristic'
    }

Local Search Improvements

1. Facility Swap Neighborhood

def facility_swap_local_search(distances, demands, initial_solution, p):
    """
    Local search with facility swap moves for p-Median

    Swap an open facility with a closed one if it improves objective

    Args:
        distances: distance matrix
        demands: customer demands
        initial_solution: starting solution
        p: number of facilities

    Returns:
        improved solution
    """
    m, n = distances.shape

    open_facilities = set(initial_solution['open_facilities'])
    closed_facilities = set(range(m)) - open_facilities

    def calculate_cost(facilities):
        """Calculate total cost for given facility set"""
        cost = 0
        for j in range(n):
            min_dist = min(distances[f][j] for f in facilities)
            cost += min_dist * demands[j]
        return cost

    current_cost = calculate_cost(open_facilities)
    improved = True

    while improved:
        improved = False
        best_swap = None
        best_cost = current_cost

        # Try all swaps
        for facility_in in closed_facilities:
            for facility_out in open_facilities:
                # Test swap
                test_facilities = (open_facilities - {facility_out}) | {facility_in}
                test_cost = calculate_cost(test_facilities)

                if test_cost < best_cost:
                    best_cost = test_cost
                    best_swap = (facility_in, facility_out)
                    improved = True

        if improved:
            facility_in, facility_out = best_swap
            open_facilities = (open_facilities - {facility_out}) | {facility_in}
            closed_facilities = (closed_facilities - {facility_in}) | {facility_out}
            current_cost = best_cost

    # Final assignments
    assignments = {}
    for j in range(n):
        best_facility = min(open_facilities,
                          key=lambda f: distances[f][j])
        assignments[j] = best_facility

    return {
        'open_facilities': list(open_facilities),
        'assignments': assignments,
        'total_distance': current_cost,
        'method': 'Facility Swap Local Search'
    }

2. Customer Reassignment

def customer_reassignment(open_facilities, distances, demands):
    """
    Improve solution by optimally reassigning customers

    Given open facilities, assign each customer to nearest facility

    Args:
        open_facilities: list of open facility indices
        distances: distance matrix
        demands: customer demands

    Returns:
        optimal assignments for given facilities
    """
    n = distances.shape[1]

    assignments = {}
    for j in range(n):
        best_facility = min(open_facilities,
                          key=lambda i: distances[i][j])
        assignments[j] = best_facility

    total_cost = sum(
        distances[assignments[j]][j] * demands[j]
        for j in range(n)
    )

    return {
        'assignments': assignments,
        'total_cost': total_cost
    }

Metaheuristics

1. Simulated Annealing for p-Median

import random
import math

def simulated_annealing_p_median(distances, demands, p,
                                 initial_temp=1000, cooling_rate=0.95,
                                 max_iterations=5000):
    """
    Simulated Annealing for p-Median Problem

    Args:
        distances: distance matrix
        demands: customer demands
        p: number of facilities
        initial_temp: starting temperature
        cooling_rate: temperature reduction factor
        max_iterations: maximum iterations

    Returns:
        best solution found
    """
    m, n = distances.shape

    def calculate_cost(facilities):
        cost = 0
        for j in range(n):
            min_dist = min(distances[f][j] for f in facilities)
            cost += min_dist * demands[j]
        return cost

    # Initial solution: random p facilities
    current_facilities = set(random.sample(range(m), p))
    current_cost = calculate_cost(current_facilities)

    best_facilities = current_facilities.copy()
    best_cost = current_cost

    temperature = initial_temp

    for iteration in range(max_iterations):
        # Generate neighbor: swap one facility
        closed = list(set(range(m)) - current_facilities)

        if not closed:
            break

        facility_out = random.choice(list(current_facilities))
        facility_in = random.choice(closed)

        neighbor_facilities = (current_facilities - {facility_out}) | {facility_in}
        neighbor_cost = calculate_cost(neighbor_facilities)

        delta = neighbor_cost - current_cost

        # Accept or reject
        if delta < 0 or random.random() < math.exp(-delta / temperature):
            current_facilities = neighbor_facilities
            current_cost = neighbor_cost

            if current_cost < best_cost:
                best_facilities = current_facilities.copy()
                best_cost = current_cost

        # Cool down
        temperature *= cooling_rate

        if temperature < 0.1:
            break

    # Final assignments
    assignments = {}
    for j in range(n):
        best_facility = min(best_facilities,
                          key=lambda f: distances[f][j])
        assignments[j] = best_facility

    return {
        'open_facilities': list(best_facilities),
        'assignments': assignments,
        'total_distance': best_cost,
        'method': 'Simulated Annealing'
    }

2. Genetic Algorithm for UFLP

def genetic_algorithm_uflp(fixed_costs, transport_costs, demands,
                          population_size=50, generations=200,
                          mutation_rate=0.15):
    """
    Genetic Algorithm for Uncapacitated Facility Location

    Args:
        fixed_costs: facility fixed costs
        transport_costs: transportation costs
        demands: customer demands
        population_size: GA population size
        generations: number of generations
        mutation_rate: mutation probability

    Returns:
        best solution
    """
    m = len(fixed_costs)
    n = len(demands)

    def decode_chromosome(chromosome):
        """Convert binary chromosome to open facilities"""
        return [i for i in range(m) if chromosome[i] == 1]

    def calculate_fitness(chromosome):
        """Calculate fitness (inverse of cost)"""
        open_facilities = decode_chromosome(chromosome)

        if not open_facilities:
            return 0  # Invalid solution

        # Fixed costs
        fixed_cost = sum(fixed_costs[i] for i in open_facilities)

        # Transportation costs
        transport_cost = 0
        for j in range(n):
            min_cost = min(transport_costs[i][j] for i in open_facilities)
            transport_cost += min_cost * demands[j]

        total_cost = fixed_cost + transport_cost
        return 1.0 / (1.0 + total_cost)

    def create_individual():
        """Create random chromosome"""
        # Ensure at least one facility open
        chromosome = [random.randint(0, 1) for _ in range(m)]
        if sum(chromosome) == 0:
            chromosome[random.randint(0, m-1)] = 1
        return chromosome

    def crossover(parent1, parent2):
        """Uniform crossover"""
        child = []
        for i in range(m):
            child.append(parent1[i] if random.random() < 0.5 else parent2[i])

        # Ensure at least one facility
        if sum(child) == 0:
            child[random.randint(0, m-1)] = 1

        return child

    def mutate(chromosome):
        """Bit flip mutation"""
        for i in range(m):
            if random.random() < mutation_rate:
                chromosome[i] = 1 - chromosome[i]

        # Ensure at least one facility
        if sum(chromosome) == 0:
            chromosome[random.randint(0, m-1)] = 1

        return chromosome

    # Initialize population
    population = [create_individual() for _ in range(population_size)]

    best_chromosome = None
    best_fitness = 0

    for generation in range(generations):
        # Evaluate fitness
        fitnesses = [calculate_fitness(ind) for ind in population]

        # Track best
        max_fitness = max(fitnesses)
        if max_fitness > best_fitness:
            best_fitness = max_fitness
            best_chromosome = population[fitnesses.index(max_fitness)].copy()

        # Selection and reproduction
        new_population = []

        # Elitism
        elite_count = int(0.1 * population_size)
        elite_indices = sorted(range(len(fitnesses)),
                             key=lambda i: fitnesses[i],
                             reverse=True)[:elite_count]
        new_population = [population[i].copy() for i in elite_indices]

        # Create offspring
        while len(new_population) < population_size:
            # Tournament selection
            parent1 = max(random.sample(list(zip(population, fitnesses)), 3),
                         key=lambda x: x[1])[0]
            parent2 = max(random.sample(list(zip(population, fitnesses)), 3),
                         key=lambda x: x[1])[0]

            child = crossover(parent1, parent2)
            child = mutate(child)

            new_population.append(child)

        population = new_population

    # Decode best solution
    open_facilities = decode_chromosome(best_chromosome)

    # Assign customers
    assignments = {}
    for j in range(n):
        best_facility = min(open_facilities,
                          key=lambda i: transport_costs[i][j])
        assignments[j] = best_facility

    # Calculate costs
    fixed_cost = sum(fixed_costs[i] for i in open_facilities)
    transport_cost = sum(
        transport_costs[assignments[j]][j] * demands[j]
        for j in range(n)
    )

    return {
        'open_facilities': open_facilities,
        'assignments': assignments,
        'total_cost': fixed_cost + transport_cost,
        'fixed_cost': fixed_cost,
        'transport_cost': transport_cost,
        'method': 'Genetic Algorithm'
    }

Complete Facility Location Solver

import time

class FacilityLocationSolver:
    """
    Comprehensive Facility Location Problem Solver

    Supports UFLP, CFLP, p-Median, and multiple solution methods
    """

    def __init__(self):
        self.fixed_costs = None
        self.transport_costs = None
        self.demands = None
        self.capacities = None
        self.problem_type = None

    def load_problem(self, fixed_costs=None, transport_costs=None,
                    demands=None, capacities=None, problem_type='UFLP'):
        """
        Load problem data

        Args:
            fixed_costs: list of facility fixed costs
            transport_costs: 2D array [facilities x customers]
            demands: customer demands
            capacities: facility capacities (for CFLP)
            problem_type: 'UFLP', 'CFLP', or 'p-Median'
        """
        self.fixed_costs = np.array(fixed_costs) if fixed_costs else None
        self.transport_costs = np.array(transport_costs)
        self.demands = np.array(demands)
        self.capacities = np.array(capacities) if capacities else None
        self.problem_type = problem_type

        print(f"Loaded {problem_type} problem:")
        print(f"  Facilities: {len(transport_costs)}")
        print(f"  Customers: {len(demands)}")
        print(f"  Total demand: {sum(demands)}")
        if capacities:
            print(f"  Total capacity: {sum(capacities)}")

    def solve_exact(self, time_limit=300):
        """
        Solve using exact MIP method

        Args:
            time_limit: time limit in seconds

        Returns:
            optimal/best solution
        """
        print(f"\nSolving {self.problem_type} with MIP (exact)...")

        if self.problem_type == 'UFLP':
            return solve_uflp(self.fixed_costs, self.transport_costs, self.demands)
        elif self.problem_type == 'CFLP':
            return solve_cflp(self.fixed_costs, self.transport_costs,
                            self.demands, self.capacities)
        else:
            raise ValueError(f"Exact solve not implemented for {self.problem_type}")

    def solve_heuristic(self, method='greedy', p=None):
        """
        Solve using heuristic method

        Args:
            method: 'greedy_add', 'greedy_drop', 'savings', 'sa', 'ga'
            p: number of facilities (for p-median)

        Returns:
            heuristic solution
        """
        print(f"\nSolving with {method} heuristic...")

        if method == 'greedy_add' and self.problem_type == 'p-Median':
            return greedy_add_p_median(self.transport_costs, self.demands, p)

        elif method == 'greedy_drop' and self.problem_type == 'p-Median':
            return greedy_drop_p_median(self.transport_costs, self.demands, p)

        elif method == 'savings' and self.problem_type == 'UFLP':
            return savings_heuristic_uflp(self.fixed_costs,
                                         self.transport_costs,
                                         self.demands)

        elif method == 'sa':
            if self.problem_type == 'p-Median':
                return simulated_annealing_p_median(
                    self.transport_costs, self.demands, p)
            else:
                raise ValueError("SA only for p-Median currently")

        elif method == 'ga' and self.problem_type == 'UFLP':
            return genetic_algorithm_uflp(self.fixed_costs,
                                         self.transport_costs,
                                         self.demands)

        else:
            raise ValueError(f"Unknown method {method} for {self.problem_type}")

    def compare_methods(self, methods, p=None):
        """
        Compare multiple solution methods

        Args:
            methods: list of method names
            p: number of facilities (for p-median)

        Returns:
            comparison dataframe
        """
        import pandas as pd

        results = []

        for method in methods:
            start_time = time.time()

            try:
                if method == 'exact':
                    solution = self.solve_exact()
                    cost_key = 'total_cost' if 'total_cost' in solution else 'total_distance'
                else:
                    solution = self.solve_heuristic(method, p)
                    cost_key = 'total_cost' if 'total_cost' in solution else 'total_distance'

                solve_time = time.time() - start_time

                results.append({
                    'Method': method,
                    'Cost': solution[cost_key],
                    'Facilities': len(solution['open_facilities']),
                    'Status': solution.get('status', 'Heuristic'),
                    'Time (s)': f"{solve_time:.3f}"
                })

            except Exception as e:
                print(f"  Error with {method}: {e}")

        df = pd.DataFrame(results)

        # Calculate gap from best
        if len(df) > 0:
            best_cost = df['Cost'].min()
            df['Gap %'] = ((df['Cost'] - best_cost) / best_cost * 100).round(2)

        return df

    def visualize_solution(self, solution, facility_coords, customer_coords):
        """
        Visualize facility location solution

        Args:
            solution: solution dictionary
            facility_coords: array of facility coordinates
            customer_coords: array of customer coordinates
        """
        import matplotlib.pyplot as plt

        plt.figure(figsize=(12, 8))

        # Plot all potential facilities (gray)
        for i, coord in enumerate(facility_coords):
            if i not in solution['open_facilities']:
                plt.plot(coord[0], coord[1], 's', color='lightgray',
                        markersize=10, alpha=0.5)

        # Plot open facilities (red)
        for i in solution['open_facilities']:
            coord = facility_coords[i]
            plt.plot(coord[0], coord[1], 's', color='red',
                    markersize=15, label='Open Facility' if i == solution['open_facilities'][0] else '')

        # Plot customers (blue)
        for j, coord in enumerate(customer_coords):
            plt.plot(coord[0], coord[1], 'o', color='blue',
                    markersize=8, alpha=0.7)

        # Draw assignments
        if 'assignments' in solution:
            for customer_id, assignment in solution['assignments'].items():
                if isinstance(assignment, int):
                    facility_id = assignment
                    fraction = 1.0
                else:
                    facility_id, fraction = assignment[0] if isinstance(assignment, list) else (assignment, 1.0)

                customer_coord = customer_coords[customer_id]
                facility_coord = facility_coords[facility_id]

                plt.plot([customer_coord[0], facility_coord[0]],
                        [customer_coord[1], facility_coord[1]],
                        'k-', alpha=0.2, linewidth=0.5)

        plt.xlabel('X Coordinate')
        plt.ylabel('Y Coordinate')
        plt.title(f'Facility Location Solution - {solution.get("method", "Unknown Method")}')
        plt.legend()
        plt.grid(True, alpha=0.3)
        plt.tight_layout()
        plt.show()


# Complete example
if __name__ == "__main__":
    print("="*70)
    print("FACILITY LOCATION PROBLEM - COMPREHENSIVE EXAMPLE")
    print("="*70)

    # Generate problem data
    np.random.seed(42)

    n_facilities = 10
    n_customers = 20

    # Generate coordinates
    facility_coords = np.random.rand(n_facilities, 2) * 100
    customer_coords = np.random.rand(n_customers, 2) * 100

    # Calculate distance matrix
    distances = np.zeros((n_facilities, n_customers))
    for i in range(n_facilities):
        for j in range(n_customers):
            distances[i][j] = np.linalg.norm(
                facility_coords[i] - customer_coords[j]
            )

    # Problem parameters
    fixed_costs = np.random.uniform(5000, 8000, n_facilities)
    transport_costs = distances * 10  # Cost per unit per distance
    demands = np.random.uniform(50, 200, n_customers)
    capacities = np.random.uniform(800, 1200, n_facilities)

    # Solve UFLP
    print("\n" + "="*70)
    print("UNCAPACITATED FACILITY LOCATION PROBLEM")
    print("="*70)

    solver = FacilityLocationSolver()
    solver.load_problem(fixed_costs, transport_costs, demands, problem_type='UFLP')

    # Compare methods
    comparison = solver.compare_methods(['exact', 'savings', 'ga'])
    print("\nMethod Comparison:")
    print(comparison.to_string(index=False))

    # Solve and visualize best solution
    best_solution = solver.solve_exact()
    print(f"\nBest Solution:")
    print(f"  Total Cost: ${best_solution['total_cost']:,.2f}")
    print(f"  Fixed Cost: ${best_solution['fixed_cost']:,.2f}")
    print(f"  Transport Cost: ${best_solution['transport_cost']:,.2f}")
    print(f"  Facilities Opened: {len(best_solution['open_facilities'])}")
    print(f"  Facility IDs: {best_solution['open_facilities']}")

    # Visualize
    solver.visualize_solution(best_solution, facility_coords, customer_coords)

Tools & Libraries

Python Libraries

Optimization:

• PuLP: MIP modeling (recommended for UFLP/CFLP)
• Pyomo: Advanced optimization modeling
• OR-Tools (Google): Fast heuristics and exact methods
• Gurobi: Commercial solver (academic license available)
• CPLEX: IBM commercial solver

Specialized:

• scikit-learn: Clustering for facility location
• scipy.optimize: General optimization
• geopy: Distance calculations with real coordinates

Visualization:

• matplotlib: Basic visualization
• plotly: Interactive maps
• folium: Geographic visualization
• seaborn: Statistical visualization

Commercial Software

• Gurobi Optimizer: State-of-art MIP solver
• IBM ILOG CPLEX: Enterprise optimization
• FICO Xpress: Optimization suite
• AMPL: Modeling language + solvers

Cloud Services

• Google Maps Platform: Distance Matrix API
• AWS Location Service: Location-based services
• Azure Maps: Microsoft mapping service

Common Challenges & Solutions

Challenge: Large-Scale Problems

Problem:

• 1000s of potential facilities
• 10,000s of customers
• Exact methods too slow

Solutions:

• Use heuristics (savings, greedy add/drop)
• Lagrangian relaxation
• Decomposition methods (Benders)
• Metaheuristics (SA, GA, Tabu Search)
• Preprocessing to reduce problem size
• Clustering customers before solving

Challenge: Capacity Constraints

Problem:

• CFLP is harder than UFLP
• May need many facilities
• Difficult to get tight bounds

Solutions:

• Valid inequalities to strengthen formulation
• Branch-and-price algorithms
• Column generation
• Specialized heuristics
• Consider multi-sourcing if allowed

Challenge: Multiple Objectives

Problem:

• Minimize cost AND maximize service
• Trade-off between # facilities and service quality
• Stakeholder preferences vary

Solutions:

• Multi-objective optimization
• Weighted sum method
• Epsilon-constraint method
• Pareto frontier analysis
• Scenario analysis with different weights

Challenge: Dynamic Demand

Problem:

• Demand changes over time
• Seasonal variations
• Facility decisions long-term

Solutions:

• Robust optimization
• Stochastic programming
• Scenario-based planning
• Flexible/modular facility design
• Rolling horizon approach

Challenge: Real-World Distances

Problem:

• Euclidean distances unrealistic
• Need actual road networks
• Traffic considerations

Solutions:

• Use Google Maps Distance Matrix API
• Pre-compute actual distances
• Consider time-of-day effects
• Use real logistics data

Output Format

Facility Location Solution Report

Problem Instance:

• Problem Type: Uncapacitated Facility Location
• Potential Facilities: 10
• Customers: 20
• Total Demand: 2,450 units
• Optimization: Minimize fixed + transportation costs

Optimal Solution:

Open Facilities:

Customer Assignments:

Cost Breakdown:

• Fixed costs: 26.9% of total
• Transportation costs: 73.1% of total
• Average cost per unit: $35.67

Questions to Ask

If you need more context:

1. What type of facility location problem? (UFLP, CFLP, p-Median, p-Center)
2. How many potential facility locations? How many customers?
3. What are the fixed costs to open facilities?
4. Are there capacity constraints on facilities?
5. What are the transportation/service costs?
6. Do you have coordinates or a distance matrix?
7. Is demand fixed or variable/uncertain?
8. Are there minimum or maximum number of facilities?
9. Single-sourcing (one facility per customer) or multi-sourcing allowed?
10. What's the planning horizon? (short-term tactical vs. long-term strategic)
11. Any special requirements? (coverage distance, service levels)
12. Budget constraints on total capital investment?

Related Skills

• hub-location-problem: For hub-and-spoke networks
• warehouse-location-optimization: Warehouse-specific considerations
• distribution-center-network: Distribution network design
• network-design: Broader supply chain network design
• set-covering-problem: Coverage-based facility location
• network-flow-optimization: Flow-based facility allocation
• optimization-modeling: General MIP formulation techniques
• multi-objective-optimization: Multiple criteria optimization
• stochastic-optimization: Uncertain demand scenarios
• inventory-routing-problem: Joint location-inventory-routing