サブロウ丸

Sabrou-mal サブロウ丸

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

Pulp ソルバー選択 / 並列計算 (python)

python pulp CBC

Pulpについて

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

install

pip install pulp

pulpの処理フロー

使用できるソルバ一覧

現在の環境で使用できるソルバーは下記のコードで確認できます.

import pulp
def AvailableSolver(cls):
    for c in cls.__subclasses__():
        try:
            AvailableSolver(c)
            if c().available():
                print(c)
        except:
            pass
AvailableSolver(pulp.LpSolver)
* 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（有料）などをインストールすると、以下のようにソルバーを変更できます 。

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ソルバーで問題が解かれます.

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に設定します.
  • cuts (Default: True)
    • cut(切除平面) generator を行なって探索をより効率的に行います.
  • preserve (Default: True)
    • モデルを分析して、冗長な制約や固定できる変数や範囲をより強く制限できる変数などを見つけます. 問題を最初に解決するには、効果がないことがわかっていない限り、これを行う価値があります.
  • strong (Default: 5)
    • 強分岐で調べる変数の数を決定します. どの変数に分岐するかを決めるために、この数までの満たされていない変数を選び、最小の上下の分岐を試みます.
  • fracGapは計算終了の判定となる, 解の上界と下界のGapのしきい値を指定できます.
  • maxSecondsは最大計算時間を指定できます. (単位は秒)
  • threadsはソルバー実行の並列数を指定できます.
  • optionsではそれ以外のCBCオプションを設定できます.
    • 例えばmaxsolオプションを使用する場合は options = ['maxsol 1']とします.
      • maxsolオプションは, 暫定解が何個見つかった時に計算を終了するかを指定できます. 最適解でなくても実行可能解が１つ欲しい場合はmaxsol 1とすればより速い時間で求めることができます.

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

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

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

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

CBC > 1. 公式サイト*2から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/配下にバージョンごとにディレクトリが作成されます.

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

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

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

CBC > 1. cbcソルバにパスを通す

maxOSの環境で話を進めます。

ソースコードからmakeでコンパイルした場合はmake installでusr/local/binにバイナリがコピーされるのでこれだけで十分です。

また、homebrewでinstallした場合は下記のようにパスを通せば, cbcコマンドを使えるようになります. homebrewでは特にインストール先を変更していない場合は/usr/local/Cellar/配下にバージョンごとにディレクトリが作成されます.[^1] なので

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

を実行してパスを張ります. (.bashrcや.zshrcに記載すると良いでしょう) これでターミナル上で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によりこのソルバが使用されます.

CBC > 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

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ソルバを指定できます.

CBC > 初期解の使用

解の探索時に実行可能解が既に分かっていると, 効率よく探索を行えて計算時間が短くなることが多いです. 初期解の情報が保存された.solファイルを用いて

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

して, 実行します*3.

GLPK

install

brew install glpk

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

SCIP

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

SCIP > install

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

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

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

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

make
make test # チェック
make install

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

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

ScipPath = scip

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

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

SCIP > 並列化への対応

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

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

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

pulpにfscipの実行パラメータファイルを置く

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の下がいいでしょう.

SCIP > pulpのsolvers.pyを書き換える

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

def initialize
initialize関数内

fscip_pathとfscip_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_pathとfscip_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の他のoption*5も指定することができます.
    • Thread: 並列スレッド数など(solver = pulp.GUROBI(Thread=10))

*1:https://www.gams.com/latest/docs/S_CBC.html#CBC_general

*2:https://projects.coin-or.org/Cbc

*3:http://miti-7.hatenablog.com/entry/2019/01/25/153134

*4:https://www.slideshare.net/tkm2261/ss-42149384

*5:https://www.gurobi.com/documentation/8.1/refman/parameters.html