Coverage for /home/martinb/.local/share/virtualenvs/camcops/lib/python3.6/site-packages/scipy/optimize/_linprog_simplex.py : 8%

Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
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.zeros(T.shape[1] - 1, dtype=np.float64)
387 else:
388 solution = np.zeros(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, tol=tol
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 Parameters
456 ----------
457 c : 1-D array
458 Coefficients of the linear objective function to be minimized.
459 c0 : float
460 Constant term in objective function due to fixed (and eliminated)
461 variables. (Purely for display.)
462 A : 2-D array
463 2-D array such that ``A @ x``, gives the values of the equality
464 constraints at ``x``.
465 b : 1-D array
466 1-D array of values representing the right hand side of each equality
467 constraint (row) in ``A``.
468 callback : callable, optional
469 If a callback function is provided, it will be called within each
470 iteration of the algorithm. The callback function must accept a single
471 `scipy.optimize.OptimizeResult` consisting of the following fields:
473 x : 1-D array
474 Current solution vector
475 fun : float
476 Current value of the objective function
477 success : bool
478 True when an algorithm has completed successfully.
479 slack : 1-D array
480 The values of the slack variables. Each slack variable
481 corresponds to an inequality constraint. If the slack is zero,
482 the corresponding constraint is active.
483 con : 1-D array
484 The (nominally zero) residuals of the equality constraints,
485 that is, ``b - A_eq @ x``
486 phase : int
487 The phase of the algorithm being executed.
488 status : int
489 An integer representing the status of the optimization::
491 0 : Algorithm proceeding nominally
492 1 : Iteration limit reached
493 2 : Problem appears to be infeasible
494 3 : Problem appears to be unbounded
495 4 : Serious numerical difficulties encountered
496 nit : int
497 The number of iterations performed.
498 message : str
499 A string descriptor of the exit status of the optimization.
500 postsolve_args : tuple
501 Data needed by _postsolve to convert the solution to the standard-form
502 problem into the solution to the original problem.
504 Options
505 -------
506 maxiter : int
507 The maximum number of iterations to perform.
508 disp : bool
509 If True, print exit status message to sys.stdout
510 tol : float
511 The tolerance which determines when a solution is "close enough" to
512 zero in Phase 1 to be considered a basic feasible solution or close
513 enough to positive to serve as an optimal solution.
514 bland : bool
515 If True, use Bland's anti-cycling rule [3]_ to choose pivots to
516 prevent cycling. If False, choose pivots which should lead to a
517 converged solution more quickly. The latter method is subject to
518 cycling (non-convergence) in rare instances.
519 unknown_options : dict
520 Optional arguments not used by this particular solver. If
521 `unknown_options` is non-empty a warning is issued listing all
522 unused options.
524 Returns
525 -------
526 x : 1-D array
527 Solution vector.
528 status : int
529 An integer representing the exit status of the optimization::
531 0 : Optimization terminated successfully
532 1 : Iteration limit reached
533 2 : Problem appears to be infeasible
534 3 : Problem appears to be unbounded
535 4 : Serious numerical difficulties encountered
537 message : str
538 A string descriptor of the exit status of the optimization.
539 iteration : int
540 The number of iterations taken to solve the problem.
542 References
543 ----------
544 .. [1] Dantzig, George B., Linear programming and extensions. Rand
545 Corporation Research Study Princeton Univ. Press, Princeton, NJ,
546 1963
547 .. [2] Hillier, S.H. and Lieberman, G.J. (1995), "Introduction to
548 Mathematical Programming", McGraw-Hill, Chapter 4.
549 .. [3] Bland, Robert G. New finite pivoting rules for the simplex method.
550 Mathematics of Operations Research (2), 1977: pp. 103-107.
553 Notes
554 -----
555 The expected problem formulation differs between the top level ``linprog``
556 module and the method specific solvers. The method specific solvers expect a
557 problem in standard form:
559 Minimize::
561 c @ x
563 Subject to::
565 A @ x == b
566 x >= 0
568 Whereas the top level ``linprog`` module expects a problem of form:
570 Minimize::
572 c @ x
574 Subject to::
576 A_ub @ x <= b_ub
577 A_eq @ x == b_eq
578 lb <= x <= ub
580 where ``lb = 0`` and ``ub = None`` unless set in ``bounds``.
582 The original problem contains equality, upper-bound and variable constraints
583 whereas the method specific solver requires equality constraints and
584 variable non-negativity.
586 ``linprog`` module converts the original problem to standard form by
587 converting the simple bounds to upper bound constraints, introducing
588 non-negative slack variables for inequality constraints, and expressing
589 unbounded variables as the difference between two non-negative variables.
590 """
591 _check_unknown_options(unknown_options)
593 status = 0
594 messages = {0: "Optimization terminated successfully.",
595 1: "Iteration limit reached.",
596 2: "Optimization failed. Unable to find a feasible"
597 " starting point.",
598 3: "Optimization failed. The problem appears to be unbounded.",
599 4: "Optimization failed. Singular matrix encountered."}
601 n, m = A.shape
603 # All constraints must have b >= 0.
604 is_negative_constraint = np.less(b, 0)
605 A[is_negative_constraint] *= -1
606 b[is_negative_constraint] *= -1
608 # As all constraints are equality constraints the artificial variables
609 # will also be basic variables.
610 av = np.arange(n) + m
611 basis = av.copy()
613 # Format the phase one tableau by adding artificial variables and stacking
614 # the constraints, the objective row and pseudo-objective row.
615 row_constraints = np.hstack((A, np.eye(n), b[:, np.newaxis]))
616 row_objective = np.hstack((c, np.zeros(n), c0))
617 row_pseudo_objective = -row_constraints.sum(axis=0)
618 row_pseudo_objective[av] = 0
619 T = np.vstack((row_constraints, row_objective, row_pseudo_objective))
621 nit1, status = _solve_simplex(T, n, basis, callback=callback,
622 postsolve_args=postsolve_args,
623 maxiter=maxiter, tol=tol, phase=1,
624 bland=bland
625 )
626 # if pseudo objective is zero, remove the last row from the tableau and
627 # proceed to phase 2
628 nit2 = nit1
629 if abs(T[-1, -1]) < tol:
630 # Remove the pseudo-objective row from the tableau
631 T = T[:-1, :]
632 # Remove the artificial variable columns from the tableau
633 T = np.delete(T, av, 1)
634 else:
635 # Failure to find a feasible starting point
636 status = 2
637 messages[status] = (
638 "Phase 1 of the simplex method failed to find a feasible "
639 "solution. The pseudo-objective function evaluates to {0:.1e} "
640 "which exceeds the required tolerance of {1} for a solution to be "
641 "considered 'close enough' to zero to be a basic solution. "
642 "Consider increasing the tolerance to be greater than {0:.1e}. "
643 "If this tolerance is unacceptably large the problem may be "
644 "infeasible.".format(abs(T[-1, -1]), tol)
645 )
647 if status == 0:
648 # Phase 2
649 nit2, status = _solve_simplex(T, n, basis, callback=callback,
650 postsolve_args=postsolve_args,
651 maxiter=maxiter, tol=tol, phase=2,
652 bland=bland, nit0=nit1
653 )
655 solution = np.zeros(n + m)
656 solution[basis[:n]] = T[:n, -1]
657 x = solution[:m]
659 return x, status, messages[status], int(nit2)