Exercise16_09.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "Exercise16.09",
      "provenance": [],
      "collapsed_sections": [],
      "authorship_tag": "ABX9TyMbKi4ncke0tcJFXdxu+Lia"
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "RhUgvN7XDOzQ"
      },
      "source": [
        "#Build model using data from file\n",
        "\n",
        "We use the `pandas` module to read in a large datafile that we then use to create a problem. The datafile contains the edges of a directed graph $G = (V,A)$ and we will use the data to create a **shortest-path problem** from the first to the last node (sorted by their index). The format of the data is: `i,j,c` where i and j are the end nodes of the directed $i \\to j$ arc, while $c$ is the cost of the arc."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "2hkhJl62CYe2",
        "outputId": "a84ae3f9-8f86-462b-f107-ce04e3f8338d"
      },
      "source": [
        "pip install xpress"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Collecting xpress\n",
            "\u001b[?25l  Downloading https://files.pythonhosted.org/packages/0b/c7/f82745387b9a0cc55d428aa92c7c680e52ed51a1ac987c9c045921599f0c/xpress-8.11.0-cp36-cp36m-manylinux1_x86_64.whl (43.8MB)\n",
            "\u001b[K     |████████████████████████████████| 43.8MB 88kB/s \n",
            "\u001b[?25hRequirement already satisfied: numpy>=1.8 in /usr/local/lib/python3.6/dist-packages (from xpress) (1.18.5)\n",
            "Installing collected packages: xpress\n",
            "Successfully installed xpress-8.11.0\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "DVzl4xZACbGC",
        "outputId": "913dcd52-f43d-4453-f11d-34d247ae870f"
      },
      "source": [
        "import xpress as xp\n",
        "import pandas as pd"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Using the Community license in this session. If you have a full Xpress license, first set the XPAUTH_PATH environment variable to the full path to your license file, xpauth.xpr, and then restart Python. If you want to use the FICO Community license and no longer want to see this message, set the XPAUTH_PATH environment variable to: /usr/local/lib/python3.6/dist-packages/xpress/license/community-xpauth.xpr\n",
            "NB: setting XPAUTH_PATH will also affect any other Xpress products installed on your system.\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "NunXLjMICf7F"
      },
      "source": [
        "#read graph data\n",
        "df = pd.read_csv('jupyter_2020-09-16_wednesday_graph_data.csv')"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "CyxE4a2JDGrb"
      },
      "source": [
        "tail = df['i']\n",
        "head = df['j']\n",
        "cost = df['c']"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "9g25lKvUFb16",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "3a2b388b-94f0-41c9-c16c-b85cc463e6a0"
      },
      "source": [
        "#V is created as a set\n",
        "V = set(tail) | set(head) #get set V of nodes\n",
        "\n",
        "#create A as a dictionary with cost as value\n",
        "A = {(tail[i],head[i]) : cost[i] for i, _ in enumerate(tail)}\n",
        "\n",
        "source = min(V)\n",
        "dest = max(V)\n",
        "\n",
        "#Let's print out V,A, source and dest\n",
        "\n",
        "print(\"Set of nodes:\", V)\n",
        "print(\"Arcs:\", A.keys())\n",
        "print(\"Find the shortest path from {0} to {1}\".format(source, dest))"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Set of nodes: {1, 2, 3, 4, 5}\n",
            "Arcs: dict_keys([(1, 2), (2, 1), (1, 3), (3, 1), (1, 4), (4, 1), (1, 5), (5, 1), (2, 3), (3, 2), (2, 4), (4, 2), (2, 5), (3, 4), (4, 3), (3, 5), (5, 3), (4, 5), (5, 4)])\n",
            "Find the shortest path from 1 to 5\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "orw2gqb7GNO0",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "339509db-4729-4907-b587-65e232ea938c"
      },
      "source": [
        "#Use the vars() functions to specify an index set for the variables\n",
        "x = xp.vars(list(A.keys()), vartype = xp.binary)\n",
        "\n",
        "#The right-hand side is also a dictionary. It specifies the\n",
        "#flow balance, i.e. the difference between outgoing and incoming\n",
        "# flow at a node. It must be 1 at the source, -1 at destination, and \n",
        "# 0 at each intermediate node.\n",
        "\n",
        "b = {i: 0 for i in V}\n",
        "b[source] = 1\n",
        "b[dest] = -1\n",
        "\n",
        "#conservation of flow: The total outgoing flow minus the total incomming flow \n",
        "#must equal the above RHS, for each node in V\n",
        "conservation = [xp.Sum(x[(i,j)] for j in V if (i,j) in A.keys()) - \\\n",
        "                xp.Sum(x[(j,i)] for j in V if (j,i) in A.keys()) == b[i] for i in V]\n",
        "\n",
        "#The objective function is the length of the path, i.e. the weighted sum of all arc variables\n",
        "obj = xp.Sum(A[(i,j)] * x[(i,j)] for (i,j) in A)\n",
        "\n",
        "p = xp.problem(x,conservation, obj)\n",
        "\n",
        "p.solve()\n",
        "\n",
        "for (i,j) in A:\n",
        "  if p.getSolution(x[(i,j)]) > 0.5:\n",
        "    print(i,j)"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "FICO Xpress v8.11.0, Community, solve started 9:01:42, Nov 20, 2020\n",
            "Heap usage: 341KB (peak 341KB, 559KB system)\n",
            "Minimizing MILP noname with these control settings:\n",
            "OUTPUTLOG = 1\n",
            "\n",
            "Original problem has:\n",
            "         5 rows           19 cols           38 elements        19 globals\n",
            "Presolved problem has:\n",
            "         5 rows           19 cols           38 elements        19 globals\n",
            "Presolve finished in 0 seconds\n",
            "Heap usage: 370KB (peak 371KB, 561KB system)\n",
            "\n",
            "Coefficient range                    original                 solved        \n",
            "  Coefficients   [min,max] : [ 1.00e+00,  1.00e+00] / [ 1.00e+00,  1.00e+00]\n",
            "  RHS and bounds [min,max] : [ 1.00e+00,  1.00e+00] / [ 1.00e+00,  1.00e+00]\n",
            "  Objective      [min,max] : [ 2.00e+00,  1.90e+01] / [ 2.00e+00,  1.90e+01]\n",
            "Autoscaling applied standard scaling\n",
            "\n",
            "Will try to keep branch and bound tree memory usage below 11.9GB\n",
            " *** Heuristic solution found:    19.000000      Time: 0 ***\n",
            " *** Heuristic solution found:     9.000000      Time: 0 ***\n",
            "Starting concurrent solve with dual\n",
            "\n",
            " Concurrent-Solve,   0s\n",
            "            Dual        \n",
            "    objective   dual inf\n",
            " D  9.0000000   .0000000\n",
            "-------- cutoff --------\n",
            "Concurrent statistics:\n",
            "      Dual: 2 simplex iterations, 0.00s\n",
            "Problem is cut off\n",
            " *** Search completed ***\n",
            "Uncrunching matrix\n",
            "Final MIP objective                   : 9.000000000000000e+00\n",
            "Final MIP bound                       : 9.000000000000000e+00\n",
            "  Solution time / primaldual integral :         0s/ 100.000000%\n",
            "  Number of solutions found / nodes   :         2 /         0\n",
            "  Max primal violation      (abs/rel) :       0.0 /       0.0\n",
            "  Max integer violation     (abs    ) :       0.0\n",
            "1 4\n",
            "4 5\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "RzcW1uNu7HjH"
      },
      "source": [
        "#Indicator Constraints:\n",
        "\n",
        "Indicator constraints are constraints that must hold if a given binary variable is equal to one, and they are not enforced otherwise. For instance, if a problem has two variables x and y and we want to constrain them to $x + 2y \\geq 4$ only if a third binary variable z is set to one, then this can be denoted as\n",
        "\n",
        "$z = 1 \\rightarrow x + 2y \\geq 4$\n",
        "\n",
        "Indicator contraints are common in Mixed Integer Linear Optimization. In the Xpress python interface, the above indicator constrain would be added to a problem as follows:\n",
        "p.addIndicator ( z == 1, x + 2*y >= 4)\n",
        "\n",
        "We can specify indicators with the `z ==0` condition too, but there are only two alternatives; constraints other that `z==0` or `z==1` as the first argument of an indicator will trigger an error.\n",
        "\n",
        "p.addIndicator (z == 0, x + 2*y > = 4)"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "9YEo1Mpx4rZb",
        "outputId": "e543f824-7df9-47ef-b2cc-f25f0a192831"
      },
      "source": [
        "x = xp.var()\n",
        "y = xp.var()\n",
        "\n",
        "p = xp.problem(x,y, xp.abs(x-y) + 2*xp.max(2*x + 4, 3*y + 2*x),\n",
        "               x+2*y >= 20,\n",
        "               x + 3*y <= 25)\n",
        "\n",
        "p.solve()\n",
        "\n",
        "print(p.getSolution(x), p.getSolution(y))"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "\n",
            "Original problem size\n",
            "   linear:    3 rows, 4 columns, 5 linear coefficients\n",
            "   nonlinear: 1 coefficients, 23 tokens\n",
            "Nonlinear presolve\n",
            "   compressing formula space (in use\\total : 24\\49)\n",
            "   removed 2 formulas and 24 formula tokens\n",
            "   checking equivalence to general constraints ABS/MIN/MAX and PWL reformulations\n",
            "   converted 1 ABS to optimizer MIP constructs\n",
            "   converted 1 MIN/MAX to optimizer MIP constructs\n",
            "   converted 4 formulas to linear constraints\n",
            "   converted objective transfer row to linear objective\n",
            "   problem is equivalent, reformulated as a MIP problem\n",
            "Presolved problem size\n",
            "   linear:    6 rows, 9 columns, 12 linear coefficients\n",
            "Problem is nonlinear presolved\n",
            "Minimizing problem using Xpress-Optimizer\n",
            "FICO Xpress v8.11.0, Community, solve started 9:46:15, Nov 20, 2020\n",
            "Heap usage: 688KB (peak 754KB, 942KB system)\n",
            "Minimizing MILP noname with these control settings:\n",
            "OUTPUTLOG = 1\n",
            "IFCHECKCONVEXITY = 0\n",
            "\n",
            "Original problem has:\n",
            "         6 rows            9 cols           12 elements         0 globals\n",
            "         2 gencons\n",
            "Presolved problem has:\n",
            "         0 rows            0 cols            0 elements         0 globals\n",
            "Presolve finished in 0 seconds\n",
            "Heap usage: 692KB (peak 754KB, 961KB system)\n",
            "Will try to keep branch and bound tree memory usage below 11.9GB\n",
            "Starting concurrent solve with dual\n",
            "\n",
            " Concurrent-Solve,   0s\n",
            "            Dual        \n",
            "    objective   dual inf\n",
            " D  75.000000   .0000000\n",
            "------- optimal --------\n",
            "Concurrent statistics:\n",
            "      Dual: 0 simplex iterations, 0.00s\n",
            "Optimal solution found\n",
            " \n",
            "   Its         Obj Value      S   Ninf  Nneg   Sum Dual Inf  Time\n",
            "     0         75.000000      D      0     0        .000000     0\n",
            "Dual solved problem\n",
            "  0 simplex iterations in 0s\n",
            "\n",
            "Final objective                       : 7.500000000000000e+01\n",
            "  Max primal violation      (abs/rel) :       0.0 /       0.0\n",
            "  Max dual violation        (abs/rel) :       0.0 /       0.0\n",
            "  Max complementarity viol. (abs/rel) :       0.0 /       0.0\n",
            "\n",
            "Starting root cutting & heuristics\n",
            " \n",
            " Its Type    BestSoln    BestBound   Sols    Add    Del     Gap     GInf   Time\n",
            "*           75.000000    75.000000      1                  0.00%       0      0\n",
            " *** Search completed ***\n",
            "Uncrunching matrix\n",
            "Final MIP objective                   : 7.500000000000000e+01\n",
            "Final MIP bound                       : 7.500000000000000e+01\n",
            "  Solution time / primaldual integral :         0s/ 97.564351%\n",
            "  Number of solutions found / nodes   :         1 /         1\n",
            "  Max primal violation      (abs/rel) :       0.0 /       0.0\n",
            "  Max integer violation     (abs    ) :       0.0\n",
            "Problem is nonlinear postsolved\n",
            "Heap usage: 837KB (peak 3392KB, 627KB system)\n",
            "Final NLP objective (feas. solution)    : 7.500000000000000e+01\n",
            "  Max validation error      (abs/rel) :      .000 /      .000\n",
            "  Observed primal integral :  100.000%\n",
            "  Total time :                  .014s \n",
            "*** Search completed ***\n",
            "10.0 5.0\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "B21Edr0gHOCB"
      },
      "source": [
        "#Manual MIP reformulation of $|x|$, max/min, and ***AND/OR***\n",
        "\n",
        "\n",
        "\n",
        "In the above problem, we defined a quantity $|x-y|$ with the `xpress.abs` operator.\n",
        "\n",
        "Below is how we would have to model the constraint $z = |x-y|$ if we did not have the `xp.abs` function. We would have to model $z = |x-y|$ by introducing a binary variable $w$ and writing the following constraints, two of which are indicator constraints:\n",
        "\n",
        "$$\n",
        "\\begin{array}{lll}\n",
        "z \\ge x-y\\\\\n",
        "z \\ge y-x\\\\\n",
        "w = 1 \\Rightarrow z \\le x-y\\\\\n",
        "w = 0 \\Rightarrow z \\le y-x\\\\\n",
        "w \\in \\{0,1\\}\n",
        "\\end{array}\n",
        "$$\n",
        "\n",
        "If, however, we only had the *inequality* constraint $z \\ge |x-y|$, then reformulation would be a lot easier: the two constraints\n",
        "\n",
        "$$\n",
        "z \\ge x - y; \\qquad z \\ge y-x\n",
        "$$\n",
        "\n",
        "would suffice."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "FSrUgucoIPpM"
      },
      "source": [
        "And/Or operator do not require extra variable: for an operation $z = x \\wedge y$, with $x$ and $y$ binary variables, one can write\n",
        "\n",
        "$$\n",
        "\\begin{array}{l}\n",
        "z \\le x\\\\\n",
        "z \\le y\\\\\n",
        "z \\ge x + y - 1\n",
        "\\end{array}\n",
        "$$\n",
        "\n",
        "and similar for an operation $z = x \\vee y$:\n",
        "\n",
        "$$\n",
        "\\begin{array}{l}\n",
        "z \\ge x\\\\\n",
        "z \\ge y\\\\\n",
        "z \\le x + y\n",
        "\\end{array}\n",
        "$$\n",
        "\n",
        "Note that $z\\le x$ above works like an implication: it is equivalent to $z = 1 \\Rightarrow x = 1$, and therefore also an implication in the opposite direction (from $x$ to $z$): $x = 0 \\Rightarrow z = 0$."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "xak8zOAaC6Yk",
        "outputId": "42249038-5ee6-48f2-ff83-02b9e188d713"
      },
      "source": [
        "x = xp.var(vartype = xp.binary)\n",
        "y = xp.var(vartype = xp.binary)\n",
        "z = xp.vars(10, vartype = xp.binary)\n",
        "\n",
        "p = xp.problem(x,y,z,\n",
        "               (x&y) + x | y | xp.Or(z[0], z[1], z[4]) == 1,\n",
        "               xp.Sum(z))\n",
        "\n",
        "p.solve()"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "\n",
            "Original problem size\n",
            "   linear:    1 rows, 13 columns, 0 linear coefficients\n",
            "   nonlinear: 1 coefficients, 16 tokens\n",
            "Nonlinear presolve\n",
            "   checking equivalence to general constraints ABS/MIN/MAX and PWL reformulations\n",
            "   converted 4 MIN/MAX to optimizer MIP constructs\n",
            "   converted 2 formulas to linear constraints\n",
            "   problem is equivalent, reformulated as a MIP problem\n",
            "Presolved problem size\n",
            "   linear:    2 rows, 18 columns, 4 linear coefficients\n",
            "Problem is nonlinear presolved\n",
            "Minimizing problem using Xpress-Optimizer\n",
            "FICO Xpress v8.11.0, Community, solve started 12:48:06, Nov 20, 2020\n",
            "Heap usage: 686KB (peak 752KB, 947KB system)\n",
            "Minimizing MILP noname with these control settings:\n",
            "OUTPUTLOG = 1\n",
            "IFCHECKCONVEXITY = 0\n",
            "\n",
            "Original problem has:\n",
            "         2 rows           18 cols            4 elements        12 globals\n",
            "         4 gencons\n",
            "Presolved problem has:\n",
            "         5 rows            5 cols           13 elements         4 globals\n",
            "LP relaxation tightened\n",
            "Presolve finished in 0 seconds\n",
            "Heap usage: 718KB (peak 752KB, 966KB system)\n",
            "\n",
            "Coefficient range                    original                 solved        \n",
            "  Coefficients   [min,max] : [ 1.00e+00,  1.00e+00] / [ 1.00e+00,  1.00e+00]\n",
            "  RHS and bounds [min,max] : [ 1.00e+00,  1.00e+00] / [ 1.00e+00,  1.00e+00]\n",
            "  Objective      [min,max] : [ 1.00e+00,  1.00e+00] / [ 1.00e+00,  1.00e+00]\n",
            "Autoscaling applied standard scaling\n",
            "\n",
            "Will try to keep branch and bound tree memory usage below 11.7GB\n",
            "Starting concurrent solve with dual\n",
            "\n",
            " Concurrent-Solve,   0s\n",
            "            Dual        \n",
            "    objective   dual inf\n",
            " D   .0000000   .0000000\n",
            "------- optimal --------\n",
            "Concurrent statistics:\n",
            "      Dual: 3 simplex iterations, 0.00s\n",
            "Optimal solution found\n",
            " \n",
            "   Its         Obj Value      S   Ninf  Nneg   Sum Dual Inf  Time\n",
            "     3           .000000      D      0     0        .000000     0\n",
            "Dual solved problem\n",
            "  3 simplex iterations in 0s\n",
            "\n",
            "Final objective                       : 0.000000000000000e-01\n",
            "  Max primal violation      (abs/rel) :       0.0 /       0.0\n",
            "  Max dual violation        (abs/rel) :       0.0 /       0.0\n",
            "  Max complementarity viol. (abs/rel) :       0.0 /       0.0\n",
            "\n",
            "Starting root cutting & heuristics\n",
            " \n",
            " Its Type    BestSoln    BestBound   Sols    Add    Del     Gap     GInf   Time\n",
            "a             .000000      .000000      1               0.0e+00        0      0\n",
            " *** Search completed ***\n",
            "Uncrunching matrix\n",
            "Final MIP objective                   : 0.000000000000000e-01\n",
            "Final MIP bound                       : 0.000000000000000e-01\n",
            "  Solution time / primaldual integral :         0s/ 98.066958%\n",
            "  Number of solutions found / nodes   :         1 /         1\n",
            "  Max primal violation      (abs/rel) :       0.0 /       0.0\n",
            "  Max integer violation     (abs    ) :       0.0\n",
            "Problem is nonlinear postsolved\n",
            "Heap usage: 840KB (peak 3443KB, 631KB system)\n",
            "Final MINLP objective (feas. solution)  : 0.000000000000000e+00\n",
            "  Max validation error      (abs/rel) :      .000 /      .000\n",
            "  Observed primal integral :  100.000%\n",
            "  Total time :                  .016s \n",
            "*** Search completed ***\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "_GC--mqJs27L"
      },
      "source": [
        "**Exercise**: Formulate a MIP to solve a SATisfaction problem defined by the clauses on three variables as defined below, where $x_i$ means a variable is included as-is and $\\bar x_i$ means its negation:\n",
        "\n",
        "$$\n",
        "\\begin{array}{lllll}\n",
        "  \\bar x_1 &\\vee&      x_2 &    &          \\\\\n",
        "           &    & \\bar x_2 &\\vee&      x_3 \\\\\n",
        "  \\bar x_1 &    &          &\\vee&      x_3 \\\\\n",
        "           &    & \\bar x_2 &\\vee& \\bar x_3 \\\\\n",
        "       x_1 &    &          &\\vee&      x_3 \\\\\n",
        "       x_1 &\\vee& \\bar x_2 &    &          \\\\\n",
        "\\end{array}\n",
        "$$"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "mmuwazL7s0pw",
        "outputId": "75d85ea0-182d-4287-b474-9446d5cef7f9"
      },
      "source": [
        "# Identify each clause with -1 if variable is negated, 1 if non-negated, 0 if not present.\n",
        "\n",
        "id_clauses = [\n",
        "              [-1,1,0],\n",
        "              [0,-1,1],\n",
        "              [-1,0,1],\n",
        "              [0,-1,-1],\n",
        "              [1,0,1],\n",
        "              [1,-1,0]]\n",
        "I = [0,1,2]\n",
        "\n",
        "x = [xp.var(vartype = xp.binary) for _ in I]\n",
        "\n",
        "clauses = [\n",
        "           [(x[j] if id[j] == 1 else (1-x[j]) if id[j] == -1 else 0)\n",
        "           for j in I]\n",
        "           for id in id_clauses\n",
        "]\n",
        "\n",
        "SATcon = [xp.Or(*clause) == 1 for clause in clauses]\n",
        "p = xp.problem(x,SATcon) #note: no objective function is necessary\n",
        "\n",
        "p.solve()\n",
        "print(p.getSolution())"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "\n",
            "Original problem size\n",
            "   linear:    6 rows, 4 columns, 0 linear coefficients\n",
            "   nonlinear: 6 coefficients, 48 tokens\n",
            "Nonlinear presolve\n",
            "   checking equivalence to general constraints ABS/MIN/MAX and PWL reformulations\n",
            "   converted 6 MIN/MAX to optimizer MIP constructs\n",
            "   converted 12 formulas to linear constraints\n",
            "   problem is equivalent, reformulated as a MIP problem\n",
            "Presolved problem size\n",
            "   linear:    12 rows, 16 columns, 18 linear coefficients\n",
            "Problem is nonlinear presolved\n",
            "Minimizing problem using Xpress-Optimizer\n",
            "FICO Xpress v8.11.0, Community, solve started 14:11:06, Nov 20, 2020\n",
            "Heap usage: 687KB (peak 753KB, 951KB system)\n",
            "Minimizing MILP noname with these control settings:\n",
            "OUTPUTLOG = 1\n",
            "IFCHECKCONVEXITY = 0\n",
            "\n",
            "Original problem has:\n",
            "        12 rows           16 cols           18 elements         3 globals\n",
            "         6 gencons\n",
            "Presolved problem has:\n",
            "         0 rows            0 cols            0 elements         0 globals\n",
            "LP relaxation tightened\n",
            "Presolve finished in 0 seconds\n",
            "Heap usage: 698KB (peak 753KB, 970KB system)\n",
            "Will try to keep branch and bound tree memory usage below 11.7GB\n",
            "Starting concurrent solve with dual\n",
            "\n",
            " Concurrent-Solve,   0s\n",
            "            Dual        \n",
            "    objective   dual inf\n",
            " D   .0000000   .0000000\n",
            "------- optimal --------\n",
            "Concurrent statistics:\n",
            "      Dual: 0 simplex iterations, 0.00s\n",
            "Optimal solution found\n",
            " \n",
            "   Its         Obj Value      S   Ninf  Nneg   Sum Dual Inf  Time\n",
            "     0           .000000      D      0     0        .000000     0\n",
            "Dual solved problem\n",
            "  0 simplex iterations in 0s\n",
            "\n",
            "Final objective                       : 0.000000000000000e-01\n",
            "  Max primal violation      (abs/rel) :       0.0 /       0.0\n",
            "  Max dual violation        (abs/rel) :       0.0 /       0.0\n",
            "  Max complementarity viol. (abs/rel) :       0.0 /       0.0\n",
            "\n",
            "Starting root cutting & heuristics\n",
            " \n",
            " Its Type    BestSoln    BestBound   Sols    Add    Del     Gap     GInf   Time\n",
            "*             .000000      .000000      1               0.0e+00        0      0\n",
            " *** Search completed ***\n",
            "Uncrunching matrix\n",
            "Final MIP objective                   : 0.000000000000000e-01\n",
            "Final MIP bound                       : 0.000000000000000e-01\n",
            "  Solution time / primaldual integral :         0s/ 97.994749%\n",
            "  Number of solutions found / nodes   :         1 /         1\n",
            "  Max primal violation      (abs/rel) :       0.0 /       0.0\n",
            "  Max integer violation     (abs    ) :       0.0\n",
            "Problem is nonlinear postsolved\n",
            "Heap usage: 837KB (peak 3398KB, 636KB system)\n",
            "Final MINLP objective (feas. solution)  : 0.000000000000000e+00\n",
            "  Max validation error      (abs/rel) :      .000 /      .000\n",
            "  Observed primal integral :  100.000%\n",
            "  Total time :              9.99E-03s \n",
            "*** Search completed ***\n",
            "[0.0, 0.0, 1.0, 1.0]\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Ku9jExIf_xIe"
      },
      "source": [
        "#Simple cutting plane algorithm\n",
        "We'll implement a very simple cutting plane algorithm. Consider the problem\n",
        "\n",
        "$$\n",
        "\\begin{array}{ll}\n",
        "\\min &\\sum_{i=1}^n\\sum_{j=1}^n q_{ij} x_i x_j + \\sum_i c_i x_i\\\\\n",
        "\\textrm{s.t.} & x_i\\in \\{0,1\\} & \\forall i=1,2,\\ldots, n.\n",
        "\\end{array}\n",
        "$$\n",
        "\n",
        "This problem is a **Binary Box-constrained Quadratic Programming problem**: it has a quadratic objective function and all variables are negative; furthermore, the only constraints are the implicit bounds on binary variables, i.e. $\\{0,1\\}$. The problem can be written in matricial form as follows:\n",
        "\n",
        "$$\n",
        "\\begin{array}{ll}\n",
        "\\min & x^T Q x + c^T x\\\\\n",
        "\\textrm{s.t.} & x \\in \\{0,1\\}^n\n",
        "\\end{array}\n",
        "$$\n",
        "\n",
        "where $Q$ is a symmetric matrix. In general, the above problem is much harder if the term $x^T Q x$ is nonconvex, i.e., if $Q$ is **not** positive semidefinite. Let's try an example."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "0RcVd0Uk_cXc",
        "outputId": "42913075-6f02-4d21-9480-4032b2651e1b"
      },
      "source": [
        "# 2BCDC\n",
        "import numpy as np\n",
        "np.random.seed(1234567)\n",
        "\n",
        "n = 30\n",
        "Q = np.random.random((n,n)) - .5\n",
        "x = xp.vars(n, vartype = xp.binary) \n",
        "c = np.random.random(n) - .5\n",
        "\n",
        "p = xp.problem(x,xp.Dot(x,Q,x) + xp.Dot(c,x))\n",
        "\n",
        "p.solve()"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "\n",
            "FICO Xpress v8.11.0, Community, solve started 14:34:23, Nov 20, 2020\n",
            "Heap usage: 356KB (peak 356KB, 636KB system)\n",
            "Minimizing MIQP noname with these control settings:\n",
            "OUTPUTLOG = 1\n",
            "\n",
            "Original problem has:\n",
            "         0 rows           30 cols            0 elements        30 globals\n",
            "       900 qobjelem\n",
            "Presolved problem has:\n",
            "       657 rows          465 cols         1527 elements        30 globals\n",
            "Presolve finished in 0 seconds\n",
            "Heap usage: 655KB (peak 899KB, 638KB system)\n",
            "\n",
            "Coefficient range                    original                 solved        \n",
            "  Coefficients   [min,max] : [      0.0,       0.0] / [ 1.00e+00,  1.00e+00]\n",
            "  RHS and bounds [min,max] : [ 1.00e+00,  1.00e+00] / [ 1.00e+00,  1.00e+00]\n",
            "  Objective      [min,max] : [ 4.33e-03,  4.81e-01] / [ 2.13e-04,  9.83e-01]\n",
            "  Quadratic      [min,max] : [ 2.13e-04,  9.83e-01] / [      0.0,       0.0]\n",
            "Autoscaling applied standard scaling\n",
            "\n",
            "Will try to keep branch and bound tree memory usage below 11.7GB\n",
            "Starting concurrent solve with dual\n",
            "\n",
            " Concurrent-Solve,   0s\n",
            "            Dual        \n",
            "    objective   dual inf\n",
            " D -34.503919   .0000000\n",
            "------- optimal --------\n",
            "Concurrent statistics:\n",
            "      Dual: 365 simplex iterations, 0.01s\n",
            "Optimal solution found\n",
            " \n",
            "   Its         Obj Value      S   Ninf  Nneg   Sum Dual Inf  Time\n",
            "   365        -34.503919      D      0     0        .000000     0\n",
            "Dual solved problem\n",
            "  365 simplex iterations in 0s\n",
            "\n",
            "Final objective                       : -3.450391856833885e+01\n",
            "  Max primal violation      (abs/rel) :       0.0 /       0.0\n",
            "  Max dual violation        (abs/rel) :       0.0 /       0.0\n",
            "  Max complementarity viol. (abs/rel) :       0.0 /       0.0\n",
            "\n",
            "Starting root cutting & heuristics\n",
            " \n",
            " Its Type    BestSoln    BestBound   Sols    Add    Del     Gap     GInf   Time\n",
            "c             .661406   -34.503919      1               35.1653        0      0\n",
            "d           -2.989883   -34.503919      2                 91.33%       0      0\n",
            "S           -3.366240   -34.503919      3                 90.24%       0      0\n",
            "   1  K     -3.366240   -24.686089      3    344      0   86.36%     252      0\n",
            "d           -5.347880   -21.336445      4                 74.94%       0      0\n",
            "   2  K     -5.347880   -21.336445      4    295    229   74.94%     331      0\n",
            "   3  K     -5.347880   -19.917226      4    274    236   73.15%     276      0\n",
            "d          -15.135307   -19.038121      5                 20.50%       0      0\n",
            "   4  K    -15.135307   -19.038121      5    257    251   20.50%     315      0\n",
            "   5  K    -15.135307   -17.791810      5    304    250   14.93%     414      0\n",
            "d          -17.598893   -17.598893      6                 -0.00%       0      0\n",
            " *** Search completed ***\n",
            "Uncrunching matrix\n",
            "Final MIP objective                   : -1.759889266754942e+01\n",
            "Final MIP bound                       : -1.759889266754932e+01\n",
            "  Solution time / primaldual integral :         0s/ 63.979947%\n",
            "  Number of solutions found / nodes   :         6 /         1\n",
            "  Max primal violation      (abs/rel) :       0.0 /       0.0\n",
            "  Max integer violation     (abs    ) :       0.0\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "032gGNbRFoex"
      },
      "source": [
        "The Xpress Optimizer can solve such problems after a *reformulation* that does the following:\n",
        "\n",
        "* Creates a (binary) variable $y_{ij}$ for each product $x_i x_j$;\n",
        "* For each variable $y_{ij}$, it adds the following constraints linking $y$ with $x$ variables:\n",
        "\n",
        "$$\n",
        "\\begin{array}{lll}\n",
        "y_{ij} &\\le& x_i \\\\\n",
        "y_{ij} &\\le& x_j \\\\\n",
        "y_{ij} &\\ge& x_i + x_j - 1\n",
        "\\end{array}\n",
        "$$\n",
        "\n",
        "After this reformulation, the problem becomes as follows:\n",
        "\n",
        "$$\n",
        "\\begin{array}{ll}\n",
        "\\min & \\sum_{i=1}^n \\sum_{j=1}^n q_{ij} y_{ij} + \\sum_{i=1}^n c_i x_i\\\\\n",
        "\\textrm{s.t.} &y_{ij} \\le x_i  &\\forall i,j\\\\\n",
        "&y_{ij} \\le x_j  &\\forall i,j\\\\\n",
        "&y_{ij} \\ge x_i + x_j - 1 &\\forall i,j\\\\\n",
        "&x_i \\in \\{0,1\\} \\forall i=1,2,\\ldots,n\\\\\n",
        "&y_{ij} \\in \\{0,1\\} \\forall i,j=1,2,\\ldots,n\\\\\n",
        "\\end{array}\n",
        "$$\n",
        "\n",
        "We can do this reformulation explicitly by introducing the variables $y_{ij}$. Below we use NumPy vectors of variables as they come in handy to write constraints."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "guz8PeWQE69U",
        "outputId": "0b340c8a-4ced-4517-d424-423b5e8c25fd"
      },
      "source": [
        "# 2BCDC\n",
        "\n",
        "import xpress as xp\n",
        "import numpy as np\n",
        "\n",
        "np.random.seed(1234567)\n",
        "\n",
        "n = 30\n",
        "Q = np.random.random((n,n)) - 0.5\n",
        "c = np.random.random(n) - 0.5\n",
        "\n",
        "x = xp.vars(n, vartype=xp.binary)\n",
        "Y = xp.vars(n, n, vartype=xp.binary)\n",
        "\n",
        "p = xp.problem(x, Y, xp.Sum(Q*Y) + xp.Dot(c,x),\n",
        "               Y == Y.transpose(),\n",
        "               Y <= x,\n",
        "               Y.transpose() <= x,\n",
        "               [Y[i,j] >= x[i] + x[j] - 1 for i in range(n) for j in range(n)])\n",
        "p.solve()\n",
        "\n"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "\n",
            "FICO Xpress v8.11.0, Community, solve started 14:38:58, Nov 20, 2020\n",
            "Heap usage: 1689KB (peak 1689KB, 1575KB system)\n",
            "Minimizing MILP noname with these control settings:\n",
            "OUTPUTLOG = 1\n",
            "\n",
            "Original problem has:\n",
            "      3600 rows          930 cols         8010 elements       930 globals\n",
            "Presolved problem has:\n",
            "      1305 rows          465 cols         3045 elements       465 globals\n",
            "LP relaxation tightened\n",
            "Presolve finished in 0 seconds\n",
            "Heap usage: 2034KB (peak 3519KB, 1577KB system)\n",
            "\n",
            "Coefficient range                    original                 solved        \n",
            "  Coefficients   [min,max] : [ 1.00e+00,  2.00e+00] / [ 1.00e+00,  1.00e+00]\n",
            "  RHS and bounds [min,max] : [ 1.00e+00,  1.00e+00] / [ 1.00e+00,  1.00e+00]\n",
            "  Objective      [min,max] : [ 7.60e-04,  5.00e-01] / [ 2.13e-04,  9.83e-01]\n",
            "Autoscaling applied standard scaling\n",
            "\n",
            "Will try to keep branch and bound tree memory usage below 11.7GB\n",
            " *** Heuristic solution found:      .000000      Time: 0 ***\n",
            " *** Heuristic solution found:     -.107342      Time: 0 ***\n",
            "Starting concurrent solve with dual\n",
            "\n",
            " Concurrent-Solve,   0s\n",
            "            Dual        \n",
            "    objective   dual inf\n",
            " D -34.503919   .0000000\n",
            "------- optimal --------\n",
            "Concurrent statistics:\n",
            "      Dual: 373 simplex iterations, 0.01s\n",
            "Optimal solution found\n",
            " \n",
            "   Its         Obj Value      S   Ninf  Nneg   Sum Dual Inf  Time\n",
            "   373        -34.503919      D      0     0        .000000     0\n",
            "Dual solved problem\n",
            "  373 simplex iterations in 0s\n",
            "\n",
            "Final objective                       : -3.450391856833882e+01\n",
            "  Max primal violation      (abs/rel) :       0.0 /       0.0\n",
            "  Max dual violation        (abs/rel) :       0.0 /       0.0\n",
            "  Max complementarity viol. (abs/rel) :       0.0 /       0.0\n",
            "\n",
            "Starting root cutting & heuristics\n",
            " \n",
            " Its Type    BestSoln    BestBound   Sols    Add    Del     Gap     GInf   Time\n",
            "a            -.923150   -34.503919      3                 97.32%       0      0\n",
            "g            -.994858   -34.503919      4                 97.12%       0      0\n",
            "q           -6.950343   -34.503919      5                 79.86%       0      0\n",
            "   1  K     -6.950343   -34.285921      5      1      0   79.73%     252      0\n",
            "   2  K     -6.950343   -34.246365      5      2      0   79.70%     253      0\n",
            "q           -7.871428   -34.103575      6                 76.92%       0      0\n",
            "   3  K     -7.871428   -34.103575      6      2      1   76.92%     252      0\n",
            "   4  K     -7.871428   -33.974131      6      3      1   76.83%     251      0\n",
            "   5  K     -7.871428   -33.939319      6      5      1   76.81%     254      0\n",
            "   6  K     -7.871428   -33.939319      6      0      4   76.81%     255      0\n",
            "Heuristic search started\n",
            "Heuristic search stopped\n",
            " \n",
            "Cuts in the matrix         : 6\n",
            "Cut elements in the matrix : 93\n",
            "\n",
            "Starting tree search.\n",
            "Deterministic mode with up to 2 running threads and up to 4 tasks.\n",
            "Heap usage: 6271KB (peak 6679KB, 1637KB system)\n",
            " \n",
            "    Node     BestSoln    BestBound   Sols Active  Depth     Gap     GInf   Time\n",
            "       1    -7.871428   -32.163133      6      2      1   75.53%     255      0\n",
            "       2    -7.871428   -31.605384      6      3      3   75.09%     275      0\n",
            "       3    -7.871428   -31.605384      6      4      4   75.09%     258      0\n",
            "       4    -7.871428   -31.605384      6      4      5   75.09%     242      0\n",
            "       5    -7.871428   -30.882686      6      5      3   74.51%     247      0\n",
            "       6    -7.871428   -29.864686      6      6      4   73.64%     219      0\n",
            "       7    -7.871428   -29.864686      6      7      5   73.64%     190      0\n",
            "       8    -7.871428   -29.864686      6      8      4   73.64%     258      0\n",
            "       9    -7.871428   -28.937935      6      9      4   72.80%     258      0\n",
            "      10    -7.871428   -27.704491      6     10      5   71.59%     266      0\n",
            "      21    -7.871428   -27.605685      6     14     12   71.49%     259      1\n",
            "q     21   -15.152430   -27.005839      7     21     12   43.89%       0      1\n",
            "      32   -15.152430   -27.005839      7     20     11   43.89%     164      1\n",
            "q     40   -15.225044   -27.005839      8     26     11   43.62%       0      1\n",
            "      42   -15.225044   -27.005839      8     26     11   43.62%     229      1\n",
            "      53   -15.225044   -26.669538      8     33      9   42.91%     273      1\n",
            "q     53   -16.979294   -26.669538      9     33      9   36.33%       0      1\n",
            "      63   -16.979294   -26.669538      9     34      7   36.33%     255      1\n",
            "      74   -16.979294   -26.669538      9     33     10   36.33%     205      1\n",
            "Elapsed time (sec): 1, estimated tree completion: 0.03809\n",
            "Heap usage: 28MB (peak 28MB, 2639KB system)\n",
            "B&B tree size: 150k total\n",
            " \n",
            "    Node     BestSoln    BestBound   Sols Active  Depth     Gap     GInf   Time\n",
            "      84   -16.979294   -26.646519      9     41      9   36.28%     219      1\n",
            "      96   -16.979294   -26.646519      9     43     10   36.28%     211      1\n",
            "     107   -16.979294   -26.646519      9     45      9   36.28%     245      1\n",
            "q    120   -17.005543   -25.997879     10     48     11   34.59%       0      1\n",
            "q    121   -17.449548   -25.997879     11     48     10   32.88%       0      1\n",
            "     208   -17.449548   -24.832347     11     53      8   29.73%     256      1\n",
            "     314   -17.449548   -23.770950     11     48     10   26.59%     199      2\n",
            "     415   -17.449548   -21.525903     11     35      8   18.94%     229      2\n",
            "q    452   -17.598893   -20.681472     12     26      9   14.91%       0      2\n",
            " *** Search completed ***\n",
            "Uncrunching matrix\n",
            "Final MIP objective                   : -1.759889266756088e+01\n",
            "Final MIP bound                       : -1.759892717147945e+01\n",
            "  Solution time / primaldual integral :         2s/ 51.390333%\n",
            "  Number of solutions found / nodes   :        12 /       491\n",
            "  Max primal violation      (abs/rel) : 2.659e-12 / 2.659e-12\n",
            "  Max integer violation     (abs    ) : 2.659e-12\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "kqR5L-OCGA_J"
      },
      "source": [
        ""
      ],
      "execution_count": null,
      "outputs": []
    }
  ]
}