Gurobi 8はPython3.5で利用可能。なので、Python3.5でどうやればGurobi 8を利用して最適化問題を解けるのかを記述する。自分の使用した環境はmacOS+Gurobi 8.0.1なのでGUROBI_HOME環境変数は/Library/gurobi801/mac64となっている。

インストール

まずはインストール。 Gurobiが環境にインストールされていることを前提。

セットアップはいとも簡単。

$ python3.5 setup.py install --user

これだけ。あとはライセンスファイルgurobi.licを~/.local/share/gurobiに設置して、GRB_LICENSE_FILE環境変数を用意するだけ。

$ GUROBI_DATA=~/.local/share/gurobi
$ mkdir -p $GUROBI_DATA
$ export GRB_LICENSE_FILE=$GUROBI_DATA

永続的にGUROBI_DATAにライセンスファイルを設置しておくなら、~/.bash_profileとか~/.zshrcとかに

export GRB_LICENSE_FILE=~/.local/share/gurobi

とか書いておくと良い。

モデルを書く

さて、準備は整ったので、モデルを書いていく。なお、今回使用する問題は次のようなSOCPとする。

$\displaystyle \min\left\{\langle \boldsymbol{c}, \boldsymbol{x} \rangle : a_i \leq x_i \leq b_i \enspace (i =0, \cdots, n), x_0 \geq \sqrt{\sum_{j = 1}^n (d_j x_j)^2} \right\}$

ここで、 $\langle \cdot, \cdot \rangle$ は内積を表していることに注意。

まず、gurobipyモジュールをインポートしておかなければ話にならない。

import gurobipy as grb

ここではgurobipyをgrbという名前で利用することにしている。そして、モデルを用意。

model = grb.Model('SOCP')

ここの`'SOCP'は何でもよいが、わかりやすさのために、モデル名とかにしておくのが無難。

次に、変数を用意していく。

x = [
    model.addVar(
        vtype=grb.GRB.CONTINUOUS,
        name="x[{}]".format(i + 1),
        lb=a[i],
        ub=b[i],
    ) \
    for i in range(n + 1)
]

lbやubはそれぞれlower boundとupper boundの略。また、vtypeは変数のタイプを表している。他は表参考。

扱いたい領域	利用するべき定数
連続値	`grb.GRB.CONTINUOUS`
整数値	`grb.GRB.INTEGER`
バイナリー値	`grb.GRB.BINARY`

変数を用意したら、モデルへ反映。

model.update()

次に目的関数を設定*1。

model.setObjective(
    grb.quicksum(c[i] * x[i] for i in range(n + 1)),
    grb.GRB.MINIMIZE
)

ここでは最小化のためにgrb.GRB.MINIMIZE定数を渡しているが、最大化したければ、grb.GRB.MAXIMIZEを渡せば良い。そしてモデルの更新。

model.update()

最後に制約を追加する。 math.sqrtはそのまま使えないので、制約式を

$\displaystyle x_0^2 \geq \sum_{j = 1}^n (d_j x_j)^2$

と変形して定義する。

model.addConstr(
    x[0] * x[0] >= grb.quicksum((d[i - 1]**2) * x[i] * x[i] for i in range(1, n + 1)),
    "q"
)

第2引数は制約式ごとの名前。忘れずに更新。

model.update()

あとは解くだけで良い。

model.optimize()

解けたかどうかはstatusで確認できる。

model.status

grb.GRB.Status.OPTIMALであれば最適値が求まっており、grb.GRB.Status.SUBOPTIMALであれば、誤差が大きすぎて最適値かどうかはわからないが、とりあえずは解が求まっている。

求まった決定変数の値はxで、変数名はvarNameで確認できる。

print("{} = {}".format(x[0].varvName, x[0].x)

求まった目的関数値はobjValで取得できる。

model.objVal

これらをクラスの中に押しこめたコードを記しておく。

import gurobipy as grb

class Socp(object):
    def __init__(self, a, b, c):
        self.__model = grb.Model('SOCP')

        self.__x = [
            self.__model.addVar(
                vtype=grb.GRB.CONTINUOUS,
                name="x[{}]".format(i + 1),
                lb=a[i],
                ub=b[i]
            ) \
             for i in range(n + 1)
        ]
        self.__model.update()

        self.__model.setObjective(
            grb.quicksum(c[i] * self.__x[i] for i in range(n + 1)),
            grb.GRB.MINIMIZE
        )
        self.__model.update()

        self.__model.addConstr(
            x[0] * x[0] >= grb.quicksum((d[i - 1]**2) * self.__x[i] * self.__x[i] for i in range(1, n + 1)),
            "q"
        )
        self.__model.update()

    def optimize(self):
        self.__model.optimize()
        if self.__model.status in [grb.GRB.Status.SUBOPTIMAL, grb.GRB.Status.OPTIMAL]:
            return self.__model.objVal

    @property
    def x(self, i):
        return self.__x[i].x

Tips

もし、出力を止めたければ

model.Params.outputFlag = 0

とすることでGurobiからの余計な出力を止めることができる。

あと、そのまま解こうとすると

Q matrix is not positive semi-definite (PSD)

と怒られることがある。具体的には

model.addConstr(
    v[i] * (t[i + 1] - t[i] - tau[i]) >= d[i],
    "q{i + 1}".format(i)
)

というコードで怒られた。どれも連続値で、tauとdは実数値の定数だった。これは2次錐制約が暗に変数の値が正となることを仮定しているため。なので、

y = [
    model.addVar(
        vtype=grb.GRB.CONTINOUS,
        name="y[{}]".fortmat(i)
    ) \
    for i in range(m)
]

というスラック変数を用意して、

model.addConstr(
    y[i] == t[i + 1] - t[i] - tau[i],
    "q'{}".format(i + 1)
)

という制約式を足すことで回避できる。

*1:quiksumを利用する際にはジェネレータ式を利用すると便利

バージョン管理された人

subversionで管理されてます

Gurobi+Python3.5で楽に最適化問題を解く

インストール

モデルを書く

Tips