Untitled
unknown
python
2 years ago
2.3 kB
11
Indexable
Never
import numpy as np def simplex(c, A, s, b): A_ext = np.hstack([A, np.eye(A.shape[0])]) c_ext = np.hstack([c, np.zeros(A.shape[0])]) base = np.arange(len(c), len(c_ext)) non_base = np.delete(np.arange(0, len(c_ext)), base) print("Initial tableau") print(A_ext) print(c_ext) print(base) print(non_base) while True: ### The loop starts here reduced_costs = np.ravel(c_ext[non_base] - c_ext[base]*A_ext[:, base]*A_ext[:, non_base]) print(f"reduced costs = {reduced_costs}") if np.all(reduced_costs <= 1.0e-6) : break entering_col = np.argmax(reduced_costs) print(f"entering col = {entering_col}") entering_limits = np.ravel(np.transpose(b*A_ext[:, base]) / (A_ext[: , base]*A_ext[:, entering_col])) if(entering_limits.any() < 0): break print(f"Entering limits = {entering_limits}") pivot_row_index = np.argmin(entering_limits) leaving_col = base[pivot_row_index] print(f"Leaving row = { pivot_row_index}") # Update base base[pivot_row_index] = entering_col # Pivot A_ext_factors = np.ravel(A_ext[:, entering_col]* (-1/A_ext[pivot_row_index, entering_col])) A_ext_factors[pivot_row_index] = 1/A_ext[pivot_row_index, entering_col] print(f"Pivot factors = {A_ext_factors}") # Pivot obj c_ext = np.ravel(c_ext - (c_ext[entering_col]/A_ext[pivot_row_index, entering_col])*A_ext[pivot_row_index, :]) # Pivot obj b_pivot_row = b[pivot_row_index] b = b + np.multiply(np.transpose(A_ext_factors),b) b[pivot_row_index] = b_pivot_row*A_ext_factors[pivot_row_index] # Pivot A matrix pivot_row = np.array(A_ext[pivot_row_index, : ]) A_ext[pivot_row_index, : ] = 0 A_ext = A_ext + np.multiply(A_ext_factors.reshape(len(A_ext_factors), 1), np.repeat(pivot_row, len(A_ext), axis=0)) A_ext[pivot_row_index, : ] = A_ext_factors[pivot_row_index]*pivot_row print(f"New cost array \t\t {c_ext}") print(f"New vector b \t\t {b}") print(f"New A_ext \t\t\n {np.around(A_ext,3)}") A = np.matrix([[2, 3, 5], [3, 4, 6], [2, 9, 7], [4, 1, 3]], dtype=np.double) c = np.array([1, 2, 6], dtype=np.double) s = np.array(['l', 'e', 'g']) b = np.array([3, 8, 9, 2], dtype=np.double) simplex(c, A, s, b)