Exercise15_09.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "Exercise15.09.ipynb",
      "provenance": [],
      "collapsed_sections": [],
      "authorship_tag": "ABX9TyP8uQC2w9gCcVpknBBC8h8+"
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "KTcac03HG0Rp"
      },
      "source": [
        "#Exercise Session 1: Overview of the Xpress python Interface\n",
        "\n",
        "###Creating a simple LP\n",
        "\n",
        ">$min\\ x + y$ \\\\\n",
        ">$ s.t.\\ 2x + 3y \\geq 6$ \\\\\n",
        ">$ 4x + 2y \\geq 7$ \\\\\n",
        ">$ x,y \\geq 0$ \\\\\n",
        "###Also:\n",
        "***Knapsack Problem***\n",
        "\n",
        "***NumPy operators and arrays of variables constraints*** \\\\\n",
        "\n",
        "***Production planning with quadratic cost***\n",
        "\n",
        "\n",
        "\n",
        "Exercise is a part of CO@Work Summer School 2020 by ZIB."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "yWKXpAXlEO5u",
        "outputId": "89eb2b6b-d174-4fff-f408-cd8c086272a8"
      },
      "source": [
        "pip install xpress"
      ],
      "execution_count": 1,
      "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 95kB/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": "j3CQoBFzC_VJ",
        "outputId": "0961a4a1-44ae-413c-83f1-d1b1d9e9754c"
      },
      "source": [
        "import xpress as xp\n",
        "print(xp.manual())\n",
        "print(xp.examples())"
      ],
      "execution_count": 2,
      "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",
            "/usr/local/lib/python3.6/dist-packages/xpress/doc/python-interface.pdf\n",
            "/usr/local/lib/python3.6/dist-packages/xpress/examples/\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "t8waEB1nEdIk"
      },
      "source": [
        "x = xp.var(name='myvar_1')\n",
        "y = xp.var(name='myvar_2') #variables by defualt are non-negative"
      ],
      "execution_count": 3,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "E7huOFkIEwMR"
      },
      "source": [
        "constr1 = 2*x + 3*y >= 6\n",
        "constr2 = 4*x + 2*y >= 7\n",
        "obj = x + y"
      ],
      "execution_count": 4,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "LahAunS5E95l"
      },
      "source": [
        "p = xp.problem(x,y, #<-- variabes(always first to be added)\n",
        "               constr1,constr2,   #<-- constraints and objectives (in no order)\n",
        "               obj)"
      ],
      "execution_count": 5,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "IXDp5tQlFayF",
        "outputId": "c1cd109d-cd14-4b06-c43e-fa25459cdcc3"
      },
      "source": [
        "p.solve()\n",
        "print(f'solution: x={p.getSolution(x)}, y ={p.getSolution(y)}')"
      ],
      "execution_count": 6,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "FICO Xpress v8.11.0, Community, solve started 8:39:37, Nov 24, 2020\n",
            "Heap usage: 338KB (peak 338KB, 593KB system)\n",
            "Minimizing LP noname with these control settings:\n",
            "OUTPUTLOG = 1\n",
            "\n",
            "Original problem has:\n",
            "         2 rows            2 cols            4 elements\n",
            "Presolved problem has:\n",
            "         2 rows            2 cols            4 elements\n",
            "Presolve finished in 0 seconds\n",
            "Heap usage: 339KB (peak 351KB, 595KB system)\n",
            "\n",
            "Coefficient range                    original                 solved        \n",
            "  Coefficients   [min,max] : [ 2.00e+00,  4.00e+00] / [ 5.00e-01,  1.50e+00]\n",
            "  RHS and bounds [min,max] : [ 6.00e+00,  7.00e+00] / [ 1.75e+00,  3.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",
            " \n",
            "   Its         Obj Value      S   Ninf  Nneg   Sum Dual Inf  Time\n",
            "     0           .000000      D      2     0        .000000     0\n",
            "     2          2.375000      D      0     0        .000000     0\n",
            "Uncrunching matrix\n",
            "Optimal solution found\n",
            "Dual solved problem\n",
            "  2 simplex iterations in 0s\n",
            "\n",
            "Final objective                       : 2.375000000000000e+00\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",
            "solution: x=1.125, y =1.25\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "MKp3O-gSFfB9"
      },
      "source": [
        "#add constraint x + 2y <= 2 to the problem and then resolve. The problem will now be infeasible"
      ],
      "execution_count": 7,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "_HH3v8wmJavJ",
        "outputId": "fd688e9d-f970-4d23-cbd8-043bcdda2578"
      },
      "source": [
        "constr3 = x + 2*y <= 2\n",
        "\n",
        "p.addConstraint(constr3)\n",
        "\n",
        "p.solve()"
      ],
      "execution_count": 8,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "FICO Xpress v8.11.0, Community, solve started 8:39:37, Nov 24, 2020\n",
            "Heap usage: 2605KB (peak 2607KB, 595KB system)\n",
            "Minimizing LP noname with these control settings:\n",
            "OUTPUTLOG = 1\n",
            "\n",
            "Original problem has:\n",
            "         3 rows            2 cols            6 elements\n",
            " \n",
            " \n",
            "The problem is infeasible due to row R3\n",
            "Presolve finished in 0 seconds\n",
            "Heap usage: 2605KB (peak 2607KB, 595KB system)\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "GOUdxSd_JkJP"
      },
      "source": [
        "#use the iisall function to understand what made the problem infeasible. If needed, type help(p.iisall) for its docstring."
      ],
      "execution_count": 9,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "EgA6MJ24MEVq",
        "outputId": "6629a612-bf45-473f-e8ca-ae24409619bd"
      },
      "source": [
        "p.iisall()"
      ],
      "execution_count": 10,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "   Name      Type    Sense    Bound \n",
            " R1          row     GE      6.000000\n",
            " R3          row     LE      2.000000\n",
            " myvar_2     column  LO       .000000\n",
            "Found IIS 1 (2 rows and 1 bounds)\n",
            " \n",
            "No more IIS can be found\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "s2wCFFXIOV-6"
      },
      "source": [
        "Now delete the constraint (using the `problem.delConstraint` method) and, just to avoid the same solve as before, change the objective function coefficient of $x$ to 4. Then re-solve and print primal and dual solution using `problem.getSolution` for the primal solution and `problem.getDual` for the dual solution. Recall that these functions can be called without arguments to obtain the whole vector, or with a variable/constraint object to get one or more values."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "3Kxngr_EOUdd",
        "outputId": "fd34babf-1a9f-4493-f357-72060b00000e"
      },
      "source": [
        "p.delConstraint(constr3)\n",
        "p.chgobj([x],[4])\n",
        "p.solve()\n",
        "\n",
        "print(f'Solution: x={p.getSolution(x)}, y={p.getSolution(y)}')\n",
        "print(f'Dual: constr1->{p.getDual(constr1)}, constr2->{p.getDual(constr2)}')"
      ],
      "execution_count": 11,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "FICO Xpress v8.11.0, Community, solve started 8:39:37, Nov 24, 2020\n",
            "Heap usage: 4956KB (peak 4956KB, 832KB system)\n",
            "Minimizing LP noname with these control settings:\n",
            "OUTPUTLOG = 1\n",
            "\n",
            "Original problem has:\n",
            "         2 rows            2 cols            4 elements\n",
            "Presolved problem has:\n",
            "         2 rows            2 cols            4 elements\n",
            "Presolve finished in 0 seconds\n",
            "Heap usage: 4956KB (peak 4956KB, 832KB system)\n",
            "\n",
            "Coefficient range                    original                 solved        \n",
            "  Coefficients   [min,max] : [ 2.00e+00,  4.00e+00] / [ 5.00e-01,  1.50e+00]\n",
            "  RHS and bounds [min,max] : [ 6.00e+00,  7.00e+00] / [ 1.75e+00,  3.00e+00]\n",
            "  Objective      [min,max] : [ 1.00e+00,  4.00e+00] / [ 1.00e+00,  4.00e+00]\n",
            "Autoscaling applied standard scaling\n",
            "\n",
            " \n",
            "   Its         Obj Value      S   Ninf  Nneg   Sum Dual Inf  Time\n",
            "     0          5.750000      D      0     1       1.000000     0\n",
            "     1          3.500000      D      0     0        .000000     0\n",
            "Uncrunching matrix\n",
            "Optimal solution found\n",
            "Dual solved problem\n",
            "  1 simplex iterations in 0s\n",
            "\n",
            "Final objective                       : 3.500000000000000e+00\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",
            "Solution: x=0.0, y=3.5\n",
            "Dual: constr1->-0.0, constr2->0.5\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "0i-bgBRyQoYW"
      },
      "source": [
        "#Knapsack problem:\n",
        "\n",
        "Formulate and solve a knapsack problem with the following vaue/weight vectors:\n",
        "\n",
        ">$ v = (12, 15, 9, 11, 8, 7, 5)$ \\\\\n",
        ">$ w = (13, 18, 9, 12, 8, 10, 4)$ \\\\\n",
        "\n",
        "with knapsack capacity $ C = 40$ \\\\\n",
        "You will need `xpress.Sum` operator. The formuation of a knapsack problem is as follows:\n",
        "\n",
        "$$\n",
        "\\begin{array}{lllrrrrrrr}\n",
        "\\max & \\sum_{i=1}^n v_i x_i\\\\\n",
        "\\textrm{s.t.} & \\sum_{i=1}^n w_i x_i \\le C\\\\\n",
        "& x_i \\in \\{0,1\\} \\forall i=1,\\ldots{},n\n",
        "\\end{array}\n",
        "$$\n",
        "\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "nI9xvx7UPZ4F",
        "outputId": "b1057e20-44dc-426f-8b27-6b9c9fd6a475"
      },
      "source": [
        "v = [12,15,9,11,8,7,5]\n",
        "w = [13,18,9,12,8,10,4]\n",
        "C = 40\n",
        "\n",
        "n = len(v)\n",
        "x = [xp.var(vartype = xp.binary) for _ in range (n)]\n",
        "\n",
        "k_con = xp.Sum(w[i]*x[i] for i in range (n)) <= C\n",
        "k_obj = xp.Sum(v[i]*x[i] for i in range (n))\n",
        "\n",
        "p = xp.problem(x,k_con,k_obj, sense = xp.maximize) # optimization sense is min by default\n",
        "\n",
        "p.solve()\n",
        "\n",
        "print(p.getSolution())"
      ],
      "execution_count": 12,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "\n",
            "FICO Xpress v8.11.0, Community, solve started 8:39:37, Nov 24, 2020\n",
            "Heap usage: 338KB (peak 338KB, 562KB system)\n",
            "Maximizing MILP noname with these control settings:\n",
            "OUTPUTLOG = 1\n",
            "\n",
            "Original problem has:\n",
            "         1 rows            7 cols            7 elements         7 globals\n",
            "Presolved problem has:\n",
            "         1 rows            7 cols            7 elements         7 globals\n",
            "Presolve finished in 0 seconds\n",
            "Heap usage: 364KB (peak 369KB, 564KB system)\n",
            "\n",
            "Coefficient range                    original                 solved        \n",
            "  Coefficients   [min,max] : [ 4.00e+00,  1.80e+01] / [ 2.50e-01,  1.12e+00]\n",
            "  RHS and bounds [min,max] : [ 1.00e+00,  4.00e+01] / [ 1.00e+00,  2.50e+00]\n",
            "  Objective      [min,max] : [ 5.00e+00,  1.50e+01] / [ 5.00e+00,  1.50e+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:      .000000      Time: 0 ***\n",
            " *** Heuristic solution found:    31.000000      Time: 0 ***\n",
            "Starting concurrent solve with dual\n",
            "\n",
            " Concurrent-Solve,   0s\n",
            "            Dual        \n",
            "    objective   dual inf\n",
            " D  39.500000   .0000000\n",
            "------- optimal --------\n",
            "Concurrent statistics:\n",
            "      Dual: 1 simplex iterations, 0.00s\n",
            "Optimal solution found\n",
            " \n",
            "   Its         Obj Value      S   Ninf  Nneg   Sum Dual Inf  Time\n",
            "     1         39.500000      D      0     0        .000000     0\n",
            "Dual solved problem\n",
            "  1 simplex iterations in 0s\n",
            "\n",
            "Final objective                       : 3.950000000000000e+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           34.000000    39.500000      3                 13.92%       0      0\n",
            "a           37.000000    38.375000      4                  3.58%       0      0\n",
            "   1  K     37.000000    38.375000      4      2      0    3.58%       1      0\n",
            " \n",
            "Cuts in the matrix         : 2\n",
            "Cut elements in the matrix : 14\n",
            "\n",
            "Performing root presolve...\n",
            "\n",
            "Reduced problem has:       3 rows       5 columns        15 elements\n",
            "Presolve dropped   :       0 rows       2 columns         6 elements\n",
            "Presolve tightened :         3 elements\n",
            "Will try to keep branch and bound tree memory usage below 11.9GB\n",
            " \n",
            "   Its         Obj Value      S   Ninf  Nneg   Sum Dual Inf  Time\n",
            "     4         37.666667      D      0     0        .000000     0\n",
            "Problem is cut off\n",
            " *** Search completed ***\n",
            "Final MIP objective                   : 3.700000000000000e+01\n",
            "Final MIP bound                       : 3.700000000000000e+01\n",
            "  Solution time / primaldual integral :         0s/ 45.660563%\n",
            "  Number of solutions found / nodes   :         4 /         0\n",
            "  Max primal violation      (abs/rel) :       0.0 /       0.0\n",
            "  Max integer violation     (abs    ) :       0.0\n",
            "[1.0, 0.0, 1.0, 1.0, -0.0, -0.0, 1.0]\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "RUGkb4jMb4Bc"
      },
      "source": [
        "### NumPy operators and arrays of variables/constraints\n",
        "\n",
        "NumPy is an essential toolbox for Python users. The Xpress Python interface can handle (multi-)arrays of floats, variables, expressions, constraints as naturally as with lists thereof."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "_r9zD5SGZvuf"
      },
      "source": [
        "import numpy as np\n",
        "\n",
        "#The folowing two statements create two equivalent NumPy arrays of variables\n",
        "\n",
        "x = np.array([xp.var(vartype = xp.binary, name = 'x({})'.format(i)) for i in range(20)], dtype = xp.npvar)\n",
        "y = xp.vars(20,     vartype = xp.binary, name = 'x')"
      ],
      "execution_count": 13,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "31V5giuis1IO"
      },
      "source": [
        "#Exercise:\n",
        "We can rewrite the knapsack problem using NumPy constructs"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "crj_yiWHcmad"
      },
      "source": [
        "import xpress as xp\n",
        "import numpy as np\n",
        "\n",
        "v = np.array([12,15,9,11,8, 7,5])\n",
        "w = np.array([13,18,9,12,8,10,4])\n",
        "C = 40\n",
        "\n",
        "n = len(v)\n",
        "x = xp.vars(n, vartype=xp.binary)"
      ],
      "execution_count": 14,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "KSaV6ZcJm8SQ",
        "outputId": "342b71ae-437d-4464-906b-c0631edb5f9a"
      },
      "source": [
        "# xp.Dot only works with numpy arrays\n",
        "k_con = xp.Dot(w,x) <= C\n",
        "k_obj = xp.Dot(v,x)\n",
        "\n",
        "p = xp.problem(x, k_con, k_obj, sense=xp.maximize)\n",
        "\n",
        "p.solve()\n",
        "\n",
        "print(p.getSolution())"
      ],
      "execution_count": 15,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "\n",
            "FICO Xpress v8.11.0, Community, solve started 8:39:38, Nov 24, 2020\n",
            "Heap usage: 338KB (peak 338KB, 623KB system)\n",
            "Maximizing MILP noname with these control settings:\n",
            "OUTPUTLOG = 1\n",
            "\n",
            "Original problem has:\n",
            "         1 rows            7 cols            7 elements         7 globals\n",
            "Presolved problem has:\n",
            "         1 rows            7 cols            7 elements         7 globals\n",
            "Presolve finished in 0 seconds\n",
            "Heap usage: 364KB (peak 369KB, 625KB system)\n",
            "\n",
            "Coefficient range                    original                 solved        \n",
            "  Coefficients   [min,max] : [ 4.00e+00,  1.80e+01] / [ 2.50e-01,  1.12e+00]\n",
            "  RHS and bounds [min,max] : [ 1.00e+00,  4.00e+01] / [ 1.00e+00,  2.50e+00]\n",
            "  Objective      [min,max] : [ 5.00e+00,  1.50e+01] / [ 5.00e+00,  1.50e+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:      .000000      Time: 0 ***\n",
            " *** Heuristic solution found:    31.000000      Time: 0 ***\n",
            "Starting concurrent solve with dual\n",
            "\n",
            " Concurrent-Solve,   0s\n",
            "            Dual        \n",
            "    objective   dual inf\n",
            " D  39.500000   .0000000\n",
            "------- optimal --------\n",
            "Concurrent statistics:\n",
            "      Dual: 1 simplex iterations, 0.00s\n",
            "Optimal solution found\n",
            " \n",
            "   Its         Obj Value      S   Ninf  Nneg   Sum Dual Inf  Time\n",
            "     1         39.500000      D      0     0        .000000     0\n",
            "Dual solved problem\n",
            "  1 simplex iterations in 0s\n",
            "\n",
            "Final objective                       : 3.950000000000000e+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           34.000000    39.500000      3                 13.92%       0      0\n",
            "a           37.000000    38.375000      4                  3.58%       0      0\n",
            "   1  K     37.000000    38.375000      4      2      0    3.58%       1      0\n",
            " \n",
            "Cuts in the matrix         : 2\n",
            "Cut elements in the matrix : 14\n",
            "\n",
            "Performing root presolve...\n",
            "\n",
            "Reduced problem has:       3 rows       5 columns        15 elements\n",
            "Presolve dropped   :       0 rows       2 columns         6 elements\n",
            "Presolve tightened :         3 elements\n",
            "Will try to keep branch and bound tree memory usage below 11.9GB\n",
            " \n",
            "   Its         Obj Value      S   Ninf  Nneg   Sum Dual Inf  Time\n",
            "     4         37.666667      D      0     0        .000000     0\n",
            "Problem is cut off\n",
            " *** Search completed ***\n",
            "Final MIP objective                   : 3.700000000000000e+01\n",
            "Final MIP bound                       : 3.700000000000000e+01\n",
            "  Solution time / primaldual integral :         0s/ 48.942165%\n",
            "  Number of solutions found / nodes   :         4 /         0\n",
            "  Max primal violation      (abs/rel) :       0.0 /       0.0\n",
            "  Max integer violation     (abs    ) :       0.0\n",
            "[1.0, 0.0, 1.0, 1.0, -0.0, -0.0, 1.0]\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "cShHXCj5n7tw"
      },
      "source": [
        "#Exercise:\n",
        "\n",
        "Create a random 20x30 numpy matrix $A$ and vectors $b$ and $c$, then create a problem with $c$ as objective function coefficient vector and constraints $Ax \\leq b.$"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "1KUdBw0OnT81",
        "outputId": "880e4a07-e621-4079-cc64-c2fe070f94d2"
      },
      "source": [
        "n = 30 #number of variables\n",
        "m = 20 #number of constraints\n",
        "\n",
        "A = np.random.random((m,n))\n",
        "b = np.random.random(m)\n",
        "c = np.random.random(n)\n",
        "\n",
        "x = xp.vars(n)\n",
        "\n",
        "objective = xp.Dot(c, x)\n",
        "consys = xp.Dot(A, x) <= b\n",
        "\n",
        "p = xp.problem(x,objective, consys, sense = xp.maximize)\n",
        "\n",
        "p.solve()"
      ],
      "execution_count": 16,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "\n",
            "FICO Xpress v8.11.0, Community, solve started 8:39:38, Nov 24, 2020\n",
            "Heap usage: 356KB (peak 356KB, 625KB system)\n",
            "Maximizing LP noname with these control settings:\n",
            "OUTPUTLOG = 1\n",
            "\n",
            "Original problem has:\n",
            "        20 rows           30 cols          600 elements\n",
            "Presolved problem has:\n",
            "        20 rows           30 cols          600 elements\n",
            "Presolve finished in 0 seconds\n",
            "Heap usage: 367KB (peak 389KB, 627KB system)\n",
            "\n",
            "Coefficient range                    original                 solved        \n",
            "  Coefficients   [min,max] : [ 5.58e-04, 10.00e-01] / [ 5.58e-04, 10.00e-01]\n",
            "  RHS and bounds [min,max] : [ 6.08e-03,  9.81e-01] / [ 6.08e-03,  9.81e-01]\n",
            "  Objective      [min,max] : [ 6.73e-02,  9.61e-01] / [ 6.73e-02,  9.61e-01]\n",
            "Autoscaling applied standard scaling\n",
            "\n",
            " \n",
            "   Its         Obj Value      S   Ninf  Nneg   Sum Dual Inf  Time\n",
            "     0           .387407      D     12     0        .000000     0\n",
            "     5           .069425      D      0     0        .000000     0\n",
            "Uncrunching matrix\n",
            "Optimal solution found\n",
            "Dual solved problem\n",
            "  5 simplex iterations in 0s\n",
            "\n",
            "Final objective                       : 6.942508063265417e-02\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"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "hEmok8vvqt9B"
      },
      "source": [
        "#Exercise:\n",
        "\n",
        "Make the objective function quadratic, with a matrix $Q$ that is PSD but not necessarily diagonal. Make the objective\n",
        "\n",
        "$(x-x_0)^TQ(x-x_0)$"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "DTtTHtKCo0-H",
        "outputId": "3cfdcb2a-aeaa-4d38-b031-7fcc1d183f9e"
      },
      "source": [
        "#The following has a strong-enough diagonal that we're sure it is PSD\n",
        "Q = np.random.random((n,n)) - .5 + 3*n*np.eye(n)\n",
        "\n",
        "x0 = 10 * np.random.random(n) #creates a point with random coordinates in [0,10]\n",
        "\n",
        "#minimized weighted distance between x and x0: (x-x0)' Q(x=x0)\n",
        "\n",
        "obj2 = xp.Dot((x-x0), Q, (x-x0))\n",
        "\n",
        "p.setObjective(obj2, sense = xp.minimize)\n",
        "p.solve()"
      ],
      "execution_count": 17,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "FICO Xpress v8.11.0, Community, solve started 8:39:38, Nov 24, 2020\n",
            "Heap usage: 2652KB (peak 2652KB, 627KB system)\n",
            "Minimizing QP noname with these control settings:\n",
            "OUTPUTLOG = 1\n",
            "\n",
            "Original problem has:\n",
            "        20 rows           30 cols          600 elements\n",
            "       900 qobjelem\n",
            "Presolved problem has:\n",
            "        20 rows           30 cols          600 elements\n",
            "       900 qobjelem\n",
            "Presolve finished in 0 seconds\n",
            "Heap usage: 2685KB (peak 2685KB, 627KB system)\n",
            "\n",
            "Coefficient range                    original                 solved        \n",
            "  Coefficients   [min,max] : [ 5.58e-04, 10.00e-01] / [ 1.31e-04,  1.99e+00]\n",
            "  RHS and bounds [min,max] : [ 6.08e-03,  9.81e-01] / [ 6.08e-03,  9.81e-01]\n",
            "  Objective      [min,max] : [ 3.44e+01,  1.61e+03] / [ 4.30e+00,  7.41e+02]\n",
            "  Quadratic      [min,max] : [ 4.25e-04,  1.81e+02] / [ 8.50e-04,  7.21e+02]\n",
            "Autoscaling applied standard scaling\n",
            "\n",
            "Barrier cache sizes : L1=32K L2=46080K\n",
            "Using AVX support\n",
            "Cores per CPU (CORESPERCPU): 2\n",
            "Barrier starts after 0 seconds, using up to 2 threads, 1 cores\n",
            "Matrix ordering - Dense cols.:     30   NZ(L):      1325   Flops:        41650\n",
            " \n",
            "  Its   P.inf      D.inf      U.inf      Primal obj.     Dual obj.      Compl.\n",
            "   0   1.81e+01   9.40e+02   0.00e+00   6.0499606e+04   6.8339408e+04   9.9e+03\n",
            "   1   8.53e+00   6.53e+02   0.00e+00   6.1213730e+04   7.0633616e+04   9.1e+03\n",
            "   2   2.84e+00   4.45e+01   0.00e+00   6.9137490e+04   7.0817343e+04   8.2e+03\n",
            "   3   1.51e+00   2.36e+01   0.00e+00   7.0372216e+04   7.1105355e+04   4.5e+03\n",
            "   4   8.10e-01   1.27e+01   0.00e+00   7.1027410e+04   7.1340105e+04   2.4e+03\n",
            "   5   4.34e-01   6.80e+00   0.00e+00   7.1388298e+04   7.1522335e+04   1.3e+03\n",
            "   6   2.31e-01   3.62e+00   0.00e+00   7.1581547e+04   7.1641047e+04   7.2e+02\n",
            "   7   1.54e-02   2.42e-01   0.00e+00   7.1787112e+04   7.1783307e+04   6.0e+01\n",
            "   8   1.32e-02   2.08e-01   0.00e+00   7.1798520e+04   7.1808484e+04   6.4e+01\n",
            "   9   1.00e-02   1.56e-01   0.00e+00   7.1801340e+04   7.1835890e+04   5.9e+01\n",
            "  10   3.96e-03   6.11e-02   0.00e+00   7.1866331e+04   7.1839353e+04   7.3e+01\n",
            "  11   7.23e-06   1.09e-04   0.00e+00   7.1871076e+04   7.1870607e+04   5.6e-01\n",
            "Barrier method finished in 0 seconds\n",
            "Uncrunching matrix\n",
            "Optimal solution found\n",
            " \n",
            "   Its         Obj Value      S   Ninf  Nneg        Sum Inf  Time\n",
            "     0       71871.07562      B      2     0        .000000     0\n",
            "Barrier solved problem\n",
            "  11 barrier iterations in 0s\n",
            "\n",
            "Final objective                       : 7.187107562149914e+04\n",
            "  Max primal violation      (abs/rel) : 5.316e-06 / 5.316e-06\n",
            "  Max dual violation        (abs/rel) :       0.0 /       0.0\n",
            "  Max complementarity viol. (abs/rel) : 9.068e-02 / 5.341e-06\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Pll-EDGMtFHQ"
      },
      "source": [
        "Obviously we can use arrays of variables in constraints. Let's consider the following problem:\n",
        "\n",
        "### Production planning with quadratic cost\n",
        "\n",
        "A factory must plan production on one machine for the upcoming $N=12$ days. If, on one day, the machine is on, it must pay a fixed cost F, regardless of how much is produced on that machine. If the machine is on, the production must range between L and U. It must also pay a cost that is C times the square of the produced amount.\n",
        "\n",
        "The demand at a given day is $w_i$, and the unsold product can be stored in a storage container of capacity $M$. Finally, the storage level at the end of the $N$ days must be the same as the beginning, which is set as $S$.\n",
        "\n",
        "**Model**. Define two classes of variables, $x$ (continuous) and $y$ (binary), indicating the amount of production in a given day and whether the machine is on or off, respectively. Also, introduce a storage variable $s$ that is also indexed by the set of periods."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "qAboo7omsw3m",
        "outputId": "4af540dc-f7fa-4fa3-f4f6-29e21f08401d"
      },
      "source": [
        "n = 12 #periods\n",
        "N = range(n) \n",
        "\n",
        "C = 2 #proportional cost\n",
        "L = 10 #minimum production level\n",
        "U = 90 #maximum production level\n",
        "F = 9000 #Fixed cost\n",
        "M = 170 #maximum storage level\n",
        "S = 25 #initial storage level\n",
        "\n",
        "#demand\n",
        "w = [60, 15, 25, 70, 70, 85, 10, 5, 65, 40, 50, 15]\n",
        "\n",
        "#variables\n",
        "x = [xp.var() for _ in N]\n",
        "y = [xp.var(vartype = xp.binary) for _ in N]\n",
        "s = [xp.var(ub=M) for _ in range(n+1)] #can assign the upper bound of M immediately\n",
        "\n",
        "#total production cost (to be minimized)\n",
        "\n",
        "objective = xp.Sum(F*y[i] + C*y[i]*x[i]**2 for i in N)\n",
        "\n",
        "#y determines bounds on x\n",
        "x_lb = [x[i] >= L*y[i] for i in N]\n",
        "x_ub = [x[i] <= U*y[i] for i in N]\n",
        "\n",
        "# storage at last period\n",
        "\n",
        "store_end = s[n] == S\n",
        "store_ini = s[0] == S\n",
        "\n",
        "# Production balance\n",
        "\n",
        "balance = [s[i] + x[i] == s[i+1] + w[i] for i in N]  # Conservation at each time period\n",
        "\n",
        "p = xp.problem(x, y, s,\n",
        "               objective,\n",
        "               x_lb, x_ub, store_end, store_ini, balance)\n",
        "\n",
        "p.solve()\n",
        "\n",
        "print('time demand machine productn storage')\n",
        "for i in N:\n",
        "    print('{0:4d} {1:6d} {2:7.0f} {3:8.4f} {4:7.4f}'.format(i, w[i],\n",
        "                                                            p.getSolution(y[i]),\n",
        "                                                            p.getSolution(x[i]),\n",
        "                                                            p.getSolution(s[i])))\n"
      ],
      "execution_count": 18,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "\n",
            "Original problem size\n",
            "   linear:    39 rows, 39 columns, 87 linear coefficients\n",
            "   nonlinear: 1 coefficients, 95 tokens\n",
            "Nonlinear presolve\n",
            "   compressing formula space (in use\\total : 96\\713)\n",
            "   removed 22 formulas and 616 formula tokens\n",
            "   linear row reductions removed 2 rows, 2 columns, 4 linear coefficients\n",
            "   bound tightening reduced 16 bounds\n",
            "   objective transfer row objtransferrow through variable objtransfercol\n",
            "Presolved problem size\n",
            "   linear:    37 rows, 37 columns, 83 linear coefficients\n",
            "   nonlinear: 1 coefficients, 95 tokens\n",
            "Problem is nonlinear presolved\n",
            "FICO Xpress v8.11.0, Community, solve started 8:39:38, Nov 24, 2020\n",
            "Control settings used:\n",
            "XSLP_POSTSOLVE = 1\n",
            "XSLP_DELETIONCONTROL = 0\n",
            "OUTPUTLOG = 1\n",
            "\n",
            "Maximum expanded nl-formula size: 95  (row 'objtransferrow')\n",
            "Total tokens: 95\n",
            "  2  parallel calculation threads\n",
            "  Jacobian: symbolic differentiation\n",
            "           1 base AD formula, 96 average complexity\n",
            "          24 in the Jacobian, 6 average complexity\n",
            "Initial point objective:  157112.50\n",
            "Absolute / relative validation:            .000  /      .904\n",
            "Minimizing problem using Xpress-SLP\n",
            "Xpress-SLP Augmentation Statistics:\n",
            "  Columns:\n",
            "          24 delta vectors\n",
            "  Rows:\n",
            "           1 nonlinear constraints\n",
            "          24 update rows\n",
            "  Coefficients:\n",
            "          25 non-constant coefficients\n",
            "\n",
            " It LP    NetObj   ValObj ErrorSum ErrorCost Validate   KKT Unconv  Ext Action T\n",
            "  1 O  65993.750  110550.      .00      .00      .00            20    0 *      0\n",
            "  2 O -49916.665  126200.      .00      .00      .00            24    0        0\n",
            "  3 O -68299.996  127350.      .00      .00      .00            24    0        0\n",
            "  4 O -76649.996  131000.      .00      .00      .00            24    0        0\n",
            "  5 O -86299.996  132750.      .00      .00      .00            24    0        0\n",
            "  6 O -88049.996  131000.      .00      .00      .00            24    0        0\n",
            "  7 O -86299.996  132750.      .00      .00      .00            24    0        0\n",
            "  8 O -88049.996  131000.      .00      .00      .00            24    0        0\n",
            "  9 O  22320.001  77098.0      .00      .00      .00            24    0 *      0\n",
            " 10 O  73265.175  84585.7      .00      .00      .00            24    0        0\n",
            " 11 O  73980.150  77718.7      .00      .00      .00            24    0        0\n",
            " 12 O  74225.764  75388.9      .00      .00      .00            24    0 *      0\n",
            " 13 O  73574.384  74348.2      .00      .00      .00            24    0 *      0\n",
            " 14 O  73219.216  73734.5      .00      .00      .00            23    0 *      0\n",
            " 15 O  73347.042  73884.6      .00      .00      .00            24    0        0\n",
            " 16 O  73087.658  73673.7      .00      .00      .00            22    0 *      0\n",
            " 17 O  73343.469  73526.9      .00      .00      .00            24    0 *      0\n",
            " 18 O  73417.923  73495.1      .00      .00      .00            24    0 *      0\n",
            " 19 O  73386.680  73422.8      .00      .00      .00            22    0 *      0\n",
            " 20 O  73377.645  73393.4      .00      .00      .00             0   10 *      0\n",
            " 21 O  73376.377  73393.8      .00      .00      .00             0   10        0\n",
            " 22 O  73362.257  73405.2      .00      .00      .00            21    0        0\n",
            " 23 O  73370.005  73384.7      .00      .00      .00             0    9 *      0\n",
            " 24 O  73380.644  73384.8      .00      .00      .00     2.56    0    9 K      0\n",
            " 25 O  73380.413  73384.5      .00      .00      .00     2.50    0    7 K*     0\n",
            " 26 O  73382.576  73383.4      .00      .00      .00      .32    0    7 K*     0\n",
            " 27 O  73383.205  73383.4      .00      .00      .00      .35    0    4 *      0\n",
            "Returning final converged solution\n",
            "\n",
            "Xpress-SLP stopped after 27 iterations. 0 unconverged items\n",
            "Problem solved using Xpress-NLP SLP\n",
            "Converged on objective function with some variables converged on extended criteria only. 0 unconverged items\n",
            "  4 step-bounded items converged on extended criteria only\n",
            "No unconverged values in active constraints\n",
            "FICO Xpress v8.11.0, Community, solve started 8:39:38, Nov 24, 2020\n",
            "Heap usage: 114KB (peak 118KB, 1545KB system)\n",
            "Minimizing MILP noname\n",
            "Original problem has:\n",
            "        61 rows           61 cols          156 elements        12 globals\n",
            "\n",
            "Coefficient range                    original                 solved        \n",
            "  Coefficients   [min,max] : [ 1.00e+00,  2.24e+04] / [ 1.32e-02,  1.41e+00]\n",
            "  RHS and bounds [min,max] : [ 5.49e-03,  1.70e+02] / [ 5.49e-03,  1.70e+02]\n",
            "  Objective      [min,max] : [ 1.00e+00,  9.00e+03] / [ 4.10e+03,  9.00e+03]\n",
            "Autoscaling applied standard scaling\n",
            "\n",
            "Will try to keep branch and bound tree memory usage below 11.8GB\n",
            " *** Search unfinished ***    Time:     0 Nodes:          0\n",
            "  Solution time / primaldual integral :         0s/ 100.000000%\n",
            "  Number of solutions found / nodes   :         0 /         0\n",
            "\n",
            "Starting nonlinear heuristics\n",
            "S     0  141700.00  73383.205   1      0     0    48.21%      0       41      0\n",
            "Heuristics stopped      .069s\n",
            "Setting tree search absolute cutoff to 141700\n",
            "\n",
            "Invoking SLP global search\n",
            "MIP strategy: SLP within MIP\n",
            "\n",
            "\n",
            "Using 2 parallel tree threads\n",
            "\n",
            "   Node  BestSoln   BoundEstim Sols Active Depth   Gap     GInf    SLPIts  Time\n",
            "      1  141700.00  73383.205   1      2     1    48.21%     12       46      0\n",
            "      2  141700.00  73383.205   1      2     2    48.21%     11       66      0\n",
            "      3  141700.00  73383.205   1      2     3    48.21%     10       86      0\n",
            "      4  141700.00  73383.205   1      2     4    48.21%      9      108      0\n",
            "      5  141700.00  73383.205   1      5     2    48.21%     11      128      0\n",
            "      6  141700.00  73383.205   1      5     3    48.21%     10      148      0\n",
            "      7  141700.00  73383.205   1      5     4    48.21%      9      169      0\n",
            "      8  141700.00  73383.205   1      8     3    48.21%     10      188      0\n",
            "      9  141700.00  73383.205   1      8     4    48.21%      9      207      0\n",
            "     10  141700.00  73383.205   1      8     4    48.21%      9      228      0\n",
            "*    21  138767.34  73383.205   2     14    13    47.11%      0      480      0\n",
            "     21  138767.34  73383.205   2     14    13    47.11%      0      481      0\n",
            "     23  138767.34  73383.205   2     14    13    47.11%      1      481      0\n",
            "*    33  137313.12  73383.205   3     14    13    46.55%      0      693      0\n",
            "     33  137313.12  73383.205   3     14    13    46.55%      0      693      0\n",
            "     35  137313.12  73383.205   3     14    13    46.55%      1      694      0\n",
            "     45  137313.12  73383.205   3     14    12    46.55%      2      890      0\n",
            "     55  137313.12  73383.205   3     34    13    46.55%      1     1535      0\n",
            "     65  137313.12  73383.205   3     34     9    46.55%      5     1535      0\n",
            "*    69  137021.87  73383.205   4     34    13    46.44%      0     1881      0\n",
            "     69  137021.87  73383.205   4     34    13    46.44%      0     1881      0\n",
            "     76  137021.87  73383.205   4     34    13    46.44%      1     1921      0\n",
            "     89  137021.87  73383.205   4     34    12    46.44%      1     2043      0\n",
            "    100  137021.87  73383.205   4     34    13    46.44%      1     2365      0\n",
            "    200  137021.87  73383.205   4    112    10    46.44%      4     3993      1\n",
            "    302  137021.87  73383.205   4    112     9    46.44%      5     5774      1\n",
            "    402  137021.87  73383.205   4    127    13    46.44%      1     7291      2\n",
            "*   434  136993.34  73383.205   5    127    13    46.43%      0     7671      2\n",
            "    439  136993.34  73383.205   5    127    13    46.43%      0     7672      2\n",
            "    507  136993.34  73383.205   5    166    11    46.43%      3     8852      2\n",
            "*   555  136974.12  73383.205   6    186    13    46.42%      0     9401      2\n",
            "    569  136974.12  73383.205   6    186    13    46.42%      0     9401      2\n",
            "    613  136974.12  73383.205   6    186    13    46.42%      1    10344      3\n",
            "    724  136974.12  73383.205   6    186    13    46.42%      1    11689      3\n",
            "    828  136974.12  73383.205   6    190    13    46.42%      1    13042      3\n",
            "    936  136974.12  73383.205   6    189    13    46.42%      1    14368      4\n",
            "   1050  136974.12  73383.205   6    194    13    46.42%      1    16149      4\n",
            "   1150  136974.12  73383.205   6    192    13    46.42%      1    17674      5\n",
            "   1253  136974.12  73383.205   6    192    12    46.42%      2    18800      5\n",
            "   1360  136974.12  73383.205   6    192    11    46.42%      3    19775      5\n",
            "   1467  136974.12  73383.205   6    187    13    46.42%      1    21377      6\n",
            "   1571  136974.12  73383.205   6    184    11    46.42%      3    22601      6\n",
            "   1689  136974.12  73383.205   6    184    13    46.42%      1    23825      6\n",
            "   1801  136974.12  73383.205   6    184    13    46.42%      1    25537      7\n",
            "   1909  136974.12  73383.205   6    180    13    46.42%      1    26578      7\n",
            "   2021  136974.12  73383.205   6    171    13    46.42%      1    28406      7\n",
            "   2124  136974.12  73383.205   6    171    13    46.42%      1    28849      8\n",
            "   2229  136974.12  73383.205   6    167    13    46.42%      1    30514      8\n",
            "   2330  136974.12  73383.205   6    167    13    46.42%      1    31794      8\n",
            "   2432  136974.12  73383.205   6    155    13    46.42%      1    33610      9\n",
            "   2536  136974.12  73383.205   6    155    13    46.42%      1    34774      9\n",
            "   2646  136974.12  73383.205   6    151    12    46.42%      2    35914      9\n",
            "   2748  136974.12  73383.205   6    144    12    46.42%      2    37265     10\n",
            "   2854  136974.12  73383.205   6    141    13    46.42%      1    38333     10\n",
            "   2954  136974.12  73383.205   6    139    13    46.42%      1    39899     10\n",
            "   3060  136974.12  73383.205   6    132    13    46.42%      1    41137     11\n",
            "   3162  136974.12  73383.205   6    120    13    46.42%      1    42101     11\n",
            "   3268  136974.12  73383.205   6    116    13    46.42%      1    43770     11\n",
            "   3372  136974.12  73383.205   6    108    13    46.42%      1    44981     12\n",
            "   3482  136974.12  73383.205   6     97    13    46.42%      1    46025     12\n",
            "   3583  136974.12  73383.205   6     91    13    46.42%      1    47219     12\n",
            "   3692  136974.12  73383.205   6     82    13    46.42%      1    48724     13\n",
            "   3793  136974.12  73383.205   6     71    12    46.42%      1    49834     13\n",
            "   3899  136974.12  73383.205   6     45    13    46.42%      1    51214     13\n",
            "   4001  136974.12  73383.205   6     26    13    46.42%      1    52348     14\n",
            "\n",
            "\n",
            "\n",
            "\n",
            "Tree search: completed\n",
            "\n",
            "Integer solution statistics:\n",
            "  NLP heuristics: 1\n",
            "  MIP heuristics: 0\n",
            "  Tree search   : 5\n",
            "\n",
            "Problem is nonlinear postsolved\n",
            "Heap usage: 1102KB (peak 3358KB, 1187KB system)\n",
            "Observed Lipschitz constant:  88200.00\n",
            "\n",
            "Final MINLP objective (feas. solution)  : 1.370344882812500e+05\n",
            "  Max validation error      (abs/rel) :      .000 /      .000\n",
            "  Observed primal integral :     .539%\n",
            "  Total time :                14.475s \n",
            "*** Search completed ***\n",
            "time demand machine productn storage\n",
            "   0     60       1  63.7500 25.0000\n",
            "   1     15       0   0.0000 28.7500\n",
            "   2     25       1  63.3438 13.7500\n",
            "   3     70       1  63.6562 52.0938\n",
            "   4     70       1  64.4375 45.7500\n",
            "   5     85       1  64.0000 40.1875\n",
            "   6     10       0   0.0000 19.1875\n",
            "   7      5       1  65.0625  9.1875\n",
            "   8     65       1  63.5000 69.2500\n",
            "   9     40       1  62.2500 67.7500\n",
            "  10     50       0   0.0000 90.0000\n",
            "  11     15       0   0.0000 40.0000\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "aw2fo0Y4xuio"
      },
      "source": [
        "#Homework:\n",
        "Rewrite the problem using NumPy for defining parameters, variables, constraints and objective."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "HcrbTHU-xkPn",
        "outputId": "34834200-eb59-44a4-a187-c6e9a6c9f7cd"
      },
      "source": [
        "#demand\n",
        "\n",
        "W = np.array(w)\n",
        "\n",
        "#variables\n",
        "\n",
        "x = xp.vars(n, name = 'prod')\n",
        "y = xp.vars(n, vartype = xp.binary, name = 'machine')\n",
        "s = xp.vars(n + 1, ub = M, name = 'storage') \n",
        "\n",
        "#total production cost\n",
        "\n",
        "objective = xp.Sum(F*y + C*x**2)\n",
        "\n",
        "# y determines bounds on x\n",
        "\n",
        "x_lb = x >= L*y\n",
        "x_ub = x <= U*y\n",
        "\n",
        "#storage at last period\n",
        "\n",
        "store_end = s[n] == S\n",
        "store_ini = s[0] == S\n",
        "\n",
        "#production balance\n",
        "balance = s[:n] + x == s[1:] + w #conservation at each time period\n",
        "\n",
        "p = xp.problem(x,y,s,\n",
        "               objective,\n",
        "               x_lb,x_ub,store_end, store_ini, balance)\n",
        "\n",
        "p.solve()\n",
        "\n",
        "print('time demand machine productn storage')\n",
        "for i in N:\n",
        "    print('{0:4d} {1:6d} {2:7.0f} {3:8.4f} {4:7.4f}'.format(i, w[i],\n",
        "                                                            p.getSolution(y[i]),\n",
        "                                                            p.getSolution(x[i]),\n",
        "                                                            p.getSolution(s[i])))"
      ],
      "execution_count": 19,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "\n",
            "FICO Xpress v8.11.0, Community, solve started 8:39:52, Nov 24, 2020\n",
            "Heap usage: 362KB (peak 362KB, 917KB system)\n",
            "Minimizing MIQP noname with these control settings:\n",
            "OUTPUTLOG = 1\n",
            "\n",
            "Original problem has:\n",
            "        38 rows           37 cols           86 elements        12 globals\n",
            "        12 qobjelem\n",
            "Presolved problem has:\n",
            "        32 rows           32 cols           74 elements        11 globals\n",
            "        12 qobjelem\n",
            "LP relaxation tightened\n",
            "Presolve finished in 0 seconds\n",
            "Heap usage: 395KB (peak 402KB, 918KB system)\n",
            "\n",
            "Coefficient range                    original                 solved        \n",
            "  Coefficients   [min,max] : [ 1.00e+00,  9.00e+01] / [ 7.81e-03,  2.81e+00]\n",
            "  RHS and bounds [min,max] : [ 1.00e+00,  1.70e+02] / [ 1.00e+00,  3.60e+02]\n",
            "  Objective      [min,max] : [ 9.00e+03,  9.00e+03] / [ 3.50e+01,  9.00e+03]\n",
            "  Quadratic      [min,max] : [ 4.00e+00,  4.00e+00] / [ 2.50e-01,  2.50e-01]\n",
            "Autoscaling applied standard scaling\n",
            "\n",
            "Will try to keep branch and bound tree memory usage below 11.8GB\n",
            "Crash basis containing 10 structural columns created\n",
            " \n",
            "   Its         Obj Value      S   Ninf  Nneg        Sum Inf  Time\n",
            "     0       77950.00000      p      8     5    2690.000000     0\n",
            "    25       58950.00000      p      0     0        .000000     0\n",
            "    25       136600.0000      p      0     0        .000000     0\n",
            " \n",
            "   Its         Obj Value      S   Nsft  Nneg       Dual Inf  Time\n",
            "    41       101405.7292     QP      0     0        .000000     0\n",
            "QP solution found\n",
            "Optimal solution found\n",
            "Primal solved problem\n",
            "  41 simplex iterations in 0s\n",
            "\n",
            "Final objective                       : 1.014057291666667e+05\n",
            "  Max primal violation      (abs/rel) : 1.421e-14 / 1.421e-14\n",
            "  Max dual violation        (abs/rel) :       0.0 /       0.0\n",
            "  Max complementarity viol. (abs/rel) : 2.511e-12 / 6.974e-16\n",
            "\n",
            "Starting root cutting & heuristics\n",
            " \n",
            " Its Type    BestSoln    BestBound   Sols    Add    Del     Gap     GInf   Time\n",
            "k         152699.9928  101405.7292      1                 33.59%       0      0\n",
            "d         137314.2856  101405.7292      2                 26.15%       0      0\n",
            "   1  K   137314.2856  101552.5426      2      5      0   26.04%      11      0\n",
            "   2  K   137314.2856  103282.3281      2      3      4   24.78%      11      0\n",
            "   3  K   137314.2856  103299.8602      2      2      2   24.77%      11      0\n",
            "   4  K   137314.2856  103299.8602      2      0      1   24.77%      11      0\n",
            "Heuristic search started\n",
            "Heuristic search stopped\n",
            " \n",
            "Cuts in the matrix         : 3\n",
            "Cut elements in the matrix : 27\n",
            "\n",
            "Starting tree search.\n",
            "Deterministic mode with up to 2 running threads and up to 4 tasks.\n",
            "Heap usage: 3224KB (peak 5177KB, 936KB system)\n",
            " \n",
            "    Node     BestSoln    BestBound   Sols Active  Depth     Gap     GInf   Time\n",
            "       1  137314.2856  105294.7815      2      2      1   23.32%      11      0\n",
            "       2  137314.2856  105294.7815      2      2      3   23.32%       9      0\n",
            "       3  137314.2856  105294.7815      2      2      4   23.32%       8      0\n",
            "       4  137314.2856  105294.7815      2      2      5   23.32%       5      0\n",
            "       5  137314.2856  108809.6560      2      5      4   20.76%       8      0\n",
            "       6  137314.2856  108809.6560      2      5      5   20.76%       7      0\n",
            "       7  137314.2856  108809.6560      2      5      6   20.76%       6      0\n",
            "       8  137314.2856  110496.6531      2      8      3   19.53%       9      0\n",
            "       9  137314.2856  110496.6531      2      8      4   19.53%       8      0\n",
            "      10  137314.2856  110496.6531      2      8      5   19.53%       7      0\n",
            "c     21  137025.0000  111899.5454      3     14     12   18.34%       0      0\n",
            "      23  137025.0000  111899.5454      3     14     12   18.34%       1      0\n",
            "      42  137025.0000  111899.5454      3     14     12   18.34%       1      0\n",
            "      58  137025.0000  112477.8316      3     31     12   17.91%       1      0\n",
            "      73  137025.0000  112477.8316      3     31      9   17.91%       2      0\n",
            "      85  137025.0000  112477.8316      3     31     11   17.91%       2      0\n",
            "      99  137025.0000  112477.8316      3     31      9   17.91%       2      0\n",
            "     113  137025.0000  112477.8316      3     31     12   17.91%       1      0\n",
            "     215  137025.0000  118987.1875      3     80     11   13.16%       2      0\n",
            "Elapsed time (sec): 0, estimated tree completion: 0.18846\n",
            "Heap usage: 5568KB (peak 5568KB, 1925KB system)\n",
            "B&B tree size: 234k total\n",
            " \n",
            "    Node     BestSoln    BestBound   Sols Active  Depth     Gap     GInf   Time\n",
            "     319  137025.0000  118987.1875      3     99     10   13.16%       2      0\n",
            "     420  137025.0000  123313.8892      3     98     11   10.01%       2      0\n",
            "     528  137025.0000  123386.7778      3     94     11    9.95%       2      0\n",
            "     635  137025.0000  123983.5937      3     84     11    9.52%       2      0\n",
            "     735  137025.0000  126162.1876      3     76     10    7.93%       3      0\n",
            "     838  137025.0000  126773.5410      3     60     12    7.48%       1      0\n",
            "     940  137025.0000  129563.8891      3     39     11    5.45%       2      0\n",
            " *** Search completed ***\n",
            "Uncrunching matrix\n",
            "Final MIP objective                   : 1.370249999549753e+05\n",
            "Final MIP bound                       : 1.370248985492461e+05\n",
            "  Solution time / primaldual integral :         0s/ 13.520596%\n",
            "  Number of solutions found / nodes   :         3 /      1035\n",
            "  Max primal violation      (abs/rel) : 3.982e-08 / 3.982e-08\n",
            "  Max integer violation     (abs    ) :       0.0\n",
            "time demand machine productn storage\n",
            "   0     60       1  63.7509 25.0000\n",
            "   1     15       1  63.7502 28.7509\n",
            "   2     25       1  63.7502 77.5011\n",
            "   3     70       0   0.0000 116.2512\n",
            "   4     70       1  63.7501 46.2512\n",
            "   5     85       1  63.7500 40.0013\n",
            "   6     10       1  63.7496 18.7513\n",
            "   7      5       0   0.0000 72.5010\n",
            "   8     65       1  63.7496 67.5010\n",
            "   9     40       0   0.0000 66.2506\n",
            "  10     50       1  63.7494 26.2506\n",
            "  11     15       0   0.0000 40.0000\n"
          ],
          "name": "stdout"
        }
      ]
    }
  ]
}