simplex_primal.py

import numpy as np

# A = np.array([[-1, -1, 1, 0, 1, 0, 0],
#               [3, -1, 0, 0, 0, 1, 0],
#               [4, 2, 0, 1, 0, 0, 1]],
#              dtype = np.double)
# b = np.array([1, 20, 100], dtype = np.double)
# c = np.array([0, 0, 0, 0, 0, 1, 1, 1], dtype = np.double)

# B = np.array([4, 5, 6])
# N = np.array([0, 1, 2, 3])
# A = np.array([[1, 1, -1, 0, 1, 0, 0],
#               [3, -1, 0, 0, 0, 1, 0],
#               [4, 2, 0, 1, 0, 0, 1]],
#              dtype = np.double)
#
# c = np.array([0, 0, 0, 0, 1, 1, 1], dtype = np.double)
# b = np.array([1, 20, 100], dtype = np.double)
# B = np.array([4, 5, 6])
# N = np.array([0, 1, 2, 3])
A = np.array([
                    [0, 1, 2, 0],
                    [2, 1, 0, 2]]
    , dtype=np.double)
b = np.array([2, 4], dtype=np.double)
c = np.array([0, 0, 4, 1], dtype=np.double)
# variable
B = np.array([3, 1])
N = np.array([ 2, 4])
B = B - 1
N = N - 1
z = np.zeros(len(c))

print(np.linalg.inv(A[:, B]) @ b)


def simplex_primal(A, b, c, B, N):
    x = np.zeros(len(c), dtype=np.double)
    x[B] = np.linalg.inv(A[:, B]) @ b
    while True:
        print("Iteration")
        # BTRAN
        print("\tBTRAN")
        y = np.linalg.solve(A[:, B].T, c[B])
        print("\t\ty = {}".format(y))
        # Pricing
        print("\tPricing")
        z[N] = c[N] - A[:, N].T @ y
        if np.all(z[N] >= 0):
            print("B = {} ist optimal. x_B = {} .STOP".format(B + 1, x[B]))
            break
        j = np.where(z[N] < 0)[0][0] if np.any(z[N] < 0) else -1
        j = N[j]
        print("\t\tz_N = {}, j = {}".format(z[N], j + 1))
        print("\t\tx_j = x_{} = {} in die Basis eintretende Variable".format(j + 1, x[j]))
        # FTRAN
        print("\tFTRAN")
        w = np.linalg.solve(A[:, B], A[:, j])
        tolerance = 1e-6
        w = np.where(np.abs(w) < tolerance, 0, w)
        print("\t\tw = {}".format(w))
        # RATIO-TEST
        print("\tRatio-test")
        if np.all(w <= 0):
            print("THE PROBLEM IS UNBOUNDED. STOP")
        valid_indices = np.where(w > 0)[0]
        ratios = x[B][valid_indices] / w[valid_indices]
        gamma = np.min(ratios)
        i = valid_indices[np.argmin(ratios)]
        print("\t\tgamma = {}".format(gamma))
        print("\t\ti = {}, gamma = {}, x_B_i = x_{} verlaesst die Basis".format(i + 1, gamma, B[i] + 1))
        # Update
        print("\tUpdate")
        x[B] = x[B] - gamma * w

        index_to_replace = np.where(N == j)[0][0]
        N[index_to_replace] = B[i]
        B[i] = j
        x[j] = gamma
        print("\t\tx[B] = {}, N = {}, B = {}".format(x[B], N + 1, B + 1))


simplex_primal(A, b, c, B, N)