Solve Network Optimization Problems with Python, NumPy, Scikit-learn, and OSRM

Introduction

Many organizations face the problem of choosing optimal locations for facilities — whether warehouses, travel hubs, corporate headquarters, or distribution centers. The goal is to minimize travel time, improve logistics efficiency, or maximize accessibility.

This is a general location-optimization problem, applicable in many contexts:

Warehouses: Minimize delivery times to retail outlets.
Airline hubs: Choose airports to reduce layover times.
Corporate headquarters: Centralize offices to reduce employee commute times.
Emergency services: Place fire stations or hospitals to minimize response times.

While the domain changes, the core problem remains: finding the point in a network that minimizes travel time to multiple destinations.

Understanding the Problem

We need a programmatic solution that:

Takes a set of destination points (outlets, offices, hubs).
Finds a candidate location for a new facility that minimizes travel time to these points.
Accounts for real-world road networks rather than straight-line distances.
Can scale to multiple destinations.

Challenges include:

Roads are not straight lines; freeways may provide faster travel than geographically closer routes.
Traffic varies by time of day.
When there are many destinations, brute-force search becomes infeasible.

Reasoning and Methodology

We considered several approaches:

1. Geometric and Travel-Time Midpoints

For two destinations, a simple strategy is to pick the travel-time midpoint along the route.
For three destinations, the geometric median (the point minimizing total Euclidean distance) is a good starting guess.

2. Hill-Climb Search

After identifying a candidate location, we can perform local optimization by exploring nearby points along the road network.
Each candidate is evaluated using actual travel time obtained from a routing engine.

3. Clustering for Multiple Destinations

For 4 or more destinations, we cluster outlets to reduce the problem’s dimensionality.
Each cluster is represented by a centroid (geometric mean), and the optimization is applied to these centroids.

Available Technologies

Several technologies can be applied to this problem:

Technology	Role	Suitability
Python	Programming language	Excellent for prototyping and data manipulation.
OSRM (Open Source Routing Machine)	Road network routing	Self-hosted, fast, accurate routing using OSM data.
Google Maps API / Mapbox	Routing services	Accurate, easy to use, but paid after free quota.
KMeans clustering (scikit-learn)	Clustering outlets	Simple, fast, Euclidean approximation suitable when travel-time table is unavailable.
Geometric Median / Weiszfeld algorithm	Centrality for 3+ points	Robust starting point for optimization.

Why we chose OSRM + Python + Hill-Climb

Free and self-hosted: No reliance on paid services.
Fast: OSRM computes routes in milliseconds.
Flexible: Allows snapping candidate points to real roads.
Scalable: Works for dozens of destinations without N² travel-time queries.

Why we didn’t choose

Travel-time-based clustering: Requires precomputed N×N tables, infeasible for large datasets.
Pure KMeans on lat/lon: Doesn’t consider roads, but we compensate using hill climb.
Complex optimization (simulated annealing, genetic algorithms): Overkill for small to medium datasets; more complex to maintain.

Definitions of Methods Employed

NumPy, short for Numerical Python, is a fundamental Python library for scientific computing. It provides support for large, multi-dimensional arrays and matrices, along with a comprehensive collection of high-level mathematical functions to operate on these arrays. Here, we use it to find the mean of a collection of coordinates.

Scikit-learn (often referred to as sklearn) is a popular, open-source Python library for machine learning. It provides a wide range of efficient tools for predictive data analysis. Here we will use its K-Means algorithm for clustering for dimensionality reduction.

OSRM is an opensource routing machine, akin to google maps. It will perform road network routing for us so that distances are not calculated just geometrically but by real-world travel-time calculations along the quickest routes on roads. This software is downloaded and installed locally. It is free and self-hosted and provides accurate routing using downloaded OSM data.

Other helpful definitions:

Route Timewise Midpoint:
Finds the point along the fastest route between two destinations where cumulative travel time is half the total.
Geometric Median:
Point that minimizes total Euclidean distance to multiple destinations, computed using the Weiszfeld algorithm.
Hill-Climb Search:
Local search method that iteratively moves the candidate point in small steps along the road network, keeping the move if it improves total travel time.
Clustering:
Groups multiple destinations into 2–3 clusters using Euclidean KMeans. Each cluster is reduced to a centroid for subsequent optimization.
Time-of-Day Multipliers:
Travel times are weighted to simulate traffic congestion based on hour of the day.

Step-by-Step Implementation

1. Prerequisites

Python 3.10+
Packages: requests, numpy, scikit-learn
OSRM installed locally with relevant map data.

Install packages:

pip install requests numpy scikit-learn

Install OSRM:

# Install OSRM backend
git clone https://github.com/Project-OSRM/osrm-backend.git
cd osrm-backend
mkdir build && cd build
cmake ..
cmake --build .

# Download OSM data (example: Dallas)
wget https://download.geofabrik.de/north-america/us/texas-latest.osm.pbf

# Prepare OSRM data
osrm-extract -p ../profiles/car.lua texas-latest.osm.pbf
osrm-contract texas-latest.osrm

# Run OSRM server
osrm-routed texas-latest.osrm

Explanatory Notes:

Python and the listed packages are used for both data manipulation and calling OSRM over HTTP.
OSRM must be installed and running locally because it provides fast routing and nearest-road snapping without relying on paid services.
Extracting and contracting the OSM data prepares the routing engine for efficient queries on your region of interest.

2. Design

Input: list of outlets (lat, lon) and current hour.
- Explanatory: This defines the problem space. Each outlet is a destination, and the hour is used to apply traffic congestion multipliers.
Determine number of outlets:
- 2 → timewise midpoint
  - Explanatory: For two outlets, we compute the point along the actual driving route where travel time is exactly halfway. This accounts for road network topology.
- 3 → geometric median
  - Explanatory: For three outlets, the geometric median provides a balanced central point minimizing total distance. This serves as a good initial guess for optimization.
- ≥4 → cluster → centroids → geometric median / midpoint.
  - Explanatory: For four or more outlets, clustering reduces complexity. We approximate groups of outlets as single points, then optimize using these representative points.
Perform hill climb:
- Snap candidate points to nearest road node.
  - Explanatory: This ensures the candidate location is accessible on the road network, not just a geographic coordinate.
- Evaluate total travel time to all outlets.
  - Explanatory: We use OSRM to compute realistic driving times, which may differ from straight-line distances.
- Apply small perturbations iteratively to find a local optimum.
  - Explanatory: Random perturbations explore nearby road nodes to minimize total travel time further. This mimics climbing toward the “best” point on a travel-time landscape.
Return best candidate point and total travel time.
- Explanatory: The algorithm outputs the recommended facility location and the cumulative weighted travel time. This allows stakeholders to evaluate efficiency gains.

3. Full Python Code

Save as warehouse_location_optimizer.py:

# warehouse_location_optimizer.py

import requests
import math
import random
import numpy as np
from datetime import datetime
from sklearn.cluster import KMeans

# -------------------------------
# CONFIGURATION
# -------------------------------

OSRM_URL = "http://localhost:5000"  # self-hosted OSRM
OSRM_TIMEOUT = 10  # seconds

# Time-of-day travel-time multipliers
TIME_MULTIPLIERS = [
    (0, 6, 1.0),
    (6, 9, 3.0),
    (9, 11, 1.5),
    (11, 14, 2.0),
    (14, 16, 1.5),
    (16, 19, 3.0),
    (19, 21, 1.5),
    (21, 24, 1.0)
]

def get_time_multiplier(hour: int) -> float:
    """Return multiplier for travel time based on the hour of the day."""
    for start, end, mult in TIME_MULTIPLIERS:
        if start <= hour < end:
            return mult
    return 1.0  # fallback

# -------------------------------
# OSRM UTILITIES
# -------------------------------

def osrm_nearest(lat, lon):
    """Snap a point to the nearest road node using OSRM /nearest."""
    url = f"{OSRM_URL}/nearest/v1/driving/{lon},{lat}?number=1"
    resp = requests.get(url, timeout=OSRM_TIMEOUT).json()
    if "waypoints" not in resp:
        # fallback if OSRM fails
        return (lat, lon)
    snapped = resp["waypoints"][0]["location"]
    return (snapped[1], snapped[0])  # lat, lon

def osrm_route_time(coords, multiplier=1.0):
    """Get travel time in seconds between coords [(lat, lon), ...] using OSRM."""
    coord_str = ";".join([f"{lon},{lat}" for lat, lon in coords])
    url = f"{OSRM_URL}/route/v1/driving/{coord_str}?overview=false"
    resp = requests.get(url, timeout=OSRM_TIMEOUT).json()
    if "routes" not in resp or len(resp["routes"]) == 0:
        return float("inf")
    # Apply time-of-day multiplier here
    return resp["routes"][0]["duration"] * multiplier

# -------------------------------
# MATHEMATICAL UTILITIES
# -------------------------------

def geometric_median(points, eps=1e-5):
    """Compute geometric median (Weiszfeld algorithm)."""
    points = np.array(points)
    guess = points.mean(axis=0)
    while True:
        dist = np.linalg.norm(points - guess, axis=1)
        if np.any(dist == 0):
            return tuple(guess)
        weights = 1 / dist
        new_guess = (points * weights[:, None]).sum(axis=0) / weights.sum()
        if np.linalg.norm(new_guess - guess) < eps:
            return tuple(new_guess)
        guess = new_guess

def route_timewise_midpoint(outlet_a, outlet_b):
    """Compute approximate OSRM travel-time midpoint between two outlets."""
    coord_str = f"{outlet_a[1]},{outlet_a[0]};{outlet_b[1]},{outlet_b[0]}"
    url = f"{OSRM_URL}/route/v1/driving/{coord_str}?overview=full&geometries=geojson&annotations=true"
    resp = requests.get(url, timeout=OSRM_TIMEOUT).json()
    route = resp["routes"][0]
    geom_coords = [(c[1], c[0]) for c in route["geometry"]["coordinates"]]

    # Use cumulative duration along the route to find halfway by travel time
    leg = route["legs"][0]
    durations = leg.get("annotation", {}).get("duration", [])
    if durations:
        cum = np.cumsum(durations)
        half = cum[-1] / 2
        idx = int(np.searchsorted(cum, half))
        prev_cum = cum[idx-1] if idx > 0 else 0
        frac = (half - prev_cum)/durations[idx] if durations[idx] > 0 else 0
        a = geom_coords[idx]
        b = geom_coords[idx+1]
        midpoint = (a[0]+(b[0]-a[0])*frac, a[1]+(b[1]-a[1])*frac)
        return midpoint
    else:
        # fallback: geometric mean
        lat = (outlet_a[0]+outlet_b[0])/2
        lon = (outlet_a[1]+outlet_b[1])/2
        return (lat, lon)

# -------------------------------
# HILL CLIMB SEARCH
# -------------------------------

def hill_climb(start, outlets, hour, steps=20, radius=0.01):
    """
    Perform hill-climb optimization:
    - start: initial candidate point (lat, lon)
    - outlets: list of outlet coordinates
    - hour: int (0-23) for time-of-day multiplier
    - steps: max number of hill climb steps
    - radius: max random perturbation in lat/lon degrees
    """
    multiplier = get_time_multiplier(hour)
    # Snap starting point to road node
    current = osrm_nearest(*start)
    best_time = sum(osrm_route_time([current, o], multiplier) for o in outlets)

    for _ in range(steps):
        improved = False
        # Explore nearby candidate points
        for _ in range(8):
            dlat = (random.random()-0.5)*radius
            dlon = (random.random()-0.5)*radius
            candidate = osrm_nearest(current[0]+dlat, current[1]+dlon)
            total_time = sum(osrm_route_time([candidate, o], multiplier) for o in outlets)
            if total_time < best_time:
                best_time = total_time
                current = candidate
                improved = True
        if not improved:
            break  # local optimum reached
    return current, best_time

# -------------------------------
# FULL PIPELINE
# -------------------------------

def find_warehouse(outlets, hour=None):
    """
    Main pipeline to determine warehouse location.
    Steps:
    1. Determine number of outlets
    2. Choose initial candidate point:
        - 2 outlets: OSRM time-wise midpoint
        - 3 outlets: geometric median
        - 4+ outlets: Euclidean KMeans clustering (2-3 clusters) → cluster centroids
    3. Hill climb starting from candidate to optimize travel time
    4. Returns best warehouse coordinates and total travel time
    """
    if hour is None:
        hour = datetime.now().hour
    n = len(outlets)

    # --- 2 outlets ---
    if n == 2:
        # Candidate = midpoint along route by travel time
        start = route_timewise_midpoint(outlets[0], outlets[1])
        return hill_climb(start, outlets, hour)

    # --- 3 outlets ---
    elif n == 3:
        # Candidate = geometric median
        start = geometric_median(outlets)
        return hill_climb(start, outlets, hour)

    # --- 4+ outlets ---
    else:
        # --- CLUSTERING STEP ---
        # Euclidean KMeans is used for clustering here:
        # - Reason: travel-time table does not exist, so we cannot cluster by road distance
        # - Euclidean coordinates provide a fast approximation
        k = 2 if n < 10 else 3
        kmeans = KMeans(n_clusters=k, n_init=10)
        labels = kmeans.fit_predict(outlets)
        cluster_means = []
        for c in range(k):
            # Compute the mean lat/lon of outlets in this cluster
            pts = np.array([outlets[i] for i in range(n) if labels[i] == c])
            cluster_means.append(tuple(pts.mean(axis=0)))

        # Initial candidate for hill climb:
        if len(cluster_means) == 2:
            start = route_timewise_midpoint(cluster_means[0], cluster_means[1])
        else:
            start = geometric_median(cluster_means)

        # Hill climb to optimize travel time
        return hill_climb(start, outlets, hour)

# -------------------------------
# EXAMPLE USAGE
# -------------------------------

if __name__=="__main__":
    #Example outlet coordinates
    outlets = [
        (31.978, -96.121),
        (32.150, -94.839),
        (32.160, -96.000),
        (39.940, -95.810)
    ]
    
    hour_of_day = datetime.now().hour
    best_location, total_time = find_warehouse(outlets, hour=hour_of_day)

    # Unpack best location
    lat, lon = best_location

    # Build Google Maps URL
    maps_url = f"https://www.google.com/maps?q={lat:.6f},{lon:.6f}"

    # Print the recommended warehouse location
    print("\nRecommended Warehouse Location:")
    print(f"Latitude: {lat:.6f}, Longitude: {lon:.6f}")
    print(f"Google Maps: {maps_url}")

Explanatory Notes:

All helper functions are defined clearly: snapping to roads, computing travel-time midpoints, geometric median, clustering, and hill climb.
The find_warehouse function integrates all steps into a single pipeline that automatically chooses the right method based on the number of outlets.
Time-of-day multipliers are applied dynamically during hill climb, making the solution context-aware.
Replace the example outlet coordinates in the code with your own outlet coordinates.
You can run the code with:
- python warehouse_location_optimizer.py

4. Did you clearly see how we used NumPy and Scikit-learn?

A Review of the use of NumPy

NumPy is primarily used for numerical operations, vector math, and efficient array manipulations. In the code, it is used in several key places:

Geometric Median Calculation
- For 3 outlets (or cluster centroids), we compute the geometric median by minimizing the sum of distances.
- NumPy arrays make it easy to:
  - Store outlet coordinates as (lat, lon) arrays
  - Compute differences between candidate points and all outlets
  - Compute Euclidean distances efficiently
  - Perform iterative optimization or convergence checks

coords = np.array(outlets)
diff = coords - candidate_point
distances = np.linalg.norm(diff, axis=1)

Hill-Climb Perturbations
- During the hill-climb search, small candidate movements are applied as vectors.
- NumPy allows you to add/subtract arrays and compute new candidate points efficiently.

perturbation = np.random.uniform(-0.001, 0.001, size=2)
candidate_point = candidate_point + perturbation

Cluster Centroid Calculations
- After clustering (for 4+ outlets), we compute centroids using NumPy’s mean over arrays.

centroid = np.mean(cluster_points, axis=0)

A review of the use of Scikit-learn

Scikit-learn is used mainly for clustering when there are 4 or more outlets. This reduces computational complexity by grouping nearby outlets together and treating the cluster centroids as “representative” outlets.

KMeans Clustering
- We use sklearn.cluster.KMeans to create 2 or 3 clusters from the list of outlets.
- This helps reduce the multi-outlet problem to a simpler 2- or 3-point problem, which we can handle with the geometric median or midpoint approach.

from sklearn.cluster import KMeans

kmeans = KMeans(n_clusters=num_clusters, random_state=42)
kmeans.fit(outlet_coords)
centroids = kmeans.cluster_centers_

Why It’s Important
- Without clustering, evaluating every possible candidate location for many outlets becomes computationally expensive.
- Clustering provides a smart approximation: optimize warehouse locations for cluster centroids instead of each individual outlet.

NumPy and Scikit-learn are super useful in Network Optimization Problems

Library	Purpose in Code
NumPy	Numerical operations, vector math, geometric median, hill-climb perturbations, centroids calculation
scikit-learn	Clustering outlets into 2–3 groups for larger datasets using KMeans, reducing complexity

5. Why there are slightly different warehouse locations with each run

1. Hill-Climb Search Has Randomness

The hill-climb algorithm we implemented doesn’t just sit on the exact geometric midpoint (or route midpoint).
It “snaps” to nearby road nodes and explores alternatives.
Depending on how that exploration starts (e.g., a random nearby node or floating-point rounding differences), the path of the search can change slightly.
Hill climbing by nature can converge to slightly different “local optima” depending on the initial seed.

2. Clustering for 4+ Outlets

When you have 4 or more outlets, the pipeline uses KMeans from scikit-learn to split them into 2 or 3 clusters.
KMeans initializes its centroids randomly (unless you fix random_state).
This means each run can cluster outlets slightly differently, shifting the starting point of the search and ultimately the warehouse location.

3. Floating Point & OSRM Routing

OSRM itself can sometimes return micro-variations in travel times because of internal approximations or changes in edge snapping.
These differences, though tiny, can affect which neighbor the hill-climb prefers at a step.

6. Making It Production Ready

We leave the code with you as a good place to start. Below are enhancements you should consider to make the code production ready.

Caching OSRM calls:
- Explanatory: Repeated route queries can be stored to reduce redundant HTTP requests, which improves performance.
Parallelization:
- Explanatory: Candidate evaluations during hill climb can run in parallel to speed up convergence for multiple destinations.
Error handling & retries:
- Explanatory: OSRM requests can fail due to network or server issues; retries ensure robustness.
Scalability:
- Explanatory: For larger datasets, dynamic clustering or grid-based preselection reduces computational load without sacrificing accuracy.
Monitoring & logging:
- Explanatory: Logging OSRM response times, candidate evaluations, and convergence ensures maintainability and easier troubleshooting.

7. Further work to be considered

The solution provided here assumes a hub in the middle star topology. The flow in such a topology has traffic always passing through the hub in the middle. It would be interesting and extremely useful to build the solution assuming other topologies. In particular it would be good to build the solution for a circular path, ring topology. Intuition would tell you that the later is a simpler problem.

Conclusion

This approach provides a practical, free, and extensible solution for optimizing facility locations based on actual travel times and road networks. While it uses simple heuristics (clustering, hill climb), it’s sufficient for city-scale deployments and can be extended for larger datasets or different facility types.

Primary technology: Python + OSRM
Key techniques: Geometric median, route midpoint, hill climb, clustering, time-of-day multipliers
Applications: Warehouses, corporate offices, travel hubs, emergency services, distribution centers

By combining mathematical optimization with real-world routing, organizations can make data-driven, cost-effective decisions for their physical locations.