{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "7fb27b941602401d91542211134fc71a",
   "metadata": {},
   "source": [
    "# 库存路径问题 (IRP)\n",
    "\n",
    "**类别：** 路径优化\n",
    "\n",
    "使用 OptAgent 的 Python 接口描述变量、约束与目标。\n",
    "\n",
    "问题与原始示例来源：[Hexaly Code Templates](https://www.hexaly.com/templates/inventory-routing-problem-irp)。\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "acae54e37e7d407bbb7b55eff062a284",
   "metadata": {},
   "source": [
    "## 问题描述\n",
    "\n",
    "**库存路径问题 (IRP)** 是一个配送问题，其中需要在给定的时间范围内将产品从一个供应商运送到多个客户。每个客户都有一个最大库存水平。供应商监控每个客户的库存并确定其补货策略，保证客户不会出现缺货（供应商管理库存策略）。从供应商到客户的运输由具有给定容量的车辆完成。\n",
    "\n",
    "### 学习要点\n",
    "\n",
    "- 添加 [列表决策变量](https://optagent.pages.dev/guide/modeling/) 来建模卡车的客户序列\n",
    "- 定义 [符号 Lambda 函数](https://optagent.pages.dev/guide/modeling/) 来计算每个时间槽的行驶距离\n",
    "- 递归地定义表达式来计算库存水平\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9a63283cbaf04dbcab1f6479b197f3a8",
   "metadata": {},
   "source": [
    "## 数据\n",
    "\n",
    "我们提供的库存路径问题 (IRP) 实例来自 [Archetti et al. instances](https://www.leandro-coelho.com/instances/inventory-routing/)。数据文件的格式如下：\n",
    "\n",
    "- 第一行：客户数、规划时间范围内的离散时间槽数以及车辆的运输容量。\n",
    "- 第二行：对于供应商：\n",
    "\n",
    "- 供应商的索引\n",
    "- x 坐标\n",
    "- y 坐标\n",
    "- 供应商库存的起始水平\n",
    "- 每个时间槽可获得的产品数量\n",
    "- 单位库存成本\n",
    "- 接下来的每一行：对每个客户\n",
    "\n",
    "- 客户的索引\n",
    "- x 坐标\n",
    "- y 坐标\n",
    "- 库存的起始水平\n",
    "- 最大库存水平\n",
    "- 最小库存水平\n",
    "- 每个时间槽消耗的产品数量\n",
    "- 单位库存成本\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8dd0d8092fe74a7c96281538738b07e2",
   "metadata": {},
   "source": [
    "## 建模思路\n",
    "\n",
    "库存路径问题 (IRP) 的 OptAgent 模型使用 列表决策变量。对每个离散时间槽，我们定义一个 list variable，表示卡车在该时间槽访问的客户。每条路线的成本取决于卡车行驶的距离。使用 符号 Lambda 函数，我们将沿途从一个客户到下一个客户的距离累加起来。\n",
    "\n",
    "我们递归地计算供应商的库存水平。事实上，时间 t 的库存水平等于时间 t-1 的库存水平，加上时间 t-1 可获得的产品数量，再减去时间 t-1 运送给客户的产品总量。供应商处的缺货约束确保供应商的库存水平始终足以运送给定时间槽内交付给客户的总数量。\n",
    "\n",
    "类似地，我们随后计算客户的库存水平。对于每个客户，时间 t 的库存水平等于时间 t-1 的库存水平，加上时间 t-1 从供应商运送给该客户的产品数量，再减去时间 t-1 消耗的产品数量。客户处的缺货约束保证每个客户的库存水平始终为正。\n",
    "\n",
    "容量约束确保车辆装载的产品总量不超过其容量。最大水平约束保证从供应商运送给每个客户的产品数量不超过该客户的最大库存水平。\n",
    "\n",
    "目标是最小化所有成本之和：供应商的总库存成本、客户的总库存成本以及总运输成本。\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8edb47106e1a46a883d545849b8ab81b",
   "metadata": {},
   "source": [
    "## Python 实现\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "10185d26023b46108eb7d9f57d49d2b3",
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "from pathlib import Path\n",
    "\n",
    "from optagent import OptModel, solve\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "def read_elem(filename):\n",
    "    return Path(filename).read_text(encoding=\"utf-8\").split()\n",
    "\n",
    "\n",
    "def main(instance_file, output_file=None, time_limit=5):\n",
    "    (\n",
    "        nb_customers,\n",
    "        horizon_length,\n",
    "        capacity,\n",
    "        start_level_supplier,\n",
    "        production_rate_supplier,\n",
    "        holding_cost_supplier,\n",
    "        start_level,\n",
    "        max_level,\n",
    "        demand_rate,\n",
    "        holding_cost,\n",
    "        dist_matrix_data,\n",
    "        dist_supplier_data,\n",
    "    ) = read_input_irp(instance_file)\n",
    "\n",
    "    model = OptModel()\n",
    "    # Quantity of product delivered at each discrete time instant of\n",
    "    # the planning time horizon to each customer\n",
    "    delivery = [\n",
    "        [model.float(0, capacity) for _ in range(nb_customers)]\n",
    "        for _ in range(horizon_length)\n",
    "    ]\n",
    "    # Sequence of customers visited at each discrete time instant of\n",
    "    # the planning time horizon\n",
    "    route = [model.list(nb_customers) for _ in range(horizon_length)]\n",
    "\n",
    "    # Customers receive products only if they are visited at least once\n",
    "    is_delivered = [\n",
    "        [model.contains(route[t], customer) for customer in range(nb_customers)] for t in range(horizon_length)\n",
    "    ]\n",
    "    # Create Hexaly arrays to be able to access them with an \"at\" operator\n",
    "    dist_matrix = model.array(dist_matrix_data)\n",
    "    dist_supplier = model.array(dist_supplier_data)\n",
    "    dist_routes = []\n",
    "    for t in range(horizon_length):\n",
    "        sequence = route[t]\n",
    "        count = model.count(sequence)\n",
    "        # Distance traveled at instant t\n",
    "        dist_lambda = model.lambda_function(lambda position: dist_matrix[sequence[position - 1], sequence[position]])\n",
    "        dist_routes.append(\n",
    "            model.iif(\n",
    "                count > 0,\n",
    "                dist_supplier[sequence[0]]\n",
    "                + model.sum(model.range(1, count), dist_lambda)\n",
    "                + dist_supplier[sequence[count - 1]],\n",
    "                0,\n",
    "            )\n",
    "        )\n",
    "\n",
    "    # Stockout constraints at the supplier\n",
    "    inventory_supplier = [start_level_supplier]\n",
    "    for t in range(horizon_length):\n",
    "        total_delivery = model.sum(delivery[t])\n",
    "        inventory_supplier.append(inventory_supplier[t] - total_delivery + production_rate_supplier)\n",
    "        model.constraint(\n",
    "            inventory_supplier[t] >= total_delivery,\n",
    "        )\n",
    "\n",
    "    # Stockout constraints at the customers\n",
    "    inventory = [[start_level[customer]] for customer in range(nb_customers)]\n",
    "    for customer in range(nb_customers):\n",
    "        for t in range(horizon_length):\n",
    "            next_inventory = inventory[customer][t] + delivery[t][customer] - demand_rate[customer]\n",
    "            inventory[customer].append(next_inventory)\n",
    "            model.constraint(\n",
    "                next_inventory >= 0,\n",
    "            )\n",
    "\n",
    "    for t in range(horizon_length):\n",
    "        model.constraint(\n",
    "            model.sum(delivery[t]) <= capacity,\n",
    "        )\n",
    "        for customer in range(nb_customers):\n",
    "            model.constraint(\n",
    "                delivery[t][customer] <= max_level[customer] - inventory[customer][t],\n",
    "            )\n",
    "            model.constraint(\n",
    "                delivery[t][customer] <= max_level[customer] * is_delivered[t][customer],\n",
    "            )\n",
    "\n",
    "    total_cost_inventory_supplier = holding_cost_supplier * model.sum(*inventory_supplier)\n",
    "    total_cost_inventory = model.sum(\n",
    "        *(\n",
    "            holding_cost[customer] * inventory[customer][t]\n",
    "            for customer in range(nb_customers)\n",
    "            for t in range(horizon_length + 1)\n",
    "        )\n",
    "    )\n",
    "    total_cost_route = model.sum(dist_routes)\n",
    "    objective = total_cost_inventory_supplier + total_cost_inventory + total_cost_route\n",
    "    model.minimize(objective)\n",
    "\n",
    "    solution = solve(model, time_limit_s=float(time_limit))\n",
    "    if not solution.feasible:\n",
    "        # With the original zero-delivery/empty-route defaults, the current\n",
    "        # OptAgent search may not find a feasible coupled list/float solution\n",
    "        # within a short run. Keep the model unchanged and report the status.\n",
    "        print(f\"Status = {solution.feasible}; no feasible solution found\")\n",
    "        return solution\n",
    "\n",
    "    values = {'objective': objective.value, 'route_cost': total_cost_route.value, **{f'route_{t}': route[t].value for t in range(horizon_length)}}\n",
    "    lines = [\n",
    "        f\"Status = {solution.feasible}; Objective = {values['objective']}; Route cost = {values['route_cost']}\"\n",
    "    ]\n",
    "    for t in range(horizon_length):\n",
    "        customers = \" \".join(str(customer + 1) for customer in values[f\"route_{t}\"])\n",
    "        lines.append(f\"Period {t + 1}: {customers}\")\n",
    "    result_text = \"\\n\".join(lines)\n",
    "    print(result_text)\n",
    "\n",
    "    if output_file is not None:\n",
    "        output_lines = [str(int(round(values[\"route_cost\"])))]\n",
    "        output_lines.extend(\n",
    "            \" \".join(str(customer + 1) for customer in values[f\"route_{t}\"]) for t in range(horizon_length)\n",
    "        )\n",
    "        Path(output_file).write_text(\"\\n\".join(output_lines) + \"\\n\", encoding=\"utf-8\")\n",
    "    return solution\n",
    "\n",
    "\n",
    "# The input files follow the \"Archetti\" format\n",
    "def read_input_irp(filename):\n",
    "    file_it = iter(read_elem(filename))\n",
    "\n",
    "    nb_customers = int(next(file_it)) - 1\n",
    "    horizon_length = int(next(file_it))\n",
    "    capacity = int(next(file_it))\n",
    "\n",
    "    x_coord = [None] * nb_customers\n",
    "    y_coord = [None] * nb_customers\n",
    "    start_level = [None] * nb_customers\n",
    "    max_level = [None] * nb_customers\n",
    "    min_level = [None] * nb_customers\n",
    "    demand_rate = [None] * nb_customers\n",
    "    holding_cost = [None] * nb_customers\n",
    "\n",
    "    next(file_it)\n",
    "    x_coord_supplier = float(next(file_it))\n",
    "    y_coord_supplier = float(next(file_it))\n",
    "    start_level_supplier = int(next(file_it))\n",
    "    production_rate_supplier = int(next(file_it))\n",
    "    holding_cost_supplier = float(next(file_it))\n",
    "    for i in range(nb_customers):\n",
    "        next(file_it)\n",
    "        x_coord[i] = float(next(file_it))\n",
    "        y_coord[i] = float(next(file_it))\n",
    "        start_level[i] = int(next(file_it))\n",
    "        max_level[i] = int(next(file_it))\n",
    "        min_level[i] = int(next(file_it))\n",
    "        demand_rate[i] = int(next(file_it))\n",
    "        holding_cost[i] = float(next(file_it))\n",
    "\n",
    "    distance_matrix = compute_distance_matrix(x_coord, y_coord)\n",
    "    distance_supplier = compute_distance_supplier(x_coord_supplier, y_coord_supplier, x_coord, y_coord)\n",
    "\n",
    "    return (\n",
    "        nb_customers,\n",
    "        horizon_length,\n",
    "        capacity,\n",
    "        start_level_supplier,\n",
    "        production_rate_supplier,\n",
    "        holding_cost_supplier,\n",
    "        start_level,\n",
    "        max_level,\n",
    "        demand_rate,\n",
    "        holding_cost,\n",
    "        distance_matrix,\n",
    "        distance_supplier,\n",
    "    )\n",
    "\n",
    "\n",
    "# Compute the distance matrix\n",
    "def compute_distance_matrix(x_coord, y_coord):\n",
    "    nb_customers = len(x_coord)\n",
    "    distance_matrix = [[None for _ in range(nb_customers)] for _ in range(nb_customers)]\n",
    "    for i in range(nb_customers):\n",
    "        distance_matrix[i][i] = 0\n",
    "        for j in range(nb_customers):\n",
    "            dist = compute_dist(x_coord[i], x_coord[j], y_coord[i], y_coord[j])\n",
    "            distance_matrix[i][j] = dist\n",
    "            distance_matrix[j][i] = dist\n",
    "    return distance_matrix\n",
    "\n",
    "\n",
    "# Compute the distances to the supplier\n",
    "def compute_distance_supplier(x_coord_supplier, y_coord_supplier, x_coord, y_coord):\n",
    "    nb_customers = len(x_coord)\n",
    "    distance_supplier = [None] * nb_customers\n",
    "    for i in range(nb_customers):\n",
    "        dist = compute_dist(x_coord_supplier, x_coord[i], y_coord_supplier, y_coord[i])\n",
    "        distance_supplier[i] = dist\n",
    "    return distance_supplier\n",
    "\n",
    "\n",
    "def compute_dist(xi, xj, yi, yj):\n",
    "    return round(math.sqrt(math.pow(xi - xj, 2) + math.pow(yi - yj, 2)))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "63e43d3d",
   "metadata": {},
   "outputs": [],
   "source": [
    "INSTANCE_DIR = Path.cwd() / \"instances\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0ade11a5",
   "metadata": {},
   "outputs": [],
   "source": [
    "solution = main(INSTANCE_DIR / \"abs1n5.dat\", time_limit=5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e2066120",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": ".venv",
   "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.14.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}