{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"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 file is: `i,j,c`, where $i$ and $j$ are the endnodes of the directed $i\\to j$ arc, while $c$ is the cost of the arc."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import pandas as pd\n",
"import xpress as xp\n",
"\n",
"df = pd.read_csv('graph_data.csv')\n",
"\n",
"tail = df['i']\n",
"head = df['j']\n",
"cost = df['c']\n",
"\n",
"# 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(\"Arc costs:\", A)\n",
"print(\"Find shortest path from {0} to {1}\".format(source, dest))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Use the vars() function 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 incoming flow\n",
"# must equal the above right-hand side, 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\n",
"# 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)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"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 \\ge 4$ only if a third, binary variable $z$ is set to one, then this can be denoted as \n",
"\n",
"$$\n",
"z=1 \\Rightarrow x + 2y \\ge 4\n",
"$$\n",
"\n",
"Indicator constraints are common in Mixed Integer Linear Optimization. In the Xpress Python interface, the above indicator constraint would be added to a problem as follows:\n",
"\n",
"```\n",
"p.addIndicator (z == 1, x + 2*y >= 4)\n",
"```\n",
"\n",
"We can specify indicators with the `z==0` condition too, but these are the only two alternatives: constraints other than `z == 0` or `z == 1` as the first argument of an indicator will trigger an error.\n",
"\n",
"```\n",
"p.addIndicator (z == 0, x + 2*y >= 4)\n",
"```\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## General functions\n",
"\n",
"The Xpress Optimizer allows for modeling specific non-linear functions so that they are automatically reformulated within a Mixed Integer format. These functions are:\n",
"\n",
"* The absolute value of a function, $|x|$;\n",
"* The maximum and minimum of a list of functions, e.g. $\\max(x+1, 2y-2,...)$;\n",
"* The logical And and Or operators on a set of functions, for example $(x \\wedge y) \\vee z$; these operators work on binary variables.\n",
"\n",
"While the C interface of the Xpress optimizer uses low-level functions for this, the Python interface allows for natural modeling using the operators `xpress.abs`, `xpress.max`, `xpress.min`, `xpress.And`, `xpress.Or`, and the binary Python operators `&` and `|`. Let's see an example of this."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import xpress as xp\n",
"\n",
"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))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Manual MIP reformulation of $|x|$, $\\max/\\min$, and AND/OR\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": {},
"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": "markdown",
"metadata": {},
"source": [
"The Xpress Python interface allows for writing\n",
"\n",
"* `xpress.abs(2*x-4)` to create the expression $|2x-4|$;\n",
"* `xpress.min(x1, x2, 2*x3-1)` for the expression $\\min(x_1, x_2, 2x_3-1)$, and similar for the maximum;\n",
"* `xpress.And(y1, y2, y3, 1-y4)` for the logical And between the binary variables $y_1$, $y_2$, $y_3$, and the negation of $y_4$, which can be modeled as $1-y_4$ as it is 0 when $y_4=1$ and viceversa, functioning therefore as a negation operator. "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import xpress as xp\n",
"\n",
"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()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"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",
"execution_count": null,
"metadata": {},
"outputs": [],
"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",
"\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(x))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Simple cutting plane algorithm\n",
"\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",
"execution_count": null,
"metadata": {},
"outputs": [],
"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)) - .5\n",
"c = np.random.random(n) - .5\n",
"x = xp.vars(n, vartype=xp.binary)\n",
"\n",
"p = xp.problem(x, xp.Dot(x,Q,x) + xp.Dot(c,x))\n",
"\n",
"p.solve()\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"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",
"execution_count": null,
"metadata": {},
"outputs": [],
"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"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"And this is where we start from with our cutting plane algorithm. The following are well-known *valid inequalities* for the above problem, each holding for each triple of distinct indices $(i,j,k)$ where $i= x[i] + x[j] - 1 for i in range(n) for j in range(n)])\n",
"\n",
"n_violated = -1\n",
"eps = 0.001 # minimum violation for inequality to be entered\n",
"p.controls.outputlog = 0 # Silence optimizer output to make sure\n",
"MAXINEQ = 4*n\n",
"\n",
"# STUDENT EXERCISE (to be completed)\n",
"\n",
"while n_violated != 0:\n",
" p.solve()\n",
" print('inequalities: {0:5d}; lower bound: {1:10.3f}'.format(n_violated, p.getObjVal()))\n",
" n_violated = 0\n",
" ineq1_triples = []\n",
" ineq2_triples = []\n",
" ineq3_triples = []\n",
" ineq4_triples = []\n",
" \n",
" # Useful to see evolution of x values from all 0.5 to different values\n",
" # print(p.getSolution([x[i] for i in range(n)]))\n",
" \n",
" for i in range(n-2):\n",
" if n_violated >= MAXINEQ:\n",
" break\n",
" for j in range(i+1,n-1):\n",
" if n_violated >= MAXINEQ:\n",
" break\n",
" for k in range(j+1,n):\n",
" xibar = p.getSolution(x[i])\n",
" xjbar = p.getSolution(x[j])\n",
" xkbar = p.getSolution(x[k])\n",
" yijbar = p.getSolution(Y[i,j])\n",
" yjkbar = p.getSolution(Y[j,k])\n",
" yikbar = p.getSolution(Y[i,k])\n",
"\n",
" if yijbar + yjkbar + yikbar + 1 < xibar + xjbar + xkbar - eps:\n",
" n_violated += 1\n",
" ineq1_triples.append((i,j,k))\n",
"\n",
" if -yijbar + yjkbar - yikbar + xibar < - eps:\n",
" n_violated += 1\n",
" ineq2_triples.append((i,j,k))\n",
" \n",
" if -yijbar - yjkbar + yikbar + xjbar < - eps:\n",
" n_violated += 1\n",
" ineq3_triples.append((i,j,k))\n",
" \n",
" if yijbar - yjkbar - yikbar + xkbar < - eps:\n",
" n_violated += 1\n",
" ineq4_triples.append((i,j,k))\n",
"\n",
" if n_violated >= MAXINEQ:\n",
" break\n",
" \n",
" if n_violated > 0:\n",
" for i,j,k in ineq1_triples:\n",
" p.addConstraint(Y[i,j] + Y[j,k] + Y[i,k] + 1 >= x[i] + x[j] + x[k])\n",
" for i,j,k in ineq2_triples:\n",
" p.addConstraint(-Y[i,j] + Y[j,k] - Y[i,k] + x[i] >= 0)\n",
" for i,j,k in ineq3_triples:\n",
" p.addConstraint(-Y[i,j] - Y[j,k] + Y[i,k] + x[j] >= 0)\n",
" for i,j,k in ineq4_triples:\n",
" p.addConstraint(Y[i,j] - Y[j,k] - Y[i,k] + x[k] >= 0)\n",
"\n",
"print('Finished. Lower bound:', p.getObjVal())"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Nonlinear optimization (trigonometric, polynomial)\n",
"\n",
"**Homework**. Create a MIP with nonlinear constraint and solve it using the Xpress solver using the general method `p.solve()`. Minimize a variant of the well-known Rosenbrock's *banana* function\n",
"\n",
"$$\n",
"f(x,y) = (a-x)^2 + b(y-x^2)^2\n",
"$$\n",
"\n",
"where $a=\\frac{5}{2}$ and $b=3$, and where both $x$ and $y$ are integer.\n",
"\n",
"Can you find an optimum for $a=4$ without solving the problem?"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# STUDENT HOMEWORK\n",
"\n",
"import xpress as xp\n",
"\n",
"x = xp.var(vartype=xp.integer)\n",
"y = xp.var(vartype=xp.integer)\n",
"\n",
"a = 4\n",
"b = 3\n",
"\n",
"p = xp.problem(x,y,\n",
" (a-x)**2 + b*(y-x**2)**2)\n",
"\n",
"p.solve()\n",
"\n",
"# The optimum for a=4 is x = 4 and y = 16, for which the \n",
"# objective function is zero as both components of the\n",
"# objective are nonnegative."
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.6"
}
},
"nbformat": 4,
"nbformat_minor": 4
}