simplex implemented

Author: Arthod

__pycache__/linear_programming.cpython-39.pyc

@@ -167,11 +167,7 @@ def add_constraint(self, expr_LHS: sympy.core.Expr, comparator: int, expr_RHS: s
         """
         assert comparator in COMPARATORS
 
-        if (comparator == EQUAL):
-            self.constraints.append(Constraint(expr_LHS, GREATER_EQUAL, expr_RHS))
-            self.constraints.append(Constraint(expr_LHS, LESS_EQUAL, expr_RHS))
-        else:
-            self.constraints.append(Constraint(expr_LHS, comparator, expr_RHS))
+        self.constraints.append(Constraint(expr_LHS, comparator, expr_RHS))
 
         # Update flags
         self._A_updated = False
@@ -228,11 +224,13 @@ def solve(self, method=SOLVER_SIMPLEX):
         if (method == SOLVER_SIMPLEX):
             lp_slacked = self.slacken_problem()
 
-            simplex(lp_slacked.A, lp_slacked.b, lp_slacked.c)
+            x_sol = simplex(lp_slacked.A, lp_slacked.b, lp_slacked.c, lp_slacked.vars_slack_amount)
 
+            # Cut slack variables
+            x_sol = x_sol[:lp_slacked.vars_slack_amount]
+            print(x_sol)
 
-
-            pass
+            return x_sol
 
         if (method == SOLVER_INTERIOR_POINT_METHOD):
             pass

linear_programming.py

@@ -3,32 +3,32 @@
 
 if __name__ == "__main__":    
 
+    """
     lp = LP.LinearProgrammingProblem()
     x1 = lp.add_variable(1)
     x2 = lp.add_variable(2)
     lp.set_objective(LP.MAX, 6 * x1 + 8 * x2)
     lp.add_constraint(5 * x1 + 10 * x2, LP.LESS_EQUAL, 60)
     lp.add_constraint(4 * x1 + 4 * x2, LP.LESS_EQUAL, 40)
-    
     """
-    lp = LP.LP()
+    lp = LP.LinearProgrammingProblem()
     x = lp.add_variable(1)
     y = lp.add_variable(2)
     z = lp.add_variable(3)
     lp.set_objective(LP.MAX, x + 2 * y + 3 * z)
     lp.add_constraint(x + y, LP.LESS_EQUAL, 20.5)
-    lp.add_constraint(y + z, LP.GREATER_EQUAL, 20.5)
-    lp.add_constraint(x + z, LP.GREATER_EQUAL, 30.5)
-    """
+    lp.add_constraint(y + z, LP.LESS_EQUAL, 20.5)
+    lp.add_constraint(x + z, LP.LESS_EQUAL, 32.5)
+    
 
 
     x = lp.solve()
 
     print("z = " + str(x @ lp.c))
     print("x = " + str(x))
-    print("iterations = " + str(lp.iterations))
-    if (len(lp.variables) == 2):
-        lp.plot_solution_path()
+    #print("iterations = " + str(lp.iterations))
+    #if (len(lp.variables) == 2):
+    #    lp.plot_solution_path()
 
 
 

main.py

@@ -12,32 +12,90 @@
 """
 
 
-def simplex(A: np.array, b: np.array, c: np.array) -> np.array:
+def simplex(A: np.array, b: np.array, c: np.array, slacks_amt: int) -> np.array:
     # https://sdu.itslearning.com/ContentArea/ContentArea.aspx?LocationID=17727&LocationType=1&ElementID=589350
-    x_current = basic_feasible_solution(A, b, c)
-    print(A, b, c)
-
-    ## Choosing a pivot
-    # Pick a non-negative column    TODO
-    positives = np.where(c > 0)
-    if (not len(positives[0])):
-        return "All coefficients negative"
-
-    x_index = positives[0][0]
-    y_index = None
-
-    val = float("inf")
-    a = A[:,x_index]
-    for i in range(b.size):
-        if (a[i] > 0):
-            val_current = b[i] / a[i]
-            if (val_current <= val):
-                y_index = i
-                val = val_current
-
-    print(y_index)
-    print(x_index)
+    #x_current = basic_feasible_solution(A, b, c)
+    N = [i for i in range(c.size - slacks_amt)]   # Not In basis
+    B = [i for i in range(slacks_amt, c.size)]    # In basis
 
+    print(A, b, c)
+    print(N, B)
+
+    while True:
+        ## Choosing a pivot
+        # Pick a non-negative column    TODO
+        positives = np.where(c > 0)
+        if (not len(positives[0])):
+            break
+
+        xx = positives[0][0]
+        yy = None
+        yy_old = None
+
+        # Picking smallest row b[i] / a[i]
+        val = float("inf")
+        a = A[:,xx]
+        for i in range(b.size):
+            if (a[i] > 0):
+                val_current = b[i] / a[i]
+                if (val_current <= val):
+                    yy = i
+                    val = val_current
+
+        # Find column from where basis was exited
+        a = A[yy,:]
+        for i in range(a.size):
+            if (a[i] == 1):
+                s = np.sum(A[:,i]) + c[i]
+                if (s == 1):
+                    yy_old = i
+                    break
+        # Remove yy_old from basis and add yy to basis
+        B.remove(yy_old)
+        B.append(yy)
+        N.append(yy_old)
+        N.remove(yy)
+
+        ## Pivot: A[xx, yy] - Peform row operations
+        # Set yy row with in column xx to 1
+        b[yy] = b[yy] * 1 / A[yy, xx]
+        A = row_mult(A, yy, 1 / A[yy, xx])
+
+        # Set column xx in c to 0
+        c = c - A[yy,:] * c[xx]
+
+        # Set other values in column xx to 0
+        for i in range(A.shape[0]):
+            if (i == yy):
+                continue
+            b[i] = b[i] - A[i, xx] * b[yy]
+            A = row_add(A, i, yy, -A[yy, xx] * A[i, xx])
+
+    # Get solution
+    x_sol = np.zeros(c.size)
+    for i in range(len(B)):
+        x_sol[B[i]] = b[i]
+    return np.array(x_sol)
+
+def row_mult(A: np.array, row_index: int, scalar: int) -> np.array:
+    # Row operation: R1 = scalar * R1
+    # Inspiration: (TODO - optimize)
+    # https://personal.math.ubc.ca/~pwalls/math-python/linear-algebra/solving-linear-systems/
+    I = np.eye(A.shape[0])
+    I[row_index, row_index] = scalar
+    return I @ A
+
+def row_add(A: np.array, row_main_index: int, row_other_index: int, scalar: int) -> np.array:
+    # Row operation: R1 = R1 + scalar * R2
+    # Inspiration: (TODO - optimize)
+    # https://personal.math.ubc.ca/~pwalls/math-python/linear-algebra/solving-linear-systems/
+    
+    I = np.eye(A.shape[0])
+    if (row_main_index == row_other_index):
+        I[row_main_index, row_main_index] = scalar + 1
+    else:
+        I[row_main_index, row_other_index] = scalar
+    return I @ A
 
 def basic_feasible_solution(A: np.array, b: np.array, c: np.array) -> np.array:
     zero_point = zero_point_solution(A, b, c)