Theory behind the solver. While its knowledge is not necessary to use the solver, it can be useful to understand how it works, especially when debugging or extending its features.

Legend

To ensure an easier reading, we will use the following notation:

• bold capital letters $\boldsymbol{A}, \boldsymbol{B}, \boldsymbol{C}$ for matrices
• subscripted bold capital letters $\boldsymbol{A}_i, \boldsymbol{B}_i, \boldsymbol{C}_i$ for column vectors
• subscripted bold capital letters $\boldsymbol{A}_{i \dots j}, \boldsymbol{B}_{i \dots j}, \boldsymbol{C}_{i \dots j}$ for submatrices going from column $i$ to column $j$
• bold lowercase letters $\boldsymbol{c}, \boldsymbol{x}, \boldsymbol{y}, \boldsymbol{z}$ for vectors
• subscripted bold lowercase letters $\boldsymbol{c}_i, \boldsymbol{x}_i, \boldsymbol{y}_i, \boldsymbol{z}_i$ for components of vectors
• non-bold lowercase letters $a, b, c, x, y, z$ for scalars

Furthermore, a subscripted letter "$^T$" will indicate the transpose of the corresponding matrix or vector and a subscripted "$^{-1}$" will indicate the inverse of the corresponding matrix.

Floating-point arithmetic model

When it comes to dealing with real numbers, the limited available memory and the desire for a fast computation forces computers to approximate their representation. Almost all modern devices operate the IEEE 754 standard implementing floating-point arithmetic. It is designed in such a way that the evaluation of an expression, denoted as $\bold{fl}(\cdot)$, where the operator $\circ = +, -, -, \star$, satisfies

$$ \bold{fl}_\epsilon(x \circ y) = (x \circ y)(1 + z) \quad |z| \le \epsilon $$

where $\epsilon$ the the unit roundoff (or machine epsilon), which represents the smallest possible value the system can distinguish from $0$ and it is inversely exponential in the number of bits of precision. It is usually in the order of $10^{-8}$ or $10^{-16}$ in single and double precision computer arithmetic. The same notation can be extended to include arbitrary linear algebra expressions within $\bold{fl}_\epsilon(\cdot)$ and it is to be understood as that each scalar operation underlying the expression is to be carried out in that same precision, in an arbitrary order.

Linear Programming

A LP problem is an optimisation problem where the objective function and the constraints are linear equalities or inequalities. The objective function is what we want to maximise or minimise, while the constraints are the conditions the solution must satisfy. The standard form of a LP problem is the following:

$$ \begin{equation*} \begin{aligned} & \max & \boldsymbol{c}^T \boldsymbol{x} \newline & \text{subject to} & \boldsymbol{A} \boldsymbol{x} \leq \boldsymbol{b} \newline & & \boldsymbol{x} \geq 0 \end{aligned} \end{equation*} $$

where $\boldsymbol{x} \in \mathbb{R}^d$ is the vector of variables to be determined, $\boldsymbol{c} \in \mathbb{R}^d$ and $\boldsymbol{b} \in \mathbb{R}^n$ are vectors of coefficients, and $\boldsymbol{A} \in \mathbb{R}^{n \times d}$ is a matrix of coefficients. It is always possible to rewrite a LP problem in standard form following these steps:

• If the problem is a minimization problem, it is sufficient to multiply the objective function by $-1$ to obtain the corresponding maximization problem.
• If some variables have no lower bound, they can be substituted with the difference of two variables, both with a lower bound of $0$ (i.e. $x = x_1 - x_2$ and $x_1, x_2 \geq 0$).
• If strict equalities exist, they can be substituted with two inequalities (i.e. $x_i = 0$ becomes $x_i \leq 0$ and $x_i \geq 0$).
• If there are some inequalities with different signs, one can be multiplied by $-1$ to ensure both have the same sign.

LP problems are usually solved via the simplex method, although some interior-point methods may be used as well. The simplex method is an iterative algorithm developed by George Dantzig in 1947. To apply the simplex method, the problem must be converted in slack form to use the simplex method. Starting with a LP problem in standard form, the slack form is obtained by introducing a slack variable for each constraint so that the inequality becomes an equality. The slack variables must also be non-negative.

$$ \begin{equation*} \begin{aligned} & \max & \boldsymbol{c}^T \boldsymbol{x} \newline & \text{subject to} & \boldsymbol{A} \boldsymbol{x} + \boldsymbol{s} = \boldsymbol{b} \newline & & \boldsymbol{x} \geq 0 \newline & & \boldsymbol{s} \ge 0 \end{aligned} \end{equation*} $$

It is possible to extract an invertible squared matrix $B$ from $A$ of dimension $n \times n$ called basis. The remaining columns of $A$ form the $N$ matrix.

$$ \begin{equation*} \boldsymbol{A} = \begin{bmatrix} \boldsymbol{B} & \boldsymbol{N} \end{bmatrix} \text{, where } \begin{cases} \boldsymbol{B} \in \mathbb{R}^{n \times n} \newline \boldsymbol{N} \in \mathbb{R}^{n \times (d - n)} \end{cases} \end{equation*} $$

The variables $x_i$ whose index $i$ corresponds to a column of $A$ also belonging to $B$ are called basic variables, while the ones that end up in $N$ non-basic variables. For commodity, we can also define $I_B, I_N$ to be the set of indexes of such basic and non-basic variables respectively. An example is shown below.

$$ \boldsymbol{A} = \begin{bmatrix} a_{1,1} & a_{1,2} & a_{1,3} & a_{1,4} \newline a_{2,1} & a_{2,2} & a_{2,3} & a_{2,4} \newline a_{3,1} & a_{3,2} & a_{3,3} & a_{3,4} \end{bmatrix} \newline \boldsymbol{B} = \begin{bmatrix} a_{1,1} & a_{1,2} & a_{1,4} \newline a_{2,1} & a_{2,2} & a_{2,4} \newline a_{3,1} & a_{3,2} & a_{3,4} \end{bmatrix} \quad \boldsymbol{N} = \begin{bmatrix} a_{1,3} \newline a_{2,3} \newline a_{3,3} \end{bmatrix} \newline I_B =1, 2, 4