Brief Introduction to MILP Solvers

Posted on In

MILP -- Mixed Integer Linear Programming

Highly recommended: CVXPY

CVXPY is an interface which is compatible with multiple solvers like CBC, CPLEX, GUROBI, etc.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
import numpy as np
import cvxpy as cp
n = 10
np.random.seed(3)
a = np.random.randint(1, 10, size=n)
b = np.random.randint(1, 10, size=n)

# Variables
x = cp.Variable(shape=n, boolean=True)
# Objective
objective = cp.Minimize(x @ a)
# Constraints
constraints = []
constraints.append(x @ b >= 10)

# Build & Solve Problem
problem = cp.Problem(objective, constraints)
problem.solve(solver=cp.CPLEX, verbose=0)

# Print Result
print(problem.status)
print("x: ", x.value)
print("Optimal value: ", problem.value)

Recommended: PULP

PULP is also an interface which is compatible with multiple solvers like CBC, CPLEX, GUROBI, etc.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
from pulp import LpProblem, LpMinimize, lpSum, LpVariable, LpBinary, LpStatus, value
from pulp.apis import PULP_CBC_CMD
import numpy as np
n = 10
np.random.seed(3)
a = np.random.randint(1, 10, size=n)
b = np.random.randint(1, 10, size=n)

# Variables
x = np.array([LpVariable('x_{0}'.format(i), cat=LpBinary) for i in range(n)])
# Objective
prob = LpProblem("PULPDemo", LpMinimize)
prob += lpSum(x @ a)
# Constraints
prob += lpSum(x @ b) >= 10

# Build & Solve Problem
prob.solve(solver=PULP_CBC_CMD(threads=4, msg=0))

# Print Result
print(LpStatus[prob.status])
v = [v.varValue for v in prob.variables()]
print("x: ", v)
print("Optimal value: ", value(prob.objective))

Give It a Try: GEATPY

GEATPY is a genetic and evolutionary algorithm toolbox. If you are not an expert in convex optimization and linear programming, why not give it a try. Genetic algorithm is always a good solution to optimization problem.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
import geatpy as ea
import numpy as np
n = 10
np.random.seed(3)
a = np.random.randint(1, 10, size=n)
b = np.random.randint(1, 10, size=n)

r = 1
@ea.Problem.single
def evalVars(Vars):
    # Objective
    f = Vars @ a
    # Constraints
    CV = np.array([10 - Vars @ b])
    return f, CV

# Build Problem
problem = ea.Problem(name='GEATPYDemo',
                        M=1, # objective dimension
                        maxormins=[1], # minimize
                        Dim=n, # variable dimension
                        varTypes=[1 for _ in range(n)], # variables are integers
                        lb=[0 for _ in range(n)], # lower bound of variables
                        ub=[1 for _ in range(n)], # upper bound of variables
                        evalVars=evalVars)
algorithm = ea.soea_SEGA_templet(problem,
                                    ea.Population(Encoding='RI', NIND=200),
                                    MAXGEN=500, # max generation 
                                    logTras=1,  # record log each logTras generation
                                    trappedValue=1e-6,
                                    maxTrappedCount=100)
# Solve Problem
res = ea.optimize(algorithm, seed=1, verbose=0, drawing=0, outputMsg=1, drawLog=0, saveFlag=0, dirName='result')

For MILFP Problems

Objective: \(maximize(\frac{x * a}{x * b})\) where \(x\) is variables and \(a\), \(b\) are matrices. This ratio can be used to represent efficiency.

For a Mixed Integer Linear Fractional Problem (MILFP), you may try to implement Charnes-Cooper tansformation which can convert a MILFP problem to a MILP problem otherwise you may try in CVXPY.

1
problem.solve(solver=cp.CPLEX, verbose=0, qcp=1)

By default CVXPY only handles DQCP problem, but it still has the ability to handle QCP problem like MIFLP problem. Warning: not all solvers support QCP problem.

A good solver is a handy tool to solve LP problem. Still it takes experience and skill to do some linear transformation so that you can solve real world problem.