GurobiとPythonで数理最適化

数理最適化ソルバーGurobiにはPython含む複数の言語のインターフェースが用意されていて、
gurobipyモジュールを使うと使い慣れたPythonがソルバーのモデリング言語になる。

例題

http://www.fujilab.dnj.ynu.ac.jp/lecture/system2.pdf

あるレストランで，手持ちの材料からハンバーグとオムレツを作って利益を最大にしたいと考
えている．手持ちの材料は，
• ひき肉 3800 [g]
• タマネギ 2100 [g]
• ケチャップ 1200 [g]
であり，それぞれの品を作るのに必要な材料の量は，
• ハンバーグ 1個あたり，ひき肉 60 [g]，タマネギ 20 [g]，ケチャップ 20 [g]
• オムレツ1個あたり，　ひき肉 40 [g]，タマネギ 30 [g]，ケチャップ 10 [g]
であるとする．(他に必要な材料は十分な量があるものとする)
販売価格は，
• ハンバーグ 400 [円/個]
• オムレツ 300 [円/個]
とする．総売上を最大にするには，それぞれハンバーグとオムレツを幾つずつ作れば良いか?

定式化

x1をハンバーグの個数、x2をオムレツの個数として、次のように定式化する。

$\begin{eqnarray} \mbox{maximize}&~& 400x_1 + 300x_2 \end{eqnarray}$
$\begin{eqnarray} \mbox{subject to} &~& 60x_1 + 40x_2 \leq 3800 \\ &~& 20x_1 + 30x_2 \leq 2100 \\ &~& 20x_1 + 10x_2 \leq 1200 \\ &~& x_i \geq 0, \\ &~& x_i \in {\bf\mbox{Z}},\:\:\: i=1,2 \end{eqnarray}$

Pythonによるモデルの記述

gurobipyモジュールのメソッドを使ってモデルに変数や目的関数、制約式などを追加していく。

from gurobipy import *

try:
    # Create a new model
    m = Model("restaurant")
    
    # Create variables
    x = [m.addVar(vtype=GRB.INTEGER, name='x'+str(i+1)) for i in range(2)]
    
    # Integrate new variables
    m.update()
    
    # Set objective
    m.setObjective(400*x[0] + 300*x[1], GRB.MAXIMIZE)
    
    # Add constraints
    m.addConstr(60*x[0] + 40*x[1] <= 3800, "ground meat")
    m.addConstr(20*x[0] + 30*x[1] <= 2100, "onion")
    m.addConstr(20*x[0] + 10*x[1] <= 1200, "ketchup")

    m.optimize()

    # OUTPUT
    for v in m.getVars():
        print v.varName, v.x
    print 'Obj:', m.objVal

except GurobiError:
    print 'Error reported'

コード内に非負制約を明示していないが、それは決定変数を宣言するaddVarメソッドのデフォルト値として設定されている。引数の数やデフォルト値は対話シェルでも確認することができる。

>>> help(m.addVar)
Help on method addVar in module gurobipy:

addVar(...) 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")

つまり、v=m.addVar() は v=m.addVar(0,inf,0,GRB.CONTINUOUS,'',None) と同義になる。
また、キーワード引数を上手く使うと楽に書ける。

決定変数の添字が定式化では1オリジンなのにソースコードでは（リストを使っているので）0オリジンで気持ち悪いという場合は辞書を使うと良い。

>>> x = {}
>>> for i in range(1, 3):
...     x[i] = m.addVar(vtype =GRB.INTEGER, name='x'+str(i))

>>> m.update()
>>> x
{1: <gurobi.Var x1>, 2: <gurobi.Var x2>}

もちろん、辞書のキーには整数だけでなく文字列やタプルなどを使うこともできる。

実行結果

Optimize a model with 3 rows, 2 columns and 6 nonzeros
Found heuristic solution: objective 24000
Presolve time: 0.00s
Presolved: 3 rows, 2 columns, 6 nonzeros
Variable types: 0 continuous, 2 integer (0 binary)

Root relaxation: objective 2.700000e+04, 3 iterations, 0.00 seconds

    Nodes    |    Current Node    |     Objective Bounds      |     Work
 Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

/*    0     0               0    27000.000000 27000.0000  0.00%     -    0s

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

Optimal solution found (tolerance 1.00e-04)
Best objective 2.700000000000e+04, best bound 2.700000000000e+04, gap 0.0%
x1 30.0
x2 50.0
Obj: 27000.0

よって、ハンバーグ30個、オムレツ50個のとき、総売上27000円で最大となる。