1""" 2Method agnostic utility functions for linear progamming 3""" 4 5import numpy as np 6import scipy.sparse as sps 7from warnings import warn 8from .optimize import OptimizeWarning 9from scipy.optimize._remove_redundancy import ( 10 _remove_redundancy_svd, _remove_redundancy_pivot_sparse, 11 _remove_redundancy_pivot_dense, _remove_redundancy_id 12 ) 13from collections import namedtuple 14 15_LPProblem = namedtuple('_LPProblem', 'c A_ub b_ub A_eq b_eq bounds x0') 16_LPProblem.__new__.__defaults__ = (None,) * 6 # make c the only required arg 17_LPProblem.__doc__ = \ 18 """ Represents a linear-programming problem. 19 20 Attributes 21 ---------- 22 c : 1D array 23 The coefficients of the linear objective function to be minimized. 24 A_ub : 2D array, optional 25 The inequality constraint matrix. Each row of ``A_ub`` specifies the 26 coefficients of a linear inequality constraint on ``x``. 27 b_ub : 1D array, optional 28 The inequality constraint vector. Each element represents an 29 upper bound on the corresponding value of ``A_ub @ x``. 30 A_eq : 2D array, optional 31 The equality constraint matrix. Each row of ``A_eq`` specifies the 32 coefficients of a linear equality constraint on ``x``. 33 b_eq : 1D array, optional 34 The equality constraint vector. Each element of ``A_eq @ x`` must equal 35 the corresponding element of ``b_eq``. 36 bounds : various valid formats, optional 37 The bounds of ``x``, as ``min`` and ``max`` pairs. 38 If bounds are specified for all N variables separately, valid formats 39 are: 40 * a 2D array (N x 2); 41 * a sequence of N sequences, each with 2 values. 42 If all variables have the same bounds, the bounds can be specified as 43 a 1-D or 2-D array or sequence with 2 scalar values. 44 If all variables have a lower bound of 0 and no upper bound, the bounds 45 parameter can be omitted (or given as None). 46 Absent lower and/or upper bounds can be specified as -numpy.inf (no 47 lower bound), numpy.inf (no upper bound) or None (both). 48 x0 : 1D array, optional 49 Guess values of the decision variables, which will be refined by 50 the optimization algorithm. This argument is currently used only by the 51 'revised simplex' method, and can only be used if `x0` represents a 52 basic feasible solution. 53 54 Notes 55 ----- 56 This namedtuple supports 2 ways of initialization: 57 >>> lp1 = _LPProblem(c=[-1, 4], A_ub=[[-3, 1], [1, 2]], b_ub=[6, 4]) 58 >>> lp2 = _LPProblem([-1, 4], [[-3, 1], [1, 2]], [6, 4]) 59 60 Note that only ``c`` is a required argument here, whereas all other arguments 61 ``A_ub``, ``b_ub``, ``A_eq``, ``b_eq``, ``bounds``, ``x0`` are optional with 62 default values of None. 63 For example, ``A_eq`` and ``b_eq`` can be set without ``A_ub`` or ``b_ub``: 64 >>> lp3 = _LPProblem(c=[-1, 4], A_eq=[[2, 1]], b_eq=[10]) 65 """ 66 67 68def _check_sparse_inputs(options, meth, A_ub, A_eq): 69 """ 70 Check the provided ``A_ub`` and ``A_eq`` matrices conform to the specified 71 optional sparsity variables. 72 73 Parameters 74 ---------- 75 A_ub : 2-D array, optional 76 2-D array such that ``A_ub @ x`` gives the values of the upper-bound 77 inequality constraints at ``x``. 78 A_eq : 2-D array, optional 79 2-D array such that ``A_eq @ x`` gives the values of the equality 80 constraints at ``x``. 81 options : dict 82 A dictionary of solver options. All methods accept the following 83 generic options: 84 85 maxiter : int 86 Maximum number of iterations to perform. 87 disp : bool 88 Set to True to print convergence messages. 89 90 For method-specific options, see :func:`show_options('linprog')`. 91 method : str, optional 92 The algorithm used to solve the standard form problem. 93 94 Returns 95 ------- 96 A_ub : 2-D array, optional 97 2-D array such that ``A_ub @ x`` gives the values of the upper-bound 98 inequality constraints at ``x``. 99 A_eq : 2-D array, optional 100 2-D array such that ``A_eq @ x`` gives the values of the equality 101 constraints at ``x``. 102 options : dict 103 A dictionary of solver options. All methods accept the following 104 generic options: 105 106 maxiter : int 107 Maximum number of iterations to perform. 108 disp : bool 109 Set to True to print convergence messages. 110 111 For method-specific options, see :func:`show_options('linprog')`. 112 """ 113 # This is an undocumented option for unit testing sparse presolve 114 _sparse_presolve = options.pop('_sparse_presolve', False) 115 if _sparse_presolve and A_eq is not None: 116 A_eq = sps.coo_matrix(A_eq) 117 if _sparse_presolve and A_ub is not None: 118 A_ub = sps.coo_matrix(A_ub) 119 120 sparse_constraint = sps.issparse(A_eq) or sps.issparse(A_ub) 121 122 preferred_methods = {"highs", "highs-ds", "highs-ipm"} 123 dense_methods = {"simplex", "revised simplex"} 124 if meth in dense_methods and sparse_constraint: 125 raise ValueError(f"Method '{meth}' does not support sparse " 126 "constraint matrices. Please consider using one of " 127 f"{preferred_methods}.") 128 129 sparse = options.get('sparse', False) 130 if not sparse and sparse_constraint and meth == 'interior-point': 131 options['sparse'] = True 132 warn("Sparse constraint matrix detected; setting 'sparse':True.", 133 OptimizeWarning, stacklevel=4) 134 return options, A_ub, A_eq 135 136 137def _format_A_constraints(A, n_x, sparse_lhs=False): 138 """Format the left hand side of the constraints to a 2-D array 139 140 Parameters 141 ---------- 142 A : 2-D array 143 2-D array such that ``A @ x`` gives the values of the upper-bound 144 (in)equality constraints at ``x``. 145 n_x : int 146 The number of variables in the linear programming problem. 147 sparse_lhs : bool 148 Whether either of `A_ub` or `A_eq` are sparse. If true return a 149 coo_matrix instead of a numpy array. 150 151 Returns 152 ------- 153 np.ndarray or sparse.coo_matrix 154 2-D array such that ``A @ x`` gives the values of the upper-bound 155 (in)equality constraints at ``x``. 156 157 """ 158 if sparse_lhs: 159 return sps.coo_matrix( 160 (0, n_x) if A is None else A, dtype=float, copy=True 161 ) 162 elif A is None: 163 return np.zeros((0, n_x), dtype=float) 164 else: 165 return np.array(A, dtype=float, copy=True) 166 167 168def _format_b_constraints(b): 169 """Format the upper bounds of the constraints to a 1-D array 170 171 Parameters 172 ---------- 173 b : 1-D array 174 1-D array of values representing the upper-bound of each (in)equality 175 constraint (row) in ``A``. 176 177 Returns 178 ------- 179 1-D np.array 180 1-D array of values representing the upper-bound of each (in)equality 181 constraint (row) in ``A``. 182 183 """ 184 if b is None: 185 return np.array([], dtype=float) 186 b = np.array(b, dtype=float, copy=True).squeeze() 187 return b if b.size != 1 else b.reshape((-1)) 188 189 190def _clean_inputs(lp): 191 """ 192 Given user inputs for a linear programming problem, return the 193 objective vector, upper bound constraints, equality constraints, 194 and simple bounds in a preferred format. 195 196 Parameters 197 ---------- 198 lp : A `scipy.optimize._linprog_util._LPProblem` consisting of the following fields: 199 200 c : 1D array 201 The coefficients of the linear objective function to be minimized. 202 A_ub : 2D array, optional 203 The inequality constraint matrix. Each row of ``A_ub`` specifies the 204 coefficients of a linear inequality constraint on ``x``. 205 b_ub : 1D array, optional 206 The inequality constraint vector. Each element represents an 207 upper bound on the corresponding value of ``A_ub @ x``. 208 A_eq : 2D array, optional 209 The equality constraint matrix. Each row of ``A_eq`` specifies the 210 coefficients of a linear equality constraint on ``x``. 211 b_eq : 1D array, optional 212 The equality constraint vector. Each element of ``A_eq @ x`` must equal 213 the corresponding element of ``b_eq``. 214 bounds : various valid formats, optional 215 The bounds of ``x``, as ``min`` and ``max`` pairs. 216 If bounds are specified for all N variables separately, valid formats are: 217 * a 2D array (2 x N or N x 2); 218 * a sequence of N sequences, each with 2 values. 219 If all variables have the same bounds, a single pair of values can 220 be specified. Valid formats are: 221 * a sequence with 2 scalar values; 222 * a sequence with a single element containing 2 scalar values. 223 If all variables have a lower bound of 0 and no upper bound, the bounds 224 parameter can be omitted (or given as None). 225 x0 : 1D array, optional 226 Guess values of the decision variables, which will be refined by 227 the optimization algorithm. This argument is currently used only by the 228 'revised simplex' method, and can only be used if `x0` represents a 229 basic feasible solution. 230 231 Returns 232 ------- 233 lp : A `scipy.optimize._linprog_util._LPProblem` consisting of the following fields: 234 235 c : 1D array 236 The coefficients of the linear objective function to be minimized. 237 A_ub : 2D array, optional 238 The inequality constraint matrix. Each row of ``A_ub`` specifies the 239 coefficients of a linear inequality constraint on ``x``. 240 b_ub : 1D array, optional 241 The inequality constraint vector. Each element represents an 242 upper bound on the corresponding value of ``A_ub @ x``. 243 A_eq : 2D array, optional 244 The equality constraint matrix. Each row of ``A_eq`` specifies the 245 coefficients of a linear equality constraint on ``x``. 246 b_eq : 1D array, optional 247 The equality constraint vector. Each element of ``A_eq @ x`` must equal 248 the corresponding element of ``b_eq``. 249 bounds : 2D array 250 The bounds of ``x``, as ``min`` and ``max`` pairs, one for each of the N 251 elements of ``x``. The N x 2 array contains lower bounds in the first 252 column and upper bounds in the 2nd. Unbounded variables have lower 253 bound -np.inf and/or upper bound np.inf. 254 x0 : 1D array, optional 255 Guess values of the decision variables, which will be refined by 256 the optimization algorithm. This argument is currently used only by the 257 'revised simplex' method, and can only be used if `x0` represents a 258 basic feasible solution. 259 260 """ 261 c, A_ub, b_ub, A_eq, b_eq, bounds, x0 = lp 262 263 if c is None: 264 raise TypeError 265 266 try: 267 c = np.array(c, dtype=np.float64, copy=True).squeeze() 268 except ValueError as e: 269 raise TypeError( 270 "Invalid input for linprog: c must be a 1-D array of numerical " 271 "coefficients") from e 272 else: 273 # If c is a single value, convert it to a 1-D array. 274 if c.size == 1: 275 c = c.reshape((-1)) 276 277 n_x = len(c) 278 if n_x == 0 or len(c.shape) != 1: 279 raise ValueError( 280 "Invalid input for linprog: c must be a 1-D array and must " 281 "not have more than one non-singleton dimension") 282 if not(np.isfinite(c).all()): 283 raise ValueError( 284 "Invalid input for linprog: c must not contain values " 285 "inf, nan, or None") 286 287 sparse_lhs = sps.issparse(A_eq) or sps.issparse(A_ub) 288 try: 289 A_ub = _format_A_constraints(A_ub, n_x, sparse_lhs=sparse_lhs) 290 except ValueError as e: 291 raise TypeError( 292 "Invalid input for linprog: A_ub must be a 2-D array " 293 "of numerical values") from e 294 else: 295 n_ub = A_ub.shape[0] 296 if len(A_ub.shape) != 2 or A_ub.shape[1] != n_x: 297 raise ValueError( 298 "Invalid input for linprog: A_ub must have exactly two " 299 "dimensions, and the number of columns in A_ub must be " 300 "equal to the size of c") 301 if (sps.issparse(A_ub) and not np.isfinite(A_ub.data).all() 302 or not sps.issparse(A_ub) and not np.isfinite(A_ub).all()): 303 raise ValueError( 304 "Invalid input for linprog: A_ub must not contain values " 305 "inf, nan, or None") 306 307 try: 308 b_ub = _format_b_constraints(b_ub) 309 except ValueError as e: 310 raise TypeError( 311 "Invalid input for linprog: b_ub must be a 1-D array of " 312 "numerical values, each representing the upper bound of an " 313 "inequality constraint (row) in A_ub") from e 314 else: 315 if b_ub.shape != (n_ub,): 316 raise ValueError( 317 "Invalid input for linprog: b_ub must be a 1-D array; b_ub " 318 "must not have more than one non-singleton dimension and " 319 "the number of rows in A_ub must equal the number of values " 320 "in b_ub") 321 if not(np.isfinite(b_ub).all()): 322 raise ValueError( 323 "Invalid input for linprog: b_ub must not contain values " 324 "inf, nan, or None") 325 326 try: 327 A_eq = _format_A_constraints(A_eq, n_x, sparse_lhs=sparse_lhs) 328 except ValueError as e: 329 raise TypeError( 330 "Invalid input for linprog: A_eq must be a 2-D array " 331 "of numerical values") from e 332 else: 333 n_eq = A_eq.shape[0] 334 if len(A_eq.shape) != 2 or A_eq.shape[1] != n_x: 335 raise ValueError( 336 "Invalid input for linprog: A_eq must have exactly two " 337 "dimensions, and the number of columns in A_eq must be " 338 "equal to the size of c") 339 340 if (sps.issparse(A_eq) and not np.isfinite(A_eq.data).all() 341 or not sps.issparse(A_eq) and not np.isfinite(A_eq).all()): 342 raise ValueError( 343 "Invalid input for linprog: A_eq must not contain values " 344 "inf, nan, or None") 345 346 try: 347 b_eq = _format_b_constraints(b_eq) 348 except ValueError as e: 349 raise TypeError( 350 "Invalid input for linprog: b_eq must be a dense, 1-D array of " 351 "numerical values, each representing the right hand side of an " 352 "equality constraint (row) in A_eq") from e 353 else: 354 if b_eq.shape != (n_eq,): 355 raise ValueError( 356 "Invalid input for linprog: b_eq must be a 1-D array; b_eq " 357 "must not have more than one non-singleton dimension and " 358 "the number of rows in A_eq must equal the number of values " 359 "in b_eq") 360 if not(np.isfinite(b_eq).all()): 361 raise ValueError( 362 "Invalid input for linprog: b_eq must not contain values " 363 "inf, nan, or None") 364 365 # x0 gives a (optional) starting solution to the solver. If x0 is None, 366 # skip the checks. Initial solution will be generated automatically. 367 if x0 is not None: 368 try: 369 x0 = np.array(x0, dtype=float, copy=True).squeeze() 370 except ValueError as e: 371 raise TypeError( 372 "Invalid input for linprog: x0 must be a 1-D array of " 373 "numerical coefficients") from e 374 if x0.ndim == 0: 375 x0 = x0.reshape((-1)) 376 if len(x0) == 0 or x0.ndim != 1: 377 raise ValueError( 378 "Invalid input for linprog: x0 should be a 1-D array; it " 379 "must not have more than one non-singleton dimension") 380 if not x0.size == c.size: 381 raise ValueError( 382 "Invalid input for linprog: x0 and c should contain the " 383 "same number of elements") 384 if not np.isfinite(x0).all(): 385 raise ValueError( 386 "Invalid input for linprog: x0 must not contain values " 387 "inf, nan, or None") 388 389 # Bounds can be one of these formats: 390 # (1) a 2-D array or sequence, with shape N x 2 391 # (2) a 1-D or 2-D sequence or array with 2 scalars 392 # (3) None (or an empty sequence or array) 393 # Unspecified bounds can be represented by None or (-)np.inf. 394 # All formats are converted into a N x 2 np.array with (-)np.inf where 395 # bounds are unspecified. 396 397 # Prepare clean bounds array 398 bounds_clean = np.zeros((n_x, 2), dtype=float) 399 400 # Convert to a numpy array. 401 # np.array(..,dtype=float) raises an error if dimensions are inconsistent 402 # or if there are invalid data types in bounds. Just add a linprog prefix 403 # to the error and re-raise. 404 # Creating at least a 2-D array simplifies the cases to distinguish below. 405 if bounds is None or np.array_equal(bounds, []) or np.array_equal(bounds, [[]]): 406 bounds = (0, np.inf) 407 try: 408 bounds_conv = np.atleast_2d(np.array(bounds, dtype=float)) 409 except ValueError as e: 410 raise ValueError( 411 "Invalid input for linprog: unable to interpret bounds, " 412 "check values and dimensions: " + e.args[0]) from e 413 except TypeError as e: 414 raise TypeError( 415 "Invalid input for linprog: unable to interpret bounds, " 416 "check values and dimensions: " + e.args[0]) from e 417 418 # Check bounds options 419 bsh = bounds_conv.shape 420 if len(bsh) > 2: 421 # Do not try to handle multidimensional bounds input 422 raise ValueError( 423 "Invalid input for linprog: provide a 2-D array for bounds, " 424 "not a {:d}-D array.".format(len(bsh))) 425 elif np.all(bsh == (n_x, 2)): 426 # Regular N x 2 array 427 bounds_clean = bounds_conv 428 elif (np.all(bsh == (2, 1)) or np.all(bsh == (1, 2))): 429 # 2 values: interpret as overall lower and upper bound 430 bounds_flat = bounds_conv.flatten() 431 bounds_clean[:, 0] = bounds_flat[0] 432 bounds_clean[:, 1] = bounds_flat[1] 433 elif np.all(bsh == (2, n_x)): 434 # Reject a 2 x N array 435 raise ValueError( 436 "Invalid input for linprog: provide a {:d} x 2 array for bounds, " 437 "not a 2 x {:d} array.".format(n_x, n_x)) 438 else: 439 raise ValueError( 440 "Invalid input for linprog: unable to interpret bounds with this " 441 "dimension tuple: {0}.".format(bsh)) 442 443 # The process above creates nan-s where the input specified None 444 # Convert the nan-s in the 1st column to -np.inf and in the 2nd column 445 # to np.inf 446 i_none = np.isnan(bounds_clean[:, 0]) 447 bounds_clean[i_none, 0] = -np.inf 448 i_none = np.isnan(bounds_clean[:, 1]) 449 bounds_clean[i_none, 1] = np.inf 450 451 return _LPProblem(c, A_ub, b_ub, A_eq, b_eq, bounds_clean, x0) 452 453 454def _presolve(lp, rr, rr_method, tol=1e-9): 455 """ 456 Given inputs for a linear programming problem in preferred format, 457 presolve the problem: identify trivial infeasibilities, redundancies, 458 and unboundedness, tighten bounds where possible, and eliminate fixed 459 variables. 460 461 Parameters 462 ---------- 463 lp : A `scipy.optimize._linprog_util._LPProblem` consisting of the following fields: 464 465 c : 1D array 466 The coefficients of the linear objective function to be minimized. 467 A_ub : 2D array, optional 468 The inequality constraint matrix. Each row of ``A_ub`` specifies the 469 coefficients of a linear inequality constraint on ``x``. 470 b_ub : 1D array, optional 471 The inequality constraint vector. Each element represents an 472 upper bound on the corresponding value of ``A_ub @ x``. 473 A_eq : 2D array, optional 474 The equality constraint matrix. Each row of ``A_eq`` specifies the 475 coefficients of a linear equality constraint on ``x``. 476 b_eq : 1D array, optional 477 The equality constraint vector. Each element of ``A_eq @ x`` must equal 478 the corresponding element of ``b_eq``. 479 bounds : 2D array 480 The bounds of ``x``, as ``min`` and ``max`` pairs, one for each of the N 481 elements of ``x``. The N x 2 array contains lower bounds in the first 482 column and upper bounds in the 2nd. Unbounded variables have lower 483 bound -np.inf and/or upper bound np.inf. 484 x0 : 1D array, optional 485 Guess values of the decision variables, which will be refined by 486 the optimization algorithm. This argument is currently used only by the 487 'revised simplex' method, and can only be used if `x0` represents a 488 basic feasible solution. 489 490 rr : bool 491 If ``True`` attempts to eliminate any redundant rows in ``A_eq``. 492 Set False if ``A_eq`` is known to be of full row rank, or if you are 493 looking for a potential speedup (at the expense of reliability). 494 rr_method : string 495 Method used to identify and remove redundant rows from the 496 equality constraint matrix after presolve. 497 tol : float 498 The tolerance which determines when a solution is "close enough" to 499 zero in Phase 1 to be considered a basic feasible solution or close 500 enough to positive to serve as an optimal solution. 501 502 Returns 503 ------- 504 lp : A `scipy.optimize._linprog_util._LPProblem` consisting of the following fields: 505 506 c : 1D array 507 The coefficients of the linear objective function to be minimized. 508 A_ub : 2D array, optional 509 The inequality constraint matrix. Each row of ``A_ub`` specifies the 510 coefficients of a linear inequality constraint on ``x``. 511 b_ub : 1D array, optional 512 The inequality constraint vector. Each element represents an 513 upper bound on the corresponding value of ``A_ub @ x``. 514 A_eq : 2D array, optional 515 The equality constraint matrix. Each row of ``A_eq`` specifies the 516 coefficients of a linear equality constraint on ``x``. 517 b_eq : 1D array, optional 518 The equality constraint vector. Each element of ``A_eq @ x`` must equal 519 the corresponding element of ``b_eq``. 520 bounds : 2D array 521 The bounds of ``x``, as ``min`` and ``max`` pairs, possibly tightened. 522 x0 : 1D array, optional 523 Guess values of the decision variables, which will be refined by 524 the optimization algorithm. This argument is currently used only by the 525 'revised simplex' method, and can only be used if `x0` represents a 526 basic feasible solution. 527 528 c0 : 1D array 529 Constant term in objective function due to fixed (and eliminated) 530 variables. 531 x : 1D array 532 Solution vector (when the solution is trivial and can be determined 533 in presolve) 534 revstack: list of functions 535 the functions in the list reverse the operations of _presolve() 536 the function signature is x_org = f(x_mod), where x_mod is the result 537 of a presolve step and x_org the value at the start of the step 538 (currently, the revstack contains only one function) 539 complete: bool 540 Whether the solution is complete (solved or determined to be infeasible 541 or unbounded in presolve) 542 status : int 543 An integer representing the exit status of the optimization:: 544 545 0 : Optimization terminated successfully 546 1 : Iteration limit reached 547 2 : Problem appears to be infeasible 548 3 : Problem appears to be unbounded 549 4 : Serious numerical difficulties encountered 550 551 message : str 552 A string descriptor of the exit status of the optimization. 553 554 References 555 ---------- 556 .. [5] Andersen, Erling D. "Finding all linearly dependent rows in 557 large-scale linear programming." Optimization Methods and Software 558 6.3 (1995): 219-227. 559 .. [8] Andersen, Erling D., and Knud D. Andersen. "Presolving in linear 560 programming." Mathematical Programming 71.2 (1995): 221-245. 561 562 """ 563 # ideas from Reference [5] by Andersen and Andersen 564 # however, unlike the reference, this is performed before converting 565 # problem to standard form 566 # There are a few advantages: 567 # * artificial variables have not been added, so matrices are smaller 568 # * bounds have not been converted to constraints yet. (It is better to 569 # do that after presolve because presolve may adjust the simple bounds.) 570 # There are many improvements that can be made, namely: 571 # * implement remaining checks from [5] 572 # * loop presolve until no additional changes are made 573 # * implement additional efficiency improvements in redundancy removal [2] 574 575 c, A_ub, b_ub, A_eq, b_eq, bounds, x0 = lp 576 577 revstack = [] # record of variables eliminated from problem 578 # constant term in cost function may be added if variables are eliminated 579 c0 = 0 580 complete = False # complete is True if detected infeasible/unbounded 581 x = np.zeros(c.shape) # this is solution vector if completed in presolve 582 583 status = 0 # all OK unless determined otherwise 584 message = "" 585 586 # Lower and upper bounds. Copy to prevent feedback. 587 lb = bounds[:, 0].copy() 588 ub = bounds[:, 1].copy() 589 590 m_eq, n = A_eq.shape 591 m_ub, n = A_ub.shape 592 593 if (rr_method is not None 594 and rr_method.lower() not in {"svd", "pivot", "id"}): 595 message = ("'" + str(rr_method) + "' is not a valid option " 596 "for redundancy removal. Valid options are 'SVD', " 597 "'pivot', and 'ID'.") 598 raise ValueError(message) 599 600 if sps.issparse(A_eq): 601 A_eq = A_eq.tocsr() 602 A_ub = A_ub.tocsr() 603 604 def where(A): 605 return A.nonzero() 606 607 vstack = sps.vstack 608 else: 609 where = np.where 610 vstack = np.vstack 611 612 # upper bounds > lower bounds 613 if np.any(ub < lb) or np.any(lb == np.inf) or np.any(ub == -np.inf): 614 status = 2 615 message = ("The problem is (trivially) infeasible since one " 616 "or more upper bounds are smaller than the corresponding " 617 "lower bounds, a lower bound is np.inf or an upper bound " 618 "is -np.inf.") 619 complete = True 620 return (_LPProblem(c, A_ub, b_ub, A_eq, b_eq, bounds, x0), 621 c0, x, revstack, complete, status, message) 622 623 # zero row in equality constraints 624 zero_row = np.array(np.sum(A_eq != 0, axis=1) == 0).flatten() 625 if np.any(zero_row): 626 if np.any( 627 np.logical_and( 628 zero_row, 629 np.abs(b_eq) > tol)): # test_zero_row_1 630 # infeasible if RHS is not zero 631 status = 2 632 message = ("The problem is (trivially) infeasible due to a row " 633 "of zeros in the equality constraint matrix with a " 634 "nonzero corresponding constraint value.") 635 complete = True 636 return (_LPProblem(c, A_ub, b_ub, A_eq, b_eq, bounds, x0), 637 c0, x, revstack, complete, status, message) 638 else: # test_zero_row_2 639 # if RHS is zero, we can eliminate this equation entirely 640 A_eq = A_eq[np.logical_not(zero_row), :] 641 b_eq = b_eq[np.logical_not(zero_row)] 642 643 # zero row in inequality constraints 644 zero_row = np.array(np.sum(A_ub != 0, axis=1) == 0).flatten() 645 if np.any(zero_row): 646 if np.any(np.logical_and(zero_row, b_ub < -tol)): # test_zero_row_1 647 # infeasible if RHS is less than zero (because LHS is zero) 648 status = 2 649 message = ("The problem is (trivially) infeasible due to a row " 650 "of zeros in the equality constraint matrix with a " 651 "nonzero corresponding constraint value.") 652 complete = True 653 return (_LPProblem(c, A_ub, b_ub, A_eq, b_eq, bounds, x0), 654 c0, x, revstack, complete, status, message) 655 else: # test_zero_row_2 656 # if LHS is >= 0, we can eliminate this constraint entirely 657 A_ub = A_ub[np.logical_not(zero_row), :] 658 b_ub = b_ub[np.logical_not(zero_row)] 659 660 # zero column in (both) constraints 661 # this indicates that a variable isn't constrained and can be removed 662 A = vstack((A_eq, A_ub)) 663 if A.shape[0] > 0: 664 zero_col = np.array(np.sum(A != 0, axis=0) == 0).flatten() 665 # variable will be at upper or lower bound, depending on objective 666 x[np.logical_and(zero_col, c < 0)] = ub[ 667 np.logical_and(zero_col, c < 0)] 668 x[np.logical_and(zero_col, c > 0)] = lb[ 669 np.logical_and(zero_col, c > 0)] 670 if np.any(np.isinf(x)): # if an unconstrained variable has no bound 671 status = 3 672 message = ("If feasible, the problem is (trivially) unbounded " 673 "due to a zero column in the constraint matrices. If " 674 "you wish to check whether the problem is infeasible, " 675 "turn presolve off.") 676 complete = True 677 return (_LPProblem(c, A_ub, b_ub, A_eq, b_eq, bounds, x0), 678 c0, x, revstack, complete, status, message) 679 # variables will equal upper/lower bounds will be removed later 680 lb[np.logical_and(zero_col, c < 0)] = ub[ 681 np.logical_and(zero_col, c < 0)] 682 ub[np.logical_and(zero_col, c > 0)] = lb[ 683 np.logical_and(zero_col, c > 0)] 684 685 # row singleton in equality constraints 686 # this fixes a variable and removes the constraint 687 singleton_row = np.array(np.sum(A_eq != 0, axis=1) == 1).flatten() 688 rows = where(singleton_row)[0] 689 cols = where(A_eq[rows, :])[1] 690 if len(rows) > 0: 691 for row, col in zip(rows, cols): 692 val = b_eq[row] / A_eq[row, col] 693 if not lb[col] - tol <= val <= ub[col] + tol: 694 # infeasible if fixed value is not within bounds 695 status = 2 696 message = ("The problem is (trivially) infeasible because a " 697 "singleton row in the equality constraints is " 698 "inconsistent with the bounds.") 699 complete = True 700 return (_LPProblem(c, A_ub, b_ub, A_eq, b_eq, bounds, x0), 701 c0, x, revstack, complete, status, message) 702 else: 703 # sets upper and lower bounds at that fixed value - variable 704 # will be removed later 705 lb[col] = val 706 ub[col] = val 707 A_eq = A_eq[np.logical_not(singleton_row), :] 708 b_eq = b_eq[np.logical_not(singleton_row)] 709 710 # row singleton in inequality constraints 711 # this indicates a simple bound and the constraint can be removed 712 # simple bounds may be adjusted here 713 # After all of the simple bound information is combined here, get_Abc will 714 # turn the simple bounds into constraints 715 singleton_row = np.array(np.sum(A_ub != 0, axis=1) == 1).flatten() 716 cols = where(A_ub[singleton_row, :])[1] 717 rows = where(singleton_row)[0] 718 if len(rows) > 0: 719 for row, col in zip(rows, cols): 720 val = b_ub[row] / A_ub[row, col] 721 if A_ub[row, col] > 0: # upper bound 722 if val < lb[col] - tol: # infeasible 723 complete = True 724 elif val < ub[col]: # new upper bound 725 ub[col] = val 726 else: # lower bound 727 if val > ub[col] + tol: # infeasible 728 complete = True 729 elif val > lb[col]: # new lower bound 730 lb[col] = val 731 if complete: 732 status = 2 733 message = ("The problem is (trivially) infeasible because a " 734 "singleton row in the upper bound constraints is " 735 "inconsistent with the bounds.") 736 return (_LPProblem(c, A_ub, b_ub, A_eq, b_eq, bounds, x0), 737 c0, x, revstack, complete, status, message) 738 A_ub = A_ub[np.logical_not(singleton_row), :] 739 b_ub = b_ub[np.logical_not(singleton_row)] 740 741 # identical bounds indicate that variable can be removed 742 i_f = np.abs(lb - ub) < tol # indices of "fixed" variables 743 i_nf = np.logical_not(i_f) # indices of "not fixed" variables 744 745 # test_bounds_equal_but_infeasible 746 if np.all(i_f): # if bounds define solution, check for consistency 747 residual = b_eq - A_eq.dot(lb) 748 slack = b_ub - A_ub.dot(lb) 749 if ((A_ub.size > 0 and np.any(slack < 0)) or 750 (A_eq.size > 0 and not np.allclose(residual, 0))): 751 status = 2 752 message = ("The problem is (trivially) infeasible because the " 753 "bounds fix all variables to values inconsistent with " 754 "the constraints") 755 complete = True 756 return (_LPProblem(c, A_ub, b_ub, A_eq, b_eq, bounds, x0), 757 c0, x, revstack, complete, status, message) 758 759 ub_mod = ub 760 lb_mod = lb 761 if np.any(i_f): 762 c0 += c[i_f].dot(lb[i_f]) 763 b_eq = b_eq - A_eq[:, i_f].dot(lb[i_f]) 764 b_ub = b_ub - A_ub[:, i_f].dot(lb[i_f]) 765 c = c[i_nf] 766 x_undo = lb[i_f] # not x[i_f], x is just zeroes 767 x = x[i_nf] 768 # user guess x0 stays separate from presolve solution x 769 if x0 is not None: 770 x0 = x0[i_nf] 771 A_eq = A_eq[:, i_nf] 772 A_ub = A_ub[:, i_nf] 773 # modify bounds 774 lb_mod = lb[i_nf] 775 ub_mod = ub[i_nf] 776 777 def rev(x_mod): 778 # Function to restore x: insert x_undo into x_mod. 779 # When elements have been removed at positions k1, k2, k3, ... 780 # then these must be replaced at (after) positions k1-1, k2-2, 781 # k3-3, ... in the modified array to recreate the original 782 i = np.flatnonzero(i_f) 783 # Number of variables to restore 784 N = len(i) 785 index_offset = np.arange(N) 786 # Create insert indices 787 insert_indices = i - index_offset 788 x_rev = np.insert(x_mod.astype(float), insert_indices, x_undo) 789 return x_rev 790 791 # Use revstack as a list of functions, currently just this one. 792 revstack.append(rev) 793 794 # no constraints indicates that problem is trivial 795 if A_eq.size == 0 and A_ub.size == 0: 796 b_eq = np.array([]) 797 b_ub = np.array([]) 798 # test_empty_constraint_1 799 if c.size == 0: 800 status = 0 801 message = ("The solution was determined in presolve as there are " 802 "no non-trivial constraints.") 803 elif (np.any(np.logical_and(c < 0, ub_mod == np.inf)) or 804 np.any(np.logical_and(c > 0, lb_mod == -np.inf))): 805 # test_no_constraints() 806 # test_unbounded_no_nontrivial_constraints_1 807 # test_unbounded_no_nontrivial_constraints_2 808 status = 3 809 message = ("The problem is (trivially) unbounded " 810 "because there are no non-trivial constraints and " 811 "a) at least one decision variable is unbounded " 812 "above and its corresponding cost is negative, or " 813 "b) at least one decision variable is unbounded below " 814 "and its corresponding cost is positive. ") 815 else: # test_empty_constraint_2 816 status = 0 817 message = ("The solution was determined in presolve as there are " 818 "no non-trivial constraints.") 819 complete = True 820 x[c < 0] = ub_mod[c < 0] 821 x[c > 0] = lb_mod[c > 0] 822 # where c is zero, set x to a finite bound or zero 823 x_zero_c = ub_mod[c == 0] 824 x_zero_c[np.isinf(x_zero_c)] = ub_mod[c == 0][np.isinf(x_zero_c)] 825 x_zero_c[np.isinf(x_zero_c)] = 0 826 x[c == 0] = x_zero_c 827 # if this is not the last step of presolve, should convert bounds back 828 # to array and return here 829 830 # Convert modified lb and ub back into N x 2 bounds 831 bounds = np.hstack((lb_mod[:, np.newaxis], ub_mod[:, np.newaxis])) 832 833 # remove redundant (linearly dependent) rows from equality constraints 834 n_rows_A = A_eq.shape[0] 835 redundancy_warning = ("A_eq does not appear to be of full row rank. To " 836 "improve performance, check the problem formulation " 837 "for redundant equality constraints.") 838 if (sps.issparse(A_eq)): 839 if rr and A_eq.size > 0: # TODO: Fast sparse rank check? 840 rr_res = _remove_redundancy_pivot_sparse(A_eq, b_eq) 841 A_eq, b_eq, status, message = rr_res 842 if A_eq.shape[0] < n_rows_A: 843 warn(redundancy_warning, OptimizeWarning, stacklevel=1) 844 if status != 0: 845 complete = True 846 return (_LPProblem(c, A_ub, b_ub, A_eq, b_eq, bounds, x0), 847 c0, x, revstack, complete, status, message) 848 849 # This is a wild guess for which redundancy removal algorithm will be 850 # faster. More testing would be good. 851 small_nullspace = 5 852 if rr and A_eq.size > 0: 853 try: # TODO: use results of first SVD in _remove_redundancy_svd 854 rank = np.linalg.matrix_rank(A_eq) 855 # oh well, we'll have to go with _remove_redundancy_pivot_dense 856 except Exception: 857 rank = 0 858 if rr and A_eq.size > 0 and rank < A_eq.shape[0]: 859 warn(redundancy_warning, OptimizeWarning, stacklevel=3) 860 dim_row_nullspace = A_eq.shape[0]-rank 861 if rr_method is None: 862 if dim_row_nullspace <= small_nullspace: 863 rr_res = _remove_redundancy_svd(A_eq, b_eq) 864 A_eq, b_eq, status, message = rr_res 865 if dim_row_nullspace > small_nullspace or status == 4: 866 rr_res = _remove_redundancy_pivot_dense(A_eq, b_eq) 867 A_eq, b_eq, status, message = rr_res 868 869 else: 870 rr_method = rr_method.lower() 871 if rr_method == "svd": 872 rr_res = _remove_redundancy_svd(A_eq, b_eq) 873 A_eq, b_eq, status, message = rr_res 874 elif rr_method == "pivot": 875 rr_res = _remove_redundancy_pivot_dense(A_eq, b_eq) 876 A_eq, b_eq, status, message = rr_res 877 elif rr_method == "id": 878 rr_res = _remove_redundancy_id(A_eq, b_eq, rank) 879 A_eq, b_eq, status, message = rr_res 880 else: # shouldn't get here; option validity checked above 881 pass 882 if A_eq.shape[0] < rank: 883 message = ("Due to numerical issues, redundant equality " 884 "constraints could not be removed automatically. " 885 "Try providing your constraint matrices as sparse " 886 "matrices to activate sparse presolve, try turning " 887 "off redundancy removal, or try turning off presolve " 888 "altogether.") 889 status = 4 890 if status != 0: 891 complete = True 892 return (_LPProblem(c, A_ub, b_ub, A_eq, b_eq, bounds, x0), 893 c0, x, revstack, complete, status, message) 894 895 896def _parse_linprog(lp, options, meth): 897 """ 898 Parse the provided linear programming problem 899 900 ``_parse_linprog`` employs two main steps ``_check_sparse_inputs`` and 901 ``_clean_inputs``. ``_check_sparse_inputs`` checks for sparsity in the 902 provided constraints (``A_ub`` and ``A_eq) and if these match the provided 903 sparsity optional values. 904 905 ``_clean inputs`` checks of the provided inputs. If no violations are 906 identified the objective vector, upper bound constraints, equality 907 constraints, and simple bounds are returned in the expected format. 908 909 Parameters 910 ---------- 911 lp : A `scipy.optimize._linprog_util._LPProblem` consisting of the following fields: 912 913 c : 1D array 914 The coefficients of the linear objective function to be minimized. 915 A_ub : 2D array, optional 916 The inequality constraint matrix. Each row of ``A_ub`` specifies the 917 coefficients of a linear inequality constraint on ``x``. 918 b_ub : 1D array, optional 919 The inequality constraint vector. Each element represents an 920 upper bound on the corresponding value of ``A_ub @ x``. 921 A_eq : 2D array, optional 922 The equality constraint matrix. Each row of ``A_eq`` specifies the 923 coefficients of a linear equality constraint on ``x``. 924 b_eq : 1D array, optional 925 The equality constraint vector. Each element of ``A_eq @ x`` must equal 926 the corresponding element of ``b_eq``. 927 bounds : various valid formats, optional 928 The bounds of ``x``, as ``min`` and ``max`` pairs. 929 If bounds are specified for all N variables separately, valid formats are: 930 * a 2D array (2 x N or N x 2); 931 * a sequence of N sequences, each with 2 values. 932 If all variables have the same bounds, a single pair of values can 933 be specified. Valid formats are: 934 * a sequence with 2 scalar values; 935 * a sequence with a single element containing 2 scalar values. 936 If all variables have a lower bound of 0 and no upper bound, the bounds 937 parameter can be omitted (or given as None). 938 x0 : 1D array, optional 939 Guess values of the decision variables, which will be refined by 940 the optimization algorithm. This argument is currently used only by the 941 'revised simplex' method, and can only be used if `x0` represents a 942 basic feasible solution. 943 944 options : dict 945 A dictionary of solver options. All methods accept the following 946 generic options: 947 948 maxiter : int 949 Maximum number of iterations to perform. 950 disp : bool 951 Set to True to print convergence messages. 952 953 For method-specific options, see :func:`show_options('linprog')`. 954 955 Returns 956 ------- 957 lp : A `scipy.optimize._linprog_util._LPProblem` consisting of the following fields: 958 959 c : 1D array 960 The coefficients of the linear objective function to be minimized. 961 A_ub : 2D array, optional 962 The inequality constraint matrix. Each row of ``A_ub`` specifies the 963 coefficients of a linear inequality constraint on ``x``. 964 b_ub : 1D array, optional 965 The inequality constraint vector. Each element represents an 966 upper bound on the corresponding value of ``A_ub @ x``. 967 A_eq : 2D array, optional 968 The equality constraint matrix. Each row of ``A_eq`` specifies the 969 coefficients of a linear equality constraint on ``x``. 970 b_eq : 1D array, optional 971 The equality constraint vector. Each element of ``A_eq @ x`` must equal 972 the corresponding element of ``b_eq``. 973 bounds : 2D array 974 The bounds of ``x``, as ``min`` and ``max`` pairs, one for each of the N 975 elements of ``x``. The N x 2 array contains lower bounds in the first 976 column and upper bounds in the 2nd. Unbounded variables have lower 977 bound -np.inf and/or upper bound np.inf. 978 x0 : 1D array, optional 979 Guess values of the decision variables, which will be refined by 980 the optimization algorithm. This argument is currently used only by the 981 'revised simplex' method, and can only be used if `x0` represents a 982 basic feasible solution. 983 984 options : dict, optional 985 A dictionary of solver options. All methods accept the following 986 generic options: 987 988 maxiter : int 989 Maximum number of iterations to perform. 990 disp : bool 991 Set to True to print convergence messages. 992 993 For method-specific options, see :func:`show_options('linprog')`. 994 995 """ 996 if options is None: 997 options = {} 998 999 solver_options = {k: v for k, v in options.items()} 1000 solver_options, A_ub, A_eq = _check_sparse_inputs(solver_options, meth, 1001 lp.A_ub, lp.A_eq) 1002 # Convert lists to numpy arrays, etc... 1003 lp = _clean_inputs(lp._replace(A_ub=A_ub, A_eq=A_eq)) 1004 return lp, solver_options 1005 1006 1007def _get_Abc(lp, c0): 1008 """ 1009 Given a linear programming problem of the form: 1010 1011 Minimize:: 1012 1013 c @ x 1014 1015 Subject to:: 1016 1017 A_ub @ x <= b_ub 1018 A_eq @ x == b_eq 1019 lb <= x <= ub 1020 1021 where ``lb = 0`` and ``ub = None`` unless set in ``bounds``. 1022 1023 Return the problem in standard form: 1024 1025 Minimize:: 1026 1027 c @ x 1028 1029 Subject to:: 1030 1031 A @ x == b 1032 x >= 0 1033 1034 by adding slack variables and making variable substitutions as necessary. 1035 1036 Parameters 1037 ---------- 1038 lp : A `scipy.optimize._linprog_util._LPProblem` consisting of the following fields: 1039 1040 c : 1D array 1041 The coefficients of the linear objective function to be minimized. 1042 A_ub : 2D array, optional 1043 The inequality constraint matrix. Each row of ``A_ub`` specifies the 1044 coefficients of a linear inequality constraint on ``x``. 1045 b_ub : 1D array, optional 1046 The inequality constraint vector. Each element represents an 1047 upper bound on the corresponding value of ``A_ub @ x``. 1048 A_eq : 2D array, optional 1049 The equality constraint matrix. Each row of ``A_eq`` specifies the 1050 coefficients of a linear equality constraint on ``x``. 1051 b_eq : 1D array, optional 1052 The equality constraint vector. Each element of ``A_eq @ x`` must equal 1053 the corresponding element of ``b_eq``. 1054 bounds : 2D array 1055 The bounds of ``x``, lower bounds in the 1st column, upper 1056 bounds in the 2nd column. The bounds are possibly tightened 1057 by the presolve procedure. 1058 x0 : 1D array, optional 1059 Guess values of the decision variables, which will be refined by 1060 the optimization algorithm. This argument is currently used only by the 1061 'revised simplex' method, and can only be used if `x0` represents a 1062 basic feasible solution. 1063 1064 c0 : float 1065 Constant term in objective function due to fixed (and eliminated) 1066 variables. 1067 1068 Returns 1069 ------- 1070 A : 2-D array 1071 2-D array such that ``A`` @ ``x``, gives the values of the equality 1072 constraints at ``x``. 1073 b : 1-D array 1074 1-D array of values representing the RHS of each equality constraint 1075 (row) in A (for standard form problem). 1076 c : 1-D array 1077 Coefficients of the linear objective function to be minimized (for 1078 standard form problem). 1079 c0 : float 1080 Constant term in objective function due to fixed (and eliminated) 1081 variables. 1082 x0 : 1-D array 1083 Starting values of the independent variables, which will be refined by 1084 the optimization algorithm 1085 1086 References 1087 ---------- 1088 .. [9] Bertsimas, Dimitris, and J. Tsitsiklis. "Introduction to linear 1089 programming." Athena Scientific 1 (1997): 997. 1090 1091 """ 1092 c, A_ub, b_ub, A_eq, b_eq, bounds, x0 = lp 1093 1094 if sps.issparse(A_eq): 1095 sparse = True 1096 A_eq = sps.csr_matrix(A_eq) 1097 A_ub = sps.csr_matrix(A_ub) 1098 1099 def hstack(blocks): 1100 return sps.hstack(blocks, format="csr") 1101 1102 def vstack(blocks): 1103 return sps.vstack(blocks, format="csr") 1104 1105 zeros = sps.csr_matrix 1106 eye = sps.eye 1107 else: 1108 sparse = False 1109 hstack = np.hstack 1110 vstack = np.vstack 1111 zeros = np.zeros 1112 eye = np.eye 1113 1114 # Variables lbs and ubs (see below) may be changed, which feeds back into 1115 # bounds, so copy. 1116 bounds = np.array(bounds, copy=True) 1117 1118 # modify problem such that all variables have only non-negativity bounds 1119 lbs = bounds[:, 0] 1120 ubs = bounds[:, 1] 1121 m_ub, n_ub = A_ub.shape 1122 1123 lb_none = np.equal(lbs, -np.inf) 1124 ub_none = np.equal(ubs, np.inf) 1125 lb_some = np.logical_not(lb_none) 1126 ub_some = np.logical_not(ub_none) 1127 1128 # unbounded below: substitute xi = -xi' (unbounded above) 1129 # if -inf <= xi <= ub, then -ub <= -xi <= inf, so swap and invert bounds 1130 l_nolb_someub = np.logical_and(lb_none, ub_some) 1131 i_nolb = np.nonzero(l_nolb_someub)[0] 1132 lbs[l_nolb_someub], ubs[l_nolb_someub] = ( 1133 -ubs[l_nolb_someub], -lbs[l_nolb_someub]) 1134 lb_none = np.equal(lbs, -np.inf) 1135 ub_none = np.equal(ubs, np.inf) 1136 lb_some = np.logical_not(lb_none) 1137 ub_some = np.logical_not(ub_none) 1138 c[i_nolb] *= -1 1139 if x0 is not None: 1140 x0[i_nolb] *= -1 1141 if len(i_nolb) > 0: 1142 if A_ub.shape[0] > 0: # sometimes needed for sparse arrays... weird 1143 A_ub[:, i_nolb] *= -1 1144 if A_eq.shape[0] > 0: 1145 A_eq[:, i_nolb] *= -1 1146 1147 # upper bound: add inequality constraint 1148 i_newub, = ub_some.nonzero() 1149 ub_newub = ubs[ub_some] 1150 n_bounds = len(i_newub) 1151 if n_bounds > 0: 1152 shape = (n_bounds, A_ub.shape[1]) 1153 if sparse: 1154 idxs = (np.arange(n_bounds), i_newub) 1155 A_ub = vstack((A_ub, sps.csr_matrix((np.ones(n_bounds), idxs), 1156 shape=shape))) 1157 else: 1158 A_ub = vstack((A_ub, np.zeros(shape))) 1159 A_ub[np.arange(m_ub, A_ub.shape[0]), i_newub] = 1 1160 b_ub = np.concatenate((b_ub, np.zeros(n_bounds))) 1161 b_ub[m_ub:] = ub_newub 1162 1163 A1 = vstack((A_ub, A_eq)) 1164 b = np.concatenate((b_ub, b_eq)) 1165 c = np.concatenate((c, np.zeros((A_ub.shape[0],)))) 1166 if x0 is not None: 1167 x0 = np.concatenate((x0, np.zeros((A_ub.shape[0],)))) 1168 # unbounded: substitute xi = xi+ + xi- 1169 l_free = np.logical_and(lb_none, ub_none) 1170 i_free = np.nonzero(l_free)[0] 1171 n_free = len(i_free) 1172 c = np.concatenate((c, np.zeros(n_free))) 1173 if x0 is not None: 1174 x0 = np.concatenate((x0, np.zeros(n_free))) 1175 A1 = hstack((A1[:, :n_ub], -A1[:, i_free])) 1176 c[n_ub:n_ub+n_free] = -c[i_free] 1177 if x0 is not None: 1178 i_free_neg = x0[i_free] < 0 1179 x0[np.arange(n_ub, A1.shape[1])[i_free_neg]] = -x0[i_free[i_free_neg]] 1180 x0[i_free[i_free_neg]] = 0 1181 1182 # add slack variables 1183 A2 = vstack([eye(A_ub.shape[0]), zeros((A_eq.shape[0], A_ub.shape[0]))]) 1184 1185 A = hstack([A1, A2]) 1186 1187 # lower bound: substitute xi = xi' + lb 1188 # now there is a constant term in objective 1189 i_shift = np.nonzero(lb_some)[0] 1190 lb_shift = lbs[lb_some].astype(float) 1191 c0 += np.sum(lb_shift * c[i_shift]) 1192 if sparse: 1193 b = b.reshape(-1, 1) 1194 A = A.tocsc() 1195 b -= (A[:, i_shift] * sps.diags(lb_shift)).sum(axis=1) 1196 b = b.ravel() 1197 else: 1198 b -= (A[:, i_shift] * lb_shift).sum(axis=1) 1199 if x0 is not None: 1200 x0[i_shift] -= lb_shift 1201 1202 return A, b, c, c0, x0 1203 1204 1205def _round_to_power_of_two(x): 1206 """ 1207 Round elements of the array to the nearest power of two. 1208 """ 1209 return 2**np.around(np.log2(x)) 1210 1211 1212def _autoscale(A, b, c, x0): 1213 """ 1214 Scales the problem according to equilibration from [12]. 1215 Also normalizes the right hand side vector by its maximum element. 1216 """ 1217 m, n = A.shape 1218 1219 C = 1 1220 R = 1 1221 1222 if A.size > 0: 1223 1224 R = np.max(np.abs(A), axis=1) 1225 if sps.issparse(A): 1226 R = R.toarray().flatten() 1227 R[R == 0] = 1 1228 R = 1/_round_to_power_of_two(R) 1229 A = sps.diags(R)*A if sps.issparse(A) else A*R.reshape(m, 1) 1230 b = b*R 1231 1232 C = np.max(np.abs(A), axis=0) 1233 if sps.issparse(A): 1234 C = C.toarray().flatten() 1235 C[C == 0] = 1 1236 C = 1/_round_to_power_of_two(C) 1237 A = A*sps.diags(C) if sps.issparse(A) else A*C 1238 c = c*C 1239 1240 b_scale = np.max(np.abs(b)) if b.size > 0 else 1 1241 if b_scale == 0: 1242 b_scale = 1. 1243 b = b/b_scale 1244 1245 if x0 is not None: 1246 x0 = x0/b_scale*(1/C) 1247 return A, b, c, x0, C, b_scale 1248 1249 1250def _unscale(x, C, b_scale): 1251 """ 1252 Converts solution to _autoscale problem -> solution to original problem. 1253 """ 1254 1255 try: 1256 n = len(C) 1257 # fails if sparse or scalar; that's OK. 1258 # this is only needed for original simplex (never sparse) 1259 except TypeError: 1260 n = len(x) 1261 1262 return x[:n]*b_scale*C 1263 1264 1265def _display_summary(message, status, fun, iteration): 1266 """ 1267 Print the termination summary of the linear program 1268 1269 Parameters 1270 ---------- 1271 message : str 1272 A string descriptor of the exit status of the optimization. 1273 status : int 1274 An integer representing the exit status of the optimization:: 1275 1276 0 : Optimization terminated successfully 1277 1 : Iteration limit reached 1278 2 : Problem appears to be infeasible 1279 3 : Problem appears to be unbounded 1280 4 : Serious numerical difficulties encountered 1281 1282 fun : float 1283 Value of the objective function. 1284 iteration : iteration 1285 The number of iterations performed. 1286 """ 1287 print(message) 1288 if status in (0, 1): 1289 print(" Current function value: {0: <12.6f}".format(fun)) 1290 print(" Iterations: {0:d}".format(iteration)) 1291 1292 1293def _postsolve(x, postsolve_args, complete=False): 1294 """ 1295 Given solution x to presolved, standard form linear program x, add 1296 fixed variables back into the problem and undo the variable substitutions 1297 to get solution to original linear program. Also, calculate the objective 1298 function value, slack in original upper bound constraints, and residuals 1299 in original equality constraints. 1300 1301 Parameters 1302 ---------- 1303 x : 1-D array 1304 Solution vector to the standard-form problem. 1305 postsolve_args : tuple 1306 Data needed by _postsolve to convert the solution to the standard-form 1307 problem into the solution to the original problem, including: 1308 1309 lp : A `scipy.optimize._linprog_util._LPProblem` consisting of the following fields: 1310 1311 c : 1D array 1312 The coefficients of the linear objective function to be minimized. 1313 A_ub : 2D array, optional 1314 The inequality constraint matrix. Each row of ``A_ub`` specifies the 1315 coefficients of a linear inequality constraint on ``x``. 1316 b_ub : 1D array, optional 1317 The inequality constraint vector. Each element represents an 1318 upper bound on the corresponding value of ``A_ub @ x``. 1319 A_eq : 2D array, optional 1320 The equality constraint matrix. Each row of ``A_eq`` specifies the 1321 coefficients of a linear equality constraint on ``x``. 1322 b_eq : 1D array, optional 1323 The equality constraint vector. Each element of ``A_eq @ x`` must equal 1324 the corresponding element of ``b_eq``. 1325 bounds : 2D array 1326 The bounds of ``x``, lower bounds in the 1st column, upper 1327 bounds in the 2nd column. The bounds are possibly tightened 1328 by the presolve procedure. 1329 x0 : 1D array, optional 1330 Guess values of the decision variables, which will be refined by 1331 the optimization algorithm. This argument is currently used only by the 1332 'revised simplex' method, and can only be used if `x0` represents a 1333 basic feasible solution. 1334 1335 revstack: list of functions 1336 the functions in the list reverse the operations of _presolve() 1337 the function signature is x_org = f(x_mod), where x_mod is the result 1338 of a presolve step and x_org the value at the start of the step 1339 complete : bool 1340 Whether the solution is was determined in presolve (``True`` if so) 1341 1342 Returns 1343 ------- 1344 x : 1-D array 1345 Solution vector to original linear programming problem 1346 fun: float 1347 optimal objective value for original problem 1348 slack : 1-D array 1349 The (non-negative) slack in the upper bound constraints, that is, 1350 ``b_ub - A_ub @ x`` 1351 con : 1-D array 1352 The (nominally zero) residuals of the equality constraints, that is, 1353 ``b - A_eq @ x`` 1354 """ 1355 # note that all the inputs are the ORIGINAL, unmodified versions 1356 # no rows, columns have been removed 1357 1358 (c, A_ub, b_ub, A_eq, b_eq, bounds, x0), revstack, C, b_scale = postsolve_args 1359 1360 x = _unscale(x, C, b_scale) 1361 1362 # Undo variable substitutions of _get_Abc() 1363 # if "complete", problem was solved in presolve; don't do anything here 1364 n_x = bounds.shape[0] 1365 if not complete and bounds is not None: # bounds are never none, probably 1366 n_unbounded = 0 1367 for i, bi in enumerate(bounds): 1368 lbi = bi[0] 1369 ubi = bi[1] 1370 if lbi == -np.inf and ubi == np.inf: 1371 n_unbounded += 1 1372 x[i] = x[i] - x[n_x + n_unbounded - 1] 1373 else: 1374 if lbi == -np.inf: 1375 x[i] = ubi - x[i] 1376 else: 1377 x[i] += lbi 1378 # all the rest of the variables were artificial 1379 x = x[:n_x] 1380 1381 # If there were variables removed from the problem, add them back into the 1382 # solution vector 1383 # Apply the functions in revstack (reverse direction) 1384 for rev in reversed(revstack): 1385 x = rev(x) 1386 1387 fun = x.dot(c) 1388 slack = b_ub - A_ub.dot(x) # report slack for ORIGINAL UB constraints 1389 # report residuals of ORIGINAL EQ constraints 1390 con = b_eq - A_eq.dot(x) 1391 1392 return x, fun, slack, con 1393 1394 1395def _check_result(x, fun, status, slack, con, bounds, tol, message): 1396 """ 1397 Check the validity of the provided solution. 1398 1399 A valid (optimal) solution satisfies all bounds, all slack variables are 1400 negative and all equality constraint residuals are strictly non-zero. 1401 Further, the lower-bounds, upper-bounds, slack and residuals contain 1402 no nan values. 1403 1404 Parameters 1405 ---------- 1406 x : 1-D array 1407 Solution vector to original linear programming problem 1408 fun: float 1409 optimal objective value for original problem 1410 status : int 1411 An integer representing the exit status of the optimization:: 1412 1413 0 : Optimization terminated successfully 1414 1 : Iteration limit reached 1415 2 : Problem appears to be infeasible 1416 3 : Problem appears to be unbounded 1417 4 : Serious numerical difficulties encountered 1418 1419 slack : 1-D array 1420 The (non-negative) slack in the upper bound constraints, that is, 1421 ``b_ub - A_ub @ x`` 1422 con : 1-D array 1423 The (nominally zero) residuals of the equality constraints, that is, 1424 ``b - A_eq @ x`` 1425 bounds : 2D array 1426 The bounds on the original variables ``x`` 1427 message : str 1428 A string descriptor of the exit status of the optimization. 1429 tol : float 1430 Termination tolerance; see [1]_ Section 4.5. 1431 1432 Returns 1433 ------- 1434 status : int 1435 An integer representing the exit status of the optimization:: 1436 1437 0 : Optimization terminated successfully 1438 1 : Iteration limit reached 1439 2 : Problem appears to be infeasible 1440 3 : Problem appears to be unbounded 1441 4 : Serious numerical difficulties encountered 1442 1443 message : str 1444 A string descriptor of the exit status of the optimization. 1445 """ 1446 # Somewhat arbitrary 1447 tol = np.sqrt(tol) * 10 1448 1449 if x is None: 1450 # HiGHS does not provide x if infeasible/unbounded 1451 if status == 0: # Observed with HiGHS Simplex Primal 1452 status = 4 1453 message = ("The solver did not provide a solution nor did it " 1454 "report a failure. Please submit a bug report.") 1455 return status, message 1456 1457 contains_nans = ( 1458 np.isnan(x).any() 1459 or np.isnan(fun) 1460 or np.isnan(slack).any() 1461 or np.isnan(con).any() 1462 ) 1463 1464 if contains_nans: 1465 is_feasible = False 1466 else: 1467 invalid_bounds = (x < bounds[:, 0] - tol).any() or (x > bounds[:, 1] + tol).any() 1468 invalid_slack = status != 3 and (slack < -tol).any() 1469 invalid_con = status != 3 and (np.abs(con) > tol).any() 1470 is_feasible = not (invalid_bounds or invalid_slack or invalid_con) 1471 1472 if status == 0 and not is_feasible: 1473 status = 4 1474 message = ("The solution does not satisfy the constraints within the " 1475 "required tolerance of " + "{:.2E}".format(tol) + ", yet " 1476 "no errors were raised and there is no certificate of " 1477 "infeasibility or unboundedness. Check whether " 1478 "the slack and constraint residuals are acceptable; " 1479 "if not, consider enabling presolve, adjusting the " 1480 "tolerance option(s), and/or using a different method. " 1481 "Please consider submitting a bug report.") 1482 elif status == 2 and is_feasible: 1483 # Occurs if the simplex method exits after phase one with a very 1484 # nearly basic feasible solution. Postsolving can make the solution 1485 # basic, however, this solution is NOT optimal 1486 status = 4 1487 message = ("The solution is feasible, but the solver did not report " 1488 "that the solution was optimal. Please try a different " 1489 "method.") 1490 1491 return status, message 1492