Coverage for /usr/lib/python3/dist-packages/scipy/optimize/_milp.py: 9%
110 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1import warnings
2import numpy as np
3from scipy.sparse import csc_array, vstack, issparse
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, format="csc")
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 if issparse(c):
78 raise ValueError("`c` must be a dense array.")
79 c = np.atleast_1d(c).astype(np.double)
80 if c.ndim != 1 or c.size == 0 or not np.all(np.isfinite(c)):
81 message = ("`c` must be a one-dimensional array of finite numbers "
82 "with at least one element.")
83 raise ValueError(message)
85 # integrality IV
86 if issparse(integrality):
87 raise ValueError("`integrality` must be a dense array.")
88 message = ("`integrality` must contain integers 0-3 and be broadcastable "
89 "to `c.shape`.")
90 if integrality is None:
91 integrality = 0
92 try:
93 integrality = np.broadcast_to(integrality, c.shape).astype(np.uint8)
94 except ValueError:
95 raise ValueError(message)
96 if integrality.min() < 0 or integrality.max() > 3:
97 raise ValueError(message)
99 # bounds IV
100 if bounds is None:
101 bounds = Bounds(0, np.inf)
102 elif not isinstance(bounds, Bounds):
103 message = ("`bounds` must be convertible into an instance of "
104 "`scipy.optimize.Bounds`.")
105 try:
106 bounds = Bounds(*bounds)
107 except TypeError as exc:
108 raise ValueError(message) from exc
110 try:
111 lb = np.broadcast_to(bounds.lb, c.shape).astype(np.double)
112 ub = np.broadcast_to(bounds.ub, c.shape).astype(np.double)
113 except (ValueError, TypeError) as exc:
114 message = ("`bounds.lb` and `bounds.ub` must contain reals and "
115 "be broadcastable to `c.shape`.")
116 raise ValueError(message) from exc
118 # constraints IV
119 if not constraints:
120 constraints = [LinearConstraint(np.empty((0, c.size)),
121 np.empty((0,)), np.empty((0,)))]
122 try:
123 A, b_l, b_u = _constraints_to_components(constraints)
124 except ValueError as exc:
125 message = ("`constraints` (or each element within `constraints`) must "
126 "be convertible into an instance of "
127 "`scipy.optimize.LinearConstraint`.")
128 raise ValueError(message) from exc
130 if A.shape != (b_l.size, c.size):
131 message = "The shape of `A` must be (len(b_l), len(c))."
132 raise ValueError(message)
133 indptr, indices, data = A.indptr, A.indices, A.data.astype(np.double)
135 # options IV
136 options = options or {}
137 supported_options = {'disp', 'presolve', 'time_limit', 'node_limit',
138 'mip_rel_gap'}
139 unsupported_options = set(options).difference(supported_options)
140 if unsupported_options:
141 message = (f"Unrecognized options detected: {unsupported_options}. "
142 "These will be passed to HiGHS verbatim.")
143 warnings.warn(message, RuntimeWarning, stacklevel=3)
144 options_iv = {'log_to_console': options.pop("disp", False),
145 'mip_max_nodes': options.pop("node_limit", None)}
146 options_iv.update(options)
148 return c, integrality, lb, ub, indptr, indices, data, b_l, b_u, options_iv
151def milp(c, *, integrality=None, bounds=None, constraints=None, options=None):
152 r"""
153 Mixed-integer linear programming
155 Solves problems of the following form:
157 .. math::
159 \min_x \ & c^T x \\
160 \mbox{such that} \ & b_l \leq A x \leq b_u,\\
161 & l \leq x \leq u, \\
162 & x_i \in \mathbb{Z}, i \in X_i
164 where :math:`x` is a vector of decision variables;
165 :math:`c`, :math:`b_l`, :math:`b_u`, :math:`l`, and :math:`u` are vectors;
166 :math:`A` is a matrix, and :math:`X_i` is the set of indices of
167 decision variables that must be integral. (In this context, a
168 variable that can assume only integer values is said to be "integral";
169 it has an "integrality" constraint.)
171 Alternatively, that's:
173 minimize::
175 c @ x
177 such that::
179 b_l <= A @ x <= b_u
180 l <= x <= u
181 Specified elements of x must be integers
183 By default, ``l = 0`` and ``u = np.inf`` unless specified with
184 ``bounds``.
186 Parameters
187 ----------
188 c : 1D dense array_like
189 The coefficients of the linear objective function to be minimized.
190 `c` is converted to a double precision array before the problem is
191 solved.
192 integrality : 1D dense array_like, optional
193 Indicates the type of integrality constraint on each decision variable.
195 ``0`` : Continuous variable; no integrality constraint.
197 ``1`` : Integer variable; decision variable must be an integer
198 within `bounds`.
200 ``2`` : Semi-continuous variable; decision variable must be within
201 `bounds` or take value ``0``.
203 ``3`` : Semi-integer variable; decision variable must be an integer
204 within `bounds` or take value ``0``.
206 By default, all variables are continuous. `integrality` is converted
207 to an array of integers before the problem is solved.
209 bounds : scipy.optimize.Bounds, optional
210 Bounds on the decision variables. Lower and upper bounds are converted
211 to double precision arrays before the problem is solved. The
212 ``keep_feasible`` parameter of the `Bounds` object is ignored. If
213 not specified, all decision variables are constrained to be
214 non-negative.
215 constraints : sequence of scipy.optimize.LinearConstraint, optional
216 Linear constraints of the optimization problem. Arguments may be
217 one of the following:
219 1. A single `LinearConstraint` object
220 2. A single tuple that can be converted to a `LinearConstraint` object
221 as ``LinearConstraint(*constraints)``
222 3. A sequence composed entirely of objects of type 1. and 2.
224 Before the problem is solved, all values are converted to double
225 precision, and the matrices of constraint coefficients are converted to
226 instances of `scipy.sparse.csc_array`. The ``keep_feasible`` parameter
227 of `LinearConstraint` objects is ignored.
228 options : dict, optional
229 A dictionary of solver options. The following keys are recognized.
231 disp : bool (default: ``False``)
232 Set to ``True`` if indicators of optimization status are to be
233 printed to the console during optimization.
234 node_limit : int, optional
235 The maximum number of nodes (linear program relaxations) to solve
236 before stopping. Default is no maximum number of nodes.
237 presolve : bool (default: ``True``)
238 Presolve attempts to identify trivial infeasibilities,
239 identify trivial unboundedness, and simplify the problem before
240 sending it to the main solver.
241 time_limit : float, optional
242 The maximum number of seconds allotted to solve the problem.
243 Default is no time limit.
244 mip_rel_gap : float, optional
245 Termination criterion for MIP solver: solver will terminate when
246 the gap between the primal objective value and the dual objective
247 bound, scaled by the primal objective value, is <= mip_rel_gap.
249 Returns
250 -------
251 res : OptimizeResult
252 An instance of :class:`scipy.optimize.OptimizeResult`. The object
253 is guaranteed to have the following attributes.
255 status : int
256 An integer representing the exit status of the algorithm.
258 ``0`` : Optimal solution found.
260 ``1`` : Iteration or time limit reached.
262 ``2`` : Problem is infeasible.
264 ``3`` : Problem is unbounded.
266 ``4`` : Other; see message for details.
268 success : bool
269 ``True`` when an optimal solution is found and ``False`` otherwise.
271 message : str
272 A string descriptor of the exit status of the algorithm.
274 The following attributes will also be present, but the values may be
275 ``None``, depending on the solution status.
277 x : ndarray
278 The values of the decision variables that minimize the
279 objective function while satisfying the constraints.
280 fun : float
281 The optimal value of the objective function ``c @ x``.
282 mip_node_count : int
283 The number of subproblems or "nodes" solved by the MILP solver.
284 mip_dual_bound : float
285 The MILP solver's final estimate of the lower bound on the optimal
286 solution.
287 mip_gap : float
288 The difference between the primal objective value and the dual
289 objective bound, scaled by the primal objective value.
291 Notes
292 -----
293 `milp` is a wrapper of the HiGHS linear optimization software [1]_. The
294 algorithm is deterministic, and it typically finds the global optimum of
295 moderately challenging mixed-integer linear programs (when it exists).
297 References
298 ----------
299 .. [1] Huangfu, Q., Galabova, I., Feldmeier, M., and Hall, J. A. J.
300 "HiGHS - high performance software for linear optimization."
301 https://highs.dev/
302 .. [2] Huangfu, Q. and Hall, J. A. J. "Parallelizing the dual revised
303 simplex method." Mathematical Programming Computation, 10 (1),
304 119-142, 2018. DOI: 10.1007/s12532-017-0130-5
306 Examples
307 --------
308 Consider the problem at
309 https://en.wikipedia.org/wiki/Integer_programming#Example, which is
310 expressed as a maximization problem of two variables. Since `milp` requires
311 that the problem be expressed as a minimization problem, the objective
312 function coefficients on the decision variables are:
314 >>> import numpy as np
315 >>> c = -np.array([0, 1])
317 Note the negative sign: we maximize the original objective function
318 by minimizing the negative of the objective function.
320 We collect the coefficients of the constraints into arrays like:
322 >>> A = np.array([[-1, 1], [3, 2], [2, 3]])
323 >>> b_u = np.array([1, 12, 12])
324 >>> b_l = np.full_like(b_u, -np.inf)
326 Because there is no lower limit on these constraints, we have defined a
327 variable ``b_l`` full of values representing negative infinity. This may
328 be unfamiliar to users of `scipy.optimize.linprog`, which only accepts
329 "less than" (or "upper bound") inequality constraints of the form
330 ``A_ub @ x <= b_u``. By accepting both ``b_l`` and ``b_u`` of constraints
331 ``b_l <= A_ub @ x <= b_u``, `milp` makes it easy to specify "greater than"
332 inequality constraints, "less than" inequality constraints, and equality
333 constraints concisely.
335 These arrays are collected into a single `LinearConstraint` object like:
337 >>> from scipy.optimize import LinearConstraint
338 >>> constraints = LinearConstraint(A, b_l, b_u)
340 The non-negativity bounds on the decision variables are enforced by
341 default, so we do not need to provide an argument for `bounds`.
343 Finally, the problem states that both decision variables must be integers:
345 >>> integrality = np.ones_like(c)
347 We solve the problem like:
349 >>> from scipy.optimize import milp
350 >>> res = milp(c=c, constraints=constraints, integrality=integrality)
351 >>> res.x
352 [1.0, 2.0]
354 Note that had we solved the relaxed problem (without integrality
355 constraints):
357 >>> res = milp(c=c, constraints=constraints) # OR:
358 >>> # from scipy.optimize import linprog; res = linprog(c, A, b_u)
359 >>> res.x
360 [1.8, 2.8]
362 we would not have obtained the correct solution by rounding to the nearest
363 integers.
365 Other examples are given :ref:`in the tutorial <tutorial-optimize_milp>`.
367 """
368 args_iv = _milp_iv(c, integrality, bounds, constraints, options)
369 c, integrality, lb, ub, indptr, indices, data, b_l, b_u, options = args_iv
371 highs_res = _highs_wrapper(c, indptr, indices, data, b_l, b_u,
372 lb, ub, integrality, options)
374 res = {}
376 # Convert to scipy-style status and message
377 highs_status = highs_res.get('status', None)
378 highs_message = highs_res.get('message', None)
379 status, message = _highs_to_scipy_status_message(highs_status,
380 highs_message)
381 res['status'] = status
382 res['message'] = message
383 res['success'] = (status == 0)
384 x = highs_res.get('x', None)
385 res['x'] = np.array(x) if x is not None else None
386 res['fun'] = highs_res.get('fun', None)
387 res['mip_node_count'] = highs_res.get('mip_node_count', None)
388 res['mip_dual_bound'] = highs_res.get('mip_dual_bound', None)
389 res['mip_gap'] = highs_res.get('mip_gap', None)
391 return OptimizeResult(res)