After finishing with the Discrete Optimization course I started to get to know the Simplex algorithm for solving linear programs. And if you learn an algorithm the best is if you try it out and do some coding (well implement the algorithm itself).
“Und I did.” 
OK, I just started to implement my version of the Simplex algorithm in Python. Well, to be honest, the algorithm itself isn’t complicated. You just have to do the pivoting over and over again — or as long as you do not get a final dictionary (or decide that you run in circles because of degeneracy). As every development, I do this in more phases too because there’s no good idea behind writing a big bang software (except of Racket where with big-bang you can create a GUI 🙂 ).
The structure of the input data for the algorithm has to be defined. What do you need for the Pivoting steps of the Simplex algorithm? Yes, a (feasible) dictionary. The best is if it is feasible because in this case you do not have any initialization phase. And what does a dictionary contain? Yes, the A matrix with the basic and non-basic variables, the b vector with the constraints and the expression to maximize.
The input will be a file. I find defining such a structure in a file is easier than any other kind (for example XML) of data-providing. And as for the first version I take a feasible dictionary, so slack-variables are included in the dictionary — AND I assume no initialization is required. Well, now I have a strange urge to try out what happens if I provide an non-feasible dictionary as input…
The structure looks as follows:
<br /> #[Line 1] m n [m: number of basic variables, n: number of non-basic variables]<br /> #[Line 2] B1 B2 ... Bm [the list of basic indices m integers]<br /> #[Line 3] N1 N2 ... Nn [the list of non-basic indices n integers]<br /> #[Line 4] b1 .. bm (m floating point numbers)<br /> #[Line 5] a11 ... a1n (first row of A matrix)<br /> #....<br /> #[Line m+4] am1 ... amn (mth row of A matrix)<br /> #[Line m+5] z0 c1 .. cn (objective coefficients (n+1 floating point numbers))<br />
Because I’m an OO developer since 10 years (and I’m doing Java since 9 years) my Python code is structured as same as Java-code and I use closures and co. currently very rarely.
However, the first step was to create a printer function which displays the dictionary. Every line of the dictionary is represented as a string line in the format:
<br /> x_i = b_i + a_ij (i in 1..m; j in 1..n)<br /> _______________________________________<br /> z = z0 + cj*mj (j in 1..n)<br />
And the line before the objective function and objective value… is hard-coded a fixed length line at the beginning but I have the plan to make it dynamic too. Anyhow, now the whole output is horrible if you look at it on the console with a big set of non-basic variables.
After this I started to write the parsing of the input-dictionary described with the format above. This was not bad because the structure of the file tells clearly what to do and how to import. However, the first version of the code is not fault-tolerant and does not do any validation on the file. So if you provide bad values for m or n you can get Exceptions or a no-good dictionary as input to the algorithm.
Third step was to implement Bland’s Rule: select the entering variable and the leaving variable. If there are more possible choices for entering variables choose the one with the lower index. If there are more possible choices for leaving variables, choose the one with the lowest index. This is no problem if you know how to select entering and leaving variables.
Now we have entering and leaving variables. Eventually it can happen that we do not have an entering variable — in this case the dictionary is final, so we do not have a leaving variable neither. If we do not obtain a leaving variable than the dictionary (and the problem too) is unbounded so the optimal value is infinity.
So, let there be an entering and a leaving variable. Now it is time to do the Pivoting step. I cut this operation into 5 (+1) steps as follows:
<br /> # 0. secure the value of the entering variable (and the objective value)<br /> # 1. swap places: entering and leaving<br /> # 2. divide new line with value of entering variable<br /> # 3. rearrange other lines of A_matrix (except the leaving): replace the entering variable with the new line (of the leaving)<br /> # 4. recalculate b (except the leaving): add the b value from the new line<br /> # 5. recalculate z: replace entering variable with the leaving one<br />
I think the steps speak for themselves but I will write some words to clarify what I do and why I do this.
0. It is easier to secure those two variables so we can update our dictionary without any regrets and can multiply without using references on the A-matrix or the b vector.
1. This is a bit fishy I guess — how do we come to let those two simply swap places? What is with coefficient of the entering variable in the line where it enters? Well, yes, but we do have two lists for the basic and non-basic variables and these lists contain the entering variable and the leaving variable — or at least the index of them. So it is not a big deal to swap them: the entering variable is entering the round of basic variables, the leaving variable leaves to get a non-basic one. And do not forget to set the value (coefficient) of the leaving variable to ‘-1.0’ (because in the step of the entering we have to subtract it from the “left side”).
2. The new line is the row where the entering variable entered. I should rename my comments 🙂 Anyhow we have to divide the values of the coefficients of this row with the coefficient of the entering variable to get it a single one on the left hand side of the equation. If we want to be mathematically more correct we have to divide both sides with this value — but for this we had had to move the coefficient with the entering variable to the left-hand side but for the basic variables we do not have a list because they have to be always multiplied by ‘1.0’. And because of the divisions we do need floating-point numbers in our A matrix and the b vector.
3. Line means in this case row too. This is the algebra part of the algorithm: replace in every line (except the one where the new variable entered) the entering variable with its new value (yes, it is the righthand side of the equation of the old row of the A matrix) and rearrange the row. So you have to multiply the coefficients of the replaced values with the coefficient of the old variable and sum the already existing ones up with the newcomers.
4. Can be done with step three under the same hut but you can change it too. If you look at my implementation I hid this step behind the calculations of 3. And as in the previous step: don’t forget to multiply with the coefficient of the entering variable for each line.
5. This is the same step as 3 and 4 but you do the calculations on the objective row. And if you do it right you get the same objective value you “secured” in the 0th step. If you are lazy (like me) you can do the algebra-part withouth the b vector and then set the value of z to the previously calculated objective value.
Now that we can find an entering and a leaving variable and do the pivoting we only need to implement the logic of the optimization algorithm:
- Test the dictionary on finality. If it is final, then exit with the objective value.
- Get the entering variable.
- Get the leaving variable based on the entering variable.
- If the leaving variable does not exist, then exit with the remark that the problem is unbounded.
- Pivot the dictionary.
- Continue with step 1.
As of degeneracy there we could add a counter that limits the algorithm to some number of pivoting steps when stalling occurs. Because of Bland’s Rule there will be no cycling so we do not have to keep the dictionaries to remember if we visited them or not.
As mentioned I plan to add the initialization phase to the algorithm. For this I plan to write a post too similar to this one. After that part is done we (you and I) can use this code for solving some discrete optimization problems. Naturally there will be validation on the input file to avoid exceptions or other misbehavior.
As always, the code is available at my GitHub repo.
 Sorry but I could not stand to quote this from The Rocky Horror Picture Show