Linear Programming (LP)

📖 Optimization API Reference

The library offers a solver for a wide range of linear optimization problems, including continuous, integer and mixed-integer problems. The LinearProgrammingProblem class uses the linprog module from SciPy (which eventually calls the HIGHS solver) for the solution, combined with the power of SymPy for being able to pose problems directly in symbolic form. The class handles all the preprocessing required to bring the problem into the form required by SciPy.

First, import some stuff. Note that the notebook requires matplotlib to be installed.

from sigmaepsilon.math.function import Function, Relation
from sigmaepsilon.math.optimize import LinearProgrammingProblem as LPP
import sympy as sy

One of the great features of the optimization module is that it handles symbolic functions pretty well. Problems can be defined using SymPy expressions, or simple strings. To solve a problem, call the solve method, which returns an instance of scipy.optimize.OptimizeResult.

x1, x2, x3, x4 = variables = sy.symbols(['x1', 'x2', 'x3', 'x4'])
obj = Function(3*x1 + 9*x3 + x2 + x4, variables=variables)
eq1 = Relation(x1 + 2*x3 + x4 - 4, variables=variables)
eq2 = Relation(x2 + x3 - x4 - 2, variables=variables)
bounds = [(0, None), (0, None), (0, None), (0, None)]
lpp = LPP(obj, [eq1, eq2], variables=variables, bounds=bounds)
solution = lpp.solve()
solution.x

solution

type(solution)

If the bounds for the variables are all the same, you can afford to input a single bound:

x1, x2, x3, x4 = variables = sy.symbols(['x1', 'x2', 'x3', 'x4'])
obj = Function(3*x1 + 9*x3 + x2 + x4, variables=variables)
eq1 = Relation(x1 + 2*x3 + x4 - 4, variables=variables)
eq2 = Relation(x2 + x3 - x4 - 2, variables=variables)
bounds = (0, None)
lpp = LPP(obj, [eq1, eq2], variables=variables, bounds=bounds)
lpp.solve().x

Furthermore, if all the bounds are in the form \(\geq\), you can omit the arguments for the parameter bounds, since it is the default value.

variables = sy.symbols(['x1', 'x2', 'x3', 'x4'])
obj = Function("3*x1 + 9*x3 + x2 + x4", variables=variables)
eq1 = Relation("x1 + 2*x3 + x4 = 4", variables=variables)
eq2 = Relation("x2 + x3 - x4 = 2", variables=variables)
lpp = LPP(obj, [eq1, eq2], variables=variables)
lpp.solve().x

2d example with matplotlib

Consider the following problem:

\[\begin{split} \begin{aligned} & \text{minimize}& \quad x_1 + x_2 &\\ & \text{subject to} & & \\ & & x_1 &\geq 1, \\ & & x_2 &\geq 1, \\ & & x_1 + x_2 &\leq 4, \\ & & x_i &\geq 0, \qquad i=1, 2. \end{aligned} \end{split}\]

Setting up the problem is as follows.

variables = sy.symbols(['x1', 'x2'])
f = Function("x1 + x2", variables=variables)
ieq1 = Relation("x1 >= 1", variables=variables)
ieq2 = Relation("x2 >= 1", variables=variables)
ieq3 = Relation("x1 + x2 <= 4", variables=variables)
bounds = [(0, None), (0, None)]  # x1 >= 0, x2 >= 0
lpp = LPP(f, [ieq1, ieq2, ieq3], variables=variables, bounds=bounds)
x = lpp.solve().x
print("The optimal solution is x1 = {0}, x2 = {1}.".format(x[0], x[1]))

Now plotting with Matplotlib, with the feasible side of inequalities visualized by hatching, the ticks being on the infeasible side.

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import patheffects

fig, ax = plt.subplots(figsize=(6, 6))

nx = 100
ny = 100
xvec = np.linspace(0.0, 4.0, nx)
yvec = np.linspace(0.0, 4.0, ny)
x1, x2 = np.meshgrid(xvec, yvec)

obj = x1 + x2
g1 = x1 - 1
g2 = x2 - 1
g3 = x1 + x2 - 4

cntr = ax.contour(x1, x2, obj, colors='black')
ax.clabel(cntr, fmt="%2.1f", use_clabeltext=True)

cg1 = ax.contour(x1, x2, g1, [0], colors='sandybrown')
cg1.set(path_effects=[patheffects.withTickedStroke(angle=-135)])

cg2 = ax.contour(x1, x2, g2, [0], colors='orangered')
cg2.set(path_effects=[patheffects.withTickedStroke(angle=-60)])

cg3 = ax.contour(x1, x2, g3, [0], colors='mediumblue')
cg3.set(path_effects=[patheffects.withTickedStroke()])

ax.plot(x[0], x[1], 'ms', markersize=15)

ax.set_xlim(0, 4)
ax.set_ylim(0, 4)

plt.show()

Linear Mixed-Integer Programs

If some of the design variables are constrained to integers, you can use the integrality parameter. As an example, consider the following problem:

\[\begin{split} \begin{aligned} & \text{minimize}& \quad 3 x_2 + 2 x_3 &\\ & \text{subject to} & & \\ & & 2 x_1 + 2 x_2 - 4 x_3 &= 5, \\ & & x_i &\geq 0, \qquad i=1, 2, 3, \\ & & x_i &\in \mathbb{Z}, \qquad i=1, 3. \end{aligned} \end{split}\]

variables = x1, x2, x3 = sy.symbols(["x1", "x2", "x3"])
obj = Function(3 * x2 + 2 * x3, variables=variables)
eq = Relation(2 * x1 + 2 * x2 - 4 * x3 - 5, op="=", variables=variables)
bounds = (0, None)
integrality = [1, 0, 1]
lpp = LPP(obj, [eq], variables=variables, bounds=bounds, integrality=integrality)
solution = lpp.solve()
solution.x

Another way to express the same integrality constraints is using assumptions when creating the SymPy variables of the problem.

x1, x3 = sy.symbols(["x1", "x3"], integer=True)
x2 = sy.symbols("x2")
variables = x1, x2, x3
obj = Function(3 * x2 + 2 * x3, variables=variables)
eq = Relation(2 * x1 + 2 * x2 - 4 * x3 - 5, op="=", variables=variables)
bounds = (0, None)
lpp = LPP(obj, [eq], variables=variables, bounds=bounds)
solution = lpp.solve()
solution.x

The solver in SciPy is capable of understanding more complex integrality constraints, colsult their API Reference for more details here. If you want to impose constraints like these, you must use the integrality parameter.