Facility Location

Solve facility location problems for optimal warehouse and distribution center placement.

Installation

pip install numpy scipy pulp scikit-learn geopy

Quick Start

from pulp import LpProblem, LpMinimize, LpVariable, lpSum

# Candidate locations for new warehouses
candidates = ['Loc_A', 'Loc_B', 'Loc_C', 'Loc_D']

# Customer demand locations
customers = {
    'C1': {'demand': 100, 'coords': (40.7, -74.0)},
    'C2': {'demand': 150, 'coords': (34.0, -118.2)},
    'C3': {'demand': 80, 'coords': (41.9, -87.6)}
}

Classic Location Problems

P-Median Problem

def p_median(candidates, customers, distances, p):
    """
    Locate p facilities to minimize total weighted distance.
    distances: dict of (candidate, customer) -> distance
    """
    prob = LpProblem("P_Median", LpMinimize)

    # y[j] = 1 if facility opened at candidate j
    y = {j: LpVariable(f"open_{j}", cat='Binary') for j in candidates}

    # x[i,j] = 1 if customer i assigned to facility j
    x = {(i, j): LpVariable(f"assign_{i}_{j}", cat='Binary')
         for i in customers for j in candidates}

    # Objective: minimize weighted distance
    prob += lpSum([customers[i]['demand'] * distances[j, i] * x[i, j]
                   for i in customers for j in candidates])

    # Exactly p facilities
    prob += lpSum([y[j] for j in candidates]) == p

    # Each customer assigned to exactly one facility
    for i in customers:
        prob += lpSum([x[i, j] for j in candidates]) == 1

    # Can only assign to open facilities
    for i in customers:
        for j in candidates:
            prob += x[i, j] <= y[j]

    prob.solve()

    return {
        'facilities': [j for j in candidates if y[j].value() > 0.5],
        'assignments': {i: [j for j in candidates if x[i, j].value() > 0.5][0]
                       for i in customers}
    }

P-Center Problem

def p_center(candidates, customers, distances, p):
    """
    Locate p facilities to minimize maximum distance to any customer.
    """
    prob = LpProblem("P_Center", LpMinimize)

    # Maximum distance variable
    z = LpVariable("max_distance", lowBound=0)

    y = {j: LpVariable(f"open_{j}", cat='Binary') for j in candidates}
    x = {(i, j): LpVariable(f"assign_{i}_{j}", cat='Binary')
         for i in customers for j in candidates}

    # Objective: minimize maximum distance
    prob += z

    # Exactly p facilities
    prob += lpSum([y[j] for j in candidates]) == p

    # Each customer assigned to one facility
    for i in customers:
        prob += lpSum([x[i, j] for j in candidates]) == 1

    # Assignment only to open facilities
    for i in customers:
        for j in candidates:
            prob += x[i, j] <= y[j]

    # z >= distance for each assignment
    for i in customers:
        for j in candidates:
            prob += z >= distances[j, i] * x[i, j]

    prob.solve()

    return {
        'max_distance': z.value(),
        'facilities': [j for j in candidates if y[j].value() > 0.5]
    }

Capacitated Facility Location

def capacitated_facility_location(candidates, customers, distances,
                                   fixed_costs, capacities):
    """
    Locate facilities with capacity constraints.
    """
    prob = LpProblem("CFLP", LpMinimize)

    y = {j: LpVariable(f"open_{j}", cat='Binary') for j in candidates}
    x = {(i, j): LpVariable(f"flow_{i}_{j}", lowBound=0)
         for i in customers for j in candidates}

    # Objective: fixed costs + transportation costs
    prob += (
        lpSum([fixed_costs[j] * y[j] for j in candidates]) +
        lpSum([distances[j, i] * x[i, j]
               for i in customers for j in candidates])
    )

    # Demand satisfaction
    for i in customers:
        prob += lpSum([x[i, j] for j in candidates]) == customers[i]['demand']

    # Capacity constraints
    for j in candidates:
        prob += lpSum([x[i, j] for i in customers]) <= capacities[j] * y[j]

    prob.solve()

    return {
        'facilities': [j for j in candidates if y[j].value() > 0.5],
        'flows': {(i, j): x[i, j].value()
                  for i in customers for j in candidates
                  if x[i, j].value() > 0},
        'total_cost': prob.objective.value()
    }

Geographic Analysis

from geopy.distance import geodesic

def calculate_distances(candidates, customers):
    """Calculate distances between all candidate-customer pairs."""
    distances = {}
    for j, cand_data in candidates.items():
        for i, cust_data in customers.items():
            cand_coords = cand_data['coords']
            cust_coords = cust_data['coords']
            distances[j, i] = geodesic(cand_coords, cust_coords).kilometers
    return distances

def find_demand_centroid(customers):
    """Find demand-weighted centroid of customers."""
    total_demand = sum(c['demand'] for c in customers.values())

    centroid_lat = sum(c['coords'][0] * c['demand']
                       for c in customers.values()) / total_demand
    centroid_lon = sum(c['coords'][1] * c['demand']
                       for c in customers.values()) / total_demand

    return (centroid_lat, centroid_lon)

def generate_candidate_grid(bounds, grid_size):
    """Generate grid of candidate locations within bounds."""
    min_lat, max_lat, min_lon, max_lon = bounds

    candidates = {}
    lat_step = (max_lat - min_lat) / grid_size
    lon_step = (max_lon - min_lon) / grid_size

    for i in range(grid_size + 1):
        for j in range(grid_size + 1):
            lat = min_lat + i * lat_step
            lon = min_lon + j * lon_step
            candidates[f"Grid_{i}_{j}"] = {'coords': (lat, lon)}

    return candidates

Coverage Models

def set_covering_location(candidates, customers, distances, coverage_radius):
    """
    Find minimum facilities to cover all customers within radius.
    """
    prob = LpProblem("SetCovering", LpMinimize)

    y = {j: LpVariable(f"open_{j}", cat='Binary') for j in candidates}

    # Objective: minimize number of facilities
    prob += lpSum([y[j] for j in candidates])

    # Coverage constraints
    for i in customers:
        # Find which candidates can cover this customer
        covering = [j for j in candidates if distances[j, i] <= coverage_radius]
        prob += lpSum([y[j] for j in covering]) >= 1

    prob.solve()

    return {
        'facilities': [j for j in candidates if y[j].value() > 0.5],
        'num_facilities': sum(y[j].value() for j in candidates)
    }

def maximal_covering_location(candidates, customers, distances,
                               coverage_radius, p):
    """
    Locate p facilities to maximize covered demand.
    """
    prob = LpProblem("MaxCovering", LpMinimize)

    y = {j: LpVariable(f"open_{j}", cat='Binary') for j in candidates}
    z = {i: LpVariable(f"covered_{i}", cat='Binary') for i in customers}

    # Objective: maximize covered demand (minimize uncovered)
    prob += lpSum([customers[i]['demand'] * (1 - z[i]) for i in customers])

    # Exactly p facilities
    prob += lpSum([y[j] for j in candidates]) == p

    # Coverage logic
    for i in customers:
        covering = [j for j in candidates if distances[j, i] <= coverage_radius]
        prob += z[i] <= lpSum([y[j] for j in covering])

    prob.solve()

    return {
        'facilities': [j for j in candidates if y[j].value() > 0.5],
        'covered_demand': sum(customers[i]['demand'] * z[i].value()
                             for i in customers)
    }

Location Scoring

def score_location(candidate, customers, competitors=None, factors=None):
    """Score a candidate location based on multiple factors."""
    score = 0

    # Distance to demand
    demand_score = sum(c['demand'] / calculate_distance(candidate, c)
                       for c in customers.values())
    score += demand_score * factors.get('demand_weight', 1.0)

    # Competition avoidance
    if competitors:
        comp_distances = [calculate_distance(candidate, comp)
                         for comp in competitors]
        min_comp_dist = min(comp_distances) if comp_distances else float('inf')
        score += min_comp_dist * factors.get('competition_weight', 0.5)

    # Infrastructure (example factors)
    if factors.get('highway_access'):
        score += factors['highway_access'] * 10
    if factors.get('labor_availability'):
        score += factors['labor_availability'] * 5

    return score