Python implementation of the simplex algorithm for solving linear programs in tabular form with  , , and  constraints and  each variable x1, x2, ... 0. See https:gist.github.comimengusf9619a568f7da5bc74eaf20169a24d98 for how to convert linear programs to simplex tableaus, and the steps taken in the simplex algorithm. Resources: https:en.wikipedia.orgwikiSimplexalgorithm https:tinyurl.comsimplex4beginners Operate on simplex tableaus  Tableaunp.array1,1,0,0,1,1,3,1,0,4,3,1,0,1,4, 2, 2 Traceback most recent call last: ... TypeError: Tableau must have type float64  Tableaunp.array1,1,0,0,1,1,3,1,0,4,3,1,0,1,4., 2, 2 Traceback most recent call last: ... ValueError: RHS must be  0  Tableaunp.array1,1,0,0,1,1,3,1,0,4,3,1,0,1,4., 2, 2 Traceback most recent call last: ... ValueError: number of artificial variables must be a natural number Max iteration number to prevent cycling Check if RHS is negative Number of decision variables x1, x2, x3... 2 if there are  or  constraints nonstandard, 1 otherwise std Number of slack variables added to make inequalities into equalities Objectives for each stage In two stage simplex, first minimise then maximise Index of current pivot row and column Does objective row only contain nonnegative values? Generate column titles for tableau of specific dimensions  Tableaunp.array1,1,0,0,1,1,3,1,0,4,3,1,0,1,4., ... 2, 0.generatecoltitles 'x1', 'x2', 's1', 's2', 'RHS'  Tableaunp.array1,1,0,0,1,1,3,1,0,4,3,1,0,1,4., ... 2, 2.generatecoltitles 'x1', 'x2', 'RHS' decision  slack Finds the pivot row and column.  Tableaunp.array2,1,0,0,0, 3,1,1,0,6, 1,2,0,1,7., ... 2, 0.findpivot 1, 0 Find entries of highest magnitude in objective rows Choice is only valid if below 0 for maximise, and above for minimise Pivot row is chosen as having the lowest quotient when elements of the pivot column divide the righthand side Slice excluding the objective rows RHS Elements of pivot column within slice Array filled with nans If element in pivot column is greater than zero, return quotient or nan otherwise Arg of minimum quotient excluding the nan values. nstages is added to compensate for earlier exclusion of objective columns Pivots on value on the intersection of pivot row and column.  Tableaunp.array2,3,0,0,0,1,3,1,0,4,3,1,0,1,4., ... 2, 2.pivot1, 0.tolist ...  doctest: NORMALIZEWHITESPACE 0.0, 3.0, 2.0, 0.0, 8.0, 1.0, 3.0, 1.0, 0.0, 4.0, 0.0, 8.0, 3.0, 1.0, 8.0 Avoid changes to original tableau Entry becomes 1 Variable in pivot column becomes basic, ie the only nonzero entry Exits first phase of the twostage method by deleting artificial rows and columns, or completes the algorithm if exiting the standard case.  Tableaunp.array ... 3, 3, 1, 1, 0, 0, 4, ... 2, 1, 0, 0, 0, 0, 0., ... 1, 2, 1, 0, 1, 0, 2, ... 2, 1, 0, 1, 0, 1, 2 ... , 2, 2.changestage.tolist ...  doctest: NORMALIZEWHITESPACE 2.0, 1.0, 0.0, 0.0, 0.0, 1.0, 2.0, 1.0, 0.0, 2.0, 2.0, 1.0, 0.0, 1.0, 2.0 Objective of original objective row remains Slice containing ids for artificial columns Delete the artificial variable columns Delete the objective row of the first stage Operate on tableau until objective function cannot be improved further. Standard linear program: Max:  x1   x2 ST:   x1  3x2  4 3x1   x2  4  Tableaunp.array1,1,0,0,0,1,3,1,0,4,3,1,0,1,4., ... 2, 0.runsimplex 'P': 2.0, 'x1': 1.0, 'x2': 1.0 Standard linear program with 3 variables: Max: 3x1   x2  3x3 ST:  2x1   x2   x3  2 x1  2x2  3x3  5 2x1  2x2   x3  6  Tableaunp.array ... 3,1,3,0,0,0,0, ... 2,1,1,1,0,0,2, ... 1,2,3,0,1,0,5, ... 2,2,1,0,0,1,6. ... ,3,0.runsimplex  doctest: ELLIPSIS 'P': 5.4, 'x1': 0.199..., 'x3': 1.6 Optimal tableau input:  Tableaunp.array ... 0, 0, 0.25, 0.25, 2, ... 0, 1, 0.375, 0.125, 1, ... 1, 0, 0.125, 0.375, 1 ... , 2, 0.runsimplex 'P': 2.0, 'x1': 1.0, 'x2': 1.0 Nonstandard:  constraints Max: 2x1  3x2   x3 ST:   x1   x2   x3  40 2x1   x2   x3  10   x2   x3  10  Tableaunp.array ... 2, 0, 0, 0, 1, 1, 0, 0, 20, ... 2, 3, 1, 0, 0, 0, 0, 0, 0, ... 1, 1, 1, 1, 0, 0, 0, 0, 40, ... 2, 1, 1, 0, 1, 0, 1, 0, 10, ... 0, 1, 1, 0, 0, 1, 0, 1, 10. ... , 3, 2.runsimplex 'P': 70.0, 'x1': 10.0, 'x2': 10.0, 'x3': 20.0 Non standard: minimisation and equalities Min: x1   x2 ST: 2x1   x2  12 6x1  5x2  40  Tableaunp.array ... 8, 6, 0, 0, 52, ... 1, 1, 0, 0, 0, ... 2, 1, 1, 0, 12, ... 6, 5, 0, 1, 40., ... , 2, 2.runsimplex 'P': 7.0, 'x1': 5.0, 'x2': 2.0 Pivot on slack variables Max: 8x1  6x2 ST:   x1  3x2  33 4x1  2x2  48 2x1  4x2  48 x1   x2  10 x1         2  Tableaunp.array ... 2, 1, 0, 0, 0, 1, 1, 0, 0, 12.0, ... 8, 6, 0, 0, 0, 0, 0, 0, 0, 0.0, ... 1, 3, 1, 0, 0, 0, 0, 0, 0, 33.0, ... 4, 2, 0, 1, 0, 0, 0, 0, 0, 60.0, ... 2, 4, 0, 0, 1, 0, 0, 0, 0, 48.0, ... 1, 1, 0, 0, 0, 1, 0, 1, 0, 10.0, ... 1, 0, 0, 0, 0, 0, 1, 0, 1, 2.0 ... , 2, 2.runsimplex  doctest: ELLIPSIS 'P': 132.0, 'x1': 12.000... 'x2': 5.999... Stop simplex algorithm from cycling. Completion of each stage removes an objective. If both stages are complete, then no objectives are left Find the values of each variable at optimal solution If there are no more negative values in objective row Delete artificial variable columns and rows. Update attributes Given the final tableau, add the corresponding values of the basic decision variables to the outputdict  Tableaunp.array ... 0,0,0.875,0.375,5, ... 0,1,0.375,0.125,1, ... 1,0,0.125,0.375,1 ... ,2, 0.interprettableau 'P': 5.0, 'x1': 1.0, 'x2': 1.0 P  RHS of final tableau Gives indices of nonzero entries in the ith column First entry in the nonzero indices If there is only one nonzero value in column, which is one
from typing import Any

import numpy as np


class Tableau:
    """Operate on simplex tableaus

    >>> Tableau(np.array([[-1,-1,0,0,1],[1,3,1,0,4],[3,1,0,1,4]]), 2, 2)
    Traceback (most recent call last):
    ...
    TypeError: Tableau must have type float64

    >>> Tableau(np.array([[-1,-1,0,0,-1],[1,3,1,0,4],[3,1,0,1,4.]]), 2, 2)
    Traceback (most recent call last):
    ...
    ValueError: RHS must be > 0

    >>> Tableau(np.array([[-1,-1,0,0,1],[1,3,1,0,4],[3,1,0,1,4.]]), -2, 2)
    Traceback (most recent call last):
    ...
    ValueError: number of (artificial) variables must be a natural number
    """

    # Max iteration number to prevent cycling
    maxiter = 100

    def __init__(
        self, tableau: np.ndarray, n_vars: int, n_artificial_vars: int
    ) -> None:
        if tableau.dtype != "float64":
            raise TypeError("Tableau must have type float64")

        # Check if RHS is negative
        if not (tableau[:, -1] >= 0).all():
            raise ValueError("RHS must be > 0")

        if n_vars < 2 or n_artificial_vars < 0:
            raise ValueError(
                "number of (artificial) variables must be a natural number"
            )

        self.tableau = tableau
        self.n_rows, n_cols = tableau.shape

        # Number of decision variables x1, x2, x3...
        self.n_vars, self.n_artificial_vars = n_vars, n_artificial_vars

        # 2 if there are >= or == constraints (nonstandard), 1 otherwise (std)
        self.n_stages = (self.n_artificial_vars > 0) + 1

        # Number of slack variables added to make inequalities into equalities
        self.n_slack = n_cols - self.n_vars - self.n_artificial_vars - 1

        # Objectives for each stage
        self.objectives = ["max"]

        # In two stage simplex, first minimise then maximise
        if self.n_artificial_vars:
            self.objectives.append("min")

        self.col_titles = self.generate_col_titles()

        # Index of current pivot row and column
        self.row_idx = None
        self.col_idx = None

        # Does objective row only contain (non)-negative values?
        self.stop_iter = False

    def generate_col_titles(self) -> list[str]:
        """Generate column titles for tableau of specific dimensions

        >>> Tableau(np.array([[-1,-1,0,0,1],[1,3,1,0,4],[3,1,0,1,4.]]),
        ... 2, 0).generate_col_titles()
        ['x1', 'x2', 's1', 's2', 'RHS']

        >>> Tableau(np.array([[-1,-1,0,0,1],[1,3,1,0,4],[3,1,0,1,4.]]),
        ... 2, 2).generate_col_titles()
        ['x1', 'x2', 'RHS']
        """
        args = (self.n_vars, self.n_slack)

        # decision | slack
        string_starts = ["x", "s"]
        titles = []
        for i in range(2):
            for j in range(args[i]):
                titles.append(string_starts[i] + str(j + 1))
        titles.append("RHS")
        return titles

    def find_pivot(self) -> tuple[Any, Any]:
        """Finds the pivot row and column.
        >>> Tableau(np.array([[-2,1,0,0,0], [3,1,1,0,6], [1,2,0,1,7.]]),
        ... 2, 0).find_pivot()
        (1, 0)
        """
        objective = self.objectives[-1]

        # Find entries of highest magnitude in objective rows
        sign = (objective == "min") - (objective == "max")
        col_idx = np.argmax(sign * self.tableau[0, :-1])

        # Choice is only valid if below 0 for maximise, and above for minimise
        if sign * self.tableau[0, col_idx] <= 0:
            self.stop_iter = True
            return 0, 0

        # Pivot row is chosen as having the lowest quotient when elements of
        # the pivot column divide the right-hand side

        # Slice excluding the objective rows
        s = slice(self.n_stages, self.n_rows)

        # RHS
        dividend = self.tableau[s, -1]

        # Elements of pivot column within slice
        divisor = self.tableau[s, col_idx]

        # Array filled with nans
        nans = np.full(self.n_rows - self.n_stages, np.nan)

        # If element in pivot column is greater than zero, return
        # quotient or nan otherwise
        quotients = np.divide(dividend, divisor, out=nans, where=divisor > 0)

        # Arg of minimum quotient excluding the nan values. n_stages is added
        # to compensate for earlier exclusion of objective columns
        row_idx = np.nanargmin(quotients) + self.n_stages
        return row_idx, col_idx

    def pivot(self, row_idx: int, col_idx: int) -> np.ndarray:
        """Pivots on value on the intersection of pivot row and column.

        >>> Tableau(np.array([[-2,-3,0,0,0],[1,3,1,0,4],[3,1,0,1,4.]]),
        ... 2, 2).pivot(1, 0).tolist()
        ... # doctest: +NORMALIZE_WHITESPACE
        [[0.0, 3.0, 2.0, 0.0, 8.0],
        [1.0, 3.0, 1.0, 0.0, 4.0],
        [0.0, -8.0, -3.0, 1.0, -8.0]]
        """
        # Avoid changes to original tableau
        piv_row = self.tableau[row_idx].copy()

        piv_val = piv_row[col_idx]

        # Entry becomes 1
        piv_row *= 1 / piv_val

        # Variable in pivot column becomes basic, ie the only non-zero entry
        for idx, coeff in enumerate(self.tableau[:, col_idx]):
            self.tableau[idx] += -coeff * piv_row
        self.tableau[row_idx] = piv_row
        return self.tableau

    def change_stage(self) -> np.ndarray:
        """Exits first phase of the two-stage method by deleting artificial
        rows and columns, or completes the algorithm if exiting the standard
        case.

        >>> Tableau(np.array([
        ... [3, 3, -1, -1, 0, 0, 4],
        ... [2, 1, 0, 0, 0, 0, 0.],
        ... [1, 2, -1, 0, 1, 0, 2],
        ... [2, 1, 0, -1, 0, 1, 2]
        ... ]), 2, 2).change_stage().tolist()
        ... # doctest: +NORMALIZE_WHITESPACE
        [[2.0, 1.0, 0.0, 0.0, 0.0],
        [1.0, 2.0, -1.0, 0.0, 2.0],
        [2.0, 1.0, 0.0, -1.0, 2.0]]
        """
        # Objective of original objective row remains
        self.objectives.pop()

        if not self.objectives:
            return self.tableau

        # Slice containing ids for artificial columns
        s = slice(-self.n_artificial_vars - 1, -1)

        # Delete the artificial variable columns
        self.tableau = np.delete(self.tableau, s, axis=1)

        # Delete the objective row of the first stage
        self.tableau = np.delete(self.tableau, 0, axis=0)

        self.n_stages = 1
        self.n_rows -= 1
        self.n_artificial_vars = 0
        self.stop_iter = False
        return self.tableau

    def run_simplex(self) -> dict[Any, Any]:
        """Operate on tableau until objective function cannot be
        improved further.

        # Standard linear program:
        Max:  x1 +  x2
        ST:   x1 + 3x2 <= 4
             3x1 +  x2 <= 4
        >>> Tableau(np.array([[-1,-1,0,0,0],[1,3,1,0,4],[3,1,0,1,4.]]),
        ... 2, 0).run_simplex()
        {'P': 2.0, 'x1': 1.0, 'x2': 1.0}

        # Standard linear program with 3 variables:
        Max: 3x1 +  x2 + 3x3
        ST:  2x1 +  x2 +  x3 ≤ 2
              x1 + 2x2 + 3x3 ≤ 5
             2x1 + 2x2 +  x3 ≤ 6
        >>> Tableau(np.array([
        ... [-3,-1,-3,0,0,0,0],
        ... [2,1,1,1,0,0,2],
        ... [1,2,3,0,1,0,5],
        ... [2,2,1,0,0,1,6.]
        ... ]),3,0).run_simplex() # doctest: +ELLIPSIS
        {'P': 5.4, 'x1': 0.199..., 'x3': 1.6}


        # Optimal tableau input:
        >>> Tableau(np.array([
        ... [0, 0, 0.25, 0.25, 2],
        ... [0, 1, 0.375, -0.125, 1],
        ... [1, 0, -0.125, 0.375, 1]
        ... ]), 2, 0).run_simplex()
        {'P': 2.0, 'x1': 1.0, 'x2': 1.0}

        # Non-standard: >= constraints
        Max: 2x1 + 3x2 +  x3
        ST:   x1 +  x2 +  x3 <= 40
             2x1 +  x2 -  x3 >= 10
                 -  x2 +  x3 >= 10
        >>> Tableau(np.array([
        ... [2, 0, 0, 0, -1, -1, 0, 0, 20],
        ... [-2, -3, -1, 0, 0, 0, 0, 0, 0],
        ... [1, 1, 1, 1, 0, 0, 0, 0, 40],
        ... [2, 1, -1, 0, -1, 0, 1, 0, 10],
        ... [0, -1, 1, 0, 0, -1, 0, 1, 10.]
        ... ]), 3, 2).run_simplex()
        {'P': 70.0, 'x1': 10.0, 'x2': 10.0, 'x3': 20.0}

        # Non standard: minimisation and equalities
        Min: x1 +  x2
        ST: 2x1 +  x2 = 12
            6x1 + 5x2 = 40
        >>> Tableau(np.array([
        ... [8, 6, 0, 0, 52],
        ... [1, 1, 0, 0, 0],
        ... [2, 1, 1, 0, 12],
        ... [6, 5, 0, 1, 40.],
        ... ]), 2, 2).run_simplex()
        {'P': 7.0, 'x1': 5.0, 'x2': 2.0}


        # Pivot on slack variables
        Max: 8x1 + 6x2
        ST:   x1 + 3x2 <= 33
             4x1 + 2x2 <= 48
             2x1 + 4x2 <= 48
              x1 +  x2 >= 10
             x1        >= 2
        >>> Tableau(np.array([
        ... [2, 1, 0, 0, 0, -1, -1, 0, 0, 12.0],
        ... [-8, -6, 0, 0, 0, 0, 0, 0, 0, 0.0],
        ... [1, 3, 1, 0, 0, 0, 0, 0, 0, 33.0],
        ... [4, 2, 0, 1, 0, 0, 0, 0, 0, 60.0],
        ... [2, 4, 0, 0, 1, 0, 0, 0, 0, 48.0],
        ... [1, 1, 0, 0, 0, -1, 0, 1, 0, 10.0],
        ... [1, 0, 0, 0, 0, 0, -1, 0, 1, 2.0]
        ... ]), 2, 2).run_simplex() # doctest: +ELLIPSIS
        {'P': 132.0, 'x1': 12.000... 'x2': 5.999...}
        """
        # Stop simplex algorithm from cycling.
        for _ in range(Tableau.maxiter):
            # Completion of each stage removes an objective. If both stages
            # are complete, then no objectives are left
            if not self.objectives:
                # Find the values of each variable at optimal solution
                return self.interpret_tableau()

            row_idx, col_idx = self.find_pivot()

            # If there are no more negative values in objective row
            if self.stop_iter:
                # Delete artificial variable columns and rows. Update attributes
                self.tableau = self.change_stage()
            else:
                self.tableau = self.pivot(row_idx, col_idx)
        return {}

    def interpret_tableau(self) -> dict[str, float]:
        """Given the final tableau, add the corresponding values of the basic
        decision variables to the `output_dict`
        >>> Tableau(np.array([
        ... [0,0,0.875,0.375,5],
        ... [0,1,0.375,-0.125,1],
        ... [1,0,-0.125,0.375,1]
        ... ]),2, 0).interpret_tableau()
        {'P': 5.0, 'x1': 1.0, 'x2': 1.0}
        """
        # P = RHS of final tableau
        output_dict = {"P": abs(self.tableau[0, -1])}

        for i in range(self.n_vars):
            # Gives indices of nonzero entries in the ith column
            nonzero = np.nonzero(self.tableau[:, i])
            n_nonzero = len(nonzero[0])

            # First entry in the nonzero indices
            nonzero_rowidx = nonzero[0][0]
            nonzero_val = self.tableau[nonzero_rowidx, i]

            # If there is only one nonzero value in column, which is one
            if n_nonzero == 1 and nonzero_val == 1:
                rhs_val = self.tableau[nonzero_rowidx, -1]
                output_dict[self.col_titles[i]] = rhs_val
        return output_dict


if __name__ == "__main__":
    import doctest

    doctest.testmod()
