freight boat transporting containers on ocean - linear programming - simplex algorithm

Playing with Mathematics and Programming – The Simplex Algorithm for Linear Programming – Two Challenges

Share.

We will give you two challenges in this article. Before we do so, we will tell you what the simplex algorithm for solving linear programming problems is. We will first show you a few simple cases.

What is the simplex algorithm for linear programming?

One-dimensional case

Let $f$ be a linear function defined as $cx$, where $c$ is a constant real number and $x$ the real variable. We impose some constraints on $x$. The first one is the requirement that $x$ be nonnegative, and the second constraint is that $ax$ be at most some fixed real number $b$, where $a$ is also a fixed real number. The function $f$ is called the objective function, and the collection of all real numbers $x$ satisfying these two constraints is called the feasible region.

A linear programming problem is to find a number in the feasible region that maximizes the objective function, meaning to find a number $x^{*}$ in the feasible region such that for any other number $y$ in that region $f\left(y\right)$ is at most $f\left(x^{*}\right)$.

The simplex algorithm find that number $x^{*}$, if it exists.

one dimension case - linear programming - simplex algorithm
Feasible region in 1 dimension

Two-dimensional case

In a two-dimensional case, the constants and variables will be vectors with two components. Let $f$ be a linear function defined as $\mathbf{c}^{\top}\mathbf{x}$, where $\mathbf{c}=\left[\begin{array}{c} c_{1}\\ c_{2} \end{array}\right],\mathbf{x}=\left[\begin{array}{c} x_{1}\\ x_{2}\end{array}\right]$, and $\bullet^{\top}$ is the transpose operator that turn columns of a matrix into rows and rows into columns. The vector $\mathbf{c}$ is constant, $\mathbf{x}$ is the variable with real components, and $\mathbf{c}^{\top}\mathbf{x}$ is a matrix multiplication. The first constraint for $\mathbf{x}$ is that each of its components be nonnegative, meaning $x_{1}$ and $x_{2}$ be nonnegative. The second constraint for $\mathbf{x}$ is that a finite number of linear inequalities hold, with unknowns $x_{1}$ and $x_{2}$. This system of linear inequalities can be written as

$\begin{array}{ccc} a_{11}x_{1}+a_{12}x_{2} & \leq & b_{1}\\ \vdots & \vdots & \vdots\\ a_{n1}x_{1}+a_{n2}x_{2} & \leq & b_{n} \end{array},$

where the coefficients $a_{ij}$ of the variables $x_{1},x_{2}$ are real constants and $b_{1},\ldots,b_{n}$ are fixed real numbers.

We can write this system of linear inequalities in a more compact way by using matrices, where $A=\left[\begin{array}{cc} a_{11} & a_{12}\\\vdots & \vdots\\ a_{n1} & a_{n2} \end{array}\right]$ is the matrix of the coefficients and $\mathbf{b}=\left[\begin{array}{c} b_{1}\\ \vdots\\ b_{n} \end{array}\right]$. The system of linear inequalities then becomes

$A\mathbf{x}\leq\mathbf{b},$

where $A\mathbf{x}$ is a product of matrices.

The feasible region is the region bounded by the lines given by $a_{i1}x_{1}+a_{i2}x_{2}=b_{i}$ for $i=1,\ldots,n$, and the feasible region contains all those lines. Geometrically, the feasible region is a $2$-dimensional polytope with the number of sides equal to the number of linear inequalities. For example, the feasible region below is a $2$-dimensional polytope with $4$ sides.

feasible region - four linear constraints - simplex algorithm
Feasible region in 2 dimensions

The simplex algorithm finds a vector $\mathbf{x}^{*}$ that maximizes the objective function.

Three-dimensional case

In the three-dimension case, the setting is similar, that is, the goal is to find a vector that maximizes the objective function under the constraints the possible vectors are within some region and their components must be nonnegative. The difference is that all vectors have $3$ components and the feasible region is bounded by half-planes.

For example, the feasible region below is bounded by $5$ half-planes, so it is a $3$-dimensional polytope with $5$ sides.

feasible region - 5 half planes constraints - simplex algorithm

When viewed along the $x$-axis, the three half-planes can be seen:

feasible region - 5 half plane constraints - x axis view - simplex algorithm
The 3 intersecting half-planes when viewed along the x-axis

When viewed along the $z$-axis, the other two half-planes can be seen:

feasible region - 5 half plane constraints - z axis view - simplex algorithm
The 2 other parallel half-planes when viewed along the z-axis

Higher-dimensional case

In the $m$-th dimensional real space, the vectors have $m$ components. If there are $n$ linear inequalities, the matrix $A$ has $n$ rows and $m$ columns, and the feasible region is bounded by $m$-th dimensional half-planes and is an $m$-th dimensional polytope with $n$ sides.

The simplex algorithm works the same as before, that is, it finds a vector in the feasible region that maximizes the objective function.

Computer implementation of the simplex algorithm

There are many computer implementations of the simplex algorithm in Python. Some of them come as packages you can install to your computer [1]. The Python library SciPy implements the simplex algorithm to solve linear programming problems [2].

Microsoft’s Excel has the Solver add-in that uses the simplex algorithm to solve linear programming problems [3,4]. The simplex algorithm has also been implemented in C++ [5], and MathWorks also implements the algorithm for its software [6].

Linear programming in applications in the world

Many problems in the world can be expressed as linear programming problems. These problems can come from several areas, such as agriculture [7], airport surface traffic [8], supply chain management [9], and renewable energy [10].

The Challenges

Before we tell you what the first challenge is, we will tell you what a group action and a rotation group are.

Group action

A group is an abstraction of the addition operation on the integers; it is a set, along with a binary operation, that has an identity element, whose each element has an inverse, and where the order of operation for any three elements is irrelevant.

A subgroup of a group is a subset of the group that contains the identity element and is closed under the same binary operation as the group.

A group action on a set captures what each element of the group does to each element of the set.

An example of a group action is the group of the integers with the usual addition acting on each element on the real line, where each integer $n$ moves to the right each real number by $n$ steps if $n$ is positive, moves to the left each real number by $n$ steps if $n$ is negative, and keeps each real number unchanged if $n=0$.

group action - integers no reals - abstract algebra

Formally, a group $G$ acts on a set $X$ if there is a function $F$ that takes pairs of the form $\left(g,x\right)$ where $g$ is an element of $G$ and $x$ is an element of $X$ and outputs elements in $X$ such that

  1. $F\left(e,x\right)=x$ for all $x$ in $X$, where $e$ is the identity element of $G$, and
  2. $F\left(g,F\left(g’,x\right)\right)=F\left(gg’,x\right)$, where $g,g’$ are elements of $G$.}

The first property means the identity element of the group does nothing to each element of the set. The second property says for any two elements of the group successively doing something on each element of the set is the same as the product of these two elements of the group to do that thing to that element of the set. The “product” is the operation of the group. That function $F$ is the group action.

Rotation group – 2-dimensional case

Imagine you have a vector in the $2$-dimensional real space that you rotate about the origin by some angle $\theta$.

rotated vector in 2d real plane

If you then rotate the result by $-\theta$, you return to the original vector. That rotation by $-\theta$ can be seen as the inverse of the rotation by $\theta$ (existence of inverses).

You could also rotate that vector by the angle $0$, meaning you keep that vector fixed; that would be your identity rotation (existence of an identity).

Also, for any three successive rotations of that vector, the order of rotations is irrelevant (associativity). When you perform two successive rotations of the vector, you can think of it as composing two functions, where the one function is the first rotation and the other function is the second rotation.

The collection of those rotations along with composition can form a group, and this group is an example of a rotation group.

In the $2$-dimensional real space, the elements of the rotation group are square matrices with determinant equal to $1$ [7]. For example, that rotation mentioned above is this matrix with $2$ rows and $2$ columns:

$\left[\begin{array}{cc} \cos\theta & \sin\theta\\ -\sin\theta & \cos\theta\end{array}\right].$

Rotation group – higher dimensional case

In the $n$-th dimensional real space, each element of the rotation group can be thought of as “rotating” each point of the real space about some reference, and a point of the real space is a vector with $n$ entries. Formally, an element of the rotation group is a matrix with $n$ rows and $n$ columns with determinant equal to $1$ [7].

Challenge #1

We are now ready to tell you what the first challenge is.

The setting is the $n$-th dimensional real space. Let $\mathbf{c}=\left[\begin{array}{c} c_{1}\\\vdots\\c_{n}\end{array}\right]$ be a constant vector, and let $f$ be the real-valued function defined as $\mathbf{c}^{\top}\mathbf{x}$, where $\mathbf{x}=\left[\begin{array}{c} x_{1}\\ \vdots\\ x_{n} \end{array}\right]$ is the variable with all the terms of $\mathbf{x}$ being nonnegative.

Let the polytope be defined as $A\mathbf{x}\leq\mathbf{b}$, where $A$ is a matrix with $k$ rows and $n$ columns and $\mathbf{b=\left[\begin{array}{c} b_{1}\\\vdots\\b_{k}\end{array}\right]}$, so the polytope has $k$ sides.

Let $\mathbb{S}$ be the subset of the polytope for which $f$ reaches a maximum. Let $\mathfrak{R}$ be the rotation group whose elements are matrices with $n$ rows and $n$ columns with determinant equal to $1$, and let the group action of the rotation group on the polytope be defined as

$\left(R,\mathbf{x}\right)\mapsto R\mathbf{x},$

where $R\mathbf{x}$ is a matrix multiplication. This map is a group action because matrix multiplication is associative. For each element $R$ of the rotation group, we write $R\mathbb{S}$ to mean all elements $R\mathbf{x}$ where $\mathbf{x}$ is an element of $\mathbb{S}$.

Is there an algorithm that can find non-identity elements $R$ of the rotation group that keep the solution subset $\mathbb{S}$ unchanged, meaning $R\mathbb{S}$ and $\mathbb{S}$ are the same? If such an algorithm exists, what are all of those non-identity elements? Do they, along with the identity element, form a subgroup of the rotation group? If these elements exist, write a Python program (or in any language of your choice) that computes these elements.

Topological spaces, compact spaces, and continuous functions

Before we tell you what the next challenge is, we need to tell you what topological spaces, compact spaces, and continuous functions are.

Topological spaces

Imagine you have the real line broken up into intervals that may be overlapping, and each number on the real line falls in one of those intervals. Given a number on the real line, you can have a notion of “closeness” to that number by finding an interval small enough that contains that number.

The collection of all those intervals is called a topology, and the real line with this topology is called a topological space.

real line - open intervals - standard topology

Formally, a topology on a set is a collection of subsets (or open subsets) of that set such that

  • the empty subset and the set is in that collection;
  • arbitrary unions of these subsets are in that collection;
  • finite intersections of these subsets are in that collection.

A topological space is a set along with a topology. A set can admit several topologies. For examples, the real line has the topology made of open intervals (the standard topology) and the topology made of singletons (the discrete topology).

Compact topological spaces

A topological space is compact if it cannot get too big, meaning the points in that space are “close” to each other, according to the topology. For example, any closed interval of real numbers is compact,

closed interval of real numbers - compact space

any circle in a $2$-dimensional real space is compact,

2d circle - compact space

any sphere in a $3$-dimensional real space is compact:

3d sphere - compact space

Formally, a compact topological space is a space with the property that any collection of open subsets that contains all of its points admits a finite subcollection of those open subsets that contains all of its points.

Continuous functions

Intuitively, continuous functions are functions such that images of close points (according to a topology) are close to each other.

More precisely, a function from a topological space to another topological space is continuous if, for each point $x$ in the domain of that function, all points close to $x$ that are also in the domain of that function admit images that are close to the image of $x$.

continuous functions - topological space

Challenge #2

We are now ready to tell you what the challenge is. The setting is the $n$-th dimensional real space with the Minkowski metric of order $p$ defined as for any two points $\mathbf{u}=\left(u_{1},\ldots,u_{n}\right)$ and $\mathbf{v}=\left(v_{1},\ldots,v_{n}\right)$

$\left(\left|u_{1}-v_{1}\right|^{p}+\cdots+\left|u_{n}-v_{n}\right|^{p}\right)^{\frac{1}{p}}.$

That metric yields a topology given by “open balls”. An open ball is a subset of the real space that has a fixed point and the distance between all other points in the ball and that fixed point is less than some positive real number, and the distance is given according to the Minkowski metric.

open ball - euclidean metric - topology

 

  1. Let $f$ be a continuous real-valued function on a compact subspace of the real space, where that subspace inherits the same topology as the real space. Is there an algorithm that can check if there are points in that compact subspace that maximizes $f$ and that can find all of those points? If such an algorithm exists, build a Python program (or in any language of your choice) that can find all of those points when $f$ is the function given by $\mathbf{c}^{\top}\mathbf{x}$ and the domain of $f$ is the unit sphere centered at the origin with radius equal to $1$.
  2. Let $\mathfrak{R}$ be the rotation group whose elements are matrices with $n$ rows and $n$ columns with determinant equal to $1$, and let the group action of the rotation group on the compact subspace be defined as $\left(R,\mathbf{x}\right)\mapsto R\mathbf{x},$ where $R\mathbf{x}$ is a matrix multiplication. This map is a group action because matrix multiplication is associative. Let $\mathbb{S}$ be the collection of all points in the compact subspace that maximize the continuous function. For each element $R$ of the rotation group, we write $R\mathbb{S}$ to mean all elements $R\mathbf{x}$ where $\mathbf{x}$ is an element of $\mathbb{S}$. Is there an algorithm that can find non-identity elements $R$ of the rotation group that keep the solution subset $\mathbb{S}$ unchanged, meaning $R\mathbb{S}$ and $\mathbb{S}$ are the same? If such an algorithm exists, what are all of those non-identity elements? Do they, along with the identity element, form a subgroup of the rotation group? If these elements exist, write a Python program (or in any language of your choice) that computes these elements.

References

[1] dantzigs-simplex-algorithm [Python package]: https://dantzigs-simplex-algorithm.readthedocs.io/en/latest/quickstart.html

[2] scipy.optimize.linprog. SciPy: https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.linprog.html

[3] Define and solve a problem by using Solver. Microsoft: https://support.microsoft.com/en-us/office/define-and-solve-a-problem-by-using-solver-5d1a388f-079d-43ac-a7eb-f63e45925040

[4] Excel Linear Programming (Using Solver and Graphical Methods). ExcelDemy: https://www.exceldemy.com/learn-excel/solver/linear-programming/

[5] Simplex. GitHub: https://github.com/chashikajw/simplex-algorithm/blob/master/Simplex.cpp

[6] linprog. MathWorks: https://www.mathworks.com/help/optim/ug/linprog.html

[7] A Mathematical Approach to Optimize Crop Allocation – A Linear Programming Model. International Information and Engineering Technology Association: https://iieta.org/journals/ijdne/paper/10.18280/ijdne.150215

[8] Optimal Airport Surface Traffic Planning Using Mixed-Integer Linear Programming. International Journal of Aerospace Engineering: https://onlinelibrary.wiley.com/doi/10.1155/2008/732828

[9] Supply Chain Process Optimization Using Linear Programming. Towards Data Science (Medium): https://towardsdatascience.com/supply-chain-process-optimization-using-linear-programming-b1511800630f

[10] Mixed-integer linear programming based optimization strategies for renewable energy communities. Energy (Elsevier): https://www.sciencedirect.com/science/article/pii/S0360544221018077

[11] Determinant. Wikipedia: https://en.wikipedia.org/wiki/Determinant