# LPduality.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

LPduality.wxmx assumes the reader has a general understanding of the simplex method
of solving for optimal solutions of "linear programming" problems (LP problems). An
introduction to simplex methods using Maxima can be found in our LPsimplex.wxmx.

1. We begin with a description of the tableau conventions in which the "z-row(s)" are
on the bottom of the tableau (rather than the top, as we used in LPsimplex.wxmx).

2. We next discuss the systems of symmetric dual LP pairs, including the Duality Theorem,
and the Complimentary Slackness Principle. Examples with and without artificial variables
are solved, using Maxima. For LP's with Step 0 tableau forms which include artificial
variables, we use B/N's "two phase simplex method", which splits the z-row into two
rows. (See References) We show how the solution to the other dual pair LP is embedded
in the final tableau values of the slack and surplus variables.

3. We next discuss systems of antisymmetric dual LP pairs. We show an explicit Maxima
matrix calculation which allows a prediction of a solution of the other antisymmetric 
dual pair LP, given the Basis variables in the optimal tableau of one of the pair.

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, 2022

## References

[B/N] Operations Research, 2nd ed, Schaum's Outlines, Richard Bronson and
  G. Naadimuthu

Sergiy Butenko, ISEN 620, A Survey of Optimization, 82 videos
https://www.youtube.com/playlist?list=PLY9yf2-4yyeQTLkCVFnGedyERuCjKW7kI

## Tableau "z-row" on the Bottom Convention

The reference (B/N) Operations Research, 2nd ed, by Richard Bronson and 
G Naadimuthu, uses tableau conventions in which the z-row is the bottom row,
instead of the top row (as was done in LPsimplex.wxmx). With minor changes, including
the use of the new function bratio(RL,ncol) [b for bottom z-row] instead of tratio(RL,ncol)
[t for top z-row], we can easily use the Bronson and Naadimuthu notation, which is
also used in many other references. In particular, we can still use tableau(RL) and
pivot1(RL, [nrow, ncol]). When using the two phase simplex method, in which the
z-row is split into two rows on the bottom, we can use the function b2ratio(RL, ncol).

### Step 0 Tableau Conventions

Let Xs be the matrix column vector of variables symbols, including slack, surplus and
artificial variables.
Xo = matrix column vector of the Step 0 Basis variable symbols, in the order these
variables appear in the constraints. 
Co = matrix column vector of the coefficients of the step 0 basis variables in the
objective, in the same order as the variable symbols in Xo.

Given the Step 0 LP:  maximize z = Cs^t . Xs, such that  As . Xs = E, with Xs >= 0,
the step 0 maximization tableau is, using these matrices:
                |           Xs^t                 |    rhs       |  Basis   
              -----------------------------------------------------------
                |            As                   |   E          |  Xo      
              -----------------------------------------------------------
                |   Co^t . As - Cs^t        | Co^t . E  |   z

Given the Step 0 LP:  minimize w = Cs^t . Xs, such that As . Xs = E, with Xs >= 0,
the step 0 minimization tableau is, using these matrices:
                |           Xs^t                 |      rhs        |  Basis          
              -----------------------------------------------------------
                |            As                   |      E           |   Xo
              -----------------------------------------------------------
                |   Cs^t - Co^t . As          | - Co^t . E  |   z

### B/N Prob. 3.1, Maximization, no Artificial Variables

maximize z = x1 + 9*x2 + x3,
  such that 
   x1 + 2*x2 + 3*x3 <= 9, 
 3*x1 + 2*x2 + 2*x3 <= 15, 
x1,x2,x3  >= 0.

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

We add nonnegative slack variables x4 and x5 to the right-hand sides of the condition
inequalities, converting them into condition equalities, and arrive at the Step 0 form:

maximize z = x1 + 9*x2 + x3 + 0*x4 + 0*x5,
such that
      x1 +  2*x2 + 3*x3 + x4       = 9,
   3*x1 + 2*x2 +  2*x3 +      x5 = 15,
with x1,x2,x3,x4,x5 >= 0.

An initial feasible solution has x1 = x2 = x3 = 0, x4 = 9, x5 = 15. 

Let G^t mean the transpose of the matrix G in the following.

If the write the Step 0 LP using matrices,
    maximize z = C^t . X such that A . X = B, with X >=0,

for this problem we then have:

X^t = (x1, x2, x3, x4, x5), C^t = (1, 2, 3, 0, 0), A = matrix ( [1,2,3,1,0], [3,2,2,0,1]), 
B^t = (9, 15), Xo^t = (x4, x5) refers to the Basis variables in the step 0 tableau,
and Co^t = (0, 0) refers to the objective coefficients of the step 0 Basis variables Xo.
The quantity Co^t . B = 0 = initial value of z, and finally, 
Co^t . A = (0, 0, 0, 0, 0) here.

The step 0 maximization tableau is, using these matrices, 
                |           X^t                |    rhs       |  Basis   
              -----------------------------------------------------------
                |            A                  |   B          |  Xo      
              -----------------------------------------------------------
                |   Co^t . A - C^t         | Co^t . B  |   z

in which Xo is the Step 0 Basis vector of variable symbols, in the order these basis 
variables appear in the step 0 constraint equations, with the same number of variables
as the vector B. C is the "cost vector", and Co is the cost vector associated with the Step 0
basis vector Xo variables, and contains the coefficients of the Xo variables in the objective,
 in the same order as the variables in Xo.

The last row is the "z-row", with the middle section containing -C^t as the main piece. 
The first term,  Co^t . A, is present just in case one or more of the basis vector Xo 
variables have nonzero objective coefficients. We want to start the simplex algorithm with 
all of the basis elements having zero z-row coefficients.

The list bL is a list of the Step 0 Basis variables, in the order they appear in the condition
equations, with z added as the last element. In this example, r1 and r2 begin with the 
coefficients coming from A . X = B, and end with the elements of B.
In this example, the last row of the tableau (the z-row) is defined by the list r3.
The list vL begins with the variable symbols in the order they appear in X, and end with
the words : rhs, Basis.

Step 0 tableau:

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

The most negative coefficient in the z-row is -9, so the pivot column is the x2 column,
ncol = 2.

Because the z-row is now the last row, we cannot use tratio(RL, ncol). Instead, use
the function bratio(RL, ncol) (b for bottom row is z-row) defined in Econ1.mac.

bratio(RL, 2);

The minimum ratio test is won by the x4 row (nrow = 1), so the pivot element has value
2,  (intersection of  row 1 and col 2). We can then use pivot1(RL, [nrow, ncol]) as usual
for the first simplex step to get the Step 1 tableau:

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

The z-row coefficients in the bottom row are now all positive, an optimum has been 
reached, with z* = 81/2, x1* = x3* = 0, x2* = 9/2, which is the same solution found by
maximize_lp.

### Minimization Example, no Artificial Variables

This minimization problem was considered in LPsimplex.wxmx in Sec. 5.2. The two
slack variables were called s1 and s2 there. C^t stands for the transpose of C.

minimize w = -8*x1 - 10*x2 - 7*x3
  s.t.  x1 + 3*x2 + 2*x3 <= 10,
         x1 + 5*x2 + x3 <= 8, 
   with 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]);

Converting from inequality conditions to equality conditions with slack variables x4, x5,

minimize w = C^t . X, s.t.  A . X = B, with X>=0,
   which translates to:  
minimize w = -8*x1 - 10*x2 - 7*x3 + 0*x4 + 0*x5,
s.t.  x1 + 3*x2 + 2*x3 + x4 = 10,
         x1 + 5*x2 + x3 + x5 = 8, 
   with x1, x2, x3, x4, x5 >= 0.

Here we have X^t = (x1, x2, x3, x4, x5), C^t = (-8, -10, -7, 0, 0), Xo^t = (x4,x5),
Co^t = (0, 0), B^t = (10, 8), A = matrix([1,3,2,1,0], [1,5,1,0,1]).

the step 0 minimization tableau is, using these matrices:

                |           X^t                 |    rhs        |  Basis          
              -----------------------------------------------------------
                |            A                   |    B           |   Xo
              -----------------------------------------------------------
                |   C^t - Co^t . A          | - Co^t . B  |   z

in which Xo is the known basic feasible solution Basis vector, with the same number
of elements as the rhs vector B, C is defined by min w =  C^t . X, and Co is contains the 
objective coefficients of the Xo variables in the same order as in Xo. The last row is 
the "z-row", with  the middle section containing C^t as the main piece; the second 
term (- Co^t . A) is present just in case one or more of the step 0 basis vector Xo 
elements have nonzero objective coefficients. We want to start the simplex algorithm 
with all of the basis elements having zero z-row coefficients.
 
The fact that we use C^t as the main part of the last row means that we are really going
to maximize z = -w with the simplex algorithm.

The list r1 elements begine with the first row of A, the last element is the first element
of B. Likewise r2 comes from A[2] and B[2,1]. r3 here is the bottom z-row, including the far
right "current value of z": - Co^t . B.

C : cvec([-8, -10, -7, 0, 0])$
Co : cvec([0, 0])$
B : cvec ([10,8])$
A : matrix([1,3,2,1,0], [1,5,1,0,1])$
"C^t - Co^t . A" = transpose(C) - transpose(Co) . A;
"- Co^t . B" = - transpose(Co) . B;

Step 0 tableau:

vL : [x1,x2,x3,x4,x5,rhs,Basis]$
bL : [x4,x5,z]$
r1 : endcons (B[1,1], A[1])$
r2 : endcons (B[2,1], A[2])$
r3 : [-8,-10,-7,0,0,0]$
RL : [r1,r2,r3]$
tableau(RL);

Choose the x2 column as the pivot column (col 2 here), use the minimum ratio test to
pick the pivot row (and hence the "pivot element"). We can either do this "by hand",

[10/3, 8/5], numer;

or use our function bratio(RL, ncol) [b for bottom z-row] to see the results of the 
minimum ratio test for column 2.

bratio(RL, 2);

Step 1 tableau:

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

z has increased from 0 to 16 going from Step 0 to Step 1.

bratio (RL, 1);

Step 2 tableau:

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

With no negative coefficients in the last row, the optimum solution maximizing z = -w is
z* = 64 with x1* = 8, x2* = x3* = 0. Hence min w = w* = -64, with x1* = 8, x2* = x3* = 0.
Recall that our method for minimizing w is to maximize z = -w.

## Symmetric Duals

We follow Ch. 4 in B/N.

Every linear program in the variables x1, x2, ..., xn has associated with it another linear 
program in the variables y1, y2, ..., ym (where m is the number of constraints in the 
original program), known as the "dual". The original program, called the "primal",
completely determines the form of its dual. 

For any matrix G, let G^t denote the transpose of G.

The dual of a (primal) linear program (LP) in the (nonstandard) matrix form is defined by:

Given: the primal LP in which X^t = (x1, x2, x3, ..., xn),
     maximize z = C^t . X  such that   A . X <= B, with X >= 0,                       (1)
 
then the corresponding dual LP in, which Y^t = (y1,y2,y3,...,ym), is:
    minimize w = B^t . Y such that A^t . Y >= C, with Y >= 0.                        (2)     

Programs (1) and (2) are "symmetrical" in that both involve nonnegative variables and
inequality constraints; they are known as the "symmetric duals" of each other. The
dual variables y1, y2,..., ym are sometimes called "shadow costs".

Butenko's Video 27, LP Duality: motivation (Operations Research Video Series) presents 
a useful approach to thinking about the relation of primal to dual.

### Duality Theorem

If an optimal solution exists to either the primal or symmetric dual program, then the
other program also has an optimal solution and the two objective functions have the
same optimal value.

In such situations, the optimal solution to the primal (dual) is found in the last row of 
the final simplex tableau for the dual  (primal), in those columns associated with the
slack or suplus variables (see Ex. 3). Since the solutions to both programs are obtained by
solving either one, it may be computationally advantageous to solve a program's dual
rather than the given program itself. 

### Example 1

Determine the symmetric dual of the program:

  maximize z = 20*x1 + 30*x2 + 40*x3 + 50*x4,
   subject to: 2*x1 + 6*x2 +   7*x3  +     x4 <= 5,
                    3*x1 + 8*x2  +     x3   + 2*x4 <= 2,
                      x1  + 5*x2  +  3*x3  +  4*x4 <= 1,
   with x1, x2, x3, x4 >= 0.

The dual LP is
  minimize w = 5*y1 + 2*y2 + y3,
  subject to: 2*y1 + 3*y2 +    y3 >= 20,
                   6*y1 + 8*y2 + 5*y3 >= 30,
                   7*y1 +    y2 + 3*y3 >= 40,
                      y1 + 2*y2 + 4*y3 >= 50,
   with y1, y2, y3 >= 0.

The primal LP has three inequality constraint conditions, so the dual has three variables.
The coefficients of the dual objective variables (y1,y2,y3) are the three right hand side
constants of the three primal constraints. The primal LP has four variables (x1,x2,x3,x4)
and hence the dual has four constraint conditions, with the right hand side constants
(20, 30, 40, 50) equal to the primal objective coefficients.

Using matrices to create the symmetric dual:

print ("primal LP")$
X : cvec([x1,x2,x3,x4])$
C : cvec ([20, 30, 40, 50])$
print ("maximize")$
z = transpose (C) . X;
B : cvec ([5,2,1])$
A : matrix ([2,6,7,1], [3,8,1,2], [1,5,3,4])$
print ("such that")$
A . X <= B;
print ("with X >= 0")$
print ("corresponding dual")$
Y : cvec ([y1, y2, y3])$
print ("minimize")$
w = transpose (B) . Y;
print ("such that")$
transpose (A) . Y >= C;
print ("with Y >= 0")$

### Example 2

Determine the symmetric dual of the program:
  maximize z = 2*x1 + x2,
   subject to: x1 + 5*x2  <= 10,
                    x1 + 3*x2  <= 6,
                   2*x1  + 2*x2  <= 8,
   with x1, x2 >= 0.

The dual LP is
  minimize w = 10*y1 + 6*y2 + 8*y3,
  subject to:    y1 + y2 +  2*y3 >= 2,
                   5*y1 + 3*y2 + 2*y3 >= 1,
   with y1, y2, y3 >= 0.

print ("primal LP")$
X : cvec([x1,x2])$
C : cvec ([2, 1])$
print ("maximize")$
z = transpose (C) . X;
B : cvec ([10,6,8])$
A : matrix ([1,5], [1,3], [2,2])$
print ("such that")$
A . X <= B;
print ("with X >= 0")$
print ("corresponding dual")$
Y : cvec ([y1, y2, y3])$
print ("minimize")$
w = transpose (B) . Y;
print ("such that")$
transpose (A) . Y >= C;
print ("with Y >= 0")$

### Example 3, Solutions to Primal and Dual LP's of Ex. 2

Show that both the primal and dual programs in Example 2 have the same optimal value.

### Primal LP Solution

The primal LP is:
maximize z = 2*x1 + x2,
   subject to: x1 + 5*x2  <= 10,
                    x1 + 3*x2  <= 6,
                   2*x1  + 2*x2  <= 8,
   with x1, x2 >= 0.

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

Adding slack variables x3, x4, x5 to the primal LP constraints,
maximize z = 2*x1 + x2 + 0*x3 + 0*x4 + 0*x5,
   subject to: x1 + 5*x2 + x3  = 10,
                    x1 + 3*x2 + x4 = 6,
                   2*x1  + 2*x2 + x5  = 8,
   with x1, x2, x3, x4, x5 >= 0.

The initial feasible solution (ifs) is: x1 = x2 = 0, x3 = 10, x4 = 6, x5 = 8, z = 0.

Using matrix notation with C^t standing for the transpose of C,
  maximize z = C^t . X, s.t. A . X = B, with X >=0, 

The step 0 maximization tableau is, using these matrices, 

                |           X^t                |    rhs       |  Basis   
              -----------------------------------------------------------
                |            A                  |   B          |  Xo      
              -----------------------------------------------------------
                |   Co^t . A - C^t         | Co^t . B  |   z

C^t = (2,1,0,0,0), X^t = (x1,x2,x3,x4,x5),Xo^t = (x3,x4,x5), Co^t = (0,0,0),B^t = (10,6,8).
Xo = matrix column vector of the initial Basis variable symbols, in the order these
variables appear in the condition equations.
Co = matrix column vector of the objective coefficients of the Xo variables, in the same 
order as the symbols in Xo.

An added step, if needed, is to calculate the two parts of the last row, using matrices.
Here it is obvious that Co^t . A = (0, 0, 0, 0, 0), Co^t . B = 0; the middle part of
the z-row is -C^t  = (-2,-1,0,0,0).

C : cvec([2,1,0,0,0])$
Co : cvec ([0,0,0])$
B : cvec ([10,6,8])$
A : matrix ([1,5,1,0,0], [1,3,0,1,0], [2,2,0,0,1])$
"Co^t . A - C^t" = transpose(Co) . A - transpose(C);
"Co^t . B" = transpose(Co) . B;

The list r1 comes from A[1] and B[1,1], likewise for list r2 and list r3, so each list
contains 6 elements. Since we have gone through the trouble of defining A and B,
we can show how this is done using endcons.

endcons ( B[1,1],  A[1]);

The list r4 contains the elements of the last row, including the current z value
Co^t . B as the last element.

The Step 0 tableau is (z-row on the bottom) [close to B/N's tableau but with Xo on the
right-hand side instead of the left-hand side]:

vL : [x1,x2,x3,x4,x5,rhs,Basis]$
bL : [x3, x4, x5, z]$
r1 : endcons (B[1,1], A[1])$
r2 : endcons (B[2,1], A[2])$
r3 : endcons (B[3,1], A[3])$
r4 : [-2, -1, 0, 0, 0, 0]$
RL : [r1,r2,r3,r4]$
tableau(RL);

We choose the x1 column (col 1 here) as the pivot column since the most negative
z-row coefficient (-2) is in col. 1.  The ratios rhs/x1-coeff here are obvious, but we
need practice with bratio(RL, ncol):

bratio (RL, 1);

Step 1 tableau using pivot1(RL, [nrow, ncol]):

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

With no negative coefficients remaining in the z-row (the last row), z* = 8 with x1* = 4
and x2* = 0, in agreement with maximize_lp.

Recall that x3, x4, x5 primal variables are all slack.
In the final primal tableau,
**the first primal slack variable x3 has z-row coefficient = 0 which implies  the 
   symmetric dual solution final tableau will have y1 = 0,
**the second primal slack variable x4 has z-row coefficient = 0 , which implies the  
   symmetric dual solution final tableau will have  y2 = 0,
**the third primal slack variable x5 has z-row coefficient = 1, which implies the 
    symmetric dual solution final tableau will have y3 = 1.

Here is a matrix method of finding the solutions to the symmetric dual LP from the primal
LP optimal tableau:

The basis variables in the optimal primal tableau are the first, third, and fourth variables 
of X,  (x1,x3,x4).
Let Cb be the matrix column vector containing the first, third, and fourth coefficients
in the primal vector C, and let Ab be the matrix containing the first, third, and fourth 
columns of the primal matrix A.
Then Y*^t = Cb^t . Ab^(-1) = transpose (Cb) . invert (Ab) = (y1, y2, y3).

Cb : part (C, [1,3,4]);

Ab : newM (A, [1,3,4]);

transpose (Cb) . invert (Ab);

which are the values of the dual LP solution  Y*^t = (y1, y2, y3).

### Dual LP Solution with B/N's Two Phase Simplex Method

The dual LP is:
minimize w = 10*y1 + 6*y2 + 8*y3,
  subject to:    y1 + y2 +  2*y3 >= 2,
                   5*y1 + 3*y2 + 2*y3 >= 1,
   with y1, y2, y3 >= 0.

minimize_lp (10*y1 + 6*y2 + 8*y3, [y1 + y2 +  2*y3 >= 2, 5*y1 + 3*y2 + 2*y3 >= 1],
       [y1, y2, y3]);

So minimize_lp has found the same solution for y1,y2,y3 as found in the primal
slack variable z-row coefficients in the optimal tableau.

To get constraint equations, subtract surplus (excess) variables y4 and y5 to get
  minimize w = 10*y1 + 6*y2 + 8*y3 + 0*y4 + 0*y5,
  subject to:    y1 + y2 +  2*y3 - y4 = 2,
                   5*y1 + 3*y2 + 2*y3 - y5 = 1,
   with y1, y2, y3, y4, y5 >= 0.

We cannot start with an initial feasible solution (ifs) y1 =  y2 = y3 = 0, since that would
require y4 = -2, y5 = -1, in contradiction with our nonnegativity assumptions for the
variables. So add artificial variables y6 and y7 to the left hand sides of the constraint
equations, and add M*y6 + M*y7 to w to get:
minimize  w = 10*y1 + 6*y2 + 8*y3 + 0*y4 + 0*y5 + M*y6 + M*y7,
  subject to:    y1 + y2 +  2*y3 - y4 + y6 = 2,
                   5*y1 + 3*y2 + 2*y3 - y5 + y7 = 1,
   with y1, y2, y3, y4, y5, y6, y7 >= 0.
Then our initial feasible solution is: y1 =  y2 = y3 = y4 = y5 = 0, y6 = 2, y7 = 1.
We seek an optimum solution in which the values of the artificial variables y6 and y7
have both been driven to zero.

Given the LP:
minimize w = C^t . X, s.t. A . X = B, with X >= 0,

the step 0 minimization tableau is, using these matrices:

                |           X^t                 |    rhs        |  Basis          
              -----------------------------------------------------------
                |            A                   |    B           |   Xo
              -----------------------------------------------------------
                |   C^t - Co^t . A          | - Co^t . B  |   z

Here, X^t = (y1,y2,y3,y4,y5,y6,y7), C^t = (10,6,8,0,0,M,M), B^t = (2,1),
Xo = matrix column vector of the initial Basis variable symbols, in the order these
variables appear in the condition equations; Xo^t = (y6,y7),
Co = matrix column vector of the objective coefficients of the initial Basis variables,
in the same order as the variable symbols in Xo; Co^t = (M,M).

The construction of this minimization tableau implies that  the simplex algorithm 
is being used to find the maximum of z = -w.

C : cvec ([10,6,8,0,0,M,M])$
Co : cvec ([M,M])$
B : cvec ([2, 1])$
A : matrix ([1,1,2,-1,0, 1, 0],[5,3,2,0,-1, 0, 1])$
"C^t - Co^t . A" = transpose(C) - transpose(Co) . A;
"- Co^t . B" = - transpose(Co) . B;

Step 0 tableau; we split the z-row into two rows.
If the original z-row has the form
  [a1 + b1*M, a2 + b2*M,...]
we set r3 to the list [a1, a2, ...] and r4 to the list [b1, b2, ...].

As long as we have already gone to the trouble of defining A and B, we can define
r1 and r2 using endcons.

Step 0 tableau, close to B/N's tableau, except Xo basis variables are on the rhs.

vL : [y1, y2, y3, y4, y5, y6, y7, rhs, Basis]$
bL : [y6, y7, z1, z2]$
r1 : endcons (B[1,1], A[1])$
r2 : endcons (B[2,1], A[2])$
r3 : [10, 6, 8, 0, 0, 0, 0, 0]$
r4 : [-6, -4, -4, 1, 1, 0, 0, -3]$
RL : [r1, r2, r3, r4]$
tableau(RL);

Looking for the most negative coefficient in the bottom row, we choose the y1 column as
the pivot column and use the minimum ratio test to choose the pivot row.

[2/1, 1/5], numer;

Since we have two z-rows, we need to use b2ratio (RL, ncol).

b2ratio (RL, 1);

So row 2 is the pivot row.

Step 1 tableau:

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

The most negative coefficient in the bottom row (excluding the "rhs" number) is -8/5
in column 3.

b2ratio (RL, 3);

The minimum ratio test is won by row 2

Step 2 tableau:

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

The only negative coefficient in the last row (excluding the "rhs" number) is in col 5
and the only positive coefficient in col 5 (excluding the last two rows) is in row 1.

Step 3 tableau:

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

Artificial variables y6 and y7 have left the Basis and are hence equal to zero. We can
delete columns 6 and 7, and then, since the bottom row will be all zeros, we can
delete the last column.

r1 : remL( RL[1], [6,7]);

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

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

vL : [y1,y2,y3,y4,y5,rhs,Basis]$
bL : [y5, y3, z]$
RL : [r1, r2, r3]$
tableau (RL);

Since there are no negative coefficients in the last row (excluding the "rhs" number),
the maximum of z = -w is z* = -8 with y1* = y2* = 0, y3* = 1, which implies that the
minimum of w = 10*y1 + 6*y2 + 8*y3 is w* = -z* = 8, with y1* = y2* = 0, y3* = 1,
in agreement with minimize_lp.

Recall that y4 and y5 are surplus (excess) variables.

In the final symmetric dual tableau,
  ** the first dual surplus variable y4 has z-row coefficient equal to 4, implies primal x1 = 4.
  **the second dual surplus variable y5 has z-row coefficient equal to 0, 
           implies primal x2 = 0.

Here is a matrix method of finding the solutions to the symmetric primal LP from the 
dual LP optimal tableau:
The basis variables in the dual optimal tableau are the third and fifth variables of Y
(y3, y5).
Let Cb be the matrix column vector conaining the third and fifth coefficients
in the dual LP vector C, and let Ab be the matrix containing the third and fifth columns 
of the dual LP matrix A.
Then X*^t = Cb^t . Ab^(-1) = transpose (Cb) . invert (Ab) gives the values of (x1, x2) in
the primal LP optimal tableau.

Cb : part (C, [3,5]);

Ab : newM (A, [3, 5]);

transpose (Cb) . invert (Ab);

If our goal is to find the optimal solution to the problem:

minimize w = 10*y1 + 6*y2 + 8*y3,
  subject to:    y1 + y2 +  2*y3 >= 2,
                   5*y1 + 3*y2 + 2*y3 >= 1,
   with y1, y2, y3 >= 0,

then, because of the need to introduce artificial variables, it is computationally 
advantageous to instead find the optimal solution to the dual problem:

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

Note that the dual of the dual is the primal. (See, for example, Butenko's video 28.)

### Complementary Slackness Principle

Given that the pair of symmetric duals have optimal solutions, then if the k'th
constraint of one system holds as an inequality - i.e., the associated slack or surplus
variable is positive for the optimal solution - the k'th component of the optimal solution of
its symmetric dual is zero.

Applying this principle to Example 3, the first and second primal slack variables x3, x4
are both greater than zero, so the first and second symmetric dual variables y1 and y2 
are equal to zero.

The second symmetric dual surplus variable y5 is greater than zero, so the second 
primal variable x2 is equal to zero.

### Example 4

This example is taken from B/N Prob. 4.4.
We take the primal LP to be:

minimize w = x1 + x2 + x3 + x4 + x5 + x6,
  subject to
    x1 +                                  x6    >= 7,
    x1 + x2                                     >= 20,
           x2 +  x3                             >= 14,
                   x3 +  x4                     >= 20,
                            x4 + x5             >= 10,
                                    x5 +  x6     >= 5, 
with all variables nonnegative.

We will find three alternative optimum solutions, each having the same minimum w value,
with the first solution given by minimize_lp.

minimize_lp (x1 + x2 + x3 + x4 + x5 + x6,
    [x1 + x6 >= 7, x1 + x2 >= 20,x2 +  x3 >= 14,x3 +  x4 >= 20,
                            x4 + x5 >= 10,x5 +  x6 >= 5],
     [x1,x2,x3,x4,x5,x6]);

To solve this LP directly requirea the introduction of six surplus variables (x7,x8,...,x12),
and then six artificial variables (x13,x14,...,x18), for a total of eighteen variables, 
and then employ the two phase simplex method.

In matrix form, minimize w = C^t . X such that A . X = B, X >= 0,
with X^t = (x1, x2,...,x18), C^t = (1,1,...,M) which has six 1's at the start, then six zeros,
then six M's at the end. The matrix A has six rows, and B^t = (7,20,14,20,10,5).

### Solving Ex 4 primal with 18 variables

A : matrix (
    [1,0,0,0,0,1,-1,0,0,0,0,0,1,0,0,0,0,0],
    [1,1,0,0,0,0,0,-1,0,0,0,0,0,1,0,0,0,0],
    [0,1,1,0,0,0,0,0,-1,0,0,0,0,0,1,0,0,0],
    [0,0,1,1,0,0,0,0,0,-1,0,0,0,0,0,1,0,0],
    [0,0,0,1,1,0,0,0,0,0,-1,0,0,0,0,0,1,0],
    [0,0,0,0,1,1,0,0,0,0,0,-1,0,0,0,0,0,1])$
B : cvec( [7,20,14,20,10,5])$
C : cvec ([1,1,1,1,1,1,0,0,0,0,0,0,M,M,M,M,M,M])$
Co : cvec ([M,M,M,M,M,M])$
"- Co^t . B " = - transpose (Co) . B;
" C^t - Co^t . A " = transpose(C) - transpose(Co) . A;

r7 : [1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0]$
r8 : [-2,-2,-2,-2,-2,-2,1,1,1,1,1,1,0,0,0,0,0,0,-76]$
r1 : endcons (B[1,1], A[1])$
r2 : endcons (B[2,1], A[2])$
r3 : endcons (B[3,1], A[3])$
r4 : endcons (B[4,1], A[4])$
r5 : endcons (B[5,1], A[5])$
r6 : endcons (B[6,1], A[6])$
RL : [r1,r2,r3,r4,r5,r6,r7,r8]$
vL : [x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,x13,x14,x15,x16,x17,x18,rhs,Basis]$
bL : [x13,x14,x15,x16,x17,x18,z1,z2]$
tableau (RL);

b2ratio(RL, 1);

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

b2ratio(RL, 2);

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

b2ratio(RL, 3);

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

b2ratio(RL, 4);

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

b2ratio (RL, 7);

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

b2ratio (RL, 5);

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

b2ratio (RL,6);

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

All the artificial variables have left the Basis and hence equal zero. We can delete their
columns of the tableau, and then delete the last row of zeros, otherwise preserving all 
other tableau elements.

colL : [13,14,15,16,17,18]$

r1 : remL (RL[1], colL)$
r2 : remL (RL[2], colL)$
r3 : remL (RL[3], colL)$
r4 : remL (RL[4], colL)$
r5 : remL (RL[5], colL)$
r6 : remL (RL[6], colL)$
r7 : remL (RL[7], colL)$
RL : [r1,r2,r3,r4,r5,r6,r7]$
vL : remL (vL, colL)$
bL : [x1,x5,x3,x7,x4,x6,z]$
tableau (RL);

Since the B/N tableau for minimization is actually maximizing z = -w, we have to
take the minimum w to be +45, with x1 = 20, x2 = 0, x3 = 14, x4 = 6, x5 = 4, x6 = 1.
This is a second optimum minimum for the primal, as compared with minimize_lp.

The z-row coefficients of the surplus variables x7, x8, ...,x12 gives the solution of the
dual: y1 = 0, y2 = 1, y3 = 0, y4 = 1, y5 = 0, y6 = 1.

The only surplus variable with a positive value in the solution of the primal is x7 = 14,
so the first constraint of the *primal* LP holds as an inequality, and the first decision
variable of the *dual* y1 = 0, an example of the complementary slackness principle.

Here is a matrix method of finding the solutions to the symmetric dual LP from the primal
LP optimal tableau:
The six basis variables in the optimal primal tableau are x1, x3, x4, x5, x6, x7.
Let Cb be the matrix column vector containing elements (1,2,4,5,6,7) of the primal 
vector C, and let Ab be the matrix containing columns (1,2,4,5,6,7) of the primal matrix A.
Then Y*^t = Cb^t . Ab^(-1) = transpose (Cb) . invert (Ab) = (y1, y2, y3, y4, y5, y6 ).

Cb : part (C, [1,3,4,5,6,7]);

Ab : newM (A, [1,3,4,5,6,7]);

transpose (Cb) . invert (Ab);

which is Y*^t = (y1,y2,y3,y4,y5,y6) for the optimal tableau of the symmetric dual LP.

### Solving Ex 4 symmetric dual with 12 variables

Rather than solving the primal LP, it is simpler (if one wishes to avoid using maximize_lp)
 to solve the symmetric dual program:
   maximize z = 7*y1 + 20*y2 + 14*y3 + 20*y4 + 10*y5 + 5*y6,
   subject to
     y1 + y2 <= 1,
     y2 + y3 <= 1,
     y3 + y4 <= 1,
     y4 + y5 <= 1,
     y5 + y6 <= 1,
     y1 + y6 <= 1,
 with all variables nonnegative.

maximize_lp (7*y1 + 20*y2 + 14*y3 + 20*y4 + 10*y5 + 5*y6,
    [y1 + y2 <= 1,y2 + y3 <= 1,y3 + y4 <= 1,y4 + y5 <= 1,y5 + y6 <= 1, y1 + y6 <= 1],
    [y1,y2,y3,y4,y5,y6]);

We introduce 6 slack variables, y7 through y12, on the left hand sides of the constraints.
The dual LP in matrix form is 
  maximize z = u^t . Y such that D . Y = v, with Y >= 0.

Given the Step 0 LP maximize z = u^t . Y, such that  D . Y = v, with Y >= 0,
the step 0 maximization tableau is, using these matrices:
                |           Y^t                 |    rhs       |  Basis   
              -----------------------------------------------------------
                |            D                   |   v           |  Yo      
              -----------------------------------------------------------
                |     uo^t . D - u^t        | uo^t . v    |   z

Y : cvec ([y1,y2,y3,y4,y5,y6,y7,y8,y9,y10,y11,y12])$
D : matrix (
    [1,1,0,0,0,0,1,0,0,0,0,0],
    [0,1,1,0,0,0,0,1,0,0,0,0],
    [0,0,1,1,0,0,0,0,1,0,0,0],
    [0,0,0,1,1,0,0,0,0,1,0,0],
    [0,0,0,0,1,1,0,0,0,0,1,0],
    [1,0,0,0,0,1,0,0,0,0,0,1])$
v : cvec( [1,1,1,1,1,1])$
u : cvec ([7,20,14,20,10,5,0,0,0,0,0,0])$
"maximize  ";
z = transpose(u) . Y;
"such that ";
 D . Y = v;
uo : cvec ([0,0,0,0,0,0])$
" uo^t . v " =  transpose (uo) . v;
" uo^t . D - u^t " = transpose(uo) . D - transpose(u);

lme is our alias for list_matrix_entries (in Econ1.mac).

Step 0 tableau 

YL : lme (Y);
vL : flatten (endcons([rhs,Basis], YL));
bL : endcons (z, rest (YL,6));
r1 : endcons (v[1,1], D[1])$
r2 : endcons (v[2,1], D[2])$
r3 : endcons (v[3,1], D[3])$
r4 : endcons (v[4,1], D[4])$
r5 : endcons (v[5,1], D[5])$
r6 : endcons (v[6,1], D[6])$
r7 : endcons (0, lme (-transpose(u)))$
RL : [r1,r2,r3,r4,r5,r6,r7]$
tableau (RL);

bratio (RL, 2);

Step 1  tableau:

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

bratio (RL, 4);

Step 2  tableau:

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

bratio (RL, 5);

Step 3  tableau:

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

bratio (RL, 6);

Step 4 tableau:

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

This is the optimized dual LP tableau with z* = 45, y1* = 0, y2* = 1, y3* = 0, y4* = 1,
 y5* = 0, and y6* = 1 as the solution to the dual, in agreement with the solution found by
maximize_lp.

A third alternative solution to the *primal* LP is given by the z-row coefficients of the slack 
variables y7, y8,...,12. Namely, w* = 45 with x1* = 20, x2* = 0, x3* = 15, x4* = 5, x5* = 5, 
x6* = 0.

Using the fact that the basis variables in the optimal dual tableau are y2,y4,y5,y6,y8,y12,
we can also use a matrix method to predict the same optimal solution for the primal
variables as is given by the dual slack variable z-row coefficients.
The step 0 matrix dual LP description is: max z = u^t . Y such that D . Y = v, with Y>=0.

BL : [2,4,5,6,8,12]$
ub : part (u, BL);

Db : newM (D, BL);

transpose (ub) . invert (Db);

which gives an optimal solution for the primal variables (x1,x2,x3,x4,x5,x6).

## Unsymmetric Duals

Let C^t indicate the transpose of the matrix C, etc.

A.       Given the primal LP:
    maximize z = C^t . X such that  A . X = B, X >= 0,                                   (1)
 the dual LP of (1) is:
   minimize  w = B^t . Y such that A^t . Y >= C, Y unrestricted in sign.         (2)

Conversely, the dual LP of (2) is (1).
-----------------------------------------------------------------
B.     Given the primal LP:
   minimize  w = C^t . X such that  A . X = B, X >= 0,                                     (3)
the dual LP of (3) is:
   maximize  z = B^t . Y such that A^t . Y <= C, Y unrestricted in sign.          (4)

Conversely, the dual of (4) is (3).
-------------------------------------------------------------------

Since the dual of a program in "standard form" is not itself in "standard form",
these duals are unsymmetric. Their forms are consistent with and a direct consequence
of the definition of symmetric duals. (See B/N Problem 4.8)

### Unsymmetric Duals Step 0 Notation

### Dual Form (1) - (2) Step 0 Notation

The path taken to arrive at a suitable Step 0 LP problem often requires increasing the
number of variables (by adding slack, surplus, and/or artificial variables). We then need
to expand the sizes of matrices: X --> Xs, Y --> Ys,  A --> As, etc.

Given the primal LP (1):
  maximize z = C^t . X such that  A . X = B, X >= 0,   
we write the Step 0 LP form as
  maximize z = Cs^t . Xs such that  As . Xs = B, Xs >= 0.            (1')

The dual LP of (1) is (2):
   minimize  w = B^t . Y such that A^t . Y >= C, Y unrestricted in sign.
We write the Step 0 LP form as   
   minimize  w = u^t . Ys such that D . Ys = v, Ys >= 0.           (2')

### Dual Form (3) - (4) Step 0 Notation

Given the primal LP (3):
     minimize  w = C^t . X such that  A . X = B, X >= 0,
we write the Step 0 LP form as
    minimize  w = Cs^t . Xs such that  As . Xs = B, Xs >= 0,         (3')

The dual LP of (3) is (4):
    maximize  z = B^t . Y such that A^t . Y <= C, Y unrestricted in sign.
We write the Step 0 LP form as
    maximize   z = u^t . Ys such that D . Ys = v, Ys >= 0.            (4')

### Example 5

Find the optimal solutions of the unsymmetric dual pair given and show that the
solution of one of the pair will predict a solution of the other.

Given the primal LP
  maximize z = x1 + 3*x2 - 2*x3 = C^t . X,
    subject to  4*x1 + 8*x2 + 6*x3 = 25,
                     7*x1 + 5*x2 + 9*x3 = 30,
           equivalent to A . X = B,
     with x1,x2,x3 >= 0.

The (unsymmetric) dual LP is
  minimize w = 25*y1 + 30*y2 = B^t . Y,
    subject to  4*y1 + 7*y2 >= 1,
                     8*y1 + 5*y2 >= 3,
                     6*y1 + 9*y2 >= -2,
            equivalent to A^t . Y >= C,
                 with y1 and y2 unrestricted in sign.

This pair is an example of the pair (1) & (2) above.

### Primal LP Solution

The primal LP is:
maximize z = x1 + 3*x2 - 2*x3 = C^t . X,
    subject to  4*x1 + 8*x2 + 6*x3 = 25,
                     7*x1 + 5*x2 + 9*x3 = 30,
      equivalent to A . X = B,
     with x1,x2,x3 >= 0.

maximize_lp (x1 + 3*x2 - 2*x3, [4*x1 + 8*x2 + 6*x3 = 25, 7*x1 + 5*x2 + 9*x3 = 30],
    [x1,x2,x3]);

float(%);

To determine the optimum tableau with an initial feasible solution x1 = x2 = x3 = 0, we 
add artificial variables x4 and x5 to the left-hand side of the constraints so that the 
initial feasible solution is x1 = x2 = x3 = 0, x4 = 25 and x5 = 30. The Step 0 LP is then:

maximize z = x1 + 3*x2 - 2*x3 - M*x4 - M*x5  = Cs^t . Xs,
    subject to  4*x1 + 8*x2 + 6*x3 + x4 = 25,
                     7*x1 + 5*x2 + 9*x3 + x5 = 30,
          equivalent to   As . Xs = B,
     with x1,x2,x3,x4,x5 >= 0 or Xs >= 0.

### Two Phase Simplex Method Primal LP Solution

We solve this Step 0 LP with artificial variables using the B/N "two phase simplex" 
method described in B/N, pp. 33,34.

Given the Step 0 LP:
maximize z = Cs^t . Xs, such that  As . Xs = B, with Xs >= 0,

the step 0 maximization tableau is then, using these matrices:

                |           Xs^t                 |    rhs        | Basis 
               ------------------------------------------------------------
                |            As                  |     B          |  Xso
               -----------------------------------------------------------
                |  Cso^t . As - Cs^t      | Cso^t . B   |   z

Xso is the known initial feasible solution Basis vector of symbols. 
Xso is defined using the Basis variable order in the constraint equations, and not 
necessarily with the order in Xs. 
Xso has the same number of variables as the rhs vector B.
Cso is the vector of objective coefficients, taken from Cs, associated with the initial basis 
vector Xso, and in the same order as Xso. 

The last row is the "z-row", with the middle section containing -Cs^t as the main piece. 
The first term,  Cso^t . As, is present just in case one or more of the basis vector Xso 
variables have nonzero objective coefficients. We want to start the simplex algorithm with 
all of the step 0 basis variables having zero z-row coefficients.

The second term of the z-row is Cso^t . B, which is the step 0 value of z.
In our problem, Xs^t = (x1,x2,x3,x4,x5), Xso^t = (x4,x5).

Cs : cvec([1,3,-2,-M,-M])$
Cso : cvec ([-M, -M])$
As : matrix( [4,8,6,1,0],[7,5,9,0,1] )$
B : cvec ([25, 30])$
"Cso^t . B" = transpose(Cso) . B;
"Cso^t . As - Cs^t" = transpose(Cso) . As - transpose(Cs);

Step 0 tableau; 
r1 and r2 are lists of the coefficients of the two Step 0 constraint equations, with 5 elements
corresponding to the variables x1,x2,x3,x4,x5, and the sixth element corresponding to
the right-hand side, coming from the vector B.
We split the z-row into two rows which supply the lists r3 and r4.
If the original z-row has the form (including the sixth element from the initial z value),
  [a1 + b1*M, a2 + b2*M,...], then we set r3 to the list [a1, a2, ...] and 
r4 to the list [b1, b2, ...].

For this problem,
Xs^t = (x1,x2,x3,x4,x5), Cs^t = (1,3,-2,-M,-M), z = Cs^t . Xs, Xso^t = (x4,x5), 
Cso^t = (-M,-M), B^t = (25,30), As = matrix( [4,8,6,1,0],[7,5,9,0,1] ),
Cso^t . B = -55*M, Cso^t . As = (-11*M, -13*M, -15*M, -M, -M),
Cso^t . As - Cs^t = (-1 - 11*M, -3 - 13*M, 2 - 15*M, 0, 0)
This last (5 element) row is separated into two 5 element bottom rows:
  (-1, -3, 2, 0, 0), with the 6'th element of r3 (rhs) = 0, the numerical part of from Cso^t . B,
  (-11, -13, -15, 0, 0) with the 6'th element of r4 (rhs) = -55, from  the symbolic part
of Cso^t . B).
We then apply the simplex method to four rows as a start, choosing the most negative
coefficient in the bottom row (row 4 here) to define the "work column" or "pivot column".

Step 0 tableau: (Note that, in general, the order of the Step 0 Basic variables in the list 
bL follows the order those variables appear in Xso)

vL : [x1,x2,x3,x4,x5,rhs,Basis]$
bL : [x4,x5,z1,z2]$
r1 : [4,8,6,1,0,25]$
r2 : [7,5,9,0,1,30]$
r3 : [-1,-3,2,0,0,0]$
r4 : [-11,-13,-15,0,0,-55]$
RL : [r1,r2,r3,r4]$
tableau(RL);

Choose the x3 column as the pivot column, having the most negative coefficient in the
bottom row (excluding the rhs number) . Use the minimum ratio test to select the pivot 
row, looking only at the r1 and r2 rows.

b2ratio (RL, 3);

Row 2 is the winning row.

Step 1 tableau:

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

Column 2 is the pivot column,

b2ratio(RL,2);

Row 1 is the winner of the minimum ratio test.

Step 2 tableau:

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

No more negative numbers in last row, x4 and x5 have left the Basis, so x4 = x5 = 0.
We can now remove columns 4 and 5. Then all numbers in last row equal 0, so 
we can remove the last row entirely.

The current value of the row 1 is given by RL[1].
We use our function remL(alist, item-numbers-list), which returns a depleted list missing
those items in alist referred to by the integers in item-numbers-list.

r1 : remL(RL[1], [4,5]);

r2 : remL(RL[2], [4,5]);

r3 : remL (RL[3], [4,5]);

Redefine the lists vL and bL, and RL is now a three element list, with elements
r1,r2,and r3.

vL : [x1,x2,x3,rhs,Basis]$
bL : [x2,x3,z]$
RL : [r1,r2,r3]$
tableau (RL);

bratio (RL, 1);

Step 3 tableau: 

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

float(RL);

With no negative coefficients in the bottom row (the z-row), an optimum solution is
z* = 70/9, x1* = 115/36, x2* = 55/36, x3* = 0  in agreement with maximize_lp.

Returning to the primal LP step 0 tableau with the knowledge that x1 and x2 are the
Basis variables of the primal solution we can predict the solution of the dual LP from 
the properties of the primal initial step 0 tableau.

display (Cs, As)$

The Basis variables in the optimal tableau are x1 and x2, elements 1 and 2 of Xs.
Let Cb be a column matrix vector whose elements are elements 1 and 2 of the vector Cs,
and let Ab be a matrix made from columns 1 and 2 of the matrix As. Then the elements
of the dual solution are given by
     Cb^t . Ab^(-1) =    transpose(Cb) . invert (Ab).

Cb : part(Cs,[1,2]);
Ab : newM (As,[1,2]);

transpose (Cb) . invert (Ab);

This predicts that a solution of the (unsymmetric) dual LP is y1 = 4/9, y2 = -1/9.

### Dual LP Solution

The (unsymmetric) dual is
  minimize w = 25*y1 + 30*y2 = B^t . Y,
    subject to  4*y1 + 7*y2 >= 1,
                     8*y1 + 5*y2 >= 3,
                     6*y1 + 9*y2 >= -2
          equivalent to A^t . Y >= C,
               with y1 and y2 unrestricted in sign.

minimize_lp (25*y1 + 30*y2, [4*y1 + 7*y2 >= 1,
                                             8*y1 + 5*y2 >= 3, 
                                              6*y1 + 9*y2 >= -2]);

The dual solution returned by minimize_lp is the same as was found using
Cb^t . Ab^(-1) on the primal lp result.

See B/N Prob. 2.6, 4.5, 4.7.
We need to take into account that y1 and y2 are unrestricted in sign before we can use
the simplex method.
Set y1 = y3 - y4 and y2 = y5 - y6 with y3,y4,y5,y6 >= 0. Then multiply the last constraint
by -1 to force a nonnegative right-hand side. This gives:

minimize w = 25*y3 - 25*y4 + 30*y5 - 30*y6,
    subject to  4*y3 - 4*y4 + 7*y5 - 7*y6 >= 1,
                     8*y3 - 8*y4 + 5*y5 - 5*y6 >= 3,
                    -6*y3 + 6*y4 - 9*y5 + 9*y6 <= 2

Convert to standard form by subtracting surplus variables y7 and y8, respectively, from
the left-hand sides of the first two constraints and add a slack variable y9 to the left-hand
side of the third constraint. Then add artificial variables y10 and y11 respectively to the
left-hand sides of the first two constraints. We finally arrive at:

minimize w = 25*y3 - 25*y4 + 30*y5 - 30*y6 +0*y7 + 0*y8 + 0*y9 + M*y10 + M*y11,

    (equivalent to minimize w = u^t . Ys),

    subject to:  4*y3 - 4*y4 + 7*y5 - 7*y6 - y7 +       y10      = 1,
                      8*y3 - 8*y4 + 5*y5 - 5*y6  -    y8 +        y11 = 3,
                     -6*y3 + 6*y4 - 9*y5 + 9*y6 +        y9             = 2,

      (equivalent to D . Ys = v),

  with y3,y4,y5,y6,y7,y8,y9,y10,y11 >= 0, (or Ys >= 0).

Then an initial feasible solution is, working top to bottom through the constraint equations:
  y3 = y4 = y5 = y6 = y7 = y8 = 0, y10 = 1, y11 = 3, y9 = 2.

### Two Phase Simplex Method Dual LP Solution

Using the Step 0 tableau matrix notation:

min w = u^t . Ys, such that D . Ys = v, with Ys >= 0, (u^t means transpose(u) )

  Ys^t = (y3, y4, y5, y6, y7, y8, y9, y10, y11), 
  u^t = (25, -25, 30, -30, 0, 0, 0, M, M), v^t = (1, 3, 2),
  D = matrix ( [4,-4,7,-7,-1,0,0,1,0], [8,-8,5,-5,0,-1,0,0,1], [-6,6,-9,9,0,0,1,0,0] ),
**  note that Yso is defined using the Basis variable's order in the constraint equations,
**    not in the order found in Ys,
Yso^t = (y10, y11, y9), 
uo^t = (M, M, 0) = transpose of the vector uo which contains the
objective coefficients of the variables in Yso, in the same order as these variables
appear in Yso (which is the order the basis variables appear in the constraint eqns).

The step 0 minimization tableau is, using these matrices:

                |           Ys^t                 |  rhs        | Basis            
           -----------------------------------------------------------
                |            D                   |   v            |  Yso
           -----------------------------------------------------------
                |   u^t - uo^t . D          | - uo^t . v  |    z

With artificial variables part of the Step 0 problem, we put the coefficients of M on the
last line, splitting the z-row into two rows.

u : cvec ([25,-25,30,-30,0,0, 0, M, M])$
uo : cvec ([M,M,0])$
v : cvec ([1, 3, 2])$
D : matrix ( [4,-4, 7,-7, -1,  0,  0,  1,  0], 
                  [8,-8, 5,-5,  0, -1,  0,  0, 1],
                 [-6, 6, -9, 9, 0,  0,  1,  0, 0]  )$
"- uo^t . v" = - transpose(uo) . v;
"u^t - uo^t . D" = transpose(u) - transpose(uo) . D;

The last output is the z-row of the tableau (ignoring the current value of z given by
 - uo^t . v = -4*M), except we need to split the z-row into two rows since we have 
artificial variables. The z-row splits into two rows which supply the lists r4 and r5 in
this problem.

If the original z-row has the form (including the sixth element from the initial z value),
  [a1 + b1*M, a2 + b2*M,...], then we set r4 to the list [a1, a2, ...] and 
r5 to the list [b1, b2, ...].

Since we have gone to the trouble of defining D and v, we can use these to define
r1,r2,and r3.

Step 0 tableau: (Note that the order of the Step 0 Basic variables in the list bL follows
the order those variables appear in Yso, ie., in the order of the Step 0 constraint equations.)

vL : [y3,y4,y5,y6,y7,y8,y9,y10,y11,rhs, Basis]$
bL : [y10,y11,y9, z1,z2]$
r1 : endcons (v[1,1], D[1])$
r2 : endcons (v[2,1], D[2])$
r3 : endcons (v[3,1], D[3])$
r4 : [25,-25,30,-30,0,0,0,0,0,0]$
r5 : [-12,12,-12,12,1,1,0,0,0,-4]$
RL : [r1,r2,r3,r4,r5]$
tableau (RL);

Choose the y3 column (column 1) as the pivot column, y3 enters the Basis.

b2ratio (RL, 1);

Row 1 wins the minimum ratio test.

Step 1 tableau:

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

Choose the y6 column (col 4) as the pivot column, y6 enters the basis. The only 
positive y6 coefficient is 9 in row 2, so y11 leaves the Basis.

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

y10 and y11 have both left the Basis and are equal to 0, so we can delete the y10 and
y11 columns, (columns 8 and 9) and then since we will have all zeros in the last row, 
we can delete the last row.

r1 : remL (RL[1], [8, 9]);

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

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

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

We preserve everything else about the last tableau, including the order of the rows.
Note that bL needs to have the order top to botton of the last tableau basis variables.

RL : [r1,r2,r3,r4]$
vL : [y3,y4,y5,y6,y7,y8,y9,rhs, Basis]$
bL : [y3, y6, y9, z]$
tableau (RL);

There are no more negative z-row coefficients (ignoring the rhs number), so this
tableau shows that z = -w has reached a maximum of -70/9, so w has reached a
minimum value +70/9. Recall that the B/N tableau simplex method for min w
actually maximizes z = -w.

Because the B/N tableau method finds min w by finding max z = -w, we have to take
the negative of the last number in the bottom row as the minimum of w, so
w* = 70/9. 

y1* = y3* - y4* = 4/9 - 0 = 4/9. 
y2* = y5* - y6* = 0 - 1/9 = -1/9.

w* = 25*y1 + 30*y2 = 100/9 - 30/9 = 70/9.

float(%);

The Basis variables in the solution vector Ys* are y3, y6, and y9, elements 1,4,and 7
in Ys.

Let ub be a column matrix vector whose elements are elements 1, 4, and 7 of the 
vector u, and let Db be a matrix made from columns 1, 4 and 7 of the matrix D. 
Then the elements of the primal solution are predicted to be:
     ub^t . Db^(-1) =    transpose(ub) . invert (Db)

ub : part (u, [1,4,7]);
Db : newM (D, [1,4,7]);

transpose (ub) . invert (Db);

This predicts a solution of the primal should be x1 = 115/36, x2 = 55/36, x3 = 0, as 
we found in our primal solution.

### Example 6

Determine the dual of the program
  minimize z = 3*x1  +    x2 + 0*x3 + 0*x4 +  M*x5 + M*x6,
    subject to:     x1  +    x2  -    x3    +              x5                 = 7,
                      2*x1  + 3*x2               -  x4   +               x6    =  8,
    with x1,x2,x3,x4,x5,x6 >= 0.
The (unsymmetric) dual is
  maximize w = 7*y1 + 8*y2,
    subject to  y1 + 2*y2 <= 3,
                     y1 + 3*y2 <= 1,
                     - y1 <= 0,
                     -y2  <= 0,
                      y1 <= M,
                      y2 <= M,
with no further restrictions on the sign of y1 and y2.
Now  -y1 <= 0 implies y1 >= 0, and -y2 <= 0 implies y2 >= 0, and y1 <= M simply 
requires y1 to be finite (a condition always presupposed), so we can simplify the dual to:
maximize w = 7*y1 + 8*y2,
    subject to  y1 + 2*y2 <= 3,
                     y1 + 3*y2 <= 1,
    with y1, y2 >= 0.

### Example 7

Find the optimal solutions of the unsymmetric dual pair given and show that the
solution of one of the pair will predict a solution of the other.

Given the primal LP
  minimize w = x1 + 2*x2 + x3 = C^t . X,
    subject to              x2 +  x3 = 1,
                     3*x1 + x2 + 3*x3 = 4,
           equivalent to A . X = B,
     with x1,x2,x3 >= 0.

The (unsymmetric) dual LP is
  maximize z = y1 + 4*y2 = B^t . Y,
    subject to  3*y2 <= 1,
                     y1 + y2 <= 2,
                     y1 + 3*y2 <= 1,
            equivalent to A^t . Y <= C,
                 with y1 and y2 unrestricted in sign.

This pair is an example of the unsymmetric dual pair (3) & (4) above.

### Primal LP Solution

Given the primal LP
  minimize w = x1 + 2*x2 + x3,
    subject to              x2 +  x3 = 1,
                     3*x1 + x2 + 3*x3 = 4,          
     with x1,x2,x3 >= 0.

minimize_lp (x1 + 2*x2 + x3, [x2 +  x3 = 1, 
                                            3*x1 + x2 + 3*x3 = 4],[x1,x2,x3]);

Add artificial variables x4 and x5 to the left-hand sides of the equality conditions.

minimize w = x1 + 2*x2 + x3 + M*x4 + M*x5 = Cs^t . Xs,
subject to
    x2 + x3 + x4 = 1,
   3*x1 + x2 + 3*x3 + x5 = 4, 

equivalent to As . Xs = B, with Xs >= 0.

An initial feasible solution is then x1 = x2 = x3 = 0, x4 = 1, x5  = 4.

Let Xs^t = (x1,x2,x3,x4,x5), Xso^t = (x4,x5),
Cs^t = (1,2,1,M,M),  Cso^t = (M,M), B^t = (1,4).

The step 0 minimization tableau is then, using these matrices:

                |           Xs^t                 |   rhs           | Basis      
             -----------------------------------------------------------
                |            As                   |   B             |  Xso
              -----------------------------------------------------------
                |   Cs^t - Cso^t . As      | - Cso^t . B  |   z

Cs : cvec ([1,2,1,M,M])$
Cso : cvec ([M,M])$
As : matrix ( [0,1,1,1,0], [3,1,3,0,1])$
B : cvec ([1,4])$
" - Cso^t . B " = - transpose (Cso) . B;
" Cs^t - Cso^t . As " = transpose (Cs) - transpose (Cso) . As;

r1 : endcons (B[1,1], As[1])$
r2 : endcons (B[2,1], As[2])$
r3 : [1,2,1,0,0,0]$
r4 : [-3,-2,-4,0,0,-5]$
RL : [r1,r2,r3,r4]$
vL : [x1,x2,x3,x4,x5,rhs,Basis]$
bL : [x4,x5,z1,z2]$
tableau (RL);

b2ratio (RL, 3);

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

b2ratio (RL, 1);

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

No more negative numbers in last row, x4 and x5 have left the Basis, so x4 = x5 = 0.
We can now remove columns 4 and 5. Then all numbers in last row equal 0, so 
we can remove the last row entirely.

The current value of the row 1 is given by RL[1].
We use our function remL(alist, item-numbers-list), which returns a depleted list missing
those items in alist referred to by the integers in item-numbers-list.

We retain the rest of the tableau. The list bL needs to be redefined to reflect the names
of the Basis variables in this last tableau.

r1 : remL (RL[1], [4,5]);

r2 : remL (RL[2], [4,5]);

r3 : remL (RL[3], [4,5]);

RL : [r1,r2,r3]$
vL : [x1,x2,x3,rhs, Basis]$
bL : [x3,x1,z]$
tableau (RL);

The optimum primal solution has x1* = 1/3, x2* = 0, x3* = 1, min w = w* = 4/3, 
the negative of z. Our minimization tableau maximizes z = -w.

The Basis variables in the optimum tableau are x1 and x3, elements 1 and 3 of Xs.

Let Cb be a column vector with elements 1 and 3 from Cs, and Ab be the matrix with
columns 1 and 3 chosen from As. We then calculate Cb^t . Ab^(-1) to get a prediction
 of values of y1 and y2 for the dual LP.

display (Cs, As)$

Cb : part(Cs,[1,3]);
Ab : newM (As,[1,3]);

transpose (Cb) . invert (Ab);

Thus the solution to the primal LP, implies a solution to the dual LP is y1 = 0, y2 = 1/3.

What does maximize_lp give as a solution of the dual LP?

maximize_lp (y1 + 4*y2,[3*y2 <= 1,y1 + y2 <= 2,y1 + 3*y2 <= 1]);

which agrees with the prediction.

### Dual LP Solution

The (unsymmetric) dual LP is
  maximize z = y1 + 4*y2,
    subject to  3*y2 <= 1,
                     y1 + y2 <= 2,
                     y1 + 3*y2 <= 1,            
                 with y1 and y2 unrestricted in sign.

Let y1 = y3 - y4, and let y2 = y5 - y6.
Add slack variables y7,y8, y9 to the left-hand sides of the constraint inequalities
to get constraint equations. We then have:

maximize z = y3 - y4 + 4*y5 - 4*y6 + 0*y7 + 0*y8 + 0*y9 = u^t . Ys,
subject to
     3*y5 - 3*y6 + y7 = 1,
   y3 - y4 + y5 - y6 + y8 = 2,
   y3 - y4 + 3*y5 - 3*y6 + y9 = 1,
   equivalent to D . Ys = v,
with y3,y4,y5,y6,y7,y8,y9 >= 0, or Ys >= 0.

The initial feasible solution is
y3 = y4 = y5 = y6 = 0, y7 = 1, y8 = 2, y9 = 1.
Ys^t = (y3,y4,y5,y6,y7,y8,y9), Yso^t = (y7,y8,y9), v^t = (1,2,1),
u^t = (1,-1,4,-4,0,0,0), uo^t = (0,0,0).

Given the LP:
maximize z = u^t . Ys, such that  D . Ys = v, with Ys >= 0,
the step 0 maximization tableau is, using these matrices:

                |           Ys^t                 |  rhs      |  Basis     
            -----------------------------------------------------------
                |            D                   |   v          |  Yso
            -----------------------------------------------------------
                |   uo^t . D - u^t          | uo^t . v   |   z

u : cvec ([1,-1,4,-4,0,0,0])$
v : cvec ([1,2,1])$
uo : cvec ([0,0,0])$
D : matrix ([0,0,3,-3,1,0,0], [1,-1,1,-1,0,1,0], [1,-1,3,-3,0,0,1])$
" uo^t . D - u^t " = transpose(uo) . D - transpose(u);
" uo^t . v " = transpose(uo) . v;

r1 : endcons (v[1,1], D[1])$
r2 : endcons (v[2,1], D[2])$
r3 : endcons (v[3,1], D[3])$
r4 : [-1,1,-4,4,0,0,0,0]$
vL : [y3,y4,y5,y6,y7,y8,y9,rhs,Basis]$
bL : [y7,y8,y9,z]$
RL : [r1,r2,r3,r4]$
tableau (RL);

bratio (RL, 3);

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

bratio (RL,1);

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

With no more negative coefficients in the bottom z-row, the optimum found is
z* = 4/3, with y3* = y4* =y6* = y7* = y9* = 0, y5* = 1/3, y8* = 5/3.
y1* = y3* - y4* = 0, y2* = y5* - y6* = 1/3, which agrees with the solution found by
maximize_lp.

The Basis variables in the optimal tableau are y3, y5, and y8, respectively elements
1, 3, and 6 in Ys.
Let ub be a matrix column vector made up of elements 1, 3 and 6 in u,
Let Db be a matrix made up of columns 1, 3 and 6 of D.
Then a predicted solution of the primal is ub^t . Db^(-1) = transpose(ub) . invert(Db).

display (u, D)$

ub : part (u, [1, 3, 6]);
Db : newM (D, [1, 3, 6]);

transpose (ub) . invert (Db);

A predicted solution of the primal is then  x1 = 1/3, x2 = 0, x3 = 1, which agrees with
our simplex solution of the primal LP.

## General Rules For Forming Duals

Butenko's Video 28 (Dual of the Dual LP and Rules for Forming the Dual) considers
general rules for forming duals, with the table:

matrix (["Primal (max )", |,"Dual (min)"], ["==========",|, "=========="],
     ["(=) Constraint", |,"Unrestricted Sign Var."],["--------------",|, "--------------"],
    [" (<=) Constraint",|, "Non-negative Var."], ["--------------", |,"--------------"],
    ["Unrestricted Sign Var.", |," (=) Constraint"],["--------------",|, "--------------"],
    ["Non-negative Var.",|, " (>=) Constraint"] );

In using this table, we associate the first dual variable y1 with the first constraint,
the dual variable y2 with the second constraint, etc. Thus if the second primal constraint
condition was an equality (=) constraint, y2 would be unrestricted in sign.

### Example 8

Given the primal LP
maximize z = 2*x1 + x2
  such that 
     x1 + x2 = 2,
     2*x1 - x2 >= 3,
    x1 - x2 <= 1,
  with  x1 >= 0, and x2 unrestricted (-urs).

A first step is to put the primal LP in "canonical form" for maximization, in which all
inequality constraints have the form (<=). We then need to multiply both sides of the
second constraint by (-1) which will reverse the sense of the inequality.

maximize z = 2*x1 + x2
  such that 
     x1 + x2 = 2,
     -2*x1 + x2 <= -3,
    x1 - x2 <= 1,
  with  x1 >= 0, and x2 unrestricted (-urs).

Because the first primal constraint is an equality constraint, the first dual variable y1
will be unrestricted in sign. The second and third dual variables (y2,y3) will be 
non-negative variables since the second and third primal constraints are (<=)
inequality constraints.

Since x1 >= 0 in the primal, the dual constraint condition associated with the x1 
coefficients will be a (>=) type of inequality constraint for minimization. Since x2 is -urs, 
the dual constraint condition associated with the x2 coefficients will be an equality (=)
type of condition.

We then have the dual LP:

  minimize w = 2*y1 - 3*y2 + y3,
    such that
              y1 - 2*y2 + y3 >= 2,
              y1  + y2  - y3 = 1,
    with y1 unrestricted, and y2,y3 >= 0.

### Example 9

Assume the primal LP is:

minimize w = 0.4*x1 + 0.5*x2,
 subject to
   0.3*x1 + 0.1*x2 <= 2.7,
  0.5*x1 + 0.5*x2 = 6,
  0.6*x1 + 0.4*x2 >= 6,
with x1, x2 >= 0.

### Finding the Dual LP

To find the dual LP, we first put the primal LP in canonical form for minimization,
with (>=) type constraints (or equality) by multiplying the first constraint condition by
(-1) which will change the sense of the inequality condition.

minimize w = 0.4*x1 + 0.5*x2,
 subject to
   - 0.3*x1 - 0.1*x2 >= - 2.7,
     0.5*x1 + 0.5*x2 = 6,
     0.6*x1 + 0.4*x2 >= 6,
with x1, x2 >= 0.

To find the dual LP we associate y1 with the first constraint of the canonical form,
y2 with the second, y3 with the third. Since the second constraint is an equality
condition, y2 will be unrestricted in sign (urs). Since the first and third constraints
for minimization are (>=) type, y1 and y3 will be non-negative variables.

Since both x1 and x2 are non-negative variables, the dual constraint conditions arising
from both the x1 and x2 coefficients will be of type (<=).

We then have for the dual LP:
maximize z = -2.7*y1 + 6*y2 + 6*y3,
  subject to
   -0.3*y1 + 0.5*y2 + 0.6*y3 <= 0.4,
   -0.1*y1 + 0.5*y2 + 0.4*y3 <= 0.5,
with y1, y3 >= 0, and y2 unrestricted in sign (urs).

### Simplex Tableau Solution to Primal LP

The primal LP is:
minimize w = 0.4*x1 + 0.5*x2,
 subject to
   0.3*x1 + 0.1*x2 <= 2.7,
  0.5*x1 + 0.5*x2 = 6,
  0.6*x1 + 0.4*x2 >= 6,
with x1, x2 >= 0.

minimize_lp (0.4*x1 + 0.5*x2, [0.3*x1 + 0.1*x2 <= 2.7,
                                              0.5*x1 + 0.5*x2 = 6,
                                              0.6*x1 + 0.4*x2 >= 6], [x1,x2]);

Convert to Step 0 form to apply the simplex method.
Add slack variable x3 to left-hand side of first constraint, artificial variable x4 to left-hand
side of the second constraint, and subtract surplus variable x5 from the lhs of 3'rd
constraint, add artificial variable x6 to lhs 3'rd constraint, to have three constraint 
equations, add M*x4 + M*x6 penalty terms to the rhs of the objective.

minimize w = 0.4*x1 + 0.5*x2 + 0*x3 + M*x4 + 0*x5 + M*x6,
  equivalent to min w = Cs^t . Xs,
 subject to
   0.3*x1 + 0.1*x2 + x3 = 2.7,
  0.5*x1 + 0.5*x2 + x4 = 6,
  0.6*x1 + 0.4*x2 - x5 + x6 = 6,
equivalent to  As . Xs  = B,
with x1, x2, x3, x4, x5, x6  >= 0,  or Xs >= 0.

An initial (Step 0) feasible solution is
x1 = x2 =  x5 = 0, x3 = 2.7, x4 = 6, x6 = 6.
Xs^t = (x1,x2,x3,x4,x5,x6), Xo^t = (x3, x4, x6),
Cs^t = (0.4, 0.5, 0, M, 0, M), Cso^t = (0, M, M)

Xo = matrix column vector of the Step 0 Basis variable symbols, in the order these
variables appear in the constraint equations.
Cso = matrix column vector of the objective coefficients of the initial Basis variables Xo,
in the same order as the symbols in Xo.
Given the Step 0 LP minimize w = Cs^t . Xs, such that As . Xs = B, with Xs >= 0,
the step 0 minimization tableau is, using these matrices:

                |           Xs^t                 |    rhs            |  Basis          
              ----------------------------------------------------------------
                |            As                   |    B              |   Xo
              ----------------------------------------------------------------
                |   Cs^t - Cso^t . As      | - Cso^t . B   |   z

Cs : cvec([0.4, 0.5, 0, M, 0, M])$
Cso : cvec ([0, M, M])$
As : matrix( [0.3, 0.1, 1, 0, 0, 0],[0.5, 0.5, 0, 1, 0, 0], [0.6, 0.4, 0, 0, -1, 1] )$
B : cvec ([2.7, 6, 6])$
"-Cso^t . B" = - transpose(Cso) . B;
" Cs^t - Cso^t . As " =  transpose(Cs) - transpose(Cso) . As;

vL : [x1, x2, x3, x4, x5, x6, rhs, Basis]$
bL : [x3, x4, x6, z1, z2]$
r1 : endcons (B[1,1], As[1])$
r2 : endcons (B[2,1], As[2])$
r3 : endcons (B[3,1], As[3])$
r4 : [0.4, 0.5, 0, 0, 0, 0, 0]$
r5 : [-1.1, -0.9, 0, 0, 1, 0, -12]$
RL : [r1,r2,r3,r4,r5]$
tableau (RL);

b2ratio (RL, 1);

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

b2ratio (RL, 2);

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

b2ratio (RL, 3);

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

Artificial variables x4 and x6 have left the Basis and are equal to zero. Remove columns
4 and 6 and preserve the remaining columns. Remove the last row.

r1 : remL (RL[1], [4, 6]);

r2 : remL (RL[2], [4, 6]);

r3 : remL (RL[3], [4, 6]);

r4 : remL (RL[4], [4, 6]);

vL : [x1, x2, x3, x5, rhs, Basis]$
bL : [x1, x3, x2, z]$
RL : [r1,r2,r3,r4]$
tableau (RL);

bratio (RL, 4);

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

No more negative coefficients in the z-row (omitting current value of z), so a primal 
optimum solution is z* = -w* = - 21/4 = -5.25, w* = 21/4 = 5.25, x1* = 15/2 = 7.5, 
x2* = 9/2 = 4.5, which agrees with the solution found by minimize_lp.

The Basis variables in the optimum tableau are (x1,x2,x5) which are elements (1,2,5) of
Xs. Let Cb be elements (1,2,5) of Cs and Ab be a matrix formed from columns (1,2,5) of
the matrix As.

Cs;

Cb : part (Cs, [1,2,5]);

As;

Ab : newM (As, [1,2,5]);

transpose (Cb) . invert (Ab);

This predicts a solution of the dual LP is y1 = -0.5, y2 = 1.1, y3 = 0.0, with
z = -2.7*y1 + 6*y2 + 6*y3 = 7.95. This prediction of the dual LP solution does not satisfy
the dual LP conditions y1, y3 >= 0, y2 unrestricted in sign (urs).

### Simplex Tableau Solution to Dual LP

The dual LP is:
maximize z = -2.7*y1 + 6*y2 + 6*y3,
  subject to
   -0.3*y1 + 0.5*y2 + 0.6*y3 <= 0.4,
   -0.1*y1 + 0.5*y2 + 0.4*y3 <= 0.5,
with y1, y3 >= 0, and y2 unrestricted in sign (urs).

maximize_lp (-2.7*y1 + 6*y2 + 6*y3, [-0.3*y1 + 0.5*y2 + 0.6*y3 <= 0.4,
   -0.1*y1 + 0.5*y2 + 0.4*y3 <= 0.5], [y1,y3]);

Since y2 is unrestricted in sign, set y2 = y4 - y5 with y4,y5 >= 0.
Add slack y6 to left-hand side of first constraint, slack y7 to left-hand side of second
constraint to convert constraints to equations.
We then have Y^t = (y1,y3,y4,y5,y6,y7), max z = u^t . Y with u^t = (-2.7, 6, 6, -6, 0, 0)
subject to D . Y = v, with Y >= 0. v^t = (0.4, 0.5), initial feasible solution:
y1 = y3 = y4 = y5 = 0, y6 = 0.4, y7 = 0.5, Yo^t = (y6,y7), uo^t = (0,0)

Given the Step 0 LP: maximize z = u^t . Y, such that  D . Y = v, with Y >= 0,
the step 0 maximization tableau is, using these matrices:
                |           Y^t              |    rhs     |  Basis   
              -----------------------------------------------------------
                |            D               |     v        |  Xso      
              -----------------------------------------------------------
                |   uo^t . D - u^t      | uo^t . v   |   z

u : cvec ([-2.7, 6, 6, -6, 0, 0])$
uo : cvec ([0,0])$
D : matrix ( [-0.3, 0.6, 0.5, -0.5, 1, 0], [-0.1, 0.4, 0.5, -0.5, 0, 1])$
Y : cvec ([y1,y3,y4,y5,y6,y7])$
v : cvec ([0.4, 0.5])$
"maximize z" = transpose(u) . Y;
"subject to";
D . Y  = v;
"uo^t . D - u^t " = transpose(uo) . D - transpose(u);
"uo^t . v " = transpose (uo) . v;

vL : [y1,y3,y4,y5,y6,y7, rhs, Basis]$
bL : [y6,y7,z]$
r1 : endcons (v[1,1], D[1])$
r2 : endcons (v[2,1], D[2])$
r3 : [2.7,-6,-6,6,0,0, 0]$
RL : [r1,r2,r3]$
tableau (RL);

bratio(RL, 2);

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

bratio (RL, 3);

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

bratio (RL, 1);

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

Optimum dual LP tableau with z = 21/4 = 5.25, y1 = 1/5 = 0.5, y2 = y4 - y5 = 11/10 = 1.1,
y3 = 0.0, agrees with maximize_lp solution.

ub : part (u, [1,3])$
Db : newM (D, [1,3])$
transpose (ub) . invert (Db);

This agrees with the optimal primal LP solution found x1 = 7.5, x2 = 4.5, and satisfies the
primal conditions x1,x2 >= 0.

## Appendix: Step 0 Tableaux

Xso is the known initial feasible solution Basis vector of symbols. 
Xso is defined using the Basis variable order in the constraint equations, and not 
necessarily with the order in Xs. 
Xso has the same number of variables as the rhs vector E.
Cso is the vector of objective coefficients, taken from Cs, associated with the initial basis 
vector Xso, and in the same order as Xso. 

Given the Step 0 LP: minimize w = Cs^t . Xs, such that As . Xs = E, with Xs >= 0,

the step 0 minimization tableau is, using these matrices:
                |           Xs^t                 |    rhs           |  Basis          
              -----------------------------------------------------------
                |            As                   |    E              |   Xso
              -----------------------------------------------------------
                |   Cs^t - Cso^t . As      | - Cso^t . E    |  z

Given the Step 0 LP: maximize z = Cs^t . Xs, such that  As . Xs = E, with Xs >= 0,

the step 0 maximization tableau is, using these matrices:
                |           Xs^t                 |    rhs          |  Basis   
              -----------------------------------------------------------
                |            As                   |     E           |  Xso      
              -----------------------------------------------------------
                |   Cso^t . As - Cs^t      | Cso^t . E   |   z