inventory-routing-problem

库存路径问题（IRP）

你是库存路径问题（IRP）和库存-配送集成优化领域的专家。你的目标是协同优化库存管理和车辆路径决策，最大限度地降低系统总成本，包括库存持有成本、路径成本以及潜在的缺货成本。

初步评估

在解决库存路径问题之前，请先了解：

1. 系统结构

   • 由供应商管理库存（VMI）还是由零售商管理？
   • 有多少客户/零售商？
   • 是单个仓库还是多个仓库？
   • 规划周期是多久（天、周）？
   • 配送频率如何？

2. 库存特征

   • 各地点的存储容量是多少？
   • 当前库存水平如何？
   • 消耗/需求率是确定的还是随机的？
   • 最低库存水平（安全库存）是多少？
   • 最高库存水平（储罐容量、货架空间）是多少？
   • 产品是否有保质期或易腐？

3. 路径约束

   • 车辆容量（重量、体积）是多少？
   • 可用车辆数量是多少？
   • 路线持续时间或距离的上限是多少？
   • 配送是否有时间窗？
   • 是否有司机班次约束？
   • 是否存在通行限制？

4. 成本结构

   • 仓库和客户处的库存持有成本是多少？
   • 运输成本是多少（按英里、按车辆、按路线）？
   • 每辆启用车辆的固定成本是多少？
   • 缺货惩罚成本是多少？
   • 每次客户访问的设置/配送费用是多少？

5. 服务要求

   • 是否必须避免缺货？
   • 每位客户的最低服务频率是多少？
   • 是否有优先级客户？
   • 是否有合同规定的配送要求？

---

IRP 基础知识

问题定义

库存路径问题（IRP）整合了两个经典问题：

1. 库存管理： 何时以及补充多少库存给每位客户
2. 车辆路径： 如何高效规划车辆路线以服务客户

关键权衡：

• 更频繁的小批量配送 → 路径成本更高，库存更低
• 较少频次的大批量配送 → 路径成本更低，库存更高

问题变体

1. 单周期 IRP

• 一次性的路径和配送决策
• 已知当前库存水平
• 在满足库存约束的条件下，最小化路径成本

2. 多周期 IRP

• 规划整个时间范围内（T 个周期）的配送
• 考虑库存动态变化
• 最贴近实际，也最复杂

3. 确定性与随机 IRP

• 确定性：消耗率已知
• 随机：需求不确定，需要安全库存

4. 海运 IRP（MIRP）

• 使用船舶而非卡车
• 容量更大，运输时间更长
• 常用于汽油/化学品配送

---

Python 实现：IRP 模型

使用 MIP 的单周期 IRP

import numpy as np
import pandas as pd
from pulp import *
from typing import List, Dict, Tuple
import matplotlib.pyplot as plt
from scipy.spatial.distance import cdist

class SinglePeriodIRP:
    """
    Single-Period Inventory Routing Problem

    Given:
    - Current inventory at each customer
    - Consumption rates
    - Vehicle capacity
    - Distance matrix

    Decide:
    - Which customers to visit
    - How much to deliver to each
    - Vehicle routes
    """

    def __init__(self, num_customers: int, customer_locations: np.ndarray,
                 depot_location: np.ndarray, current_inventory: np.ndarray,
                 consumption_rates: np.ndarray, max_inventory: np.ndarray,
                 vehicle_capacity: float, num_vehicles: int,
                 holding_cost: float = 1.0, routing_cost_per_km: float = 1.0):
        """
        Parameters:
        -----------
        num_customers : int
            Number of customer locations
        customer_locations : ndarray
            (n x 2) array of customer coordinates
        depot_location : ndarray
            (2,) depot coordinates
        current_inventory : ndarray
            Current inventory level at each customer
        consumption_rates : ndarray
            Daily consumption at each customer
        max_inventory : ndarray
            Maximum storage capacity at each customer
        vehicle_capacity : float
            Vehicle capacity (units)
        num_vehicles : int
            Number of vehicles available
        holding_cost : float
            Inventory holding cost per unit per day
        routing_cost_per_km : float
            Cost per kilometer traveled
        """
        self.n = num_customers
        self.customer_locations = customer_locations
        self.depot = depot_location
        self.I = current_inventory
        self.d = consumption_rates
        self.C = max_inventory
        self.Q = vehicle_capacity
        self.K = num_vehicles
        self.h = holding_cost
        self.c_routing = routing_cost_per_km

        # Calculate distance matrix
        all_locations = np.vstack([depot_location, customer_locations])
        self.dist_matrix = cdist(all_locations, all_locations, metric='euclidean')

    def solve_mip(self, time_until_next_delivery: int = 1) -> Dict:
        """
        Solve single-period IRP using Mixed-Integer Programming

        Decision variables:
        - x[i,j,k]: binary, 1 if vehicle k travels from i to j
        - y[i,k]: binary, 1 if customer i is visited by vehicle k
        - q[i]: quantity delivered to customer i
        """

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

        # Nodes: 0 = depot, 1..n = customers
        nodes = range(self.n + 1)
        customers = range(1, self.n + 1)
        vehicles = range(self.K)

        # Decision variables
        # Routing variables
        x = {}
        for i in nodes:
            for j in nodes:
                for k in vehicles:
                    if i != j:
                        x[i, j, k] = LpVariable(f"x_{i}_{j}_{k}", cat='Binary')

        # Visit variables
        y = {}
        for i in customers:
            for k in vehicles:
                y[i, k] = LpVariable(f"y_{i}_{k}", cat='Binary')

        # Delivery quantities
        q = {i: LpVariable(f"q_{i}", lowBound=0) for i in customers}

        # Objective: Minimize routing cost + inventory holding cost
        routing_cost = lpSum([
            self.dist_matrix[i, j] * self.c_routing * x[i, j, k]
            for i in nodes for j in nodes for k in vehicles if i != j
        ])

        # Inventory after delivery
        inventory_after = {i: self.I[i - 1] + q[i] for i in customers}
        holding_cost = self.h * lpSum([inventory_after[i] for i in customers])

        prob += routing_cost + holding_cost

        # Constraints

        # 1. Each customer visited at most once
        for i in customers:
            prob += lpSum([y[i, k] for k in vehicles]) <= 1

        # 2. Visit variable linking
        for i in customers:
            for k in vehicles:
                prob += lpSum([x[j, i, k] for j in nodes if j != i]) == y[i, k]
                prob += lpSum([x[i, j, k] for j in nodes if j != i]) == y[i, k]

        # 3. Vehicle starts and ends at depot
        for k in vehicles:
            prob += lpSum([x[0, j, k] for j in customers]) <= 1
            prob += lpSum([x[i, 0, k] for i in customers]) <= 1
            prob += (lpSum([x[0, j, k] for j in customers]) ==
                    lpSum([x[i, 0, k] for i in customers]))

        # 4. Flow conservation
        for k in vehicles:
            for j in customers:
                prob += (lpSum([x[i, j, k] for i in nodes if i != j]) ==
                        lpSum([x[j, i, k] for i in nodes if i != j]))

        # 5. Vehicle capacity
        for k in vehicles:
            prob += lpSum([q[i] * y[i, k] for i in customers]) <= self.Q

        # 6. Delivery quantity constraints
        for i in customers:
            # Don't deliver more than capacity minus current inventory
            prob += q[i] <= (self.C[i - 1] - self.I[i - 1]) * lpSum([y[i, k]
                                                                      for k in vehicles])

            # If visited, deliver enough to avoid stockout until next delivery
            min_delivery = max(0, time_until_next_delivery * self.d[i - 1] - self.I[i - 1])
            prob += q[i] >= min_delivery * lpSum([y[i, k] for k in vehicles])

        # 7. Subtour elimination (MTZ formulation)
        u = {i: LpVariable(f"u_{i}", lowBound=0, upBound=self.n) for i in customers}

        for i in customers:
            for j in customers:
                for k in vehicles:
                    if i != j:
                        prob += u[i] - u[j] + self.n * x[i, j, k] <= self.n - 1

        # Solve
        prob.solve(PULP_CBC_CMD(msg=0))

        # Extract solution
        routes = self._extract_routes(x, vehicles, nodes)
        deliveries = {i: q[i].varValue if q[i].varValue else 0 for i in customers}

        total_distance = sum(
            self.dist_matrix[i, j] * x[i, j, k].varValue
            for i in nodes for j in nodes for k in vehicles
            if i != j and x[i, j, k].varValue > 0.5
        )

        total_delivery = sum(deliveries.values())

        return {
            'status': LpStatus[prob.status],
            'routes': routes,
            'deliveries': deliveries,
            'total_cost': value(prob.objective),
            'routing_cost': self.c_routing * total_distance,
            'holding_cost': self.h * sum(self.I[i - 1] + deliveries[i]
                                         for i in customers),
            'total_distance': total_distance,
            'total_delivery': total_delivery,
            'vehicles_used': len([r for r in routes if len(r) > 2])
        }

    def _extract_routes(self, x, vehicles, nodes):
        """Extract route sequences from solution"""

        routes = []

        for k in vehicles:
            route = [0]  # Start at depot
            current = 0

            while True:
                next_node = None
                for j in nodes:
                    if j != current and (current, j, k) in x:
                        if x[current, j, k].varValue > 0.5:
                            next_node = j
                            break

                if next_node is None or next_node == 0:
                    if len(route) > 1:
                        route.append(0)  # Return to depot
                        routes.append(route)
                    break

                route.append(next_node)
                current = next_node

        return routes

    def plot_solution(self, solution: Dict):
        """Visualize IRP solution"""

        fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(16, 6))

        # Plot 1: Routes
        colors = plt.cm.tab10(np.linspace(0, 1, len(solution['routes'])))

        # Plot depot
        ax1.plot(self.depot[0], self.depot[1], 'rs', markersize=15,
                label='Depot', zorder=5)

        # Plot customers
        for i in range(self.n):
            ax1.plot(self.customer_locations[i, 0],
                    self.customer_locations[i, 1],
                    'bo', markersize=10, zorder=3)
            ax1.text(self.customer_locations[i, 0],
                    self.customer_locations[i, 1],
                    f'  {i+1}', fontsize=9)

        # Plot routes
        for route_idx, route in enumerate(solution['routes']):
            if len(route) > 2:
                route_coords = np.vstack([
                    self.depot if node == 0
                    else self.customer_locations[node - 1]
                    for node in route
                ])

                ax1.plot(route_coords[:, 0], route_coords[:, 1],
                        'o-', color=colors[route_idx], linewidth=2,
                        markersize=8, label=f'Route {route_idx + 1}',
                        alpha=0.7)

        ax1.set_xlabel('X Coordinate')
        ax1.set_ylabel('Y Coordinate')
        ax1.set_title('Vehicle Routes', fontweight='bold')
        ax1.legend()
        ax1.grid(True, alpha=0.3)

        # Plot 2: Inventory levels
        customer_ids = np.arange(1, self.n + 1)
        current_inv = self.I
        deliveries = [solution['deliveries'][i] for i in customer_ids]
        final_inv = current_inv + deliveries
        capacity = self.C

        x_pos = np.arange(self.n)
        width = 0.35

        ax2.bar(x_pos - width/2, current_inv, width, label='Current Inventory',
               alpha=0.7, color='orange')
        ax2.bar(x_pos + width/2, final_inv, width, label='After Delivery',
               alpha=0.7, color='green')
        ax2.plot(x_pos, capacity, 'r--', linewidth=2, label='Max Capacity')

        # Mark delivered customers
        delivered_customers = [i for i in customer_ids if deliveries[i - 1] > 0]
        if delivered_customers:
            ax2.scatter([c - 1 for c in delivered_customers],
                       [final_inv[c - 1] for c in delivered_customers],
                       s=200, marker='*', color='red', zorder=5,
                       label='Delivered')

        ax2.set_xlabel('Customer ID')
        ax2.set_ylabel('Inventory Level (units)')
        ax2.set_title('Inventory Levels Before and After', fontweight='bold')
        ax2.set_xticks(x_pos)
        ax2.set_xticklabels(customer_ids)
        ax2.legend()
        ax2.grid(True, alpha=0.3, axis='y')

        plt.tight_layout()
        return plt


# Example Usage
def example_single_period_irp():
    """Example: Single-period IRP with 8 customers"""

    print("\n" + "=" * 70)
    print("INVENTORY ROUTING PROBLEM (IRP): SINGLE-PERIOD")
    print("=" * 70)

    np.random.seed(42)

    # Problem setup
    num_customers = 8
    depot = np.array([50, 50])

    # Customer locations (random)
    customer_locations = np.random.rand(num_customers, 2) * 100

    # Current inventory (random, 20-80% of capacity)
    max_inventory = np.random.randint(80, 150, num_customers)
    current_inventory = max_inventory * np.random.uniform(0.2, 0.8, num_customers)

    # Consumption rates (units per day)
    consumption_rates = np.random.uniform(5, 20, num_customers)

    # Days until will stock out if not replenished
    days_until_stockout = current_inventory / consumption_rates

    print("\nProblem Data:")
    print(f"  Number of Customers: {num_customers}")
    print(f"  Vehicle Capacity: 200 units")
    print(f"  Number of Vehicles: 3")
    print(f"  Planning Period: 1 day")

    print("\n  Customer Inventory Status:")
    print(f"\n  {'Customer':<12} {'Current':<12} {'Capacity':<12} {'Usage/Day':<12} "
          f"{'Days to Stockout'}")
    print("  " + "-" * 65)

    for i in range(num_customers):
        print(f"  {i+1:<12} {current_inventory[i]:<12.0f} "
              f"{max_inventory[i]:<12.0f} {consumption_rates[i]:<12.1f} "
              f"{days_until_stockout[i]:<.1f}")

    # Create and solve IRP
    irp = SinglePeriodIRP(
        num_customers=num_customers,
        customer_locations=customer_locations,
        depot_location=depot,
        current_inventory=current_inventory,
        consumption_rates=consumption_rates,
        max_inventory=max_inventory,
        vehicle_capacity=200,
        num_vehicles=3,
        holding_cost=1.0,
        routing_cost_per_km=2.0
    )

    print("\nSolving IRP with MIP...")
    solution = irp.solve_mip(time_until_next_delivery=3)  # Plan for 3 days

    print(f"\n{'=' * 70}")
    print("OPTIMAL SOLUTION")
    print("=" * 70)

    print(f"\n{'Status:':<30} {solution['status']}")
    print(f"{'Total Cost:':<30} ${solution['total_cost']:,.2f}")
    print(f"{'Routing Cost:':<30} ${solution['routing_cost']:,.2f}")
    print(f"{'Holding Cost:':<30} ${solution['holding_cost']:,.2f}")
    print(f"{'Total Distance:':<30} {solution['total_distance']:.1f} km")
    print(f"{'Vehicles Used:':<30} {solution['vehicles_used']}")
    print(f"{'Total Delivered:':<30} {solution['total_delivery']:.0f} units")

    print("\n  Vehicle Routes and Deliveries:")
    for route_idx, route in enumerate(solution['routes']):
        if len(route) > 2:
            print(f"\n  Route {route_idx + 1}: ", end='')
            print(" → ".join([f"Depot" if node == 0 else f"Cust {node}"
                             for node in route]))

            route_delivery = sum(solution['deliveries'][node]
                               for node in route if node > 0)
            print(f"    Total delivery on route: {route_delivery:.0f} units")

            for node in route:
                if node > 0:
                    delivery = solution['deliveries'][node]
                    if delivery > 0:
                        print(f"      Customer {node}: Deliver {delivery:.0f} units")

    # Plot solution
    irp.plot_solution(solution)
    plt.savefig('/tmp/irp_single_period.png', dpi=300, bbox_inches='tight')
    print(f"\nSolution plot saved to /tmp/irp_single_period.png")

    return irp, solution


if __name__ == "__main__":
    example_single_period_irp()

---

多周期 IRP

滚动时域方法

class MultiPeriodIRP:
    """
    Multi-Period IRP using rolling horizon approach

    Solve single-period IRP repeatedly, updating inventory levels
    """

    def __init__(self, single_period_irp: SinglePeriodIRP,
                 num_periods: int, delivery_frequency: int = 2):
        """
        Parameters:
        -----------
        single_period_irp : SinglePeriodIRP
            Single-period IRP model
        num_periods : int
            Number of periods to plan
        delivery_frequency : int
            Minimum periods between deliveries to same customer
        """
        self.irp = single_period_irp
        self.T = num_periods
        self.freq = delivery_frequency

        # Track inventory over time
        self.inventory_history = np.zeros((num_periods + 1, self.irp.n))
        self.inventory_history[0] = self.irp.I

        # Track deliveries
        self.delivery_history = []

        # Track routes
        self.route_history = []

    def solve_rolling_horizon(self) -> Dict:
        """Solve multi-period IRP using rolling horizon"""

        total_cost = 0
        total_distance = 0
        total_delivered = 0

        for t in range(self.T):
            print(f"  Period {t+1}/{self.T}...", end='')

            # Update current inventory in IRP model
            self.irp.I = self.inventory_history[t]

            # Solve single-period IRP
            solution = self.irp.solve_mip(time_until_next_delivery=self.freq)

            # Record solution
            self.route_history.append(solution['routes'])
            self.delivery_history.append(solution['deliveries'])

            # Update costs
            total_cost += solution['total_cost']
            total_distance += solution['total_distance']
            total_delivered += solution['total_delivery']

            # Update inventory for next period
            for i in range(1, self.irp.n + 1):
                delivered = solution['deliveries'][i]
                consumed = self.irp.d[i - 1]
                self.inventory_history[t + 1, i - 1] = (
                    self.inventory_history[t, i - 1] + delivered - consumed
                )

            print(f" Cost: ${solution['total_cost']:.2f}")

        return {
            'total_cost': total_cost,
            'total_distance': total_distance,
            'total_delivered': total_delivered,
            'avg_cost_per_period': total_cost / self.T,
            'inventory_history': self.inventory_history,
            'delivery_history': self.delivery_history,
            'route_history': self.route_history
        }


# Example: Multi-period
def example_multi_period_irp():
    """Example: 7-day multi-period IRP"""

    print("\n" + "=" * 70)
    print("MULTI-PERIOD INVENTORY ROUTING PROBLEM")
    print("=" * 70)

    np.random.seed(42)

    # Setup (smaller problem for multi-period)
    num_customers = 5
    depot = np.array([50, 50])
    customer_locations = np.random.rand(num_customers, 2) * 100

    max_inventory = np.array([100, 120, 80, 150, 100])
    current_inventory = np.array([80, 90, 60, 100, 70])
    consumption_rates = np.array([10, 15, 8, 20, 12])

    # Create single-period IRP
    irp = SinglePeriodIRP(
        num_customers=num_customers,
        customer_locations=customer_locations,
        depot_location=depot,
        current_inventory=current_inventory,
        consumption_rates=consumption_rates,
        max_inventory=max_inventory,
        vehicle_capacity=150,
        num_vehicles=2,
        holding_cost=0.5,
        routing_cost_per_km=2.0
    )

    # Create multi-period wrapper
    multi_irp = MultiPeriodIRP(irp, num_periods=7, delivery_frequency=2)

    print("\nSolving 7-day multi-period IRP...")
    print("  Using rolling horizon approach")

    solution = multi_irp.solve_rolling_horizon()

    print(f"\n{'=' * 70}")
    print("MULTI-PERIOD RESULTS (7 days)")
    print("=" * 70)

    print(f"\n{'Total Cost (7 days):':<30} ${solution['total_cost']:,.2f}")
    print(f"{'Average Cost per Day:':<30} ${solution['avg_cost_per_period']:,.2f}")
    print(f"{'Total Distance:':<30} {solution['total_distance']:.1f} km")
    print(f"{'Total Delivered:':<30} {solution['total_delivered']:.0f} units")

    # Plot inventory evolution
    fig, axes = plt.subplots(num_customers, 1, figsize=(12, 10))

    for i in range(num_customers):
        axes[i].plot(range(8), solution['inventory_history'][:, i],
                    marker='o', linewidth=2, color='blue')
        axes[i].axhline(y=max_inventory[i], color='red', linestyle='--',
                       label='Capacity')
        axes[i].axhline(y=0, color='black', linestyle='-', linewidth=0.5)

        # Mark delivery days
        for t in range(7):
            if multi_irp.delivery_history[t][i + 1] > 0:
                axes[i].plot(t, solution['inventory_history'][t, i],
                           'g^', markersize=12, label='Delivery' if t == 0 else '')

        axes[i].set_ylabel(f'Cust {i+1}\nInventory')
        axes[i].grid(True, alpha=0.3)
        if i == 0:
            axes[i].legend()

    axes[-1].set_xlabel('Day')
    plt.suptitle('Inventory Evolution Over 7 Days', fontsize=14, fontweight='bold')
    plt.tight_layout()
    plt.savefig('/tmp/irp_multi_period.png', dpi=300, bbox_inches='tight')

    print(f"\nInventory evolution plot saved to /tmp/irp_multi_period.png")

    return multi_irp, solution


if __name__ == "__main__":
    example_multi_period_irp()

---

工具和库

Python 库

• pulp、pyomo：MIP 建模
• ortools：Google OR-Tools 路径规划
• numpy、scipy：数值计算

商业软件

• Blue Yonder TMS：支持 VMI 的运输管理
• Manhattan Associates：集成库存管理的 WMS/TMS
• SAP TM + EWM：集成运输和仓库管理
• Oracle Transportation Management：结合库存的路线优化
• Descartes：考虑库存因素的路径规划

---

常见挑战与解决方案

挑战：问题规模和复杂性

问题： 客户和周期较多时会出现组合爆炸 解决方案：

• 使用滚动时域方法
• 先聚类、后规划路线的启发式方法
• 分解方法
• 限制优化时间，并使用效果良好的启发式方法

挑战：需求不确定性

问题： 消耗率具有随机性 解决方案：

• 在客户处设置安全库存
• 使用需求场景进行鲁棒优化
• 频繁重新规划
• 在仓库进行风险池化

挑战：时间窗口与服务要求

问题： 客户有配送时间窗口和最低服务频率要求 解决方案：

• 在 MIP 中添加时间窗口约束
• 多目标优化（成本与服务）
• 对违约行为设置惩罚成本
• 基于合同的服务水平协议

挑战：异构车队

问题： 车辆类型各不相同（容量、成本） 解决方案：

• 在模型中按类型索引车辆
• 针对不同类型设置路线成本
• 优先使用成本较低的车辆

---

相关技能

• vehicle-routing-problem：纯路线优化
• route-optimization：运输规划
• inventory-optimization：库存管理
• multi-echelon-inventory：网络库存
• network-design：战略配送网络
• fleet-management：车队运营
• demand-forecasting：消耗率预测