EVOLUTION-MANAGER

Edit File: _linprog_simplex.py

"""Simplex method for linear programming The *simplex* method uses a traditional, full-tableau implementation of Dantzig's simplex algorithm [1]_, [2]_ (*not* the Nelder-Mead simplex). This algorithm is included for backwards compatibility and educational purposes. .. versionadded:: 0.15.0 Warnings -------- The simplex method may encounter numerical difficulties when pivot values are close to the specified tolerance. If encountered try remove any redundant constraints, change the pivot strategy to Bland's rule or increase the tolerance value. Alternatively, more robust methods maybe be used. See :ref:`'interior-point' <optimize.linprog-interior-point>` and :ref:`'revised simplex' <optimize.linprog-revised_simplex>`. 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. """ 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=1e-9, 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: # ma.mask is sometimes 0d return True, np.nonzero(np.logical_not(np.atleast_1d(ma.mask)))[0][0] return True, np.ma.nonzero(ma == ma.min())[0][0] def _pivot_row(T, basis, pivcol, phase, tol=1e-9, 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-9): """ 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, stacklevel=5) def _solve_simplex(T, n, basis, callback, postsolve_args, maxiter=1000, tol=1e-9, phase=2, bland=False, nit0=0, ): """ 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 callback : callable, optional If a callback function is provided, it will be called within each iteration of the algorithm. The callback must accept a `scipy.optimize.OptimizeResult` consisting of the following fields: x : 1D array Current solution vector fun : float Current value of the objective function success : bool True only when a phase has completed successfully. This will be False for most iterations. 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. postsolve_args : tuple Data needed by _postsolve to convert the solution to the standard-form problem into the solution to the original problem. maxiter : int The maximum number of iterations to perform before aborting 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. 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. 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. nit0 : int The initial iteration number used to keep an accurate iteration total in a two-phase problem. 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 status = 0 message = '' complete = False if phase == 1: m = T.shape[1]-2 elif phase == 2: m = T.shape[1]-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, tol) 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[:] = 0 solution[basis[:n]] = T[:n, -1] x = solution[:m] x, fun, slack, con, _, _ = _postsolve( x, postsolve_args, 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, tol) nit += 1 return nit, status def _linprog_simplex(c, c0, A, b, callback, postsolve_args, maxiter=1000, tol=1e-9, disp=False, bland=False, **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 If a callback function is provided, it will be called within each iteration of the algorithm. The callback function must accept a single `scipy.optimize.OptimizeResult` consisting of the following fields: x : 1D array Current solution vector fun : float Current value of the objective function success : bool True when an algorithm has completed successfully. 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 algorithm being executed. status : int An integer representing the status of the optimization:: 0 : Algorithm proceeding nominally 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. postsolve_args : tuple Data needed by _postsolve to convert the solution to the standard-form problem into the solution to the original problem. 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. unknown_options : dict Optional arguments not used by this particular solver. If `unknown_options` is non-empty a warning is issued listing all unused options. 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, callback=callback, postsolve_args=postsolve_args, maxiter=maxiter, tol=tol, phase=1, bland=bland ) # if pseudo objective is zero, remove the last row from the tableau and # proceed to phase 2 nit2 = nit1 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 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, callback=callback, postsolve_args=postsolve_args, maxiter=maxiter, tol=tol, phase=2, bland=bland, nit0=nit1 ) solution = np.zeros(n + m) solution[basis[:n]] = T[:n, -1] x = solution[:m] return x, status, messages[status], int(nit2)