{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "7fb27b941602401d91542211134fc71a",
   "metadata": {},
   "source": [
    "# 带约束的露天矿坑极限问题 (CPIT)\n",
    "\n",
    "**类别：** 装箱\n",
    "\n",
    "使用 OptAgent 的 Python 接口描述变量、约束与目标。\n",
    "\n",
    "问题与原始示例来源：[Hexaly Code Templates](https://www.hexaly.com/templates/constrained-pit-limit-cpit)。\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "acae54e37e7d407bbb7b55eff062a284",
   "metadata": {},
   "source": [
    "## 问题描述\n",
    "\n",
    "在带约束的露天矿坑极限问题 (CPIT) 中,我们考虑一组可在若干时间段内从矿坑中开采的块。开采一个块需要消耗一定数量的资源。对于每种资源,在某个时间段内的资源消耗总和不能超过某个上限。每个块必须在其所有前驱块被开采的同一时间段或更晚的时间段开采。一个块的开采利润取决于该块和开采时间段。目标函数是最大化从块开采中获得的利润。\n",
    "\n",
    "### 学到的要点\n",
    "\n",
    "- 添加 [集合决策变量](https://optagent.pages.dev/guide/modeling/) 来建模每个时间段的开采情况\n",
    "- 通过 [`find`](https://optagent.pages.dev/api/model-reference/#optagent.OptModel.find) 算子获取每个块开采时间段的索引\n",
    "- 定义一个 [lambda 函数](https://optagent.pages.dev/guide/modeling/) 来对每个时间段的利润求和\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9a63283cbaf04dbcab1f6479b197f3a8",
   "metadata": {},
   "source": [
    "## 数据\n",
    "\n",
    "所提供的实例来自 [Espinoza 等人](https://link.springer.com/article/10.1007/s10479-012-1258-3),其格式如下:\n",
    "\n",
    "- 一个 .cpit 文件存储与带约束的露天矿坑极限问题 (CPIT) 相关的数据:\n",
    "\n",
    "- 块的数量\n",
    "- 时间段的数量\n",
    "- 资源的数量\n",
    "- 折现率(用于计算块开采的利润)\n",
    "- 对每种资源和每个时间段,该时间段内该资源消耗的上下界\n",
    "- 对每个块,其利润\n",
    "- 对每个块和每种资源,该块的开采对该资源的消耗(仅当非零时写入文件)\n",
    "- 一个 .prec 文件存储每个块的前驱块。\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8dd0d8092fe74a7c96281538738b07e2",
   "metadata": {},
   "source": [
    "## 建模思路\n",
    "\n",
    "带约束的露天矿坑极限问题 (CPIT) 的 OptAgent 模型使用 [集合决策变量](https://optagent.pages.dev/guide/modeling/)。每个集合表示在某个时间段内开采的块。我们添加一个虚拟时间段来表示保持未开采的块。借助 [`partition`](https://optagent.pages.dev/api/model-reference/#optagent.OptModel.partition) 算子,我们确保每个块最多被开采一次。\n",
    "\n",
    "对于给定类型的资源,一个 lambda 函数对所有被开采的块应用 [`sum`](https://optagent.pages.dev/api/model-reference/#optagent.OptModel.sum) 算子来计算该时间段内该资源的总消耗量。注意,该求和中项的数量在搜索过程中是变化的,集合的大小也随之变化。然后我们可以根据为问题定义的上下界对该数量施加约束。\n",
    "\n",
    "然后,使用 [`find`](https://optagent.pages.dev/api/model-reference/#optagent.OptModel.find) 算子,我们可以获取每个块被选定的开采时间段的索引。这使我们能够写出块之间的优先关系约束:每个块必须在与其前驱块的同一时间段或更晚的时间段被开采。\n",
    "\n",
    "一个时间段的利润通过一个 lambda 函数对该时间段内被开采的块应用 [`sum`](https://optagent.pages.dev/api/model-reference/#optagent.OptModel.sum) 算子来计算。同样地,该求和中项的数量在搜索过程中是变化的,集合的大小也随之变化。我们最大化所有时间段的总利润。\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "72eea5119410473aa328ad9291626812",
   "metadata": {},
   "source": [
    "## Python 实现\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8edb47106e1a46a883d545849b8ab81b",
   "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",
    "    lines = Path(instance_file).read_text().splitlines()\n",
    "    index = 2\n",
    "    nb_blocks = int(lines[index].split()[1])\n",
    "    index += 1\n",
    "    nb_periods = int(lines[index].split()[1])\n",
    "    index += 1\n",
    "    nb_resources = int(lines[index].split()[1])\n",
    "    index += 1\n",
    "    discount_rate = float(lines[index].split()[1])\n",
    "    index += 2\n",
    "    lower_bounds = [[None] * nb_periods for _ in range(nb_resources)]\n",
    "    upper_bounds = [[None] * nb_periods for _ in range(nb_resources)]\n",
    "    for _ in range(nb_resources * nb_periods):\n",
    "        resource, period, relation, value = (\n",
    "            lines[index].split()[0],\n",
    "            lines[index].split()[1],\n",
    "            lines[index].split()[2],\n",
    "            lines[index].split()[-1],\n",
    "        )\n",
    "        index += 1\n",
    "        target = lower_bounds if relation == \"G\" else upper_bounds\n",
    "        target[int(resource)][int(period)] = float(value)\n",
    "    index += 1\n",
    "    profit = [0.0] * nb_blocks\n",
    "    for _ in range(nb_blocks):\n",
    "        block, value = lines[index].split()\n",
    "        index += 1\n",
    "        profit[int(block)] = float(value)\n",
    "    index += 2\n",
    "    resource_use = [[0.0] * nb_blocks for _ in range(nb_resources)]\n",
    "    while lines[index] != \"EOF\":\n",
    "        block, resource, value = lines[index].split()\n",
    "        resource_use[int(resource)][int(block)] = float(value)\n",
    "        index += 1\n",
    "    precedence = {}\n",
    "    prec_file = Path(instance_file).parent / \"prec\" / f\"{Path(instance_file).stem}.prec\"\n",
    "    for line in prec_file.read_text().splitlines():\n",
    "        values = [int(value) for value in line.split()]\n",
    "        if values[1]:\n",
    "            precedence[values[0]] = values[2:]\n",
    "    discounted_profit = [[value / (1 + discount_rate) ** period for value in profit] for period in range(nb_periods)]\n",
    "    return nb_blocks, nb_periods, lower_bounds, upper_bounds, discounted_profit, resource_use, precedence\n",
    "\n",
    "\n",
    "def main(instance_file, time_limit=60):\n",
    "    nb_blocks, nb_periods, lower_bounds, upper_bounds, discounted_profit, resource_use_data, precedence = read_instance(\n",
    "        instance_file\n",
    "    )\n",
    "    model = OptModel()\n",
    "    # 最后一个集合( period_nb_periods )是 虚拟时间段 ,表示\"该矿块选择不开采\"\n",
    "    period_sets = [\n",
    "        model.set(nb_blocks)\n",
    "        for period in range(nb_periods + 1)\n",
    "    ]\n",
    "    periods = model.array(period_sets)\n",
    "    model.constraint(model.partition(period_sets))\n",
    "    resource_use = [model.array(values) for values in resource_use_data]\n",
    "    for resource, (lower, upper) in enumerate(zip(lower_bounds, upper_bounds)):\n",
    "        for period in range(nb_periods):\n",
    "            consumption = model.sum(\n",
    "                period_sets[period], model.lambda_function(lambda block: resource_use[resource][block])\n",
    "            )\n",
    "            if lower[period] is not None:\n",
    "                model.constraint(consumption >= lower[period])\n",
    "            if upper[period] is not None:\n",
    "                model.constraint(consumption <= upper[period])\n",
    "    extraction_period = [model.find(periods, block) for block in range(nb_blocks)]\n",
    "    for block, predecessors in precedence.items():\n",
    "        for predecessor in predecessors:\n",
    "            model.constraint(extraction_period[block] >= extraction_period[predecessor])\n",
    "    profit_per_period = []\n",
    "    for period in range(nb_periods):\n",
    "        profit = model.array(discounted_profit[period])\n",
    "        profit_per_period.append(model.sum(period_sets[period], model.lambda_function(lambda block: profit[block])))\n",
    "    objective = model.sum(profit_per_period)\n",
    "    model.maximize(objective)\n",
    "    solution = solve(model, time_limit_s=float(time_limit))\n",
    "    values = {'discounted_profit': objective.value}\n",
    "    print(f\"Discounted profit = {values['discounted_profit']}; Status = {solution.feasible}\")\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 / \"zuck_small.cpit\", time_limit=5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "769cd5cc",
   "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
}