サブロウ丸

主にプログラミングと数学

【MacOS】Pulp ソルバー選択 / 並列計算

Pulpについて

Pulp線形計画問題を解く Python パッケージです 1 .

install

pip install pulp

使い方

PuLPによるモデル作成方法 2 に分かりやすくまとめられています.

pulpの処理フロー f:id:inarizuuuushi:20190228183552p:plain

ソルバー選択

現在の環境で使用できるソルバーはpulp.pulpTestAll()で確認できます[^3].

>>> import pulp
>>> pulp.pulpTestAll()
* Solver <class 'pulp.solvers.PULP_CBC_CMD'> passed.
Solver <class 'pulp.solvers.CPLEX_PY'> unavailable
* Solver <class 'pulp.solvers.COIN_CMD'> passed.
Solver <class 'pulp.solvers.COINMP_DLL'> unavailable
Solver <class 'pulp.solvers.GLPK_CMD'> unavailable
...

単にPuLPをインストールしただけだと、使えるソルバーはCBCだけです。 別途、GLPK(無料)、Gurobi(有料)、CPLEX(有料)などをインストールすると、以下のようにソルバーを変更できます [^3] 。

from pulp import LpProblem, PulpSolverError
solver1 = pulp.GLPK_CMD()
solver2 = pulp.GUROBI_CMD()
solver3 = pulp.CPLEX_CMD()
for solver in [solver1, solver2, solver3]:
    try:
        m = LpProblem()
        m.solve(solver)
    except PulpSolverError:
        print(solver.path, 'はインストールされていません')




CBC

pulpのデフォルトソルバーで, solve()を実行するとCBCソルバーで問題が解かれます.

オプション(並列計算など)

また, PULP_CBC_CMDを使用すれば以下のようなオプションを使用することができます.

PULP_CBC_CMD(path=None, keepFiles=0, mip=1, msg=0,
cuts=None, presolve=None, dual=None, strong=None,
options=[], fracGap=None, maxSeconds=None, threads=None)

例えば4スレッドで並列計算させたい場合は

m = LpProblem()
solver = pulp.PULP_CBC_CMD(threads=4)
m.solve(solver)

のようにオプションを指定します.

  • keepFilesを1(True)に設定すると, 定義した問題が記録されたmpsファイルとその解.solファイルが削除されずに残ります.
  • 解きたい問題に整数変数が含まれる場合はmipをTrueに設定します. 整数変数が含まれない場合はFalseに設定して解くことができます.

  • cuts (Default: True)

    • cut generator を行なって, 探索をより効率的に行います.
  • preserve (Default: True)

    • モデルを分析して、冗長な制約, いくつかの変数を固定する制約、範囲に変換できる制約などを見つけます. 問題を最初に解決するには、効果がないことがわかっていない限り、これを行う価値があります.
  • strong (Default: 5)

    • 強分岐で調べる変数の数を決定します. どの変数に分岐するかを決めるために、コードはこの数までの満たされていない変数を選び、最小の上下の分岐を試みます.
  • fracGapは計算終了の判定となる, 解の上界と下界のGapのしきい値を指定できます.

  • maxSecondsは最大計算時間を指定できます. (単位は秒)

  • threadsはソルバー実行の並列数を指定できます.

  • optionsではそれ以外のCBCオプションを設定できます.

    • 例えばmaxsolオプションを使用する場合は options = ['maxsol 1']とします.

      • maxsolオプションは, 暫定解が何個見つかった時に計算を終了するかを指定できます. 最適解でなくても実行可能解が1つ欲しい場合はmaxsol 1とすればより速い時間で求めることができます.

cuts, presolve, strongは特に変更する必要はないでしょう. (より早く問題を解くためにチューニングする場合を除いて)

CBCで使用できるオプションについてはGAMS3 が参考になります. (上記の説明もこのサイトを参考にしました)



自分でインストールしたCBCソルバーを使用する

pulpと一緒にインストールされるCBCはバージョンが古い場合があり, 最新版のCBCソルバを用いて問題を解きたい場合は以下の処理を行う必要があります.

1. 公式サイト4からCBCをインストールする

maxOSにhomebrewを使ってCBCをインストールする場合は

brew tap coin-or-tools/coinor
brew install cbc --with-parallel --with-suite-sparse

を実行します. (--with-parallelをつけることで並列計算が可能になります. --with-suite-sparseはおまけです)

特にインストール先を変更していない場合は /usr/local/Cellar/配下にバージョンごとにディレクトリが作成されます.5







solver = COIN_CMD()で使用されるCBCソルバは(下記「ソルバ指定」節参照), PuLPの設定ファイル(pulp.cfg.*)のCbcPathで指定されたソルバが使用されます (デフォルトではターミナル上でcbcと打ったときに使用されるソルバ). なので, 方法としては以下の2通りあります.

2-1. 自分でインストールしたcbcソルバにパスを通す
2-2. 設定ファイルを書き換える

なので, 2-1か2-2のどちらかを行えば十分です.

2-1. cbcソルバにパスを通す

maxOSの環境に合わせて話を進めます. /usr/local/binにパスを通せば, cbcコマンドを使えるようになります.

オペレーティングシステム(OS)によって考え方は異なっているが、一般に/usr/local以下には、自分でコンパイルしたものを置いておくために使われる。 /usr/binや/binはシステムの標準のコマンドが置かれるが、/usr/local/binにはそれ以外のものを置く。[^1]

brewでは特にインストール先を変更していない場合は /usr/local/Cellar/配下にバージョンごとにディレクトリが作成されます.[^2] なので

ln -s /usr/local/Cellar/cbc/2.9.9_1(自分がインストールしたバージョン)/bin/cbc /usr/local/bin/cbc

を実行して, パスを張ります. これでターミナル上でcbcを実行すれば, インストールしたcbcが使用されます.

$ cbc
Welcome to the CBC MILP Solver 
Version: 2.9.9 
Build Date: Mar  1 2019 

CoinSolver takes input from arguments ( - switches to stdin)
Enter ? for list of commands or help
Coin:

これで, COIN_CMDにより, このソルバが使用されます.

2-2. PuLPの設定ファイルを書き換える

pulp.__path__によって確認できるPulPモジュールの保存先には, pulp.cfg.*ファイルがあり, maxOSの場合はpulp.cfg.osxに使用するソルバのパスなどの情報が書かれています. そのcbcに関するパスを下のように書き換えます.

# CbcPath = cbc
CbcPath = /usr/local/Cellar/cbc/2.9.9_1(自分がインストールしたバージョン)/bin/cbc

こうすれば, COIN_CMDによりリンク先のCBCソルバーが使用されることになります.




ソルバ指定

from pulp import LpProblem
m = LpProblem()

# COIN CBCソルバ
solver = pulp.COIN_CMD()
m.solve(solver)

あとは, 上記のようにsolverとしてCOIN_CMDを指定すれば, 自分でインストールしたCBCを使用することができます.

solver = pulp.PULP_CBC_CMD()とすると, PuLPと一緒にインストールされるCBCソルバを指定できます.



初期解の使用

解の探索時に実行可能解が既に分かっていると, 効率よく探索を行えて計算時間が短くなることが多いです.

初期解の情報が保存された.solファイルを用いて

solver = pulp.PULP_CBC_CMD(options=["mips", "initialSolution.sol"])

して, 実行します6.




GLPK

install

brew install glpk

インストールすると, 何もしなくてもpulpで使用できるようになります. ただCBCよりも遅い場合が多いため, 使用することはほぼなさそう.7




SCIP

SCIPとはドイツのZIBという研究機関が開発しているMIP(mixed integer programming)とMINLP(mixed integer nonlinear programmin)ソルバーです. 非商用利用(Academic License)であれば, 無料で使用することができます.

f:id:inarizuuuushi:20190309125049p:plain:w300

SCIPのダウンロード

公式サイトからDownloadからSCIP Optimization Suiteをダウンロードします.

README.mdに従って, makeを行います.

CMake

mkdir build
cd build
cmake ..
make
make test # チェック

cmakeがない場合はbrew install cmake でcmakeをインストールできます

Make

make
make test # チェック

pulpの設定ファイルへscipのパスを追記

scipのパスを通す

ln -s scipoptsuite-6.0.1(ダウンロードしたscipパス)/scip/bin/scip /usr/local/bin/scip

pulp.cfg.osxに下を上書きします . pulp.cfg.osx ファイルは
>> import pulp
>> pulp.__path__
で確認できるフォルダにあります

ScipPath = scip

PulpCbcPathの下に書き込むのがいいでしょう. これでscippulpで使えるようになります.

m = LpProblem()
solver = pulp.SCIP()
m.solve(solver)
  • デフォルトで使用できるオプション
    • msg (Default: False) Trueにすると実行ログが表示されます
    • keepFiles (Default: False) Trueにすると.lpファイルと.solファイルが保存されます.

並列化への対応

scipを並列化versionで使用するには, さらに細かい設定をする必要があります.

# scipoptsuite ディレクトリ上
make PARASCIP=true
make ug

を実行するとfscipが使用できるようになります(バイナリファイルは/ug/bin/にできます)

パスを通す

ln -s scipoptsuite-6.0.1(ダウンロードしたscipパス)/ug/bin/fscip /usr/local/bin/fscip

pulpfscipの実行パラメータファイルを置く

fscipを実行するにはパラメータ設定ファイルを指定する必要があります. これはpulpフォルダのsolverdir以下に置くのがいいでしょう.

# pulpフォルダ上
mkdir -p solverdir/scip/fscip
cp scipoptsuite-6.0.1(ダウンロードしたscipパス)/ug/setting/default.set solverdir/scip/fscip/default.set 

pulpフォルダは
>> import pulp
>>pulp.__path__
で確認できます

pulpの設定ファイル(pulp.cfg.osx)に追記する

FscipPath = fscip
FscipParamPath = %(here)s/solverdir/scip/fscip/default.set

PulpCbcPathもしくはScipPathの下がいいでしょう.

pulpsolvers.pyを書き換える

fscipscipの実行コマンド, ならびに最適化結果ファイルのフォーマットが異なるので 調節する必要があります. solvers.pydef initialize周辺とclass SCIP_CMDを以下のように書き換えます.

def initialize
initialize関数内

fscip_pathfscip_param_pathを追加.

try:
    scip_path = config.get("locations", "ScipPath")
except configparser.Error:
    scip_path = 'scip

↓↓↓↓↓↓↓↓

try:
    scip_path = config.get("locations", "ScipPath")
    fscip_path = config.get("locations", "FscipPath")
    fscip_param_path = config.get("locations", "FscipParamPath")
except configparser.Error:
    scip_path = 'scip

initialize返り値

fscip_path, fscip_para_pathを追加

return cplex_dll_path, ilm_cplex_license, ilm_cplex_license_signature,\
    coinMP_path, gurobi_path, cbc_path, glpk_path, pulp_cbc_path, scip_path

↓↓↓↓↓↓↓↓

return cplex_dll_path, ilm_cplex_license, ilm_cplex_license_signature,\
    coinMP_path, gurobi_path, cbc_path, glpk_path, pulp_cbc_path,\
    scip_path, fscip_path, fscip_param_path

if name == 'main'の何行か下

fscip_pathfscip_param_pathを追加

cplex_dll_path, ilm_cplex_license, ilm_cplex_license_signature, coinMP_path,\
        gurobi_path, cbc_path, glpk_path, pulp_cbc_path, scip_path\
        initialize(config_filename, operating_system, arch)

↓↓↓↓↓↓↓↓

cplex_dll_path, ilm_cplex_license, ilm_cplex_license_signature, coinMP_path,\
        gurobi_path, cbc_path, glpk_path, pulp_cbc_path,\
        scip_path, fscip_path, fscip_param_path = \
        initialize(config_filename, operating_system, arch)

class SCIP_CMD

class SCIP_CMD(LpSolver_CMD):
    """The SCIP optimization solver"""

    SCIP_STATUSES = {
        'unknown': LpStatusUndefined,
        'user interrupt': LpStatusNotSolved,
        'node limit reached': LpStatusNotSolved,
        'total node limit reached': LpStatusNotSolved,
        'stall node limit reached': LpStatusNotSolved,
        'time limit reached': LpStatusNotSolved,
        'memory limit reached': LpStatusNotSolved,
        'gap limit reached': LpStatusNotSolved,
        'solution limit reached': LpStatusNotSolved,
        'solution improvement limit reached': LpStatusNotSolved,
        'restart limit reached': LpStatusNotSolved,
        'optimal solution found': LpStatusOptimal,
        'infeasible':   LpStatusInfeasible,
        'unbounded': LpStatusUnbounded,
        'infeasible or unbounded': LpStatusNotSolved,
    }

    FSCIP_STATUSES = {
        'objective value': LpStatusOptimal, 
        'No Solution': LpStatusInfeasible
    }

    def defaultPath(self):
        return self.executableExtension(scip_path)
    
    def __init__(self, keepFiles = 0, msg = 0,
            threads = None):
        if threads is None or threads == 1:
            path = scip_path
        else:
            path = fscip_path
            self.fscip_param_path = fscip_param_path
        LpSolver_CMD.__init__(self, path, keepFiles, msg)
        self.threads = threads

    def available(self):
        """True if the solver is available"""
        return self.executable(self.path)

    def actualSolve(self, lp):
        """Solve a well formulated lp problem"""
        print(self.path)
        if not self.executable(self.path):
            raise PulpSolverError("PuLP: cannot execute "+self.path)

        # TODO: should we use tempfile instead?
        if not self.keepFiles:
            uuid = uuid4().hex
            tmpLp = os.path.join(self.tmpDir, "%s-pulp.lp" % uuid)
            tmpSol = os.path.join(self.tmpDir, "%s-pulp.sol" % uuid)
        else:
            tmpLp = lp.name + "-pulp.lp"
            tmpSol = lp.name + "-pulp.sol"

        lp.writeLP(tmpLp)
        if self.path == scip_path:
            proc = [
                self.path, '-c', 'read "%s"' % tmpLp, '-c', 'optimize',
                '-c', 'write solution "%s"' % tmpSol, '-c', 'quit'
            ]
        else:
            proc = [
                self.path, self.fscip_param_path, tmpLp,
                '-sth', '%d' % self.threads, '-fsol', tmpSol
            ]
        proc.extend(self.options)
        if not self.msg:
            proc.append('-q')
        print(' '.join(proc))
        if os.path.exists(tmpSol):
            os.remove(tmpSOl)
        self.solution_time = clock()
        subprocess.check_call(proc, stdout=sys.stdout, stderr=sys.stderr)
        self.solution_time += clock()

        if not os.path.exists(tmpSol):
            raise PulpSolverError("PuLP: Error while executing "+self.path)

        if self.path == scip_path:
            lp.status, values = self.scip_readsol(tmpSol)
        elif self.path == fscip_path:
            lp.status, values = self.fscip_readsol(tmpSol)

        # Make sure to add back in any 0-valued variables SCIP leaves out.
        finalVals = {}
        for v in lp.variables():
            finalVals[v.name] = values.get(v.name, 0.0)

        lp.assignVarsVals(finalVals)

        if not self.keepFiles:
            for f in (tmpLp, tmpSol):
                try:
                    os.remove(f)
                except:
                    pass

        return lp.status

    def scip_readsol(self, filename):
        """Read a SCIP solution file"""
        with open(filename) as f:
            # First line must containt 'solution status: <something>'
            try:
                line = f.readline()
                comps = line.split(': ')
                assert comps[0] == 'solution status'
                assert len(comps) == 2
            except:
                raise
                raise PulpSolverError("Can't read SCIP solver output: %r" % line)

            status = SCIP_CMD.SCIP_STATUSES.get(comps[1].strip(), LpStatusUndefined)
            if not status == LpStatusOptimal:
                return status, dict()
            # Look for an objective value. If we can't find one, stop.
            try:
                line = f.readline()
                comps = line.split(': ')
                assert comps[0] == 'objective value'
                assert len(comps) == 2
                float(comps[1].strip())
            except:
                raise PulpSolverError("Can't read SCIP solver output: %r" % line)

            # Parse the variable values.
            values = {}
            for line in f:
                try:
                    comps = line.split()
                    values[comps[0]] = float(comps[1])
                except:
                    raise PulpSolverError("Can't read SCIP solver output: %r" % line)

        return status, values

    def fscip_readsol(self, filename):
        """Read a SCIP solution file"""
        with open(filename) as f:
            try:
                line = f.readline()
                line = f.readline()
                comps = line.split(': ')
            except:
                raise 'Problem is unknonw states' 
            status = SCIP_CMD.FSCIP_STATUSES.get(comps[0].strip(), LpStatusUndefined)
            # Look for an objective value. If we can't find one, stop.
            if comps[0] == 'objective value':
                float(comps[1].strip())

            # Parse the variable values.
            values = {}
            for line in f:
                try:
                    comps = line.split()
                    values[comps[0]] = float(comps[1])
                except:
                    raise PulpSolverError("Can't read SCIP solver output: %r" % line)

        return status, values
SCIP = SCIP_CMD

これで並列化バージョンのSCIPを使用することができます. ただ, SCIPが出力する.solファイルの形式がcbcなどと違うため, 上のコードでは拾いきれていないエラーがあるかもしれません.




GUROBI

gurobipyを用いて実行されます.

GUROBI(mip = True, msg = True, timeLimit = None, epgap = None, **solverParams)
m = LpProblem()
solver = pulp.GUROBI(Threads=4)
m.solve(solver)
  • timeLimitは最大計算時間です.
  • またGUROBIの他のoption8も指定することができます.
    • Thread: 並列スレッド数など(solver = pulp.GUROBI(Thread=10))