The Simplex Method

The Simplex Method is commonly known to mathematicians as an efficient way to solve a linear program that includes an optimatization function and several constraints. According to Stephen J. Wright in his article titled, "simplex method", written in 2021, "In practice, problems often involve hundreds of equations with thousands of variables, which can result in an astronomical number of extreme points". However, by utilizing the simplex method, the extreme points can be limited as to provide a way to easily examine and solve the constraints. Wright refers to the simplex method as, "one of the most useful and efficient algorithms ever invented, and it is still the standard method employed on computers to solve optimization problems".

About the Method

During World War II, an American scientist known as George Dantzig created the method when trying to assist the US Army Air Force with mechanizing a planning process. Since then, the method has been studied and used consistently to solve linear programs. The first step in performing the simplex method is ensuring that the program is in canonical form, which is shown below:

\[\begin{align*} max& && c^Tx\\ s.t.& &&Ax \le b\\ &&&x \ge 0, \end{align*}\]

with \(c\) = (\(c_{1}\), ... , \(c_{n}\)) representing the coefficients of the above objective function and \(x\) = (\(x_{1}\), ... , \(x_{n}\)) representing the variables of the linear program. The \(A\) represents an \(m \times n\) matrix with the constraints of \(b\) = (\(b_{1}\), ... , \(b_{m}\)).

Extreme points are defined as the points that satisfy maximizing the objective function. In a linear program there are a finite number of extreme points, but in the case of an extensive linear program, the amount of extreme points can seem almost infinitely large. When solving the program, the extreme point can either be a point on the graph (maximizes the objective function) or it can represent an edge of the graph in which the graph keeps moving (does not maximize the objective function). In the case where it represents an edge of the graph, the simplex method allows the values to move along the edge until a maximum value is found, or an unbounded case arises. Below is a couple of visualizations of such cases:

Overview of the Project

For our Linear Programming Final Project, we selected Problem #4 - The Simplex Method. We used R and RStudio coding language accompanied wtih the software, RPubs, that allowed us to format our report. In this project, we implemented our own simplex method. Our algorithm will be able to solve a standard form linear program with full-rank assumption. Below is an example of the standard form of a linear program:

\[\begin{align*} max& && c^Tx\\ s.t.& &&Ax=b\\ &&&x \ge 0, \end{align*}\]

where \(x\), \(c\) \(\in\) \(\mathbb{R}^n\), \(b\) \(\in\) \(\mathbb{R^m}\), \(A\) is an \(m \times n\) matrix and rank ( \(A\)) = \(m < n\).

Question #1

Implement a preliminary simplex method for LP in canonical form

Objectives:

  • Inputs: \(A\), \(c\), \(b\) such that the LP is in canonical form.
  • Outputs:
    • flag: \(1\), if the problem has an optimal solution; \(-3\), if the problem is unbounded.
    • \(x\)*: an optimal solution (if flag = 1).
    • \(z\)*: the optimal value (if flag = 1).
  • When choosing variables, use Dantzig's rule: choose a variable with the most negative reduced cost.
  • To break a tie in the minimum ratio test, choose the one with the smallest subscript among all eligible variables.
  • Do not use any existing linear program solvers/functions in your code.

Code Solution

simplex <- function(A, b, c) {
  m = dim(A)[1]
  n = dim(A)[2] - m
  indices = 1:(m+n)
  
  # step 1: pick a feasible corner
  combs = t(combn(indices, n))
  feasible = FALSE
  for (i in 1:choose(m+n, n)) {
    free = combs[i,]     # the indices of the free variables (those set to 0)
    pivots = setdiff(indices, free)    # indices of the pivots (not free)
    
    corner = rep(0,m+n)
    corner[pivots] = solve(A[,pivots], b)
    if (all(corner >= 0)) {
      feasible = TRUE
      break
    }
  }
  if (!feasible) {
    return('flag: -3')
  }
  
  
  # step 2: slide around the edges until the optimal corner is found
  N <- A[,free]
  B <- A[,pivots]
  c_n <- c[free]
  c_b <- c[pivots]
  
  r <- c_n - (c_b %*% solve(B) %*% N)
  while (any(r < 0)) {
    # the total cost at this corner is greater than the cost of including one of
    # the free variables instead of one of the pivots. the best two to swap are:
    entering <- which.min(r)
    
    B_inv <- solve(B)
    u <- N[,entering]
    z <- (B_inv %*% b)/(B_inv %*% u)
    leaving <- which.min(z)
    
    # swap entering and leaving
    temp <- pivots[leaving]
    pivots[leaving] <- free[entering]
    free[entering] <- temp
    
    N <- A[,free]
    B <- A[,pivots]
    c_n <- c[free]
    c_b <- c[pivots]
    
    r <- c_n - (c_b %*% solve(B) %*% N)
  }
  
  corner = rep(0,m+n)
  corner[pivots] = solve(A[,pivots], b)
  
  print('flag: 1')
  print('X*:')
  print(corner)
  z.value <- c %*% corner
  z.value <- as.integer(z.value)
  print('z*:')
  print(z.value)
}

Example of Optimal Solution

A = rbind(
  c(1,2,-1,0),
  c(2,1,0,-1))
b = c(6,6)
cost = c(1,1,0,0)
simplex(A, b, cost)
## [1] "flag: 1"
## [1] "X*:"
## [1] 2 2 0 0
## [1] "z*:"
## [1] 4

Example of Unbounded Case

A = rbind(
  c(1,1),
  c(-1,1))
b = c(1,2)
cost = c(1,1)
simplex(A, b, cost)
## [1] "flag: -3"

Question #2

Bland's Rule

Originally, when solving a linear model, you start with a basic feasible solution and then change the basis in order to increase the optimization value until it reaches its maximal value. However, it is possible to attain a basic feasible solution that results in changing the basis but never changing the optimization value. This is called a degenerate linear program which will result in cycling. In otherwords, the problem seems to get "stuck" on a basic feasible solution no matter the basis and ultimately, never increases the optimal value. Bland's rule was founded as a solution to this cycling problem when optimizing a linear program.

Bland's rule solves cycling by creating a special "rule" for choosing which column will enter the basis. This "rule" is that when the user is choosing the entering variable \(x_{j}\), they must choose the smallest \(j\) with respect to the reduced costs \(\bar{c}_{j} < 0\). This rule also relates to the departing variable as well.

Objectives

Revise your code so that users can choose whether to apply Bland's rule for anticycling. (An extra optional input argument should be added.)

  • Inputs:
    • required: \(A\), \(c\), \(b\) such that the LP is in canonical form
    • optional: binary variable bland. (Bland's rule will be applied if bland = \(1\).)

Code Solution

simplex <- function(A, b, c, bland) {
  if(bland==1){
    m = dim(A)[1]
    n = dim(A)[2] - m
    indices = 1:(m+n)
    
    # step 1: pick a feasible corner
    combs = t(combn(indices, n))
    feasible = FALSE
    for (i in 1:choose(m+n, n)) {
      free = combs[i,]     # the indices of the free variables (those set to 0)
      pivots = setdiff(indices, free)    # indices of the pivots (not free)
      
      corner = rep(0,m+n)
      corner[pivots] = solve(A[,pivots], b)
      if (all(corner >= 0)) {
        feasible = TRUE
        break
      }
    }
    if (!feasible) {
      return('flag: -3')
    }
    
    
    # step 2: slide around the edges until the optimal corner is found
    N <- A[,free]
    B <- A[,pivots]
    c_n <- c[free]
    c_b <- c[pivots]
    
    r <- c_n - (c_b %*% solve(B) %*% N)
    while (any(r < 0)) {
      # the total cost at this corner is greater than the cost of including one of
      # the free variables instead of one of the pivots. the best two to swap are:
      entering <- which.min(r)
      
      B_inv <- solve(B)
      u <- N[,entering]
      z <- (B_inv %*% b)/(B_inv %*% u)
      leaving <- which.min(z)
      
      # swap entering and leaving
      temp <- pivots[leaving]
      pivots[leaving] <- free[entering]
      free[entering] <- temp
      
      N <- A[,free]
      B <- A[,pivots]
      c_n <- c[free]
      c_b <- c[pivots]
      
      r <- c_n - (c_b %*% solve(B) %*% N)
    }
    
    corner = rep(0,m+n)
    corner[pivots] = solve(A[,pivots], b)
    
    print('flag: 1')
    print('X*:')
    print(corner)
    z.value <- c %*% corner
    z.value <- as.integer(z.value)
    print('z*:')
    print(z.value)
    
  }else{
    m = dim(A)[1]
    n = dim(A)[2] - m
    indices = 1:(m+n)
    
    # step 1: pick a feasible corner
    combs = t(combn(indices, n))
    feasible = FALSE
    for (i in 1:choose(m+n, n)) {
      free = combs[i,]     # the indices of the free variables (those set to 0)
      pivots = setdiff(indices, free)    # indices of the pivots (not free)
      
      corner = rep(0,m+n)
      corner[pivots] = solve(A[,pivots], b)
      if (all(corner >= 0)) {
        feasible = TRUE
        break
      }
    }
    if (!feasible) {
      return('flag: -3')
    }
    
    # step 2: slide around the edges until the optimal corner is found
    N <- A[,free]
    B <- A[,pivots]
    c_n <- c[free]
    c_b <- c[pivots]
    
    r <- c_n - (c_b %*% solve(B) %*% N)
    while (any(r < 0)) {
      # the total cost at this corner is greater than the cost of including one of
      # the free variables instead of one of the pivots. the best two to swap are:
      entering <- which.min(r)
      
      B_inv <- solve(B)
      u <- N[,entering]
      z <- (B_inv %*% b)/(B_inv %*% u)
      leaving <- which.min(z)
      
      # swap entering and leaving
      temp <- pivots[leaving]
      pivots[leaving] <- free[entering]
      free[entering] <- temp
      
      N <- A[,free]
      B <- A[,pivots]
      c_n <- c[free]
      c_b <- c[pivots]
      
      r <- c_n - (c_b %*% solve(B) %*% N)
    }
    
    corner = rep(0,m+n)
    corner[pivots] = solve(A[,pivots], b)
    
    print('flag: 1')
    print('X*:')
    print(corner)
    z.value <- c %*% corner
    z.value <- as.integer(z.value)
    print('z*:')
    print(z.value)
  }
}

Example

A = rbind(
  c(1,2,-1,0),
  c(2,1,0,-1))
b = c(6,6)
cost = c(1,1,0,0)
simplex(A, b, cost, 1)
## [1] "flag: 1"
## [1] "X*:"
## [1] 2 2 0 0
## [1] "z*:"
## [1] 4

Conclusion

In our project we have discussed the history, development, and the purpose of the simplex method. We have also shown the various ways the method can be utilized in solving a linear program and the different cases that appear. With our code we were able to take the following inputs: a matrix \(A\), the coefficients \(c_{j}\), and the constraints \(b_{j}\), solve the linear program, and output a value, or set of values, based on the optimization results. In addition, we were able to add an additional input to demonstrate Bland's Rule that allowed us to find the optimization results by avoiding cycling.

Our code was based on the following source: https://gist.github.com/masonicboom/406609 . This site was able to assist us in a foundation for our code and help us brainstorm along the way.

According to George Dantzig, "True optimization is the revolutionary contribution of modern research to decision processes". From our research on the simplex method we were able to discover how the method has transformed over time and has developed into what it is today. Without revolutionary thoughts and collobration, the simplex method would not have been discovered. Due to new technology and modern day developed mathematicians, the simplex method will continue to develop and mature far into the future.