Untitled

mail@pastecode.io avatar
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)