{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 容量约束车辆路径问题(CVRP)\n",
    "\n",
    "**类别：** 路径优化\n",
    "\n",
    "使用 OptAgent 的 Python 接口描述变量、约束与目标。\n",
    "\n",
    "问题与原始示例来源：[Hexaly Code Templates](https://www.hexaly.com/templates/capacitated-vehicle-routing-problem-cvrp)。\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 问题描述\n",
    "\n",
    "**在容量约束车辆路径问题(Capacitated Vehicle Routing Problem, CVRP)**中,一组具有相同容量的配送车辆必须为具有单一商品已知需求的客户提供服务。车辆从一个共同的配送中心出发并返回。每个客户必须恰好由一辆车服务,且每辆车服务的客户总需求不得超过其容量。目标是最小化总行驶距离,同时最小化所使用的车辆数量。\n",
    "\n",
    "### 学习要点\n",
    "\n",
    "- 添加 列表决策变量 以建模每辆卡车的客户序列\n",
    "- 通过 `partition` 约束保证每个客户恰好由一辆卡车服务\n",
    "- 定义 lambda 函数 来计算行驶距离\n",
    "- 通过 `count` 算子统计使用的卡车数\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 数据\n",
    "\n",
    "所提供的容量约束车辆路径问题(CVRP)算例来自 [Augerat 等人的 Set A 数据集](http://neo.lcc.uma.es/vrp/vrp-instances/capacitated-vrp-instances/)。它们遵循 [TSPLib 格式](http://comopt.ifi.uni-heidelberg.de/software/TSPLIB95/DOC.PS):\n",
    "\n",
    "- 节点数量由关键字 DIMENSION 指定(其中包括一个配送中心,因此客户数量为节点数减 1)。\n",
    "- 卡车容量由关键字 CAPACITY 指定。\n",
    "- 边的类型由 EDGE_WEIGHT_TYPE 指定。请注意,在我们的模型中仅接受 EUC_2D 这一种边类型。\n",
    "- 在关键字 NODE_COORD_SECTION 之后:每个节点的 ID 以及其 x、y 坐标。\n",
    "- 在关键字 DEMAND_SECTION 之后:每个节点的 ID 及其需求。\n",
    "- 配送中心列在关键字 DEPOT_SECTION 之后。请注意,在我们的模型中仅接受一个配送中心。\n",
    "\n",
    "可用车辆数量等于客户数量。\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4d1e725a",
   "metadata": {},
   "source": [
    "## 建模思路\n",
    "\n",
    "容量约束车辆路径问题(CVRP)的 OptAgent 模型使用 列表决策变量。对每辆卡车,我们定义一个列表变量,表示其所访问的客户序列。通过对所有列表施加 partition 约束,我们确保每个客户恰好由一辆卡车服务。\n",
    "\n",
    "当一辆卡车至少访问一个客户时,它才被车队所使用。借助 count 算子(返回列表中的元素数量),我们可以检查每辆卡车是否被使用,从而计算车队中使用的卡车总数。\n",
    "\n",
    "我们可以使用需求数组上的 at 算子来访问序列中每个客户的需求。每辆卡车所配送的总量通过一个 lambda 函数计算,该函数将 `sum` 算子应用于所有被访问的客户。请注意,该 sum 中的项数以及列表的大小在搜索过程中会变化。\n",
    "\n",
    "从一个客户到下一个客户所行驶的距离同样使用二维距离矩阵上的 at 算子来访问。我们使用另一个 lambda 函数 计算每辆卡车的总行驶距离,该函数对序列中相邻客户之间的距离求和。\n",
    "\n",
    "最后依次声明两个最小化目标:先最小化使用的车辆数,再最小化总行驶距离。OptAgent 按目标声明顺序执行字典序优化。\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "970623e3",
   "metadata": {},
   "source": [
    "## Python 实现\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "018bb97e",
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "from pathlib import Path\n",
    "\n",
    "from optagent import OptModel, solve\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "def read_cvrp_instance(filename):\n",
    "    \"\"\"Read an Augerat-format CVRP instance and return a dict of arrays.\"\"\"\n",
    "    tokens = Path(filename).read_text(encoding=\"utf-8\").split()\n",
    "    nb_nodes = 0\n",
    "    truck_capacity = 0\n",
    "    depot_x = 0\n",
    "    depot_y = 0\n",
    "    index = 0\n",
    "    while index < len(tokens):\n",
    "        token = tokens[index]\n",
    "        if token == \"DIMENSION\":\n",
    "            nb_nodes = int(tokens[index + 2])\n",
    "            index += 3\n",
    "            continue\n",
    "        if token == \"CAPACITY\":\n",
    "            truck_capacity = int(tokens[index + 2])\n",
    "            index += 3\n",
    "            continue\n",
    "        if token == \"EDGE_WEIGHT_TYPE\":\n",
    "            if tokens[index + 2] != \"EUC_2D\":\n",
    "                raise ValueError(f\"Only EUC_2D is supported, got {tokens[index + 2]}\")\n",
    "            index += 3\n",
    "            continue\n",
    "        if token == \"NODE_COORD_SECTION\":\n",
    "            index += 1\n",
    "            break\n",
    "        index += 1\n",
    "\n",
    "    nb_customers = nb_nodes - 1\n",
    "    customers_x = [0] * nb_customers\n",
    "    customers_y = [0] * nb_customers\n",
    "    for _ in range(nb_nodes):\n",
    "        node_id = int(tokens[index])\n",
    "        index += 1\n",
    "        x = int(tokens[index])\n",
    "        index += 1\n",
    "        y = int(tokens[index])\n",
    "        index += 1\n",
    "        if node_id == 1:\n",
    "            depot_x = x\n",
    "            depot_y = y\n",
    "        else:\n",
    "            customers_x[node_id - 2] = x\n",
    "            customers_y[node_id - 2] = y\n",
    "\n",
    "    while index < len(tokens) and tokens[index] != \"DEMAND_SECTION\":\n",
    "        index += 1\n",
    "    index += 1\n",
    "    demands = [0] * nb_customers\n",
    "    for _ in range(nb_nodes):\n",
    "        node_id = int(tokens[index])\n",
    "        index += 1\n",
    "        demand = int(tokens[index])\n",
    "        index += 1\n",
    "        if node_id != 1:\n",
    "            demands[node_id - 2] = demand\n",
    "\n",
    "    distance_matrix = compute_distance_matrix(customers_x, customers_y)\n",
    "    distance_depots = compute_distance_depots(\n",
    "        depot_x, depot_y, customers_x, customers_y\n",
    "    )\n",
    "    return {\n",
    "        \"nb_customers\": nb_customers,\n",
    "        \"nb_trucks\": nb_customers,\n",
    "        \"truck_capacity\": truck_capacity,\n",
    "        \"distance_matrix\": distance_matrix,\n",
    "        \"distance_depots\": distance_depots,\n",
    "        \"demands\": demands,\n",
    "    }\n",
    "\n",
    "\n",
    "def compute_dist(xi, xj, yi, yj):\n",
    "    exact_dist = math.sqrt(math.pow(xi - xj, 2) + math.pow(yi - yj, 2))\n",
    "    return int(math.floor(exact_dist + 0.5))\n",
    "\n",
    "\n",
    "def compute_distance_matrix(customers_x, customers_y):\n",
    "    nb_customers = len(customers_x)\n",
    "    distance_matrix = [\n",
    "        [None for _ in range(nb_customers)] for _ in range(nb_customers)\n",
    "    ]\n",
    "    for i in range(nb_customers):\n",
    "        distance_matrix[i][i] = 0\n",
    "        for j in range(nb_customers):\n",
    "            dist = compute_dist(\n",
    "                customers_x[i], customers_x[j], customers_y[i], customers_y[j]\n",
    "            )\n",
    "            distance_matrix[i][j] = dist\n",
    "            distance_matrix[j][i] = dist\n",
    "    return distance_matrix\n",
    "\n",
    "\n",
    "def compute_distance_depots(depot_x, depot_y, customers_x, customers_y):\n",
    "    nb_customers = len(customers_x)\n",
    "    distance_depots = [None] * nb_customers\n",
    "    for i in range(nb_customers):\n",
    "        dist = compute_dist(depot_x, customers_x[i], depot_y, customers_y[i])\n",
    "        distance_depots[i] = dist\n",
    "    return distance_depots\n",
    "\n",
    "\n",
    "def main(instance_file, output_file=None, time_limit=10):\n",
    "    data = read_cvrp_instance(instance_file)\n",
    "    nb_customers = data[\"nb_customers\"]\n",
    "    nb_trucks = data[\"nb_trucks\"]\n",
    "    truck_capacity = data[\"truck_capacity\"]\n",
    "\n",
    "    model = OptModel()\n",
    "    customers_sequences = [\n",
    "        model.list(nb_customers) for k in range(nb_trucks)\n",
    "    ]\n",
    "    model.constraint(model.partition(customers_sequences))\n",
    "\n",
    "    demands = model.array(data[\"demands\"])\n",
    "    dist_matrix = model.array(data[\"distance_matrix\"])\n",
    "    dist_depot = model.array(data[\"distance_depots\"])\n",
    "    trucks_used = [model.count(customers_sequences[k]) > 0 for k in range(nb_trucks)]\n",
    "\n",
    "    dist_routes = []\n",
    "    for k in range(nb_trucks):\n",
    "        sequence = customers_sequences[k]\n",
    "        count = model.count(sequence)\n",
    "\n",
    "        demand_lambda = model.lambda_function(lambda customer: demands[customer])\n",
    "        route_quantity = model.sum(sequence, demand_lambda)\n",
    "        model.constraint(route_quantity <= truck_capacity)\n",
    "\n",
    "        dist_lambda = model.lambda_function(\n",
    "            lambda i: dist_matrix[sequence[i - 1], sequence[i]]\n",
    "        )\n",
    "        dist_routes.append(\n",
    "            model.sum(model.range(1, count), dist_lambda)\n",
    "            + model.iif(\n",
    "                count > 0,\n",
    "                dist_depot[sequence[0]]\n",
    "                + dist_depot[sequence[count - 1]],\n",
    "                0,\n",
    "            )\n",
    "        )\n",
    "\n",
    "    nb_trucks_used = model.sum(*trucks_used)\n",
    "    total_distance = model.sum(*dist_routes)\n",
    "\n",
    "    # Lexicographic objectives, in the same order as the original Hexaly model.\n",
    "    model.minimize(nb_trucks_used)\n",
    "    model.minimize(total_distance)\n",
    "\n",
    "    solution = solve(model, time_limit_s=float(time_limit))\n",
    "    result_values = {'trucks_used': nb_trucks_used.value, 'total_distance': total_distance.value, **{f'truck_{k}': sequence.value for k, sequence in enumerate(customers_sequences)}}\n",
    "\n",
    "    lines = [\n",
    "        f\"Trucks used = {result_values['trucks_used']}; \"\n",
    "        f\"Total distance = {result_values['total_distance']}; \"\n",
    "        f\"Status = {solution.feasible}\"\n",
    "    ]\n",
    "    for truck in range(nb_trucks):\n",
    "        sequence = result_values[f\"truck_{truck}\"]\n",
    "        if sequence:\n",
    "            customers = \" \".join(str(customer + 2) for customer in sequence)\n",
    "            lines.append(f\"Truck {truck + 1}: {customers}\")\n",
    "\n",
    "    result_text = \"\\n\".join(lines)\n",
    "    print(result_text)\n",
    "    if output_file is not None:\n",
    "        Path(output_file).write_text(result_text + \"\\n\", encoding=\"utf-8\")\n",
    "    return solution"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "44207fd2",
   "metadata": {},
   "source": [
    "## 本地运行\n",
    "\n",
    "Notebook 直接调用 `main` 并显式传入实例路径。以下三个代码格相互独立,可以按需要单独运行;调整 `time_limit` 可以控制每个实例的求解时间。\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d8e19edb",
   "metadata": {},
   "outputs": [],
   "source": [
    "INSTANCE_DIR = Path.cwd() / \"instances\"\n",
    "print(\"Instances:\", INSTANCE_DIR)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "96692e94",
   "metadata": {},
   "outputs": [],
   "source": [
    "solution_a_n32_k5 = main(\n",
    "    INSTANCE_DIR / \"A-n32-k5.vrp\",\n",
    "    time_limit=1,\n",
    ")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7f8ab63f",
   "metadata": {},
   "outputs": [],
   "source": [
    "solution_a_n33_k5 = main(\n",
    "    INSTANCE_DIR / \"A-n33-k5.vrp\",\n",
    "    time_limit=10,\n",
    ")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b52c1ddb",
   "metadata": {},
   "outputs": [],
   "source": [
    "solution_a_n33_k6 = main(\n",
    "    INSTANCE_DIR / \"A-n33-k6.vrp\",\n",
    "    time_limit=10,\n",
    ")\n"
   ]
  }
 ],
 "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
}