Hide keyboard shortcuts

Hot-keys on this page

r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

1"""Revised simplex method for linear programming 

2 

3The *revised simplex* method uses the method described in [1]_, except 

4that a factorization [2]_ of the basis matrix, rather than its inverse, 

5is efficiently maintained and used to solve the linear systems at each 

6iteration of the algorithm. 

7 

8.. versionadded:: 1.3.0 

9 

10References 

11---------- 

12.. [1] Bertsimas, Dimitris, and J. Tsitsiklis. "Introduction to linear 

13 programming." Athena Scientific 1 (1997): 997. 

14.. [2] Bartels, Richard H. "A stabilization of the simplex method." 

15 Journal in Numerische Mathematik 16.5 (1971): 414-434. 

16 

17""" 

18# Author: Matt Haberland 

19 

20import numpy as np 

21from scipy.linalg import solve 

22from .optimize import _check_unknown_options 

23from ._bglu_dense import LU 

24from ._bglu_dense import BGLU as BGLU 

25from scipy.linalg import LinAlgError 

26from numpy.linalg.linalg import LinAlgError as LinAlgError2 

27from ._linprog_util import _postsolve 

28from .optimize import OptimizeResult 

29 

30 

31def _phase_one(A, b, x0, callback, postsolve_args, maxiter, tol, disp, 

32 maxupdate, mast, pivot): 

33 """ 

34 The purpose of phase one is to find an initial basic feasible solution 

35 (BFS) to the original problem. 

36 

37 Generates an auxiliary problem with a trivial BFS and an objective that 

38 minimizes infeasibility of the original problem. Solves the auxiliary 

39 problem using the main simplex routine (phase two). This either yields 

40 a BFS to the original problem or determines that the original problem is 

41 infeasible. If feasible, phase one detects redundant rows in the original 

42 constraint matrix and removes them, then chooses additional indices as 

43 necessary to complete a basis/BFS for the original problem. 

44 """ 

45 

46 m, n = A.shape 

47 status = 0 

48 

49 # generate auxiliary problem to get initial BFS 

50 A, b, c, basis, x, status = _generate_auxiliary_problem(A, b, x0, tol) 

51 

52 if status == 6: 

53 residual = c.dot(x) 

54 iter_k = 0 

55 return x, basis, A, b, residual, status, iter_k 

56 

57 # solve auxiliary problem 

58 phase_one_n = n 

59 iter_k = 0 

60 x, basis, status, iter_k = _phase_two(c, A, x, basis, callback, 

61 postsolve_args, 

62 maxiter, tol, disp, 

63 maxupdate, mast, pivot, 

64 iter_k, phase_one_n) 

65 

66 # check for infeasibility 

67 residual = c.dot(x) 

68 if status == 0 and residual > tol: 

69 status = 2 

70 

71 # drive artificial variables out of basis 

72 # TODO: test redundant row removal better 

73 # TODO: make solve more efficient with BGLU? This could take a while. 

74 keep_rows = np.ones(m, dtype=bool) 

75 for basis_column in basis[basis >= n]: 

76 B = A[:, basis] 

77 try: 

78 basis_finder = np.abs(solve(B, A)) # inefficient 

79 pertinent_row = np.argmax(basis_finder[:, basis_column]) 

80 eligible_columns = np.ones(n, dtype=bool) 

81 eligible_columns[basis[basis < n]] = 0 

82 eligible_column_indices = np.where(eligible_columns)[0] 

83 index = np.argmax(basis_finder[:, :n] 

84 [pertinent_row, eligible_columns]) 

85 new_basis_column = eligible_column_indices[index] 

86 if basis_finder[pertinent_row, new_basis_column] < tol: 

87 keep_rows[pertinent_row] = False 

88 else: 

89 basis[basis == basis_column] = new_basis_column 

90 except (LinAlgError, LinAlgError2): 

91 status = 4 

92 

93 # form solution to original problem 

94 A = A[keep_rows, :n] 

95 basis = basis[keep_rows] 

96 x = x[:n] 

97 m = A.shape[0] 

98 return x, basis, A, b, residual, status, iter_k 

99 

100 

101def _get_more_basis_columns(A, basis): 

102 """ 

103 Called when the auxiliary problem terminates with artificial columns in 

104 the basis, which must be removed and replaced with non-artificial 

105 columns. Finds additional columns that do not make the matrix singular. 

106 """ 

107 m, n = A.shape 

108 

109 # options for inclusion are those that aren't already in the basis 

110 a = np.arange(m+n) 

111 bl = np.zeros(len(a), dtype=bool) 

112 bl[basis] = 1 

113 options = a[~bl] 

114 options = options[options < n] # and they have to be non-artificial 

115 

116 # form basis matrix 

117 B = np.zeros((m, m)) 

118 B[:, 0:len(basis)] = A[:, basis] 

119 

120 if (basis.size > 0 and 

121 np.linalg.matrix_rank(B[:, :len(basis)]) < len(basis)): 

122 raise Exception("Basis has dependent columns") 

123 

124 rank = 0 # just enter the loop 

125 for i in range(n): # somewhat arbitrary, but we need another way out 

126 # permute the options, and take as many as needed 

127 new_basis = np.random.permutation(options)[:m-len(basis)] 

128 B[:, len(basis):] = A[:, new_basis] # update the basis matrix 

129 rank = np.linalg.matrix_rank(B) # check the rank 

130 if rank == m: 

131 break 

132 

133 return np.concatenate((basis, new_basis)) 

134 

135 

136def _generate_auxiliary_problem(A, b, x0, tol): 

137 """ 

138 Modifies original problem to create an auxiliary problem with a trivial 

139 initial basic feasible solution and an objective that minimizes 

140 infeasibility in the original problem. 

141 

142 Conceptually, this is done by stacking an identity matrix on the right of 

143 the original constraint matrix, adding artificial variables to correspond 

144 with each of these new columns, and generating a cost vector that is all 

145 zeros except for ones corresponding with each of the new variables. 

146 

147 A initial basic feasible solution is trivial: all variables are zero 

148 except for the artificial variables, which are set equal to the 

149 corresponding element of the right hand side `b`. 

150 

151 Runnning the simplex method on this auxiliary problem drives all of the 

152 artificial variables - and thus the cost - to zero if the original problem 

153 is feasible. The original problem is declared infeasible otherwise. 

154 

155 Much of the complexity below is to improve efficiency by using singleton 

156 columns in the original problem where possible, thus generating artificial 

157 variables only as necessary, and using an initial 'guess' basic feasible 

158 solution. 

159 """ 

160 status = 0 

161 m, n = A.shape 

162 

163 if x0 is not None: 

164 x = x0 

165 else: 

166 x = np.zeros(n) 

167 

168 r = b - A@x # residual; this must be all zeros for feasibility 

169 

170 A[r < 0] = -A[r < 0] # express problem with RHS positive for trivial BFS 

171 b[r < 0] = -b[r < 0] # to the auxiliary problem 

172 r[r < 0] *= -1 

173 

174 # Rows which we will need to find a trivial way to zero. 

175 # This should just be the rows where there is a nonzero residual. 

176 # But then we would not necessarily have a column singleton in every row. 

177 # This makes it difficult to find an initial basis. 

178 if x0 is None: 

179 nonzero_constraints = np.arange(m) 

180 else: 

181 nonzero_constraints = np.where(r > tol)[0] 

182 

183 # these are (at least some of) the initial basis columns 

184 basis = np.where(np.abs(x) > tol)[0] 

185 

186 if len(nonzero_constraints) == 0 and len(basis) <= m: # already a BFS 

187 c = np.zeros(n) 

188 basis = _get_more_basis_columns(A, basis) 

189 return A, b, c, basis, x, status 

190 elif (len(nonzero_constraints) > m - len(basis) or 

191 np.any(x < 0)): # can't get trivial BFS 

192 c = np.zeros(n) 

193 status = 6 

194 return A, b, c, basis, x, status 

195 

196 # chooses existing columns appropriate for inclusion in initial basis 

197 cols, rows = _select_singleton_columns(A, r) 

198 

199 # find the rows we need to zero that we _can_ zero with column singletons 

200 i_tofix = np.isin(rows, nonzero_constraints) 

201 # these columns can't already be in the basis, though 

202 # we are going to add them to the basis and change the corresponding x val 

203 i_notinbasis = np.logical_not(np.isin(cols, basis)) 

204 i_fix_without_aux = np.logical_and(i_tofix, i_notinbasis) 

205 rows = rows[i_fix_without_aux] 

206 cols = cols[i_fix_without_aux] 

207 

208 # indices of the rows we can only zero with auxiliary variable 

209 # these rows will get a one in each auxiliary column 

210 arows = nonzero_constraints[np.logical_not( 

211 np.isin(nonzero_constraints, rows))] 

212 n_aux = len(arows) 

213 acols = n + np.arange(n_aux) # indices of auxiliary columns 

214 

215 basis_ng = np.concatenate((cols, acols)) # basis columns not from guess 

216 basis_ng_rows = np.concatenate((rows, arows)) # rows we need to zero 

217 

218 # add auxiliary singleton columns 

219 A = np.hstack((A, np.zeros((m, n_aux)))) 

220 A[arows, acols] = 1 

221 

222 # generate initial BFS 

223 x = np.concatenate((x, np.zeros(n_aux))) 

224 x[basis_ng] = r[basis_ng_rows]/A[basis_ng_rows, basis_ng] 

225 

226 # generate costs to minimize infeasibility 

227 c = np.zeros(n_aux + n) 

228 c[acols] = 1 

229 

230 # basis columns correspond with nonzeros in guess, those with column 

231 # singletons we used to zero remaining constraints, and any additional 

232 # columns to get a full set (m columns) 

233 basis = np.concatenate((basis, basis_ng)) 

234 basis = _get_more_basis_columns(A, basis) # add columns as needed 

235 

236 return A, b, c, basis, x, status 

237 

238 

239def _select_singleton_columns(A, b): 

240 """ 

241 Finds singleton columns for which the singleton entry is of the same sign 

242 as the right-hand side; these columns are eligible for inclusion in an 

243 initial basis. Determines the rows in which the singleton entries are 

244 located. For each of these rows, returns the indices of the one singleton 

245 column and its corresponding row. 

246 """ 

247 # find indices of all singleton columns and corresponding row indicies 

248 column_indices = np.nonzero(np.sum(np.abs(A) != 0, axis=0) == 1)[0] 

249 columns = A[:, column_indices] # array of singleton columns 

250 row_indices = np.zeros(len(column_indices), dtype=int) 

251 nonzero_rows, nonzero_columns = np.nonzero(columns) 

252 row_indices[nonzero_columns] = nonzero_rows # corresponding row indicies 

253 

254 # keep only singletons with entries that have same sign as RHS 

255 # this is necessary because all elements of BFS must be non-negative 

256 same_sign = A[row_indices, column_indices]*b[row_indices] >= 0 

257 column_indices = column_indices[same_sign][::-1] 

258 row_indices = row_indices[same_sign][::-1] 

259 # Reversing the order so that steps below select rightmost columns 

260 # for initial basis, which will tend to be slack variables. (If the 

261 # guess corresponds with a basic feasible solution but a constraint 

262 # is not satisfied with the corresponding slack variable zero, the slack 

263 # variable must be basic.) 

264 

265 # for each row, keep rightmost singleton column with an entry in that row 

266 unique_row_indices, first_columns = np.unique(row_indices, 

267 return_index=True) 

268 return column_indices[first_columns], unique_row_indices 

269 

270 

271def _find_nonzero_rows(A, tol): 

272 """ 

273 Returns logical array indicating the locations of rows with at least 

274 one nonzero element. 

275 """ 

276 return np.any(np.abs(A) > tol, axis=1) 

277 

278 

279def _select_enter_pivot(c_hat, bl, a, rule="bland", tol=1e-12): 

280 """ 

281 Selects a pivot to enter the basis. Currently Bland's rule - the smallest 

282 index that has a negative reduced cost - is the default. 

283 """ 

284 if rule.lower() == "mrc": # index with minimum reduced cost 

285 return a[~bl][np.argmin(c_hat)] 

286 else: # smallest index w/ negative reduced cost 

287 return a[~bl][c_hat < -tol][0] 

288 

289 

290def _display_iter(phase, iteration, slack, con, fun): 

291 """ 

292 Print indicators of optimization status to the console. 

293 """ 

294 header = True if not iteration % 20 else False 

295 

296 if header: 

297 print("Phase", 

298 "Iteration", 

299 "Minimum Slack ", 

300 "Constraint Residual", 

301 "Objective ") 

302 

303 # :<X.Y left aligns Y digits in X digit spaces 

304 fmt = '{0:<6}{1:<10}{2:<20.13}{3:<20.13}{4:<20.13}' 

305 try: 

306 slack = np.min(slack) 

307 except ValueError: 

308 slack = "NA" 

309 print(fmt.format(phase, iteration, slack, np.linalg.norm(con), fun)) 

310 

311 

312def _phase_two(c, A, x, b, callback, postsolve_args, maxiter, tol, disp, 

313 maxupdate, mast, pivot, iteration=0, phase_one_n=None): 

314 """ 

315 The heart of the simplex method. Beginning with a basic feasible solution, 

316 moves to adjacent basic feasible solutions successively lower reduced cost. 

317 Terminates when there are no basic feasible solutions with lower reduced 

318 cost or if the problem is determined to be unbounded. 

319 

320 This implementation follows the revised simplex method based on LU 

321 decomposition. Rather than maintaining a tableau or an inverse of the 

322 basis matrix, we keep a factorization of the basis matrix that allows 

323 efficient solution of linear systems while avoiding stability issues 

324 associated with inverted matrices. 

325 """ 

326 m, n = A.shape 

327 status = 0 

328 a = np.arange(n) # indices of columns of A 

329 ab = np.arange(m) # indices of columns of B 

330 if maxupdate: 

331 # basis matrix factorization object; similar to B = A[:, b] 

332 B = BGLU(A, b, maxupdate, mast) 

333 else: 

334 B = LU(A, b) 

335 

336 for iteration in range(iteration, iteration + maxiter): 

337 

338 if disp or callback is not None: 

339 if phase_one_n is not None: 

340 phase = 1 

341 x_postsolve = x[:phase_one_n] 

342 else: 

343 phase = 2 

344 x_postsolve = x 

345 x_o, fun, slack, con, _ = _postsolve(x_postsolve, 

346 postsolve_args, 

347 tol=tol, copy=True) 

348 

349 if callback is not None: 

350 res = OptimizeResult({'x': x_o, 'fun': fun, 'slack': slack, 

351 'con': con, 'nit': iteration, 

352 'phase': phase, 'complete': False, 

353 'status': 0, 'message': "", 

354 'success': False}) 

355 callback(res) 

356 else: 

357 _display_iter(phase, iteration, slack, con, fun) 

358 

359 bl = np.zeros(len(a), dtype=bool) 

360 bl[b] = 1 

361 

362 xb = x[b] # basic variables 

363 cb = c[b] # basic costs 

364 

365 try: 

366 v = B.solve(cb, transposed=True) # similar to v = solve(B.T, cb) 

367 except LinAlgError: 

368 status = 4 

369 break 

370 

371 # TODO: cythonize? 

372 c_hat = c - v.dot(A) # reduced cost 

373 c_hat = c_hat[~bl] 

374 # Above is much faster than: 

375 # N = A[:, ~bl] # slow! 

376 # c_hat = c[~bl] - v.T.dot(N) 

377 # Can we perform the multiplication only on the nonbasic columns? 

378 

379 if np.all(c_hat >= -tol): # all reduced costs positive -> terminate 

380 break 

381 

382 j = _select_enter_pivot(c_hat, bl, a, rule=pivot, tol=tol) 

383 u = B.solve(A[:, j]) # similar to u = solve(B, A[:, j]) 

384 

385 i = u > tol # if none of the u are positive, unbounded 

386 if not np.any(i): 

387 status = 3 

388 break 

389 

390 th = xb[i]/u[i] 

391 l = np.argmin(th) # implicitly selects smallest subscript 

392 th_star = th[l] # step size 

393 

394 x[b] = x[b] - th_star*u # take step 

395 x[j] = th_star 

396 B.update(ab[i][l], j) # modify basis 

397 b = B.b # similar to b[ab[i][l]] = j 

398 else: 

399 status = 1 

400 

401 return x, b, status, iteration 

402 

403 

404def _linprog_rs(c, c0, A, b, x0, callback, postsolve_args, 

405 maxiter=5000, tol=1e-12, disp=False, 

406 maxupdate=10, mast=False, pivot="mrc", 

407 **unknown_options): 

408 """ 

409 Solve the following linear programming problem via a two-phase 

410 revised simplex algorithm.:: 

411 

412 minimize: c @ x 

413 

414 subject to: A @ x == b 

415 0 <= x < oo 

416 

417 Parameters 

418 ---------- 

419 c : 1-D array 

420 Coefficients of the linear objective function to be minimized. 

421 c0 : float 

422 Constant term in objective function due to fixed (and eliminated) 

423 variables. (Currently unused.) 

424 A : 2-D array 

425 2-D array which, when matrix-multiplied by ``x``, gives the values of 

426 the equality constraints at ``x``. 

427 b : 1-D array 

428 1-D array of values representing the RHS of each equality constraint 

429 (row) in ``A_eq``. 

430 x0 : 1-D array, optional 

431 Starting values of the independent variables, which will be refined by 

432 the optimization algorithm. For the revised simplex method, these must 

433 correspond with a basic feasible solution. 

434 callback : callable, optional 

435 If a callback function is provided, it will be called within each 

436 iteration of the algorithm. The callback function must accept a single 

437 `scipy.optimize.OptimizeResult` consisting of the following fields: 

438 

439 x : 1-D array 

440 Current solution vector. 

441 fun : float 

442 Current value of the objective function ``c @ x``. 

443 success : bool 

444 True only when an algorithm has completed successfully, 

445 so this is always False as the callback function is called 

446 only while the algorithm is still iterating. 

447 slack : 1-D array 

448 The values of the slack variables. Each slack variable 

449 corresponds to an inequality constraint. If the slack is zero, 

450 the corresponding constraint is active. 

451 con : 1-D array 

452 The (nominally zero) residuals of the equality constraints, 

453 that is, ``b - A_eq @ x``. 

454 phase : int 

455 The phase of the algorithm being executed. 

456 status : int 

457 For revised simplex, this is always 0 because if a different 

458 status is detected, the algorithm terminates. 

459 nit : int 

460 The number of iterations performed. 

461 message : str 

462 A string descriptor of the exit status of the optimization. 

463 postsolve_args : tuple 

464 Data needed by _postsolve to convert the solution to the standard-form 

465 problem into the solution to the original problem. 

466 

467 Options 

468 ------- 

469 maxiter : int 

470 The maximum number of iterations to perform in either phase. 

471 tol : float 

472 The tolerance which determines when a solution is "close enough" to 

473 zero in Phase 1 to be considered a basic feasible solution or close 

474 enough to positive to serve as an optimal solution. 

475 disp : bool 

476 Set to ``True`` if indicators of optimization status are to be printed 

477 to the console each iteration. 

478 maxupdate : int 

479 The maximum number of updates performed on the LU factorization. 

480 After this many updates is reached, the basis matrix is factorized 

481 from scratch. 

482 mast : bool 

483 Minimize Amortized Solve Time. If enabled, the average time to solve 

484 a linear system using the basis factorization is measured. Typically, 

485 the average solve time will decrease with each successive solve after 

486 initial factorization, as factorization takes much more time than the 

487 solve operation (and updates). Eventually, however, the updated 

488 factorization becomes sufficiently complex that the average solve time 

489 begins to increase. When this is detected, the basis is refactorized 

490 from scratch. Enable this option to maximize speed at the risk of 

491 nondeterministic behavior. Ignored if ``maxupdate`` is 0. 

492 pivot : "mrc" or "bland" 

493 Pivot rule: Minimum Reduced Cost (default) or Bland's rule. Choose 

494 Bland's rule if iteration limit is reached and cycling is suspected. 

495 unknown_options : dict 

496 Optional arguments not used by this particular solver. If 

497 `unknown_options` is non-empty a warning is issued listing all 

498 unused options. 

499 

500 Returns 

501 ------- 

502 x : 1-D array 

503 Solution vector. 

504 status : int 

505 An integer representing the exit status of the optimization:: 

506 

507 0 : Optimization terminated successfully 

508 1 : Iteration limit reached 

509 2 : Problem appears to be infeasible 

510 3 : Problem appears to be unbounded 

511 4 : Numerical difficulties encountered 

512 5 : No constraints; turn presolve on 

513 6 : Guess x0 cannot be converted to a basic feasible solution 

514 

515 message : str 

516 A string descriptor of the exit status of the optimization. 

517 iteration : int 

518 The number of iterations taken to solve the problem. 

519 """ 

520 

521 _check_unknown_options(unknown_options) 

522 

523 messages = ["Optimization terminated successfully.", 

524 "Iteration limit reached.", 

525 "The problem appears infeasible, as the phase one auxiliary " 

526 "problem terminated successfully with a residual of {0:.1e}, " 

527 "greater than the tolerance {1} required for the solution to " 

528 "be considered feasible. Consider increasing the tolerance to " 

529 "be greater than {0:.1e}. If this tolerance is unnaceptably " 

530 "large, the problem is likely infeasible.", 

531 "The problem is unbounded, as the simplex algorithm found " 

532 "a basic feasible solution from which there is a direction " 

533 "with negative reduced cost in which all decision variables " 

534 "increase.", 

535 "Numerical difficulties encountered; consider trying " 

536 "method='interior-point'.", 

537 "Problems with no constraints are trivially solved; please " 

538 "turn presolve on.", 

539 "The guess x0 cannot be converted to a basic feasible " 

540 "solution. " 

541 ] 

542 

543 if A.size == 0: # address test_unbounded_below_no_presolve_corrected 

544 return np.zeros(c.shape), 5, messages[5], 0 

545 

546 x, basis, A, b, residual, status, iteration = ( 

547 _phase_one(A, b, x0, callback, postsolve_args, 

548 maxiter, tol, disp, maxupdate, mast, pivot)) 

549 

550 if status == 0: 

551 x, basis, status, iteration = _phase_two(c, A, x, basis, callback, 

552 postsolve_args, 

553 maxiter, tol, disp, 

554 maxupdate, mast, pivot, 

555 iteration) 

556 

557 return x, status, messages[status].format(residual, tol), iteration