| """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 : 2-D array |
| A 2-D 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(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 : 2-D array |
| A 2-D 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 : 2-D array |
| A 2-D 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 : 1-D 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] |
|
|
| |
| if np.isclose(pivval, tol, atol=0, rtol=1e4): |
| message = ( |
| f"The pivot operation produces a pivot value of:{pivval: .1e}, " |
| "which is only slightly greater than the specified " |
| f"tolerance{tol: .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.") |
| 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 : 2-D array |
| A 2-D 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 : 1-D 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 : 1-D 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 : 1-D 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 : 1-D 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: |
| |
| |
| |
| |
| |
| |
| |
| |
| 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.empty(T.shape[1] - 1, dtype=np.float64) |
| else: |
| solution = np.empty(max(T.shape[1] - 1, max(basis[:m]) + 1), |
| dtype=np.float64) |
|
|
| while not complete: |
| |
| pivcol_found, pivcol = _pivot_col(T, tol, bland) |
| if not pivcol_found: |
| pivcol = np.nan |
| pivrow = np.nan |
| status = 0 |
| complete = True |
| else: |
| |
| 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 |
| ) |
| 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: |
| |
| 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 |
| |
| User-facing documentation is in _linprog_doc.py. |
| |
| Parameters |
| ---------- |
| c : 1-D 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 : 2-D array |
| 2-D array such that ``A @ x``, gives the values of the equality |
| constraints at ``x``. |
| b : 1-D array |
| 1-D 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 : 1-D array |
| Current solution vector |
| fun : float |
| Current value of the objective function |
| success : bool |
| True when an algorithm has completed successfully. |
| slack : 1-D 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 : 1-D 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 : 1-D 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 |
|
|
| |
| is_negative_constraint = np.less(b, 0) |
| A[is_negative_constraint] *= -1 |
| b[is_negative_constraint] *= -1 |
|
|
| |
| |
| av = np.arange(n) + m |
| basis = av.copy() |
|
|
| |
| |
| 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 |
| ) |
| |
| |
| nit2 = nit1 |
| if abs(T[-1, -1]) < tol: |
| |
| T = T[:-1, :] |
| |
| T = np.delete(T, av, 1) |
| else: |
| |
| 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: |
| |
| 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) |
|
|
