Pythonで使える 線形計画(LP)ソルバ と モデラの 一覧

Pythonには使線形計画問題(LP)を扱える最適化アプリケーションは大きく以下の2種類に分類されます。

- Solver (ソルバー);
  - 問題を解くアルゴリズムを内包したアプリケーション
- Modeler (モデラー);
  - 最適化問題をプログラムしやすくするアプリケーション。 ユーザーとSolverの橋渡しを行う
  - ユーザーがモデラーで問題を記述し、モデラーがソルバーを呼び出す

使用可能なSolverとModelerの組み合わせと、Modelerの簡単な実装例をまとめました。

## 履歴

- 2020/05/26: 初稿
- 2024/09/13: モデラーにPython-MIPを追加。ソルバーにHiGHSを追加。


# モデラーとソルバのまとめ

行がソルバーで、列方向がモデラーです。

| LP Solver / Modeler                   | Python-MIP | PuLP | Scipy | CVXOPT | PICOS | Pyomo |
| -------------------                   | ---------- | ---- | ----- | ------ | ----- | ----- |
| GLPK                                  |            | ok   |       | ok     | ok    | ok    |
| [CBC](https://github.com/coin-or/Cbc) | ok         | ok   |       |        |       | ok    |
| [HiGHS](https://highs.dev)            | ok         | ok   |       |        |       |       |
| SCIP                                  | ok         | ok?  |       |        | ok    | ok    |
| CPLEX                                 | ok         | ok   |       |        | ok    | ok    |
| MOSEK                                 |            |      |       | ok     | ok    |       |
| GUROBI                                | ok         | ok   |       |        | ok    | ok    |
| CVXOPT                                |            |      |       | ok     | ok    |       |
| SCIPY                                 |            |      | ok    |        |       |       |


## ソルバーのインストール

mac環境を想定

- glpk

非商用

```sh
brew install glpk
pip install swiglpk
```

- mosek

商用、アカデミックライセンスあり

```sh
pip install mosek
```

and get academic license from https://www.mosek.com/products/academic-licenses/ if you can.


## Modeler Install

- Python-MIP、PuLp、Scipy、CVXOPT、PICOS、Pyomo

```sh
pip install mip pulp scipy cvxopt picos pyomo mip
```

# LP (Modeling)

次の問題をそれぞれのモデラーを使って解いてみます

```
min -x+4y
s.t. -3x +  y ≤ 6
      -x - 2y ≥ -4
            y ≥ -3
```

## Python-MIP

```python
solver = mip.HIGHS
# solver = mip.CBC
m = mip.Model(solver_name=solver)
x = m.add_var()
y = m.add_var()
m.objective = -3 * x + y
m.add_constr(-x -2 * y >= -4)
m.add_constr(y >= -3)
m.optimize()

print("obj", m.objective.x, "x", x.x, "y", y.x)
```

### 使用可能なソルバー

```python
# solvers
CBC = "CBC"
CPX = "CPX"  # we plan to support CPLEX in the future
CPLEX = "CPX"  # we plan to support CPLEX in the future
GRB = "GRB"
GUROBI = "GRB"
HIGHS = "HiGHS"
SCIP = "SCIP"  # we plan to support SCIP in the future
```


https://github.com/coin-or/python-mip/blob/0ccb81115543e737ab74a4f1309891ce5650c8d5/mip/constants.py#L26-L33


## PuLP


```python
from pulp import LpVariable、LpProblem, value

x = LpVariable('x')
y = LpVariable('y', lowBound=-3)  # y ≥ -3

prob = LpProblem()
prob += -x + 4*y  # 目的関数
prob += -3*x + y <= 6
prob += -x - 2*y >= -4

prob。solve()
print(value(prob.objective), value(x), value(y))
```


### 使用可能なソルバーの確認方法

CBC、CPLEX、GLPK、GUROBI、XPRESS、(CHOCO)

使用可能なソルバーの表示.
(参考; http://www.logopt.com/download/OptForPuzzle.pdf)

```python
from pulp import LpSolver
def AvailableSolver(cls):
    for c in cls.__subclasses__():
        try:
            AvailableSolver(c)
            if c().available():
                print(c)
        except:
            pass
AvailableSolver(LpSolver)
```

```python
from pulp import PULP_CBC_CMD # GLPK, SCIP
solver = PULP_CBC_CMD()  # solver=GLPK(), solver=SCIP()
prob.solve(solver=solver)
```


## Scipy

min c^T x, s.t. Ax <= b

定式化を以下の形に変更

```
min -x+4y
s.t. -3x +  y ≤ 6
       x + 2y ≤ 4
           -y ≤ 3
```


```python
from scipy.optimize import linprog

c = [-1, 4]
A = [[-3, 1], [1, 2]]
b = [6, 4]
bounds =[
    (None, None),  # -∞ ≤ x ≤ ∞
    (-3, None)     # -3 ≤ y ≤ ∞
]

res = linprog(c, A_ub=A, b_ub=b, bounds=bounds)
print(res)
```


## CVXOPT

min c^T x s.t. Gx <= h、Ax = b

定式化を以下の形に変更

```
min -x+4y
s.t. -3x +  y ≤ 6
       x + 2y ≤ 4
      0x + -y ≤ 3
```


```python
from cvxopt import matrix, solvers

c = matrix([-1., 4.])
G = matrix([[-3., 1., 0.], [1., 2., -1.]])
h = matrix([6., 4., 3.])

sol = solvers.lp(c, G, h)
print(sol['primal objective'], sol['x'])
```


### 使用可能なソルバー

conelp(CVXOPT)、glpk、mosek

https://cvxopt.org/userguide/coneprog.html#s-lpsolver


```python
sol = solvers.lp(c, G, h, solver='mosek') # solver='glpk'
print(sol['primal objective'], sol['x'])
```



## CVXOPT (modeling)


```python
from cvxopt import matrix
from cvxopt.modeling import variable, op

x = variable()
y = variable()

c1 = ( -3*x + y <= 6 )
c2 = ( -x - 2*y >= -4 )
c3 = ( y >= -3 )

lp = op(-x+4*y, [c1,c2,c3])

lp.solve()
print(lp.objective.value(), x.value, y.value)
```



### 使用可能なソルバー

conelp(CVXOPT)、glpk、mosek


```python
sol = solvers.lp(c, G, h, solver='mosek') # solver='glpk'
print(sol['primal objective'], sol['x'])
```



## PICOS (solver CVXOPT)


```python
import picos as pic

prob = pic.Problem()

x = prob.add_variable('x')
y = prob.add_variable('y')

prob.set_objective('min', -x+4*y)
prob.add_constraint(-3*x + y <= 6)
prob.add_constraint(-x - 2*y >= -4)
prob.add_constraint(y >= -3)

sol = prob.solve()
print(sol.value, x.value, y.value)
```


    /usr/local/lib/python3.7/site-packages/ipykernel_launcher.py:5: DeprecationWarning: Problem.add_variable is deprecated: Variables can now be created independent of problems, and do not need to be added to any problem explicitly.
      """
    /usr/local/lib/python3.7/site-packages/ipykernel_launcher.py:6: DeprecationWarning: Problem.add_variable is deprecated: Variables can now be created independent of problems, and do not need to be added to any problem explicitly.



### Available Solvers

- None to let PICOS choose.
- "cplex" for CPLEXSolver.
- "cvxopt" for CVXOPTSolver.
- "ecos" for ECOSSolver.
- "glpk" for GLPKSolver.
- "gurobi" for GurobiSolver.
- "mosek" for MOSEKSolver.
- "mskfsn" for MOSEKFusionSolver.
- "scip" for SCIPSolver.
- "smcp" for SMCPSolver.

https://picos-api.gitlab.io/picos/api/picos.modeling.options.html#option-solver


```python
sol = prob.solve(solver='glpk')
print(sol.value, x.value, y.value)
```

## Pyomo


```python
import pyomo.environ as pyo

model = pyo.ConcreteModel()

model.x = pyo.Var()
model.y = pyo.Var()

model.OBJ = pyo.Objective(expr=-model.x + 4*model.y)

# add constraints
model.Constraint1 = pyo.Constraint(expr=-3*model.x+model.y <= 6)
model.Constraint2 = pyo.Constraint(expr=-model.x-2*model.y >= -4)
model.Constraint3 = pyo.Constraint(expr=model.y >= -3)

# add constraints
# def rule(model, i):
#     A = [[-3, 1], [1, 2], [0, -1]]
#     b = [6, 4, 3]
#     return A[i][0]*model.x + A[i][1]*model.y <= b[i]

# model.const_set = pyo.Set(initialize=[0, 1, 2])
# model.CONSTRAINTS = pyo.Constraint(model.const_set, rule=rule)

opt = pyo.SolverFactory('glpk')
res = opt.solve(model)
print(model.OBJ(), model.x(), model.y())
```


### 使用可能なソルバー

```
Serial Solver Interfaces
------------------------
The serial, pyro and phpyro solver managers support the following
solver interfaces:

    asl                  + Interface for solvers using the AMPL Solver
                           Library
    baron                  The BARON MINLP solver
    bilevel_blp_global   + Global solver for continuous bilevel linear
                           problems
    bilevel_blp_local    + Local solver for continuous bilevel linear
                           problems
    bilevel_bqp          + Global solver for bilevel quadratic
                           problems
    bilevel_ld           + Solver for bilevel problems using linear
                           duality
    cbc                    The CBC LP/MIP solver
    conopt                 The CONOPT NLP solver
    contrib.gjh            Interface to the AMPL GJH "solver"
    cplex                  The CPLEX LP/MIP solver
    cplex_direct           Direct python interface to CPLEX
    cplex_persistent       Persistent python interface to CPLEX
    gams                   The GAMS modeling language
    gdpbb                * Branch and Bound based GDP Solver
    gdpopt               * The GDPopt decomposition-based Generalized
                           Disjunctive Programming (GDP) solver
    glpk                 * The GLPK LP/MIP solver
    gurobi                 The GUROBI LP/MIP solver
    gurobi_direct          Direct python interface to Gurobi
    gurobi_persistent      Persistent python interface to Gurobi
    ipopt                  The Ipopt NLP solver
    mindtpy              * MindtPy: Mixed-Integer Nonlinear
                           Decomposition Toolbox in Pyomo
    mosek                * Direct python interface to Mosek
    mpec_minlp           + MPEC solver transforms to a MINLP
    mpec_nlp             + MPEC solver that optimizes a nonlinear
                           transformation
    multistart           * MultiStart solver for NLPs
    path                   Nonlinear MCP solver
    pico                   The PICO LP/MIP solver
    ps                   * Pyomo's simple pattern search optimizer
    py                   + Direct python solver interfaces
    scip                   The SCIP LP/MIP solver
    trustregion          * Trust region filter method for black
                           box/glass box optimization
    xpress                 The XPRESS LP/MIP solver
```

https://pyomo.readthedocs.io/en/stable/solving_pyomo_models.html


## その他とまとめ

1. Scipy は 制約を行列の形で入力しないといけないので手間
2. PuLPとCVXOPTはSolverも一緒にinstallされるのでモデラをインストールするだけで、すぐに問題を解くことができます
3. CVXOPTとPICOSは線形計画以外の問題(二次計画など、CVXOPTはそちらがメイン)を解くことができるので、汎用性が高い
4. Pyomoでの実装はクセが強い