Coverage for /usr/lib/python3/dist-packages/scipy/optimize/_linprog_highs.py: 14%
77 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1"""HiGHS Linear Optimization Methods
3Interface to HiGHS linear optimization software.
4https://highs.dev/
6.. versionadded:: 1.5.0
8References
9----------
10.. [1] Q. Huangfu and J.A.J. Hall. "Parallelizing the dual revised simplex
11 method." Mathematical Programming Computation, 10 (1), 119-142,
12 2018. DOI: 10.1007/s12532-017-0130-5
14"""
16import inspect
17import numpy as np
18from ._optimize import OptimizeWarning, OptimizeResult
19from warnings import warn
20from ._highs._highs_wrapper import _highs_wrapper
21from ._highs._highs_constants import (
22 CONST_INF,
23 MESSAGE_LEVEL_NONE,
24 HIGHS_OBJECTIVE_SENSE_MINIMIZE,
26 MODEL_STATUS_NOTSET,
27 MODEL_STATUS_LOAD_ERROR,
28 MODEL_STATUS_MODEL_ERROR,
29 MODEL_STATUS_PRESOLVE_ERROR,
30 MODEL_STATUS_SOLVE_ERROR,
31 MODEL_STATUS_POSTSOLVE_ERROR,
32 MODEL_STATUS_MODEL_EMPTY,
33 MODEL_STATUS_OPTIMAL,
34 MODEL_STATUS_INFEASIBLE,
35 MODEL_STATUS_UNBOUNDED_OR_INFEASIBLE,
36 MODEL_STATUS_UNBOUNDED,
37 MODEL_STATUS_REACHED_DUAL_OBJECTIVE_VALUE_UPPER_BOUND
38 as MODEL_STATUS_RDOVUB,
39 MODEL_STATUS_REACHED_OBJECTIVE_TARGET,
40 MODEL_STATUS_REACHED_TIME_LIMIT,
41 MODEL_STATUS_REACHED_ITERATION_LIMIT,
43 HIGHS_SIMPLEX_STRATEGY_DUAL,
45 HIGHS_SIMPLEX_CRASH_STRATEGY_OFF,
47 HIGHS_SIMPLEX_EDGE_WEIGHT_STRATEGY_CHOOSE,
48 HIGHS_SIMPLEX_EDGE_WEIGHT_STRATEGY_DANTZIG,
49 HIGHS_SIMPLEX_EDGE_WEIGHT_STRATEGY_DEVEX,
50 HIGHS_SIMPLEX_EDGE_WEIGHT_STRATEGY_STEEPEST_EDGE,
51)
52from scipy.sparse import csc_matrix, vstack, issparse
55def _highs_to_scipy_status_message(highs_status, highs_message):
56 """Converts HiGHS status number/message to SciPy status number/message"""
58 scipy_statuses_messages = {
59 None: (4, "HiGHS did not provide a status code. "),
60 MODEL_STATUS_NOTSET: (4, ""),
61 MODEL_STATUS_LOAD_ERROR: (4, ""),
62 MODEL_STATUS_MODEL_ERROR: (2, ""),
63 MODEL_STATUS_PRESOLVE_ERROR: (4, ""),
64 MODEL_STATUS_SOLVE_ERROR: (4, ""),
65 MODEL_STATUS_POSTSOLVE_ERROR: (4, ""),
66 MODEL_STATUS_MODEL_EMPTY: (4, ""),
67 MODEL_STATUS_RDOVUB: (4, ""),
68 MODEL_STATUS_REACHED_OBJECTIVE_TARGET: (4, ""),
69 MODEL_STATUS_OPTIMAL: (0, "Optimization terminated successfully. "),
70 MODEL_STATUS_REACHED_TIME_LIMIT: (1, "Time limit reached. "),
71 MODEL_STATUS_REACHED_ITERATION_LIMIT: (1, "Iteration limit reached. "),
72 MODEL_STATUS_INFEASIBLE: (2, "The problem is infeasible. "),
73 MODEL_STATUS_UNBOUNDED: (3, "The problem is unbounded. "),
74 MODEL_STATUS_UNBOUNDED_OR_INFEASIBLE: (4, "The problem is unbounded "
75 "or infeasible. ")}
76 unrecognized = (4, "The HiGHS status code was not recognized. ")
77 scipy_status, scipy_message = (
78 scipy_statuses_messages.get(highs_status, unrecognized))
79 scipy_message = (f"{scipy_message}"
80 f"(HiGHS Status {highs_status}: {highs_message})")
81 return scipy_status, scipy_message
84def _replace_inf(x):
85 # Replace `np.inf` with CONST_INF
86 infs = np.isinf(x)
87 with np.errstate(invalid="ignore"):
88 x[infs] = np.sign(x[infs])*CONST_INF
89 return x
92def _convert_to_highs_enum(option, option_str, choices):
93 # If option is in the choices we can look it up, if not use
94 # the default value taken from function signature and warn:
95 try:
96 return choices[option.lower()]
97 except AttributeError:
98 return choices[option]
99 except KeyError:
100 sig = inspect.signature(_linprog_highs)
101 default_str = sig.parameters[option_str].default
102 warn(f"Option {option_str} is {option}, but only values in "
103 f"{set(choices.keys())} are allowed. Using default: "
104 f"{default_str}.",
105 OptimizeWarning, stacklevel=3)
106 return choices[default_str]
109def _linprog_highs(lp, solver, time_limit=None, presolve=True,
110 disp=False, maxiter=None,
111 dual_feasibility_tolerance=None,
112 primal_feasibility_tolerance=None,
113 ipm_optimality_tolerance=None,
114 simplex_dual_edge_weight_strategy=None,
115 mip_rel_gap=None,
116 mip_max_nodes=None,
117 **unknown_options):
118 r"""
119 Solve the following linear programming problem using one of the HiGHS
120 solvers:
122 User-facing documentation is in _linprog_doc.py.
124 Parameters
125 ----------
126 lp : _LPProblem
127 A ``scipy.optimize._linprog_util._LPProblem`` ``namedtuple``.
128 solver : "ipm" or "simplex" or None
129 Which HiGHS solver to use. If ``None``, "simplex" will be used.
131 Options
132 -------
133 maxiter : int
134 The maximum number of iterations to perform in either phase. For
135 ``solver='ipm'``, this does not include the number of crossover
136 iterations. Default is the largest possible value for an ``int``
137 on the platform.
138 disp : bool
139 Set to ``True`` if indicators of optimization status are to be printed
140 to the console each iteration; default ``False``.
141 time_limit : float
142 The maximum time in seconds allotted to solve the problem; default is
143 the largest possible value for a ``double`` on the platform.
144 presolve : bool
145 Presolve attempts to identify trivial infeasibilities,
146 identify trivial unboundedness, and simplify the problem before
147 sending it to the main solver. It is generally recommended
148 to keep the default setting ``True``; set to ``False`` if presolve is
149 to be disabled.
150 dual_feasibility_tolerance : double
151 Dual feasibility tolerance. Default is 1e-07.
152 The minimum of this and ``primal_feasibility_tolerance``
153 is used for the feasibility tolerance when ``solver='ipm'``.
154 primal_feasibility_tolerance : double
155 Primal feasibility tolerance. Default is 1e-07.
156 The minimum of this and ``dual_feasibility_tolerance``
157 is used for the feasibility tolerance when ``solver='ipm'``.
158 ipm_optimality_tolerance : double
159 Optimality tolerance for ``solver='ipm'``. Default is 1e-08.
160 Minimum possible value is 1e-12 and must be smaller than the largest
161 possible value for a ``double`` on the platform.
162 simplex_dual_edge_weight_strategy : str (default: None)
163 Strategy for simplex dual edge weights. The default, ``None``,
164 automatically selects one of the following.
166 ``'dantzig'`` uses Dantzig's original strategy of choosing the most
167 negative reduced cost.
169 ``'devex'`` uses the strategy described in [15]_.
171 ``steepest`` uses the exact steepest edge strategy as described in
172 [16]_.
174 ``'steepest-devex'`` begins with the exact steepest edge strategy
175 until the computation is too costly or inexact and then switches to
176 the devex method.
178 Curently, using ``None`` always selects ``'steepest-devex'``, but this
179 may change as new options become available.
181 mip_max_nodes : int
182 The maximum number of nodes allotted to solve the problem; default is
183 the largest possible value for a ``HighsInt`` on the platform.
184 Ignored if not using the MIP solver.
185 unknown_options : dict
186 Optional arguments not used by this particular solver. If
187 ``unknown_options`` is non-empty, a warning is issued listing all
188 unused options.
190 Returns
191 -------
192 sol : dict
193 A dictionary consisting of the fields:
195 x : 1D array
196 The values of the decision variables that minimizes the
197 objective function while satisfying the constraints.
198 fun : float
199 The optimal value of the objective function ``c @ x``.
200 slack : 1D array
201 The (nominally positive) values of the slack,
202 ``b_ub - A_ub @ x``.
203 con : 1D array
204 The (nominally zero) residuals of the equality constraints,
205 ``b_eq - A_eq @ x``.
206 success : bool
207 ``True`` when the algorithm succeeds in finding an optimal
208 solution.
209 status : int
210 An integer representing the exit status of the algorithm.
212 ``0`` : Optimization terminated successfully.
214 ``1`` : Iteration or time limit reached.
216 ``2`` : Problem appears to be infeasible.
218 ``3`` : Problem appears to be unbounded.
220 ``4`` : The HiGHS solver ran into a problem.
222 message : str
223 A string descriptor of the exit status of the algorithm.
224 nit : int
225 The total number of iterations performed.
226 For ``solver='simplex'``, this includes iterations in all
227 phases. For ``solver='ipm'``, this does not include
228 crossover iterations.
229 crossover_nit : int
230 The number of primal/dual pushes performed during the
231 crossover routine for ``solver='ipm'``. This is ``0``
232 for ``solver='simplex'``.
233 ineqlin : OptimizeResult
234 Solution and sensitivity information corresponding to the
235 inequality constraints, `b_ub`. A dictionary consisting of the
236 fields:
238 residual : np.ndnarray
239 The (nominally positive) values of the slack variables,
240 ``b_ub - A_ub @ x``. This quantity is also commonly
241 referred to as "slack".
243 marginals : np.ndarray
244 The sensitivity (partial derivative) of the objective
245 function with respect to the right-hand side of the
246 inequality constraints, `b_ub`.
248 eqlin : OptimizeResult
249 Solution and sensitivity information corresponding to the
250 equality constraints, `b_eq`. A dictionary consisting of the
251 fields:
253 residual : np.ndarray
254 The (nominally zero) residuals of the equality constraints,
255 ``b_eq - A_eq @ x``.
257 marginals : np.ndarray
258 The sensitivity (partial derivative) of the objective
259 function with respect to the right-hand side of the
260 equality constraints, `b_eq`.
262 lower, upper : OptimizeResult
263 Solution and sensitivity information corresponding to the
264 lower and upper bounds on decision variables, `bounds`.
266 residual : np.ndarray
267 The (nominally positive) values of the quantity
268 ``x - lb`` (lower) or ``ub - x`` (upper).
270 marginals : np.ndarray
271 The sensitivity (partial derivative) of the objective
272 function with respect to the lower and upper
273 `bounds`.
275 mip_node_count : int
276 The number of subproblems or "nodes" solved by the MILP
277 solver. Only present when `integrality` is not `None`.
279 mip_dual_bound : float
280 The MILP solver's final estimate of the lower bound on the
281 optimal solution. Only present when `integrality` is not
282 `None`.
284 mip_gap : float
285 The difference between the final objective function value
286 and the final dual bound, scaled by the final objective
287 function value. Only present when `integrality` is not
288 `None`.
290 Notes
291 -----
292 The result fields `ineqlin`, `eqlin`, `lower`, and `upper` all contain
293 `marginals`, or partial derivatives of the objective function with respect
294 to the right-hand side of each constraint. These partial derivatives are
295 also referred to as "Lagrange multipliers", "dual values", and
296 "shadow prices". The sign convention of `marginals` is opposite that
297 of Lagrange multipliers produced by many nonlinear solvers.
299 References
300 ----------
301 .. [15] Harris, Paula MJ. "Pivot selection methods of the Devex LP code."
302 Mathematical programming 5.1 (1973): 1-28.
303 .. [16] Goldfarb, Donald, and John Ker Reid. "A practicable steepest-edge
304 simplex algorithm." Mathematical Programming 12.1 (1977): 361-371.
305 """
306 if unknown_options:
307 message = (f"Unrecognized options detected: {unknown_options}. "
308 "These will be passed to HiGHS verbatim.")
309 warn(message, OptimizeWarning, stacklevel=3)
311 # Map options to HiGHS enum values
312 simplex_dual_edge_weight_strategy_enum = _convert_to_highs_enum(
313 simplex_dual_edge_weight_strategy,
314 'simplex_dual_edge_weight_strategy',
315 choices={'dantzig': HIGHS_SIMPLEX_EDGE_WEIGHT_STRATEGY_DANTZIG,
316 'devex': HIGHS_SIMPLEX_EDGE_WEIGHT_STRATEGY_DEVEX,
317 'steepest-devex': HIGHS_SIMPLEX_EDGE_WEIGHT_STRATEGY_CHOOSE,
318 'steepest':
319 HIGHS_SIMPLEX_EDGE_WEIGHT_STRATEGY_STEEPEST_EDGE,
320 None: None})
322 c, A_ub, b_ub, A_eq, b_eq, bounds, x0, integrality = lp
324 lb, ub = bounds.T.copy() # separate bounds, copy->C-cntgs
325 # highs_wrapper solves LHS <= A*x <= RHS, not equality constraints
326 with np.errstate(invalid="ignore"):
327 lhs_ub = -np.ones_like(b_ub)*np.inf # LHS of UB constraints is -inf
328 rhs_ub = b_ub # RHS of UB constraints is b_ub
329 lhs_eq = b_eq # Equality constaint is inequality
330 rhs_eq = b_eq # constraint with LHS=RHS
331 lhs = np.concatenate((lhs_ub, lhs_eq))
332 rhs = np.concatenate((rhs_ub, rhs_eq))
334 if issparse(A_ub) or issparse(A_eq):
335 A = vstack((A_ub, A_eq))
336 else:
337 A = np.vstack((A_ub, A_eq))
338 A = csc_matrix(A)
340 options = {
341 'presolve': presolve,
342 'sense': HIGHS_OBJECTIVE_SENSE_MINIMIZE,
343 'solver': solver,
344 'time_limit': time_limit,
345 'highs_debug_level': MESSAGE_LEVEL_NONE,
346 'dual_feasibility_tolerance': dual_feasibility_tolerance,
347 'ipm_optimality_tolerance': ipm_optimality_tolerance,
348 'log_to_console': disp,
349 'mip_max_nodes': mip_max_nodes,
350 'output_flag': disp,
351 'primal_feasibility_tolerance': primal_feasibility_tolerance,
352 'simplex_dual_edge_weight_strategy':
353 simplex_dual_edge_weight_strategy_enum,
354 'simplex_strategy': HIGHS_SIMPLEX_STRATEGY_DUAL,
355 'simplex_crash_strategy': HIGHS_SIMPLEX_CRASH_STRATEGY_OFF,
356 'ipm_iteration_limit': maxiter,
357 'simplex_iteration_limit': maxiter,
358 'mip_rel_gap': mip_rel_gap,
359 }
360 options.update(unknown_options)
362 # np.inf doesn't work; use very large constant
363 rhs = _replace_inf(rhs)
364 lhs = _replace_inf(lhs)
365 lb = _replace_inf(lb)
366 ub = _replace_inf(ub)
368 if integrality is None or np.sum(integrality) == 0:
369 integrality = np.empty(0)
370 else:
371 integrality = np.array(integrality)
373 res = _highs_wrapper(c, A.indptr, A.indices, A.data, lhs, rhs,
374 lb, ub, integrality.astype(np.uint8), options)
376 # HiGHS represents constraints as lhs/rhs, so
377 # Ax + s = b => Ax = b - s
378 # and we need to split up s by A_ub and A_eq
379 if 'slack' in res:
380 slack = res['slack']
381 con = np.array(slack[len(b_ub):])
382 slack = np.array(slack[:len(b_ub)])
383 else:
384 slack, con = None, None
386 # lagrange multipliers for equalities/inequalities and upper/lower bounds
387 if 'lambda' in res:
388 lamda = res['lambda']
389 marg_ineqlin = np.array(lamda[:len(b_ub)])
390 marg_eqlin = np.array(lamda[len(b_ub):])
391 marg_upper = np.array(res['marg_bnds'][1, :])
392 marg_lower = np.array(res['marg_bnds'][0, :])
393 else:
394 marg_ineqlin, marg_eqlin = None, None
395 marg_upper, marg_lower = None, None
397 # this needs to be updated if we start choosing the solver intelligently
399 # Convert to scipy-style status and message
400 highs_status = res.get('status', None)
401 highs_message = res.get('message', None)
402 status, message = _highs_to_scipy_status_message(highs_status,
403 highs_message)
405 x = np.array(res['x']) if 'x' in res else None
406 sol = {'x': x,
407 'slack': slack,
408 'con': con,
409 'ineqlin': OptimizeResult({
410 'residual': slack,
411 'marginals': marg_ineqlin,
412 }),
413 'eqlin': OptimizeResult({
414 'residual': con,
415 'marginals': marg_eqlin,
416 }),
417 'lower': OptimizeResult({
418 'residual': None if x is None else x - lb,
419 'marginals': marg_lower,
420 }),
421 'upper': OptimizeResult({
422 'residual': None if x is None else ub - x,
423 'marginals': marg_upper
424 }),
425 'fun': res.get('fun'),
426 'status': status,
427 'success': res['status'] == MODEL_STATUS_OPTIMAL,
428 'message': message,
429 'nit': res.get('simplex_nit', 0) or res.get('ipm_nit', 0),
430 'crossover_nit': res.get('crossover_nit'),
431 }
433 if np.any(x) and integrality is not None:
434 sol.update({
435 'mip_node_count': res.get('mip_node_count', 0),
436 'mip_dual_bound': res.get('mip_dual_bound', 0.0),
437 'mip_gap': res.get('mip_gap', 0.0),
438 })
440 return sol