import numpy as np5 Linear algebra and optimization
In this chapter we will see how Python can be used to perform linear algebra tasks and solve (integer) linear optimization problems.
Linear algebra forms the basis of many statistical tasks like principal component analysis (PCA) that can reduce the dimensionality of the data and may improve computational efficiency.
Applications of (integer) linear optimization range from scheduling, routing, transportation and food aid problems to solving Sudoku puzzles.
Let us first import NumPy again. Later on we will also use its linear algebra package (or module) linalg.
5.1 Linear algebra
In this section we will often refer to one-dimensional arrays as mathematical row or column vectors, and to two-dimensional arrays as matrices.
5.1.1 Matrix multiplications
Recall from Section 4.1.1 that for an m \times n matrix A and m \times 1 column vector x, the command A*x gives a matrix in which every column of A gets pointwise multiplied by the column vector x, that is,
\begin{align*} A*x = & \left[ \begin{matrix} a_{00} & a_{01} & \dots & a_{0(n-1)} \\ a_{10} & a_{11} & \dots & a_{1(n-1)} \\ \vdots & \vdots & \ddots & \vdots \\ a_{(m-1)0} & a_{(m-1)1} & \dots & a_{(m-1)(n-1)} \end{matrix} \right]*\left[ \begin{matrix} x_{0} \\ x_{1} \\ \vdots \\ x_{(m-1)} \\ \end{matrix} \right] \\ = & \left[ \begin{matrix} a_{00} \cdot x_0 & a_{01} \cdot x_0& \dots & a_{0(n-1)} \cdot x_0\\ a_{10} \cdot x_1 & a_{11} \cdot x_1& \dots & a_{1(n-1)} \cdot x_1\\ \vdots & \vdots& \ddots & \vdots \\ a_{(m-1)0} \cdot x_{m-1} & a_{(m-1)1} \cdot x_{m-1} & \dots & a_{(m-1)(n-1)} \cdot x_{m-1} \end{matrix} \right] \end{align*}
A = np.array([
[1,2],
[3,4]])
x = np.array([1,5])
x = x[:,None] # Turn x into (2,1) shaped column vector
print('A = \n', A)
print('x = \n', x)
print('A*x = \n', A*x)A =
[[1 2]
[3 4]]
x =
[[1]
[5]]
A*x =
[[ 1 2]
[15 20]]In the linear algebra sense, the multiplication Ax gives an m \times 1 column vector defined by
\begin{align*} Ax = & \left[ \begin{matrix} a_{00} & a_{01} & \dots & a_{0(n-1)} \\ a_{10} & a_{11} & \dots & a_{1(n-1)} \\ \vdots & \vdots & \ddots & \vdots\\ a_{(m-1)0} & a_{(m-1)1} & \dots & a_{(m-1)(n-1)} \end{matrix} \right]\left[ \begin{matrix} x_{0} \\ x_{1} \\ \vdots \\ x_{(n-1)} \\ \end{matrix} \right] \\ = & \left[ \begin{matrix} a_{00}x_0 + a_{01}x_1 + \dots + a_{0(n-1)}x_{n-1}\\ a_{10}x_0 + a_{11}x_1 + \dots + a_{1(n-1)}x_{n-1}\\ \vdots \\ a_{(m-1)0}x_0 + a_{(m-1)1}x_1 + \dots + a_{(m-1)(n-1)}x_{n-1}\\ \end{matrix} \right] \end{align*}
This can be achieved in Python with the @ (‘at’) operator, that is, by writing A @ x (the spaces are not needed; they are included here for readability).
print('A = \n', A)
print('x = \n', x)
print('Ax = \n', A @ x)A =
[[1 2]
[3 4]]
x =
[[1]
[5]]
Ax =
[[11]
[23]]If we talk about matrix multiplication, it will always be clear from context whether we mean A*x, which is element-wise multiplication with broadcasting, or A @ x, which is multiplication in the linear algebra sense.
We can also use @ to multiply a k \times m matrix and m \times n matrix with each other in the linear algebra sense.
A = np.array([
[1,2],
[3,4],
[3,7]])
B = np.array([
[2,5,1,7],
[3,4,-1,8]])
# Note that k = 3, m = 2, n = 4
print('AB = \n',A @ B)AB =
[[ 8 13 -1 23]
[18 31 -1 53]
[27 43 -4 77]]A special case of this is if we multiply a column vector
x = \left[\begin{matrix} x_0\\ x_1 \\ \vdots \\x_{n-1} \end{matrix} \right]
with its transpose (being a row vector) resulting in a matrix that contains x_ix_j at position (i,j).
x = np.array([
[1],
[2],
[4]])
print('xx^T = \n', x @ x.T)xx^T =
[[ 1 2 4]
[ 2 4 8]
[ 4 8 16]]5.1.2 Matrix properties
Using the linalg subpackage from NumPy, we can compute various well-known properties of a matrix (that you should recall from your Linear Algebra course).
np.linalg.matrix_rank(A): Computes the rank of matrix A.np.linalg.det(A): Computes the determinant of square A.np.linalg.eig(A): Computes eigenvalues and (right) eigenvectors of square matrix A, i.e., values \lambda that satisfy Av = \lambda v for some (eigen)vector v \neq 0.np.linalg.inv(A): Computes the inverse matrix A^{-1} of a square, invertible matrix A.np.linalg.norm(A): Computes the Frobenius norm of a matrix or a vector.
Finally, there is also the trace function np.trace(A) in NumPy that computes the trace of A, i.e., the sum of the elements on the diagonal. This function is implemented directly in NumPy, so you don’t have to go to the subpackage linalg first.
Let us look at some examples of these commands.
A = np.array([
[1,2],
[3,4]])# Rank of A
print('Rank of A: ', np.linalg.matrix_rank(A))Rank of A: 2# Determinant of A
print('Determinant of A: ', np.linalg.det(A))Determinant of A: -2.0000000000000004# Trace of A
print('Trace of A: ', np.trace(A))Trace of A: 5# Inverse of A
print('Inverse of A: \n', np.linalg.inv(A))Inverse of A:
[[-2. 1. ]
[ 1.5 -0.5]]The function linalg.eig() is special in that it does not output one, but two arrays.
The first array is one-dimensional and contains the eigenvalues [\lambda_0,\dots,\lambda_{n-1}] of A; the second array is two-dimensional, and contains the corresponding eigenvectors v^0,\dots,v^{n-1} as columns of the matrix. The eigenvalues and eigenvectors are paired in the sense that Av^i = \lambda_iv^i for i = 0,\dots,n-1.
# Eigenvalues and eigenvectors of A
lambdas, V = np.linalg.eig(A)
print('Eigenvalues of A: \n', lambdas)
print('Matrix with eigenvectors of A: \n', V)Eigenvalues of A:
[-0.37228132 5.37228132]
Matrix with eigenvectors of A:
[[-0.82456484 -0.41597356]
[ 0.56576746 -0.90937671]]You can store only the eigenvalues or the eigenvectors in a variable by writing _ in the position of the other output.
# Only eigenvalues of A
lambdas, _ = np.linalg.eig(A)
# Only eigenvectors of A
_, V = np.linalg.eig(A)If you want to recall how you can handle multiple outputs, and do not assign the ones that your are not interested in, have a look at Appendix B.
Let us check that the eigenvalues and eigenvectors indeed satisfy Av^i = \lambda_iv^i for i = 0,1. Although the output should be zero, due to numerical inaccuracies the output can also be a number very close to zero, typically in the order of 1e-16.
v0 = V[:,0]
v1 = V[:,1]
print('Verify first eigenvector:', A @ v0 - lambdas[0]*v0)
print('Verify second eigenvector:', A @ v1 - lambdas[1]*v1)Verify first eigenvector: [0.00000000e+00 5.55111512e-17]
Verify second eigenvector: [0. 0.]The function linalg.norm(A) computes the Frobenius norm of a matrix A = (a_{ij})_{i = 0,\dots,m-1, j = 1,\dots,n-1} \in \mathbb{R}^{m \times n}, that is defined as
||A||_{\text{F}} = \sqrt{\sum_{i=0}^{m-1} \sum_{j = 0}^{n-1} a_{ij}^2 }.
If you apply the function to a row or column vector x = [x_0,\dots,x_{n-1}], you get the usual L^2-norm defined as
||x||_{2} = \sqrt{\sum_{i=0}^{n-1} x_i^2 }.
You can also use this function in a vectorized manner. For example, if you want to compute the L^2-norm of every row in a matrix, you can use the axis keyword argument.
A = np.array([
[1,2],
[3,4],
[5,6]])
row_norms = np.linalg.norm(A,axis=1)
print(row_norms)[2.23606798 5. 7.81024968]5.1.3 Equation solving
For a given m \times n matrix A = (a_{ij}) and m \times 1 column vector b = [b_0,\dots,b_{m-1}]^T, perhaps the most important question in linear algebra is to compute an x = [x_0,\dots,x_{n-1}]^T that satisfies Ax = b, i.e., the linear systems of equations \begin{array}{rrrrcrrrl} a_{00}x_0 &+& a_{01}x_1 &+& \dots &+& a_{0(n-1)}x_{n-1} &=& b_0 \\ a_{10}x_0 &+& a_{11}x_1 &+& \dots &+& a_{1(n-1)}x_{n-1} &=& b_1 \\ &&&& \vdots &&&& \\ a_{(m-1)0}x_0 &+& a_{(m-1)1}x_1 &+& \dots &+& a_{(m-1)(n-1)}x_{n-1} &=& b_{m-1} \end{array}.
If A is a square, invertible matrix then this can be done with linalg.solve(A,b). This gives the unique solution x = A^{-1}b.
A = np.array([
[1,2],
[3,4]])
b = np.array([[1],[4]])
print('Solution to Ax = b: \n', np.linalg.solve(A,b))Solution to Ax = b:
[[ 2. ]
[-0.5]]In fact, you can compute x directly as x = A^{-1}b using the inverse of A that we have seen earlier. For computational efficiency, it is recommended, though, to use np.linalg.solve(A,b) instead to solve Ax = b.
print('Solution to Ax = b: \n', np.linalg.inv(A) @ b)Solution to Ax = b:
[[ 2. ]
[-0.5]]If m > n, the system is overdetermined and may have no (unique) solution. This can happen when there are more constraints than variables, i.e., when m > n. In this case, we can find an approximately optimal x using the least squares method. It finds an x that solves the (non-linear) problem
\min_x ||Ax - b||_2^2.
Note that ||Ax - b||_2^2 = 0 if and only if Ax = b. The least squares method can be executed with linalg.lstsq(A,b); we will see this method in the next chapter, when we consider non-linear optimization problems. This function can also be used to solve underdetermined systems, where n > m.
5.2 Linear optimization
In this section we will see how to solve small- and larger-scale optimization problems using Gurobi which is a state-of-the-art software program for solving optimization problems. As a student, you can get a free academic license to use this software.
Below you will find the instructions to install Gurobi if you want to use it on your personal device; the university computers are equipped with Gurobi already and the program should be ready to use on them.
5.2.1 Basics
All the relevant functionality needed to use Gurobi via Python is contained in the gurobipy package. For our purposes, we only need two classes from this package:
- The
Modelclass that we will use to build optimization problems - The
GRBclass that contains various “constants” that we use to include words such as maximize and minimize, as well as binary and integer.
In Python you can include the modules by adding the following line at the top of a script.
from gurobipy import Model, GRBExample: Duplo problem
To illustrate how to implement optimization problems in Python, we consider the Duplo toy problem that you have seen in the course Linear Optimization (35B108).
In this illustrative problem, we have 8 small and 6 big (Duplo) bricks available to build two products: chairs and tables. A chair requires 2 small and 1 large brick to build, and a table 2 small and 2 large bricks.
Building one chair yields a profit of 15 euros, whereas building a table yields a profit of 20 euros. The goal is to decide the number of chairs and tables to build so that our total profit is maximized.
Below we have given a linear optimization formulation that solves the Duplo problem based on the given input data (the numbers) and the (for now possibly fractional) decision variables
- x_1: Number of chairs to build, and
- x_2: Number of tables to build.
\begin{array}{llll} \text{maximize}& 15x_1 + 20x_2 && \text{(profit)} \\ \text{subject to} &x_1 + 2x_2 \leq 6 &&\text{(large bricks)} \\ &2x_1 + 2x_2 \leq 8 &&\text{(small bricks)} \\ &x_1, x_2 \geq 0 \end{array}
Convince yourself that this formulation is indeed correct based on the problem description.
This problem can be implemented in Python and solved with Gurobi using the Python code below that will be explained step-by-step.
from gurobipy import Model, GRB
# Initialize Gurobi model.
duplo = Model(name='Duplo problem')
# Declare the two decision variables.
x1 = duplo.addVar(name='chairs')
x2 = duplo.addVar(name='tables')
# Specify the objective function.
duplo.setObjective(15*x1 + 20*x2, sense=GRB.MAXIMIZE)
# Add the resource constraints on bricks.
duplo.addConstr(x1 + 2*x2 <= 6, name='big-bricks')
duplo.addConstr(2*x1 + 2*x2 <= 8, name='small-bricks')
# Optimize the model
duplo.optimize()
# Inspect the objective value
print('Objective value:', duplo.ObjVal)
# Inspect the optimal solution (x1,x2)
for var in duplo.getVars():
print(var.VarName, "=", var.X)The implementation of the linear optimization problem follows the concept of object-oriented programming. We first create an object from the Model() class after which we use various methods to add the decision variables (which are themselves also objects), objective function and constraints to the object. Finally, we use a method to optimize the problem.
After optimizing, we can inspect the solution by accessing various attributes of the Model() object and decision variable objects.
Model object
To start, we will create a variable (or object) duplo that will contain all the information of the problem, such as the decision variables, objective function and constraints.
In programming terms, we create an object from the Model() class. We can give the instance a name by adding the name keyword argument, with value 'Duplo problem' in our case. We add this name with quotations so that Python knows this is plain text.
# Initialize Gurobi model.
duplo = Model(name='Duplo problem')You can think of the Model class object duplo as a large box that we are going to fill with decision variables, an objective function, and constraints. We will also reserve a small space in this box to store information about the optimal solution to the linear optimization problem, once we have optimized the problem.
The concept of having “objects” is in fact what the programming language Python is centered around. This is known as object-oriented programming, which you will learn about more later in your first full programming course. You can have a look, for example, here already if you are interested in this concept.
Decision variables
To add decision variables, an objective function, and constraints, to our Model object we use so-called methods which are Python functions.
To create a decision variable, we can use the addVar() method. As an input argument, we include a name for the decision variables.
# Declare the two decision variables.
x1 = duplo.addVar(name='chairs')
x2 = duplo.addVar(name='tables')If your Model object is not called
duplo, but for examplemodel_Duplo, then you should usemodel_Duplo.addVar(name='chairs')instead ofduplo.addVar(name='chairs'). The same applies to all later methods that are introduced. Never call a modelModel(with capital M) because this spelling is reserved for theModel()class in Gurobi.
When creating a variable, you can specify more properties. For example, we might have a lower and upper bound for a decision variable (such as a nonnegativity constraint). You can use the keyword arguments lb and ub, respectively, for this.
To see all these properties, you can check out the addVar documentation by typing help(Model.addVar) in Python. Looking up the documentation in this way can be done for every method.
help(Model.addVar)Help on cython_function_or_method in module gurobipy._model:
addVar(self, lb=0.0, ub=1e+100, obj=0.0, vtype='C', name='', column=None)
ROUTINE:
addVar(lb, ub, obj, vtype, name, column)
PURPOSE:
Add a variable to the model.
ARGUMENTS:
lb (float): Lower bound (default is zero)
ub (float): Upper bound (default is infinite)
obj (float): Objective coefficient (default is zero)
vtype (string): Variable type (default is GRB.CONTINUOUS)
name (string): Variable name (default is no name)
column (Column): Initial coefficients for column (default is None)
RETURN VALUE:
The created Var object.
EXAMPLE:
v = model.addVar(ub=2.0, name="NewVar")
If we wanted to create a variable named test variable, having lower bound 3 and upper bound 7, we can do this with x = duplo.addVar(lb=3,ub=7,name='test variable') using in particular the lb and ub keyword arguments.
As you can see above, all the keyword arguments have default values, meaning that if we do not specify them, Gurobi uses the specified default value for them. For example, by default a variable is continuous (vtype keyword argument) and non-negative (lb keyword argument), so we do not have to specify the nonnegativity constraint x_1,x_2 \geq 0 ourselves for the Duplo example.
Also note that ub is set to ifinite which roughly speaking indicates that the decision variable has no upper bound value by default.
NoteInteger and binary linear optimization problems
In practice, decision variables are often restricted to be binary, i.e., in the domain \{0,1\}, or integer-valued, i.e., in the domain \mathbb{Z}= \{\dots,-3,-2,-1,0,1,2,3,\dots\}, as opposed to being continuous, i.e., with domain \mathbb{R}.
You can force these conditions using the vtype keyword argument. As an example, if we wanted to model the number of chairs to be binary or integer, we could set, respectively
x1 = duplo.addVar(name='chairs', vtype=GRB.BINARY)x1 = duplo.addVar(name='chairs', vtype=GRB.INTEGER)
Objective function
Now that we have our decision variables, we can use them to define the objective function and the constraints.
We use the setObjective() method to do so by referring to the variables created above. The sense keyword argument specifies whether we want to maximize or minimize the objective function. For this, we use GRB.MAXIMIZE or GRB.MINIMIZE, respectively.
# Specify the objective function.
duplo.setObjective(15*x1 + 20*x2, sense=GRB.MAXIMIZE)Constraints
The next step is to declare the constraints using the addConstr() method. Also here, we can specify the constraints referring to the variables x1 and x2. We give the constraints a name as well.
# Add the resource constraints on bricks.
duplo.addConstr(x1 + 2*x2 <= 6, name='big-bricks')
duplo.addConstr(2*x1 + 2*x2 <= 8, name='small-bricks')<gurobi.Constr *Awaiting Model Update*>
NoteOther constraint types
If you instead want to model a \geq or = constraint, you can do this with >= or ==, respectively, instead of <=.
duplo.addConstr(x1 + 2*x2 >= 6, name='big-bricks')duplo.addConstr(x1 + 2*x2 == 6, name='big-bricks')
Now that the model is completely specified, we are ready to compute the optimal solution. We can do this using the optimize() method applied to our model duplo.
# Optimize the model
duplo.optimize()Gurobi Optimizer version 13.0.1 build v13.0.1rc0 (win64 - Windows 11+.0 (26200.2))
CPU model: 12th Gen Intel(R) Core(TM) i7-1265U, instruction set [SSE2|AVX|AVX2]
Thread count: 10 physical cores, 12 logical processors, using up to 12 threads
Optimize a model with 2 rows, 2 columns and 4 nonzeros (Max)
Model fingerprint: 0xadc88607
Model has 2 linear objective coefficients
Coefficient statistics:
Matrix range [1e+00, 2e+00]
Objective range [2e+01, 2e+01]
Bounds range [0e+00, 0e+00]
RHS range [6e+00, 8e+00]
Presolve time: 0.01s
Presolved: 2 rows, 2 columns, 4 nonzeros
Iteration Objective Primal Inf. Dual Inf. Time
0 3.5000000e+31 3.500000e+30 3.500000e+01 0s
2 7.0000000e+01 0.000000e+00 0.000000e+00 0s
Solved in 2 iterations and 0.01 seconds (0.00 work units)
Optimal objective 7.000000000e+01We can see some output from Gurobi, and the last line tells us that Gurobi found an optimal solution with objective value 70.
Recalling our linear optimization problem object as being a “box” filled with decision variables, an objective function and constraints, the box now also contains a small space where the optimal solution is stored.
Let us have a look at how to access information about the optimal solution. The objective function value of the optimal solution, that was also displayed in the output when optimizating the model with duplo.optimize() is stored in the ObjVal attribute. An attribute is a piece of information about an object in Python.
# Print text 'Objective value' and the duplo.ObjVal attribute value
print('Objective value:', duplo.ObjVal)Objective value: 70.0The optimal value of a decision variable can be obtained by accessing the X attribute of a variable.
print('x1 = ', x1.X)
print('x2 = ', x2.X)x1 = 2.0
x2 = 2.0If the model has many variables, printing variables in this way is a cumbersome approach. It is then easier to iterate over all variables of the model using the getVars() method. This method creates a list of all the decision variables of the model.
print(duplo.getVars())[<gurobi.Var chairs (value 2.0)>, <gurobi.Var tables (value 2.0)>]As you can see above, we already see the optimal values of the decision variables, but also a lot of unnecessary information.
Let us print these values in a nicer format by iterating over the list duplo.getVars() using a for-loop with loop variable var; you can choose a name other than var for the loop variable if you want. We print for every variable var its name, stored in the attribute VarName, and its value, stored in the attribute X. In between, we print the = symbol in plain text (hence, the quotation marks).
Recall that if you want to print multiple data types (such as variable values and text) in a print() statement, you should separate them with commas.
for var in duplo.getVars():
print(var.VarName, "=", var.X)chairs = 2.0
tables = 2.05.2.2 Beyond the basics
Writing out the objective function and constraints explicitly makes it easy to check if the problem is correctly implemented, but this is not scalable. In practice, problems can have thousands of decision variables and constraints, in which case a more efficient implementation is needed.
To illustrate useful commands for implementing large-scale optimization problems, we will develop general code for implementing a standard-form problem with m constraints and n decision variables: \min \mathbf{c'x} subject to \mathbf{Ax} = \mathbf{b}, \mathbf{x} \geq \mathbf{0} where \mathbf{c} \in \mathbb{R}^n, \mathbf{b} \in \mathbb{R}^m, \mathbf{A} \in \mathbb{R}^{m \times n}, and x = [x_1,\dots,x_{n}].
More explicitly, the problem is as follows:
\begin{array}{rrrcrrlll} \min & c_1 x_1 &+& \cdots& + & c_{n}x_{n}& & & \\ {\rm s.t.} & a_{11}x_1& +& \cdots& +& a_{1n} x_{n} &= &b_1 &\\ & & & \vdots & & & & \\ & a_{m1}x_1& +& \cdots& +& a_{mn} x_{n} &= &b_{m}& \\ & x_1, & & \cdots& & x_{n} & \geq & 0 & \\ \end{array}
Suppose we are given the following input data in NumPy arrays.
import numpy as np
# Objective function coefficients
c = np.array([3, 1, 7, 2, -1]) # n = 5
# Matrix A (m = 3, n = 5)
A = np.array([
[2, 1, 0, 3, -1],
[1, 0, 4, 0, 2],
[0, 3, -2, 1, 1]
])
# Right-hand side vector b
b = np.array([10, 12, 5]) # m = 3Below we have given an implementation that solves this standard-form problem. The implementation is explained below.
from gurobipy import Model, GRB, quicksum
# We define m and n here for convenience later
m = len(b)
n = len(c)
# Create Model object
standard_form = Model("Standard form problem")
# Add decision variables
x = standard_form.addVars(n,name=[f"Var{i+1}" for i in range(n)])
# Create objective function
standard_form.setObjective(quicksum(x[j]*c[j] for j in range(n)), sense=GRB.MINIMIZE)
# Create constraints
for i in range(m):
standard_form.addConstr(quicksum(A[i][j]*x[j] for j in range(n)) == b[i])
# Optimize the model
standard_form.optimize()
# Print optimal values of decision variables
for var in standard_form.getVars():
print(var.VarName, "=", var.X)Decision variables
Instead of creating decision variables one-by-one with addVar(), it is also possible to instantiate a whole collection of variables using its “plural” version addVars().
# We define m and n here for convenience later
m = len(b) # m = 3 for the example data
n = len(c) # n = 5 for the example data
# Create Model object
standard_form = Model("Standard form problem")
# Add decision variables
x =standard_form.addVars(n, name=["Var1","Var2","Var3","Var4","Var5"]) With the syntax x = standard_form.addVars(n) we tell Python to create n (= 5) decision variables. You can use the i-th decision variable to define constraints and the objective function by indexing x at position i-1, i.e., by using x[i-1]. Recall that Python starts counting from zero when indexing data objects such as lists.
The name keyword argument is optional, but if used, the number of elements in the list should match the value of n. In the code above we hard-coded the list for n = 5. This can be avoided using list comprehension as illustrated below.
# Add decision variables
x = standard_form.addVars(n, name=[f"Var{i+1}" for i in range(n)])Objective function
We continue by creating the objective function by indexing the variables in x and the array c containing the objective function coefficients. If we write out the objective function completely, it would look like this.
# Create objective function
standard_form.setObjective(x[0]*c[0] + x[1]*c[1] + x[2]*c[2] + x[3]*c[3] + x[4]*c[4])We can implement it more compactly by using the quicksum() function from gurobipy.
# Create objective function
standard_form.setObjective(quicksum(x[j]*c[j] for j in range(n))) The above command intuitively carries out the following steps:
- Create the list
[x[0]*c[0], x[1]*c[1], ..., x[n-1]*c[n-1]]with list comprehension. - Use the
quicksum()function to add up the elements in this list. - Set the resulting quantity as the objective function.
Constraints
We continue by adding the equality constraints. We use a for-loop to loop over the equalities a_{i1}x_1 + \dots + a_{in}x_n = b_i for i = 1,\dots,m (note that Python starts counting at 0, though).
# Create constraints
for i in range(m):
standard_form.addConstr(quicksum(A[i][j]*x[j] for j in range(n)) == b[i])For a fixed i, the list [A[i][j]*x[j] for j in range(n)] contains the coefficients of the i+1-th row in A multiplied by the variables in x, i.e., [a_{i1}x_1, ..., a_{in}x_n]. Applying quicksum() to this list then computes the expression a_{i1}x_1 + \dots + a_{in}x_n. This should equal b_i so we model the constraint by quicksum([A[i][j]*x[j] for j in range(n)]) == b[i].
Having created the decision variables, objective and constraints, we can then optimize the problem and inspect the values of the optimal decision variables, just like before. Recall that the constraints x_i \geq 0 do not have to be added explicitly, because when creating the decision variables, their default range is [0,\infty).
Putting both input data and implementation together, we obtain the following code.
import numpy as np
from gurobipy import Model, GRB, quicksum
# Objective function coefficients
c = np.array([3, 1, 7, 2, -1]) # n = 5
# Matrix A (m = 3, n = 5)
A = np.array([
[2, 1, 0, 3, -1],
[1, 0, 4, 0, 2],
[0, 3, -2, 1, 1]
])
# Right-hand side vector b
b = np.array([10, 12, 5]) # m = 3
# We define m and n here for convenience later
m = 3
n = 5
# Create Model object
standard_form = Model("Standard form problem")
# Suppress Gurobi output message (you do not need to know this command)
standard_form.setParam('OutputFlag',0)
# Add decision variables
x = standard_form.addVars(n,name=[f"Var{i+1}" for i in range(n)])
# Create objective function
standard_form.setObjective(quicksum(x[j]*c[j] for j in range(n)), sense=GRB.MINIMIZE)
# Create constraints
for i in range(m):
standard_form.addConstr(quicksum(A[i][j]*x[j] for j in range(n)) == b[i])
# Optimize the model
standard_form.optimize()
# Print optimal values of decision variables
for var in standard_form.getVars():
print(var.VarName, "=", var.X)Var1 = 4.75
Var2 = 0.0
Var3 = 0.0
Var4 = 1.3749999999999998
Var5 = 3.6255.2.3 Double-indexed decision variables
In some cases, linear optimization problems have decision variables with multiple indices, for example, x_{ij} for i = 1,\dots,k and j = 1, \dots, \ell. This means we would have k \cdot \ell decision variables in total. We will see such an example in the exercises.
In this case, the decision variables can be defined in a two-dimensional object.
k = 3
ell = 5
# Create Model object
double = Model("Double-indexed model")
# Add decision variables
x = double.addVars(k, ell) You can also add a name keyword argument to addVars(). Instead of giving a custom name to every x_{ij}, you can also enter one word (string) as the value of name.
x = double.addVars(k, ell, name="combinations") When printed, the variables in this case will be displayed as combinations[0,0], combinations[0,1], ..., combinations[k-1,ell-1].
To use variable x_{ij} in a constraint or the objective function, you can access it through x[i-1,j-1]. These variables can also be used in combination with list comprehension-like approach and quicksum().
For example, if the objective function is
\sum_{i = 1,\dots,k} \sum_{j=1,\dots,\ell} c_{ij}x_{ij}
for given coefficients c_{ij}, then we can set it in our Gurobi model using a nested list comprehension-like approach.
# Objective function coefficient array c
c = np.array([
[2, 1, 0, 3, -1],
[1, 0, 4, 0, 2],
[0, 3, -2, 1, 1]
])
# Objective function expression
z = quicksum(quicksum(c[i,j]*x[i,j] for j in range(ell)) for i in range(k))
# Set objective in model
double.setObjective(z, sense=GRB.MINIMIZE)To see how the double application of quicksum() works, note that for a fixed i, the expression y_i = quicksum(c[i,j]*x[i,j] for j in range(ell)) computes the inner summation \sum_{j=1,\dots,\ell} c_{ij}x_{ij}. Next summing over i then means quicksum(y_i for i in range(k)) which then corresponds to the double summation that is the objective function.
It is also possible to set constraints in this way, either with single or double list comprehension. You will get to practice with this in the exercises.
5.2.4 Remarks
Python has various other packages and functions for implementing optimization problems.
The linprog function of SciPy is typically convenient for small linear optimization problems, especially for problems whose input data (such as the constraint matrices) you can define explicitly.
The pulp package is more useful when the input data is defined implicitly, e.g., using different sets of inequalities defined through indices. Furthermore, for larger-scale problems, as you will typically encounter in practice, it is better to use state-of-the-art solvers such as Gurobi, CPLEX, or MOSEK. PuLP is able to be coupled to such solvers, that is, you can define a problem in PuLP and then have it solved with the external solver software.
Many of these solvers can also handle other problems such as mixed-integer linear optimization problems, in which there is a mixture of continuous, integer and binary variables, and various nonlinear optimization problems. We will see more of those in the next chapter.