Linear programming with Python and Julia
“True optimization is the revolutionary contribution of modern research to decision processes.” — George Dantzig
I was intrigued by the concept of optimization when I attended the course Operations Research (OR) during my undergraduate studies in Mechanical Engineering half a decade ago. The main reason this course was fascinating to me was that it dealt with solving real-world problems such as optimizing the workflow in a factory, supply chain management, scheduling flights in an airport, travelling salesman problem, etc. Operations Research deals with how to make decisions efficiently through the use of different mathematical techniques or algorithms. In a real-world setting, this could mean maximizing (profit, yield) or minimizing (losses, risks) the given expression while satisfying the constraints such as costs, time, and resource allocation.
There are several applications of OR in domains such as energy system optimization, supply chain management, logistics and inventory management, routing and pathfinding problems, etc. [1]. When I was an undergrad student, solving optimization problems in the OR course manually with hand and a calculator was a daunting task. A human error in a single step meant all the following steps are wrong, and one has to redo the entire procedure from the start. Thanks to the advancement in the programming techniques, now there are open-source tools such as Google OR-Tools, and different packages in Python, and Julia, which facilitate solving the optimization problem in a matter of seconds. One just has to define the problem in a framework understandable to the given package.
Besides the Google OR-Tools, some open-source packages available for solving optimization problems in Python are scipy.optimize, PuLP and Pyomo. In Julia, there is a similar package embedded within the language called JuMP. In this post, I am going to solve a simple linear optimization problem first using the Pyomo package in Python, replicate it in Julia using JuMP, and share my experience. The code used in this post is available in this GitHub repository. Let’s get started:
Problem Statement
For this post, I consider a problem of determining an optimal product mix in a factory based on the processing time of the products in different machines, availability of the machines, and the unit profit of each product.
A company manufactures two products: X and Y. To manufacture each product, it has to go through three machines: A, B, and C. Manufacturing X require 3 hours in machine A, 9 hours in machine B, and 2 hours in machine C. Similarly, manufacturing product Y require 2, 4, and 10 hours in machines A, B, and C respectively. The availability of each of the machines A, B, and C during a manufacturing period are 66, 180, and 200 hours respectively. The profit per product X is USD 90 and that per product Y is USD 75. How many units of X and Y should be produced during a production period to maximize profit?
Assuming x and y as the units of X and Y to be produced respectively, x and y are our decision variables. The problem could be represented in algebraic form as follows:
Profit = 90x+75yObjective: maximize 90x+75y subject to:3x+2y≤669x+4y≤1802x+10y≤200x, y≥0
Solver
Solvers embed powerful algorithms to solve optimization problems and help improve decision-making around planning, allocation, and scheduling of the resources under constraints to meet the given objectives. Based on the problem class, one needs to select the suitable solver to solve the optimization problem. The following table shows the examples of some solvers which have both Python and Julia wrappers available.
It is important to note that not all solvers are open-access. For e.g. IPOPT and GLPK are open-access solvers, while CPLEX and GUROBI require commercial licenses.
The problem used in this post is an example of linear programming since both the objective and constraints are linear. In this case, we can use open-access solvers either GLPK or IPOPT to solve the problem. If the problem was non-linear (for e.g. with quadratic constraints), then only the IPOPT solver would have been usable in that case, and not GLPK based on the problem class.
Solving in Python with Pyomo package
Having some basic hands-on experience of working with the Pyomo package, I find it pretty intuitive to work with this package in Python. Pyomo allows two strategies for declaring models. When a mathematical model is defined by using symbols that rely on unspecified parameter values (e.g. ax+by=c), it is referred to as an Abstract Model. Passing data through an Abstract Model creates an instance of the model (also referred to as a Concrete Model; e.g. 2x+3y = 5) which could be passed through the solver, as depicted in the flowchart below.
The declaration of modelling components such as variables is quite straightforward in Pyomo. The objective function and constraints are declared in the form of expressions. While declaring an objective, it is important to provide the sense: whether to minimize or maximize the given expression.
After declaring the necessary components, the model is passed through a solver (here GLPK) for solving the linear optimization problem. When the solver status is “ok” and the termination condition “optimal”, it means that the optimization problem has been solved successfully.
For the given problem, the model determines that 10 items of X and 18 items of Y are optimal to maximize the profit of the company, which is $2250 in this case. This problem could also be solved graphically by plotting the model variables and constraints in the graph. I have plotted them using matplotlib, with the code in the following gist:
As shown in the plot above, the red marker at (10, 18) represents the optimal product strategy. I define the space where all the red, green, and blue shades overlap each other as “Feasibility Space”, wherein, the model meets all the constraints.
Solving in Julia with JuMP
Launched in 2012, Julia is an open-source language whose community is growing in recent years. Thanks to the strong performance and flexibility of the language, Julia has been deemed apt for computationally intensive tasks [3]. At the time of writing, there are over 4000 packages in the General registry of the language.
As a newbie with the Julia language, I wanted to give it a try by replicating the same linear problem. In this section, I am going to describe how I did it. First I opened Julia’s command-line REPL (read-eval-print-loop) and added all the required packages using Pkg, which is the built-in package manager of Julia.
using Pkg #In-built package manager of JuliaPkg.add(“LinearAlgebra”) #Analogous to numpy package in PythonPkg.add(“Plots”) #Add Plots.jl framework for plottingPkg.add(“JuMP”) #Add mathematical optimization packagePkg.add(“GLPK”) #Add solverPkg.add(“IJulia”) #Backend for Jupyter interactive environmentPkg.add(“PyPlot”) #Julia interface to matplotlib.pyplot in Python
Installing these packages was pretty convenient and each package was installed in my system within a matter of seconds. After installing the packages, they are imported using the syntax: using <package> in Julia which is similar to import <package>in Python.
After the packages are called, the model is declared and the optimizer is set from GLPK. As shown in the gist above, declaring variables, constraints and objectives are much easier in Julia because one can declare directly using the real algebraic expression, which makes the script shorter and more comprehensible.
As shown in the screenshot above, I get the same results with JuMP as with Pyomo. Next, I plot the linear problem graphically in Julia using PyPlot.
The PyPlot module in Julia provides a Julia interface to the matplotlib plotting library for Python, and specifically to the matplotlib.pyplot module. Therefore, it is a pre-requisite to have the matplotlib library installed in one’s system to use PyPlot [4]. I noticed that the parameters to be passed using PyPlot in Julia are basically the same as for matplotlib.pyplot in Python, and that only the nomenclature is different.
Conclusion
In this post, I demonstrated an example of solving a simple linear optimization problem first using the Pyomo in Python and then using JuMP in Julia. These open-source packages as well as the solvers are a real boon to solve complicated optimization problems in the domain of Operations Research, which otherwise would cost a lot of time and resources. Moreover, they also facilitate improving speed, precision, and flexibility in solving optimization problems (for e.g. including sensitivity analysis).
Although a newbie with the Julia language, with my experience in Python, I found it relatively easy to understand the syntax, which is user-friendly, and formulate the given problem in Julia. Julia has been deemed to combine the interactivity and syntax of scripting languages such as Python, and the speed of compiled languages such as C [3]. It might be slow when a function is keyed in initially, however, the subsequent runs are supposedly faster.
References:
[1] Vasegaard, A. (2021). Why Operations Research is awesome — An Introduction?
[2] Hart et al., (2021). Pyomo - Optimization Modeling in Python.
[3] Perkel, J.M. (2021). Julia: come for the syntax, stay for the speed
[4] Johnson, S.G. (2021). The PyPlot module for Julia.
