diff --git a/cflp_cuopt_milp/README.md b/cflp_cuopt_milp/README.md new file mode 100644 index 0000000..efa189a --- /dev/null +++ b/cflp_cuopt_milp/README.md @@ -0,0 +1,107 @@ +# Capacitated Facility Location Optimization + +A Mixed Integer Linear Program (MILP) example for optimizing distribution center (DC) placement and customer assignment to minimize total logistics cost. + +## Problem Overview + +A logistics company needs to determine: +- Which distribution centers (DCs) to open from a set of candidate locations +- How to assign customers to open DCs +- How to balance fixed operating costs with transportation costs + +**Goal**: Minimize total annualized logistics cost while meeting demand and respecting capacity constraints + +## Model Assumptions + +**Important**: This model uses the following assumptions: + +1. **Single assignment**: Each customer is assigned to exactly one DC +2. **Capacity limits**: Each DC has a maximum pallet-handling capacity +3. **Fixed costs**: Opening a DC incurs a fixed annual operating cost +4. **Known demand**: Customer demand is deterministic and known +5. **Euclidean distances**: Transportation costs are proportional to straight-line distance + +This formulation works well for: +- Distribution network design +- Warehouse location planning +- Supply chain optimization +- Retail store placement + +**Note**: Extensions include multi-period capacity expansion and fractional (LP relaxed) assignments. + +## Notebook Contents + +### Setup & Data +- 5 candidate distribution centers +- 20 customers with varying demand +- Synthetic 2D coordinates for visualization +- Transportation cost proportional to distance ($0.05/pallet-km) +- Fixed DC operating costs (~$80,000-$120,000) + +### Optimization +- **Decision Variables**: + - DC open/close binaries (y_i) + - Customer assignment binaries (x_ij) +- **Objective**: Minimize fixed costs + transportation costs +- **Constraints**: + - Assignment: Each customer assigned to exactly one DC + - Capacity: DC load cannot exceed capacity + - Linking: Customers can only be assigned to open DCs + +### Results & Analysis +- Optimal DC selection and network design +- Customer-to-DC assignment matrix +- DC utilization rates +- Cost breakdown (fixed vs. transportation) +- Network visualization with assignment arcs + +### Extensions +- **LP Relaxation**: Fractional assignments for lower bound analysis +- **Multi-Period Expansion**: Time-phased DC opening decisions with discounting + +## Installation + +**Requirements**: +- NVIDIA GPU with CUDA 12 or 13 support +- Python 3.9+ + +**Install cuOpt** (choose one based on your CUDA version): + +```bash +# For CUDA 12 +pip install --upgrade --extra-index-url=https://pypi.nvidia.com cuopt-cu12 + +# For CUDA 13 +pip install --upgrade --extra-index-url=https://pypi.nvidia.com cuopt-cu13 +``` + +**Install visualization dependencies**: + +```bash +pip install matplotlib seaborn +``` + +## Quick Start + +```bash +jupyter notebook cflp_cuopt_milp.ipynb +``` + +The notebook includes GPU detection and will guide you through any missing dependencies. + +## Possible Extensions + +**Multi-Period Planning** (included in notebook): +- Open DCs over multiple time periods +- Time-discounted costs +- Time-varying demand with growth rates +- Capacity expansion decisions + +**Additional**: +- Multiple DC capacity tiers +- Product-specific assignment +- Stochastic demand scenarios +- Service level constraints (max distance) +- Multi-echelon supply chain +- Inventory considerations + diff --git a/cflp_cuopt_milp/cflp_cuopt_milp.ipynb b/cflp_cuopt_milp/cflp_cuopt_milp.ipynb new file mode 100644 index 0000000..73ad121 --- /dev/null +++ b/cflp_cuopt_milp/cflp_cuopt_milp.ipynb @@ -0,0 +1,1366 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Capacitated Facility Location with NVIDIA cuOpt (MILP)\n", + "\n", + "This notebook formulates and solves a capacitated facility location problem (CFLP) using the NVIDIA cuOpt MILP Python API.\n", + "\n", + "## Problem Description\n", + "\n", + "A logistics company operates a distribution network. They need to determine:\n", + "- Which distribution centers (DCs) to open from a set of candidate locations\n", + "- How to assign customers to open DCs to fulfill their demand\n", + "\n", + "The goal is to minimize total annual logistics cost while:\n", + "- Meeting all customer demand\n", + "- Respecting DC capacity constraints\n", + "- Balancing fixed operating costs with transportation costs\n", + "\n", + "### Model Assumptions\n", + "\n", + "**Important**: This model makes the following simplifying assumptions:\n", + "1. **Single assignment**: Each customer is assigned to exactly one DC\n", + "2. **Capacity limits**: Each DC has a maximum pallet-handling capacity per week\n", + "3. **Fixed costs**: Opening a DC incurs a fixed annual operating cost\n", + "4. **Known demand**: Customer demand is deterministic and known in advance\n", + "5. **Euclidean distances**: Transportation costs are proportional to straight-line distance\n", + "\n", + "This formulation is well-suited for:\n", + "- Distribution network design\n", + "- Warehouse location planning\n", + "- Supply chain optimization\n", + "- Retail store placement decisions\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Environment Setup\n", + "\n", + "First, let's check if we have a GPU available and install necessary dependencies.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import subprocess\n", + "import html\n", + "from IPython.display import display, HTML\n", + "\n", + "def check_gpu():\n", + " try:\n", + " result = subprocess.run([\"nvidia-smi\"], capture_output=True, text=True, timeout=5)\n", + " result.check_returncode()\n", + " lines = result.stdout.splitlines()\n", + " gpu_info = lines[2] if len(lines) > 2 else \"GPU detected\"\n", + " gpu_info_escaped = html.escape(gpu_info)\n", + " display(HTML(f\"\"\"\n", + "
\n", + "

✅ GPU is enabled

\n", + " 
{gpu_info_escaped}
\n", + "
\n", + " \"\"\"))\n", + " return True\n", + " except (subprocess.CalledProcessError, subprocess.TimeoutExpired, FileNotFoundError, IndexError) as e:\n", + " display(HTML(\"\"\"\n", + "
\n", + "

⚠️ GPU not detected!

\n", + "

This notebook requires a GPU runtime.

\n", + " \n", + "

If running in Google Colab:

\n", + "
    \n", + "
  1. Click on Runtime → Change runtime type
    2. \n", + "
  2. Set Hardware accelerator to GPU
    3. \n", + "
  3. Then click Save and Runtime → Restart runtime.
    4. \n", + "
\n", + " \n", + "

If running in Docker:

\n", + "
    \n", + "
  1. Ensure you have NVIDIA Docker runtime installed (nvidia-docker2)
    2. \n", + "
  2. Run container with GPU support: docker run --gpus all ...
    3. \n", + "
  3. Or use: docker run --runtime=nvidia ... for older Docker versions
    4. \n", + "
  4. Verify GPU access: docker run --gpus all nvidia/cuda:12.0.0-base-ubuntu22.04 nvidia-smi
    5. \n", + "
\n", + " \n", + "

Additional resources:

\n", + " \n", + "
\n", + " \"\"\"))\n", + " return False\n", + "\n", + "check_gpu()\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Check if cuOpt is installed, if not provide installation instructions\n", + "try:\n", + " import cuopt\n", + " print(\"✓ cuOpt is installed\")\n", + "except ImportError:\n", + " print(\"⚠️ cuOpt is not installed!\")\n", + " print(\"\\nTo install cuOpt, uncomment and run ONE of the following commands:\")\n", + " print(\" For CUDA 12: %pip install --upgrade --extra-index-url=https://pypi.nvidia.com cuopt-cu12\")\n", + " print(\" For CUDA 13: %pip install --upgrade --extra-index-url=https://pypi.nvidia.com cuopt-cu13\")\n", + " print(\"\\nThen restart the kernel and run again.\")\n", + " raise ImportError(\"cuOpt is required. Please install it using the instructions above.\")\n", + "\n", + "# Uncomment ONE of the following lines to install cuOpt (requires GPU with CUDA):\n", + "# %pip install --upgrade --extra-index-url=https://pypi.nvidia.com cuopt-cu12 nvidia-nvjitlink-cu12\n", + "# %pip install --upgrade --extra-index-url=https://pypi.nvidia.com cuopt-cu13 nvidia-nvjitlink-cu13\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%pip install matplotlib seaborn\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Import Required Libraries\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import pandas as pd\n", + "import time\n", + "import matplotlib.pyplot as plt\n", + "import seaborn as sns\n", + "\n", + "from cuopt.linear_programming.problem import Problem, VType, sense, LinearExpression\n", + "from cuopt.linear_programming.solver_settings import SolverSettings\n", + "\n", + "# Set style for better visualizations\n", + "sns.set_style(\"whitegrid\")\n", + "plt.rcParams['figure.figsize'] = (12, 6)\n", + "\n", + "print(\"✓ All imports successful\")\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Problem Statement\n", + "\n", + "- We have `num_dcs` candidate distribution centers (DCs) with weekly pallet-handling capacities and fixed operating costs.\n", + "- We have `num_customers` customers, each with a weekly pallet demand and a location.\n", + "- Transportation cost is proportional to distance between DC and customer.\n", + "- We must:\n", + " - Decide which DCs to open.\n", + " - Assign every customer's demand to exactly one open DC.\n", + "- Objective: minimize total annualized logistics cost = fixed DC cost + transportation cost.\n", + "- Constraints:\n", + " - Each customer's demand is fully assigned to a single DC.\n", + " - The sum of assigned demand at each DC does not exceed that DC's capacity.\n", + " - No customer can be assigned to a closed DC.\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Problem Data Setup\n", + "\n", + "Generate a synthetic instance with candidate DC locations and customer locations.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# -----------------------------\n", + "# Synthetic instance generation\n", + "# -----------------------------\n", + "\n", + "rng = np.random.default_rng(42)\n", + "\n", + "# Sets\n", + "num_dcs = 5\n", + "num_customers = 20\n", + "\n", + "I = list(range(num_dcs))\n", + "J = list(range(num_customers))\n", + "\n", + "# DC locations on a 2D grid (for visualization)\n", + "dc_coords = rng.uniform(low=0.0, high=1000.0, size=(num_dcs, 2))\n", + "\n", + "# Customer locations\n", + "cust_coords = rng.uniform(low=0.0, high=1000.0, size=(num_customers, 2))\n", + "\n", + "# Euclidean distances and cost per pallet\n", + "dist_matrix = np.linalg.norm(\n", + " dc_coords[:, None, :] - cust_coords[None, :, :],\n", + " axis=2\n", + ")\n", + "\n", + "alpha = 0.05 # $ per pallet-km\n", + "unit_cost = alpha * dist_matrix\n", + "\n", + "# Weekly customer demand (pallets)\n", + "demand = rng.integers(low=80, high=250, size=num_customers).astype(float)\n", + "\n", + "# DC capacity: scaled so that total capacity > total demand\n", + "total_demand = demand.sum()\n", + "avg_capacity = total_demand * 1.4 / num_dcs\n", + "capacity = avg_capacity * rng.uniform(low=0.8, high=1.2, size=num_dcs)\n", + "capacity = capacity.astype(float)\n", + "\n", + "# Fixed opening costs (proportional to capacity + noise)\n", + "base_fc = 80000.0\n", + "fixed_cost = base_fc + 20.0 * capacity + rng.normal(loc=0.0, scale=5000.0, size=num_dcs)\n", + "fixed_cost = np.maximum(fixed_cost, 0.0).astype(float) # Ensure non-negative costs\n", + "\n", + "# Allowed arcs: all DC–customer pairs\n", + "A = [(i, j) for i in I for j in J]\n", + "\n", + "print(f\"Problem Size: {num_dcs} DCs, {num_customers} customers\")\n", + "print(f\"Total demand: {total_demand:.0f} pallets/week\")\n", + "print(f\"\\nDC capacities: {capacity.round(0)}\")\n", + "print(f\"Fixed costs: ${fixed_cost.round(0)}\")\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Visualize problem data\n", + "fig, axes = plt.subplots(1, 2, figsize=(15, 5))\n", + "\n", + "# Plot 1: DC and Customer locations\n", + "ax1 = axes[0]\n", + "ax1.scatter(cust_coords[:, 0], cust_coords[:, 1], c=\"tab:blue\", s=demand/2, alpha=0.6, label=\"Customers\")\n", + "ax1.scatter(dc_coords[:, 0], dc_coords[:, 1], c=\"tab:red\", marker=\"s\", s=100, label=\"DCs\")\n", + "\n", + "for i, (xi, yi) in enumerate(dc_coords):\n", + " ax1.text(xi + 15, yi + 15, f\"DC {i}\", fontsize=9, color=\"red\")\n", + "\n", + "ax1.set_xlabel(\"X coordinate\")\n", + "ax1.set_ylabel(\"Y coordinate\")\n", + "ax1.set_title(\"DC and Customer Locations\\n(Customer size = demand)\", fontsize=14, fontweight=\"bold\")\n", + "ax1.legend()\n", + "ax1.grid(True, alpha=0.3)\n", + "\n", + "# Plot 2: DC capacity and fixed cost\n", + "ax2 = axes[1]\n", + "x_pos = np.arange(num_dcs)\n", + "width = 0.35\n", + "\n", + "bars1 = ax2.bar(x_pos - width/2, capacity, width, label='Capacity (pallets)', color='tab:blue', alpha=0.7)\n", + "ax2_twin = ax2.twinx()\n", + "bars2 = ax2_twin.bar(x_pos + width/2, fixed_cost/1000, width, label='Fixed Cost ($K)', color='tab:orange', alpha=0.7)\n", + "\n", + "ax2.set_xlabel('DC')\n", + "ax2.set_ylabel('Capacity (pallets/week)', color='tab:blue')\n", + "ax2_twin.set_ylabel('Fixed Cost ($K)', color='tab:orange')\n", + "ax2.set_title('DC Capacity and Fixed Costs', fontsize=14, fontweight='bold')\n", + "ax2.set_xticks(x_pos)\n", + "ax2.set_xticklabels([f'DC {i}' for i in I])\n", + "ax2.legend(loc='upper left')\n", + "ax2_twin.legend(loc='upper right')\n", + "\n", + "plt.tight_layout()\n", + "plt.show()\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Problem Formulation\n", + "\n", + "### Decision Variables\n", + "- `y[i]`: Binary, 1 if DC i is open, 0 otherwise\n", + "- `x[i,j]`: Binary, 1 if customer j is assigned to DC i, 0 otherwise\n", + "\n", + "### Objective Function\n", + "Minimize: Fixed DC costs + Transportation costs\n", + "\n", + "$$\\min \\sum_i f_i y_i + \\sum_{i,j} c_{ij} d_j x_{ij}$$\n", + "\n", + "### Constraints\n", + "1. **Assignment**: Each customer assigned to exactly one DC: $\\sum_i x_{ij} = 1$ for each customer $j$\n", + "2. **Capacity/Linking**: DC load cannot exceed capacity: $\\sum_j d_j x_{ij} \\le u_i y_i$ for each DC $i$\n", + "3. **Binary**: $x_{ij} \\in \\{0,1\\}$, $y_i \\in \\{0,1\\}$\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# -----------------------------\n", + "# Problem and variable creation\n", + "# -----------------------------\n", + "\n", + "prob = Problem(\"CFLP\")\n", + "\n", + "# Decision variables\n", + "y = {} # DC open/close binaries\n", + "x = {} # assignment binaries\n", + "\n", + "# y_i: 1 if DC i is open\n", + "for i in I:\n", + " y[i] = prob.addVariable(\n", + " name=f\"y_{i}\",\n", + " lb=0.0,\n", + " ub=1.0,\n", + " vtype=VType.INTEGER\n", + " )\n", + "\n", + "# x_ij: 1 if customer j is assigned to DC i\n", + "for (i, j) in A:\n", + " x[i, j] = prob.addVariable(\n", + " name=f\"x_{i}_{j}\",\n", + " lb=0.0,\n", + " ub=1.0,\n", + " vtype=VType.INTEGER\n", + " )\n", + "\n", + "print(f\"Number of variables: {prob.NumVariables}\")\n", + "print(f\" - DC open/close (y): {num_dcs}\")\n", + "print(f\" - Assignment (x): {num_dcs * num_customers}\")\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# -----------------------------\n", + "# Objective function\n", + "# -----------------------------\n", + "\n", + "obj = LinearExpression([], [], 0.0)\n", + "\n", + "# Fixed facility costs\n", + "for i in I:\n", + " obj += float(fixed_cost[i]) * y[i]\n", + "\n", + "# Transportation costs: c_ij * d_j * x_ij\n", + "for (i, j) in A:\n", + " flow_cost = float(unit_cost[i, j] * demand[j])\n", + " obj += flow_cost * x[i, j]\n", + "\n", + "prob.setObjective(obj, sense=sense.MINIMIZE)\n", + "\n", + "print(\"Objective set: Minimize (Fixed DC Costs + Transportation Costs)\")\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# -----------------------------\n", + "# Constraints\n", + "# -----------------------------\n", + "\n", + "# 1) Each customer assigned to exactly one DC\n", + "for j in J:\n", + " expr = LinearExpression([], [], 0.0)\n", + " for i in I:\n", + " expr += x[i, j]\n", + " prob.addConstraint(expr == 1.0, name=f\"assign_{j}\")\n", + "\n", + "print(f\"Assignment constraints added: {len(J)}\")\n", + "\n", + "# 2) DC capacity and logical linking: sum_j d_j x_ij <= capacity[i] * y_i\n", + "for i in I:\n", + " expr = LinearExpression([], [], 0.0)\n", + " for j in J:\n", + " expr += float(demand[j]) * x[i, j]\n", + " expr -= float(capacity[i]) * y[i]\n", + " prob.addConstraint(expr <= 0.0, name=f\"capacity_{i}\")\n", + "\n", + "print(f\"Capacity constraints added: {len(I)}\")\n", + "print(f\"\\nTotal constraints: {prob.NumConstraints}\")\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Solve the MILP with cuOpt\n", + "\n", + "We now configure a few basic MILP parameters and call `solve()`.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# -----------------------------\n", + "# Solve\n", + "# -----------------------------\n", + "\n", + "settings = SolverSettings()\n", + "\n", + "# Example tuning parameters (adjust according to your environment)\n", + "settings.set_parameter(\"time_limit\", \"300.0\") # seconds\n", + "settings.set_parameter(\"mip_relative_gap\", \"1e-4\") # relative optimality gap\n", + "\n", + "print(\"Starting optimization...\")\n", + "print(f\"Problem size: {prob.NumVariables} variables, {prob.NumConstraints} constraints\")\n", + "print()\n", + "\n", + "start_time = time.time()\n", + "prob.solve(settings=settings)\n", + "solve_time = time.time() - start_time\n", + "\n", + "print(f\"\\nOptimization completed in {solve_time:.3f} seconds\")\n", + "print(f\"Solve status: {prob.Status}\")\n", + "print(f\"Best objective value: ${prob.ObjValue:,.2f}\")\n", + "print(f\"Solve time (internal): {prob.SolveTime:.3f}s\")\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Extract and Analyze the Solution\n", + "\n", + "We inspect which DCs are opened, how customers are assigned, and the resulting loads on each DC.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# -----------------------------\n", + "# Solution extraction\n", + "# -----------------------------\n", + "\n", + "# Status codes: 1 = Optimal, 2 = Feasible, others indicate issues\n", + "# See cuOpt documentation for full status code reference\n", + "if prob.Status not in (1, 2):\n", + " status_msg = {0: \"Unknown\", 3: \"Infeasible\", 4: \"Unbounded\", 5: \"Time Limit\"}.get(prob.Status, \"Unknown\")\n", + " raise RuntimeError(f\"Solver did not find a solution: Status = {prob.Status} ({status_msg})\")\n", + "\n", + "# Open DCs\n", + "open_dcs = [i for i in I if y[i].getValue() > 0.5]\n", + "print(\"Open DCs:\", open_dcs)\n", + "\n", + "# Assignment matrix\n", + "assign_matrix = np.zeros((num_dcs, num_customers))\n", + "for (i, j) in A:\n", + " val = x[i, j].getValue()\n", + " assign_matrix[i, j] = val\n", + "\n", + "# For each customer j, identify assigned DC i\n", + "assigned_dc = assign_matrix.argmax(axis=0)\n", + "\n", + "# Load per DC\n", + "dc_load = assign_matrix @ demand\n", + "\n", + "print(f\"\\nDC loads: {dc_load.round(0)}\")\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Create summary DataFrames\n", + "dc_df = pd.DataFrame({\n", + " \"DC\": I,\n", + " \"Open\": [\"Yes\" if i in open_dcs else \"No\" for i in I],\n", + " \"Capacity\": capacity.round(0),\n", + " \"Load\": dc_load.round(0),\n", + " \"Utilization %\": (dc_load / capacity * 100).round(1),\n", + " \"Fixed Cost ($)\": fixed_cost.round(0)\n", + "})\n", + "\n", + "cust_df = pd.DataFrame({\n", + " \"Customer\": J,\n", + " \"Demand\": demand.round(0),\n", + " \"Assigned DC\": assigned_dc,\n", + " \"Transport Cost ($)\": [unit_cost[assigned_dc[j], j] * demand[j] for j in J]\n", + "})\n", + "\n", + "print(\"\\n\" + \"=\"*70)\n", + "print(\" \" * 20 + \"DC SUMMARY\")\n", + "print(\"=\"*70)\n", + "display(dc_df)\n", + "\n", + "print(\"\\n\" + \"=\"*70)\n", + "print(\" \" * 20 + \"CUSTOMER ASSIGNMENTS\")\n", + "print(\"=\"*70)\n", + "display(cust_df)\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Calculate cost breakdown\n", + "total_fixed_cost = sum(fixed_cost[i] for i in open_dcs)\n", + "total_transport_cost = sum(unit_cost[assigned_dc[j], j] * demand[j] for j in J)\n", + "\n", + "print(\"\\n\" + \"=\"*70)\n", + "print(\" \" * 20 + \"COST BREAKDOWN\")\n", + "print(\"=\"*70)\n", + "print(f\"\\n📊 Cost Summary:\")\n", + "print(f\" • Total Fixed Costs: ${total_fixed_cost:,.2f}\")\n", + "print(f\" • Total Transport Costs: ${total_transport_cost:,.2f}\")\n", + "print(f\" • Total Logistics Cost: ${prob.ObjValue:,.2f}\")\n", + "if prob.ObjValue > 0:\n", + " print(f\" • Fixed Cost %: {total_fixed_cost/prob.ObjValue*100:.1f}%\")\n", + " print(f\" • Transport Cost %: {total_transport_cost/prob.ObjValue*100:.1f}%\")\n", + "\n", + "print(f\"\\n🏭 Network Design:\")\n", + "print(f\" • DCs Opened: {len(open_dcs)} of {num_dcs}\")\n", + "print(f\" • Customers Served: {num_customers}\")\n", + "print(f\" • Total Demand Met: {total_demand:,.0f} pallets/week\")\n", + "\n", + "print(f\"\\n📦 Utilization:\")\n", + "if open_dcs:\n", + " avg_util = np.mean([dc_load[i]/capacity[i]*100 for i in open_dcs])\n", + " print(f\" • Average DC Utilization: {avg_util:.1f}%\")\n", + " for i in open_dcs:\n", + " print(f\" • DC {i}: {dc_load[i]/capacity[i]*100:.1f}% ({dc_load[i]:.0f}/{capacity[i]:.0f})\")\n", + "else:\n", + " print(\" • No DCs opened (unexpected - check solver status)\")\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Visualize the Network Design\n", + "\n", + "We now plot DCs, customers, and assignment arcs in a 2D plane (using the synthetic coordinates defined above).\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Create comprehensive visualizations\n", + "fig, axes = plt.subplots(2, 2, figsize=(16, 12))\n", + "\n", + "# Plot 1: Network design with assignments\n", + "ax1 = axes[0, 0]\n", + "\n", + "# Draw assignment lines\n", + "for j in J:\n", + " i = assigned_dc[j]\n", + " x1, y1 = dc_coords[i]\n", + " x2, y2 = cust_coords[j]\n", + " ax1.plot([x1, x2], [y1, y2],\n", + " color=\"tab:green\" if i in open_dcs else \"lightgray\",\n", + " linewidth=0.7, alpha=0.7)\n", + "\n", + "# Plot customers\n", + "ax1.scatter(cust_coords[:, 0], cust_coords[:, 1], c=\"tab:blue\", s=demand/2, alpha=0.6, label=\"Customers\")\n", + "for j, (xj, yj) in enumerate(cust_coords):\n", + " ax1.text(xj + 5, yj + 5, str(j), fontsize=8, color=\"blue\")\n", + "\n", + "# Plot DCs\n", + "colors = [\"tab:red\" if i in open_dcs else \"gray\" for i in I]\n", + "ax1.scatter(dc_coords[:, 0], dc_coords[:, 1], c=colors, marker=\"s\", s=80, label=\"DCs\")\n", + "\n", + "for i, (xi, yi) in enumerate(dc_coords):\n", + " label = f\"DC {i}\"\n", + " if i in open_dcs:\n", + " label += \" ✓\"\n", + " ax1.text(xi + 5, yi + 5, label, fontsize=8, color=\"red\" if i in open_dcs else \"gray\")\n", + "\n", + "ax1.set_xlabel(\"X coordinate\")\n", + "ax1.set_ylabel(\"Y coordinate\")\n", + "ax1.set_title(\"CFLP Network Design Solution\", fontsize=14, fontweight=\"bold\")\n", + "ax1.legend()\n", + "ax1.grid(True, alpha=0.3)\n", + "\n", + "# Plot 2: DC utilization\n", + "ax2 = axes[0, 1]\n", + "dc_colors = ['tab:green' if i in open_dcs else 'lightgray' for i in I]\n", + "bars = ax2.bar([f'DC {i}' for i in I], dc_load, color=dc_colors, alpha=0.7, label='Load')\n", + "ax2.bar([f'DC {i}' for i in I], capacity - dc_load, bottom=dc_load, color='white', edgecolor='gray', alpha=0.5, label='Unused')\n", + "\n", + "# Add capacity line markers\n", + "for i, bar in enumerate(bars):\n", + " ax2.hlines(capacity[i], bar.get_x(), bar.get_x() + bar.get_width(), colors='red', linestyles='--', linewidth=2)\n", + "\n", + "ax2.set_ylabel('Pallets/week')\n", + "ax2.set_title('DC Load vs Capacity', fontsize=14, fontweight='bold')\n", + "ax2.legend()\n", + "ax2.grid(True, alpha=0.3, axis='y')\n", + "\n", + "# Plot 3: Cost breakdown pie chart\n", + "ax3 = axes[1, 0]\n", + "cost_labels = ['Fixed Costs', 'Transport Costs']\n", + "cost_values = [total_fixed_cost, total_transport_cost]\n", + "colors_pie = ['tab:red', 'tab:blue']\n", + "ax3.pie(cost_values, labels=cost_labels, autopct='%1.1f%%', colors=colors_pie, startangle=90)\n", + "ax3.set_title('Cost Breakdown', fontsize=14, fontweight='bold')\n", + "\n", + "# Plot 4: Customers per DC\n", + "ax4 = axes[1, 1]\n", + "customers_per_dc = [sum(1 for j in J if assigned_dc[j] == i) for i in I]\n", + "demand_per_dc = dc_load\n", + "\n", + "x_pos = np.arange(num_dcs)\n", + "width = 0.35\n", + "bars1 = ax4.bar(x_pos - width/2, customers_per_dc, width, label='# Customers', color='tab:blue', alpha=0.7)\n", + "ax4_twin = ax4.twinx()\n", + "bars2 = ax4_twin.bar(x_pos + width/2, demand_per_dc, width, label='Demand (pallets)', color='tab:orange', alpha=0.7)\n", + "\n", + "ax4.set_xlabel('DC')\n", + "ax4.set_ylabel('Number of Customers', color='tab:blue')\n", + "ax4_twin.set_ylabel('Demand (pallets)', color='tab:orange')\n", + "ax4.set_title('Customer and Demand Distribution', fontsize=14, fontweight='bold')\n", + "ax4.set_xticks(x_pos)\n", + "ax4.set_xticklabels([f'DC {i}' for i in I])\n", + "ax4.legend(loc='upper left')\n", + "ax4_twin.legend(loc='upper right')\n", + "\n", + "plt.tight_layout()\n", + "plt.show()\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Extension: Semi-Relaxation (Fractional Customer Assignments)\n", + "\n", + "The following section demonstrates a **semi-relaxation** of the MILP:\n", + "- **DC open/close decisions (y)**: Remain INTEGER (binary)\n", + "- **Customer assignment (x)**: Changed to CONTINUOUS (fractional allowed)\n", + "\n", + "This allows customers to be partially assigned to multiple DCs, which can provide a lower bound on the optimal MILP solution or be used when fractional assignments are acceptable (e.g., splitting orders across DCs).\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# -----------------------------\n", + "# Semi-Relaxation: Fractional assignments allowed\n", + "# -----------------------------\n", + "\n", + "# Create a new problem for the semi-relaxed version\n", + "prob_lp = Problem(\"CFLP_SemiRelaxed\")\n", + "\n", + "# Decision variables\n", + "y_lp = {} # DC open/close - remain INTEGER (binary)\n", + "x_lp = {} # assignment - now CONTINUOUS (fractional allowed)\n", + "\n", + "# y_i: 1 if DC i is open (still binary/integer - facility decisions remain discrete)\n", + "for i in I:\n", + " y_lp[i] = prob_lp.addVariable(\n", + " name=f\"y_lp_{i}\",\n", + " lb=0.0,\n", + " ub=1.0,\n", + " vtype=VType.INTEGER # Still INTEGER (DC open/close decision)\n", + " )\n", + "\n", + "# x_ij: Fractional assignment (0 <= x_ij <= 1), now CONTINUOUS\n", + "for (i, j) in A:\n", + " x_lp[i, j] = prob_lp.addVariable(\n", + " name=f\"x_lp_{i}_{j}\",\n", + " lb=0.0,\n", + " ub=1.0,\n", + " vtype=VType.CONTINUOUS # Changed to CONTINUOUS for fractional assignments\n", + " )\n", + "\n", + "print(\"Semi-Relaxed Problem - Number of variables:\", prob_lp.NumVariables)\n", + "print(\"Note: y_i remain INTEGER (facility decisions), x_ij are CONTINUOUS (fractional assignments)\")\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# -----------------------------\n", + "# Objective function for relaxed LP\n", + "# -----------------------------\n", + "\n", + "obj_lp = LinearExpression([], [], 0.0)\n", + "\n", + "# Fixed facility costs (same as MILP)\n", + "for i in I:\n", + " obj_lp += float(fixed_cost[i]) * y_lp[i]\n", + "\n", + "# Transportation costs: c_ij * d_j * x_ij (same as MILP)\n", + "for (i, j) in A:\n", + " flow_cost = float(unit_cost[i, j] * demand[j])\n", + " obj_lp += flow_cost * x_lp[i, j]\n", + "\n", + "prob_lp.setObjective(obj_lp, sense=sense.MINIMIZE)\n", + "\n", + "print(\"Semi-relaxed objective set (minimize total cost).\")\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# -----------------------------\n", + "# Constraints for relaxed LP\n", + "# -----------------------------\n", + "\n", + "# 1) Each customer's demand must be fully assigned (sum of fractional assignments = 1)\n", + "for j in J:\n", + " expr_lp = LinearExpression([], [], 0.0)\n", + " for i in I:\n", + " expr_lp += x_lp[i, j]\n", + " prob_lp.addConstraint(expr_lp == 1.0, name=f\"assign_lp_{j}\")\n", + "\n", + "print(f\"Assignment constraints added (fractional allowed): {len(J)}\")\n", + "\n", + "# 2) DC capacity and logical linking: sum_j d_j x_ij <= capacity[i] * y_i\n", + "for i in I:\n", + " expr_lp = LinearExpression([], [], 0.0)\n", + " for j in J:\n", + " expr_lp += float(demand[j]) * x_lp[i, j]\n", + " expr_lp -= float(capacity[i]) * y_lp[i]\n", + " prob_lp.addConstraint(expr_lp <= 0.0, name=f\"capacity_lp_{i}\")\n", + "\n", + "print(f\"Capacity constraints added: {len(I)}\")\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# -----------------------------\n", + "# Solve the semi-relaxed problem\n", + "# -----------------------------\n", + "\n", + "settings_lp = SolverSettings()\n", + "settings_lp.set_parameter(\"time_limit\", \"300.0\")\n", + "\n", + "prob_lp.solve(settings=settings_lp)\n", + "\n", + "print(f\"Semi-Relaxed Solve status: {prob_lp.Status}\")\n", + "print(f\"Semi-Relaxed Objective value: ${prob_lp.ObjValue:,.2f}\")\n", + "print(f\"Solve time: {prob_lp.SolveTime:.3f}s\")\n", + "\n", + "# Compare with MILP solution (only if both solved successfully)\n", + "if prob.Status in (1, 2) and prob_lp.Status in (1, 2):\n", + " gap = prob.ObjValue - prob_lp.ObjValue\n", + " print(f\"\\nComparison:\")\n", + " print(f\" MILP (full integer) objective: ${prob.ObjValue:,.2f}\")\n", + " print(f\" Semi-relaxed objective: ${prob_lp.ObjValue:,.2f}\")\n", + " if prob.ObjValue > 0:\n", + " print(f\" Integrality gap: ${gap:,.2f} ({gap/prob.ObjValue*100:.2f}%)\")\n", + " else:\n", + " print(f\" Integrality gap: ${gap:,.2f}\")\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# -----------------------------\n", + "# Analyze fractional assignment solution\n", + "# -----------------------------\n", + "\n", + "if prob_lp.Status not in (1, 2):\n", + " status_msg = {0: \"Unknown\", 3: \"Infeasible\", 4: \"Unbounded\", 5: \"Time Limit\"}.get(prob_lp.Status, \"Unknown\")\n", + " raise RuntimeError(f\"Semi-relaxed solver did not find a solution: Status = {prob_lp.Status} ({status_msg})\")\n", + "\n", + "# Open DCs (still integer)\n", + "open_dcs_lp = [i for i in I if y_lp[i].getValue() > 0.5]\n", + "print(\"Open DCs (relaxed LP):\", open_dcs_lp)\n", + "\n", + "# Fractional assignment matrix\n", + "assign_matrix_lp = np.zeros((num_dcs, num_customers))\n", + "for (i, j) in A:\n", + " val = x_lp[i, j].getValue()\n", + " assign_matrix_lp[i, j] = val\n", + "\n", + "# For each customer, show fractional assignments\n", + "print(\"\\nFractional assignments (customer -> DC: fraction):\")\n", + "for j in J:\n", + " assignments = [(i, assign_matrix_lp[i, j]) for i in I if assign_matrix_lp[i, j] > 1e-6]\n", + " if len(assignments) > 1: # Only show customers with split assignments\n", + " assignment_str = \", \".join([f\"DC{i}: {frac:.3f}\" for i, frac in assignments])\n", + " print(f\" Customer {j}: {assignment_str}\")\n", + "\n", + "# Load per DC\n", + "dc_load_lp = assign_matrix_lp @ demand\n", + "print(f\"\\nDC loads (relaxed LP): {dc_load_lp.round(0)}\")\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Extension: Multi-Period Capacity Expansion\n", + "\n", + "This section extends the CFLP to a multi-period setting where we decide **when** to open DCs (now vs defer) and handle time-varying demand. Key features:\n", + "\n", + "- **Multiple time periods**: Plan over T periods (e.g., quarters or years)\n", + "- **Deferred opening decisions**: DCs can be opened in any period, not just period 0\n", + "- **Time-discounted costs**: Future costs are discounted using a discount factor\n", + "- **Time-varying demand**: Customer demand may change over periods\n", + "- **Capacity constraints**: DCs can only serve customers after they're opened\n", + "- **Objective**: Minimize total discounted cost over all periods\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# -----------------------------\n", + "# Multi-period problem setup\n", + "# -----------------------------\n", + "\n", + "# Number of time periods\n", + "num_periods = 3 # e.g., periods 0, 1, 2 (now, next year, year after)\n", + "T = list(range(num_periods))\n", + "\n", + "# Discount factor (e.g., 0.95 means $1 next period is worth $0.95 today)\n", + "discount_factor = 0.95\n", + "\n", + "# Time-varying demand: demand may grow over time\n", + "demand_growth_rate = 1.1 # 10% growth per period\n", + "demand_period = {}\n", + "for t in T:\n", + " demand_period[t] = demand * (demand_growth_rate ** t)\n", + "\n", + "# Time-varying fixed costs: opening costs may change over time\n", + "cost_inflation = 1.05 # 5% cost increase per period\n", + "fixed_cost_period = {}\n", + "for t in T:\n", + " fixed_cost_period[t] = fixed_cost * (cost_inflation ** t)\n", + "\n", + "# Transportation costs may also vary over time\n", + "unit_cost_period = {}\n", + "for t in T:\n", + " unit_cost_period[t] = unit_cost * (1.0 + 0.02 * t) # 2% increase per period\n", + "\n", + "print(f\"Multi-period setup:\")\n", + "print(f\" Number of periods: {num_periods}\")\n", + "print(f\" Discount factor: {discount_factor}\")\n", + "print(f\" Demand growth rate: {(demand_growth_rate - 1) * 100:.1f}% per period\")\n", + "print(f\" Cost inflation: {(cost_inflation - 1) * 100:.1f}% per period\")\n", + "print(f\"\\nTotal demand by period:\")\n", + "for t in T:\n", + " total_demand_t = demand_period[t].sum()\n", + " print(f\" Period {t}: {total_demand_t:,.1f} pallets\")\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# -----------------------------\n", + "# Multi-period problem and variable creation\n", + "# -----------------------------\n", + "\n", + "prob_mp = Problem(\"CFLP_MultiPeriod\")\n", + "\n", + "# Decision variables\n", + "y_mp = {} # y_mp[i][t]: 1 if DC i is open at the START of period t (state variable)\n", + "z_mp = {} # z_mp[i][t]: 1 if DC i opens specifically in period t (decision variable)\n", + "x_mp = {} # x_mp[i][j][t]: 1 if customer j is assigned to DC i in period t\n", + "\n", + "# y_mp[i][t]: State variable tracking if DC i is open at start of period t\n", + "for i in I:\n", + " y_mp[i] = {}\n", + " for t in T:\n", + " y_mp[i][t] = prob_mp.addVariable(\n", + " name=f\"y_{i}_{t}\",\n", + " lb=0.0,\n", + " ub=1.0,\n", + " vtype=VType.INTEGER\n", + " )\n", + "\n", + "# z_mp[i][t]: 1 if DC i opens specifically in period t\n", + "for i in I:\n", + " z_mp[i] = {}\n", + " for t in T:\n", + " z_mp[i][t] = prob_mp.addVariable(\n", + " name=f\"z_{i}_{t}\",\n", + " lb=0.0,\n", + " ub=1.0,\n", + " vtype=VType.INTEGER\n", + " )\n", + "\n", + "# x_mp[i][j][t]: assignment in period t\n", + "for i in I:\n", + " x_mp[i] = {}\n", + " for j in J:\n", + " x_mp[i][j] = {}\n", + " for t in T:\n", + " x_mp[i][j][t] = prob_mp.addVariable(\n", + " name=f\"x_{i}_{j}_{t}\",\n", + " lb=0.0,\n", + " ub=1.0,\n", + " vtype=VType.INTEGER\n", + " )\n", + "\n", + "print(f\"Multi-period - Number of variables: {prob_mp.NumVariables}\")\n", + "print(f\" DC opening state variables (y): {num_dcs * num_periods}\")\n", + "print(f\" DC opening decision variables (z): {num_dcs * num_periods}\")\n", + "print(f\" Assignment variables (x): {num_dcs * num_customers * num_periods}\")\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# -----------------------------\n", + "# Multi-period constraints\n", + "# -----------------------------\n", + "\n", + "# 1) DC opening logic constraints\n", + "for i in I:\n", + " # Period 0: y_mp[i][0] = z_mp[i][0]\n", + " expr = LinearExpression([], [], 0.0)\n", + " expr += y_mp[i][0]\n", + " expr -= z_mp[i][0]\n", + " prob_mp.addConstraint(expr == 0.0, name=f\"open_logic_{i}_0\")\n", + " \n", + " # Period t > 0: State transition constraints\n", + " for t in T[1:]:\n", + " # Monotonicity: once open, stays open\n", + " expr = LinearExpression([], [], 0.0)\n", + " expr += y_mp[i][t]\n", + " expr -= y_mp[i][t-1]\n", + " prob_mp.addConstraint(expr >= 0.0, name=f\"monotone_{i}_{t}\")\n", + " \n", + " # If opens in period t, must be open\n", + " expr = LinearExpression([], [], 0.0)\n", + " expr += y_mp[i][t]\n", + " expr -= z_mp[i][t]\n", + " prob_mp.addConstraint(expr >= 0.0, name=f\"open_if_opens_{i}_{t}\")\n", + " \n", + " # Can only be open if was open before OR opens now\n", + " expr = LinearExpression([], [], 0.0)\n", + " expr += y_mp[i][t]\n", + " expr -= y_mp[i][t-1]\n", + " expr -= z_mp[i][t]\n", + " prob_mp.addConstraint(expr <= 0.0, name=f\"open_logic_{i}_{t}\")\n", + "\n", + "# 2) Each DC can open at most once\n", + "for i in I:\n", + " expr = LinearExpression([], [], 0.0)\n", + " for t in T:\n", + " expr += z_mp[i][t]\n", + " prob_mp.addConstraint(expr <= 1.0, name=f\"open_once_{i}\")\n", + "\n", + "print(\"DC opening logic constraints added\")\n", + "\n", + "# 3) Assignment constraints: Each customer fully assigned in each period\n", + "for j in J:\n", + " for t in T:\n", + " expr = LinearExpression([], [], 0.0)\n", + " for i in I:\n", + " expr += x_mp[i][j][t]\n", + " prob_mp.addConstraint(expr == 1.0, name=f\"assign_{j}_{t}\")\n", + "\n", + "print(f\"Assignment constraints added: {num_customers * num_periods}\")\n", + "\n", + "# 4) Capacity constraints\n", + "for i in I:\n", + " for t in T:\n", + " expr = LinearExpression([], [], 0.0)\n", + " demand_t = demand_period[t]\n", + " for j in J:\n", + " d_jt = demand_t[j]\n", + " expr += float(d_jt) * x_mp[i][j][t]\n", + " expr -= float(capacity[i]) * y_mp[i][t]\n", + " prob_mp.addConstraint(expr <= 0.0, name=f\"capacity_{i}_{t}\")\n", + "\n", + "print(f\"Capacity constraints added: {num_dcs * num_periods}\")\n", + "print(f\"\\nTotal constraints: {prob_mp.NumConstraints}\")\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# -----------------------------\n", + "# Multi-period objective function\n", + "# -----------------------------\n", + "\n", + "obj_mp = LinearExpression([], [], 0.0)\n", + "\n", + "# Fixed opening costs (discounted)\n", + "for i in I:\n", + " for t in T:\n", + " discount = discount_factor ** t\n", + " cost_t = fixed_cost_period[t][i]\n", + " obj_mp += float(cost_t * discount) * z_mp[i][t]\n", + "\n", + "# Transportation costs (discounted)\n", + "for i in I:\n", + " for j in J:\n", + " for t in T:\n", + " discount = discount_factor ** t\n", + " demand_t = demand_period[t]\n", + " d_jt = demand_t[j]\n", + " cost_t = unit_cost_period[t][i, j]\n", + " flow_cost = float(cost_t * d_jt * discount)\n", + " obj_mp += flow_cost * x_mp[i][j][t]\n", + "\n", + "prob_mp.setObjective(obj_mp, sense=sense.MINIMIZE)\n", + "\n", + "print(\"Multi-period objective set (minimize total discounted cost).\")\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# -----------------------------\n", + "# Solve multi-period problem\n", + "# -----------------------------\n", + "\n", + "settings_mp = SolverSettings()\n", + "settings_mp.set_parameter(\"time_limit\", \"600.0\") # Longer time limit\n", + "settings_mp.set_parameter(\"mip_relative_gap\", \"1e-4\")\n", + "\n", + "print(\"Starting multi-period optimization...\")\n", + "print(f\"Problem size: {prob_mp.NumVariables} variables, {prob_mp.NumConstraints} constraints\")\n", + "print()\n", + "\n", + "start_time = time.time()\n", + "prob_mp.solve(settings=settings_mp)\n", + "solve_time_mp = time.time() - start_time\n", + "\n", + "print(f\"\\nOptimization completed in {solve_time_mp:.3f} seconds\")\n", + "print(f\"Multi-period Solve status: {prob_mp.Status}\")\n", + "print(f\"Multi-period Objective value (discounted): ${prob_mp.ObjValue:,.2f}\")\n", + "print(f\"Solve time (internal): {prob_mp.SolveTime:.3f}s\")\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# -----------------------------\n", + "# Analyze multi-period solution\n", + "# -----------------------------\n", + "\n", + "if prob_mp.Status not in (1, 2):\n", + " status_msg = {0: \"Unknown\", 3: \"Infeasible\", 4: \"Unbounded\", 5: \"Time Limit\"}.get(prob_mp.Status, \"Unknown\")\n", + " raise RuntimeError(f\"Multi-period solver did not find a solution: Status = {prob_mp.Status} ({status_msg})\")\n", + "\n", + "# Extract opening decisions by period\n", + "opening_schedule = {}\n", + "dc_open_period = {}\n", + "\n", + "for t in T:\n", + " opening_schedule[t] = []\n", + " for i in I:\n", + " if z_mp[i][t].getValue() > 0.5:\n", + " opening_schedule[t].append(i)\n", + " dc_open_period[i] = t\n", + "\n", + "never_open = [i for i in I if i not in dc_open_period]\n", + "\n", + "print(\"DC Opening Schedule:\")\n", + "for t in T:\n", + " if opening_schedule[t]:\n", + " print(f\" Period {t}: Open DCs {opening_schedule[t]}\")\n", + " else:\n", + " print(f\" Period {t}: No DCs opened\")\n", + "if never_open:\n", + " print(f\" Never opened: DCs {never_open}\")\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Calculate costs by period\n", + "costs_by_period = {}\n", + "for t in T:\n", + " opening_cost = 0.0\n", + " transport_cost = 0.0\n", + " \n", + " # Opening costs\n", + " for i in opening_schedule[t]:\n", + " cost_t = fixed_cost_period[t][i]\n", + " opening_cost += cost_t\n", + " \n", + " # Transportation costs\n", + " for i in I:\n", + " for j in J:\n", + " if x_mp[i][j][t].getValue() > 0.5:\n", + " demand_t = demand_period[t]\n", + " d_jt = demand_t[j]\n", + " cost_t = unit_cost_period[t][i, j]\n", + " transport_cost += cost_t * d_jt\n", + " \n", + " discount = discount_factor ** t\n", + " costs_by_period[t] = {\n", + " 'opening_cost': opening_cost,\n", + " 'transport_cost': transport_cost,\n", + " 'total_cost': opening_cost + transport_cost,\n", + " 'discounted_cost': (opening_cost + transport_cost) * discount\n", + " }\n", + "\n", + "# Create summary DataFrame\n", + "summary_data = []\n", + "for t in T:\n", + " summary_data.append({\n", + " 'Period': t,\n", + " 'DCs Opened': len(opening_schedule[t]),\n", + " 'Opening Cost ($)': costs_by_period[t]['opening_cost'],\n", + " 'Transport Cost ($)': costs_by_period[t]['transport_cost'],\n", + " 'Total Cost ($)': costs_by_period[t]['total_cost'],\n", + " 'Discounted Cost ($)': costs_by_period[t]['discounted_cost'],\n", + " 'Total Demand': demand_period[t].sum()\n", + " })\n", + "\n", + "summary_df = pd.DataFrame(summary_data)\n", + "print(\"\\nMulti-Period Cost Summary:\")\n", + "display(summary_df.round(2))\n", + "\n", + "print(f\"\\nTotal discounted cost: ${prob_mp.ObjValue:,.2f}\")\n", + "print(f\"Total undiscounted cost: ${summary_df['Total Cost ($)'].sum():,.2f}\")\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# -----------------------------\n", + "# Visualize multi-period solution\n", + "# -----------------------------\n", + "\n", + "fig, axes = plt.subplots(1, 2, figsize=(14, 5))\n", + "\n", + "# Plot 1: Opening schedule timeline\n", + "ax1 = axes[0]\n", + "for i in I:\n", + " if i in dc_open_period:\n", + " t_open = dc_open_period[i]\n", + " ax1.barh(i, num_periods - t_open, left=t_open, alpha=0.6, color='tab:green')\n", + " ax1.text(t_open + 0.1, i, f'DC {i} opens', va='center', fontsize=9)\n", + " else:\n", + " ax1.barh(i, 0.1, left=0, alpha=0.3, color='gray')\n", + " ax1.text(0.15, i, f'DC {i} (never opened)', va='center', fontsize=9, color='gray')\n", + "\n", + "ax1.set_xlabel('Period')\n", + "ax1.set_ylabel('DC')\n", + "ax1.set_title('DC Opening Schedule', fontsize=14, fontweight='bold')\n", + "ax1.set_yticks(I)\n", + "ax1.set_yticklabels([f'DC {i}' for i in I])\n", + "ax1.set_xticks(T)\n", + "ax1.set_xlim(-0.2, num_periods)\n", + "ax1.grid(True, alpha=0.3)\n", + "\n", + "# Plot 2: Costs by period\n", + "ax2 = axes[1]\n", + "periods = summary_df['Period']\n", + "width = 0.35\n", + "x = np.arange(len(periods))\n", + "\n", + "ax2.bar(x - width/2, summary_df['Opening Cost ($)'], width, label='Opening Cost', alpha=0.8, color='tab:red')\n", + "ax2.bar(x + width/2, summary_df['Transport Cost ($)'], width, label='Transport Cost', alpha=0.8, color='tab:blue')\n", + "\n", + "ax2.set_xlabel('Period')\n", + "ax2.set_ylabel('Cost ($)')\n", + "ax2.set_title('Costs by Period (Undiscounted)', fontsize=14, fontweight='bold')\n", + "ax2.set_xticks(x)\n", + "ax2.set_xticklabels([f'Period {t}' for t in periods])\n", + "ax2.legend()\n", + "ax2.grid(True, alpha=0.3, axis='y')\n", + "\n", + "plt.tight_layout()\n", + "plt.show()\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Summary and Business Insights\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "print(\"\\n\" + \"=\"*70)\n", + "print(\" \" * 20 + \"BUSINESS INSIGHTS\")\n", + "print(\"=\"*70 + \"\\n\")\n", + "\n", + "print(\"📊 Single-Period MILP Results:\")\n", + "print(f\" • Optimal Cost: ${prob.ObjValue:,.2f}\")\n", + "print(f\" • DCs Opened: {len(open_dcs)} of {num_dcs}\")\n", + "print(f\" • Solve Time: {solve_time:.3f} seconds\")\n", + "print()\n", + "\n", + "print(\"📈 Multi-Period Results:\")\n", + "print(f\" • Total Discounted Cost: ${prob_mp.ObjValue:,.2f}\")\n", + "print(f\" • Planning Horizon: {num_periods} periods\")\n", + "print(f\" • DCs Eventually Opened: {len(dc_open_period)} of {num_dcs}\")\n", + "print(f\" • Solve Time: {solve_time_mp:.3f} seconds\")\n", + "print()\n", + "\n", + "print(\"🔄 Semi-Relaxation Insights:\")\n", + "if prob_lp.Status in (1, 2) and prob.Status in (1, 2):\n", + " gap = prob.ObjValue - prob_lp.ObjValue\n", + " print(f\" • Semi-Relaxed Objective: ${prob_lp.ObjValue:,.2f}\")\n", + " if prob.ObjValue > 0:\n", + " print(f\" • Integrality Gap: ${gap:,.2f} ({gap/prob.ObjValue*100:.2f}%)\")\n", + " else:\n", + " print(f\" • Integrality Gap: ${gap:,.2f}\")\n", + "else:\n", + " print(\" • Solution not available\")\n", + "print()\n", + "\n", + "print(\"⚡ Performance:\")\n", + "print(f\" • GPU Acceleration: Enabled ✓\")\n", + "print(f\" • Single-period problem: {prob.NumVariables} variables, {prob.NumConstraints} constraints\")\n", + "print(f\" • Multi-period problem: {prob_mp.NumVariables} variables, {prob_mp.NumConstraints} constraints\")\n", + "print(\"\\n\" + \"=\"*70)\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Conclusion\n", + "\n", + "This notebook demonstrated how to use the cuOpt linear programming API to solve a Capacitated Facility Location Problem (CFLP).\n", + "\n", + "### What We Achieved:\n", + "1. **Optimal DC selection** to minimize total logistics cost\n", + "2. **Customer-to-DC assignment** respecting capacity constraints\n", + "3. **LP relaxation** for lower bound analysis\n", + "4. **Multi-period extension** for capacity expansion planning\n", + "\n", + "### Key Advantages of cuOpt:\n", + "- **GPU Acceleration**: Significantly faster than CPU-based solvers\n", + "- **Scalability**: Handles problems with thousands of variables and constraints\n", + "- **Integration**: Easy to integrate into existing Python workflows\n", + "- **Real-time**: Fast enough for operational decision-making\n", + "\n", + "### Possible Extensions:\n", + "\n", + "**High Priority - Service Levels:**\n", + "- Maximum distance constraints (customers must be within X km of assigned DC)\n", + "- Multiple service tiers with different SLAs\n", + "- Backup DC assignments for redundancy\n", + "\n", + "**Additional Enhancements:**\n", + "- Multiple product types with different handling requirements\n", + "- DC capacity expansion options (modular capacity)\n", + "- Stochastic demand scenarios\n", + "- Multi-echelon supply chain (plants → DCs → customers)\n", + "- Inventory considerations and safety stock\n", + "- Seasonal demand patterns\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "SPDX-FileCopyrightText: Copyright (c) 2025 NVIDIA CORPORATION & AFFILIATES. All rights reserved.\n", + "\n", + "SPDX-License-Identifier: Apache-2.0\n", + "\n", + "Licensed under the Apache License, Version 2.0 (the \"License\");\n", + "you may not use this file except in compliance with the License.\n", + "You may obtain a copy of the License at\n", + "\n", + "http://www.apache.org/licenses/LICENSE-2.0\n", + "\n", + "Unless required by applicable law or agreed to in writing, software\n", + "distributed under the License is distributed on an \"AS IS\" BASIS,\n", + "WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.\n", + "See the License for the specific language governing permissions and\n", + "limitations under the License.\n" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "cuopt_fresh", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.12.12" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +}