Introduce 08-network

08_network/models/integer_lp/Test_model.ipynb

@@ -0,0 +1,244 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "id": "5533a1a5-19e4-4206-a7f5-6fbeb897cdda",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import jijmodeling as jm\n",
+    "\n",
+    "def build_ip_formulation() -> jm.Problem:\n",
+    "    \"\"\"Create integer programming formulation for the arc-based flow problem.\n",
+    "\n",
+    "    This model aligns with the ZPL (0-based) specification, formulating a\n",
+    "    multi-commodity flow problem with integer scaling and big-M constraints.\n",
+    "\n",
+    "    Sets:\n",
+    "        - N = {0..n−1}\n",
+    "        - A = {(i, j) | i ≠ j}\n",
+    "        - T = {(k, i, j) | i ≠ j and k ≠ j}\n",
+    "\n",
+    "    Parameters:\n",
+    "        - n (int): Number of nodes.\n",
+    "        - t (ndarray, shape (n, n)): Demand matrix, with zero diagonal.\n",
+    "        - M (int): Big-M constant for capacity constraints.\n",
+    "        - intscale (int): Integer scaling factor.\n",
+    "\n",
+    "    Variables:\n",
+    "        - x[i, j] ∈ {0, 1}: Binary arc selection variable.\n",
+    "        - f[k, i, j] ∈ ℤ, 0..intscale·M: Flow of commodity k on arc (i, j).\n",
+    "        - z ∈ ℤ, 0..intscale·M: Global upper bound on flow.\n",
+    "\n",
+    "    Objective:\n",
+    "        - Minimize z.\n",
+    "\n",
+    "    Constraints:\n",
+    "        - c1: ∀ i ∈ N: Σ_{j ≠ i} x[i, j] = 2 (out-degree = 2).\n",
+    "        - c2: ∀ j ∈ N: Σ_{i ≠ j} x[i, j] = 2 (in-degree = 2).\n",
+    "        - c11: ∀ (k, i), k ≠ i:\n",
+    "                  Σ_{j ≠ i} f[k, j, i] − Σ_{j ≠ i, j ≠ k} f[k, i, j]\n",
+    "                  = t[k, i]·intscale\n",
+    "        - c14: ∀ (k, i, j), i ≠ j, k ≠ j:\n",
+    "                  f[k, i, j] ≤ M·intscale·x[i, j]\n",
+    "        - c100: ∀ (i, j), i ≠ j:\n",
+    "                  Σ_{k ≠ j} f[k, i, j] ≤ z\n",
+    "\n",
+    "    Returns:\n",
+    "        jm.Problem: JijModeling problem instance with all variables,\n",
+    "        objective, and constraints defined.\n",
+    "    \"\"\"\n",
+    "    # ---- Placeholders ----\n",
+    "    n = jm.Placeholder(\"n\")\n",
+    "    t = jm.Placeholder(\"t\", ndim=2)   # (n,n)\n",
+    "    M = jm.Placeholder(\"M\")\n",
+    "    intscale = jm.Placeholder(\"intscale\")\n",
+    "\n",
+    "    # ---- Indices ----\n",
+    "    i = jm.Element(\"i\", belong_to=(0, n))\n",
+    "    j = jm.Element(\"j\", belong_to=(0, n))\n",
+    "    k = jm.Element(\"k\", belong_to=(0, n))\n",
+    "\n",
+    "    # ---- Vars ----\n",
+    "    x = jm.BinaryVar(\"x\", shape=(n, n), description=\"arc i->j selected\")\n",
+    "    f = jm.IntegerVar(\n",
+    "        \"f\", shape=(n, n, n),\n",
+    "        lower_bound=0, upper_bound=intscale * M,\n",
+    "        description=\"flow of commodity k on arc i->j\"\n",
+    "    )\n",
+    "    z = jm.IntegerVar(\"z\", lower_bound=0, upper_bound=intscale * M)\n",
+    "\n",
+    "    # ---- Problem & Objective ----\n",
+    "    problem = jm.Problem(\"ip_formulation\", sense=jm.ProblemSense.MINIMIZE)\n",
+    "    problem += z\n",
+    "\n",
+    "    # c1: ∀ i ∈ N :  Σ_{j ≠ i} x[i,j] = 2\n",
+    "    problem += jm.Constraint(\n",
+    "        \"c1_outdeg_eq_2\",\n",
+    "        jm.sum([(j, j != i)], x[i, j]) == 2,\n",
+    "        forall=[i]\n",
+    "    )\n",
+    "    \n",
+    "    # c2: ∀ j ∈ N :  Σ_{i ≠ j} x[i,j] = 2\n",
+    "    problem += jm.Constraint(\n",
+    "        \"c2_indeg_eq_2\",\n",
+    "        jm.sum([(i, i != j)], x[i, j]) == 2,\n",
+    "        forall=[j]\n",
+    "    )\n",
+    "\n",
+    "    # c11: flow balance\n",
+    "    problem += jm.Constraint(\n",
+    "        \"c11_flow_balance\",\n",
+    "        jm.sum([(j, j != i)], f[k, j, i])\n",
+    "        - jm.sum([(j, (j != i) & (j != k))], f[k, i, j])\n",
+    "        == t[k, i] * intscale,\n",
+    "        forall=[k, (i, k != i)]\n",
+    "    )\n",
+    "\n",
+    "    # c14: capacity bound\n",
+    "    problem += jm.Constraint(\n",
+    "        \"c14_capacity_by_x\",\n",
+    "        f[k, i, j] <= M * intscale * x[i, j],\n",
+    "        forall=[k, i, (j, (i != j) & (k != j))]\n",
+    "    )\n",
+    "\n",
+    "    # c100: z upper bound on flow\n",
+    "    problem += jm.Constraint(\n",
+    "        \"c100_z_upper_bounds_flow\",\n",
+    "        jm.sum([(k, k != j)], f[k, i, j]) <= z,\n",
+    "        forall=[i, (j, i != j)]\n",
+    "    )\n",
+    "\n",
+    "    return problem\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "id": "cc77aecb-7444-4ada-a553-3b6b9b8b23d5",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/latex": [
+       "$$\\begin{array}{cccc}\\text{Problem:} & \\text{ip\\_formulation} & & \\\\& & \\min \\quad \\displaystyle z & \\\\\\text{{s.t.}} & & & \\\\ & \\text{c100\\_z\\_upper\\_bounds\\_flow} & \\displaystyle \\sum_{\\substack{k = 0\\\\k \\neq j}}^{n - 1} f_{k, i, j} \\leq z & \\forall i \\in \\left\\{0,\\ldots,n - 1\\right\\} \\forall j \\in \\left\\{j \\in \\left\\{0,\\ldots,n - 1\\right\\} \\mid i \\neq j \\right\\} \\\\ & \\text{c11\\_flow\\_balance} & \\displaystyle \\sum_{\\substack{j = 0\\\\j \\neq i}}^{n - 1} f_{k, j, i} - \\sum_{\\substack{j = 0\\\\j \\neq i \\land j \\neq k}}^{n - 1} f_{k, i, j} = t_{k, i} \\cdot intscale & \\forall k \\in \\left\\{0,\\ldots,n - 1\\right\\} \\forall i \\in \\left\\{i \\in \\left\\{0,\\ldots,n - 1\\right\\} \\mid k \\neq i \\right\\} \\\\ & \\text{c14\\_capacity\\_by\\_x} & \\displaystyle f_{k, i, j} \\leq M \\cdot intscale \\cdot x_{i, j} & \\forall k \\in \\left\\{0,\\ldots,n - 1\\right\\} \\forall i \\in \\left\\{0,\\ldots,n - 1\\right\\} \\forall j \\in \\left\\{j \\in \\left\\{0,\\ldots,n - 1\\right\\} \\mid i \\neq j \\land k \\neq j \\right\\} \\\\ & \\text{c1\\_outdeg\\_eq\\_2} & \\displaystyle \\sum_{\\substack{j = 0\\\\j \\neq i}}^{n - 1} x_{i, j} = 2 & \\forall i \\in \\left\\{0,\\ldots,n - 1\\right\\} \\\\ & \\text{c2\\_indeg\\_eq\\_2} & \\displaystyle \\sum_{\\substack{i = 0\\\\i \\neq j}}^{n - 1} x_{i, j} = 2 & \\forall j \\in \\left\\{0,\\ldots,n - 1\\right\\} \\\\\\text{{where}} & & & \\\\& x & 2\\text{-dim binary variable}& \\text{arc i->j selected}\\\\& z & 0\\text{-dim integer variable}\\\\ & & \\text{lower bound: }0 & \\\\ & & \\text{upper bound: }intscale \\cdot M & \\\\& f & 3\\text{-dim integer variable}& \\text{flow of commodity k on arc i->j}\\\\ & & \\text{lower bound: }0 & \\\\ & & \\text{upper bound: }intscale \\cdot M & \\\\\\end{array}$$"
+      ],
+      "text/plain": [
+       "<jijmodeling.Problem at 0x115904600>"
+      ]
+     },
+     "execution_count": 2,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "problem = build_ip_formulation()\n",
+    "problem"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "id": "5931df18-ab64-4a0b-9828-1ff77eb60008",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# ---- 準備資料 ----\n",
+    "\n",
+    "def cut_matrix(t: list[list[int]], n: int) -> list[list[int]]:\n",
+    "    if not (1 <= n <= 24):\n",
+    "        raise ValueError(\"n must be between 1 and 24.\")\n",
+    "    return [row[:n] for row in t[:n]]\n",
+    "    \n",
+    "n_val = 6\n",
+    "your_t_0based = [\n",
+    "    [0, 24, 43, 23, 21, 41, 61, 21, 20, 0, 0, 0, 0, 0, 40, 0, 0, 0, 0, 0, 24, 19, 23, 0],\n",
+    "    [0, 0, 0, 0, 21, 18, 39, 23, 0, 0, 40, 19, 64, 19, 17, 64, 80, 0, 0, 22, 24, 18, 19, 0],\n",
+    "    [16, 0, 0, 20, 44, 0, 42, 22, 20, 0, 0, 23, 21, 0, 40, 21, 38, 97, 22, 17, 20, 37, 17, 20],\n",
+    "    [42, 20, 18, 0, 40, 57, 43, 42, 0, 0, 19, 17, 0, 39, 0, 0, 22, 41, 0, 0, 17, 42, 40, 43],\n",
+    "    [0, 0, 0, 40, 0, 83, 60, 0, 44, 0, 37, 60, 0, 40, 0, 0, 40, 19, 0, 39, 41, 17, 0, 0],\n",
+    "    [60, 18, 0, 37, 18, 0, 21, 41, 23, 39, 0, 63, 60, 39, 0, 19, 0, 16, 16, 0, 0, 40, 0, 16],\n",
+    "    [22, 0, 0, 0, 0, 0, 0, 61, 36, 80, 96, 19, 19, 41, 16, 0, 0, 0, 22, 0, 43, 0, 44, 22],\n",
+    "    [0, 0, 20, 19, 17, 20, 40, 0, 0, 60, 61, 0, 20, 62, 20, 0, 0, 38, 0, 0, 0, 0, 24, 22],\n",
+    "    [21, 22, 38, 0, 44, 20, 40, 39, 0, 36, 22, 21, 19, 39, 19, 0, 0, 21, 24, 16, 23, 21, 37, 0],\n",
+    "    [0, 24, 23, 39, 20, 0, 0, 41, 0, 0, 0, 22, 0, 0, 44, 42, 22, 42, 22, 19, 20, 58, 18, 0],\n",
+    "    [60, 57, 0, 0, 16, 0, 16, 37, 0, 0, 0, 44, 63, 0, 18, 0, 17, 18, 0, 0, 0, 100, 24, 23],\n",
+    "    [0, 44, 44, 0, 23, 17, 39, 21, 0, 17, 40, 0, 24, 78, 17, 24, 20, 18, 0, 24, 24, 0, 20, 0],\n",
+    "    [23, 16, 0, 0, 0, 23, 0, 0, 0, 0, 43, 58, 0, 0, 24, 60, 0, 0, 19, 0, 21, 0, 20, 0],\n",
+    "    [44, 20, 0, 19, 21, 0, 39, 19, 0, 0, 0, 39, 22, 0, 0, 64, 24, 22, 0, 39, 0, 43, 42, 16],\n",
+    "    [0, 60, 37, 18, 0, 0, 0, 20, 0, 41, 43, 16, 43, 24, 0, 0, 18, 18, 0, 44, 20, 0, 21, 37],\n",
+    "    [0, 0, 23, 39, 0, 24, 40, 0, 37, 0, 40, 20, 44, 43, 0, 0, 0, 0, 0, 16, 0, 59, 0, 0],\n",
+    "    [0, 42, 0, 0, 23, 24, 38, 19, 36, 0, 20, 60, 57, 0, 23, 40, 0, 0, 0, 16, 42, 0, 23, 0],\n",
+    "    [0, 41, 36, 43, 23, 41, 17, 0, 38, 0, 0, 21, 21, 17, 16, 16, 39, 0, 22, 0, 21, 23, 16, 23],\n",
+    "    [17, 0, 23, 23, 20, 0, 17, 58, 17, 0, 20, 0, 17, 24, 0, 0, 17, 42, 0, 58, 19, 22, 0, 24],\n",
+    "    [42, 0, 16, 0, 43, 0, 24, 36, 0, 16, 24, 41, 41, 0, 24, 0, 0, 0, 0, 0, 56, 38, 63, 19],\n",
+    "    [37, 23, 23, 0, 42, 16, 23, 76, 23, 0, 0, 24, 20, 41, 20, 24, 40, 23, 0, 0, 0, 39, 20, 0],\n",
+    "    [43, 20, 17, 17, 0, 20, 19, 0, 80, 0, 0, 0, 40, 0, 40, 16, 19, 0, 0, 0, 18, 0, 17, 0],\n",
+    "    [18, 20, 44, 40, 21, 18, 0, 20, 0, 0, 16, 24, 0, 0, 19, 18, 0, 17, 23, 0, 23, 44, 0, 42],\n",
+    "    [44, 0, 0, 0, 62, 0, 17, 41, 0, 0, 0, 63, 0, 37, 22, 0, 20, 0, 0, 19, 57, 18, 38, 0],\n",
+    "]\n",
+    "\n",
+    "instance_data = {\n",
+    "    \"n\": n_val,\n",
+    "    \"t\": cut_matrix(your_t_0based, n_val),\n",
+    "    \"M\": 1000,\n",
+    "    \"intscale\": 1000,\n",
+    "}\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "id": "40cc7ee7-ff07-47b3-9d5b-1c98a9677a50",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "objective=101000.0, feasible=True\n"
+     ]
+    }
+   ],
+   "source": [
+    "import ommx_pyscipopt_adapter as scip_ad\n",
+    "problem = build_ip_formulation()\n",
+    "ommx_instance = jm.Interpreter(instance_data).eval_problem(problem)\n",
+    "solution = scip_ad.OMMXPySCIPOptAdapter.solve(ommx_instance)\n",
+    "\n",
+    "print(f\"objective={solution.objective}, feasible={solution.feasible}\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "47de8dd9-3867-40b6-841b-1a53d8680907",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "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.10.6"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}

08_network/models/integer_lp/model.py

@@ -0,0 +1,104 @@
+import jijmodeling as jm
+
+def build_ip_formulation() -> jm.Problem:
+    """Create integer programming formulation for the arc-based flow problem.
+
+    This model aligns with the ZPL (0-based) specification, formulating a
+    multi-commodity flow problem with integer scaling and big-M constraints.
+
+    Sets:
+        - N = {0..n−1}
+        - A = {(i, j) | i ≠ j}
+        - T = {(k, i, j) | i ≠ j and k ≠ j}
+
+    Parameters:
+        - n (int): Number of nodes.
+        - t (ndarray, shape (n, n)): Demand matrix, with zero diagonal.
+        - M (int): Big-M constant for capacity constraints.
+        - intscale (int): Integer scaling factor.
+
+    Variables:
+        - x[i, j] ∈ {0, 1}: Binary arc selection variable.
+        - f[k, i, j] ∈ ℤ, 0..intscale·M: Flow of commodity k on arc (i, j).
+        - z ∈ ℤ, 0..intscale·M: Global upper bound on flow.
+
+    Objective:
+        - Minimize z.
+
+    Constraints:
+        - c1: ∀ i ∈ N: Σ_{j ≠ i} x[i, j] = 2 (out-degree = 2).
+        - c2: ∀ j ∈ N: Σ_{i ≠ j} x[i, j] = 2 (in-degree = 2).
+        - c11: ∀ (k, i), k ≠ i:
+                  Σ_{j ≠ i} f[k, j, i] − Σ_{j ≠ i, j ≠ k} f[k, i, j]
+                  = t[k, i]·intscale
+        - c14: ∀ (k, i, j), i ≠ j, k ≠ j:
+                  f[k, i, j] ≤ M·intscale·x[i, j]
+        - c100: ∀ (i, j), i ≠ j:
+                  Σ_{k ≠ j} f[k, i, j] ≤ z
+
+    Returns:
+        jm.Problem: JijModeling problem instance with all variables,
+        objective, and constraints defined.
+    """
+    # ---- Placeholders ----
+    n = jm.Placeholder("n")
+    t = jm.Placeholder("t", ndim=2)   # (n,n)
+    M = jm.Placeholder("M")
+    intscale = jm.Placeholder("intscale")
+
+    # ---- Indices ----
+    i = jm.Element("i", belong_to=(0, n))
+    j = jm.Element("j", belong_to=(0, n))
+    k = jm.Element("k", belong_to=(0, n))
+
+    # ---- Vars ----
+    x = jm.BinaryVar("x", shape=(n, n), description="arc i->j selected")
+    f = jm.IntegerVar(
+        "f", shape=(n, n, n),
+        lower_bound=0, upper_bound=intscale * M,
+        description="flow of commodity k on arc i->j"
+    )
+    z = jm.IntegerVar("z", lower_bound=0, upper_bound=intscale * M)
+
+    # ---- Problem & Objective ----
+    problem = jm.Problem("ip_formulation", sense=jm.ProblemSense.MINIMIZE)
+    problem += z
+
+    # c1: ∀ i ∈ N :  Σ_{j ≠ i} x[i,j] = 2
+    problem += jm.Constraint(
+        "c1_outdeg_eq_2",
+        jm.sum([(j, j != i)], x[i, j]) == 2,
+        forall=[i]
+    )
+
+    # c2: ∀ j ∈ N :  Σ_{i ≠ j} x[i,j] = 2
+    problem += jm.Constraint(
+        "c2_indeg_eq_2",
+        jm.sum([(i, i != j)], x[i, j]) == 2,
+        forall=[j]
+    )
+
+    # c11: flow balance
+    problem += jm.Constraint(
+        "c11_flow_balance",
+        jm.sum([(j, j != i)], f[k, j, i])
+        - jm.sum([(j, (j != i) & (j != k))], f[k, i, j])
+        == t[k, i] * intscale,
+        forall=[k, (i, k != i)]
+    )
+
+    # c14: capacity bound
+    problem += jm.Constraint(
+        "c14_capacity_by_x",
+        f[k, i, j] <= M * intscale * x[i, j],
+        forall=[k, i, (j, (i != j) & (k != j))]
+    )
+
+    # c100: z upper bound on flow
+    problem += jm.Constraint(
+        "c100_z_upper_bounds_flow",
+        jm.sum([(k, k != j)], f[k, i, j]) <= z,
+        forall=[i, (j, i != j)]
+    )
+
+    return problem