# LPsimplex.wxmx

load(draw)$  set_draw_defaults(line_width=2, grid = [2,2], point_type = filled_circle,
          head_type = 'nofilled, head_angle = 20, head_length = 0.5,
             background_color = light_gray,  draw_realpart=false)$     
 fpprintprec : 5$  ratprint : false$

load (simplex);

load ("Econ1.mac");

## Preface

LPsimplex.wxmx uses Maxima to introduce the simplex method of solving "Linear
Programming" problems - optimization problems in which both the objective function
and the constraint conditions are linear functions of the input variables (decision
variables).

Our tableaux follow (most closely) the tableaux conventions used by Sergiy Butenko
of Texas A&M University in his YouTube videos (see references).

1. We begin with graphical solutions to two variable LP's.

2. We next demonstrate the powerful maxima functions:
   maximize_lp (obj, condL, posL) and minimize (obj, condL, posL)
from Maxima's simplex package

3. We next demonstrate Maxima functions depending on row operations to generate the 
next optimization tableau, functions tableau(RowL) and pivot1(RowL, [row, col]). The
latter function executes one simplex step based on selection of the pivot row and pivot
column. Both maximization and minimization simplex optimization examples are worked
through in detail.

4. Examples of degeneracy and cycling are then discussed.

5. The two-phase simplex method of handling handling optimization problems in which
artificial variables are needed is discussed next. This includes an unbounded example.

6. the "Big M" method of handling optimization problems which need artificial variables
is next discussed. This includes an unbounded example.

7. The matrix simplex method is demonstrated with the Maxima functions Mdefine(),
Mdisplay(), Mtableau(), Mratio(col), and Mpivot(Enter, Leave). Examples of using
the matrix method to arrive at the optimized solution are worked out in detail. The 
matrix method is simple to use. We include an example of the use of toMTab(rows),
which converts the display produced by tableau(rows) into a matrix tableau which is
used in the matrix method.

A code file Econ1.mac is available in the Economic Analysis with Maxima section of
my CSULB webpage. This code file defines the Maxima functions used in this worksheet.

Inside your worksheet, use load ("Econ1.mac");

This worksheet is one of a number of wxMaxima files available in the section
    Economic Analysis with Maxima
on my CSULB webpage.

We have inserted this worksheet in the middle of Dowling's Ch. 13 worksheets,
since it seems appropriate to discuss the solution of linear optimization (LP) problems
before the solution of nonlinear optimization problems.

Edwin L. (Ted) Woollett
http://home.csulb.edu/~woollett/
Nov. 14, 2021

## References

Operations Research, 2nd ed., Richard Bronson and Govindasami Naadimuthu,
 Schaum's Outline Series, McGraw-Hill, 1982, Ch. 1, 2, 3.

Sergiy Butenko YouTube videos, 
ISEN 320, Operations Research I,
https://www.youtube.com/
               watch?v=J6ocZbsqGqE&list=PLY9yf2-4yyeT0Jh80w5yuevm0o9lAAUJg


A Simple LP Model: The Heavenly Pouch Example  =   Video 5
   Problem description
   Summary of the problem data
   Introducing the decision variables
   Formulating the objective function
   Stating the constraints
   Some terminology
   Summary of the model
   Definition of a linear function

Graphical Solution for Two Variable LP's = Video 13
   Plotting the constraints
   Plotting the objective (iso-profit line) and finding the optimal solution
   LP with alternative optima 
   A minimization LP with unbounded feasible region
   An unbounded LP
   An infeasible LP
   Classification of LPs

The Simplex Method Part 1  =   Video 15, 
   Converting an LP into the standard form
   Dictionary representation of an LP in the standard form
   The initial basic feasible solution 
   Improving the solution: Selecting an entering variable
   Ensuring the feasibility of the new solution: The ratio test to determine the leaving 
             variable
   Performing the first pivot
   Performing the second step of the simplex method
   Summary of the steps for the considered example
   Detecting the optimality and reporting the solution

The Simplex Method Part 2  =   Video 16, 
   Solving the Heavenly Pouch example using the simplex method
  Tableau format: Setting up the initial tableau
  Comparing the dictionary and tableau formats
  Performing a simplex pivot in the tableau format
  Summary of the simplex step in the tableau format
  Proceeding to step 2 of the simplex; optimal solution
  Termination: detecting optimality
  Termination: detecting unboundedness

The Simplex Method, Part II (continued)  =  Video 17, 
    Intro
    Termination: Degeneracy and cycling
    Avoiding cycling: Bland's rule
    Potential pitfalls of the simplex method
    Issues with initialization of the simplex method

The Simplex Method, Part III, Video 1 =  Video 19, 
   Intro
   Introducing an artificial variable
   Auxiliary and Big-M problems
   The idea of the two-phase simplex method
   The idea of the Big-M method

The Simplex Method, Part III, Video 2, 
                      The Two-Phase Simplex and Big-M Method = Video 20

The Two-phase Simplex Method: An Example = Video 21

An Example of the Big-M Method = Video 22
 
Alternative Optima and the Simplex Method = Video 24

The Simplex Method in the Matrix Form  =  Video 39
=================================
Sergiy Butenko, ISEN 620, A Survey of Optimization, 82 videos
https://www.youtube.com/playlist?list=PLY9yf2-4yyeQTLkCVFnGedyERuCjKW7kI

## Graphical Solution for Two Variable LP's

We follow Butenko's Video 13.
In Video 5 Butenko's motivates the "Heavenly Pouch" example, in which two types of
baby slings (carriers)  are considered for manufacture and sale. 

### Graphical Example 1, Heavenly Pouch LP,
       Maximize Profit, Scenario 1 

Type 1 carrier is a cheaper solid color (inside and out), let x1 equal the number of units 
to manufacture. This is called a non-reversible baby carrier (NR) which is to sell for $23 
per carrier while costing $8 each to manufacture using 2 yards of solid color fabric (and 
no printed color fabric). The total amount of solid color fabric on hand (for both types of 
carriers) is 900 yd. Marketing projections are that at most 350 NR carriers can be sold.

Type 2 carrier is the more expensive model, let x2 be the number of units to manufacture.
This is called a reversible baby carrier (R) with 2 yd solid color fabric needed for one side
and 2 yd printed color fabric for the other side, designed to sell for $35 per carrier, 
costing $10 per carrier to manufacture, and subject to the fabric available restriction that
only 600 yd of printed color fabric are on hand.

A general carrier project budget restriction is that at most $4000 can be spent on 
manufacturing costs including both types.

The economic goal is to maximize the manufacturer's profit = revenue - costs
 = (23*x1 + 35*x2) - (8*x1 + 10*x2) = 15*x1 + 25*x2, subject to four constraints.

1. The solid color fabric availability constraint implies
   2*x1 + 2*x2 <= 900,  or x1 + x2 <= 450.

2. The printed color fabric availability constraint implies
   2*x2 <= 600, or x2 <= 300.

3. The project budget constraint implies
   8*x1 + 10*x2 <= 4000,  or  4*x1 + 5*x2 <= 2000.

4. The marketing projections for NR carriers implies
   x1 <= 350.

In summary, our LP problem (Linear Programming problem) is:

  Maximize  z  = 15*x1 + 25*x2,
   subject to:
    x1 + x2 <= 450,
    x2 <= 300,
    4*x1 + 5*x2 <= 2000,
       x1 <= 350,
   and the decision variables x1 and x2, being numbers of units to manufacture of each
 type, are inherently non-negative numbers, so we need to add:
    x1, x2 >= 0.

Graphical Solution
1. Plot the constraints.

Let x1 >= 0 label the horizontal "x-axis" and x2 >= 0 label the vertical "y-axis". The 
optimum value of z = 15*x1 + 25*x2 must then lie in the first quadrant of the (x1, x2) 
plane.

xmax : 500$
ymax : 500$
wxdraw2d (xlabel = "x1", ylabel = "x2", yrange = [0, ymax], key_pos = bottom_left,
    title = "Constraints",
    key = "x1 + x2 <= 450", explicit (450 - x1, x1, 0, xmax),
    color = red, key = "x2 <= 300", explicit (300, x1, 0, xmax),
    color = brown, key = "4*x1 + 5*x2 <= 2000", explicit (400 - 4*x1/5, x1, 0, xmax),
    color = magenta, key = "x1 <= 350", parametric (350, yy,  yy, 0, ymax))$

2. Add "iso-profit lines" z = constant.

For example, add a line corresponding to z  = 15*x1 + 25*x2 = 3000.
The symbols xmax and ymax remain bound to 500 respectively.

[xmax, ymax];
wxdraw2d (xlabel = "x1", ylabel = "x2", yrange = [0, ymax], 
    title = "z = profit = 15*x1 + 25*x2",
    key = "x1 + x2 <= 450", explicit (450 - x1, x1, 0, xmax),
    color = red, key = "x2 <= 300", explicit (300, x1, 0, 350),
    color = brown, key = "4*x1 + 5*x2 <= 2000", explicit (400 - 4*x1/5, x1, 0, xmax),
    color = magenta, key = "x1 <= 350", parametric (350, yy,  yy, 0, 300),
color  = black, key = " z = 5000", 
         implicit (15*x1 + 25*x2 = 5000, x1, 0, xmax, x2, 0, ymax),
color = yellow, key = "z = 9000", 
         implicit (15*x1 + 25*x2 = 9000, x1, 0, xmax, x2, 0, ymax))$
    

We see from the two iso-profit lines that as we increase the profit z, the maximum value
inside the constraints is the corner intersection of the lines x2 = 300 and 
4*x1+5*x2 = 2000, which is (x1 = 125, x2 = 300) which gives 
profit = 15(125) + 25(300) = $9,375.

soln : solve ([x2 = 300, 4*x1 + 5*x2 = 2000], [x1, x2]);

at (15*x1 + 25*x2, soln[1]);

### Graphical Example 1, Maximize Profit, Scenario 2 LP

If we modify the above baby carrier project by increasing the sales price of the 
non-reversible baby carrier by $5, keeping all the same constraints, we want to 
maximize z = 20*x1 + 25*x2. We can reuse most of the above plot elements.

wxdraw2d (xlabel = "x1", ylabel = "x2", yrange = [0, ymax], 
    title = "z = profit = 20*x1 + 25*x2",
    key = "x1 + x2 <= 450", explicit (450 - x1, x1, 0, xmax),
    color = red, key = "x2 <= 300", explicit (300, x1, 0, 350),
    color = brown, key = "4*x1 + 5*x2 <= 2000", explicit (400 - 4*x1/5, x1, 0, xmax),
    color = magenta, key = "x1 <= 350", parametric (350, yy,  yy, 0, 300),
   color  = black, key = " z = 5000", 
            implicit (20*x1 + 25*x2 = 5000, x1, 0, xmax, x2, 0, ymax),
   color = yellow, key = "z = 9800", 
             implicit (20*x1 + 25*x2 = 9800, x1, 0, xmax, x2, 0, ymax))$

The yellow line is starting to leave the allowed solution region at the intersection of the
lines x1 + x2 = 450 and 4*x1 + 5*x2 = 2000.

soln : solve ([x1 + x2 = 450, 4*x1 + 5*x2 = 2000], [x1, x2]);

at (20*x1 + 25*x2, soln[1]);

The profit for the project has increased by $625.

### Graphical Example 2, Minimize Cost

A second Butenko Example (22:53 of Video 13) follows:

A retail store is planning an Ad campaign to increase customer traffic at both the online
web site and the physical store. Advertising options include magazine ads, with the 
estimate that $1,000 worth of magazine ads will yield 100 new visitors to the store and 
500 new visitors to the web site.

A second campaign option will be online advertising, with the projection that $1,000 
invested in online ads will yield 50 new store visitors and 1,000 new web site visitors.

The goals are to bring at least 500 new visitors to the physical store and at least 
5,000 new visitors to the web site.

Let x1 be the amount to be spent on ads in magazines (in units of $1,000 ) and likewise
x2 the amount to be spent on online ads. Then we want to minimize the advertising
expense while meeting the new visitor goals. Then with the data:
     if x1  = 1 --->100 store and 500 website visitors,
     if x2 = 1 ----> 50  store and 1000 website visitors,
we can set up the LP problem as:

Minimize  x1 + x2
  subject to
     100*x1 + 50*x2 >= 500,  (store goal)
     500*x1 + 1000*x2 >= 5000,  (website goal)
    with x1, x2 >= 0.

Dividing the first constraint by 50, and the second constraint by 500, the constraint 
conditions become
   2*x1 + x2 >= 10,
   x1 + 2*x2 >= 10

xmax : 10$
ymax : 10$
wxdraw2d (xlabel = "x1", ylabel = "x2", yrange = [0, ymax],
    key = "2*x1 + x2 >= 10", explicit (10 - 2*x1, x1, 0, xmax),
    color = red, key = "x1 + 2*x2 >= 10", explicit (5 - x1/2, x1, 0, xmax))$

Now add some iso-cost lines x1 + x2  = const.

wxdraw2d (xlabel = "x1", ylabel = "x2", yrange = [0, ymax],
    key = "2*x1 + x2 >= 10", explicit (10 - 2*x1, x1, 0, xmax),
    color = red, key = "x1 + 2*x2 >= 10", explicit (5 - x1/2, x1, 0, xmax),
    color = black, key  = "x1 + x2 = 13", implicit (x1 + x2 = 13,x1, 0, xmax,x2,0,ymax),
    color  = yellow, key = "x1 + x2 = 7", implicit (x1 + x2 = 7,x1, 0, xmax,x2,0,ymax))$

soln : solve ([2*x1 + x2 = 10, x1 + 2*x2 = 10], [x1, x2]);

at (x1 + x2, soln[1]);

float(%);

We can add the minimum cost line x1 + x2 = 20/3.

wxdraw2d (xlabel = "x1", ylabel = "x2", yrange = [0, ymax],
    key = "2*x1 + x2 >= 10", explicit (10 - 2*x1, x1, 0, xmax),
    color = red, key = "x1 + 2*x2 >= 10", explicit (5 - x1/2, x1, 0, xmax),
    color = black, key  = "x1 + x2 = 12", implicit (x1 + x2 = 12,x1, 0, xmax,x2,0,ymax),
    color  = yellow, key = "x1 + x2 = 7", implicit (x1 + x2 = 7,x1, 0, xmax,x2,0,ymax),
    color = dark_turquoise, key = "x1 + x2 = 20/3",
      implicit (x1 + x2 = 20/3,x1, 0, xmax, x2, 0, ymax),
    line_width = 1, color = black, key = "", explicit (10/3, x1, 0, 10/3),
    parametric (10/3, yy, yy, 0, 10/3))$

## Classification of LP's

An LP is called "feasible" if it has at least one feasible solution (agrees with all 
constraints) and "infeasible" otherwise.  (Older language is "unfeasible".)

An LP is called "optimal" if it has an optimal solution.

An LP is called "unbounded" it is feasible and its objective function is not bounded
(from above for a maximization problem and from below for a minimization problem) in
the feasible region (the region in which all the constraints and non-negativity conditions
are met).

## tableau (RowL), pivot1 (RowL, [row,col]), 
      simplex package

The Simplex method is an approach to solving linear programming models by hand 
using slack and/or excess variables, tableaus, and pivot elements as a means to finding 
the optimal solution of an optimization problem. 

A linear programming problem (LP problem) goal is achieving the best outcome 
given a linear objective function subject to linear constraint relations. 

Most LP problems can be solved using a solver such as Maxima's simplex package
functions maximize_lp (obj, condL, posL) or minimize_lp (obj, condL, posL). We will 
provide many examples of invoking both of these functions, as well as solving problems
 "by hand" via the simplex algorithm using our tableau and pivot1 functions defined in
 Econ1.mac.

### maximization example

A simple maximization LP example is:

maximize z = x1 + x2,
subject to:   x1 + 5*x2 <= 5,
                  2*x1 + x2 <= 4,
                       x1, x2 >= 0.

Both the objective function z and the constraint conditions are linear in the decision
variables x1 and x2. Both of the constraint inequalities have the form lhs <= pos-constant.
We require both of the decision variables x1 and x2 to be non-negative. Hence the
LP is already expressed in "standard form".

### Using maximize_lp (obj, condL, posL)

Using the Maxima simplex package function maximize_lp(obj, condL, posL), we get
(we loaded the simplex package at the top of this file):

maximize_lp ( x1 + x2, [x1  + 5*x2 <= 5, 2*x1 + x2 <= 4], [x1, x2]);

The returned solution is z = 7/3 using x1 = 5/3 and x2 = 2/3.

### "Hand Calculation" with Simplex Method & Solve

We convert <= constraints into equality constraints using two non-negative slack 
variables s1 and s2, 

maximize z = x1 + x2,
  subject to:
                   x1 + 5*x2 + s1         = 5,
                2*x1 +    x2         + s2 = 4, 
     with         x1, x2, s1, s2 >= 0.

Rewrite the constraint equations in the form:
    s1 = 5 - x1 - 5*x2,
    s2 = 4 - 2*x1 - x2.

z : x1 + x2;
eqn1 : s1 = 5 - x1 - 5*x2;
eqn2 : s2 = 4 - 2*x1 - x2;

We start with the basic feasible solution x1 = 0, x2 = 0 (z = 0), s1 = 5, s2 = 4.
Since z = x1 + x2, we can increase z by increasing x1 or x2 from zero to some
positive value as long as s1 and s2 remain non-negative. Since in the objective
expression (z) both x1 and x2 have the same coefficient (+1), we can increase z equally 
fast choosing either variable. If the coefficients were both negative, we could not 
increase z by increasing the decision variables x1 and x2 larger than zero. If only one
 coefficient was  positive, then that variable would be chosen for increase.
 
Let's arbitrarily choose x1, which then becomes the "pivot variable" and we say "x1 enters 
the basis". Basic variables are those which are greater than zero. With our basic feasible 
solution, s1 and s2 are basic variables (s1 = 5, s2 = 4) and x1 and x2 are non-basic 
variables (x1 = x2 = 0).

With x1 chosen as the pivot variable, we look at which of the constraint equations will 
limit the increase of x1 the most. In order for a constraint equation to limit the increase 
of x1, the coefficient of x1 in our original statement must have been positive, and after 
solving the constraint equations for the basic variables, as above, the coefficient of x1 
then is a negative number.

We then choose the "winning constraint equation" using the "minimum ratio test", which
involves the ratio: constant/coeff, which for our case depends of the ratios (5/1,  4/2).

Leaving x2 = 0, as we increase x1 we have s1 = 5 - x1, and s1 becomes equal to zero 
when x = 5 (ie., 5/1). We have s2 = 4 - 2*x1, which becomes zero when 4 - 2*x1 = 0, 
or x1 = 4/2 = 2. So we cannot satisfy both of the constraint equations (keeping both s1
and s2 non-negative) if we increase x1 by a value larger than 2. The constraint equation
which presents the tightest constraint is the "winning equation" and  as x1 --> 2, s1 --> 3
and s2 --> 0. We say "s2 leaves the basis". We solve the "winning constraint equation" 
for x1 as a function of s2:

solve (eqn2, x1), expand;

Use this result to define a "new eqn2":

eqn2 : %[1];

Use this new eqn2 to define a "new eqn1":

eqn1 : at (eqn1, eqn2), expand;

Since both x2 and s2 are non-basic variables, equal to zero, at the end of this first step,
x1 = 2 and s1 = 3.

Express the objective z in terms of non-basic variables x2 and s2 (both equal to zero)
by eliminating x1 (using our new eqn2) in favor of x2 and s2:

z : at (z, eqn2), expand;

Since x2 and s2 are equal to zero at the end of this step, we have increased z to the
 value 2. 

We are now ready for Step 2. Do any of the variables (x2 and s2) in z (now) have positive
coefficients? We see that only x2 has a positive coefficient (1/2), so we increase x2 from
its present zero value as large as possible (leaving s2 = 0) and without violating our new
constraint equations eqn1 and eqn2. Since s2 = 0, eqn1 says s1 = - 9*x2/2 + 3, so s1
will decrease from 3 to zero when x2 = 2/3.

Since x2 = 0, eqn2 says x1 = -x2/2 + 2 which equals zero when x2 = 2/(1/2) = 4. Hence 
the strongest constraint (minimum ratio) test is won by eqn1 (x2 increases to 2/3, s1 -> 0).
Hence x2 enters the basis with value 2/3 and s1 leaves the basis with value 0.

Use the winning eqn1 to solve for x2 as a function of s1 and s2.

solve (eqn1, x2), expand;

We replace the old eqn1 by the new eqn1, which expresses x2 as a function of s1
and s2.

eqn1 : %[1];

Since both s1 and s2 equal zero at this stage, x2 = 2/3.

Rewrite eqn2 by using our new eqn1 to express x2 as a function of s1 and s2.

eqn2;

eqn2 : at (eqn2, eqn1), expand;

Rewrite the objective z by using our new eqn1 to express z as a function of s1 and s2.

z;

z : at (z, eqn1), expand;

At the end of this Step, both s1 and s2 are non-basic variables and are both equal to 
zero, so we have increased z to 7/3.

Since both coefficients (of s1 and s2) are negative numbers, we cannot increase z by
increasing either s1 or s2 to a positive value (greater than zero). 

This "hand calculation" using solve thus has produced the same solution as 
as produced by maximum_lp: z = 7/3, x2 = 2/3, x1 = 5/3.

### "Hand Calculation" with Simplex Method & tableau (RowL)

Our tableau notation is quite similar to that used by Butenko (see References).

The "Simplex method" is a technique for solving linear programs by hand. To solve a 
linear programming model using the Simplex method the following steps are necessary: 

1. Standard form
2. Introducing slack and/or excess (surplus) variables
3. Creating the tableau
4. Choosing the Pivot element (and checking for optimality)
5. Creating a new tableau
6. Checking for optimality
7. Repeating Steps 4,5,6 until optimality is reached

Working this problem "by hand" with the simplex method, we need two non-negative 
slack variables s1, s2 to convert the constraint inequalities into constraint equations. 

If, instead, we have a problem in which at least one of the constraint inequalities has
the form  g(x1, x2) >= b1, with b1 > 0, we introduce a positive "excess variable"
(or "surplus variable") e1 > 0 by writing, for example,
   g (x1, x2) - e1 = b1,
with  x1, x2, e1, b1 >= 0.

In this simple introductory problem we only need slack variables s1 and s2 to write the
constraints as equalities. We then have the LP problem:

maximize z = x1 + x2
  subject to:
                   x1 + 5*x2 + s1         = 5,
                2*x1 +    x2         + s2 = 4,
              x1, x2, s1, s2 >= 0.

For Step 3 we create the tableau (French for "table") form. The virtue of the tableau form
 is  that attention is focused on the coefficients of the variables, their size and sign, 
which govern the crucial decisions: which variable is the pivot variable (entering the 
basis with a positive value) and which variable leaves the basis with value 0.
 
A Simplex tableau is used to perform row operations on the linear programming model 
as well as to check a solution for optimality. The tableau consists of the coefficients
corresponding to the linear constraint variables and the coefficients of the objective 
function. 

Butenko chooses the top numerical row (row 1) of the tableau to have coefficients which 
correspond to first rewriting the definition of z as  z + something = 0, or in our case, 
z - x1 - x2  + 0*s1 + 0*s2 = 0.

Thus:
z -    x1 -      x2  + 0*s1  + 0*s2    =  0 
        x1 +  5*x2 +    s1   + 0*s2    =  5,
     2*x1 +     x2 +  0*s1  +    s2    =  4,


In the definition of the tableau form, used by Sergiy Butenko's Operations Research
videos on YouTube (video 16), an extra column on the left side is included, labeled as 
the z column of coefficients, which we show using the Maxima matrix function:

In the tableau below, the top row of the tableau states what each column represents. 
More precisely, the x1 column of numbers are the coefficients of x1 in the numerical
objective equation (z - x1 - x2 = 0) and in the numerical constraint equations, and 
likewise for the x2, s1, s2 columns.

kill(z)$
matrix ( [z, x1, x2, s1, s2, rhs, basis],
            [1, -1, -1, 0, 0, 0, z],
            [0, 1, 5, 1, 0, 5, s1],
            [0, 2, 1, 0, 1, 4, s2] );

The first column ends up playing no role in the solution, so our tableau function ignores
 that first column in applying the simplex method of solution and our future tableaus will
 ignore that "z column" of Butenko.

We use our tableau (RowL) function which uses global lists vL and bL,  the first  list (vL)
holding the column headings, and the second list bL holding [z, s1, s2,...], in which all 
except z are the current basic variables.

For our example, RowL = [r1, r2, r3], where r1, r2, r3 are lists of coefficients reading left to 
right across the table. Our convention is that r1 is the list of the z eqn coefficients (after 
bringing over all terms to the left hand side). r2 is a list of the coefficients of the variables 
in the first constraint equation, etc. We often call row 1 (r1) the z-row.

kill (z, x1, x2, s1, s2)$
vL : [x1, x2, s1, s2, rhs, Basis];
bL : [z, s1, s2];
r1 : [ -1, -1, 0, 0, 0];
r2 : [ 1, 5, 1, 0, 5];
r3 : [ 2, 1, 0, 1, 4];

tableau ([r1, r2, r3]);

Our initial feasible solution, corresponding to x1 = x2 = 0, z = 0, s1 = 5, s2 = 4 is 
represented by these three row lists, r1, r2, r3, of coefficients. To be a valid starting 
tableau for the Simplex algorithm, the basic variables (here s1 and s2) must have zero 
coefficients in the z-row (r1), which is true here, and as a consequence the s1 column
 of coefficients has a value of 1 in the s1 row, with the remaining s1 column coefficients
 equal to 0. Likewise for the s2 column of coefficients. x1 and x2 are non-basic variables
 with value 0.

Step 4: Identify the Pivot Element in the tableau.

The pivot element is at the intersection of the pivot column and the pivot row.

First choose the pivot column by finding the most negative number in the top row of
numbers.  Since z = x1 + x2 is our initial definition of z, we can increase z by increasing 
either x1 or x2 to a value greater than zero, leaving the other decision variable equal to
zero. 

Using the tableau form, we look in the numerical top row (row 1) for the largest 
NEGATIVE coefficient to decide which of the two variables (x1, x2) to increase from 
value 0, leaving the other at value 0. 

If there are no negative coefficients in the objective row of the tableau, we cannot 
increase z from the ifs z = 0, x1 = 0, x2 = 0, and we have found the optimal solution.

Since both x1 and x2 have the same positive coeffient (+1), we arbitrarily start by 
increasing x1, leaving x2 = 0. The "pivot column" is the x1 column. x1 "enters the 
basis" with a positive value.

To determine the "pivot row", we use the "minimum ratio test".

Looking at the x1 column of coefficients, we look only at the constraint rows. A candidate 
for a pivot row must contain a positive x1 coefficient. In our example, both x1 column 
constraint row coefficients (1 and 2) are positive, so we choose between these two rows 
using our min. ratio test. For row 2, we divide the rhs constant by 1, ratio  = 5/1 = 5. For 
row 3, we divide the rhs by 2, ratio = 4/2 = 2. Since 2 < 5, row 3 is the winner of the test 
and becomes the "pivot row".  Since row 3 is the "s2 row", s2 leaves the basis with 
value 0.

Hence the pivot element is row 3, column 1, represented by the list [3,1], and has
the value 2.

You can use our function tratio (rowL, ncol) to see the ratios calculated for you.

tratio ([r1, r2, r3], 1);

Step 5. Create the New Tableau

We again display the initial feasible solution (ifs) tableau of coefficients.

tableau ([r1,r2,r3]);

Applying the simplex algorithm, we divide the pivot row (row 3, r3) by 2, so the pivot 
element becomes equal to 1. We *redefine* the basic label list bL so row 3 becomes 
the x1 row (as we interchange the roles of s2 and x1).

r3 : r3/2;
bL : [z, s1, x1];
tableau ([r1,r2,r3]);

The next step of the simplex algorithm is to add multiples of the pivot row (the third row)
to the other rows so that all other coefficients in the pivot column (x1) equal zero. We 
can accomplish this step by adding numerical row 3 to numerical row 1 and adding
(-1)*row3 to row2. At the completion of these steps, we say that "x1 has entered the
basis" (has a value greater than zero), and "s2 leaves the basis" (s2 becomes equal
to zero), and we have indicated this in the tableau basis column (the right most column).

r1 : r1 + r3;
r2 : r2 - r3;

tableau ([r1, r2, r3]);

Step 6. Check for Optimality

We still have a negative coefficient in the objective row (the z-row) in the x2 column, 
so we can increase x2 from its present value 0 to a value greater than zero to increase 
the value of z above its present value (z = 2). The x2 column becomes the pivot column
(x2 becomes the "pivot variable").

The first numerical row says
  z + 0*x1 - x2/2 + 0*s1 + s2/2 = 2
but x2 = 0, and s2 = 0 (s2 left the basis), so this reduces to z = 2, which can be
read off more quickly from the values in the rhs column.

The third numerical row says
  x1 + x2/2 + 0*s1 + s2/2 = 2,
but with x2 and s2 both equal to zero, this reduces to x1 = 2.

Likewise, the second numerical (first constraint) row says s1 = 3.

tratio ([r1,r2,r3], 2);

Step 7. Repeat Steps 4, 5, 6.

We have two positive candidates (9/2, 1/2) in the x2 column for the pivot element. 
Applying the minimum ratio test, 3/(9/2) = 2/3 and  2/(1/2) = 4, so row 2 (the s1 row) is 
the winner, and s1 leaves the basis with value 0 while x2 enters the basis with a 
positive value.

The pivot element is numerical row 2 and column 2, with value 9/2. 

We multiply numerical row 2 by 2/9 to get the pivot element equal to 1. We switch out
the s1 label for the x2 label in the basis label list bL.

r2 : (2/9)*r2;
bL : [z, x2, x1];
tableau ([r1, r2, r3]);

We then add 1/2 of the pivot row (row 2, r2)  to the first row (z-row, r1), and add
(-1/2) of the pivot row to the third row to get all elements of the x2 column equal to
zero except the pivot element. This results in (x2 enters and s1 leaves the basis):

r1 : r1 + r2/2;
r3 : r3 - r2/2;

tableau ([r1,r2,r3]);

There are no more negative coefficients in the objective row (the z-row, r1), so the new 
tableau represents the optimal solution: x1 = 5/3, x2 = 2/3, z = 7/3 = 2.33. The second 
pivot operation increased z by 0.33.

This "hand calculation" using the simplex algorithm result agrees with the return 
values from maximize_lp above.

### Using RL : pivot1(RL, [nrow, ncol]), tratio(RL, ncol)

The Maxima function  pivot1 (RowL, [nrow, ncol]), defined in Econ1.mac, shows 
the result of one pivot operation based on a selected pivot element from the constraint 
rows. nrow is the pivot row number (greater than 1 since row 1, r1, the z-row, cannot be
the pivot row).

The lists vL and bL need to be defined prior to using pivot1 for a succession of
simplex pivot operations. Because at this point the symbols x1, x2, s1, s2, rhs, Basis, 
and z are not bound to any other expression or value, we can omit including double 
quotes ("), as in "x1". ("x1" is a Maxima string.)
 
Both of the lists vL and bL are then global lists. The function pivot1 preserves vL and 
changes bL each step as a basic variable is swapped with a non-basic variable. 

We go back to the beginning of the previous problem, so we need to redefine vL, 
bL, r1, r2, r3 corresponding to our initial feasible solution tableau, in which s1 and s2
are the basic variables (greater than zero), and x1 and x2 are the non-basic variables 
and hence equal to zero.

Step 0 tableau; initial feasible solution (ifs):

vL : [x1, x2, s1, s2, rhs, Basis]$
bL : [z, s1, s2]$
r1 : [ -1, -1, 0, 0, 0]$
r2 : [ 1, 5, 1, 0, 5]$
r3 : [ 2, 1, 0, 1, 4]$
RL : [r1, r2, r3]$
tableau ( RL);

tratio ( RL, 1);

Our pivot column is the x1 column, our pivot row is the s2 row, using the minimum ratio 
test, (5/1 = 5, 4/2 = 2),  so our pivot element is defined by the [nrow, ncol] = [3,1]. 
(Our row numbers refer only to the set of numerical rows). The value of the pivot 
element is 2.

Here we define the lists vL and bL. 

Here we call pivot1 (RowL, [row, col]) using the tableau rows defined via r1, r2,
and r3 assignments above, (RowL = [r1,r2,r3]) and implicitly using the global lists vL 
and bL. pivot1 then carries out the first Simplex algorithm step, which involves dividing
 the pivot row by the value of the pivot element, and then adding multiples of the pivot
row to the other constraint rows in such a way that all elements of the pivot variables 
column equal zero, except for the pivot element.

The pivot element list has the form  [nrow, ncol], which here is [3, 1], standing for row 3
and column 1.

Step 1 tableau:

RL : pivot1 ([r1,r2,r3],[3,1])$

We have redefined RL as the output list produced by pivot1 - this new RL can then
be used for further calls to tratio and pivot1

 Using the minimum ratio test (3/(9/2) = 2/3, 2/(1/2) = 4,  row 2 is the winner, x2 enters 
the basis with a positive value, s1 leaves the basis with value 0.

tratio(RL, 2);

Step 2 tableau:

RL : pivot1(RL, [2,2] )$

This Step 2 tableau has no negative coefficients in the objective row, so optimization
has been achieved with the variable values and z values indicated in the two right
hand columns: z = 7/3, x2 = 2/3, x1 = 5/3, agreeing with maximize_lp.

Using our Maxima function pivot1 thus can reduce your work in hand calculations. You
still need to identify the most negative coefficient in the objective row to determine the
pivot column, and the constraint row which produces the minimum ratio (using positive 
pivot column coefficients) to select the pivot row, and hence the pivot element [row,col]
at each step. You can use tratio (RL, ncol) to help with the arithmetic.

### Minimization Problem using the tableau, pivot1 method

Here is an minimization LP example:

Minimize  w = -8*x1 - 10*x2 - 7*x3
subject to
   x1 + 3*x2 + 2*x3 <= 10,
  -x1 - 5*x2 - x3 >= -8,
   x1, x2, x3 >= 0. 

### Using mimimize_lp (obj, condL, posL)

minimize_lp (-8*x1 - 10*x2 - 7*x3, [x1 + 3*x2 + 2*x3 <= 10, -x1 - 5*x2 - x3 >= -8],
    [x1, x2, x3] );

### Converting to a Maximization Problem

"Hand calculation" using the simplex method involves several steps:

Step 1: Convert problem to "standard form", which means
  (1) Must be written as a maximization problem,
  (2) All constraint conditions must be either (for decision variables = x1,x2,x3 case)
           f(x1,x2,x3) <= positive-constant,  or
           f(x1,x2,x3) >= positive-constant.
  (3)  All decision variables must be non-negative.

To transform a minimization LP into a maximization LP, multiply both sides of
the objective function by (-1). Thus
   Minimize  w = -z  = -8*x1 - 10*x2 - 7*x3 is equivalent to
    Maximize (-1) * [ -z = -8*x1 - 10*x2 - 7*x3],  
   or:    Maximize   z = 8*x1 + 10*x2 + 7*x3.

To transform g(x1,x2,x3 >= -A, where A > 0, into a type <= condition, multiply both 
sides of the given condition by (-1):
  (-1)* [ -x1 - 5*x2 - x3 >= -8 ] then becomes
     x1 + 5*x2 + x3 <= 8.

### Standard Form

We now have the maximization linear problem (LP) in the "standard form":

  maximize z = 8*x1 + 10*x2 + 7*x3
  subject to:
    x1 + 3*x2 + 2*x3 <= 10,
   x1 + 5*x2 + x3 <= 8,
   x1,x2,x3 >= 0.

### Using maximize_lp (obj, condL, posL)

maximize_lp (8*x1 + 10*x2 + 7*x3, [x1 + 3*x2 + 2*x3 <= 10, x1 + 5*x2 + x3 <= 8],
                       [x1, x2, x3] );

Since the maximum of z is 64, the minimum of (-z) is -64.

### Introduce Slack Variables

Step 2: Introduce slack variables and/or excess variables.

Slack and excess variables are additional positive variables that are introduced into the
linear constraints of a linear program to transform them from inequality constraints into
equality constraints. 

If the model is in standard form, the slack variables will always have a coefficient of  (+1)
and the excess variables (surplus variables) will have a coefficient of (-1).

Slack and/or excess variables are needed in the constraints to transform them into
solvable equalities with one definite answer. In this example, we need two slack variables
and we call the slack variables s1 and s2.

maximize z = 8*x1 + 10*x2 + 7*x3
  subject to:
    x1 + 3*x2 + 2*x3  + s1 = 10,
   x1 + 5*x2 + x3  + s2 = 8,
   x1, x2, x3, s1, s2 >= 0.

### Using tableau(RowL) and pivot1(RowL, [row,col])

Step 3: Setting up the Tableau

A Simplex tableau is used to perform row operations on the linear programming model 
as well as to check a solution for optimality. The tableau consists of the coefficients
corresponding to the linear constraint variables and the coefficients of the objective 
function. In the tableau below, the top row of the tableau states what each 
column represents. 

The first row (of numbers) in the tableau represents the objective function decision 
variable coefficients after the decision variables are brought over the the left hand side,
 as in    z - 8*x1 - 10*x2 - 7*x3 = 0.
 
Then we have two rows to represent the linear constraint variable coefficients from the 
linear programming model. 

vL : [x1,x2,x3,s1,s2,rhs,Basis]$
bL : [z,s1,s2]$
r1 : [-8,-10,-7,0,0,0]$
r2 : [1,3,2,1,0,10]$
r3 : [1,5,1,0,1,8]$
RL : [r1, r2, r3]$
tableau (RL);

This tableau represents an initial feasible solution (ifs) with x1, x2, x3 non-basic variables
with value 0, and with basic variables s1 = 10 and s2 = 8.

The most negative coefficient in the z-row (Row 1, r1) is (-10) so x2 is the pivot variable,
x2 enters the basis with a positive value. To find the pivot row we use the minimum ratio
test for x2 column coefficients in the constraint row which are positive: (10/3 = 3.33,
8/5 = 1.6) so the s2 row is the winner, s2 leaves the basis with value 0, and the pivot
element is 5 in row3, col 2.

tratio(RL, 2);

For our first use of pivot1 on a specific problem, we need to enter the rows r1,r2,r3 the
same way we entered them using tableau. pivot1 then takes one step based on the 
chosen pivot element  [row, col].

Step 1 tableau:

RL : pivot1 ( RL, [3, 2])$

z has increased from 0 to 16 in the first step. The solution is not optimal, since there are
negative coefficients in the z-row, the most negative being (-6) so x1 is the pivot variable,
x1 enters the basis with a positive value. There are two positive coefficients in the
constraint rows of the x1 column (2/5, 1/5). Using the minimum ratio test to select the
pivot row, the ratios are  (26/5)/(2/5) = 26/2 = 13 and (8/5)/(1/5) = 8 = min, so row 3 
is the winner, x2 leaves the basis with value 0. The pivot element is [3,1] with value 1/5.

tratio (RL, 1);

For our second step using pivot1, we use as input row list RL the output of the 
previous step.

RL : pivot1 (RL, [3, 1])$

There are no more negative coefficients in the z-row, so the optimum solution has been
reached with x1 = 8, x2 = 0, x3 = 0,  and z = 64.

Since the maximum of z is 64, the minimum of w = -z is -64.

## Degeneracy and Cycling

### A Cycling Example

We follow Butenko, Video 17 here. Consider the LP:
Maximize z = 5*x1 + 4*x2 - 20*x3 - 2*x4
subject to:
  x1/4 - x2/8 + 12*x3 + 10*x4 <= 0,
  x1/10 + x2/20 + x3/20 + x4/5 <= 0,
  with x1, x2, x3, x4 >= 0.

maximize_lp (5*x1 + 4*x2 - 20*x3 - 2*x4, [x1/4 - x2/8 + 12*x3 + 10*x4 <= 0,
            x1/10 + x2/20 + x3/20 + x4/5 <= 0], [ x1, x2, x3, x4] );

The simplex package function maximize_lp is unable to calculate an optimum and
returns the initial feasible solution with all decision variables equal to zero and z = 0.

If we introduce non-negative slack variables x5 and x6 in the two constraint conditions, 
we construct the Step 0 tableau:

vL : [x1, x2, x3, x4, x5, x6, rhs, Basis]$
bL : [z, x5, x6 ]$
r1 : [-5,-4, 20, 2, 0, 0, 0]$
r2 : [1/4, -1/8,12, 10, 1, 0, 0]$
r3 : [1/10,1/20,1/20,1/5, 0, 1, 0]$
RL : [r1, r2, r3]$
tableau (RL);

A peculiar feature of this tableau is that all of the right-hand sides are equal to zero.
The tableau represents an initial feasible solution in which the four input variables are all
equal to zero, but also the basic variables x5 and x6 are equal to zero.

The most negative coefficient in the z-row is -5, so we let x1 be the pivot variable, x1 
enters the basis, and we would normally use the minimum ratio test to choose the pivot 
row. Since 0/(1/4) = 0 and 0/(1/10) = 0, there is no "winner", so we arbitrarily choose 
row 2 with the pivot element [2, 1] with value 1/4.

Next is the Step 1 tableau:

RL : pivot1 ( RL, [2, 1])$

x1 has entered the basis, but with a zero value, x5 has left the basis with a zero value, 
the new solution has all the variables equal to zero and z = 0, which is the same 
solution we started with.

An iteration of the Simplex method, which results in a new basis with a solution 
identical to the previous feasible solution is called a "degenerate iteration"; the 
phenomenon is referred to as "degeneracy".

Hopefully, a few of these degenerate iterations will be followed by an iteration which is 
not degenerate. If we are unable to get out of this degenerate trap, the result is 
called "cycling".

The only negative coefficient in the z-row is -13/2 so x2 enters the basis as the pivot 
variable, and x6 leaves the basis with a zero value since the only positive coefficient in 
the x2 column of the constraint rows is 1/10.

Next is the Step 2 tableau:

RL : pivot1(RL, [3, 2])$

Next is the Step 3 tableau:

RL :  pivot1 (RL, [2, 3])$

tratio (RL, 4);

Next is the Step 4 tableau:

RL : pivot1 (RL, [3, 4])$

Next is the Step 5 tableau:

RL : pivot1 (RL, [2, 5])$

Next is the Step 6 tableau:

RL : pivot1 (RL, [3,6])$

The Step 6 tableau is identical to the Step 0 tableau we started with. Each of the 
iterations is degenerate. Further iterations would simply cycle through the same set 
of tableaus with no  solution which improves on z = 0. 

### Avoiding Cycling using Bland's Rule

Several methods are used to avoid cycliing.

Bland's Rule:

Let the variables be ordered in some way which reflects an ordering, for example,
x1, x2, x3, x4, x5, ...
Then, when there are several candidates for the choice of the entering variable (the pivot
variable) or the leaving variable, choose the variable that appears earlier in the ordering 
method.

If Bland's Rule is used to select the entering and leaving variables in the simplex 
algorithm, then cycling never occurs.

In our cycling example, we followed Bland's rule with the Step 0 tableau: both the x1 and
x2 variables are candidates for the pivot variable with z-row coefficients -5 and -4 
respectively. We chose x1. Then there were two candidates for the pivot row, the x5 
row and the x6 row, and we arbitrarily chose the x5 row.

In the resulting Step 1 tableau there were no choices available for the entering variable
other than x2 and no choices for the leaving variable other than x6.

In the resulting Step 2 tableau we chose x3 as the entering variable rather than x4 or x5 
and there was no choice in the leaving variable.

In the resulting Step 3 tableau we chose x4 as the entering variable rather than x5, and 
chose x2 as the leaving variable rather than x3.

In the resulting Step 4 tableau there was no choice in the entering variable and no 
choice in the leaving variable.

In the resulting Step 5 tableau we chose the most negative coefficient in the z-row (-10) 
to make x6 the entering variable, with the alternative being x1 (-6) and x2 (-9/2), which 
did not follow Bland's rule.

If we follow Bland's rule with the Step 5 tableau, what happens?

We need to define bL, r1, r2, r3 to imply the Step 5 tableau.

bL : [z, x5, x4 ]$
r1 : [-6,-9/2, 39/2, 0, 0,-10, 0]$
r2 : [-19/4, -21/8, 19/2, 0, 1, -50, 0]$
r3 : [1/2,1/4,1/4,1, 0, 5, 0]$
RL : [r1,r2,r3]$
tableau (RL);

We follow Bland's rule and choose x1 as the entering variable. There is no choice in the
leaving variable.

RL : pivot1 (RL, [3, 1])$

There is only one negative coefficient in the z-row so no choice in the entering variable.
There is only one positive coefficient in the constraint rows of the entering variable, so 
again there is no choice on choosing the leaving variable.

RL : pivot1 (RL, [3, 2])$

The solution has not changed compared to the initial feasible solution at Step 0, but this
last tableau proves its optimality, since there are no negative coefficients in the z-row, 
which implies an optimal solution. The optimal solution is x1 = x2 = x3 = x4 = 0, z = 0.

### Potential Issues with the Simplex Algorithm

Initialization: Can we always find an initial starting solution?

Iteration: Can we always choose an entering variable, a leaving variable, and construct
 the next tableau using pivoting? 

If there are no negative coefficients in the z-row, then there is no way we can increase 
the objective z. Given an entering variable, if we cannot find a leaving variable, then 
there is no solution since we have an unbounded LP.

Termination: Is it possible that an endless sequence of solutions will be generated and
 optimal tableau is never found?

If the simplex method avoids cycling, it must terminate by either finding an optimal
solution or else finding that the LP is unbounded.

## Artificial Variables, Two-Phase Simplex Method

We follow Butenko's example in his YouTube video 20: 
   The Simplex Method Part III, video 2 -
            The two phase simplex method and the big M method.

The example considered is:

Maximize z = x1 - 2*x2 + 3*x3,
subject to
   - 2*x1 + 3*x2 + 4*x3 >= 12,
    3*x1 + 2*x2 + x3 >= 6,
    x1 + x2 + x3 <= 9,
   with x1, x2, x3 >= 0.

maximize_lp (x1 - 2*x2 + 3*x3, [- 2*x1 + 3*x2 + 4*x3 >= 12,3*x1 + 2*x2 + x3 >= 6,
                                                   x1 + x2 + x3 <= 9], [x1, x2, x3]);

We achieve the optimum value of z = 27 by using x1 = 0, x2 = 0, and x3 = 9.

### Reduction to Standard Form with Artificial Variables

When do we need Artificial Variables?

 Introduce an artificial variable in the cases:

    1. A constraint equation which is an equality constraint (without a slack or excess
 variable) so that we are sure to have a starting nonnegative basic variable 
corresponding to that given original constraint equation,

    2.  A <= type constraint in which the right-hand side is negative, since a slack 
variable would start with a negative value when the input variables start at value zero. 
The slack variable then cannot be part of the initial basis, and we need to introduce 
also an artificial variable. 
With the steps:  f(x1,x2,x3) <= -10 =>  f(x1,x2,x3) + s = -10  => f(x1,x2,x3) + s - a = -10,
  we can start with  x1 = x2 = x3 = s = 0 and a = 10.

    3.  A  >= type constraint with a positive right-hand side. Again the goal is to
achieve a standard form equality constraint in which we have a basic variable which is
nonnegative when the input variables start at a zero value.  With the steps
f(x1,x2,x3) >= 10  =>   f(x1,x2,x3) - e = 10  =>  f(x1,x2,x3) - e + a = 10, we can start with
 x1 = x2 = x3 = e = 0 and a = 10.

In our example, we need to add two excess variables (e1 and e2) to convert the first 
two conditions into equations and add one slack variable (s3) to convert the third 
condition into an equation.
    - 2*x1 + 3*x2 + 4*x3 - e1 = 12,
    3*x1 + 2*x2 + x3 - e2 = 6,
    x1 + x2 + x3 + s3 = 9,
    with x1, x2, x3, e1, e2, s3 >= 0.

We then need to add artificial variables a1 and a2 to the left hand side of the first two
constraint equations to arrive at

Maximize z = x1 - 2*x2 + 3*x3,
 subject to
 - 2*x1 + 3*x2 + 4*x3 - e1 + a1 = 12,
    3*x1 + 2*x2 + x3 - e2 + a2 = 6,
    x1 + x2 + x3 + s3 = 9,
    with x1, x2, x3, e1, e2, s3, a1, a2 >= 0.

### Phase I

Phase I determines if the given LP has a feasible solution by solving the
"auxiliary problem" (A).

The objective of the auxiliary problem (A) is to minimize the sum of the artificial variables.
We want to drive them all down to the value zero, to make sure we can get rid of them.

1. If the optimal objective value of (A) is 0, we can conclude that the LP problem is
feasible. One needs to make sure that all the artificial variables are non-basic (ie.,
equal zero) and then one can proceed to Phase II.

2. If the optimal objective value of (A) is NOT equal to zero, we can conclude that
the given LP problem is not feasible, and we are done.

Using minimize_lp (a1 + a2, condL, posL):

We can use minimize_lp (a1 + a2, condL, posL) to establish feasibility of the given
example LP:

minimize_lp (a1 + a2, [- 2*x1 + 3*x2 + 4*x3 -e1 + a1 = 12,
                                    3*x1 + 2*x2 + x3 - e2 + a2 = 6,
                                 x1 + x2 + x3 + s3 = 9], 
                               [x1, x2, x3, e1, e2, s3, a1, a2]);

Since the optimum value is zero for the sum of the artificial variables, this establishes
that the given LP is feasible.

Initial Tableau for Phase I:

Working Phase I "by hand", we convert the Phase I problem from
   minimize (a1 + a2)
to
   maximize (- a1 - a2),
so letting z again represent the Phase I objective, 
  maximize   z  = - a1 - a2
subject to
   - 2*x1 + 3*x2 + 4*x3 -e1 + a1 = 12,
  3*x1 + 2*x2 + x3 - e2 + a2 = 6,
      x1 + x2 + x3 + s3 = 9
with x1, x2, x3, e1, e2, s3, a1, a2 >= 0.

Construct the initial tableau for Phase I, and in standard form

  z  + 0*x1 + 0*x2 + 0*x3 + 0*e1 + 0*e2 + 0*s3  +  a1  +    a2  = 0, 
      - 2*x1 + 3*x2 + 4*x3  -    e1 +  0*e2 + 0*s3 +   a1  + 0*a2 = 12,
        3*x1 + 2*x2 +    x3  + 0*e1     - e2 + 0*s3 + 0*a1  + a2   =   6,
           x1 +    x2 +    x3  +  0*e1 + 0*e2 +   s3  + 0*a1 + 0*a2 =  9

The 9 numerical tableau columns are:
x1, x2, x3, e1, e2, s3, a1, a2, rhs.

vL : [x1, x2, x3, e1, e2, s3, a1, a2, rhs, Basis]$
bL : [z, a1, a2, s3]$
r1 : [0, 0, 0, 0, 0, 0, 1, 1, 0]$
r2 : [-2, 3, 4, -1, 0, 0, 1, 0, 12]$
r3 : [3, 2, 1, 0, -1, 0, 0, 1, 6]$
r4 : [1, 1, 1, 0, 0, 1, 0, 0, 9]$
tableau ([r1,r2, r3, r4]);

Step 0 Tableau for Phase I:

The basic variables are a1, a2, and s3. A general feature of the Simplex method is that
the z-row coefficients of the basic variables must be zero. If the add (-1)*row2 to row 1 
and add (-1)*row3 to row 1 our final z-row will have zero coefficients in z-row.
We can do this in one step.

r1 : r1 - r2 - r3$
RL : [r1,r2,r3,r4]$
tableau (RL);

which is Step 0 Tableau for Phase I.

Step 1 Tableau for Phase I:

The coefficients of x2 and x3 in the z-row are equally negative (-5) and we arbitrarily
choose x2 as the pivot variable (x2 enters the basis) and look for the appropriate pivot
row by using the minimum ratio test for the positive x2 coefficients in the three constraint
rows.

tratio (RL, 2);

For row2, 12/3 = 4, for row3, 6/2 = 3, for row4, 9/1 = 9, so row3 has the minimum
ratio and becomes the pivot row, (a2 leaves the basis - becomes equal to zero). 
When using pivot1(RowList, [row,col]), our pivot element in then [3, 2].

Our first use of pivot1 produces the Step 1 Tableau. Recall that the first step of the 
Simplex algorithm is to divide the pivot row by the value of the pivot element (2) to
achieve a pivot element value equal to 1, and then add multiples of the pivot row to
the other rows to result in all other pivot column elements to becoming equal to zero. 
pivot1 shows the overall result of the pivot operation step.

RL : pivot1 ( RL, [3, 2] )$

Step 2 Tableau for Phase I:

The z-row (r1) still has two negative coefficients, so the Step 1 tableau is not optimal. 
The  most negative coefficient is (-5/2) in column 3, x3 is our pivot variable, col 3 is 
our pivot column.

tratio (RL, 3);

 
All the column 3 constraint row elements are positive and we use the minimum ratio 
test to find row 2 has the minimum ratio (1.2), so row 2 is the pivot row, the pivot 
element is [2, 3] with the value 5/2. We can then find the Step 2 tableau.

RL : pivot1 (RL, [2, 3] )$

There are no negative coefficients in the z-row, so we cannot drive our objective function
any larger than zero, and the Step 2 tableau is the optimal solution for Phase I, in which 
both a1 and a2 are zero, and our objective function z = -a1 -a2 = 0. Recall the auxiliary
problem was to minimize (a1 + a2) which we translated into maximize (-a1 - a2) for
application of the simplex algorithm.

Since the optimal objective value of the auxiliary problem (A) is 0, we can conclude that 
the given LP problem is feasible.  We can proceed  to Phase II.

### Phase II

Phase II is performed only if Phase I concludes the given LP is feasible. Phase II either
finds an optimal solution or finds that the solution is unbounded.

1. Remove the artificial variable columns from the optimal tableau of (A), and replace the
z-row of the resulting tableau with the objective of the original LP.

2. Eliminate all nonzero basic coefficients in the z-row to determine the initial feasible
solution tableau for application of the Simplex algorithm.

3. Apply the Simplex algorithm to either find the optimal solution or conclude the LP
problem is unbounded.

Initial Tableau for Phase II:

The current list RL defines the rows of the final tableau of Phase I. 

We remove the 7th and 8th column (the coefficients of a1 and a2) from the tableau, 
which means the 7th and 8th elements of each of our four rows and our column
label list vL. Each row will then contain six coefficients, corresponding to the
variables x1,x2,x3,e1,e2,s3, and the last (and seventh) element is the "rhs" value.

For row1 (the "z-row") we return to our original objective definition z = x1 - 2*x2 + 3*x3,
or z - x1 + 2*x2  - 3*x3 = 0. Redefining r1 we have coefficients equal to zero for e1, e2, 
s3, and rhs = 0. This leaves a seven element list.

r1 : [-1, 2, -3, 0, 0, 0, 0];

For vL and the constraint rows r2, r3, r4 we can use our function
      remL (aL, vL) 
to remove from list aL the elements whose numbers are in the list vL. For example, 
remL (aL, [7, 8]) returns a depleted list with elements seven and eight missing. 
remL is defined in Econ1.mac.

vL : remL (vL, [7, 8]);

r2 : remL (RL[2], [7, 8] );

r3 : remL (RL[3], [7, 8] );

r4 : remL (RL[4], [7, 8] );

We leave the list bL which as the Basis labels unchanged from the last tableau of
Phase I:

bL;

tableau ( [r1, r2, r3, r4] );

Step 0 Tableau for Phase II:

To arrive at Step 0 tableau for Phase II, we need to make sure that all basic variables
 (nonzero variables) have zero coefficients in the z-row (row 1), since  basic variables
 have all zeros in their column, except for a 1 in their row. As an example, see the 
column and row corresponding to the variable s3. Since x2 and x3 are basic variables
 in this tableau, we redefine r1 by adding -2*r3  and 3*r2 to r1.

r1 : r1 - 2*r3 + 3*r2$
RL : [r1,r2, r3, r4]$
tableau (RL);

Step 1 Tableau for Phase II:

The  most negative coefficient in the z-row is - 72/5, so x1 is the pivot variables, x1 enters
the basis. Looking at the positive coefficients in the x1 column,  (12/5)/(14/5) = 12/14, 
and (27/5)/(4/5) = 27/4, so row 3 wins the minimum ratio test, and x2 leaves the basis.
The pivot element is [3,1] with value 14/5. We can then compute the Step 1 tableau.

tratio (RL, 1);

RL : pivot1 ( RL, [3, 1] )$

The most negative z-row coefficient is -5/7 in the fifth column, so e2 enters the Basis.
The only positive e2 column constraint row coefficient is 3/7 in the fourth row.

tratio (RL, 5);

Step 2 Tableau for Phase II:

RL : pivot1 (RL, [4, 5])$

The only negative coefficient in the z-row is - 1/3 in the fourth column, so e1 enters the
Basis.

tratio (RL, 4);

Step 3 Tableau for Phase II:

RL : pivot1 (RL, [3, 4])$

The Step 3 tableau is optimum since there are no negative coefficients in the z-row.
The optimum solution is z = 27 with x1 = x2 = 0, and x3 = 9, in agreement with
the maximize_lp result.

## Artificial Variables, Big M Method

In Butenko's Simplex Part III, Video 2 on Youtube, (10:55) the same problem is solved 
using the Big M method, instead of the two phase method.

The example considered is (same as with the two phase method):

Maximize z = x1 - 2*x2 + 3*x3,
subject to
   - 2*x1 + 3*x2 + 4*x3 >= 12,
    3*x1 + 2*x2 + x3 >= 6,
    x1 + x2 + x3 <= 9,
   with x1, x2, x3 >= 0.

Refer to Setting up the Auxiliary Problem, When Do We Need Artificial Variables section
above. In the Big M method, since we have to introduce two artificial variables a1 and 
a2 in the constraints, we rephrase the problem as:

Maximize  x1 - 2*x2 + 3*x3 - M*a1 - M*a2,
 
in which instead of minimizing (a1 + a2) explicitly, we subtract a large (M) multiple of 
the values of a1 and a2 from the original objective as a penalty and an impediment to
maximization. The artificial variables vanish from the problem whenever the original LP
is feasible. We essentially treat M as an arbitrarily large positive number at each step.
 

Bringing all terms to the left, and with x1,x2,x3,e1,e2,s3,a1,a2 all non-negative, and
M a large positive number (larger than any other number involved in the problem), 
the standard form is:

z  -   x1 + 2*x2 -  3*x3  + 0*e1 + 0*e2 + 0*s3  +  M*a1 +  M*a2  =  0,
   - 2*x1 + 3*x2 + 4*x3  -  e1   +  0*e2 + 0*s3  +     a1  +   0*a2  = 12,
     3*x1 + 2*x2 +    x3  + 0*e1 -     e2 + 0*s3  +  0*a1  +      a2   =  6,
        x1 +    x2 +    x3  + 0*e1 + 0*e2 +    s3   +  0*a1 +   0*a2   =  9

Initial Tableau:

We next put the standard form in the tableau format, redefining our starting row r1,
with constraint rows r2, r3, r4 the same as Phase I initial tableau of the two phase 
method.

vL : [x1, x2, x3, e1, e2, s3, a1, a2, rhs, Basis]$
bL : [z, a1, a2, s3]$
r1 : [-1, 2, -3, 0, 0, 0, M, M, 0]$
r2 : [-2, 3, 4, -1, 0, 0, 1, 0, 12]$
r3 : [3, 2, 1, 0, -1, 0, 0, 1, 6]$
r4 : [1, 1, 1, 0, 0, 1, 0, 0, 9]$
tableau ( [r1, r2, r3, r4] );

The Simplex algorithm is used as usual. However when we compare two non-basic
variable coefficients which both involve M, the coefficient with a greater multiplier for M
is considered the larger coefficient.

After solving the big-M LP, we check if all artificial variables are equal to zero in the 
optimal solution. If they are all equal to zero, the corresponding solution for the 
decision variables is optimal for the given original LP.

If at least one of the artificial variables in the optimal solution is nonzero, the original LP
is not feasible.

Step 0 Tableau:

We first need to check there are no basic variable coefficients in the z-row (row 1). 
Since a1 and a2 are basic variables in the initial tableau we need to add to r1 the 
quantity  -M*r2 -M*r3

r1 - M*r2 - M*r3;

r1 : %;

r1 : r1 - M*r2 - M*r3$
RL : [r1, r2, r3, r4]$
tableau (RL);

Step 1 Tableau:

The most negative coefficient in the z-row (row 1) is -5*M - 3, so x3 is the pivot variable,
x3 enters the basis. There are three positive coefficients (4,1,1) in the x3 constraint rows,
with corresponding ratios (3, 6, 90) so row 2 is the winner of the minimum ratio test, and
row 2 is the pivot row, the pivot element is row 2, col 3 with value 4. a1 leaves the basis
(becomes equal to 0). Using pivot1 we get the Step 1 tableau.

tratio (RL, 3);

RL : pivot1 (RL, [2,3])$

In the z-row (r1) the x1 column has the most negative coefficient (-7/2)*M.

tratio (RL, 1);

Step 2 Tableau:

RL : pivot1 (RL, [3,1])$

tratio (RL, 5);

Step 3 Tableau:

RL : pivot1 (RL, [4,5])$

tratio (RL, 4);

Step 4 Tableau:

RL : pivot1 (RL, [3,4])$

There are no negative coefficients in the z-row, both a1 and a2 are non-basic variables 
and therefore equal zero, so we have the optimum solution to the original LP 
with x1 = x2 = 0,  x3 = 9, z = 27.

## An Unbounded Example with the Two Phase Method

We follow Butenko's video 21:
https://www.youtube.com/
       watch?v=_wnqe5_CLU0&list=PLY9yf2-4yyeT0Jh80w5yuevm0o9lAAUJg&index=21

The example is:
  Maximize 6*x1 + x2,
   subject to:
    - x1 + 3*x2 <= 6,
     x1  - 3*x2 = 6,
     x1 + x2 >= 1,  with
     x1, x2 >= 0.

maximize_lp (6*x1 + x2, [ - x1 + 3*x2 <= 6,  x1  - 3*x2 = 6, x1 + x2 >= 1], [x1, x2] );

### Standard Form

The standard form is
 Maximize  6*x1 + x2
subject to
    - x1 + 3*x2 + s1 = 6,
     x1  - 3*x2 = 6,
     x1 + x2 - e3 = 1,
with x1, x2, s1, e3 >= 0.

### Phase I: the Auxiliary Problem

If we assume an initial feasible solution in which x1 = x2  = 0, we see that we 
need to add an artificial variable a2 to the second constraint equation and add
an artificial variable a3 to the third constraint equation. We then  have
the auxiliary problem:
  Minimize (a2 + a3), 
subject to 
   - x1 + 3*x2 + s1 = 6,
     x1  - 3*x2 + a2 = 6,
     x1 + x2 - e3 + a3 = 1,
with x1, x2, s1, e3, a2, a3  >= 0.

Initializing Phase I:

We then Maximize z = (- a2 - a3), which we write as z + a2 + a3 = 0 to form a
coefficients tableau.

vL : [x1, x2, s1, e3, a2, a3, rhs, BV]$
bL : [z, s1, a2, a3]$
r1 : [0,0,0,0,1,1,0]$
r2 : [-1,3,1,0,0,0,6]$
r3 : [1,-3,0,0,1,0,6]$
r4 : [1,1,0,-1,0,1,1]$
tableau ([r1,r2,r3,r4]);

Step 0 Tableau:

Correct row 1 so that all elements in the a2 column (except the a2 row) are zero and
likewise for the a3 column.

r1 : r1 - r3 - r4$
RL : [r1, r2, r3, r4]$
tableau (RL);

This is the Step 0 tableau:

tratio (RL, 1);

Step 1 Tableau:

RL : pivot1 (RL, [4,1])$

Step 2 Tableau:

RL : pivot1 (RL, [3,4])$

There are no more negative coefficients in the z-row (row 1), a2 and a3 are non-basic
variables and therefore zero, z = 0 and the original LP is therefore feasible. We eliminate
the a2 and a3 columns from the tableau and replace the z-row (row 1) by a row 
corresponding to the given LP with z - 6*x1 - x2 = 0.

### Phase II

Remove column 5 and 6 elements from r2,r3,r4,vL, and redefine r1.

r1 : [-6,-1, 0, 0, 0]$
r2 : remL (RL[2], [5, 6])$
r3 : remL (RL[3], [5, 6])$
r4 : remL (RL[4], [5, 6])$
vL : remL (vL, [5, 6])$
tableau ([r1, r2, r3, r4] );

x1 is a basic variable, so correct the z-row (row 1) to a coefficient 0 for x1.

r1 : r1 + 6*r4;

We then have the initial feasible tableau for Phase II, which is the Step 0 tableau.

tableau ([r1,r2,r3,r4]);

The only negative coefficient in the z-row is -19 in the x2 column. There are no positive
coefficients in the column x2 constraint rows. This means we can increase x2 from
zero to positive infinity without any constraints. Therefore the problem is unbounded.

Let's be more explicit and unravel this tableau. 

The four rows say: z = 36 + 19*x2, s1 = 12, e3 = 5 + 4*x2, and x1 = 6 + 3*x2, with 
current values x1 = 6, x2 = 0, s1 = 12, and e3 = 5. We can now increase x2 from its 
current zero value to plus infinity, increasing z without bound. s1 will remain equal to 
12, and both e3 and x1 will become arbitrarily large and positive without bound. 

## An Unbounded Example with the Big M Method

This example is from Butenko's video 22:
https://www.youtube.com/
      watch?v=4qdWAETaFrU&list=PLY9yf2-4yyeT0Jh80w5yuevm0o9lAAUJg&index=22

The example is the same as in the previous section where we used the two phase
method to deal with artificial variables.
  Maximize 6*x1 + x2,
   subject to:
    - x1 + 3*x2 <= 6,
     x1  - 3*x2 = 6,
     x1 + x2 >= 1,  with
     x1, x2 >= 0.

The standard form is

 Maximize  6*x1 + x2
subject to
    - x1 + 3*x2 + s1 = 6,
     x1  - 3*x2 = 6,
     x1 + x2 - e3 = 1,
with x1, x2, s1, e3 >= 0.

We again need two artificial variables a2 and a3 in the second and third constraint 
equation. In the Big M method we subtract M*a2 and M*a3 from the objective function 
and seek a solution in which a2 and a3 become equal to zero.

 Maximize  6*x1 + x2 - M*a2 - M*a3
subject to
    - x1 + 3*x2 + s1 = 6,
     x1  - 3*x2 + a2 = 6,
     x1 + x2 - e3 + a3 = 1,
with x1, x2, s1, e3, a2, a3 >= 0.

Our z-row corresponds to the equation z - 6*x1 - x2 + M*a2 + M*a3 = 0.

vL : [x1, x2, s1, e3, a2, a3, rhs, BV]$
bL : [z, s1, a2, a3]$
r1 : [-6, -1, 0, 0, M, M, 0]$
r2 : [-1,3,1,0,0,0,6]$
r3 : [1,-3,0,0,1,0,6]$
r4 : [1,1,0,-1,0,1,1]$
tableau ([r1, r2, r3, r4]);

Step 0 tableau:

r1 : r1 - M*r3 - M*r4$
RL : [r1, r2, r3, r4]$
tableau (RL);

tratio (RL, 1);

RL : pivot1 (RL, [4,1])$

RL : pivot1 (RL, [3,4])$

The only negative coefficient in the z-row is -19, and there are no positive coefficients in
the x2 column constraint rows. This means we can increase x2 from zero to plus infinity
without meeting any constraints. The LP problem is unbounded.

## Alternative Optima and the Simplex Method

### The Heavenly Pouch LP

Maximize z  = 15*x1 + 25*x2
  subject to
  x1 + x2 <= 450,  --->  x1 + x2  + s1 = 450,
         x2 <= 300,  --->  x2 + s2 = 300,
 4*x1 + 5*x2 <= 2000,  --->  4*x1 + 5*x2 + s3 = 2000,
x1 <= 350,  --->  x1 + s4 = 350,
 with  x1, x2, s1, s2, s3, s4 >=0.

maximize_lp (15*x1 + 25*x2, [x1 + x2 <= 450, x2 <= 300,4*x1 + 5*x2 <= 2000,
        x1 <= 350], [x1, x2] );

maxima_lp finds the optimum z = 9375 with x1 = 125, and x2 = 300.

Butenko looks at the z-row of the optimum tableau to see if any non-basic variables have
a zero coefficient in the z-row (a sign that it is possible there is more than one optimum
solution).

Using the solution found by maximize_lp, we can solve for the final tableau values of the
four slack variables.
  s1 = 450 - x1 - x2 = 25,
  s2 = 300 - x2 = 0,
  s3 = 2000 - 4*x1 - 5*x2 = 0,
  s4 = 350 - x1 = 225.

Hence s1 and s4 are basic variables (> 0) and s2 and s3 are non-basic variables (= 0) 
in the optimum tableau.

Do either s2 and s3 have zero coefficients in the z-row of the optimum tableau? 

Step 0 tableau:

vL : [x1, x2, s1, s2, s3, s4, rhs, Basis]$
bL : [z, s1, s2, s3, s4]$
r1 : [-15,-25, 0, 0, 0, 0, 0]$
r2 : [1,1,1,0,0,0,450]$
r3 : [0,1,0,1,0,0, 300]$
r4 : [4,5,0,0,1,0,2000]$
r5 : [1,0, 0, 0, 0, 1, 350]$
RL : [r1, r2, r3, r4, r5]$
tableau (RL);

tratio (RL, 2);

RL : pivot1 (RL, [3, 2])$

tratio (RL, 1);

RL : pivot1 (RL, [4, 1])$

There are no more negative coefficients in the z-row (Row 1), so this is the optimum 
solution. Both of the non-basic variables (s2, s3) have nonzero coefficients in the 
z-row, a sign that this LP has a unique optimum solution. 

Since z = 9375 - (25/4)*s2 - (15/4)*s3, with both s2 and s3 equal to zero, at  the end of
this step, and since s2 and s3 cannot be negative (to satisfy the given conditions), 
there is no way to increase the value of z.

### A Modified Heavenly Pouch LP

If we change the objective to
  Max z = 20*x1 + 25*x2, with the same four constraints,

maximize_lp (20*x1 + 25*x2, [x1 + x2 <= 450, x2 <= 300,4*x1 + 5*x2 <= 2000,
        x1 <= 350], [x1, x2] );

vL : [x1, x2, s1, s2, s3, s4, rhs, Basis]$
bL : [z, s1, s2, s3, s4]$
r1 : [-20,-25, 0, 0, 0, 0, 0]$
r2 : [1,1,1,0,0,0,450]$
r3 : [0,1,0,1,0,0, 300]$
r4 : [4,5,0,0,1,0,2000]$
r5 : [1,0, 0, 0, 0, 1, 350]$
RL : [r1,r2,r3,r4,r5]$
tableau (RL);

tratio (RL, 2);

RL : pivot1 (RL, [3, 2])$

tratio (RL, 1);

RL : pivot1 (RL, [4,1])$

This optimum tableau has z = 10,000, x1 = 125, x2 = 300, whereas maxima_lp found the
optimum z = 10,000 for a different set: x1 = 250 and x2 = 200. Of the two non-basic 
variables, (s2 and s3) s2 has a zero coefficient in the z-row, indicating the possibility of 
alternate optimum solutions. An alternate solution was returned by maxima_lp, in which 
z = 10,000 but x1 = 250 and x2 = 200.

Because there are only two decision variables (x1 and x2), one can plot the constraint 
conditions with a 2d plot and visually identify cases in which more than one optimum 
case exists. See Butenko's video.

With s2 non-basic and equal to zero, we can try using s2 as the pivot variable, and 
as usual, use the minimum ratio test to select the pivot row: s2 enters the basis with 
a positive value and s1 leaves the basis with zero value, and our pivot element is 
row 2 and col 4:

tratio (RL, 4);

RL : pivot1 (RL, [2, 4])$

We have arrived at the solution returned by maxima_lp: x1 = 250, x2 = 200, z = 10,000.
Can we continue to play this game? The non-basic variable s1 has a zero coefficient in 
the z-row, so we let s1 be the pivot variable:

tratio (RL, 3);

RL : pivot1(RL, [2, 3])$

We have returned to our first simplex algorithm solution x1 = 125, x2 = 300, z = 10,000.

The presence of a zero coefficient of a non-basic variable in the z-row does not 
guarantee an alternative optimum solution. For example, if we seek to use the minimum 
ratio test for the positive coefficients in the pivot column constraint rows, and the rhs 
value is equal to zero, instead of finding another optimum solution you arrive at a case 
of degeneracy in which the next pivot results in the same solution that you start with. 

See Butenko, Video 17 for a  discussion with examples of degeneracy and cycling. One 
will change the basic variables in a degenerate iteration and it may be possible to get 
out of a degenerate trap. However the usual alternative (and what we have here) is
 "cycling." 

## The Matrix Form of the Simplex Method 

### Using Mdefine, Mdisplay, Mtableau, Mratio, Mpivot

We use a maximization example from Butenko's Video 42.

Maximize 4*x1 + 3*x2 + 5*x3,
subject to
  x1 + 2*x2 + 2*x3 <= 4,
  3*x1 + 4*x3 <= 6,
  2*x1 + x2 + 4*x3 <= 8,
   with  x1, x2, x3 >= 0.

We first use maximize_lp (objective, condL, posL).

maximize_lp (4*x1 + 3*x2 + 5*x3, [x1 + 2*x2 + 2*x3 <= 4, 3*x1 + 4*x3 <= 6,
  2*x1 + x2 + 4*x3 <= 8], [x1,x2,x3]);

Adding slack variables x4, x5, x6, this becomes:

Maximize 4*x1 + 3*x2 + 5*x3
subject to
  x1 + 2*x2 + 2*x3 + x4 = 4,
  3*x1 + 4*x3 + x5 = 6,
  2*x1 + x2 + 4*x3 + x6 = 8,
   with  x1, x2, x3, x4, x5, x6 >= 0.

### Define Initial Global Entities: X, c, A, b, NV, BV

The matrix method functions use the following global matrices:
       c, X, A, b,  N, B, Xn, Xb, cN, cB,
and the global lists NV, BV.

The first four (c, X, A, b) remain the same throughout the optimum solution process.

c, X, and b are matrix column vectors. A is a matrix of coefficients of the
left-hand sides of the constraint equations.  b is a matrix column vector of the
right-hand side constants from the constraint equations. X is a matrix column vector whose
elements are the variables,

The remaining eight global entities (NV, BV, N, B, Xn, Xb, cN, cB) are redefined in 
each simplex algorithm step (each pivot).
 
 We seek to optimize z = transpose (c) . X subject to:
         A . X = b, 
         with    X >= 0.

With N a matrix of the columns of A corresponding to the non-basic variables
 (in order left to right) and with B a matrix of the columns of A corresponding to the
 basic variables (in order left to right), the matrix equation A . X = b turns into the
 matrix equation  N . Xn + B . Xb = b, which can then be solved for Xb (the basic variables)
 in terms of the non-basic variables Xn:
       Xb = - invert(B) . N . Xn + invert(B) . b

NV is an ordered list of the numbers associated with the non-basic variables for a
particular tableau, and BV is an ordered list of the numbers associated with the basic 
variables of the same particular tableau.

cB is a matrix row vector with elements taken from the row vector transpose(c) which 
correspond to the variables in Xb. Likewise cN is a matrix row vector with elements taken 
from the row vector transpose(c) which correspond to the variables in Xn. Then we have 
the identity
     z = transpose(c) . X = cB . Xb + cN . Xn.
Using Xb = - invert(B) . N . Xn + invert(B) . b, we can express z entirely in terms of the 
variables in Xn to get:
      z = cB . invert(B) . b  + (cN - cB . invert(B) . N) . Xn.
The first term is some number (a constant) and the second term is a function
of the non-basic variables at a given step. Since all non-basic variables, at a given step,
are equal to zero, the second term is zero, and the current value of z is the first term.

The association of variables names with variable numbers (defined by our definition of
the matrix column vector X), is (for this problem):
    [x1,x2,x3,x4,x5,x6]  <==>  [1,2,3,4,5,6].

We convert
Maximize z = 4*x1 + 3*x2 + 5*x3
subject to
  x1 + 2*x2 + 2*x3 + x4 = 4,
  3*x1 + 4*x3 + x5 = 6,
  2*x1 + x2 + 4*x3 + x6 = 8,
   with  x1, x2, x3, x4, x5, x6 >= 0,

into matrix form:

X : cvec ([x1,x2,x3,x4,x5,x6]);
c : cvec ([4,3,5,0,0,0]);
A : matrix ( [1,2,2,1,0,0], [3,0,4,0,1,0], [2,1,4,0,0,1] );
b : cvec ([4,6,8]);
NV : [1,2,3];
BV : [4,5,6];

The above definitions of X, c, A, and b will remain unaltered throughout the optimization
process. The definitions of NV and BV will change, swapping a basic and non-basic 
variable to produce the next simplex step (pivot).

### Mdefine() to define global N, B, Xb, Xn, cB, cN

Mdefine() uses global X, c, A, b, NV, and BV to define global N, B, Xn, Xb, cN, and cB.
Unless details is set to true (default is false), the results are not displayed.

Here is a look at the Maxima code for Mdefine():

display2d:false$
fundef (Mdefine);
display2d:true$

Here we call Mdefine():

Mdefine()$

### Mdisplay() to see NV, BV, N, B, Xb, Xn, cB, cN values

Since Mdefine(), by default, does not display the current definitions of N, B, Xn, Xb,
cN, and cB, you can use Mdisplay()$ to have the current values of these global
variables displayed, and as a bonus, a reminder of NV and BV.

Mdisplay()$

cN and cB are each matrix row vectors.

### Mtableau() to see current Matrix tableau

Once the lists NV and BV have been defined, and the values of N, B, Xn, Xb, cN,
and cB have been defined, the function Mtableau() creates a "matrix tableau" of the 
current state of the optimization process.

Mtableau() creates a matrix display based on the following organization.
ident(m) produces an (m x m) identity matrix in Maxima. The code for Mtableau()
uses length (Xb) to define m. Note that use of the Maxima function matrix provides
a method of making a table in which Maxima formats the width of each column to
accommodate the width of the various column elements.

matrix ( ["transpose(Xb)", "transpose(Xn)", "  ", "  " ],
    ["Xb as row","Xn as row","rhs", "Basis"],
    [0, "- (cN - cB . invert(B) . N)", "cB . invert(B) . b","z"],
    ["ident(m)", "invert(B) . N", "invert(B) . b","Xb"]);

The first numerical row is called the z-row in the following.

All the basic variables, whose symbols are in transpose(Xb) (in the first row) have 
matrix columns on the far left, all together, in the variable order established with the 
definition of the matrix column vector X.

All the non-basic variables, whose symbols are in transpose(Xn) (in the first row), have
matrix columns all together in the middle, with the variable order established by the
 vector X.

In each step of the simplex algorithm, one non-basic variable and one basic variable
trade roles, and appear in the opposite set of columns, maintaining the established
order set by the constant vector X.

Here is a look at the Maxima code of Mtableau():

display2d:false$
fundef (Mtableau);
display2d:true$

fundef (zrow);

fundef (coefNV);

Here is the Step 0 tableau, the initial feasible solution tableau:

Mtableau();

The first numerical row, the z-row, is meant to be read as
   z - 4*x1 - 3*x2 - 5*x3 = 0.  (since x1 = x2 = x3 = 0)
The three constraint condition rows can be read off in a manner not too different
from our previous work with tableau.
The first constraint row says:
   x1 + 2*x2 + 2*x3 + x4 = 4  (with x4 = 4 since x1 = x2 = x3 = 0 ), and so on.

### Mratio(colNum) to see ratios for the minimum ratio test

colNum (in the function Mratio) is the column number of the matrix of coefficients 
invert(B) . N (the middle matrix in the lower structure), and is usually the same 
number as the column number of the row matrix zrow() which has the largest negative 
coefficient. The x3 variable (column 3 of the middle matrix) has the coefficient (-5), which
is the most negative coefficient in the z-row. The x3 variable will be the "pivot variable",
and "x3 enters the Basis" with a positive value.

Mratio(3);

The "minimum ratio test" is won by the variable x5, which is the fifth variable in X (variable
number 5).

### Mpivot (n_Enter, n_Leave) for one simplex step

In Mpivot(n_Enter, n_Leave), n_Enter is the number of the variable entering the
Basis and n_Leave is the number of the variable leaving the basis, based on the
choice of the user. The simplex algorithm says make this choice based on the
non-basic variable with the most negative coefficient in the z-row and the basic variable
for which the ratio test produces the smallest value. The "variable number" depends on
the order of the variable in the constant vector X.

Here is a look at the Maxima code for Mpivot:

display2d:false$
fundef(Mpivot);
display2d:true$

Step 1 tableau with the third variable in X entering and the fifth variable in X leaving:

Mpivot (3, 5);

The first simplex pivot operation has caused z to increase from 0 to 15/2 = 7.5. The 
variable x3 has increased from 0 to 3/2 = 1.5.

The most negative coefficient in the z-row is in column 2 of the middle structure.

Mratio(2);

Since the most negative coefficient in the z-row corresponds to the second variable
in the matrix vector X, and the  minimum ratio test is won by the row corresponding to
the fourth variable in X, we have n_Enter = 2, n_Leave = 4.

Step 2 tableau:

Mpivot (2, 4);

The second simplex pivot operation has caused z to increase from 7.5 to 9.

At this stage, there is only one negative coefficient in the z-row, corresponding to the
first column of the middle structure, and in that first column there is only one positive 
coefficient (3/4) corresponding to the variable x3, the third variable in X.

Mratio(1);

Step 3 tableau:

Mpivot (1,3);

In this matrix tableau, there are no negative coefficients in the z-row, z has reached the
value 11, and the variables have the values:
  x1 = 2, x2 = 1, x3 = 0, x4 = 0, x5 = 0, x6 = 3.
This is the optimum solution, in agreement with maximize_lp.

### The Heavenly Pouch Example

We repeat the Heavenly Pouch example using the matrix method. This example
shows the importance of keeping track of the variable number as defined by the
matrix column vector X, since instead of using slack variables x3, x4, x5, x6, we use
variable names s1, s2, s3, s4.

There are two input variables (decision variables) x1 and x2. There are four
slack variables, which are called s1, s2, s3, and s4. s1 is variable 3, s2 is variable 4,
s3 is variable 5, and s4 is variable 6 (as will be seen from our definition of the matrix
column vector X below).

Maximize z  = 15*x1 + 25*x2
  subject to
  x1 + x2 <= 450,  --->  x1 + x2  + s1 = 450,
         x2 <= 300,  --->  x2 + s2 = 300,
 4*x1 + 5*x2 <= 2000,  --->  4*x1 + 5*x2 + s3 = 2000,
x1 <= 350,  --->  x1 + s4 = 350,
x1, x2, s1, s2, s3, s4 >=0.

### Using maximize_lp (obj, condL, posL)

maximize_lp (15*x1 + 25*x2, [x1 + x2 <= 450, x2 <= 300,4*x1 + 5*x2 <= 2000,
        x1 <= 350], [x1, x2] );

### Matrix Method:  X, c, A, b, NV, BV

The correspondence between variable name and variable number, defined by our
definition of the matrix column vector X, is:
   [x1, x2, s1, s2, s3, s4]  <==> [1, 2, 3, 4, 5, 6].

The entities X, c, A, and b remain the same throughout the optimization process.
The entities NV, BV, N, B, Xn, Xb, cN, and cB change at each simplex step.
NV is an ordered list of the numbers associated with the non-basic variables for a
 particular tableau, and BV is an ordered list of the numbers associated with the basic 
variables of the same particular tableau.

The matrices are related to the specific problem via the matrix problem
statement:
maximize  z = transpose(c) . X, subject to A . X = b, with X >= 0.

For the initial feasible solution, x1 and x2 are equal to zero, so the list NV is [1,2].
the basic variables are s1,s2,s3,and s4, corresponding to the third, fourth, fifth and
sixth variables, so the list BV is [3,4,5,6], with s1 = 450, s2 = 300, s3 = 2000, s4 = 350.

X : cvec ([x1,x2,s1,s2,s3,s4])$
c : cvec ([15, 25, 0,0,0,0])$
b : cvec ([450, 300, 2000, 350])$
A : matrix ([1,1,1,0,0,0], [0,1,0,1,0,0], [4,5,0,0,1,0], [1,0,0,0,0,1] )$
NV : [1,2]$
BV : [3,4,5,6]$
z = transpose(c) . X;
A . X = b;
Mdefine()$
Mtableau();

The coefficient of x2 in the z-row is the most negative, so x2 enters the Basis. In this
matrix tableau the x2 column happens to be column 2 of the middle structure, so we
calculate the ratios using Mratio(col-num = 2) here.

Mratio(2);

The minimum ratio corresponds to variable s2 which is variable number 4. s2 = var. 4
leaves the Basis.

Step 1 tableau:  use Mpivot (n_Enter, n_Leave) = Mpivot (2, 4):

Mpivot (2, 4);

Mratio (1);

The minimum ratio appears in the s3 row, and s3 is variable number 5.

Step 2 tableau:

Mpivot (1, 5);

This is the optimum matrix tableau, z = 9375, x1 = 125, x2 = 300, no more negative
coefficients in the z-row. s1 = 25, s2 = 0, s3 = 0, s4 = 225.

### Minimization Problem using the Matrix Method

Here is the  minimization LP example we discussed above using the tableau,
pivot1 method:

Minimize  w = -8*x1 - 10*x2 - 7*x3
subject to
   x1 + 3*x2 + 2*x3 <= 10,
  -x1 - 5*x2 - x3 >= -8,
   x1, x2, x3 >= 0. 

minimize_lp (-8*x1 - 10*x2 - 7*x3, [x1 + 3*x2 + 2*x3 <= 10, -x1 - 5*x2 - x3 >= -8],
                                                                          [x1, x2, x3] );

"Hand calculation" using the simplex method involves several steps:

Step 1: Convert problem to "standard form", which means
  (1) Must be written as a maximization problem,
  (2) All constraint conditions must be either (for decision variables = x1,x2,x3 case)
           f(x1,x2,x3) <= positive-constant,  or
           f(x1,x2,x3) >= positive-constant.
  (3)  All decision variables must be non-negative.

To transform a minimization LP into a maximization LP, multiply both sides of
the objective function by (-1). Thus
   Minimize  w = -8*x1 - 10*x2 - 7*x3 is equivalent to
    Maximize (-1) * [ w = -8*x1 - 10*x2 - 7*x3],  or
letting z = -w,
    Maximize   z = 8*x1 + 10*x2 + 7*x3.

To transform g(x1,x2,x3 >= -A, where A > 0, into a type <= condition, multiply both 
sides of the given condition by (-1):
  (-1)* [ -x1 - 5*x2 - x3 >= -8 ] then becomes
     x1 + 5*x2 + x3 <= 8.

We now have the maximization linear problem (LP) in the "standard form":
  maximize z = 8*x1 + 10*x2 + 7*x3
  subject to:
    x1 + 3*x2 + 2*x3 <= 10,
   x1 + 5*x2 + x3 <= 8,
    with  x1,x2,x3 >= 0.
or, introducing slack variables x4 and x5,
subject to:
  x1 + 3*x2 + 2*x3 + x4 = 10,
   x1 + 5*x2 + x3 + x5 = 8,
    with  x1,x2,x3,x4,x5 >= 0.

Using now the matrix method of solution, we first define X, c, A, b, and starting values
of NV and BV.

X : cvec ([x1,x2,x3,x4,x5 ])$
c : cvec ([8,10,7,0,0 ])$
A : matrix ( [1,3,2,1,0], [1,5,1,0,1] )$
b : cvec ([10,8])$
NV : [1,2,3]$
BV : [4,5 ]$
z = transpose(c) . X;
A . X = b;
Mdefine()$
Mtableau();

Mratio(2), numer;

Mpivot (2, 5);

Mratio(1);

Mpivot (1, 2);

z* = 64 for x1* = 8, x2* = 0, x3* = 0.
So the minimum of w = -z is -64 with x1 = 8, x2 = 0, x3 = 0, in agreement with
minimize_lp.

### BN Prob. 3.1: tableau to Mtableau with toMTab(rows)

B/N solved problem 3.1, p. 34:

Maximize z = x1 + 9*x2 + x3
 subject to
  x1 + 2*x2 + 3*x3 <= 9 ==> x1 + 2*x2 + 3*x3 + x4 = 9,
  3*x1 + 2*x2 + 2*x3 <= 15  ==> 3*x1 + 2*x2 + 2*x3 + x5 = 15,
                                                  x1,x2,x3,x4,x5 >= 0.

B/N refers to the Schaum's Outlines text "Operations Research, 2nd ed" by Richard
Bronson and G. Naadimuthu.

maximize_lp (x1 + 9*x2 + x3, [x1 + 2*x2 + 3*x3 <= 9, 3*x1 + 2*x2 + 2*x3 <= 15],
      [x1, x2, x3]);

We will work this problem using the non-matrix tableau method and display the
corresponding matrix tableau, allowing comparisons at each step. This is a simple
problem with two slack variables, x4 and x5.

We illustrate the function toMTab(Rows) which shows the matrix tableau corresponding
to any simplex step set of tableau rows. Before calling toMTab(Rows) for the first time,
you should define the global matrix column vector c such that the original objective
can be written z = transpose(c) . X, where X will be created from global list vL by the
function toMTab, and X is the matrix column vector whose elements are the symbols of
the variables.

Then when called, toMTab(Rows) will define the matrices X, A, b, and the lists NV, BV,
before calling Mdisplay() and then Mtableau().

In the non-matrix method, the list r1 corresponds to the equation
z - x1 - 9*x2 - x3 + 0*x4 + 0*x5 = 0. The lists r2 and r3 correspond to the two constraint
equations. Each of these lists will have six elements, corresponding to five variables
and the sixth entry is the "rhs".

vL : [x1, x2, x3, x4, x5, rhs, Basis]$
bL : [z, x4, x5]$
r1 : [-1, -9, -1, 0, 0, 0]$
c : cvec (-r1)$
r2 : [ 1, 2, 3, 1, 0, 9]$
r3 : [ 3, 2, 2, 0, 1, 15]$
RL : [r1, r2, r3]$
tableau (RL);
toMTab(RL);

Continuing with the tableau method:

tratio (RL, 2);

The tableau pivot element is row 2, col 2, value = 2.
We define a new set of rows, RL, as the output of pivot1.

RL : pivot1(RL, [2,2])$

Since the z-row contains no negative coefficients, this is the optimum solution 
with z = 81/2, x2 = 9/2, x5 = 6, x1 = x3 = x4 = 0. This solution agrees with maximize_lp.

toMTab (RL);

### BN Prob 3.2: Minimization with Matrix Big M Method

We work this same problem (BN Prob. 3.2) later with the matrix two phase method.

minimize w = 80*x1 + 60*x2
  subject to
   0.2*x1 + 0.32*x2 <= 0.25,
   x1       +         x2  = 1, 
    x1,x2 >= 0.

minimize_lp likes whole number fractions better than decimals.

rat(0.32);

minimize_lp (80*x1 + 60*x2, [x1/5 + 8*x2/25 <= 1/4, x1 + x2 = 1], [x1, x2]);

float(%);

### Convert to a Maximization Problem in Standard Form

1. Convert to a maximization problem with z = - w:
  maximize z = - 80*x1 - 60*x2
  subject to
   0.2*x1 + 0.32*x2 <= 0.25,
   x1       +         x2  = 1, 
    x1,x2 >= 0.
2. Write in standard form. Use slack variable x3 in the first constraint condition.
   We need an artificial variable x4 in the second constraint in order to have a second
non-zero basic variable (x4 = 1) in the initial feasible solution (x1 = 0, x2 = 0, x3 = 1/4,
 x4 = 1).
   maximize z = - 80*x1 - 60*x2 - M*x4,
  subject to
   0.2*x1 + 0.32*x2 + x3 = 0.25,
   x1       +         x2  + x4 = 1, 
    x1,x2,x3,x4 >= 0.

To be able to compare with the output of minimize_lp, we will use:
0.2 = 2/10 = 1/5, 0.32 = 8/25.

rat(0.32);

The initial feasible solution is x1 = x2 =  0, x3 = 1/4, x4 = 1.

X : cvec ([x1,x2,x3,x4 ])$
c : cvec ([-80,-60,0,-M ])$
A : matrix ( [1/5,8/25,1,0 ], [1,1,0,1] )$
b : cvec ([1/4,1])$
NV : [1,2 ]$
BV : [3, 4 ]$
z = transpose(c) . X;
A . X = b;
Mdefine()$
Mtableau();

The coefficient of x2, (60 - M), is more negative than the coefficient of x1, (80 - M), in 
the z-row - remember that we treat M as an arbitrarily large positive number.
We then use Mratio (colNum) (here the second column of the Xn matrix) to perform the
ratio test for the second column (which happens to be the x2 column):

Mratio(2);

We then call Mpivot (n_Enter, n_Leave) with n_Enter = 2 (variable x2, the second
variable in X, and with n_Leave = 3 (variable x3, the third variable in X):

Mpivot(2, 3);

Mratio(1);

Mpivot (1, 4);

The large positive parameter M has disappeared the current value of z and from the
values of the decision variables x1 and x3. The values of both the slack variable x3 
and the artificial variable x4 are zero. This is the optimum solution.

Since z* = -215/3, the minimum value of w = -z is w* = 215/3 with x1* = 7/12, x2* = 5/12
in agreement with minimize_lp.

### BN Prob. 3.2: Two Phase Matrix Simplex Method

Repeating the statement of BN Prob. 3.2:

minimize w = 80*x1 + 60*x2
  subject to
   0.2*x1 + 0.32*x2 <= 0.25,
   x1       +         x2  = 1, 
    x1,x2 >= 0.

rat(0.32);

minimize_lp (80*x1 + 60*x2, [x1/5 + 8*x2/25 <= 1/4, x1 + x2 = 1], [x1, x2]);

float(%);

1. Convert to a maximization problem
  maximize z = - w = - 80*x1 - 60*x2, 
  subject to
   0.2*x1 + 0.32*x2 <= 0.25,
   x1       +         x2  = 1, 
    x1,x2 >= 0.
2. Write in standard form:
   We need an artificial variable x4 in the second constraint in order to have a second
non-zero basic variable (x4 = 1) in the initial feasible solution (x1 = 0, x2 = 0, x3 = 1/4,
 x4 = 1) - ie., as many basic variables as there are conditions.

   maximize z = - 80*x1 - 60*x2, 
  subject to
   0.2*x1 + 0.32*x2 + x3 = 0.25,
   x1       +         x2  + x4 = 1, 
    x1,x2,x3,x4 >= 0.

With x4 an artificial variable, Phase 1 seeks to minimize (drive to zero) the sum of
all artificial variables in the condition equations. Here we have only one artificial variable
x4, so we seek to minimize x4 (drive to zero). We then convert this Phase I objective
into: we seek to maximize ( - x4):

Auxiliary Problem:
maximize z = -x4, 
subject to
  x1/5 + 8*x2/25 + x3 = 1/4,
  x1   +  x2  + x4 = 1,
  with    x1, x2, x3, x4 >= 0.

### Phase I Matrix Method, Mdefine, Mtableau, Mratio, Mpivot

X : cvec ([x1,x2,x3,x4 ])$
c : cvec ([ 0, 0, 0, -1 ])$
A : matrix ( [1/5,8/25,1,0 ], [1,1,0,1] )$
b : cvec ([1/4,1])$
NV : [1,2 ]$
BV : [3, 4 ]$
Mdefine()$
z = transpose(c) . X;
A . X = b;

Step 0 matrix tableau:

Mtableau();

Choose column 1 as the pivot column:

Mratio(1);

Call Mpivot(Enter, Leave) with n_Enter = 1 (var x1 enters Basis) and with n_Leave = 4
 (var x4 leaves the Basis)  (the order of the variables is defined by the order in the
vector X): 

Mpivot(1, 4);

There are no more negative coefficients in the z-row, z = 0, and x4 is not in the Basis
and hence equals zero, so we have established that the maximization problem is
a feasible problem to solve, and we go on to Phase II.

### Phase II Matrix Method

From the final matrix tableau of Phase I we eliminate the artificial variable (x4) column,
work now with the three remaining variables x1, x2, x3, and rewrite the z-row in the 
original maximization form
  max z = - 80*x1 - 60*x2 + 0*x3,
with the condition equations from reading off the condition rows of the above matrix
tableau (ignoring the x4 column):
    x1 + x2 = 1,
    x3 + (3/25)*x2 = 1/20
to define the matrices A and b.

X : cvec([x1,x2,x3])$
c : cvec ([-80, -60, 0])$
A : matrix ( [1, 1, 0], [0, 3/25, 1])$
b : cvec ([1, 1/20])$
NV : [2]$
BV : [1, 3]$
z = transpose(c) . X;
A . X = b;
Mdefine()$
Mtableau();

Mratio (1);

Mpivot (2, 3);

This matrix tableau shows no more negative coefficients in the z-row, and we have
the optimum solution with z = - 215/3, with x1 = 7/12, and x2 = 5/12. Hence the minimum
of w (= - z) is 215/3 with the same final values of x1 and x2.