Toolchain Tuesday No. 5
TL;DR: Part of a series of posts about tools, services, and packages that I use in day-to-day operations to boost efficiency and free up time for the things that really matter. Use at your own risk - happy to answer questions. For the full, continuously expanding list so far see here.
This is the fifth installment of a series of posts; the full list is expanding over time. This time around will be about modeling languages for optimization problems.
Python Libraries
CVXOPT
Low-level Python interface for convex optimization.
Learning curve: ⭐️⭐️⭐️⭐️
Usefulness: ⭐️⭐️⭐️
Site: https://cvxopt.org
CVXOPT is basically a Python interface to various optimization solvers, providing an intermediate, relatively low-level, matrix-based interface. This is in contrast to some of the modeling languages from below that provide a higher level of abstraction; basically these modeling languages generate the matrix structure by transcribing the statements of the modeling language. Nonetheless, CVXOPT is a great tool to solve optimization problems in Python.
From https://cvxopt.org:
CVXOPT is a free software package for convex optimization based on the Python programming language. It can be used with the interactive Python interpreter, on the command line by executing Python scripts, or integrated in other software via Python extension modules. Its main purpose is to make the development of software for convex optimization applications straightforward by building on Python’s extensive standard library and on the strengths of Python as a high-level programming language.
Here is sample code from https://cvxopt.org to give you an idea:
# Risk-return trade-off.
from math import sqrt
from cvxopt import matrix
from cvxopt.blas import dot
from cvxopt.solvers import qp, options
n = 4
S = matrix( [[ 4e-2, 6e-3, -4e-3, 0.0 ],
[ 6e-3, 1e-2, 0.0, 0.0 ],
[-4e-3, 0.0, 2.5e-3, 0.0 ],
[ 0.0, 0.0, 0.0, 0.0 ]] )
pbar = matrix([.12, .10, .07, .03])
G = matrix(0.0, (n,n))
G[::n+1] = -1.0
h = matrix(0.0, (n,1))
A = matrix(1.0, (1,n))
b = matrix(1.0)
N = 100
mus = [ 10**(5.0*t/N-1.0) for t in range(N) ]
options['show_progress'] = False
xs = [ qp(mu*S, -pbar, G, h, A, b)['x'] for mu in mus ]
returns = [ dot(pbar,x) for x in xs ]
risks = [ sqrt(dot(x, S*x)) for x in xs ]
Pyomo
Pyomo is a python-based open-source optimization modeling language supporting a wide range of optimization paradigms and solvers.
Learning curve: ⭐️⭐️⭐️
Usefulness: ⭐️⭐️⭐️⭐️
Site: http://www.pyomo.org/
Pyomo is a Python-based open-source optimization modeling language. It supports a variety of different optimization paradigms and integrates with a wide range of solvers including BARON, CBC, CPLEX, Gurobi, and glpsol; check the Pyomo Manual. Another great resource with examples is the Pyomo Cookbook.
What sets Pyomo apart from MathProg and CVXOPT is that it is a relatively high-level modeling language (compared to CVXOPT) while being written in Python (compared to MathProg) allowing for easy integration with a plethora of other packages.
From http://www.pyomo.org/:
A core capability of Pyomo is modeling structured optimization applications. Pyomo can be used to define general symbolic problems, create specific problem instances, and solve these instances using commercial and open-source solvers. Pyomo’s modeling objects are embedded within a full-featured high-level programming language providing a rich set of supporting libraries, which distinguishes Pyomo from other algebraic modeling languages like AMPL, AIMMS and GAMS.
Supported problem types include:
- Linear programming
- Quadratic programming
- Nonlinear programming
- Mixed-integer linear programming
- Mixed-integer quadratic programming
- Mixed-integer nonlinear programming
- Stochastic programming
- Generalized disjunctive programming
- Differential algebraic equations
- Bilevel programming
- Mathematical programs with equilibrium constraints
Here is an example of a model:
from pyomo.environ import *
model = ConcreteModel()
# declare decision variables
model.y = Var(domain=NonNegativeReals)
# declare objective
model.profit = Objective(
expr = 30*model.y,
sense = maximize)
# declare constraints
model.laborA = Constraint(expr = model.y <= 80)
model.laborB = Constraint(expr = model.y <= 100)
# solve
SolverFactory('glpk').solve(model).write()
MathProg
MathProg (aka GMPL) is a modeling language for Mixed-Integer Linear Programs.
Learning curve: ⭐️⭐️
Usefulness: ⭐️⭐️⭐️
Site: https://www.gnu.org/software/glpk/
MathProg, also known as GMPL is a modeling language for Mixed-Integer Linear Programs (MILP). It is included with glpk, the GNU Linear Programming Kit and it supports reading and writing from data sources (databases via ODBC or JDBC) and, apart from glpsol which is glpk’s own MILP solver it supports various other solvers such as CPLEX, Gurobi, or SCIP through LP-files. From the [GLPK wikibook], which is also a great resource for GMPL:
GNU MathProg is a high-level language for creating mathematical programming models. MathProg is specific to GLPK, but resembles a subset of AMPL. MathProg can also be referred to as GMPL (GNU Mathematical Programming Language), the two terms being interchangeable.
MathProg is in particular great for fast prototyping. Unfortunately, it does not directly support other IP solvers through its built-in interface but requires to go the route via LP-files. As a consequence the relatively powerful output post-processing of MathProg cannot be used in that case and glpk’s own solver glpsol that integrates with MathProg natively is only able to handle smaller to midsize problems. This limits the use case (therefore the ⭐️⭐️⭐️-rating on usefulness). Probably the natural course of things is to ‘graduate’ to Pyomo over time. Despite all of this, the actual MathProg language very useful. I have used it very often for the actual modeling task, to separate (supporting) code from model. I then generate an LP file from that model and its datasources that I then solve with e.g., Gurobi. Finally, I parse the output back in, mostly using Python; there are other Python packages for glpsol that handle this parsing for you if you do not want to implement it yourself (relatively easy though). Pyomo for example also goes through the LP-file route to interface with glpsol.
See here on how to obtain glpk. To give you an idea of the syntax, check out this example that also includes some solution post-processing (similar to what e.g., OPL can do):
# A TRANSPORTATION PROBLEM
#
# This problem finds a least cost shipping schedule that meets
# requirements at markets and supplies at factories.
#
# References:
# Dantzig G B, "Linear Programming and Extensions."
# Princeton University Press, Princeton, New Jersey, 1963,
# Chapter 3-3.
set I;
/* canning plants */
set J;
/* markets */
param a{i in I};
/* capacity of plant i in cases */
param b{j in J};
/* demand at market j in cases */
param d{i in I, j in J};
/* distance in thousands of miles */
param f;
/* freight in dollars per case per thousand miles */
param c{i in I, j in J} := f * d[i,j] / 1000;
/* transport cost in thousands of dollars per case */
var x{i in I, j in J} >= 0;
/* shipment quantities in cases */
minimize cost: sum{i in I, j in J} c[i,j] * x[i,j];
/* total transportation costs in thousands of dollars */
s.t. supply{i in I}: sum{j in J} x[i,j] <= a[i];
/* observe supply limit at plant i */
s.t. demand{j in J}: sum{i in I} x[i,j] >= b[j];
/* satisfy demand at market j */
solve;
# Report / Result Section (Optional)
printf '#################################\n';
printf 'Transportation Problem / LP Model Result\n';
printf '\n';
printf 'Minimum Cost = %.2f\n', cost;
printf '\n';
printf '\n';
printf 'Variables (i.e. shipment quantities in cases ) \n';
printf 'Shipment quantities in cases\n';
printf 'Canning Plants Markets Solution (Cases) \n';
printf{i in I, j in J}:'%14s %10s %11s\n',i,j, x[i,j];
printf '\n';
printf 'Constraints\n';
printf '\n';
printf 'Observe supply limit at plant i\n';
printf 'Canning Plants Solution Sign Required\n';
for {i in I} {
printf '%14s %10.2f <= %.3f\n', i, sum {j in J} x[i,j], a[i];
}
printf '\n';
printf 'Satisfy demand at market j\n';
printf 'Market Solution Sign Required\n';
for {j in J} {
printf '%5s %10.2f >= %.3f\n', j, sum {i in I} x[i,j], b[j];
}
data;
set I := Seattle San-Diego;
set J := New-York Chicago Topeka;
param a := Seattle 350
San-Diego 600;
param b := New-York 325
Chicago 300
Topeka 275;
param d : New-York Chicago Topeka :=
Seattle 2.5 1.7 1.8
San-Diego 2.5 1.8 1.4 ;
param f := 90;
end;