MILP: Waterborne Urban Logistics in Amsterdam

MILP: Waterborne Urban Logistics in Amsterdam#

Problem and model recap#

Amsterdam’s historic center is served by many Hotels, Restaurants and Cafés (HoReCa) that require frequent deliveries of food supplies. Traditional truck-based logistics contributes to congestion, emissions, and stress on fragile quay walls and bridges. To reduce these impacts, the city explores a logistics system that uses urban waterways for bulk transport and light vehicles (LVs) for last-mile delivery. In this system, goods are transported by vessel from a consolidation hub to Transfer Points (TPs) along the canal network. From these TPs, supplies are delivered by LVs to nearby HoReCa demand zones.

Your task is to design this network optimally by answering two key questions:

  1. Which transfer points should be established?

  2. How should HoReCa demand be assigned to the selected transfer points?

Each TP has a fixed establishment cost, a limited handling capacity, and a geographic location. Each demand zone has a known weekly demand for supplies, and a specific geographic location. Serving demand through a TP incurs transport costs (water + road) and and associated greenhouse gas emissions.

The challenge is to find a configuration of open TPs and demand assignments that:

  • satisfies all HoReCa demand,

  • respects TP capacities and planning constraints,

  • stays within budget and siting limits,

  • and minimises total system costs.

This forms a location-allocation MILP problem where strategic decisions (which TPs to open) interact with operational decisions (how demand flows are distributed).

Parameters#

  • \(d_i\): demand at zone \(i\)

  • \({cap}_j\): capacity of TP \(j\)

  • \(f_j\): fixed cost of establishing TP \(j\)

  • \(t^{\text{water}}_{0j}\): water distance from consolidation hub to TP \(j\)

  • \(t^{\text{road}}_{ij}\): LV road distance from TP \(j\) to zone \(i\)

  • \(c^{water}\): unit water transport cost

  • \(c^{road}\): unit road transport cost

  • \(e^{water}\): unit emission factor for vessels

  • \(e^{road}\): unit emission factor for LVs

  • \(D_{max}\): maximum delivery distance

  • \(\delta_{ij}\): if TP \(j\) is allowed to serve zone \(i\) (\( t^{\text{road}}_{ij} \leq D_{max}\))

  • \(B\): investment budget

  • \(K\): maximum number of TPs allowed

  • \(\lambda \ge 0\): emission penalty factor

Decision variables#

  • \(x_j \in \{0,1\}\): 1 if TP \(j\) is opened; 0 otherwise

  • \(y_{ij} \geq 0\): flow from TP \(j\) to demand zone \(i\)

Model#

\[\begin{split}\begin{align*} min \quad & Z = {\sum_{j \in J} f_j\,x_j}+{\sum_{(i,j)\in A}{CT}_{ij}.y_{ij}}+\lambda{\sum_{(i,j)\in A}{Em}_{ij}.y_{ij}}\\ \text{s.t.} \quad \\ & \sum_{j \in J} y_{ij} = d_i \qquad \forall i \in I\\ & \sum_{i \in I} y_{ij} \leq {cap}_jx_j \qquad \forall j \in J\\ & \sum_{j \in J} f_j x_j \leq B\\ & \sum_{j \in J} x_j \leq K\\ & y_{ij} \leq M \delta_{ij} \qquad \forall i \in I,\forall j \in J\\ & x_{j} \in \{0,1\}\\ & y_{ij} \geq 0 \end{align*}\end{split}\]

Data description#

In this assignment you will work with pre-processed data for Amsterdam’s historic city center. The raw inputs (exact locations and demands of individual HoReCa businesses, and candidate Transfer Points) have been processed into two main components: HoReCa demand zones and a Smaller set of candidate TPs.

You’ll work with two instances:

  • a small instance with 15 demand zones and 10 TPs

  • a larger instance with 300 demand zones and 36 TPs

You will receive separate data for each instance.

HoReCa demand zones#

The historic center contains more than 1500 individual HoReCa entities. To keep the optimisation model tractable and interpretable, these are grouped into HoReCa demand zones using spatial clustering:

  • Each HoReCa entity has a known longitude, latitude, and estimated weekly demand (derived from Pourmohammadzia and Koningsveld (2024)). Entities are clustered into \(n\) groups (e.g. \(n=15\) for the small instance). For each cluster, the centroid is taken as the demand zone locationa and the sum of demands of all entities in that cluster gives the zone’s total weekly demand (\(d_i\)).

#

Candidate Transfer Points (TPs)#

The Municipality of Amsterdam (MoA) has identified 36 potential locations for HoReCa Transfer Points (TPs) along the canals, and classified them into three categories: Spacious, Moderate, and Poor. For the optimisation experiments, we select:

  • Small instance: 10 TPs (a subset of the 36, spatially distributed)

  • Large instance: 36 TPs

For each TP, we define:

  • Its geographic location (longitude, latitude)

  • A capacity (\(cap_j\)) (maximum weekly throughput), which depends on its category

  • A fixed establishment cost (\(f_j\)), which depends on its category

Distance-based feasibility and service limits#

We define a central consolidation hub (outside the historic centre) as the origin of all vessel flows. Using the coordinates of the hub, TPs, and demand zones, we compute:

  • \(t^{water}_{0j}\): water distance (in km) from the hub to TP \(j\)

  • \(t^{road}_{ij}\) : road distance (in km) from TP \(j\) to demand zone \(i\)

We also impose a maximum last-mile distance for LVs: \(D_{max} = 2\)

We encode service feasibility with: \( \delta_{ij} = \begin{cases} 1 & \text{if } t^{road}_{ij} \le D_{max} \\ 0 & \text{otherwise} \end{cases} \)

Global parameters#

In addition to zone and TP data, we use a few global parameters that are common to all instances.

Transport cost parameters#

We model cost per unit of demand per kilometre. These are derived from estimated cost-per-vehicle-km and typical vehicle capacities. For this assignment, we use \(c^{road} = 0.03\) and \(c^{water} = 0.01\) €/unit/km.

Emission parameters#

Similarly, we define emission factors per unit of demand per kilometre. For this assignment, we use \(e^{road} = 0.045\) and \(e^{water} = 0.015\) CO₂/unit/km.

These values are consistent with literature values for diesel trucks, LVs, and inland vessels, scaled down to “per unit” (pallet-equivalent).

Planning limits#

For each instance, we also specify Budget and Maximum number of TPs. Suggested values are:

  • Small instance: \(B = 1000\) and \(K = 10\)

  • Large instance: \(B = 2000\) and \(K = 20\)

Finally, we will experiment with different values of the emission penalty factor: \(\lambda \in \{1, 5, 10, 15\} \)

Task 1: (small instance)

Run the cells below to:

  • Import the required Python packages
  • Load the required data
  • Transform the data into the inputs required for the model

Note: You are not expected to fully understand all details of the data processing code at this stage. These cells make use of Pandas and DataFrames, which will be explained in more depth later in the course. For now, focus on how the processed data becomes structured inputs for the optimization.

# Task 1 – small instance
# Step 1: Import required packages

import math
import pandas as pd
import os
from urllib.request import urlretrieve

Note: If you do not have pandas, numpy, and math already installed, you can simply use: 

pip install YOUR_PACKAGE_NAME

# Task 1 – small instance
# Step 2: Load the small-instance data (demand zones + TPs)

def findfile(fname):
    if not os.path.isfile(fname):
        print(f"Downloading {fname}...")
        urlretrieve('http://files.mude.citg.tudelft.nl/'+fname, fname)

findfile("GA_2_5_data.xlsx")

data_file = "GA_2_5_data.xlsx"

# Read only the small sheets
df_sml_dem = pd.read_excel(data_file, sheet_name="Small_Demand")
df_sml_tp  = pd.read_excel(data_file, sheet_name="Small_TP")

# Standardise column names
df_sml_dem = df_sml_dem.rename(columns={
    "number": "zone_id",
    "latitude": "lat",
    "longitude": "lon",
    "demand": "demand"
})

df_sml_tp = df_sml_tp.rename(columns={
    "number": "tp_id",
    "latitude": "lat",
    "longitude": "lon",
    "capacity": "capacity",
    "establishment_cost": "fixed_cost"
})

print("Small demand zones:", df_sml_dem.shape[0])
print("Small TPs:", df_sml_tp.shape[0])

display(df_sml_dem.head())
display(df_sml_tp.head())

# Task 1 – small instance
# Step 3: Create index sets and basic parameters for the small instance

# Index sets
I_sml = df_sml_dem["zone_id"].tolist()  # demand zones
J_sml = df_sml_tp["tp_id"].tolist()     # transfer points

# Demand, capacity, fixed cost
d_sml   = dict(zip(df_sml_dem["zone_id"], df_sml_dem["demand"]))
cap_sml = dict(zip(df_sml_tp["tp_id"], df_sml_tp["capacity"]))
f_sml   = dict(zip(df_sml_tp["tp_id"], df_sml_tp["fixed_cost"]))

print("\nTotal demand (small):", sum(d_sml.values()))
print("Total TP capacity (small):", sum(cap_sml.values()))

# Task 1 – small instance
# Step 4: Define hub location and distance function

# Consolidation hub coordinates (given)
hub_lat = 52.405141
hub_lon = 4.872129

def haversine_km(lon1, lat1, lon2, lat2):
    """
    Compute great-circle distance between two points on Earth (km).
    Inputs and outputs are in decimal degrees and kilometers.
    """
    # convert degrees to radians
    lon1_rad, lat1_rad = math.radians(lon1), math.radians(lat1)
    lon2_rad, lat2_rad = math.radians(lon2), math.radians(lat2)

    dlon = lon2_rad - lon1_rad
    dlat = lat2_rad - lat1_rad

    a = math.sin(dlat/2.0)**2 + math.cos(lat1_rad)*math.cos(lat2_rad)*math.sin(dlon/2.0)**2
    c = 2 * math.asin(math.sqrt(a))
    R = 6371.0  # Earth radius in km
    return R * c

# Task 1 – small instance
# Step 5: Compute water distances t0j_water from hub to each TP (small instance)

t0j_water_sml = {}   # key = tp_id, value = distance (km)

for _, row in df_sml_tp.iterrows():
    j = row["tp_id"]
    dist_km = haversine_km(hub_lon, hub_lat, row["lon"], row["lat"])
    t0j_water_sml[j] = dist_km

print("Example t0j_water_sml (first 5 TPs):")
for j in list(t0j_water_sml.keys())[:5]:
    print(f"  TP {j}: {t0j_water_sml[j]:.2f} km from hub")

# Task 1 – small instance
# Step 6: Compute road distances tij_road between TPs and demand zones (small instance)

tij_road_sml = {}   # key = (i,j), value = distance (km)

for _, dz in df_sml_dem.iterrows():
    i = dz["zone_id"]
    for _, tp in df_sml_tp.iterrows():
        j = tp["tp_id"]
        dist_km = haversine_km(tp["lon"], tp["lat"], dz["lon"], dz["lat"])
        tij_road_sml[(i, j)] = dist_km

print("Example tij_road_sml (first 5 pairs):")
for key in list(tij_road_sml.keys())[:5]:
    i, j = key
    print(f"  Zone {i} – TP {j}: {tij_road_sml[key]:.2f} km")

# Task 1 – small instance
# Step 7: Define global cost/emission parameters and compute CT_ij and Em_ij

# Global transport cost parameters (€/unit/km)
c_road  = 0.03   # LEV last-mile
c_water = 0.01   # vessel hub -> TP

# Global emission parameters (kg CO2 / unit / km)
e_road  = 0.045
e_water = 0.015

# Maximum LV service distance (km) for small instance
D_max_sml = 2

CT_sml = {}   # transport cost per unit
Em_sml = {}   # emissions per unit

for i in I_sml:
    for j in J_sml:
        t_water = t0j_water_sml[j]
        t_road  = tij_road_sml[(i, j)]
        CT_sml[(i, j)] = c_water * t_water + c_road * t_road
        Em_sml[(i, j)] = e_water * t_water + e_road * t_road

print("Example CT_sml and Em_sml (first 5 pairs):")
for key in list(CT_sml.keys())[:5]:
    i, j = key
    print(f"  (i={i}, j={j})  CT={CT_sml[key]:.3f} €/unit  Em={Em_sml[key]:.3f} kg CO2/unit")

# Task 1 – small instance
# Step 8: Compute delta_ij (feasibility) and big-M for small instance

delta_sml = {}   # key = (i,j), value = 0 or 1

for i in I_sml:
    for j in J_sml:
        if tij_road_sml[(i, j)] <= D_max_sml:
            delta_sml[(i, j)] = 1
        else:
            delta_sml[(i, j)] = 0

# Big-M: upper bound on any y_ij. A safe choice is total demand.
M_sml = sum(d_sml.values())

print("D_max (small):", D_max_sml, "km")
print("Big-M (small):", M_sml)
print("Number of feasible (i,j) pairs (delta=1):",
      sum(delta_sml[(i, j)] for i in I_sml for j in J_sml))

Building the MILP model in Pyomo (small instance)#

You now have all the inputs for the small instance (10 TPs, 15 demand zones). Start by implementing the MILP model in Pyomo for \(\lambda = 1\):

Task 2: Implement the Pyomo model for the small instance (λ = 1)

Translate the mathematical formulation into Pyomo code:
  • Create a ConcreteModel for the small instance
  • Define sets, parameters, decision variables
  • Define the objective function (with λ = 1)
  • Add all constraints: demand, capacity, budget, maximum TPs, service feasibility
# Task 2: Build the Pyomo model for the small instance (λ = 1)

from pyomo.environ import (
    ConcreteModel, Set, Param, Var, Objective, Constraint,
    Binary, NonNegativeReals, minimize
)

# Planning parameters for the small instance
B = 1000   # budget
K = 10     # maximum number of TPs
Lambda = 1  # emission penalty factor 

# 2.1. Create the model object
ConcreteModel()

# 2.2. Define index sets I and J based on I_sml and J_sml
m.I = Set(initialize=I_sml)
m.J = Set(initialize=J_sml)

# 2.3. Define parameters using the dictionaries from Task 1: demand, capacity, fixed cost, transport cost, emissions, delta (0/1 feasibility), B and K
m.d     = Param(m.I, initialize=d_sml)
m.cap   = Param(m.J, initialize=cap_sml)
m.f     = Param(m.J, initialize=f_sml)
m.CT    = Param(m.I, m.J, initialize=CT_sml)
m.Em    = Param(m.I, m.J, initialize=Em_sml)
m.delta = Param(m.I, m.J, initialize=delta_sml)

m.B = B
m.K = K
m.M = M_sml
m.Lambda = Lambda

# 2.4: Define decision variables:
# Hint: You can add decision variables in the following form
#   m.x = Var(m.J, domain=Binary)
### YOUR CODE HERE ###

# 2.5: Define the objective function for λ = 1:
### YOUR CODE HERE ###

# 2.6: Add constraints: (Demand fulfillment, TP capacity, Budget constraint, Maximum number of TPs, Service distance feasibility (big-M))
# Hint: You can define each constraint set as a separate function + Constraint, e.g.
#   def demand_rule(m, i):
#       return sum(m.y[i,j] for j in m.J) == m.d[i]
#   m.demand = Constraint(m.I, rule=demand_rule)

# Demand fulfillment
### YOUR CODE HERE ###

# TP capacity
### YOUR CODE HERE ###

# Budget
### YOUR CODE HERE ###

# Maximum number of TPs
### YOUR CODE HERE ###

# Service feasibility (big-M)
### YOUR CODE HERE ###

Solving the MILP (small instance, λ = 1)#

You now have:

  • a fully defined Pyomo model m for the small instance, and

  • all parameters loaded from the processed data.

Task 3: Solve and analyse the small-instance solution (λ = 1)

Do the following:
  • Call the HiGHS solver via SolverFactory("highs") and solve the model m
  • Check and print the solver status and termination condition
  • Identify which TPs are opened (values of x[j])
  • For each open TP, compute how much demand it serves (sum of y[i,j] over all zones)
  • Verify that: all demand is satisfied and no TP exceeds its capacity
  • Compute and report: total fixed cost, total fixed cost, total emissions, total objective value
# Task 3: Solve and analyse the small-instance model (λ = 1)

from pyomo.environ import SolverFactory, value
import time

start_time = time.time()
# 3.1. Solve the model m, with tee=True
solver = SolverFactory("highs")
results = solver.solve(m, tee=True)

solve_time_small = time.time() - start_time
print(f"\nSolve time (small instance, λ = 1): {solve_time_small:.3f} seconds")

# 3.2. Print solver status and termination condition
print(results.solver.status)
print(results.solver.termination_condition)

# 3.3. Identify which TPs are opened
open_tps_small = [j for j in m.J if value(m.x[j]) > 0.5]
print("\nOpen TPs in small instance:", open_tps_small)

# 3.4. For each open TP j, compute total flow served: sum_i y[i,j] and compare it to cap[j]
print("\nFlow through each open TP:")
for j in open_tps_small:
    flow_j = sum(value(m.y[i, j]) for i in m.I)
    print(f"  TP {j}: total flow = {flow_j:.2f} units")

# 3.5. Check demand satisfaction for each zone:
print("\nDemand satisfaction per zone:")
for i in m.I:
    served_i = sum(value(m.y[i, j]) for j in m.J)
    print(f"  Zone {i}: demand = {value(m.d[i]):.2f}, served = {served_i:.2f}")

# 3.6. Compute cost components and emissions:
total_fixed_cost = sum(value(m.f[j]) * value(m.x[j]) for j in m.J)
total_transport_cost = sum(value(m.CT[i, j]) * value(m.y[i, j]) for i in m.I for j in m.J)
total_emissions = sum(value(m.Em[i, j]) * value(m.y[i, j]) for i in m.I for j in m.J)
total_objective = value(m.obj)

print("\nCost and emissions summary (λ = 1):")
print(f"  Total fixed cost      = {total_fixed_cost:.2f}")
print(f"  Total transport cost  = {total_transport_cost:.2f}")
print(f"  Total emissions       = {total_emissions:.2f} kg CO2")
print(f"  Total emission penalty = {(total_emissions*Lambda):.2f}")
print(f"  Objective value Z     = {total_objective:.2f}")

Adding emission penalty and exploring different λ values#

So far, you solved the small-instance model with \(\lambda=1\). For larger \(\lambda\), emitting more becomes more expensive, which may push the model to: change which TPs are opened, change how demand is assigned, and accept higher monetary costs in exchange for lower emissions.

In this task, you will solve the same small-instance problem for several \(\lambda\) values and compare the results.

Task 4: Experiment with different λ values (small instance)

Do the following:
  • Choose the set of λ values, e.g. [1, 5, 10, 15]
  • For each λ:
    • build a model with that λ
    • solve it with HiGHS
    • record:
      • which TPs are opened
      • total fixed cost
      • total transport cost
      • total emissions and emission penalty
      • objective value
  • Compare how the solution changes when λ increases.
Hint: To keep your code organised, it is convenient to write a small function that builds a model for a given λ. 
lambda_values = [1, 5, 10, 15]
results_lambda = []

def build_model(lambda_value):
    """
    Build and return a Pyomo model m for a given λ,
    using the small-instance data prepared in Task 1.
    """
    m = ConcreteModel()

    # Sets
    m.I = Set(initialize=I_sml)
    m.J = Set(initialize=J_sml)

    # Parameters
    m.d     = Param(m.I, initialize=d_sml)
    m.cap   = Param(m.J, initialize=cap_sml)
    m.f     = Param(m.J, initialize=f_sml)
    m.CT    = Param(m.I, m.J, initialize=CT_sml)
    m.Em    = Param(m.I, m.J, initialize=Em_sml)
    m.delta = Param(m.I, m.J, initialize=delta_sml)

    # Planning parameters (use same B, K as before)
    m.B = B
    m.K = K
    m.M = M_sml
    m.lambda_value = lambda_value

    # Decision variables
    ### YOUR CODE HERE ###

    # Objective
    ### YOUR CODE HERE ###

    # Constraints
    ### YOUR CODE HERE ###

    # Demand fulfillment
    ### YOUR CODE HERE ###

    # TP capacity
    ### YOUR CODE HERE ###

    # Budget
    ### YOUR CODE HERE ###

    # Maximum number of TPs
    ### YOUR CODE HERE ###

    # Service feasibility (big-M)
    ### YOUR CODE HERE ###

    return m

# Loop over λ values
solver = SolverFactory("highs")

for lam in lambda_values:
    m_lam = build_model(lam)
    res = solver.solve(m_lam, tee=False)

    # Identify open TPs
    open_tps = [j for j in m_lam.J if value(m_lam.x[j]) > 0.5]

    # Cost components
    total_fixed_cost = sum(value(m_lam.f[j]) * value(m_lam.x[j]) for j in m_lam.J)
    total_transport_cost = sum(value(m_lam.CT[i, j]) * value(m_lam.y[i, j]) for i in m_lam.I for j in m_lam.J)
    total_emissions = sum(value(m_lam.Em[i, j]) * value(m_lam.y[i, j]) for i in m_lam.I for j in m_lam.J)
    total_obj = value(m_lam.obj)

    
 
    results_lambda.append({
        "lambda": lam,
        "open_tps": open_tps,
        "fixed_cost": total_fixed_cost,
        "transport_cost": total_transport_cost,
        "emissions": total_emissions,
        "objective": total_obj
    }) 

# Print summary of results
for r in results_lambda:
    print(f"λ = {r['lambda']}:")
    print(f"  Open TPs         = {r['open_tps']}")
    print(f"  Fixed cost       = {r['fixed_cost']:.2f}")
    print(f"  Transport cost   = {r['transport_cost']:.2f}")
    print(f"  Emissions        = {r['emissions']:.2f} kg CO2") 
    print(f"  Emission penalty = {r['emissions']*r['lambda']:.2f}")     
    print(f"  Objective value  = {r['objective']:.2f}\n")

Scaling up: large instance and runtime comparison#

So far, you have built and solved the small instance (10 TPs, 15 demand zones) and recorded its solve time.

We now consider a larger instance with: 36 candidate TPs and 300 HoReCa demand zones. The model structure is exactly the same, but the number of decision variables and constraints increases and the MILP becomes harder to solve.

The main goal of this task is to see how the solver runtime changes when the problem size increases. We will not re-experiment with different λ values here; we use a single λ = 1.

Task 5: Solve the large instance and compare runtimes

Do the following:
  • Load and compute the large-instance data (similar to task 1)
  • Build a new Pyomo model for the large instance:
    • Use the same structure as before (same component names m.I, m.J, m.x, m.y, etc.)
    • Use a single λ value (λ = 1)
  • Solve the large-instance model with HiGHS and record the solve time
  • Compare:
    • small vs large solve time
    • How many TPs are opened in each case
    • Any noticeable differences in model behaviour
# Task 5 – large instance
# Step 1: Load the large-instance data (demand zones + TPs)

data_file = "GA_2_5_data.xlsx"

df_lg_dem = pd.read_excel(data_file, sheet_name="Large_Demand")
df_lg_tp  = pd.read_excel(data_file, sheet_name="Large_TP")

df_lg_dem = df_lg_dem.rename(columns={
    "number": "zone_id",
    "latitude": "lat",
    "longitude": "lon",
    "demand": "demand"
})

df_lg_tp = df_lg_tp.rename(columns={
    "number": "tp_id",
    "latitude": "lat",
    "longitude": "lon",
    "capacity": "capacity",
    "establishment_cost": "fixed_cost"
})

I_large = df_lg_dem["zone_id"].tolist()
J_large = df_lg_tp["tp_id"].tolist()

d_large   = dict(zip(df_lg_dem["zone_id"], df_lg_dem["demand"]))
cap_large = dict(zip(df_lg_tp["tp_id"], df_lg_tp["capacity"]))
f_large   = dict(zip(df_lg_tp["tp_id"], df_lg_tp["fixed_cost"]))

print("Large instance – |I| =", len(I_large), ", |J| =", len(J_large))

# Task 5 – large instance
# Step 2: Compute distances (water + road), CT_ij, Em_ij, delta_ij, M for large instance

# Water distances hub -> TP
t0j_water_large = {}
for _, row in df_lg_tp.iterrows():
    j = row["tp_id"]
    dist_km = haversine_km(hub_lon, hub_lat, row["lon"], row["lat"])
    t0j_water_large[j] = dist_km

# Road distances TP -> demand zone
tij_road_large = {}
for _, dz in df_lg_dem.iterrows():
    i = dz["zone_id"]
    for _, tp in df_lg_tp.iterrows():
        j = tp["tp_id"]
        dist_km = haversine_km(tp["lon"], tp["lat"], dz["lon"], dz["lat"])
        tij_road_large[(i, j)] = dist_km

        
D_max_large = 2

CT_large = {}
Em_large = {}
delta_large = {}

for i in I_large:
    for j in J_large:
        t_water = t0j_water_large[j]
        t_road  = tij_road_large[(i, j)]
        CT_large[(i, j)] = c_water * t_water + c_road * t_road
        Em_large[(i, j)] = e_water * t_water + e_road * t_road
        delta_large[(i, j)] = 1 if t_road <= D_max_large else 0

M_large = sum(d_large.values())

print("D_max (large):", D_max_large, "km")
print("Big-M (large):", M_large)
print("Feasible (i,j) pairs (delta=1):",
      sum(delta_large[(i, j)] for i in I_large for j in J_large))

# Task 5 – large instance
# Step 3: Build and solve the large-instance model

B_large = 2000   
K_large = 20 
lambda_large = 1


m_large = ConcreteModel()

# Sets
m_large.I = Set(initialize=I_large)
m_large.J = Set(initialize=J_large)

# Parameters
m_large.d     = Param(m_large.I, initialize=d_large)
m_large.cap   = Param(m_large.J, initialize=cap_large)
m_large.f     = Param(m_large.J, initialize=f_large)
m_large.CT    = Param(m_large.I, m_large.J, initialize=CT_large)
m_large.Em    = Param(m_large.I, m_large.J, initialize=Em_large)
m_large.delta = Param(m_large.I, m_large.J, initialize=delta_large)

m_large.B = B_large
m_large.K = K_large
m_large.M = M_large

# Variables
### YOUR CODE HERE ###

# Objective
### YOUR CODE HERE ###

# Constraints
### YOUR CODE HERE ###

start_time_large = time.time()
# Solve and measure runtime
solver_large = SolverFactory("highs")
results_large = solver_large.solve(m_large, tee=False)

solve_time_large = time.time() - start_time_large
print(f"\nSolve time (large instance, λ = {lambda_large}): {solve_time_large:.3f} seconds")

# Task 5 – large instance
# Step 4: Basic comparison with small instance

open_tps_large = [j for j in m_large.J if value(m_large.x[j]) > 0.5]
print("\nOpen TPs in large instance:", open_tps_large)

total_obj_large = value(m_large.obj)
print(f"Objective value (large instance) = {total_obj_large:.2f}")

print("\nRuntime comparison:")
print(f"  Small instance solve time: {solve_time_small:.3f} seconds")
print(f"  Large instance solve time: {solve_time_large:.3f} seconds")

By Nadia Pourmohammadzia and Gonçalo Homen de Almeida Correia, Delft University of Technology. CC BY 4.0, more info on the Credits page of Workbook

Contents