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