{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Exercise session 1: Overview of the Xpress Python Interface\n",
"\n",
"## Creating a simple LP\n",
"\n",
"First, create the simple linear program:\n",
"\n",
"$$\n",
"\\begin{array}{lll}\n",
" \\min & x + y\\\\\n",
" \\textrm{s.t.} & 2x + 3y \\ge 6\\\\\n",
" & 4x + 2y \\ge 7\\\\\n",
" & x,y \\ge 0\n",
"\\end{array}\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Location of manual and examples can be found by calling appropriate functions:\n",
"\n",
"import xpress\n",
"print(xpress.manual())\n",
"print(xpress.examples())"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import xpress as xp\n",
"\n",
"x = xp.var(name='myvar_1') # variables by default non-negative\n",
"y = xp.var(name='myvar_2')\n",
"\n",
"constr1 = 2*x + 3*y >= 6\n",
"constr2 = 4*x + 2*y >= 7\n",
"\n",
"obj = x + y\n",
"\n",
"p = xp.problem(x, y, # <-- variables (always first to be added)\n",
" constr1, constr2, # <-- then constraints and objective (in no order)\n",
" obj)\n",
"\n",
"# Uncomment the following if you want to turn off logging from Xpress\n",
"# p.controls.outputlog = 0\n",
"\n",
"p.solve()\n",
"\n",
"print(f'solution: x={p.getSolution(x)}, y={p.getSolution(y)}')\n",
"\n",
"# Same effect with this call that uses names instead of objects:\n",
"# print(f'solution: x={p.getSolution('myvar_1')}, y={p.getSolution('myvar_2')}')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Add the constraint $x + 2y \\le 2$ to the problem, then re-solve. The problem should now be infeasible."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# STUDENT EXERCISE\n",
"\n",
"bad_constr = x + 2*y <= 2\n",
"p.addConstraint(bad_constr)\n",
"p.solve()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Use the `iisall` function to understand what made the problem infeasible. If needed, type `help(p.iisall)` for its docstring."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"p.iisall()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"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",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"p.delConstraint(bad_constr)\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)}')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Knapsack problem\n",
"\n",
"Formulate and solve a knapsack problem with the following value/weight vectors:\n",
"\n",
"$$\n",
"\\begin{array}{lllrrrrrrr}\n",
"v& =& (12,&15,& 9,&11,& 8,& 7,&5)\\\\\n",
"w& =& (13, &18,& 9,& 12,& 8,& 10,& 4)\n",
"\\end{array}\n",
"$$\n",
"\n",
"and with knapsack capacity $C=40$. You will need the `xpress.Sum` operator. The formulation 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"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import xpress as xp\n",
"\n",
"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())"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Reading and solving a MIP\n",
"\n",
"Let's try to solve something bigger, like a file from the MIPLIB 2017 collection. Specifically, `roi2alpha3n4`. It should be available in the directory with this exercise session; load it and solve it as a Python problem."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"filename = 'mipinst.mps.gz'\n",
"p = xp.problem()\n",
"p.read(filename)\n",
"p.controls.maxtime = -30\n",
"p.solve()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Add a new constraint: the norm of the variable vector is at most 10. Even if the problem was read from a file, you can get the vector of all variables with `problem.getVariable`.\n",
"\n",
"Note that before the constraint is added we'll have to postsolve the problem (using the `problem.postsolve()` function) as, for reasons of performance, the problem remains in a presolved state after the `solve()` call."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"p.postsolve()\n",
"v = np.array(p.getVariable(), dtype=xp.npvar)\n",
"\n",
"p.addConstraint(xp.Dot(v,v) <= 100)\n",
"\n",
"# Or equivalently:\n",
"#\n",
"# v = p.getVariable()\n",
"# p.addConstraint(xp.Sum(i**2 for i in v) <= 100)\n",
"#\n",
"# Note that this form does not require importing numpy\n",
"\n",
"p.controls.maxtime=-30\n",
"p.solve()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"With a real MIP at hand, how about we try to see the structure of the coefficient structure from a graphical standpoint? Use `matplotlib` to draw a density chart."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"import scipy.sparse as sparse\n",
"\n",
"p.postsolve()\n",
"mstart, mclind, dmatval = [], [], []\n",
"p.getrows(mstart, mclind, dmatval, p.attributes.elems, 0, p.attributes.rows - 1)\n",
"\n",
"# Suppose the matrix is as follows:\n",
"#\n",
"# A = [0 0 2 9 0 7\n",
"# 1 0 0 5 1 9\n",
"# 0 0 8 5 0 3]\n",
"#\n",
"# Then we want a sparse representation where we store only information\n",
"# on the nonzero coefficients. This comes with three lists (I'm using\n",
"# spaces liberally to show information of each row):\n",
"#\n",
"# coeff = [2 9 7 1 5 1 9 8 5 3] <-- all nonzeros in the order they appear (first row, then column)\n",
"# indices = [2 3 5 0 3 4 5 2 3 5]\n",
"# mstart = [0 3 7 10]\n",
"#\n",
"# Here, mstart contains an entry for each row i: the position, within coeff and indices,\n",
"# of the first element of row i. This is obviously 0 for the first row. The difference\n",
"# between two consecutive elements of mstart is the number of corresponding nonzeros in a row.\n",
"# There are m+1 elements in mstart if the problem has m rows. The (m+1)-st contains the \n",
"# total number of elements in the problem.\n",
"\n",
"# mclind contains variable objects, not indices.\n",
"# Extract the index for all of them with getIndex()\n",
"mclindnum = [p.getIndex(v) for v in mclind]\n",
"A = sparse.csr_matrix ((dmatval, mclindnum, mstart))\n",
"\n",
"# plt.matshow(A.toarray())\n",
"plt.spy(A, markersize=.2)\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"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",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"\n",
"# The following two statements create two equivalent NumPy arrays of variables.\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')\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can also rewrite the knapsack problem using NumPy constructs."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"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)\n",
"\n",
"# xp.Dot only works with numpy arrays (of numbers and/or Xpress objects: expressions and variables)\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())"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Exercise: Create a random 20x30 numpy matrix $A$ and vectors ${\\bf b}$ and ${\\bf c}$, then create a problem with c as objective function coefficient vector and constraints $A{\\bf x} \\le {\\bf b}$."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"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()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's spice it up: make the objective function quadratic, with a matrix $Q$ that is PSD but not necessarily diagonal. Make the objective\n",
"\n",
"$$\n",
"(x-x_0)^T Q (x-x_0)\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# The following has a strong-enough diagonal that we're sure it is PSD (in fact, it's diagonally dominant)\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",
"# Minimize weighted (with Q) 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()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"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",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import xpress as xp\n",
"\n",
"n = 12 # number of periods\n",
"N = range(n)\n",
"\n",
"# Input parameters\n",
"\n",
"C = 2 # proportional cost\n",
"F = 9000 # fixed cost\n",
"L = 10 # min. production level\n",
"U = 90 # max. production level\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",
"\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 upper bound of M immediately\n",
"\n",
"# Total production cost (to be minimized)\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",
"\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])))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Homework**: Rewrite the problem using NumPy for defining parameters, variables, constraints, and objective."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# STUDENT EXERCISE\n",
"\n",
"import xpress as xp\n",
"import numpy as np\n",
"\n",
"# Demand\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') # Can assign upper bound of M immediately\n",
"\n",
"# Total production cost (to be minimized)\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",
"\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])))"
]
}
],
"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
}