{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 二次分配问题 (QAP)\n",
    "\n",
    "**类别：** 选址\n",
    "\n",
    "使用 OptAgent 的 Python 接口描述变量、约束与目标。\n",
    "\n",
    "问题与原始示例来源：[Hexaly Code Templates](https://www.hexaly.com/templates/quadratic-assignment-problem-qap)。\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 问题描述\n",
    "\n",
    "**在 Quadratic Assignment Problem (QAP) 中**，需要将 n 个设施分配到 n 个位置。问题数据包含每对位置之间的距离以及每对设施之间的流量（或权重），即它们之间运输的物料量。问题目标是将每个设施分配到一个位置，使距离与对应流量乘积之和最小。该问题类似于指派问题（Assignment Problem），不同之处在于其目标函数由二次不等式表示，因此得名。更多细节，请参阅 [QAPLIB 网页](https://coral.ise.lehigh.edu/data-sets/qaplib/qaplib-problem-instances-and-solutions/)。\n",
    "\n",
    "### 学习要点\n",
    "\n",
    "- 添加 [列表决策变量](https://optagent.pages.dev/guide/modeling/) 来建模设施的排列\n",
    "- 使用 [‘count’ operator](https://optagent.pages.dev/guide/modeling/) 约束 list 中的元素数量\n",
    "- 使用 [‘at’ operator](https://optagent.pages.dev/guide/modeling/) 访问 list 中的元素以定义目标函数\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 数据\n",
    "\n",
    "我们提供来自 [QAPLIB](http://anjos.mgi.polymtl.ca/qaplib/) 的实例。数据文件的格式如下：\n",
    "\n",
    "- 点的数量\n",
    "- 矩阵 A：每对位置之间的距离\n",
    "- 矩阵 B：每对设施之间的流量\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 建模思路\n",
    "\n",
    "Quadratic Assignment Problem (QAP) 的 OptAgent 模型仅使用一个 [列表决策变量](https://optagent.pages.dev/guide/modeling/)。该 list 表示设施的一个排列：list 中第 i 个位置的元素表示分配到位置 i 的设施索引。\n",
    "\n",
    "使用 **count** 算子，我们将 list 的大小约束为等于设施总数。从而确保所有设施都被分配到一个位置。\n",
    "\n",
    "最后，我们计算目标函数。使用 **at** 算子，我们可以访问分配到位置 i 的设施（p[i]，其中 p 是 list 变量）。然后我们可以轻松地获取分配到位置 i 和 j 的设施之间的流量（B[p[i]][p[j]]）。将该数量乘以位置 i 和 j 之间的距离（A[i][j]），即可得到与这两个设施相关联的成本。\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Python 实现\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Copyright (c) Hexaly. Permission is hereby granted to use, copy,\n",
    "# and modify this code for applications developed with Hexaly.\n",
    "from pathlib import Path\n",
    "\n",
    "from optagent import OptModel, solve\n",
    "\n",
    "\n",
    "def read_integers(filename):\n",
    "    return [int(elem) for elem in Path(filename).read_text(encoding=\"utf-8\").split()]\n",
    "\n",
    "\n",
    "def read_instance(filename):\n",
    "    file_it = iter(read_integers(filename))\n",
    "\n",
    "    # Number of points\n",
    "    n = next(file_it)\n",
    "\n",
    "    # Distance between locations\n",
    "    distances = [[next(file_it) for _ in range(n)] for _ in range(n)]\n",
    "    # Flow between factories\n",
    "    flows = [[next(file_it) for _ in range(n)] for _ in range(n)]\n",
    "    return n, distances, flows\n",
    "\n",
    "\n",
    "def main(instance_file, output_file=None, time_limit=30):\n",
    "    n, distances, flows = read_instance(instance_file)\n",
    "    model = OptModel()\n",
    "\n",
    "    # Permutation such that p[i] is the facility on the location i\n",
    "    permutation = model.list(n)\n",
    "\n",
    "    # The list must be complete\n",
    "    model.constraint(model.count(permutation) == n)\n",
    "\n",
    "    # Create the flow matrix as a model array for expression-based indexing.\n",
    "    flow_array = model.array(flows)\n",
    "\n",
    "    # Minimize the sum of product distance*flow\n",
    "    objective = model.sum(\n",
    "        distances[i][j] * flow_array[permutation[i], permutation[j]]\n",
    "        for i in range(n)\n",
    "        for j in range(n)\n",
    "    )\n",
    "    model.minimize(objective)\n",
    "\n",
    "    solution = solve(model, time_limit_s=float(time_limit))\n",
    "    print(\n",
    "        f\"Facilities = {n}; Objective = {objective.value}; \"\n",
    "        f\"Status = {solution.feasible}\"\n",
    "    )\n",
    "\n",
    "    #\n",
    "    # Write the solution in a file with the following format:\n",
    "    #  - n objValue\n",
    "    #  - permutation p\n",
    "    #\n",
    "    if output_file is not None:\n",
    "        result_text = (\n",
    "            f\"{n} {objective.value}\\n\"\n",
    "            + \" \".join(str(facility) for facility in permutation.value)\n",
    "            + \"\\n\"\n",
    "        )\n",
    "        Path(output_file).write_text(result_text, encoding=\"utf-8\")\n",
    "    return solution\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 实例调用\n",
    "\n",
    "下面使用仓库提供的 QAPLIB 实例运行模型。\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "INSTANCE_DIR = Path.cwd() / \"instances\"\n",
    "print(\"Instances:\", INSTANCE_DIR)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "solution_qap = main(INSTANCE_DIR / \"nug30.dat\", time_limit=1)\n",
    "solution_qap.feasible\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python",
   "version": "3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}