{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 混流问题\n",
    "\n",
    "**类别：** 非线性优化\n",
    "\n",
    "使用 OptAgent 的 Python 接口描述变量、约束与目标。\n",
    "\n",
    "问题与原始示例来源：[Hexaly Code Templates](https://www.hexaly.com/templates/pooling)。\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 问题描述\n",
    "\n",
    "**在混流问题** 中，具有特定属性浓度的原材料流被混合以生产最终产品。最终产品的属性浓度必须在容差区间内。事实上，这些属性是材料中重要的主要成分。流既可以直接从源节点发送到目标节点，也可以在中间池中混合后再发送到目标节点。此外，必须满足流图以及不同节点的容量约束，且原材料供应不得超过。目标是最大化运营利润。收益来自销售最终产品，而成本对应于需要购买的原材料。\n",
    "\n",
    "换句话说，混流问题结合了最小成本流问题和混合问题的特点。\n",
    "\n",
    "### 建模要点\n",
    "\n",
    "- 使用 `json` 标准库 读取输入文件\n",
    "- 添加 [浮点决策变量](https://optagent.pages.dev/guide/modeling/) 来建模流\n",
    "- 使用 [非线性算子](https://optagent.pages.dev/guide/modeling/)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 数据\n",
    "\n",
    "混流问题实例来自此 [GitHub 仓库](https://github.com/cog-imperial/pooling-network/tree/main/pooling_network/instances/data)，并采用 JSON 格式：\n",
    "\n",
    "- “components”：表征原材料的对象数组，包含：\n",
    "\n",
    "- “name”：源节点的名称\n",
    "- “lower”：最小流出量（未使用，等于 0）\n",
    "- “upper”：最大流出量（供应量）\n",
    "- “price”：一个单位原材料的价格\n",
    "- “quality”：包含原材料中每种属性浓度的对象\n",
    "- “products”：表征最终产品的对象数组，包含：\n",
    "\n",
    "- “name”：目标节点的名称\n",
    "- “lower”：最小流入量（需求量）\n",
    "- “upper”：最大流入量（容量）\n",
    "- “price”：一个单位最终产品的销售价格\n",
    "- “quality_lower”：包含最终产品中每种属性最小浓度的对象\n",
    "- “quality_upper”：包含最终产品中每种属性最大浓度的对象\n",
    "- “pool_size”：包含池容量的对象\n",
    "- “component_to_product_bound”：表征源节点和目标节点之间边的对象数组，包含：\n",
    "\n",
    "- “component”：源节点 c 的名称\n",
    "- “product”：目标节点 p 的名称\n",
    "- “bound”：边 (c, p) 上的最大流量\n",
    "- “cost”：边 (c, p) 上每单位流量的成本\n",
    "- “component_to_pool_fraction”：表征源节点和池节点之间边的对象数组，包含：\n",
    "\n",
    "- “component”：源节点 c 的名称\n",
    "- “pool”：池节点 o 的名称\n",
    "- “fraction”：来自源节点 c 的池 p 流入量的最大比例\n",
    "- “cost”：边 (c, o) 上每单位流量的成本\n",
    "- “pool_to_product_bound”：表征池节点和目标节点之间边的对象数组，包含：\n",
    "\n",
    "- “pool”：池节点 o 的名称\n",
    "- “product”：目标节点 p 的名称\n",
    "- “bound”：边 (o, p) 上的最大流量\n",
    "- “cost”：边 (o, p) 上每单位流量的成本\n",
    "\n",
    "原材料和最终产品的价格可以通过边上的 “price” 或 “cost” 字段表示。“cost” 表示用于 “randstd” 实例，而 “price” 表示用于其他实例。\n",
    "\n",
    "Python 实现使用标准库 `json` 读取实例。\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 建模思路\n",
    "\n",
    "混流问题的 OptAgent 模型对应于问题的 Q 形式，同时使用流量比例和流量量。我们首先声明三个浮点决策变量数组，分别表示：\n",
    "\n",
    "- 从每个源节点 c 到每个目标节点 p 的流量；\n",
    "- 从每个池 o 到每个目标节点 p 的流量；\n",
    "- 每个池 o 中来自每个源节点 c 的流入量比例。\n",
    "\n",
    "通过将后两个量相乘，我们计算出从每个源节点 c 经过每个池 o 到每个目标节点 p 的流量。\n",
    "\n",
    "然后我们可以开始编写约束。首先，我们确保每个池中来自源节点的流入比例之和等于 1。然后，每个目标节点的总流入量必须满足需求同时不超出产品容量。类似地，每个源节点的总流出量不得超过原材料的供应量。最后，每个最终产品的每种属性浓度必须在一定的容差区间内。为了强制执行该约束，我们计算直接来自源节点和来自池的属性量，以及目标节点的总流入量。\n",
    "\n",
    "最后，我们计算目标函数的值。总利润等于销售产品获得的总收入减去初始组件的成本再减去除流量成本。\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Python 实现\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pathlib import Path\n",
    "import json\n",
    "\n",
    "from optagent import OptModel, solve\n",
    "\n",
    "\n",
    "class PoolingInstance:\n",
    "\n",
    "    #\n",
    "    # Read instance data\n",
    "    #\n",
    "    def __init__(self, instance_file):\n",
    "        with Path(instance_file).open(encoding=\"utf-8\") as problem:\n",
    "            problem = json.load(problem)\n",
    "\n",
    "            self.nbComponents = len(problem[\"components\"])\n",
    "            self.nbProducts = len(problem[\"products\"])\n",
    "            self.nbAttributes = len(problem[\"components\"][0][\"quality\"])\n",
    "            self.nbPools = len(problem[\"pool_size\"])\n",
    "\n",
    "            # Components\n",
    "            self.componentPrices = [problem[\"components\"][c][\"price\"]\n",
    "                                    for c in range(self.nbComponents)]\n",
    "            self.componentSupplies = [problem[\"components\"][c][\"upper\"]\n",
    "                                      for c in range(self.nbComponents)]\n",
    "            self.componentQuality = [list(problem[\"components\"][c][\"quality\"].values())\n",
    "                                     for c in range(self.nbComponents)]\n",
    "            self.componentNames = [problem[\"components\"][c][\"name\"]\n",
    "                                   for c in range(self.nbComponents)]\n",
    "\n",
    "            componentsIdx = {}\n",
    "            for c in range(self.nbComponents):\n",
    "                componentsIdx[problem[\"components\"][c][\"name\"]] = c\n",
    "\n",
    "            # Final products (blendings)\n",
    "            self.productPrices = [problem[\"products\"][p][\"price\"]\n",
    "                                  for p in range(self.nbProducts)]\n",
    "            self.productCapacities = [problem[\"products\"][p][\"upper\"]\n",
    "                                      for p in range(self.nbProducts)]\n",
    "            self.demand = [problem[\"products\"][p][\"lower\"]\n",
    "                           for p in range(self.nbProducts)]\n",
    "            self.productNames = [problem[\"products\"][p][\"name\"]\n",
    "                                 for p in range(self.nbProducts)]\n",
    "\n",
    "            productIdx = {}\n",
    "            for p in range(self.nbProducts):\n",
    "                productIdx[problem[\"products\"][p][\"name\"]] = p\n",
    "\n",
    "            self.minTolerance = [[0 for _ in range(self.nbAttributes)]\n",
    "                if (problem[\"products\"][p][\"quality_lower\"] == None)\n",
    "                else list(problem[\"products\"][p][\"quality_lower\"].values())\n",
    "                for p in range(self.nbProducts)]\n",
    "            self.maxTolerance = [list(problem[\"products\"][p][\"quality_upper\"].values())\n",
    "                                 for p in range(self.nbProducts)]\n",
    "\n",
    "            # Intermediate pools\n",
    "            self.poolNames = list(problem[\"pool_size\"].keys())\n",
    "            self.poolCapacities = [problem[\"pool_size\"][o] for o in self.poolNames]\n",
    "            poolIdx = {}\n",
    "            for o in range(self.nbPools):\n",
    "                poolIdx[self.poolNames[o]] = o\n",
    "\n",
    "            # Flow graph\n",
    "\n",
    "            # Edges from the components to the products\n",
    "            self.upperBoundComponentToProduct = [[0 for _ in range(self.nbProducts)]\n",
    "                                                 for _ in range(self.nbComponents)]\n",
    "            self.costComponentToProduct = [[0 for _ in range(self.nbProducts)]\n",
    "                                           for _ in range(self.nbComponents)]\n",
    "            # Edges from the components to the pools\n",
    "            self.upperBoundFractionComponentToPool = [[0 for _ in range(self.nbPools)]\n",
    "                                                      for _ in range(self.nbComponents)]\n",
    "            self.costComponentToPool = [[0 for _ in range(self.nbPools)]\n",
    "                                        for _ in range(self.nbComponents)]\n",
    "            # Edges from the pools to the products\n",
    "            self.upperBoundPoolToProduct = [[0 for _ in range(self.nbProducts)]\n",
    "                                            for _ in range(self.nbPools)]\n",
    "            self.costPoolToProduct = [[0 for _ in range(self.nbProducts)]\n",
    "                                      for _ in range(self.nbPools)]\n",
    "\n",
    "            # Bound and cost on the edges\n",
    "            for edge in problem[\"component_to_product_bound\"]:\n",
    "                self.upperBoundComponentToProduct[componentsIdx[edge[\"component\"]]] \\\n",
    "                    [productIdx[edge[\"product\"]]] = edge[\"bound\"]\n",
    "                if len(edge) > 3:\n",
    "                    self.costComponentToProduct[componentsIdx[edge[\"component\"]]] \\\n",
    "                        [productIdx[edge[\"product\"]]] = edge[\"cost\"]\n",
    "\n",
    "            for edge in problem[\"component_to_pool_fraction\"]:\n",
    "                self.upperBoundFractionComponentToPool[componentsIdx[edge[\"component\"]]] \\\n",
    "                    [poolIdx[edge[\"pool\"]]] = edge[\"fraction\"]\n",
    "                if len(edge) > 3:\n",
    "                    self.costComponentToPool[componentsIdx[edge[\"component\"]]] \\\n",
    "                        [poolIdx[edge[\"pool\"]]] = edge[\"cost\"]\n",
    "\n",
    "            for edge in problem[\"pool_to_product_bound\"]:\n",
    "                self.upperBoundPoolToProduct[poolIdx[edge[\"pool\"]]] \\\n",
    "                    [productIdx[edge[\"product\"]]] = edge[\"bound\"]\n",
    "                if len(edge) > 3:\n",
    "                    self.costPoolToProduct[poolIdx[edge[\"pool\"]]] \\\n",
    "                        [productIdx[edge[\"product\"]]] = edge[\"cost\"]\n",
    "\n",
    "\n",
    "class _OptAgentRunner:\n",
    "    def __init__(self, time_limit):\n",
    "        self.model = OptModel()\n",
    "        self.time_limit = time_limit\n",
    "\n",
    "    def __enter__(self):\n",
    "        return self\n",
    "\n",
    "    def __exit__(self, exc_type, exc_value, traceback):\n",
    "        return False\n",
    "\n",
    "    def solve(self):\n",
    "        return solve(self.model, time_limit_s=float(self.time_limit))\n",
    "\n",
    "\n",
    "def main(instance_file, output_file=None, time_limit=20):\n",
    "    data = PoolingInstance(instance_file)\n",
    "\n",
    "    with _OptAgentRunner(time_limit) as optimizer:\n",
    "        #\n",
    "        # Declare the optimization model\n",
    "        #\n",
    "        model = optimizer.model\n",
    "\n",
    "        # Zero-capacity arcs are constants, shared by every derived expression.\n",
    "        zero = model.constant(0.0)\n",
    "\n",
    "        def bounded_flow(upper):\n",
    "            return model.float(0, upper) if upper > 0 else zero\n",
    "\n",
    "        # Decision variables\n",
    "\n",
    "        # Flow from the components to the products\n",
    "        flowComponentToProduct = [[bounded_flow(data.upperBoundComponentToProduct[c][p])\n",
    "            for p in range(data.nbProducts)] for c in range(data.nbComponents)]\n",
    "\n",
    "        # Fraction of the total flow in pool o coming from the component c\n",
    "        fractionComponentToPool = [[bounded_flow(data.upperBoundFractionComponentToPool[c][o])\n",
    "                                    for o in range(data.nbPools)]\n",
    "                                   for c in range(data.nbComponents)]\n",
    "\n",
    "        # Flow from the pools to the products\n",
    "        flowPoolToProduct = [\n",
    "            [bounded_flow(data.upperBoundPoolToProduct[o][p])\n",
    "             for p in range(data.nbProducts)] for o in range(data.nbPools)]\n",
    "\n",
    "        # Flow from the components to the products and passing by the pools\n",
    "        flowComponentToProductByPool = [\n",
    "            [[(fractionComponentToPool[c][o] * flowPoolToProduct[o][p]\n",
    "               if data.upperBoundFractionComponentToPool[c][o] > 0\n",
    "               and data.upperBoundPoolToProduct[o][p] > 0 else zero)\n",
    "              for p in range(data.nbProducts)] for o in range(data.nbPools)]\n",
    "            for c in range(data.nbComponents)]\n",
    "\n",
    "        # Constraints\n",
    "\n",
    "        # Proportion\n",
    "        for o in range(data.nbPools):\n",
    "            proportion = model.sum(fractionComponentToPool[c][o]\n",
    "                                   for c in range(data.nbComponents))\n",
    "            model.constraint(proportion == 1)\n",
    "\n",
    "        # Component supply\n",
    "        for c in range(data.nbComponents):\n",
    "            flowToProducts = model.sum(flowComponentToProduct[c][p]\n",
    "                                       for p in range(data.nbProducts))\n",
    "            flowToPools = model.sum(flowComponentToProductByPool[c][o][p]\n",
    "                for p in range(data.nbProducts) for o in range(data.nbPools))\n",
    "            totalOutFlow = model.sum(flowToPools, flowToProducts)\n",
    "            model.constraint(totalOutFlow <= data.componentSupplies[c])\n",
    "\n",
    "        # Pool capacity (bounds on edges)\n",
    "        for c in range(data.nbComponents):\n",
    "            for o in range(data.nbPools):\n",
    "                flowComponentToPool = model.sum(flowComponentToProductByPool[c][o][p]\n",
    "                                                for p in range(data.nbProducts))\n",
    "                edgeCapacity = model.prod(data.poolCapacities[o],\n",
    "                                          fractionComponentToPool[c][o])\n",
    "                model.constraint(flowComponentToPool <= edgeCapacity)\n",
    "\n",
    "        # Product capacity\n",
    "        for p in range(data.nbProducts):\n",
    "            flowFromPools = model.sum(flowPoolToProduct[o][p] for o in range(data.nbPools))\n",
    "            flowFromComponents = model.sum(flowComponentToProduct[c][p]\n",
    "                                           for c in range(data.nbComponents))\n",
    "            totalInFlow = model.sum(flowFromComponents, flowFromPools)\n",
    "            model.constraint(totalInFlow <= data.productCapacities[p])\n",
    "            model.constraint(totalInFlow >= data.demand[p])\n",
    "\n",
    "        # Product tolerance\n",
    "        for p in range(data.nbProducts):\n",
    "            for k in range(data.nbAttributes):\n",
    "                # Attribute from the components\n",
    "                attributeFromComponents = model.sum(\n",
    "                    data.componentQuality[c][k] * flowComponentToProduct[c][p]\n",
    "                    for c in range(data.nbComponents)\n",
    "                    if data.upperBoundComponentToProduct[c][p] > 0 and data.componentQuality[c][k] != 0)\n",
    "\n",
    "                # Attribute from the pools\n",
    "                attributeFromPools = model.sum(\n",
    "                    data.componentQuality[c][k] * flowComponentToProductByPool[c][o][p]\n",
    "                    for o in range(data.nbPools) for c in range(data.nbComponents)\n",
    "                    if data.upperBoundFractionComponentToPool[c][o] > 0\n",
    "                    and data.upperBoundPoolToProduct[o][p] > 0 and data.componentQuality[c][k] != 0)\n",
    "\n",
    "                # Total flow in the blending\n",
    "                totalFlowIn = model.sum(flowComponentToProduct[c][p]\n",
    "                                        for c in range(data.nbComponents)) \\\n",
    "                    + model.sum(flowPoolToProduct[o][p] for o in range(data.nbPools))\n",
    "\n",
    "                totalAttributeIn = model.sum(attributeFromComponents, attributeFromPools)\n",
    "                model.constraint(totalAttributeIn >= data.minTolerance[p][k] * totalFlowIn)\n",
    "                model.constraint(totalAttributeIn <= data.maxTolerance[p][k] * totalFlowIn)\n",
    "\n",
    "        # Objective function\n",
    "\n",
    "        # Cost of the flows from the components directly to the products\n",
    "        directFlowCost = model.sum(\n",
    "            data.costComponentToProduct[c][p] * flowComponentToProduct[c][p]\n",
    "            for c in range(data.nbComponents) for p in range(data.nbProducts)\n",
    "            if data.upperBoundComponentToProduct[c][p] > 0 and data.costComponentToProduct[c][p] != 0)\n",
    "\n",
    "        # Cost of the flows from the components to the products passing by the pools\n",
    "        undirectFlowCost = model.sum(\n",
    "            (data.costComponentToPool[c][o] + data.costPoolToProduct[o][p]) *\n",
    "            flowComponentToProductByPool[c][o][p]\n",
    "            for c in range(data.nbComponents)\n",
    "            for o in range(data.nbPools) for p in range(data.nbProducts)\n",
    "            if data.upperBoundFractionComponentToPool[c][o] > 0 and data.upperBoundPoolToProduct[o][p] > 0\n",
    "            and data.costComponentToPool[c][o] + data.costPoolToProduct[o][p] != 0)\n",
    "\n",
    "        # Gain of selling the final products\n",
    "        productsGain = model.sum((model.sum(flowComponentToProduct[c][p]\n",
    "                                            for c in range(data.nbComponents))\n",
    "                                  + model.sum(flowPoolToProduct[o][p]\n",
    "                                              for o in range(data.nbPools)))\n",
    "                                 * data.productPrices[p] for p in range(data.nbProducts))\n",
    "\n",
    "        # Cost of buying the components\n",
    "        componentsCost = model.sum(\n",
    "            (model.sum(flowComponentToProduct[c][p]\n",
    "                       for p in range(data.nbProducts))\n",
    "             + model.sum(flowComponentToProductByPool[c][o][p]\n",
    "                         for p in range(data.nbProducts)\n",
    "                         for o in range(data.nbPools)))\n",
    "            * data.componentPrices[c] for c in range(data.nbComponents))\n",
    "\n",
    "        profit = productsGain - componentsCost - (directFlowCost + undirectFlowCost)\n",
    "\n",
    "        # Maximize the total profit\n",
    "        model.maximize(profit)\n",
    "\n",
    "\n",
    "        # Solve with OptAgent\n",
    "        solution = optimizer.solve()\n",
    "        if not solution.feasible:\n",
    "            print(f\"No feasible solution found; Status = {solution.feasible}\")\n",
    "            return solution\n",
    "\n",
    "        #\n",
    "        # Write the solution\n",
    "        #\n",
    "        if output_file is not None:\n",
    "            with open(output_file, 'w', encoding=\"utf-8\") as f:\n",
    "                component_to_poduct = []\n",
    "                component_to_pool_fraction = []\n",
    "                pool_to_product = []\n",
    "\n",
    "                # Solution flows from the components to the products\n",
    "                for c in range(data.nbComponents):\n",
    "                    for p in range(data.nbProducts):\n",
    "                        component_to_poduct.append(\n",
    "                            {\"component\": data.componentNames[c],\n",
    "                             \"product\": data.productNames[p],\n",
    "                             \"flow\": flowComponentToProduct[c][p].value})\n",
    "\n",
    "                # Solution fraction of the inflow at pool o coming from the component c\n",
    "                for c in range(data.nbComponents):\n",
    "                    for o in range(data.nbPools):\n",
    "                        component_to_pool_fraction.append(\n",
    "                            {\"component\": data.componentNames[c],\n",
    "                             \"pool\": data.poolNames[o],\n",
    "                             \"flow\": fractionComponentToPool[c][o].value})\n",
    "\n",
    "                # Solution flows from the pools to the products\n",
    "                for o in range(data.nbPools):\n",
    "                    for p in range(data.nbProducts):\n",
    "                        pool_to_product.append(\n",
    "                            {\"pool\": data.poolNames[o],\n",
    "                             \"product\": data.productNames[p],\n",
    "                             \"flow\": flowPoolToProduct[o][p].value})\n",
    "\n",
    "                json.dump({\"objective\": profit.value, \"solution\":\n",
    "                    {\"component_to_pool_fraction\": component_to_pool_fraction,\n",
    "                     \"component_to_product\": component_to_poduct,\n",
    "                     \"pool_to_product\": pool_to_product}}, f)\n",
    "        print(f\"Profit = {profit.value}; Status = {solution.feasible}\")\n",
    "        return solution\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 本地运行\n",
    "\n",
    "在 notebook 所在目录执行以下 cell，即可调用一个混流实例。\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_pooling = main(INSTANCE_DIR / \"haverly1.json\", time_limit=1)\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python",
   "version": "3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}