# 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;
```