A brief introduction to Linear Programming

• What is Linear Programming?
• Solve LP with Graphing
• Standard Form of LP
• Simplex Method
• Simplex Algorithm with Python from Scratch

What is Linear Programming?

Let’s find your answer to the following problem. Suppose you are making two kinds of jam (call jam#1 and jam#2) using strawberry and grape. You need 1packages of grape and 2package of strawberry for 1kg of jam#1 and 3 packages of grape and a package of strawberry for 1kg of jam#2. If you have 7packages of grape and 4packages of strawberry, and with each kg of jam you made, you earn 1 point. How many kgs of jam#1 jam#2 should you make to earn the maximum point? Other materials and the remaining fruit packages are not necessary to be considered.

This kind of problem, maximum or minimum a cost under a specific restriction, is not only commonly found in daily life but also in economy, science, and engineering.

Now, let x and y be the weight of jam#1 and jam#2 in kg. Your problem can be rewritten as follow. Here, max. refers to maximize and s.t. means subject to. This kind of problem is called Linear Programming. In the following part, we call max. part “objective function” and s.t. part “constrain” or “restriction” or “requirement”.

A Linear Programming (LP) is a type of mathematical optimization, a selection of variables that make a function takes maximum or minimum value, which to be optimized function and all requirements are represented by a linear relationship.

Solve LP with Graphing

We can solve this problem with simple arithmetic learned in elementary school. In the following graph, the blue lines are the boundary of inequality expression of the restriction. The green part is the Feasible Region where all the requirements are met. To find the optimal answer to the problem, we plot some of the destination function contours with red lines. Among these contours, those which intercept with the feasible region at the highest value is the optimal answer.

From the graph, your point, x+y, is getting the largest value while all the requirements are satisfied when x=1 and y=2. Means, you can earn a maximum point by making 1kg of jam#1 and 2 kg of jam#2. Alright, the problem is solved.

Now, let’s say you want to make another kind of jam (call jam#3) which you need 1package of grape and strawberry to make 1kg of it. To deal with this in maths, here, you need another variable, z, for jam#3 and your LP is:

Are you thinking about graphing these inequality restrictions? If you are going to solve it with a method we did just before this, now you may need a three-dimension graph. The following 3D graphic is an example. Can you find the fittest amount of each kind of jam to be made which meets all the requirement while maximizing your point?

Do not worry if you cannot get any answer to the question. In general, it is difficult for a human to understand data in high-dimension space. However, we cannot give up on solving our problem here. In this post, I will introduce a well-known method in solving the LP, Simplex Method in the following parts.

Standard Form of LP

Before discussing on Simplex Method, let’s have a look at a standard form of linear programming.

Some characteristics of a standard form of LP are:
1. It is about maximization, not minimization.
2. Each variable is non-negative.
3. Other constraints are in a linear combination of the variables with an equation(≤).

Of course, it does exist LP which is asked to minimize function or containing the opposite inequality signs in requirement expressions. In that case, some techniques are used to transform those LP into standard forms. Here I am introducing some examples.

Simplex Method

I will not dive into deep maths on this topic, but I will share the method’s concept and calculation examples.

Idea

The method name is derived from a geometric term Simplex. Actually, simplex is not actually used in this algorithm, but simplicial cones are where the idea works. The feasible region, saying in a geometric term, is a convex polytope. The simplicial cones are corners of the polytope. Let’s use the earlier jam-making example again here.

In this case, red points are the vertex of the convex polytope made by the constrains. These extreme points are called basic feasible solutions.

In the simplex method, if an LP in standard form has a maximum objective value, we can reach the optimal point by continue going along edges of the convex polytope to extreme points with bigger and bigger objective values. The LP can be concluded that it has no solution when an unbounded edge is visited.

How does it work?

Simplex method (or simplex algorithm) takes two steps to define the solution of an LP. In step1, a starting extreme point is found, and a basic feasible solution is found or feasible region is empty as a result. In step2, a simplex method is applied using a basic feasible solution from PhaseI as a starting point. The answer to the LP is said found when an optimum basic feasible solution is reached. Otherwise, the objective function is unbounded above.

Let’s do some maths!

Here, we solve the remaining problem of jam making example with Simplex Method. First, we make an initial dictionary (tableaus) of this LP as the following format where a, b are new variables (call slack variables) to make inequality constraints transform to equality. To make it simple, we start a basic feasible solution by setting all variables on the right-hand side to 0. Now we got an initial basic feasible solution (x,y,z,a,b) = (0,0,0,7,10) an objective function value of 0.

Focus on objective function p = x+y+z, the value of p is increasing when moving x,y,z from 0 to greater values. Taking the simplest way, we fix y and z at 0 and increase x value to make p bigger. Since a=7-x and b=10–2x are non-negative, therefore, x must not exceed 5. Hence, we increase x from 0 to 5, and at that time, b became 0. Our dictionary is updated by exchanging b and x in the second equality and use a new express of x to update the first equation. The process of moving basic and non-basic variables in the dictionary is called a pivot operation. The new dictionary is shown in the following figure. Ending this process, we are now moving extreme point explain in the previous section from (x,y,z)=(0,0,0) to (x,y,z)=(5,0,0). Now the objective function value is p=x+y+z=5.

Similarly, if we keep increase y,z, then the value of p keeps increasing. Fixing the value of z at 0 and increase y value. a=2-(5/2)y and x=5-(1/2)y are non-negative, therefore, y is less than or equal to 4/5. We move y from 0 to 4/5 where a becomes 0. We update the dictionary by exchanging a and y of the first equality. The updated dictionary is shown in the following figure. In this step, we are moving extreme point from (x,y,z)=(5,0,0) to (x,y,z)=(23/5, 4/5, 0) and the objective function value is increase to 27/5.

With the same procedure, we increase z to make p greater. y=4/5 — (1/5)z and x=23/5 — (2/5)z are non-negative, so the maximum value of z we can move to is 4 and at that time y becomes 0. The proceeded dictionary is as follows. In this step we are moving extreme point from (x,y,z)=(23/5, 4/5, 0) to (x,y,z)=(4,0,3).

Now, p=7–2y-a. Since a and y are non-negative, there is no variable can be move from zero to make p greater. Therefore, (x,y,z)=(3,0,4) is the optimum basic feasible solution where the objective function value is 7, and the problem is solved.

How was it? Using the simplex method, you do not need to do a three-dimension graphing as shown earlier. Before starting your work with the simplex method make sure your LP is in a standard form and your initial dictionary is feasible. In case your initial dictionary is not feasible a technique of using an artificial variable is popular. Anyway, we will not discuss this method here.

In the following section, I am talking about how to code the simplex algorithm to solve an LP with python. If you have no idea about programming or you prefer using open-source solver, here is the end of this post.

Simplex Algorithm with Python from Scratch

There are some open-source solvers where you can easily access to solve a linear programming problem. However, in this section, I will share an idea on how to code the algorithm from scratch and also to see LP in another form.

Mathematical Preparation

To make it simple in coding with python, I will transform LP into a matrix and vector expression.

Let matrix A be coefficients of equations in the constrain, vector b be the right-hand side of constraining equality, and vector c be the coefficient of all variables in the subjective function. vector x_N and x_B are non-basic and basic variables vectors respectively. Now we decompose matrix A and coefficient vector c into a non-basic and a basic part as shown in the above figure.

It is also can be rewritten as the following form to make a dictionary (tableaus) of our LP. Here, matrix B and N, vector c_N, c_B are new decomposed matrixes, vectors after a pivot operation.

For example, when process pivot operation on x1 and x5, our B, N,c_B,c_N are now decomposed to:

In this form, we can easily select a non-basic variable for pivot operation just by defining a positive element of the coefficient vector of x_N in the new dictionary. Also, to stop the repetition, the following contains are considered:

Sample codes

Select variable for pivot operation

Non-basic variables

def selectNindex(x,n_variable,x_N):
  pos_exist = 0
  for i in range(n_variable):
    if x[i]>0:
      pos_exist = 1
      return i,x_N[i]
    else:
      continue
  return -1

Basic variables

def selectBindex(x,y,n_constrain,x_B):
  res = x/y
  print(res)
  min_index = 0
  for i in range(n_constrain):
    if (res[i]<res[min_index]):
      min_index = i
    else:
      continue
return min_index,x_B[min_index]

Decompose coefficient vectors and matrixs

coefficient vector

def updateC(b,n,n_constrain,n_variable,c):
  cTB = np.zeros(n_constrain)
  cTN = np.zeros(n_variable)
  for idx,i in enumerate(b-1):
    cTB[idx] = c[i]
  for idx,i in enumerate(n-1):
    cTN[idx] = c[i]
  return cTB,cTN

coefficient matrix

def updateBN(b,n,n_constrain,n_variable,A):
  B = np.zeros((n_constrain,n_constrain))
  N = np.zeros((n_variable,n_constrain))
  AT = np.transpose(A)
  
  for idx,i in enumerate(b-1):
    B[idx] = AT[i]
  for idx,i in enumerate(n-1):
    N[idx] = AT[i]
  return np.transpose(B),np.transpose(N)

Calculate elements of new dictionary

def nextStep(B,N,b,cB,cN,n_constrain):
  Binv = np.linalg.inv(B)
  nb = np.dot(Binv,b)
  nN = np.dot(Binv,N)
  nB = np.identity(n_constrain)
  ncN = (cN-np.dot(np.dot(cB,Binv),N))
  Z = np.dot(cB,np.dot(Binv,b))
  return nB,nN,nb,ncN,Z,(nb>=0)

Check optimal state

def checkMoving(c_N,n_variable):
  pos_var = 0
  for i in range(n_variable):
    if c_N[i]>0:
      pos_var += 1
    else:
      pos_var += 0
  if pos_var > 0:
    return True
  else:
    return False

solve LP:

def solve(n_variable,n_constrain,x_B,x_N,b,c,A):
  x = True
  i = 0
  c_B,c_N = updateC(x_B,x_N,n_constrain,n_variable,c)
  B,N = updateBN(x_B,x_N,n_constrain,n_variable,A)
  idx_n,x_n_idx = selectNindex(c_N,n_variable,x_N)
  while x and (i<10):
    i += 1
    idx_b,x_b_idx = selectBindex(b,N[:,idx_n],n_constrain,x_B)
    #print("choose non-base :#",idx_n, " at x",x_n_idx)
    #print("choose base     :#",idx_b, " at x",x_b_idx)
    print("Pivote: x",x_n_idx," & x",x_b_idx)
    x_B[idx_b] = x_n_idx
    x_N[idx_n] = x_b_idx
    print("x_B : ",x_B)
    print("x_N : ",x_N)
    #print("updating C ...")
    c_B,c_N = updateC(x_B,x_N,n_constrain,n_variable,c)
    #print("c_B:",c_B)
    #print("c_N:",c_N)
    #print("finished!")
    #print("updating B,N ...")
    B,N = updateBN(x_B,x_N,n_constrain,n_variable,A)
    #print("B:")
    #print(B)
    #print("N:")
    #print(N)
    #print("finished!")
    #print("new tableaus ...")
    nB,nN,nb,nc_N,Z,feasible=nextStep(B,N,b,c_B,c_N,n_constrain)
    print("#",i,"##################")
    print(feasible)
    if feasible.all():
      print("Feasible")
    else:
      print("infeasible")
      break
    print("b    : ",nb)
    print("N    : ")
    print(nN)
    print("Z    : ",Z)
    print("c_N  : ",nc_N)
    print("########################")
    x = checkMoving(nc_N,n_variable)
    if x:
      idx_n,x_n_idx = selectNindex(nc_N,n_variable,x_N)
      print("Next>>")