dual_simplex.py

import numpy as np

# # constants
A = np.array([[1, 2, -2, 0],
              [1, 1, -1, -1]], dtype = np.double)
b = np.array([1, 1], dtype = np.double)
c = np.array([1, -1, 1, 0], dtype = np.double)
# variable
B = np.array([2, 4])
N = np.array([1, 3])

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

#initialization
B = B - 1
N = N - 1
x = np.zeros(len(c), dtype = np.double)
z = np.zeros(len(c), dtype = np.double)
z[N] = c[N] - A[:, N].T @ np.linalg.inv(A[:, B]).T @ c[B]
while(True):
    print("Iteration")
    # BTRAN
    print("\tBTRAN:")
    x_basis = np.linalg.solve(A[:, B], b)
    print("\tB = {}".format(B + 1))
    print("\t A[:, B] = {}".format(A[:, B]))
    print("\t\tx_B = x_{} = {}".format(B + 1, x_basis))
    x[B] = x_basis
    print()
    # Pricing
    print("\tPricing:")
    if np.all(x[B] >= 0):
        x[N] = 0
        print("OPTIMAL SOLUTION IS x = {}. END".format(x))
        break
    i = -1
    for index, value in enumerate(x[B]):
        if value < 0:
            i = index
            break
    print("\t\ti = {}, x_B_i = {} < 0, B_i = {} verlaesst die Basis B".format(i + 1, x[B][i], B[i] + 1))
    print()
    # FTRAN
    print("\tFTRAN")
    ith_unit_vector = np.zeros(len(B))
    ith_unit_vector[i] = 1
    print("\t\te_i = {}".format(ith_unit_vector))
    w = np.linalg.solve(A[:, B].T, ith_unit_vector)
    print("\t\tw = {}".format(w))
    alpha_N = -A[:, N].T @ w
    print("\t\talpha_N = alpha_{} = {}".format(N + 1, alpha_N))
    print()
    #Ratio-Test
    print("\tRatio-Test")
    if np.all(alpha_N <= 0):
        print("THE PROBLEM IS UNBOUNDED. END")
        break
    valid_indices = np.where(alpha_N > 0)[0]
    ratios = z[N][valid_indices] / alpha_N[valid_indices]
    gamma = np.min(ratios)
    print("\t\tz[N] = {}".format(z[N]))
    j = valid_indices[np.argmin(ratios)]
    print("\t\tj = {}, gamma = {}, N_j = {} tritt in die Basis ein".format(j + 1, gamma, N[j] + 1))
    print()
    #Update
    print("\tUPDATE")
    z[N] = z[N] - gamma * alpha_N
    z[B[i]] = gamma
    tmp = B[i]
    B[i] = N[j]
    N[j] = tmp
    print("\t\tB = {}, N = {}, z = {}".format(B + 1, N + 1, z))
    print()