""" Simplex method for solving linear programming problems """ import numpy as np from warnings import warn from .optimize import OptimizeResult, OptimizeWarning, _check_unknown_options from ._linprog_util import _postsolve def _pivot_col(T, tol=1.0E-12, bland=False): """ Given a linear programming simplex tableau, determine the column of the variable to enter the basis. Parameters ---------- T : 2D array A 2D array representing the simplex tableau, T, corresponding to the linear programming problem. It should have the form: [[A[0, 0], A[0, 1], ..., A[0, n_total], b[0]], [A[1, 0], A[1, 1], ..., A[1, n_total], b[1]], . . . [A[m, 0], A[m, 1], ..., A[m, n_total], b[m]], [c[0], c[1], ..., c[n_total], 0]] for a Phase 2 problem, or the form: [[A[0, 0], A[0, 1], ..., A[0, n_total], b[0]], [A[1, 0], A[1, 1], ..., A[1, n_total], b[1]], . . . [A[m, 0], A[m, 1], ..., A[m, n_total], b[m]], [c[0], c[1], ..., c[n_total], 0], [c'[0], c'[1], ..., c'[n_total], 0]] for a Phase 1 problem (a problem in which a basic feasible solution is sought prior to maximizing the actual objective. ``T`` is modified in place by ``_solve_simplex``. tol : float Elements in the objective row larger than -tol will not be considered for pivoting. Nominally this value is zero, but numerical issues cause a tolerance about zero to be necessary. bland : bool If True, use Bland's rule for selection of the column (select the first column with a negative coefficient in the objective row, regardless of magnitude). Returns ------- status: bool True if a suitable pivot column was found, otherwise False. A return of False indicates that the linear programming simplex algorithm is complete. col: int The index of the column of the pivot element. If status is False, col will be returned as nan. """ ma = np.ma.masked_where(T[-1, :-1] >= -tol, T[-1, :-1], copy=False) if ma.count() == 0: return False, np.nan if bland: return True, np.nonzero(ma.mask == False)[0][0] return True, np.ma.nonzero(ma == ma.min())[0][0] def _pivot_row(T, basis, pivcol, phase, tol=1.0E-12, bland=False): """ Given a linear programming simplex tableau, determine the row for the pivot operation. Parameters ---------- T : 2D array A 2D array representing the simplex tableau, T, corresponding to the linear programming problem. It should have the form: [[A[0, 0], A[0, 1], ..., A[0, n_total], b[0]], [A[1, 0], A[1, 1], ..., A[1, n_total], b[1]], . . . [A[m, 0], A[m, 1], ..., A[m, n_total], b[m]], [c[0], c[1], ..., c[n_total], 0]] for a Phase 2 problem, or the form: [[A[0, 0], A[0, 1], ..., A[0, n_total], b[0]], [A[1, 0], A[1, 1], ..., A[1, n_total], b[1]], . . . [A[m, 0], A[m, 1], ..., A[m, n_total], b[m]], [c[0], c[1], ..., c[n_total], 0], [c'[0], c'[1], ..., c'[n_total], 0]] for a Phase 1 problem (a Problem in which a basic feasible solution is sought prior to maximizing the actual objective. ``T`` is modified in place by ``_solve_simplex``. basis : array A list of the current basic variables. pivcol : int The index of the pivot column. phase : int The phase of the simplex algorithm (1 or 2). tol : float Elements in the pivot column smaller than tol will not be considered for pivoting. Nominally this value is zero, but numerical issues cause a tolerance about zero to be necessary. bland : bool If True, use Bland's rule for selection of the row (if more than one row can be used, choose the one with the lowest variable index). Returns ------- status: bool True if a suitable pivot row was found, otherwise False. A return of False indicates that the linear programming problem is unbounded. row: int The index of the row of the pivot element. If status is False, row will be returned as nan. """ if phase == 1: k = 2 else: k = 1 ma = np.ma.masked_where(T[:-k, pivcol] <= tol, T[:-k, pivcol], copy=False) if ma.count() == 0: return False, np.nan mb = np.ma.masked_where(T[:-k, pivcol] <= tol, T[:-k, -1], copy=False) q = mb / ma min_rows = np.ma.nonzero(q == q.min())[0] if bland: return True, min_rows[np.argmin(np.take(basis, min_rows))] return True, min_rows[0] def _apply_pivot(T, basis, pivrow, pivcol, tol=1e-12): """ Pivot the simplex tableau inplace on the element given by (pivrow, pivol). The entering variable corresponds to the column given by pivcol forcing the variable basis[pivrow] to leave the basis. Parameters ---------- T : 2D array A 2D array representing the simplex tableau, T, corresponding to the linear programming problem. It should have the form: [[A[0, 0], A[0, 1], ..., A[0, n_total], b[0]], [A[1, 0], A[1, 1], ..., A[1, n_total], b[1]], . . . [A[m, 0], A[m, 1], ..., A[m, n_total], b[m]], [c[0], c[1], ..., c[n_total], 0]] for a Phase 2 problem, or the form: [[A[0, 0], A[0, 1], ..., A[0, n_total], b[0]], [A[1, 0], A[1, 1], ..., A[1, n_total], b[1]], . . . [A[m, 0], A[m, 1], ..., A[m, n_total], b[m]], [c[0], c[1], ..., c[n_total], 0], [c'[0], c'[1], ..., c'[n_total], 0]] for a Phase 1 problem (a problem in which a basic feasible solution is sought prior to maximizing the actual objective. ``T`` is modified in place by ``_solve_simplex``. basis : 1D array An array of the indices of the basic variables, such that basis[i] contains the column corresponding to the basic variable for row i. Basis is modified in place by _apply_pivot. pivrow : int Row index of the pivot. pivcol : int Column index of the pivot. """ basis[pivrow] = pivcol pivval = T[pivrow, pivcol] T[pivrow] = T[pivrow] / pivval for irow in range(T.shape[0]): if irow != pivrow: T[irow] = T[irow] - T[pivrow] * T[irow, pivcol] # The selected pivot should never lead to a pivot value less than the tol. if np.isclose(pivval, tol, atol=0, rtol=1e4): message = ( "The pivot operation produces a pivot value of:{0: .1e}, " "which is only slightly greater than the specified " "tolerance{1: .1e}. This may lead to issues regarding the " "numerical stability of the simplex method. " "Removing redundant constraints, changing the pivot strategy " "via Bland's rule or increasing the tolerance may " "help reduce the issue.".format(pivval, tol)) warn(message, OptimizeWarning) def _solve_simplex(T, n, basis, maxiter=1000, phase=2, status=0, message='', callback=None, tol=1.0E-12, nit0=0, bland=False, _T_o=None): """ Solve a linear programming problem in "standard form" using the Simplex Method. Linear Programming is intended to solve the following problem form: Minimize:: c @ x Subject to:: A @ x == b x >= 0 Parameters ---------- T : 2D array A 2D array representing the simplex tableau, T, corresponding to the linear programming problem. It should have the form: [[A[0, 0], A[0, 1], ..., A[0, n_total], b[0]], [A[1, 0], A[1, 1], ..., A[1, n_total], b[1]], . . . [A[m, 0], A[m, 1], ..., A[m, n_total], b[m]], [c[0], c[1], ..., c[n_total], 0]] for a Phase 2 problem, or the form: [[A[0, 0], A[0, 1], ..., A[0, n_total], b[0]], [A[1, 0], A[1, 1], ..., A[1, n_total], b[1]], . . . [A[m, 0], A[m, 1], ..., A[m, n_total], b[m]], [c[0], c[1], ..., c[n_total], 0], [c'[0], c'[1], ..., c'[n_total], 0]] for a Phase 1 problem (a problem in which a basic feasible solution is sought prior to maximizing the actual objective. ``T`` is modified in place by ``_solve_simplex``. n : int The number of true variables in the problem. basis : 1D array An array of the indices of the basic variables, such that basis[i] contains the column corresponding to the basic variable for row i. Basis is modified in place by _solve_simplex maxiter : int The maximum number of iterations to perform before aborting the optimization. phase : int The phase of the optimization being executed. In phase 1 a basic feasible solution is sought and the T has an additional row representing an alternate objective function. callback : callable, optional (simplex only) If a callback function is provided, it will be called within each iteration of the simplex algorithm. The callback must require a `scipy.optimize.OptimizeResult` consisting of the following fields: x : 1D array The independent variable vector which optimizes the linear programming problem. fun : float Value of the objective function. success : bool True if the algorithm succeeded in finding an optimal solution. slack : 1D array The values of the slack variables. Each slack variable corresponds to an inequality constraint. If the slack is zero, the corresponding constraint is active. con : 1D array The (nominally zero) residuals of the equality constraints, that is, ``b - A_eq @ x`` phase : int The phase of the optimization being executed. In phase 1 a basic feasible solution is sought and the T has an additional row representing an alternate objective function. status : int An integer representing the exit status of the optimization:: 0 : Optimization terminated successfully 1 : Iteration limit reached 2 : Problem appears to be infeasible 3 : Problem appears to be unbounded 4 : Serious numerical difficulties encountered nit : int The number of iterations performed. message : str A string descriptor of the exit status of the optimization. tol : float The tolerance which determines when a solution is "close enough" to zero in Phase 1 to be considered a basic feasible solution or close enough to positive to serve as an optimal solution. nit0 : int The initial iteration number used to keep an accurate iteration total in a two-phase problem. bland : bool If True, choose pivots using Bland's rule [3]_. In problems which fail to converge due to cycling, using Bland's rule can provide convergence at the expense of a less optimal path about the simplex. Returns ------- nit : int The number of iterations. Used to keep an accurate iteration total in the two-phase problem. status : int An integer representing the exit status of the optimization:: 0 : Optimization terminated successfully 1 : Iteration limit reached 2 : Problem appears to be infeasible 3 : Problem appears to be unbounded 4 : Serious numerical difficulties encountered """ nit = nit0 complete = False if phase == 1: m = T.shape[0]-2 elif phase == 2: m = T.shape[0]-1 else: raise ValueError("Argument 'phase' to _solve_simplex must be 1 or 2") if phase == 2: # Check if any artificial variables are still in the basis. # If yes, check if any coefficients from this row and a column # corresponding to one of the non-artificial variable is non-zero. # If found, pivot at this term. If not, start phase 2. # Do this for all artificial variables in the basis. # Ref: "An Introduction to Linear Programming and Game Theory" # by Paul R. Thie, Gerard E. Keough, 3rd Ed, # Chapter 3.7 Redundant Systems (pag 102) for pivrow in [row for row in range(basis.size) if basis[row] > T.shape[1] - 2]: non_zero_row = [col for col in range(T.shape[1] - 1) if abs(T[pivrow, col]) > tol] if len(non_zero_row) > 0: pivcol = non_zero_row[0] _apply_pivot(T, basis, pivrow, pivcol) nit += 1 if len(basis[:m]) == 0: solution = np.zeros(T.shape[1] - 1, dtype=np.float64) else: solution = np.zeros(max(T.shape[1] - 1, max(basis[:m]) + 1), dtype=np.float64) while not complete: # Find the pivot column pivcol_found, pivcol = _pivot_col(T, tol, bland) if not pivcol_found: pivcol = np.nan pivrow = np.nan status = 0 complete = True else: # Find the pivot row pivrow_found, pivrow = _pivot_row(T, basis, pivcol, phase, tol, bland) if not pivrow_found: status = 3 complete = True if callback is not None: solution[basis[:n]] = T[:n, -1] x = solution[:m] c, A_ub, b_ub, A_eq, b_eq, bounds, undo = _T_o x, fun, slack, con, _, _ = _postsolve( x, c, A_ub, b_ub, A_eq, b_eq, bounds, undo=undo, tol=tol ) res = OptimizeResult({ 'x': x, 'fun': fun, 'slack': slack, 'con': con, 'status': status, 'message': message, 'nit': nit, 'success': status == 0 and complete, 'phase': phase, 'complete': complete, }) callback(res) if not complete: if nit >= maxiter: # Iteration limit exceeded status = 1 complete = True else: _apply_pivot(T, basis, pivrow, pivcol) nit += 1 return nit, status def _linprog_simplex(c, c0, A, b, maxiter=1000, disp=False, callback=None, tol=1.0E-12, bland=False, _T_o=None, **unknown_options): """ Minimize a linear objective function subject to linear equality and non-negativity constraints using the two phase simplex method. Linear programming is intended to solve problems of the following form: Minimize:: c @ x Subject to:: A @ x == b x >= 0 Parameters ---------- c : 1D array Coefficients of the linear objective function to be minimized. c0 : float Constant term in objective function due to fixed (and eliminated) variables. (Purely for display.) A : 2D array 2D array such that ``A @ x``, gives the values of the equality constraints at ``x``. b : 1D array 1D array of values representing the right hand side of each equality constraint (row) in ``A``. callback : callable, optional (simplex only) If a callback function is provided, it will be called within each iteration of the simplex algorithm. The callback must require a `scipy.optimize.OptimizeResult` consisting of the following fields: x : 1D array The independent variable vector which optimizes the linear programming problem. fun : float Value of the objective function. success : bool True if the algorithm succeeded in finding an optimal solution. slack : 1D array The values of the slack variables. Each slack variable corresponds to an inequality constraint. If the slack is zero, the corresponding constraint is active. con : 1D array The (nominally zero) residuals of the equality constraints, that is, ``b - A_eq @ x`` phase : int The phase of the optimization being executed. In phase 1 a basic feasible solution is sought and the T has an additional row representing an alternate objective function. status : int An integer representing the exit status of the optimization:: 0 : Optimization terminated successfully 1 : Iteration limit reached 2 : Problem appears to be infeasible 3 : Problem appears to be unbounded 4 : Serious numerical difficulties encountered nit : int The number of iterations performed. message : str A string descriptor of the exit status of the optimization. Options ------- maxiter : int The maximum number of iterations to perform. disp : bool If True, print exit status message to sys.stdout tol : float The tolerance which determines when a solution is "close enough" to zero in Phase 1 to be considered a basic feasible solution or close enough to positive to serve as an optimal solution. bland : bool If True, use Bland's anti-cycling rule [3]_ to choose pivots to prevent cycling. If False, choose pivots which should lead to a converged solution more quickly. The latter method is subject to cycling (non-convergence) in rare instances. Returns ------- x : 1D array Solution vector. status : int An integer representing the exit status of the optimization:: 0 : Optimization terminated successfully 1 : Iteration limit reached 2 : Problem appears to be infeasible 3 : Problem appears to be unbounded 4 : Serious numerical difficulties encountered message : str A string descriptor of the exit status of the optimization. iteration : int The number of iterations taken to solve the problem. References ---------- .. [1] Dantzig, George B., Linear programming and extensions. Rand Corporation Research Study Princeton Univ. Press, Princeton, NJ, 1963 .. [2] Hillier, S.H. and Lieberman, G.J. (1995), "Introduction to Mathematical Programming", McGraw-Hill, Chapter 4. .. [3] Bland, Robert G. New finite pivoting rules for the simplex method. Mathematics of Operations Research (2), 1977: pp. 103-107. Notes ----- The expected problem formulation differs between the top level ``linprog`` module and the method specific solvers. The method specific solvers expect a problem in standard form: Minimize:: c @ x Subject to:: A @ x == b x >= 0 Whereas the top level ``linprog`` module expects a problem of form: Minimize:: c @ x Subject to:: A_ub @ x <= b_ub A_eq @ x == b_eq lb <= x <= ub where ``lb = 0`` and ``ub = None`` unless set in ``bounds``. The original problem contains equality, upper-bound and variable constraints whereas the method specific solver requires equality constraints and variable non-negativity. ``linprog`` module converts the original problem to standard form by converting the simple bounds to upper bound constraints, introducing non-negative slack variables for inequality constraints, and expressing unbounded variables as the difference between two non-negative variables. """ _check_unknown_options(unknown_options) status = 0 messages = {0: "Optimization terminated successfully.", 1: "Iteration limit reached.", 2: "Optimization failed. Unable to find a feasible" " starting point.", 3: "Optimization failed. The problem appears to be unbounded.", 4: "Optimization failed. Singular matrix encountered."} n, m = A.shape # All constraints must have b >= 0. is_negative_constraint = np.less(b, 0) A[is_negative_constraint] *= -1 b[is_negative_constraint] *= -1 # As all constraints are equality constraints the artificial variables # will also be basic variables. av = np.arange(n) + m basis = av.copy() # Format the phase one tableau by adding artificial variables and stacking # the constraints, the objective row and pseudo-objective row. row_constraints = np.hstack((A, np.eye(n), b[:, np.newaxis])) row_objective = np.hstack((c, np.zeros(n), c0)) row_pseudo_objective = -row_constraints.sum(axis=0) row_pseudo_objective[av] = 0 T = np.vstack((row_constraints, row_objective, row_pseudo_objective)) nit1, status = _solve_simplex(T, n, basis, phase=1, callback=callback, maxiter=maxiter, tol=tol, bland=bland, _T_o=_T_o) # if pseudo objective is zero, remove the last row from the tableau and # proceed to phase 2 if abs(T[-1, -1]) < tol: # Remove the pseudo-objective row from the tableau T = T[:-1, :] # Remove the artificial variable columns from the tableau T = np.delete(T, av, 1) else: # Failure to find a feasible starting point status = 2 nit2 = nit1 messages[status] = ( "Phase 1 of the simplex method failed to find a feasible " "solution. The pseudo-objective function evaluates to {0:.1e} " "which exceeds the required tolerance of {1} for a solution to be " "considered 'close enough' to zero to be a basic solution. " "Consider increasing the tolerance to be greater than {0:.1e}. " "If this tolerance is unacceptably large the problem may be " "infeasible.".format(abs(T[-1, -1]), tol) ) if status == 0: # Phase 2 nit2, status = _solve_simplex(T, n, basis, maxiter=maxiter, phase=2, callback=callback, tol=tol, nit0=nit1, bland=bland, _T_o=_T_o) solution = np.zeros(n + m) solution[basis[:n]] = T[:n, -1] x = solution[:m] return x, status, messages[status], int(nit2)