{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "7fb27b941602401d91542211134fc71a",
   "metadata": {},
   "source": [
    "# 集群车辆路径问题 (cluVRP)\n",
    "\n",
    "**类别：** 路径优化\n",
    "\n",
    "使用 OptAgent 的 Python 接口描述变量、约束与目标。\n",
    "\n",
    "问题与原始示例来源：[Hexaly Code Templates](https://www.hexaly.com/templates/clustered-vehicle-routing-cluvrp)。\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "acae54e37e7d407bbb7b55eff062a284",
   "metadata": {},
   "source": [
    "## 问题描述\n",
    "\n",
    "**在集群车辆路径问题 (cluVRP)** 中,一组具有相同载货能力的运输车辆必须为已知需求的客户集群提供单一商品的配送服务。集群在事先已知,由彼此相邻的客户组成。所有客户必须被访问恰好一次。车辆从一个共同的配送中心出发并最终返回配送中心,每辆车服务的总需求量不能超过其载货能力。同一集群内的客户必须被一起服务。换句话说,当车辆访问某集群中的一个客户时,必须在离开该集群之前访问该集群中的所有其他客户。目标是最小化总行驶距离。\n",
    "\n",
    "### 学到的要点\n",
    "\n",
    "- 添加 [列表决策变量](https://optagent.pages.dev/guide/modeling/) 来建模每辆卡车访问集群的顺序以及每个集群内客户的访问顺序\n",
    "- 添加 [partition](https://optagent.pages.dev/guide/modeling/) 约束以确保所有集群都被访问\n",
    "- 定义 [lambda 函数](https://optagent.pages.dev/guide/modeling/) 来计算行驶距离\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9a63283cbaf04dbcab1f6479b197f3a8",
   "metadata": {},
   "source": [
    "## 数据\n",
    "\n",
    "我们提供的**集群车辆路径问题 (cluVRP)** 实例来自论文 [Exact Algorithms for the Clustered Vehicle Routing Problem](https://www.researchgate.net/publication/260941176_Exact_Algorithms_for_the_Clustered_Vehicle_Routing_Problem/)。它们遵循 [TSPLib 格式](http://comopt.ifi.uni-heidelberg.de/software/TSPLIB95/DOC.PS),具体如下:\n",
    "\n",
    "- 节点数量跟在关键字 DIMENSION 之后(由于有一个仓库,客户数量等于节点数减 1)。\n",
    "- 卡车载货能力跟在关键字 CAPACITY 之后。\n",
    "- 卡车数量跟在关键字 VEHICLES 之后。\n",
    "- 集群数量跟在关键字 GVRP_CAPACITY 之后。\n",
    "- 关键字 NODE_COORD_SECTION 之后:对每个节点,给出其 ID 和 x、y 坐标。\n",
    "- 关键字 GVRP_SET_SECTION 之后:对每个集群,给出其 ID 以及属于该集群的节点(以值 -1 结束)。\n",
    "- 关键字 DEMAND_SECTION 之后:对每个集群,给出其 ID 和总需求(该集群内客户需求的总和)。\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8dd0d8092fe74a7c96281538738b07e2",
   "metadata": {},
   "source": [
    "## 建模思路\n",
    "\n",
    "集群车辆路径问题 (cluVRP) 的 OptAgent 模型使用列表决策变量。对于每辆卡车,我们定义一个列表变量来表示它访问的集群序列(`truckSequences`)。通过在所有这些列表上使用 [**partition**](https://optagent.pages.dev/guide/modeling/) 约束,我们确保每个集群恰好由一辆卡车服务。然后我们使用第二组列表决策变量来建模卡车访问每个集群内客户的顺序(`clustersSequences`)。为了确保我们访问所有客户,我们将这些列表的大小约束为等于集群中的客户数。\n",
    "\n",
    "每辆卡车交付的总数量通过 [**lambda 函数**](https://optagent.pages.dev/guide/modeling/) 对所有访问过的集群应用 `sum` 算子来计算。注意,该求和中项的数量以及列表的大小在搜索过程中是变化的。我们将该数量约束为小于卡车的载货能力。\n",
    "\n",
    "然后我们计算每个集群内部的行驶距离。使用另一个 [**lambda 函数**](https://optagent.pages.dev/guide/modeling/),我们沿路径将从一个客户到下一个客户的距离累加起来。我们还记录每个集群中访问的第一个和最后一个节点(分别为 `initialNodes` 和 `endNodes`)。\n",
    "\n",
    "然后我们可以通过累加以下各项来计算每辆卡车的行驶距离:\n",
    "\n",
    "- 在每个访问过的集群内行驶的距离\n",
    "- 从一个集群的最后一个客户到下一个集群的第一个客户的距离之和\n",
    "- 从配送中心到第一个集群的第一个客户的距离\n",
    "- 从最后一个集群的最后一个客户返回配送中心的距离。\n",
    "\n",
    "最后,我们最小化总行驶距离。\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "72eea5119410473aa328ad9291626812",
   "metadata": {},
   "source": [
    "## Python 实现\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8edb47106e1a46a883d545849b8ab81b",
   "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 read_input_cvrp(filename):\n",
    "    file_it = iter(read_elem(filename))\n",
    "\n",
    "    while True:\n",
    "        token = next(file_it)\n",
    "        if token == \"DIMENSION:\":\n",
    "            nb_nodes = int(next(file_it))\n",
    "            nb_customers = nb_nodes - 1\n",
    "        elif token == \"VEHICLES:\":\n",
    "            nb_trucks = int(next(file_it))\n",
    "        elif token == \"GVRP_SETS:\":\n",
    "            nb_clusters = int(next(file_it))\n",
    "        elif token == \"CAPACITY:\":\n",
    "            truck_capacity = int(next(file_it))\n",
    "        elif token == \"NODE_COORD_SECTION\":\n",
    "            break\n",
    "\n",
    "    customers_x = [None] * nb_customers\n",
    "    customers_y = [None] * nb_customers\n",
    "    depot_x = 0\n",
    "    depot_y = 0\n",
    "    for node in range(nb_nodes):\n",
    "        node_id = int(next(file_it))\n",
    "        if node_id != node + 1:\n",
    "            raise ValueError(\"node identifiers must be consecutive\")\n",
    "        x = int(float(next(file_it)))\n",
    "        y = int(float(next(file_it)))\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",
    "    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",
    "\n",
    "    if next(file_it) != \"GVRP_SET_SECTION\":\n",
    "        raise ValueError(\"missing GVRP_SET_SECTION\")\n",
    "    clusters_data = [None] * nb_clusters\n",
    "    for cluster_index in range(nb_clusters):\n",
    "        if int(next(file_it)) != cluster_index + 1:\n",
    "            raise ValueError(\"cluster identifiers must be consecutive\")\n",
    "        cluster = []\n",
    "        value = int(next(file_it))\n",
    "        while value != -1:\n",
    "            cluster.append(value - 2)\n",
    "            value = int(next(file_it))\n",
    "        clusters_data[cluster_index] = cluster\n",
    "\n",
    "    if next(file_it) != \"DEMAND_SECTION\":\n",
    "        raise ValueError(\"missing DEMAND_SECTION\")\n",
    "    demands = [None] * nb_clusters\n",
    "    for cluster_index in range(nb_clusters):\n",
    "        if int(next(file_it)) != cluster_index + 1:\n",
    "            raise ValueError(\"demand identifiers must be consecutive\")\n",
    "        demands[cluster_index] = int(next(file_it))\n",
    "\n",
    "    return (\n",
    "        nb_customers,\n",
    "        nb_trucks,\n",
    "        nb_clusters,\n",
    "        truck_capacity,\n",
    "        distance_matrix,\n",
    "        distance_depots,\n",
    "        demands,\n",
    "        clusters_data,\n",
    "    )\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",
    "        distance_depots[i] = compute_dist(\n",
    "            depot_x, customers_x[i], depot_y, customers_y[i]\n",
    "        )\n",
    "    return distance_depots\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 main(instance_file, output_file=None, time_limit=20):\n",
    "    (\n",
    "        nb_customers,\n",
    "        nb_trucks,\n",
    "        nb_clusters,\n",
    "        truck_capacity,\n",
    "        matrix_data,\n",
    "        depot_data,\n",
    "        demands_data,\n",
    "        clusters_data,\n",
    "    ) = read_input_cvrp(instance_file)\n",
    "    model = OptModel()\n",
    "    matrix, depot_distances = model.array(matrix_data), model.array(depot_data)\n",
    "    demands, clusters = model.array(demands_data), model.array(clusters_data)\n",
    "\n",
    "    # 决定 cluster 内部的访问顺序（每个cluster内部有哪些客户都是确定的，即数量一致）\n",
    "    cluster_sequences = []\n",
    "    for cluster_id, cluster in enumerate(clusters_data):\n",
    "        sequence = model.list(\n",
    "            len(cluster),\n",
    "        )\n",
    "        model.constraint(sequence.count() == len(cluster))\n",
    "        cluster_sequences.append(sequence)\n",
    "\n",
    "    cluster_distances, initial_nodes, end_nodes = [], [], []\n",
    "    for cluster_id, sequence in enumerate(cluster_sequences):\n",
    "        count = sequence.count()\n",
    "        customers = clusters[cluster_id]\n",
    "        distance_lambda = model.lambda_function(\n",
    "            lambda position: matrix[\n",
    "                customers[sequence[position - 1]],\n",
    "                customers[sequence[position]],\n",
    "            ]\n",
    "        )\n",
    "        cluster_distances.append(\n",
    "            model.sum(model.range(1, count), distance_lambda)\n",
    "        )\n",
    "        initial_nodes.append(customers[sequence[0]])\n",
    "        end_nodes.append(customers[sequence[count - 1]])\n",
    "\n",
    "    # 存放每个集群内部的行驶距离\n",
    "    cluster_distances_array = model.array(cluster_distances)\n",
    "    initial_nodes_array, end_nodes_array = model.array(initial_nodes), model.array(end_nodes)\n",
    "    \n",
    "    # 分配truck所负责的cluster\n",
    "    truck_sequences = [\n",
    "        model.list(nb_clusters)\n",
    "        for truck in range(nb_trucks)\n",
    "    ]\n",
    "    model.constraint(model.partition(truck_sequences))\n",
    "\n",
    "    # 各个部分的距离计算\n",
    "    route_distances = []\n",
    "    for sequence in truck_sequences:\n",
    "        count = sequence.count()\n",
    "        demand_lambda = model.lambda_function(lambda cluster: demands[cluster])\n",
    "        route_quantity = model.sum(sequence, demand_lambda)\n",
    "        model.constraint(route_quantity <= truck_capacity)\n",
    "\n",
    "        route_distance_lambda = model.lambda_function(\n",
    "            lambda position: cluster_distances_array[sequence[position]]\n",
    "            + matrix[\n",
    "                end_nodes_array[sequence[position - 1]],\n",
    "                initial_nodes_array[sequence[position]],\n",
    "            ]\n",
    "        )\n",
    "        route_distances.append(\n",
    "            model.sum(model.range(1, count), route_distance_lambda)\n",
    "            + model.iif(\n",
    "                count > 0,\n",
    "                cluster_distances_array[sequence[0]]\n",
    "                + depot_distances[initial_nodes_array[sequence[0]]]\n",
    "                + depot_distances[end_nodes_array[sequence[count - 1]]],\n",
    "                0,\n",
    "            )\n",
    "        )\n",
    "    total_distance = model.sum(route_distances)\n",
    "    model.minimize(total_distance)\n",
    "    solution = solve(model, time_limit_s=float(time_limit))\n",
    "    values = {'total_distance': total_distance.value, **{f'truck_{truck}': route.value for truck, route in enumerate(truck_sequences)}, **{f'cluster_{cluster}': route.value for cluster, route in enumerate(cluster_sequences)}}\n",
    "    print(f\"Total distance = {values['total_distance']}; Status = {solution.feasible}\")\n",
    "    route_lines = []\n",
    "    for truck in range(nb_trucks):\n",
    "        customers = [\n",
    "            clusters_data[cluster][customer] + 2\n",
    "            for cluster in values[f\"truck_{truck}\"]\n",
    "            for customer in values[f\"cluster_{cluster}\"]\n",
    "        ]\n",
    "        route_line = \" \".join(map(str, customers))\n",
    "        route_lines.append(route_line)\n",
    "        print(f\"Truck {truck + 1}: {route_line}\")\n",
    "    if output_file is not None:\n",
    "        output_text = \"\\n\".join(\n",
    "            [str(int(values[\"total_distance\"])), *route_lines]\n",
    "        )\n",
    "        Path(output_file).write_text(output_text + \"\\n\", encoding=\"utf-8\")\n",
    "    return solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "optagent-example",
   "metadata": {},
   "outputs": [],
   "source": [
    "INSTANCE_DIR = Path.cwd() / \"instances\"\n",
    "solution = main(INSTANCE_DIR / \"A-n32-k5-C11-V2.gvrp\", time_limit=5)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": ".venv",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python",
   "version": "3.14.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
