Untitled
unknown
python
3 years ago
2.3 kB
17
Indexable
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)Editor is loading...