simplex-algorithm-implementation.ipynb

{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.10.13","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[],"dockerImageVersionId":30698,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"import numpy as np\nfrom numpy.linalg import inv\nfrom numpy.matlib import matrix  \n\nnp.set_printoptions(precision = 3, threshold = 10, edgeitems = 4, linewidth = 120)  \n\nTRUNC_THRSHLD = 10**(-12)\n\ndef LPP(A: matrix, b: np.array, c: np.array, \\\n            _EPSILON : float) -> (int, np.array, float, np.array):\n    \"\"\"\n    executes simplex method: phase 1 and 2\n    @param A: constraint matrix\n    @param b: independent terms in constraints\n    @param c: costs vector\n    \"\"\"\n    __m, __n = A.shape[0], A.shape[1]  # no. of rows, columns of A\n\n    if __n < __m: raise ValueError(\"Incompatible dimensions \"\n                         \"(no. of variables : {} > {} : no.of constraints\".format(n, m))\n    \n    if b.shape != (__m,): raise \\\n        ValueError(\"Incompatible dimensions: c_j has shape {}, expected {}.\".format(b.shape, (m,)))\n    \n    if c.shape != (__n,): raise \\\n        ValueError(\"Incompatible dimensions: c has shape {}, expected {}.\".format(c.shape, (n,)))\n\n    \"Check full rank matrix\"\n    if not np.linalg.matrix_rank(A) == __m:\n        # Remove ld rows:\n        A = A[[i for i in range(__m) if not np.array_equal(np.linalg.qr(A)[1][i, :], np.zeros(__n))], :]\n        __m = A.shape[0]  # Update no. of rows\n\n    A[[i for i in range(m) if b[i] < 0]] *= -1  # Change sign of constraints\n    b = np.abs(b)\n\n    A_I = matrix(np.concatenate((A, np.identity(__m)), axis = 1))  # Phase 1 constraint matrix\n    x_I = np.concatenate((np.zeros(__n), b))  # Phase 1 variable vector\n    c_I = np.concatenate((np.zeros(__n), np.ones(__m)))  # Phase 1 c_j vector\n    basic_I = set(range(__n, __n + __m))  # Phase 1 basic variable set\n\n    print(\"Executing phase 1\")\n    ext_code_I, x_init, basic_init, z_I, _, ITERS_I = simplex_core(A_I, c_I, x_I, basic_I)\n    # ^ Exit code, initial BFS, basis, z_I, and no. of iterations\n    print(\"Phase I terminated.\")\n\n    assert ext_code_I == 0  # assert phase I has optimal solution\n    if z_I > 0:\n        print(\"\\n\")\n        print_boxed(\"Unfeasible problem (z_I = {:.6g} > 0).\".format(z_I))\n        print(\"{} iterations in phase I.\".format(ITERS_I), end = '\\n\\n')\n        return 2, None, None, None\n    \n    # If artificial variable in the basis for the initial BFS, exit:\n    if any(j not in range(__n) for j in basic_init): raise NotImplementedError(\"Artificial variables in basis\")\n\n    x_init = x_init[:__n]  # Get initial BFS for original problem\n\n    print(\"Found initial BFS at x = \\n{}.\\n\".format(x_init))\n\n    \"\"\"Phase II execution\"\"\"\n    print(\"Executing phase 2\")\n    ext_code, x, basic, z, d, ITERS_II = simplex_core(A, c, x_init, basic_init)\n    print(\"Phase II terminated.\\n\")\n\n    if ext_code == 0:\n        print_boxed(\"Found optimal solution at x =\\n{}.\\n\\n\".format(x) +\n                    \"Basic indexes: {}\\n\".format(basic) +\n                    \"Nonbasic indexes: {}\\n\\n\".format(set(range(n)) - basic) +\n                    \"Optimal cost: {}.\".format(z))\n    \n    elif ext_code == 1: print_boxed(\n        \"Unbounded problem. Found feasible direction d =\\n{}\\nfrom x =\\n{}.\".format(d, x)\n    )\n\n    print(\"{} iterations in phase I, {} iterations in phase II ({} total).\".format(\n        ITERS_I, ITERS_II, ITERS_I + ITERS_II), end = '\\n\\n')\n\n    return ext_code, x, z, d","metadata":{"execution":{"iopub.status.busy":"2024-04-24T11:21:12.605620Z","iopub.execute_input":"2024-04-24T11:21:12.606096Z","iopub.status.idle":"2024-04-24T11:21:12.634835Z","shell.execute_reply.started":"2024-04-24T11:21:12.606050Z","shell.execute_reply":"2024-04-24T11:21:12.633557Z"},"trusted":true},"execution_count":2,"outputs":[]},{"cell_type":"code","source":"MAX_ITERS : int = 100000007 \n\ndef simplex_core(A: matrix, c: np.array, x: np.array, basic: set) \\\n        -> (int, np.array, set, float, np.array):\n    \"\"\"\n    @param A: constraint matrix\n    @param c: costs vector\n    @param x: initial BFS\n    @param basic: initial basic index set\n    \n    @return: a tuple consisting of the exit code, the value of x, basic index set,\n    optimal cost, and BFS corresponding to feasible direction\n    \"\"\"\n    __m, __n = A.shape[0], A.shape[1]  \n\n    assert c.shape == (__n,) and x.shape == (__n,)\n    \n    assert isinstance(basic, set) and len(basic) == __m and \\\n           all(i in range(__n) for i in basic)  # Make sure that basic is a valid base\n\n    B, N = list(basic), set(range(n)) - basic  # Basic /nonbasic index lists\n    del basic\n    B_inv = inv(A[:, B])  # Calculate inverse of basic matrix (`A[:, B]`)\n\n    z = np.dot(c, x) # Objective function\n\n    iteration : int = 1 \n    while iteration <= MAX_ITERS: \n        r_q, q, p, theta, d = None, None, None, None, None\n        print(\"\\tIteration no. {}:\".format(iteration), end='')\n\n        \"\"\"Optimality test\"\"\"\n        prices = c[B] * B_inv  # Store product for efficiency\n\n       # Bland rule\n        optimum = True\n        for q in N:  # Read in lexicographical index order\n            r_q = np.ndarray.item(c[q] - prices * A[:, q])\n            if r_q < 0:\n                optimum = False\n                break \n        \n        # Minimal reduced cost rule\n        #r_q, q = min([(np.ndarray.item(c[q] - prices * A[:, q]), q) \\ \n        #              for q in N], key=(lambda tup: tup[0]))\n        #optimum = (r_q >= 0)\n\n        if optimum:\n            print(\"\\tfound optimum\")\n            return 0, x, set(B), z, None, iteration  # Found optimal solution\n\n\n        \"\"\"Feasible basic direction\"\"\"\n        d = np.zeros(__n)\n        for i in range(__m): d[B[i]] = trunc(np.ndarray.item(-B_inv[i, :] * A[:, q]))\n        d[q] = 1\n\n\n        \"\"\"Max Step Length\"\"\"\n        # List of tuples of \"candidate\" theta an corresponding index in basic list:\n        neg = [(-x[B[i]] / d[B[i]], i) for i in range(__m) if d[B[i]] < 0]\n\n        if len(neg) == 0:\n            print(\"\\tidentified unbounded problem\")\n            return 1, x, set(B),  None, d, iteration  # Flag problem as unbounded and return direction\n\n        # Get theta and index (in basis) of exiting basic variable:\n        theta, p = min(neg, key=(lambda tup: tup[0]))\n\n\n        \"\"\"Variable updates\"\"\"\n        x = np.array([trunc(var) for var in (x + theta * d)])  # Update all variables\n        assert x[B[p]] == 0\n\n        z = trunc(z + theta * r_q)  # Update objective function value\n\n        # Update inverse:\n        for i in set(range(m)) - {p}:\n            B_inv[i, :] -= d[B[i]]/d[B[p]] * B_inv[p, :]\n        B_inv[p, :] /= -d[B[p]]\n\n        N = N - {q} | {B[p]}  # Update nonbasic index set\n        B[p] = q  # Update basic index list\n\n        \"\"\"Print status update\"\"\"\n        print(\n            \"\\tq = {:>2} \\trq = {:>9.2f} \\tB[p] = {:>2d} \\ttheta* = {:>5.4f} \\tz = {:<9.2f}\"\n                .format(q + 1, r_q, B[p] + 1, theta, z)\n        )\n        iteration += 1\n    raise TimeoutError(\"Iterations maxed out (probably due to an endless loop)\")","metadata":{"execution":{"iopub.status.busy":"2024-04-24T11:21:17.420439Z","iopub.execute_input":"2024-04-24T11:21:17.420850Z","iopub.status.idle":"2024-04-24T11:21:17.442627Z","shell.execute_reply.started":"2024-04-24T11:21:17.420821Z","shell.execute_reply":"2024-04-24T11:21:17.441343Z"},"trusted":true},"execution_count":3,"outputs":[]},{"cell_type":"code","source":"def print_boxed(msg: str) -> None:\n    lines = msg.splitlines()\n    max_len = max(len(line) for line in lines)\n\n    if max_len > 100: raise ValueError(\"Overfull box\")\n\n    print('-' * (max_len + 4))\n    for line in lines: print('| ' + line + ' ' * (max_len - len(line)) + ' |')\n    print('-' * (max_len + 4))\n\ndef trunc(x: float) -> float: return x if abs(x) >= TRUNC_THRSHLD else 0","metadata":{"execution":{"iopub.status.busy":"2024-04-24T11:21:21.475653Z","iopub.execute_input":"2024-04-24T11:21:21.476687Z","iopub.status.idle":"2024-04-24T11:21:21.486850Z","shell.execute_reply.started":"2024-04-24T11:21:21.476638Z","shell.execute_reply":"2024-04-24T11:21:21.485088Z"},"trusted":true},"execution_count":4,"outputs":[]},{"cell_type":"code","source":"","metadata":{"execution":{"iopub.status.busy":"2024-04-24T11:22:21.129812Z","iopub.execute_input":"2024-04-24T11:22:21.130858Z","iopub.status.idle":"2024-04-24T11:22:21.147107Z","shell.execute_reply.started":"2024-04-24T11:22:21.130820Z","shell.execute_reply":"2024-04-24T11:22:21.145639Z"},"trusted":true},"execution_count":6,"outputs":[{"name":"stdout","text":"Executing phase 1\n\tIteration no. 1:\tq =  1 \trq =     -1.00 \tB[p] =  1 \ttheta* = 3.0000 \tz = 4.00     \n\tIteration no. 2:\tq =  2 \trq =     -2.50 \tB[p] =  2 \ttheta* = 1.6000 \tz = 0.00     \n\tIteration no. 3:\tfound optimum\nPhase I terminated.\nFound initial BFS at x = \n[0.6 1.6 0.  0. ].\n\nExecuting phase 2\n\tIteration no. 1:\tfound optimum\nPhase II terminated.\n\n---------------------------------\n| Found optimal solution at x = |\n| [0.6 1.6 0.  0. ].            |\n|                               |\n| Basic indexes: {0, 1}         |\n| Nonbasic indexes: {2, 3}      |\n|                               |\n| Optimal cost: -5.4.           |\n---------------------------------\n3 iterations in phase I, 1 iterations in phase II (4 total).\n\n","output_type":"stream"},{"execution_count":6,"output_type":"execute_result","data":{"text/plain":"(0, array([0.6, 1.6, 0. , 0. ]), -5.4, None)"},"metadata":{}}]}]}