Simplex: Convert to standard form


Creation date: 2017-08-12

Tags: simplex, python, constraint-programming, linear-programming, optimization

How to solve all different kinds of linear problems?

In the last blog article I explained the Simplex algorithm. A small recap:

We can solve maximization linear problems if every constraint is a \(\leq\) constraint and the right hand side is non negative. For example:

\[ \begin{array}{rrcrcrcr} \text{max} & 45x_1 & + & 60x_2 & + & 50x_3 \\ \text{subject to} & 20x_1 & + & 10x_2 & + & 10x_3 & \leq & 2400 \\ & 12x_1 & + & 28x_2 & + & 16x_3 & \leq & 2400 \\ & 15x_1 & + & 6x_2 & + & 16x_3 & \leq & 2400 \\ & 10x_1 & + & 15x_2 & & & \leq & 2400 \\ & x_1 & & & & & \leq & 100 \\ & & & x_2 & & & \leq & 40 \\ & & & & & x_3 & \leq & 60 \\ \end{array} \]

The general standard form is:

\[ \begin{array}{rrcrl} \text{max} & \sum_{j=1}^n c_jx_j \\ \text{subject to} & \sum_{j=1}^n a_{ij}x_j & \leq & b_j & (i=1,\dots,n) \\ & x_j & \geq & 0 & (j=1,\dots,n) \\ & b_i & \geq & 0 & (i=1,\dots,m) \end{array} \]

We can also express it in matrix form which might be quite nice because we always work with matrices:

\[ \begin{array}{rrcrl} \text{max} & c^Tx \\ \text{subject to} & Ax & \leq & b \\ & x & \geq & 0 \end{array} \] \[ \begin{array}{rcl} c^T & = & (c_1,\dots,c_n) \\ x & = & \begin{pmatrix} x_1 \\ \vdots \\ x_n \end{pmatrix} \\ A & = & \begin{pmatrix} a_{11} & \dots & a_{1n} \\ \vdots & \ddots & \vdots \\ a_{m1} & \dots & a_{mn} \end{pmatrix} \\ b & = & \begin{pmatrix} b_1 \\ \vdots \\ b_m \end{pmatrix} \\ \end{array} \]

Therefore, we have a list of restrictions:

Our goal in this article is to remove all of this restrictions to be able to solve a lot of different linear programs.

Let's use a diet example to explain what we want to achieve. Yes they are called diet problems.

We have a list of ingredients with different nutrients:

\[ \begin{array}{l|c|c|c|c} & \text{Energy (kcal)} & \text{Proteins (g)} & \text{Calcium (mg)} & \text{Price (cent)} \\ \hline Oats & 110 & 4 & 2 & 25 \\ Chicken& 205 & 32 & 12 & 130 \\ Egg & 160 & 13 & 54 & 85 \\ Milk & 160 & 8 & 285 & 70 \\ Cake & 420 & 4 & 22 & 95 \\ Bean & 260 & 14 & 80 & 98 \\ \end{array} \]

Oats might stand for a portion of 25g and the same is true for the other ingredients. From now on we only use "portion" as the unit.

Of course we want to minimize the price for our diet, right? ;) We need to take 2000kcal, 55g proteins and 800mg calcium every day.

This gives us the following linear program in matrix form:

\[ \begin{array}{rl} \text{min} & (25,130,85,70,95,98)x \\ \text{subject to} & \begin{pmatrix} 110 & 205 & 160 & 160 & 420 & 260 \\ 4 & 32 & 13 & 8 & 4 & 14 \\ 2 & 12 & 54 & 285 & 22 & 80 \\ \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \\ x_3 \\ x_4 \\ x_5 \\ x_6 \end{pmatrix} \geq \begin{pmatrix} 2000 \\ 55 \\ 800 \end{pmatrix} \end{array} \]

Here we have a minimization problem with \(\geq\) constraints. Our program can only solve maximization problems at the moment so let's reformulate the problem.

We can simply maximize the negative cost function so instead of

\[ \begin{array}{rl} \text{min} & \begin{pmatrix}25 & 130 & 85 & 70 & 95 & 98\end{pmatrix}x \\ \end{array} \]

we can use

\[ \begin{array}{rl} \text{max} & -\begin{pmatrix}25 & 130 & 85 & 70 & 95 & 98\end{pmatrix}x \\ \end{array} \]

The next restriction was:

Okay we can just use the same idea again. Just multiply everything with \(-1\)...

\[ \begin{array}{rl} \text{max} & -\begin{pmatrix}25 & 130 & 85 & 70 & 95 & 98\end{pmatrix}x \\ \text{subject to} & \begin{pmatrix} -110 & -205 & -160 & -160 & -420 & -260 \\ -4 & -32 & -13 & -8 & -4 & -14 \\ -2 & -12 & -54 & -285 & -22 & -80 \\ \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \\ x_3 \\ x_4 \\ x_5 \\ x_6 \end{pmatrix} \leq \begin{pmatrix} -2000 \\ -55 \\ -800 \end{pmatrix} \end{array} \]

Okay we just moved from restriction #2 to restriction #3:

Actually, why are we not able to use negative values?

We need a basic feasible solution and the basic feasible solution (BFS) needs to have non negative values. The normal way to find a BFS was to add slack variables and assign the values of the right hand side to them. This isn't possible with a negative right hand side.

First of all we add another variable called \(x_0\) (new column) which can be set high enough to fulfill all the constraints. Every other value is set to 0. To save some space in our matrix we use the same notation as last time.

Instead of the one above we can write our tableau:

\[ \begin{array}{rrrrrr|r} -110 & -205 & -160 & -160 & -420 & -260 & -2000\\ -4 & -32 & -13 & -8 & -4 & -14 & -55\\ -2 & -12 & -54 & -285 & -22 & -80 & -800\\ \hline 25 & 130 & 85 & 70 & 95 & 98 & 0 \end{array} \]

Let's add \(x_0\) as our last column:

\[ \begin{array}{rrrrrrr|r} -110 & -205 & -160 & -160 & -420 & -260 & -1 & -2000\\ -4 & -32 & -13 & -8 & -4 & -14 & -1 & -55\\ -2 & -12 & -54 & -285 & -22 & -80 & -1 & -800\\ \hline 25 & 130 & 85 & 70 & 95 & 98 & 0 & 0 \end{array} \]

In this example we would be able to set \(x_1 = \dots = x_6 = 0\) and \(x_0 = 2000\) to fulfill all of our constraints.

But we don't really want to solve this problem because \(x_0\) isn't in our objective function and only \(x_0\) would be set. That doesn't make any sense.

Instead we want no \(x_0\) at all. How can we achieve that without just going back to the step before? At the moment we are able to find a BFS but for a wrong problem. In the step before we couldn't find one but had our correct LP.

Well, we can simply try to get rid of \(x_0\) by solving the tableau with a different objective. We want to remove \(x_0\) therefore we can try to minimze \(x_0\) or because minimization isn't our standard form: maximize \(-x_0\).

The tableau will look like this now:

\[ \begin{array}{rrrrrrr|r} -110 & -205 & -160 & -160 & -420 & -260 & -1 & -2000\\ -4 & -32 & -13 & -8 & -4 & -14 & -1 & -55\\ -2 & -12 & -54 & -285 & -22 & -80 & -1 & -800\\ \hline 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 \end{array} \]

Let's add our slack variables for this problem:

\[ \begin{array}{rrrrrrrrrr|r} -110 & -205 & -160 & -160 & -420 & -260 & -1 & 1 & 0 & 0 & -2000\\ -4 & -32 & -13 & -8 & -4 & -14 & -1 & 0 & 1 & 0 & -55\\ -2 & -12 & -54 & -285 & -22 & -80 & -1 & 0 & 0 & 1 & -800\\ \hline 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \end{array} \]

Normally we looked at negative values in our objective row (bottom row) to start our optimization step. Here it is a bit different. Our first step is to get only non negative b values. To achieve that we use Gaussian elimination. We want to have one 1 in the \(x_0\) column and the rest values for that column should be 0. Because our goal is to have non negative values for the whole right hand side we choose to get a 1 in the first row where the right hand side has it's smallest value.

We do the same thing as normal like having an entering and a leaving variable.

After doing this step we get the following tableau:

\[ \begin{array}{rrrrrrrrrr|r} 110 & 205 & 160 & 160 & 420 & 260 & 1 & -1 & -0 & -0 & 2000 \\ 106 & 173 & 147 & 152 & 416 & 246 & 0 & -1 & 1 & 0 & 1945 \\ 108 & 193 & 106 & -125 & 398 & 180 & 0 & -1 & 0 & 1 & 1200 \\ \hline -110 & -205 & -160 & -160 & -420 & -260 & 0 & 1 & 0 & 0 & -2000 \\ \end{array} \]

Let's have a short look on our new tableau. Our goal was to maximize \(-x_0\) and at the moment we have the objective \(-2000\). Which also means that we have \(x_0 = 2000\). We can validate this by having a look on our first row. Now we have a tableau which we can solve using our normal simplex algorithm because we have only positive b values. We have negative values in our objective row which means we are able to improve our solution and as described we have our BFS.

If we run the simplex algorithm on this tableau (without adding extra slack variables because we have a BFS already) then we should be able to get an objective of 0 which would mean that we can ignore \(x_0\) and have a BFS for our real problem or we aren't able to get an objective of 0. Then we know that we can't solve the problem.

\[ \begin{array}{rrrrrrrrrr|r} 0.1 & 1& 0.35 & 0& 0& 0.33 & 0.03 & 0& -0.03 & 0& 0.66 \\ -0.01 & 0& 0.16 & 1& 0& 0.24 & 0 & 0& 0 & -0& 2.51 \\ 0.22 & 0& 0.15 & 0& 1& 0.37 & -0.01 & -0& 0.02 & 0& 3.49 \\ \hline 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \end{array} \]

First of all we have our objective set to 0 which is exactly what we aimed for. For better visibility we will write this tableau in a different form to understand how we can use it to achieve our real goal to get a good diet. We first add a row displaying the name of the variables in our tableau:

\[ \begin{array}{rrrrrrrrrr|r} x_1 &x_2 & x_3 & x_4 & x_5 & x_6 & x_0 & s_1 & s_2 & s_3 & b \\ \hline 0.1 & 1& 0.35 & 0& 0& 0.33 & 0.03 & 0& -0.03 & 0& 0.66 \\ -0.01 & 0& 0.16 & 1& 0& 0.24 & 0 & 0& 0 & -0& 2.51 \\ 0.22 & 0& 0.15 & 0& 1& 0.37 & -0.01 & -0& 0.02 & 0& 3.49 \\ \hline 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \end{array} \]

Therefore our first constraint can be written as:

\[ 0.1x_1 + 1x_2 + 0.35x_3 + 0x_4 + 0x_5 + 0.33x_6 + 0.03x_0 + 0s_1 - 0.03s_2 + 0s_3 = 0.66 \]

Let's remove the zero entries to make it shorter.

\[ 0.1x_1 + x_2 + 0.35x_3 + 0.33x_6 + 0.03x_0 - 0.03s_2 = 0.66 \]

The variable which is in our basis and isn't zero is \(x_2\) so let's rewrite it to have \(x_2\) on on side and the rest on the other side.

\[ x_2 = 0.66 - 0.1x_1 - 0.35x_3 - 0.33x_6 - 0.03x_0 + 0.03s_2 \]

If we do the same thing for the other constraints we have this at the end:

\[ \begin{array}{rcl} x_2 & = & 0.66 - 0.1x_1 - 0.35x_3 - 0.33x_6 - 0.03x_0 + 0.03s_2 \\ x_4 & = & 2.51 + 0.01x_1 - 0.16x_3 - 0.24x_6 \\ x_5 & = & 3.49 - 0.22x_1 - 0.15x_3 - 0.37x_6 + 0.01x_0 - 0.02s_2 \\ \end{array} \]

Before we really build our new objective function we should remove the \(x_0\) values because we don't want to have that variable anymore. It's not really part of our problem, it just helped us to get to this stage but now bye bye \(x_0\).

\[ \begin{array}{rcl} x_2 & = & 0.66 - 0.1x_1 - 0.35x_3 - 0.33x_6 + 0.03s_2 \\ x_4 & = & 2.51 + 0.01x_1 - 0.16x_3 - 0.24x_6 \\ x_5 & = & 3.49 - 0.22x_1 - 0.15x_3 - 0.37x_6 - 0.02s_2 \\ \end{array} \]

Now let's have a look back to our objective function from the beginning...

\[ 25x_1+130x_2+85x_3+70x_4+95x_5+98x_6 \]

To get our new objective function we simply need to replace the terms \(x_2,x_4 \text{ and } x_5\) with the new expression of that terms mentioned above.

This gives us:

\[ -7.88x_1 + 14.29x_3 + 3.35x_6 + 0.2s_1 + 2.78s_2 + 0.06s_3 \]

We also get our objective value \(-592.07\). Our simplex algorithm can maximize this value and at the end we have to take the absolute value to get our solution.

If we run our simplex algorith on this tableau we get:

\[ \begin{array}{rrrrrrrrrr|r} x_1 &x_2 & x_3 & x_4 & x_5 & x_6 & s_1 & s_2 & s_3 & b \\ \hline 1 & 9.9 & 3.4 & 0 & 0 & 3.3 & 0 & -0.3 & 0 & 6.47 \\ -0 & 0.1 & 0.2 & 1 & 0 & 0.3 & 0 & -0 & -0 & 2.60 \\ 0 & -2.2 & -0.6 & 0 & 1 & -0.3 & -0 & 0.1 & -0 & 2.08 \\ \hline 0 & 77.8 & 41.3 & 0 & 0 & 29.2 & 0.2 & 0.2 & 0.1 & -541.1 \end{array} \]

We read the tableau as always but use the absolute value of our obj. Therefore our diet plan is to eat 6.47 portion of oats, 2.6 poritions of milk and a bit more than 2 cake portions. You might think that this isn't the best option because well portions of cake and there is so much stuff on the list we don't eat at all... We can change that by adding extra constraints. I personally like oats but too many of them is just not tasty.

Therefore I suppose we can add the following constraints:

# oats
m.add_constraint(x1 <= 4)
# chicken
m.add_constraint(x2 <= 3)
# egg
m.add_constraint(x3 <= 2)
# milk
m.add_constraint(x4 <= 8)
# cake
m.add_constraint(x5 <= 1)
# bean
m.add_constraint(x6 <= 2)

It's pretty simple to read: We don't want to eat more than 4 portions of oats and not more than one cake portion and so on. This gives us our new diet:

oat 4
milk 3.875
cake 1
bean 2

That looks much better to me. Now you might say it doesn't make much sense to use 3.875 milk portions and you want to have only whole numbers. The simplex algorithm generates real number all the time and in some cases you might want to have only integers or part of the values should be integers. I'll write another blog article about that kind of stuff.

Hope you enjoyed this one!

I changed the structure of my repository and added the new functionality but it's nothing special so I don't explain the changes in this blog post.

You can check out the new version on Github.

Next post: The Dual problem.

If you're interested in another series check out my posts on constraint programming: Solving Sudokus again but also much more, faster and in depth: Constraint Programming Solver from Scratch.

If you've enjoyed this and some other posts on my blog I would appreciate a small donation either via PayPal or Patreon whereas the latter is a monthly donation you can choose PayPal for a one time thing. For the monthly donation you'll be able to read a new post earlier than everyone else and it's easier to stay in contact.



Want to be updated? Consider subscribing on Patreon for free
Subscribe to RSS