サブロウ丸

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

[ Python ] LP modeler と Solver 一覧

Pythonで使線形計画問題(LP)を扱える最適化アプリケーションをいくつか使用することができます. アプリケーションは大きく2種類に分類されます.

  • Solver (ソルバー); 問題を解くアルゴリズムを内包したアプリケーション
  • Modeler (モデラー); 最適化問題をプログラムしやすくするアプリケーション. ユーザーとSolverの橋渡しを行う

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

LP (Modeling & Solver summary)

LP Solver / Modeler PuLP Scipy CVXOPT PICOS Pyomo
GLPK ok ok ok ok
CBC ok ok
SCIP ok? ok ok
CPLEX ok ok ok
MOSEK ok ok
GUROBI ok ok ok
CVXOPT ok ok
SCIPY ok

Solver Install

  • glpk
brew install glpk
pip install swiglpk
  • mosek
pip install mosek

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

Modeler Install

pip install pulp
  • Scipy
pip install scipy
  • CVXOPT
pip install cvxopt
  • PICOS
pip install picos
  • Pyomo
pip install pyomo
  • まとめて
pip install pulp scipy cvxopt picos pyomo

LP (Modeling)

以下の問題を解くためサンプルコードを載せています.

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

PuLP

from pulp import LpVariable, LpProblem, value

x = LpVariable('x')
y = LpVariable('y', lowBound=-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))

Available Solvers

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

from pulp import LpSolver
def AvailableSolver(cls):
    for c in cls.__subclasses__():
        try:
            AvailableSolver(c)
            if c().available():
                print(c)
        except:
            pass
AvailableSolver(LpSolver)
<class 'pulp.apis.coin_api.PULP_CBC_CMD'>
<class 'pulp.apis.glpk_api.GLPK_CMD'>
<class 'pulp.apis.choco_api.PULP_CHOCO_CMD'>
<class 'pulp.apis.scip_api.SCIP_CMD'>
<class 'pulp.apis.mosek_api.MOSEK'>
from pulp import PULP_CBC_CMD # GLPK, SCIP
solver = PULP_CBC_CMD()  # solver=GLPK(), solver=SCIP()
prob.solve(solver=solver)

Scipy

min cT x, s.t. Ax <= b

from scipy.optimize import linprog

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

res = linprog(c, A_ub=A, b_ub=b, bounds=bounds)
print(res)
     con: array([], dtype=float64)
     fun: -21.99999984082497
 message: 'Optimization terminated successfully.'
     nit: 6
   slack: array([3.89999997e+01, 8.46872172e-08])
  status: 0
 success: True
       x: array([ 9.99999989, -2.99999999])

CVXOPT

min cT x s.t. Gx <= h, Ax = b

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'])
     pcost       dcost       gap    pres   dres   k/t
 0:  7.7458e+00 -2.2000e+01  2e+01  0e+00  2e+00  1e+00
 1: -3.3531e+00 -9.3245e+00  5e+00  3e-16  4e-01  1e+00
 2: -2.2577e+00 -1.0334e+01  1e+01  3e-16  5e-01  2e+00
 3: -1.3020e+01 -1.8365e+01  1e+01  9e-16  4e-01  2e+00
 4: -1.9962e+01 -2.0107e+01  2e+00  2e-15  5e-02  9e-01
 5: -2.1977e+01 -2.1973e+01  3e-02  1e-15  8e-04  2e-02
 6: -2.2000e+01 -2.2000e+01  3e-04  2e-16  8e-06  2e-04
 7: -2.2000e+01 -2.2000e+01  3e-06  1e-15  8e-08  2e-06
Optimal solution found.
-21.99999773873395 [ 1.00e+01]
[-3.00e+00]

Available Solvers

conelp(CVXOPT), glpk, mosek

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

sol = solvers.lp(c, G, h, solver='mosek') # solver='glpk'
print(sol['primal objective'], sol['x'])
Problem
  Name                   :
  Objective sense        : min
  Type                   : LO (linear optimization problem)
  Constraints            : 3
  Cones                  : 0
  Scalar variables       : 2
  Matrix variables       : 0
  Integer variables      : 0

Optimizer started.
Presolve started.
Linear dependency checker started.
Linear dependency checker terminated.
Eliminator started.
Freed constraints in eliminator : 1
Eliminator terminated.
Eliminator started.
Freed constraints in eliminator : 0
Eliminator terminated.
Eliminator - tries                  : 2                 time                   : 0.00
Lin. dep.  - tries                  : 1                 time                   : 0.00
Lin. dep.  - number                 : 0
Presolve terminated. Time: 0.00
Optimizer terminated. Time: 0.01



Interior-point solution summary
  Problem status  : PRIMAL_AND_DUAL_FEASIBLE
  Solution status : OPTIMAL
  Primal.  obj: -2.2000000000e+01   nrm: 3e+01    Viol.  con: 0e+00    var: 0e+00
  Dual.    obj: -2.2000000000e+01   nrm: 6e+00    Viol.  con: 0e+00    var: 0e+00

Basic solution summary
  Problem status  : PRIMAL_AND_DUAL_FEASIBLE
  Solution status : OPTIMAL
  Primal.  obj: -2.2000000000e+01   nrm: 3e+01    Viol.  con: 0e+00    var: 0e+00
  Dual.    obj: -2.2000000000e+01   nrm: 6e+00    Viol.  con: 0e+00    var: 0e+00
-22.0 [ 1.00e+01]
[-3.00e+00]

CVXOPT (modeling)

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)
     pcost       dcost       gap    pres   dres   k/t
 0:  7.7458e+00 -2.2000e+01  2e+01  0e+00  2e+00  1e+00
 1: -3.3531e+00 -9.3245e+00  5e+00  3e-16  4e-01  1e+00
 2: -2.2577e+00 -1.0334e+01  1e+01  3e-16  5e-01  2e+00
 3: -1.3020e+01 -1.8365e+01  1e+01  9e-16  4e-01  2e+00
 4: -1.9962e+01 -2.0107e+01  2e+00  2e-15  5e-02  9e-01
 5: -2.1977e+01 -2.1973e+01  3e-02  1e-15  8e-04  2e-02
 6: -2.2000e+01 -2.2000e+01  3e-04  2e-16  8e-06  2e-04
 7: -2.2000e+01 -2.2000e+01  3e-06  1e-15  8e-08  2e-06
Optimal solution found.
[-2.20e+01]
 [ 1.00e+01]
 [-3.00e+00]

Available Solvers

conelp(CVXOPT), glpk, mosek

sol = solvers.lp(c, G, h, solver='mosek') # solver='glpk'
print(sol['primal objective'], sol['x'])
Problem
  Name                   :
  Objective sense        : min
  Type                   : LO (linear optimization problem)
  Constraints            : 3
  Cones                  : 0
  Scalar variables       : 2
  Matrix variables       : 0
  Integer variables      : 0

Optimizer started.
Presolve started.
Linear dependency checker started.
Linear dependency checker terminated.
Eliminator started.
Freed constraints in eliminator : 1
Eliminator terminated.
Eliminator started.
Freed constraints in eliminator : 0
Eliminator terminated.
Eliminator - tries                  : 2                 time                   : 0.00
Lin. dep.  - tries                  : 1                 time                   : 0.00
Lin. dep.  - number                 : 0
Presolve terminated. Time: 0.00
Optimizer terminated. Time: 0.01



Interior-point solution summary
  Problem status  : PRIMAL_AND_DUAL_FEASIBLE
  Solution status : OPTIMAL
  Primal.  obj: -2.2000000000e+01   nrm: 3e+01    Viol.  con: 0e+00    var: 0e+00
  Dual.    obj: -2.2000000000e+01   nrm: 6e+00    Viol.  con: 0e+00    var: 0e+00

Basic solution summary
  Problem status  : PRIMAL_AND_DUAL_FEASIBLE
  Solution status : OPTIMAL
  Primal.  obj: -2.2000000000e+01   nrm: 3e+01    Viol.  con: 0e+00    var: 0e+00
  Dual.    obj: -2.2000000000e+01   nrm: 6e+00    Viol.  con: 0e+00    var: 0e+00
-22.0 [ 1.00e+01]
[-3.00e+00]

PICOS (solver CVXOPT)

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

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

Pyomo

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())

Available Solvers

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