Assumptions
- $A \in M_{m\times n}(\Bbb R)$ ($m\le n$) has rank $n$ and basis matrix $B$.
- $x_B$ denotes the basic solution.
- $c_B$ denotes the reduced objective function so that $c^T x=c_B^T x_B$. (So the order/arrangement of basic variables is very important.)
The initial simplex tableau for your problem(, which is already in canoncial form,) will be
\begin{array}{cc|c}
A & I & b \\ \hline
-c^T & 0 & 0
\end{array}
We multiply the $m$ constraints by $B^{-1}$. Note that $\mathbf{e}_1,\dots,\mathbf{e}_n$ can be found in $[B^{-1}A\;B^{-1}]$ since $B$ is chosen from columns of $[A\;I]$. We add a linear combination of the first $m$ row to the last row.
\begin{array}{cc|c}
B^{-1}A & B^{-1} & B^{-1}b \\ \hline
c_B^T B^{-1}A-c^T & c_B^T B^{-1} & c_B^T B^{-1}b \tag{*} \label{t1}
\end{array}
Why is $c_B^T$ chosen? We write $x_B$ as $(x_{B_1},\dots,x_{B_n})^T$. Choose an index $i$ in $\{1,\dots,n\}$, and suppose that $x_j = x_{B_i}$
Case 1: $1\le j \le n$. In the $j$-th column of the last row,
$$(c_B^T B^{-1}A-c^T)_j = c_B^T B^{-1} \mathbf{a}_j -c_j = c_B^T \mathbf{e}_i - c_j = c_{B_i} - c_j = 0.$$
Case 2: $n+1 \le j \le m+n$. Note that $c_{B_i}= 0$. In the $j$-th column of the last row,
$$(c_B^T B^{-1})_{j - n} = c_B^T \mathbf{e}_i = c_{B_i} = 0.$$
Therefore, the entries corresponding to the basic variables in the last row in tableau $\eqref{t1}$ will be zero. Then, $\eqref{t1}$ is a simplex tableau that you can compute from the given optimal solution.
Computational results for your example
Here's the GNU Octave code for finding the optimal tableau. The syntax should be the same in MATLAB except for the solver glpk. You may test the code at Octave Online.
c=[60 100 80]'; A=[144 192 240; 100 150 120; 1 1 1];
b=[120000 60000 500]';
[x_max,z_max]=glpk(c,A,b,[0 0 0]',[],"UUU","CCC",-1)
x=[x_max; b-A*x_max] % Include slack variables
A1=[A eye(3)]; % Part of initial tableau
c1=[c' zeros(1,3)]';
B=A1(:,find(x > 0)')
c_B=c1(find(x > 0));
disp('Initial tableau');
[A1 b; -c1' 0]
disp('Final tableau');
[B^-1*A B^-1 B^-1*b; c_B'*B^-1*A-c' c_B'*B^-1 c_B'*B^-1*b]
The computer output
x_max =
0.00000
400.00000
0.00000
z_max = 4.0000e+04
x =
0.0000e+00
4.0000e+02
0.0000e+00
4.3200e+04
-7.2760e-12
1.0000e+02
B =
192 1 0
150 0 0
1 0 1
Initial tableau
ans =
144 192 240 1 0 0 120000
100 150 120 0 1 0 60000
1 1 1 0 0 1 500
-60 -100 -80 -0 -0 -0 0
Final tableau
ans =
Columns 1 through 6:
6.6667e-01 1.0000e+00 8.0000e-01 0.0000e+00 6.6667e-03 -0.0000e+00
1.6000e+01 0.0000e+00 8.6400e+01 1.0000e+00 -1.2800e+00 0.0000e+00
3.3333e-01 1.1102e-16 2.0000e-01 0.0000e+00 -6.6667e-03 1.0000e+00
6.6667e+00 0.0000e+00 0.0000e+00 0.0000e+00 6.6667e-01 0.0000e+00
Column 7:
4.0000e+02
4.3200e+04
1.0000e+02
4.0000e+04
The optimal tableau:
\begin{equation*}
\begin{array}{ccccccc|l}
& x_1 & x_2 & x_3 & s_1 & s_2 & s_3 & \\ \hline
x_2 & \frac23 & 1 & \frac45 & 0 & \frac{1}{150} & 0 & 400 \\
s_1 & 16 & 0 & \frac{432}{5} & 1 & -\frac{32}{25} & 0 & 43200 \\
s_3 & \frac13 & 0 & \frac15 & 0 & -\frac{1}{150} & 1 & 100 \\ \hline
& \frac{20}{3} & 0 & 0 & 0 & \frac23 & 0 & 40000 \\
\end{array}
\end{equation*}
The proposed answer is not right.
The tableau has a reduced cost for a non-basic variable equal to $0$. This means that there are alternative optimal solutions. As it is not possible to go to another extreme point (the column $y_2$ is $\le 0$), there are optimal solutions which are obtained as:
$x^* + \lambda(-y_2,e_2)$, where $(-y_2,e_2)$ is an extreme direction, so the multiple solutions are:
$\left(\frac 23,0,\frac43\right)' + \lambda(1,1,0)' = \left(\frac23+\lambda, \lambda, \frac43\right)'$ for $\lambda \ge 0$.
So, as an example, you can observe that when $\lambda=1$ you obtain the solution $\left(\frac 53,1,\frac 43\right)$ which is feasible and also optimal.
Best Answer
In this problem, the non-uniqueness in the simplex method comes from the substitution $y = m-n$: a single value of $y$ can be expressed as $m-n$ in many ways.
To see that this is the only reason for non-uniqueness, we can parametrize the solutions found by the simplex method and find all the possible solutions.
The bottom row of your tableau actually corresponds to the equation $z = 55 - 2a - d$. So we know that we obtain the optimal value of $z=55$ exactly when $a=d=0$.
To make this substitution, we delete the $a$ and $d$ columns from the tableau, and get the system of equations $$ \begin{bmatrix} 0 & 1 & -1 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 1 \\ 1 & 0 & 0 & 0 & 0 \end{bmatrix} \begin{bmatrix} x \\ m \\ n \\ b \\ c\end{bmatrix} = \begin{bmatrix}5 \\ 15 \\ 80 \\ 15\end{bmatrix} $$ which tells us that the optimal solutions are those feasible solutions which have $m-n=5$, $b = 15$, $c = 80$, and $x=15$ (and $a=d=0$).
Since $y = m-n=5$ is fixed, the simplex method confirms that actually there's only one solution $(x,y) = (15,5)$ after we undo this substitution and return to the original formulation of the LP.