Implementation in Python

Implementation in Python#

In this Notebook, we show how you can use Gurobi with Python API to solve the clay and sand extraction problem.

Note

Gurobi cannot be loaded in this online book, so download this notebook to work on it with Gurobi locally installed. Instruction on how to do that are in the README.md and PA_2_4_A_gurobilicious.ipynb of week 2.4.

Problem description#

  • A company extracts sand and clay from a site which when sold gives a profit of 57 and 60 monetary units per thousand units of product. For this extraction, sand needs a manpower of 50 men x hours to extract 1000 units of product 1 whilst clay needs 13 men x hours for 1000 units.

  • 4h of backhoe work are needed to extract one thousand units of sand and 5h of backhoe work per thousand units of clay. The number of hours needed of truck transport is 8h and 4h respectively for sand and clay for each 1000 units of product transported.

  • The company has a work schedule of 40 hours per week for the men but also for the equipment (truck and backhoe). There are 5 men who can be used interchangeably between the transport of the two products. There is only one truck and one backhoe.

  • What should be the extraction plan of this company in a week? What should their objective be? What decisions need to be made?

Gurobi tutorial#

  • In this notebook, we are going to build a simple model using gurobipy.

  • You can turn to the built-in function help() to see how different commands work, for example:

help(model.addVar())

Import gurobipy package#

After preparing the environment for gurobi+python, you should be able to run the following code to import the gurobipy package.

import gurobipy as gp

Construct the model#

Initialization#

  • First, initialize an empty model.

model = gp.Model()

Define attributes#

  • Now let’s store some parameter values we need for constructing the model.

  • The multydict() function in gurobipy allows us to store these values conveniently.

gurobipy.multidict()


help(gp.multidict)

Help on cython_function_or_method in module gurobipy:

multidict(data)
    ROUTINE:
      multidict(data)
    PURPOSE:
      Split a single dictionary into multiple dictionaries.
    ARGUMENTS:
      data: A dictionary that maps each key to a list of 'n' values.
    RETURN VALUE:
      A list of the shared keys, followed by individual tupledicts.
    EXAMPLE:
      (keys, dict1, dict2) = multidict( {
               'key1': [1, 2],
               'key2': [1, 3],
               'key3': [1, 4] } )

Profit

Manpower

Backhoe work

Truck transport

Sand

57

50

4

8

Clay

60

13

5

4

 
# define attributes
material, profit, manpower, backhoe_work, truck_trans = gp.multidict({
    'sand' :[57, 50, 4, 8], 'clay' :[60, 13, 5, 4]
})
men_equipment, work_hour = gp.multidict({
    'men' :[200], 'truck' :[40], 'backhoe' :[40]
})

Add desicion variables#

  • This step is to add desicion variables. You can use these two functions to add decision variables to the model.

Model.addVar() # add a single variable
Model.addVars() # add multiple variables

You can define the type of the variables. Such as:

GRB.CONTINUOUS  # continuous variable
GRB.INTEGER # integer variable

  • Use the built-in function in the following cells to see detailed information about Model.addVar() and variable types in class GRB

help(gp.Model.addVar)

Help on cython_function_or_method in module gurobipy:

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")

  • Use following cell to add decision variables x into the model. \(x_0\) and \(x_1\) represent the amount (ton/week) of sand and clay that will be extracted (notice that the index starts from 0).

  • Update the model after adding variables.

x = model.addVars(2, vtype = gp.GRB.CONTINUOUS, name = 'x')  # x0, x1
model.update() #Remember to update the model after adding decision variables

Add constraints#

You can use the following functions to add constrains to the model.

Model.addConstr()  # add a single constrain
Model.addConstrs()  # add multiple constrin

Besides, for convenience, you can use a function gurobipy.quicksum() to calculate sums and build expressions.

help(gp.quicksum)

Help on cython_function_or_method in module gurobipy:

quicksum(list)
    ROUTINE:
      quicksum(list)
    
    PURPOSE:
      A quicker version of the Python built-in 'sum' function for building
      Gurobi expressions.
    
    ARGUMENTS:
      list: A list of terms.
    
    RETURN VALUE:
      An expression that represents the sum of the input arguments.
    
    EXAMPLE:
      expr = quicksum([x, y, z])
      expr = quicksum([1.0, 2*y, 3*z*z])

Turn to built-in function help() to see detailed information like following.

help(gp.Model.addConstr)

Help on cython_function_or_method in module gurobipy:

addConstr(self, lhs, sense=None, rhs=None, name='')
    ROUTINE:
      addConstr(tc, name)
    
    PURPOSE:
      Add a constraint to the model.
    
    ARGUMENTS:
      tc (TempConstr): The constraint to add
      name (string): Constraint name (default is no name)
    
    RETURN VALUE:
      Depending on the data of 'tc':
        - A Constr object if tc arose from a linear expression
        - A QConstr object if tc arose from a quadratic expression
        - An MConstr object if tc arose from a linear matrix expression
        - A GenConstr object if tc arose form a general constraint expression
    
    EXAMPLE:
      c = model.addConstr(x + y <= 1)
      c = model.addConstr(x + y + z == [1, 2])
      c = model.addConstr(x*x + y*y <= 1)
      c = model.addConstr(z == and_(y, x, w))
      c = model.addConstr(z == min_(x, y))
      c = model.addConstr((w == 1) >> (x + y <= 1))
      c = model.addConstr(A @ x <= b)

In this step, you need to add the following constraints.

\[\begin{split} \begin{align} & 8 x_0 + 4 x_1 \leq 40 \\ & 4 x_0 + 5 x_1 \leq 40 \\ & 50 x_0 + 13 x_1 \leq 200 \\ \end{align} \end{split}\]

Note that non-negativity constraints are not always necessary in Gurobi since by default, variables are non-negative in Gurobi.

model.addConstr(gp.quicksum(manpower[material[i]] * x[i] for i in range(2)) <= work_hour['men'], name = 'manpower limit')
model.addConstr(gp.quicksum(backhoe_work[material[i]] * x[i] for i in range(2)) <= work_hour['backhoe'], name = 'backhoe limit')
model.addConstr(gp.quicksum(truck_trans[material[i]] * x[i] for i in range(2)) <= work_hour['truck'], name = 'truck limit')
model.update()

Set objective function#

  • You can use the following command to set the objective

Model.setObjective()

  • Use built-in function help() to see detailed information

help(gp.Model.setObjective)

Help on cython_function_or_method in module gurobipy:

setObjective(self, expr, sense=None)
    ROUTINE:
      setObjective(expr, sense=None)
    
    PURPOSE:
      Set the model objective equal to a LinExpr or QuadExpr
    
    ARGUMENTS:
      expr: The desired objective function.  The objective can be
            a linear expression (LinExpr) or a quadratic expression (QuadExpr).
            This routine will replace the 'Obj' attribute on model variables
            with the corresponding values from the supplied expression.
      sense (optional): Objective sense.  By default, this method uses the
            modelSense model attribute to determine the sense.  Use
            GRB.MINIMIZE or GRB.MAXIMIZE to ignore modelSense and pick a
            specific sense.
    
    RETURN VALUE:
      None.
    
    EXAMPLE:
      model.setObjective(x + y)
      model.setObjective(x + y + 2*z, GRB.MAXIMIZE)

This step is to set the objective function shown in the following:

\( Max \quad L = 57 x_0 + 60 x_1 \)

model.setObjective(gp.quicksum(profit[material[i]] * x[i] for i in range(2)), gp.GRB.MAXIMIZE)

Solve#

Use following code to solve the model, and you can also show the optimization process log and limit the calculation time.

You can also print the objective function value and optimal decision variables values.

You can choose the algorithm used, or let Gurobi decide which algorithm is the best. See the following links for detailed information:

https://www.gurobi.com/documentation/9.5/refman/method.html

model.Params.LogToConsole = True  # show the optimization process log
model.Params.TimeLimit = 100  # limit the calculation time to 100s
model.Params.MultiObjMethod = 0  # default is -1 (let Gurobi decide)
model.optimize()

Set parameter TimeLimit to value 100
Set parameter MultiObjMethod to value 0
Gurobi Optimizer version 11.0.0 build v11.0.0rc2 (win64 - Windows 10.0 (19045.2))

CPU model: 11th Gen Intel(R) Core(TM) i7-1185G7 @ 3.00GHz, instruction set [SSE2|AVX|AVX2|AVX512]
Thread count: 4 physical cores, 8 logical processors, using up to 8 threads

Optimize a model with 3 rows, 2 columns and 6 nonzeros
Model fingerprint: 0xd942f567
Coefficient statistics:
  Matrix range     [4e+00, 5e+01]
  Objective range  [6e+01, 6e+01]
  Bounds range     [0e+00, 0e+00]
  RHS range        [4e+01, 2e+02]
Presolve time: 0.00s
Presolved: 3 rows, 2 columns, 6 nonzeros

Iteration    Objective       Primal Inf.    Dual Inf.      Time
       0    1.4625000e+31   3.609375e+30   1.462500e+01      0s
       2    4.9500000e+02   0.000000e+00   0.000000e+00      0s

Solved in 2 iterations and 0.01 seconds (0.00 work units)
Optimal objective  4.950000000e+02

print("Optimal Objective function Value", model.objVal)

Optimal Objective function Value 495.0

for var in model.getVars():
    print(f"{var.varName}: {round(var.X, 3)}")  # print the optimal decision variable values.

x[0]: 1.667
x[1]: 6.667

# see the model
model.display()

Maximize
  57.0 x[0] + 60.0 x[1]
Subject To
  manpower limit: 50.0 x[0] + 13.0 x[1] <= 200
  backhoe limit: 4.0 x[0] + 5.0 x[1] <= 40
  truck limit: 8.0 x[0] + 4.0 x[1] <= 40

C:\Users\tomvanwoudenbe\AppData\Local\Temp\ipykernel_8432\17538946.py:2: DeprecationWarning: Model.display() is deprecated
  model.display()

Exercise: The sand and clay extraction problem#

Now that you have seen how to model and solve a simple problem in gurobi, let’s try to make some changes in the problem. Imagine in the meantime the company has purchased another truck and has hired an extra (full time) worker to increase its profit. So you have to model and solve a very similar problem with minor changes to reflect the new employee and the new truck.

You need to make changes when defining attributes.

Solution

The correct change is:

men_equipment, work_hour = gurobipy.multidict({'men' :[240], 'truck' :[80], 'backhoe' :[40]})

Which leads to an Optimal Objective Value of 510.91 with x[0] equal to 3.434 and x[1] equal to 5.253.

Attribution

This chapter is written by Gonçalo Homem de Almeida Correia, Maria Nogal Macho, Jie Gao and Bahman Ahmadi. Find out more here.

Contents