-
- home
-
- PDF (letter size)
-
- PDF (legal size)
-
- nma_matlab.m source code
Basic Matlab implementation of the Simplex matrix
algorithm
June 22, 2020 Compiled on June 22, 2020 at 1:53pm
Contents
List of Figures
List of Tables
1 Introduction
This is a description of a Matlab function called nma_simplex.m that implements the matrix
based simplex algorithm for solving standard form linear programming problem. It supports
phase one and phase two.
The function solves (returns the optimal solution \(x^{\ast }\) of the standard linear programming problem
given by\[ \min _x J(x) = c^T x \] Subject to \begin{align*} Ax &= b\\ x & \geq 0 \end{align*}
The constraints have to be in standard form (equality), which results after adding any needed
surplus and/or slack variables. The above is equivalent to Matlab’s \(A_{eq},b_{eq}\) used with the standard
command linprog.
This function returns the final tableau, which contains the final solution. It can print all of the
intermediate tableau generated and the basic feasible solutions generated during the process by
passing an extra flag argument. These are generated as it runs through the simplex
algorithm.
The final tableau contains the optimal solution \(x^{\ast }\) which can be read directly from the tableau.
Examples below illustrate how to call this function and how to read the solution from the final
tableau.
The tableau printed on the screen have this format
The optimal \(x^{\ast }\) is read directly by looking at the columns in \(A\) that make up the identity matrix.
With the debug flag set, the optimal \(x^{\ast }\) is also displayed on the screen.
For example, if the final tableau was the following
\(\begin{pmatrix} 1 & 0 & 0.4 & -0.2 & 1.4\\ 0 & 1 & -0.2 & 0.6 & 3.8\\ 0 & 0 & 0.6 & 2.2 & 24.6 \end{pmatrix} \)
Then \(x_1=1.4\) and \(x_2=3.8\). This means the optimal solution is \[ x^{\ast }= \begin{pmatrix} 1.4\\ 3.8\\ 0\\ 0 \end{pmatrix} \]
The function accepts \(A,b,c\) and a fourth parameter which is a flag ( true or false). If the flag is true,
then each tableau is printed as the algorithm searches for the optimal \(x\) solution, and it also prints
each \(x\) found at each step.
The following are few example showing how to use this function to solve linear programming
problems, and comparing the answer to Matlab’s linprog to verify they are the same.
2 Examples
2.1 Example 1
A vendor selling rings and bracelets. A ring has 3 oz. gold., 1 oz. silver. Bracelet has 1 oz. gold, 2
oz. silver. Profit on a ring is $4 and the profit on bracelet is $5. Initially we have 8 oz. gold and 9
0z silver. How many rings and bracelets to produce to maximize profit? Note: we need integer
LP to solve this. But for now we can ignore this to illustrate the use of this function. \begin{align*} x_{1} & =\text{number of rings}\\ x_{2} & =\text{number of bracelets}\\ J\left (x\right ) & =4x_{1}+5x_{2} \end{align*}
Since we want to maximize \(J(x)\), then we change the sign \[ J(x) = -4x_1 - 5x_2 \] With \(x_{i}\geq 0\). The constraints are \(3x_1+x_2 \leq 8\) and \(x_1+2 x_2\leq 9\). We start
by converting to standard linear programming model by adding slack variables. This results in
standard model\[ \min _{x} c^T x =\min _x \begin{pmatrix} -4 & -5 & 0 & 0 \end{pmatrix}\begin{pmatrix} x_{1}\\ x_{2}\\ x_{3}\\ x_{4}\end{pmatrix} \] subject to \begin{align*} Ax & =b \\\begin{pmatrix} 3 & 1 & 1 & 0\\ 1 & 2 & 0 & 1 \end{pmatrix}\begin{pmatrix} x_{1}\\ x_{2}\\ x_{3}\\ x_{4}\end{pmatrix} & =\begin{pmatrix} 8\\ 9 \end{pmatrix} \end{align*}
Here is the call and result returned which is the final tableau
Which returns
ans = 1.0000 0 0.4000 -0.2000 1.4000
0 1.0000 -0.2000 0.6000 3.8000
0 0 0.6000 2.2000 24.6000
solution_vector =
1.4000
3.8000
0
0
We see from above that the optimal \(x\) is \begin{align*} \begin{pmatrix} x_{1}\\ x_{2}\\ x_{3}\\ x_{4}\end{pmatrix} = & \begin{pmatrix} 1.4000\\ 3.8000\\ 0\\ 0\end{pmatrix} \end{align*}
Since \( J(x) = 4 x_1 + 5x_2 \) then we calculate \(J(x)=4(1.4)+5(3.8)\) which gives the optimal objective function of \(24.4\) dollars. This mean the
vendor should make \(1.4\) rings and \(3.8\) bracelets for maximum profit.
The same solution using Matlab’s linprog is
which gives
Optimization terminated.X = 1.4000 3.8000 0 0FVAL = -24.6000
EXITFLAG =
1
OUTPUT =
iterations: 0
algorithm: 'simplex'
cgiterations: []
message: 'Optimization terminated.'
constrviolation: 8.8818e-16
firstorderopt: 8.8818e-16
Which is the same.
2.2 Example 2
Given \begin{align*} A_{eq} &= \begin{pmatrix} 2 & 2 & 1 & 0 & 0\\ 2 & 2 & 0 & -1& 0\\ -1.5 & 1 & 0 & 0 & -1 \end{pmatrix} \end{align*}
And \begin{align*} b_{eq} &= \begin{pmatrix} 10&4&0 \end{pmatrix} \end{align*}
And \begin{align*} c^T &= \begin{pmatrix} \frac{1}{30}&\frac{1}{15}&0&0&0 \end{pmatrix} \end{align*}
The following gives the solution
With the output printed on the console as
tab = 0 0 1.0000 1.0000 0 6.0000
1.0000 0 0 -0.2000 0.4000 0.8000
0 1.0000 0 -0.3000 -0.4000 1.2000
0.0333 0.0667 0 0 0 0
solution_vector =
0.8000
1.2000
6.0000
0
0
From the above, we see that the solution is \begin{align*} \begin{pmatrix} 0.8\\ 1.2\\ 6\\ 0\\ 0 \end{pmatrix} \end{align*}
Given the optimal solution, the optimal objective function is now known.
To see each step and each \(x\) solution found, set the last argument to true.
Here is the result
>>>>Current tableau [phase one]
2.0000 2.0000 1.0000 0 0 1.0000 0 0 10.0000
2.0000 2.0000 0 -1.0000 0 0 1.0000 0 4.0000
-1.5000 1.0000 0 0 -1.0000 0 0 1.0000 0
0 0 0 0 0 1.0000 1.0000 1.0000 0
***********************
Current tableau [phase one]
2.0000 2.0000 1.0000 0 0 1.0000 0 0 10.0000
2.0000 2.0000 0 -1.0000 0 0 1.0000 0 4.0000
-1.5000 1.0000 0 0 -1.0000 0 0 1.0000 0
-2.5000 -5.0000 -1.0000 1.0000 1.0000 0 0 0 0
pivot row is 3
current basic feasible solution is
0
0
0
0
0
10
4
0
***********************
Current tableau [phase one]
5.0000 0 1.0000 0 2.0000 1.0000 0 -2.0000 10.0000
5.0000 0 0 -1.0000 2.0000 0 1.0000 -2.0000 4.0000
-1.5000 1.0000 0 0 -1.0000 0 0 1.0000 0
-10.0000 0 -1.0000 1.0000 -4.0000 0 0 5.0000 0
pivot row is 2
current basic feasible solution is
0.8000
1.2000
0
0
0
6.0000
0
0
***********************
Current tableau [phase one]
0 0 1.0000 1.0000 0 1.0000 -1.0000 0 6.0000
1.0000 0 0 -0.2000 0.4000 0 0.2000 -0.4000 0.8000
0 1.0000 0 -0.3000 -0.4000 0 0.3000 0.4000 1.2000
0 0 -1.0000 -1.0000 0 0 2.0000 1.0000 8.0000
pivot row is 1
current basic feasible solution is
0.8000
1.2000
6.0000
0
0
0
0
0
***********************
Current tableau [phase one]
0 0 1.0000 1.0000 0 1.0000 -1.0000 0 6.0000
1.0000 0 0 -0.2000 0.4000 0 0.2000 -0.4000 0.8000
0 1.0000 0 -0.3000 -0.4000 0 0.3000 0.4000 1.2000
0 0 0 0 0 1.0000 1.0000 1.0000 14.0000
***********************
Current tableau [phase two]
0 0 1.0000 1.0000 0 6.0000
1.0000 0 0 -0.2000 0.4000 0.8000
0 1.0000 0 -0.3000 -0.4000 1.2000
0.0333 0.0667 0 0 0 0
Using Matlab linprog
Which gives
Optimization terminated.X = 0.8000 1.2000 6.0000 0 0FVAL = 0.1067
EXITFLAG =
1
OUTPUT =
iterations: 0
algorithm: 'simplex'
cgiterations: []
message: 'Optimization terminated.'
constrviolation: 0
firstorderopt: 0
Which is the same solution.
2.3 Example 3
minimize \(2 x_1 + 3 x_2\) subject to \begin{align*} 4 x_1+2 x_2 &\geq 12\\ x_1+4 x_2 &\geq 6 \\ x_i \geq 0 \end{align*}
We convert the problem to standard form, which results in minimize \(2 x_1 + 3 x_2\) subject to \begin{align*} 4 x_1+2 x_2 -x_3 &= 12\\ x_1+4 x_2 -x_4 &=6 \end{align*}
with \(x_i \geq 0\). Therefore \begin{align*} A =& \begin{pmatrix} 4&2&-1&0\\ 1&4&0&-1 \end{pmatrix} \end{align*}
And \begin{align*} b =& \begin{pmatrix} 12\\ 6 \end{pmatrix} \end{align*}
And \(c^T = \begin{pmatrix} 2&3&0&0 \end{pmatrix}\). To solve this using nma_simplex the commands are
And the solution is
>>>>Current tableau [phase one] 4 2 -1 0 1 0 12
1 4 0 -1 0 1 6
0 0 0 0 1 1 0
***********************
Current tableau [phase one]
4 2 -1 0 1 0 12
1 4 0 -1 0 1 6
-5 -6 1 1 0 0 0
pivot row is 2
current basic feasible solution is
0
1.5000
0
0
9.0000
0
***********************
Current tableau [phase one]
3.5000 0 -1.0000 0.5000 1.0000 -0.5000 9.0000
0.2500 1.0000 0 -0.2500 0 0.2500 1.5000
-3.5000 0 1.0000 -0.5000 0 1.5000 9.0000
pivot row is 1
current basic feasible solution is
2.5714
0.8571
0
0
0
0
***********************
Current tableau [phase one]
1.0000 0 -0.2857 0.1429 0.2857 -0.1429 2.5714
0 1.0000 0.0714 -0.2857 -0.0714 0.2857 0.8571
0 0 0 0 1.0000 1.0000 18.0000
***********************
Current tableau [phase two]
1.0000 0 -0.2857 0.1429 2.5714
0 1.0000 0.0714 -0.2857 0.8571
2.0000 3.0000 0 0 0
solution_vector =
2.5714
0.8571
0
0
We see from the last tableau that \(x_1=2.5714\) and \(x_2=0.8671\). The optimal \(x\) is also printed in the display since the the
debug flag is true.
The optimal objective function is \begin{align*} J(x)&=2(2.5714)+3(0.8671) \\ &=7.7441 \end{align*}
2.4 Example 4
maximize \(3 x_1 + 5 x_2\) subject to
\begin{align*} x_1 + 5 x_2 & \leq 40 \\ 2 x_1 + x_2 & \leq 20 \\ x_1 + x_2 & \leq 12 \\ x_i \geq 0 \end{align*}
Introducing slack and surplus variables and converting to standard form gives
minimize \(- 3 x_1 - 5 x_2\) subject to \begin{align*} x_1 + 5 x_2 + x_3 &= 40 \\ 2 x_1 + x_2 + x_4 &= 20 \\ x_1 + x_2 + x_5 &= 12 \\ x_i \geq 0 \end{align*}
Therefore \begin{align*} A =& \begin{pmatrix} 1&5&1&0&0\\ 2&1&0&1&0\\ 1&1&0&0&1 \end{pmatrix} \end{align*}
And \begin{align*} b =& \begin{pmatrix} 40\\ 20\\ 12 \end{pmatrix} \end{align*}
And \(c^T = \begin{pmatrix} -2&-5&0&0&0 \end{pmatrix}\). To solve this using nma_simplex
And the solution is
>>>>Current tableau [phase one] 1 5 1 0 0 1 0 0 40
2 1 0 1 0 0 1 0 20
1 1 0 0 1 0 0 1 12
0 0 0 0 0 1 1 1 0
***********************
Current tableau [phase one]
1 5 1 0 0 1 0 0 40
2 1 0 1 0 0 1 0 20
1 1 0 0 1 0 0 1 12
-4 -7 -1 -1 -1 0 0 0 0
pivot row is 1
current basic feasible solution is
0
8
0
0
0
0
12
4
***********************
Current tableau [phase one]
0.2000 1.0000 0.2000 0 0 0.2000 0 0 8.0000
1.8000 0 -0.2000 1.0000 0 -0.2000 1.0000 0 12.0000
0.8000 0 -0.2000 0 1.0000 -0.2000 0 1.0000 4.0000
-2.6000 0 0.4000 -1.0000 -1.0000 1.4000 0 0 56.0000
pivot row is 3
current basic feasible solution is
5
7
0
0
0
0
3
0
***********************
Current tableau [phase one]
0 1.0000 0.2500 0 -0.2500 0.2500 0 -0.2500 7.0000
0 0 0.2500 1.0000 -2.2500 0.2500 1.0000 -2.2500 3.0000
1.0000 0 -0.2500 0 1.2500 -0.2500 0 1.2500 5.0000
0 0 -0.2500 -1.0000 2.2500 0.7500 0 3.2500 69.0000
pivot row is 2
current basic feasible solution is
5
7
0
3
0
0
0
0
***********************
Current tableau [phase one]
0 1.0000 0.2500 0 -0.2500 0.2500 0 -0.2500 7.0000
0 0 0.2500 1.0000 -2.2500 0.2500 1.0000 -2.2500 3.0000
1.0000 0 -0.2500 0 1.2500 -0.2500 0 1.2500 5.0000
0 0 0.0000 0 -0.0000 1.0000 1.0000 1.0000 72.0000
pivot row is 3
current basic feasible solution is
0
8
0
12
4
0
0
0
***********************
Current tableau [phase one]
0.2000 1.0000 0.2000 0 0 0.2000 0 0 8.0000
1.8000 0 -0.2000 1.0000 0 -0.2000 1.0000 0 12.0000
0.8000 0 -0.2000 0 1.0000 -0.2000 0 1.0000 4.0000
0.0000 0 0.0000 0 0 1.0000 1.0000 1.0000 72.0000
***********************
Current tableau [phase two]
0.2000 1.0000 0.2000 0 0 8.0000
1.8000 0 -0.2000 1.0000 0 12.0000
0.8000 0 -0.2000 0 1.0000 4.0000
-3.0000 -5.0000 0 0 0 0
pivot row is 1
current basic feasible solution is
0
8
0
12
4
***********************
Current tableau [phase two]
0.2000 1.0000 0.2000 0 0 8.0000
1.8000 0 -0.2000 1.0000 0 12.0000
0.8000 0 -0.2000 0 1.0000 4.0000
-2.0000 0 1.0000 0 0 40.0000
pivot row is 3
current basic feasible solution is
5
7
0
3
0
***********************
Current tableau [phase two]
0 1.0000 0.2500 0 -0.2500 7.0000
0 0 0.2500 1.0000 -2.2500 3.0000
1.0000 0 -0.2500 0 1.2500 5.0000
0 0 0.5000 0 2.5000 50.0000
solution_vector =
5
7
0
3
0
We see that \(x_1=5\) and \(x_2=7\) and \(x_4=3\) with all others zero. We always read the solution from the identity matrix
inside the final tableau. All other \(x_i\) are zero. Hence the optimal solution is \[ \begin{pmatrix} 5\\ 7\\ 0\\ 3\\ 0 \end{pmatrix} \] And the corresponding
optimal objective function is \(3 x_1+5 x_2=3(5)+5(7)=50\)
Result using Matlab linprog is
A=[1,5,1,0,0; 2,1,0,1,0; 1,1,0,0,1];b=[40,20,12];c=[-3,-5,0,0,0];
options = optimset('LargeScale','off','Simplex','on');
[X,FVAL,EXITFLAG,OUTPUT]=linprog(c,[],[], A,b,zeros(size(c)),[],[],options)
Optimization terminated.
X =
5
7
0
3
0
FVAL =
-50
EXITFLAG =
1
OUTPUT =
iterations: 0
algorithm: 'simplex'
cgiterations: []
message: 'Optimization terminated.'
constrviolation: 0
firstorderopt: 0
2.5 Example 5
maximize \(2 x_1 + 5 x_2\) subject to
\begin{align*} x_1 & \leq 4 \\ x_2 & \leq 6 \\ x_1 + x_2 & \leq 8 \\ x_i \geq 0 \end{align*}
Introducing slack and surplus variables and converting to standard form we now have the
following problem
minimize \(- 2 x_1 - 5 x_2\) subject to \begin{align*} x_1 + x_3 &= 4 \\ x_2 + x_4 &= 6 \\ x_1 + x_2 + x_5 &= 8 \\ x_i \geq 0 \end{align*}
Therefore \begin{align*} A =& \begin{pmatrix} 1&0&1&0&0\\ 0&1&0&1&0\\ 1&1&0&0&1 \end{pmatrix} \end{align*}
And \begin{align*} b =& \begin{pmatrix} 4\\ 6\\ 8 \end{pmatrix} \end{align*}
And \(c^T = \begin{pmatrix} -2&-5&0&0&0 \end{pmatrix}\). To solve this using nma_simplex
And the solution is
>>>>Current tableau [phase one] 1 0 1 0 0 1 0 0 4
0 1 0 1 0 0 1 0 6
1 1 0 0 1 0 0 1 8
0 0 0 0 0 1 1 1 0
***********************
Current tableau [phase one]
1 0 1 0 0 1 0 0 4
0 1 0 1 0 0 1 0 6
1 1 0 0 1 0 0 1 8
-2 -2 -1 -1 -1 0 0 0 0
pivot row is 1
current basic feasible solution is
4
0
0
0
0
0
6
4
***********************
Current tableau [phase one]
1 0 1 0 0 1 0 0 4
0 1 0 1 0 0 1 0 6
0 1 -1 0 1 -1 0 1 4
0 -2 1 -1 -1 2 0 0 8
pivot row is 3
current basic feasible solution is
4
4
0
0
0
0
2
0
***********************
Current tableau [phase one]
1 0 1 0 0 1 0 0 4
0 0 1 1 -1 1 1 -1 2
0 1 -1 0 1 -1 0 1 4
0 0 -1 -1 1 0 0 2 16
pivot row is 2
current basic feasible solution is
2
6
2
0
0
0
0
0
***********************
Current tableau [phase one]
1 0 0 -1 1 0 -1 1 2
0 0 1 1 -1 1 1 -1 2
0 1 0 1 0 0 1 0 6
0 0 0 0 0 1 1 1 18
***********************
Current tableau [phase two]
1 0 0 -1 1 2
0 0 1 1 -1 2
0 1 0 1 0 6
-2 -5 0 0 0 0
pivot row is 3
current basic feasible solution is
0
6
2
0
0
***********************
Current tableau [phase two]
1 0 0 -1 1 2
0 0 1 1 -1 2
0 1 0 1 0 6
-2 0 0 5 0 30
pivot row is 1
current basic feasible solution is
2
6
2
0
0
***********************
Current tableau [phase two]
1 0 0 -1 1 2
0 0 1 1 -1 2
0 1 0 1 0 6
0 0 0 3 2 34
solution_vector =
2
6
2
0
0
We see that \(x_1=2\) and \(x_2=6\) and \(x_3=2\) with all others zero. We always read the solution from the identity matrix
inside the final tableau. All other \(x_i\) are zero. Hence the optimal solution is \[ \begin{pmatrix} 2\\ 6\\ 2\\ 0\\ 0 \end{pmatrix} \] And the corresponding
optimal objective function is \(2 x_1+5 x_2=2(2)+5(6)=34\)
Result using Matlab’s linprog
>> A=[1,0,1,0,0; 0,1,0,1,0; 1,1,0,0,1];b=[4,6,8];c=[-2,-5,0,0,0];
>> options = optimset('LargeScale','off','Simplex','on');
>> [X,FVAL,EXITFLAG,OUTPUT]=linprog(c,[],[], A,b,zeros(size(c)),[],[],options)
Optimization terminated.
X =
2
6
2
0
0
FVAL =
-34
EXITFLAG =
1
OUTPUT =
iterations: 0
algorithm: 'simplex'
cgiterations: []
message: 'Optimization terminated.'
constrviolation: 0
firstorderopt: 0
3 Appendix
3.1 References
-
Lecture notes, ECE 719 optimal systems, Univ. Wisconsin, Madison spring 2016 by
Professor B. Ross Barmish
-
An introduction to Optimization. Edwin Chong, Stanislaw Zak. Wiley publication,
NY 1996.