Python Mixed Integer Linear Programming

Python Mixed Integer Linear Programming

Pulp is a python modeling interface that hooks up to solvers like CBC(open source), CPLEX (commercial), Gurobi(commercial), XPRESS-MP(commercial) and YALMIP(open source).

You can also use Pyomo to model the optimization problem and then call an external solver, namely CPLEX, Gurobi GLPK and the AMPL solver library.

You can also call GLPK from GLPK/Python, PyGLPK or PyMathProg.

Yet another modelling language is CMPL, which has a python interface for MIP solvers (for linear programs only).

All the above solvers solve Mixed Integer Linear Programs, while some of them (CPLEX, GUROBI and XRESS-MP for sure) can solve Mixed Integer Quadratic Programs and Quadratically constrained quadratic programs (and also conic programs but this probably goes beyond the scope of this question).

MIP refers to Mixed integer programs, but it is commonly used to refer to linear programs only. To make the terminology more precise, one should always refer to MILP or MINLP (Mixed integer non-linear programming).

Note that CPLEX and GUROBI have their own python APIs as well, but they (and also) XPRESS-MP are commercial products, but free for academic research. CyLP is similar to Pulp above but interfaces with the COIN-OR solvers CBC and CGL and CLP.

Note that there is a big difference in the performance of commercial and free solvers: the latter are falling behind the former by a large margin. SCIP is perhaps the best non-commercial solver (see below for an update). Its python interface, PySCIPOpt, is here.

Also, have a look at this SO question.

Finally, if you are interested at a simple constraint solver (not optimization) then have a look at python-constraint.

UPDATES

More solvers and python interfaces that fell into my radar:

Update: MIPCL links appear to be broken.

MIPCL, which appears to be the fastest non-commercial MIP solver, has a python interface that has quite good documentation. Note, however, that the Python API does not include the advanced functionality that comes together with the native MIPCLShell. I particularly like the MIPCL-PY manual, which demonstrates an array of models used in Operations Management, on top of some small-scale implementations. It is a very interesting introductory manual in its own right, regardless of which solver/API one may want to make use of.

Google Optimization Tools, which include a multitude of functionalities, such as

• A constraint programming solver and a linear programming (not MIP) solver
• An interface for MIP solvers (supports CBC, CLP, GLOP, GLPK, Gurobi, CPLEX, and SCIP)
• Specialized algorithms for graphs, for the Travelling Salesman Problem, the Vehicle Routing problem and for Bin packing & Knapsack problems

It has extensive documentation of several traditional OR problems and simple implementations. I could not find a complete Python API documentation, although there exist some examples here. It is somewhat unclear to me how other solvers hook up on the interface and whether methods of these solvers are available.

CVXOPT, an open-source package for convex optimization, which interfaces to GLPK (open source) and MOSEK
(commercial). It is versatile, as it can tackle many problem classes (notably linear, second-order, semidefinite, convex nonlinear). The only disadvantage is that it modeling complex problems may be cumbersome, as the user needs to pass the data in a "Matlab-y" fashion (i.e., to specify the matrix, rhs vectors, etc). However, it can be called from the modeling interfaces PICOS and...

CVXPY, a python-embedded optimization language for convex optimization problems, which contains CVXOPT as a default solver, but it can hook up to the usual MIP solvers.

Thanks to RedPanda for pointing out that CVXOPT/CVXPY support MIP solvers as well.

For a very comprehensive article on optimization modeling capabilities of packages and object-oriented languages (not restricted to Python), check this article.

mixed integer programming optimization

The problem is currently unbounded (see Objective: -1.E+15).

Use m.Intermediate() instead of m.MV(). An MV (Manipulated Variable) is a degree of freedom that the optimizer can use to achieve an optimal objective among all of the feasible solutions. Because tempo_b1, tempo_b2, and tempo_total all have equations associated with solving them, they need to either be:

• Regular variables with m.Var() and a corresponding m.Equation() definition
• Intermediate variables with m.Intermediate() to define the variable and equation with one line.

Here is the solution to the simple Mixed Integer Linear Programming (MINLP) optimization problem.

 ----------------------------------------------------------------
 APMonitor, Version 1.0.1
 APMonitor Optimization Suite
 ----------------------------------------------------------------
 
 --------- APM Model Size ------------
 Each time step contains
   Objects      :            0
   Constants    :            0
   Variables    :            7
   Intermediates:            2
   Connections  :            0
   Equations    :            6
   Residuals    :            4
 
 Number of state variables:              7
 Number of total equations: -            3
 Number of slack variables: -            0
 ---------------------------------------
 Degrees of freedom       :              4
 
 ----------------------------------------------
 Steady State Optimization with APOPT Solver
 ----------------------------------------------
Iter: 1 I:  0 Tm: 0.00 NLPi: 4 Dpth: 0 Lvs: 0 Obj: -3.00E+03 Gap: 0.00E+00
 Successful solution
 
 ---------------------------------------------------
 Solver         :  APOPT (v1.0)
 Solution time  :   1.529999999911524E-002 sec
 Objective      :   -3000.00000000000     
 Successful solution
 ---------------------------------------------------

Braco_gas_1:[0.0]
Braco_gas_2:[1.0]
Braco_eh_1:[0.0]
Braco_eh_2:[1.0]
vazao_gas:[150.0]
vazao_eh:[150.0]
volume_gas:1000
volume_eh:2000
tempo_b1:[1000.0]
tempo_b2:[2000.0]
tempo_total:[3000.0]

The complete script:

from gekko import GEKKO
import numpy as np

m = GEKKO(remote=False)

# binary variables
braco_gas_1 = m.Var(integer=True, lb=0,ub=1)
braco_eh_1 = m.Var(integer=True, lb=0,ub=1)

braco_gas_2 = m.Var(integer=True, lb=0,ub=1)
braco_eh_2 = m.Var(integer=True, lb=0,ub=1)

# continuous variables
vazao_gas = m.Var(value=100,lb=0,ub=150)
vazao_eh = m.Var(value=100,lb=0,ub=150)

#constants
volume_gas = 1000
volume_eh = 2000

# I want to see each parcel of the objective
tempo_b1 = m.Intermediate(vazao_gas*braco_gas_1 + volume_gas*braco_gas_2)
tempo_b2 = m.Intermediate(vazao_eh*braco_eh_1 + volume_eh*braco_eh_2)

# that is supposed to be the objective
tempo_total = m.Var(lb=0, ub = 4000)
m.Equation(tempo_total==tempo_b1+tempo_b2)

# Only of binary variable of each group can be true
m.Equation (braco_gas_1+braco_gas_2 == 1)
m.Equation (braco_eh_1+braco_eh_2 == 1)

# I want to maximize the objective
m.Maximize(tempo_b1+tempo_b2)

m.options.SOLVER = 1
m.solve()    # solve

print('Braco_gas_1:'+str(braco_gas_1.value))
print('Braco_gas_2:'+str(braco_gas_2.value))
print('Braco_eh_1:'+str(braco_eh_1.value))
print('Braco_eh_2:'+str(braco_eh_2.value))
print('vazao_gas:'+str(vazao_gas.value))
print('vazao_eh:'+str(vazao_eh.value))
print('volume_gas:'+str(volume_gas))
print('volume_eh:'+str(volume_eh))
print('tempo_b1:'+str(tempo_b1.value))
print('tempo_b2:'+str(tempo_b2.value))

print('tempo_total:'+str(tempo_total.value))

Additional tutorials are available in the documentation or in the 18 example problems (with videos).

Integer division in Mixed Integer Linear Programming

model.addRange(val - divVal*res, 0, 0.99999, name="Range")

I would prefer to use this above mentioned Range Constraint only. Incorporating tighter bounds (there is only integer in the given range, which we require) directly into the model can not only improve the numerical behavior, but it also speed up the optimization process (because gurobi use branch and bound algorithm to get solutions)
https://www.gurobi.com/documentation/9.1/refman/improving_ranges_for_varia.html

Optimality - A small change in model can easily calculate an optimal result, adding res in a Minimisation type objective function or negative of it in the maximisation function would shrink its value at lower side if divVal*res will become integer. Gurobi does not provide Less than constraint. Moreover, An integrality restriction on a variable is considered satisfied in Gurobi when the variable's value is less than IntFeasTol from the nearest integer value. Default value of IntFeasTol tolerance is 1e-5, and it could be further reduced up to 1e-9 for a better result. However, making multi-objective model, add extra complexity to the model. I would not like to recommend it.

model.addRange(val - divVal*res, 0, 1, name="Range")

model.setObjective(res, GRB.MINIMIZE)

---

Related Topics