Optimization modeling in Python: Pulp, Gurobi and CPLEX
I have been involved in design, development and implementation of operations research (OR) and optimization models such as Linear Program (LP), Mixed Integer Linear Program (MILP) and Quadratic Program (QP) for more than a decade. In the past 4 years, I have realized the importance of OR solutions, i.e., software solutions that are based on optimization models. In the past, we used to model a real-world optimization problem with LP/MILP packages such as GAMS, AMPL, OPL, etc in isolation, solve it with an optimization solver (such as CPLEX, Gurobi, Mosek, Xpress, etc.) and give the optimal result to managers and decision makers. Thus, optimization models were very well adopted for strategical/ tactical decisions rather than operational decisions.
However, these days, many in industry want to plan and make optimal decisions regularly as a part of their hourly, daily or weekly operations. Recent computational advances have provided required infrastructures for us to incorporate optimization models in analytic software solutions. This means that todays OR practitioners need to design, model and implement robust software engines that are based on LP/MILP models. They need to utilize a programming language such as C++, Java, C#, Python, etc. for that purpose.
A good and popular programming language recommended by many in OR and Data Science community is Python. It is easy, flexible yet powerful, and has great libraries for Machine Learning, Optimization and Statistical modeling. In this blog, I plan to show how one can use Python and write OR models (LP/MILP).
Many optimization solvers (commercial and open-source) have Python interfaces for modeling LP, MILP, and QPs. I select CPLEX and Gurobi since they are among leading commercial solvers, and PuLp which is a powerful open source modeling package in Python. I provide a side by side tutorial for each of these packages, and I hope it will help you to easily translate your model from one to another. Here, we use “gurobipy” which is Gurobi’s Python API, “docplex” which is IBM Decision Optimization CPLEX Modeling for Python, and “pulp” which is an LP/MILP modeler written in python. I assume a reader of this blog is familiar with Python, i.e., knows how to install and use python packages and use python data structures like lists, tuples and dictionaries, and also knows what linear programming, mixed integer programming and constraints optimization are.
Now lets dive in to optimization modeling with Gurobi, CPLEX and PuLp. An optimization model is a problem which has an objective (or a set of objectives in multi-objective programming), a set of constraints and a set of decision variables. The following is a simple optimization model.
In the above optimization example, n, m, a, c, l, u and b are input parameters and given. In order to write Python code, we set these parameters as follows:
import randomn = 10
m = 5
set_I = range(1, n+1)
set_J = range(1, m+1)
c = {(i,j): random.normalvariate(0,1) for i in set_I for j in set_J}
a = {(i,j): random.normalvariate(0,5) for i in set_I for j in set_J}
l = {(i,j): random.randint(0,10) for i in set_I for j in set_J}
u = {(i,j): random.randint(10,20) for i in set_I for j in set_J}
b = {j: random.randint(0,30) for j in set_J}
Now, its time to implement our OR model in Python. When we want to code up an optimization model, we put a place holder for that model (like a blank canvas), then add its elements (decision variables and constraints) to it. Here is how:
- Gurobi
import gurobiby as grbopt_model = grb.Model(name="MIP Model")
- CPLEX
import docplex.mp.model as cpxopt_model = cpx.Model(name="MIP Model")
- PuLp
import pulp as plpopt_model = plp.LpProblem(name="MIP Model")
After this step, we have a Model Object named opt_model. Next, we need to add decision variables. It is standard to store decision variables in python dictionaries (or Pandas Series) where dictionary keys are decision variables indexes and values are decision variable objects. A decision variable is defined with three main properties: its type (continuous, binary or integer), its lower bound (0 by default), and its upper bound (infinity by default). For the above example, we can define decision variables as:
- Gurobi
# if x is Continuous
x_vars ={(i,j):opt_model.addVar(vtype=grb.GRB.CONTINUOUS,
lb=l[i,j],
ub= u[i,j],
name="x_{0}_{1}".format(i,j))
for i in set_I for j in set_J}# if x is Binary
x_vars = {(i,j):opt_model.addVar(vtype=grb.GRB.BINARY,
name="x_{0}_{1}".format(i,j))
for i in set_I for j in set_J}# if x is Integer
x_vars ={(i,j):opt_model.addVar(vtype=grb.GRB.INTEGER,
lb=l[i,j],
ub= u[i,j],
name="x_{0}_{1}".format(i,j))
for i in set_I for j in set_J}
- CPLEX
# if x is Continuous
x_vars =
{(i,j): opt_model.continuous_var(lb=l[i,j], ub= u[i,j],
name="x_{0}_{1}".format(i,j))
for i in set_I for j in set_J}# if x is Binary
x_vars =
{(i,j): opt_model.binary_var(name="x_{0}_{1}".format(i,j))
for i in set_I for j in set_J}# if x is Integer
x_vars =
{(i,j): opt_model.integer_var(lb=l[i,j], ub= u[i,j],
name="x_{0}_{1}".format(i,j))
for i in set_I for j in set_J}
- PuLp
# if x is Continuous
x_vars = {(i,j):
plp.LpVariable(cat=plp.LpContinuous,
lowBound=l[i,j], upBound=u[i,j],
name="x_{0}_{1}".format(i,j))
for i in set_I for j in set_J}# if x is Binary
x_vars = {(i,j):
plp.LpVariable(cat=plp.LpBinary, name="x_{0}_{1}".format(i,j))
for i in set_I for j in set_J}# if x is Integer
x_vars = {(i,j):
plp.LpVariable(cat=plp.LpInteger,
lowBound=l[i,j], upBound= u[i,j],
name="x_{0}_{1}".format(i,j))
for i in set_I for j in set_J}
After setting decision variables and adding them to our model, it is time to set constraints. Any constraint has three parts: left hand side (normally a linear combination of decision variables), right hand side (usually a numeric value) and sense (Less than or equal, Equal or Greater than or equal). To set up any constraints, we need to set each part. Here is how:
- Gurobi
# <= constraints
constraints = {j :
opt_model.addConstr(
lhs=grb.quicksum(a[i,j] * x_vars[i,j] for i in set_I),
sense=grb.GRB.LESS_EQUAL,
rhs=b[j],
name="constraint_{0}".format(j))
for j in set_J}# >= constraints
constraints = {j :
opt_model.addConstr(
lhs=grb.quicksum(a[i,j] *x_vars[i,j] for i in set_I),
sense=grb.GRB.GREATER_EQUAL,
rhs=b[j],
name="constraint_{0}".format(j))
for j in set_J}# == constraints
constraints = {j :
opt_model.addConstr(
lhs=grb.quicksum(a[i,j] * x_vars[i,j] for i in set_I),
sense=grb.GRB.EQUAL,
rhs=b[j],
name="constraint_{0}".format(j))
for j in set_J}
- CPLEX
# <= constraints
constraints = {j :
opt_model.add_constraint(
ct=opt_model.sum(a[i,j] * x_vars[i,j] for i in set_I) <= b[j],
ctname="constraint_{0}".format(j))
for j in set_J}# >= constraints
constraints = {j :
opt_model.add_constraint(
ct=opt_model.sum(a[i,j] * x_vars[i,j] for i in set_I) >= b[j],
ctname="constraint_{0}".format(j))
for j in set_J}# == constraints
constraints = {j :
opt_model.add_constraint(
ct=opt_model.sum(a[i,j] * x_vars[i,j] for i in set_I) == b[j],
ctname="constraint_{0}".format(j))
for j in set_J}
- PuLp
# Less than equal constraints
constraints = {j :
plp.LpConstraint(
e=plp.lpSum(a[i,j] * x_vars[i,j] for i in set_I),
sense=plp.plp.LpConstraintLE,
rhs=b[j],
name="constraint_{0}".format(j))
for j in set_J}# >= constraints
constraints = {j :
plp.LpConstraint(
e=plp.lpSum(a[i,j] * x_vars[i,j] for i in set_I),
sense=plp.LpConstraintGE,
rhs=b[j],
name="constraint_{0}".format(j))
for j in set_J}# == constraints
constraints = {j :
plp.LpConstraint(
e=plp.lpSum(a[i,j] * x_vars[i,j] for i in set_I),
sense=plp.LpConstraintEQ,
rhs=b[j],
name="constraint_{0}".format(j))
for j in set_J}
Next step is defining objective which is a linear expression. Here is how we can define an objective:
- Gurobi
objective = grb.quicksum(x_vars[i,j] * c[i,j]
for i in set_I
for j in set_J)- CPLEX
objective = opt_model.sum(x_vars[i,j] * c[i,j]
for i in set_I
for j in set_J)- PuLp
objective = plp.sum(x_vars[i,j] * c[i,j]
for i in set_I
for j in set_J)Now, we need to add objective to our model:
- Gurobi
# for maximization
opt_model.ModelSense = grb.GRB.MAXIMIZE# for minimization
opt_model.ModelSense = grb.GRB.MINIMIZEopt_model.setObjective(objective)
- CPLEX
# for maximization
opt_model.maximize(objective)# for minimization
opt_model.minimize(objective)
- PuLp
# for maximization
opt_model.sense = plp.LpMaximize# for minimization
opt_model.sense = plp.LpMinimizeopt_model.setObjective(objective)
Finally, we call solver to solve our optimization model. In PuLp, the default solver is CBC, but it can go with other solvers as well. Here is the final step in solving our model:
- Gurobi
opt_model.optimize()- CPLEX
# solving with local cplex
opt_model.solve()# solving with cplex cloud
opt_model.solve(url="your_cplex_cloud_url", key="your_api_key")
- PuLp
# solving with CBC
opt_model.solve()
# solving with Glpk
opt_model.solve(solver = GLPK_CMD())Now we are done. We just need to get results and post process them. I’ve found Pandas package as a good data processing library. If the problem is solved to optimality, we can get and process results as:
import pandas as pdopt_df = pd.DataFrame.from_dict(x_vars, orient="index",
columns = ["variable_object"])opt_df.index =
pd.MultiIndex.from_tuples(opt_df.index,
names=["column_i", "column_j"])opt_df.reset_index(inplace=True)# Gurobi
opt_df["solution_value"] =
opt_df["variable_object"].apply(lambda item: item.X)# CPLEX
opt_df["solution_value"] =
opt_df["variable_object"].apply(lambda item: item.solution_value)# PuLp
opt_df["solution_value"] =
opt_df["variable_object"].apply(lambda item: item.varValue)opt_df.drop(columns=["variable_object"], inplace=True)
opt_df.to_csv("./optimization_solution.csv")
“opt_df” is a pandas dataframe that holds optimal values of each decision variable. We can also save the result in a csv file as shown above.
I hope you’ve found this blog useful. We just covered high level modeling in Python, however, all of above packages contain useful functions and data structures that should be considered when you write production ready code. For example, in Gurobi you can add a set of variables at once using opt_model.addVars(), in CPLEX you can do that with opt_model.continuous_var_dict(), opt_model.binary_var_dict(), or opt_model.integer_var_dict(). Similarly, in PuLp that can be achieved with plp.LpVariable.dicts().
