JunPyo Park

UNIST FE Lab

Contact me

2023 © JunPyo Park

최적화 문제 풀이 Solver GUROBI 사용하기

2021, May 08    

최적화(Optimization) 문제 풀이를 위한 GUROBI solver 사용하기

Written by JunPyo Park

• 최적화(Optimization) 문제를 풀기 위한 솔버 중 하나인 GUROBI 사용법을 간단하게 정리해 보았습니다.
• 사용 후기: 규모가 큰 문제를 풀려면 라이센스가 필요합니다. 학교 다니시면 아카데믹 라이센스 무료로 이용가능합니다.

Reference

• GUROBI Python Documentation
• GUROBI Tutorial Slide
• Introduction to Mathematical Optimization Modeling

---

Table of Contents

• 1  Installation
• 2  Simple MIP(MIxed Integer Linear Programmin) Example
  • 2.1  Objective (최적 문제의 목적식)
  • 2.2  Constraints (제약식)
  • 2.3  전체 코드
  • 2.4  모델(최적화 문제) 생성
  • 2.5  모델에 변수 추가 하기
  • 2.6  목적식(objective) 세팅
  • 2.7  제약식(constraints) 추가
  • 2.8  모델 최적화
  • 2.9  결과 분석
• 3  Consulting Company Problem
  • 3.1  Problem Description
  • 3.2  Data
  • 3.3  Decision Variables
  • 3.4  Job constraints
  • 3.5  Resource constraints
  • 3.6  Objective function

Installation

pip install gurobipy

Simple MIP(MIxed Integer Linear Programmin) Example

Objective (최적 문제의 목적식)

$\text{Maximize}$
$x + y + 2 z$

Constraints (제약식)

$\text{subject to}$
$x + 2 y + 3 z <= 4$
$x + y >= 1$
$x, y, z$ $\text{binary}$

x,y,z 가 0,1 값을 가지는 binary 형태이고 목적식과 제약식이 모두 linear 하기에 위 문제는 MIP라고 불립니다.

전체 코드

밑에 부분별 설명이 있습니다.

import gurobipy as gp
from gurobipy import GRB

try:

    # Create a new model
    m = gp.Model("mip1")

    # Create variables
    x = m.addVar(vtype=GRB.BINARY, name="x")
    y = m.addVar(vtype=GRB.BINARY, name="y")
    z = m.addVar(vtype=GRB.BINARY, name="z")

    # Set objective
    m.setObjective(x + y + 2 * z, GRB.MAXIMIZE)

    # Add constraint: x + 2 y + 3 z <= 4
    # m.addConstr(x + 2 * y + 3 * z <= 4, "c0")
    constr = gp.LinExpr()
    constr += x
    constr += 2 * y
    constr += 3 * z
    constr
    m.addConstr(constr <= 4, "c0")

    # Add constraint: x + y >= 1
    m.addConstr(x + y >= 1, "c1")

    # Optimize model
    m.optimize()

    for v in m.getVars():
        print('%s %g' % (v.varName, v.x))

    print('Obj: %g' % m.objVal)

except gp.GurobiError as e:
    print('Error code ' + str(e.errno) + ': ' + str(e))

except AttributeError:
    print('Encountered an attribute error')

Gurobi Optimizer version 9.1.2 build v9.1.2rc0 (win64)
Thread count: 8 physical cores, 16 logical processors, using up to 16 threads
Optimize a model with 2 rows, 3 columns and 5 nonzeros
Model fingerprint: 0x98886187
Variable types: 0 continuous, 3 integer (3 binary)
Coefficient statistics:
  Matrix range     [1e+00, 3e+00]
  Objective range  [1e+00, 2e+00]
  Bounds range     [1e+00, 1e+00]
  RHS range        [1e+00, 4e+00]
Found heuristic solution: objective 2.0000000
Presolve removed 2 rows and 3 columns
Presolve time: 0.00s
Presolve: All rows and columns removed

Explored 0 nodes (0 simplex iterations) in 0.00 seconds
Thread count was 1 (of 16 available processors)

Solution count 2: 3 2 

Optimal solution found (tolerance 1.00e-04)
Best objective 3.000000000000e+00, best bound 3.000000000000e+00, gap 0.0000%
x 1
y 0
z 1
Obj: 3

모델(최적화 문제) 생성

# Create a new model
m = gp.Model("mip1")

모델에 하나의 최적화 문제 배정. 모델은 변수들과 constraint들 그리고 최적화 문제와 관련된 attribute(변수 범위, 변수 타입, 목적식의 계수 등)들로 구성

모델에 변수 추가 하기

# Create binary variable
x = m.addVar(vtype=GRB.BINARY, name="x")

model.addVar() 메서드를 통해 모델에 변수를 추가 가능합니다. 어레이 같은 구조를 통해 여러 변수를 한번에 추가하려면 addVars() 메서드를 사용하면 됩니다.

help(m.addVar)

Help on method addVar in module gurobipy:

addVar(lb=0.0, ub=1e+100, obj=0.0, vtype='C', name='', column=None) method of gurobipy.Model instance
    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")

목적식(objective) 세팅

Maximize
$x + y + 2 z$

# Set objective
m.setObjective(x + y + 2 * z, GRB.MAXIMIZE)

식이 복잡할 경우 아래와 같이 단계적으로(incrementally) 수식 표현 가능

obj = gp.LinExpr()
obj += x
obj += y
obj += 2*z
m.setObjective(obj, GRB.MAXIMIZE)

제약식(constraints) 추가

# Add constraint: x + 2 y + 3 z <= 4
m.addConstr(x + 2 * y + 3 * z <= 4, "c0")

• model.addConstr() 메서드를 통해 제약식 추가 가능
• 두번째 arg("c0")는 제약식의 이름
• 식이 복잡할 경우 위의 2.3과 같이 incrementally 추가 가능
• addConstrs() 를 사용해서 여러 제약식 한번에 추가 가능

모델 최적화

# Optimize model
m.optimize()

추가된 목적식과 제약식을 바탕으로 최적화를 진행합니다.

결과 분석

for v in m.getVars():
    print('%s %g' % (v.varName, v.x))

• variable.varName() 메서드를 통해 사전에 정의한 변수 이름에 접근 가능
• v.x: 변수 v의 최적 value

print('Obj: %g' % m.objVal) # optimal value

---

Consulting Company Problem

이거 보면서 풀어보면 감 잡힙니다. 번역은 귀찮아서....

Problem Description

Consider a consulting company that has three open positions: Tester, Java Developer, and Architect. The three top candidates (resources) for the positions are: Carlos, Joe, and Monika. The consulting company administered competency tests to each candidate in order to assess their ability to perform each of the jobs. The results of these tests are called matching scores. Assume that only one candidate can be assigned to a job, and at most one job can be assigned to a candidate.

The problem is to determine an assignment of resources and jobs such that each job is fulfilled, each resource is assigned to at most one job, and the total matching scores of the assignments is maximized.

Data

The list $R$ contains the names of the three resources: Carlos, Joe, and Monika.

The list $J$ contains the names of the job positions: Tester, Java Developer, and Architect.

$r \in R$: index and set of resources. The resource $r$ belongs to the set of resources $R$.

$j \in J$: index and set of jobs. The job $j$ belongs to the set of jobs $J$.

For each resource $r$ and job $j$, there is a corresponding matching score $s$. The matching score $s$ can only take values between 0 and 100. That is, $s_{r,j} \in [0, 100]$ for all resources $r \in R$ and jobs $j \in J$.

# Resource and job sets
R = ['Carlos', 'Joe', 'Monika']
J = ['Tester', 'JavaDeveloper', 'Architect']

help(gp.multidict)

Help on built-in function multidict in module gurobipy:

multidict(...)
    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] } )

# Matching score data
combinations, scores = gp.multidict({
    ('Carlos', 'Tester'): 53,
    ('Carlos', 'JavaDeveloper'): 27,
    ('Carlos', 'Architect'): 13,
    ('Joe', 'Tester'): 80,
    ('Joe', 'JavaDeveloper'): 47,
    ('Joe', 'Architect'): 67,
    ('Monika', 'Tester'): 53,
    ('Monika', 'JavaDeveloper'): 73,
    ('Monika', 'Architect'): 47
})

The following constructor creates an empty Model object “m”. We specify the model name by passing the string "RAP" as an argument. The Model object “m” holds a single optimization problem. It consists of a set of variables, a set of constraints, and the objective function.

# Declare and initialize model
m = gp.Model('RAP')

Decision Variables

To solve this assignment problem, we need to identify which resource is assigned to which job. We introduce a decision variable for each possible assignment of resources to jobs. Therefore, we have 9 decision variables.

To simplify the mathematical notation of the model formulation, we define the following indices for resources and jobs:

For example, $x_{2,1}$ is the decision variable associated with assigning the resource Joe to the job Tester. Therefore, decision variable $x_{r,j}$ equals 1 if resource $r \in R$ is assigned to job $j \in J$, and 0 otherwise.

The Model.addVars() method creates the decision variables for a Model object. This method returns a Gurobi tupledict object that contains the newly created variables. We supply the combinations object as the first argument to specify the variable indices. The name keyword is used to specify a name for the newly created decision variables. By default, variables are assumed to be non-negative.

combinations

<gurobi.tuplelist (9 tuples, 2 values each):
 ( Carlos , Tester        )
 ( Carlos , JavaDeveloper )
 ( Carlos , Architect     )
 ( Joe    , Tester        )
 ( Joe    , JavaDeveloper )
 ( Joe    , Architect     )
 ( Monika , Tester        )
 ( Monika , JavaDeveloper )
 ( Monika , Architect     )
>

# Create decision variables for the RAP model
x = m.addVars(combinations, name="assign")

Job constraints

We now discuss the constraints associated with the jobs. These constraints need to ensure that each job is filled by exactly one resource.

The job constraint for the Tester position requires that resource 1 (Carlos), resource 2 (Joe), or resource 3 (Monika) is assigned to this job. This corresponds to the following constraint.

Constraint (Tester=1)

$$ x_{1,1} + x_{2,1} + x_{3,1} = 1 $$

Similarly, the constraints for the Java Developer and Architect positions can be defined as follows.

Constraint (Java Developer = 2)

$$ x_{1,2} + x_{2,2} + x_{3,2} = 1 $$

Constraint (Architect = 3)

$$ x_{1,3} + x_{2,3} + x_{3,3} = 1 $$

The job constraints are defined by the columns of the following table.

In general, the constraint for the job Tester can defined as follows.

$$ x_{1,1} + x_{2,1} + x_{3,1} = \sum_{r=1}^{3 } x_{r,1} = \sum_{r \in R} x_{r,1} = 1 $$

All of the job constraints can be defined in a similarly succinct manner. For each job $j \in J$, take the summation of the decision variables over all the resources. We can write the corresponding job constraint as follows.

$$ \sum_{r \in R} x_{r,j} = 1 $$

The Model.addConstrs() method of the Gurobi/Python API defines the job constraints of the Model object “m”. This method returns a Gurobi tupledict object that contains the job constraints. The first argument of this method, "x.sum(‘*’, j)", is the sum method and defines the LHS of the jobs constraints as follows: For each job $j$ in the set of jobs $J$, take the summation of the decision variables over all the resources. The $==$ defines an equality constraint, and the number "1" is the RHS of the constraints. These constraints are saying that exactly one resource should be assigned to each job. The second argument is the name of this type of constraints.

# Create job constraints
jobs = m.addConstrs((x.sum('*',j) == 1 for j in J), name='job')

Resource constraints

The constraints for the resources need to ensure that at most one job is assigned to each resource. That is, it is possible that not all the resources are assigned.

For example, we want a constraint that requires Carlos to be assigned to at most one of the jobs: either job 1 (Tester), job 2 (Java Developer ), or job 3 (Architect). We can write this constraint as follows.

Constraint (Carlos=1)

$$ x_{1, 1} + x_{1, 2} + x_{1, 3} \leq 1. $$

This constraint is less or equal than 1 to allow the possibility that Carlos is not assigned to any job. Similarly, the constraints for the resources Joe and Monika can be defined as follows:

Constraint (Joe=2)

$$ x_{2, 1} + x_{2, 2} + x_{2, 3} \leq 1. $$

Constraint (Monika=3)

$$ x_{3, 1} + x_{3, 2} + x_{3, 3} \leq 1. $$

Observe that the resource constraints are defined by the rows of the following table.

The constraint for the resource Carlos can be defined as follows.

$$ x_{1, 1} + x_{1, 2} + x_{1, 3} = \sum_{j=1}^{3 } x_{1,j} = \sum_{j \in J} x_{1,j} \leq 1. $$

Again, each of these constraints can be written in a succinct manner. For each resource $r \in R$, take the summation of the decision variables over all the jobs. We can write the corresponding resource constraint as follows.

$$ \sum_{j \in J} x_{r,j} \leq 1. $$

The Model.addConstrs() method of the Gurobi/Python API defines the resource constraints of the Model object “m”. The first argument of this method, "x.sum(r, ‘*’)", is the sum method and defines the LHS of the resource constraints as follows: For each resource $r$ in the set of resources $R$, take the summation of the decision variables over all the jobs. The $<=$ defines a less or equal constraints, and the number “1” is the RHS of the constraints. These constraints are saying that each resource can be assigned to at most 1 job. The second argument is the name of this type of constraints.

# Create resource constraints
resources = m.addConstrs((x.sum(r,'*') <= 1 for r in R), name='resource')

Objective function

The objective function is to maximize the total matching score of the assignments that satisfy the job and resource constraints.

For the Tester job, the matching score is $53x_{1,1}$, if resource Carlos is assigned, or $80x_{2,1}$, if resource Joe is assigned, or $53x_{3,1}$, if resource Monika is assigned. Consequently, the matching score for the Tester job is as follows, where only one term in this summation will be nonzero.

$$ 53x_{1,1} + 80x_{2,1} + 53x_{3,1}. $$

Similarly, the matching scores for the Java Developer and Architect jobs are defined as follows. The matching score for the Java Developer job is:

$$ 27x_{1, 2} + 47x_{2, 2} + 73x_{3, 2}. $$

The matching score for the Architect job is:

$$ 13x_{1, 3} + 67x_{2, 3} + 47x_{3, 3}. $$

The total matching score is the summation of each cell in the following table.

The goal is to maximize the total matching score of the assignments. Therefore, the objective function is defined as follows.

\begin{equation} \text{Maximize} \quad (53x_{1,1} + 80x_{2,1} + 53x_{3,1}) \; + \end{equation}\begin{equation} \quad (27x_{1, 2} + 47x_{2, 2} + 73x_{3, 2}) \; + \end{equation}\begin{equation} \quad (13x_{1, 3} + 67x_{2, 3} + 47x_{3, 3}). \end{equation}

Each term in parenthesis in the objective function can be expressed as follows.

\begin{equation} (53x_{1,1} + 80x_{2,1} + 53x_{3,1}) = \sum_{r \in R} s_{r,1}x_{r,1}. \end{equation}\begin{equation} (27x_{1, 2} + 47x_{2, 2} + 73x_{3, 2}) = \sum_{r \in R} s_{r,2}x_{r,2}. \end{equation}\begin{equation} (13x_{1, 3} + 67x_{2, 3} + 47x_{3, 3}) = \sum_{r \in R} s_{r,3}x_{r,3}. \end{equation}

Hence, the objective function can be concisely written as:

\begin{equation} \text{Maximize} \quad \sum_{j \in J} \sum_{r \in R} s_{r,j}x_{r,j}. \end{equation}

The Model.setObjective() method of the Gurobi/Python API defines the objective function of the Model object “m”. The objective expression is specified in the first argument of this method. Notice that both the matching score parameters “score” and the assignment decision variables “x” are defined over the “combinations” keys. Therefore, we use the method “x.prod(score)” to obtain the summation of the elementwise multiplication of the "score" matrix and the "x" variable matrix. The second argument, GRB.MAXIMIZE, is the optimization "sense." In this case, we want to maximize the total matching scores of all assignments.

scores

{('Carlos', 'Tester'): 53,
 ('Carlos', 'JavaDeveloper'): 27,
 ('Carlos', 'Architect'): 13,
 ('Joe', 'Tester'): 80,
 ('Joe', 'JavaDeveloper'): 47,
 ('Joe', 'Architect'): 67,
 ('Monika', 'Tester'): 53,
 ('Monika', 'JavaDeveloper'): 73,
 ('Monika', 'Architect'): 47}

x.prod(scores)

<gurobi.LinExpr: 53.0 assign[Carlos,Tester] + 27.0 assign[Carlos,JavaDeveloper] + 13.0 assign[Carlos,Architect] + 80.0 assign[Joe,Tester] + 47.0 assign[Joe,JavaDeveloper] + 67.0 assign[Joe,Architect] + 53.0 assign[Monika,Tester] + 73.0 assign[Monika,JavaDeveloper] + 47.0 assign[Monika,Architect]>

# Objective: maximize total matching score of all assignments
m.setObjective(x.prod(scores), GRB.MAXIMIZE)

# Save model for inspection
m.write('RAP.lp')

m.display()

Maximize
   <gurobi.LinExpr: 53.0 assign[Carlos,Tester] + 27.0 assign[Carlos,JavaDeveloper] + 13.0 assign[Carlos,Architect] + 80.0 assign[Joe,Tester] + 47.0 assign[Joe,JavaDeveloper] + 67.0 assign[Joe,Architect] + 53.0 assign[Monika,Tester] + 73.0 assign[Monika,JavaDeveloper] + 47.0 assign[Monika,Architect]>
Subject To
   job[Tester] : <gurobi.LinExpr: assign[Carlos,Tester] + assign[Joe,Tester] + assign[Monika,Tester]> = 1.0
   job[JavaDeveloper] : <gurobi.LinExpr: assign[Carlos,JavaDeveloper] + assign[Joe,JavaDeveloper] + assign[Monika,JavaDeveloper]> = 1.0
   job[Architect] : <gurobi.LinExpr: assign[Carlos,Architect] + assign[Joe,Architect] + assign[Monika,Architect]> = 1.0
   resource[Carlos] : <gurobi.LinExpr: assign[Carlos,Tester] + assign[Carlos,JavaDeveloper] + assign[Carlos,Architect]> <= 1.0
   resource[Joe] : <gurobi.LinExpr: assign[Joe,Tester] + assign[Joe,JavaDeveloper] + assign[Joe,Architect]> <= 1.0
   resource[Monika] : <gurobi.LinExpr: assign[Monika,Tester] + assign[Monika,JavaDeveloper] + assign[Monika,Architect]> <= 1.0

# Run optimization engine
m.optimize()

Gurobi Optimizer version 9.1.2 build v9.1.2rc0 (win64)
Thread count: 8 physical cores, 16 logical processors, using up to 16 threads
Optimize a model with 6 rows, 9 columns and 18 nonzeros
Model fingerprint: 0xb343b6eb
Coefficient statistics:
  Matrix range     [1e+00, 1e+00]
  Objective range  [1e+01, 8e+01]
  Bounds range     [0e+00, 0e+00]
  RHS range        [1e+00, 1e+00]
Presolve time: 0.00s
Presolved: 6 rows, 9 columns, 18 nonzeros

Iteration    Objective       Primal Inf.    Dual Inf.      Time
       0    4.6000000e+32   1.800000e+31   4.600000e+02      0s
       5    1.9300000e+02   0.000000e+00   0.000000e+00      0s

Solved in 5 iterations and 0.01 seconds
Optimal objective  1.930000000e+02

# Display optimal values of decision variables
for v in m.getVars():
    if v.x > 1e-6:
        print(v.varName, v.x)

# Display optimal total matching score
print('Total matching score: ', m.objVal)

assign[Carlos,Tester] 1.0
assign[Joe,Architect] 1.0
assign[Monika,JavaDeveloper] 1.0
Total matching score:  193.0

Twitter Facebook Google+

# Optimization # MIP # Mixed Integer Programming # Linear Programming # Solver