In [2]:
%qtconsole
%load_ext autoreload
%autoreload 2
In [4]:
import numpy as np

The simplex algorithm

Implementing the algorithm seen during the lectures

In [1]:
def simplexit(A,b,c, init, debug=True, expectNullOptimum=False):
    """
    Implements the standard matrix form of the simplex algorithm

    Hugues Talbot 2018, from an initial version in R ca. 2009.
    """
    iter = 0
    dim=A.shape[0]
    nbvar=A.shape[1]
    lastx=0
    lastxl=0
    sol=None
    IBVmem=None # to detect cycles

    vars=np.arange(0,nbvar) # indices start at 0
    IBV=init # Initial basis vector
    while (True):
        if debug:
            print("*** Iteration: ", iter)
        NBV = np.setdiff1d(vars,IBV) # Set complement -> non basis variables
        if debug:
            print("    IBV = ", IBV)
            print("    NBV = ", NBV)

        B=A[:,IBV]
        E=A[:,NBV]
        cb = c[IBV]
        ce = c[NBV]

        if debug:
            print("    B = \n", B)
            print("    E = \n", E)

        try:
            iB = np.linalg.inv(B)
        except:
            print("B is not invertible, error")
            break

        if debug:
            print("    iB = \n", iB)

        bbar = iB.dot(b)
        if (np.min(bbar) < -1e-7):
            print("bbar is not positive, not a feasible solution")
            break

        if debug:
            print("    bbar= ", bbar)

        cebart = ce.T - cb.T.dot(iB).dot(E)

        if debug:
            print("  reduced costs", cebart)

        if (np.min(cebart) >= 0):
            print(" optimum found!")
            sol = (IBV, bbar)
            print("IBV = ", IBV)
            print("Bbar= ", bbar)
            if ((np.min(cebart) > 0) or expectNullOptimum):
                break
            if (IBVmem is None):
                IBVmem = IBV
            else:
                if (IBV == IBVmem):
                    print("Cycle detected\n")
                    break

        l = np.argmin(cebart) # may not be unique
        P = iB.dot(A[:,NBV[l]])
        if (debug):
            print("P = ", P)
        if (np.max(P) <= 0):
            print("Solution in unbounded!")
            break
        ratios = bbar/P ## elementwise division
        if (debug):
            print("Ratios = ", ratios)
        ## Ignore negative members f P
        negP = P < 0
        ratios[negP] = np.inf ## this way they will not be selected as the min
        nanP = np.isnan(ratios)
        ratios[nanP] = np.inf ## nan are eliminated as well
        if (debug):
            print("Fixed ratios = ", ratios)
        xl = np.min(ratios)
        j = np.argmin(ratios)
        lastx = xl
        lastxl = NBV[l]
        if (debug):
            print("Variable (", IBV[j],") in the basis replaced by variable (",NBV[l],")\n")
            print("xl = ",xl)

        IBV[j] = NBV[l]
        if (debug):
            print("IBV = ", IBV)

        iter += 1
        ## end while
        
    return sol
In [2]:
## exercice 2, TD1 (production de A,B,C liƩe)
def exampleII(debug=True):
    # A <- matrix(data =c(-1, 0, 1, 2, -1, 2, 0, 1, 3, 1, 0, 0, 0, 1, 0, 0, 0, 1), ncol=6, nrow=3)
    A = np.array([[-1,  2, 0, 1, 0, 0], 
                  [ 0, -1, 1, 0, 1, 0], 
                  [ 1,  2, 3, 0, 0, 1]])
    b = np.array([  0,  0, 35])
    C = np.array([-10,-36,-44, 0, 0, 0])
    IBV = np.array([3,4,5]) # index of the IBV
    sol = simplexit(A,b,C,IBV,debug)
    
    # calcul optimum
    z = C[sol[0]].dot(sol[1])
    return(z)
In [6]:
exampleII(debug=True)
*** Iteration:  0
    IBV =  [3 4 5]
    NBV =  [0 1 2]
    B = 
 [[1 0 0]
 [0 1 0]
 [0 0 1]]
    E = 
 [[-1  2  0]
 [ 0 -1  1]
 [ 1  2  3]]
    iB = 
 [[1. 0. 0.]
 [0. 1. 0.]
 [0. 0. 1.]]
    bbar=  [ 0.  0. 35.]
  reduced costs [-10. -36. -44.]
P =  [0. 1. 3.]
Ratios =  [        nan  0.         11.66666667]
Fixed ratios =  [        inf  0.         11.66666667]
Variable ( 4 ) in the basis replaced by variable ( 2 )

xl =  0.0
IBV =  [3 2 5]
*** Iteration:  1
    IBV =  [3 2 5]
    NBV =  [0 1 4]
    B = 
 [[1 0 0]
 [0 1 0]
 [0 3 1]]
    E = 
 [[-1  2  0]
 [ 0 -1  1]
 [ 1  2  0]]
    iB = 
 [[ 1.  0.  0.]
 [ 0.  1.  0.]
 [-0. -3.  1.]]
    bbar=  [ 0.  0. 35.]
  reduced costs [-10. -80.  44.]
P =  [ 2. -1.  5.]
Ratios =  [ 0. -0.  7.]
Fixed ratios =  [ 0. inf  7.]
Variable ( 3 ) in the basis replaced by variable ( 1 )

xl =  0.0
IBV =  [1 2 5]
*** Iteration:  2
    IBV =  [1 2 5]
    NBV =  [0 3 4]
    B = 
 [[ 2  0  0]
 [-1  1  0]
 [ 2  3  1]]
    E = 
 [[-1  1  0]
 [ 0  0  1]
 [ 1  0  0]]
    iB = 
 [[ 0.5  0.   0. ]
 [ 0.5  1.   0. ]
 [-2.5 -3.   1. ]]
    bbar=  [ 0.  0. 35.]
  reduced costs [-50.  40.  44.]
P =  [-0.5 -0.5  3.5]
Ratios =  [-0. -0. 10.]
Fixed ratios =  [inf inf 10.]
Variable ( 5 ) in the basis replaced by variable ( 0 )

xl =  10.0
IBV =  [1 2 0]
*** Iteration:  3
    IBV =  [1 2 0]
    NBV =  [3 4 5]
    B = 
 [[ 2  0 -1]
 [-1  1  0]
 [ 2  3  1]]
    E = 
 [[1 0 0]
 [0 1 0]
 [0 0 1]]
    iB = 
 [[ 0.14285714 -0.42857143  0.14285714]
 [ 0.14285714  0.57142857  0.14285714]
 [-0.71428571 -0.85714286  0.28571429]]
    bbar=  [ 5.  5. 10.]
  reduced costs [ 4.28571429  1.14285714 14.28571429]
 optimum found!
IBV =  [1 2 0]
Bbar=  [ 5.  5. 10.]
/usr/local/Cellar/ipython/7.3.0/libexec/vendor/lib/python3.7/site-packages/ipykernel_launcher.py:77: RuntimeWarning: invalid value encountered in true_divide
Out[6]:
-500.0
In [7]:
## Bakery problem
def bakery2(debug=True):
    # A <- matrix(data =c(-1, 0, 1, 2, -1, 2, 0, 1, 3, 1, 0, 0, 0, 1, 0, 0, 0, 1), ncol=6, nrow=3)
    A = np.array([[ 5, -1, 1, 1, 0], 
                  [ 1,  0, 0, 0, 1]])
    b = np.array([30, 5])
    C = np.array([-30, 4, -4, 0, 0])
    IBV = np.array([3,4]) # index of the IBV
    sol = simplexit(A,b,C,IBV,debug)
    
    # calcul optimum
    z = C[sol[0]].dot(sol[1])
    return(z)
In [8]:
bakery2(debug=False)
 optimum found!
IBV =  [2 0]
Bbar=  [5. 5.]
Solution in unbounded!
/Users/talboth/Library/Python/3.6/lib/python/site-packages/ipykernel_launcher.py:77: RuntimeWarning: divide by zero encountered in true_divide
Out[8]:
-170.0

Alloy problem

In [9]:
## Alloy problem
def alloys(debug=True):
    A = np.array([[ 10, 10, 40, 60, 30, 30, 30, 50, 20, 1, 0, 0], 
                  [ 10, 30, 50, 30, 30, 40, 20, 40, 30, 0, 1, 0],
                  [ 80, 60, 10, 10, 40, 30, 50, 10, 50, 0, 0, 1]])
    b = np.array([30, 30, 40])
    C = 100*np.array([4.1, 4.3, 5.8, 6.0, 7.6, 7.5, 7.3, 6.9, 7.3, 100, 100, 100])
    IBV = np.array([9,10,11]) # index of the IBV
    sol = simplexit(A,b,C,IBV,debug)
    
    # calcul optimum
    z = C[sol[0]].dot(sol[1])
    return(z)
In [10]:
alloys(debug=False)
 optimum found!
IBV =  [3 1 0]
Bbar=  [4.00000000e-01 6.00000000e-01 1.11022302e-16]
Out[10]:
498.0

Collisions

In the dual, if the solution is unbounded, this means that there is no solution in the primal

In [11]:
## Alloy problem
def collisions(debug=True):
    A = np.array([[ 1,  1, -2, -1, -1, 2, 1, 0], 
                  [ 1, -2,  1,  1, -1, 1, 0, 1]])
    b = np.array([0, 0])
    C = -1*np.array([5, -10, -4, 6, -10, 10, 0, 0])
    IBV = np.array([6,7]) # index of the IBV
    sol = simplexit(A,b,C,IBV,debug)
    
    # calcul optimum not needed
    # z = C[sol[0]].dot(sol[1])
    return
In [12]:
collisions(debug=False)
Solution in unbounded!

Unbounded solution in the dual, so no solution in the primal, so no collision !

In [14]:
## Young couple
def youngcouple(debug=False):
    A = np.array([[ 1, 1, 1, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0],
                  [ 0, 0, 0, 0, 1, 1, 1, 1, 0, 1, 0, 0, 0, 0],
                  [ 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0],
                  [ 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0],
                  [ 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0],
                  [ 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1]
                 ])
    b = np.array([2, 2, 1, 1, 1, 1])
    C = np.array([4.5, 7.8, 3.6, 2.9, 4.9, 7.2, 4.3, 3.1, 100, 100, 100, 100, 100, 100])
    IBV = np.array([8,9,10,11,12,13])
    sol = simplexit(A,b,C,IBV,debug)
    z= C[sol[0]].dot(sol[1])
    return(z)
In [15]:
print(youngcouple())
 optimum found!
IBV =  [ 2  5  4 11  6  3]
Bbar=  [1. 1. 1. 0. 0. 1.]
18.6
/Users/talboth/Library/Python/3.6/lib/python/site-packages/ipykernel_launcher.py:77: RuntimeWarning: divide by zero encountered in true_divide
/Users/talboth/Library/Python/3.6/lib/python/site-packages/ipykernel_launcher.py:77: RuntimeWarning: invalid value encountered in true_divide
In [ ]:
## tighter problem formulation