{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 设施选址问题 (FLP)\n",
    "\n",
    "**类别：** 选址\n",
    "\n",
    "使用 OptAgent 的 Python 接口描述变量、约束与目标。\n",
    "\n",
    "问题与原始示例来源：[Hexaly Code Templates](https://www.hexaly.com/templates/facility-location-problem-flp)。\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 问题描述\n",
    "\n",
    "**设施选址问题 (FLP)** 的定义如下。给定一组地点以及每对地点之间的运输成本,从中选取 p 个地点作为设施,以最小化运输成本。一个地点的运输成本等于其到最近设施的距离。因此,目标是为设施提供最优布局以最小化运输成本。该问题也称为 P-Median 问题。\n",
    "\n",
    "### 学到的要点\n",
    "\n",
    "- 添加布尔决策变量来建模一个地点是否被选为设施\n",
    "- 使用 `iif` 和 `min` 算子计算每个地点到最近设施的运输成本\n",
    "- 使用 OptAgent 的 `ModelBuilder` 建模、`solve` 求解并读取解中的表达式值\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 数据\n",
    "\n",
    "我们提供的数据文件来自 [OR-LIB](http://people.brunel.ac.uk/~mastjjb/jeb/orlib/pmedinfo.html)。其格式如下:\n",
    "\n",
    "- 地点数量\n",
    "- 原始实例中的边数\n",
    "- 要选择的设施数量\n",
    "- 距离矩阵,使用 Floyd 算法从原始文件计算得出。\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 建模思路\n",
    "\n",
    "设施选址问题 (FLP) 的 OptAgent 模型使用布尔决策变量来表示每个地点是否被选为设施。我们将这些布尔值的和约束为不超过 p,以确保设施数量最多为 p。\n",
    "\n",
    "然后我们需要计算运输成本。注意,在模型中无需为每个地点定义其最近设施,只需计算该地点与其最近设施之间的距离即可。利用三元算子,我们可以计算地点 i 与设施 j 之间的运输成本:若 j 为设施则等于 i 与 j 的距离,否则为有限的大惩罚值 `2 * max_distance`。利用对所有可能设施的 **min** 算子,我们可以计算地点 i 与其最近设施之间的距离,对所有 i 均如此。\n",
    "\n",
    "最后,需要最小化的目标即为这些运输成本之和。\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Python 实现\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pathlib import Path\n",
    "\n",
    "from optagent import OptModel, solve\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "def read_instance(instance_file):\n",
    "    file_it = iter(int(value) for value in Path(instance_file).read_text().split())\n",
    "    n = next(file_it)\n",
    "    next(file_it)  # Skip the number of edges.\n",
    "    p = next(file_it)\n",
    "    distances = [[next(file_it) for _ in range(n)] for _ in range(n)]\n",
    "    max_distance = max(max(row) for row in distances)\n",
    "    return n, p, max_distance, distances\n",
    "\n",
    "\n",
    "def main(instance_file, output_file=None, time_limit=20):\n",
    "    n, p, max_distance, distances = read_instance(instance_file)\n",
    "\n",
    "    model = OptModel()\n",
    "    # One variable for each location: 1 if facility, 0 otherwise\n",
    "    is_open = [model.bool() for _ in range(n)]\n",
    "    opened_locations = model.sum(is_open)\n",
    "    # Keep the original Hexaly formulation: select no more than p facilities.\n",
    "    # The finite 2 * max_distance penalty makes fewer (or no) open facilities feasible.\n",
    "    # Nonnegative distances make extra facilities non-worsening, but a time-limited\n",
    "    # heuristic solve may return such a feasible solution before finding a better one.\n",
    "    model.constraint(opened_locations <= p)\n",
    "\n",
    "    # Costs between location i and j is w[i][j] if j is a facility or 2 * wmax if not\n",
    "    location_costs = []\n",
    "    for location in range(n):\n",
    "        # 这里相当于直接根据布尔变量取值来进行不同罚分的设定\n",
    "        candidate_costs = [\n",
    "            model.iif(is_open[facility], distances[location][facility], 2 * max_distance)\n",
    "            for facility in range(n)\n",
    "        ]\n",
    "        # Cost between location i and the closest facility\n",
    "        location_costs.append(model.min(candidate_costs))\n",
    "\n",
    "    total_cost = model.sum(location_costs)\n",
    "    model.minimize(total_cost)\n",
    "\n",
    "    solution = solve(model, time_limit_s=float(time_limit))\n",
    "    values = {'total_cost': total_cost.value, **{f'is_open_{i}': variable.value for i, variable in enumerate(is_open)}}\n",
    "    facilities = [i for i in range(n) if values[f\"is_open_{i}\"]]\n",
    "    result_text = (\n",
    "        f\"Total cost = {values['total_cost']}; Facilities = {facilities}; \"\n",
    "        f\"Status = {solution.feasible}\"\n",
    "    )\n",
    "    print(result_text)\n",
    "    if output_file is not None:\n",
    "        Path(output_file).write_text(\n",
    "            f\"{values['total_cost']}\\n{' '.join(map(str, facilities))}\\n\",\n",
    "            encoding=\"utf-8\",\n",
    "        )\n",
    "    return solution\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 本地运行\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "INSTANCE_DIR = Path.cwd() / \"instances\"\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 在当前建模下，什么都不做也是可行解\n",
    "solution = main(INSTANCE_DIR / \"pmed1.in\", time_limit=5)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b1568a20",
   "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
}