Mathematical programming with Pyomo

Author: Jérémie Decock (2019)

What is mathematical programming

Mathematical programming aims at finding the solutions (vectors) which minimize (or maximize) an objective function, the feasible solution space being defined by a set of constraints (equations and inequations).

Very basic example of linear program (i.e. a mathematical optimization problem with a linear objective function and linear constraints):

$$ \begin{align} \min_{x_1,x_2} & \quad -x_1 + 4 x_2 \\ \text{s.t.} & \quad -3 x_1 + x_2 \leq 6 \\ & \quad -x_1 - 2 x_2 \geq -4 \\ & \quad x_1 \leq 10 \end{align} $$

The first line is the objective function and the following inequations are the constraints that define the solution space.

The optimal solution of this problem is $\boldsymbol{x}^* = \begin{pmatrix}10 \\ -3\end{pmatrix}$ with a cost equals to -22.


In [ ]:
%matplotlib inline

import matplotlib.pyplot as plt
import seaborn as sns

In [ ]:
sns.set_context('talk')

In [ ]:
def intersect(A, b):
    # Ax = b   <=>   x = A^{-1} b
    inv_A = np.linalg.inv(A)
    return np.dot(inv_A, b)

fig, ax = plt.subplots(figsize=(16, 10))

a1 = [-3,  1]
a2 = [-1, -2]
a3 = [ 1,  0]

a4 = [ 1,  0]
a5 = [ 1,  0]
a6 = [ 0,  1]
a7 = [ 0,  1]

b1 =  6
b2 = -4
b3 = 10

b4 = -3
b5 = 12
b6 = -10
b7 = 50

A = np.array([a1, a2])
b = np.array([[b1], [b2]])
p12 = intersect(A, b)

A = np.array([a1, a3])
b = np.array([[b1], [b3]])
p13 = intersect(A, b)

A = np.array([a2, a3])
b = np.array([[b2], [b3]])
p23 = intersect(A, b)

A = np.array([a1, a4])
b = np.array([[b1], [b4]])
p14 = intersect(A, b)

A = np.array([a1, a5])
b = np.array([[b1], [b5]])
p15 = intersect(A, b)

A = np.array([a2, a4])
b = np.array([[b2], [b4]])
p24 = intersect(A, b)

A = np.array([a2, a5])
b = np.array([[b2], [b5]])
p25 = intersect(A, b)

A = np.array([a3, a6])
b = np.array([[b3], [b6]])
p36 = intersect(A, b)

A = np.array([a3, a7])
b = np.array([[b3], [b7]])
p37 = intersect(A, b)

ax.plot([p14[0], p15[0]], [p14[1], p15[1]], ":k")
ax.plot([p24[0], p25[0]], [p24[1], p25[1]], ":k")
ax.plot([p36[0], p37[0]], [p36[1], p37[1]], ":k")

ax.plot([p12[0], p13[0]], [p12[1], p13[1]], "o-k")
ax.plot([p12[0], p23[0]], [p12[1], p23[1]], "o-k")
ax.plot([p13[0], p23[0]], [p13[1], p23[1]], "o-k")

plt.fill_between([float(p12[0]), b3], np.array([p12[1], p13[1]]).flatten(), np.array([p12[1], p23[1]]).flatten(), facecolor='red', alpha=0.5)

ax.plot([p23[0]], [p23[1]], "or")

# CONTOURS ###

delta = 0.1
x1 = np.arange(b4, b5, delta)
x2 = np.arange(b6, b7, delta)
X1, X2 = np.meshgrid(x1, x2)

Z = -X1 + 4*X2

cs = plt.contour(x1, x2, Z,
                 alpha=0.5,
                 label="TC")
ax.clabel(cs, inline=False, fontsize=12)

ax.set_xlabel("x1")
ax.set_ylabel("x2");

Here, the solution space defined by the three constraints is depicted by the red area. Colored lines are isocontours of the objective function. The red dot is the optimal solution (the point in the red area that minimize the objective function).

Note: a linear program can be fully defined with one matrix and two vectors (respectively $\color{\orange}{A}$, $\color{\green}{b}$ and $\color{\red}{c}$ in the following example):

$$ \begin{align} \min_{\boldsymbol{x}} & \quad \color{\red}{\boldsymbol{c}}^{\top} \boldsymbol{x} \\ \text{s.t.} & \quad \color{\orange}{\boldsymbol{A}} \boldsymbol{x} \leq \color{\green}{\boldsymbol{b}} \end{align} $$$$ \color{\red}{ \boldsymbol{c} = \begin{pmatrix} c_1 \\ c_2 \\ \vdots \end{pmatrix} } \quad \color{\orange}{ \boldsymbol{A} = \begin{pmatrix} A_{1,1} & A_{1,2} & \cdots \\ A_{2,1} & A_{2,2} & \cdots\\ \vdots & \vdots & \ddots \end{pmatrix} } \quad \color{\green}{ \boldsymbol{b} = \begin{pmatrix} b_1 \\ b_2 \\ \vdots \end{pmatrix} } $$

Kind of problems solved in mathematical programming:

  • Linear programming (LP): linear objective function and constraints, continuous solution space
  • Quadratic programming (QP): quadratic objective function and constraints, continuous solution space
  • Nonlinear programming: nonlinear objective function and constraints, continuous solution space
  • Mixed-integer linear programming (MILP): linear objective function and constraints, locally discrete or continuous solution space
  • Mixed-integer quadratic programming: quadratic objective function and constraints, locally discrete or continuous solution space
  • Mixed-integer nonlinear programming: nonlinear objective function and constraints, locally discrete or continuous solution space
  • Stochastic programming
  • Generalized disjunctive programming
  • Differential algebraic equations
  • Bilevel programming
  • Mathematical programs with equilibrium constraints

The two main components involved in mathematical programming are:

  • Algebraic modeling tools and Algebraic modeling language: languages or libraries used to easily define optimization problems (optional but very often used)
  • Solvers: they implement algorithms used to find the optimal solution to optimization problems

Problems are defined with Algebraic modeling tools (in a specific Algebraic modeling language). Then the Algebraic modeling tools submit this problem to a selected solver and get back the optimal solution(s). Algebraic modeling tools are usually compatible with several solvers (thanks to common formats or libraries).

Algebraic modeling language (AML)

See https://en.wikipedia.org/wiki/Algebraic_modeling_language

Examples of commercial AMLs:

Examples of open source AMLs:

Solvers

Examples of commercial solvers:

Examples of open source solvers:

What is Pyomo

Pyomo is a Python-based open-source optimization modeling language for mathematical programming.

See http://www.pyomo.org/about

Getting started with Pyomo

TODO Installation instructions (Pyomo + open source solvers) can be find on this page or this one.

Installation instructions (Pyomo + open source solvers) can be find on this page.

As an alternative, can be used to work directly in the cloud and avoid local installation (check for instance this page).

Pyomo's official documentation is there and very nice tutorials are available there.