_linprog.py 19 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490
  1. """
  2. A top-level linear programming interface. Currently this interface solves
  3. linear programming problems via the Simplex and Interior-Point methods.
  4. .. versionadded:: 0.15.0
  5. Functions
  6. ---------
  7. .. autosummary::
  8. :toctree: generated/
  9. linprog
  10. linprog_verbose_callback
  11. linprog_terse_callback
  12. """
  13. from __future__ import division, print_function, absolute_import
  14. import numpy as np
  15. from .optimize import OptimizeResult
  16. from ._linprog_ip import _linprog_ip
  17. from ._linprog_simplex import _linprog_simplex
  18. from ._linprog_util import (
  19. _parse_linprog, _presolve, _get_Abc, _postprocess
  20. )
  21. __all__ = ['linprog', 'linprog_verbose_callback', 'linprog_terse_callback']
  22. __docformat__ = "restructuredtext en"
  23. def linprog_verbose_callback(res):
  24. """
  25. A sample callback function demonstrating the linprog callback interface.
  26. This callback produces detailed output to sys.stdout before each iteration
  27. and after the final iteration of the simplex algorithm.
  28. Parameters
  29. ----------
  30. res : A `scipy.optimize.OptimizeResult` consisting of the following fields:
  31. x : 1D array
  32. The independent variable vector which optimizes the linear
  33. programming problem.
  34. fun : float
  35. Value of the objective function.
  36. success : bool
  37. True if the algorithm succeeded in finding an optimal solution.
  38. slack : 1D array
  39. The values of the slack variables. Each slack variable corresponds
  40. to an inequality constraint. If the slack is zero, then the
  41. corresponding constraint is active.
  42. con : 1D array
  43. The (nominally zero) residuals of the equality constraints, that is,
  44. ``b - A_eq @ x``
  45. phase : int
  46. The phase of the optimization being executed. In phase 1 a basic
  47. feasible solution is sought and the T has an additional row
  48. representing an alternate objective function.
  49. status : int
  50. An integer representing the exit status of the optimization::
  51. 0 : Optimization terminated successfully
  52. 1 : Iteration limit reached
  53. 2 : Problem appears to be infeasible
  54. 3 : Problem appears to be unbounded
  55. 4 : Serious numerical difficulties encountered
  56. nit : int
  57. The number of iterations performed.
  58. message : str
  59. A string descriptor of the exit status of the optimization.
  60. """
  61. x = res['x']
  62. fun = res['fun']
  63. success = res['success']
  64. phase = res['phase']
  65. status = res['status']
  66. nit = res['nit']
  67. message = res['message']
  68. complete = res['complete']
  69. saved_printoptions = np.get_printoptions()
  70. np.set_printoptions(linewidth=500,
  71. formatter={'float': lambda x: "{0: 12.4f}".format(x)})
  72. if status:
  73. print('--------- Simplex Early Exit -------\n'.format(nit))
  74. print('The simplex method exited early with status {0:d}'.format(status))
  75. print(message)
  76. elif complete:
  77. print('--------- Simplex Complete --------\n')
  78. print('Iterations required: {}'.format(nit))
  79. else:
  80. print('--------- Iteration {0:d} ---------\n'.format(nit))
  81. if nit > 0:
  82. if phase == 1:
  83. print('Current Pseudo-Objective Value:')
  84. else:
  85. print('Current Objective Value:')
  86. print('f = ', fun)
  87. print()
  88. print('Current Solution Vector:')
  89. print('x = ', x)
  90. print()
  91. np.set_printoptions(**saved_printoptions)
  92. def linprog_terse_callback(res):
  93. """
  94. A sample callback function demonstrating the linprog callback interface.
  95. This callback produces brief output to sys.stdout before each iteration
  96. and after the final iteration of the simplex algorithm.
  97. Parameters
  98. ----------
  99. res : A `scipy.optimize.OptimizeResult` consisting of the following fields:
  100. x : 1D array
  101. The independent variable vector which optimizes the linear
  102. programming problem.
  103. fun : float
  104. Value of the objective function.
  105. success : bool
  106. True if the algorithm succeeded in finding an optimal solution.
  107. slack : 1D array
  108. The values of the slack variables. Each slack variable corresponds
  109. to an inequality constraint. If the slack is zero, then the
  110. corresponding constraint is active.
  111. con : 1D array
  112. The (nominally zero) residuals of the equality constraints, that is,
  113. ``b - A_eq @ x``
  114. phase : int
  115. The phase of the optimization being executed. In phase 1 a basic
  116. feasible solution is sought and the T has an additional row
  117. representing an alternate objective function.
  118. status : int
  119. An integer representing the exit status of the optimization::
  120. 0 : Optimization terminated successfully
  121. 1 : Iteration limit reached
  122. 2 : Problem appears to be infeasible
  123. 3 : Problem appears to be unbounded
  124. 4 : Serious numerical difficulties encountered
  125. nit : int
  126. The number of iterations performed.
  127. message : str
  128. A string descriptor of the exit status of the optimization.
  129. """
  130. nit = res['nit']
  131. x = res['x']
  132. if nit == 0:
  133. print("Iter: X:")
  134. print("{0: <5d} ".format(nit), end="")
  135. print(x)
  136. def linprog(c, A_ub=None, b_ub=None, A_eq=None, b_eq=None,
  137. bounds=None, method='simplex', callback=None,
  138. options=None):
  139. """
  140. Minimize a linear objective function subject to linear
  141. equality and inequality constraints. Linear Programming is intended to
  142. solve the following problem form:
  143. Minimize::
  144. c @ x
  145. Subject to::
  146. A_ub @ x <= b_ub
  147. A_eq @ x == b_eq
  148. lb <= x <= ub
  149. where ``lb = 0`` and ``ub = None`` unless set in ``bounds``.
  150. Parameters
  151. ----------
  152. c : 1D array
  153. Coefficients of the linear objective function to be minimized.
  154. A_ub : 2D array, optional
  155. 2D array such that ``A_ub @ x`` gives the values of the upper-bound
  156. inequality constraints at ``x``.
  157. b_ub : 1D array, optional
  158. 1D array of values representing the upper-bound of each inequality
  159. constraint (row) in ``A_ub``.
  160. A_eq : 2D, optional
  161. 2D array such that ``A_eq @ x`` gives the values of the equality
  162. constraints at ``x``.
  163. b_eq : 1D array, optional
  164. 1D array of values representing the RHS of each equality constraint
  165. (row) in ``A_eq``.
  166. bounds : sequence, optional
  167. ``(min, max)`` pairs for each element in ``x``, defining
  168. the bounds on that parameter. Use None for one of ``min`` or
  169. ``max`` when there is no bound in that direction. By default
  170. bounds are ``(0, None)`` (non-negative).
  171. If a sequence containing a single tuple is provided, then ``min`` and
  172. ``max`` will be applied to all variables in the problem.
  173. method : str, optional
  174. Type of solver. :ref:`'simplex' <optimize.linprog-simplex>`
  175. and :ref:`'interior-point' <optimize.linprog-interior-point>`
  176. are supported.
  177. callback : callable, optional (simplex only)
  178. If a callback function is provided, it will be called within each
  179. iteration of the simplex algorithm. The callback must require a
  180. `scipy.optimize.OptimizeResult` consisting of the following fields:
  181. x : 1D array
  182. The independent variable vector which optimizes the linear
  183. programming problem.
  184. fun : float
  185. Value of the objective function.
  186. success : bool
  187. True if the algorithm succeeded in finding an optimal solution.
  188. slack : 1D array
  189. The values of the slack variables. Each slack variable
  190. corresponds to an inequality constraint. If the slack is zero,
  191. the corresponding constraint is active.
  192. con : 1D array
  193. The (nominally zero) residuals of the equality constraints
  194. that is, ``b - A_eq @ x``
  195. phase : int
  196. The phase of the optimization being executed. In phase 1 a basic
  197. feasible solution is sought and the T has an additional row
  198. representing an alternate objective function.
  199. status : int
  200. An integer representing the exit status of the optimization::
  201. 0 : Optimization terminated successfully
  202. 1 : Iteration limit reached
  203. 2 : Problem appears to be infeasible
  204. 3 : Problem appears to be unbounded
  205. 4 : Serious numerical difficulties encountered
  206. nit : int
  207. The number of iterations performed.
  208. message : str
  209. A string descriptor of the exit status of the optimization.
  210. options : dict, optional
  211. A dictionary of solver options. All methods accept the following
  212. generic options:
  213. maxiter : int
  214. Maximum number of iterations to perform.
  215. disp : bool
  216. Set to True to print convergence messages.
  217. For method-specific options, see :func:`show_options('linprog')`.
  218. Returns
  219. -------
  220. res : OptimizeResult
  221. A :class:`scipy.optimize.OptimizeResult` consisting of the fields:
  222. x : 1D array
  223. The independent variable vector which optimizes the linear
  224. programming problem.
  225. fun : float
  226. Value of the objective function.
  227. slack : 1D array
  228. The values of the slack variables. Each slack variable
  229. corresponds to an inequality constraint. If the slack is zero,
  230. then the corresponding constraint is active.
  231. con : 1D array
  232. The (nominally zero) residuals of the equality constraints,
  233. that is, ``b - A_eq @ x``
  234. success : bool
  235. Returns True if the algorithm succeeded in finding an optimal
  236. solution.
  237. status : int
  238. An integer representing the exit status of the optimization::
  239. 0 : Optimization terminated successfully
  240. 1 : Iteration limit reached
  241. 2 : Problem appears to be infeasible
  242. 3 : Problem appears to be unbounded
  243. 4 : Serious numerical difficulties encountered
  244. nit : int
  245. The number of iterations performed.
  246. message : str
  247. A string descriptor of the exit status of the optimization.
  248. See Also
  249. --------
  250. show_options : Additional options accepted by the solvers
  251. Notes
  252. -----
  253. This section describes the available solvers that can be selected by the
  254. 'method' parameter. The default method
  255. is :ref:`Simplex <optimize.linprog-simplex>`.
  256. :ref:`Interior point <optimize.linprog-interior-point>` is also available.
  257. Method *simplex* uses the simplex algorithm (as it relates to linear
  258. programming, NOT the Nelder-Mead simplex) [1]_, [2]_. This algorithm
  259. should be reasonably reliable and fast for small problems.
  260. .. versionadded:: 0.15.0
  261. Method *interior-point* uses the primal-dual path following algorithm
  262. as outlined in [4]_. This algorithm is intended to provide a faster
  263. and more reliable alternative to *simplex*, especially for large,
  264. sparse problems. Note, however, that the solution returned may be slightly
  265. less accurate than that of the simplex method and may not correspond with a
  266. vertex of the polytope defined by the constraints.
  267. Before applying either method a presolve procedure based on [8]_ attempts to
  268. identify trivial infeasibilities, trivial unboundedness, and potential
  269. problem simplifications. Specifically, it checks for:
  270. - rows of zeros in ``A_eq`` or ``A_ub``, representing trivial constraints;
  271. - columns of zeros in ``A_eq`` `and` ``A_ub``, representing unconstrained
  272. variables;
  273. - column singletons in ``A_eq``, representing fixed variables; and
  274. - column singletons in ``A_ub``, representing simple bounds.
  275. If presolve reveals that the problem is unbounded (e.g. an unconstrained
  276. and unbounded variable has negative cost) or infeasible (e.g. a row of
  277. zeros in ``A_eq`` corresponds with a nonzero in ``b_eq``), the solver
  278. terminates with the appropriate status code. Note that presolve terminates
  279. as soon as any sign of unboundedness is detected; consequently, a problem
  280. may be reported as unbounded when in reality the problem is infeasible
  281. (but infeasibility has not been detected yet). Therefore, if the output
  282. message states that unboundedness is detected in presolve and it is
  283. necessary to know whether the problem is actually infeasible, set option
  284. ``presolve=False``.
  285. If neither infeasibility nor unboundedness are detected in a single pass
  286. of the presolve check, bounds are tightened where possible and fixed
  287. variables are removed from the problem. Then, linearly dependent rows
  288. of the ``A_eq`` matrix are removed, (unless they represent an
  289. infeasibility) to avoid numerical difficulties in the primary solve
  290. routine. Note that rows that are nearly linearly dependent (within a
  291. prescribed tolerance) may also be removed, which can change the optimal
  292. solution in rare cases. If this is a concern, eliminate redundancy from
  293. your problem formulation and run with option ``rr=False`` or
  294. ``presolve=False``.
  295. Several potential improvements can be made here: additional presolve
  296. checks outlined in [8]_ should be implemented, the presolve routine should
  297. be run multiple times (until no further simplifications can be made), and
  298. more of the efficiency improvements from [5]_ should be implemented in the
  299. redundancy removal routines.
  300. After presolve, the problem is transformed to standard form by converting
  301. the (tightened) simple bounds to upper bound constraints, introducing
  302. non-negative slack variables for inequality constraints, and expressing
  303. unbounded variables as the difference between two non-negative variables.
  304. References
  305. ----------
  306. .. [1] Dantzig, George B., Linear programming and extensions. Rand
  307. Corporation Research Study Princeton Univ. Press, Princeton, NJ,
  308. 1963
  309. .. [2] Hillier, S.H. and Lieberman, G.J. (1995), "Introduction to
  310. Mathematical Programming", McGraw-Hill, Chapter 4.
  311. .. [3] Bland, Robert G. New finite pivoting rules for the simplex method.
  312. Mathematics of Operations Research (2), 1977: pp. 103-107.
  313. .. [4] Andersen, Erling D., and Knud D. Andersen. "The MOSEK interior point
  314. optimizer for linear programming: an implementation of the
  315. homogeneous algorithm." High performance optimization. Springer US,
  316. 2000. 197-232.
  317. .. [5] Andersen, Erling D. "Finding all linearly dependent rows in
  318. large-scale linear programming." Optimization Methods and Software
  319. 6.3 (1995): 219-227.
  320. .. [6] Freund, Robert M. "Primal-Dual Interior-Point Methods for Linear
  321. Programming based on Newton's Method." Unpublished Course Notes,
  322. March 2004. Available 2/25/2017 at
  323. https://ocw.mit.edu/courses/sloan-school-of-management/15-084j-nonlinear-programming-spring-2004/lecture-notes/lec14_int_pt_mthd.pdf
  324. .. [7] Fourer, Robert. "Solving Linear Programs by Interior-Point Methods."
  325. Unpublished Course Notes, August 26, 2005. Available 2/25/2017 at
  326. http://www.4er.org/CourseNotes/Book%20B/B-III.pdf
  327. .. [8] Andersen, Erling D., and Knud D. Andersen. "Presolving in linear
  328. programming." Mathematical Programming 71.2 (1995): 221-245.
  329. .. [9] Bertsimas, Dimitris, and J. Tsitsiklis. "Introduction to linear
  330. programming." Athena Scientific 1 (1997): 997.
  331. .. [10] Andersen, Erling D., et al. Implementation of interior point
  332. methods for large scale linear programming. HEC/Universite de
  333. Geneve, 1996.
  334. Examples
  335. --------
  336. Consider the following problem:
  337. Minimize::
  338. f = -1x[0] + 4x[1]
  339. Subject to::
  340. -3x[0] + 1x[1] <= 6
  341. 1x[0] + 2x[1] <= 4
  342. x[1] >= -3
  343. -inf <= x[0] <= inf
  344. This problem deviates from the standard linear programming problem.
  345. In standard form, linear programming problems assume the variables x are
  346. non-negative. Since the problem variables don't have the standard bounds of
  347. ``(0, None)``, the variable bounds must be set using ``bounds`` explicitly.
  348. There are two upper-bound constraints, which can be expressed as
  349. dot(A_ub, x) <= b_ub
  350. The input for this problem is as follows:
  351. >>> c = [-1, 4]
  352. >>> A = [[-3, 1], [1, 2]]
  353. >>> b = [6, 4]
  354. >>> x0_bounds = (None, None)
  355. >>> x1_bounds = (-3, None)
  356. >>> from scipy.optimize import linprog
  357. >>> res = linprog(c, A_ub=A, b_ub=b, bounds=(x0_bounds, x1_bounds),
  358. ... options={"disp": True})
  359. Optimization terminated successfully.
  360. Current function value: -22.000000
  361. Iterations: 5 # may vary
  362. >>> print(res)
  363. con: array([], dtype=float64)
  364. fun: -22.0
  365. message: 'Optimization terminated successfully.'
  366. nit: 5 # may vary
  367. slack: array([39., 0.]) # may vary
  368. status: 0
  369. success: True
  370. x: array([10., -3.])
  371. """
  372. meth = method.lower()
  373. default_tol = 1e-12 if meth == 'simplex' else 1e-9
  374. c, A_ub, b_ub, A_eq, b_eq, bounds, solver_options = _parse_linprog(
  375. c, A_ub, b_ub, A_eq, b_eq, bounds, options)
  376. tol = solver_options.get('tol', default_tol)
  377. iteration = 0
  378. complete = False # will become True if solved in presolve
  379. undo = []
  380. # Keep the original arrays to calculate slack/residuals for original
  381. # problem.
  382. c_o, A_ub_o, b_ub_o, A_eq_o, b_eq_o = c.copy(
  383. ), A_ub.copy(), b_ub.copy(), A_eq.copy(), b_eq.copy()
  384. # Solve trivial problem, eliminate variables, tighten bounds, etc...
  385. c0 = 0 # we might get a constant term in the objective
  386. if solver_options.pop('presolve', True):
  387. rr = solver_options.pop('rr', True)
  388. (c, c0, A_ub, b_ub, A_eq, b_eq, bounds, x, undo, complete, status,
  389. message) = _presolve(c, A_ub, b_ub, A_eq, b_eq, bounds, rr, tol)
  390. if not complete:
  391. A, b, c, c0 = _get_Abc(c, c0, A_ub, b_ub, A_eq, b_eq, bounds, undo)
  392. T_o = (c_o, A_ub_o, b_ub_o, A_eq_o, b_eq_o, bounds, undo)
  393. if meth == 'simplex':
  394. x, status, message, iteration = _linprog_simplex(
  395. c, c0=c0, A=A, b=b, callback=callback, _T_o=T_o, **solver_options)
  396. elif meth == 'interior-point':
  397. x, status, message, iteration = _linprog_ip(
  398. c, c0=c0, A=A, b=b, callback=callback, **solver_options)
  399. else:
  400. raise ValueError('Unknown solver %s' % method)
  401. # Eliminate artificial variables, re-introduce presolved variables, etc...
  402. # need modified bounds here to translate variables appropriately
  403. disp = solver_options.get('disp', False)
  404. x, fun, slack, con, status, message = _postprocess(
  405. x, c_o, A_ub_o, b_ub_o, A_eq_o, b_eq_o, bounds,
  406. complete, undo, status, message, tol, iteration, disp)
  407. sol = {
  408. 'x': x,
  409. 'fun': fun,
  410. 'slack': slack,
  411. 'con': con,
  412. 'status': status,
  413. 'message': message,
  414. 'nit': iteration,
  415. 'success': status == 0}
  416. return OptimizeResult(sol)