Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-packages/scipy/optimize/_milp.py: 9%
106 statements
« prev ^ index » next coverage.py v7.3.2, created at 2023-12-12 06:31 +0000
1import warnings
2import numpy as np
3from scipy.sparse import csc_array, vstack
4from ._highs._highs_wrapper import _highs_wrapper # type: ignore[import]
5from ._constraints import LinearConstraint, Bounds
6from ._optimize import OptimizeResult
7from ._linprog_highs import _highs_to_scipy_status_message
10def _constraints_to_components(constraints):
11 """
12 Convert sequence of constraints to a single set of components A, b_l, b_u.
14 `constraints` could be
16 1. A LinearConstraint
17 2. A tuple representing a LinearConstraint
18 3. An invalid object
19 4. A sequence of composed entirely of objects of type 1/2
20 5. A sequence containing at least one object of type 3
22 We want to accept 1, 2, and 4 and reject 3 and 5.
23 """
24 message = ("`constraints` (or each element within `constraints`) must be "
25 "convertible into an instance of "
26 "`scipy.optimize.LinearConstraint`.")
27 As = []
28 b_ls = []
29 b_us = []
31 # Accept case 1 by standardizing as case 4
32 if isinstance(constraints, LinearConstraint):
33 constraints = [constraints]
34 else:
35 # Reject case 3
36 try:
37 iter(constraints)
38 except TypeError as exc:
39 raise ValueError(message) from exc
41 # Accept case 2 by standardizing as case 4
42 if len(constraints) == 3:
43 # argument could be a single tuple representing a LinearConstraint
44 try:
45 constraints = [LinearConstraint(*constraints)]
46 except (TypeError, ValueError, np.VisibleDeprecationWarning):
47 # argument was not a tuple representing a LinearConstraint
48 pass
50 # Address cases 4/5
51 for constraint in constraints:
52 # if it's not a LinearConstraint or something that represents a
53 # LinearConstraint at this point, it's invalid
54 if not isinstance(constraint, LinearConstraint):
55 try:
56 constraint = LinearConstraint(*constraint)
57 except TypeError as exc:
58 raise ValueError(message) from exc
59 As.append(csc_array(constraint.A))
60 b_ls.append(np.atleast_1d(constraint.lb).astype(np.double))
61 b_us.append(np.atleast_1d(constraint.ub).astype(np.double))
63 if len(As) > 1:
64 A = vstack(As)
65 b_l = np.concatenate(b_ls)
66 b_u = np.concatenate(b_us)
67 else: # avoid unnecessary copying
68 A = As[0]
69 b_l = b_ls[0]
70 b_u = b_us[0]
72 return A, b_l, b_u
75def _milp_iv(c, integrality, bounds, constraints, options):
76 # objective IV
77 c = np.atleast_1d(c).astype(np.double)
78 if c.ndim != 1 or c.size == 0 or not np.all(np.isfinite(c)):
79 message = ("`c` must be a one-dimensional array of finite numbers "
80 "with at least one element.")
81 raise ValueError(message)
83 # integrality IV
84 message = ("`integrality` must contain integers 0-3 and be broadcastable "
85 "to `c.shape`.")
86 if integrality is None:
87 integrality = 0
88 try:
89 integrality = np.broadcast_to(integrality, c.shape).astype(np.uint8)
90 except ValueError:
91 raise ValueError(message)
92 if integrality.min() < 0 or integrality.max() > 3:
93 raise ValueError(message)
95 # bounds IV
96 if bounds is None:
97 bounds = Bounds(0, np.inf)
98 elif not isinstance(bounds, Bounds):
99 message = ("`bounds` must be convertible into an instance of "
100 "`scipy.optimize.Bounds`.")
101 try:
102 bounds = Bounds(*bounds)
103 except TypeError as exc:
104 raise ValueError(message) from exc
106 try:
107 lb = np.broadcast_to(bounds.lb, c.shape).astype(np.double)
108 ub = np.broadcast_to(bounds.ub, c.shape).astype(np.double)
109 except (ValueError, TypeError) as exc:
110 message = ("`bounds.lb` and `bounds.ub` must contain reals and "
111 "be broadcastable to `c.shape`.")
112 raise ValueError(message) from exc
114 # constraints IV
115 if not constraints:
116 constraints = [LinearConstraint(np.empty((0, c.size)),
117 np.empty((0,)), np.empty((0,)))]
118 try:
119 A, b_l, b_u = _constraints_to_components(constraints)
120 except ValueError as exc:
121 message = ("`constraints` (or each element within `constraints`) must "
122 "be convertible into an instance of "
123 "`scipy.optimize.LinearConstraint`.")
124 raise ValueError(message) from exc
126 if A.shape != (b_l.size, c.size):
127 message = "The shape of `A` must be (len(b_l), len(c))."
128 raise ValueError(message)
129 indptr, indices, data = A.indptr, A.indices, A.data.astype(np.double)
131 # options IV
132 options = options or {}
133 supported_options = {'disp', 'presolve', 'time_limit', 'node_limit',
134 'mip_rel_gap'}
135 unsupported_options = set(options).difference(supported_options)
136 if unsupported_options:
137 message = (f"Unrecognized options detected: {unsupported_options}. "
138 "These will be passed to HiGHS verbatim.")
139 warnings.warn(message, RuntimeWarning, stacklevel=3)
140 options_iv = {'log_to_console': options.pop("disp", False),
141 'mip_max_nodes': options.pop("node_limit", None)}
142 options_iv.update(options)
144 return c, integrality, lb, ub, indptr, indices, data, b_l, b_u, options_iv
147def milp(c, *, integrality=None, bounds=None, constraints=None, options=None):
148 r"""
149 Mixed-integer linear programming
151 Solves problems of the following form:
153 .. math::
155 \min_x \ & c^T x \\
156 \mbox{such that} \ & b_l \leq A x \leq b_u,\\
157 & l \leq x \leq u, \\
158 & x_i \in \mathbb{Z}, i \in X_i
160 where :math:`x` is a vector of decision variables;
161 :math:`c`, :math:`b_l`, :math:`b_u`, :math:`l`, and :math:`u` are vectors;
162 :math:`A` is a matrix, and :math:`X_i` is the set of indices of
163 decision variables that must be integral. (In this context, a
164 variable that can assume only integer values is said to be "integral";
165 it has an "integrality" constraint.)
167 Alternatively, that's:
169 minimize::
171 c @ x
173 such that::
175 b_l <= A @ x <= b_u
176 l <= x <= u
177 Specified elements of x must be integers
179 By default, ``l = 0`` and ``u = np.inf`` unless specified with
180 ``bounds``.
182 Parameters
183 ----------
184 c : 1D array_like
185 The coefficients of the linear objective function to be minimized.
186 `c` is converted to a double precision array before the problem is
187 solved.
188 integrality : 1D array_like, optional
189 Indicates the type of integrality constraint on each decision variable.
191 ``0`` : Continuous variable; no integrality constraint.
193 ``1`` : Integer variable; decision variable must be an integer
194 within `bounds`.
196 ``2`` : Semi-continuous variable; decision variable must be within
197 `bounds` or take value ``0``.
199 ``3`` : Semi-integer variable; decision variable must be an integer
200 within `bounds` or take value ``0``.
202 By default, all variables are continuous. `integrality` is converted
203 to an array of integers before the problem is solved.
205 bounds : scipy.optimize.Bounds, optional
206 Bounds on the decision variables. Lower and upper bounds are converted
207 to double precision arrays before the problem is solved. The
208 ``keep_feasible`` parameter of the `Bounds` object is ignored. If
209 not specified, all decision variables are constrained to be
210 non-negative.
211 constraints : sequence of scipy.optimize.LinearConstraint, optional
212 Linear constraints of the optimization problem. Arguments may be
213 one of the following:
215 1. A single `LinearConstraint` object
216 2. A single tuple that can be converted to a `LinearConstraint` object
217 as ``LinearConstraint(*constraints)``
218 3. A sequence composed entirely of objects of type 1. and 2.
220 Before the problem is solved, all values are converted to double
221 precision, and the matrices of constraint coefficients are converted to
222 instances of `scipy.sparse.csc_array`. The ``keep_feasible`` parameter
223 of `LinearConstraint` objects is ignored.
224 options : dict, optional
225 A dictionary of solver options. The following keys are recognized.
227 disp : bool (default: ``False``)
228 Set to ``True`` if indicators of optimization status are to be
229 printed to the console during optimization.
230 node_limit : int, optional
231 The maximum number of nodes (linear program relaxations) to solve
232 before stopping. Default is no maximum number of nodes.
233 presolve : bool (default: ``True``)
234 Presolve attempts to identify trivial infeasibilities,
235 identify trivial unboundedness, and simplify the problem before
236 sending it to the main solver.
237 time_limit : float, optional
238 The maximum number of seconds allotted to solve the problem.
239 Default is no time limit.
240 mip_rel_gap : float, optional
241 Termination criterion for MIP solver: solver will terminate when
242 the gap between the primal objective value and the dual objective
243 bound, scaled by the primal objective value, is <= mip_rel_gap.
245 Returns
246 -------
247 res : OptimizeResult
248 An instance of :class:`scipy.optimize.OptimizeResult`. The object
249 is guaranteed to have the following attributes.
251 status : int
252 An integer representing the exit status of the algorithm.
254 ``0`` : Optimal solution found.
256 ``1`` : Iteration or time limit reached.
258 ``2`` : Problem is infeasible.
260 ``3`` : Problem is unbounded.
262 ``4`` : Other; see message for details.
264 success : bool
265 ``True`` when an optimal solution is found and ``False`` otherwise.
267 message : str
268 A string descriptor of the exit status of the algorithm.
270 The following attributes will also be present, but the values may be
271 ``None``, depending on the solution status.
273 x : ndarray
274 The values of the decision variables that minimize the
275 objective function while satisfying the constraints.
276 fun : float
277 The optimal value of the objective function ``c @ x``.
278 mip_node_count : int
279 The number of subproblems or "nodes" solved by the MILP solver.
280 mip_dual_bound : float
281 The MILP solver's final estimate of the lower bound on the optimal
282 solution.
283 mip_gap : float
284 The difference between the primal objective value and the dual
285 objective bound, scaled by the primal objective value.
287 Notes
288 -----
289 `milp` is a wrapper of the HiGHS linear optimization software [1]_. The
290 algorithm is deterministic, and it typically finds the global optimum of
291 moderately challenging mixed-integer linear programs (when it exists).
293 References
294 ----------
295 .. [1] Huangfu, Q., Galabova, I., Feldmeier, M., and Hall, J. A. J.
296 "HiGHS - high performance software for linear optimization."
297 https://highs.dev/
298 .. [2] Huangfu, Q. and Hall, J. A. J. "Parallelizing the dual revised
299 simplex method." Mathematical Programming Computation, 10 (1),
300 119-142, 2018. DOI: 10.1007/s12532-017-0130-5
302 Examples
303 --------
304 Consider the problem at
305 https://en.wikipedia.org/wiki/Integer_programming#Example, which is
306 expressed as a maximization problem of two variables. Since `milp` requires
307 that the problem be expressed as a minimization problem, the objective
308 function coefficients on the decision variables are:
310 >>> import numpy as np
311 >>> c = -np.array([0, 1])
313 Note the negative sign: we maximize the original objective function
314 by minimizing the negative of the objective function.
316 We collect the coefficients of the constraints into arrays like:
318 >>> A = np.array([[-1, 1], [3, 2], [2, 3]])
319 >>> b_u = np.array([1, 12, 12])
320 >>> b_l = np.full_like(b_u, -np.inf)
322 Because there is no lower limit on these constraints, we have defined a
323 variable ``b_l`` full of values representing negative infinity. This may
324 be unfamiliar to users of `scipy.optimize.linprog`, which only accepts
325 "less than" (or "upper bound") inequality constraints of the form
326 ``A_ub @ x <= b_u``. By accepting both ``b_l`` and ``b_u`` of constraints
327 ``b_l <= A_ub @ x <= b_u``, `milp` makes it easy to specify "greater than"
328 inequality constraints, "less than" inequality constraints, and equality
329 constraints concisely.
331 These arrays are collected into a single `LinearConstraint` object like:
333 >>> from scipy.optimize import LinearConstraint
334 >>> constraints = LinearConstraint(A, b_l, b_u)
336 The non-negativity bounds on the decision variables are enforced by
337 default, so we do not need to provide an argument for `bounds`.
339 Finally, the problem states that both decision variables must be integers:
341 >>> integrality = np.ones_like(c)
343 We solve the problem like:
345 >>> from scipy.optimize import milp
346 >>> res = milp(c=c, constraints=constraints, integrality=integrality)
347 >>> res.x
348 [1.0, 2.0]
350 Note that had we solved the relaxed problem (without integrality
351 constraints):
353 >>> res = milp(c=c, constraints=constraints) # OR:
354 >>> # from scipy.optimize import linprog; res = linprog(c, A, b_u)
355 >>> res.x
356 [1.8, 2.8]
358 we would not have obtained the correct solution by rounding to the nearest
359 integers.
361 Other examples are given :ref:`in the tutorial <tutorial-optimize_milp>`.
363 """
364 args_iv = _milp_iv(c, integrality, bounds, constraints, options)
365 c, integrality, lb, ub, indptr, indices, data, b_l, b_u, options = args_iv
367 highs_res = _highs_wrapper(c, indptr, indices, data, b_l, b_u,
368 lb, ub, integrality, options)
370 res = {}
372 # Convert to scipy-style status and message
373 highs_status = highs_res.get('status', None)
374 highs_message = highs_res.get('message', None)
375 status, message = _highs_to_scipy_status_message(highs_status,
376 highs_message)
377 res['status'] = status
378 res['message'] = message
379 res['success'] = (status == 0)
380 x = highs_res.get('x', None)
381 res['x'] = np.array(x) if x is not None else None
382 res['fun'] = highs_res.get('fun', None)
383 res['mip_node_count'] = highs_res.get('mip_node_count', None)
384 res['mip_dual_bound'] = highs_res.get('mip_dual_bound', None)
385 res['mip_gap'] = highs_res.get('mip_gap', None)
387 return OptimizeResult(res)