Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-packages/scipy/optimize/_linprog_simplex.py: 8%
107 statements
« prev ^ index » next coverage.py v7.3.2, created at 2023-12-12 06:31 +0000
1"""Simplex method for linear programming
3The *simplex* method uses a traditional, full-tableau implementation of
4Dantzig's simplex algorithm [1]_, [2]_ (*not* the Nelder-Mead simplex).
5This algorithm is included for backwards compatibility and educational
6purposes.
8 .. versionadded:: 0.15.0
10Warnings
11--------
13The simplex method may encounter numerical difficulties when pivot
14values are close to the specified tolerance. If encountered try
15remove any redundant constraints, change the pivot strategy to Bland's
16rule or increase the tolerance value.
18Alternatively, more robust methods maybe be used. See
19:ref:`'interior-point' <optimize.linprog-interior-point>` and
20:ref:`'revised simplex' <optimize.linprog-revised_simplex>`.
22References
23----------
24.. [1] Dantzig, George B., Linear programming and extensions. Rand
25 Corporation Research Study Princeton Univ. Press, Princeton, NJ,
26 1963
27.. [2] Hillier, S.H. and Lieberman, G.J. (1995), "Introduction to
28 Mathematical Programming", McGraw-Hill, Chapter 4.
29"""
31import numpy as np
32from warnings import warn
33from ._optimize import OptimizeResult, OptimizeWarning, _check_unknown_options
34from ._linprog_util import _postsolve
37def _pivot_col(T, tol=1e-9, bland=False):
38 """
39 Given a linear programming simplex tableau, determine the column
40 of the variable to enter the basis.
42 Parameters
43 ----------
44 T : 2-D array
45 A 2-D array representing the simplex tableau, T, corresponding to the
46 linear programming problem. It should have the form:
48 [[A[0, 0], A[0, 1], ..., A[0, n_total], b[0]],
49 [A[1, 0], A[1, 1], ..., A[1, n_total], b[1]],
50 .
51 .
52 .
53 [A[m, 0], A[m, 1], ..., A[m, n_total], b[m]],
54 [c[0], c[1], ..., c[n_total], 0]]
56 for a Phase 2 problem, or the form:
58 [[A[0, 0], A[0, 1], ..., A[0, n_total], b[0]],
59 [A[1, 0], A[1, 1], ..., A[1, n_total], b[1]],
60 .
61 .
62 .
63 [A[m, 0], A[m, 1], ..., A[m, n_total], b[m]],
64 [c[0], c[1], ..., c[n_total], 0],
65 [c'[0], c'[1], ..., c'[n_total], 0]]
67 for a Phase 1 problem (a problem in which a basic feasible solution is
68 sought prior to maximizing the actual objective. ``T`` is modified in
69 place by ``_solve_simplex``.
70 tol : float
71 Elements in the objective row larger than -tol will not be considered
72 for pivoting. Nominally this value is zero, but numerical issues
73 cause a tolerance about zero to be necessary.
74 bland : bool
75 If True, use Bland's rule for selection of the column (select the
76 first column with a negative coefficient in the objective row,
77 regardless of magnitude).
79 Returns
80 -------
81 status: bool
82 True if a suitable pivot column was found, otherwise False.
83 A return of False indicates that the linear programming simplex
84 algorithm is complete.
85 col: int
86 The index of the column of the pivot element.
87 If status is False, col will be returned as nan.
88 """
89 ma = np.ma.masked_where(T[-1, :-1] >= -tol, T[-1, :-1], copy=False)
90 if ma.count() == 0:
91 return False, np.nan
92 if bland:
93 # ma.mask is sometimes 0d
94 return True, np.nonzero(np.logical_not(np.atleast_1d(ma.mask)))[0][0]
95 return True, np.ma.nonzero(ma == ma.min())[0][0]
98def _pivot_row(T, basis, pivcol, phase, tol=1e-9, bland=False):
99 """
100 Given a linear programming simplex tableau, determine the row for the
101 pivot operation.
103 Parameters
104 ----------
105 T : 2-D array
106 A 2-D array representing the simplex tableau, T, corresponding to the
107 linear programming problem. It should have the form:
109 [[A[0, 0], A[0, 1], ..., A[0, n_total], b[0]],
110 [A[1, 0], A[1, 1], ..., A[1, n_total], b[1]],
111 .
112 .
113 .
114 [A[m, 0], A[m, 1], ..., A[m, n_total], b[m]],
115 [c[0], c[1], ..., c[n_total], 0]]
117 for a Phase 2 problem, or the form:
119 [[A[0, 0], A[0, 1], ..., A[0, n_total], b[0]],
120 [A[1, 0], A[1, 1], ..., A[1, n_total], b[1]],
121 .
122 .
123 .
124 [A[m, 0], A[m, 1], ..., A[m, n_total], b[m]],
125 [c[0], c[1], ..., c[n_total], 0],
126 [c'[0], c'[1], ..., c'[n_total], 0]]
128 for a Phase 1 problem (a Problem in which a basic feasible solution is
129 sought prior to maximizing the actual objective. ``T`` is modified in
130 place by ``_solve_simplex``.
131 basis : array
132 A list of the current basic variables.
133 pivcol : int
134 The index of the pivot column.
135 phase : int
136 The phase of the simplex algorithm (1 or 2).
137 tol : float
138 Elements in the pivot column smaller than tol will not be considered
139 for pivoting. Nominally this value is zero, but numerical issues
140 cause a tolerance about zero to be necessary.
141 bland : bool
142 If True, use Bland's rule for selection of the row (if more than one
143 row can be used, choose the one with the lowest variable index).
145 Returns
146 -------
147 status: bool
148 True if a suitable pivot row was found, otherwise False. A return
149 of False indicates that the linear programming problem is unbounded.
150 row: int
151 The index of the row of the pivot element. If status is False, row
152 will be returned as nan.
153 """
154 if phase == 1:
155 k = 2
156 else:
157 k = 1
158 ma = np.ma.masked_where(T[:-k, pivcol] <= tol, T[:-k, pivcol], copy=False)
159 if ma.count() == 0:
160 return False, np.nan
161 mb = np.ma.masked_where(T[:-k, pivcol] <= tol, T[:-k, -1], copy=False)
162 q = mb / ma
163 min_rows = np.ma.nonzero(q == q.min())[0]
164 if bland:
165 return True, min_rows[np.argmin(np.take(basis, min_rows))]
166 return True, min_rows[0]
169def _apply_pivot(T, basis, pivrow, pivcol, tol=1e-9):
170 """
171 Pivot the simplex tableau inplace on the element given by (pivrow, pivol).
172 The entering variable corresponds to the column given by pivcol forcing
173 the variable basis[pivrow] to leave the basis.
175 Parameters
176 ----------
177 T : 2-D array
178 A 2-D array representing the simplex tableau, T, corresponding to the
179 linear programming problem. It should have the form:
181 [[A[0, 0], A[0, 1], ..., A[0, n_total], b[0]],
182 [A[1, 0], A[1, 1], ..., A[1, n_total], b[1]],
183 .
184 .
185 .
186 [A[m, 0], A[m, 1], ..., A[m, n_total], b[m]],
187 [c[0], c[1], ..., c[n_total], 0]]
189 for a Phase 2 problem, or the form:
191 [[A[0, 0], A[0, 1], ..., A[0, n_total], b[0]],
192 [A[1, 0], A[1, 1], ..., A[1, n_total], b[1]],
193 .
194 .
195 .
196 [A[m, 0], A[m, 1], ..., A[m, n_total], b[m]],
197 [c[0], c[1], ..., c[n_total], 0],
198 [c'[0], c'[1], ..., c'[n_total], 0]]
200 for a Phase 1 problem (a problem in which a basic feasible solution is
201 sought prior to maximizing the actual objective. ``T`` is modified in
202 place by ``_solve_simplex``.
203 basis : 1-D array
204 An array of the indices of the basic variables, such that basis[i]
205 contains the column corresponding to the basic variable for row i.
206 Basis is modified in place by _apply_pivot.
207 pivrow : int
208 Row index of the pivot.
209 pivcol : int
210 Column index of the pivot.
211 """
212 basis[pivrow] = pivcol
213 pivval = T[pivrow, pivcol]
214 T[pivrow] = T[pivrow] / pivval
215 for irow in range(T.shape[0]):
216 if irow != pivrow:
217 T[irow] = T[irow] - T[pivrow] * T[irow, pivcol]
219 # The selected pivot should never lead to a pivot value less than the tol.
220 if np.isclose(pivval, tol, atol=0, rtol=1e4):
221 message = (
222 "The pivot operation produces a pivot value of:{0: .1e}, "
223 "which is only slightly greater than the specified "
224 "tolerance{1: .1e}. This may lead to issues regarding the "
225 "numerical stability of the simplex method. "
226 "Removing redundant constraints, changing the pivot strategy "
227 "via Bland's rule or increasing the tolerance may "
228 "help reduce the issue.".format(pivval, tol))
229 warn(message, OptimizeWarning, stacklevel=5)
232def _solve_simplex(T, n, basis, callback, postsolve_args,
233 maxiter=1000, tol=1e-9, phase=2, bland=False, nit0=0,
234 ):
235 """
236 Solve a linear programming problem in "standard form" using the Simplex
237 Method. Linear Programming is intended to solve the following problem form:
239 Minimize::
241 c @ x
243 Subject to::
245 A @ x == b
246 x >= 0
248 Parameters
249 ----------
250 T : 2-D array
251 A 2-D array representing the simplex tableau, T, corresponding to the
252 linear programming problem. It should have the form:
254 [[A[0, 0], A[0, 1], ..., A[0, n_total], b[0]],
255 [A[1, 0], A[1, 1], ..., A[1, n_total], b[1]],
256 .
257 .
258 .
259 [A[m, 0], A[m, 1], ..., A[m, n_total], b[m]],
260 [c[0], c[1], ..., c[n_total], 0]]
262 for a Phase 2 problem, or the form:
264 [[A[0, 0], A[0, 1], ..., A[0, n_total], b[0]],
265 [A[1, 0], A[1, 1], ..., A[1, n_total], b[1]],
266 .
267 .
268 .
269 [A[m, 0], A[m, 1], ..., A[m, n_total], b[m]],
270 [c[0], c[1], ..., c[n_total], 0],
271 [c'[0], c'[1], ..., c'[n_total], 0]]
273 for a Phase 1 problem (a problem in which a basic feasible solution is
274 sought prior to maximizing the actual objective. ``T`` is modified in
275 place by ``_solve_simplex``.
276 n : int
277 The number of true variables in the problem.
278 basis : 1-D array
279 An array of the indices of the basic variables, such that basis[i]
280 contains the column corresponding to the basic variable for row i.
281 Basis is modified in place by _solve_simplex
282 callback : callable, optional
283 If a callback function is provided, it will be called within each
284 iteration of the algorithm. The callback must accept a
285 `scipy.optimize.OptimizeResult` consisting of the following fields:
287 x : 1-D array
288 Current solution vector
289 fun : float
290 Current value of the objective function
291 success : bool
292 True only when a phase has completed successfully. This
293 will be False for most iterations.
294 slack : 1-D array
295 The values of the slack variables. Each slack variable
296 corresponds to an inequality constraint. If the slack is zero,
297 the corresponding constraint is active.
298 con : 1-D array
299 The (nominally zero) residuals of the equality constraints,
300 that is, ``b - A_eq @ x``
301 phase : int
302 The phase of the optimization being executed. In phase 1 a basic
303 feasible solution is sought and the T has an additional row
304 representing an alternate objective function.
305 status : int
306 An integer representing the exit status of the optimization::
308 0 : Optimization terminated successfully
309 1 : Iteration limit reached
310 2 : Problem appears to be infeasible
311 3 : Problem appears to be unbounded
312 4 : Serious numerical difficulties encountered
314 nit : int
315 The number of iterations performed.
316 message : str
317 A string descriptor of the exit status of the optimization.
318 postsolve_args : tuple
319 Data needed by _postsolve to convert the solution to the standard-form
320 problem into the solution to the original problem.
321 maxiter : int
322 The maximum number of iterations to perform before aborting the
323 optimization.
324 tol : float
325 The tolerance which determines when a solution is "close enough" to
326 zero in Phase 1 to be considered a basic feasible solution or close
327 enough to positive to serve as an optimal solution.
328 phase : int
329 The phase of the optimization being executed. In phase 1 a basic
330 feasible solution is sought and the T has an additional row
331 representing an alternate objective function.
332 bland : bool
333 If True, choose pivots using Bland's rule [3]_. In problems which
334 fail to converge due to cycling, using Bland's rule can provide
335 convergence at the expense of a less optimal path about the simplex.
336 nit0 : int
337 The initial iteration number used to keep an accurate iteration total
338 in a two-phase problem.
340 Returns
341 -------
342 nit : int
343 The number of iterations. Used to keep an accurate iteration total
344 in the two-phase problem.
345 status : int
346 An integer representing the exit status of the optimization::
348 0 : Optimization terminated successfully
349 1 : Iteration limit reached
350 2 : Problem appears to be infeasible
351 3 : Problem appears to be unbounded
352 4 : Serious numerical difficulties encountered
354 """
355 nit = nit0
356 status = 0
357 message = ''
358 complete = False
360 if phase == 1:
361 m = T.shape[1]-2
362 elif phase == 2:
363 m = T.shape[1]-1
364 else:
365 raise ValueError("Argument 'phase' to _solve_simplex must be 1 or 2")
367 if phase == 2:
368 # Check if any artificial variables are still in the basis.
369 # If yes, check if any coefficients from this row and a column
370 # corresponding to one of the non-artificial variable is non-zero.
371 # If found, pivot at this term. If not, start phase 2.
372 # Do this for all artificial variables in the basis.
373 # Ref: "An Introduction to Linear Programming and Game Theory"
374 # by Paul R. Thie, Gerard E. Keough, 3rd Ed,
375 # Chapter 3.7 Redundant Systems (pag 102)
376 for pivrow in [row for row in range(basis.size)
377 if basis[row] > T.shape[1] - 2]:
378 non_zero_row = [col for col in range(T.shape[1] - 1)
379 if abs(T[pivrow, col]) > tol]
380 if len(non_zero_row) > 0:
381 pivcol = non_zero_row[0]
382 _apply_pivot(T, basis, pivrow, pivcol, tol)
383 nit += 1
385 if len(basis[:m]) == 0:
386 solution = np.empty(T.shape[1] - 1, dtype=np.float64)
387 else:
388 solution = np.empty(max(T.shape[1] - 1, max(basis[:m]) + 1),
389 dtype=np.float64)
391 while not complete:
392 # Find the pivot column
393 pivcol_found, pivcol = _pivot_col(T, tol, bland)
394 if not pivcol_found:
395 pivcol = np.nan
396 pivrow = np.nan
397 status = 0
398 complete = True
399 else:
400 # Find the pivot row
401 pivrow_found, pivrow = _pivot_row(T, basis, pivcol, phase, tol, bland)
402 if not pivrow_found:
403 status = 3
404 complete = True
406 if callback is not None:
407 solution[:] = 0
408 solution[basis[:n]] = T[:n, -1]
409 x = solution[:m]
410 x, fun, slack, con = _postsolve(
411 x, postsolve_args
412 )
413 res = OptimizeResult({
414 'x': x,
415 'fun': fun,
416 'slack': slack,
417 'con': con,
418 'status': status,
419 'message': message,
420 'nit': nit,
421 'success': status == 0 and complete,
422 'phase': phase,
423 'complete': complete,
424 })
425 callback(res)
427 if not complete:
428 if nit >= maxiter:
429 # Iteration limit exceeded
430 status = 1
431 complete = True
432 else:
433 _apply_pivot(T, basis, pivrow, pivcol, tol)
434 nit += 1
435 return nit, status
438def _linprog_simplex(c, c0, A, b, callback, postsolve_args,
439 maxiter=1000, tol=1e-9, disp=False, bland=False,
440 **unknown_options):
441 """
442 Minimize a linear objective function subject to linear equality and
443 non-negativity constraints using the two phase simplex method.
444 Linear programming is intended to solve problems of the following form:
446 Minimize::
448 c @ x
450 Subject to::
452 A @ x == b
453 x >= 0
455 User-facing documentation is in _linprog_doc.py.
457 Parameters
458 ----------
459 c : 1-D array
460 Coefficients of the linear objective function to be minimized.
461 c0 : float
462 Constant term in objective function due to fixed (and eliminated)
463 variables. (Purely for display.)
464 A : 2-D array
465 2-D array such that ``A @ x``, gives the values of the equality
466 constraints at ``x``.
467 b : 1-D array
468 1-D array of values representing the right hand side of each equality
469 constraint (row) in ``A``.
470 callback : callable, optional
471 If a callback function is provided, it will be called within each
472 iteration of the algorithm. The callback function must accept a single
473 `scipy.optimize.OptimizeResult` consisting of the following fields:
475 x : 1-D array
476 Current solution vector
477 fun : float
478 Current value of the objective function
479 success : bool
480 True when an algorithm has completed successfully.
481 slack : 1-D array
482 The values of the slack variables. Each slack variable
483 corresponds to an inequality constraint. If the slack is zero,
484 the corresponding constraint is active.
485 con : 1-D array
486 The (nominally zero) residuals of the equality constraints,
487 that is, ``b - A_eq @ x``
488 phase : int
489 The phase of the algorithm being executed.
490 status : int
491 An integer representing the status of the optimization::
493 0 : Algorithm proceeding nominally
494 1 : Iteration limit reached
495 2 : Problem appears to be infeasible
496 3 : Problem appears to be unbounded
497 4 : Serious numerical difficulties encountered
498 nit : int
499 The number of iterations performed.
500 message : str
501 A string descriptor of the exit status of the optimization.
502 postsolve_args : tuple
503 Data needed by _postsolve to convert the solution to the standard-form
504 problem into the solution to the original problem.
506 Options
507 -------
508 maxiter : int
509 The maximum number of iterations to perform.
510 disp : bool
511 If True, print exit status message to sys.stdout
512 tol : float
513 The tolerance which determines when a solution is "close enough" to
514 zero in Phase 1 to be considered a basic feasible solution or close
515 enough to positive to serve as an optimal solution.
516 bland : bool
517 If True, use Bland's anti-cycling rule [3]_ to choose pivots to
518 prevent cycling. If False, choose pivots which should lead to a
519 converged solution more quickly. The latter method is subject to
520 cycling (non-convergence) in rare instances.
521 unknown_options : dict
522 Optional arguments not used by this particular solver. If
523 `unknown_options` is non-empty a warning is issued listing all
524 unused options.
526 Returns
527 -------
528 x : 1-D array
529 Solution vector.
530 status : int
531 An integer representing the exit status of the optimization::
533 0 : Optimization terminated successfully
534 1 : Iteration limit reached
535 2 : Problem appears to be infeasible
536 3 : Problem appears to be unbounded
537 4 : Serious numerical difficulties encountered
539 message : str
540 A string descriptor of the exit status of the optimization.
541 iteration : int
542 The number of iterations taken to solve the problem.
544 References
545 ----------
546 .. [1] Dantzig, George B., Linear programming and extensions. Rand
547 Corporation Research Study Princeton Univ. Press, Princeton, NJ,
548 1963
549 .. [2] Hillier, S.H. and Lieberman, G.J. (1995), "Introduction to
550 Mathematical Programming", McGraw-Hill, Chapter 4.
551 .. [3] Bland, Robert G. New finite pivoting rules for the simplex method.
552 Mathematics of Operations Research (2), 1977: pp. 103-107.
555 Notes
556 -----
557 The expected problem formulation differs between the top level ``linprog``
558 module and the method specific solvers. The method specific solvers expect a
559 problem in standard form:
561 Minimize::
563 c @ x
565 Subject to::
567 A @ x == b
568 x >= 0
570 Whereas the top level ``linprog`` module expects a problem of form:
572 Minimize::
574 c @ x
576 Subject to::
578 A_ub @ x <= b_ub
579 A_eq @ x == b_eq
580 lb <= x <= ub
582 where ``lb = 0`` and ``ub = None`` unless set in ``bounds``.
584 The original problem contains equality, upper-bound and variable constraints
585 whereas the method specific solver requires equality constraints and
586 variable non-negativity.
588 ``linprog`` module converts the original problem to standard form by
589 converting the simple bounds to upper bound constraints, introducing
590 non-negative slack variables for inequality constraints, and expressing
591 unbounded variables as the difference between two non-negative variables.
592 """
593 _check_unknown_options(unknown_options)
595 status = 0
596 messages = {0: "Optimization terminated successfully.",
597 1: "Iteration limit reached.",
598 2: "Optimization failed. Unable to find a feasible"
599 " starting point.",
600 3: "Optimization failed. The problem appears to be unbounded.",
601 4: "Optimization failed. Singular matrix encountered."}
603 n, m = A.shape
605 # All constraints must have b >= 0.
606 is_negative_constraint = np.less(b, 0)
607 A[is_negative_constraint] *= -1
608 b[is_negative_constraint] *= -1
610 # As all constraints are equality constraints the artificial variables
611 # will also be basic variables.
612 av = np.arange(n) + m
613 basis = av.copy()
615 # Format the phase one tableau by adding artificial variables and stacking
616 # the constraints, the objective row and pseudo-objective row.
617 row_constraints = np.hstack((A, np.eye(n), b[:, np.newaxis]))
618 row_objective = np.hstack((c, np.zeros(n), c0))
619 row_pseudo_objective = -row_constraints.sum(axis=0)
620 row_pseudo_objective[av] = 0
621 T = np.vstack((row_constraints, row_objective, row_pseudo_objective))
623 nit1, status = _solve_simplex(T, n, basis, callback=callback,
624 postsolve_args=postsolve_args,
625 maxiter=maxiter, tol=tol, phase=1,
626 bland=bland
627 )
628 # if pseudo objective is zero, remove the last row from the tableau and
629 # proceed to phase 2
630 nit2 = nit1
631 if abs(T[-1, -1]) < tol:
632 # Remove the pseudo-objective row from the tableau
633 T = T[:-1, :]
634 # Remove the artificial variable columns from the tableau
635 T = np.delete(T, av, 1)
636 else:
637 # Failure to find a feasible starting point
638 status = 2
639 messages[status] = (
640 "Phase 1 of the simplex method failed to find a feasible "
641 "solution. The pseudo-objective function evaluates to {0:.1e} "
642 "which exceeds the required tolerance of {1} for a solution to be "
643 "considered 'close enough' to zero to be a basic solution. "
644 "Consider increasing the tolerance to be greater than {0:.1e}. "
645 "If this tolerance is unacceptably large the problem may be "
646 "infeasible.".format(abs(T[-1, -1]), tol)
647 )
649 if status == 0:
650 # Phase 2
651 nit2, status = _solve_simplex(T, n, basis, callback=callback,
652 postsolve_args=postsolve_args,
653 maxiter=maxiter, tol=tol, phase=2,
654 bland=bland, nit0=nit1
655 )
657 solution = np.zeros(n + m)
658 solution[basis[:n]] = T[:n, -1]
659 x = solution[:m]
661 return x, status, messages[status], int(nit2)